From 1b951f4d89aaab86be19b8c0ca6583bb844cc2ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pawe=C5=82=20Magnuszewski?= <47724273+pawelmagnu@users.noreply.github.com> Date: Fri, 8 Nov 2024 18:35:21 +0100 Subject: [PATCH] comparison against trixi via 2d advection (#454) --- .github/workflows/tests+pypi.yml | 4 +- .../trixi_comparison/__init__.py | 7 + .../advection_comparison.ipynb | 659 ++++++++++++++++++ examples/docs/pympdata_examples_landing.md | 3 +- examples/setup.py | 1 + 5 files changed, 672 insertions(+), 2 deletions(-) create mode 100644 examples/PyMPDATA_examples/trixi_comparison/__init__.py create mode 100644 examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb diff --git a/.github/workflows/tests+pypi.yml b/.github/workflows/tests+pypi.yml index 21573c02..9adbf8e4 100644 --- a/.github/workflows/tests+pypi.yml +++ b/.github/workflows/tests+pypi.yml @@ -176,7 +176,7 @@ jobs: - uses: actions/checkout@v2 with: submodules: recursive - fetch-depth: 0 + fetch-depth: 0 - uses: actions/setup-python@v5.2.0 with: @@ -195,6 +195,8 @@ jobs: sudo make install cd ../../../ rm -rf libmpdataxx + - uses: julia-actions/setup-julia@v2 + - run: julia --version # https://github.com/numba/numba/issues/6350#issuecomment-728174860 - if: matrix.platform == 'ubuntu-latest' diff --git a/examples/PyMPDATA_examples/trixi_comparison/__init__.py b/examples/PyMPDATA_examples/trixi_comparison/__init__.py new file mode 100644 index 00000000..c58d4bbd --- /dev/null +++ b/examples/PyMPDATA_examples/trixi_comparison/__init__.py @@ -0,0 +1,7 @@ +""" +This example uses a basic 2D advection test case to compare PyMPDATA +solution against Trixi.jl (Julia DG code) + +advection_comparison.ipynb: +.. include:: ./advection_comparison.ipynb.badges.md +""" diff --git a/examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb b/examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb new file mode 100644 index 00000000..1e95b0d9 --- /dev/null +++ b/examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb @@ -0,0 +1,659 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e333839d", + "metadata": {}, + "source": [ + "[![preview notebook](https://img.shields.io/static/v1?label=render%20on&logo=github&color=87ce3e&message=GitHub)](https://github.com/open-atmos/PyMPDATA/blob/main/examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb)\n", + "[![launch on mybinder.org](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/open-atmos/PyMPDATA.git/main?urlpath=lab/tree/examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb)\n", + "[![launch on Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/open-atmos/PyMPDATA/blob/main/examples/PyMPDATA_examples/trixi_comparison/advection_comparison.ipynb)" + ] + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": [ + "# Introduction\n", + "Trixi.jl is a numerical simulation framework for conservation laws written in Julia. It is based on the Discontinuous Galerkin (DG) method and for the purpose of this comparison, we will use the StructuredMesh for data representation.\n", + "\n", + "This notebook compares the results of a simple advection equation solved in 2D by PyMPDATA and Trixi.jl.\n", + "The general flow of the notebook is as follows:\n", + "1. define the advection equation and the common settings for both PyMPDATA and Trixi.jl in the JSON file;\n", + "2. run the simulation in Trixi.jl and save the results;\n", + "3. use Trixi2Vtk to convert the results to a vtk file;\n", + "4. reshape the results from Trixi.jl to match the shape of the results from PyMPDATA;\n", + "5. run the simulation in PyMPDATA for a bigger nx and ny, to account for the polynomial degree in Trixi.jl;\n", + "6. compare the results from PyMPDATA and Trixi.jl;\n", + "7. assert that the results are close to each other, this is to ensure that the implementation of PyMPDATA is correct.\n", + "\n", + "To run the notebook, Julia and the following Julia packages are required:\n", + "- JSON\n", + "- Trixi\n", + "- OrdinaryDiffEq\n", + "- Trixi2Vtk\n", + "- Pkg" + ], + "id": "2448bffa3ee6d9ff" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:10.871429Z", + "start_time": "2024-10-30T19:29:10.857440Z" + } + }, + "cell_type": "code", + "source": [ + "import sys\n", + "if 'google.colab' in sys.modules:\n", + " !pip --quiet install open-atmos-jupyter-utils\n", + " from open_atmos_jupyter_utils import pip_install_on_colab\n", + " pip_install_on_colab('PyMPDATA-examples')" + ], + "id": "6127e1a5c94b8ece", + "outputs": [], + "execution_count": 1 + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "if 'google.colab' in sys.modules:\n", + " JULIA_URL = \"https://julialang-s3.julialang.org/bin/linux/x64/1.11/julia-1.11.1-linux-x86_64.tar.gz\"\n", + " !wget -nv $JULIA_URL -O /tmp/julia.tar.gz\n", + " !tar -x -f /tmp/julia.tar.gz -C /usr/local --strip-components 1\n", + " !rm /tmp/julia.tar.gz" + ], + "id": "6c2e6e6520f308da" + }, + { + "cell_type": "markdown", + "id": "0f162ce9-5704-4464-8b67-be8c86ecabc8", + "metadata": {}, + "source": [ + "## common settings" + ] + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:10.887433Z", + "start_time": "2024-10-30T19:29:10.873430Z" + } + }, + "cell_type": "code", + "source": [ + "SETUP = {\n", + " \"nx\": 32,\n", + " \"ny\": 32,\n", + " \"ux\": 0.25,\n", + " \"uy\": 0.25,\n", + " \"dt\": 0.025,\n", + " \"tmax\": 2.0,\n", + " \"polydeg\": 2,\n", + " \"omega\": 3.141592,\n", + " \"min_x\": -1.0,\n", + " \"min_y\": -1.0,\n", + " \"max_x\": 1.0,\n", + " \"max_y\": 1.0\n", + "}\n", + "\n", + "assert SETUP[\"nx\"] == SETUP[\"ny\"]\n", + "\n", + "import json\n", + "import subprocess\n", + "with open('setup.json', 'w', encoding='UTF-8') as f:\n", + " json.dump(SETUP, f)" + ], + "id": "dff76910f0610a2d", + "outputs": [], + "execution_count": 2 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "## Trixi.jl", + "id": "52cd27020f7efe9" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:10.903438Z", + "start_time": "2024-10-30T19:29:10.888434Z" + } + }, + "cell_type": "code", + "source": [ + "%%writefile trixi.jl\n", + "import Pkg\n", + "Pkg.add([\"JSON\", \"Trixi\", \"OrdinaryDiffEq\", \"Trixi2Vtk\"])\n", + "using JSON\n", + "using Trixi\n", + "using OrdinaryDiffEq\n", + "using Trixi2Vtk\n", + "\n", + "setup = JSON.parsefile(\"./setup.json\")\n", + "\n", + "advection_velocity = (setup[\"ux\"], setup[\"uy\"])\n", + "equations = LinearScalarAdvectionEquation2D(advection_velocity)\n", + "solver = DGSEM(polydeg = setup[\"polydeg\"])\n", + "\n", + "function initial_condition(x, t, equations::LinearScalarAdvectionEquation2D)\n", + " return SVector(sin(setup[\"omega\"]*sum(x)) + 1)\n", + "end\n", + "\n", + "cells_per_dimension = (setup[\"nx\"], setup[\"ny\"])\n", + "coordinates_min = (setup[\"min_x\"], setup[\"min_y\"])\n", + "coordinates_max = (setup[\"max_x\"], setup[\"max_y\"])\n", + "\n", + "mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)\n", + "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)\n", + "\n", + "tspan = (0.0, setup[\"tmax\"])\n", + "ode = semidiscretize(semi, tspan);\n", + "\n", + "summary_callback = SummaryCallback()\n", + "save_solution = SaveSolutionCallback(save_initial_solution = false, interval=100)\n", + "\n", + "stepsize_callback = StepsizeCallback(cfl = 1.6)\n", + "\n", + "callbacks = CallbackSet(summary_callback, save_solution, stepsize_callback)\n", + "\n", + "time_int_tol = 1e-6\n", + "sol = solve(ode, CarpenterKennedy2N54();\n", + " abstol = time_int_tol,\n", + " reltol = time_int_tol,\n", + " dt = setup[\"dt\"],\n", + " ode_default_options()..., callback = callbacks);\n", + "\n", + "summary_callback()" + ], + "id": "6586bff9a39d588f", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Overwriting trixi.jl\n" + ] + } + ], + "execution_count": 3 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:51.466613Z", + "start_time": "2024-10-30T19:29:10.905438Z" + } + }, + "cell_type": "code", + "source": "subprocess.run([\"julia\", \"trixi.jl\"], check=True)", + "id": "56fb8302adfc01e7", + "outputs": [ + { + "data": { + "text/plain": [ + "CompletedProcess(args=['julia', 'trixi.jl'], returncode=0)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "execution_count": 4 + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "## PyMPDATA", + "id": "a30cc2b4961f1be7" + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:52.484809Z", + "start_time": "2024-10-30T19:29:51.469600Z" + } + }, + "cell_type": "code", + "source": [ + "import numpy as np\n", + "import meshio\n", + "from open_atmos_jupyter_utils import show_plot\n", + "import matplotlib.pyplot as plt\n", + "from PyMPDATA import Solver, ScalarField, VectorField, Stepper, Options\n", + "from PyMPDATA.boundary_conditions import Periodic\n", + "import os" + ], + "id": "9aaadc4a5234804a", + "outputs": [], + "execution_count": 5 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:52.500805Z", + "start_time": "2024-10-30T19:29:52.485809Z" + } + }, + "cell_type": "code", + "source": [ + "dt = SETUP[\"dt\"]\n", + "tmax = SETUP[\"tmax\"]\n", + "nt = int(tmax / dt)\n", + "\n", + "nx = SETUP[\"nx\"] * SETUP[\"polydeg\"] + 1\n", + "ny = SETUP[\"ny\"] * SETUP[\"polydeg\"] + 1\n", + "ux = SETUP[\"ux\"]\n", + "uy = SETUP[\"uy\"]\n", + "omega = SETUP[\"omega\"]\n", + "\n", + "min_x, min_y = SETUP[\"min_x\"], SETUP[\"min_y\"]\n", + "max_x, max_y = SETUP[\"max_x\"], SETUP[\"max_y\"]\n", + "dx_temp = (max_x - min_x) / (nx - 1)\n", + "dy_temp = (max_y - min_y) / (ny - 1)\n", + "min_x, max_x = min_x - dx_temp/2, max_x + dx_temp/2\n", + "min_y, max_y = min_y - dy_temp/2, max_y + dy_temp/2\n", + "dx = (max_x - min_x) / nx\n", + "dy = (max_y - min_y) / ny\n", + "Cx = ux * dt / dx\n", + "Cy = uy * dt / dy" + ], + "id": "9a0f60b51e32ce3e", + "outputs": [], + "execution_count": 6 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:52.516818Z", + "start_time": "2024-10-30T19:29:52.502806Z" + } + }, + "cell_type": "code", + "source": [ + "opt = Options(n_iters=3)\n", + "boundary_conditions = (Periodic(), Periodic())" + ], + "id": "64e3274fa4ac14a6", + "outputs": [], + "execution_count": 7 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:52.532819Z", + "start_time": "2024-10-30T19:29:52.518818Z" + } + }, + "cell_type": "code", + "source": [ + "def initial_condition():\n", + " return np.array(\n", + " [\n", + " np.sin(omega*(x+y)) + 1 for x in np.linspace(min_x, max_x, nx)\n", + " for y in np.linspace(min_y, max_y, ny)\n", + " ],\n", + " dtype=float\n", + ").reshape((nx, ny))\n", + "\n", + "advectee = ScalarField(data=initial_condition(), halo=opt.n_halo, boundary_conditions=boundary_conditions)" + ], + "id": "cab790be5c425ea5", + "outputs": [], + "execution_count": 8 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:52.547815Z", + "start_time": "2024-10-30T19:29:52.534814Z" + } + }, + "cell_type": "code", + "source": [ + "field_x = np.full((nx+1, ny), Cx, dtype=opt.dtype)\n", + "field_y = np.full((nx, ny+1), Cy, dtype=opt.dtype)\n", + "\n", + "advector = VectorField(\n", + " data=(field_x, field_y),\n", + " halo=opt.n_halo,\n", + " boundary_conditions=(boundary_conditions[0], Periodic())\n", + ")" + ], + "id": "b454a74473b8f900", + "outputs": [], + "execution_count": 9 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:52.988941Z", + "start_time": "2024-10-30T19:29:52.548817Z" + } + }, + "cell_type": "code", + "source": [ + "stepper = Stepper(options=opt, n_dims=2)\n", + "solver = Solver(stepper=stepper, advector=advector, advectee=advectee)" + ], + "id": "ce055b2d3a61a491", + "outputs": [], + "execution_count": 10 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:29:53.003936Z", + "start_time": "2024-10-30T19:29:52.989935Z" + } + }, + "cell_type": "code", + "source": [ + "vmin = np.min(solver.advectee.get())\n", + "vmax = np.max(solver.advectee.get())" + ], + "id": "5290dbc73bd36d29", + "outputs": [], + "execution_count": 11 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:27.249334Z", + "start_time": "2024-10-30T19:29:53.004932Z" + } + }, + "cell_type": "code", + "source": "_ = solver.advance(n_steps=nt)", + "id": "b29adee15e8ff545", + "outputs": [], + "execution_count": 12 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:27.264339Z", + "start_time": "2024-10-30T19:30:27.250335Z" + } + }, + "cell_type": "code", + "source": "pympdata_result_state = solver.advectee.get().copy()", + "id": "f59fd725765b765e", + "outputs": [], + "execution_count": 13 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:28.048648Z", + "start_time": "2024-10-30T19:30:27.267338Z" + } + }, + "cell_type": "code", + "source": [ + "plt.imshow(pympdata_result_state, cmap='viridis', vmin=vmin, vmax=vmax)\n", + "plt.colorbar()\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')\n", + "plt.title('PyMDATA solution')\n", + "show_plot(inline_format='png')" + ], + "id": "a041cd5f2c2dbaa", + "outputs": [ + { + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAHHCAYAAADaqqCfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABydElEQVR4nO3deXxU5d3//9dM9oUkhiWBCkir4soiFEzV26WpyG35aqUu1FbE7VcLVk1729K7dWsrtt5VtEVorYq2IGqrVquiiIJfv4ILllZqbwqKgkLCYrPvM+f3R8ro5HxOzGWGTJi8n4/HPB5w5eTkTGYSLs55v68T8jzPQ0RERORjwsk+ABEREel7NEEQERERH00QRERExEcTBBEREfHRBEFERER8NEEQERERH00QRERExEcTBBEREfHRBEFERER8NEEQEWcnnXQSJ510UkL3+e677xIKhVi8eHFC9ysin44mCNKnLV68mFAoFHtkZ2dz6KGHMmfOHKqqqnq0v5deesn3cc/zGD58OKFQiC9/+ctxH/v4caSnp1NcXMyECRO48soreeutt7r8unfeeSehUIjJkyfHjR900EFx+w16fPwfzaeeeopQKMSwYcOIRqPO34NkW7p0KfPnz0/2YYjIJ0hP9gGIdMeNN97IqFGjaG5u5qWXXmLhwoU89dRTbNiwgdzcXOf9ZWdns3TpUo4//vi48dWrV/P++++TlZVlft6XvvQlLrjgAjzPo6amhr/+9a/cd9993HnnnfzsZz+joqLC/LwlS5Zw0EEH8eqrr7J582YOPvhgAObPn099fX1su6eeeooHHniA2267jUGDBsXGv/CFL/j29e677/L8889TXl7u/PyTaenSpWzYsIGrrroqbnzkyJE0NTWRkZGRnAMTkXieSB927733eoD32muvxY1XVFR4gLd06dJPtb+zzjrLGzRokNfW1hb38UsvvdSbMGGCN3LkSO/000+P+xjgzZ4927fP3bt3e2VlZR7gPfnkk76Pv/POOx7gPfLII97gwYO966+/PvD4brnlFg/wtmzZYn68vr7ey8vL8+644w5v/Pjx3oUXXtiNZ514J554onfiiSd+qs89/fTTvZEjRyb0eEQk8XSJQfZLp5xyCgBbtmzhnXfeIRQKcdttt/m2e/nllwmFQjzwwANx4zNmzGDPnj2sWLEiNtba2sof/vAHvva1rzkdy8CBA1m2bBnp6en89Kc/9X18yZIlHHDAAZx++ul89atfZcmSJU77/7hHH32UpqYmzj77bM477zweeeQRmpubu/W5mzZtYvr06ZSWlpKdnc2BBx7IeeedR01NTWyb9vZ2fvzjH/O5z32OrKwsDjroIH7wgx/Q0tLS5b73Xrp5991348ZXrVpFKBRi1apVQEd24cknn+S9996LXT456KCDgOAMwvPPP88JJ5xAXl4eRUVFnHHGGfzjH/+I2+b6668nFAqxefNmLrzwQoqKiigsLGTWrFk0NjZ26/sjIvE0QZD90ttvvw10/OP82c9+luOOO878h3fJkiUMGDCAM844I278oIMOoqysLG7i8PTTT1NTU8N5553nfDwjRozgxBNPZO3atdTW1vqO4ayzziIzM5MZM2awadMmXnvtNeevsXdfJ598MqWlpZx33nnU1dXxxBNPfOLntba2MmXKFNauXcsVV1zBggULuOyyy3jnnXeorq6ObXfJJZdw7bXXcswxx3Dbbbdx4oknMm/evE/1PbH893//N+PGjWPQoEH87ne/43e/+12XeYTnnnuOKVOmsHPnTq6//noqKip4+eWXOe6443yTEYBzzjmHuro65s2bxznnnMPixYu54YYbEnLsIv1Osk9hiHRl7yWB5557ztu1a5e3bds2b9myZd7AgQO9nJwc7/333/c8z/N+/etfe4D3j3/8I/a5ra2t3qBBg7yZM2f69vfaa695v/rVr7wBAwZ4jY2Nnud53tlnn+2dfPLJnud5TpcY9rryyis9wPvrX/8aG3v99dc9wFuxYoXneZ4XjUa9Aw880LvyyivNfXR1iaGqqspLT0/37rrrrtjYF77wBe+MM84IPKa9/vKXv3iA9/DDDwdus379eg/wLrnkkrjx7373ux7gPf/887GxzpcY9n5fOx/3Cy+84AHeCy+8EBsLusSwZcsWD/Duvffe2Ni4ceO8IUOGeHv27ImN/fWvf/XC4bB3wQUXxMauu+46D/AuuuiiuH1+5Stf8QYOHBj4nEUkmM4gyH6hvLycwYMHM3z4cM477zzy8/N59NFH+cxnPgN0/M8xOzs77izCM888w+7du/n6179u7vOcc86hqamJP//5z9TV1fHnP//Z+fLCx+Xn5wNQV1cXG1uyZAklJSWcfPLJQEcT4txzz2XZsmVEIhGn/S9btoxwOMz06dNjYzNmzODpp5/mX//6V5efW1hYCHR8T4JOuT/11FMAvqDld77zHQCefPJJp+PtqR07drB+/XouvPBCiouLY+NjxozhS1/6Uux4P+6b3/xm3N9POOEE9uzZ4zurIyKfTBME2S8sWLCAFStW8MILL/DWW2/xzjvvMGXKlNjHi4qKmDZtGkuXLo2NLVmyhM985jOxvEJngwcPpry8nKVLl/LII48QiUT46le/+qmPcW8bYcCAAQBEIhGWLVvGySefzJYtW9i8eTObN29m8uTJVFVVsXLlSqf9//73v2fSpEns2bMntq/x48fT2trKww8/3OXnjho1ioqKCn77298yaNAgpkyZwoIFC+LyB++99x7hcDjWsNirtLSUoqIi3nvvPafj7am9X2/06NG+jx1++OHs3r2bhoaGuPERI0bE/f2AAw4A+MQJlIj4qeYo+4VJkyYxceLELre54IILePjhh3n55Zc5+uijefzxx/nWt75FOBw8D/7a177GpZdeSmVlJVOnTqWoqOhTH+OGDRtIS0tj1KhRQEe4bseOHSxbtoxly5b5tl+yZAmnnnpqt/b98dzCIYccYu7rsssu63Ifv/jFL7jwwgv505/+xLPPPsu3v/1t5s2bx9q1aznwwANj24VCoW4d08cFfY7rWZKeSktLM8c9z+vV4xBJBZogSMo47bTTGDx4MEuWLGHy5Mk0NjbyjW98o8vP+cpXvsL/9//9f6xdu5YHH3zwU3/trVu3snr1asrKymJnEJYsWcKQIUNYsGCBb/tHHnmERx99lEWLFpGTk/OJ+1+yZAkZGRn87ne/8/0j+NJLL3HHHXewdetW3/+gOzv66KM5+uij+eEPfxgL+y1atIif/OQnjBw5kmg0yqZNmzj88MNjn1NVVUV1dTUjR44M3O/e/6l/PPAImGcdujsB2fv1Nm7c6PvY//7v/zJo0CDy8vK6tS8RcacJgqSM9PR0ZsyYwdKlS/nHP/7B0UcfzZgxY7r8nPz8fBYuXMi7777LtGnTPtXX/fDDD5kxYwaRSIT//u//BqCpqYlHHnmEs88+27xsMWzYMB544AEef/xxzj333E/8GkuWLOGEE04wty0rK+OOO+7ggQce4Hvf+575+bW1teTm5pKe/tGP/NFHH004HI5VGP/zP/+TH/zgB8yfP59f//rXse1uvfVWAE4//fTA4/vc5z4HwIsvvsi4ceOAjrMHv/nNb3zb5uXlxV3aCDJ06FDGjRvHfffdx9y5c2NndzZs2MCzzz4bmC0RkcTQBEFSygUXXMAdd9zBCy+8wM9+9rNufc7MmTO7vf9//vOf/P73v8fzPGpra/nrX//Kww8/TH19PbfeeiunnXYaAI8//jh1dXX8n//zf8z9HHvssbGzHZ80QXjllVfYvHkzc+bMMT/+mc98hmOOOYYlS5YEThCef/555syZw9lnn82hhx5Ke3t77GzE3tDj2LFjmTlzJr/5zW+orq7mxBNP5NVXX+W+++7jzDPPjAUtLUceeSTHHnssc+fO5cMPP6S4uJhly5bR3t7u23bChAk8+OCDVFRU8PnPf578/PzAydktt9zC1KlTKSsr4+KLL6apqYlf/vKXFBYWcv3113f5fRORHkp2jUKkK0ErKXblyCOP9MLhcKwC+Wn2F1Rz3PsIh8NeUVGRN378eO/KK6/0/v73v8dtO23aNC87O9traGgI/BoXXnihl5GR4e3evTs2ZtUcr7jiCg/w3n777cB9XX/99b6K5ce988473kUXXeR97nOf87Kzs73i4mLv5JNP9p577rm47dra2rwbbrjBGzVqlJeRkeENHz7cmzt3rtfc3By3nbWS4ttvv+2Vl5d7WVlZXklJifeDH/zAW7Fiha/mWF9f733ta1/zioqKPCBWebRqjp7nec8995x33HHHeTk5OV5BQYE3bdo076233orbZm/NcdeuXXHjQfVLEflkIc9TekdSy/jx4ykuLnZuCYiIyEdUc5SU8vrrr7N+/XouuOCCZB+KiMh+TWcQJCVs2LCBdevW8Ytf/ILdu3fzzjvvkJ2dnezDEhHZb+kMgqSEP/zhD8yaNYu2tjYeeOABTQ5ERHpIZxBERETER2cQRERExEcTBBEREfHZLxZKWrBgAbfccguVlZWMHTuWX/7yl0yaNKlbnxuNRtm+fTsDBgz4VGvMi4hIcnmeR11dHcOGDevy3io90dzcTGtra0L2lZmZmRo5qCSuwdAty5Yt8zIzM7177rnH+/vf/+5deumlXlFRkVdVVdWtz9+2bVvcAjd66KGHHnrsn49t27btk39nmpqavNIhaQk7ztLSUq+pqWmfHGtv6vMhxcmTJ/P5z3+eX/3qV0DHGYHhw4dzxRVX8P3vf/8TP7+mpoaioiIOvfRa0jLjZ3QZDf6nbo0BZNRH7fEG/1Ky6Q1t5rbhBnt2GmpqNse9phb/YKsxBnjNAePt/rvpedGAl9yzn2NChPyz/lCGfQIrlG7fkS+UlWWP5/hn6l62va2XZ8/qI7kZvrG2fPv42vLs42vL8z/Htnz7rFVbwD2G2vP8r017vv16RXPtOyWm59nvv7w8/3vkgNxGc9shOfXmeElWbbfGAIak2/dbGJLm3744zT6OorD9M1MQ9n9fc0L+1xAgK2S/jmnGe9JVi2d/r5s9/++F+qj9etVG7ffTh1H/Tbx2tw8wt90dyTfHq9oKfWM7mv1jADtb7H3vacw1x2sa/T9LLY2Z5rbU269NWr3/NUhv8L+20ZZm3rn9RqqrqykstI+/J2prayksLOS9dQdRMKBn74vauigjJ7xLTU0NBQUFCTrC5OjTlxhaW1tZt24dc+fOjY2Fw2HKy8tZs2aN+TktLS2xm88A1NXVAZCWmU1aVvwbOq3N/4s3rdX+ZZyeYf/jmZ5uTBACbjkbTrP/sQiF7a/pWeOhgG0Dx/3HF7Qt9PIEIeAXd+B42P7lEwr7JwNeWsAEIWA8lO7ft5duH4eXYb++0Uz/c4xm2q951D4Motn+1yZsjAGQY/+DE861jy/N+D2fnmfvIyPH/oc5K9v/iz47y/4+5QR8//KMCWB+mv1LOT/gdPIAY4KQG/APflbAeGImCPY+Mo3xUMDEPBq199FijDe2269tTrv9vc5q879emcZ7HSAjYDwtZL9Z0/BPEMJewAQhYk8Qwu3+55jWHnwpeF9fJs4fECJ/QM++RpTUuZTdp0OKu3fvJhKJUFJSEjdeUlJCZWWl+Tnz5s2jsLAw9hg+fHhvHKqIiOznIl40IY9U0acnCJ/G3LlzqampiT22bduW7EMSEZH9QBQvIY9U0acvMQwaNIi0tDSqqqrixquqqigtLTU/Jysri6yAa9UiIiLSPX16gpCZmcmECRNYuXIlZ555JtARUly5ciVz5sxx2ldGg+fLHGTU+2d6mUFhxHr/tXyA9Hp/SCncYAcGQ41BYUR7nBb/foLDiPbxmYHERJwCC7h+GzKuDQNg5DJCmQGZgqyA8YDakJdjZBCCwoh59r6tQGJbvn29tzXffu5WILHNzo7RFhA8tAKJ0YCcQHq+HZDLz7PfT8V5/iBgSU6due3QbDtgWJrlHy8JCCOWBowPTGvwjQWFEQsD3k+5RiAxUWFE6xRxO/ZrYIURAeqMQGJNQBhxjxFGBNjV7g+47QoIKe5oKzLHtzf7x6ua7X3sarSTszUN9vE1Nxg/S0FhxDr7Ncio97++Gf63B2H7117CRYn2OInV8z30HX16ggBQUVHBzJkzmThxIpMmTWL+/Pk0NDQwa9asZB+aiIikkIjnEelhsa+nn9+X9PkJwrnnnsuuXbu49tprqaysZNy4cSxfvtwXXBQREZHE6fMTBIA5c+Y4X1IQERFxkYiQoUKKIiIiKSaKR0QThJh+M0HIaPB8iyBZgUSXMCLYgcREhBHBDiQ6hRE7PmCPW1yCh0HbBix4Q4Y/vJSIMCLYgUSXMGLHuD9A5hJG7Bi3xrofRgQ7kJiIMCLYgUSXMCLYgUSXMCLYgUSXMCLYgcREhBHBDiS6hBHBDiS6hBE7xv1hQpcwItiBxISEEcEMJLqEEcEOJFrh8XDAAnayb/WbCYKIiEhXdIkhniYIIiIiqMXQWcqtpCgiIiI9pzMIIiIidNyurucLJaUOTRBERESASAJaDD39/L6k30wQMuqjvls2W40Fl7YC2I2FRLQVwG4seBE7Ne0iFHA7ams5ZAi4xWrQthkBb6kMfxI6qMXg0lYAu7Hg0lYAu7Hg0lboGPf/YmjP635bAezGQiLaCmA3FlzaCmA3FlzaCmA3FlzaCpCYWzW7LJ/s0lYAqDbu5/1hxH7juCyf7NJWALuxkIi2AtiNBZe2Qsf2/p8Ps8XQ1jv/6Ea8jkdP95EqlEEQERERn35zBkFERKQryiDE0wRBREQEiBIiQsAdaR32kSp0iUFERER8dAZBREQEiHodj57uI1X0mwlCRkM76enx6WSrseDSVoCAxkIC2goQ0FgIWqXLahpgNxZC6UFNAzutbN5fIbDFYO/D2ndgWyHXHne5v4JLW6Fje//3z6WtAHZjIZrf/bYC2I2FRLQVwG4suLQVwG4suLQVwG4s7Mu2Qotnf69d7q/g0lYA2BXx31+hsq3Q3Nbl/goubQWwGwuJaCuA3VhwaSsEjVv3yGlv650r+5EEXGLo6ef3JbrEICIiIj795gyCiIhIV3QGIZ4mCCIiIkDUCxH1ethi6OHn9yW6xCAiIiI+/eYMQnpDG+mdgnVWINEpjAhmIDEhYUSwA4kOYUSwA4mhzICQUmZASNEKNQYGHe1xz9h3NCCMGM21j6NtQMB4nn+e6xJG7Bi3xgLCiAHj1vLJLmFEsAOJiQgjgh1IdAkjgh1IdAkjAmSE/O9V1zBixDOCbA5LJ4Pb8sl7ovbyxLva/WFEgMp2fyBxZ5u9rcvyyS5hRIDmeuNnrCEgEOoQRgQ7kOgSRgQ7kGgtgR8K+N2ZaLrEEK/fTBBERES6EiFMpIcn1nt+t5y+QxMEERERwEtABsFTBkFERERSmc4giIiIoAxCZ5ogiIiIABEvTMTrYQZBSy3vf8INrYTT4md2VmPBpa0AdmMhsK3guki31VgISnsHtRisxkJOdve3BUj379sLaCsEtRii2f7xSI6dcm/PD2grBDQTWo0WQ2BbwV6pljZjmeTAtkLQ8sl5/sbCgPwmc9sDcu1xq7GQiLYC2I0Fl7YC2I0Fl7ZCV+MWq60AdmOhMWo3RhoC9hG0fLLVWHBpKwDsaC3yjVW12PtwWT7Zqa0AUO//uUur73lbASCjzv/zkdHQ/bZCx9f0/660lsAnYr+2sm/1mwmCiIhIV6KEiPYwmhcldU4haIIgIiKCMgidqcUgIiIiPjqDICIiQqJCirrEICIiklI6Mgg9vFlTCl1i6DcThFBTM6Fw/MzObCw4tBXAbiwE3lvBldFYcLnnAgBZ/nRzKNtOPHtZAS0Go5ngZdjHEc0MGM/y76M9zz5m694KYLcVwG4sWPdW6Nh3wP0VjHHr3gpgtxXAvr+CS1sB7MZCItoKYDcWXNoKYDcWskIBa/s73l/B4nJ/hUS0FcBuLOxqt5sGQfdXsBoLLm0FsBsLzQ0BP6NB91cwGgvObYWg+ysYjQWXtgLYjQXrHjnhiP0+lX2r30wQREREuhJNwL0YUqnFoJCiiIgIH2UQevpw8eKLLzJt2jSGDRtGKBTiscce63L7Cy+8kFAo5HsceeSRsW2uv/5638cPO+ww5++HJggiIiJ0nEFIxMNFQ0MDY8eOZcGCBd3a/vbbb2fHjh2xx7Zt2yguLubss8+O2+7II4+M2+6ll15yOi7QJQYREZGkmTp1KlOnTu329oWFhRQWfrSC52OPPca//vUvZs2aFbddeno6paWlPTo2nUEQEREBIl4oIQ+A2trauEdLQAC+p+6++27Ky8sZOXJk3PimTZsYNmwYn/3sZzn//PPZunWr8777zRkEr6kFr1OLwWosuLQVIKCxENSDte6tAIH3V7AaC6GAex0ENROscS+oxZBjJ6StxkI0wz7moBZDe45/vD03oJUQ1FbIs79/7UYI3GoldDVuNRbS8rvfVgAozmv0jQ3OqTe3HZJttxiGZPrHg9oKQ9LtfbjcX8GlrQB2YyERbYU2z24rBI03Rv3jdVH7OKqj9vv9w4hdddljjFe12fdcqAy4v8KuZv8+9jTlmtvWNtr3Rmlu9P88egFthfSg+ys0+F/fdP/btGPboLZCwLjVWAhqK6Q12D9LVmPBukdOKLpv/nHtLJKAkGLk3yHF4cOHx41fd911XH/99T3ad2fbt2/n6aefZunSpXHjkydPZvHixYwePZodO3Zwww03cMIJJ7BhwwYGDAi4IY2h30wQREREesu2bdsoKPhoApllVM576r777qOoqIgzzzwzbvzjlyzGjBnD5MmTGTlyJA899BAXX3xxt/evCYKIiAgQ9cJEe7iSYvTfZ5ALCgriJgiJ5nke99xzD9/4xjfIDLoT778VFRVx6KGHsnnzZqevoQyCiIgIH11i6OmjN6xevZrNmzd364xAfX09b7/9NkOHDnX6GpogiIiIJEl9fT3r169n/fr1AGzZsoX169fHQoVz587lggsu8H3e3XffzeTJkznqqKN8H/vud7/L6tWreffdd3n55Zf5yle+QlpaGjNmzHA6tv5ziaG1BUKdllo2AolOYUSwA4kJCCOCHUgMBVzHCmXbQScv1z8ezbf3Ec22g2nRTP9xB4UUI8a2AO05/u9Je07AksrG0skA7UHLJ+f7XwNrDIKXT04b4A9R5eXawSgrjAh2ILE0u9bcdmimHTwcmvEv31hQGHFwmj3usnyySxgR9l0gscWzQ2x1UfvnsdoIJH4YtX8GrNAhwE5jSWWAHa1FvrHtLXZIsaopIKTY5E/OVhtLJwM0NQSEhuv9r0F6nf27Imj55HQjr5pZ1/MwYsf2/tfGWjoZ7DAi2IFEawl8L9o7Sy1HIdZC6Mk+XLz++uucfPLJsb9XVFQAMHPmTBYvXsyOHTt8DYSamhr++Mc/cvvtt5v7fP/995kxYwZ79uxh8ODBHH/88axdu5bBgwc7HVv/mSCIiIh04dMsdGTtw8VJJ52E18UdIBcvXuwbKywspLExoI4CLFu2zOkYgugSg4iIiPjoDIKIiAh8qnspWPtIFZogiIiIAFFCROlpBqFnn9+XJHWq80l3sfI8j2uvvZahQ4eSk5NDeXk5mzZtSs7BiohISkvG3Rz7sqSeQdh7F6uLLrqIs846y/fxn//859xxxx3cd999jBo1ih/96EdMmTKFt956i+yA1H4Qr7kFr3OLwWgseFHHe3kbjYWgVgJBbYX07i+fHNhWyLMT0tEB/u3b8+1FNazlkAGimf7nGM2wZ8kRY1uA9myjxZBrb9sW1FYIWj7ZaCxE8+22Qka+nYa2lk8eGNRWyLaXTx6a7W8mlGbZbYVhGdX2vtP8rYegpZOLA9oKAxyWT05EWyHi2bntduzXoNnz/9zVGUsnA9RE7ffknqj//b4roJVQ2W43EKy2AkCVsXzyjoC2QtDyydX1/vGWhoDFbIy2AkCGsXxyekBbIcN+S5rNhMzAJZUD2goNdpPEaiy4tBXAbixgjXm902KQeEmdIHR1FyvP85g/fz4//OEPOeOMMwC4//77KSkp4bHHHuO8887rzUMVEZEUl5h7MaTOGYQ++0y2bNlCZWUl5eXlsbHCwkImT57MmjVrAj+vpaXFdxctERGRTxL1Qgl5pIo+O0GorKwEoKSkJG68pKQk9jHLvHnzYvfLLiws9N1RS0RERD5Zn50gfFpz586lpqYm9ti2bVuyD0lERPYD0QTch6GnCy31JX225lhaWgpAVVVV3A0mqqqqGDduXODnZWVl7ZPbaoqISGpLzN0cNUHY50aNGkVpaSkrV66MTQhqa2t55ZVXuPzyy53357VH8ELxaVyzsRCQyA5iNRYCWwlBt+QMvL+Cf9y6twLYbQWA9gH+r9laYB9fe3bQ/RWsFoO5KZGApxgx78Vgbxt0H4X2fPu18Yz7K2QGtBUKjLYCwMBcf1NgSI59r4OhgfdXqPaNlabbLYbB6fY+Bob9zYnCsL2+vUtbASAj5H+vut5bwbqPgjUG0Bh4fwX/61sT8IbaE7VbAlZjwbWtsKPZbibsbBrgP45G/70VAGob7J+71nr/D0KowW5kpBttBYCMOv/rm2EXWgLvo5DRYLUYAho+xr0VANIaun9/hcC2gnHfG8BsLHit/p9dTy2GpEjqBKG+vj7u/tR772JVXFzMiBEjuOqqq/jJT37CIYccEqs5Dhs2jDPPPDN5By0iIikpQohIDxc66unn9yVJnSB80l2srrnmGhoaGrjsssuorq7m+OOPZ/ny5c5rIIiIiHwSXWKIl9QJwifdxSoUCnHjjTdy44039uJRiYiISJ/NIIiIiPSmCD2/RGAnPPZP/WaC4EU931LLLoHEoOWTrUBiYBgxx740YoURATxjPJpvbxu0fLIVSGzLs0+BtRlBQoBIVvdDitGgkKLx1NuDlk7ODQgjBiyfnJnnDzAV5jeZ2xbn2MsnD831hwaHZAWEFI0wItiBRJcwItiBxESEEYPGg5ZJjmK/NlYg0SWMCHYgcWfEXl97V8QOEla1+QOJO1oDQorN9viuJvtr7m70ByPrAsKIbQHLJ4eM5ZMzApZJdlk+2QodQlfLJ/tfr6Clk9Pq7SBguNEetwKJLmFECAgkWkvgBwRhE02XGOL1mwmCiIhIV3S753ip80xEREQkYXQGQUREBPAIEe1hBsFTzVFERCS16BJDvNR5JiIiIpIw/ecMghcFutFaCFp6NqDFQIaRJM+00+VB7QYvK2A8xz8ezbb33Z5jH5+1fHJQW6E91x63mgnuSy37U9aRbDt57eXYr1NGjp2Wz8vxJ6cLs+3UdHGW3R4oNtawHZRhtxiK04x4ecB4UdhuUwxwWD4526GV0NW4i2jAz0qL50+YNwasZVLn2b9eqqP+NbarA5ZU/rDdbhp82O5f+vjDVns55OpWe03vmha7mdDY7G8KtTcH/Ew3BSyf3OR/HdOMsY5tzWHSm/zf14xG+3ud3mi/XulN/vR/WqPdYgg32e/JUHPAMsctxniL3WKwmgkAXsR/3NYS+F2tl5NIibhdcyrd7rn/TBBERES6sPeOjD3dR6pInWciIiIiCaMzCCIiIugSQ2eaIIiIiABRwkR7eGK9p5/fl6TOMxEREZGE0RmETkIB696HQgHjaf45lnV/BgDSg5oQ9vZehn/7aKY9p4tmBtxHwRi37q3QsQ/78KxmQvC9GOy0cTTDPx7NspPXoSx73fWMTDsJnZPpT1/npdvJ6wEZdruhIN0/PiAcsG2aPZ4X9ie4c0P2MQcURsgwWjT7tq3Q/XsuALQZ27cFBMybPfv4mj3/m6cxat9jpDHgTVnb7m8gNARUaBra7PGmVvtN3NbmP26vxf65S2uxX8iw8fZLCygDWNsGbZ/Wan+z01rtn6Vwi/91DLXa70naAsZb7XaD1+Yf9wL24UXs95M5bt0fxOG+OT0R8UJEeniJoKef35dogiAiIoIyCJ1pgiAiIgJ4Cbibo6eVFEVERCSV6QyCiIgIECFEpIc3W+rp5/cl/XuCYC2r7LrUsjUeEFL0HMKIANEM/7FYYx3jAcFDI4cVGDB0GLdCh13uwwovZtjBo7R0ezwrIyDwl+EPS+UGhBTz0+ylYHON8QEBYcTckL2PbCOQmB0KWDY6IPSagf+9kG6MJUrQksqRgPBiq7HkrUsYEaDBCCRaYxAcUmwyAomN7fa2ze32z11be8DPXat/PNRm/9yF2oJCiv7x4DBiUPDQPx4O2DbcFhBSNMZDbXZgMBQUMAxYJtkMNQaEEQPHzUCi8Rx7banlnmcIjJWi91u6xCAiIiI+/fsMgoiIyL9FExBS7Onn9yWp80xERER6IEooIQ8XL774ItOmTWPYsGGEQiEee+yxLrdftWoVoVDI96isrIzbbsGCBRx00EFkZ2czefJkXn31VddvhyYIIiIiydLQ0MDYsWNZsGCB0+dt3LiRHTt2xB5DhgyJfezBBx+koqKC6667jjfeeIOxY8cyZcoUdu7c6fQ1dIlBRESE5KykOHXqVKZOner8dYYMGUJRUZH5sVtvvZVLL72UWbNmAbBo0SKefPJJ7rnnHr7//e93+2v06wmCtayytXQyQCiggRDKMJLaAdsGjUczA9LUxngkYKlla0nljnFrv/bhBaxUazYQAtsKWXaE1zP2Ecq0k9eZWfbSrtaSymAvq+yypDLYyypbSyd3NW4tq5ztsKRyx7j/NU8LatY4ihiJ8UhAOrwtYGnbFmPzoBZD0PLJDcYbsD7iXzoZoKY9xxyvazeaEAFLKjcHLKnc2hrQKmr1f7+tVgIEL59sjbssqdwxbrQYgtoKxpLKELCscsDSyc5LLRvtBi9iH58XFO3vpXZCdyUyg1BbWxs3npWVRVaW/TPxaYwbN46WlhaOOuoorr/+eo477jgAWltbWbduHXPnzo1tGw6HKS8vZ82aNU5fQ5cYREREEmz48OEUFhbGHvPmzUvIfocOHcqiRYv44x//yB//+EeGDx/OSSedxBtvvAHA7t27iUQilJSUxH1eSUmJL6fwSfr1GQQREZG9oiTgXgz/Dilu27aNgoKC2Hiizh6MHj2a0aNHx/7+hS98gbfffpvbbruN3/3udwn5GntpgiAiIgJ4n6KFYO0DoKCgIG6CsC9NmjSJl156CYBBgwaRlpZGVVVV3DZVVVWUlpY67VeXGERERPjobo49ffS29evXM3ToUAAyMzOZMGECK1eu/Oh5RaOsXLmSsrIyp/3qDIKIiEiS1NfXs3nz5tjft2zZwvr16ykuLmbEiBHMnTuXDz74gPvvvx+A+fPnM2rUKI488kiam5v57W9/y/PPP8+zzz4b20dFRQUzZ85k4sSJTJo0ifnz59PQ0BBrNXRX/5kghML++yxY91GwWgkAGQERf2N7L9PeRzQ7oMWQZY+35/iPrz3Hnp22B8TlI8b2AYFxIjlB91cwWgzWvRWw2woAoWx/yjoz205H5wa0GAZk2u2BA7IafWNFGU3mtoXp/m0BitKMfYTtbQcExNHzwv4Ed3ZAAyE7ZL/m4QTc6MVqKwC0438NWjw7ud4YkC5v8PzHXRe131C1AeMfRvL9Y+155rY1bXaLobrFP17bbH+9xmb7ZzfSHNBMava/ZmnNAS2GpoBxoyyT3mx/T4PG05r9r2N6o91WCLfYr2Oo2f+zFGoJaCW02D9fQfdi8Iz7K1hjHR+w35Mm8z4lIQJuD5JQyVhJ8fXXX+fkk0+O/b2iogKAmTNnsnjxYnbs2MHWrVtjH29tbeU73/kOH3zwAbm5uYwZM4bnnnsubh/nnnsuu3bt4tprr6WyspJx48axfPlyX3Dxk/SfCYKIiEgXEnGJwPXzTzrpJLwu6p6LFy+O+/s111zDNddc84n7nTNnDnPmzHE6ls6UQRAREREfnUEQERGBT3UvBWsfqUITBBEREZJziaEv6zcThFBGOqFOwbBQpj+8FMqyA01B416Of/GLaK69IEYkxw4vtucFhBRz/VeA2nPsq0LtuQHhRSPj1Z5nX++KZAeEFLOMgFGGHToKXD7ZCCTm5dihqKIcO2A4KLvBHs+s942VZNQaW8Lg9DpzvDjNv4+isH0cA4wllQFyjXBVVkAYMR17eWKXZZVdwogAzUYgsSFgH3VR+/iqjeDhHiN0CLCr3e6A724b4B9rCdhHsz1e3ex/Y9c32T93bU32z12o0X6O6Y3+1yC9wf75SrffIqQ3+n+WMoyxjm3t18AKJKY12QHDcKMdnA01+3/GPGMM6GIJ5oBxK5AYFEYMusZuBRLNn4Fwr4QUJV6/mSCIiIh0RWcQ4mmCICIigiYInanFICIiIj46gyAiIoLOIHSmCYKIiAgdOcie36wpdfSbCUIoPc3fYjCaCaFse7lWq60A4BmNhWhuQFsh3x5vy7Ov9Fjjrfn2m7fNDnvTlu9/u7bn2kljLyegmZBlpKnTA9oKAcskW8snB7UVBmf7GwUAJVl2M2FIpn98UHpAiyHNHh+Y5m9IFIYDlnwO269Bbsj/+ga1GHq7rQBQF/VvXxPQVtgTtZc4tpoJle2F5rY72+wWQ1WLf3xXQIthT1OuOV7b6P85bWkIWA69IeA1qLdfg/R6/+ubYa+6TUa9/c9BRoPRYqi3X8eMBvv1Sm/wv/8C2wpNAcskW42FoCWVW+19uyy17NRWALOxEMrwv14hz4OA8kUi6QxCPGUQRERExKffnEEQERHpis4gxNMEQUREBE0QOkvqJYZ58+bx+c9/ngEDBjBkyBDOPPNMNm7cGLdNc3Mzs2fPZuDAgeTn5zN9+nSqqqqSdMQiIiL9Q1InCKtXr2b27NmsXbuWFStW0NbWxqmnnkpDw0eBsauvvponnniChx9+mNWrV7N9+3bOOuusJB61iIikor1nEHr6SBVJvcSwfPnyuL8vXryYIUOGsG7dOv7jP/6Dmpoa7r77bpYuXcopp5wCwL333svhhx/O2rVrOfbYY7v9tUJZWYTC8Slnq7EQ2FbIs9sNkTx/crptQEBbId+ej7UGthj8b7T2oLZCwP0V2vP9yWkv306/Z+TYqf2MTH+KOSvDTjbnZAYk/zP9EeSgeysEtRVKs2rs7dP946XGGNhtBYCisD/B7dJWALux4NJWALuxkIi2AtiNBZe2AtiNhaC2wvbmInO8qtl/L4ZdjXnmtjUN9vE11xs/p0FthTr7Ncgw2goAGcZbxKWtAJBpNBYC2wr1AfdXaPD/zIQam81tvSZ7HGM8IW2FIA5tBQhoLKRbLYZor7QYPC+E18N/4Hv6+X1Jn2ox1NR0/FIvLi4GYN26dbS1tVFeXh7b5rDDDmPEiBGsWbMmKccoIiLSH/SZkGI0GuWqq67iuOOO46ijjgKgsrKSzMxMioqK4rYtKSmhsrLS3E9LSwstH+v51tba/xsVERH5uCihHi+U1NPP70v6zBmE2bNns2HDBpYtW9aj/cybN4/CwsLYY/jw4Qk6QhERSWXKIMTrExOEOXPm8Oc//5kXXniBAw88MDZeWlpKa2sr1dXVcdtXVVVRWlpq7mvu3LnU1NTEHtu2bduXhy4iIpKSkjpB8DyPOXPm8Oijj/L8888zatSouI9PmDCBjIwMVq5cGRvbuHEjW7dupayszNxnVlYWBQUFcQ8REZFPsjek2NNHqkhqBmH27NksXbqUP/3pTwwYMCCWKygsLCQnJ4fCwkIuvvhiKioqKC4upqCggCuuuIKysjKnBgNAKCebUDg++Ww1FlzaCgBt+f5vYdC9FQLbCgH3V2g3gt3WvRUA2gPGvTx/Ajkzz04x5+XYMWGrmZCbYSev89LtfR+Q5V/MflCmfc8F694KYLcVwG4suLQVAAqNxoJLWwH23f0VEtFWALux4NJWALux4NJWALux4NRWALOxkIi2Qsf21n0Uut9W6Nje/5q5tBXAbiwEthUc7q+QkLZCEIe2AgQ0FrL9r3koGoKA1yuRtFBSvKROEBYuXAjASSedFDd+7733cuGFFwJw2223EQ6HmT59Oi0tLUyZMoU777yzl49URERSnWqO8ZI6QfCC7vz1MdnZ2SxYsIAFCxb0whGJiIgI9KGao4iISDJ5CbjEoDMIIiIiKcYDunFi+xP3kSr6zQTBy87CS+sUUjQCiS5hxI5xfyCsNWBJ5aAwYlvA8sntxvLJQSHFaMDyyZn5/pBSYX6TuW1hth2AsoKHuQFhxAEZ9j6KMvxfsyTDDiMOSg9Yatlh+eTigDCiy/LJvR1GBDuQmIgwItiBRJcwItiBRJcwItiBRJcwItiBROcwYl33l092CSOCHUh0CSNCQCAxKIzYHDBuBBITEkYEM5DoEkaEgEBijv/9EYr2iUZ+v9NvJggiIiJdiRIipJUUYzRBEBERQS2GznTeRkRERHx0BkFERISORY5CWigpRhMEEREROhoMPW4xpFCNod9MELy8bF+LwWosuLQVwG4sBLYV7LA3bUZbAewWQ9RYOhkgw2grABTk+ZPQxTn+ZY8Bio3lkMFuJuSn2anpgnQ7kV2Y7t/34PQ6c9vBaXaLwWX5ZJe2AtiNBZe2AtiNBZe2AtiNhUS0FcBuLLi0FcBuLLi0FSCgseDQVgC7sZCItgJAZp3/dcxo6H5bAezGglNbAczGgktbAXp/+WSXtgLYjQVzCfwEFS/ETb+ZIIiIiHRFIcV4miCIiIigCUJnajGIiIjw0d0ce/pw8eKLLzJt2jSGDRtGKBTiscce63L7Rx55hC996UsMHjyYgoICysrKeOaZZ+K2uf766wmFQnGPww47zPXboQmCiIhIsjQ0NDB27Nhu35DwxRdf5Etf+hJPPfUU69at4+STT2batGn85S9/idvuyCOPZMeOHbHHSy+95HxsusQgIiJCcloMU6dOZerUqd3efv78+XF/v+mmm/jTn/7EE088wfjx42Pj6enplJaWuh1MJ/1mghDJzSCUHt9asBoLLm2Fju39p5OC7q3g0lYAu7GQNsBOTecbbQWAgbn+aPfQXLslUBwQA7eaCbkBLYYB4YB7MaT5WwzFafXmti5tBYBCo7Hg0laAfXd/BZe2AtiNhUS0FcBuLLi0FcBuLDi1FcBsLLi0FcBuLCSirdCxvf81c2krgN1YcGkrgN1Y6CttBbAbCy5tBQhoLBj3yPEivXNdv2OC0NMMQoIOppui0Sh1dXUUFxfHjW/atIlhw4aRnZ1NWVkZ8+bNY8SIEU777jcTBBERkd5SWxv/H7GsrCyysgImzT3wP//zP9TX13POOefExiZPnszixYsZPXo0O3bs4IYbbuCEE05gw4YNDBgQ0Lc3KIMgIiLCRy2Gnj4Ahg8fTmFhYewxb968hB/v0qVLueGGG3jooYcYMmRIbHzq1KmcffbZjBkzhilTpvDUU09RXV3NQw895LR/nUEQEREBvH8/eroPgG3btlFQ8NFlvUSfPVi2bBmXXHIJDz/8MOXl5V1uW1RUxKGHHsrmzZudvobOIIiIiCRYQUFB3CORE4QHHniAWbNm8cADD3D66ad/4vb19fW8/fbbDB061Onr9JszCG356XidQjVWINEljNgxbo0FhBEDxoOWT07L9wej8nLtQNPAPHuZ5CE5/uWMh2TZSxwPyrDHreDhgDQ7cJUXto+vKOw/vqJwk7ltYdgOhLksn9zbYUSwA4kuYUSwA4mJCCOCHUh0CSOCHUh0CSOCHUh0CSOCHUhMRBgR7ECiSxgRAgKJDmFEsAOJfSWMCHYg0SWMCHYg0VoCP9Juv4aJloyFkurr6+P+Z79lyxbWr19PcXExI0aMYO7cuXzwwQfcf//9QMdlhZkzZ3L77bczefJkKisrAcjJyaGwsON3wne/+12mTZvGyJEj2b59O9dddx1paWnMmDHD6dh0BkFERAQ+usbQ04eD119/nfHjx8cqihUVFYwfP55rr70WgB07drB169bY9r/5zW9ob29n9uzZDB06NPa48sorY9u8//77zJgxg9GjR3POOecwcOBA1q5dy+DBg52Ord+cQRAREelSAs4g4Pj5J510El4X3cjFixfH/X3VqlWfuM9ly5Y5HUMQnUEQERERH51BEBERITkrKfZlmiCIiIiguzl21m8mCG15aXgZ8clxq7Hg0lboGPdPFwPbCvl2Ajk9r/vLJxcHtBUGZ9vLFg/N9i+rPDSz2tw2aOnjAqOxkBuyk9dBLYYBxjLJA0J2itylrQB2Y8GlrQB2Y8GlrQB2Y8GlrQB2YyERbQWwGwsubQUIaCw4tBXAbiy4tBU6tvePJ6KtAHZjwamtAGZjwaWtAPuusZCItgLYjQWXtgLYjQVrCfz29n7zT1Wfou+6iIgIdAQMezmk2JdpgiAiIoIyCJ2pxSAiIiI+OoMgIiICib0ZQwrQBEFERAS1GDrrNxOEtrww0cz4KypWY8GlrQDQnucfD7q3gktbAezGwuCcoLZCjT1uNBZK0+1tg1oMVjMhO6CBkBswnhf2J8xzQz1vK8C+u7+CS1sB7MaCS1sB7MZCItoKYDcWnNoKYDYWXNoKYDcWXNoKYDcWEtFWALux4NJWALux0NttBbAbC4loK4DdWHBpK4DdWLDukdPeZreBZN/qNxMEERGRT5RClwh6ShMEERERdImhM00QREREQCHFTlRzFBEREZ9+cwahLT9ENDPUaczarvthRLCXT07Pt0NRA/KbzPEDcu1xK5BYaiydDFCaZQcPrUDi4HR7H0Vh+zis4GF2yA77ZQecXcs2QlRBocPeDiOCHUh0CSOCHUh0CSOCHUhMRBgR7ECiSxgR7ECiSxgR7ECiSxixY3v/65WIMCIEBBIdwohgBxJ7O4wIdiAxEWFEsAOJLmHEjnH/z5K1BH6ktbf+Lxv696On+0gN/WaCICIi0iVdYoijSwwiIiLiozMIIiIioDMInWiCICIiArqbYye6xCAiIiI+/eYMQlseRDuFca3GQntAiyFw+WSjsRC0dHJQW6Ekp84cH5LtHx+aabcVhmVUm+NWY2Fg2L+EM8CAsJ0CzzUmxBkByyRnBKSss41mQjp2G8ClrQB2Y8GlrQB2Y8GlrQB2Y8GlrQB2YyERbQUIaCw4tBXAbiy4tBU6tvePu7QVwG4sJKStAGZjwaWtAPuuseDSVgC7sZCItgLYjQWXtgLYjQVrCfxIa+/8r1y3e47XbyYIIiIiXVIGIY4uMYiIiIiPziCIiIiAQoqdJPUMwsKFCxkzZgwFBQUUFBRQVlbG008/Hft4c3Mzs2fPZuDAgeTn5zN9+nSqqqqSeMQiIpKqQl5iHqkiqROEAw88kJtvvpl169bx+uuvc8opp3DGGWfw97//HYCrr76aJ554gocffpjVq1ezfft2zjrrrGQesoiIpCovQY8UkdRLDNOmTYv7+09/+lMWLlzI2rVrOfDAA7n77rtZunQpp5xyCgD33nsvhx9+OGvXruXYY491+lrteR7R7PhXzmosuLQVwG4sFOfZLYGgtsLQbLuZMCTTaDFk/MvcdnCafX8Fq7FQGNBWGBDufjMhI6CBkBGyx8PG+uSJaCuA3VhwaSuA3VhwaSuA3VhwaSuA3VhISFsBzMaCS1sB7MaCS1sB7MaCS1sB7MZCItoKYDcWerutAHZjwaWtAHZjIRFtBbAbCy5thY7t/e8z6x45Efulkn2sz4QUI5EIy5Yto6GhgbKyMtatW0dbWxvl5eWxbQ477DBGjBjBmjVrknikIiKSkvZmEHr6SBFJDym++eablJWV0dzcTH5+Po8++ihHHHEE69evJzMzk6KiorjtS0pKqKysDNxfS0sLLR/7n0Ftrf0/axERkTiqOcZJ+hmE0aNHs379el555RUuv/xyZs6cyVtvvfWp9zdv3jwKCwtjj+HDhyfwaEVERPqHpE8QMjMzOfjgg5kwYQLz5s1j7Nix3H777ZSWltLa2kp1dXXc9lVVVZSWlgbub+7cudTU1MQe27Zt28fPQEREUoJCinGSfomhs2g0SktLCxMmTCAjI4OVK1cyffp0ADZu3MjWrVspKysL/PysrCyysvwhnPZ8j3CnkKIVSHQJI4IdSHQNI5Zm2eMl6f7xIen2vgem2evdWoHEoDBidkDA0AoeJmqZZItLGBHsQKJLGBHsQKJLGBHsQKJLGBHsQGIiwohgBxJdwohgBxJdwogd2/tfL5cwItiBxESEEcEOJPZ2GBHsQKJLGBHsQGIiwogd4/6fJZcwYse4NWaEx9N76V9dXWKIk9QJwty5c5k6dSojRoygrq6OpUuXsmrVKp555hkKCwu5+OKLqaiooLi4mIKCAq644grKysqcGwwiIiLiJqkThJ07d3LBBRewY8cOCgsLGTNmDM888wxf+tKXALjtttsIh8NMnz6dlpYWpkyZwp133pnMQxYRkVSllRTjJHWCcPfdd3f58ezsbBYsWMCCBQt66YhERKS/SsRKiFpJUURERFKaJggiIiKQlBbDiy++yLRp0xg2bBihUIjHHnvsEz9n1apVHHPMMWRlZXHwwQezePFi3zYLFizgoIMOIjs7m8mTJ/Pqq6+6HRif4hLDzJkzufjii/mP//gP5y+WTNHcCOTEJ5GtxoJLWwHsxkIi2gpgNxYGp9kthuJwqzluNRZyQxnmtkHLJAeN91Qi2gpgNxZc2gpgNxZc2gpgNxZc2gpgNxYS0VYAu7Hg0lbo2N4/7tJWALux4NJWgIDGQgLaCrDvGgsubQWwGwsubQWwGwuJaCuA3VhwaSt0jPvfT+YS+L3VYkiChoYGxo4dy0UXXdStew1t2bKF008/nW9+85ssWbKElStXcskllzB06FCmTJkCwIMPPkhFRQWLFi1i8uTJzJ8/nylTprBx40aGDBnS7WNzPoNQU1NDeXk5hxxyCDfddBMffPCB6y5ERET6nBAJuJuj49ecOnUqP/nJT/jKV77Sre0XLVrEqFGj+MUvfsHhhx/OnDlz+OpXv8ptt90W2+bWW2/l0ksvZdasWRxxxBEsWrSI3Nxc7rnnHqdjc54gPPbYY3zwwQdcfvnlPPjggxx00EFMnTqVP/zhD7S12V1mERGR/qS2tjbu0RJwhsvVmjVr4u5RBDBlypTYPYpaW1tZt25d3DbhcJjy8nLn+xh9qgzC4MGDqaio4K9//SuvvPIKBx98MN/4xjcYNmwYV199NZs2bfo0uxUREUmeBN6safjw4XHL/s+bNy8hh1hZWUlJSUncWElJCbW1tTQ1NbF7924ikYi5TVf3MbL0qOa4Y8cOVqxYwYoVK0hLS+M///M/efPNNzniiCP4+c9/ztVXX92T3YuIiPSeBK6kuG3bNgoKPsorWSv89nXOE4S2tjYef/xx7r33Xp599lnGjBnDVVddxde+9rXYN+PRRx/loosu0gRBRET6pYKCgrgJQqKUlpZSVVUVN1ZVVUVBQQE5OTmkpaWRlpZmbtPVfYwszhOEoUOHEo1GmTFjBq+++irjxo3zbXPyySf7btOcbOl5bYRz49O4VmPBpa0AdmPBta1QGjBu3V+hyKGtAHZjobfbCmA3FhLRVgC7seDSVgC7seDSVgC7seDSVoCAxkIC2gpgNxZc2gpgNxZc2gpgNxac2gpgNhb6SlsB7MaCS1sB7MaCS1sB7MZCItoKHdv732cubQUIaCwY98iJhhPwWnXHfnAvhrKyMp566qm4sRUrVsTuUZSZmcmECRNYuXIlZ555JtBxj6OVK1cyZ84cp6/lPEG47bbbOPvss8nOtt+QAEVFRWzZssV11yIiIkmTjJUU6+vr2bx5c+zvW7ZsYf369RQXFzNixAjmzp3LBx98wP333w/AN7/5TX71q19xzTXXcNFFF/H888/z0EMP8eSTT8b2UVFRwcyZM5k4cSKTJk1i/vz5NDQ0MGvWLKdjc54gfOMb33D9FBERETG8/vrrnHzyybG/V1RUAB1rDi1evJgdO3awdevW2MdHjRrFk08+ydVXX83tt9/OgQceyG9/+9vYGggA5557Lrt27eLaa6+lsrKScePGsXz5cl9w8ZP0uds9i4iIJEUSLjGcdNJJeF7wJ1mrJJ500kn85S9/6XK/c+bMcb6k0JkmCCIiIrBfZBB6k+7FICIiIj795gxCXl4LabnxY1ZjwaWtAHZjIRFtBbAbC4UObQWwmwm93VYAu7GQiLYC2I0Fl7YC2I0Fl7YC2I0Fp7YCmI2FRLQVwG4suLQVOrb3v2YubQWwGwsubQWwGwt9pa0AdmPBpa0AdmPBpa0AdmMhEW2FjnFrrPttBbAbC9Y9cqLh3lmlV7d7jtdvJggiIiJd+thKiD3aR4rQBEFERASUQehEGQQRERHx0RkEERERlEHorN9MEA7IbSS9UyDGCiS6hBHBDiQmIowIdiAxKIyYFQpYjtcIVwUFCYNEjXNmUex9RAL6vC1GILEh4DhcwohgBxJdwohgBxJdwohgBxJdwohgBxITEUbs2N4/7hJGBDuQ6BJGhIBAokMYEexAYkLCiGAGEl3CiGAHEl3CiGAHEl3CiB3j/p+lRIQRO8b97yeXMCLYgURrCfxIKDG3Sv5EusQQR5cYRERExKffnEEQERHpUgIuMaTSGQRNEERERECXGDrRJQYRERHx0RkEERER0BmETvrNBGFITj0ZOfFtAaux4NJWALuxkIi2AtiNBZe2QhCrldAxbifa2zx/AjkSsI+2gGZCo9FuqEtAWwHsxoJLWwHsxoJLWwECGgsObQWwGwuJaCuA3VhwaSuA3VhwaiuA2VhwaStA7y+f7NJWALux4NJWALux4NJWALuxkIi2AtiNBZe2AtiNBWsJ/HZ6p8WgmmM8XWIQERERH00QRERExKffXGIQERHpkjIIcTRBEBERQRmEznSJQURERHz6zRmEkqxasrLjWwFWY8GlrQB2YyERbQWwGwsubQWwGwjWGNj3SwBoM86ZtQbec8E+jgbP/1yqo3Z626WtAHZjwaWtAHZjwamtAGZjwaWtAHZjIRFthY7t/a+vS1sB7MaCS1sB7MZCX2krgN1YcGkrgN1YcGkrgN1YcGkrdGzvf58loq0AdmPBpa0AdmPBukdOW9T+nbpPpNAZgJ7qNxMEERGRLimDEEeXGERERMRHZxBERERQSLEzTRBERERAlxg66TcThJKsWrKz4p+uFUh0CSOCHUhMRBgR3AKJkYAljq1AYqNnB4ms5ZA79uEfa/bssFTQeJ0RSNwTsdNSLmFEsAOJLmFEsAOJLmFEsAOJLmFEsAOJiQgjgh1IdAkjQkAg0SGMCHYgsa+EEcEOJLqEEcEOJLqEETvG/T9LLmHEjnFrrOdhRLADiS5hRLADidYS+C3t9u8s2bf6zQRBRESkK7rEEE8TBBEREdAlhk7UYhAREREfnUEQEREBnUHoRBMEERERlEHorN9MEIak15DTKbVsNRZc2gpgNxaC2goZITvh77p8sqUdO2lsNRbqovY7uM5YDhnsZkKzZz/Hxqid6q41WgxBSyq7tBXAbiy4tBUgoLHg0FYAu7Hg0lbo2N4/noi2AtiNBae2ApiNBZe2Auy7xkIi2gpgNxZc2gpgNxZc2gpgNxZc2god4/73UyLaCmA3FlzaCmA3Fqwl8Jvb7PdSwukMQhxlEERERMSn35xBEBER6ZLOIMTRBEFERARlEDrrM5cYbr75ZkKhEFdddVVsrLm5mdmzZzNw4EDy8/OZPn06VVVVyTtIERGRfqJPTBBee+01fv3rXzNmzJi48auvvponnniChx9+mNWrV7N9+3bOOuusJB2liIikNC9BD0cLFizgoIMOIjs7m8mTJ/Pqq68GbnvSSScRCoV8j9NPPz22zYUXXuj7+GmnneZ8XEm/xFBfX8/555/PXXfdxU9+8pPYeE1NDXfffTdLly7llFNOAeDee+/l8MMPZ+3atRx77LFOX2dIWi156fFpYaux4NJWALuxENRWCBp3Yd1bAaDZs1O+VmOhJmo3EKqjAQl/o7HQENBWaIja681/aNx3YXebfb+Eqha7xeByfwWntgKYjQWXtgLYjQWXtgLYjYVEtBXAbiy4tBXAbiz0dlsB7MZCItoKYDcWXNoKYDcWXNoKHdv732cubQWwGwuJaCuA3VhwaSuA3Viw7pHTlN47LYZkXGJ48MEHqaioYNGiRUyePJn58+czZcoUNm7cyJAhQ3zbP/LII7S2fvTv1J49exg7dixnn3123HannXYa9957b+zvWVkBv/u6kPQzCLNnz+b000+nvLw8bnzdunW0tbXFjR922GGMGDGCNWvWBO6vpaWF2trauIeIiEhfdOutt3LppZcya9YsjjjiCBYtWkRubi733HOPuX1xcTGlpaWxx4oVK8jNzfVNELKysuK2O+CAA5yPLakThGXLlvHGG28wb94838cqKyvJzMykqKgobrykpITKysrAfc6bN4/CwsLYY/jw4Yk+bBERSUUJvMTQ+T+qLcZZudbWVtatWxf3H+FwOEx5eXmX/xH+uLvvvpvzzjuPvLz4M6mrVq1iyJAhjB49mssvv5w9e/Z0+9sQOxbnz0iQbdu2ceWVV7JkyRKys+3Td5/G3LlzqampiT22bduWsH2LiEgKS+AEYfjw4XH/WbX+I7x7924ikQglJSVx45/0H+G9Xn31VTZs2MAll1wSN37aaadx//33s3LlSn72s5+xevVqpk6dSsTxsl/SMgjr1q1j586dHHPMMbGxSCTCiy++yK9+9SueeeYZWltbqa6ujjuLUFVVRWlpaeB+s7KyPtW1FhERkUTZtm0bBQUf5an2xb9Ld999N0cffTSTJk2KGz/vvPNifz766KMZM2YMn/vc51i1ahVf/OIXu73/pE0QvvjFL/Lmm2/Gjc2aNYvDDjuM733vewwfPpyMjAxWrlzJ9OnTAdi4cSNbt26lrKzM+esVpzWSnxZ/wsQKJLqEEQGyQka4LQFLJ4MdSGwxlk4GqIvaM0MrkLjTCAwCVEdzzXFr+eSgkGJ9xD4b9GG7P0i4u8U+jl1B4w7LJ7uEEcEOJLqEEcEOJLqEETu294exEhFGhIBAokMYEexAYm+HEcEOJCYijAh2INEljNgx7g8kuoQRO8atse6HEcEOJCYijAh2INEljAh2INFaAr8hLQHvsW4I/fvR030AFBQUxE0QLIMGDSItLc1X3/+k/wgDNDQ0sGzZMm688cZPPKbPfvazDBo0iM2bN+8fE4QBAwZw1FFHxY3l5eUxcODA2PjFF19MRUUFxcXFFBQUcMUVV1BWVubcYBAREflEvbySYmZmJhMmTGDlypWceeaZAESjUVauXMmcOXO6/NyHH36YlpYWvv71r3/i13n//ffZs2cPQ4cO7f7B0Qdqjl257bbbCIfDTJ8+nZaWFqZMmcKdd96Z7MMSEZEUlIyaY0VFBTNnzmTixIlMmjSJ+fPn09DQwKxZswC44IIL+MxnPuPLMNx9992ceeaZDBw4MG68vr6eG264genTp1NaWsrbb7/NNddcw8EHH8yUKVOcjq1PTRBWrVoV9/fs7GwWLFjAggULknNAIiIi+9C5557Lrl27uPbaa6msrGTcuHEsX748FlzcunUr4XD85amNGzfy0ksv8eyzz/r2l5aWxt/+9jfuu+8+qqurGTZsGKeeeio//vGPnXMQfWqCICIikjRJulnTnDlzAi8pdP6PM8Do0aPxPPsL5eTk8Mwzz7gfhEETBBERkb1S6GZLPdVvJghF4VbyO52msRoLLm0FcGssRDw7ud6OndC1Wgx1UXvJ0ZqovYzrHqOZsCtiJ2s/bLfbA43G8snWGEBNu50Yr2nzj+9qtr/enia7TeG0fLJDWwHsxoJLW6Fje/+4S1sB7MZCQtoKYDYWXNoKsO8aCy5tBbAbC4loK4DdWHBpK4DdWHBpK3SM+99PLm0FsBsLiWgrgN1YcGkrgN1YsJbAz0qzf45k3+o3EwQREZGu6HbP8TRBEBERgaRlEPqqpN+sSURERPoenUEQERFBlxg60wRBREQEdImhk34zQSgIhxjQqbVgNRaS0VZo9uzEeKNxf4XqqH0ce6J2gntXu7+xUNVWaG5r3S8BoLbdn/Zuitgthrp2OzFe3eI/vupm+5hrG+10ucv9FVzaCmA3FlzaCmA3FlzaCmA3FhLRVgC7sdDbbQWwGwsubQWwGwuJaCuA3VhwaSt0bO9/n7m0FcBuLLi0FcBuLCSirQB2Y8GlrQB2Y8G6R056WC2GZOg3EwQREZGu6BJDPE0QREREQJcYOtEEQUREBDRB6EQ1RxEREfHpN2cQckIZ5HYKTVmBRJcwYpCgMGJj1A4SNQSEGuuMQOKHUTtwZYURASrb/YHEHa0BIcVWO6TYYAQSG9vtgFdDmz1e2+w/7vomO1TW0mDvw2X5ZJcwItiBRJcwYsf2/sCfSxgR7EBiIsKIYAcSezuMCHYg0SWMCHYgMRFhxI5xfyDRJYzYMW6NdT+MCHYg0SWMCHYgMRFhRLADiS5hRLADidYS+GFjbF9QBiFev5kgiIiIdEmXGOLoEoOIiIj46AyCiIgIEPI8Ql7PTgH09PP7Ek0QREREQJcYOtElBhEREfHpN2cQskLpZHVKWyeisdDi+VPFQUsnB7UVaqL2Mq7VUX9Se0/EXq/VaisA7Ggt8o8129tWt9qJcauZ0Nxuv3WaW/3LVwM0Nvv30dZkb+vSVgC7seDSVujY3j/u0lYAu7Hg0laAgMZCAtoKsO8aCy5tBbAbCy5tBbAbC4loK4DdWHBpK3SM+99PLm0FsBsLLm0FsBsLiWgrgN1YcGkrgN1YsJbAbw/1zlLLajHE6zcTBBERkS7pEkMcXWIQERERH51BEBERQZcYOtMEQUREBHSJoRNNEERERNAZhM76zQQhLRTuVmshEtA0CLq/gtVYqIva2wa1FfZE7QT3h0ZjYWfAPRestgLAjmb/9rua7Oh1TYu9ln2T0Uxoaw9IgLfab6lIs7H+fqO9j7T67rcVwG4suLQVwG4suLQVwG4sOLUVwGws9JW2AtiNBZe2AtiNBZe2AtiNhUS0FTq297/PXNoKYDcWXNoKYDcWXNoKYDcWEtFWALux4NJWALuxYN0jJ6uXWgwSr99MEERERLqkSwxxNEEQERH5t1S6RNBTqjmKiIiIj84giIiIAHhex6On+0gRmiCIiIigFkNn/XqCYDUWXNoKYDcWXNsKuwKaCdZ9F4LaClUt9j52Ng3wje1uzDW3bWy2k+Rtbf7nE221n6PXGpB+b/aPpzfa26Y7tBUAMur9P5EubYWO7f2vr0tbAezGgktbAezGQl9pK4DdWHBpK4DdWHBpK4DdWEhEW6Fj3BrrflsB7MaCS1sB7MaCS1sB7MZCItoKYDcWXNoKYDcWrLZZmr1b2cf69QRBREQkRi2GOJogiIiIAKFox6On+0gVajGIiIiIj84giIiIgC4xdNKvJwhWINEljAh2INE1jLir3R8kBKhqK/SNbW/xjwHsaAoIOjbm+cbqGuxAWHuzHSTyWvwnmkJt9smncKudJkpr9o+nNwSEEe3VZM0wYtB4Zn33w4hgBxJdwogQEEh0CCOCHUhMSBgRzECiSxgR7ECiSxgR7ECiSxixY9z/c5eIMGLHuP/95BJGBDuQ6BJGBDuQ6BJGBDuQmIgwItiBRJcwItiBxGRKVothwYIF3HLLLVRWVjJ27Fh++ctfMmnSJHPbxYsXM2vWrLixrKwsmps/en95nsd1113HXXfdRXV1NccddxwLFy7kkEMOcTquvvXqiIiIJMvedRB6+nDw4IMPUlFRwXXXXccbb7zB2LFjmTJlCjt37gz8nIKCAnbs2BF7vPfee3Ef//nPf84dd9zBokWLeOWVV8jLy2PKlClxk4ju0ARBREQkSW699VYuvfRSZs2axRFHHMGiRYvIzc3lnnvuCfycUChEaWlp7FFSUhL7mOd5zJ8/nx/+8IecccYZjBkzhvvvv5/t27fz2GOPOR2bJggiIiJ8dImhpw+A2trauEeLcbmxtbWVdevWUV5eHhsLh8OUl5ezZs2awOOsr69n5MiRDB8+nDPOOIO///3vsY9t2bKFysrKuH0WFhYyefLkLvdp0QRBREQEPgop9vQBDB8+nMLCwthj3rx5vi+3e/duIpFI3BkAgJKSEiorK81DHD16NPfccw9/+tOf+P3vf080GuULX/gC77//PkDs81z2GaRfhxRFRET2hW3btlFQ8FF4PCvLDu26Kisro6ysLPb3L3zhCxx++OH8+te/5sc//nFCvsZe/WaC0OK10eLFnzCxGgsubQWA6qj/RQ9qK1S22w2EnW0B2xvLJ1cFtRWa7OWTa43GQluDnRgPNdnPMa3Fn1YOtQW0FezAM2lNRouhyd42sK0QtHyy0VjIqO9+WwHsxoJTWwHMxoJLWwF6f/lkl7YC2I0Fl7YC2I0Fl7YC2I2FRLQVwG4suLQVwG4suLQVwG4suLQVwG4sJKKtAHZjYX9pKwRJZIuhoKAgboJgGTRoEGlpaVRVVcWNV1VVUVpa2q2vl5GRwfjx49m8eTNA7POqqqoYOnRo3D7HjRvXzWfRYf941URERPa1Xm4xZGZmMmHCBFauXBkbi0ajrFy5Mu4sQVcikQhvvvlmbDIwatQoSktL4/ZZW1vLK6+80u197tVvziCIiIj0NRUVFcycOZOJEycyadIk5s+fT0NDQ2ytgwsuuIDPfOYzsQzDjTfeyLHHHsvBBx9MdXU1t9xyC++99x6XXHIJ0NFwuOqqq/jJT37CIYccwqhRo/jRj37EsGHDOPPMM52OTRMEERERkrNQ0rnnnsuuXbu49tprqaysZNy4cSxfvjwWMty6dSvh8Ecn+//1r39x6aWXUllZyQEHHMCECRN4+eWXOeKII2LbXHPNNTQ0NHDZZZdRXV3N8ccfz/Lly8nOti/7BUnqJYbrr7+eUCgU9zjssMNiH29ubmb27NkMHDiQ/Px8pk+f7rtWIyIikhAJbDG4mDNnDu+99x4tLS288sorTJ48OfaxVatWsXjx4tjfb7vttti2lZWVPPnkk4wfPz5uf6FQiBtvvJHKykqam5t57rnnOPTQQ52PK+kZhCOPPDJuRaiXXnop9rGrr76aJ554gocffpjVq1ezfft2zjrrrCQerYiISP+Q9EsM6enpZlqzpqaGu+++m6VLl3LKKacAcO+993L44Yezdu1ajj32WKev0+y1k9mpxWA1FlzaCgC7Iv6UalBbYUdrkTleZbQVAHY1++PXu5r891YAqK63Wwyt9f7EeKjeftnTjaYBgBVuDrznQlCLwQj+pzf2vK0AdmPBpa0AdmPBpa0AdmOhr7QVwG4suLQVwG4suLQVwG4suLQVOrb3v/8S0VYAu7Hg0lYAu7Hg0lYAu7Hg0lYAu7GQiLYC2I2FRLQVWjz/97rF6517KCfrXgx9VdLPIGzatIlhw4bx2c9+lvPPP5+tW7cCsG7dOtra2uJWgzrssMMYMWJEl6tBtbS0+FawEhER+URRLzGPFJHUCcLkyZNZvHgxy5cvZ+HChWzZsoUTTjiBuro6KisryczMpKioKO5zPmk1qHnz5sWtXjV8+PB9/CxERCQlJCmD0Fcl9RLD1KlTY38eM2YMkydPZuTIkTz00EPkBJzi/CRz586loqIi9vfa2lpNEkRERBwl/RLDxxUVFXHooYeyefNmSktLaW1tpbq6Om6bT1phKisrK7aCVXdWshIREQEIkYCbNSX7SSRQ0kOKH1dfX8/bb7/NN77xDSZMmEBGRgYrV65k+vTpAGzcuJGtW7c6rwYFUB+NEOp0bcgKJO6J2mcuPozYCajKNn8gMWjp5KAwYlXzAHPcWj65usE+vpag5ZMb/M8xoz4gYBgQUrSChwE5p8Dx9Gb/ebeMoJBiUBixofvLJ7uEESEgkOgQRgQ7kNhXwohgBxJdwohgBxJdwogd4/73pEsYsWPcGut5GBHsQKJLGBHsQKJLGBHsQKJLGBHsQGIiwoiw7wKJ1hL4zb0UUnRdCTFwHykiqROE7373u0ybNo2RI0eyfft2rrvuOtLS0pgxYwaFhYVcfPHFVFRUUFxcTEFBAVdccQVlZWXODQYRERFxk9QJwvvvv8+MGTPYs2cPgwcP5vjjj2ft2rUMHjwY6FgQIhwOM336dFpaWpgyZQp33nlnMg9ZRERSlGqO8ZI6QVi2bFmXH8/OzmbBggUsWLCgl45IRET6rUS0EFJogtCnQooiIiLSN/SpkKKIiEiyhDyPUA9Dhj39/L6k30wQaqNpRKPxJ0ysxsKu9oBlj9vtpsGOtiLf2PZm/xgEtxV2NdrLJ9c2+hPjTQ12upyg5ZPr/SeJ0gNaDOlN9q6tZkJaq/1DELTUstViSG/seVsB7MaCU1sBzMaCS1sB9l1jIRFtBbAbCy5tBbAbCy5tBbAbCy5thY5x//spEW0FsBsLLm0FsBsLLm0FsBsLLm0FsBsLfaWtAHZjwVoCvz7aSy2G6L8fPd1HitAlBhEREfHpN2cQREREuqJLDPE0QRAREQG1GDrRBEFERAS0kmInyiCIiIiIT785g/BhNIeWTi0Gq7Hg0lYAu7Hg2laoCbi/QnOjPzHuBbQVMoy2AkBGnT/FnFFvbkp6U/ebCcEthoDxZn+0N70xIF3e0P22AtiNBZe2AtiNhd5uK4DdWEhEWwHsxoJLWwHsxoJLW6Fje/970qWtAHZjIRFtBbAbCy5tBbAbCy5tBbAbCy5tBbAbC32lrQB2Y8G6R059tHdugaSVFOP1mwmCiIhIl3SJIY4uMYiIiIiPziCIiIgAoWjHo6f7SBWaIIiIiIAuMXTSbyYIu9sH0NgeH36xAokuYUSwA4nOYcQGOxDmNfhfnvQ6OxAWtHxyhrEya0aD/QbOaOx+8DAcEEYMt9nTZyuQmNYUEEZstINYTssnO4QRwQ4k9nYYEexAYiLCiGAHEl3CiB3j/vefSxixY9wa634YEexAYiLCiGAHEl3CiGAHEl3CiGAHEl3CiGAHEhMRRgQ7kOgSRgQ7kGgtgd/QW0stS5x+M0EQERHpkhZKiqMJgoiICFpquTO1GERERMRHZxBERERAIcVONEEQERGBjvxAT/OQqTM/6D8ThN2RfHLa45+u1VhwaSuA3VgIbCvU2+lyjLYCQLqxfHJGUFshYPnkjHr/uzXTGANIb7R/MtJa/eNBbYVwi51WDrf4082BbYWmgKaBw/LJLm0F2HeNBZe2AtiNhUS0FcBuLLi0FcBuLLi0FTrG/e8/l7YC2I2FRLQVwG4suLQVwG4suLQVwG4suLQVoPeXT3ZpK4DdWLCWwG9sT8DPZzcogxBPGQQRERHx6TdnEERERLrkkYAMQkKOpE/QBEFERAQUUuxElxhERETER2cQREREoKPBYGdu3faRIvrNBKGqrZCstvgEsNVYcGkrgN1YCGwr1AckjY22AkBGg/+dmm4Hoc22Qsc+/OMZ9QHJ8KagBoJ/PKjFEGq1WwKhZn8SOhTUNAgYJ6DF4LX6U+C93VYAu7Hg0lYAu7GQiLYC2I0Fl7ZCx/b+96RLWwHsxoJLWwHsxkIi2gpgNxZc2gpgNxZc2gpgNxb6SlsB7MaCS1sB7MaCdY+cpoh9DImWrBbDggULuOWWW6isrGTs2LH88pe/ZNKkSea2d911F/fffz8bNmwAYMKECdx0001x21944YXcd999cZ83ZcoUli9f7nRcusQgIiKSJA8++CAVFRVcd911vPHGG4wdO5YpU6awc+dOc/tVq1YxY8YMXnjhBdasWcPw4cM59dRT+eCDD+K2O+2009ixY0fs8cADDzgfmyYIIiIi8FFIsacPB7feeiuXXnops2bN4ogjjmDRokXk5uZyzz33mNsvWbKEb33rW4wbN47DDjuM3/72t0SjUVauXBm3XVZWFqWlpbHHAQcc4Pzt0ARBREQEen2C0Nrayrp16ygvL4+NhcNhysvLWbNmTbf20djYSFtbG8XFxXHjq1atYsiQIYwePZrLL7+cPXv2dPu49uo3GQQREZHeUltbG/f3rKwssrLi80S7d+8mEolQUlISN15SUsL//u//duvrfO9732PYsGFxk4zTTjuNs846i1GjRvH222/zgx/8gKlTp7JmzRrS0uyciKXfTBB2NBeSmR4f3rICiS5hRIDmBiMQFrB0cmAYMWD55HQjc5VZZ89Og5ZPzjQCiRkNduAnrTEgYGgED0NtAWG/VjvoFGrxjweGEY2lk8EOI4IdSOztMCLYgUSXMCLYgcREhBE7xv2/GFzCiB3j1lj3w4hgBxJdwohgBxITEUYEO5DoEkYEO5DoEkYEO5DYV8KIYAcSXcKIHeP+38HWEvgtbfYxJ1wC10EYPnx43PB1113H9ddf37N9d3LzzTezbNkyVq1aRXb2R78nzjvvvNifjz76aMaMGcPnPvc5Vq1axRe/+MVu77/fTBBERES6lMCa47Zt2ygo+Ghi1PnsAcCgQYNIS0ujqqoqbryqqorS0tIuv8z//M//cPPNN/Pcc88xZsyYLrf97Gc/y6BBg9i8ebPTBEEZBBERET6qOfb0AVBQUBD3sCYImZmZTJgwIS5guDdwWFZWFnicP//5z/nxj3/M8uXLmThx4ic+r/fff589e/YwdOhQp++HJggiIiJJUlFRwV133cV9993HP/7xDy6//HIaGhqYNWsWABdccAFz586Nbf+zn/2MH/3oR9xzzz0cdNBBVFZWUllZSX19xy196+vr+a//+i/Wrl3Lu+++y8qVKznjjDM4+OCDmTJlitOx6RKDiIgIJOVeDOeeey67du3i2muvpbKyknHjxrF8+fJYcHHr1q2Ewx/9X37hwoW0trby1a9+NW4/ezMOaWlp/O1vf+O+++6jurqaYcOGceqpp/LjH//YPIvRFU0QREREAKIehHo4QYi6f/6cOXOYM2eO+bFVq1bF/f3dd9/tcl85OTk888wzzsdg6TcThJ0tA8jo1GKwGgtObQWAen8COa3Ora2Q4bB8cuCSyvX20scZ9f5kclq93QYINwUkhdusFkPA0qcB457VTAhoPLi0FWDfNRZc2gpgNxZc2gpgNxYS0VYAu7Hg0lboGPe//1zaCmA3FlzaCmA3FhLRVgC7seDSVgC7seDSVoB911hIRFsB7MaCS1sB7MaCtQR+a7P9fZZ9q99MEERERLqk2z3H0QRBREQEgARMEEidCYJaDCIiIuKjMwgiIiKgSwydaIIgIiIC/24g9H6Loa/qNxOEPY25pIXik+NWY8GlrQB2YyERbYWg8cygtkLQ/RUa/CnmcKOdCA4FJYWNtkFQoyCwmWBtH7C+em+3FcBuLLi0FcBuLLi0FcBuLCSirdCxvf996dJWALux4NJWALux4NJWALuxkIi2AtiNBZe2AtiNhd5uK4DdWEhEWwHsxoJLWwHsxoJ1j5y2FrUYkqHfTBBERES65EU7Hj3dR4rQBEFERASUQehEEwQRERFQBqGTpNccP/jgA77+9a8zcOBAcnJyOProo3n99ddjH/c8j2uvvZahQ4eSk5NDeXk5mzZtSuIRi4iIpL6kThD+9a9/cdxxx5GRkcHTTz/NW2+9xS9+8QsOOOCA2DY///nPueOOO1i0aBGvvPIKeXl5TJkyheZme2lWERGRT2XvJYaePlJEUi8x/OxnP2P48OHce++9sbFRo0bF/ux5HvPnz+eHP/whZ5xxBgD3338/JSUlPPbYY5x33nnd/lo1jdmkEZ8cNxsLDm0FsBsLgW2FuoC2QoM9bjUWrHsrAKTX2ynmcIP/HgihxoDJVUBS2LPaBkH3XHBpIAS0Enq7rQB2Y8GlrQB2Y8GlrQB2YyERbYWOcWus+20FsBsLLm0FsBsLLm0FsBsLiWgrgN1YcGkrgN1Y6O22AtiNhUS0FTrG/W0Dl7YC2I0F6x45kcZe+qfKIwEZhIQcSZ+Q1DMIjz/+OBMnTuTss89myJAhjB8/nrvuuiv28S1btlBZWUl5eXlsrLCwkMmTJ7NmzRpzny0tLdTW1sY9RERExE1SJwjvvPMOCxcu5JBDDuGZZ57h8ssv59vf/jb33XcfAJWVlQCx+2LvVVJSEvtYZ/PmzaOwsDD2GD58+L59EiIikhp0iSFOUicI0WiUY445hptuuonx48dz2WWXcemll7Jo0aJPvc+5c+dSU1MTe2zbti2BRywiIikrGk3MI0UkdYIwdOhQjjjiiLixww8/nK1btwJQWloKQFVVVdw2VVVVsY91lpWVRUFBQdxDRERE3CQ1pHjcccexcePGuLF//vOfjBw5EugILJaWlrJy5UrGjRsHQG1tLa+88gqXX36509dqacwk7HUKhhmBRJcwItiBxMClkx3CiB378QePXMKIYAcSvWZ7W1rscc8KJAYGDO3nYgYPg1YcCzpFF7JfAyuQ6BJGBDuQ6BJGBDuQ6BJG7Bj3B8gSEUbsGPd/X13CiGAHEl3CiGAHEl3CiGAHEhMRRgQ7kOgSRoR9F0h0CSOCHUhMRBgR7ECiSxgR7ECitQR+pDHgZz/RtFBSnKROEK6++mq+8IUvcNNNN3HOOefw6quv8pvf/Ibf/OY3AIRCIa666ip+8pOfcMghhzBq1Ch+9KMfMWzYMM4888xkHrqIiKQaTRDiJHWC8PnPf55HH32UuXPncuONNzJq1Cjmz5/P+eefH9vmmmuuoaGhgcsuu4zq6mqOP/54li9fTna2XR8TERGRnkv6Ustf/vKX+fKXvxz48VAoxI033siNN97Yi0clIiL9jpZajpP0CYKIiEhf4HlRvB7ejbGnn9+XaIIgIiICHfmBnp4BUAZhP1SfAZH4JLLVWHBpK3Rs738zBLUYXNoKAGkN/hSzS1sBAhoLTQHbJmKZ5KAfLmtWnYC2AtiNBZe2AtiNBZe2AtiNBZe2AtiNhUS0FcBuLLi0FcBuLLi0FcBuLLi0FcBuLCSirQB2Y6G32wpgNxZc2gpgNxYS0VYAu7Hg0lYAu7FgLYEfbUqd/5XvT/rPBEFERKQrXgIyCDqDICIikmKiUQj18GxFCmUQkrqSooiIiPRNOoMgIiICusTQiSYIIiIigBeN4vXwEoNqjvuhtPow4fb4KypWY8GlrRA07tpWcLm/QmBbIaCZYDUWvFY71b1P76NgSUBbAezGgktbAezGgktbAezGgktboWN7//ckEW0FsBsLLm0FsBsLLm0FsBsLLm0FsBsLiWgrgN1Y6O22AtiNBZe2AtiNhUS0FcBuLLi0FcBuLFj3yKHJbm/IvtVvJggiIiJd0iWGOJogiIiIQMciSSFNEPZSi0FERER8dAZBREQE/v2//56ug5A6ZxD6zQQhvSFEWnt8WMkKJLqEEcEOJCYijAh2INEljAh2IDFwSeVELJMcxAokJiCMCHYg0SWMCHYg0SWM2DHuD5C5hBE7xq2xnocRwQ4kuoQRwQ4kuoQRwQ4kuoQRwQ4kJiKMCPsukOgSRgQ7kOgSRuwY9wcJExFGBDuQ6BRGBDOQaC2BH2runZPdXtTD6+ElBk8TBBERkRTjRen5GYTUqTkqgyAiIiI+OoMgIiKCLjF0pgmCiIgI6BJDJyk/Qdg7m4u2+MNYYSMbGG61Z3/htoCgWJv/zRAKCAESCQgpRuwV4EJR/wF6UXtbvIDVEY1xz7NDUYEz30SEFLECZAEhxYB9hwJ+8EJR/75DUXvfAU8dL+LfR6Td/nrt7faPTXubP1QWabWPI9JqB+oixnsymh4QRgwaD9tPMhr2v/8iITsg24493ma8/1raA1YIbLN/DprS/eMNafYxZ6XZr0F62D8eDggptgcsnZsVMJ4WsMCnixbjvdoc8P6tjwaN+w+kIWDbxnb7+9cU8X+vW9rs16u12f4d0tZij0ca/T8HkUb7mxdtCvhH01gh0QokRps7fn/v6/+dt9PW43WS2rG/v/ujkJdK50MM77//PsOHD0/2YYiISA9t27aNAw88MOH7bW5uZtSoUVRWViZkf6WlpWzZsoXsbHu59v1Fyk8QotEo27dvZ8CAAdTV1TF8+HC2bdtGQYFdC9rf1dbWpvRzTPXnB3qOqULPMXE8z6Ouro5hw4YRDu+bbH1zczOtAfepcZWZmbnfTw6gH1xiCIfDsRln6N9d/IKCgpT9gd0r1Z9jqj8/0HNMFXqOiVFYWLhP95+dnZ0S/6gnkmqOIiIi4qMJgoiIiPj0qwlCVlYW1113HVlZ9nK7qSDVn2OqPz/Qc0wVeo6yv0v5kKKIiIi461dnEERERKR7NEEQERERH00QRERExEcTBBEREfHpNxOEBQsWcNBBB5Gdnc3kyZN59dVXk31In9qLL77ItGnTGDZsGKFQiMceeyzu457nce211zJ06FBycnIoLy9n06ZNyTnYT2nevHl8/vOfZ8CAAQwZMoQzzzyTjRs3xm3T3NzM7NmzGThwIPn5+UyfPp2qqqokHbG7hQsXMmbMmNgiM2VlZTz99NOxj+/vz6+zm2++mVAoxFVXXRUb29+f4/XXX08oFIp7HHbYYbGP7+/Pb68PPviAr3/96wwcOJCcnByOPvpoXn/99djHU+F3jvj1iwnCgw8+SEVFBddddx1vvPEGY8eOZcqUKezcuTPZh/apNDQ0MHbsWBYsWGB+/Oc//zl33HEHixYt4pVXXiEvL48pU6bQ3Oy/YVVftXr1ambPns3atWtZsWIFbW1tnHrqqTQ0NMS2ufrqq3niiSd4+OGHWb16Ndu3b+ess85K4lG7OfDAA7n55ptZt24dr7/+OqeccgpnnHEGf//734H9//l93Guvvcavf/1rxowZEzeeCs/xyCOPZMeOHbHHSy+9FPtYKjy/f/3rXxx33HFkZGTw9NNP89Zbb/GLX/yCAw44ILZNKvzOEYPXD0yaNMmbPXt27O+RSMQbNmyYN2/evCQeVWIA3qOPPhr7ezQa9UpLS71bbrklNlZdXe1lZWV5DzzwQBKOMDF27tzpAd7q1as9z+t4ThkZGd7DDz8c2+Yf//iHB3hr1qxJ1mH22AEHHOD99re/TannV1dX5x1yyCHeihUrvBNPPNG78sorPc9Ljdfwuuuu88aOHWt+LBWen+d53ve+9z3v+OOPD/x4qv7OEc9L+TMIra2trFu3jvLy8thYOBymvLycNWvWJPHI9o0tW7ZQWVkZ93wLCwuZPHnyfv18a2pqACguLgZg3bp1tLW1xT3Pww47jBEjRuyXzzMSibBs2TIaGhooKytLqec3e/ZsTj/99LjnAqnzGm7atIlhw4bx2c9+lvPPP5+tW7cCqfP8Hn/8cSZOnMjZZ5/NkCFDGD9+PHfddVfs46n6O0f6wSWG3bt3E4lEKCkpiRsvKSlJ2K09+5K9zymVnm80GuWqq67iuOOO46ijjgI6nmdmZiZFRUVx2+5vz/PNN98kPz+frKwsvvnNb/Loo49yxBFHpMzzW7ZsGW+88Qbz5s3zfSwVnuPkyZNZvHgxy5cvZ+HChWzZsoUTTjiBurq6lHh+AO+88w4LFy7kkEMO4ZlnnuHyyy/n29/+Nvfddx+Qmr9zpEPK381R9n+zZ89mw4YNcdd2U8Xo0aNZv349NTU1/OEPf2DmzJmsXr062YeVENu2bePKK69kxYoVKXuXvKlTp8b+PGbMGCZPnszIkSN56KGHyMnJSeKRJU40GmXixIncdNNNAIwfP54NGzawaNEiZs6cmeSjk30p5c8gDBo0iLS0NF9yuKqqitLS0iQd1b6z9zmlyvOdM2cOf/7zn3nhhRdit+2GjufZ2tpKdXV13Pb72/PMzMzk4IMPZsKECcybN4+xY8dy++23p8TzW7duHTt37uSYY44hPT2d9PR0Vq9ezR133EF6ejolJSX7/XPsrKioiEMPPZTNmzenxGsIMHToUI444oi4scMPPzx2KSXVfufIR1J+gpCZmcmECRNYuXJlbCwajbJy5UrKysqSeGT7xqhRoygtLY17vrW1tbzyyiv71fP1PI85c+bw6KOP8vzzzzNq1Ki4j0+YMIGMjIy457lx40a2bt26Xz3PzqLRKC0tLSnx/L74xS/y5ptvsn79+thj4sSJnH/++bE/7+/PsbP6+nrefvtthg4dmhKvIcBxxx3nqxj/85//ZOTIkUDq/M4RQ7JTkr1h2bJlXlZWlrd48WLvrbfe8i677DKvqKjIq6ysTPahfSp1dXXeX/7yF+8vf/mLB3i33nqr95e//MV77733PM/zvJtvvtkrKiry/vSnP3l/+9vfvDPOOMMbNWqU19TUlOQj777LL7/cKyws9FatWuXt2LEj9mhsbIxt881vftMbMWKE9/zzz3uvv/66V1ZW5pWVlSXxqN18//vf91avXu1t2bLF+9vf/uZ9//vf90KhkPfss896nrf/Pz/Lx1sMnrf/P8fvfOc73qpVq7wtW7Z4/+///T+vvLzcGzRokLdz507P8/b/5+d5nvfqq6966enp3k9/+lNv06ZN3pIlS7zc3Fzv97//fWybVPidI379YoLgeZ73y1/+0hsxYoSXmZnpTZo0yVu7dm2yD+lTe+GFFzzA95g5c6bneR21ox/96EdeSUmJl5WV5X3xi1/0Nm7cmNyDdmQ9P8C79957Y9s0NTV53/rWt7wDDjjAy83N9b7yla94O3bsSN5BO7rooou8kSNHepmZmd7gwYO9L37xi7HJgeft/8/P0nmCsL8/x3PPPdcbOnSol5mZ6X3mM5/xzj33XG/z5s2xj+/vz2+vJ554wjvqqKO8rKws77DDDvN+85vfxH08FX7niJ9u9ywiIiI+KZ9BEBEREXeaIIiIiIiPJggiIiLiowmCiIiI+GiCICIiIj6aIIiIiIiPJggiIiLiowmCiIiI+GiCICIiIj6aIIiIiIiPJggi+6Fdu3ZRWlrKTTfdFBt7+eWXyczMjLurnojIp6V7MYjsp5566inOPPNMXn75ZUaPHs24ceM444wzuPXWW5N9aCKSAjRBENmPzZ49m+eee46JEyfy5ptv8tprr5GVlZXswxKRFKAJgsh+rKmpiaOOOopt27axbt06jj766GQfkoikCGUQRPZjb7/9Ntu3bycajfLuu+8m+3BEJIXoDILIfqq1tZVJkyYxbtw4Ro8ezfz583nzzTcZMmRIsg9NRFKAJggi+6n/+q//4g9/+AN//etfyc/P58QTT6SwsJA///nPyT40EUkBusQgsh9atWoV8+fP53e/+x0FBQWEw2F+97vf8X//7/9l4cKFyT48EUkBOoMgIiIiPjqDICIiIj6aIIiIiIiPJggiIiLiowmCiIiI+GiCICIiIj6aIIiIiIiPJggiIiLiowmCiIiI+GiCICIiIj6aIIiIiIiPJggiIiLiowmCiIiI+Pz/81IjcTmSjboAAAAASUVORK5CYII=" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "HBox(children=(HTML(value=\".\\\\tmpjzrlzgss.pdf
\"), HTML(val…" + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "8d95a28c505d40e7b9820b2a48f1366a" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 14 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:28.064651Z", + "start_time": "2024-10-30T19:30:28.050633Z" + } + }, + "cell_type": "code", + "source": "solution_filename = [f for f in os.listdir(\"./out\") if \"solution\" in f][0]", + "id": "c58862c96237cccf", + "outputs": [], + "execution_count": 15 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:28.079639Z", + "start_time": "2024-10-30T19:30:28.065644Z" + } + }, + "cell_type": "code", + "source": [ + "%%writefile to_vtk.jl\n", + "using Trixi2Vtk\n", + "trixi2vtk(joinpath(\"out\", ARGS[1]))" + ], + "id": "8d35a730ab764f86", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Overwriting to_vtk.jl\n" + ] + } + ], + "execution_count": 16 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.084503Z", + "start_time": "2024-10-30T19:30:28.080639Z" + } + }, + "cell_type": "code", + "source": "subprocess.run([\"julia\", \"to_vtk.jl\", solution_filename], check=True)", + "id": "960cee65b6c01544", + "outputs": [ + { + "data": { + "text/plain": [ + "CompletedProcess(args=['julia', 'to_vtk.jl', 'solution_000030.h5'], returncode=0)" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "execution_count": 17 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.100519Z", + "start_time": "2024-10-30T19:30:40.085510Z" + } + }, + "cell_type": "code", + "source": [ + "try:\n", + " vtu_filename = [f for f in os.listdir(\"./\") if \"vtu\" in f and \"celldata\" not in f][0]\n", + " mesh = meshio.read(vtu_filename)\n", + " trixi_points = ((mesh.points[:,:2] + 1)*SETUP[\"nx\"]*SETUP[\"polydeg\"]/2).round().astype(np.int16)\n", + " assert trixi_points.shape[0] == SETUP[\"nx\"]**2 * (SETUP[\"polydeg\"] + 1)**2\n", + "except Exception as e:\n", + " e.args += (list(os.walk(os.path.curdir)),)\n", + " raise e" + ], + "id": "451911db51e18682", + "outputs": [], + "execution_count": 18 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.116510Z", + "start_time": "2024-10-30T19:30:40.101521Z" + } + }, + "cell_type": "code", + "source": [ + "try:\n", + " trixi_output = np.zeros_like(pympdata_result_state)\n", + " for i in range(trixi_points.shape[0]):\n", + " trixi_output[trixi_points[i][0], trixi_points[i][1]] = mesh.point_data['scalar'][i][0]\n", + "except Exception as e:\n", + " e.args += (list(mesh.point_data.keys()),)\n", + " e.args += (list(mesh.points.shape),)\n", + " raise e" + ], + "id": "58595cff705f196c", + "outputs": [], + "execution_count": 19 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.684640Z", + "start_time": "2024-10-30T19:30:40.117511Z" + } + }, + "cell_type": "code", + "source": [ + "plt.imshow(trixi_output, cmap='viridis', vmin=vmin, vmax=vmax)\n", + "plt.colorbar()\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')\n", + "plt.title(\"Trixi solution\")\n", + "show_plot(inline_format='png')" + ], + "id": "a126dbabdf719ee3", + "outputs": [ + { + "data": { + "text/plain": [ + "
" + ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAHHCAYAAADaqqCfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABi/klEQVR4nO3de3xV1Z3//9dJyIUACXJLSA0XRwVRuQgS46WipkbGcUSpVUcr4u0rghUyrZaZCnipIFZFLYJaBXSqeGmltSqKUXAcATXKVKtDRVEikKD+5JIAScjZvz9OPXqyPydkkU1OOHk/H4/zeMjKzsk+ubmy9/u9VsjzPA8RERGR70lJ9AmIiIhI26MJgoiIiPhogiAiIiI+miCIiIiIjyYIIiIi4qMJgoiIiPhogiAiIiI+miCIiIiIjyYIIiIi4qMJgoijUaNGMWrUKOf3W758OaFQiOXLlwd+TvHMmDGDUCgU+PP269ePSy+9NPDnFZG2QxMEaXdCoVCzHq35P/K26M0332TGjBls3bo10aciIgnQIdEnINLaHnvssZh/P/rooyxbtsw3fsQRR5jv//LLL+/Tx/3hD3/Irl27SE9P36f3b21vvvkmN910E5deeildu3aNedvatWtJSdHfFyLJTBMEaXcuvvjimH+vWrWKZcuW+cYb27lzJ1lZWfv8P/iUlBQyMzP36X3bmoyMjESfgojsZ/oTQMQwatQojjrqKMrLy/nhD39IVlYW//Ef/xF92/czCOPGjSMzM5OPPvoo5jlKSko46KCD2LRpE9D8DMKOHTuYPHky/fr1IyMjg169evGjH/2Id999N+a4p59+muHDh9OxY0d69OjBxRdfzMaNG5t87s8++4xQKMTChQt9bwuFQsyYMQOIZBd+8YtfANC/f//obZfPPvsMsDMIn376Keeddx7dunUjKyuL4447jueffz7mmG8/B0899RS//vWvOfjgg8nMzOS0005j3bp1TZ67iLQuXUEQiePrr79m9OjRXHDBBVx88cXk5uaax91zzz28+uqrjBs3jpUrV5KamsoDDzzAyy+/zGOPPUZ+fr7Tx7366qt55plnmDRpEoMGDeLrr7/mjTfe4KOPPuKYY44BYOHChYwfP55jjz2WmTNnUlVVxT333MP//M//8N577/luCbg699xz+fvf/84TTzzB3XffTY8ePQDo2bOneXxVVRXHH388O3fu5Gc/+xndu3dn0aJF/Ou//ivPPPMM55xzTszxs2bNIiUlhZ///Ods27aN2bNnc9FFF7F69eoWnbeIBMgTaecmTpzoNf5ROPnkkz3Amz9/vu/4k08+2Tv55JNjxl566SUP8G699Vbv008/9Tp37uyNGTMm5pjXXnvNA7zXXnutyfPJycnxJk6cGPftdXV1Xq9evbyjjjrK27VrV3T8L3/5iwd406ZNi45Nnz495rWtX7/eA7wFCxb4nhfwpk+fHv33HXfc4QHe+vXrfcf27dvXGzduXPTfkydP9gDvv//7v6NjO3bs8Pr37+/169fPa2ho8Dzvu8/BEUcc4dXW1kaPveeeezzAe//99+O+bhFpXbrFIBJHRkYG48ePb9axp59+Ov/v//0/br75Zs4991wyMzN54IEH9unjdu3aldWrV0dvTTT2zjvvsGXLFq655pqYTMOZZ57JwIEDfZf1W8MLL7zAyJEjOfHEE6NjnTt35qqrruKzzz7jww8/jDl+/PjxMVmOk046CYjcphCRtkETBJE4fvCDHzgFEn/zm9/QrVs31qxZw7333kuvXr326ePOnj2bDz74gIKCAkaOHMmMGTNi/sf5+eefAzBgwADf+w4cODD69tb0+eefm+fzbROk8Tn16dMn5t8HHXQQAN98881+OkMRcaUJgkgcHTt2dDr+vffeY8uWLQC8//77+/xxf/KTn/Dpp59y3333kZ+fzx133MGRRx7Jiy++uM/P+a14iyY1NDS0+LldpKammuOe57XqeYhIfJogiASgpqaG8ePHM2jQIK666ipmz57N22+/vc/P17t3b6655hqWLFnC+vXr6d69O7/+9a8B6Nu3LxBZi6CxtWvXRt9u+fYv9caLH1lXHVxWYOzbt695Pv/3f/8Xc84icuDQBEEkADfccAMbNmxg0aJF3HXXXfTr149x48ZRW1vr9DwNDQ1s27YtZqxXr17k5+dHn2vEiBH06tWL+fPnxzz/iy++yEcffcSZZ54Z9/mzs7Pp0aMHr7/+esz4/fff7zu2U6dOgH8yYfnnf/5n3nrrLVauXBkdq6mp4cEHH6Rfv34MGjRor88hIm2Lao4iLfTqq69y//33M3369GgNccGCBYwaNYobb7yR2bNnN/u5duzYwcEHH8yPf/xjhgwZQufOnXnllVd4++23ufPOOwFIS0vj9ttvZ/z48Zx88slceOGF0Zpjv379mDJlSpMf44orrmDWrFlcccUVjBgxgtdff52///3vvuOGDx8OwH/+539ywQUXkJaWxllnnRWdOHzfL3/5S5544glGjx7Nz372M7p168aiRYtYv349f/jDH7TqosgBSBMEkRbYsWMHl112GcOGDeM///M/o+MnnXQS1113HXfeeSfnnnsuxx13XLOeLysri2uuuYaXX36ZP/7xj4TDYQ499FDuv/9+JkyYED3u0ksvJSsri1mzZnHDDTfQqVMnzjnnHG6//fa9roEwbdo0vvzyS5555hmeeuopRo8ezYsvvugLVR577LHccsstzJ8/n6VLlxIOh1m/fr05QcjNzeXNN9/khhtu4L777mP37t0MHjyY5557rskrGiLSdoU8pYJERESkEV33ExERER9NEERERMRHEwQRERHx0QRBREREfDRBEBERER9NEERERMTngFgHYe7cudxxxx1UVlYyZMgQ7rvvPkaOHNms9w2Hw2zatIkuXbo4LR0rIiJtg+d57Nixg/z8/P226Nbu3bupq6sL5LnS09Njdlo9YCV2t+m9W7x4sZeenu498sgj3t/+9jfvyiuv9Lp27epVVVU16/0rKio8QA899NBDjwP8UVFRsV/+P7Nr1y4vr1dqYOeZl5fn7dq1a7+ca2tq8wslFRYWcuyxx/Lb3/4WiFwRKCgo4Nprr+WXv/zlXt9/27ZtdO3alYNn/IqURjO6DtX+KwoddtpXGdKq7edPq/F/+tJqwnGewx7vUFNvjqfu3OMbS9lpr+0f2rXbHPd2+2fEXm2cY+vs8/Dq/edBOIDd/1LsHf1CafaFrVB6mj2e4Z+phzLtbZq9jvasPpyV4RtryLLPY08n+zzqO/v/sqnPsv/aqe9sf5/V+xcpZE8n+0c03ni4k/21Se3k//p27mR/P+V03GWO98ry/yDkZuwwj81N32aO90rb7hvrkWo/R49U+wcvJ8X/fZ0T5w/LjiH7+ywjZH8d0+Ic76Le838Naj3752uXcSzANuPXxbaw/X39VUPnOONdfGNb6rPNY6vqcuzxWv9zAHy1y//N+s3OLPPYmhr/zxfAnhr/1yClxv/5D+/ezRc33crWrVvJybHPsyW2b99OTk4On5f3I7tLy65QbN8Rpu/wz9i2bRvZ2fbn+kDRpm8x1NXVUV5eztSpU6NjKSkpFBcXx2wK8321tbUxG9js2BH5xZOSmembIKTs8f+STm2wf3GnxrnylFrv/yXdoS7ORCAtzngH+xdSaqr/F0qc/6cSSrH/Z+EZ3+teyD4PL84dGM+6NRMK4DJfnF/EoVCcCULI/uUYSvGPh1LsX0heqj0eNsZDHez/gZBmj3tp/s9JON3+PIXT7U922Di9cKb9tU2JM05H+384KVn+z3eq/fucDllxJrlZ/h+E9Ez785EZZ0LX0ZgAZsX5GeiUan/+OhuXmeP9Xs+K872aEWc8LYDv7Xrj765a64cR6BDnb7QG40uwJ2w/x84G+/PXcY//c51Zb39d0uviTJg62D93HYyfsVTsn68Uz56Yp4SNCUKc1wJuO4zui85dQnTu0rKPESZ5bmW36ZDiV199RUNDA7m5uTHjubm5VFZWmu8zc+ZMcnJyoo+CgoLWOFURETnANXjhQB7Jok1PEPbF1KlT2bZtW/RRUVGR6FMSEZEDQBgvkEeyaNO3GHr06EFqaipVVVUx41VVVeTl5Znvk5GRQUaGfZlLREREmqdNTxDS09MZPnw4ZWVljBkzBoiEFMvKypg0aZLTc3WoDvkyB2k1/ntFaXZWygwjAqQbwcO4YcTqOGHEOCFFK5AY2hkvjGiHzbzd/uO9OFUeM4wI+y2QGD+MGCdrEKc2FMr0Twi9rOaHEQEajODhns7NDyMC1HUyQorxwoh27ot6I3gYN4zYOU4YsbP9/dSlk/97oWuWHUbMzbJ/EPIy/QHD3nHCiLlp9njPDv7n6BknpNjVCCNGxv1jWW0kjAh2IHFnnGO3xrkivdUIJH5phA4BvtwTJ3hY7w/0bY4TRqzcbT/Hl7vsAOT/V+MPsOyosX/uGqrtr4EVSOxg/F4O726d+/phwrT0BkHLn6HtaNMTBIDS0lLGjRvHiBEjGDlyJHPmzKGmpobx48cn+tRERCSJNHgeDS0s9rX0/duSNj9BOP/88/nyyy+ZNm0alZWVDB06lKVLl/qCiyIiIhKcNj9BAJg0aZLzLQUREREXQYQMFVIUERFJMmE8GjRBiGo3E4QOO0O+RZCsQKJLGBHsQGIQYUSwA4kuYUSwA4mtHUYEO5AYRBgR7ECiSxgR7ECiSxgxcrw/SFVv57vMMCLYgcQgwohgBxJdwohgBxJdwohgBxJdwohgBxLbShgR7ECiSxgR7ECiSxgR7EBiEGFEgGojkOgSRgR7NVsrPN5QmzyLDx1I2s0EQUREpCm6xRBLEwQRERHUYmgs6VZSFBERkZbTFQQREREg/I9HS58jWWiCICIiAjQE0GJo6fu3Je1mgpBW7d+y2WosuLQVwG4sBNFWALux4NJWgDiNhVZuK4DdWAiirQB2Y8GlrQB2Y8GlrQD28skubQWwGwtBtBXAbiy4tBXAbiy4tBXAbiy4tBXAbiwE0VYAu7Hg0lYAu7Hg0lYAu7Hg0lYAu7Hg2lZwWT7Zpa0AzV/uPk7JJXANXuTR0udIFsogiIiIiE+7uYIgIiLSFGUQYmmCICIiAoQJ0UDLFmUKt/D92xLdYhAREREfXUEQEREBwl7k0dLnSBbtZoKQVuORWh/7lbMaCy5tBbAbC0G0FSLjxl4MLm0F2G+NBZe2AtiNhSDaCmA3FlzaCmA3FuK2FRz2V3BpK4DdWAiirQB2Y8GlrQB2Y8GlrRAZ94+5tBWg9fdXcGkrRMb9PwcubQWwGwsubQWwGwsueyuA2/4KLm0FiLMfTrX/ZyalrnX+r9sQwC2Glr5/W6JbDCIiIuLTbq4giIiINEVXEGJpgiAiIgKEvRBhr4Uthha+f1uiWwwiIiLi026uIKTVhOlQF5sosgKJLmFEsAOJQYQRwQ4ktnYYEexAoksYMTLuDxgGEUYEO5DoEkaMHO+f9VtLJ4Pb8skuYUSwA4lBhBHBDiS6hBHBDiS6hBHBDiS2lTAi2IFElzAi2IFElzAi2IFElzAi2IFEl6WTwW35ZJcwYuR4Y7n7Gv8ne0996yw/pFsMsdrNBEFERKQpDaTQ0MIL6wH8qdZmaIIgIiICeAFkEDxlEERERCSZ6QqCiIgIyiA0pgmCiIgI0OCl0OC1MIOgpZYPPGnVYTqkxSZhrcaCS1sB7MZCEG0FiNNYaOW2AtiNBZe2AtiNhSDaCmA3FlzaCpFxY8yhrQB2Y8GlrQB2YyGItgLYjQWXtgLYjQWXtgLYjYW20lYAu7Hg0lYAu7Hg0lYAu7Hg0lYAe/nkINoKYDcWXNoK0Pzl7kOt1GKQWO1mgiAiItKUMCHCLYzmhUmeSwiaIIiIiKAMQmNqMYiIiIiPriCIiIgQVEhRtxhERESSSiSD0MLNmpLoFkO7mSB0qKmnQ4fYlK7VWHBpK4DdWAikrQD7rbHg0lYAu7Hg0lYAu7EQRFsB7MZC3LaCw/4KLm0FsBsLLm0FsBsLQbQVwG4suLQVIuP+MZe2Auy/xkIQbYXIuP/nwKWtAHZjwaWtAHZjIV5bwWV/hSDaCmA3FlzaCuCwH84e+2sr+1e7mSCIiIg0JRzAXgzJ1GJQSFFERITvMggtfbh4/fXXOeuss8jPzycUCrFkyZImj7/00ksJhUK+x5FHHhk9ZsaMGb63Dxw40PnzoQmCiIgIkSsIQTxc1NTUMGTIEObOndus4++55x42b94cfVRUVNCtWzfOO++8mOOOPPLImOPeeOMNp/MC3WIQERFJmNGjRzN69OhmH5+Tk0NOznd5liVLlvDNN98wfvz4mOM6dOhAXl5ei85NVxBERESABi8UyANg+/btMY/aWjsA31IPP/wwxcXF9O3bN2b8448/Jj8/n0MOOYSLLrqIDRs2OD93u7mCkLpzD6mpsUlYq7Hg0laIjBt7MbSRtgLYjQWXtkJk3N9AcGkrgN1YCKKtEDnen7K29lYAt/0VXNoKYDcWXNoKYDcWgmgrgN1YcGkrgN1YaO22AtiNhSDaCmA3FlzaCmA3FlzaCmA3Fqy9FcBtf4Ug2gqR4/0/My5thci5NG8/HK8hzu/OgDUEEFJs+EdIsaCgIGZ8+vTpzJgxo0XP3dimTZt48cUXefzxx2PGCwsLWbhwIQMGDGDz5s3cdNNNnHTSSXzwwQd06RKnymVoNxMEERGR1lJRUUF29neTwowM+4+nlli0aBFdu3ZlzJgxMePfv2UxePBgCgsL6du3L0899RSXX355s59fEwQREREg7KUQbuFKiuF/rKSYnZ0dM0EImud5PPLII/z0pz8lPc5V4W917dqVww8/nHXr1jl9DGUQRERE+O4WQ0sfrWHFihWsW7euWVcEqqur+eSTT+jdu7fTx9AEQUREJEGqq6tZs2YNa9asAWD9+vWsWbMmGiqcOnUql1xyie/9Hn74YQoLCznqqKN8b/v5z3/OihUr+Oyzz3jzzTc555xzSE1N5cILL3Q6t3ZziyFlZ60vx2cFEl3CiGAHEttKGBHsQKJLGBHsQKJLGBHsQGIQYUSwl092CSOCHUh0CSOCHUh0CSOCHUgMIowIdiDRJYwIdiCxtcOIYAcSgwgjgh1IdAkjgh1IdAkjgr18sksYEexAYhBhRLADiS5hRIiz3H2N/+crpWH/NAAaC0O0hdCS53DxzjvvcMopp0T/XVpaCsC4ceNYuHAhmzdv9jUQtm3bxh/+8Afuuece8zm/+OILLrzwQr7++mt69uzJiSeeyKpVq+jZs6fTubWbCYKIiEhT9mWhI+s5XIwaNQqviR0gFy5c6BvLyclh586dcd9n8eLFTucQj24xiIiIiI+uIIiIiMA+7aVgPUey0ARBREQECBMiTEszCC17/7YkoVOdve1i5Xke06ZNo3fv3nTs2JHi4mI+/vjjxJysiIgktUTs5tiWJfQKwre7WF122WWce+65vrfPnj2be++9l0WLFtG/f39uvPFGSkpK+PDDD8mMk8SPJ7RrN6GU2CCI1VhwaStAnMZCEG0FMBsLLm0FsBsLLm0FsBsLLm0FsBsLzm0Fh+WTXdoKYDcWXNoKYDcWXNoKYDcWgmgrRMb9Yy5tBdh/jQWXtgLYjYUg2gpgNxZc2gpgNxbitRVclk92aSuA3VgIoq0AdmPBpa0AcZa732WMhVunxSCxEjpBaGoXK8/zmDNnDr/61a84++yzAXj00UfJzc1lyZIlXHDBBa15qiIikuSC2Yshea4gtNlXsn79eiorKykuLo6O5eTkUFhYyMqVK+O+X21trW8XLRERkb0Je6FAHsmizU4QKisrAcjNzY0Zz83Njb7NMnPmzOh+2Tk5Ob4dtURERGTv2uwEYV9NnTqVbdu2RR8VFRWJPiURETkAhAPYh6GlCy21JW225piXlwdAVVVVzAYTVVVVDB06NO77ZWRk7JdtNUVEJLkFs5ujJgj7Xf/+/cnLy6OsrCw6Idi+fTurV69mwoQJzs/n7a6j8dfNaiw4tRWg1fdXcGkrRMb9kyWXtgLYjQWXtgLYjQWXvRXAbX8Fl7YC2I0Fl7YC2I0Fl7YC2I2FINoKYDcWWrutAHZjwaWtEBn3/xwE0VYAu7Hg0lYAu7Fg7a0AbvsruLQVwG4sBNFWiJyL/+vo0laA5u+H44Xt73XZvxI6Qaiuro7Zn/rbXay6detGnz59mDx5MrfeeiuHHXZYtOaYn5/PmDFjEnfSIiKSlBoI0dDChY5a+v5tSUInCHvbxer666+npqaGq666iq1bt3LiiSeydOlS5zUQRERE9ka3GGIldIKwt12sQqEQN998MzfffHMrnpWIiIi02QyCiIhIa2qg5bcIAlpHt01oNxMEr3Y3Xig2bGMFEttKGBHsQKJLGBHsQKJLGBHsQKJLGDFyvP+HzmXpZHBbPtkljAh2INEljAh2INEljAh2IDGIMCLYgcTWDiOCHUh0CSOCHUgMIowIdiDRJYwI9vLJLmFEsAOJLmHEyPH+n5kgwohgBxJdwojQ/OXuPa91Qoq6xRCr3UwQREREmqLtnmMlzysRERGRwOgKgoiICOARItzCDIKnmqOIiEhy0S2GWMnzSkRERCQw7eYKgldXT+NdOM3GQhtpK4DdWHBpK4DdWHBpK4DdWHBpK4C9fHIQbQWwGwsubQWwGwsubQWwGwsubQWwGwtBtBVg/zUWXNoKYDcWXNoKYDcWgmgrgN1YiNdWcFk+2aWtAHZjwaWtAHZjIYi2AtiNBZe2QmS8ecvde3G+x4IWxHbNybTdc7uZIIiIiDTl2x0ZW/ocySJ5XomIiIgERlcQRERE0C2GxjRBEBERAcKkEG7hhfWWvn9bkjyvRERERALTbq4gePV78EKNLv3sp8ZCEG2FyLi/geDSVgC7seDSVgC7sRC3reCwv0IQbQWwGwsubQWwGwsubQWwGwsubYXIuH+srbQVwG4suLQVIuP+nwOXtgLYjYUg2gpgNxasvRXAbX8Fl7YC2I0Fl7YC2I2FINoKYDcWXNoK0Pz9cDwvzh45AWvwQjS08BZBS9+/LWk3EwQREZGmKIMQSxMEERERwAtgN0dPKymKiIhIMtMVBBEREaCBEA0t3Gyppe/flrSfCUK4AUItuGDisHxyEGFEsAOJLmFEsAOJLmHEyPH+b3hr6WRwWz45iDAi2IFElzAi2IFElzAi2IFElzAi2IHEthJGBDuQ6BJGBDuQ6BJGBDuQGEQYEezlk13CiGAHEl3CiJHj/T8zLmHEyHn4v45BhBHBDiS6hBHBYbn7ON+nQQt7Lc8QhO1fgQck3WIQERERn/ZzBUFERKQJ4QBCii19/7YkeV6JiIhIC4QJBfJw8frrr3PWWWeRn59PKBRiyZIlTR6/fPlyQqGQ71FZWRlz3Ny5c+nXrx+ZmZkUFhby1ltvuX46NEEQERFJlJqaGoYMGcLcuXOd3m/t2rVs3rw5+ujVq1f0bU8++SSlpaVMnz6dd999lyFDhlBSUsKWLVucPoZuMYiIiJCYlRRHjx7N6NGjnT9Or1696Nq1q/m2u+66iyuvvJLx48cDMH/+fJ5//nkeeeQRfvnLXzb7Y2iC0JhDWwHsxkIQbQWwGwsubQWwGwsubYXIuDHm0FYAu7EQRFsB7MaCS1sB7MaCS1sB7MaCS1sB7MZCEG0FsBsLLm0FsBsLLm0FsBsLLm0FsBsLrm0Fl+WTXdoKYDcWXNoKYDcWXNoKYDcWgmgrRMaNpZZd2goQzHL3AQoyg7B9e+zvj4yMDDIy7P8H7IuhQ4dSW1vLUUcdxYwZMzjhhBMAqKuro7y8nKlTp0aPTUlJobi4mJUrVzp9DN1iEBERCVhBQQE5OTnRx8yZMwN53t69ezN//nz+8Ic/8Ic//IGCggJGjRrFu+++C8BXX31FQ0MDubm5Me+Xm5vryynsja4giIiI8I+QYkvXQfhHSLGiooLs7O+udAV19WDAgAEMGDAg+u/jjz+eTz75hLvvvpvHHnsskI/xLU0QREREAG8fWgjWcwBkZ2fHTBD2p5EjR/LGG28A0KNHD1JTU6mqqoo5pqqqiry8PKfn1S0GERERvtvNsaWP1rZmzRp69+4NQHp6OsOHD6esrOy71xUOU1ZWRlFRkdPz6gqCiIhIglRXV7Nu3brov9evX8+aNWvo1q0bffr0YerUqWzcuJFHH30UgDlz5tC/f3+OPPJIdu/eze9+9zteffVVXn755ehzlJaWMm7cOEaMGMHIkSOZM2cONTU10VZDc7XvCYLRWHBpK4DdWAiirQB2Y8GlrQB2YyFuW8FhfwWXtgLYjYUg2gpgNxZc2gpgNxZc2gqRcf+YS1sBWn9/BZe2QmTc/3Pg0lYAu7Hg0lYAu7HgsrcCuO2v4NJWALux4NJWALux4NJWALuxEERbAezGwoHSVognESspvvPOO5xyyinRf5eWlgIwbtw4Fi5cyObNm9mwYUP07XV1dfz7v/87GzduJCsri8GDB/PKK6/EPMf555/Pl19+ybRp06isrGTo0KEsXbrUF1zcm/Y9QRAREfmHIG4RuL7/qFGj8Lz4OzwtXLgw5t/XX389119//V6fd9KkSUyaNMnpXBpTBkFERER8dAVBREQE9mkvBes5koUmCCIiIiTmFkNb1n4mCCmp0CgAZgUSXcKIkXF/wDCIMCLYgUSXMGLkeP83q7V0Mrgtn+wSRgQ7kBhEGBHsQKJLGBHsQKJLGBHsQGJbCSOCHUh0CSOCHUh0CSOCHUh0CSOCHUh0WToZ3JZPdgkjRo73/8y4hBEj5+H/OrqEEcEOJAYRRoQ4gcQgwojWcvdeGOJ8r8r+034mCCIiIk3QFYRYmiCIiIigCUJjajGIiIiIj64giIiIoCsIjWmCICIiAni0vKYYf8mjA0+7mSCE0joQCsW+XKux4NJWALuxEERbAezGgktbAezlk13aCmA3FlzaCmA3FoJoK4DdWHBpK4DdWHBpK4DdWGgrbQWwGwsubQWwGwsubQWwGwsubQWwl08Ooq0AdmPBpa0AdmPBpa0AdmPBpa0AdmMhkLYC7LfGgtku8zywX3qgdAUhljIIIiIi4tNuriCIiIg0RVcQYmmCICIigiYIjSX0FsPMmTM59thj6dKlC7169WLMmDGsXbs25pjdu3czceJEunfvTufOnRk7dixVVVUJOmMREZH2IaEThBUrVjBx4kRWrVrFsmXLqK+v5/TTT6empiZ6zJQpU3juued4+umnWbFiBZs2beLcc89N4FmLiEgy+vYKQksfySKhtxiWLl0a8++FCxfSq1cvysvL+eEPf8i2bdt4+OGHefzxxzn11FMBWLBgAUcccQSrVq3iuOOOa/bHCqWnEQrFJratxoJLWwHsxkIQbQWwGwtx2woO+yu4tBXAbiy4tBXAbiwE0VYAu7Hg0laIjPvHXNoKsP8aC0G0FSLj/saCS1sB7MaCS1sB7MZCvLaCy/4KQbQVwG4suLQVwG4suLQVwG4suLQVIuPGXgxtpK0Azd8PJ+TRKi0GzwvhtfB/8C19/7akTbUYtm2L/A+gW7duAJSXl1NfX09xcXH0mIEDB9KnTx9WrlyZkHMUERFpD9pMSDEcDjN58mROOOEEjjrqKAAqKytJT0+na9euMcfm5uZSWVlpPk9tbS21td9NNbdvt/9KFRER+b4woRYvlNTS929L2swVhIkTJ/LBBx+wePHiFj3PzJkzycnJiT4KCgoCOkMREUlmyiDEahMThEmTJvGXv/yF1157jYMPPjg6npeXR11dHVu3bo05vqqqiry8PPO5pk6dyrZt26KPioqK/XnqIiIiSSmhEwTP85g0aRLPPvssr776Kv379495+/Dhw0lLS6OsrCw6tnbtWjZs2EBRUZH5nBkZGWRnZ8c8RERE9ubbkGJLH8kioRmEiRMn8vjjj/OnP/2JLl26RHMFOTk5dOzYkZycHC6//HJKS0vp1q0b2dnZXHvttRQVFTk1GABCGZmEUhq3GPwNBJe2AtiNhSDaCpHj/d9o1t4K4La/gktbAezGgktbAezGQhBtBbAbCy5tBbAbC63dVgC7sRBEWwHsxoJLWwHsxoJLWwHsxoK1twK47a8QRFshcrz/Z8alrRA5F//X0aWtAHZjwaWtAHZjoa20FaD5++GEwikQ5+sVJC2UFCuhE4R58+YBMGrUqJjxBQsWcOmllwJw9913k5KSwtixY6mtraWkpIT777+/lc9URESSnWqOsRI6QfC8vW+MmZmZydy5c5k7d24rnJGIiIhAG6o5ioiIJJIXwC0GXUEQERFJMh7QjAvbe32OZNFuJgihzHRCKbFBQyuQ6BJGBDuQGEQYMTJujDmEEcEOJLqEEcEOJLqEEcEOJAYRRgQ7kOgSRgQ7kNjaYUSwA4lBhBHBDiS6hBHBDiS6hBHBXj7ZJYwIdiAxiDAi2IFElzAi2IFElzAi2IFElzAixAkkBhFGBDOQ6BJGhOYvdx+K8zMg+1e7mSCIiIg0JUyIkFZSjNIEQUREBLUYGmsTKymKiIhI26IrCCIiIkQWOQppoaQoTRBERESINBha3GJIohpDu5kgeB0z8VJj07FWY8GlrQB2Y8G5reCwfLJLWwHsxoJLWwHsxoJLWwHsxkIQbYXIuH/Mpa0A+6+x4NJWALuxEERbAezGgktbAezGQry2gsvyyS5tBbAbC0G0FcBuLLi0FcBuLLi0FSLjxlLLLm0FaPXlk13aCpHx5i137zUkz1/lB5J2M0EQERFpikKKsTRBEBERQROExtRiEBER4bvdHFv6cPH6669z1llnkZ+fTygUYsmSJU0e/8c//pEf/ehH9OzZk+zsbIqKinjppZdijpkxYwahUCjmMXDgQNdPhyYIIiIiiVJTU8OQIUOavSHh66+/zo9+9CNeeOEFysvLOeWUUzjrrLN47733Yo478sgj2bx5c/TxxhtvOJ+bbjGIiIiQmBbD6NGjGT16dLOPnzNnTsy/b7vtNv70pz/x3HPPMWzYsOh4hw4dyMvLczuZRtrNBCGclUG4UYvBaiy4tBXAbiy47K0AbvsruLQVwG4suLQVwG4suLQVwG4sBNFWALux0NptBbAbCy5thci4PwUeRFsB7MaCS1sB7MaCtbcCuO2v4NJWALuxEERbIXIu/q+jS1sB7MaCS1sB7MZCW2krgN1YcGkrQPP3wwlq+4i9iUwQWppBCOhkmikcDrNjxw66desWM/7xxx+Tn59PZmYmRUVFzJw5kz59+jg9d7uZIIiIiLSW7dtj/yjKyMggI8OeKLXEb37zG6qrq/nJT34SHSssLGThwoUMGDCAzZs3c9NNN3HSSSfxwQcf0KVLnF69QRkEERERvmsxtPQBUFBQQE5OTvQxc+bMwM/38ccf56abbuKpp56iV69e0fHRo0dz3nnnMXjwYEpKSnjhhRfYunUrTz31lNPz6wqCiIgI4P3j0dLnAKioqCA7+7tbeEFfPVi8eDFXXHEFTz/9NMXFxU0e27VrVw4//HDWrVvn9DF0BUFERCRg2dnZMY8gJwhPPPEE48eP54knnuDMM8/c6/HV1dV88skn9O7d2+njtJsrCA1ZHQh1iA1NWYFElzBi5Hh/oMVl6WRwWz7ZJYwIdiDRJYwIdiDRJYwIdiAxiDAi2IHE1g4jgh1IdAkjgh1IDCKMCHYg0SWMCPbyyS5hRLADiS5hxMjx/p+ZIMKIYAcSXcKIYAcSXcKIECeQ2EbCiGAHEl3CiND85e4b9rROSjERCyVVV1fH/GW/fv161qxZQ7du3ejTpw9Tp05l48aNPProo0DktsK4ceO45557KCwspLKyEoCOHTuSkxP52f/5z3/OWWedRd++fdm0aRPTp08nNTWVCy+80OncdAVBREQEvrvH0NKHg3feeYdhw4ZFK4qlpaUMGzaMadOmAbB582Y2bNgQPf7BBx9kz549TJw4kd69e0cf1113XfSYL774ggsvvJABAwbwk5/8hO7du7Nq1Sp69uzpdG7t5gqCiIhIkwK4goDj+48aNQqviW7kwoULY/69fPnyvT7n4sWLnc4hHl1BEBERER9dQRARESExKym2ZZogiIiIoN0cG2s3E4Q9ndIgLTYdazUWXNoKkXFjLIC2AtiNBZe2AtiNBZe2AtiNBZe2AtiNhSDaCrD/GgsubQWwGwsubQWwGwtBtBXAbizEayu4LJ/s0lYAu7Hg0lYAu7EQRFsB7MaCS1shMm4stezSVoD91lgIoq0QGfc3EFzaCtD85e731LfSWssSo91MEERERJrkhZxDhuZzJAlNEERERFAGoTG1GERERMRHVxBEREQg2M0YkoAmCCIiIqjF0Fi7mSDUd07BS4u9o2I1FuK2FRz2VwiirQB2Y8GlrQB2Y8GlrQB2Y8GlrRAZ94+1lbYC2I0Fl7ZCZNyfAndpK4DdWAiirQB2Y8HaWwHc9ldwaSuA3VhwaSuA3VgIoq0AdmPBpa0AdmOhtdsKYDcWgmgrgN1YcGkrQPP3w9lTr7vhidBuJggiIiJ7lUS3CFpKEwQRERF0i6ExTRBERERAIcVGdGNHREREfNrNFYT6rBTC6bHzISuQaC2dDG7LJwcRRgQ7kOgSRgQ7kOgSRgQ7kOgSRgQ7kNhWwohgBxJdwohgBxJdwohgBxKDCCOCvXyySxgR7ECiSxgxcrz/Z8YljBg5D//XMYgwItiBRJcwIsQJJLZyGBHsQGIQYUSwA4kuYUSIs9y9MdZQ11p/y4b+8WjpcySHdjNBEBERaZJuMcTQLQYRERHx0RUEERER0BWERjRBEBERAe3m2IhuMYiIiIhPu7mCUN85RDg9dmZnLZ/s0lYAu7EQRFsB7MaCS1sB7MaCS1sB7MaCS1sB7MZCEG0FsBsLLm0FsBsLLm0FsBsLLm0FsBsLrm0Fl+WTXdoKYDcWXNoKYDcWXNoKYDcWgmgrRMaNpZZd2gqw3xoLLm0FsBsLQbQVwG4suLQVIM5y953832MNaa3zV7m2e47VbiYIIiIiTVIGIYZuMYiIiIiPriCIiIiAQoqNJPQKwrx58xg8eDDZ2dlkZ2dTVFTEiy++GH377t27mThxIt27d6dz586MHTuWqqqqBJ6xiIgkq5AXzCNZJHSCcPDBBzNr1izKy8t55513OPXUUzn77LP529/+BsCUKVN47rnnePrpp1mxYgWbNm3i3HPPTeQpi4hIsvICeiSJhN5iOOuss2L+/etf/5p58+axatUqDj74YB5++GEef/xxTj31VAAWLFjAEUccwapVqzjuuOOcPlZ9JwhnNB7zfyVd2gpgNxaCaCuA3VhwaSuA3VhwaStExv1jLm0FaP39FVzaCpFxfwrcpa0AdmPBpa0AdmPBZW8FcNtfwaWtAHZjwaWtAHZjwaWtAHZjIYi2AtiNhdZuK4DdWHBpK0TG/Q2EINoKYDcWXNoKkeON/XCMH7sG+0so+1mbCSk2NDSwePFiampqKCoqory8nPr6eoqLi6PHDBw4kD59+rBy5coEnqmIiCSlbzMILX0kiYSHFN9//32KiorYvXs3nTt35tlnn2XQoEGsWbOG9PR0unbtGnN8bm4ulZWVcZ+vtraW2trvppvbt9udfxERkRiqOcZI+BWEAQMGsGbNGlavXs2ECRMYN24cH3744T4/38yZM8nJyYk+CgoKAjxbERGR9iHhE4T09HQOPfRQhg8fzsyZMxkyZAj33HMPeXl51NXVsXXr1pjjq6qqyMvLi/t8U6dOZdu2bdFHRUXFfn4FIiKSFBRSjJHwWwyNhcNhamtrGT58OGlpaZSVlTF27FgA1q5dy4YNGygqKor7/hkZGWRk+MM2ezp5hDM935jv4zuEEcEOJAYRRgQ7kOgSRgQ7kOgSRgQ7kNhWwohgBxJdwohgBxJdwohgBxJdwohgBxJdlk4Gt+WTXcKIkeP9PzMuYcTIefi/ji5hRLADiUGEESFOILGVw4hgBxJdwohgBxKDCCOCHUh0CSNC85e7D6e20v91dYshRkInCFOnTmX06NH06dOHHTt28Pjjj7N8+XJeeuklcnJyuPzyyyktLaVbt25kZ2dz7bXXUlRU5NxgEBERETcJnSBs2bKFSy65hM2bN5OTk8PgwYN56aWX+NGPfgTA3XffTUpKCmPHjqW2tpaSkhLuv//+RJ6yiIgkK62kGCOhE4SHH364ybdnZmYyd+5c5s6d20pnJCIi7VUQKyFqJUURERFJapogiIiIQEJaDK+//jpnnXUW+fn5hEIhlixZstf3Wb58OccccwwZGRkceuihLFy40HfM3Llz6devH5mZmRQWFvLWW2+5nRj7cIth3LhxXH755fzwhz90/mCJtKeTR0qjFoPVWHBpK4DdWAiirQB2Y8GlrQB2Y8GlrQB2Y6GttBXAbiy4tBXAbiy4tBXAbiy4tBXAXj45iLYC2I0Fl7YC2I0Fl7YC2I0Fl7YC2I2FQNoKsN8aCy5tBbAbCy5tBbAbC0G0FcBuLMRtK9g/BvZy952NFkOHJLpu30hNTQ1Dhgzhsssua9ZeQ+vXr+fMM8/k6quv5ve//z1lZWVcccUV9O7dm5KSEgCefPJJSktLmT9/PoWFhcyZM4eSkhLWrl1Lr169mn1uzlcQtm3bRnFxMYcddhi33XYbGzdudH0KERGRNidEALs5On7M0aNHc+utt3LOOec06/j58+fTv39/7rzzTo444ggmTZrEj3/8Y+6+++7oMXfddRdXXnkl48ePZ9CgQcyfP5+srCweeeQRp3NzniAsWbKEjRs3MmHCBJ588kn69evH6NGjeeaZZ6ivt/9iEBERaU+2b98e8/j+FgAtsXLlypg9igBKSkqiexTV1dVRXl4ec0xKSgrFxcXO+xjtUwahZ8+elJaW8r//+7+sXr2aQw89lJ/+9Kfk5+czZcoUPv744315WhERkcQJcLOmgoKCmGX/Z86cGcgpVlZWkpubGzOWm5vL9u3b2bVrF1999RUNDQ3mMU3tY2RpUc1x8+bNLFu2jGXLlpGamso///M/8/777zNo0CBmz57NlClTWvL0IiIirSfAlRQrKirIzv4ul2St8NvWOU8Q6uvr+fOf/8yCBQt4+eWXGTx4MJMnT+bf/u3fop+MZ599lssuu0wTBBERaZeys7NjJghBycvLo6qqKmasqqqK7OxsOnbsSGpqKqmpqeYxTe1jZHGeIPTu3ZtwOMyFF17IW2+9xdChQ33HnHLKKb5tmhMt3KkBOsamk63GgktbAezGQhBtBbAbCy5thci4f8ylrQD7r7EQRFshMu5Pgbu0FcBuLLi0FcBuLMRrK7jsrxBEWwHsxoJLWwHsxoJLWwHsxoJLWyEybuzF0EbaCmA3FlzaCpFx/1+cLm0FsBsLQbQVIsf7v8+svRXAbitAnP1wOvm/LuGUAL5WzXEA7MVQVFTECy+8EDO2bNmy6B5F6enpDB8+nLKyMsaMGQNE9jgqKytj0qRJTh/LeYJw9913c95555EZ55saoGvXrqxfv971qUVERBImESspVldXs27duui/169fz5o1a+jWrRt9+vRh6tSpbNy4kUcffRSAq6++mt/+9rdcf/31XHbZZbz66qs89dRTPP/889HnKC0tZdy4cYwYMYKRI0cyZ84campqGD9+vNO5OU8QfvrTn7q+i4iIiBjeeecdTjnllOi/S0tLgciaQwsXLmTz5s1s2LAh+vb+/fvz/PPPM2XKFO655x4OPvhgfve730XXQAA4//zz+fLLL5k2bRqVlZUMHTqUpUuX+oKLe9PmtnsWERFJiATcYhg1ahSeF/+drFUSR40axXvvvdfk806aNMn5lkJjmiCIiIjAAZFBaE3ai0FERER82s0VhNRO9aRkxSaOrcaCS1sB7MZCEG0FsBsLLm0FsBsLrd1WALuxEERbAezGgktbAezGgktbAezGgrW3ArjtrxBEWyFyvP9PG5e2QuRc/F9Hl7YC2I0Fl7YC2I2FttJWALux4NJWALux4NJWALuxEERbITJujDm0FaD5++GEUlpnlV5t9xyr3UwQREREmvS9lRBb9BxJQhMEERERUAahEWUQRERExEdXEERERFAGobF2M0Ho3KmW1EYZMiuQ6BJGBDuQGEQYEexAoksYEexAYmuHEcEOJAYRRgQ7kOgSRgQ7kOgSRgR7+WSXMCLYgcQgwohgBxJdwohgBxJdwohgBxJdwogQJ5AYRBgRzECiSxgR7ECiSxgR7ECiSxgR7ECicxjRYflklzAiNH+5+4ZQMFsl75VuMcTQLQYRERHxaTdXEERERJoUwC2GZLqCoAmCiIgI6BZDI7rFICIiIj66giAiIgK6gtBIu5kg5HTcRYes2MS21VhwaSuA3VgIoq0QGfePubQVYP81FlzaCmA3FoJoK4DdWHBpK4DdWIjXVnBZPtmlrQB2YyGItgLYjQWXtgLYjQWXtkJk3Fhq2aWtAK2+fLJLWyEy7m8guLQVwG4suLQVwG4suCydDG7LJ7u0FcBuLHTrtNP/sWidFoNqjrF0i0FERER8NEEQERERn3Zzi0FERKRJyiDE0ARBREQEZRAa0y0GERER8Wk3VxB6ZVWTlhWblLYaCy5tBbAbC0G0FcBuLLR2WwHsxoJLWyEy7k+BB9FWALux4NJWALuxYO2tAG77K7i0FcBuLATRVoici//r6NJWALux4NJWALux0FbaCmA3FlzaCmA3FlzaCmA3FlzaCpHj/d9nLnsrgNv+Ci5tBbD3w+nZsdp/bmH7d+d+kURXAFqq3UwQREREmqQMQgzdYhAREREfXUEQERFBIcXGNEEQEREB3WJopN1MEHIzdpCeGRv6sQKJLmFEsAOJQYQRwQ4ktnYYEexAoksYEexAYhBhRLADiS5hRLCXT3YJI4IdSHQJI0aO9/92CSKMCHYg0SWMCHYg0SWMCHECiW0kjAh2INEljAh2INEljAh2INEljBgZN8YCCCOCHUh0CSNC85e7r9tjf0/L/tVuJggiIiJN0S2GWJogiIiIgG4xNKIWg4iIiPjoCoKIiAjoCkIjmiCIiIigDEJj7WaCkJu+jcz02LSw1VhwaSuA3VgIoq0A+6+x4NJWALux4NJWALuxEERbAezGQry2gsvyyS5tBbAbCy5tBbAbC0G0FcBuLLi0FSLjxlLLLm0F2G+NhSDaCpFxfwPBpa0AdmPBpa0AdmMhblvBYfnkINoKYDcWXNoK0Pzl7nent1KLQVcQYiiDICIiIj7t5gqCiIhIk3QFIYYmCCIiIiiD0FibucUwa9YsQqEQkydPjo7t3r2biRMn0r17dzp37szYsWOpqqpK3EmKiIi0E21igvD222/zwAMPMHjw4JjxKVOm8Nxzz/H000+zYsUKNm3axLnnnpugsxQRkaTmBfRwNHfuXPr160dmZiaFhYW89dZbcY8dNWoUoVDI9zjzzDOjx1x66aW+t59xxhnO55XwWwzV1dVcdNFFPPTQQ9x6663R8W3btvHwww/z+OOPc+qppwKwYMECjjjiCFatWsVxxx3n9HF6pW2nY6OUs9VYcGkrRMb9Y22lrQB2Y8GlrRAZ96fAXdoKYDcWgmgrgN1YsPZWALf9FVzaCmA3FlzaCmA3FoJoK4DdWHBpK4DdWGjttgLYjYUg2gpgNxZc2gpgNxZc2gqR4/3fZ9beCuC2v0IQbQWwGwsubQVo/n44u9LifI8FLBG3GJ588klKS0uZP38+hYWFzJkzh5KSEtauXUuvXr18x//xj3+k7ns/h19//TVDhgzhvPPOiznujDPOYMGCBdF/Z2TY38NNSfgVhIkTJ3LmmWdSXFwcM15eXk59fX3M+MCBA+nTpw8rV66M+3y1tbVs37495iEiItIW3XXXXVx55ZWMHz+eQYMGMX/+fLKysnjkkUfM47t160ZeXl70sWzZMrKysnwThIyMjJjjDjroIOdzS+gEYfHixbz77rvMnDnT97bKykrS09Pp2rVrzHhubi6VlZVxn3PmzJnk5OREHwUFBUGftoiIJKMAbzE0/kO1ttZ/ta6uro7y8vKYP4RTUlIoLi5u8g/h73v44Ye54IIL6NSpU8z48uXL6dWrFwMGDGDChAl8/fXXzf40RM/F+T0CUlFRwXXXXcfvf/97MuNcAtwXU6dOZdu2bdFHRUVFYM8tIiJJLMAJQkFBQcwfq9Yfwl999RUNDQ3k5ubGjO/tD+FvvfXWW3zwwQdcccUVMeNnnHEGjz76KGVlZdx+++2sWLGC0aNH09DgdtsvYRmE8vJytmzZwjHHHBMda2ho4PXXX+e3v/0tL730EnV1dWzdujXmKkJVVRV5eXlxnzcjI2Of7rWIiIgEpaKiguzs7/JU++P/Sw8//DBHH300I0eOjBm/4IILov999NFHM3jwYP7pn/6J5cuXc9pppzX7+RM2QTjttNN4//33Y8bGjx/PwIEDueGGGygoKCAtLY2ysjLGjh0LwNq1a9mwYQNFRUXOH69H6g6yOsQGnqxAoksYEexAYlsJI4IdSHQJI4IdSHQJI4IdSAwijAj28skuYUSwA4kuYcTI8f50kksYMXIe/q9jEGFEsAOJLmFEiBNIbOUwItiBxCDCiGAHEl3CiGAHEl3CiGAvn+wSRgQ7kBhEGBHsQKJLGBHsQGJeB//YztQAvseaIfSPR0ufAyA7OztmgmDp0aMHqampvvr+3v4QBqipqWHx4sXcfPPNez2nQw45hB49erBu3boDY4LQpUsXjjrqqJixTp060b179+j45ZdfTmlpKd26dSM7O5trr72WoqIi5waDiIjIXrXySorp6ekMHz6csrIyxowZA0A4HKasrIxJkyY1+b5PP/00tbW1XHzxxXv9OF988QVff/01vXv3bv7J0QZqjk25++67SUlJYezYsdTW1lJSUsL999+f6NMSEZEklIiaY2lpKePGjWPEiBGMHDmSOXPmUFNTw/jx4wG45JJL+MEPfuDLMDz88MOMGTOG7t27x4xXV1dz0003MXbsWPLy8vjkk0+4/vrrOfTQQykpKXE6tzY1QVi+fHnMvzMzM5k7dy5z585NzAmJiIjsR+effz5ffvkl06ZNo7KykqFDh7J06dJocHHDhg2kpMTenlq7di1vvPEGL7/8su/5UlNT+etf/8qiRYvYunUr+fn5nH766dxyyy3OOYg2NUEQERFJmARt1jRp0qS4txQa/+EMMGDAADzP/kAdO3bkpZdecj8JgyYIIiIi30qizZZaqt1MEHqkVtMpNfYyjdVYcGkrgN1YCKKtAHZjwaWtAHZjwaWtAHZjwaWtAHZjwbWt4LJ8sktbAezGgktbAezGgktbAezGQhBthci4sdSyS1sB9ltjwaWtAHZjIYi2AtiNBZe2AtiNhbhtBYflk13aCmA3FoJoK4DdWHBpK4C93H231GrfWEZqnOqV7FftZoIgIiLSFG33HEsTBBEREUhYBqGtSvhmTSIiItL26AqCiIgIusXQmCYIIiIioFsMjbSbCUJOSh2dUxq3GPzHubQVoPX3V3BpK0TG/Slwl7YC2I0Fl7YC2I0Fl70VwG1/BZe2AtiNBZe2AtiNBZe2AtiNhSDaCmA3Flq7rQB2Y8GlrRAZ9zcQgmgrgN1YcGkrRI73f59ZeyuA2/4KLm0FsBsLQbQVwG4suLQVoPn74XRIUYshEdrNBEFERKQpusUQSxMEERER0C2GRjRBEBERAU0QGlHNUURERHzazRWEnBTo0mg6ZAUS20oYEexAoksYEexAoksYEexAoksYEexAosvSyeC2fLJLGDFyvH/a7xJGjJyH/+voEkYEO5AYRBgR4gQSWzmMCHYg0SWMCHYgMYgwItiBRJcwYmTcGHMII4IdSHQJI4IdSAwijAh2INEljAjNX+4+tZX+lFUGIVa7mSCIiIg0SbcYYugWg4iIiPjoCoKIiAgQ8jxCXssuAbT0/dsSTRBERERAtxga0S0GERER8Wk3VxA6hlLJCsXOh6zGQltpK4DdWHBpK4DdWHBpK4DdWHBpK4C9fHIQbQWwGwsubQWwGwsubQWwGwsubQWwGwuBtBVgvzUWXNoKYDcWXNoKYDcWgmgrgN1YiNtWcFg+2aWtAHZjwaWtAHZjIYi2AtiNBZe2QmTcP2a1y/aE7M9/0NRiiNVuJggiIiJN0i2GGLrFICIiIj66giAiIoJuMTSmCYKIiAjoFkMjmiCIiIigKwiNtZsJQkYojYxGLYb91VgIoq0QGfenwF3aCmA3FlzaCmA3FuK1FVz2VwiirQB2Y8GlrQB2Y8GlrQB2Y8GlrRAZN/ZiaCNtBbAbCy5thci4v4Hg0lYAu7EQRFshcrz/+8zaWwHc9ldwaSuA3VhwaSuA3VgIoq0AdmPBpa0Azd8PJyMU55ek7FftZoIgIiLSJN1iiKEJgoiIyD8k0y2CllLNUURERHx0BUFERATA8yKPlj5HktAEQUREBLUYGms3E4S0UCppoX2/o+Kyv0IQbQWwGwsubQWwGwsubQWwGwvW3grgtr9CEG2FyPH+n0iXtkLkXPxfR5e2AtiNBZe2AtiNhbbSVgC7seDSVgC7seDSVgC7sRBEWwHs/RVc2gpgNxZc2gpgNxZc2gpgNxaCaCuA3VhwaStA8/fDSWulvRgkVruZIIiIiDRJLYYYmiCIiIgAoXDk0dLnSBZqMYiIiIiPriCIiIiAbjE0oglCIy5hRLADiUGEEcEOJLqEEcEOJLqEEcFePtkljAh2IDGIMCLYgUSXMCLYgUSXMCLYgUSXMCLECSQGEUYEM5DoEkYEO5DoEkYEO5DoEkYEO5DoHEZ0WD7ZJYwIdiDRJYwIdiDRJYwIdiAxiDBiZNw/5hJGhGCWuw9SoloMc+fO5Y477qCyspIhQ4Zw3333MXLkSPPYhQsXMn78+JixjIwMdn/vd43neUyfPp2HHnqIrVu3csIJJzBv3jwOO+wwp/PSLQYRERH4bh2Elj4cPPnkk5SWljJ9+nTeffddhgwZQklJCVu2bIn7PtnZ2WzevDn6+Pzzz2PePnv2bO69917mz5/P6tWr6dSpEyUlJTGTiObQBEFERCRB7rrrLq688krGjx/PoEGDmD9/PllZWTzyyCNx3ycUCpGXlxd95ObmRt/meR5z5szhV7/6FWeffTaDBw/m0UcfZdOmTSxZssTp3DRBEBER4btbDC19AGzfvj3mUVvrvwVZV1dHeXk5xcXF0bGUlBSKi4tZuXJl3POsrq6mb9++FBQUcPbZZ/O3v/0t+rb169dTWVkZ85w5OTkUFhY2+ZwWTRBERETgu5BiSx9AQUEBOTk50cfMmTN9H+6rr76ioaEh5goAQG5uLpWVleYpDhgwgEceeYQ//elP/Nd//RfhcJjjjz+eL774AiD6fi7PGY9CiiIiIgGrqKggO/u7kHhGhh3mdVVUVERRUVH038cffzxHHHEEDzzwALfccksgH+Nb7XqCYDUWXNoKYDcWgmgrgN1YcGkrgN1YiNdWcFk+2aWtAHZjIYi2AtiNBZe2AtiNBZe2QmTcWGrZpa0Arb58sktbITLu/yXn0lYAu7Hg0lYAu7HgsnQyuC2f7NJWALux4NJWALux4NJWALuxEERbAezGwoHSVognyBZDdnZ2zATB0qNHD1JTU6mqqooZr6qqIi8vr1kfLy0tjWHDhrFu3TqA6PtVVVXRu3fvmOccOnRoM19FhG4xiIiIQKu3GNLT0xk+fDhlZWXRsXA4TFlZWcxVgqY0NDTw/vvvRycD/fv3Jy8vL+Y5t2/fzurVq5v9nN9q11cQREREEqm0tJRx48YxYsQIRo4cyZw5c6ipqYmudXDJJZfwgx/8IJphuPnmmznuuOM49NBD2bp1K3fccQeff/45V1xxBRBpOEyePJlbb72Vww47jP79+3PjjTeSn5/PmDFjnM5NEwQRERESs1DS+eefz5dffsm0adOorKxk6NChLF26NBoy3LBhAykp313s/+abb7jyyiuprKzkoIMOYvjw4bz55psMGjQoesz1119PTU0NV111FVu3buXEE09k6dKlZMa5dRhPQm8xzJgxg1AoFPMYOHBg9O27d+9m4sSJdO/enc6dOzN27FjfvRoREZFABNhicDFp0iQ+//xzamtrWb16NYWFhdG3LV++nIULF0b/fffdd0ePrays5Pnnn2fYsGExzxcKhbj55puprKxk9+7dvPLKKxx++OHO55XwDMKRRx4ZsyLUG2+8EX3blClTeO6553j66adZsWIFmzZt4txzz03g2YqIiLQPCb/F0KFDBzOtuW3bNh5++GEef/xxTj31VAAWLFjAEUccwapVqzjuuOOcPk6910B9o/CI1VhwaStExv0p8CDaCmA3FlzaCmA3Fqy9FcBtfwWXtgLYjYUg2gqRc/F/HV3aCmA3FlzaCmA3FtpKWwHsxoJLWwHsxoJLWwHsxoJLWyFyvP/7zGVvBXDbX8GlrQB2Y8GlrQB2Y8GlrQB2YyGItgLYjYUg2gpWu6zea509lBO1F0NblfArCB9//DH5+fkccsghXHTRRWzYsAGA8vJy6uvrY1aDGjhwIH369GlyNaja2lrfClYiIiJ7FfaCeSSJhE4QCgsLWbhwIUuXLmXevHmsX7+ek046iR07dlBZWUl6ejpdu3aNeZ+9rQY1c+bMmNWrCgoK9vOrEBGRpJCgDEJbldBbDKNHj47+9+DBgyksLKRv37489dRTdOzYcZ+ec+rUqZSWlkb/vX37dk0SREREHCX8FsP3de3alcMPP5x169aRl5dHXV0dW7dujTlmbytMZWRkRFewas5KViIiIgAhAtisKdEvIkAJDyl+X3V1NZ988gk//elPGT58OGlpaZSVlTF27FgA1q5dy4YNG5xXg4JIILHWi50PWYFElzAi2IHEIMKIYAcSXcKIYC+f7BJGBDuQ6BJGjBzvv+4WRBgR7ECiSxgR7ECiSxgR4gQS20gYEexAoksYEexAoksYEexAoksYEezlk4MII4IdSHQJI4IdSHQJI4IdSHQJI4IdSAwijAj7L5BohcdrWymk6LoSYtznSBIJnSD8/Oc/56yzzqJv375s2rSJ6dOnk5qayoUXXkhOTg6XX345paWldOvWjezsbK699lqKioqcGwwiIiLiJqEThC+++IILL7yQr7/+mp49e3LiiSeyatUqevbsCUQWhEhJSWHs2LHU1tZSUlLC/fffn8hTFhGRJKWaY6yEThAWL17c5NszMzOZO3cuc+fObaUzEhGRdiuIFkISTRDaVEhRRERE2oY2FVIUERFJlJDnEWphyLCl79+WtJsJwi6vgQ6NvnBWY8GlrQB2YyGItgLYjYV4bQWX5ZNd2gpgNxZc2gpgNxaCaCuA3VhwaStExo2lll3aCrDfGgtBtBUi4/4GgktbAezGgktbAezGQty2gsPyyUG0FcBuLLi0FcBuLLi0FcBuLLi0FSLj/rG20laA5i93v6u1Wgzhfzxa+hxJQrcYRERExKfdXEEQERFpim4xxNIEQUREBNRiaEQTBBEREdBKio0ogyAiIiI+7eYKwrYwNDRKl1qNBZe2AtiNhSDaCmA3Fqy9FcBtfwWXtgLYjQWXtgLYjYUg2gpgNxZc2gpgNxZau60AdmMhiLYC2I0Fl7YC2I0Fl7ZC5Hj/95m1twK47a8QRFsB7MaCS1sB7MaCS1sB7MaCS1sB7MZCW2krQPP3w4nzayVwWkkxVruZIIiIiDRJtxhi6BaDiIiI+OgKgoiICBAKRx4tfY5koQmCiIgI6BZDI+1mgrAtnM6ecOwdFSuQ6BJGBDuQGEQYEezlk13CiGAHEl3CiJHj/d/wLmHEyHn4Q0pBhBHBDiS6hBEhTiCxlcOIYAcSgwgjgh1IdAkjgh1IdAkjRsaNMYcwItiBxCDCiGAHEl3CiGAHEl3CiGAHEl3CiGAHEoMII4IdSHQJI0Lzl7uvCSfRn+UHkHYzQRAREWmSFkqKoQmCiIgIWmq5MbUYRERExEdXEEREREAhxUY0QRAREYFIfqClecjkmR+0nwnCVw2d2dkQm961GgsubQWwGwuubQWX5ZNd2gpgNxZc2gpgNxZc2gpgNxaCaCtExo2lll3aCrDfGgsubQWwGwtBtBXAbiy4tBXAbizEbSs4LJ/s0lYAu7EQRFsB7MaCS1sB7MaCS1shMu4fc2krQOsvn+zSVoiMN2+5+50NAfx8NoMyCLGUQRARERGfdnMFQUREpEkeAWQQAjmTNkETBBEREVBIsRHdYhAREREfXUEQERGBSIPBzty6PUeSaDcThK8autBxT+zLtRoLLm0FsBsLLnsrgNv+Ci5tBbAbCy5tBbAbCy5tBbAbC0G0FcBuLLR2WwHsxoJLWyEy7m8gBNFWALux4NJWiBzv/z6z9lYAt/0VXNoKYDcWgmgrgN1YcGkrgN1YcGkrgN1YaCttBbAbCy5tBWj+fji7GuL8PAcsUS2GuXPncscdd1BZWcmQIUO47777GDlypHnsQw89xKOPPsoHH3wAwPDhw7nttttijr/00ktZtGhRzPuVlJSwdOlSp/PSLQYREZEEefLJJyktLWX69Om8++67DBkyhJKSErZs2WIev3z5ci688EJee+01Vq5cSUFBAaeffjobN26MOe6MM85g8+bN0ccTTzzhfG6aIIiIiMB3IcWWPhzcddddXHnllYwfP55BgwYxf/58srKyeOSRR8zjf//733PNNdcwdOhQBg4cyO9+9zvC4TBlZWUxx2VkZJCXlxd9HHTQQc6fDk0QREREoNUnCHV1dZSXl1NcXBwdS0lJobi4mJUrVzbrOXbu3El9fT3dunWLGV++fDm9evViwIABTJgwga+//rrZ5/WtdpNBEBERaS3bt8fmUzIyMsjIiM0NffXVVzQ0NJCbmxsznpuby//93/816+PccMMN5Ofnx0wyzjjjDM4991z69+/PJ598wn/8x38wevRoVq5cSWpq8/Mq7WaCsKU+m8z62ICPFUh0CSOCHUh0WToZ3JZPdgkjRo73z2ZdwoiR8/CHlFzCiGAHEoMII0KcQGIrhxHBDiS6hBHBDiQGEUYEO5DoEkYEe/lklzAi2IFElzAi2IHEIMKIYAcSXcKIYAcSXcKIYAcS20oYEexAoksYEZq/3P3uevvcAhfgOggFBQUxw9OnT2fGjBkte+5GZs2axeLFi1m+fDmZ3/tdc8EFF0T/++ijj2bw4MH80z/9E8uXL+e0005r9vO3mwmCiIhIkwKsOVZUVJCd/d0EqPHVA4AePXqQmppKVVVVzHhVVRV5eXlNfpjf/OY3zJo1i1deeYXBgwc3eewhhxxCjx49WLdundMEQRkEERERvqs5tvQBkJ2dHfOwJgjp6ekMHz48JmD4beCwqKgo7nnOnj2bW265haVLlzJixIi9vq4vvviCr7/+mt69ezt9PjRBEBERSZDS0lIeeughFi1axEcffcSECROoqalh/PjxAFxyySVMnTo1evztt9/OjTfeyCOPPEK/fv2orKyksrKS6upqAKqrq/nFL37BqlWr+OyzzygrK+Pss8/m0EMPpaSkxOncdItBREQEErIXw/nnn8+XX37JtGnTqKysZOjQoSxdujQaXNywYQMpKd/9LT9v3jzq6ur48Y9/HPM832YcUlNT+etf/8qiRYvYunUr+fn5nH766dxyyy3mVYymaIIgIiICEPYg1MIJQtj9/SdNmsSkSZPMty1fvjzm35999lmTz9WxY0deeukl53OwtJsJQlVdDul1salgq7Hg0lYAe/nkINoKYDcWXNoKYDcWXNoKYDcWXNoKYDcWAmkrwH5rLLi0FcBuLLi0FcBuLATRVgC7sRC3reCwfLJLWwHsxoJLWwHsxkIQbQWwGwsubYXIuH/Mpa0A+6+xEERbITLu/zlwaStA85e7r6trpRaDxGg3EwQREZEmabvnGJogiIiIABDABIHkmSCoxSAiIiI+uoIgIiICusXQiCYIIiIi8I8GQuu3GNqqdjNBqKrtQlqH2NSt1ViI11Zw2V8hiLYC2I0Fl7YC2I0Fl7YC2I0Fl7ZCZNzYi6GNtBXAbiy4tBUi4/4GgktbAezGQhBthcjx/u8za28FcNtfwaWtAHZjwaWtAHZjIYi2AtiNBZe2AtiNhdZuK4DdWAiirQB2Y8GlrQDN3w+nvtb+/Mv+1W4mCCIiIk3ywpFHS58jSWiCICIiAsogNKIJgoiICCiD0EjCa44bN27k4osvpnv37nTs2JGjjz6ad955J/p2z/OYNm0avXv3pmPHjhQXF/Pxxx8n8IxFRESSX0InCN988w0nnHACaWlpvPjii3z44YfceeedHHTQQdFjZs+ezb333sv8+fNZvXo1nTp1oqSkhN1xlukVERHZJ9/eYmjpI0kk9BbD7bffTkFBAQsWLIiO9e/fP/rfnucxZ84cfvWrX3H22WcD8Oijj5Kbm8uSJUu44IILmv2xvtrViQ4psclxq7Fg7a0AbvsrBNFWiBzv/0ZzaStEzsWfYnZpK4DdWHBpK4DdWGgrbQWwGwsubQWwGwsubQWwGwtBtBUi48aYQ1sB7MaCS1sB7MaCS1sB7MZCEG0FsBsLLm0FsBsLrd1WALuxEERbAezGgktbAezGQtVO/8fbE+d3U+A8AsggBHImbUJCryD8+c9/ZsSIEZx33nn06tWLYcOG8dBDD0Xfvn79eiorKykuLo6O5eTkUFhYyMqVK83nrK2tZfv27TEPERERcZPQCcKnn37KvHnzOOyww3jppZeYMGECP/vZz1i0aBEAlZWVANF9sb+Vm5sbfVtjM2fOJCcnJ/ooKCjYvy9CRESSg24xxEjoBCEcDnPMMcdw2223MWzYMK666iquvPJK5s+fv8/POXXqVLZt2xZ9VFRUBHjGIiKStMLhYB5JIqEThN69ezNo0KCYsSOOOIINGzYAkJeXB0BVVVXMMVVVVdG3NZaRkUF2dnbMQ0RERNwkNKR4wgknsHbt2pixv//97/Tt2xeIBBbz8vIoKytj6NChAGzfvp3Vq1czYcIEp4/1zc4sUokNjFnLJ7uEEcEOJAYRRgQ7kOgSRgQ7kOgSRgQ7kOgSRoQ4gcQgwohgBhJdwohgBxJdwohgBxJdwohgBxKdw4gOyye7hBHBDiS6hBHBDiS6hBHBDiQGEUaMjPvHXMKIsP8CiS5hRLADiUGEEcEOJLqEEcFe7n7rzo6+sYZdrfS3rBZKipHQCcKUKVM4/vjjue222/jJT37CW2+9xYMPPsiDDz4IQCgUYvLkydx6660cdthh9O/fnxtvvJH8/HzGjBmTyFMXEZFkowlCjIROEI499lieffZZpk6dys0330z//v2ZM2cOF110UfSY66+/npqaGq666iq2bt3KiSeeyNKlS8mMU0ETERGRlkv4Usv/8i//wr/8y7/EfXsoFOLmm2/m5ptvbsWzEhGRdkdLLcdI+ARBRESkLfC8MF4Ld2Ns6fu3JZogiIiIQCQ/0NIrAMogHHhqajJI8WJzC1ZjwaWtAHZjIYi2AtiNBZe2AtiNBZe2QmTcWGrZpa0Arb58sktbITLubyC4tBXAbiy4tBXAbiy4LJ0Mbssnu7QVwG4suLQVwG4suLQVwG4sBNFWALux0NptBbAbCy5thci4/+cgiLYC2I0Fl7YC2MvdW+2y8E7z3WU/azcTBBERkSZ5AWQQdAVBREQkyYTDEGphhiCJMggJXUlRRERE2iZdQRAREQHdYmhEEwQRERHAC4fxWniLQTXHA9CemjRSwrFJZKux4NJWALuxEERbIXIu/hSzS1sB7MaCS1sB7MZCW2krgN1YcGkrgN1YcGkrgN1YcGkrRI73f5+57K0AbvsruLQVwG4suLQVwG4suLQVwG4sBNFWALux0NptBbAbCy5tBbAbC0G0FcBuLLi0FaD5++GEdwW0b4s4aTcTBBERkSbpFkMMTRBEREQgskhSSBOEb6nFICIiIj66giAiIgL/+Ou/pesgJM8VhHYzQUipSSWlITZo1MFaJtkhjAiQVu3/ZkivaXkYEexAYkqNHR4L7YoXPPSPu4QRIU4gsY2EEcEOJLqEEcEOJLqEEcEOJNbFCSnWd4qzfLIRSIwbRuwcJ4zYqfnLJ8cLI3brZK9t27NjtW/MJYwIdiAxr4N9bLdU/8cDO5AYRBgR9l8g0SWMCHYg0SWMCHYgMYgwIkDVTv/H3Lqzo3msFUaEOMvdVxuf/90t/5o0hxf28Fp4i8HTBEFERCTJeGFafgUheWqOyiCIiIiIj64giIiIoFsMjWmCICIiArrF0EjSTxC+nc2FjWBeeLc/KNZQa4fH4izSRkqdf7a4p97+BgnFGWePHV7yGvzhwJSGOCsmhuOEFMPGKohenDBinBCV5xkhxTjBKidxfpBCcWbg8Sb2obD/Tlm81VK9Bvvra2UuG/bYr3FPfbxx/3k01Nl38RrS7POwvrzh1DhhxA5xxlPs8wul+L++DSH7+2YP9ni98f1UF+f7d3e6Pb4rzf/9tDPVPueM1DjB3hT/eGqcG6Z7QvbnOiPON0lanONd1Bvf27Vxvt93xRm3Ms01YfvYnQ3252+X8Ttkd739damrs8fra+3fF3uM1VsbdtlfhLCdebVXSDQCid/+/t7ff53vob7F6yTtwf48HohCXjJdDzF88cUXFBQUJPo0RESkhSoqKjj44IMDf97du3fTv39/KisrA3m+vLw81q9fT2acJd8PFEk/QQiHw2zatIkuXbqwY8cOCgoKqKioIDvbru4c6LZv357UrzHZXx/oNSYLvcbgeJ7Hjh07yM/PJyVl/2Trd+/eTV2cqrer9PT0A35yAO3gFkNKSkp0xhn6x6XD7OzspP2B/Vayv8Zkf32g15gs9BqDkZNjr9MQlMzMzKT4n3qQVHMUERERH00QRERExKddTRAyMjKYPn06GRn2MrzJINlfY7K/PtBrTBZ6jXKgS/qQooiIiLhrV1cQREREpHk0QRAREREfTRBERETERxMEERER8Wk3E4S5c+fSr18/MjMzKSws5K233kr0Ke2z119/nbPOOov8/HxCoRBLliyJebvneUybNo3evXvTsWNHiouL+fjjjxNzsvto5syZHHvssXTp0oVevXoxZswY1q5dG3PM7t27mThxIt27d6dz586MHTuWqqqqBJ2xu3nz5jF48ODoIjNFRUW8+OKL0bcf6K+vsVmzZhEKhZg8eXJ07EB/jTNmzCAUCsU8Bg4cGH37gf76vrVx40YuvvhiunfvTseOHTn66KN55513om9Pht854tcuJghPPvkkpaWlTJ8+nXfffZchQ4ZQUlLCli1bEn1q+6SmpoYhQ4Ywd+5c8+2zZ8/m3nvvZf78+axevZpOnTpRUlLCbmPDqrZqxYoVTJw4kVWrVrFs2TLq6+s5/fTTqampiR4zZcoUnnvuOZ5++mlWrFjBpk2bOPfccxN41m4OPvhgZs2aRXl5Oe+88w6nnnoqZ599Nn/729+AA//1fd/bb7/NAw88wODBg2PGk+E1HnnkkWzevDn6eOONN6JvS4bX980333DCCSeQlpbGiy++yIcffsidd97JQQcdFD0mGX7niMFrB0aOHOlNnDgx+u+GhgYvPz/fmzlzZgLPKhiA9+yzz0b/HQ6Hvby8PO+OO+6Ijm3dutXLyMjwnnjiiQScYTC2bNniAd6KFSs8z4u8prS0NO/pp5+OHvPRRx95gLdy5cpEnWaLHXTQQd7vfve7pHp9O3bs8A477DBv2bJl3sknn+xdd911nuclx9dw+vTp3pAhQ8y3JcPr8zzPu+GGG7wTTzwx7tuT9XeOeF7SX0Goq6ujvLyc4uLi6FhKSgrFxcWsXLkygWe2f6xfv57KysqY15uTk0NhYeEB/Xq3bdsGQLdu3QAoLy+nvr4+5nUOHDiQPn36HJCvs6GhgcWLF1NTU0NRUVFSvb6JEydy5plnxrwWSJ6v4ccff0x+fj6HHHIIF110ERs2bACS5/X9+c9/ZsSIEZx33nn06tWLYcOG8dBDD0Xfnqy/c6Qd3GL46quvaGhoIDc3N2Y8Nzc3sK0925JvX1Myvd5wOMzkyZM54YQTOOqoo4DI60xPT6dr164xxx5or/P999+nc+fOZGRkcPXVV/Pss88yaNCgpHl9ixcv5t1332XmzJm+tyXDaywsLGThwoUsXbqUefPmsX79ek466SR27NiRFK8P4NNPP2XevHkcdthhvPTSS0yYMIGf/exnLFq0CEjO3zkSkfS7OcqBb+LEiXzwwQcx93aTxYABA1izZg3btm3jmWeeYdy4caxYsSLRpxWIiooKrrvuOpYtW5a0u+SNHj06+t+DBw+msLCQvn378tRTT9GxY8cEnllwwuEwI0aM4LbbbgNg2LBhfPDBB8yfP59x48Yl+Oxkf0r6Kwg9evQgNTXVlxyuqqoiLy8vQWe1/3z7mpLl9U6aNIm//OUvvPbaa9FtuyHyOuvq6ti6dWvM8Qfa60xPT+fQQw9l+PDhzJw5kyFDhnDPPfckxesrLy9ny5YtHHPMMXTo0IEOHTqwYsUK7r33Xjp06EBubu4B/xob69q1K4cffjjr1q1Liq8hQO/evRk0aFDM2BFHHBG9lZJsv3PkO0k/QUhPT2f48OGUlZVFx8LhMGVlZRQVFSXwzPaP/v37k5eXF/N6t2/fzurVqw+o1+t5HpMmTeLZZ5/l1VdfpX///jFvHz58OGlpaTGvc+3atWzYsOGAep2NhcNhamtrk+L1nXbaabz//vusWbMm+hgxYgQXXXRR9L8P9NfYWHV1NZ988gm9e/dOiq8hwAknnOCrGP/973+nb9++QPL8zhFDolOSrWHx4sVeRkaGt3DhQu/DDz/0rrrqKq9r165eZWVlok9tn+zYscN77733vPfee88DvLvuust77733vM8//9zzPM+bNWuW17VrV+9Pf/qT99e//tU7++yzvf79+3u7du1K8Jk334QJE7ycnBxv+fLl3ubNm6OPnTt3Ro+5+uqrvT59+nivvvqq984773hFRUVeUVFRAs/azS9/+UtvxYoV3vr1672//vWv3i9/+UsvFAp5L7/8sud5B/7rs3y/xeB5B/5r/Pd//3dv+fLl3vr1673/+Z//8YqLi70ePXp4W7Zs8TzvwH99nud5b731ltehQwfv17/+tffxxx97v//9772srCzvv/7rv6LHJMPvHPFrFxMEz/O8++67z+vTp4+Xnp7ujRw50lu1alWiT2mfvfbaax7ge4wbN87zvEjt6MYbb/Ryc3O9jIwM77TTTvPWrl2b2JN2ZL0+wFuwYEH0mF27dnnXXHONd9BBB3lZWVneOeec423evDlxJ+3osssu8/r27eulp6d7PXv29E477bTo5MDzDvzXZ2k8QTjQX+P555/v9e7d20tPT/d+8IMfeOeff763bt266NsP9Nf3reeee8476qijvIyMDG/gwIHegw8+GPP2ZPidI37a7llERER8kj6DICIiIu40QRAREREfTRBERETERxMEERER8dEEQURERHw0QRAREREfTRBERETERxMEERER8dEEQURERHw0QRAREREfTRBEDkBffvkleXl53HbbbdGxN998k/T09Jhd9URE9pX2YhA5QL3wwguMGTOGN998kwEDBjB06FDOPvts7rrrrkSfmogkAU0QRA5gEydO5JVXXmHEiBG8//77vP3222RkZCT6tEQkCWiCIHIA27VrF0cddRQVFRWUl5dz9NFHJ/qURCRJKIMgcgD75JNP2LRpE+FwmM8++yzRpyMiSURXEEQOUHV1dYwcOZKhQ4cyYMAA5syZw/vvv0+vXr0SfWoikgQ0QRA5QP3iF7/gmWee4X//93/p3LkzJ598Mjk5OfzlL39J9KmJSBLQLQaRA9Dy5cuZM2cOjz32GNnZ2aSkpPDYY4/x3//938ybNy/RpyciSUBXEERERMRHVxBERETERxMEERER8dEEQURERHw0QRAREREfTRBERETERxMEERER8dEEQURERHw0QRAREREfTRBERETERxMEERER8dEEQURERHw0QRARERGf/x980ZDKrSf/iQAAAABJRU5ErkJggg==" + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "HBox(children=(HTML(value=\".\\\\tmp__kuzitk.pdf
\"), HTML(val…" + ], + "application/vnd.jupyter.widget-view+json": { + "version_major": 2, + "version_minor": 0, + "model_id": "c75b8c33ee774761ad1c8d8fbe1eee33" + } + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "execution_count": 20 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.700641Z", + "start_time": "2024-10-30T19:30:40.685639Z" + } + }, + "cell_type": "code", + "source": [ + "residual = pympdata_result_state - trixi_output\n", + "rmse = np.sqrt(np.mean(residual**2))\n", + "mse = np.mean(residual**2)\n", + "max_diff = np.max(np.abs(residual))\n", + "min_diff = np.min(np.abs(residual))" + ], + "id": "ec2fffd2627288b4", + "outputs": [], + "execution_count": 21 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.716644Z", + "start_time": "2024-10-30T19:30:40.702644Z" + } + }, + "cell_type": "code", + "source": [ + "assert np.allclose(rmse, 6.94e-2, 0.1)\n", + "assert np.allclose(mse, 4.81e-3, 0.1)\n", + "assert np.allclose(max_diff, 0.285, 0.1)\n", + "assert np.allclose(min_diff, 2.69e-5, 0.1)" + ], + "id": "a7496a6f898495a3", + "outputs": [], + "execution_count": 22 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2024-10-30T19:30:40.732655Z", + "start_time": "2024-10-30T19:30:40.718646Z" + } + }, + "cell_type": "code", + "source": "", + "id": "119bd01509fa1ceb", + "outputs": [], + "execution_count": 22 + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/docs/pympdata_examples_landing.md b/examples/docs/pympdata_examples_landing.md index ed60f6fb..160d53e4 100644 --- a/examples/docs/pympdata_examples_landing.md +++ b/examples/docs/pympdata_examples_landing.md @@ -28,7 +28,8 @@ The examples are grouped by the dimensionality of the computational grid. | Spectral-spatial advection, particle population condensational growth in a vertical column of air, time dependent flow | `PyMPDATA_examples.Shipway_and_Hill_2012`
spectral-spatial | | shallow-water equations
$$\begin{eqnarray} \partial_t h + \partial_x (uh) + \partial_y (vh) &=& 0~ \\\ \partial_t (uh) + \partial_x ( uuh) + \partial_y (vuh) &=& - h \partial_x h~ \\\ \partial_t (vh) + \partial_x ( uvh) + \partial_y (vvh) &=& - h \partial_y h~ \end{eqnarray}$$ | `PyMPDATA_examples.Jarecka_et_al_2015`* | | advection equation, solid body rotation | `PyMPDATA_examples.Molenkamp_test_as_in_Jaruga_et_al_2015_Fig_12`* | -| advection equation, coordinate transformation, spherical coordinates, see also examples in [PyMPDATA-MPI](https://pypi.org/project/PyMPDATA-MPI/) $$ \partial_t (G \psi) + \nabla \cdot (Gu \psi) = 0 $$ | `PyMPDATA_examples.Williamson_and_Rasch_1989_as_in_Jaruga_et_al_2015_Fig_14`
mpi-gif | +| advection equation, coordinate transformation, spherical coordinates, see also examples in [PyMPDATA-MPI](https://pypi.org/project/PyMPDATA-MPI/) $$ \partial_t (G \psi) + \nabla \cdot (Gu \psi) = 0 $$ | `PyMPDATA_examples.Williamson_and_Rasch_1989_as_in_Jaruga_et_al_2015_Fig_14`
mpi-gif | +| advection equation, comparison against DG solution using [Trixi.jl](https://trixi-framework.github.io/) | `PyMPDATA_examples.trixi_comparison` | ## in 3D | tags | link | diff --git a/examples/setup.py b/examples/setup.py index 6fbb48c6..fcaad0c1 100644 --- a/examples/setup.py +++ b/examples/setup.py @@ -41,6 +41,7 @@ def get_long_description(): "joblib", "sympy", "imageio", + "meshio", ], author="https://github.com/open-atmos/PyMPDATA/graphs/contributors", license="GPL-3.0",