notebooks/08_firas_plot.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "sys.path.append(\"../\")\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.pylab as pylab\n",
    "from matplotlib import gridspec\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "from scipy.stats import chi2\n",
    "from tqdm import *\n",
    "\n",
    "from grf.grf import FIRAS\n",
    "from grf.units import *\n",
    "from grf.pk_interp import PowerSpectrumGridInterpolator\n",
    "\n",
    "from IPython.display import set_matplotlib_formats\n",
    "set_matplotlib_formats('retina')\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# FIRAS data checks and plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:855: MatplotlibDeprecationWarning: \n",
      "examples.directory is deprecated; in the future, examples will be found relative to the 'datapath' directory.\n",
      "  \"found relative to the 'datapath' directory.\".format(key))\n",
      "/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/__init__.py:846: MatplotlibDeprecationWarning: \n",
      "The text.latex.unicode rcparam was deprecated in Matplotlib 2.2 and will be removed in 3.1.\n",
      "  \"2.2\", name=key, obj_type=\"rcparam\", addendum=addendum)\n"
     ]
    }
   ],
   "source": [
    "# Plot parameters\n",
    "\n",
    "from plot_params import params\n",
    "pylab.rcParams.update(params)\n",
    "\n",
    "cols_default = plt.rcParams['axes.prop_cycle'].by_key()['color']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Where to save plots\n",
    "plots_dir = \"../paper/draft-letter/plots/\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FIRAS data checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pspec = PowerSpectrumGridInterpolator(\"nonlin_matter_bj\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "../grf/grf.py:633: FutureWarning: read_table is deprecated, use read_csv instead.\n",
      "  names=['freq', 'I', 'residual', 'uncert', 'poles'])\n"
     ]
    }
   ],
   "source": [
    "firas = FIRAS(pspec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "chi2_null = minimize(firas.chi2_FIRAS,x0=[2.725],args=(0, np.ones_like(firas.omega_FIRAS)), method='Powell')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(2.72501704)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi2_null.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 is 48.19672388442093\n",
      "chi2 / ndof is 1.147541044867165\n"
     ]
    }
   ],
   "source": [
    "ndof = len(firas.d) - 1\n",
    "\n",
    "print('chi2 is', chi2_null.fun)\n",
    "print('chi2 / ndof is', chi2_null.fun / ndof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "58.12403768086803"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi2.ppf(q=0.95, df=ndof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.01859729, 1.01219083, 1.00665943, 0.96947179, 0.99214751,\n",
       "       0.99333799, 1.01169952, 0.99478948, 0.97820668, 0.98441236,\n",
       "       1.39763349, 0.95627042, 1.03539677, 0.931326  , 0.90965366,\n",
       "       1.04189886, 0.91012744, 0.96399769, 0.94892646, 1.01679604,\n",
       "       1.01600407, 0.99166357, 0.91752175, 1.1235758 , 0.91454281,\n",
       "       1.02221851, 1.00215006, 0.99539556, 1.0102167 , 0.98502826,\n",
       "       0.97942204, 0.92540096, 0.97359247, 1.01144426, 0.97696511,\n",
       "       0.9837021 , 1.01099055, 0.94466418, 1.00737463, 0.98468032,\n",
       "       1.00800534, 1.00442703, 1.00053248])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Amazing!\n",
    "(firas.d - firas.B_CMB(firas.omega_FIRAS, chi2_null.x)  / (1e6 * Jy)) * 1e3 / firas.resid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "../grf/grf.py:170: RuntimeWarning: invalid value encountered in sqrt\n",
      "  z_ary, self.omega(omega_0, z_ary, evolve_z))), nan=1e-18 * eV)\n",
      "../grf/grf.py:176: RuntimeWarning: invalid value encountered in sqrt\n",
      "  z, self.omega(omega_0, z, evolve_z))), nan=1e-18 * eV) - m_A\n",
      "../grf/grf.py:366: RuntimeWarning: invalid value encountered in log\n",
      "  sigma_LN_sq_ary / 2.) ** 2 / (2 * sigma_LN_sq_ary)) * (1. / (m_Ap ** 2 + (b - 1) * m_A_sq))\n",
      "../grf/grf.py:369: RuntimeWarning: invalid value encountered in log\n",
      "  ** 2 / (2 * sigma_LN_sq_ary)) / np.abs(self.dz_dt(z_mesh)) * (m_Ap ** 2 / (m_Ap ** 2 + (b - 1) * m_A_sq))\n"
     ]
    }
   ],
   "source": [
    "# Get disappearance probabilities for benchmark point\n",
    "p1 = 1 - firas.P_tot_perturb(firas.omega_FIRAS, 6e-6, 4e-15 * eV)[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FIRAS blackbody plot (Fig. 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/smsharma/anaconda3/lib/python3.7/site-packages/matplotlib/figure.py:2369: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n",
      "  warnings.warn(\"This figure includes Axes that are not compatible \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAFnCAYAAAAPGfZrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnX98U9X9/1+nECnlR9NShQ4Rm4LDwZhLi6KoH4F0bjp0alpEnW4KqeCc8wepzDmY+rGm+NnHfRUkQRB/WxJ007mpCeLwBypN8AOiTElAp0ClNAGk/Gjp+f6Rey83yU2aNDe9Sft+Ph730d57ftz3OffmvO85533Om3HOQRAEQRBakKe1AARBEETfhZQQQRAEoRmkhAiCIAjNICVEEARBaAYpIYIgCEIzSAkRmsMYM2otA0EQ2kBKKMtgjBkYYzbGmJ8xFmSMeRljbsaYWQi3MsYscdKahbhexphT+N+tFJ8xZhTi+RljXDi8ssMvHDbGmF6WziQLD8pkFI+g7LAnUV4jAGc38hbrwxynHG6FPKPLKZbRzRgzpfickoov1LP4HJzCYRLktqZyTyG/pN+POM84odxC/mJcMX9Td98XgugSzjkdWXIAsALgALwAjAphXiHcopDWLYSZoq6bZHnqFdJZxHCFMLMQ5o9OK0sXVEhnFNK4kyizM7o8XeQtymRNNk1UuD+ODByAPcnnZIj3HOLEi36WdqUyZOr9EOpMDEv4TATZ3PHk6+77kotH9G+Jjswc1BPKEoSvVxsAH+e8gnPuk4dzzhsANMZJa0dY2dRyzj1R6TwAqhFWDM5UZOKcuwAEEG5MFXtfcdL5hHsaEsUTvphNAFanKFMdgPJk0yRBnfDXkmQPR4xf20U88XlGP8taAD7lJMqk834INAIIATAxxhSfi6wHE0pFNpkM3XpfshGhLmxay9EXICWUBQgvvDh0VR8vntDQRDQQwnCW+INXbMyFxkFsgMwpiifeL6lGX2zgohvJOFgArOacp9rouQAUp5gmLpzzgOw0mfmpGjFuF/NZRsRXxF0OVYqk835EIaaNpzwtqcgVh5TelyympusohBqQEsoORCUSEhRGIqIVzUzhr6+LxlzsIc1MECcCofETG9kue1GCApJ/PXbVoNUmEUfM28gYswGS0tiYTLpk85adeuJGhNQjcSD8xQ8k7g0FAOjFeZWosNUIK9NkSOf9kOOIyi+aSUl+PCiS6vuSrQjvA/WCeghSQtlBlfA3kDBWGBsiG0rxR99V2tao+AmRGQyEoDDMpxBfjxPDVACkL/N48U0IN6rJNnoRyjNR3snCGNMLcojlrE5CHlFxisrTkmAiXoxjBOAWJvP9okFCVA8sEem8HxLCR4oLYcUYoYiEekg0nJeQVN8X+X0ZY3bBoMEuGm/IwkQDDLdgdCEa3HiVDDvYCeMcp+yvQeGe7qh4ViHMIpRDj3BPVzS8MMnSBoVnaWOMWYQ8vMKHklEmc1B2T5s8nVhnUQYmelk9+IVhdrFMNkFOf5JDxrmD1pNSdHAgPJHLATgzlRbhxomHH3nEdWlCHycaVzvCP0S/kM6gkJ+YTpwoT6kMUDBISCJvW4L8kjVM4LIyihPwQaUyKuRhhDCxj3AjJeYX18AAJybrow9/MvdM9/2QyWCVlSHGQAMyowycMNRIZJiQ0vsSRy6Dwvtogsx4QiZvhNGNTA6n7Jr4jhvl5ZK/E7J0Vtk1sX71crmgYHwRJZNbyE80/rFFvRvBqHQ2RL3Hsrj+KJlE4xVnVLnjvuO5emguAB0cSNJyKU5adzJpZS91dOMT19pJCBd/oOY46aJ/4F02lMIPL+4PKTpvIb4T6iih6PKLDUgQXVh0CTKYZedupTzjpDUJjZD4rFNR2N1+P4T0khKKys8onBuiGsZklFBK70ucuGJjbo9qaOX/x1UIsnsZZfGin69RXhax7hWejTWZe0aFc9k1s/z9UXoXZXVni7ouvg8GhbjuOGXOeetD8aDhuOygSfib0JoMkIaQ5PHE4aOuJurFNMkOAYmIY+PJjPFHzEcIww1Kcz4LcGJ+okt4eBgp3QnzeHl7IMzdIMFktGjJxyPnZMS6MSgNkcjnmjjnHs55Hee8AkCFcDlZI5F03g8lRAOFBcLfWiQweEiRpN8XHh769CHc4IrDlW4k/46K8SohG2YWhq5swrCXOIxbLnseEfkLz6Y7w7vS0C3n3MVTN7CJgCsPz0YPD4v3UM0wR2tICWUH4g+2y0YG4YZDHk8cxzd2sUhQbCRTnTAWG0B01bhxzkOc82rZpUqEv7qjSdkKS1AWajWU0Yg/7KoEcSwAWuULNBE5eV2nkEaxroXGN1mjBHk+3Xk/lO4vWkuahWeqT7cBlZH0+yLIUoFwvTsQbnBNUH5nkiUgKHv5wXjYLL7bxClLq8K17qJW/eccpISyAKGB9QDSmp9EmLhs0ldo0MRexQKlBLJJ6ADnPOkeiECl7F6p9qLqEPXVKZiIN3UjL6jYUEYTY7Sh0OjUAqjg4TU60oET1nEmhY+A4gQm8a1I8os/nfcjAeJ74Ia6lmxJvy/CJL9V6InUyhSSvgvT9+h7SfUjv3/UvYz8hNFJjEIRepCJPuLU6Innutl6RiAllCVwzsWvQYtgHRPxgxCGttwA5iikrUX4R2iNthgSGkE7wg1eRXTaRESZXCt96SdKa0f4x94UFZS0WXYPotQ4SQ2zMNTmUVKCglIXryt9BDgVLNHEtV1Jf52n835AeehG6lUmqbS6pBvvix7AgiiF3yTIFD0MFbEmS6hTPQAH5zwgPJs6KFv+OXGiDmqFa9HvoPS8ZcrTIMQ1IcXFxTgxxCtHHI2gbY3kaD0pRUfkgfCLL1pu+XHiS9WGrifOzUJcMY14KG3zY0Sk5RkXzr2y60Ehr4hJ4yTSSVZgCmVLZJCglLdbuKa4hQpObEkTLY+7C3ntUflYhTheoc5MsjoKCoc3jrxBnDBu8Ir1Lcgu7h7glZUlok4z9X4oyR/1LCPeDZmcYnnE+pDXRUrvSxdlEc3j7VF/jVHlFWWxyeJ4oWD8IOQp1rP4/ittZeUWDtG6zxAVxxxVv3pZWtGwQ6xTo4IcYn2J97Ah1rLQKMtLrE+D7F0Uy20X8pQ/m7i/iVw7mFA4gsg47MRi05R6VUTfRegl+SFsV6S1PIT69NdaAKJPYQFQprUQBEFkDzQnRPQIwjh9E8+ccQFBEDkIKSGip6gF7cdFpAA7saUSEF6LlY7pNpGl0JwQ0SNEmcgSBEEAICVEEARBaAgNxxEEQRCakTPWcYwxO4/aekOY7BZXuxfzqN0AugonCIIgtCUnhuOEVf82znm57JoFQCsXNpQU4kiKpqtwhXu0AdAB2JvRwhAEQfQ9TgbQzjkviA7IeiUk7l6MWCXkjV68Jr/WVbjCfdrz8vL6Dx8+PCX52traUFAQU68po1Y+Bw8exJAhQ9LOB8i+smVbPmrVdbaVK9vyoXrumXwyWc/Nzc3o7Ozs4JzrYhJovWVDVwcEPx2Q+QlBHH80EPxsdBUe5z67SktLeapUVFSknCaT+YQfqTpkW9myLR+16jrbypVt+VA990w+mazn0tJSDmAXV2h7s9owQdiwUGlzxUoob6MeEsK6CicIgiCygKxWQgjP4SitsNdD2f9Ga5Lhihw8eBCMMcVj0aJFimksFovi9VRRKx81ybayZVs+apFt5cq2fNQi28qVbfmowaJFi+D1emPaz927dwOA4lhf1s4JMcbM/IRRgR7hHYzLxTBEzREJ1/0Ir8zXJwrnClvXM8Z2lZaWlu7atSszBeohGGPI1mfa26C67hmonnuGTNbz9773PezevXs35/x70WFZ2RMSds5NtLo+3v5jxUmGEwRBEFlAtq4TMiHsnEo8H4awl0orwsqpCcoKRY8TTtS6CicIgiA0JiuVEI9ddGpA2IFVg+xaK2NMHzVnJHpY7DK8t7Jw4UKtRegzUF33DFTPPYNW9Zy1c0JyBCXk5ZwXya5ZAYR45OJUyM4Thivco1fMCcnhnOPQoUPYu3cv9u7di2+//Rb79+9HWVkZfvSjH2HQoEFai0gQRB8g0ZxQVvaE5AhGCLUID8/ZADRyzn2c8wbGmEXY7l0PRCqYrsJ7G21tbVi9ejXWr18vKZy9e/fiyJEjivHz8vIwbtw4VFRUwGg0orKyEmVlZcjLy8ppQoIgeitKi4f64oFuLlbVmj179vAHH3yQjxs3jpeWlsY9Ti4o4CU6XcxxckGBFGfcuHH8mmuu4c888ww/cuSI1kUjFPD7/dxisXAA3Gazcbvdzu12O7dYLNztdkvx3G43NxqN3G63d+s+brebGwwG7nQ6UwpLFbvdzvV6PQ8Gg93Ow2QycavVyu12O7darVyv13ObzcZtNhs3m83carWmLWdP4PV6I+ROVCeJwhPlI9ZTMBjkdrs94p3hPPL9slqt3O/3S2EWi4UbDIZuvVOJFqtq3vhny5FrSuiTTz7ht956Kz/ttNMiFMwZBYP4tOJh/MrhI/jcykq+ZMkS3tjYyNetW8c/+eQTvnPnTv7GG2/w//7v/+ZXXXUVLy8vV1RaZ511Fl+2bBn/7rvvtC4qEUUwGIxZ3S5e83q90jWxwekuFoslrqJJFJYqRqOx20rI7/dHlNFut3OTyRQRJ5066ClEpSDidDq50WiMGx/h3V8iDpvN1mU+opIWFXU8WaLfL8553PjJQEqolyih48ePc4/Hw6urq2OUxowZM/hrr73GOzo6Usqzvb2db926lT84bTqvGVHKzygYFNFTMhQV8f/5n/9J60uVUJd4jYTRaIxoKMSeUnexWq1xFU2isFRJRwlFy2A2m2May+iv/WxE7F2KiM9YqV78fn9MucXn3FU+yTwzpfcrXUWes9v2ECdobm7G5Zdfjl/+8pd49913AYTndS677DK8+uqreOWVV3DJJZegX79+KeXbv39//OAHP8CCtR407t6FT0JBLHn2WYw75xzoSkpwOD8fDz/8MCZNmoQHHngAzc3NmSgekSahUAiBQABms1lrUXoUk8kUce7z+WKuRZ9nIyaTCU6nUzoPBAIAAL1eHxO3uLg44jm7XC7U1NSklI/Pl7yTY4fDkdFdGUgJ5QAvzbwaFxrK8cEbb6C9pQX9g0H8LK8fVk2bhmXLlqGiQnFj8G6h0+lQU1ODt99+G8uXL8eECRMAAIcOHcLSpUtxzjnnYNGiRTh8+LBq9yS6h8fjgcfjgcvlwpw5c+B0OmEwGBLG9/l8cDgcaGhoiAjz+XxoaGiQ8gyFIlcy1NbWoqKiAh7Pic1GNm7cKMVXys/hcMDj8cDhcEiNYfS9XC5XxL1cLhfKy8tRWxt2HRYKhVBeXo66ujrFMskbV1ERG43GuHWQzcjlbmxshNVqVYwXXebW1taIa4nyCQQC8Hg8MBgMqKur61IZ1dXVwWazpVyWlFDqHvXFA1k6HPfaa6/xsrIyadjtrrvu4vv37++x+3d2dvK33nqLTysriximmzhkCH/yhxP5lrnzekwWIozScInb7Y6ZhI4ejjMYDFK42WyWhqmCwWDMvIE4ye90OmPmGeRxou/PeXi4KHpeRsw/+l7RcnHOpXuLJDvs53Q6Y+4rYrPZuNPp5G63m7vdbm6z2bjX6+Umk0lx+M5gMHCv18u9Xq80zOl0Orndbk84NxIMBrnVak14dDW0FQwG45YjGqvVGncos6t8/H5/xNCdPB0ASc7oYd7ukGg4LutNtPsqnHM8+uijeOihhwCEh94WLlyI2bNnQ7aTRMZhjGHq1KmYGgjgo48+wgMPPICmpibsBfCH1n24dfRp+H57O3S6WDchucgf//hHbN26tcfvO378eNx3333dTm8ymeB2u6UekRJerxd6vV7qlYg9kNWrV6Oy8sTm8uKXb11dHQKBAKZPn461a9fG5FdcfGJTksrKSlRVVQEI92bE/+VxPR4PAoFAxL2A2KEii8WC8vJy2Gy2lHo2brdbMW5VVRXsdrvUSwyFQnA6nbBardKXfnSvw2AwSHlVVlbCbDZL6cvLyyPOo8uSbs+hrq4Obrc7qbgejyfu/ZTyCYVCUn0bDIaIHmo04hCf0+lMWOZ0ISWUhRw9ehR33XUX1qxZAwAYPHgwHn/8cUyfPl1Tuc4++2y8/PLLsNvtsNlsaG9vx1/+8hd4PB48+uijGDdunKbyqcHWrVuxYcMGrcXoFuXl5XC5XHHD6+vrMWzYsJQaE7FRra6uTtgwyhXJvn37khc6DmazGQ6HI2b+IxEejwd2uz3imjjcJC+vXq9HdXW1dG40GuHz+WA0GhEIBLoc0jSZTHHjhEIh1NfXJ5SzvLw87hxLQ0ODpFTkCiOeLPIPga7y8Xg8qKurg9frTSifiFxZWa1W1NbWJq0cU4GUUJbR0tKCG2+8EU1N4S3uRo0ahaeeeiprGvh+/fph3rx5mDp1Km699VZ8+umn2Lp1Ky6++GLU1dWhtrY2ZeOIbGL8+PE5e1+n0xm3wRbng8RGRJxLEBtVpcY7FArBYDBIk93RE9Stra0R+YthtbW10pyOPK7JZEIoFIq5V/T8k5iH2INJBnE+KNoIoampSVFhyOPV1tbCZrPBbrcjFAop9qY8Hg8qKysRCoUSzsGm0xNyuVwwm81S47969WqpTgOBAIqLiyOUks/nU1RS8fKprKyMkE2Mlww2m036yFHd+EVpjK4vHsiCOaFPP/2UT5o0SZr/ueyyy3hLS4umMiXi2LFjvL6+no8cOVKS+fLLL+c7duzQWrReS6LFqvJ5FHEuw2Qycb/fz4PBIDebzRHzIvIFruI8idvt5k6nk7/88svcaDRys9ksrcUR78k5l+KK8ZXmVcTFkDabLWLRoxjm9Xql+ReLxRIzt2GxWJKqE7vdzs1msyRf9L2U5kXEOGL5xTjiOit5GovFErNoUy3zdLk8iFr3I5+vMZvNMXNJSvNTXeUjLmSNnneTp1darCrO5cnfgVSgdUI5oIS2bNnCx44dKzXmv/3tb/nRo0c1kycVmpqa+HnnnSfJXl5ezl9//XWtxSJyHLUaelERy5ErX865pMjlaUSilZCSQiASQ+uEspx3b7wJ1ZMnI7hzJ9pb9uGaAfmY+eFGfP6727UWLSkqKirwyA/Go4rlob2lBfu/+grXzbgM88vKsWXuPK3FI3KI2tpayYxYrWEft9sNl8sFl8slmZSbTCb4fD7YbDZpKFEcvhJNyn0+H3w+H5qamqS04jxVNnkzzXVyYhftnkCrXbQ7Ojowa9YsvPfeewCAe++9F3Pnzu1RGdTE7XZj7ty5aGtrAwDMmjULDz30UK+xniMyi7hGSW6dRuQ+iXbRJiUkwBjbVVhYWDpmzBgAYTPRnvjaWbhwIZYvXw4AuOKKK/DYY4/1qAl2Jvj0009x/fXXQ1ToU6ZMwfLlyxNa+hAE0ftwOBxwOBzYvHkz2tvb93POYxoBUkICWvSEXC4Xfvvb3wIIW0e98sorGDhwYI/dP5M0Nzfj17/+NT7++GMAYTPPZ555BmVlZRpLRhBET5OoJ0RzQhqxZcsWzJ8/HwBQVFSElStX9hoFBADDhw/HSy+9hBkzZgAIm5hecskleP/99zWWjCCIbCKpdUKMsccBKK+KSg4GYB/nPHcnO1Rk3759uPHGG3H06FHk5eXBbrdj1KhRWoulOvn5+Xj88cdRXl6ORx55BPv378fVV1+NhoYGXH311VqLRxBEFpDsYtVhnPOadG7EGFudTvreQkdHB26++WZ88803AMKGCOeff77GUmWOvLw8WK1WGAwG3HnnnWhvb8cdd9yBPXv24Lbbbsv5+S+CINIj2eG4jSrcS408cp77779fsoS78sor+4yp57i31uGPQwuRHwqhvaUF/221Yu7o08mEmyD6OEn1hDjni9O9kRp55Doul0uyhBs/fjwWL17cZ3oCE5YuwQQA/7V9O6qrq9Hc3IxXOtpxcskwjO/sRF4eTU8SRF8krV8+Y+x0dcTo/Xz22We92hAhWcaMGYO//vWv0hzYihUrYLVacfz4cY0lIwhCC9L9/Mx+l4VZwOa58zDnvCn47ptvcHzfPtyKPHxdczU+mXeL1qJpwujRo/HXv/5V2ljy+eefx29/+1u0t7drLBlBED1NwuE4wZgg3sIOJoQ9obZQvY2tF16AwN/+Ct3AfMybNw83/+EPWoukOaWlpXj55Zcxc+ZMbNu2DS+//DKOHDmCxx9/HCeddJLW4hEE0UMkXKzKGCsEYAGg5KSEAbByzm/OkGw9SqYWqwaDQZx//vkIBoMoLS3F+vXrMWjQIFXvkcsEg0HMmjULmzdvBgBMnToVK1asQH5+vsaSEQShFt1erMo53w8gwDnfoXAEACi7cCQk6uvrEQwGAQB/+tOfSAFFUVRUFOHZc926dbjuuutw6NAhjSUjCKIn6HJOiHO+JkFYrM9fQsLn8+G5554DAPzXf/0XLr30Uo0lyk6GDh2KF198UVov9f777+OXv/wlDh8+rLFkBEFkGrKLzRDHjx/HggULwDmHTqfDgw8+2GfMsbtD4K75mPflf/DD9na0t7Tgnb//HTNOHQVfba8Y7SUIIg4pKyHGWGIH6gQA4Omnn8aWLVsAALfccgtt3NkFE5YuwUUfbsArX3+N6VdeCV1JCT4dcBL+t+0Qjh07prV4BEFkiO70hOhzvgv27t0r+XIfNWoUbr31Vo0lyh0GDBiAJ598Eueddx4AYO3ataitrSXzbYLopXRHCZHvhy544IEHcODAAen/vrgoNR0GDhyIp556SjJWeOONN/Cb3/wGHR0dGktGEITa0JyQynz44YdwOsNGg1VVVaiqqtJYotxk0KBBePbZZ3HWWWcBAF599VXcfvvttLMCQfQySAmpSHt7OxYsWAAg7Mbg/vvv11ii3Gbo0KF44YUXMH78eADAmjVrYLVa0dnZqbFkBEGoBSkhGW1tbaisrERlZSUcDkfK6VesWIFt27YBAG677TacdtppaovY5ygsLERjYyPGjRsHAHjhhRdwzz33gDwCE0T243A4UFlZiZaWFgAoUIqTsntvxthDnPO7VZAvq0h3x4Q9e/bgggsuwKFDh2AwGPDWW2/R9jMqsnfvXlx55ZXw+/0AAIvFgoULF5LZO0HkAGq79/arIFOvY/HixdIq/wcffJAUkMqcfPLJuL98DIYd/A7tLS1Y8uCDuK3MgA3nTumzG8ESRG8gZSXEOV+eCUFyGc/11+O5xx5De0sLzmrvgK5uATWOGeCiVU/i9W2fYfTEidCVlMB17Cg+nX0jJixdorVoBEF0k2Tde8fAGBsKAJzzA+qJk5u8OWIE8oYNQx6Ah15/HRMnTtRapF7LyJEj0djYiF/84hfYt28f7r33XgwZMgQ1NWl5nycIQiPSMUx4C+Edtvs033zzDRobGwGETbJJAWWe8vJyvPDCCxg6dCgA4I477sA///lPjaUiCKI7pKOE6jnnDysFMMbOSiPfnOLRRx+VVvPffvvtGkvTd5gwYQKeeeYZ5Ofno7OzEzfffDPeeecdrcUiCCJF0lFCfsbYtDhhM9PIN2fYtWsXXnjhBQDAtGnTpIWVRM8wadIkrFixAjqdDu3t7fj1r38Nr9ertVgEQaRAOkroCQAOxthxxtgXsmM7AKtK8mU1S5YsoV6QxkydOhWPPfYY8vLy0NbWhmuvvRafffaZ1mIRBJEk6S5WrQYwBsBPoo64Poh6C3v27MGzzz4LIOwrqKKiQmOJ+i4zZszA4sWLAQAHDhzA1VdfjZ07d2orFEEQSZGOEprDOd8Ux+Nqr3f3IO8F3XnnnRpLQ8yaNQsLFy4EEF7YWlNTgz179mgsFUEQXdFtJcQ53xR9jTF2FmNsqFJYb6K5uVnqBV1wwQXSbs+Etkz5v824UncS2ltasOPjj/GzM76PNyadTeu1CCKL6bYSYowtY4w9LiieQsbYFwBqAVgSGCz0CpYuXYqjR48CCJsHE9nBhKVL8OjOHbDMnw9dSQmaBw/CX04uwekNNq1FIwgiDukMx7k553M55x8DWABgrXD+MIAidcTLPr799ls8/fTTAIApU6bgnHPO0VgiQg5jDPfffz+uuOIKAMDHH3+MG2+8UfpoIAgiu0hHCYVk/08H8GacsF7FsmXLqBeU5eTl5eGRRx6ByWQCALz77ruYO3cuOcUjiCwkHSXEAYAxVgjAyDl/KTqst9HS0oKnnnoKAHDuuefi3HPP1VgiIh46nQ4OhwOTJ08GALz++uuYP38++SIiiCwjHSVULsz92ABIm5oyxmanLVWWsmzZMhw+fBgArQvKBfLz87Fq1SpMmDABANDY2Ij77ruPfBERRBaRjnXccgDlALyc85sF44SHEF43ZFBLwGyhtbUVq1atAgCcffbZmDJlirYCEUkxdOhQPP/88zAYwq+kw+HA//t//09jqQiCEOn2LtpApFsHzvl+AL3O2Z2I3W5HW1sbgPC6IHKmljuUlJTgxRdfxOWXX47du3fDZrOhsLAQv/rVr7QWjSD6POTeW4aSe+9P5t2Ct86ejGX33Y/2lhYYDh9Bv/l1tPYkxzj11FPR2NiI4uJiAMA999yDNWt6/cYeBKEpGXHv3VtJ5N571apV+P3vfw8AWLlyJX7605/2tHiECnwy7xZsevc93PPFv3H4+HHkMYa7DeWomjqVHOMRRAZR2703AIAx1uu35gGAzs5OrFy5EgAwatQoVFVVaSwR0V0mLF2CX27+GKvXrcPgkSPRb9gwPNJ2CMFZV2stGkH0WdIZjqvrzZZwIuvXr8f27dsBAL/+9a/Rr18/jSUi0mXy5MkRLiB+9atfoampSWuxCKJPko4SquWcP8EYm84Yu7K3OrJbsWIFAGDgwIG45pprNJaGUItp06ZJLiAOHz6M6667Dlu3btVaLILoc6Rrog3O+Vphoep+QRldqZp0GhMIBLB27VoAQE1NjeROmugdzJgxAw8/HHYOLLqA8Pv9GktFEH0L1azjOOc7ADAAvxec283O9d6ROBcEADfeeKOGkhCZ4uqrr8Z9990HANi3bx9qamrw9ddfaywVQfQd0jFMOEv4ezpj7CHGWCsAC4A6zvlYzvkT4WB2JWMs57oQBw4cQGNjI4Cw07qxY8fZbBCvAAAgAElEQVRqLBGRKWbPno358+cDAHbv3o2amho0NzdrLBVB9A3S6Qk5GWNvAPAhvFdcGef8Ys75WjGC4PTuJQCmNOXscVavXo1Dhw4BCDdSRO/md7/7HW6++WYAwM6dOzFr1iwEg0GNpSKI3k86SmgYADvnvJhzvkDYMSGGXPQtdPz4cWkorqysDFOnTtVYIiLTMMZwxe49mAqG9pYWbHn3Xfy0zADPpHNoYTJBZJC0TLSjds6OR00a99CEt956Czt37gQA3HTTTcjLo40l+gI/fHwpnv7PV6iZMwe6khL8p2Ag/veUEoy2PaS1aATRa0nbOi6JeDcnqayyhieeeAIAMHjwYFRXV2ssDdGT9OvXD3/5y1/w85//HACwadMmXHvttdLQLEEQ6pK0EmKM3ZVJQbKFzz//HO+88w6AsOXUkCFDNJaI6Gn69++PJUuW4OKLLwYANDU14frrr5fceBAEoR6p9IQaGGMbuzjeYIxdkTFpewBxcSpjjMyy+zA6nQ52ux3Tp08HAGzYsAG/+tWvyE04QahM0huYMsa2A7AnEXUmgGWCiXZaMMYsAPQI+y0Kcc7rFMJbhdNizrkjlfCouLtGjBhROnDgQBw5cgRVVVWSF1Wi73L06FHccMMNWL9+PYDwTgsrV67ESSedpLFkBJE7JNrANBUl9BDnPCl/QYyx+ZzzxSnKGZ2HRa40GGM2AAbOebUYDqCVc+4Szs2QKZquwhXut6uwsLC0oCC82/iLL76ICy+8MJ0iEL2Ew4cP45e//CXef/99AMBPfvITLF++HDqdTmPJCCI3UGsX7WR6QSK+FOLGwBjT40QPRqQegFkIA8J717nEQOH/Wln8rsJjEJ3WjR07FhdccEH3C0D0KgYOHIinnnoKkyZNAgC8+eabmDdvHjo6OjSWjCByn6SVkLAtT0KErXpOB1CYhkxA2D24M+r+ITFMUERKLsSNjDF9V+Hxbnr8+HEAYbNs8pxKyBk0aBCeffZZGI1GAMBrr72GefPmob29XWPJCCK3UXsBzN0AbEizJ8Q59wGokF9jjIlKJQCgErE9JQAICWFdhcdl6NChMJvNqYpM9AG+rLsbv2vei9PaDqO9pQUvr1qFK0aeik21N2stGkHkLKoqIc75GM75TM75ThXyilZktQAcQo9Ij7BCiUZUPF2FK9PZiQsPteH/pldhyfgJYIxJx6JFi1IrANHrmLB0CUwbP8SbX32Js6uqoCspga9/P/zv4TaymiMIAIsWLYpoN8Vj9+7dAKC43iUn3HsLvSA757xKODcDsHHOy6Pi+RFWVvpE4Zxzj8I9dhUUFJRu27YNo0aNylRRiF7CgQMHcO2118Lr9QIALrroIqxcuRL5+fkaS0YQ2YcqhgmpuPPOwMJWGwD51gVKvRwAKE4yXJHCwkJSQERSDB06FC+88ALOOeccAMDbb79NC1oJohukuk7I2WXEMGbOuSq+DwTT7HqZYYJoPbeDc14UFZcDEK/FDZfnJQvbVVpaWrpr1y41xCb6CG1tbbjhhhvw3nvvAQDOPfdcPP300xg0aJDGkhFE9qCWiTYQtjgrT+JI2ONIFmGtj11UGoLlm0E4b1WwdAtwzkNdhashG0EAQEFBAZ5++mlpTdmGDRtwzTXX4ODBgxpLRhC5QUpzQoyx6QjPt/gSmWyrtFjVhPBiU59wrkfYaZ6Dcx5ijFkR3kVBvjgVsvOE4Qr3o54Q0W2OHj2Km266CW+99RYAoKKiAs899xy5hCcIqLRjQkQixsoA/BjhRv4thfDCeP6FkszfAMCvEBSSD7EJiiWAsGKMt21P3PCouKSEiLQ4duwY5syZA7fbDQCYMGECnn/+eZSUlGgsGUFoSyIl1L87GQq9oB2A1DsqRLh3tFMI77YCEtIHAHS5WjSRUkkmnCDU5KSTTsLv8gci2HEcG0JBbHr7bUwrM+C+sWdgzORzMGHpEq1FJIisI+11QpzztYK/IMYYm5aLnlQJQi1+bF+Gl3Z9g+tuuQW6khK0DBmMP/Hj0N32W61FI4isRM3FqoUIe1H1MMYaVcyXIHKK/v37489//jPmzp0LAGhubsYvfvELaU0RQRAnSFsJCfvFNQFYCyAIoJxzPjNtyQgih2GM4d5778Uf/vAHAMD+/ftRU1ODt99+W1vBCCLL6JYSYoydxRhbxhg7jvAi0nrO+TDO+YJkNjoliL7CvHnz8PDDDyMvLw+HDx/GDTfcgL/97W9ai0UQWUMqOyYMFXo92xHu9YQAjOGcX8w5XxMV93RVpSSIHOaaa66Bw+GATqdDe3s75s2bh1WrVmktFkFkBan0hIII93qsQq/n7gS9Hlv6ohFE7+GSSy7Bc889h0GDBoFzjt///vdYvHgxcmHvRoLIJKlu2yMqF3kiFnVeBOBuzvkwVSTsIWidENETrKmZifl/fxUHBYd4FxUPw62jT0dxhZFMuIlei1rrhOyc8+XJRGSM5ZQCIoie4qrVjTjL78d1112HL7/8Eu8B6PjeCKz87we0Fo0gNCEVz6pJb8PDOb+7e+IQRO+nvLwcf//731FZGfav+OGHH+LnP/85duwgmx6i75GUElJjASotYiWIEwwbNgxOpxOXXXYZAGDHjh249NJL8eGHH2osGUH0LMn2hGpVuJcaeRBEr2HAgAFYunQpbrvtNgBAKBRCTU0NXnrpJY0lI4ieI9k5ofI0d0FgCLuBIAhCRl5eHurq6nD66adj/vz5aG9vx29+8xvs2LEDd9xxBxjrcgtFgshpcsK9d0/AGNtVWFhYOmbMGACAxWKBxWLRWCqiL/Huu+9i9uzZOHDgAADgyiuvxOLFizFw4ECNJSOI7uFwOOBwOLB582a0t7fv55xH+3gjJSRCJtpENvDaNdei7tVX0Hz0KACgrKAACwxjMJZ24SZyGDU9qxIEkUEuff45vL1jBy6cMQO6khJ8XVCAe460Yd/MGq1FI4iMQEqIILKMkpISvPjii5gzZw6AsMHCNddcg0cffZR2WCB6HaSECCIL0el0+NOf/oQlS5YgPz8fnHPU19dj9uzZ+O6777QWjyBUo9tKiDFWr6YgBEHEcsUVV+C1117D6NGjAQD//Oc/cckll2D79u0aS0YQ6pBOT6iOMTZbNUkIglDkzDPPxOuvv45p08Lrvbdv345LLrkEr7/+usaSEUT6pKOEajnnTzDGpjPGrmSMnaWaVARBRFBYWIinnnoKt99+OwDgu+++w4033oh7770XRwVLOoLIRbqthMTNTDnnaznnLwHYLyijK1WTjiAIiX79+uFnO3bCWlgEXTCI9pYWLLPZMOWU4Xjtmmu1Fo8guoVqhgmCbyEG4PeMsS8EB3jUOyIIFZmwdAl+99lWrN+2Def85CfQlZRg16AC3PrBBjz77LNkPUfkHOkYJpwl/D2dMfYQY6wVgAVAHed8LOf8iXAwu5IxNlQleQmCADB69Gi8/PLLuO2228AYw5EjR2C1WjFnzhyEQiGtxSOIpEmnJ+RkjL0BwIewU7sywdX3WjEC53yTMFRnSlNOgiCi0Ol0qKurw+rVqzFixAgAwD/+8Q9Mnz4dGzZs0Fg6gkiOdJTQMIQd3RVzzhdwzvcrRSIXDgSRWaZMmYK1a9fi4osvBgDs3r0b1dXVaGhoQHt7u8bSEURi0lFCS4VeDgCAMVbGGLtKmAs6XRaP9hshiAxTVFSElStX4qGHHsKAAQPQ2dmJRx55BD/72c+wZcsWrcUjiLiko4Qiej6c8x2c8zXCXJBRdv1mubIiCCIzMMZg/OBD1BcW4XuH2tDe0oL/W78epspK/G7yZDLlJrKSdJRQIjOcSWnkSxBEN5mwdAmu/tiHDS17cXd9PQaOGIF+w4Zh9VdfoaqqCl6vV2sRCSKCpJWQMNw2nzH2JmPsCwALBFNs+bGdMbYPwMbMiUwQRFfodDrceeedeOONNzBx4kQA4Z0WLrvsMixcuBBtbW0aS0gQYbrlT4gxZgfQCsARHSasF8o5yJ8Q0Vvp6OiA3W7H4sWLcezYMQBhE++HH34YU6ZM0Vg6oi+QCX9CVgABYR4o4khPVIIg1KZ///645ZZb4PF4UFlZCQD48ssvUV1djdtvvx179+7VWEKiL9MtJcQ53y9u29ObaGtrQ2VlJSorK+FwxHTyCCKnGTNmDO7/3kj8amABWGsr2lta8OySJag87TQsvOgiMucmVMfhcKCyshItLS0AUKAUJyPuvRlj9ZzzBapnnEFoOI7oS3z11Vf44x//iDfffFO6dsYZZ+CBBx7A+eefr6FkRG9EleE4xtgb8oWnohGCwtGK8HAdQRBZymmnnYZVq1bh2WefRVlZGQDg888/R01NDebMmYOvv/5aYwmJvkIqw3EeAAHZeQBAJQCDwrFGLQEJgsgc06ZNw7p163DPPfegoCA8WvLaa6/hwgsvxJ///GccOXJEYwmJ3k63h+MYY2XxDBEShWUrNBxH9HX27NmD+++/Hy+//LJ0bcSIEbjzzjsxc+ZM9O/fX0PpiFwm0XCcqnNCws7aAc75AdUy7SFICRFEmOevMuPPHg++PHxiLVHpgHzMOf983PLa35GXp5oHGKKPkAkTbTDGljHGHmeMncUYKxQWsNYCsNCmpQSRu1yzxoUP9rXgsWeewek/+hF0JSVoGTIY9f/3MX76059i3bp15LeIUI10PmncnPO5nPOPASwAsFY4fxhAkTriEQShBf3798fMmTPx3nvv4b777kNxcTEA4JNPPsG1114Ls9lMWwARqpCOEpJ7zpoO4M04YQRB5CgnnXQSZs+ejQ8++AB33XUXBg8eDADYsGEDZsyYgeuuuw4ffPAB9YyIbpP2BqaMsUIAxqidsumNJIhexODBg3HHHXfggw8+wBUjTwXfF17s+sbq1Zhx4YW4cNgwOGZchs7OTq1FJXKMdJRQuTD3YwMg7Z7AGJudtlQEQWQlxcXFWLLxI3h3BDDnrrsweORI6EpK4M/PxyJvE0wmE9asWYOOjg6tRSVyhLSs4xhjcwCAc75c6BGJuyRsF/wK5QxkHUcQqbN371488cQTWLVqFQ4ePChdHzVqFG6++WbMmjUL+fn5GkpIZAM9ZqKdy5ASIojuc+DAATzzzDNwOBwRG6IWFRVh1qxZuP7663HaaadpKCGhJT2uhGjvOILomxw9ehT/c8kleP6jj9Ac4cmV4byyMtz2+FJcdNFFtNaoj5ERJcQYOx1AA4AfRwcBKOOc9+tWxhpBSogg1KOjowP/+Mc/8OSTT+LDDz+MCBs9ejSuv/56zJo1C3q9XiMJiZ4kU0poGQA3wnvIyU2yGYCHOOc13cpYI0gJEURm+Oyzz7Bq1SqsWbMmwqPrgAEDcPnll6OmpgaTJ0+m3lEvJlNKaDrnfG2csB9zzjd1K2ONICVEEJnlwIEDcLlcePLJJ+H3+yPCRo4cCbPZjOrqahgMBo0kJDKFFkpoaK7tH0dKiCB6hi1z5+G99e/gH3u/xUf7QzgutEFMp0P/wkJUVFTAbDbj8ssvp+G6XkIm54QMnPO3FMLIMIEgiC4JBoP429/+BpfLBZ/PFxGm0+lgMpnw85//HCaTCUOGDNFISiJdMqWEtgMoBlCI2Dmhwlw0TCgsLCwdM2YMAMBiscBisWgsFUH0HbZv3w6XywWXy4Xoj0GdTocLL7wQl156KS6++GIUFdH2lLmAw+GAw+HA5s2b0d7evp9zHtO1TUcJvYnwrtmt0UEAHGSYQBBEd9g8dx42vPMO1u3bhw2hENqOh3dfEIfr+vXrh3PPPReXXHIJfvazn2H48OEaS0x0RaZ6QnGND8ipHUEQatDe3o733nsPr732Gl5//XXs27cvJs7EiRMxdepUTJ06FUajkZzvZSEZXawqzA2Bc76TMVbIOd+fVoYaQUqIILKb48eP46OPPsKTv7kV/9q2Da3txyLCmU6H4lGjcOGFF2Lq1Km46KKLUFpaqpG0hJxM9YQKAaxFeE7IwzmfK1y/CmFfQ2QdRxBERujs7MTHH3+MN954A2+//Ta2bNmiGO/MM8/Eeeedh3PPPReTJ0+W/CIRPUumlFA9gNWc803R5tqMsdm0gSlBED3F3r178a9//Qvr1q3D22+/jb07d4K3t8fEKzv5ZJhuuAGTJ0/G5MmTccopp2ggbd8jU0roKs75GuH/aXJTbXlYrkBKiCB6B8ePH8fmzZuxbt06rF+/Hps2bUK7gkICAIPBgEmTJsFoNMJoNOL73/8+zSllgEwpoStFR3bUEyIIIls5cuQIvF4vPvjgA2zYsAFNTU04duyYYtz8/HxMnDgRRqMRZ511FoxGI0aOHAnGWA9L3btIpITSUfnlsh6QpMkEp3a0zJkgiKwgPz8fU6ZMwZQpUwAAx44dg2vWNfjI68Un3x3EF4cOoe34cQBAh06Hjz76CB999JGUvri4GD/4wQ8wfvx46RgzZgx0Op0m5eltpOvUbjWA6QivFQoBMCA8TzRXHfF6jkz3hP7zxAp8vWJlzPVTb7oRo2bflJF7EgTRNZ2dnQgEAvD5fNi0aRM2bdqE/3vvPXREuKIII65V0ul0OOOMMzB+/HiMGzcOY8eOxdixY3HqqafSRqwKZNpEuwyASTj15Nr6IJGeGo7bOu83AIDxSx/L6H1yiUAgQJtWElnFkSNHsHXrVvh8PnzyySfYunUrvvjii7hzSyIDBgxAeXm5pJTGjBmDsWPHYvTo0SgoKOgh6bOPTA3HAQAEpbM83XyIvovdbofNZtNaDFVwuVwIBAKwWCy0+WYOk5+fj4qKClRUVEjX2tvbsX37dmzduhXrGhqwzR/AjsNtONjRIcXp0Onw6dGj+PTTT2PyPOWUUzB69GicfvrpGD16NMrKyqTzoqKiPjvv1G0lxBg7i3P+sey8DIARQBHCPaKd6YtHpENDQwMAQK/Xo7i4GAaDAR6PJ6KBDIVCcDgcMBqNAMK9kuLiYpjNZimfQCAAu92OhoYG2Gw2Ka3X60VVVZUUNxAIwGazweFwRMTz+/1wOBxYu3atdB+RUCiE8vLyLtM3NDSAcw6fz4fGxkY0NDTAarWitrYWAOLKV1dXF7eX5fF44PP5YLVaY8J8Ph88Ho+Ul8FgkJRLIsxmM8rLy2E2m5NWQtQTzA10Oh3OPPNMnHnmmdI7zznHvn37sH37dnzxxRfSsWXtWjSHQhHpv2lpwa4vvsDGwsKYvAcOHIiRI0fGPUaMGIH8/PweKWdPk4513F2c84fjhEmWc7lCbxqOC4VCqKiogNPpjGj0A4EAysvLEQwGodfrEQqFMH36dKxduzaiwXQ4HPB6vbDb7RF5FhUVIfp9YYzB7XbDZDIljCeOucuVGxBWlNFKUSm9y+WC0WiEwWBQjKN0zefzoaKiQipvND6fD9XV1TG+bUKhEObMmQOn0xkhJwBFhRWNWPfJKBafz4dAIBBTL0Tuc+jQIfj9fvj9fuzcuRNffvml9Pebzz9XXMckzjkpMXToUIwYMQLDhw/HiBEjcMopp0jnJ598MkpKSlBSUoKhQ4dmXa8qo8NxcZgEQHMlxBiz4MQGq8Wcc4eW8vQUc+bMQW1tbUyvw2AwRHzJV1dXo7a2NqaBtlgsqKiogMfjkZRLPIxGI3w+X9x4oVAIer0eBoMhZqt+ANi3b1/CHoOY3mQyoampKaUeg1j+pqamGPlCoRAMBgOMRmNMOQOBAAKBQER8i8WC1atXJ33vZAiFQqivr8fMmTNVzZfIDgYNGoSJEydi4sSJMWGHDx+OUEq7du3CN998Ix0tLS3o2L8/QlHta2nBvkAAnyZQVEC4x1ZcXCwppZKSEhQXF6OoqCjuMXDgQM0UV9JKSBhuMwOoAlAGoJgxVhsdDeHhuDmqSdhNRAXEOXcJ52bGmKUvKCKXyxV3jqWurk5q9D0eT9x4lZWVsNvtCZWQz+dDKBRKOEQlNvB6vT7ma9/n82HSpEkJy+JwOGC1WiVFlAoulwsWi0UxXSAQgNFoRG1tbUw5jUYjQqGQpKRF+eOVUxy6E9OFooZhPB6PdM9QKCT1ppqamhAKheB2u9Ha2gqTySQNmSrFJ3oPAwcOxLhx4zBu3DjF8CNHjmD37t2SUtqzZw+am5ulY8+ePfjm3/9WtODr0OnQ3N6O5ubmpOXR6XQoLCzE0KFDY/6Kx+DBg6VjyJAhGDRoEIYMGYLBgwdj0KBBGDRoEPr1S92DT9JKSDBAWAxgMWPMjnAPI6ZBzyLruFrOuTSryDl3Mca8UJC5NyH2NuLtkSX2JMR48XoW5eXlaGpqirkuNpChUAiNjY1wOp2KPRmHwyF96e/YofxKNDY2xlWC8vSpNMKifE6nE+Xl5RFDikqYTCZUV1dLPS4Rr9eL+vp61NbWSsNl8uE5EXHozuv1Stfq6uoi4tTW1sLr9cJkMqG2tlZSzCaTCU6nM2JeLVF8ou+Qn5+PsrIylJWVxY3T2dmJYDCI5uZm7N27Fy0tLRGH/x//QMu33yLY3o6DHcdxtDO8Fkoc8pP3tNoBtO3ejd1IPCTYFQMGDEBBQQEGDRoU8Tf6w0xOd4fjrABqskjhRMAY0yO8ZikaI2NMzzmPXyM5jjgE1draqqgcxMZWbogQPWwHhI0BlBSUvDE0m82oqqpCbW1tTC9H7DXI84hu6BPRXYeConwmkwlVVVUAYudxfD4fmpqaJCVbWVmJ1atXR9xTr9fDZrPBZrNJikbsNclZvXo1KisrI65Fl9Hr9UrKsbW1NWaoL5pU4xN9k7y8PAwbNgzDhg1TjvC//xtxeuzYMQSDQQSDQYRCoYj/W1tbcfDgQezfvx8HDhzAgQMHsNvrw4ED+3Gwo0NywQ4oKzGRdgCHdDoEo5TYUYUem0i3lJDgriGuWXYWuPeuRKyzPSC8oLYSgKdnxelZTCYTXC6XYg/C4/FICsNkMknDSErxkjGbrq6uRl1dXdyJdfl1+b1dLldScyHpeLetqqqC2+2OqYfW1taIfCsrK1FdXS1d8/l80jwWEFYqy5cvx/Tp0xXv09qq9KqFEY0/REOFjRs3xo3r8/lgMBiSjk8QqXDSSSdh+PDh3XICePToUXz33XcRx8GDB3Ho0CHp/7a2tojj0KFD2PWvf+Fgyz60cI5inU7Rr0Y6JtqnA2gA8OPoIITnjLRUQnpEuhwXid9aADh48GDcybmFCxdi0aJF6UvWAzidTlRUVMBkMkUomOgusRgvek1LQ0MDTCZTUhZbbrdbUYkpsXHjRilP+f+JSGetjV6vl3o78l5YdD0YjUa0trbC5/NJ/9fX10cMvwUCAcUhsZqampjekTx/j8cjmcfLw+RzZaISE3s9ieIThBYMGDAAAwYMiN/rEli0aBHq6+sVw1ra279Tup6OddzdABoB1COywWcAHkojX7VIufUaMmQIDh48mAlZehS9Xg+/34+6ujpprUtxcXGMcYBer5fmPkQDgdbWVhgMhojeg7hOCAgrKPk6nOLiYilMXOcTHc/v98PlcklzJaFQSPFllt+nrq4OVVVVMQ1vtCxmsxnFxcVwOMJTfQ6HAzU1NZIhgdPphMPhkBp2UQZRAQAn5pHmzJmDBQsWwGAwYObMmVK61tZWhEIhxZ6hOGzX0NAgGSaI8i9fvhwmkwmNjY1wuVzQ6/Worq6OMN+ura2FzWaTTNCLi4sTxieIbGbRokWKH+uCibZi45rOOqGInbOjwuK6/u4JGGMmAHbOeXnU9SCAas55zHBcb1onlO04HA7JEowgiN5PonVCmdppz991lIzSBEDJPEwvhBEaEs/ogSCIvkc6SsjPGJsWJ0zL+SAI1m+tgpWcnIBWlnH/eWIFNpw7BQc2bcKBTZuw4dwp2HDuFPzniRVaiKMZoVCoy7VBBEH0HdIZjtuOcG+jELFzQoWc89RXLakIY8wKICQuThUWryLeYlVyakcQBJEZMrVtTwDh3ROiLc4YsmBBKOe8gTFmEeaH9MI1zeUiCIIgTpCOEqqLt1iVMVandL2nIaVDEASR3SQ9JySsC5LowvotPU95BEEQRJ8gFcOE6M1K1YpLEDmD6H9JjritDm2vk134fD64XC64XC6tRSESkIoSqmOMvZHMgfDecgTR65AvlhWprq5GRUWFoqsKQjvq6+thNptp/70sJxUlFAKwKcmj124QSqiP6B1VbTweT4R7ZpFAIICGhgZ4PB40NDQk3OE3GnFXAzkLFiyA0+lM2jGdy+VCRUUFysvLpd0aRBwOB4qKimJ24tYal8sFh8MBl8sVI7MaqP2sHA5HhDdcWpeWxXDOkzoAPJSJuNlyANhVWlrKiZ7F6/Xy8GuoLm63O27eJpNJ+t/v93OLxZJS3vL0nHNus9l4MBhMKR+32831er1imN1uT0meTON0OrnNZuOch+vLaDSqmn8mnpXVauUWi0V6LsFgUDV5idQpLS3lAHZxhba32+uEehu0TkgbXC4X5syZg2AwmJH8GWMRLr8DgQCqq6sj/P8UFRVJ94+e7wHCvpnkvRxxd+5oGhoaYDAYku4RFRUVwel0RuyP5/F4UFlZmdbGrWpTXl4e4wI9E6j5rOrq6jBp0iSYzWbJNxU5B9QOLdx7E0SXuFwumM1mzJmT2BFvXV0dFixYENEwiz6BUnX14PP5FB3+BQKBGPfnySAO+0TviJ0MFoslxqtrvN26tUI+zyXuMp6IbHlWVVVVEfNA2aTUiUgytXccoRGhUEjasVp0O+1wOODxeKTzhoYGzSdqxYYkGRYsWBChqETjgO74GlJy9ldcXJz0XIPL5UJTU5NkcVVTUyPVud/vj5FJrG/5LuIitbW1EZZboVAorkfcrnA4HNKhJoFAAHq9Hi6XCwaDQZqfiUe2PCuTyYRQKASPx9OlC3pCW6gn1MtoamqSehczZ86EyWRCIBBAVVUVvF4vjEYjAoFAXKd3okvtZBzaiXETUV5ertgA+Hy+pIetRKdy1dXVWLBgAex2e5duuzOF2WyOcYchnvth+uYAACAASURBVEeXR3RtIQ7dBQIBOByOCK+zRqNRuiZ3+pcKtbW1KC8vjzvclM5zEv0siXJZLBaUlZXFHT7Npmcl1kc29SyJWEgJpcn6LV/inU++AgDcM+uCiLDHXvkI+w8dxcSyUzBj8vel6182h/DsW1sAANdN+yFGDz/xtffqB//G5h3fonDQAPzmsrNTlkdUOpWVldLQidiIyP0AVVdXK6ZfvXp1zJeu6NMmGtGXTqp0x0GbXq9HbW0tpk+fntb8kdKXdDxX6Oki+gSSD2nJ5zcASC7DLRZLSkN5IuK6JbvdHtdFRnefE3BCUcrzCoVCCXuyufisCO0gJZQmLQfasO0/LYph/l2taDlwGCcXFkRcP3S0XUpz6Gikj/bdrd9h239aUDJ0YLdl8ng8qKqqks7dbneE0vF4PIpfp6LyStY8OJ0v7NWrV0fk05WPoVAoBLvdjrVr10oNd3cQPadGkwkT3n379kU04krzKTU1NaitrZUMElJF7PkmGm5K5zkZDIYYRSA6SUx0v1x7VoSGKJnM9cUD3TTR/tfmnfyB59fzB55fHxP26N8+5A88v56/smFbxPWde4JSmp17Ik1HX9mwjT/w/Hr+6N8+TFkWEbPZzL1er3RuMBik/71er2Sy6na7I9KJpsHy+G63Oyae2qALE+1gMMjNZrN0Lpahu3nLTYy7Y6KdLF6vN8acW/5cRMxmc1yzZ7/fn/Aefr8/om4ygclkkkycuzLRztVnRWSWRCbamjf+2XIA2FVYWMgrKip4RUVF1q3VSAV5wxfdKPj9fm61WrnT6YxI43Q6ud1u5263mxsMBqnRyWQ9BINBbrPZOAButVrjNrhWqzVmnYfX600om9vt5larVcpbrkj9fj+32Wzc6XRyq9WqTmESyGGz2bjb7eZOp1NxvYoYRwmz2dyljHa7XXp2SkouXYLBILdardxutys+Czm5/KwI9bHb7byiooLrdDqOsGudmLaX1gkJ9PV1QvJ5GnFSuStzXKJn6M4cGkFkE1q49yZyBJ/Ph6qqKmk+IhQKIRQKobGxUWPJCJFUthUiiFyDekICfb0nRGQn2biDAkGkCu2YQBA5Cg3DEb0dGo4jCIIgNIOUEEEQBKEZpIQIgiAIzSAlRBAEQWgGKSGCIAhCM0gJEQRBEJpBSoggCILQDFJCvYxFixZpLUKfgeq6Z6B67hm0qmfaMUGgt+yYwBgDPdOegeq6Z6B67hkyWc+0dxxBEASRlZASIgiCIDSDlFCaOByOrMpHTbKtbNmWj1pkW7myLR+1yLZyZVs+apGyPEpOhvrigW56Vq2oqEg5TSbzQRdeSlMh28qWbfmoVdfZVq5sy4fquWfyyWQ9J/KsSoYJAoyx9ry8vP7Dhw9PKV1LSwtKSkrSvr9a+ezevRulpaVp5wNkX9myLR+16jrbypVt+VA990w+mazn5uZmdHZ2dnDOddHxSQkJMMbaAOgA7E0xaQGANhVEUCufIQAOqpAPkH1ly7Z81KrrbCtXtuVD9dwz+WSynk8G0M45L4iOTEqIIAiC0AwyTCAIgiA0Q1UlxBg7izE2VM08CYIgiN5Lt4fjGGPLAHAAdgA7ADQB8ADwA/Bxzt9SS8ieII05IYIgCCIxceeE+qeRqZtzvgYAGGMPAVjLOZ8rnF+VRr5aoROs41IyD2lra0NBQUy9poxa+Rw8eBBDhgxJOx8g+8qWbfmoVdfZVq5sy4fquWfyyWQ9C9ZxygmU7LaTOQBMl/2/EcCVSmG5coDWCcWQbWXLtnzUqutsK1e25UP13DP5RNfzV8uf4O9PPi/m+Gr5EynLk2idUDo9IQ4AjLFCAEbO+UvRYQRBEERuMmr2TRg1+yZsnfcbAMD4pY9l5D7pGCaUM8amAbABWC5eZIzNTlsqjWhra0NlZSUqKyuT3nrCYrGocm+18lGTbCtbtuWjFtlWrmzLRy2yrVzZlo9ayOVxOByorKxES0sLEF4/FENa64QYY3MAgHO+XOgRLRCCtnPOn+h2xhpArhyIVKG67hmonnuGePWsRk8okSuHdIbjwDlfLvt/P4C7AWmIjiAIgiASkqnFqndnIlPGmIUxZhaOpPugQnx7JmQiCIIguk9SSogxtowxti/JoxWAVW1BBaXTyjl3cc5dAFq7UkSMMSNjzAagCkBxV/fozpxQtrFw4UKtRegzUF33DFTPPUMm6lm1OSFhHZAb4QWpXUYH4OCc16QgazIyeDnnFV1di5PWDGAm57w6QZxeMSdEEAShJtkyJ2TnnO9I9oaMsbqkpUsuPz0Ag0KQkTGm55yH1LwfQRAE0TMkNRyXrAJijP1Y2C1BbcOESgCtCtdDQhhBEASRg6hqmMA538TDW/morRj0CCucaJQUE0EQBJEjpKWEGGPzGWNvMMY2yo4vAHQ5T9MN9BnIM4KDBw+CMaZ4LFq0KNO3JwiCyGkWLVqk2H7u3r0bCDvNi6Hb64QYY/OFf+9GuOcjGi0YEN5NW03izfl0afGWCkOGDMHBg2o5JSUIguhbLFq0SPGDXTBMUGxc01msGhIXqzLGijjnm4Trm4TtfNR05dAEZYWjR3IWewRBEEQWks5w3D7Z/0WZdGYnWL+1ClZycgJkGUcQBJG7pKOEGAAwxuoRHn5zMsaGCMqoSg3horADkNYeCQtVbbJzI2PMmYH7EgRBEBmi28NxnPM1wgam5Qi7blgMYL/wv6rrhIT7NQjb9pggGClwzuXbGhQDMMnXDTHGDABqAZgAGITdEzYKOy7EIO6YAIR3gs223WkJ4j9PrMDXK1bGXD/1phsxavZNGkhE5Co98S45HA44HI7M7aKtmCFjZQD2cc4PqJpxhqEdE4hcItM+Xoi+Q1fvUqZ3TFB9A1NhYevyLiMSBEEQfZ5uDccJrhoSWaUpbbFDEARBEBGk0xPaj7ChwE9kx80AHMjRrXR6wy7aBEEQ2UIyu2h3qyfEOd/PGJsjWxsksgPAWsbYlQCiw7KegoICNDXRsiOCIDJLXzEwEQ28hDmhNqU46VjHJVIytHaHIAgiDqNm34RRs28iAxOk6d47AUaou2MCkWH6ypcZQfQVcuU3nY5hwg6E1wRFBAnX5qQpF9HD0JcZ0VPkSuOY6+TKbzqdOSEPFJQN53x/2lJpRHcXq6b7o6IfJdGXyIXGkX6TXZNMHSWzWDWd4bj6XFY4SnTXMCHdH1Uu/CiJ7CD6h7/h3CkAqHFUm3R/k31BiSVTR1oaJhAEkQHEHz6RHplWEvRhmTwZMUxgjDVyzmdmIu+eRo2XtS98FfUFknmOWj/rbHhfta6DZGToDUoiG+pZDZJSQoIDu5ouI57A2D1xupTDghMuvYujNjBVLY0cNV7W3vDCp0su/GDUaLi0ftbZ8L5qXQfZIkOm6S1lTLYnNAzhnRDkEyYG4fBFxa0C8FD6okUiKhNxB2zGmJkxZkmkVFJN097SIo2xN5/1I/zi8aVqFyNnSVeJqPGD6UqGbJCRIIgTqGmY0Bg9B8QYK+OcL1aIK+6YoDa1nPMK8YRz7mKMeRFWjqqkYTodzt3wnmoC5xJdTXhnQwPdlQzZIGO2kws9UqL3oJphQhwjhES7IqhqNSd4VFXaFNUo9x+UbpreTLJDTZnKn8gOSFET2UY6hgkGxN8VoTCNfJWoxIl5HTkhIcyjUppeS6YbH2rc+iZtR9rR0dmJ/nl5KMjXSdePHOtAc7+B6ASD4Wg7Bg44Ebbx8ZXY/urr6Nd5HKP275F63d+78UZ8fd5PsNfjQd6ra3DyoaCUZl9BIb77eQ30F5yPC384GgWy/N7Z8iX2HmhDydACXPjD0dL1g4ePwT3wVADAwN1BGEqLpLBN23djk38PWgpG49K2LyPK9Nxbm3Gg7Ri+f+owDJdd/8/e/Xj5vW3gAMznn4mRJUOlsH9u3I5Pv/wWQwcNwJyfSYMv6ATwwuAzMLTxXZw//jScP+E0Kexj/x68MPj7AIAF3x2BfnC+FPaU+2N8tnELhn66CVP9H0rXvztpIDwX34D8U0/FL6aciR+Xj5DC1m4K4O3NOwEA998wLaJMzw8eC/7kWlSMLcVV5/9Auu7f1Qr70PEAgN9F1VHjv7bCt3038gafgeu++1y6zjlH3Ypw81llNEDuIOijf38D5/pPw2W6+nwUDxmIrui2UzvG2FUIN+YPcs4Pyq7fBSDEOX+iWxkr38sMYIF8aE247kd4yC1GoaSahjG2a1D//qWHOjoUZZg1Ziz0ej1OufSSuHJ++9o/ACBunHTD1SDde3S3DPt9PhzwxXaohxp/jEJjanYsma5HNfJX61l2dnIc7+xEJwcG6PpJ1znn2Lu/Dcc7OzF44EkYMnCAFLb/0BF88e5H6ATDj6efB13/E+k++2ovDrQdRd6XOzC646AkX9vRdnzw2dfo5BzjR5+M0uIhUhlaJ56Nf3+9D+Acl583LuJZ+kaOx7H+OpzyXSvOGj1Mepbfhg5hvfcLAMBFlWegpPDEdEDT57vw1d79GDRAh4srx0SU6eX3twEAzhxVgmFbNkp1+GVzCN7tuwEAUwYcxrGPT7xLnw4fiwP5gzAsvz+mXfAj6fqhI8fwz/fCDeKkCadj9HD9iTRf7cW2//z/9u48OIrrTuD49+nW6BZIIMQpgYFw2BbYjo0xSRDxuSlvgJDE6yRrjEicaze1NsbJbpFKbdmi/Ed2t8obyRvnqpDE4Nh72A4BbAfiOI4F2CbYOLbEoRgZJEtC6J7RvP2je2a65xAjzYx6Rvw+VSqY7umenp7u/vV779fvdeAdHOKmwbO23+k3Te/RP+RmTnkRs/7yun8bznf38fvjZwBYvXQ2ZUV5Id/JlZ3JLZbvdP7Z5ziYM4O0nGwWz5rK4tll/nmnz3Xz2p9PAXDbqiW2IH7w2Gk6evopK3Sx8Mxx/zb0DQ6z93AzACvmVzBnWrH/d+pYdg0nWjsA+PSqxbZteDV7Gu5cF3PKi1ixIBA22rv7+J35O61ZuSDsd8oa6OPaofP+fTTa7/Tci3+gfcgo27QefZGR4UEserXWBQSJdXjvauCCUsp3y1IM7NRaPzre9Y6i+NJviW0ZV1YWvW532Hn+O/wdOyIuf/x8x6jviTTfX5WVYR6Ee/cBianKutQ2xrp8rOsfzaX2U7z2Yyy/04W+QT7sGeBE5yBzPL22dfym6T26eweZU17E9R+Z5Z/e3NbFT/e9zpBnhLpbV9juRv/z/5r4w1utlBe7eLjuZv90r1fzxUefAeDmVYtsF52Dx07zWNYiAP7h/o22C8u//uIQJ1o7mDH7r3yh9x3/9nVc6OcfG/YCcNetNdy0bI5/P7xz11d5+mXjovMv/3QnaWnK2E/3fY3/KFxGX2k5q5bM4su3B0ZweftMO50NzwLwlbrbWDKn3D+v8bnDHPrzGUoLctnxlVts+7hl59MA3HnDQhbtesz/Oxw6dhrP80YO1HfrPkl5cZ7/nHxm9SZOtHawcOYUvvP5m/zr6ujp54NHfolCc+9d61i1JFAKefZP7/L8n95l8P33ue/CGyy1/E7pPz9IV+8A11wxg+W7G/zb0Hy2k8Lnj4CCLbfUUD2j1L/MU79/m9eb2yjJz+Vb66/3Tz9+vgNXwSIKli1lzfK5fGz5XP+8N1rOsesHTwKw/TtbbCWhXS8c40z7BWaXF3GlZRu6ewcpf/4wALdes4Clc8v9v1P3V7fx6om/AoRswxzXHLJXrWbx7DJuXlHtn9fafoGynUY1+r0PbGZmWaB0t+9IMydaP2To4EvcOnDGf6xorSn7XyPwXLuwkrx/e9i/fZvuaeeF108CcHft0xS6jJsjs03IX1ixiuk5Ia31TqVUA4Hxg5oS1ItCpPab0gjTx7WMdrsdeQL9Uu0x0t5iuNR+Gm+7lmfES9/gMF27fk7bE4H9/LM7Pk+nq4iZq65j/Te/5F9/58UBdj75Mv3Dbj63ZimzLAHlqYZfsfeDYcip5huHfmI7nvZ5Z9J2toMFp47BO4f8y7xfWM67128g3eWib3DYtm2ZGcaQX0PuEdv0tDRFeppixKtxe7z2ZdIDw4R5Ruw1HWVFLi705VN43v452VkZXF09nYz0NKYU2qtQrqicwq3XzCc9LQ2Nxugi0nDT4FnKamupKLXf4M6cWsiG3mbS8TKrzF47/6nrF/LxK+eRlRE6nNnOe9eRpiA/N4tTuwLTP7p4JlfNryBdKbKz0m3LPLhpFQrlD44+UwtdPNhtBK4lS+zHxe3XLuD2axf4A5nVP98VCGTHdwemV88o5ZHNtYBxTr5iOSdnmH8zN98DBAIAwN9fPMGSv7N/Tut//ZD+Hz7Bnebrt9f9Egic05//xLKw21Ccn8P9G1eFbDPAqiWzWLVkVth5d/SfZsnf3h8yfVZZERv7jJKVNQABrKupZl1NNcf3/tQ2XSnF1z51bWD7LPMWzy6zlfaiEfPDqmbQORDrei6hifDBo5jII7yOeZnMqVNDsuOi6SYl0V2pXKq9ZSKCVKK/Yzy+g9er6R9y0z/kprw4cPevtWb3obc4+2oTJQd/yxUdp/zz3imbx29rbiXd5WLnvZu4fkvgsw49+TLHTp2ns6SY9ZbPyUxP4/0PjZu6iwP2i3nF6hvIevHPAFz90ou2tovsHx0guyCf6RvWU/hcL2D8lh909XJy/xtkZ2aQn5tlW991CyupnFJgW4/Ptz+3mrS0NEoLcmzTr11YyUNdh0lDU1H6Jdu8utuM2unj9/23bXpBbpbt7tlq6dxyls4tDzvvyuEPWVJTHTK9wJXNYrdRQeK7G/aZXpLP9JKQRQCoKM0POz0zI91WrWiVnhYazGI9XhOdLSo9XwQkqseEh7XW2+O1Pq11t1KqM0xWW0ukLLfxLBNONAeL0wfURCQFjLW0NtaTPtJ3GPHa7/K7ewf53bFT9PQNsWb5XGaXB+6yf3bgDfYfPYkrO5OGb97hn66U4sDRk/TrQmq/8QCFTz/u/4yc5g848NQrACGlEF8d/ZDHXgrJzc7kmitmkJOVwYwp9hJAzfwKppfkk52VQVbQRfN7X/y4/4J5/Lmf+KdPL8mPeHe7bN40ls2bFnbegsopIdPicbOQ6Av4RIj1eHX6nL6cRB2ElFJ7gXqt9Qvm6/eAcPczCiM7Lm5ByNSA0WtDo/n5dUC9ZftqMBIRNka7TLDx9qJ9KclwUibaeE7aC32DtLb30NU7yEcXVdrudL//9B9563Q7CyqncP/GG/zTB4bd7Dn0NgDzppfYglC+2TjfP+RmZMRLuqVaqigvm4z0NH/1ls/MqQWsv3ExeTmZIdVQm2++mq23rQi5A89IT+Mbd14X9jtNK8lnWkn4u/lwd+zxFo+LZ6zrSIULeCps42QQ71609wMtltctwFZC06AVoz9AOi5m+1OdUqoWM+EgqOeDUqDWWvKJYhmbRA3vfTkd8J4RLxmWi7/Xq/nFS8fo6h3kyqrprLakqL7Zco5Gs7F5QWUp0y0Xb7fHy8Cwh55+W3YNRS6j6ik3KwPPiL2UtHzeNFzZGRTkZocMdFW/uRalzAZ1Sx17WVEed96wKOx3yQ1TBZZo0v43OVwON57RiGsv2mF6R9iqtT4Z7r1KqW1j2M6ojRZAzJTrkJLZWPuKm6zicVKMjHhp7+mns2eAiin5lOQHSg7P/OEEvzVTRx/7+u3+6WlpioPHztA/5CY/J8sWhKzZQN29g7YgdPX86ZQX54W0EeRmZ/DEtz4Vtn1gQWUpCyrD5534AlCyS3TVqlwcJ0Yq3Hgmy7EQS4p2SABSSl2F0eYSNjgJ54zlpHij5Ryt7RfIzcpg7dWBTifOdff5H1Kru22FLaBAoJF+2D1CluWZluml+fQODIc0rldVlPDQZ2+kOD+Xwad+xSuWzLS8++7mCsxsI0tKqVIqYgP1pSTLSeckuThePmJt95qo32HcQUgp9QOMobwbMIb6bsKosmtWSh3xtR0J5414vbg9XnKyAj93c1sXz/7pL3T2DHDf31xjyybbd6SZN1rOMaus0BaErE8/d10csH3G/BmlfGz5XEoLc81kgkCg+O7dHwu7XXk5WYF0zi2bmbclsReYVLgAp4KJygZ16vMni1Rp24slO26f1vopAKXUI8ABrfVXzNfrR10ySSUqMSHRevqHaOu8yNuZJf60WJ+HfnSA9zsusmb5HO65+Wr/9MFhN6+9cxaA9gt9tiDkCzYDw/beI3KyMvjy7Ssozs8JyQobLY1XTC5OB3OnPz9ZpEIwTvTw3tY057XAwxHmpYxEJSaM1ciIl0G3h7ycwDMjg8MeXsitpFdl4j15zpa2+9Kbp9h98C3c+dU82BU8sgZ4tQ4puUwtdDGtJI8pBS5bIgHAxtUf4bNrltq6EfGxPnUuhHBOKgTjhA7vjVEVh1KqCKjRWv86eF6qGfaMsOvFY0wtdPFJSzuE1prW9h5yszPIz8myZU0NuT0Mu0dQSpGXk2lrAP9rew8Dw25cOZlUTgk8jXy+u4/fvXma/iE3t6ystqX0NjzbxMvHW5kxpcD2dPapJ37M72/8AgCubTvobX3Tf8cz8uorDHcY3Q21vf0XBix3RDcuv4negWFbKjMYqcSPbvlk2P1QEPRwoYifVLh7FWIixRKEqs0L7meAx30TlVL3Yk/lThluj5fnX3uPmVMLbUFoyD3Ct39sNHHdvXa5bd7+oyf55UvGE/KN37zDFqC+/8wfOdfVx0cXVfJVSzcXXb0D/M8f3wGgZv50WxDKysxAY1Sx+fjueHb9+7NkpKdRvX4N1193hX/+mns2Ud3RQ3FeDjOmbLaVbMJ34iGckgp3r0JMpFiy4x5XSm0BDpv/D35ANeUSEzweD+fPtXH2vWM0epr9bUL9Q4FOTYMzfa29kAenAfuemHcHPc/iys5EKePf4GddrqqaRqErm6K80NLIY1+/LWyqcWlBblRdpgshUsdkKDVH0yY07qEc/CtQai6A1vqUUqooQR2YJpxS6mxFRUXF2bNnGfF6bU+3Dw57ePPkOQaHPVTPKLFVrZ36oJt3z36IV8Paq+bZSiFvnWnH7fFSkp9jqw7zeo19HtzhohA+Mi5TcpOHigOiOVbNNqE2rfWM4HmxjCdUhNFxaRGw35IZ92nzdc+4VuwQaxASwilycROpJtYgFEub0IPAFq31UaXUWt9ErfWvzXahuA1qJ8TlQtqMxOUmlh4Vm7TWvuENg4tTXcFvFkIIIYLFEoSsgSe4cSPCaCHjZ3ZEusH8i/opUvP9DdG81/ew6sqVK2lslC7nhBAiFo2NjaxcuTJhD6tWK6U+YXbP4w9IZlXceIbijsgMOp1a6z3m6w1KqbrROic1h3bYZG7LaCOw+iXLw6pCCDEZRPOw6rhLQmav2l9WSn0INCilXjP/v0Jr/eh41xvBVl8AMj97D8YwEqNt3xGt9TZgX5y3RQghRJzENLKq1vozSql5QK05aX+8e9BWShUDVWFm1YQZNVUIIUQKiXl4bzPoWHtMWGtMjlsv2isJHTgPjP7pVmL03C2EECIFxX28Ya31AeLbJlRM+A5RwwWmmFy8eBGlVNi/HTt2xPvjhBBiUtmxY0fY62dbWxtAQbhlxhSElFKFSqmHlVKPKaXmRHjPVcC6MW/96OKa6BBJQUEBWuuwfxKEhBBidDt27Ah7/ayoqAC4GG6ZqKvjzB4SDpsvS4FNSql5Wusec4C7tQTabjZGWEctEM3Q391aa986IrX5RJXxJoQQInmNpU1oC0aW2gEApdQDwHYzcaAb2IlRRXYkUnKC1no/Y2/DaSJ8wCk25wkhhEhRY6mOm+oLQABa653ACmCb1vpBrfXjWuun4p0dZ2a/dZrBzqol3plxk+FhVak2nDiyryeG7OeJkYj9HM3DqlF3YKqUut98NmjUaYlglrq6fQ+n+npMsLyuAbZbqvCsy24ANoWbF/S+SdGBqVKKaH9TERvZ1xND9vPEGO9+jrUD01iz40LHksZIYIhxvTZmqQulVK0ZVPwByFQK1FpLS0qpKqVUPcYYR7VKqXrfskIIIZLDWNqEwoXISGFzO/YB7mI2Whc9ZltTSdC0FqJLghBCCOGQsZSEHlJK7bX+AfXB08zpDyRoe5NOvNqOkrENKtm+W7KtJ16S7Xsl23riJdm+V7KtJ17Guj1jCUKdwNGgvwNhph0FWsa0FUliPIkJk/VAguT7bsm2nnhJtu+VbOuJl2T7Xsm2nnixbk+8ExMe0Vo/GOV7JyRhIZ6UUu60tLSMadOmjWm5jo4Opk6dGvPnx2s9bW1tvgfDYpZs3y3Z1hOvfZ1s3yvZ1iP7eWLWM9797LlwAYCMoqKI23Pu3Dm8Xq9Ha50ZvPy4h/eebJRS/UAm0D7GRV1A2C7KHVpPARGeTB6HZPtuybaeeO3rZPteybYe2c8Ts56472dXWnqBKz0tv8vjIQPFoHckeOw5CUJCCCGcE/cOTIUQQohoSRASQgjhGAlCQgghHCNBSAghhGMkCE0CZhdFtZd+pxCpR47vyS3m4b1FUqjF6L3CN6zGfozezVPyoeFkY/Y5uE5rvTXMvDoCo/yWjta9lBjdKPtZju84Mo/ZYqAao2PobWHmT9gxLUFocujUWpcopYrjPbzF5czsnX0TxgkbMqaV72TVWu8xX29QStVJIBqbS+1n5PiOm+Dj0+zYebdvlAEnjmmpjptE5ASNL631EfMucV+Et2z1nazm+/cAIaUlMboo9rPvfXJ8x8AsSXYGTX4Y2GAZgWDCj2kJQkKMg3nSVoWZVRNmAEYhkkEVsNs6wRLYq5w6pqU6bpIIarit8Y3BJBJmJaF3lWC0Waxk7MPYi1HI8R07rfURpdQK6zSllC/otODQMS1BaHI4glGP2w3GXbpSqj64wVHEla+RPFi4k1jERo7vONFaBw9EuhVo1Fp3WxI/giX0mJbquElAa91irS8363EvmzGdHCTVbhNAju/EMEtBNUHZiBN+TEsQmry6zawjkRiRGsnDJ7NujAAAA4pJREFUZXeJ+JPjO3b1wEbLa0eOaQlCKU4pVaOU6oowW7KJEqeJ8CdnsTlPxIEc34mhlKoHtgRlHDpyTEsQSn3dQPDDZsVgVGM4skWXAfPk7QyTNdQiqcRxJcd3nJnPAjUEtbFVOXVMSxBKcRFOxO0EnbgiIRqAz/hemCd3vXObM/nI8R1fZpZhk2+/mgHH2kPChB/TMqjdJGA5kLoxuuJolqf2Y2c23G7F6DamCmgEXrM+zGeepC2YT/vLfh+7S+1nOb7jw9zPzWFmdWutSyzvm9BjWoKQEEIIx0h1nBBCCMdIEBJCCOEYCUJCCCEcI0FICCGEYyQICSGEcIwEISGEEI6RICSEEMIxEoSEEEI4RoKQEEIIx0gQEkII4RgJQkIIIRwjQUgIIYRjJAgJIQD/uDK7nd4OcXnJcHoDhBDOU0ptwBhGocrpbRGXFykJCeEApdRhpdRuc5Axx5lj9+wJN08p9YC5rQ0TvFniMiAlISGc0aK13uj0RkRDa70TQKrqRCJISUgIIYRjpCQkxCRmDtVcHWm+1nrbBG6OECEkCAkxiWmtG53eBiFGI0FIiBiYiQXbgFqttQqaVwes0FpvjXJddeZ/O4FS4EmtdbdSqgaox8hcW0cgg22jOR1zWrE5f5vWunv830qIiSNtQkLEQGu9H9gNtISZvQ04HM16lFIPAMVa60at9R6zBFNnfsYRYCtGYKrSWu83P3cf0GAut9/McDtMIDBFzRJMq5RSdUopSdUWE0JprZ3eBiFSmpk11mJtXzEv6ru11iWRlvFlxymlioEua0nKDALNvmkR3lMDHA6aVgvUa61XxPVLBm2zEPEiJSEhYleLUSqxqge2RLn8SsBWfaa1biE0oSC4iq2b0BJYZ5SfKURSkCAkRAzMEksx0GSZ1gD4qsfGzQxEQkxqEoSEiE0tcMSXCKCUqge6x5j63AQUm1VuQlxWJAgJEZt1wH6lVJVZAto31mdvzAC2E9hunZ4sXfoIkUgShISIjS9Q1Gqtt5pZa2NmBq5ms5+2DWaHoi3gr/J7HKO0VG/2du1P2zZLX/6kBKDGN02IZCfZcUI4IBUzzVJxm0Xyk5KQEEIIx0gQEkII4RgJQkIIIRwjQUgIZ1Ql06B2o/ENaoc8CCsSQBIThBBCOEZKQkIIIRwjQUgIIYRjJAgJIYRwjAQhIYQQjpEgJIQQwjEShIQQQjjm/wGoNJJidklxygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 443.077x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gs = gridspec.GridSpec(2, 1, height_ratios=[2.8, 1], hspace=0.04)\n",
    "\n",
    "# Top plot\n",
    "\n",
    "ax0 = plt.subplot(gs[0])\n",
    "\n",
    "line2 = ax0.errorbar(firas.df['freq'], firas.d, yerr=firas.df['uncert'] / 1e3, c=cols_default[0], fmt='none', capsize=3, capthick=1, elinewidth=1.5, label='COBE/FIRAS data')\n",
    "line1, = ax0.plot(firas.df['freq'], firas.B_CMB(firas.omega_FIRAS, chi2_null.x)   / (1e6 * Jy), c='k', ls='-', label='Blackbody $T_{\\mathrm{CMB}} = 2.725\\,\\mathrm{K}$')\n",
    "\n",
    "ax0.set_xlim(np.min(firas.df['freq']), np.max(firas.df['freq']))\n",
    "ax0.set_ylim(-50, 420)\n",
    "\n",
    "ax0.set_ylabel(r\"Intensity\\,[MJy\\,sr$^{-1}$]\")\n",
    "line3, = ax0.plot([],[], c=cols_default[1], ls=':', label=r'$m_{A^\\prime}=4\\times 10^{-15}$\\,eV, $\\epsilon=6\\times 10^{-6}$')\n",
    "\n",
    "ax0.set_ylim(-100, 410)\n",
    "first_legend = plt.legend(handles=[line1], loc='upper right', fontsize=14)\n",
    "\n",
    "# Add the legend manually to the current Axes.\n",
    "plt.gca().add_artist(first_legend)\n",
    "\n",
    "# Create another legend for the second line.\n",
    "plt.legend(handles=[line2, line3], loc='lower left', fontsize=14, handlelength=1.8)\n",
    "\n",
    "# Bottom plot\n",
    "\n",
    "ax1 = plt.subplot(gs[1])\n",
    "\n",
    "ax1.errorbar(firas.df['freq'], firas.d - firas.B_CMB(firas.omega_FIRAS, chi2_null.x)   / (1e6 * Jy), yerr=firas.df['uncert'] / 1e3, c=cols_default[0], fmt='none', capsize=3, capthick=1, elinewidth=1.5)\n",
    "ax1.plot(firas.df['freq'], p1 * firas.B_CMB(firas.omega_FIRAS, chi2_null.x)   / (1e6 * Jy) - firas.B_CMB(firas.omega_FIRAS, chi2_null.x)   / (1e6 * Jy), c=cols_default[1], ls=':')\n",
    "ax1.axhline(0., c='k', ls='-', lw=0.5)\n",
    "ax1.set_xlim(np.min(firas.df['freq']), np.max(firas.df['freq']))\n",
    "ax1.set_ylim(-0.15, 0.15)\n",
    "\n",
    "ax1.set_ylabel(\"Residuals\", labelpad=0.1)\n",
    "ax1.set_xlabel(r\"$\\nu\\,[\\mathrm{cm}^{-1}]$\")\n",
    "\n",
    "plt.suptitle(r\"\\bf{COBE/FIRAS CMB spectrum}\", y=0.94, fontsize=18)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig(plots_dir + \"firas_bb.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}