{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Impact of top-side TEC on SAR range delay\n", "\n", "Dataset: Chile Sen asc track 149\n", "\n", "+ Figure. Time-series of GIM/SUB/TOP_TEC and its predicted range delay for X/C/S/L-band SAR" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Go to directory /Users/yunjunz/data/geolocation/ChileSenAT149/offset_comp/BoxLR\n" ] } ], "source": [ "%matplotlib inline\n", "import os\n", "import h5py\n", "from pprint import pprint\n", "import numpy as np\n", "import datetime as dt\n", "import numpy.polynomial.polynomial as poly\n", "from matplotlib import pyplot as plt, ticker, dates as mdates, colorbar\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "\n", "from mintpy.objects import timeseries\n", "from mintpy.utils import ptime, readfile, utils as ut, plot as pp\n", "from mintpy.simulation import iono\n", "plt.rcParams.update({'font.size': 12})\n", "speed_of_light = 299792458 # speed of light in m / s\n", "\n", "# dir\n", "proj_dir = os.path.expanduser('~/data/geolocation/ChileSenAT149')\n", "proj_name = os.path.basename(proj_dir)\n", "\n", "work_dir = os.path.join(proj_dir, 'offset_comp/BoxLR')\n", "os.chdir(work_dir)\n", "print('Go to directory', work_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure - Topside Ionospheric Impact on SAR\n", "\n", "### 1. Read data" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "tec_ipp min/max/span: 6.6 / 59.1 / 52.5 TECU\n", "tec_sub_ipp min/max/span: 0.6 / 49.2 / 48.6 TECU\n", "tec_sub_tpp min/max/span: 0.6 / 53.3 / 52.7 TECU\n", "tec_top_tpp min/max/span: 2.9 / 11.5 / 8.6 TECU\n" ] } ], "source": [ "geom_file = os.path.join(proj_dir, 'mintpy_offset/inputs/geometryRadar.h5')\n", "tec_file = os.path.join(proj_dir, 'mintpy_offset/inputs/TECsub.h5')\n", "\n", "# geometry\n", "inc_angle = 42 # use the median value of NISAR\n", "inc_angle_iono = iono.incidence_angle_ground2iono(inc_angle)\n", "\n", "# TEC\n", "tec_dict = {}\n", "with h5py.File(tec_file, 'r') as f:\n", " dnames = [x for x in f.keys() if x.startswith('tec_')]\n", " for dname in dnames:\n", " tec_dict[dname] = f[dname][:]\n", "\n", "# time\n", "date_list = timeseries(tec_file).get_date_list()\n", "dates = ptime.date_list2vector(date_list)[0]\n", "\n", "# stats\n", "for key, value in tec_dict.items():\n", " print('{:<11} min/max/span: {:4.1f} / {:4.1f} / {:4.1f} TECU'.format(key, np.nanmin(value), np.nanmax(value), np.nanmax(value)-np.nanmin(value)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Plot" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "number of tec_top_tpp observations: 93\n", "min / max: 3 / 11 TECU\n", "mean / std: 6.6 / 1.8 TECU\n", "L mean / std: 1.7 m / 0.462 m -> 0.07 pixel for 24 MHz\n", "L mean / std: 1.7 m / 0.462 m -> 0.14 pixel for 44 MHz\n", "L mean / std: 1.7 m / 0.462 m -> 0.25 pixel for 80 MHz\n", "S mean / std: 30.6 cm / 7.825 cm -> 0.039 pixel for 75.00 MHz\n", "C mean / std: 11.5 cm / 3.033 cm -> 0.013 pixel for 64.35 MHz\n", "X mean / std: 3.8 cm / 1.020 cm -> 0.007 pixel for 109.89 MHz\n", "Ka mean / std: 2.8 mm / 0.773 mm -> 0.001 pixel for 200.00 MHz\n", "save figure to file /Users/yunjunz/data/geolocation/ChileSenAT149/offset_comp/topTEC_pred.pdf\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxIAAADwCAYAAACUoFLfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACz+0lEQVR4nOydd3xT5frAv6cts0DL3rQUKUlLBJXlRONCSx1XRb2oxQXBONDr79LrIjiwehVnNA6Eet3rqqWKV4ly4aoMB4Y2YbS0jLJHoey27++P08Q0Tdq0mW3f7+dzPs15z3ve87znnDTv877PUIQQSCQSiUQikUgkEkljiIm0ABKJRCKRSCQSiaT5IRUJiUQikUgkEolE0mikIiGRSCQSiUQikUgajVQkJBKJRCKRSCQSSaORioREIpFIJBKJRCJpNFKRkEgkEolEIpFIJI0mLtIC1IeiKJlAJnB7pGWRSCQSiSTYCCGUSMsQTBRFWSSEmBBpOYKJ7FPLoDX2ORxE9YqEECJPCDG15rPX7fbbbw/6sVC129Rjp512WpNkbU59jMSzqu++RlMfm9N7HE33NNqeRyD9aMp9bS33pjm/q8DhyP7KBp8BXZSTMSUITAmmSMsSRHpEWoAQ0BL71BCtsc8hJ6pXJPwhMzMz6MdC1W4g8tSHr3ObUx8j8axC0WY0HYvUNaNJlmh6HvL7H13POBRtNvZYbKdubZKz8wUwuyQnw9QYGaOVLQdEGabyfpGWQyKRhAelZlYkqlEURTQHOUPFqFGjWLVqVaTFaHHI+xp85D0NDfK+Bp9ouKeKohwWQsRHVIggoyjKKiHEqEjLEUz86ZMuV9cNmAdcBOwG/mHLsr3npV4WcDcwFDgAvAc8YMuyVQZd8Hpoic+pIVpjn8NBVJs2SVSmTp0aaRFaJPK+Bh95T0ODvK/BR95TSZAxA8eB3sBk4FVdri7dS72OwAxUM5uxwPnA/WGSUSIJOnJFQiKRSCSSCODPikRydv47qIPNeGA78HRJTsab4ZCvKbTEWd+G+qTL1cUD+4DhtizbupqyfwFbbVm27Pra1uXq7gPOs2XZmm7f2ARa4nNqiNbY53AgVyQkEolEIolengSSS3IyugCXAY8nZ+efFmGZJLVJBaqcSkQNqwFvKxKenAMUhESqeugTF9fPrtEKu0ZrCve1JS2LqHa2dgv/KpFIJBJJq6MkJ8N9kClqtiHAL5GRqPUR1y2uny5X524WMduWZTO57XcCyj1OKwc619euLld3MzAKuC0YcjaG7ZWVZVqHXTrFSwImqhUJIUQekKcoiswjIZFIJJIWhVvUJideozclZ+e/AkwBOgC/AV+FRUAJAJV7K8tsWbb6Bt0VQBePsi7AQV8n6HJ1VwA5wAW2LNvugIWUSCKENG2SSCQSSYujuLiY9PR04uLiSE9Pp7i4ONIi1aGqYu+JkpwMxW0zeatXkpNxB+rs9tnAZ8CxMIopaZh1QJwuVzfUrWwEPkyWdLm6CcAbQKYty2YLg3wSIDG+Zz+zwSrMBqsp0rK0JKSztUQikUhaHOnp6TgcDqqrq4mJiUGj0VBQEHZT9HppSvjX5Ox8C1BYkpPxYojECoiBCTFlm+/t3BeYjancFGl5goGf4V8/QDU7uw0YibpqdIYty1bgUU8PfAxcacuy/Tc0EjdMa3Q8bo19DgdRbdokkUgkEklTWLt2LdXV1QBUV1ezdu3aCEsUNOJQfSSiklackO4O4C1gJ7AHmG7LshXocnWDgEIgzZZl2wQ8DCQAX+lydc5zl9qybJdEQGaJJGCkIiGRSCSSFsewYcMoLCwEICYmhmHDhkVYosaTnJ3fC9ADC4EjwAXA9cBfIymXpC62LNte4Aov5ZtQnbGd++eFUSyJJORIRUIikUgkLY68vDy0Wi3Hjx8nNTWVvLy8SIvUFAQwHbCg+jSWAjNKcjK+iKhUklaBXaOtk61b67DXydZdU/deYCZqQIBPgelah/2Y2/HrgFnAINR8KFO0DvvS0PZAEg6kIiGRSPymuLiYzMxM1q5dy7Bhw8jLyyMlJSXSYkkkdUhJSaFz584MGTKEu+++u1m+pyU5GbuA8ZGWQ9Jqcc/WPRLIt2u0q7UOey2/D7tGezGQjbp6Vgb8G5hdU4Zdo70QeAq4FlgB9A2T/JIwENVRmxRFyVQU5fVIyyGRSFQyMzOx2+1UVVXhcDjIzJRpXiTRSVVVFfv372fatGl88sknkRbHK87wr8nZ+aZIyyKRuGPXaOOBq4CHtQ57hdZhXwZ8CdzopXoWME/rsBdoHfZ9wGOo4YqdzAYe1TrsP2sd9mqtw75V67BvDXEXJGEiqlckZB4JiSS6WLt2Lc4Iai3MgVXSwtizZw+JiYlceeWVzJgxg4MHD9K5c735wcJOTfjXtpGWQ1I/bgnpPBPRtWRSgSqtw+6ZrdvbClk68IVHvd52jbY7sB816d6Xdo12A9Ae+Bz4P63DfiQEckvCTFSvSEgkkugiNTXV9bm5OrBKWge7du2iZ8+edO3albPOOov8/PxIiyRpptQkpFNakhLRJy6un12jFW6byaNKY7J1e9Z1fu6MahbVBrgaNQ/KSOAU4KGAOiCJGqQiIZFI/ObZZ58FQFEUNBpNc3VglbQCnIoEwPjx45k6dWpUJ6eTSMLJ9srKMq3DrrhtJo8qjcnW7VnX+fkgarQxgJe0Dvs2rcO+G5gLXBpQB5qATEgXGoKuSCiKcqeiKKsURTmmKMoCj2PnK4riUBTlsKIo3yuKkhTs6zdXmkMWVonk6NGjANx5550UFBQ0SwdWSevAXZGYP38+Bw8elL49YWBAF6UfpgSBKcEUaVkkAbEOiLNrtP5k6y6oOeZeb4fWYd9T4zOxBTUCWUTZf2hXmdGiV4wWvSnSsrQkQuEjUQY8DlyMGgYMAEVRegCfoWZ9zEN1xvkQGBcCGZoVQgjOO+88Nm3aBOD6oYu2LKwSyYYNGxg4cCBbtmyJtCgSSb24KxIbNmxwlUvfntDSihPStSi0Dvshu0b7GfCoXaN1Zuu+HDjDS/W3gQV2jfZdYBuq2dICt+PzgbvsGu0i4AQwAzU3iqQFEPQVCSHEZ0KIz1EzO7rzF6BACPGxEOIoYAJGKIqiCbYM0YznysNrr73G2LFjXUoEyB86SfSyfv16zjvvPKlISKKeXbt20atXL0BNTqcoChBdvj0yapMkyrkDdUJ4J/A+am6IArtGO8iu0VbYNdpBAFqHfRHwNPA9aq6TUtScEU4eA1airnLYgd+AJ8LWC0lICWfUpnRUT34AhBCHFEUpqil3hFGOiJKZmYnD4aC6uprCwkJmzJjBO++8wyOPPOIqj6YfOonEnfXr15OVlcW3334baVEkknrZtWuX6/9oXl4eF154IcXFxVHl2yOjNkmiGa3D7jVbt9Zhr5Wtu6ZsLqrvg7d2TqAqJXcEX0pJpAmnItEJ2OVR5isCAIqiTAWmhlqocLN27Vqqq6td+ydOnOCqq67ilFNOITMzk8LCQlJSUqLmh04icWf9+vWcc8457N69mxMnTtCmTZtIiySReGXXrl2cddZZgJqcbsOGDQwYMIDPP/9c+vZIJBKvmA3WOtm8jRa912zeZoO1TjZvo0V/rDHtmA3WWagWOhcaLfrvgt6hMBDOqE2NiQCAEOJ1IcQoIcSokEsWRoYO/dNvyX3lISUlhYKCAq655hoeeeQR+UMniToOHz7Mnj17SE5OplevXmzbti3SIkkkPnH3kQA10tiECRNYtGhRBKWSAGBK6IYp4d+YEg5hSijFlPBXP86x1jhxR3X+K0mzxz2b92TgVbPBml6nksHqzOZ9PpAMpKAm3vO7HbPBOgQ1LG6z/jENpyJRy6tfUZR4YAjeIwC0WM477zwSEhKIjY31usR+5pln8r///S9C0kkkvtmwYQMpKSnExMQwYMAA6SchiWo8FQlAKhLRQ51BFqaEOoM1F6aEyUQwga4zIZ0uV2eKlAyS0GM2WF3ZvI0WfYXRom8wm7fRoi8wWvS1snk3op2XUVc0joeiP+EiFOFf4xRFaQ/EArGKorRXFCUO+DcwXFGUq2qOPwL8IYRoNf4RxcXFfPjhhxQUFFBZWek1fKZUJCTRyvr1610ralKRkEQ7O3furKNIXHDBBSxdutQVxlgSAUwJrkEWpvIKTOX1DdbAlJCA6rj797DJ6EFLTEgn8UoqUGW06D2zeXtTcmv5/dZ87m02WLv7047ZYL0GOG606L8KlvCRIhQrEg+hJiDJBm6o+fyQEGIX6j+PJ4B9wFjguhBcP2r5v//7P+677z769+/vs86IESMoKSlh//794RNMIvEDqUhImgvV1dXs3buXHj161Crv2rUrOp2OpUuXRkiy2rTSqE2pQBWmcn8GawBzgFeB7aEWrDXhltnaFGlZwoVbQjrhIzFdsLJ519uO2WDthPpez2hsH6KRoC8VCiFMqI4j3o59B7SqcK9OrFYrv/32G++++2699dq0acPo0aP56aefuOSSSwK+bnFxMZmZmaxdu5Zhw4aRl5cn/S8kTWL9+vWMHTsWkIqEJLrZt28fnTp1qhMMoLi4mPXr13PxxRej1Woj/v+wJUZtciWk+5PZmMpNbvv+D9ZMCaOAM4F7gAHBlbR1U5PZulXl+6hJSFdfn4OVzbuhdmYD/zJa9Bv9kTvaCaePRKulsrKSe+65h2eeeYb27ds3WL8h8yb3XBRpaWnYbDZ2797Nli1b2LBhA2vWrGHVqlUsW7YMvV6P3W6XGV39QGYXrx+5IiFpLnjzjwA1/Pbu3bsRQsj/hyGiJiGd4raZPKr4N1gzJcQArwD3YCqvDJW8Eokb64A4s8Ha5GzeRot+jx/tnA/cbTZYt5sN1u3AQOAjs8E6M0j9CCtRrUgoipKpKMrrkZYjUF577TV69uzJlVde6Vf9M844o15FIjMz06Uc2O12RowYgUajYdy4cUyYMIHrrrsOg8FAdnY2paWlCKFODvmT6K41D6bd76scZNTFXZFwz27dmt+Z5kZreVa+FIm1a9c26v+hk9Zy38LEOiAOU0JDg7UuwCjgQ0wJ21ETmgFswZRwdujFlLQ2jBb9IeAz4FGzwRpvNljPRM3m/S8v1d8GbjUbrGlmg7Urbtm8/WjnfGA4arbwkUAZMA01CEGzQ3H+U41mFEURzUFOb+zZswetVsvixYvR6XR+nbN//34GDhzI3r17vcbpj4uLo6qqyrUfGxtLZaX3CZv09HRXojtFUdBqtRQU+A6U5V4/JiYGjUZTb/2WRGxsbK0cH/Xd19bGwYMH6d27NxUVFcTExFBaWspZZ53F5s2bW/U709xIT0/HbrcjhGjRz+qzzz7j7bff5vPPP69V7vn/MC4ujj/++AONpn6L21C944qiHBZCxAfcUBShKMqqBsO2mxI+AARwG+pA6ivgDEzlBW51FNSoTk4GAitQTZx2YSoPW6Qbv/rUzGiJfWoIf/pck//hLeBCYA+QbbTo3zMbrIOAQiDNaNFvqql7H7XzSBg88kjUacfHNUuA25prHgmEEFG/qWI2T+68805xxx13NPq84cOHixUrVng9dtJJJwnUf8IiJiZGpKWl+WynqKhIpKWliZiYGJGQkCCKiorqvW5sbKyrbUDExsY2WvbmyNatW0VsbKxQFMWv+9ra+PXXX8Xw4cNd+8eOHRNt2rQRJ06caLXvTKhwfmdjY2NFWlpag9/ZxtBanpXFYhG33XZbnXLPe5uTkyN69eolrFar13b27t0rPvjgA9f/hWDfN+CQiILf2GBuwKoG683q0k3M6vK5mNXlkJjVZZOY1eWvNeWDxKwuFWJWl0FezkkWs7oIMatLXFT2qZltLbFPss+R2WRilxCyZs0aPvzwQ+x2e6PPdfpJjB49us6xCy64gL1791JeXu5yoPaFM9FdWVkZw4cPZ+DAgfVet0+fPmzduhWonTCvJVNdXc2UKVMwGo189913FBYWkpycLLOLu+Fu1gTQtm1bunfvzo4dOxg2bJhrlhugW7durplbSePJzMx0zX4XFhai1WqpqqoKSrCEYcOGUVhYCHj/freU4Ay+TJuc/w/dGTNmDNdccw3t2rVjx44dJCUl8Ze//IWff/6Z1atXM378ePr06cOOHTtc73Ww/i86ozYBs0tyMkxBabQ5YCrfC1zhpXwTqjO2t3NKACWEUkkkLQqzwfpfP6seNVr0FzX1OvKXPkQIIbjnnnt45JFH6N69e6PPP/PMM/nxxx/rlFdWVvLFF1+wdOlSn7kovNGvXz9SU1NZsmRJvfVOOukkOnTogKIoXhPmtURefvllDhw4wLPPPktBQQEPPvggV155ZbMcQIWC4uJi7rrrLj7//PNa9uFOh+svv/ySuLg4YmNjSU1NpX///lx77bUcOnQowpI3T9auXVvLxO748eNB89txN/Xx9v12KjHN3U/IlyLhjfPOO4/ExETKysqoqqqiuLiYt956iwcffJAdO3aQl5fHsmXL0Gg0PhOJNpWaqE1Kq1IimiEyIZ2kmTIamNfA9hZqOoamE+klET+Xo0Rz47PPPhPp6enixIkTTTq/qKhI9OvXT1RXV9cq//LLL8Xpp5/epDaffPJJYTQafR7fsWOHSEhIEEuWLBHp6elNukZzo6CgQHTv3l2sW7euVln//v1FVVVVBCWLHtLS0rya0l1xxRXik08+EX/88YdITk52vatHjx4VU6ZMESNGjBAlJSUBXTuUZj7RikajqWVGQxBNaoqLi0Xfvn1Fu3btxOHDh+scbymmT9dff73417/+5Xf9SPWbBkybkmYubJc0c+G8pJkLS5NmLjyYNHPhb0kzF15S3zmR3miB5iOyTy1ja219fnna4sV+1vtPINeRKxIh4OjRo/ztb3/jhRdeIC6uadZjgwcPprq6mtLS0lrl8+bN49Zbb21Sm1dccQWff/658wtVh48//piMjAxOP/10Nm7cyIEDB5p0nebC8ePHmTx5MnPmzKlltpOWlkb37t2jJmlVKPEnGo17ZBv3SDfOFYm8vDwyMzNRFNXqoF27drz11ltkZWUxbtw4li1b1mT5omGGPNwRey699FISEhKIjY2lbdu2tUzE6ktm6Q/r1q0jLS2NIUOGsG7dujrHhw0b5nqOzdm0sTErEqD223mfo6zfccBmYDyQADwMfJScnZ8cSaEkkuaIW0I6U6RlCQdGi/58P+s12awJkCsSoWDOnDniiiuuCLidq666Srzzzjuu/bKyMpGYmCgOHjzY5DaHDRsmVq5c6fXYGWecIRYuXOj6vHjx4iZfJ5pxznIriiI6deokNmzYUKfOk08+KaZNmxYB6YJLQzP6Tkd8amZi27ZtW6duSkqK1xWJnJwccf/994tx48aJ//znP16vv2jRItGrVy8xZ86cJq0sRMMMufs9CrUT/oEDB0TPnj2F3W4XQtR+fsnJyaJ79+5i3rx5TW7/xRdfFNOnTxdXX321eO+99+ocLyoqEt26dROASEpKiooVoKasSo0YMUL88ssvIb1GMKAJztZJMxf+kTRz4VWNPS9c24AuSlmNU7Qp0rIEa6MFzmT3iYsrKxymEYXDNC3mObXG5xgNmyv8q6Iob/upexwTQtwekPbiJ4qiZAKZwO1OOaOdrVu3MmLECJYvX86QIUMCauu5555j/fr1vPLKKwA89dRTrF+/njfffLPJbWZnZ9OmTRsee+yxWuUbN25kzJgxlJWV0aZNG+677z569epFdnZ2QH2IRvwJgVlSUsLo0aNd96O50lDYSs9Qwk7c6z7wwANYLBYOHDhQywH33Xff5a233uLXX39lx44dtG3rPUHvunXr0Ol0nDhxot577g2tVovD4agjUzhpTLjlQHn66af59ddf+eCDD7weX7t2LRMmTOCKK67gP//5T6Odou+66y6GDBnCvn37qK6urvN/AODcc89l+/bt3H777fztb38LuE+B4hmytU2bNg06n/fr14/ly5c3GFwi0jQ2/Gtydn5voBQYWZKT4QidZE2nJYYVlX1qGbTGPjsxG6wjgOdQwy07AxoogDBa9N5/vP3E3bTpWqDIj+3aQC7YGIQQeUKIqeG6XjD4xz/+wdSpUwNWIqB2hmshBPPmzeO2224LqE2neZMnH3zwAddcc41r0Dx27FiWL18e0LWikTVr1lBYWOicnfCZlCo5OZmhQ4fy7bffhlvEoOLuuOutr+6mLO641/3ll1+YN29eHef+AQMGYLVaufDCC30qEQCpqalUVVU1eM+9cf/99wNE1Pk/XOY+hw8fZu7cuTz44IP1yvLTTz9hsVgoLCx0JaXUarV+mV6tX7+e1NRU0tLSXNGb3BFC8Mcff3DzzTfz66+/BqVfgeL+DgshGnQ+F0Kwe/fuRpk2RQpn1Ca3zeSrbnJ2fhvgXSA3WpUIiUQStbwP/A84B9DWbJqav4HhtuSzwc+lIUe4l01oJqZNP/30k+jXr584cOBAUNo7fvy4iI+PF/v37xdLliwRaWlpdZyvG0tVVZXo06ePWL9+vausurpapKeni6VLl7rKNm7cKPr27Rvw9aKFffv2iXvuuUf07NlT9OnTxy9TlZdeeklMnjw5zJIGF1+O0k6KiorEwIEDhaIoom3btq77oiiKSEtLE0eOHBGdOnUS+/btq3OeM59Jv379GjQFcZqSNXTPPbn33ntFSkqKyMrKaky3g4rzHgFi4MCBITN7ee6558SVV17pV11Pk6/6nrE7gwcPFuvWrRN//PGH0Gq1dY5v2rRJ9O7dW6xevVpoNJom9yWYuL/Dnps3U7d9+/aJzp07R0DSxoOfpk1JMxfGJM1c+EHSzIVfJc1c2MafcyK10QLNR2SfWsbWGvvs3F6etnjvy9MWK6Fo27UiIYQ4yU/Fo/4UoK0QpzPm6aefDqiOfsGgTZs2nHbaafz888+8+eab3HbbbV5njxtDTEwMl112GV988YWrzGazcfDgQc444wxXWVJSEpWVlWzZsiWg60Wa6upqcnNz0Wq1HD58mIKCAv73v//5FcrxmmuuYeHChRw+fDjMUgePvLw8OnbsCHgP95mSksKDDz7Irbfeit1ud92X2NhY/v3vf7Ns2TJ0Oh2JiYm1zsvMzHTNfG/fvr1BJ+i8vDx69+7d6JUFq9XKLbfcQklJiX8dDgEpKSkuE7/bbrstJGGBjx49yj//+U8efvhhv+q7Owe7U99qz7FjxygrKyM5OZnU1FQ2btzI8eO1kwOvXr2aESNGoNVqKS0tpaKiovGdCTLPPfccbdq0cTmfO/G1OtRYR+toJzk7X0EN09gbuKokJ+NEhEWSSCTNj1zgr6FoWEZtCgITJ050JZ3zZ1DVGM4880y++uorvvzyS2688cagtOlp3vTee+9x/fXX1xqYKIrS7M2bfv31V8466yxeeeUVvvzyS15//XV69uzpSkrVUB6O3r17M2bMGBYuXBhmyYNHSkoKPXv2pFu3bj77WlZWRr9+/WrdlzFjxuBwOPjPf/7DRRfVDejQkMmUNzk++OADTj/9dL9zn+zevZuNGzdy9dVXN1qRCHakpa1bt5KWlsYff/wRUDu+mDdvHqeeeiqnnHKKX/Xz8vJcSp97ZCdFUXyaXhUVFZGUlESbNm1o164dgwYNYsOGDbXqOBWJNm3aMHz4cFavXh1Yx4LAtm3buOqqq6isrMRut9OrV696FdKWpkgAr6KaH2SW5GQcibQwEomkWZIDPGY2WAvMBqvVfQu04VqKhKIomxVF2eSxFSmKYlUUJSwO1s2NnTt31srq2xj7b38YMmQIL774IgcOHGD8+PFBCT2p1+ux2Wzs3LmT6upq3n//fSZPnlynXnNVJPbs2cP06dO59NJLufXWW/npp5+8Zgj3h+uvv573338/yBKGj6qqKsrKyqioqODIEe9jEKci4c5tt93GG2+8wbfffsuFF15Y55ymhMtMT0+noKDA9V1piO+//56zzz6bwYMHs23bNq8Ozr4UhszMTOx2e9DCxpaVlTFhwgSfikRTFZfi4mLS0tK48847KSgo8Ps8d6XPuZIUExNDx44dfa72rFu3jtTUVNe+Nz8JpyIBcOqpp0aFn8TKlSsZM2YMoPbbbDZz+eWX+1RIW5IikZydnwRMQ3WQ3J6cnV9Rs9X9hy2RSCS++QTYiDox8a7HFhge9mPjvWwXAFOBQuD/ImTXJqKRZcuWiQEDBogePXqELDzksGHD/LZ/bgyTJk0Sb775pli6dKkYPny41zr/+c9/xNlnnx2U64WaoqIiodVqRUxMjIiNjRU33nij2Lt3b8Dt7t+/X3Tp0qWOj0BzYfPmzaJv374iOTnZa5hbIYS49NJLxZdfflmrzGazud5prVZbxy+gqeEye/XqJbZs2eJXXYPBIJ599lkhhBADBw4UGzdurFPHV2hW95C2BCFs7EUXXSTy8vJEhw4dREVFhd9yOPF2v6qrq8XQoUOb5DviDbvdLuLi4lzteD6Tp59+Wtx7772u/QceeEDMnj27Vp3U1FTxxx9/CCGEeO2118SUKVOaLE+wGD16dC3/rdWrV/u8T0VFRaJv374uH59oCF9bHzQh/Gu0b7RAO/S4bnFlwxcMF8MXDDdFWhb5nGSfm7K9PG3xwZenLW4bkvvaiAeQCqyL0MMX0UR1dbV49tlnRa9evcTChQtDGoM8VHH033vvPTFx4kRhMBjEnDlzvNbZt2+fiI+Pb3J27nDSkFNxUykqKhKdO3f2OTiLdv73v/+JsWPHijPPPFMsWbLEax1vMfeb6hzdEOedd5745ptv/Ko7dOhQ8fvvvwshhDj77LPF999/X6eOt+/Hhg0bRGxsbFDlT09PF6tXrxannHKKWL58eYNyOGVxvjPuioaiKCI+Pl707NkzqJmrG3pmt912m3j11Vdd+++884649tprXfuHDh0S7du3F8ePHxdCCLFy5Upx8sknN1meYHD06FHRsWPHWsqbU87Kyso69UP13oaK2E7djifNXCiSZi40iSgYbARja4mDNdmnlrElxvcse3naYvHytMWmSMsS7u3laYu/enna4pGhaLtxlWF/ODuOmkPi9WhSJPbv3y+uvPJKMWrUKK8zpMEmVMmwfv/9d1e7Q4cO9TlA1mg0rsFcNBPsGWgnoRyYhCMJ1vvvvy+uueYaMWnSJPH+++97rdOzZ0+xbdu2WmWhUmDvvPNOMXfuXCGE7/4XFRWJoUOHCsBVfuONN4r58+fXac8zoo9T7h49eojBgwfXaiMQunbtKnbt2iWmTJkiXn/99QblcG7OpIee5TExMWLz5s1B/X439MzOOeecWkkmf/31V6HT6Vz7y5cvFyNHjnTtHzlyRHTo0EEcOXKkyTIFyooVK8SIESPqlA8cOFAUFxfXKY+GBIaNoSWuSMiEdM1ja4l9kn32vb08bbH55WmLd748bfFrL09b/Kj7FmjbfjtbK4oyGghrCB8RZXkkfv/9d0aNGkXfvn1ZtmwZycnJIb+mu1NlMOPo//Wvf3U5yxYVFfm0IW8OfhJHjx4lNjY2JLH+165d6/wHFFT/FyEEF154YVDt+L2xadMmBg0aRL9+/di6dWud48ePH2ffvn11bMqb4gPhD04/CVD9GBwOR53+Z2ZmupyAneXJycleHa7nzJnjiujjTFQGsHfvXjp06EBqairvvvtuQJGWjhw5wuHDh+nevTsnn3yyVz+JvLw8YmJiiI2NrVUuhODw4cMMHTq01v3UaDQMGDAgqN9vz0hObdu2ZefOna59Tx+Jtm3bsmbNGpdfx3fffefyjwBo3749qamp2Gy2JssUKCtWrPDq35Samsq6devqlA8YMMD1OZQ5PyS+2XJAlGEqVzCVmyIti0QicdERyAfaAgPdtgH1neQXHtraLV62acCTwHZgSoS0SBFp5s2bJ3r06CHee++9SIsSFPyduXvllVfELbfcEmbpGofJZBIXXnhhSGb3gzljfPjwYbFw4UIxffp0MWjQoKCatfjCaDSKF154Qfzzn/+sZR/vpLS0VPTv379OeahWS/773/+KcePGCSF8v4PeyufNm+c1l8T06dNddv7ezrv99tvF888/H5DMGzZsEMnJyUIIIb777jtxzjnn1Klz9OhR0aZNG1FVVeX1nQnH6pPnNe644w6RlJQkVq9eLcrLy0XHjh1FVVWVq76nOWDXrl1dq0XO9hITE0Ni1ufv/cjKyvK6AjR9+nTx4osv1inPyMgQffr0Cel9Dia0wBUJWuCsb0vsU5+4uLLCYRpROExjirQs8jk2783zJn/vZfsWNYb1RRF8+CJSHDp0SEyZMkWkpaWJwsLCiMkRbPwdIP/yyy8iPT09zNL5z/r160X37t1FaWlpSNp3DnhowATMF6WlpeKVV14RGRkZonPnzuKcc84RTz31lFizZk1YBpyZmZni3//+t3jvvffEpEmT6hz/8ccfxZgxYwK6RmPYs2eP6NKli6iurq41kHU6xwqhOvy6D3DT0tLE4sWLxfjx42u1dfjwYdGtWzfXs/d2P999912/k7z54r///a8488wzhRBC7Ny5UyQmJtZJ1FhcXCwGDhwohAiPyZq/vPvuu6Jr164ufwx3ebz5dXi+h06zPmfCwmD1yV+TQa1WK3777bc65c8995y48847a5Xt3r1bJCQkBCXIQriQikTz2GSfWsbW2vr88rTFHYJZz+d9jXRH/Xz4IhKsXbtW6HQ6MXnyZK+RWpoz/g52nNm1y8vLwyxhw1RXV4uLLrpIPP300yG/1qmnnipWrFjRYL3KykqxbNkykZ2dLXQ6nejevbu44YYbxPvvv19ngOPuhNu5c+c6TrnB8MlwOlIvWbLENRh255NPPhFXXHFFQNdoLH369BGbNm0SX3/9tYiLixOxsbEiLi5OrFu3TgghxP/93/+JxMTEWu9mUVGRSEpKqtXOO++8Iy6++GLXvrd3esuWLaJ79+61ZuIbi9PPxFN+d5YtW+ZaaYk2nL4inu+U+7vmTZloavZsf/BnRbS8vFzEx8e7nL/dyc/Pr/XshRDi+eefb3aZ6KUi0Tw2f/o0fMHwbsMXDP/38AXDDw1fMLx0+ILhf/VRb/jwBcO/Gb5g+O7hC4aLaO5TS9taW59fnrb4gJ/19gZynTh3MydFUdKBDCHE054mUIqi/B3IE0LYPY+1RD755BOmT5/OY489xrRp0wLOKB1tOOPQN0SbNm0YMWIEq1atQq/Xh0Ey//n4448pKytjxowZIb9W37592bZtm2u/uLiYzMxM1q5dy0knnYTBYGDlypV88803DBgwgIyMDCwWC2PHjq1jN+/E+QyOHDlCcnIyF154Ya08AsHwyXD6SOzfv5+ysrI6x73lkAg1Tj+Jn376ibvvvptnn32W0047jdLSUoYMGcKnn37KN99848odAKovhzOXRFyc+m9r3rx5TJ8+3VXH1zudmJhIYWEhw4cPb5K8nvdoxIgR/PHHHwwcONBVtnXrVvr379+k9kPNpk2bXJ/d36m8vDzXO+z0LXGvM2zYMBwOh8uXylsbTty/D8OGDSMvL69ev5SUlBTWr18P+PZl+OWXX1zJ8TwZOnSoy0fCee3CwkKSkpIoLi4OSfbxUBDbqVub5Ox8AcwuyckwRVoeSUCYgeOoGchHAvm6XN1qW5bN85/SCeAj4BXg83AKKGl1tDcbrG/7Ua/uP9lG4Ols/Qiw2Ufd0prjLZrjx48zY8YM/v73v7No0SIMBkOLUyIaSzQ6XB84cID77ruPV1991etAI9h4KhLuCc/Wrl3Lww8/zDnnnMNvv/3G77//zhNPPMEZZ5zhU4lwp0OHDsTExNRJRhaos+jBgwc5duwY3bt3p3///pSVlTlnZVxs3bo1YorERx99xKRJkwC48cYbeeedd/j+++/p1KlTHQfbtm3b0rt3b7ZsUeM9FBUVsWbNGi677LIGrzd+/HiWLFnSZHndlYTi4mJWrVpFZmZmrcRzkVDI/MWX47x7Uru0tLQ6dTyzZ9cXzMCX47wvZs+e7WrLl5P5ihUraimT7iQnJ1NWVsaxY8dc30WAzZs3hyRoQaioqth7oiQnQ5FKRPNGl6uLB64CHrZl2SpsWbZlwJfAjZ51bVm2tbYs2zyg4Zk8iSQwngCK/NhyArlInMf+6UCWj7qfA88EcrHGoihKJmoI2LCwefNmJk2aRM+ePfnll1/o2rVruC4d1YwdOzZqsjs7Zx/tdjsJCQlhG7x5KhLu0ZxAjewzbdq0Jre/a9euOmV9+vQJKIrP5s2bGTRoEIqi0KFDBzp27MiePXvo0aOHq05ZWVnYV5rS09N54403OHLkiGugePrpp/O3v/2Nt99+m969e7Nx48Y6s8rJycls3LiR5ORk5s+fz+TJk2nXrl2D19NqtTzwwAPcc889fs2We1JWVsapp54KqAPmvXv3IoRwDZgLCgqiekXCfeXB2X9/6riv8Lh/77p27epqo7i4mIyMDBwOh6stf1bSnFnKp06dyquvvuq1zsqVK7nqqqu8HmvTpg1JSUkUFRWFLLKaROIkrltcP12uzn0WZrYty2Zy208FqmxZNvdQYqtRk/pKJBHBaNHPDsd1PFckugFV3ioC1UBYR9YijOFfv/nmG0aPHs2VV17J559/LpUIN5wrEp6z2eFGCMEFF1yA3W5HCMGBAwfCNvvoqUi4z8gGI8yk56xx7969ufLKKwMy0XCaNTnp379/nRCwkZhJT0xMZNWqVWzatInhw4dTXFzMLbfcQnV1NUIIdu7c6fW5OkPAVlVVsWDBAm699Va/rvf6669z4MCBJofZdb9HvgatkVjZ8Rf3lYeCggKv71RDdZzHi4uLEUK4lNGLLrqolhIB/n0fNm3axIgRI+qswrnjK/SrE2cI2KSkpEZdWyJpLJV7K8tsWTbFbTN5VOkElHuUlQOdwyKgxC8S43v2MxuswmywmiItS0vCU5HYCJzho+4ZQElIpYkAVVVVmEwmbrnlFj788EP+/ve/14rFLoGkpCSqq6tdZiXhZu3atcyaNYvU1FQ2btwYkdlHT0UiLy/Plb8gGPk9PPMJvPDCC6xatSqgNj0VCW+5JCKhSDzyiGoh6T6r7/4cfT1XpyLxzTff0K9fP799HgL1O3FfbfBlJlRWVha1KxLBJDk5mdNPP53U1FTi4uIoKiqqU2fw4MENfh82bdrEueee61OR2L59OxUVFQwZMsRnG05FIj09nV69egU9146kaQzoovTDlCAwJZgiLUsYqQC6eJR1AQ5GQJagYddou9k12n/bNdpDdo221K7R/rWeuvfaNdrtdo223K7RvmXXaNu5HXvHrtFus2u0B+wa7Tq7RntbeHpQm/2HdpUZLXrFaNGbInH9lorniPkN4E1FUU5zL1QU5VTUDNOvhUuwcLBr1y4mTJjAkiVL+OWXXxg/Xq5CemPjxo0cPnyY5OTkWnbhoWTr1q0uJ9xzzz2XAwcO8N5773m15Q4HnopEUlISiqJQUVHhc5a3MXjOCGdkZGCz2Thx4kST2/S2IuHpcB0Jkxxnsjmo7djb0HMdPHgwJSUlzJs3z+/VCFAH/572/cXFxaSnp7uSsfl6p4UQtZQtp8IHMGTIENegNZpXJIKN3W5nx44dtRy0Qb238fHxPPfccw1+HzZt2sQ555zD5s2bXWZO7qxcuZLRo0fX65+WmprKV199xU8//cT69evrXXGRhI9WmpBuHRCny9UNdSsbQfP3g3B3IJ8MvGrXaNM9K9k12ouBbOB8IBlIAdzNap4EkrUOexfgMuBxu0Z7mmc7kuZJLUVCCPEi8DWwXFGUjYqi/KgoykZgObBICPFSJIQMBT/++COnnnoqY8aM4dtvv6VPnz6RFilqyczM5NChQ1RXVzdoGuLvAM1bfY1Gw5w5czjvvPPQ6XQUFhby9NNPs2XLFp577jlGjx4dskzfDeGpSGzbto1u3brRvn37kFyvU6dOJCcn+xVZyxcNmTZVVFRw/PhxEhMTAxG10XhTGvx5rsnJyaxatYrFixdz3XXX+X29vLw8Bg8eDOBq291Zvr53ury8nDZt2tCpUyfgT4Xv/PPP56WXXiIlJcWlbLSGFQmA0tLSOmXO53bFFVewcePGBtvYtGkTQ4cOpVevXl5XOutztAb1/8aTTz7pcqLfvXt3I3oQPTijNiVn55siLYuk6diybIeAz4BHdbm6eF2u7kzgcuBfnnV1uTpFl6trj5phGF2urr0uV9ews1eYsWu0LgdyrcNeoXXYfTqQo/rWztM67AVah30f8BgwxXmwpvxYza4z7LPv5UZJ0DEbrN1D1banszVCiLsVRXkRVbPsDuwBFgshNnjWbQqKovwAjAOc01BbhRBhmVZ2Ogw6HA4UReG1115r1Mxma8XTLrywsJCOHTvW2Tp06MDq1as5dOgQQC1nVF9MnDgRh8OBEIK1a9eSk5NDbm4ul1xyiddBur9ha4NNnz592LlzJ9XV1cTExFBaWkpycnJIrzlq1ChWrVrFyJEjGx1eE7wrEr/++qtrf9u2bfTr1y/sUckacuz1RUxMDIWFhSiKwhlnnOG303RKSgobNmxAq9Xy1ltvkZKS4reDrq+VBmdo1Isvvpjy8nJiY2Pp3Ll1mEO7h4V1Rl1yPrtnnnmGkpKSes8XQrjezZSUFIqKiup8l1auXMkdd9zhs43MzEw2b1YDDO7du7fB/zPRSk3UpraRlkMSFO4A3gJ2oo6bptuybAW6XN0goBBIs2XZNgFJqGbkTo6gRsVMDqewfeLi+tk12loO5FqH3eS2nwpUaR12fxzI04EvPOr1tmu03bUO+x4Au0b7Cqpy0QH4Dfgq4E5IGsNms8H6Lapy+6XRoj8erIbrKBIANUpDHcVBUZQeQohgTP3cKYR4MwjtNApnrHFQByVz586VioQfeBs4rFixgsOHD3PkyBEOHz7s2tzNw6qrq7Hb7fz++++MHDmyVpsHDx7kzTffdIVtdHL48GGuvPLKcHSrUbRt25aEhAR27dpF7969KSkpCZsicdttt7lm0IUQ2O12tFotVVVV9SoVpaWldRQJ95n+SIUsbaoy6MwZ4RkxyR8UReHaa6/lww8/ZOzYsfTv39+VX6E+Ezlf98j5naivTkulvihQgwcP5n//+1+95+/fv5+YmBgSEhJISUmhuLiY888/33VcCOEybfLF2rVrXfktZKQmSTRgy7LtBa7wUr4J1RnbuV8CRDym/PbKyjKtw17fP67GOJB71nV+7oyqVKF12O+wa7R3oUYHPRc4hiScJAHXAzOB180G6yfA20aLflmgDXsmpNsrhOjmtr9YCHG+W5Vi6joUNRv8ceiU1MXbwCE+Pp74+Pg6dTUaTS2lo2fPnlx66aUMHz6coqIiSktLSUxMpLKykosvvpjBgwdTWlrqqh/NEVec5k3hVCQWLFgA1F4VEkJw/Lg6meBrQF1VVUVZWRkDBgxwlXk6W0dzyFJvOBOQQdO+v9deey0XXXQRzzzzDL169aKqqoqtW7fW8nXwxNc90mg0fPHFF/XWaanUpwgOHjy4QdMm95WyIUOG1DF/LC4uJj4+nr59+/psw3NyI5r/b0i8YEroBswDLgJ2A//AVP6el3pZwN3AUOAA8B7wAKbyuo41kmDTGAdyz7rOz7Xqah32KmCZXaO9AZgOvBgcUYOH2WCt824aLfq676Za917UgXkH4FNgutGiP9ZQO2aDdRyq+ddpqJFSfwDuNlr02zyvESyMFv0u1Pv9otlgHYZqovYvs8EqgHeAeUaLvq7dqh94Olt7ZvY6xWM/WFr0k4qi7FYU5X+KopzrrYKiKFMVRVmlKEpgoWvc8MehU1IXf8JHOvG0df/xxx8pKirijz/+oLi4mKqqKvbu3UuvXr348MMP+e677yLi89AU3P0kSkpKaoWdDAUjRozAbrdz9OhRn34MvrIMa7Vajh8/zqmnnuoaqHk6Wze3mfRAv79paWl07dqVZ555hj179lBSUsL48eN59dVXfb7TckWicThzfdQXKtp9pcxp2uROQ2Ffoe7/mWj+vyHxSh0nXkwJdZx4gY7ADKAHMBbV5Pr+MMnY2lkHxNk1Wn8cyAtqjrnX2+E0a/JCHNHrI1Hn3TQbrHXeTbPB2pCDeX3tdEUNYJSMulJwEJgf/K74pE/N1gU1IV1/4DezwZrdlMY8TZsaShQQjEQCM1HtBY8D1wF5iqKMFELU+jURQryOeqNRFCUoCQz8ScwkCQxfs5XuzpBCCNfgNlI+D03BXZEoLS31mSwrWHTs2JGhQ4eyePFiKisrGTp0KMXFxcTGxnLixAmEELWiEDnf7ZiYGFe0J/cVi169erF//36OHTtGu3btmt0AONDvb3FxMTt27GDmzJn07duXTZs2cdJJJ7Fhw4ZapjXulJWVkZqaWqd84MCB7Nu3j4MHD7a6FYn6cObf2bdvH926dfNax31Fwmna5E5DjtbO85rL/w2JB6YEpxPvcEzlFcAyTAlOJ97aAxlTuXu2wq2YEt4FzguXqE7cEtJ5JqJrsWgd9kN2jfYz4NGacK0jUR3IvaUIeBtYYNdo3wW2AQ8BCwDsGm0vQA8sRPUHuQDVxMZnKNlIYTZYXe+m0aKvAJaZDVbv72aNg7nRoi+oOfcx4F0gu6F2jBb91x7XfRlYEsKuUaPE3ICq1FQAucDJRot+q5v8f9CELNdefSRCiRBiudturqIo1wOXAiGPCCV/fCJHSzBF8FyRCLVpU3FxMZs2bWLixIl07dqVRYsWuQZeTp+Jzp0714pCJISoFZbTfcUiJiaGPn36sG3bNqqrq5k3bx4VFRXMnz+/0dmeI0Gg39/MzEz27FEnyHbs2EFmZiY33XQT69ev93nO1q1bOe+8uuOWmJgYVx6DrVu3Nsv3ORQoiuIyb6pPkXCu5g0ZMqTOisTKlSt59NFHQy5rNOCM2gTMLsnJMEVanjCRClRhKm9KFuhziEBI1ZqEdM1n1iV41HEg1zrsBXaN1uVArnXYN2kd9kV2jfZp4Hv+NPOZVdOGQDVjsqBawZQCM7QO+xeEGWdCOrei2R45JVKBKqNF32QH85roSIMa0Q6E573+L/A+cLXRol/hedBo0ZeYDdbnm9KwpyLRXlGUt9324z32QxGiTBAFjkeS0NISVoP69u3LunXrqK6urhMRKRRkZmZSXq76rJWXl7tWFpwD6oqKClJTU9Hr9V5DckJtE6Di4mJ27drFkCFDiIuLa9DPoqXhzUH3pJNO4qeffvJav7i4mEWLFvHll1+6zGfcla1hw4axdu1aysrK0Ov1YelDc8CpSJx2mvcw8Zs2bXIFX+jevTtVVVXs27ePffv2MXHiROx2O0ajkfz8/KhXbgOlJUZtciWk+5PZHjklmpYF2pRwMzAKiEgys9aI1mH36kCuddhrOZDXlM0F5nqpuwv/lMSQU5OQLhwO5n63YzZYTwYeQV3tCSV9G4rUZLToH2lKw56KxBMe+3Ma2G8UiqIkoto5LkEN/3otqiY2I5B2JdFPS1gN6tu3L0uWLGHHjh0kJCTQsWPHkF6voRClnTp1IjY2to4SoSgKbdq0qRXVCVTF5OjRowAuJcJX2y0Rb6tiTtMmb0ycONF1v7wpW872pGlTbRpyuHZXwhVFcflJZGVlud7DdevWtQrltiVSk5CuvsFa47NAmxKuQDW5uABTefNMGiJpDgTLwdyvdswG60moudvuMVr0S5sos18YLfrjZoO1NzAG1edIcTv2ViBtezpbrxNCzK5vC+RiqM7cjwO7UL3Y7wKuEEK0/FGMpNnjNG0Kh1kT+Odc7J4kz4lWq8Vut9dxjvelLDRXU7PG4s1B96STTqKoqMi1UuFOQ1He3FckmpOvSahpjCIBf0ZukiFdvZOcnX9ncnb+quTs/GPJ2fkLIi1PEFgHxGFK8C8LtClhAvAGkImp3BZ68SStmHVAnNlgbbKDudGi3+NPO2aDNQn4DnjMaNHXSVwYbMwG6xWojtWPAq+hjr9fw3uCwUbhqUi8FmiD9SGE2CWEGC2E6CyESBRCjBNCfBvKa0okwSLcioQ/kWk8lY20tDSfkbXc6yqKQtu2bVtV1Btv0cfi4+Pp2rVrrbC4TuLj413J+rwpWxqNhsLCQnbu3FlvqNLWRn2KxIkTJ9i5c2ctxcvp9+Meqri1KLd+UoY6ARfQrGHUYCp3ZYHGlBCPKcFnFmhMCXpUB9arMJXXseuWSIKJ0aJ3vZtmgzXebLD6fjdVB/NbzQZrmtlg7Yqbg3lD7ZgN1v6AFTAbLXpLiLvl5HHgZqNFfwpwqObvVOCXQBv2VCSkr4JE4gN3RSLUoV/Bv7C7jQmD6V7X16pFa2To0KF1HK63b98OqMqXr3ubmppKQUEBXbt2pU0bz8jZrZf6FImtW7fSp08f4uL+tKodMmSIK2P5wIEDW5Vy6w8lORmfleRkfE5NYq8Wwh2oTrk7UR1Ap2MqL8CUMAhTQgWmBOeS1cNAAvBVTXkFpoSvfbQpaQTOzNZ2jdYUaVmijDrvptGiLzAbrIPMBmuF2WAdBGC06BcBTgfz0pptVkPt1By7DTVc7KyaNivMBmtFiPs1yGjRf+xRlgvcFGjDinu8b0VRDgMZ1KNQCCGsgV7UXxRFyQQygdvri0sukYSLxMRELrroIs4991zuuOOOSIsjCQK33norY8eOZerUqa6yp556ivXr1/Pmm2/We27//v3p3bs3v/76a6jFbDYcOnSIHj16cOjQIdcKmJP//ve/PPDAAyxbpiZTLS4u5rzzzmPTpk107tyZ33//vVUptXGdu58YYHzbXQv1Gb0pOTv/cWBASU7GlHDI1lQURVklhBgVaTmCiexTy6A19tmJ2WDdAJxptOh3mA3W31AVnd3Az0aLvnsgbXs6W7dDzcTnS5EQqFpUWBBC5KHmmbg9XNeUSOqjb9++/PTTT0yZMiXSokiChKfDtRCCN998k3/9q36z1eLiYvbv38+2bdtIT09vFiF0w0F8fDydO3dm+/btHD16tFa0tltuuaWWf0RmZiZbtmwBVAWktTlYt8SoTRKJJCp5AzgLNTTvc6grKdXAs4E27GnadEgIkSKEGOxjk7+SklZN37592bJlS1h8JCThwd20qbi4mMGDB7NhwwZuvfXWOsnS3MnMzOTIkSMIIVxRnSQqTvOmzMxMHA4HVVVVOBwOnn766VqKhHSwlkgkktBjtOifMlr0n9Z8fhs1Z8ZpRov+4UDb9lQkJBJJPTidasPhIyEJD+4rEpmZma5wug0pBw2F523NOBUJT0Vh586d/POf/yQ9PZ3i4mK/IpNJJJHGmdlal6szRVoWSdNxJqQzG6ymSMsSaYwW/SajRW8PRluePhIHhRD1J4WJAIqiCOkjIYkG7r//fnJzc9m1a1ekRZEEiYqKCnr16kVFRQVt2rSpFQo2NjaWyspKr+elp6fXykuh0WhalVmOL4qLixkzZgx79+7F1/9t5/3ylqiyNZmHKYpyWAgRX1+d5Oz8OFQz5FnAAOB2oLIkJ8P7ixlhBibElG2+t3Nf6iaia7a0RNv6ltinhmhtfTYbrJtRXRLqxWjRB5Rdt5aPRDQqERJJtFBcXMw777zD7t27pU18C6JTp04kJCRQVFREbGwsQgiEEA3OkLeEbO2hIDMzs5YS0bZtW6qqqqiqqnLVca7gtIRElWHgIWpHg7kBmA2YIiJNA/iRkE4ikYSHG8JxEdeKhKIoPwghzm3wBEVZLIQ4P9SCeVxTrkhIIk56ejp2u901yJQz0C2D4uJihg8fzpEjR4iPj6dv375s3LixVc6QB4O4uLhaSoNzVUeu4NTFLWqTz2hNzY2WOOsr+9QyaI19DgfuKxJjFUW5mYZzSYTtIbiFf5VIIo60iW+ZOJ2mAY4cOULbtm19mjNJGmbYsGG1FAbnqo5cwamLjNokkUgawmyw/gv/TJR85oQwG6ztgEeA64HuRos+wWywXgSkGi36lwORz12RWI5/iSl+DuSCjUGGf5VEE74GSJLmjbtCKBXEwPGlMEgzJolEImkSG9w+9wCygDzUJHiDUCfccxto4zmgPzAZcCZ1LKgpD44i4Y9Zk0TSmpEzqi0TqSAGF6kwSCTBRZer8yvCpi3LVt1wLUlzw2jRz3Z+Nhus3wAZRot+qVvZWahZ4OvjSuAko0V/yGywVte0u9VssPYPVD7PhHQSicQHcoDUMpEKokTSgjAl+OfUZCr3nSQm+qjED9MWINbfBvvExfWza7QCmK112E1NFSwU2DVavT/1tA67NdSyRCHjqGsZtBw4vYHzjuMx5jcbrD2BPYEKJBUJiUTSqpEKoqSlk5yd79eMdklORkuY0d6AOuhWqD349tz3e9AdBQx2+5wBXA08iWrakgTMRM1Y7DfbKyvLtA57tEbXmuex3x/12e0BuqM+yy1Aa4yE8Rswx2ywPmK06I+YDdYOqFHcfm/gvI+BXLPBei+A2WDtCzwPfBCoQC0iId0//vEPnn/+eb/q/uUvf2HRokWhFaiVsnTp0nrNQqZMmcJDDz0URokkEokkeont1K1Ncna+SM7ON4X4UpXACT+25o+pPAZTeSym8hjgNtSBkgZoX/P3PeDWxjQZ6YR0tixbqXMD7gP+YsuyfWvLsq2zZdm+Ba4B7o+EbKFA67APdm7AG8BLQNcaxacr8GJNeaNoIQnppgBnAuVmg3UHUA6cheo3UR8PACWADUgE1gNlqEpIQDR7RWLXrl28/fbbTJs2za/62dnZPPjggyGWKnro1KmTa4uJiaFDhw6u/XfffTeo1zr77LOD7qi6adOmWn1QFIX4+HjX/tKlS5kyZQpt27atVW/EiBGuNo4fP47JZGLo0KHEx8eTnJzMLbfcQklJSVBljRTJycm0bduW3bt31yofOXIkiqI0m34uW7aMM844g4SEBLp168aZZ57JypUr/XqHk5OT6dChA507dyYxMZEzzjgDi8VSK7lcMPnggw/QarXEx8czZMgQli5dGpS6TeXll19m1KhRtGvXjilTptRbt6SkhEsvvZSuXbvSp08f7rzzThmlKkLURG1SwhD6dTDq7G0KcBewBJgAaGv+fg/cGYwLDeii9MOUIDAlmILRXoA8BtyGqXw9pvLjmMrXA9OAxxvTSOXeyjJblk2xZdlMoRCykSQAHT3KOtaUt0TuBbK1DvthgJq//0BVqBrF/kO7yowWvWK06E1BlTCMGC36EqNFfwYwBLgM1e/hDKNFv7GB844bLfoZRou+E9Ab6Gy06O81WvTHA5WplmmToihnApcJIWZ6VlQUJQf4XAgRtqhN/rBgwQIuvfRSOnTo4Ff9MWPGcODAAVatWsWoUS0/nHBFRYXrc3JyMm+++SYXXHBBBCVqHIMGDarVB0VRWL16NSeddJKrbN68efz973/n8ce9/zZcffXVbNmyhffee49TTjmFQ4cO8c4777B48WJuvbVRE1NRy+DBg3n//fe56667ALDZbK6QppHEZDLV+uuLAwcOMHHiRF599VUmTZrE8ePHWbp0Ke3atfP7Hc7Ly+OCCy6gvLycJUuWcM8997B8+XLmz58fzC7x7bffMnPmTD788EPGjBnDtm3bglLXE3/vHUC/fv146KGH+Oabbxp87nfccQe9evVi27Zt7N+/nwsvvJBXXnmFu+++22/ZJM2LkpyMUufn5Oz8+4BRJTkZ+2uK1iVn568CVgGvBnqtKEtIFwMkA3a3siSal1mTJ7nAd7pc3fPAZmAgcDcNR+1prhwCxgD/cysbDRyOjDhRwzFgFxBnNlhTAIwWfS2/H2e5DzqbDaqLied5jcVzReIB4L8+6v4AhHUqX1GUTEVRXq+vztdff8348eNd+/v27WPixIn07NmTrl27MnHiRLZs2VLrnHPPPZf8/PzQCN1MOHbsGDNmzKBfv37069ePGTNmcOzYMQB++OEHBgwYwJw5c+jRowfJycm1Vi+++uor0tLS6Ny5M/379+eZZ56pdZ6T3377jVNPPZXOnTtz7bXXcvTo0VoyLFy4kJEjR7pmkP/444+g9/O7777j22+/5YsvvmD06NHExcWRkJCA0WhsMUoEwI033sjbb7/t2s/NzeWmm2pHcy4rK+Oqq66iZ8+eDB48mBdffNF1LCcnhyFDhtC5c2fS0tL497//Xevc5ORknnnmGU4++WQSEhK8Ps9AWLduHQDXX389sbGxdOjQgYsuuoiTTz650W0lJCRw2WWX8eGHH5Kbm8uaNWuCJifArFmzeOSRRxg3bhwxMTH079+f/v29B75oTN1A+Mtf/sIVV1xB9+7dG6y7ceNGJk2aRPv27enTpw8TJkyQPiKti9Y0o/0cYMWUMAdTwnRMCXOAxTXlzZW/o5r2XAvMBa5DDd/590gKFUIeBhbZNdr37BrtU3aN9j1gEWrG91aH2WCdYDZYtwLbUf2BnNt6L9Wd5c6/6z3qr/dxXqPwVCRGoj4gb3wHnBboBRuDECJPCDG1vjo2m62WXX51dTU333wzpaWlbNq0iQ4dOnDnnbVXbLVaLatXrw6N0M2EJ554gp9//pnff/+d1atXs2LFiloz+tu3b2f37t1s3bqV3Nxcpk6d6jJbuvXWW3nttdc4ePAga9asQa+vG2Dh+PHjXHHFFdx4443s3buXa665hk8//dMX7Ndff+WWW27htddeY8+ePUybNo3LLrvMpcwEi++++44xY8YwcODAoLYbbYwbN44DBw5gt9upqqriww8/5IYbbnAdr66uJjMzkxEjRrB161YWL17M888/zzfffAPgMrkpLy9n1qxZ3HDDDXVmzz/66CMWLVrExo0b+eOPP1iwYEHQ5E9NTSU2NpasrCy+/vpr9u3bF3CbY8aMYcCAAT5NiSZOnEhiYqLXbeLEiV7PqaqqYtWqVezatYuTTjqJAQMGcOedd3pdBWhM3XByzz338MEHH3D48GG2bt3K119/zYQJEyIqkySs5ALfJWfnT03Ozr8kOTt/KvANLXFG21T+T+BmVFOOy4A+wC2Yyp+OqFwBYMuyVduybBZblu18W5ZNa8uy6Wv2qxo+u/mhddj/BYxFXVXqAjiAcTXlrREzqslevNGij3Hb6qyyOcuNFr27v9AwAvAX8oZn1KYuQFvA2y9dG6BzoBcMNvv376dz5z/F6t69O1dddZVr/8EHH+S8886rdU7nzp3Zv39/uESMSt59911eeuklevXqBagzp9OmTeOxxx5z1Xnsscdo164d48ePJyMjg48++oiHH36YNm3aUFhYyIgRI+jatStdu3at0/7PP//MiRMnmDFjBoqicPXVVzN37lzX8TfeeINp06YxduxYALKyspgzZw4///xzrRUmf3nmmWd4+eU/c6pcfvnl5ObmsmfPHvr27dvo9pojzlWJ8ePHo9Foas18r1y5kl27dvHII48AaqSi22+/nQ8++ICLL76Ya665xlX32muv5cknn2TFihVcfvnlrvK7776bfv1Ui4XMzEx+//13r3JMnDiRZcuWAbhWLZzBEM466ywWLlxY55wuXbqwbNkynnrqKW6//Xa2b9/OpZdeyhtvvEHv3r2bfE/69evH3r17vR7zJkdD7NixgxMnTvDJJ5+wdOlS2rRpw+WXX87jjz/OE0880eS6Tppy7xrL+PHjeeONN+jSpQtVVVVkZWVxxRVXBNyupNnwd9QZyWuBfsA21BntRjuvNgtM5YvwPUHaLNHl6i5Cnfjt5F5uy7I9EhGBQozWYS8ECiMtR5TQFXjNaNH7Ew7YnceAoUaL3jm+X282WKcB64AFgQjkqUg4gIuAL7zUvajmeFTRtWtXDh486No/fPgw9957L4sWLXLNah48eJCqqipiY2Nd+4mJiZEQN2ooKysjKSnJtZ+UlERZWZlrv2vXrsTHx3s9/umnn/L444+TnZ3NySefTE5ODqefXjuEcVlZGf3790dRlFptOCktLSU3N5eXXnrJVXb8+PFaMjSG+++/36uPRPfu3V1mMy2dG2+8kXPOOYeNGzfWMWsqLS2lrKys1ntfVVXF2WefDcDbb7/N3LlzXY7ZFRUVdZy3+/Tp4/rcsWNHn8/KfbDbGDt/rVbrWuVwOBzccMMNzJgxg/fff7/Bc32xdetWunXr1uTzPXH6Yt11110uBfW+++7zqhw0pq6Tpt47f6murubiiy9m2rRp/Pjjj1RUVHDLLbcwc+ZMnn662U7SNlucUZuA2WFwuAZcIV4tNVvLxpTQFjXKzUg8Bt2Yym+qe0L0o8vVvQxMQnWQd/cTaOzAsllg12i7oUakGonHM9Q67OdEQqYIMw91le2tRp4XMn8hT0XiOeA1RVFiUR2rqxVFiQGuQF1OabSXfKg5+eSTWbduHaNHjwbg2WefZe3atSxfvpw+ffrw+++/c8oppyDEn98xu91eK6pPa6Rfv36UlpaSnp4OqNGRnLPNoPqaHDp0yKVMbNq0ieHDhwMwevRovvjiC06cOMHLL7/MpEmT2Lx5c632+/bty9atWxFCuJSJTZs2MWTIEAAGDhzIgw8+GPIIWhdccAEvvPACW7ZsqeW/0RJJSkpi8ODBfPXVV8ybVzsM98CBAxk8eDDr19c1hywtLeX2229n8eLFnH766cTGxjJy5Mha35lwo9FomDJlCq+99lqT21i5ciVbt27lrLPO8nr8kksu8Wn2dPbZZ/P111/XKe/atSsDBgyopSD7ojF1w8XevXvZvHkzd955J+3ataNdu3bcfPPNPPTQQ1KRiAA1UZvahvu6ydn5Xme0S3IyWtqMdi4wAsgDdkRYlmBxPTDSlmXb3GDNlsF7QDvgI6SDNagJ6e42G6zZqH4SLowWfX2K1XOA1WywzudPJ/0pBMFfqJYiIYR4T1GUPqhfvnaKouwGegBHgVlCiKZPDYaISy+9lCVLljB58mRAXW3o0KEDiYmJ7N27l9mz64bIXbJkCe+88064RY0qrr/+eh5//HFGjx6Noig8+uijtWzqQTV3mjNnDsuXL2fhwoXMnj2b48eP8/HHHzNx4kQSEhLo0qWLa6XHndNPP524uDhefPFFjEYjX375JStWrHCZmd1+++1ceeWVXHDBBYwZM4bDhw/zww8/cM4559QyVQuUCy64gAsvvJArr7wSi8XCiBEjOHLkCO+++y5t27bllltuCdq1ooF58+axb98+4uPja4X0HDNmDF26dOGpp57i7rvvpm3bttjtdo4cOUJ8fDyKotCzZ08A5s+fHzQHZX9n0x0OB/n5+Vx77bUMGDCAzZs38/777zNu3LhGX/PAgQP897//5Z577uGGG25Ap9N5redNUfCHm2++mZdeeokJEybQpk0bnn/+eZ8+FY2p60ljViIqKyuprKykqqqKqqoqjh49SlxcHHFxteeKevToweDBg3n11Ve5//77qaioIDc3t9VPrLQmkrPzW9OM9gRgMKby/ZEWJIjsAfYH2kg0Z7b24Aygp9ZhD64DZfPlzZqtURgt+n+aDVYbas6RU1BNGm8xWvSBm/0JIepsqL4SFwN/rfnbxVu9cG2qmN7ZtWuX6N+/vzh8+LAQQoitW7eK8ePHi/j4eDF06FBhsVgEIE6cOCGEEGLFihVi5MiRPttrySQlJYlvv/1WCCHEkSNHxF133SX69Okj+vTpI+666y5x5MgRIYQQ33//vejfv794/PHHRffu3cXAgQPF22+/LYQQ4tixY+Liiy8WiYmJonPnzmLUqFFi6dKltc5zsnLlSjFy5EjRqVMnMWnSJDFp0iTx4IMPuo5//fXXYtSoUSIhIUH06dNHXH311eLAgQP19gEQ69evr1WWlZUl2rRpI+Lj411b9+7dXcePHTsmHnnkETFkyBDRsWNHMWjQIHHrrbeK0tLSAO5m9OD+XN05ceKEAMTGjRuFEOp347rrrhO9e/cWiYmJYuzYsa7zHnjgAdG1a1fRvXt3ce+994pzzjlHvPHGGz6vMWvWLDF58mSv8kyYMKHWs3DfJkyY4PWcLVu2iGuuuUb069dPdOzYUfTr109MnTpVlJeX+9XXpKQk0b59e9GpUyfRpUsXMW7cOPHyyy+LysrK+m9eEzh+/LiYPn26SEhIEL1796713ZkwYYJ44okn/KrrjabcOyHU54E6EHRts2bN8irTb7/9JsaPHy8SExNF9+7dxdVXXy127NgR4F2RNAXgkAjz72nSzIV7kmYuHBiq9oFV4e6Tz21Wl9ViVpfegbYT1y2ubPiC4WL4guGmSPdp+ILh04YvGP7l8AXDTx++YHiK+9Zsn1M9W+EwzbLCYZohwWgrMb5n2cvTFouXpy2O+HNsSZsixJ+TEIqi9BFC1FoqcUdRlNOEEL8ErL00EkVRhLucnjzwwAP06tWLGTNmNNjWVVddxa233sqll14aRAlbFj/88AM33HBDnbC5EolEIgkeiqIcFkLEN1wzeCRn568DTivJyTjYYOUmMDAhpmzzvZ37ArMxlZtCcQ2/MSX8DXUG9gU8TZtM5VZ/m1EUZZUQIioST+lydb6ybApbls1ve/do6lN92DXaR1HNuebjYcqjddgb5SfQXPrcEGaDtTdqbo0egMt21mjRN9ZvIih4+kisQ12NAEBRlPVCiKFux793Px5qFEXJBDIbqjdnzhy/23QPQSqRSCQSSSvjWeDd5Oz8J/EYXJfkZASUmAqiLiGdM/a75yBBoGb5bnbYsmyeYftbOmcDW4ALPcoFjXc4bvaYDdYrgHdQ8z+kAwXAcGAZEbofnoqEp1dgjwaOhxQhRB6QpyjK7eG8rkQikUgkoSYSUZv4M3u1p7OOoHlnfK6LqXxwpEWQBIbWYT+v4VqtiseBm40W/cdmg3Wf0aI/xWyw3oyqVEQET0XC036ooX1JC+Tcc8+VZk0SiUQSYiIRtakkJ6O1zWi3KHS5ujjgDmA8HqYttixbiw6HatdoFdz6q3XYfZl5tWQGGS36jz3KclHNvu5v6GSzwRoD9DZa9NsaqusvnoqERCKRSCQSSfPHlNAFMOFl0I2pfFBkhAqY5wA98DrwBPAgMB01a3GLw67R9kdNmHgOkOhxuGWtoPnHTrPB2tto0e8ASswG6+nAbhq4F2aDNRF4BbgaOAHEmw3Wy4AxRov+oUAE8lQkOiqK8l+3/c5u+wrQIZCLASiK0g01ocZFqJ3/hxDivUDblUgkEomkpZGcnV/nN7MkJ6PJv5nJ2fk+Z7RLcjLCM6NtSqjTJ0zl3vtkSrgXmIk6/vgUmI6p3N9QoK8AA4BHUe3KbwD+r6adoKLL1dXpky3L5rVPulxdnT7Zsmz+9ukvwOm2LNsmXa5uti3L9oIuV/cN8Bqq0hQ0apLB1eqT1mH32ie7RlunT86QrY1pxwsW1DDF5wNLUBUKE/BV03pVP2aDtY6sRoveq6xmg7VOn40W/TF/2jEbrOej5mcbBCwHphgt+lI/RHwDOKvmes+h+i5Xo/o+1YcF2IeahM6ZJfynmvMCUiQ8lzhvRe24c7vN7fObNfuBYgaOA72BycCriqJEzLarOfD6669HWoQWibyvwUfe09Ag72vwaUb3tM5vZnJ2fiC/mc8B04D/AqehDkh6AX5HMQoCdfqEKaFun0wJFwPZqIPIZFQH6brJoXxzEXAVpvIvgKqav9cCNwYivA/q9EmXq6vTJ12uLtA+dURNKAZwRJer62jLsjlQcwMEmzp9smu0dfpk12gb6pNf7fjgDOAWrcP+OyC0Dvtq1LHq35rQH3+oI6vZYK0jq9lgbXSfne2YDdYewGfAw0A3YBXwoT/CGS36p4wW/ac1n98GUoHTjBb9ww2cej5wd41Jk6g5fxfqdz8gPBWJvkKI3Pq2QC6mKEo8cBXwsBCiQgixDPiSAL7UeXl5QT8WqnabeqyhHzxf5zanPkbiWdV3X6Opj83pPY6mexqJa4aqH025r63l3rSUd9Ubydn5rt/MkpyMipKcjIB/M1FntC8pycl4Aais+XsFEB6nVlOCq0+YyiswldfXpyxgHqbyAkzl+4DHULPx+ksMUF7zuQJTQiJqIq6Tmii9V3S5OlefbFm2CluWrcE+2bJsBbYsW1P6ZAdG13xeBZh0ubqHgK1Nld/rRTRaV5+0DnuF1mFvsE9ah71A67DX6lMj2/FGFeDMrrrfrtH2BA4B/ZvWM9+YDVaXrEaLvsJo0TfYZ6NFX2C06Gv12Y92/gIUGC36j40W/VHUFZYRZoNV46ecsWaD9UyzwXoNaobqdX6cVo5HACWzwToI9fsQEJ6KxAOBNtgAqUCVEMK906sJwNs82n6cwvkj09C5zamPkXhWoWgzmo5F6prRJEs0PQ/5/Y+uZxyKNkMgSypQVZKTEbTfTDxmtJOz8zuW5GSEakbbG6moqwP+9Cm95ph7vd6YErr7ea3VqCZcAEtRZ4lfxb+BV2NIBapsWbYm90mXq/O3T/fw58D6PuBU1DD5UxslccOkAlVah73JfbJrtN0b2Y43lgPOxF/foM7cf4aqRAWbVKDKaNE3uc9mg9XV53raqXWu0aI/BBT5uE4tzAbryaihXz9GNdP7GFhvNlhHNHDqm8CnZoP1PCCmxrciF9XkKSA8E9IdFEJ0DrRRnxdTlLOBj4UQfdzKbgcmCyHO9ag7lT+/GKeFSiaJRCKRSCJF0syF7ru1wsAmZ+efDXxckpPRx63sdmBySU7GuU25XnJ2/o/AjJKcjBXJ2fl5qDPcB2ra1DalTXfSesZuKDR2GuJWVDsxnSnhbOBjTOV93MpuByZjKj+3VmOmhCLAiKl8Uc1+G1RzkcGYyksaFMaUkAIomMqLMCX0BHKATjUyFdZ/8p+0799+w9A5Q2v1yZZlc/VJl6s7G/jYlmXr41Z2OzDZlmWr1Sddrq4IMNqybItq9l19smXZGu5TkDipXbsNeYNTavVJ67CbnDt2jfZs4GOtw97Hrex2YLLWYT/XvS27RlsEGLUO+6Ka/T+fkzpj7lc73rBrtIlAjNZh32vXaDugRibqBDyvddgbNZvet2vShoeunV+rz0aL3tVns8F6NvCx0aLv41Z2OzDZaNHXktVssBYBRqNFv6hmv06ffbVjNljnAbuMFn222/H/AW8YLfoF9fXBbLCuAt4H5hotemE2WBXg3pq2fY6Va+rdgzquTgI2ofrVvGC06AOKyOrpbB2nKMrN1JMvQggRSMKLCuomtOsCHPRynddRoxJIJBKJRNIa8fs3sxHcg2ouAuqM9qtAZ4I0o124q6ohs6HG9MmzrvOzf/03lRe7fd6FalvfaI5uPRo9fQJ0uTo9arbnfkAZ8IEty7bY3/MBNhw7Fq4+BfQOax32/W6fj6CaEDWJbftKo6XPgdyTVOB55+C/Rpl4gQYc7WvqP1+zBRVPRaINcFM99QPNJLgOVVkZKoRYX1M2AjUzn0QikUgkkj9ZB8QlZ+cPLcnJCMpvZklOxkq3z+uBCwITsdGsA+IwJQzFVN5Qnwpqjn3kVm8HpvI9fl/NlHALHoNu4C1M5cHMi7UOiNPl6obasmxN6pMty+ZXn3S5uvtQnXznA7+hRv15T5ere9qWZWsock9jWAfE2TXaoVqHvUl90jrse+wa7dFGtFMHu0bbFjWqkOczfELrsB9tQr/qYx0QZzZYhxot+ib12WjR7zEbrEcbaKcA1ccCcPlUDPFxHU++Ai4D/u1Wlgnk13eS2WDV+zh0DNjiZ8Qor3iaNh0QQnhqSUFFUZQPUBWS24CRqDflDCGEVCYkEolEInEjOTvf629mSU5Gk38zk7Pz68xol+RkNGpGOyBMCV77hKm8wKPeBGABat6EbagRplZgKs/GH0wJTwOXo87ClqKadNwN5GEq/3vA/XBDl6vz2idblq3Ao57XPtmybH71SZer2wpcbMuyrXErSwe+tWXZ+gXekz+xa7Re+6R12As86nntk9Zhz25MOz5kmAcMQ82Z4XyG/wA2aB32WwLtoydmg9WrrEaLvsCjntc+O82V6mvHbLD2BDYAt6AqALOB8UaLfpwPmf7FnwmhO6AqEr+g+joNRDX//8Jo0U+qp18bUb/vAHsAp0/OTqAP8AdwnZvi4zeRyHB5B+qN2Ilq5zVdKhESiUQikXilzm9mgErEfagzuntRBzF7gPeSs/P/FgRZ/aVOn9TITAmDMCVUYEpQk8WpvhFPo8bKL63ZZjXiOlOA8zGVv4qp/CtM5a+ihoS9OWg9+ZM6fbJl2Qp0ubpBulxdhS5XNwigxjcikD6BOgh1p5g/B5rBpE6ftA57gV2jHWTXaCvsGu0ggBrfiPr65LUdP2W4Apioddi/1jrshVqH/euasisC7Jsv6shaM/gfZDZYK2oiHVHjG9GoPjuVkZqwq1ehKkf7gLHAdfXItAHVGbsIWAPMQXU8L6z5O4eGVzPmAS8CiUaLvh9qcr8XUJ2tE4GVqHlXGo3nisQaIcTwpjQkkUgkEokkuknOzt8KXFySk7HGrSwd+LYkJyOoM9oRR3XWPhVTeblbWSLwC6byIb5OizZ0uTr3Sd9bgXNRbeK3oM5IPwwssWXZ3gy7cCHGrtEWABdqHfYyt7L+wH+0DrvMQeYnZoN1F9DXaNFXupW1AcqMFn3PGvOqLUaLvmtj267lIyGVCIlEIpFIWjzhmtEOP2qkJifPA59hSsjhz0H3/6Em5WtOVPLn83EGw7neo+yvqCE+mz12jdbdnv9fwCK7RvsSfz5DI/B2JGSLJswGa77Ros/ws/oh1PwjP7mVnYaaNRzU7NhNwtPZWiKRSCQSSQsiOTvffUbbBMxLzs43UXtGu7HmNdHKBtQBtnv0Sc9ke3rg5bBJFDiDIy1AmJnnpcwzz9k04KkwyBLNnN2Iuo8A/zEbrF+i+lYMQHXSvqvm+PnAJ00RopZpk0QikUgkkpZFcnZ+NXVntPEoEyU5GbFhFUwSELpc3QBblm1LpOWQRAazwXrQaNH7nfvNbLCmofpm9EN1EP/EaNH7nU/FF3JFQiKRSCSSlk1rm9GuiynhTEzl/4u0GEGmkLr5CFosdo32eq3D/n6k5YgipjWmco3SELDi4IlLkVAUJaW+ik6EEMUN1wouiqLIZROJRCKRtDiEED4TwAaLkpwMrzHik7PzB5TkZAR1RltRlEVCiAnBbDNIfE0TB91R3Kcm00z79BpqBKQm0Uz77BOjRf+ev3XNBms31KzgI1Ezg7u3c04gcrjbTW4A1tf89bU1Or5ssBBCeN1uv/32oB8LVbtNPXbaaac1Sdbm1MdIPKv67ms09bE5vcfRdE+j7XkE0o+m3NfWcm+a87vKn46OkSLos5NAjxC0GQwCUdhkn6KDgJTuxPieJ5sNVmE2WE1BkifsmA3WLmaDda7ZYP3FbLCWmg3WTc6tgVPfA04H8lB9UNy3gHCtSAghIpFTImAyMzODfixU7QYiT334Orc59TESzyoUbUbTsUhdM5pkiabnIb//0fWMQ9FmqJ6xJCg0OXNvNKDL1Z0JXGbLss10ltmybJ1rjuUAn9uybD9HSr5gY9do/6J12D/zKF5ac6wN8IjWYX+4MW3uP7SrrCaHQnPmFVRH6UeBd4AbUCORfdrAeWcAPY0W/bFgC9QsnK0VRRHNQc5QMWrUKFatWhVpMVoc8r4GH3lPQ4O8r8EnGu6poiiHhRDxkbp+cnb+wZKcjM7J2fmxwKySnIxHAm1TUZRVQohRQRCv6ZgSTgOOYSpfU7PfEzX5Vjpq+Mv7MZVX+NtcNPRJl6vLB16xZdnyvRybABhtWTa/tdVo6FN92DXaUmAFYNQ67Dvdys9CDXO7WeuwX9iYNqO9z/5gNlh3AlqjRb/HbLDuN1r0iWaDtT+QZ7ToT63nvGVAltGiLwq2TF6drRVFiUPNyjcedfnLtZwkhAjIlkrSeKZOnRppEVok8r4GH3lPQ4O8r8FH3lMoyclwRnyJAx5EDRHZEngemI2aBRjUgWdf4HXU/AtPo45xgoYuV9cN1UzkImA38A9blq2ODbsuV6cAj6Fm1+4E/IaqBDSUmXgksMjHse+At5omedSSDvwTKLBrtH8D/o363K4BZmod9oBNcpopMYAzwWKF2WBNRI3ANLSB86zAIrPBOh/Y7n7AaNEH9O74Mmd6DtUb/L+oCSs+BXrVCBI2FEXJVBTl9XBeMxqRP3ihQd7X4CPvaWiQ9zX4yHtah6A4fQ/oovTDlCAwJZiC0V4T0VJjBlOTyfoS4AZM5WZURSIUdmZm4DjQG5gMvKrL1XnLvHwNcAtqDoBuqCsk//Kj/S5AWx/H2gB+hwFtDmgd9gqtwz4duBp4FnWw3BcY3oqVCIDVqJP8oL7jZuBVwNHAeWej5o25ELjRbbshUIF8hX/9C3C6EGKToiizhRAvKIryDarHvCnQi/qLECIPyFMU5fZwXTPUFBcXk5mZydq1axk2bBh5eXmkpPgVMEsikUgkklARFPvhLQdEGabySNuhx6EO6gHGAdsxla8DwFS+uUa5CBq6XF08anz+4bYsWwWwTJer+xJ1oJbtUX0wsMyWZSuuOfcd4F4/LuNAXe34wsuxi2h4IFmLPnFx/ewarQBmax12U2PODRd2jbY76qT2CeB3IA1IxWNGvTVgNljPNlr0S4Hb+VPpvxt4EkgE6nW2Nlr0nkkZg4YvRaIjauY7gCOKonQUQjgURTklVIK0FjIzM3E4HFRXV+NwOMjMzKSgoKEVTYlEIpFImk5ydr6+nsO+ZrqbKwWoM/8fAdehmv6omBL686dpiF/EdYvrp8vVuStas21ZNpPbfipQZcuyrXMrc585ducD4Fpdri4V2Ahk4dtkyZ3ngNd0ubpYVMfqal2uLga4AnVW+j5/+wOwvbKyTOuwR1rh84ldo70OeBFYCKRpHfb9do32BuATu0b7GfB3rcN+IKJChpcvzAbrJUaLfrmzwGjR7wJuMxusc1EVWb8wG6wKbiuQRou+OhDBfCkSdmA0qqPLKsCkKMoBYGsgF5PA2rVrqa5Wn1l1dTVr166NsEQSiUQiiQSxnbq1Sc7OF8DskpwMU4gv15A5SEPhI5sTM4E8TAkWoAo4y+3YtUCjEtNV7q0ss2XZ6ht0d6KuclKOd3OjbagmKWtrZNsM1KfkAWDLsr2ny9X1AXKBdrpc3W5UH9ajwCxblq2lJWp7Episddi/dRZoHfZ37Brtf4CXUUMXD4iUcBHgDmCh2WC9yGjR/+YsNBusZmAC3pVW3Or1R71v56CuYLgTUEZ7X4rEPagvOKha7quoXwhpVBogqamp2O12ABRFYdiwYRGWSCKRSCSRoKpi74mSnIywrAaU5GS0nuzWpvJlmBIGoa4UrMNUftDtaD7qqkAwqaBusrsuwEEvdWehTtQORDXRuQGw6nJ16bYsW715RWxZtrm6XN2bqPkAugN7gJ9sWbaWODM/XOuwH/IsrIngNMmu0V7e2AYT43v2MxusAphttOhNgYsYPowW/Qdmg7Ud8I3ZYNUbLfo1ZoN1Hqrvw3ijRd9QYkkLat6a84ElqAqFCfgqUNlaRPjXf/zjH/Tu3ZsZM2bwww8/cMMNN7Bli/d7et9995GamorBYAiVuPXy9ttvM23aNI4dO0bPnj356aefotpH4txzz+WGG27gtttuC0n7DT0vg8FA//79efjhRoWLdqEoCuvXr+ekk04KREyJRCIJOpEO/xoKGgyxaUpohxoL/wJUZ+MNwAOYyr8Oj4SNp6E+1fhI7APSbVm29TVlbwNltixbtkfdhcC3tizbC25l+4ELbFm2sMUjbgmhUBtLS+iz2WCdihqR7GdURfl8o0XfoM+I2WDdAwwyWvSH3MLGdgN+NFr0mkBk8hq1SVGUbEVRRnuUjVEU5e+BXCwU7Nq1yzU494f/+7//44knnuD48eMNVw4BixYt4umnn+add97h7LPPDosSsWzZMs444wwSEhLo1q0bZ555JitXrgz5dYOBxWJxKRE//PADAwYEZyUzPT2dTp060alTJ2JjY2nfvr1rf86cOSxYsIDY2FhXmXMrKytztfHee+8xatQoOnXqRN++fbnkkktYtmxZUOSTSCSSFkocqjnPeCABeBj4CFNCciSFCgRblu0Q8BnwqC5XF1+TPO5yvEdjWglco8vV9dbl6mJ0ubobUaMubQifxJLmhtlg1ZsNVj3qe/I/4DzUMMJpbsfqowqorPm832yw9gQOAf0Dla0+06aXPMoKgc9R4/hGDQsWLODSSy+lQ4cOftXv27cvGo2GL7/8kquvvjrE0tVm//795Ofn89JLL7Fv3z5mzpzZ8EkBcuDAASZOnMirr77KpEmTOH78OEuXLqVdu3Yhv3ZDVFZWNlwpRLg7uHtbdVmwYAGnn366T8Vg7ty55OTkYLFYuPjii2nbti2LFi3iiy++4KyzzvJ6jkQikbR6TOWHqB39cSGmhI2ooeZLIiFSkLgDNZfDTlSTo+m2LFuBLlc3CHX8lGbLsm0CnkINp/87EI86MLzKlmXbH05hm0PUJkktPH2c9qH6kTgRQH0z08uBS1HzcXwDfAgcQfWDDghfikRb1HBb7hwH2gd6wcagKEomDcR7/vrrr7nlllvqlM+ZM4e5c+fSqVMnnnjiCSZPnuw6du6555Kfnx92ReLDDz/kwgsvpHv37nTr1o0jR45QVlZGv36hC5ywbp0aROL6668HoEOHDlx00UWu4yaTiQ0bNvDOO+8AUFJSwuDBgzlx4gRxcerrUVRUxJgxY1i7di3nnnsu8+fPp1u3bl6vd+zYMWbOnMlHH30EwKRJk3jqqado166dy4zprrvu4rnnnuPCCy/k1ltvBXw/rylTpjBgwAD+8Y9/cMkll3Ds2DE6derk6tuWLVu45557sNvtdOjQgauuuoq5c+fStm3ozI7Ly8t55JFHmD9/Pn/5y19c5ZmZmWRmhiI8uUQikbRQTAm9UU00mnX4QluWbS9qBCXP8k2oztjO/aOAsWaLGNEetUlSG6NFH6iP0438aYU0A/gbqu/z8wG261OR+AVVu3a/gAH4NdALNgZ/8kjYbLY6Dsvbt29n9+7dbN26lZ9//plLL72UUaNGuepptVo+/fTTkMrujQULFvDQQw8Bqu3+mDFjWLFiBVdccUXIrpmamkpsbCxZWVlcd911jBs3jq5duzaqjbfffptvvvmGwYMHc9NNN3H33Xe7FA9PnnjiCX7++Wd+//13FEXh8ssv5/HHH+exxx4D1Gezd+9eSktLqa6uZvny5Q0+L4D4+Hi+/vrrOv4U27Zt47nnnmPUqFFs2bKFSy65hFdeeYUZM2Y0/mb5yU8//cTRo0e58sorQ3YNiUQiaY64EtL9yWxM5SavlU0JbYB3gVxM5Y3KgyCRhIIav4FaGcqNFn2dDOU1de9FjRDWATVx83SjRX+sxim6jh+Q0aKPiB+Q2WCNBV6gJmCS0aI/AjwerPZ9Zba+F/i7oii/KIrykaIov6LerLuDdeFgsX//fjp3rhth7bHHHqNdu3aMHz+ejIwM1ww5QOfOndm/f38YpQS73U5JSQkXX3yxq2zs2LEsX768nrMCp0uXLixbtgxFUbj99tvp2bMnl112GTt27PC7jRtvvJHhw4cTHx/PY489xkcffURVVZXXuu+++y6PPPIIvXr1omfPnsyaNYt//etPM9GYmBhmz55Nu3btapmj1fe86uO0005j3LhxxMXFkZyczLRp01iyZInffauPn3/+mcTERNc2ZMgQAPbs2UOPHj1cKzYSiUQSKpKz8+9Mzs5flZydfyw5O3+BW/m45Oz8b5Oz8/cmZ+fvSs7O/zg5O79vBEUFXAnpFLfN5LWiKSEG1YfgOHBnGEWUSOqjToZys8FaJ0O52WC9GDXZ4PlAMqpZ0eyaw179gMwGa3KIZfeK0aKvQlWMAsoX4QuvioQQogB1qfGfqI5BTwPDhBCFoRAiELp27crBgwfrlMXH/xkIIykpqZaT7MGDB0lMTAyXiIC6GnHjjTfWGnyOGTMm5IoEqCswCxYsYMuWLaxZs4aysrJGzdgPHDjQ9TkpKYkTJ06we/duDAZDLQdlgLKyMpKSkmrVd7/3PXv2pH372hZyDT2v+li3bh0TJ06kT58+dOnShQceeIDdu3f73bf6GDduHPv373dtRUVFAHTv3p3du3dH1MdDIpG0GspQZw/f8ijvCryOOohJQg01Oj+skjUVU4KCOuvbG7gKU7mnKbVEEnbMBqszQ/nDRou+wmjRLwOcGco9yQLmGS36AqNFvw/V8XkKgNGiP2S06E1Gi77EaNFXGy36hajJB08LS0e88xww22ywtgl2wz6nVIUQFQQ/1nLQOfnkk1m3bh2jR/8ZZGrfvn0cOnTINTjdtGkTw4cPdx232+2MGDEibDJWVlbyr3/9i++++65W+ZgxY1i1ahVVVVXExgaUD8RvNBoNU6ZM4bXXXgNUk6HDh/8MXb19e90oYps3b3Z93rRpE23atKFHjx5YLBYsFkutuv369aO0tJT09HRXfXcfEEVR8KSh51XfudOnT+eUU07h/fffp3Pnzjz//PN88skn9d6DQDn99NNp3749n3/+edj9bCQSSeuiJCfjM4Dk7PxRuCXgKsnJqGUmkZyd/zJqfPjmwKuAFrgAU/mRSAsjaR245ZFw4plPIhWoMlr0/mQoTwe+8KjX22ywdjda9HvcK5oN1mjwA7oL6APcZzZYd6E6ZwNgtOgHBdKwa0VCUZRFbp+XKoryX29bIBcLBZdeeqlXU5ZZs2a5IhQtXLiQa665xnVsyZIlXHLJJWGT8T//+Q8DBw4kLS2tVnn37t3p3bu3K0FdKHA4HDz77LMuv4LNmzfz/vvvM27cOABGjhzJf//7XzZt2kR5eTlPPvlknTbeeecdCgsLOXz4MI888ghXX321T8Xn+uuv5/HHH2fXrl3s3r2bRx99lBtuuKFBOet7Xk569+7Nnj17KC//M4HowYMH6dKlC506dcLhcPDqq6/6dV8CISEhgUcffRSj0cjnn3/O4cOHOXHiBF9//TV//3vURUiWSCRRijOztdtmCqC5c2gODsumhCRgGjAS2I4poaJmm1z/iRJJYOw/tKvMaNErbpvJo0pjMpR71nV+rlW3ZgXgXSDXaNFH0g/oBlSfjYtrPt/otgWE+4rE226f3wy04XBx0003MXLkSI4cOeKyue/Tpw9du3alX79+dOzYEYvFgkaj5tvYtm0bhYWFIXVw9mTBggXcfPPNXo+NHTuWFStWeJ2BDwadO3dm+fLlzJ07l/3795OYmMjEiRP55z//CcCFF17Itddey8knn0yPHj2YOXMmX375Za02brzxRqZMmYLD4WD8+PH1DtYfeughDhw4wMknnwzANddc43Iw90V9z8sdjUbD9ddfT0pKClVVVRQWFvLMM88wdepUnn76aU455RSuvfZarFZrY2+TV3766SdXhCgn33//PaNHj+a+++6jd+/ePP7440yePJnOnTtz2mmn8eCDDwbl2hKJpOUTrMzWydn5JwOPoOYuiG5M5aVA3eXlKCauW1w/Xa5OALNtWTZTpOUJBjL8q1cak6Hcs67zs6uu2WCNGj8go0UfstVKr5mtFUUZK4SoY7yvKMoYIcSKUAnji4YyWz/wwAP06tXLL7v/v/3tbwwZMoQ77rgjiBL6Zu/evaSkpLBx40av0ZJeeukl1qxZ4zI1kkgkEknroDGZrZOz8x8HBpTkZEzxKD8J1aQpuyQnw1sCtLDSErIHeyL71DJoqM81PhL7gHSjRb++puxtoMxo0Wd71H0P2Gi06B+s2dcD7xkt+j41+wqqX1MycGlNpKQWiS8fiW+pq5UBLEINZRVVOB19/eHZZ58NoSR1ef/997nkkkt8hlwdM2YM8+Z55hmRSCQSiaR+krPzk4DvgMeiQYmQSJozRov+kNlg/Qx41Gyw3oZqfnc5cIaX6m8DC8wG67vANuAhYIHbcZcfUEtWIsBDkVAUJQZ1yVFRVM9W9+XHIfyZXjss+JOQLtqZP38+TzzxhM/jI0eOZP369bWcjSUSiUQiAUjOzo9D/a2OBWKTs/Pbo/4W9wasgLkkJ8NSTxNhxS2PhO/8ERJJ9FInQ7nRoi8wG6yuDOVGi36T0aJfZDZYnwa+5888ErMAzAar0w/oGLDdbHCZW08zWvTvhrU3YaCWaZOiKNWontwKbh7dNVQDTwghTGGT7k+56jVtilZsNhuXXnopJSUl9UZlGjt2LP/85z8555xzwiidRCKRSCKJP6ZNNQ7YszyKZ6P+RpuAQ+4HSnIyOhFBWqLJjOxTy6A19tmJ2WC932jRP+Ol/D6jRT83kLY9TZsGoyoRS1AjQDgRwC4hRItengk2CxYs4KabbmowtKvT4VoqEhKJRCJxpyQnw4SqMHhjto9yiUQicecRoI4igWqSFTxFQghRqihKLGrijO1CiGOBNN6aOXHiBO+88w5Lly5tsO7YsWP5/PPPQy+URFIPxcXFZGZmsnbtWoYNG0ZeXh4pKSmRFksikUgkkoBxyyPhmT+ixVLjBA4QazZYz6O2y0IK3iNSNYo6ztZCiCpFUQbjI+u1xD++/vprhg4dSmpqaoN1x4wZwz/+8Y8wSCWR+CYzMxOHw0F1dTUOh4PMzEwKCqI/LL1EIpFIGkdrDP9ak0eiX8M1WxTOaD7tUX0/nAhgO2qiuoDwFf71FlTTplnAFtz8JYQQ1X41rChDARvwiRDihpqy8wEzMAhYDkwRQpT60Vaz85G48sorycjI4LbbbmuwrhCCHj16sGbNGvr27RsG6SSSusTFxVFVVeXaj42NpbIyrPEVJJJWRWPCvzYXWqIduuxTy6A19tmJ2WB922jR3xSKtn2tOrwJ3AQUoybSOIEaJeJEI9o2AyudO4qi9AA+Ax5GDSG7Cviw8SJHP7t27eL7779n0qRJftVXFIUxY8awfHmd1B0SSdgYNmwYarA2iImJYdiwYRGWSCJp2bhltjZFWhaJRNJycVcizAZrjPsWaNu+8kgMDqRRRVGuA/YDPwIn1RT/BSgQQnxcU8cE7FYURSOEiGTa8KDz7rvvkpmZSZcu3lJxeMfpcO2ZcVvarUvCRV5eHqNGjWLfvn0MGTKEvLy8SIskkbRogpXZWiIJBXaNthuqacxFwG7gH1qH/T0v9bKAu4GhwAHgPeABrcMul7SjBLPBeirqBP/JqGZO8GeE1vojAjWAV01ECFFaY3K0GTju3PfTDKkL8CjwN49D6cBqt2scAopqyr21M1VRlFWKoqzyryvRgRCC+fPnc/PNNzfqPF8rEhkZGdjtdqqqqlx26xJJKEhJSWH06NGAmv9EKqwSiUQCcd3i+ulydUKXqzNFWpYwY0a1SukNTAZetWu03sZsHYEZQA9gLHA+cH+YZJT4Ry5qzotRqE7WKaiLBgH/0HtdkVAUJRF4Bbga1ZwpXlGUy4AxQoiHGmjzMWCeEGKz00yihk7ALo+65UBnb40IIV4HXq+Rp9k4SPz++++Ul5dz7rnnNuq8MWPGsGrVKqqqqoiNjeXQoUO88sorOBx/LtZUV1ezdu3aIEsskfzJxo0bGTFiBCUlJZx55pmRFkcikTQzWmJCusq9lWW2LFurctK1a7TxwFXAcK3DXgEss2u0XwI3AtnudbUO+6tuu1vtGu27wHlhE1biD0nAg0aLPujjaV+2URbUQX4SqjYK8BNwbX2NKYoyErgAeM7L4QrA09anC0EIPRVNzJ8/n6ysLGJiGmd2duDAAY4cOULbtm3p3bs3SUlJrFq1iiFDhki7dUlYqKqqorS0lPHjx1Na2uDio0QikdRhywFRhqlcaSlKREvFGbXJbTN5VEkFqrQO+zq3stX4sCLx4BxAhvyLLv6NaqIWdHz5SJwP9BNCnHCuBgghdimK0quB9s4FkoFNNYPfTkCsoihpqMpJlrOioijxwBBa0Mt27Ngx3n///SY5TWdmZnLihOrLvnPnToYMGcKHH37o8pEoLCxk4MCB0m5dEjK2bt1Kjx490Gg0/P7775EWRyKRSCQhYntlZZnWYa9vlaUT6oSyOz6tSJzYNdqbUc1nGg5ZKQkn7YF/mw3WZahhX10EGs3J17R5OaqtmwtFUQYB2xpo73VU5WBkzWYB8oGLUbWh4YqiXKUoSnvULHt/tCRH64ULF5Kent4k23JPk6WSkhJAtVsvKCjg7bffZvDgwdJuXRIyNm7cyODBg0lOTna9fxKJRCJplTTaisSu0V4B5ACXaB323aETrWk4E9KZDVZTpGWJAIXAU8D/UP2T3beA8LUi8SbwqaIoDwIxiqKcDsxBVQx8IoQ4DBx27iuKUgEcFULsqtm/CngZeAc1j8R1gXYgGnCuGtjtdvr27UtxcXGjB/zDhg1zJQPzZsJ0/fXXM2vWLH788UfOOOOMYIovkQC43tukpCRp2iSRSCStm3VAnF2jHap12NfXlI3AhxWJXaOdALwBZGgddluYZGwUrTQhHQBGi352qNqulZBOUZQYIUS1otol3QNMRfWT2AS8BrwQicxw0Z6QLj09vZYSoNFoGp0R2J8wrxaLhYULF7Jw4cJgii+RAPDwww8TGxvL//3f/9GjRw8OHz6MR8AEiUQSRGRCuuZBa+2TXaP9ADU86G2oViZfAWdoHfYCj3p64GPgSq3D/t/QSBw4LfE5NgazwdoWGIZqceT6cTda9NZA2vVckdiqKMq/gLeFEM8DzwfSeGth7dq1VFerCb+bGlnJacJUH1OmTOHRRx/l999/Z+TIkU0RVSLxycaNG7nooouIj4+nU6dO7Ny5k969e0daLIlEIpFEhjuAt4CdwB5gutZhL7BrtINQTWXStA77JtREwwnAV3aN1nnuUq3DfkkEZJZ4wWywnoWq7LVDNVE7gOrvspkAQ8B6KhIG4AZgpaIodmAB8J4QIups3aKJgQMHumzKQxlZqX379tx3333k5OTwwQcfhOQaktZLcXExgweruSiTkpIoKSmRioREIpG0UrQO+17gCi/lm1CdsZ37MtRr9PMc8LTRon/ObLDuM1r03cwG6yO4uSM0lVrO1kKIL4QQ1wB9UU2ZJgFbFEX5ssZJuk2gF2wMiqJkKoryejiv2ViOHDmCoigMGDCA2NhYNBpNSCMrTZs2jcWLF7Nu3bqGK0uCSnFxMenp6cTFxZGenk5xcXGkRQoq7r49ycnJ0k9CIgkxsZ26tUnOzhfJ2fmmSMsSLFx5JEwJpkjLIvGNW/hXU6RlkYSFVOAFj7Ic4N5AG/aV2Xq/EOI1IcRZgBZYharNNBS1KagIIfKEEFPDec3GYjKZOO2009i8eTOVlZUUFBSENLJS586dufPOO3n66adD0n5LHywHQmZmJg6Ho0VmGT98+DDl5eX07dsXQDpcSyRhoKpi74mSnAylJCfDFGlZgoXMI9E8qAn/qmgddlOkZZGEhXL+jMK1zWywpgFdcVtZair1Zk1TFKUdMBo15XlvICo98SPFihUrWLBgAWazOazXzczMZP78+SEZ7DujT7XEwXKgBMMXJlrZuHEjSUlJrkSKTtMmiUQikUgkzZ7PgEtrPs8Dvgd+QfWbCAivioSiKGfVmBTtAB4HfgZShRDSDq6Go0ePcvPNN/PCCy/Qq1dDefqCy4033ogQIiSDfYfDgTNCVksbLAfKsGHDWmyWcc+QxdK0SSKRSFTiusX10+XqhC5XZ4q0LJKm05rzSBgt+hlGi/69ms/PAlcDt6NGZw2IWs7WiqKYgBuBbqhaSoYQ4n+BXqQl8uijjzJs2DCuvfbasF977dq1tQb7drudr776igkTJrhmlJvCjh07aoX7bGmD5UDJy8vjlFNO4cCBA6SkpLSoLOMbN26spUjIFQmJRCJRqdxbWWbLsrXK/AMtidacR8ITo0W/NFhteUZtGgc8CHwuhDgarIu0NFatWsW8efNYvXp1ROLseyav69OnDw8//DD33HMPd955J1OmTCEhIaHR7d55553cfvvtfPzxx+zduzfkjuPNjZSUFLRaLZs2bWL27NktKsu4e8Qm+NNHQgghc0lIog5/8u5IJBJJa8ZssC5FzQNSL0aL/pxAruMZtWmCEOIDqUT45tixY9x8883MnTuXPn36RESGvLw8NBqNK0rU0qVLWbVqFbm5ufz8888MHjyYO++8E4fD4Xebn332GTabjeeee47HH3+cW2+9NeSO482RjRs3ctVVV/HLL79EWpSg4mnalJiYSGxsLHv37pUO+M2clvj8WnLgA4lEIgkSb6L6QzS0BYYQIuo3Vczo4KGHHhKXXXaZqK6ujrQoPtm6dat4+OGHRe/evcVZZ50lBg4cKGJjY0VaWpooKiqqU3/Pnj2ib9++YtmyZUIIIb755htx3nnnhVvsqOfQoUOiXbt24uuvvxbjx4+PtDgBUVRUJNLS0kRMTIxo27atAERKSkqt9+Pkk08Wv/zyi0hLSxOKoghAxMTEiLS0tAhKLmkszufckp5fbGysQJ1pE4CIjY2NtEhNAjgkouA3NpgbsCrSMsg+Nbz1iYsrKxymEYXDNKZIyyKfY/PelJqbG9UoiiKiQc5ff/2VCRMmsHr1aleYzGjm2LFjpKSkUFZWBoCiKGg0GgoLC2vVmzJlCl26dOHFF18EoKioiPPPP1/ayHtQWFjIlVdeyY8//sjgwYPZv39/QD4pkSQ9Pd1lHuckJiYGjUbjyrB++eWXM2XKFK6++upa9WJjY6msrAy7zJKmERcXR1VVlWu/JTy/9PR01/8xz/e2OaEoymEhRHyk5QgmiqKsEkKMirQcwUT2qWXQGvtsNlhPA44ZLfo1Nfs9geeB4cBPwP1Gi74ikGtE9SgomhLSHT9+nJtvvplnnnmmWSgRAO3atWPHjh2ufSEEdrudDz/80DUwXLRoEUuWLGHOnDmueoMGDWLbtm0cP3487DJHMxs3bmTw4MF0796d7t27s379+kiL1GTcQ9k6cY/SVVxczLJly/jLX/5SR9mQDvjNC/fn1VKe32effeZS4gcOHCh9uaIImZBOIokqngfc7fDfRE1O9zqqMhFwUrKoViREFCWke/LJJxk4cCA33nhjpEVpFMOGDXP94MbExDBo0CDmzp1LWloaAwcO5JJLLkEIwc6dO13ntGnThv79+8sVCQ+cigTAaaed1qz9JE466aQ6Ze6DzMzMTPbt2weoK1lxcXGuFS05aGtefPbZZy6H+Zby/DZv3sxpp53GbbfdRnZ2drP15WqJma1lQjqJJKrQAksBzAZrInAJMNlo0ZuB64GAHcyiWpGIFlavXo3ZbOa1115rdhFsPB2zv//+e37++WcOHz7Mli1bAPVH2dNZcciQIc3KKTMcDqUlJSUkJycDMGrUKFatWhX0a4SL2bNn06FDB2JiYmjbtq3r/XAOMt1DDDvtIIcNGyYd8Jsh5eXlpKam0qZNG/74448W8fw++eQTrrrqKoYMGUJRUVGkxWkyLTGzdWtFl6vrpsvV/VuXqzuky9WV6nJ1f/XjHGtNfgrPCJqSENBK80jEAU7zknHAdqNFvw7AaNFvBhIDvYBUJBrgxIkTTJkyhaeeeor+/ftHWpxGk5KSQkFBAZWVla5BoKIoLr8J8J54rjn9QB87doxzzjmHwsLCkEZxaUkrEnv37mXy5MlUVVVx7NixWu8H1F3JGjp0qCscrKR5sWzZMvR6PT169GDbtm2RFidgqqqq+Pzzz1uEIiFpUZhRB2y9gcnAq7pcXbqvyrpc3WTqhuCXhJCaPBKK0aI3RVqWMFIAXFPz+TrgO+cBs8HaHygP9AJSkWiAp556it69ezNlypRIixJUPAeKnnbTKSkpUfkD7b7ykJqaitFoJCkpia1bt7rqhCojt6ci8dtvv9XxM2gurFmzhuHDh/s87rmSlZ+fT6dOnWr53EiaB0uXLuXss89m0KBBbNq0KdLiBMz//vc/evfuzUknnSQVCUlUoMvVxQNXAQ/bsmwVtizbMuBL1AS/3uonALOAv4dPSkkrZSbwmtlg3QtkAE+5HbsWCDjptFQk6sFms/HCCy/wxhtvNDuTpobwHCh62k1Hq2mTe/z49evX88EHH/DDDz+QlpbmUowURQmJQ6m7IrF//36OHj1K27Ztm2Vs/oYUCW8rWYMHD2bjxo1hlFISKEIIli1bxllnncXAgQPZvHlzpEUKiOLiYq688kpsNhvp6enExsZSVFQkV8qaG6aEOzElrMKUcAxTwoJIi9MQcd3i+tWYIDk3k0eVVKDKlmVb51a2GvC1IjEHeBXYHnxp/aNPXFw/u0Yr7BqtKVIySEKP0aJfBgwCLgRSjBa9+yxrPnBvoNeQy2o+qKys5Oabb3Y5Wbc0nANFX0TrTJ9ntKHy8nKXIuRUMjp37hx0h9L9+/dTWVlJ9+7dAVWhOXHiBIDLlKq5hJ8UQjSoSHgjOTmZkpISTj/99BBJ1jKIpqzLDoeD+Ph4Bg4c2CIUiczMTPbu3QuoffvrX/9K+/bt2blzJ717946wdJJGUAY8DlwMdIiwLA1SubeyzJZl61dPlU7UNREpBzp7VtTl6kYBZwL3AAOCJmQj2V5ZWaZ12Ovrk6SFYLToDwJ1bLE9lIomIxUJH/zzn/+kW7du3HrrrZEWJSKkpKRQXFysJhuJktWY3bt3ExcXR3V1NUKIWiZZTsVoz549DBkyhC5dugT12k5Ha+e9cDedCpUpVajYvn07MTEx9OrVq1HnyRUJ/3AqtNXV1RFXMpctW8bZZ58NqGGdm9vKmSfevnennXYaRUVFLVqRSM7Ovw7VFGYQ6iz2lJKcjKWRlSoATOWfqX8TRhHBwXQQqQA8f3S6AAfdC3S5uhjgFeAeW5atUperC5N4En8xG6zdULM9XwTsBv5htOjf81H3XlTToQ7Ap8B0o0V/rObYncAUQAe8b7Top4Rc+AgR1aZNkcojUVBQwNy5c1ukSZO/dOnShY4dO0aNTfy+ffu46KKLuPnmm9FqtT5Nsrp3785ll13GggULgnp9d7MmaNjHJJpxrkY09t12rkhI6sd91SzSSqbTPwJoESsSPXv2dH12fu+idfU0WCRn51+Iatd8M+oM9zlA1GqEbnkkRCvKJ7EOiNPl6oa6lY1AdXR1pwswCvhQl6vbDqysKd+iy9WdHXoxJX5Qx2nebLDWMVEzG6wXA9nA+UAykALMdqviXHV7K8TyRpyoXpEQQuQBeYqi3B6ua1ZWVnLLLbfw+OOPk5SUFK7LRiXOH+g+ffo0XDmEHDhwgAkTJnDuuefy7LPPNjgAnj59OjfddBP33Xdf0DJPeyoSTlMqu91OYmJis4rN3xSzJlBXJD799NMQSNSySEpKcs38R1rJXLp0KdnZ2UDLUCR0Oh2xsbFs377dZTY2f/78Fq1IoA5OHi3Jyfi5Zn9rfZUjTU0eiVZlMmPLsh3S5eo+Ax7V5epuA0YClwNneFQtB9zvzUBgBXAasCsMokrqwWywOp3mh9dke15mNlidTvPZHtWzgHlGi76g5tzHgHed9YwW/Wc15S1l1c0nUb0iEQnmzp1Lp06dmDo1KvLgRRSneVMkqaio4NJLL2XUqFF+KREA48aNo2PHjixevDhocngqEk5TKmdI1B49egTtWqFmzZo16HSNX1IP5YpEOPKAhItJkyZFRQK/LVu2cPDgQbRaLdD8FYmqqipWrlzJqlWragUBaMkrEsnZ+bGoM9g9k7PzNyRn529Jzs5/OTk7P+r9Clohd6CauOwE3gem27JsBbpc3SBdrq5Cl6sbZMuyCVuWbbtz40/lYYcty3bcV8OS4OCWR0L4yCeRClQ58yzU4MtpPr3mmHu93maDtXtQhW4GRPWKRLhxOBw8/fTTrFy5stWaNLkT6R/ow4cPk5mZiUaj4aWXXvL7mSiKwvTp07FYLFx44YVBkaWkpITzzz+/TvnAgQO58MILmT9/Pvfcc09QrhVq1qxZ0yTfn6SkJDZv3kxVVRWxsbFBlSma/AoC5ddff+Xll18mOzubP/74I+j3yl+c0Zqc35vevXu7oo21b98+IjIFwi+//MKAAQPqrJAOGTKE118PrQVsqBzonZmt3YpmeySn6w20Aa4GzgZOAF8ADwEPBiyAJGjYsmx7gSu8lG9Cdcb2dk4JEJHBhjNqEzBb67CbIiFDuKnJIxEUp3kvdZ2fOwN7mixkM0SuSNRQVVXFzTffzKOPPlpr5rk1E0lF4ujRo1x55ZUMGDCA1157rdEmSpMnT+b777+vlV8iEDxXJNyZMWMGL730ElVVVUG5Viiprq6msLCQ9HSfeZJ80r59e7p16xaSpGbR5FcQCEePHuXHH3/k2muvpV+/fvz2228Rk8XdPwJUM6v+/fu7Mto3N7777jsuuOCCOuXh+D/lHnY6mAkv3TJbKz4yXB+p+ftSSU7GtpKcjN3AXODSoAgQKUwJcZgS2gOxQCymhPaYEuTEZhipidqktBYlwk/8cpr3Udf52VvdFk1QFQlFUdopijJPUZRSRVEOKorym6Iol7gdP19RFIeiKIcVRfleUZSocUJ4/vnnad++PQaDIdKiRA2RMm06fvw4kyZNIiEhgfnz5zdpRrdz585ce+21zJs3LyBZnCY3BQUFXHfddV7vx7hx4+jWrRv5+fkBXSsclJaW0rVrVxISEpp0fqgiN7n7EUTaryAQfvrpJ9LT00lMTOS8887jhx9+iJgsS5cu5ayzzqpV1pzNm3wpEn379uXgwYMcPBi63+9IKbolORn7gC1AS0uU8RCqkpQN3FDz+aGISiSR1DjNmw3WhpzmqSkb4VFvh9Gib1WrERB806Y4YDMwHtiEOmvykaIoOlTt7TPgNiAPeAz4EBgXZBkazbp163jyySdZvnx50JxzWwKRWJGorKzkr3/9KzExMbz77rvExTX9FTUYDEycOJEHHnigye04ZyJBfU+8mdwoisKMGTN4/vnnueyyy5osbzhoqqO1E6efhPtMdzDIy8tDq9Vy/PjxiPoVBMq3337rMqc777zzeOutt7j//vvDLse+ffvYuHEjp556aq3y5prd+vDhw6xYsYJzzjmnzjFFUUhJSaGoqIiRI0eG5PoDBgygtLQUiIiiOx+4Kzk7fxGqadMMYGE4BQg6pnITYIqwFJIGsGu0dUKhah32OqFQ7RrtcOBZVKfx7lqHvVnahhst+kNmg/Uz4FGzwVqf0zzA28ACs8H6LrANVRFe4DxoNljjUMfEsUCs2WBtD1QaLfrKkHYiAgR11CyEOCSEMAkhSoQQ1UKIhcBG1JfrL0CBEOJjIcRR1H8iIxRF0QRThsZSVVXFLbfcwqxZsxgyZEgkRYk6+vbtS3l5ORUVFWG5XlVVFTfddBOHDx/mww8/pE2bNgG1N2LECAYOHBjQSoG/M5FXX301a9eu5Y8//mjytUJNcXExt912G4sWLWqyQ3OoViS6du1Ku3btaN++PcuXL49YArdAcZ81Hz9+PMuWLaOyMvy/Gz/++CNjxoyp8x1qjisSxcXFaLVaDh06xLhx47y+t6Gc9Dh27BgxMTEMHDjQZ9jpEPMYapjQdYAd+A14IpwCSFotdUKh2jVab3axJ4CPgJaQeKuO07zRoi8wG6yDzAZrhdlgHQRgtOgXAU8D3wOlNdsst3Zaz6qbECJkG+rLdxTQAC8Ar3ocXwNc5Uc7IlQ899xz4uyzzxZVVVUhu0ZzRqvVij/++CPk16mqqhJTpkwR559/vjh8+HDQ2n377bfFxRdf3KRzq6urRWJiokA1KxAxMTEiLS3NZ/0nnnhC3HLLLU0VNeSkpaX53RdfvPHGG2LKlClBl23RokXi3HPPFRqNJizvWyjYu3ev6Ny5szh69KirTKfTiZ9//jnsssycOVPMmjWrTvkrr7wipk6dGnZ5AiEtLU0oiuLzvS0qKhLdunUTiqKItLQ0UVRUFNTrP/roo+Kyyy4LaptOgEMihL/BkdiAVZGWQfYp8D4VDtPEFw7THC8cpkl1K/tX4TBNTj3nnFQ4TCMi3bfW9ByjYQuZHY+iKG1QY+rmCiEcNM4bHkVRpiqKskpRlFWhknHDhg08/vjjzJs3T5o0+SDU5k3uYT8/+ugjnn/+eTp0CF5kw2uuuYZffvmlSbPvOTk59O3bF41G49dM5NSpU/nss/9v79zjm6rSvf9bTYGWAr2BtlCg9yYtAUEYwRlEUI7nHU94UdAOOFpkHAatio4onBE1njOvoMO8WLBjwZFDvXCRo86HCuqAIKKIUrkY2hTa0EoZCnXkWm6lyXP+yOXstElz6U52kj7fz2d92qy91trrWXu12c/ez+UDNDc3d2XKAUOObNwZGRloaGiQPVzrN998g7Fjx4Zt9uyjR49i+PDhaGlpwahRoxzrMXHiROzYsSPo82nvaG0nHN9IHD582H4T4HLf6nQ6nDlzBkQkqyM0YP2OKCkpwfLly2UbM9KRJKTTKz0XuYhOih6oLdeStlyrV3ouQSQXgFlTY/QmFCrTjQnI3bMQIgrA27C+EnvUVu2LNzyIaBURjSai0YGYo8ViwW9+8xssWrQIOTk5njt0UwKtSNiTuhERrly5gsLCQlnHj4mJQVFREVauXOlTv7Vr16KsrAzbtm2D0Wh0ilvvjv79+2P69Ok+nytYXHfddY7f/bXzTk9PR319vexRbPbs2eNSkQiX/BI6nQ7/+Mc/OtzMKqFIXL58GQcOHMBNN93U4Vg4+khI/z+72rftFY3q6mpZ9gsRobi4GAsXLuz2yUl9wZaQTtj8ICKCttNtJwxFBmEoMuiVnotc2MO/Soq+XROfHv6GA5I8Enql5xJRBODVkYDVOWwHgFhJ/RwAX0k+xwG4BEDtxZgkNytWrKCbb76Z2traZB87kigpKaFHHnkkYOOrVCqHuQ0AUqlUsp/jyJEjNGDAACeTE1eYTCbKz8+nqKgoUqlUtGXLFp/PZTAYKDU1la5evervdAPGmDFjKC0tjVQqld8mIDU1NU7XS47rZrFYKCkpiZqammjp0qU0b948xzH79UAXzLGCgbt9bDd3CuZ+2LlzJ40ZM8blsdOnT1O/fv2CNhc5WLt2LcXExLjdt9I9Ii1d3S/r1q0jrVZLra2tXRXBLWDTprAo3VGm6jz1yOo89aV2dU9V56krOunDpk3dsATijcTrADQAdER0WVL/IYBhQohpQogYAM8D+J6sZk9B5ejRo9Dr9Vi9erViyaLChUC/kZDjKbkncnJyMGLECLz//vudtpMmRbNYLH5F2xk2bBgKCgrw3nvv+TvdgGA3R6qrq/Pq7Yo77r777g51QoguXbfa2lr069cPKSkpHd5IhEt+idTUVMfv0n2cmJiI7Oxs7N27N2hzcWfWBAAJCQkwm804d679g8bQpb6+Ho888ojbfVtRUeEwP5TSlf1y9uxZ/P73v0dZWVmXgz4wTJhyBEC0Ua3xJhQq042RO4/EUAC/gzVk1kkhRIut3EdEPwKYBmu0iTMAbgLwKznP7wm7mURWVhaio6P5C8ILsrKyAmpOkpGRgbS0tIBHQ7Fnuu4M6U0rEfl9EzJv3jwsW7bM/gQkJHj33Xdx7733olevXl0ax9WaREVF4YMPPui0X2cmSnv27HGY4WRkZDgd82TWEirk5uZi4MCBLvdxZ+ZNnky3fDHtsrd97rnnsHHjRpdthRAYMmRIWPlJ7N69Gzff7Cr6opXMzExUVVWhra0N+fn5Dn+3ruyXRYsWQafTdXpeObBntk5fuFkf0BMxjI9oaowXYQ3Z/x9GtSbOqNb8HNZQqG+3b2tUa4RRrYkB0NP2Ocao1nTty4YJH5R+JeJNgUymTZ6ifzAduXz5MvXq1ctrEzC7eZA35jNNTU2UkJAga5Qmd7S2ttLAgQPJYDC4bXP99dfLYhZhNpspOzubdu3a5e90ZcN+PQBQenp6lyPauDI1mjBhAr333ntu+7S2tlJaWprbtX344Ydp2bJlRER05swZiouLI4vFQkREGzduJAABi8gjBydPnqSEhAS6cOGCy+MVFRU0adIkl8c8/U/yxbTLm7Ymk4ni4uIcx0NxPaWYzWZKTEykpqYmr9pL93tWVpZf8n377beUkpJCP/30k899fQVs2hQWpbvKVJ2nTqrOU/+tOk99sTpPfaw6Tz3TVj+kOk/dUp2nHmL7nF6dp6Z2pUFpGbvDdQyFovgEvLz4JAfBsMePRNLS0qi+vt6rtr7c+JSUlND9998v0yw98/zzz9Ojjz7q8tjJkycpMTGRsrKyuuRDYGfFihU0ffp0v/vLhdzKsytFccOGDXTrrbd2aHvs2DFatGgRpaamuvSpkN70ZWRkONY7ISGBmpubiYho+fLldMcdd1BqamqX5h1Ili5d2mlI3LNnz1KfPn1c+ui0t+1v/z/Jl/9Z3rSV7gcA1LNnT1n2e6CoqqqizMxMn/vNnDmT3nzzTZ/7Xbt2jUaOHElvv/22z339gRWJ8CgsU2SU7ihzMEq3inmal5cny2vv7kZmZqbX5k2+2LSvW7cOM2bMkGWO3vDQQw9h7dq1LhPs/eEPf8Ds2bO77ENgZ9asWdi+fbsjG65SeAqd6StSMxL7Gg0fPhy7du1ymN+sXr0aU6ZMwYgRI3D27Fls3brVyeTE7lMhzRr+ww8/OCIdZWZmOvwkKisrMW3aNFy9ehVNTU1dmrvc2E2J5s+fjx07drj9G4mPj4dGo8GePXuc6k+fPg2VSgUhrElgXf1PysvLcxwH0Ome9MYMTLofAKC1tVW2yFuBwJNZkzvGjRuHr7/+2uv29mvZq1cvHD58GOPGjfP5nAzDMN2RkFYkhBA6IcQqucaTOuUpkJ00bPHW4ZqIkJiY6FTnTlmrr69HXV2dIwtwMBg8eDDGjx+P9evXO9Xv3bsXW7ZswXPPPSfbufr06YNZs2bhtddek21Mb5Ha1duVOiBwyvO0adNgsVhgNptRXV2N4uJiTJkyBY2NjVixYgUKCgocf3tRUVHo1asXNm3a5FbplDpcV1ZWYsyYMRg1ahT27dsn+9y7gj10MQA0NjZ2eiPe3k+CiDB79mz8+te/hkajAQAMHDiww/+kiooKxMfHIyoqCikpKbh48aJbheXuu+9G3759O/3/Jn2YIiVUHdm/+uor/PznP/e5n6+KhDTQwpUrVzBlyhSfz8kw4YQk/Kte6bkwYY7Sr0S8KQhgZmvGM3/84x9pwYIFnbaxWCy0cOFCysvLo9zcXFKpVKRSqWjz5s0u27/00ks0d+7cQEy3Uz7++GMaNWqUwwbfbDbT2LFjafXq1bKf6+jRo5SUlOTWdj5QtA+HGWjzFV/Mb8xmM2k0Gho0aJDbUJ3z58+nxYsX04ULF6h3797U2tpKCxYsoBdffFH2uXcFX+T++OOP6ZZbbnF8LikpodGjRzvCwq5Zs4Z0Op3LviNHjqSvv/6aiIhef/11Gjx4MNXW1jqOm0wmUqvVBIAyMzM7vcZS07SePXuGvM9YTk6OX5nOW1tbqU+fPnT27Fmv2itl9ooING1K6ydO0Av9iF7op1d6LnIVRKBJTCTKxDIrtK5KT8DLi0+Mcqxdu5buuecet8ctFgstWLCARowYQT/++KOjftmyZTR16lSXfbRaLe3cuVP2uXrCbDZTZmYmffvtt0RE9NZbb9GYMWPIbDYH5Hx33303lZaWBmRsdwT7psjXXA9Dhgzp4DMhVXJKS0tpzpw5tHPnTho7diwREb333ns0ZcqUgMrhK3b/Dm/kvnDhAsXFxdGlS5eosrKSBgwYQHV1dY7j58+fp/j4eIdviJ2rV69SbGwsXbx40VG3cuVKSk1Npezs7C4pBCaTiQYPHhyyjuzNzc0UHx/vd66fW265hT799FMi8hwEQqPRyJZ/whciUZGIxJs1likySkLcgBOv/e4zeu13n+mVnkskFcUn4NUkWZFQlG+++YZGjRrl8pjFYqGnn36abrjhBvrnP//pdOzy5cs0aNAg2rt3r1P9oUOHKC0tLWA37554+eWX6cEHH6Tz58/TwIEDHU97A8HOnTspNzc34LJKb5TcPekPxrm9uSH1pOhs2bKFJk+eTH/+858dzvEmk4kGDRoUMBn8oaKignr06OGV3CaTiWJjYykqKop69OhBy5cv79Bm5syZtGLFCqe6gwcPklqt7tDWlQO7P4pjbW0tDR061Ov2weRvf/sb3XHHHX73X7BgAen1eiLyrOyuWLGC4uLigu54zopEeJTopOgTw9YMo2FrhumVngtfJ5Y51EpI+0gwoYHdR8L2h+iAiPD0009j27Zt+Oyzz5CcnOx0PCYmBosWLerge7B+/XoUFha6tNUOBhMnTkR5eTni4+Nx8eJFp6R4cjN+/HjExcXhk08+Ccj40two1dXVMJvNAICePXsGzRfIlQN2Z3gKemB37q+srMTo0aMBWP0mLl68iFOnTgVGCD/Yu3cvHn/8ca/k1ul0uHLlCiwWC9ra2lzmNLn//vvx9tvOIdr379+PkSNHdmjb3Nzs8jy++sFkZGSgubkZFy5c8LpPsPDXP8KO1E/C7v8AuPYH+eijj7Bq1SpZAi0wkUfb6bYThiKDMBQZ9ErPhWFCDVYkGI8kJSWBiHDmzBlHHRHhqaeewueff45t27YhKSnJZd/Zs2ejpqYGX375paPfunXr8KtfBTUXoROzZs2CxWIBEeHChQsBjVYjhMATTzyBV199NSDjSyMfSTGbzSF7U+Qp6MHQoUPR2NiIb775xqFICCFcOlz7krBNbj788EPcddddXrWVRksicp3s8Pbbb8exY8ecjh04cAA33HBDh7ZSZUwI4bfiqFKpoNFoUF1d7XUfOens+vkbsclOamoqtm7dCpVK1WnQgcbGRuzdu9fra8kwDMP8L6xIMB4RQjhFbiIiPPnkk9i1axe2bt3qVokArE/GX3jhBSxatAhEhMrKSgghcOONNwZr+h2Q3qgFI1pNYWEhDAYDqqqqPLb19cZYGvnITqiHNvb0BiMmJgb9+/fHyZMnoVarHfU33ngjvvvuO6e2v/zlL2E0Gr0KYepubf1RRkwmE06dOoWxY8d6JbM3oaejo6MxY8YMvPPOO446d4qEVBnTaDQwGo1+K47Dhg3DoUOHfOojF3ZFuP31u3r1Kg4cOODIeO4PDz74ICwWi+PvIzo6GkKIDspWeXk5CgsLERsb2zVh/IAzWzNKwVGbGNlQ2rbKmwL2kVCc6dOn0/r168lisdBjjz1GY8aMoTNnznjV99q1a5Sbm0tbt26lJ598kp577rnATtYDvjoHy4Fer6c5c+bIPrf20Y/QznE5HLH7E7SXpaSkhPr27UsqlYpyc3PpoYce8sk/wN3a+rMfli5dSr/97W99kskbP5J9+/ZReno6mc1mslgslJCQQKdOnfL6PP7w8ssv0xNPPBHQc7jDnb/MV1995dYvqytjJyQkOPkrmc1mysjIoMrKyi6dy1/APhJhUVimyCjdUeagrKvSE/Dy4hOjHCaTiZKTk0kIQYmJiTR8+HCvlQg7r776KsXExBAAysrKUvRG11fnYDk4efIkJSQkdHBIb4+vEZfGjh1LaWlpIZ2d2FfcRUPKzs52WpukpCTKycnxWglwt7ae1ly6X7Kzsx3zGDJkiOzrbbFYqKCggL744guqr6+ngQMHyjq+KzZv3kyTJ08O+HlcIc20bY8eRUT0pz/9yW0Wel/Gbr83cnNz6eDBg44227dvp+HDhzvCQQcbViTCo7BMkVG6o8zBKCFt2iR3QjrGP3Q6HU6fPg0iq59Ea2srEhISfBpj1apVuHLlCgBrMjols+j66hwsB9dffz2mTp2KN954o9N2ubm5jt/tGaAB1+Y3dXV1MJlMsmXjDhXcmZ7Zk9TZOXfuHD755BOo1WoIIZCamuowWZGuV05ODrKzsx2O6IDz2mZlZTmN237NpY7sdXV1qKurAwAcP35c9n0shHA4Xe/fv9+lWZPcFBQUBMW0ydUerqiocCSxTEpKcly/rjpaA659ccaPH48vvvjC0Wb16tWYPXu2U/Zwpmuk9RMDoY8n6OP1Ss+FYZggoLQm46UWSYxyyJGXQKmET6HE/v37adCgQdTa2uq2zbPPPkt9+vShqKgoio2NdTzxdvV0dcGCBTR//vxgTT9o+GOCVFJS4mRqJH3S7aqoVCoyGo1ERPT4449TYmKiIydDcnKyI6Giu/6B3MeNjY2UlJREzzzzDD377LOyj98ei8VCffv29fi2zFvcvfHTaDQu812MGTOGli5d6vhssVhowIABdOzYMVnmI6W8vJzuvfdeIiI6e/YsxcfHyya3P4DfSIRFYZkio3AeicAUxSfg1SRZkVAUOXwKlPBLCEUmTJhA69at61DfPjux0WikpKQkx82Uu5tapc3EAoG7G9HOTNIOHDhAubm5js+dKQEqlYomTpxI69ato6tXr1JKSgodOnSIiIhyc3M9Kg/2Eqh9bDKZqHfv3gSA0tLSgnJ9x44dK1uCSOnfuhCCevbs6ZRpXXodLl26RL1796aWlhZKSUmh2tpaqq2tpbS0NFnm0p6jR49SSkoKWSwWev311ztNtBkMWJEIj8IyRUbpjjIHo4S0aRMTGngK1xmsMSIBd6FgdTqdw4SnoaEB06ZNw9SpU7Fx40YA1pCorlDaTCwQuDM968wkTavV4scff0RTUxMAICUlxeXY9ohJDz/8MMrKyvDBBx9Ao9GgoKAAAByRydyRnZ2N7OzsgO5jnU6Hy5cvAwBOnDgRlOtbUFDgVVQxb5BGEiMitLa2uo0stm/fPmg0GsTFxUGn02HTpk2ymDW5Iz09HSqVCiaTyWHWxDAMw/gPKxKMR+TwKVDCLyEU0el0aG5uxp49e5zqpXkG7H4BhYWF2LBhA4gIycnJSE1NhUqlcuoXjPC14UBUVJST/XtaWhrS0tKgUqlc3vxrtVp8+eWXmDFjBmprax0hX93lZ8jPz4fJZEJtbS1qa2sDuo9d7YVAI2cI2IyMDI9tcnJyUFFRgT179jhC6E6ZMiXgioQQArfccgv+8pe/oKmpCZMnTw7IeRgm1OHwr4xcsCLBMEFEpVLh8ccfR0lJiVO99Am6/Wlteno6vvvuO0RHR+P777/Hzp070dbWhvz8fI85CbojEyZMwM6dO1FbW4v6+nqYTCa0tbW5vPmfNm2a4ym59Km/nPkZ/MWbnBNyk5ycjDVr1siS2C83NxfXXXcdVCoVevbs6SRLfn4+dDodnn/+eWRmZmLPnj0YN24cAOC2227D/v378emnn3YpEV1nHD16FJ9//jmWLVuGK1eu4IcffgjIeRgm1DnZ1nZCU2MUmhqjXum5MGGO0rZV3hSwjwQTQZw9e5YSExOpsbHRUecqjKs0DKo0NKYS4WvDgcrKSsrPz6dnnnnGoxN6KDv/K3F9pb4hrnw/PM3JfjwqKoqio6PJYDC47ffGG29QYWEhEVnzoNTV1Tna9u3blwCQRqMJiNyh5qsF9pEIixKdFH1i2JphNGzNML3Sc+HrxDKHWhG2xQ1phBAUDvNkGG+ZN28eevfujcWLF8NkMmHcuHFobGxEr169HG2io6OdQpaqVCq0tbUpMd2woLa2Fmq1GhaLBVlZWfj73//u9i1CQUEBampqYLFYEBUVBbVaLZuPQDjiaa8VFBTAaDSCiFyul3Q9hRDQaDRu1/PkyZPQaDT47rvvcNNNN6G5uRlCCI/nCIacwSa6b/K1tOK3egB4sWHJnXrFJiIjQohKIhqt9DzkhGWKDLqjzMGATZsYRgEee+wx/PWvf8WlS5ewZs0a3HfffU5KBKCMiUs4M3XqVIe5kicndHb+d6a9b0j7vWa/wQesfhvV1dVOZlDtHaw78+tISUmBWq3GK6+8gnHjxjlyOATDNyTU/qbMLaevNSy5U0SKEsEwTPcjpBUJTkjHRCrZ2dkYN24cysvLsWbNGpfRY/hm1zfcJbJzBTv/O2Pfa0IIEBFqamqQk5ODnJwcREVFWV9ft0vaZjabUVNTA51Oh7y8PMdxb27Qf/GLX2DlypX46KOPHMpIMG7y+W8q8HBCOiZUSYgbMLB07nYqnbtdr/RcIgk2bWIYhXjnnXfwwAMPgIiQn5+PioqKbn9D2xXYXKnrSM2LpAgh0KNHD5jNZifTIMBqHlRVVeUIoZuXl+dxL2dnZztC7dqvVUVFhSMMsjdjRAJCiEtEFKf0POQkEs1HWKbIoDvKHAxC+o0Ew0Qyixcvdtyw2Z/sMv7DT5u7jtS8SAoRwWw2u40aVl9fj9GjR3v9hqehocHxu/3tEb8lYpjgweFfGblgRYJhFMIXUxzGM3wj2nWk5kVSpKZGUjOo5ORkVFRUYO3atZg5c6Zf5wkFXwWG6W5w+FdGLliRYBiF4JspJtSQvtVxl8XbrrAdPHgQPXv2RP/+/bFp0ybce++9fp2H3x4xDMOEL0H3kRBCJAF4E8C/APgngH8norUe+nRrH4lVq1Zhzpw5Sk8j4lB6XY8ePRpxNuFKr2mkEqrr+rOf/QyHDx/G+fPnw87PJxTW1BsfifSFmzt8ZzYsubPT70wl8coOXR/fQSboz4W1TNpybQeZDEUGlzJpy7VPAlgAIBbA+wAeNhQZrso7687xRiajWtNBJk2N0aVMRrWmg0yaGmNQZfKENzKXzt3eQebiskkuZS6du72DzMVlk676Ok64o8QbiVIArQCuB3AfgNeFEAX+DtbZkyx/jwVqXH+PrVrVeeAqd33DSUYlrlVn6xoMOaSmOEuWLHF7AxZO+1jpNVX6nIGSw591DcbaNDU14fz58wCsfj633nprQM8nZ99Q26ud0OE7M33hZr+/M0OEDjJBHx9xMmnLtR1k0pZr7wCwEMBtANIBZAJ4MXjT9IkOMhnVmg4yGdWacJLJEx1kLp27vYPMpXO3e5LZq3EigaAqEkKIOADTADxHRC1E9CWATQDu93fMUPtyUuBLhhWJAKxrKMkYTvs41OYSStcjEv/+m5qaHL9bLBYcP348oOeTu6/cY8o9l/SFmx3fmQ1L7mxpWHJnl78zFUcf75AJ+nMt0J8Le5m05VqHTIYiQ4uhyNCZTEUA3jQUGaoMRYYzAP4TwKygTdZLjGqNQyZNjbFFU2P0KJOmxlilqTGGrEyeKJ273SFzcdmkluKySR5lLi6bVFVcNslJZh/HCXuCatokhBgJYDcRxUrq5gOYQES6dm3nALC/d74xaJNkGIZhmCAxdMFH0o9OGa7TF24eCWB3w5I7YyV18wFMaFhyZ0iGecsfoKqrLu6TJal6EfpzescnffxIALuhPxcrqZsPYAL050JSpphBMXU5L+U4yWQoMujtH7Tl2pEAdhuKDLGSuvkAJhiKDE4yacu1BwG8ZCgybLB97g/gRwD9DUWGnwIohhPZvXrVVWRkOskkdbw2qjUjAezW1BhjJXXzAUzQ1BidZDKqNQcBvKSpMW6wfXbIpKkxBk0mT6QmDq1bVPhfTjIXl03S2z+Uzt0+EsDu4rJJsZK6+QAmFJdNcpK5dO72gwBeKi6btMH22SEzgCHejhMJRAf5fH0AnGtXdw5A3/YNiWgVAE5GxzAMw3RXvP7ODBWqfzRne2gSdjJd+ccVOWVq39b+e18AQbvprrt6NeJk8kTTmR+CJXPY7fGuEGwfiRYA/drV9QNwIcjzYBiGYZhQJxK/M7u7TO3b2n8PNfkjUSZPyCVzJO5xtwRbkTgCIFoIkSOpGwGA088yDMMwjDNHAESnL9wcSd+ZRwBEQx8fcTJpy7XeyFRlOyZtdyqYZk1ecgRAtFGt8VumUDJr8pIjAKJL5273W+biskk/+ThO2KNE+Nf1AAjAQwBuALAFwM1EFJELzDAMwzD+kr5ws8vvzIYld4bvd6Y+3qVM0J8LW5m05VqXMhmKDFXt2v0rgDUAJgFogjVs6LeGIsPCYM7XG4xqjUuZNDXGqnbtXMqkqTGGnEyeKJ273aXMxWWTqtq1cylzcdmkhb6MEwkoEf71EVhj7jYDWAfg4e6iRAghegkh3hRC/CCEuCCE2C+E+D+S47cJIWqEEJeEEDuEEEMlxyba6s4JIRrcjD9PCFEvhLgohDAKIXKDIJbiBGpdhRBDhBAt7QoJIZ4KoniKEMi9KoS4QQixy3b8uBDi+SCJpSgBXtObhRDf2sb9XgjxiyCJpThdXNenhRCHbP3qhRBPtxs73dbnkm2M24Mpm40O35lhrURY6SBTOCsRNjrIZCgyVGnLtUO05doWbbl2CAAYigyfAHgFwA4AP9jKCwrN2RMdZNLUGKuMas0Qo1rTYlRrhgCApsYYTjJ5ooPMxWWTqkrnbh9SOnd7S+nc7UMAoLhskieZXY4TPDGCCBFxCVIBEAdAD2vM4SgA/warzVw6rJ7+5wDcAyAGwJ8A7JH0/RmsocPmAGhwMfZDAL4HkA9AAMgCkKS0zOG+ru3OkwHADCBdaZnDeU0BVAP4fwBUtn3aBGCK0jKH65oCSII14dE9tjX9NYAzABKVljkM1vUZAKNgDTySB+vNwK8kx78G8P9hvSGYBuAsgAFKy8yFCxcuoVKCbtrEOCOE+B7WJCbJAGYR0c22+jhYbw5GElGNpP3tAP5KROmSuihYvwBnEdFnQZx+yCLHuroY8wUAtxLRxEDOPVSRa02FEJcAjCaiatvnjQD2EdHioAgSQsj09/9vAF4mogJJ3RFb3ZtBESTE8HVdJf2Ww2ry+5jtja4BQH8iumA7vgvAu0RUFiRRGIZhQholTJsYG0KI6wHkwuqAUwDgoP0YEV0EYLLVeyLNVoYJIRptr+hftCkY3Q4Z17U9DwAol2OO4YbMa/oqgAeEED2EEHkAxgHYJuuEwwAZ11TYSvu6YfLMNLzwd12FEALAePyvQ2QBgKN2JcLGQVd9GYZhuivd8kYzFBBC9ADwLoBy25OxrsQdTrP9/BcAWgATAcwA8Bt5Zhs+yLyu0nHHw5rq/r/lmGc4EYA1/QjAdACXAdQAeJOI9so03bBA5jXdDWCgEGKGTTkrgtVkrLeccw4Huriueli/E//L9rlbxYJnGIbxB1YkFMD2puBtAK0AHrVVdyXu8GXbz1eI6CwRNQBYCeCXXZ9t+BCAdZVSBOB9Imrp0iTDDLnXVAiRBOATAP8Bq836YAB3CCEekWvOoY7ca0pEPwH4vwB+D+AUgH+F9Q3PcZmmHBZ0ZV2FEI/C+sbxTiK66ktfhmGY7gwrEkHG9vr8TVifbk8jomu2Q04xiW22vFnwLu7wYVi/PLutw0uA1tXeJxZWZ81uZdYUoDXNBGAmoreIqI2IjgNYj26i9AZqnxLRTiIaQ0RJsDpl5wH4Vs65hzJdWVchxGwACwHcZtuPkPTNFEJI30BEbCx4hmEYf2BFIvi8DkADQEdElyX1H8Lq4zBNCBED4HkA39sdAoUQUbb6HtaPIkYI0RMAiOgSgA0AnhFC9BVCpAH4LawmJN0F2ddVwl2wRmvZEWghQoxArOkRW91MW7sUAIWQ2LFHOAHZp0KIkTazpn4AlgI4TkSfBkuoEMDfdb0PwEsAJhPRUemARHQEwAEAL9jW+y4Aw2GNF88wDMMAHP41mAXAUFjfGlyB9bW5vdxnO347rDbjlwF8DkmYUQC32vpKy+eS4/1gfbJ7AUAjrF+YQmmZw31dbW0+BfCfSssZKWsKawKfvbDam58E8AaA3krLHOZrus62nudgfahwndLyhsm61gO41q5fmeR4uq3PZVjf/N6utLxcuHDhEkqFw78yDMMwDMMwDOMzbNrEMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzPsCLBMAzDMAzDMIzP/A/1SSKfG76cLQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ysteps_dict = {\n", " 'tec_top_tpp' : [1, 0.2, 0.1, 0.03, 0.002],\n", " 'tec_sub_tpp' : [6, 1.0, 0.4, 0.10, 0.01],\n", "}\n", "rg_bw_dict = {\n", " 'S' : 75e6, #NISAR-S\n", " 'C' : 64.35e6, #Sentinel-1\n", " 'X' : 109.89e6, #TSX\n", " 'Ka': 200e6, #SWOT\n", "}\n", "unit_dict = {'L':'m', 'S':'cm', 'C':'cm', 'X':'cm', 'Ka':'mm'}\n", "ylabels = [f'{x}-band' for x in iono.SAR_BAND.keys()]\n", "ylabels[-1] += '\\nSlant range delay [m]'\n", "colors = ['C0', 'C1', 'C2', 'C3', 'C4']\n", "fig, axs = plt.subplots(nrows=2, ncols=1, figsize=[12, 3.5], sharex=True)\n", "for i, (ax, tec_type, ymax, ystep) in enumerate(zip(axs, ysteps_dict.keys(), [13, 58], [1, 5])):\n", " ts_tec = tec_dict[tec_type]\n", " # omit nan values\n", " flag = ~np.isnan(ts_tec)\n", " ts_tec = np.array(ts_tec)[flag]\n", " x = np.array(dates)[flag]\n", " if tec_type == 'tec_top_tpp':\n", " print('number of {} observations: {}'.format(tec_type, np.sum(flag)))\n", " print('min / max: {:.0f} / {:.0f} TECU'.format(np.nanmin(ts_tec), np.nanmax(ts_tec)))\n", " print('mean / std: {:.1f} / {:.1f} TECU'.format(np.nanmean(ts_tec), np.nanstd(ts_tec))) \n", "\n", " # plot\n", " ax.plot(x, ts_tec, '-o', label=tec_type, color='k', ms=4, lw=1)\n", " ax.tick_params(which='both', axis='x', direction='in', top=True, bottom=True)\n", " ax.set_ylim(0, ymax)\n", " ax.yaxis.set_minor_locator(ticker.MultipleLocator(ystep))\n", "\n", " # predicted range delay\n", " ysteps = ysteps_dict[tec_type]\n", " yaxis_locs = [1.0, 1.12, 1.24, 1.36, 1.50]\n", " for j, (bname, c, yaxis_loc) in enumerate(zip(iono.SAR_BAND.keys(), colors, yaxis_locs)):\n", " # calculate\n", " ts_delay = iono.vtec2range_delay(ts_tec, inc_angle=42, freq=iono.SAR_BAND[bname])\n", " if tec_type == 'tec_top_tpp':\n", " # mis-registration in meter\n", " dmean, dstd = np.nanmean(ts_delay), np.nanstd(ts_delay)\n", " dunit = unit_dict[bname]\n", " scale = 1 if dunit == 'm' else (100 if dunit == 'cm' else 1000)\n", " msg = f'{bname:2s} mean / std: {dmean*scale:5.1f} {dunit:2s} / {dstd*scale:6.3f} {dunit:2s}'\n", " # mis-registration in pixel\n", " if bname == 'L':\n", " for rg_bw in [24e6, 44e6, 80e6]:\n", " rg_delay_pix = dstd / (speed_of_light / rg_bw / 2)\n", " print(f'{msg} -> {rg_delay_pix:6.2f} pixel for {rg_bw/1e6:.0f} MHz')\n", " else:\n", " rg_delay_pix = dstd / (speed_of_light / rg_bw_dict[bname] / 2)\n", " print(f'{msg} -> {rg_delay_pix:6.3f} pixel for {rg_bw_dict[bname]/1e6:.2f} MHz')\n", "\n", " # add y-axis\n", " ax2 = ax.twinx()\n", " ax2.spines['right'].set_position(('axes', yaxis_loc))\n", " ax2.tick_params(which='both', axis='y', colors=c, width=1)\n", " ax2.yaxis.set_major_locator(ticker.MultipleLocator(ysteps[j]))\n", " ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())\n", "\n", " # ylim sync btw. ax and ax2\n", " ratio = ((np.max(ts_delay) - np.min(ts_delay)) / (np.max(ts_tec) - np.min(ts_tec)))\n", " ax2.set_ylim(np.array(ax.get_ylim()) * ratio)\n", " #ax2.plot(x, ts_delay, 'o', mfc='none') # test\n", "\n", "# axis format\n", "pp.auto_adjust_xaxis_date(ax, dates, every_year=1, buffer_year=None)\n", "fig.tight_layout()\n", "fig.subplots_adjust(hspace=0.1)\n", "\n", "# label\n", "top_tec_mean = np.nanmean(tec_dict['tec_top_tpp'])\n", "top_tec_std = np.nanstd(tec_dict['tec_top_tpp'])\n", "for ax, num, label, yloc in zip(axs, ['(a)', '(b)'], ['Topside TEC', 'Sub-orbital TEC'], [0.09, 0.8]):\n", " ax.annotate(num, xy=(0.015, yloc), xycoords='axes fraction', ha='left')\n", " ax.annotate(label, xy=(0.180, yloc), xycoords='axes fraction', ha='left', color='k')\n", "axs[0].annotate(r'Mean $\\pm$ STD = {:.1f} $\\pm$ {:.1f}'.format(top_tec_mean, top_tec_std), xy=(0.45, 0.09), xycoords='axes fraction', ha='left')\n", "plt.annotate('Vertical TEC [TECU]', xy=(0.002, 0.55), xycoords='figure fraction', va='center', rotation='vertical')\n", "for ylabel, color, xloc in zip(ylabels, colors, [0.593, 0.662, 0.723, 0.793, 0.88]):\n", " plt.annotate(ylabel, xy=(xloc, 0.54), color=color, xycoords='figure fraction', ha='center', va='center', rotation='vertical')\n", "\n", "# output\n", "out_fig = os.path.join(proj_dir, 'offset_comp/topTEC_pred.pdf')\n", "print('save figure to file', out_fig)\n", "#plt.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure - Percentage of TOP TEC" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "percentage min / median / max: 8 / 31 / 91\n", "34.3021240234375 - 0.3429257273674011·x¹ + 0.002435453934594989·x² +\n", "4.556585190584883e-05·x³ - 4.718176285223308e-07·x⁴ +\n", "1.4302660167331283e-09·x⁵ - 1.3914712481202796e-12·x⁶\n", "write data to file: /Users/yunjunz/Papers/2021_Geolocation/figs_src/data/top_tec_perc_s1a.txt\n", "Manually copy over the text file to ~/tools/dev/tools/data to be used by iono_tec.py.\n" ] } ], "source": [ "# calc top/total TEC percentage time-series\n", "perc = tec_dict['tec_top_tpp'] / (tec_dict['tec_top_tpp'] + tec_dict['tec_sub_tpp']) * 100\n", "print('percentage min / median / max: {:.0f} / {:.0f} / {:.0f}'.format(np.nanmin(perc), np.nanmedian(perc), np.nanmax(perc)))\n", "\n", "# fit\n", "flag = ~np.isnan(perc)\n", "doys = np.array([x.timetuple().tm_yday for x in dates])\n", "perc_ts = np.array(sorted(zip(doys[flag], perc[flag])), dtype=np.float32)\n", "ffit = poly.Polynomial(poly.polyfit(perc_ts[:,0], perc_ts[:,1], deg=6, full=True)[0])\n", "print(ffit)\n", "\n", "overwrite_txt = True\n", "txt_file = os.path.expanduser('~/Papers/2021_Geolocation/figs_src/data/top_tec_perc_s1.txt')\n", "if overwrite_txt and not os.path.isfile(txt_file):\n", " # calc percentage pred for each day-of-year\n", " x = np.linspace(1, 366, num=366)\n", " top_perc_pred = ffit(x) / 100.\n", " # write text file\n", " header = 'poly fit as a function of day-of-year:'\n", " header += '\\n{}'.format(str(ffit))\n", " data = np.hstack((x.reshape(-1,1), top_perc_pred.reshape(-1,1), 1 - top_perc_pred.reshape(-1,1)))\n", " np.savetxt(txt_file, data, fmt='%s', delimiter='\\t', header=header)\n", " print('write data to file: {}'.format(txt_file))\n", " print('Manually copy over the text file to ~/tools/dev/tools/data to be used by iono_tec.py.')" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "save figure to file /Users/yunjunz/data/geolocation/ChileSenAT149/offset_comp/topTEC_perc.pdf\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAADMCAYAAAA8nNe2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABO/0lEQVR4nO2deZyN9f7A359ZGMtgZN+XhKQkIjdtpCypCGlBFPd220t1cRGpK5ekIvplSUjCzdIiJZFbltwSoqxjzTI0M2aY5fP743nOOOfM2WbmzJkzM9/363VeM+f7fb7Lc85zns/z/Xw/i6gqBoPBYDAUNBEFPQGDwWAwGMAIJIPBYDCECUYgGQwGgyEsMALJYDAYDGGBEUgGg8FgCAuMQDIYDAZDWGAEksFgMBjCgpAJJBF5TEQ2i8h5EZntVtdBRH4VkXMiskZE6jrViYiMF5FT9us1EZFQzdtgMBgMoSGUK6QjwMvATOdCEakELAH+CVQENgMLnQ4ZDNwFXAVcCXQDhuT/dA0Gg8EQSkImkFR1iar+BzjlVtUD2K6qi1Q1FRgNXCUiTez6/sBEVT2kqoeBicCA0MzaYDAYDKEiHPaQmgE/Od6oajKwxy7PVm//3wyDwWAwFCmiCnoCQFnghFvZWSDWqf6sW11ZERF1C8QnIoOxVHyUKVPmmiZNmmAwGAyFmS1btpxU1cre6m+77Vo9deqsx7otW3Z/oaq359vkgkw4CKQkoJxbWTkg0Ut9OSDJXRgBqOoMYAZAq1atdPPmzcGfrcFgMIQQETngq/7UqbP8sHG6x7qoyJsr5cuk8olwUNltxzJYAEBEygAN7fJs9fb/2zEYDAZDkSKUZt9RIhIDRAKRIhIjIlHAUuAKEelp148EflbVX+2m7wPPiEhNEakBPAvMDtW8DQaDwRAaQrlCGgGkAC8CD9j/j1DVE0BPYByQALQB7nVqNx1YDmwDfgFW2mUGg8FgKEKEbA9JVUdjmXR7qlsNeLRAsPeKnrdfBoPBYCiihMMeksFgMBgMRiAZDAaDITwwAslgMBgMYYERSAaDwWAIC4xAMhgMhmKKiJQUkfdE5ICIJIrIVhHp7FTvKxPDUBH5xW63T0SGuvVdz25zzu6jo7/5GIFkMLihezahi0ahs5+w/u7ZVNBTMgADBgxgxIgRAR1br149Vq9e7bV+1apV3HXXXXme0/79+xER0tPT/R77zDPP8M477+R5zCATBcQDNwLlsbIufGQLE3+ZGAToB8QBtwOPiYizy84CYCtwCTAc+FhEvIZAAiOQDAYXdM8m2PAhJCdYBckJsOFDI5SKGMOGDePFF18M6ZhDhw5l3LhxXLhwIaTj+kJVk1V1tKruV9VMVV0B7AOuwU8mBlV9TVV/VNV0Vd0FfAL8BUBELgNaAqNUNUVVF2P5kvb0NR8jkAwGZ35cARlprmUZaVa5oUiwadMmzp49S9u2bUM6bvXq1WnSpAnLli0L6bg5QUSqApdhhWfzl4nBuZ0A7bkY1q0ZsFdVE50O85upwQgkg8EZx8oo0HKDC/Xq1WPChAlceeWVlClThkGDBnH8+HE6d+5MbGwsHTt2JCHh4me5bNkymjVrRoUKFbjpppvYuXNnVt3WrVtp2bIlsbGx9OnTh9TUVJexVqxYQYsWLahQoQLt2rXj559/DmiOn332GTfeeKNLmYgwZcoUGjRoQKVKlRg6dCiZmZkAZGZm8vLLL1O3bl2qVKlCv379OHs2e3TtRYsWcc0117iUTZw40UU1eNNNN7Fy5cqA5hkkKtmZuh2vwd4OFJFoYB4wxw7d5p5pAVwzMTgzGkuezLLf56RtFkYgGQzOlInLWXkBIyIhewXK4sWL+fLLL9m9ezfLly+nc+fOvPLKK5w8eZLMzEymTJkCwO7du+nbty+TJ0/mxIkTdOnShTvuuIMLFy5w4cIF7rrrLh588EFOnz5Nr169WLx4cdYYP/74IwMHDmT69OmcOnWKIUOG0L17d86fP+93ftu2baNx48bZypcuXcrmzZv58ccf+eSTT5g500puPXv2bGbPns2aNWvYu3cvSUlJPPbYY9nad+/enX379rkI1Q8++IAHH3ww633Tpk356aefsrXNR06qaiun1wxPB4lIBDAXuAA4Ts5fJgZH28ew9pK6qur5nLR1xwgkg8GZlt0gMtq1LDLaKjcExOOPP07VqlWpWbMm7du3p02bNlx99dWULFmSu+++m61btwKwcOFCunbtyq233kp0dDTPPfccKSkpbNiwge+//560tDSeeuopoqOjueeee2jdunXWGO+++y5DhgyhTZs2REZG0r9/f0qWLMn333/vd35nzpwhNjb7g/oLL7xAxYoVqVOnDk899RQLFiwAYN68eTzzzDM0aNCAsmXL8uqrr/Lhhx9mM2QoWbIkffr04YMPPgBg+/bt7N+/n27dLl47sbGxnDlzJsefaX5iq9veA6oCPVXVobP2l4kBERmIFZ+0g6oecup2O9BARJw/aL+ZGrzGshORWwI6G8hQ1bUBHmswhDXSsDUK1p5RcoK1MmrZDWnY2l9Tg03VqlWz/i9VqlS290lJSQAcOXKEunWzrIiJiIigdu3aHD58mMjISGrWrOmyMnM+9sCBA8yZM4c333wzq+zChQscOXLE7/zi4uJITMz+oF67dm2XsRx9uc+zbt26pKenc/z48Wx99O/fn759+/Lyyy8zd+5cevfuTcmSJbPqExMTqVChgt85hphpQFOgo6qmOJUvBSaISE+soNYumRhE5H7gFeBmVd3r3KGq7haR/wGjRGQE0Bm4Ej9GDb6Cq34JHMAy7fNFZSx9ocFQJJCGraGQCCAPeSoLDTVq1GDbtm1Z71WV+Pj4LEF0+PBhVDVLKB08eJCGDRsClvAYPnw4w4cPz/G4V155Jbt3785WHh8fT7NmzbLGqlGjRtY8Dxy4mCPv4MGDREVFUbVqVQ4dOuTSR9u2bSlRogTr1q1j/vz5zJ8/36V+586dXHXVVYQLtl/REOA8cMzpAWCIqs6zhdFbwAfAD7hmYngZy6R7k1O7D1T1r/b/92KlCkoADgL32NkdvOJLZXdOVRuoan1fLyDNRx8Gg8Hgkd69e7Ny5Uq++uor0tLSmDhxIiVLlqRdu3Zcd911REVFMWXKFNLT01myZAkbN27MavvII4/wzjvv8MMPP6CqJCcns3LlSo8rH3e6dOnC2rXZlToTJkwgISGB+Ph43njjDfr06QNA3759ef3119m3bx9JSUkMGzaMPn36EBXl+Xm+X79+PPbYY0RFRXH99de71K1du5bOnTt7bFcQqOoBVRVVjVHVsk6veXb9alVtoqqlVPUmVd3v1La+qka7tfurU/1+u00pVW1sZ3XwiS+BdFeA5+RzCWYwGAyeaNy4MR988AGPP/44lSpVYvny5SxfvpwSJUpQokQJlixZwuzZs4mLi2PhwoX06NEjq22rVq149913eeyxx4iLi+PSSy9l9uzZAY3bsmVLypcvzw8//OBSfuedd3LNNdfQokULunbtyqBBgwAYOHAgDz74IDfccAP169cnJibGRVXozoMPPsgvv/ziYswAcPToUXbs2BEUh9yiihTmJb8vWrVqpZs3by7oaRgMhjBk1apVTJ06lf/85z+AZa3422+/cemll+a575SUFKpUqcKPP/5Io0aNssqfffZZGjZsyKOPPpqj/kRki6q28lbfqlVj/WGj55ylUZE3+2wbbuQoQZ+dQnwi0BzYC/xDVX1aTRgMBkO40alTJzp16pQvfU+bNo3WrVu7CCOwfJIMvslpxti3gY+BUcAtWHGNrgj2pAwGg6EwUq9ePVQ1a+VlyBk+BZKILAGeVNV4uygO+I+qJovIGWBsPs/PYDAY8p1gbV3s378/KP0UV/w5xk7Eivz6oohEAVOBnSKyHsvB6V/5PUGDwWAwFA98CiRV/Q64HkgHNgAnsLxtnwGaq6pRihoMBoMhKPjdQ1LVDODfIvIh8AaQCjyrqsfye3IGg8FgKD7420O6HBgP1MdS0T2FFf7hcxGZBbypqpn5PUmDwWAweOHUGSLeX1LQswgK/vaQPgQ+x3J+XQe8o6orgbZYIYP+m7/TMxgMBkNxwZ9Aqg7MtrMBzgWqAahqqqqOAB701Tgn2ClzPxWRBBE5JiJv2YYUPvO6GwwGg6Fo4E8gvYYVOG8e1grpNedKVc0eoTD3TAX+wBKCLbByvD8aQF53g8FgMBQBfO4hqeoEEZkL1AF+V9XT+TiX+sBbdu72YyLyOVa626y87gAiMho4KSJNHGHQDQaDwVD48ZugT1WPqerGfBZGYFnw3SsipUWkJlb+DIdQCjSv+2BHqt4TJ3xGOTcYDAZDmOFVIInIXm91bsf9FqS5rMUSMn8Ch7BUc/8hB7nZVXWGI1Vv5cqVgzQtg8FgMIQCXyq7miIyJoA+qvo/xDd2PvcvgOlAOywhNBPL5DxXudkNBoPBULjwJZDmA7V91Dv4MAjzqGiP9ZaqngfO235OLwNTgP6OAz3ldTcYDAZD4cerQFLVh0I1CVU9KSL7gL+JyL+xVkj9sfaOfOZ1NxgMBkPRwK9RQwjpAdyOFS/vd6z4eU/bOdh7AuOwcrO3wTWvu8FgMBhygYiUFJH3ROSAiCSKyFYR6exU79UHVERutsvOish+D323E5GNdr8/i8j17se4EzYCSVX/Z+dfj1PVSqraS1X/sOu85nU3GAwGQ66JAuKx/D7LY/l7fmQHKvDnA5qMtdc/1L1TEakILAMmABWwfFiXi0icr8mEjUAyGAwGQ2hR1WRVHa2q+1U1U1VXAPuAa3DyAbX9Q0cDV4lIE7vtRlWdi5U93J12wHG7bYaqfoCl/erhaz5GIBkMBkPRpZLDN9N+DfZ1sIhUBS7DMhoL2AfUU1f2y73MZ4ZxnwJJRO4Qkele6qY76xoNBoPBEHacdPhm2q8Z3g4UkWhgHjDHNhoL2AfUAxuAGiLSV0SiRaQ/lnV0aV+N/K2QngE+8FI3Fw+6Q4PBYDAULmxf0LnABeAxuzjXPqCqegq4E0uGHMcyWFuNFfTAK/4S9F2uquu81H1HYEs3g8FgMIQpIiLAe1hBDrqoappdtZ08+ICq6lqgtd02Ckvd5zPLuL8VUikR8bY8KwuUCmRiBoPBYAhbpgFNgTtUNcWpfClwhYj0FJEY3HxARSTCLo+23kqMiJRwNBaRq211XTng38AhVf3C10T8CaStwD1e6noA//PT3mAwGAxhiu1XNAQr5c8xEUmyX/cH4AN6A5ACfIqVESIFWOVU/zxwEsusvDpwt7/5+FPZvYJlkx4HLAaO2h33xLJN7+NvAIPBYDCEJ6p6gOzWcM71q4EmXuq+8dO2b07n4y8f0hciMghL7zfBqSoeeFhVV3luaTAYDAZDzvC3QkJVPwY+FpHGwCXAKTulucGQY3TPJvhxBSQnQJk4aNkNadi6oKdlMBjCAH9+SG0c/9tCaKuzMBIRvzpBg8GB7tkEGz60hBFYfzd8aJUbDIZij78V0pe42qEfxopp5GAOliWGweCfH1dARpprWUaaVZ6LVdLOFdtY9/rXJB47S2y18rR/+haadmsepMkaDIZQ408geQr94Ou9weAdx8oo0HIf7FyxjVUjV5Ceagm4xKNnWTVyBYARSgZDIcWf2bfm8L3B4J0yXgL9eiv3wbrXv84SRg7SU9NY9/rXuZmZwWAIA0xwVUPoaNkNIqNdyyKjrfIcknjMPcSW73KDwRD++FPZlRWRg07vyzu9F/wEyjMYnJGGra0ldRCs7GKrlSfxaHbhE1utfN4nWkQwFo3Fg/Onlf0fpfk/sBDgTyDdHJJZGIoN0rB1rgwY3Gn/9C0ue0gAUTHRtH/6ljz3XRTIsmh0GJE4LBrBCCVD2OJPIPVR1UdDMhODIQc4DBeMlZ0XgmzRaDCEAn8C6QHACCRDWNK0W3MjgLwRRItGgyFU5NTs22AIGWYPJA+UifMsfHJh0WgwhAp/AqmEiIzxdYCqjgzifAwGwOyB5JmW3Vw/P8i1RaPBECoCWSHV9lFv/JAM+YPZA8kTwbRoNBhChT+BlKqqD4VkJgaDM2YPJM8Ey6LRYAgV/hxjQ7qHJCL3ishOEUkWkT0i0t4u7yAiv4rIORFZYyeVMhRlghjVwWAojOieTeiiUVxTr8o1BT2XUOFPIK0LySwAEbkVGA88BMRiZSPcKyKVgCVYCQErApuBhaGal6GACGJUB4OhsJEtMn4xwV+Cvi6hmgjwEjBGVb+33x8GEJHBwHZVXWS/Hw2cFJEmjtzuhqKHpz2Q1JjLOTH2P2ScmElk5YrE9b+T2JuvLeipGgzBx9MeajEgLGLZiUgk0AqoLCK/i8ghEXlLREoBzYCfHMeqajKwxy5372ewiGwWkc0nTpwI1fQN+YQ0bI30egkZMIWkSl05PnMDGX+cBoWMP05zaso8EtdsLOhpGgzBJ0QrIxEpKSLvicgBEUkUka0i0tmp3ut2iYjcbJedFZH9HvpuISLr7PpDIuLXIjssBBJQFYgG7gHaAy2Aq4ERQFnAPWjZWSy1nguqOkNVW6lqq8qVK+frhA2hJWHOJ+j5Cy5lev4CCXM+KaAZGQz5SOj2SqOAeOBGoDzW1shHIlIvgO2SZGAmMNRL3/OBb+22NwJ/E5HuviYTLgIpxf77pqoeVdWTwCSgC5CEa5JA7PeJIZyfoYDJOHE6R+UGQ6HG0x5qPqCqyao6WlX3q2qmqq4A9gHXAD2wt0tUNRUYDVwlIk3sthtVdS6w10v39YB5qpqhqnuA9XjQbDnjdQ9JRAKKUqmqeU5Ao6oJInIIz35N24H+TvMqAzS0yw3FhMjKFS11nYdyg6Go4bKHmjcqichmp/czVHWG13FFqgKXYd1f/4bbdomIOLZLAtm/nwz0E5F/Ag2A64DXfDXwZdTwXgADqj1QMJgFPC4inwNpwFPACqwU6RNEpCewEhgJ/GwMGooXcf3v5NSUeS5qOylZgrj+dxbgrAyG/MPhR7al95gteejmpKq2Cmg8kWhgHjBHVX8VkbKA+2a8x+0SL6wA3geeAyKxjNY2+WrgVSCpav0ABw0WY4FKwG4gFfgIGKeqqbYwegv4APgBuDfEczMUMA5ruoQ5n5Bx4nSurOxMbDyDwTMiEgHMBS4Aj9nFud4uEZGKwOd2X/OBasDHInJcVad6a+cvUkPIUNU0rMji2aKLq+pqoEnIJ2UIK2JvvjbXZt4mNp7B4BkRESyNWFWgi30vhrxtlzQAMlT1ffv9IRH5EMsuwKtACsioQUTKicgkEdlimwcedLwCaW8wFDi+YuMZDMWbaUBT4A5VTXEqXwpcISI9RSQGt+0SEYmwy6OttxIjIiXstrvtsvvs46oBfXDak/JEoCukqUAtYAyW2uwBLFO/xQG2NxgKlnyIjWdUgIbCju1XNAQ4DxyzFksADFHVeX62S24A1ji9TwHWAjep6p8i0gMr+s40u245MM7XfAIVSJ2Apqp6SkQyVPUT23JjOfB6gH0YDPlCQIIhyPmBjArQECwS12zM095oXlDVA/iIWepru0RVv/HT9msgRz+GQP2QIrjonJokIhWAo8ClORnMYAg22WJ+OQTDHjdjnmDHxjMqQEMQSFyzkVNT5pkIJDaBrpB+wvK0/Qor4OrbWBYYu/NpXoZCRoGprwLMmxT0/EAmPYYhCPiKQFIc4zQGKpAe4eLS7AngVaAC0C8f5mQoZBSo+ioHgiGo+YFMinCDE7l9IPMXgSQjIyOo8wx3AlXZVbZDP6CqJ1T1YVXtgxVnzlDcKUj1VUHlTTLpMQw2AauNPeAt0oijfPz48UGbZ2EgUIH0pZfyz4M1EUMhpiDVVwUkGKRha2h370XBVyYO2t1rDBqKI3l4IIvrfydSsoRLmSMCyebNmxk1alQwZxr2+FTZ2d67Yv0rgqtFRUMgPR/nZigsFKD6Kuh7Qzkc26QIN+TlgcxbBJKIa5txf8uWpKcXr1usvz2kdC4GPHX/ZDLxY1NuKCa07Oa6hwQhVV8ZwWAoUPL4QOYpAslf//pXdu8ufjZj/gRSfaxV0VosJyjBElAKnHDz6jUUUwpylWIwFDhBfiBbvnw506dPD/j4MxeiWbK/eq7GCjf8pTA/YP9bF7JUeFWB46qamc9zMxQizCrFUFwJ5gPZkSNHeOihh7Le9+zZk8WLi09AnIDMvkUkFsv36F67TZodKO8JVXXP5mowGAzFimA8kGVmZtK/f39OnToFQM2aNZk+fXqxEkiBWtm9CZQBrgBKAc2B0sCUfJqXwWAwFCsmTZrE6tWrARAR5s6dyyWXXFLAswotgTrG3g40UNVz9vvdIvIQsCd/pmUwGAzFhy1btjBs2LCs9y+++CI333xzAc6oYAhUIKUClYEDTmWVsCLEGgxhRUZGBidOnOD48eNZr9OnT3P+/PmsV3p6OiVLlqR06dJZr8qVK1OzZk1q1KhBlSpViIyMzNG4BRkkM1AKwxyLG0lJSfTt25e0NMsoonXr1rz00ksFPKuCIVCB9H/AlyIyCUso1QWeBrzmZjcY8pvMzEx+++03Nm3axM6dO9m1axe//vorjRMzeKrBNdSIKcufqUlsOniEWGlBXHR5UtPO8tmJ1Wz9c5vPviMjI6lbty7NmjWjWbNmXHHFFVx55ZVcfvnlHgWVI0imIy6ZI0gmEDY3/MIwx6KOpxBDT77yDr/99hsAZcuWZcGCBURHR/vpqWgSqEAaBxwB7gNq2P+/BszMp3kZDNlITU3lu+++Y82aNfzwww9s2rSJs2ddbWq6V7+UMVfcSOko6wddq1QsjzdqzI+nSnAoWahYogK9qncH8CmUMjIy2Lt3L3v37mX58uVZ5eXKleO6667j+uuv5/rrr6ddu3aUKFGiUATJLAxzLMp4ivmYvm4eKTvWZR0zdepUGjZsWEAzLHgCFUjXqupM3ASQiFwLFM846YZ8R1XZsWMHK1eu5Msvv2T9+vWkpqb6bPNc4zZZwshBVAQ0i0vkUHJpAEpElOC+y+6he9cenDt3jnPnzpGUlMTx48c5fPgwR44cybJ0cufPP//kj/8eIeX302yet5avMpZzvHYCL5Ys6TExjLfgmQWBv0CegWCSEuYBDyGGoshkXM92LPj+N+677z4eeOCBAppceBCoQPoSKOeh/HPAc3RAgyEXZGZmsnnzZpYsWcKSJUuyVBneqFy5Mtdeey0tWrSgcePG1FywweNxpSNd3ebkHIwePdprvykpKezevZvt27ezfft2fvnlFzZu3Ej1c5XpVb07JSKs+GMVospT+nApTtQ4TJWYmGz9eAueWRBEVq5o5d3xUB4IJilhHvESSqjOJbHUq1ePqVOn4pSxtVhiYtkZwoLdu3czd+5c5s6dy4EDB7we17hxYzp27MgNN9zAtddeS926dV1+xAe/+tXjTfdchquHQ2y18j7nU6pUKa666iquuuqqrDJVZdpNk0g5kexybImIEvySEEvbyudcVmfnyST+sipUPX+ekiVL+hwvFMT1v9NlDwkuBvIMiABzTxm84CXEUPypRObPn0/58r6vyfxAREoCU4GOWIuL34FhqvqZXd8Bywe1DlYK8wGOgAkicjMwEmgJJKhqPad+6wA73IYrAzynqhO9zcfEsjMUGAkJCSxcuJD333+f//73vx6PKVOmDJ07d6Zz58507NiROnXq+OzT0003PVPYnhCb9T4qJpr2T9+S4/mKCCknkz3WJZy/hOE7vuDZRq2pEVOWI6lJ/HvXDyz7bAaVJ49j4MCBDB48mAYNGuR43GDhLZBnwPtHJilh3vAQYij5fBpbqEHP664rqFlFAfFYCVgPAl2Aj0SkOVYS1iXAw8ByYCywEGhrt03G2sZZAAxz7lRVD+KUnkhE6mMJO59evjmNZZc1HiaWnSGXbNmyhbfffpsFCxZ43BOKi4vjzjvv5O677+bWW2+lVKlSAfft6aZ74cqWnF15CM6dJbZaedo/fQtNuzXP1dxjq5Un8Wj24CTlqldgxjfrWLRoES/On8+GDRdVhydOnGD8+PGMHz+ezp0789xzz3HzzTcXiHrGUyDPgDFJCfOENGzN6YQEktYtpGaF0hw8lcj8X//kxelvFdicVDUZGO1UtEJE9gHXAJcA21V1EYCIjAZOikgTVf1VVTcCG0WkYwBD9QO+VdX9vg4SVfVVH1JEpBGwDfhYVR+wy7wuGX1Rt1wdfaL2w3m+ARmCQ2pqKh999BFvv/02Gzdmt4OJioqia9eu9OvXj65du4aFissTO1dsY9XIFaSnXnzKjYqJptOYbi7X2L59+5g7dy7vvvsuhw4dytbPNddcw/PPP0/Pnj1z7O9UUGTbQwIriKjJAxUQ6enp3Hzzzaxfvx6AatWq8b///Y+qVav6bCciW1S1lbf62qVq6lP1h3ise27nKJ9tPYxVFcu1pwXwN6CEqv7Nqf4XYJSqLnYq6wj8n7PKzkO/e4Cxqjrb1/iBGjWEireBrDSLIlIJ30tGr2SmZYBC4tGzrBppJcoqjkKpoB0hT548ydSpU3nrrbc4ceJEtvoWLVowcOBA7r33XipXrhyyeeUWxzW07vWvSTzmfcVVv359Ro4cybBhw/j000955513+Pzzz3E8AG7ZsoU+ffrQoEEDnnvuOf5SrQ3fv7XOZ58FjYnqnjdGjBiRJYwiIiL48MMP/QqjIFBJRDY7vZ+hqh79R0UkGpgHzFHVX0WkLOD+oz0LxGZr7AMRaY8VlPtjf8eGjUASkXuBM8AG4FK7uAc+loyB9p2emsa6178Oux94flOQjpB79+5l0qRJzJw5k5QUV81uiRIl6N27N3//+99p06ZNgVkW5VZYN+3WPOBrKSoqiu7du9O9e3d+//13Jk2axKxZs7JUlXv37uXdF6dxtsYRosUyiAjnhygT1T13rFixwiUd+dixY7nxxhtDMfTJQFZItgHbXOAC8JhdnER26+pyQGIO59AfWKyqSf4ODDS4ar4iIuWAMcCzblXNgJ8cb2x95x67PEckHit+Qcl9OULmF5s3b6Z37940atSIt99+20UY1a5dm1dffZX4+Hjmzp1L27ZtC1QYnZoyz7LI04vCOnFN/rnVXXrppUydOpUDBw4wYsQI4uKsvZfOlTtmCSMHjoeogiJxzUYODhjOvq5/4+CA4fn6uRR1Dhw4QP/+/bPe33777bz44osFOCNXbAvq97BWMT1V1aGT3Q5c5XRcGSzr6u056LsU0AuYE8jx/sy+vwNWAp+q6v8CnUQuGAu8p6rxbjeoHC0ZRWQwMBigVkwNlzp/Zr5FkWA4QgbK999/z5gxY/jss8+y1bVo0YKhQ4fSq1cvjyFRCkKtWJBRC6pUqcLYsWN54YUXmD59Ovp/nh8cE4+eRVVDLrRNiKGcsXPFNq8q3HPnznH33Xdz+rT1m6tVqxZz584lIiIs1gIOpgFNgY5uhmpLgQki0hNLDowEfnZop+xVVQkg2norMUCmqjr/sO7G0nytCWQi/lR2z2GZAb4nIlWwHGE/Bb4MZPkVCCLSAssG/moP1TlaMtq60RlgbfQ5ynNr5put/0LmpZ5XR8hA+O677xgzZgyrVq3KVtepUyeGDh1Khw4dvN5UC+rm500opx8/zcRmY4Kyj+PrRgVW3LJnn32W6StfJ+lY9kv69IUzdOzYkYkTJ9KiRYtczyOnmBBDgeNu5OKsbm3S9QoeeeQRtm7dCkB0dDQLFy6kUqVKfvt0vm7KRJbON+9qEakLDMEKlH3M6Xc6RFXn2cLoLeADLKOye52a34CroEnBssi+yamsP/C+Bmg95y9j7H+B/wL/FJFqWMLpfmCGiPyEJZw+zcl+jgduAuoBB+0PoywQKSKXA+9gnRCQsyVjRHQkCHm6sWRmHkfZC5yHNIHDm4lwmL0WAi/1PDtC+mDt2rWMGTOGr792VSuJCH369OHFF190cSr1Rm5vfqsWbWb6mE/541ACVWrFMWRkFzr1CtiYyKuwPpcRERRjGF83Kvf+bnimYzbLvQuZF6wgsHu20bJlSwYMGMDLL79MjRquK//8IJQr68LOute/dvne4KK69dNdq5g/f35W+Ztvvkm7du189ufpuqkQXb5u8GduYVsse12Cq+pqoImXum98tbWPuS0n8wl43aiqx1R1pqreA1THcoqtDiwWkedzMqgbM7CETAv79Q7W8vA2rCXjFSLS014OuiwZfVH5sio8u30kg796Mg/CaBdZGTaiFdpeRWa9WhcPcnipu6F7NqGLRqGzn7D+7tmU7ZhQEHvztVzyxP1EVqkIApFVKnLJE/fn+ilXVfn666+58cYbuemmm1yEUUREBA888AA7duxgwYIFAQkjyN3Nb9WizYx/4iOOxyegCsfjExj/xEesWrTZaxt34vrfiZQs4VKWnomLA21e9nF83ajcadqtOZ3GdCO2enkQKFO1LAnNUvg52XJ0V1VmzZpFo0aNGDNmTDYjkWDjbQUdTmGQAiW/98K87U0nHj3L889fvC0OHjyYIUM8m2Y74+m6ESSs9Hv5Sa6s7FQ1HWuptgYYapsL5go76Z8j8R8ikgSkquoJ+72vJWO+Ya2MXOOfERUFLZvBfie/EjdHwXCL9+XNETIn6kdV5auvvmL06NF89913LnWRkZE8+OCDDBs2jEaNGuV4frlRK04f8ynnU1x/tOdT0pg+5tOAV0nuDrTJaRFsT4jNCsDqILfGMF5vVF7KPVnuPbTzYYYOHcrKlSsBaz9i1KhRzJ49m8mTJ3PHHXfky/5Sfq6sQ0ko1MHeHKXPZvxJZqZ1/7juuuuYMiWw5NrF0fjKmaBIXierjGD0NdrhFGu/X62qTVS1lKre5M/TN3h4yT1YprTbezcvdV/xvsKELKHprn50W8mpKqtXr6Z9+/bceuutLsIoKiqKQYMGsWvXrqyn99zgaaXi7+b3xyHPoWq8lXsj9uZrqTN7HPVXTmNDxmXZhBHk3hjGW7uc9Ne0aVNWrFjB6tWrXVac+/bt484776Rbt27s2RP8pM3BXlkXFKGwMm3/9C1Exbg+j6dpOiuOWXuq1atXZ/HixQE7ehdH4ytnwsYPKfwoiUehlHzu4v+R0VZ8Kpf6QhDvy0+QTIdqbvTo0VmOfA7urtWYkS1uonyGEJlekbiDpyAP+VtyE1+tSq04jsdn/zyr1Mp9CJv2T9/iMQJDbo1hgtlfhw4d2LJlC8P/PomvP9xNtJbhvCay8cvvaNasGc8//zwvvvgipUtnF6i5JU8hhsKEUOyFuTtKn5NUlsavYOuf2yhRogRLliyhevXqAffn6bpRNNNHkyJFkRZIebGKExrYe0gXr4X01AzSN/xKKUC89VcY4n35EJoOQbRu3TqXqujoaP7dZxB3JEbBBSvObrBUIDm9+Q0Z2YXxT3zkorYrWSqaISO75HoOgUZgKKj+vlqylU2fnKAEZUEgRspxaamO/J6ymrFjxzJ37lwmT57MnXcWLrVafhIKK1OwvusmXa9gwIABvP/++1nlc+bMoW1bv0FlsvUFrtfNmd1n/YZKKyoEHMvO3idqC9RQ1YW2xZvDWTXsqFW5oh6YOICITKcg5TmMu5WZeZwL53YRFZNB4rHzfDd5H7s+PeExdpmDwhDvSxeN8iiUjiVdoMbj013KoqOjGThwIP/4xz+QUTM8/8CrVKTO7NAGfs+rlV1ho2fzMR5XhamZf7Il8WLezDvuuIM33niD+vXrB3V8fybs4Yj7HhJY6uD8UD+OGTOGUaNGZb0fN24cw4YN89EicEIZy66gCWiFZIciX4alw6qFFU/uRiyT7D75Nrs8EBcT4SqMIMe5WyIiqvL+nR9m27T0FYqoUMT78hIG/9l5Fy3AnAVR3bqW1em+MDIH7tSrVZEWQO542x+LiSjHJZdcwl9KxPFc4zbUSCvL4UEj2XBlHe4ZPzIoQWpzYsIeTuQ53UaAvP/++y7CaNCgQfzjH/8I6hjFhUBVdtOAkao6V0Qcv4y1wLv5M628UzLSi71GDvdyfFlLeYswEO7xvrT+Nfy4ZQs1jm+latkSHDyVyPDFG1jw/W9ERUUxcOBAhg0bliWIHIRKBWLIjrd9s6q141g0biGJ0xcRrZbFXc2YssTtOM6T7W6l12uj6NChQ57G9mXCHs4CCfJ/L2z58uUMHDgw633Hjh2ZNm1asc/8mlsCFUjNsMyuwU7Yp6rJdpyisCQt07MqMrNUeXIS7N+bWWejGhS68CppaWl8+OGHjB8/nu3bXX2Lo6KieOSRRxg2bBj16tXz2L6omAMXRnztm6UtWZoljByUjorm/nJ1uKFjR/r27cvEiRNztLnuTE5N2IsL33zzDb169SIjIwOAK664go8//thjeKz8JCEtnaVHwshoKg8Eava9HythUxYici1WBsCwJKpcJS642aYkn0/jsVmfM2vWrCwfAX94MuuMiommWVxi0E1K88uJ79y5c7z11ls0atSIfv36uQijEiVKMHjwYHbP+w8jLlRG/z7e69hFxRw4p4Qq0Kgvh+pOvVrxwpTeVK0dh4i1MnphSm869WrlVWVaI8ZK2LlgwQKaNGnCm2++SXq6e+Jn/wTDhL2osWnTJu644w7On7cscRs0aMCqVasKJA15USIgowYR6YYVDfYdrIjc44C/Ao+oavYgZmFAq1atdNPCaaT9sJTI84kuaimwEqRNmjSJG264wU9Pnjd0Y6ZNvZjc3RmB+iun5Xi++bEBGx8fz9SpU3n33Xc5deqUS13ZsmX529/+xlNPPUXsrkMh2/wtbIRqYzwvxjAHBwz3qEpNkEyu+dQ19c3VV1/NtGnTaNOmTcBzCzQpYXEgcc1Gjr+7CElIzEpTv1nOsX79+qAbkjjwZ9RQNqqqtih7n8e6785OLnpGDaq6QkQ6YyXKWwvUBXqo6pb8nFxekYatKWH71fx3wQK+jV+dVbdlyxZuvPFGOnTowJgxY3zGmPLkRX/w48D3U7ztNTmbpcekCqWrpJEcf7FdbgJaqirr169nypQpLF26NEud4KBSpUo8+eST/P3vf89Kf3Bw2NteV3tl60h4G2jkMyELNOrHN8wX3lSpjZ64n9XP9ObRRx9l9+7dAGzdupXrrruOwYMH88orr1Cxov/9v2CbsBdWEtds5MTk94lMywARapWK5dXmNxHxQOd8E0bFjbBKYR5MWrVqpZs3u8Y2S05O5l//+hcTJkzIWmo7uO222xg5cqTf4IcOAn1y9nZc1YHtiDmz0eUmlJmunNqUQXK803cS4Irr5MmTzJs3j5kzZ/Lzzz9nq69Xrx7PPPMMgwYNyuZAua/r3zyu9srUESr/pVRYm7DnN94+m9yuhL2hs5/wXK5w4OM0vxZivlJ4nD9/nn//+9+8/PLLWYkBwXo4mTBhAv379zeb8AHwe9/niDyb3cslULeH3JrOmxUSICJjAulAVUcGbzr5S5kyZRg7diyDBg1i7NixzJkzJ2sF8cUXX/DFF1/Qpk0bnnrqKXr27OlzczJQk1JvT9hRhzdAjOudLiJKiLsykuT4i3p+XxZs6enpfPHFF8yaNYtly5aRlpY9gtMtt9zCE088Qbdu3YiM9GzO4c16Lq5FdK6f2osKIbMs9OJQnXFOXRIIgmejGV/WZCVLlmT48OHcd999PP7441mx8U6ePMlDDz3Ee++9x9SpU2nePPQrnsKS0mX9+vVUP5MEHgR3IG4PhdV0PtR4XSGJyCyntzFAT2ATcACoA1yLlZa2b35PMjd4WiG589tvvzF27FjmzZuXzcihRo0a9OvXj/79+9Okicfo6wHh7Qm7bq8oj0+lqsqJP2pQof+VRFYqjV6IIDKmCRERVQHLUu6bb75h8eLFLF26lD/++CNbH6VKleKBBx7g8ccfD+gm420VV+dO9RpbXgYEFiyysBOMPaRAbrqe9pA8rZgjq1Sk9j/vyvVNXFVZtmwZTzzxBAcPHrzYb2QkTz/9NKNGjaJs2bIB9ZVXCoUT+Z5NnFv/ETGZqVxIziRxm7pqMAhshTSjwxserXVjq5dn8FdP+mxbnFZIgRo1fAgsUtXFTmU9gF6FWSA52LVrF+PHj2fevHlcuHAhW/21115L79696d69e46DiHrbcK7VvQRRMdk/+/SateEvVxMRc3Hxqip8//1p3vu/ZSxbtiwr+6Q7bdu25aGHHqJPnz45tvbxpPIpe3Kl1zBI0uulHPVfmMlLRtuc3HSdBVd6spLwc0a2m1/p2kJc62iio9Rnf/7UQ8nJyYwdO5aJEye6WN7VqlWLyZMn06NHj6Cq8Qrj9aV7NpH27QdEy8XP2v0hIdCHk4nNxnhV/T673beSyQgk94NEzgIVVTXDqSwSOK2qYWnnmBOB5OD48eNMnz6dqVOncvz4cY/HNGnShNtvv5327dvzl7/8hapVq/rsMyd7SERGk9m7u5V7yY0DB47RsEF22V+tWjUefPBBBgwYwOWXXx7gmQZGYXiCDXe8hWnyd9P19iBTo2sUJcp4EBRO/eXEKm7Hjh08+uijrF271qX89ttv56233qJhHgLnOiiMK/CMjAzOznyGOA+/xfRU4dDyCzl6ODErpMAI1A/pd+DvbmWPAsGPfV+AVK1alZEjRxIfH88nn3xCjx49su0j/frrr0yePJmePXtSrVo1GjZsyF133cXw4cOZN28e69evZ8+ePZw7Z0UF9+S7U/GxviS2uZmDtdqREhGDKpw+D//6dj+ZkRmepkbt2lWy/q9VqxZPPvkk3377LYcOHeK1114LujACO39Tu3svBoYtE2eEUU7JZfR3T2k5MlSILu3lFu7UX06SA15++eWsWbOG999/nypVLl5jn3/+Oc2aNWPMmDEuhhC5wds+akaql3Mp4EDEZ86c4Y477qB8lGdfxagYpf7KadSZPS7glbI3f8bcRpMvqgQaqeFhYKmdGfYwUBNIB3rk18QKkujoaLp370737t05deoUn3zyCcuWLWPVqlXZsnXu3buXvXv38skn2R1iS5QoQUxMDDExMURFRZGamkpKSgqp77+Kt5Vpn2G9qFu3WrbyP/44y0svvcRtt91G69atiYgITRLJcA+DFM6sWrSZq85GUKW8hxubn5uuJ6OZjdszqPxnKuX89JeTyAoO1V6NCsf57Z/3UbZMBgdP/ckw22dv1KhRzJo1i5dffpm+ffvm6rrztumf8L80z1ac7ildQsj27du56667+P333zl4S3/qVSqX/aBcCMxwNZ0XkZLAVKAjUBFr8TFMVT+z6zsAb2PZDfwADLDTniMiN2Nl8W4JJKhqPQ/9Pwk8BVQBDgJ3qupur/PJYbTv67DSlh8F/hvMxHzBJjcqO3+kpKSwZs0a1q1bx7p169i0aZPHPae8cO+9HZjx7nOULh2TVaYqRMhFwwZD+ONIs96+6Xmev1cp5bzYyaXac0aHN6hR4Tidbk/CeeGeli5E3/xAVn+Bqoccqr2G9ROz9ZmSlsHDM7/MciQHaNGiBePHj6dTp045mrc39WNeDTSCiaoyc+ZMnnzySZKTLdPuvm0bMeuR2ynhLIMLQGWdnyo7O2vDUGA2lsDoAiwAmgNJWFqwh4HlwFigvaq2tdteCzQGSmEJsXpufT8MPIGV5Xsn0ABLcHk1SyxWfkjBJjU1lZ07d7J9+3a2b9/Orl27OHLkCEePHuXo0aMezbABri7XnK5VO1E+KpaUiPOcrHeW6tfXpkmTJjRp0oRGjcqD7McKrl4SoYERRgXAvHkLGD58BAcPxlOnTm3GjXuZ++8PzIbHOV1Ex5aZDOmmVImDk39GUKXb/bm6oTkLkPY3nCO2XCaJiREk1riRWvfcne04f3tIDsH18JDTHlddiZlRNHj+/WxRPjp27Mirr75Kq1aBbU2EMg1EoDgbkGSUKs+U9Xt59q35WfWlS5dm9uzZ3NOyXoELzFDvIYnIz8BLwCVYK6J2dnkZ4CRwtar+6nR8R+D/nAWSiERgWWQPUNWvAh3blx/STlVtav8fj2cbEVS1TqCDFTViYmK4+uqrufrqq7PVqSoXLlwgNTWV1NRU0tLSKFWqFPFr9rH2lVWkp1qWTaUzY2hwJJZOV3R1W77nLhCmITjMm7eAwYP/mrUXeODAQQYP/itAQELJOV3E6h8jWP2j9b8IrHsydzc0Z7XP/83wrvYJVD3kUOHFlvO8VxIbkc6ePXuYMGECkyZNylJXr169mtWrV9OlSxdGjhzpNwxRqNJABIq7sU5kylkGX1GOzW0bseD732jSpAkfffTRRZeJYqSyFpGqwGXAduBvwE+OOjug9h6sYNu/eu4hi1r26woRmY21xfM+8JKq9wy4vvyQrlfV9fb/N3rrQFXXeqsrSEKxQsoNebG28UVxS1iX39Sr15ADBw5mK69btw779/u35fGWUK9q7TgWbwsPX3J/KyRny70jR47w0ksv8d5772ULR9WpUyeGDx9O+/bt8y3iQzAdaL1ZPu4/+Sf/3leGCRMmBDUdfF7J4wrpANaqxsEMVZ3h6Vh7W+YzYI+qDhGR94ATqvqi0zHfAe+q6mynMk8rpHbAd8CnwP1ABWAVMEFVvaYt8rpD6RBG9v9rvb28tTd4Jj9C+Tv2K47HJ6AKx+MTGP/ER6xaFH4CubBw8GB8jsrdGTKyCyVLuVpV5TXNerBxWH6t+7Y02bTLbsYFNWrUYPr06fzyyy/ce++9LoJn1apV3HjjjbRu3Zq5c+cGfV81a0XjECLJCbDhQ5do6IFy5swZNMmzhWPdSuV4++23w0oYBYGTqtrK6eVNGEUAc4ELwGN2cRLgbtVRDkgMYFyH9ddrqnpGVfcD07H2qLwSkMmMiDwjIi3s/9uKyEER2Ssi1wXS3nCR/AjlP33Mpy55cgDOp6Qxfcynue6zuFOnTu0clbvjK11EqNJZ+CM+5Tw7zqexdVs0qz4vy5+JkZZe3od5f5MmTViwYAE7duzggQcecLG627JlC/369aNu3bq89NJLHDhwIDgT9RV4NkASExMZN24c9evX58CpPz0eIwVsbl5QiPV08R5QFejpZKy2HbjK6bgyQEO73B+7sIRbjowUArXhfBrYZ///KjAJKwXF5JwMVlhZtWgzPZuPoX3c0/RsPiZPK4/88Efwlt7aW7nBP+PGvZztSbl06dKMG/dywH106tWKxdtGsi7hdRZvG5kljE5NmWdZnTnFqAu1UHKsqrcfTuCLP/7kX1+l0OOVaFaX6Yf0esmvOqxJkybMnTuXX3/9lYcffpiYmItWoceOHWP06NHUq1ePW265hTlz5pCUlBTQvNyF9a5Ji72uaALJ/nzgwAFeeOEF6taty4gRIzhz5gzDF28g+bybgPNgbr5zxTZmdHiDic3GMKPDG+xcsS2gc5g3bwH16jUkIqIE9eo1ZN68BQG1K0CmAU2BO1TV2a9lKdYeUE8RicEy8f7ZYdAgIhF2ebT1VmJEpASAqp4DFgLPi0isiNQCHgF8PkUEKpDKq+pZEYnFkphvqup7WCZ/RZpgq8OadmtOpzHdiK1eHsTaO8prXpkqtTw/2XkrLwiCKdRDwf3392XGjHeoW7cOIkLdunWYMeOdgK3svOErnUUoCdaqulGjRrz77rvEx8czbtw4atSo4VK/Zs0aBgwYQNWqVenRowdz5szhxIkTHvvyJKwjVn3FucScOdAmJSXxwQcf0LlzZxo0aMBrr71GQsJF4bXppPK/Uk1QHw7fDkvFxKNnQS8GQ/UnlBzGMAcOHERVs4xhwlUoiUhdYAjQAjgmIkn2635VPYEVw3QckAC0wTLhdnADlmruUyw/pRSsfSIHj2Gp/Y4A/wXmAzN9zifA0EHbsaRbM6Crqt4lIuWAfap6id8O/Pefa+csbwTLqCFcN6edY4OllS7D7O3pfHvk4g+3ZKnoLBVRQeMQ6u7pt8NlfqEkVOks/NE+7mk8/fRFYF3C67nu98KFCyxdupRZs2bx5ZdfeszMHBERwdVXX0379u254YYbaNu2LdWqVSP+oRGe/ZVqRlLt2ggXPylnf6C0tDR+/vlnvv32Wz777DO+/fbbbOllwMrqOmLECB588EGioi4aGHsyCNr/zne5Mj7KqzGMJ4pT6KBAIzUMBT7G0gn2tMu6AcHSM0QB8cCNXHTO+khEHM5ZS3B1zloItA3S2D4JR3WYu19HdHIygxpGUjq2LF/sPh92Vna+nsbDZY6hImTpLPxQpVacxwetnK6qPQVN7dOnD3369OHw4cPMmzePOXPmsGPHjqw2mZmZbNmyhS1btjB58mQALrnkEja2vsejlV76oQxWnS530ffqzwhWp5XkP6OnsHv3bn7++Wef4Y06dOjAk08+SZcuXdj92Q5m3vZ2ljl8+evqMeO977KuT4cGpEuc54jn/oyP8moMU9wJNGPsp0ANt+JF9ivPqGoyMNqpaIWI7AOuwXLO2q6qiwBEZDRwUkSaODtn5RfB+uEGE09qn4iMDHpUTKZUjcrEVilP7VIl8zRGMM3Iw1GoFxTesrvG9b8zpPMYMrKLx1VrTqwA3R+M3HM21axZk+eff57nn3+eXbt2sWyZFa1+w4YN2VZOp06d4nBqErVKxWYb51xGBLt2xrBrp7VPdfrCGV7Z86rPuTVv3py+ffvSt29f6tWrB3jOSXRmyU9UBg45tT2fkkZqnJVzxx1/xkd16tT2uEIK1BimuBNwYCoRaSQiw0XkbREZDtTLr9BBbs5ZzXBzzsIKZ9HMQ7vBIrJZRDZ701PnlHA03/UWGyyG9Bzpu70R7H2zwrDHFSo8BdstiIgFvqwAAyUn+2GNGzdm6NChrFu3jlOnTrFy5UpeeOEF2rVrl5V/6d+7fuBcultQ2EzYnnBRSKVlpvNrYhrtyj3JNbEDqRRtbWPXrVuX3r17M336dA4cOMDPP//MP/7xjyxhBJ6DzkYCzWJLZZvvtoSkXBkfBcMYpjgT0ApJRO4DZgArscJBNAdeFJEhqjrfZ+McYjtnzQPmqOqvIlIWcJcuZ4Fsj1K2jf0MsPaQgjEfxw80nJxOval9zmVcfL5wRHfOjbFEsFVswXgaL0r4yu4aSjr1apWn69jbg1H68dNMbDbGa4SIChUq0KVLF7p0sb7/zMxM4uPj+eWXX9ixdjONdx2nbFompzWDr8+lcj6lAqVQUkhj+5/nOX4+ChGIkXJcEdeNR8dNpdfD/q1UvanbSkVmfy5Pq1SWTmO65TgYqsPoJbchp4o7gRo17MUyJPjWqaw9MNdThNdcT8ZyzpqP5Xx1p6qmicgbQLSqPup03DZgtHPCQHfCNVJDMPAUGyw9E348VZ4ydSJcdO3lchE3LT82vE0kiaKHt6CpyekRfHHIjr0YJXyevJav9n2b55tzXg2MvEVJOZeRyRfHL5aHm8FNcTJqCFRlF4tltufM90CZYE0kn5yziiTuap8UjcoSRp1uT6Jc+UxEsMLB5MKjPT9UbJ58cgyFG085m9xVbKQr18qVQTGBzutepDcfwMt6XZ0n1aUheAQqkCYBr9hOUIhIKSzb9ElBnEuunLOKK7E3X0ud2eOov3Ia8vfBHMsoT/sbzuGWTzDHHu0QnvtmhvDD/cEoOT2CH0+V51Cy6x5KXPRFQ4Bz584xfPiIXI3n60EpECdWbz6A6U3PsSVxDt+dfZstiXM4ceE3D6MYQkGgKrt4oBqWB0UCEAcIVl6kLHIb+dt2ztqPlW8h3alqiKrOs4P3vQXU5aIf0n5ffcaVqqlXxPQuNuqhnSu20fjEu3iLbZnTlNBGxWbIKd5UYpZV3JtZ70WEzMycx7vz5s82eNBfOPbpjoBStrvjHtUdLCOEYDhBB4vipLILVCB5jfbtTDgFW3X+ksJNJ5xfeIti7By12WDILzzlYbqQeYFFR1ew9c+LGva8OIn6c2Jt3DQ1aw81KTmK2M59fe6h5ocja7DxJ5BERL3bp6UXPYFUGHF/aijoyAqhwD3PC1AgGS4NxRdHSvTEY2eJiI1i4Z5P+OHkj1n1+bH6mNhsDKgljNwz3/q7/iMiSuDpHpjbVVx+UJwEUqBm39HACOBBLAfZI1ihysepanh8a34oDk6Y0rC1FZUmDFJCG4onTbs1d1GTVZlXn2P5bAIdW608iUfP+t5D9fIbMI6s4UWgoYNeA64F/orlh1QX+CeWefbT+TO14FJcnDClYetileHSEN7cf39fjwLIeSUVqI+PN9o/fQurRq7wmvnWV1TwceNe9riHZBxZC4ZABVIv4CpVPWW/3yUiP2JFUAh7gWQsxAyG8MFTCJ9VIy1L0NwIJUebpP2ziC2bnv0AH3mOjCNreBGoUcNh4EongYSIVMIyv3aPcRcWFDcrO4OhsODNGs9fJG1/FNU9VLOHZCMifVV1AVYQ1eUi8hJWNO66WHtKH+X/FHNHw2bVWbc592H0DQZD/uAthI+/SNr+MHuohR9/KrvpwALgeSwB9DYXjRoWAEbRajAYcoTDCMFTeV4xe6iFG3+RGgRAVS+o6khVvVRVS9t//6mq2bNgGQwGgw+8hfDxF0nbUPTxt0KKFJGbsQWTJ1T16+BOyWAwFGUcRgjBsrIzFB38CaSSWAFPvQkkBRoEdUYGF4JpHmswhAvu/kqGgkFESgJTgY5AReB3YJiqfmbXd8DaqqnDxbBtB+y6m7Fii7YEEtwzP4jIfqxg2Rl20QZV7eRrPv4EUrKqGoFTQATbPNZgMBjciALigRuxDNa6AB+JSHMgCVgCPAwsB8YCC4G2dttkYCaWPcEwL/3foaqrA51MwBljDaHHU4ZLR+I9g8FgyCuqmqyqo1V1v6pmquoKYB9wDdAD2K6qi1Q1FRgNXCUiTey2G1V1LrA3WPMJyKjBUDDkl3mswWAweEJEqgKXYeWba4YV/ACwhBewxy4PlHkickJEVonIVf4O9imQVDVbmnBD6PBmBhsM81iDwVAsqCQim51eg70daMcsnQfMsfPNlQXcn37PYiVsDYT7gXpYfqtrgC9EpIKvBkZlF8YY81iDwZBHTqpqK6fXDE8HiUgEVsDsC8BjdnESVrxSZ8oBiYEMrKrfqWqKqp5T1VeBM0B7X20CjWVnKACMeazBYMhvRESwrKmrAl1U1bFxvR3o73RcGaChXZ4bFD/bQEYghTnGPNZgMOQz04CmQEdVTXEqXwpMEJGewEosE++fbXWeY1VVAoi23koMkKmqF0SkDlAb2ISliXscqAR852siRmVnMBgMxRQRqQsMAVoAx0QkyX7dr6ongJ7AOCABaAPc69T8BiAF+BTLTykFWGXXxWIJugTgMHA70Nk5QLcnzArJYDAYiim2k6uvSDyrgSZe6r7x1lZVtwNX5nQ+ZoVkMBgMhrDACCSDwWAwhAVGIBkMBoMhLDACyWAwGAxhQaEQSCJSUUSWikiyiBwQkfsKek4Gg8FgCC6FxcrubSwP4qpY5okrReQn25LDYDAYDEWAsF8h2d7BPYF/qmqSqq4HlgEPFuzMDAaDwRBMCsMK6TIgQ1V3O5X9hJW/wwU7cKAjeOB5EfklBPOrBJwsAmOEahxzLuE7TqjGKkqfWyjGqOun/gtIr+SlLlTXTVAoDAIp4IizduDAGQAisllVW+X35EIxjjmX8BynKJ1LKMcqSp9bKL8bb6jq7QU5fjAJe5UdeYw4azAYDIbCQWEQSLuBKBFp5FR2FbmPOGswGAyGMCTsBZKdpXAJMEZEyojIX4A7sXJ3+MJj3o98IBTjmHMJz3GK0rmEcqyi9LmF8rsp8oiqFvQc/CIiFYGZwK3AKeBFVZ1fsLMyGAwGQzApFAKpqCIiA4CHVfX6gp6LoWAQkdnAIVUdUdBzMRgKmrBX2QUTEflGRB7Ox/73i0iKU06RJBGpEeT+L4hIJbfy/4mIiki9YI3l1Pc3IpIgIiWD2GfIz8PuP1+//1CNlR/fiZdxrheRDSJyVkROi8h3ItI6n8YaICLbROSciBwTkWkiUiGAdioil/o5Zr+IHLd9Gh1lD4vIN3mfucsYKSKSKCJn7M/tr3YSO0OAmA8r+NyhqmWdXkeC3P8+oK/jjYg0B0rlpiMR8Wn2bwuG9liph7sHeYygnUdxIhjfSYDjlANWAG8CFYGawEvA+XwY61lgPDAUKA+0xfK9+VJESgRpmCjgySD15Y07VDUWa+7/Al7ASg1uCJBiKZBEJE5EVojICftJc4WI1HKq/0ZExtpPhIkissr9aT4HY5UXkfdE5KiIHBaRl0Uk0vUQedN+Cv1VRDr46XIu0M/pfX/gfafOuorIVhH5U0TiRWS0U109+4lykIgcBL72M1Y/4Htgtj2Oo5/ZIvKOiHxpfz5rxco86ahXEfm7iPwG/JYP57FSRB537kxEfhaRu/ycj+PYASKy3q0s60nbPr+37XESReQHEWkYSN85HSsXePtOXFZk7uOKSCcR2WVfZ1Pt78zXCu4yAFVdoKoZqpqiqqtU9We7v4EistP+/Xzh4ft/QkT2ishJEZngbaVgC76XgMdV9XNVTVPV/UBvrBv7AyISKSLDRGSP/X1sEZHaIvKt3c1PYmkj+vg4nwnAc55WXSLSTkQ22Z/NJhFpZ5ffKyKb3Y59WkSW+RgHVT2rqsuAPkB/EblCREqKyL9F5KBYq7V3RCTrAUxE7hRLQ/CnfZ5FxrcoJxRLgYR13rOwLnhH6t233I65D3gIqIKVN/65XI41B0gHLgWuBjoBzjeCNsBeLI/vUcASsYw4vPE9UE5EmtqCrQ/wgVN9MtZNqwLQFfibhxv1jUBT4DY/c+8HzLNft4lIVae6+4Gx9rz/Zx/jzF32uV2eD+cxB3jAcaCIXIX1BP+pn/PJCX2xbpRxwO9YaZzDAV/fiUfsh6mPgX8AlwC7gHZ+mu0GMkRkjoh0FpE4p/7uAoYBPYDKwDpggVv7u4FWQEssq9iBXsZpB8RgWdJmoapJwGdYhkzPYH0fXbB8EAcC51T1Bvvwq2xtxEIf57MZ+Aa337H9W1sJTMH6bCZhxcq8BCtEWWNxdTm5DwjIoEpVNwKHsFa047GEfAuse0FNYKQ9h2uxHsaGYl3vNwD7AxmjqFEsBZKqnlLVxap6TlUTsW427qGIZqnqblVNAT7CupAC4T9i6ZDPiMhnQGfgKVVNVtU/gNdxzUv/BzDZfjJciHWz6OpnDMfq4lbgV6yc9Y5z+0ZVt6lqpv00u8DDuY2255PibQARuR5LYH+kqluAPVg/RgcrVfVbVT0PDAeuE5HaTvWvquppX2Pk4Tw+ARo53SgeBBaq6gUfY+WUJaq6UVXTsW7+LYLYd64I4DvxRhdgu6ousc9nCnDMVwNV/RO4Hks1+C5wQkSW2QJwCNb3u9Pu7xWghfMqCRhvf/8Hgck4qWfdqASctPtx56hd/zAwQlV3qcVPqnoqgPN2ZyTwuIhUdirrCvymqnNVNV1VF2Bdi3eo6jmsa60vgH29NcESVIFyBEvl+QjwtP2ZJGJ9Zo77wCBgpqp+aV/vh1X111ycX6GnWAokESktItPFSmXxJ/AtUEFcVWnOP9hzWCGMAuEuVa2gqhWwVjzRwFGHkAKmY626HBxWV1PHA4A/Q4i5WDeiATipuexzayMia8RSR54F/or1o3YmPoDz6A+sUlVHLKz5OKmInPuwn2ZPu807kDFydR62EPwIS50TgXXD8OeXllNy+/3nJ/6+E2/UwPX7Uqwnd5/YAmeAqtYCrrD7mYwlFN9wuqZPA4L11O/A+fv3dU2fBCqJ573G6nZ9bSzhmydU9ResfbEXnYpr2PNz5gAXz2U+F4XpfcB/bEEVKDWx9q9KA1ucPrPPsVaXEKTzKwoUS4EEPAs0BtqoajmsJTJYP6pgEo+1CVzJIaRUtZyqNnM6pqaIOI9bB+upyiuqegDLKKALbqoOrB/QMqC2qpYH3iH7efm09bd1272BG8WyeDoGPA1cZavHwPoROY4vi/UU6Dxvv/4EeTyPOVhqww5Y6pv/+hvPiWSsG4Rj/tVy0DanBGWsAL4Tl3EA53GOAs57pOL8PhDsJ/bZWIIpHhjidE1XUNVSqrrBqYnzatnXNf1frN9ID+dCsSziOgNf2ePlag/PA6OwVisOgXOE7MFL63Bxtb4KS2C2wBJMAfs/imWRWBP4D9a2QDOnz6u8qjoecoJ5foWa4iqQYrEukDO2DnlUfgyiqkexLuiJIlJORCJEpKGIOKvQqgBPiEi0iPTC2tsJZC9kEHCLHcnCmVjgtKqm2rrp3CQzvAvIwNr/aWG/mmLtFTgMEbqIZRZcAmsv6QdVDWRV5E6uzsMWQJnARHK+OvoJaCYiLUQkBhidi3mHeqy78P2d/A/oYa/+L8X6XB2sBJqLyF32SuTvuAqsbIhIExF5VmxjH1sd2xdr7+8d4B8i0syuK29fu84MFct4qDaWdZvH/R1VPYu1V/emiNxu/w7qAYuwVnFzgf8DxopII7G40t7jATgONPB1Lm7j/W7P5Qm76FPgMhG5T0SixDKMuBxrJYWtSvwYyyiiIvClvzHs33o34EPgA1X9CUvt+bqIVLGPqSkijj3c94CHRKSDfY+oKSJNAj2nokRxFEiKpXYohaUO+B5r+Zxf9MMyitgBJGBd3NWd6n8AGtlzGQfcE4h+XFX3qOpmD1WPYoVZSsTSmX+Uizn3x9pDO6iqxxwvLMOP+7FUEPOxBPlp4Bq7PMfk8TzeB5rjagwRwJC6GxgDrMayAlzvu0muCeZY/r6T17GSWB7HWj1mGZnYKr5ewGtYkU4ux9rk92XCnYhllPKDiCRj/U5+AZ5V1aVYm/Qf2irvX7BWM858AmzBEpQr8WH+rKqvYRlJ/Bv4E+s3EQ90sNWzk7C+/1V2/XtcdBEYDcyxVWG9fZyPM2OAMvbYp4BuWFqTU8DzQDcntShY13pHYJGXvS4Hy+3rNR5rX3USlmEUWCbgvwPf25/ZaiwtjcP44SGs7/AssBb/KSeKJMUqUoOI/AiMUdX/FPRcCjMSJtEFRKQfMFgDjHQRyu8/nK81e9/tEHC/qq7Jh/4VaGSvRgyGgCk2KyRbvdAU2FrQczHkHREpjbWKCii4ZSi//3C81kTkNhGpIFZ0h2FY+3HfF/C0DAYXioVAEpHxWMv9F+yNdEMhxta9n8BST/ndZA7l9x/G19p1WJZcJ4E7sKxBfZnkGwwhp1ip7AwGg8EQvhSLFZLBYDAYwh8jkAwGg8EQFhiBZDAYDIawwAgkg8FgMIQFRiAZDAaDISwwAslgMBgMYcH/A3anTecxyrrkAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# prepare data\n", "years = sorted(list(set([d.year for d in dates])))\n", "num_year = len(years)\n", "year0 = 2016\n", "\n", "# plot\n", "fig, ax = plt.subplots(figsize=[6, 3])\n", "\n", "# plot - model fit\n", "#ax.axhline(np.nanmedian(perc), c='k', ls='--', lw=3)\n", "x = ptime.date_list2vector(ptime.get_date_range('20160101', '20161231'))[0]\n", "top_perc_pred = ffit(np.linspace(1, 366, num=366))\n", "ax.plot(x, top_perc_pred, c='k', ls='-', lw=3, label='model (poly)')\n", "\n", "cmap = plt.get_cmap('magma', lut=num_year)\n", "for i, year in enumerate(years):\n", " flag = np.array([d.year == year for d in dates], dtype=np.bool_)\n", " xs = [x.replace(year=year0) for x in np.array(dates)[flag]]\n", " ax.plot(xs, perc[flag], 'o', c=cmap(i))\n", "\n", "# axis format\n", "ax.set_ylabel('Topside / total TEC [%]')\n", "ax.set_ylim(0, 100)\n", "ax.set_xlim(dt.datetime(year0, 1, 1), dt.datetime(year0, 12, 31))\n", "# centering labels beetween ticks (https://matplotlib.org/stable/gallery/ticks_and_spines/centered_ticklabels.html)\n", "ax.xaxis.set_major_locator(mdates.MonthLocator())\n", "ax.xaxis.set_minor_locator(mdates.MonthLocator(bymonthday=16))\n", "ax.xaxis.set_major_formatter(ticker.NullFormatter())\n", "ax.xaxis.set_minor_formatter(mdates.DateFormatter('%b'))\n", "for tick in ax.xaxis.get_minor_ticks():\n", " tick.tick1line.set_markersize(0)\n", " tick.tick2line.set_markersize(0)\n", " tick.label1.set_horizontalalignment('center')\n", "ax.legend(frameon=False)\n", "\n", "# colorbar\n", "cax = make_axes_locatable(ax).append_axes('right', 0.1, pad=0.1, axes_class=plt.Axes)\n", "cbar = colorbar.ColorbarBase(cax, cmap=cmap, ticks=((np.arange(num_year) + 0.5) / num_year))\n", "cbar.ax.set_yticklabels(years)\n", "\n", "fig.tight_layout()\n", "# output\n", "out_fig = os.path.join(proj_dir, 'offset_comp/topTEC_perc.pdf')\n", "print('save figure to file', out_fig)\n", "#plt.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)\n", "plt.show()" ] }, { "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 }