{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to load and use Cotrending Basis Vectors for Kepler, K2 and TESS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cotrending Basis Vectors (CBVs) are generated in the PDC component of the Kepler/K2/TESS pipeline and are used to remove systematic trends in light curves. They are built from the most common systematic trends observed in each PDC Unit of Work (Quarter for Kepler, Campaign for K2 and Sector for TESS). Each Kepler and K2 module output and each CCD in TESS has its own set of CBVs. You can read an introduction to the CBVs in [Demystifying Kepler Data](https://arxiv.org/pdf/1207.3093.pdf) or to greater detail in the [Kepler Data Processing Handbook](https://archive.stsci.edu/kepler/manuals/KSCI-19081-003-KDPH.pdf). The same basic method to generate CBVs is used for all three missions.\n", "\n", "This tutorial provides examples of how to load CBVs from MAST and set them up in a design matrix for use to remove systematic trends. Please consult the [DesignMatrix page](https://docs.lightkurve.org/reference/api/lightkurve.correctors.DesignMatrix.html#lightkurve.correctors.DesignMatrix) in the API docs for the full details on that class. A convenient tool has been created called [CBVCorrector](https://docs.lightkurve.org/reference/api/lightkurve.correctors.CBVCorrector.html?highlight=cbvcorrector) which will utilize the CBVs for the custom removal of systematic trends. See the [CBVCorrector HOWTO page](https://docs.lightkurve.org/tutorials/2-creating-light-curves/2-3-how-to-use-cbvcorrector.html) for example usage." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Cotrending Basis Vector Types" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are three basic types of CBVs: \n", "- **Single-Scale** contains all systematic trends combined in a single set of basis vectors. \n", "- **Multi-Scale** contains systematic trends in specific wavelet-based band passes. There are usually three sets of multi-scale basis vectors in three bands.\n", "- **Spike** contains only short impulsive spike systematics.\n", "\n", "There are two different correction methods in PDC: Single-Scale and Multi-Scale. Single-Scale performs the correction in a single bandpass. Multi-Scale performs the correction in three separate wavelet-based bandpasses. Both corrections are performed in PDC but we can only export a single PDC light curve for each target. So, PDC must choose which of the two on a per-target basis. Generally speaking, single-scale performs better at preserving longer period signals. But at periods close to transiting planet durations multi-scale performs better at preserving signals. PDC therefore mostly chooses multi-scale for use within the planet finding pipeline and for the archive. You can find in the light curve FITS header which PDC method was chosen (keyword “PDCMETHD”). Additionally, a seperate correction is alway performed to remove short impulsive systematic spikes.\n", "\n", "For an individual's research needs, the mission supplied PDC lightcurves might not be ideal and so the CBVs are provided to the user to perform their own correction. All three CBV types are provided at MAST for TESS, however only Single-Scale is provided at MAST for Kepler and K2. For TESS, Cotrending Basis Vectors are currently only supplied at a 2-minute cadence. For Kepler/K2 only for the 30-minute target cadence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Obtaining the CBVs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Two tools are available to automatically download the CBVs relevent to your target of study: `load_tess_cbvs` and `load_kepler_cbvs`. Here is an example loading in the TESS Multi-Scale Band 2 CBVs for TESS target TIC 99180739, which happens to be on camera 1 and CCD 1. Note that you do not need to have a lightcurve object already loaded. One can directly request whichever CBVs one desires." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:20.721272Z", "iopub.status.busy": "2023-11-03T14:13:20.721086Z", "iopub.status.idle": "2023-11-03T14:13:26.138120Z", "shell.execute_reply": "2023-11-03T14:13:26.137657Z" } }, "outputs": [], "source": [ "from lightkurve import search_lightcurve\n", "from lightkurve.correctors import load_tess_cbvs, load_kepler_cbvs\n", "import numpy as np\n", "lc = search_lightcurve('TIC 99180739', author='SPOC', sector=10).download(flux_column='sap_flux')\n", "cbvs = load_tess_cbvs(sector=lc.sector, camera=lc.camera, ccd=lc.ccd, cbv_type='MultiScale', band=2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This TessCotrendingBasisVectors object contains 8 Multi-Scale band 2 CBVs for Sector 10, Camera 1, CCD 1. A CBV object contains only one type of CBVs. To obtain all types you must request each seperately. The basis vectors themselves are containted in an astropy.timeseries.TimeSeries Table:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:26.140528Z", "iopub.status.busy": "2023-11-03T14:13:26.140292Z", "iopub.status.idle": "2023-11-03T14:13:26.145019Z", "shell.execute_reply": "2023-11-03T14:13:26.144676Z" } }, "outputs": [ { "data": { "text/html": [ "
TessCotrendingBasisVectors length=3\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
timeCADENCENOGAPVECTOR_1VECTOR_2VECTOR_3VECTOR_4VECTOR_5VECTOR_6VECTOR_7VECTOR_8
Timeint32boolfloat32float32float32float32float32float32float32float32
1569.43311737673732462240-0.000260.00163-0.000160.000140.00005-0.00150-0.00004-0.00007
1569.4345062650322462250-0.000360.00158-0.00025-0.000070.00006-0.00152-0.00006-0.00007
1569.4358951533282462260-0.000450.00152-0.00034-0.000280.00007-0.00153-0.00008-0.00008
" ], "text/plain": [ "TESS CBVs, Sector.Camera.CCD : 10.1.1, CBVType.Band: MultiScale.2, nCBVs : 8" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Show the head of the CBV table\n", "cbvs[1:4]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All data contained in the MAST CBV FITS files is downloaded. The `time`, `CADENCENO` and `GAP` columns can be used to select which cadences you desire. Extra information is in the `cbvs.meta` dict. The CBVS are interpolated accross the gaps and so should be used in gapped cadences with extreme caution. Below we will plot the first 4 CBVs. Note: _CBVs use 1-based indexing!_" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:26.147073Z", "iopub.status.busy": "2023-11-03T14:13:26.146926Z", "iopub.status.idle": "2023-11-03T14:13:26.373725Z", "shell.execute_reply": "2023-11-03T14:13:26.373366Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqUAAAGQCAYAAACXhxmxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACR0ElEQVR4nOzdd3RT9f/H8WdWk6Z77xZa9ipLhmwFRXGLgKjgxq24FbfinqD4VQE3Mhw4cDEEQQVk71FG6d4rzc69vz8q+VlZBdqmhffjnB5oc3PvJ7m5977uZ0WjqqqKEEIIIYQQPqT1dQGEEEIIIYSQUCqEEEIIIXxOQqkQQgghhPA5CaVCCCGEEMLnJJQKIYQQQgifk1AqhBBCCCF8TkKpEEIIIYTwOQmlQgghhBDC5ySUCiGEEEIIn5NQKk5bTz75JGvWrPF1MYQQzdhdd93FokWLfF0MZs2axXPPPXfExwsLC7niiiuorq5u8LLccMMN/PXXXw2y7qlTp/LLL780yLqF7+l9XQBR/6644grv/51OJ1qtFr2+Zld36NCBp59+mgsvvBA/Pz+02v+/L+natSuTJk0CYMmSJXz99dcUFhai1+tJTEzkxhtvpE2bNnV6/L88Hg/fffcdixYtoqCgALPZTMuWLbnssstIT09n0aJFTJ06FT8/PwD8/f3p06cPN910E5WVlVx//fW8+uqrtG7dutZ6f/nlF+bOncsHH3xQ67Ucy6ZNm6ioqKBnz57ev33zzTf8/PPPlJWVYTAYaNmyJXfeeScxMTF1Xu/hvPHGGwQGBnLTTTed1Hr+a9GiRSxYsIDs7GxMJhMJCQlccMEF9O/fv16305i2bt3K3Llz2blzJwBRUVEMHjyYiy66CIPBABz9dW/evJlHH30Uk8mERqPBaDSSmprKeeedR58+fU64XJmZmcyYMYOMjAyqqqr44osvCAwMrLXMTz/9xNy5c6mqqqJz587ceeedhIeHH3Z9mzZtYvbs2ezZsweNRsPs2bOPuv3ly5fz3XffsXfvXhISEpgyZcpxv4Zjvbf/PifodDpatmzJTTfdRGpqKk899RSxsbHccssttdZptVq55pprGD16NPPmzQNAURScTicmk8m73FNPPUXHjh2Pu8zHa9asWcyZM8d7HgkMDOTss8/mqquuQqPRNPj2T8SFF16I0Wjkk08+wWw2e//+9NNPs2bNGh599FH69u17Qut96623SE1NBSA6Otq7j6Bm33300UesWrUKq9VKYGAg7du358EHHzz5F3USvv76a3777TfvdWLAgAGMGzfOe/yPGjWKBx98kLPOOsv7N3HqkFB6Cvr3ieeRRx6hT58+XHzxxYcs98orr3hPWP+2detW3n//fR5//HE6dOiAw+Fgy5Yt3hPAsR4/nNdee439+/dz66230q5dOwDWr1/PH3/8QXp6OgApKSnei21JSQlPPfUUX3/9NaNHj6ZHjx4sXLjwkFC6cOFChg4delyBFGDBggUMHTrU+/tvv/3GDz/8wBNPPEFKSgoWi4X169c3iQuZ2+323lQc9NFHH7Fs2TJuvfVWunTpgsFgYNu2bfz8889NIpR6PB50Ot1xPWf16tW8+uqrXHXVVdx7772EhISQlZXFl19+SVlZGdHR0XV63QEBAd6QZ7VaWb16NVOmTOHAgQOMGjXqhF6PTqejf//+XHDBBTz77LOHPL5x40Y+/vhjnn76aVJSUnjvvfd47bXXmDx58mHXZzKZGDp0KIMHD2bmzJnH3H5QUBAXXXQRubm5/PHHH8dd/rq8t/D/5wS3282sWbN48cUXef/99xk2bBhvv/02N9xwQ63jfNmyZYSHh3PFFVd439vNmzczefLkYwbthnLGGWfw2GOPAZCbm8sjjzxCYmIigwcP9kl56iIyMpLly5dz7rnnAlBaWsrOnTsJDQ1tsG1Onz6d0tJS3nzzTcLCwiguLubvv/9usO3VlaIo3HnnnaSmplJeXs7kyZP54osvGDduHAAxMTHEx8fzxx9/NOl9Kk6MhFJxiJ07d5KWluat2TCZTLVqFI/1+H9t3ryZv/76i2nTphEXF+f9e69evejVq9dhnxMREUH37t05cOAAAMOGDeOtt97ixhtv9NaCZGVlsXv3bu+d/dKlS5k1axbl5eX4+/tz3nnnMWbMmEPW7Xa7WbduHVdddVWt15Senk5KSgpQU8MyYMCAWs/bsGEDn3zyCbm5uYSHhzN+/Hh69+4N1JxIf/jhB3788UdKS0sJDQ1lwoQJ5OTksGzZMjQaDb/++itRUVFMmzYNq9XKzJkzWb16NQC9e/fmhhtuwGQyUVBQwI033sjdd9/N3LlzsdlsfPrpp95y5Ofn88033/D888/Xqn3q3LkznTt3Bmqa6qZOncrevXvxeDy0b9+eW265xVvr+8Ybb6DT6aiurmbdunVERUXx0EMPsW3bNubMmYPL5WLs2LGMGDHCu/7ff/+defPmUVRURFxcHDfffDPt27cHam5+2rRpw969e9m+fTsPPvgger2eTz/9lNzcXIxGI3369OGGG27AaDQesk9UVeX999/n8ssvr3UDlZSUxMSJE+v8uv/LbDYzePBgtFotb775Jueddx5BQUGHXfZoEhMTSUxMpKCg4LCPL1q0iMGDB9O2bVsAxo0bx7XXXkt+fj6xsbGHLN+mTRvatGnD5s2b67T9rl27erdzvOry3v6XXq9n8ODBzJs3D5fLRe/evXn33XdZuXJlreNi8eLFDB069Ig3b/v27eOhhx7i448/xt/fH6i54bzxxhuZPn06ubm5TJ48mauuuop58+ahqirDhw9n7Nix3nUe7bg7lvj4eNq3b+89jwB8+OGHLF++HIvFQmRkJGPHjvXe0BwM1Ndddx2zZ8/G4XAwbNgwrrvuOu/zf/jhB7766iscDgfDhw+vUzmOZejQoSxatMgbSpcsWUL//v1Zu3atd5lZs2axd+9eb+AGGDNmDJMmTTrk83/vvfcC8MADD6DVarniiisYNGgQN954o7eWf+fOnVx++eWEhYUBNcH4vPPO867jSOe0Hj16sG7dujof23B8+3DkyJHe/0dGRjJkyBD+/PPPWsukp6ezatUqCaWnIOlTKg7Rvn17tm3bxscff8ymTZuwWq3H9fh/rV+/njZt2tQKpMdSVFTEunXr6NChA1ATYA0GQ61+SosWLaJr165ER0djt9t58803ueuuu5g7dy7vvPMO3bt3P+y6c3NzcTgcJCQk1HpNK1asYM6cOWzbtg2n01nrOfv27ePFF19k/PjxzJo1i9tvv53XX3+d7OxsoOZC9d1333H//fczZ84cnnvuOaKjo7nooosYNGgQ559/PvPmzWPatGkAfPDBB+Tl5fH2228zdepUsrOzmT59eq1trlq1itdff50PPvig1t83bNhAeHj4UZtDVVXl4osv5sMPP2TmzJkYjUbefvvtWsusWLGCiy++mNmzZ9O6dWuee+458vPz+eCDD3jwwQeZPn06ZWVlAKxZs4aZM2dy9913M2vWLK644gqeffZZKisrvetbvHgx11xzDfPmzSM9PR2j0cgdd9zBrFmzePnll9m8eTPz588/4j4pKChg4MCBR3xNdXndR9K3b188Ho+36fq/li5dyp133nnc6z1o//79tGzZ0vt7WFgYoaGh7N+//4TXWV/q8t7+l9PpZMmSJbRq1QqDwYBer2fIkCEsXLjQu8yBAwfYvXt3rRaH/2rZsiUJCQm1aneXLFlC165diYiIAMBms7Fnzx7ef/99nn/+eRYtWsSSJUuAYx93x5KVlcX27du9N08Hy/T666/zxRdfMGbMGF5//XXy8/O9j9tsNrKysnjvvfd46aWXWLBggffmYePGjXz66afeoA01XTuOZsyYMWzduvWoy3Tt2pXi4mKysrKAmnPb0d7XY3n99deBmprvefPmHbaFoH379syePZuff/6Zffv2oapqrcePdE4DjuvYPtl9uGXLFlq0aFHrb0lJSezbt69OzxfNi4TS09jDDz/MmDFjvD9ffPEFUHOyeuqpp8jNzeXll1/mqquu4oUXXqCioqJOj/9XRUWF9wJ0NJmZmYwZM4bRo0dz/fXX4+/vz5AhQ4Ca5tOzzjrLW1Pk8Xj47bffGDZsmPf5Op2OrKwsb/+oI/VvtVgsGI3GWs3LgwYN4u6772bHjh0888wzjB07lqlTp2K32wH4+eefOfvss0lPT0er1dKxY0fOOOMMVqxYAdT0Jxw7diytWrVCo9EQHR1NUlLSYbevKApLly5l3LhxBAcHExISwrhx41iyZAmKoniXu/LKKwkMDKzVN6+u72dMTAw9e/bEz88Ps9nMqFGj2Lp1a631n3HGGXTo0AGdTseAAQMoLCxk7NixGAwG0tPTCQgI8F5wFyxYwGWXXUarVq3QarWceeaZJCYm1hooNmjQINq0aePty9mxY0fS0tLQ6XTExsYyfPhwtmzZctjyHvzsHO111fVzdDgGg4Hg4GAsFsthHx88eDBTp049oXUD2O32Q/qYBgQEYLPZTnid9aUu7+1BB88Jo0aN4ueff2b8+PHex4YNG8bGjRspKioCarrOdO/e/ZjrHTZsGIsXL/b+frB29SBFUbj22msxmUwkJSUxYsQIfvvtN+DYx93hrFmzxvsabrvtNtq2beutaYaafR0aGopOp2PgwIEkJiayY8cO7+OqqnL11Vfj5+dHUlIS7du3JyMjA6jprjBo0CDatWuHwWBg7Nixhxyf/zV79uxj3khptVqGDBnC4sWL2b59Ozqd7ojnr/oyYcIEzjvvPBYvXsx9993HNddcUytYHu2cdjzH9onsw4N++eUXtm/ffkioNpvNRzyWRfMmzfensRdffPGwfUqhpnnkYF/Pffv28eabb/L+++/zwAMP1OnxfwsODq7TXfG/+5TabDY+++wznnzySV555RWgponrjjvuoLCwkP379+N2u71NQCaTiSeeeIJvvvmGjz76iJSUFK6++mq6dOlyyHYCAwNxOByH9Hvs168f/fr1Q1VVtm3bxmuvvcbcuXMZN24chYWFbNq0qdbF1ePxeAcmFBYWEh8ff8zXCDUhwe12e2sdAGJjY3G5XLVqHqOiog77/ODgYEpKSo65jffff59t27Z5R9u6XC5sNhsBAQEAtfqrGY1G/P39azW/GY1Gb6gqKCjgk08+YdasWd7H3W53rXL8t7y7du3ik08+ITMz0/t+JyYmHvE1QU3T7pFq1Ovyuo/k4Hv73+BYX0wm0yGjmq1Wq7fJ2pfq8t4edPCc4PF42LJlCy+88AIvvvgiLVq0IDk5mdatW7NkyRJGjhzJ0qVLufXWW4+5/YEDBzJz5kzy8/MpLy+nsrKyVrcdPz+/Wp/F6Oho734+1nF3OD179vQ2cVdVVTFt2jTefPNN77lp/vz5/Prrr5SUlKDRaLDZbLWOO7PZXCto/vs4KCkpqdVUrtfrvc3fJ2vo0KE88sgjlJeXn1QtaV0ZDAYuvfRSLr30UlwuFytWrGDKlCkkJyfTvXv3o57TjufYPpF9CDWtF5999hnPPPPMIQMGD1Y8iFOPhFJxTC1btmTo0KFHnIbjWI93796d+fPnH7F/3eH4+/tzzjnn8N1331FZWUlwcDBJSUm0bduWxYsXs3fvXoYMGVJr0MXBoOx2u/nxxx+9HeT/OwgqPj4eo9FITk4OycnJh2xbo9HQsWNH+vXr521+jYyM5MILL+Taa689bHmjo6PJzc31DuL6t/9uPyQkBL1eT2FhofeCVlBQ4K3NO1gTdaR+el27duV///sf27Zt83Zv+K+PP/4Yh8PBm2++SUhICHv37uXuu+8+pImurqKiorjwwgtr9Tn7r/+W99VXX+Xss8/msccew2Qy8e2339a6MP1bQkIC0dHR/P7774wePfqwy9TldR/JX3/9hV6vP+z+qQ8tWrSo1ZxYXl5OWVnZIc2OvlCX9/a/dDod6enpxMXFsX79eu/rGDZsGF999RXJycmoqnrEPuH/FhgYSN++fVmyZAmlpaUMHjy41nHrdDopLy/3BtOioiJv7euxjrtjCQoKYsiQIbz66qtAzSDNL774gsmTJ5OamopWq+Wuu+6q83ERERFBYWGh93e32+3t4nKy4uPjiY2NZdmyZYcd/GYymXA4HN7f7Xb7UbtOHc8gTYPBwJAhQ/j222/JzMyke/fuRz2nHc+xfSL7cOnSpXzwwQc888wztbrFHJSVlXXYv4vmT5rvxSH++usvlixZ4m32y8/PZ9myZd6T07Ee/6/OnTvTt29fnnvuObZu3YrL5cLtdrN27Vrefffdwz7H4XCwaNEiwsPDaw1MOeecc/jll19Ys2YN55xzjvfvZWVl/PXXX1itVnQ6HWaz+Yijv/V6Pd26das1yGTRokWsXLnS2ySUmZnJqlWrvH3Rhg8fzuLFi9m0aRMejweXy8WOHTu8fcCGDx/O7Nmz2bt3L6qqUlhY6H0sNDSU/Px874VPq9UyaNAgPv30U6qqqqisrOTTTz9lyJAhdZpFIC4ujksvvZRXX32Vv//+G7vdjsfjYevWrbz22mtATU2C0WgkICCAyspKb9eMEzVixAi+/vprMjIyUFUVu93Ohg0bKC4uPuJzDtZmmEwmsrKy+Omnn464rEajYcKECXz11Vd8//333pqrnJwcpkyZQmFhYZ1e93/ZbDaWL1/Oe++9x+jRo0+4dkVVVZxOJy6XC6ipeXU6nd59OnToUH777Td27dqF3W7nk08+oVOnTke8CTs4bdLB9TmdzkP6Mf+bx+PB6XTi8XgOKQvUDIJ55JFHDvvcury3h3u9W7Zs4cCBA97BfwADBgygrKyM6dOnM2TIkENmhTiSg034K1asqNXlBmqOh08++QSHw0F2djYLFizwDmA51nF3LNXV1Sxbtsz7Gmw2G1qtluDgYFRVZeHChcfsE/pvAwcOZNmyZezcuROXy8Xs2bO9XXzqwz333MMLL7xw2NrXtLQ0du7cSVZWFk6nk08++eSowTM0NJS8vLwjPv7FF1+wfft2b03nqlWryMrK8p7Hj3ZOO55j+3j34bJly3j//fd56qmnSEtLO+wymzZt4owzzjjiNkXzJTWlp7GDIzMPSkpK4vXXXycoKIgff/yRGTNm4HQ6CQoKonfv3t4pOY71+OHcd999fPvtt7zzzjsUFhbWmqf0oMzMTO8cq3q9nlatWvHkk0/WOvH279+f999/n9TU1FoXS1VV+e6773jrrbdQFIWEhAQefvjhI4a8ESNG8NFHH3lHlwcEBDB//nzeeustPB4PoaGhDBw40Fu+tLQ07r//fj777DOysrLQaDSkpqZy/fXXAzVzAiqKwksvvURpaSnh4eFMmDCBpKQkzjnnHF566SWuvPJKoqKimDp1KjfffDPTp0/ntttuA2pG3x9c1+HcdtttjBo1ynuxvvbaa0lISGDWrFlkZWVhMplITEzkwgsvBOCqq67ijTfe4MorryQiIoJLLrmElStXHnH9x9KrVy+cTidTp0711uq2bt36kHkr/+32229n+vTpfPTRR6SlpTFgwABWrVrlffydd97xLndwG08++SRz587l888/B2pqaIcMGeK9SB/rdUNNELniiitqzVN6xx13HHWux6VLlzJ37lzvQLT/Kiws5MYbb/T+fvCzPn36dGJiYkhPT2f8+PE8//zzWCwWOnXqxH333XfE9W/dupVHH33U+/jll18OwPfffw/A3Llz2bp1K08//TRQM2XZW2+9VWv56OhoZsyYAdTULv57MM9/1eW9hf8/J2g0GsLDw7n++utrDRg0m83079+fRYsW1bopPJbOnTuj1WqJiYk5pIbL39/fOyeqqqqce+65nHXWWcCxj7vD7be///7bex4xGAx07NjRuy+6d+9Ov379uPPOO721g0d73/6ra9eu3j70TqeT4cOH1zoPHc4VV1xR53la4+LijtjFIj09nXPPPZcHH3wQo9HIlVdeedTuIVdddRXvv/8+U6dO5fLLLz9koJtWq2XatGkUFBR4982dd97pfT+Odk471rH9b8e7Dz/55BOsVmut4+PgrCVQcyxmZ2c3ianvRP3TqCfanidEM/fEE09w0UUXHXU6KyGagzvuuIPnn3/e23+0KZo0aRJ9+/blggsu8P7N13Oaiubn7bffpnXr1t7ps8SpRUKpEEKIBrVjxw6eeOIJPvzwQ+9AO5BQKoSoTZrvhRBCNJgnn3ySnTt3cvPNN9cKpEII8V9SUyqEEEIIIXxORt8LIYQQQgifk1AqhBBCCCF8TkKpEEIIIYTwOQmlQgghhBDC5ySUCiGEEEIIn5NQKoQQQgghfE5CqRBCCCGE8DkJpUIIIYQQwucklAohhBBCCJ+r09eMKopCaWkp/v7+aDSahi6TEEIIIYQ4Raiqis1mIzw8HK32yPWhdQqlpaWlXHfddfVWOCGEEEIIcXr58MMPiYyMPOLjdQql/v7+3pWZzeb6KVkDURSFoqIioqKijprGRdMm+/HUIPux+ZN9eGqQ/XhqaK770Wq1ct1113nz5JHUKZQebLI3m83NIpT6+/tjNpub1Q4Ttcl+PDXIfmz+ZB+eGmQ/nhqa+348VhfQ5veKhBBCCCHEKUdCqRBCCCGE8DkJpUIIIYQQwucklAohhBBCCJ+TUCqEEEIIIXxOQqkQQgghhPA5CaVCCCGEEMLnJJQKIYQQQgifk1AqhBBCCCF8TkKpEEIIIYTwOQmlQogG89VXX/Hrr7/6uhhCCCGaAQmlQoh6kZGRQXV1tfd3RVFYvHgxy5YtQ1VVH5ZMCCFEcyChVIjTwOrVq1m3bl2DrX/FihW89957PPDAAzgcDgA2bdpEeno6qamp7Nmzp8G2LYQQ4tQgoVSI08D06dOZOXNmg63/yy+/5Omnn+aiiy7iyy+/BOCnn35i+PDhdOnSpUEDsRBCiFODhFIhTnEHDhygZcuWmEwmnE5nva+/srISf39/zGYz55xzDn/++SelpaXk5eWRkpJC+/bt2bhxY71vVwghxKlFQqkQp7jffvuNIUOGkJSURE5OTr2vf926dfTo0QMArVbL1VdfzU033cTw4cMBCAwMxGKxSL9SUSeqqspnRYjTlIRSIZohVVX54osv+Oyzz455AV+zZg09e/YkKSmJAwcO1HtZMjIyaNu2rff3vn378vbbb3Peeed5/5aUlERWVtYx12WxWCgsLKz3MormYceOHUyYMIFbbrmFbdu2+bo4QohGJqFUiGbozz//JD8/H4B33333iMvt37+fmJgY9Hp9nYPh8dqzZw+pqam1/hYXF4dGo/H+3rVrV9avX4+iKCxfvvyI3Qiee+45HnjgAWw2W72XUzRtdrudt956i1deeYXXXnuNKVOmUFlZ6etiNUlut5uPP/6YH374Abvd7uviCFFvJJQ2oAULFrBo0SJfF0OcYhRF4fPPP+fmm2/m6quvxuPxMHfu3FrL2O12ZsyYweTJkxk3bhwAsbGxFBQU1Ht5rFYrAQEBR13mYCh9++23Wbp0Ke+8884hy6xZs4a4uDguvvhili9fXu/lFE3bjBkzuOqqqwgJCSEwMJCbb76ZDz74wNfFanI8Hg/PPPMMERERANx///18+OGHPi6VEPVDQmkDWbx4MRs2bODrr7/2TpEjRH1YvHgxZ555pjcI3n777RQUFPDaa6/hdruprq7m4YcfpmXLlrz99tskJycDEBUVRVFRUb2WxW63YzQaj7lcZGSkt3b08ccfp7i4mP3793sfd7lczJw5k+uvv54hQ4awbNmyei2naJo8Hg/z589n8uTJlJSU0L9/f+9j3bt3p6SkpEFq95srVVV55ZVX6NWrFxdccAEXXHABU6dOpby8nDVr1vi6eEKcNL2vC3AqcDgcbN26lbVr17Jjxw5UVSUkJIRJkybx2WefsXnzZnr27OnrYopTgMfj4ZtvvuH111/3/k2r1XLnnXeyaNEi7rrrLjweD7feeitdu3at9VyTyVTvN0j79++nZcuWdVr2+eef9/7/1ltvZcqUKTzxxBP8+OOP/Pbbb1x99dUEBQUBUF1djdPpxM/Pr17LK5qW2bNn43K5uO6664iNjT3k8QkTJvC///2PyZMn+6B0Tc+MGTNITk7mggsu8P5No9Fw9dVX89Zbb8l1RjR7EkpPgN1uZ/Pmzaxbt44dO3ag1Wrp2LEjvXr14tprr8VgMHiXTU9PZ8OGDXKyEPViwYIFDB48GJPJdMhjQ4cO5cwzzwTAbDY3Snn27t17SH/SukhMTKR79+489dRTnHvuuUydOhW9/v9PRz169GDdunX06dOnPosrmpjVq1fz2muv1dr3/5aSkkJwcDCbN2+mc+fOjVy6pmXjxo3k5OTwxBNPHPJYVFQUTqeTiooKQkJCfFA6IeqHhNI68Hg8bN++nVWrVrFlyxb0ej2dOnWiX79+3HDDDUc8oQKkpaXx1VdfNWJpxanKbrfz008/MXXq1CMuU5cwqqpqrUFIJ2Pv3r3eqZ+O15gxYxgzZsxhHxs4cCCzZ8+WUHoKq6iowGw2H/X8CTW1pQ8//DBPPfXUYWtTTweqqvLhhx/y5JNPHvHYPffcc/n111+54oorGrl0TU9GRgazZs2irKwMrVaLTqfDZrORmJjIgw8+WG/nv+Zi4sSJdOzYkRtvvNHXRTkmCaVHkJeXx6pVq1i9ejUWi4X27dvTp08fxo8ff8yT6L8FBwdTVVXVgCUVTUFlZSWBgYFotQ3XTXvhwoUMHz78uD5//xUaGkp5eTlhYWH1UqbMzExvn9X6lJKSQlZWVr0GaNG0rFq1it69ex9zudDQUB577DGeeuopnn76aWJiYhqhdE3Lli1baNmy5VGP2wEDBnD//fczcuRIdu3aRWpqaq1Wu9PFihUrmD9/Pvfccw+JiYl4PB5cLhcmk4np06ezZMkSzj77bF8Xs9EUFhYSEhLChg0b8Hg86HQ6XxfpqCSU/qO6upq1a9eycuVKDhw4QFxcHL179+ahhx466eYQnU6H2+0+qTAhmiZVVZk5cyYbN24kISGBhx56qMG2tXDhQl5++eWTWkd0dDSFhYX1FkpdLleD9fvs2bMnv/zyywnXxIqm7Y8//uCOO+6o07KJiYneYPrGG28ctvvKqWzevHlMmDDhqMv4+fnRtWtXbrjhBhITEwkJCeG+++5rpBI2DS6Xi88++4ypU6d6A7lOp/MGsauvvpp7772XIUOGNGgFQlOybds2OnfuTHFxMZs3bz5krEFTc3rslcNQVZX9+/fz2WefMXHiRJ5++mny8/MZPXo0U6dOZdKkSQwdOrRe+uckJCQ0yDfpCN/76KOPUFWVKVOmUF5eTnFxcYNsZ+fOnSQnJ5/0xbg+R+DXdeT9ibr66qv54YcfyMvLa7BtCN9QFIWSkhKioqLq/JzExETGjRt31O4rvmaz2ep90v/i4mLsdjsJCQnHXPb666/n/fff55lnniE3N5fy8vJ6LUtTt3btWs4888wj1hCbTCbOPPNMfvvtt0Yume/s2LGD9u3bM2jQoGYxq8lpFUqdTicrV67k1Vdf5Y477mDu3LmkpaXx4osv8vLLLzNq1ChSUlLqvbkwNTWVffv21es6he/9+OOPFBcXc8MNNwA1A41+//33BtnWN998w6WXXnrS64mNjfVOun+ytm3bRvv27etlXYej1+uZOHEiM2fObLBtCN84eKE8Xn379kWv1/PTTz81QKlOTklJCffffz8ff/wxP//8c72td+bMmVx11VV1Xv5gi9yFF17IggUL6q0czcGff/5Jv379jrrMZZddxtdff33afJVtRkYGrVq1om3btuzatavJv+5TPpQWFRUxf/58Hn74YR566CF27drF5Zdfzttvv82DDz5I3759G7S2B6Bly5YSSk8Rqqqyfft2Jk+ezIYNG7j33nu9NzHdu3dnw4YN9b7N6upq8vPzSUtLO+l11ecE+ps2bWrwEdFpaWkUFhbicrkadDuicf3++++15iQ9HnfeeSebNm1i4sSJfPnll3g8Hu9jVqu10b/hSFEUVq1axaRJk7jvvvuYPHky3377ba1yHaSqKt999x3r168/5nrdbjczZ85Eq9WSnp5+3OXq378/f/zxR5MPIfWpLrOBmM1munXrxp9//tlIpfIdVVW9Xaw0Gg0tW7Zkz549vi7WUZ1ynRwVRWHnzp188803bNu2jfDwcPr168ekSZO8cyA2thYtWvD111/7ZNuifuTk5PDVV1+xY8cO2rRpw9ixYw+ZnzMsLKxBmstmz57NJZdcUi/rio2NrZfmcFVVWbdu3XHV4Jyobt26sWHDBs4444wG35ZoeKqqsnHjRm6++eYTer5er+ehhx7C4XCwYMECHnroIZ588kn+/PNPvvvuO1RVZdSoUQwePLh+C/4fFRUVTJ8+nczMTLp06cKLL75IaGgoAP369WP58uW1yqCqKu+++y5+fn4sX74cjUZzxP59hYWFPPXUU1x44YVcd911J1Q+vV5Peno6a9eupVWrVsydO5f09HSio6NPaH1NXV5eHrGxsXVq6Rw1ahSPPfYYZ5555ik1kFJRFF555RUuv/xyWrVqRWFhYa39PWrUKN5++21eeumlJvu6T4lQarfbWblyJb///jsFBQW0aNGCYcOGcfPNNzeJkWahoaHyHc5NVH5+Pu+88w7l5eWYzWZuv/1272jyg31Ef/nlF/bv38+1117LnXfeedSDOSoq6pATwclYtGgRe/fu5frrr6+X9QUGBlJdXX3S65k1axbdu3dvlNG9/fr1Y8GCBRJKmzGr1cqUKVMoKirCaDTSr1+/kx5oYjQaueyyy2jbti0PP/ww7du356233kJVVV577TW2bdvGhAkTGuQakJGRwauvvsott9xy2GB5ySWXMGnSJG8oLSoq4oMPPiAuLo7rrruO6upq7r//fp588kliY2NrzTJRVlbGk08+yaOPPkpSUtJJlfOSSy7hkUcewc/Pj2uuuYYZM2ZgNBqb/GCXE/HXX39552k+luDgYNq0acO6devo0aNHA5es8fz4448EBgby7rvv8tprr7F9+/Za3WSSk5M5//zzcTgcTXawYJMOpb/99hvR0dF07NjxkMfKyspYvnw5f/zxBy6Xi169enHTTTcRExPjDQVNaXSdRqNBUZQmVaaG9O9mA1/Zvn07u3btokWLFrRt2xZFUfjuu+/4888/0el0KIqCn58ft9xyC2lpaWRnZzNt2jQsFgsA4eHhREVF0aNHD26//fY6bbNz585s2rSJoUOHnnT5v/vuOzZs2MDTTz/dpO5qly9fzr59+5g0aVKjbK9Vq1ZkZGQ0yrZE/bPb7UyaNImrr76a9PR08vPz6zRop646duzIO++8U+tvDz/8MD/++CP33HMPd9xxB23btq237a1YsYJ58+bx/PPPEx4efthlAgMDSUlJYfny5Wzbto2MjAyuueYaunTpAkBAQACTJk3iueeeQ6/XoygKAQEBDBw4kO+//5777rvvpAMpQExMDJMnTyY0NBSj0UhcXBwffPABP/30E6NGjTrhLkHFxcV8//33WK1WunXrRq9evXw+u8zq1at57LHH6rz82LFjeeqpp+jWrVuTvy5XVVWxdetW/vzzTwIDA+nXr98hFR92u50FCxYwZcoUXnjhBbKzs9m+fTtDhgyptVxDtyCcrCYdSjt06MDrr7/OSy+9BNSMrPv999/Zt28fwcHB9O/fn8cee6xWs7yiKL4q7lHFxMRQUFBAXFycr4vSINxuNxkZGezYsYONGzdSXFyM0WjE4/Ewbtw4unXr1uBlODijQkZGBr/++isxMTF0796dv//+m88++wyAYcOG8frrr6PX6w+ZAzMxMbHWV2GeiC5duvDtt9+edChduHAh69ev5/HHH6/3E+bJTFFWWlrKF198wZtvvtloQVmj0RAXF0dubi7x8fHY7XYsFguRkZGNsv3TiaqqOBwOrFYrVqsVm81GdXU1eXl5+Pn54XA48Hg8qKrq/bFYLGRkZFBdXc0999xzSOB86aWXGD16tLdGKjExsVFey/nnn0+vXr2YMmUKZrOZ66+/vtaF/GBfy39/jvfv349Wq/W2luzYsYP33nvP+7jH4yElJYVXXnnlmDfct9xyC9OmTaNLly6Hnc4pMTGRt99+2/t7aWkpy5cv5/HHH6/X68TBdSmKgtlsZvLkyWRmZjJlyhTGjh173C0QGzdu5L333mP8+PGEhoayatUqPv30U3r37s3IkSMJDAyst7LXlcPhwOFwHNe2w8PDGTx4MB9//PEJd5FoSKWlpbz99tsUFxcTFBRE+/btufjiiyktLeW9997D398fp9OJy+XinHPOYd++fYwePRqDwcBFF13E/Pnz2bVrFzfddJOvX8px0ah16AVttVoZPXo0c+bMabSvLzzotddeIykpifXr1xMXF8eIESNITU094gVRUZQmWVM6Z84ckpKS6ty8cLIsFgsff/yxt4ZJr9fTv39/zj777BM+adhsNgoLC6mqqiInJ4cDBw5QVFTknWLo4Ai/Ll26eE/+5eXlzJgxg+LiYtq2bUtJSQkOh4OkpCSSk5Np27atdzLszMxMtm/fTkxMDG3btqWqquqo+1FVVbKzs9myZQubNm3iwIEDpKam0rZtW3r16uWTvlOKonDPPfcwZcqUE3q+zWZj8eLF/PHHHzz77LMNUvvwxhtvcMkll9T5O+sPUlWVRx99lPHjx9OuXbs6Pae+jsfffvuNoqIihg8fziOPPEJISAgBAQHccccdx5y2TVVVFixYQHl5OWPGjPF5jU5jUVWViooKiouLKS4upqKigoqKCiorK6moqKC8vNzbKnCQRqPBZDJhNpvx9/f3/rjdbmJiYvD390en06HRaLzn4ICAAFJTU3G5XDz77LO89dZb3vd4wYIFZGdnH3OOzYa2bds2Zs6cSYsWLbjgggtYsWKFdxDQyJEjSU9PZ+bMmdhsNvR6PSUlJej1ekJCQpgwYcJxTV3VVP33WLTb7TzyyCOMGDGCjRs3Ul5ezn333UdoaChutxun03nI9X7p0qX88MMPPPXUU7WuI4qisGLFCmbPns2ll17KsGHDTqqsK1eu5KOPPsJkMvHQQw8dM6QvXbqUvLw8rrzyyuPajqqqvPPOO+zbt4/w8HDcbrf3hn38+PG0aNHiJF7FiXO5XNxzzz3cc889tG7dutZjB/ej2Wz27p/58+fj8XgYPXo0UPO67rzzTlJSUnjggQcavfyHU9cc2eRDqcvlYsGCBaSnp9fpItpUQ+nWrVtZunQpt99+O3a7vd77c2RkZPDbb7/Rvn17XC4Xc+bM4YYbbqBnz55oNBpsNhtLlixhyZIlaLVaEhMTCQ8PJywsjNTUVNq0aeOtQTs40lmj0VBVVcXq1av5/fffURSFxMREAgMDiY+PJzk5mejoaMLDw495oc/Ly6OgoIDIyEgMBgNZWVns37+fXbt2UVhYiKqqJCcn06FDB/Lz89m+fTtVVVVHfZ88Hg9JSUl06tSJzp07k5CQ0CSaue+77z6ef/7545rVoby8nBdffBFFUejevTsjR45ssPB0cDqd884776jLbdiwgXnz5mGxWNDpdBiNRjp16nRcg5vq63h0OBzceeed6PV67rrrLtq1a8emTZt49913ue222446C8D8+fPJzs6mZcuWrFq1qsl1hzheHo+HsrIyb9j8909JSQlut9u7bEhICJGRkURERBAWFkZwcDAhISHen4CAgGPul+PZh/Pnz8dutzNmzBiysrJ45ZVXvC0TTcHKlSv5448/6NGjBwMGDMDlcjF37lwyMjIYOXKkt3ndZrPhdrt9Nji2IRxuP1qtVj7//HPvVFtTpkxh4sSJvPXWW9hsNh5//HFvMPvrr7+YP38+kydPPuL+dLlcvPfee1gsFu69997D1iYvXbqUDz/8kJYtWzJx4sRDbip37tzJtGnTePnllykpKeH555+ndevWFBUVkZqayrXXXlvrc+h2u7n77rt54YUXCA4OPqH35mDLgF6vR6/XU1lZybPPPstDDz10xGDqcrn4/PPPWbVqFR07dmTw4MFkZ2dTXl7O8OHDvYPe6sLtdmOxWLzPmTp1Kh07duSss846ZNm6Ho+KotS6efS1UyaUHq+mGkpVVeW2224jPDycoqIihg4dyqhRo+pl3bm5uTz33HPcfPPNbN++HY/Hw+WXX46/v/9hl7fb7eTl5VFaWkpZWRm7d+9m165dABgMBu8JR1VVAgIC6N69OwMGDGjUE3RT3Y91MX36dHr27HnEwQRut5vt27cTFRVFbGysd9DExIkT67Xv25Hs37+fefPm1bqDLi0t5Z133qG6uprY2Fj27t1LWloaY8aMISYmBofDQVlZ2XF/93h97sfi4mIURalVA26xWHjwwQd58sknD/v1kzk5Obzyyiu89tpr6HQ6PvvsM7RaLWPHjj2psjQEp9NJSUnJYcNmWVmZt7lZq9USHh7uDZuRkZHen4iIiHoffHY8+/BgS8HEiRN55ZVXeOKJJ07b76tvauqyH/fu3cu7777L7bffTkBAAE899RRvvvkm2dnZvPHGG7zyyit1utlevnw5s2fPZtKkScTHx3v//sMPP7Bq1SqeeOIJdu3axdtvv01KSop3gHLHjh358ssva81k4Ha72bt3L3Fxcfzyyy9kZWUxceJEoOba99JLL3HZZZcxaNCgk3+T/qW4uJjHH3+c8ePH07t3b2+4y83N5cMPP6SgoIDhw4dz3nnnsXHjRlatWkVycjL+/v58/fXXdOzYkeHDh5OcnHzUYKgoCpMmTcLhcBAXF0eXLl1YvXo1jz/++BGXb47XRgmlTXCHFRUV4XA4SEhI4OGHH+bmm28+6bknXS4XEydO5JFHHqnXwQO+1pT347Hs3LmT+fPnc/bZZ3vnGezfvz99+/Zl/fr1/Pzzz3Tr1o38/HxKSkoICQlh4sSJRxw0Ud9UVeX222/3DoDIyMjgtdde8w6uyM/P9zbVnqzG2I+ZmZm8/vrrvPbaa7VqcNxuNw888AD33Xeftx+jqqo8+eSTDBw4kD59+uBwOAgLC2vQz1h1dfVhg2ZxcbG3+Vyj0WAwGLxh878/oaGhPptJ5Hj3YVZWFtOmTWPMmDEnNL+maBgnciwuXLiQ33//naKiIp577rnj6sedk5PD22+/jc1m8w70TU1N5Y477vB+lh0OB6WlpURHR7N79242bdrEsGHDjvo1yO+99x6hoaGkpKTw2Wef8fDDDzdYP+WKigrmzp3L+vXrOeecc3C73fz+++/ce++9R23aV1XV28J44MABwsPD6datG+edd94hoX7q1Km0aNGCCy+8kLVr17J27VrGjx9/xPDfXK+NEkqb+A4rKiri2Wef5Y033kCn07Fs2TJmz57N5Zdf7h0kU5fBKC+//DK9evVq8iPqjldz2Y9HMm3aNGw2GxMmTMBkMrFs2TI2bNhA27ZtGT58uM+bM3fs2MEnn3yCzWYjPDycO+6446gXghPVWPvx559/ZuXKlQwdOtTbb3LVqlWMGDHikP5tDoeDjz76iKysLAwGA+3bt6+3Vot/mzx5MsXFxZjN5sMGzcjISAIDA5tM89qRNPdjUdQ40f24adMmkpOTj6s5+r9UVcXtdtdLLb6qqkybNg273c7tt9/eKFMbuVwufvrpJzQaDeeff/5x3yCWlpby559/8sMPP3gHllVWVlJaWkpCQgK33HJLndfVXI9HCaXNYIfNnz8fq9VK27ZtmT17Ns8++yzTp09nx44d6HQ6bDYbF198MSNGjDjs82fNmkVFRQW33nprI5e84TWn/SiOrDH346ZNm9ixYwdhYWEEBQXRpk2bRqt9PpXJsXhqkP3oey6Xi82bN2MwGAgODsZoNPq0S1RjqmuObBq9z09TF198Mc888wzbtm3jmWeewWQycccdd3gf93g8vP/++0yaNIno6GhKSkowm83ExsaydetW2rVrd1x3WEKcyrp06eIdqCKEEE2NwWCge/fuvi5Gkyah1Ic0Gg1PPvnkER/X6XTceuutZGdn43K5CA8Px2azkZuby4UXXkhEREQjllYIIYQQouFIKG0G/t2JOyQkREazCiGEEOKU03w6JAghhBBCiFOWhFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nN7XBRBCCHFi3B4Fi9NDld1NleP/fywON06Pisuj4FZq/nX9869HUdFrteh1GvRaDTqNxvt/f4OOIKO+1k+AnxZFVX39UoUQpwEJpUII0USoqkqRxUl2hY3scjsFFgdFFidF1TX/VtjdtZbXazUEGg8NkoFGPUa9FoNWg8mgxaA1YNBpMOi06LQaPIqK++CPR8Gjqrg8KmVWFwfKbP8EWw9VDjeVdhdFldUY/A6gQYMKBBn1RAf60S46kA4xQXyxPoeschtzx/XA7CeXFSHEiZGzRwNRVZVfdxZh9tMxIDXC18URQjQBbo/CgXIbGcXV7CmxkllqI7vC5g2bGiAq0I/EEH8SQkzEBhnpFBtEVKCRqAA/gk16NBpNo5ZZURQKCwuJjo5Gq9WiqioWh4cCi4Ot+VX8vreECX1T2JxXyax1OdzYJ6VRyyeEOHVIKK1HlXYXa7MrWJVZxpKMYgalRbA5rwqPojK4VaSviyeEaES/7ixka34Ve0qsZJXb/mk215ASbiYtwkxaRAADUyNIDDER4m/wdXHrTKPREGTSE2TS0yoygIs7xQLQKjKAa2dvkFAqhDhhEkpPkMujsDmvktUHyll9oJz8KgdBRj09k0Lo2yKMuwa0xOynp8LmYsxnaxmUFtHoNRxCCN+ptLvpEhfMpZ3jSAr1R6c9tY//QKMem8uDR1FP+dcqhGgYEkrrQFVV9pVavQF0Z5EFnUZD57ggeieHMfn8dsQFmw773BB/A/1ahPPDtgIu7BjbyCUXQvjKyPR4Xxeh0bWLDmR3kYV2MUG+LooQohmSUHoYJdVO/s4qZ1VmGRtyK3G4FVIjzPRKDuXmPsm0iQpEexw1AfcMTOXimX8zrE0UJoOuAUsufGF/qZUlu4tJjw+mR1Kor4sjhM90TwhhXU6FhFIhxAk57UOp3eVhQ24lqzLL+DurnFKri3CzgV7JoQxvF81DZ7U66SAZaNTz0FlpnD99Fe2jgxjbPYF+LcPr6RWI+lJucwFg1Gsx6bW1ultU2Fx8tSkPFegYG4TV6WHhriLWZVeQEu7P0NZRvLJ0D6PS47msS1yjlHdfiZUym5POccEYdDLlsPC9bgkhzFx9gLHdE31dlFOWoqjMWH2A5XtLaRHuT9uoQJLD/OnfMly6iP2H1elmfU4luZV2YgKNdIgJJDLQ6OtiiaM4rUKpoqjsLLJ4m+H3lVoxGbR0SwihV1Io1/RMJNzs1yDbPqdtNMPaRLG/1MZtX2/incs6kxoR0CDbOp3kV9pZk12Bw+3B9c+8jNkVdlZlluFWVLQaDfEhJjSA06MQHWgkIcREXqWdclvNiOcKu4sSq5Nwfz90Wg0Ot4L1n75xPZNCKLO62Ftq5eruifjpNSzaVYRRr2Vklzgmn9fOW2t+aedYRkxfRf+W4UQH1e+JT1VVPv47my35lQQZ9azNriAq0I+4YBNP/LwTg05LvxbhXN0jgdh/dSVRVZWNuZU8/vNOwvwNfDSm63HV8gtxPNpGB7KzqNrXxThl7S6ycM+3WxnRPoY3Lu7I/lIrO4ssLNxVxItLMnjrkk5YHG7+91cmrSICuKN/i9Oude77rfl89HcWNpeC/z/X96RQf/7KLGPm6gMcKLdx/+A0zm8f4+uiisM4pUNpXqWd1QfKWXWgjC15VSiqStvoQHolhfLAkDRSwvwb9c5So9HQMsLM1Es78/CC7cwd17PRtt0cuRWVIouDYqubHYVVbCuwkFdpx+6qmVex0OIgKsBI3xZhmA06DDoNJoOOM1uEce+gVIx6HR5FJa/SjkYDBq2W/CoHORU2BqSGE+ZvQIOGAKPusDcjLo/ChpxKAvx0dIg9dnOkQaflpREduOfbrXx+Vbd6+2x5FJW7vtlCUqiJG3snU+lwc9/gNAKN/3/4Wp1uft9byk3zNpEUaiK/yoHTrXhD+fRR6cxal82Xm/IY1fX06+soGodOq0FRVRRFlZufeuRRVN78fS+/7y1h2mWdSQk3AxAR4OftMpRZauXRH3cQ6m/gvkGpbMit5Pzpq7j8n4F24WY9qs3G5opi/s6qoNDioGt8CP1Tw0mLMHvPVx5F5cftBaw6UM7V3RO8XTFsLg+bcivZlFfJjkILDreC7p8vX9BpNUQG+NEpNoiOsUEkh/ofdf/bXB4MWg36o7TwuD0Ke0ut7PrnJqdVhJnWUYFHHERXZHFw33fbSAgx8cmV3QgwHj7eONwerp+9kQqbmyu7Jxz9jReNTqOqx/6qDqvVyujRo5kzZw5ms7kxynXcLA43a/7pB/rHngJcGj3xwSZ6J4fRKzmUjrFBTaqJ8/7vtnJBh5hTaqooVVVZtKuYdTkVDE6LID0+mL8yy/h9TwkRAX50jA0iIcREYZWD7Ao7Of/8lFprJgV3K2qtdbldTuLDAokKNNI2KpCOsUHEB5swGWomAA81GZrkhe+1pXvIrrBxZotwtuZXUWl3ExdsJDHEn6RQE4mh/kQG+BHgp6sVXD2Kytrscr7fWsCW/CpahpvpEhfMgu0FXNQxhmt6Jh1z26qqsru4mvhgU63QCjVdEMZ/sZ751/eq99d8JP+d41I0P8e7D+/9dis39k6u042cOLZNuZU8tGAbl3eO44beycd1s2t1ulm8u5j8KkfNebe4nLTYcHqnhBEbZGRddgUr9pWSWWZDpeb84fQonNMmmv6p4XywMpPcSjtajQaTXkuXuGC6xAfTPjoQf4MOj6riUWp+auattbA1v4p9pVYC/HQMTosg0KhnU16lN1xqAJNBS4XNTeuoABJDTGwtqKLM6kKj0RBhNlBocWLQaUiNMNeM4dDAzsJqNuVV0is5lMs7x2F3K6zMLGNfqZWCKgcWp5vnhrerU79+t0fh3u+2YXd7GN42mnbRgbSOCvBmhDKrk6JqJ60iAprcNeZkz6k3zd3IlEs74d/INeh1zZFNOpRml9uICPA75M1zexS25Fd5m+HzquwE+unpkRhCz6QQko1OUhPjmvRFsMzqZPSna/n5pj5N7kN/JG6PQkZxNUXVTsptLspsLvaVWDH71eyfhbuK6NcinAGpEfyWUczmvEp6p4QxpFUkZVYnW/MtNX17gmqa0BNDTMSHmAj3NxDib6h109Dcw8yiXUVkldvoFBtMiL+evMqaAJ5Vbie73EaJ1UW1083Bo0+jqQmlXRNCuLBDDN0SQthXamVDTgUdY4PoFBdcL+W6+vN1vHxBB+JDDj9bxMnwKDW117FBRu+Fs7nvR3H8+/CbzXkUWhxM6Nui4Qt3CtpRUMW3WwvYnFdJqdVFQoiJJ89pQ2Ko/0mtt7GPxQqbi+X7SrG5PHSKDaLNYWo5dxZaKLI4aB8TRESAH4qiUmJ1EmH2O+x1UVVVVh0o54dtBfgbtPRJDqNVZAAxQcYT6qawLb+KvzLL2F5Qxe7iajz/VIyEmAxEBPixvaCK8WckMaZrvLdWV1VVluwu5rc9JQT66bikU2yjDuw73v1YaXcx/osNxAQZaRsVQKHFyQsj2jdCSWs7JULpbxnFvPdXJp+N7cb+Mhufrc1mTVY5Wo2GTv9Mx3RGUmitC2xzugi+9ftedFoNt/dr0aQ6qLs8CnuKqym1umgbHUiwSc/Hf2fxxfpcuiYEExtkJMRkIMxsoEWYGavLjUeBvi3C6u3uqzntx+ZkwbYCdhRauG9wWr2ud/GuIp5fnEFSqIkKu5v/jeyCn07De39lkuTv4co+bWQ/NlPHeywWWRzcM38rn1/dvRFK17x4FJXvtuZTanUxIDWc1pEB3nN/QZWDZ37dRZXDzfW9kkiPDyasHsc4yDn1+NldHj78O4svN+Yxums8XeKDefm3DDrEBHF55zgq7C5mb8glu9xGSpiZ7AobSaH+jOuZyBlJoSd9XXd5FBRVxaj//+vq0fZjuc3Fx39n0ScljN4pYaiqytjP1nHXgJYEGfUs3VPz7Wu+aDWua45s0n1Kh7SKpMLm4oIZq0kO8+eaHok8MaxNs6lZPJbb+7Xg2YW7uWDGajRA14QQruqeQHsfTKeyOa+SmasPsLOwGqNeS2qEmTB/A//7K5Nym4sLO8bw8829m1QXCHH8zmkbxZTl++o1lL77537WZlUw/7ozCDLp2ZJXyfVzNqDTaLixdxKz/t6Px5DNuDOS622boumKCjRicboprHLU+4C/pkxVVbYXWNhaUMWg1Ihar93icPPh6izmb8nngg7RJIf58/aK/ewutmDQalFUFZNBxwOD0+idEubDVyH+zWTQceuZLbixdzKfr83hu635vHlxR5LD/j9UDW4VicXhprjaSWKIid3F1Xy4Oosnft6Jn06LQachLthEm6gAws1+6P8ZTJsWYaZfy3Dcisrve0qwuxXOaRuFQafF4nAzZfk+lmQUE2LS4/T8f91hiElPSqCG9BQPDrdCmc2FW1HJLKv5+uIbeiUzZfk+eh8oo6TaxZktwunboma2n/pqcWtITbqm9EQ017tBRVFZdaCMz9bmePvjdEsI4bLOJ9804PIorMosIzHUH71Ww/K9pSzaXURBlQOtRoNHVUmLCOD6Xkl0TQipp1d0cprrfmwO7vh6M7ed2aJe+vwt3FnE5+uy+XBM18PWCiiKQl5+AQ/9lsc5baMY1iaKCLMffnrZp83FiRyLm3IrefCHbRh0Woa1ieSOfi1PmcqEg9wehU15lfyxr4y/MssotTppHxNI59hgvt6cR6/kMC7sEMOXm3JZk1XB9b2SuLxL3FEH9zQkOaf6htOtkFtpZ2ehhQq7G5dHwU+vZUteFasOlKHRwOC0SPx0Gn7cXkhEgB9lVhe392vBhR1jDjmvFlXZ+WNnFqUePwKNBsL8DRh0/ww0+yd0qqrK+yszUVWY0DelSbTEnhLN9yfiVDnwLA43a7PL+XxdDiXVTu4fnOa92/m3gioHX27MRVEhIsCAR1Eptbow6rWckRTK6qxyvtyYx4DUcIosTpwehX4twjm7dSRJYSfXR6khnSr7sSlam1XO1BX7mDm6K25FpczmotzmIjHEdMQRq4eTWWrlxrkb+e6GXkfstnFwP4ZFRDJlxX4yiqvJKLYyMj2OW/5zstyWX8U7f+xnRIfoRpuupdrhZkt+ldROHcXJHIseReV/f+5ndVY500elnxItLQVVDl5aksG2gip6JIbQr2U4fVLCas3goaoqP2wr4K/MMs5rF90k5hCVc2rTp6oqFXY3ISb9ET8vzXU/nhLN96ezQKOeQWmRDEqLJLfCzmvL9vDC4gzuHZTKoLQINBoNX23K5b2/MrmjX0tMBi3F1U6Meh2dYk1UOz18vTmPjrFB0uwuaumRFMqgtAgu/vBv/HQawvz9CPHXk1lmo9rppltCCLFBRjKKaybnTwzxp01UAG2iAjHptewurmZHoYXf95YwY1TXOvUjNui0PDCkFVDTKjB1xT4u/3gN9wxIJS7YyGdrc9hWUMXT57bl9WV7KbW6uLpHw0/AvqPQwo/bCyWUNhCdVsPt/VsSszGXCfM2MWN0us/D2fGqdrj5dVcRm3Ir2Vdqpczm4oHBabx+cccjPkej0XBhx1j5amlxXDQaDaH+Bl8Xw6cklDYD8SEmXruoI0UWB68v28tLv2Vg0utoExXAghuPHDgv6iQnRHF41/VK5rpeh/bxVBSV9TkVlFidXNAhhjB/A9kVdnYVWVi2pwSHW6F1VADntYvmsaGtT6gpUqvVcPfAVC7oEMPsDTkUVDk5t20UT53bBo1Gw3tXdGH8F+sJ9TdwQYeGrTHdV2qlZXjTbv05FYxMj2dPiZUpy/dx98BUoKZWaPHuYtZkldM+JojuCSEEm2pmqsirqpnFoV10YKOFWLdHodrpwaOqWBxulu4p4YdtBTjcCue2jWZEhxgSQkzEBdf/zBVCiBoSSpuRqEAjL4xoj8VRM5VQkEl2n6hfWq3mkHn+wsx+dG6ADvJpkQFMGtrmkL/rtBpmju7KqE/WoAFGNGAw3VdqpVdyaIOtX/y/B4ekccOcjXy0OouYID+mrthPj8QQzm0bxeb8Kn7eUeid0zcu2ER2hY3NeVWMaB/N0DZRrNhXyoq9pfRJCePGPsn1MtOHw+3hlx1FzNmQS7ndRZi/Ab225ks4+rcM538juzTYt/wJIQ4lqaYZ+u+k6EKcavz0WmZd3Z0nft7J3I25vH9Fl1rTotSXPSVWxnSVb3VpDBqNhv+N7MK0P/dTaHHw6dhuRATUBL7+qRGHfY7Lo7BgWwGfrsmmT0oYUy/rxK87i7hk5t+YDFoMOi0xgUbS44M5UG5jbXY5Oo0Gs5+O1pEBuDwq5XYXdpeCVgNajQaNBlQVCiw1Az2Hto7k1Ys6SA2oEE2ApBshRJNk9tPz6kUd+WFbAbd9tZnpo+q/P2J2uY3EUAkjjcVPr+Wef5rv68Kg03JJ5zgu6Rzn/duYbgmM6ZaAoqi4FZWcCjsbcytIjw/mmXPbotVqsDjc7CmpxqjTEuJvwKTXolLTPeXgF8dFBhx+gnYhhO9IKBVCNGkXdIhhb0k1D/6wnbHdEii0ODAZtHRLCEGn0VBiddaaN/B4NbeBN6KGVqvBT6uhZYSZlhG193+gUU96fNOY3k4IUXcSSoUQTd5dA1KZtzGXb7bkER1opNrp4d0/a+bhM/vpqHK4efXCDrQ4zKClA2VWPvo7mxCTni7xwXRLCCHU30BBlYOowNNncnchhGjqJJQKIZqFK9LjuSI9/rCP7S6ycP/323C4FTTUDJZqGW4mp8KOXqvhpj7JONwKa7LKeX7Rbh4b1poii5N+LWQqKCGEaCoklAohmr3WUYF8Ob6n93eXR2FviZUwf0Otr3s8t100t53Zgss+WoPF6WbeuJ6HW50QQggfkFAqhDjlGHRa2kYHHvaxAKOe2dd0Z0+JlfgQGeQkhBBNhYRSIcRpJ8zsR0+Zf1IIIZoU+e5JIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvicTJ4vhBBCiCbFZrNRWlqKoii+LkqToqoqJSUlOJ1ONBqNr4vjpdVqCQgIIDg4GK32xOs7JZQKIYQQosnYsWMHM2bMwO12+7ooTY6qqiiKglarbVKh9KC0tDSuvPJKIiIiTuj5EkqFEEII0STYbDZmzJhBq1atOPfcc9HpdL4uUpPjdrvR65tWfFMUhZKSEhYsWMArr7zCs88+i8FgOO71NK1XJYQQQojTVmlpKW63m3PPPZcWLVr4ujhNjqqq3lDa1GpKU1JSCA0NZerUqRQVFREfH3/c65CBTkIIIYRoEg72IZUa0ubpYO2ox+M5oedLKBVCCCGEED4noVQIIYQQQvic9CkVQgghhDiGH3/8kV9++YXMzExGjRrF2LFjfV2kRvX222+zevVq7HY70dHRjBs3jl69etXrNiSUCiGEEEIcQ1hYGGPHjmXZsmW+LopPXHzxxUyYMAGDwcCuXbt4/PHH+eCDDwgODq63bUgoFUIIIYQ4hr59+wKwZs0aH5fEN5KSkrz/12g0uN1uSkpKJJQKIYQQQojGNW3aNBYvXozT6aRnz571Pm2XhFIhhBBCCHFMt912GxMmTGDLli1kZmbW+1ypEkqFEEII0Wy88MILFBYWnvDzo6OjeeSRR+qxRI3n1r++IttSgVZ7YmEwwRzCu2deflJl0Ol0pKen8+233xIfH0/Pnj1Pan3/JqFUCCGEEM1Gcw2U9eHdvpc3mW908ng85Obm1us6ZZ5SIYQQQohj8Hg8OJ1OFEXx/v9Ev7mouamurmbp0qXYbDY8Hg8rVqxg8+bNdOrUqV63IzWlQgghhBDHMGfOHL744gvv73PnzuXuu+9m6NChPixV49BoNPz666/873//Q1VV4uLiuP/++0lNTa3X7UgoFUIIIYQ4hrFjx552E+YfZDabef755xt8O9J8L4QQQgghfE5CqRBCCCGE8DkJpUIIIYQQwucklAohhBBCCJ+TUCqEEEIIIXxOQqkQQgghhPA5CaVCCCGEEMLnJJQKIYQQQgifk1AqhBBCCCF8TkKpEEIIIcQxuFwu3nrrLa677jpGjRrF/fffz44dO3xdrEa1Y8cOLrroIubMmdMg65dQKoQQQghxDB6Ph+joaF566SVmz57NRRddxDPPPIPNZvN10RqFoihMnz6d1q1bN9g2JJQKIYQQQhyDyWTiyiuvJDo6Gq1Wy8CBAzEYDOTk5Pi6aI3il19+oU2bNiQmJjbYNiSUCiGEEEIcp9zcXKqqqoiLi/N1URpcZWUl3377LWPHjm3Q7egbdO1CCCGEEPXo1i83kVNhP+HnJ4SYeHdkl5Mqg8Ph4LXXXuOKK64gICDgpNZ1PPI+vhVXaTYazYnVKRrCEoi79t3jft6nn37KxRdfTGBg4Altt64klAohhBCi2TjZQHmy3G43L774InFxcYwZM6ZRtx03/l3cbjd6vR6NRtMo29yzZw+7d+/mlltuafBtSSgVQgghhKgDRVF4/fXX0Wg0TJw4sdGCoS9t2bKFnJwcrr32WgCsVis6nY68vDzuueeeet2WhFIhhBBCiDp45513KC0t5ZlnnkGn0/m6OI3i3HPPZeDAgd7f33//fWJiYhg5cmS9b0tCqRBCCCHEMRQWFvLrr7/i5+fHVVdd5f37U089RceOHX1YsoZlMpkwmUze3/38/PD392+Q/qUSSoUQQgghjiE6Oprvv//e18XwuYkTJzbYumVKKCGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCNAlabU0s8Xg8Pi6JOBEulwvghKfLktH3QgghhGgSwsPD0ev1/PLLL5x77rmnzVygx+PgNzo1JYqiUFJSwg8//IDRaCQqKuqE1tO0XpUQQpwitpcX0C4k+rT4xhch6ou/vz833HADM2bMYMeOHb4uTpOjqiqKoqDVapvkuSUtLY3bb78dg8FwQs+XUCqEEA3g68zNrCjYx8SOAxkW36ZJXkCEaIratWvHM888Q2lpKYqi+Lo4TYqqqpSUlBAREdGkzilarZbAwECCgoK8XTBOhIRSIcRpq8plZ0NJLgNiU+t93ZPSh1LmsPLWtuW8uXU593QcIOFUiDry9/cnISHB18VochRFwc/Pj+jo6JMKf03VqfeKhBCijvZUlvBr7q4GW3+Y0cxT3c7l80Fj+bNwPyMWzmBRA25PCCGaMwmlQojTVq61kgRzcINv59/h9Le8PVy86EPWFmc3+HaFEKI5kVAqhDht5VgriDeHNNr2woxmJvc4j3fPvIwPdq3kmt9nsaeyuNG2L4QQTZn0KW0gdreLccu/IMIYwLtnXu7r4gghDiPHWkG3iMbvtxZvDuF/Z45kR3khj637mUhTAI+nDyXaP6jRyyKEEE2FhNJ64FI8ZFQWs628gIzKYuweNyuLMrmn40A+zlhDga2KGLnYCNHk5ForiW+E5vsjaRcazReDr2ZlYSbXrZjDmdEtuLfjIPz1JzadihBCNGcSSo+D3e1iZ2UR28sL2FZewPaKQmxuF3qtltbBkXQIjeHM6BaY9X7c2KY3CQEh5FRXsCx/D6NadvV18YUQ/5FvqyK2Cdww9olO4YehN/B15mZGLJzOLe36ckWLdBmp30B2VxShAm1C/n+C71KHlV0VRXQNj8fUDG4KVFWVz4c45UgoPQyr28m2g8GzvJCdFYU4FDcmnYG2IVF0CI3hkuROPBQSRYDBeNR19Y1OYcau1RJKhWiCFFVFq2kaXes1Gg2Xt+jCiMT2vLltORcsmsEz3YbTIzLR10VrlhweNzN2reLX3F14VAUNGsKNZnKtlbQOjgQgo7KYATGpbCnPx6V4SA+P45kNC0kJDOPi5I4oqso+Syl51kqK7NUUO6pxKm6CDCbC/PzJt1Xh8LhRVBWdVkOHkBhaBIVjdTsptFkotNf8aNDQNzqFcxPakmkpY0dFIcEGE6lB4bQOjqJlUDgGrQ6b20WR3QKA3ePmp+wdrCjch93jQoOG3lHJtA6OZN7+TVhcDgDOSWiLQaMl11ZJgjmEQbGpdAiNaTKfayGOx2kdShVVIdNSxqbSPDaV5bG5LI9qtxN/nYEOoTF0DI1lbGo3WgdHnvCdc6vgSDKqTq+BDPnWSnZVFtEiMJykgFDv3bxHUSh1Wgn3M6Orw/xqVrcTl+LB+K9lVVWlzGkj11pBudNOy8Bw4s3BUmMgjpuiKjTFj41Jb+DhLmcx3tqTJ9b/grpD5dnuw4nzYTeDpsLudvHOjj9YUbCfdiHRXNf6DFoHR/J3cRZf7F2Pn1ZHl/B4qt1OZu9dz3WtezF38DX46fQoqkKJw0qE0ewNbA6PmzXFWdzUtnetLlbbywtYmLsLvUZLalA4/WNaEmUKINIYgEGrw+J2UOawEW0K9F4b3IqH7eWFZFWX4683EG0KJNo/kEhjAB5VYUXBPr7J3ELLoHDOimtFpcvBnspiFuXuZr+lDLfqwajVE+MfhAbQabQMjW/NDW16EWgw4vynW1hGZQlTe19CrDkYq9vJioJ9qMCg2DQyq8v4OGMN28sL0WhAg4YoUyAJ5mAijGaslmpCyoJxKQp5tpquK70ik+kWkYBRd1rHAdFEaFRVVY+1kNVqZfTo0cyZMwez2dwY5TphiqJQWFh4yMSylU47m8tqwuem0jwOVJehQUOLwDC6hMfTJSyOjmExBBlM9V6mEQuns2DYjfW+3sZS6rDyd3EWi3N3c6C6nDj/ILqGx9MqOJJiezUrCvexo6IQDRoUVSXM6E+n0Fgyq8s4YClH+8+VX6vREOrnT561kqSAUJIDQylz2LC4nQTq/QjxM1HisJJVXY5L8aBzK2DQ41YUwoz+lDtt3tqOeHMwoX7+7K0qId9WBYCKSowpiJZB4YQYTESaAmgXEk3bkGhvHz2PoqCiotf+//cpn0gzmEdRsLgdBBmMUiNxFEc6HpuC7Opynt2wiPf6jfR1UY5qXXE2j6//mX7RLbm348BGb1puCvtQVVXmH9jC1G0ruLXdmVyU3JENpTl8uPtvMi1ldA6LY1yrHqgqbCrLA+CKFl3wO82DlqIqFNmryamuoMhuoaSsjOCQYEx6A7H+QeRYK1lVlMm6khxciof2ITEkB4byR8F+qt1OTDo9HUJjSAwIRUNNbf7BM6VHVfGoCoF6I72ikmgdHCnnwkbQFI7HE1HXHNmkj1hFrfl6seP5oHsUhQxLKcus+Wwpz2drWT4OxU2wwUTnsDi6hMcxIrF9rRq8hhZkMFHptBPsV/+B92QoqsLuymJyqisodlRTZK+mzGEj3OiP1e1idfEBLC4n4UYzPSITuTqtO2lBEeRYK1lfksPvBXsJ8/NnfKuex91cdMBSRoGtijCjmUC9Hxa3kwqnnXCjP0kBoWjReA88NFDpdBBq9D/qOlVVpcBWxX5LGZUuO4V2C1/u38TOyiLsHhdQU/tQ6bRzcXInOoXF8sGuVZQ7bQD0i24BwMqiTA6eegP0fnQMiyHWPxi7x0Wl08HakmycipswPzOVLjuKqmLU6WkVFMHguDTOimsttQ4nwO528VPODloHR9IpLK7e1lvptJNnq8Sk03PAUk6508aFyR3ZXJZH57DYettOQ+kemcgPQ2/gy/2bOH/hdO7rNIgRSR18XawGt6JgH/P2bWRPVQkqKn2jWrBg2I3eG8zeUSn0jko55Hmdw+vvs9PcaTVaYvyDiPEPqgkz+tphplNYHOcmtAVqrgc7KgrZV1XKuLSehBr9sbldbCsvINdagUrNOfZgLZZeq0Wn0VDutDNt+59kVBWjQUNSQCitgyPRa7VUOO01Py47To+b7hGJnB3fio6hsT5p3cq1VrCuJIdOobG0CApv9O03hDKHFaNOj1nv5/2boipUOh0EGPxQVJUSRzWx/kGoKkzZvoJb2vZtsoMpm3RN6Ue7/8ZPq2NsWvfDPm5xOdhYmsuG0lzWl+SQa61Eq9GQoDfTOyGNruHxtA+N8fmb//T6XzkvsR29opIbfdtuxUO5046iKrgUhc1lefxVlMmGklzcqoc2wVEkBYQSZQokyhRAmNGfEnvNh/yMyCSfBemGvBv0KApf7t/EfkspY1O7kxQYilvxsKroAAC9opIx/FOTanE52FZeQIGtCn+9gQC9H+nh8bVOAFATqHZXFrMwdxe/5WfgrzPQMjCcxH9O0B1DY4nxD6TIXo1b9dQE76OE+FKHlVxrRa1wpqoqTsXjDbyZllK+ztxCkMFI57BYOobGEniEPs5Wt5Nyp61R5uQ8eEpZX5LD29v/YE9ZIWfGp3Jd616Y9Qb+Ls5idVEWudYK2oVGMyAmlb1VJXy0+29GtujC+tJcAK5s2ZVqt5M8WxV9opLpEh4P1Lw3n2SsYWNpLr2ikhmR2J7kwDBUVSXTUsbmsjx2VhSxz1LKvqpSAg1+JAeEYfe4iPUPpshuoW1IFMWOaobGteHMmBYN/p7Ul2qXgxc2LWFHRSEv9RxB2j99IxtSY9fMlDmsPPD3D5j1ftze/kypgasnjbEfFVUhp7qm+5aKSrDBRIifiWCDCb1Wy9ribBbl7mZ7RQG6f/ZpuNFMjCmQGP8gov0DCTaY2F1ZxOayfG//2tEtu3Jlajc8isL60hyiTIGkBR3+u98PWMpYUbCPPwr3s7eqBL1W623Fi/EPolt4PKuLs3CrChPa9qF/dEs0mppBj0X2atKCIo54Hm0sVS47eypL8KgqiqqwvGAfi3J3o9GAqkK3gCiCAgNZVrCXMKMZm9uFS/EQbw4m11oJQNg/lUsAYX7+5NkqcSsKY1O7Mb71GY3+muqaI5t0KLW7XZy/cDo/nXMTBq2W5QX7+PbAVvZWleBSPN6A0DU8nm4RCcT5B6OqapOr2p63byPVbifXNuAHochu4X87/iLXWsmIpPakh8XzyZ41/Ja3h+TAUHQaLVo0dAyLoU9UCunh8U26Nq+5NlEcVOWyk11dwYHqcnZXFLGlPJ9Cm4Vo/0B0Gi2ZljIMWh19o1NQVZWVRQdweNxoNTUnzyCDkXCjGY+qMKX3JeypKuHuVd9i1hswaHWEG81UOu1c36YX1S4nW8rz2FKWj9Xt8p64NBpoGRhBVnU5KioaNAyIacm9nQbV62v1KAq7Kov4cv8mlhfsw6jToao1I5tvbt2bELvCLqx8mbkJh8dNz8gkekUmkRQQyuayPFYWHSDSFMDVqd29TdMbS3NZmLOLED8TUaZAFuftZk9lCWa9AUVVGdeqJ2dGp7CiYD8/5ewg31aJqkJKYBidw2JpGxJNalD4YcO/qqqMXfY5W8ryWX/xxFpdOZqLjMpiHlqzgA6hMTzS5axDbpLqU0Mfiwdr6LaU5fN3cRYbS/N4ptu59Ik+tBZUnLimeE5VVIUyh40CexUFNgsFtioqnHZaB0fSOTyOKFMgDo+bD3ev5rsD29BrtZwRmUSBrYrdlcXEm4O9AbLYXk2F005SQCj9Y1rSL6YFLQPDj1gju7+qlM/2rOPv4iwUFOL8g4k0BZBRWUy12+ntduZRVLqE19QoF9urvWEXarqMHUxQBq2OpIBQ2odG0zsymQ6hMawpyebrzM3k26poGxxFm5AoiuwWcqorqHI7cHg8KKpCYkAo6eFxhPuZmbNvI/m2KrqEx6H7pwy9IpMZGt8aP50em8vJT7s2EBwaysDYVG8XFZfiodBmIdY/qE5jNhrbKRFKAX7K3s4Lm5Zg1OnpH92SUS3T/2kaOPyFpCkeeFvK8vg0Yy0vnXHBSa/L4XGzt6qE3ZXF7K4sIqu6guzqchRV5bb2Z5IWFMH8A1vZWpbPZSmdOS+xXbMcBNQU92N9s7ld/3QVqKmdPRgsFFXxBqmfs3cwedNiIoxm3jtzJDH+QZTYqylz2mh1jFoyp8dNVnU5ceZgzHo/VFXlyfW/oNNoeaLrMDQaDRVOG5VOB1GmgDr3VcyoLObbA1v4PX8fCgpatLQKjmB4QjuGxreu9Xmrz/1od7twq0q91GIcHB3d3JvwfsjaxutblnF7+35cltK5QY71hjwW91eVcuMf80gPj6NLWBzp4fGkh8c3y3NWU3cqnlPzrZU4FDcAIQb/Y3bxOhGqqrKuJIeFubuI8Q/kzOgWtAmOOuQz6lI8ZFrK2FZewKqiA2wrL6BbRDyXJHciJTCMXRVF7KwsIsYUREJAMMEGk7di6IClnE1luRTYLFya0om2IdFHLE9z3Y+nTCiFmp1tqGNtRlPcYQ6PmzFLP+Obs6894XUcsJTx0JoF2DwuWgdH0jo4itbBkSQHhBLrH3TMqamam6a4H08VL21awuK8DPQaLcF+NVPbFDuqcXhqTu4He42pas3gtDbBUdg8Lg5Ul6GoKskBNdPlDIlrdczadtmPDc/udvHqlmX8XZzFiz3Pp31oTL2uv6H2YbnDxmVLPuKjAWNIDgyrt/WKw5Nj8dTQXPfjKTHQ6aC6BtKmyqjT4/znbu5E7K8q5foVc3jvzJG0/tdkz0KciIe6nMWDnYfUqTbKpXjYU1mCv15Pojm0STYLne5MegOPdR1KpqWUh9f8SGJACE90HdYgM4nUF5vbxfUr5vBizxESSIUQXnKFaSQmnQHbP52Oj0eR3cL1K+bw4YDREkhFvalr86hBq6NdaDQpgeESSJu4lMBwvhh8NcPi23Dxog/5fM866tAQ1uj2VBZz0eKZTGjX1yeDP4UQTZdcZRpJ25AodlUW1Xl5q9vJh7tXM+q3T5nS5xJSApt33zchROM4J6EtP59zE7nWCkYsnMHGf2YzaAoyKouZ8OdXzOg3yjsVkRBCHNQsmu9PBV3C4lhXkk36P9PaHEmhrYoXN//G9vICRrfsyoJhNzToyFohxKnHT6fngc5DGJvanUfW/kionz9PdzuHMKPvvvxEVVUmrv6W6f2ukCZ7IcRhSU1pIxkW34ZfcnYddZkFWdu4+vcvGJnShZ/OuYlrW58hgVQIccISAkL4ZOCVjGzRmVFLP+XdHX/iURSflOXD3X8zJLZVs5/xQAjRcCSUNpIIUwD+Oj33rv6On7K3H/L4ysJMZu7+m++HXt+sJvMWQjR9A2PT+HnYTeg0WoYv/ICleRmNuv0DljLm7NvAXR36N+p2hRDNi4TSRjS93yhubNOLWXvXsyj3/2tNSx1WHl37IzP6jWrSE9oLIZovnVbLzW37MHfwNXx7YCtjl31OpqW0QbdZ5rDy4qYlXLdiDu+eeXmz/KICIUTjkVDaiHRaLR1CY/lf38t5YdMSyhxWVFXljpXf8FLPEQ0y8a8QQvxbmNHMG70v5vH0odyz6jueXP8LVrez3rczd98Gxiz9jC5hcSw892ZSgyLqfRtCiFOLVMv5QIDByOTu53HdijnE+gfRNyqFM2RqFCFEI2ofGsPXZ43nh6xtXLBwBje37cPoll3r5duUNpfmMW//JhYMu0FqR4UQdSY1pT7SJzqFJ9KHcU1aD+6UflZCCB/QaDRcmNyRn865iQPV5VywaAZ/Fx046fXO2L2aJ7sOk0AqhDguEkp9qHtkIv1iWvq6GEKI05xRp+fBzkOY0W8UH2b8zdXLZrGnsviE1qWqKlvL8+kYGlvPpRRCnOqk+V4IIQQAseZgpvW9nJ0VhUxa9zMxpkAe7XLWca1jfUkO3cIT6qUbgBDi9CI1pUIIIWppGxLN7MFXM6plOuNXzOGtPSvrPBhq/oEtXJrSqYFLKIQ4FUkoFUIIcVj9Ylqy4OzraR0QwcWLP+LdHX/i8LiP+pyVRQfoLQM3hRAnQEKpEEKII9JoNJwf25qfht1AgN6P8xdO5+1tK7C7XYcsm1FZTGpQOFqNXFqEEMdPzhxCCCGOSa/VMa5VT3455ybCjGZGLJrB61uWUe6weZf5OnMzl6Z09mEphRDNmYRSIYQQdabX6rgqrTu/nnMzLQLDufy3j1mYswtVVVmcu5uz4lr5uohCiGZKRt8LIYQ4bjqtlstadOa8xHbc8MdcZu1dR9/oFAwyN6kQ4gRJKBVCCHHC/PUGPh1wJcvy9zAoNs3XxRFCNGMSSoUQQpwUnVbLWfGtfV0MIUQzJ31KhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSSoUQQgghhM9JKBVCCCGEED4noVQIIYQQQvichFIhhBBCCOFzEkqFEEIIIYTPSSgVQgghhBA+J6FUCCGEEEL4nIRSIYQQQgjhcxJKhRBCCCGEz0koFUIIIYQQPiehVAghhBBC+JyEUiGEEEII4XMSShuQI2cb9qzNvi6GEEIIIUSTJ6G0AeV/cR/5n93p62IIIYQQQjR5EkobiCNvJ4aweLTGAF8XRQghhBCiyZNQWo9UxYOnugyA0l/eIHzY3WhNQSiOah+XTAghhBCiadP7ugDNneKwUr11IVXrv8NVloMuIByPpQRDeCKm5C4YIlvgLNqPKbGjr4sqhPiHYrfgKs3CGN/e10URQgjxDwmlJ8BdWYRlww9UbfoR1e0ksOMwoi55EkNEMgCqqnqXNYQn4S7NAgmlQjQJquIh+39X4akqIvH2uRjCE31dJCGEEEgorRNVVXHkbMWy/nusu5ajNYcS1PUC4q/7AF1A6CHLazQa7/91QVG4q4oasbSnPlXxgKqi0R3546sqCrY9K1E9LnT+IaDTozUFYQhPQqNt2F4rqqJQumgqqstOQKdz8FQUUL1tMY6crehD44kY8RDG2NYNWoZjqVr3HegMBHYZXuvz2pypbheushw8lmI8VcUojmrMbQeiD46qtVzx988T3ONSdEGRVK75iohz7vZRiU8dqqqi2KvQ+Qf7uihCiGZMQukRqG4n1TuWYdnwPY68nRgTOhLU9QIizrsfjd5Q5/XogyKxH9jYgCU9fXhslZQtnoZl60I0Wj2m5HRCB1xP2bLp2DPXYQhPImbMK2iNAeTOvAlDRAq6wAgUWwWqx41iq8BVmvX/K1Q8BPcaRXDvMWj9TLjKcrFnrscQ2eKkuluULZ6G6nJgatGDqr+/QhcUSciZ12BM6owzfxfF3z6L6rKDTo9iq0QfGkfkhZPwi2px8m9SHTgL9lC+/EOMSZ0pX/YB0aNe8nlIPhmK007xD89j27sav+hW6IOi0AVHoTGYyHnvakL7jSPkzKtQVZXSX97AXVFA5EWPodirqPjjU5BQekSq24WrNBt39k6qizQo1nI81WUo1WWoqoIhNB60Oir+/AxPVRFJE3/AEJ7g62ILIZopCaX/4q4swrLpRyybfkKxWzC3G0TYWbfhF9f2hGuTpKb0UB5LKZbNv+CxlqE1BaPRaKjesRR3WQ5odQBo9H7Yqy04/Aze916j9yO03zgizn8QjVZL9Y7fKf31LYJ7jyL2ylex7V9HzrtjURU3UZc8RUD7wUcth+K0UbHyC7KmXIpGo0EfEoupZU8sm3/GVbgHc/sh+EWnEdBuMLrA8Dq9NlXxULVuPskP/IpGqyWw49Bajxvj2pIw4ZN/BsRp0AWE4sjdTsHs+zFEJBN+zt34RaYc93t6PCpWzyF82J0EdDgLZ+FeCr96DH1QJFGXPo27qoiirx7HU12GMb49Eec/gD4sAVfhXlTFjTGubYOW7XhZM1ZS+OWjhA+7i+jLnjnk8dB+4yj69lmy3xmNx1ZJYJfhxF4zFY1Gg84/GMVe5YNSNy0eawWOvB04c7fjyN2Os2hvzU0ToNEZ0IXE4dKaccQkoQ8IxxCehC45HdDgLstBcVpJvPNL7PvXUr58JlEXP+7bFySEaLZO+1DqKsmicu3XVG9ZiNYcQlCX84kd9w76wIh6Wb8uKBJPVXG9rKs5c5XlUrXmKyxbFqI1mAjsOgJDZEsUeyWq20XkRY/VCmNuh5Wi0gpiYmLQHqG5PaDdQALaDfT+7t+iO8kP/AKKUqfabK2fP2EDryds4PWHPKY4rNj2rsJZsIec98fhF9eOiHMnHrMWyLLpJwK7nHfMLgK6gDDv/43x7Um680usGX9R9NXjKC4bYYNuJKDjsAbpamDf+zcRw+8DwC86lcRbZ1G9czk5749HFxhO9Mjn8YtJw7ZvDQVf3Idir8IQlYriqMYvOo2oix+r9zIdi6qqFM57FI+lmMgRD6EPT6L422dwlWaTdOdXtd7Pf9Po9ERf9jTOov3ozKGHdLfR+ofgsVWeFs3OqqriKt6Pfd8abPvX4sjZBqqC1j8YY3x7jHHtCOl3DX7RaWj9/L3PUxSFwsJCwqOjj3gsApjbDqLkx1ca46UIIY7BWbgX687fCR1wra+LclxOy1DqLNpH5d9fYd2+BH1ILEE9LyPszlvQGoz1vi19YCQey6kXSlVVxbZnJfYDGzHGtsEvpjWqx4Vir8JVlot1x28483f/0//Tgy4gnKCel5N4+5xaF7wj0RpMaDSVx10ujVbnrW09GVqjmYD2QwhoP4SwwTdi3bWCgrkPojqqCT7jCoJ6XHLYOWjLl39E3LX/O6Ftmlv1xdyqL+6KAsp/n0nxjy8T1OV8wofdiUbvd7IvyUt1Ow75rAe0HUBA2wG1/ubfsieJt8+p9beC2Q9QtnQ6YYNvrLfy1EXlytlo/fwJHX4fhd88iaeigLCzbyN65OQ6Pf9IXSP84trizN+Ff8ue9VjapsFjKcW6+w+su//EkbMVVA+GyBb4t+hJSJ+xGBM6HLVf9vHSaLVo/MwodgtaU2C9rRdA9bhBqztii5WzcC9Va79BFxSFMbETqtuJq3APjvxdaLQ6/GJa1wxE1enR6o0YEzuh0elxFuyhetsizO3P8nZhcZVmo7qd+EWn1utrEPXHXVmIzhxar+fFY3GVHECxWzAmdDip9XhslSiOavQhsfXWn99TXUb5io8J7nkZhohkVLeLvI9vRecfgiGyxTFbDZuS0yaUOvJ3UfX3V1TvXIYhPJngMy4n4py7GvxDrdEbUD2uBt1GY3Dk7aTomydRFQ/6wAicBRmYUnthTuuDbf86KlfPRWMwoTUGoguOJqTv1Rjj2x9X/9umzNymP+Y2/fFYK6j8+0uyp43BLzqN6CteROtnAsBVmoPWYEIfFHlS29KHxBB54SNEjHiQylVzyJp6OYm3z/Nu52S4q4rRBZ54+aJHvUTehzeh9fMn5MyrTro8deEqzaF8xUckT/wBjd5A4i2f19u6jXHtcOTuOCVCqavkANady7HuWoGzeB+6wAjMrfsR0ns0xoSOjXIsmlufiXX3nwR2PqfW390VBVSu+QpTix6Y03p7/664HNj3/Y3iqMaU3BV9SAyq24ll0094bJUYwpOo3rYYW8ZfoNEQPfplTAmdKP/jE6rWfo3WGFhzk2UOJaTPlXiqy6ha8zUagwm/mFYEdbsIFA/Ogt1Ydy1HVTwodgv2eY+AqmCISCaw87mULHgRj7UCUNH6mdHoDChOK+FD78SUnI7GGHDcAcJjraBs6fsY4zsQmH5+sxlQ6LaU1GopdJXl4irNwr9Fj2PexKiKgquo5rN3sFXCXVFA5eq52A9sxNSiByF9x6Iz17RQVPz5Oaakzpjb9K9T2RSXg6Kvn8BZmIGnuoyQvlcROuDamsqIBmTbv47CuQ+hD0tEFxhB9MjJKLYKFIcVnTkErX8wqsuOxs+MRqOpqVhY8RHVWxcBGsztBxPc/RIqVs3Gvm8N+tB43BV5oNER0G4QoQNvQGsOoXrbYipXfoEhKpWwIROOeC1RnVZKfn4d3A5UtwNbxkpCB15PzgfXEn/DDMp/n0nogOsI6HA2eR/eRED7wageN1lvXYJfTCtixr7RZD+PTTqU2g9sRB8cgz409oSe78jZRuXfX2Ld/Qd+Ma0I7nk5Eec/UK+1A6eCmtkFtmEIS/D2caxa/z1+0WkEdhmOdedySn55g/gbZqILDMdTVYQ+LNHbtBzk4/I3Jp05hLBBNxA26AaqNv5I9tuXkzDhc3QBoZQufofQQfVXg6jR6mpO4IER5H96O/E3zDjpddoPbMSUnH4SZdISd9375H04AcVlJ2zQDSddpqNRHFbyPryJ2KvebJBQ5RfXlqo1X9f7ehuaqig4crdj3bUc264VuC0lGMKTMLcdQMSIhzBEpvjkomNuO5Cqdd/WCqVVG3+i9JfXCTv7dir/+pzSX98i5MyrsWz6CWfBbsxtBqA1BVHx1yzclYUABKWfjy44GtveVQR0OJuoy57BU1VMwRf34a4qJLjnSJLumo/qcQIc85vz/NN6HfXx4DNG4q4sQqPVefuPO4v2U77iI0oXTUVxWP9/YdWD1hxK8BlXENjlvMO2sNmzt5D/2V1EnDsR+4ENlC2ZRsT5D9XqbnRwVhdn3k5MLbrjF9Xy6G9uPXMWZ1I4+wE81nIiRjxEQPshFMy+H2f+TgyRLYm96k0qV8+l4q/PMbXoQeGXjxLc/RLChkzwVuaoioLHUoy7PA/L5p+p3rIQv/j2KNVleKzlqIoHXUAYwb1GEXXGSKw7l5Pz3tX/hH4bof3HU7HyC8p/n4kppTu2/WvB48LcbhCBnc9FH5aIqyQTe+YG7PvXYs34k4hz7iZm9Euobhdly6aT9eZFaHQGoq94EWN8u+N+H9zl+aiK54hds5zFmRTMeYCk2+ehCwyncu18st8eiT4oCo0xAMVagWKvRKM31vRR12jRmoII7T+eiPPuBzRYNv9M2dL3CUwfQdSlT3uPTdXjxrJxAbnTr0Vx2fFP60PUJU/hyN1O7vTraj4TGi3Owj3og2OIvPhxVBWqPr6WmOF3YwiLByDqkqfQaLWYW/cj54Px+LfoQUifMQDowxKxZ26geucygnpcitYU9E+APnaLpS9o1H9PqnkEVquV0aNHM2fOHMxmc2OUC4CKlbNR3Q5C+4+v0/KqqmLLXE/+sk/RFmzBlNCRoJ6XY259ZoPfSR3NgTcuJHni9z7b/pF4rBVYd/9J2eJ3MESm4LGU4rGW4xedRlCPS3AWZFC9+RcMUS2JHjm5Ts3u9eVgP7boY/Rj8zXbvrUUzn2I0IHXU7n2G5LumNcg2yn8+omaPn99x57Uekp+fh1jUjqBHc8+qfWoikLO/8YSecGjmJK7HHG5k9mPquIhe9oYws+6lYAOZ51UeY9YPpcDd2kWfjGtGmT99UVVFBxZm6jeuQzbrj9QHBb84ttjbjMAc+t+h0x7VZ+OZx8erI1JvvcHACxbF1O+bDrxN37orel35O7AsvlnAjoOxZTYqcHK3ZBqav/mYdn0I7rgaIxx7dCHJaA1BeHM3Y51z0rir/8AfXC0d/mSX97EnrkOfXA0WmMgzsI9GBM6YEpKx5rxF+7yHNDo0Gi1GCKSibxwEvqQmJMua/WO33FkbcQvoRPlVgem0l3Y96xEddqIHv0ShrAECr64H3vWRkL7X0vYkJupWDmbssXT8G/Vh+grXkCj1aEqHir+/JzyPz7GGNceV0kmAPrgaPQhcfi37kdg53PqdK31VJfX6t/tyNmGs2hfTYuFTo91+29YNv2Mu7IAQ0QypuRumFK61dRaH6ZSyV1RQPa7V5J4y6w6V2Kpqkrx989jP7ABrcEf1eMk/Jy70YfEYtu7GtvuP3AW7QfFTey4aY02O8q/2bO3oNEb8YtOxZG7g+LvnkNxO9Cd/SBx7XvV6ZzqLNxL7ozr0ZqCSLr72wafEvFI6pojm3QodZVkUfTdc8Rf994Rl1FVFfu+NVSu+Qr7/rUYE7vgTjuLuB7notM3jRrRphZK7ZkbKPzyUbT+wfi3PIOQftd4T55NRXMJpVAzxVLVxgWEDri2wQbMqB43B964kLhx75xUX7ec98YRM/b1k+5iAOAqzyNv5k0kTfz+iLVyJ7Mfi759DkN4YrPrqF8fVI8b+4ENWHcsw5rxJ6rLjjGxM+a2gzC36deoA7OOdx8eePNiEm+bg0arJfO180i+53u0xsa7bjQ2V3kersK9uMpyUOxV6EPjCOw8/IgXf4+lFMVhQR+WcNgApyqemnP03IeIGPEQgZ2GnXDZyld8QvXWRQT3uRL7gY1UlRYQ2eVsAtr0O2RwoOKw1mk/qYqCq3g/hvDERu3TeSyOnG3kf3EvSXfNr1NXp+IfXkT1uL0DN53FmZQv/QCPtQz/Fj3xb9Mfv+hUn1ZoHc6JnFMtWxbWdI9pwJvXY6lrjmwaqe0IDBH/fBvSfxycGL1yzVc4sjZhatGT4DOuIPqKF2pG6RYW+uxu4LA0WlRFaRJlqlzzDeXLZxI/4dN6m2HgdOcXk0bEOXc16DY0Oj1x498l7+NbSL7n+xNuyvZYiuslkAIYQuPwT+tD1bpvCe5xSb2s8yBPdTm2jD+JvHdBva63qVJcDuyZ67DuXI4t4y9UjwtTcteaaemGTKj3gUMNyT+1F7a9q7EfWE9o/2tP6UAKNceBITSuzsvrAsOPOsWcRqvDv2UPku75joIv7sOy+ZeaGSeOI1B4LKWU/T4TZ/4u4m/6EI1WR0Dnc/EUFhJ4hDBT1/2k0Wqb5CAwY0IHIs65h5x3ryR23Dvepu3DqVzzDa7iTGLHT/P+zS8yheiRzzVGURvdydzYNLYmHUrhnylbLKXoAsOp+OsLKv76HFAxpXQntN94jEmda9XS1KHit9FpTYEoDkuj1m7YD2zEkbcDV3Em9sx13r4uxoROJN35VZO6wxV14xfVgtAB11H4zZPEXPH8UZdVHNVYNv2MI38n4UPvqJmT02lHYzj5wVL/FjH8Xg68eSFB6efX62eqbMm7hJ19e5PtjH+yXCVZWDP+wrZnJc78nWh0Bkwp3TG36U/40DuadZALPmMk+Z/fA6pC8n0/+bo4zZbWaCbu2nep3raEvI8moAuMJGzIBEwtuh/xuHCVZJH36R1o/cwEdhlOxHn3nbLH0OEEdhmOLiSGvI8m4N/yDCLOux+t0YzqduKuLAKNBuuOZVSunkvi7XNOq/emuWjyoTSox6WU//kZHksJir2KpHu+bXLV6cei8w9BsVY0Sig9OBWE1mjG1PIM/NN6EzZkAjpzSINvWzS8kN6jsWz8EVvmevxTugE1zX3O/F24q4qxZ67Hun0JAAGdh2OMa0/+J7eTMOFTHDlbMCZ2rtfyaI1mQvuNp3TJu8f1dZ2q4iHv49tqWkJ0BkyJnQjuPQZTUmdUt4vqbYuJGPFQvZbVF1RVxVVyAMeBDTWDNbI3ozgsGMKT8G/Vl7DBN+EX27ZJtKLUF2NcW8KH3YkpsdMp9bp8JaDDWQR0OAtH7nbKl39IwZwHiBnzKv4tutdazlm4l9yZNxJ//fQmWZPZWPxTupF0z3dUrfuWrLcuRqP3Q6Pz8/bPNUSnkXDbbKmYaaKafCgN7nEpuTNvxNxmQKPPjVhftOZQPNZyDBFJDbodVVXJ++R2AjoNI6T36AbdlvCd6FEvkjfzRpIm/oDqdpL7/jj0oXHoQ+MxxrcnbNDntUYjV29bjHXPKuz71+Lfoke9lyfkzKs58OpwQs68+ohdQhzZW8j/+BbMbfoTdekz5H9+N+Y2/Qntdw2q24U9axPF3z2HqeUZaPQGgnuPbvaBJnfmzbjLc/8ZqNGVwC7nEXH+A826FrSugtLP93URTjnG+PbEjH4Zj6WUnA+uJajrBTjyd+HI2ohGbwSNloSbPm7w60xzoNFoCO5xSb13KxINr8mHUo3eQMLNH/u6GCdFZw5BsVU0+HaKf3gBY1xbCaSnOENoHOZ2Qyj95U2sO5cROngCQennHXH5qIseI/fDm9BodcTf/Gm9l0ej1RJ9xfPkf3IHCbfOOmyTWNH8p0i8fR62vavJfHEwYYMnENrvmprn6w34t+xBwm2zKfvtPVyFe4ge9WK9l7OxxV33njQPinqnCwwn6c6vaiZL73UF5itf9XWRhKg3TT6Ungq0/iH/TMzcMFTFQ/H3z6PYKokc/XKDbUc0HRHn3EXJT68ReeEk/FPPOOqy+tBYAjoOQ6kuQ+ffMLPK+rfsiSmlG0XfPIXWFEj15l8IP3ciAV3Ow52zCX1IHIbwBAzhlxLc89LDrkOj0RB+1i0NUj5fkEAqGopGb2i2LYdCHI2E0kagM4fWa01p+R+fUbHio5pBKxoNqsdF8BlXEHrRY3IhPE1o9H5EXvhInZdv6NkBACLOf4DK1XPRGgMIH3YX2VMvQx+RjH3pNBKveb3Bty+EEKJ5k1DaCLT+IbiK99XLuirXfI1t9x8kP/Brs+9zJ04tGo2mVteR+Bs/JOe9a9Al9cIvsoXvCiaEEKJZkFDaCHTmEBz10HyvuByULnqb5Pt+kkAqmjx9cDRJ9/1EYWGhr4sihBCiGZBk0wh05lA89dB8X7VuPsFnXHHY71oWQgghhGjOJJQ2Au0/85SerIq/Zp30d58LIYQQQjRFEkobgS4gFI+17KTW4Szciz4oUibBF0IIIcQpSUJpI9D4mVEc1pNaR/nvMwgdeEM9lUgIIYQQommRUNoITnaaJtXjxrZ3Nf6t+tZTiYQQQgghmhYJpc1A1bpvCex6gcxBKoQQQohTloTSJk5VVcqWfUBo/2t9XRQhhBBCiAYjobSRaLRaVMVz3M8rX/o+Ae2HyAAnIYQQQpzSJJQ2Eq059LinhSpb+gG2vX8Tcd4DDVQqIYQQQoimQUJpI9EFhOOpLq3z8pati7DtWUXcde/LtzcJIYQQ4pQnaaeR6ALC8VjqFkoVu4Xi758n9qo3JZAKIYQQ4rQgiaeR6ALC6lxTWvjlJCIvfBStKbCBSyWEEEII0TRIKG0kusAI3FXFx1yuauOPqB4XgR2HNkKphBBCCCGaBgmljcQQnoi7POeoy5T//iEVf3xKzJhXG6lUQgghhBBNg4TSRmKISMZVcuCIj5ct/QD7gQ0k3PI5WqO5EUsmhBBCCOF7el8X4HShD43HXXb4mtLKNV9jy/iLuOuny8AmIYQQQpyWJAE1Eo1Oj+pxAWDP3oKzIAMAy9bFVPw1i7hr35NAKoQQQojTltSUNiKtKRjL5l8oXTgVfXA0rtJsdEGRJEz4FI3e4OviCSGEEEL4jITSRhTSdyxF858h+b4f0ZlDUBUPGq3O18USQgghhPA5CaWNKKjbhQR1u9D7uwRSIYQQQoga0olRCCGEEEL4nIRSIYQQQgjhc6dcKHW5XHz33Xe4XC5fF0WcBNmPpwbZj82f7MNTg+zHU8Opvh9PyVD6/fffn7I77HQh+/HUIPux+ZN9eGqQ/XhqONX34ykXSoX4v/buPSiq8v8D+JsVFxBBQFgUNUFLUIQESwhKjTRvwBQ2jppllGYzlqZTk7dmopkuyjhkOZaNGgPqeJlyuqjR4Cqo5CWduCk33QHTFUi8cBPZZX9/MHu+7Y8Fd+Eczh55v/5yz+VznvN8fODDOc85S0RERMrDopSIiIiIZMeilIiIiIhkJ1tRevjwYbkO3SNStVvK/lBqbKkotT+UGlsqSu0PpcaWilL7Q6mxpaLU/lDi73Qpyd1uFqV2UuJ/YKXGlopS+0OpsaWi1P5QamypKLU/lBpbKkrtDyX+TpeS3O3m7XsiIiIikh2LUiIiIiKSHYtSIiIiIpKdsy0bmUwmAEBTU5NoB25raxM1npk5ZlNTE1Qq8WtuqdotVVylxpYyj0rsD6XGZh6VH1upP1MZ2xLz2LuxlfgzFZC+3eZ6sjNOpodtAeDff/9FcnKyOC0jIiIioj7nhx9+gK+vb6frbSpK29raUFdXBzc3Nzg5OYnaQCIiIiJ6dJlMJjQ3N8PHx6fLK7w2FaVERERERFLig05EREREJDsWpUREREQkOxalRERERCQ7m14J1Ruam5vx008/oaysDGVlZWhoaMDKlSsxbdo0i+0KCwuxbt06qzFSU1MREhKCtLQ0aLXaTo+Vnp6OwYMHC59bW1uxZ88eHD9+HA0NDQgMDMSiRYsQEREhzsn1IXLl0ZZ4ZBsxc/hfN27cwO7du3Hp0iXU19fDz88PU6ZMwcsvvwxXV1dhO45HcciZR45H8UiRx4qKCmRmZuLy5csAgODgYCQnJ2PUqFEW+3EsikeuPCptLDpMUXrv3j3s27cPfn5+CAoKQmFhYZfbJyQk4IknnrBYNnToUADArFmzMGHCBIt1JpMJ27Ztg0ajsShIAeCrr77C6dOnkZiYiICAABw7dgwpKSn47LPPEBoa2vOT60PkzOPD4pFtxMyhWW1tLVavXg13d3fMmTMHHh4eKCkpwd69e3HlyhVs2LBB2JbjURxy59HWmNQ1sfNYUVGBjz76CL6+vliwYAFMJhMOHz6MtWvXYvPmzRg+fLiwLceieOTMoy3xHIXDFKU+Pj7IyMiAt7c3ysvLsXr16i63Dw0NRWxsrNV1ISEhHar/4uJitLS0YOrUqRbLy8rKkJubi+TkZCQlJQEA4uLi8O677yI9PR2pqandP6k+SK482hKPbCNmDs2OHz+OxsZGbNy4ESNHjgQAzJw5EyaTCVqtFg0NDRg4cCDHo4jkzKM9MalrYudxz549UKvVSE1NhaenJwBg6tSpeOedd5CRkSFcVeNYFJdcebQ1nqNwmDml/fv3h7e3t137NDU1wWg02rRtTk4OnJycMGXKFIvlp0+fhkqlwsyZM4VlarUa06dPR0lJCWpra+1qU18nVx67G486kiKH5m/z8PLyslju7e0NlUoFZ+f2v485HsUjZx7tiUldEzuPxcXFePLJJ4VCBmgvmEJDQ3H+/Hk0NzcD4FgUm1x5tDWeo3CYK6X22rJlC5qbm6FSqRAaGork5OQOl6bNDAYDTp06hZCQEPj7+1usu3r1KoYNG4YBAwZYLB8zZgwAQKfTwc/PT5qTINHy2J14JA5b+jwsLAw//vgjvvnmGyxcuFC47Xv06FHEx8cLcxE5HuUjZh7tiUnielift7a2wsXFpcN+Li4uMBgMqKysREhICMeizMTKo63xHIXiilJnZ2fExMTgqaeegqenJ6qqqnDo0CGsWbMGmzZtwujRozvsc/HiRdTX11u95VtXV2f1rxfzslu3bol+DiR+HrsTj3rGnj6fOHEiFi1ahAMHDuDs2bPC8nnz5uG1114TPnM89j4p8sjx2Pts7fPhw4ejtLQURqMR/fr1A9Be4JSVlQH43xjjWJSH2HlU2lhUXFE6duxYjB07VvgcFRWF2NhYvPfee8jIyEBKSkqHfXJycuDs7Ixnn322w7oHDx6gf//+HZar1WphPYlP7Dx2Jx71jL19rtFoMH78eMTExMDDwwN//fUXDh48CG9vb8THxwPgeJSDFHnkeOx9tvb57NmzsW3bNnz99deYO3cuTCYT9u/fj9u3bwP43xjjWJSH2HlU2lhUXFFqTUBAAKKjo5GXl2fxVwPQ/hqGs2fPIiIiwmLuhZlarUZra2uH5eaEmgcgSa8nebQ3Hkmjsz7Pzc3F1q1bsX37dvj6+gIAYmJi0NbWhvT0dEyePBmenp4cjw6ip3m0JyZJx1qfz5o1C7W1tTh06JDwyr3HH38cSUlJOHDggDAFg2PRcfQkj7bGcxSPRFEKAL6+vjAYDGhpabGYA3PmzJkun9b28fGxehvC/NeGtdcOkXS6m0d745F0rPX5kSNHMHr0aKGQMYuKisKxY8dw9epVTJgwgePRgfQkj/bEJGlZ6/PXX38dSUlJqKyshLu7OwIDA5GRkQEAGDZsGAD+bnQ03c2jPfEcgcM8fd9TN2/ehFqt7vDXwYkTJ+Dm5oZJkyZZ3S8oKAjXr18Xnio1Ky0tFdZT7+luHu2NR9Kx1ud37txBW1tbh20NBgMACE+Ecjw6jp7k0Z6YJK3O+nzgwIEIDQ1FYGAgAODvv/+Gr6+v8H5LjkXH0t082htPboorSu/evdthmU6nw7lz5xAREQGVSmWxbX5+PqKjozvt+NjYWLS1teH3338XlrW2tiI7OxvBwcF8ulAiYufRnngkDnv6PCAgAFeuXMH169ctts/NzYVKpRJ+oHI89j4p8sjx2Pt60ucnT55EeXk5EhMThe04FuUhdh6VNhYd6vb9b7/9hsbGRuGWwblz54R/x8fHw93dHZs2bYJarUZISAi8vLxQVVWFrKwsuLi4YPHixRbxTp48CaPR2OUt3+DgYMTGxiIjIwN3797F0KFDodVqUVNTgxUrVkh2ro8yOfJoTzx6OLFzmJSUhAsXLmDNmjXCNwGdP38eFy5cwIsvvijcCuR4FJdceeR4FJeYeSwqKsK+ffsQEREBDw8PlJaWIjs7G5GRkUhMTBS241gUnxx5VNpYdDKZTCa5G2H21ltvoaamxuq6HTt2wN/fH7/88gtycnKg1+vR1NSEQYMGITw8HAsWLEBAQIDFPh988AGqq6uRnp7e5UTeBw8eYPfu3Thx4oTF9/tGRkaKen59hRx5tCcePZzYOQTavyFm7969uHr1Kurr6+Hv74+4uDjMnTvXIq8cj+KRK48cj+ISM496vR7ffvstrly5gubmZiF/L730Uoen7TkWxSVHHpU2Fh2qKCUiIiKivsmxJhMQERERUZ/EopSIiIiIZMeilIiIiIhkx6KUiIiIiGTHopSIiIiIZMeilIiIiIhkx6KUiIiIiGTHopSIiIiIZOdQXzNKRI+ehIQEu7bXaDTYuXMn1q5di6KiIuGbTpQgOzsbW7ZsET67urri4MGDwufq6mosWbIE48ePxxdffCHacdPS0qDVavH5558jLCzMpn0KCwuxbt06xMXFYdWqVcLyn3/+GTt27BA+m/NBRCQ1FqVEJKm4uLgOyy5fvgy9Xo+goCAEBQVZrPP09OytpknGfF5qtVqUeOavJ/z1119FideVESNGCDnTarWSH4+IyIxFKRFJ6r9X4czS0tKg1+sRHR2NhQsXdrpfS0sLBg8eLHUTRdfVeUlh8eLFeOWVV+Dn59fjWJGRkcJ3m7MoJaLexKKUiBySRqORuwmK4ePjAx8fH7mbQUTUIyxKicghdTanNCEhARqNBt9//z0OHjwIrVaLW7duQaPRYO7cuZg2bRoAID8/H/v370dFRQVUKhUmTZqEJUuWWJ0eYDQakZWVBa1Wi6qqKhiNRgwbNgwvvPAC4uPj0a9fP9HPr6mpCbt378aff/6Ju3fvwt/fHzNmzEBiYiJUqvZnUM3zPv977mb/nevZ1ZzSyspKZGZmoqioCG1tbQgKCsK8efNEm1pARCQWFqVEpEgbN25EQUEBwsLCMGTIEBQVFQkPGbm5uSE1NRXBwcGIjIxESUkJjh8/jurqanz55ZdwcnIS4rS0tODTTz9FQUEBPDw8EBwcDLVajbKyMuzYsUMoDM2FohhaW1uxfv166PV6hIeHw2AwID8/Hzt37oROpxOmPHh5eSEuLg55eXm4f/++xfxcW+belpeXY/369WhubsbIkSMxcuRI3LhxAykpKZg1a5Zo50NEJAYWpUSkODU1NXBzc8P27dsxaNAgAEBBQQHWr1+PzMxMoeh7+umnAbRflfzwww9x6dIlFBYWIjw8XIi1a9cuFBQU4LnnnsPy5cvh7u4u7JOamoqzZ88iKytL1CKutLQUgYGBFu3X6/VYs2YNtFotoqOj8cwzz2DEiBFYtWoVioqKcP/+favzcztjMpmQlpaG5uZmzJ8/H6+++qqw7vDhw/juu+9EOx8iIjHwPaVEpEhLly4VCjoACA8Px6hRo1BXV4eJEycKBSkADBgwADNmzAAAFBUVCcvv3LmDP/74A76+vli5cqVQkJr3WbFiBZydnXHkyBHR2//mm29atH/o0KGYP38+gPaisacKCwtx7do1DBkyRIhrNmfOHAQHB/f4GEREYmJRSkSK4+zsjPHjx3dYPmTIEABAREREp+vq6uqEZYWFhTAYDJg4cSJcXFw67OPt7Y2AgABUVlaipaVFrObDw8PDahsnT54MoP2VWW1tbT06RnFxMQAgNjbW6pxY87GIiBwFb98TkeJ4eXlZLbRcXV0BwOprpNzc3AC0z+c0q6mpAQBkZWUhKyury2M2NDRYLVy7o7NXN7m7u8Pd3R2NjY1oaGjo0TtbzcV3Z8fi2w2IyNGwKCUixXnYQ0e2PpRkvho5atQoBAYGdrmtszN/XBIRSYk/ZYmoz/L19QUAjBs3DsuWLeu149bW1lpd3tTUhMbGRqjVaov5rd1hfm9pZ8cyXyUmInIUnFNKRH1WeHg4VCoVzp07B4PB0GvHra+vR35+foflubm5AICQkBCL6Qnmq7RGo9HmY4wbNw4AkJeXZ3V+6smTJ+1qMxGR1FiUElGfNXjwYEyfPh01NTVITU3F7du3O2xz48YNnD59WvRj79q1C/fu3RM+37x5E/v27QPQ/nT8f5mvev7zzz82xw8LC8Pw4cOh1+uxf/9+i3VHjx5FSUlJd5tORCQJ3r4noj5t6dKlqK6uRl5eHi5evIigoCD4+fmhpaUFVVVV0Ov1iIqKQmxsrGjHDA4OhsFgwNtvv43w8HAYjUbk5+ejpaUFU6dORUxMjMX2UVFRKCoqwoYNGxAeHg4XFxd4enrijTfe6PQYKpUK77//PjZs2IC9e/ciLy8Pjz32GPR6PSoqKjB79mxJXnVFRNRdLEqJqE9zcXHBJ598gpycHBw7dgw6nQ7l5eXw9PSERqPB888/L/rrk/r374+UlBRkZGTgzJkzuHfvnsXXjP5/CQkJaGhoQG5uLvLy8mAwGKDRaLosSoH24jc1NRWZmZkoLi7GzZs3ERgYiI8//hiurq4sSonIoTiZTCaT3I0gInoUZGdnY8uWLViwYAEWLlzYa8fduHEjTp06hc2bN2PMmDGixU1ISIBGo8HOnTtFi0lE1BleKSUiEtmZM2dQXV0NtVqN5cuXS3oso9EInU4HJycn4QsCeuLixYvIyckRoWVERPZhUUpEJDKdTgedTgdXV1dJi9KtW7fi8uXLuH79OiIjI3v0sn2za9euQavVitA6IiL78PY9EZFCmb/TPiIiAsuWLYOXl5e8DSIi6gEWpUREREQkO76nlIiIiIhkx6KUiIiIiGTHopSIiIiIZMeilIiIiIhkx6KUiIiIiGTHopSIiIiIZMeilIiIiIhkx6KUiIiIiGTHopSIiIiIZPd/5EQaADwdJtwAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cbvs.plot(cbv_indices=np.arange(1,5));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " As further examples, below we will request both Kepler and K2 CBVs." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:26.375717Z", "iopub.status.busy": "2023-11-03T14:13:26.375546Z", "iopub.status.idle": "2023-11-03T14:13:27.201241Z", "shell.execute_reply": "2023-11-03T14:13:27.200805Z" } }, "outputs": [], "source": [ "cbvsKepler = load_kepler_cbvs(mission='Kepler', quarter=8, module=16, output=4)\n", "cbvsK2 = load_kepler_cbvs(mission='K2', campaign=15, channel=24)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Converting the CBVs to a Lightkurve `DesignMatrix`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can access individual CBVs directly by referencing the column in the table (i.e. `cbvs['VECTOR_#']`) however a better way to use the CBVs for correcting your light curves is to convert the CBVs to a LightKurve DesignMatrix:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:27.203584Z", "iopub.status.busy": "2023-11-03T14:13:27.203275Z", "iopub.status.idle": "2023-11-03T14:13:27.207922Z", "shell.execute_reply": "2023-11-03T14:13:27.207613Z" } }, "outputs": [ { "data": { "text/plain": [ "10.1.1.SingleScale DesignMatrix (18900, 8)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cbv_designmatrix = cbvs.to_designmatrix(cbv_indices=np.arange(1,9), name='10.1.1.SingleScale')\n", "cbv_designmatrix" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:27.209728Z", "iopub.status.busy": "2023-11-03T14:13:27.209582Z", "iopub.status.idle": "2023-11-03T14:13:27.231674Z", "shell.execute_reply": "2023-11-03T14:13:27.231354Z" } }, "outputs": [], "source": [ "cbv_designmatrix.plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aligning the CBVs with your light curve" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The CBVs obtained from MAST may not have the same cadences as the your target. A method called `align` will allow you to align your CBVs to the cadence numbers in your target lightcurve object. The method will use the cadence number (i.e. [lc.cadenceno](https://docs.lightkurve.org/reference/api/lightkurve.KeplerTargetPixelFile.cadenceno.html?highlight=cadenceno)) to perform the alignment. Only cadence numbers that exist in both the CBVs and the light curve will have values in the returned CBVs. All cadence numbers that exist in the light curve but not in the CBVs will have NaNs returned for the CBVs on those cadences and the GAP set to True." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:27.233448Z", "iopub.status.busy": "2023-11-03T14:13:27.233349Z", "iopub.status.idle": "2023-11-03T14:13:27.237037Z", "shell.execute_reply": "2023-11-03T14:13:27.236781Z" } }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Take a cut of the LC loaded above\n", "lc_short = lc[501:1501]\n", "# Take a different cut of the CBVs\n", "cbvs_short = cbvs[0:1000]\n", "# These cuts do not overlap\n", "np.all(lc_short.cadenceno == cbvs_short.cadenceno)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:27.238549Z", "iopub.status.busy": "2023-11-03T14:13:27.238464Z", "iopub.status.idle": "2023-11-03T14:13:28.670565Z", "shell.execute_reply": "2023-11-03T14:13:28.670262Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "The MultiScale CBVs do not appear to be well aligned to the light curve. Consider using \"interpolate_cbvs=True\"\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Align the cuts\n", "cbvs_aligned = cbvs_short.align(lc_short)\n", "# They now fully overlap\n", "np.all(lc_short.cadenceno == cbvs_aligned.cadenceno)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Interpolating the CBVs to an arbitrary light curve" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above `align` method will only keep cadences that line up exactly between the CBVs and the light curve based on the cadence numbers. What if you have a light curve with cadence times not exactly lining up with the CBV cadences? For example, what if you want to use the 2-minute CBVs for cotrending against the 30-minute FFIs? A more general method is `interpolate`, which uses PCHIP interpolation to generate CBVs at the cadence times of an arbitrary light curve. If the light curve has cadences past either end of the cadences in the CBVs then one must extrapolate. An optional argument, `extrapolate`, can be used to also extrapolate the CBV values to the light curve cadences. If `extrapolate=False` then the exterior values are set to zeros, which will probably result is a very poor fit." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:28.672514Z", "iopub.status.busy": "2023-11-03T14:13:28.672390Z", "iopub.status.idle": "2023-11-03T14:13:37.641020Z", "shell.execute_reply": "2023-11-03T14:13:37.640636Z" } }, "outputs": [], "source": [ "# Get a light curve from a TESS FFI\n", "from lightkurve import search_tesscut\n", "search_result = search_tesscut('HAT-P-11', sector=14)\n", "tpf = search_result.download(cutout_size=20)\n", "target_mask = tpf.create_threshold_mask(threshold=15, reference_pixel='center')\n", "ffi_lc = tpf.to_lightcurve(aperture_mask=target_mask)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:37.643519Z", "iopub.status.busy": "2023-11-03T14:13:37.643163Z", "iopub.status.idle": "2023-11-03T14:13:37.863848Z", "shell.execute_reply": "2023-11-03T14:13:37.863387Z" } }, "outputs": [], "source": [ "# Get the Single-Scale CBVs for this light curve\n", "cbvs = load_tess_cbvs(sector=ffi_lc.sector, camera=ffi_lc.camera, ccd=ffi_lc.ccd, cbv_type='SingleScale')" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:37.866761Z", "iopub.status.busy": "2023-11-03T14:13:37.866301Z", "iopub.status.idle": "2023-11-03T14:13:37.894273Z", "shell.execute_reply": "2023-11-03T14:13:37.893858Z" } }, "outputs": [], "source": [ "# Interpolate the CBVs to the FFI cadence times\n", "cbvs_interpolated = cbvs.interpolate(ffi_lc, extrapolate=False)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:37.897016Z", "iopub.status.busy": "2023-11-03T14:13:37.896762Z", "iopub.status.idle": "2023-11-03T14:13:37.899976Z", "shell.execute_reply": "2023-11-03T14:13:37.899579Z" } }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# All cadence times agree\n", "np.all(ffi_lc.time.value == cbvs_interpolated.time.value)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2023-11-03T14:13:37.901796Z", "iopub.status.busy": "2023-11-03T14:13:37.901680Z", "iopub.status.idle": "2023-11-03T14:13:37.923353Z", "shell.execute_reply": "2023-11-03T14:13:37.923023Z" } }, "outputs": [], "source": [ "# The CBVs are now interpolated to this FFI-derived light curve\n", "import matplotlib.pyplot as plt\n", "_, ax = plt.subplots(2, sharex=True)\n", "ax[0].plot(ffi_lc.time.value, ffi_lc.flux, '.b')\n", "ax[0].set_title('HAT-P-11')\n", "cbvs_interpolated.plot(cbv_indices=np.arange(1,4), ax=ax[1]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Some notes about interpolating Cotrending Basis Vectors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Currently, only TESS 2-minute and Kepler/K2 30-minute CBVs are archived on MAST\n", "- 20-second CBVs will be available for the TESS extended mission beginning with Sector 27. Each set of 20-second CBVs will be for the entire field of view.\n", "- FFI CBVs are also being developed and will begin to be exported soon.\n", "- The CBVs are generated to account for systematics at a specific cadence. They will not necessarily properly represent systematics at a different cadence, but in some cases can still be beneficial.\n", "- Please remain conscious of the Nyquist frequency and aliasing when interpolating CBVs." ] } ], "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.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }