jason-neal/eniric

View on GitHub
docs/Notebooks/spectral_gradient.ipynb

Summary

Maintainability
Test Coverage
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing the Numerial Gradient\n",
    "\n",
    "This notebook explores the assess the differences and validity in calculating the spectral slope using finite differences or np.gradient().\n",
    "\n",
    "In *Figueira et al. 2016* the slope is calculated as `dy/dx = np.diff(flux)/ np.diff(wav)`. The `np.diff()` function shrinks the array by 1 which can be significant when slicing the wavelengths into many chunks for the telluric masking as each loses 1 pixel.\n",
    "\n",
    "The `np.gradient` function does not drop the end pixel.\n",
    "\n",
    "From the numpy documentation     \n",
    "    - \"The gradient is computed using second order accurate central differences in the interior points and either first or second order accurate one-sides (forward or backwards) differences at the boundaries. The returned gradient hence has the same shape as the input array.\"\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\jason\\miniconda3\\envs\\rtd\\lib\\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": [
    "import matplotlib\n",
    "\n",
    "matplotlib.rcParams[\"text.usetex\"] = False\n",
    "matplotlib.rcParams[\"text.latex.unicode\"] = True\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "from eniric.utilities import load_aces_spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from eniric.precision import slope, slope_grad\n",
    "\n",
    "def slope(wavelength, flux):\n",
    "    \"\"\"Original version used to calculate the slope. [Looses one value of array].\"\"\"\n",
    "    delta_flux = np.diff(flux)\n",
    "    delta_lambda = np.diff(wavelength)\n",
    "\n",
    "    return delta_flux / delta_lambda\n",
    "\n",
    "\n",
    "def slope_grad(wavelength, flux):\n",
    "    \"\"\"Slope using gradient.\"\"\"\n",
    "    return np.gradient(flux, wavelength)  # Yes they should be opposite order\n",
    "\n",
    "\n",
    "def slope_adjusted(wavelength, flux):\n",
    "    \"\"\"Slope that adjust the wave and flux values to match diff.\"\"\"\n",
    "    derivf_over_lambda = np.diff(flux) / np.diff(wavelength)\n",
    "    wav_new = (wavelength[:-1] + wavelength[1:]) / 2\n",
    "    flux_new = (flux[:-1] + flux[1:]) / 2\n",
    "    return derivf_over_lambda, wav_new, flux_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load spectrum\n",
    "wl_, flux_ = load_aces_spectrum([3900, 4.5, 0.0, 0])  # M0 star\n",
    "wl_ = wl_ * 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Here we plot a small section of the spectrum and the slopes from the gradient and dy/dx methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd809X6wPHPSZq06aZTaIGyZ7FAWQIKskTAgQIORLgi158DRcUrioobRcV59eJ1K8oQmQpeliIqUKBsiowCbbG0pXukaXJ+f6QtHWlpQ9qk7Xm/Xnm1+Sbf5Gna5sk55znnCCkliqIoiuJqNM4OQFEURVFsUQlKURRFcUkqQSmKoiguSSUoRVEUxSWpBKUoiqK4JJWgFEVRFJekEpSiKIriklSCUhRFUVySSlCKoiiKS3JzdgCOEBQUJCMiIpwdhqIoSpOxe/fuVCllcF0+R6NIUBEREcTExDg7DEVRlCZDCHG6rp+jXhOUEOJTYCxwXkrZ3cbtAngHuB7IA6ZKKffUVTwr9yayYEMcSRn5tPA3MHtUJ27qGVZn5ymKoig1V98tqM+B94Evq7h9NNCh+NIP+LD4q8Ot3JvInBUHyDeZAUjMyGfOigMA1SYbe88rOVclNkVRlJqp1wQlpfxVCBFRzV1uBL6U1iXW/xRC+Ashmkspzzk6lgUb4kqTTIl8k5lHl8by2vqj6LQadFqBTqtB76bBTWP9PvZsBsYiS6Xznll5kMSMfLz0Wjz1bhj0WjyLv/fUa/njRCoLN/5Veq5KbIqiOENDej9xtTGoMOBsmesJxccqJSghxAxgBkCrVq1q/URJGfk2j1skDGofhMlswWSRmIos1u/NkkKzpVJyKpFtLGLBhrhaxZBvMjN7+T7W7EvC31NPgJeOZl56mnlaLwFeevacvsDbm/6iwFS7xNaQ/ggVRak9e/7HL6cHyBlcLUEJG8dsblglpVwELAKIjo6u9aZWLfwNJNpIUmH+BhZMuLLK8wbO31zFeR5semwIeYVm8gqLir8Wf280M/1L20UcJrPkXGYBR85lcSGvsDQRVSffZObplQdIzTHSwt9Acz8PWvgbCPZ2R6MRDe6PUFGU2rnU/7jFIknPK+R8tpGUbGPp1/e3/GWz52jBhjiXfG9wtQSVALQscz0cSKqLJ5o9qlO5XzCAQadl9qhOdp7XGQ+dFg+dlgAvfaXzwqpJiD8+PLj0en6hmfS8Qusl18TkT3bYjCPXaOaldUfKHdNpBaG+HpzPMlJortwNuWDDUZf8I1SUpsyeltBr64/aTDSzl+9j/k9HSc0xUmSp+ef2qnqUnM3VEtRq4EEhxHdYiyMy62L8CS62JGr7h2HveTVNiAa9FoPeQAt/A1BdYvPgx5lXk5iRz7nMfJIyC0jKyOdcRj4rY23n9MSMAm798Hc6hHrTLtibDqE+dAjxprmfB9YCStU1qCj1qbqW0OjIKzidlsfJlBxOpuZyMiW39PuMPJPNxzOZJYM7BBHs406IjzvBPh6E+LoT7O1OsI87Ixf+avP9pOT9xtWI+tzyXQjxLTAECAKSgecAHYCU8qPiMvP3geuwlplPk1JecoJTdHS0bAjzoBzRZwzWxPbq+Mgqz62qG9JTr6V7mB/Hz+dwIbew9LiXXkv7EG90WkHs2cxyn7wu9VyKotivqv9VrUYgpaRsIyjEx502QV60DfbmxwNJZOYXVTovzN/A9ievrfL57Hk/qYoQYreUMrpWJ9VSfVfx3X6J2yXwQD2FU+9u6hlW6z8Ce1psVbXWXrn54h9hWo6R4+dz+Ot8DseLL7+fSKVir0C+ycy81YfoHuZHu2Cv0paWoij2MZktHEjMZOepCzaTE4DZIpk5rANtg7xoG+xFmyAvfDx0pbf3axNg1xCFvT1AzlKvLai60lBaUPXJntZamyfX2a5IKRbopadvm4DSS+crfNFqVNegopSw9X9wXfcr2Hsmg52nLrAzPo09pzNKE4ubRtgcK7pUS6iq56rP/7n6aEGpBKWUqqq7IdTXncdGdOLPU2nsPHWBhHTrfXw93OgTEYCXXsuGw8nlSvBV16DS1NjqPtMIa2myWYIQ0PkKX/q1CaBfmwCiIwLYfjzVYV1u9a3RdfEprq2qrsE5o7twU88wJvaxFlgmZuSz69QFdpxKY8epC5xMya30WK5cuqoojlZktvDi2sOVKussErzd3XjntiiiWwfg56krd3tD63KrbypBKaVq+s8S5m8grMx4WlVdg4kZ+VzILbRZdq8ojcGp1FyWxZzl+z0JpJUpPCor11jEsC6hVT6GPWPTTYVdCUoIESKlPF/hWCcpZe2WUlBcjj3/LFVNegbo/8omro+8gsn9W9O7dTNVZKE0eHmFRazbf45lMQnsjL+ARsDQTiGYzBnlqmNLuGoJd0NgbwtqmxDiGSnlUgAhxGPAPUBXh0WmNBhVdQ0+PKw9f2cZ+X53Aitjk+h8hQ939m/NzT3D8HZXjXfFdVUsQHh8ZEdaBXqyLCaBNfuSyC000ybIiyeu68QtvcIJ9fWosoT7UpV1StXsKpIQQjTHusxQARAKHAEek1LmODa8mlFFEs5XXUVRXmERq2OT+HrHaQ4mZuGl13JTzzAm929N3N/Zqv9dcSm2Eo3Auuaap17LmMjmTOzTkmgbPQLOrqyrTy5dxSeEeACYA1iA26WU2x0ZWG2oBNUwSCnZl5DJ13+eZs2+JIxFFoSAsn+CDaWCSWm8qqpm9Tfo+O3Ja1Xrv1h9JCiNPScJIf6HdSmi7lg3F1wohHjDkYEpjY8QgqiW/rwx4Up2PDUMP4MbFT8flVT/KYozWCyyyvHUzHyTSk71zK4EBXwgpZwipcyQUh4ErgIyHRiX0sj5e+rJsrFUC7juwpVK4xafmsvtH/9Z5e2q2KH+2ZWgpJQrK1wvklK+6JiQlKaiqn94L3c3isyX3nZEURzBbJH8d9tJrnvnVw4nZTEpOhyDrvxboyp2cA57u/iyhRBZxZcCIYRZCKFaUEqtzB7VCYNOW+6YViPIMRZx1yc7Sck2Oikypan4KzmbWz/6nZfWHWFguyB+fvRqXrv1Sl4d34MwfwMC67w/NS7qHHZ1qEopfcpeF0LcBPR1SERKk1HVxOAii+TpHw4w9r1t/PvOXvRuHeDkSJXGxmS2sOjXk7yz8S+83LW8c1sUN1zZorQqT02edQ0OW4tPCPGnlLK/Qx6sllQVX+NzOCmL//tmN4np+cwd04W7r4pQk3wVhziYmMkTy/dz+FwWY3o05/kbuhHk7e7ssBocl12LTwgxvsxVDRBNFVuzK4o9urbwZfWDg3hsaSzz1hxmz5kMXh0fiZeqolJqoey8pOZ+HnRr4cvmuBSaeer5aHJvrut+hbNDVKph73/7uDLfFwHxwI2XHY2ilOFn0LHormg+/OUEb/4cx5FzWXw4uTftQ7ydHZrSAFSccJuUWUBSZgF9Wjfj47uj8fdUa0S6OnvHoKY5OhBFsUWjETwwtD1RLf2Z+e1ebnz/NxZMuJLCIkuTmbGv2GfBhrhKq4uDNVGp5NQw1CpBCSHeo5quPCnlzMuOSFFsGNg+iLUzB3H/N3u4/5s95TZ6S8zIZ86KAwAqSSmlqppPp+bZNRy1LTOPAXZXc1GUOtPcz8CSGQPw0msr7UKqVqBQKgr0tt1KUhNuG47advF9I6W0Pf1fUeqB3k1DXmHlbhtQn4yVi86k5ZFnLCpd5LWEmnDbsNS2BbWz5Jvi7j5FqXdVfQJWn4wVgBxjEfd+GYPOTctT13dRE24bsNq2oMpORBnoyEAUpaZs7T/l7qZRn4wVLBbJI9/Fcjwlhy+m9WVQhyDuvbqts8NS7FTbBKXmOilOV3EFCiGsWyGoOS3KGz/HsfFIMvPGdWVQhyBnh6NcptomqM5CiP1YW1Ltir+n+LqUUvZwaHSKUoWyS9H8eiyFKZ/u5NUfj/D8jd2dHJniLKtiE/n31hPc3rcld18V4exwFAeobYLqUidRKMpluLpjMPcMasMnv51iSKcQhnYOcXZISj3bdzaDJ5bvp29EAM/f0F0ti9VI1KpIQkp5urpLXQWpKJcye1QnOl/hw+zl+0jNUaugNyXJWQXc+2UMQd7ufDi5F3o3e7e5U1yN+k0qjYKHTss7t/Ukq6CIJ5bvx1GLICuurcBkZsaXMeQYi/jv3dEEqkVfGxWVoJRGo9MVPswZ3ZnNR8/z9Y4zzg5HqWNSSv71/X72JWTy1sQoujT3dXZIioOpBKU0KlOviuCajsG8tPYwx89nOzscpQ59+MsJVsUm8diIjqqCs5GqVYISQhwQQuyv6lJXQSpKTQkhWDChB17ubsz8NhZjke1VJ5SGbePhZBZsiGNsj+Y8eG17Z4ej1JHatqDGYt1qY33x5c7iy4/AcseGpij2CfHx4PVbenD4XBZv/XzM2eEoDhb3dzYPf7eX7i38WHDrlapirxGrVZl5SaWeEGKglLLsShJPCiG2Ay84MjhFsdfwrqHc2a8V//n1JFd3DGZgezVpsyEru/GgRiPw1GlYNKU3Br3W2aEpdcjeMSgvIcSgkitCiKsAL8eEpCiOMXdMV9oGe/HY0n1k5BU6OxzFTiUbDyZm5CMBs0VSaJbsOHnB2aEpdczeBHUP8IEQIl4IcQr4N/CPmpwohLhOCBEnhDguhHjSxu2thBBbhBB7i8e2rrczRqWJM+i1vHtbT9JyjcxZcUCVnjdQtjYeNBZvWKk0bnYlKCnlbinllUAPIEpKGSWl3HOp84QQWuADYDTQFbhdCNG1wt3mAkullD2B27AmP0WxS/cwPx4b2YmfDv7Nst0Jzg5HsYPaeLDpsitBCSFChRCfAEuklJlCiK5CiHtqcGpf4LiU8qSUshD4Drixwn0kUDKhwQ9IsidGRSlx7+C29G8bwPOrDxGfmuvscJRaUturNF32dvF9DmwAWhRfPwY8UoPzwoCzZa4nFB8rax4wWQiRgLU68CFbDySEmCGEiBFCxKSkpNQ8cqXJ0WoEb02MQqsRTPl0B1fN30SbJ9cxcP5mVu5NdHZ4yiWM6VF5jpPaeLBpsDdBBUkplwIWgOJddmsy4cRWPWjFgYHbgc+llOHA9cBXQohKcUopF0kpo6WU0cHBwbWLXmlyWvgbuLlnGGcu5JOUUYAEEjPymbPigEpSLsxktrDpyHmCvPW08PNQGw82MbVdzbxErhAikOLkIoToD2TW4LwEoGWZ6+FU7sK7B7gOQEr5hxDCAwgCztsZq6IAsPFI5T+hfJOZBRvi1Judi/r6z9OcSMnlv1OiGd411NnhKPXM3hbUY8BqrHtCbQe+BGbW4LxdQAchRBshhB5rEcTqCvc5AwwDEEJ0ATwA1YenXDY12N6wpOcW8vbGvxjcIYhhXdQWKk2RXS0oKeVuIcQ1QCes3XZxUkpTDc4rEkI8iHX8Sgt8KqU8JIR4AYiRUq7Gmvw+FkLMwtpCmypVfbDiAC38DSTaSEZqsN01Ldx4jOwCE8+M7apWi2ii7EpQQogTwAIp5Udljq2VUo691LlSyh+xFj+UPfZsme8PAwMrnqcol2v2qE7MWXGg3JwaD51GDba7oGPJ2Xyz4wx39mtNx1AfZ4ejOIm9XXwmYKgQ4rPirjqoXI2nKC7lpp5hvDo+krAyLaaRXULV+JOLkVLy4trDeOm1zBrR0dnhKE5kb4LKk1JOAo4A24QQralcjacoLuemnmFsf/Ja4ueP4drOIWyJSyFN7cDrUjYdOc+2v1J5ZHhHArz0lz5BabTsTVACQEr5OvAU1jGlcEcFpSj14anrO5NnMrNwo1rx3FUUFll4+ccjtAv24q4BrZ0djuJk9iaosmNGm4BRwPsOiUhR6kn7EB8m92vF4h1nOJasNjd0BV/8Hs+p1Fzmju2KTqv2U23qarthYefibxOFEL1KLkAgsNbh0SlKHXt4eEe83N145ccjzg6lyUvNMfLupr8Y0imYoZ1UWblS+yq+x4B7gTdt3CaBay87IkWpRwFeemZe24GXfzzCL8dSuKajWpXEWd78+Rj5JjNzx1RcP1ppqmq7YeG9xV+H1k04ilL/plzVmq93nObldYcZ2G4wbqprqd4dTspiya4z3H1VBO1DvJ0djuIiapWghBDjq7tdSrni8sJRlPrn7qZlzujO3Pf1HpbEnOXOfmpwvj5JKXlh7SH8DDoeGabKypWLatvFN66a2ySgEpTSII3qdgV9IwJ46+dj3HBlC3w8dM4OqcnYcOhv/jx5gRdv7Iafp3rdlYtq28U3ra4CURRnEkIwd2wXbnh/Ox9sOcGToztf+iTlshmLzLz84xE6hnpze99Wzg5HcTH2rmaOEGIM0A3rYq4ASClfcERQiuIMPcL9Gd8rjE9/O8Wd/VrRMsDT2SE1ep/+Fs/ZC/l8dU9fNfanVGLvjrofAZOwbiYogAmA6rhXGrzZozqh0cBr6486O5RG73x2Ae9v/ovhXUIY3EFVTyqV2fuR5Sop5RQgXUr5PDCA8vs8KUqD1NzPwIyr27F2/zl2n77g7HAatTc2xFFotvC0KitXqmBvgirZsyBPCNEC6+KxbRwTkqI4133XtCXEx50X1h7BYlFLTDrSyr2JDJy/mTZPrmNpTAKD2gfSJsjL2WEpLsreBLVWCOEPLAD2APHAd44KSlGcyVPvxuxRndh3NoM1+ytu+KzYa+XeROasOEBiRn7pytJ/nLzAyr2JTo1LcV12JSgp5YtSygwp5fdYx546SymfcWxoiuI8t/QKp1sLX15fH0dBmf2jFPst2BBXbi8ugAKThQUb4pwUkeLq7C2S0AohbhBCzAQeAO4RQjzq2NAUxXk0GsHcMV1JzMjnk99OOTucRiHJxm7G1R1XFHu7+NYAU7EuEutT5qIojcaAdoGM7BrKv7cc53x2gbPDafBa+HtUcdxg87ii2DsPKlxK2cOhkSiKC5pzfReufWMrQxZsJb/QTAt/A7NHdVK78NphQNtAlu8pP95k0GmZPaqTkyJSXJ29LaifhBAjHRqJorigfWcz0GgEeYVmJJCYkc+cFQfUwH4tnUjJYe2Bc3S6wocwfw8EEOZv4NXxkSrZK1WytwX1J/CDEEKDtcRcAFJK6euwyBTFBSzYEIe5Qql5vsnMgg1x6o21hkxmC48uicVDp+XLf/Ql1Nd2V5+iVGRvgnoT6+TcA1JKNVFEabTUwP7l+2DLcfYlZPLBHb1UclJqxd4uvr+Agyo5KY1dVQP4amC/ZmLPZvDe5uPc3DOMMT2aOzscpYGxtwV1DtgqhPgJMJYclFK+5ZCoFMVFzB7ViTkrDlSav/PPa9TCKZeSX2jm0SWxhPi4M++Gbs4OR2mA7G1BnQI2AXpUmbnSiN3UM4xXx0cS5m9AACE+7rhpYO3+vzGZLc4Oz6W9+tMRTqbm8uaEK/EzqH2elNqrdQtKCKEFvKWUs+sgHkVxOTf1DCtXELFybyKPLIll/k9HeWasWujUlq1x5/nyj9PcM6gNV7UPcnY4SgNV6xaUlNIM9KqDWBSlQbipZxhTr4rgk99OsWafWquvovTcQp5Yvp8OId5qjpNyWewdg4oVQqwGlgG5JQellGrLd6VJeOr6LhxIzORf3++n0xU+dAxVPdwAUkrmrjxIel4hn07tg4dO6+yQlAbM3jGoACANuBYYV3wZ66igFMXV6d00/PvOXnjq3fjnV7vJKjA5OySXsCo2iXUHzvHI8I50D/NzdjhKA2dXC0pKOc3RgShKQxPq68EHd/Tkjv/u4LGl+/jP5N5oNMLZYTlNUkY+z6w6SO/WzbjvmnbODkdpBOxdzTxcCPGDEOK8ECJZCPG9ECLc0cEpiqvr1zaQp67vwv8OJ/PhLyecHY7TWCySx5ftw2yRvDXxSrRNOFErjmPvGNRnwGJgQvH1ycXHRjgiKEcwmUwkJCRQUKBWob5cHh4ehIeHo9OpUmFb/jEwgr1n0nnz5zh6hPsxuEOws0Oqd5/9Hs/vJ9KYPz6S1oFqh1zFMexNUMFSys/KXP9cCPGIIwJylISEBHx8fIiIiEAI9WnOXlJK0tLSSEhIoE0bNTnVFiEEr93Sg2PJ2cz8di9rHhpEeDNPZ4dVb44lZ/Pa+qMM7xLCpD4tnR2O0ojYWySRKoSYXLxxoVYIMRlr0cQlCSGuE0LECSGOCyGerOI+E4UQh4UQh4QQi+0JsKCggMDAQJWcLpMQgsDAQNUSvQQvdzc+mtybIrPk/m/2NPpdeFfuTWTg/M20eXIdY97dhk4jeHV8D/X/pjiUvQnqH8BE4G+syx7dWnysWsWTfD8ARgNdgduFEF0r3KcDMAcYKKXsBtjdMlP/LI6hXseaaRvszZsTr2R/QibzVh9ydjh1ZuXeROasOEBiRj4SMJklJrOF7cdTnR2a0sjYlaCklGeklDdIKYOllCFSypuklKdrcGpf4LiU8qSUshD4Drixwn3uBT6QUqYXP9d5e2JUFGcY2e0KHhjaju92neW7nWecHU6dWLAhrtLahIVmyYINcU6KSGmsajUGJYR4tpqbpZTyxUs8RBhwtsz1BKBfhft0LH6u7YAWmCelXG8jlhnADIBWrVpd4mkvbeXeRBZsiCMpI9+hu6a+/PLLLF68GK1Wi0aj4T//+Q/9+lX8ke3zyiuv8NRTTznksRTHeXREJ/YnZPL0Dwd483/HSM02NqqdeNUWJEp9qW0LKtfGBeAe4F81ON9WX1HFLTvcgA7AEOB24L9CCP9KJ0m5SEoZLaWMDg6+vKqpil0Wjto19Y8//mDt2rXs2bOH/fv3s3HjRlq2dNwg8iuvvGLzuJQSi0UtZOosWo1gVLcrsEhIyTY2qp14d59Or3Kul9qCRHG0WrWgpJRvlnwvhPABHgamYe2qe7Oq88pIAMq+Q4cDFRczSwD+lFKagFNCiDisCWtXbWIt6/k1hziclFXl7XvPZFBYYWXqfJOZJ5bv59squmm6tvDluXHVbyFw7tw5goKCcHd3ByAoyLpoZkREBJMmTWLLli0ALF68mPbt25OSksJ9993HmTPW53z77bcZOHAgOTk5PPTQQ8TExCCE4LnnnmPXrl3k5+cTFRVFt27dePnllxk9ejRDhw7ljz/+YOXKlXTr1o2cnBwAli9fztq1a/n888+ZOnUqBoOBo0ePcvr0aT777DO++OIL/vjjD/r168fnn39+6RdVqdaHW09U+uTVkHfiNVskH/1ygrf+dww/gxu5RjPGoov/MwadVq27pzhcrceghBABQoiXgP1YE1wvKeW/ajhWtAvoIIRoI4TQA7cBqyvcZyUwtPi5grB2+Z2sbZy1UTE5Xep4TY0cOZKzZ8/SsWNH7r//fn755ZfS23x9fdm5cycPPvggjzxirQN5+OGHmTVrFrt27eL7779n+vTpALz44ov4+flx4MAB9u/fz7XXXsv8+fMxGAzExsbyzTffABAXF8eUKVPYu3cvrVu3rja29PR0Nm/ezMKFCxk3bhyzZs3i0KFDHDhwgNjY2Mv6uZXG1Q2WnFXAXZ/sYMGGOK6PbM7W2UN57ZYepVuQhPkbeHV8ZINMvIprq+0Y1AJgPLAIiJRS5tTmfCllkRDiQWAD1vGlT6WUh4QQLwAxUsrVxbeNFEIcBszAbClljUrYq3Kpls7A+ZtJtPHGEeZvYMk/B9j9vN7e3uzevZtt27axZcsWJk2axPz58wG4/fbbS7/OmjULgI0bN3L48OHS87OyssjOzmbjxo189913pcebNWtm8/lat25N//79axTbuHHjEEIQGRlJaGgokZGRAHTr1o34+HiioqJq/wMrpVr4G2z+TWmE4McD5xjd/YoGUR256Ugyjy/bR4HJwuu39mBC73CEEJW2IFGUulDbibqPYd1Bdy7wdJl/MIG1SML3Ug8gpfwR+LHCsWfLfC+BR4sv9cLWrqmO6rLQarUMGTKEIUOGEBkZyRdffAGUL90u+d5isfDHH39gMJTvy5dS1ujNzMur/Az+sudUnMdU0u2o0WhKvy+5XlRUVJMfTamGrb8pvZuGQE8d93+zh6vaBTLvhm4uuwq6scjM/J+O8tn2eLo09+W923vSPsTb2WEpTUytuviklBoppUFK6SOl9C1z8alJcnJVFXdNdVSXRVxcHH/99Vfp9djY2NKutyVLlpR+HTDA2kobOXIk77//frn72zqenp4OgE6nw2SqehXt0NBQjhw5gsVi4Ycffrisn0WpHVt/U6/f0oNt/7qWF27sxqGkLEa/s40X1hx2uZXQT6TkcPMHv/PZ9nimDYzgh/uvUslJcQp7lzpqdOqiy6KkuCEjIwM3Nzfat2/PokWLWLt2LUajkX79+mGxWPj2228BePfdd3nggQfo0aMHRUVFXH311Xz00UfMnTuXBx54gO7du6PVannuuecYP348M2bMoEePHvTq1YuXX3650vPPnz+fsWPH0rJlS7p3715aMKHUj6r+pqYMiGBsjxYs2BDHZ7+fYvW+RJ64rjO39gp3ymroZadY+HnqyC0owtvDjU/ujmZYl9B6j0dRSghrj1rDFh0dLWNiYsodO3LkCF26dHFSRNWLiIggJiamtKqvIXDl17MhO5CQyXOrD7LnTAZRLf15/oZunErNrZM5ebaUTLEo2xWpEfDsuK5MvUqtvahUTQixW0oZXZfPoVpQiuJEkeF+LL/vKn7Ym8irPx3lxg+2o9UIzBbrB8eS+VOAQ5OUlJL4tDzmrTlUaVUIi4SPfz2lEpTidCpBOUF8fLyzQ1BciEYjuKV3OCO7hTLg1U3kGMsnjHyTmVd/OsK4K1tUu89SdauhGIvMHEzMJCY+nd2n09lzJp3UnMIqH6shlsMrjU+jTlA1rX5TqtcYuoEbAh8PHblG26ugJ2cZ6fLMeloGGIgI9KJVoCcRgV60Lv66+/QF5q682BpKzMjnieX7WL0vkaz8IvYnZlJYPLE2ItCTazqG0Lt1M97eeIzz2cZKz6dWhVBcQaNNUB4eHqSlpaktNy5TyX5QHh4ezg6lSahq/pS/Qcekvi05nZpHfFouv59Iq9Q1V1GhWbL5aAo9W/lz94DW9G4dQO/WzQj2uTitwFOvrbMpFopyuRptggoPDychIYGUlBRnh9Lgleyoq9S9qubkzbuhW7kxKCklKdlFrxlpAAAgAElEQVRG4tPyOJ2Wy+zl+20+ngB+uH9glc9X8pj1VZShKLXRaBOUTqdTO8AqDU5NE4YQghBfD0J8PejbJoC3N/5ls+VVk646tSqE4qoabYJSlIbKnoRRl6uhKIqzqASlKI2A6qpTGiOVoBSlkVBddUpj0yhWkhBCpACngSAg1cnhuBr1mtimXpfK1GtSmXpNKit5TVpLKS9vt9hLaBQJqoQQIqaul95oaNRrYpt6XSpTr0ll6jWprD5fk1pvWKgoiqIo9UElKEVRFMUlNbYEtcjZAbgg9ZrYpl6XytRrUpl6TSqrt9ekUY1BKYqiKI1HY2tBKYqiKI2ESlCKoiiKS1IJSlEURXFJKkEpiqIoLkklKEVRFMUlqQSlKIqiuCSVoBRFURSXpBKUoiiK4pJUglIURVFckkpQiqIoiktSCUpRFEVxSSpBKYqiKC6pUWz5HhQUJCMiIpwdhqIoSpOxe/fu1LreUbdRJKiIiAhiYmKcHYaiKEqTIYQ4XdfPobr4FEVRFJfUKFpQiqK4vuiX/kdqTmGl40HeemLmjoD9S2HTC5CZAH7hMOxZ6DHRCZEqrkIlKEVR6oWt5FR6fP9SWDMTTPnWg5lnrddBJakmrNEmKJPJREJCAgUFBc4Opcny8PAgPDwcnU7n7FAUV7fphYvJqYQp33pcJagmq9EmqISEBHx8fIiIiEAI4exwmhwpJWlpaSQkJNCmTRtnh6O4NGltMdmSmVC/oSgupdEWSRQUFBAYGKiSk5MIIQgMDFQtWAUpZbnrwi0LQ6v/ILTZhIvzfKpbUPXJfuF1HJ3iyhptggJUcnIy9fo3bWaL5Ivf47n5379jLDKXHtcHbkLrGU+vkI/ZqJ9NX81RiJwAOkP5B9AZrIUSSpPVaLv4FEVxnv0JGTz9w0EOJGYyuEMQ2QVF+HSeC6Ko9D7H/M/Tx785WNyY6jObWWNHotlcoYovrDf8+gYMfgzUB54mRyUoalD+aietVktkZGTp9ZUrVxIfH8+NN95YOi4TFBTExo0bmTdvHh9//DHBwcHk5uYSGRnJSy+9RNeuXe1+fkWpb/mFZl796Qhf/XmaYG933ru9J2N7NEcIweYblvHcj9PYVpRemmyuanEVUR73sWDTFmK75jP37u9p798erUZrfcAtr8Avr4G0kNJnGrN/nc0b17xBkCHIiT+lUl9UguIS5a+XwWAwEBsbW+5YfHw8gwcPZu3atZXuP2vWLB5//HEAlixZwrXXXsuBAwcIDq7T1UQUpVaq+0C346nh+B1fyV7vL/EznUdsDgdhnc+UbTGxoygDhECn0VFkKSLMK4z7BkTxy98ric1bwa1rPsVb582VIVfSM7gn0wY9ij79NGx5mY8yYtmTcZAP933IM/2fccJPrtS3JpOgJv3nj0rHxvZozl0DIi557oXcQv7v693lji355wBHhWbTpEmTWLduHYsXL+bhhx+u0+dSlNqo7gOd9uAyHi14H1F2PtMP90FRITFeeqRWx+hWw7kn8h6WHVtGan4qQggW3/IcM5cPYv3xP+nRJYvk3ON8m/YtM3rMoHdeDIVtWkHGAQCWxi1ladxS9Fo9uyfvthmL0jg0mQTlDPn5+URFRQHQpk0bfvjhBwC2bdtWenzChAk8/fTTNs/v1asXR48erZ9gFcURNr1wMTkVk9KM2PwCEx8/xqiIUfi5+wEwt//c0vtoNBrevmUYDyxuRnx8HqsfGojEhBCC9ePXc9e6ySTmJoEAvdAyIuI6Hu/zeL3+aEr9azIJ6nJaPAFeervOt9XFB1TZxVdRxfJcRXF5FeYtJWu1PBQazJy0dHpCaXKyxU2r4d3be5JfaMbdTQtYx6GCPYO5Knwgy48tBykplEUUZiepcagmoFGXmTd0e/fupUuXLs4OQ1Fqrsy8pVSNhulXhHBa54bWK6RGp7u7afH31GMsMjNrSSx/nkwD4ELBBSZ2msiiq9/AC8Fv53dTlLi3Tn4ExXU0mRZUdYK89VUO+jrL999/z88//8ybb77ptBgUpdYGPQrrZpGu0XBv8xCS3bR8mJpFj1G1+zvOM5o5kJjJ/w4n892M/rw99O3S277T+5L3/T24HVsPYT0d/RMoLsRpCUoI0RL4ErgCsACLpJTvCCECgCVABBAPTJRSptdlLJdTSu5ICxcu5OuvvyY3N5fu3buzefNmVcGnuJwqP9B56eDkFrKFln+2aM5ZLXyQI+g96s1ar6fXzEvPV/f05dYP/+DuT3ey7L4BtA32BiAifABM/xU8A1lzYg1DWw7FW+/tkJ9NcS3CWeMcQojmQHMp5R4hhA+wG7gJmApckFLOF0I8CTSTUv6ruseKjo6WFTcsPHLkiOoecwHq99C4/HoshTd/juOju3rT3K/Cyg/7l8KKezENe5ZnSeX6NtczOHzwZT3fyZQchr35C7bepQL8M5AtFtDbJPl3pgl9VpLapqMeCSF2Symj6/I5nNaCklKeA84Vf58thDgChAE3AkOK7/YFsBWoNkEpilL3MvNMzF6+D18PHc08K3R/ZyaS9+PjFLSMJmDgI7xaMtH2MrUN9raZnAAuZPjzXtfePJW+k6fcC3gdiUZt09GouESRhBAiAugJ7ABCi5NXSRKzOboqhJghhIgRQsSkpKTUV6iK0mQ9u/ogaTmFvDUxCg9dmQRksVCw8j5mNjNwTzMDJiz1FtO4Ezt49EI6G7y9eD2gmTWZlWzToTR4Tk9QQghv4HvgESllVk3Pk1IuklJGSymj1TiNotSttfuTWBWbxMxhHYgMt5aKp+SlMHX9VJK2L+TRvCPs9NDzjyv/iU5Tj/t/ZSYwNTObyZlZfOPnw2G9vvS40vA5tYpPCKHDmpy+kVKuKD6cLIRoLqU8VzxOdd55ESqKIqV1VfIrW/pz/5B2pcc/2v8Re5L3MM2ymyRPA8/1f5Zx7cbVb3B+4YjMs8y+kMGI3Hy6FRaWHlcaPmdW8QngE+CIlPKtMjetBu4G5hd/XeWE8BRFKSaE4Kt7+pGRZ8JNq6H3170pNF+s4ksq7od5ddd8bu00oX6DG/YsrJmJxpRPL6MRgD893DG16sHllWcorsCZXXwDgbuAa4UQscWX67EmphFCiL+AEcXXFUVxgh0n08gxFuGh03KFnwcA68ev5/o21+NR3JWn0+gY02YMG27ZUCcxVDcfcZVlIIx7F/xaAgLpF877Ic15LGsf+/d+UifxKPXHaQlKSvmblFJIKXtIKaOKLz9KKdOklMOklB2Kv16ol4D2L4WF3WGev/Xr/qWX/ZBarZaoqKjSS3x8PFu3bsXPz6/02PDhwwGYN28eYWFhREVF0aFDB8aPH8/hw4fLPd6tt97KyZMnq3y+IUOGULHcHiAmJoaZM62VTUajkeHDhxMVFcWSJUt4++23ycvLq9XPtXXrVsaOHQvA2rVree6552p1vtIwnE7LZdrnu5i3+lC548GewRgLMjCaC9GjochShJfeq86WHoqZO4L4+WPKXeJeuo5B7YPILzRbq/VmHYR5GYhZh3h7/CqCcOOBvW9x6vhPdRKTUj+cXiThEvYvtZamZp4FpPXrmpmXnaRK1uIruURERADWtfhKjm3cuLH0/rNmzSI2Npa//vqLSZMmce2111JSoXjo0CHMZjNt27atdRzR0dG8++67gHX5JJPJRGxsLJMmTbIrQZU1ZswYVq9efVmPobges0Xy2NJ9aDWCR0d0LHfbwdSDbD73Bx29w1k84mMmdppIWn5avcbn7qbly3/05ba+rQAoMl+sHAzyj+A/oz5FIwT3/TqbI/Fbmbp+Kqn5qfUao3L5mk6C+mxM1ZdVD1pLU8sy5cNPxdOvctMqn1PHJk2axMiRI1m8eDEA33zzDTfeeCMAZrOZqVOn0r17dyIjI1m4cGHpecuWLaNv37507NiRbdu2ARdbPOfPn2fy5MnExsYSFRXFO++8Q1JSEkOHDmXo0KEA/PzzzwwYMIBevXoxYcIEcnJyAFi/fj2dO3dm0KBBrFixovT5hBAMGTKkRovfKg3Hol9PEnM6nRdu7EYL/4sTclPzU3l400M092rOorGL6dSiL3P7zy23FFF90Wismx5uOXqekW//SnJWQeltLVv05t+DXiNDwPObZrIneQ8f7vuw3mNULk/TSVDVMRttH8+/vN7Fku02oqKiuPnmm0uPl2y3ERUVxcsvv1zl+WW329i+fTu9e/cGIDY2lsTERA4ePMiBAweYNm1a6TlFRUXs3LmTt99+m+eff77c44WEhPDf//63tAX38MMP06JFC7Zs2cKWLVtITU3lpZdeYuPGjezZs4fo6GjeeustCgoKuPfee1mzZg3btm3j77//Lve40dHRpclQafgOJ2Xx1v/iGN39Cm6KCis9XmguZNb6e8jKS+GddrcR4BHgxCgvCvZx5+/MAqZ9tosc48Ut5afseJY8jeCQm0QiWRq3lMgvIun9dW8nRqvURtNZLHbauqpvW9i9uHuvAr+W1q9egdWfXwVHbrdx7ty50nX52rZty8mTJ3nooYcYM2YMI0eOLL3f+PHjAejduzfx8fG1ivfPP//k8OHDDBw4EIDCwkIGDBjA0aNHadOmDR06dABg8uTJLFq0qPS8kJAQkpKSavVciuvy8XBjRNdQXropElG8NTvAG3++QmzWSRYUaOnczXVWaege5se/7+zFPV/E8H9f7+bTqX3QaTWsH7+eN2LeYPOZzRSYCxAIBoUN4oWBahJvQ6FaUGAtVdVVWFdMZ7Aed6Ky220YDAYKCqxdGM2aNWPfvn0MGTKEDz74gOnTp5ee4+7uDlgLNIqKiio/aDWklIwYMaJ0fOzw4cN88om1EqrsG1VFBQUFGAyGKm9XGpaWAZ78+87eBHiVr567+e9TPH4hg+vGfgx6LydFZ9uQTiG8cnN3tv2VytM/HEBKSbBnMF46L4xmIzqNDonk98TfOJx2+NIPqLiEptOCqk7Jml2bXrDOQHeBBScrbrfRpUsXjh8/TkREBKmpqej1em655RbatWvH1KlT7X4eHx8fsrOzCQoKon///jzwwAMcP36c9u3bk5eXR0JCAp07d+bUqVOcOHGCdu3a8e2335Z7jGPHjtG9e/fL+XEVJ4l+6X9VbjUTM3cE7F9K2uYXCMw4SxegS6froWXf+g+0Bib1aUViej7peSYsErTi4j5SEzrcypc/P8RmUxoPbX6Ix6MfZ3KXydV+8FKcTyWoEj0mOn1xyeq22xgzZgxbt25l+PDhJCYmMm3aNCwWa+XSq6++avdzzpgxg9GjR9O8eXO2bNnC559/zu23346xeNLjSy+9RMeOHVm0aBFjxowhKCiIQYMGcfDgwdLH2LJly2XFoDiPreRUenz/Us79+Ai3hfhzl/RlemYWnNxirW510YVYZxVXHAohKCyylCveePm2n8kryuep357i9V2v4651Z2In1/w5FCunbbfhSE1hu438/HyGDh3K9u3b0Wods1K0IyQnJ3PHHXewadMmm7c3tt9DYxPxZNVjq0dC/8XdXibO6txYnPQ3bUzFXcZ+La3zjlzYqdRcpn62k5du6s7gDuXX6rQk7ubrVVO5ZcyHeLUe5KQIGz6X2W5DCDFBSrnsUseUumMwGHj++edJTEykVatWzg6n1JkzZ9Suv42S5Dl9Hkf1nryfnHIxOUGDWIg1yFtPQno+d32ys9Jtnbxy2eADfHMb+e7ezPMo5MEiT1oOVftIuZqaFknMqeExpQ6NGjXKpZITQJ8+fYiKinJ2GIqD6QN+5SdvL2amZ3J1fkH5GxvAQqw+HjrMFtu9Q3G5XtBnOhRmc8aYynaDB3f4WNiz4VGHrCCjOE61CUoIMVoI8R4QJoR4t8zlc6B2JWKKojQYFlMzbgq4kntyK4xRuUB1q0Ps/BiAToUmvklKxt9sYXqwP2u2qRJ0V3KpFlQSEAMUYN2SveSyGhhVt6EpilJXvtt5hie/329jIVYz92tX8XRRHC+O/Qpxw3ulC7Hi19K6MGsj6AaTZbopWxcV8fW5ZHoWGHnKW7D4iHX1lpL9rtQSSc5T7RiUlHIfsE8IsVhKaaqnmBRFqUNLY84y54cDXN0hmO1PXou7m7XoJqcwh6nrp9LF5MY4vEFKl6hurQspmmBCLBe3mvOzWPjo7/MsDPBncNwv0OWO0v2uPtz3Ic/0f8aJ0TZdNS0z7yuEmAe0Lj5HAFJKWfuVSxVFcZoVexL41/f7GdQ+iP/c1bs0OSXnJnPL6vFkF+ZwxciPIbQ3aBrvPP73uJ058kM8xcUuTJ3OwBNhI+mdtZPCLyJLjy+NW8rSuKXotXp2T97tjHCbrJr+BX4CvAUMAvoA0cVfGxVHN+mFEDz22GOl19944w3mzZvnkMdWlNpaFZvI48v2MaBtIB9PicZDd3G6wkNrbiezMIuogC70ad4XNK4zlcFeVe0jFeStZ86Tz/JLp2dIlEFYpMDkHWbtvrzlY9bf+jNXh11den8NGvqE9uGn8WrrjvpW0xZUppSy0f92HN2kd3d3Z8WKFcyZM4egoLrZK0dRaqqZp55rOgbz7zt7lyanirvj7rlwiMgvIhtFayFm7ohqbx99x0zOZ83g472JzLi6LQjBrvgL/N/XsWR7GdE3s3YVmbGwK3kX1y++iWXjv6RNwh6XWnWmMatpgtoihFgArABKl/6WUu6pk6jqwLT10yodGxUxits630bvr3pTaLn4T1rSpHcTbuydspf0gnQe3fpouXM/u+6zSz6nm5sbM2bMYOHChZVWLZ86dSpjx47l1ltvBcDb25ucnBy2bt3Kc889R2hoKLGxsYwfP57IyEjeeecd8vPzWblyZenyRh4eHhw6dIjk5GTeeustxo4dy+DBg3nvvfdKS78HDhzIhx9+SI8ePWr9mimNQ1JGPi38DVzdMZjBHYLKLe+zvt3d3HfgA0646zEDHloPhrUaxuN9HndewPUoxNeDf17TDoD03EKmfLKTfJMZD78cCtP7Y8roi87/TzwMp+lviSfs0zGQf4Hf9BqC9G50Ltk7DlSSqgM17eLrh7Vb7xXgzeLLG3UVVH1bedNKAtwDEFj/cQWCAI8AHo1+9BJnXtoDDzzAN998Q2ZmZo3P2bdvH++88w4HDhzgq6++4tixY+zcuZPp06fz3nvvld4vPj6eX375hXXr1nHfffdRUFDA9OnT+fzzzwHrGnlGo1ElpyZsw6G/GbJgKz8fsm6RUpKc8kx5xG1/g+ANz3ClewgWBHqtHqPZWKe747qyZl56PrnbujBCQeJdGJNvwmJsgTF5PJnxj1CQMAV9XirSXMiCgGZMCGvOlOYh/KQXmDZdLE9X1X+OU6MWlJRyaF0HUteqa/G09GnJsNbDWH5sOXqtHpPZxPDWw7mr610ANPNoVqMWky2+vr5MmTKFd999t8Yrfvfp04fmzZsD0K5du9LtNCIjI9myZUvp/SZOnIhGo6FDhw60bduWo0ePMmHCBF588UUWLFjAp59+elkLySoN28bDyTy4eA/dw/wY0C6w9HhqfioPrrmdczmJ/NjuWi60CGeiVygTOk5g2bFlTfqN9ar2VSVmwUZLb6TFjAC+PPc3K729WeLrzRMhQQQVmXls9weM7XW/qv5zoJoudRSKtfXUQko5WgjRFRggpfykTqOrR6WrHtfBP+kjjzxCr169ym0s6ObmVrrYq5SSwsKLXYwlW2YAaDSa0usajabcFhoVV2IWQuDp6cmIESNYtWoVS5cupeIahUrTsPloMv/3zW66Nvfli3/0xcdDB0B8Zjz/9+NdpBZcYIE2DK/bFvN2ma1m5vaf66yQG4RESyDhmlT8LJK7s7K5Kyub3w0efOvrwzMHPmTOwY9K71syVIB048DUvdZVKtTYVa3UdAzqc+Az4Oni68eAJVir+xqFsqseO/qfNCAggIkTJ/LJJ5/wj3/8A4CIiAh2797NxIkTWbVqFSZT7aeZLVu2jLvvvptTp05x8uRJOnXqBMD06dMZN24cgwcPJiDANXY9VepWVdtmJKTn4+uhg/1L2bf1BR70tiCATws8iZy+qvI+aEq1kvs+QcjuZ9BL61C8BhhUJBjU+0luWZPIoeDjuPnuR4jiD58WN4py2/PNun9y875VeBbmWx+owthVSl4Ks3+dzRvXvNEku1erUtMxqCAp5VLAAiClLALMdRZVI/TYY4+RmnqxVXbvvffyyy+/0LdvX3bs2IGXV+03gOvUqRPXXHMNo0eP5qOPPsLDwwOw7qbr6+tbrsWmNG5VbZuRlmvdNoM1M1miycPHYuHrpGQiLyTA0Uvv6twUVVee3nvsP9Hf/D7SryUSQYFXCxj3LsZe09htikaa3QEJFo31i8kPnUcSb6ZsR2uyJqcfvL1Y5OfLLq2FguKxq7LdgspFNdpuQwixFbgF+J+UspcQoj/wmpTymjqOr0aawnYbFVWsAiwrKSmJIUOGcPToUTROnmzZ2H8PrqK6bTMOhz6BZ2YCRgG5QkNAcddyQ9g2o6EwFpnpNHc9HmFfIYt8iqv/diLcstEkTuBXz+kEFb/uc4IDWetd/IFUSrC1aWID6BZ0me02gEexrr/XTgixHQgGKr8zKk735Zdf8vTTT/PWW285PTkp9SOroKruYQv64J+5w9vMV1kCHylxl5aLNzeAbTMaipIVOQoS7yo9Zky+qfR7H0NzyE0E4NWUNJ5MS2evuztbDZ787OVJrhYsQiAtOjQSAs0mHlo2htbJR2ltyae73o0uZboFo1c3IzWnEOGWhUeLbylIvANp9im3E3J1ia2hdCnWtIpvjxDiGqAT1rlrcWptPucqKSWvaMqUKUyZMqV+g1Gc5mBiJg8stjUdsYgrWnxOrt9xehW4YbDVU9IAts1oLNxHzUOunokosnbz+Vks9M+3sCJzEmnBKeia7QSLFkQR7sZmtDGnczYjnj+8PDD6eHJLVg7z0i5gMeVz15/P0to/iLae/iS655HlmYJvyCoyz91Jao4RuW8JYu3DYLI93gWOX5SgrlSboIQQ10opNwshxle4qaMQAinlijqM7bJJKStVuin1pzHs1uzKDiZmMv7D3wnwvDhmItyy8Aj7muZcIMUzh+6pETwzfHL5NyxoPNtmuJAgb73NscAgbz30GGOdZVmmVeMx7FlWL/bCw+0rTOn9SrsFi9yy2Zz4BCfd7wQh+dvt4rJTORqBt7mI/d4ZIC7OrbT4H8THfw7SoiVqbxFezQPxlha8LdbL3ZnZDNv0Ar33vWZzUQJXXTnkUi2oa4DNwDgbt0msK0vUCSHEdcA7gBb4r5Ryfm3O9/DwIC0tjcDAQJWknEBKSVpaWmnhhuI4JR+8ujb35f4h7ZgyIIKRC38hNacQn8D1SMMZUpHcmhzAJuO9iCtvsI5zuOhYRmNxqaWVbK4Mv3hdld2Cwj8ckXmWFkUX69F8LZJX0zT0LngK75CV4HMMNGY0Fg3G7B4Upg5lVsA8cjQacjSCHI2GXI3GWg2XmcA7E9dy/8b7kVg/PLr6yiGX2m7jueKv9VoOJoTQAh8AI4AEYJcQYrWU8nBNHyM8PJyEhARSUlJs3n4uMx+zBYSwILS5SLMXUmrQaqC5nwEKc6EgEyxFoHEDDz/Qe9l9XlPk4eFBeLjqRnKkw0lZPLPqIO/f0ZPmfgYeGd4RABnxJD7mQkrarBLB8tB09NrngRsa7bYZjZkY9qy1a65Cyzdg3EvIxV4Umn3RCQtY3DALM9LsjqUwlIekL1w4W/kB/VoyKGwQt3a8leXHlqPT6lx+5ZBLdfFVu9aPlPItx4ZTqi9wXEp5sjiO74AbgRonKJ1OR5s2baq8fXRx1ZN76A/omu3ElNkLY/I4kFpO3FaAdu0jlbtExr3L6JVe5c9L74sx+WYA4u/IhbWV/6AY9+5lD2oqTZuUku92nWXe6kP4GXSkZButH4gALBaebtaXj85tJd3NjQLh+p+MFatLdQsCtt8XFq9DuOWU6xoUbtnW+1eR2Eq6dOtyUQJHq7bMXAjxXPG3nbBur7G6+Po44Fcp5fQ6CUqIW4HrSh5fCHEX0E9K+WCZ+8wAZgC0atWq9+nTp2v1HN0/64nQVL1r/eicXF5PSQNgZMsWFCLQIflbq7VZFiql4M58M77GbLwtFnoYjfQ0FmIG9vuH8mTGneSY/cgKjEH676+c2Gz9QV0qsd2QrpJaE5BrLGLuyoP8sDeRwR2CWDgpiiBv6+oiMjeNT1ZM4B1LCgFoSceCTqvDZDYxodMElx4AV+xX1cTs+vzAWx9l5jWdB/UzcIuUMrv4ug+wTEp5XZ0EJcQEYFSFBNVXSvmQrfvbmgd1KW3mfot7yI+4+RxCaExIqcVcEIo5txOPalfSrrCQEXnWhPFagD8FQmASgsOEkuxRQK6bGbMAnQWQWowWAz7aLHKLS7unZ2TycHommRoNg1pX3c2ls0CfwiL8i4z4my34WSz4Wcz0yzfS3hDCTanTSHMTZATElktsN2h+412vz2wmtZI/xIZSSqpcVNUbD8CjIzryYNBeNJutbzxGv3Dm+XuyVpPPaN+O5Pu2ILTCmnplV0hRFEdypXlQrYCy/zWFQITDo7koAWhZ5no4kOTIJ5BFvtZZ36IIaXEDYcaS35LClFHcF7oZ8i724f7rQob1G7+WRCQ/U9q9h0VLoTBjyojGmHwzMQGP456XQJ4oWRcdDBYLbyXn8a5lMOl+x8kyZCGFRGsRdM01MDzDjf8Fp3DGXU+mRku21prg5qWk0T4riTMdv8AkSlZZB33ADvQBO9gmJcTnk+imZYOXJ4FmCwFmM4Gb5nLvWiOJWX64h65G12wPAz96EmPyzao7sQGoKjkBzAzeW66l/S/3fDZpBA/6X8mMG74qVwyk1tRTGoOaJqivgJ1CiB+wVu/dDHxZZ1HBLqCDEKINkAjcBtzh6Cexqw93cdXnGa57HtbMxLvMeXqdgRHXvcq9i71wd/8BnWEnWNwoEgBt+rsAABq7SURBVGZ2myL5PedmjuofxyPPmn+LgGyNBncpwbcF7eNHkxj6J9mGTKSQCAk+RRpeTrVunxCn17MwoFmFn+x1fMIuXitJaoUSNn7bjvDTO2hbkIcebM6RUC0vF7XphXJ/k9Mys7g+J5eRF2Jtr0agKA1cTSfqviyEWI91y3eAaVLKvXUVlJSySAjxILABa5n5p1LKQ458jiBvPak2yjsvNTgZtPp/dp1X3aCmR3Fiw5SPG9DMYrEmxOHz2LnYC/eC9NLEJoWZ1Ow+tM/bApoChubl82f8WS5otaRpNZzWePOUuAmz4Rx676NIjRkhQSMlZo1gVuEJaB7EqoQk2pqKWO9lncke/vsLbNi0jYDcQmJ9c9D6HlAtL1eTmcBWg4Gj7jruy8jiSmNJa0utCKE0TjUagyq9sxAhQOnEFinlmboIqrbsGYOqb/YOakY8uc7m+l4jz7WpcgwqYrHXxW5IqQVhto5dpVzHOq/7SNK5cU1ePnpgqY83X/n6kKhzK+1KrEgjBTFBw9DFfgNFxkrPp5KUY2QVmOgx72ebtwWSwYPhL7DQW0e3wkI+P5eMe8m/rlpTT3ECVyqSuAHrLrotgPNYx6SOSim71WVwNdUQEpS9qk1sVVTxVZXUChLv4lToEwgba7CdlUEMZxZeIeso9DyH0JjRWgS+ZgvpOoGHxUKU0UjffCP9CgroUfLpvcybo+oavDwPLt7D2v3nKh0P1yQS3fxtNvlqGZVXwIvnUy4uXaQ+JChO4kpFEi8C/YGNUsqeQoihwO11F5ZSokaz022ocnb6sOdsjq81H/MKxu+8oTAMnVeidR8bYSY5JxqR05E7fBexy+DO/7d35+FVlNcDx7/nZiMhG4Gwh1VKBVQUZNG6UFEQUeoGtKhgrRtad+uCP0WsttTaUkFL3agLilFRS7EIAm6AIGtAFmULqxAJSQgh+/n9MZNwE7IByb03N+fzPHly7ztzZ955nyQnM/POOc8nxHNWbi6v790PwH+K0umQlkK3pt3qTX6vQPXokFP56oefyDxSUPpYQZs959OmzavMjwplbMdh3BZ3GrLgKbvMahqEmgaoAlU9ICIeEfGo6kIRmVinPTMn7EQe/gs9fTjMqPg+WW72GVyfH8Ijnh856PFw0J1pmCvCE80SKPxkVJn9BHp+r0CyYW8WM5bt4PHLu9M6PpI1T1wCwFPfPMV7m7ZzQauN9CwK5dc972HwGTc5HzpjhB97bIzv1DRAZYhINPAlMF1E9uNMODMB6IRygrkqO/NqPORJCubeT5PiXGcSBxAR0ognIm/gid1H0OjP0PAMEFD1EJHVhaFZLSE7DbYutMkVFVi2LZ2bXv+W6IhQbruwM63iIun1Vi/yi7ySeUZHkAyEr33xaIAypoGoacGgYUAOcC8wB9hCxQlkTT1WVSXRJv2uI+xXk8GtJEpcEivOeJJ7VnUnc39v8g53RRFUQwClrWYxLCSZlZO7w4e3OtPZ0aPT2lOSfXpsgeaz9fu4/tWlJMZE8P7t55SmLZpz1Rz6R3d0CtkBjUIiuKzjZXx69af+7K4xflHtGZSbuPVjVR2IU/L99TrvlfGLmp55lcz16w0s7H+YAX/9/JhLg5tCD3FzbG/yo3dy4eEc7jmYQecC96S74IhzRtVAz6I+WrWb+99bQ/fWsUwbczZN3bRFALnLX2NF1hbweAj3hJFXlB/QyTyNqUvVBihVLRKRHBGJU9XM6tY3DUvHZk7y3IouDeZKPuOa38Vr8bFcFdWKK7MPM/ZgJs2Lihp0NdfW8ZEM6JrIpJFnEh1x9FcwMy+TsXs/pdjj4dL2g7jp9JsDPpmnMXWppvegcoG1IjIPOFzSqKp31UmvTHDQcG4sjuWanXt4KT6WGbEx9MjL45pDhyEkHLL3Q3Rzf/fSJ1SVpdvS6depKX06JtCnY8LRhUWFsGsZ/83Zxu78DF4eNI3eLZ3Zu5ayyDRkNQ1Qs90voLTkjOVWMdUad+hKng55mYfSM7gu6xAtCosgJIyPY2LJ2vQuI3reSmZeZtA9P1XZ82txkWHOTD3vB7MjYiEvi9+MXUr/y9+nU3wnP/TYmMBTXT2oYUBbVX3Bfb8MSMQJUg/VffdMfVDZtPaEqDA2Jgzmvj1FTGj8Aa0L9iNxSXDR4yxKX8b/1r3E9O2zaUsYK7NTg+r5qcqSvmYeKXCCk9ezaG9FFNNXIuny4xo6NdD7csZUpLp6UIuAkaq6032/GvglEA1MU9WLfNLLagRzJon6Lr+wmIlzNvLq19uYMKw7N/TvULrsrDfPoqC44JjPBMPzUx0enl3psu0tHnJnNcLH0Y15LLEpI7MOMa4w2lIWmXrDF5kkqptmHl4SnFxfq2q6m4OvYdYxN8clPNTD/w3txoxb+jGqb3sA0g87ZxefXv0pl3a8lFAJcdb1hHNZ0sDgn1LtThBZ1iiC8c0S6Hsklz8cONigJ44YU5HqAlSZOg7eFW1xLvUZUyP9OjUlxCP8lJ3HoElfMu7DtcSEJRAdFk2RFhPmCaOguIDGG2bTbNELzsSBYBWdyNawUO5pnkj7ggL+tj+NMHAeYjbGlKpuksRSEblZVV/2bhSRW4FlddctE6ziIsO46sw2/OvLrcz4didhrdajhX3JzuhDo/glLIpM4bGN/2Z86iJCu18FS6bUuwwUBUXFlS6LIhcUXo2LI1yVF/alEVusR+uNGWNKVXcPqjnwEZAHrHSbewERwK9UdV+d97AG7B5U/bNw435u/Pe3x7SHJXxJoxafMPTwEf64P42QMgvrR+buv8zZyIufb6lw2T8iX2GYLqTgF/ew+7v36XCwfgVfY0r4PZu5qu4HzhGRXwIlpTVmq+qCuuyUCX4Dfl7x808F6ecjUsR/m39KSLMEJvyUfvQ6dD3IQLFkywH++cUWRvROYuI1p5dZVrxuJv/6bDkZp44lfuB4Ogwc748uGlNv1LSi7gLAgpLxifwDA7g3LJkXm8QRCjzuHaQCeCJBRk4+9yWvpkPTxjx+ebeyCzN38bcvHub1JvE0b9+Tq/3TRWPqlZo+qGuMT92msRQezCQzpNw8ngCdSKCqPPrhWtIO5TFz7Dk0dlMYpeWk8eAXD3Dej1t4PTqCke2HcFXXa/3cW2PqBwtQJiDJRY9zp/swqwDpHg9NJBQJ0IkEqtC1RSw9k+I5vW18afvUlKms3L+KFR7lgpjOPHT+04hYEhZjasIClPGbyjJQRIR60NOudXJpzZ9A+qE9jGjbmotb9uPBkvYA4/EIdw/sUvq+fF0ngC8ObaHv233r/UPIxvhKlbP46gubxRc8Js//gefmfc9Tv+rB9f2cB3tVlYnfTmT6hunc2PI87k0ajHS7ws89deQXFnP7WysYc24Hzuty9NHAtJw0/rr0zyzYMY9clIiQCAa2G8gDZz8QNPkGTcMWCJkkjPGpOwacwoCuibyzdAeF7vNEIsJDZz/EiJ+NYNqPXzH527+gxZU/a+RLf//se+Zv3M/hvLIPFsdGxLIrazt5qoR7wsi3uk7GHDcLUCageDzCpJFn8t5t/Qn1miAhIjza71Gu7nAZL5PJB5tnkpaTxpg5Y/xWL2nxlp+Y+sUWRp6dxOAerUrbVZUnFz9JSsb3XND2At6+7B2Gdx3OgSMH/NJPY+oruwdlAk5cZBgAR/KLePOb7dz0i06EeASPeHj8/GdISujCoBZ9mTT396zMXO+XLOgZOfnc9+4aOjRtzP8NLTulfNqKSczaOos7TruV285ysoNZXSdjjp8FKBOwPtuwj2c+2Uh2biH3XdIVAI94eHHNi0xaOal0veRNySRvSvZpFvTk5Tv5KbvslHKAz1MXMGndaww+kset7Yf4pC/GBCu/XOITkWdFZKOIpIjIhyIS77XsERHZLCKbRGSQP/pnAsPlZ7RmeO+2PL9gMws2Hs2qNeeqOQzpcGnpf1ehEsJlHS/zaRb0m8/rxMd3nltmSnluYS5PfvUwp+bnMaHPo0hTKzxozMnw1z2oeUAPVT0d+B54BEBEugEjcdIqDQZeFJGQSrdigt6EYT3o1iqWe99dw870HAASoxJpHB5NEYJHoVCL2J2V6pMJCKkHDrPjQA4iQvfWcWWWNdr3HVN3pvJ8XG8izxpd530xJtj55RKfqs71evsNcI37ehgwQ1XzgG0ishnoAyzxcRdNgGgUFsLU63oxdPJXPDwzhem/6wdAem46w7sO58r47tz/9aOsPrCOdzZM59enjqrV/VdWur1ZdDjLH7sYUpIpmD+BLwsOcFFuPl3DY+CKF8AexjXmpAXCPajfAu+6r9vgBKwSu9w204C1axrF1Ot70TY+qrRt0oCj96BmHc7ivqUT+NfySQztfAUx4TG1tu/KSrf/lJ0PKcnorLt4OjaSDxKa8c7uH+lRkAM/zAvohLbG1Bd1dolPRD4TkXUVfA3zWmccUAhML2mqYFMVPkksIreIyHIRWZ6Wllb7B2ACyjmdm9GuaRSqysYfs8osC+s1mufaDObN1K3EbFvku07Nn8DbkSF8EBvNzRmZ9MjPh6I8J+O6Meak1VmAUtWBqtqjgq+PAURkNDAUGKVH01nsApK8NtMW2FPJ9l9S1d6q2jsx0Yr7NhQvfr6FYVMWsW53Zpn28CHPkZTYHf3wFp5f/BQfb/64zvuyOD+NvyQ0YcDhHO486NWfAM64bkx94pdLfCIyGHgIuEBVc7wW/Qd4W0T+BrQGumCVe42XkWcn8dzcTQyd/PUxy05vfAsf9N3I2r3LeOX7dwn58HaGhibUTTFAzxEebJ7IKfkF/DntQNn/9AI047ox9Y2/ZvFNAWKAeSKyWkSmAqjqd0AysB6YA9yhqkV+6qMJQE2jIyiuJH1kyuEmhLU8nec3fUvv3DzGJSYwp/AAzLoLUpKPe18rUg9Wuqw/W5nQ4VdMTj9ElHc+Syvdbkyt8UuAUtVTVDVJVXu6X7d5LXtaVTuraldV/Z8/+mfqsfkTiMw/wpR9afTMy+PhxKbMC+O47wst25bO9a8uJTy0/K9IEQmNNjEt/Fkuysml9ZBJEJcEiPO9HpSkN6a+CIRZfMbUHvf+T5QqL/6Yxh0tEikSOa77QitS07lx2jJaxjVixi39aB7TyCk8+OWDtItpx6wts9h35mTad7kUGsVZQDKmjliAMkFF49oimTsBaKzKaz/udy4TRDUlIzeD+EbxVX5+1Y6DjH7tW5rHNuKdm53gBG7hwX0rWLFvBTd0u4H2p42s4yMxxlg2cxNU3oy6AQ2NLH3vARAPS4qzGZR8IV+nLqg0C7qqMn7WehKiQ/nHqA7ERim93urFaa+fRvKm5NLnHd5Y/wa93urls2MypqGyAGXqnWbR4RW2R4WHMH57d54OuZ386DaU3hca9gLdzhhN+9wj3L3wbsZ/8QdW7lvJlFVTAEjNSuWZpc8wdv5YClr+ieyWDzJq7uWs2reKOVfNoXfcKYg7ESIiJMLnef+Maaisoq4JKku3HuCuGas4nFfE1w8NID7qaDDr9caZ5GvhMZ8JlTDQUDo3aUe72Ha0jWlLUkwS57U5j5Zpm5kwewzvN44kLCScguICru16rc/LexgTaHxRUdfuQZmg0rdTUz656zxWpB4sDU55hUVEhIYw55q5PPP1Yyzcs5gigXA89G0xgEXLzyUqpAkvXPkLWsQ2Orqx1CUwfTjpLRIZ3nko13a/gfe+f89vBRKNaWgsQJmg0zQ6gku6twRgzrq9PPPJRqb85kxOb5tIk5g2FIsQjlCgRUT9MJ9Onr5MuqV/2eC0YylMvwZiWzNpxGyIaQFY4UFjfMkClAlqiTGNKCwq5up/LiY8xENR4nq0sC/ZGX1o1+QT0sK2klYQSfvds+GNCc509OhEOJLpZIQYPas0OBljfMvuQZmgl5GTzwPvpfDZhn0VLFWu8Czi+ahpUHjEq13g0onQ91ZfddOYesUX96BsFp8JevFR4bx8Q2XTwoU/hCaXC04ACosn13XXjDFVsABlGgSpooBga6lk0oNlJTfGryxAmQZvj1ZSKt6ykhvjV0FxD0pE0oBUoBlgc4DLsjFxhbc8pcLrfPFk05b9xSJH/2FTpXhHlqb+lKPpvuuh39nPyrFsTI5VMibtVbVOi/EFRYAqISLL6/qmXX1jY1IxG5dj2Zgcy8bkWL4cE7vEZ4wxJiBZgDLGGBOQgi1AveTvDgQgG5OK2bgcy8bkWDYmx/LZmATVPShjjDHBI9jOoIwxxgQJC1DGGGMCkl8DlIgkichCEdkgIt+JyN1u+7MislFEUkTkQxGJ9/rMIyKyWUQ2icggr/btIrJWRFaLSIWJ+UTkQXf5ahFZJyJFIpLgLosXkffd/W4Qkf51ffyVCbBxudftwzoReUdEGlW0jbrmhzGJE5FZIrLG3d+NXstGi8gP7tfoujzuqgTKmIhITxFZ4raliMiIuj72ygTKmHgtjxWR3SIypa6OuTqBNCYi0k5E5rp9WS8iHarsvKr67QtoBZzlvo4Bvge6AZcAoW77RGCi+7obsAaIADoCW4AQd9l2oNlx7PtyYIHX+9eB37mvw4H4hj4uQBtgGxDpvk8GxjSEMQEe9dpWIpDu/lwkAFvd703c100a+Jj8DOjitrcG9vrr9ydQxsRr+T+At4Ep/hiPQBsT4HPgYvd1NBBV1bb8egalqntVdaX7+hCwAWijqnNVS0uffgOU5JwZBsxQ1TxV3QZsBvqc4O5/DbwDzn85wPnAq25f8lU14wS3e9ICZVxcoUCkiIQCUcCeE9zuSfHDmCgQIyKC84uUDhQCg4B5qpquqgeBecDgkzy8ExIoY6Kq36vqD24/9gD7cf4w+VygjAmAiPQCWgBzT/KwTkqgjImIdMMJiPPcvmSrak5VGwqYe1Duqd6ZwNJyi34L/M993QbY6bVsl9sGzqDMFZEVInJLNfuKwvmj8oHb1AlIA6aJyCoReUVEGp/godQqf46Lqu4G/grswPmvOFNV/frLBj4bkynAqTgBeS1wt6oWV7Ndv/HzmHj3ow/OWdWWEzqQWuTPMRERD/Ac8OBJHkat8vPPyc+ADBGZ6f6dfVZEQqrqb0AEKBGJxvmjeI+qZnm1j8P5b2R6SVMFHy+ZJ3+uqp4FXArcISLnV7HLy4FFqqV51kKBs4B/quqZwGHg4RM9ntri73ERkSY4/011xLl001hErjuJQzppPhyTQcBqnOPuCUxxz7Sr2q5fBMCYlOyvFfAmcGP5wOVrATAmY4FPVHVnBZ/xiwAYk1DgPOAB4GycE4MxVfXZ7wFKRMJwBm26qs70ah8NDAVGqXvBEieSJ3l9vC3uJSf30gKquh/4kKpPSUdS9jLWLmCXqpb8V/E+TsDymwAZl4HANlVNU9UCYCZwzskc18nw8ZjcCMxUx2ace3E/r2q7/hAgY1JymXw28JiqflN7R3j8AmRM+gN3ish2nKsQN4jIn2vtII9TgIzJLmCVqm51Ly1+RHV/Z9VPN+7csRDgDWBSufbBwHogsVx7d8revNsKhACNgRh3ncbAYmBwJfuMw7km2rhc+1dAV/f1eODZhj4uQF/gO5x7T4IzkeT3DWFMgH8C493XLYDdOFmcE9xfuCbu1zYgoYGPSTgwH+c/c7/8zgTamJRbZwz+nSQREGPibmNNyf6AacAdVfbdzz9Mv8A5dUzBOSVcDQzBuSm306ttqtdnxuFc394EXOq2dXIPfA3OH9RxXuvfBtxW7odlRgV96Qksd/vyEX6amRWA4/IksBFYh3P5JqIhjAnO5Ym5ONfQ1wHXea33W3e/m3EuZzWIn5PKxgS4Dijw2t9qoGdDHpNyfRqDfwNUwIwJcLHbj7XAv/Ga8VjRl6U6MsYYE5D8fg/KGGOMqYgFKGOMMQHJApQxxpiAZAHKGGNMQLIAZYwxJiBZgDJBR0T+LiL3eL3/VERe8Xr/nIjcV8v7zK7N7bnb7CkiQ7zejxeRB2rwORGRBd5ZHk6iD+Ei8qWbi9EYn7IAZYLRYtyMF25OtGY4Dx+WOAdY5Id+Ha+eOM+rHK8hwBr1SmdzolQ1H+chXL+V0DANlwUoE4wWcTQlU3echwUPiUgTEYnASWS5SkSiRWS+iKx0a9wMAxCRiSIytmRj7pnL/e7rB0XkW3Fq6DxZ0c4rWkdEOrg1cF4Wp0bOXBGJdJed7a67xE2guU5EwoEJwAhxau+UBIhuIvK5iGwVkbsqOf5RwMc12O/n7tnml+46Z7uJPH8QkT96be8jd5vG+JQFKBN01MkXVigi7XAC1RKc7M39gd5AintmkAtcqU7yywHAcyIiwAzKnjEMB94TkUuALjj5x3oCvcony6xmnS7AC6raHcgArnbbp+E8hd8fKHKPIR94HHhXVXuq6rvuuj/HScbZB3jCzbFW3rnACq/3le0XIF9Vzwem4gS1O4AewBgRaequsw4nuacxPmUBygSrkrOokgC1xOv9YncdAZ4RkRTgM5ySAi1UdRXQXERai8gZwEFV3YFT4O0SYBWwEidYdCm336rW2aaqq93XK4AO4lQxjVHVkj69Xc1xzVanTs9POHWXWlSwToI6dX9KHLNfr2X/cb+vBb5Tp3ZQHk7+tSQAVS0C8kUkppq+GVOr7ManCVYl96FOwzkD2AncD2QBr7nrjMIprNdLVQvczNMlJe3fB64BWuKcUYET0P6kqv+qYr8VriNOHZ48r6YiIJKKSxtUpfw2KvodLhQRjx4teVHRfstvr7jcesXlth2Bc8ZpjM/YGZQJVotwygikq2qROjWu4nEu8y1x14kD9rvBaQDQ3uvzM3DKj1yDE6wAPgV+69bVQUTaiEjzcvutyTql1KnKe0hE+rlNI70WH8Ip0X28NuEk9qwV7qW+kpIrxviMBSgTrNbizN77plxbpnt5DJwCbb1FZDnO2dTGkhVV9Tuc4LBbVfe6bXNxLsEtEZG1OIGrTACpyToVuAl4SUSW4JxRZbrtC3EmRXhPkqiJ2cCFx7F+dQYAn9Ti9oypEctmboyfiUi0qma7rx8GWqnq3SexvVbAG6p6cS31bybwiKpuqo3tGVNTdg/KGP+7TEQewfl9TKWaMtjVUdW97rTy2JN9Fsqd7v6RBSfjD3YGZYwxJiDZPShjjDEByQKUMcaYgGQByhhjTECyAGWMMSYgWYAyxhgTkP4fbmmMVBACuVkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mask = ((wl_ / 1000) > 2.20576) & ((wl_ / 1000) < 2.20586)\n",
    "wl, flux = wl_[mask], flux_[mask]\n",
    "flux = flux / max(flux)  # Normalize maximum to 1\n",
    "\n",
    "f_slope = slope(wl, flux)\n",
    "f_grad = slope_grad(wl, flux)\n",
    "f_adj, new_wl, new_flux = slope_adjusted(wl, flux)\n",
    "\n",
    "# plt.figure(figsize=(15,10))\n",
    "ax1 = plt.subplot(211)\n",
    "plt.plot(wl, flux, \"o-\", label=\"Spectrum\")\n",
    "# plt.plot(new_wl, new_flux, \"+--\", label=\"Diff adjusted\")\n",
    "plt.ylabel(\"Normalized Flux\")\n",
    "plt.ticklabel_format(axis=\"x\", style=\"plain\")\n",
    "ax1.ticklabel_format\n",
    "plt.legend()\n",
    "\n",
    "ax2 = plt.subplot(212, sharex=ax1)\n",
    "plt.plot(wl[:-1], f_slope, \"s--\", label=\"FFD\")\n",
    "plt.plot(new_wl, f_adj, \"o-.\", label=\"FFD(shifted)\")\n",
    "plt.plot(wl, f_grad, \"*--\", label=\"Numpy\")\n",
    "plt.ylim([-30, 32])\n",
    "plt.xlabel(\"Wavelength (nm)\")\n",
    "plt.ylabel(\"Gradient\")\n",
    "ax2.ticklabel_format(useOffset=False)\n",
    "ax1.tick_params(labelbottom=False)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"spectra_gradient_example1.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an wavelength offset between the dy/dx (blue) and gradient (green) due to the wavelenght points. The wavelength and flux are adjusted to the center of each difference (orange).\n",
    "\n",
    "We later show that the difference difference between the blue and orange creates a change in precision of around 0.1%. (e.g. the wavelength difference is squared)\n",
    "\n",
    "\n",
    "The gradient is less sharp then the dy/dx method and so produces a lower precision (higher rms).\n",
    "We also show that the difference in RV between the gradient and dy/dx method is between $2-7$% in the given bands!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xd4k1X7wPHvSZruSQdCC5Q9y6wMEZmCCCgiQxAVFPH9vQ7Eia8DERCcDEERB7hQQAURFGWUIYJQaNmrQIWWVVoKLV1pc35/JC0daZu2SZO253Ndudo8T/I8d0LJnfOcc+4jpJQoiqIoiqPR2DsARVEURTFHJShFURTFIakEpSiKojgklaAURVEUh6QSlKIoiuKQVIJSFEVRHJJKUIqiKIpDUglKURRFcUgqQSmKoigOycneAVhLQECADA0NtXcYiqIoNcLevXuvSCkDbXmOapOgQkNDiYyMtHcYiqIoNYIQ4l9bn8Nml/iEEF8KIS4LIQ4Vs18IIeYLIWKEEAeEEB3z7XtECHHSdHvEVjECrI6Kp/vszTScso7uszezOirelqdTFEVRLGTLPqilwF0l7B8INDXdJgKfAAghagFTgS5AZ2CqEMLPFgGujornlZ8PEp+cjgTik9N55eeDKkkpiqI4AJslKCnlNiCphIfcC3wtjXYBvkKIOsAAYIOUMklKeRXYQMmJrtze++M46fqcAtvS9Tm898dxW5xOURRFKQN7juILBs7lux9n2lbc9iKEEBOFEJFCiMiEhIQyB3A+Ob1M2xVFUZTKY88EJcxskyVsL7pRysVSynApZXhgYNkHk9T1dSvTdkVRFKXy2DNBxQH18t0PAc6XsN3qXhzQHDedtsA2J43gxQHNbXE6RVEUpQzsmaDWAA+bRvN1Ba5JKS8AfwD9hRB+psER/U3brG5oh2BmDQsj2NcNAbg6adBpBf1a1bbF6RRFUZQysNk8KCHE90AvIEAIEYdxZJ4OQEq5CPgNuBuIAdKA8aZ9SUKI6cAe06HeklKWNNiiQoZ2CGZoB2MX14G4ZO5ZsIOv/o7lyd5NbHVKRVEUxQI2S1BSytGl7JfAk8Xs+xL40hZxlaRtiC+9mwfy+fbTjLstFA+XajOPWVEUpcpRtfgKeaZvU66m6fl6p80nSSuKoiglUAmqkA71/bijWSCfbT/Njcxse4ejKIpiVk2ogqMSlBmT+jYl6UYW3+5SrShFUWyvrMmmplTBUZ0sZnRq4EePpgEs3naah7o1wN1ZvU2KothGbrLJrWqTm2wA7mlXl+R0PQkpmVxJNd4SUjKZu/FksVVwcgd9VQelfvIKIYKklJcLbWsupazW9YAm9W3K8EU7WfbPWSb0aGTvcBRFqaaKK7n23IpoXli5n2yD2ToFZlW3KjiWNA22CyFel1KuABBCPA88BrSyaWR2Fh5ai+5N/Fm09TQPdmmAm7O29CcpiqKUUXFJxSDhPz0bEejlQoCn8Rbo5UKgpwt3z99GfHJGkedI4IlvInlxQHOaBHnZOHLbs6QPqhfwkBBipRBiG9AMY5Xxam9S32ZcSc1k2e6z9g5FUZRq6GDcNYS54m5AsK8bL93VgvHdGzKkXV26NfanSZAnPu46XhzQokgVHFedhoFtbmFHTCL952zjxZX7ia/iLapSE5SpusN6oBsQirECeaqN43IInRvWolsjfxZtPUVGoSa4oihKRfxx+CIjPv0bb1cdLk4FP4rddNoSS64VroIT7OvG7GFt+WRsJ7a91JtHuzfkl+jz9H5vC9PXHiHpRpaNX41tCON82RIeIMQG4ALwDMa6eF8C26SUL9g+PMuFh4dLW6you+t0Ig8s3sXUIa0Y372h1Y+vKErNIqXks+2nmfX7MdqG+PLZw534OyaR9/44zvnkdOr6uvHigOYVHuwQn5zOvI0n+HFvHO7OTky8oxGP3d6QDUcuWeVcQoi9UsrwCgVZ2jksSFBDpZSr8913Al6RUk63ZWBlZasEBTDq053EJt5g64u9cdWpvihFUcpHn2Pg9dWH+GHPOQaF1eGDke1s/pkSczmF9/84wfrDF/Fw1pKVY0Cfc/Nz302nZdawsDInqcpIUJZc4ltd6H62oyUnW5vUtymXrmeyfM+50h+sKIpixrV0PeOW7OaHPed4qncTPhrdoVK+8DYJ8mLRQ51Y9d/b0BtkgeQEjr1Ia6kJSgiRIoS4brplCCFyhBDXKiM4R9GtsT+3hvrxyZZTZGarvihFUcrmbGIawz7ewe4zSbw/oh0vDGiORlPM6Agb6VDfD322wew+Rx2ebkkLyktK6W26uQL3AwttH5rjEEIwqW8zLl7PYIVqRSmKUgaRsUkM/XgHV1Kz+OaxLgzvFGK3WKraIq1lLnVkuuTXxwaxOLTuTfzp1MCPj1UrSlEUC/0SHc+Yz/7Bx03Hqv/eRtdG/naNx9wiraWNGLQnSypJDMt3VwOEU8wS7NWZsRXVlIe/3M2Pe+N4sEsDe4ekKIqDWR0VnzdCztPViZSMbDo3rMWnYzvh5+Fs7/DyBkJYe8SgrVhSSWJIvt+zgVjgXptE4+B6NA2gQ31fPo44xYhO9XB2UrV2FUUxKlxTLyUjG60QjOwU4hDJKVf+RVodnSV9UOPz3R6XUs4sXJuvpshtRcUnp/PTvjh7h6MoigMxV1MvR0rmbDxpp4iqvmJbUEKIjyjhUp6U8hmbROTgejYLpF09XxZGxDC8Uwg6rWpFKYpS/Eg4Rx0hVxWU9OkaCewt4VYjGVtRTYi7ms7PqhWlKIrJLT6uZrc76gi5qqCkPqjvpJRqSVkzejcPop6fG/9bdYgpPx10+I5GRVFsL9TfnQvXClYYd+QRclVBSS2o3bm/mC73KSa/RJ/n0vVMcgyyWq9mqSiKZfafS2bXmSR6NQsoUMC1PCWElJtKakHln+bc3daBVCXv/XGcrJyCM7Kr42qWiqKUzmCQvLHmMP4eLnw0piNerjp7h1RtlNSCqnFznSylOkMVRcm1cu859p9L5n93t1DJycpKakG1EEIcwNiSamz6HdN9KaVsa/PoHFRdXzezC4GpzlBFqVmupel5Z/1xwhv4cZ+6emJ1JSWolpUWRRXz4oDmBSbkAWg1QnWGKkoN8+GG4ySnZTHt3s6I4pbGVcqt2AQlpfy3MgOpSgqXC3Fz1pKelUO7er52jsy+8pd5KcvIxvI+T1Hs6cj563yz61/Gdm1A67o+9g6nWrKk1JFiRv5yIZdTMuj93hZm/36UTx+y6fpdDqtwmZfckY1AicmmvM9TFHuSUjJ1zSF83Z157s5m9g6n2lIJygqCvFz5T8/GfLDhBLvPJNG5YS17h1TpzJV5Sdfn8Oaaw6Rk6MnKkWTnGNDnGMjKkehzDGTnGFi2+6zZ56kRkYojWx0dz57Yq8weFoavu+PU2atuVIKykgk9GvHdP2eZue4Iq/7bvdIXI7O34kYwJqfref2Xw0W267QCnVZDWpb5pUvUiEjFUaVk6Hn7t2O0C/FhZHg9e4dTrZVUi+8gJdfiq7Gj+MxxczbOGH9+5X5+PXCee9vXnG//8cnp6LSaInPDAG7xdmHtMz3QaTV5SclJI/I6lLvP3mx2RKQQ8Pn204zuXB8PF/U9SnEc8zed5EpqJp8/HF7jvohWtpLmQQ3GuNTGetPtQdPtN+BH24dW9dzXIZjWdb15d/1xMvQ1Y1HDNfvPc9fcbYBEpy34n9VNp2XKwJYEeLrg46bD3dkJnVZTYLSTuQXUnLUaGgZ4MGPdUW6bvZkPN5zg6o2syng5ilKik5dSWLIjllHh9Wr8oKjKUGyCklL+axrJ111K+ZKU8qDpNgUYUHkhVh0ajeDVQS2JT05nyY5Ye4djU9cz9ExeHs0z30fRNMiTDc/15L3h7cpc5mVoh2BmDQsr8Lx3h7dl0/O9+Pm/t9G5YS3mbzrJbbM389avR9SlP8VujAMjDuPurOrrVRYhZckFI4QQ0cBTUsq/TPdvAz6WUravhPgsFh4eLiMjI+0dBgATvtrDP6eT2PJiL/w9XewdjtXtPpPE5OXRXLyewdN9mvBU7yY42XDZkZOXUvhk6yl+iT6PRsDQ9sE0qe3J13//q4amK5Vm3YELPLlsH9Pvbc1D3ULtHY7dCSH2SiltOmzZkgTVCfgS8MHYJ3UNeFRKua/UgwtxFzAP0AKfSylnF9rfwHTsQCAJGCuljDPtewcYZHrodCnl8pLO5UgJKuZyKgPmbmNM5/pMH9rG3uFYjT7HwNyNJ/hkyynq1XJnzqj2dKzvV2nnj7uaxufbz/DtrliyC3V3uem0qjCnYjNpWdn0/WArfu7O/Pr07WhV31OlJChLVtTdK6VsB7QF2ksp21uYnLTAQmAg0AoYLYRoVehh7wNfmwZcvAXMMj13ENARaA90AV4UQnhb/rLsq0mQJw92qc+y3WeJuZxq73Cs4nRCKsM/+ZuFEacY3imEdc/0qNTkBBDi586b97QmwKvouju5Q9MVxRYWbI7hwrUM3rq3tUpOlajU4VFCiNrA20BdKeVAU5LpJqX8opSndgZipJSnTcf5AbgXOJLvMa2AyabfI4DV+bZvNa1HlS2E2A/cBayw7GXZ36S+TVm1L57Zvx/l80dutXc4ZVKwsoMrtzUJYO3+C7joNCwa25G72tSxa3yXCq25kys+OZ0dMVe4rbG/KjujWM3phFQ+236aYR2CCQ+teXMc7cmSjoOlwB9AXdP9E8CzFjwvGDiX736caVt++4H7Tb/fB3gJIfxN2wcKIdyFEAFAb6BKTTjw93Thv72bsPHoZf4+dcXe4Vgst7JDfHK6aa2rDFZGxlHPz431k+6we3KC4ovyagQ8+Pk/DPvkbzYfu0Rpl68VpTRSSqb9egQXJy1T7m5h73BqHEsSVICUcgVgADC1aiwZQ23uK2zhT4wXgJ5CiCigJxAPZEsp/8Q4nP1v4HtgJ1BkdV8hxEQhRKQQIjIhIcGCkCrX+O6hBPu6MXPdUQyGqvFhaa4iBEBqVnaxS1pXNnND0910Wt69vy0zhrbh8vVMHl0ayeCP/mL9oQtV5r3Pb3VUPN1nb6bhlHV0n71ZLYZZyfLe/1d+Y+uJBPq1DCLIzKVlxbYsmQF5w9SqkQBCiK4YB0qUJo6CrZ4Q4Hz+B0gpzwPDTMf1BO6XUl4z7ZsJzDTtWwacLHwCKeViYDEYB0lYEFOlctVpeemu5kz6IZpVUfHc3ynE3iGVqrhh3BeSzV9Ws4fCxXoLj+IbdWs9VkfF8/GWU/zn2300q+3Jk72bMLhtXX7df75SC9OWpxCuqk9oX4Xff4D1hy+yOipevf+VzNJRfPOBNsAhjCPuRkgp95fyPCeMlwP7YmwZ7QHGSCkP53tMAJAkpTQIIWYCOVLKN0wDLHyllIlCiLbAMowDNIq0onI50ii+/AwGyX0f7+DS9UwiXuiFm7O29CfZUafpG0g0Myk22NeNHVP62CGi8ssxSNYeOM/CiBhOXEolwNOZa+l69Dk3/+ZtOfrP3Add7vkGt61D4o0sLl3P4PL1TC6lZHDpeiaXr2ewKiqezMLDFKma/wZVUXHVTdT7X1BljOIrtQUlpdwrhOgJNMd42e64lFJvwfOyhRBPYey/0gJfSikPCyHeAiKllGuAXsAsIYQEtgFPmp6uA7abOrqvYxx+XmxycmQajeC1wa0YsWgnn28/zdN9m9o7pGKtP3SRq2lZCAH5v7e46armxEStRnBv+2CGtK3Ln0cu8vT3UQWSE+SO/jtmUaumLC2hrGwDb/921Gwh3OdWRPPcimgKX3kUAgI8XcwmJzC2pC5fzyDIW11qsiW1YrbjsKQFdQp4T0q5KN+2tVLKwbYOriwctQWV6z/f7GXbyQS2vNjLIa9lrztwgWd+iKJtiA8jw0NYsPlUtZsE23DKumKLS7au602jQE8aBXjQKNCDRgGeNAz0wNPFqcSWUO/mQcQkpHIq93b5BqcSUjmblEZOCX1fT/dpQpC3K7W9XKjt7Uptb1cCPJ1x0mqK/QYPxqTbu3kgI8Lr0adFEDobTpCuqVQLyjKOMlH3GMZRdWnAE1LKLCFElJSygy0DKytHT1CxV27Q78OtjAivx6xhYfYOp4BfouOZvDyajvX9WDL+VrxcdfYOySaK++DxcNES3qAWp6+kEnc1vUDrMcjLheR0PVlmWjUaQYFWkLNWQ2iAO02CPGkc6Mm3u/7lalrRiw2lfdAVlxAn39mUpBt6ftoXR0JKJgGezgzrGMLI8Ho0CfJUCz9ayQ+7zzLF1OeXS00EL8ohLvEBaVLKUUKIlzBedhtJCVXOFfNCAzx4qFsDluyIZdPRSySkZDrEh8iPe+N46cf93Bpaiy/H3VqtK4e/OKC52Q/+mUNvfvBk6HM4m5TG6YRUTiXc4MyVG/y4N87s8QwSXhnYIi8hhfi5FSj51DjQ0+z5SrtcWtogkBf6N2PL8QRWRJ7jy7/OsHjbaUL93YlPTs+7hKkGVpTftXTjl4pATxeupDrG/9OaypJPIwEgpXxXCLEXY5+Smq1WDk2DPAG4nJIJ2P9DZPke4zfF2xr78/nDtzr8AI6KKu2DH4wjL5vV9qJZba+8bTtPJRZ7yeeJno0rdL6Snlvc45y0Gvq1qk2/VrVJSMlkVVQc764/TrbBXP+aWvixLDL0OXy2/TQ9mgbwzWNd7B1OjWdJgnoj9xcp5SYhxADgEduFVH0tjDhVZJu9PkS+3fUvr60+xB3NAln8UCdcddU7OeUq6YO/OMW1vCwZOFKe85VFoJcLE+9ozKzfjpndrzr2y+aH3We5kprFU72b2DsUhZIXLGwhpTwGxAshOhbavda2YVVPjjI6aMmOM0z79Qh9WwSx8MGONSY5lVdFWkKVpa6vm9lWXh0HmVxdFWRlG/h022luDfWjSyN/e4ejUHIL6nngceADM/skoIazlJEjfIh8tu00M387Sv9WtVkwpiPOTmoUmCVs3RKqKHOtPABfdx2Z2Tm4OKkvIaX5eV8cF65lMPt+tVi4oyhpwcLHTT97m7mp5FQO5kr0AOi0govFFEC1po+3xDDzt6MMCqvDwgdVcqpOzC38OCI8hCMXUnhqWRT6HPNzqxSj7BwDH285RdsQH+5oGmDvcBSTki7xDSvpiVLKn60fTvVm7lJR/1ZBrIiMY/BH21k4pqNVLy3kH3bs6epESkY297Sry4cj29l0gUHFPsy18trU9WHqmsNM+iGK+Q90UP/uxfj1wHnOJqXx2qBOqhK+AynpEt+QEvZJQCWocjD3ITKmSwOe+GYvYz7/h1fvbsn47qEV/k9SeC5NSkY2WiHo1SxQfUjVII/cFoo+x8CMdUdx0uxnzqj2aj2jQgwGycKIU7S4xYt+LWvbOxwln2ITlJRyfGUGUpM1re3F6qe68/yK/by19gj745KZPaxthYZ9v7P+WJH+iBwp+WDDCYZVgaK1ivVM6NEIfY7knfXHcNII3hvRTiWpfP44fJGYy6nMH90BjXpfHIpFszJNK9y2BvJ686WUb9kqqJrI21XHp2M78cnWU7z/53GOX0xh8UPh1Pd3t/gYUkr2x13j+3/OcqGYPi017Lhm+r9ejcnOMfDBhhPotBpmDQtTH8YY/898tDmGhgEeDAqz/1pnSkGWrKi7CHDHuGjg58BwYLeN46qRNBrBk72b0CbYh2e+j2LwR9uZN7oDvZsHlfi86xl6fomKZ9nucxy9cB13Zy3uzlrSsoqu61TcYn9K9fd036bocwzM3xyDk1YwY2ibGt/fEnH8MkcuXOe94W1Vq9IBWdKCuk1K2VYIcUBKOU0I8QGq/8mmejYL5NenbueJb/fy6NI9PNevGSG+bry/4UTe4IoX+jejQYAH3/9zlrUHLpCuz6F1XW9m3teGe9rVZdPRy+WeXKpUX5PvbEZWjmTR1lPotBqmDmlVY5NUbusp2NfNoacQ1GSWJKjca0JpQoi6QCLQ0HYhKQD1/d35+f9u45WfD/DBhhMFCpPGJ6fz3Ir9SMDDWcvQDsGM6VyfsBCfvOdXhcmlSuUTQvDyXc3JzjHw+V9ncNIIXh3UskYmqZ2nEok6m8z0oW1UVXgHZUmCWiuE8AXeA/ZhHMH3uU2jUgBwc9YyZ1R7Io5f5lp6weWwJODrpuOvKX3wLKbAq6NPLlXsQwhjUso2SD7/6wxnrtzg2MXrnE/OqFFfZD7aHEOQlwsj1KAhh2XJgoXTTb/+JIRYC7jmLsuu2J4Qguvp5tdqvJauLzY5KUpJhBBMHdKKE5eus+nY5bzt9i5gXFn2/pvEztOJvDaopSr15cAsGSShBQYBobmPF0IgpfzQtqEpuYorkaQGPCgVIYQgNjGtyPaaUAV9weYYank4M6ZLfXuHopTAkguvvwLjAH/AK99NqSTmSiSpAQ+KNVxIrnnTEQ7FXyPieAKP3d4Qd2d1BcKRWfKvEyKlVNUT7UgNeFBspSa2zhdsjsHL1YmHujWwdyhKKSxJUL8LIfpLKf+0eTRKsdSAB8UWzFVB1wh4/s6mdozKdk5cSmH94Ys806cJ3q46e4ejlMKSBLULWCWE0AB6jCvsSimlt00jUxTF5gq3zn3ddVxN03PGTN9UdfBxRAzuzlrGd1czZaoCSxLUB0A34KCUUpb2YEVRqpbCrfMXV+5nQUQMXRr6c3s1Wnoi9soN1uw/z+M9GuHn4WzvcBQLWDJI4iRwSCUnRakZpt3bmiaBnjy7PJrLKbZfp6yyfLLFWD3jsR6q9VRVWJKgLgBbhBCvCCGey73ZOjBFUezD3dmJhQ92JDVTz+Tl0eQYqv530/jkdH7aF8cDt9YjyKvyVrBWKsaSBHUG2AQ4o4aZK0qN0Ky2F9Puac2OmEQ+joixdzjltjoqnu6zN9N99mayDZLQAA97h6SUQYl9UKZJup5SyhcrKR5FURzEyPB6/H0qkTkbT9C5YS2rrvZcGQov2gnw7vrj+Lk7qxGxVUSJLSgpZQ7QsZJiURTFgQghmHlfGA38PXjmhygSUzPtHVKZvPfH8SKLduZWyVCqBksu8UULIdYIIR4SQgzLvdk8MkVR7M7TxYkFYzpwNU3P8yv3Y6hC/VHFVcOozlUyqhtLElQtjEts9AGGmG6DbRmUoiiOo3VdH14f1JItxxP4bPtpe4djsVt8zA+GqM5VMqobS6qZj6+MQBRFcVxjuzZg5+lE3vvjOOGhtejUwM/eIZUq2NeVC9cKDpNXNSyrllJbUEKIECHEKiHEZSHEJSHET0IItYCKotQgQghmDWtLHV9Xnvk+iuS0LHuHVKL1hy4Q+W8yd7WuTbCvGwII9nVj1rAwNUCiCrGkksQSYBkwwnR/rGnbnbYKylr0ej1xcXFkZFSfyYb25OrqSkhICDqdqmFWE/m46VgwuiPDF/3Niz8eYPFDnRxyJd4rqZm8uuoQbYK9+WhMR7VabhVmSYIKlFIuyXd/qRDiWVsFZE1xcXF4eXkRGhrqkP+RqhIpJYmJicTFxdGwoZqJX1O1q+fLy3e1YMa6oyz9O9bhatpJKXl11UFSMrL5fmR7lZyqOEsS1BUhxFjge9P90RgHTTi8jIwMlZysRAiBv78/CQkJ9g5FsbPHbm/IrtOJTF97hI+3nOJKSqbDLAGzOjqePw5f4pWBLWhWW9UTqOos+XrxKDASuIix7NFw07ZSCSHuEkIcF0LECCGmmNnfQAixSQhxQAixJX/flhDiXSHEYSHEUSHEfFHOLKOSk/Wo91IB499BnxZBSAkJKZlIbi4Vvzoq3m5xXbiWzhu/HKZTAz8m9GhktzgU6yk1QUkpz0op75FSBkopg6SUQ6WU/5b2PFMVioXAQKAVMFoI0arQw94HvjYtiPgWMMv03NuA7kBboA1wK9CzDK9LURQbWhhxisIzouw5CVZKycs/HSQ7R/LBiHZoNerLVHVQ7CU+IcQbJTxPSimnl3LszkCMlPK06Xg/APcCR/I9phUw2fR7BLA69/iAK8b6fwLQAZdKOV+FrY6Kt8mqtTNnzmTZsmVotVo0Gg2ffvopXbp0sULE8Pbbb/O///3PKsdSFEs52iTY73efY9uJBKbf21rV26tGSmpB3TBzA3gMeNmCYwcD5/LdjzNty28/cL/p9/sALyGEv5RyJ8aEdcF0+0NKebTwCYQQE4UQkUKIyIr2jeTW7YpPTrfqJYudO3eydu1a9u3bx4EDB9i4cSP16tWr0DHze/vtt81ul1JiMBisdh5Fya+4ya7FTY61pbOJacxYd4TbmwTwYBe1jHt1UmwLSkr5Qe7vQggvYBIwHvgB4yKGpTHXxi58VeAFYIEQYhywDYgHsoUQTYCWQG6f1AYhxB1Sym2FYlwMLAYIDw8vsQbLtF8Pc+T89WL3R51NJiun4Ad6uj6Hl348wPe7z5p9Tqu63kwd0rqk03LhwgUCAgJwcXEBICDAuABcaGgoo0aNIiIiAoBly5bRpEkTEhIS+M9//sPZs8Zzzp07l+7du5OamsrTTz9NZGQkQgimTp3Knj17SE9Pp3379rRu3ZqZM2cycOBAevfuzc6dO1m9ejWtW7cmNTUVgB9//JG1a9eydOlSxo0bh5ubG8eOHePff/9lyZIlfPXVV+zcuZMuXbqwdOnSEl+XUrOZWyoewFkruJ6hr7Tl1A0GyQs/7kcrBO8Mb4tGXdqrVkrsgxJC1BJCzAAOYExmHaWUL0spL1tw7Dggf1MhBDif/wFSyvNSymFSyg7Aq6Zt1zC2pnZJKVOllKnA70BXS19UeRROTqVtt1T//v05d+4czZo147///S9bt27N2+ft7c3u3bt56qmnePZZ48j9SZMmMXnyZPbs2cNPP/3EhAkTAJg+fTo+Pj4cPHiQAwcO0KdPH2bPno2bmxvR0dF89913ABw/fpyHH36YqKgoGjQo+dvk1atX2bx5M3PmzGHIkCFMnjyZw4cPc/DgQaKjoyv0upXqbWiHYGYNCyswCXZ89wacv5bBg5/9U2kTeZf8HcvuM0m8MaQVwaqEUbVTUh/Ue8AwjC2UMFOiKIs9QFMhREOMLaMHgDGFzhEAJEkpDcArwJemXWeBx4UQszC2xHoCc8t4/gJKa+l0n72ZeDPXz4N93Vj+RLdyn9fT05O9e/eyfft2IiIiGDVqFLNnzwZg9Oj0r9LvAAAgAElEQVTReT8nTzZ2xW3cuJEjR252012/fp2UlBQ2btzIDz/8kLfdz898qZkGDRrQtatluXzIkCEIIQgLC6N27dqEhYUB0Lp1a2JjY2nfvn3ZX7BSYxReKh7g9iaB/N+3+3hg8S6+ndCFAE8Xm50/5nIq764/Rr+WQQzvpIrbVEcltaCeB+oCrwHnhRDXTbcUIUTx18pMpJTZwFPAH8BRYIWU8rAQ4i0hxD2mh/UCjgshTgC1gZmm7T8Cp4CDGPup9kspfy37y7PciwOa46bTFthmrbpdWq2WXr16MW3aNBYsWMBPP/0EFBy2nfu7wWBg586dREdHEx0dTXx8PF5eXkgpLRrm7eFRsIM4/3MKV9TIveyo0Wjyfs+9n52dXcZXqSjQt2VtvhgXTmziDR5YvIvL121TxSU7x8DzK/fj5qzl7WFhagpENVVsgpJSaqSUblJKLymld76bl5TS25KDSyl/k1I2k1I2llLONG17Q0q5xvT7j1LKpqbHTJBSZpq250gpn5BStpRStpJS2nyJeXOXLKxRt+v48eOcPHky7350dHTepbfly5fn/ezWzdhK69+/PwsWLCjweHPbr169CoBOp0Ov1xd7/tq1a3P06FEMBgOrVq2q0GtRFEv0aBrI0vGdOZ+czshPd9pkZN+irafYfy6ZGUPbqCXcqzFVBySfoR2C2TGlD2dmD2LHlD5WGWKemprKI488QqtWrWjbti1HjhzhzTffBCAzM5MuXbowb9485syZA8D8+fOJjIykbdu2tGrVikWLFgHw2muvcfXqVdq0aUO7du3yBldMnDiRtm3b8uCDD5o9/+zZsxk8eDB9+vShTp06FX49imKJro38+eaxziSmZjHy052cS0qz2rEPn7/GvE0nGdy2DoPb1rXacRXHI6SsOguQlSQ8PFxGRkYW2Hb06FFatmxpp4hKFhoaSmRkZN6ovqrCkd9TxfEciEvmoS924+6s5bsJXWgU6Fmu4+Sfo6jVCNx0Gra91Ac/D2crR6xYSgixV0oZbstzqBaUoig20zbEl+8f70pWtoFRi3dx8lJKmY9ReI5itkGSmS3ZekLVhazuVIKyk9jY2CrXelKU8mhV15sfJnZFAKMW72JhRAzdZ2+m4ZR1dJ+92exkeCklF69l8NfJK0xdc6jIfKusHIPdyioplceSauZVmqWj35TSVZfLwUrla1rbi+VPdGPowr8KJJb45HRe/ukAh89fI8DThZjLqZy8nMqpy6mkZJY8ktReZZWUylOtE5SrqyuJiYn4+/urJFVBuetBubqqEVNK+TQM8MBVp+VaesHEk5lt4LPtZwAI9HKhaZAn93UMpmmQJ02CvJi8PJqLZoarF1duSak+qnWCCgkJIS4uTq1hZCW5K+oqSnldvp5pdrsAot/oj4970RJJUwa2KFJWyVpzFBXHVq0TlE6nU6u/KooDqevrZrZiS11fN7PJCcib7mGLlQYUx1atE5SiKI7FXJFZS1pD5soqKdWfSlCKolQa1RpSykIlKEVRKpVqDSmWqjaVJIQQCUDuUvQBwBU7huNo1PtRlHpPilLvSVHqPSkq9z1pIKUMtOWJqk2Cyk8IEWnrEhxViXo/ilLvSVHqPSlKvSdFVeZ7oipJKIqiKA5JJShFURTFIVXXBLXY3gE4GPV+FKXek6LUe1KUek+KqrT3pFr2QSmKoihVX3VtQSmKoihVnEpQiqIoikNSCUpRFEVxSCpBKYqiKA5JJShFURTFIakEpSiKojgklaAURVEUh6QSlKIoiuKQVIJSFEVRHJLDJighhKsQYrcQYr8Q4rAQYpq9Y1IURVEqj8OWOhJCCMBDSpkqhNABfwGTpJS77ByaoiiKUgkcdkVdacycqaa7OtPNMbOpoiiKYnUOm6AAhBBaYC/QBFgopfyn0P6JwEQADw+PTi1atKj8IBVFUWqgvXv3XlEr6gJCCF9gFfC0lPKQuceEh4fLyMjIyg1MURSlhhJC7LX1yroOO0giPyllMrAFuMvOoSiKoiiVxGETlBAi0NRyQgjhBvQDjtk3KkVRFKWyOHIfVB3gK1M/lAZYIaVca+eYFEVRlErisAlKSnkA6GDvOBTHptfriYuLIyMjw96h1Eiurq6EhISg0+nsHYpSDTlsglIUS8TFxeHl5UVoaCjGqXNKZZFSkpiYSFxcHA0bNrR3OEo15LB9UIpiiYyMDPz9/VVysgMhBP7+/qr1qtiMSlBKlaeSk/2o916xJZWgFEVRFIek+qCUGiN8xgaupGYV2R7g6Uzka3eW+7harZawsLC8+6tXryY2NpZ77703r28mICCAjRs38uabb/LZZ58RGBjIjRs3CAsLY8aMGbRq1arc51eU6kolKKXGMJecStpuKTc3N6Kjowtsi42NpUePHqxdW3RmxOTJk3nhhRcAWL58OX369OHgwYMEBtq0aoyiVDkqQSnVyqhPdxbZNrhtHR7qFlrqc5NuZPF/3+4tsG35E92sFZpZo0aNYt26dSxbtoxJkybZ9FyKUtWoPihFqaD09HTat29P+/btue+++/K2b9++PW/7zJkzi31+x44dOXZMFUlRlMJUC0qpVirS4qnl4Vyu55u7xAcUe4mvsKpQsFlR7EG1oBTFzqKiomjZsqW9w1AUh6MSlFJjBHg6l2l7Zfjpp5/4888/GT16tN1iUBRHpS7xKTVGRYaSW9OcOXP49ttvuXHjBm3atGHz5s1qBJ+imOGwCxYKIeoBXwO3AAZgsZRyXnGPVwsW1kxHjx5Vl8fsTP0b1EyVsWChI7egsoHnpZT7hBBewF4hxAYp5RF7B6YoiqLYnsP2QUkpL0gp95l+TwGOAsH2jUpRFEWpLA6boPITQoRiXBvqn0LbJwohIoUQkQkJCfYITVEURbERh09QQghP4CfgWSnl9fz7pJSLpZThUspw1cmsKIpSvTh0ghJC6DAmp++klD/bOx5FURSl8jhsghLGhWa+AI5KKT+0dzyKoihK5XLYBAV0Bx4C+gghok23u+0dlFLFHVgBc9rAm77GnwdWVPiQWq02r+Ze+/btiY2NZcuWLfj4+ORt69evHwBvvvkmwcHBtG/fnqZNmzJs2DCOHCk4MHX48OGcPn262PP16tULc1MqIiMjeeaZZwDIzMykX79+tG/fnuXLlzN37lzS0tLK9Lq2bNnC4MGDAVi7di1Tp04t0/PtLXzGBkKnrCtyC5+xwfgAG/wtKNblsMPMpZR/AWq5TsV6DqyAX58Bfbrx/rVzxvsAbUeW+7DWXG7j8OHD5OTk0KhRozLHER4eTni4cVpKVFQUer0+L67Q0FDGjh2Lu7t7mY8LMGjQIF5//XVefvnlch+jspW4vIqN/hYU63LkFpSilN2SQcXffnnq5gdSLn06/P6y8fcbiUWfY2OjRo2if//+LFu2DIDvvvuOe++9F4CcnBzGjRtHmzZtCAsLY86cOXnPW7lyJZ07d6ZZs2Zs374duNniuXz5MmPHjiU6Opr27dszb948zp8/T+/evenduzcAf/75J926daNjx46MGDGC1NRUANavX0+LFi24/fbb+fnnm92+Qgh69eplUfHbKmHTW+b/Fja9ZZ94FLNUglJqjpxM89vTkyp0WGsut7Fjxw46deoEQHR0NPHx8Rw6dIiDBw8yfvz4vOdkZ2eze/du5s6dy7Rp0wocLygoiM8//5wePXoQHR3NpEmTqFu3LhEREURERHDlyhVmzJjBxo0b2bdvH+Hh4Xz44YdkZGTw+OOP8+uvv7J9+3YuXrxY4Ljh4eF5ybDKuxZXtu2KXTjsJT5FKZfx64rfN6eN8VJOYT71jD89/Et+fjGsudzGhQsX8uryNWrUiNOnT/P0008zaNAg+vfvn/e4YcOGAdCpUydiY2PLFO+uXbs4cuQI3bt3ByArK4tu3bpx7NgxGjZsSNOmTQEYO3YsixcvznteUFAQ58+fL9O57CUzO6fE/Skut+CVeaHoDqGByC+h/YMkZF3nxW0v8n7P9wlwC7BRpEpJVAtKqTn6vgE6t4LbdG7G7XaUf7kNNzc3MjIyAPDz82P//v306tWLhQsXMmHChLznuLi4AMYBGtnZ2WU6n5SSO++8k+joaKKjozly5AhffPEFYLyUV5yMjAzc3NyK3e8o/k28wfBPbq6sLJyu41b/U4Q2hdYilgW6ebyZMoR0ClWx1zqDT31YOxnmtWPBH0+y79I+Ptn/SSW/AiWXSlBKzdF2JAyZb2oxCePPIfPt2ileeLmNli1bEhMTA8CVK1cwGAzcf//9TJ8+nX379pX7PF5eXqSkpADQtWtXduzYkXeetLQ0Tpw4QYsWLThz5gynTp0C4Pvvvy9wjBMnTtCmTZtyx1BZImOv8m/iDbxdjReInP03oXWPpU3QUn52nsqt2pM89sBIvg6YTJQIZI+rK6sC63Fp4CyYFEWHRg0JC9Dx8/WjSCQrjq8g7KswOn3byc6vrOZRl/iUmqXtSLuP0ippuY1BgwaxZcsW+vXrR3x8POPHj8dgMAAwa9ascp9z4sSJDBw4kDp16hAREcHSpUsZPXo0mZnGfrkZM2bQrFkzFi9ezKBBgwgICOD222/n0KFDeceIiIioUAy2lJmdw6H4a3RqUIv7O4XQu0UQd67qhlfOzZF8sb7xhPvWQadxou6pacT7xLPAyw0wtgofSvcmc91R5MXHyPH9Ea3zVRASpKRlZjZdE9pDVhocW2scTHEtDnxCjC1wNfLPJhx2uY2yUstt1EzVbamH9PR0evfuzY4dO9BqtfYOJ8+lS5cYM2YMmzZtKrLP3v8G/ybe4Mll+zh1+QbbXupNoJfx8mdCWgLvbvsff17chUGAExoGNLyLsa0eYsmhJYR4hRDiFUKwZzD1POux8p8UFkbEkmOQuNRehc5vN0gtQmQTlKVlaVwq9YZNh3WTC44A1LnZvSVuD1ViuQ0hxAgp5crStimKUjo3NzemTZtGfHw89evXt3c4ec6ePcsHH3xg7zCKWHfgAlN+OoBGI5g/ukNeckJKAg6u4vi57RicdTgJLTnSgIezJ20C2vBBr6Kv5bk7YXjHUO54LwLhlIr+ahf0yZ3R+f5DvC6Fu7JG8Mu2NxHoaZz/ibnD02tYgqoM1rjE9wpQOBmZ26YoigUGDBhg7xCKuPXWW+0dQgFSSt745TDf7PqX9vV8WTCmAyF+pgnE+gxYNZFv4zZzxt+Plr7NmN7jbVaeWMmV9CslHre+v/EYGfEP5W3LvHRz6sAM50yO1L2F168kMeRGvsocani6TZQ7QQkhBgJ3A8FCiPn5dnljXGxQURSlwsJnbDBbFcJNp+HxHg15OfggTkufKNAnpEfwS53G9LulAx/0+hCN0PBa19cqHMu7GS68KDL5X1AA+66nMCXpKi4S43kVq6tIC+o8EAncA+zNtz0FmFyRoBRFUXIVV7IoXW/g1XqH4NdJRUoW6YbMZ2mLu9EIDRphvcHKQX3e4Itfn+Ejz0y+9PXhsIsLH1+6TEDjvlY7h3JTuf/lpJT7pZRfAU2klF/lu/0spbxqxRgVRVHMK1SyKEmjYbaXK+mb3sLT2RN3XdnrBgZ4Ohe77/Nr4TgNmc9kgzcfXUzAX+jwqdUMor6Ggz+W6yUoxbNGH1RnIcSbQAPT8QQgpZRlr3aZjxDiS2AwcFlK6fiTLxRFqXz5+n6ygGdrB3DE2ZmhFy7RopyHjHztziLbsnMMPPNDFDPWHcVlaBcemnyIXkAvgKw0kr+7j+83Pc8EjY7khrepChRWYo227xfAh8DtwK1AuOlnRS0F7rLCcRSlgIS0BMatH1dqh7mlhBA8//zzeffff/993nzzTascWymFqe9HAtMCahHl6sqMK0m0cKtt1dM4aTXMHdWBfi2DeH31IQ7EJd/c6ezOxi4P8bGfN4/ueIk5ES+pChRWYo0EdU1K+buU8rKUMjH3VtGDSim3ARWr4qkoZiw6sMiqHyAuLi78/PPPXLlinYSnlEFrY03CL3y8WePlyX+vJnNXlrRJ+SpnJw0LxnTkw5HtCAv2KbBveKuxOAknol2c+fVKpKpAYSXWuMQXIYR4D/gZyCsXLaUsf10WRSmn8evHF9k2IHQAD7R4gE7fdCLLcLPDfcXxFaw4vgIn4UTUw1FczbjKc1ueK/DcJXctKfWcTk5OTJw4kTlz5hSpWj5u3DgGDx7M8OHDAfD09CQ1NZUtW7YwdepUateuTXR0NMOGDSMsLIx58+aRnp7O6tWrady4MePGjcPV1ZXDhw9z6dIlPvzwQwYPHkyPHj346KOPaN++PQDdu3fnk08+oW3btmV+zxxdgKez2YESdT2Ao2u47lmbb3ydGZh6g/9IHxhiu8oOrjotwzoaW23HLl4nLimdfq2MrbU/h//Jm9tfZfuFnUgBzsKJO0MH8MKtL9gklprAGgmqi+ln/hnFEuhjhWOXSAgxEZgIONSkRsUxrR66mrHrxnI18yoSiUDg5+rHhLAJpT+5FE8++SRt27blpZdesvg5+/fv5+jRo9SqVYtGjRoxYcIEdu/ezbx58/joo4+YO3cuYFz8cOvWrZw6dYrevXsTExPDhAkTWLp0KXPnzuXEiRNkZmZWy+QkpTTbJwRAxCzYegbvh39hWVBT/F39EU6ulRbbO78f46+YK3z2cDi9mgcR6B5Ibe8QuCjQSIleZuPh5Kb6oSqgwglKStnbGoGU89yLgcVgLHVkrzgcUXFzRwI8nY3/4Q+sqJb1xEpq8dTzqkffBn358cSPOGud0efo6degHw+1Mk7K9HP1s6jFZI63tzcPP/ww8+fPt7ji96233kqdOnUAaNy4cd5yGmFhYUREROQ9buTIkWg0Gpo2bUqjRo04duwYI0aMYPr06bz33nt8+eWXjBs3rlxxO7KkG1mMXryL1we34vamhT7kE0+R8Pc8fmvenYcb9iS4hCrstjJ3VAdGf7aLJ77Zy5Lxt3Jb4wCSMpIY2XwkI+rewcqTq0jISOJY0jFa1CrvkI2arcJ9UEKI2kKIL4QQv5vutxJCPFbx0JSKsGi562vnAHlzuesDKyo3SDvI/QBZdvcyRjYfSWJ6hbtL8zz77LN88cUX3LhxI2+bk5NTXrFXKSVZWTf/XXKXzADQaDR59zUaTYElNAovgSGEwN3dnTvvvJNffvmFFStWMGbMGKu9Dkfxzu/HOJWQerN8US4pyVj3PJOC/FiYk0Bcqn2qOPi46/jmsc7Ur+XOhK8iiYxNYm7vubzW9TWa17+D1/rOoZFvI8aufYC/j6jCOuVhjUt8S4ElwKum+yeA5RhH95WbEOJ7jKM4A4QQccBUKWWFjqmYlLTctakVlZCWUC2Hys7tPTfvd2tUFsivVq1ajBw5ki+++IJHH30UgNDQUPbu3cvIkSP55Zdf0Ov1ZT7uypUreeSRRzhz5gynT5+mefPmAEyYMIEhQ4bQo0cPatWqZdXXYm+RsUksjzzHE3c0ovktXgX2ycOreD3lAIc8PZhzx2zqedWzU5Tg7+nCdxO6MGrxLsZ89g9ZOYYC+321HjQIzeKZPTOYnxTDbXu+q3ZXLWzJGqP4AqSUKwADgJQyGyh5OUsLSClHSynrSCl1UsoQlZysRxa73PU5OLMdDDlWH+lWUzz//PMFRvM9/vjjbN26lc6dO/PPP//g4eFR5mM2b96cnj17MnDgQBYtWoSrq7GfpVOnTnh7exdYCr46yM4x8NrqQ9T1ceWZvk3ztiekJTDut4f4cPtrrPf0YFKHZ+hb3/4VHIK8XVk+sWuR5ASQnFOHy7FP0cA1gKdjvuPvrARq2lWLirBGC+qGEMIf48AIhBBdgWtWOK5STgkpmSXuP2+oRbDG3KUtQact/yFLc/N7S+5IN2etM3vH7jXzHCU1NTXv99q1a5OWllbg/q5du/Lu566n1KtXL3r16pW3fcuWLXm/F97XvXt35syZU+S858+fx2AwFFgKvjr44/Aljl1MYdHYTni43PyIWnRgEfsSotnnpuGeOt151AqDW6wlyLv4wRkXchrw+4VLTHDXMyXQn/XnzuMupaqCbgFrtKCeA9YAjYUQO4CvgaetcFylHLaeSGDgvG159/Mvd32XZjeepJHU7RWyhPG6fpoQRLs4s8Lbl5hus+hxfgi++pt/FhopaZziTdN/TX0cB1bAnDbwpq/xp/oGaBdff/01Xbp0YebMmWg01Wth7LvDbuGHiV0Z0No4fLvTt50I+yqMFcdXIAEpBGsu7CD8O5suRWRVvslxfH7xMgsvJhiTUy5VBb1E1hjFt08I0RNojrHM0XEpZdkvtCsVdiohlXFLdtMsyIscg+Rqmj5vueteteexKOkIn2geIKPDBN7Q9+XolYPEasBg6oRP3e6LzK6Pq/dpnJyOIhAYgH89rjEq7Uc4UM94WaJQYU5AfQu0kaVLl5rd/vDDD/Pwww9XbjCVIDE1E39PF7o28s/btn7Yet7d8y6bz24iy6DHVetK3/p9q9b8Ip8QfK+dw9c0SGaZlyf1srPp4Rxo58AcW0WW2+gjpdwshBhWaFczIQRSyp8rGJtiodTMbDxdnGgc6MnCMR3p0yKI7ss7F1juOtInlTCf+jhrIwm5MZjIzARaNuhNf/8WtKzVkgaeTTl0VsOkH6IBDfqrXY2LtflvBfeT/KLvxZRNbyH16RQYU+YAlymklEVGuimVw5orcm89kcAT30Ty7WNdCA+9Oegj0D2QuEvRZBn06IQTmTmZeDh7VK3BO33fyPtypwdWe3kS46xjbp07ucPesTmwirSgegKbgSFm9kmMlSUUG5JS8vO+eKb9epilj3amY30/7g4zzqtZP/Q3ZqyfwObUMyAEGgTdg2/nre5v4e/qz5DGRf/ZGteCST9EF1ys7fxoMpGkAOg/YbGvN7E6Hc8nXSUgt1PYjpcpXF1dSUxMxN/fXyWpSialJDExMW/QRkVk6HOY+ssh6vq4ERZSsIzQqeRTHE6/RIjWnbkDv2LlyR+tVkfRmoqreFHLQ3fzC9ymt9Bdi+OzFHjc28CzFzYwd/8X3NFOzcwxp9wJSko51fSzeg0hqiJSMvS8tvoQv0Sfp3PDWtySv5M29TIuP08kMvsUaLXoNDqyDdnU8axTzm+dpg9+nxDgGn94uLPF3Y3/Xr3G6OspOHnY75tsSEgIcXFxJCQk2C2GmszV1ZWQkIov1rdo6yliE9P49rEuuDhp87bnGHJ44+838Hbx5pt7VxPgFsBr/tadHmAtxVa8AAwGyX7fO+kw2ZiofIDPLh1i4trRPBs1hznZWbRqObxaTu2oiIpc4nuupP1Syg/Le2yloOKqQgA8d2czngqIQrPkcWNLxsMf9JnsdNGQEuBDz+CePN3xaYuWuy7NH7dMZOLptxmQeoHZ/n686+/Hz14eTE+4Spt9X0PHyu8T0el0NGzYsNLPq1hP7JUbfLzlFEPa1S1SMeK7nW9zIOEAszq/WqU/tL/46wyzfj/K+yPa5dXy86ndhsX3/sgTa0YQt/0dtl2LyZva8XrX1+0csWOoyCW+3NlzzTEur7HGdH8IsM3sM5RyKS45ATwTGFVw4MKNK4BgQI+ptOs4mls8bgEsn5Ra3GUKFycNT+xvzIxGL/DgjSV8cimezQEhvFcrEOHSCNY8DUmnSej6f7z418vqW6Bise0xV3Bx0vDaoJYFd+Rk43V4DXdjYFDDQfYJzkrGdKlPxPHLPLdiP2lZOYzt2gAAn4DmnNTpmO3nBef+BNTUjvxERTs5hRB/AvdLKVNM972AlVLKSl3LKTw8XEZGRlbmKStN6JR1xe6Lrf2yqWSRsePvg1q+3J6WTleXIJh8yGoxSCmZt+kkczeeZEi7unw0ugMA2YZsnCTw2wu8c2YVB33rcECmMaL5CPUtULHY1RtZ+HkUWsl250L4438w8htodY99ArOiDH0OT363j03HLvPKwBY80bMxYJyA/H7k+2w8u5GsnCyc0NC/4V28eOuLDv0lTwixV0pp07H+1pioWx/jYpa5soBQKxxXsUS+AQqLfL35yscbZynpetW6AxeEEDzbrxmh/h40CrxZDcFJY/wT6nRtO1k+3iCNdejUt0ClNKmZ2cRcTqV9Pd8iyWn1waXIf95naJM7ES3NjcOqelx1WhY91Ilnl0fz4YYT3B1Wh3q13Al0D8RD54E+R48GyMbA0SuH8XHxKfWY1Z01EtQ3wG4hxCqMX+LvwzhZV6kMPiFw7RwrvDz52M+Xe1JSefrqNfCxTX2yoR2C836fv+kktzX2Jzy0FuuHref9yPfZEPsHepkDEtr7NWdO/0U2iaMy1dTK8LZgyXt5IeItZnlL2rpqGdqoF1Sj0Zk6rYb5D3Tg2MXr1Kvlnrc9t4jx/Q0HM/XvNzh6/QxPb3qa93u+j6ezpx0jti9rTNSdKYRYj3HJd4DxUsqoih5XMUrLyi52n45s6PsGmza8wEx/b3qkpfPmlSSEzs0mK4rml5qZzeqoeBZExPDe8Lbc2z4YD50H2dKAs3AiS+o5mHwcQdX/cLGoMryZycvha/y4kpqFcLqOa93vyYgfg8zxqtGJrbT3Uv76DNNqeSJxMf4tR8wAz6Bq9b5oNYLWdY2toxV7zvHGmkNk6I0jAD8nHngCnc8e/pY/sWzlfUy8fKFG/Y3kZ5UaKVLKvcD3wCogUQihVg+0knfXHy9mj2Sx2wI4s43tzXvRJgfev5yIzqceDJlv8z9iTxcnfvq/22hfz5dJP0Tz0aaTN5eyGPwDIxsOolPtcPzd/JGZqcSlVM+SLjkbpxVbGV6TeolAruLqvwGteyzOARsBteRJsTa9xRoXDTvc3Xg2KZng7JybE8GrqbNJaWToixaZ1V8L55V4Fx49uRuunSOnhv6NVLgFJYS4B/gAqAtcxtgndQxobYVj3wXMA7TA51LK2RU9ZlUipSQrx8D47qFMHVLo7TywAn7eDQFDeOO2p0jPTsddV/ZK2RXh5+HMN6OsY1cAACAASURBVI91ZspPB/lgwwkg91tgDJjmx/fc8wuTgt7ibVc9L90wMOJyHKIKfRP8N/Hm2k4FW0KedBAxaK4XTLw5wEZ3N86I62S0eAcpbv4nc671D861/kFnANYkQHbJS57UNCnX43mnXl06ZmTwQMrNArzVuV7dCwOasyAixsweQY/MJJw0cEWr4fFbgpiclMwd6TXrb8QafVDTga7ARillByFEb2B0RQ8qhNACCzF+6sUBe4QQa6SURyp67KpCCMHb94VhMBQaaXn9Auf/eJFXGzTmrTb3Uk9o8Kjk5JTLxUnLhyPbsSoq3uz+uBuCO1zr8Ef6CaZ7uPBPYC2mJsbh7cA1/K5n6MnJkfh5OHM64WaCcvbfhJN7LG1qf8bjKRfIcUnmA50vsc46GmXpee5qMhpgWoA/KVoNTno3tOSQ5ZRlTFQGQfM0Hdec9CzxcGLADS11swutTFONP4xL4+UdzOyEBOrrswte2vGp+ETgqqiuMK44IAGdhKdrB/JK4lUeqEF/I9a4xKeXUiYCGiGERkoZAbS3wnE7AzFSytNSyizgB+BeKxy3SvhkyykOxhlXLdFo8vXjSMnVNU/yhJ87x511pBtKXlqjMpRUYigHLQFJsXxy8RLPJl1lk4c7I+vW4aDIcahLN1JK/jmdyHMrouk8cyOLt58G4I5mgXg2fw2vllNwrvUPCEmsz2VeDdHyRqA/3/v5cd5Jh5NxtRkEsCwhmd3tpnA1ZiopKR0wIJAGJ/QCog3N+FcG82EtPwbUC+bBOrX5xtuLlNz3MF9VjoS0BMatH+eQZX2sTmRB3ze4I1sQmm81YSqhP9VRnZfGv4XAHANLL1zi9vQMZgbU4oPadTDIopcFqyNrtKCShRCeGCfnfieEuAwU37NvuWDgXL77cUCX/A8QQkwEJgLUr199ur3+OnmFd9YfI6lHwwJ1yRLSEnh+3UNk3DjNeTd3Fvf7mGZ+zewYqYWuxaEBHruWQqeMTKYEBpCs1VR6a6G4EWQezlqCvF0J+//2zju+qirb4991b256IAmB0IuoICKDgAji6AAq1rGCPMvYHQvP7owNZHSeM6BjmWGeimMbxwpSH0oRdHSUKi0gvQgJLSSQXm7uXe+PcxJukptA6j1J9vfzySf37LPPPvusT3LX2Xuv/VuZC3jC8xl/cWfg3dCRzLaPMWvvEtr7iskA3AhFLgG/i5KCbkTlXMGKx2/HveFz29lmQ+vOdC+dvvxkHhKWi/fI2ZbwbvwKJCyH/N3j+O7avcxfOokFkeG8nBjP5bl5gLBt0M0kFBwmKbJNuaSRzWFP2UsLthAXGUZOYfmvB3Hn0OGkSczMvZ6rr/hriwscqYrJJWN4LeZdxFtAtCqvHUxnUpsE3msVR8z8e7nn4jeaVYRjMOrDQV0JFAAPAzdiyUzVx6txMMuXm+tS1anAVLA26tbDPUNOTqGX33++npPaxvDoRb3KnXt95YusyUuFiAhePe8lBiYPDFEva0Z2RDKtig4A0L+omDmp+7B2vQjzFz/B4KGP4UMbXIesqgiyvGIfV7m/5/6odwjzFbLd42FqVCGLUl6mRIQBnhi6x3VgRc4uwt0evHi5YcBZjB9yl9VAvzFBv0STYsM5HCi8e/CqsvIuZ91D4sZ8/pH2BgX+/STEdYKRE/jT/vms/Gx4uXaaw56yrzcfYsrX27l9WA8mXNGn3LnfLXmQRXtKOKOoGM4JbsvmTFXqLXP853JG69bcUfwBruw0wlp35qmzH+a0zR8xcvN3MDyHdH9Rs9bvq5ODsteJZqvqBVgp39+vl15ZpAKBm3k6A/vqsX1H8j/zNrE/q4Dp955DpMcSzRz4r4EUl6bOsN+YHvr3w4T/p2l8YT2Tcw2vRL2L22cFBYQDhEWSmdCVCXvmErv3S05v15/V6WsaZbRQMewb4Bb5mEy8tAPyXMJ/oqIYm53DdRpLzwdX8NDXDzGm41mMPnX0CesaViceqqpsbHMxj2/pRWJMOE//6jSuOqMTT3UZwOdbpjFjyyfkY03jhEsYF3Yf1bTyHwVwOLeIx6evo3f7OH53cfmXriV7lvDl3iWM6/dbTu5/X4h6GFqC/Z2oKu//sJs/zhP+lTCF1+8eSJ+OrRDgmkG3QnYaxZ5Ibpv9X+zJ2dtsRtkVqQ+poznAzapar2neRSQM2AqMBNKAlcANqroxWP3mIHW0bGcGY6cu457ze/LEJb3LytPz03lp0X+zJDOFQperXMI2p7w1VbcB84M7zua09Pmw+Dk0KxV/q064L3gWzhjNgA/OtDb2VqA+Rwv7swp44YvNzF1nvd9EJM/Ek7AC75Gz6JsbR1zCd6TEFHFFbh4TMzJRoEiESFVAYOLReulHMDakZfH0rA2s23uUs3skMunafnRPiuG5pc8xfes0VEFQRrvbMP6KDyC+aU1lqyp3vr+K77YfZs64YfRu36rsXFZ+BlfPuITEuM58fMWneFyeEPbUmfz4cyb3fbiarAIvs+4/Zr9yL60BNOYou6lIHRUCKSKyCCgLeVLVB+rSqKqWiMg4YAFWmPk7VTmn5sLAbglMvKIPYweX/xLKLMxknTeTIpeLcHe4IxO2VTdaAKDDGA71+DXnTvqamFw3TxafxnUKC65bxKSlz/NV6jfWXg9gVPdRPDH4CTZlbKJNVBvaRbercX/yi0vYn1VIz7axxESEsWp3JrG9nkFcx9Y/whNXsDURUOW/CkoYnZ0DWHPLkaUvbg0cQda3U2tm3nsOn67ay2tfbcMlwqA/LiK39U9oyRA0uy994+dyyPMzS94cSp9uv6L9/rWQldYk1mi+357B4s2HGH95n3LOKT0/nbtmXkWWN48pXS4xzqkKBnZL5P/++5d8tHwPvZLjyspLlVsW//wVRf5jjqpbXDeyirKajUxSfYygbgk4LG1MVLU+p/uOS1MfQRUU+4gKd1cqzy/OY+y8saTl7uOSHpdwc5+by6aYXh3+agh6Wje2Hszh6ZkprNx9hMHdE9l2KIf8uM/wJKzApS5UfBQfOZs+uecR1W82KUe30q9EGZmdxUhXa7oNn1CtQsO0e85hzbypDNk1hWQySE3ozKb+17IhuhWfrFxIT18WO2JzKXIJohBXEM+htNvZNiamvCIEWBFkjbDpuZTiEj/hYa6g4sAd5ADuk1/FL34mZGRycV5+SPpYG77bls6wnknlolGfXziOafu+4YrwZP7nhsWh61wTY29mPuNnb+CFq8/gH5teZPrW6XjcHrw+L2fGn0K77INMvmY2EtOGmdtm0iWuCwOTBzZIMs/GGEHV2kGJyJVAZ1X9u328AmiL5aR+r6rT6q2XJ0BTdlCLNx3kyRkpfHjn2ZwS8JYE8Ozsscw8spG3fjmZs3teEqIe1i9+vzL9x1Re+HITR/O9RHb6AC2JKxfpNvVQKp1b72NxOCyODGNjRAQA1+UWkJE2mjn+cwOm6gZTdPAqxJPBpf6NvBj2LgvjPExOTCDHbe2k8Cj0Ki6iW5GPL+KiUXWj4sd7ZDBxedc7SnqoKvX6z6If5LVkF+sjI7giJ48nMzKJU7V0F+tRub4+KC7xk3okn5PalteRc8LUVFPmu23p3PPBj+QX+4io8H/T0bObhYdX44nvSMmZNzNq5z855BK6+ZRrulzAr899pmzWJT0/vc7BFU53UN8DY1V1r328FhgBxALvqurIeuvlCdBUHdTR/GIueuVbEmPCmT1uWLlsovN3zefxbx/nroguPHD9vGYXUpqZV8yA5xcFPdeRwyyMfZbYkiMA7He7WRITRSevj0fbJVHsCm6Lvx72MjxnP6sjIpgXG02fomL6FBdzcokfz4XP81DeRpJiO5YLdnDaSLQqB7Uz4gb8AlPjWzM1vhXJJT6m7dtPKz8Nuk5WGybN38w7/9nFV4+cX04Ude2htdz15W8oVD+IOHI91elsP5TLBS//O+i5/rKdWa1fhsKjFIiwKCaaz+NiWB0ZiRvh6aHjGX3qaJ5f9jzTtkyrU1ocp69BhZc6J5v/qGomkCkioZE1aIJMnLORzLxi3rn1rEqprl9f9zr92vbj3ovfa3bOCSCxYv6fAPaRRLT3aNlmgw4+HzdmW/I3X+5NY3xSW5ZFh+MXQVTp4i0hMut0+hZ8DcCAoiIGFAVuYhYYej+BruhEkzg6hX2aRGc5zH1HsxhWUMAPUZG08iuERaCqDTKNUxuW7czgjX/vYMzALuWc09YjW3l40b34/T7E5cLj0PVUp3Nyu6rVzdfqyRAeDYVHiVLl17l5/Do3j12eMGa26cCflv+J55Ye2wXk9C0MdVGSSAg8UNVxAYdt69Bui2H+hgPMWruPcSNOpm+n8oua7lXv8G5kb/5y7qQWu4Dsig8eoFBcksg3Refhw4X6w/DjYmfuMNYcvpW2sR2DN9YM5HIml4xBPVEA/KKomHuPZoM7go1D7uKmL29i59GdIVefyCrw8uhn6+iWGF1uv1NKegq3fXkLrqJsBhDJmFPH8NGlHzGm1xgyCjJC0tdmS/b+SkU9vCU8ciCVBdcu4NIel5blcYt0R3JZj8tYcO2Cxu7lCVEXB7VcRO6qWCgivwVW1KHdFsOinw7Sp0Mr7h9+crnyHzZ/TsnC8SQe3k77uE5VXN0CGDnBCgIIxBPF5JIxZQoN+bvvw3vkbDQsDz+uKq9pSnI5SbHBR5Y/RI9gfo+nyAxLRhFr7enKKWSddhl7svdw/eyreXz26DL1iVAwYfYGDmQX8urYM4mJODZBs3z/MuKKC3g/PYe3rpzOM0PH0yuxF88MecZxU6xNnqpexlp3LkuO6PP7HBsRHEhdpvgeBmaJyA3AartsIBABXFXXjjVHqtorNPRPi8sW6Zf9+w/cEyc8GOPhjpMvapZTe4FUtYs+KTYc+l1mHVQIXPhhTgKFaceuCVRoqOoaJ0e5VaS6kP3XvurOvet7ckan1rxx80A6xUdxDpDvzacYPz8WW6ORUEzd+PxKu7gIHrnwVPp3iQe7X9GeaO7ofSPXb19OXP8RkNijUfrTYhk5IXhEqv2SVpoWpyabzkNFfYSZj+BYao2Nqrqkzr2qBU0hSKKqxW+A3TfkkTnvQa5tl0Arv59P9h0gKizS8SHEhsZn0U8HeeTTtYSHuZhywwCG9mxjbeZe9RILd82nBD9hCqM6DuOxX/4xZG/H83fN588r/szbo96mZ3zPkPShueKELM+OjuJzGk3dQe1M/j3jIgtZHhnJR/sO0MvrtU44MITYEHp2pOdy9z9XkXa0gO9+N4K2cRG2+sR0POLG6/cyOreQEWfeCd1/ybBOwxqsLz6/8uhna7lpSDcGdU8EYMa2GUz8YSJntv0FU7J9xA1/BpL7HKclQ1PC6VF8hnrkQz3Kd9EJPHU485hzghadH8hQNT3bxjLr/mGs2n2EtnHWHrGMgoCpmw3/5PCOhbyb8jbLt3/IBV1G8LvBT9AhtkO93D/YG/ystftIig3nwWsOMHnlZIYVK6+s/D+iAJJ6QfKz9XJvQ8vBOKhGorqRqhsfA4r93JSVXT6TKDSL6DNDwxAX6WF4b0sGqt/EBWQXBmY0PgcPg/lD1L8Y0moZb+oS/rNnMb89ksVvaE14Had8qlKGPyqrmLzyIy7ML2LSwYOUxZ8ufx3anWamqw01oj4SFhpOgBmrg2ecBT+/dc/l9MI8fp+VXz7HSBOLPjOEjuzCyinYvITxVMGt3Nnzaubs3ce5+QW8lhjPt97D1iL6+s/qvR8luafxWK6PyYHOCY6lszcYaoBxUI3EJWe0JyaI1l5kh+lM7+BHL30Jrvy7teZUGkJsAiQM9cGOJXTwlfDKocN8sO8AI/MLwFvAN98+x4E8K0/Xie6fSs8pKjcbIGHZRHV9g/C2X9DVvYs/u9/llvS04FMzZrraUEMcOcUnIqOBicBpwGBVdXb0QzUcyi4kOiKM2IgwNj53cblzc1Pe46nVq7n1F/ci/QOS3xkM9YhmpZaNzPsXWVNzRQLPRvspmHk5d/e+iX3F2dVm7911OI83vtnBjDWp/OuOY4mtw9t8RVj0bsJidjNOZ3FlThGEx0JxbqU2zHS1oaY40kEBG4BrgDdD3ZG6UOj1cdc/V+F2CZ/fe06ZFE16fjoPLB7HjoyfGOCO5e5+d4e4p4bmTE5ARuNSIhQ+ynFxaXwRr218u6y8dP+Ux+Xh9r63k5aVzeq96fyceRSXq4TRgx6mS2I0sb2fQsRfrs3n2ybyfBs3KQOfqnYfjsFwojjSQanqJsAx2mK1QVUZP2sD61KzmHpzebn719e+zobMn/CIm0lDJpbJjhgMDYFv+Hi8Cx/F4y8sK/O6Iuk4fAJd5x7Gm/g+aVF5IEK4X+mY245Deipvrn8Tl99FtProFOsjQtw80W0DreIHEHtkBFGt5pPpduMXAX8Y3py+xOReBf2utW7ShDdLG5xBk/5mFJG7gbsBunZ1VqbRfy3fw7QfU3lgxMlcdHp7oHKqAS8+LvzhMcKXP+VIoUZD06E6RY6EITdBdDi6+A+QlUaGuy1zku7k9n5jSPloHhHRZ+CJWoHLL3jFz2B28mjGOqKi+6P71uMKSIjH3ocgLJJlD7/CcwsLmb7/e8LdHrx4uXHwqYwfYjunfmOMQzLUmZA5KBH5Cmgf5NTTqjr7RNpQ1anAVLA26tZj9+rEjz9n8tzcjQzv1ZaHLji1rHzWlbOY8s0TLDm8jkJX+VQDBkNdOG5G435jENthJAG3+AMDHSxdw9K8QrM9h+jsi+C21IUIFf6t1AdfTYR+Y8j0RDQZyRxD0yRkDkpVLwjVvRua5FaRXNSnPS9cc0ZZFtGtR7by6OJxxGfuoSjc0ySEGg3NF3dAPq3CtJvLPhcdvIoi4A/AbWELg1+cvQ+gnMhrU0tdYmgamDDzesTr8+P3K50Tovn7jQNoHWXtBJm7Yy43zruB3Jz9uMTFmJMuN6kGDM6nGlVsg6ExcOQalIhcDfwNK6/UPBFZq6qjQtyt4zJh9kbSc4p48+aBuF1Cka+ISSsmMW3rNAZJNC+mHSDp5jnQbShg3joNDuc4qtgGQ0PjyBGUqs5U1c6qGqGqyU3BOX28Yg8fr9jDKcmxZdMn69PX8/m2z7mjVR/e2rmZpIteKHNOBoMTqCr3lJW6ZIy1WdxsHjeECKNmXg+s2XOE699cxtknJfLebYPZn5dG5zhrGmT3mvfpPvsB6H+jpRTRhEPnDQaDoRSjZu5QqsrFsnHfUf6+9m+8k/IP3sv20z8jle4A8d3gspeNczIYDIYa4MgpPqcTzDmJO4eCNq/zVspbXJWbT+/MVECtn9yDsGlOo/fTYDAYmjLGQdWQykKZb+KO/YnoHn/DHfUzz+eUMDE9ncjAqdOSQqPkbDAYDDXETPGdAIVeH//emm79bEkvKw9vs5iw6N0kylHyVMjefR9XyePBGzFKzgaDwVAjWrSDKl1LkrBsIjt+TGHaDagvjqTYcN67bTCtts2k65qXiMhKJdGVxP7wESR16EF2+49Bjo2QsqOPABDd/X/haKfgzsjsHTEYDIYa0aIdVOlaUnibxbijdxORPBdv9plkeTJ4bNprjJKNPJyViQATOoaR6lkKLCVS/Xj8Sr7LjU9A/G6Kc/pSdPByuN5l9o4YDAZDPdCiHVRsr2cQ17FMpJ7W6/G0Xg9All/5ueDYEt34jCNE+pWueHjRexdfxO1CE9aA342KD/VF0iaqDfSzNdGMkrPBYDDUiRbtoPJ2/I6Idl8QFrcRcXlRv5uS/JMoPng5612PlBPKPKegNFVBMZMmTqTo64dIiiovlPnqcNs5GSVng8FgqDMt2kFpSSvUFwFSgvrDQHxocSL+4mQkuTNk7a18kb2WZIQyDQaDoWFp8WHmpakG8nffh/fI2UiYnap65ARr7SgQs5ZkMBgMjUaLHkElxYZzuEKqgdJy+l1mFZq1JIPBYAgJRovPYDAYDDWmMbT4HDnFJyIvishmEVkvIjNFJD7UfTIYDAZD4+JIBwUsAvqqaj9gK/BkiPtjMBgMhkbGkQ5KVReqaukGpWWAkWEwGAyGFkZTCJK4Hfg02AkRuRu42z7MFZEt9uck4HAj9K2pYOxRGWOTyhibVMbYpDKlNunW0DcKWZCEiHwFtA9y6mlVnW3XeRoYBFyjNeioiKxq6MW7poSxR2WMTSpjbFIZY5PKNKZNQjaCUtULqjsvIrcAlwMja+KcDAaDwdA8cOQUn4hcDPweOF9V80PdH4PBYDA0Po4MkgCmAHHAIhFZKyJv1PD6qQ3Qp6aMsUdljE0qY2xSGWOTyjSaTZrNRl2DwWAwNC+cOoIyGAwGQwvHOCiDwWAwOBJHOCgR6SIiX4vIJhHZKCIP2uWj7WO/iAyqcM2TIrJdRLaIyKiA8t0ikmKvXQUV5xOR1iIyV0TW2e3fFnDuFhHZZv/c0lDPfDycYhMR6S8iS+2y9SJyfUM+d1U4xR4B51uJSJqITGmI5z0RnGQTEekqIgvtvvwkIt0b5qmrx2E2mWyXbRKRv4qINNRzV0cIbJIglkTdehFZISJ9A85dbLe5XUSeOG7nVTXkP0AHYID9OQ5L3qgPcBrQC/gGGBRQvw+wDogAegA7ALd9bjeQdJz7PQVMsj+3BTKBcCAR2Gn/TrA/J7Rwm5wKnGKXdwT2A/Et1R4B518DPgKmtPT/G/v4G+BC+3MsEN2SbQKcA3wPuO2fpcCvWohNXgSetT/3Bhbbn912WyfZNloH9KmuLUeMoFR1v6qutj/nAJuATqq6SVW3BLnkSuATVS1S1V3AdmBwTW4JxNlvNLFYf1QlwChgkapmquoRLE3Ai2v9YHXAKTZR1a2qus3uxz7gENY/YqPiFHsAiMhAIBlYWOsHqgecYhMR6QOEqeoiuy+5GqLtIU6xiV0eifVFHAF4gIO1fKw6EQKb9AEW2/fbDHQXkWS7je2qulNVi4FP7HtViSMcVCD21MCZwPJqqnUCAtPdptplYP1hLBSRH8WSQgrGFKy3h31ACvCgqvqP027ICLFNAvsxGOsfbkcNH6FeCaU9RMQF/AV4vNYP0ACE+G/kVOCoiMwQkTViZSNw1/ph6olQ2kRVlwJfY8047AcWqOqmWj5KvdFINlkHXGPfbzCWJFLn47QbFEdt1BWRWOBz4CFVza6uapCy0nj5Yaq6T0TaYe2j2qyq31aoOwpYC4wAetr1vjtOuyEh1DYpvaeIdAA+AG6p6Lgak1DbA/gN8IWq7g3RkkIlHGCTMOCXWF98e7C0M28F3q7lI9UZB9ikHZbjKhW6XiQi5wW5vtFoRJv8GXhNRNZiOe01WKPKGn+/OmYEJSIeLON9qKozjlM9FegScNwZ6w2mdBoKVT0EzCT40PQ2YIZabAd2Yc2VVtluKHCITRCRVsA84BlVXVb7J6obDrHHUGCciOwGXgJ+IyJ/rvVD1RGH2CQVWGNP3ZQAs4ABtX+quuEQm1wNLLOnO3OBL4EhtX+qutGYNlHVbFW9TVX7Y73QtcWyS42/Xx3hoOz527eBTar68glcMgcYKyIRItIDOAVYISIxIhJntxkDXARsCHL9HmCkXS8Za6FwJ7AAuMiOQkmwr19Qt6erHU6xiYiEY/0h/lNVp9X1uWqLU+yhqjeqaldV7Q48hmWX40cjNQBOsQmwEkgQkdK1yRHAT7V/strjIJvsAc4XkTDbOZyPtfbT6DS2TUQk3v7eALgT+NYesa0EThGRHvb5sfa9qkZDEFVS8Qc4F2uotx5ruLwWuBTrLSQVKMJaYFwQcM3TWGshW4BL7LKTsOY/1wEbsZTRS+vfA9xjf+6ItcCdYhv4poB6t2MtCm4HbmvpNgFuArwBfVgL9G+p9qjQp1sJbRSfY2wCXGj3IwV4j4CIx5ZoE6yItTexnNJPwMst6O9kKLAN2AzMICAS2r7vVrvtp4/XdyN1ZDAYDAZH4ogpPoPBYDAYKmIclMFgMBgciXFQBoPBYHAkxkEZDAaDwZEYB2UwGAwGR2IclKHZISKviMhDAccLROQfAcd/EZFH6vmeufXZnt1mfxG5NOB4oog8dgLXiYgssTdY17UP4SLyrYg4SnXG0DIwDsrQHPkBS00asbTzkoDTA86XKk07nf5Y+0ZqyqXAOq1ezuaEUEvUczEQkjQrhpaNcVCG5sj32A4KyzFtAHJshZAILI20NSISKyKLRWS1WDlurgQQkUkicl9pY/bI5VH78+MislKsXDd/CHbzYHVEpLtY+XjeEisHz0IRibLPnWXXXSqW0OoGe6f9c8D1YuXeKXUQfUTkGxHZKSIPVPH8NwKzT+C+39ijzW/tOmeJJfi6TUT+GNDeLLtNg6FRMQ7K0OxQSy+sRES6YjmqpVjqzUOBQcB6e2RQCFytqgOA4cBfbFmYTyg/YhgDTBORi7BkXwZjjW4Gish5gfc+Tp1TgL+r6unAUeBau/xdrF34QwGf/QzFwATgU1Xtr6qf2nV7YwmUDgaetWV0KjIM+DHguKr7AhSr6nnAG1hO7X6gL3CriLSx62wAzgpyH4OhQTEOytBcKR1FlTqopQHHP9h1BHhBRNYDX2FJ/yer6hqgnYh0FJFfAEdUdQ+W9thFWOrMq7GcxSkV7ltdnV2qutb+/CNWnpx4IE5VS/v00XGea55aeXoOY+XmSg5SJ1GtvD+lVLpvwLlSLbQUYKNauYOKsPTkugCoqg8oLtVhMxgaC7PwaWiulK5DnYE1AtgLPApkA+/YdW7EUloeqKpesRTKI+1z04HrgPZYIyqwHNqfVPXNau4btI5YeXiKAop8QBTBUxBUR8U2gv0Pl4iIS4+lRQl234rt+SvU81doOwJrxGkwNBpmBGVornwPXA5kqqpPVTOBeKxpvqV2ndbAIds5DcdKrFbKJ1hqy9dhOSuwlO1vFyuvDiLSSay8OIGcSJ0y1MrcnCMipakYxgaczsFK0V1TtmAJe9YL9lRfuqp666tNg+FEMA7K0FxJj6vpRAAAAOBJREFUwYreW1ahLMueHgP4EBgkIquwRlObSyuq6kYs55CmqvvtsoVYU3BLRSQFy3GVcyAnUicIdwBTRWQp1ogqyy7/GisoIjBI4kSYB/yqBvWPx3Dgi3psz2A4IYyaucEQYkQkVq2kdojIE0AHVX2wDu11wMpTdWE99W8G8KSqbqmP9gyGE8WsQRkMoecyEXkS6//xZ6w8U7VGVffbYeWt6roXyg53n2WckyEUmBGUwWAwGByJWYMyGAwGgyMxDspgMBgMjsQ4KIPBYDA4EuOgDAaDweBIjIMyGAwGgyP5f3LYTejT4DtJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Another example\n",
    "mask = ((wl_ / 1000) > 2.10580) & ((wl_ / 1000) < 2.1059)\n",
    "\n",
    "wl, flux = wl_[mask], flux_[mask]\n",
    "\n",
    "flux = flux / max(flux)\n",
    "\n",
    "f_slope = slope(wl, flux)\n",
    "f_grad = slope_grad(wl, flux)\n",
    "f_adj, new_wl, new_flux = slope_adjusted(wl, flux)\n",
    "\n",
    "# plt.figure(figsize=(15,10))\n",
    "ax1 = plt.subplot(211)\n",
    "plt.plot(wl, flux, \"o-\", label=\"Spectrum\")\n",
    "# plt.plot(new_wl, new_flux, \"+--\", label=\"Diff adjusted\")\n",
    "plt.ylabel(\"Normalized Flux\")\n",
    "\n",
    "plt.legend()\n",
    "ax2 = plt.subplot(212, sharex=ax1)\n",
    "plt.plot(wl[:-1], f_slope, \"s--\", label=\"FFD\")\n",
    "plt.plot(new_wl, f_adj, \"o-.\", label=\"FFD(shifted)\")\n",
    "plt.plot(wl, f_grad, \"*--\", label=\"Numpy\")\n",
    "plt.ylim([-2, 3])\n",
    "plt.xlabel(\"Wavelength (nm)\")\n",
    "plt.ylabel(\"Gradient\")\n",
    "ax2.ticklabel_format(useOffset=False)\n",
    "ax1.tick_params(labelbottom=False)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"spectra_gradient_example2.pdf\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gradient method is ~7 times slower than the difference.\n",
    "# %timeit np.diff(flux_) / np.diff(wl_)\n",
    "# %timeit np.gradient(wl_, flux_, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions to calcualte wis, and q and RV from the slopes.\n",
    "def slope_wis(wavelength, flux):\n",
    "    derivf_over_lambda = slope(wavelength, flux)\n",
    "    wis = np.sqrt(\n",
    "        np.nansum(wavelength[:-1] ** 2.0 * derivf_over_lambda ** 2.0 / flux[:-1])\n",
    "    )\n",
    "    return wis\n",
    "\n",
    "\n",
    "def grad_wis(wavelength, flux):\n",
    "    derivf_over_lambda = slope_grad(wavelength, flux)\n",
    "\n",
    "    wis = np.sqrt(np.nansum(wavelength ** 2.0 * derivf_over_lambda ** 2.0 / flux))\n",
    "    return wis\n",
    "\n",
    "\n",
    "def slope_adjusted_wis(wavelength, flux):\n",
    "    derivf_over_lambda, wav_new, flux_new = slope_adjusted(wavelength, flux)\n",
    "\n",
    "    wis = np.sqrt(np.nansum(wav_new ** 2.0 * derivf_over_lambda ** 2.0 / flux_new))\n",
    "    return wis\n",
    "\n",
    "\n",
    "def q(wavelength, flux):\n",
    "    \"Quality\"\n",
    "    return slope_wis(wavelength, flux) / np.sqrt(np.nansum(flux))\n",
    "\n",
    "\n",
    "def grad_q(wavelength, flux):\n",
    "    \"Quality\"\n",
    "    return grad_wis(wavelength, flux) / np.sqrt(np.nansum(flux))\n",
    "\n",
    "\n",
    "from astropy.constants import c\n",
    "\n",
    "\n",
    "def slope_rv(wavelength, flux):\n",
    "    return c.value / slope_wis(wavelength, flux)\n",
    "\n",
    "\n",
    "def grad_rv(wavelength, flux):\n",
    "    return c.value / grad_wis(wavelength, flux)\n",
    "\n",
    "\n",
    "def slope_adjusted_rv(wavelength, flux):\n",
    "    return c.value / slope_adjusted_wis(wavelength, flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Numerical Gradient Effect on Bands\n",
    "Need to calcualte the relative difference on the same wavelength and flux sections. They are all normalized differently due to the maximum in the given range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Band    wl_min wl_max   dy/dx   gradient Q(dy/dx) Q(grad) Q(frac) RV(dy/dx) RV_adj RV(grad) RV(frac_grad) RV(frac_adj)\n",
      "VIS      0.38  0.78  1.86e+07  1.78e+07  40206.0  38342.5 -0.046     16.1    16.2    16.9   0.049   0.006\n",
      "CARM_VIS 0.52  0.96  1.43e+07  1.36e+07  24999.0  23763.9 -0.049     20.9    21.0    22.0   0.052   0.003\n",
      "Z        0.83  0.93  3.90e+06  3.80e+06  12713.0  12404.2 -0.024     76.9    77.0    78.8   0.025   0.001\n",
      "Y        1.00  1.10  3.83e+06  3.58e+06  17489.7  16341.9 -0.066     78.3    78.5    83.8   0.070   0.002\n",
      "J        1.17  1.33  2.01e+06  1.92e+06   7230.6   6902.7 -0.045    149.3   149.4   156.4   0.047   0.001\n",
      "H        1.50  1.75  2.51e+06  2.45e+06   9131.7   8909.1 -0.024    119.4   119.5   122.3   0.025   0.001\n",
      "K        2.07  2.35  1.95e+06  1.90e+06   8069.9   7852.8 -0.027    153.4   153.7   157.7   0.028   0.002\n",
      "CARM_NIR 0.96  1.71  6.51e+06  6.24e+06  11235.1  10778.1 -0.041     46.1    46.2    48.0   0.042   0.001\n",
      "NIR      0.83  2.35  8.13e+06  7.85e+06  10766.4  10392.9 -0.035     36.9    36.9    38.2   0.036   0.001\n"
     ]
    }
   ],
   "source": [
    "import eniric\n",
    "from eniric.utilities import band_limits\n",
    "\n",
    "wl_, flux_ = load_aces_spectrum([3900, 4.5, 0.0, 0])  # M0 star\n",
    "wl_ = wl_ * 1000\n",
    "\n",
    "bands = [\"VIS\", \"CARMENES_VIS\", \"Z\", \"Y\", \"J\", \"H\", \"K\", \"CARMENES_NIR\", \"NIR\"]\n",
    "\n",
    "print(\n",
    "    \"Band    wl_min wl_max   dy/dx   gradient Q(dy/dx) Q(grad)\"\n",
    "    \" Q(frac) RV(dy/dx) RV_adj RV(grad) RV(frac_grad) RV(frac_adj)\"\n",
    ")\n",
    "\n",
    "for band in bands:\n",
    "    wl_min, wl_max = band_limits(band)\n",
    "    mask = ((wl_ / 1000) > wl_min) & ((wl_ / 1000) < wl_max)\n",
    "    wl, flux = wl_[mask], flux_[mask]\n",
    "    flux = flux / max(flux)\n",
    "\n",
    "    s = slope_wis(wl, flux)\n",
    "    sa = slope_adjusted_wis(wl, flux)\n",
    "    g = grad_wis(wl, flux)\n",
    "\n",
    "    # Quality\n",
    "    qs = q(wl, flux)\n",
    "    qg = grad_q(wl, flux)\n",
    "\n",
    "    qfraction = (qg - qs) / qs\n",
    "\n",
    "    # RV_rms\n",
    "    rvs = slope_rv(wl, flux)\n",
    "    rvg = grad_rv(wl, flux)\n",
    "    rvsa = slope_adjusted_rv(wl, flux)\n",
    "\n",
    "    rv_fraction = (rvg - rvs) / rvs\n",
    "    adj_rv_fraction = (rvsa - rvs) / rvs\n",
    "\n",
    "    if \"CARMENES\" in band:\n",
    "        band = \"CARM\" + band[-4:]\n",
    "    print(\n",
    "        (\n",
    "            \"{0:8} {1:1.02f}  {2:1.02f}  {3:7.2e}  {4:7.2e}  {5:7.1f}  {6:7.1f} \"\n",
    "            \"{7:-6.03f}  {8:7.01f} {9:7.01f} {10:7.01f}  {11:-6.03f}  {12:-6.03f}\"\n",
    "        ).format(\n",
    "            band,\n",
    "            wl_min,\n",
    "            wl_max,\n",
    "            s,\n",
    "            g,\n",
    "            qs,\n",
    "            qg,\n",
    "            qfraction,\n",
    "            rvs,\n",
    "            rvsa,\n",
    "            rvg,\n",
    "            rv_fraction,\n",
    "            adj_rv_fraction,\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Changing the RV precision calculation from using dy/dx to using `np.gradient` gives a $2.5-7$% **decrease in RV precision** (Increase in RV error) alone.\n",
    "The gradient function returns the exact same size section of wavelength. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adjusting the wavelength values to the center of each diff changes the RV precsion by only ~0.1% (last column)\n",
    "\n",
    "Gradient is larger change by ~30x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### The current version of eniric now uses np.gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}