Home.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0.2\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import readligo as rl\n",
    "from gwpy.table import EventTable\n",
    "from gwpy.timeseries import TimeSeries\n",
    "\n",
    "import matplotlib.mlab as mlab\n",
    "from scipy.interpolate import interp1d\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import gwpy\n",
    "print(gwpy.__version__) # 2.0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['GW150914',\n",
       " 'GW151012',\n",
       " 'GW151226',\n",
       " 'GW170104',\n",
       " 'GW170608',\n",
       " 'GW170729',\n",
       " 'GW170809',\n",
       " 'GW170814',\n",
       " 'GW170817',\n",
       " 'GW170818',\n",
       " 'GW170823']"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eventtable = EventTable.fetch_open_data(\"GWTC-1-confident\")\n",
    "eventtable['commonName'].tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Select a GW events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBH\n",
    "event_name = 'GW150914'\n",
    "# event_name = 'GW151012'\n",
    "# event_name = 'GW151226'\n",
    "# event_name = 'GW170104'\n",
    "# event_name = 'GW170608'\n",
    "# event_name = 'GW170729'\n",
    "# event_name = 'GW170809'\n",
    "# event_name = 'GW170814'\n",
    "# event_name = 'GW170818'\n",
    "# event_name = 'GW170823'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Two ways to load data.\n",
    "\n",
    "- Way 1\n",
    "\n",
    "Download from GWOSC host..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tevent = eventtable[eventtable['commonName'] == event_name]['GPS'].tolist()[0]\n",
    "sample_rates = 8192 # [Hz]\n",
    "tag = 'CLN' #'C00'\n",
    "\n",
    "# Download 32sec open data from GWOSC (It may takes a few minutes...)\n",
    "data_H1 = TimeSeries.fetch_open_data(ifo='H1', \n",
    "                                     start=tevent-16,\n",
    "                                     end=tevent+16,\n",
    "                                     sample_rate=sample_rates*2, # 4096Hz or 16384Hz\n",
    "                                     tag=tag)\n",
    "strain_H1 = data_H1.value[::2]  # 16384Hz -> 8192Hz\n",
    "data_L1 = TimeSeries.fetch_open_data(ifo='L1', \n",
    "                                     start=tevent-16,\n",
    "                                     end=tevent+16,\n",
    "                                     sample_rate=sample_rates*2, # 4096Hz or 16384Hz\n",
    "                                     tag=tag)\n",
    "strain_L1 = data_L1.value[::2]\n",
    "try:\n",
    "    data_V1 = TimeSeries.fetch_open_data(ifo='V1', \n",
    "                                         start=tevent-16,\n",
    "                                         end=tevent+16,\n",
    "                                         sample_rate=sample_rates*2, # 4096Hz or 16384Hz\n",
    "                                         tag=tag)\n",
    "    strain_V1 = data_V1.value[::2]\n",
    "except ValueError:\n",
    "    print('No V1 data for {}'.format(event_name))\n",
    "\n",
    "# both H1 and L1 (and V1) will have the same time vector, so:\n",
    "time = data_H1.times.value[::2]    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Way 2\n",
    "\n",
    "Use local downloaded data. (`./hdf5`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/herb/anaconda3/lib/python3.7/site-packages/h5py/_hl/dataset.py:313: H5pyDeprecationWarning: dataset.value has been deprecated. Use dataset[()] instead.\n",
      "  \"Use dataset[()] instead.\", H5pyDeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "tevent = eventtable[eventtable['commonName'] == event_name]['GPS'].tolist()[0]\n",
    "sample_rates = 8192 # [Hz]\n",
    "\n",
    "if event_name == 'GW150914':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1126259447-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1126259447-32.hdf5'\n",
    "elif event_name == 'GW151012':\n",
    "    fn_H1 = './hdf5/-H1_GWOSC_16KHZ_R1-1128678885-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1128678885-32.hdf5'\n",
    "elif event_name == 'GW151226':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1135136335-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1135136335-32.hdf5'\n",
    "elif event_name == 'GW170104':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1167559921-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1167559921-32.hdf5'\n",
    "elif event_name == 'GW170608':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1180922479-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1180922479-32.hdf5'\n",
    "elif event_name == 'GW170823':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1187529241-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1187529241-32.hdf5'\n",
    "elif event_name == 'GW170729':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1185389792-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1185389792-32.hdf5'\n",
    "    fn_V1 = './hdf5/V-V1_GWOSC_16KHZ_R1-1185389792-32.hdf5'\n",
    "elif event_name == 'GW170809':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1186302504-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1186302504-32.hdf5'\n",
    "    fn_V1 = './hdf5/V-V1_GWOSC_16KHZ_R1-1186302504-32.hdf5'\n",
    "elif event_name == 'GW170814':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1186741846-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1186741846-32.hdf5'\n",
    "    fn_V1 = './hdf5/V-V1_GWOSC_16KHZ_R1-1186741846-32.hdf5'\n",
    "elif event_name == 'GW170818':\n",
    "    fn_H1 = './hdf5/H-H1_GWOSC_16KHZ_R1-1187058312-32.hdf5'\n",
    "    fn_L1 = './hdf5/L-L1_GWOSC_16KHZ_R1-1187058312-32.hdf5'\n",
    "    fn_V1 = './hdf5/V-V1_GWOSC_16KHZ_R1-1187058312-32.hdf5'\n",
    "\n",
    "\n",
    "# read in data from H1 and L1, if available:\n",
    "strain_H1, time_H1, chan_dict_H1 = rl.loaddata(fn_H1, 'H1')\n",
    "strain_L1, time_L1, chan_dict_L1 = rl.loaddata(fn_L1, 'L1')\n",
    "strain_H1 = strain_H1[::2]  # 16384Hz -> 8192Hz\n",
    "strain_L1 = strain_L1[::2]\n",
    "try:\n",
    "    strain_V1, time_V1, chan_dict_V1 = rl.loaddata(fn_V1, 'V1')\n",
    "    strain_V1 = strain_V1[::2]\n",
    "except NameError:\n",
    "    pass\n",
    "\n",
    "time_H1 = time_H1[::2]\n",
    "# both H1 and L1 will have the same time vector, so:\n",
    "time = time_H1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydeZxO1R/HP98Zg7Eb+5aySyIS2ZeEEkqFX0SKylIUqVQa0UKLSAohEhVZSijZ92xlG0vW7Pu+jJnz++PeZ+buz93vMzPn/XrNa+4999xzv/d57nO/Z/kuxBgDh8PhcDhKooIWgMPhcDiRCVcQHA6Hw9GEKwgOh8PhaMIVBIfD4XA04QqCw+FwOJpwBcHhcDgcTbiCiGCIKAsRnSSiRQZ18hFRxNoqE9ElIiqpUb6fiKprlDcgoq2S/RgiGkxEO4kogYj2ENEEIiosqRNFRK8Q0Q5JnRFElN2irCWJ6EUb91iUiDZaPU+jHc3PSlGnARG1cHotryGipkS0TvzedhPRkpDcRNSBiDYr6m8nolmSfSKi00RUTdyPIaLnpc+GWP4qER0Xr7OTiBZKjg0Qn4U9RDSLiAqK5ZWJaDER7SKibUTUTtFmQSKaT0Sv6txbQyK65fQzSgtwBRHZPAngLwBViKhC0MIExFwA+QHcwxirAOAuAEcBLCai0PM7HkB9AHUldW4AWENEsRaudTuAHlYFZIwdZYyplJ1HNATwkE/XsgURdQAwAsDTjLGKjLFyAN4GMJWIagP4A0BlIooT65cEkBlAfSKKEZupAiARwGYiehbAIQAfaVyuBoCnxOtUZIw1E9tsCuG7rM0YKwvgJIDPxHPmAhjNGCsP4GEAI4joTvG8+QC2A2ikc2/FAAwEEG3z40lbMMb4n8k/ANUBfA1gkLhfD8AUADMAlPDgemsBtAUwBsBIsSybKMM+ALsBvA+AQVD2BwHUkZw/DUBXANlFOQ9A+KEtAVBIrMMA9IHwozgE4H3J+d0A7ASQAGAFgLvF8gIApgPYJsrRWyyPFuXZC2APgOEArgAoqXFv+wFU1yhvAGCruN0SwD8ASFEnCkBJ8XrVICiMWI22FgN4UaP8UfF+EwBsBdAMQBPx87sh3nNOUcauogzPAigPYKl4z/8B+AYAibJcEtseBGAChJfgHgCbAVTUkEH3sxKPjRCvvx/ARgAVALwB4BSAcwBmie08B2CX2Ma/AB7T+UzXABgn1t0HoK3keB8AW8TnYyaAvGJ5K/Hz2SN+Nj0l9zhSrLtA476OAbhXQ46CAHKL2xsAtBa3XwDwAYA/ATQSy14FMFHv2ZCU/QvgB/G+fgdQVSz/DMDHknpVAZyF0Nm4BSBKcmwmgNcU7U4E8KqiLAbAQgDlALCg30d+/AUugCc3JTyIORVlZRT7HQFMk+yPgvBC3gDgYYO2SyJVQbwLoBiADgAecfkeqokvgiwA6gK4ACAHgKEA5ogPaxSA0aGHFUA8gC/E7VgAJyC86F6EoGSiILzQlkF8cUJQEIMln9sV8Z7qiS+MYuKxdgA2itu/AXhA3M4FYDWAphAUyjoAOcRjr4ntl9S4PzMKYiiAjyTH3oLw8t4J4DiA2hB6id/rfIavSr9jxbWbitt1Q3WgeAGJ9RYByCPufwehVxy678sAKkKtII4CKC7ufwJgnIYMup8VBMU4F0Am8djk0Ocgtj9K3I6DoGDyi/tdAazT+UyTJPf8CIA94vaTEF6GoWv1AzBVfFb2QPzdAGgM4IREhksQeufKa5UN1RP3K0m+s4MS2YcCGCFuzwVQR/y+Qvc5D8ATes+G5BlfAOG3EgWgOwTFnU28p/6SugWQ+jv5F0B3iXwnAHyquJaWgvgCwJOh342bv/dI/UuvU0wNAMwhoqwAQETdIPR4QERx4nzx8FBlImoGoBRjrBaEXtNoyfSFEQsg9KT7QejluUkvAD8yxm4AWAVBWXQE0BzC8DiRMZYMYJjknEkAHieiaAhD5wWMsUuMsTEAPgfQWax/F4RRRYixAMAYOwmhh1oYwucwjTF2RDz2A2Osujiv3xTAKCLaCeElVwzCy7o5gG8YY5fFdj+DME1glzjp+YyxIUycSoCg6GLE/1aZBmA8EY0EUBTC6ECPUYyx8+J2JwB7iagnhHuLhfxzDPEbY+w/cXsrhM9Tie5nxRj7FUBvAO2JaAiE6Q7VdRhjZyGMausR0RsQXpB66y77GWN/aMjUFsL3uVX8PrtDUNzJAO4GUI6I+gN4XdH2IsbYao3rKL+z7ZLvbDqETg4g9PYbElEWCL37tRCUQgsiygThefoDBjDGrjHGmjPGNjHGkhljYyEorqpilWRJdel2Gwi/k90QOnlzxfN0IaKnhEuyH43qpTfSpYJgjP0E8eUtfrFdIfSUwBg7y4T54vaSU6pBfBgZY0chzFfeRkQvEdE08W+kxqW6Q+hZPQnjl4wlxLnZ9hAe4mMQeqT5AfSE0MtkWucxxvZDmDppgtSeIYjodQhK4CyE6Y9fFadekzYj/o+G5IdORNFEdB+EFzJBWBMI/fDvgziSUcimKacFVkF4eckgorshjHYA4cXSUGet4WEIU2MyGGNvQpi7XibW2WGwoH1Usj0fwrOUAKEXfVjnHK3PUyUGdD4r8Zn9BcJ0148AvtVqQFyX+hvCy/5PAEN0rmUkUzSAeMl3WRnAI0RUAMKzdBeE72GAor2j0GYbgBxEVEnjWBPJ9moAdwB4DMByxlgSYywBghJqC2GkcF6jjRSIqJy43iHlFoT1jPMA8knKC4llgKAsHmKMlWOMPQHhWTpgdC0ATwBoHloMF6+/k4juCHNemiZdKggAYIwNA3AdwtRKK0lPTYvMkP+ArgDIwhgbyRjrIP69pHHeYggv3kEQFJJbdAWwjTGWjzFWhDFWBMK8ZzkAFwH0EXtegDB/LGUChHnpOyG8AAFhUXMSY2wOhB9JXYT/7v8A8FTI8gPAUxCG4ZchvDBeAwAiygZhyqkSBMXzPBHlFs/pBaGXb5epAK6KFknZxOvdAWGkBABgjG2EMA00lYjyiHViiOgDCD/8SdIGReuYfQAKMMZmijLGAcgL4cWpXHxMFs/LBEFZvc8Y+xPCVFBx2P8NGX1WzQDMFjs6ByG8WEPXkcpYG8BBcYS4UTzPqjzzALwY+uwA9IWw1lFZ3P9EHCk0VbQt7ZGnwBi7IrbxfciwgogyE9FHEEY7oXqJEJ7P9yDvsPwGYfppvgnZrwMYGVJG4sJ0HgiGHQshPL9xYt0eSP2NjgbwtHhOVQgGDrONLsQYa8MYKytRpBC395uQM82SKWgBvIKI6kHo/cyC8ALvZVA9EUBWyX4shN6bCsbYQQhz/WCMfQdhXto1iIggLNrJemyMsaNENBVCD+ssgJ1EdBnCy1HKzxB+AJ8xxkI9xc8BfExEfSH0emcCeImIxuvJwRibT0SjAawgwYz2GIAu4uFOEKbh9kDosX3NGFtGRCsAlAawiYiuQVgMv2JwuwuJSDoFtUKUPSQDI6LmECxgtotyJEPoVe+TnNcFwktptfj5RUN4wdRijF1V3BcTzRd/lEwjDmaM/SfKkkOceqiuOO8WEX0IYDkRnYCw+DwHQq+9m8E96jEZ+p/VOAhTYG0hzI8vAPACEX0BYUpvChEVBfA8gK6iwjslfi5xRNSJMTbFpBzfQlB264koGcJaQVcIz/92APuI6CQEQ4wTRDQwXIOMsYniyPc7IsoPQbFshTBtJv1cF0LovCyUlP0G4bcaVkEwxg4RUXcIMwWZIIz8W4pKagERfQVgrfjc7EHqKP8FCJ/vGxBGw88yxs6Fu15GhFLfIekHIroHwjzzwxAWGr8DsI8x9pakTgMALzDGOogvoR6MsVZEVARCD7ksYywpAPE5HA4nIkivU0zZIFhA/CsutnWGMFTXhDG2AMAhIloJYajZmysHDoeT0UmXIwgOh8PhOCe9jiA4HA6H4xCuIDgcDoejSbqyYho2bBjLkydP+IoaZM2aFdevX3dZosiG33PGIKPdc0a7X8DZPR87duz0oEGDCmgdS1cKIk+ePOjevbutcxMSElChQsaKh8fvOWOQ0e45o90v4Oye4+PjdQ14+BQTh8PhcDThCoLD4XA4mnAFweFwOBxNuILgcDgcjiZcQXA4HA5Hk8CsmIhoFISQy5kgJOCZpzi+EULkUkAI/asVTZXD4XA4HhGIgpAm6BGjUq4movli3CSIsf3PMsZUuQA02uoOIS8DunTpgvr169uS6fTp00hISLB1blqF33PGIKPdc0a7X8C7ew5qBCFL0COGE74NqUk7KgPIS0SLIcS/H8gYW6vVkJhFaiwAjB07ltm1Bea20xkDfs/pH9/u98QJ4Phx4M8/gRo1gHr1vL+mDl7dc1AKIjOEnAYhrkDIvRziJoREJdOIqBqAGURUmvHIghwOJxKYOxdo3Vpelg5fT0EtUodL0LMLwE8AwBjbBOAMhLzHHA6HEyxnzqiVAwBs2+a/LB4TlILYBDE/rZigpyDk+X17Q8yoRkSFIKSEPOGzjBwOhyPn6lUgf37tY5UrAwvczDwcPIEoCK0EPQDyE9ESscpoANWIaDmEkUQ3MYcth8PhBMewYcbHW7TwRw6fCMzMlTGmlSO6kXjsCoC2/krE4XA4YYiPD1oCX+GOchwOh2OGzz8PWgLf4QqCw+FwwnHzJtCnT9BS+A5XEBwOJ1jWrRN8CSKZli3N150+3Ts5fIYrCA6HExx79wK1agEPPBC0JMb88Yf5uh06eCeHz3AFwUl77NoFVKsGXLwYvi4ncjlwAChbNnU/0kcRGRCuIDhpjwoVgM2bgdy5gaefDloajh0YA+64Q172wAPWeupBsmKFsef0xx/7J4uHcAXBSVso7dCnTAFOnQpGFo59unXTLn/wQeDWLX9lCUdSkrqsbl3hf7Nm2uf07++dPD7CFYRb7N8vWDpwvOOvv4ABA9TlBQv6LwvHPmPGAN98o3880r7Pr77SP/bLL8A77wCxsepj//3nnUw+wRWEUw4fBtq3B0qVAv73v6ClMc+nnwKTJwMnTwYtiTmuXwfuu0//+Pr1/snCsc+tW0CPHsZ1zp0DkpP9kccMvRQ+vV26pG7HxAjOcw0bqs97800vpfIFriCcctttwA8/CNszZwJNmwoPeKTz6qvCqKdQIeDFF4UF30ieqjlyxPh4zZpCEDVOZGP2xR8dDVy75q0sZrhwQV325ZfqMqnSCLFypevi+A1XEE6YMEFdtmgREBfnvyxmSUwEjh6Vl331lbDgG2lD+xAXLwJlyoSvlz8/0Lix9/Jw7LNxo/m6//zjnRxmef99+f5992lPJz35pHqtYv9+7+TyCa4g7NK0KfDss/rH163zTxazJCcDxYsDxdJY5PTcuc3XXbIksqYnOHJq1zZfNxIsmq5fl+/37atfN0rjdfrkk+7K4zNcQdhl0SLj47VqRd6Lqm/f8GsOder4I4tZNm+2fs7s2e7LwXFOaCrWLG+/HXwSnpEj5fvhXvhKJfHTT+7K4zNcQdjhhMnUFNHRwpROpKB82LVYvTr8IqKfVKtm/Zy2bYErV9yXheOM9u2tnxNp36PWKEHKpEm+iOEXXEHY4bHHzNddvNg7OaxAZL7umDHA9997J4tZnPQec+RwTw6Oc/bt0z+2VjPdvECjRu7L4iWdOqnLIs2vwwJcQVjl4kWhl60kXz7t+koTOb85d85eFMqnnnJfFquE661VrGh83OxIj+M9pUtrly9dKligde+ufXzDBs9ECsvu3fL9adPMnac0b5071x15AoArCKtoLZgyBhw8qD1fvncv8O+/3sulR9u29uPYK38gkUByMrBli2ACuWOHUBYdrV03Eg0FOHIaNBD+jxkDrFmjXSeoHvhnn8n3zU6RFS0q30/DDrRcQVjh2DF1WeihyZ4dqFpVeyhtxkTTK/7+O3ydNm20y8uXd1cWK2iFN3jjDWGqrEoVIGtWoezcOUFZPP+8un64EQgnWKSj7qgowbAjb151verV/ZNJilQxhTojZnjhBfn+r7+6I08A8F+QFSpUUJcpY67ccUdqnJZIQOtFCwijnsqVBVPCWbO0nX+CRMsGXmmTDgB58gjerFrhEB55JLKMBDIq2bNrl58+rS5buhTIlk1e9s8/+s+xl4wfn7qdK5f585Qj2qlT3ZEnALiCsIIyvPS+fdpWNloL08ePeyNTOLQ8QUPmt489lhqH/8UXtc/XG/Z7jfJz1VvjCUekhpC+cEGYlrTSM02rXL2qLhsyRLvu3Xdrm2JPnuyuTOFQmqjrTWPqcfasfD8tRFfQgCsIs2j9kLWGw4DQo33jDXlZkSLuyxQOLRPB0aP1LZpmzFCXBbl+IkWrt6lkzx51mV7U0CB57TVh5HP77UClSkFL4y16lmjK34cUrRFH167uyGMW5e9dbxSkh3Ktsl07Z/IEBFcQZlH+kPPlE37ketix+XYbpaln69bGPg5t26otMDp1SjvJ2rXWeiIpomZSkuCLMnx40JL4h9bU5bp1kb8+9PbbqdtFigA5c1o7X3l/adRgIsK/pQgmXGC7u+9Wl/k5H641Z2vGw/i999Rlfidrf+01+b5Z80JAe6QXCesQZ84ICuzll4OWxF+0zLyNovKG0DIl15ou9QrpbyVcoEgzpNHsh1xBmOHSJXWZGcezVavk+376FtidWtHr2fkZ8kDZw7YyGqtYEXj0UXnZjRvOZXICY0IMrAMHtI8TWQtil1Y4fFhddvmyuXPvv189+itZ0rlMdrDiZCpF6aeTBqMNcwVhhvh4+f7XX5s7TxmYzK+4LIwBEyfKy5YvN3++lmmsX1YkykizDz5ovY2ff5bv6zlh+cHRo0DmzOqgb0ruvTfywko45bbb1GVW5vKVQSX9GkGE+67MorRe2rXLnXZ9hCsIMygfGCt+DV984a4sZtCKQ1+vnvnztabH/JrLV8qpFVLdKlamqNymWDHzjl45cqSfxEdaUyrNmztv1+wIxAlSX6YsWey3c8898v1IC4RpAq4gzLB9u3zfSs6BsmXdlcUMSrNCO6aqyjaUCea9QuloaDc0ec+e8n0rIyi3sONBW7Nm+hhJLFmiLitXzno7I0bI95XezV4gHX1rWcZlILiCMMPSpanblStbO1fpjWzGXNMpyimxWrWst6GVFMVrlC9UJyagQ4fK9zt2tN+WXez2PnPkECzOnHLqVHD+N1re+R9/bL0dZRIrP/xGpHKWKOH99SIYriDCofTotdobL1kSGDgwdb9AAecyhcMt5zblVJVR8nY3UL5Qt22z35bSDl1rwdRLZs1ydr4bAd4KFkz1vzl1KlinvD//FPyDrKI0UJg+3R15/EL53EVyWl8NuIIIhzIOjFWHGUDfa9QLlA+kk2iyyjlTPW/rSEW52K4cWXnFtWvGIeEbNAAeekhYwDZKonP77drmnlaJjxeURaVKwKZNwjy+35Zdt99u7zyitKcUpBQvLt8fNSoYOWzCFUQ4vIgkuXOn+22GUFqOtGrlrL1585ydbxYv0ksqzQzffdf9a2hhFHuHMWHKct48oXf/5JP6PeuDBwUlbUfJ16+vXV69uuD0lTWrdwYUWspHGV/JCkovZC/Dv3jtM/Pdd9627zJcQRihdJBx8rKVzoHfeaf9dqyitKSwykMPyfeNEr84QWnO6sZQPFMmddm1a87bDYdVH5RjxwRloMfo0ebbOn9e6HWvWBG+bu/e5tu1QiiEt5TChd1r30pea6t4EbRSatSyf7/77XsIVxBGKNcflPb1VtCKROo2ShPJ8eOB/PndvYaWCa0XuCE3ETBnjrzM6x6iMtQzIEz7JSbqWzXlyyeM/C5f1p/C1HLWunlTHlRu/Xr9+GB6eOGgpwwrMWCA+9fwCmnUALc6cmk4LzVXEEYoe89WIzpKUVpDeDF1pVzYdCvAmXS6onPn4BPJW0E56vPaTFLLibJ4cWE0E26RNnt245FEfDzQt6+wfe2asKgfHS0oj7x5BRNZq9x7LzB4sPXz9FBGQQWADz903m4Q1kRu+OAAQFycfN8oUGGEwRVEUHgRI15q2nn2rP0QAUqUjmYffeROu3q4bQrslye1MqhhuXLWHbvy5dNXwO++K/gF1KypntM/f97adaQMGuSevb+V6TArVKki3/ejk2JH4ZrBDYXpE1xB6KHsCRlFQTWL1CKqSxfn7RmhjOTqBGUKRbd7QMuWyfft5n7Qo1kz+b4XsfkPH1YHNezTx57VG5A6UtDCC29rO8pMi5deku+7ZbXz/ffyfS8MDrZscb9NPdwIAOgDXEHooew1jxzpvE1lAnY3c9UqF1/t2JwHRcOG3ravdDpTDvndQCstplYaVLN8+ql3BgF6WA1prUSZJAdwr6efM6d8Ud3NabEQLVumbhspaDsojQYiKQy9AVxB6KH0vHWy/qCHm45L0ikHP5zxvMKLiJ3R0dqWNW6xb5/a6mryZOc5D/wKbyJFL+KsGbSmHt2MYPzkk/J9rUx1TpD26t32vFc6gWqt1UQgXEFo8dtv8v1Bg9xrWzqX6tQEVY/ff3e/Tb+czLyKmaRUmloh3O0i7XmG6NTJnbbdCJVhJUf6HXfYjwU1bJi6zM3RmlLheplbQyv/vJt07uxOO4wBU6Z45qEdmIIgolFEtJaINhDRw1aPe4rSvNXN6RovXt7KHL5Vq7p/jXfeMb6mXY4dk+9rhYh2g6eflu9bSUJvxM2basdHraRLdilUyFyCHT2mTBHs8F9/XRjpMCakkTWy98+Rwx1zZrcD3SlNn9000VX26J049mmhzLHu1mcTEyM82174bwDQ8CTyHiJqBqAUY6wWERUFsJqI5jPGks0cV7TVHUB3AOjSpQvq63mQhuH06dNISEgQdnLlkseAuftuIHTMDaRtu9Hu2LGpbUZHm25Tds9m+Pjj1HWUl192Z2Q1YoT7n4cWZcsCkybh9H//ISF0PTeudf26Ol5Q27bu3se33woeuFbyg9evL4xQ8+QRvufOnQUP55BcjRoJJr96qTBXr7bmi7Jjh/pzuHXL/e9Teo38+TXbt/xcA4LhgtfP4aRJwIIFqftr1lj3W1HyxBMAgNN33WX9ns3AGPP9D8AbAPpI9tcDuN3scb2/r7/+mtll586dqTtCPyv1z22kbU+c6G57W7eaPk12z2b47Tf5tXbssCioBtL2smVz3p4RJ0+yne3bp14vOdlZexcuqJ+Vxx5zR1YlycmMtW2rvl7oL1MmxhYvFupeucJYUlLKqbrf86lTjEVH67e5YIE52UaOVJ8bHe3whnUw8du0/FwzxliLFqlt9u7tUEgD3H63iO3sbN/edhPvvvvuBqbzTg1qiikzAKnZzRUAWSwc9w8v8jFLTf9ef91ZW/Pny/fvustZe0a0aCHf18o85wRlNjm3Ua5DOLVUUeacANRmnm5BBMyYAbz5pjBdoVzAvnJFGBUAwvSImQXy/PmNLemaN1evx2mhdc9GQQid4IWJMiD/HSlT1kYqPpjKBqUgEgFklezHArhh4bh3KH8wr73m/jWkUVFPnHDW1syZqdtOQ0xbxakJo3JhTRmi22uUjm1WSEpSB17r3Fk/SJ5bDB0qzL0vXixMz/XsKUxTZM5sr72oKON1sYfDLP/pBZ+791578oQjTx75vhdz7yFFG+ls3uz5JYJSEJsANAEAIioCoCCAwxaOe4fSikaZsMQNlL075UKtFb75JnVbK0mLlzz7bLDnu4FdXxStF+CkSe55r4fj9tuFdaAvvtD2P7BC06ZAjRr6x42smvSstdx2dtRDaxSXUXj7bfm+B97lgSgIxtgCAIeIaCWA2QB6A8hPREv0jjPGknwR7pln5Pte+D8QAf37p+5/+6371/CKJk1St69dcxYd9ZdfnMtjFWX4BDu+BgcOqL1utSLHpiWMQmjnyKH98tEL7/Hxx+568itRZmncu9dZe246rIZDaS7+wQf221I+gx4EogzMzJUx1osxVpcxVpMxNo8xdoIx1kjvuC9CnTol93B0amFghHTqym7oCrcdhcygzNvglr24MtyGVzRrJsyth7Cz7qFlXul3xjq3iY42zuCnzAvCmP7vw4t1OynSFMAAsGqVs/bc9IkJx+OPy/fffFOYrrSKVrwyD0av3FFOijLGzZQp3l3LjZHJK6+kbjvJHGcF5UN46JA77Xo9dx+CSO3YZnWKT/kjP3rU3XwHQWGUA/yRR+T7emtzX37pzahbivKzduqVLDXn7dfPWVvh0AohbidnvFa0BA/C63AFIUXp3Rhugc4Jyt7XuHHWzr95Ux5aWhmaPNLxa8SghXIR0koYDq0591De5/SA0UiISAhB8eCDwjSSFtIpSL949VX3QlcMH+5OO0YoU/kqY7TZQRmGxCW4gpBy5oy/15N6yFoNSa00MVVGLPUSZRIeq1Nd5897H6DPCGUvbs8e8wt8SqsdK2Es0gLFi6tTfEqZOtU4PWy5cu7LpIV0JHvunDBVk1aYPNn9NpXpdV2CK4gQI0bI972IFqlE+UOzYoUgDX3hRmA4KyiT8FgNaa2MbOnnImEIpYIyo2C1pmCWLHFFnIhi+nR1aIhIQzkCHT/eXjt25v+d4iQIY1KS+1kiDeAKIsSFC/J96fy+V+TKJc/6puyZGyF12XcrMJwTrPzQlFN5QYQm//NP+f4ffxibi544oY6+O25c2rde0mPGDOvn+DkCr1dPfW07VjxSPyKldZRXaC0mG1mRSenRQ/05L17sXCYduIIAtJOa2E30YhXpFIVZD84gej1KpP4XAPDWW+bP9cob1gpRUer5cq10oSG0Xh5uReSMRKz2cjds8CbPhhXsZNWTTqd54RSrhzLybe3a4X/XZ84IcdeUeOjYxxUE4F14BDM0bizfNxPuWnlOECj9RdJQGsUUlGlf33xTnansyBHBrFU5whw0KG0lZbKDlYVfIwsor1D6Dzk183Qrh7sZtBwtw0V6rl3bG1kM4ApCCz9fwMoEOWZCRUuVSLdu7spjFq0f4/bt4c9T2tMrFY2fFCqkLouPB+6/XzB3XLBAWLTV+jF7kfIy0iAyN210551A1qzh67mNMoT7xInWzr91yz1ZrKK1mG9kiXjkCLB7t7pcKweHi3AFoYXfC3TSJEKLFhnXVaYq/PPj2nIAACAASURBVOIL9+UxS4cO8n0z4YaVPggTJrgnj1usXQt88ok6OGGItDhasktcnHHmwzvuCNZkWYrVZFPKtMJ+UqyYtfrFi2uXSyMyeABXEFr43TtUusxPn65ft0QJ+b7dIG1uoFSkjz9unNnKzhyx19hZDxkwwH05IpmKFYV4TFo99I4dfbWqMcSqgpCOQDwyEzVEy+di9Wr5/o0b+pnzlAv1HsAVhJJOnfxboJYiNVvt0EHbt0AZ+XX2bG9lCoeWpZdRcEOlc6CXjohmyZNHGDGY5YY/QYUjjmzZgC5d1IYUHvdgwyKNAHzxov2wGW44q1lF6/dTp07qNGxiojB1N3Kk9vl6kXRdhCsIJR6l7guL8oemVFKzZqlDDLRu7a1M4XDqe/Hrr+7I4RRlAD8jghyxRQI//yxMKV29Kvjt5MwZrDzK0Bhly5o7Tzmf73aKUTPo/X5athQ+W2XoHyVepeeVwBWEEi+jUFq9bigv8K+/Ao895q88TlCa6x05ojYTjTT/AaNAdSEiwTw3EqhfH4iNDVoKAeXC7okT5hwvI93IICoqIrz0uYKIJJQWPvXqCZYkykBpgPFcv59o9XLGjZNPxRQvru6xdeniqViWqVRJsNgZOFBuNBDi1Cl1shpO8FSrpk40tXVr+POkC9S9e7srkxXsZq/zIPeDFlxBSLHiyewFVgLuRcrCoFaimRdfFOZOL1/WP+/9972TyS5xccCQIUKmrh9/FKb5atUSnotI+bw5apQ5Fn791dirWunT4iSzoFN+/jm4a5uAKwgpacXxqUePoCWQc/y4dvkjjwg+BUq2btUOVxwpEAFPPCEouDVr1LGnOJGF1ghAz/InOVk9EvQrC2AaJMImggOid+/ImUK4dct4fr5MmWB9H7TQcjgD1IldQihNdTkcJ2gt9kpjlUkJIothOBITrXVOlXHEPISPIADBHrl9e+3ert9ERwPz5+sf37MnMns8Q4ear6ucM+Zw3Gb/fu1y5bRnUEYpUqwabPgY6YErCADIksW/SI5maN5cOyvXv//6L4tZzGbi+vRTb+XgZEy0AicOGiRPzblvn+DYJ2XSJE/FMs3AgUFLoAlXEJHKV1+lbo8cKUw9lSoVnDzhyJxZO0+ulKpVgb59/ZGHk7HQ8vIePDh1VuDmTaB0aXWdtm29lcssfqUMtghXEJHKc88Jpmz//CM8PF7n+XWDfPmMj2/e7I8cnIwHEbBwobp8717hWJYs/stkhcKFzYWi8dnSkiuISKdy5chcc9BDL3BbJE3hcdInuXJZq6+3ThEUuXMDP/2kf7xVK98t6riC4LjLPfdol//1l79ycDIeVkKmAOrYYJHA44/rH/Mh9pISriA47pIzp9zLc/RoYe0kiFg3nIyF1ZF2pFrTaY0iChUKJO4VVxAcb7h2Tcjx3KOHYH2VFtZQOGkfo2jCUg4c8FQMRzz+uOAbUaOG4PDHmL4zqsdwRzmON2TNGkyWMU7G5uBB4IMPwtdTZnKMNDJlCh/N1Qf4CILD4aQfsmYFHnjAuM7Fi/7Ikg7gCoLD4aQvjOJ8jRwZfA6LNARXEBwOJ/2hF3MpQh3SIhWuIDgcTvqjZUtgxozU/ZIlge3b05ZPUQTAFQSHw0mftG0LPPigsJ2QANx5Z7DypEG4FROHw0m/zJsnmIxyizpbcAXB4XDSL5kyRV7+8zQEn2LicDgcjiZcQXA4HA5HE9NjLyIqAqCwtIwxxuM3czgcTjrFlIIgoqEAXgZwCkAoEhsDoJGBg8PhcDjpAbMjiKcA3MEYO+WlMBwOh8OJHMyuQVzmyoHD4XAyFmYVxDgiepeIihJR3tCf3YsS0SgiWktEG4joYZ06G4loifg30u61OBwOh2MPs1NMn4n/35GUMQCWg/wTUTMApRhjtYioKIDVRDSfMZYsqRML4CxjrKmJ9roD6A4AXbp0Qf369a2KBAA4ffo0EhISbJ0biTDGkMySER2l/xWlt3s2A7/n9E9Gu1/Au3s2pSAYY26aw1YD8IfY7lEiOgngNgAHJHUqA8hLRIshKKKBjLG1OrKNBTAWAMaOHcsqVKhgS6iEhATYPTcSGbZqGAYsGoAT/U6gYHbtJCrp7Z7NwO85/ZPR7hfw7p4NFQQRdWaMfUtEfbWOM8Y+0yoPQ2YAZyX7VwBkUdS5CeATxtg0IqoGYAYRlWZMmsuSY8SARQMAACPXjcSQxkMClobD4aRFwo0gyon/K2scC/uyJqJuAJ5WFC8CIA2MEgvghqLOLgDbAIAxtomIzgAoBuC/cNfkABSfGrFy6Iqh2H5qO2a1mxWgRBwOJy1iqCAYYwPF/12Vx4ioUbjGGWPjAIxTnNccQA8An4vOdwUBHFac2htADIChRFQIQByAE+Gux5ErhxCzE2YHIAmHw0nrmFpbIKL7iWgaEc0V/34FMM3OBRljCwAcIqKVAGYD6M0YSyKiQkS0RKw2GkA1IloO4CcA3RhjiXauxxH4aOVHQYvA4UQUSclJoHgCxROSkpOCFiciMbv4PB7Abgie0/MgTAl9YveijLFejLG6jLGajLF5YtkJxlgjcfsKY6wtY6y++LfY7rUyEscuHdM9Nnj5YB8l4XAimxu3biDTe6kTKO1mtANf4lRjVkHEMMYGAZgOYB0Ez+pWnknFsUXRT4vqHruaeBVXbl7xURoOJ3IZtHSQbH/mzplY+O/CgKSJXMwqiCgiygIgAcBdjLHrAPJ7JxbHC45cOhK0CBxORPDRKvWU64HzB/wXJMIxqyAmAlgq/r1LRJMAXPdGJI4dEpPCL9HweVYOB0hO9cmV8eK8F32WJPIxpSAYY0MBdBfjMT0L4CCA9l4KxrHG+yveD1vnpx0/+SCJMbeSb2Hx/sW4fPNy0KJwMihDluv7Bd1KvuWjJJGPWSumfxhjWwGAMbaMMTaIMbbLW9E4Vth2alvYOsp5V7+5lngNMe/FoMnkJnj0h0cDlYWTcRm3aZzusefmPuejJJGP2SmmtUTUxlNJOI6YsWOGqXo7Tu3wWBJtbibdRLb3s6XsL9q3CNO3TQ9EFk7GJSk5Cf9d1Pe3/fbvb32UJvIxqyByAZgmRlj9hYjmENEcLwXjeMOkLZMCuW7hjwuryjrM7IC/j/8dgDScjEpisnyt7p367+jU5ADmo7n+Jv5x0jiHLhzy/ZqMMZy7fk7z2KZjm1ClcBWfJeJkVLaf3C7bj28UjwdLP4i6E+sGJFFkY3YEUYQxNln6B6CIl4JxzLNg7wJV2fCmw8EGqR1/ftj+gx8iyThz7Yzusa5zVVFcOBzP6LtQHXe0ciGtUHMcIIyCIKLyRPQogF5E9Kjk73EAffwRkROOFlNbqMq63iO8eJd3We63OCoemfaI4fEbt5SxGjkcb1hxaEXKdui3kStLLlmdtf9pZhaIWH7c/iNWHlzpSdvhRhA1IQTOywugl7jdG8ALAMZ4IhHHMS/XfBlxsXEAgHol6wUsTfgfXNahWQ2PczhusOnYJtm+3m/j/m/u90Mc12g3ox3+PPCnJ20bKghxOqkxgN8BDBa3ZwLYB4Av90cA56+fV5WNaD5Ctr+081LZ/rlr2usBXsCd8ziRwu4zu3WPPVDqAR8lSTuYXYO4COBpIroLQF8IeRnGeyYVxzTKXpEWNYvXlO0/9P1DXomj4uSVk6bqzd0112NJOBkdaWdl8dPy+J9z2nOjTC3MKog6ALoBaAZgNGNsMIDinknFsU3m6MyqsqyZ5FM4W45v8UscjN8k70ewQQyLOi1S1Ws9vbVfInEyKEksVUEofxPZYrLJ9vec2eOLTE7x2vPbdK5pxlgygIYAQqueMV4IxLGGMq5MizLqBWsAKJS9UMr29VvXfQtt/M5StZ357Xlu16zLwxxwvOSL9V+kbN9X7D7V8Walm6Vsp5VQMPP3zE/Z9mI616yC2EBEKwBUALCJiCYC+Mt1aTiW2Xt2r2x/9EOjNetVK1JNtr/nrP89pBeqvwAAKB1XWtWDA/g0E8db/jqa+sqKjopWHc+fLTVA9efrPvdFJqecuJKaaFPpBOgGZhVENwipQ1swoeu5DgAPWhIBKCNQFstVTLPe583lD7wfvfWjl47K9muXqJ2yfaiP2mFv6YGlXovE4egi7VyllZAb0hA7Wp0up5iN5npZtGjaK+5/xRhL19lnRq0bhQmbJwQthmuEzF5DENS5q93m0zWfyvY73t0xZbtA9gJoW7Gt7Pio9aM8l4mTMZGGw+9fu79mndxZc/sljmt4neTI9BpERuOlBS/h2bnPRnQawombJ8r2e9/XW7duvmz5ZPuHLx72RCYpZeLKyPaJ5EpJui7CSV+0md4Gg5YEGz1YyssLXk7ZHtZ0mKlz9PJGZCS4glDw79l/QfGpL7Lvt34foDTGKPM7dK7S2fS5zb5rFr6SQ6TTX0s6L1Edf7fhu6qySFbInPAkJiWixdQWmLNrDgYvH4z95/YHLRIAYFbCLFP1Xq6ZqkjW/bfOK3HSDFxBSLiVfAtlRsl7vR1ndcS/Z/8NSCJjlAtt1YtWN6zf494eXopjiNKMEBCmmZSYSXzEiVye+OkJWWyw45ePByhNKmbleLjswynb129FdtLMzcc2e34NriAkxLynbbmrVBqRQqYos8F4BT5p9olHkoSnRtEamuV5suaR7U/dOtUPcVxj79m9WHN4jSpKaEZlzi65w1ntCbWx6tCqgKQRkI5K21Qwn9bGC6sgN3lz8Zsp2y1Ka5u3O4UrCJNEooXN7ITZlup7YeVgFuX6Q4jZ7eT3sPP0Thw8f9APkVyh7KiyqD2hNu4acxconnx1Qow0pFOzUrSiDfvJzaSbKds5MucwrNvw9oYp273n66/pRQKXblxK2a5RTLsD5hSuIETil8UbHm/0bSOfJLFHkzuaWD7HKDaNU8Zt1E/rKEVLcZT/orzb4rjOsUvHNF+I93x9D85eOxuARJHLkBVDAnWCvJKYanB5d8G7DevGRKfOInj5+3DKjVs3sOpw6shMrwPmFK4gABw4f8BUvcX7F4evFBDz/jfP8jn/XdBPveiUMRtSg/0+WkE//3T5fGplcCMpssN/t5rWCkU/Lap7PN+wfLrH0itd5xjn9ej1Wy+fJFFz5WaqguhTK31kKdBLwOU2XEEAKDeqnKl6S/arLXGCQulWnyVTFlPnyWzAPXSFkPZoZj45U7deoRzapq6Rmor0l12/4Jfdv4StV+STImkmXIMbTNwy0fD41xu/9kkSNd1+6ZayLR0hmGHD0Q1ui+MKfo3IuIKA+cWoISuGeCyJeez2snNmzpmybZTpzSnSKLPhhr8V8ldQlUlDCEQKV25eQavprUzVPX75eESuW3nBkYtHTNU7cTmY79SqM1m+2NQR4Ji/IjPtzX8XU0f/b9d/27PrcAVhkUjJfmZ3NCN1mNtxaodb4siw6suglfXuveXvuSWOazz242OW6mcUn47in5kL7Hw18arHkrjDhu6po4brSZFp6ipNajS40WDPrsMVhAa7eu3CD49r526OlOxnLae1TNmWBhkLx/PVn/dCHBlSq5E78twRtr6WP8TKQ96kULTD9VvXQfGE3//93dJ5raa3inhbeqdoLcg/WelJzbpR5P/rZtfpXSnb8Q2NDVFCFMlRJGX7WuI112VKS3AFIWFRp0W48uYVlMtXDk9WehJskHYPcP2R9T5LZkyVQlVM19WKYuk2UquR7Jmz224nUnrgI9eNtH1ukU+KhK8UIIwxrDm8BhRPeGbOM5YVmjLfx8KOC/HD4z9gxTMrMLDeQNmxICKkbjy2MWW7ZrGaBjVTka5TSDs7GRGuIAAs67IM3at1R5NSTVQev1qONTXHm3vQ/ILB/ovUi4VUqdXItpPbTJ3zQZMPVGWfrAnOsU+Kk1GAVkrYSOKbzd+g9gQhyu6kLZPQcFJDS+crFcSDpR8EANS9rS6GNJav2X229jP7gtokmlI7RGZ9BaQjnXl7rFsHeo00B0SpvKU8vRZXEADql6yPIjm1e3pfPvSlZrk0OqTfuLkO0n5Ge9faCvHbnt9Stv982lwyda0Imz9s157m85ufd/7s6Pz3lkXeekoIqYUPAKw7ss60BRljzHJeEb9jM0lHzMqIxkY0vqOxF+K4ws7TO1O2H6tgbV3MKlxBhKFIziLY23uvqvzecfcGII2Acn7eqhXD3PapiXnWHXE/INkL815I2Tb7Q9Oa+ooEE8Njl47h7xPOTG7fWfpORJlIhzh99bRmedWvq5oa+bw0/yXZ/ownZqjqzGonD5LnVNla5Zk5z9g6TxrGJlKmOkNIrcH619EOXe4WXEGYoHRcaVXZPyf+CezBkfYgAHl4ADNULVzVRWm8heIJZ656Z44bjk6zOrnSTuPJjSPOL6LAcLVxQIi8H+U1PDeZJeOLv1JTeMZExaDtnW1V9ZRTtFtO+BuKxO5nLs2X4pW1n12GrU4NV14we0FPr8UVhEmmt52uKhu3yVw4CbdxGiNG6peg14uMJNb+tzaQ6x6+cBh/7jc3RWYGZfa/IDFjcmo0imgwqYFs/9s25jKwfffPd6bqBY10HeKfE/8EKEmwcAVhknZ3tVOVRUIwr4+bfmz5HOWiq5shmS9cv2D73F87/KpZ7mQR3gld5xqHj7BKJL0c60+sH7ZO3o/yygLCSVFOc3ao3EG3ncI5ClsTziWcWCBJF9j/9/P/3BAnTcIVhANuJt0MPOuUXg5qI3JnkadWdDP9qNQEWJlSNBzKrHchgggXnZiUiEX7Fvl+XT84deWUzPzTiFwf5lLFIGsy2VpgSKW5q1/k/CA1asD656yZplcrUs1tcVxBGmJn3CPez2BwBWGBbS+qTTajB0er4iL5idWcEIC2Y5pbSKcu9Bym9KhVvJZm+YerPnQkkx28ykux89TO8JU85vVFr1uq/97y91Ji/+w5s0elMLpV66Z1Wgr/qyzvge85Y83yyS7SEYRX4bD9ptY3qb8Rvd+Lm3AFYYFKBStplk/+e7JvMigXxrWioVrFzcBfXoVTOHXllCft6vHCry+ErfP4nY+j490dceH1C9jZ09yL/84v78TGo+Z6714xYcsES/WXHliKmPdiMPnvyaj0pfo38H4T4yyASvPS5lObW7q+XZqX8ec6fiK17HPjtx+OwBQEEXUkomkGx0cR0Voi2kBED+vV85umpZqqyrrO7YrX/njNl+uP/mu0bL9yocqO23TTmUsahtiOlVffWn01ywt+7K21hpSk5KSwwRCL5yqOn574CVMenYJcWXKhQv4KuPj6RfSqET6s9b3j7sW5a/6Ea3aTzrM7awa2NBPqpV2l1DW8fef2+RKNNJSoqESuEo7bWnN4jeM23MZqZFo7WJ+fcAgRxQH4A0BRAEt16jQDUIoxVouIigJYTUTzGVNP+BNRdwDdAaBLly6oXz/84psWp0+fRkJCQth6n9/7OQZfVQfHOrzvsKnznXL92HW0L5zq3Gb3ms+XfB4lYkqgfeH2GDJrCAY1GOSKfCs2rUiRL+ZcjGX5upfsjmMHjmke27x1M2JjYh3JZ+Z7/v3f32WfsZKeNXoiLjZOs52epXri0QKPhrVw++b3b9CyfEvDOoCQq6RAtgKOQpZI7/nv438b3ptVKhWoZOo7HlB+AOhM6lrXth3bPMtwePr0aezcuTPlPs3KqOTNim9i37l9AIB1f69D3ivGpr9ec/2W/m/f7PvLKr4rCMbYWQDViagBAL1xfDUISgSMsaNEdBLAbQAOaLQ3FsBYABg7diyrUEEdOtoMCQkJMHtuvcv10PO3nqrySicr4bU6r+H01dMYvmo4KheqjDvy3IFGd7iXjW7YrmGYfjzV5HZaBd1BmCF98vdB/Mz4lLbstqNk+g+psk18dqKtl4C0DSk7Vu3A3y84c1oz8z1X/KGi7rFeNXqhbrW6hudXQAU0+c14IXf68emodGcltCibmks4MSkR/f/oj6erPI2hK4bKnMpalGmB3576TaupsEjv2eje7PBBpQ9M/W7OXz+P6T+nfq/fVvwWmaMzuypLiISEBOQomiPl2e6Qv4Pp37aUctfL4f05wvTZ9OPT0adlsMmGEk4n6P72rby/rBCpaxCZAUjDKF4BYC4jjg/0qNFDs/ztJW8jy5AsKPZpMYxYNwLPzn0WjSc3dtUaRjo0f7PumwY1jdHKweA2dnuIdW/TfgH/c+IfHLukPbpwC6O4S6MfGo1RD40y1Y4ZRfbQ9w/h4o2LuHTjEnrO64nMQzLj83Wfo/rY6iqP4/l75ztOZeqFo97rdc0teOfJmke27/X3WHlM6tSrmYjCWgRlXq1HxdGpyt1pR8ksnioIIupGRCsUf8YmDwKJAKRvl1gAkZGIQaRsXFnTdZtOaYoftrkTV2jKP1NStvvV7udKm5FGuTj9DH9GqT7doPX01qqyUS1GgQ1iuh0DLe4uZJz7OETuD3Mj14e58OUG7ZhfUvINy2c6OY8WboceH9FshKX69xdPzWHQZ6F3vXHGmGxd7d2G79pqJ2gTdiPMPl9O8VRBMMbGMcbqKf7MGO9uAtAEAIioCICCAA57KatVrMY/aj/T/aB4Xg3R7eKWpVG4Xrqe85YbaOV86HWfvXzKI5vbDxOuxwNTHrB9bvPv1FY9maIyYWJr43ShenSu2tlS/dolaqdsz06YbeuaZgitG4Swu5gbyQrCLyJmiomIChHREgBgjC0AcIiIVgKYDaA3Yyw4ZwMNOt7d0fI5VhPOKFFaBZnNQ22GebudhzV+d+m7KdstyrTQrxiGbDHZ8Or9r+oef+PPN2y3bcTes+qgjE7oXbO3amrFKQmn7S1E3ky6qXKOy5s1LxLfTkSXql1wfoB1Szar96YV0t0LpC/2GkXt+z8oz42UbJIxUd5bL4UITEEwxpYxxjpI9k8wxhpJ9nsxxuoyxmoyxiIuKDsRWQ561+y7Zo4CzylfDtJY93aQxpJXBgC0w4UbqWE2vnvMWViJOiXq6B4b/ddoLD+oTlPqlDbT1bk/vmjxhUZN8/zzgvtxfOzc++rDq1Vle3qnOqzlzpobN98yH5rCzovXD7NMANh+anvKdvfq3W23U6WwPBHXhyv9d9gMMXbj2JTtzlWsjdycEDEjiLTInPZzLJ+Tf3h+21Fg7/zyTtm+NOieHcrkLZOy7UZkWmnPzY6Ht5Tbct9meFwZLM4NpC8WAHjizifQ8z61tZoVSuQugSzR7tpXNJjUwLL3vvQFE0LpwBYTHYPDfc3N5C7suNDS9bVwOqLWQxqePXuMffNgJWeuBRdV+PlfU1MFf/iAf4qKKwgH2A21O2eXdcWixK5lhh5OlQ0ATNuWanbn9KVYvWj1sHV+3a0d3M8OysBurcu3djwKClE8V3FX2pGS6b1MWPef+Vwe0u8GAHb02KH5nRfPVRybum8ybGtAnQHIG2vPJ0CaH2TLce9Df5fLp2/wYJUxG8bYPnf9kfXo9VsvW8EsNx2Tfx96Mcu8gCsIB2TNlBWT21gPs/HoD486vvYvHX5x3IbUA1b5AnGKG+sjXap2MTz+yLRHcOjCIcfXAYAsQ+TytijTwjUjgCWdvUkW1PZHc8EQt51QxxCrWEDfH+KeIvcYtmfWtFWLVuVapWxfvHHRdjt6KDM9muloGPFKrVdStu14fzPGUGNcDdQcXxOj/xqNGuOsT80dvXTU8jluwRWEQzpVsZdQxurDduD8Adm+XlwoK5TNl2qqq+ylOKFYTusRZrUY21I9LaLk5QUvO77ON5u+UZU9V+05x+2GKJG7hCsKXcmRS+ZMXmcmzJTth5u+A7Szw4VwsvAuDdUxdMVQ2+3oMXz1cFfbe6fBO47On7hloix+0p6zezD578m4dOOSaYMDaXDEaW3d7ciFgysIF/ij0x+Wz9GaEzbi8IWIsvJVIZ0TL5SjkCttxkTHYFEnYydDN8wlv9ksVxCFcxTWTIHqhJblWlpel/ms2Wdh69hJZbqxe/hggW3vbIsLr1+QTWVObjPZlrWTlOpFnPXowzFz58zwlSygfA6smL6OXj8az859VlXeeXZn5PowFyqOrhjWZPvopaP4bG3qc1CpgPOOoRW4gnCBB0pZt03XCtVhRP1J9mJM+YXUy9fpArWUJqXC5x4YtMRZHKk1/8kDsVnNHWCW0/2tZe/rU6sP2CCmStsppfFk45zfc3fNVZWZCa4HALmy5EpZjP665dfoVKUTcmfNHeYsY5RhZ7af3K5T0x7SkfCZ15wvKiuf5Z7zzP1uZ+2chV7zw/vPzNtjbKBZ7FP5aLxEbueBB63AFYRL2Jmv9io0tl2U01hWkIZwdmp+q+T1OsZz3oOXD7blJZyYlIj8w+Qvy3n/m+fZj9DKyzXx7dSpmFntZhnWvZZ4TbOcMabyDF/xzArTMgDCNOSF1y+Ezflgl7vG3OVJu4Dg5+EUpYL4auNXps577MfHTNXrMLMDFu5dqGlFqJWB0G2/mnBwBeESdobe/5tpL5Vhvdvq2TovHE/+ZC3BjxRpz+3xOx93Q5wUXqsTPpR67NBY9F2gHSpcj5NXTqpMFx8q+5ClNqxixjz06CtHVS8mqc+KkmzvZ9N8wShHRgBwe57bwwupIFeWXK5YuYVQ+g+55YA28E955jo3ZLYzGtZT2Ho0n9ocH6+Wpw4+cP4AOs2yt77pJlxBuISdMNRzds0JPC+ANNGNW2aHejkd7JI3Ni9eqx1eSYxYNwJfrDfv2NbvD/9jWT1Y+kFMbjMZgxvKQ8ZPbD0RiW8nIvmdZBTJWUR1XrhF7vdXyJP2nL12VmUGnDk6sycmt1b58fEfZfunr1qbetPi0o1LeH+lceIiu0hjSAFCx8KIrzaYG2VIeW3Ra6B4wqUbl8AYwx2fq83Y32v0nuV2ncIVhIus7qr2Vg1Hla+qhK2jNN0Ll8HLCrmy5ErZLh1X2pU23exthvio6Ue48uaVsL4nvef3DjvdxBjDhM0TMH2bPKz4sVe9jTAaolOVTni7wdtgtA0bTwAAGwRJREFUgxhO9DuBm2/dRJeqXZApKpPuZ3dngTuxo8cO3TbfWvKWbD/fsHz4YKU8tMWJfiecC+8CSYqoOW4kD9p1Zpds/7l73LNCU34nxy8fN6yvNXIzS64PcyFqsPZr+anKT9lu1y5cQbjI/SXux9U3ra0rHL4Y3jpJOfzUC4dtB6nysRvnxy+yxWTDiX4nsPKZlYb1YofGGk5bPDDlAU3rksI5CjuW0SoFsxc0HYLCyHcBACieQPGE77d+r3nc7/lrPZRm0NI8y3ZRTrGNa2UmJqg99AxMriZeBcUTftrxkyfXvSOvu86xZuAKwmViY2Jx9jVrcfvDmbodu+xdz1a5cKocrZjBjTAdVqhzWx2wQcbXzDo0KyieVN7Ge87skdmVhyid153RUyTw1M/+9zStkDNLTpmnfbgeeTjWHF6D+8bf51QsXZRGEisPaXdQ7HhJmyWo0R9XEB6QNzYvHiz9oOn6sxKMrVTczBmtJE/WPLLcFsqQE2ZYd8R8yAc3OT/gPB4p94hhnVrf1EqZctp6Yiu+36bdu179rPXpwSCQrhlZ4clK9g0QvGB8q/GyfSe5Kh6ZZvwMOOWR8ur2tSwQT1zRf4mzQQxsEENCT3ujdLOmyW7DFYRHzG0/VxYt04jdZ3YbHpcmCXqm6jOO5NIiZ5acKdtGD7keUic5PxfScmfNjbkd1Hb+SmKHxoLiCXd/pZ9kxW5cLb+pkL+CrZdFEPPXRijD5Q9bNczS+XvO7EmZUlNaoj1f/Xmds9wjfmm8bH/VoVW452vtECXdq6VGlC2fvzwuvG59pBFFwbyquYLwiCyZsqBMXBks7xI+NLOVkANevIClU0R3j7GeqUq66OiG7blVTvYztioJx9DG7od88JJ1z1kfsbUq3yp8pQAZtDS8s+P2k9tB8YSCwwui6tf6ofb9WEsatnoYKJ6QmJSIa4nXMGDRAN26I1vIE0flypILL977oulreeW4aQauIDymXsl6ONjnYNh6EzdrZ/VSLlBLe/tuMaBO6sN9JfGK5fPXH0l9gL2wYApHgewFHIVwaFepnYvSeE+pvKVw/FXz8/aDGjjzNPeLPWf0R9yL9i1Kcao7dfVUxDiZ3jXmLmR7PxtWHV6lW0crcOWXD3+JQtnNhaSpUcx+0iOncAXhA2aCo3Wd21Wz/J0l8mBhbucWAIB2dzl7Qfb/o3/KdlBD4Q3dN4SvpIObmfn8wq14V0ES31A+TVPuC3Vo7jNXz4DiCU2nNDXV5vuNvfGF0FvDCTc9XKu4voXWqq76SiVEh7s6hK3jJVxB+ESVQuH9HZSBwJJZMq7dkntlRvrLLMjF0MtvXLZ8ztz2cyPCecwrHq3gPLS8V/yvsjqSwNYTW1Piel2/dR0Fhhew1Kbd/OHheOAOe7nA1zyr7xNROq502PhcA+sNNDzuNVxB+ISZcM9KCyKjIbeXOMnNrMxS5ifZM2fH9YHWrGG0LFTSCuGMIGKiYvBzu599ksY6ZeLKqMru/upu5BuWD3N3zUXs0FgwWDOh9mIKFvBuxJYvWz5seV47gsHe3ntdCevvBK4gfKJE7hJ46b6XDOs8M0duoTR+k9wU8LGK5gKAOeXV3181XXfmDnfDKzslS6Ys+KDJB+ErpgPKxJVBzWI1dY+bCRceqSiDDJrBaDrHKeHMqbUwaxhQNGdRzXI7cbPchisIH9F7EEIoQz98vEa+QD31samuyxRCOs1iZeQyf+98L8RxxOt1X0f/2v3D1jObfzmS0TPzLZKjCHrU6OGzNNbR6z3b4c+n/3StLSV2jC++e9RcytoC2Qvg7GtnkfROEs68dibFZ8LtnCR24ArCR1qWaxm2zudrP8fNpJv4bI2695c1U1YvxAIgj12z87R5ZyypiWvlgpVdlckJw5qGt6tPD2sPBbMXxK231bGM/ur2VyAWZVbJFpPNcRtxsXG4/MZlV9pyEyvTXXlj8yKKogKdotWCKwgfqVSwEm68ZRzauM/CPij6SVG88vsrsvJ7i97rpWiqAGpmkY56fu/0u1vieMqiTotcSSYTKURHRaNXDfnibLFc7qR99Ro3lPTyLsuRPXN2F6Qx5ucnza/npIfOB8AVhO9kjs6M7DHGD7PSMxSQeyt7AcFeb1MaIsELE1wnjH9kPKa1nYZ9L+1LKSufrzyalGoScT01p4x6aBTYIIYJrSZgQzf7Jr9+ExsTa5jrwgx+LeQ2K9PMdF07CcQiEfdyQ3JMc6jvIeQbls/SOV6bj75y/ysYvHxw+IoGuJlq1A2erZYasfX4q8ex5r81KM/KByiR9zxzj/uhWLxmYceFKDuqbPiKIp3u7oSu93TFzzt/9jzHtZSYKHNRdwFg/lORtzZnBz6CCAA7vVczWdWcoIzquuZw+Jj2ysxZkbCopkehHIXQpkKbNDEvn9EoE1cGE1pNMFW3Z42emPzoZDS8vSFGthiJzlU7eyxdKmY7QB82+RDl8qmd/tIikdXl4+jit4fypC2TcH+J+w3rPPS9PD1npI0gOGmHZ+55RjeaQIgN3TagelF7I4Zp06bh0iXjsPpmeCfHO2ETHEVvjsbYHWMdX8sKiYmJWL48fNy3nDlzokMH897Z/BcdEJu6b0K1sdWCFkOXsZvG4utHvjass/TAUtk+VxAcJ/Sv3R/DVw/XPPZXt79sKwcAuHTpErp37x6+Yhi2ntiKG0nGhial85ZG3lh/g1Zeu3YNsbHh0x6PHWtNcfEppoC4p8g9KJDNXBgBvxa8prWd5uj8oOIwcdIHT9z5hGZ516pdPbfiM0u4ZzxzdGbflYOX8C5fgPz9wt8o+qmx8xwAfPWw9STodmhaylxANEBtVRVJPhCctEmNYjUwq90sNCvdDPvO7cPNpJu4p4h2joWgKB1XGttObtM9flfBu3yUxnu4ggiQfNnMWTL5teBlZZF547GNsn0+euC4QZsKbQD4YLpq1lhBkU43a6asKJqzKI5eOqqqGkVR6e53kL7uJo2ROTozRjYfaVgnNlMs6txWxxd5YjPJ5zB/2PaDbt03/3xTtv/lw196IhOHE2lIrRA3rtmIgT2EiKsV8lUAALzxxhv49ttvU+osW7YMdeqY/w3PmTMH//77r+n6ffv2xX///We6vhW4ggiY3jV7Gx4/0OeAP4JACHQntfVuP7O9bt0/98vj3tQuUdszuTicSCJrpqwpyiBEFEUhW2Z1qI969eqhY8eOSE5OVh3TY/bs2di713xE5c8++wzFi3vjuc2nmCKApqWa4o99f2gey5Ull6+y5I3Ni5NXnKXw5HDSO7ExqaPtPFnzoFoRbYvEFStW4ODBg2jfXruz1b17d+zatQvXrl1Djx49QERYsGABNm3ahBw5cuCtt95Co0aNsHLlSvz000945plncPLkSZw/fx59+vTBc889h0aNGmHEiBGYPXs2jh07hiNHjmD//v0YOnQoWre2HhVXCh9BRAALOi7Ag6UfVJXXva2upwH6tBjRbIRs//RVdUKTE5dP+CUOhxORREdFo2xcWZTJWwZLlyxFo0aNUv6+//57U21cvnwZS5YswezZs7Fs2TJcvXoVnTt3RvPmzTFs2DDZtNTvv/+OXbt2oV27dli5ciW++uorTJyoTlN88+ZN/PLLL5gyZQpGjjSevjYDVxARQBRFYXa72aryININdqgsv+aQ5UNUdT5Z84lsP+isVxxOEOTOmhuZojOhcePGWLJkScrf//6nzpSnRY4cOTBq1Ci8/PLLeOKJJ5CYmKhZr1UrIa9EsWLFsHr1avTv3x9Tp05FUpI6PluNGkL+6ri4OFy96jxvN1cQEUJsTCwuvH4hZf/7x77Hi/e+GKBEAp+v+1xVpnRm8nsajMNxDGPm/jzk2LFj2LJlCyZPnoy5c+fi888/x4ULF0TxUq8dcoAbPnw4qlatiuHDh+PRR7VTyUZFuftK52sQEYT0RavsyQfJhqMbDB2VWpd3Ns/J4WREChUqhH379qF27dqIjo5G69atkTt3blSqVAmvvPIKcuaU55N4+OGH0bdvX8yYMQM1atTAmTNnsH79ek9l5Aoiwvilwy8okatEoDJkj8mOK4lXUvZHrR+Fb9sIZnvbT25X1S+fP31HSOVw9GjQoAEaNGggK/vgA3nK25IlS2LNGnXwy6ioKM3QF/369UO/fv0AAEuWLEkpb9q0KbZtS3XSe/fdd1PqXLt2DYMGDQp7TasENsVERB2JSDe2AxFtJKIl4p/z1ZY0QstyLVGlcJVAZXjlfnmyosl/T8b56+dxK/kW7hqTvjxFORyOPr6PIIgoDsAfAIoCWKpTJxbAWcZY2NgPRNQdQHcA6NKlC+rXr29LrtOnTyMhIcHWuWkVvXtuW6At9hSW56V+8VthPaR9YbW5Xlr63Pj3nP7Rut/ExERcu3ZN54y0z61bt0zdX2JioqVnwXcFwRg7C6A6ETUA8IJOtcoA8hLRYgAMwEDG2Fqd9sYCGAsAY8eOZRUqVNCqFpaEhATYPTetYnTP02dO1yzXYloFZ0H+/IR/z+kfrftdvny5qWinaRWz0VxjYmIsPQuRugZxE8AnjLFpRFQNwAwiKs2Yx2YFHA6Hw0nB0zUIIupGRCsUf91MnLoLwE8AwBjbBOAMgLSRhT2dYNbEdvHTiz2WhMPhBIWnIwjG2DgA42yc2htADIChRFQIQBwA7r7rI92rd8eYDWPC1qt7W10fpOFw3IXizUVzZYMy9qRFxDjKEVEhIgrZdI0GUI2IlkMYSXRjjGm7GXI8oUiOIqbqxESbT+TO4aQ3li1bppvC89y5c6hRowZ27dplqq2DBw9i5syZpq+9cOFCjBtnp/9tnsDWIBhjywAsk+yfANBI3L4CoG1AonEAFMpRKGyd79uaiznD4WQ0Pv74Y4wdOxaHDx82fc6BAwcwY8YMtG1r7tXXrFkzu+KZJlIXqTkRwJz2c9B6ur6XdMPbG/onDIeThgg5uzVq1Ejz+Jw5czBs2DBERUWhSJEiGDVqFPr06YOjR4/i7bffRpkyZbBixQocO3YMPXv2xJ49ezB16lQAggf2nDlzMGXKFCQkJOCFF15A+/btUa5cOezZswdly5aV5aNwQsRMMXEij1blW+key5E5h4+ScDjpi2nTpqFfv35YsWIF2rRpg+TkZIwYMQKNGzfGe++9BwBYt24dfvzxR7Ro0QJHjx7F2rVrsXbtWuzevVuVUGj37t0YMWIEVq9ejS1btlgauRjBRxAcQ3JlyYWLNy6qyn/t8GsA0nA46YOPP/4Yn332GSZNmoS4uDg0bar2CW7SpAmyZ88OAChYsCD69esHIkJiYqIqkmupUqWQN29eAO5FcgW4guCE4fnqz6uit+aLzYcGtzfQOYPDiXyCtk4aP348Bg4ciLi4OAwaNAjTpk1DlSpVNKO4btmyBXPnzsWyZctw/fp1zJ6tTg3gdhTXEFxBcAwZ2nioSkEc7uvO8JXDSQ8sXrxYttbQq1evsAvNVatWRbNmzZA9e3ZER0dj0qRJSEpKwpo1a/DWW2+hbNmyKXXLlCkDxhgaNWqEYsWKoVatWpgwYQIqVark2T2F4AqCY0hMdAxuvX0LG49txF9H/kLP+3oGLRKHEzE0aNAAJ07ou2hJo7FKadOmDdq0aaMqP3jwoKosR44cWL58uaEcS5cuDXtNO3AFwQlLdFQ07it2H+4rdl/QonA4HB/hVkwcDofD0YQrCA6Hw+FowhUEh8PhcDThaxAcDifdkzNnTs30numFxMRExMSEj4umzHMdDq4gOBxOukcvoF56waukUHyKicPhcDiacAXB4XA4HE24guBwOByOJulqDeLYsWOn4+Pj1a6IJli1alX+OnXqnHZbpkiG33PGIKPdc0a7X8DxPZfUO0DS4FAZGSLawBi7N2g5/ITfc8Ygo91zRrtfwLt75lNMHA6Hw9GEKwgOh8PhaMIVRCrp14tGH37PGYOMds8Z7X4Bj+6Zr0FwOBwORxM+guBwOByOJlxBcDgcDkcTriA4HA6HowlXEBwOh8PRJMMqCCIaRURriWgDET2scTyaiCYS0Wrxr5FWO2mJcPcs1ulKRJuIaAsRDfRbRrcxc89ivUzifet6laYVTDzb2YloOhEtJ6I1RHR/EHK6gYl7NfX9pyWM7sn175YxluH+ADQDME/cLgrgAIAoRZ0WAH4St8sA+CtouX2453IAfoPQcSAAbyvrpKU/M/csqTsAwF4AJYOW24fv+SUAg8XtygDWBy23F/dq5ftPK38m7tnV7zajjiCqAfgDABhjRwGcBHCbos45ALmJKBpAQQBpPbaLmXt+GMB+AHMBLAewlTGW7KeQLmPmnkFEZQGUArDCV+m8wcw97wAwQdy+CCDaN+ncJdy9mvr+0xjh7snV7zZdBeuzQGYAZyX7VwBkUdTZDCAGQAKAwgCe8kc0zzBzz0UgPIBNAOQAsJqIVjLG0qpyDHvPREQAPgDwIoBh/onmGWHvmTG2CACIqASAKQDe8U06dwl3r2ae+bSG4T25/d2mewVBRN0APK0oXgQgq2Q/FsANRZ3XAaxmjA0UP+zVRPQnY+yKd9K6g4N7vgbgN8bYVQBXiWgjgIpIAz1rB/fcDcCvjLFTgq5IOzi4ZxBRZwBdALzCGNvglYwekwjjew13PC0S9p5c/W6DnlMLaB6vOYC54nYRAPsARCvqDAHQR9zOAeAggGxBy+7xPTcBsADC+kMsgO0ACgUtu8f3PA3AUgBLABwDsAZA3qBl9/ie2wCYrCxPa3/h7tXMZ5HW/kzcs6vfbeA3HOAH/QWAlQDWAXhYLCsEYIm4nQ/AfAi9578AdApaZq/vWdwfCmCTeN8Z4p4ldScijS9Sm7lnAHMgTKEuEf8mBC2zW/eq8TyrPou0/md0z25/tzwWE4fD4XA0yahWTBwOh8MJA1cQHA6Hw9GEKwgOh8PhaMIVBIfD4XA04QqCw+FwOJpwBcHhcDgcTbiC4NiGiF4moiLidhsiGuPhtWoTkWbeXSKKJaJFRBRrsc33iegZdyQMe61cRDSBiEZpHCtJRBuJqK2kLJqIPiCiXUT0rxhZOFZy/C0i2iseG0NEUWJ5NSJaRUS7iehvImouOecXItpHRDvFv/fE8tuIaI54rV1E9JKGjNnE9jpLysqIEUN3idF/7zcju87n05CIbln7VDmeE7TTB/9Lu38QAvtV9+E6UQC2AChsUOdZAAOD/kx0ZBsHIajaFQCjFMe2ATgF4CaAtpLybhACJmaFEH9nLlKjdHYDsAFANvHYKgBPQgidcwDAQ2K92gAuA8gh7mtGA4DgSf6quJ0PQsC3Zoo6kyB41ncW90PXaifutwSwI5zsOp9PMQgB6FjQ3xX/k//xEQTHFkS0HsIPeyYRPU1EnYnoF/HYEiJ6h4hWENEhIhpGRIPEfAtHiOhZSTt9xN7nASKaSUR5NS7XHMBBxthxIqpIQn6OXUS0h4j6i3WmA3iRiFTxxYjoBUnv+P/tnW/Ml2MUxz9f/xcWtUppKqTIimqZmRdUyJqENCsT07xIi5GkF2b+bNlEkvmTJRuqaWWKF80sTbZ4gVBqrVYLY6SSovh6cU6PX7/dT8/zFIVdn1f3fV9/7nNdz7PrnOtc9++cjyX1zeezJd2T105ZvkiZH6tpPzYt7jU5pt75vF3G3v88LfPxVXNle6zt9lQEA7R9vu12xCJfS1dguu3dtn8D3gXOyrLbgYdt/5Jl1xG/mu0ArLf9dva9glA8HTOemIBFubuYI6ld9tcfeCPb/EBEEBhWM/47iaCVK2vkGwxstT0v75cAI5oh+35IOpaIPjquqrxwZCkKonBQ2B4AbCGs3lcqqvQDLgMuASYCP9ruS0RNnQIg6UagD9DfdlciDtIzFX1dBazI6/uIsAI9gIuAoZJaO4Iobsr3NZAL0NPA0GwzlcYj87ax3YtYMCdIOl3SpSnvINs9U77ZWX8OMMv2+cAFwE2SBjfSd4uwPcX2ghxDa+BWYHEW9wDOViSFWQNMBnbY3mJ7YM3YxwDfEvF6OhNW/c1Ab8DA81l1JTBW0lGSOhPhG07LPi4mYnRNrROxB7AhFc3n2XfbZshez5PAS7bXtnCKCoeB/30018IR41Xbe4HNknYSFibAWuDkvL6eWNBXKSKpHg1U5Z84i7+s1zeB59IiXg7cZHtblm0GugPL9jW0vUfSfGCxpDezzf2NyPxCtvlO0vfEInkN8LrtLVk2D5gn6UTCiu6mv6LAtiLcOkubmJtmo8gY9gTwgu3Xa95zMZE85hjC+r+XCDCJpDZECPO+wBDbvxPK98OaficBGxX5TsYATxG7hI3AQqCTpPb57qG2rf2j3bYCrgYG275F0hBgoaQutn8+gOy1YxtFuJXmH9IkFf4xioIo/FP8XnPtmvva4F9HAw/ZfhFA0vFUJ3RpUBq2F0laSiyQg4BPJQ1xhDU2Ff/TtkenQrmU8I8/SOwS6tlVJ/M+GRue54Laj/DTC7jQ9u4s61A37oMmXWXPAt2AYba/qin+Gphte1fWfY1070i6iMgDMAMYl0qa3NnsTLcTwF5iro4i8gmMdoR5R9J0QlFcTuw8Pkjl0BG4UlLblGGV7eUAtt+RtB3oKemTA8heywigl6TVNeNeTZyhbGj5rBX+boqLqXAomEPLWLWEODc4Je/vJtwl9awnzjuQ9CowypEYZTKwCuiZ9Tpl3QYknSppE7DH9mvAeOC8dD01h6XAqLSmIdxT09JK/oBweSGpFZGutVcz+22Ku4gw81dULLALgOvSJXQsscv5SNJJwHzC7Tdjn3JI2gIz6uZ6se09hAttUo7jDGAkMNf2XNtn2D7X9rnEzmKy7WnE3+5MSedluwHEznBdE7I3YPta291r+ievi3L4l1B2EIVD4T3i0LMxl01TzAG6ACsl/QGsBm6rqPcWMCGvHwFm5eGyiLDHC3L30Z265Ea2t0p6FHg/P6P8A7gjXU9NCpiW8UxguSQTOSPGZPHNwExJ6wiL/Hnby6p7ajEjCYv9yxo5F9p+AHiAOAfZCOwgvkJ6nLD4OwDz68Y20PZcSecAH0v6FfiMOOyGcLm9nOM4jvjiaDUHwJFgaRThVgL4CbjB9jZJjcouaRww3Pagg5iTwmGmhPsu/CeQtIKwjL9ppHw00Mf2xKryQqHQcoqCKPwnkNSfsPzHVpSdQLg8htveftiFKxT+pxQFUSgUCoVKyiF1oVAoFCopCqJQKBQKlRQFUSgUCoVKioIoFAqFQiVFQRQKhUKhkj8BUg9BEC8XFkoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the time sample interval (uniformly sampled!)\n",
    "dt = time[1] - time[0]\n",
    "assert 1/dt == sample_rates\n",
    "\n",
    "# Plot\n",
    "deltat = 0.5\n",
    "shiftt = 0.3\n",
    "\n",
    "indxt = np.where((time >= tevent-deltat-shiftt) & (time < tevent+deltat-shiftt))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(time[indxt]-tevent,strain_H1[indxt],'r',label='H1 strain')\n",
    "plt.plot(time[indxt]-tevent,strain_L1[indxt],'g',label='L1 strain')\n",
    "try:\n",
    "    plt.plot(time[indxt]-tevent,strain_V1[indxt],'b',label='V1 strain')\n",
    "except NameError:\n",
    "    pass\n",
    "plt.xlabel('time (s) since '+str(tevent))\n",
    "plt.ylabel('strain')\n",
    "plt.legend(loc='lower right')\n",
    "plt.title('Advanced LIGO strain data near '+event_name)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to whiten data\n",
    "def whiten(strain, interp_psd, dt):\n",
    "    Nt = len(strain)\n",
    "    freqs = np.fft.rfftfreq(Nt, dt)\n",
    "    freqs1 = np.linspace(0,2048.,Nt//2+1)\n",
    "\n",
    "    # whitening: transform to freq domain, divide by asd, then transform back, \n",
    "    # taking care to get normalization right.\n",
    "    hf = np.fft.rfft(strain)\n",
    "    norm = 1./np.sqrt(1./(dt*2))\n",
    "    white_hf = hf / np.sqrt(interp_psd(freqs)) * norm\n",
    "    white_ht = np.fft.irfft(white_hf, n=Nt)\n",
    "    return white_ht\n",
    "\n",
    "# number of sample for the fast fourier transform:\n",
    "NFFT = 2* sample_rates\n",
    "Pxx_H1, freqs = mlab.psd(strain_H1, Fs = sample_rates, NFFT = NFFT)\n",
    "Pxx_L1, freqs = mlab.psd(strain_L1, Fs = sample_rates, NFFT = NFFT)\n",
    "try:\n",
    "    Pxx_V1, freqs = mlab.psd(strain_V1, Fs = sample_rates, NFFT = NFFT)\n",
    "except NameError:\n",
    "    pass\n",
    "\n",
    "# We will use interpolations of the ASDs computed above for whitening:\n",
    "psd_H1 = interp1d(freqs, Pxx_H1)\n",
    "psd_L1 = interp1d(freqs, Pxx_L1)\n",
    "try:\n",
    "    psd_V1 = interp1d(freqs, Pxx_V1)\n",
    "except NameError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now whiten the data from H1 and L1, and the template (use H1 PSD):\n",
    "strain_H1_whiten = whiten(strain_H1,psd_H1,dt)\n",
    "strain_L1_whiten = whiten(strain_L1,psd_L1,dt)\n",
    "try:\n",
    "    strain_V1_whiten = whiten(strain_V1,psd_V1,dt)\n",
    "except NameError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Advanced LIGO strain data near GW150914')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydd3wURfvAv5OQUEIoQSDUIE26VFFRuoJiAbFEfRVQwYYF9KdiQywvYgO7BgUEeYkKCtiwAAKKoDSlhSagdOmEGsj8/ti93N7d7t3e3d5dkpvv53Ofu5udnXlmyzwzz8w8I6SUKBQKhSL+SIi1AAqFQqGIDUoBKBQKRZyiFIBCoVDEKUoBKBQKRZyiFIBCoVDEKUoBKBQKRZyiFECMEEKUFELsEUL86CdOJSFEoZ2nK4Q4IoTIMAnfLIRoYxLeSQix0vA/SQjxrBBirRAiRwixQQgxTgiRboiTIIQYKoRYY4gzRgiREqSsGUKIu0MoY3UhxNJgzzNJx/RaecXpJIS4LNy8Io0Q4hIhxGL9vq0XQsx1yS2EuFEIsdwr/mohxBeG/0IIsVcI0Vr/nySEuNP4bOjhDwkhdun5rBVCfGc49qj+LGwQQnwhhKiihzcXQswRQqwTQqwSQtzglWYVIcS3QoiHLMrWWQhxOtxrVFRQCiB2XA/8DpwrhGgUa2FixEzgLKCVlLIR0AzYAcwRQriezQ+AjsBFhjgngV+FEKWDyKsOcE+wAkopd0gpfZRZhOgMXB6lvEJCCHEjMAa4VUrZWErZEHgKmCyEuBD4AWguhEjT42cAyUBHIUSSnsy5QB6wXAhxO/A3MMoku3bAzXo+jaWUPfQ0L0G7lxdKKRsAe4DR+jkzgbellOcAvYAxQogm+nnfAquBLhZlqwE8ASSGeHmKHlJK9YnBB1gE9AXeBd7Qw8oA7wN/AeuB/wISTVFvBToYzp8C3AakAJOALWgv0lygqh5HAg+iPfR/A/81nD8QWAvkAAuAFnp4ZSAbWKXLcZ8enqjLsxHYALwMHAUyTMq2GWhjEt4JWKn/vgL4ExBecRKADD2/1mgKobRJWnOAu03C++jlzQFWAj2Abvr1O6mXOVWX8TZdhtuBc4Cf9DJvAz4EhC7LET3t4cA4tEpuA7AcaGwig+W10o+N0fPfDCwFGgHDgH+BA8AXejp3AOv0NDYB11hc01+BsXrcv4C+huMPAiv052MaUFEPv0q/Phv0a3OvoYxv6HFnmZRrJ9DWRI4qQHn99xLgav33XcBIYDbQRQ97CBhv9WwYwjYBn+jl+h5oqYePBl4xxGsJ7EdrTJwGEgzHpgGPeKU7HnjIKywJ+A5oCMhY1w/R+sRcgHj8oFVsB4CSwEXAIaAs8AIwQ38YE4C3XQ8jMAJ4S/9dGtiNVpHdjaZEEtAqrHnoFSOaAnhW/11Fr4RqABfrFUIN/dgNwFL99zdAd/13OWAhcAmawlgMlNWPPaKnn2FSPjsK4AVglOHYk2iV81pgF3AhWivvfxbX8CFgikXel+i/L3LF8a5g9Hg/AhX0/x+jtWpd5c4FGuOrAHYANfX/rwJjTWSwvFZoim8mUEI/NtF1HfT039R/p6EpkLP0/7cBiy2u6RlDma8ENui/r0er7Fx5PQxM1p+VDUB9PbwrsNsgwxG01rV3Xg1c8fT/TQ33bKtB9heAMfrvmUAH/X65yvk1cJ3Vs2F4xmehvSsJwCA0xVxGL9P/GeJWxv2ebAIGGeTbDbzmlZeZAngLuN713sS6jojWpwSKWDAY+FRKeVII8QuaMvgP0BMYJqXMAxBCvITbbDEBzezxAFrXdpaU8gjwrm5C6gc0QTORGO3jWQBSyj1CiH+BdLTW3xQp5Xb92CfAJ7pd/RLgbCGE6/wyaJVxC+BDKWWuHj4aeD6Ma5AG7HP9kVI+70pPCLEbTQkK81P9MgX4QAgxA/gZrXVvxZtSyoP671uAC4QQ96JVOqXRruMxr3O+kVJu03+vROvFedMTi2slpfxKt3Vn6vetC5rS90BKuV8fR+mqx7saz/tqZLOU8geDTK4xlL5ole9K/X4mAvlSynwhRAugixCiD9o9N6b9o5RyoUk+aWimG5eMq9GUJEKIUWiNDNBa668LIUqitc4Xod3rT4UQT6A9TzdblMWV9nG06+giSwgxRE8PIN9wzPi7N/CqEOJh4A80BXTEX15CiJu1LOWn/uIVR9QYQJTRbaOZwLVCiJ1oLcqzgHvRWommg75Sys1opo1uuFt2CCEeQ6vk96OZJ77yOvW4MRn9OxHDiyyESBRCnIdW4Qo0m3xjKWVj4Dz0noiXbOEOTv+CVvF4oFdMropkEdDZwtbfC8105YGU8nE02/E8Pc4aPwPGOwy/v0VrZeegtYL/sTjH7Hr6iIHFtdIrmy/RzFGfAh+ZJaBX+n+gVeaz8a9srWRKBEYY7mVz4EohRGW0Z6kZ2n141Cu9HZizCigrhGhqcqyb4fdC4GzgGmC+lPKMlDIHTcn0RWvpHzRJowAhREN9vMHIabTxhINAJUN4VT0MNGVwuZSyoZTyOrRnaYu/vIDrgJ6uwWY9/7VCiLMDnFfkUQog+twGrJJSVpJSVpNSVkOzOzYEDgMP6i0n0Oy3Rsah2YWboFVwoA0aTpBSzkB7CS4i8H39AbjZNXMCrTX2mt5i/QXNZIEQogyaSagpmmK5UwhRXj9nMForPVQmA8f0GT1l9PzORuvpACClXIpmppkshKigx0kSQoxEe7EnGBPUZ5f8BVSWUk7TZUwDKqJVjN6De/n6eSXQlNF/pZSz0Uw1NQn9/fB3rXoA06WUn6GZTboZ8jHKeCGwVUr5Lto4QY8Q5PkauNt17YAhaGMNzfX/r+ot/Uu80ja2qAuQUh7V0/ifa+KCECJZb/23McTLQ3s+n8OzQfINmnnoWxuynwDecCkbfeC3AtrEie/Qnt80Pe49aOYi0Bort+rntESbQDDdX0ZSyt5SygYGRYn+e7MNOYs0ygQURYTWD78LrxaXlHKHEGIyWgtpP7BWCJGLVvkZ+RztAR/tMngCrwOv6N3jf9AGve4XQnxgJYeU8lshxNvAAqFNM90J9NcP3wK8LYTYgNbiel9KOU8IsQCoBywTQhxHG2w+6qe43wkh8gz/F+iyu2SQQoieaDNIVuty5KO1iv8ynNcfrdJZqF+/RLQK5HwppYd5Rk/zITRTg6tCe1ZKuU2XpawQYj2Gyko/77QQ4kVgvm5+Wo5mlnkezZ4fLBOxvlZj0UxUfdHs07OAu4QQb6GNG0wSQlQH7gRu0xXav/p1SRNC3CKlnGRTjo/QlNlvQoh8NFv9bWi9j9XAX0KIPcBUYLdunvGLlHK83nP9WAhxFpriWAnch+d1/Q6tcfKdIewbNGUYUAFIKf8WQgwCsnUFvQe4QldCs4QQ7wGL9OdmA25T311o13cYWm/2dinlgUD5xSvCXY8oFAqFIp5QJiCFQqGIU5QCUCgUijhFKQCFQqGIU5QCUCgUijilyMwCeumll2SFChUCRzShVKlSnDhxwmGJCjeqzPGBKnN8EE6Zd+7cuXf48OGVzY4VGQVQoUIFBg0aFNK5OTk5NGoUX/7WVJnjA1Xm+CCcMo8YMWKr1TFlAlIoFIo4RSkAhUKhiFOUAlAoFIo4RSkAhUKhiFOUAlAoFIo4RSkAhUKhiFOUAlAoFIo4RSkAhULhCCt2rWDRtkWxFkMRBEVmIZhCoSjctHq/FQByuHIxX1RQPQCFQqGIU2KqAIQQJYQQy4QQGbGUQ6FQKOKRWPcAHgLKxVgGhUKhiEtiNgYghGgA1EXbK9YqziBgEED//v3p2LFjSHnt3buXnJyckM4tqhTqMksJJ05A6dKOJluoyxwhClOZM9MzASIuT2Eqc7SIVJljogD0zb1HAncDL1nFk1JmAVkAWVlZMlRveMp7YCHjscdg1CjYtw/S0hxLtlCXOUIUpjJnf5INwJRGUyKaT2Eqc7SIVJljZQIaCHwlpfw3RvkrYsm0adr3vn2xlUOhiHNiZQLqAlQTQvQDGgGNhBCXSykPxEgeRTSRapqgQlEYiIkCkFLe6PothBgPPKMq/zhEiFhLoFDENTFfCCalHBBrGRQKhSIeifU0UIVCoVDECKUAFApFZHjkEfj661hL4Tx//w2nTsVaCkdQCkARfdQgcHTYvRs+/DB2+b/8MlxxRezyjwRHj0JGBgwaFGtJHEEpAEXsUIPAkaVPH7jjDti6NdaSRIfDh+H48cjmceyY9l1MejZKAShih+oJRJbdu7Xv06djK0e0KF8e6tePtRRFCqUAFApFkWfNv2s4mgTs2BFrUYoUSgEUZU6cgPz8WEsROsoEpHCAvDN5NH2nKddeH2tJih5KARRVTp3SnKk99FCsJQmeImr6mbZmGodOHIq1GMFTRK+3XU7nayaun+rEVo6iiFIARZWTJ7XvDz6IrRxxwvp967n2s2vpP6N/rEWxT5z0sCTFW8FFEqUAFAobHD11FIAtB7eElsDIkVqF7JpF4gBHTh7h1JniMR9dERuUAijqFOXufRFqoYbdynzjDe37kHMmpHIvluOVha8EjliUnxEbCBx6jiZNgr177cU1u6ZSwltvwcGDzsgTBZQCKKoUocrTLnM2zyHvTF6sxfBLyJVNhCrhvHw/16sYPiOAdi2zswumtzpiAtqyBW69Fa67zjN8+XL702gXLoT77oOBA8OXJ0ooBaCIPiaV4aJti+g2sRtPznmSf4/+y5IdS2IgWGDCrmyKa6UcTVatghtv5MjLL3Dtp9ey5+ie8NN0jakZp5GuWgWtW8OTT/rGN7uPJ05o3/v3hy9PlFAKALTl8n37xlqKuGZ3rrZoKWdfDm3HtqXd2HYxlshhzHoA27drFcnPP0cvz0AcPw5nzgR/3vz5wZ8TKvo4yoQDc5m2dhrPzXsuMvns3Kl9L11qL34RNLUpBQDacvnPP4+1FKFRBB86M4wt678P/R1DSfwTtr3Z2HL86Sft+913w0szGPLz4Z133K1Vb8qU8TWDBGLWLOjUKXzZigtFqJenFIAidhShF0WGq2hjoajNru/UqXDvvfDUU9bnffFFcPn8809w8R1GOPQcvXwhJGWu9z0Q6Xu3ZIl2r1avjmw+JigFYMTBKXpRowhVogWYvFCOzeQo7MT6fh05on0XFjv1/v2wJzwbvo9y9tfDsWLGDB65FE4ba8Ro3atPPtG+v/kmOvkZUArAyIgRsZYgeMJtnaSmwqhRzsgSBoV9MY9PK1PK4K59MTHVOU6lSlC1KsycGfo4wtSpnv8NPZy+n/bljcVvBE7j0Uetj82eHZpcRQClAIzk5sZaAvs41TrJzYXHHnMmLbvYld3unGwn2bFDG5z1wqeV2bkzJITw+sSiB+CA8vn70N+UeaEMa/5d407Tla4Tyu3qq4MeR3BdSXnksO9BvYfz+drPeWDWA2EKB2ze7Pnfah2Aznljz+Ph7x8OP98IE3cKIF/mk3uqEFT0UtraVUhKGZT9ed6Wefy2/bdwJIs8dk1AR49GQRgvatSAmjUtDxf0BMxaq6dOaT7pzSgsYwAhMm3NNI6fPs7YpWPh33815ffWW46lHxTelbFORFVrkC61f9/xO6/++mpwecTgGYk7BfDYj4+ROjI19krglVegZEm/rdxjecdIeDaBZ+c96w70Nj14PTSdP+pM+w/a25NhxYpgJA6dL77QKiOX/dmFoYIq7CYgW3Tvrvmk94dZpVzUzEN/67O0JkyIft7ffQfr1ul/PK+bxz8LJeEYTvbkYjguFHcKYNKfkwDNj0pMmThR+3bNNTbB5XnyvaXvuQMvvVRrfYXy0Bw+7Dn416pV8GmEwrO6Atu4Ufv2U+FFczB46pqp7D3moJlpwQLrY2ZljvWAcLT44gvrnlGw9OzpE2T6NM2dG1Y2t824jaxdIe765eR9veQSd10RAeJOAbgotC3OY8f8+/j/8cfQ065WTRt0M2HJjiVMWDEh9LSLGLtyd3HdZ9dxzYTLnBn7CWWT8CLU8jd9X+zKf801MGCAswIB3kYfx6rdd95h/Irx3LlpTHDnReJ+/vgj9OvnfLo68aEApITFi4HotDBP55/mrJfO4uM/P9YC/v4b/vor8ImnTkFKCjz4YETkuubKY5z1iPmxdmPbMWBGJF7SwsnJ09rS/60blkCvXuEn+IjFhbVDsC3GP/6AtLTgzpk8Ge65J7hzTBCn80Jr4W7ZEnbegXCs+r33XqdSKvTEhwJYvhzOPx+mTy8IshxYPXEC8kwcbAkBvXvbyu7IySPsO76P+769TwvIyIB69axP2L4dvvrK7Y9k/HhNxiAe6RW7VvDTlp/8xvmiMewrYztJ2+SeymXmupnBnxioIomgiWTTgU3uP064MVi50l48m2XauH8jFUdVNHc//dprcOCAfdkA/vMfzxXHobZWs8YGjFLQ8AkC8QwM7RGCPIC0cUkPHD/AnM1zPAPXrClcvbDVqzWZWrZ01GusP+JDAbgGWje5X/pdubvMlUDp0tCihXk6M2YEla1p+q6Hznjsggvgyist07HTa2n1fiu6fNQlKPlC4aJxF/n0FO6YeQdXZ1/N2n/XRjx/p+g2sVtQ8cM2GUrJxy3g7T8/tBX9w2UfcvDEQaasnBJcPhs3wuuvWx8PV6kaTV0WlectX9ziG2ijoh19gfUxKWXA5+tEEiyobX6s1/960W1iN47l6Ys9586Fpk0hKyugXEBw1y3UazxxIgwfrvXwvv8+tDSCJD4UgAltx7b1sHm/3wZaVfxU+5OTE5lMrR4611L6SLZGRo92//79d9qNbUeJp03ijRzpN5lf/vnFZ6zA1Zq2nFlVmFpZYWKpjI0v/YEDmkvgY8f4bftvfLfxOwBuuQYGz3+UXbm7PM+Vks0HNps3GPbsDizUc8/Bbj3exRdrJkTjFFoHr//BUpA2qwujOgSpFIOcRulN1tIsmrzTxHTbx51l3b873mZ+/so9Wg/tTL7u6M41k2j5cnsCeF9D/X/embyIzyg8lhS5tONDARhvnmHA76etPxX8vutKWJFkMSPkzz9tZbPm3zXuhTKYrB5dr/sZWbbMfyvB65jfF83Oy/333zB0qPv/oEEs2bGEM2Z3//HHA6eH9iKdOO253D5ghRBsy6gQKg5bld6IEdpWnWPH0v6D9vSc3JNRbY4XHB6/XDPxua7H/NJ7qPtGXT764yMAZm2c5d7pa4yf1ryLp592DxQaZ9tEYMrp4hpwIO8wj10Cb5+9tyDNPSkBTrRrIrPA5R58VAdYV0mr9F2l+75+WEkHh35N11U4gxghSH4+mdSRqc6l73V/pjSDlCdwxuW1CfGhAHa5W1zicAjTP20OEjZ9pylN32lqL81wK7dgKlM/ra8V6fBVw+Czv2HqDZR+obQmSrAD615lv2rKVfT5pA8Af2zz43r33Xe1ckfBPjr7r9l8vymIbrjxfrjKZxj4fOwit7J8fM7jHvsdrC6plWfxtsX8/PfPXDb5Ml5b9Jq9vFy4Wvx2nyu7z8+ePdrmKzpGe/u6VG3MavGpzVT9P3vJefNrTQIrDwOzGkB2c6ge7iLbUN8//bxF6Q5sXLRuHbz9tt8ornfT5S7daeJDAURhBoIpZoPJ3th4EcOauTRtGnSztne3uguuvCmEZNdO8wmzu2L5eGI+d14B+09pFd+X678sOLbl6DZ3RO9r49pW0bhphxUXXwzvv296aOcR99oLsyu7qgp0n9SdHh/7jkoGdS/efNPykMss5M2/R/+1n74/TBTS4ZLQZ+/bPNMZ8i16Mv+WgUFXumdJccUVcOONyFyt4bS2stcJK1dy/v7QfUldeAecf0fIpweNRGomXteMqFDt9bkmq9RNnv99X2pm5afmPEXfT732HDnvPBg82FOGKK8NiQ8F4CLYi/vttyEtYPnfyv9pP8weEvC7ktcU7xW0VmmZce21EVeAwbrj/ajBMbLawpNfDQktQ0OZP1n1Cav2rNJMDMaewc8/w113aVv0eTF1zVSfMCMjDC5pjq/5k7mb5/pVbl+s/YLq7eZzKtF+EQBmH1tFzaFwVLh7aN4mpqXVYcFWrwVm/haVGY953Zd328L0E8sZ0Rmyy5hPS/6/S2FsG8hepbf6XeNTZmtTpLSeMGHCqJ9HIUb4PiubK9pOwhkCjHN5I56BgYufsBV3WzlgrXuwuvaiG+CPP3h+wfN8vtZrzxGnFseFQXwpAC/8tua2bIHLL9f2CQ2SpTstzBgmFaWUkpOBKo4oOqm78HZC9iVkd1DQFUvOnuM3nh0yp2XS/N3mWkU0xmThjsFfzR0z72D0r6N944A2PmNyf+588ly6TuzKxv0bC8L6fNKHvte749w/6352ljzFaxfAhJb2ZX9k7xS2l4OcZHdF4K1opjWBjhM6+nWZvKusSWAApexSOlsPbvUwL0Ry1OWx2TacDtpZLxMDPtj0mcd/s6v7zYHfqDUUZqa4e7HHkuHAv+HvlxCphatxpwBst1VdNtUNGzx91kjJx5emc+KD9zyim7UQ7cxPztoylVJPwdby+LTs/JpUwhlD8HPur7Xgvsn/CSo522aRV02cY23c6Bvmwnt+vvfsrA0b7OWr8+HyDxn6/VDzHovFJiirdZPH4ZPuSnp6znQ+b+Ibd1h3GNAb9/X1UwkbX2hbd/Kaa9y/vdK9/jrMewAFPnOsqfN6HdJfTfcrX1Tp4txU5gO93SY8j2f0a4OLh/fe83UnHSJLj2rP4+/VPcOf3Gxv6u/xEr5hkTYIFX8FEKpHSYtKctbGWdzSYTePT73bI7zrxK62k96YcJA9JbUW2KfbtIHGjcaFnV4vuJB4dCu9Mc48MrL32F7EM5hOnfOH3BhcxWqbj10LhAzXdts206gA3GIyn9zIOef4BJ3ucQlihGB4Z3joUthuNkHDsOdtoGrOW4nn5Qcx+BegFe6qlKQwKgMLiXZ7DgIaW/0HStnL2xhilY+VxGaNERGijvDbsHGZRUaNCtsePrCkexDfo7z79nlGnDzZMo2Zvo8YAPnBzME46Z4BZlp2IZjWGMo8CStKHbSfsAMUfwXg1cL0eGjz82HFCvNXwWJh1sET2g3a6dXtDrQKV8tce2oaVPuMmleYbD33tbnzKSmgybvNeOxH3y709AZnLGceLd6mub94+cLAovmwZo3H9Nen5z5tar818tz851i9R9vWLvdULvnS3KeRx/UOp8VneJnyBcyqD2tXaL6Snu0Mr10It/ZBWwhocNUgHjBxtWGzsvF5gUPtiR07VuCxcnz5zboI9muVM9Y1teU5Iy+ynTwDZgxg1KCmHjPovJlV9QjTGwVOq3cmPGj04SYlfPSR/x34bOxREejK7ysdWDYX4hnz8Ktv9Py/6dBmbukDt18VRNrfuQf8xy0f5xtBSr5poP1cutGPQ8EIUPwVgD9+W0z7t1tR12y/iEgMnBpezrwEk8d32TKPv8YKYW2lfEb9YphtobvkXVnF+jWwbE0GqGgkaAvWzj1XC1i6lJELAg+cfbPhGy4afxFHTx0ldWQqj/7guctSbjIcPXW0oFfmZPd2RTosrqmbQwzkJQL33w8vv1wQ9pnNmbrgW8m4FhQVYDbF1swc483ChT4t6GD3fvDA21xmMmh70Fgh2ljs+Eaqu2f58p++G9dvSD1Fn0z/aRwuCTMawevnu8PknNnQvz889JD5ScbZcw6tBfFerPVsJ63Sl0HmcdvcB/n4XMj3qjlP51tPtTbeZ+PaowIM+ctTJ23L4gTFXwH4qew+Kr2e32rCFn+zELzOf32xYWGOnWmewWB8EA8eRD7/nHXcxo31c0LIJ4DnSp+xi7ZtrdcSfPaZR0V+6swpDh/Xekkf/6pPw9TLlfo4VHyxQvjXrWlTGOI5g8jVJTcdd/HymzOvjvt3QfRA5hqL4yv3rGTbYS8zlus+Blj9GrIC9CerK2+LKbAFBLkn8J7jobnNNm2Fu8w8Vr0L4zRff+NVJax7J4EYbux4BjmWZMbEP6xdNvvcrb//9ij7yZ3bGNfaffi1C2DZA9dT9/W6TLY/ySokYqIAhBApQohsIcR8IcSvQgg/XkBiy9yzcnmijdsut3j7YvfB/zNf/fLdJvM53kHx2GPwrjbQvNNgx/Z+HZalW7uONg58edgs/bmbtkzLguuvR+zw3NNAntErviNHkFKyqtyJArnzpHHKYxAcP+7532zGD752aQkwa1YwOXmeH6CWbvFBG99Af6YNnbX5Fgt7fl1oHu6q9Fu3hgkTzOW6917Iy2NNZZj3+lCmNNPDQ2xF7y4Lu4NYpGUXaWEaDJb/lfQ/yG1bwQbhpsJqwsPxPIOdP1DGGRmaa3ad/dPd/p4W14CHekCbO2Hzwc225QqVWPUAbgfWSyk7AoMAG+vdQ8R7QDXI07t22cp/W1rMw583zzR4xxHzhUq5p495VMQS+POQYSzA+KKePMnKqr5pvNBRM6W4mH6Oxct96pSHCWh+hnk0MzxS/OOPwCe4dojyYlcqjF8xnubd13vkb3fwcPDl8Pz857U/JquxjbMtXEnmeC1UsjMTyx9/+E6QCUy+tszqg9bWUf6XlMPas7wCDx9CWig2pETu2cOR1b6+a1bpz8muie/Q4D5oei90HgA3XQsLa4Ugv86ZBEgPcYWvGa77brUIDeCtc0/S7TvDykSH3YGUf7E8f3tv2mYzj33H9jFvh7mC3nxwM7lnjpse88DfpAfArD0RSUwmHkWFNYBrdPYwYDoTXggxCE1B0L9/fzp27Bh8TgcOsLdZM3IyMyEtjcsSMzlsNmtCJ8fLppmZ7hnu+t9MQE4SBbbUzHRfY2jJ05CTkwOZmZxKhJGJi+HdTFwxJw+F7vrsn5M9ICcjAzIz2Vc+me0Nq/DF+e64LlbfBC/1gJtWesrjzbSVX9K8TnMy0zNpkARHy7vjftcJMtOtJ6ynVoS1mbC/NFR68UXIzOS6ar52T4Ab34PWxyBDdzOdlJDE9q07C67Hke1HyEzPZPN17vzLdtR+128NOanWZdh3Fexbu5qcKjlQtSpkel6ND86HTP1eHu4Fzco287kPtUr73jsjJSrrxytVKki/RRMoUcUzXsqhFJ+0czLN08wpdQ4b7slkdgt87p8/6ic2JrFrpjXz5j0AACAASURBVHmal6bxy4LJ/PhGJkMXaootM91d5pxMWFQT2taHtobz9l8BK/fs4sDtnummdNSeTVeZcvTnuPG5vvlbldMfgycO5obbMimR73vue8fWcElmJlSp4nNNf7kRqkjIydTdeq5bBwkJNJQNC+Ka3WdveQHaN4KqJubduXd4vtciLY3MdHfN++kD5uUdNWOUZb7bN28HSpCZnkmN8zzPr9/a/b/JwZIsvSUTAZTRLaFHkgNf38SjouAeOUlMFICU8kcAIUQtYBJg5pcSKWUWkAWQlZUlGzWyMeXAm9OnyVm1ikbZ2dCuHd9uy+YfP9u2Tsn2/J/dyDPc9T9zJTwxDZiidd+yP/E6ESh/Aiakv8uChdksrgnZDTyPl1sP2eW03yV3wYRS6ZCdrc1IOKJ/TKh4HJ71ksebzPRMzq18Ltm7srl8PVy8CLL1lmJ2NSCA+fTcTTC4Fyz8AC7YBlPrwymLp6XkaXBNky+TUJKRZe4je5cmYHpGOtm7smn3C2TrM6cuXqCVu+IyaPS1dRlcTGk0xcMfjYtvqsE2/V5m14DM3MyCfI28FOBaTcoG2rQpyOPP62Cq10Bxl8kzye7madqZkm2e5pwz07lp0wmyq/ge80fFhKupOWcG2Sabtk35vj63NPiJJYfWMGQ21DgM2Wdr9zl7VzZTsuHb893X2EX/7+GGit+yttYhj3veeT40em1KwXM7ZWYreOQRRq3IJtur12RVzkAcPgAj5kJ2Hc/wOSXKc1/2IUhIIPtpc3NQwXt4+DBHRjzBstPLCu6tq8xWuM79LRnm1PU9fv7PkK2bVSd/AgkNGpJdxd0Tz66I6ftR5UgVW07Z6i+GbEMjs9JSyNabuIl/wuO6XV8+o33vLAvZJnIa6VPvCkKq/wIQs0FgIUQ/YCIwVEoZ4uabwfHV6TV+K39vfrbwLQ6aQ6r9AaaZHSoFW/t0peNt2hJ7f3zUEku7dsgYTFQbg9xA6rlOoZ3HiZNwww0Ff4N1ExEM4Zp3gsKGXd/FnsQTjAlxVCvzOosDGzcijOsXbJb9n/KwNtnXed5d3rOcH33UJ064BDKs7ChnoxDffEO5rzu4XVM4wLMGVx+t7oSTCc6MSQTLgz0hL8bTcGI1CNwb6AZ0l1IuCRTfCfIS4Mrj9lbkubjYwre4i0kt4LPVn/mdH9/sApv+xoPArlV0+levFPy+54rg8thtaEmKZ6xb/97yHEuGX077Dl7tDWMnMqv1BE4x5nzI5RSHS0Y0m7DITYZjJ7QuoQS+tOnB9V9/132N+QLCaBEpH0CfN/Z/fL/hmvyZDn+lODybzwurMa/Xz4dXLgzv3QiXWOmfAUBz4EchxFwhhMnqCIcQgsPJkGxqZLJPlsXgzJhfTNwbGMj1U6mE2oC12/qbcU54+dg91zvODSat2JEX+4a92w7fgVATEp9N9PGXJIHt5WwIZ4MhPaHS0acpP8w6TqirXp0i9XFYfVqbXLAz1b5C97titWkQCyIiwekzgeOEQN8bYFO0HcwZeL6T539/t+Dx7tAi/K2aQyYmCkBKebWUspWUsov+CdDWDp2Fe5czOpSVsF7cabIweGwbWLhzse+BCHOoFLZWYOYF6Z0yWhjr0t9q2DvnaLLmaM3lLyVY9xaBOIX/yuhIFHoH7+6wt+XoqiDGFoI1k01oFVz8UPI+lq8tdrIjWqh6t77Z4k4r7F4jp9f9FAJiNQsoaqw87MfZWJB4LxdfHeQgnzfhNCr7ZELuC/7jHPIz28ku4UwjBBybxldJN1H/VgOG/mrhATOCPBTihuWRIJgrGtTV9+eYz0Fy808EjqRzOgpN1P3J9noie04dCBzJBKcbK05S7BVAYSarbeA4/jDd0tFh3m0X3vlisa9r6UCD5/7YnQIN7g9DoGKAWcv6hMWb/LR9H4XQoAE8E4pE5vhTPqsra9s7BiIa4/wXXep/bn64GNfzRHplb7AUf1cQQfJtfV93roWVYDwSxowNvk7vgqqUosyGIGY9BbsBTCSZlwFDewaOZ8TboaHT+DM/XXAHzPRjxmx8r9bjDqexUBReD7tEagiq+PcAgnwKLg/OFb7j/BhgPrCRT2M8hmeH1xr49yHzbCe/h32I9NTPhvfDtavtxX3FgbElp+gZwHO2GY9b7xTqwbgQxwVOJkK3fubHjiX5P9e1ojvUvYYV9ij2PYCw9tONAZcEsQGZ2cC0GV+HsOl7tPgryHUGXwSY4mdGsK2nFTZXvcZq2qhTrcGtFezFuzPIKcQu/kiPwXaPBpZWCxynqHAmPzIzpoq9Aojx7D1FIeDsB+GLIBZRbjRZiWvGshhVMHZmgNlh7tnu39Maw0GLSQOnQzR1+Xv3ovFeHgzDfFTY+HqLAw4mTVAmoEJEMBtYFCZOFvKnaGsFuCYYpzw2+aGe82naYUUEFM+1N8Alm5xN02pgGsz9SimsidRiyGJ/G4qSCcjHS6FCEUWcVmg7HFqop4gcxV4BKBQKhcIcpQAKER8XsjnCCoWieFP8FUAEvVE6zWuFaFqhQqEo/sSBAoi1AAqFQlE4Kf4KQKFQKBSmxIECUF0AhUKhMCMOFIBCoVAozCj2CiCSWxIqFApFUabYKwCFQqFQmKMUgEKhUMQpxV8BKBOQQqFQmFL8FYBCoVAoTFEKQKFQKOIUpQAUCoUiTokDBaDGABQKhcKM4q8AVP2vUCgUphR/BaBQKBQKU4q9AjiRfyrWIigUCkWhpNgrALUpvEKhUJhT7BWAQqFQKMwp9gpAStUHUCgUCjOKvQJQKBQKhTlKASgUCkWcUuwVgFTDwAqFQmFKCbsRhRDVgHRjmJRyueMSOc3WrWoxmEKhUJhgSwEIIV4AHgD+xT2zUgL1IiSXcyxZAu2ax1oKhUKhKHTY7QHcDJwtpfw3ksIoFAqFInrYHQPILaqVv1QbwigUCoUpdhXAWCHEM0KI6kKIiq5PRCVzClX/KxQKhSl2TUCj9e+nDWESSAw1YyHEm0A7XYbhUsqvQ00rQE6RSVahUCiKOLYUgJTS0emiQogeQF0p5flCiOrAQiHEt1LKfCfz0TJzPEWFQqGILhGaze5XAQgh+kkpPxJCDDGVScrRZuE2aA38oKexQwixB6gNbPHKfxAwCKB///507Ngx6IwqNuhCs7JVyEzPDFHUokmzss1UmeMAVeb4oFlKM3JychxPN1APoKH+bTaPMhydlAzsN/w/CpT0yUDKLCALICsrSzZq1CjojL4ZPZdVqU3J3pUdoqhFk8z0TFXmOECVOT7ITM/kiUZPOJ6uXwUgpXxC/77N+5gQoksY+eYBpQz/SwMnw0jPEqlMQAqFoqgTCxOQCyHEBcD9QIoelAC0xWtlcBAsA+4BXtdXGFcB/gkxrQAoDaBQKBRm2B3c/QBYj7by92u01vqroWYqpZwF/C2E+BmYDtwnpTwTanp+81L1v0KhUJhidxpokpRyuBDiNLAY+AhtEPflUDOWUg4O9VyFQqFQhI/dHkCCEKIkkAM0k1KeAM6KnFjOoXoACoVCYY7dHsB44CfgKuBXIUR34ESkhHISWbVKrEVQKBSKsIhUO9ZWD0BK+QIwSPcHdDuwFSgSE3Fl2bKxFkGhUCgKJXZnAf0ppWwBIKWcB8yLqFQOoraDUSgUCnPsjgEsEkL0jqgkCoVCoYgqdscAygFThBBrgB1APoCU8upICeYUaktIhUKhMMeuAvhG/xQ58qVSAAqFQmGGXQVQTUo5yhgghHg0AvI4juoBKBQKhTmBvIGeAzQBBgsh1hsOJQIPAqNMTyxEqHUACoVCYU6gHkB7oD9QERiMezpqPvBu5MRyENUBUCgUClP8zgKSUk6UUnYFvgee1X9PA/5CcwdR6FEmIIVCoTDH7jTQw8CtQohmwBBgG5qDuEJPekK5WIugUCgUhRK7CqADMBDoAbwtpXwWqBkxqRykbVJGrEVQKBSKQontvX71/Xo7A/P1oKRICOQ0rZOVAlAoFAoz7E4DXSKEWIC2AcwyIcR44PfIiaVQKBSKSGNXAQwErgEWSimlEGIxMClyYjmIlGpTMIVCUaQpfToylZgtBSClzAUmGv6/FxFpFAqFQuFD9WOJEUnX9hiAQqFQKIoXxV8BCGX/USgURZt8EZn1TMVfASgUCkUR569yZyKSbvFXAMobqEKhUJhS/BWAQqFQKExRCkChUCjilOKvANQgsELhQfXDsZZAUVgo/gpAAUDF47GWQFFYqL8/1hIoCgtKAcQJ9dRLHzfcvixwnO6bIi+HovBT/BVAgFlAJSIzu0qhiBnN9gSOU/FE5OVQFH6KvwIIgBohUBQ3Aq0ZEsD39aIiSrEn/Qh8PTnWUoRO8VcAAQaBI7TArtBx1rFYS6CINZcdrFzw+1CpGApSjCh7CtrssD5+doWzoydMCBR/BRAAp3oAe0c5lJAJVXPDT+PjzwPHyRAVw8ojIR+uWBf6+QOWh5V9oeH+RbHNP00N+BcaEhOcceJWrXQVR9LxJu4VQIJDPYCk/NDPbbjX//H3vww9bReVDJWC1bhHGZEcVh77XoIO//iP02qn9bHEMK6hUzT6F/a/GGspwuPmlebhrsaOBMoVsjGAxgeTWDTWN3zp+87n9dVkOD/Ac+oUCcKpKlb5AgoNKem0xfqwUyagcNIZsMIZGewiLbo9IkYjIuULU2VUobzpAGmZU75hVoq562ZnRQoWy0ZNCfcmfrcVst6WRNJ+u294az8NhlCpchQaBWh0OcXTHZ+OTkYhUvwVAP5b+U5VeQklS4Z87qM/OyREmKgBcaBaNdNg72sjn4FBS82TuHodHBzpqFSOIBLd5ojEQjb2FU1xpIDSp6OT11XnXOVIOk3KN3AkHW+KvwIQwrR13lZvbTi1QEqc0yj0c50RwYMU6btls8vEYvWyRaMH4LoXycJ3L6JItPacItgeXvmTkZHDxcz/hXe+U6ZPI+1MWvB2seqVRooXf3QuLTPRm5fKQA6XpJZMDSqt0nkm6UvocFbr0IQLQPFXABZMmK59lzV07e/9LbA93oqE7E/CF8ohBILbzrTwCT/4Ihz+r7/zokdJ4VZQrsq179ooCmCBHSV4jgPmgxaVmoZ1fiApR8+Csw94BZbQlW7jRiQ+8lhY+Zsx6+PQz23U6GLA3NQWCon5UH+f+TEhoZxDCjp7qsWBhNCqVjPTVNlTICLk0iZuFYBZC+itb+CG1SGmV8dzute3YbwM4VI5pbJpi6rsKUj184KF0gO4fL3ngGKghqVxINJfWGHDeGWSHFg8WKZE6bDOD9QjeXARNP7XM0ye01A7t0pVWlVrFVb+ZgQ7+2jqJ9BvBXw3CT7uNwOA2oe8Ii3zXdY89PyhttKPRoOmzU6Le1G2bEjpmSUVyanqcasArHjmJ5gxJfjzvCvPBhatj+gR/OMvknzNRlZ8pS9+MSoaOw9qYV93YdXSMsrtSMUSZiLG0xt5VfT8YzHFRR8DkEiub3p9yHm3dMhU13et1hO/dBPWppJWoSmqaD5nCXffHb3MHKb4K4DMzIBREvPhyT8rAFrPwNU9DGZaonfFUfsQdEhryU/9fiqYymb0x/PxNEPk5OCmX4Y7s8DK3tqh+eW20wi5/kr1fdFDeVk/PN0rqPgVbLZOL6p1UcA4RnHbbwtKjALMelsv/hBaWj7jDTVr2s4zFJoYFM6g/XVZMGAB699wJGlHxwIsn6v/3OxcJkBCannfvIuIF+KoKwAhRIoQIlsIMV8I8asQ4oKIZlje9+YYkQJOPwvPrUgLK5tE4bngIykffu70EZ3qdKL1Tjj2PFxlWCSVYhzsCXIG0T3t7uHirX6Ot72n4Hf9fdqMFTuM6fl6UHK4SEmwL39pNDv0o+l9C8JCeVXMxjj88cs4e0v237z8TZg92yfcSsZ544MSwy+P/AIX+bmvRmrYcOkcqYHVDIOZRgAX1b6IBnrj5pXv/J/bvrz/sQ8nW+6Wxb//ARCCnDeDT7Nd9XY+YQkmOVUu41513b1udxodS6HZbvv51Kvo9tXx4Uwi5tbedypG5LkdWC+lzBRCNAc+BM4ziyiEGAQMAujfvz8dO3YMKcOUJheQeQIu/BsW1tbC9l4FmelQqRzkZAKdOsG8eQCcqKAdE9L+S7Ru3TrW5mQyorP2PycT2L8fcnIKeiF16mnpAtDd/TunTzKfHDjFF40tEjfEBTiw7QBZxzMZke4bdehCKHv99fzzzyIyEzNJc5UPoG1bWLKEG8xnOrJpwyYy0wP3mADOXKLJVD8JnrjkLnJ2rmLnlw2pvGUemfmbqHEY+i+HFzrpWW+HJTWgToXSZC06Dq2vZ3W+NhevdB4cT4J/rvEspz9yTtdg3JZMvq/vDmtWtpmH/PfuqcPbVbYAsO9KOFMq8H39a8NfUL26z3UoeRpO6m9LlbKGa4qvzK5j/spSs2x1MtI9HfKsy4SxxzG9r0buXAIlu2jpNyvbjBoNMzm7nHYsIR9y9GeuTRMoX8l9XsPkhqSmp5JRKoOcnBzb99qbm6tAQq2qbMrbTf3yFQvyA0ir6b/cHSu24uzS2iKEnEzgvPPgnHMK3pPudaFliju+mZyHDxwOKHtCPqRVgFYpvscObz9Mzg03APafNxctq7SkXr77vuVkQn56LTLJpFtGF2ZvnQvA/7X6P+26AG+2exPWfQwHNjHiXP/pVy0LjcrCoNaDyFqWBUCzFrA3N7cgPScRMsp75gohugMbpZRbhBAZwOdSyjaBzsvKypKDBg0KKc+x341l0KJBDFsAmTc8R6lhT3FGQJPB2oyOnLeA/Hw4fhyaNGEeW+k8QDMBnfHTR7r3N3hbV11yuNSmnD6j/38GWL4cWrYs0N4PfX4nr/2p2YOmT4HeN2oDdWsmpvJ9lSP0uMU8H1dcF955Gdn1MlTNlax65y2uXXsfb3yr2Vi1E63Pc6UrRmiyNi9Zm5Un/7Ys+9eTodfN0HMDfPux+xka+fVjPL5kFI/+lsyL35wqyOvHj6B7P+h6KI3Zo/fD118jftfMOJWOwb4ymuzz6kBWG5hd1zO/addPo++n7l6DPPEYiye9yPkD3XEy0zPJ3pXtjrPsKprUmsnayrD6bdh8Y0+uSJvFZRvgW4tp1XK4VhbXdXDRcC+sP0v73Xw3/Pmu+5j39XT1uKyuM8D56W1ZtGtJ0OcZ4w2+HPZdlUn6F9mM0fvR696EhnslfPEFt0+4hnGG2YO9GvTi6w1f0zGjI/P6z/Mpoz++yIY+ep0r/7yGdy+pwD27x3Hn/rq89/qmgmf81Qvg4R7W6QyvdzsjNn2opfN8CcgzdIWFoPG9kONuPHs8ky4W915M++ntTdN33dsSZ7R9D4xpuTgy7Ahl9TEHq2vteia9GdByAONXuLt98hngmWdg+HA4fhzxUpkCuT249FL44YeA9/bcXfBHOiy/czmt3m9VkEfOwoU0uiA0Y8mIESOWDh8+vK3ZsaibgKSUP+qVfy1gEhDVpXItStam4T6T7qEQUKYMNLC/4OKV721GbK2/hWW05sjwn9yHbA0W16zh9/BthokSrpZticQkct4yVP4Gfhqvtc790TU1OBNLAeV0k1sbT53u0+I2dGkLuv3LlnL9arjQZAzzmsbXMKLNw/7T9IO85hro1s0zPwPrBq/j5wHWK/J+nAj3LTbI21CbUUNWln0hHCbjoPadeONNBWENXc9Tnz688VkuEw0+oGSI86zGfAu9jY3PadOgVi0ARN0g3Yo6YMrwN5ZhnMBhFuvSepdSNjnwDB3j9PDLNkCmhXsND0qEZlB5dk5IpzlCRBWAEGKgEGKB12egEKIfMBEYKqX8OpIyGJGaUP4jffIJPPUUAOf5WdjSfDeUslhNeNZRV4b6C7d4MZx0j9Sleg/aSf+vZkK5Cn5F/nCmSeCFF8KwYe7yGsrdaSuMn+EZ/dubv/UMSLTnxMrqaoY04Fjdv6IzG0D2i9e9LpOoucA084zasFJDOtTuYJlUrcMwbIH2+7KNQAtdQVbwf2+8ubyB/YF2u4hq1U3DU5JTuGXk1/TW11ekp2j2jpZVWwaVflmzqcNVqgIg69U1OehmkGcnBw+LQ4QHSv2OJzzxhO10Uk9C97/0NM2ea1c5/M2i82Np6WFopEV7klxEFYCUcqyU8mLjB/gX6AZ0l1IuCZCEI3jcNK8FGj4XPC2toKWYdAZurWren7V6uE48B9tf0//U1gccSpSwnuljY5rbZfvTGLLI98Ez7u2ampzqKVdSEvz3v+4BTb3FZkZf0ZSe9Xt6hMlgK9sAFMhVUvdDXK5cwbFvJkO/daWonGLSX7fiKhtL7A0vnQA6V2zNO1/BW3vcveFFY2Fi74m+p3p34YFqubDtVXjBd4zYlLYmDYjGZ1kN9ESI9u2pqy8Ia1K5CYtuX8TLl75sGT3cmUKX/OX5v9UuP5FDND/b7ckIf7Prnn8e3nnH+two1cTbXjVvaLruQ+MS+oBdislghgPEYhroAKA58KMQYq4QYlxUc9e1tb8b3Ogsza3Dnde9yDu3TaNuRa2Vk5yYzOfXe/pVPvewp6Gw5BlIPgO89x5UquRxrHU1zRTU1Dhvu7ZWMSf7WVxUQgpe+9G3Rf7rh/Dpp9rvUiUsHLzXqaN9+2ltTap+r29gmvmsqM6b3WsAzGhTXTP9nCc9R5oLzDUNG8CUKdDB3dpuuwMmzK0Q0HOiR+V0wQXIhQv9xvc5PyGBu5dA6ml3Pu23wy3nWgy+mFDjiO5Hx1B5mSkQgDkfaTZ5DxkKyhBcRfvO5b6VldEEtmTgEibmmCiXSpXgJreJqH3N9iQnWleMux7excZ/+tiWy1thtNhtf9aZGcMc9Isl/DR6ABg40PKQvwWTTlLjiOf/rC/hYjJodFYjjj1+jBWPbdGetSCnitslFmMAV0spW0kpu+if2yKdZ7VUrTJqvx1b3c6qZasih0tuuvxRUpJTyGyqjX493fHpAmXgYt6iRqy8246BEG5ucTPrBq+j50bfY523wMgOw3mrycO+BzGXu/YhuG6NVzSrzP2Uu3Sb8y2PedNzI/Ta4P7vrUcvrXcp/wz5hz6Yt3RFYgmftRnRmDEdsQadELSvaT4gmXpKs8mvvXctnet0BiiofMslB9fDKleynN/jbaq34ZZJf2oTGbxJDzzV5ZVLXiHvqTyqpFSh3mn/eYWDTDXY302eyaY2trMMlkqlK5kf8GOz/yqnNa39bPRixeYHNrNsUOBNmXOH5bLjFd/w9tthfomBlCxRktJJpf0qayco/gvBgHpp9di+/kptICsEu2P/lv0pm1yWm5rf5NP9LH+6BM2qNPM9ycKE0rBSQ9NwATzW/iEqNbaYECUlfdbC1D7Z5sdD4IUWQ7i9eT9TM1Sfxv5bgf56UDXLmS9ECpYrD2ibYFi9iFVSPDfJqJlQwdPm7OEIMEwV8Pjj7sF88DIvac9UcmIy24b4rgxrdFYjKpTSxgraVG/DS91fYlKP9wqOzx+nzZLyh9nCIp97UKIElAptq6/EhERKJPhWiALhYaN24dcM08fw7HivcTGY/kI1AdnFdV9m3zqbpzo+5Tduh1qe4z+1v5zPIze+FXSedSrUseVmIyU5hWoObPQULnGhAACq5+kvhvcYgA190KBSA44MO8LZFd3+fvye9s47tlYgB8vnn0DfxteYHrN8If28ZI83uZMPrplgeszVYvUmc5VXQACfJ4sYyFc9J2pmH3+4uriffopsqynBQPboemn1+POuPwv+l6QE739liCCl5/01U/5vv+1fLhcvvABLl8KmTTB/vmmU2uVrU6Oc10B2X23qqrEs/9fh/zyU18V/Q7cAewh4XAspYcSIyHnQNFyn/OH51LSx6MyDzw1m0jFjPA4FmnbeeC/UPhhkfgZckmccgovTtMHutNJpPhW8N/3O7ccRo6PElBQ466yCvxH1VnrFFZ7/o7iKOG4UAPm6XwfXGECIyfg8wGYP9N13+/cG+NFHBa6Pb2t5m99Kuu5+r4AAL1AkB68+z/ZcBQrARRYvln6d24ua9Gp/C7KSRTcctFlX3+tzaq+7DnpepiXhk6TvXWtetbkNyfXr4jKFdOniPnDPPabxLalbFy6+2CPIr/KdOtU8jnGW1ZeGnWXGmmyLhcnOUo3dJrZYuB0IdbC4Xpr/aaNl8mDrGOvjDfdCSpL1gGiJfPjsU22F9ugmQ1l19ypqlQ8wFqDjPdupRVVtppdxBb8QIjRngP7e265dQ0jQGeJHAbhugEMvS1gV7a23UuuwNlh2daOrDYkKnxerwGVEALldJgbT9Q0O0cdsIWKg9L0efO/KqlV6K3j2WW01aAACVjp+bLpSABkZ8NdfWmveKYzrGWxUigVxjD0nYwvwjjsssvFNO5Z+9SyV3k5rT3Hz+s+j37n9AiduMXjbZA+seyvwPrvXrtEGV5NFCZpW0VxPuEyvvc/pbXqOmfmrceXGHB+6n5t/2ufpvbaQOzQMhvhUAAMGuIOdShe0dNsEXNRsL60g+eGWH3j9W5M57k7aWQ1puWYy3djsRvO4NhTP6adOs8R7ongYiJSy8IPBo1r9+r4v69ln217jECz+7OIuX1HBttZbpbfiPy3+43dnKSecvPmYmIzcd5/1ed7l8TPg3DGjo73y/229Aj0oDOU4u+LZ5A7L5a62d5lG/U+L/5iGl0qt6DEjLuRr7e899DajKhNQBHDdgIQEGDcOsW6983mMGwdLQqjQHLjhdSrU4f7Fkc3DSO1DkP90Preee6v/iH4e/MSERNOpnwNaDaBymcrcsi/AwjAvqqZUhe7d+W/X//LLbb/AyJFQLvz1DHP7zfUNNBkE9sfbl7/N4HaD6dUgOC+mNcrVYFKfSdbTfG3gmuYa8vaEzz0HD1vMTnOS//4XXnstcDy7eD17KckplgooKdG+K3THMTRIo038KACvMYBQCXU5vV28Z9AU5t6m39ac65j+EgZz3epWrMue/9tDnVMmzlj80LWuZksdqGUdGgAAF6tJREFUdvEwLqx1ISQnI10t0hsteio2sBoQB2w/T1XLVuXNy990rqI55xzbA5Mt01sih0saVPJ0c9K5TmfKl9Rcd3jcS+8ylS8PL1svHnOMYcNgyBBAMxf9cIvWm5t2/TSfqBcc0FrN4W7O8+ddf/LWZe7ZPreugJe6v+QTzzVI36+lDROWXZYuhTlzPE2XHTrA7bc7l0cA4kcBjBypeR7UV/mmlda6dVeE2BEIuz3dqBFco8/oMbRUOtTu4PvAX2TwUR+sSSfI+LNvnc3bl9ucGeMPq41VwrhygcwHpgvJCrZBbBJyvoWSFi20qamEPgg8t9/cwD24GNExoyPd63YH4JxKvuNDpf15aXRh49lvXrU5957nXgj50XRtlpY3dQ9oY3YX1b6IuQGm7Jpi5mOsdWvPCQkAP/8MVar4xo0Q8aMAmjTRfPLo8/MrlanEjqE7eDnITTgc8566dq3mVMuI/iJf0/ga9/TGhg21QdJwTTg2z+96dlfuaRfkzJgIUTVPm0Pu7V4gGArzxhz9W/anYqmKfuP4VZgOu+sIixYtgjLfpB8JHMcOz/zk52CE1hmYOSsMyJgx8HXU3J7ZJn4UgAnVUqtRIohdv4xEbSZAyeSIDVpGhTBewpp5pdk62u17J1BlGVW8Bhh71Oth6RLCivFXj2f/o97zfIMRIfyHMJg0BrZ2u07wOe+PPwrMN4H4vf9C/ngvcDw7tN8GG+8zWVpf2ChVCi533hFguMS1AgBgwQJYsyZwPB2XHfVJ11qgKO+nEDSuBVYuW/j8+VAxChWp9xhAiNep9iHd9w4wqE1o+0FEmhIJJZj1n1lcUCuym9uB505RLhyfBWThIyfryiwfJ3mB8pZSUrt8bUZ1H1UQ1rZ6W6oc9XNSEEgReG1BUeDFbi+ydNDSqOerFMBFF3ksqglEuZLlkMOle058tBTA4MHad7A+x2vVgvHj4YsvtP8XXwyVg/C6GSpWYwDBmGS84iYmJPLNTd+w6PZFtpO4o5U2r756qrnL5KLEskHLWHyHe6pXnQp1AHz8UwWD6eD8hReGnJ4ZWx/cyiMdHnE0TVtPUQTfzXnjYdEZ59yYPXrRowXOIqNJLLaEVITCK6/AqFGhmYP693dcHNs4/BJe1uCyoOIPuWAIQy7wNU28c/k7BV5fwyLSPu0N6Xv7mMlslkmF3Ar0bN3T+7Sw8rHDFQ2vYPC3g7mjtfnCtVDTdZQIKoCOWwGCm6ZcGFEKoDBg50EVInDr/2AYTlScxokX3+U+omrV8NPy4u52d4eXwBtvaDvIXRacQgqWRpWslZQQgrMrnh2TSjajQobpnglRo2B6ceA4CmuUCShcQvS+aIrhRS5ZQpsBUzUliMqvfHntE4hRo7QZJK4Na0LBhusGIKR1AAW4vErOmOE/XiyoVQsmT3b2/uvctErw9E8w56bveb7r846nbyTae4I7RaGY21VEr50R1QMIlYcegldfjZgjp4aVGjLuqnFcec6VzifeuzccDtbFoxc5Zo6BDHgNAhcEB/PqXnUV5OYG3A1pROcRXNHwCvB2VFdEmTw9EU6fhjqdIUorVJ0YSDbDVMG4no3q9sZlXD7xjW5OKuZpVVdJfwvBgq2gk5I8N6m3Ytky2LYNfvtN+x9ED+ynfj8VqvEo1QMIFTst7TAZ0GoAZ5U5K3DECHH/efczr/88APo27hvcyU6ZJWxshfd0p6djMoBWHMiokAFoK5WjRokSWu9pwQJb0RtUasC7vd7ls8/cYe+vzOC1WXDxVj8nBqsA/vkncMMGtP0zrjQ0zILIp1OdTj4rsmOJUgChcu212vf118dWjgjy+mWv0zGjI3K4ZOr1U2MtTtR57dLXuLPNnbEWI6I8dMFDTL9huo+Cb5AW4UrqppuCMkHe1fYuj6mjFfNKMGSR2xR0d1t9TKeHYQ/vYBVA1ar2TZsQVadtkUKZgEKlceNiYQOMOGGuA4glZrOHokIUr1ViQqKnS3KdP+76g1Nnwt8YN2ID1F7X6J1e7/BOL33f5O7dYfbsYlFBRxqlAAoDCxbApEnarJLCTmIinLHhgctqDEC9lOY88YR71lMhoHRSaUonlQ47nZgo/gkTYPRozbGawi9KARQGWrf23G+2MPPHHzBvXuB4XgrA5eX0vOrnRUqyos3zhtk+SkkGxp9iqVlTm6ChCIhSAIrgaNpU+wRJ86rNWXX3KmcWXxV2li2DIw55O1P457nnYi1BkUYpAEXUcG3PV+xp1SpwHEV43HSTpmgt/BYp7KFmASkig8UYgELhCEOHwvHjEVklbpti8GyrHkBRZ8gQz6lvhQWlABROM3astqczaM9XBFZhh0QRHrNRCqCo4+Qeqk5ShF8KRSHlDt3xnJ3FWgpbKAWgiAw33KBtw3mHf2+RChOKQe/p9la38+HyD6OW35QpUzgSyYF3fftNsrLcYampWnipUp7hESAvL4/58+cHjJeamsqNQex/rRSAIjLUrg37Q9/tKq5ZtAg++0zzTVNESUqIruxHjhxh0KAIbhi0ZIn23batO2z7dti5U/NpZNOvUagcP36c0qUDr8vIClIRKQWgUBQ2itK6EAsqp2ibDpUvFXmfWVGhefNYSxARlAIoqsybB1OmxFoKhcKUJzs+Sa1ytchslhlrUZyhZMlYSxARlAIoqnTsqH0UikJIcmIyA9vEcI6+3UkIRXicxQnUOgCFQqGwYN68eT6DqsOGDeOjjz5yx/n1VzoE4XdoxowZbNq0yXb8IUOGsG3bNtvxg0H1ABQKhSJELh44kC27d1MzI8P2OdOnT6dUqVLUq1fPVvzRo0dz/PjxUEX0i1IACoVCESILxo5la34+mQ88YHp80KBBrFu3juPHj3PPPfcghGDWrFksW7aMsmXL8uSTT9KlSxd+/vlnPvvsMwYMGMCePXs4ePAgDz74IHfccQddunRhzJgxTJ8+nZ07d7J9+3Y2b97MCy+8wNVX+7ryDgZlAlIoFAo/zJkzhy5duhR8/ve//9k6Lzc3l7lz5zJ9+nTmzZvHsWPH6NevHz179uSll17yMBt9//33rFu3jhtuuIGff/6Z9957j/Hjx/ukeerUKb788ksmTZrEG2+8EXbZVA9AoVAo/NC1a1emGGbcDRs2zNZ5ZcuW5c033+SBBx5g//79XHLJJabxrrrqKgBq1KjB5MmTWbJkCbm5uZwx2XejXbt2AKSlpXHs2DGf48GiFIBCoSh+FILZPTt37mTFihVMnDiR/Px86tevT//+/QHPjXJcC7xefvllWrZsyYABA/j+++9Zvny5T5oJCc4abWKmAIQQJYDfgD5SSn9bOysUCkWRo2rVqvz1119ceOGFJCYmcvXVV1O+fHmaNm3K0KFDSU1N9Yjfq1cvhgwZwtSpU2nXrh379u3jt99+i6yQUsqYfIBHgY1Ahp3477//vgyVtWvXhnxuUUWVOT5QZdYIp34ImW3bpPz9dym3b494VseOHbMVz+w6PPPMM0ukRb0akx6AEKIBUBdYECDeIGAQQP/+/ekY4sKnvXv3khNnHgRVmeMDVWaNvLy8iE2VtCQpCdLStO8I53369Glb5cvLywvqeYi6AhDaruAjgbuBl/zFlVJmAVkAWVlZslGj0LYTzMnJIdRziyqqzPGBKrPG/PnzbTlLc5T9+7VPqVIQ4bztOoNLSkoK6nmIqAIQQgwEbvUK/hj4Skr5r1A+4xUKhSJmRFQBSCnHAmONYUKIKUA1IUQ/oBHQSAhxuZTyQCRlUSgUCoUnUTcBSSkLHGsIIcYDz6jKX6FQKKJPTNcBSCkHxDJ/hUJRPBEj7JmX5fDYrxeIJcoVhEKhUFhg5g3UxYGDB2l3662s27jRVlpbt25l2rRptvP+7rvvGDt2bOCIYaAUgEKhUATJK6+8Qvsrr2TVX3/ZPmfLli1MnTrVdvwePXowcGBk91RQriAUCoUiSB5++GEevvFGulh445wxYwYvvfQSCQkJVKtWjTfffJMHH3yQHTt28NRTT1G/fn0WLFjAzp07uffee9mwYQOTJ08GtBXEM2bMYNKkSeTk5HDXXXeRmZlJw4YN2bBhAw0aNPDYjyAcVA9AoVAoHGbKlCk8/PDDLFiwgN69e5Ofn8+YMWPo2rUrzz33HACLFy/m008/5bLLLmPHjh0sWrSIRYsWsX79ep8NY9avX8+YMWNYuHAhK1as4J9//nFETtUDUCgUCod55ZVXGD16NBMmTCAtLc3UE2i3bt1ISUkBoEqVKjz88MMIIcjLy/PxBFq3bl0qVqwIOOcJFJQCUCgUxZBYz+754IMPeOKJJ0hLS2P48OFMmTKFc88919QL6IoVK5g5cybz5s3jxIkTTJ8+3Sc9p72AulAKQKFQKPzg2hDGxeDBg+nbt69743kTjwYtW7akR48epKSkkJiYyIQJEzhz5gy//vorTz75JA0aNCiIW79+faSUdOnShRo1anD++eczbtw4mjZtGvGyKQWgUCgUFnTq1Indu3ebH6xalbnTp0P16j6HevfuTe/evX3Ct2719XxftmxZ5s+f71eOn376qeD33Llz/QsdBEoBKBQKRSgkJkLNmrGWIizULCCFQqGIU5QCUCgUijhFKQCFQqGIU9QYgEKhKPKkpqaSlZUVazEiRl5eHklJSQHjee8zHAilABQKRZHHymFbcSFSO78pE5BCoVDEKUoBKBQKRZyiFIBCoVDEKUVmDGDnzp17R4wY4buMzga//PLLWR06dNjrtEyFGVXm+ECVOT4Is8wZVgeE0TlRcUUIsURK2TbWckQTVeb4QJU5PohUmZUJSKFQKOIUpQAUCoUiTokXBVB8V4hYo8ocH6gyxwcRKXNcjAEoFAqFwpd46QEoFAqFwgulABQKhSJOUQpAoVAo4hSlABQKhSJOKbYKQAjxphBikRBiiRCil8nxRCHEeCHEQv3TxSydokSgMutxbhNCLBNCrBBCPBFtGZ3GTpn1eCX0cluuiiwK2HiuU4QQ2UKI+UKIX4UQF8RCTiewUVZb974o4a9MEbm3Uspi9wF6AF/rv6sDW4AErziXAZ/pv+sDv8da7iiUuSHwDZriF8BT3nGK0sdOmQ1xHwU2AhmxljvC9/h+4Fn9d3Pgt1jLHYmyBnPvi8rHRpkdv7fFtQfQGvgBQEq5A9gD1PaKcwAoL4RIBKoARd23iJ0y9wI2AzOB+cBKKWV+NIV0GDtlRgjRAKgLLIiqdM5jp7xrgHH678NAYtSkc5ZAZbV174sYgcrk+L0tMs7ggiQZ2G/4fxQo6RVnOZAE5ADpwM3RES1i2ClzNbSHrBtQFlgohPhZSllUlV/AMgshBDASuBt4KXqiRYSA5ZVS/ggghKgFTAKejpp0zhKorHae96KG3zJF4t4WeQUghBgI3OoV/CNQyvC/NHDSK85jwEIp5RP6BV0ohJj9/+2df8xWZRnHP1/IkjfKlCGizh8pCeK0wOmwuVVqgrFESZkDArN31ciwzAxpueaPpluWopk/EnFTfiymKGobNTICF2LTQBEIRRzpYmX+QE3Ub3/cF/Dw+Ly/5H14eH2uz3a2+5z7Pve5rvu8731d5zrPuW7bW+onbfewCzq/CTxk+w3gDUmPA0PoAZ7xLujcCiy0vbnYgp7BLuiLpEnAZOAHtlfUS8Y6s5X2de2ovifSoU7dfm8bHfeqUyxtJHB/lAcCzwK9q9pcCVwU5b7A80BLo2Wvs86nAL+nxP/7AE8BAxote511ng38CVgMvAg8CuzbaNnrqO8Y4K7q4z1t60jXzoxFT9s6oXO339uGK13HwbwR+AvwV+ArcWwAsDjK/YCHKd7vY8DERstcb51j/yrgb6F3U+hc0XYmPfglcGf0BRZQwpuLY7uj0TJ3l641/pbfNxY9fWtP53rc28wFlCRJ0qR8WH8FlCRJknRAGoAkSZImJQ1AkiRJk5IGIEmSpElJA5AkSdKkpAFIkiRpUtIAJG0iaaqkgVEeI+nmOl7rJEk11z2V1EfSHyT16WKfV0s6v3sk7PBan5R0h6QZNeoOlfS4pLEVx3pL+rmkNZLWR2baPhX1P5H0j6i7WVKvOD5M0lJJayU9KWlkxTkPSHpW0urYrojjh0haENdaI+l7NWRsif4mVRw7MrJOronssSM6I3sb4/MFSe90bVSTutPoDx9y23M3SuK44bvhOr2AJ4AD2mlzATC90WPShmy3URJ3bQFmVNWtAjYDbwNjK463UhLy7U3JAXM/OzI9tgIrgJaoWwqcS0ndsgE4I9qdBLwO9I39ml+zU76EvjjK/ShJxU6vanMn5cvwSbG/7VrjYn808HRHsrcxPgdRkpy50fcqt523fAJIaiJpOeUfd76kr0uaJOmBqFss6aeSlkjaKOlaSZdHvv1Nki6o6Oei8B43SJovad8alxsJPG/7JUlDVNZnWCNpnaRLos0c4DuS3pe/StK3K7zbFZKGxfGZki6OskOWp0LmqyvObw2P+ZnQ6dg43j/yr68Kz/rCWmNlu9X2/tRINmf7GNv9KZN4JYcB19t+y/bbwB+BI6Lum8AVtt+IurMpX34OANbbfij6XkYxLAMjn5WA++LpYJak/tHf8cDv4px/U76AP7NC/+9SkiIur5DvNOBl23Nj/0HgnE7IvhOS9qJksJxSqz5pLGkAkprYPgHYRPFa76rRZDjwReDzwCXAf2wPo2TdnA4g6VzgOOB424dR8vDcWKOvkcCyKP+I8un7UcCJwGhJ+7gk6dsY19tOTDA3AKPjnGtoO7PrfraHUibEqZIOknRyyHuq7cEh38xoPwu43fYxwGeB8ySd1kbfXcL2dNvzQ4d9gPOBhVF9FHCkysIfzwDTgNdsb7J9SoXuk4GXKDljDqZ45ROBYwEDt0TT5UCrpF6SDqakGDgg+hhByRF1TZWIRwHPhSFZFX3364Ts1fwS+K3ttV0comQ30OOzgSYN427b7wAvSNpC8RAB1gKfiPJYyoS9UiUTZ2+g1voDR7DD+1wA/CY82iXAebZfiboXgEHAI9tOtL1V0jxgoaQFcc6P25D51jjnX5I2UybBrwKzbW+KurnAXEkfp3jBh2tHFtEWSthlUQdj02lUVn36BXCr7dkV1xlBWSDkIxTv/YeUBIZI2o+S4noYMMr2uxTj+mhFv5cCG1TWu5gM/Iri5W8A7gUOlLR/XHu0bWvnbKktwBnAabYnSRoF3CvpUNuvtyN7pW7jKWGfebs0SEndSAOQfFDerSi7Yr8yuVRv4Ge2bwOQ9DFqL9qx3SjYvk/SIsoEeCrwpKRRLqlvTY2/WdsTwmCcTIlPX07x8qt5s0rmbTJuPx4T5nBKnFzA52y/FXUDqvT+wEQo69fA4cCZttdUVP8TmGn7zWh7DxF+kXQiJRf8DGBKGGHiyWRLhIUA3qGMVS9KTvkJLmnAkXQ9xRB8ifLksDQm/4HA6ZL6hQwrbS8BsP2wpFeBwZKeaEf2Ss4BhkpaXaH3aso7jOe6PmpJd5MhoKQ9zK6tOvQgJW7/qdj/PiWcUc16yvsGJN0NjHdZ/GIasBIYHO0OjLbbkbSvpI3AVtv3ABcCR0doqDMsAsaHNwwlfHRdeLlLKSEpJLVQltMc2sl+O+IiShryL9eYQOcDZ0fIZi/KU8pjkvoC8yhhuRnbJv+gHzCjaqwX2t5KCXFdGnocAowD5tieY/sQ20NsD6E8GUyzfR3l3n1a0tFx3gmUJ7t1Hci+HdtjbA+q6J8o5+S/h5BPAEl7LKa8VGwrpNIRs4BDgeWS3gNWA9+o0e4BYGqUrwRuj5e3oqTGnR9PD4OoWrzG9suSrgL+HD8zfA/4VoSGOhQwPNubgCWSTFkzYHJUTwRukrSO4lHfYvuR2j11mXEUj/vpCjnvtX0ZcBnlPcQG4DXKr3iupXjsA4B5VbqdYnuOpM8AKyT9D/g75WUylJDYnaHHRym/2FlNO7gsoDOeEvYB+C/wNduvSGpTdklTgLNsn/oBxiTZzWQ66GSPQNIyimf7Yhv1E4DjbF9Sqz5Jkq6TBiDZI5B0PMVzb61RtzclJHGW7Vd3u3BJ8iElDUCSJEmTki+BkyRJmpQ0AEmSJE1KGoAkSZImJQ1AkiRJk5IGIEmSpEn5P+sZX1lFTrmNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "deltat = 0.5\n",
    "shiftt = 0.3\n",
    "\n",
    "indxt = np.where((time >= tevent-deltat-shiftt) & (time < tevent+deltat-shiftt))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(time[indxt]-tevent,strain_H1_whiten[indxt],'r',label='H1 strain')\n",
    "plt.plot(time[indxt]-tevent,strain_L1_whiten[indxt],'g',label='L1 strain')\n",
    "try:\n",
    "    plt.plot(time[indxt]-tevent,strain_V1_whiten[indxt],'b',label='V1 strain')\n",
    "except NameError:\n",
    "    pass\n",
    "plt.xlabel('time (s) since '+str(tevent))\n",
    "plt.ylabel('strain')\n",
    "plt.legend(loc='lower right')\n",
    "plt.title('Advanced LIGO strain data near '+event_name)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save GW event\n",
    "np.savetxt('./{}_H1_{}Hz_{}sec_whiten.txt'.format(event_name, sample_rates, indxt[0].size/sample_rates), np.asarray([time[indxt],strain_H1_whiten[indxt]]).T)\n",
    "np.savetxt('./{}_L1_{}Hz_{}sec_whiten.txt'.format(event_name, sample_rates, indxt[0].size/sample_rates), np.asarray([time[indxt],strain_L1_whiten[indxt]]).T)\n",
    "try:\n",
    "    np.savetxt('./{}_V1_{}Hz_{}sec_whiten.txt'.format(event_name, sample_rates, indxt[0].size/sample_rates), np.asarray([time[indxt],strain_V1_whiten[indxt]]).T)\n",
    "except NameError:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "# Save noises\n",
    "deltat = 0.5\n",
    "shiftt = 0.3\n",
    "for i in range(1,4+1):\n",
    "    indxt = np.where((time >= tevent-deltat-(i+0)-shiftt) & (time < tevent+deltat-(i+0)-shiftt))\n",
    "\n",
    "    # plt.figure()\n",
    "    # plt.plot(time[indxt]-tevent,strain_H1_whiten[indxt],'r',label='H1 strain')\n",
    "    # plt.plot(time[indxt]-tevent,strain_L1_whiten[indxt],'g',label='L1 strain')\n",
    "    # plt.xlabel('time (s) since '+str(tevent))\n",
    "    # plt.ylabel('strain')\n",
    "    # plt.legend(loc='lower right')\n",
    "    # plt.title('Advanced LIGO strain data near '+eventname)\n",
    "    print(i)\n",
    "    np.savetxt('./{}noise_{}_H1_{}Hz_{}sec_whiten.txt'.format(event_name, i, sample_rates, indxt[0].size/sample_rates), np.asarray([time[indxt],strain_H1_whiten[indxt]]).T)\n",
    "    np.savetxt('./{}noise_{}_L1_{}Hz_{}sec_whiten.txt'.format(event_name, i, sample_rates, indxt[0].size/sample_rates), np.asarray([time[indxt],strain_L1_whiten[indxt]]).T)\n",
    "    try:\n",
    "        np.savetxt('./{}noise_{}_V1_{}Hz_{}sec_whiten.txt'.format(event_name, i, sample_rates, indxt[0].size/sample_rates), np.asarray([time[indxt],strain_V1_whiten[indxt]]).T)\n",
    "    except NameError:\n",
    "        pass        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}