brandonshensley/BB_Validation.ipynb Secret

## BB_Validation.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8e66400a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysm3\n",
    "import pysm3.units as u\n",
    "import healpy as hp\n",
    "import numpy as np\n",
    "import pymaster as nmt\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=True)\n",
    "res_dpi = 300\n",
    "ext = 'pdf'\n",
    "import cmocean\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bedb034a",
   "metadata": {},
   "source": [
    "This is the Keck Array/BICEP2 mask found here: http://bicepkeck.org/bk14_2015_release.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e41012ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAFxCAYAAABA772BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArU0lEQVR4nO3dW3Ic2Z3f8d//nMysQgFgs9mtmyfGE9bEvPlJnglvQB3hBUihFdjawUx4BQ4tQTMrkFsbcLR2MJYeHH5w2FZ7JuzRSK0mGyRQt7ycvx9OFlAAARIkkajb9xMBolCXrEIRLHx5zqlMc3cBAADgYYVNPwAAAIB9RGQBAAAMgMgCAAAYAJEFAAAwACILAABgAEQWAADAAIgsAHcys9+a2ee3nP8zM3vj/l/M7KmZ/frm6fd4DN83s9/e9fjeZ5v32TYAfKhi0w8AwNb7/i3n/VDS2WPcubt/KenPd23bAMBIFoC3+YWZ/Wj1hZn9QNKvNvh4AGAnEFkA3uaXkn6y9vVPJP1i/Qpm9h/M7Iv+44dv2lh/ne/3p39uZn/dn/6BmX1uZn9tZr/ur/f05lRjf5tfm9nPJT1bO//a7W6538/7QFx9/dtbtn3zvt/4WO/x3AE4YEQWgDfqp9TWpwx/6O6/WX3Rh8tn7v6Zu38m6We3Rc6aLyT9aO3rv+o//0TS3/fb+jeSfirpZ+s37APu++7+b9z9p7c8hltv1/tFfx+r6/9m/cI7tvGmx3otNAHgJiILwH38wsx+dMdU4U91PTh+pbxm6y6/lPRZP0L0W0lP+/N/KOlPJD3rR4l+Jukvb9z2M0nrI0gv+s8/ecvt5O6/1FUw3RZJt23jTY+VKVMAb8TCdwD38Uvl8PhS9xvBeXbXBe7+ZR8tq9GkT9am8RaS/sbdLwPmLaNi667d7g6r+/6hu//NLdt+bRt3PVZ3P7vn4wJwoBjJAvBWa1OGP1ifKux9rjyatfIjSf/5LZv8jaSf9EHzhaT/qDwy9Iv1ba3WQ635QtKP+8ue6irm3na7lZ/31/vylsvu2sZdjxUA3ojIAnBfP1eOjGtW8dEvGP+18mjQ2Vu29Qv167z62/9I0hd9wK0W0N9cD7W67pf9/axG1vS22635laS/1i2jcW/Yxq2P9S3fHwDI3N+4P0EAAAC8B0ayAAAABkBkAQAADIDIAgAAGACRBQAAMAAiCwAAYAD33Rkpb0EEAAB4nd11AXt8B3DNZ+HHm34IO++LxLGjAdx/P1mMZAE7jnjaPsQYsBfuHMkisoA9QUTtHyIM2AlEFrBriCa8K6IM2AgiC9hWxBSGRnwBgyKygG1AUGFbEF7AgyGygCERT9hXxBjwVkQW8JCIKhwqogt4DZEFfAiiCrgd0QUQWcAbEVHAMIgwHAAiC1hHVAGbQXRhDxFZOEzEFLAbiC/sMCILh4GoAvYD0YUdQmRhfxFWwH4juLDliCzsPmIKwDriC1uCyMJuIqwA3AfBhQ0isrA7CCsAH4LgwiMjsrCdCCoAj4HwwoCILGwPwgrAJhFceGBEFjaPuAKwTYgtPBAiC4+LoAKwiwgvvAciC8MjrADsE4IL90Rk4eERVQAOCdGFOxBZeDjEFYBDRmzhBiILH4awAoDXEVwQkYV3RVQBwLsjug4SkYX7Ia4A4MMRWweFyMKbEVcA8PCIrYNAZOE6ogoAHh/RtZeILGTEFQBsHrG1V4isQ0VUAcD2I7p22p2RFR7zUeBxEVgAsBt4vd5PjGTtGf6hAsDuY2RrpzBduM8IKwDYXwTX1iOy9hFxBQCHg9jaWqzJ2jcEFgAcFl73dw8jWTuCf1wAgJsY3doKjGTtMgILAHAbfj9sN0aythj/eAAA98Wo1saw8H1XEFYAgA9FcD0qImvbEVcAgIdGbD0KImsbEVYAgMdCcA2Ghe/bhsACADwmfu88PkayHhk/5ACATWNU60ExXbhpxBUAYNsQWw+C6cJNIrAAANuI30/DYiRrQPzwAgB2BaNa743pwsdCWAEAdh3B9U6YLnwMBBYAYB/w++xhEFkPhB9IAMA+4ffah2O68APwAwgAOBRMId6J6cKHRmABAA4Jv/feHSNZ74gfMgDAoWNU6xreXfihiCsAAK4jtiQxXfhhCCwAAF7H78c3YyTrDfjhAQDgfg54VIuRrHdFYAEAcH/83nwdkXULflAAAHh3/P68junCNfxwAADwMA5o+pB3F74JcQUAwDAOILZYkwUAAPCYDnYki9ErAAAe156OajGStY7AAgDg8R3a79+DGsk6tL9cAAC21R6NajGSRWABALA9DuH38kFE1iH8RQIAsGv2/ffzXk8X7vtfHgAA+2KHpw8Pb7qQwAIAYHfs4+/tvY0sAACATdqr6cJ9rGAAAA7RDk0fHt50IQAAwCbtTWQxigUAwP7Yh9/rOz9duA9/CQAA4G5bPnXIdCEAAMBj2tmRLEawAAA4LFs6orVfI1kEFgAAh2fXfv/vXGTt2hMMAAAezi51wE5F1i49sQAAYBi70gM7E1m78oQCAIDh7UIX7ERk7cITCQAAHte298HWR9a2P4EAAGBztrkTtjqytvmJAwAA22Fbe2FrI2tbnzAAALB9trEbtjKytvGJAgAA223b+mHrImvbniAAALA7tqkjtiqytumJAQAAu2lbemJrImtbnhAAALD7tqErtiayAAAA9om5+32ud68rvY9tKE0AALC/vkifD7l5u+sCRrIAAAAGsNHIYhQLAAAMbVO9sbHIIrAAAMBj2UR3bCSyCCwAAPDYHrs/Hj2yCCwAALApj9khLHwHAAAYwKNGFqNYAABg0x6rRx4tsggsAACwLR6jS5guBAAAGMCjRBajWAAAYNsM3SeDRxaBBQAAttWQncJ0IQAAwAAGjSxGsQAAwLYbqlcGiywCCwAA7IohumWQyCKwAADArnnofmFNFgAAwACILAAAgAE8eGQxVQgAAHbVQ3bMg0YWgQUAAHbdQ/UM04UAAAADILIAAAAG8GCRxVQhAADYFw/RNYxkAQAADOBBIotRLAAAsG8+tG8+OLIILAAAsK8+pHOYLgQAABgAkQUAADCAD4ospgoBAMC+e9/eYSQLAABgAO8dWYxiAQCAQ/E+3cNIFgAAwACILAAAgAG8V2QxVQgAAA7Nu/YPI1kAAAADILIAAAAGQGQBAAAM4J0ji/VYAADgUL1LBzGSBQAAMAAiCwAAYADvFFlMFQIAgEN33x5iJAsAAGAARBYAAMAA7h1ZTBUCAABk9+kiRrIAAAAGQGQBAAAMgMgCAAAYAJEFAAAwACILAABgAEQWAADAAIgsAACAARBZAAAAAyCyAAAABkBkAQAADIDIAgAAGACRBQAAMAAiCwAAYABEFgAAwACILAAAgAEQWQAAAAMgsgAAAAZAZAEAAAyAyAIAABgAkQUAADAAIgsAAGAARBYAAMAAiCwAAIABEFkAAAADILIAAAAGQGQBAAAMgMgCAAAYAJEFAAAwgHtH1hfp8yEfBwAAwM64TxcxkgUAADAAIgsAAGAA7xRZTBkCAIBDd98eYiQLAABgAEQWAADAAN45spgyBAAAh+pdOoiRLAAAgAEQWQAAAAMgsgAAAAbwXpHFuiwAAHBo3rV/GMkCAAAYAJEFAAAwgPeOLKYMAQDAoXif7mEkCwAAYAAfFFmMZgEAgH33vr3DSBYAAMAAiCwAAIABfHBkMWUIAAD21Yd0zoOMZBFaAABg33xo3zBdCAAAMIAHiyxGswAAwL54iK5hJAsAAGAARBYAAMAAHjSymDIEAAC77qF65sFHsggtAACwqx6yY5guBAAAGACRBQAAMIBBIospQwAAsGseul8GG8kitAAAwK4YolsGnS4ktAAAwLYbqldYkwUAADCAwSOL0SwAALCthuyURxnJIrQAAMC2GbpPmC4EAAAYwKNFFqNZAABgWzxGlzzqSBahBQAANu2xeoTpQgAAgAE8emQxmgUAADblMTtkIyNZhBYAAHhsj90fG5suJLQAAMBj2UR3bHRNFqEFAACGtqneYOE7AADAAMzd73O9e13pQ3wWfjz0XQAAgAPySCNYdtcFjGQBAAAMYGsii/VZAADgoWxDV2xNZEnb8YQAAIDdti09sVWRJW3PEwMAAHbPNnXE1kWWtF1PEAAA2A3b1g9bGVnS9j1RAABge21jN2xtZEnb+YQBAIDtsq29sNWRJW3vEwcAADZvmzth6yNL2u4nEAAAbMa298FORJa0/U8kAAB4PLvQBTsTWdJuPKEAAGBYu9IDOxVZ0u48sQAA4OHtUgfsXGRJu/UEAwCAh7Frv//N3e9zvXtdaRM+Cz/e9EMAAAAD2vK4srsu2MmRLAAAgG238yNZK4xoAQCwX7Z8BGtl/0eyduQvAgAA3MM+/F4vNv0AAGAnmEl2y/9LPa2d3vpBfwCPaG+mC29i+hDAe7N+9N+CLPRxFUy2Oj+sxVbKkeXuUtfl011HcAHvYUdHr+6cLmQkCwBuWsWVJMWY4yqEHF8hXMWW+riSZCnl/40ml8wlJUILOHB7O5IlMZoF4B31U4IWYx65ivEqrGK4HNHSWmTJPYeVJ3nTSu55JCutPneb+36AHbKjo1jSG0ay9jqyVogtAG8VoizGHFNlKSuKHFRFkUeuVsEVw/XpwtVraJekpskjW3UjdZ28beV1Le8jjJEt4HU7HFcrhx1ZEqEF4A3MZFWVR65ilFVlH1Sxj60glYU89lOGRbweV+6yLuX1WZefO3mX5NOpvGkZ1QJusQeBJRFZGaEF4DWWpwVtNJKVRQ6rsswjV6u4Cv1IVgzyEKRo8qIfzUp5PZY6l/XThVY3l7Hli4W8bqSmUVouGc0CensSWBKRdR2xBeBy/VVZyIpCNpnIiigVhXxc5cCKQV4V8mjyMsqDyYMpVVEyya1f4+6SJVdokpRcoc2BZW2SzWtZ3cjnC6XnL+Rtu+nvHNioPYqrFSLrJkILOHA3RrDs+FgqorwspFElL4K8CEpVXpuVqpAjK5pSkWPrclPukkuhdVmXP0LdydqksGhki0Y2Wyh9/Vypbpg2xMHaw8CSiKw3I7iAA7O2yN2OJ7KylJ9MpKqUV4W6o7KPLFM3ivIodVWQB/WRpXza7Orl1aXQuSytYit/Lqad4rJTmDWK//y10nSmdHHBtCEOxp6G1Tr2kwUAKxYsv3Mw5ulBlaVUlUqjUj6KSuOoVAalaOrGOa5SafnzKrJiH1pB+SXWJXOTJclaU+gjy6OUFkFFYQoXxzJ32XzOtCFwABjJWsOIFnAAzGRFebUW6/QkB9ZHx0rjQmkU1RzndVipNHWVKcWryMrBpX59Vo4tmfIC+H59VmhdoZWslaqpKy6TinlS9Xyu+M1U6Xe/Z9oQe+8ARrBWmC58F8QWsKdW67COjmRVmacJT4/lR5Xaj4/UjaK6UVAzCUqFlApTqqQUJS/WImttJCuVkoKUCs+jWpJibbIuR1Y5leLSVcyloz+2ql7WKr/8vdL5hXy5ZEQLe+eA4mqF6cJ38UX6nNAC9pHlvbZbUeRpwqLIu2goo7oyqBsFdZWpG1kfWf0IVtTlOiyFq9Ne9HG1mjqM+f+jq9Et6/KIVxyZUikViyjzSsWTk36xvBNZ2CsHGFhvRGTdgdAC9oyZLFjeo3vsQ6sslEal0qhQqq4Cqx3nUSyPq6nBHFTe/3/1ciQrurqR5IXLo8tLv1qf1Vhem9VI3ofacmFKMar85ESF+v/+LpZMG2IvEFivY7rwHogtYA+EqFCVec/uHz3J+8KqSjXPJmonUcuPC7XjHFjtsV0brZJW04N+fbowutI4SYVLZVIx6hRiUoxJ7qauDWouKlkdFOam0fOg0ErHv0savUwaf71Q8Q9/UJrO5IulvG141yF2DnHFdOEHYVQL2H35HYV5D+6KQSqiUlWoGwWlMqgdmbpKSlX+7KuQMknmSv3Qk5vkpV9ND5YuqzrFMuloslQRkiTJ3WQj1zS42jqqKws1tam8MNVP8o5QpbGOz58qlKX85SulWcqH3yG0sCMIrDcjsu5p9YNEbAE7yq4O6uxFzHtwH0WlUVCq8pqpbmTqRlIa5bVWsj60JCn0wRX64IouL64Cqxo1GpetxkWrUZHXWTVdVAyueVlqGV1tHeQxKC5NbvmjOp+oiKYQTNa2+ZiHTcsUIrYacXU/TBe+J2IL2C1WFLKjI4XjifzZR0rHI7XHpebfrtRMTMunpvZI6sZSe5KuRqpWbdZPCHjwPD0YXVYkVeNGVdXqeFTrO5MLnZRLPS3nCpaUPOisOdJZfaRXy7FeTCeazyrpq5HK86BQS5M/uKpzV/Wq0+QfzqSmlU3nSt+c5QNMDzWFaHYZnlYWMutH+qR8zMXVfSaXd50smDz153litO1AEVe3YrrwoTGFCOwYCzKzHC5FUCpC3kVD/87ByxGskSuNPMdUUN75VVD+r2Y/dWhlkkWXBVdVtZpUjSZlo0lR60m50NNyptI6ldZpFBodF0tNiomCuV6EpPN5oSZIxdRUf7R692JUsTiVNUnlN5UsJdlyKV8EedvmwPmQ0S2zq+chmKyq8s5YU8oHxraQ322ZOim5tHrXY0pSa/0mVpFlTGseIALr3TGS9UAILmDLhagwHik8OVX67idqT0dqTgtdfK9Qc2JaPpO6I1c3SdJpI4uuEFwWkswkT6tIcRVFXtxexk4fHS00iq2eVAv9q+PnelZM9WlxrnFoVFqbR7O6ic7TWP979m3Nu1L/48V3dD4bazkvpW8qFTNTMTONnud9apVz19HXrYppq/hyIVssZctGalv5YpmnFNs2R1LT6ObruFnem71VldR10miU31VZlvn4jDHIxyNZ1+X4cs+f207WtDms6kaeUr79ZXD1u5xI6Sr8GNXaa4TVvTCSNTRGtoDd4dHyzkVjf+zBIKXK1Y3yOquyunqXYBH7hez9bU1SDElV0akISUdFo5NyqdNiqdO40Cg0ehpnGodapTo1ipqEpc7TkS5GYz2vT/Tdk3ONi1Yv4kQzl5pRXoBvbZC5qZmbpEL6pNDorFSxmKi4aGRNlz9mi/yq7i51SV7XsrLsQ6mVQpBPxnn6LyX58Vgp9r8HQl7oL0nWJq0Gp9QmWZF3b6G2kyXP95GS3MJlTJl0OYpl6uQpSM76sX1EYH04IusBEVrA9vOmkczkIUfHaqejHiRFl6qksmpVxBxSVdEqmqvuoiSpCEnBXGXsLkewjmOt0zIH1mlY6DTMNQ6NKuX4WHipsTWaViNNQq3Gg6rQqoyd/hiSZsVYXRE1j4WKWd55qUKQXOoqU6yjYl2qmCeFOil0xwp1Hz3J8/SeJI9BaRwV6qRUBnmwHG8pj1T5+v+3bXUga1eoU46qJij0o2CSZF3K04lSHj3zfF+XX8eYQ2u1YUa09gaB9TCIrAe2/oNJcAFbJOV37dnIZU0n+ep4hP3e3UcuHyUV41ZPj+c6KWuNilYnxVJF6NSmqCRTFfLUWWlJR7FfbxVqTeJS3yle6klc6ElYaBIaja1TKVfSXLUHPQkLvSrH+l75jf7QfqSvm1N9dXKqPy5O9Go51sv5WPN5pWUdtVhE2SJq/h2TzBQXQWEZ5UGKtRTaHIeh1tVolHQtpFbnW9LlUFxoXNYpH8C6ycdYjHU+vmJoC8VlIWuSwjzmg1k3nbSs+whr82Zak7pO1oeWVtOG2GmE1cMjsgbEbh+ALeNJvlxefhma1eL2/E5CK5OKstO4aHVcLjUpGj0t5yrC9emwqHS5P6xnxVSj0Ki0TuPQaGyNgiVFuUq5yj56SkuSliqtU1JQZZ1Ow0KldTqOtc6qI/2xONHFuNK8LjWdj9Q2Uc2olCS1nSlOQ37HY+qPjSjleOoDKtamVHo+SHWdF9SbS2GZv8fQSKHJD6iYu9xyYJpLHoJinUeqQrC8jaaTmymklAO1iLIm73oij7R5XrOFnUZcDYfIegRMIwLbwbtOlvI0W1x2ak6jrPP8brmY48TMNYqtxrHVcbHUJ9WFRtYqyVRaDorGo0ahkSSdhoXGoVFQ0nFYKiipUlLoyydKiv07+0JIit5o4XNNwjLfzlzPiqn+UDzRcbHU8+WxZqNKs/FS58tKzUlU20bVdSF/aurqICWTgkt1v/Ou6FJraoLL2iAl5WMnNn2MTfLhfaw1hX4UTDLF4ApN7qTQmMwl64LU5mnIJCmoy0EVgtSlq3cpYi8QWMMish4Jo1rAdvAuKb6aKU0qFfMkt0KhlSyZrMijU5Oi1nFR61vVhf6k+kbHYXkZUJKU+p1ndR40DnUetbJWY2sU5Qr9PN36BFqQNDbT2KSxzdS59N041Z8WLzT1SrM00lk30TSN9KI71qwb6VU7VuNR5+1Yyy6/XM/aSkXodNGMVMVObQq6qEeKIWlaVwrmarqgRV2qaaK6JqpbRqkJirOg0JjCMi/492CySnmUq8jrv8r1tVr992CLNr9j8bbASqzD2kXE1eMgsh4ZsQVskOcda6pLiudLtcel4tLVTkxWm1ITZZNabYoq+ynC47DUaZjrSVgoWFKlTgsvFSyp8Xg5ulWpu5wmXOlkauSXu0hY7dd0bDlyanfFUGvsnY6t1mmY6zwd6Wmc6TyNtShLzdLo8n5etMcqrdNFN1Iw17QdqfGgNAmatpUWo0KzttKizVOMIbja6Fq0Ie9+YXVYoOjyYPJ+qjBFSUmKqX8DgF1f23XzOZSUd/OQkhRMYsZwZxBXj4vI2hCmEIEN8SSfzaTJWKFNspRDw/q1TU0TFSxp2RWKVVJU0tiafv1UjqiJGiXPodT1u8hZH8G6Flr9Tkw7d1V9aK1iqzLrR8E6TbzTzDsdW6Mums7TWJ1M0zRSUlDtUd8pX6rxqNoLLVOppszvAvy6OdGzKuir5amKkFTFTou2kHunrguy4JeHBrp8Gi73YN+fEVaB1b8LsX/3pa2iqv/82r4V2VfWziCwHh+RtUG8ExHYAHf5Yil78VJlETX+pFKqguLMlEZRjVU6b8YqQtKsq9T1i9Sfhloj60ehJEWZFp7USWo8T6t1ytHUre2bsJGp8RxeqY+v1TaCpGCmUpJMOlbSapKx81pd/1XneTsLj2o8qJPpRXei2qPO05Ek6bwbq7C8e4k2BTVdUN0Watsgb4LUmay9isnLdySu9mS/9vyYu6z1fMfJL9exXS5y71I/KsgBrbcdYbVZRNaWYBoReDy+XObj9Z2dKy6fqHplCl1QmAd1o6SXi7FOy4XOmok6N029UlC/iF2mkeWXztKSZt71YSU1nheL3wwtqR/x8rzOaSG/2t7aNOLaQJOCmaLyCFgwqZRLyve18KjSWi08Tws2HrVMheZdqUVbatkVaruY12S1UWpDXnPWmZTy+jPdFlv919blESxL3kdZP+WZ+tDquj6uErtu2FLE1XYgsrYMo1vA8Lzr1L26UOiSjv5pKrcTVWdB1phqizo7yccZXHSlfjv+jhov9K14rm+FpaK5TiwqyBQtaOSdGu/UqNPSk2p35T1D+OXo1konuxz5kvopxn4UKJpUKi9vimuPdTWa1bg0S4UWXmjhpf7YPtFZN9HLbqL/u3imV+1I/zT9SPOm1KwuNZuO84L3Nq83C/07DfMi/7V9ZyVdnh86z5d1eSTL2pT3k9W0smUjr2v5ss6R1bRME24Zwmr7EFlbjNEtYCCe68IXS4Vlo9HzpapvRaXKVEyDmjbo5fRIR2Wjl+2R/mAf6S9Gv9fLlBStUePd2mhWH1wyJW/6pUw5sPLs4PWF8JLUuSmaX0ZXlCu51CiPZjX99VeRtpoiPEtHWnipaRrpj+0TzVKlF+2xzpojndVHmjelpstKbRuV+ilCrUaw1katLF19Xl+PJu8DK+W9vVuT+kP5tFKTj52oppF36cMOVo0HRVxtr/D2q2DT+AcEDGD1TsMXLxXnjWLtql664lJK56XaJurri2P9fvFE37QT/a75WFMvdZYKLb1Vq05dP1UWLahQVGlBpQVF5VGpaLq20H21ID6u7aK9c+tjy/p3I1q//ipo4VGzVGjqhc7SWK/S+PJg0+fd+DKwXjVjzdtSy6ZQ1wW1TZT3U4NajVh5H1pptaB97bnoYyt0efcNoUkKdXc5iqW6kTeNvG6u1mFhK/D7YbsxkrUjbv5DYnQLeACpUzp7qZBcp6cj1R+VshTVHkU17UgX41L/q/hUsyeVJrHW75qP9WfV14rVHzROSz0LQSdhpNKiogWd2FidJ5WhVVJS40mL1UJ25d4J/bqtlVVwrY9ySVKjoMbD5fTgwkuddRPN0kjnaaznzbHOmommXaXzeqR5U6puo9o2KF2uuTKZX4+q/k2R/aiVZG0+TE+xdMU6qTprFOpO4eVMNl9Kbas0ncnrOu8Gom2YItwgomq3MJK1o/iHBjyQrpOfn6v86lzFvNPR86TRmak8D4rnUfNlpa+mJ/o/s0/0sjvSP9af6p/aJ5p6oT8m18xrNX59VGtkhUZWamSFxhZUmqnq30UYTSr7j/UF77eNcklSp6DO80fjhZp+9w3LVCrJ1KZwdYDmnplySZnLg8uj8uegy9qypMvjFhaLpGKeVJ53CrNGYVbL5kv5spYvFnktVtMSWBvG6/7uYSRrh7FmC/hw3raSBYVlreoPU6XyVHERVZ2Zls+k+flITV3on6pa0VwflXN9FGeS8o5KxzZTsqXKfiRLyqElTxpZoWiWI0yuRklReSekUt4V1eWo1vr6LbtaHJ88KCl/dApq/GpZfHdLXIXg8uT5sDuFpMb62FodNifvyiEu87Eb41IqZ0lhmRRnrcJsKasb+WKRpwnbNi9yZw3WxhBXu4vI2gO8IxH4MN7USr//SnZ0pKO6USo/UX0SFJdBM6vUjVz/qGd6tRjruKo170r9z/K7+pej55pWv9dxWOpZWOhbcarSgiZWqbQcQ1FBI8u7Wug8Kcn7qcROQUGNOiX3HFbyPsDy7hokaWzN5VTi2Go1IaqT6SjWajwo9u+CjCGpS0F1iGqCy70/zmFwWZcPpxNnpmKeA2v8TVI5TYqLpNGLpcJ0KdWN9PJC3rZKF9O13TQwevXYCKv9QGTtGYILeD+pbhRCkD3/Rke/n8j/5Ej6WmonQfVTqXlV6Rs3zcalJuXHSpOgZSoUlfQvym/6PbO/0sQ6pbC8HMUqlNdrSXmEa5VeRX+qVFTq3+6X37XoavoYi2qVzK4iKzR5VMuDRqHVSVzqlQcVllSGTjEkxWBqFPPuqzqTunyswmJuKi+kYuEqp67qPKmYdorLTmG6lM0WUtMqLRbyxZKpwQ0grPYPkbXHOHQP8A5Sp7RY5ncDTpc6+p2kf3Gk0Yt8IL9UBnWp1FLSy2U+YPNH1VyldWq80CjkHS8ch6WSZhpZrbGZJlaqXAWVXU31XYZX/2fnKR+ouR/VGispyLWwpLFadQqXh+AZhUaldVpaoTLkYyZKUnJT00V1XVBaFFJjirOg8twUa6k6d8Wlq5r2C9yXreK0ls0W8nk/PTifsxf3DSCw9pO9dhyq2/GvbU8QXcBbmClMJrLJRHYy0eL7n2r+aanZt4PqJ1J35Ko/7aQiqTqt9Z2n5zqplqpCp784/Uqflhf6KM70reI8TyPGC0W5TkOjiXm/AN76Pbqbwtr7j1bvSOzkmvU7NT1LlaZeaZFK/b79SLM00ov2RF81p3rZHOm8Gemfp080r0u9fDlRmhWyOqh6ERWXUjGVRi9d5TRp/LxVXHaK01rh1UyaL+RNo3R+wbqrR0ZU7ZW7DqdOZB0qYgt4gxBlZaFweiJ965nq757q1Z+N1ExMaSQtPnU1T5K8Shp9vFBRJD07nunpeK5vjS70tJzp29W5giX9WfW1JOlpmOlpnKlS0sg6JVmeWlR+p6GUDxW46NdZnaVKyYOep2MtUqWzbqJ/bp5qkUpddCP9v9lTzdpK3yyO9PzsRF0bpOcjhdpUTE3luVRe5IXtRy86xWVS9XyusGhl07n87GX/jsE2756BkatHQVztJSILb0Z0AbcwUzg6kh2N5X/ybTXPJmpOCk2/F9UembqxtPzYlSpXd5Rkx62KUavjo1qTUa1x0erj0UxV6HRc1DoulhqFVpNQq7RO0ZJKy6NHk7DUeRqrSYWSTOfdWMtU6EU90aIrtehKPZ9PtGgKzZeVlvNSaVrKalP1TR61qs76NVczV/WyU3nRKi5axa9fSYul0stX7IrhkRFVB4HIwv0QW8ANISpUpez0VHr2kdLJSLM/PVZzFJRK0+w7plRI7YmrG7m64yRVScVRq1h0Oj5aqox5zdTT8VxV7DSOjQpLCuYqQ46saVupCp0umpGK0Omb5USS9GI6UXJT2wXNpyN5E2TTqDgPKqam0EqjF67QSJPnnULtiouk6qupbLGUzqfyi6l8uWSt1SMirg4KkYV3Q2wBa1brtKpKNqrUfe9Tdcellp9WWp5G1U/60DqWurGrG0vdJC9k10mT97wek2LZKQRXWXYyc43LVk0bNSpbzetSyU0pmVIKcpeaZaHUBmkRZU0+rmJo8r6uRt/kPbZX565yllSdtarOlgqv5rJlncNqOsvTgW276WfwYBBXB4nIwochuoDear3WZCIbVfKTidLHJ/Jomn93rK4K6iqpPukPkzMydUd5D+vtUQ4jj1KqXKExpbLfMWmdj3cTllcHbi5m+YDN5YUUm/6dgS/zFGDxfCpb1tKyVnp13h9XsGPx+iMjqiAiCw+J4MLBM5PFKMWYR7cmR7Kqko8rpY8mak8qdeOoVJrasSmVphRN5q52bPIghUZKlSTPOwdNZT7MTTHPh78plq64TCpmnUKd8uFuLubSq4t8DMHZLE8BJiesHhlhhRuILDw8YgvQ5ciWmcnGI9nRkVQUUlkonYylEOTR1E2qvB7KTB6vXpOtc8ldllyh7mRtktokc5ddzOSzhSTlqKrzvrhYuL4ZxBXuQGRheEQXsMb6191+p6MW7PL0JU95JGrt6/yZl9xtQFThnogsPC6CC8AuIqzwHogsbB7hBWCbEFR4IEQWtgexBWCTiCs8MCIL24ngAvAYCCsMiMjC7iC8AHwIggqPjMjCbiK4ANwHYYUNIrKw+wguAOsIK2wJIgv7i/gC9hsxhS1HZOEwEFzAfiCssEOILBwmogvYDUQVdhiRBawjvoDNIKawh4gs4E2ILmAYRBUOAJEFfAgiDLgdEQUQWcCDIrpwqIgq4DVEFjAkogv7iqgC3orIArYBMYZtQTwBD4bIArYV4YWhEVTAoIgsYNcQX3hXxBSwEUQWsO+Isv1DNAE7gcgCDh0Rtn2IKGAvEFkA7ocY+3DEE3BQiCwAAIAB3BlZxYduAAAAAK8Lm34AAAAA+4jIAgAAGACRBQAAMAAiCwAAYABEFgAAwACILAAAgAEQWQAAAAMgsgAAAAZAZAEAAAyAyAIAABgAkQUAADAAIgsAAGAARBYAAMAAiCwAAIABEFkAAAADILIAAAAGQGQBAAAMgMgCAAAYQLHpB4DD8Kl912vVV2eYyS5PX/5x7XK98fIb59mNE3bjyjdu/uZtZH7X7e6xfb9le3fdj8zy9d94nVu2f8/78DvOf+N5esfv4cb5b73Pt933u97vrZf5ez9+3fwbue1H4dqPzvXr2y2PZf1H1G5cdvM269tb/zFdnb9++3z59fOv3X7te3rT5de3fWN7dv2y6/d/47Ib38fN+7y62O/c/s3v5er0+rbXtmHS2ivKa9dfXbZ+DZP06/+2/C/u/u8EDITIwqOoVevfhs9kYfWKHvJpC1KwHF0hXL3ir04Hk61d57XLVx+Xl/eDs+vXX79ckpvlMdy1yy/Pk14/b3X66jeC/PL81fX77UpS0OX1V5fnz7p8bKuvr12uW84z9fel/rHfuGz9s26eZ9e+fv3y12977fy1867ue21799j2a9vX3be/Ot+vXefmbS4vf22bvnadta9vXG43r5uv0V92/XP/1y2Zy8z7H42166zuau3ycHn6Kp5W5wVdv07+cegvu7y8P712+eqycNtll+ena1/ny5Pi5Xn58tifXt0+X766rP+8ulz5dLB0ebv128e120VLl9fPt833na+T+uvk+5Kk2J+3uk2U97dZbduvtn15On/kx95fR1I0Kcou/wlHMwVZf571XweF/i81Wv4qfu9/fSpgQEwXAgAADIDIAgAAGACRBQAAMAAiCwAAYABEFgAAwACILAAAgAEQWQAAAAMgsgAAAAZAZAEAAAyAyAIAABiAub921DTgwZnZf5e02PTjAIA1Y3f/15t+ENhfHLsQj2Xh7n+56QcBACtm9l83/Riw35guBAAAGACRBQAAMAAiC4/lbzf9AADgBl6XMCgWvgMAAAyAkSwAAIAB8O5CDM7Mfibp+5K+dPe/2fTjAbDb3vaaYmY/70+erS6/z3lm9gNJn0v6TX/+v3f3sztu+3l/3gt3/+lDfn/YH4xkYVBm9kNJz939x5Ke9y9iAPBe3vaaYmY/kvTrVfiY2Q/ue16/iV+6+4/7j7M7bvsflAPsx2uPCXgNkYWhfaar/xV+KYkXIwAf4m2vKX8labX/q7+X9JfvcJ4kfd/MfrYWXa9dz93/1t3P+vN+K+nph31L2FdEFob2VNKL/vSZpE829kgA7IOnevNrynow/VV//fued9Z//Z8k/Z2Zff+O6637TNKv3uP7wAFgTRaGdibpWX/6qaTnG3skAPbBmd7wmuLuvzSzL8zsz5XXbf32Hc77UtJqCvAXkn7o7n9783qr++rXav10bVQLuIaRLAztC0nrw+6/ecN1AeBt3vqa4u6frS2I/8/3Pe/G+q4/V56OvPW2/eL7n/VhBtyKkSwMyt1/ZWaf9f/jO3N3htUBvLe7XlPM7HN3/7GZPZX0d/3Vf9EvXr/veWdr7xr8sr+v267318prwX5gZpL0c3f/5XDfNXYVOyMFAAAYANOFAAAAAyCyAAAABkBkAQAADIDIAgAAGACRBQAAMAAiCwAAYABEFgAAwACILAAAgAH8f0UcoXuiHn+GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 612x388.8 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mask_512_BICEP = hp.read_map('bk14_mask_gal_n0512.fits', verbose=False)\n",
    "idx = np.where((mask_512_BICEP < 0) | (~np.isfinite(mask_512_BICEP)))\n",
    "mask_512_BICEP[idx] = 0\n",
    "hp.mollview(mask_512_BICEP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "edb43470",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/a6997a166010916ce04ada58b4872b86/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/1a77be884471968effb02a9e53c2a236/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/2233e1624d4c24ef8100824b8c088c68/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/ca73bb01bc98b80a5f50b6066aa9c8fe/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/2233e1624d4c24ef8100824b8c088c68/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/ca73bb01bc98b80a5f50b6066aa9c8fe/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/a6997a166010916ce04ada58b4872b86/contents\n",
      "No physical unit associated with file /Users/bhensley/.astropy/cache/download/url/1a77be884471968effb02a9e53c2a236/contents\n"
     ]
    }
   ],
   "source": [
    "nside = 512\n",
    "sky_0 = pysm3.Sky(nside=nside, preset_strings=[\"d9\", \"s4\", \"f1\", \"a1\", \"co1\"])\n",
    "sky_1 = pysm3.Sky(nside=nside, preset_strings=[\"d10\", \"s5\", \"f1\", \"a1\", \"co3\"])\n",
    "sky_2 = pysm3.Sky(nside=nside, preset_strings=[\"d12\", \"s7\", \"f1\", \"a2\", \"co3\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b0147ecb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n"
     ]
    }
   ],
   "source": [
    "map353_0 = sky_0.get_emission(353*u.GHz).to(u.uK_CMB, equivalencies=u.cmb_equivalencies(353*u.GHz))\n",
    "map353_1 = sky_1.get_emission(353*u.GHz).to(u.uK_CMB, equivalencies=u.cmb_equivalencies(353*u.GHz))\n",
    "map353_2 = sky_2.get_emission(353*u.GHz).to(u.uK_CMB, equivalencies=u.cmb_equivalencies(353*u.GHz))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9c4997a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ClBB(mask, map1):\n",
    "    b = nmt.NmtBin.from_nside_linear(hp.get_nside(mask), 20)\n",
    "    f_2 = nmt.NmtField(mask, map1, purify_b=True)\n",
    "    cl_22 = nmt.compute_full_master(f_2, f_2, b)\n",
    "    ell_arr = b.get_effective_ells()\n",
    "    return (ell_arr, cl_22[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b27f3dfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "(ells, clBB_353_0_BICEP) = ClBB(mask_512_BICEP, map353_0[1:,:])\n",
    "(ells, clBB_353_1_BICEP) = ClBB(mask_512_BICEP, map353_1[1:,:])\n",
    "(ells, clBB_353_2_BICEP) = ClBB(mask_512_BICEP, map353_2[1:,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9eb33599",
   "metadata": {},
   "outputs": [],
   "source": [
    "template_353 = hp.read_map('COM_CompMap_IQU-thermaldust-gnilc-varres_2048_R3.00.fits', field=(0,1,2))\n",
    "template_353 = hp.ud_grade(template_353, 512)\n",
    "template_353 *= 1.e6 # to uK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "62b6d8a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "(ells_512, clBB_353_tmp_BICEP) = ClBB(mask_512_BICEP, template_353[1:,:])\n",
    "Dlbb_tmp_BICEP = ells*(ells+1.)*clBB_353_tmp_BICEP/(2.*np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "733e77f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "npipe_353_A = hp.read_map('npipe6v20A_353_map.fits', field=(0,1,2))\n",
    "npipe_353_B = hp.read_map('npipe6v20B_353_map.fits', field=(0,1,2))\n",
    "npipe_353_A = hp.ud_grade(npipe_353_A, 512)\n",
    "npipe_353_A *= 1.e6 # to uK\n",
    "npipe_353_B = hp.ud_grade(npipe_353_B, 512)\n",
    "npipe_353_B *= 1.e6 # to uK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f85e3112",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ClBB_cross(mask, map1, map2):\n",
    "    b = nmt.NmtBin.from_nside_linear(hp.get_nside(mask), 20)\n",
    "    f_2_A = nmt.NmtField(mask, map1, purify_b=True)\n",
    "    f_2_B = nmt.NmtField(mask, map2, purify_b=True)\n",
    "    cl_22 = nmt.compute_full_master(f_2_A, f_2_B, b)\n",
    "    ell_arr = b.get_effective_ells()\n",
    "    return (ell_arr, cl_22[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ec0e1461",
   "metadata": {},
   "outputs": [],
   "source": [
    "(ells_npipe, clBB_353_npipe_BICEP) = ClBB_cross(mask_512_BICEP, npipe_353_A[1:,:], npipe_353_B[1:,:])\n",
    "Dlbb_npipe_BICEP = ells_npipe*(ells_npipe+1.)*clBB_353_npipe_BICEP/(2.*np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a5f11ca1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAADRCAYAAAAZrWyfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABuYUlEQVR4nO3dd3iUVfbA8e+dyaT3RnqjhE4IodiQakcsgN1dFcGy7uqqoGtZu4L6W3dXUcHeaSqLnSIWmvTeQ0gP6b3O3N8f804YQiBtJjMJ9/M8eTLzzjvvnISEObn33HOFlBJFURRFURTFfnSODkBRFEVRFKW7UwmXoiiKoiiKnamES1EURVEUxc5UwqUoiqIoimJnLo4OoC2EEAmAPzABWCKlTHVsRIqiKIqiKC3raiNcyUAqsBKY4uBYFEVRFEVRWsVpEi4hxAQhxIomx6Zox2cBSCmXSClL0Ea4HBCmoiiKoihKmzlNwiWlXGl9Xwgxxep4iRBignbckmwVdXqQiqIoiqIo7eDMNVzDgYXa7VQgWQgBMFu7vwI1yqUoiqIoShfgzAmXf5P7Qdpo18pmzgXA3d1d6vX6xvvBwcGEhITYJ7oOys/Pd5rYSktL8fPz65ave+DAAerr6xk4cKBdrm/rr8GZfi4U56J+NtqnO37f8vPzKSgoAKCqqqpWSunu4JCUVnDmhKsECGzLE/R6PZWVlfaJxsZSUlLYvHmzo8MAYMaMGcyfP79bvu727du58cYb7fa9tvXX4Ew/F4pzUT8b7dPdv29CCKOjY1Bax2lquJqxiROjXAmYpxDPyGg0MmPGDJYvX27PuLqdSZMmddvXTUpKwtPT027Xd9T3TlGUs9vy5cuZMWMGgL6lcxXn4DQJl1Ykn2JVLL8ESLAUyzctqm+OXq9n/vz56k2wjbpzwrVy5UrKysrsdn31s6YoiiNMmjTJMrquRri6CKeZUtQSrCVNjs1tyzWCg4NtGpM9aX+ZKHb23HPPoS226BLUz4VyOupno33Ogu9bvqMDUFpHSCkdHYPN9O7dW44dO5ZJkyapkQcFgDFjxgCwZs0ah8ahKIpiS8uXL2f58uUsWLDgsJSyt6PjUVrWrRKulJQU2Z2LI5W2UwmXoijdmRBii5QyxdFxKC1zmhouRVEURVGU7qrFGi4hhC9tbM+gKZJS2q9auRmlpaXMmDFDTSkqiqIo3ZplShHo/CaKSru0OKUohHgRc4uGtlYep0gpH21vYO2hphSVpg4cOABAYmKigyNRFEWxPTWl2HW0ZpXiZinll229sBCiPaNiimJTKtFSFEVRnEGLCZeUcmlL5wgh4qSUadrta4BAKeWCjoenKB1jaYKrppgVRVEUR+pwHy4hxMPATCHECmALkAIECiGGSSnv7uj120LVcClNvfrqq4BKuBRF6V5UDVfX0+G2EEKIO6WUC4QQ4wF/y4iY5bgtgmwtVcOlNKXaQiiK0p2pGq6uwxZtIQoBpJSrgCNWx480f7qiKIqiKMrZxRZb+yQIIeK02ylCiBLtdjKw2gbXVxRFURRF6dJskXBdD/TkRNsIy9DmMOAVG1y/1VQNl6IoinI2UDVcXY8tariGSim3tfa4PakaLqWpjIwMAKKjox0ciaIoiu2pGq6uwxYjXKfL2LrPJo1Kl6USLUVRFMUZ2CLhWiKEsBTIC8yJlgDiAbWDueJQCxcuBOC6665zcCSKoijK2cwWCddM7XOqlPKo5aAQ4k4bXFtROuTNN98EVMKlKIqiOFaHEy6tHQRCiHghRJJ2bLvqNK8oiqIoimJmiz5cAEgpj2qJ1nYhxFvaptedyrJK0bKdi6IoiqJ0R8uXL2fGjBmgVil2GR1epdh4IfPo1l2Ya7eWOGKES61SVJpSneYVRenO1CrFrsMWeyk+BEwEtgIvWW1iPU5KqRqfKoqiKIpy1rNF0fxFwBztdoIQIgHzKsVZqE7zioMtWbLE0SEoiqIoik0SrtmnaXxaZINrK0qHBAcHOzoERVEURbFJ0fwEIcSbQohx0LhacVxnd5lXlOZ88MEHfPDBB44OQ1EURTnL2SLh2go8YqnX0lYrrrYkYJ1JrVJUmlIJl6Io3ZFapdj1nHGVohDCF5gGrLQqhh8PSEuCdbrieCHEeEuPrs6iVikqTalVioqidGdqlWLX0VIN1z+0zxcJIQqllHdLKVcJIQ5xYtse/9M8V2XdiqIoiqIotJxwrbDqJO+ntYBYAJRandNTCBFnGQHTzo0Deto4VkVRFKUT1NU0oNfr0Bts1hvb7jL2FRHR2x+9S9eJWTm7tJRw+QshrpFSfimlLAVeEUJci9XolZTyZSHET0IICaRiTrSklPJi+4WtKIqi2MvX/7eNwAgvJvy5v6NDaZWygmr+9+/tJF8SyzlXOcff+nlHy3D3NuAX4uHoUBQn0VLCtRJzDVcjKeVSIURJk2MXCSGGAimYu8x3au2WopzOd9995+gQFKVLKSuoJj+9nLLCaqRJInTC0SG1yNhgAmDn6gwGj43Cy8/NwRHBkjnmeuJ73+r09WOKkzrj2KuUsrS5LXqaS6iklNuklAtUsqU4E09PTzw9PR0dhqJ0Gel7CgGorWygILPCwdG0TUOdiS3fpTk6DEVpVqsmu4UQ01txzkMdD0dRbGvevHnMmzfP0WEoSpdxbE8RHj4GADL2d63+1Z5+ruz5LZvS/GpHh2IXhVkV1FbVOzoMpZ1aW104+0wPanVdc850jqI4wqJFi1i0aJGjw1CULsFYbyLzQDE9k0MJCPMka3+xo0Nqk6ETY9DpBX8sT3V0KHbxxbN/8NX/qZ7iXVVrE65VzY1gCSHihBA/YU621E+BoihKF5Z9uISGWiOxA4KI6htI9uGSxvqorsDLz43B46I4uCmvy02HtlZhN/26zgatSriklHdJKV+xnloUQjyMeVXiFillL2C8nWJsNdVpXlEUpf2O7SlE5yKITAwgKjGAhjoTeUfLHB1Wmwy9KBY3Dxc2Ljvi6FDsSnWa73ra1LBESvmOEOJOrfHpNKCnlPJR7bHSMz/b/vz8/Jg/fz6TJk1ydCiK4hTy0srY+tMxGuqMjg5F6QLSdxcS2ScAg5ueiD7+CAGZXayOy93LwNCLYkjbVUjO4RJHh2M3kyZNYv78+XByX0zFibW2aD7OcltbtbhUSjlcSnm0uXMURXGshnoj6748zNI5m1n/5RGWzNlCcW6lo8NSnFhZQTXFuVXEDggCzIlLSIwPmQe6Vh0XwOBx0Xj6ubL+6yOcafs6RelMrR3hmml9R0r5iBBiurbXYrPnKIozWLNmzVm3j2JuaimLnt/Etp/S6XduOBffOZDKkloWv7iZg3/kOjo8xUlZ2kHEDAhsPBbVN4C81DLqahocFVa7GFz1DL8sjpzDpRzbXejocBQFaMMqRSGE0foDmA8Ua/dNwCz7hakoSksa6oysXXqYL1/eQn2tkUl/HcLYW/rRa1go1z0+nOAob1a8t5efP9mvphiVUxzbU4RvsDv+PU70rYtKDMRkkuQc7nqzVv3Oj8A3xIMNy1KRJjXKpTheaxOu+UBgk4+AJvdPaZCqKI72yiuv8Morrzg6DLvLOVLKwuc3sX1FOv3Pj+CGJ0cS0z+o8XHvAHeu+vtQki+OZe/v2SyZs1lNMSqNLO0gYgYEIcSJzvJhvfzQuYguV8cFoNfrGHllPIWZFRzakufocBSl1QnXHK3r/Ok+SlB9uBQn9M033/DNN984Ogy7qa8z8vuSQ3z5yhaM9Sau/FsSY27qi6vHqbt26fQ6zrm6J1f8ZQiVJXUsenEzBzaqKUbl5HYQ1gyuesIT/LpkHRdA72E9CIryZuOy1C7V3kLpnlrbFuLomR4XQowH1JitonSi7MMlLHzuD3aszGDgBZFc/+QIovsFtvi82IFBXPf4cEKivVn5/l5+/nifmmI8y1m3g2gqqm8ABRkVVFfUOSCyjhE6wajJCZQV1LBvbbajw1HOcm1qC2EhhPhJCHFICPGmEOIaYBOQbNvQFEVpTn2dkd8WHeSrV7diMkom35/EhTcm4ure0l70J3gHuHPVA0NJviSWvWtz1BTjWS59dyGRvf0xuOlPeSyqrzmJzzpQ0slR2UbswCDCe/mx6ds06tUfFooDtSvhklJeJKXsDSwBLgKOAgm2DExRlFNlHyrmi2f/YOfqTAaNjuT6J0Y0viG2lU6v45yrenLFfUOoLFVTjGerskJzO4iYJtOJFqGxPhjc9V2yjgtACME5V/WkqqyOnaszHB2OchZrV8JlIaVcpXWhD8KcdCmKU/Hw8MDDw8PRYXRYfa2RXxce5KtXt4GUXPXAUEbf0LZRrdOJHRDEdY+NaJxiXP3xPjUScBZJ32NOpGIHNp9w6fQ6Inr7O3UdV0uttsJ7+RM3KIhtP6VTU6k2f1Yco71TiklNenABdMpvoxBihhBiQnOPmaqqMNV1vToDxX6+//57vv/+e0eH0SFZB4r54tmN7Po5k0Fjo7j+iZHN1tp0hHeAG1c9MJRhl8ayb10OS17aTFGOmmI8GxzbXYhP0MntIJqKSgyg9Hg15UU1nRiZbY2c3JPa6ga2/XTM0aEoZ6n2/nk8EbhOCCGBlcARYBiw2laBncFmTjN9WZd6lIPDUnDr1w+PIUPwGDwYj6QhGKKiTlrqrCi2Ik0m6tLSqN6xE2NJCcLVgHB1RefqirD+MFjfN5z6uOVDd+JvoLqaBtZ/dYTdv2ThG+LB1Q8OJaK3bRMtazq9jlGTexLRy58V7+9l8YubuPDGRBJHhiHr6pDV1Zhqa5E1NZhqapE11ebPtTWYamq04zXImlpMtTXI6hpMtTXovLzwnzwZQ2Sk3WJX2sfSDqLvyLAz/h9pmbbO3F9Mv3PDOys8mwqO8qbPiB7sXJ3J4LHRePm7OTok5SzTroRLSvky8LIQwg+YgDkBSxFCbMKcgC2UUm5vyzW1UavZUsqJVsemACVAspRybkvXMMREE/inW6nevoOSJUso/vhjAPQBAeYEbMhgPIYMwX3QIPQ+Pm0JT+minn32WQCeeOIJm1zPWFpK9c5dVO/YYf7YuRNTqQ2bQrq4IFxdKQlIZE/stdQY/Ikt30Jizh/U/kNHmpawNSZ1huYSt5OTPvT6U5IgaUmUrO/XnEicRhhd2dljMqs+MLH3xXfpc2gRelPbp2KEqyuyvp6C19/Ae+xYAm++Cc9Ro9QfQE6isR3EaaYTLYIivPDwMZB5oKjLJlwAI65I4PDm42z6Lo0xNyY6OhzlLNOhAhBtw+ql2gdWCdh1wPY2XmulEGK25b6WbFmOJwghJkgpV57pGnpfX0Ifesh8vYYGag8donrHzsY3xgrLFi9C4NozAY/BJ0bB3Hr1Qrh0vB5GcS6rVq0C2pdwSaOR2sOHqd6uJVfbt1OXmmp+UAjcevfG96KL8EgagseQIbj06GEeCbL6MDXerjd/rj/z43U1DezM60FqeSje+ioucP+DQM/jyLqwE8+rqcVUVk6D5Rr1VtfXPjCdpueQwYDOzQ3h7m7+7OGOzs3dfN/LC31wMDo3Nzzc3RnrlsW+anf2cw6V8cM4P7EI/wAXdB7uCDd3hLsbOnd3dO7m5ws3t8bbOu2+0Omoz86m+IuFlCxeTPqqVbj27EnATTfid+Vk9N5e7f2n7XTVDdWklqQS4R1BgLv9Rho705naQVgTOvM5mfuLkVI6b8LcQlh+IR4MOD+CPb9lkzQhGv/Q00+jKoqttZhhCCGSWjtaZZ2ACSHGSSk7MsU4HFio3U7F3HZiJeaELkgIsVJruNp83C4uuPfrh3u/fgRcfx0AxrIyqneZRydqduykYvVqSr/80ny+hwceAwfiMWQw7kOG4DF4CIYeoR0IX+lqGgoLzQn69u3mn5FduzBVVQEnRkn9rpyER1IS7gMHovf2Pun51ZUVGHx9MRgM7Xr9jH1F/PzxfsorahgyIZqRVyZgcL3ijM+pqKsgpzKH3Moccitzya3MJacyh7yyHArKsikqP44wmkgI6UtSzAiGR44iOTQZT0Pr3mgigN57C1n5/l5+2ufBhTcm0ndU20Y4DBERhP79AYLvvYey776n+JNPyHvmWfL/71/4XX01ATfegFt8fJuuaW91xjoOFR9iT+EedhfsZk/hHo6UHMEojeiEjuTQZMbFjGNczDgivbvuVOmZ2kE0FZUYwOHNxynOrSIwvOskyk0NuyyOfetz+GP5US66Y4Cjw1HOIq0Z0rleCJHajmtfRMdquvyb3A8CONPUYk32Hl7/ywQ+2FACwJ1/up4Z9/69sS5G7+uL93nn4X3eeWjXoj4jQ5seMo+EFX74EdSbp05cwsPNI2DadKT7gAHo3N078CUpzkLW1VFz4MCJ0asdO6jP0JaMu7jgnpiI39VXN45eGaKjG/+ql1JyMK+CI798Rei+D4ib/hHB/n7sWDqHQUfmc1uPT+kTE8GogFL6hfkSGd/vpNqspuqqG1j75WH2/paNfw9PrnloGOE9/ag31pNZntOYRDX3uaK+4qRruQgXQj1DCfMKo39UMmGeYeiEji15W/jowKe8t+9DXIQLA4IHMCJsBMPDhpMUmoSHy+lXcsb0N69i/OndPaz6YB9ZB4oZfX1iq96krenc3PC/+ir8rppMzY4dFH3yKcVffEHxxx/jdf75BNx0I96jRyP0bbtuR9Wb6jlScoQ9BXvYU2j+OFh8kAaTecNmfzd/BgQN4MKoC+kT2IdDxYdYnb6auZvmMnfTXPoG9mVctDn56hPQx3lHf5qwtIPof35Eq863ruPqygmXl58bQ8ZHs+X7Ywy9KIaQ6K5XXjJ//nzmz59vuRvsyFiU1mtNwrUQ8xRhW0hOjE61VwnmPRpbzejixZDxU9j8+l1UV5bjOjeaN585wPro6QyO8Ga8biuxQ8cRFGr+i1QIgWtMDK4xMfhNmgSAqbaWmr17qdm5szERK//xR/MLuLjg3qcPHklDcNcSMde4uC7zH+zZrD4396TkqmbPHmRtLQAuoaF4JCURcP315n/b/v3RNdNK4ti+LRR++xTPVE1le1UQ43S7ecF1DwUZhwj2T8G31zlsriqn0uTBJxuO0Zu38NRvpu7xY7jpdOz4eTENRhORwyfTw9eNopoi9u44xt4vi2koh/pBx9mWuJXv9meTuyWX/Op8ZJMNHALcAgjzCiPaJ5oRYSMI9wonzCuMMK8wwr3CCfYIRq9rPmGpqq9ie/52NuVu4o/cP3hv93ss2LUAg87AoOBBjAgfwYiwEQwOGYyb/uSCYi9/Nybfn8Smb9PY/H0aeWnlXHLnQAIj2v7GK4TAIymJyKQkesyeRfGiRZR8sZDMu+/BEB1NwA034H/tNej9/Np87ZYYTUaOlh5tTKz2FO7hQNEBao3mnwUfgw/9g/pzS/9bGBg0kAHBA4jwijjpd/ySuEu4b+h9HCs7xs/pP7M6YzVv7niTeTvmEekdaR75ih7H0NChp/23cAYttYNoyi/EA58gdzL3FzF4bJQ9Q7O7oRfFsvvXLDZ8ncqk+4Y4Opw2mzFjBjNmzABACFHg4HCUVhKypQYmnUgIscJSNG9Vw7VECDEDSG2phislJUVu3rwZgKqKUnb9+B4ba2L4sTAUU95evjfMYk3/Zxkz7a8U5KZz+Nv/YBpyA336DiTY+/QrVhoKCqjeudP8hr1zJzU7dzZONen8/MyjYFotmMegQej9/W3y/VDax1RTQ83evVRv38EtLzyPsbSU14LMfwQKV1fcBwwwj1paaq/CTl6hJaXkSH4l2/YdIGHTM9QNuZlzJk4l8/AuDJ9MZmns4wQPmsg58YFEBzWfcJTVVrBl188cy9mJb99+5FTmsG/du+SLBnbrg3DX1XJO+uX0O34OxR65/NzzM0r9c09JoKzvh3mFnXEkqq0q6yvZmre1MQHbV7QPkzThpndjSMgQhocNZ0TYCAYFD8KgPzFNmrGviBXv7aG+1sjo6xNtUkQt6+spX7GCok8/o3rLFoS7O36TJhFw8024J7avuNkkTaSXpZ9Irgr2sK9oH9UN1QB4uHjQP6g/A4IGmD+CBxDtE41OtL1bTkF1AWsy1rA6fTUbcjZQb6onwC2AC6MvZHzMeEaFj8LdxblGx7+dt5PCrApuee6cVv/RuPrjfaRuy+f2Vy5Ap3OePzSLcir5/OmNXDR9AL1TerTqOVt/PMb6r45w41MjCQiz/YjdG3eZJ3jufWucXa8phNgipUyx2YsoduM0VeJagpUihJgipVyiJVqzLD23Wkq2AEpLS5kxYwaTJk1i0qRJjLz2AUYCfwVqalLYv7s3faP6AJBzYDMj0t/hmsPRbJe5XOx7jJmG7wi4ai7xvfsjTabGaSCX4GB8xo3DZ5z5h1wajdQeOWIeKdESsYLff2/svucaF3dSLZh7Yh9EO+t6lDOTUlKfmXlSYXvN/v3QYJ4O+k9cvDYlbE6w3Pv2Na/ca0Z9XS07/3s9P1T2ZkHVhbhSz0r3A2SUHwcgMmEAPHmYe6ymB8vrytl2fBtb8raQWpraON1XWmu1cnEdCATBQUH46f05t3AIiTuH417nhgzayOT7ruQe/0vIeO1ajngPoGDwnQwJ92dgDze8vLybhmkzXgYvLoi6gAuiLgCgrK6MrXlb+SP3DzblbuKN7W/wBm/g4eJBUkgSI8LNU5ADEgdw3eMjWPHuHlZ/tI/sg8WMvqHtU4zWhMGA72WX4XvZZdTs20fRp59S+r//UbJ4MZ4pKQTcfDM+E8Y3u7BFSklhTSFZFVlklWexv2g/ewr3sLdwb+OUq5vejb6Bfbmm9zWNCVasb6zNRqCCPYKZ0mcKU/pMobK+kt+zfmd1+mpWHVvF14e/xsPFg/MizmNczDhGR43Gz832o3dt0dp2EE1F9Q1g39oc8tPL6RHXtBWjmUmaOFxyGA+9B9G+0bYK2eaCIs2/W7XVDQ6OpH2WL1/O8uXLARz7w6S0mlONcHWU9QhXa1SUl7A7r5Zd2ZXI/d9yWc4bNNy+grjoaP5Y/AoRe9/h9Z5v0ys2muRgE4lRwXj5+Dd7LWNFBTW7d59YFbljB8bCQgCEm9uJURWtNUXTURWldYwVldTs3nXS9KCxyDw1Ijw9tYUPQ/AYmoTH4MG4BJ9a3iClJK2wig2phUStexyjwZsx984DYPeLF3LYdxTVw+/lnIQgYoM8T/p3KqopYlveNjbnbWZL3hYOFB/AJE246FyI841rHJUK9w6nh2ePxtuhHqGYauH3JYfZvy6HgHAvxt/aj6BoT1xcXDA2NLDjtWtYVdOXNyouxEADO92mszF2BmNufx5pMnF4x2/E9BuOm3vnrKwqqSlhc97mxgTscMlhADxdPEnukczw0BGE7uvPsTWVBIR5cfGdAwiKsF2CaCwpoWTpUoo/+5z6rCxkSCDFl47k8Og4julLyK7IJqsii5zKnMYpQQCDzkBiQCIDgs2JVf+g/vT074mLrvP/vqw31rMpbxOr01fzc/rPHK8+jl7oSQlLaaz7CvMKO+3zpZRUN1RTVldGaW0pZXVl5o/asuZv15XRYGogMSCRgcEDGRwymF7+vU752jP2F/G/17Zz2T2DiR/c+hKgqrI63p/1O6OuSmDYJXGNMaaXp7MxZyMbczayKXcTxbXmPthDQ4dyVa+ruDjuYrwM9qv7as8I17HdhXzz+g6unT2MsHjb5yxqhEtpqlslXIP6J8nNmzfh5tnx0aQdqxdRtXURf6+7i5yyWh53+Zgb9KsxzU7Hx9OdQ9t+wSQhbvD5uLmc+leylJL6rGxqdu5orAWr2bvXvGQfcAkJOakWzGPAAHReXbcQ1R6kyUTd0aMnJVe1hw41tjxwjY/HIympcfSqudYejz76KAAvvvgiABvee5jKjF3cUX0fAP/yeI/g4BAuuPetZmPIq8xjS96Wxo8jpUcAGqfdUnqkMKzHMAaFDDrjdF/argLWfLKfqvJ6hl4Uw/DL43AxND+6UlBRy56jmbj+8RZuvceQPPpyMg/vIuqT83msYQa7wiYzMkxwkW4zcedcTUh4TBu+q+1XWF3IprxNbMoxT0GmlaUB0LsiiQsP3oCL0ZX+kwO5cEJSm6blTNJEQXVBYwJlnUhlV2STU5bFwEO1XLJZMiRNUq+HLQPd2XVhNPTvRaR3JBHeEUR6RxLuFU6cb9xJU6DOwiRN7CnYw+qM1axOX01qqXktkiUprKivoKyujPLa8pOSqQZ5+hEYgcDH1QdfV1983XzxdfVFINhXtI+S2hIA3PXu9A/qz8DggQwKGcTg4MGkfl/Jrl8ymf7q6DaPTH7+zEYM3uA2+bg5ycrdSG6leQ/OUM9QRoWPYkTYCAqqC/j68NeklaXh4eLBxNiJXNXrKob1GNauadszUQmXSri6gm6VcMWEJMpZ17yBW6Ak6fxexPQPIiTWp8O1BsfLa0jfsYbqzN1ccL25z9eOORfhXpnFFcaXSQzz4TaP34kOC2HEFdNPe52TVsZpRfn16enmB3U63HomoPP0Ar3evFJLrzdPa7roETq91XGddl+H0Ls03hcuetDpEXod6F3Mn3VWz2m8ptU1XPQnXavxOSddS7vGKfd15gRHZ772yTG7mD83d1+vxygFxXnVHM+oJP9YOXnHyqkpq8XdpQG3+jIMpbnoc9NwLT+Oa20p7i4N+PeJxm9wIp5Dz1wrJ6Uko6iaDamFPH3TaLzr8/nlQAnuri5s+OgJRMF+Dp/7MqN6BpMQ7HXS6sPMiky25G1hc655BCuzIhMwT78NDR3KsB7DSOmRwoCgAa16U6+prGft4kPs35BLYIQX4//Uj9DY5qdizqSirJiD65axriaetcfd8M/6mTfFS2wc/REjx00mdfdGjv/2HqVJM0nsnXjKyJw9HK86zqbcTWzK3cSOY3vot208kWV9ONpjO6bzs0mJSmZE2Aji/eLNCVXliYSqaVJV36SpaqB7YGMiFeEVYf7sHUFEgQn3ZWuoXLYcU1UV7oMHE3jzTfhccgm600wVO6ujpUdZnb6a1RmrSS9LP5E4WSVPLd32Nng3m7xYfpZ35e9iV4H5Y1/hPupM5j/4btjxOC7eEr9rSxkUMoiBwQPxdT39z2VpbSmbcjexMWcjpT+7E5ExgPeHP4KPhzcjwkYwMmwkI8NHEusbe0o95I78HXx9+Gt+SPuByvpKIr0jmdxrMpN7TibCu3UrJFtyNiZclinFBQsWHJZS9rbZiyh2060SriGDhsr5zy8lfU8hx9PLQYKblwvR/QKJ6R9ETP9Am23nkJtxhCNHU/m9KpqdmSX8I/NuKtx6MPLRHwDY/O8byPPuR+3Q2xkc5U9CkCc6/an/MTYUFTUW4tfs24+srUWaTNDQYP5sNCKNRqTJCEYTmIzIhhP3pdFoPqeZ5zQeNzp2I2KT0FHl2YMyn1jKfWIo84mlwjsSqTMnLC71lfhUZOBWU0ydqy91bv7UegRSrz91xEjoBJ6+rnj5ueLp54aXvxtevgY8/Vzx8ncnc/9qona+yNV1sygWvpgW3kuIroIfNu4jJDjkpGtJKUktTTUnWNoU4fEqc72Wv5s/yaHJpISZR7D6BPRp85TU0Z0FrPl0P9Xl9SRfHMPwy+LRG2zzl73JaCTj8C4CI3vi4+3D5m8XMPCPR7mg9jXyCWCq+yb+5LaGwNs+IyLMNm9qLckpz2HlVzsoXWegwquQ73q9Q7FnDgJxymrLIPegEwmVVVIV6R1JuHd4i4sDjBUVlH69jOJPP6Xu6FH0QUH4T5tKwHXXYQg7/RTd2azeWM/BkoNsP7KHkveC2Z/4C2sCv2x8PM43jkHBgxgUMohBwYMorS1lY655mnBf4T4kEg8XD8YYryB23fkMvSOQUSmDWz1aVd1Qzap0c03bxpyNCAQjwkcwuedkJsRO6NCCEHsnXFJKdhfsZl32OibGTSTBr9nd5U6iRriUptqVcAkhpksp37FDPB1iXcNVXVFHxr4i0vcUkb63iOoy8192QZHexPQPJHpAIBE9/W34BmiirKwI/4BgpMnEnjnj+Lm6F6/WXY3AxDq3v3Ew7mYuvO0ZAHKOHSAsuvcZ+zPZipTSPA2nJWCywWhO3CxJmSWRM5rA2GB13+pxo3UyZ0IaG8BklfAZjWAyYao3Ul5uoqBYUFSqo7BUT1G5C0aT+a9eF52JQK9aAtxrCPCoItCtCk99DcJoROftbV5sMHAQem8vjA0mqsrqqCyppaq0jsrSWipLaqksq6Mov4qiwmrqy2oQDc18DwW4+Rj415L70bvo+M9Tn+Dpa6DMpYh0Uyr7a3azrWITuaYsEJIQj5DG6cFhPYaR4J/Q7mmPmsp6fl90iAMbcwmK9GLcre0b1Wqr+vo6DuZXszOzFLFrMYNylhD/8C94urmyY9UX+ITEkDD4XLvHkbG/iBXv7aWuuh6fsZUUxqYS5h3WmFyFe7WcULWWNJmoXL+e4k8+Ne8kodPhGhuLITISQ1QkrpGR2u0oDJGR6AMCzvrayd2/ZvHLZwe48amRuASa2F2wm90Fu9lZsJNd+bsorClsPNdF58KQkCGNI1iDggdhqhO8+/dfGXZpHCOvbDnxaE52RTbLjixj2eFlZFVk4W3w5uK4i7mq11UMCRnS5n8jeyVc1Q3VfH/0exYeWMjewr0A6IWeqX2mck/SPWfccUAlXEpT7a0inQk4XcLVdJVin+Fh9BkeZl7FlFWhJV+F7FidwbYV6bi46ohMDGgc/fIL9Wj3f8Y6vQ7/AK31gE7HwEfX0M8kuTi/gt1Hszm2+UIMoeZR36K8DMLfH8EccRu7o25gWLg7F7rsIi55IgHBPcyjVwibJWNCCLBM99nkimZSSsoLazh+rJzjx8o4fqyM/GPl1NWYR9T0Bh0h0d4MSPYlNNaH0Dhf/EM9Ea2c4tW76PAJdMcn0J26ujpcXV0pyE6jbsGV/F/d1SwxXUhCSBVzPT6BEQ8QFTOUylItOSuppbykBpMwUV9bz44NqbjUWUY3g4llDLGMAZ3E3ceAr78HXoVueGW5UeynZ79/rnkkzd8NLz83PLwNrYo7dXs+v3x2gJqKelIujyPl0jj0LvZPqgEMBlcGRLgyIMIPRjwIPAiYR8P81j5HdoMP7xyZz4MT+xB0hjYoHRXdN5DrHhvOivf2kvWTZMjoi7nQTvvWCZ2usZlxXWYmpV9+Se2hw9RnZTW7z6Xw9MQ1MgJDxMmJmCEyEteoSHR+ft0+ITu2uxCfIHf8e5innc+JOIdzIs4BzL/TuZW57C7cjaeLJ0NDh566K4EHhMb5krm/qN0JV4R3BHcPuZuZg2eyJW8LXx/+mu+OfsfSQ0uJ841jTPQYBocMZlDwoDMuKrCXtNI0Fh5YyLIjyyivK6eXfy8eG/kY50Wcx4d7P2TRwUV8m/otMwbP4MZ+N+Kq7/zpbLVKsetp7wjXYWAW5uakRZi34UmQUj5q0+jaqLWrFOtqGsg6WELGnkKO7S2iLN/cl8c32J2Y/kFE9w8kqm8Aru72WdVUVlLAvhUf8lt9Iqvy/fA5volFhqf5LeV1LrjiFvau+4b+P93ETQ1PslM/kHP1e3na9Do5l77L0JFj2L/xB1j5NAsC/k6RRyyJDfuZWPYVgVe9RHzPRA7v+J2CDZ+xLfpPmDyC6FFzlOjSzfSaeCdBQcFkH93H8cPbKI8ajaubB571hXjW5hPVNwU3VzeqK0qpq6nC1ScEV4MevZZoVJbUkpdWRn66JcEqp6bCXHuj0wuCIr0JjdOSq1hfAsObn0ZtjdzCEjakV7LxcB537b2Z1KAxjP3LPKTJyOZ/X09hz2uIHX45iT1O1OhVN1SzK39XY4H7jvwdHJp3CIALZ13IsOAUhngNo5dbX9xrvaksqaOqtNY8claq3S6po6by1E2adTqBp2Ua088VLz83vPwt991w9zawY1UGhzblERTlzfhb+xES4zwdrEuLCvhg9Q7+u7WWENda/q/PPlKm/B2Dq/16Q5lMkpXv7+XQpjzuen1MpyWe1ozl5dRnZ1OfmUl9Vhb1WVnUZZo/12dmYqo4uVO/zsurSSIWgatVUqb3tf9IpT0Z602889Bv9B0Z1qEkeMPXR9j6UzrTX70AVw/b/D9ZWV/JT2k/8b8j/2NH/o7Gur5Qz1AGBw9uLPjvH9T/lCTQFiNcDaYG1mSsYeGBhWzI2YCLzoWJMRO5ru91JIcmn5SIHyk5wqubX+W3rN+I9I7kgWEPcFHsRSedo0a4lKba+5tSLKX8UggxFOipHWuxT5azcHV3IX5wcONy6JLjVWTsNU897t+Yy+5fs9DpBGE9/YgZYK7/Co7ybvXITEt8/YMZOfVBRgIPATXVyezfPYT+CeaOx/7hPVkfPZ1zg4bSR9+DgPIqMo4Px99P6wgtdDToXKkxCgoq6oisLiCkej+VNeal8SXpe0jKXsRDaSPIlCHcoF/FFMO7HBx2DUFBwWRu/JoR+18iueYtivDlDv23PGH4lGN37ic2MpwdX77MkMMfMrlyHkFGV84x5eHd4Em1ydL4X+Ljks0hrygqIlwY5LKeFLGGAQ99g7tBz5blb3Jw6x/8EPMgbi46+pf9RkhdJuf/6VkA9q77ltKCLPJjr8BVryOoYj+e1DJg1MUA7HnxQjKqDNxf/wC+7i5c4D+awLhh2peuZ/gDiwHzPoLrctY2Jli7CnbRYGpAIOgb2JcpfaYw7PNhDA0dSpBH67ppg/lNqbLMehqzThs5MydmpfnV5BwuPSUx0+kEw6+IZ9glsQ5JLs7ELzCYv00Zz+Wjy/n1i1cZeeA1dm4/j6QRF9rtNXU60a5O9Lak9/FBn5h42uapxtJScxKWlUW9JRHLyqI+I4PKDRuQWoNjC52vr3l3ivh4XONicYuP127HNbs7gbPJPlJCQ62RmFZ2lz+dqL4BbPnhGNmHSohrQ1uJM/EyeHF176u5uvfV1BnrOFB0gJ0FO9mZv5NdBbtYmW5+i9ELPb0DejMoeBCDQwYzOHgwvjKkhaufXklNCV/t+IwlB5dwvOo4YV5h3Df0Pq7pfQ3BHs1/bT39ezJvwjzWZa/jlc2v8NAvD5EUksTDwx9mcMjgdseidG8dHeHaKqVMs3VQ7dW7d285duzYxinF9jA2mMg5UkrG3kKO7SmiMNP8F7CHj4Ho/ubkK7pfIJ6+zr8iSkpJndFEXU01dVWl+AaEYjAYKM7PoSDrMKW+idSZdBhzU2lIS8PDP5mS7GqyDuZTY/WHv96zHl+3bHqPHUdUT3/Sdy+Bw9+zIHg2tQ0mxpR+zbDq9Qx8ZDUueh0b3nmAwMzVTNO9TG29iX/yNmN02wl7+igAW/41heDiHVxY9xoArxv+Q399BglPmWskNn7xIrXClcDzp9Mv3LdxhK2kpoQtx0+0aNhftN/cA0u40D+4f+MKwqTQpDOuuLKVhnojVaV1jXVmgRFedulYbWtSSvbt2ED/JPM00uZlbxDa7zxi+iTZ/LU2f5/GxmWpDhvh6ggpJcaSkpMSsbrMDOqPHaM2LY2G7JyTzncJDz+RhMVpiVh8PIbwsE7fH/J01i45xM417WsHYa2h3sg7f/+NgRdEcv60zlkgV1RTxO6C3ezI38Gu/F3sLthNeX05AJF1CUza8jdyztlMVWwuAoFO6BDC/FmH1W2hQyBwywrCd9UAlg36NzneqZwXcR7TEqcxOmp0mxbIGE1Gvj78Nf/d9l8Kawq5NP5S7k++n68e2g+oVYrKCe1NuDZLKVOEEOOBqZj3TlwspezIZtUd1tbGp61RWVrbWHyfsbeocVQjJMaHmP6BxAwIpEeCH/p2Tp11tvpaI/kZ5RxPK2usvSo9Xt34uG+wO6GxvoTE+tAj1peQGB+bTBk0NDTgovXIKi06TkVFBdXuIdQ2mBCFR6C+kv7JF5z0nONVx9mat7VxBaGl8aab3o3BIYMbC9wHBw8+tc5Ec//99wPw2muvdfhr6I6qK8qofmUAq41DOTDqJe4b3xtfd9v1sOrKCVdLTNXV1B07Rt3Ro9QePUpdWhp1R9OoO3r0pKlK4eqKa2xsYwLmGh+HW1wcrvHxdtkv8kw+e3ojXn6uTL5/aIevtey1bVSX13H9EyNtEFnbmaSJtNI0dhbsZM+hQ/gvG8aupB/ICz+ESZowSRMSeeK2lJg4cTu0IIELdt5IzaT9TLngcmJ9YzsUT2V9Je/tfo8P93yIlJLb174MqClF5YT2vpOuFELESSlXAauEEPHAFMChCZc9ePm50XdUOH1HhWMySfLTy8nYW0j6niK2/pTOlh+OYXDXE5UYQMwAc/G9b7BzTC001BspzKxsLGg/fqyc4pxKyw5EePm7ERrrQ99R4Y11V+7e9mkY6WLVkNQvMBS/wNATD0Ykm/sGlWee1GQ0vdzco8xSvHtZ/GWkhJl7YLW2SHX79u22/DK6HQ9vXyrvXs/en4/w/u9H2bBlK0/0P86wq/6KvpltdJQTdB4euPfti3vfvicdl1JiLCg4JQmrPXCA8lWrTmrTog8MNDfwHTKEwJtvwhBhvxYeZYXVFOdU0v+8ju99CeZpxQ1fp1JVVueQEX+d0JHgn0CCfwKjPSv5fNlGHkx5sG01XDt3cFP/mwjz7Xji62Xw4r6h9zG1z1T+s/U/sNZ8/KqvryK5R3LjH4iOWASgOId2/Y8qpXxEG91K0+4fBV62YVxOSacT9IjzpUecLymXxVNbVU/mgeLG1Y9Hd5g3bffv4WluPdE/kMjEAAyu9p9OMBpNFGVXmkeu0s0jWEVZlZhM5uzKw8dAaKwvCUNDCI01F7Z7+dlvpVpbZFVk8dS6p9iQswEAPzc/kkOTmZY4jZQeKSQGJrZ7WxavhgYe2b8fSkuhk0cTuorgHlE8eX0U11xQyv7PZjNw15ekD7uS+Phejg6tSxJC4BISgktICF4jRpz0mKyvpy4jk7q0oydGxo6mUfTRRxR99BF+l19O0PQ7cOtt+xmi9D3mLbBiBnSsfssiMtHcEiHrYHGrk5yzQZhXGC9c8AJvfGoefwj3Duf7o9+z+KC59jTSO5JhPYaRHGpOwpo2i1W6rza9iwkhfKWUZQDa6NZZzc3TQM+hofQcGoqUkuJcS/F9IXt+z2bnz5noXXSE9/Izt54YEEhghFeHf7lMJklJblXjqNXxY2UUZFZgrNe2vPFwITTWh6SJMY3tGLwD3Jzul1pKyeKDi3l186sA3J98PxdEXUAv/1422/rj3MJCLigshOXL4eabbXLN7mpgpB8DHpzH4b3T6a0lWxs+f4G486YSFqNKRGxBGAy4JcTjlhB/0vH67GyKPvyQ4kWLKV22DO9x4wi6czqeQzs+9WdxbHchPoHuBITZZi/O0BgfXN31ZO5XCdeZvDnhTYwmIweLD7Ilbwtbj2/l96zf+d+R/wHmJsDJPZJJDEikp39PEvwTiPaJxqBzvu2plI5pVcIlhPADVgHJQggJzJFS/sPqsWnALEcX7jXtw9WZhBAEhnsRGO7FkPHRNNQZyT5cQvpec/3Xui8Ps+5L8zSeufg+kOh+gbh7nfmXSkpJaX61tv2Nuc9Vfno59bXmaQkXNz2hMT4MvDCycVrQL6T9/cQ6S3ZFNk+ue5KNORsZFT6Kp8992mbbfFi7LEcrbH7vPZVwtYLQ6eg90FwOcjzrKEP2/4sFe1Mxjp7FXRf2xKO9o7XdZ0MLuzBERNDj0UcJuusuij/9jOKPP+bY6tV4pqQQdOd0vEaP7tDvtLHeROaBYhJHhtns/wadXkdEnwAy9xfZ5HrdmV6np19QP/oF9ePm/jcjpSStLM2cgOVtZevxraw4tqLxfBedC7E+sST4J5iTML+Ek7rbLzm4hP9t+R97y/aC6sPVZbR2hGs+8CLm1g89gWlCiB+llBdLKUuBBUKIt+0VZGv5+fkxf/58R4cBgIurXmuoGgRToLyopnH06+j2fPavy0EIcwNBc/F9EKGxPlSW1plHrtLMI1f56eXUVpk3r9W76AiO9qbvOSdqrvzDPDu8V2Rnajqq9cSoJ5jaZ6rtEsQJE2DVicHXgZbmsWvXgvVrjB8PK7tMJxOHCI2MJ/eOdRz9NZ+vVh1izx+ruCvJhWGX3tEpOyScjVwCAgj5y70E3X4bJUuWUPj+B2TMvAu3xESCpk/H99JLTtmgvTUs7SBiBwS2fHIbRPUNIG1nAWUF1U5Tu9oVCCGI94sn3i+eKX2mAFBVX8XRsqOklqSSWprKkZIjHCw+yKr0VZikefbiLv4NwNPrnybQPZDz+5zPFraUnvaFFKfS2t/czVLKpdrtrcBWIYS/EOJNYLY2zaj+hj0Dn0B3+p8fQf/zIzAZTeSllZOuFd9v+i6NTd+modMLTEbzt1GnEwRGetFzWCihMeZpwcAIry6zGrI52RXZ/HPdP9mQs4GR4SN55txnbD+q9dhjsH49aP2TDCbzf1TU1Z04x9MTHn/ctq/bTYXF9OZfN/fmhqNFlH7+HtGbtlJy3lQCTrNxeFNOPtDqtHSengTeeisB119P6bffUfjuO2Q//DD5//43gbffhv8116Bzb33T2vTdhehcRGPdla1E9TVfL/NAMf1VwtUhngZPBgQNYEDQgJOO1xprOVZ2jNSSVI6sNx9bdtUy4n3jEULwby0JU5xfu5chSSlLgLuFEA8JIZa2dL5ygk6vI7ynH+E9/Rg5KYGainoy9hWRl1aGX4gHobG+BEV54WJwjt49HSWlZMmhJby6+VWklLYf1bI2dix88w1ccUVj0nUST0/49lsYM8b2r92NjYgPxDhrIceO7KOHvz8mo5E/PvoHvS79C8Fh0Y4Or9sSrq74X30VfpOvpGLNGgrfnk/eM89S8Pob5oTsxhta1f3+2J4iInr523z3jMBwLzx9XcncX0z/8zpnk3Sbc/KhAje9G30C+tAnoA9vaI0AWrN5tuJ8WjtcslIIcacQ4pAQIsn6ASnlK0Ay2HSbvnZpMDU4OoR2cfc20Ht4D86f2ptBY6LoEe/bbZKtnIocZq6YyTPrn2Fg0EC+nPwl0xKn2bfGbOxYWLgQmo4AuLubj6tkq130Li4kJA4CIHXnWoalLeDVN99i/q9HqGswOTi67k3odPiMG0fsF58T+/FHuA8cQP5rr3F47DjyXn6Z+rzjp31ueVENxTmVxHawu3yzcQnzqFnmgWLa09PRoRz+jqWcbVr1546UcpsQIhVIlVJub+bxpUKInqc+s3MdLj7M1FlTufWCWzu9aF45mZSSpYeW8srmVzBJE4+PfJypiVNttvqwRSUl4OKCEajX6XAHcHExH1c6rNfQ0Rzz+Y3jv1Xw+Xf7ObruS64fEc3gsdc5/YKNrkwIgefw4cQMH07N/v0UvvMuRe9/QPFHH+N31WQCb78dt/iTV0Ae210I2K4dRFNRfQM4tCmPouxKgiK97fIayqnU5tVdT6vHl7Xi+NO2gtB6cTmUQW9g/4D9HAw/SL2xHoNeLat1hJyKHJ5a/xTrstcxMmwkT5/3NJHekZ0bxLvvQlUVR7y9eTs+nlddXGDHDrVa0YZiew3g3V7w84Hj+C56CbffK6m/YAquBtUwtTO49+1L5CsvE/K3v1L0/vuULP2SkiVL8brgfFwCAkGnAyE4WDEYT50PNfPmkqMX2oIHYX5cJ8wJstBp55tH0xA6PEeOwPu881qMo7GOa3+xSrg6kWU1/oIFC1TRfBfR4v+MWtuHR4ECYIkz7Z3YVIJfAlcnXs0Hez5gc+5m5l44l2gfVV/SWaSUfHnoS17e/DImaeKxkY8xLXFa541qWfPzg5dfZuayZUghzCsXX3sNfvut82Pp5sYmhlI3awW52em4Glyorixnx+f/pN81/3B0aGcF1+howp58kuB776Xo448p/+FH6o6kgpQY0XG85/mEl+6iYtdqMEkwmUBK8xSgydTsfWkyUTh/Pn6TJ9Pj0UfQn2GRhG+QB74hHmQeKGbIePX/raKcTosJlzay9Yil35YQIgE4AiyyNEF1FkIIHh/1OKPCR/HkuieZtnwa/zznn1wSf4mjQ+v2citzeWrdU6zNXsvwsOE8c+4zRPlEOS6gr78GQP7P3FwQvR4efND8odicq5s7MfF9ADi4/huGZ7zHjNdjOTf8QgCks1cmdwMuQUGE3n8/odr+oQAZ+4swvradwY/dRvyQWa2+lqmujsK33qJg/gIq1q4l7J9P4jtx4mnPj0oM4PDmPExGE7ouvJJaUeyp1b8ZUspSKeUCKeWjmKcWZwohXhRC2G5nThuZEDuBxZMW09O/Jw//+jBPrXuK6obqlp+otJllVOvqZVez9fhWHhv5GO9c9I5jky0rSUlJJCUlOTqMs8qQCTdw+IbfKe8xgh/35AFQVpTv4KjOTu1tB6FzdSXkr38lfvEiXEJDyLrvr2Q+8AANhYXNnh/VN4C6GiPH08ttEbbShQkhJmgfc6yOFQshtjQ5dsp5Ta7jL4SYJYSYYnXuLO2xZMv1tPMmCCGOCCFm2P8rbL92/SkipTwqpXxZS76KhRAvaclXkm3Da79I70jev+R9pg+azpeHvuSGb27gUPEhR4fVreRW5nL3qrv557p/0i+oH0uvXMr1fa93zBTiabz22mu89tprjg7jrJPYdyALZ4zimgHmet7dn6npRUfoaDsI9379iF+4kJD776di5SpSL7+C0uXfnLIiMSrxRB2XI6gRVOcghEgGJkopV2LemcbSv2KqlHKYlHJ2C+dZW4y5jGmJdt5mzI3XkVJuBVKBhVLKEu3xEmCRPb++jurwO6OUcpuU8hEt+eqpJV8P2SC2DjPoDPwt+W+8NfEtSmpLuOHbG1h8cHHXW77sZE4a1crbyqMjHuWdi95R9XLKSYQQ9I4OA6Dvjc3+EavYkaUdREdXJwqDgeC7ZhL/1ZcYYmPIfvhhMu+5l/q8vMZzPHxcCYr0dljCpTgHKeVWKeVsIYQ/5q4GqdpD/tZJ1RnOA8yjX9p5qVbPKQFa3NHGalRsihBiRUvndyabDkVIKZdKKR8BFtjyuq1l2UtRWyrb6NyIc1ly5RKG9RjGM+uf4aFfHqKszqnKz7oM61GtvoF9WTppKTf2u9GpRrWs3XzzzdysViU6XGBgIPV1NexcvdDRoTgdY4OJmsp6m1/X0g4i1kbtINx69SLus88InT2byvXrSb1iEiVLlzb+ARvVN4DcI6U01Blt8nrKmS1fvpwZM2aAc7aFSME84mQRCBQ1swVg0/MsEpo7ro1sWbNMNU4A/LVjW6WUS4DrgJltDdyezvguKYTwFUJMb+tFtUL7TmfZS7G5HlzBHsG8OeFNHhj2AKvTVzNt+TR25u90QJRdk5SSrw59xTXLrmkc1Xr34neJ9nXuUa3MzEwyMzMdHYYCbF34PIN/ncHPv61xdChOw9hgYvl/t/P+7N9Z89kBygpsV2uavqcQ70A3AsI9bXZNodcTdNufSVj2Ne59+5Lz2ONk3DGd+qwsovoGYGwwkZPquC4FZ1MPuEmTJln2Dna6thDaFJ+/EGKKdn++NkJVYjnW3HlWNmNOugAQQiRo9VpbtFExi5VSypVWU4pIKVdqCdgpI2eO1tKwxD+AAK0+6xohxGEhxI9CiJb3knBCOqHj9oG388GlHwDwp+//xHu732vcGFRpXm5lLvesuocn1z1Jn8A+Tj+qpTinpKmP8mLA08z8sZrNaUWODsfhpJSs+ewAWQdKiOkfxL512Xz65AZWfbSPkrxmtqVqA2ODicz9xcQODLZLEuIaG0vMhx8Q9s8nqd6+ndRJV+KxbQU6nVDTimcxLSmyFK6XAIFCiBlWCVXh6c6zvo6lRssyDaklTisw7+tc0opQZmtTlk0TOYdq6R1zk5TyZcyFaIFSyl6Ye3I9avfI7GhIyBAWTVrE2Jix/GvLv7hn5T0UVBc4OiynYz2qtTl3M4+MeIT3Ln7P6Ue1FCckwc3dk7um30OkvwcvfrSMrNQ9jo7Kobb9lM7+dTmkXBbH5fcM5pZnz2HgmEgObcrjs6c28NO7eyjMrmjXtXMOl1BfayR2QGDLJ7eT0OkIuOEGEpb/D4+hQyl64Rn86nPJ2JFrt9dUnN7bmBOlCYC/lHI+5vyhxKoua8lpzjuJlHIqYKnFmoJ5C8Ej0Fh0nwBcZ1mlqN2fpq16XKEdu87eX3BbtLR0RULj1j6B2u2tQoj4Mz/N+fm6+vLqha+y+OBi5m6ay9TlU3nh/Bc4J+IcR4fmFPIq83h6/dP8lvUbyaHJPHves8T4xjg6LKWLC/By5f1bh6KbN5Piz4KJeOz3s2oayCJ1Wz7rvz5Cr5RQRkwy/3fqHeDOBdP6MOySOLavTGf3L1kc2pRHwtAQUi6NIyTGp9XXP7anqF3tINrDEBlJ9DsLKP3yK45+uIWjuhBy5n9A2B23IPTdY09YpXW0kSjLNN5K7ViJ5bbVsVPOO8315p7m+FZgmNWhlUBzP+ynvbYjtJRw/UMb0lvJyXuqd4sxYyEE0xKnkRSaxMO/PMzMFTOZPmg69yTdg4vu7NyeRErJsiPLmPvHXOpN9Twy4hFu6HtDl50+POcclUA7UnPJVFyoH/uumIebT9hZmWwdP1bGivf20CPOl/G39jvle+Dp68q51/Qi+aJYdvycwc7VmaRuyyd2UBApl8YRltByjXT6nsIOtYNoKyEE/tdew4CooRx9/yiHPvye+o2/E/XmPHSurp0Sg6I4u5beRV8EtgHXY+42/6MQ4kXg9C2Hu6A+AX34/PLPubr31SzYtYDbfriN7IpsR4fV6fIq8/jL6r/wxNon6B3Qm6VXLuWmfjd12WQL4MUXX+TFF190dBhKE/1SxtM7cQAAu9YsRZrOjjrKiuJavpu3E3cfA5fdPRgX19OPALl7Gxg5KYFbXziXkZMTyEstY+ncLSx7bRtZB4pP296mvKiGouyOt4Noj6jkWFwMOmovuYXKtWvJfngW0qhWLSoKtJBwaW0eVml9ti6SUl6MebSrZ+eE13k8DZ48fe7TzB09l0Mlh5iyfAorjznVaKTdSClZdngZV//vav7I+YNZw2fx/iXvqylExe52rlnCoDW388PCeY4Oxe7qa418O28HdTVGrrh3CJ6+rRv5cfNwIeXSOG594VzOm9KLouxKvv7XNr56dSvpewpPSbxs3Q6iLfQGHeG9/CiQIYTOmkX5jz+S++yzqvehotCOPlxSylXAbDvE4hQujb+UxVcsJsYnhgfWPMBzG56j1ljr6LDs5njVcf6y+i88vvZxevn3YsmVS7il/y1delTL2rXXXsu1117r6DCU0xg0+ho+jnyCe3fE8vW2LEeHYzfSJFnx3h4KMyu4aPoAgiK923wNg5uepAkx3PLcOYy+vg/lhTUs/+8Olry0mdTt+UiTOamxRzuItojqG0hRdiXu195I0J3TKfliIQX//a9DYlEUZ9KuCX4p5VFbB2ILlsankyZNarYXV2tF+0bz8aUf8++t/+bDvR+y7fg2Xr7wZRL8mtt9oGuSUrI8dTkv/fES9cZ6Zg2fxY19b0Sv615FroWn2ftNcQ5Cp2PabQ/wzbt/8NySdfQyhjIw5UJHh2Vz678+wtEdBZw/rTdxg4I7dC0XVz2DxkTR//wIDmzMZcv3aXz/1i6CIr1JviSGzP3F9BnRw2H1cVF9zbXLWQeK6f33v9NQVETBvDfRBwQSeItqQmwry5cvtzT5dsbGp0oz7NL41FHO1Pi0rQx6Aw8Nf4g3xr9BflU+139zPV8f/rpbDI0frzrOfavv47HfHztpVKu7JVtK1+DmouftW4bxmvsCenzzJ8rLu9cuEHvXZrPtp3QGjo5k8Fjbbequd9HR/7wIbnp6FBNu64/JaGLFu3uprzU6pH7LIjjaBzdPFzL3FyOEIPzpp/GeMJ6855+n9JtvHRZXU119uYYzNz5VmndWNT5tj9FRo1ly5RIGBQ/iibVP8Ojvj1JZX+nosNpFSsnyI8u5atlVbMjZwMMpD/P+xe8T6xvr6NCUbq6lP1P8PV2Jv+lfpF74X3x8us9/L1kHivnl0wNE9w/kgut622XUSafXkTgyjOufHMnFdw5k8NgoYuzYf6vFeHSCgDAvyotqABAuLkS++iqew4eT/cgjVPz2m8NiUxRHOisbn7ZVqGco8yfO596ke/n+6PdMWz6NPYVdq2ljflU+f139V/7x+z/o6deTJZOWcOuAW9WoluI0IhMGMHKseXT6wOZV1NZ0rNu6o5XkVfH927vwC/Xg4ukD0OntWxep0wl6DQvlguv64GJw7O9107xS5+ZG1Lw3cOvdm8y//o3q7dsdEpfiGEKIZCHECiHErCbHJwghjmiNTFu6xpymz29ynWY3qtYap0443XM7U0v/AzQ2PgWOare3Yt7n6Kyi1+m5a8hdvHfxe9Qaa7n5u5v5eO/HTj/FaD2qtT5nPQ+lPMQHl3xAnF+co0PrFOPHj2f8+PGODkNpg+y0A8Qvn8qvCx52+t+v06mprOebN3YgdILL7x2Cm6fB0SE5nN7Hh5gF83EJCSFj5l3UHj7s6JCUTqLlDXNovvN7ajObUjdn4Rmu37iXojWrvRxXYtXt3lE63PhUCHGNdjNVSrnd6niclDLNRnE6jWE9hrFk0hKeWPcEczfNZWPORp4971kC3O3f0bmt8qvyeWbDM6zJWMOQkCE8e96zxPt1+U0C2uSJJ55wdAhKG0XEJfJD32d4eEcP7lh1iPsn9HF0SG1ibDDxw9u7KC+q4ar7h+IX4uHokJyGS3AwMe++Q9qNN5J+x3TiPv8MQ0RExy/cNfNyh4l75Ns1zRxelPbS5fPiHvnWE/iumcc/SHvp8g/iHvk2GFhi/UDaS5ePaeVLrxRCJGs71vhzots8ANr+ipuBFMt2P9rI1FbMW/tYzvMHZmjHE5rbGkgznBOJWqp2DYf1e7JF41MBBGHOHsdZ1Xf52zpYZ+Hv7s9/xv6HR0Y8wrrsdUz53xQ25W5ydFiNpJR8k/qNeVQr2zyq9eElH551yZbSdV18/b1clNyH/67cz5pVzlNo3RIpJb98doCsgyWMu6Uf4b38HR2S03GNjibmnXcwVVWRfsd0GorURuZnkbeBmdrtBKDxH19LrDZbbVw9Q0vAtmojVNaJ0qPASu249RY/Tfk3ue+41SS0MMIlpVyq3VxlOaYNyc2wPkcI4Yf5G5AAxAshFmLOJrfbOmBnIYTgpn43kRyazMO/Psz0n6Yzc/BMZg6e6dC6qILqAp5Z/ww/Z/zM4JDBPHfec2d1onXppZcC8P333zs4EqUthBC8eM0gzk1/i3N/XcjRmF+J7z3Q0WG1aNuKdPZpG1InjgxzdDhOyz0xkeg355F+x3QyZswk5oMP0Ht7OTqss8aZRqTSXrq8CjjT4wVnevxMpJSp2qxZcyZyYuQslROJWXN9P5Mxb1CdjDmJO50SwHErSJpoVRWn9ciVllGe9A2QUpZq3einSSmHA9uklF/aPlzn0y+oHwuvWMjl8Zfz5o43ueOnO8itzO30OCyjWpO/nszarLU8OOxBPrrko7M62QKorq6murra0WEo7eDqomPC7c+wftAzxPUa4OhwWpS6PZ/1Xx2h17BQRlxxdv/etYZnSgqR//oXNfv2kXnfXzDV1Tk6JKVzLBZCzGmmbmsr5kEbtM+btA/LVKJ14rQCGmvDTpqWbGITJ0a5EizPc5TWLpspBtK0KcVrWmp8qnWjP2t4Gbx44YIXeP7859lbuJepy6fyS8Yvnfb6BdUF3P/z/Tz626PE+cax+MrF/Hngn9UKRKXL8wsM4cIpf0EIQVbqHorzcxwdUrPy08tZ8d4eQmN9Gf+nfghdV+/y1Dl8xo0l/LnnqFq/gexZs9W+i92UNjM2W/u8yOqhaUCCVtc1G0jWzkmWUs6VUs61PI458ZoohPDXjlvOTbE83nT0TEq5RHv+BO2+Q/fra22n+UBgvLZasZEQYigQIKVcbfPIuqAre17JoOBBzPp1Fn9Z/Rdu7nczDwx7AFd96/ZMayspJd8f/Z4X/niB6vpq/j7s79zaX7V6UJxQB4uaa6orcfloEofdejP44e9wc3Gen/GK4lq+fWMH7l4GLrt70Bk3pFZO5X/1VRiLizk+dy65/n6E/fOfDuuSr9hHMzVYs7Xj84H5VufN1W6ubObYVmBuM8ctmt3juZnzHKa1I1x+TZMtaGwXscVqpeJZL94vnk8u+4Qb+97IJ/s+4ebvbuZY2TGbv45lVGv2b7OJ9Yll8aTF3DbwNpVsKc7FRu+b7h5epI16lofLpjJryU6naRdRX2vkuzd3Uldj5PJ7h+Dl5+bokLqkoNtvI2j6Hdq+i687OhxFsYvWjnCddhNBKWWpEGKVEGKcvUe6tKWglj4aK6WUJfZ8vfZy07vx6MhHGRk+kifWPsG05dN48pwnuTzh8g5fu+mo1gPDHuDW/rfiomvXtpjd3hVXXOHoEBQbGXnJTUx1O8zLPx5guFs6N1892aHxWDakLsgo57J7BhMc1fYNqZUTQh58kIbiYgrmzcNUU0PIvfeg81KF9Er30dp36TP+naolXf4dD6dFj0opZ1v14HCaocLmjIsZR7/Afsz+bTaP/PYIG3I28OiIR/E0eLbregXVBTy/4XlWpq9kUPAgnj3vWXr6NzuKqmgeeughR4eg2NA9Y3rik/otN+94ki2+Cxg2fprDYtmwTNuQemrHNqQuzs+hvraa0KjT/l17VrDsuyh0Ooree4+yb74h9OGH8L3iCjXFqHQLrZ1STG3FtGGHll4215q/mZb8/lYPd4lMI9w7nPcufo8Zg2ew7PAyrvvmOg4UHWjTNSyjWlcvu5pfMn/h/uT7+ejSj1SypZx1hBBcf/NdrEl4iP7nOW6Ea9+6bLb+mM6A0ZEMHte2DakL8zLZ+NNCnvh6Nxf96xd+/8+fcVkwmqoa896D0mSyR8hdgnBxIfzZZ4n9/DNcQkLIfngWx268ierdXWsrNUVpTqsSLq0f1/VCiLFnOM2/I4E0XT1wmpb8JVanHOnI63UmF50L9w29jwUXLaCivoIbv72RL/Z/0ao6lMLqQv6+5u/M+nUWUd5RLJ60mDsG3aGmEFtpzJgxjBkzxtFhKDbk6ubGmFufwMPdjfLSIrJS93Xq62cdKGbNJweI7hfQqg2p83LS2bz8bapqagE4uPxfpKydyU9bDxLm50HV8L+Qev6reLq7A7DvpdH8+N/7+O1QPkaTc9SqtUdHyuw8hw4lbvEiwp9/jrr0dNKmTiXniSdoKCy0XYCK0sna8q59J7BKCHEEeEFKucPygBBius0ja74l/9uWRAyrlQ1dxcjwkSyZtITH1j7G8xufZ0POBp4+92n83PxOOVdKyY9pP/L8xueprK/kb8l/488D/qwSLUWxkjbvWvzqcim8fxNBfvavoSrJq+L7+dqG1HcORN/MhtQlBXkc+n0x31YPYnWGiaSSFfzH9Q22hfRh6KixxE6YwZGiSawddB4uhpP3WKyrraXcpycbCt15/90/iPI18FLIj8SPn05kQj+7f302Y4MZQKHT4X/ttfhcdBEFb8yj6JNPKPvhR4LvvYfAm25CGNT+lErX0urt67XmpinAFmCbEMIohDgkhDACw6SUr9g4Nv8m94OklKlSyiXaR0nTJ+Tn55OSktL4MX++8+VkQR5BzBs/jweHPcgvGb8wdflUth/fftI5hdWFPPjLgzz868NEekey6IpFTB80XSVbSpdlr1WFholP8M+GPzPjs53U1Nu/h9MP83chOLEhtZSSo1l5/LH0X+zf+QcA+ZkHGb79MaoOrqZvmA/DJ17H4au/Y3DKBYB5r8g+yWNOSbbAPHo38r4Pmf3YS7x+41AuCsxlVOZ7pO7ZCEBFWTGlhZ3fWNmR9D4+9HhkNgn/W4ZHUhLHX5pD6uSrqPh9rU2u7yQLXttk/vz5je9zQPsLCDuJECJZCLHCqjzIcnyCEOKI1kerpWvMafr8JtdptqnpmR7rbG1+B9d6WszVvkEBmPc+KrV5ZO1oyR8SEsLmzZvtEIpt6YSOPw/8M8N6DOPhXx/mzz/8mXuT7uX2gbezIn0FL2x4gYr6CjWqpXR5wlZ9IU6jb8o4prj2497PtvJ/ny7j0VsnI+zYGqUwu5KhF0Wz7+fX2VTqz8fH46kuL2a729P8VjyTvoNHED9gFEcNP/FS3xR0+vbF4m7Qc8XgCK4YfCv5uWMY5md+T9397Vsk73+Zf8R9yriRyVyYGIKhmVG27sgtIYHo+W9TsWYNeS++RMb06XiPG0ePR2aDvh0lxF24Dn/GjBnMmGHeYU8IUeDgcFqkbVY9B5jDqYvdUpvpOt+chZzoUtD0+iuFEDPP8Fhz2wN1una/k7fyG9QRbW7JX1payowZM5g0aRKTJk2yZ2w2MShkEIsnLebp9U/zn23/YcnBJWRXZtM/qD/PnfccvQN6OzpERXF6lw8OpyzDwDV/3Mq6T3dy3i1PAebi81qjpLrOSHVRJrK+tnFabu+67yiurCbdbziVtQ3EH1uMh7sH5075KwAb5s0gq9aDLzyuo7LWyNOlj1HjGw/8Cb1eR/yOeeTohjGq92OMTOhNRtA6LkzoC4CLwUD8gJE2+/pCwmIab4cNvZhfK8r5MdOFzw5sZrbn/xgS6sKoGa+jOwu62wsh8Bk7Fq/zzqPoww8pfPMtUi+/Av0NdwF9HR1ep1q+fDnLly8HOLUmpSVP+a0BPuCp0g94ys+A+f31HZ4q/YSn/DyB74A3eap0IU/5+QHLgP/wVOmXPOUXjHnPw1d5qnQ5T/mF8VRpa4ddV2pd5bdq3QZO2pZH26x6M5CiNUW1bGq9lRNb/FhaRM3QjidYznV2bUq4tP0Ur8O8yaT1P3IJ5n+wRVLKsvYEotVmpQghplimDYUQs9rSkt/Pz88ppxHPxMfVh5dHv8w54ecwb8c8/jr0r9w28DY1qmUj06Y5rm2A0nmuv3QCvxX+jV4T7wRg6yuTcCnL5Mq65wD4yPAiwYZaIp/cBIBc8yKGmloerfsnAItdl6B38wHMCZdLdT6uDd4YvHVE+LtSIgbgHhQFBwAhcLnnd64KDudqXeeOLsX1TSaubzIbjCbWHMjH9YfPoayiMdna9tMnRA8eTbBVktYd6VxdCb7zTvwmTyb/1f8ja+kSGP44denHYFioo8PrFJaBhQULFthjhsle3sbcZX4m5oGUxoRLS6xWaslYoJZ8AWzVRqmKODHC9SiwUDv3TJtXO5VWv6sLIe7E/A1aKKVc0MzjQ4F/CCEK2lPPpe15tKTJMafus2UrQgiu7XMt1/a51tGhdDv33HOPo0NQOoHQ6Rh965ON9+t7XUJxSSGz4hLxNOgxlj9ArduJESD/69/CR+rZEByHp5seT5eLcXE58d9hyoNfAXBinDwFgJ2/m3s7B4ZG2vXraYlBr2Ni/x7Q/xNMRnMbidKifAau/Suf/zaBfUOf4N6xPYkKaF/Pv67CEBpKxJyXYPVmWFRG4fwF1A3rgWtM9044O+yp0jFWt+sB6/tVTe6XNrlf0OR+q4sKpZSpTfc7tDKREzlAKuakDLRtgJpIBlZopU3dK+ESQjwMvH2m0Sttm59tQgg/IcRDdiiib1FXm1JU7K+qqgoAT8/u/cajnGzkVfcCML7xSPxJj0cmDOjUeOxJp9Vw+QWGkHbDSrJ2VbB0SyY7tqzjqR6/Ez/tBULCoh0cpX25JfYF/gCTiYw7ZxD7xee4BAQ4Oiy76tCUomMtFkLMsWpibrGVE6NeCZjLisCcXKVyck33CmisDbO+hlNrbR+ul1s7VaitZuz0ZAtOTCmqZEuxuOyyy7jsssscHYaidIq4vsk8OnU0ax4ew/SEEnoV/ky9yTyyZzLafxWnowXdOZ36nBwy77kXk9ZItruaNGmSpYTG6acUtdKg2drnRVYPTQMStLqu2UCydk6ylHKuNsuVoI1kJQMThRD+2nHLuSmWx5sbPbMuV7Lzl9kiVSikKIrSzUT4e3DNHbMpL7+LQB/zAMiOVy+n1LsXvW98mUh/DwdHaB9uCQmEzJ1L1gMPkD37ESL/9X+ITq6zU06l1WBb12HP1o7Px6qnplUZ0cpmjm3FaoVjMyVHzW690ly5kqPY7CdRCDHOVtdqL8uUojbMqiiKYltdrGeTj5ZsNdTXUefRg7U5MObln3nsy53kZKU5Njg78b3kYkJnzaL8xx85/rJDJls6xfLlyy2tIbralOJZy5YjXBOB1Ta8Xpt1xVWKitKtdf9OBV2Ci8GVkfd9SFRJNVU/HyZ7y7cE7ZjL1vEfkTz6CkeHZ3OBf/4T9VlZFL3/PoaICAJvudnRIdlcF12leFZrbdH8ZmDomU7B/Lffo7YISlEUxVl15Rwy0t+D568eRG6KB5uWpzMoeQwAh7b/hm9wOD2ierXrulJKiqvqSSuspOzQOspyq/EIdtxKTiEEPR59hPqcHPJeeAFDRDg+48e3/ERFsaPWjnDNBoq0lYjNEkK8ZJuQFMV2/vznPzs6BEVxOmHRvQi7Zx5gTpbE8r9R1NDAG0M/5e6xvQj3O7XGS0pJUWUdWccOUXvoV/qMvxU/H282Lv03ibvmcl7Nf6jGnbv1/yOuNgmDwbwydNPX82ioLiX+0vsJ68TaMaHXE/nKyxz705/JevAhYj/8AI8hQzrt9RWlqVYlXFLKVa2o0XL4XkWqLYTSlEq4FOXMhBB43fo5X/22jc/+yGDJpjTejPyR/lMfJyQkjD3rv8d15RP8teE+9tUGM1n3O/92nceuuBQGJY3AO7w3B7In8o8B8USERxLnm8K2T9Mbm7HqD/+ALMtl1I7BpMQGcGfEUYYOO4fQqNO1Y7IdnYcH0W/OI+36G8i4+x7ivvi82/To6sJtIc5ara7hklKesT5LSrmq4+F0jKrhUpoqKDBvMxYc7PT7uyqKw4THJnJPbCKTiqr49tsvOffwZ2zadiEhF03FzdOHaoMvl/cOYGpUf3r59ibD/Tr6xvcHYMC5l8G5l2G9mdEOXWbj7eSH/kdqVi4P7i/nx53pnLP1IXYeGk3oQwsBKMzLJKhHlN2+NpegIKLnv82x62/oVj26VA1X16PaQijd2pQp5tYra9ascWwgitIFRAd6ctctN1OQcwF9XdwA6DXkfBiymkEduG5CZBj3RYZx3/jeZBz8nnCtN1hO+iF6vDuc1/3ux3PEn7h0UFiz05kd5RYfT9Sb80j/821k3nMvMe+/Z/PXUE5P65M1B1hh3c5B66P1NjC1pf2Ztc2vC5vbgcbS50tKObHJcX9ObAc0XOv15TCqQYmiKPbXxdopnO2Cw2MJDAmzy7Wj+yTRs6+5lsrN3ZON8XexhQE8881eHpjzOgeeG0lObjYAu9d+wx//uYXnlm7kn8t288GnH7Fx3p0Ul5gHdfas/YaNb/+Fj34/DMDGLZvZtPT/qK6pBaAgN53Mw7spqqjFMCSJiLlzqd6+nezZjyBN6oeys2jJ1BzMezE3ldpSsqVZeIbrr8S8p3NT04BArReXZXNsh+lWI1yqhktRFKXrCAyN5Jw/v8Q5QGp+Bbt/KUTsM1KjdYmvyjtCQtGvrCq5mmK8mWraxjXyWyrqXiAAKD3yB8nZX/DgkYu4AW8KDq7nOs//UnHpvXgAB5e9TEr2p/Sp/QgQzHZbysQkSfmPP1Ln3xfoQ0lBDuE9u14ZVHtruAZ9OOg1IMnG4Wzf9add97fivJVaV3nLljyp1g9qCdFmIEVrimrZ1Hor5k7zlvP8gRna8QTLuU01OZ6Ag/dd7FYJl6rhUhTnIrpyDwWlUyWEeJMw5Vbg1sZjI675G1zzN35uPHIR8C98tXvn3voM8AzLsir44tk/GH/1dPJjbifU3RWAHufezLajg/ln0ADKaxoIzUoiN0jQa2BPDn/zIwzpw/FvniVi8Af4eRg67Wu1hS5aw/U25q4HMzmxbyLQmFit1JKxQKvRqK1SypVCiCJOTA8+CizUzm0xidK2/CmSUqa2dK49dSjhEkJMl1K+Y6tgFEVRFKW93N09CIsKbbzfc9BIeg4aaVXQ3xsAaTSSlfU8mGB1ej/e+GgzH90+AneDvtNj7mytHImyCyllanP7HWomcmILnlTMSRlo2wA1kQys0GrDWjNqNUVKObPl0+yrozVcze5dpCjO4u677+buu+92dBhKd6KG7bo8odcTPONOAC7evwGvX1ew4N23MDY0ODiys8JiIcScZuq2tmIe9UL7vEn7sEwlBlqduwIaa8POOGolhJhiKbTXiusdpqNTiup/HsWpXXddczWaiqKc7XSu5lWYHgP68/D6L4h0KWK9XwXn39A5C9nyc9JJ2/wdW1xHkFFlYFSMFxP7heLu6d0pr9+ZLKsIhfmPlUWc2JVmGpCg1XXNFkLM0s5JtkqSZmkjWcnARCHEfCnlXKtz0aYbk4UQCdbThtrrzhFCWF7PoasUO5pwOdUyD1U0rzSVkZEBQHR0tIMjURTFXmQH3olCH3wQ46sFZK1fT+Ik20/afL7qDw5UeqHL2c7NeXMov+hfDBk1nsx9fzB8y2zm1D7JfteBGDf/xIUun5F76xrieia2eN2u1PhUW0W40urQbO34fGC+1XmWlg8rmzm2FZjbzHGLU/7xtNd1mpm4NiVcQoiHOTGsJ4DxWoZpGekqkFI6bHt2VTSvNHXLLbcAqg+XoijN07m5ET5vHpn33EPBE0/SUFdDbohg6MSbbHL99auWscrlAob7GSh3j8DgYq4T650ynmMRP/NJbF9c3TzYvcnAri3VjIoz15mt+/Q5ZH01sZMfIyrA85TrdtGi+bNamxIuKeXL1veFEC9JKdWG1YqinJFTDYW3k+zIMIrSOdpZ5KJzdydq3jwy77mX4meeJ3p4GT/7D2Ts8LbtvXg8K41Dm37ivKtOtHt66qG/8++AQMyDE9Maj3v7BuDte6Lj/eCR42DkiR30DLlbKS0p4vw5P3NuzyBmxuZC4/pMpSvqaNG8+h9IURRF6fLMSdcbeIwayfFNvix67Ss2pha2+vnLtmfx64K/M3jbkxTkH288HhgYhGjHQovhD35J4l+X8cCEPpQU5TNq7fQ2X0NxLqrTvKIoiqJgTrpi334bt3PO5S9bFrH2qX9wdM/GMz6ntOg4j360kr99sZ1lwTMovflHgkNCz/ic1ooO8eNvE3rz7cOXc+TSTxqPp+7fzs6XxrNh05ljU5xLRxMutUpRURRFOVUXnf/QubkR/9Y8DKNGcekfv1L+8n2nnU6ur6+j9L8XMvHwczx0UR/ev/sionq3bRqyNYROR/9RlzTeL885THBNOn9ZesTmr6XYT0dXKb5okygUxU4efPBBR4egKGed9kyhOROdmxu93n6TQ7f9CcOGHZQsXkzAtBP1V7V1tbi5umEwuJI37EFiohIZl9S70+IbMnYKdedfxRN78rhqTqe9rNJBHRrhklI61eoIS1sIbamsoqgWIYrNdfFcQmklnZsbvd//EK8LR5P75D/Z8vdbqKooYe+eneS9OJR1P5r3Uh5+xXR6JV3Q6fH9+MP3fDvvaegCbSG0nlmzhBBThBDF2ucZQgibpItCiAlCiBUtnDPFFq/VEd2qhsvSFkK9wSoWBw4c4MCBA44OQ1GULkjn5kbUf/8LSQPx/G4z3953N9d+mkqm6IGfj49DY5s0aZKlDZJTDXycRqqUcq6UcgnmPQ2XaD24Ntni4lq/rZLTPa5tdj3RFq/VEd1q82pFaWrmTPP2WaoPl6Io7aFzdSXxo0/ZduONDFy/nX+M7UX/V37A38vN0aG1y76+/V4Dkmx82e399u+7/3QPaolWc1Zqm1ZbtvVJxdwUdQ7mzvKpmBOpqdrxFO2x2ZbztWTrJFqHeevHE4AUIcQEbSNsf2CG5XW15M/uWjXCJYQYKoQYL4SIs3M8iqJ0R92hh1U3+BKU9tG5ujL0s8/wuGISKT8voXz2gxgrKh0dVnfwKLBSS4qGaZ8Dtc9LgOu021uAFO12kZRypZYknbJxtbY5tiWJsnS038rJydlJr2vnr7FRa0e4AqSUq4QQ8UKIO4EVUso0y4NCiHhAWh9TFEXp6sXTimKhc3Ul9uU5FA8eSN6cuaRdfx3Rr7+Oa1yco0NrkzONRDlAMrBC2yvRkjxZb0ZtuV1kdazE+vGm+ydqt+dro1hteV27a20NlxBCvAn4SSkXaAcaW+JKKY9yYkdvRVEURel2hBAE3norMe++g7GgkKNTp1Hxyy+ODqsrWwEnRqBa+Rx/q9snJVsAQohkbZqyWVqS1Z7X7bBWJVxSylXAI8D1QohNmIfpkoUQ1q1vg+wQn6IoiqI4Fa9Ro4hbsgRDVBQZd91NwfwFauunFggh/LWVgoHaCsUEbQPqZK3mKkVLhpK1acEJVrcncqLoPVBLqmYAM7VrWz8vAfMoWALmETDL6kTL7dSmr9sZXz+0oWheawHxCIAQYjzmEa0UIcRMzDt7d9qwnKK01uOPP+7oEBRF6YZcoyKJ++xTch5/gvz/+z9q9u4l4vnn0Hl5OTo0pySlLMFcl7WkyfG5TU7tqX2er32AllhpUrWRqa1W19hq9TzrEauZVufMbuF17a5dqxS1Ea9VlvtCiKGqfktxRhMmTHB0CIqiOLGOjEvpPDyIeOVl3Pv35/irr5KWmkrUG6/jGh1ts/iUE7QRqeSmdVtdRWtXKZ5xi3Ip5bam57f0HEXpDNu3b2f79u2ODkNRFLtqR9pko/UcQgiC7rid6Pnzqc/L4+iUqVSsXWubiysn0VYn9uyKyRa0vmh+onWR/Jlo040TpJRl7Q+rfVSneaWp+++/n/vvv9/RYZz1ukN5Szf4EhQ78j7/POIXL8IQGkrGnTPs/nrLly9nxowZ0AU6zStmrZpSlFIu1XpxvQUcxjx3mop5qWYg5uK0Ydrnt7Upx05n6TSvKIqinH2ErYat2sk1Joa4Lz4n+9F/QJ19X8uybdmCBQu6Qqd5hbYVzW8D7tJ6bk0ALsK8PLMEKACWaO0hFEVRlLNcdxjVbA+dlxeR/34N7v4ZgIbiYlwCAhwblOIU2lw0ryVVC+wQi6IoiqJ0edYNf4s/+4yQe+91YDSKs1Bb+yiKoiiKnRR/8imm6mpHh+FQWt+sFUKIOdp9fyHE20KIOdpjW4QQs4QQUyyfrZ63RTvPv5XnzhJCTNAe33KaeKzPOW2TVO3cCUKIFdrtBCHEYu32lDM9rzlqax+lW3vhhRccHYKiKGcxY3ExJV9+SeBNNzk6FIeRUm7Vkq3ZQohk7f7bmHtqlQghUjHvbbgVQAhRLIRYqZ2XCizU+nhtbcW51o+VCCH8teeiHZsCbLXsq2hJAs8Q+0qt36hl26Cp2rZBE2nSU6wlrU24LFv7vC2lXCCEiBNCjJNSrtaCOCqEuAZIa8uLK4q9nXvuuY4OQVGUs5hHUhJF771PwHXXIVza1frSpt64a/VrQJKNL7v93rfG3d+K82YCizEvsiuxcQyNtIRogpSyuYQoFZgjhEjVEqgXtedMAOZg3kkngZM3u7ZcN9nqnBQhxISm55yJ2tpH6dbWrVvHunXrHB2GoihnqaA7p1OflUXZDz86OhSH0xKchacZVUrRpu9mAHdaj0q18dwJwAzMHRSai2Er5qRvsRDiCObkCi1xKtJ6fc2nmd1ztOeWWPZgbEuyBWprH6Wb+8c//gHAmjVrHBuI0vVpy+6EYzsPKF2M99ixuCYkUPjOO/heftlJBfWO0MqRKLuRUs7VaquOAIusHtpsmQpshTOda5leTADzqJT1udoU43xgvnbO25zYp7HE6jqptu5o39rGpyeRUq6SUr4spbxLSjkcWKTqtxRFOYVKTpSznNDpCLrjdmr376fyd9WBXjMV89Sc3UgpU7WpxaYjXdO045YRtxKrx/ytbrcq2dKmGVulXQlXU0239lEURVEUxcx30iRcQkMpfPddR4fiEFp91GztM01qp5IxT+tdZ0mErJ530mOtODcZmGBZgYh5z+fmkqYJVlOS1rNzgdpqxxloG19brqutUGy8jXkEbMpprt8sx1fwtZH2jWjz3KmiKIqiOILO1ZXAP/2J4y+/TPWu3XgMGujokDqV9n69ssmxudrNrZiL6Jt7XtPHWjq3Z5PDpxTNa9OJp5OqXWer1flNr2u5PfsM12mWTUa4OtlmTh72UxRFURSn5n/dNHQ+PhS+846jQ1GaoY2+WUav7KLTEy7rJmJWx6Zox8/YgKwl+fn5HQuuE6k9HzvHa6+9xgUXXODoMFpN/Vwop6N+NtrHWb5vem9vAm64gfKffqIuLc2Wlw625cXOVtrqxJ62LJJvqtMTrmb6WkyxOl5imePVkjDrD/+Wrl1QUGCPkO3CWf4TAPOu8931dZOSkvj+++/tdn1bfw3O9HOhOBf1s9E+zvR9C7z1FoTBQOF779vysiG2vJhiP84wpTicE0VnqZiL3pBSLmnyUaKdMwEY3poETGmd7pxwrVy5krKyMrtd31Hfuy7nLN3IWOkCOvFn0yU4GL+rr6b0669p6EIzMoptCOmALd2FECuklBO1229j7mC/VRvdmiilbHMxmnatGsBodSgfcNZhr2CcJzY/oLQbv649v9e2/hqc6edCcS7qZ6N9uuP3LZgTI1t6KaW7I4NRWscZVimWcJqOsG2lfugURVEURXFGzjCluIkTqw4TgBWnP1VRFEVRFKXrccQqxSmYtwSyFMsvARKsGqKp/lqKoiiKonQrDqnhsjUhxBxL3ZeWyJUAyVaN1RxOCFGMeVHASmePtauz+r4mWJrcOdv32tJ52VLLqB07JUZni1uxH20h0ATt7vAz/T+hfi5OZfmjHas64O72vWvpva4rf21nA2eYUuwQ7ZfMskllsy0mnMRUKeWwJr8szhprl6V9Hy07EaRq2zQ43fe6Ne1RnDFuxa6mAYHaqD9CiBnq56J1tC1XJmrfE8s2LN3qe9fSe11X/trOFl064bLsZ2R1qNkWE07Cv0kHW2eOtSvbDCy27LmlbcvQFb7XzcXYFeJWbERKOd9q25EEzFuhqJ+LVpBSbpVSztZGCVO15pXd5nvXyve6Lvm1nU26dMLFqbt5+zd5PKgTY2lJIFCktcEA5461y9L6tb0NLObEnlf+TU5zxu+1f5P7Qac5pnRz2ptrkfZ/m3+Th9XPxZmlYJ5Sg+71vWvNe11zxxQn0mUTLiHEhGYK7EuwUYsJW9P+ei3BPNRrmWd3yli7Mu17u1JK2dPqfgnO/70u4dQYmzumdH9TpJQztdslqJ+LVtPeE/zP8Hvf3DGn1ob3uuaOKU7EGfpwtVeRNkftj3mVYzJO2mJCCDED81+sS4BC7bBTxtoNJFhqYIAXMdfFdIXvdXMx+jdzTOnGhBBTrAqgJ6B+LlpFCDEHOKJNyZZgTjy6y/eute91/s0cU5xIlx3h0ubsV2L+xfLXjjlri4lFWBUxWrYrwjlj7erma8XGE4Bp2sii032vW9MexRnjVuxH+3eeI4TYIoTYAurnog3exrxIZgLgf7rf+674vWvte11X/NrONt2iLYSiKIqiKIoz67IjXIqiKIqiKF2FSrgURVEURVHsTCVciqIoiqIodqYSLkVRFEVRFDtTCZeiKIqiKIqddeU+XIqiOAmrfduKtO2UFEVRFCtqhEtRlA7R9q+batUrSFEURWlC9eFSFKVDtC7fC9H2sGuy55uiKIqCGuFSFKXjJmjTiBNUsqUoitI8lXApitJuWu3WZiFEAqC2ElEURTkNlXApitIRycARzJuGq9EtRVGU01AJl6IoHdETKFEb5SqKopyZKppXFKXdhBBbpJTDHB2HoiiKs1MjXIqitIsQIhmrui2tjktRFEVphhrhUhSl3bQkawJQBGxVdVyKoijNUwmXoiiKoiiKnakpRUVRFEVRFDtTCZeiKIqiKIqdqYRLURRFURTFzlTCpSiKoiiKYmcq4VIURVEURbEzlXApiqIoiqLYmUq4FEVRFEVR7EwlXIqiKIqiKHb2/+loeVHqSGfXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set up plot area\n",
    "fig, ax = plt.subplots(1,1,figsize=(8., 3.))\n",
    "\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel(r'$\\ell$',fontsize=14)\n",
    "ax.set_ylabel(r'$\\ell(\\ell+1)C_\\ell^{BB}/2\\pi\\ [\\mu{\\rm K}_{\\rm CMB}^2]$',fontsize=14)\n",
    "ax.axis([40., 400.,  0.1, 100.])\n",
    "ax.xaxis.set_major_formatter(matplotlib.ticker.NullFormatter())\n",
    "ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())\n",
    "ax.set_xticks([40, 50, 80, 100, 200, 300, 400])\n",
    "ax.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())\n",
    "ax.tick_params(axis='both', which='both', labelsize=10,\n",
    "                 bottom=True, top=True, left=True, right=True,\n",
    "                 direction='in')\n",
    "\n",
    "ax.plot(ells_512, ells_512*(ells_512+1.)*clBB_353_0_BICEP/(2.*np.pi), label=r'${\\rm Model}\\ 0$', linestyle='--')\n",
    "ax.plot(ells_512, ells_512*(ells_512+1.)*clBB_353_1_BICEP/(2.*np.pi), label=r'${\\rm Model}\\ 1$', linestyle=':')\n",
    "ax.plot(ells_512, ells_512*(ells_512+1.)*clBB_353_2_BICEP/(2.*np.pi), label=r'${\\rm Model}\\ 2$')\n",
    "ax.plot(ells_512, ells_512*(ells_512+1.)*clBB_353_tmp_BICEP/(2.*np.pi), label=r'${\\rm Template}$')\n",
    "ax.plot(ells_npipe, ells_npipe*(ells_npipe+1.)*clBB_353_npipe_BICEP/(2.*np.pi), label=r'${\\rm NPIPE\\ Split}$')\n",
    "\n",
    "ax.vlines(80.,0.,1000.,linestyle='--',color='k')\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "order = [0,1,2,3,4]\n",
    "leg = ax.legend([handles[idx] for idx in order],[labels[idx] for idx in order],\n",
    "                     loc=(1.02,0), prop={'size': 10}, frameon=False, numpoints=1,\n",
    "               title=r'$353\\ {\\rm GHz}$')\n",
    "\n",
    "ax.plot(80., 4.4, color='r', marker='*', markersize=10)\n",
    "plt.savefig('dust_models_BB_BICEP.' + ext,format=ext,dpi=res_dpi,bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f4d9b54",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}