{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solid Earth tides - time series and frequency aliasing" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Go to directory /Users/yunjunz/Papers/2021_Geolocation/figs_src/SET\n" ] } ], "source": [ "%matplotlib inline\n", "import os\n", "import datetime as dt\n", "import numpy as np\n", "from scipy import signal\n", "from matplotlib import pyplot as plt, ticker, dates as mdates\n", "import pysolid\n", "from mintpy.utils import utils as ut, isce_utils\n", "plt.rcParams.update({'font.size': 12})\n", "\n", "work_dir = os.path.expanduser('~/Papers/2022_Geolocation/figs_src/SET')\n", "os.chdir(work_dir)\n", "print('Go to directory', work_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fig:SET - Calculate" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PYSOLID: calculate solid Earth tides in east/north/up direction\n", "PYSOLID: lot/lon: 34.0/-118.0 degree\n", "PYSOLID: start UTC: 2019-12-01T00:00:00\n", "PYSOLID: end UTC: 2021-02-01T00:00:00\n", "PYSOLID: time step: 1800 seconds\n", "Done.\n" ] } ], "source": [ "## SET time-series\n", "# inputs\n", "lat, lon = 34.0, -118.0 # Los Angles, CA\n", "step_sec = 60 * 30\n", "dt0 = dt.datetime(2019,12,1,0,0,0)\n", "dt1 = dt.datetime(2021,2,1,0,0,0)\n", "\n", "# run\n", "(dt_out,\n", " tide_e,\n", " tide_n,\n", " tide_u) = pysolid.calc_solid_earth_tides_point(lat, lon, dt0, dt1,\n", " step_sec=step_sec,\n", " display=False,\n", " verbose=False)\n", "\n", "# tide: ENU2LOS\n", "inc_angle = 42 * np.pi / 180.\n", "az_angle = 100 * np.pi / 180.\n", "tide_los = ut.enu2los(tide_e, tide_n, tide_u, inc_angle=inc_angle, az_angle=az_angle)\n", "\n", "# SET time-series at SAR acquisitions\n", "dt_sen12 = np.array([dt.datetime(2019, 9, 3, 13, 30, 0) + dt.timedelta(days=12) * i for i in range(-200,200)])\n", "dt_sen06 = np.array([dt.datetime(2019, 9, 3, 13, 30, 0) + dt.timedelta(days=6) * i for i in range(-400,400)])\n", "tide_sen12 = np.array([tide_los[dt_out==d] for d in dt_sen12 if np.sum(dt_out==d) > 0], dtype=np.float32)\n", "tide_sen06 = np.array([tide_los[dt_out==d] for d in dt_sen06 if np.sum(dt_out==d) > 0], dtype=np.float32)\n", "# remove zero values\n", "dt_sen12 = dt_sen12[np.multiply(dt_sen12>=dt0, dt_sen12<=dt1)]\n", "dt_sen06 = dt_sen06[np.multiply(dt_sen06>=dt0, dt_sen06<=dt1)]\n", "\n", "print('Done.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fig:SET - Plot time-series of SET up components" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "save figure to file: /Users/yunjunz/Papers/2021_Geolocation/figs_src/SET/SET_TS.pdf\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAACrCAYAAACkAK8AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACWL0lEQVR4nO1dd3gU1fp+z5bUTSWhJCC9S1NUBEVQFBCxYQUFe7sq9p/lKth7uXbhKqIIil0UBRSUJkhTeq8JmEJI3yRbzu+PzTc5c/bM7mxIQuKd93nywM7OzsyZOfO9Xz+Mcw4LFixYsGDBQuOG7VhfgAULFixYsGAhPCzCtmDBggULFpoALMK2YMGCBQsWmgAswrZgwYIFCxaaACzCtmDBggULFpoALMK2YMGCBQsWmgAswrZgwYIFCxaaACzCtmDBggULFpoAHEZfMMb2mzyGm3PetY6ux4IFCxYsWLCggCFhA0gDMDLM7xmAb+vucuoeaWlpvF27dsf6MixYsGDBggUl1qxZk885Tw+3XyjC/oxz/lu4AzDGPo/oyhoIjLHRAEZ36tQJq1evPtaXY8GCBQsWLCjBGCtkjE0BMIdzPsdwv396L/H+/ftzi7AtWLBgwUJjBWNsDee8f7j9QlnYqoMmAnCJ2zjnByO8NgsWLFiwYMFChDBF2IyxYQCmAGiLQNyawAHY6+G6LFiwYMGCBQsCzJZ1vQ/gGQBJAJzCX1Q9XZcFCxYsWLDwj8Rvv4VND1PCrEs8BsA0zrmvVmexYMHCUcHv9yM/Px+FhYXw+azX0IKFpgS73Y7k5GSkpaXBZrPh7LPPRnFxMWJiYiI6jlnCfhXAA4yx5/g/PUvNgoVGiKysLDDG0K5dOzidTjDGwv/IggULxxycc3g8HuTk5CArKwvHHXccEhMTUVJSEjFhm3WJfwngRgBFjLHd4l+kF2/BgoXIUVZWhszMTERFRVlk3cTBOUdFRYVuW1ZWVtA2C/8MMMYQFRWFzMxMlJWVAQCcTie8Xm/ExzJrYX8BYAmAzwG4Iz6LBQsWjho2m9VJ+J8Azjk2bdqEE088UdtWUlKC5OTkY3dRFuoNfr8fxcXFuufrcDjqlbDbA+jHOfdHfAYLFixYsKDB7/dbytf/ELxeL/bt21cnhG121nwL4MyIj27BggULDYjJkyfjqquuavDfRgLOeZMJazDGsHPnzgY953vvvYe77roLALB3714wxmpFbg2N9evXY+DAgUHbVQqaw+GAx+OJ+BxmCTsawHeMsXmMsY/Ev4jP2EBgjI1mjE0pKio61pdiwcI/HkuXLsXAgQORlJSE1NRUDBo0CKtWrTrWl9UooSJsxhjkfF4zJHXo0CGcf/75yMjIAGMMe/fu1X1/3333oXPnzkhISEC3bt3w0UeNVmQDAKqqqvDUU0/h/vvvr/NjFxYW4rrrrkPLli2RkJCALl264Pnnnw/a75prroHD4cDBg/qeYJMnT4bT6YTL5UJycjIGDhyI33//Xfu+d+/eSE5Oxpw5+s6iRoQtPd8kxtiU6pbahjBL2JsAPA9gOYBd0l+jBOd8Duf8pqSkpGN9KRYsHFN8sy4bg55biPYP/oBBzy3EN+uy6/T4xcXFOO+883DHHXegoKAA2dnZmDRpEqKjo+v0PP9kEIGLQvzPP/8M+zubzYYRI0bgyy+/VH4fHx+POXPmoKioCNOnT8fEiROxfPnyOrnm+sC3336Lbt26ITMzs86Pfffdd6O0tBRbtmxBUVERvvvuO3Ts2FG3T1lZGb788kskJSXhk08+CTrG5ZdfjtLSUuTn52Po0KG49NJLdd+PGzcO7733nm6birAVSWdFnPObQvURB0wSNuf8caM/M7+3YMHCscE367Lx0FcbkF3oBgeQXejGQ19tqFPS3r59OwDgyiuvhN1uR2xsLM455xz07t0bALBr1y6ceeaZaNasGdLS0jBu3DgUFhZqv2/Xrh1efPFF9O7dG/Hx8bj++uuRk5ODkSNHIiEhAcOGDcORI0cA1LhIp0yZgoyMDLRq1Qovv/yy4bWtWLECAwcORHJyMvr06YNff/1V+27Pnj0444wzkJCQgLPPPhv5+fkhx/n999+jb9++mnW1fv16AMBzzz2HSy65RLfvxIkTceeddwIAtm3bhuuuuw6tWrVCZmYmHnvsMfj9gXSgDz/8EIMGDcKzzz6L4447Do8++ihSU1OxYcMG7Vi5ubmIjY1FXl5e0DW1aNECt912G0466STlNT/++OPo1q0bbDYbTjnlFJx++uk6q1DGiy++iFatWiEjIwMffPCB7rsffvgB/fr1Q2JiItq0aYPJkydr340aNQpvvPGGbv/evXvjm2++Aeccd999N5o3b46kpCT07t0bGzduVJ7/xx9/xBlnnBG0/YMPPgh63n///Tfi4uJw+PBhbb81a9YgPT1d6W5etWoVxo4di5SUFNhsNnTr1i3ouX355ZdITk7GY489hunTpxveJ4fDgXHjxiE7O1v3XIYMGYJffvkFlZWV2jaVR6Vek84YYw8C+IVzvkrYdjKAIZzzFyI+qwULFo4Kj8/ZhM0Hi8Put25/Iap8+lxRt8eHB75Yj1l/hF7yvkdGIiaN7hn2HF26dIHdbseECRNwxRVXYMCAAUhJSdG+55zjoYcewuDBg1FcXIwxY8Zg8uTJeO2117R9vvzySyxYsABerxf9+vXDunXr8P7776NHjx4YOXIkXn/9dUyaNEnbf9GiRdixYwd2796NM888E3369MGwYcN015WdnY1Ro0bh448/xogRI/DLL79gzJgx2Lp1K9LT0zF27FiceuqpmD9/PlauXIlRo0bhggsuUI5x7dq1uO666zBnzhz0798fM2bMwPnnn49t27bhyiuvxBNPPIHi4mIkJibC5/Nh9uzZ+PrrrwEAN998M7p06YKdO3eirKwMo0aNgtPpRL9+/QAAK1euxJlnnomdO3ciLi4ORUVFmDFjhma9zZo1C8OGDUN6etjVF0PC7XZj1apVuO2225Tf//TTT3jppZfwyy+/oH379rjxxht138fHx+Ojjz5Cz549sXHjRpx99tno27cvLrzwQkyYMAEvv/wy7rjjDgDAX3/9hezsbJx77rmYP38+Fi9ejO3btyMpKQlbt241zIjfsGEDRo4MXtXZ6HkPGTIEs2fPxq233goAmDFjBq644go4nc6gYwwYMACPPPIIjhw5gtNOOw2dO3cO2mf69Om48sorccUVV+Dee+/F2rVrccIJJwTtV1VVhY8++gjNmjXTzfXMzEw4nU5s27ZNU1hVqO+ks4kANkvbNgO4K+IzWrBgocEgk3W47bVBYmIili5dCsYYbrzxRqSnp+P8889HTk4OAKBTp044++yzER0djfT0dNxzzz1BrRnvuOMOtGjRApmZmTj99NNxyimnoF+/foiOjsZFF12EdevW6fafNGkS4uPj0atXL1x77bWYNWtW0HXNmDED5557Ls4991ytu1T//v0xd+5c7N+/H6tWrcKTTz6J6OhoDB48GKNHG4cPp06diptvvhmnnHKKppxER0djxYoVaNu2LU444QR88803AICFCxciLi4OAwYMQE5ODpYvX46XX34Z8fHxaN68OW6//XbMmzdPO3ZGRgbGjx+veScmTJiAmTNnalb4xx9/jKuvvrpWz0bELbfcgj59+mD48OHK72fPno1rr70Wxx9/POLj43UWNBCwHnv16gWbzYbevXvjyiuv1J7jBRdcgB07dmDHjh3aNV9++eWIioqC0+lESUkJtm7dCs45unfvjlatWimvobCwEAkJCUHbjZ73hAkTMGPGDACAz+fDrFmzDO/VG2+8gXHjxuHNN99Ejx490KlTJ/z444/a9/v378eiRYswduxYtGjRAmeddVaQlT179mwkJycjNjYWU6dOxRdffAGHQ2/3JiQk6DxIDW5hI9AzXPYxVCHQstSCBQtHiVtvvRXvvPOO6f3NWL4AMOi5hcguDG6dkJkci89uPtX0+cKhe/fu+PDDDwEAW7duxVVXXYW77roLs2bNQm5uLu68804sWbIEJSUl8Pv9OqsECLh2CbGxsUGfS0tLdfu3adNG+3/btm11LmTCvn378Pnnn+uSgDweD4YOHYqDBw8iJSUF8fHxuuMcOHBAOb59+/Zh+vTpOrdvVVWVlpg0duxYzJo1C+PHj8fMmTMxduxY7Xder1d3vX6/H82bN1eOBQBOOeUUxMfHY+3atUhMTMTOnTtx/vnnK6/LLO6//35s3LgRixYtMsxQP3jwoK42vG3btrrvV65ciQcffBAbN25EVVUVKisrNS9AdHQ0LrvsMsyYMQOTJk3CrFmz8MUXXwAAzjzzTNx+++3417/+hf379+Oiiy7CSy+9hMTExKBrSElJQUlJSdB2o+d9wQUX4JZbbsHu3bs1C/7kk09Wji82NhYPP/wwHn74YRQXF+O5557DpZdeiv379yM1NRUff/wxunfvjr59+wIIxKPvvfdevPTSS5rFTmPMz8/HmDFjsGbNGgwZMkR3HjM19fWdJb4GgOxHuQXA2ojPaMGChSC8++679XLc+4d3RaxTv6BerNOO+4d3rZfzAUC3bt1wzTXXaHHKhx56CIwxrF+/HsXFxZgxY0ZQRnSkEIl1//79yMjICNqnTZs2uPrqq1FYWKj9lZWV4cEHH0SrVq1w5MgRrfMUHccIbdq0wSOPPKI7Vnl5Oa688koAwKWXXopff/0VWVlZ+PrrrzXCbtOmDaKionDw4EHtdzk5OZo1DgQSzuQs8fHjx+PHH3/Exx9/jEsuuSTiFpYiJk2ahB9//BHz589XkiShVatWQfdVxNixY3H++efjwIEDKCoqwi233KK75gkTJuCTTz7BL7/8gri4OJx6ao1CeOedd2LNmjXYtGkTtm/fjhdffFF5Db1799ZyIkQYPe+YmBhcdtll+OSTTyLyRCQmJuLhhx9GWVkZ9uzZAwD46KOPsHv3brRs2RItW7bEPffcg/z8fJ0VTkhLS8N7772HyZMn49ChQ9r2gwcPorKyEl271rxfdWlhmyXsuxHoJb6GMTabMbYWwP8BuDPiM1qwYCEspk+fjpkzZx71cS7sl4lnL+6FzORYMAQs62cv7oUL+9VdFu7WrVvx8ssvIysrC0BAuM6aNQsDBgwAELA4qBQmOzvbUFhHgieffBLl5eXYtGkTpk2bhssvvzxon6uuugpz5szBvHnz4PP5UFFRoZFq27Zt0b9/f0yaNAlVVVVYunRpUDmOiBtvvBHvvvsuVq5cCc45ysrK8MMPP2jWYHp6OoYMGYJrr70W7du3R/fu3QEESPDUU0/FAw88gOLiYvj9fuzevTuo5E0m7HHjxuHXX3/FjBkzMH78+JD3oqKiQktyqqys1LU4ffbZZzFz5kwsWLAAzZo1C3mcyy67DB9++CE2b96M8vJyPP64Pqe4pKQEqampiImJwR9//BE0P0899VTYbDbce++9OuJctWoVVq5cCY/Hg/j4eMTExMBuV6/KPHLkSF1iICHU8x4/fjw+/PBDfPfddyHr6J988kmsWrUKVVVVqKiowH/+8x8kJyeja9eu+P3337Fr1y788ccf+PPPP/Hnn39i48aNGDt2rGHyWbdu3TB8+HC88EJNGtevv/6K/v37h62QsNvtWsgjEpjNEt8EoAuAFwGsAvACgK6cczmubcGChTCYPXs2cnNzQ+6zY8cO7NpVN1WTF/bLxLIHz8Se50Zh2YNn1ilZA4GY3cqVKzVX7oABA3D88cdr2byTJk3C2rVrkZSUhFGjRuHiiy8+6nOeccYZ6NSpE8466yzcd999OOecc4L2adOmDb799ls888wzSE9PR5s2bfDiiy9qgnLmzJlYuXIlUlNT8fjjj4ckxv79+2Pq1Km4/fbbkZKSgk6dOmkhAMLYsWPx888/a9Y14emnn0ZVVRV69OiBlJQUjB07NigjXSbszMxMdOvWDYwxnH766QAClprKMxEbGwuXywUgQCKxsbHadw8//DD279+Pzp07w+VyweVy4ZlnnlGOceTIkbjrrrtw5plnolOnTjjzTH2vrLfffhuPPfYYEhIS8MQTT+Cyyy4LOsb48eOxYcMGHXEWFxfjxhtvREpKCtq2bYtmzZrhvvvu08YkYtSoUdi4cWNQDfTgwYMNn/egQYNgs9lwwgknoF27dsqxAYF7fO211yItLQ0ZGRlYsGABfvjhB7hcLkyfPh0XXHABevXqpVnYLVu2xMSJE/H999+joKBAecz7778fU6ZM0d7nTz75RDm/ZQvbZrPVirC1SfBP/TvxxBN5fSMvL4+vWLGi3s9j4Z+BQYMG8cWLF+u2BV7FGjz00EP8qaee0j5v3ry5Qa6tsWPPnj0cAPd4PMf6Ukxjw4YNvLy8XPtcVlbGN27cqNtn9+7dPC8vT/vs8Xj4+eefzx955BFtW05ODs/Jyan/Cz4KTJ8+nQ8aNMj0/uvXr+c+n0/7XFVVxf/973/ziRMnats2btzIKysrQx5n6NChfOrUqRFfb11i/fr1fMCAAXzVqlW67UeOHOHbt2/nnNe8xyNGjOBz587V9gGwmpvgM0MLmzFmXISm329a5GrCPwvLly/HU089dawvw0ITgZn4VWVl5VHFLS00HphpQypb2Hv27MHChQtx/fXXa9s8Hk+jbtFZXl6Ot99+GzfddJPp33i9Xp2l6ff7cemll+pK/nw+X8ich1WrVmHt2rXKsEhDolevXlqNu3y9dWVhh3KJX8IYG8oYOzPUH4CLIj7rPwxVVVWIiorSbYtk0lr434LX69WVgqiEkcfjCSoXsdB0IT5jHqaX+KOPPop+/fphwoQJaN++venfHUvMmzcP6enpaNGiRVBIIBRULVlV5Cbv43YHKh8mTJiAYcOG4bXXXlOWgx0LyGNSvd+1JexQEiEPwAchvif8HfFZGwDVPVlHd+rUqc6P/dBDD+HZZ5/VPns8niDCnjp1KqZMmVLn57bwz4OqdaFKkFkIdEVravfFzLMUSerJJ5/Eo48+ii1btuj2acyEPXz4cF3GvQperxcHDhzQKSFmyI0xFkRumzZtQv/+/UN2IztWMPO8FYSdxBibAmAOD9Ge1NDC5py345y3N/HXLdIBNQR4PfYSf+6553SfLWvIQiRwOBzw+XzaZ5/PZ4qww7XOtHDs4fP5gsqhVKSkIt6mpohECp/PF1RjbcbCNkPqjQlmnreCsOuul7iF0FC5bIDGP7EsHBvY7XZdLNLn8wWVuahcZvJKTP9ENPV3xufzaX3PCbXxljRma7q2MONJMnIfi9v9fn+jvjdGz1t8n+sjhm3BJGw2W1AyiN1u11lRFiwQ5JV6VIJMfqHj4+NRVlaGqqqqJk9qobBz505dHTGAJvUemQ1vmFleU0ZTD5MYWZrhrFHZJc45D7rHjQnyc/L7/eCcIzs7W+usVx8xbAsmoZpAlAlsucotVFRU6DK+5SxxlYUtv/StW7fGgw8+iIyMjEadKRwpcnNzdW06Dx48iPLycl1OyL59+4LaZDZWVFVVIT8/Xxd/zsnJQXl5uTYHKioqUFhYqHu+BQUFcDgcWr2vx+NBbm6u7jhHjhyBzWbT9aluSqB7I8rKnJwceL1e7Xkb7VNWVqbVl/t8PuTl5enuTUVFBZxOp2FDloZEbm4ubDabJvtLS0tRWFiIPn36IC0tDYBVh92gddiQamanT5/Or7rqKt22+Ph4XlJSUufnttD0EB8fr/t8wQUX8K+++kr7fOTIEZ6YmKjb58EHH+RPP/20bps875o6/H5/0Jj69OnD165dq9vWlMb9+++/85NPPlm3bdiwYXz+/Pna50WLFvHBgwfr9rnrrrv4K6+8on3evn0779Spk26fhx9+WFeb39hxyy23cL/fr31evXo179evn26fk046SdfDYs2aNbxv3766fUaOHMl/+OEH7fPBgwd5ixYtdPuMGTOGz549uy4v3xS8Xq+uxp5zztu3b8937dqlff7www/5+PHjdftcccUVfObMmdpnHG0dtgjGmLLNC2MsdLumfwDmzZuH1atXh9zH7/cHaXayFZWVlaUtAWfhfwty9qwZC1uxwP0/DhUVFUEtHJ1Op25RBK/XG+S9WrZsWaNNwFNVjMjWlEpeyPuoEhFVobfGjKlTpwYlV6rmeVVVVch95PCiah/VPW0IfPHFF7j55pt121TPO1zIyyzMBgKCFhdljDkBHHv/Qz3j22+/xcqVK0Puo3q5ZKHsdrvx888/18s1Wmi88Hg8YclYJWxUpSy0b1PFY489plt1q6ysTLdaFhD83pSWlmptNwnPPvus1qCiscGIaCMV4EakzqUY9tGu4lWfMEvGsqJSG8L2er3HhLD9fr9ujEDgOYnX22CEzRhbwhhbDCCGMbZY/AOwDcDyiM8Y+nzRjLH3GWP7GGMljLF1jLGRwvdnMca2MsbKGWOLGGP1HtiKiooK0moZY0EPJJyFXV5ejri4ONPnveaaa2p3wRYaFcyQkkrIA8EZs9HR0doiD00RH3zwga4ns+qdkC1sVVMi3oiTjsySsVnCFrepks5CLVhyrBEVFaWbr2YIWyVLzRB2rWPCR4mqqqogL5FKCakrwg6XEfVfAAzASQDeF7ZzADkAFkZ8xvDXcwDAGQD2AzgXwGzGWC8ApQC+AnADgDkAngTwGYABdXwNOqhuLAlcmjRmLOxIE9CmT58etLgA/weWevzToSIluazLyMKWUdsl+RoLYmNjdQLc6/Vq6wwTZMJWCefG3PfADJmYIWyfzwen0xm0b1N6/6Ojo3XWp9G4w/UkUN0b+R4fK2W2srIyiLBrEwIxi5BqKud8Ouf8QwD9qv9Pfx9xzudxziNfgTv0+co455M553s5537O+fcA9gA4EcDFADZxzj/nnFcAmAygD2OsXhu3qG6sGTeOWSsKAN577z28//77yu8IGzZs0NbetXBs8dBDD+GXX37RbcvOzlbu6/P5gsjFjAaugmxp7Ny5Ew888EAkl35MoSJjedxmCLsxx/ePJj4tC3kibBGyhd2YIRN2fbrE5XMBCPpcH1B5e8y6xGtTrmhKTeWcb2WMnQOgLwCX9N1jEZ/VJBhjLRBY1nMTgFsB/CWct4wxtgtATwBb6/Eagm6s6oGoCFsUPKGSIg4fPhzUAUiGx+PBjh07Ir18C/WAnTt3ol+/frptrVu3VgpTs8LZTPxNFlz5+fn47bffIr38Y4bakLFZa2rQoEFYtmxZPVx1ZKjt8zaysM0IddU9OhaQPYBRUVFhCVtlYdeGsGX3OwC0bNnScFnMuoJZg06eE/W6HjZj7E0AMxCwdNsIf60jPqNJVCe1fQJgOud8KwKKQpG0WxGAoI7vjLGbGGOrGWOr8/Lyjuo6juaBhBM8Zr4jVFZWBrkPLTQMZE3drEUM1N7iUkHV0rSxuoYB4LbbbtN9lpVYI8+UnCUuj1EmAiCwYl5jgJnnbTaGrbKwZZd4TExMUKOZY4Hdu3dj+PDhum1mCLuuYtgqJU7uOFdfkBV1MwadglfSiLOq/5SrR5l9268E0JdzfsDk/kcFxpgNwMcAqgDcXr25FECitGsigCDTlHM+BcAUAOjfv/9R+ZBUmZlmLexwkyzcd6LG6na7I0pas1B3aNWqFQ4fPqx9ttvtpl2TtY3JqSArgY25MY/H48GUKVPw9ttva9vMWtjh9lERdmNBbZPOVN2xZBmigiop9ligtLQUhw4d0m0zk3RWVzHsY9VZ0oz1bDJLPJ9z3j/c+cymWh4GUGhy36MCCzDU+wBaABgjxMk3Aegj7BcPoGP19jrDHXfcofusijXUxuWhemji71UTMFyWeXl5uYkR1Q4lJSVBGur48ePr7XyE++67T/d59uzZWLx4cb2fNxRkt1ok8SezFpfD4QipBHDOgwT4sSplMQOjDPBwLVnNEHZjzpavbVmX6jhmYtiynDhWkLv5AcEkamQ910UM+1gSdl14VMzCLGG/DOATxtipjLEO4l/EZwyPdwB0BzCac+4Wtn8N4HjG2BjGWAyAxwCsr3aX1xnefPNN3WfVjTXr8ginOYb6vaq8Rc5GbN++PYqLi02MKnLMmDEDkydP1m37+OOPdZ/LysqwaVOd6kt4+eWXdZ8XLFiAbdu26bbRWrgNgcrKSmXSWDjh4Pf74fV6TVsMDocj5Avs8/kQFRXVZFzibrdbayVJMOsSDyec5eMcSzz00EO6z7W1sGWQSzzcPGsslQO1Ta40I0ubImFH6jUwC7OE/Q6A8wAsA7BT+KvTLKjquuqbEUhu+5sxVlr9N45zngdgDICnARwBcAqAK+ry/KobaGQNmbGwjTIFhw4dqttXZVmZnaR1RdgbNmzQfa6oqAgrVNatW4cbb7yx1ue8//77sWDBgpD7qAR/enp6rc8ZKVSNO2RlSgQpE9OmTcODDz5oWgNXWVMiVITdmC1sVQjHjPVsJvfDjKv4jz/+wNy5c2t7+aYhL7Vr5nmbyVkwsrDlGHZjyZhXzUWzCWWRytLazon6gBlFpcHKugicc5vBX51KC875Ps4545zHcM5dwt8n1d//zDnvxjmP5ZwP4Zzvrcvzq2KC8guiEpyRaoW//vorAOCJJ57QWi+qXF3htLTY2Ng6Szjp3bu37rPqmmQcbV343r17wyaGqASB3OqzPqGqFVaVkBB69uyJQ4cOweFwICcnp84ygo0IuzFb2LKL1Axhm7GwVXXs8jxcsWIFfvzxx6MeR6QwGwIJp2iZUeKA4HtRUFCAGTNm1OLKjw5GxBXuWZohdbPEfywUF/laqMzrWLvEAQCMsTaMsXptVHIsYUYAql4kMzFs1SR77rnnUFVVZZiJHm6SxsbG1pt72MjtKJK4yvqNxFXp8XjCZr5HRUWFPSY1xq8LXHzxxUHNPSJxyZaVlYExpik8Zt2AtbWwG0vlgFyWqHq2MmEbKbrhLGz53VA9o4Zykcpln3UdwzbjEhf32b9/P1566aVIhlAnUPV8N+slNBPDjjSTnO6xuG3NmjV4+umnazE6Y8jXqwplNLhLnDF2HGNsGQL1zj9Xb7uEMfbfiM/YiGGGQFQvktnieaPaQaPEtnDWVG0Ju7CwsFaxsejoaJ1FX1FREUTYLVq0MDzmO++8o/usajspQ1VfKWP69OmYNGlSyH3MYsWKFbqFJcxYDiLI+ibCVglnlZCqTQxbNScuv/xy5OY2/Jo8nTp10nlL6tP9aebdaKhMcjkBrrZlXTLMWtiqNQvkd7I+ICvIZmPY4axns95KlbwVr8nn8yEuLk43B/Ly8uo8gdWIsI+pSxzAewB+QKDmmVTkBQDOjviMjRiRELYZV1W4SUZCRdUj2IwbSFWDOWBAeAfI6NGjwy5oAgS/lLIQVAnKUC5uVV1uODeWmRhdSUmJruzqaCDfU9UYQzU9EJUwv99fpzHs6OjosNbnX3/9VWf3IhI4nU5dqKK2rk2zLvFwhN1QFraKsOsq6cxMXFZ+hxqqbWvHjh11n1XPwCwZ14USJ8Pr9SImJkb3u7qYE2KZoupaVMq3KnRY34R9MoDnOOd+BPqIg3NeBCAp4jM2YphxMaosHRXMlDSQNaay0M1opyrCVhGxTHjkqg0F1b1QvThHIxxCJW8RzJB6XcZyZS+CmWQa8TnRM6XVtho6hp2YmIiiIrm/UP1DJq6jcX9G6hKnexPu2K+99lqdt/ZUjVtVY11XSWciqNRPJuyGCJPs2bMnyKI1E8MO54k0q+iFK4P0er2mFNxI8a9//Uv3WVa+zGb31zdh5wDoJG5gjPVAYIGOJotzzjlHV8tsRjs1ImxV0li4lzSU2642rkHVBN61axeGDBmi2ya70lW/M2O1GGnVKoJVncMMYcsC3Kj9p3xvdu/eHfK4RoiJiQkSwGYEEQkQeqbhYtgqwg71AqsEgRFh11epXyjUpg1lbRsOqeZhTExMWOKaNGkSCgsLAQDfrMvGoOcWov2DP2DQcwvxzTp1L/hwkMfdkElnfr8/qHFKQxG2GcW2NjLMrKIXzmhSWdgNcW+M3mWVB7U+CfslAN8zxq4F4GCMXYnASlnPR3zGBgJjbDRjbEooa2Pz5s26phhmLWwzLnEzmuPRErZMeKq1lw8fPhxkhcuErXrZVGRgZgUyOXPd5/MhKysLlZWVQfFqsxa2GQ1cfiFkl50KO3bswPz580NevxmXuOiurk+XuBkLO1QGe33CjIVtprzJjIUtuxfNErbL5UJZWRm+WZeNh77agOxCNziA7EI3/u/L9fhg6W7kFFfgcGklSio8qPD44PNz7Nixw9CaM2NFGiVNhXveZrwu0dHRun1+P+jF7p7XHLUiEg6yDDGr2NZGmTHyqIQyVlQWtpmcmaOF0btswiWexBibwhgbHer4Zhf/+IAxVgDgJgSWv5wA4FHO+Tdmfn8swDmfA2BO//79DQuF4+PjdXG3o7Gw5QdixiVeF4QdTkgVFBQgNTVVty0uLi6IsOnFp3OYtbBVbnq3240vvvgCQ4YMgdvtxkUXXYRVq1Yps4bDubtVSgK9iHR9qpig0+lULn0n4o8//sAPP/yAc845R3f94QhbZSmSpUPEJbrEzbhI6yrpTDWnGsKyMNs3ui5i2DJUhK1SvuPi4lBeXo4X52XB7dG/v5VeP574fgue+H5L0PE59yPauRNRdhscdgan3QanjcHpsAHnPorrP98FV1w2HDaGvJw0lCcMxLXT/gjsZ7fhr7h+OJCXiP1fbYDTzrA5PxGFcT3x8vxtiHY64LAzbPK2gA0cn6zcB6fdhj9zgYP2lli6pwRZ/hw47TYc8iXAzjj+OlAIp90Gn6cS9qSW+Lu4ErnFFViwOQczdwKeqEAH5+xCNx76KtBf4cJ+mSHvYaQwo/SbDYGYUXjkfcIpuKJMI9TXeyDKfrO5B7IMAFDEOVf2DxdhOvBXTc7fmN2/KSAuLi6IsMM9UDPWEOdcm2Qvvvgihg0bppyITqdTSdji7wlGwi0cYVdWVgbVw8bGxgaFAsh9FIqwzSw+QVbzV199heTkZLRu3RpxcXHKfc10rJItLq/Xi9jYWN31qYSzy+VCaWmpjrDPPfdcXTMNVRKKqv+xGcuByJSEFP1fvKcjRozATz/9FKTcGcW95BhhbQk7IyMDR7sITjiY6RttNmtYpcyEqvlXWVMq5ZsI+2ChcWXFmLZV+GPVGky49jpU+fzweDmeee55XPuvf8HujIbH54fHx/HbkqXoe8oAHPjrb7RwdURMbBQ8vkB5oc/mQH5pVfW+fhTak1FSGY2cLTnw+Pwoc7vgc8bj7d92o+YRHwcAWPn1RrqjALpi6aK/AfxdvS2wkvDct4RVyXpfhzsWFAELfqneoL9Hbo8PL87bVueELSu2qvfETMcvM3NCTspVvQsqr4tqTshy4tVXX8Udd9xRZzkwZj2wte1fYfoqGWOnA+iH4OU1n6nVmRsBZAvbTPKSGVeVSOqLFy9Gt27dTGmX4u/NdLWSCc+IsGUrU7ZsVfEeI405nEucxkQWZ1lZGeLj45XHMxvDDhevVF1HQkICSkpK0KxZM22b3EjD6J6Guw8qUqHnRUJKJG/6/bx585RjNPOSq4SUak6p+myLZWpAYH6tWbMGJ510krZt6dKl6NmzJ1JSUrRt9Oxk+P1+/PnnnzjhhBMMz1vbGLYM0VIyejd9Pp8pl3hcXBx25pbAZmPw+YNd3JnJsTg13Y29RVtww+k1XZefu/oH3HTKE8jIyAAQKGdMurI/lr9aifmPj8WjD1+EHj16AABef30lduzYgTfuuEr7/cSJE9EhswMmTpxYvc/rOHDgAJ58/Ek4owJKwBNPPQ0/bLj73vvg8fnxzXdzsGffAQw+YyjadegIj8+Pqe9Pgx8MV4wdB4+Po7C4BO+8NwVnDjsHnbp0wSMa2euRXViOoqIiJCXVXY6wSoaEk29G8f1IPSpmYtg0J8JZ2M8//zyuuOIKtGrVKuQ5zcJM0tnRJD6arcN+A8AXAAYj0Oeb/rrV+swNjAMHDgT1vVZZmnURw5YtLq/Xa6oeV/V7cZvKmlW5i8VjqnqQq4heFnhGmetmSN3n82kuH4/HownccAs9qCayHJ8mwpaXYJSvg6ypUFB5H1Qx1Ei636ks7NrW3IrnEePkhNrWH//9998477zzdNsefvhh/PXXX7ptmZlqq6yiogKnnXaabpvs7amt+1OGGeGsUuKUllp6J0xaUowYhw3RDqlroNOO+4d3Vf5OlhPl5eWaIqNS3szEaSkEYrcxxDjtiGJ+xDAvWiTGoHVKHFKdPjSP9qFtkh192iSjf7tUtLKXINNejLO6t8CI41virM7JyKjKwsCWwLhT2iIzWV1/zSqKsWfPHsP7VxuYeU9qk1BmpMTJ74KZOWHGwpaNNjMIRbhm+EElW83C7K/GAejHOb+Ec3618Ff/yzcdJbb+XYJv1mVj7ty5+M9//qP7TiYMcq2GeyBmtDtRgMsWF00+o1WfzFpTKpd4XFycbkxGyV7hLGwVzC7MQF2P/H6/RqZGLvFIyy5U1hQQTKJyhzTVC6TyPpjxIqiaRpixsEONqT5j2HI3roqKiqA+36rrlBM233vvPbz33num8htq6xKXYTZeqZoTIn5YfwiHul2CeAcwd+LpeH5Mb7RKjAbnfmQmx+LZi3vhwn6Zhpa5HK+lfcxakbVJMjRbm0/nv394VzhY8P7cGYvfdtZ+bWjOeVAnOzOeKLMhEDNKnOwSNzsnwlnYZhR7EeHi06IyFmqf+ibsAwAa53p2YeDx+fHQVxuw7rA9qAOQLNTJGqyLF0kmbHqRxUlr5BJXncOIHOXrj4+P1wlsVWZkqESuUIjEJS4SNp2P9iVLzuilD1dfKVvYKsjEpSJn1ZKAqtKpUKRC5EzP62gsbLNKoHhtZghbtiJU98LMkpU7d+5EYWFhrQnbyP0pPu9QypARjOYvxT6nLt6Nf81ci1h3Hu7p5UfbZvG4sF8mvrm+F9zTrseyB8/UYrwqj5SceS/OCRr3woULcejQIdONU8wI9UgVlQv7ZWJEymHE+MvBEHDxP3Zed0RXHsGrq8vxyvxtylBAOGzfvh2nnnqqbptK2TYTwzZD6rUJD8kwmyUeqYUty0AZZt/l+ibs6wFMZYxdyhgbLP7V6qwNDLfHh2WlaUHCQFUWRQlNRjBrYYsCnIQ6TWASUvJk5Zxj6tSpynOo3CiqLHHZwo7EJR7Owla5xOklnTVrlm4fEpa0j/iy9u3bV9tXfunlBhgyxKSzUJAT+lQkpdpmxuona/7IkSPo27ev7iWlZ0yCW7a4jEIgdbVal0zYHo8HLpdLlyCkGreZhWSo7aV43uefD1R21tb9GYm3KtQ+Kgvbz4HJ323C03O34NxeLdE5+ydEoWbeq9p4qoQ6zdMtW7bgrbfe0iW00bjffPNNLF++XCmMa0PYNCciJaUuMcU4270Ye54bhWUPnonrTuuAbllzcXKaD68v3Inrp69CUXlkS5MWFxcH3SeV8l6bqoC66qsvw6yFLYc7QoGuLRw/mKmfD+dZMoJZwj4RwEgEltn8RPhr+GVhaokSvwOLYk/DdR+uwuTvNmHasj0oSWiLgyU+VFSXeJhxrUUqXImsRBcp/V6erKWlpbj77rsNCVuGyiUuW9iRuMQj1RxFwh47dqx2TSoLW35ZVZnwZlzzZq4VUFvYsjUdSjiL5zPKaygpKUFRUVGQR8XIwjZayORorApVKEC0lCsqKpCYmKg7r8qzoOqaJ4PGJj73Bx98UPuuNhZ2uGzfSISz7rlxGxZ5OmH67/tw4+nt8eaVJyDaYdPdh0gJOysrC19//bVu/DRuKnOShbHR6k11YWGrFBUVKSXExeC85kV4+qLjsWxnPka/uRSbD5pvrlNcXBy0xKyZihEz4QIjC9uMQRSpMqO6N6p3UuzNIR8znHwy89yOxiVuNkv8GQCjOec/1+osjQBO+JBsq8TBQjdW7D6M8iof0Hwotu4A/vPYT2iZGAMXPCjMHIp3F+9Bl4xUtG0Whwq/DTG24PiT/NBUMRYSWpR0JmZQi8KdQORgRkgBaktZZWGHIyUzFjbFy1RuMHHsdGzZwhaFnNhcRE4oExuQqBCJS1wmKdmqVHktQrnpZTKhJSTF5yXHsEXiMirjq63FpepBL88pIuxw3gYzhE0eqbqMYUfirTKCTNj5pZX4qaor8nk8Jo/ugWsGtdddP8HtdsPlcumutaqqStkvgK5DHr8cApItbGofqlJU5DGJz9Lv9wfVEavGbYaUAi79StxwSlt0b5WIW2eswcXvLMPzY3rjgr7hy71USqsqrCYrKqpxy/dGpcSZVV4jie/TNRqVoYpo1qyZYfdHVXKuCKM6bHpXKZ+kvgm7DEDdLnPSgIh12nESdqO9LR+P33U1OOfIL63Crfc/iowuvdG+18nYX1COP3dmoSiuNd74bR+AfdW/PhFR8GLtm0txXGocoj3FyIpqgy35XnQrcqNFQsBSkYvnSRDRgxKFuSzcCaI1apaww7nEVZPUTKcoGWJ5jThOGgNNQBJgZMXb7fYgC1tu/xnqHKp9amthh2qiQlC5xEk4yASkUrBkC1v8nVEGt1k3mhmvi7y9oqICCQkJQYmIqjitGQvb5/MpBbj4ndfrrbWFLcOshU3EtSuvFNdOW4UjPBZnOnfimkGBbHjOuZKwyfsgErZRzgc9a1XSGb3nqhCIrIzVVdKZ6n0JZ0WecFwKvr/jdJx017uY+Kkffx4oxMPndofTbkwgRgpaqNCR6vplBZlIXVbiYmNj62ROqMhVln1m1ioQjxku18foudHc16oDaukSN0vYjwF4jTH2BADd2n3VC4I0WjjtNjx7cS/8+e0KUFMBxhjSE6KRzkrQLaYYt53dBQDwxRc78csvv+D/HnkMbns89h0ux38//QYlPBpJsS2xPqsIWUfK4bd3x5o1VXhhzUJEO2yI8R2PBFYJNmcT2qbGIYFVotyRgEqvOunMyMIWrTIzhK2KYcsucRpvKJhx9YS6JlFIiIRdVVWlJGwjcjBD2OQGDGdhq2LYMTExIRtwAMYJeSoCUilYRha2KmeBYDZRxUxiIKB/3mYtbHkuEQkBwYIulIX9zjvv4MCBA+jZs6epGHZtku1k0DPamOPGK98vh50xjIjahnSbvsdCbGysbt6Ul5cjMTExrIUqVj7I3iKZsOV7Jb4L4rhVHhVZ6Tdzb8yUs8nnT0+Ixr4P78Xj36zHB8v2YFN2Md4c1w/NE/RhknD3JJRL3KwyJhO/isRVv6uNS1wFWYkT573qmGZDh/L1i3O/IVziH1T/e7OwjSGwctfRLX9ST6juyTq6U6dOuLBfJv78NngfVVlXTEwMomwc7VolonurRKxwBLoMTb7+WgDA0uW/48ufFqFNt35I79AD+w+XYd6ytSjh0fj0jwNCy8M+wB6OmJTh2JFtQ5w/Bg5bNOZtyYOzeQcUlgWvgx2phS3/XmVhm4GqDluGysUk/p5eaLp++pesbPGFNrKwzZCSKulMFpLkTpRJiixN0YJSuZSNCFt1LSolTCRsVUKaatx1kWijgtvtDrKwVTHsUF3zTj31VMybN08ZwyaILvGysrIGtbB9Ph/+drbCV0sK0To1Hh9eexKmvbFMt09lZSVcLldQuCg+Pj6oF4CRFWaz2bT5IRI2lWqqPB4qayqU8iaOyayFHSkpBS7Mh8dG90CfNkn4vy/XY/QbS/H2uBNxYtuUoN+rnnc4l7hZ5VuVdW30nmiX7vcjJiYm4rIuFRHL90ZUZqdOnYouXbrgjDPOMDymakwqZULMYzAg7CTG2BQAc6rbaithlubbV/91EP7oc6ME53wO5/ymUN19VFni4YiLcT/Sojm6Jvlx9YC2eGRUD5wZtRMXRG/C5ieG48CbV+GtiztiaOx+nJpYhFRfATx+hj2eBGywdcDD3++C96z7cNqrK7CpywS8vN6G36o6YK0nE99vyoe9ZRfkFFfAKQmp3b5m+Lyit66pvzwBa0vYkbjEVfvIFja9vOQeFa0SWtFKJOz9+/fjoYceqjOXOHkaZJJKTEwMW7pkFG80IlqVhS0nnckVA4QDBw7ghhtuqLPMWBXMWtihCHvv3r0oLy/X9glF2LLCIu9T1zFszjl+2uvBb97OaJdkx5e3DkTbZsHd2aqqquByuYJKHuV5AqjXBVDNZaCGXFSE/fzzzxtaU+FWa6pt0plK4TAKxXDOcUHfTHx92yBEO+y4/L3fMWPFPuVCGqoQiHxecZyRWNgyjMpdIzm2WWVGnvfiu7F69Wps2VLTW/5oLGxRCSE5IaGIc35TKLIGTBI253wf53wfAvXYVfS5eluThVFZV20FCGMM/rJCHN8iFj3jy3Cq6zBOqNqAK9OzcVX8Rlxm/wPTruiCuHUzccfpbZBStg/RNo487sIGXys8tygbbNi9uHTGDvycPAqTV3pxw/TVmPDBSiz1tEMZorXVhR76agN2efWLehi5xMPBbNIZvSQLFy7UNdVQxfREC1sV5xXPdfjwYcydO9eUBqtybaoUF1W2fEJCQthwQSiXuOpa5Mz/UBa2TFhVVVVYuHChLrZlhKMl7HAx7FBtbmmVKznznyBm/IvLiqrK4WpD2IYKk59j8neb8NVuoFN0Ce47MQqp8epkRSJsMxa2DHkuq8q6VIT94IMPhiTscFZkuCYdRqQUrnyV4vn0LnRvlYh3LmwH38FN+Pc3G/HAF+u1yhk6T7h1BWSY8RAYxXJl93FtKwfCKf9A8L0RcxjIE7Vt2zb88MMPppQAo+cmj6ley7oYY8mMsZkAKgDsrN52PmPsqVqd9Rgh3ERWuVtlmIktqYS0z+dDjJ2jY2o0XEd2YNyJzdHh8Arc0MmNMVF/4aroNZg2pi0cy97DvwakoYM/G2kxHPsLyvDb9nz4pUfl9viw1ttat83Iwg7Xu9aM5ii6xCdPnox169ZpxyUBJo+brBISauIqW+L9q6io0O67GZe4WQtbRdiyhS3fGxVhG2WkimNRPW/ZwpaFNJ2rtrFcxphq1R8dVK5gMxZ2VVWV5nYM1XNftqxVhE0lf2YXPQk3biBQtnXLjDWY/vs+DGnpxUXNC2DjxscxqqCgPvehIHsWQiWdyTAi7Lp0iavcvuK1yNUYlZWVQd6mWAcHFr+LO8/qjM/XZOHSd39H1pFAfbJqIZVwyVqhSJWuxcDSDLJGKQfEzLEJtQ0XiO8GeSu2bduGd999NyKXuMrCrosYttlfvQugCEBbACQFfwdwea3O2kggu6XMxnIjcePR5BNLGFRJZ3bG0dJlR9ThnTivayJOtO3FDV19mH/3GTBKkSqD3pogoVQfSWfiC0AkRccla0w8jmiV0MstJm/JpCjvIyJcL3EZRi1aZdewfG9UjRFq4xIPZWHL7naR8M1Y2KoEnVDzVeVZMOrwJrvESYmSV7UTCVtMLqQ5ISosQE1THaMxRWopubkDP1V1xc9bcjB5dA+cm1GJ2NjwSpzsmamqqgqaJ6p3he4xEVcohUWGGZe4EdHXJumMrjfU7yoqKpCUlKSbE+Xl5YiPi8U9Z3dB7hdPYE9+GUa88it+25pjyiUuw+hdFmVeKJe4vI8qkc9MuEA8v5leFmLfCmpJK77TR+MSF8dU34R9FoA7OeeHEEg0A+c8D0DzWp31GMEscclCUUSoyUL7hhPSsnCn61IlMQFAhkFT/1imnzhGFnY4K8xIUTHSauUxiYRNpE4vtIrU5N+LbnPVSy6/rLWxsFWxXPnZqkID4VziqrH5/X5dzaXctlQet9l4peo6wglnlStY1URGlcBIhC0rCirClmO5RsJYNaZwCrJ4b3bnlWJuVXcc4bF476oTcc2g9tqcCKfEyV4Foxi2DNHCpnkuE7ZoqYnvmxmXuMpFejRWZLhuXG63G0lJSToLW2wiw7PXY/b1J6A4NwvXTl+NZYdjYZNK1WrrGVFZzzLM3Jv6yhInzxLnXOuCJoaCzBC26p1Ujak2q3aZJewiAGniBsbYcQAORXzGYwQV8coTIdJYruq7cHFL2X0qk7mKsO8f3hV2SCQGoILbsd+XrG0ja0gWQMosUQFGFrbRalEqwjUiZfE7I5e4aLnUp0tcdgPKCpyqBtdMDFulhIn3zUzdfW1j2HLsWQaRkpxsJVtMqjwAcd1xo9p+sWSPLGtVZ7tIxmS0z+q9Bbj4neXwcDtGRG3DOT1batejmr/iO68ac20I24yFLVuIRoQtKvjyPpEocSoLO9SYyMI2Iuz4+Hg0j2OIXfomTmsbj18LU7DedRJKK73aeWvriRTvU6gYtkjqqvunIkWVJ07cx2yWeFxcnDbGiooK7TrNyCcjg05V1lWbNbHNEvZ/AXzJGBsKwMYYOxXAdARc5U0C8s1RxfFUvcRVCQ9GLxK9yKFKeWQBLltcqrKuC/tlYqBjL+JRqTX1f3hEJ6TbK7HQ0wkfLtujXb8q6SwcYZsJBYjEpVo0hbR6sWGMOKZwLnHRQg+nOcuuTUD/sqpCA0ZJZ/J9kIW6+CymTZsW5HUwsrCBGktLTkijY8hJa7VRFEVrSpWgo7KwzcTQxDmhiutTxr/YIEdW0MT3RGVNyJaSKsOZxr1kbynG/nclUuKicG7UlqAa63AruKmerdmks1BZ4mKSJZ0/nIVI+4Syws1akSpFJdz77na7kZycHOQSJ8ImyzIuyo67Tk7Eaa5cHLK3wIVvLcOFV99cZ8QVSpmJxG1OkHsFmEkQk71HYs4HZdcTV4j32+i4jcUl/jyA2QDeAuBEoC77WwD/CfWjxgwja+poYtiyAFaV8oRzl8sWGwmQjo4CXBqzXmvqf2nfFrikWRba2Aoxec5mPDFnMyqrgolK1eFJdd3hLBTxmuTGJyQM6QWQx0sCzMhCly1sn8+H66+/PuS1hrOwo6KiwibaqH4nC3Xxmm677TbduM0qaLK3hQSL7JExU+YjWwwiUcg15oDawqbfhroPYsWEkYUtWloiYctzmDEGr9erCamFCxfi3XffDRqTSpD5/RyLc514cVkBemUm4ctbByLRpn+GZoSzEWHHxcWFDAWJ1yh6G+Q6bFFxkq3IcKSk2sdMAxFRQRYRLq/ByMKmJVdpyUlaEKZ3VB7O4OtRUFaFtWnDMG/T36ZyXuRENRq37BJX9UKQ9wl1HKPzm/G6hEpEpHeE4tqi7Gnbtq2hl9WoDjucEmcGZsu6OOf8Nc55D855POe8e/XnyJ3wxxCyAJctbCMXk7wPCdeSkhKlxSXGK0PVnqpimrLFNmLECGRnZyuvI9Zpx1DnTlwzsB0+WLYHX+Wkwh4VXPIk1mGqhKJZjVkkbPFll2PQqhwA8fcqi022sKdNm6a8HiMLW2VNyR2/VBa2bJmHImxa4EH8TtU4RXymocg8lEtcFcIJJ5xV16/KEhe7cskoLS1F586ddd4mOSYqWpoqwqbx0nvicDhQWVmpCan9+/djxYoV2vUbCWefn+MP73H4Zq8Np2RG45MbTtHKtuT3zowSJ98bVc4HZd4DwP/93//pxqiysMUV+MIR9saNG7FkyRLtd6EsLlWfbRk0b2SICrrqGZOFbeQSp3lOlrbH40GmoxRzbh8Ez+Es3PXFZmyL6oIqjxdVVVWG3pPaJJQZ7aM6TjjCDifTKF4tglzi4u9oxUPxmIWFhTpZQDDylqlc4rWB2bKuBxljJ0nbTmaMPVCrsx4DyJOKCDtcWY8M8YEMGjQIW7du1b4TLS6aZHK5BwleVZmPysI+cOAAjhw5oh2LHjqdi4Fj8vk98dh5PbCzIg5P/16CQrexi0xlhRnF7ulFGjBggI64VDW74Tq0EZmrPAuqOHdsbGzIfuOhrCIj4pI7fqmSzkK5xGXCFl3i4SxscW7I4zayno2uQ4T4LFReI1XDECDYpUj3oqSkBBUVFUdF2KpSN0pMA/Tz0ciaclf5cMuMNdjia4FhrYE7+8cjxlkjuGX3p5mkM9nrIl+LeB8A4IUXXtDdB3HcRn0H5DGJwnnZsmWYMSOwwGE4i4sxFrb0i+67DPF5qeYEJZ2Jc0IkbOpESPOdxp0e70D+7EdwQa90bEIbfJ6Thhtum4gvvvgi6BpChQLMhAvC7SPeYyOFIVwCq2olw1BlkKLVbrQcrVFCWUOXdU0EsFnathnAXbU66zGC+LBU60QD5hsCkHA2cpEaPZBILWyxjaeoTMgv63Wntcfw+P3Ye8SDDw6kYldeqfadKJTMhgLECbdy5UrdNalKgOTmKDJEog2XdEbjrqioUFqatYlX+nw+ZcmbOFaj35GlqLKwxeel8qgQGRNhqBQ0swlG8j5iDNno+lUxbCCYqMQSPfLIkKVh5BIPRdiAvk2taGGLx1MRtps7cMXUFfh5Sw5OduzDZV2c4GGIS9U/QVRwKelMVoSM7oPq2LKFLTeM0Uo0DQg7JiZGmz9mLK7aWtji/VV5ElUucTGGTSEvmgdk7VdWViLaYcNDZx2HwbFZOFAZg5XJQ7CrQN1m2Kj+PFQIJJJ9QlnhZi1soxCS1+vV5IJI2CJxqyzsUC78hoxhRwGQ1dcqAOpu8Y0AjLHRjLEp1JFLbiig0q4A8w0BRAtZTiIK9UBEUgtlYZNwFt3P8kR2OBw6AdOGFeDtS7ui0s8w5p3l+GNPYF1XMUmMBJdskaoyKmXPAJGp7HKjiUwvgMpSlBUV8YVQucSNlns0comLUBEX3QeRsBljQYSnmhN03824xOVnSmSsujeigkfkNnPmTN25xZCCShCI8zWUADKbiEjzQ84SD2dhy6WJdN9obFSrTedWuY99Ph+KeQzmVnXHtr+L8e5VJ6KHIzeIAM2uNhdOmZH3kWG32+F2uzUhLRO2KmcDMBbO4jtlxuIyY2GT4iCPSSRs2e1Lq5SFs7BVzXQolntCkhsXJe0DmB3/3ZuMXb5m2piys7O1+Vqb+HSkMWwjwjbyGopKnKzMhHOJ0zyTZQHhKLwGSYyxKdVrYBjCLGGvAXCbtO0WAGtN/r7BoeolLk5A0poiXW5NFM5EhCIpyS5xo9+rLE3RtUb7iIQt7q9K6vB6vTipQzouTd6L1PgoXPXfldjtS9UIA6ghJZVLUTVun68maYgmrSzsiWDFF1p1nOjoaJ37XLTqRAKg47nd7lotBWqWsIFg74PqdyRAZCVC5dJWJRLSOWQLmciNSM3v9+Pqq6/WxmYmWzZcDJuuWyZn8Xey8kbzWnSJezwenYJFVpzKM0MlS+KSgqJLXBwL3a+3334ba/YVYnbBcfBwO2bdOADDq8u2ZGtKNe/DlTcZzYlQJVAulwuFhYXanJQbpxh5lEghpfdGFR4zQ9hmEqvCzQnVOvBiu9p9+/bhrbfe0hE2WdiqkifR0kxjpRhUtgStoiqxxNMBKzzHobyiCq1bt9aecW3i02azxMMRttGcEDkglEtclXQmErcoCzjnePnll00pXwb71F0vcQB3A3iAMbaGMTabMbYWwP8BuNPk7xsFxBWcQrnEQ8VHxZdUzpY24xIXXaRGAly0ZmXCNnKJE6KiohDrK8NXtw5E3zbJWOzpiN9y9eOOJIYtunJVSWMimdPEDuUSp/suW5py4xUjDdZsvLI2hG30O9HrISsqocq6AOgsTHrp5ZI30SUujjucUKIs8VDhDhp3KJe4OG7Rwg5X1hVKeaIYrMolLiqbtM/9r8/ErbM3I9bux7lRW9DvuJpVo0TiUs17ozCJmWcbqoIiJiYGJSUlyqQzozAHjYnmucfj0RSYo7WwVeEh2csmj0mVqyPGsLdu3Yqvv/5at3gHvcuqpiJyLDfO5sMlzQ6ip/1vbPakY/y01bDHp9RZfDoUYYdziatkmpFHinOOc889V9muVixbpWPSPRLDIvfdd5+hd7VBXeKc800AugB4EcAqAC8A6Mo5l+PajRpRUTWrREXiEjdqiycvEylaTOFc4mLcV/xOdrEaEbZR/IpesuS4KHx8w8lobzuM7w/Y8fYfR+D1+ZVWZCjBq4pPq4grHGHT70UtXTyHTHxi/AwAioqKMGPGDMNrFQWBGcKm46qEupEbTx63yiWu6ipnt9t1sWESLqJHRhw/XaMqg1rl/gzn9lVl1Br9joSYXNalUjSMYscEElKyS1wui5q7qwJpF/wfujaPw9WtcoPKtmprTZl1iRsRdnR0NMrKyrRkNSIYUUE1ar9JrmjyLIjZ5+K9oTGFyxJXyRQjpV0ct8qTKFrYZWVlcLlcut+T9UjHkcOI4v32+/1w2G04yXkA57iysOXvErS85j/YnFuhXS8lMdK9EZ9luBi2GeI38jaqKlbkhDxSZsrLy7Fo0SKlSxyA1rVQVGLpfRcrBI5GUTED0zTPOS/lnH/KOX+x+t/S8L9qHKAJJ5KfKhkD0PfINco+NCpvEgVwOMIm4vT7/UGuRiKHcBa2irBFzTPaYcdg526MbGvDz3srcP301SgqD1ZUxNixDFVCWKSETa0qjfYRCUA+BwmaAwcO4Jlnnjkqa0okQ7KSVL9TufGMCNsoy1v+PT1vr7emHlmMQdJx5JyFcIIgXAzbCLK7mH4XziVOOQsUZgmV2EP3jQhMHLfT6YSfAzO3evDpNg/Kt/+O/1zYEa7gx6YkbPHd5NVLtoZyiRvdm1Ax7JiYGJSVlRmSgSovga6X5IsYCpAtbDOWpmiFy+M2kgFy0lmoGLbcV568PBTDNrKwKY9HlHXd48vx8fg+4FUVeGJpMRYfAjgH7r//fkybNk0bU6TWs5nENJXXJVxVBckixpi2AJHcTMdoHQPxvbPb7ZqCE+p6G7Ksy8EYu5Mx9iVj7DfG2GL6q9VZGxgkIERhaCZLXPUiicJZJAByH5slbNFyVTXSUFmaogA30q5lMAZc3j0W1/SIwtKd+bj/x4Moh7HgWrBggW7pTFEohWqxSd6GUC5xGksolzi5nMK548y4eGWFS3xeNK5IXOKyAJNd4qIVKsJms+mIUCQueoZicxm6RpU1Ih6fiF62Is20RxBJWPS6hEo6IytRtDBD5RPQXCaFRbSwmSMaizydMH9vFYa1Zjj87fOw8dq5PwG1Z8yMhR0qhh0dHY3S0tIgUhS9RURc8vVSja8RYdOYZs6cqZQXjLGgxCqbzRZkmYerw6YYtjgnqCeBx+MJInSPx6OtaicrzaKFLScW0v/bp0bj0Ed3o1e6A7O2+7HU0x4+1FTT1EVCmdnjqCoAjEreKIZfWVlpuIKbaFCI883hcMDtdmsGWGPIEn8VwM0AFgM4EcCXCCz8sbBWZ21gkBAShaEZl7iswT766KM6l7j48G02myaUVC+gaEWLFrbK0lS5n+WuVkbatQpOpxMDW3K8P6E/DpV68dpGJ4749QuK0LHvu+8+7N69W9uucomL7j0xXkcvuUzY5E4S3Ygql7joxpLjriLRq0gpUgvbiLBpThQVFeH111/XxigrEWRpiha2x+MxFC6iq1lMvhJdw6SoGVUFGBFVOCtShugaHDRoEPLy8rT7pbKwVWVdpKCF6nglKioiYRdV+rE6aRAO+JNxTe94XNSeIzY2EC+mfUpLS3Hw4EHtOPJ9kBUjuc2vfG8icYlTWCMmJkZJ2OEsbJIFMmGLz5Cud9y4caZi2PRbcU5E4hIXQZ4csrDF76uqqpCQkKBziQPQFDsxj0W+x/R/XlmGO/tF4YKOTuzyp2FZ7AAcqbJpY2qIGLZ4j0UYlbwRYYvzXoYRYdvtgWoCsXoIANavX6/10JBzFurbJX4xgJGc8/8A8Fb/eyGAobU6awODJpocww6XdCa+EJxzPPXUUzriEh8+vaRGFjbFOkSBKCedqWLYokCRY5qyIDGyrOh8Q7o2x6MD4wEGzK3qjiU78rR9yHqntoQEul45JiRaCkSmRgKMrtdoHyOXuLxghdjaVEZdELZoaRYUFOCVV17Rxig/CzEDWPzOSHCEconTPnIIxKwVIVuR4RYVEJPV8vLykJ2drRE93RexDlt2G9NzIneqSEoiyDVMlovNZkORPxqf5LRCiS0RQ507MbprgjbfKyoqNLfvggULcMsttyjvQzhFNVR+gipBS+VKp2dRWlqqm0sU3jHyFtH1ikqs0fOOJEtcRdhG90JFSuK4xe6HKgvb5XJpBCTOdyJ3mvdksYrxaarf5n4/xnSNwVnO7ShnsZh1+Dhk+xI1ou3Xr59p93FtrXCjxT5UyowYGpDfc7G0UuwlIY9bdonffffdWL16ddD1NkTjlDgAB6r/72aMxXHOtwLoV6uzHgUYY6mMsa8ZY2WMsX2MsbHhfkOkY9YlrkrsMiJTUcsiwjaaiDabTRf/M0qeEa1a2fUSiYVN1yGSX0asH5NOjYOLVeLaaaswe9UB3bGTkpJ0LnHZwlZprDTxzSadyfuICTZGiW2iR0KlqKhcw6qSPTl+ZUT0YnYtCVejLHHxOyMBZOQSly3s6OhobQ3ecIQte10iiWFryYnJycjNzQ3pEhdDMXTdooVNRKsat0hcWZXRmFvVA5V+hsG+dWhrL9S9U263WzuOGJM2Y1WKoHtVW5e4qDyqYtg034zmuxzDpuuNlLBVzz9caScQTEpyDBuoUeLlGDYpamVlZbpYrtPpRFlZmc7ClM8lE5fdbkcbexHO8q1BnM2DBZ4uWPh3YP8///zTtEtcdW/CGS9GkD1SooXtcrmC5LoImveigSGPW1QwxKTkBo1hA9gCgFqTrgYwmTH2bwDBTa7rH28h0LSlBYBxAN5hjPUM9QMSZHLSmSzcZMEvCgfRMlSRKXUCoolkNBFVgpsgWqGUjCRb2F6vFw899JBScKlWalJZkS0SozEyagtO7dgMD3y5Hms9mZpwcLlcKC8v17RKkSjF84kucVE4mbGwjRKV5Bh2ZWVlUPxMFcsFzFnYcqc71b2hOSFaJiKphMoSVylhdI9ULnFVDNvpdGoEaNYlrlJUQjXcEH+XnJxs6BIXm6KIvcdlC5vmvQx6J2JiYrAm14/p+xIRzby4yLULLR360jVx3LIiYsb9KYLug1FiHUG0ELdu3YrFixfr7o3KJS6HgkIRNuVayMl2dG/MJFapXOLh8ljEcavqsIGadcxlT6MYwhBdw06nE6WlpTr5JJ/LKJabaKvEmMTdaG8rwNwsB55dkg8WFat8lmKZlDhu1b2JNJ+HxiEqMzQnKioqkJaWpj1vUT6JOUZiLwo5hi27xI3K+Boihj0RAI3gHgAnABgN4KZanbWWYIzFAxgD4NHqrPWlAL4DcHWo3xkRtjyRZUFBpMw515VwqCaraE0YWVrhCJtit2LNqoqwn3vuOVNaJQkuVZw2ivnxwTUn4bL+rbHel4FH5mxHpUffHIQ0RzFrVIZsYRtlS4eysMXxizFs0XITtWF6EWbNmqV9NkPYIuh4MmHTOcT5IZJKVVWVjsAoPh/KwqbnLrvEVTFsUlTIwvZ6vcjOzg6ZJS4rHKq4rAzax+Vy4ciRI1rpkmhhUytXuSSJrpMIO5SFXVFRgbLWAzBrXywyY304N2oL4vxlunGLhE33RhR8onA2Q9gkUI3CHar7sHjxYnzwwQeaFUsxaJWFTddE8131nRgmChUCoe0UChChyiSXY8BU460aE2Ac+hOTDun6bDabZinKrnSysF0ul04+iUoFxXJlGWiz2WDzezHYuRsXt+dYdbASrca/gqySmnfhu+++w9dffx10THreqjUfZFKXKwdEEEmGimGnpaWhpKTEMBFRdonLMWzZJS6SdIOWdXHOV3HO11b/fwfnfBjn/BTO+ZJanbX26ALAxznfLmz7C4DOwmaM3cQYW80YW52Xl6eLYUfSOEV8wURtWiRT0cIORcaAvrxHZY35fDW9w1WELWuV4QhbNW5Rq3TabXh+TG/0c2Thu/V/4/vStrDHBhJOyLKsqKhAfHy8zl0pj4msCSNSp7GJddhi9ilBjmG73e6g+y/i+uuvR2lpadB9CkXYoudAJjfx3hgRtkhOdP2ULFYbC1sW4DJhHzp0CH379g3pGpQtbFW9uTh+MVs2Pj4eBQUFmmVOZS3i8op0L8Xn5XA4gpQ5+VzMZseHG8qQm3kausRX4KYulYhhXt07JHooZAWNzkv7VFZWmrawHQ5HkGtTnm+iR83r9WrlYWQxJSYmKpPOaPzieyrCKOlMzHEQQyCkqITrCiaHSVRjAkInnd15552635EXj+53eXm57j0Vx0tucvG9lO+J2+1GfHy87lnSfGQMOPs4Ox4ZmAhbjAtP/+HBmpzAWP766y+sWbNGu3+RWNi0jyrJjCAaL6p743a7kZ6eruUsiPNG9CyJvfONks7EcYuym2L3Bi7xNOKs6j+lMWxI2IyxM838Gf2+nuACUCRtKwKQIG7gnE/hnPfnnPdPT08P6RJXNcmQXS0kyOQHIiediWRUGwtbdBsbWdhk8ZhxA6ksbDmm5fV60cdxCK9c1hs5vnisST4DOaV6V08owhYtbCPCEuuwRcuD9pVdrSpyVBF2YmIiioqKtDh9pBa2yiUeysKm+2FU5qN6pmTBiIQtxjRlAS4SgcPhQEFBAai9ropwROFixsKmcdHcTUhIQGFhoc4lTklHRCZGwlm0sB2Omo5bjDF4uQ0f7Y7CL/t9aFW0Gecl/62ttkVKg3hv5eft8Xi0+0/WVEZGRq1d4kag+07PlX5XXl6OhISEoKQz8XdGCqqoxMoWtuwSp9g97VNWVoY9e/Zo+4SKYRNC5WrILvGpU6eipKREt39iYiKKi4s1wiYCF983srBlMpbHrSJs0Qtgt9vRLdWGQ9PvQqs44LXVZVjjaQ1mq2nnqnKJq8oZQyXkyfdENSfkpLP09HTNo6IK2RlZ2GIoQFSsSfmme+D3+7XYvYIf8omzqv+mBF0AQlvY75v4+2+I39cHSgEkStsSAZQo9tWgIuxQSSg+nw+PPfaY7mUVXeLiA5GTzswSdqikM/ElF7Vsh8OB0tJSLc5rZM3SZKVziQJcfIE552jXrh0A4OIT2uBs53ZU2mIwI6clDvN40xa2KiNWvh7R3Su7/eUYkSoJSeURoV7PJGgjIWwipKOxsAmiu1z8ThRSdE3i3BCfoZGFTVm3RlBZ2KEIm0hWdomLSWdUh0ukrLKeSbFQucTL/Xb8VNUVm47YcFFbH7q6N8HjCc4kFwUveVRUChop0QUFBaZCQaJLXCV4VSDCJgu7vLwcLpcrpEs8nEdFlSUuE7aYZAgAa9aswTXXXAMguMOiHMOWIWbH0z5EoIRu3bph27Zt2v6MMSQmJqKkpEQjZSPCLi0t1cjY6J5UVFQgLi5O+Z5QfNrv98NXchh3Hu/H2R1iscHXCj+5O8Lt18u8iRMnGiZbqixs2fsgQlRmvV4vrr32Wp0sFC1sInWVW10Vw6Zxy7F7sYy43su6OOftTfx1qNVZa4/tAByMsc7Ctj4ANoX6kYqwabKqtFOv14vp06cjKysLDodD5/IQYz5y0pnK7SlCzBKnfcRJQRPCyNXmcDhQUlKiq49VQU4YE10zooVdUVGB/Px87Xet7CUY4dgEO/fjx8qu2OtJ1AibNG8ZsnCSrWYat9hAQ6WwUOMZulYibCOXOGMMLpcLRUVFQYQtPiMZ9LzJWjCK94nkTcJBdNuKFqLsNRHHTg0wVIStyhIXY7liTNBoLKos8VCEXV5ejpiYGO13sbGxKCws1MZKJE6JRbL1TBAtbNElXuSPxtyq7jjC43BDN2BAak1DCllBEwWvyiVO1yQmctI+ogCUBSu5xMNZ2OKYvF6vZvGGcomLSphZC1v0qMhZ4uQ1U1llIlSkJI+bji+O2+126zxqXbt21b3zQEDxLSkpQVRUlNaKVVaQxRi2kYUteuRUyrvYDMZut8NXVYHbTkrBQMceHPTGYQE7AYf9cdo+r7/+eshkSzlbXuV9EGWh6Fn66aefsHfvXp2CTklnqrwhOg9VL6hi2GJSLFBD2KHazC5ZsgQff/xx0PiMYDpVjTFmZ4wNYoxdyhgbyBirnYpwFOCclwH4CsATjLF4xtggABcACDlikbhI4zGqVaUbm5CQgCNHjgQRtgg5jlGbpDNZaIRyLZOFLZbbhIJoYdO4Sau02WwoKyvTVughNHN6cBZfixR7Jb4vbIVfs7mWOR7KJS7H3uV96MVTucSBGtc/CTuKn9I+clY/5xxxcXFKC9sMZMK+5JJLdAJYFFh0TXFxcUGeBhqTaGGLhET7qJLOVFnioktc1NyNIGdCh7Kw8/Pz0bNnTyQlJemUBzFeXFVVpSUWGYUA6Lxy0lmu34W5Vd1Rxe0YHrUV/ZrblMQlwsh7IT5vEo4JCQkoKioKa2ma8TT4/cHLdJKSILrEZQtbzhIPZ2FTCAvQZweLpYKyoquSTaqkMzPjJrKg33Xp0gW//PKLbizitRythe12u4MsbFGxoHdBnDddHPm4NOUAOBjmVnXHwj1lEcewSU7TNlmZkS3sjIwMZGdn6xTE1NRUzdNA7nnRyKPzUEJmqHI2oGbtClJMVVnie/bswfz585XPUwWzrUl7A9gB4HMA9wP4AsAOxlhf02eqO9wGIBZALoBZAG7lgcVJDEHEJVvYKtAkTkhI0JJxxBiFal+gxsKWCVuerEYWJqCPYRslrYmrB9FLYzQWEsSyS5wsrJKSEsTFxQX9zumrwMXJ+9Ex1o25OS5sjz8e5QaCm8hMtrDFe0MvEpGR7BIH9IQtu9XEpg0iYmNjte5YspASobo/MmH//vvvOHTokK6MTCzrkuO14thE61klXMJZ2KKrUHaJi54dFVQxbLnhBWH37t0oLCxEamqqzkIQyVEkbCLRUO5PUtS2lcdhnqcbopkPo6K2oLmtLGhuhErEJAVNHLeYdEYJckTYJDjLy8sN45WhXOLicpKE0tJSuFyukIQN1Lh2Q5WziUlnKtewqKCJ77tRLb2Y5UxjUikcFOaRiYtImDGmc5EDAUswIyNDs7DF0BUdw0wMm+SkkUscqCF1ypOgfZo73DjLuwrNbaV4eWkevjkQDdj0JbJz587Fk08+qR1HtrDFbbLSTOOhedOiRQvk5OTo7l1qaqrOJU73RiRsmq8iYdO4w+U4qVziRvFyI5i1sD9AoP45k3N+MoBMAG8iEMduUHDOCzjnF3LO4znnx3HOZ4b7jcolHg6yhU3Zk4rrAaB3iRtp3rVJOhMhW9ikTISK28iELbqGyL0uw+v1IsZhwyUtj6BvXCG2oTX+cPSGlwdr/maSzkQLm4hA9hDQNZLlJlqz4gshggib7kNOTg7uvPNOneXg8/mUQrC8vFxH2BkZGTh06JD2vdg4ha6b3GFGFnao6gDVcxeFkRzDJmtKNW6gxq1mFMMmj4qIoqIiREdHIy0tTbteXr1SEQlxv9+vlXOJLnHVMyUlarOvJWZnudCMlWOkczOS7DW90EMRF+0jK2jy86b5Q3XRTqcTdrsdq1evxvnnnx90TLrXosCUFR4iHxElJSVISEjQCFtsYkQQ564RGYsucVlRoTAcKWj0/oh5J7Jh4PF4IorTqkIBRvkXADBy5Ei8+eabhi5xqiqgygEjr4vqWdLv5X1k4q+qqoLLCZzt3IbLeqdiZUE0Wox9FgXumnckNzcXO3fuBGDeJS7eGzFnpXnz5sjNzdV9T4QtykuZsEnOyxa2UdKZGDJVucTFWm0zMEvYXQC8xquPXP3vfwB0DvmrRoCqqipNoIo3X9TAfD4fDh8+rPtdXFwciouLdS5xo1gioHeDiW5P2QoPlXRGSTjhCJtiRHK8Sn7wpO2Kk0Is6ykuLkZsbKwyucLhcMDpsGNA1AGc7spFXnQGbpq1EW6uf1FFr0EoK0q0sEnoyIRNglCOjYqxSxEiYTudTuzbtw8LFy7UHdPj8SitKdnCTk9PR15eTatWMZ5N1y0qHuLYQj1TcR/ZJU6WlijAVUlnTqe+a1t+fr52rlBJZ7Lbt6SkBAsWLMAJJ5ygE0rl5eU6Tws9S5Gw5VIUh8MBd0UFPtlSic3R3dEz2YfhUVsRa6sRmFQ5IYZLVFUZsgCXFRUx3l5aWqqNu6KiQteVjyAqqrI1SgmQBQUFSElJ0f2utLRUR9iiEimOm7apLGx6lkahH/G5q1ziKsJu27atdhwzhK3yLDidThw5ciRoKU2SB1RpQeNWWdh0bbWNYQN6UhctbDq2jQG3DWyFS1oVISq9HV74y4ZZG0vxeUVvTN6ahhXpo7DLm6oR4MyZM3UWtpG1KiozXq8X6enpQRZ2SkoKBg8erCRsMeeJCFss2VK5xMXjkBJHz08MjdQHYc8FIKuyowH8YPpMxwi5ubmawCZBJMLpdGLnzp047bTTdNtjY2ODCFtl6YjWHE3y2pZ1ATXkYJRRW1paiuTkZM1V53Q68fLLL2Pu3LlBLjKRqEgxEc9TXFyMlJSUoHsivgAVFRU4Lb0KbbMXYGtOGeZWdUeRv0ZxEWOQNCGJKKi5iBzDpn3El54EocoKV2nOjDHtGZFwzsnJQWpqqrYPkahMSEANSZklbCIuOUuZxkSCwKisS+USlxPzVFniqlDMaaedhu3bt+tcw0CN5k7noTppQmlpKTIyMnTJaowxlJWV6e4PCSkxhi0nOfqZHSudvfH99jJklm7FhE5eOFiNUkj3Ro5hq2qNw3lUSBBTb2+am0YhnXCWpsfjwfDhw9G6dWvdd126dEGLFi1055fHLc5Lo5wNMRRgFN4Sk87CEXZ5ebkuFGDWwhbnlxFhi89DdInLLUtlD1gkMWx5H1XliSgP7HY7urvcKPnqMXA/x+ebS1CGaAAMlY54LPe2w7d/BhaGGTdunPY8SBZR4yfCvffeiwULFugSMuPi4oKMt6ioKHz33XdwOmt6QNC7IBo9FMKhZyuOW3yWojIt5oqIpE4KpFmYJWw7gE8ZY8sZY58xxpYD+AyAnTH2Ef2ZPmsDwmazaS4/IgoRDocD+fn5QS+JaL2FImyCqFWHc42Gc5+GsrBLSkqQnJysCS6Hw4ENGzZg37592n40AWWiOuuss7T96FgtW7YMslJEoqR4U1TOZnx8XX94uB1zq3ogx+/SjhPOeiaBK45NrkUUY8Cie4leRCIX8fnFxMToCDs3N1dH2PSyqSxsIlJyH6enp2suMsYYKitrltkjIaVy5YvjlpUQ8V7Qc6fxqyxN2eIyCgW4XC6tnajKiqR7KRMxuXxpHxJE8jwR20+ScBaJq4I78NpfHHnOVrhjUCtk5KxAlDO02zxUXobsRpXJkEiKLGya98XFxbrrpnuqKmck0DO48MILceqpp+q++/DDD9GxY0fNwqasXzFXJJyFTWMyCgWIFq1I2HRv6NpFiLkaKiuSxk0GA82Dyy+/HBs3bgQQIKOCggIdYas6pBFhFxcXIzExUEErKngimanGHY6wbTab0iUu7+Pz+eAsz4cNwWTmgx0vztumOyfNid27d2PAgAG6/bOysrBnzx4tFEiyZP369Yb3QSzzk8NLIuGL76nsWZC9q5S4Ku4TaYmXWcLeCOAZAPMAbK7+9xkEyql2CX+NDowxFBQUaC+22OEHCDycw4cPa80pCKKFTQ/E7FKCMhmLcW5xHyNSVmnu9NKoLGxyZcmQiSouLk6bpETYLVq0QHFxse53YtYlTcTy8nL0b9cM50ZtQTTzYF5VV+zxpWqkJBO2w+HQXEiiZRrKJS4mnRFhi64uuYOX7BLPycnRuTmpIYTKwiZBTs+CMuHp2B6PR1tmkK7JyCVOAkBWQkTFRHaJyyWFcpa4nHQmQqydVglwGlNBQQGaNWumbSeXLz3/srIyxMTEYOjQobpwD3miyNL6qzAKc9jJeLugBz6r6IOvK3viQClHp5xFuOrkzJAuUrLyw1lccgxbfM/o3oi9vWn+ynHo2267DaWlpWFdw8cdd5ySdOi6ibDFDHr6fTgLW0xANXING1nY8vkA6EIB9Lw3b94cdEw5V6N169Ya6RpZ2CLEGLZI2PR7sdwynEtcbDMrJ94a9XYQPTOUz1PiVz+j7EI3/JwhKSkJhw8f1uREfn5+EAE3a9YMhw4dClLievXqpTy2GBpQZe3LXiN6T0PNc5Gw5X4D8vMOBbOtSR8382f6rA0Iim2SwOacY+bMmbrvDx8+jISEBN2DiYuL0/rK2u12rZGCKpkHCN3pTHSbG9Vq0yQTY9gqK4wImwSXmPEtH0smqh49emgvOmnRKsIWXUz0cpHgSbRVYqRzC9JYGX7zdMTP2QxeiYxlNyJptbKFrSJsmdxEAS53qktOTtYs7JiYGBw8eFBnYScnJ6OwsDBIcVm1ahU6duyoOy9jDJMmTdLdQ3oBZZe4bGGLySRGFna4Gn3ROo6EsFUlQGRF5uTkID09Xdsu1uQ6nTWLOchlJaKFfdCZidl7HXCzWAAMbkShEk6M7hiF2MPbNVITn7dY+kSEbTbpjMYttxSlLHHRJS5b2IwxfPfdd8jKytISEVUWdm5ublD8WnzXRMKWw0pmLWyjcjY5IdJMDJtkEd2bwsJCnHLKKUHHlGunExMTMW3aNG3cMmHLcycqKkpTdjp16qQkbJWFLRItPUuVMkMJeaKFbeR1oaYy8VDLWwD4srIXXCech7wjRSHDJCJhi2EeI8hGA8kgsRqDGjyJctJo3HRMImy5o1+dEzZjbChjrH31/1syxqYzxj5gjLU0faYGBmNsNGNsSmVlJfbv36/r7jVhwgSUlZUBqCFsu92u09ZVLnEVYVO8NhKXuJigJDdwoeMYxbdFwnY6nVr7RFXiArmnCC6XC9dff712LLKwZZe4aNnSy6WzbG0+nBO1De1th/HVLj/+bj0EzFZDxjJhi9tkC5vc3HLylEiOosZN11FWVoaUlBTtPiQlJeHgwYM6gZSUlIR///vfWL16te4+bNu2TctZIDdgVFQUJk+eDKBGkBkRtkoZofsWysIORdhATU9yOflKfqHFkkMVSLj+9NNPQVYEjY2ev7z0IudcZxkeTO8Pj19WChgWZ3l17nJRSMnWoypbWrxHFH4Qxy17sqqqqpCUlKQp0SqXOACkpaXh4MGDOmvK6/XqPGoyYau6uJGFJSt7oieI5rIqLm9m3OQ1kl3iRoRNxKHyLDz33HP4/vvvDcN2Zixs0SW+YsWKoMYpKgtbtECJDGuTJS4eh5SZqKgo9PTvQrRDur/w4cbT2yOeVYGdeBmm5XfAZ5tKUAX1nBAVX6fTiaKiImXugwh6DmI5ndjvgEq46DvyLIh19yLE7HJRURGUrCTG2BTG2OhQ12XWJf42AHqDXgHgBMABKPudNgZwzudwzm9KTk7W1dfu3r0bY8aMwe7duwHUEDa1ayTQDRZdw0aEbbfbQ1rPBJWlRROCkhrEGAsdW7S+5Ri2y+XSFsAQsWLFCuTl5WnCho7x5ptv6o7VsmVLpYUthgJEwtY0TsYx2Lkbozo4UdK8D/49PwtubyAurLKwxbIuFRmLrirxZTGysEWPBykuZIXQWJOTk7F8+XLMnz9fdx+ysrK0hCOn04nCwsKgF5gxhpiYGE0bV5GT+LzEeyM+bzGGTYqaKjGNfi/Ge0m4qWKaFK4h3HPPPdr/RW9FixYtdGMS9yELW3xOIon4/Rxep65Nv4bDbn9IJUZMUDLrKjRyiQMBj0pSUlJICxsIdn96vV506tRJN26VhS3eG9HCpmsUfy+6xB2O4LatovKuimHTMcXkpVAWV0xMjNa4yGjcu3btwrZt2wwJOyoqKiKXuHy9KgtbNW5VL3ECGTbh1ieg9yU2NhbpZXvxwNDWiEclAI54VGKgYy8eGdUD50ZvReKq/yLFX4RP/irEs5viMT/XBWdS86BjijX2RUVFOoVHVVpFY83IyMDQoUN14RWS4aIMI0VF9ugQVPOcjlX9vIs45zdxzueon071cUJ9KSCTc76fMeYAMByBZTVvBTDQ5O+PGaKiopCVlaV93rRpEyZOnIj//jfQBt3pdCI/Px/NmzcPyi4koRPKJU6EG+olJahIXSauUJamHMO22+0aUYmIj4/Hww8/jC+//FJnHchZo2Zc4jJhi/2QGQOu7BGH5B1zsfpAKXKPH4sjbr+yVtnIJR6OsMUYtkzY9II4nU6NqEUFpmPHjpgwYQI2bdqkW0BD9DyQG018gek+id2YZJc4WZGi+1Z0iZP7T45hy8qc2FCB7i1p7kYucSo5pOfAOcerr76qe7bioiGq509zSc4E1nIFmnXCeW8sDTxkBdLi7GG9DmYIG6hJ4iEBLrvEgWDCJoVTJq6UlBTk5ubC6QyUcRYUFODIkSNafgIRdnJysvYb2SIiRVVFfrJLnNy38nwPZ2GTAKe2raS0i+VUhKioqCDCli3slJQU5OTkhLSw5aQz1T4yYcuJjKKFLSfIicQVyt2timHLVji1hnW73RjZPQ2Xx23E+Kg/cGnMenR0FGj7JlXlo9eRpZhycTv0SvFjoycd2zpfiWWedtiVV2PIFBUVITk5WbOwxfunSk6k93XIkCG4+eabdfuoCJvGJB9HrphQucRDJTPLMEvYxYyxFgDOALCZc053wrzz/RiBkpEIPXv2RK9evbRYJ7nW2rRpo0xMoliY2+1GQkKC7oFwzjXhKiaaGMWnRdewXN4kWm9yXa2KsGnBACJsUTN3uVxo3749Nm/eHNTxh0ACr3nz5ujevbvuOxI2dC0Ui6Nrkt3dMdlr8PKFneGNTcVNX+xAIY9TWtgqlzh9J74QYq9e8d7IhE0vCI27ffv2SEpK0lz8LpcLr7/+OnJzc3WWpgxVUhqgd4kbeQZEa9DIJS4KN3luiL+ncYttKKm8RrQ4qeMXjZusB9Ea2r9/P1q1amU4ZlX3OIfDgb+9cbjru72wDb0DJZUexB9agyhJSti5D+P7JetKvkTCkgnbiLjErG7RKhOVVgK5xEULWxVDTk5ORm5urhavPHDgAE455RRNaTeysEXIFrYIlYVt1EzHjPIuKjP0vGULmxYIoX1ULvHY2FgUFRVF7BIXlTixrEuGuNQqPXc5yU1OOjMat5kscSJsUlTExiMiKFzQtYUL13RjOKNkETIr9mK3rxmGvfIbFlV1RJ4vDkVFRUhMTITD4UBhYaHu/qlyHUTDRN5HLs8iBSM2Nta0hS26xOuDsN8AsArAJwh0PAOAQQC2mj7TMcTdd9+t/X/GjBm6SUsv4AMPPIDbb78dgN5FKVrYoqVJCOUSlxtXkAAXM4rlloOyAJcnBrVMJOJq1qwZevfurUs4SUhIQHl5Od5/v6YRnUo4l5SUoG3btvjhhx+084kkQdcrZmerJrLP58PgzmlwLn4DnAOrkgZjf2WNZU/ET5ambGHTMY0UFfE6aNxUsiQS9s6dO3WETbjllls0i4o6NYmQu16Ri0xcq5jCI+KckD0ioqVF7j+RsOn3MqnT78XnL1rY1CqUEBMTo8VyASAvLw9du3bFkSNHAAQE7759+3SELVry4pwgC3vb3yWYW9IGvzn7Y39hJUp/ex+/3DMESTt+wg194xDtLQPAkWCrwon+rTi7c7IuO14kTpHUyMI2Ssah5y2Xdfl8Ph2ZkIUtxrBlAe7z+dCsWTOdpXno0CHccccdeOONNwAEyO/QoUOmk862bNmi+15lYcshIJovRp3OaB9VtrERYZeXl2tywqhDIS3OQc9WjkHLhC3LJzGGLUNcgpOUB5I5ctJZuCU4VdnS8j7l5eU6wqY5IUNucOUtysEJfCfGRP2Jfw3phIP+RPzg6Ykhk7/Eyr1FmktcVNBV1QRyaELlEheNLXqWRo1bxBh2Q2SJPw9gGIBBnPNPqzdnA7jB9JmOIV555RXD78jtSS4poIa4RJe47BqmJAlRuBLxyO5ugjjp6DvRmlJZ2LJLnCxect+3bt0aM2bM0LlWWrdujdTUVFx33XXaucklRBCzkQn0AtLLTucWCVtlYWsWcd5efHjV8Yj2leGTrGTs8KYB0FuRNBaVpRpJDFtF2DabDUlJSUEu/ueff14bp+p72cKm5yIKcdnqVxF2uKQzUUFTzQnZwhY9OyJhM8Z0487NzcXxxx+veZJcLhd27typI2zKHBdRWVmJEn8U7vnsT4z4z2JkVcaiS8VWLPm/M7Hi4xcQ5QjkIwxpF4++B77Ec31LcV36XmR4DurGbcbCVlmaYjtOce1oet5ibLGyshKJiYlaKEjlNi8qKkKrVq1063sfOnQIHTp0QFpamnZvDh06pJW3qSASdrdu3XTfqSxsVTMdoyxxMYYtu8RpWzgLWxXDpvFTZrfqfZfd3XJWuarHP3kQExISUFxcrCVF0rWqLGwzvcTDlXWZIWyPx6MtDSp6H+Li4hADD+4b3hWXRv+FE+37UWxz4ar3V+KSKatwyNESscL9U3l0ZMIW9xFd4iT7ZS+hfP+MYtiA8UJUKpherYtzvp1zvkv6vMH0mRopSHMUQS+baGlSEoSYfCUKV5GgZcIWM8DFEiDZwpZdrCriIsKWu5OJL/qQIUO05DJCUVGRLqZplNhDL6ToEldZ2KLFIbq7M5Jj0fPQPHRw+bDM2x7rPBlaIhrdNzMucRKMRjHs4uJiJCUlBb1YiYmJynaVBLICRMgWtugRIEtMtH5FASISR7jGKaKCpnoOpAyJMWwiMxq36Eamcefl5ekIOyUlBZs3b9YRNnl+CHkllcjJHIyHl1Xihw2HcNPgDvhX2xy0LduKuCiHlqglE69IxnIIQxyT3HCG5rmqmkHlUaF5J7ohyeUoW9h0TCIpChEAAcIW7wOFkEQhSU10xOs3Wk7WjIUtu8RDZcfLLnGVYkMWtkxK8r0U3/HCwkLd+y4q+QSxmx8A3XtOEO/3Cy+8oCuDFC1sMaFMbg5C95jGbcYlriJsWYl2u93aClt0jaWlpUhNTdVyFqKYH8c7/saSB4bi2Yt7oaTCg4LuF+PBxWX4bNV+VHrVRKuysGkfUvBD5eHQmEWDTOUSjxSGv2KMbRH+f4Axtl/1V6uzNiKoBDzFI+jFV7mGZRenKIxFwpZfZDneKVvYctKZKMjIKiYrWITYFAMI1tpUhC1PUrEnuJFLXFVyJZdswePGDV2q0Nmej798mfhoO4PXrw8zyOOWk85UngWRsAsKCpCamhrkUlK5xEVcfvnl2oo/BHHxC7o3Xq8XQ4YM0Wr2jQjbTAxbtLDlntzhXOLy/RetUppbOTk56NWrl0bYDocDBQUFQRZWVVUViis8eGneNpzx4iKUtToBZ7WPxW/3D8VDI7sjIdqujOWJlpNKsZQ9SSKpGdWvizATAhHb/opKnKjM0Bz3eDwaERUWFurc36qqCpm4aLxGMWyxrEtUTsRnKrYmNSIluRsayRAZMmEXFxdrcx+oIeN7771Xc5XLFjZZoqHGrfLCiCGwe++9VzcX6ZmI76nsNSCISqhR4xQxXECepVAxbLfbjWbNmmkucZovKiMsxmnHlScfh1/uHQL7immIi3Lg/77cgDNe+BVrSpPg9qqzxAmivKR3IRxhAzUyoa4sbOPqceBG4f9XmT5iE4PKRUqETf+qXMNiRi9NKJmwVckoooUtCnTR3aoS4PRCcM6VFrYcs5Jx+PBhneBSuYHExi9inFa8bprIottWTijzer2IctgxOOYAYsvLseLv4xDT6UIUlXu0MYZyJxkJcJGwr7jiCk2TFuNyaWlpuOKKKwzvQ2Jioq4hBICgFp50b2w2mya0xWdhFMM2CoGIlpOKsOWcBXntY3Hc4hwjgeL1enHyySdj3Lhx2nHFFpUAYI+KwVa0xuAXFqGw3IPzerfCL6/djTvHv4WWSTXNVIwIm1ZqCqVYyvePhLOooMmWrep5q0IgZM3RuUlpjYuLQ3Z2NoBA61U5qeq3337TnVNF2HLfbFXSmVhrLFc8GFnYRjFssujJayC6xFVxWpGwKYbdvHlzLWchPj4eWVlZOq/akSNHdISdlJSEwsJC3XFVLnGVpSnnJ5AMJPmoykdRERfNH7kZkwyy3sO5xMvLy9GsWTPNwo6Li0N5ebnOmxBkNNkY3NuXY9rHr2KPOwbv/LoLC3YnwpE6HK/M34YJA9uhmUu9SptI2LGxsUEJw3IsXCzNDUXYkcCQsDnnS4X//1arozcBqLQxeiDR0dHK+DQQbGGHcokTQlnYZlziBJWFDYTW1ORMaVWiBVnYw4YNQ/PmzVFaWqrFEgkqK1CVUOZwOOB0OtCD78OpfbrgvXWZGPPucvhikoNc4kZJZ6Fc4pSXIAvp9PR0XccyM5DbMKrCBWZi2EYucfHeqFa9UsWwRWVGHLfb7UZ6errO8tm9ezfi4uLw73//Wzsu9Rnw+vz4fE0W3stuheL4TJzROhn3D++K4zOT0P6R/UpFRTVuWXlV5RyIK8GJXhdxvhCioqI0iy+US1xlzZFwpmsi9+eoUaNw3HHH4ZJLLtHOM3jwYN14XC5XkKUpr8RH4xUFPYFkgDg2o6SzcHXYojVqZGET0ckx7ObNm2uZ7/Hx8VojKEJ+fr6uy53K8yQrqioLWyZskdQplERyI5SlSe+NmCWuKoEiWVNZWamzsFWE7Xa7teZJdrtdu0+ivCorK0NGRkbQuF0uFwa3ScHgLul4/r+f4bP1BXh9YTSmLNmNK046Dm5brGHSmegSF99TlddS3KcuXOKGhM0Ye8LMATjnj9XqzI0EycnJmDNHX6tOgoCEnpFrWLaGSNDKmdCEUBa2GZc4EChLU8WjbrvtNnTp0sVwnFdeeaVu4qpeThIypKnv27dPmYxBloeYCQ1Ap7BQTMnr9eKMtrH44qOpyI27FvYR/4ft+ZU6MpYnu5j8FhMTo41XtZ65uKBFbSEm6wDGxCWGJMTnGyrpTBUmidQlLnZ4E+PHJEhVCUitWmXgh/WH8PL8bdidX4Z2LiB95xxMf+5dbR85dh/KwhZLHOncqnkqx7nJwpbLV6Kjo7UlbVVeG1LQRMIWLWyypkSFh5rHfP7550aPGikpKfjPf/6j20YhIAKd36jtq+gSF99bSkKl78LVYZuxsCk7nso4o6KiUFhYiPT0dM1TEBsbG0TYdrtdp6Cr3PvFxcW6d0elqMorvsmxXMrroTK6UK5hkpNy90QCLX5D41bFsMVnIsbLo6OjtTlBqwQC6nXPaYliQockGzoeWojZb32A937bjRkr9sF74i2Y8lcFmncuQZcWCTr5dNdddwXVYYdKOiNFhbLEjeL7ZhDKJd5G+H8MgDEIlHbtA3AcgJMBfGn6TI0UdrsdgwYN0m0jwiaCUxG2WHJVGwtbTDoTLWxRgMsF+gCwYcMGeDyeIOI688wzQ47z6aef1n2mvswi5KxPlVZL1yQnnQF6wiaNW7NID+/Gp7cOxFlPfIV/fbkD0a62ukx0VWMCEnqUKKQi7HvuuUdpCUUCWWDTGETIyVdmLGyRoFUWtthcRTyH7BIXY7lEbkblIJxzLNmRjxfmbcXG7GJ0aeHC1PH90dzzN9ZknKzbVxZmKkWFzvfJJ5+gY8eO2Lhxo6EnyOPxaN2k5Bi26GLknCvHJIdi5Bi22ONcFM6RxAAZY5gwYYJum9x+VGUpicqEkYVNigkRdlRUVFBYjGBkYcuJdHLb1ri4OJSWlmoZ80BwqSUAbNy4Mei9OHTokO5zenp6UPWBDNltLj53SpgTLexQLvFQhM051xQ+GndiYqJmjdK7JCedxcbGaiEGsrAvuugiDB8+3PDeUAtgAsn1Ts0T8OKlfXD32V0w9NansNI5COe8uhjDurdA//gamfD888/rxhROUZFj2LKiEglCucSvpf8zxj4FcCXn/Eth28UALq31mesZ1T1ZR4ttCc1ClTWsEpyycFVZU2Zd4qKmTvvITearx6V0Z0cKVRKKHMszImzR/SkqGHR9ssZN/3ZukYDybx/H8fd+gM29L8e768qwLjsBhW3G4rIZ2+BrfQKAGsWBFqugkhLR4iK8+OKLR3UfAOCOO+7QCSuVYiAStmhhm0k6E+cECRcicZEkQ1nYMrmpLId1+4/ghZ+24ffdh9E6JRavXNYHF/TNhN3GALRAnz59dPvL5WwqgUNzgnqSi89UtLDJVSo/d3Itixa6uOSmOCaa76qcETFu6HQ6ERcXh7KysoisEyOoEjJVHhbRilNZ2Kr4thwCIYRyidMiHECNokLziM4vxuLT09N1i4HQGGSILnIAePfdd8PWAJOiSRDfk6eeegrt27fHp59+qo1ffifE8leRsMU2nvQuiHKOCJvmvWpdalnRovBIUlKS9jwTExOD3hO5A6AcFs1IjkXUxu/w7r2XYlVhHD5cvhc/l3sR3eUS/LotF2d0SQ8aUzhFRSRsOZO8+p4mMcamAJgTqj2pWdNkJIBx0rZvAUwz+fsGR/Wg5/Tv3//GsDtLEMtCAOgEiFx6JGpQtOABlT1E4hIXNXXZwpZjZJFYFEZISEgIcqOpLK5QFra4dKZ4bUTqooWtkVJpAf479ngMmvQFFu1tA8AGMCCnxAPHqePxzbpsTXCT5Wuz2dCvXz9deVhd4vXXX9d9pvioCDEmG87CFpUZ2cIWyZisJqMYtph8RcQVFxeneRzIjb89pwQvzduG+Ztz0Cw+CpNH98CVpxyHaEfoxBa6pwQxHkwgtydBtDBFjwgRLhE2KShkEctWmMrCJkVNnDd0b2Th3Lx5c8ydOzdIyawNYmJigsYoz3tSRtLS0nDddddpc0G8VhVhy0ocEQ7V0otuc1LyRUVFbMYUFRWlPSORsEeMGIERI0bUatwyVq9erftsVIoHIGjFMLoXonzyeAJd7FSlX/TekPIjyjkiN4rpkzVPiVyAMWGLePbZZ4OUkl279CtBk0Egwu12o2VKAu7q3xE3De6AZz79DbPWJeGaaavQvVUibh3SER5vsEtclk+qsi56tiTPqu9XEef8JuWNFmCWsHcC+BcAUbLdhka6BvbRQo6LqpJgZAtbLOMSLS45GYUmm1HSmUhuJMDk8ou6QGxsLP744w/dNrm8ySgJRY5hG8W5VRq3z+dDQlw0uFOv9QIAc0bj+Z+2BrnEmzVrhi+++AKcc63zVH1C1Q0tlIUtx7CJZGXCFsu6VBa6qLCJyXYiub3wwgtwOBy44YYbUOxz4oXZf+HrdVmIj3Lg3rO74LrT2iM+2txrvXbtWt3nlJQULfOYoBKKZDXTXBYtZHrudE/I7S9a2Eb5IAB090YVryQwxnDSSSeZGmc4fPrpp0ENQIwSMpOTk/Hvf/8bpaWlmrCl+S6GQOj65TCJ3IxH5RI3CgXQveWcK625usCJJ56o+3zjjTeGbHFLYyJZICfNVVZWIjU1VedZkkNvYuzeaKlQWohHtLDlUlan06lbBId+J0NcIx4AWrZsqauwAPTJrHFRDpzXJR7z334d9785G+/+tgt3zloHZ+crsSLPgdN4jaynZyI31lK5xFVlnuFglrBvAPA1Y+wBBDqcZQLwArjY9JmaEMgdQ1AlFYiJNmIdtvjQIinrUiWd0cuqImz5xYoUjDH07NlTt022sF0ul2aFFxQUaNckWtgqtzk1cVFa2NWfeVwKVH6CQ0UVSLxoEqYuz0KzLifCGVUzbsZYkFuvPkCrI4kQrV/RwhbHJj5TUQCLgluOT8v3RnQfq1ziTqcT+aWVOH/yRxj11kqAAdef1h63DumE1HjzPYlVyMjI0JZfJRDhEsglK0KOq6vmBB1HFR9WlYXJsVxZcahLyEJdVKwJctKcOEb6jshEJC7Zwgb05Wzi2gOiS1xOwJRJLiMjA0899VQd3QFjDB06NGgbdY0jyEsEixBd+nQf5OdtJEPE8YoWNt0/uboDAF566aUIRxjoDCkTvVx94nQ64a2qxKX922DMCa2xYEsO7nzne8zckYT5b6xETL/RqPTbdLKfPARGLnFVr45wMLUn53wdY6wzgAEAMgAcAvA757zuVbxGgI0bN+pqGIlExRdJJmy6+XLHLvFhqJLOQlnYJAhUDeJl11VdoG/fvrpabTExjbaTNSG6xI0sbHqBRcKhe8LLDoO59C8+ALiiHai0O/GfhTvBznkA/Z/+BYM6peH0zuk4vXMa2qSGXse2LqBabo8UjV27diEuLi5kDFsmbPE72cJWucRVhF1VVYWSCg+mLtmD95fshtsTh8v6Z+LOszojI7luiKx58+aYOHGibtuQIUN0n8V34NFHHwUAHUmpcjdE0JjEZDGVABctdKBG8G/atOnoBxoGRslXRkoFkQu5j8XfqzrbUYWHmLMgeulULnFZBjRr1gwXXHBB3Q7cJA4ePKj7HCqPQB4joKiNNrDMw1nYVN2xdWvdL2nRs2fPoOz4GsuYYXjPlmi2bhquuf9pLM6PQf7p43HFzB3wHX8ecksqlGW/MmGrxh0Opqm9mpyXRHT0JgqRrIEaqyImJgadO3cGUKP50r9iDFu0mIw0SVG4yxa2rLnLL3194bPPPlOOW4SYdCaPiUCETWNLTk7WmjbQPXGv+AzNRtyBCm/NhHUyjqcuPB4XnTACh0srsWxnPpbsyMOSHfn4cePfAID2afE4vXOAwAd0SEVCTN0vGNelSxfcdJM+nETlLTQ3jFziKgtb/k7OeQhH2NzmQGW7QRj8wiIcKfdgVK9WuOecLuiYbrxUYl1h0aJFQdvIIn3iCX3lJ8XwQrn5iNDkmltRONO4RTf4t99+C8YYevToUSfjCofevXvrPquUb1mxFuOrBJnAVfJBnC9kYRM5UIxeDEX98ssvdT/gCCDHhK+77jpkZmYCAL78Ul84RM9SnhPyPZHrveU5IdZYk+z5+eefERsbi65du9bd4KqxatUq3WdVkq6nqgp9M+Jw4wUnIaZVZ1z9zEdYUDkcpz2/CJm8AwbbPDqPGpV1yS7xSHB09TD/IyBL0+VyYfny5QD0ySAkXFWEHcrCFl2kYuyS9qHGBpHGOeoKNptNK2EgkAARE21kwq6oqNCVo8lNGxhjuGXkiTjhol544ts/caQCyEyJw/3Du+LCfplYtmwZUuOjMLpPBkb3yQDnHLvySrFkRz6W7MjH56uz8NHv++CwMZxwXEqAwLuko1dmUnVG9NGhTZs2uPLKK3XbVAuEhIphy4Qtrp0tli6pXOKaFwIMsT3PwtWf7gA74RIcn5mE+4d3Re/WyUc9xtpiwIABQX0LTjvtNAA15TFer9dwvtI9E4mM3hsCkdOYMWMwatQoAMD5559fH8MxxJ9//qn7LJOxOD6RsFVQrdgHBPeOV1nYFAoQLfxwJZwNjTPOOEP7/8UXB6Kkl112GYAaBSWUDKOYPdWyA8EhELKwRUtdLsetT6jyO8TnXpWzCy9c2A2d+w3E2Fc/x2crvfjoMENakg09q5xaMq7sErcIux6g6iQk9jYWE07kjl1GhC0mI4klEaIAX758OTIyMoIScBoSDzzwgO7ziy++iJYtW6Jt27bo2rUr7HY77rjjDt0+/fv3R2pqqk67lJPFqJ833/sHPv74Y3z99dfadwMHDtTtyxhDp+YJ6NQ8AdcOao9Krw9r9xVq1vcrP2/Hywu2IynWidM6peG0zmk4vXMaWqfUnftcdokaES49fzmWS89b/k4lwD0eD5YfcOOW75Yg7dyJSHdFY9fMx/HxpmV1Np7awmazBZXJLFkScLxRkk2/fv2C4qsXXnghgJrQwoABA3Dcccfp9tm/P7A0ASX2RUVFRbRWcF1Cft9UzYoIcgw7FMJZ2DJhi1nwx0Jpry3IW0fj7d27N9q0aWO4v3y/5RAANU25+uqrj7rvQm2QkJCgW/0Q0MuEDz74IFANUZCNZy/uhUPzp8J93KlY5c/E27uc2DztDyC9k9Y7HtAn25mFRdgmEBcXF5QpKceWVBaXaGGLCUniZ3Kpk+AXX8q2bdsCACZOnBjkojtWGDNmDADgvvvu07bJy5d+9913APQJN0ZuKzl+aQbRDjtO7dgMp3ZshgdGIOA+33UYS7YHCPyHDYEGER1E93nHZnCZzJ42gspSkhuIiGV4quctNweRy/i2FPjBz74fk+YfQKfmLuR+9TRmrZmHVo/8dVTX3hA4//zzkZmZiZYtW6Jly5a670ghozF37txZCy8RSKCnpaVh9OjRDXPRJqHK1SCQ8q5yiYugWL+s0It5DVRnTjkfcqexpgaXyxW01C9BVMwp214uX+3Xrx+AQJjq66+/Rvv27ev/ohWw2Wx47bXXdNtEpeLaa6+F2+2uyRJ3F2FEqwq4136LzDOuwB9ZRUi48DFcP2szYjuejK/XZOGDvPYoie6M+Ioq2GITU81cR8QSjDF2NoDjAazknC+P9PdNETabDTt37tRtI6IlwqGEEzELVrXWtRzTlJuMqEBupqaGUMlHhBNOOCEoZyBSNHNF4/w+GTi/2n2+M5fc53mYvToL08l93jYFg6sJ/PgI3edygo8R4ao6dgF6C1t2ift8Pvx1oBDlp9yAR37OAWIS8Njw9phwRne8Gz1OW9+6scNMLbBR5y8RnTp10vVFbww477zz0Lx5c+V3YoJRqGQ7WmqTrHDZwlaVhlLzIOrc1dQwadIkQ+/gsmUBjxHdj7vvvjso65tKD2NiYkK2Xj4WkN32YsKxlh3vLsG57ex4dcSZaHnqRci/ZCKaX/IY7vn8L3A4AQaUIRqOxPS2Zs4ZkrAZY7MA/MI5/2/15/8D8ASA9QCeYozdwjn/uHbDbVqQJx1p1aQdUxxGdImHygg2ahbRlEFuT4IZwlZZWkcDxhg6t0hA5xYJuO60gPt8zb4jGoG/NH87Xpq/HclxTgzqlIbBndNwWud0ZIbJtP7mm290n8V4tFg5QBa27BKnOlT5uxw3Q16XC3HBW8uAxFa4Y1BLPDlhAi56cgfsNoZ//etfABBksTYVyEliDocjKLxUm6YfDY1Q13j11VejRYsWIUt0KF4pusyJ6GluiGWAYq1xfHw8fvrpp7odUANB5cYnzyGB3N0dOnTQtokhssaKxx57TFdmKmZ9i53PoqOjEeO0o2jNHCxa8TU63PMZ7HH6bmtgzFS8I5yFPQjAxMDxmA3AfQDGcs6/ZIyNBPAcgEZJ2EfTmtQMxAxuxpiuoYRI2FSfKPeNFglbdIkDQGqqKe9Io4P8kjkcjiCr5I033mjIS0K0w46BHdMwsGMa/m9ENxwurcTSnfkagf+wvtp9nh6PwdWlYwM6NAvbfISetxjbJAtb5RIHAsoEfVfid2JLwon4bVExkNIedw3rjOevPQdj/+93PFbpjjhM0FixceNG3ec+ffoEKaY//vhjQ15SnWPy5MkAQnsPKL4vVnyIVSRiiaRI2F9//bVuEY9/AmhRJcJll10WVI0iK/+NEbfccovus2jU0TOU3fxOhx22WPViRXXRmjSZc55b/f9+CCwC8k31558AzArz+2OGo2lNagaXXnop0tPTNVenirDlkg5xHWsS7vR70Q1G6/s2dcg9ugHg9ttvP0ZXE0AzVzQu6JuJC/pmgnOOHbmlWLw9D0t35uPTVfvx4fK9cNoD2eeDuwQI/PiMJNgk9zkt/edwOHDuuecCCB3DJnjtMXjttwOYtTMZ/rhEjOoch2VTn8NdLy3Bk+UlWg/1phy3FCE/f7pX/wScffbZus+hPEpkRcsNRERhrkpsGzBgQB1f9bGHbHUbhRqaMoiwVYvT+Irz4UgKHnNdtCbNZ4y145zvBTAUgWYpJIHiAUSW4vYPwgsvvAAgEKPp3r279rLJTRPEMh+xEQm9rNTBh0o3AHU7vaaIxp7VyhhDlxYJ6NIiATec3gEVHh/W7juCxdXW94vztuHFeduQornP03Fa5zRkJMeiZ8+eWLo0sGQ8lTkZJZYBAIuKxWs/b0fiuFcxc/VB9Ev2wL9+Dm648k4sfSNQ30lem5dffrlB6u4tHB3mz5+v+xyKsOVGSUBwwiXJgvpoRWyhYbB9+3YANWWAqnK2I79Nx3FjHoTbI9An56bqu8IR9n8B/MAYmwdgPACxfmcwgPpv7NzIESqblV5gMZZNri5ykY8fPx4tW7aEy+XSak4tHBvEOO0Y2CkNAzul4cGR3ZBf3bxl8fYAgX9f7T7v1NxVnX2ehlNi4nXu82/WZWNDu8tw76poNN98AJUZ/fDB0j3IvPm/eO3nHajc9xd++s892Lj8Z3y+uiSoAYfT6cTdd999TMZv4egQKoYtKu8kwI26eZ122mno2LFjg1yzhboF5eNcf/316NGjB956660gwi7f8huevfgTPP71WhypZEh0+HCgOG+fmeOHJGzO+TOMsWwA/QFM5JyLLvB0AC9HMph/OvLy8nSfqdmGkebNGMNtt92mfX711Vfr/RotmEea5D7fnlOKJTvysHhHPmb9sR/TlgXc5ye2TcHpndPR45J78NBX6+GJCsSocsu8QP+xeOL7zajK3Yt5L96GE9qORtfZj2OzYG3Jne0sNE1cfPHFhnXYlHwkWtjy86awWnJy8lFXTlg4tjjnnHMAAP/+97+1pWlFXNgvE56dy/HZZ59hxIgRuNFdXGDmuGHLujjn0wFMN9huQYDcFP+9995D165dtX7QAHD48OFjcWkWjhKMMXRtmYCuLWvc52v2HcHiHXlYsj0fL87bBnQ8E/D45R+iWXwU1n72CPp++rAWyySLq127dnjnnXcAANOmTWv0YQQLxpAX0xHh8XjgcrmCei3YbDYcf/zxAAJJmqGai1hoelBVFwwbNgxATQht5MiRpo8XlrAZY+0ATAZwNoA0APkAfgEwmXO+2/iXFmi9WNEN1lQzwC3oEeO0Y1CnNAzqlIaHRgJ5JZU46emflfsWlNUsUXrrrbfq1kqPi4vTFti45pprGuDKLTQEpkyZAgDIzQ3k7FISUo8ePYIyotevXw8AaNeuHRYuXNiwF2qhwbFgwQIANVVE1IfdDMLVYXcHsAzACgCPILBKVysAlwFYzRgbxDn/n49jyzj55JN1n8UVjiz8M5GeEI3M5FhkF7qDvstIjsXe6v+//fbbAGpeVgv/TNx4Y6A4hep0qSZX1QRJnAeWh+Wfi7/+0ncrHDJkCCItOw43O54D8Bbn/FzO+Yec83nV/54L4G0AL0R0NgMwxqIZY+8zxvYxxkoYY+uq67zFfc5ijG1ljJUzxhYxxkx1hjkWWLlype6zqhe5vP6qhaaP+4d3RaxTn6sQ67Tj/uHBbVlDdbaz8M+DpbRbkNtLt2rVCieddFJExwjnEh8MYILBdy8D2BPR2UJfxwEAZwDYD+BcALMZY70453sZY2kAvgJwA4A5AJ4E8BkC63M3enTv3l230hMAvPyyla/3T8OF/QKurRfnbcPBQjcykmO1FchkJCYmBnV8svDPRYcOHYKWZ/zqq6+O0dVYaKpgoVaXYYwVA2jFOS9TfOcCcIhzrm7bcrQXxth6AI9Xd1W7CcA1nPOB1d/FIxBL78c5D7l6ef/+/fnq1avr4xItWDCNK664Ap9++qlum6qpggULFv73wBhbwznvH26/cC7xVQCuNfjuGgD1woSMsRYAugDYVL2pJwAtAFCtQOyq3m7BQqOHTNZAcBcwCxYsWAiFcC7xRwHMY4x1BfAFapLOLkXAVV7nS8gwxpwAPgEwXbCeXQDypF2LACit+2qL/CYAQWvuWrBgwYIFC40MaYwx0QCewjmfIu8UrnHKcsbYOQCeB3ArAha5H8DvAEaYXV6TMfYrAvFpFZZxzk+r3s+GwGIiVQDEptOlABKl3yUCKIEC1QOdAgRc4mau0YIFCxYsWDhGyDfjEjfTOOV3AIMZY7EAUgEc4ZyXM8bsjLEnOOePmTjGkHD7sIB/8H0ALQCcyzkXUyo3QUh+q45hd0SNy9yCBQsWLFj4R8N00R/n3M05z+acU9W/A4Ha7LrCOwC6AxjNOZeLWb8GcDxjbAxjLAbAYwDWh0s4s2DBggULFv4pONoq/TrJmqmuqb4ZQF8AfzPGSqv/xgEA5zwPwBgATwM4AuAUAFfUxbktWLBgwYKFpoCQZV0hf8hYNIByznmjXgeQMVYEYIe0OQmBpDWjz3W1T30dt7Gfm1rYHotzm9nHOnfjOLc8T/5Xxm2d23ifupgTTXHcnTnnSQgHWmRb9QfgzBB/IwD4Qv2+MfwhkG0Xclt97fM/fO7Vx+rcTeDeWOc2mCf/Q+O2zl2Pc6KJjjtom+ovXNLZ+2G+3x/m+8aAOSa21dc+/6vnVqExXZ917sZx7nDHqM9z/6/e88Z+7nC/qc9zN6Z7rkStXeIWLBiBMbaamyhRsPC/DWueWJBhzYnQsJaGsVAfCCr4t2BBAWueWJBhzYkQsCxsCxYsWLBgoQnAsrAtWLBgwYKFJgCLsC1YsGDBgoUmAIuw/wfBGItmjL3PGNvHGCthjK1jjI0Uvj+LMbaVMVbOGFtU3diGvrufMbax+nd7GGP3S8duV/2b8upjDAtzLU8yxjYwxryMscmK78dWX2cZY+wbxlhqHdwCCybQVOYJY6wVY+w7xthBxhhnjLWrmztgQUYTmhOjGGNLGWOFjLG/GWNTGWP1shR0Q8Ii7P9NOAAcQGBBliQEVmWbXf3CpAH4qnpbKgJLqH4m/JYBGA8gBYFa/NsZY2LXuVkA1gFohkDr2i8YY+khrmUngAcA/CB/wRjrCeA9AFcj0GO+HMDbkQ7WQq3RJOYJAgsS/YRAN0QL9YumMieSADwFIAOBltetAbwY0UgbIaykMwsAAMbYegCPI/CyXMM5H1i9PR6BzkP9uKJ3O2PsdQTm0R2MsS4ANgBI45yXVH+/BMAnnPN3w5x/BoCdnPPJwrZnALTjnI+t/twRwBYAzej4FhoWjXGeCN85AHgAtOec7639KC1EgsY8J4R9LgbwOOe8V23G2FhgWdgWwBhrAaALAquf9QTwF33HOS8DsKt6u/w7BuB01Kya1hPAbolM/1L91iTka9mFwNKrXWp5PAtHgUY8TywcIzShOTEY/4DVHcMur2nhnw3GmBPAJwCmc863MsZcAPKk3YoAqOI/kxFQ+qZVf3YhuEduEYDMWl6e0fGafCyqqaGRzxMLxwBNZU4wxs5GYHnmU472WMcaloX9PwzGmA3AxwhYrbdXby4FkCjtmghA54JmjN2OQDxqFOe80sxvGWObWM1KbKebuERT12KhftEE5omFBkZTmROMsQEAZgK4hHO+3ezvGisswv4fRbVL6n0EkrnGcM491V9tAtBH2C8eQEcI7iTG2HUAHgRwFuc8SzjsJgAdpGzMPvRbznlPzrmr+m+JicuUr6UDgGgATf7FaypoIvPEQgOiqcwJxlg/AN8BuI5z/kuEw2yUsAj7fxfvIJA9OZpz7ha2fw3geMbYGMZYDIDHAKynpBEWWKP8GQBnc853iwes1mD/BDCJMRbDGLsIQG8AXxpdBGPMWX0eGwBH9e/s1V9/AmA0Y+z06pf/CQBfWQlnDYqmME9Q/V109cfo6s8W6geNfk4wxo5HoHLgDs65qYU1mgTMLOll/f2z/gC0BcABVCDgiqK/cdXfDwOwFYAbwK8IZGrTb/cgkIkr/u5d4ft21b9xA9gGYFiYa/mw+lrEv2uE78cisCpcGYBvAaQe6/v3v/LXxOaJ/B0/1vfvn/jXVOYEArFxv3SuTcf6/h3tn1XWZcGCBQsWLDQBWC5xCxYsWLBgoQnAImwLFixYsGChCcAibAsWLFiwYKEJwCJsCxYsWLBgoQnAImwLFixYsGChCcAibAsWLFiwYKEJwCJsCxYsWLBgoQnAImwLFixYsGChCcAibAsWLFiwYKEJ4P8BM3xXmHqNrRQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot\n", "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=[7, 2.5])\n", "ax.plot(dt_out, tide_u*100, lw=0.5, c='k')\n", "ax.plot(dt_sen12, tide_sen12*100, 'o-', c='C0', label='Sampled every 12 days (by SAR)')\n", "#ax.plot(dt_sen06, tide_sen06*100, 'o-', c='C1', mfc='none', label='Sampled every 6 days')\n", "\n", "# axis format\n", "ax.tick_params(which='both', direction='out', bottom=True, top=False, left=True, right=True)\n", "ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))\n", "ax.xaxis.set_major_locator(mdates.MonthLocator())\n", "ax.xaxis.set_minor_locator(mdates.DayLocator())\n", "ax.yaxis.set_minor_locator(ticker.AutoMinorLocator())\n", "ax.set_ylabel('LOS displacement [cm]')\n", "ax.set_xlim(dt.datetime(2020,9,10), dt.datetime(2020,12,15))\n", "#ax.set_ylim(-22, 40)\n", "ax.legend(loc='upper right', ncol=2)\n", "fig.tight_layout()\n", "\n", "# output\n", "out_fig = os.path.abspath('SET_TS.pdf')\n", "print('save figure to file:', out_fig)\n", "fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=600)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fig:SET - Plot power spectrum in semi-diurnal and diurnal periods" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "save figure to file: /Users/yunjunz/Papers/2021_Geolocation/figs_src/SET_PSD.pdf\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAC+CAYAAADdq82HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwIklEQVR4nO3dfXwU1dnw8d8VkqAkgICABBBCokVFDBAj3iUY0dZWi3i39bVFaSvw1Cq+gOhdraK1Vp8ACrdWRaFURaqtVUHxuVUUkeptDBIF3yoRUIlGRAwCBVGu54+ZDZtks5nd7GR2N9f385nPbs7MnLnm5ezJvJ0jqooxxhhj0lNG0AEYY4wxxj9W0RtjjDFpzCp6Y4wxJo1ZRW+MMcakMavojTHGmDRmFb0xxhiTxqyiN8YYY9KYVfQ+E5FJQceQaOm2Tum2PpCa65SKMfvFtsV+ti32i3dbWEUfBxEZG8PknnZMLHn6MW2qrJNP28lz4Wnv6+SzVv+gJ2JdkiSPtNgWCTq2bFvsZxV9G/LjhzGWPP2YNlXWya/t5Eee6bhOyS4R65IsebRWMqxHMmwHaOfbQrw0gSsiQ1X1zTaIJ1AHH3ywDhw4sMXp6urq6Nq1q6c8t2zZQs+ePROapx/Tpso6+ZGn1/Xxa/mpsk6rV6/+XFW9ZRqH3NxcHTx4cKvyiGVbJnMesew/v2JIRB6JiMG2xX6rV6/eoaqdY50v0+N0y0WkBngAWKSqn8S6oFQwcOBAKisrgw7DmKQkIpv8zH/w4MFW/oyJQkTei2c+r5fu+wDXAccB74vIMyLycxHpFM9CjTHGGNM2PFX0qvqNqj6hqmcCfYFHgOlArYjcLyLf9TNIY4wxxsQnpofxRCQXOAM4B+gH/BV4H1gkIncmPDpj2sC2bdsQEQYMGNAg/eOPP6ZTp0507949oMiMMab1PN2jF5HTgPHAD4F/AvcBj6vqbnf8ncCHwG98itMY31RVVdG3b1+2bdvG9u3b6dKlCwC//e1v6devH3379g04QmOMiZ/Xh/FuAf4CXB7pQTxV/UJELktkYMa0laqqKoYNG8bnn3/O22+/zciRI3n99dd5+eWXKS0tpVu3bkGHaIwxcfN66f73qjqzcSUvIj8NfVfV+xIamTFtZM2aNRQVFTF06FDWrVsHwLRp0/jjH//I22+/TVFRUbABGmNMK3it6JurxOclKpBYiEhHEZkvIptE5CsRWSMiPwwbf5KIvCsiu0TkBREZEC0/075VVVXVV/RvvfUWS5YsYffu3fz4xz9m7dq1DBs2jFdeeYXjjz+eE044gXPPPZe9e/cGHbYxxngStaIXkUEiMgjIEJH80N/ucDKwu23CbCIT+Ag4AegK/A54REQGisjBwD/ctO5AJfBwQHGaJLdnzx7eeeed+oq+qqqKq6++mtmzZ/Pee++xb98+jjjiCAYMGMDzzz/Piy++yKBBg3jiiSeCDt0YYzxp6R79ekABAaobjfsUmOFDTC1S1Z2Nlv2kiGwARgA9gLdU9W8AIjID+FxEBqvqu20dq0lu69at48ADD2TQoEH06NGDFStWcPbZZzNy5EgWLVrEkCFDyMzMJC8vr36ezMxMMjKs9WhjTGqIWtGragaAiLyoqie0TUixE5HewOHAW8CvgTdC41R1p4hUA0cBVtGbBtasWcMxxxyDiHDQQQfxwgsvMGTIEGD/Jf1wGzZs4Omnn+aaa64JIFpjjImd1wZzkrmSzwIWAX9xz9hzgbpGk9UBEdsHFpFJIlIpIpVbtmzxN1iTdBpX5mVlZRx88MHA/of0QrZv384FF1zAAw88QHZ2dhtHmhQODpUVd0hEr2JW/ozxLq4y2GynNiLy/1T1B+73l3Au4TehqqPjjbi1RCQDeAjoAoxT1b0iMgfIUtWLwqZbC8xQ1Uej5VdcXKzW1nZ6qa6upqamhry8PAoKCuLO55tvvmHcuHFMnTqVMWPGJDDC1CEiq1W12K/8rfwZE128ZTDaGf39Yd/vA+Y3MwRCRMRdfm/gJ6oaegz6LeCYsOlygAI33bQTFRUVlJaWUlhYyOjRoyksLKS0tJSKioq48lu8eDGvvvoqN954I2VlZTz8sD3faYxJDZ66qU1GInI3UAScrKo7wtJ74jxE+EvgKeAG4ARVHdlSnnZGkR4qKiooKyujR48eXHrppRQVFVFVVcWcOXPYunUrK1asoKSkJOgwU46d0RsTrHjLoNf+6M8FqlT1HRH5Ds77898CFwXxJLv7XvxGYA/wTdioyaq6yH317w5gAPAqMEFVN7aUr/3QpIfS0lI2btxIZWUlvXv3rk+vra2luLiY/Px8Vq5cGWCEqckqemOC5cel+3A3AV+432cCrwErgT/FusBEUNVNqiqqeoCq5oYNi9zxz6nqYFU9UFXLvFTyJj1UV1ezatUqLr300gaVPEDv3r2ZMmUKL730EtXVjd8WNcaY9OS1ou+pqrUicgAwCrgGuBHn0rkxSaOmpgagwdPy8+btb8AxlB6azhhj0p3Xin6LiBTi9F73mqruAQ7AaUjHmKQRatimqqqqPm3SpP1voITSwxvAMcaYdOa5UxtgNc5T7uVu2kmENUxjTDIoKChg1KhRzJkzh9raWgBmzJgBOPfo586dS2lpaatetTPGmFTiqZtaVV0oIo+433e5ya8C5/gVmDHxmjVrFmVlZRQXFzNlyhRuuOEGcnJymDt3Llu3buXRR6M2p2CMMWnFc4PdbgWfG9bRTS7QybfIjIlTSUkJK1asID8/n+nTpwMwffp08vPz7dU6Y0y74+mMXkR+gHPZvk+jUQp0SHRQxrRWSUkJK1eupLq6msLCQtavX2+X640x7ZLXM/o7ce7T56hqRtiQFpW8iIwVkXl1dY2byDeprqCggMrKSqvkE6OriMwTkbGJzNTKnzGexVUGvVb03YB7VPXfsceV/FR1qapO6tq1a9ChGJPM6lR1kqouTWSmVv6M8SyuMui1op8P/CL2mIwJXnGxb425GWNM0vN0jx4YCUwRkauBT8NHBNl7nTHGGGOi81rR3+cOxhhjjEkhXt+j/4vfgRjjl+uvvz7oEIwxJjCe7tGLY6KIPC8ib7ppo0XkLH/DM6b1Qi3jGWNMe+T10v2NwPeA24G73bSPgduAR6LNKCJe+wPdrarf9zitMZ7l5eVZJzbGmHbLa0U/ARimqp+LyF1u2gZgkId5jwX+TwvTCDDHYyzGxOSTTz4JOgRjjAmM14q+A7DD/a7uZ25YWjQve7nHLyLneYzFGGOMMR55fY9+GTBbRDqCc88ep6W8Fl/aV9WTvCzALtsbvwwfPjzoEIwxJjBeK/orgDygDuiKcyY/ALjKp7iMSZjVq1cHHYIxxgTGU0WvqttV9Qycyn0kUKCq/6mqX8WyMPep/QsjpD8VSz6JZm1tp7dJkyYFHUK6sLbujQlWXGVQVDXyCBGv/wTs87wwkd1ANfACcKmqfuumb1fVLl7z8UtxcbFWVlYGHYZJMBGhuePceCciq1XVt/aErfwZE128ZTBaZf4NsNfDEIuvca4IDASeE5HubrrEmI8xxhhjPIhW0efjvD43CLgEeBH4AXCE+/kCcHGsC3Qv948F/heoFJGh7H+S3xhjjDEJ1Ozrdaq6KfRdRK4AilX1SzfpXyJSCVQCd0WYvTni5q3Af4nIG8BzwAExxm2MZ5s3bw46BGOMCYzXp+67Ap0apXVy02Pxq/A/VPWvwCnAH2LMxxjP7Kl7Y0x75rXBnL/g3FO/HfgI6A9McdM9U9UmzeWq6hpgTSz5GBOL008/3R7GM8a0W14r+unAeuBsnPfpPwHuAO71MrOIvEQL9+GtX3tjjDEm8bx2U7sPpzObu1uathnhfdkLcCdwUZx5GWOMMcYjr2f0rdK4rXsRmW193Ju2cs899wQdgjHGBMbrw3jGpCxrGc8Y055ZRW/SntMHkzHGtE9tculeRMY0Xq6InEhYi3iq+nxbxGKMMca0J81W9CJyo5cMVPU6D5PNb/T3VmBBeDY4LfAFwu0gYGxhYWFQIRiTCrqKyDxgqaq22EW1V1b+jPEsrjIY7Yy+f+tjcqhqfqLy8oO7wZYWFxdPDDoWk3g/+tGPgg4hXdSpasIfeLDyZ4xncZXBaE3g/qJ18RiTHJYuTdjJpzHGpJyYHsYTkc4iki8ig0KDX4EZkyhjxya0+3RjjEkpnh7GE5EjgUXAMTj304X9Ld118Cc0YxLjySefDDoEY4wJjNcz+j/hdEvbHdgOdAPuAS7wMrOIfCQi80TkDBHJiStSY4wxxsTMa0V/DHCV202tqGodcCXwe4/zlwCvAuOBjSLyrIhcLiKHxxqwMcYYY7zz+h79biAL2At8LiKHAtuAHl5mVtVPcF6xmy8imcBo4FTgcRHJBpa5wwuquie2VTAmOuu5zhjTnnk9o38JOMv9/nfgaeBFIOZGblT1G1V9XlWnqeqRwMnAe8Al7mBMQs2bNy/oEIwxJjAS69mOiGQA5wGdgftVdacfgQWhuLhYKysrgw7DJJiI2Fl9AojIalUt9it/K3/GRBdvGWzxjF5EOojIChHpCE6Xtar6oKreFa2Sd+f7sfsAXmZY+pmxBmmMMcaY+LRY0avqt0C+l2kbuR8YDhQBq0Qk1L7lr2PMxxhjjDFx8lp53wDcJSID3DP1jNAQZZ48Vb1WVWcA5wL3iUhZ68L1h4iMFZF5dXV1QYdifLBkyZKgQ0gXXd3XZBPaApGVP2M8i6sMerpHLyL73K/hEwugqhqxwRwReRk4MfQUvYh0BhYDJaraK5Yg24rdI0xPNTU15OXlBR1GyrN79MYEK94y6PX1ung6pbkCp2GdTwFU9SsRGYdzdm9Mm+nbt689jGeMabe8Xro/U1U3NR6AnzQ3g6r+r6p+2ijtW1V9sDUBG2OMMcY7r2f01wEzI6RfC8z2ujAR6QpMAYYBueHjVPX7XvMxxhhjjDdRK3oRGeN+7SAiJ+Lclw8ZBHwV4/L+htMJzmPAv2Oc15i4TJxo3ZwbY9qvls7o57ufBwALwtIVqCX2luxGAj1UdW+M8xkTN2sZzxjTnkWt6FU1H0BE7lfV8xOwvFXAEcCbCcjLGE9GjBjB6tWrgw7DGGMC4fVhvNki0j88QUT6i8gxMS5vArBARO4UkevChxjzQUQuFpFKEdkjIgsbjTtJRN4VkV0i8oKIDIg1f5M+Xn/99aBDMCYhtm3bhoiwadMmwOmw6brrrmPgwIG88cYbAUdnkpXXh/EeBE5vlJYNPAAMjWF5fwD6AxuBLmHp8bz7VAPcBJwCHBhKFJGDgX8AFwJLcbrSfRjntoExxqSsqqoqunXrxoABA9i5cyfnn38+tbW1VFRU0KtXUjZPYpKA14r+UFX9IDxBVatFZGCMyzsHONzttrZVVPUfACJSDPQLG/Vj4C1V/Zs7fgZO17qDVfXd1i7XpJ4+ffoEHYIxCVFVVUVRUREffvgh48aNY9iwYSxevJjs7OygQzNJzOul+49FZHh4gvt3TYzL+wCnT3s/HQXUX8NyO96pdtNNO1RTE+thakxyWrNmDXv37mXkyJGMHz+eBQsWWCVvWuS1or8NeEJELhGRU0XkEpxX5Dy/Q+96AFgiIueKyJjwIcZ8oskFGjeaXYfTrW4TIjLJvddfuWXLlgSGYZLFjBkzgg4hXRwcKivuMKm1GVr5i01VVRVvvfUWQ4cO5YorrqhPr6uro6SkhNzcXNatWxdghMZncZVBz/3Ru93L/grnHvtHwH2q+vdYIhSRDc2MUlUdFEteYXneBPRT1Qnu33OALFW9KGyatcAMVX00Wl7W1nZ6sv7oE8Paug/Wnj17yM3N5ZlnnmHChAlcfvnlXHbZZQDs3buXL7/8kiuvvJJp06YxZMiQYIM1vvC7rXvce95/i3UBjfKIp838WL0FXBD6Q0RygAI33RhjUtK6devo0KEDpaWlPPbYY5xwwgkMHTqUMWPGkJWVRc+ePYMO0SQpT5fuxTFRRJaLyJtu2mgROcvf8KLGlCkiB+C0tNdBRA4QkUycWwpDROQn7vjrgDftQTxjTCpbs2YNQ4YMITMzk+HDh3PnnXdy1llnsWFDcxdKjXF4vUd/I85l+3uBQ920j4Gr/AjKo2txmtG9Gvi5+/1aVd2C09nOH4BtwHE4T/ubdsouB5t0EHriPuT888/nvPPO44wzzmDnzp3BBWaSntdL9xOAYar6uYjc5aZtwGnvPhCqOgOY0cy454DBbRmPMca0RnV1NTU1NeTl5VFQUNBk/B133NEkbe7cuW0RmklxXs/oOwA73O+hp5pyw9KMSVrFxb49P2ZMq1VUVFBaWkphYSGjR4+msLCQ0tJSKioqYsrn1FNP5ZlnnmHixIksXLjQn2BNSvJ6Rr8Mpxncy8G5Z4/T4tzSWBYmIrer6mUxRWiMMWmqoqKCsrIyevToQXl5OUVFRVRVVTFnzhzKyspYsWIFJSUlnvJatmyZz9GaVOW1or8CuB/nffQsnDP5Z4BmO7oRkc+AD4F97O/e9jAR+Q9V9XbkGmNMGps6dSo9evSgsrKS3r17A3DyySczfvx4iouLmTZtGitXrgw4SpPqPF26V9XtqnoGzoN4I4ECVf1PVY3WH/2lOG3az1LVY1X1WOAVq+RNW7v++uuDDsGYJqqrq1m1ahWXXnppfSUf0rt3b6ZMmcJLL71EdXV1QBGadOH1Hj0ichDwPaAMOElEukWbXlUXA2cCmSLyuIicx/4ze2PajLWMZ5JRqGnm8Cfpw4XSrQln01pe36Mfg3N2PgU4FrgE2CAiJ0WbTx2LcDqayQSqWhOsMfHIy8sLOgRjmggdl1VVVU3SwtPt+DWt5fWM/g5gkqoep6pnqepIYCJwp5eZVXWfqt6vqlfHG6gx8frkk1Z3lmhMwhUUFDBq1CjmzJlDbW0tsP9Yra2tZe7cuZSWlkZ81c6YWHit6POAxu3EPwYc4nVBInKEiNwsIk+IyPPu580icoTXPPwiImNFZF5dXeO+cIwxYbqKyDwRGZvITNtz+Zs1axZbt26luLiY8vJyAMrLyykuLmbr1q3MnDkz4AiTy86dO7n22mspKCigc+fOHHnkkdxzzz1Bh9WW4iqDXiv6+4HfNEr7tZveIhE5F3gFp9/4lcBDwItAX+BlETnbYxy+UNWlqjqpa9euQYZhfDJ8+PCWJzJe1KnqJFWN6bXalrTn8ldSUsKKFSvIz89n+vTpAEyfPp38/PyYXq1rD7Zt28aoUaPYsGEDy5cvZ/v27dx777387ne/Y/78+UGH11biKoOeeq8TkVU4TcnWAptxKuhewKvsb0AHVR3dzPwbgJ+r6j8jjPsusEhVB8YSuB+s9yxjmme91/mrpZbx2rvx48ezZcsWnn76aZymXBy33HILjzzyCK+//nqA0bUNv3uvu9cd4tUTaG4vrAEObkXexkQ1adIk5s2bF3QYxkRVUFDArbfeasdqBBs3buShhx7itddea1DJg7PdrGOf6Ly+R/8XL0OULJ4FFohIg39T3b/vdccb44t7723N/6jGtB07ViN75pln6N+/f8TbcJs3b6Znz57U1dVRUlJCbm4u69atCyDK5OX19bpzQw/Nich3RORF94E6rx3H/NL9fFtEdopIjYjswOkjXsLGG2OMMQ1s2bKFfv36RRz32GOPccopp9CpUyeeeuopfvrTn7ZxdMnP66X7m4D/cL/PBF7DaQb3T8CYlmZW1W3AuSLSCTic/R3i/EtVd8UatDHGmPYjPz+fTZs2sW/fPjIy9p+fPvvss6xevZqFCxeSlZVFz549A4wyeXl96r6nqtaKyAHAKOAanD7qi2JZmKruUtUqVV3lflolb3y3efPmoEMwxhM7ViM7/fTTAbj22mvZtWsXe/bs4cEHH+Scc85hwYIF5OfnBxxhcvNa0W8RkULgh8BrqroHOIAENGkrIh1E5LrW5mNMc1avXh10CMZ40l6P1erq6qjt+ufm5rJ8+XLefPNNBgwYwCGHHMKiRYt48sknOeuss9o42tTj9dL974HVwLdA6J33k4A3EhTD9ThXCIxJuNNPPx0vr5EaE7T2dqxWVFQwdepUVq1aVZ82atQoZs2a1aQNgcMPP5wnn3yyrUNMC54qelVdKCKPuN9Dl9tfBc7xMr+ILGhtDMYYY9JHRUUFZWVl9OjRg/LycoqKiqiqqmLOnDmUlZXF1WDQqaeeSlVVFe+99x6TJ09mwoQJ/gSfYjw1mNPqhYjsBuYDX0QY3QG4SlU7+B5IC9p7gx2JtGjRImbPnk11dTUZGRkMHTqURYsW0bdv3zaPRUTa1VmSX6zBHP+1p2O1tLSUjRs3UllZ2aCb3traWoqLi8nPz2flypUBRph8/G4wp7XWAv+jqksaj3Af8Au0sxu33eCxhYWFQYaRNhYuXMhNN93Eww8/zPDhw9m6dSuPP/443bpF7dnYN+2sLWw/dRWRecDSRDaDa+Vvv/ZyrFZXV7Nq1SrKy8sbVPIAvXv3ZsqUKUyfPp3q6mprJbChuMqg5/7oW2lhlGXtBW5oozgias9tbfth/vz5TJ48mREjRiAiHHzwwVx44YV06tQpkHgmTZoUyHLTkLV177P2cqzW1NQAUFRUFHF8KD00nakXVxlssaIXkQwRGSMi2fFGpqp3qurjzYz7VlUDrehNYh144IEsWLCARx55hM8//zzocJo0mWlMsmovx2peXh4AVVVV9Wnh6x5KD01nWqfFil5V9wFPqOrX8SxARH7vcTqr7NPE/fffzymnnMK0adPo3bs3Y8eO5bPPPuOVV17h+OOP54QTTuDcc89l7969QYdqjAlAQUEBo0aNYs6cOdTW1jYYV1tby9y5cyktLbXL9gni9dL9ShEZGecyLhORfBEZFG0ApsSZv0kyhxxyCLfffjsffvghFRUVvPnmm9x6660MGDCA559/nhdffJFBgwbxxBNPBB2qMSYgs2bNYuvWrRQXF1NeXg5AeXk5xcXFbN26lZkzZwYcYfrw+jDeJuBpEXkC+IiGXdO21NhNDrCelhvX2e0xFpNCRowYwdFHH83OnTsbXIbLzMxs0JSln370ox+1yXKMaa32dKyWlJSwYsUKpk2bxvTp0wGYPn06paWlPProozG/Wmea57WiPxB43P0euWeBZqhqWz3wZ5LALbfcQmlpKcXFzhsgixcvZsWKFSxfvrx+mg0bNvD0009zzTXXtElMS5cm9NkxY3zT3o7VkpISVq5cSXV1NTU1NeTl5dnleh94bTDnF34HYtLD9u3b+cUvfkFNTQ05OTkMHz6c5cuXc9xxx9WPv+CCC3jggQfIzo77+c6YjB07tt39gJrU1F6P1YKCAi677LJ2ue5twXODOW43tT8FeqvqxSLyHaCjqr7pZ4BtyRrsaFlr/vP+5ptvGDduHFOnTmXMmBY7PUyY9tQIiZ+swRz/tedjtT2vu1fxlkGv/dGfCawE+gLnu8mdgdmxLtCkpoqKCkpLSyksLGT06NEUFhZSWlpKRUWF5zwWL17Mq6++yo033khZWRkPP/ywjxEbY4wB7/fobwS+p6pVIhLq1OYN4Bh/wjLJJFFtUo8fP57x48e3QcTGGGNCvFb0vdjfU52Gfdp1lnZg6tSp9OjRo0Gb1CeffDLjx4+nuLiYadOmJXWb1HY50KSK9nystud195vXJ+JXA41Pxc4BvF+3TWIiMlZE5tXV1QUdStIJtUl96aWXNtsmdbR+pJPBvHnzgg4hXXQVkXlu2/QJY+Vvv3iO1W3btiEi5Obm0qlTJ/Ly8rj99tsTH5zPrJx6ElcZ9FrRTwFuEpEXgRwR+R+cPuovjzHIpGRtbTcvUpvUM2bMqP+eCm1ST548OegQ0oW1de+zeI7VqqoqevbsyY4dO9i1axd33XUXl19+OR9//LEPEfrHyqkn/rR1D6Cq7wKDgTuBa4E/A0er6vsxh2lSSqQ2qW+4YX9rxdYmtTHBqqqq4thjj63/O/Qq69dfx9VquUlDXp+6H6qqu1T1EVUtV9W/quoOv4MzwbM2qY1JbmvWrKl/GPbLL7/kmmuuYcSIEeTn5wccmUkWXi/dPykiW0XkcRG5XESGS3vpZsk0aZP6hhtuSKk2qZcsWRJ0CMZ4Es+xWlVVRXl5Od27d6/vGnrp0qVs376dkpIScnNzWbdunQ/RJpaVU/94bRnvULfjmdHACcDFQA8RWaWq7adx5nYqUpvUQMq0ST1ixIigQzDGk1iP1T179vDOO++wYcMG+vVr2Dr53r17eeqpp7jyyisTGaJvrJz6x+vrdajqByKSCWS7ww9wXrsz7UB4m9SFhYWsX78+ZS7X9+3b117dMSkh1mN13bp15OTkNKnkAbKysujZs2ciw/OVlVP/eL1H/1cR+Qi4HxgELAIGqmpyn8qZhAtV7qlSyRuTztasWcNRRx0VdBgmyXm9R18MfIvTaM4bQJWqfuVbVMYYY1pUVVXFkCFDgg7DJDmv9+gLReQQnPvzo4GrReRAYKWqXuhngCb5TJw4MegQYpJq8Zr2K/xY9dKB1B133NFWofnOyql/PPdeByAiRcCJQJn7+ZWq9vUlsgBY71nGNM96r2sbFRUVTJ06lVWrVtWnjRo1ilmzZsX84Oupp55KVVUVAwYMYPLkyUyYMCHB0Zq25HfvdUtE5AvgCWA4sBQYkU6VvPEu1Z6OTbV4Tfs1ePBgysrK2LhxI+Xl5Tz77LOUl5ezceNGysrKYuotEmDZsmXU1NTwyiuvJH0lb+XUP57O6EVkAvCiqm7wPaIAuO0Gjy0sLJz4/vvW2F9LUq3f6FSLN1mJyHrgBWBpIpvBtfK3n4jQr1+/Bh1IgdM4VXFxMfn5+UndgVRrWDltWbxl0GsTuAuBj0RktIic6356fjUv2Vlb28Z4Ym3d+yjUMVQqdyBlfOdfW/ciMhh4B3gIp4Obh4B3ReSImMM0Ka9Pnz5BhxCTVIvXtE+ROpCC/b26pUIHUq1h5dQ/Xl+v+xMwD+ivqseraj/gbjfdtDOp9kOTavGa9ilSB1IAkyZNapCerh1IWTn1j9eKvgiYrQ1voNzuppt2Jryb2lSQavGa9qmgoID+/fs36UBqxowZ7aIDKSun/vH6MN46YIqqPh+WdiJwh6qmTbNM9nqPN6n20EyqxZus7PU6/4kIHTt2pEuXLvzyl7/kpJNO4vvf/z79+vVj69atrFixIun7loiXldOWxVsGvT5Q91tgiYg8CWwCBgCnAT+PdYEm9XhpuKO9Cm2bvXv3kpWVZdvIxKW6uprnnnsOcDqq2bJlC7feeiu33norAPn5+SnRgZRJTl5bxlsiIsOAs4E8YB1wnar+y8/gTLD+8Y9/MGPGDNauXVufNmrUqAAjik2oEvZDRUUFF110EatXr24yLt7GTUz7E6mM5efnc8UVV7B582bmz5/Pli1bmDlzZloeT+EnEY3/DuIf5nT9xz3qpXsR6QRcCwwBXgf+qKp72ii2NnfooYfq1KlT6d27N7169eKzzz6jtra2VX8DCZk2EeOa+954nvXr1/PnP/+ZDz74oH7bHH744QwZMoSXX36ZL774grlz53LYYYd5XkYivmdlZbF3794Wp/3iiy+4++67G/x4FhcXM27cOAoLCz3lFW3866+/zsyZMxtcZuzWrRuqyo4dO+jUqRO7d+/mtttu4+ijj26ST6S8o6UlcvrWjBs9evS/VPU7jctNogwaNEjvu+++uLZVsm3flqZZu3Yts2fPbvCqXGZmJqeddhqrV6+uv0w/YMAAjjrqKPr27cttt93W7HK8frYUc6jCDa/sEv2Zl5fHG2+80eQfnJycHHbu3Fn/99ChQ5k8eTJHH320p7haM83777/PvffeG/Ef96FDh3L99ddzzDHHNPknoLnlRRqXiLTCwsL3VHVwC0WpKVVtdgAWAGuB/4tzFv/f0aZPpgHoDjwG7MS53XCeh3nUhqZDx44dA48h3mHQoEHatWtX3/LPzMzUX/3qV9qrVy8VkcDX1+/B5zIb+PoFMQwaNEgB7dq1q3bs2FGXLVum/fr106KiIh01alTg8bXF+nfo0EEzMjI0MzOzfnsENXTt2lW/+93vamZmZuDbJtIQV9lqoeB9AvRxv/cHNgRdgcfwo7EYeBjIBUYBdcBR9kMT+5CRkaEdOnSo/x5KD6Ul6xAeX+iflURXxjk5OdqxY8f6H4VQ/sm+beIdrKJP7JCRkaEHHXSQAnrzzTcroKWlpXrxxRcroN27d6+fVkR8+WdSRBqUaz+WF2kZ2dnZeskll2h2dnb9ckREs7KyNDMzU7t3767Z2dlR4/ISe7RpQumh6TMyMrRjx47asWNH7dOnj+bm5jbZDkOHDm02r0jLiTR9a9L8qOi3N/r7i6ArcI8/GDnA18DhYWkPALe09EMTvkPDv2dkZHj+OysrS7OysurHNf7b67SJGpeTkxPxe+N5QpVV6GAN/T18+HDNzs5uUnkdeOCBmpWVpRkZGc0uo7XfQ9sqVJDCv0eaFtBDDz1U8/LyGhTQzp0715/ZZ2ZmRsyrpWWFF+IuXboooEcccUT9D1Xnzp3rx/fv37/BD0h2dnb9dKG8o6VlZ2c32Dfh6bFOn4hxoW3td0Xf3H6ItO4t7bsgtm9oO0WbJny6nj17alFRkQKan59fP/6www6rP75CaYcccogOGzZMO3TooDk5OZ4/W4onVL4jDfEsL9JnaBmhZR566KHar18/PfbYYxWoL5uhz0MOOUQ//fRTHThwYMT4Q3F5iT3aNE899ZQCmpeXVx9XaFz48vv06aOA9urVS7t06aJ9+vRpktfatWubpOXl5TWZvjVpwNd+VPS7cHqpG+MO2xv9PSboSr2ZuIcB/26UNg2nfeCoPzTNDSNHjozp71jySvXhxBNP9C3v8G3V3PfQEPrBBHTSpElNxp922mkRYw7l1dKywvOPtPzw/KNN5zUtGY8hvyt6L/s50vhY5kmG7Rs6VsrKyhTQvn37tjjP1Vdf3eDYbbzuzX22doh1eY0/Q+vYeLjqqqvqv5999tkNtsuvf/1rXb9+vadt2JppFi9e3Oy45cuXN0k777zzFNDJkyd7mj70OxQ+fWvSgM/8qOg3AhuiDB8EXak3E3cp8GmjtInAigjTTgIq3SHaQbEtxr9jySsVh6/DDz4fl7PNw/fQsDXse22Eab9qJuZdHpcVnv83EZa/w+O2iRR7rMdErNMnaqgMGyYloKxGKn8t7edI42OZJxm2byj/0HEUOl6/DZtmr/u5z/38kIbHVON1b/y5I8ryYxm8Lq+5z+Z+Hz4M+x6aNlRGPwDebSGujzzE3tI01VHGbYqQVtNof7U0fW2E6VuT9g1xlMHAK2WfKvphwK5GaVNp4Yw+hvznxTBtpQ95JnzaVFknn/L0tD62Tv4OscTs57okQx7psi0SFINti1ZuC69N4KaafwGZInJYWNoxwFsJyj+hvXfFkacf06bKOvm1nfzIMx3XKdklYl2SJY/WSob1SIbtAO18W3hqAjcVichfcS51XIjTJv8y4D9UNVGVvdc4KtXHZkODkG7rlG7rA6m5TqkYs19sW+xn22K/eLdFup7RA1wEHIhzf2gx8Ou2ruRd8wJYpt/SbZ3SbX0gNdcpFWP2i22L/Wxb7BfXtkjbM3pjjDHGpPcZvTHGGNPuWUVvjDHGpDGr6GMkIheLSKWI7BGRhWHp2SLydxHZ6LbQVdZCPitEZLeI7HCH93wOPVosza3TSBF5VkS+EJEtIvI3EekTJZ/uIvKYiOwUkU0icl6brEDTOBK1Pqmwj45007e5w3MicmSUfALfRyLSUUTmu8v/SkTWiMgP3XEx7aNUF21bNJruevd35eQg4vRbS9tBRDqJyJ9E5HMRqRORlUHG6ycP2+IsEXnHHfe2iJzRUp5W0ceuBrgJp8OfxlYBPwc+9ZjXxaqa6w6+9QrmQXPr1A3n4Y+BwACcxiz+HCWfO3Eai+kN/Ay4S0SOSnSwHiRqfSD591EN8FOcTpwOBpYAf42STzLso0ychkxOALoCvwMeEZGBxLePUlm0bQGAiBTg7ONPggiwjbS0HebhHONHuJ+XBxBjW2l2W4hIX+BB4AqgC3Al8JCI9IqaY2tf4G+vA86P7sJmxn0MlLUw/wrgwqDXw+s6ueOHA181My6u/gWSdX1ScR+5PxC/oVFjUcm8j8LieBP4Saz7KB2HxtsCeBo4Fael0pODjq+ttwPwHZzm17sEHVMSbIvjaNQMLrAFOD7a/HZGH6w/upei/tnSpf4kMZrmGx06HPhWVf8VlvYGEMQZvVfR1ickJfaRiHwJ7Ab+G7i5mcmSch+JSG+c2CLtCy/7KG003hYiciZORybLAg2sjTXaDsfhNC97g1sW14rITwINsA012haVwDsicrqIdHAv2+/B+UegWZm+R2macxXwNs4Z1jnAUhEpUtXqYMOKTESGAtcB45qZJBenK+BwdUBnP+OKl4f1gRTaR6p6kIjkABfg/ChGknT7SESygEXAX1T13UbjvOyjtNF4W4hILs4/bd8PNrK2FWE7/BgYAjwK5AHHA0+JyNuq+k6AofouUvkQkfuBh4ADcH6bzlTVndHysTP6gKjqq6r6laruUdW/AP/EuTyXdESkEOfy4aWq+lIzk+3AuWcUrgvOPdak4nF9UmofAbiF/W7g/mbu2SXVPhKRDJxbB18DFzca52kfpYtmtsUNwAOquiGwwNpYM9vh3zgd/Nykql+r6ovAC6T5P0CRtoX7MOb/BcqAbJz7+PeJSFG0vKyiTx4KSNBBNCYiA4DngN+r6gNRJvW7f4GEiGF9IknKfdRIBtAJ6BthXNLsIxERYD7OQ4E/UdW9YeNas49STpRtcRIwRUQ+FZFPgf44D2VdFVCovoqyHaJelk5HUbZFEbBSVStVdZ+qvga8CkR9G8Mq+hiJSKaIHAB0ADqIyAEikumO6+iOA8h2xzWpGETkIBE5JTSviPwM517k/7TZijSMJ+I6uU94Pg/cqap3R8vDPZv8B3CjiOSIyHdxLrm2+Q91ItYnhfbR90RkmHu/rgswG6fLzyaXNJNpHwF34TxBPVZV/x1KjGUfpZGI2wKnoh+C8+NehPOGxWScNyfSUXPbYSVOl7b/5R7z38U5ow2kLLaR5rbFa0Bp6AxeRIbhdMse/Z+hoJ8mTLUBmEHTPodnuOM2Rhg30B33W+Bp93tPd4d9BXwJ/C/wvWRbJ+B69vdrXT+EzVe/Tu7f3YHHgZ04BfO8VF2fFNpHZ+L0270D5+nbZcDQZN5HOK/NKc7Dg+H74mct7aN0G6JtiwjTbiRNn7pvaTvgPDD6invcvg38Z9AxB7gtLgbWu79NHwBTW8rT2ro3xhhj0phdujfGGGPSmFX0xhhjTBqzit4YY4xJY1bRG2OMMWnMKnpjjDEmjVlFb4wxxqQxq+iNMcaYNGYVvUkoEXlaRC6Ic96NblvOkcatEJHdIrLSy/SJJiIni8gOEdnXVss0Jh5WBk1jVtGbUGH9t1uIakXkz27PWTFT1R+q0wGMHy5W1dE+5R2Vqj6nqrk4rckZk1BWBltmZTB+VtGbkLFuIRoOHAtcG8vM4kjL4ynUl4ExPrMy2Awrg62TlgeFiZ+qbsbpHnQIgIiMFJGXReRLEXlDRMpC07qX8v4gIv8EdgGD3LQL3fEZInKtiGwSkc9E5H4R6Ro2/3h33FYRuSbOkItE5E0RqRORh8M6FUJEJorIehH5QkSWiEiemz5QRDT8x6NR3BNE5J8icpuIfIHTprwxbcLKoJXBRLOK3jQgIv1x+lxfI05PYk8BN+F0hjINeFREeobNMh6YBHQGNjXKboI7nAgMAnKBO9zlHInTQ9N4IA/oAfSLI+SzgB8A+cBQd3mIyBjgj+74Pm5sf40h3+NwOozoBfwhjriMiYuVwXpWBhPEKnoT8riIfAmsAl4EbgZ+DixT1WXq9H38LFCJ8yMUslBV31LVbzSsT3HXz4DZqvqBqu4A/gs4x/0v/qfAk6q6UlX3AL8D9sUR91xVrVHVL4ClON15hpa9QFVfd/P/L+B4ERnoMd8aVf1vd73+3fLkxrSalcGGrAwmiN33MCFnqOpz4QkiMgA4U0TGhiVnAS+E/f1RlDzzaHiGsQnnmOvtjqufV1V3isjWOOL+NOz7Ljff0LJfD8t/h5t/X2Czh3yjrZcxfrAy2JCVwQSxit5E8xHwgKpOjDJNtH6Oa3D6Vg45FPgGqAU+AY4IjRCRTjiXDhOlwbJFJMfNfzNOn9YAnYDt7vdDGs1v/TebZGBl0LSaXbo30TwIjBWRU0Skg4gcICJlIuL1Pt5i4HIRyRfnVaGbgYdV9Rvg78CPRGSUiGQDN5LY4/Eh4BciUiQiHd1lv6qqG1V1C86Pzc/d9folUJDAZRuTKFYGTatZRW+apaofAeOA3wJbcM4ursT7cbMAeABYCWwAdgOXuHm/BfwG58fgE2Ab8HECY1+Oc8/xUTf/AuCcsEkm4qzLVuAo4OVELduYRLEyaBJBVO3qiEl+IvIMcDxQqaonBrD8k3B+sDoCp6rqCy3MYkxasTKYuqyiN8YYY9KYXbo3xhhj0phV9MYYY0was4reGGOMSWNW0RtjjDFpzCp6Y4wxJo1ZRW+MMcakMavojTHGmDT2/wHaPXKfzOuWLAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "## calc PSD\n", "fs = 1. / step_sec\n", "freq, psd = signal.periodogram(tide_los, fs=fs, scaling='density')\n", "# get rid of zero in the first element\n", "psd = psd[1:] / 1e4\n", "freq = freq[1:]\n", "period = 1./3600./freq # period (hour)\n", "\n", "## plot\n", "fig, axs = plt.subplots(nrows=1, ncols=2, figsize=[7.2, 2.8], sharey=True)\n", "for ax in axs:\n", " ax.plot(period, psd, 'ko', ms=8, mfc='none', mew=1.5)\n", "\n", "# axis format\n", "for ax in axs:\n", " ax.tick_params(which='both', direction='out', bottom=True, top=True, left=True, right=True)\n", " ax.xaxis.set_minor_locator(ticker.AutoMinorLocator())\n", " ax.set_xlabel('Period [hour]')\n", "axs[0].set_xlim(11.3, 13.2)\n", "axs[1].set_xlim(22.0, 28.0)\n", "ax = axs[0]\n", "ax.set_ylabel('Power spectral density\\n'+r'[$10^4$ '+r'm$^2/$ Hz]')\n", "ax.set_ylim(0, ymax=axs[0].get_ylim()[1] * 1.1)\n", "ax.yaxis.set_major_locator(ticker.FixedLocator([0,10,20]))\n", "ax.yaxis.set_minor_locator(ticker.AutoMinorLocator())\n", "fig.tight_layout()\n", "\n", "# Tidal constituents\n", "for ax in axs: pysolid.add_tidal_constituents(ax, period, psd, min_psd=0.4)\n", "\n", "# output\n", "out_fig = 'SET_PSD.pdf'\n", "print('save figure to file:', os.path.abspath(out_fig))\n", "fig.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=600)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Frequency Aliasing" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Signal frequency [mHz]: 0.02236427796891005\n", "Sampling rate [mHz]: 0.0019290123456790122\n", "Alias frequency [mHz]: 0.0007838701792380956 -> 14.77 days\n" ] } ], "source": [ "t_m2 = 12.4206012 * 3600 # seconds, period of M2 tide\n", "t_sen = 6 * 24 * 3600 # seconds, sampling interval of Sentinel-1 and NISAR\n", "f_m2 = 1000. / t_m2 # mHz\n", "f_sen = 1000. / t_sen # mHz\n", "print(f'Signal frequency [mHz]: {f_m2}')\n", "print(f'Sampling rate [mHz]: {f_sen}')\n", "\n", "N = np.rint(f_m2 / f_sen)\n", "f_alias = np.abs(N * f_sen - f_m2)\n", "t_alias = (1000. / f_alias) / (24 * 3600)\n", "print(f'Alias frequency [mHz]: {f_alias} -> {t_alias:.2f} days')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Impact of SAR Acquuisition Time Variation" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "max LOS SET rate: 1.9 mm/min\n", "extract metadata from ISCE/topsStack xml file: /Users/yunjunz/data/geolocation/ChileSenAT149/reference/IW2.xml\n", "This is the Open Source version of ISCE.\n", "Some of the workflows depend on a separate licensed package.\n", "To obtain the licensed package, please make a request for ISCE\n", "through the website: https://download.jpl.nasa.gov/ops/request/index.cfm.\n", "Alternatively, if you are a member, or can become a member of WinSAR\n", "you may be able to obtain access to a version of the licensed sofware at\n", "https://winsar.unavco.org/software/isce\n", "Satellite speed on orbit / ground: 7590.2 / 6831.6 m/s\n", "Time to image 100 km: 14.6\n", "SET ramp in the along-track direction: 0.5 mm / 100km\n" ] } ], "source": [ "tide_los_max_rate = np.max(np.abs(np.diff(tide_los) / step_sec)) * 1000 * 60\n", "print('max LOS SET rate: {:.1f} mm/min'.format(tide_los_max_rate))\n", "\n", "xml_file = os.path.expanduser('~/data/geolocation/ChileSenAT149/reference/IW2.xml')\n", "meta, burst = isce_utils.extract_isce_metadata(xml_file)\n", "Vs = float(meta['satelliteSpeed'])\n", "Re = float(meta['earthRadius'])\n", "Hs = float(meta['altitude'])\n", "Vg = Vs * Re / (Re + Hs)\n", "print('Satellite speed on orbit / ground: {:.1f} / {:.1f} m/s'.format(Vs, Vg))\n", "\n", "t = 100e3 / Vg\n", "ramp = tide_los_max_rate * (t / 60)\n", "print('Time to image 100 km: {:.1f}'.format(t))\n", "print('SET ramp in the along-track direction: {:.1f} mm / 100km'.format(ramp))" ] }, { "cell_type": "code", "execution_count": null, "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.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }