Examples/gmimtools_Example_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a4d57f1c",
   "metadata": {},
   "source": [
    "Example: GMIM tools\n",
    "=================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7eb3bf59",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import pyrotd\n",
    "import gmimtools"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12b4c7df",
   "metadata": {},
   "source": [
    "# Input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ec4197b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = './Data/1989_LomaPrieta_TreasureIsland/'\n",
    "path = './Data/1994_Northridge_RRS/'\n",
    "path = './Data/2023_Pazarcik_TK_0125/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b36c9271",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaling = 9.81 # To convert acceleration units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a5072fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = os.listdir(path)\n",
    "file  = files[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50a5c4c4",
   "metadata": {},
   "source": [
    "# Read acceleration time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b37797b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "[accel,dt] = gmimtools.read_AT2(path+file,scaling)\n",
    "gmtime     = np.linspace(dt,dt*len(accel),len(accel))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8c3ad538",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArEAAADZCAYAAAA34mgzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABn40lEQVR4nO2dd3gU1dfHv7ubZNNJQkIivUPoTelNehEQUESkCIgUqYKgIL0jRaqKgoUqP0BBEUFAEAVpiYCEHggQIEBCOim78/6Rd8Y7szOzs5vdZEPO53l4mMxOOTt3du53zj33HB3HcRwIgiAIgiAIogChz28DCIIgCIIgCMJWSMQSBEEQBEEQBQ4SsQRBEARBEESBg0QsQRAEQRAEUeAgEUsQBEEQBEEUOEjEEgRBEARBEAUOErEEQRAEQRBEgYNELEEQBEEQBFHgcMtvA/ISs9mM2NhY+Pn5QafT5bc5BEEQBEEQhASO45CcnIzixYtDr1f2txYqERsbG4tSpUrltxkEQRAEQRCEFe7cuYOSJUsqfl6oRKyfnx8A4NatWwgMDMxnawo3JpMJV65cQZUqVWAwGPLbnEINtYXrQG3hWlB7uA7UFq5DXrRFUlISSpUqJeg2JQqViOVDCPz9/eHv75/P1hRuTCYTfH194e/vTw+kfIbawnWgtnAtqD1cB2oL1yEv28Ja6CdN7CIIgiAIgiAKHCRiCYIgCIIgiAIHiViCIAjCKk+ePEFERER+m0EQBCFAIpYgCIKwSr169VCvXj2cOXMmv00hCIIAQCKWIAiC0EBMTAwA4JdffslnSwiCIHIoUCL22LFjeOWVV1C8eHHodDr88MMP+W0SQRDEc09SUpKwHBwcnI+WEARB/EeBErGpqamoXbs2Vq9end+mEARBFBpSUlKEZbXqOQRBEHlJgcoT26lTJ3Tq1Cm/zSAIgihUPHv2TFjOzs7OR0sIgiD+wy4Re+vWLfzxxx+4desW0tLSEBISgrp166Jx48bw9PR0tI0EQRBEPpKRkSEsWxOxJpMJACghPUEQTscmEbtlyxasXLkSp06dQrFixVCiRAl4eXkhPj4eN27cgKenJ/r164fJkyejTJkyzrJZMxkZGaKHLx/XZTKZhActkT/w15/aIf+htnAdXLUtUlNTheWsrCxV+4YNG4bt27fj0qVLKFGiRF6Y5zRctT0KI9QWrkNetIXWY2sWsfXq1YNer8egQYPw/fffo3Tp0qLPMzIycOLECWzbtg0NGjTA2rVr8dprr9lmtYNZsGABZs2aZbH++vXr8PX1zQeLCClXr17NbxOI/4fawnVwtba4cuWKsBwbG4uoqCjZ7R4/fowNGzYAAJYtW4ahQ4fmiX3OxtXaozBDbeE6OLMt2Dh8NTSL2Dlz5qBLly6KnxuNRrRq1QqtWrXC3LlzER0drfXQTuPDDz/EhAkThL+TkpJQqlQpVKxYEYGBgfloGWEymXD16lVUrlyZhh3zGWoL18FV2+Lhw4fCclBQEMLDwy224TgO7u7uwt8Gg0F2u4KEq7ZHYYTawnXIi7ZgM6KooVnEqglYKcHBwS6RhsVoNMJoNFqsNxgM9CNwEagtXAdqC9fB1doiKytLWOY4TtY2s9ks+luv17vUd8gNrtYehRlqC9fBmW2h9bgFKjtBSkoKrl+/LvwdHR2NyMhIBAUFWYQ3EARBEI5By8QuORFLEAThTGx+yqxduxZt27bF66+/jsOHD4s+e/z4McqXL+8w46ScOXMGdevWRd26dQEAEyZMQN26dTF9+nSnnZMgCKKwoyXFFsdxor/JW0YQhLOxScSuXLkSkyZNQtWqVWE0GtG5c2csWLBA+NxkMuH27dsON5KnVatW4DjO4t/XX3/ttHMSBEEUdlhPrNKsYfLEEgSR19gUTvD5559j/fr1ePPNNwEAI0eORI8ePZCeno7Zs2c7xUCCIAgif9HiiSURSxBEXmOTiI2OjkaTJk2Evxs3bozDhw+jTZs2yMrKwrhx4xxtH0EQBJHPaPHESsMJSMQSBOFsbBKxwcHBuHPnDsqWLSusq169Og4fPoyXX34Z9+7dc7R9BEEQRD5jjyeWktITBOFsbHpVbtasGXbu3Gmxvlq1ajh06BD279/vMMMIgiCI/MFkMuHMmTNCai17RCyblosgCMIZ2CRip0yZgtq1a8t+Vr16dRw5coQyBRAEQRRwFi1ahBdffBFTp04FYF84AYlYgiCcjU3hBLVq1UKtWrUUP69evTqqV6+ea6MIgiCI/IMXr0uWLMHixYtFnlipx1VpfWZmpvMMJAiCgB15YgFg//79OH78uPD3mjVrUKdOHbz55ptISEhwmHEEQRBE/sN6YrWKWFs9sYmJiYiKirLdOIIgCi12idhJkyYJdW0vXLiA999/H507d8bNmzcxYcIEhxpIEARB5C+sJ1YaNqC03lYR26xZM1SrVg0XLlyw3UCCIAoldpWdjY6ORrVq1QAAO3fuRNeuXTF//nycO3cOnTt3dqiBBEEQRP5iTziBrSL24sWLAIAffvgBNWvWtNFCgiAKI3Z5Yj08PJCWlgYA+O2339C+fXsAQFBQkOChJQiCIJ4P7Akn4GNio6OjsWjRIs19g06ns9NKgiAKG3Z5Yps1a4YJEyagadOmOHXqFLZv3w4AuHr1KkqWLOlQAwmCIIj8RYsnVimcoFmzZoiNjUVUVJRiiXA244FSCi+CIAgpdnliV69eDTc3N/zvf//DunXrUKJECQDAL7/8go4dOzrUQIIgCCJ/yU04QWxsLABg7969isdnhSuJWIIgtGKTJ/bAgQNo3bo1SpcujZ9++sni8+XLlzvMMIIgCMI1SE5OFpaVJnZZi4lVS7nFbqsllpbjOAo7IAjCNk/s8OHDERISgj59+mDr1q1ITEx0ll0EQRCEnezZsweNGzfGtWvXcnUcd3d3AEBKSoqwTipWt27dihIlSuDkyZOi9VLRyh5Dii0idsaMGShXrhwePnyobjxBEM89NonYmzdv4tixY6hZsyaWL1+O0NBQtGnTBitXrsStW7ecZCJBEARhC927d8fJkycxcOBAm/dlRaqbW85gXWpqquznAPDmm28iNjYWr7/+umi9nBh98uSJ7DltEbGzZ8/G7du3sXjxYtXtCIJ4/rE5JrZWrVqYNm0aTp06hZs3b+K1117D/v37ER4ejtq1a2P69Ok4c+aMM2wlCIIgbODBgwc278N6UHkRywpLpZhYKXJi9Ny5c7LbsnGwbPytVjsJgiic2DWxi6d48eIYPnw49u3bh8ePH+Pjjz/GrVu30LFjR8yfP99RNhIEQRB2YDAYbN6HTafF789mD8iNiPXy8rK6rVYRy9qxYcMGTJs2TTFelyCI55NciVgWHx8f9O7dG99++y3i4uLwzjvvOOrQBEEQhB3o9bY/4lkRy+/PilitQlHOU8rH2EpRErEnTpxAyZIlhTSOLKxNQ4YMwbx583DkyBFNthEE8XxgV55YADh16hR+//13xMXFid6IdTodli5dipCQEIcYSBAEQdhHbj2x/LOdHe7PjSeWFZ5K27IitkOHDkhOTsYbb7yBPn36iPaRs+P27duabJOye/dubN++HevXr4efn59dxyAIIu+xS8TOnz8f06ZNQ5UqVRAaGipKdUJpTwiCIFyD3HpiedHpqHACpRyw7Lbp6enCMpvaSwrvEWbtsTdOtmfPngCAcuXKYcGCBXYdgyCIvMcuEfvpp59iw4YNGDRokIPNIQiCIBxFbkUsLzpzI2JDQ0OFdFhKIlZuYpdWQcpuZ4+IvXfvnrBMabsIomBhV0ysXq9H06ZNHW0LQRAE4UByG07Ai0t7wgl4QckKSy2eWF7Efvzxxzbbm1sRW7RoUZv3Jwgi/7BLxI4fPx5r1qxxtC0EQRCEA3EFT6zc8ZS2Bf4TsVorQLLClT2XVtjvVrJkSZv3t4XZs2ejbt26iI+Pd+p5CKKwYFc4wcSJE9GlSxdUqFAB1apVs5hxumvXLocYRxAEQdiPNRF78OBBXLlyBSNGjBC8tqwQ5DgOZrNZJD61ZieQE7FKAlhOxBoMBk0laNnjs0UZtMJOJLNHBNvCjBkzAOQI9Dlz5jj1XARRGLBLxI4ePRpHjhxB69atUbRoUZrMRRAE4YKohRNcvXoV7du3BwBUrlxZWJYKuezsbJFwtcUTazKZRJ5OW7ITWAuF4G1i7VWbCKZEbsMR7IGdvEYQhP3YJWK//fZb7Ny5E126dHG0PQRBEISDUPPEjh07Vlg+f/68ooiV/m1NxOp0OnAch8zMTAtRaIuI1eoccaSIdbYntiCRmpqKDh06oE2bNpg1a1Z+m0MQstgVExsUFIQKFSo42haCIIjngq1bt6JixYr4559/8vzcrNdUTcTu379fWGbFX25FrNFoBJAjWKXVt5T2ZcVuWloaAOsiVs4Tm5KSorqPHOz+WsIXHEFBqCy2f/9+/Pnnn5g9e3aBsJconNglYmfOnIkZM2YIDxuCIAjiP958803cuHHDIkF/XnDr1i1h2cPDQ9M+d+/eFZZzK2LZc0pjVLV4YtPT03Hjxg3Nk9Jy64llhba1cIJffvkFpUqVwr59+2w+T0GDbZMnT57koyUEoYxd4QQrV67EjRs3EBoairJly1pM7Dp37pxDjCMIgijIxMXF5fk5N2zYICyridgmTZrgr7/+AgA8evRIWG9NxFrzyvGeWMBSVCqJWKl4rFixIgIDAy22Y8+dHzGxc+fOxd27d9GlS5fn3jvJxu3eu3cPwcHB+WgNQchjl4jt0aOHg80gCIJ4/siriUIsL730krCslNIKAKpUqSKIWNbr5qhwAgBITEzUtK/cdZLzxMqJYHZfZ4tY/nrxZGVl4eLFi6hdu7ZN6cwKggBm286erA8EkRfYJWL5NCEEQRCEMvk9UUga42k2m3H69GnUq1dP9JkjRSw7MmdPOAGPXEysnCjXGhN79+5dGAwGvPDCC4r72/rSMXDgQGzduhXLly/HuHHjVLdlv3tBELFJSUnCsjS2mSBcBbtiYrVQEH6kBEHkjh07dqBy5cqIiIiwab87d+6gUaNG+Oqrr5xkmWug5gl1BGlpaTh8+LCiIJWKw3HjxqFRo0ZYv36900Ssm9t/vhGpqJTuy/cTcuJRTvCy6/jYXy3hBOnp6ShVqhSKFy9ucdzciNitW7cCyAkzYLl8+TI6d+4s8txu2rRJWLZ2Dfft24dWrVrh5s2bNtnjSFhnFaUEI1wVzSI2PDwcW7Zssfojv3btGkaMGIFFixbl2jiCIFyb119/HdeuXUOvXr2Edffv38dnn32mOgQ5f/58/P333xg6dGhemPncMnToULRp0wbz588X1rHCWfq8XrVqFQBg1KhRIuHK7mOLiJVzVuj1eiHHq1TE8gLSZDKhYcOGaNCgAbKzs2X7FTYsgRfGrJ1HjhyxsE9JbD18+FBYltrEXgd7wz+koQSvv/46fvnlF1F59rVr1wrL1kRsly5dcPToUQwcONAue3KL9OWLPLHyZGRkoEWLFnj99dfz25RCi+ZwgjVr1mDy5MkYNWoU2rdvjwYNGqB48eLw9PREQkICLl26hOPHj+PSpUt47733MHLkSGfaTRCECxEbGyssd+3aFefOncOlS5ewcuVK2e2vX7+eV6Y91/CewPnz5wueMzVPLIu9nlhWuMp5mnkRazKZFMMJ7t27h1OnTgHIuXf483t4eAhCkg1L4EML5ESmFk8q631NTk5GkSJFZL+Do0RsTEyMqp1aC0bcuHHDLntyC5utAiiYntiIiAjMmjULc+fORY0aNZxyjpMnT+KPP/4AAMTHxyMoKMgp5yGU0eyJffnll3H69Gn8/PPPCAsLw5YtW/Dee++hX79+mDlzJq5du4YBAwbg7t27WLhwIfz9/Z1pN0EQDoTjODx+/Nju/dnOn89OsnnzZsXtrVVjIuRRCtPy9PQUlrWKMkeEE8iJZNYTKxWx/L7sS096erpgJ+vRZ1M4ZmVlgeM47Ny50+J8bAhBVlaWcI4rV67g119/Fc7BI/UqOsITK72f2faQO69WEasUQ5xbDh06hO3btyt+znquAfuuS1xcXL5OCGvfvj1+/PFH1KxZU9WTnJmZiZUrV+Ly5cvCuqVLl6JUqVJWX7b5FzEAiI6Ozr3RhM3YHBPbpEkTfPrpp4iIiEBCQgKePXuGu3fvYu/evXjvvfcQEBDgBDMJgnAmy5cvxwsvvCB0+rYiJ674dYsWLYK/v78o8b9S5/zkyRPMmzcP9+7ds+n8MTExeP/993Hx4kWb9itILFq0CHq9XvblgPda7tmzB7t27RLWaxWx9oYTKE3IUgsnSEtLw/r164V1qampwnG8vb2F9dI85CaTySI7wD///COqPMbaVLVqVXTs2BEnT54UHcsZIlbqifXy8rLYhj2v1gl/WsQux3E2F2lo27Yt3njjDURGRsp+npCQIPrb1gmKCQkJKFu2LGrVqmXTfo6EfSlXKzqyZs0ajB07FuHh4cK6iRMn4u7du5gyZYrifmlpafjggw+Ev6XeayJvcNrELoIgCg4bN24Ex3H4/PPPHXZM3vs1ZcoUJCcnY8SIEcJnSjO116xZg2nTpqFRo0Y2neujjz7CsmXLMGTIEJv2S05Oxvjx43HixAmb9rMFtUmuWj1y2dnZQof61ltvyX6enp6O7t27Y8+ePcJ6Z4QT5MYTazKZsHDhQlEuW1bEsuJPOoSdlZWF0qVLC3+XK1cOn3zyicX5pfZGRESIjpWRkQGz2SwIbGeIWDlPLPvdNm7cqOm41iZIZ2ZmQq/Xw8PDA3fu3NF0TJZ///1Xdr1UxNp6XU6ePIn09HTcvHlT8z3On3f79u2KEyLtnTCulv7s2LFjip+phVGwXljA0ntN5A0kYgmikMN2DPwwb3R0NFq0aIEffvhB83HYhPmApcfrxIkTQtoeJSH0zTffALDdq/HLL78AsOxYrLFz506sWLECTZo00bT94cOHUaZMGfz000+az6E0pNqtWzdUr15dU7whWzGJjefkycrKkj0OLz4eP36MTp06CesNBoNImLAvFbaIWKX8rvxELLmys9IXpaysLFkRKyUrKwtVq1YV/s7OzpYVIJmZmaL0UL6+vhZe0H79+iEwMBC3b9+2W8Syvxtr4QSPHj0SDVcDwPTp0/Hjjz9qPocc7OS37777TnVbHlYgKsXcSkXsuHHjRCEg1mCvoy2lgCdNmoQ33ngDkydPlv2sZcuWdsUJq/3GpG3OXnM204YUafjAgwcPNNly6dIlXL16VdO2QI44njBhgsX9Q+RAIpYgCjmsB4EPB5o8eTL++OMPvPrqq5qP8/vvv1usk3YQXbt2BSD2jLAiT86DpQW1zkaNqKgoYVmLgBk6dChiYmLwyiuvKG4jFR5JSUm4efMmOnTogN9++01Yv3fvXly+fBlHjx61el62E65Tp47F56wQlK4HcgTA/v37hfVubm52hxOw38+aJ1Z6TU0mk0X+V9Z2Dw8PxXhp6XfMysqSvV8yMzNF4Sju7u6i/Xbs2IFt27YhOzsba9euFX1mS7EE9rtZ88QeOHDAYv85c+ZYLUss58W8desW5s2bh6dPn4rWa52Hwgp6aWowHqmIBYAJEyZYrDOZTNi5cyfu3LmD/v37o0+fPujbt69IpEkLXqjBp9zjBXlaWhrat2+PadOmYfny5UhISMCaNWs0HcvHx0dYtkXEstdVLXZfGoqhpTRvcnIyGjRogNq1a2ueLDd48GAsX74cHTt21LR9YcO+Jz9BEM8NbGlUXlAqDU1mZ2crCka5uLnbt2+L/uZn8rKkpqYKJUalsbIXL15EWFiY1ZKXAQEBwvdIS0sTxVaqwXr+4uPjERYWJvqc4ziMHj0aISEhmDFjBu7fvy/aXm42srRTTExMxNixY3Hw4EEcOHAAHMeJROOMGTOsdlBqMZ1AjqCTE+H8OqnXyBYRKz0f20bWYmLlvLhSEZudnS0cx93dHW5ubrIx01IRm52dbVHynD8n6/1LT08XeSxXrFghLBsMBtF3v3z5MubPn4+PPvrI4rjSlxP25UstJjY+Pl7Ro5qRkQGO42QLOwDyIrZnz56IiIjA2bNnRet9fX1ljyGFbU+59svOzsa0adMs1st5vTdt2oRBgwZZrN+2bZuwbM/kLv4l4Oeff8bBgwdx8OBBkX3WMJvNovNK46ul2/JkZmaKPKpKQvPZs2dC5pWgoCDEx8cLfy9btkxR/EZHRwvHjI2NRYUKFax+l3379gGwfJYSOZAnliAKOWxHFh8fD0B+4tX3338PPz8/YYa4VDTJebGUvBOsuGE7G7ZDuXbtGmrWrKn6oN+3bx+OHj0qEilSD5Ua7LnlsjPcvHkTa9aswcyZM3Hx4kWRAGC9uCxS0ZeUlGQhItnvryUEgu1MlTzGcuuVYmI9PT1Fn6mFE6gNtzraE+vu7i4rTPljSUWsh4eH7Hbsvfjs2TPF62A0Gi0+mzp1quy20nZlhZFUxLLfoW3btorfCVCPW5YTsXxhEWmoj9bJV+z38PT0tPit8/l3pfAvmixq8aQ8agJy3759GDZsmIVY5K+n3POD9bDKkZmZiWHDhonWqWUnYO/HxMRE0bC9kgBnU6hVrlxZWF65ciVOnjypeC72JXjkyJGq2xLasNsTazabcf36dcTFxVn80Fq0aJFrwwiCyBtYocEPI8p5O/ihzz59+iA7O9si1o2NQ+RR6ljZc/IdxZYtW3Dt2jVh/fHjx4Xjynmr4uPj0aVLFwA5k3xYO4oXLy57XimszXIdJtv5SQu4yA25ApbDp4mJiRaeOGmneu7cOdStW1fRTlYIaLmm7Dqp5xewFG+2hBMkJiaid+/e6N27N6pVq2ZxTlbESvfVImKVPP1y4QRyIlbqic3IyFAUih4eHppn9kuFFtsm0lhe9jtERESoTjhMT0+X/R6Aekys9DOtw9Psvffs2TNcvXpVNDOfFVosfn5+FuvKlClj9XxqIpb//YaFhWH27NnC+tu3b+PBgwey97o1sb57926LSoBq14a9HomJibhw4YLwt5KIZW146623RGI0JiZGVOSChb22Bw4cEEZmCPuxyxN78uRJVKxYEeHh4WjRogVatWol/GvdurWjbSQIwonIeWLVhux4z41UxCYnJ6sOubJIReypU6fQr18/0Tas90quE+JtBSCKgZSK6fj4eFHHJLWZR07EstdGOnFNKdaPL4eqZA9g2RHXr19fmGl/5coVi3NpEbFy6zmOg8lksmhPqSfWlnCCR48eYefOnejbt6/VcAKpsFYKJ2ALHEi9lrznTc4TK+fhlHpi1USsj4+P3SKWvf8jIyMxatQo4dpJPbPs74AdNZA7LostM/u1eGKPHTtmEYfLhhMBynHBcs8EOe+slLS0NBw5cgSrVq1SFGzXrl2zOP7MmTNlJ3gqTRT766+/EBUVJevdVPPEst83MTFRdHyl5xf/eyxbtqzFNVDLEqF18lde8OTJE3Ts2NGhGWnyA7tE7PDhw9GgQQNcvHgR8fHxSEhIEP6xHQtBEK5DbGwswsPDhRKljx8/RkpKikhoZGRkID09XTTEqJTTVfqAT0pKko0HlYPdLi0tTbbjYbeROw57fnZbqWjs1asXateuLQzbHzp0CBMnTkRiYqJVEcsKDGkOXTlxKnccOdvlOtWZM2fi33//Re/evUWz8AGIYmZt8cQCOULckZ5Y6bGlsNkJ5MIJ5I6hFk5gq4jNyMgQ3RtqIlZ6TDWkYlMqSNauXSt465WO2a1bNwub1WKObRGx7HG+//57jB8/3uLYLVu2FIqR8Ej7bPbaNWzYUFiWu9+1XLu0tDS8/PLLGDNmDPbu3Su7jaenp8XvJDo6WjaFWnJyMnbu3IkaNWoIccE//PADmjZtimrVqolinnnUXhRY0ZqYmCj6/tZErLe3t0V4g1qOa7mQJUd7YjmOw4YNG0SjWnKsWrUKv/76K4YPH16gvcF2idhr165h/vz5CA8PR0BAAIoUKSL6RxCE67Fnzx5cvnwZU6dORWxsLMqVK4f27dtbdEQJCQmKMasscuEEUsGjFJ8q9cTKDSGznas1Ecty6tQpDBs2TOhMfv/9d3Ach23btuHx48do27Ytli5dipdfflnUMct1MGpDoUoCXfqdZ82aZdFJyAlD3mMlPbb0PEpiVU3cSsWQXq/XHBOrJmKVUmypeWKlsBO7PDw8RPeCXq8XJvlIBafJZFIMoWDb1VkiVk6sXLp0CYB6LLLUSyu9vmzVMns8sb/++iv69OmDFStW4I033rC6n/R3xP89YsQI0f5yHlotGT3Y3xArrNgXJw8PD4vfDZvJgyU5ORlTp07Fv//+iwEDBgCA1Swqap5YNRGr5PXlxbiciFVLDyj3zFB6GQYsBW5GRgbi4uIsin6wbN26FUOGDBHF6srBFoYpyGXA7RKxDRs2LNBfmiAKI2xnOXfuXKSkpODEiRMWXroff/xRlItRaXhRTsRKOws5Efv06VNNIpkd5pR7+Ct1MFOnTsX69evRrl07UQealJQkEsbnzp2z6om1JmLPnj2Lr776StSZDx48WLSd3PCikjBkZzXz55aKJX5fqVdT6TryJVtZzGaz0zyxatkJbPXEuru7C/GicoKTvQfKly8vnHPSpEki+5VE5Zw5c2TTX6WlpeHPP/+E2WyGyWTCw4cPLe4FudhR/jxK4Tienp4WokUqBNncsbZ4yPjfHuu1Z4+lJDiVRKyvry8GDRokhB/IeTO1emKtnVev11tcF6UMI0lJSbhy5QoA7eLLXk/s06dP8eeff4q2P3fuHJYvXw4gZ5RAaqeaiJUTrGqpuaR2JyYmon379mjatKniBDxW4KqVLb5586aw/Pfffytu5+rYJWJHjx6N999/H19//TXOnj2L8+fPi/4RxPPIrVu3bC6H6iqcP39e5NlgOzdpJzNx4kTR34cPH7aImzp27JjdntiwsDDRdUxNTZUVpOwMYFtELE9UVJQo7ZDBYLDoFFghYi2cQEpiYiIaNGiAoUOHonr16gDUh/R5TCaTomeI3Z+3Ryoy+GssFRBqoRtSMWQymRRFLG8DP1mJP5+ct1yuU7aWnUCKXIotHg8PD1URy1+bX3/9FSEhIbLnXLp0qWIhAKV2GDJkCJo1awYPDw+MGzcOJUqUECYa8sjd2/y51TyxUtReEuwRsfx1kCKNfeXhr+GtW7fQp08fIWzGx8cHAQEBWLJkCQD534IWTyx7j7Dfh33uyE0ULVq0qOzx2O3UimOwKP2O+TLIPFIRCwDNmjUT/c2KR6PRaOGJlXtpTU5ORrdu3fD9999bfMY+dziOE9kjtWXFihVCCd3jx4/L3h/s70eaFYXHbDYLLwIALMK5Tpw4genTp9uVHi2vsUvE9urVC1FRURg8eDBefPFF1KlTB3Xr1hX+J4iCxoYNG9C8eXPFt+jk5GSUK1cOJUuWVH27Zfnpp5+EdDg8aWlpOHbsmMNikDIyMrB06VLVCjBPnz5FvXr1RFWm2Oo7UvEmFbUDBgzA8OHDRetatmyJL774QrQuOTlZk4iVmzgkJ4bYB7BSfllrsNdZr9db7MPGNX799dcWLylyXqS+ffsCEHfOvEdIyWvNdqKjRo1S9E69//77wjI/QUx6bfjcolKhxHuHpMjlkJV6Ys1ms3Ct+PbhO2f+b7kZ9D179rRYp5adgD0Pa5+9nli+PT09PYUJU3KiUJpTVQ5+hn5AQIAozylfxEOaO1bupYFvVzVPrBRrIjYhIUFTWAF/nPr161t8lpqaqujx46/h/Pnz8f333wvDzHz7857GzMxMi2efFhF76NAhYZlte/Y3kZqaavHbUXrOsllBrKXb4lF6WZE+D6ZOnSr7LGXtloYdSW2Qe+Z9//33ivHAT548EY7/1VdfwcfHR0ifJhX2CxYsEJYTExMxefJkBAYGigqZsM+ws2fPolOnTihTpoxowumpU6dEzyB24uvOnTvRpEkTzJkzB0uWLMGSJUswePBg1VGp/MQuERsdHW3x7+bNm8L/zmTt2rUoV64cPD09Ub9+fdnOjSC0YDabsWnTJty4cQNTpkzB8ePH8eGHH8puy8e6AcDBgwcRERGBevXqYePGjfjnn3+QkpIieiBev34dr7zyCurVq4cJEyZg9uzZiI6OxpgxY9CyZUtRuqYHDx5g9+7dMJvNiIuLQ0ZGBr788kscPnzYYhIG23GazWasXbsWEydORJUqVTBhwgQ8fPgQX3zxBXQ6nXCO6OhoVeEt7Yy1Vr/ik3DzlYK0hhNIycjIkBWxbPycXGUhW8pZAjkdrjXhK00ZJffgfvnllwHIe27Z7zFmzBhhmQ1j+PzzzzV1CHxde+m14QWsVNApDQlmZmZaeKKknlh+HfCfGOLFC/+3Wq5TFrXsBLaGE1jzxPL3gKenp2g7vo1sgc+OoTSMrVbZiefw4cOy2SB4lCqMKWE2mxEaGir7siCF/+2xv+dy5cohISEBpUqVkq30Bvx3j0tjUHlxxno7pfcRb3vx4sUVrxv/nOC/j/S8QE4c5/jx40X7KT07WE8nXwRFmvFBitJLvvQZkpqaKlvWlhWGrEeb4zgLEZueni76no8fP1aM7wVyqhi6ubnhk08+wTvvvAMAePPNNwV7lLh16xaWLFmCp0+fCuEz8fHx+N///idsM2HCBOzfvx8xMTHYsGGDsJ7fhi/Wwnplly1bJizPmjULH3zwATZu3CjrRXYF7BKxZcqUUf3nLLZv345x48Zh6tSpiIiIQPPmzdGpUyfRsCNBADkPSGvei927d6N///6oWLGikNJo06ZNePbsGWbOnCkqB8pWS+nUqRPatGmDiIgIDB48GHXq1IGfnx8qVaqElJQUcBwnEr3Lly/HjBkzUL58eSF/ISuWx48fj549e6JMmTIIDQ1Fw4YN8c4776BNmzaoX78+jh07hi+++AIdO3ZEQECA4AEdMmSIqBTk8uXLUbp0abz77rsAgClTpiA+Pt5qWhepiNU6OZO/vi+88AIA9XACuao+PEoiVpqHle/0UlJScPv2bZuHuqQvGnJI7ZB6h8LDwwXRLr2uUpG8dOlSIZ2U9F7UYvujR49w5coVWQEzcuRIxeFWKXIiVi5OlhdeSp5YrSKWzU4gvR+WL18ue97chhMYjUbRdlo8hFJ4saalIhSg7HWPi4uzKZxg48aNqufJysoShf8o8eeff2LdunU4ceKEsC47OxuRkZGKOY0BsTebhW9/dr2SiH3nnXdQrFgxqzay+0tf5KRCkxeYTZo0UTweP7lPyZs9cOBAAFAMc9T6IsyKPLYvMJvNsuKd/Z5vv/22yLMvxWQywWw2i+K4+ZdANfvYkTV+u7Vr14q2YUfc2GcQf4/wzoGHDx8Kz2qlymBscY2srCyMGDECzZs3x9GjR+36vTkKuyt23bhxA6NHj0bbtm3Rrl07jBkzRvYNxpEsW7YMQ4YMwdChQxEeHo4VK1agVKlSWLdunVPPS+QvJpMJKSkpSEtLU/UochyHQ4cOoWrVqvDx8YHBYMA333wjO3RvNpvx888/yx7nlVdewaxZs9CqVSth0oc0b6dcp3D//n2sWrUKer0e3bt31/Tdfv75Z+EBx4cy8DFPPC1btsS7774rxKq9++674DgOX3/9tcXxpA+TUaNGWeQtlSLtjLXUAGfhCwvIeWJ5gaxWBlZJxErZtGkTOI5D//79UbZsWdn0O2pcv35dU3UsnsuXL4s6ASBHLPEiXypiU1NTRfkj3dzcFL+3FhH75MkTTJ48Wfj7tddeE0SxNJm7Gr///jsePnwoWic3vOpIEavkiZVblxtPLN8xG41GUThBbkSs1spXSvdscnKyTZ7Y7777Dr6+vsjKyrIpGwEP78G8ffs2Ro4cKfrs6dOnFs8uKdZELJshQknEenh4iK55WFiY7EuW1lKwLGqFS549e6bYXrVr1xZGo/gRLin8s89aWIKSiM3Ozpbdl/1urNjUCn+91Z4T7Pe5desWFi9erBrDz3q2+WdXrVq1hOt7+fJlmM1mRafHjz/+iKVLlwLIeRYfP34cJ06cQKtWrfDqq68iKSlJ9v7lOA579uwRwnIcjV0Vu3799Vd069YNderUQdOmTcFxHP766y9Ur14de/fuRbt27RxtJzIzM3H27FlMmTJFtL59+/aK6SYyMjJEDc0/dFq3bi287SvVrGbX07J9y9Y+S0lJga+vryjlDOu5MhgMePz4Mc6fPy8IK09PT9SsWRMtW7bEzZs34e7ujtDQUDx8+BC//PKLRccyaNAgDBo0CIGBgareCBZ26KdDhw7o0aOHRYlHJeRqritx69YtdO3aVfP2LNI0PUps27ZN1QsAKE8I0grviU1JSbHwHPAPTjUBlJ6erml27Lvvvgt/f3+hLaTCzBrnzp0TwjOKFCkifO+qVavC09MTkZGRAHJemv744w/Zwi3u7u7CZDHpwz4xMVG4/7y9vWEymeDt7S3bEWnxAH366ado1aqV8Le/v79dsdTvvfeexTq5Dp33YPIviryo4wWvVhEL/Hd/yolJqXhhxSkrgIEcccSfV618rJubm7Dd5cuXNb0USeFFsNbfg9I5EhISVCuEyZGamorSpUur5vZMTU2VFcFKE7mAHJFmbaQyJSUFJpPJYkieLUnr5eWFZ8+eISUlBXPnzkVWVhamT58u3EcGg8GiJLOHh4eFwEtNTRWOqeTJlhIWFiYsly1bVvRSnp6eriiGP/74YwQFBUGn04HjODx+/Fh0LOC/51PJkiVhNBoVPbaXL18WPKbs9czIyIBOp0P16tVx9epVod2TkpIQFBSkmtpLDaPRCJPJpPkaAcDkyZNRs2ZNADm6qHbt2sKkPCAnrIG/9nx/GhAQgFKlSiE2NhaNGzfGgwcPZB1Fnp6eePbsGT766CMMGTLEIqPHvn37UKRIEdSqVQu1atXChQsX0KNHD4wbNw7vv/++EMpQt25dZGZmwsfHB1euXIFOp0PdunWRmJiIZ8+ewWQyCf+0vkzaJWKnTJmC8ePHY+HChRbrJ0+e7BQRyzdAaGioaH1oaKjim8OCBQswa9Ysi/VSTxdRcHj27BlOnz6N06dPK25TokQJiwk6WgWsHFoFrK2wpVLzE6UykwaDQdMkNlbcSCtj8ZMg1DwKsbGxFh5PJaZPn65pO2uwItbDwwMjR47EsGHD4OPjg6ioKEUvr8lkUvRUnz9/XvCm63Q6REVFKQo/6Ux3JdhO1ZFDdnIem0uXLomEDC+Y+Y5eSUDzIoEnLS1N6MzlziNd9/DhQ0HUx8bGir6n2WxGVFQUgJzwHyUP5507dwTB8Omnn8puYw2tIxBeXl5IT09X3P7ixYuKHjG1GPEHDx6oZveZOHGiIFJYOnTooPryLDdiA/z3+46Li0NUVJSFaHj06JFw7fn7+M8//8THH38MAGjatKnw+05ISECfPn2wbt06tG7dWvgddOzYUTTpKDY2VjimtWT8PKxdzZo1E4nY1NRUIXbcYDBgy5YtQkqwxMREXLlyBX5+fkhKSsLJkydRvnx54bvcvXsXnTp1ApAzGrB+/Xq0aNFC1obIyEhERUXh8ePHFjmuo6Ki8OWXXyIjIwM9evTA06dPhXvgl19+sThWx44d0bx5cxQtWhS3bt3Cli1bLF40TCYToqKiBA9ww4YN0atXL3zwwQeq14p//laoUAF9+/YVvKVATghBVFQUsrOzhWdffHw8ypQpIzgR+CwF/v7+ope0o0ePomHDhsjMzMR3330nOHqk7ctmqPrnn38s9Jd0ojMAxVRhWrFLxEZFRckG+Q4ePFi2WoYjkXrz5Gqq83z44YeimMGkpCSUKlUKW7duFb0hsg9gpWXpOQvScn6fX26Zn8QUEhIieG2kM8mzsrLg7e2NatWqgeM4VKlSBQ8fPsTp06exe/duBAYGwsfHB4mJiUKA/8iRI4Ukz99//70QIK+FyZMniyZcyTF16lQkJydj/fr1MBgMNk8u0sKGDRswbNgwzbF5cnzyyScWqbKUUBrClKbCUqJly5b45ptvAFiGDfDXR21IUHp+Pz8/eHt7C57WihUrCrP/tXZ81ihZsqTQcYSFhQm1zt3d3REeHq4YBhAYGKg4QaZ79+7CfoGBgUIxGDmBvmPHDkXbunXrhj179gAQC1drycttQe7eKl++vMhTyDsMeM+o0pCrj4+P6Hfg7+8v2K3lHvb39xee4RUrVhTFZPv5+QmT4uS8hTx169ZV9UhqoWLFipq2CwsLQ3R0tKJQVYsp559TSpQtW1bxs88++0x2fe3atVWPqZRmKTQ0FLGxsdDpdHj8+LGFwGjSpAkqVKgAIKcd4uLiRPHKBoNB8AyXKlUKQ4YMQffu3dG4cWPhXuFHaXj0ej0iIyNRt25dxZhu6ctz+fLlcebMGVy6dAnZ2dnYtGmT8FlmZiZKly4NIMd7yYr82rVrIzw8HCEhIUhKSkLv3r0B5MSNDh06FFu2bBG2vXnzJho1aiR6IVu1ahXGjx+P7OxsmEwmhIeHiwoEAMCwYcOErBZAzr389OlThIaGIjw83GIIvU+fPvjyyy9Fk+VGjRqFUqVKibZ78OABLly4IIw8FytWDPXq1ZO9XnJUqlQJtWrVwh9//IF9+/ahW7duiImJQc2aNYUXHp1Oh4YNGyI8PFzQc/ykuRIlSkCv1+Pp06do2bIl6tevjxEjRmDdunXYuHEjEhIS4O7ujs2bN8PPzw/Hjh3D1atXMW3aNDx9+hTdu3fH3r17kZ2djRIlSmDjxo3Q6XSIjIxEiRIlYDabERQUBC8vL1y/fh0hISHw9vaGm5sbDAYDDAYDMjIy0KZNG6vf1S4RGxISgsjISFSqVEm0PjIyUlNwtz0EBwfDYDBYeF3j4uIsvLM8bJwUS4cOHTTVfCacB/+mGR4eLho+tEZoaChq1aqFIUOGWN22b9++Qhzh9evXMWLECLi7u4uGQngbgJyO8MKFC9i7dy/69u0r6ykNCQnB3Llz8emnn4LjONy7d8/iAcTTpUsXIe526dKl+PPPP7Fr1y6L7fz8/HDp0iW89tprGDlyJPr374+SJUviyy+/tDojtG3btggODoanpyeKFi2K2NhYLF++HMWKFcOaNWsUOzAWpWHR0NBQTSK2Vq1agndKOhTLe1HU8jmy3qlSpUphxYoV+PTTTwURW7p0aYcXV2GfGQEBAYL4fPr0qezMfR6j0ag6oYof2ixTpgwMBoNIDL/77rua6pT36dNHELGsOJQKAkfg7u4OjuMEscn/r9PpBIFiLSbW19dXZKfBYFCc2CUHm+7Ly8tLJKSNRiPeeecdzJs3D6VKlRJNWmLx8/OTHWq3BalId3d3l70PihUrhujoaMWh3tdee03xHCkpKUJ4ktFolE03ZwsTJkxQfH42atQIJ0+eVBxu50VsWloaunXrZvE5+9LE38dsCE98fLxwfTw9PeHp6YkOHTqIjsHmaQZyvOm7d+8GYDkJiUc6+lOiRAnUr18fderUsYgvffbsmXDPenh4iPRHpUqVYDAYLPr6kSNHYsKECRbX2t3dHZ6ensLLybBhw1CpUiV07NgRaWlpMBgMwrX09/fH1q1b0aFDB9H1ZzN6GAwGi2dr165dLa6JkhZhHTB+fn6qjgApvF4CLEf8+LLjAQEB8PDwQHBwMHx8fJCamorRo0cDABo0aIDevXvj888/x9q1a2EwGNCoUSOsW7dOeNlp2rQpihQpAoPBgDZt2qBNmzbo27cvnjx5ggoVKgje75deegkBAQEAcvorKWzIFIvWkCC7ROw777yDYcOG4ebNm2jSpAl0Oh2OHz+ORYsWifIcOhIPDw/Ur18fBw8eFJWYO3jwoOZJNEThg+9Mq1SpgsOHDwMA9u/fj7feeguff/45qlatih9++AEnT55Er1694Obmhho1agDImTAhzb3JChidToeSJUti/PjxWLduneihuHHjRgwaNAj//vsvfvnlF4waNQo9e/ZEs2bNYDKZRDNRedHKdtDt2rVDu3btsHXrVowZMwZr1qwRPluxYgXGjRuHQYMGqc5svnr1KsqUKSN4AocOHYpmzZpZZApQE7FaCAwMhLe3N9LT0y3qsPOopcDhQz38/PwE7yjrbSlRooQmO2yB7ex0Op3IvilTpigO3cslN1c7Prut1k4oJCREVvCWLFlS0/624O7uDrPZjOzsbGRnZwteKKPRKIpF5beVQ9opsym2tExUkmYnkOaJ5bNBPHv2THa43dPTEwaDQTHeVCvSF62wsDDZxPXSuEpbSE5OxtKlS9G2bVu4u7sLmUR4bI3z5if+bdiwwaJSnPSlRxrbz38PuRytfDlXHv7aSIuDsOWC5VCb0MlmcFGDFWFSwcdxnPACxb9gbtq0CW5ubsJ9w6eRYpEK2J07dwIQh7qwMb18OBR/rjJlyqBz584Wx+W350tIT506VfQ5bxOLt7e31dAtHx8fxfuuTp06Qjw/D/udK1eubBHyA4j7suLFi4tGuXr37o1u3bqJXm6koSxynuGAgABBsPr7+6N9+/aK38lR2JWd4OOPP8b06dOxatUqtGzZEi1atMDq1asxc+ZMi0ZzJBMmTMCXX36JDRs2ICoqCuPHj0dMTIxFInaCUKNjx4549OiRUKO8e/fuWLBggUV+1GXLliEqKgrXr1/HX3/9hWXLlsnWIl+2bJmQWuvTTz9Fly5dhLis6tWrY+LEiTAajShbtizGjx+PcePGiUSaWoEQvV6P1atXC8P169atw9ixY8FxnNXUPG5ubiIRvnz5ctlzKYlYuQeu0nZ8JydNBM6j5iXjhS8r+KQPWDX8/Pw02cnCitgLFy6IROzy5ctx5swZ2f20TqjjhR3biYeGhmra38vLy0IYAvJejNzi7u4uCM7s7GzBM2g0Gi0yDCiJWKmol07Osoa1FFv8vSMt/8nD33vWcoVaQ3qPKnnI1Ib8pUjnX/Tr1w/ly5fHqFGjZMMOlOLT5Thx4oRwH8uJRenv4pNPPhEJFT60QeqprV+/PmbPni1ax19jdiQ0LS1NlJ1ADrUXvtWrVwvLer0eHTt2hJubm+gFH4AgigD5ZxIb2w7kXGP++SvdX45vv/1WMQ+vVMTyYl/pmcO3Q2pqqmyuYjn7dTqdSMBKrz2Q8zwJDAzE7NmzRRlLgJwYcD4umIcVsR4eHjCbzbhw4YJIVLLPWOkLj5z4rFWrluh7q6U+y0vsErE6nQ7jx4/H3bt3kZiYiMTERNy9exdjx45VjE91BH369MGKFSswe/Zs1KlTB8eOHcO+ffucmpuWeD7Rep9WrVoVFSpUQOPGjTF+/HjFQgB8pz1mzBj89NNPqsPnbm5u6NevH3r37o0WLVpoitvt378/4uPjbX5hY4cEfX19ZcWR0gQVrUUP/P39hYe3UviD0WiUrSQE/OeJZTti9mGpFuvYrVs3UYclR61atSzWsbGJb7/9toUAUirRyU98sBbmwdvPCoSiRYta3BfffPONRZJ3Hx8fi3aaNWsWDAYDLly4gJCQEJtEIsvKlStFf7OiUSpipSEBSkJFaqutIvbBgweqKbb4v9mQkvLlywvL/H3jaE+sklCxZUKm1AvIerPkRLfWCY6AOBZWKmIPHjxo0S5BQUGiF0L++0knXf72228WfSp/fFZkp6enWxWxap5Ynq5du+L06dPYt28f4uPjLZ4TrAiVeyZJRawUay/j7PH5GFZ+Vr+tIpb1xNpiCyso+YIHcsf9+OOPsXDhQtG9ExwcbHGfyb2A1ahRAw0bNpQ9JytiAwICZJ0OBoMBv//+O0aPHo1Tp06hSpUqst8lr7E7TyyPn5+fXZ4Qexk5ciRu3bqFjIwMnD17VnE2IUG4Ojt27MDRo0c1Peh1Op1dcdx16tTB3r17hTd1td9qp06dhAd6pUqVNA/je3h4CAKAf3hLQxGMRqNiNTRexLJeG7azknvwN27cGEDOsJc1zxgfgye1md9vyJAhmr14fM5NtbhH4L/rzBbMCAwMtOgc+vfvj2XLlomEEZvGi4e3r0aNGoiLi9MUWyuHNPTK3d1dEAGPHz+WFbH2hBNofQECcrJ/8MO0cp5Y6XmbNWsm2oa/9xwtYpUEhy2eWNYmaTU4uXtu7Nixdh2bfYZUr14dbdu2Ff2edDodihQpIvqOrIhl18t5iOXCCbSIWC2hNzVq1EC9evWg0+ng5+dncSxrxVf4+1fpN2yLiG3ZsiU4jhMmxfL28yNtWj2xShN+lfZjf0Ny84qk14Bt75CQEIu+QS6EAoDICcKG+rAvN2qx9/Xq1cPKlSttmmTmbDSL2Hr16gmdTd26dVGvXj3FfwRBuA5du3YVOlC1B3rjxo2RkJCA+/fv49SpUxZxcUCOuJVLVyYV4nwJTx6j0YhevXqpppZij8F68qQiZuLEidi2bRt27dqFfv36YdSoUejSpQtmzJghe1w5r3j9+vVx8eJFPHnyBIGBgaqii41Vbt68ueJ2LHyn1LJlS2Fdw4YNRbbwHTe7Hd+BSDt/aQetRbDJbSP1kLL2NG3aVMjSwIYT8MJW6RrlNpyARc0Tq7SNo8IJ1EZPWGzJgsC2gXQitFz72FLsQG5CEQBhoikrjAICAqDX60UjUHxfnZqaKnjVJk2aJDtKxV8b1huuRcRquabSZwd7LDc3N9nRIxZrnlhreX/Vwg34c/OVwXgRq2QT/1tQGt1Smvi+b98+6PV6VKlSBXq93uIlWSpi2QmHQUFBFve+ksODFaust5jP8ABonwvhKmh+Xe7evbtwobp37+7UsAGCIJyDWkfPCwN+AoHcw91oNIrW8zGe0o5ImkqI90CqeYLZY7DPF3boGMgR4qVLlxYevAEBAfjpp59w79492bzQrPezSJEi2LhxI1588UUA6p4inU6HuLg4BAcH49atW1i7dq0wexewrFLEwn/PjRs34oMPPsDEiRPh6+urGDKxfPlyBAUFYciQITCbzRadpLSDZtsxMDAQvr6+FpOQvL29LeyTikupZ5jfnvXE8sLKWeEELFpFbG49sXw2Dek6FiXBYS1NFouHhwd27tyJVatWiSZnAuq/RWnGB/bccnHn7H3Ffw/23mY9c7t378b58+eFTAJ8IQDA0lMvPSaLFhErpxMaNGggijlXE7F8sQI1rHlireUIV/P0stcwNTVVaBNrnli2StqxY8dw8+ZN6PV6xXNVq1YNp0+fFka/li5dimfPnmHv3r0ALJ0P0mwgUtQKk/D3Fito2dGFgiZiNXtiZ8yYITTQzJkzMWPGDMV/BEEUDNjhS7kHnzSGUtpZKeUQlXZ6fAcjHVJlYY/B5mJt3ry5kGAdUI6zU/L6sEKtfv36ouwmLNI0Q4MGDRIES5kyZbBo0SJRhoAtW7bAy8vLwsMG/CfsypUrhx07dgixaKwXj30ZCAgIwNKlSwWPmFI4AY90mPq7776zsEHuOkknlnl6egoJ36Xnk3petU7sYrMT8AwbNkx2XylawgmURKwtnli5e8VoNIoEk1KOSqmA6d69O+bPn4+zZ89abOvh4YGePXviyJEjFuE5avZevXpVWGYFqNxkIUBexLL3EOuZ69GjB6ZPny5qNz4rgjVxxpKenm41/ZpcijJfX19RFiM1ESvnUeTzvfJY88T26NFDdj2P0tA7kOMJ5o/LZnGwFhO7atUqADm/r+bNm2PgwIHo37+/qh316tUTBGSpUqUwZswY4TNrIRW2sHz5ctSuXRvz5s0T1rEiVutok6tgV0xs+fLlZd3lT58+tfCaEAThurAPcLmOaPTo0appf3ixwh6HT2LNwnfYasP27LBf7969sWbNGpw7dw46nU5UqUspgb6SuGVFrFoeznHjxgnLH330EZYtW6a4LQD06tULycnJGDhwoMVnSp0cmyZHzcMkFYZqotbd3V32fHLXQyoujUajbAycVNDx55FDiye2SZMmaNCggez+0nNIU2xpDSew1RMrxcPDw6Lk7a5du9CnTx9hGL1q1aoW4rNMmTL48MMPLULp5MQ8i5qIDQ0Nxfnz5zFnzhwsXrxYWK9FZPLLbLvICTVWoPFiU2mYXO56paWlCb8npRdIuXtGr9eL7m/pvlJPrJQvv/wSv/32G6pWrQrAuifWWvy6tbhddnKX1phYHmuhEGqwwlUqYvmXQrZgg9ay0EOHDkVkZKRIq1WoUEFoB3tLoecXdonYW7duyeY0y8jIEMrNEQTh+rAPWaVOQG1CGe/ZY49z4sQJRRGrBpurUq/XY+TIkUJKMFb8KhUi4As+SGHFmFIKMACiGdnz5s2zmpoHyBGFcp2gUkoxVsQqVQQELDtJNU+skoiVi3+WemKNRqOsrUajUSi4wJ5Hi61yIlYqXAAI9dSl55B6YqXi1FmeWLZIA3+ebt26Ydq0aShbtizu37+Pc+fOydojh4eHh+qLCmvvnDlzhOXp06dDr9ejZs2amDZtmug+5L2KUnHEfh/eHvZ6S8ua8kjbREmcWROxSvd7jx49UK1aNdELol6vF9kvtcGaJ9bX1xdt2rQRzs1ndFB6iVGyTa/XY/To0VbDFeRErLWYWKW/bYEd7pc6B1etWoWLFy+ib9++wjotzysl/Pz8cPToURw6dKjAZXuyqdgB+1D79ddfRW8HJpMJhw4dcpl68ARByMMm1mYfskqdgJr3lP+MFS0VK1a0SL6tpZqStQkYISEhePTokezwPU9UVBTi4+MFL40UaUJ4lvLly+PAgQM2J7OXdmhVq1ZVPD8bb6Y2xC4VoNIJP6wAMhqNsuKjRIkSOHfunGidVFy6ubkpilipBz434QQGg0G0P+sFlNpjLSbWzc1NdE/y9ssdb9myZaLS49J9pHZLRSwLf19I44zZ7dauXYuRI0cq2sPCtiFb8lY67Mwe5+WXX8axY8csSuSyHkD+5Yi9Ly9fvixrg4+Pjyhm1JZwAi0i1t/fX8iMwpek1+l0mkWs2lA/PyKzf/9+AOovMQEBAXj69Kmo4MPp06c1TURnRazWmFjpvvZQokQJbNiwAenp6RZxqh4eHqhevbpo3caNG9G1a1erpdOV4OcJFDRsErH8W6BOp7MYQnN3d0fZsmWxdOlShxlHEITjGTlyJFatWoWOHTuKOhMtkzOUvIfsNnq93i5PrDRfpZQff/wRkZGRFjFxLCEhIQgJCREmj/AxrPPnz8f58+etVhRs166dVTulsNfwzJkzqFOnjuIwMiuQ1Up0SztJaTyltDSrXKfKzjjmkROxSvGh0rZWErHStpXzxEoranl5ecm+HEmFpJbsBLyAkrvH3nrrLWzZssWieIXUI823hZqIVVrP/s2W0LRFxLI5N6UvUexxPDw8ZGMWpS81gNiTJ1fuGrCeBYNH7h5JSUkRJsfZUvI3NTU1V55YngoVKohGftWu99mzZ/HTTz9hyJAhqF27NpKSklRfhll4W1NSUjTHxPJoSZ+oxttvv6152y5duiAhIUFzkZrnBZvCCcxmM8xmM0qXLo24uDjhb7PZjIyMDFy5cqXAxVMQRGFj4cKF2LZtG7Zs2SLqTLTGVAEQCgy89dZbAP5LqcVPYFITsXyO04ULF4q2mTZtmuo5GzdujBEjRmia+b5582a8/fbbQqnhDz/8EFu3brUpf6lW2GvI1iyXIzw8XFhWmwUs7YjYJOWA+HqePHlStrPUEk7A14uXYjQaLcI2lO4PaWJ3JRHL7u/v76/YFvam2FJKKSZ3Hlag37lzR4h5lQpoOXQ6nYWw5JETk0qw16hMmTLYunUrfvjhBwvPPivklK4Z+xLJTx6sWrUqXn31Vbz99tuKkxnZ32mzZs0Uh9ZZccbfixEREeA4Dnq9XjaMRwo/OtG7d2/R8ez1xEoFntr1Ll++PMaMGQMfHx+cO3cOUVFRmvPb5yYmVil+31nwqdQKE3Y90aOjox1tB0EQeYS3t7cgQrV4Yln4Tm716tXo2rWrcBze88nnqFQTscOGDUOXLl1QvHhxLF68GPHx8di8ebNsSV97qVy5smzMpTNgBZW1WcQvvfQSGjRoAE9PT9ViEmzHzr8osLBtde/ePVnxIScA5Dyxcp2e0WhE//79sW7dOmGdkoiVesu0hBPkVsSyYkHNE6skYtk5HWzGCS2eWN4uuXK8rA3Wfk8hISF4/fXX4e3tjaCgIMX7/8UXX8TIkSNF1ffk8Pf3R1JSkpDBQK/XK3pgedgXGLWYSvaFiPda8y85oaGhmp4dR44cwZEjR9C7d28cO3ZMWC8VhFo9sdKRDK0T+2z1VNoSEysdvVCq3EU4DrvdEqmpqTh69ChiYmIsYoTY1BAEQbgurFjS4onlH9LBwcEW4ootFyn17kn/5gXcgQMHEB8fb9cwvqvACkhr3h0PDw+cPn3a6jH1ej1atGiBY8eOyQ4fsx32u+++K/qsU6dO6Nmzp6KoY3Fzc5NNsG80GjFz5kxNInbAgAGiUplKE7vY/f38/BSPpyWcgL1v1WJipeKZp127drh06ZJFzldbRKzcduw1t+b11+l02L59u+o2/HbSHLNy/P3333j06JFQzU4L1ip18aiVgVbLXsISFhYmTERq2LChEO4gnUfDXls1m7RUQHME/L3Gen6V7JJOeFcqekA4DrtEbEREBDp37oy0tDSkpqYiKCgIjx8/hre3N4oVK0YiliAKCLZ6YrUOVWmNiZXWSS+I1K9fH97e3lZDCWzl66+/xvHjx616qPkJMzt37sTZs2cxd+5cRYEkbT81ESsd5pUTdY0aNRJSU/EduJaYWD8/P4d5Yi9evChsK0XJE1upUiXcuXPHwtOXWxHLrlfLPuEM1CYUKsG+XKp5KFkRy3qu7cXX1xdXrlyB2Wy2aDd2BEEpEwmQE4LRrFkzoQpgbssOKyE3OUtJxHbr1g2tW7fGkSNHnGILYYldwRPjx4/HK6+8gvj4eHh5eeHkyZO4ffs26tevj08++cTRNhIE4STYTlyLJ9bRIvZ5oGjRorh+/TrOnz/v0OOWK1cO/fv3l22X8uXLY+jQoZg2bZogRHr27Il58+YJnmHpfv3797cY6nd3d5cVW9nZ2aoTmHj4FGjsfaElJtbb21tRxEqFpDURe+3aNQC2hRMYjUaULFnSQqBI88QqIc3TK7e+ICCtZqcEK2Jr164t+kw6EqAVX19fWeFsMBjQu3dvBAcHK8by8jRp0kRYdta1l1bCA5RDL4oUKYLDhw8LOsharmki99jliY2MjMTnn38Og8EAg8GAjIwMlC9fHosXL8bAgQPRs2dPR9tJEIQT0Or14LFHxOr1eqdMqHIl5AoGOBOdTof169erbsOKKzZG2M3NTfCaurm5yYrYmJgYizaTE3V8/k+DwSDcP1piYj09PTXN/tdS7EDNPmnKLrVt+e3l7JCiJGLz0xNrD1pFLPuckBat+PTTTx1u1/fff4+srCyr3lVWTDpLxD548MBinbWsAxMmTMDAgQNtKlFM2Iddnlh3d3fhbT80NFRIpFykSBHFpMoEQbgerNhk0/wooXW4nO0crSUTJ5wDK646duwom9NXScQmJCRAr9eLXlqkonbFihXCZCNrnlhpTKynp6fFfcGHtkg9m9bKzn7wwQfCtlKUXqDkivUA2kXs06dPhWVWaGlJR+dKsNdMLZygePHiWL58OZYvX47SpUtj7dq1AHIErS1ZTbQizQChBCtilaqG5RY5b6q1Z5pOpyMBm0fYJWLr1q0r5N1r3bo1pk+fjs2bN2PcuHGoWbOmQw0kCMK5XLx4EZs3b1adLd+oUSMA2ocO2Q6lIHTmzyNKHkJWqLm5uWHo0KHQ6XTo1auXsJ5/CVE6BiAWw1rCCdiZ2kajUVQCeNSoUcKENy0Tu1gRwU8WktrH34NyIlYp9ZFWERsXF2d1u4Jw39+7d09YtiYCx40bJ3jehw4dinXr1mHz5s3ONM8qbDhUbnOyKtGhQwdkZWUJWR/Ya0bkP3aN8c2fP19INTFnzhwMHDgQI0aMQMWKFbFx40aHGkgQhHOpWrWq1Q533759iIyMFCVzV4MVFHIThwjno5TLVCrUSpUqhZSUFHh5eWHjxo2YPXs2Vq1aJXyekZEhLLOwwpUVrUrhBNu2bRP+9vT0RJ06deDp6Yly5cph9erVIptYu+UqdrEilP9uUk/s1KlTZe0GgFq1alms448tZ4cazvBE5hXsUHnnzp017+fu7o7hw4c7wySbyAtPLJBzX/z2228wmUzPfWhUQcPm1uA4DiEhIULJs5CQEOzbt8/hhhEE4ToEBgaidevWmrenEIL8R64gAGAZTgD858UaPHiwqDQv22HLCVMeqSdW2tEbDAbRy4ynpyeCg4Nx9epVi5ybWiZ2sV5c/l6TpnHjBQ57vK1bt8LLywsvvfQS5GDvW62z3ZViMQuCJ/bVV1/FhQsX0Lp1a9XiG64KK2K1Fi+wF2k1OcI1sEvEVqpUCf/++6/msm0EQRBE3sKKK3bSjjScQA1WQPKFEXgxaks4gV6vF4k63ja+OIbSOT08PCzEpJLnU1oAYciQIYLdPC+99BLKly8vuz8gFp5aRaySB7AgiNjx48dDr9djwIAB+W2KXbAili2zSxQebI6J1ev1qFSpEiXxJQiCcGFYEct6qWwZMmc/NxgMIk8lK1S1hBO0aNFC+FtNIGqJiZWLf9TpdGjevDmKFy+Ox48fC55ZVmRaE+3shC+1a8PHiAPKGTsKgogNCAjAjBkzLAoOFBRYEasW0088v9g1sWvx4sWYNGmSkGSaIAiCcC2UZp7LhRMoIQ0nUIqDZdfLiVg3Nzd89tlnom2UkHpipeLZ3d0dw4YNQ5s2bUSxtEBOadObN2+K8r8qCXg52JAHNRG7f/9+2eMTeUtoaCj69euH/v37F1ghTuQOuwI83nrrLaSlpaF27drw8PCwGE6Jj493iHEEQRCEfUgrZPHkNpyAR0nE8vnDWTw9PUWeMjUvpVxGBHd3d6G8ubu7O/z8/PDbb79Z7Ct3bltELGsXW4VMSpEiRbB48WLcv39fseocxYU7H51Oh02bNuW3GUQ+YpeI5cscEgRBEK4J64llJ0/ZK2LVPLFSQSsVkkajUXQuNYEnF+7AC1ipTVqQ5qdVw5ZMGpMmTVL9nEQsQTgfu0TswIEDHW0HQRAE4UBYEcuKN1vCCbTGxGrJE8vua0s4gdrnWmC/uzQTghRHxrFqrW5HEIT92P0ru3HjBqZNm4a+ffsKiZ/379+Pf//912HGEQRBEPbBhnmxSeFt8cRKt7XXE1u6dGnNdksndkmxVcT27t0bbdu2xbRp06wKS0fmNCYRSxDOx65f2dGjR1GzZk38/fff2LVrF1JSUgAA58+fx4wZMxxqIEEQBRNWOBF5j4eHBzZv3owNGzagbNmywnpnhBOoeWK/++47C0Fnqye2Y8eOmm2WEhYWhoMHD2LOnDlWt3WEiOW/m1JBBYIgHIddInbKlCmYO3cuDh48KHpTbt26NU6cOOEw4wiCKLh89dVX8PHxwY4dO/LblELLm2++ibffflu0Lq/DCeS8qWXKlFE8J1uNi88ywIYEOLNCliPCCf766y8MGTIEX3zxhQMsIghCDbtiYi9cuIAtW7ZYrA8JCaH8sQRBAMipBpSUlETDqi6Gs8MJpBkC2NjcPXv24NSpU3jzzTcVz1msWDFhmU8Nllci1hGe2EaNGonyyBIE4TzsErEBAQG4f/++RV62iIgISjhMEIQACVjXwxnhBOx2ap7YV155Ba+88orqOUuWLIldu3YhNDRUtqSsUplXR+DImFiCIJyPXT3Mm2++icmTJ+PBgwfQ6XQwm834888/MXHixAJbvo4gCKIw4IxwAlao6vV60d/2iM5XX30VTZo0Ef5mJ6lpLQdrDyRiCaJgYZeInTdvHkqXLo0SJUogJSUF1apVQ4sWLdCkSRNMmzbN0TYSBEEQDsIZxQ6kIpatmOUI0cmWmXWmJ7YglIolCOI/7AoncHd3x+bNmzF79mxERETAbDajbt26qFSpkqPtIwiCIByIM8rOssLSGSKWPLEEQchhl4jlqVChAipUqOAoWwiCIAgn44yYWDVPrCM8p3nliSURSxAFC80idsKECZoPumzZMruMIQiCIJxLbsIJtMbEOlp05pUnlsIJCKJgoVnERkREaNqO6kUTBEG4LkqpseTQGk4gzU7Ailg+TVZuYIUrhRMQBMGjWcQeOXLEmXYQBEEQeYAzwgnY4+j1eoSEhAgpGENCQnJtM5tii/XKOprMzEynHZsgCMeTqySO169fx6+//or09HQANBRDEATh6jgjnEDqiXV3d8eFCxcQFRXlkOIEbNlcZ4rYF154wWnHJgjC8dglYp88eYI2bdqgcuXK6Ny5M+7fvw8AGDp0KN5//32HGkgQBEE4jtzkidXqiQVySsY6ahJWy5YtMWrUKKxevdohx1Nix44daNeuHY4fP+7U8xAE4RjsErHjx4+Hu7s7YmJiRLFPffr0wf79+x1mHEEQBOFYpN5VNbTGxMqJWEdiMBiwevVqjBo1yuHHZqlWrRoOHDiApk2bOvU8BEE4BrtSbB04cAC//vorSpYsKVpfqVIl3L592yGGEQRBEI6H9Y5aG+q3J5zAmjAmCIJwFHa9Mqempoo8sDyPHz92ag4/giAIInewz2hrz2tWkHp6eoq8rOwyux31AQRB5BV2idgWLVrg22+/Ff7W6XQwm81YsmQJWrdu7TDjCIIgCMfCikxr6arYybpGo1FTii0SsQRB5BV2jfssWbIErVq1wpkzZ5CZmYkPPvgA//77L+Lj4/Hnn3862kaCIAjCQdjiiWXzphqNRsVwAtYT68w8rgRBECx2eWKrVauG8+fP46WXXkK7du2QmpqKnj17IiIigsrQEgRBuDC5EbFaJnaRJ5YgiLzC7gj8sLAwzJo1y5G2EARBEE6GFZ/WBKfJZBKW9Xq9poldJGIJgsgr7PLEbty4ETt27LBYv2PHDnzzzTe5NoogCIJwPrZ4YqWQJ5YgiPzGLhG7cOFCBAcHW6wvVqwY5s+fn2ujCIIgCOfAelOtxa9KP2cnepGIJQgiv7FLxN6+fVuoi81SpkwZxMTE5NoogiAIwvlYy+n64osviv5mPbNK4QQ0sYsgiLzCLhFbrFgxnD9/3mL9P//8g6JFi+baKIIgCMI5qIUISHn11Vfx+eef49y5cwDIE0sQhGthl4h94403MGbMGBw5cgQmkwkmkwmHDx/G2LFj8cYbbzjaRgDAvHnz0KRJE3h7eyMgIMAp5yAIgnjesUVk6nQ6DBs2DHXr1gVAIpYgCNfCruwEc+fOxe3bt9GmTRvh4WU2mzFgwACnxcRmZmbitddeQ+PGjfHVV1855RwEQRDPO4MHD8ZXX32FTp062byvUjgBiViCIPIDu0Ssh4cHtm/fjrlz5yIyMhJeXl6oWbMmypQp42j7BPh0Xl9//bXTzkEQBPG84+/vjwsXLti1r5Inll1PMbEEQeQVdueJBYBKlSqhUqVKjrLF4WRkZCAjI0P4OykpCQCEEAgi/+CvP7VD/kNt4Tq4eluwYpXjOMHOrKwsYb2bm5vL2m8rrt4ehQlqC9chL9pC67HtErG9e/dGgwYNMGXKFNH6JUuW4NSpU7I5ZPODBQsWyBZkuH79Onx9ffPBIkLK1atX89sE4v+htnAdXLUtnj17JixHRUUJy/fv3xeWb968Kars9Tzgqu1RGKG2cB2c2RYpKSmatrNLxB49ehQzZsywWN+xY0d88sknmo8zc+ZMq1W/Tp8+jQYNGthsIwB8+OGHmDBhgvB3UlISSpUqhYoVKyIwMNCuYxKOwWQy4erVq6hcubJoWJLIe6gtXAdXbws2lVZ4eLiwzD5Pq1evnqc2ORNXb4/CBLWF65AXbcGPnFvDLhGbkpIiG/fk7u6u+cQA8N5771nNZlC2bFlbzRMwGo2ykwwMBgP9CFwEagvXgdrCdXDVtlCKiWXDtlzR7tziqu1RGKG2cB2c2RZaj2uXiK1Rowa2b9+O6dOni9Zv27YN1apV03yc4OBg2cpfBEEQhOvBiliWoKCgPLaEIAjCThH78ccfo1evXrhx4wZefvllAMChQ4ewdetWp8XDxsTEID4+HjExMTCZTIiMjAQAVKxYkeJbCYIg8gClWNdevXph3rx56NevXx5bRBBEYcYuEdutWzf88MMPmD9/Pv73v//By8sLtWrVwm+//YaWLVs62kYAwPTp0/HNN98If/PJt48cOYJWrVo55ZwEQRDEfygN8YWFheHevXvQ6XR5bBFBEIUZu1NsdenSBV26dLFYHxkZiTp16uTGJlm+/vpryhFLEASRj6jFqZGAJQgir3FIHpTExESsXbsW9erVQ/369R1xSIIgCMLFoAk1BEG4ErkSsYcPH0a/fv3wwgsvYNWqVejcuTPOnDnjKNsIgiAIF4JELEEQroTN4QR3797F119/jQ0bNiA1NRWvv/46srKysHPnTpsyExAEQRAFCwoZIAjClbDJE9u5c2dUq1YNly5dwqpVqxAbG4tVq1Y5yzaCIAjChSARSxCEK2GTJ/bAgQMYM2YMRowYgUqVKjnLJoIgCMIFIRFLEIQrYZMn9o8//kBycjIaNGiAhg0bYvXq1Xj06JGzbCMIgiBcCKU8sQRBEPmBTU+kxo0bY/369bh//z7effddbNu2DSVKlIDZbMbBgweRnJzsLDsJgiCIfIY8sQRBuBJ2vVZ7e3tj8ODBOH78OC5cuID3338fCxcuRLFixdCtWzdH20gQBEG4ACRiCYJwJXI9NlSlShUsXrwYd+/exdatWx1hE0EQBOGCkIglCMKVcFiAk8FgQI8ePbBnzx5HHZIgCIJwIapWrZrfJhAEQQjYXXaWIAiCKFx88sknyMjIwLBhw/LbFIIgCBKxBEEQhDaKFSuG7du357cZBEEQABwYTkAQBEEQBEEQeQWJWIIgCIIgCKLAQSKWIAiCIAiCKHCQiCUIgiAIgiAKHIVqYhfHcQCApKQkGAyGfLamcGMymZCSkkJt4QJQW7gO1BauBbWH60Bt4TrkRVskJSUB+E+3KVGoROyTJ08AAGXLls1fQwiCIAiCIAhVkpOTUaRIEcXPC5WIDQoKAgDExMSoXhTC+SQlJaFUqVK4c+cO/P3989ucQg21hetAbeFaUHu4DtQWrkNetAXHcUhOTkbx4sVVtytUIlavzwkBLlKkCP0IXAR/f39qCxeB2sJ1oLZwLag9XAdqC9fB2W2hxdlIE7sIgiAIgiCIAgeJWIIgCIIgCKLAUahErNFoxIwZM2A0GvPblEIPtYXrQG3hOlBbuBbUHq4DtYXr4EptoeOs5S8gCIIgCIIgCBejUHliCYIgCIIgiOcDErEEQRAEQRBEgYNELEEQBEEQBFHgIBFLEARBEARBFDgKjYhdu3YtypUrB09PT9SvXx9//PFHfpv03LNgwQK8+OKL8PPzQ7FixdCjRw9cuXJFtA3HcZg5cyaKFy8OLy8vtGrVCv/++28+WVx4WLBgAXQ6HcaNGyeso7bIW+7du4e33noLRYsWhbe3N+rUqYOzZ88Kn1N75A3Z2dmYNm0aypUrBy8vL5QvXx6zZ8+G2WwWtqG2cA7Hjh3DK6+8guLFi0On0+GHH34Qfa7lumdkZGD06NEIDg6Gj48PunXrhrt37+bht3g+UGuLrKwsTJ48GTVr1oSPjw+KFy+OAQMGIDY2VnSM/GiLQiFit2/fjnHjxmHq1KmIiIhA8+bN0alTJ8TExOS3ac81R48exahRo3Dy5EkcPHgQ2dnZaN++PVJTU4VtFi9ejGXLlmH16tU4ffo0wsLC0K5dOyQnJ+ej5c83p0+fxhdffIFatWqJ1lNb5B0JCQlo2rQp3N3d8csvv+DSpUtYunQpAgIChG2oPfKGRYsW4bPPPsPq1asRFRWFxYsXY8mSJVi1apWwDbWFc0hNTUXt2rWxevVq2c+1XPdx48Zh9+7d2LZtG44fP46UlBR07doVJpMpr77Gc4FaW6SlpeHcuXP4+OOPce7cOezatQtXr15Ft27dRNvlS1twhYCXXnqJGz58uGhd1apVuSlTpuSTRYWTuLg4DgB39OhRjuM4zmw2c2FhYdzChQuFbZ49e8YVKVKE++yzz/LLzOea5ORkrlKlStzBgwe5li1bcmPHjuU4jtoir5k8eTLXrFkzxc+pPfKOLl26cIMHDxat69mzJ/fWW29xHEdtkVcA4Hbv3i38reW6P336lHN3d+e2bdsmbHPv3j1Or9dz+/fvzzPbnzekbSHHqVOnOADc7du3OY7Lv7Z47j2xmZmZOHv2LNq3by9a3759e/z111/5ZFXhJDExEQAQFBQEAIiOjsaDBw9EbWM0GtGyZUtqGycxatQodOnSBW3bthWtp7bIW/bs2YMGDRrgtddeQ7FixVC3bl2sX79e+JzaI+9o1qwZDh06hKtXrwIA/vnnHxw/fhydO3cGQG2RX2i57mfPnkVWVpZom+LFi6NGjRrUNk4mMTEROp1OGD3Kr7Zwc9qRXYTHjx/DZDIhNDRUtD40NBQPHjzIJ6sKHxzHYcKECWjWrBlq1KgBAML1l2ub27dv57mNzzvbtm3D2bNncebMGYvPqC3ylps3b2LdunWYMGECPvroI5w6dQpjxoyB0WjEgAEDqD3ykMmTJyMxMRFVq1aFwWCAyWTCvHnz0LdvXwD028gvtFz3Bw8ewMPDA4GBgRbbUP/uPJ49e4YpU6bgzTffhL+/P4D8a4vnXsTy6HQ60d8cx1msI5zHe++9h/Pnz+P48eMWn1HbOJ87d+5g7NixOHDgADw9PRW3o7bIG8xmMxo0aID58+cDAOrWrYt///0X69atw4ABA4TtqD2cz/bt27Fp0yZs2bIF1atXR2RkJMaNG4fixYtj4MCBwnbUFvmDPded2sZ5ZGVl4Y033oDZbMbatWutbu/stnjuwwmCg4NhMBgs3gTi4uIs3vAI5zB69Gjs2bMHR44cQcmSJYX1YWFhAEBtkwecPXsWcXFxqF+/Ptzc3ODm5oajR49i5cqVcHNzE643tUXe8MILL6BatWqideHh4cJkU/pt5B2TJk3ClClT8MYbb6BmzZro378/xo8fjwULFgCgtsgvtFz3sLAwZGZmIiEhQXEbwnFkZWXh9ddfR3R0NA4ePCh4YYH8a4vnXsR6eHigfv36OHjwoGj9wYMH0aRJk3yyqnDAcRzee+897Nq1C4cPH0a5cuVEn5crVw5hYWGitsnMzMTRo0epbRxMmzZtcOHCBURGRgr/GjRogH79+iEyMhLly5entshDmjZtapFu7urVqyhTpgwA+m3kJWlpadDrxV2hwWAQUmxRW+QPWq57/fr14e7uLtrm/v37uHjxIrWNg+EF7LVr1/Dbb7+haNGios/zrS2cNmXMhdi2bRvn7u7OffXVV9ylS5e4cePGcT4+PtytW7fy27TnmhEjRnBFihThfv/9d+7+/fvCv7S0NGGbhQsXckWKFOF27drFXbhwgevbty/3wgsvcElJSfloeeGAzU7AcdQWecmpU6c4Nzc3bt68edy1a9e4zZs3c97e3tymTZuEbag98oaBAwdyJUqU4H766ScuOjqa27VrFxccHMx98MEHwjbUFs4hOTmZi4iI4CIiIjgA3LJly7iIiAhhxruW6z58+HCuZMmS3G+//cadO3eOe/nll7natWtz2dnZ+fW1CiRqbZGVlcV169aNK1myJBcZGSnqzzMyMoRj5EdbFAoRy3Ect2bNGq5MmTKch4cHV69ePSHNE+E8AMj+27hxo7CN2WzmZsyYwYWFhXFGo5Fr0aIFd+HChfwzuhAhFbHUFnnL3r17uRo1anBGo5GrWrUq98UXX4g+p/bIG5KSkrixY8dypUuX5jw9Pbny5ctzU6dOFXXO1BbO4ciRI7J9xMCBAzmO03bd09PTuffee48LCgrivLy8uK5du3IxMTH58G0KNmptER0drdifHzlyRDhGfrSFjuM4znl+XoIgCIIgCIJwPM99TCxBEARBEATx/EEiliAIgiAIgihwkIglCIIgCIIgChwkYgmCIAiCIIgCB4lYgiAIgiAIosBBIpYgCIIgCIIocJCIJQiCIAiCIAocJGIJgiAIgiCIAgeJWIIgiHxg5syZqFOnTr6d/+OPP8awYcM0bTtx4kSMGTPGyRYRBEHYBlXsIgiCcDA6nU7184EDB2L16tXIyMhA0aJF88iq/3j48CEqVaqE8+fPo2zZsla3j4uLQ4UKFXD+/HmUK1fO+QYSBEFogEQsQRCEg3nw4IGwvH37dkyfPh1XrlwR1nl5eaFIkSL5YRoAYP78+Th69Ch+/fVXzfv06tULFStWxKJFi5xoGUEQhHYonIAgCMLBhIWFCf+KFCkCnU5nsU4aTjBo0CD06NED8+fPR2hoKAICAjBr1ixkZ2dj0qRJCAoKQsmSJbFhwwbRue7du4c+ffogMDAQRYsWRffu3XHr1i1V+7Zt24Zu3bqJ1v3vf/9DzZo14eXlhaJFi6Jt27ZITU0VPu/WrRu2bt2a62tDEAThKEjEEgRBuAiHDx9GbGwsjh07hmXLlmHmzJno2rUrAgMD8ffff2P48OEYPnw47ty5AwBIS0tD69at4evri2PHjuH48ePw9fVFx44dkZmZKXuOhIQEXLx4EQ0aNBDW3b9/H3379sXgwYMRFRWF33//HT179gQ7UPfSSy/hzp07uH37tnMvAkEQhEZIxBIEQbgIQUFBWLlyJapUqYLBgwejSpUqSEtLw0cffYRKlSrhww8/hIeHB/78808AOR5VvV6PL7/8EjVr1kR4eDg2btyImJgY/P7777LnuH37NjiOQ/HixYV19+/fR3Z2Nnr27ImyZcuiZs2aGDlyJHx9fYVtSpQoAQBWvbwEQRB5hVt+G0AQBEHkUL16dej1//kWQkNDUaNGDeFvg8GAokWLIi4uDgBw9uxZXL9+HX5+fqLjPHv2DDdu3JA9R3p6OgDA09NTWFe7dm20adMGNWvWRIcOHdC+fXv07t0bgYGBwjZeXl4Acry/BEEQrgCJWIIgCBfB3d1d9LdOp5NdZzabAQBmsxn169fH5s2bLY4VEhIie47g4GAAOWEF/DYGgwEHDx7EX3/9hQMHDmDVqlWYOnUq/v77byEbQXx8vOpxCYIg8hoKJyAIgiig1KtXD9euXUOxYsVQsWJF0T+l7AcVKlSAv78/Ll26JFqv0+nQtGlTzJo1CxEREfDw8MDu3buFzy9evAh3d3dUr17dqd+JIAhCKyRiCYIgCij9+vVDcHAwunfvjj/++APR0dE4evQoxo4di7t378ruo9fr0bZtWxw/flxY9/fff2P+/Pk4c+YMYmJisGvXLjx69Ajh4eHCNn/88QeaN28uhBUQBEHkNyRiCYIgCije3t44duwYSpcujZ49eyI8PByDBw9Geno6/P39FfcbNmwYtm3bJoQl+Pv749ixY+jcuTMqV66MadOmYenSpejUqZOwz9atW/HOO+84/TsRBEFohYodEARBFDI4jkOjRo0wbtw49O3b1+r2P//8MyZNmoTz58/DzY2mUhAE4RqQJ5YgCKKQodPp8MUXXyA7O1vT9qmpqdi4cSMJWIIgXAryxBIEQRAEQRAFDvLEEgRBEARBEAUOErEEQRAEQRBEgYNELEEQBEEQBFHgIBFLEARBEARBFDhIxBIEQRAEQRAFDhKxBEEQBEEQRIGDRCxBEARBEARR4CARSxAEQRAEQRQ4SMQSBEEQBEEQBY7/A08TE5CL7NW9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,2))\n",
    "ax.plot(gmtime,accel,'k')\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Acceleration (m/s2)')\n",
    "ax.set_xlim([0,gmtime[-1]])\n",
    "plt.grid(which='both',color='lightgray')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d802e3b8",
   "metadata": {},
   "source": [
    "# Compute ground motion intensity measures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d8591bd",
   "metadata": {},
   "source": [
    "## PGA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "6e70edd4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGA = 1.286 m/s2\n"
     ]
    }
   ],
   "source": [
    "PGA = gmimtools.get_PGA_cy(accel)\n",
    "print('PGA = {:.3f} m/s2'.format(PGA))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77beffe9",
   "metadata": {},
   "source": [
    "## PGV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c7ec3216",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PGV = 0.336 m/s\n"
     ]
    }
   ],
   "source": [
    "PGV = gmimtools.get_PGV_cy(accel,dt)\n",
    "print('PGV = {:.3f} m/s'.format(PGV))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cdf2c7c",
   "metadata": {},
   "source": [
    "## Ia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6eb52485",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ia = 0.717 m/s\n"
     ]
    }
   ],
   "source": [
    "[Ia,Ia_hist] = gmimtools.get_Ia_cy(accel,dt)\n",
    "print('Ia = {:.3f} m/s'.format(Ia))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24a5a351",
   "metadata": {},
   "source": [
    "## Ia times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c03718aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "percentages = np.array([0.05,0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "93206c7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "Ia_times = gmimtools.get_Ia_times_cy(accel,dt,percentages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4f15a320",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAD/CAYAAADi4l4lAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAnUlEQVR4nO3deVxU9f4/8Newy74PmwgqKgqCiZWiGaIs6jUrb179lpZlGZopZUpaLqVo3ZTrmlbaT9O0xbpUphIqrqUioAaJC4jKIILsyDZzfn8gc0O2YZwVXs/HYx4yn/M5Z94z74g3h3M+b5EgCAKIiIiIiPSQgbYDICIiIiJSFotZIiIiItJbLGaJiIiISG+xmCUiIiIivcViloiIiIj0FotZIiIiItJbLGaJiIiISG8ZaTsATZPJZMjNzYWVlRVEIpG2wyEiIiKiBwiCgLKyMri5ucHAoPVzr52umM3NzUXXrl21HQYRERERteHGjRvw8PBodU6nK2atrKwAANnZ2bCzs9NyNCSVSnHp0iX07t0bhoaG2g6nU2MudAdzoTuYC93BXOgWdeejtLQUXbt2lddtrel0xWzDpQXW1tawtrbWcjQklUphaWkJa2tr/s9Jy5gL3cFc6A7mQncwF7pFU/lQ5JJQ3gBGRERERHqLxSwRERER6S0Ws0RERESkt1jMEhEREZHeYjFLRERERHqr061mQERERNTZCYKA2tpalJeXo6CgAMXFxbh79y6ysrJQWVkJmUwGqVTa4r91dXW4c+cObG1tIQiCfLylfWQyWbviq62tVXgui1kiIiIiPdFQgJaVleHOnTu4e/cuqqurUVVVhaqqKmRnZ+PWrVsoLi5GeXk5KioqUFlZicrKyiZfS6VSbb8dlWAxS0RERKRhMpkMd+7cQVFREWpqauSP6upq1NTUoLy8HHfu3EFeXh4yMjLwxx9/IC8vD9XV1SqPxcrKCnZ2drC1tUW3bt1gY2MDQ0NDGBgYtPivSCRCcXExnJycYGRk1Op8Q0NDiEQihdaMbVBVVYX58+crNJfFLBEREZGKlZWV4fvvv0d+fr78jGhFRQVKS0vx119/4fz58+36U/qDTExM4ODgAEdHR5iZmcHMzAympqZwdXVF9+7dYWtrC0tLS1hYWMDCwgLm5ubyfx/82sTEpN2vL5VKkZGRAV9fX7V1AGMxS0RERKRB9+7dQ2pqKuLj47FlyxbcvXu31fkikQi2trYwNTWFiYlJo4e5uTmcnJzg7OyMXr16ITAwED4+PrCysoKlpaVSBWhHxWKWiIiISEm1tbXYsGED4uLikJOTA0EQ5Nt69OiB4OBg+dlRCwsLWFpaolu3bggKCoKHhweMjY21GH3HwGKWiIiIqAWCICA1NRUSiQRubm7o27cvsrOzUVBQgMuXL2PlypX466+/5PPFYjEGDhyIV199FWPHjlXLn+CpMRazRERERM1ITEzE7NnRSE8/Lx8zNDSFVNr4JixHR0d88MEHePrpp+Hs7NyuG53o4bGYJSIiInpAYmIiwsMjIJMNAfArAH8AFyCVLgdwAq6uLvD09ER4eDjmzp0LW1tbrcbbmbGYJSIiIvobQRAQFTUbMtlgCMJvABqua3UHEAoDg5GwsyvGqVOneBZWB7CYJSIiok5HJpPhxo0bKCsrQ2FhIRISEpCSkoKSkhLk5eXh6tWrqD8j++ANWsaQyWKQnh6JtLQ0BAYGaj54aoTFLBEREXVYFRUVOHbsGG7cuIFr164hLy8PmZmZOH/+PMrLy9vY27/V8dzcXBazOoDFLBEREek9QRBQWFiIvLw8SCQSeeesLVu2oLCwsNl9TExMYGNjAysrKwwYMACjRo2Co6Mj8vPzERUVBeAC6i8teNAFAICbm5va3g8pjsUsERER6bXFixdj1apVLbZ69fT0hJ+fH7y9veHm5gYvLy8EBASgV69eza7zKggC1q//FBkZsRCEUDS+1KAWBgax8PUNQEBAgHreELULi1kiIiLSS1KpFB988AGWLVsmH3NwcICLiwtcXV3h6uqK0aNHY8KECTAyUrzkEYlEWLduDcLCwiEIIyGTxaBhNQMDg1iIRCexdu0B3vylI1jMEhERkUo92GggICBA6cJPEARkZ2fj3Llz+P3335GdnY27d++ipKQEly9fRmlpKQBg+fLleOutt2BqaqqS9zBixAgcOLD//jqzkfJxX98ArF17ACNGjFDJ69DDYzFLREREKtNco4G+fftj3bo1ChWAdXV1OHv2LP773/9i9+7dyMvLQ1VVVYvz7ezs8O9//xsvvfSSys+UhoaG4uLFVKSlpSE3N/ehC3NSDxazREREpBItNRrIyIhFWFg4DhzYj0ceeQS3bt1CQUEBTp8+jQMHDiA3Nxfl5eUoKytDeXk5pFJpo+MaGxvD398fDg4OePLJJ+Hl5QUbGxu4u7vDz8+vXZcQtJdIJEJgYCBXLdBhLGaJiIjooQmCgNmzoyGTDWnSaEAQQiGVjkBk5D9QW3uvzWPZ29tj8ODBeOGFFzBo0CC4u7vDyMgIGRkZ8PX1haGhoVrfC+kXFrNERESktHPnziE5ORlJSUn3Ly1ovtEAsBC1tfXXnjo6OsLBwQFeXl546qmn4OvrC0tLS1hZWcHKygpisbhJwfrg2VqiBixmiYiIqN0KCgrw2muvYe/evQ9sab3RwO7duzFx4kS1xkadC4tZIiIiUkhxcTGSkpLw119/Yd26dbh16xYMDQ0REhICZ2dn7Nq1C201Gujdu7cmQ6ZOgMUsERERtUgmkyE5ORnffPMNvvzySxQUFMi39erVC9988w0CAgLuL8d1kY0GSOPaXcwKgoCkpCQcO3YM2dnZqKyshJOTEwYMGICRI0eia9eu6oiTiIiINKy4uBghISFITU2Vj3l6emLw4MEYMGAAXn/9dVhbWwNgowHSHgNFJ967dw8rVqxA165dERkZiV9++QXFxcUwNDTElStXsHjxYnh7e2P06NH4/fff1RkzERERKUEQBKSkpGDfvn1ITU2FIAgtzr137x7+9a9/yQvZ5557Dt9//z0uX76M3bt3Y/78+fJCtkFDo4E+fYoBRALwABAJX98SHDzIRgOkHgqfme3Vqxcee+wxfPrppwgPD2+2l/H169exa9cuTJw4EYsWLcL06dNVGiwREREpp7VmBo6OjtiyZQtu3ryJiooKVFRU4NKlS7h79y66dOmCpKQkDBo0SKHXYaMB0jSFi9lff/0Vfn5+rc7p1q0bYmJi8NZbb+H69esPHRwRERE9vNaaGYwaFQ4jIwPU1NQ02a9bt27Ytm2bwoVsAzYaIE1SuJhtq5D9OxMTE/j4+CgVEBEREamOIAiYOfPNFpsZCMII1NScwhNPPIHJkyfD0tISFhYWcHJywuOPP84GBaTzlFrNYP/+/bC0tMTQoUMBABs2bMBnn32Gvn37YsOGDbCzs1P4WBs3bsTHH38MiUSCfv36IS4uDsOGDWtxfnV1NZYtW4avvvoKeXl58PDwwMKFCzFt2jRl3goREZHeEwQBOTk5uHLlCk6dOoWSkhKUlZWhrKwMJ06cuP/X0pabGQCRWLlyJQYPHqzx2IkellLF7Lx587Bq1SoAwIULF/DWW28hOjoahw4dQnR0NLZt26bQcfbs2YM5c+Zg48aNCA4OxubNmxEZGYn09HR4eno2u89zzz2H27dv44svvkDPnj2Rn5+Puro6Zd4GERGR3qirq0NRURFKSkqQm5uLy5cvo6SkBLdu3cL333+vwOV9rTczKCoqUmm8RJqiVDGblZWFvn37AgC+//57jB07FitWrMC5c+cwevRohY+zevVqvPzyy3jllVcAAHFxcThw4AA2bdqE2NjYJvP379+PpKQkXLt2Dfb29gAALy+vVl+juroa1dXV8uelpaUA6tvisTWe9jXkgLnQPuZCdzAXukMbuaitrUVmZiZOnDiBH3/8EcXFxSgrK8P169dRWVnZ4n7GxsYQi8UYMmQI3N3d5a1hpVIpFixYgLaaGYjFYp3+b47fF7pF3floz3GVKmZNTEzk31C//fYbpkyZAgCwt7eXF4ttqampQXJy8v1vsP8JCwvDyZMnm90nPj4eQUFB+Oijj7Bjxw5YWFhg3Lhx+OCDD9ClS5dm94mNjcXSpUubjF+5cgWWlpYKxUrql5mZqe0Q6D7mQncwF7pD1bmorKzErVu3kJ2djYyMDOTm5qKiogISiQTXrl1DbW1ti/uam5vD1tYW3t7esLW1haWlJQYNGoThw4fDzMysyXxBEPDZZ9tw7VrzzQxEolh0794bJiYmyMjIUOn7VAd+X+gWdeWjvLxc4blKFbNDhw5FdHQ0goODcfr0aezZswdA/Rvy8PBQ6BgFBQWQSqUQi8WNxsViMfLy8prd59q1azh+/DjMzMzwww8/oKCgAFFRUbh79y62bt3a7D4xMTGIjo6WPy8tLUXXrl3Rs2fPdl3bS+ohlUqRmZmJXr168SYDLWMudAdzoTtUnQtBEDB//nysX7++2dUDGlhZWaFfv34YNWoUHnnkEVhYWMDd3R09evSAkVH7f3Rv2rQOkZGjW2xmsGnTPvlfXHUVvy90i7rzoejJUUDJYnb9+vWIiorCd999h02bNsHdvf7PFr/++isiIiLadawH150TBKHFtehkMhlEIhF27twJGxsbAPWXKkyYMAEbNmxo9uysqakpTE1Nm4wbGhrym0GHMB+6g7nQHcyF7mguF/XtW1MhkUhaXEtVEAQUFRXh+vXruHbtGn755Rf5fSV2dnbw8vLCwIED0bt3b9jY2EAsFqN///7o1q2bStdlHTVqFA4c2H9/ndlI+bivbwDWrtWvZgb8vtAt6spHe47ZrmL24MGDCAkJgaenJ37++ecm29esWaPwsRwdHWFoaNjkLGx+fn6Ts7UNXF1d4e7uLi9kAcDX1xeCIODmzZtcDoyIiDSitQYEI0aMQElJCeLi4rB+/XoUFBQ02X/dunWYOXOmRhsJsJkBdVQKt7MFgBkzZsDJyQkTJ07E119/jZKSEqVf2MTEBAMHDkRCQkKj8YSEBAwZMqTZfYKDg5Gbm9voOorMzEwYGBgofHkDERHRw2hoQJCRYYv65a5uAvgVGRm2CAsLx4EDBzBq1CgsWbJEXsiKxWI89thjeP755/Htt99i1qxZWikiG5oZjB49GoGBgSxkqUNoVzF77do1HD16FP7+/lizZg3EYjFCQ0Oxdu1aZGdnt/vFo6Oj8fnnn2Pr1q3IyMjA3LlzkZOTgxkzZgCov9614eYyAJg8eTIcHBzw0ksvIT09HUePHsW8efMwbdq0Fm8AIyIiUhVBEDB7dvTfGhBEoH6FgAgIwm+QSh9HRMRYnDlzBkD9GdiKigrk5eXh999/x44dOzBhwgRtvgWiDqfd18z2798f/fv3x6JFi5Cbm4v4+HjEx8dj/vz56NWrF5566imMGzcOQUFBbR5r4sSJKCwsxLJlyyCRSODn54d9+/ahW7duAACJRIKcnBz5fEtLSyQkJOCNN95AUFAQHBwc8Nxzz+HDDz9s79sgIiJqt9TU1PuXFrTegMDKygqbN2/GpEmTNB4jUWej1A1gDdzc3DBjxgzMmDEDFRUV+PXXXxEfH4+IiAhER0fj3XffbfMYUVFRiIqKanbbl19+2WSsT58+TS5NICIi0gSJRHL/q9YbEOzatQtjx47VSExEnd1DFbN/Z2FhgQkTJmDChAmQyWQoLCxU1aGJiIg0SiaToaqqCsXFxcjOzkZhYSHu3Lkjv3ygrQYEvI+DSHOULmZPnz6NI0eOID8/HzKZTD4uEonwySefwMnJSSUBEhERaVJmZiZCQ0Nx8+bNFmYYAlgOoGkDAgODWPj6BiAgIEDtcRJRPaWK2RUrVmDRokXo3bs3xGJxo7sheWckERHps1WrVskLWZFIBAsLCzg7O8PV1RUuLi6QSqWIj/8ZQPMNCNauPcCfhUQapFQx+5///Adbt27Fiy++qOJwiIiI1KetRge3b9/Gjh07AEB+D8ilS5fg6+vbaBH3/60zq98NCIg6gnYtzSXfycAAwcHBqo6FiIhIbRITE+HnF4hHHnkEY8aMwYABA+DnF4hDhw4BAD788EO4uLigtrYWFhYWCAsLg4FB8z8mGxoQpKSk4JdffkFKSgouXEhhIUukBUoVs3PnzsWGDRtUHQsREZFatNXoYMeOHVi6dCkAwNvbG1u3bm22FfrfsQEBkW5Q6jKDt99+G2PGjEGPHj3Qt29fGBs3Xmtv7969KgmOiIjoYTVtdNDwM8sdghAKqTQUL7/8Gurq6hAeHo79+/drM1wiaielitk33ngDhw8fRkhICBwcHPjbKBER6ay2Gx28i9raSBgbG+ODDz7QfIBE9FCUKma3b9+O77//HmPGjFF1PERERCqlaKODDRs2YNCgQRqJiYhUR6lrZu3t7dGjRw9Vx0JERKRyrq6u97+60MKM+nEWskT6SalidsmSJVi8eDEqKytVHQ8REZFKnDlzBsuXL0d0dDT+1+ig9oFZ9Y0O+vVjowMifaXUZQZr167F1atXIRaL4eXl1eQGsHPnzqkkOCIiovZKTU3F/PnzcfDgwQe2nISBARsdEHU0ShWz48ePV3EYRERE7fNgA4QePXpg6dKlWLduHWpqamBkZITw8HD84x//QEBAAEpLSzF37jw2OiDqYJQqZhcvXqzqOIiIiBT2vw5c5+Vj1tb2KC29C6C+qcGaNWvg79/4pq+LF0chLS0Nubm5zXYAIyL9o1QxqwhBEPg/CCIiUrmGBggy2RDUL7dVf8lAaelyACcQE7MAH374YbPduxoaHQQGBmo2aCJSG4VvAPP19cWuXbtQU1PT6rzLly/j9ddfx6pVqx46OCIior9r2gAhAoD7/X8PAQjGf//7C0+mEHUiCp+Z3bBhA+bPn4+ZM2ciLCwMQUFBcHNzg5mZGYqKipCeno7jx48jPT0ds2bNQlRUlDrjJiKiTqjtBggLkZ4eibS0NJ59JeokFC5mR4wYgTNnzuDkyZPYs2cPdu3ahezsbNy7dw+Ojo4YMGAApkyZgueffx62trZqDJmIiDorRRsg5Obmspgl6iTafc3skCFDMGTIEHXEQkRE1Cpra+v7X11A/eUFD6pvgODm5qapkIhIy5RqmkBERKRpdXV1mDJlCtgAgYj+jsUsERHphRMnTiArKwuAFCLRSYhEIwHsB3ALwH4YGIy83wBhNW8AI+pEWMwSEZFeiI+PBwBMmTIFCQkH4etbDCASgAeASPj6luDgQTZAIOps1LbOLBER0cP4e4cvFxcXfP/99wCAcePGITQ0FBcvprIBAhGxmCUiIt3TXIcvwBAWFhaIjKxvR8sGCEQEqOAyg3v37qG0tLTRg4iISFkNHb4yMmxRv57szfv/DkZlZRVOnTql1fiISLcoVcxWVlZi1qxZcHZ2hqWlJezs7Bo9iIiIlNFWhy+RKBizZ0dDEATtBkpEOkOpYnbevHk4dOgQNm7cCFNTU3z++edYunQp3NzcsH37dlXHSEREnURDhy9BiEFzHb5kshikp59HWlqaNsIjIh2kVDH7008/YePGjZgwYQKMjIwwbNgwLFq0CCtWrMDOnTtVHSMREXUS7enwRUQEKFnM3r17F97e3gDqu7HcvXsXADB06FAcPXpUddEREVGn4urqev+rCy3MYIcvImpMqWK2e/fuyM7OBgD07dsX33zzDYD6M7a2traqio2IiDqZHj16wNzcGuzwRUSKUqqYfemll+TXK8XExMivnZ07dy7mzZun0gCJiKjzWLZsGSorSwGcgEgUCnb4IqK2KLXO7Ny5c+Vfh4SE4K+//sLZs2fRo0cP/rZMREQK+3tjBDc3N3z33XcAgPfffw/fffcj0tMj5XN9fQOwdi07fBFRYyppmuDp6QlPT09VHIqIiDqJlhojGBgYYN68eViyZAk7fBFRm9pVzK5du1ahebNnz1YqGCIi6hwaGiPIZENQ3xDBH/U3dy2HTHYSf/zxB0JDQ9nhi4ja1K5ids2aNW3OEYlELGaJiKhFTRsjNKwn6w4gFCJRKGbPjsbFi6k8E0tEbWpXMZuVlaWuOIiIqJNoaIxQf0a2aWMEQXgX6emRSEtL41lZImqTUqsZEBERKWvPnj33v2JjBCJ6eCxmiYhIY3JycrB69er7z9gYgYgeHotZIiLSiK+++gphYWGora1Fly5WEIliwcYIRPSwWMwSEZHa/fnnn3jhhRdw6dIlmJmZYf36OBgYnISBwUiwMQIRPQytF7MbN26Et7c3zMzMMHDgQBw7dkyh/U6cOAEjIyPeHEBEpKMEQUBKSgr27dsnv7TA0tISp0+fxrRp03DgwH706VMMIBKAB4BI+PqW4OBBNkYgIsUp3TTh6tWr2LZtG65evYr//Oc/cHZ2xv79+9G1a1f069dPoWPs2bMHc+bMwcaNGxEcHIzNmzcjMjIS6enprTZhKCkpwZQpUxAaGorbt28r+xaIiEhNWmqI8MYbb8Dfv/4Gr9DQUFy8mMrGCET0UJQqZpOSkhAZGYng4GAcPXoUy5cvh7OzM86fP4/PP/9c3o6wLatXr8bLL7+MV155BQAQFxeHAwcOYNOmTYiNjW1xv9deew2TJ0+GoaEhfvzxx1Zfo7q6GtXV1fLnpaWlAACpVAqpVKpQnKQ+DTlgLrSPudAd+p6LxMREjB49ptmGCB999DGefPJJhIaGyuf7+/vLC1yZTKaNkFuk77noSJgL3aLufLTnuCJBEIT2vsDgwYPxz3/+E9HR0bCyskJaWhq6d++OM2fOYPz48bh161abx6ipqYG5uTm+/fZbPP300/LxN998E6mpqUhKSmp2v23btmHjxo04deoUPvzwQ/z4449ITU1t8XWWLFmCpUuXNhk/deoULC0t236zRESkMEEQ8PTT/8S1a+IHGiIAQC1EopHo3v02fvjhW56BJaIWlZeXY/DgwSgpKYG1tXWrc5U6M3vhwgXs2rWrybiTkxMKCwsVOkZBQQGkUinEYnGjcbFYjLy8vGb3uXz5MhYsWIBjx47ByEix0GNiYhAdHS1/Xlpaiq5du6Jnz56ws7NT6BikPlKpFJmZmejVqxcMDQ21HU6nxlzoDn3ORUpKCq5evQQgDs03RIjB1auRqKmp0Yt7HvQ5Fx0Nc6Fb1J2Phr+kK0KpYtbW1hYSiQTe3t6NxlNSUuDu7t6uYz34m7kgCM3+ti6VSjF58mQsXboUvXr1Uvj4pqamMDU1bTJuaGjIbwYdwnzoDuZCd+hjLvLz8+9/1XpDhNu3b+vVe9PHXHRUzIVuUVc+2nNMpVYzmDx5MubPn4+8vDyIRCLIZDKcOHECb7/9NqZMmaLQMRwdHWFoaNjkLGx+fn6Ts7UAUFZWhrNnz2LWrFkwMjKCkZERli1bhrS0NBgZGeHQoUPKvBUiIlIhV1fX+1+xIQIRaYZSxezy5cvh6ekJd3d3lJeXo2/fvnjiiScwZMgQLFq0SKFjmJiYYODAgUhISGg0npCQgCFDhjSZb21tjQsXLiA1NVX+mDFjBnr37o3U1FQ89thjyrwVIiJSocDAQPTt258NEYhIY5S6zMDY2Bg7d+7EsmXLkJKSAplMhgEDBsDHx6ddx4mOjsYLL7yAoKAgDB48GFu2bEFOTg5mzJgBoP5611u3bmH79u0wMDCAn59fo/2dnZ1hZmbWZJyIiLRDJBJh3bo1GDkyDMAIAAvRsJqBgUHs/YYIB3jzFxGpjNJLcw0fPhw9evRAjx49lH7xiRMnorCwEMuWLYNEIoGfnx/27duHbt26AQAkEglycnKUPj4REWleSEgIfH17Iz39FOobItTz9Q3A2rVsiEBEqqXUZQajRo2Cp6cnFixYgIsXLz5UAFFRUcjOzkZ1dTWSk5PxxBNPyLd9+eWXOHLkSIv7LlmypNVluYiISLMSExPh5xeI9PR0APXrRHbt6oXNmzfjwoUUFrJEpHJKFbO5ubl45513cOzYMfTv3x/9+/fHRx99hJs3b6o6PiIi0hOJiYkID49ARoYt6psl3ATwK27e9ERU1EzeqEtEaqFUMevo6IhZs2bhxIkTuHr1KiZOnIjt27fDy8uLv3UTEXVCgiBg9uxoyGRD7jdLiADgDiACgvAbBGEIZs+OhhJ9eoiIWqVUMft33t7eWLBgAVauXAl/f/8WO3cREVHHlZqaivT08xCEGDTXLEEmi0F6+nmkpaVpIzwi6sAeqpg9ceIEoqKi4OrqismTJ6Nfv374+eefVRUbERHpCYlEcv+r1psl5ObmaiQeIuo8lFrN4N1338XXX3+N3NxcjBw5EnFxcRg/fjzMzc1VHR8REemBxs0SmusEyWYJRKQeShWzR44cwdtvv42JEyfC0dFR1TEREZGeaWiWkJERC0EIReNLDeqbJfj6slkCEameUsXsyZMnVR0HERHpsYZmCaNGhUMQ2CyBiDRH4WI2Pj4ekZGRMDY2Rnx8fKtzx40b99CBERGRfgkJCcELL/wf/t//+wpslkBEmqJwMTt+/Hjk5eXB2dkZ48ePb3GeSCSCVCpVRWxERKQnEhMTMXt2NNLTz8vHunb1wqJFMZg+fTrPyBKR2ii8moFMJoOzs7P865YeLGSJiDoXNksgIm1Sammu7du3o7q6usl4TU0Ntm/f/tBBERGRfmCzBCLSNqWK2ZdeegklJSVNxsvKyvDSSy89dFBERKQf2CyBiLRNqWJWEIRmr3+6efMmbGxsHjooIiLSD2yWQETa1q6luQYMGACRSASRSITQ0FAYGf1vd6lUiqysLERERKg8SCIi0k1slkBE2tauYrZhFYPU1FSEh4fD0tJSvs3ExAReXl549tlnVRogERHpLjZLICJta1cxu3jxYgCAl5cXJk6cCDMzM7UERURE+qGhWUJYWDikUjZLICLNU6oD2NSpU1UdBxER6QlBEJCamgqJRAI3NzeEhITgp5/iMWbMUxAENksgIs1SqpiVSqVYs2YNvvnmG+Tk5KCmpqbR9rt376okOCIi0i3NNUfo27c/IiJGQhBqYW9vj+3bt8Pd3R0BAQE8I0tEaqfUagZLly7F6tWr8dxzz6GkpATR0dF45plnYGBggCVLlqg4RCIi0gUtNUfIyLDF6tVrAADTp0/HmDFjEBgYyEKWiDRCqWJ2586d+Oyzz/D222/DyMgIkyZNwueff473338fv//+u6pjJCIiLWurOQIQDJHIGO+99552AyWiTkepYjYvLw/+/vVrB1paWsobKIwdOxa//PKL6qIjIiKd0FZzBGAhBKEWly9f1kJ0RNSZKVXMenh4yBfK7tmzJw4ePAgAOHPmDExNTVUXHRER6QQ2RyAiXaVUMfv0008jMTERAPDmm2/ivffeg4+PD6ZMmYJp06apNEAiItK+xs0RmsPmCESkHUqtZrBy5Ur51xMmTICHhwdOnjyJnj17Yty4cSoLjoiIdAObIxCRrlKqmH3Q448/jscff1wVhyIiIh3E5ghEpKsULmbj4+MVPijPzhIRdTwjRozAgQP7MW7cM6isZHMEItINChez48ePV2ieSCSCVCpVNh4iItJhI0aMgLOzPbKzSzFz5ky8/PLLXFOWiLRK4RvAZDKZQg8WskREHVNiYiL8/AKRnZ0NANiwYQOef/5FHD58WLuBEVGnptRqBkRE1Lm01v0rLCxcvsINEZGmKXUD2LJly1rd/v777ysVDBER6Z6m3b8aVjJwhyCEQhBGYvbsaFy8mMrLDYhI45QqZn/44YdGz2tra5GVlQUjIyP06NGDxSwRUQfS0P2r/oxs0+5fMlkM0tMjkZaWhsDAQM0HSESdmlLFbEpKSpOx0tJSvPjii3j66acfOigiItId7en+xWKWiDRNZdfMWltbY9myZXjvvfdUdUgiItIB7P5FRLpMpTeAFRcXo6SkRJWHJCIiLWvo/iUSxQKofWBrffevfv3Y/YuItEOpywzWrl3b6LkgCJBIJNixYwciIiJUEhgREekGdv8iIl2mVDG7Zs2aRs8NDAzg5OSEqVOnIiYmRiWBERGR7ggJCcGcOW/ik0/iALD7FxHpDqWK2aysLFXHQUREOioxMRGzZ0ffX9GgXteuXli0KAbTp0/nGVki0io2TSAioha11Czh5k1PREXNxKFDh7QbIBF1ekqdma2qqsK6detw+PBh5OfnQyaTNdp+7tw5lQRHRETaw2YJRKQPlCpmp02bhoSEBEyYMAGPPvoo/ydGRNQBsVkCEekDpYrZX375Bfv27UNwcPBDB7Bx40Z8/PHHkEgk6NevH+Li4jBs2LBm5+7duxebNm1Camoqqqur0a9fPyxZsgTh4eEPHQcRETXGZglEpA+UumbW3d0dVlZWD/3ie/bswZw5c7Bw4UKkpKRg2LBhiIyMRE5OTrPzjx49ilGjRmHfvn1ITk5GSEgI/vGPfzTbkYyIiB4OmyUQkT5Qqpj95JNPMH/+fFy/fv2hXnz16tV4+eWX8corr8DX1xdxcXHo2rUrNm3a1Oz8uLg4vPPOOxg0aBB8fHywYsUK+Pj44KeffnqoOIiIqCk2SyAifaDUZQZBQUGoqqpC9+7dYW5uDmPjxtdS3b17t81j1NTUIDk5GQsWLGg0HhYWhpMnTyoUh0wmQ1lZGezt7VucU11djerqavnz0tJSAIBUKoVUKlXodUh9GnLAXGgfc6E7dCkXcXH/RkTEaAhC880S1qzZ1+Qm4I5El3LR2TEXukXd+WjPcZUqZidNmoRbt25hxYoVEIvFSt0AVlBQAKlUCrFY3GhcLBYjLy9PoWN88sknqKiowHPPPdfinNjYWCxdurTJ+JUrV2Bpadm+oEltMjMztR0C3cdc6A5dyIWLiwvGjx+HvXv/i783S/D27o2YmE/h6uqKjIwM7QWoIbqQC6rHXOgWdeWjvLxc4blKFbMnT57EqVOnVPKnpQcLYUEQFCqOv/76ayxZsgT//e9/4ezs3OK8mJgYREdHy5+Xlpaia9eu6NmzJ+zs7JQPnFRCKpUiMzMTvXr1gqGhobbD6dSYC92hK7lITEzE3LnzmjRLePfd+XjllVc6xUo2upILYi50jbrz0fCXdEUoVcz26dMH9+7dU2ZXOUdHRxgaGjY5C5ufn9/kbO2D9uzZg5dffhnffvstRo4c2epcU1NTmJqaNhk3NDTkN4MOYT50B3OhO7SZi8TERIwePQYy2RDUL81Vf3nBzZuxmDXrDfj4+CA0NFQrsWkDvy90B3OhW9SVj/YcU6kbwFauXIm33noLR44cQWFhIUpLSxs9FGFiYoKBAwciISGh0XhCQgKGDBnS4n5ff/01XnzxRezatQtjxoxRJnwiImpF02YJEQDcAURAEH6DIAzB7NnREARBy5ESESl5ZjYiIgIAmvxW3nCJgKIX7UZHR+OFF15AUFAQBg8ejC1btiAnJwczZswAUH+JwK1bt7B9+3YA9YXslClT8J///AePP/64/Kxuly5dYGNjo8xbISKiB7BZAhHpE6WK2cOHD6vkxSdOnIjCwkIsW7YMEokEfn5+2LdvH7p16wagfsHuv685u3nzZtTV1WHmzJmYOXOmfHzq1Kn48ssvVRITEVFnx2YJRKRPlCpmhw8frrIAoqKiEBUV1ey2BwvUI0eOqOx1iYioeY2bJbg3M4PNEohIdyhVzB49erTV7U888YRSwRARkfY1NEtIT18OIBSNLzWob5bg68tmCUSkG5QqZp988skmY39fooULGhMR6S+RSITY2A/w1FPjATTfLGHt2gOdYmkuItJ9Sq1mUFRU1OiRn5+P/fv3Y9CgQTh48KCqYyQiIjURBAEpKSnYt28fUlNT5SsUGBgYABBgYnIO9c0SPABEwte3BAcPHsCIESO0GDUR0f8odWa2uZUDRo0aBVNTU8ydOxfJyckPHRgREalXYmIiZs+ObtQUoW/f/pg+/SVs2rQJAPDii8/j9ddfR25uLtzc3BAQEMAzskSkU5QqZlvi5OSES5cuqfKQRESkBomJiQgPj2jSFCE9fQXmzo0GIMDc3Bxz5syBr68vVy0gIp2lVDF7/vz5Rs8FQYBEIsHKlSt5QwARkY5r2hSh4QYvd9Tf8DUCNjYZOHPmFHx8fLQXKBGRApQqZgMDAyESiZp0f3n88cexdetWlQRGRETq0VZTBGAhSkoiUVFRofngiIjaSaliNisrq9FzAwMDODk5wczMTCVBERGR+rApAhF1JEoVsw0duoiISP+wKQIRdSTtWprr0KFD6Nu3L0pLS5tsKykpQb9+/XDs2DGVBUdERKoXGBiInj37AFgOoPaBrfVNEfr1Y1MEItIP7Spm4+LiMH36dFhbWzfZZmNjg9deew2rV69WWXBERKR6IpEIkyc/B+AE6psi7AdwC8B+GBiMvN8UYTWX4CIivdCuYjYtLQ0REREtbg8LC+Mas0REOqyhSUJmZiYAAba2f4FNEYhIn7Xrmtnbt2/D2PjBO1//djAjI9y5c+ehgyIiItVrrkmCiYkJNm/eDA8PDzZFICK91K4zs+7u7rhw4UKL28+fP/+3GwuIiEhXNDRJyMiwRf2SXDcB/Io7d3ogKmomTE1N5csuEhHpk3YVs6NHj8b777+PqqqqJtvu3buHxYsXY+zYsSoLjoiIHl7TJgkRqF/FIAKCkAhBGILZs6ObrB1ORKQP2nWZwaJFi7B371706tULs2bNQu/evSESiZCRkYENGzZAKpVi4cKF6oqViIiU0FaTBJksBunpkUhLS+O6skSkd9pVzIrFYpw8eRKvv/46YmJi5L/Fi0QihIeHY+PGjRCLxWoJlIiIlLNq1ar7X7FJAhF1PO1umtCtWzfs27cPRUVFuHLlCgRBgI+PD+zs7NQRHxERKamyshJxcXHYs2fP/RE2SSCijkepDmAAYGdnh0GDBqkyFiIiUqFXX30VO3fuBADY2TmhuDgWghCKxpca1DdJ8PVlkwQi0k/tugGMiIj0w7Vr17B7924AwPz58/H11ztgYHASBgYjwSYJRNSRKH1mloiIlCMIAlJTUyGRSFSytqsgCLh8+TJu376N7Oxs/Prrr/j2228hlUoREhKClStXAgAOHNh/f53ZSPm+vr4BWLuWTRKISH+xmCUi0qDmGhf07dsf69atUbigFAQBVVVVKCsrw4kTJ7BmzRocO3asyTwvLy+sW7dO/jw0NBQXL6YiLS0Nubm5bJJARB0Ci1kiIg1paFwgkw1B/TJZ/gAuICMjFmFh4ThwYD9CQ0Ob7JednY2MjAzcvXsXR48exVdffYXKyspGc0QiEXr06AGxWIzg4GA8/fTTePTRR2FgYNBkXmBgIFctIKIOg8UsEZEGNG1c0HATljsEIRQyWSimTJmGd96JRllZGa5cuYKamhr89ddfSElJafG4NjY2mD59OmbOnAkvLy9NvBUiIp3CYpaISM0EQcDJkydbbVwgCO8iNzcSc+bMabK/kZERfH194ezsDHd3d0ydOhVBQUGwsLCAoaGhBt4BEZHuYjFLRKQmSUlJmDNnDtLT01FTU3N/tPXGBcOGDYOPjw/q6urg5+cHDw8PhISEwMXFRSMxExHpGxazREQqIJPJkJiYiD/++ANZWVm4du0ajhw50szM1hsXrF27Fv7+/sjIyICvry/PvBIRtYHFLBFRG+rq6nDnzh1IJBJIJBJcu3YNP/30EzIzM1FRUYHKysomN2Q1iIiIwIYNG+Ds7IzHHgtGRkbbjQtkMplG3hcRUUfAYpaIOr3q6mqcOHEC6enpKCoqwu3bt3Ht2jV58Xrnzh2FCkxLS0s8/fTT6NmzJ7y9veHr64tHHnlEvqLAunVrEBYWDkEYCZksBg2rGRgYxN5vXHCAy2QREbUTi1ki6pDq6upQUFCA/Px83LlzB7dv38bly5dx/vx5FBcXo7y8HIWFhSgoKEBJSUmbxzMwMIBYLIarqytcXFwwfPhwDBs2DFZWVrCwsIC5uTns7OxgYmLS4jFGjBjBxgVERCrGYpaINO5hOmAJgoA7d+6gvLwcFRUVjR43b97Ejh07cOnSJRQVFbUrJltbWzz55JNwcnKCo6MjunfvDnd3d7i4uMDV1RVOTk4quX6VjQuIiFSLxSwRaVRbHbCqq6tx48YN7NmzB1lZWSgqKsLdu3fl/xYWFrZ4feqDRCIRHB0d4eTkBGdnZ3h5eSEgIABisRgWFhawt7eXb7ezs2vSYEBd2LiAiEh1WMwSkca03AFrBUJDR6FLF1Pcu3dPoWOZm5vDwsKi0cPS0hIjR47E2LFj4eLiAnt7e64GQETUwbGYJSKNaKsDFjAC9+6dAgAYGhpi+PDhCAkJgb29Pezt7WFnZyf/19PTs9VrU4mIqPNgMUtEapeUlIR///vfrXbAAhYCiERSUhKCg4N5RpWIiBTCYpaI1Grv3r149tln/zbSeges8vJyFrJERKQwzdztQESdjiAIOHToEF577TUA+NuyUxda2KN+3M3NTf3BERFRh8Ezs0SkEnV1dbh06RKOHz+O5ORknD17FikpKQCAPn364Oeff0ZQ0OMKdcAiIiJSFItZIlKaIAi4dOkSTp06hYULF0IikTTabmZmhqeeegrLly9Hly5d2AGLiIhUjsUsEbWpuroahYWF2L9/P+Lj45GXl4f8/Hzcvn270Zqv5ubmGDBgAJ588kn4+Phg1KhRjS4bYAcsIiJSNa0Xsxs3bsTHH38MiUSCfv36IS4uDsOGDWtxflJSEqKjo/Hnn3/Czc0N77zzDmbMmKHBiIk6rtu3b2Pfvn24evUqSktLUVJSggsXLuDmzZst7mNubo5evXrh2Wefxbx582Bqatrqa7ADFhERqZJWi9k9e/Zgzpw52LhxI4KDg7F582ZERkYiPT0dnp6eTeZnZWVh9OjRmD59Or766iucOHECUVFRcHJyeuBuaSJqS11dHSQSCXJycpCbm4sdO3bg559/hiAIzc43NDSEo6Mjnn/+eQwdOhTOzs4Qi8Xw8PBos4B9EDtgERGRqmi1mF29ejVefvllvPLKKwCAuLg4HDhwAJs2bUJsbGyT+Z9++ik8PT0RFxcHAPD19cXZs2fx73//u93FbHx8PCwtLeXPm/sBrokxbb62LsQjk8kgkUjg6uoKAwMDvfksZDJZux4P7tPa8/bMbWnf5sZqa2tRVVWFqqoqlJWV4datW5BKpU3eX3BwMEaNGgVbW1tYW1vLW8Da2tpqrN0rERGRorRWzNbU1CA5ORkLFixoNB4WFoaTJ082u8+pU6cQFhbWaCw8PBxffPEFamtrYWz84ELs9df6VVdXy5+XlJQAAF588cWHfAdE+s/Q0BDu7u5wdnaGr68vRo8ejZEjRza7zmvD9w6pn1QqRXl5OYqKirjmrpYxF7qDudAt6s5HaWkpgJZPNv2d1orZgoICSKVSiMXiRuNisRh5eXnN7pOXl9fs/Lq6OhQUFMDV1bXJPrGxsVi6dKnqAifqQKRSKXJycpCTk4OzZ89ix44d2g6JiIhIrqysDDY2Nq3O0foNYA/e9CEIQqs3gjQ3v7nxBjExMYiOjpY/Ly4uRrdu3ZCTk9Pmh0PqV1paiq5du+LGjRuwtrbWdjidGnOhO5gL3cFc6A7mQreoOx+CIKCsrEyhRjpaK2YdHR1haGjY5Cxsfn5+k7OvDVxcXJqdb2RkBAcHh2b3MTU1bfbmFBsbG34z6BBra2vmQ0cwF7qDudAdzIXuYC50izrzoehJR63dzWFiYoKBAwciISGh0XhCQgKGDBnS7D6DBw9uMv/gwYMICgpq9npZIiIiIurYtHprcnR0ND7//HNs3boVGRkZmDt3LnJycuTrxsbExGDKlCny+TNmzMD169cRHR2NjIwMbN26FV988QXefvttbb0FIiIiItIirV4zO3HiRBQWFmLZsmWQSCTw8/PDvn370K1bNwCQr4HZwNvbG/v27cPcuXOxYcMGuLm5Ye3ate1alsvU1BSLFy9u97qYpB7Mh+5gLnQHc6E7mAvdwVzoFl3Kh0hQZM0DIiIiIiIdxBXQiYiIiEhvsZglIiIiIr3FYpaIiIiI9BaLWSIiIiLSW52umN24cSO8vb1hZmaGgQMH4tixY9oOqcOLjY3FoEGDYGVlBWdnZ4wfPx6XLl1qNEcQBCxZsgRubm7o0qULnnzySfz5559airjziI2NhUgkwpw5c+RjzIXm3Lp1C88//zwcHBxgbm6OwMBAJCcny7czF5pRV1eHRYsWwdvbG126dEH37t2xbNkyyGQy+RzmQn2OHj2Kf/zjH3Bzc4NIJMKPP/7YaLsin311dTXeeOMNODo6wsLCAuPGjcPNmzc1+C46htZyUVtbi/nz58Pf3x8WFhZwc3PDlClTkJub2+gY2shFpypm9+zZgzlz5mDhwoVISUnBsGHDEBkZ2Wj5L1K9pKQkzJw5E7///jsSEhJQV1eHsLAwVFRUyOd89NFHWL16NdavX48zZ87AxcUFo0aNQllZmRYj79jOnDmDLVu2oH///o3GmQvNKCoqQnBwMIyNjfHrr78iPT0dn3zyCWxtbeVzmAvNWLVqFT799FOsX78eGRkZ+Oijj/Dxxx9j3bp18jnMhfpUVFQgICAA69evb3a7Ip/9nDlz8MMPP2D37t04fvw4ysvLMXbsWEilUk29jQ6htVxUVlbi3LlzeO+993Du3Dns3bsXmZmZGDduXKN5WsmF0Ik8+uijwowZMxqN9enTR1iwYIGWIuqc8vPzBQBCUlKSIAiCIJPJBBcXF2HlypXyOVVVVYKNjY3w6aefaivMDq2srEzw8fEREhIShOHDhwtvvvmmIAjMhSbNnz9fGDp0aIvbmQvNGTNmjDBt2rRGY88884zw/PPPC4LAXGgSAOGHH36QP1fksy8uLhaMjY2F3bt3y+fcunVLMDAwEPbv36+x2DuaB3PRnNOnTwsAhOvXrwuCoL1cdJozszU1NUhOTkZYWFij8bCwMJw8eVJLUXVOJSUlAAB7e3sAQFZWFvLy8hrlxtTUFMOHD2du1GTmzJkYM2YMRo4c2WicudCc+Ph4BAUF4Z///CecnZ0xYMAAfPbZZ/LtzIXmDB06FImJicjMzAQApKWl4fjx4xg9ejQA5kKbFPnsk5OTUVtb22iOm5sb/Pz8mB81KykpgUgkkv9FSVu50GoHME0qKCiAVCqFWCxuNC4Wi5GXl6elqDofQRAQHR2NoUOHws/PDwDkn39zubl+/brGY+zodu/ejeTkZJw9e7bJNuZCc65du4ZNmzYhOjoa7777Lk6fPo3Zs2fD1NQUU6ZMYS40aP78+SgpKUGfPn1gaGgIqVSK5cuXY9KkSQD4faFNinz2eXl5MDExgZ2dXZM5/PmuPlVVVViwYAEmT54Ma2trANrLRacpZhuIRKJGzwVBaDJG6jNr1iycP38ex48fb7KNuVG/Gzdu4M0338TBgwdhZmbW4jzmQv1kMhmCgoKwYsUKAMCAAQPw559/YtOmTZgyZYp8HnOhfnv27MFXX32FXbt2oV+/fkhNTcWcOXPg5uaGqVOnyucxF9qjzGfP/KhPbW0t/vWvf0Emk2Hjxo1tzld3LjrNZQaOjo4wNDRs8ptBfn5+k9/4SD3eeOMNxMfH4/Dhw/Dw8JCPu7i4AABzowHJycnIz8/HwIEDYWRkBCMjIyQlJWHt2rUwMjKSf97Mhfq5urqib9++jcZ8fX3lN6Ty+0Jz5s2bhwULFuBf//oX/P398cILL2Du3LmIjY0FwFxokyKfvYuLC2pqalBUVNTiHFKd2tpaPPfcc8jKykJCQoL8rCygvVx0mmLWxMQEAwcOREJCQqPxhIQEDBkyREtRdQ6CIGDWrFnYu3cvDh06BG9v70bbvb294eLi0ig3NTU1SEpKYm5ULDQ0FBcuXEBqaqr8ERQUhP/7v/9DamoqunfvzlxoSHBwcJMl6jIzM9GtWzcA/L7QpMrKShgYNP5xaGhoKF+ai7nQHkU++4EDB8LY2LjRHIlEgosXLzI/KtZQyF6+fBm//fYbHBwcGm3XWi7UdmuZDtq9e7dgbGwsfPHFF0J6erowZ84cwcLCQsjOztZ2aB3a66+/LtjY2AhHjhwRJBKJ/FFZWSmfs3LlSsHGxkbYu3evcOHCBWHSpEmCq6urUFpaqsXIO4e/r2YgCMyFppw+fVowMjISli9fLly+fFnYuXOnYG5uLnz11VfyOcyFZkydOlVwd3cXfv75ZyErK0vYu3ev4OjoKLzzzjvyOcyF+pSVlQkpKSlCSkqKAEBYvXq1kJKSIr9DXpHPfsaMGYKHh4fw22+/CefOnRNGjBghBAQECHV1ddp6W3qptVzU1tYK48aNEzw8PITU1NRGP8+rq6vlx9BGLjpVMSsIgrBhwwahW7dugomJifDII4/Il4ci9QHQ7GPbtm3yOTKZTFi8eLHg4uIimJqaCk888YRw4cIF7QXdiTxYzDIXmvPTTz8Jfn5+gqmpqdCnTx9hy5YtjbYzF5pRWloqvPnmm4Knp6dgZmYmdO/eXVi4cGGjH9DMhfocPny42Z8RU6dOFQRBsc/+3r17wqxZswR7e3uhS5cuwtixY4WcnBwtvBv91lousrKyWvx5fvjwYfkxtJELkSAIgvrO+xIRERERqU+nuWaWiIiIiDoeFrNEREREpLdYzBIRERGR3mIxS0RERER6i8UsEREREektFrNEREREpLdYzBIRERGR3mIxS0RERER6i8UsEZGWLFmyBIGBgVp7/ffeew+vvvqqQnPffvttzJ49W80RERG1HzuAERGpgUgkanX71KlTsX79elRXV8PBwUFDUf3P7du34ePjg/Pnz8PLy6vN+fn5+ejRowfOnz8Pb29v9QdIRKQgFrNERGqQl5cn/3rPnj14//33cenSJflYly5dYGNjo43QAAArVqxAUlISDhw4oPA+zz77LHr27IlVq1apMTIiovbhZQZERGrg4uIif9jY2EAkEjUZe/AygxdffBHjx4/HihUrIBaLYWtri6VLl6Kurg7z5s2Dvb09PDw8sHXr1kavdevWLUycOBF2dnZwcHDAU089hezs7Fbj2717N8aNG9do7LvvvoO/vz+6dOkCBwcHjBw5EhUVFfLt48aNw9dff/3Qnw0RkSqxmCUi0iGHDh1Cbm4ujh49itWrV2PJkiUYO3Ys7Ozs8Mcff2DGjBmYMWMGbty4AQCorKxESEgILC0tcfToURw/fhyWlpaIiIhATU1Ns69RVFSEixcvIigoSD4mkUgwadIkTJs2DRkZGThy5AieeeYZ/P2Pd48++ihu3LiB69evq/dDICJqBxazREQ6xN7eHmvXrkXv3r0xbdo09O7dG5WVlXj33Xfh4+ODmJgYmJiY4MSJEwDqz7AaGBjg888/h7+/P3x9fbFt2zbk5OTgyJEjzb7G9evXIQgC3Nzc5GMSiQR1dXV45pln4OXlBX9/f0RFRcHS0lI+x93dHQDaPOtLRKRJRtoOgIiI/qdfv34wMPjfeQaxWAw/Pz/5c0NDQzg4OCA/Px8AkJycjCtXrsDKyqrRcaqqqnD16tVmX+PevXsAADMzM/lYQEAAQkND4e/vj/DwcISFhWHChAmws7OTz+nSpQuA+rPBRES6gsUsEZEOMTY2bvRcJBI1OyaTyQAAMpkMAwcOxM6dO5scy8nJqdnXcHR0BFB/uUHDHENDQyQkJODkyZM4ePAg1q1bh4ULF+KPP/6Qr15w9+7dVo9LRKQNvMyAiEiPPfLII7h8+TKcnZ3Rs2fPRo+WVkvo0aMHrK2tkZ6e3mhcJBIhODgYS5cuRUpKCkxMTPDDDz/It1+8eBHGxsbo16+fWt8TEVF7sJglItJj//d//wdHR0c89dRTOHbsGLKyspCUlIQ333wTN2/ebHYfAwMDjBw5EsePH5eP/fHHH1ixYgXOnj2LnJwc7N27F3fu3IGvr698zrFjxzBs2DD55QZERLqAxSwRkR4zNzfH0aNH4enpiWeeeQa+vr6YNm0a7t27B2tr6xb3e/XVV7F792755QrW1tY4evQoRo8ejV69emHRokX45JNPEBkZKd/n66+/xvTp09X+noiI2oNNE4iIOiFBEPD4449jzpw5mDRpUpvzf/nlF8ybNw/nz5+HkRFvtyAi3cEzs0REnZBIJMKWLVtQV1en0PyKigps27aNhSwR6RyemSUiIiIivcUzs0RERESkt1jMEhEREZHeYjFLRERERHqLxSwRERER6S0Ws0RERESkt1jMEhEREZHeYjFLRERERHqLxSwRERER6S0Ws0RERESkt/4/NXHOS0AzJwcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 800x250 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,2.5))\n",
    "ax.plot(gmtime,Ia_hist,'k')\n",
    "ax.scatter(Ia_times,Ia*percentages,zorder=2,s=35,facecolor='b',edgecolor='k')\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Cumulative Ia (m/s)')\n",
    "ax.set_xlim([0,gmtime[-1]])\n",
    "ax.set_ylim([0,None])\n",
    "plt.grid(which='both',color='lightgray')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "338aaf0e",
   "metadata": {},
   "source": [
    "## CAV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ea9242ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAV = 13.717 m/s\n"
     ]
    }
   ],
   "source": [
    "CAV = gmimtools.get_CAV_cy(accel,dt)\n",
    "print('CAV = {:.3f} m/s'.format(CAV))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1b09c2f",
   "metadata": {},
   "source": [
    "## CAV5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "12c7c279",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAV5 = 12.915 m/s\n"
     ]
    }
   ],
   "source": [
    "CAV5 = gmimtools.get_CAVn_cy(accel,dt,0.05)\n",
    "print('CAV5 = {:.3f} m/s'.format(CAV5))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5af43cd7",
   "metadata": {},
   "source": [
    "## CAVstd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c68b10bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAVstd = 11.544 m/s\n"
     ]
    }
   ],
   "source": [
    "CAVstd = gmimtools.get_CAVstd_cy(accel,dt)\n",
    "print('CAVstd = {:.3f} m/s'.format(CAVstd))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a50ab5e8",
   "metadata": {},
   "source": [
    "## CAVdp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69591a54",
   "metadata": {},
   "source": [
    "#### Read three-component record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "061ff619",
   "metadata": {},
   "outputs": [],
   "source": [
    "[accel_H1,dt] = gmimtools.read_AT2(path+files[0],scaling)\n",
    "[accel_H2,dt] = gmimtools.read_AT2(path+files[1],scaling)\n",
    "[accel_V,dt]  = gmimtools.read_AT2(path+files[2],scaling)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "600dadd4",
   "metadata": {},
   "source": [
    "#### Compute PSA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "55ed76ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "osc_damping = 0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "972f60ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "freqs1 = np.logspace(0,np.log10(2),30)\n",
    "freqs2 = np.logspace(np.log10(2),1,30)\n",
    "freqs  = np.unique(np.concatenate([freqs1,freqs2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "ebc56c49",
   "metadata": {},
   "outputs": [],
   "source": [
    "PSA_H1 = pyrotd.calc_spec_accels(dt,accel_H1,freqs,osc_damping).spec_accel\n",
    "PSA_H2 = pyrotd.calc_spec_accels(dt,accel_H2,freqs,osc_damping).spec_accel\n",
    "PSA_V  = pyrotd.calc_spec_accels(dt,accel_V,freqs,osc_damping).spec_accel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f2c1291",
   "metadata": {},
   "source": [
    "#### Compute PSV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "eef00c95",
   "metadata": {},
   "outputs": [],
   "source": [
    "w_freq = 2*np.pi*freqs\n",
    "PSV_H1 = PSA_H1/w_freq\n",
    "PSV_H2 = PSA_H2/w_freq\n",
    "PSV_V  = PSA_V/w_freq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29644b72",
   "metadata": {},
   "source": [
    "#### Compute CAVstd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "80a74456",
   "metadata": {},
   "outputs": [],
   "source": [
    "CAVstd_H1 = gmimtools.get_CAVstd_cy(accel_H1,dt)\n",
    "CAVstd_H2 = gmimtools.get_CAVstd_cy(accel_H2,dt)\n",
    "CAVstd_V  = gmimtools.get_CAVstd_cy(accel_V,dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad96aa28",
   "metadata": {},
   "source": [
    "#### Compute CAVdp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "3cdbad64",
   "metadata": {},
   "outputs": [],
   "source": [
    "PSA_max = max(np.max(PSA_H1[freqs>=2]),np.max(PSA_H2[freqs>=2]),np.max(PSA_V[freqs>=2]))\n",
    "PSV_max = max(np.max(PSA_H1[freqs<=2]),np.max(PSA_H2[freqs<=2]),np.max(PSA_V[freqs<=2]))\n",
    "CAVstd_max = max(CAVstd_H1,CAVstd_H2,CAVstd_V)\n",
    "CAVstd_gm  = (CAVstd_H1*CAVstd_H2)**0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b2e5d3d1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAVdp = 10.629 m/s\n"
     ]
    }
   ],
   "source": [
    "CAVdp = gmimtools.get_CAVdp_cy(CAVstd_gm,CAVstd_max,PSA_max,PSV_max)\n",
    "print('CAVdp = {:.3f} m/s'.format(CAVdp))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5592578",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "6138c156",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>GMIMs</th>\n",
       "      <th>Unit</th>\n",
       "      <th>Value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>PGA</td>\n",
       "      <td>m/s2</td>\n",
       "      <td>1.286023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>PGV</td>\n",
       "      <td>m/s</td>\n",
       "      <td>0.335510</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Ia</td>\n",
       "      <td>m/s</td>\n",
       "      <td>0.717027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CAV</td>\n",
       "      <td>m/s</td>\n",
       "      <td>13.717308</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>CAV5</td>\n",
       "      <td>m/s</td>\n",
       "      <td>12.914516</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>CAVstd</td>\n",
       "      <td>m/s</td>\n",
       "      <td>11.544224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>CAVdp</td>\n",
       "      <td>m/s</td>\n",
       "      <td>10.629479</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    GMIMs  Unit      Value\n",
       "0     PGA  m/s2   1.286023\n",
       "1     PGV   m/s   0.335510\n",
       "2      Ia   m/s   0.717027\n",
       "3     CAV   m/s  13.717308\n",
       "4    CAV5   m/s  12.914516\n",
       "5  CAVstd   m/s  11.544224\n",
       "6   CAVdp   m/s  10.629479"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table = pd.DataFrame({\n",
    "    'GMIMs': ['PGA','PGV','Ia','CAV','CAV5','CAVstd','CAVdp'],\n",
    "    'Unit':  ['m/s2','m/s','m/s','m/s','m/s','m/s','m/s'],\n",
    "    'Value': [PGA,PGV,Ia,CAV,CAV5,CAVstd,CAVdp]\n",
    "})\n",
    "table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90da0304",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}