{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quijote $C_{\\ell}$ of the LOS pecululiar velocity covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import healpy\n",
    "from h5py import File\n",
    "\n",
    "import csiborgtools\n",
    "import healpy as hp\n",
    "from scipy.special import spherical_jn\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.integrate import quad\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "k, Pk = np.load(\"/mnt/users/rstiskalek/BayesianBulkFlows/data/pk_fiducial.npy\").T\n",
    "loglog_Pk = interp1d(np.log(k), np.log(Pk))\n",
    "\n",
    "def Pk(k):\n",
    "    return np.exp(loglog_Pk(np.log(k)))\n",
    "\n",
    "\n",
    "def djn(n, x):\n",
    "    return (n / x) * spherical_jn(n, x) - spherical_jn(n + 1, x)\n",
    "\n",
    "\n",
    "def C_ell_theory_equal_radius(ell, r):\n",
    "    dDdtau = 35.5\n",
    "    \n",
    "    norm = dDdtau**2 \n",
    "    integrand = lambda k: norm * Pk(k) * djn(ell, k * r)**2\n",
    "    return quad(integrand, 0.0001, 10, limit=1000, limlst=5)[0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Num. simulations:   50\n",
      "Num. observers:     8\n",
      "Num. radii:         5\n",
      "Radii:              [100.  112.5 125.  137.5 150. ]\n",
      "Dividing the dipole by a factor of 9:\n"
     ]
    }
   ],
   "source": [
    "num_simulations = None\n",
    "load_skymaps = False\n",
    "\n",
    "with File(\"/mnt/extraspace/rstiskalek/BBF/Quijote_Cell/C_ell_fiducial.h5\", 'r') as f:\n",
    "    num_observers = f.attrs[\"num_observers\"]\n",
    "    num_radii = f.attrs[\"num_radii\"]\n",
    "    npix_per_map = f.attrs[\"npix_per_map\"]\n",
    "    if num_simulations is None:\n",
    "        num_simulations = f.attrs[\"num_simulations\"]\n",
    "\n",
    "    ell = f[\"ell\"][:]\n",
    "    radii = f[\"radii\"][:]\n",
    "\n",
    "    Cell = []\n",
    "    if load_skymaps:\n",
    "        skymap = []\n",
    "\n",
    "    for n in range(num_simulations):\n",
    "        grp_simulation = f[f\"nsim_{n}\"]\n",
    "\n",
    "        for i in range(num_observers):\n",
    "            grp_observer = grp_simulation[f\"observer_{i}\"]\n",
    "\n",
    "            Cell_ = np.full((num_radii, len(ell)), np.nan)\n",
    "            if load_skymaps:\n",
    "                skymap_ = np.full((num_radii, npix_per_map), np.nan)\n",
    "            else:\n",
    "                skymap_ = None\n",
    "\n",
    "            for j in range(num_radii):\n",
    "                Cell_[j] = grp_observer[f\"C_ell_{j}\"][...]\n",
    "                if load_skymaps:\n",
    "                    skymap_[j] = grp_observer[f\"skymap_{j}\"][...]\n",
    "\n",
    "\n",
    "            Cell.append(Cell_)\n",
    "            if load_skymaps:\n",
    "                skymap.append(skymap_)\n",
    "\n",
    "\n",
    "Cell = np.stack(Cell)\n",
    "if load_skymaps:\n",
    "    skymap = np.stack(skymap)\n",
    "\n",
    "\n",
    "print(f\"Num. simulations:   {num_simulations}\")\n",
    "print(f\"Num. observers:     {num_observers}\")\n",
    "print(f\"Num. radii:         {num_radii}\")\n",
    "print(f\"Radii:              {radii}\")\n",
    "\n",
    "\n",
    "\n",
    "print(\"Dividing the dipole by a factor of 9:\")\n",
    "Cell[..., 1] /= 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# healpy.mollview(skymap[2, 4], title=None, cmap=\"coolwarm\", unit=r\"$V_{\\rm rad} ~ [{\\rm km/s}]$\", )\n",
    "# plt.savefig(\"../../plots/skymap.png\", dpi=450)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1113107/3790696556.py:17: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
      "  the requested tolerance from being achieved.  The error may be \n",
      "  underestimated.\n",
      "  return quad(integrand, 0.0001, 10, limit=1000, limlst=5)[0]\n"
     ]
    }
   ],
   "source": [
    "k = -1\n",
    "\n",
    "\n",
    "ytheory = [C_ell_theory_equal_radius(i, radii[k]) for i in ell]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAHKCAYAAAANcD9RAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAACDeElEQVR4nO3dd3TV9f3H8ee9N7k3e+9F2HsThqA4kKFFce89quIeP7VWsbUVrbV1kLqtW3GBoy5EQGXIBtlhhRCy97xJ7v3+/giJRFaAO0LyepyTc5r7Xe9LYu6rn2kyDMNAREREpIMze7sAERERkbZAoUhEREQEhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIREREBFAoEhEREQEUikREREQAhSIRERERQKFIpN169NFHMZlMFBYWersUj2t67/t64403MJlM7Ny584judSTXHei5x6q19+zIP28RV1EoEmlDNm7cyM0330y3bt0ICAggODiY0aNH88ILL9DQ0ODt8li0aBGPPvoopaWlR3xtW39v0jorVqxg4sSJhISEEBwczPjx41m9erW3yxJxCYUikTbipZdeYtCgQcyaNYtzzjmHZ599lmnTpuHj48Mtt9zCM8884+0SWbRoEX/5y1+OOBS1hfd2xRVXUFNTQ6dOnTxyXXu0cuVKxowZw/bt25k2bRqPPPIIGRkZjB07ls2bN3u7PJFj5uPtAkQE3nvvPW6++WbOPPNMPvjgAwIDA5uP3Xvvvbz33ntERUV5scKj11bem8ViwWKxeOy69ujhhx/G39+fxYsXExkZCcDll19Ojx49+NOf/sQnn3zi5QpFjo1aikS8LDc3l6lTp9KnTx8+/vjjFqGhyaWXXsr48eOP6v6FhYVceOGFhISEEBkZyR133EFtbW3z8auvvprU1NT9rvv9WJZHH32U++67D4DOnTtjMpkOO9bG3e8N4OeffyYtLQ0/Pz+6du3KSy+9dMDzfj826OOPP8ZkMrFgwYL9zn3ppZcwmUysW7fuoGOKWvtcgOzsbK699lpiY2Ox2Wz07duX119//ajfy6Ec6uc9b948TCYTs2bN2u+69957D5PJxOLFiw96759++olx48Y1ByKA+Ph4xo4dy5dffkllZeUR1yvSlqilSMTLnn32WUpLS3n//fex2Wwuv/+FF15Iamoq06dPZ8mSJTz33HOUlJTw1ltvHdF9zj33XLZs2cL777/Pv//97+bWnejo6INe4+739uuvvzJ+/Hiio6N59NFHaWhoYNq0acTGxh722jPPPJOgoCA+/PBDxo4d2+LYzJkz6du3L/369WP58uXH9Ny8vDxGjhyJyWTi1ltvJTo6mq+//prrrruO8vJy7rzzzmN+L/s61M/75JNPJjk5mXfffZdzzjmnxXXvvvsuXbt2ZdSoUQe9t91ux9/ff7/XAwICqKurY926dYwcOfKI6hVpUwwR8aqEhAQjOTnZcDgcLr3vtGnTDMA466yzWrx+yy23GICxZs0awzAM46qrrjI6dep00Ov39dRTTxmAsWPHjlbV4K731mTKlCmGn5+fkZmZ2fzahg0bDIvFsl/t//3vf/er/ZJLLjFiYmKMhoaG5tdycnIMs9ls/PWvfz3odUfy3Ouuu86Ij483CgsLW7x+8cUXG6GhoUZ1dfUR3/NAWvvzfvDBBw2bzWaUlpY2n5Ofn2/4+PgY06ZNO+Qz+vfvb/To0aPFv5fdbjdSUlIMwPj4448PW6dIW6buMxEv2rZtG3v27OHMM8/EbHbPf45Tp05t8f1tt90GwFdffeWW5zVx93tzOBx8++23TJkyhZSUlObXe/fuzYQJE1p1j4suuoj8/Hzmz5/f/NrHH3+M0+nkoosuOubnGobBJ598wuTJkzEMg8LCwuavCRMmUFZWxsqVK13yXpoc7ud95ZVXYrfb+fjjj5vPmTlzJg0NDVx++eWHvPctt9zCli1buO6669iwYQPr1q3jyiuvJCcnB4CampojqlWkrVEoEvGirKwsALp06dLqa6655hrS0tJwOp2tOr979+4tvu/atStms/mI1+s5Ukf63l544QWGDBmCr68vjz766GHPLygooKamZr/3B9CzZ89WPXPixImEhoYyc+bM5tdmzpzJoEGD6NGjxzE/t6CggNLSUl5++WWio6NbfF1zzTUA5Ofnu+S9NDncz7tXr16kpaXx7rvvNp/z7rvvMnLkSLp163bIe99000386U9/4r333qNv377079+fbdu28X//938ABAUFHVGtIm2NQpGIFzWtz1NfX9+q8xctWkRpaSkmk6nV1/ze7xcCPNjCgA6H46ju3+RI31t8fDyPPvoo55133jE990jYbDamTJnCrFmzaGhoIDs7m4ULFx60lehINQXXyy+/nDlz5hzwa/To0S551sEc6Od75ZVXsmDBAnbv3s22bdtYsmTJYVuJmvz9738nLy+Pn376ibVr17Js2bLm93mwIClyvNBAaxEvavp/5uvWrWvV+Z999hlnn302M2bMaPXA5YyMDDp37tz8/datW3E6nc0zzsLDww+47lBmZuZ+rx3Jas1H+t6mTJkCtL5bLzo6Gn9/fzIyMvY7diRr5lx00UW8+eabzJ07l40bN2IYxiFD0ZE8Nzo6muDgYBwOB+PGjTvoPR0Oh0veCxz+5w1w8cUXc/fdd/P+++9TU1ODr6/vEQXB8PBwxowZ0/z9999/T1JSEr169TqiWkXaGrUUiXhRamoqaWlpfPTRR6xcuXK/4xUVFTz99NPN38+fPx+Hw8Hpp5/e6mekp6e3+P75558HYNKkSUBj90pZWRlr165tPicnJ+eA07abptS3ZvHGI31vR8pisTBhwgRmz57Nrl27ml/fuHEj3377bavvM27cOCIiIpg5cyYzZ85k+PDhLULFsTzXYrFw3nnn8cknnxwwHBYUFLj0vcDhf94AUVFRTJo0iXfeeYd3332XiRMnHvVaUTNnzmTZsmXceeedbhsXJ+IxXh7oLdLhrV692ggNDTX8/PyMW265xXjhhReM559/3rjmmmuM0NBQ48wzzzQMwzAaGhqMpKQk47LLLjMyMjIOe9+m2Uj9+/c3Jk+ebKSnpxuXX365ARiXXnpp83mFhYVGYGCg0aVLF+OZZ54xHn/8cSM5OdkYMmTIfrOeli5dagDGGWecYbz11lvG+++/b1RWVh7ze9vXH//4x8POgmqyZs0aw8/Pz0hJSTGeeOIJ429/+5sRGxtrDBgwoFWzz5pcf/31RlBQkGEymYynn376sNcdyXNzc3ONTp06GQEBAcYdd9xhvPTSS8b06dONCy64wAgPDz+qex5Ia3/eTT7++GMDMABj5syZh72/YRjGggULjNNOO8148sknjVdffdW4/vrrDYvFYkycONGor69v1T1E2jKFIpE2YOvWrcY111xjJCYmGr6+vkZsbKwxcuRI49FHHzW2bt1qGIZhZGdnGykpKcbVV1/dqns2fUhu2LDBOP/8843g4GAjPDzcuPXWW42ampoW53733XdGv379DKvVavTs2dN45513Djgl3zAM47HHHjMSExMNs9ncqun5rXlv+zqSUGQYjR/UQ4cONaxWq9GlSxfjxRdfPGDthwpFc+bMMQDDZDIZWVlZrbqutc81DMPIy8szpk6daiQnJxu+vr5GXFyccdpppxkvv/zyUd/z947k520YjVPpw8PDjdDQ0AMeP5CtW7ca48ePN6KiogybzWb06tXLmD59umG321t1vUhbZzIMw/Bgw5SIHKUtW7YwaNAgNm3a1GLadntz0003ERcX16oZaHL0GhoaSEhIYPLkybz22mveLkekTVAHsMhxYuXKlfTo0YPk5GTmzp3r7XJcrqGhgdraWhwOR4v/Le4xe/ZsCgoKuPLKK71dikiboVAkchyoq6vjq6++Yvz48Zx88smUl5d7uySX+9vf/oa/vz+vvvoqf//73/H39+ftt9/2dlntzi+//MIrr7zC3XffzeDBg/fb4kSkI1P3mYhIB3L11VfzzjvvMGjQIN544w369evn7ZJE2gyFIhERERHUfSYiIiICKBSJiIiIANrmo5nT6WTPnj0EBwcf0VYGIiIi4j2GYVBRUUFCQsIxr6quULTXnj17SE5O9nYZIiIichSysrJISko6pnt0+FCUnp5Oenp6847eWVlZhISEeLkqERERaY3y8nKSk5MJDg4+5ntp9tle5eXlhIaGUlZWplAkIiJynHDl57cGWouIiIigUCQiIiICKBSJiIiIABpoLSIiXuBwOKivr/d2GXIc8PX1xWKxeORZCkUiIuIxhmGQm5tLaWmpt0uR40hYWBhxcXFuX0dQoUhERDymKRDFxMQQEBCgxXLlkAzDoLq6mvz8fADi4+Pd+jyFIhER8QiHw9EciCIjI71djhwn/P39AcjPzycmJsatXWkaaC0iIh7RNIYoICDAy5XI8abpd8bd49AUikRExKPUZSZHylO/MwpFIiIiIigUiYiIeMzOnTsxmUysXr36gN+Ld3X4UJSenk6fPn1IS0vzdikiItKGZWVlce2115KQkIDVaqVTp07ccccdFBUVtfoeycnJ5OTk0K9fvwN+3xonn3wyd95555GWL63Q4UPR1KlT2bBhA8uWLfN2KSIi0kZt376dYcOGkZGRwfvvv8/WrVt58cUXmTt3LqNGjaK4uLhV97FYLMTFxeHj43PA78W7OnwoEhEROZypU6ditVr57rvvGDt2LCkpKUyaNInvv/+e7OxsHnroIaBxQPDs2bNbXBsWFsYbb7wBtK77bMGCBQwfPhybzUZ8fDwPPPAADQ0NAFx99dUsWLCAZ599FpPJhMlkYufOnQCsW7eOSZMmERQURGxsLFdccQWFhYXu/GdpdxSKRETEawzDoLquwStfhmG0qsbi4mK+/fZbbrnlluY1c5rExcVx2WWXMXPmzFbf71Cys7M544wzSEtLY82aNbzwwgu89tpr/O1vfwPg2WefZdSoUdxwww3k5OSQk5NDcnIypaWlnHrqqQwePJjly5fzzTffkJeXx4UXXnjMNXUkaq/7HVf8UouISOvU1Dvo88i3Xnn2hr9OIMB6+I/BjIwMDMOgd+/eBzzeu3dvSkpKKCgoOOaa/vOf/5CcnMyMGTMwmUz06tWLPXv2cP/99/PII48QGhqK1WolICCAuLi45utmzJjB4MGDefzxx5tfe/3110lOTmbLli306NHjmGvrCNRS9Ds5ZbXeLkFERNqgw/2fZqvVeszP2LhxI6NGjWqxLs/o0aOprKxk9+7dB71uzZo1zJs3j6CgoOavXr16AbBt27ZjrqujUEvR7+wsqqJHkoHZrMXFRETczd/Xwoa/TvDas1ujW7dumEwmNm7cyDnnnLPf8Y0bNxIdHU1YWBgmk2m/8OTuVZgBKisrmTx5Mk8++eR+x9y9X1h7olD0O/Z6J7tLakiJ1DL0IiLuZjKZWtWF5U2RkZGcfvrp/Oc//+Guu+5qMa4oNzeXd999l6lTpwIQHR1NTk5O8/GMjAyqq6tb/azevXvzySefYBhGc2vRwoULCQ4OJikpCWhskXI4HC2uGzJkCJ988gmpqamayXYM1H12ADuKqmhwOL1dhoiItBEzZszAbrczYcIEfvzxR7Kysvjmm284/fTT6dGjB4888ggAp556KjNmzGDVqlUsX76cm266CV9f31Y/55ZbbiErK4vbbruNTZs28dlnnzFt2jTuvvtuzObGj+zU1FR++eUXdu7cSWFhIU6nk6lTp1JcXMwll1zCsmXL2LZtG99++y3XXHPNfgFKDk6h6ADqG5xkFrc+2YuISPvWvXt3li1bRpcuXbjwwgvp1KkTkyZNokePHixcuJCgoCAAnn76aZKTkznxxBO59NJLuffee49oA9zExES++uorli5dysCBA7npppu47rrr+POf/9x8zr333ovFYqFPnz5ER0eza9cuEhISWLhwIQ6Hg/Hjx9O/f3/uvPNOwsLCmsOUHJ7J0HQrAMrLywkNDWX2LxkEBgVjsZgY3TUKq49+mUREXKG2tpYdO3bQuXNn/Pz8vF3OMZs2bRr/+te/mDNnDiNHjjyqe2zevJlevXqRkZFBt27dXFxh+3Go352mz++ysjJCQkKO6Tn6xD8Ih8NgZ1GVt8sQEZE26i9/+QvPPfccS5Yswek88iEXxcXFfPzxx4SEhJCcnOyGCuVIaTTWIewuqSYlIgC/Vs5QEBGRjuWaa6456muvu+46VqxYwQsvvIDNZnNhVXK0FIoOwemE7QVV9Ek4tuY4ERGR35s1a5a3S5DfUffZYeSU1VBd1+DtMkRERMTNFIoOwzBgW77GFomIiLR3CkWtkFdeS3mt+1ckFREREe9RKGqlbfmV3i5BRERE3EihqJWKKusoqarzdhkiIiLiJh0+FKWnp9OnTx/S0tIOe+62ArUWiYiItFcdPhRNnTqVDRs2sGzZssOeW1pdT0GF3QNViYjI8cRkMjF79mxvl3FErr76aqZMmeLtMtqUDh+KjpRai0REOp7DBYicnBwmTZrkuYKOwM6dOzGZTKxevdrbpbR5CkVHqLK2gdyyWm+XISIibUhcXJzXV6U2DIOGhra/rl5dXdsdn6tQdBS2F1SifXRFRKTJvt1nTS0zn376KaeccgoBAQEMHDiQxYsXt7jm559/5sQTT8Tf35/k5GRuv/12qqp+Wxfv7bffZtiwYQQHBxMXF8ell15Kfn5+8/H58+djMpn4+uuvGTp0KDabjZ9//nm/2jp37gzA4MGDMZlMnHzyyS2O//Of/yQ+Pp7IyEimTp1Kff1vS9DY7XbuvfdeEhMTCQwMZMSIEcyfP7/F9Z988gl9+/bFZrORmprK008/3eJ4amoqjz32GFdeeSUhISHceOONnHrqqdx6660tzisoKMBqtTJ37txD/2O7kULRUaiuc5BdWuPtMkRE2o2qqqqDftXW1rb63Jqamlad6wkPPfQQ9957L6tXr6ZHjx5ccsklzS0527ZtY+LEiZx33nmsXbuWmTNn8vPPP7cICvX19Tz22GOsWbOG2bNns3PnTq6++ur9nvPAAw/wxBNPsHHjRgYMGLDf8aVLlwLw/fffk5OTw6efftp8bN68eWzbto158+bx5ptv8sYbb/DGG280H7/11ltZvHgxH3zwAWvXruWCCy5g4sSJZGRkALBixQouvPBCLr74Yn799VceffRRHn744Rb3gMbgNXDgQFatWsXDDz/M9ddfz3vvvYfd/ts43XfeeYfExEROPfXUI/63dhlDDMMwjLKyMgMwZv+SYcxZn3vYrx+35BsNDqe3yxYROW7U1NQYGzZsMGpqavY7Bhz064wzzmhxbkBAwEHPHTt2bItzo6KiDnjekbrqqquMs88++6DHAWPWrFmGYRjGjh07DMB49dVXm4+vX7/eAIyNGzcahmEY1113nXHjjTe2uMdPP/1kmM3mA/77GIZhLFu2zACMiooKwzAMY968eY2fW7NnH7L2pnpWrVq133vq1KmT0dDQ0PzaBRdcYFx00UWGYRhGZmamYbFYjOzs7BbXnXbaacaDDz5oGIZhXHrppcbpp5/e4vh9991n9OnTp/n7Tp06GVOmTGlxTk1NjREeHm7MnDmz+bUBAwYYjz766AHfw6F+d5o+v8vKyg72T9Bqaik6SvZ6J7tLqr1dhoiItFH7ttrEx8cDNHd/rVmzhjfeeIOgoKDmrwkTJuB0OtmxYwfQ2AozefJkUlJSCA4OZuzYsQDs2rWrxXOGDRt21DX27dsXi8XSos6mGn/99VccDgc9evRoUeeCBQvYtm0bABs3bmT06NEt7jl69GgyMjJwOBwHrdHPz48rrriC119/HYCVK1eybt26A7aEeZKPV59+nNtRWEVCmD++FmVLEZFjUVl58Jm9+35oAy3G1fye2dzy7/HOnTuPqa5j4evr2/y/TSYTAE6nE2h8v3/84x+5/fbb97suJSWFqqoqJkyYwIQJE3j33XeJjo5m165dTJgwYb+ByoGBgS6psanOfWu0WCysWLFiv59BUFDQET3nQDVef/31DBo0iN27d/Pf//6XU089lU6dOh3hO3AthaJj0OAwyCyqplvMkf1yiIhIS0fywe6ucz1pyJAhbNiwgW7duh3w+K+//kpRURFPPPEEycnJACxfvvyonmW1WgFatNy0xuDBg3E4HOTn53PiiSce8JzevXuzcOHCFq8tXLiQHj167Bekfq9///4MGzaMV155hffee48ZM2YcUX3uoFB0jLJKqkmO8Mfmc+gfvoiIHN/Kysr2W+snMjKyObQcifvvv5+RI0dy6623cv311xMYGMiGDRuYM2cOM2bMICUlBavVyvPPP89NN93EunXreOyxx46q7piYGPz9/fnmm29ISkrCz8+P0NDQw17Xo0cPLrvsMq688kqefvppBg8eTEFBAXPnzmXAgAGceeaZ3HPPPaSlpfHYY49x0UUXsXjxYmbMmMF//vOfVtV2/fXXc+uttxIYGMg555xzVO/PldTvc4wcDoOdhRpbJCLS3s2fP5/Bgwe3+PrLX/5yVPcaMGAACxYsYMuWLZx44okMHjyYRx55hISEBACio6N54403+Oijj+jTpw9PPPEE//znP4/qWT4+Pjz33HO89NJLJCQkcPbZZ7f62v/+979ceeWV3HPPPfTs2ZMpU6awbNkyUlJSgMYWrw8//JAPPviAfv368cgjj/DXv/611WODLrnkEnx8fLjkkkvw8/M7mrfnUibD0II7AOXl5YSGhjL7lwwCg4KP6FqzGU7oGoWfr1qLREQOpra2lh07dtC5c+c28QEo3rdz5066du3KsmXLGDJkyEHPO9TvTtPnd1lZGSEhIcdUj1qKXMDp1PYfIiIirVVfX09ubi5//vOfGTly5CEDkScpFLlIblktVfa2v7y6iIiIty1cuJD4+HiWLVvGiy++6O1ymmmgtYsYRmNr0YCkMG+XIiIi0qadfPLJbXK7LLUUuVB+uZ2ymvrDnygiIiJtjkKRi2lskYjIobXFFgJp2zz1O6NQ5GLFlXUUV9Ud/kQRkQ6mafXk6motYyJHpul35vcrcLuaxhS5wbaCSiICI7xdhohIm2KxWAgLC2vepiMgIKB5+wuRAzEMg+rqavLz8wkLCzvsKtnHSqHIDcqq68mvqCUmWOtwiIjsKy4uDjj0/mUivxcWFtb8u+NOCkVusi2/iuggm/5fkIjIPkwmE/Hx8cTExFBfr4kpcni+vr5ubyFqolDkJlX2BnLLa4kP9fd2KSIibY7FYvHYB51Ia2mgtRttL6jC6dQsCxERkeOBQpEb1dQ5yC6t8XYZIiIi0goKRW62o7AKh1qLRERE2jyFIjera3Cyq1hrcoiIiLR17WqgdWpqKiEhIZjNZsLDw5k3b563SwIgs6iKpHB/fC3KoCIiIm1VuwpFAIsWLSIoKMjbZbTQ4DDILKqmW0zbqktERER+o6YLD8kqrsbe4PB2GSIiInIQbSYU/fjjj0yePJmEhARMJhOzZ8/e75z09HRSU1Px8/NjxIgRLF26tMVxk8nE2LFjSUtL49133/VQ5a3jcBrsKKzydhkiIiJyEG0mFFVVVTFw4EDS09MPeHzmzJncfffdTJs2jZUrVzJw4EAmTJjQYqn4n3/+mRUrVvD555/z+OOPs3btWk+V3yp7SmuoqVNrkYiISFvUZkLRpEmT+Nvf/sY555xzwOP/+te/uOGGG7jmmmvo06cPL774IgEBAbz++uvN5yQmJgIQHx/PGWecwcqVKw/6PLvdTnl5eYsvd3M6GzeLFRERkbanzYSiQ6mrq2PFihWMGzeu+TWz2cy4ceNYvHgx0NjSVFFRAUBlZSU//PADffv2Peg9p0+fTmhoaPNXcnKye9/EXnnltVTaGzzyLBEREWm94yIUFRYW4nA4iI2NbfF6bGwsubm5AOTl5TFmzBgGDhzIyJEjufLKK0lLSzvoPR988EHKysqav7Kystz6HpoYBmzLV2uRiIhIW9NupuR36dKFNWvWtPp8m82GzWZzY0UHV1Bhp6y6ntAAX688X0RERPZ3XLQURUVFYbFYyMvLa/F6Xl4ecXFxXqrq2GzV2CIREZE25bgIRVarlaFDhzJ37tzm15xOJ3PnzmXUqFFerOzolVTVUVRp93YZIiIislebCUWVlZWsXr2a1atXA7Bjxw5Wr17Nrl27ALj77rt55ZVXePPNN9m4cSM333wzVVVVXHPNNcf03PT0dPr06dM8/mjlogXHdL8jsa1A6xaJiIi0FSbDMNrEFu7z58/nlFNO2e/1q666ijfeeAOAGTNm8NRTT5Gbm8ugQYN47rnnGDFihEueX15eTmhoKEHBobzwyRziElNcct/DGZAUSkyIn0eeJSIi0t40fX6XlZUREhJyTPdqM6HI25r+UQG69x3AM29/jtXm/rASYLMwqkskJpPJ7c8SERFpb1wZitpM91lbERwaTsb6tcz4+0MeeV613UFOWa1HniUiIiIHp1D0O/c89i9MJhNff/Iu33z6nkeeub2gCodTDXYiIiLe1OFD0e8HWg8aMYarbrsfgOcee5CtG391ew219Q5WZ5XQ4HC6/VkiIiJyYB0+FE2dOpUNGzawbNmy5tcuueF2Ro49nbjEZHx9rR6po6SqntVZpQpGIiIiXtJuVrR2JbPZzP1PzMBkNhMYFOyx55ZW17Mqq5RByWH4Wjp8XhUREfEoffIeRFBIaItAVJiX45HnllXXs2pXKfVqMRIREfEohaLDMAyDj994kSsnjGD5wvkeeWZ5TT0rM0sUjERERDxIoegwTCYTWTu2Ul9fx/T/u4W8PVkeeW5FbQMrMkuoa1AwEhER8YQOH4p+P/vsQKb+6W907zuA8tJiHrvrBurqPLNnWeXeYGRvcHjkeSIiIh2ZVrTeq2lFzNm/ZBxwcHVu9i5uuWACFWUlnHnhldw57R8eqy3Q5sOQTmHYfCwee6aIiMjxQCtae0FcYgoPPpmOyWTifx++xXezZ3rs2VX2xhaj2nq1GImIiLiLQtHv1B2iqyrtxFO54pZ7AHj2r/eTk5XpqbKotjtYqWAkIiLiNgpFvzP9601szq046PHLbrqbUadM4Lo7/0RcUooHK4PqOodajERERNxEY4r2auqTTL7zQ8y2AE7pGc15Q5Lw891/HI9hGF7d1d7famFISjj+Vo0xEhGRjk1jitzohK6RAMzbXMCjX6xnw57y/c7ZNxBVVVbw8/dfeaw+gJq9LUY1dWoxEhERcRWFot+5ZHgK95zeg6ggK4WVdfzr+y28tXgn1XUN+51bWV7GbRdP4rG7rmfl4h89WmdtvYPlmcUHrEtERESOXIcPRQdap6h3fAiPTu7LqT1jAPgxo5Bpn6/n1+yyFtcGhYTSd3AaTqeTx++7ifycbI/Wbq93siKzhCq7gpGIiMix0piivQ62TtGWvAr+u2gnBRWNCzae0DWSi4YlE2hr3EvXXlvDnZefxdaNv9JrwBCefnMWVqvNo7VbfcwM6RROkE37+4qISMeiMUUe1CM2mEcn9+H03rGYgEXbinjk8/WszioFwObnzyPPvEpQSCib1q7kxSenebzGugYnKzNLqFSLkYiIyFFTKGoFm4+Fi9KSeWBSL+JC/SirqWfGvK28/ON2KmrriU/qxANPzADgiw/e4PsvPvZ4jXUNjV1pFbX1Hn+2iIhIe6BQdAS6Rgcx7Q99mNg3DpMJlu4s5pHP17M8s5gRY0/nspvuAuCdF/5FQ73nw0n93mBUrmAkIiJyxDQI5Qj5WsycPzSJoZ3CeWPRTrJLa3hxwXaGdirhomvuwOl0cs5l1+Pj6+uV+hocBiszSxicEk6ov3dqEBEROR6ppegodY4K5M9n9uYPA+KxmEysyCzhr//bRN/JNxIWGeXV2hocBqt2lVBWrRYjERGR1lIoOga+FjNTBiXy0Jm9SQ73p9LewCs/7yB93jZKq+v44ctP+fStl71SW4PDYGVWCaXVdV55voiIyPGmw3efpaenk56ejsNx9KtDp0QE8NCZvflmXS5frM1h9e5SVs9YTOab92E2m+ncow+DR45xYdWt43AYrMoqZXByGGEBVo8/X0RE5HiidYr2Otg6RUcqu6SG/y7awY7CKoq+eoaqdXMJCY/kpU++Jyo23oUVt57FbGJgchgRgQpGIiLSvmidojYsMdyfByf15vyhScRMvAXf6FTKS4q47+arqauze6Umh9NgTVYpRZXeeb6IiMjxQKHIDSxmE5P6xfOXc4cw9Lq/YbIFsnvzGm6+/S4KvRRMHE6DNbsVjERERA5GociN4kP9+ctlp3HOnY8DsOunT7lrejo/bMrH6YVeS6cT1uwu9VowExERacsUitzMbDZx85UXcdZVUwGoKczmvaW7+Od3m8krr/V4PU4nrN1dSn6F558tIiLSlmmg9V6uGmh9MA6Hg3Urf6E4uBsfr9yNvcGJ1WJmyuAExvWKxWw2ufyZh2I2Q7/EUGKC/Tz6XBEREVfSQOvjkMViYWDaCZzSK4a/nNWXntF+2BscfLh8N09+u4mcshqP1uN0wq+7y7zSWiUiItIWKRR5gbO8gMw376NH/nz8fS1sK6jisf9tZF12mUfrMAxYl11GbpmCkYiIiEKRFyxfOJ/Nv67ih7ef4aKkCnrFBVPX4OT5H7ayZHuRR2sxDFi/p8zjLVUiIiJtjUKRF0w6/zJOm3w+ToeD5x6ayhX9gxmeGoHDMHj15x3M2ZDn0XoMAzbsKSe7VMFIREQ6rg4fitLT0+nTpw9paWkee6bJZOLOaf+gc4/elBQVMP2+m7h6ZBLjescAMHN5Fp+s3I0nx8AbBmzcU05WcbXHnikiItKWaPbZXu6efXYguzO3M/XCCVRXVnDO5Tdw8wN/5et1uXy6KhuAMd2iuGJkJywenpnWIzaYlMgAjz5TRETkaGj2WTuR1KkL909/HoBZ77zCVx+9wxn947lyVCdMJvh5ayEvzN9GXYPTo3Vtyasgs6jKo88UERHxNoUiLzvh1Ilcd+dDJHfpzoix4wA4qXs0t4ztio/ZxOrdpfz7+y1U1zV4tK6MvEp2FCoYiYhIx6Hus7280X22r5rqKvwDAlu8tiWvgud/2EpNvYPEMH/uGtedsADP7nTfOTqQrtFBHn2miIhIa6n7rB3aNxD98OWnzP96Nj1ig/m/iT0J9fclu7SGJ77ZRK6HF1vcUVDF1vxKjz5TRETEGxSK2phfVyzhiQem8sQDt7J43nckhwfwwMRexATbKKys48lvNrHTw+N9dhZWkZFX4dFnioiIeJpCURvTZ1Aap/7hPBwNDTx21/WsWLSA6GAbD0zsRUpEABW1DTz17WY27Cn3aF2ZRdVszlUwEhGR9kuhqI2xWCzc97dnOPH0M6mvr2PabVezdvliQvx9uW98T3rFBWNvcPLsDxks21ns0dqyiqvZmOPZMCYiIuIpCkVtkMXHhwf/8QLDTzoNe20Nf775cjatXYm/1cIdp3VnWKdwHE6Dl3/czrxN+R6tLbukhg17yj26sKSIiIgnKBS1Ub5WK4/8+1UGjRhDTXUVD/7xEvL37MbXYubGE7twSs9oDODdpbv4bHW2R0PKntIa1isYiYhIO6NQ1IbZ/Pz56/Nv0mdQGhPPvYTo+EQAzGYTlw5P4ayBCQB8sTaHd37ZhdPpuZCSW1arYCQiIu2Kj7cLkEPzDwzkH699iNXmh8n023YfJpOJswYmEOznw3u/7GLBlgIq7Q1cP6YzvhbPZN3cslqchkG/hFDMHt6KRERExNU6fEuRNzaEPVI2P//mQFRnr+XZv/4f+TmN+6Od0jOGP57UBR+ziRWZJTw7N4OaOofHassvt/NrdplHW6lERETcQSta7+XtFa1b65lH7+N/H71NYqcu/OvN2URExwCwMaecGfO2Ym9wkhIRwB2ndSfU39djdUUGWRmYFKYWIxER8SitaN2BXfrHO4lNSCI7czv3X38BZSVFAPSOD+G+CT0J9vNhV3E1T3yziYIKu8fqKqqsY/XuUhxqMRIRkeOUQtFxJiY+kX+89jGRMXHs3LqZB2+8hMryMgBSIwN5YGIvooKsFFTYeeKbTWQVV3ustuLKOlZnlSgYiYjIcUmh6DiUkJLKP177iLCISDI2rOWhmy+jpqpx64/YED8emNiLpHB/ymrq+ce3mz26EnVJVT2rdpXQ4HB67JkiIiKuoFB0nErp0p0nXvmQ4JAwNqxezuP/d3PzsbAAK/83oSc9YoOoqXfw7++3sHJXicdqK62uZ1VWKfUKRiIichxRKDqOde3Vl8dfeo/ouEQuueH2FscCrD7ceVoPBieH0eA0eGHBNn7MKPBYbWXV9azapWAkIiLHD80+2+t4mX12IHV1dqxW2wGPOZwGby/J5OethQBMGZTAmf3jW6x55E5Bfj4MSQnH6qP8LSIirqfZZ9LCvoEoY8Na0h9/CKezsYXGYjZx1ahOnNE/DoDZq/fw/rIsnB7KwpW1DazILMHe4Lm1k0RERI6GVrRuR6qrKnnwxospKymmvq6OO6b9A5PJhMlk4tzBSYT4+fLBsix+2JRPZW0D145OxccDq19X2RuD0ZCUcPx8LW5/noiIyNFQS1E7EhAYxG1/no7ZbOZ/H73NC0880mJvsnG9Y7lhTGcsJhNLdxbz/A9bqa33TAtOtd3ByswSjz1PRETkSCkUtTNjJ57NPY/9G4BZ77zCf597osXxEV0iue3Ublh9zKzPKefpOVuoqK33SG3VdQ5WKBiJiEgbpVDUDo2fchG3/Xk6AO+//Czvvfxsi+P9EkO59/QeBNl82FFYxRPfbCK/otYjtdXUOVi+s8Sj+7OJiIi0hkJRO3XWJddw472PAPDfZ6cz/+vZLY53iQ7i/ok9iQi0kldu5/GvNpGR75lFHmvrHSzPLKa6rsEjzxMREWkNhaJ27IJrbuHKqfcyZNRJjBh7+n7H40P9+dOkXnSKDKDS3sDT323hl+1FHqnNXu9kRWaJx7ruREREDkfrFO11PK9TdCiGYeBoaMDH1/eg59jrHby6cAerdpUCcPbABP4wwDNrGVnMJnrHhxAX6uf2Z4mISPujdYqk1UwmU3MgMgyDt//zT3767ssW59h8Ldw8tisT+sQC8NmaPby+cKdHVqN2OA3WZZexJa8C5XMREfEmrVPUgfz43Re8lf5PfHx8edTPjxEnjWs+ZjaZuGBYMjEhfrz7SyaLtxdRWGln6sndCPJz/6/JrqJqKmob6J8YqtWvRUTEK/Tp04GMGXcmJ0+aQkNDPX+54zpWLfl5v3PG9ojmjtO64+9rISO/kse/3khuuWdmppVU1bF0RzHlGmckIiJe0OFDUXp6On369CEtLc3bpbidxWLh/unPc8KpE6mvs/PIrVeyftWy/c7rmxDKA5N6ERloJb/CzvSvNrI513Mz01bsLCGnrMYjzxMREWmigdZ7tdeB1gdSV2fnkalXsWLRfAKCgnnq9Y/p0XfgfueV1dSTPm8r2wurmvdQO6FrlMfqTIkMoHtMkMc2rxURkeOPBlrLMbFabTz63Ov0HzaS6soKHrrpUqqrKvc7L9Tfl3vH92RYp3AcToPXF+5k9upsjw2I3lVUzcpdpdQ1uH/At4iIiEJRB+XnH8Df/vMOfQalceXU+wgIDDrgeVYfMzee1IUz+sUB8OXaHF75aYdHZqaBxhmJiIjnqPtsr47Ufbavhvr6FmsYlRYXEhoeecAuq58zCnl7SSYOw6BrdCC3ntKNYL+Dr3/kSmYz9I4PIT7U3yPPExGR44O6z8Rl9g1E5aXF3Hn5WTx+383U2fefcTamexR3nd6dAKuFbQVV/P2rjewp9cyAaKcT1meXszlX6xmJiIh7KBRJsw2rl5ObvYv5X8/mvmsvoLS4cL9zesWF8OCkXkQH2SisrGP615vYmFPusRqziqtZuatE44xERMTlFIqk2ciTxzP9pQ8ICgllw+pl3H7pmezanrHfefGh/vzpjF50jQ6kpt7BM99n8HPG/gHKXUqq6jXOSEREXE6hSFoYPHIMz777JXFJKeRkZXLHZX844CKPwX6NM9OGp0bgMAzeWLyTT1buxumhrq3aegfLdxZ7rPtORETaP4Ui2U9Kl+48//5X9BmURmV5GQ/+8WIW/fDtfuf5WszccGJnJg+IB+Drdbm89ON2j3VtOZ2wYY/GGYmIiGsoFMkBhUVE8dTrH3HypClExyXQe+CQA55nMpk4e1Ai145OxWI2sSKzhKe+20xZjee6tprGGdkbHB57poiItD+akr9XR52SfzhOp5PSokIiomOaX3M0NGDx2X+T2C15FaTP20pVnYPIQCu3n9adxDDPTaG3+ZoZkBRGqL9nlgkQERHv05R88Riz2dwiEH0z633uuvJsSgoL9ju3R2wwfzqjN7HBNoqq6nji602s31PmsVrt9U5WZGqckYiIHB2FImm1muoqXv/342xcs4LbLz2DzK2b9zsnNsSPByf1pkdsEDX1Dp6dm8GCLfsHKHdpGme0Kbccp1ONoCIi0noKRdJq/gGB/Ovtz0hI6Uxudha3X/YHVixasN95QX4+3DWuB6O6ROI04O0lmXy4PMujIWV3cY3GGYmIyBFRKJIjktSpC8+//z/6D23cTPZPN13KVx+9s995vhYz145O5exBCQB8tyGPFxZsw17vuZBSWt24npEnB32LiMjxS6FIjlhIWARPvDqT0yafj9Ph4N+P3ssr//zrftPiTSYTkwckcMOJnfExm1iVVco/vttMaXWdx2ptGmeUrXFGIiJyGApFbUBogC+DUsKIC/XDYt5/I9a2yGq1cf/057ly6n0A+AcGHXATWYARnSO5Z3wPgmw+ZBZV8/evNpJVXO2xWp1O2KhxRiIichiakr+Xt6bkm0wwvHNE827zDQ4n+RV2cspqKa2u43j46axbuZS+g9MOGoqa5FfU8tzcreSW12LzMfPHk7owICnMM0XuFRbgS/+kUGw+Fo8+V0RE3ENT8tuRTpEBzYEIwMdiJiHMn6GdwhndLYpuMUEE2vZfE6gt6TdkeHMgqqmu4tE7rmXHlo37nRcT7MeDk3rRKy4Ye4OT5+dt5fuNeR5djbp5nFG1xhmJiEhLCkVeFGC10Dkq6KDH/XwtpEYFMqprJCO6RJASGYDVp23/yP777BMs/P4r7rx8Mst++mG/44E2H+48rTtjukVhGPDBsiyenrOF/Ipaj9Vor3eyYlexR7vwRESk7Wvbn7DtXM+44FaPIQr286VHbDAndo9icBsef3T5zXczMO0Eqqsq+fMtl/PFB2/sd46PxcxVozpxcVoyVouZTbkVPPr5Br5dn4vDQ2N+nE7YnFvBmqxS6h2e2atNRETaNo0p2svTY4riQv3olxh6TPdoGn+UW15LSVXbGX9UX1fHvx+9lzmffQjAeVf+kRvufQSLZf9xPPkVtby9OJONuRVAY3fiVaNSSYkI8Fi9fr4W+iWGEBZg9dgzRUTENVw5pkihaC9PhiJfHzOjukS6tCustt5BXnktOWW1VNY2uOy+R8swDN57+VneeO4JAEadMoEHn/wP/oGBBzx34bYiPlyeRXWdA7MJJvaNY/LABHwtnmnMNJmgS3QQqZEBhx0wLiIibYcGWh/nuscEuXxskJ+vhU6RgYzs0jbGH5lMJi7745089M8X8bXa2Lh2BeVlxQc9d0y3KB47ux9DO4XjNOCrdbk8+sV6tuRVeKRew4Bt+ZWs3FVKrQcXmBQRkbZDLUV7eaqlKDzQytBO4W67/74Mw6C4qo6csloKKuweG6/zextWL8dkNtN7wJBWnb9yVwnv/rKreSXqsT2iOX9IEv5Wz0yj9/Ux0zchhKggm0eeJyIiR0/dZ27giVBkNsPILpEEWD0/xb7B4aSgsnH9I2+PP1o87zvMZhMjxp5+0HOq6xr4eMVufswoBCA8wJfLRnRiUHKYh6psHN/UNToIcxsc0C4iIo3UfXYI1dXVdOrUiXvvvdfbpeync1SQVwIRNM74ig/1Z0hK4/pH3WODCPLzfC2Z27bw+H1/5JFbr2L2u68e9LwAqw9Xjkrl3vE9iAm2UVJdz4x5W3lxwTaP7WWWWVTN8swSaurUnSYi0hG0u1D097//nZEjR3q7jP0E2nzo5MEZVYfy+/FHMSGe6yZKTOnMyWecg9PpJP3xP/PcYw9QV2c/6Pm94kJ4dHJfJvaNw2yC5ZklPPLZOhZuK/TIoo/lNfUs2VFEbpnn1lESERHvaFehKCMjg02bNjFp0iRvl7Kf3vHBbbIbJtjPlwFJYQxMDsPP1/1jdnx8fbn7L09zwz0PA/DFB29w+yVnsmt7xkGvsfqYOX9oEg+d0ZvkcH+q6hz8d+FOnvk+g8LKgwcqV3E4DNZll7FhT7nXxmWJiIj7tZlQ9OOPPzJ58mQSEhIwmUzMnj17v3PS09NJTU3Fz8+PESNGsHTp0hbH7733XqZPn+6hilsvMdy/za+BEx1sY+TeWWvunpFuMpm48NqpPJb+FiFhEWzbtI5bLhjPVx+/e8jrOkUG8tCZvTlvSCK+FhPrc8p55PP1zNmQ55GNXveU1vDLjiIqarVFiIhIe9RmQlFVVRUDBw4kPT39gMdnzpzJ3XffzbRp01i5ciUDBw5kwoQJ5OfnA/DZZ5/Ro0cPevTo4cmyD8vma6ZbzMG38mhLfCxmesQGk9Y5gmAPjDcaefJ4Xp49jyEnjMVeW0Nu9q7D12g2M6lfPI9O7kuP2CDqGpzMXJ7FE99sIrukxu01V9sdLNupLUJERNqjNjn7zGQyMWvWLKZMmdL82ogRI0hLS2PGjBkAOJ1OkpOTue2223jggQd48MEHeeedd7BYLFRWVlJfX88999zDI488csBn2O127Pbful7Ky8tJTk52+eyzAUmhxIT4uex+nmIYBlnFNWwrqHR7l5HT6WTuFx9zyhnn4OPbuDlufV0dvtZDt645DYOfMgr5eMVuauodWMwmzugXxxn94z2y6GNMiI3e8SEeW2BSRET21+Fmn9XV1bFixQrGjRvX/JrZbGbcuHEsXrwYgOnTp5OVlcXOnTv55z//yQ033HDQQNR0fmhoaPNXcnKyy+uOCrYdl4EIGoNpSmQAo7pGEhXs3oHYZrOZ08++sEUguvvKs3n1X3+jvq7u4NeZTIztEc1fz+7LoOQwHE6DL9bm8NcvN7A1v9KtNQPkl9v5ZXsxZdXqThMRaQ+Oi1BUWFiIw+EgNja2xeuxsbHk5uYe1T0ffPBBysrKmr+ysrJcUWozi8VErzj376Hmbn6+FgYlh9E/KdRjK2QvWTCHTb+uYuZrM7jj8j+we+e2Q54fHmBl6slduemkLgT7+ZBTVsuT32zivV92uX116tp6B8szi9lRWOXW54iIiPsd9cCRtLS0A+4RZRgGJpNpv0HQnnT11Vcf9hybzYbN5r4WkK5RQR6ZzeUpsSF+RARa2Zpf6faxOyeefibTnn2dfz1yDxnr13Lz+eOY+qe/M+GcSw66L5nJZGJYagS94kP4cHkWi7YV8cPmfFbvLuWKkZ3of4yb7x5K0xYhJdV19E0IwebTfn7uIiIdyVGHoo8//tiVdRxSVFQUFouFvLy8Fq/n5eURFxfnsTpaK8Tfl+QIf2+X4XK+FjO940OID/VjY04FVXb3bTw7ZtwZ9Oo/mCcfuJXVSxfy9MN3s/SnH7jr0X8SHBp20OuCbD5cO7ozIztH8taSnRRW1vHs3AxGdI7g4rRkgv183VZzcWUdv2wvpm9CCJHaIkRE5Lhz1P0hnTp1OuSXK1mtVoYOHcrcuXObX3M6ncydO5dRo0a59FnHymRqXJOoPe+0HhZgZUTnCLrGBGF2Y49aVGw8T7z6Idfd9RAWHx9++u5Lnn7k7lZd2ychhL9M7svpfWIxmeCXHcU8/Nl6ftle5NZFH+sanKzaVcrW/AqPLC4pIiKu02a6zyorK9m6dWvz9zt27GD16tVERESQkpLC3XffzVVXXcWwYcMYPnw4zzzzDFVVVVxzzTVH+xaAxrWP0tPTcThcM/YkJSLAra0RbYXZbKJzVCCxITY25lRQUnXwAdHHwmKxcPH1tzF45Ik8+5f/44a7H271tTZfCxcNSyYtNZw3F2WSXVrDKz/vYMmOYi4fkeLW1pydhdWUVNfTLyHUYxvZiojIsTnqKfmZmZmHPH6krUXz58/nlFNO2e/1q666ijfeeAOAGTNm8NRTT5Gbm8ugQYN47rnnGDFixBE952BcsSGsv9XCyC6RWNrgytXullNWw5a8SuobnG57RlPgbvLlh28xeMSJJHbqfNhrGxxOvlmfy5drc2hwGth8zJw3JImTe0ZjdmOrno/FRO/4EGKP01mIIiJtnSun5LtknaKvv/6arVu3ctttt5Gbm0tJSQm9e/c+1tt6lCtC0eCUsA49lqSuwUlGfgU5pe7fJ2zVkp+5//oL8PMP4NaHHuf0sy9sVZflntIa3lqcydaCxin7veKCufqEVKLc/HNLDPenR2xwhwzMIiLu1KbWKbr33nv54IMPmleitlgsrZr91d7Ehfp16EAEjXuU9U0IZUincALc3GWUmJJKv6EjqKmu4qmH7uDx+26msrzssNclhPnzfxN7cunwFKw+ZjblVjDt8/X8uKXArWOAsktqWLqjmEo3Dk4XEZFjc8yhaO7cubz55pv4+zfOtoqOjqa2tmPtKO5jMdE99vjYysMTIgKtjOwSSWpUoNsGYsckJPHU659wzR0PYrZYmP/1bG467zTWrfjlsNeaTSZO7RXDtD/0oXtMEPYGJ28tyeTZuRmUVLtnbBRAlb2BZTuK2V2iLUJERNqiY/7I8vX1xel0NnddFBcXY3bnlCQXS09Pp0+fPqSlpR31PbrHBmttmt8xm010iwlieOdIwgLcM/DcYrFw6Y138Mw7XxCf3Im8Pbu55+pz+OCV51t1fWyIH/eN78mFw5LwMZtYt6ecaZ+vZ/E2981QczgNNuVUsHJXiVuXNBARkSN3zOnl9ttv56KLLqKwsJDHHnuMk046if/7v/9zRW0eMXXqVDZs2MCyZcuO6vrwQF8Sw9rfmkSuEmTz2buoYjA+FveMp+k9YAgvfjKX08++EKfTSUR0TKuvNZtNjO8TxyN/6ENqZADVdQ5eW7iD9PnbKKtx3/YdxZV1/LKjiIy8Choc7hucLiIirdfqgdYbNmzg/fff55577iEsLKzFsU2bNjF37lycTiennXYaffr0cUetbnU0A63NZhjROZJAm/t3lG8P7A0OtuRWklfuvu7VdSt+oe+Q4c0tlyVFBYRHRrfqWofT4Jv1uXy+Zg8Op0GQzYfLRqSQlhrhtnqhcSxW99gg4kMVrkVEjpRXZp9dccUVVFZWMmvWrP2O1dbWsnPnTnr16nVMxXjT0YSiLtGBdInWWKIjVVhpZ1NOhdv3JSsvLeaP557GwOGjue3P01v9c80qqeb1n3eQtXc7k2GdwrlsRIrb158KC/ClR1wwIR1gnSsREVfxyuyzJUuWcPvttx/wmJ+fHzfccAPTp08/pmKOJwE2C6mRgd4u47gUFWRjVNdIOkUG4M6Fv1f/spDiwnzmfvExN513GhtWL2/VdcnhATx0Rm/+MCAeswmWZ5Yw7fP1rNpV4r5igdLqepbtKGZjTjl1blzvSUREDqzVoWj37t1069btoMdvuukmPv/8c5cUdTzoEx+CWWvOHDWL2UT32GBGdY0kJsQ9SxmcNGEy/3pzNnGJyeTu3sVdV57NOy/8q1Wrl/tYzEwZlMifJvUmIdSP8toG0udv47Wfd1Bd574B0obROH1/0bZCsoqrtVWIiIgHtToURUREkJOTc9Djw4cPb7FNR3uWEOZPWIDV22W0CwFWHwYkhZHWOYLwQNf/m/YdnMaLn8zltD+ch9Ph4M0Z/+Deq88lb09Wq65PjQrk4T/0YWLfOEwmWLy9iGmfr2dd9uHXRDoWDQ6DzbkV/LKjmFI3LhMgIiK/aXUoOumkk5q32zjgjczm43J9oiOdkt80KFZcK9Tfl6GdwhmcEkawn2sHrgcGh/DAk+k88GQ6AYFBrFv5C++88K9WX+9rMXP+0CQemNiL2GAbJdX1PDM3g7cW73T7uKjK2gaW7yxhXXaZ258lItLRtXqg9YoVKxg5ciTp6enceOON+x1/5513eOqpp1izZo3Li/SE1g607p8Uqn2sPCC3rJbtBZVU17k2CORkZfLqv//GXY/+k6CQUGD/PdUOxd7g4NOV2czdlA9AVJCVq09IpVfcsQ3uaw2LxUSXqECSwwPUdSsispfX9j575ZVXuPnmmznllFOYOnUqQ4YMISgoiJ9++olbbrmF22+/nfvvv/+YCvKW1oSiyCArg1PCPVxZx+V0GmSX1rCjsMptA48Nw2DabVfTvU9/LrruVqy21gXeTbnlvLFoJ4WVjV1bp/aK4bzBidh83b+IZ4DNQs/Y4A6/rYyICHh5Q9iff/6Zu+++m+XLlzf/v2vDMBg/fjyff/45VuvxOdbmcKHIYjYxqmskfh740JOWHE6DXcXV7CyqwuFw7cDjVUt+5v+uOx+AhJTO3P7wEww9YWyrrq2td/DRit0s2FIAQEywjWtHd6ZbjGe6V6ODbfSIDcbfzfvMiYi0ZV4NRU02bdrEypUrqa6upl+/fowcOfKYCvG2w4WiHrHBpEQGeKEyaVLX4GRnURW7S6pxuqjhyDAMfvz2c/7zxCMUF+QBMHbiWdx0/1+Jiolr1T3WZZfx5uKdlFTXYwLG94llyuBEfC3u3+7GbIZOkYGkRgZiUZeaiHRAbSIUtTeHCkXBfj4M7xzR6nEn4l619Q62FVSSW1aLq357qyoreGvGU8x+91WcTicBgUFcdev/cfal12LxOfzA7+q6Bj5YlsWibUUAxIf6ce3oznSO8sxaVn6+FnrEBhGj8W4i0sEoFLnBwUKRyQRpnSO0ynAbVGlvYFt+JQUVdpfdc+vGdTz32P1sXLOC5M7dePHTuVitrR+7szqrlLcW76S8tgGzCc7oF88fBsTj44FWI4CIICs9Y4O19YyIdBgKRS6Unp5Oeno6DoeDLVu27BeKUiID6BHbuu0hxDtKq+vYml9JabVrNnB1Op1888l7JHbqzMDhowFoqK+nuqqSkLDDD7SvrG3gvaW7WLqzGIDkcH+uHd2Z5AjPdL+azZAUHkCXqECPhTEREW9RKHKDA7UU+flaGNU1UmM1jhMFFXa25ldSZXf9itMfvfECM199nhvueYTxUy5qVVfq8p3FvPPLLirtDVjMJiYPiGdSv3iP/T5Zfcx0iwkiIUwbzYpI++WVvc86ol7xwQpEx5HoYBsju0TQNzHEpbMEDcPgp+++pKykmH/++U7uueocdm7ddNjrhqVG8Jez+jI4OQyH02D26j1M/3oje0prXFbbodQ1ONmwp5xlO4spr3VNK5qISHumlqK9ft9SFBviR/+kUG+XJUfJ6TTYXVLDjqIq6l2wxlFDfT2fvv0yb//nn9TW1GDx8eG8K27k8pvvwT/w0IOpDcNgyY5i3l+6i+o6Bz5mE2cPSmBc71iPzFCDxrFxCWH+dI0Owuqj/y8kIu2Hus/cYN9QFBoawqiukdh8tP7L8a7B4SSzuJpdRdU4nMf+q56/Zzf/efIRFn7/FQDRcYnc/8TzDEw74bDXllTX8ebinazLLgcaV8OeMiiR4Z0jMHtoZqOPxUSXqCCSwv21KraItAsKRW6wbyga2j2BpHCtSdSe2Bsc7CisYk9pjUvWOFoy/zvSH3+IgrwcXvpkLp269WzVdYZhsGh7EbNWZlNa09illRIRwHlDEumb4LmWSZuvmS7RQSSE+mmpCRE5rikUuUHTP+q8tTs4uX+qt8sRN6mua2B7QRW5Zce+eXFtTTW/rviFtDGnNL+2YtECBgwbhe9hVna3NziYuzGfr9flUrN3o9c+8SGcPyTJo4uEBlgtdI0J0n5+InLcUihyg6Z/1D35RcRHR3i7HHGz8tp6NuVUUF7jugHIWzf+ytQLJ5CY2pXb/zydQSPGHPaaitp6/vdrDvM2FzR3743oHME5gxOJ8uDeZsF+PnSNCfLoM0VEXEGzz1woPT2dPn36kJaWBqBF7zqIED9f0lLDSY0KxFW9RyWFBYSER5C1PYP7rj2fJ+6fSnFB/iGvCfbz5eK0FP52dj+GpzaG8V92FPPn2euYuSyLylrXLy9wIBW1DazeVcqKzGJKq+s88kwRkbZGLUV7uTJpyvGlpKqO9XvKqd3bjXUsKspK+e9zT/DlzDcxDIPA4BCuuf0B/nDRVVgshx+4v7Ooik9W7GZjbgUA/r4Wzugfx2m9Yj06aywq2EbX6ECCtZK7iLRx6j5zA4Wijq3e4WRTTgV55cc+1ghg86+rePax+8lYvxaAIaNO4olXZrZqULNhGKzfU84nK3eTVdK4plF4gC9nD0rkhC6RHp01FhfqR5foQAKsakEVkbZJocgNFIoEYE9pDZvzKnA4jv0/C4fDwZcz3+T1Z6dz3Z1/4qxLrjmi652GwS/bi5m1OpviqsYurcQwf84bkkj/xFCPzRprWuOoc1SgSxfFFBFxBYUiN1AokibVdQ2syy532SDsprFGTd1nKxYtoLgwn9P+cB5m8+G7xOodTn7YlM//fs2huq6xi69HbBDnD02iS1SQS2psDYvZRFK4P6lRgR5bdFJE5HAUitxAoUj2ZRgG2wqqyCyqwpX/hdTWVHPD2WPJzc6iS8++XHfXn0gbc2qrWn2q7A18tS6HuRvzadg7U21Yp3DOGZzo0Sn1PhYTnSIDSYkI0DY4IuJ1CkVuoFAkB1JaXce6bNcMwgaor6vjkzdf4v1Xn6O6snEw9YC0UVx315/pM3Boq+5RVGnnszV7WLytCAOwmEyM7RHNHwbEE+LvuYHRVh8znaMCSQzT6tgi4j0KRW6gUCQH4+pB2ADlpcW8/8pzfPbef6mvswMw+rRJXHfXQyR37taqe2SVVPPJyt3N24bYfMxM7BfH+N6x2Dw49sfP10KX6EDitTq2iHiBQpEbKBTJ4bhyEHaT/Jxs3v7PP/lu9kycTidPvvYRQ0aeeET32JRbzscrdrOzqBqAUH9fzhqYwJhuUR7t3gq0+dA1OpAYrY4tIh6kUOQGCkXSGtV1DazfU05ZtetWwgbI3LaFn777kstvvrv5tRWLFtC9T39Cwg6/wrrTMFi+s4RZq7IpqGxseYoL8ePcIYkMTg7zaAtOiL8vXaMDidTq2CLiAQpFbqBQJK3lrkHY+yopKuCqiSMxmc1cdN2tnHP59fgHBB72ugaHk/lbCvhybQ6V9sbVsLtGB3L+0CS6xwS7p9iDCA+00i0miFAPjnMSkY5HociF0tPTSU9Px+FwsGXLFoUiaTVXD8Le146MjTxx/61s37wegIioGC6/+W4mnXcZPr6HDxk1dQ6+WZ/LnA151DmcAAxODuPcIYnEh/q7vN5DiQ620TUmiCBtoSMibqBQ5AZqKZKj4Y5B2E2cTifzv5rNGzOeJCcrE4CE5FSuvv0Bxk48q1VrHJVW1/H5mj38tLUQwwCzCcZ0i+LsQYkebcExmSA2xI9uMUFaAFJEXEqhyA0UiuRY5JTVsCnXtYOwm9TX1fHVx+/wzov/orSoEB8fX974ehGxCcmtvsee0ho+XZnN6t2lAPj5mjmjXzyn94n16EKMZjMkhwdoAUgRcRmFIjdQKJJj5a5B2E1qqqr45K2XsNtrue7OPzW/vmfXThJSUlt1jy15FXy4PKt5plpkoJXzhiSRlhru0cHYPhYTnaMCSQ4P0BpHInJMFIrcQKFIXMEwDLYXVrGz0H2DsPe1deOv3HLBeEafdgbX3PEAKV26H/aapj3VPl21m5K9Aa5rdCAXpSV7dNsQaFzjqGtMoMfHOYlI+6FQ5AYKReJK7hyEva9P336Fl/4xDafTidlsZvw5F3PlLfcSHZdw2GvtDQ6+W5/H1+tzqWtoHIw9onME5w5O9Ph0+mA/H7rHBhMRaPXoc0Xk+KdQ5AYKReJq9Q4nm3MryC1z/SDsfe3cuon/PvsEi374BgCrzY8pl17LRdffRkhY+GGvL6muY9aq7OZtQ3wtJib0iWNivziPD4qOCLLSPSaIYD9N4xeR1lEocgOFInEXdw7C3tf6Vct47d9/59cVSwBISu3Ka1/81KpZagCZRVXMXJ7FlrxKoHFl7CmDEhjdNcqj435MJogL9aNrtGaqicjhKRS5gUKRuFNNnYN1e8rcNgi7iWEYLP1pLq//+3EmnX8pUy67Hmic3u90OA67xpFhGKzKKuWjFbspqGhcGTs53J+L0pLpFefZ/y7MZkiJCKBTpGaqicjBKRS5gUKRuJsnB2E7nU4MpxOLT+OCiQu++Yz/Pvck19x+PyeOn3zY1qN6h5MfNuXz5docavaOixqUFMb5w5KI8/DeZr4+ZjpHBpIU7q+ZaiKyH4UiN1AoEk8pra4jI7/S7a1G+7rz8rNYv2opAN169+fSG+9g9LgzDhuOKmrr+WJNDvO35OM0wGIycUqvaP4wIMHjK1T7Wy10jQ4iLlQbzorIbxSK3EChSDytyt5ATlkNOWW12Oudbn1WdVUln7z5Ih/99wVqqqsASO7SnYuvu5VTzzz3sN1qe0pr+GjFbn7NLgMg0Gph8sAETu4ZjU8rxyy5Soi/L91jggjXTDURQaHILRSKxFsMw6Cwso6cshoKK+043ZiPykqKmP3ua8x+9zUqyxsDztiJZ/Hnp19u1fXr95Tx4fLdZJfWABAbYuOCockMTAr16OKPAFHBNrppTzWRDk+hyIW0Iay0JXUNTvLKa8kuraGytsFtz6mqrOB/H77Fx2+8yH1/f5a0E09tfh3DIDD44P8NOJwGP28tZPbqbCr21tg7LpgLhyWTHBHgtpoPxGSC+FB/ukQHaqaaSAelUOQGaimStqaitp6cslpyymqpb3BP81GdvRZfq625lefNGf/g07df4exLr+Wcy68nPDL6oNfW1Dn43685fL8xjwangYnGzWanDPbsZrMAFrOJ5IgAUiMD8NFMNZEORaHIDRSKpK1yOg0KK+3sKaulqNLutplrhmFw95VTWLfyFwBsfv5MOu9SLrj6ZmISkg56XUGFnU9W7mZ5ZknjdT5mzugfz+m9Y7H6eDag+PqY6RIVSGKYZqqJdBQKRW6gUCTHA3uDg9yyxu61arvrtxBxOp0snvct77/8LJvXrQbA4uPDuMnnc9F1t5LcudtBr92aX8nM5VnsKGwcyB0RaOW8IYkMT43w+HijAKuFbjFBxHh4+QAR8TyFIjdQKJLjTVl1PXvKasgrr6XBxatlG4bBqiU/8f7Lz7J66UIAJp53Kff89V+HvM5pGCzdUcwnK3/bbLZLVONms12jPbvZLEBoQONMtbAAzVQTaa8UitxAoUiOV06nQX6FnT1lNZRU1bm8e23DmhV88Mpz3HjfNJI6dQEga8dWSosL6T905AGvsTc4mLMhj6/X5WLfOx5qeGoE5w3x/GazANF7Z6oFaqaaSLujUOQGCkXSHtTWOxoHZ5fWUF3n+u61Jn+/94/M//oz+g0ZwSU33k7amFMP2EVWunez2UV7N5u1mE2M7BzBhL5xJIT5u62+A9GeaiLtk0KRGygUSXtTUlXHnrIa8svtOJyu+8/cMAyef+wBvvn0ferr6wDo2qsfl9x4O2PGnYnFsn/g2FVUzYcrstiUW9H82uDkMCb2i/N4t5rZDIlhAaRGBWDzUTgSOd4pFLmBQpG0Vw0OJ/kVdnLKaiipct3WIoX5uXzy5ot8OfMtamuqAUhK7crVt9/P2AlnHfCabQWVfLMul1VZpc2v9YgN4ox+8fRNCPHogOymafydIgO04azIcUyhyA0UiqQjqK5rYHdJDdklNS5rPSovLW5cJfud16goL+XGex/hgmtuOeQ1e0pr+GZ9Lr9sL8ax909Qcrg/E/vFMaxTBBYPTqf3sZjoFBlISkSAR58rIq6hUOQGCkXSkdgbHOwqqma3C8NRdVUlX330DmdccDkBgY1dYkvmf8f2LRs56+KrCQoJ3e+a4qo65mzI48eMguYB2VFBVib0jWN01yiPrnNk9THTWWsciRx3FIrcQKFIOqK6Bie7iqvJKqnG4YZp/TeffzrbNq0jICiYsy65hnOvuOGAq2RX2huYtzmfuRvzqbQ3bh0S7OfDuN6xnNIzmgCr52aN+VstdIkOJC7Ez+PrK4nIkVMocgOFIunI6h1Osoqr2VVc7bI1j5xOJ/O+msUHrzzHzq2bAfC12jj1zHOYctl1dOvdf79r7A0OFm4t4tv1uRRVNQ7i9vM1M7Z7NKf3ifXoekOBNh+6xgQSE6wFIEXaMoUiN1AoEmkclJ1VUsOu4mqX7bfmdDpZMv873n/5WTb9uqr59fOu/CM33f+XA9fhdLJsZwnfrMslu7QGAB+ziVFdIpnQL444D65UHeLvS9foQK+sryQih6dQ5AYKRSK/cTgNdpdUk1lUTZ2LwpFhGGxcs4JZ77zKT3O+5JF/v8oJp04EoLK8DKfTQUhYxH7X/JpdxtfrcsnIrwTABAxJCWdivzg6RwW6pLbWCA+00i06iNAAz252KyKHplDkBgpFIvtzOA32lNaws6gKe71rwhFAYV4O4VExzWsavf2ffzLztXRO+8N5TLnsOjr36L3fNRn5FXyzLpc1u8uaX+sVF8ykfnH0iffcdP7oYBtdY4II0urYIm2CQpEbKBSJHJzTabCnrIadhdXU1rt+pew/3XQpy376ofn7QSPGcM7l1zNi7On7LQaZXVrDN+tyWbrjt+n8KREBTOoXx9CUcI/MHDOZIDakcXVsf6sWgBTxJoUiN1AoEjk8p9Mgp7yWzMIql24jYhgG61b+wqx3XmXh91/hdDa2SsUlpXDB1Tdz1iXX7HdNUaWd7zbk8dPWwuYuvphgGxP6xnFC10iPLMhoNkN8qD+dowK1dYiIlygUuYFCkUjrGYZBbnktOwqrqLa7tuUof89uPv/gDb76+F0qykqYdN5l3P3Xpw96fkVtPT9syueHTflU7Q1qIX4+nN4nlrE9PDOdv3F1bH86RQZqdWwRD1MocqH09HTS09NxOBxs2bJFoUjkCBiGQX6Fne0FVVTtXV/IVWprqvnhf7PoO2gYnbr1BCBjw1pe+/fjnHP5daSdeBpm828BxF7v4KethXy3Po/i6sbp/P6+Fsb2iGZc7xiPTOdvWh07OdwfH4UjEY9QKHIDtRSJHJv8ilp2FFRRUevacLSvf/zpduZ89iEACSmdmXLZdYyfchGBQcHN5zQ4nSzdUcw363LZU1YLNE7nP6FrJKf3iSU+1N9t9TWx+phJjQwkKVyrY4u4m0KRGygUibhGQYWdHYVVlNe4bvPZJjlZmXz2/ut88+n7VFWUAxAQGMT4cy7m7EuvJalTl+ZznYbB2t1lfL0uh20FVc2vD0gMZXzfWHrGBrt9xpqfr4XO0YEkhGp1bBF3UShyA4UiEdcqqmwMR6XVrg9HNVVVfP/FR8x69zWytmcAEB2XyDtzlrXoUoPGLr6M/Eq+25DHmqxSmv7gJYf7M75PHGmp4W7v6rL5mkkKDyAhzA+bjwZki7iSQpEbKBSJuEdJVR3bC6so2btthys5nU5WLv6R2e++St/Bw7nkhtsBcDQ08N1nH3LyxLPxD/xtgce88lq+35jHwq1F1DkaZ6yF+ftyaq8YxvaIJtDNaw+ZzRAT7EdSuL9HtywRac8UitxAoUjEvUqr69hZVE1RpR13/NUxDKO5i+qn777kr3ddT2BwCBPPvYSzL72W+KROzedW1jawIKOAHzblU7a3m8/qY2ZM1yjG9YnxyH5nwX4+JEUEEBfih0XjjkSOmkKRGygUiXhGdV0DWcU17CmrweGizWd/7+fvv+KVpx9jz64dAJhMJkaMPZ0zL7ictDGnYvFpbBGqdzhZurOYORvy2F3SuMeaCRiUEsb4PrF0iw5y+1ggH4uJxDB/ksIDtBCkyFFQKHIDhSIRz6p3ONlTWsPukhpqXLgQZBOn08myn39g9juvsXzhvObXI2PieHnWDy32WTMMg405FXy3MZd12eXNr3eOCmR8n1iGpIS7vTXHZIKIQCvJEQFEafNZkVZTKHIDhSIR7zAMg4IKO7uKq90yKBtg1/YMvv7kXeZ89iHxyak8//5Xzcc2rF5Ot979sNoau8z2lNYwZ0Mei7cX0eBs/PMYEWhlXO8YTuwW7ZHWnACrhaTwAOLD/LQYpMhhKBS5gUKRiPeV19azq6ia/IpanK7bf7ZZfV0dRQW5xCWmND6vtISLTxmEn38A4yafz6TzL6Vz98bNaMtq6pm/OZ95mwuo3LswpZ+vmRO7RzOuVwyRHmjNsZhNxIb4kRzhT7Cfr9ufJ3I8UihyA4Uikbajtt7B7pIasktrqG9wQzraa+PalTx21w0U5GY3v9ZrwBAmnXcZJ086m4DAIOoanCzZXsR3G/PI3bsYpNkEQzuFM75PHJ2jAg92e5cKC/AlKTyAmGCbFoQU2YdCkRsoFIm0PU0b0O4qqnb5NiJNHA4HKxct4KuP32Xx/G9xNOxtFfIP4OF/v8LwE09rrMUwWJddxpwNeWzMrWi+vntMEOP7xDIwKcwjYcXqYyYx3J/EMH9tQiuCQpFbKBSJtG1FlY3jjooqXb/eUZOSwgLmfP4hX3/yHjm7M3nv+5VERMcAjeOSwiIiCQmLYFdxNXM25LF0ZzGOveOOYoJtjOsdy+iukdg8EFZMJogOtpEcHkB4oNY8ko5LocgNFIpEjg/VdQ3sKq4mp7S2OZC4mmEY7Nq2pXkjWoD7rjmP9auWMXrcGUw6/zIGDR9NWW0D8zblM39LAdV7Z9AFWC2c3COaU3t5ZhNagECbD8kR/sSH+mvNI+lwFIrcQKFI5PhS73CSXdI4pb+23vVT+vdlr63hrivPJmP92ubX4pJSmHjupUyYchHBETEs3FbEnI15FFTYgcZB0sNTIxjfJ5bkiAC31tfEx2IiPtSf5Ah/AqzuXZ1bpK1QKHIDhSKR45NhGOTvndJf5qYp/U0yNqzl60/e44f/fdq8Ia3ZbOai62/j2jsexOk0WLO7lO825JGRX9l8Xe+4YCb0jaNvQojHNoaNCLKSFO5PdJBNm9FKu6ZQ5AYKRSLHv7LqerJKqskrr3XLViJNamuq+WnO//j6k3f5dfkS/m/685x+1gUAVJSVUl5agj0wmu835LM8s5imXr7EMH8m9I1leGqE2zehbWLzNRMf6kdcqD9Bbt7bTcQbFIrcQKFIpP1onNJfze6SGhrctJVIk6wdW4mOS8DPv7GL7KM3XuDlp/7CoOGjmXTeZfQedRoLtpfxU0Yh9r3LC4QH+DKudywndo/yaDdXiL/v3oCkRSGl/VAocgOFIpH2x+E0yCmrYVdxNdV29447ajLj73/iiw/ewLl39cngkDBOOfMcThg/hVy/ZH7YXNC8Ca2/r4WTukdxWu9YIjw4g8xshqggG3Ghfupek+OeQpEbKBSJtG+FlXYyi6opqXLflP4m+TnZfDvrA7759D3yc35bGLJT1x7M+HguyzPL+HZDLjl7F4O0mEwM7xzBhL6xJIV7ZlB2E6uPmbhQP+JD/bRqthyXFIrcQKFIpGNo2krE3eOOoHFhyFVLfmLulx+z8PuvGX7Safz56ZeBxsUgX331VXKCupNVH9R8Tb+EECb0jaNXXLDHW3CC/HxICPUnNtSGzUcLQ8rxQaHIDRSKRDqW2noHu4qryS6tweHmcUfQODi7sryMqNh4ADK3beH6s04CoFv/YQT3PZnCmMGY/UMBSIkIYELfWIZ1ivD42kMmE0QG2UgI9SMqSNuKSNumUOQGCkUiHVPTekdZJdXY6923z9rvbf51Fa88/Rhrly+m6c+wxeJDfN8R1Hc+AWuXEZitfkQGWjm9TyxjukV5ZVsPH4tpb/eaP6H+6l6TtkehyA0UikQ6NqfTILe8lkw37rN2IAW5e5j/zWf88OWnbN34a/PrZ979NJnBfamobaylaaXs03rHei2cBNgsJIT6Exfqp33XpM1QKHIDhSIRaeLJQdn7yty2hXlfzWLpT3N55p0vwOzL4u1FvPnS85TkZhHYZyyBKX0Z1TWa8X1iSQjz92h9TUwmCA+0khDqT3SwTVuLiFcpFLmBQpGI/J4nB2UfjGEYXDVpFDlZOwGwhEQT2HssgX3GMmzwQCb2i6d7TJDXptVbLCZig/1ICPPz2F5vIvtSKHIDhSIRORhPD8rel2EYrFryEz/871N+mvM/qisrmo/5RqUQNGgSAydczIS+sQxJDvfqoOgAq4W4UD8SwvzVvSYeo1B0AKWlpYwbN46GhgYaGhq44447uOGGG1p9vUKRiByOtwZlN7HX1rD0x7n88L9ZLFkwh4b6OoIHjidi4u0ARAVZGZNo5fSh3b0+pT480JeEMH9ig/00e03cSqHoABwOB3a7nYCAAKqqqujXrx/Lly8nMjKyVdcrFIlIa3lrUPa+KsvL+HnuV8R06sEucxzzNuVTvHMDue/cR1CXIYw8/SyuufQCYiLDvVJfEx+LifhQfxLCtDikuIdC0WEUFxczZMgQli9fTlRUVKuuUSgSkaPhrUHZv2evd/D00/9k3tv/bn7N5GMjeeBoxo4/kymT/0BIaKgXK2zcey0x3J/YYJvHNsSV9s+Vn99t5rfyxx9/ZPLkySQkJGAymZg9e/Z+56Snp5Oamoqfnx8jRoxg6dKlLY6XlpYycOBAkpKSuO+++1odiEREjlZUkI2hncIZ3iWCuFA/vLWNmM3Xwp8euJ/XvlzIaZfdin9UEkaDnV0rfuDt6fdw/pi+PPvRD6zfU0aD0/NdfwDlNfVs3FPOT1sL2bCnvHkPOJG2os20FH399dcsXLiQoUOHcu655zJr1iymTJnSfHzmzJlceeWVvPjii4wYMYJnnnmGjz76iM2bNxMTE9PiXnl5eZx77rl8+umnxMbGtur5aikSEVfw5qDsfTmdTn5YuJT/fT6LLUt/oKGmgqSpb2EymQm0WvDf8i1JoX5MOWcKSSmpXqszyM+HxLDGtY981XokR6Hdd5+ZTKb9QtGIESNIS0tjxowZQON/8MnJydx222088MAD+93jlltu4dRTT+X8888/4DPsdjt2u735+/LycpKTkxWKRMQlvD0oe18Op8GKzZlsLjOxclcJ5TV1ZP/nahyVxQCEJXVn+CkTmHL2FLr16uuV6f0Ws4noYBtJ4f6a2i9HpF12nx1KXV0dK1asYNy4cc2vmc1mxo0bx+LFi4HG1qGKisapqmVlZfz444/07NnzoPecPn06oaGhzV/JycnufRMi0qH4WsykRgUyumsU/RJDCQvw3iBji9nE8N6pXDGyE/88fyB3ntKFoWddS2DqQDCZKd2dwXdvz+CW88dx7snDeO6Zf9Hg8GyQczgNcstqWb6zhEXbCtlVVE1dg3fDpHQ8Pt4uoDUKCwtxOBz7dYXFxsayadMmADIzM7nxxhsxDAPDMLjtttvo37//Qe/54IMPcvfddzd/39RSJCLiSmZz495hcaF+VNTWs7ukhtzyWq91rVnMJvqnRNH/oXtxOu9hVcYuvvjyS9YunEPF1hVUFmazYHUGGR+uYVByGAMTAmnYvY6hI8fga/VMC0613cGWvAq2FlQQHeRHYrg/EYFqPRL3Oy5CUWsMHz6c1atXt/p8m82GzWZzX0EiIr8T7OdL73hfuscEkVNWy+6SGq9N6YfGwDa0ZyeG9pyK865bWLcrjy//9zXZpihq6h0s3l7ED99/R/5H0/D1D6LfyJOZ9IezGXnSqfgHBLq9PqcT8sprySuvJcBqIT6scWq/t9dgkvbruAhFUVFRWCwW8vLyWryel5dHXFycl6oSETk6PhYzyREBJEcEUFJVx+6SGgoqa/HSpDCgMSANSI1jwNRrcBoG2/IrWZ5ZwvcZ1VgCw6mvKmHVvC9ZNe9LLL5Weg4dw4QzJzP29DMJDHb/OMzqOgfb8ivZXlBJVJCNxHB/IgOtXtveRNqn4yIUWa1Whg4dyty5c5sHXzudTubOncutt97q3eJERI5BeKCV8EAr9oYg9pTWkl1SQ229w6s1mU0muscG0z02mIvS7iRj6rV888NP/DLva0o2LKShNJcNS35gw5IfWFcbybgxafRLCMWCE4uPez9WDAMKKuwUVNjx87UQH+ZHorYVERdpM6GosrKSrVu3Nn+/Y8cOVq9eTUREBCkpKdx9991cddVVDBs2jOHDh/PMM89QVVXFNddcc0zPTU9PJz09HYfDu3+ERKRjs/lY6BwVSGpkAAWVdnaX1FBc6d0FIaExIPWMC6XnpX/AecmZbMuvZM7CZSye+zVluzezuT6cLfO3YfMxUzvnGZylOZw68UzGnn4mKV26u7W22noHOwqq2FlYRUSglcRwf6KDbGo9kqPWZqbkz58/n1NOOWW/16+66ireeOMNAGbMmMFTTz1Fbm4ugwYN4rnnnmPEiBEueb7WKRKRtqa6roHskhqyS2to8OKaRwfiNAx2FFaxPLOEFZklFFXUsPv5y3DWVjafE53UmdGnnM6osePoN3QEVqv7x3FafcwkhPkRG6JtRTqKdr9OkTcoFIlIW9W019rukhrK2+Aq0IZhsKOoip9WZ/Dj3G/I//VnajPXgPO3QeRd+w/jP+9/gdmDrTj+VgtRQTaigqyEB1i1MW07pVDkBgpFInI8KK+tZ3dxDXnltTicbe/Pt2EY7CyqZtHGXcyfN5eCDb9Qu305QQMn0m3SNYzsHMnAWCtP3Hoxw04YS9qJp9J38HC3T/e3WExEBlqJDrYRGWjD6nNcLNMnraBQ5AYKRSJyPKl3OMkprWV3STXVdW1zTKTTMNicW8HirQUs35FPHY3dWdWbF1Ew+/Hm8wICgxgyaizDTzqNtDGnEBUb79a6TCYI9fclKshGdLCNQFubGV4rR0GhyIX2HWi9ZcsWhSIROe4UV9Wxu6Saggo7bfUvur3BwZqsMhZvL2Lt9myqt62gZscKaravwFld1uLcB55M57Q/nOex2gKsFqKCbUQF2QgP8NVA7eOMQpEbqKVIRI53tfUOsktr2FNa4/X91g6lrKaeZTuLWbK9iB2FldTlbqVm+wrsO5ZTu2cL096awwmD+2I2mfj+i49ZPO9bhp94GmljTiUiOubwDzgGPhbT3nFINiKDrNqk9jigUOQGCkUi0l4YhkFBhZ2skhpKqrw/rf9Q9pTWsGRHEUu2F1NcVYejpgKLfzARgVZGdo5gyUsPsOrn75vP795nQGNAOvFUeg0YgsXivvWJTCYIC7ASHWQjKthKgFXdbG2RQpEbKBSJSHtUU+cgp6yG3LLaNjv2CBrHH2XkVbJkexHLM0uo2buApT1nC9Y9q6nbsYLcbetbXBMWGcW736/wyFR/gACbhZi93Wyh/upmaysUitxAoUhE2ruy6npyy2vJLa+lvg3vQF/vcLImq5TF24tYl12OY+/HlFFdQnjxBpy7VrN99SK69OjNv96a3Xzdvx65h+i4eIafeBrd+w7EbHZf15evj5nIQCsxwTYiAq34qJvNaxSK3EChSEQ6CsMwKKysI7es1ut7rh1ORW09y3aWsGR7EdsLq5pft1kMeocZnDakFz1jg6koK+aCE/vR9JEWEhbB4JFjGDziRIaMOon45E5uq9Fs/q2bLTrYpi1HPEyhyA0UikSkI2pwOMmvsJNTVktpdV2bnb0GkFtey5LtRSzZXkThPlughAf4MjjOD9OOJWSs+JEVixZQXVXZ4tqzLrmG2/483SN1hgb4khjmT1yInxaM9ACFIhfSlHwRkUa19Q5yy2rJKaulyt5w+Au8xDAMtuZXsnjv+KN9x0olh/uTlhJCRPUuMlYtZtWSn9mwZjl3PPIPJp13KQCZWzcz/f6pDBl1IoNHnkT/oSPw8w9weZ0+FhMJYf4khvlrLSQ3UihyA7UUiYj8pqK2ntyyxvFHbXl6f73DydrdZSzZXsTa7LIWq3ynRAQwICmUHuE+pEYFERAYCMCsd17lP9P/3Hyer6+VPoOGMXjkiQwedRI9+w7E4uPaEBMe6EtiWAAxwTa1HrmYQpEbKBSJiOzPMAyKq+rIKauloNKOo41tTLuvytoGlmU2rn+0vaCKfSsN9vOhX0IoA5JCSfKrZ+OKRaxa8iMrF/9Ifk52i/s8/uJ7pJ14KgB19lp8rTaXzTTz9TGTEOpHYri/pvi7iEKRGygUiYgcmsPZuP5RTlkNxVVte/xRWU0967LLWJtdxoY95c1T/AEsJhNdYwIZkBhG/8QQjPJcVi35iVWLf2T9qmW88fVi/AMaW5Ve+edfmfu/Txk8cgxDRp7E4JEnumwbkvBAK0nh/kQHqfXoWCgUuYFCkYhI69kbHOSVNQakitq2O/4IoMHpZGt+JWt3N4ak3LLaFsejgqwMSAxr7GqLDcLq89vssdsvPZONa1a0OD+5S3eGjGyc1TZi7OnHvICk1cfcPPbI36qZa0dKocgNFIpERI5Olb2BnLJacstqqa1vuwtENimosLN2dylrs8vYnFtBwz7jkKw+ZnrHBTMgKYz+iaEEWhysX72cVYt/ZNWSn9iyfk3ztP/wyGhmLljb3LW2a3sGcUkpR72YpMnUsvVIi0O2jkKRGygUiYgcu5K944/yK2ppaMPjj5rY6x1szK1g7e5Sfs0uo6S6vsXxpHB/BiSGMiApjC5RgVSWl7Jm2SJWLf6RgKBgrr+7ccC2YRhcOLY/NVVV9B0ynIFpJzBg2Ch69huEr9V6xHXZfH9rPdK6R4emUORCmpIvIuJ6TqdBYZWd/HJ7mx+g3cQwDLJKapoD0u8HawfZfOibEMKApFD6JoQStM80+8K8HG65YDwlRQUt7mnz86f3wKGMm3w+E865+IhrMpkgMshGYpg/UUFWtR4dgEKRG6ilSETEPY7HgASNyxKs21POr7vLWLenrMV6SCYTdIsOon9i44y2xDB/AHZu3cSapQtZu2wxa5cvpqykGIALr53KDfc8DEBVZQWfvPkiA9NOoNeAIdj8/FtVj5+vhYQwPxLUetSCQpEbKBSJiLjf8RqQHE6DbQVNg7VL2VPacrB2RKC1OSD1igvG5mPBMAwyt21m7bLF9B44lO59BgCw9Ke5PHTTZUDjGkm9BgxhQNooBqSdQJ+BQw+7kKTJBFFBNpLC/YkM8sxmuG2ZQpEbKBSJiHjW8RqQAIoq7c2z2TblllO/T+02HzNDO4UzqkskPWOD95tuv3HtSma9/Qprli2iuCCvxTEfH18e/Md/OGnC5FbV4W+1kBjmT3yYHzafjtl6pFDkBgpFIiLe43QaFFXVkVfe9heJ/D17g4PNuRXNIam4quW+bCM6R3JC10gSwlp2kxmGQfauHaxduoi1yxu72wpy9/DKZ/NJ7dYLgG8+fY+vPn6XAcNGMSBtFH0HDycwKHi/GsxmiA7yIzbERliAFauP2b1vug1RKHIDhSIRkbZh34BUWGk/LmaxNTEMg60FlSzZXsyyncUtxiGlRARwQtdIhqdGEOLve8Brc3fvIi4ppXlA9eP33cy8r2Y1n2M2m+nWuz8DhzfObhtywtgDLgEQaPMhPNCX8AAr4e08JCkUuYFCkYhI23M8B6SmfdkWbyvi1+wyHHs/bs0m6JcQyqiukQxMCjtkYMnfs5s1yxaxZllja1JOVmbzMbPZzKeLNhEY3PiZtWt7BuGR0QSHhu13n31DUliAb7vqalMocgOFIhGRtu14DkgVtfUs21nCom2F7Cyqbn7d39fCsE7hjOoaSbeYIMyHmXKfn5Pd2NW2bDGV5aU88sxrzcfuuGwyG1YvI7VbT/oOGUG/IcPpN2Q4sQnJ+03lD7BZiAi0touQpFDkQlqnSETk+NMUkPIraimoOL4CUk5ZDYu3F7Fke3GL8UdRQVZGdo5kZNdI4kL8juieDoeDW84/ne1bNux3LDImjhNOncjtDz9x0OsDbBbCA6xEBB5/IUmhyA3UUiQicnxyOg2Kq/cO0j6OApLTMNiSV8GS7cUszyymtt7ZfKxLVCCjukSSlhpBkJ/PIe7SUmlxIetXLWPdyl9Yt3IpGRvW4mhoYPS4M3j02deBxrFLj993E5269qTfkOH07D+4eQPcJsdTSFIocgOFIhGR419TQCqqrKOspp5Kez1O5+Gv8zZ7g4M1WWUs2l7Ihj3lNG3HZjGbGJDYOP6of2IovpYjGzBdW1PN5l9X4Wvzo8/AoQDkZGVy5cQRzedYfHzo1qsf/YaMoO/eLrfwyOgW92kKSeEBVsID21ZIUihyA4UiEZH2x+k0qKhtoKymvvmrrW9aW1ZTzy87GrvXdhX/Nv4owGpheGoEo7pG0iUq8Ki3/CgvLWbeV7ObW5QKcve0OH7eVX/kpv/7CwB19lry9uwmKbVri+e1pZCkUOQGCkUiIh2DvcFBWU095XtDUnlNAw5n2/wozC5pGn9URGnNb5vVxgTbGNUlkpFdIokOPrZVrfP37GbdyqWNXW6rlnL1bfdzwqkTAVi55Cfuv+4CwiIi6Tt4+N6WpBF069WvxUa3AVYLEUGN3W0RAVZ8jrBF61goFLmBQpGISMdkGAaV9patSdX2ttWa5HQabMqtYPH2IlbsKqGu4bc+we4xQYzqEsmw1HACrK0ff9Qa33z6Hs899iD1dfYWr9v8/OnVfzDX3vVQc7dcE5MJQv19iQyyERFoJcTPx60b2SoUuYFCkYiINKl3OJtbkpq+2soA7tp6B6t2lbJoeyGbcipoqsrHbKJ/Yij9E0PplxhKRKD1kPdprbo6O1s3/No8eHvdyqVUlJUA8MLH39Otdz8Afv7+K5YvnE+fgUPpPXBoc5ebr4+ZiAArEUFWIgOtLt/MVqHIDRSKRETkUKrr9mlNqq6n0t6Atz9Bi6vq+GVHEYu3F+23SW1imD/9EkPonxhKt+ggl3VpOZ1OsnZsZf2qZUw452IslsaQ89RDd/Dd7JnN5wWHhtN74BD6DBxG70HD6D9kBL5WKwE2C1F7W5HCA6xYzMfWiqRQ5AYKRSIiciQcToOK2patSfZ670x1MwyDrJIa1mSV8mt2GTsKq9j3w93mY6Z3fGNA6pcQQmTQsY1DOpCVi39k+cJ5bFi9gi3r1+zX5TZr8WaCQkIByNiwFj//AFK6dCUswEZkoJXIICvBfvtvf3I4CkVuoFAkIiLHqrbe0bioZHktJdV1XlsOoLK2gQ055fyaXca6PWVU1Da0OB4f6ke/xFD6J4TSPTboiKf6H059XR3bNq1jw5oVbFyznIqyMp545YPm43dfNYVfly/ZrzVpwKChJMVGELl30HZrZrUpFLmBQpGIiLhSvcNJQYWdPC8HJKdhkFVc3RiQssvZVljZotvP6mOmV1ww/RMaxyId62y2wzEMgwdvvIS1yxfv15pkMpkYOHw0T73+MQBBfj5EBVmJCLQR5u+L+QBdbQpFLqRtPkRExN2aAlJ+hZ3iKrtXF5SssjewsbkVqZyyfab6A8SG2PZ2s4XSMy7Y5a1ITX7fmrRxzQry9uxm5NjTeew/bzef98dzTyMqNo6+g4YxfMRITjxhJJ3iowiyNc60UyhyA7UUiYiIJ9Q7nBRW2skr935AahqLtG5vN9vW/Er2XbLJajHTIy6ouRUp9gj3ZDtShfm51FZXkZTaFWjcAPeycb+f8m8itXsv+g9OY9KZkxl3ykl0T45VKHIlhSIREfG0BoeTgko7+eV2irwckKBxht3GnIrmkFRS3bIVKSbYRr+EUPolhtAzLtjtK1k31NezdeOvbFy7kg2rl7NxzXLy9uxuPj7l8uu5+rb7mTKiu0KRKykUiYiINzU4nBRWNm5s2xYCkmEYZJfWsC67nHV7ysjIr2yx8reP2UTPuGD6JTSujRQbYnPrIo1NigrymgPSsDGn0rPfIIUiV1MoEhGRtmLfgFRcVdcmtiGprXewMaecdXsaxyMVV9W1OB4dbGNoSjjDUsPpFBHgkYAEUFVZoVDkagpFIiLSFjU4nBRV7W1BqmwbAckwDHLKapun/GfkVdKwT11RQVaGdgpnWKcIUiPdG5AUitxAoUhERNo6h9PYO0i77QQkaGxF+jW7jOWZJfy6u4w6x299f5GBjQFpaKdwukQFujwgKRS5gUKRiIgcTxxOg6K9s9gKK+1tJiDZ6x38uqeMFZklrN1dhn2fzWsjAqwM6RTGsE4RdIkOxOyCgKRQ5AYKRSIicrxqCkj5FXZKquu8tt3I79kbHKzfU87ynSWs2V3aIiCFB/gyJKWxBalbTNBRBySFIjdQKBIRkfaioraekqp6CqvslFXXt4lWpHqHk3V7u9jW7C6ldp/gFurvy9C9Aal7TNABV64+GIUiN1AoEhGR9sjpNCiprqO4qo6iqjoqf7cPmjfUO5ys31POiswSVmeVUlPvaD4W4ufDkL2z2LrHBGM5TEBSKHIDhSIREekI7A2OxlakyrbR1VbvcLIxp5zlewNSdd1vASl4b0AamhJOz7gDBySFIjdQKBIRkY6o0t5AcWVdm+hqa3A42ZRbwfLMElbtKqFqn4AUZPNhcHIYw1IbA5KPuXFPNoUiF9KGsCIiIo3aUldbg9PJ5twKVmSWsHJXKZX232oJtFoYnBLOsE7hJAfB+Sf0VChyJbUUiYiItNRWutocToMteY0tSCt3lVCxT1izGbVs+cf5Lvn89jnWQkVERKR9svlYiAu1EBfqB3ivq81iNtE7PoTe8SFcNjyFLfmNLUgrMksoLXNdUFNL0V5qKRIREWm9pq62kuo6Ciu909XmdBqs3ZHDfWcNUUuRiIiIeIfZbCIyyEZkkI1uMVBT5yC3vJbcslqq7J4JSGazie6xwS67n0KRiIiIHDN/q4XOUYF0jgqkoraevPJacsvs1O6zBlFbp1AkIiIiLhXs50uwny/dYoIpqaojt7yWvPJaGhxte8SOQpGIiIi4TXiglfBAKz1jgymqqiO3rLZNbWC7L4UiERERcTuz2UR0sI3oYBsNDicFlXZyymopqaqjrUz5UigSERERj/KxmIkP9Sc+1B97g4P8cju55bWUVdd7ty6vPl1EREQ6NJuPheSIAJIjAppnsOWU1VBt9/wAbYUiERERaRP2ncFWXltPXlktueW1HltJW6FIRERE2pwQP19C/HzpFhNESXU9uWW15Fe4dwabQpGIiIi0WSaTiYhAKxGBVno5gymssjfPYHO6uAFJoUhERESOC2aziZhgP2KC/WhwOMmvsLM1u9Zl91coEhERkeOOj8VMQpg/QeZwl93T7LI7iYiIiBzHFIpEREREUCgSERERARSKSE9Pp0+fPqSlpXm7FBEREfEik2G0lR1HvKu8vJzQ0FDKysoICQnxdjkiIiLSCq78/O7wLUUiIiIioFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAigUiYiIiAAKRSIiIiKAQpGIiIgIoFAkIiIiAoCPtwtoK5p2OykvL/dyJSIiItJaTZ/brti1TKFor6KiIgCSk5O9XImIiIgcqaKiIkJDQ4/pHgpFe0VERACwa9euY/5HPZ6Ul5eTnJxMVlZWh9oIV+9b77sj0PvW++4IysrKSElJaf4cPxYKRXuZzY3Dq0JDQzvUL1OTkJAQve8ORO+7Y9H77lg66vtu+hw/pnu4oA4RERGR455CkYiIiAgKRc1sNhvTpk3DZrN5uxSP0vvW++4I9L71vjsCve9jf98mwxVz2ERERESOc2opEhEREUGhSERERARQKBIREREBFIpEREREAIWiZunp6aSmpuLn58eIESNYunSpt0tyq+nTp5OWlkZwcDAxMTFMmTKFzZs3e7ssj3riiScwmUzceeed3i7FI7Kzs7n88suJjIzE39+f/v37s3z5cm+X5VYOh4OHH36Yzp074+/vT9euXXnsscdcskdSW/Ljjz8yefJkEhISMJlMzJ49u8VxwzB45JFHiI+Px9/fn3HjxpGRkeGdYl3oUO+7vr6e+++/n/79+xMYGEhCQgJXXnkle/bs8V7BLnK4n/e+brrpJkwmE88884zH6nOX1rzvjRs3ctZZZxEaGkpgYCBpaWns2rWr1c9QKAJmzpzJ3XffzbRp01i5ciUDBw5kwoQJ5Ofne7s0t1mwYAFTp05lyZIlzJkzh/r6esaPH09VVZW3S/OIZcuW8dJLLzFgwABvl+IRJSUljB49Gl9fX77++ms2bNjA008/TXh4uLdLc6snn3ySF154gRkzZrBx40aefPJJ/vGPf/D88897uzSXqqqqYuDAgaSnpx/w+D/+8Q+ee+45XnzxRX755RcCAwOZMGECtbW1Hq7UtQ71vqurq1m5ciUPP/wwK1eu5NNPP2Xz5s2cddZZXqjUtQ73824ya9YslixZQkJCgocqc6/Dve9t27YxZswYevXqxfz581m7di0PP/wwfn5+rX+IIcbw4cONqVOnNn/vcDiMhIQEY/r06V6syrPy8/MNwFiwYIG3S3G7iooKo3v37sacOXOMsWPHGnfccYe3S3K7+++/3xgzZoy3y/C4M88807j22mtbvHbuuecal112mZcqcj/AmDVrVvP3TqfTiIuLM5566qnm10pLSw2bzWa8//77XqjQPX7/vg9k6dKlBmBkZmZ6pigPONj73r17t5GYmGisW7fO6NSpk/Hvf//b47W504He90UXXWRcfvnlx3TfDt9SVFdXx4oVKxg3blzza2azmXHjxrF48WIvVuZZZWVlAC7ZUK+tmzp1KmeeeWaLn3l79/nnnzNs2DAuuOACYmJiGDx4MK+88oq3y3K7E044gblz57JlyxYA1qxZw88//8ykSZO8XJnn7Nixg9zc3Ba/76GhoYwYMaJD/Y2Dxr9zJpOJsLAwb5fiVk6nkyuuuIL77ruPvn37erscj3A6nfzvf/+jR48eTJgwgZiYGEaMGHHIrsUD6fChqLCwEIfDQWxsbIvXY2Njyc3N9VJVnuV0OrnzzjsZPXo0/fr183Y5bvXBBx+wcuVKpk+f7u1SPGr79u288MILdO/enW+//Zabb76Z22+/nTfffNPbpbnVAw88wMUXX0yvXr3w9fVl8ODB3HnnnVx22WXeLs1jmv6OdeS/cQC1tbXcf//9XHLJJe1+s9Qnn3wSHx8fbr/9dm+X4jH5+flUVlbyxBNPMHHiRL777jvOOecczj33XBYsWNDq+/i4sUY5TkydOpV169bx888/e7sUt8rKyuKOO+5gzpw5R9bH3A44nU6GDRvG448/DsDgwYNZt24dL774IldddZWXq3OfDz/8kHfffZf33nuPvn37snr1au68804SEhLa9fuWlurr67nwwgsxDIMXXnjB2+W41YoVK3j22WdZuXIlJpPJ2+V4jNPpBODss8/mrrvuAmDQoEEsWrSIF198kbFjx7bqPh2+pSgqKgqLxUJeXl6L1/Py8oiLi/NSVZ5z66238uWXXzJv3jySkpK8XY5brVixgvz8fIYMGYKPjw8+Pj4sWLCA5557Dh8fHxwOh7dLdJv4+Hj69OnT4rXevXsf0ayM49F9993X3FrUv39/rrjiCu66664O1VLY9Heso/6NawpEmZmZzJkzp923Ev3000/k5+eTkpLS/HcuMzOTe+65h9TUVG+X5zZRUVH4+Pgc89+5Dh+KrFYrQ4cOZe7cuc2vOZ1O5s6dy6hRo7xYmXsZhsGtt97KrFmz+OGHH+jcubO3S3K70047jV9//ZXVq1c3fw0bNozLLruM1atXY7FYvF2i24wePXq/JRe2bNlCp06dvFSRZ1RXV2M2t/wzZ7FYmv9fZUfQuXNn4uLiWvyNKy8v55dffmnXf+Pgt0CUkZHB999/T2RkpLdLcrsrrriCtWvXtvg7l5CQwH333ce3337r7fLcxmq1kpaWdsx/59R9Btx9991cddVVDBs2jOHDh/PMM89QVVXFNddc4+3S3Gbq1Km89957fPbZZwQHBzePLQgNDcXf39/L1blHcHDwfmOmAgMDiYyMbPdjqe666y5OOOEEHn/8cS688EKWLl3Kyy+/zMsvv+zt0txq8uTJ/P3vfyclJYW+ffuyatUq/vWvf3Httdd6uzSXqqysZOvWrc3f79ixg9WrVxMREUFKSgp33nknf/vb3+jevTudO3fm4YcfJiEhgSlTpnivaBc41PuOj4/n/PPPZ+XKlXz55Zc4HI7mv3MRERFYrVZvlX3MDvfz/n348/X1JS4ujp49e3q6VJc63Pu+7777uOiiizjppJM45ZRT+Oabb/jiiy+YP39+6x9yTHPX2pHnn3/eSElJMaxWqzF8+HBjyZIl3i7JrYADfv33v//1dmke1VGm5BuGYXzxxRdGv379DJvNZvTq1ct4+eWXvV2S25WXlxt33HGHkZKSYvj5+RldunQxHnroIcNut3u7NJeaN2/eAf97vuqqqwzDaJyW//DDDxuxsbGGzWYzTjvtNGPz5s3eLdoFDvW+d+zYcdC/c/PmzfN26cfkcD/v32svU/Jb875fe+01o1u3boafn58xcOBAY/bs2Uf0DJNhtLOlXUVERESOQocfUyQiIiICCkUiIiIigEKRiIiICKBQJCIiIgIoFImIiIgACkUiIiIigEKRiIiICKBQJCIiIgIoFIlIB7B7927OP/98Hn74YW+XIiJtmEKRiLR7d911F927d+ejjz7ydiki0oYpFIlIu1ZWVsb8+fMZM2YMCQkJ3i5HRNowhSIRadd++OEHTjzxRBYsWMDo0aO9XY6ItGEKRSLSrv3000+MHDmSL7/8knPOOcfb5YhIG6ZQJCLt2vLlyykpKcHPz48hQ4Z4uxwRacN8vF2AiIg77dy5k4KCAv7yl794uxQRaeNMhmEY3i5CRMRdbDYbvXr1YtWqVZjNahwXkYPTXwgRadd8fX15+umnFYhE5LD0V0JE2q0333yTqqoqbDYbS5Ys4csvv/R2SSLShmlMkYi0S7W1tXz66ae8+eab3HjjjaSmpvLWW295uywRacM0pkhEREQEdZ+JiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICKBSJiIiIAApFIiIiIoBCkYiIiAigUCQiIiICwP8DUJpSq0pm5M8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "plt.title(r\"$C_{\\ell}$ but $C_{1}$ divided by 9\")\n",
    "ylow, ymed, yhigh = np.percentile(Cell[:, k, :], [16, 50, 84], axis=0)\n",
    "plt.plot(ell, ymed, label=\"Quijote\")\n",
    "plt.fill_between(ell, ylow, yhigh, alpha=0.3)\n",
    "plt.plot(ell, ytheory, label=\"Linear theory\", color=\"black\", linestyle=\"--\")\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(r\"$\\ell$\")\n",
    "plt.ylabel(r\"$C_{\\rm ell}$\")\n",
    "plt.xlim(ell.min(), ell.max())\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "plt.savefig(\"../../plots/C_ell.png\", dpi=450)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "fig, axs = plt.subplots(len(radii), 1, figsize=(5, 4 * len(radii)), sharex=True)\n",
    "\n",
    "for k in range(len(radii)):\n",
    "\n",
    "    ylow, ymed, yhigh = np.percentile(C_ell[:, k], [16, 50, 84], axis=0)\n",
    "    axs[k].plot(ell, ymed, label=rf\"$R = {radii[k]} ~ [\\mathrm{{Mpc}} / h]$\")\n",
    "    axs[k].fill_between(ell, ylow, yhigh, alpha=0.1)\n",
    "\n",
    "    Cell_theory = [C_ell_theory_equal_radius(ell_, radii[k]) for ell_ in ell]\n",
    "    axs[k].plot(ell, Cell_theory, color=\"black\", label=\"Theory\")\n",
    "\n",
    "    axs[k].set_yscale(\"log\")\n",
    "    axs[k].legend()\n",
    "\n",
    "    axs[k].set_ylabel(r\"$C_{\\ell}$\")\n",
    "    axs[k].set_xlabel(r\"$\\ell$\")\n",
    "\n",
    "axs[-1].set_xlim(0, ell.max())\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../../plots/C_ell_fiducial.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(skymap[4], unit=r\"$V_{\\rm rad} ~ [\\mathrm{km} / \\mathrm{s}]$\",  cmap=\"coolwarm\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 4\n",
    "\n",
    "Cell_theory = [C_ell_theory_equal_radius(ell_, radii[k]) for ell_ in ell]\n",
    "\n",
    "plt.figure()\n",
    "plt.title(f\"R = {radii[k]}\")\n",
    "\n",
    "\n",
    "for n in range(50):\n",
    "    plt.plot(ell, C_ell[n, k], alpha=0.1, color=\"black\")\n",
    "\n",
    "\n",
    "plt.plot(ell, C_ell[0, k], color=\"red\", label=\"Simulated\")\n",
    "plt.plot(ell, Cell_theory, color=\"blue\", label=\"Theory\")\n",
    "\n",
    "plt.yscale(\"log\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_csiborg",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}