get_gps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bG3GpKz_mO0t"
   },
   "source": [
    "# Download of GPS time series\n",
    "This notebook collects GPS data given a lat-long box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "T618TOwh8Y9K"
   },
   "outputs": [],
   "source": [
    "# import modules\n",
    "import requests, zipfile, io, gzip, glob, os\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collect information about the GPS site location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 309,
     "resources": {
      "http://localhost:8080/nbextensions/google.colab/files.js": {
       "data": "Ly8gQ29weXJpZ2h0IDIwMTcgR29vZ2xlIExMQwovLwovLyBMaWNlbnNlZCB1bmRlciB0aGUgQXBhY2hlIExpY2Vuc2UsIFZlcnNpb24gMi4wICh0aGUgIkxpY2Vuc2UiKTsKLy8geW91IG1heSBub3QgdXNlIHRoaXMgZmlsZSBleGNlcHQgaW4gY29tcGxpYW5jZSB3aXRoIHRoZSBMaWNlbnNlLgovLyBZb3UgbWF5IG9idGFpbiBhIGNvcHkgb2YgdGhlIExpY2Vuc2UgYXQKLy8KLy8gICAgICBodHRwOi8vd3d3LmFwYWNoZS5vcmcvbGljZW5zZXMvTElDRU5TRS0yLjAKLy8KLy8gVW5sZXNzIHJlcXVpcmVkIGJ5IGFwcGxpY2FibGUgbGF3IG9yIGFncmVlZCB0byBpbiB3cml0aW5nLCBzb2Z0d2FyZQovLyBkaXN0cmlidXRlZCB1bmRlciB0aGUgTGljZW5zZSBpcyBkaXN0cmlidXRlZCBvbiBhbiAiQVMgSVMiIEJBU0lTLAovLyBXSVRIT1VUIFdBUlJBTlRJRVMgT1IgQ09ORElUSU9OUyBPRiBBTlkgS0lORCwgZWl0aGVyIGV4cHJlc3Mgb3IgaW1wbGllZC4KLy8gU2VlIHRoZSBMaWNlbnNlIGZvciB0aGUgc3BlY2lmaWMgbGFuZ3VhZ2UgZ292ZXJuaW5nIHBlcm1pc3Npb25zIGFuZAovLyBsaW1pdGF0aW9ucyB1bmRlciB0aGUgTGljZW5zZS4KCi8qKgogKiBAZmlsZW92ZXJ2aWV3IEhlbHBlcnMgZm9yIGdvb2dsZS5jb2xhYiBQeXRob24gbW9kdWxlLgogKi8KKGZ1bmN0aW9uKHNjb3BlKSB7CmZ1bmN0aW9uIHNwYW4odGV4dCwgc3R5bGVBdHRyaWJ1dGVzID0ge30pIHsKICBjb25zdCBlbGVtZW50ID0gZG9jdW1lbnQuY3JlYXRlRWxlbWVudCgnc3BhbicpOwogIGVsZW1lbnQudGV4dENvbnRlbnQgPSB0ZXh0OwogIGZvciAoY29uc3Qga2V5IG9mIE9iamVjdC5rZXlzKHN0eWxlQXR0cmlidXRlcykpIHsKICAgIGVsZW1lbnQuc3R5bGVba2V5XSA9IHN0eWxlQXR0cmlidXRlc1trZXldOwogIH0KICByZXR1cm4gZWxlbWVudDsKfQoKLy8gTWF4IG51bWJlciBvZiBieXRlcyB3aGljaCB3aWxsIGJlIHVwbG9hZGVkIGF0IGEgdGltZS4KY29uc3QgTUFYX1BBWUxPQURfU0laRSA9IDEwMCAqIDEwMjQ7Ci8vIE1heCBhbW91bnQgb2YgdGltZSB0byBibG9jayB3YWl0aW5nIGZvciB0aGUgdXNlci4KY29uc3QgRklMRV9DSEFOR0VfVElNRU9VVF9NUyA9IDMwICogMTAwMDsKCmZ1bmN0aW9uIF91cGxvYWRGaWxlcyhpbnB1dElkLCBvdXRwdXRJZCkgewogIGNvbnN0IHN0ZXBzID0gdXBsb2FkRmlsZXNTdGVwKGlucHV0SWQsIG91dHB1dElkKTsKICBjb25zdCBvdXRwdXRFbGVtZW50ID0gZG9jdW1lbnQuZ2V0RWxlbWVudEJ5SWQob3V0cHV0SWQpOwogIC8vIENhY2hlIHN0ZXBzIG9uIHRoZSBvdXRwdXRFbGVtZW50IHRvIG1ha2UgaXQgYXZhaWxhYmxlIGZvciB0aGUgbmV4dCBjYWxsCiAgLy8gdG8gdXBsb2FkRmlsZXNDb250aW51ZSBmcm9tIFB5dGhvbi4KICBvdXRwdXRFbGVtZW50LnN0ZXBzID0gc3RlcHM7CgogIHJldHVybiBfdXBsb2FkRmlsZXNDb250aW51ZShvdXRwdXRJZCk7Cn0KCi8vIFRoaXMgaXMgcm91Z2hseSBhbiBhc3luYyBnZW5lcmF0b3IgKG5vdCBzdXBwb3J0ZWQgaW4gdGhlIGJyb3dzZXIgeWV0KSwKLy8gd2hlcmUgdGhlcmUgYXJlIG11bHRpcGxlIGFzeW5jaHJvbm91cyBzdGVwcyBhbmQgdGhlIFB5dGhvbiBzaWRlIGlzIGdvaW5nCi8vIHRvIHBvbGwgZm9yIGNvbXBsZXRpb24gb2YgZWFjaCBzdGVwLgovLyBUaGlzIHVzZXMgYSBQcm9taXNlIHRvIGJsb2NrIHRoZSBweXRob24gc2lkZSBvbiBjb21wbGV0aW9uIG9mIGVhY2ggc3RlcCwKLy8gdGhlbiBwYXNzZXMgdGhlIHJlc3VsdCBvZiB0aGUgcHJldmlvdXMgc3RlcCBhcyB0aGUgaW5wdXQgdG8gdGhlIG5leHQgc3RlcC4KZnVuY3Rpb24gX3VwbG9hZEZpbGVzQ29udGludWUob3V0cHV0SWQpIHsKICBjb25zdCBvdXRwdXRFbGVtZW50ID0gZG9jdW1lbnQuZ2V0RWxlbWVudEJ5SWQob3V0cHV0SWQpOwogIGNvbnN0IHN0ZXBzID0gb3V0cHV0RWxlbWVudC5zdGVwczsKCiAgY29uc3QgbmV4dCA9IHN0ZXBzLm5leHQob3V0cHV0RWxlbWVudC5sYXN0UHJvbWlzZVZhbHVlKTsKICByZXR1cm4gUHJvbWlzZS5yZXNvbHZlKG5leHQudmFsdWUucHJvbWlzZSkudGhlbigodmFsdWUpID0+IHsKICAgIC8vIENhY2hlIHRoZSBsYXN0IHByb21pc2UgdmFsdWUgdG8gbWFrZSBpdCBhdmFpbGFibGUgdG8gdGhlIG5leHQKICAgIC8vIHN0ZXAgb2YgdGhlIGdlbmVyYXRvci4KICAgIG91dHB1dEVsZW1lbnQubGFzdFByb21pc2VWYWx1ZSA9IHZhbHVlOwogICAgcmV0dXJuIG5leHQudmFsdWUucmVzcG9uc2U7CiAgfSk7Cn0KCi8qKgogKiBHZW5lcmF0b3IgZnVuY3Rpb24gd2hpY2ggaXMgY2FsbGVkIGJldHdlZW4gZWFjaCBhc3luYyBzdGVwIG9mIHRoZSB1cGxvYWQKICogcHJvY2Vzcy4KICogQHBhcmFtIHtzdHJpbmd9IGlucHV0SWQgRWxlbWVudCBJRCBvZiB0aGUgaW5wdXQgZmlsZSBwaWNrZXIgZWxlbWVudC4KICogQHBhcmFtIHtzdHJpbmd9IG91dHB1dElkIEVsZW1lbnQgSUQgb2YgdGhlIG91dHB1dCBkaXNwbGF5LgogKiBAcmV0dXJuIHshSXRlcmFibGU8IU9iamVjdD59IEl0ZXJhYmxlIG9mIG5leHQgc3RlcHMuCiAqLwpmdW5jdGlvbiogdXBsb2FkRmlsZXNTdGVwKGlucHV0SWQsIG91dHB1dElkKSB7CiAgY29uc3QgaW5wdXRFbGVtZW50ID0gZG9jdW1lbnQuZ2V0RWxlbWVudEJ5SWQoaW5wdXRJZCk7CiAgaW5wdXRFbGVtZW50LmRpc2FibGVkID0gZmFsc2U7CgogIGNvbnN0IG91dHB1dEVsZW1lbnQgPSBkb2N1bWVudC5nZXRFbGVtZW50QnlJZChvdXRwdXRJZCk7CiAgb3V0cHV0RWxlbWVudC5pbm5lckhUTUwgPSAnJzsKCiAgY29uc3QgcGlja2VkUHJvbWlzZSA9IG5ldyBQcm9taXNlKChyZXNvbHZlKSA9PiB7CiAgICBpbnB1dEVsZW1lbnQuYWRkRXZlbnRMaXN0ZW5lcignY2hhbmdlJywgKGUpID0+IHsKICAgICAgcmVzb2x2ZShlLnRhcmdldC5maWxlcyk7CiAgICB9KTsKICB9KTsKCiAgY29uc3QgY2FuY2VsID0gZG9jdW1lbnQuY3JlYXRlRWxlbWVudCgnYnV0dG9uJyk7CiAgaW5wdXRFbGVtZW50LnBhcmVudEVsZW1lbnQuYXBwZW5kQ2hpbGQoY2FuY2VsKTsKICBjYW5jZWwudGV4dENvbnRlbnQgPSAnQ2FuY2VsIHVwbG9hZCc7CiAgY29uc3QgY2FuY2VsUHJvbWlzZSA9IG5ldyBQcm9taXNlKChyZXNvbHZlKSA9PiB7CiAgICBjYW5jZWwub25jbGljayA9ICgpID0+IHsKICAgICAgcmVzb2x2ZShudWxsKTsKICAgIH07CiAgfSk7CgogIC8vIENhbmNlbCB1cGxvYWQgaWYgdXNlciBoYXNuJ3QgcGlja2VkIGFueXRoaW5nIGluIHRpbWVvdXQuCiAgY29uc3QgdGltZW91dFByb21pc2UgPSBuZXcgUHJvbWlzZSgocmVzb2x2ZSkgPT4gewogICAgc2V0VGltZW91dCgoKSA9PiB7CiAgICAgIHJlc29sdmUobnVsbCk7CiAgICB9LCBGSUxFX0NIQU5HRV9USU1FT1VUX01TKTsKICB9KTsKCiAgLy8gV2FpdCBmb3IgdGhlIHVzZXIgdG8gcGljayB0aGUgZmlsZXMuCiAgY29uc3QgZmlsZXMgPSB5aWVsZCB7CiAgICBwcm9taXNlOiBQcm9taXNlLnJhY2UoW3BpY2tlZFByb21pc2UsIHRpbWVvdXRQcm9taXNlLCBjYW5jZWxQcm9taXNlXSksCiAgICByZXNwb25zZTogewogICAgICBhY3Rpb246ICdzdGFydGluZycsCiAgICB9CiAgfTsKCiAgaWYgKCFmaWxlcykgewogICAgcmV0dXJuIHsKICAgICAgcmVzcG9uc2U6IHsKICAgICAgICBhY3Rpb246ICdjb21wbGV0ZScsCiAgICAgIH0KICAgIH07CiAgfQoKICBjYW5jZWwucmVtb3ZlKCk7CgogIC8vIERpc2FibGUgdGhlIGlucHV0IGVsZW1lbnQgc2luY2UgZnVydGhlciBwaWNrcyBhcmUgbm90IGFsbG93ZWQuCiAgaW5wdXRFbGVtZW50LmRpc2FibGVkID0gdHJ1ZTsKCiAgZm9yIChjb25zdCBmaWxlIG9mIGZpbGVzKSB7CiAgICBjb25zdCBsaSA9IGRvY3VtZW50LmNyZWF0ZUVsZW1lbnQoJ2xpJyk7CiAgICBsaS5hcHBlbmQoc3BhbihmaWxlLm5hbWUsIHtmb250V2VpZ2h0OiAnYm9sZCd9KSk7CiAgICBsaS5hcHBlbmQoc3BhbigKICAgICAgICBgKCR7ZmlsZS50eXBlIHx8ICduL2EnfSkgLSAke2ZpbGUuc2l6ZX0gYnl0ZXMsIGAgKwogICAgICAgIGBsYXN0IG1vZGlmaWVkOiAkewogICAgICAgICAgICBmaWxlLmxhc3RNb2RpZmllZERhdGUgPyBmaWxlLmxhc3RNb2RpZmllZERhdGUudG9Mb2NhbGVEYXRlU3RyaW5nKCkgOgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAnbi9hJ30gLSBgKSk7CiAgICBjb25zdCBwZXJjZW50ID0gc3BhbignMCUgZG9uZScpOwogICAgbGkuYXBwZW5kQ2hpbGQocGVyY2VudCk7CgogICAgb3V0cHV0RWxlbWVudC5hcHBlbmRDaGlsZChsaSk7CgogICAgY29uc3QgZmlsZURhdGFQcm9taXNlID0gbmV3IFByb21pc2UoKHJlc29sdmUpID0+IHsKICAgICAgY29uc3QgcmVhZGVyID0gbmV3IEZpbGVSZWFkZXIoKTsKICAgICAgcmVhZGVyLm9ubG9hZCA9IChlKSA9PiB7CiAgICAgICAgcmVzb2x2ZShlLnRhcmdldC5yZXN1bHQpOwogICAgICB9OwogICAgICByZWFkZXIucmVhZEFzQXJyYXlCdWZmZXIoZmlsZSk7CiAgICB9KTsKICAgIC8vIFdhaXQgZm9yIHRoZSBkYXRhIHRvIGJlIHJlYWR5LgogICAgbGV0IGZpbGVEYXRhID0geWllbGQgewogICAgICBwcm9taXNlOiBmaWxlRGF0YVByb21pc2UsCiAgICAgIHJlc3BvbnNlOiB7CiAgICAgICAgYWN0aW9uOiAnY29udGludWUnLAogICAgICB9CiAgICB9OwoKICAgIC8vIFVzZSBhIGNodW5rZWQgc2VuZGluZyB0byBhdm9pZCBtZXNzYWdlIHNpemUgbGltaXRzLiBTZWUgYi82MjExNTY2MC4KICAgIGxldCBwb3NpdGlvbiA9IDA7CiAgICB3aGlsZSAocG9zaXRpb24gPCBmaWxlRGF0YS5ieXRlTGVuZ3RoKSB7CiAgICAgIGNvbnN0IGxlbmd0aCA9IE1hdGgubWluKGZpbGVEYXRhLmJ5dGVMZW5ndGggLSBwb3NpdGlvbiwgTUFYX1BBWUxPQURfU0laRSk7CiAgICAgIGNvbnN0IGNodW5rID0gbmV3IFVpbnQ4QXJyYXkoZmlsZURhdGEsIHBvc2l0aW9uLCBsZW5ndGgpOwogICAgICBwb3NpdGlvbiArPSBsZW5ndGg7CgogICAgICBjb25zdCBiYXNlNjQgPSBidG9hKFN0cmluZy5mcm9tQ2hhckNvZGUuYXBwbHkobnVsbCwgY2h1bmspKTsKICAgICAgeWllbGQgewogICAgICAgIHJlc3BvbnNlOiB7CiAgICAgICAgICBhY3Rpb246ICdhcHBlbmQnLAogICAgICAgICAgZmlsZTogZmlsZS5uYW1lLAogICAgICAgICAgZGF0YTogYmFzZTY0LAogICAgICAgIH0sCiAgICAgIH07CiAgICAgIHBlcmNlbnQudGV4dENvbnRlbnQgPQogICAgICAgICAgYCR7TWF0aC5yb3VuZCgocG9zaXRpb24gLyBmaWxlRGF0YS5ieXRlTGVuZ3RoKSAqIDEwMCl9JSBkb25lYDsKICAgIH0KICB9CgogIC8vIEFsbCBkb25lLgogIHlpZWxkIHsKICAgIHJlc3BvbnNlOiB7CiAgICAgIGFjdGlvbjogJ2NvbXBsZXRlJywKICAgIH0KICB9Owp9CgpzY29wZS5nb29nbGUgPSBzY29wZS5nb29nbGUgfHwge307CnNjb3BlLmdvb2dsZS5jb2xhYiA9IHNjb3BlLmdvb2dsZS5jb2xhYiB8fCB7fTsKc2NvcGUuZ29vZ2xlLmNvbGFiLl9maWxlcyA9IHsKICBfdXBsb2FkRmlsZXMsCiAgX3VwbG9hZEZpbGVzQ29udGludWUsCn07Cn0pKHNlbGYpOwo=",
       "headers": [
        [
         "content-type",
         "application/javascript"
        ]
       ],
       "ok": true,
       "status": 200,
       "status_text": ""
      }
     }
    },
    "id": "5FyBzHzEm7de",
    "outputId": "4369d988-9eae-4e68-cdcf-7d3460df32f1"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     station   latitude   longitude    elevation\n",
      "0       AB01  52.209503 -174.204758    25.483989\n",
      "1       AB02  52.970604 -168.854671   192.803423\n",
      "2       AB04  63.656863 -170.567443   136.547210\n",
      "3       AB06  54.885321 -163.423455   500.432692\n",
      "4       AB07  55.349275 -160.476760    89.571502\n",
      "...      ...        ...         ...          ...\n",
      "1805    WWAY -81.577000  -28.402700  1195.388000\n",
      "1806    YBRA  45.122100 -109.267400  2072.518000\n",
      "1807    YESX  28.378336 -108.919639  1506.487003\n",
      "1808    YMER  77.432895  -24.326329  1070.068249\n",
      "1809    YUMX  32.027692 -115.199203   -30.200000\n",
      "\n",
      "[1810 rows x 4 columns]\n"
     ]
    }
   ],
   "source": [
    "# read csv file\n",
    "gps_loc = pd.read_csv(\"gps_locations.csv\",usecols=['station','latitude','longitude','elevation'],\n",
    "    dtype={'latitude': np.float64, 'longitude': np.float64,'elevation':np.float64})\n",
    "print(gps_loc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick a station\n",
    "stalist = \"AB01\"\n",
    "# index of the station chosen\n",
    "ista=gps_loc[gps_loc[\"station\"]==stalist].index.values.astype('int64')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 71
    },
    "id": "JAvnzb8_DnOm",
    "outputId": "71cbe680-d972-4552-b310-687dcfde57a4"
   },
   "outputs": [],
   "source": [
    "# OR select within a lat long box\n",
    "latmin = 35.\n",
    "latmax = 36.\n",
    "longmin = -120.\n",
    "longmax = -119.\n",
    "ik=gps_loc[ (gps_loc[\"latitude\"]>=latmin )& (gps_loc[\"latitude\"]<latmax ) & (gps_loc[\"longitude\"]>=longmin )& (gps_loc[\"longitude\"]<longmax )].index.values.astype('int64')\n",
    "# "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['BVPP', 'P522', 'P537', 'P543', 'P544', 'P545', 'P547', 'P563', 'P564', 'P565', 'P809', 'P810']\n",
      "12\n",
      "BVPP\n"
     ]
    }
   ],
   "source": [
    "stalist = gps_loc[\"station\"].iloc[ik].values.tolist()\n",
    "print(stalist)\n",
    "print(len(stalist))\n",
    "print(stalist[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-mT97In6fp9p"
   },
   "source": [
    "## download data\n",
    "Build URLs and download from the FTP server."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "NLao9dfvmy-T",
    "outputId": "a7c9e478-77f9-4341-9a2e-cd103efce4bf"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Downloading and processing BVPP\n",
      "writing to diskgps_BVPP.csv\n",
      "Downloading and processing P522\n",
      "writing to diskgps_P522.csv\n",
      "Downloading and processing P537\n"
     ]
    }
   ],
   "source": [
    "# initialize nested lists\n",
    "df=[]\n",
    "for ista,sta in enumerate(stalist):\n",
    "    dfsta=[]\n",
    "    print(\"Downloading and processing \" + sta )\n",
    "    \n",
    "    \n",
    "    # what's the lat long el\n",
    "    iista=gps_loc[gps_loc[\"station\"]==sta].index.values.astype('int64')[0]\n",
    "    lat=gps_loc[\"latitude\"].iloc[iista]\n",
    "    lon=gps_loc[\"longitude\"].iloc[iista]\n",
    "    el=gps_loc[\"elevation\"].iloc[iista]\n",
    "    \n",
    "    \n",
    "    file_url=\"http://geodesy.unr.edu/gps_timeseries/tenv/IGS14/\"+ sta + \".tenv\"\n",
    "    r = requests.get(file_url).text.splitlines()  # download, read text, split lines into a list\n",
    "    ue=[];un=[];uv=[];se=[];sn=[];sv=[];date=[];date_year=[]\n",
    "    for iday in r:  # this loops through the days of data\n",
    "        crap=iday.split()\n",
    "        if len(crap)<10:continue\n",
    "        date.append((crap[1]))\n",
    "        date_year.append(float(crap[2]))\n",
    "        ue.append(float(crap[6])*1000)\n",
    "        un.append(float(crap[7])*1000)\n",
    "        uv.append(float(crap[8])*1000)\n",
    "#     #             # errors\n",
    "#         se.append(float(crap[10])*1000)\n",
    "#         sn.append(float(crap[11])*1000)\n",
    "#         sv.append(float(crap[12])*1000)\n",
    "\n",
    "          # make dataframe\n",
    "        crap={'station':sta,'latitude':lat,'longitude':lon,'elevation':el,'date':date,'date_year':date_year,'east':ue,'north':un,'up':uv}\n",
    "        if len(df)==0:\n",
    "            df = pd.DataFrame(crap, columns = ['station', 'latitude', 'longitude','elevation', 'date','date_year','east','north','up'])\n",
    "        else:\n",
    "            df=pd.concat([df,pd.DataFrame(crap, columns = ['station', 'date','date_year','east','north','up'])])\n",
    "        if len(dfsta)==0:\n",
    "            dfsta = pd.DataFrame(crap, columns = ['station', 'latitude', 'longitude','elevation', 'date','date_year','east','north','up'])\n",
    "        else:\n",
    "            dfstadf=pd.concat([dfsta,pd.DataFrame(crap, columns = ['station', 'latitude', 'longitude','elevation','date','date_year','east','north','up'])])\n",
    "\n",
    "    print(\"writing to disk\"  + \"gps_\"+sta+\".csv\")   \n",
    "    dfsta.to_csv(\"gps_\"+sta+\".csv\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot time series\n",
    "fig,ax=plt.subplots(3,1,figsize=(11,8),sharex=True)\n",
    "for ista,sta in enumerate(stalist):\n",
    "    ax[0].plot(df['date_year'][df['station']==sta],df['east'][df['station']==sta]);ax[0].grid(True);ax[0].set_ylabel('East (mm)')\n",
    "    ax[1].plot(df['date_year'][df['station']==sta],df['north'][df['station']==sta]);ax[1].grid(True);ax[1].set_ylabel('North (mm)')\n",
    "    ax[2].plot(df['date_year'][df['station']==sta],df['up'][df['station']==sta]);ax[2].grid(True);ax[2].set_ylabel('Up (mm)')\n",
    "    ax[2].set_xlabel('Time (years)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 391
    },
    "id": "DcLbwDKqxHaF",
    "outputId": "f8d45657-3637-48e4-9891-c9fedc6ff4f8"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(32.0, 39.0)"
      ]
     },
     "execution_count": 45,
     "metadata": {
      "tags": []
     },
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAFlCAYAAAD226FvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd1gUV9sG8HuoAhaKqIiiYFdM7J0Ido0aNXaTqCm8RqPRL9GosSRGRRNbYmJiib2X2DVgw6iJwd6NICrYFSsgKOz9/bGwsmwFFtT1+V3XuYDZmTOzh5l55sycOUchCSGEEEJYD5sXvQFCCCGEsCwJ7kIIIYSVkeAuhBBCWBkJ7kIIIYSVkeAuhBBCWBkJ7kIIIYSVsXvRG5CZq6sry5Yt+6I3w+olJCTAxcXlRW+GVZMyzn1SxnlDyjn3HTly5C5JT0vl99IF96JFi+Lw4cMvejOsXnh4OAIDA1/0Zlg1KePcJ2WcN6Scc5+iKFcsmZ/clhdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjAR3IYQQwspIcBdCCCGsjMngrihKPkVRIhRFOaEoyhlFUb5Nm75PUZTjaem6oigbDCzfW1GUyLTU29JfQAghhBDa7MyYJxlAE5LxiqLYA9ivKMp2kgHpMyiKsg7AxswLKoriDmAsgFoACOCIoiibSN63zOYLIYQQIjOTNXeqxaf9aZ+WmP65oigFATQBoK/m3hLADpL30gL6DgCtcrzVQgghhDDIrGfuiqLYKopyHMBtqIP1vxk+7gBgF8lHehb1BhCb4e+radOEEEIIkUvMuS0PkqkAqimK4gpgvaIo/iRPp33cA8C8nGyEoijBAIIBwNPTE+Hh4TnJTpghPj5eyjmXSRnnPinjvCHl/OoxK7inI/lAUZQ9UN9aP60oSmEAdQB0NLDINQCBGf4uASBcT75zAMwBgAoVKjAwMDDzLMLCwsPDIeWcu6SMc5+Ucd6Qcn71mNNa3jOtxg5FUZwANAdwPu3jzgC2kEwysHgogBaKorgpiuIGoEXaNCGEEELkEnOeuXsB2KMoykkAh6B+5r4l7bPuAFZknFlRlFqKoswDAJL3AHyXttwhAOPSpgkhhBAil5i8LU/yJIDqBj4L1DPtMICPM/w9H8D87G+iEEIIIbJCeqgTQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKSHAXQgghrIwEdyGEEMLKmAzuiqLkUxQlQlGUE4qinFEU5du06YqiKBMURbmgKMo5RVEGGVg+VVGU42lpk6W/gBBCCCG02ZkxTzKAJiTjFUWxB7BfUZTtACoBKAmgIkmVoihFDCz/hGQ1C22vEEIIIUwwGdxJEkB82p/2aYkAPgXQk6Qqbb7bubWRQgghhDCfoo7dJmZSFFsARwCUBfALya8URYkDMA1ARwB3AAwiGaln2RQAxwGkAJhEcoOeeYIBBAOAp6dnzdWrV2f/GwmzxMfHI3/+/C96M6yalHHukzLOG1LOuS8oKOgIyVqWys+c2/IgmQqgmqIorgDWK4riD8ARQBLJWoqidAIwH0CAnsVLkbymKIofgN2KopwieTFT/nMAzAGAChUqMDAwMPvfSJglPDwcUs65S8o490kZ5w0p51dPllrLk3wAYA+AVgCuAvgj7aP1AN4wsMy1tJ/RAMIBVM/mtgohhBDCDOa0lvdMq7FDURQnAM0BnAewAUBQ2myNAVzQs6yboiiOab8XBtAQwFnLbLoQQggh9DHntrwXgEVpz91tAKwmuUVRlP0AlimKMgTqBncfA4CiKLUA9CP5MdQt6mcriqJKW3YSSQnuQgghRC4yp7X8Sei5lZ52i/5tPdMPIy3Qk/wbQNWcb6YQQgghzCU91AkhhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuAshhBBWRoK7EEIIYWUkuItcd/DgQfzwww84dOgQUlJSXvTmCCGE1ZPgLnJd3bp1ce3aNdSpUwfu7u54++23JdgDOHnyJG7cuPGiN0MIYYXsXvQGCOunKAqmTZuGuLg4LF26FNu2bcO2bdsAAAUKFEBAQAACAwMRGBiI6tWrw87u9dgt3d3d0bBhQwBAgwYN0LBhQzRo0AD+/v6wtbV9wVsnhHiVvR5nUfHC2djYYP78+bh3754msAPA48ePX9tgX6JECezatQsBAQFYtmwZli1bBkBdBvXq1dME/Lp166JgwYIveGuFEK8S6zxripeSvb091qxZg+bNm+Pvv//WO4++YN+tWzf8/PPPcHR0zPE2PHv2DFOnToWHhwfq16+PSpUqvdBasq+vL3bt2oW33noLt2/fBqAugx07dmDHjh0A1BdGVatW1ardly5dGoqivLDtFkK83OSZu8hTzs7O2LJlC/z9/c2av2fPnvjhhx8sEtgB9QVGr169MHHiRFStWhVubm5o3rw5Ro8eja1btyIuLs4i61myZAliYmKwcuVKXLt2zei8FSpUwK+//ooCBQro/VylUuHEiRP49ddf8d5778HPzw/e3t4YNGgQnj17ZpHtFUJYF6m5izzn5uaG0NBQNGzYEJcvXzY43zfffIOxY8dafP0lS5bE7t270bhxY8TGxmLnzp3YuXOn5vNy5cqhfv36qFevHurXrw9/f/8sPxp47733sHjxYvTo0QOAuoYeEBCARo0aISAgAHfu3IGPjw/y58+PZ8+ewdXVNUvraNq0Kb755hvY29tnabuEEK8Jki9VKl++PEXu27Nnz4veBEZGRrJIkSIEYDB16dKFMTExubZ+Ly8vo+sHQBcXFwYGBnL48OHcuHEjb926ZVb+e/bsYbdu3fTm6erqSkdHR7q4uLBw4cJ0dHRk0aJFTW4LAI4aNSpXyuNV9DLsx68DKefcB+AwLRhL5ba8eGHKli2LP//802hjsTVr1qBixYoICQlBcnKyxde/e/duFClSxOh8CQkJCA8Px6RJk/D555/j3LlzZq9j8eLFaNOmjc70Bw8eIDk5GQkJCbh79y6Sk5Nx69Yts/IcP348ChUqhODgYPz1119ISkoye3uEEK8HCe7ihapevTo2bdqk9Uy9ZcuWyJ8/v+bvxMREjBw5Ev7+/ti+fbtF11+xYkXs3LkTHh4eJuf95JNPcPLkSTRu3Njs/B0cHLB27dosLWOOR48eYe7cuWjcuDGcnJxQoUIFxMbGWnQdQohXlwR38cI1btwYq1atgo2NenccMGAALly4gPfee09rvqioKLRp0wYdOnTApUuXLLb+qlWrIiwsDK6urgbncXR0hI+PT7Za1js5OWHTpk2oXbt2TjbTIEVR0Lp1a3h7e+dK/kKIV48Ed/FSeOeddzBv3jwAgIeHB7y8vLBkyRL89ddfeOONN7Tm3bhxIypVqoRvvvkGT548scj6a9SogdDQUIMt1pOTkzF69GiUL18eixYtgkqlylL+BQsWxPbt21GlShVLbK4Wkvjxxx9RtWpV7N692+L5CyFePRLcxUujb9+++P7777VukQcEBODIkSOYOXMmChUqpJmenJyMb7/9FpUrV8aGDRugbo+SM3Xq1MH27dvh4uJicJ5r166hT58+qFWrFvbs2ZOl/D08PBAWFgY/P7+cbqpeZ8+eRdOmTVGwYEG0a9cOkyZNwvbt2zXvzwshXh8S3MVLZejQoShXrpzWNDs7O3z22We4cOECPvroI63PLl++jI4dO6JNmza4cOFCjtffsGFDbN68Gfny5dNMW7FiBerUqaM137Fjx9CkSRO0a9cO58+fN5jfkydPsHfvXjx69AgAULx4cezcuRPFixfP8bYa8vjxY2zZsgUjRozA6NGjER8fn2vrEkK8nEwGd0VR8imKEqEoyglFUc4oivJt2nRFUZQJiqJcUBTlnKIogwws31tRlMi01NvSX0BYn/Rn75kVKVIE8+bNw8GDB1GzZk2tz/7880/4+/tjxIgROQ5mQUFB2LhxIxwcHACoG/j9888/WL58OUqVKqU1b3qHPAMGDMCdO3d08nJyckJUVBTc3NxQsWJFvPfee9iwYQMmTZoEd3d3zXxTpkzRuYCwhNjYWOzatQupqalZWu7atWs4ffq0xbdHCJFHTL0rB0ABkD/td3sA/wKoB6AvgMUAbNI+K6JnWXcA0Wk/3dJ+dzO2PnnPPW+86u+tnjx5kuPGjaO7u7vOe+De3t5cuXIlVSpVjtaxefNm2tvbMzk5WTPtyZMnnDRpEgsWLKiz3gIFCjAkJIRPnjwhqV3Gv//+OxVF0VnGxsaGADhs2DDu27ePc+bMMevd+6ymQoUKsX79+mzVqhWbNm3KgIAA1q1bl9WrV+eAAQMYGRnJ+fPns0+fPvTz86OiKDx69GiOyi8vvOr78atCyjn3wcLvuWdtZsAZwFEAdQFEAChrYv4eAGZn+Hs2gB7GlpHgnjde9YP10aNHrFOnDj08PFi6dGm9AS0oKIinT5/O0Xo2bdqkd/rt27c5YMAA2tra6qzXx8eHy5Yt0ynj+fPn6w3wmYN95cqV+cYbb+jNOzeCfrFixXSmf/TRRzkqt7zyqu/Hrwop59xn6eCu0IyGSIqi2AI4AqAsgF9IfqUoShyAaQA6ArgDYBDJyEzLfQkgH8nxaX+PBvCE5JRM8wUDCAYAT0/PmqtXrza5TSJn4uPjtd4lfxWlpqbiv//+M9livmjRovDy8sqVAWKSkpJw9epVPHz4UOczHx8fODk5aZVzXFyc0S53M1MUxSKNBbNCURSULl0aycnJKFKkSK4NrKNSqUAyR/lbw378KpByzn1BQUFHSNayWIZZuRIA4ApgDwB/APEAvkib3gnAPj3zfwlgVIa/RwP40tg6pOaeN6zlSvzmzZssX768yRqqk5MTg4ODefHixSzfrj979iwHDx7MX3/9lXv37uWdO3d05tm9ezerV6+utc4pU6YQADt27MgLFy5o5l20aJHJGjwAtmrVijdu3OCsWbNYqlSpXKm5OxqYHpT20xfgZkdHcuBA8sEDxsXF8f79+2aVW2JiIs+dO8ft27dz1qxZHDZsGLt06cLatWuzcOHCLFq0KG/evJml/0Vme/bs4Z07dxgWFsZJkyaxW7duDA4OzvEjGaHNWs4XLzNYuOaepdEwSD5QFGUPgFYArgL4I+2j9QAW6FnkGoDADH+XABCelXUKYUzRokWxc+dONGrUCDExMQbne/LkCebMmYM5c+bA09MT77zzDmbMmGH0tbd0lSpVwltvvYXu3bvj6dOnAABPT09UqVIFlStX1qStW7ciLCwMX3/9tdZIcOvXr8fmzZvRv39/jBkzBh988AEURUHv3r2N1srj4uKwatUqBAQE4Pz585g/fz6++uori7Z+19ehrxeAflBfxb8BYEtyMqbMnIlzixYhn6srTp06pV42ORkxMTG4fPkyLl26hMuXL2v9fvPmTaPrDg0NRdGiRdV/PHwIFCwIGBnGliSuXr2Ko0eP4tixYzh27BiaNWuGoKAgzTyOjo44efKkDIcrhKnoD8ATgGva704A9gFoC2ASgA/TpgcCOKRnWXcAl6BuTOeW9ru7sfVJzT1v5PWVeGpqKvft28fExMRcyT8yMlLvs2N96S0bG15ZvjzL69i1axfz589vNG93d3c2aNCAtWvX5tSpU3U+d3V15ZQpU5iUlMQlS5ZoGtSZSs7OzgwMDOTgwYNZrVq1XH0Ob2/ks927d2vKIyIigm3bts3WOj7++GPt2vXcueTx45o/U1NT+d9//3HFihUcNmwYmzdvTg8PD5180u+OpKfvvvsuR/uR0E9q7rkPed2gDuqL92MATgI4DWAMn9+i3wrgFIB/ALyZNr0WgHkZlv8QQFRa6mtqfRLc88aLOFj/+OMPuri4sHnz5vz+++95/PhxpqamWiz/U6dO6W09nzEFA1QB6tvM2RAREaE3yOhLmQNPxlTLx4fbZ8/msmXLzA7wmVN2l8tu6gow4e5dnTI5cuQIO3TokOX8ChUqxAYNGvDjjz9mlL8/k0aPZmpqKhcsWGBytEBDZdywYUP26tWLQ4cO5Y9TpvDCmjU53a0EJbjnhTwP7nmdJLjnjRd1sM6cOVPrZFykSBH26tWLCxcu5LVr13Kc/6FDh1igQAGjAaGHgwNj5s7N9jrOnTvHEiVKmBV4OnTowNWrVzNfvnya6Z4AxwLcamPDuIULuXz58iwHal9fX27fvp0VKlTQ+7k3wMq5EODz2dvT39+f7du35zcffshLS5ZoyuXEiRPs0qWLWe0JMiZHgPEAU+rU0eQVFxfHadOmsWzZslkK7unJFeB5QH2K+/zzHO1TQoJ7XpDgLiziRR6sQ4cONXiyrlKlCgcPHsxt27YxPj5e7/IqlYqLFy/mxo0beePGDZ3P//rrLzo5ORkNCk5OThwzZozBdZhy5coVg4E1Pc2YMYN///03STI6OtrgeO2+vr6sW7euJig2bdrUZEBs164dIyMjef/+fc1t+vQLiFoAE9IC2zAT+ThkM8hnXAdDQrTK5vTp0+zRo4fZQT4/wPbFi6tvy2dqCJeamsqwsDB26NBB7wWQoeBeKn3bALJ+/Wz9j8VzEtxznwR3YREv8mBNTU1l9+7dTZ70HRwcGBQUxIkTJ/Lw4cNat/CvXr2qqdX5+PiwS5cu/OGHH/jXX38xISGBoaGhdHBw0OTl7++vN+B7e3tz6dKlvHHjBpctW5al1vS3b99mzZo1TdYqmzVrxr179/LSpUsG38nPnJydnU3OY2Njw08++YT9+vWjq6sru3fvziZNmrB5hsA220Qe9gAHAayfxeDeNmPwDA7WWz7nzp3j+++/rxOUS5QowS1btnDy5Mns3bs3a9WqxW7dupks75iYGI4aNUrrIiljcC9atCj79OnDFi1asFrlypyeLx8TfX3JzZvN2zGFQRLcc58Ed2ERL/pgTUpKYuPGjc0OJn5+fty3b59WHrGxsfTz89OZ19bWltWqVWOLFi00gSUkJIQxMTHs2bOn3vzr1q3LIUOGEFA/KnjnnXcYEhLC8PBwo7X7hw8fMigoyGhwT09vvfUWly1bpnebc5IcHBxYpkwZrWkejo5sUrgwOwQGGm0jMC1DkG5pZB22isJVy5fz559/5pAhQ9j17bc5r2hRqjp2JKOijP6vIyMj2bdvX61OeaJMLGNMcnIyV65cybfeekurjO3s7Hjr1q1s5ysMe9Hni9eBBHdhES/DwXr//n1WrlzZeM3S3p6jRo0y2Mr+ypUrZtWGHR0dGRQUxOHDhzMkJETnnfT0pG+6ra0ta9SowQEDBnDJkiWMiorSqt0/efJEp0FZr169OG3aNIPrMKeL2W7duvGtt97SCuINGjQwWlZZvTCYniG4tzYxb2WAhzw9mbp0abb+39HR0QwODqa9vT2nTZuWrTwy2759O/v37695i8FS+QptL8P5wtpJcBcW8bIcrFeuXDEa6IoVK8b169cbzePSpUv08fExO6AVLFiQCxYs4KJFi1i8ePFsBckiRYqwffv2DAkJ4Z49e/jgwQP26dNH8/mZM2e4c+dOjh49moUKFTJ4wWFsHa1ateKBAwfYsGFDAurbziqVips2bTL79r6p5ADwS4Dts7Gsvb09165dm63/+ZIMDfFyIn0/fvToEWfNmsWePXtKBza54GU5X1gzCe7CIl6mg/XYsWMmW7i3b9+eV65cMZhHVFQUvb29TQakoKAgrXzi4+M5evRordbs2UklSpTgyZMn+X//938EwFOnTmnK+P79+xw3bhzd3Nyylfe7777LN954g2XKlNFs99OnT/nTTz/pHcDG3JTTV+nmz5+f4/99Tl+FzLwfq1QqCe654GU6X1grCe7CIl62gzUsLIx2dnZGg4mLiwunTp3KZ8+e6c3jwoULJm93u7m5sWPHjpw0aRJDQ0N5+/ZtkurapDmN/PSl7t27827a+98qlYoTJ07kiRMndMr40aNHnDRpEj09PbO8DltbW1aoUIHXr1/XyvPevXscMmSIybLTl8aMGcN58+aZ7JjHUPLz8+MXvXrxUKbb9MnJyYyNjeWhQ4e4ZcsWzps3j+PHj+fAgQPZpUsXBgQEsHz58ixZsqRWt7zZ8bLtx9ZKyjn3SXAXFvEyHqyLFi3SBI6mTZty3rx5emu71apV47///qs3j/Pnzxt85cxQSu9MZfjw4ZwwYQKrVq1q9rJly5bVG6BSU1MNlnF8fDynTZtmsEc9Dw8PtmrVSu9teycnJ44YMUKnf/fIyEh27NgxywG6VKlSnDlzJgcMGJDl99PtAb4BsAnA2h4erFy5stkd/ADghg0bcrzPvIz7sTWScs59EtyFRbysB+v48eMJgO+//z5J8tatW3z//fd1AoOiKOzfvz8fPHigk8eZM2fM7uHMUKpTpw4LFy5s1rx2dnb8/PPPGRcXp7Udpso4MTGRM2fONNghTokSJdiwYUO9t8/d3Nw4efJkJiQkaOW5d+9e1qhRI8vft3Tp0pw8eTJr1aqVo3IzN33xxRc521HMLGNhGVLOuU+Cu7CIl/VgValUDA4O5tChQ7Wm79q1i+XKldMJEsWKFeOqVat0nrOeOnVKKzi3adMmW0HI0dHR7Bqtq6srp06dyqSkJJLml3FSUhJnz57NkiVL6s23aNGiBvuTL168OGfPns2nT59q8ktNTeXixYv1NhY0ldzd3dmhQwfNe/aF9MyjWCC429jYsGrVquzbty9//vlnHjx4kE+ePMny/vKy7sfWRso590lwFxbxMh+sz54946FDh3SmP3nyhGPHjtXqnCY9tWrVilFRUVrvOZ84cULT1/yOHTv466+/5kmt1M/Pj2vWrMlyGT99+pTTp0832Lueh4eHwe5Yy5Urx1WrVmk1UEtISOC4cePM6hDHYBAGaJfp4qZWrVqsX78+S/n40C7T/LZQd/2anXWVtLXlmmLFqAoLM7vMXub92JpIOec+Ce7CIl7lg/X8+fN6O47Jly8fO3TowBo1avC7777jyZMneeTIEbq5uXH//v0kyd9//z3Lz5aNpdIAqxv47LfffuPBgwez/P3i4+MN1tQrAOxmb0//TJ3WpKcaNWowNDRU607GtWvX+OGHH+b6RU2hQoXoDPAuwK8A5stiOYcD6lMSQF66ZFZZvcr7cXalpKRo3anJC69jOec1Ce7CIl71gzW9f3l9z8Uzjgzn6+vLHj16cOnSpZpW9lkZatVgjTZt+UNpwSgUoH/FilrzpPee1qNHD16+fDlL3+/Jkyds1qyZ0W0oWLAgXVxc9H4WFBSkc2Fx7NgxzTN1WwsH9hkzZlClUvHQoUOc/cknTB04kFcjIhgcHKzVMx0ANmnShHPmzOH//vc/1qxZU9OvwPb0wG5vT169alY5WXo/fhlfo0tISOCCBQvYtm1bFi1alPnz5+esWbPydBte9fPFq0CCu7AIazlY7969y48++sisAOTu7s4PPviA69at4+LFi7VeH8tql7CVK1fm//73P26xsVEfRm5uTLl9m3PmzNG01s/YNaqjoyO/+uorvQ0ADUlISGBgYGCOgm7Hjh15/PhxTpo0iSEhIVy1ahWnTp3KVh4edLNgcC9QoAC9vb1ZuXJlVqtWjbVq1WK9evW4dOlSnj9/np07d9bMW6JECa0gmpSUxEOHDnHujz8yeuBA8tw5s8sop/vx1atXuWLFCn766ads1KiRTqPIF+HOnTucOnUqg4KCDA5hrO+xVW6ylvPFy0yCu7AIaztY//rrL5Nd2WZMjo6OrFmzpqZWOXv2bJ49ezZLLc3t7OzYt3dv3vrtN/LOHc22PHr0iKNGjdLb/WzhwoX5yy+/GHxXP7PHjx9reqgzlby9vfV2amNjY8NP27fX6Vq3oIsLHbLxfry5qU+fPlptACIiItikSRMC4NmzZy3yf8/KfqxSqRgVFcX58+ezT58+Ohd0M2bMsMg2ZdXFixc5evRo1qlTx6w+B9q0aZPn22ht54uXkQR3YRHWeLAmJydzwoQJ2eptzs/PjyEhIbx48SLDw8Oz1L1r/vz5OW7cOJ0BZsLCwvjBBx/oXaZixYrcvHmzWbeBHz58yDp16hjdBn9/fz5+/Jj379/nyJEjtRrRfQXwKcA7ZcqwQ+vWuRbM9ZVLjx49uG7dOs0reyqVimFhYdy1a5dZ/88rV67w4MGD3Lx5M1NSUnTmMbYfp6am8tSpU5w1axa7d+9u9O2BChUq5Nlz7FOnTnHQoEH09/c32QWxvmRvb8/8+fPT3d2dXl5e9PX1ZeXKlVmzZk3+888/ubLN1ni+eNlIcBcWYc0Ha1RUFFu2bJmlE6aNjY3Wu+oqlYorVqygq6ur2Xl4eXlxzpw5mlp5ehkfOXLE4O31Jk2a8NixYya/07179wwOdpOefH19Na8FXr9+nf3796ednR1DAfWhDjDl+HEOGjQozwJ8enJ2duakhg1548gRJiQkMCoqivv27ePq1as5Y8YMfvXVV/zggw/YvHlz+vv763SGM3nyZL3lknE/Tn/LYurUqXznnfBRjMUAACAASURBVHcM3tLWl3x8fNi6dWt269aNn3zyCb/44gtu2bIlZzuiHiqVitu3b2fjxo112iJYIk2fPt3i20xa9/niZSHBXViEtR+sKpWKK1euNNgLXMZUvXp1njPwnPfZs2ccPnx4lk6wlSpV4saNG7XKWKVSccOGDSxfvrzO/IqisE+fPrxqohHZnTt3zOo9r169ejxw4ABJ9YXO9EaN+LRoUbJPH/LxY5Lk9OnTLfrWgKk0K+3i4jHAsllctkOHDgbvcGQs4+joaI4aNcpig+pk7mvBknbv3s0mTZrkuGFneuoDcAvAh7a2pLu7xbfX2s8XLwMJ7sIiXpeD9f79++zXr5/WidAB4DCA76X9bWNjw+HDh2s6n9Hn6NGjdHFxyVJtS9+rcOkDvujrptXZ2Zljx441On78rVu3WDFDq/zChQsbvIDp3Lnz83HT9QzQsm7dOpOPMEqUKKH1nbNb29yI53cPamVhOVtbWy5fvpw3b97UWx769uPU1FTu27eP//vf/7I8WI87wC8AfgDwRrFi5Pz55JQpWm0qLOnhw4dcs2YNu3XrZrB/A3NSxqF7qSgW387X5XzxIklwFxbxuh2s4eHhmpNnSIYTYdcMJ8jKlSsbbYX833//8ezZs2YPMJPeWr5Lly6MjIzUyuv+/fv88ssvdTrk6QFwfYECfHb+vMHtuHbtmqYzm9q1a/Off/4xGBjs7Ow4cOBAg63A//77b7O62U2v5VeoUIEhISFZ7r+/PMC1AD9D9l/D8/b2Zrt27Th27Fhu3LiRV69eNbkfJyUl8Y8//mCnTp30dn6UnoKCgtizZ0+2a9eOl5ydmZoxUFaqZHQdlpKSksIDBw4wODjY6OOgX3/9VXdhLy8J7q84Ce7CIl7Hg/Xy5cssUKAAx2UI7p0z3Ra1tbXl119/bbQWT6r7rw8MDGQzM4J7epD97LPPtHrQI9Utpbt27aoJgJoTdFCQ0fXHxMTQ19eXgYGBJNV3Fow9Y7a1tWXNmjU5efJk7tq1SyvYR0ZGGuz5Tl+QHzBgAM+cOcPZs2ebtVxZqGvu2wEWS5tWD+BCgKWyGejtoL77sn3KFEYOG2bW///evXucM2cOAwICdPLr2rXr8xmHDXv+fwDUf78A0dHRHDRokM7/1cvLS3f/lOD+ypPgLizidT1YDx48yHz29gwGuOmTT3js2DG++eabOif7qlWr8siRIybz+zcsjMf692d7PT3mpQd3hwzTChQowPHjx+vcej9w4ADbV6/OREVRH5Zduphc96VLl9i/f3/N36dOncrSgDk+Pj585513+Oeff/LOnTusX7++5jN/f3/Onz/f4OuFNjY27Ny5Mw8cOMA1a9ZoOsdxBjgYYMEM8/6SIVD+26QJAwICWBVgUlr6qVEjuhl4/m/osUHFtPz2TJnCxzY2Wd4PLl26xAkTJmgecdja2jImJkb94b592sH9r7+ynL+l3bt3jwMHDqS7uzudnJx0a+8S3F95EtyFRbzOB+uSJUsIgJs3byapfuVqzJgxOs+T7ezsOHbsWCYnJ5uV78GDB7V6lZsyZQo/Angf4Kq0wJex9jV37lyt991VKhVjDxwgFy0i9bz2pU/m17fOnTtnckz7jLXwsWPHarYhMTGR7777rjp4Vqyo2aZt27Zp3k/Xlxo1asT169dz586dbBUUxEsAVQDX1ajBKlWq8Evg+W3u+fNJkjdv3uTeefPI5ctJlYrbWrTg+EwXBQA4f/58RkdHc+3atRw5ciRbtWpFT09PegG8lRbcr3l6mlVW+qhUKh4+fJiDBw9+3tL82TN1ozRA/dPMPgksKTExkceOHdPb6VFKSgovXryoPVGC+ytPgruwiNf9YB02bBhPnz6tNe3IkSOsVKmSTvCqVq0ajx8/bnbeu3fvZv369TllyhT+laEG6KcnMFauXJmbNm2yaLenkZGRBkeY03c7Ojo6WrNsamoqv/jiC5YqVUon38OHD7NHjx4GG9WVL1+es2fP5sUzZ8g//yTTvlNsbCxjJkwgjdwJeZqQwDMLF/L4H39w2LBhmhp79+7ddeZVqVSMjY3l4rFjGbp0KRMt1Kuc1v/gvffU/7f33rNI3oYkJCTwyJEjXLx4MYcPH8727duzbNmytLGxYePGjbU6ATJKgvsrT4K7sIjX/WA1NPhGfHy83gBvb2/PcePGmd3RSfqrb9OKF2eqgwMZFMR9W7eydu3aujVogMOrVOHFEycs9v0uXbpEX19fswK8g4MDv/jiC967d0+z/LJlyzS/7969m2vXruXevXt55swZHjlyhJ9//rnB3tQ8PT05btw43slBC/Nr166xX79+9PLyMhrgcm0/XrVKfXpcudIi2T1+/JgRERFcuHAhhw4dyrfffpu+vr4GX0d0cXHRrZ0bI8H9lSfBXViEHKyGJSYmGhyVrVy5cmY9iyfVZaxSqbReo1KpVFy1apVW16eb007KqcWKkWk9uVlCTEwMy5UrZ1aAB9R978+YMUPnMcTDhw/ZsWNHnVv67u7u9PT0NNgKvZy9PX9o0oSxMTEmGygaEhkZyYcPHxr8PNf24wcPSGdn9c9sCgsLY+vWrenj42P2/yA9DRkyhI8ePTJ/ZRLcX3kS3IVFyMFq3P379412GOPn58dvv/2WBw4cMFibN1bGycnJmvfdT2e4dc8MtWdT0l+dun//vsF5rl+/rnMnwthrVuUA/uLmxrhTp7TyUalUnDx5crY6XXmzcGFev37d7O+VFbm6H4eE5DiLU6dO8cMPPzT6Gp6xVLhwYdauXZsrVqwwviIJ7q88Swd3GwghdLi6uiI0NBS+vr56P4+OjsbYsWPRsGFDuLm5oVWrVpg8eTIiIiKQkpJiMn8HBwcMHDgQFy9exNGPPkJqYCDw+++Aq6vZ22hra4uHDx/Cy8sLZcqUQZcuXRASEoLQ0FDcvn0bAODl5YXw8HBUrVpVs9zWrVsxe/ZseHt7a+VnA+A0gP7378O5VSutzxRFwbBhw7Bjxw4ULlzY7G0MBLC3WTN4eXmZvcxLY8iQHGfh7++P33//HTExMRgzZkyWyg4A7t69i0OHDiEhISHH2yJeM5a8UrBEkpp73pArcfNcuHCBnp6eZte0mjZtqmmglldlHBoaqveVsRIlSrB9+/YcO3YslyxZQn9/fwLg0aNHSarHjJ8+fbqmExsbqAeYIcAogG3bttXb731MTIzJgWwA0Nnenns7dKAq07v9lvSq7ceJiYmcPXu2Vi+DxlIPqN80UCmK8Rq51NxfeZDb8sIS5GA13+HDh00OxWlnZ8eQkBCtFtd5Wca7du3SGgnO2HY2aNCAI0aM4Jo1a3jx4kU+evSIEyZMYKFChVgB6h78imdYpmvXrjyfqce8pKQknW59DaWqVavyxx9/zJWx0l/V/Tg1NZVbt25l06ZNdcqrYMGCPHjwIFevXs1N3bs/D9qA4QwluL/yJLgLi5CDNWt27dpl8rmpvb09p06dqhmaNK/LeO/evXRxcTH7LoOrqytXr16tWf7evXv8+uuv9eZhY2PDPn368NKlS1rrXLBggdlD7Do6OrJHjx7ctWuX+a94mWAN+/Hx48f5wQcf0N7eXlNWoaGh6g+XLZPg/pqQ4C4sQg7WrFuzZo1ZI6nVqlWLFy5ceCFlvH//fhYoUMDkNgYEBPDKlSt687h16xaHDBmid6xxe3t79u/fX6uB3NGjR3VeuytbtqzRux1+fn6cMGECr127pv+LmPnevzXtx9euXePIkSPp5ubGPn36qCdKcH9tWDq4S4M6IczUuXNn/PLLLybnO3z4MAIqVULctWt5sFXaGjZsiLCwMBQsWNDofFWrVoWjo6Pez4oUKYJp06YhKioKwcHBsLOz03z27NkzzJo1C35+fhg6dCju3r2L6tWr4/Dhw2jdurVmvpYtW+LGjRv4/fffUb9+fZ11REdH4+uvv0bJkiXRvn17bNy4Ec+ePXs+w/r1Wfzmr77ixYtjwoQJiI2NRbNmzdS1LyGyy5JXCpZIUnPPG3Ilnn1jx47V1EC7d+/OGTNm6L2V/fOUKQwdNeqFbOOhQ4eMvvIGqIeYHTlypNFX6Uj1mPDvv/++3rsWQc7OXD9wIEn1c+SxY8dSURS+++67WnmcPn2aQ4YM0TvUbXoqVqwYhw8fzqhTp0gfH7O64H1Z9uPHjx9z586dXLp0qWUzlpr7awNyW15Yghys2adSqfjpp58SAAcPHqzpmKZDhw5awSp94Jh3331Xcwv84cOHnDBhAjds2MDbt2/n6naaGikuPbm6ujIkJIQJJjrQOX36NDt16qRZrhMy9BmfYSCTrVu3sl27dnrzSEpK4qpVq9i8eXOD29MmLc8UMwZseVH78dWrV7lq1SoOGjSINWrU0HTJu3HjRsuuSIL7ayPPgzuAfAAiAJwAcAbAt2nTFwK4BOB4WqpmYPnUDPNsMrU+Ce55Qw5WdQcjffv25bx58xgVFZWl/t1TUlLYuXNnfvfddyTJf//9l/nz56ePjw8LpD1rzjjkq7OzMydOnMikpCTu379f07K9fPny7NOnD+fOncuzZ89arKFZuuPHj2uN1/7555+zffv2BmvOv/zyi8mBcg4fPsxWrVoxOGPQGTdOa557ZnTGEx0dzdGjR9Pb21trO0YCjLe3J2fMMJlHXuzHqampPHnyJH/99Vf26tWLpUuX1lt+LVq0sOgYASQluL9GXkRwVwDkT/vdHsC/AOqlBffOZiwfn5UNkuCeN+RgVVu+fLnmdnOJEiX43nvvmR3sk5KSePjwYc3fYWFhWi2eMwb39FS+fHmGhoZy586dehusubu78+233+bEiRMZHh5usjZtjoxDwc5PG5Xt77//ZmBgoN4g5evry8WLF2ta/RtyYMcORr/3HjlkSI66aU1JSeHWrVvZsWNH2tnZEQC3//GHWY3qcmM/TkhIYHh4OMePH8/WrVubfLyRnsaNG8ewsDBGREQwMjKSd+7c0Rr1L1skuL82XuhteQDOAI4CqCvB/dUmB+tzM2fO1HuyzmqwJ8mVK1dqLhb0Bff0VK9ePf7000+aYGYo2dnZsXbt2hw8eDDXrFljuHW5CWfPnmWxYsW4ZMkSzTSVSsWwsDDNWOyZU5UqVbh+/XrL10aNuHHjBmfMmGHywiKdJfdjlUrFZcuWmT1krrlp5syZ2d8oCe6vjRcS3AHYpt1WjwcwOW3aQgD/ATgJYDoARwPLpgA4DOAggA6m1iXBPW/IwaotYyM5Q8ncYJ9+sWAsuKcnc2uF6al48eKMiIjI1nf877//+Jee59gqlYrr1q3TOxoeANatW5e7du3K1jpzW27sx8nJyVy/fj07dOhg8uIrc3IE6JT2uzPAlgDD/u//yPDw5yvISmc+EtxfG5YO7oo6T/MoiuIKYD2AgQDiANwE4ABgDoCLJMfpWcab5DVFUfwA7AbQlOTFTPMEAwgGAE9Pz5qrV682e5tE9sTHxyN//vwvejNeKrGxsZo+2c1RqFAh+Pn5wcZG943S69evw9bWFlevXjWZj6IoMOc49PDwQMmSJWFra2v2NmZVXFwcrl+/jqdPn+p8VrBgQXh7e8PZ2TnX1p9Vub0fp6Sk4N69e4iLi0NiYqJZy1QBcBfAYwCV0ic6OAAVKwI3bgAqFVC6tHkbcO8ecOnS879r1tQ/38mTQMZXCQ3Nl01yvsh9QUFBR0jWsliGWb0aADAGwJeZpgUC2GLGsgth4la+1NzzhlyJ60pNTWW3bt1M1s4cHBz43XffGW14plKpuGjRIovd2i1RogT/+eefPCmHpKQk/vTTT5rn9JlTx44defr06TzZFlPycj8+ffo0hw0bxuLFi+uUib29PWNiYnjr1i0+bNeOBPjMyel5bdreXp0AculSpqam8sSJE6aHwpWa+2sDed2JjaIonmk1diiK4gSgOYDziqJ4pU1TAHSAekCpzMu6KYrimPZ7YQANAZw1tU4hXgQbGxssXrwYLVq0MDpf7969MWzYMDg4OBicR1EU+Pj4oHPnzhbZtqtXr6J+/fro1asXYmJiLJKnIY6Ojhg4cCCio6MxYcIEFCpUSOvz9evXo2rVqujduzcuZaxVWrkqVapg8uTJiImJwZ9//okePXogX758ANSd+yxatAhFihRBwV69AAB2T55ollU9e4aTz57hJ0VBpxUr4OnpiREjRhjsSEiIHDMV/QG8AeAY1M/WTwMYkzZ9N4BTadOW4nmL+loA5qX93iBtnhNpPz8ytT6puecNuRI37PHjx6xbt67RmnSFChWe9/9twJ49e5iUlMQmTZpYrAYPgPny5ePo0aP5+PHjPCmPuLg4Dh8+nE5OTjrb0tbWln/07Gnw+wcEBPCdd97hp59+yvHjx3PBggUMDQ3l6dOnef/+/Rw31nvR+/GDBw84d+5cNmrUiN7e3nz69Cl5/z5TbW15AuCPADsCdM9Ubra2tjxz5ozpFUjN/bUB6cRGWIIcrMbdvXvXYAOzjKlTp068fPmy3jzSy/jRo0esWbOm1nLmDiNbyMGBY0aP1htYvby8uGDBAou/G2/I9evX2b9/f02HLT0zBp1p0/Quc+jQIZ332DMmZ2dnli1blo0bN2aPHj345Zdfcvr06Vy9ejX379/Pu3fvGt2ml2U/Tk1N5ZYtWzhx4kR27NiR7iYa4r3//vtMSkoyfXEjwf21IcFdWIQcrKbFxMSwZMmSmhNyz549WaVKFZ0TtZOTE8eNG8cnT55oLZ+xjG/dusVy5coRAF1cXHj//n0OGjSINjY2+u8MpP2sCJC//sqYmBj26tVL77w1a9bU2wo+t0RERNDb25uFAXYBOAzgb+3bMywsjJGRkTptEa5fv27yToi+FBAQ8FIH9wsXLvDHH39UB3MzegLUl+zs7Ojq6soSJUqwYsWKvHr1qvZKJLi/Niwd3GXgGCEMKFmyJMLCwuDh4QEAaNWqFY4dO4bp06drDczy5MkTjBkzBlWqVMGWLVv05lWkSBGEhYWhePHi8PDwgKurK3788UccPXoUjRo10pn/AoDKUHcsAXt7lCxZEkuXLsXBgwd1BmI5cuQI3nrrLXTp0iVPnoHXrl0bkZGRqNGsGdYA+B5Av02b0KJFC5QrVw5OTk4oVaoUgoKC8NFHH+H3339Hv379tAaWMaVnz57YsWOHpuxfRt7e3ihatChSU1Px+PHj7GWSkoIHDx7g9tWreHr+PPJduwbEx6s/y9j6XYissuSVgiWS1NzzhlyJmy8iIoIuLi7cunWrZtqNGzfYu3dvvbWxt99+m5GRkXrL+NSpU2zatKnWNJVKxSVLlrBYsWI6eRV0cuKSJUu0bt+qVCquWLGCPj4+OvM7ODhw2LBhfPjwYY6+85IlS9i0aVP279+fP/30E0NDQ3nlyhWtRwCpqakcNmxYtmqsxlL+/Pm5bNkyk93gki/PfhwXF8c5c+YwMDDQrGGB09PXACcCbJaxdj5/PvnDD+Snn0rN/TUCuS0vLEEO1qzZsWMHjx49qjN9//79rFatmt4gu2zZMr3dx966dUvvOh4+fMghQ4ZonmlnTAEBATxx4oTW/ImJiRw/frzeEemKFCnCOXPmmN3Tmz7z58/X2RZnZ2dWq1aN3bt35zfffMMVK1Zw5MiRBh8vZE4uLi4cNWoUCxUqZHLeokWLctSoUYyNjTW4jVndj1UqFZOSkhgXF8crV67w7NmzPHToEM+ePZvtcsosNjaWU6ZMYY2iRQ1+ty+//JJr167ltpEjSYCpNjbPg7OtrfrnggUS3F8jEtyFRcjBajkpKSmcNWsW3dzctE7gU6ZMoY+PD9euXZulVuGnTp1i48aNdQKCra0tBw0apDNE6/Xr19m3b1+9NcY33niDO3fuzPZ327RpE/Ply2dW4Danxrpy5UqS5Pnz51m+fHmz8rW1tWWnTp24c+dOnXLcs2cPz58/z6+++or9+/dn7969+e6777JVq1YMCAhg9erVWb58eRYvXpyFChXS2+NcgQIFePLkyWyXkUF79/IcwDEAy2a6+AkKClLPk5JCenhoB/D0AH3rlgT314gEd2ERcrBa3p07dxgcHKy3b/lmzZplqXaoUqm4fPlyvf2cFylShAsXLtRpJX/kyBG9FwUA2L59e164cCFb32v//v06Fy45Sa1bt+bOnTt57949tmzZUuuzwMBAlilTxuCyFStW5E8//cQHaQPVpO/H+/bt4xtvvJHlbbG1tTX5SiNJs8aW1/HsGenmRgJUde3KiIgIDh48WPP45dixY+r53n9fN7jXq6f+TIL7a0OCu7AIOVhzT0REBOvUqaPTt7ydnR2HDh3KR48emZ3Xo0eP+OWXX+qtcTZo0OB5gEiT3k+8n5+fzvz29vYcMmSIWcOxZnb69Gmjr7SlpzJlyrBChQpmBdZq1apx4cKFHDx4sGba2LFjmZqaym3btvHtt982eDfAxcWF/fr147Zt2zTb+OzZM86cOdOsW/7pqW3btoyOjjb+5R8+JDMMuJMlvXqpT7NLl2ompaSkcOfOndy8ebN6wqpVusF9wgT1ZxLcXxsS3IVFyMGau1JTU7l27VqtsdTTk5eXF5ctW5alW/VnzpzR2xmOjY0NBwwYoBOwk5KS+P3337NgwYI6y3h4eHDmzJnqDley4MqVK6xYsaLRYOnt7c2QkBA2a9ZMZ7qDg4PeZYoXL84uXbrQwcGBgwYN0lpndHQ0hw0bRg8PD73LTpkyhQEBAVy5cqWmAd6tEyfYt2/fLNXgS5Uqxd69e3PhwoW6/RasWUN26JClstJYuZK0sTE+WMyDB6SdnXYQT39MIMH9tSHBXViEHKy5b8+ePbx37x4/++wzvQ3O3nrrLZ1GcsaoVCquWrVKbw3a09OTv//+u86t+lu3brFfv3466+8LcF3x4lTdvp2l73T37l3Wq1fPZLB0d3fXGko2ODiYN27c4Ndff23wnfB8+fKxZs2aemvST5484cKFC1m7dm2d4J7+e7FixThmzBgmdO5MJiXx77//ZvXq1bMU5NOTr68vP/zwQy5evJjxnTuTzs5kpn4MzHL/Ptmkien5mjZ9HpxLlXo+lr0E99eGBHdhEXKw5r6MZXz8+HE2atRIJ4gYaiRnzOPHj/nVV1/R3t5eJ7+6devy8OHDOsucPHlSU5tukjFYvP9+lr9XfHw827RpY1aQTK+tv/vuu1rL//LLLyxbtqzeZWxsbNi5c2eDA+VERESwT58+dHR01Hn04QgwHuDNxYtJGm7sCKhfuTO1/TYAb6eXVYZHAFlizhC9M2Y8/58MGPB8urHgfvTo84sAQ8H93DlSzxsbWSXni9wnwV1YhBysuS9zGatUKi5dulTv++yenp5cMHeu0a5kFy1axO+//57h4eF8/Pgxz507x+bNm2sFotEAz9nYMF7PM2KVSsXNmzezc8mSTEkPBJ99lq3v9vTpU37wwQeadXt4eHD9+vVatfXMAXvAgAFat7xTUlK4fv16vRc96alBgwZct26d3lf67t69y5UrV9LX11czfyWAB1xdyXXrtOa9c+cOP/74Y61n+F27duXt27e5du1aDhgwgJUrV9atwQOcAfBmUBC5dm22ysosFy8+D85//vl8urHg3q8fWbcuuX27bnA/d079vL9MmecXADkg54vcJ8FdWIQcrLnPUBk/fPiQX3zxhVYjubEA7wOM+/RTg/klJSVpasw2Njb09/dn3759+b///Y/FihWjZ8Yg8MYbBvN5+vQpr27dSi5aROagX3qVSqXpxMbNzU0zLSwsjIGBgXqDtZ2dHXv37s1z585p5fXvv/+yW7duet/xB9QN9WbOnMn4+Hit5fbs2cOUlBRu2bKFrVu3pqIoXLZsmcFt/vfffzUXID4+Pjqf37x5k6tWrWK/fv007QscHByYmJiY7XIyW+XKZP78ZMZhYI0F961btQN6xvnS35sfONAimybni9wnwV1YhBysuc9UGWdsJPcg/aTs7Gx0mSdPnrBVq1Z6A2BBW1ueTTvJXwwMZHR0dI5HXTPH1KlTWbBgQZ3pBw4cYNu2bfVuq6IofPfdd3UeIVy6dIlDhgwxeMvczc2N0z75hPfSGqhlLuOoqCiTY6SnpKRwzpw59PDw4LVr14zOe/36dYaFhZlRChYwbBjZqZP2NGPBPTGRzDhmvL5kzmt+ZpDzRe6T4C4sQg7W3GdOGatUKq5Zs4bb3d3Vh2O3biaXefLkidbteK1n+ABrZPi7cOHCbNOmDb/55htu374910aQW7NmjcHPjh8/zu7duxvsxa5ly5bcu3ev1jIPHjzgDz/8wBIlSmjNOzstaD2rXp1MSsrRfnz37l1ev34928tb3L596q5nMzLVoK59e8OBPfNdgByQ80Xuk+AuLEIO1tyXlTJWqVRk5hHBjEhMTGTTpk0NPqvOnLy8vPhnxme5L8CFCxf48ccf620ICICNGjXitm3btO42PH36lMuWLWONGjUIgGczBq/7961rP05J0X1lzlRwnzvXcHDP0Igxp6yqnF9Slg7uMiqcEC8BRVEAb2+z53dycsKmTZsQFBRkct4uXbrg1KlTaNmyZU42McfKlSuHuXPnIjo6GoMHD4aTk5PW5/v370ebNm1Qo0YNrFmzBqmpqbC3t0fPnj1x+PBh7NmzB6HVqyPF2xsYORIoUOAFfZNcYmsLuLtnbZm33zb8Wbt2Odse8UqT4C7EK8rZ2RmbN29G48aNjc53/fp1XLhwIY+2yrQSJUpg+vTpuHLlCr7++msUKlRI6/Pjx49jZteuWFe0KFLThrBVFAWBgYEYfPQo7K5eBSZMUAfD152XF1C7tu50RQGyMMSusD4S3IV4hbm4uGDLli0ICAgwOM+BAwfQoEEDdO7c+aUK8p6enhg/fjyuXLmCkJAQeHp6AgA8APwFoGtcHGw7d36h2/hKaNtWd1q9ekCRInm/LeKlIcFdiJfUvXv3kJqaanK+/PnzY+vWrWjYds6tBQAAIABJREFUsKHR+datW4cqVapgwIABuHXrlqU2M8cKFSqE4cOH4/Lly5g5cya8vL2RlP5hplv31kalUiE+Pj5nmei7/S635F97EtyFeEk9evQI1atXR2BgIIYOHYrVq1fj0qVL6pawmRQoUADbtm1D/fr1NdMmTZqEMWPGwNnZWTMtJSUFs2bNQtmyZfHdd98hISEhT76LOZydnfHZZ5/hSHQ0Tq1aBUydCvz554veLItITU1FZGQkNm3ahEmTJuGDDz5A7dq14eHhgcjIyJxlXq0a4OCgPU2C+2vP7kVvgBBCv9KlS2PTpk0IDAzE3r17NdMLFy6MWrVqoXbt2ppUrFgxFCxYENu3b0eLFi0QERGBIkWKoG/fvujXrx++/fZbzJs3T3MnID4+HmPGjMGsWbMwbtw49O3bF3Z2z08Hx48fR4ECBVC6dGnY5vGzbQcHB9Tu2jVP12kpT58+RWRkJM6dO4ezZ8/i7NmzOHfuHP777z8kJyfrzP/NN9+gevXqOVupogCursDt28+nVamSszzFK0+CuxAvsdKlSyM8PByBgYG4cuUKAODu3bv4888/8WeGWm2JEiU0gX7EiBEYN26cpjW6l5cXfvvtN3z++ecYMWIENm7cqFnu5s2bCA4OxvTp0zFp0iS0a9cOiqLA3d0dLVq0wJUrV1C+fHlUqlRJK5UvXx6Ojo55WxgvoejoaMyfP18TzCMjI816lAIA1apVw8iRIy2zIW5uz4O7ra064IvXmgR3IV5y6QG+cePGiImJ0TvP1atXcfXqVaxfv14z7YcffkCTJk1QJK1hVaVKlbBhwwbs27cPQ4cOxb///quZ99y5c3jnnXcQEBCAH374AXXr1tW8mnbo0CGcPHlSa302Njbw8/PTCvgVK1ZEpUqVdFq/WzM/Pz906tQJ06ZNQ1RUlNmBHVDfHfHx8YGvry/8/Px0fnp7e8PoPZNWrYBevYD33wcKFnw+Pf1Oy8mTwCefAFu3AoULZ+v7iVeXou/53YtUoUIF/vfffy96M6xeem1Q5B5Ll/GlS5cQGBhoMMBn1L59e8ydO1cT2DMjiXXr1mHEiBGIiorS+bxLly6YOHEiihUrhk6dOmHHjh1mb6eXlxfq16+PefPmwc3NzezlsiM7ZZyamorQ0FAkJCTg6dOnePbsmd6fxj5zcnLCtGnTkD9/fk2+sbGxmDlzJubMmYOHDx/m+LvZ29vjM3d3TMvY+DHj+bp1a3WbhKJF1TX1mzfV0+3sgDfeAI4eBQoVAh48yPG2yPki9ymKcoRkLYtlaMkecSyRpIe6vCE9TuW+3CjjixcvsmTJkkZ7owsMDDTZv3q6p0+f8ueff6anp6feQV4+++wzxsbGslu3bmb3hvfmm28yKirK4t9dn+yW8ZkzZ9iwYUOzv1PG5ODgwN27dxvM+9GjR5w+fTpLlSqVrfwzpv6Ze53LaONGw73Tpad27bJVPpnJ+SL3QXqoE+L15efnh/DwcJQsWdLgPOHh4ahUqRKWL18OlUplND97e3sMGDAAUVFRGDVqlFavcSkpKfj5559RuXJlVKlSBcHBwSa3r27duvj7779RpkwZ87/UC1C5cmX89ddfmD17dpYeIyiKgiVLlhjtGbBAgQIYPHgwoqKisHr1atStW1fvfN999x0mTJiAjz/+GE2aNIGvr69O40Wjt+XbtlXX0o0ZMMD458J6WfJKwRJJau55Q67ELWfZsmX84osvuH79et6+fVszPTfLOCoqSmdQFX2pevXqWRrV7Nq1a/z444/1DvLi5eXF9u3bm1xn3bp1uWnTpjwZkc4SZXzjxg2z70zY29uzXbt2/PnnnxkZGWlW/iqVivv372enTp20xpMfN26czrxPnz5ldHQ0d+7cyblz53Jl5oFhMnvzTcO1dnt7i4zlTsr5Ii9ABo4RliAHq+WkpqayT58+mpN2uXLl2LdvX65Zs4bnzp3LtSCXOcD7+PjQ399fb1Bq1qwZjxw5YnbeZ86cYbt27fTmVaxYMa0gZSi9+eabXL16NVNSUnLl+5OW3Y+3bduW5Vvpfn5+/PTTT7l+/Xo+fPjQ5DqioqL42Wef0dnZmZUqVTK9b5gaOGbSJMPBvWbNbJaELjlf5D4J7sIi5GC1rJSUFHbp0kXrxD9lyhQC6mFX27dvz++//54HDhww+3m4OSIjI+nt7a0J4CkpKVywYIHBWn2PHj148eJFs/MPDw9n7dq1nz9vBrgY4EKA9ra2muk2NjZ0d3fXu86KFSty8eLFfPbsmcW+dzpL78fx8fEcOnQobTN8N3OTnZ0dAwICOH78eEZERBgdXjcuLo4hISGMjY3VtxHPfzcU3FUqMiGBfPTIcHCfPt1i5SLni9wnwV1YhByslpecnMy3335bJ7hnTo6OjmzUqBG/+uorbt68mXGZh/nMovQA37ZtW820xMREfv/993R1ddV7a3nQoEFajxCMUalUXLVqFf38/OiZIXiE+fnRxcWFAOjk5MRHjx7x+++/Z5EiRfR+b19fX86ePduiFze5tR8fO3ZM66IGUDekGz16NOvWrWvWnQsPDw/27NqVu6dOzdrKJ04k+/QhIyN1g7tKRW7apK6VR0er5y9RQjewKwqZmGix8pDzRe6T4C4sQg7W3JGYmMigoCCjwT1zwJg0aVKOb91HRkZy8ODBOtPv3bvHoUOH0tHRUWfdBQoU4Hfffcf4jDVFI5KTk/nL1Kn8x92dKnt78ttv+e+//9LDw4MODg5aZTBz5kyDdw+8vb05Y8YMJiQk5Og7k7m7H6ekpPCnn35i/vz5NdseExNDUl3zXrVqFT/88EPNnZPMyQXgmfRg+8kn5q/42DH1Mra2us/Ua9RQ/6xS5fn8wcG6wd3Hx6JlIeeL3CfBXViEHKy55/Hjx6xXr57J4F6rVi2ePn3aYut9+vSpwc9iYmLYt29fvTXOYsWK8bfffjO6fEYqlYq8c0fz99mzZ+nr66szX3JyMufOncsyZcro/f6enp6cNGmSWc+qDcmL/Tg2NpYdOnQgAG7YsEHnc5VKxTNnznDatGls2bIl8+XLRwAskjHYVqv2/+2dd1gVRxfG36ELCDZUUAGxgL2hYjR2RROJxBJNjL1GjRrBYCzRqJjYSOwl+mkSKxYUUYwFe2IvRI1GUYiixgIWlH7P9weXGy63Aru0nN/zzCN3dnbm7HF33t2pxheoUBBVqaIp2FnD5Mn/pv/zT83jY8ZIcOX/wvWF/LC4M5LAD6u8xMXF0YoVK/SK++eff07x8fGSl61QKGj79u0UFhZGMTExaq0Cf/zxB3Xr1k2rPTVr1qSdO3fmqhVBa9+xktTUVNq4cSPVqlVLa7mlS5emGTNm5Kp7Ij/v4127dtGaNWsMpktMTKSDBw/SpC++oJUVKlBazZpEu3blrLDPPtMv7qdPq6cvWVL9eGaTvURwfSE/LO6MJPDDKj+HDh0iDw8PvQJfrlw5WrFiheSDzR4+fEh16tQhAGRvb08tW7akUaNG0fLly+nEiRO0d+9eat68uVabvLy86MSJE5LaQ5Qxq2Dnzp1Ur149reXa2tpSQEAAPX782Og88/s+zo/pfUREtG+fbmEvV44o+wyELl3+PW5vL7k5XF/ID4s7Iwn8sMrP0aNH6cGDB1S1alWDfe9169alQ4cOSVr+s2fPyNPTU2eZlStXpkaNGlHp0qW1Hu/WrRv98ccfktpERBQfH0+NGzcmMzMzreVaWFjQyJEj9bYGZFJs7+PERCJra+3iPmCAZvqsq9VlGVgpFcXWz4UIqcWdV6hjGBmpVKkSjhw5AicnJ1Xcd999h/bt26ulu3btGjp16oTu3bvnfX9vJWXLlsWRI0fQunVrrccfPHiAy5cvIz4+XuvxsLAwNGjQAEOGDMH9+/clsQkASpUqhbNnz2LYsGFaj6ekpGD16tWoUqUKypYti04dOmDcuHEICgrCzp07cfHiRcTFxUlmT6HEygro1En7MW17tWddrY5XpWMAw1/uAKwAnANwFcB1AN8o4zcAuAfgijI01HH+QAC3lWGgofL4yz1/4Ddx+cnq4xs3blC5cuUIAJ08eZIUCgXt3r1b62Azc3NzmjhxomT98W/evKGuXbsabD3IDNrmeFtZWdHMzz+nuLg4SWwiymji/v77742aVpY1VKhQgX7//XciKub38dq1ml/t5uZEugYgNmwo6ap0WSnWfi4koAC+3JMBtCeiBgAaAugihPBSHptERA2V4Ur2E4UQZQDMANAcQDMAM4QQ8m4TxTCFkFq1auHQoUMoVaoU7OzsIIRA9+7dcf36dSxYsAB2WbbsTE1NRVBQEGrUqIHVq1cjLS0tT2VbW1tj9+7d6N27t8G0o0aNQmxsLMaNGwdzc3NV/JikJAQsXQqTzp2BV6/yZE8mQghMmDABe/bsgY2NjVHnNGzYEOfOnYOXl5fhxEWd997TjGvbVn1716z07ZuxGxzv5c4AhsVd+VKRoPxprgxkZP7eAA4RURwRxQM4BKBLrixlmCJOw4YNER4eDgcHB1WcpaUl/P39cfv2bYwcORImJv8+ks+ePcOoUaPQuHFjHDlyJE9lW1hYYMuWLRg8eLDedGvWrMGkSZMwZswY/Pnnn/j4448BAL7IaMKzv3ABePo0T7Zkx8fHBydPnkSlSpUMph08eLDObWyLHY6OQNOm6nHamuQzGT0amDBBXpuYIoNR+7kLIUwBXARQHcByIgoQQmwA0AIZX/ZHAEwmouRs5/kDsCKiOcrf0wEkEtHCbOlGABgBAA4ODk2Cg4Pzel2MARISEtT2omakJzc+TkxMxP379/H69WuNY6VKlULlypVhaWmZJ7vu37+PJ0+eGExXtmxZODo6Ij09HQnR0SifnJyxP7irK2Ai/XCd1NRU3LlzB2/fvtWbzszMDA4ODihfvjySkpKK93386BHw8OG/v+vVAyws8t0Mri/kp127dgW3nzuAUgCOAqgLwBGAAGAJ4CcAX2tJ7w9gWpbf0wH46yuD+9zzB+5Dk5/c+lihUFBISIjO/nh/f3968eJFru1SKBQ0ffp0jbnm2vraTU1NaeDAgRk7oCUm5rpMY3n9+rVRO89BOap+/fr1dPnyZdntKjAyV6sD1Fely2e4vpAfFORoeSJ6oRT3LkT0SGlTMoD1yOhTz04sgKwbT1dWxjEMowMhBHx9fXH9+nXMnz8fJUuWVB1LTU3FwoULVf3x6enpucp/1qxZmD9/viquQ4cO+OuvvzBkyBC1PcXT09Px008/wcPDA4M/+wxRUVF5uzgD2NraYteuXZg4caJavKenJ2rVqqUWl5KSgufPn6NRo0Zo164dQkNDc+WPQk2DBkDlyhl/62uSZ5hsGBR3IYSDEKKU8u8SADoBuCmEcFTGCWR0yV3TcvqvADoLIUorB9J1VsYxDGMAS0tLTJo0Cbdv38aIESMgsgyUevr0qao/PiIiAgCwZcsW9OrVC0FBQTh79ixSUlL05j9p0iSsWrVKla+bmxvWrVuHv/76C4MHD9YQ+Q0bNsDd3R1DhgyRVeRNTU2xaNEirFq1SmVDo0aNcP36dRw4cADe3t4a5xw7dgzdu3eHu7s7li5dqrVbo0giRMY0N4DFnckRBvvchRD1kdHsboqMl4FgIpolhIgA4ICMpvkrAEYRUYIQwlP59zDl+UMATFFmF0hE6/WV5+7uTrdu3crLNTFGcOzYMbRt27agzSjWSO3jq1evYsKECTh27JjGMV9fXyxYsACbNm3CzJkzAQBWVlZo2rQpWrZsiZYtW6JFixYoW7asxrlbtmxBaGgotmzZohYfFRWFOXPm4JdfftH4IjY1NcXAgQMxdepUuLm5SXaN2Tl48CB69+6N9957T82+GzduYPHixXB3d4efn5/Gefb29hg2bBjGjh0LV1dX2ezLCwcPHsTvv/+OtLQ0pKenq0L237WjozH4+HGM7dkTaURqxwIDA1G7dm3ZbeX6Qn6EEAXX554fgfvc8wfuQ5MfOXysUCho165d5ObmprU/3s/Pj8aNG6ezn9rDw4OGDBlC69ato5s3b6qWU33w4IHOMm/fvk0DBw4kExMTjfzMzMxo6NChdFfitcyzcu3aNfr888+1Hjt8+DAFBgaSo6Oj1us1MTGhXr160enTp/Nv6VgjSUtLo6CgILK2ttY7tsASoGVa4hcuXJhvtnJ9IT/g5WcZKeCHVX7k9HFSUhLNnz+fSpYsqbnIS7ly1FvHuvHZQ9myZcnHx4e+++47OnHiBKWnp2st7969e3T+/HkaMGCATpEfNmwY3bt3T5br1bX2fqaPk5OTaePGjdSkSROd19q0aVPavHmz0bvf5Rd3796lzp076/1/Kpftd0BAQL7ayPWF/LC4M5LAD6v85IePHz9+TMOHD1et8lYBoEiAUoWgrXqELnsoX7487dKzc1l0dDR5eHhQ+fLlqWnTplS9enWtK8uZmZnR8OHDZRP57GT3sUKhoJMnT1KPHj20voQAGfvJf/vtt7nahU4uFAoF/fzzz1SmTBmj/r8aNWpEU6dOpVOnTkm+6ZA2uL6QHxZ3RhL4YZWf/PTx5cuXqU2bNtQwc9oUQIquXemjjz4yKBR9+vShp1n2Z9fF06dPde4klz2YmpqSr68vRUZGynrd+nx89+5d+uKLL7S2bgAgJysrWuftTYp8EEdj+eeff+jjjz82+qUMyJjG2KdPH9qwYUOOdtTLCVxfyA+LOyMJ/LDKT0FsRxq6aRNdqlePyNmZaM8eSklJIR8fH73i0LFjR6NFOCEhIUfr1AOgEiVKUJs2bfT26+cWY3z88uVL+uGHH9R257MD6Fnmi9D770tuV14JCwujKlWq5MjPmaFJkyY0bdo0On36NKVl3xo2l3B9IT8s7owk8MMqPwXq4yx954mJidSxY0e9gmBiYkIjR46kf/75x2DWKSkp1L9/f6PFxgqg6SYmdHLsWMkvMyc+TktLo5CQEGrdujWVz9LCQQ0bSm6XFLx69YrGjh2r0f0xdepU8vX1JVtbW4O+L1OmDPXt25d++ukno/5vdcH1hfywuDOSwA+r/BQmHyckJFCrVq3UKn5tfeZ2dnY0b948SkpK0ptfeno6TZo0yShxX5xVSPfulfS6cuvjixcu0I4PPyTy8SE6f15Sm6Tm9OnTVKtWLZU/d+7cSUQZgwiPHDlC/v7+VKdOHf3dJACNBeiqpyfRw4c5tqEw3cvFFanFnfdzZ5j/ADY2NggLC4On57/TaNeuXYuJEyeq7f726tUrBAQEoHbt2ti5c2fGF4AWTExMMH/+fCxcuFDr8aykZv2RpayCpHGTJui5axcQGgp4Sje1WA7eeecdXL58GTNnzoS5uTmuXMnYgNPCwgLt27fHggULcO3aNcTExGDVqlXo3r27xi577QAsBVD/wgXgq6/y/yKYfIfFnWH+I9jb2+PAgQOoW7cuAMDZ2RmLFi3C9evX0b17d7W0d+/eRa9evdC2bVtcunRJZ55+fn745ZdfYGZmpjNNAIAxAAZUrIg1MTFITk7WmZbRjqWlJWbMmIHLly/D2tpaaxpnZ2eMHDkSu3fvxvPnz3H48GH4+fmhdu3aeIksL1nlyuWX2UxBImUzgBSBm+XzB25mk5/C6uPHjx9TzZo16ezZs2rxR44cofr162s06QohaNCgQRQbG6szz/DwcI3FWHRN63JycqJFixbR69ev83wthdXHhY3o6GiKWLGCaMeOXJ3PfpYfcLM8wzB5oUKFCjhy5AiqVKmiFt++fXtcunQJP/74o9qe6USEDRs2oGbNmpgzZw4SExM18uzSpQsiIiLUlre9cOECVq5cqbH868OHD+Hn5wcXFxfMmjULcXFx0l4go4GLiwvaffYZ0LNnQZvC5BMs7gzzH6Ry5cpwdHTUiDc1NcWwYcNw+/ZtTJ48WW3v+Ddv3mD69Olwd3fHli1bNPrjmzdvjtOnT8PZ2RlAxr7ro0aNwu3bt/HLL79orIEeFxeHGTNmwMXFBV9++SUePXokw5UyzH8TgxvH5De8cUz+wBtByE9x8PG9e/cQEBCA7du3axzz8vLC999/Dy8vL7X42NhYdOnSBXv27FHbVEahUCA0NBSBgYG4cOGCRn6WlpYYMmQIJk2ahKpVqxpln5w+TklJQVBQEB48eAArK6tch4oVK+odk1AUKA73cmGHN45hJIH70KRFoVBQaGgoRUZGqhYOKU4+PnHihM512z/55BP6+++/1dLHx8fTixcvtOalUCjo0KFD1K5dO+3TtkxNqX///nT9+nWt58fFxdGJEyfor7/+oiNHjsi6Iczr169p6NChuVpMBgANHjxY53r9RYnidC8XVsB97gxT+BBCwM3NDe3bt0fp0qXRsWNHPHz4EPv27cOzZ88K2rw88+677+LcuXNYv369RnP+5s2b4e7ujq+//hoJCQkAgFKlSsHe3l5rXkIIdOzYEREREfjtt9/gk22f8vT0dPzyyy+oU6cOevToofGVX6pUKURGRqJ27dqq0eOurq5o3rw5PvjgA4wYMQLTp0/HsmXLsH37dpw4cQK3bt3CixcvdE7t04WtrS3Wrl2LXbt2ad0uVx+DBw/G2rVrYWLC1SyT/3Cz/H8UbmaTh6tXr6Jdu3aIj4/HwoUL4e/vDwCoXr06WrRoAS8vL3h5eaF+/fpFtqn2hx9+wPTp0/H27VsoFAq1Y05OTpg7dy769++fI1GLjIzEd999h23btmnkCQCdOnXClClT0KZNGwghAAAnT57EpUuXMGHCBKPLsbS0RIUKFTB69GgEBAQYfR6QMRBw8ODBOHjwoMG0Li4uOHz4MKpXr671OBHh9OnTSElJgRBCLQAwOs7d3V3nS5SUcH0hP9wsz0gCN7PJx/nz58nOzo4WLlyos7nW2tqaWrduTQEBARQSEkKPHj0qaLNzxLZt28jMzEzn9ZUtW5Y+++wz2rt3L925c8foNc5v375Nw4cPJ3Nzc635tmjRgvbu3atqij948CA1a9YsR03lY8eOzfWa6+np6bR48WKytLQ0qqwmTZrQvHnztO53HxkZSY0bN851k7+vr2++7AhHxPVFfgBefpaRAn5Y5eW3336joKCgHFXWixYtKmizc8TBgwfJxsbGqGuzsrKiBg0aUN++fSkoKMiguN6/f5+++OILjbnzmaF+/fq088cf6ejRo5SYmEhDhgwxyo4xY8ZI0kd/7do1rWsC6Auenp40f/58te1wU1JSKDAwkCwsLHKUV6dOnQwuESwlXF/ID4s7Iwn8sMrP3r17dYpT1lClShU6fPhwQZubK86ePUtly5ZVu55ayFjHfCxAn2W71jp16tCff/5pdP5Pnz6l6dOnU6lSpVR5mAG0G6AUgI5u3EhEGYP0Vq5cqfOLP2vo2LEjBQcHU3Jycp6uPSkpifz9/TXyd3Bw0LmXfGZo2rQpLViwQCX0f/zxB3l6ehol7BYWFhQSEpJn+3MC1xfyI7W4c5/7fxTuQ5OfY8eOIT09Hd26dUNSUpLOdOvWrcPgwYNV/ahFjT///BOdO3fGgwcPAACDAfxPeSwRQOZiqR999BH+97//aax7bgyvXr3CqlWrEBQUhDf//IPXyvhjy5ah7ZgxqnSnT59Gr1698PjxY4N5VqhQAYMHD8bw4cPVpuxlhYiwdu1anD17FpaWlmpT3DJ/R0dHY/369Xjx4gUAwNvbGyNHjsSZM2dw9OhRXLhwAfrq2WbNmuGjjz6Cr68vtm/fjhkzZiAlJcWg/SVLlkSXLl3wwQcfoGvXrjke8JcTuL6QH+5zZySB38TlJ9PH4eHhBptdW7VqpbEcbFHi77//Jg8PDwJAg4GMRkGA3ma5xgoVKtD8+fPp1atXuS7n7du3tGLJEgqrUIEUNjZ0dOtWjTSxsbHUokULNf8KIfSOEejcuTPt2LGDUlJSNPJLS0ujefPmGdUqoC20a9eOlixZQu3btzf4Rd+8eXPy9/enhg0b5qgMExMTat26NS1cuJBu3bqVa//qgusL+QE3yzNSwA+r/GT1cWhoqF5xyQx9+/ZV65OVE0n7bH/7jZ7+/Tc1a9ZMp7i7A+QEUOnSpenrr7+m+Pj4vJWZmKjzPk5KSqKRI0eqyjY1NaXY2FgKDAwkV1dXnf6vWLEiTZkyRev/waVLl9S2XjUmeHt7U2JioiqPx48f08qVK6ldu3YGhd7Z2VnjnilRogQ5ODgYLNfd3Z0mTZpEJ0+ezPXgwaxwfSE/LO6MJPDDKj/Zfbxjxw4yNTVVVcBDhw6l8uXLa+1TnTRpUt7FzwA3b96k5s2b04cffkhz5syh8PBwevLkSe4yW7aMqGJFSpo3j5Z4eKiJuwdAGwF6DVCJbF+bVatWpSFDhtD3339PYWFhdOvWLa1fz7owdB+vWbNG1Wry/PlzIsoY8R4eHk4ffvih2v9H9i/9Ll260K5du9RGpL99+5bGjh1rlLDb2NjQkSNHdC5i8/jxY1qxYgW1bdtWr9DbZfnb1taWUlNT6fTp0xQQEGDUy0bZsmVpwIABtGPHjly3mnB9IT8s7owk8MMqP9p8vHnzZlVF/vvvv9OrV69o2rRpVKJECa2V8pIlS3Ikdjnl1q1b5OTkpPHF+OGHH1JgYCAdOHCAnj59ajij6GiVoCtMTVV/E0Dpyn9PZhkUZ6iJuVq1auTt7U1jx46lxYsX0/79++mvv/7S8IUx9/Hvv/9OTk5OdPPmTY1jsbGxNHv2bHJ2dtZpj6OjI02bNo2io6NV5+3fv58qVqxo1PVUqlSJJkyYQKdPn9Yp9I8ePaLly5dT27ZtSQjx7/8FQGkALQDISnnfZLWDKGP6YFBQELVt21bny0rWF0dvb29avny5xqqC+uD6Qn5Y3BlJ4IdVfnT5eMOGDSSEoMjISFXc/fv3acCAAWoVe2aoUaMGhYSEyLbM6l9//UWVKlXSKwouLi7Uo0cPmjt3Lv3666/07NkzzYzq1VMT9ew7lw6PAAAgAElEQVRBsWoVHT58mNq2bWuUKGoLpqam5ObmRkOHDqVEPc3y2Xn06BE9ePBA5/G0tDTav38/de/eXedXtBCC3nvvPdq9ezelpqbSkydPyNfXN0f2V6pUicaPH2+U0Ldp04acAEpW+u9my5bUokULCgsL03kdz58/p02bNlGfPn3Izs5Opx0tAPoNoFdjxxK9fWvQf1xfyI/U4s6j5f+j8OhX+dHn4x9//BEdO3bU2CDl0qVL8Pf3x9GjRzXOad26NRYtWgRPT+kG1GZy584dtGvXTjXi3RiaNm2K3bt3w8nJKSNi6lRg7lzdJ8TGAsq0p06dQmBgIA4cOJAre8uWLYuqVati1KhRuHbtGpydneHi4gJnZ2d4eHjA1tY2V/lmmBmLdevWYe3atbh//77WNJUqVcLQoUMxdOhQHDp0COPHj8ebN29Ux8uUKYO3b9/qnSVhbW0NNzc31KpVC87OzrC0tISFhYVaaNasGZ6cOgUfU1Ng5EikW1jg8ePHqFSpksHrSElJwcmTJxEaGorQ0FBER0erjoUCUC36e+EC0KSJ3ry4vpAfHi3PSAK/icuPIR/r+hJXKBS0d+9e1ejz7KFfv34aTbNScOfOHapSpYpRX6De3t6aq+r9/rvuL/cmTbSWef78eerevXuuv+SzrwLo6+tLr1+/lsQfaWlpFBYWRj4+Pjq/5k1MTKj7e+/RxoULqXnz5qr4Zs2aUWJiIs2fP59KliyZq2sbOnSoZJvOKBQKioyMpMDAQGrevDkFZv6/VKxIFBNj8HyuL+QH3CzPSAE/rPKTVx+npqbSihUrtI6OtrS0pMmTJ+vceS23REVF6e1/BjIWYElISNA8OT2dyMFBu7jPmKG33MjISOrZs2eexD0gIEC2Hdj+/vtvmjFjhkb3RRWAbgKUYmJCqd9/TzNnziRTU1MqXbq06tzXr1/TsGHDcnRdtWrVohMnTsh2PY8ePaJn+/cTGTnAjusL+WFxZySBH1b5kcrHL1++pK+++oqsrKw0RMDBwYGWL18u6aC7u3fvkouLi17xcXJyotWrV2uubT5okHZxv3DBqLL/+OMPql69eo7FvUqVKhQeHq4hhk+fPqX+/ftT165dadCgQRQQEEBBQUG0adMmOnz4MP3xxx/0zz//GC2iqampFBoaSu+//z4JIahl1mvs1YuIMgbwVatWTW1cQmpqKi1YsIDKlCmTI5F3cnKi8ePH06lTpwp061iuL+SHxZ2RBH5Y5UdqH8fExNCnn36qVQQ8PDwoNDRUskF39+7d0zsfPDPUrFmTgoOD/y13505NYXdyIsqBXQqFgsaPH5/jL3cgY3738uXL1ZrmX716RT169NCbj4mJCVWoUIHq169PnTp1ok8//ZT8/Pxo/vz5tG/fPq12xsTE0NxJk+hk5cpEVasShYerjr1+/Zpevnyplj4iIkJjZkJOQkczMzo+dGiOfCkVXF/ID4s7Iwn8sMqPXD6+cOECtWnTRqsAtGvXji5evChJOdHR0VS1alW1VoImTZpoLbdJkyZ06NAhotevSWFhoSbuiuHDc1X+2rVrDU7t0rXznr29Pfn5+akWo1EoFDR37lytsxH0hapVq2qdQqeBkYIbHx9P/fr1y5W4OyFjvf6j/ftLsjBNTuD6Qn5Y3BlJ4IdVfuT0sUKhoN27d1PNmjW1CkH//v1zNI9ZFzExMeTm5kYAqEGDBpSenk7BwcFUo0YNreV26NCBznt5qYl7H2tratq0KX366ac0a9Ys2rZtG12+fFl7v3029u3bR9bW1mRlZUUBAQEaG/Ho21Y384v8ww8/pGPHjpFCoaDw8HC1TWj0hcaNG8u2Fe+2bdv0TlUzFMqXL08jRoyggwcPyroOQiZcX8gPizsjCfywyk9++DglJYWWLl2qsTMbkLHN6pQpU/K0ljtRxhz8atWqUbNmzdTKXb16tc5m5l7IGGiWKITaqnTZQ5UqVahDhw40evRo+uGHHyg8PJyioqLUvkzPnTtHFSpUIIVCQc+fP6epU6eq1nlfuHAhuRgpiA0bNqT169fT9evXqV69egbTd+jQgUJDQ2XbM/3BgwfUqVOnXAt8ZihTpgwNHjyYwsLCZNsGlusL+ZFa3Hme+38UnrcqP/np45cvX2Lu3LlYvHgxkpOT1Y59+eWXmDdvXp7yj42NxYwZM7B27Vq1+Ldv32Lp0qX47rvvVLuiZWIKYFDlyrjp4oLTp0/nqLwBAwZg3bp1MDMzAwBERUXB1dUVpqamADKut1evXujSpQsq+vujDIB51arheFSUwbwdHBwwZMgQ3Lx5E3v27DGYvnLlyqo57VWqVDGY/uHDhxg3bhxiYmJgY2OjCra2tmq/bWxsUKJECZw5cwabN29GWlqawbwNYWdnBx8fH/Tq1Qve3t4oUaJEnvMEuL7ID3ieOyMJ/CYuPwXh43v37tHHH39MAMgaoE42NhR/4EDGNLU8oq/5Ny4ujiZPnqx1GV0rKyudc/azB1NTU1q0aJHRAwN37NhBP4wfrxrMdubMGerRo4dRfevm5ubUuHFjo/vhhRBUr149mjx5MkVERNCNGzfo6dOnWkexJyQk0PDhw/P8VZ6XYGNjQ71796Y9K1fmeaAl1xfyA26WZ6SAH1b5KUgfnz17llq3bElvrK2JGjXKt3JjY2Np5MiRZKpFMA1tewtk7IqXk0Vojh49qnVw2c2bN2n48OFGlQlAbfc1a2traty4sdEiamJiQuXKlaNatWrRu+++Szt27FDZsX37dqP7+DNfbrK+bDRp0oS++eYb6tKli8Z4A2PCRiCjmv/kkzz9v3J9IT8s7owk8MMqPwXtY4VCQemffkr09df5XvZfwcH0UQ5ELWsoU6YMzZgxQ/v69dkw5OOHDx9SQEBAjgevPX78mH788UejtunNGvz9/TW+5GNiYujdd9816vy+fftSWFgYVatWTfXikP1FwsrKyuB2sZnhVaa429rm5b+zwO/l/wJSi7sJDCCEsBJCnBNCXBVCXBdCfJPt+BIhRIKOc12FEIlCiCvKsMpQeQzDSIMQAiY+PkC3bvledo2ePbFt7lxcuHABnTt3ztG5cXFx+Oabb+Di4oKJEyciNjY213Y4Ojriu+++w/379zF//nw4OjoadZ6rqyvOnTuHTZs2wdXV1ahz7OzsYGNjo2Gvs7MzIiIi8M0338DERH+Vu3XrVnTv3h0uLi5o06YNrKysUKdOHdVxhUKBpKQkKBQKo2z6HsBLAHtr1MCvv/6K1NRUo85jigGG1B+AAGCr/NscwFkAXsrfngB+AZCg41xXANdy8rbBX+75A7+Jy0+h8HFCgiT97bkiS3P5kSNHqGnTpjq/MO3t7XXOoTc3N6ehQ4fSrVu3NIrIqY+TkpJo3bp15O7ursq/IUCuOuzyAmibhwf16NAhR8303bp1o71792p0GZw8edLg8r5Zg5mZGXXp0oVatmyZq1YQAGSa5e+hNjZ0fuBAohzOACgU93IxB/n95a4sN/PL3FwZSAhhCmABgC8N5cEwTAFhYwMY+FqUDeXIdgBo3749zp49i507d8Ld3V0j6cuXL2FnZ4clS5agQ4cOasdSU1Oxbt06eHh4oHfv3rh06VKuTbK0tMSQIUNw48YNhISEwMvLC80B3AEQAqBtll36qgE4DeCjmzexw8wMM2bMMKoMhUKBsLAw+Pj4wMXFBV9//TWWL1+Ow4cPw9HREefOnUOvXr00ztPml7S0NBw4cACnT5+GpaUlhBAGy3///fcxaNAglCtXDgCQroz/BMDaN2/g+dNPQGCgUdfCFF2MmgqnFPKLAKoDWE5EAUKI8QBMiOh7IUQCEWnssSiEcAVwHcBfAF4BmEZEJ7WkGwFgBAA4ODg0CQ4Ozv0VMUaRkJCQp20xGcOwj3Xz7NkzPHz4UK2Z2BpAJQCW5csjvUwZPHr8WGN6XSZ2dnaqJva8+jjh9WuYx8XB0twcqFgRicnJePLkCRKeP0edzPrRzg6oUQMvX77EvXv3kJ6erj9TLQghMls0YWFhASGE2rRFKysrWFlZqZreU1JScn1NJiYmKF26NGxsbJCYmIgXL16gZGoqVK8uTk6AkV0UAN/L+UG7du0KbiocgFIAjgJoDeAUADNlvK5meUsAZZV/NwFwH4CdvjK4WT5/4GY2+WEf6+ft27e0YMEC1Wjyk5mDvwCiqCgiIrpx4wYNHDhQ58C2lStX0p49e2TZVOXJkyd07NtvMwYkZtl9Lyoqiho0aKCywcXFhY4fP079+/cnS0tLnc3jHgAdASgYIMccNKvndMnc7KFu3boUFBREoaGh9Pv48URLl6p1mRgD38vyg4IeLQ/gawAzADwGEK0MCgB3jDj3GABPfWlY3PMHfljlh31sHPHx8TR58mTaYmqaUSWVK0f0zz9qaaKjo+nzzz/XmEefufxsnTp16Oeff86XpViJiN68eUMDBw4kIGMef+Y88vj4eFq+fDk1atRIQ2TnZHl5GZQHsdYX9L1cWFpa0scff0xHjhzJ8csQ38vyk+/iDsABQCnl3yUAnATQLVsaXV/uDgBMlX+7AYgFUEZfeSzu+QM/rPLDPs4ZT548IcXevUTx8XrTTJs2TfW1n31teVdXV1q2bBm9fftWdnsfPXpEM2bMIDMzM3qR5cs+k4sXL9Lo0aOpZMmSBID6AJQIUDxALWUS96xB33S5atWq0dy5c+nhw4dGXSvfy/IjtbgbM9LGEcBRIUQkgPMADhFRmK7EQogPhBCzlD9bA4gUQlwBsAPAKCKKM6JMhmH+Yzg4OEB06waUKqU3zezZsxETE4P58+fD3Nxc7Xh0dDTGjh0LV1dXfPvtt3j58qVs9pYpUwb3799HWloaqlWrBk9PT/To0QMTJkzA999/j3v37mHQoEGIjIzEzz//jEetW6MWgBrIGKiXV2xtbREUFITOnTtrHWi3RKHAegDv2NhoHIuKisKUKVNQpUoV+Pr6IiwsTJLlb5lChJRvClIE/nLPH/hNXH7Yx/ITERFBq1atUu1clz3Y2dlRQECAbLu7KRQKmj9/vsF+cUtLS6pRowY1bdqUKleuTFZWVpJ8nbu5udHy5cvp1q1b1LBhQ2rVqpVqfMKfWboBrp89S1988YXWDYYyQ6VKlWjatGmqbXKzwvey/KCg+9zlDizu+QM/rPLDPpafTB+npqbS5s2bqX79+lqFa4CZGR0cNEg2O0JCQnK1PKxUoVy5cjRz5ky6d+8ePXv2jFavXk1zPTwoydGRaNo01VoHSUlJtHXrVurYsaPOvIQQ1KlTJwoODqbk5GQ1PzPyweLOSAI/rPLDPpaf7D5WKBS0b98+atWqlUqsRmX5gqVly2Sz5dKlS1SpUiWjxLhVq1ZUvnx5yUXe2tqaxo0bp/XrOztRUVE0depUndv2VgJonZUVPf7lF76X8wGpxb2AVrdgGIaRHiEE3nvvPZw8eRInT57Ee++9B7VeeRmXX23UqBHOnTsHT0/DU5Xr1auH4OBg1K9fX1Ib3r59iyVLlqB69ero168frly5ojOtm5sb5syZg5iYGISGhuKDDz5QbakLAFsBDElKQoX+/YE8zLlnCgYWd4ZhiiWtWrXCvn37MPTiRUT7+wNLlwKffy5rmU5OTjh+/LjWFeiysnLlSrRt2xYVK1Y06mUAAFq2bIlSegYbZiU9PR2bN29Go0aN4O3tjYiIiIymWi2YmZnBx8cHe/bswd9//43AwEBUrVoVbzITmJsDRqyMxxQuWNwZhinWNGjcGK4LFgBjx6otiSsX1tbW2LZtG6ZOnWow7cGDB3HhwgWULl3aYFpLS0usWbMGbdq0QadOnWBnZ2eUPQcPHkSHDh3QtGlTBAcH6x0V7+TkhClTpuDOnTuw2LcPCbNnA9evZwg8U6RgcWcYhpEYExMTzJkzBz///DMsLCxU8e+++y4++eQTmJmZqaWPj483mGdERAQ++ugjPHz4EC1atMCtW7fwww8/oGqW9fD1cfHiRfTp0wfu7u5YsWIFEhMT9drf7r33YDttGlCjhlH5M4ULFneGYRiZ6N+/P44cOaLaxKV8+fLYtGkTYmJiMG3aNFV8Trh9+zZmzZqF2rVr49GjRzh69Ch27tyJli1bGnX+3bt3MWbMGLi4uGD27Nl4/vx5jm1gCj8s7gzDMDLSqlUrnD17FrVr18arV68AZDR/z549G/fv38e6desMDqwrXbo0rKys1OLi4+Mxb948VK9eHTt27MD333+Ps2fPom/fvmoD43Tx9OlTfP3113B2dsb48eMRExOT+4tkCh0s7gzDMDLj5uaG3377De+8845avJWVFYYMGYIrV67g6NGjaN26tdbz4+PjkZSUpPVYWloatmzZgmbNmmHixIno1asX7ty5g0mTJsHe3l6VzhvAzwAaZDs/c4R9bzc37GjfXtYZBUz+weLOMAyTD9jb22PmzJlajwkh0LZtWxw/fhz79++HjZYlY43h9OnT6NWrF9q1awdHR0fcuHEDS5YsQbVq1eAIoD+AywAO9O2L9u3bq85rBOCMQoFeR48Co0blqmymcMHizjAMU4jo2rUrIiMj4erqajCtnZ0d2rdvD2tra7X46OhoTJw4ER4eHrh79y4OHDiAHlu24NoHH0CMGAHvJUtw5MgRnD9/Hh999BFshPhXDN680SiHKXqwuDMMwxQyMpvx69Wrpzfdq1evEBERgVKlSqFr165wcnJSO/769Wv88MMPcHd3x4bt2/Fi0iTQqlWAgwMAwNPTE9u2bcOG27dxZuJEYOJEYMUK2a6LyT9Y3BmmiJDCq4T9p3B0dMTx48c1+um18fDhQ4SHh6NNmzbYunUrmjdvrnZcoVBg165dePfdd9GsWTNs3rwZqVn61qtVqwavRYuARYuAMmUkvxYm/2FxZ5giwl9//YV3330Xo0ePxrZt2/Do0aOCNomRmdKlS+PgwYPo0qWLzjROAASATgBsQkIQdesWDq5bh99++w29e/eGiYl6NX/hwgX069cPVatWxXfffacawc8UL1jcGaaIULduXcyfPx8//fQT+vbtCycnJ9SsWRPDhw/Hxo0b8ffffxe0iYwM2NjYYM+ePejbt68qrm3btggJCUGjRo3QGEAkgHUAfkxKQqMZM7CtcWMcOHAAq1atwt27d+Hn56exol1sbCxmzpypcxQ+U7RhcWeYIkSLFi2wa9cumCuXA719+zbWrl2L/v37w8XFBVWrVsWgQYOwfv16REVF6VxPvKBJTEzEypUrcezYMSQkJBS0OYUeCwsLbNy4EZ999pkqztfXFxcvXsToHTtQXQhUUcZ3BbAzJQWzZs2Cq6srVq9ejcmTJ+PB/ftYvHgx3NzcYAZgCYB91aqh/OzZwNWr/xZWSO8ZJmeYGU7CMExhwtvbGxs3bkTfvn01xDs6OhrR0dH46aefAACVKlVC69at0aZNG7Rp0wbu7u4QhWATkBIlSqB+/frw9vZGYmIi6tatCy8vLzRv3hxeXl7w8PDQaE7+r2Nqaorly5ejXLlyOH78OICMKXRde/YEvf8+EBYGAEgAcEx5zuvXr/Htt99i8eLFCOzRA8MtLDDm6FHsO3MGH/TpA9y4kRHatwdq1gTWrAEsLIAsLxFMEUXK/WOlCLyfe/7A+zPLj9w+Xr16tdH7fAshaPz48ZSQkKA3z/v379OXX35J27ZtowcPHshqPxHRsWPHqESJEhr22tnZUceOHWnatGkUFhZGT58+1Xq+1D5OS0sjhUIhaZ5yoHHdq1er9qx/0qoVtW/fXsOnJgA9ASjVxIQS+vX7d497gKhXL6KKFTP+PnvWcHmM5EDi/dwLXMyzBxb3/IEfVvnJDx/PnTvXoLDXqlWLfvvtN6PzDA0NJTMzMwJArq6u1K9fP1qxYgVdvXqV0tLSJL+Gw4cPk5WVlcHrqFatGn3yySe0ZMkSOnfuHCUnJ0vu49TUVBo7diy1bt2axo0bR//73//o0qVLlJSUJGk5kvPgwb9CvW4dERGdOnWKvL291Xy4IaugawsVKxKlp2tkz/WF/LC4M5LAD6v85IePFQoF+fv76xXFDh060J9//pmjfLdu3UpCCK1f1F26dKHZs2dTRESEwZYAYwkPDycLCwujWyKcnJzo8OHDsvg4LS2Nhg4dqlaemZkZ1a9fnwYMGEBBQUEUERFBz58/l7zsPNG4cUaV/vixWvTZs2epW7duBIB6GhL3oUO1Zs31hfywuDOSwA+r/OSXjxUKBQ0ZMkSvGJqZmZGfnx+9fPnS6HzXrVtnUGTNzMzI09OTJkyYQNu3b6eHDx/m+jr27t1L5ubmBsvs168fxcXFEZF8Pk5PT6fRo0cbtMXZ2Zl8fHxo+vTptHPnToqKiiq4Zv0ZM4iaNdN5+NKlS9TPx4eS9Yj7g+XLtZ7L9YX8sLgzksAPq/zkp49TU1OpR48eBsWoQoUKtGHDBkrX0vSqjcWLFxv9NQ2AGjduTNHR0bm+jpCQEFWXgLZga2urZr+cPlYoFOTn55ej6zc3N6c1a9bIZpNeLlwgmj3bYLJXXl5ahT0RIFshaEDv3nTj4kW1c7i+kB8Wd0YS+GGVn/z2cVJSEnXo0EElNDNnzqQxY8aQiYmJhgh5eXnR+fPnjcrXmH59c3Nzmj17NqWkpOT5OoKDg7XanDU0bNhQtmb5rCgUCpo2bZpRwl6tWjWjfSqTsUQxMYbTLV6sVdz3AdQCoBiAXtjbE50+rTqF6wv5kVrcea4JwxQTLC0tERISgqZNmwLIWFJ02bJluHTpksZWomfOnEGzZs0wfPhwPH36VG++X331Fb766iu9aerWrQtfX1/V/Pu80Lt3b/zyyy96p+xduXIFHTt2xJ07d3Dt2rU8l6kLIQRmz56NwMBAg2n9/f3RqFEj2WwxiBCAs7PhdD4+WqP3m5igBQBnAPYvXwK//y6peUz+wuLOMMWIkiVLIjw8HLVr10bp0qUBAA0aNMCxY8ewZcsWVKpUSZWWiLB27VrUqFEDS5YsQVpams58AwMD8fnnn+s8fvnyZTRo0ACjR4/Gs2fP8nwdn3zyCdavX68m8CtWrIBzNvF6+fIlGjRogBEjRsi6HO+UKVMQFBSkN81nn32GBg0aYNeuXRnNooWVqlWBOnU0ogNOnkSZnj3xtEKFjONZtoRliiBSNgNIEbhZPn/gZjb5KUgfP3jwgG7fvq0Rn5CQQFOmTNE6Mr1u3boUERGhM8/09HQaPHiwweZpe3t7CgoKouTk5Dxfx5o1a1T5JicnU2JiIs2fP5/s7e0JAC1cuFB13MbGhmbOnCnZCH5trFixwuixB/v37y+8c+YDAtSb5Rs2/PdYampGyALXF/ID7nNnpIAfVvkpzD6+c+cO+fj4aBWm3r17U4yOvtu0tDTq06ePKm2LFi1owoQJWgfB1ahRg0JDQ/MscMuXLyczMzO1uKdPn9K4ceNo0aJFGuU6OjrS2rVrZZmTT5QxiyDrNMGePXuSs7OzVl++8847el+YCoxTp9TFffp0vckL871cXGBxZySBH1b5KQo+3r9/P9WsWVNDlEqUKEGzZs2ixMREjXNSUlJULwa9e/cmIqKbN2+q5lJnDx07dqTIyMg82blOuTBLdg4cOEA9e/bUWm7dunUpPDxclq/nTZs2kampKQGgjRs3UlJSEi1btowqVqyo1Zb27dvnaCEh2UlLIypX7l9xP3dOb/KicC8XdVjcGUngh1V+ioqPk5OTad68eWRra6shSlWrVqWQkBANgUxMTKQOHTrQ0GyLnvz6669Up04djXxMTExo1KhR9OTJE0ltz/TxqVOnqHnz5lqFtVOnTnTlyhVJyyUi2rFjB5mbm9PWrVtVcW/evKEFCxZQ2bJltdry/vvv06VLlyS3JVcMGJAhATpWpctKUbmXizIs7owk8MMqP0XNx7GxsfTpp59qFaXOnTtrrHL3+vVrCg4O1sgnNTWVVqxYoVXg7O3tadGiRZL0xxOp+1ihUFBwcDC5ublplCuEoEGDBtH9+/clKTeTvXv30t69ezXiX716RbNnz1aNDcgeevbsSdeuXZPUlhwTHJwhATpWpctKUbuXiyIs7owk8MMqP0XVx6dOnaJGjRppCFJOV7mLi4ujiRMnau2Pr169Ou3ZsyfPTebafJyUlERBQUFUunRprd0N06ZNo1evXuWp3Kzou4a4uDiaOnUq2djYaH3h6Nevn9aBj/nCy5dE5uZEu3cbTFpU7+WiBIs7Iwn8sMpPUfZxWloarV69WuvXd05Xubt165bOwXsdOnTIU3+8Ph/HxcWRn5+f1pkB5cuXp5UrV1JqtlHhcvHkyRPy8/PTukGOqakpDRs2TOcgRlnx8SEyYnZBUb6Xiwos7owk8MMqP8XBx8+fP5dklTsiooMHD+rsjx85ciT9888/ObbPGB/fvXuX+vbtq/XlwsPDg/bu3ZtvU9YePHhAo0eP1rqGvoWFBY0dOzZP6/PnmBs3jEpWHO7lwg6LOyMJ/LDKT3Hy8dWrV6l169aaA8QA2tS9u9H5pKam0sqVK6lcuXIaednZ2dHChQtz1B+fEx+fOXOGWrVqpVXk27ZtSxcuXDA6r7xy7949GjJkiGrEfdZgZWVF/v7+OvewLwiK071cWGFxZySBH1b5KW4+VigUtHXrVqpcuTIBoH6Z06gAokWLcpRXfHy83v743bt3G/U1nVMfKxQKCgkJoRo1aqiVaQ6QP0D7e/c2OHJcSm7dukUff/yx1u11bW1tadq0aRQfH59v9uiiuN3LhRGpxd3g8rNCCCshxDkhxFUhxHUhxDfZji8RQiToOf8rIcQdIcQtIYS3ofIYhimcCCHQp08f3Lx5E1OnToWdqem/B9++zVFepUqVwqJFi3D9+nV88MEHasfu3LkDX19fdOzYEZGRkVKYrkIIAV9fX1y/fh1Lly5F2bJlAQBTASwA0HX7dmDjRknL1EfNmjWxefNmXL16FR9++NcvR38AAAhGSURBVKHasYSEBFyYMwcnK1TA2+PH880mpnhgzNryyQDaE1EDAA0BdBFCeAGAEMITQGldJwohagPoC6AOgC4AVgghTHWlZxim8GNjY4M5c+bA7/p1PJoyBZg7F5g0KVd51axZE3v27MGhQ4dQt25dtWMRERFo1KgRRo4ciSdPnkhhugpzc3OMHTsWUVFRmDx5MhRZX1QsLCQtyxjq1auHXbt24fz58+jatSsAwBRAOACflBRYDxyY7zYxRRuD4q5sMcj8MjdXBlKK9AIAX+o5vTuArUSUTET3ANwB0CyPNjMMUwio5u4Ox8BA4KuvAEvLPOXVsWNHXL58GatWrUK5cuVU8c4KBXqsWYO4OnWAGzfyarIG9vb2+PbbbzH0zh28WLQI2LMH6NtX8nKMxdPTE/v378epU6fQpk0bqLbgcXAoMJuYoolRu8IJIUyFEFcAPAFwiIjOAhgLIJSI9G3FVAnA/Sy/HyjjGIZh1DAzM8PIkSNx+/Zt+Pn5wdzcHF0AeAPwePYM2LlTtrKdXV1RauJEIFsXQUHRsmVLHD56FC/OnQP+9z/g5MmCNokpYoiMfnwjEwtRCkAIgBkA5gJoS0RpQogEIrLVkn4ZgDNEtFH5ex2AcCLakS3dCAAjAMDBwaFJcHBwbq+HMZKEhATY2mr8lzESwj7OG8nJyXgWHQ3HN29gIgTg5gbY26ulYR/nD+xn+WnXrt1FIvKUKj+znCQmohdCiKMA2gGoDuCOcr9layHEHSKqnu2UWABVsvyurIzLnu8aAGsAwN3dndq2bZsTs5hccOzYMbCf5YV9LA3pT5/CND0dqFhR4xj7OH9gPxc9jBkt76D8YocQogSATgAuElFFInIlIlcAb7UIOwCEAugrhLAUQlQFUAPAOenMZximuGPq4KBV2BmG0Y0xX+6OAH5SDqAzARBMRGG6EgshPgDgSURfE9F1IUQwgBsA0gCMIaJ0KQxnGIZhGEY7BsWdiCIBNDKQxjbL36HI+GLP/B0IIDAPNjIMwzAMkwOMGi3PMAzDMEzRgcWdYRiGYYoZLO4MwzAMU8xgcWcYhmGYYgaLO8MwDMMUM1jcGYZhGKaYweLOMAzDMMUMFneGYRiGKWawuDMMwzBMMYPFnWEYhmGKGSzuDMMwDFPMYHFnGIZhmGIGizvDMAzDFDNY3BmGYRimmMHizjAMwzDFDBZ3hmEYhilmsLgzDMMwTDGDxZ1hGIZhihks7gzDMAxTzGBxZxiGYZhiBos7wzAMwxQzWNwZhmEYppjB4s4wDMMwxQwWd4ZhGIYpZggiKmgb1BBCvAZwq6Dt+A9QDsCzgjaimMM+lh/2cf7AfpYfdyIqKVVmZlJlJCG3iMizoI0o7gghLrCf5YV9LD/s4/yB/Sw/QogLUubHzfIMwzAMU8xgcWcYhmGYYkZhFPc1BW3AfwT2s/ywj+WHfZw/sJ/lR1IfF7oBdQzDMAzD5I3C+OXOMAzDMEweKFBxF0L0FkJcF0IohBCeWeI7CSEuCiH+UP7bPsuxA0KIq8rzVgkhTAvG+qJBTn0shLAWQuwTQtxUnvddwVlfNMjlfRwohLgvhEgoGKuLHrn0cxNl/B0hxBIhhCgY64sGenxcVghxVAiRIIRYlu2cPkKISOV58/Lf6qJFLn38sfI+jlRqYDlD5RT0l/s1AD0AnMgW/wyADxHVAzAQwC9Zjn1ERA0A1AXgAKB3fhhahMmNjxcSkQeARgBaCiG65oulRZfc+HgvgGb5Y16xITd+XglgOIAaytAlH+wsyujycRKA6QD8s0YKIcoCWACgAxHVAVBRCNEhPwwtwuTUx2YAFgNoR0T1AUQCGGuokAKd505EfwJA9pdpIrqc5ed1ACWEEJZElExEr5TxZgAsAPCgAT3kwsdvARxVpkkRQlwCUDmfzC2S5PI+PqPtHEY3OfUzgDIA7LL4+mcAvgDC88XgIogeH78BcEoIUT3bKW4AbhPRU+XvwwB6Ajgis6lFllz4WCiDjRDiOQA7AHcMlVPQX+7G0BPAJSJKzowQQvwK4AmA1wB2FJRhxQgNHwOAEKIUAB/wgyoFWn3MSE5WP1cC8CDLsQfKOEY67gBwF0K4Kr8wfQFUKWCbihVElArgMwB/AHgIoDaAdYbOk/3LXQhxGEBFLYemEtEeA+fWATAPQOes8UTkLYSwArAJQHsAhyQyt0gih4+VD+oWAEuI6K5UthZV5PAxown7WX7y4uPsEFG8EOIzANsAKAD8BqBa3q0s2kjpYyGEOTLEvRGAuwCWAvgKwBx958ku7kTUMTfnCSEqAwgBMICIorTkmySE2AOgO/7j4i6Tj9cgo7nth7zaVxyQ6z5m1JHYz7FQ71KqrIz7T5NbH+vJby8yxpBACDECQLqU+RdFJPZxQ2WeUQAghAgGMNnQSYWyWV7ZHLwPwGQiOp0l3lYI4aj82wzA+wBuFoyVRRtdPlYemwPAHsCEgrCtuKDPx4x06PIzET0C8EoI4aUcJT8AQI6+mhjDCCHKK/8tDWA0gLUFa1GxIxZAbSGEg/J3JwB/GjyLiAosAPgQGf1gyQD+AfCrMn4agDcArmQJ5QFUAHAeGaMFryGjecKsIK+hsIdc+LgyMgYp/pklflhBX0dhDjn1sfLYfOU5CuW/Mwv6Ogp7yKWfPZV1RRSAZVAu3MUhZz5WHosGEAcgQZmmtjJ+C4AbytC3oK+hsIdc+niUsk6OREYrSVlD5fAKdQzDMAxTzCiUzfIMwzAMw+QeFneGYRiGKWawuDMMwzBMMYPFnWEYhmGKGSzuDMMwDFPMYHFnGIZhmGIGizvDMAzDFDNY3BmGYRimmPF/DJ0tZ6ZIy54AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light",
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make a map with the vectors\n",
    "fig0=plt.figure(figsize=(8,6))\n",
    "for ista,sta in enumerate(stalist):\n",
    "    plt.quiver(gps_loc['longitude'][gps_loc['station']==sta],gps_loc['latitude'][gps_loc['station']==sta],Ve[ista]/100,Vn[ista]/100,scale=1,units='xy')\n",
    "    plt.quiver(gps_loc['longitude'][gps_loc['station']==sta],gps_loc['latitude'][gps_loc['station']==sta],0,Vu[ista]/50,scale=1,color='r',units='xy')\n",
    "plt.grid(True)\n",
    "plt.axis('equal')\n",
    "plt.xlim(-123,-118)\n",
    "plt.ylim(32,39)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "avsDxHQwzmQj"
   },
   "outputs": [],
   "source": [
    "# Next activities could be fitting a sinisoidal functions\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "bG3GpKz_mO0t"
   ],
   "name": "EPS55_PS4.ipynb",
   "provenance": []
  },
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}