06aa_soil_column_FDM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)\n",
    "\n",
    "Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.\n",
    "Held at Technische Universität Bergakademie Freiberg.\n",
    "\n",
    "Comments to:\n",
    "\n",
    "*Prof. Dr. Thomas Nagel  \n",
    "Chair of Soil Mechanics and Foundation Engineering  \n",
    "Geotechnical Institute  \n",
    "Technische Universität Bergakademie Freiberg.*\n",
    "\n",
    "https://tu-freiberg.de/en/soilmechanics\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Some plot settings\n",
    "import plot_functions.plot_settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 6 - Finite differences: soil column under gravity and top load\n",
    "\n",
    "![soil_column.png](Images/soil_column.png)\n",
    "\n",
    "## Governing differential equation\n",
    "\n",
    "Consider the equilibrium conditions on an infinitesimal element of the soil column:\n",
    "\n",
    "$$\n",
    "    \\uparrow: \\quad F_z - \\frac{\\partial F_z}{\\partial z}\\text{d}z - F_z - \\varrho g A \\text{d}z = 0\n",
    "$$\n",
    "\n",
    "The vertical force is determined by $F_z = \\sigma_{zz}A = E_\\text{s} A \\epsilon_{zz} = -E_\\text{s} A u_{z,z}$. Therefore, the equilibrium conditions read:\n",
    "\n",
    "$$\n",
    "    0 = \\left[ \\frac{\\partial}{\\partial z} \\left(E_\\text{s} A \\frac{\\partial u_z}{\\partial z}\\right) - \\varrho g A \\right]\\text{d}z\n",
    "$$\n",
    "\n",
    "Considering, that the equation has to hold for an arbitrary infinitesimal non-zero $\\text{d}z$, we find the ordinary differential equation\n",
    "\n",
    "$$\n",
    "    0 = \\frac{\\partial}{\\partial z} \\left(E_\\text{s} A \\frac{\\partial u_z}{\\partial z}\\right) - \\varrho g A\n",
    "$$\n",
    "\n",
    "Nota bene: compare this to the differential equation for a rod ('Zugstab'): $(EAu')' + n = 0$. While the cross section of a rod can vary along its length, the cross-sectional area of a soil column is considered a constant, simplifying the equation further:\n",
    "\n",
    "$$\n",
    "    0 = \\frac{\\partial}{\\partial z} \\left(E_\\text{s} \\frac{\\partial u_z}{\\partial z}\\right) - \\varrho g\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's put this to test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite difference discretization\n",
    "\n",
    "We assume $E_\\text{s} = \\text{const.}$ and arrive at\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial^2 u_z}{\\partial z^2} = \\frac{\\gamma}{E_\\text{s}}\n",
    "$$\n",
    "\n",
    "In the finite difference method, we introduce a grid made up, in the one-dimensional case, of a series of points. The differentials are the evaluated by finite differences between values at these points:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial u_z}{\\partial z} \\approx \\frac{\\Delta u_z}{\\Delta z} = \\frac{u_i - u_{i-1}}{z_i - z_{i-1}} = u_i'\n",
    "$$\n",
    "\n",
    "where here, backward differences were chosen (viz. forward and central differences). We also assume a constant $\\Delta z$ here.\n",
    "\n",
    "A second derivative can likewise be approximated (now using forward differences):\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial^2 u_z}{\\partial z^2} \\approx u_i'' = \\frac{u_{i+1}' - u_i'}{\\Delta z}\n",
    "$$\n",
    "\n",
    "Substitution of the first derivatives yields:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial^2 u_z}{\\partial z^2} \\approx u_i'' = \\frac{u_{i+1} - 2 u_i + u_{i-1}}{\\Delta z^2}\n",
    "$$\n",
    "\n",
    "Our differential equation has now transformed into a finite difference equation evaluated at a series of points $i$ as:\n",
    "\n",
    "$$\n",
    "u_{i+1} - 2 u_i + u_{i-1} = \\Delta z^2 \\frac{\\gamma_i}{E_{\\text{s},i}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We collect all solutions $u_i$ and $z_i$ in a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = 10.\n",
    "number_of_points = 10\n",
    "Delta_z = H/(number_of_points - 1)\n",
    "z = np.linspace(0,H,number_of_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right-hand side is identical in every point, unless the stiffness or the specific weight changes with depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.00019523, 0.00019523, 0.00019523, 0.00019523, 0.00019523,\n",
       "       0.00019523, 0.00019523, 0.00019523, 0.00019523, 0.00019523])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma = 2600. * (1.-0.38) * 9.81 # N/m³\n",
    "Es = 1.e8 #Pa\n",
    "RHS = np.ones(number_of_points) * gamma/Es * Delta_z**2\n",
    "RHS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left hand side can be expressed via the finite difference matrix $A_{ij}$, so that the equation system finally can be written as $A_{ij} u_{j} = b_i$ for a system with $n=10$ points:\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "    -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
    "    1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
    "    0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
    "    0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0\\\\\n",
    "    0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\\\\n",
    "    0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0\\\\\n",
    "    0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0\\\\\n",
    "    0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 & 0\\\\\n",
    "    0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1\\\\\n",
    "    0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2\\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "    u_0 \\\\ u_1 \\\\ u_2 \\\\ u_3 \\\\ u_4 \\\\ u_5 \\\\ u_6 \\\\ u_7 \\\\ u_8 \\\\ u_9\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\Delta z^2\n",
    "\\begin{pmatrix}\n",
    "    \\frac{\\gamma_0}{E_{\\text{s},0}} \\\\ \\frac{\\gamma_1}{E_{\\text{s},1}} \\\\ \\frac{\\gamma_3}{E_{\\text{s},3}} \\\\ \\frac{\\gamma_4}{E_{\\text{s},4}} \\\\ \\frac{\\gamma_5}{E_{\\text{s},5}} \\\\ \\frac{\\gamma_6}{E_{\\text{s},6}} \\\\ \\frac{\\gamma_7}{E_{\\text{s},7}} \\\\ \\frac{\\gamma_8}{E_{\\text{s},8}} \\\\ \\frac{\\gamma_9}{E_{\\text{s},9}}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Let's build this matrix now ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0., -2.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0., -2.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0., -2.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0., -2.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0., -2.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0., -2.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -2.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -2.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -2.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Central_Diag = np.diag([-2.]*number_of_points,0)\n",
    "Central_Diag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Upper_Diag = np.diag([1.]*(number_of_points-1),1)\n",
    "Lower_Diag = np.diag([1.]*(number_of_points-1),-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Central_Diag + Upper_Diag + Lower_Diag\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have everything we need: we have the load vector (right-hand side), the system matrix (vertex connectivity and second derivative), and can solve for the vector of unknown vertex displacements. Let's give it a shot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.linalg.solve(A,RHS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAHvCAYAAACIbWPJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABQ/klEQVR4nO3deXzcV33v//eZGe3SaLU2S7IkJ463LLackJDNaqGllEJuA1z2S1sIFMrWAJdfaUtbLi1tEqDQXihd4LaBUvZSCoEWpIQskFi241i2Y8eSbI1Wa9+lWc7vjxmNNZJsS7JG35nR6/l4+CHrO2dGHx+Pvnr7+CzGWisAAAAA8eNyugAAAAAg1RG6AQAAgDgjdAMAAABxRugGAAAA4ozQDQAAAMSZx+kC4u1lL3uZfeSRR5wuAwAAAKnPXOqBlB/pHhgYcLoEAAAAbHIpH7oBAAAApxG6AQAAgDgjdAMAAABxRugGAAAA4ozQDQAAAMQZoRsAAACIM0I3AAAAEGeEbgAAACDOCN0AAABAnBG6AQAAgDgjdAMAAABxRugGAAAA4szx0G2MqTDGPGKMsU7XAgAAAMSDo6HbGPM/JD0lafsV2uUaY/7GGPO8MeaEMebHxpg9G1MlAAAAcHWcHun+iKSXSnriCu2+IWmfpH3W2t2SfiGp2RizNc71AQAAAFfN6dB9u7X2zOUaGGNeKullkv7IWjsVufxxSW5JfxDn+gAAAICr5mjottYGVtDsXkl+SY8veN6cwqPj98apNAAAACShOX9QE9N+p8tYwuN0AStwg6TuSNBeqF3SK4wxpdbafgfqAgAAgEOstRocnVFHz5jau0fV0T2m9p4xdV2Y0Kvu2q7f/o3EWv6XDKG7RNL4MtfHIh+LJcWEbmPMfZLuk6Sampq4FgcAAID4mpkL6Hzv+MWA3TOmju6xZUe0XUaamFo8Vuu8ZAjdl2Iu9YC19ouSvihJBw4cYCtCAACAJGCtVf/wtDq6R9UeCdYdPaPqHpiUXSbR5WWnqa4yX7UV3vCvSq9qyr3KSHNvfPFXkAyhe0BS5TLX8yIfBzewFgAAAKyDqRm/zveOq33B9JBzvWOamlm65M/tMqoqy1VtRb5qK8MBu67SqyJvpoy55DhsQkmG0H1M0gFjTPqied11kvqYzw0AAJC4QiGr3qHJ8JzryMh1R8+Yegenlm1fkJcRHbmuq/SqtiJf1WW5SvMk3uj1aiRD6P62pHdIerGkZkkyxqRHPv+ac2UBAABgoYlpv871jMVMDznXO6aZueCSth63SzVledGR6/npIYV5mQ5UHn8JH7qttT82xvxI0seNMb8a2av7o5JCkv7c2eoAAAA2n2AwpO6ByfCCxgWLGy8MTy/bvjg/c0GwzlddhVdbS3PlcTt9ZMzGcTR0G2MeUPhEyprI50cjD92yaCrJayR9UtJRY0xQkk/SQWtt1waWCwAAsOmMTc6Fp4QsmB5yvndcc4HQkrbpHpdqKryqWzByva3cq/zcDAcqTyyOhm5r7YdW2G5c0rvjXA4AAMCmFQiG1NU/EZkWcnF6yNDYzLLtSwuzYhY21lZ4VbklV25Xcixs3GgJP70EAAAA62t4fCayHd/FqSGdfeMKBJfuy5eZ7ta2+YWNkekh2yq8ys1Kc6Dy5EXoBgAASFH+QFCdfRPq6BkNTw2JBO2Ridll21cU5yxZ2FhelCMXo9dXjdANAACQ5Ky1Ghqbicy5HosciT4qX/+EQqGlo9fZmZ4lCxtryvOUncnodbwQugEAAJLIrD+o871jC6aHhBc3jk8tPRLdGGnrllzVVi5c3Jiv0sKspDlUJlUQugEAABKQtVYXhqfDwXrB7iE9AxNaZvBauVmRI9EXTA+pKc9TZjpxLxHwtwAAAOCw6dmAzkVGr+cXNnb0LH8kustlVFOWu+DUxnzVVnhVnJ88R6JvRoRuAACADRIKWfUNTV3c9zoy/7pncHLZ9vm56apbtC1fdVme0tOS+0j0zYjQDQAAEAeT0/7oiPX81nznei51JLpRVWme6iq90b2v6yq8KsjLYPQ6RRC6AQAArkIwZNUzMHFx15DIwsb+SxyJXuTNUG1FfiRghxc2bt2SqzTP5jkSfTMidAMAAKzQ+NRcdDu++ekh53vHNedfOnqd5nGppjwvZt51bQVHom9WhG4AAIBFAsGQui5MLDm1cXB0+SPRSwqyIuHaG52DXVmSI7eb0WuEEboBAMCmNjI+G17YOL/ndfeYzveNKxAMLWmbke7WtvK82OkhFV7lZqc7UDmSCaEbAABsCv5AUL7+ieipjfOj1yPjyx+JXlaUfXFqSGRhY1lxjtwciY41IHQDAICUMn8k+vzCxvmA7eufUHCZU2WyMtzhHUMi00NqK/K1rYIj0bG+CN0AACBpzfqD6uwdV0fPaHTP6/buMY1PzS1pa4xUWZIT2fP64vSQ0sJsuRi9RpwRugEAQMKz1urCyPSS0evuC8sfiZ6TlRYeuY5syVdX6VVNWZ4yM4g+cAbvPAAAkFBm5o9Enx+5jhwuMzntX9LWZaTqstwl00NKCjgSHYmF0A0AABwRCln1D09FFzZ29IyqvXtMvYOTssuMXudlp4dD9YJt+arL8pTBkehIAoRuAAAQd1MzC45E7754NPr0bGBJW7fLqDpyqMzC3UMKORIdSYzQDQAA1k0wZNU7OBlzamNHz5j6hqaWbV+YlxE9Cn1+ekhVaR5HoiPlELoBAMCaTEzNRXcMmZ8ecq53XLNzS49E97gXHokenh6yrcKrgjyORMfmQOgGAACXFZw/Er1n7OKpjT1jGhiZXrZ9SX5mdOR6PmRXbsmVhyPRsYkRugEAQNToxOyCHUPCJzae7x2XP7D0SPT0tPkj0SOLGyNBO48j0YElCN0AAGxC/kBIvv7xRQsbRzU0tvyR6KVF2eE9rxcE7HKORAdWjNANAEAKs9ZqeHw2Eqwvntro6x9XILj8kejbymMXNm4r9yoniyPRgatB6AYAIEXM+YM63zceM3Ld0TOm0Ynlj0SvKMlZcGpjePSaI9GB+CB0AwCQZKy1GhiZiYbq+TnYXRcmFFrmTPScTE/MwsbayOh1FkeiAxuG7zYAABLYzFxA53vHIzuGXAzZE5c4Er2qNPfivOvIqY1bCrI4VAZwGKEbAIAEYK1V39DUolMbR9U9cKkj0dOiu4XMh+yaci9HogMJitANAMAGm5rx61zPeMzCxssdiV5VlqvayKj1/OLGIm8mo9dAEiF0AwAQJ6HIkejtC0auO3rG1Du4/JHoBbkZMcG6tiJf1WW5SvMweg0kO0I3AADrYGJq7uLUkPmQ3Tt26SPRy/KiAXt+ekhhXqYDlQPYCIRuAABWIRgMqXtgMrJjyGg0ZF8YXv5I9CJvZmTUOrz3dV2FV1tLORId2GwI3QAAXMLoxOyShY3ne8c1t9yR6B6Xaub3vF6wLV9+boYDlQNINIRuAMCmFwyG1Nk/oY7u8Mj1/BzsobGZZdtvKcyKbsc3Pz2ksiRHbkavAVwCoRsAsOnM+oM6fX5YrW2Dam0b1KmOIc0sM/c6M92tbZFQXReZHrKtwqtcjkQHsEqEbgBAypuc9utkx1A0ZJ/pHFEgGDtFpKwoW/VbY/e9Li/K4Uh0AOuC0A0ASDkj47NqbR/UibZBHW8bVEf3qBaejm6MVFfp1Z66Yu3ZXqw9dcUq9LJzCID4IXQDAJJe/9CUWtvDo9jHzw6q68JEzONul9F1NQXaW1+sPfXF2lVbpNzsdIeqBbAZEboBAEnFWitf/0R0qkhr++CS7frS09zaua1Qe+uLtbu+WNfVFCozgx95AJzDHQgAkNCCIav27tFoyD7RPqjRibmYNjlZadpdVxSdLrJ9a4HSPOwkAiBxELoBAAnFHwjq9PkRnWgPz8c+1TGkqZlATJvCvAztri+OThfZVu5lwSOAhEboBgA4amrGr1PnhqOLHk+fH5Y/sHRnkT0LQnZFSY6MIWQDSB6EbgDAhhqbnNOJyKLH1rZBne0aVWjh1iKSasrztKc+vKvInvpilRRkOVQtAKwPQjcAIK4GR6d1/OxgdHeR873jMY+7XEbXVheEQ3Z9sXbXFcubw84iAFILoRsAsG6steoZmAxv3RcZye4bmoppk+Zx6bpthdFR7J21RcpiZxEAKY67HABgzUIhq3O9Y9GQfaJtUMPjszFtsjI82lVXFN6+r65YO2oKlOZxO1QxADiD0A0AWDF/IKSzXSNqjUwXOdE+pMlpf0wbb056dNHj7vpi1VXmy83OIgA2OUI3AOCSZuYCen7BziKnzg1rzh+MabOlMCtm0WNVaS47iwDAIoRuAEDUxNScTnQMRUP2C50jCi7aWaSqNDe66HFPXbFKi7IdqhYAkgehGwA2seGxmfCuIpHpIh09Y7ILMrbLSNur8qOj2LvrilWQl+FcwQCQpAjdALBJWGvVNzQV3R+7tW1Q3QOTMW08bpeurS7Q3u3hgL2rtkg5WWkOVQwAqYPQDQApKhSy6uwfjwnZg6MzMW0y093aWVsUnS6yo6ZQGWnsLAIA643QDQApIhgM6WzXqE60D+r42fDOIuNTczFt8rLTtDsyVWRPfbHqt+bL43Y5VDEAbB6EbgBIUnP+oE6fH46OYp86N6Tp2didRYq8mdpbX6w928Mhu7o0Ty627wOADUfoBoAkMTXj18mOoWjIPn1+RIFgKKZNRUlOOGRHfpUVZbN9HwAkAEI3ACSo0YnZcMBuD4fs9q5RLdy9zxiptsIbPYRmT32xiryZzhUMALgkQjcAJIj+4ano/tgn2gfV2TcR87jbZXRdTUF4+77txdpdW6Tc7HSHqgUArAahGwAcYK1V14UJtc6H7LZB9Q9Px7RJT3Nr57bC6FSR62oKlZnBbRsAkhF3bwDYAMGQVUf3aHS6yIm2IY1MzMa0ycn0aFddcXRO9vaqAqV52FkEAFIBoRsA4sAfCOpM50h00ePJjiFNzQRi2hTkZUSPUt+7vVg15V652VkEAFISoRsA1sH0bEDPnxuKTBUZ0vPnhjQXiN1ZpKwoOzpVZG99sSpKcthZBAA2CUI3AKzB+NRczKLHF3yjCi3cWkRSTXleeNFjJGiXFGQ5VC0AwGmEbgBYgcHR6Zjj1M/1jsc87nIZXVNdEN6+r65Yu+uKlJ+b4VC1AIBEQ+gGgEWsteoZnIyOZLe2Dap3cCqmTZrHpR01hdE9snduK1R2ZppDFQMAEh2hG8CmFwpZne8bV+vZgeh0kaGx2J1FsjI82lVXFJ0usqOmQGket0MVAwCSDaEbwKYTCIZ01je/s8iQTrQPamLaH9PGm5MenYu9p75YdRVeud1s3wcAWJukCN3GmAOSPi6pVlJA0rSkB6y133CyLgDJYWYuoNPnh9XaNqTWtgGdOjes2blgTJuSgqzo/th76otVVZrLziIAgHWT8KHbGFMr6SeSvifpN6y1AWPMOyV93RjzSmvtfzhaIICEMzHt18n2i4seX/CNKBCM3Vlk65Zc7d0eXvS4t75YpUXZDlULANgMEj50S3q5JK+kT1lrA5Jkrf2CMeYvJb1BEqEb2OSGx2d0om1Ix9sGdKJtSO09o7ILMrYxUv3W/Ogo9u66IhXmZTpXMABg00mG0D1/hFu0VhP+P1+XJFYxAZuMtVb9w9NqbRuIThfpujAZ08bjNrq2ujAasnfVFikni51FAADOSYbQ/TVJH5D0h8aYN0qakvT/ScqQ9IXlnmCMuU/SfZJUU1OzQWUCiAdrrTr7xtXaPqTWs4NqbRvQwOhMTJvMdLd2bivSnu3zO4sUKiONf5MDABKHsdZeuZXDjDGVkr4k6W5JE5JGJf22tfbRKz33wIED9tChQ3GuEMB6CQZDausejY5in2gf0tjkXEyb3Ky0mJ1F6rfmy8POIgAA511yBX7Cj3QbY65TeCHlDyQVSZqR9FpJ3zbGvMla+0Mn6wNwdeb8QZ3pHNHxtgG1nh3UqXNDmp6N3VmkyJsZPYRmb32xqsvy5HKxswgAIHkkfOhWeKvAAknvs9ZOR659zRjzPyX9P2NM5fwCSwCJb2rGr1Mdw+FFj+1Dev7csALBUEybipKc6CE0e+qLVV6czfZ9AICklgyh+3pJvgWBe95pSfdIqpN0ZqOLArAyoxOzOtE+GJ0u0tY1qtCinUVqK7wx00WKvOwsAgBILckQuvsl3WSM8Swa0d4myUoadqYsAMu5MDyt1gV7ZHf2jcc87nYZXVddcHFnkboi5WWnO1QtAAAbIxlC9+ckfUPSnxljPmqttcaYRkm/KenfrLUDzpYHbF7WWnUPTOr42UGdaB/U8bZB9Q9NxbRJ97i0s7YoegjNddsKlZmRDLceAADWT8L/5LPWftMY8zJJH5F0whgTlBSS9FFJn3W0OGCTsdaqo2dMx89GRrLbBzUyPhvTJifTo12R+dh764u1vapAaR52FgEAbG4JH7olyVr7I0k/croOYLPqvjCh5sM+Nbf41DMYexBNQV5GzKLHbRVeudlZBACAGEkRugFsvNGJWT1+tEtNh316/tzFpRNF3gzdtKM0PF1ke7EqS3LYWQQAgCsgdAOImvUH9XRrr5pbfGo51adgZJuRrAy3bru+Uo0NVbr+mi2MZAMAsEqEbmCTC4WsjrcNqLnFpyeOdWtqJrxJkMtldGBXmQ7ur9KL9pSz+BEAgKvAT1FgkzrXM6amlk49etingdGZ6PVrqgvU2FClO2/aqsI89ssGAGA9ELqBTWRwdFqPHelSc4tPbd2j0eulRdlq3F+lu/dXqbosz8EKAQBITYRuIMVNzfj18+M9amrx6diZC9HTIHOy0nTnTVt1cH+VdtUWycU8bQAA4obQDaSgYDCkI6cvqLnFp6eO92jOH5QkedwuvWh3mRobqnRgV5nSPG6HKwUAYHMgdAMpwlqrF3wjam7x6bEjXRqZuHhoze66IjU2VOv2Gys5ch0AAAcQuoEk1zc0pebDnWpu8cnXPxG9vnVLrhoPVOnufVUqL85xsEIAAEDoBpLQxNScHn+2W82HfWptG4xez89N1137qtTYUKVrqgo4tAYAgARB6AaShD8Q1KGTfWpq8emZE30KBEOSpPQ0t27dW67GhmrdtGOLPG6Xw5UCAIDFCN1AArPW6kT7kJoP+/T40S5NTPslScZIN127RY0HqnTr3gplZ6Y5XCkAALgcQjeQgHz942pu8anpsE/9Q1PR63WVXjU2VOuufVtVnJ/lYIUAAGA1CN1AghgZn9VjR31qbvHpTOdI9HpxfqYO7q/SwYZq1VZ4nSsQAACsGaEbcNDMXEC/ON6r5sM+HX6+X6HIyTVZGR7dfkOlGg9UaU99idwcXAMAQFIjdAMbLBiyeu6FC2pq8emp57o1PRs+uMbtMrp5d5kaG6p1y55yZaRxcA0AAKmC0A1skPbuUTW1+PToYZ+Gxmai16+rKdTBhirdedNW5edmOFghAACIF0I3EEcDI9N69LBPzYd96ugZi14vL87Wwf3VamyoUuWWXAcrBAAAG4HQDayzqRm/njzWraYWn547OyAbnqatvOw03XHTVv1SQ7Wu21bIwTUAAGwihG5gHQSCIR1+vl/NLT794niP5gLhg2vSPC7dsrtcBxuq1LCzTGkeDq4BAGAzInQDa2St1enzw2pu8emxo10am5yLPrZ3e7EO7q/W7TdWKjeLg2sAANjsCN3AKvUMTKr5sE/NLZ3qHpiMXq8uy1VjQ7Xu3lel0qJsBysEAACJhtANrMDY5Jwef7ZLzS0+newYil4vyMvQ3fuqdLChStu35jNPGwAALIvQDVzCnD+oZ070qamlUy2n+hQIhldEZqS7ddv1FWrcX60bry2R2808bQAAcHmEbmCBUMiqtX1QzS0+PfFslyZnApIkl5H2X1eqgw1VunVvhbIy+NYBAAArR3IAJJ3vHQvP0z7s04Xh6ej17VX5Ori/Wnft26oib6aDFQIAgGRG6MamNTw2o0ePdKn5cKfO+kaj17cUZung/iod3F+lmnKvgxUCAIBUQejGpjIzG9BTx3vU3OLT0dP9CkUOrsnJ9Oj2G7fqYEOV9tQVy+ViQSQAAFg/hG6kvGAwpGfPDKjpcKd+/lyPZuaCkiSP2+jmnWVqbKjWzbvLlJ7mdrhSAACQqgjdSEnWWrV1jaqpxafHjvg0PD4bfWxXbZEONlTpjhu3ypuT7mCVAABgsyB0I6X0D0/p0cM+NbX41Nk3Hr1eUZKjxoZqHdxfpYqSHAcrBAAAmxGhG0lvYtqvJ57tVvPhTh0/Oxi97s1J1103hedp76gp5OAaAADgGEI3kpI/EFLLqT41t/j09Ile+QMhSVK6x6UX7a3QwYYq7b+uVB4OrgEAAAmA0I2kYa3VqY5hNR3u1ONHuzQ+5ZckGSPdcE2JGhuqdNv1lcrJSnO4UgAAgFiEbiS87gsTamrxqflwp3oHp6LXt5XnqbGhWnftq9KWwiwHKwQAALg8QjcS0ujErH52tEvNLT49f344er3Im6m791epsaFKtRVe5mkDAICkQOhGwpj1B/X08V41He7U4VP9CkZOrsnKcOu26yvV2FCl66/ZIjcH1wAAgCRD6IajQiGr584OqLnFpyeOdWt6NiBJcrmMGnaWqrGhWi/aU67MDN6qAAAgeZFk4IhzPWNqaunUo4d9GhidiV6/prpAjQ1VuvOmrSrMy3SwQgAAgPVD6MaGGRyd1qOHu9R8uFPt3WPR66WFWToYObimuizPwQoBAADig9CNuJqa8eup53rU3OLTsy9ckA1P01ZOVpruuLFSjQ3V2lVbJBfztAEAQAojdGPdBYMhHTl9QU0tnfr58V7N+YOSJI/b6Obd5WpsqNKBXWVK87gdrhQAAGBjELqxbs76RvSTQ5362ZEujUzMRq/vritSY0O1br+xUnnZ6Q5WCAAA4AxCN66aPxDUl//zhL73WFv02tYtOWpsqNbd+6tUXpzjYHUAAADOI3TjqnT2jeuBhw+pvXtMbpfRy26r1S/fXK1rqgo4uAYAACCC0I01sdbqx784py9+97jm/EGVF2frQ286oB01hU6XBgAAkHAI3Vi18ak5/c03jurJYz2SpMaGKr3zN29Qdmaaw5UBAAAkJkI3VuX42QE99JUWDYzOKCvDo3fde4MONlQ7XRYAAEBCI3RjRYLBkP71v57XN/77tEJWuq6mUB98UwOLJAEAAFaA0I0r6hua0oMPH9Kpc8MyRnrtS3bo9b9ynTxul9OlAQAAJAVCNy7rsSM+/e03n9XUTEDF+Zn6/Tfs1w3XbHG6LAAAgKRC6MaypmcD+rvvHNNPnumUJN26t1zvee0+eXM43AYAAGC1CN1Y4kznsB58uEXdA5NK97j0tlft1ctuq2XfbQAAgDUidCMqFLL67qMv6F9+eFKBoFVthVcfelODasq9TpcGAACQ1AjdkCQNjc3o0189rKNnLkiSXnFHnX7rFXuUnuZ2uDIAAIDkR+iGnj7Rq7/+2hGNTc7Jm5Ou971un27ZXe50WQAAACmD0L2JzfmD+tL3W/X9x9slSTddu0UfeMN+FXkzHa4MAAAgtRC6N6nzvWN64OEWdfSMye0yesvLd+meu6+Ry8ViSQAAgPVG6N5krLV65KkO/cO/H9dcIKTKkhx98E0Nura60OnSAAAAUhahexMZm5zT575+RD8/3itJ+uWbq/WO/3GDsjJ4GwAAAMQTaWuTeO6FAT301RYNjs4oO9Ojd7/6Rt21r8rpsgAAADYFQneKCwRD+uqPTumbPz0ja6Wd2wr1wTcdUFlRttOlAQAAbBqE7hTWOzipBx9u0fPnh+Uy0mtfukOvf+l1crtdTpcGAACwqRC6U1RzS6f+77eOaXo2oJL8TN3/xgbt3V7idFkAAACbEqE7xUzN+PWFbx9TU4tPknTb9RV6z2tvUl52usOVAQAAbF6E7hRy+vywHny4RT2Dk0pPc+u+e/bqV160Tcaw9zYAAICTkiZ0G2PulfQ+STmSCiUNSfpra+2/OFpYAgiFrL7VdEZfeeSUgiGr+sp8ffBNDaouy3O6NAAAAChJQrcx5gOS3izpldZanzEmTdL/k/TLkjZ16B4cndanvnpYx14YkCS98q56vfXXdyvN43a4MgAAAMxL+NBtjKmV9ElJd1hrfZJkrfUbYz4oqdLJ2pz2i+M9+ut/O6rxqTnl56br/a/brwO7ypwuCwAAAIskfOhWeIR7xFr7zMKL1tpuSd3OlOSsWX9Q//S94/rBkx2SpH07tugDr9+vQm+ms4UBAABgWckQul8sqSMyp/v9krYoPJ/7H6y1/7TcE4wx90m6T5Jqamo2qMyNca5nTA88fEjnesflcRv9r1/frVfeuV0uF4slAQAAEpWx1jpdw2UZY45LqpX0nKT/Ialf0r2S/lXSx6y1n7jc8w8cOGAPHToU7zLjzlqrHzzRrn/8j1b5AyFt3ZKrD72pQdurCpwuDQAAAGGXHAVNhpHuTIV3LPmQtbY3cu0bxpjXSfoDY8ynrbVTzpUXf6MTs/rc14/qF63hP/5Lb6nRffdcr8yMZPjrAwAAQDKktvHIx6OLrh+R9JuSdktK/qHsS3j2zAV96quHNTQ2o5xMj979mpt0501bnS4LAAAAq5AMofuUpJskuRZdD0Y+Lr6eEgLBkB7+4Ul9u/kFWSvtqi3SB9/YoNKibKdLAwAAwColQ+j+D0mvk3SDpMcXXN8raVpSqxNFxVPPwKQeePiQznSOyGWk1//KdXrtS3bI7U7Jf18AAACkvGQI3f+m8K4l/8cY8wpr7YQx5k5Jr5b0Z9baSUerW2c/PdSpL3z7WU3PBrWlMEsffGODdtcVO10WAAAArkLCh25rbdAY8zJJfymp1RgzI2lW0u9Za//e2erWz9SMX5//1jE1H/ZJku64sVLvfs1Nys1Kc7gyAAAAXK2ED92SZK0dkvR2p+uIl1PnhvTgwy3qG5pSRrpb77jner3klhoZw97bAAAAqSApQneqCoasvvXTM/rKj04pFLLaXpWvD76xQVWleU6XBgAAgHVE6HbIwMi0PvXVw3ru7IAk6Z67t+stL9+lNI/b4coAAACw3gjdDnjquR597utHND7lV0Fehj7wuv3av7PU6bIAAAAQJ4TuDTQzF9A/fq9VjzzVIUlq2Fmq979uvwryMpwtDAAAAHFF6N4g7d2jeuDhFnX2jcvjdum3XrFbv3FnPYslAQAANgFCd5xZa/X9x9v1pe+3yh8Iqao0Vx9+8wHVVeY7XRoAAAA2CKE7jkYnZvWZrx3RoZN9kqRfvXWb3vaqvcpMp9sBAAA2E9JfnDx7+oIe+mqLhsdnlZuVpt977U26/YZKp8sCAACAAwjdcTA149cnvvy0pmcD2lNfrPvf0KAthVlOlwUAAACHELrjYHI6oOnZgAryMvSJ371dbheLJQEAADYzl9MFpDKPyxC4AQAAQOgGAAAA4o3QDQAAAMTZuoZuY8y/rOfrAQAAAKlg1QspjTHXSLpbUrkk96KHf3k9igIAAABSyapCtzHmXZI+J+lSqwPtVVcEAAAApJjVTi/5sKR3S9pirXUt/iXp2PqXCAAAACS31U4vuWCt/cJlHn/z1RQDAAAApKLVjnR/xxhz/WUef8fVFAMAAACkolWNdFtr/9wY85fGmO2SXpA0tajJvZLes17FAQAAAKlgLQspP3SZJiykBAAAABZZ7fSSP1B4CkkJCykBAACAlVntQsp+a+3fX+ZxFlICAAAAi6x2pPu/WEgJAAAArM5qR7pnJX3XGHNELKQEAAAAVmS1ofsPIx/rLvE4CykBAACARVY7veTZ5RZQspASAAAAuLTVhu4/vsLjTC0BAAAAFllV6LbW/scVHn/86soBAAAAUs9qR7oBAAAArBKhGwAAAIgzQjcAAAAQZ4RuAAAAIM6uGLqNMbcbY8xGFAMAAACkopWMdP+XpAvGmIeNMa8zxhTEuSYAAAAgpawkdBdL+m1Jk5IekNRvjHnUGPMhY8yuuFYHAAAApIArhm5r7bS19nvW2ndYa6slvUjh0e97JT1njGkzxnzWGPMrxpj0eBcMAAAAJJtVL6S01h6x1v4fa+2tkiok/Vnk49clDRpjvmOM+e11rhMAAABIWle1e4m19oK19svW2tdIKpF0j6QOSR+++tIAAACA1LBuWwZaawPW2p9Yaz9grd25Xq8LAAAAJDv26QYAAADijNANAAAAxBmhGwAAAIizVYVuY8zfGmPuilcxAAAAQCpa7Uj370r6kTHmVxY/YIzJNcZcsz5lAQAAAKljLdNLHpH0HWNM46LreyQ9f/UlAQAAAKllLaH7LyQ9JOk/jDF3rnM9AAAAQMrxrOVJ1to/jhz5/p/GmF+11j61znUBAAAAKWNNoVuSrLUfMcakSfqhMeal61gTAAAAkFLWHLolyVp7f2TE+0eSPro+JQEAAACpZbWh+w5JZxZesNa+xxjjkfTZdasKAAAASCGrCt3W2icvcf13jTF+SW9Zl6oAAACAFLJuJ1Jaa99rrS1Yr9dLZmmecLcOjc/qu4++oFDIOlwRAAAAnMQx8HFQkJehexuvUShk9Y/fa9Wf/sPPNTw243RZAAAAcAihO07e+oo9+sPfukV52ek6/Hy/3vtQsw6d7HO6LAAAADiA0B1HL9pboc998KBuuKZEIxOz+tN/+Ln+/t+fkz8QdLo0AAAAbCBCd5wV52fp4+94sf7Xr++W22X0vcfadP9fP6bOvnGnSwMAAMAGIXRvAJfL6NW/dK3+6j13qqI4R+3dY3r/px/Vj37eIWtZZAkAAJDqCN0baEdNoT7z+3frlw5Ua84f1N9841l98p+f0fjUnNOlAQAAII4I3RssOzNNH3j9ft3/xgZlZXj05LEevffBJh0/O+B0aQAAAIgTQrdDDu6v0mfvP6jrago1MDqjj37+CT38yEkFgyGnSwMAAMA6I3Q7qLw4R5/8vTv02pfskJX0b/91Wh/528fVNzTldGkAAABYR4Ruh3ncLr3513bpE++8XcX5mTp1bljvfahJjx3xOV0aAAAA1gmhO0Fcf02JPnt/o27dW66pmYAeeLhFn/naYU3N+J0uDQAAAFeJ0J1AvDnp+oO33qJ3vfpGpae59ZNnOvX+Tz+qM53DTpcGAACAq5CUodsY8zNjjDXG1Dpdy3ozxujXbqvVp99/l2orvOoZmNSHP/czfbvpjEIh9vQGAABIRkkXuo0x90q6w+k64q2m3KuH3neXfuPOegWCVl/6/gl97ItPaWhsxunSAAAAsEpJFbqNMemS/kLSD5yuZSOkp7l13z3X649/50Xy5qTr6JkLes+DTXr6RK/TpQEAAGAVkip0S3q3pEOSnnG6kI108+5yfe6DjbppxxaNTc7p4//4C/3dd45pzh90ujQAAACsQNKEbmNMkaQPSfoDp2txQpE3U3/69tv0W6/YI4/b6PuPt+v+v35M53vHnC4NAAAAV5A0oVvSH0t62FrbcaWGxpj7jDGHjDGHLly4EP/KNojLZfSbjdfor95zpypLctTRM6YPfPpR/fDJdlnLIksAAIBElRSh2xhzjaTXSvrEStpba79orT1grT2wZcuW+BbngGurC/WZ3z+ol9xco7lASP/3W8f0519+WmOTc06XBgAAgGUkReiW9FeSPmmtHXW6kESRleHR+163Tx96U4OyMz36+fFevfehJh17IXVG9gEAAFJFwoduY8ydkvZK+rzTtSSiu/ZV6bP3N2rntkINjs7oD7/wpP75BycUCIacLg0AAAARCR+6Jb1UklvSM8aYo8aYo5LeGXnsB5FrL3esugRQVpStT777Dv3Pl+6QkfSNn5zRR/7mcfUOTjpdGgAAACSZZFyAZ4z5E0kfk1R3pYWVBw4csIcOHdqIshLC8bMDeuirhzUwMq2sDI/ede8NOthQ7XRZAAAAm4G51APJMNKNVdi7vUSfu/+gXnxDhaZnA3roq4f1qa+2aGrG73RpAAAAm1ZShW5jzMuXm17iXEWJKTc7XR95y836vdfcpPQ0t5pafHrfp5p1+vyw06UBAABsSkk5vWQ1Ntv0ksU6+8b14MMtauseldtl9MaX7dS9jdfK5brk/34AAABgbZhesllVl+XpwffdqVfeVa9gyOqff3BSf/R3T2pwdNrp0gAAADYNQvcmkOZx6+2vul4fe9utKsjN0LEXBvSeB5v18+M9TpcGAACwKRC6N5EDu8r02Q8e1P7rSjU+NadPfOlpff5bz2rWH3S6NAAAgJRG6N5kCvMy9bG33arfeeUeedxGP3iyQ7//mUd1rmfM6dIAAABSFqF7E3K5jO65+xo9+N67tHVLrs73jusDn3lU//l4m1J9YS0AAIATCN2b2PaqAn3mA3frpbfUyB8I6QvfeU7/55+e1ujErNOlAQAApBRC9yaXmeHRe//nPv3vtxxQTqZHT5/o1XsfatKzpy84XRoAAEDKIHRDknTHjVv12fsbtau2SENjs/qjLz6pL3+/VYFgyOnSAAAAkh6hG1GlRdn6i3fdrjf8ynUykr7V9II+/LmfqXtgwunSAAAAkhqhGzHcbpde/6s79RfvvkNbCrN0pnNE7/9Us3566DyLLAEAANaI0I1l7a4r1mfvb9QdN1ZqejaoT//rET30lcOanPY7XRoAAEDSIXTjknKz0vThNx/Qe197kzLS3Xr0iE/v+1SzTp0bcro0AACApELoxmUZY/TSF23TZz5wt7ZX5atvaEr/+28e17/99/MKhphuAgAAsBKEbqxIVWmeHnjPnbrn7u0Khawe/uEp/dEXntTAyLTTpQEAACQ8QjdWLM3j1u+8cq/+9O23qSAvQ8+dHdB7HmzSU891O10aAABAQiN0Y9X27yzV5+5vVMPOUk1M+/XnX35Gf/vNZzUzF3C6NAAAgIRE6MaaFORl6GNvu1Vvf9VeedwuPfJUh37/M4/q2TMXFGKuNwAAQAyT6nsvHzhwwB46dMjpMlJae/eo/upfDsnXHz5Ep6QgSwf3V6mxoUo15V6HqwMAANgw5pIPELqxHmbmAvrWT1/QTw6d14Xhi4sr67fmq7GhSnftq1KRN9PBCgEAAOKO0I2NEQpZnWgfVFOLT08826XJmfA8b5eRbrx2iw42VOu26yuUleFxuFIAAIB1R+jGxpvzB/XMyT41HepUy6k+BYLh91pGulu37a1QY0O1bry2RG43SwsAAEBKIHTDWWOTc3ri2S41tfh0suPiiZYFeRm6a99WNTZUa/vWfBlzyfcqAABAoiN0I3H0DEzq0SM+NR3qVPfAZPR6dVmuDu6v1sH9VSotynawQgAAgDUhdCPxWGt1pnNETS2deuxIl8Ym56KP7akvVmNDlW6/catys9IcrBIAAGDFCN1IbIFgSEee71dTi0+/ON6juUBIkuRxu3TLnjI1NlSrYWeZ0jzM/wYAAAmL0I3kMTXj15PHetTU0qnnzg5o/i2al52mO24Mz//eWVvI/G8AAJBoCN1ITgMj03rsiE9NLT519IxFr5cXZ+vu/VVqbKjW1i25DlYIAAAQRehG8mvvHlVzi0/Nh30aGpuJXt9RU6DGhmrdedNW5edmOFghAADY5AjdSB3BkNXxFwb005ZOPfVct6Zng5Ikl8to/3Wl+qWGat2yt1wZaW6HKwUAAJsMoRupaWYuoKdbe9XU4tPh5/sVCoXfz1kZHr34hvABPHu3l8jtYv43AACIO0I3Ut/I+Kx+drRLTS2dOtM5Er1enJ+pu/dVqfFAtWorvM4VCAAAUh2hG5uLr39czS0+NR32qX9oKnq9tsKrxoZq3b1/q4rzsxysEAAApCBCNzYna61OdgypqcWnx492aWLaL0kyRrrhmhI1NlTrtusrlJ3JATwAAOCqEboBfyCoQyf71dTSqWdO9CkQDB/Ak57m1q17y9XYUK2bdmyRx80BPAAAYE0I3cBCE1NzevzZbjUf9qm1bTB6PT83XXftq9LB/VW6trqAA3gAAMBqELqBS+kbmlLz4U41HfKp68JE9PrWLblqbKjS3furVF6c42CFAAAgSRC6gSux1uqsb1RNLZ167EiXRiZmo4/trivSwYZq3XFjpfKy0x2sEgAAJDBCN7AawWBIR89cUNMhn5463qM5f/gAHo/bpZt3l+ng/irdvLtMaR4O4AEAAFGEbmCtpmb8+vnxXjW1dOrYmQuKnL+jnKw03XFjpRobqrWrtkguDuABAGCzI3QD62FwdFqPHelSc4tPbd2j0eulRdk6uD+8ALO6LM/BCgEAgIMI3cB6O9czpqaWTj162KeB0Zno9WuqC9S4v0p37tuqwrxMBysEAAAbjNANxEsoZNXaNqimlk49caxbUzMBSZLLZbRvxxYdbKjWrXvLlZnucbhSAAAQZ4RuYCPM+oN65kSvmg751HKqT8HIBPCsDLduu75SB/dX6YZrt8jN/G8AAFIRoRvYaKMTs3r8aJeaDvv0/Lnh6PUib4bu2lelxoZq1VV6OYAHAIDUQegGnNQ9MKHmFp+aW3zqGZyMXt9WnqeDDdW6e1+VthRmOVghAABYB4RuIBFYa/X8+WE1HerUz452a3xqTpJkjHT99hId3F+lF99QqZysNIcrBQAAa0DoBhKNPxDS4VN9ajrs09OtvfIHQpKkdI9Lt+wpV2NDtfbvLJXH7XK4UgAAsEKEbiCRTU779cSxbjW3+PTc2YHo9bzsdN21b6sONlTpuppC5n8DAJDYCN1AsugfntKjh31qavGps288er2iJEeN+6t0d0OVKktyHawQAABcAqEbSDbWWrV3XzyAZ3h8NvrYzm2FOthQrTtv2ipvTrqDVQIAgAUI3UAyC4asnj1zQc0tnXrquR7NzAUlSW6X0YFdZTrYUKVbdpcrPc3tcKUAAGxqhG4gVczMBvTz4z1qOuzT0ef7FTl/R9mZHt1+Q6UaG6q1p75YLg7gAQBgoxG6gVQ0PDajx452qbmlUy/4RqPXSwqydHB/lRobqlRT7nWwQgAANhVCN5DqOvvGo/O/+4eno9frt+arsaFKd+2rUpE308EKAQBIeYRuYLMIhaxOdgypqaVTjx/t0uRMQJLkMtKN127RwYZq3XZ9hbIyPA5XCgBAyiF0A5vRnD+oZ072qbmlU4dO9ikQDH+/Z6S7ddveCjU2VOvGa0vk5gAeAADWA6Eb2OzGp+b0+NEuNbX4dLJjKHq9IC9Dd+3bqsaGam3fms8BPAAArB2hG8BFvYOTaj7sU3NLp7ouTEavV5flqrGhWnfvq1JpUbaDFQIAkJQI3QCWstbqTOeImlo69bOjXRqdmIs+tqe+WI0N1br9xkrlZqU5WCUAAEmD0A3g8gLBkI4836/mFp9+frxHc4GQJCnN49LNu8vU2FCthp1lSvMw/xsAgEsgdANYuakZv5481qPmw5069sKA5m8Tedlpuv3Grdp/Xal21xUpPzfD2UIBAEgshG4AazMwMq3HjvjU1OJTR89YzGPVZXnaW1+s3fXF2ltfrJKCLIeqBAAgIRC6AVy99u5RPXmsRyfaB3WqYyg6BWVeaVG29tYXa0/kV2VJDruhAAA2E0I3gPXlD4T0QueIWtsH1do2qJPtg9GDeOYV5GVoT93FEL6twiu3ixAOAEhZhG4A8RUMWZ3rGdPxtgGdaBtSa9ugRiZmY9rkZHq0KxLC99YXa3tVAQszAQCphNANYGNZa9U9MKnjZwfV2jag1vYh9Q9NxbRJT3Nr57bC8Eh4XbGu21aoTI6nBwAkr+QN3caYmyS9W9J+SR5JaZL+W9LHrbUXrvR8QjeQOC4MT0eno7S2DaqzbzzmcbfL6Jqqguh0lN11RcrNTneoWgAAVi2pQ/cpSa2S3mKtnTTGbJX0E0kuSTdaa6cv93xCN5C4RidmdaJ9UK1tQ2ptG1Bb16hCC25Jxkjbyr3REL6nvlhF3kznCgYA4PKSPnS/wlr7woJrvyPpHyS92lr7rcs9n9ANJI+pGb9OdQyH54W3D+n5c8MKBGN3SKkoyYkuzty7vVhlRdnskAIASBSX/IGUDJMnb7DWzi261h35WLjRxQCIn+zMNO3fWar9O0slSXP+oM50jkQXZ57sGFTPwKR6Bib138+clyQVeTNj9gqvLsuTix1SAAAJJuFHupdjjHmfpE9L2mmtPb3M4/dJuk+SampqGs6dO7fBFQKIh2AwpPbuMR1viyzObBvS+FTsv8nzstO0e8E2hfVb8+Vxs0MKAGBDJO/0ksWMMW5JhyU9ba19+5XaM70ESF2hkJWvf1yt7UNqjeySMjA6E9MmM92tnduKtGd7OITvqClURprboYoBACkupUL3n0j6DUl3W2snrtSe0A1sHtZa9Q9PR0fBW9sG1HVhMqaNx210bXVhdCR8V22RcrLSHKoYAJBiUiN0G2N+S9KHJB201vav5DmEbmBzGx6f0Ym2oei88PaeUS287bmMVFuZH50XvqeuWAV5Gc4VDABIZskfuo0xb5b0vyW9xFrbu9LnEboBLDQx7depjiEdPxveIeVM57ACwdj74NYtudq7vVi768KLM0uLsh2qFgCQZJI7dBtj3iTpI1oQuI0xr5BUaa394uWeS+gGcDkzcwGdPj8cnY5y6tywZueCMW22FGZFtyncU1+sqtJctikEACwneUO3MeaNCu/J/UeSFo5w3ympx1r7J5d7PqEbwGoEgiGd9Y1ETs0cUmv7oCan/TFt8nPTY3ZIqavMl5ttCgEASR66h3Tp/bj/lNANIJ5CIavzfeNqPTsQ3iWlbUBDY7MxbbIyPNpVVxQdDd9RU6A0DzukAMAmlLyh+2oRugGsJ2utegYndaJtUMfbBnWibUg9g7E7pKR5XNpRUxhdnLlzW6GyM9khBQA2AUI3AMTL4Oh0dIeU1rZBnesdj3nc5TKq3xreIWVPfXiBpjcn3aFqAQBxROgGgI0yPjWnE22D0ekoL/hGFQrF3mtryvNiFmeWFGQ5VC0AYB0RugHAKdOzAT1/big6HeX5c0OaC4Ri2pQVZUcD+N76YlWU5LBDCgAkH0I3ACQKfyCoFzpHo9NRTnYMaWomENOmIC8jGsD31BdrW7lXLnZIAYBER+gGgEQVDFl1dI+qtX1QrZHR8JGJ2B1ScrLStKu2KBrCt1cVKM3jcqhiAMAlELoBIFlYa9V1YSKyV3j4V//wdEyb9DS3dm4rjE5JuW5boTLTPQ5VDACIIHQDQDLrH56KWZzZ2TcR87jbZXRNdUF4ceb2Yu2uLVJuNjukAMAGI3QDQCoZnZjVifb5vcIH1dY1qoUbpBgjbSv3RvcK31NfrCJvpnMFA8DmQOgGgFQ2NePXqY7h6OLM0+dHFAjG7pBSUZITnRO+p75YZUXZ7JACAOuL0A0Am8mcP6jT54fDizPPhndImZkLxrQpzs+MTkfZU1es6rI8dkgBgKtD6AaAzSwYDKmte3TB4swhjU/NxbTJy07T7gUH9mzfmi+3mx1SAGAVCN0AgItCIavO/nGdaAvPC29tG9Tg6ExMm8x0t3bWFkVD+I6aQmWkuR2qGACSAqEbAHBp1lr1DU2FF2eeHdSJ9kF1XZiMaeNxu3RtdYH2bi/W7rpi7aotUk5WmkMVA0BCInQDAFZneGwmemBPa9ugOnrGtPBHhstIdVvzw/PC68NBvCAvw7mCAcB5hG4AwNWZmPbr5IIQfqZzRMFQ7M+QqtLc6HSUPfXFKi3MdqhaAHAEoRsAsL5m5gLhHVLODqq1fVAnO4Y154/dIWVLYVY4gEdGw6tKc9mmEEAqI3QDAOLLHwjpbNdIdHHmifYhTU77Y9rk56Zrd11x9NCeusp8udmmEEDqIHQDADZWKGR1rncsZoeU4fHZmDbZmR7trC2KHtpzbXWB0jzskAIgaRG6AQDOstaqZ3AyOh2ltW1QvYNTMW3SPS7t2FYYnY6ys7ZIWRkehyoGgFUjdAMAEs/g6LRa56ejtA3qXO94zOMul9H2rfnRhZm764rlzUl3qFoAuCJCNwAg8Y1Nzulk+8XpKGe7RhVatENKTXledHHm3u3FKs7PcqhaAFiC0A0ASD7TswGd6hiKTkc5fW5Yc4FQTJvy4uzo4sw99cWqKMlhhxQATiF0AwCSnz8Q1JnOkehe4SfahzQ9G4hpU5iXEbNX+LZyr1zskAJgYxC6AQCpJxiy6ugeDYfwyGj46MRcTJucrDTtrisKL87cXqxrqgrkcbscqhhAiiN0AwBSn7VWvv4JnWi/uDizf3g6pk1GulvX1RRG9wq/bluhMtPZIQXAuiB0AwA2p/7hqZi9wn39EzGPe9xG26sKtKu2SHWV+aqr9KqqNJf9wgGsBaEbAABJGhmf1Yn2i9NR2rtGtWiDFLldRlWluaqtCIfw2kqvaiu8KvJmskgTwOUQugEAWM7UjF8nO4Z0pnNEHd1j6ugZVffApJb78ZiXnR4N4XUVXtVW5Ku6PE8ZaYyKA5BE6CZ0AwBWbmYuoPO942qPhPCOnjF1dI9pYtq/pK3LSJVbclVb4VVdZX50VHxLQRaj4sDmQ+gGAOBqWGs1ODqjjp4xtXePqqN7TO09Y+q6MLHkAB9Jysn0qLYyPzwiHgni28q9yuRYeyCVEboBAIiHOX9QnX3zo+LhkfH27jGNTc4taWuMVF6cc3FUvMKrukqvSguz2UscSA2EbgAANoq1ViPjs2rvGVNH92jk45h8/eMKBJf+3M3KcGtbeez0lNoKr7Iz0xyoHsBVIHQDAOA0fyCkrgsT0ekp8yPjQ2Ozy7YvLcqOTk+pqwgH8vLiHLkZFQcSFaEbAIBENToxG50jPr9w83zvuPyB0JK26WlubSvPi05PmR8Zz8tOd6ByAIsQugEASCbBYEjdA5PhUfGeseic8YGR6WXbl+RnhhduLpiesnVLrtwceQ9sJEI3AACpYGJqLjpHvCM6Mj6uOX9wSds0j0vVZXmRIH5xJ5X83AwHKgc2BUI3AACpKhiy6h2cjExRuThfvG9oatn2hXkZS6anVJXmKc3DqDhwlQjdAABsNpPTfp3rHYse7jM/Mj49u3RU3O0yqi7Lizlts7bSq8K8DA75AVaO0A0AAKRQyKp/eCpmX/GO7jH1DE5quUjgzUmPTk+ZHxmvKctTepp744sHEh+hGwAAXNr0bEDnF4yKz+8xPjkTWNLW5TLauiU35rTN2op8lRRkMiqOzY7QDQAAVsdaqwsj00ump3T1Tyi0THzIzUqLCeF1kVHxzAzPxhcPOIPQDQAA1sesP6jO3vHwkffzI+PdYxqfmlvS1hipsiQnOkd8fjvDsqJsRsWRigjdAAAgfqy1GhqbWTQqPqbOvnEFlxkWz8rwROeIzy/c3FaRp+zMNAeqB9YNoRsAAGw8fyAoX//ExYWb3eHR8ZHx2WXblxdnR6en1FZ6VVfpVXlRjlwuRsWRFAjdAAAgcQyPz+hcZDS8vTs8On6+b1yBYGhJ24x0t2rLvTHTU2or85Wbxag4Eg6hGwAAJLZAMKSuCxMx01Pau0c1ODqzbPsthVnREF4XGRmvLMmR280hP3AMoRsAACSnsck5neu5eNpme8+YzveMaS6wdFQ83eNSTXnekoWb+bkZDlSOTYjQDQAAUkcwZNUzMBGzcLO9e1T9w9PLti/yZi5YtOlVXWW+tpbmysOoONYXoRsAAKS+yWl/zNSUjp4xnesZ08xccElbj9uouixvycLNwrxMBypHiiB0AwCAzSkUsuobmgrvK9598dTNnsHJZdsX5GZEtzOcHxWvLstVmse9wZUjCRG6AQAAFpqeDehc7/zhPqPREfKpmcCSti6XUVVp7sWFm5XhEzeLvJkc8oOFCN0AAABXYq3VheHp6PSU+RM3ewYmtMwZP8rLTotZtFlX6VV1WZ4y0z0bXzwSAaEbAABgrWbmAursG4/unjI/Oj4x7V/S1mWkipLc6MLNusp81VZ4taUwi1Hx1EfoBgAAWE/WWg2Nzah90fQUX/+EQssMi2dnemIO96mr9GpbuVdZGYyKpxBCNwAAwEbwB4Lq7Ju4uHCzO7zH+OjE3LLtK4pzYqan1Fbkq6woWy4Xo+JJiNANAADgpOHxmci0lDF19IRHxjv7xhUILs1imelubVuwaHN+hDwnK82ByrEKhG4AAIBEEwiG1NU/EZ2eMj9ffGhsZtn2pYVZMXuK11Z4VVGSKzej4omC0A0AAJAsRidmda53LDo9paNnVOd6x+UPhJa0Tfe4VFMROW2z0qu6inxtq/DKm5PuQOWbHqEbAAAgmQWDIXUPTEbniM8v3LwwPL1s++L8zJh9xWsrvdq6JVcet2uDK99UCN0AAACpaGLar3MLjr3v6B5TR++YZueCS9p63C7VlOXFTE+prchXQV6GA5WnJEI3AADAZhEKWfUOTcZMT+noGVPv4NSy7QvyMpactllVmqs0j3uDK096hG4AAIDNbmrGr3M94+HtDOdHxXvGND0bWNLW7TKqKs1VbUU4hM9va1jkzeSQn0sjdAMAAGApa636hqaic8TnR8a7Bya1XEzMy06PhvC6yPSU6vI8ZaQxKi5CN6EbAABgNWZmAzrfNx6zr3h795gmp/1L2rqMVLklN2bRZm2FV1sKsjbbqDihGwAAAFfHWquBkZloCA/vpDKmrgsTCoWWZsqcrLRwEK+4OD1lW7lXmRkeB6rfEIRuAAAAxMecP6jzfePROeIdPaNq7x7T2OTckrbGSOXFOdEj78Oj416VFmbLlfyH/CR36DbGlEr6tKQDkUvPSXq/tdZ3pecSugEAADaetVbD47PROeLzCzd9/eMKBJfmz6wMT3QHldrodoZeZWemOVD9miVv6DbGpEt6RtJpSa+XZCX9k6QXS9pnrZ243PMJ3QAAAInDHwjJ1z++YNFmeI/x4fHZZduXFWVHg3hdRXi+eHlxjtyJOSqe1KH77ZK+KGm7tbYtcq1cUpekj1hrH7jc8wndAAAAiW9kfDZ8yE/PxYWb53vH5Q+ElrTNSHdrW3ledHrK/E4qudnpDlQeI6lD9yOSdllrty26/pykSWvtrZd7PqEbAAAgOQWDIXVdmIhuZxg+7GdUA6Mzy7YvKchSbYVXt+wp16/dVruxxYZdMnQnw9LRGxSeWrJYu6RfXu4Jxpj7JN0nSTU1NfGrDAAAAHHjdrtUU+5VTblXd+27eH18ai5mekp4ZHxcAyPTGhiZVnF+pnNFX0IyhO4SSS3LXB+TlG2MybLWTi98wFr7RYWnpOjAgQOJPZQPAACAVcnLTtf120t0/faS6LVgyKp3cFId3WMqKSB0r6eEnD0PAACAjed2GW3dkqutW3KdLmVZLqcLWIEBSXnLXM+TNLV4lBsAAABINMkQuo9Jql3mep3C+3UDAAAACS0ZQve3JW0zxtTOXzDGlEnaJelbThUFAAAArFQyhO4vKzyi/ZfGGI8xxiXpkwrvXvJ5JwsDAAAAViLhQ7e1dk7SSyUFJZ2QdFKSV9IvXek0SgAAACARJMXuJdbaPklvcLoOAAAAYC0SfqQbAAAASHaEbgAAACDOCN0AAABAnBG6AQAAgDgjdAMAAABxRugGAAAA4ozQDQAAAMQZoRsAAACIM0I3AAAAEGeEbgAAACDOCN0AAABAnBlrrdM1xJUx5oKkc1doViJpYAPKSTX02+rRZ2tDv60N/bY29Nva0G9rQ7+tTaL224C19mXLPZDyoXsljDGHrLUHnK4j2dBvq0efrQ39tjb029rQb2tDv60N/bY2ydhvTC8BAAAA4ozQDQAAAMQZoTvsi04XkKTot9Wjz9aGflsb+m1t6Le1od/Whn5bm6TrN+Z0AwAAAHHGSDcAAAAQZ4RuAAAAIM4I3QAAXIIx5uXGmF5jzJedriWZ0G+rR5+tTTL1W0qFbmPM+40xJ4wxx4wxh40x96zwea8yxnzdGPNs5NcZY8yXjTHVy7RtMMY8aow5box53hjzoDEmc5l2240x3zfGnDTGnDbG/KMxpnAd/pjrbi39ZoxJM8a83hjzY2NMizGm1RhzyBjzv4wxZlHbDmPM0WV+DRpjfrag3UFjzNAy7Z6Mwx/7qm1Av701ciNZ3B/fXOZ1eb+F2xljzK8bY/498trPRb6nP2CMSV/0mrzflrbn/nbxeVnGmL+V9AlJZZdpx/0t9nkr7beUu7/Fs8+4ty37vBW91yJtE+PeZq1NiV+SPqLwyUTbI5+/VJJf0q+t4LmHJH1ZUkbk8zJJxyWdkuRZ0O5aSWOS3hf5vEDSc5L+ddHrFUvySfqUJCMpQ9J/S3pCksvpvlqPfpN0q6SQpNcvuPZqSVbSxxe17bjEa5yS9LsLPj8o6ctO90kC9dtbJf3JCmrh/XbxWnnk2od1caH47ZJmJP3Lotfk/Rbblvtb7HPfLOkBSemR/lr2vcL9bc39llL3t3j3Gfe2q3qvJcy9zfEOX6e/tAJJk5L+bNH1/5TUuoLn/52kqkXX3h/5S7xhwbWvSOqYf8NHrr0m0u7mBdc+Eaknd8G1myPtXuN0f61Hvyn8w/zxZa7/LPLmXthH9yzT7k5JE5K8C64lxY1iA/vtrVrZDyXebxd/CJVL8i3T7l8UDp9beL9d8v3G/S223cIBl8v9QL9nmWub8v62yn5LmfvbRvQZ97areq8lzL0tVaaXvExStqSmRdd/Kmm3MWbn5Z5srX2Htda36LJX4Y4elCRjjEfSqyQ9aiN/Ewu+hiTdu+DavZIOWWsnFlw7pPAPuYXtnHY1/fYLSY3LXO+WlCMpbf6Ctfa7y7R7m6R/s9aOrabgBLEh/bYKvN8i/Wat7ZVUf4l2RlL+KmtOBHHvN+5vS1lrAyv5ItzfYq2031YhGd5vce8z7m1LraTfEu3eliqh+4bIx/ZF19sXPX5FkXlTjZLeI+lj1tquyEP1Cv+Qivka1tpBSePzXyMyR+jaZdpZSedWU8sGWHO/2TD/Mg/tkPSUtXbuUs81xuQr/K/Mv1vu+caY70Tmqp0yxvy9Mab2kn8CZ2xkv91ijHkkMn/vOWPMp40xJfMP8n5b2m+XeO/tkHReUtvi67zfJHF/Wzeb+f62Bqlyf9uQPuPetiYJdW9LldA9/006vuj6/ChD8UpexBjz+5IuKPxfG5+S9Ocr+BrzX2f+axQp3K9XapcI1qXf5hljblH4jfkHV2j6RkmnrbVPL7ruV/i/yT5qrb1R0t2SaiQ9k2A3i43qtxlJAUnvsNZeL+keSXdJesoYUxBpw/vtCu83E14Q/WuS/sBaG1rwEO+3K3+N+a+z6e9vq8D9bWVS6f7myHuNe9tVfY35r7Oh97aEDN3GmJcYY+wKfjVf6aVW83WttZ+y1pZIepGk10lqWm5161V8nVXVs1pO9Vvka+dK+kdJf2itfewKzd+mZY5vtdY+Ya29w1p7IvJ5n6TfUnje10dXW9NKJWq/WWu/Zq19pbX2XOTzs5LeKekaSe9eycuvtp7VSNR+W9TOI+lLkv6ftfYrCx/j/XbF79PVfp2Ufb+tEve3FUjk+1ui9tmiGrm3Xb0Nv7d51uuF1tmTknatoN1U5ONA5GOeInOwF3yuRdeuyFr7nDHmnZE63qXwqPfCr7HYwq87pPBc8Cu1iwdH+s2Etyv6lqQfW2v/4gpt90m6TtLDK3lta223MaZT4QVh8ZLw/bZAi8KjGPP9wfvt0u2MwuFnSOHv4yvaxO837m/rgPvbVUuU+1tC9xn3tlVJqHtbQoZua+2UwtstrdSxyMdahVeozqtb9PgSkTdvprV2etFDz0Y+NkQ+tim8qrV20fOLFf4LORapfcYYc3qZdkbSNkk/vMKfZc02st/mRX6Qf1vSCWvt/Sv4mm+X9LXlFhgZY4okjVprg4seCkpyr+C11yRR+80Ys8Vae2GZh6wi/cH77ZLtjKQvKPy9+fpl3lO832Jxf1sf3N9WKJHvb4naZxL3tjVIqHtbQk4vWYNHFP6X08FF1xsV/iETfRMYY7JNeKHLvG2Sji7zmrWRj4NSdJXs9yTdHflLWPg1pPAPNC34/c3GmJwF1xoU3hHlW1f+42yYq+m3hSNnZ6y1H1hw/e+MMRWLv5gxJkvSG7TMf71GfFvSLYueU6Lw39HhFf6ZNsJG9dszy/TjXoX3JF3YH7zfFvRT5Pvz8wrP0Xt95HtXxpiPGWMaFrwk77dIv3F/W9pvq8X9bdVS6f62IX3GvW31/ZZw9zabAHs1rscvhTdYvyCpPvL5S7TMBuuSno+0y4l8Xqvwv6x/Txf3+i2W9F8Kvxn2Lnju/Abr74l8nq/wv5KW22C9S9KDurjB+o+VYJv5X2W/pSv8Rm6V9KZFv05Kql3ma71F0rOXqaVZ4X9N5kc+z5L0DYX/lbr3av6cydhvCv/L/0sK/0/M/PvqMUl9kip4vy3tt8if//OR/njronaPSTrI++2S7zfub5F+W+Z1rK6w77G4v62q35Ri97d495m4t13Ney1h7m2Od/g6/+W9X9KJSGce0fKHFjQr/ANo/vTJdEn3SfqJwqdQHlP4v0e+ImnPMs8/IOnRyGuclvTQ/E1jUbtrJH1f4R9spxVevFTodB+tY7+9MvImv9Sv2mVe4zFJ775MHbcrfBN+LlJLl6R/l3ST033kRL8pvCr965H+OC6pM/K+rOP9tny/Kbwrx+XaHeT9dunvU3F/y1hwLUvh/wU9GumrocjvP3mJr7Pp72+r6Tel4P0tnn0m7m1X9T2qBLm3zY/sAgAAAIiTVJnTDQAAACQsQjcAAAAQZ4RuAAAAIM4I3QAAAECcEboBAACAOCN0AwAAAHFG6AYAAADijNANAAAAxBmhGwAAAIgzQjcAQMaYbGPMtDHmTqdrAYBUROgGAEjSSyRNS3rS6UIAIBURugEAkvQKSY9Ya4NOFwIAqYjQDQBJwhiTboyZM8bcv+h6tTHGGmPuvYqXf7mk76+ghlcYY54xxpw2xjxnjPmxMebNxphXGmOORur4uDHmr4wxR4wxXcaY34s896PGmEPGmHPGmHdfRa0AkHQI3QCQPPZKSpN0ZNH1myIfj67lRY0x+ySVS3rkCu1eLem7kj5hrd0h6UZJz0n6nLX2e9ba+TreJOkr1tp9kv5U0ueMMQ9I+r619oCkT0au7VhLvQCQjAjdAJA89kU+Hl3m+riktjW+7iskPWWtHbpUA2OMkfSQpGZr7XclyVobkvQxSROLmh+x1j4b+f23Ih+3L7j2DUlG0t1rrBcAkg6hGwCSxz5J55YJx/skPWuttWt83VfoylNLrpNUI+mZhRettRPW2qpFbV9Y8PuhZa4NRj6Wr7JOAEhahG4ASB77tXRqiRQO3UfX8oLGmFJJB3Tl0F0S+XjJ0fAFpuZ/s+AfAstdc6+wTABIeh6nCwAAXFlkesf1kn606Pp2Sdu0IHQbY/5Z0g2RTwsUDsz5l9iZ5OWSOq21rVcoYSDysXC1tQMAGOkGgGSRJylXUs+i6x+KfDw6f8Fa+5bIosZ7JI1KuvcyWwGuZGqJJD0v6bzCo+JRxpgiY8zPjTEFK3gNANi0CN0AkASstWMKB+7XGmMqjTH1xphPSXqNpICk4wvbG2N2KRym32Wt/dGSFwy3SZP0Uq0gdEemhNwvqdEY8xuR53sk/bmks9bakbX+2QBgMyB0A0DyeIukKoV3Kfle5ON3JJ2y1s7ONzLGNEj6pqQ3WWufuMzr3aXwNMPmlXxxa+03Jf2mpD82xpxReH75nKS3G2PuNMYcjTR9pzHmb40x16/g2ldW8rUBINmZtS92BwAkGmPM3ZL+RuEpJaev0PbTkuqtta/akOIAYBMjdANACjHGXJA0o4vb8gUl3Wqt9S/T9rSkB6y1f7+BJQLApkToBgAAAOKMOd0AAABAnBG6AQAAgDgjdAMAAABxRugGAAAA4ozQDQAAAMQZoRsAAACIM0I3AAAAEGeEbgAAACDO/n84khLpuAiAoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(u*100,z)\n",
    "plt.xlabel('$u_z$ / cm')\n",
    "plt.ylabel('$z$ / m');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That doesn't look plausible. What's missing are the boundary conditions: \n",
    "\n",
    "$$\n",
    "    u (z = 0) = 0 \\qquad u'(z=H) = 0\n",
    "$$\n",
    "\n",
    "They can be integrated easily by manipulating the first and last equations (rows).\n",
    "\n",
    "Let's start with the lower boundary condition, $u(z=0) = u_0 = 0$. Thus, we want the first equation to yield $u_0 = 0$ which can be achieved by manipulating the RHS and System matrix such that the first equation reads\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "1 & 0 & \\dots & 0 \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "u_0\n",
    "\\\\\n",
    "\\vdots\n",
    "\\\\\n",
    "u_9\n",
    "\\end{pmatrix}\n",
    "= \n",
    "0\n",
    "$$\n",
    "\n",
    "Nota bene: if a non-zero displacement were given as the boundary condition, the RHS entry would have to be set to that given value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lower bc\n",
    "A[0,:] = 0.\n",
    "A[0,0] = 1.\n",
    "RHS[0] = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we get to the upper boundary condition, $u'(z=H) = 0$. In other words, we want the slope of the $u_z - z$ curve to be zero at the upper end. Because of \n",
    "\n",
    "$$\n",
    "        0 = u'(z=H) \\approx \\frac{u_9 - u_8}{\\Delta z}\n",
    "$$\n",
    "\n",
    "this means that the displacements of the two final vertices should be equal. Thus, we want the final equation to yield $u_8 - u_9 = 0$ which can be achieved by manipulating the RHS and System matrix such that the first equation reads\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "0 & \\dots & 0 & 1 & -1\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "u_0\n",
    "\\\\\n",
    "\\vdots\n",
    "\\\\\n",
    "u_7\n",
    "\\\\\n",
    "u_8\n",
    "\\\\\n",
    "u_9\n",
    "\\end{pmatrix}\n",
    "= \n",
    "0\n",
    "$$\n",
    "\n",
    "Nota bene: if a non-zero displacement were given as the boundary condition, the RHS entry would have to be set to that given value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#upper bc\n",
    "A[-1,-1] = -1.\n",
    "RHS[-1] = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check whether the manipulations were successful:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.],\n",
       "       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -1.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.00019523, 0.00019523, 0.00019523, 0.00019523,\n",
       "       0.00019523, 0.00019523, 0.00019523, 0.00019523, 0.        ])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RHS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's solve the system again with the newly implemented boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAHvCAYAAACIbWPJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABIJklEQVR4nO3dd3xcV4H3/89RteTeLdtykROXxN1yikuKjUMSktDDEgIsLLD8aNvobKE8S2f32QQeWB7gYalLQk1CepzixE7i3uISd9mSbbkXSVY7vz9G1hrjLo2uRvq8X695TXTu3PFXJ+PR11dn7g0xRiRJkiSlT1bSASRJkqSOztItSZIkpZmlW5IkSUozS7ckSZKUZpZuSZIkKc1ykg6QbjfffHN89NFHk44hSZKkji+cbUOHP9K9b9++pCNIkiSpk+vwpVuSJElKmqVbkiRJSjNLtyRJkpRmlm5JkiQpzSzdkiRJUppZuiVJkqQ0s3RLkiRJaWbpliRJktLM0i1JkiSlmaVbkiRJSjNLtyRJkpRmlm5JkiQpzRIv3SGEohDCoyGEmHQWSZIkKR0SLd0hhDcCi4BR53lctxDCt0MIG0IIr4QQHg8hXNk2KSVJkqSWSfpI96eBecAL53nc/cAUYEqM8QrgJeCZEMKQNOeTJEmSWizp0j0zxvjquR4QQpgH3Az8U4yxqmn4S0A28Nk055MkSZJaLNHSHWOsv4CHvRmoA54/Zb9aUkfH35ymaC0SY6ShMXWTJEmSkj7SfSEmAuVNRftUW4GBIYQBCWQ6p32HanjDJx7gff/r8aSjSJIkqR3IhNLdDzh6hvEjTfd9T98QQvhACGFJCGFJZWVlWsNJkiRJ55MJpftswtk2xBi/H2MsjTGW9u/fvy0zSZIkSX8mE0r3PqD7GcZPju1vwyySJEnSRcuE0r0KGBxCyDttfCSwJ8a4N4FMkiRJ0gXLhNL9WyAXmHFyoKmAzwB+k1QoSZIk6UK1+9IdY3wceAz4UgihsGn4c0Aj8OXEgkmSJEkXKCfJPzyE8A1SV6Qc1vT1iqZNV512isC3Al8FVoQQGoCdwA0xxl1tGFeSJEm6JImW7hjjJy7wcUeBD6c5jiRJkpQW7X55iSRJkpTpLN2SJElSmlm6JUmSpDSzdEuSJElpZumWJEmS0szSLUmSJKWZpVuSJElKM0u3JEmSlGaWbkmSJCnNLN2SJElSmlm6JUmSpDSzdEuSJElpZumWJEmS0szSLUmSJKWZpVuSJElKM0u3JEmSlGaWbkmSJCnNLN2SJElSmlm6JUmSpDSzdEuSJElpZumWJEmS0szSLUmSJKWZpVuSJElKM0u3JEmSlGaWbkmSJCnNLN2SJElSmlm6JUmSpDSzdEuSJElpZumWJEmS0szSLUmSJKWZpVuSJElKM0u3JEmSlGaWbkmSJCnNLN2SJElSmlm6JUmSpDSzdKdRQ2NMOoIkSZLaAUt3GvTqnk/XglwOHj3B1vLDSceRJElSwizdaZCbk8UNU4cC8MTLOxJOI0mSpKRZutNk3lXDAHhmaRl19Q0Jp5EkSVKSLN1pMmpoL0oG9+RoVR0vrtmddBxJkiQlyNKdRvOuTh3tfuKl7QknkSRJUpIs3Wl0/dSh5OZkseLVSvYerEo6jiRJkhJi6U6j7oV5XDu+iBjhqcVlSceRJElSQizdafaapg9UPrl4B42et1uSJKlTsnSn2aTL+zOgdwF7D1SxetO+pONIkiQpAZbuNMvKCrxmetMHKj1ntyRJUqdk6W4Dc6cPIwRYuLqcY1W1SceRJElSG7N0t4EBfQqZdHl/6uobeXbZzqTjSJIkqY1ZutvIyStUPrHYJSaSJEmdjaW7jVwzvohuBbls3nmYLbsOJx1HkiRJbcjS3UbycrO5YepQAJ542StUSpIkdSaW7jY07+rhADyzdCe1dQ0Jp5EkSVJbsXS3oZIhPSkZ0pNj1XW8tGZ30nEkSZLURizdbeympg9UPu4SE0mSpE7D0t3Grp86lNycLFa+WsneA1VJx5EkSVIbsHS3sW6FeVw7oYgY4SlPHyhJktQpWLoTcNNVqQ9UPrl4B42NMeE0kiRJSjdLdwImXNaPAX0K2XuwmlWbKpOOI0mSpDSzdCcgKyvwmulNV6h82SUmkiRJHZ2lOyFzpxcTAixaXcGxqtqk40iSJCmNLN0JGdC7kMmX96euvpFnlu1MOo4kSZLSyNKdoHlNH6h0iYkkSVLHZulO0DUTBtG9MJctuw6zeeehpONIkiQpTSzdCcrNyeb6qUMBeNKj3ZIkSR2WpTthN12dWmLyzLKd1NY1JJxGkiRJ6WDpTtjIwT25bGhPjlXX8eKaiqTjSJIkKQ0s3e3Aa05+oPIll5hIkiR1RBlRukMIpSGER0II60IIq0MIL4cQ3pp0rtZy/dSh5OVksXJTJXsOVCUdR5IkSa2s3ZfuEMII4ClgHzAhxjgB+BFwXwjh9iSztZZuBblcO2EwMcJTiz3aLUmS1NG0+9IN3Ar0AP4txlgPEGP8HnAEuCvJYK1p3tWpy8I/uXgHjY0x4TSSJElqTZlQuuub7nNODoQQAqns2YkkSoMJo/oxsE8hlQerWflqZdJxJEmS1IoyoXT/N7Ae+McQQrcQQhbwWSAf+N6ZdgghfCCEsCSEsKSyMjMKbFZW4DVXNR3t9pzdkiRJHUq7L90xxiPAXKALqXXde4H3AvNijPPPss/3Y4ylMcbS/v37t13YFppTWkwIsGhNBUerapOOI0mSpFbS7kt3CGEM8DKwHegDDAA+B/w2hHBLktla24DehUwZPYC6+kaeWboz6TiSJElqJe2+dANfAnoBfxNjrIoxNsYY/xt4DvivEELOOffOMM0fqHSJiSRJUoeRCaV7ArAzxlh92vhGoD8wsu0jpc/VVw6ie2EeW8oPs2nnoaTjSJIkqRVkQuneCxSd4Yj2cCACB9s+Uvrk5mRz47ShgEe7JUmSOopMKN33kjpP9xebThVICOFG4E3Ar2KM+5IMlw4nz2LyzLKdnKhrSDiNJEmSWqrdl+4Y46+Bm4FrgVdCCGuA/yD1Ycq/TDBa2owc3JPLintxvLqOF1dXJB1HkiRJLZQRH0KMMT4GPJZ0jrY076phbCo7xO+f3cSMiYPJzWn3/z6SJEnSWdjk2qkbpg6lX88ubNp5mB8/tDbpOJIkSWoBS3c7Vdgll0+9ezo52YEHFmxhwYpdSUeSJEnSJbJ0t2Njh/fhvbePB+De+5ZTtudowokkSZJ0KSzd7dxts0Zy3eQhVJ9o4Ks/WUzNifqkI0mSJOkiWbrbuRACH7lzMsUDu7Fj91G+ff9KYoxJx5IkSdJFsHRngIL8HD7z7qvokpfNs8t38vDCbUlHkiRJ0kWwdGeI4oHd+eidkwH4wR9Ws3FHh7oQpyRJUodm6c4g100Zym0zR1LfEPnqTxZz5Hht0pEkSZJ0ASzdGea9d4xnzLDeVB6s5lu/WEpjo+u7JUmS2jtLd4bJzcnik+8qpXthHsvW7+VXT25MOpIkSZLOw9KdgQb0LuTjd08jBPjl4+tZtmFv0pEkSZJ0DpbuDDV1zADePm8MMcI3f7aUyoPVSUeSJEnSWVi6M9jb5o1h6pgBHK2q5Ws/WUxdfWPSkSRJknQGlu4MlpUV+Pu7ptKvVwEbdhzkRw+uSTqSJEmSzsDSneF6dsvnM++eTk524KHnt/Lc8p1JR5IkSdJpLN0dwOhhvXnf6ycAcO99KyjbczThRJIkSTqVpbuDuHXGCK6fMpSa2ga+8l8vU32iPulIkiRJamLp7iBCCHzkrZMoHtidsj3H+PZ9K4jRC+dIkiS1B5buDqRLfg6fefd0CvKzeW7FLv74wtakI0mSJAlLd4dTPLA7H71zCgA/fGAN67cfSDiRJEmSLN0d0OzJQ7h9dgn1DZGv/WQJh4+dSDqSJElSp2bp7qDec9uVjB3em32HqvnWz5fS0Oj6bkmSpKRYujuo3JwsPvWu6fTomsfyjZX86okNSUeSJEnqtCzdHVi/XgV8/B3TCAH++4kNLF2/J+lIkiRJnZKlu4ObMmYAd712LDHCt36+lL0Hq5KOJEmS1OlYujuBO+eOZtrYARytquNrP1lMXX1D0pEkSZI6FUt3J5CVFfj7u6bRv3cBG3cc4gd/WJN0JEmSpE7F0t1J9Oiax6ffNZ2c7CweXriNZ5btTDqSJElSp2Hp7kRGD+vNB94wHoBv37+C7buPJJxIkiSpc7B0dzI3XzuCG6YN5URtA1/58WKqauqSjiRJktThWbo7mRACH37zJIYN6s6uymPce98KYvTCOZIkSelk6e6EuuTn8Jl3T6cgP4fnV5bz4PNbko4kSZLUoVm6O6mhA7rzN2+bAsCPHljL+m0HEk4kSZLUcVm6O7GZkwZzx3UlNDRGvvaTxRw+diLpSJIkSR2SpbuTe89tVzJuRB/2Ha7hmz9bSkOj67slSZJam6W7k8vJzuJT7yqlZ7c8VrxayS8fX590JEmSpA7H0i369izgE+8oJSvAr57YyJJ1e5KOJEmS1KFYugXApNH9uevmsQB86+dL2XOgKuFEkiRJHYelW83eOmc0peMGcqy6jq/+ZDF19Q1JR5IkSeoQLN1qlpUV+Pu7pjKgTyGbyg7xf3+/JulIkiRJHYKlW3+ie2Een3nXdHKys3hk0TaeXlqWdCRJkqSMZ+nWn7msuBd//cYJAHz7/pVsrziScCJJkqTMZunWGb32muHMKS2mtq6Br/zXy1TV1CUdSZIkKWNZunVGIQT+vzdPZERRD3ZVHucLP3iRI8drk44lSZKUkSzdOqsueTl85i+n07dnF17ZeoBP3PMc5fuOJR1LkiQp41i6dU6D+3XjW39zHSWDe1K+7zifuGcB67cdSDqWJElSRrF067z69izgKx+eybSxAzhyvJbPfvcFnl+5K+lYkiRJGcPSrQtS2CWXf3rv1dxy7Qjq6hv52k+W8Jv5rxJjTDqaJElSu2fp1gXLzs7i/3vzRN5z2xUA/PiPr/Dd36yioaEx4WSSJEntm6VbFyWEwJtuvJxPvauU3JzUBXS+9KOXPKWgJEnSOVi6dUlmTRrCv35wJj265rF0/V4+/Z3n2X+4OulYkiRJ7ZKlW5ds3Mg+fPNj1zGkf1e2lh/hH/7jObaWH046liRJUrtj6VaLFPXrytc/eh1XjOzD/sM1fOrbC1i6fk/SsSRJktoVS7darEfXPL701zO4bsoQqk808MUfvsSji7YlHUuSJKndsHSrVeTlZvMPd03jzteMprEx8p1fr+THD62lsdFTCkqSJFm61WqysgLvvGUcH71zMllZgd88vYlv/GwJtXUNSUeTJElKlKVbre6mq4fz+fddQ2GXHJ5fWc4/fm8hh4+dSDqWJElSYizdSospYwbwtY/Mpl+vAtZtO8An7l1AeeWxpGNJkiQlwtKttBlR1INvfmw2JUN6UrHvOB+/ZwFrt+xPOpYkSVKbs3Qrrfr2LOCrH55F6biBHK2q5R+/t5AFy3clHUuSJKlNWbqVdgX5Ofzje67i1hkjqG9o5Os/W8L9T20kRs9sIkmSOgdLt9pEdnYWH3zTRP7qjisJAX7y8Dq+ff9K6hsak44mSZKUdpZutZkQAm+4/jI+/a7p5OVk8fhL2/niD16kqqYu6WiSJElpZelWm5sxcTD/+qGZ9OyWx/KNlXzq289TebA66ViSJElpY+lWIsYO78M3P3YdQ/p3Y1vFET5+z3Ns3nko6ViSJElpYelWYgb17co3Pjab8aP6cuBIDZ/+zvMsWbcn6ViSJEmtLmNKdwjhzSGE50IIS0MIW0IIS0II70w6l1qme2EeX/zAtdwwdSg1tQ186Ycv8vDCrUnHkiRJalUZUbpDCH8HfA64K8Y4DRgDbATmJhpMrSI3J5u/v2sqb5s3msYI3/3NKn704FoaGz2loCRJ6hhykg5wPiGEEcBXgVkxxp0AMca6EMLHgcFJZlPrCSFw983jGNSnkG/fv5LfPbOJvQeq+Lu7ppKfm510PEmSpBbJhCPd7wQOxRgXnzoYYyyPMS5JKJPS5DVXDecL77+Wwi45vLCqnM999wUOHzuRdCxJkqQWyYTSPQPY1rSme0EIYX0IYWEI4b1n2yGE8IGmNd9LKisr2zCqWsOk0f35+kdn0793ARu2H+Tj9zzHzr1Hk44lSZJ0yUJ7vxR3CGENMAJYDbwR2Au8Gfgl8C8xxn891/6lpaVxyRIPiGeiA0dq+NIPX2TTzsN0K8jlc++5ivGj+iUdS5Ik6WzC2TZkwpHuLkBX4BMxxt0xxsYY4/3AH4DPhhAKk42ndOnTowtf+dAsrrpiEMeq6/in/1zEM8t2Jh1LkiTpomVC6T65rmDFaePLgULgijZNozbVJT+Hz77nKm6bNZL6hka+9fOl/OrJDbT339BIkiSdKhNK9/qm+9OzNpxlXB1Mdlbgr984kfe/fjwhwM8eWc+9962gvqEx6WiSJEkXJBMK64NN9xNPGx8PVANr2zaOknLHdaP4zLuvIi83myde3sEX/u+LHK+uSzqWJEnSeWVC6f4VsBj4XyGEbgAhhNnAW4B/jTEeTzKc2ta1E4r4yodm0qtbPitereST317A3oNVSceSJEk6p3ZfumOMDcDNwKvA2hDCBuA7wEfOd+YSdUyjh/XmGx+bTfHAbuzYfZSP/8dzbCo7lHQsSZKks2r3pwxsKU8Z2HEdq6rlyz9ezOrN+8jPy+aTd5dy1ZWDko4lSZI6r4w+ZaB0Rt0K8/jCB67lxmlDOVHbwL/+v5f44/Nbko4lSZL0Zyzdymi5OVn83dunctdNY2iM8L3freYHf1hDQ2PH/g2OJEnKLJZuZbwQAm9/7Vj+9i+mkJMd+MNzm/naTxZTU1ufdDRJkiTA0q0OZO70YXzhA9fStUsOi1ZX8LnvvsDBozVJx5IkSbJ0q2OZeFl/vv7R2QzoXcDGHYf4+D0LKNtz9Pw7SpIkpZGlWx3OsEE9+ObHruPy4l7sPVDFJ+5dwOpN+5KOJUmSOjFLtzqk3j268OUPzeTqKwdxvLqOf/7+Qp5eWpZ0LEmS1ElZutVhdcnL4TN/eRV3XFdCfUPk336xjF8+tp6Ofm56SZLU/li61aFlZwXe//oJfOANE8gK8IvHN/CP31tIxb7jSUeTJEmdiKVbncLts0v43Huvpme3PFZt2sdHvvk0v392k+fzliRJbcLLwKtTOXzsBD/4wxqeWbYTgNHDevGxO6cwvKhHwskkSVIHcNbLwFu61SktfmU3/+fXK9l3uIac7MBb547mrXNHk5vjL38kSdIls3RLp6uqqePHD73CI4u2ATBsUHf+5m1TGD2sd7LBJElSpjpr6W7Vw3ohhJ+25vNJ6VTYJZcPvWUSX/7QTIr6dWXH7qN84p7n+OEDa7yEvCRJalUXfaQ7hHAZcD0wCMg+bfMHY4yDWylbq/BIty5ETW09v3xsA79/dhONEQb1LeSjd05m4mX9k44mSZIyR+ssLwkhfAi49xxPGGOMpxfxRFm6dTFeLTvIPb9awbaKIwC89prhvOe2K+lakJtwMkmSlAFabXnJJ4EPA/1jjFmn34BVLUkpJe3y4t78299ez903jyUnO4vHXtzOh74+n5fX7k46miRJymAXe6R7cYxx+jm2j48xrmmVZK3EI926VDt2H+Ge+1awYftBAK6bPIT3v2ECvbrnJ5xMkiS1U612pPt3IYQJ59j+1xf5fFK7NWxQD772kdm8//Xjyc/L5rkVu/jQ1+fzzNIyLyUvSZIuyqV8kPJrwChgE1B12mY/SKkOaff+43zn/pWseLUSgNJxA/nQmyfRv3dBwskkSVI70qofpPz2OR7iBynVYcUYefLlHfzwgTUcr6mnID+H99x2Ba+9ZgRZWWf9OyZJkjqPVivdO4EvAL+JMR44w/blMcYplxQxTSzdam37D1fzvd+u4sU1qQ9Xjh/Vl4++dTKD+3dLOJkkSUpYq5XuZTHGqefY7gcp1SnEGHlhVTn/+dvVHDp2grycLN5x81hef90osrO9lLwkSZ1Uq32Q8gk/SClBCIFZk4bwnU/OYU5pMbX1jfy/h17h4/c8x9byw0nHkyRJ7czFHun+IvAOYDl+kFJqtnT9Hr7z65VUHqwmOyvwljmX87Z5o8nNaVcfcZAkSenVastLGs/zED9IqU6rqqaOnzy8jj++sBWA4oHd+NidUxg7ok/CySRJUhtptdJ9zg9K+kFKCdZu2c+99y1nV+VxQoDbZ5XwzlvG0SU/J+lokiQpvVptTfc/n2f7Ry/y+aQO58qSvtzzDzfyljmXE0LggQVb+PA3n2bFxr1JR5MkSQm56IvjZBqPdCtJm3Ye4t5frWBL04cr5101jPfeMZ5uBbkJJ5MkSWnQOstLMpGlW0mrb2jkt09v4pePb6C+oZE+PfL54Jsmce2EoqSjSZKk1mXplpJWtuco9963gnXbUteVmjlpMH/9xgn07t4l4WSSJKmVWLql9qCxMfLwwq381x9foaa2gW4Fubz/DeO5cVoxIXgpeUmSMpylW2pP9hyo4jv3r2D5xkoApo4dwIffMokBvQsTTiZJklrg0kt3CGEmsDBmaDu3dKu9ijEyf0kZP/jDGo5V11GQn827b72CW2aMJCvLo96SJGWgFpXuKlJXnnwUeAh4NMZ4qDXTpZOlW+3dwSM1fO93q1i4qgKAK0b24aN3TmbogO4JJ5MkSRepRaW7AJgHvA64FRgILCJVwB+KMa5rvZytz9KtTPHCqnK+99tVHDp6gtycLN5+0xjedMNlZGdf7On0JUlSQlpvTXcIYQqpAn4bUArsoKmAA8/EGGsvPWfrs3QrkxyrquWHD6zlycU7ABg1tCcfu3MKJUN6JpxMkiRdgPR8kDKE0J9UAX8dqaPh2cCTwIMxxh9d8hO3Iku3MtGyDXv5zv0r2HuwmqyswJtvvIy/mDeGvNzspKNJkqSzS//ZS0IIOcD1pI6A3xJjHNsqT9xClm5lquoT9fz0kXU89PwWYoQh/bvxsbdN5oqRfZOOJkmSzsxTBkqZ6pWt+7n3vhXs3HuMEOB1M0fyrluvoCA/J+lokiTpT1m6pUxWW9fAr57cyG/mv0pDY2RA7wI+/NbJTB0zIOlokiTpf1i6pY5gy67D3HPfcjbvPAzA3OnF/NUd4+lemJdwMkmShKXb0q2Oo6Ghkd89u5lfPLaeuvpGenXP54NvmsjMiYOTjiZJUmd31tJ9UScADiF8J4RwXcvzSLpU2dlZvGXO5dz78Ru5sqQvh46e4Kv/tZgv//hlDhypSTqeJEk6g4s60h1CaAROAK+PMT5+2rZuwKAY46bWjdgyHulWR9bYGHlk0Tb+649rqT7RQNeCXN53x3jmTi8mBC8lL0lSG2udI91NHgV+F0K48bTxK4ENl/B8ki5RVlbgdTNH8u1PzGHq2AEcr67jP361nH/+/iL2HKhKOp4kSWpyKaX7K8C3gAdDCLNbOY+kSzCgdyGff981/P1dU+lemMuKjZV85BvzeXDBFhobO/bnNiRJygSXUrqJMf4z8G3gjyGEa1s3kqRLEULgxmnF/J9PzmXWpMHU1Dbw/d+v5tPfeZ6yPUeTjidJUqd2SaUbIMb4aeD/Ao+EEKa3XiRJLdGrez6fetd0PvuXV9GnRz7rth3gY996hvue3Eh9Q2PS8SRJ6pQuuXQDxBj/Afgp8BgwtVUSSWoV104o4jufnMu8q4ZR39DITx9Zx9//72fZtPNQ0tEkSep0LvbsJTOAdTHGg6eNfxd4H5AVY8xu3Ygt49lLJFi5sZJ771/BngNVZGUF3nj9KN7+2rHk57arv66SJGW69F8cJ4RwD/CuGGOvVnnCVmLpllJqTtTzs0fX88CCzcQIQ/p35cNvmcyEy/olHU2SpI7CK1JKSlm//QD3/GpF84crp44dwDtvHsdlxb2SDSZJUuazdEv6H3X1Dfzm6U389ulXqT7RAKTWgL/j5rEMH9Qj4XSSJGUsS7ekP3f42Al+8/Qm/vj8FmrrGwkBrp86lLtuGktRv65Jx5MkKdNYuiWd3f7D1dz35EYef2k79Q2R7KzAa64axl/MG0O/XgVJx5MkKVNYuiWd354DVfzy8fU8vaSMxgi5OVncMmMEb50zml7d85OOJ0lSe2fplnThyvYc5RePref5leUAdMnL5vbZJbzphsvoVpiXcDpJktotS7eki7dl12F+9ug6Fr+yB4CuBbm88YZR3DF7FAX5OQmnkySp3bF0S7p067cd4KePrGPVpn0A9OyWx1vmjObWGSPI8wI7kiSdZOmW1HIrX63kp4+sY8P21EVp+/bswtvmjWHeVcPIyc5KOJ0kSYmzdEtqHTFGFq/bw88eWcfW8iMADOpbyNtvGsv1U4eSnXXW9xtJkjo6S7ek1tXYGHlhVTk/f3Q9uyqPAVA8sDvvuHksMyYUEYLlW5LU6Vi6JaVHQ0MjTy/dyS+f2MDeA1UAjBrak7tvHse0sQMs35KkzsTSLSm96uobefyl7dz35AYOHDkBwBUj+3D3LeOYMKpfwukkSWoTHat0hxAWALOAkTHGbed6rKVbals1tfU8/MI2fj3/VY5W1QIweXR/3nnLOEYP651wOkmS0qrjlO4QwpuBXzd9aemW2qmqmjr+8NwWfv/sJqpq6gG4+spB3H3LOEYU9Ug4nSRJadExSncIIQ9YA7wK3IqlW2r3jhyv5bdPv8qDz2+ltq6BEGD25CG847VjGdy/W9LxJElqTR2mdP8dMB3YCPwLlm4pYxw8UsN9T23k0UXbqW9oJCsrMLe0mL+4aQwDehcmHU+SpNaQ+aU7hNCH1FHuGcBfYumWMtLeA1X89xMbeGpJGY2NkZzsLG6+djh3zh1N7x5dko4nSVJLdIjS/b+B2hjjJ0MIn+ccpTuE8AHgAwDDhg2btn379jZMKulC7Ko8xi8eW8+CFbuIEfLzsrl9VglvuvEyuhfmJR1PkqRLkdmlO4RwGfAcMC7GePh8pftUHumW2rdtFUf42SPreGntbgAKu+Twhusv4/XXlVDYJTfhdJIkXZSML92/BZ6JMd7T9PXnsXRLHcrGHQf56SPrWLGxEoDuhXm8Zc7lvG7WSPJzsxNOJ0nSBcnc0h1CmA38ELgyxljXNPZ5LN1Sh7R68z5++vA61m07AECfHvm8bd4Y5l01nNycrITTSZJ0Thldur8IvAM4esrwIGAgsA6oBT4bY3z4TPtbuqXME2Nk6fq9/OzRdWzeeRiAAX0KueumMdwwdSjZ2ZZvSVK7lLml+0w80i11DjFGFq6u4OePrqNszzEAhg7oxl2vHcvMiYPJyjrre5skSUk46w+mnLZMIUkXI4TAzImDuWZ8Ec8u28kvH1/Pzr3H+PpPl1AyuCd33zKW0nEDCcHyLUlq3zLqSHcI4Vbgy5y2vCTGOPls+3ikW+o46hsaeeLlHfzqiQ3sP1wDwNjhvXnnreOYeFn/hNNJktTBlpdcDEu31PGcqGvgkYXb+PX8jRw+VgvApMv7cfct4xg7vE/C6SRJnZilW1LHU32ingcWbOZ3T2/ieE09AFddMYi7bxnLyME9E04nSeqELN2SOq5jVbX89plNPLhgCzW1DQDMnjyEu147hqEDuiecTpLUiVi6JXV8h46e4P75G3lk4Tbq6hvJCjCndBh/cdMYBvYpTDqeJKnjs3RL6jwqD1bzqyc38OTLO2hojORkB157zQjufM1o+vToknQ8SVLHZemW1PlU7DvOLx5fz7PLdhIj5OVmc9vMkbx5zuX06JqXdDxJUsdj6ZbUeW3ffYSfP7qeRasrACjIz+EN14/iDdePorBLbsLpJEkdiKVbkjaVHeKnj65j2fq9AHQvzOXNN17O62aNpEue1wqTJLWYpVuSTlq7ZT8/fWQda7fsB6B393zufM1oXnvNcHJzshNOJ0nKYJZuSTpVjJHlGyv56SPr2FR2CID+vQt4+7wxzCktJjs7K9mAkqRMZOmWpDOJMfLimt38/NF1bN99FIAh/bty12vHMmvSELKyzvr+KUnS6SzdknQuDY2RBSt28YvH1lOx7zgAI4p6cPfNY7nqykGEYPmWJJ2XpVuSLkR9QyNPLS7jv5/YwL5D1QCMHtaLd94yjkmX97d8S5LOxdItSRejtq6BR1/cxv1PvsqhYycAmDCqH++8ZRzjRvZJOJ0kqZ2ydEvSpag5Uc+Dz2/ht09v4lh1HQCl4wZy981jGTW0V7LhJEntjaVbklriWHUdv392Ew88t5nqEw0AzJw4mHfcPJbigd0TTidJaics3ZLUGg4fO8Gv57/Kwy9spba+kawAN0wr5u03jWFQ365Jx5MkJcvSLUmtaf/han715EYef3E7DY2R7KzATVcP5003Xmb5lqTOy9ItSemwe/9xfvn4Bp5ZWkZjhKwAV105iDtmj2L8qL6e7USSOhdLtySlU9meo/x6/qs8t3wX9Q2NQOo837fPLuH6qUPJz/Xy8pLUCVi6JaktHDxaw6MLt/Hwom0cOpo61WD3wjxuvnY4r5s5kr49CxJOKElKI0u3JLWluvoGFqwo58EFm9m08zAA2VmBmRMHc/t1JYwd7rm+JakDsnRLUhJijKzfdpAHFmxm4eoKGhtT77mjh/Xi9lklzJw0hNycrIRTSpJaiaVbkpJWebCahxdu5bEXt3G0KnWhnT498rllxkhuvmYEvbrnJ5xQktRClm5Jai9qaut5dtkuHlywme27jwKQk53F9VOHcPusEq90KUmZy9ItSe1NjJFVm/bx4IItvPzKbk6+HV9Z0pfbZ5dwzZWDyM526YkkZRBLtyS1ZxX7jvPQC1t48uUdVNXUA9C/dwG3zRzJvKuH070wL+GEkqQLYOmWpExQVVPH/CVlPLhgC+X7jgOQn5fNjdOKuX3WSIYN6pFwQknSOVi6JSmTNDZGlm3YywPPbWb5xsrm8cmj+3PH7BKmjR1IVpZXu5SkdsbSLUmZqmzPUR5csIX5S8s4UdsAwOB+XbltVglzpxdT2CU34YSSpCaWbknKdMeqann8pR388YUt7D1YDUBBfg7zrhrGbbNKKOrXNeGEktTpWbolqaNoaGjkpbW7eWDBFtZu2Q9ACDB93CDumF3CxMv7EYJLTyQpAZZuSeqItuw6zIMLtvDs8p3U1TcCMHxQd26fXcL1U4fSJS8n4YSS1KlYuiWpIzt09ASPvbiNhxdu5cCREwB0L8zlpquH87qZJfTvXZBwQknqFCzdktQZ1NU38sKqch5csJmNOw4BkJUVuHZCEXfMLmHciD4uPZGk9LF0S1Jns377AR5csIUXVpbT0Jh6r79saE9un13C7MlDyM3JTjihJHU4lm5J6qz2H67m4YXbeHTRNo4crwWgV/d8brl2BLdcO4LePboknFCSOgxLtyR1difqGnhu2U4eWLCFbRVHAMjJDsyePIQ7Zo/isuJeyQaUpMxn6ZYkpcQYWbN5Pw8s2MzLa3fTtPKEcSP6cPvsEmZMKCI7OyvZkJKUmSzdkqQ/t3v/cf74wlaeeGk7x2vqAejXswu3zhzJa68ZQY+ueQknlKSMYumWJJ1d9Yl65i8p48EFW9hVeQyAvNxsbpw2lNtnlTC8qEfCCSUpI1i6JUnn19gYWbGxkgcWbGbp+r3N4xMv68cds0sovWIQ2VmeclCSzsLSLUm6ODv3HuWh57fy1OId1NQ2ADCobyG3zSrhNdOH0bUgN+GEktTuWLolSZfmWHUdT768nYee38qeA1UAFORnM3f6MG6fVcLg/t0STihJ7YalW5LUMg2NkcWv7ObBBVtYtWlf83jpuIHcPruEKaP7e7VLSZ2dpVuS1Hq2lh/mwQVbeHbZTmrrGwEoHtiN22eVcOO0Yrrk5yScUJISYemWJLW+w8dO8NiL23l44Vb2H64BoGtBLjddPZzbZo5kQJ/ChBNKUpuydEuS0qe+oZFFqyp4YMFm1m8/CEBWgKvHF3HH7BKuLOnr0hNJnYGlW5LUNjbuOMiDC7bw/Mpd1DekfsaUDO7J7bNLuG7KEPJysxNOKElpY+mWJLWtA0dqeGThNh5dtI1Dx04A0LNbHjdfM4JbZoygb8+ChBNKUquzdEuSklFb18CCFbt4YMEWtuw6DEB2VmDWpCHccV0Jo4f1TjihJLUaS7ckKVkxRl7ZeoAHFmzmxdUVNDb9+BkzvDd3zC5hxsTB5GRnJRtSklrG0i1Jaj/2Hqzi4Re28tiL2zlWXQdAnx5duHXmCG6+ZgQ9u+UnnFCSLomlW5LU/tScqOfpZTt5cMEWyvYcBSA3J4sbpg7l9tkljBzcM+GEknRRLN2SpPYrxsiKjZU8sGALS9btaR4fP6ovd8wu4aori8jO8pSDkto9S7ckKTOUVx7joRe28uTLO6g+UQ/AgD6F3DZzJPOuHk63gtyEE0rSWVm6JUmZpaqmjidf3sFDz2+lYv9xAPLzsplTWszts0ooHtg94YSS9Gcs3ZKkzNTQGFm6bg8PLNjMylf3NY9PHTOA22eXMHXMALJceiKpfbB0S5Iy3/bdR3hwwRaeXrqT2roGAIb078bts0YyZ/owCvJzEk4oqZOzdEuSOo4jx2t5/KXt/PGFrew7VA1A1y45zLt6OK+bOZJBfbsmnFBSJ2XpliR1PA0Njby4ZjcPLNjMK1sPABACXHXFIO64roQJo/oRgktPJLUZS7ckqWPbVHaIB5/fwnPLd1Hf0AjAiKIe3DpzJLMnD/GsJ5LagqVbktQ5HDxaw6MLt/Hwom0cOnoCgLycLK4eX8Sc0mKmjO5Ptpebl5Qelm5JUudSV9/A8yvLefLlHazevI+TP+769Mjn+qnFzC0tZnhRj2RDSupoLN2SpM5r78Eqnl5axvzFZZTvO948PmpoT+aUFnP9lKH07JafYEJJHYSlW5KkGCMbth/kqSVlLFi+k+M1qSteZmcFSscNZO70YkrHDSI3x+Unki6JpVuSpFPV1jXw0trdzF9SxrL1e2hs+nHYvTCP66cMYe70YYwa2tOzn0i6GJlbukMIk4EPA1OBHCAXeBL4Uoyx8nz7W7olSedz8EgNzyzbyfwlZWyrONI8PmxQd+aWFnP91KH07VmQYEJJGSKjS/d6YC3wrhjj8RDCEOApIAuYFGOsPtf+lm5J0oWKMbJl12HmLynjmWU7OXK8FoCsAJNHD2BOaTHXTCgiPzc74aSS2qmML923xRg3nTL2V8APgLfEGH9zrv0t3ZKkS1Hf0MjSdXt4akkZi1/ZTX1D6udlYZccZk0awpzSYq4Y2cflJ5JOddY3hJy2THGJJsYYa08bK2+6793WYSRJnUNOdurc3lePL+LI8VoWLN/JU0vKeLXsEI+/tJ3HX9pOUd+u3FhazJzSYgb2KUw6sqR2rN0f6T6TEMLfAP8OjI0xbjzD9g8AHwAYNmzYtO3bt7dxQklSR7Vj9xHmLynj6aU7OXCkpnl8/Ki+zC0dxoyJRRR28eqXUieVuctLThdCyAaWAS/HGN9/vse7vESSlA4NjZGVr1Yyf3EZi1aXU1ufuvR8fl42MyYUMbd0GBMu60dWlstPpE6kQ5XuzwO3A9fHGI+d7/GWbklSulXV1PH8ynLmLylj7Zb9zeP9ehVw47ShzCktZuiA7gkmlNRGOkbpDiG8B/gEcEOMce+F7GPpliS1pYp9x3l6aRlPLSlj74Gq5vExw3szp7SY6yYPoVthXoIJJaVR5pfuEMI7gU8Br4kx7r7Q/SzdkqQkNDZG1m7dz/zFZbywahfVJxqAkx/QHMTc0mKmjhlAdrZXv5Q6kMwu3SGEu4FPc0rhDiHcBgyOMX7/XPtauiVJSas5Uc+iNRXMX1zGyk2VnPzR26t7PjdMTS0/GTm4Z7IhJbWGzC3dIYR3kDon9z8Bpx7hng1UxBg/f679Ld2SpPak8mA1zywr46nFZeyq/J+PJpUM7smc6cVcP2UovbrnJ5hQUgtkdOk+wNnPx/0FS7ckKRPFGNm44yBPLSljwfJdHKuuAyA7KzBt7EDmTC/mqisGkpvj1S+lDJK5pbulLN2SpPaurr6Bl9fu4aklO1i6fi+Njamfzd0KcrluyhDmTh/G5cW9vPql1P5ZuiVJygQHj9bw7LJdzF+yg63lR5rHhw7oxpzSYm6cVky/XgUJJpR0DpZuSZIyzdbywzy1uIxnl+3k0LETAIQAky7vz9zSYq6ZUESXvJyEU0o6haVbkqRMVd/QyLINe5m/uIyX1u6mviF19cuC/BxmTRrMnNJirhjZ16tfSsmzdEuS1BEcraplwYpdzF9cxoYdB5vHB/YpZE5pMXNKixnUt2uCCaVOzdItSVJHU7bnKE8vLePpJWXsO1zTPH5lSV/mlBYza9JgCrvkJphQ6nQs3ZIkdVQNjZHVmyp5akkZi1ZXcKI2dfXLvNxsrh1fxJzpxUy6vD/ZLj+R0s3SLUlSZ1BVU8fCVeU8taSMNZv3N4/37dmFG6YOZe70YRQP7J5gQqlDs3RLktTZ7N5/nKeX7mT+kh3s3l/VPH55cS/mlhYze8pQenTNSzCh1OFYuiVJ6qxijLyy9QDzl5Tx/MpdVNXUA5CTHZh+xSDmlhYzbdxAcrKzEk4qZTxLtyRJghN1Dby4uoL5S8pYsXEvTRe/pGe3PK6fMpQ5pcWUDOnp1S+lS2PpliRJf2r/4WqeWbqTp5aUUbbnaPP4iKIezCkt5oapQ+ndo0uCCaWMY+mWJElnFmNk085DzF9cxrPLd3G0qhaArKzA1DEDmFNazNVXDiIvNzvhpFK7Z+mWJEnnV1ffyJJ1u3lqcRlL1u2hoWn9SdeCXGZPHsLc0mLGDO/t8hPpzCzdkiTp4hw+doJnl+9k/pIyNu883Dw+pH9Xbiwt5sZpxQzoXZhgQqndsXRLkqRLt73iCE8tKeOZpWUcPHoCgBBgwqh+zJ1ezLUTBlOQn5NwSilxlm5JktRyDQ2NLN9YyfwlZby4poK6+kYAuuRlM2PiYOZOL2Z8ST+yvPqlOidLtyRJal3Hqut4fsUu5i8pY922A83jA3oXcOO0YuZML2Zwv24JJpTanKVbkiSlT3nlMeYvKWP+0jIqD1Y3j48b0Yc5pcXMmjyEbgW5CSaU2oSlW5IkpV9jY2TNln08tbiMhavKqaltACA3J4trxhcxp7SYKaP7k+3VL9UxWbolSVLbqj5Rz6LV5Ty1uIzVm/dxsnL07p7PDdOKmVtazPCiHsmGlFqXpVuSJCVn78Eqnl5axvzFZZTvO948PmpoT+aUFnP9lKH07JafYEKpVVi6JUlS8mKMbNh+kKeWlLFg+U6O19QDkJ0VKB03kLnTiykdN4jcHJefKCNZuiVJUvtSW9fAS2t3M39JGcvW76Hp4pd0L8zj+ilDmDO9mMuG9vLql8oklm5JktR+HTxSwzPLUle/3FZxpHm8eGB35pYWc8O0ofTtWZBgQumCWLolSVL7F2Nky67DzF9SxjPLdnLkeC0AWQEmjx7AnNJirplQRH5udsJJpTOydEuSpMxSV9/I0vV7mL+kjMWv7Ka+IdVZCrvkMGvSEGZOGszEy/qR4+kH1X5YuiVJUuY6cryW55bv5KklZWwqO9Q83q0gl6vHD2LGxMFMGd2f3ByPgCtRlm5JktQx7Nh9hOeW72Lh6nLK9hxrHi/Iz2H6FQOZMXEw08YOoEteToIp1UlZuiVJUsdTtucoC1eVs3BVBVvKDzeP5+VmM23sAGZMHMxVVwyksIuXoFebsHRLkqSOrWLf8VQBX13Oxh2HmsdzsrOYMqY/MyYM5urxg+hemJdcSHV0lm5JktR5VB6sZtHqchauruCVrfubL0GfnRWYcFk/ZkwczDXjB9G7e5dkg6qjsXRLkqTO6eCRGl5cU8HCVRWs2ryPxqar8GQFuKKkLzMmDGbGxCLPA67WYOmWJEk6cryWl9ZUsHB1BSs27m0+DSHAmOG9mwv4oL5dE0ypDGbpliRJOtXx6joWv7KbhasrWLpuD7X1jc3bSob0ZObEVAEfOqB7gimVYSzdkiRJZ1Nzop6l6/fywqpylqzbTfWJhuZtwwZ1bz4CPqKoByGctVdJlm5JkqQLUVvXwPINe1m4uoKX1u7meHVd87aifl2ZMaGIGRMHc3lxLwu4TmfpliRJulh19Y2s3rSPhavLeXFNBYeP1TZvG9C7gGubjoCPHd6HrCwLuCzdSceQJEkZrqGhkbVb97NwVQWLVpdz4MiJ5m19euRzzfjUEfDxJX3Jzs5KMKkSZOmWJElqLY2NkQ3bD7JwdTkLV5Wz92B187YeXfO4+spBzJw0mImX9Sc3xwLeiVi6JUmS0iHGyOadh1m4upwXVpZTvu9487auXXK46spBzJg4mCljBpCfm51gUrUBS7ckSVK6xRjZvvto6nL0q8rZvvto87YuedmUjhvIjImDKR03kIL8nASTKk0s3ZIkSW1tV+Wx5gK+aefh5vG8nCymjh3AjImDueqKQXQtyE0wpVqRpVuSJClJu/cfZ9HqChauKmf99oPN4znZgUmX92fGxMFcfeUgenbLTzClWsjSLUmS1F7sP1zdVMArWLtlH41NdSwrKzC+pC8zJw3m2vFF9O7RJdmguliWbkmSpPbo0NETvLQ2VcBXvlpJQ1MDDwHGjejDjImDuXZCEQN6FyacVBfA0i1JktTeHauq5aW1u1m4qoLlG/dSV9/YvO3y4l7MmJi6GM/gft0STKlzsHRLkiRlkqqaOpas28PCVRUsWb+HE7UNzdtGDu6RKuATihg2qEeCKXUaS7ckSVKmqqmtZ/mGvbywsoKXX9lN9Yn65m1DB3RrLuAlQ3oSgpejT5ClW5IkqSOoq29gxcZKFq6q4KW1FRytqmveNqhvITMmpJagjB7W2wLe9izdkiRJHU19QyNrNu9j4aoKFq2u4NCxE83b+vXswrVNR8DHjexLdpYFvA1YuiVJkjqyhsbIuq37Wbi6gkWrytl3uKZ5W6/u+VwzvoiZE4sYP6ofOdlZCSbt0CzdkiRJnUVjY+TVsoMsXFXBwtXl7N5f1byte2EuV19ZxIyJRUwe3Z/cnOwEk3Y4lm5JkqTOKMbIll2HWdh0Ncyde481byvsksP0cYOYMbGIqWMH0CUvJ8GkHYKlW5IkSbBj95HmAr61/EjzeH5eNqVjBzJjYhGl4wZS2CU3wZQZy9ItSZKkP1W+7xiLmpagbNxxqHk8NyeLKaMHMGNiEVdfOYhuhXnJhcwslm5JkiSd3d6DVSxqOgK+btsBTlbE7KzAxMv6MWPiYK4ZX0Sv7vnJBm3fLN2SJEm6MAeO1PDimlQBX715P42Nqb6YFeDKkn7MmFjEtROK6NuzIOGk7Y6lW5IkSRfv8LETvLR2NwtXlbPy1UrqG/6nO44d3jt1NcyJgxnYpzDBlO2GpVuSJEktc6y6jsWvpAr4svV7qa1vbN42amhPZjYV8CH9uyWYMlGWbkmSJLWe6hP1LF2/hxdWlrNk3R5qahuatw0f1L35CPjwQd070+XoLd2SJElKjxN1DSzfsJeFq8p5ee1ujtfUN28b3K9rUwEv4rKhvTp6Abd0S5IkKf3q6htZtamShasqeHFNBUeO1zZvG9C7IFXAJwxmzPDeZGV1uAJu6ZYkSVLbamhoZM2W/SxcVc6Layo4cORE87Y+Pbpw7YTU5eivHNmX7OysBJO2Gku3JEmSktPYGFm//QALmy7GU3mwunlbz255XDO+iBkTBjPhsn7k5mRsAbd0S5IkqX2IMbJp5yEWrqrghVXlVOw73ryta0EuV185iBkTipgyZgB5udkJJr1olm5JkiS1PzFGtlUcaT4CvmP30eZtBfnZlI4bxMyJg5k2dgBd8nMSTHpBMrt0hxAGAP8OlDYNrQb+Nsa483z7WrolSZIyx869R5sL+Oadh5vH83KzmTZ2ADMmFDH9ikF0LchNMOVZZW7pDiHkAYuBjcDbgQj8CJgBTIkxHjvX/pZuSZKkzLR7//HmAr5h+8Hm8ZzsLCaP7s+MCUVcPb6IHl3zEkz5JzK6dL8f+D4wKsa4pWlsELAL+HSM8Rvn2t/SLUmSlPn2Hapm0epUAX9ly34amypsVlZg4qh+zJhYxDXji+jdo0uSMTO6dD8KjIsxDj9tfDVwPMZ4zbn2t3RLkiR1LIeOnuDFNRUsXFXOqk37aGhq4CHAFSP7Mre0mHlXDz/Ps6TFWUt3u1+NDkwktbTkdFuBuWfaIYTwAeADAMOGDUtfMkmSJLW5Xt3zufnaEdx87QiOVtXy0prdLFxdzvINlazdsp/B/bomVbrPKhOOdNcCj8UYbz9t/GfAO4DCGGP1GXfGI92SJEmdRVVNHYtf2cOQ/t24rLhXEhEy+kj32XS464ZKkiTp0hV2yeX6qUOTjnFGmXC5n31A9zOMdweqznWUW5IkSWoPMqF0rwJGnGF8JKnzdUuSJEntWiaU7t8Cw0MII04OhBAGAuOA3yQVSpIkSbpQmVC6f0zqiPbXQgg5IYQs4Kukzl7y3SSDSZIkSRei3ZfuGGMtMA9oAF4B1gE9gDnnuxqlJEmS1B5kxNlLYox7gLuSziFJkiRdinZ/pFuSJEnKdJZuSZIkKc0s3ZIkSVKaWbolSZKkNLN0S5IkSWlm6ZYkSZLSzNItSZIkpZmlW5IkSUozS7ckSZKUZpZuSZIkKc0s3ZIkSVKahRhj0hnSKoRQCWxPOEY/YF/CGToq5zZ9nNv0cF7Tx7lND+c1fZzb9ElqbvfFGG8+04YOX7rbgxDCkhhjadI5OiLnNn2c2/RwXtPHuU0P5zV9nNv0aY9z6/ISSZIkKc0s3ZIkSVKaWbrbxveTDtCBObfp49ymh/OaPs5tejiv6ePcpk+7m1vXdEuSJElp5pFuSZIkKc0s3ZIkSVKaWbolKQOEEMaGEOpCCM8knUVSMkIIt4YQdocQfpx0Fl08S3eahBD+NoTwSghhVQhhWQjhDRe4X20IYcUZbkdDCD9Nc+x271LntWnfbiGErzbttzKEsCOE8GAIoSSNkTNGC16zn2+ay9Nfs/ekOXLGaMnr9hTfAHJaOVpGa8Fr9vUhhPua3gdWhhBeDSH8OIRQnObIGeNS5jaEkBtCeHsI4fEQwtIQwtoQwpIQwrtDCKENYrd7LXjNFoQQvgP8KzAwrSHbsRDCgBDCz0MIG5puvw4hDL3AfXNDCF8KIawPIawJISwMIcxKd+Y/EWP01so34NOkroI0qunreUAdcMsF7LvtDGMFwMEL2b8j31o4r12Al4AvAzlNYyOAvcANSX9vSd9aOLefB/4y6e+hvd5aMrenPMccYCWwA3gm6e+pPdxa+JpdAvwYyG/6eiCwBlh/8v2hM98udW6Ba4BG4O2njL0FiMCXkv6+kr618DX7TlL/8M5rms8fJ/39JDB/eU3vg/eTOgCRDfwX8CrQ7QL2/x6wEejf9PX7gGpgcpt9D0lPYke7Ab2A48AXTxv/I7D2AvZ/wxnG3glsA7KS/v4yeF4/0/RDNZw2PuPkX8DOemuFubV0p2lumx6bBSwHbmp6H3gm6e8r6VsrvGb/Exh62tjfNpWZiUl/f5k6t02l+/kzjC8Ajpz+/tuZbq3wms055b87a+l+f9P3XnLK2CCgAfjEefYdQ+ofhO89bXwt8Me2+h5cXtL6bgYKgadPG58PXBFCGHuunWOMvz/D8PuAH8YYG1slYWZq0bySmsM/xqa/ZSfFGBfGGCtbL2ZGaunc6uxaY27fDeyJMT7e2uEyWEvfZ/86xrjztOEepH6g72+1lJmpJXP7EnDjGcbLga5AbqskzEwtfc3WpytYBnkzsCPGuOXkQIxxN/BK07ZzeSMQOPP83xRC6NaaQc/G0t36Jjbdbz1tfOtp2y9ICGE0qaOxP2xhrkx3yfMaQhgAlACVIYSvhRCWN63hvD+EMD4NWTNNa7xmbw4hPNO0hnNZCOGLIYTC1ouYsVo0t01z+AXgE62cK9O12vtsSLkR+CjwLzHGXa2QL5Nd8tzGlLozbBoNLIox1rZCvkzVqt2gk5rIn88fTWMTLmDfRlJL9E7fNwe4osXpLoClu/X1a7o/etr4kab7vhf5fCeP0Ja3KFXma8m8Dmu6/3zT468CpgD1wCKP5Lb4NVtF6temd8YYrwT+CrgbeCqE0JmPbEHL5/YTwOMxxtWtmirztcr7bAjh74FKUr/i/zdSn/no7Fr1Z1gI4SpSheezLcyV6Vq7G3RG/fjz+YPUHBaGEArOs29VjLHhDPtCG82/pfs8QgivCSHEC7g9c76nuoQ/Oxd4F+3wUqYt1cbz2qXpvgL4coyxLsZ4DPgIqQ9mfOrSv5P2p61fszHGr8cY/yrGuLfp6+Wk5vQa4M4WfTPtTFvObQihCPgQ8E+tkb09S+p9Nsb4bzHGfsDVwF8AT4cQupxnt4yS8M+wbqR+S/uPMcbnLiV/e5XkvOrPtGQO23T+Pf3U+S0Exl3A46qa7vc13XfnT9cGdm+6v5j1grcDJ4BHL2KfTNGW83ryX8YrT13THWPcH0IoA6ZfQI5MkuRr9qSXmu6vAX5+Cfu3V205t/8L+G6MseKiEmamRF+zMcbVIYQPNuX4EKmj3h1FInMbQsgDfkPqNzVfuZB9Mkx7eJ/tbPbxP/N1qu6kjmJXn2ffwhBC9mlHu9t0/i3d5xFjrCJ1GqkLtarpfgSpMw2cNPK07RfifcAPOuIHKNt4Xl8l9enmM/1m52zjGautX7MhhP5n+DDqyTe17IvI0e618dzOAurCn57HdzDQL4SwAqiPMZZeRJZ2qy3nNYQQgC5n+AG9sul+2kXkaPeS+BnWVLh/C7wSY/yHi/izM0bC3aCzWgWcaTnoSOB8S/BWAW8Hivnz+a8H1rVCvvPqUGWjnXiU1L9sbzht/EZSb0DNf0lDCIUhhJ5nepKQOtn7XOBHacqZaS55XpveHOcDE5p+4J58XHdS670XpzF3Jmjpa3Z7COH0cn2yuCxrzaAZqCWv2zExxvExxsknb6TOArGk6esOUbgvUUtes8OBFWd4zhFN9539iGOL3g9OOcL9aozx704Z/8+mJVOdVat0g07ut8DwEMKIkwMhhIGkfuPwm1MfGEIYGEI4teP+jtTZiW447TlvJPXbmDOtFW99bXVuws50I3UC/EqaziUJvIYznAAf2ND0uK5neI5/Bv6Q9PfSnm4tmVdSS0hqgI82fR2A75BaejI26e8t6VsL5zaSWgqR3fT1cFKncFoPFCb9vSV9a433g1Mesw3P092ieSVVriOpz3SEprG+wBOkStH4pL+3pG8tmNs84AFS5z6++7TbOmBE0t9bJs7rGZ6ns56nO4/UEetfkVqpkQX8P067OA4wk9RvW7972v7fa5rbfk1fv5c2vjiOy0vSIMb41RBCDfBQCKGe1P/8t8YYHzntoRWkfq3xJ+ffbDoa+x5SPxTUpCXzGmNcHEK4CfhyCOFvSS17WAPMjKccYeisWviafQdwF7Ci6Yh3IamjOv8UU79l6NRa+n4AEEL4Hqn18acuL3kgxvjPaQ3fjrVgXsuBvwbeBnwwhNBI6of5UmB6jHFtm3wD7VgL5vZmUp9FAvhpm4TNIC15LwipM3MsOuUxdzS9DzwaY/x0epO3DzHG2hDCPODfSR3YiaR+js+JqZMjnHQMOExqHk/1UeBfgBdCCHWkDrrdFGNcke7sJ538V74kSZKkNHFNtyRJkpRmlm5JkiQpzSzdkiRJUppZuiVJkqQ0s3RLkiRJaWbpliRJktLM0i1JkiSlmaVbkiRJSjNLtyRJkpRmlm5JEiGEwhBCdQhhdtJZJKkjsnRLkgBeA1QDC5MOIkkdkaVbkgRwG/BojLEh6SCS1BFZuiUpQ4QQ8kIItSGEfzhtvDiEEEMIb27B098KPHQBGW4LISwOIWwMIawOITweQnhnCOGOEMKKphxfCiF8PYSwPISwK4TwkaZ9PxdCWBJC2B5C+HALskpSxrF0S1LmGA/kAstPG5/cdL/iUp40hDAFGAQ8ep7HvQX4PfCvMcbRwCRgNXBvjPGBGOPJHHcDP48xTgG+ANwbQvgG8FCMsRT4atPY6EvJK0mZyNItSZljStP9ijOMHwW2XOLz3gYsijEeONsDQggB+BbwTIzx9wAxxkbgX4Bjpz18eYxxZdN//6bpftQpY/cDAbj+EvNKUsaxdEtS5pgCbD9DOZ4CrIwxxkt83ts4/9KSMcAwYPGpgzHGYzHGoac9dtMp/33gDGP7m+4HXWROScpYlm5JyhxT+fOlJZAq3Ssu5QlDCAOAUs5fuvs13Z/1aPgpqk7+xyn/EDjTWPYFxpSkjJeTdABJ0vk1Le+YADx22vgoYDinlO4Qwk+AiU1f9iJVmHue5cwktwJlMca154mwr+m+98VmlyR5pFuSMkV3oBtQcdr4J5ruV5wciDG+q+lDjW8ADgNvPsepAC9kaQnABmAHqaPizUIIfUIIL4YQel3Ac0hSp2XplqQMEGM8Qqpw3xlCGBxCKAkh/BvwVqAeWHPq40MI40iV6Q/FGB/7sydMPSYXmMcFlO6mJSH/ANwYQri9af8c4MvA5hjjoUv93iSpM7B0S1LmeBcwlNRZSh5ouv8dsD7GeOLkg0II04BfA3fHGF84x/NdR2qZ4TMX8ofHGH8NvAn45xDCq6TWl9cC7w8hzA4hrGh66AdDCN8JIUy4gLGfX8ifLUmZLlz6h90lSe1NCOF64NuklpRsPM9j/x0oiTG+vk3CSVInZumWpA4khFAJ1PA/p+VrAK6JMdad4bEbgW/EGP9vG0aUpE7J0i1JkiSlmWu6JUmSpDSzdEuSJElpZumWJEmS0szSLUmSJKWZpVuSJElKM0u3JEmSlGaWbkmSJCnNLN2SJElSmv3/W2XiyyRkfn0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "u = np.linalg.solve(A,RHS)\n",
    "plt.plot(u*100,z)\n",
    "plt.xlabel('$u_z$ / cm')\n",
    "plt.ylabel('$z$ / m');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks plausible. Let's go ahead and perform a basic exploration of the numerical properties of the scheme.\n",
    "\n",
    "We first define a simulation run and then vary the number of grid points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_sim(number_of_points):\n",
    "    #discretization\n",
    "    Delta_z = H/(number_of_points - 1)\n",
    "    z = np.linspace(0,H,number_of_points)\n",
    "    #RHS\n",
    "    RHS = np.ones(number_of_points) * gamma/Es * Delta_z**2\n",
    "    #FD Matrix\n",
    "    A = np.diag([-2.]*number_of_points,0) + np.diag([1.]*(number_of_points-1),1) + np.diag([1.]*(number_of_points-1),-1)\n",
    "    #lower bc\n",
    "    A[0,:] = 0.\n",
    "    A[0,0] = 1.\n",
    "    RHS[0] = 0\n",
    "    #upper bc\n",
    "    A[-1,-1] = -1.\n",
    "    RHS[-1] = 0.\n",
    "    #solution\n",
    "    return np.linalg.solve(A,RHS), z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the numerical result against the analytical solution\n",
    "\n",
    "$$\n",
    "u_z(z) = \\frac{\\varrho gH^2}{2E_\\text{s}} \\left[ \\left(\\frac{z}{H}\\right)^2 - 2 \\frac{z}{H} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytical solution\n",
    "reference_z = np.linspace(0,H,100)\n",
    "reference_solution = ((reference_z/H)**2 - 2.*reference_z/H) * gamma * H**2 / (2.*Es)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4cc2f5b655b84aceb2cc53a4cd3f38e3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='nodes', min=5), Output()), _dom_classes=('widget-intera…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Interactive playground\n",
    "\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "#Compute reference solution with 100 cells\n",
    "\n",
    "@interact(num_nodes=widgets.IntSlider(min=5, max=100, value=10, step=1, description='nodes'))\n",
    "def plot(num_nodes=10):\n",
    "    fig,ax = plt.subplots()\n",
    "    ax.set_ylabel('$z$ / m')\n",
    "    ax.set_xlabel('$u_z$ / cm')\n",
    "    solution, z = run_sim(num_nodes)\n",
    "    ax.plot(reference_solution*1e2, reference_z, '--', color='k', label='reference solution')\n",
    "    ax.plot(solution*1e2,z,label='solution')\n",
    "    ax.set_ylim([0,10])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convergence of this simple scheme is rather low. A comparatively simple finite element scheme converges much faster, as we shall see in one of the following exercises.\n",
    "\n",
    "If we plot the convergence rate of the surface settlements with respect to the grid resolution, we observe a linear behaviour:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0sAAAIrCAYAAAA+xwKYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABXr0lEQVR4nO3de5yVZb3//9eHk+OMioeE1PAUoQ6EJ8LC1G+S6RcYFTPJdKcMqfmVPG6Ppbm3/gQ1U5NKSYZMKVEDYwbSNEtN1MQUBBJNTckMPCTSyCiH6/fHGmocFziHNXOvmXk9H4/1WHNfs9a93rNjj7y5rvu6I6WEJEmSJOmDumUdQJIkSZKKkWVJkiRJkvKwLEmSJElSHpYlSZIkScrDsiRJkiRJefTIOoBa7rDDDkv33HNP1jEkSZKkji7yDTqz1IG98cYbWUeQJEmSOi3LkiRJkiTlYVmSJEmSpDwsS5IkSZKUh2VJkiRJkvKwLEmSJElSHpYlSZIkScrDsiRJkiRJeViWJEmSJCkPy5IkSZIk5WFZkiRJkqQ8LEuSJEmSlIdlSZIkSZLysCxJkiRJUh6WJUmSJEnKw7IkSZIkSXlYliRJkiQpD8uSJEmSJOVhWVKLLVoEgwblniVJkqTOxrKkFqmthREjYPFiGDkydyxJkiR1JpYltUhlJSxfDinBsmUwblzWiSRJkqTCsiyp2aqqYPZsqKvLHdfVQXV1blySJEnqLCKllHUGtdCQIUPSvHnz2v1z+/bNzSo11qdPbpZJkiRJ6mAi36AzS2q2CROgrOyDY6WlMHFiNnkkSZKktmBZKhIR0TMi/jsiaiNiUNZ5NqayMrepQ0lJ7rikBCoqYOzYbHNJkiRJhWRZ2oiI2C4i7omI9lireDLwCFDaDp/ValVVuWV3EblleVOmZJ1IkiRJKizL0gZExGjgUeCTH/G6PhExLSKW1D/uiohPNPfzUko/TCk92tK87a2sDObMgfLy3GYPjZflSZIkSR2dZWnDLgAOITfbk1dE9ALuA3oBA4FyoBb4XURs1h4hszRwICxcmHuWJEmSOpseWQcoYvunlNZE5N0YY70TgMHA6JTSGoCIOB94FTgVuLp+7FGg3wbO8dmU0t8KllqSJElSQViWNmB9+fkIXwZeSSm92OB9/4iIxfXfu7p+7HNtk1KSJElSW3EZXusMBl7KM/4S8Ol2ziJJkiSpgCxLrfMxYGWe8XeA0ojYtKkniogDI2JS/eFFETFmA687OSLmRcS8119/vfmJJUmSJDWJy/DaxkYvdMonpfQQ8BAw/iNeNxmYDDBkyJD22NJckiRJ6pKcWWqdN4DN84xvDrybUlrVznkkSZIkFYhlqXUWADvnGd8FeKZ9o0iSJEkqJMtS68wAdoqIndcPRERfYA/gl1mF6mwWLYJBg3LPkiRJUnuxLLXOT8nNIF0ZET0iohswkdxueD/OMlhnUVsLI0bA4sUwcmTuWJIkSWoPlqUNiIirI+Jp4PD646frH73Wvyal9D5wCLAWWAz8GdgCODil9K/2T935VFbC8uWQEixbBuPGZZ1IkiRJXYW74W1ASuncJr5uGfC1No7TJVVVwezZUFeXO66rg+rq3HhlZbbZJEmS1PlFSu4+3VENGTIkzZs3L+sYbaZv39ysUmN9+uRmmSRJkqQCyXvrH5fhdUARURERk1esWJF1lDY1YQKUlX1wrLQUJk7MJo8kSZK6FstSB5RSqk4pndy7d++so7Spysrcpg4lJbnjkhKoqICxY7PNJUmSpK7BsqSiVlWVW3YXkVuWN2VK1okkSZLUVViWVNTKymDOHCgvz2320HhZniRJktRW3A1PRW/gQFi4MOsUkiRJ6mqcWZIkSZKkPCxLkiRJkpSHZUmSJEmS8rAsSZIkSVIeliVJkiRJysOy1AFFREVETF6xYkXWUSRJkqROy7LUAaWUqlNKJ/fu3TvrKJIkSVKnZVmSJEmSpDwsS5IkSZKUh2VJkiRJkvKwLEmSJElSHpYlSZIkScrDsiRJkiRJeViWJEmSJCkPy5IkSZIk5WFZkiRJkqQ8LEsdUERURMTkFStWZB1FDSxaBIMG5Z4lSZLU8UVKKesMaqEhQ4akefPmZR1DQG0tlJfD0qWw4465wlRWlnUqSZIkNVHkG3RmSSqAykpYvhxSgmXLYNy4rBNJkiSptSxLUitVVcHs2VBXlzuuq4Pq6ty4JEmSOi6X4XVgLsMrDn375maVGuvTJzfLJEmSpKLnMjypLUyY8OHrk0pLYeLEbPJIkiSpMCxLUitVVsLIkVBSkjsuKYGKChg7NttckiRJah3LklQAVVW5ZXcRuWV5U6ZknUiSJEmtZVmSCqCsDObMyW0fPnu224ZLkiR1Bj2yDiB1FgMHwsKFWaeQJElSoTizJEmSJEl5WJYkSZIkKQ/LkiRJkiTlYVnqgCKiIiImr1ixIusokiRJUqdlWeqAUkrVKaWTe/funXUUSZIkqdOyLEmSJElSHpYlSZIkScrDsiRJkiRJeViWJEmSJCkPy5IkSZIk5WFZkjqJRYtg0KDcsyRJklrPsiR1ArW1MGIELF4MI0fmjiVJktQ6liWpE6ishOXLISVYtgzGjcs6kSRJUsdnWZI6uKoqmD0b6upyx3V1UF2dG5ckSVLLRUop6wxqoSFDhqR58+ZlHUMZ69s3N6vUWJ8+uVkmSZIkfaTIN+jMktTBTZgAZWUfHCsthYkTs8kjSZLUWViWpA6usjK3qUNJSe64pAQqKmDs2GxzSZIkdXSWJakTqKrKLbuLyC3LmzIl60SSJEkdn2VJ6gTKymDOHCgvz2320HhZniRJkpqvR9YB1HwRUQFU9O/fP+soKiIDB8LChVmnkCRJ6jycWeqAUkrVKaWTe/funXUUSZIkqdOyLEmSJElSHpYlSZIkScrDsiRJkiRJeViWJEmSJCkPy5IkSZIk5WFZkiRJkqQ8LEuSJEmSlIdlSZIkSZLysCxJkiRJUh6WJUmSJEnKw7IkSZIkSXlYliRJkiQpD8uSpDazaBEMGpR7liRJ6mgsS5LaRG0tjBgBixfDyJG5Y0mSpI7EsiSpTVRWwvLlkBIsWwbjxmWdSJIkqXksS5IKrqoKZs+GurrccV0dVFfnxiVJkjqKSCllnUHNFBEVQEX//v1Pev7557OOI31I3765WaXG+vTJzTJJkiQVmcg36MxSB5RSqk4pndy7d++so0h5TZgAZWUfHCsthYkTs8kjSZLUEpYlSQVXWZnb1KGkJHdcUgIVFTB2bLa5JEmSmsOyJKlNVFXllt1F5JblTZmSdSJJkqTmsSxJahNlZTBnDpSX5zZ7aLwsT5Ikqdj1yDqApM5r4EBYuDDrFJIkSS3jzJIkSZIk5WFZkiRJkqQ8LEuSJEmSlIdlSZIkSZLysCxJkiRJUh6WJUmSJEnKw7IkSZIkSXlYliRJkiQpD8uSJEmSJOVhWZIkSZKkPCxLkiRJkpSHZUlSp7doEQwalHuWJElqKsuSpE6tthZGjIDFi2HkyNyxJElSU1iWJHVqlZWwfDmkBMuWwbhxWSeSJEkdhWVJUqdVVQWzZ0NdXe64rg6qq3PjkiRJHyVSSllnUAsNGTIkzZs3L+sYUtHq2zc3q9RYnz65WSZJkqR6kW/QmSVJndaECVBW9sGx0lKYODGbPJIkqWOxLEnqtCorc5s6lJTkjktKoKICxo7NNpckSeoYLEuSOrWqqtyyu4jcsrwpU7JOJEmSOgrLUgcUERURMXnFihVZR5GKXlkZzJkD5eW5zR4aL8uTJEnaEDd46MDc4EGSJEkqCDd4kCRJkqSmsixJkiRJUh6WJUmSJEnKI7OyFBHbZvXZkiRJkvRRemT42csiYhnwNDC/wePZlNK6DHNJkiRJUqZl6UZgH+Ag4FBg/bZ870XEYj5YoOanlN7OIqQkSZKkrimzspRS+n8AEdEdGEiuMF0JvArsSq5IQX2JioilKaWd2z+pJEmSpK4oy5klAFJKa4EFEVFLriydn1KaERHbA4PqH58GyjOMKUmSJKmLybwsbUhK6e/A34HfZJ1FkiRJUtfj1uGSJEmSlIdlSZIkSZLyyPI+S+MiYq+IKNqlgJJUaIsWwaBBuWdJklTcsiwqPyG30937EfEM8FL98U4RUZZSqs0wmyQVXG0tjBgBS5fCyJG5wlRWlnUqSZK0IVkuw/sCcBbw8/ochwMBXAOsiIgXI6I6IiZExHERsWeGWSWp1SorYflySAmWLYNx47JOJEmSNibL+yw9CDy4/rj+fku7AXs2eOwLjFz/FqB7O8eUpIKoqoLZs6GuLndcVwfV1bnxyspss0mSpPwipZR1ho2KiG2BvYDBKaVrMo5TVIYMGZLmzZuXdQxJTdC3b25WqbE+fXKzTJIkKVORb7Dod8NLKb2eUrrPoiSpI5sw4cPXJ5WWwsSJ2eSRJEkfrejLkiR1BpWVuU0dSkpyxyUlUFEBY8dmm0uSJG2YZUmS2klVVW7ZXURuWd6UKVknkiRJG2NZkqR2UlYGc+ZAeXluswe3DZckqbgVzQ1h629OOwyYn1JaERE9gc+tP842nSQVxsCBsHBh1ikkSVJTFNPMUm/gd+S2CwfYutGxJEmSJLWbYipL8OEt+/Ju4SdJkiRJba3YypIkSZIkFQXLkiRJkiTlYVmSJEmSpDwsS5IkSZKUh2VJkiRJkvKwLEmSJElSHpYlSZIkScrDsiRJkiRJefTIOkADK4AvAE/XH7/V6FiSJEmS2k3RlKWU0hrgwQbHqxse6z8iogKo6N+/f9ZRJEmSpE7LZXgdUEqpOqV0cu/evbOOIqmDW7QIBg3KPUuSpA+yLElSF1VbCyNGwOLFMHJk7liSJP1HUZSliCjJOoMkdTWVlbB8OaQEy5bBuHFZJ5IkqbgURVkCVkbEfRHhujJJagdVVTB7NtTV5Y7r6qC6OjcuSZJyiqUsdSe3891vI2Lrxt+MiP0i4jvtH0uSOqcLL/zwsrt3382NS5KknGIpSwAXA72B30XEto2+1x/4n/aPJEmd04QJUFb2wbHSUpg4MZs8kiQVo2IqSy8DBwKbAA9GxMczziNJnVZlZW5Th5L6K0ZLSqCiAsaOzTaXJEnFpJjKEiml14CDgLXAQxHxiYwjSVKnVVUFffpABPTtC1OmZJ1IkqTiUlRlCSCltIxcYfoX8HBE7JJxJEnqlMrKYM4cKC/PbfbQeFmeJEldXdGVJYCU0lvkNnxYDjwE7J5tIknqnAYOhIULc8+SJOmDemQdoF5/cjNJ/5ZSWhERXwR+DXwbSFkEkyRJktQ1FcXMUkrpxZTS8jzjK4EvAXcD3ltekiRJUrsplpmlDUopvQsclXUOSZIkSV1LUcwsSZIkSVKxsSxJkiRJUh4FX4YXEWsLfc5GUkqp6JcPSpIkSerY2qJ0RBucsz3PL0mSJEmFL0spJZf2SZIkSerwLDaSJEmSlIdlSZIkSZLysCxJkiRJUh5FUZbaYQc9SVKRW7QIBg3KPUuSVAyKoizhDneS1KXV1sKIEbB4MYwcmTuWJClrxVKWUtYBJEnZqayE5cshJVi2DMaNyzqRJEkb2To8IoY15QQppbmFiyNJ6mqqqmD2bKiryx3X1UF1dW68sjLbbJKkri1Syj+pExEvNeH9KaW0a6tDRKxNKXVv7Xm6miFDhqR58+ZlHUOSWqVv39ysUmN9+uRmmSRJagd5Lwva4MxSSmmXtssiSVLOhAlw+ukfvE6ptBQmTswukyRJsJGy1FhE9AB2BbahQfNyGZ4kqTUqK+Hee2HWrNwSvJISqKiAsWOzTiZJ6uqaVJYi4vPAdGBTYAvgHWBzYCm5AiVJUotVVUF5OSxdmluWN2VK1okkSWr6bnjXAxNTSlsDK+ufvwtMbrNkkqQuo6wM5szJFabZs3PHkiRlbYMbPHzgRRHvAL1TSiki/plS2ioiegIvp5S2b3UIN3hoETd4kCRJkgoi7wYPTZ1ZegP4WP3Xf4uIfYHtyS3LkyRJkqROp6llaTJwQP3X1wIPAs8AN7ZFKEmSJEnKWpM2eEgpTWzwdVVE/BYoSyktbrNkkiRJkpShJm8d3lBK6eVCB5EkSZKkYtKkZXgRcUFEfKbR2Gci4vy2iSVJkiRJ2WrqNUtnAIsajS0GzixoGkmSJEkqEk0tSz2BtY3G1gK9ChtHkiRJkopDU8vSXOCiRmMXAI8VNo4kSZIkFYembvDwLaA6Ik4EXgZ2Bt4BRrVNLEmSJEnKVlO3Dn85IvYC9gP6AUuBP6aUGi/NkyRJkqROYYPL8CJivwZfDwM+CwTwt/rn/erHJUmSJKnT2djM0hRgUP3X0zbwmgTsWtBEkiRlbNEiGDMGpk+HgQOzTiNJysoGy1JKaVCDr3dpnziSJGWrthZGjIClS2HkyFxxKivLOpUkKQtNvSnt2A2Mn1DYOJIkZauyEpYvh5Rg2TIYNy7rRJKkrDR16/DrNzD+/UIFkSQpa1VVMHs21NXljuvqoLo6Ny5J6no2uhtegw0cukXE58ht7LDerkBtWwWTJKm9XXhhbhleQ+++mxuvrMwmkyQpOx+1dfj6jR02BX7eYDwBy4Az2iKUJElZmDABTj/9g4WptBQmTswukyQpOxstS+s3doiIn6SUTmqfSJIkZaOyEu69F2bNyi3BKymBigoYm/fKXUlSZ9fUa5aejoht2jSJJElFoKoK+vSBCOjbF6ZMyTqRJCkrTS1LhwB/jYjZEfG1iChty1CSJGWlrAzmzIHy8txmD24bLkld10ddswRASunIiNgSOBr4BvCjiKgGbgPuSymta7uInV9EbA1cDbwLrAN2Bv47pfR8lrkkqasaOBAWLsw6hSQpa02dWSKl9HZK6eaU0sHAAcCewK+B1yLiyvq/8HcqEbFdRNwTEamNP2pHoC6l9K2U0hnAfcDNbfyZkiRJkjaiyWUpIkoiYkxEzAIeAhYBo4CDga2A37RNxGxExGjgUeCTH/G6PhExLSKW1D/uiohPNOezUkpPp5ROazD0ErBd81NLkiRJKpQmlaWIuBX4B3AyMBPYMaV0bEppTkppEXAqsFvbxczEBeSu1XpkQy+IiF7kZoF6AQOBcnL3nvpdRGzWis+uAG5sxfslSZIktVKTrlkCFgMXpZSW5vtmSmltROxQuFhFYf+U0pqI2NhrTgAGA6NTSmsAIuJ84FVyBfLq+rFHgX4bOMdnU0p/W38QESOB3sB1rf0BJEmSJLXcBstSRAxrcPgg0C8iPvQX/pTS3PrndwofLzvry89H+DLwSkrpxQbv+0dELK7/3tX1Y59rymfWF6UjgP9y0wxJkiQpWxubWZrWhPcnYNcCZemIBgPP5Rl/CRjenBNFxFeA/YFTUkopIq6v3+xBkiRJUgY2WJZSSru0Z5AO6mPAk3nG3wFKI2LTlNKqjzpJRAwGfgG8AXy1fulfb+BDZSkiTiZ37Rg77rhjy5NLkiRJ2qimXrNERPQgN4u0DfDvC3nWL8PTB2z0QqfGUkoLaPo9ryYDkwGGDBnS1luaS5IkSV1Wk/6CHhGfB6YDmwJbkJs52RxYStdehvcGuf87NLY58G5TZpUkSZIkFaem3mfpemBiSmlrYGX983epn+HowhYAO+cZ3wV4pn2jSJIkSSqkppalTwGTGo1dDZxe2Dgdzgxgp4jYef1ARPQF9gB+mVUoSZIkSa3X1LL0BrnNDAD+FhH7AtuTW5bXlf2U3AzSlRHRIyK6ARPJ7Yb34yyDSZI6vkWLYNCg3LMkqf01tSxNBg6o//pacvddega4sS1CFYOIuDoingYOrz9+uv7Ra/1rUkrvA4cAa8nduPfP5K7pOjil9K/2Ty1J6ixqa2HECFi8GEaOzB1LktpXpNT8DdUiYiegLKW0uCAhItamlLoX4lxdyZAhQ9K8efOyjiFJagNjxsCsWVBXByUlcMQRcPvtWaeSpE4r727WTZ1Z+oCU0suFKkqSJOmDqqpg9uxcUYLcc3V1blyS1H6aPLMUEcOBPYHNGo6nlP631SGcWWqWiKgAKvr373/S888/n3UcSVKB9e0Ly5d/eLxPH1i2rP3zSFIX0PKZpYj4EXArsC/Qr8HjE4VKp6ZLKVWnlE7u3bt31lEkSW1gwgQoK/vgWGkpTJyYTR5J6qqadFNa4KvAwJTSa20ZJiKGNeV1KaW5bZlDkqQsVVbCvfd+8JqligoYOzbrZJLUtTS1LL0GvNeWQepNa8JrErBrWweRJClLVVVQXg5Ll+aW5U2ZknUiSep6mlqWfgLMjIjrgA+sli7kLE9KaZdCnUuSpI6srAzmzMntijd9+oeX5UmS2l6TNniIiJc28K2UUmr1LI8bPLSMW4dLkiRJBZF3g4cmzSxlNeMTET3ILbnbhgY/gNcsSZIkSWprTV2G1+7FJSI+D0wHNgW2AN4BNgeW4jVLkiRJktpYk8pSc4tLRHwJeCCltKYV2a4HJqaUboiIf6aUto6Ii4B1rTinJEmSJDVJk+6zxH+Ky9bAyvrn7wKTN/D6e4BfRcQmDQcj4sBmZPsUMKnR2NXA6c04hyRJkiS1SFPLUnOLSy0wH5gTEaUNxmuake0N4GP1X/8tIvYFtic3u9WlRURFRExesWJF1lEkSZKkTqupZam5xSWllC4CHgbujYj1G57m3WViAyYDB9R/fS3wIPAMcFMzztEppZSqU0on9+7dO+sokiRJUqfV1A0e1heXGfynuKwDfriB1wdASunSiFgFPBARh5K7oWyTpJQmNvi6KiJ+C5SllBY39RySJEmS1FJN3Tq8ucXlvgavvzIi3gN+19TPA4iIC4DfppSeqD/PyxHxmYg4P6V0ZVPPI0mSJEkt0dRleB+QUnp5YzM8KaWjGh1fB9wIvNeMjzkDWNRobDFwZjPOIUmSJEkt0uSZnoYi4o/krh9asP45pfTGxt6TUrqJ5l1v1BNY22hsLdCrGeeQJEmSpBZpUVkCXgcOBcZSfx1SRCzjwwXq6VZkmwtcRG6L8vUuAB5rxTklSZIkqUlaVJZSSiMBImIbYG9gL2AouQL1xQYv7d6KbN8CqiPiROBlYGdyN8Md1YpzSpKkFli0CMaMgenTYeDArNNIUvto6cwSACmlN4H76x/UbxH+I3LF5tpWnvvliNgL2A/oBywF/phSarw0T5IktaHaWhgxApYuhZEjc8WprOyj3ydJHV2LNnjYkJRSLXAiueuKtm/u+yNivwZfDwM+S24b8r/VP+9XPy5JktpJZSUsXw4pwbJlMG5c1okkqX20amYpn5RSioiZwGnkZpmaYwowqP7raRv6CGDXFsaTJEnNUFUFs2dDXV3uuK4Oqqtz45WV2WaTpLYWKTX5PrH/eVPEFOBpcps5zE8pvd3o+5cDZ6aUNmvi+damlFpzfVOXEhEVQEX//v1Pev7557OOI0nqxPr2zc0qNdanT26WSZI6icg32NJleIcC1wMPAG9GxF8j4u6IuCIirgfOIrcjXotFxNgNjJ/QmvN2Biml6pTSyb179846iiSpk5sw4cPXJ5WWwsSJ+V8vSZ1Ji2aWACJiK2DPRo/dgU2A+cCJKaUmFaZ8M0sR8U5KaYs8r30zpbRNi0J3MkOGDEnz5s3LOoYkqZMbMwZmzcotwSspgSOOgNtvzzqVJBVU3pmlFl+zlFL6J/D7+sd/PiWie2t2rGuwgUO3iPgcHwy+K1Db0nNLkqTmq6qC8vLcbnh9+8KUKVknkqT20RYbPLR2a+/1GztsCvy84amBZcAZrTy/JElqhrIymDPnP/dZcttwSV1FwctSa6WUdgGIiJ+klE7KOo8kScrdiHbhwqxTSFL7Kuh9lgrs6Yjw2iRJkiRJmSjmsnQI8NeImB0RX4uI0qwDSZIkSeo6irYspZSOBPoBM4FvAH+PiFsj4tCIKNrckiRJkjqHgpSOiDg4IrYrxLkaSim9nVK6OaV0MHAAue3Jfw28FhFXRsTWhf5MSZIkSYLCzSzdB1QU6Fz/FhElETEmImYBDwGLgFHAwcBWwG8K/ZmSJEmSBIXbDS/vTZxadcKIW8kVsCeB24DjUkorG3z/VODtQn+uJEmSJEERbh3ewGLgopTS0nzfTCmtjYgd2jmTJEmSpC6iqMpSRAxrcPgg0C8i+jV+XUppbv3zO+2VrZhERAVQ0b9//6yjSJIkSZ1WUZUlYFoTXpOAXds6SDFLKVUD1UOGDPGmvZIkSVIbKaqylFLaJesMkiRJkgRFVpYai4ge5GaRtqHBJhLrl+FJkiRJUlsp2pu7RsTngZeBx4CHgRpy1zHdlmUuSZKUvUWLYNCg3LMktZWiLUvA9cDElNLWwMr65+8Ck7ONJUmSslRbCyNGwOLFMHJk7liS2kIxl6VPAZMajV0NnJ5BFkmSVCQqK2H5ckgJli2DceOyTiSpsypUWfoCUF2gc633BvCx+q//FhH7AtsDmxb4cyRJUgdRVQWzZ0NdXe64rg6qq3PjklRoBSlLKaUHU0qvFeJcDUwGDqj/+lpy1ys9A9xY4M+RJEkdxIUXfnjZ3bvv5sYlqdCKdhleSmliSmlG/ddVwEDgsyklfx1KktRFTZgAZWUfHCsthYkTs8kjqXMr2rLUWErp5ZTS4qxzSJKk7FRW5jZ1KCnJHZeUQEUFjB2bbS5JnVOx32dpOLAnsFnD8ZTS/2aTSJIkZa2qCsrLYelS6NsXpkzJOpGkzqpoZ5Yi4kfArcC+QL8Gj09kmUuSJGWrrAzmzMkVptmzP7wsT5IKZYMzSxExrCknSCnNLVycD/gqMLANNo6QJEkd3MCBsHBh1ikkdXYbW4Y3rQnvT8CuBcrS2GvAe210bkmS2tyECRPo0aMH5557btZRJEktsMGylFLapT2D5PETYGZEXAcsa/iNNpzNkiSpIObOnctFF10EwB/+8Aduvvlmtt1224xTSZKaI1JKTXthRA9ys0jbALF+vBDFJSLWppS6Nxp7aQMvTymltprN6hAiogKo6N+//0nPP/981nEkSXmsXLmSAw44gPnz5wOwxRZbcMUVV3DKKafQo0dR768kSV1R5B1sSlmKiM8D04FNgS2Ad4DNgaWFKC75ypI+2pAhQ9K8efOyjiFJ2oDVq1dz4IEH8thjj/17bPDgwUyaNIkDDjhgI++UJLWzvGWpqbvhXQ9MTCltDaysf/4uMLlA4fKKiB4RMSAiPhcRw9Y/2vIzJUkqlJ49ezJnzhx22203ADbddFMWLFjAgQceyHHHHcff//73jBNKkjamqWXpU8CkRmNXA6cXNs5/1M9mvQw8BjwM1AAPArdt4PVfql8qKElS0dhqq62oqalh6623ZtWqVXz2s5+lpKSEn//85+y2225cffXVvP/++1nHlCTl0dSy9Abwsfqv/xYR+wLbk1uW11aaO5t1D/CriNik4WBEHNiGGSVJ+kj9+/dnxowZ9OzZk8cee4yLL76Y0aNH869//YvzzjuPwYMHc++992YdU5LUSFPL0mRg/eLqa8nN8DwD3NQWoeo1dzarFpgPzImI0gbjNW2QTZKkZjnooIO46abcfzYvueQS/t//+3/cc889DBgwgCVLlnDYYYcxevRoXnppQ/sbSZLaW5PKUkppYkppRv3XVcBA4LMppQvaMFtzZ7NSSukickv27o2I9ffzznuxliRJ7W3s2LGcf/75rF27lqOPPpqddtqJZ555hquuuorNNtuMu+++m/Lycv7nf/6HVatWZR1Xkrq8JpWliLggIj6z/jil9DJQFhHnt1myDc9m3bihmPXZLiU3m/RARGxJ7sa5kiQVhSuuuILRo0ezYsUKRo0axcqVKzn33HNZsmQJxx13HHV1dVx66aWUl5czc+ZMmnqLD0lS4TV1Gd4ZwKJGY4uBMwuapoGNzGZduIG33NfgvVcCvwB+x0ZuvCtJUnvr1q0bt956K/vssw8vvPACRx11FO+99x7bb789t912Gw899BCDBw/mr3/9K0cddRSHHXYYS5YsyTq2NmDRIhg0KPcsqfNpalnqCaxtNLYW6FXYOBuWUno5pbR4I98/qtHxdeRmod5r42iSJDVLWVkZs2bNYvvtt+ehhx7ilFNO+fcM0gEHHMCTTz7JpEmT2HLLLfnNb37Dpz/9ac477zxWrlyZcXI1VFsLI0bA4sUwcmTuWFLn0tSyNBe4qNHYBeS29W4XEfHHiJgSEWdExMER8bGPek9K6aaU0lbtkU+SpObYYYcdqK6uprS0lFtuuYWrrrrq39/r0aMHp512Gs899xwnnXQSa9as4eqrr2a33XZj2rRpLs0rEpWVsHw5pATLlsG4cVknklRo0ZRfuBGxE1AN9CZ376OdgJXAqJTSX1sdImJtSqn7R7xmNrAnuU0e1odeRu46pgXrn1NKT7c2T0cxZMiQNG/evKxjSJJaYebMmXz5y18mpcQvf/lLjjrqqA+95oknnuBb3/oWjz/+OJCbfbrhhhvYc8892zuu6lVVwemnf3A2qbQUbrghV6IkdTh5N4VrUlkCiIhuwH5AP2Ap8MeUUuOleS1L1oSy1OC12wB7A3sBQ4FDgc3Wf7+p5+kMLEuS1DlcddVVnH/++Wy66aY8/PDD7Lvvvh96zbp167jllls4//zzef311+nWrRunnnoql112GVtt5SKK9ta3b25WqbE+fXKzTJI6nLxlqanL8EgprUspPZpSuqP+uSBFqblSSm+mlO5PKX0vpXQMuZmm24A/AF/OIpMkSa1x7rnnMnbsWFatWsXhhx/Oq6+++qHXdOvWjbFjx/Lcc89xxhlnEBH88Ic/ZMCAAfzkJz9h7dpM/rPcZU2YAGVlHxwrLYWJE7PJI6ltbHBmKSL2Syk9Xv/1sA2dIKU0t9UhmjGztIH3B7nrqm5NKf2otXk6CmeWJKnzeP/99znkkEN46KGH2HvvvXn44Ycpa/y38QaeeeYZTj/9dH7/+98DMGTIECZNmsR+++3XTok1ZgzMmgV1dVBSAkccAbffnnUqSS3UvGV4EbEwpTSo/usN3U48pZR2bXWyVpal+nOcB5yQUhrY2jwdhWVJkjqXN998k/32248XXniBI488kl/+8pd067bhRSApJe644w7OOeecf89GjR07lokTJ9KnT5/2it1l1dZCeTksXQo77pjbPnwj/VZScWveMrwGRak7UArskVLapdGj1UWpqep3wvtWRBxUf7PZxrYgt/GEJEkd0jbbbENNTQ1bbrkld999Nxdd1Hgj2g+KCMaMGcOzzz7LhRdeSM+ePZk6dSoDBgzg+uuvZ82aNe2UvGsqK4M5c3KFafZsi5LUGTV1N7yFwMEppTyXMhYgRNN2w/sbH9wJbynwNLmb45YB3yC3G97n2iJjMXJmSZI6p9/+9rccdthhrFmzhqqqKsaOHduk9z3//POcccYZ/PrXvwZg0KBB3HDDDfyf//N/2jCtJHUKrdrgYSpQHRFjImL/iBi2/lG4fBuXUvoEsA0wHDgbeIDcznxnAKcBzwInt1ceSZLayvDhw/nRj3KX4J5yyik8+OCDTXrfpz71KWbPns2sWbPYddddWbhwIV/4whf46le/yt/+9re2jCxJnVJTZ5aK+pqliOie1e58WYiICqCif//+Jz3//PNZx5EktZFzzjmH73//+2y99dY89thjfOpTn2rye+vq6vje977HFVdcwapVqygtLeU73/kOZ599NptsskkbppakDql191lqS4XY4KErchmeJHVua9eu5cgjj6SmpoYBAwbw2GOPNfueSi+//DL//d//zV133QVA//79uf766xkxYkRbRJakjqp191mSJEntq3v37vz85z9n8ODBPPfccxx99NGsXr26WefYaaeduPPOO7n//vvZY489+Mtf/sLIkSM5/PDDeeGFF9oouSR1Dk0uSxGxW0RcGBE/qD/ePSIGt100SZK0+eabU11dTd++fXnggQcYP348LVkVMnz4cObPn88111zz73MOHDiQiy++mHfffbcNkktSx9ekshQRY4AHgR2AE+uHy4Dr2iSVJEn6tx133JFZs2ZRUlLC5MmTue6661p0np49e3L22WezZMkSvv71r/Pee+9x+eWXs8cee3DXXXe1qIRJUmfW1Jml/wW+mFIaD6zfSGEB4MySJEntYOjQodxyyy1AbuOH6urqFp9ru+2245ZbbuGRRx5h77335pVXXuErX/kKhxxyCIsXLy5UZEnq8JpalrYmdz8j+M99jgLwbneSJLWTY445hssuu4yUEsceeyzz589v1fmGDRvGE088wY9//GO23nprfvvb37Lnnntyzjnn8M477xQotSR1XE0tS48ApzcaOwl4uLBxJEnSxnz729/muOOOo7a2loqKCl577bVWna979+5885vf5LnnnuOb3/wma9eu5fvf/z4DBgzgZz/7GevWrStQcknqeJp6n6UdgGqgFNgV+DOwDhiVUnq11SHcOrxF3Dpckrqmuro6hg8fzty5c/nMZz7Dgw8+yKabblqQc//pT39i/PjxPProo0Bu9umGG25gn332Kcj5JalItXzr8PpCtC9wAnAc8E1gSCGKkiRJap6SkhJmzpzJzjvvzBNPPMEJJ5xQsBmgffbZhz/84Q/ccsst9O3bl7lz5zJkyBBOPfVU3nzzzYJ8hiR1FM25z9KXgIOATwP/F/huRPxvm6SSJEkb1adPH2pqathiiy248847ufTSSwt27m7duvH1r3+dJUuWcPbZZ9O9e3duvPFGBgwYwI033sjatWs/+iRqd4sWwaBBuWdJhdHUrcNvBn4ClAPbNnpIkqQMDBw4kDvuuINu3bpx2WWXcdtttxX0/L179+aaa65h/vz5DB8+nLfeeotTTz2Vz3zmM8ydO7egn6XWqa2FESNg8WIYOTJ3LKn1mjqz9GVg75TSiSmlUxs+2jKcJEnauEMPPZQf/OAHAIwbN45HHnmk4J9RXl7Offfdx5133km/fv146qmn2H///TnhhBP4xz/+UfDPU/NVVsLy5ZASLFsG48ZlnUjqHJpalpYAhblyVJIkFdRpp53G+PHjef/99znyyCN58cUXC/4ZEcHRRx/Nn//8Z77zne+wySab8LOf/YwBAwbw/e9/n9WrVxf8M9U0VVUwezbU1eWO6+qgujo3Lql1mrob3q7AjcBvgOUNv5dS+lmrQ7gbXou4G54kab01a9YwatQo7r33XsrLy5k7dy69e/dus8974YUXOOuss/59c9w99tiDG264geHDh7fZZyq/vn1zs0qN9emTm2WS1CQt3w0P+AbweeBI4L8aPI4vRDJJktQ6PXr0YPr06ZSXl7N48WLGjBnDmjVtd+/4T37yk8yaNYvZs2fTv39//vznP/PFL36Rr3zlK7zyyitt9rn6sAkToKzsg2OlpTBxYjZ5pM6kqWVpPPDplNLnU0qHNHh8qS3DSZKkpuvduzc1NTVsu+223HvvvZx11llt/pkjRoxg4cKFXHHFFZSWlnLXXXex++67c/nll1O3fl2Y2lRlZW5Th5KS3HFJCVRUwNix2eaSOoOmlqWlgDdXkCSpyO2yyy7cfffd9OrVi0mTJjFp0qQ2/8xNNtmECy+8kGeffZYxY8awatUqLr74YgYOHEh1dTVNWfKv1qmqyi27i8gty5syJetEUufQ1LI0BfhVRIyOiGENH20ZTpIkNd+wYcOoqr+6/4wzzuCee+5pl8/t168ft99+Ow888AADBw7kxRdf5PDDD2fUqFE8//zz7ZKhqyorgzlzoLw8t9lD42V5klqmqRs8vLSBb6WU0q6tDuEGDy3iBg+SpI255JJLuOyyy9h8882ZO3cugwYNarfPXr16NT/60Y+45JJLeOedd+jVqxfnnHMOF110EZtttlm75ZCkJsq7wUOTylJbsyy1jGVJkrQx69at49hjj+WOO+5g55135vHHH6dPnz7tmmHZsmVceOGFTJ06FYAddtiBa665hmOOOYaIvH83kaQstGo3PEmS1MF069aNn/70pwwdOpS//vWvjB49ut03Xejbty9VVVU8+uij7Lvvvrz66qt89atf5eCDD2bhwoXtmkWSmqupy/C+vqHveZ+l7DizJElqin/84x8MHTqUpUuXctxxx3HrrbdmMquzdu1aqqqquPDCC3nzzTfp3r0748eP59JLL2XLLbds9zyS1ECrZpb+q9HjXOBmYIMlSpIkFYePf/zj1NTUsNlmmzFt2jQuv/zyTHJ0796dk046ieeee47TTjuNlBLXX389AwYMYOrUqaxbty6TXJK0IS2+ZikivgZ8NqV0eqtDOLPULBFRAVT079//JHcXkiQ1VU1NDYcffjgpJW6//XbGjBmTaZ758+czfvx4/vCHPwCw3377MWnSJIYMGZJpLkldUsGvWbodZ5YykVKqTimd3Lt376yjSJI6kFGjRnHNNdcAcOKJJ/L4449nmmfPPffkoYceYtq0aWy33XY8/vjjDB06lJNOOonXX38902ySBE0sSxGxfaNHf+A7wGttG0+SJBXSmWeeycknn0xdXR1HHHEEr7zySqZ5IoKvfe1rLFmyhHPPPZfu3btz8803M2DAAH74wx+yZs2aTPNJ6tqaOrP0N2Bp/fPfgKeBQ4ET2iaWJElqCxHBpEmTGD58OMuWLaOiooKVK1dmHYvNN9+cq666imeeeYYvfelLvP3224wfP559992Xhx9+OOt4krqoJpWllFK3lFL3+uduKaXNUkr7p5T+2NYBJUlSYfXs2ZM777yT3XbbjQULFnDssceydu3arGMBsPvuu3PPPfcwc+ZMdtppJxYsWMCBBx7Icccdx9///ves40nqYpq8wUNEDAf2BD5w2+2U0v+2OoQbPLSIW4dLklrjL3/5C/vttx9vvfUWZ511Ft///vezjvQBq1at4sorr+TKK6+krq6OzTbbjEsuuYQzzjiDXr16ZR1PUufS8g0eIuJHwK3AvkC/Ro98r/9SRPRoWU5JktQe+vfvz4wZM+jZsyfXXnstN910U9aRPmDTTTfl0ksvZfHixRx55JH861//4rzzzmPw4MHce++9WceT1AU09ZqlrwL7ppSOSymd1PCxgdffA/wqIjZpOBgRB7YmrCRJKqyDDjro3yXptNNO4/7778840YftsssuzJw5k3vuuYcBAwawZMkSDjvsMEaPHs1LL72UdTzlsWgRDBqUe5Y6sqaWpdeA95px3lpgPjAnIkobjNc04xySJKkdjB07lvPPP5+1a9dy9NFH8+yzz2YdKa9DDz2UZ555hiuvvJKysjLuvvtuysvL+Z//+R9WrVqVdTzVq62FESNg8WIYOTJ3LHVUTS1LPwFmRsToiBjW8LGB16eU0kXAw8C9EVFWP553LaAkScrWFVdcwejRo1mxYgWjRo3izTffzDpSXr169eK8885jyZIlfO1rX6Ouro5LL72U8vJyZs6cSVOvxVbbqayE5cshJVi2DMaNyzqR1HJN2uAhIjY0x51SSrvmef3KlNLm9V+fDxxFbqvxV1JKW+R5vRs8tIAbPEiSCqm2tpYDDzyQP/3pTxx44IH85je/YZNNNvnoN2booYce4lvf+hYLFiwA4Etf+hI/+MEP2G233TJO1jVVVcHpp39wNqm0FG64IVeipCKWd1KnybvhNeuTImaklI5qcHwmuXsy7ZZSKs3zestSC1iWJEmF9uqrrzJ06FD+/ve/c8IJJzB16lQiinthyJo1a7jxxhu5+OKLefvtt+nZsydnnnkmF198MZtvvnnW8bqUvn1zs0qN9emTm2WSiljLd8NrroZFqf74OuBGmnfdkyRJamc77LAD1dXVlJaWcsstt3DVVVdlHekj9ejRg/Hjx/Pcc8/xjW98gzVr1nD11Vez2267MW3aNJfmtaMJE6Cs7INjpaUwcWI2eaTWatHMUkT8EXgGWLD+OaX0RotDOLPUIs4sSZLaysyZMznqqNy/ff7yl7/899cdwRNPPMH48eP54x//CMABBxzADTfcwJ577plxsq5hzBiYNQvq6qCkBI44Am6/PetU0kcq6MzS6+SuQboWuA9YFhF/j4h7I+LqiPh6ROzVwnNLkqSMjR49mon10wHHH388Tz75ZMaJmu4zn/kMjz76KFVVVWy77bY8/PDD7LPPPowfP55//vOfWcfr9KqqcsvuInLL8qZMyTqR1HKtumYpIrYB9gb2AoaSK1Cbrf9+U2eLnFlqGWeWJEltKaXEuHHjmDp1Kttvvz1//OMf2WGHHbKO1Sxvv/023/3ud/nhD3/I2rVr+djHPsYVV1xBZWUl3bv7V4+2smhRboZp+nQYODDrNFKTtP0GD/VbhP8I2Bm4NqV0dxPfZ1lqAcuSJKmtvf/++xxyyCE89NBD7L333jz88MOUNb4opQN45pln+Na3vsWDDz4IwJAhQ5g0aRL77bdfxskkFYm23+AhpVQLnAj0ArYv5LklSVL769WrFzNmzOCTn/wkTz31FMcffzzr1q3LOlazffrTn+Z3v/sdt99+OzvssAPz5s3js5/9LJWVlSzPt32bJNEGu+Gl3FTVTOC0Qp9bkiS1v2222Yaamhq23HJL7r77bi666KKsI7VIRDBmzBieffZZLrjgAnr27MnUqVMZMGAA119/PWvWrMk6oqQi06KyFBFTIuJbEXFQRGyZ5yVbADu1KpkkSSoau+++O3fddRc9evTgyiuvZOrUqVlHarHNNtuMCRMmsHDhQv7v//2/rFixgjPPPJO9996b3//+91nHk1REWrp1+N/ILbNb/+alwNPAYqAM+Aa57cQ/18Tzec1SC3jNkiSpvf3kJz/h5JNPpmfPntx3330cdNBBWUdqlZQSNTU1nHHGGbz00ksAjBkzhu9973t84hOfyDidpHZUuGuWUkqfALYBhgNnAw8A/YAzyC2/exY4uUUxJUlS0TrppJM4++yzWb16NUcddRTPP/981pFaJSKoqKhg8eLF/O///i+bbrop06dPZ7fddmPChAm89957WUeUlKGC7oYHEBHdU0prm/medSmlgl8/1dk5syRJysLatWs58sgjqampYcCAATz22GNstdVWWccqiJdffplzzjmHX/7ylwD079+f66+/nhEjRmScTFIba/vd8ACaW5Tq32NRkiSpg+jevTs///nPGTx4MM899xxHH300q1evzjpWQey0007cdddd3Hfffey+++785S9/YeTIkRx++OG88MILWceT1M4sKZIkqdk233xzqqur6du3Lw888ADjx4+n0KtVsvTFL36R+fPn873vfe/fP+vAgQO5+OKLeffdd7OOJ6mdWJYkSVKL7LjjjsyaNYuSkhImT57Mddddl3WkgurVqxfnnHMOS5Ys4etf/zrvvfcel19+OXvssQd33XVXpyqHkvIreFmKiLVt/PAmCJIkFYmhQ4dyyy23AHDOOedQXV2dcaLC22677bjlllv4wx/+wF577cUrr7zCV77yFQ455BAWL16cdTxJbagtZpaijR/OhkmSVESOOeYYLrvsMlJKHHvsscyfPz/rSG1i//33Z968efz4xz9m66235re//S177rkn55xzDu+8807W8ZTHokUwaFDuWWqJgu+Gp/bjbniSpGKRUuK//uu/mDZtGv369ePxxx9nu+22yzpWm3nzzTf5zne+w0033URKib59+3LVVVdx/PHH062b/65bDGprobwcli6FHXfMFaaysqxTqYi1z254kiSp64kIbr75ZoYNG8bSpUs54ogjWLVqVdax2sw222zDj3/8Y+bNm8fnPvc5li1bxgknnMABBxzAn/70p6zjCaishOXLISVYtgzGjcs6kToiy5IkSSqIkpISZs6cyc4778wTTzzBCSecwLp167KO1ab22Wcf/vCHP3DLLbfQt29f5s6dy5AhQzj11FN58803s47XZVVVwezZUFeXO66rg+rq3LjUHJYlSZJUMH369KGmpoYtttiCO++8k0svvTTrSG2uW7dufP3rX2fJkiWcffbZdOvWjRtvvJEBAwZw4403snZts29BqVa68MLcMryG3n03Ny41h2VJkiQV1MCBA7njjjvo1q0bl112GbfddlvWkdpF7969ueaaa5g/fz4HH3wwb731Fqeeeiqf+cxnmDt3btbxupQJEz58fVJpKUycmE0edVyWJUmSVHCHHnooP/jBDwAYN24cjzzySMaJ2s/AgQO5//77ufPOO+nXrx9PPfUU+++/PyeccAL/+Mc/so7XJVRWwsiRUFKSOy4pgYoKGDs221zqeCxLHVBEVETE5BUrVmQdRZKkDTrttNMYP34877//PkceeSQvvvhi1pHaTURw9NFH8+c//5nvfOc79OrVi5/97GcMGDCA73//+6xevTrriJ1eVRX06QMR0LcvTJmSdSJ1RG4d3oG5dbgkqditWbOGUaNGce+991JeXs7cuXPp3bt31rHa3V/+8hfOOussampqANhjjz34wQ9+wBe/+MWMk3VuixbBmDEwfToMHJh1GhU5tw6XJEntq0ePHkyfPp3y8nIWL17MmDFjWLNmTdax2l3//v2prq6mpqaG/v378+c//5lDDjmEr3zlK7zyyitZx+u0Bg6EhQstSmo5y5IkSWpTvXv3pqamhm233ZZ7772Xs846K+tImRk5ciQLFy7kiiuuoLS0lLvuuovdd9+dyy+/nLr1+1xLKhqWJUmS1OZ22WUX7r77bnr16sWkSZOYNGlS1pEys8kmm3DhhRfy7LPPcswxx7Bq1SouvvhiBg4cSHV1NW19icRrr73GYYcdRkTeVUeSGrAsSZKkdjFs2DCq6u8KesYZZ3DPPfdknChb/fr1Y/r06TzwwAMMHDiQF198kcMPP5xRo0bx/PPPt8lnzpw5k8997nO88MILbXJ+qbOxLEmSpHZz3HHHcfHFF7Nu3TqOOeYYFi5cmHWkzH3hC1/gqaee4rrrrmOLLbZgzpw5DBo0iIsuuoh//etfBf2siRMnct9997H//vsX9LxSZ2VZkiRJ7erSSy/lmGOOYeXKlVRUVLB8+fKsI2WuZ8+enHHGGTz33HOMHTuW999/nwkTJrD77rszffr0gi3Ne+SRR/jUpz5VkHNJXYFlSZIktatu3brx05/+lKFDh/LXv/6V0aNHu7lBvb59+1JVVcWjjz7Kvvvuy6uvvspXv/pVDj744ILMwvXo0aMAKaWuw7IkSZLa3aabbsqvfvUr+vXrx9y5c/nGN77R5hsbdCSf/exnefzxx5k8eTLbbLMNv//979lrr70488wzefvtt7OOJ3UZliVJkpSJj3/849TU1LDZZpsxbdo0Lr/88qwjFZXu3btz0kkn8dxzz3HaaaeRUuL6669nwIABTJ06lXXr1mUdUer0LEuSJCkzgwcP5he/+AURwSWXXML06dOzjlR0tt56ayZNmsSTTz7J5z//eV5//XUqKysZNmwY8+bNyzqe1KlZliRJUqZGjRrFNddcA8CJJ57I448/nnGi4rTXXnvx0EMPcdttt7Hddtvx+OOPM3ToUE466SRef/31rONJnZJlSZIkZe7MM8/k5JNPpq6ujiOOOIJXXnkl60hFKSI47rjjWLJkCeeeey7du3fn5ptvZsCAAfzwhz9kzZo1WUeUOhXLkiRJylxEMGnSJIYPH86yZcuoqKhg5cqVWccqWptvvjlXXXUVzzzzDIcccghvv/0248ePZ9999+Xhhx/OOp7yWLQIBg3KPavjsCxJkqSi0LNnT+6880522203FixYwLHHHsvatWuzjlXUdt99d+69915mzJjBTjvtxIIFCzjwwAP54Q9/mHU0NVBbCyNGwOLFMHJk7lgdg2VJkiQVja222oqamhq23nprZs+ezbnnnpt1pKIXEYwePZrFixfz3e9+l2222YYjjjgi72vPPfdc9tprL2bNmgXkroPaa6+9eP/999szcpdTWQnLl0NKsGwZjBuXdSI1VXhPg45ryJAhyV1wJEmd0YMPPsghhxzC6tWrufHGGznllFOyjtRh1NbWUlZWlnUM1auqgtNP/+BsUmkp3HBDrkSpaETeQctSx2VZkiR1ZlOnTqWyspLu3btzzz338MUvfjHrSFKz9e2bm1VqrE+f3CyTikbesuQyPEmSVJTGjh3L+eefz9q1azn66KN59tlns44kNduECdB4oq+0FCZOzCaPmseyJEmSitYVV1zB6NGjWbFiBaNGjeLNN9/MOpLULJWVuU0dSkpyxyUlUFEBY8dmm0tNY1mSJElFq1u3btx6663ss88+vPDCCxx11FG89957WceSmqWqKrfsLiK3LG/KlKwTqaksS5IkqaiVlZUxa9Ystt9+ex566CFOOeUUvOZaHUlZGcyZA+XlMHv2h5flqXhZliRJUtHbYYcdqK6uprS0lFtuuYWrrroq60hSswwcCAsX5p7VcViWJElSh7DPPvtw2223AXDBBRcwY8aMjBNJ6uwsS5IkqcMYPXo0E+u3ETv++ON58sknM04kqTOzLEmSpA7lvPPOY+zYsaxatYrDDz+cV199NetIkjopy5IkSepQIoIbb7yRAw88kL///e9UVFRQW1ubdSxJnZBlSZIkdTi9evVixowZfPKTn+Spp57i+OOPZ926dVnHktTJWJYkSVKHtM0221BTU8OWW27J3XffzUUXXZR1JEmdjGVJkiR1WLvvvjt33XUXPXr04Morr2Tq1KlZR5LUiViWJElShzZ8+HB+9KMfAXDKKafw4IMPZpxIUmdhWZIkSR3eSSedxNlnn83q1as56qijeP7557OOJKkTsCxJkqRO4aqrrmLUqFG89dZbjBo1in/+859ZR5LUwVmWJElSp9C9e3d+/vOfM3jwYJ577jmOPvpoVq9enXUsqegsWgSDBuWetXGWJUmS1GlsvvnmVFdX07dvXx544AHGjx9PSinrWFLRqK2FESNg8WIYOTJ3rA2zLEmSpE5lxx13ZNasWZSUlDB58mSuu+66rCNJRaOyEpYvh5Rg2TIYNy7rRMXNsiRJkjqdoUOHcssttwBwzjnnUF1dnXEiKXtVVTB7NtTV5Y7r6qC6Ojeu/CxLkiSpUzrmmGO47LLLSClx7LHHMn/+/KwjSZm68MIPL7t7993cuPKzLEmSpE7r29/+Nscddxy1tbVUVFTw2muvZR1JysyECVBW9sGx0lKYODGbPB2BZUmSJHVaEcHNN9/MsGHDWLp0KUcccQSrVq3KOpaUicrK3KYOJSW545ISqKiAsWOzzVXMLEuSJKlTKykpYebMmey888488cQTnHDCCaxbty7rWFImqqqgTx+IgL59YcqUrBMVN8uSJEnq9Pr06UNNTQ1bbLEFd955J5deemnWkaRMlJXBnDlQXp7b7KHxsjx9kGVJkiR1CQMHDuSOO+6gW7duXHbZZUybNi3rSFImBg6EhQtzz9o4y5IkSeoyDj30UH7wgx8AUFlZySOPPJJxIknFzLIkSZK6lNNOO43x48fz/vvvc+SRR/Liiy9mHUlSkbIsSZKkLufaa6/l0EMP5Y033qCiooIVK1ZkHUlSEbIsSZKkLqdHjx5Mnz6d8vJyFi9ezJgxY1izZk3WsSQVGcuSJEnqknr37k1NTQ3bbrst9957L2eddVbWkSQVGcuSJEnqsnbZZRfuvvtuevXqxaRJk5g0aVLWkSQVEcuSJEnq0oYNG0ZVVRUAZ5xxBvfcc0/GiSQVC8uSJEnq8o477jguvvhi1q1bxzHHHMPChQuzjiSpCFiWikBEXBsRkyPimoi4PyK+nHUmSZK6mksvvZRjjjmGlStXUlFRwfLly7OOJCljlqUNiIjtIuKeiEjt8HGrU0onp5TOAf4/YHI7fKYkSWqgW7du/PSnP2Xo0KH89a9/ZfTo0dTV1WUdS1KGLEt5RMRo4FHgkx/xuj4RMS0iltQ/7oqITzT381JK5zU43AN4qrnnkCRJrbfpppvyq1/9in79+jF37ly+8Y1vkFJ7/LuppGJkWcrvAuAQ4JENvSAiegH3Ab2AgUA5UAv8LiI2a+4HRsTeETED+Hr9Q5IkZeDjH/84NTU1bLbZZkybNo3LL78860hSp7VoEQwalHsuRpal/PZPKT3/Ea85ARgMnJ9SWpNSWgucD+wKnLr+RRHxaET8bQOPf89CpZSeSikdBVwM/CEiygr/Y0mSpKYYPHgwv/jFL4gILrnkEqZPn551JKnTqa2FESNg8WIYOTJ3XGwsS3mklJpyC+8vA6+klF5s8L5/AIvrv7d+7HMppU9s4PG3iOjecCYqpXQfsDkwpHA/kSRJaq5Ro0ZxzTXXAHDiiSfy+OOPZ5xI6lwqK2H5ckgJli2DceOyTvRhlqWWGwy8lGf8JeDTzThPPxps6BAR25MrS/nOLUmS2tGZZ57JySefTF1dHUcccQSvvPJK1pGkTqGqCmbPhvV7qNTVQXV1bryYWJZa7mPAyjzj7wClEbFpE8/zFtA9IqZGxPeBKuDElFLe38YRcXJEzIuIea+//nqLgkuSpKaJCCZNmsTw4cNZtmwZFRUVrFyZ7z//kprjwgs/vOzu3Xdz48XEslR40ZwXp5TeSSmNSSmNTSmdnVI6LKV0+0ZePzmlNCSlNGTbbbdtfVpJkrRRPXv25M4772S33XZjwYIFfO1rX2Pt2rVZx5I6tAkToKzRFfqlpTBxYjZ5NsSy1HJvkFsu19jmwLsppVXtnEeSJLWRrbbaipqaGrbeemtqamo499xzs44kdWiVlblNHUpKcsclJVBRAWPHZpurMctSyy0Ads4zvgvwTPtGkSRJba1///7MmDGDnj17cu2113LTTTdlHUnq0KqqoE8fiIC+fWHKlKwTfZhlqeVmADtFxM7rByKiL7mbyv4yq1CSJKntHHTQQf8uSaeddhr3339/xomkjqusDObMgfLy3GYPjZflFQPLUsv9lNwM0pUR0SMiugETye1i9+Msg0mSpLYzduxYzj//fNauXcvRRx/Ns88+m3UkqcMaOBAWLsw9FyPLUh4RcXVEPA0cXn/8dP2j1/rXpJTeBw4B1pK7t9KfgS2Ag1NK/2r/1JIkqb1cccUVjB49mhUrVjBq1CjefPPNrCNJagORUso6g1poyJAhad68eVnHkCSpS6qtreXAAw/kT3/6EwceeCC/+c1v2GSTTbKOJall8u5o7cySJElSC5SVlTFr1iy23357HnroIU455RT8R2ipc7EsdUARURERk1esWJF1FEmSurQddtiB6upqSktLueWWW7jqqquyjiSpgCxLHVBKqTqldHLv3r2zjiJJUpe3zz77cNtttwFwwQUXMGPGjIwTSSoUy5IkSVIrjR49mokTJwJw/PHH8+STT2acSFIhWJYkSZIK4LzzzmPs2LGsWrWKww8/nFdffTXrSJJaybIkSZJUABHBjTfeyIEHHsjf//53KioqqK2tzTqWpFawLEmSJBVIr169mDFjBp/85Cd56qmnOP7441m3bl3WsSS1kGVJkiSpgLbZZhtqamrYcsstufvuu7nooouyjiSphSxLkiRJBbb77rtz11130b17d6688kqmTp2adSRJLWBZkiRJagPDhw/nRz/6EQCnnHIKDz74YMaJJDWXZUmSJKmNnHzyyZx99tmsXr2ao446ir/85S9ZR5LUDJalDigiKiJi8ooVK7KOIkmSPsJVV13FqFGjeOuttxg5ciT//Oc/s44kqYksSx1QSqk6pXRy7969s44iSZI+Qvfu3fn5z3/O4MGDee655zj66KNZvXp11rEkNYFlSZIkqY1tvvnmVFdX07dvXx544AHGjx9PSinrWJI+gmVJkiSpHey4447MmjWLkpISJk+ezHXXXZd1JEkfwbIkSZLUToYOHcott9wCwDnnnEN1dXXGiSRtjGVJkiSpHR1zzDFcdtllpJQ49thjmT9/ftaRJG2AZUmSJKmdffvb3+a4446jtraWiooKXnvttawjScrDsiRJktTOIoKbb76ZYcOGsXTpUo444ghWrVqVdSxJjViWJEmSMlBSUsLMmTPZeeedeeKJJzjhhBNYt25d1rEkNWBZkiRJykifPn2oqalhiy224M477+TSSy/NOpKkBixLkiRJGRo4cCB33HEH3bp147LLLmPatGlZR5JUz7LUAUVERURMXrFiRdZRJElSARx66KH84Ac/AKCyspJHHnkk40SSwLLUIaWUqlNKJ/fu3TvrKJIkqUBOO+00xo8fz/vvv8+RRx7Jiy++mHUkqcuzLEmSJBWJa6+9lkMPPZQ33niDiooKXEUiZcuyJEmSVCR69OjB9OnTKS8vZ/HixYwZM4Y1a9ZkHUvqsixLkiRJRaR3797U1NSw7bbbcu+993LWWWdlHUnqsixLkiRJRWaXXXbh7rvvplevXkyaNIlJkyZlHUnqkixLkiRJRWjYsGFUVVUBcMYZZ3DPPfdknEjqeixLkiRJReq4447j4osvZt26dRxzzDEsXLgw60hSl2JZkiRJKmKXXnopxxxzDCtXrqSiooLly5dnHUnqMixLkiRJRaxbt2789Kc/ZejQofz1r39l9OjR1NXVNfs8Tz/9NCeddBL77rsve+65J+Xl5Zx++um8/vrrbZBa6hwsS5IkSUVu00035Ve/+hX9+vVj7ty5fOMb3yCl1KxzfPWrX+Wtt97ioYceYv78+dx333385je/Yf/992fVqlVtlFzq2CxLkiRJHcDHP/5xampq2GyzzZg2bRqXX355s89x5ZVXUlZWBsAOO+zAueeey/PPP8+cOXMKHVfqFCxLkiRJHcTgwYP5xS9+QURwySWXMH369Ca/d8GCBfTv3/8DY9tvvz0A//znPwuaU+osLEsdUERURMTkFStWZB1FkiS1s1GjRnHNNdcAcOKJJ/L444836X29evX60Nhzzz1HRHDggQcWNKPUWURz17uqeAwZMiTNmzcv6xiSJKmdpZT45je/yeTJk+nbty9//OMf2XHHHZt1jrVr17LPPvswdOhQfvKTn7RRUqnDiHyDzixJkiR1MBHBpEmTGD58OMuWLaOiooKVK1c26xyXXXYZPXr04Nprr22jlFLHZ1mSJEnqgHr27Mmdd97JbrvtxoIFC/ja177G2rVrm/TeqVOncscdd/DrX/+azTbbrI2TSh2XZUmSJKmD2mqrraipqWHrrbempqaGc8899yPfc+utt3LNNdfwwAMP0KdPn3ZIKXVcliVJkqQOrH///syYMYOePXty7bXXctNNN23wtbfddhtXXnkl999/Px//+McBqKmpYfLkye0VV+pQemQdQJIkSa1z0EEHcdNNN1FZWclpp53GJz/5Sb74xS9+4DXTpk3jpJNO4rLLLuP+++//9/jDDz/Mdttt196RpQ7BsiRJktQJjB07liVLlnDjjTcS8eGNvb71rW9RV1eXd6ned7/73faIKHU4bh3egbl1uCRJamjdunUsXbqUnXbaKesoUkfj1uGSJEmdWbdu3SxKUgFZliRJkiQpD8uSJEmSJOVhWZIkSZKkPCxLkiRJkpSHZUmSJEmS8rAsSZIkSVIelqUOKCIqImLyihUrso4iSZIkdVqWpQ4opVSdUjq5d+/eWUeRJEmSOi3LkiRJkiTlYVmSJEmSpDwsS5IkSZKUh2VJkiRJkvKwLEmSJElSHpYlSZIkScrDsiRJkiRJeViWJEmSJCkPy5IkSZIk5WFZkiRJkqQ8LEuSJEmSlEeklLLOoBaKiNeBt4EVrTxV71ac42PAG638fBVOa/637Eg6ws+Zdcb2+vy2+pxCnrcQ5/L3ZOeQ9f9ftpeO8HNmnbE9P9/fkx+tGH5PvpFSOqzxoGWpg4uIySmlk7M6R0TMSykNac3nq3AK8eehI+gIP2fWGdvr89vqcwp5Xn9Par2s//+yvXSEnzPrjO35+f6ebNJ7i/b3pMvwOr7qIjmHikNX+d+yI/ycWWdsr89vq88p5Hn9Pan1usr/jh3h58w6Y3t+vr8nOzBnltQqxfwvAZJUDPw9KUkbV8y/J51ZUmtNzjqAJBU5f09K0sYV7e9JZ5YkSZIkKQ9nliRJkiQpD8uSJEmSJOXRI+sAkiR1VRGxNXA18C6wDtgZ+O+U0vNZ5pKkYhER1wJlwEpgT+DHKaVfttfnO7OkNhURPSPivyOiNiIGZZ1HkgolIraLiHsiojUX/+4I1KWUvpVSOgO4D7i5MAklKTsF+h0JsDqldHJK6Rzg/6OdN4OwLCmvAv4BPxl4BCgtQCxJKgoRMRp4FPjkR7yuT0RMi4gl9Y+7IuIT67+fUno6pXRag7e8BGzXNqklqX0U6nckQErpvAaHewBPFT7xhlmW9CEF/gP+w5TSo22ZV5IycAFwCLl/DMorInqRmynqBQwEyoFa4HcRsdkG3lYB3FjYqJLU7gr6OzIi9o6IGcDX6x/txrKkfNrqLwGS1Fns34Trik4ABgPnp5TWpJTWAucDuwKnNn5xRIwEegPXFTirJLW3gv6OTCk9lVI6CrgY+ENElLVF6HwsS8qn4H8JkKTOJKW0pgkv+zLwSkrpxQbv+wewuP57/1ZflI4A/iultK6QWSWpvRXqd2REdG/4j/AppfuAzYEhhU28YZYlfUih/xIgSV3UYHLXIDX2EvDp9QcR8RVys/mnpJTWRMT17ZRPkrLUlN+R/WiwoUNEbE+uLOV7X5tw63C11GDguTzjLwHD2zmLJBWjjwFP5hl/ByiNiE2BTwG/AN4AvhoRkFuKd0Z7hZSkjDTld+RbQPeImAr8k9xlHyemlF5pr5CWJbXUR/4BTymtiogDgWPqv3dRRPwqpTS93VJKUvGJ9V+klBbgf4slqaGGvyPfAcZkmMVf0Cq4aHiQUnoIeAgYn00cScrMG+SWizS2OfBuSmlVO+eRpGLSIX5Hes2SWqpD/AGXpAwtAHbOM74L8Ez7RpGkotMhfkdaltRSHeIPuCRlaAawU0TsvH4gIvqSu6niL7MKJUlFokP8jrQsqaU6xB9wScrQT8n949GVEdEjIroBE8lthPPjLINJUhH4KR3gd6RlSS31UzrAH3BJagsRcXVEPA0cXn/8dP2j1/rXpJTeJ7cl+Fpyt1X4M7AFcHBK6V/tn1qS2kdn+h0ZKaWsM6jIRMTV5P7w7ghsBcyv/9bQ+j/Y61/XF7iW3I3BErAQODOltLR9E0uSJEmFZ1mSJEmSpDxchidJkiRJeViWJEmSJCkPy5IkSZIk5WFZkiRJkqQ8LEuSJEmSlIdlSZIkSZLysCxJkiRJUh6WJUmSJEnKw7IkSZIkSXlYliRJkiQpD8uSJEmSJOVhWZIkqYGIWBQRP4mI/4qIpyJiVUQsjogvZJ1NktS+IqWUdQZJkopCRJQA/wKWAn8GbgZWA1cDZSmlfhnGkyS1sx5ZB5AkqYh8GugOPJlSOnr9YET0A34YEZumlFZllk6S1K5chidJ0n/sXf98SaPxjwHvFLooRcTbEbFzIc8pSSocy5IkSf+xF/C3lNLiRuN7AwvaP44kKUuWJUmS/mNv4MkNjP9p/UFEpIi4ICIej4iXImJsg+8dEhHzImJBRDwSEfs0+N4REfHniHg0IiY2/ICI+ExEPFD/3qci4iv145tGxPT6TSbmR8RvCv5TS5Ly8polSZKAiOhG7pqlXzca3wrYCXiq0VtWp5T2i4hy4I8RcSuwFXA7MDyl9HREjARmRsSngC2BqcABKaVFEXES0Lv+M7YEfgKMSCn9PSI+BvwpIuYCQ4EtU0rl9a/dug1+fElSHs4sSZKUMwAo48MzS+uvY/pTo/FpAPVL9hLwceCzwMKU0tP135vd4NyfBRaklBbVj00F1tR/PQzYBZgTEU8D99eP7wbMB8oj4saI+GqD90iS2pgzS5Ik5exV/9y4FO0NvAc0vo6prsHXa/nPf1M3dk+Oxt+LBs+LUkrD8r0pIvYADga+BFwVEXumlP65kc+RJBWAM0uSJAEppdtTSpFSeq3R+DUppZKUUlNmdB4DBkXEYICI+L/kCtJzwKPA4IgYWP/aseS2KQeYC+wSEYesP1FE7BURvSLiE7kYaRZwdv17vN+TJLUDZ5YkSSqQlNLrEXEsMDUiepC7we1RKaX3gdcjohL4ZUT8E5gDrKh/3z/rr2/6XkR8D+gJvAIcSe46qokREeT+kfOnKSV35pOkdhApbWy1gCRJkiR1TS7DkyRJkqQ8LEuSJEmSlIdlSZIkSZLysCxJkiRJUh6WJUmSJEnKw7IkSZIkSXlYliRJkiQpD8uSJEmSJOXx/wNTDgLFOUt6jgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#HIDDEN\n",
    "fig, ax = plt.subplots()\n",
    "disc_n = np.logspace(1,3,20)\n",
    "analytical_top = reference_solution[-1]\n",
    "for nn in disc_n:\n",
    "    numerical, z = run_sim(int(nn))\n",
    "    numerical_top = numerical[-1]\n",
    "    relerr = np.abs((numerical_top - analytical_top)/analytical_top)\n",
    "    ax.plot(int(nn),relerr,ls='',marker='d',color='blue')\n",
    "ax.plot([1e2,3e2],[5e-3,5e-3/3],ls='-',color='black')\n",
    "ax.text(3.1e2,5e-3/3,'1')\n",
    "ax.plot([1e2,3e2],[5e-3,5e-3/9],ls='-',color='black')\n",
    "ax.text(3.1e2,5e-3/9,'2')\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel('$n_\\\\mathrm{nodes}$')\n",
    "ax.set_ylabel('$\\\\left| \\\\frac{u_z^\\\\mathrm{numerical} - u_z^\\\\mathrm{analytical}}{u_z^\\\\mathrm{analytical}} \\\\right|_{z=H}$',size=24)\n",
    "fig.tight_layout();"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}