06aa_soil_column_FDM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "<div style=\"display: flex; align-items: center;\">\n",
    "    <img src=\"https://github.com/nagelt/Teaching_Scripts/raw/9d9e29ecca4b04eaf7397938eacbf116d37ddc93/Images/TUBAF_Logo_blau.png\" width=\"500\" height=\"auto\" height=\"auto\" style=\"margin-right: 100px;\" />\n",
    "    <div>\n",
    "        <p><strong>Prof. Dr. Thomas Nagel</strong></p>\n",
    "        <p>Chair of Soil Mechanics and Foundation Engineering<br>Geotechnical Institute<br>Technische Universität Bergakademie Freiberg.</p>\n",
    "        <p><a href=\"https://tu-freiberg.de/en/soilmechanics\">https://tu-freiberg.de/en/soilmechanics</a></p>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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": 19,
   "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": 20,
   "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": 20,
     "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": 21,
   "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": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Central_Diag = np.diag([-2.]*number_of_points,0)\n",
    "Central_Diag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 23,
   "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": 23,
     "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": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.999999999999996\n"
     ]
    }
   ],
   "source": [
    "u = np.linalg.solve(A,RHS)\n",
    "print(np.linalg.det(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAIaCAYAAAAqfyubAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABoWklEQVR4nO3dd3hUZd7/8c+k90ISEkhCaAlVivSOrthF14Krguiqq6vurq7t8Xncn+XZR3dXVl3Fte6CAnZde5cqEBSUriShhISakDLp9f79kWRgSJmElJnMvF/XlWuXc58z883xJPnk5D7f22KMMQIAAAA8jJezCwAAAACcgSAMAAAAj0QQBgAAgEciCAMAAMAjEYQBAADgkQjCAAAA8EgEYQAAAHgkgrADxhhZrVbRbhkAAMC9EIQdKCoqUnh4uIqKipxdCgAAADoQQRgAAAAeiSAMAAAAj0QQBgAAgEciCAMAAMAjEYQBAADgkQjCAAAA8EgEYQAAAHgkgjAAAAA8EkEYAAAAHokgDAAAAI9EEAYAAIBHIggDAADAIxGEAQAA4JG6RRAuLy/X//zP/8jPz0/XXXedw/0/++wznXnmmQoPD1dYWJimTp2qd955p/MLBQAAQLfh8kF4xYoVGjFihJ577jlVVVU53P+pp57S+eefr+TkZKWlpWnv3r0aP368rrjiCj300EOdXzAAAAC6BYsxxji7iOa8/vrr+u1vf6tHH31UgYGB+vWvf6358+dr8eLFTe6/fft2jRo1Sqeffro2bNggi8ViGzvjjDO0cuVKrV+/XhMnTmx1DVarVeHh4SosLFRYWFh7PyUAAAC4CJe+I9yvXz/t3LlTt956q12obc4TTzyhmpoa3XTTTY32v/nmmyVJjz/+eKfUCgAAgO7Fx9kFtKQtd24l6aOPPpIkTZkypdHY1KlTJUmffvqpqqur5ePj0p86AAAAOplL3xFuiwMHDig3N1dS3Z3kk8XHx8vPz0/l5eVKS0vr6vIAAAA8QllljTZl5unV9fvkwjNwJbn4HeG2OHz4sCQpICBAgYGBjcYtFosiIiJ09OhRHT58WEOHDm3ydSoqKlRRUWH7t9Vq7ZyCAQAAurmqmlrtOlykLdkF2ppVqC3ZBUo/Wqya2roAPCMlRklRwU6usnluE4RLS0slSX5+fs3u4+/vb7dvUx577DE9/PDDHVscAABAN1dba7Qnt1hbsgq1NbtAW7ILtfOQVZXVtc0esyW7kCDcFYKCgiRJlZWVze7TcKe3Yd+m3H///frjH/9o+7fValViYmIHVQkAAOD6jDHKzi/T1uyG0Fug7QesKq6obvE4by+LknuGaGRChEYkhmtMUmQXVXxq3CYIx8XFSapbfKOsrKzR9AhjjAoKCiRJvXr1avZ1/P39bXeOAQAAPEFOUYXtLu/W7AJtyy7UsZLmby426BcdrBEJ4RqREKGRCeEa1jtcgX7eXVBxx3CbIBwfH6/o6Gjl5uZq7969jeYAHzhwQJWVlQoMDFRycrKTqgQAAHAua3mVtmUX2ub1bs0u0MHCcofH9QoPOCH0Rui0+HCFB/l2QcWdx22CsCRddNFFWrRokdatW9coCK9du1aSdN5559E6DQAAeITyqhrtOFhom9e7NbtQe3JLHB4XGeRru8s7on6aQ8/QgC6ouGu5VSK866679Oqrr+rll1/WDTfcYLeoxosvvihJuueee5xVHgAAQKdp6OBwfF5vodKOFNk6ODQn2M9bw+PDNTIxQiMSwjUyIUIJkYGtWsysu3OrIDxs2DAtWLBAd955p2699VY9/PDD8vHx0aOPPqrly5frwQcfbPMiHQAAAK6mroNDie0u75bsAu08aFVFCx0cJMnP20tDeofZ7vSOTAhX/5gQeXu5f+htisW4eKfjln4bWbRoka677rpG2z/99FM9/vjj2rRpk4wxGjFihO68805dfvnlbX5/q9Wq8PBwFRYWKiwsrM3HAwAAtIcxRgcKymyBd2tWobYfKFSRgw4OXhYpJTbUbl7voLhQ+fm4zXpq7ebyQdjZCMIAAKAr5RbXd3A4YV5vazo49I0KqpvPm1A3zWFY7zAF+bnVH/87HGcHAADASazlVdqeXWhrW7Y1u1AHCsocHhcXFmALvCMSwjUiPqLbd3BwBoIwAABAF6jr4GC1m9e7J8dxB4eIkzo4jEwIV88w9+vg4AwEYQAAgA5WVVOrtCMndHDIquvgUO2gg0NQQweHE+b1JvbwjA4OzkAQBgAAaIfaWqO9x0rs5vXuaG0Hh16hdvN6B3hwBwdnIAgDAAC0kjFGBwvLtTXLfjni1nRwSO5Z38EhsW56w+C4MDo4OBlBGAAAoBnHiiuOty2rD765xY47OCTVd3BomOIwrHeYgv2JXa6G/yIAAACSisqrtO1Aod283tZ0cIgN87dfjjghXBFBfl1QMdqLIAwAADxOeVWNdh6yamvWCR0cckvkaHWFiCBfnRZftwxxw7zeWDo4dFsEYQAA4Naqa2q164QODluzC7XrcCs7OPQOt5vX26dHEB0c3AhBGAAAuI1T7eDg623RkF5hdssRD+xJBwd3RxAGAADdUpMdHA4Uqqi8dR0cTks43q93cK9Q+ft4d1HlcBUEYQAA0C2cageHPj2C6ubz1s/rHR4fTgcHSCIIAwAAF3SqHRx6hp7QwSExQiPiwxUZTAcHNI0gDAAAnOpUOziEB/ra3emlgwPaiiAMAAC6THVNrdKOFNfd5a2/29uaDg6Bvt46LZ4ODuhYBGEAANApamuN9h0rsZvXu+Ngocqr2t7BYUBMsHy8WY4YHYsgDAAA2s0Yo0OF5XZ3erdmO+7gYLFIyT1D7FZmo4MDugpBGAAAtFleSWXdXd76Xr1bsguVW1zh8Dg6OMCVcOUBAIAWFVdUa9sJd3m3ZBcoO58ODuj+CMIAAMCmvKpGPx2y2s3r3Z1T7LCDQ1iAj0YmRtjN640Lp4MDXBtBGAAAD1VdU6v0o/YdHH4+1LoODsPjwzSioW1ZQoSSoujggO6HIAwAgAcwxmjfsVLb4hRbswu0vZUdHAbHhR2f15sYroExIXRwgFsgCAMA4GaMMTpsLbcF3obV2ayt6OAwMKa+g0NifQeHuFAF+NLBAe6JIAwAQDeXV1JpF3i3ZBcqp8hxB4fEHoF2bcuGx4crhA4O8CBc7QAAdCPFFdXafqDQbl5vVp7jDg4xof62wNvwQFsPOjjAwxGEAQBwUe3p4HBi4B2ZGK64sAAeZgNOQhAGAMAFNNXBYdfhIlXVtJx6A3y9NLx3uN283r50cABahSAMAEAXa6qDw46DVpVV1bR4nI+XRYN7hdrN603uSQcH4FQRhAEA6ETt6eAwICbEbjniIb3C6OAAdCCCMAAAHSi/pNI2n7ctHRwSIgNtgbeug0OYQgN8u6BiwHMRhAEAOEUlFdXadgodHKJDTujgkBiuEfHhigrx74KKAZyIIAwAQCtUVNfop0NFdvN6M1rRwSE0wOd494b6/+0VTgcHwBUQhAEAOElNrVH60SJtzTretuznw9ZWdXAY1jvcbl5v36hgeXkRegFXRBAGAHg0Y4wyj5XazevdfqB1HRwGxdl3cEiJpYMD0J0QhAEAHuVwYXl96G0IvoUqLKtq8RiLReofHXz8YbbECA2lgwPQ7RGEAQBuK7+kUlsPFGpr1vGH2Y62ooNDfESgbXGKEQnhOi0+nA4OgBsiCAMA3Mb+Y6ValZ6jDXuOaWt2ofbnlTo8JjrEzxZ4G+740sEB8AwEYQBAt1VcUa31u49pdVqOVqfnKPNYy8E31N9Hp53YwSExQr3p4AB4LIIwAKDbqK012n6wUGvSc7UqLUc/ZOarurbpTg7+Pl4aHk8HBwDNIwgDAFzaEWu5VqflaE16rr7NyFVeSWWT+/l4WTQmKVLTU2I0LTlaQ3qFyZcODgBaQBAGALiU8qoafb8vzxZ+fz5c1Oy+faOC6oNvjCYNiFKIPz/WALQe3zEAAE5ljFH60eL6eb652rDnmCqqa5vcN8TfR5MHRGl6SoymJ8eoT1RQF1cLwJ0QhAEAXS6/pFLfZuRqTXrdXd9DheVN7mexSCPiw+uCb0qMRiVGMN0BQIchCAMAOl1VTa02ZxXY7vpuzS6QaWa14riwAE1Ljtb0lBhNHRityGC/ri0WgMcgCAMAOkVWXqlWpeVodVqO1u8+pqKK6ib38/fx0oT+UZpeH36Te4bQzgxAlyAIAwA6RHFFtVJ3H9Pq9Lrwu6+Fnr6DYkM1PSVa05JjNL5fD5YqBuAUBGEAwCmprTXacdBqC74/7M9XVU3T8x0ig3w1NTnGdtc3Niygi6sFgMYIwgCAVjtqLdfq9LqH3L5Nz9WxFnr6np4UaQu+w3uHs5AFAJdDEAYANKu8qkYb9+Xb7vq21NM3KSpI05PrFrOYNCBKoQG+XVgpALQdQRgAYGOMUcbRYq1Oz9XqtBxt2HtM5VXN9/SdZOvpG62kqOAurhYA2ocgDAAerqC0vqdvWq5Wp+e02NP3tPhwTU+u6+k7ug89fQF0bwRhAPAw1U309K1tpqdvbJi/ptUH36kDo9WDnr4A3AhBGAA8QFZeqW2e77qMlnv6ju/Xw3bXNyWWnr4A3BdBGADcUElFtVL3HLPd9d2bW9LsvimxIbbgS09fAJ6EIAwAbqC21mjnIatWpeVoTXqONmU239M3IshXUwdG1z/kFqO4cHr6AvBMBGEA6KaOFpXbHnBz2NO3T6SmNfT0jQ+XNz19AYAgDADdRXlVjTZl5mt1Wo5WOejp26dHkG0J48n09AWAJhGEAcBFGWO0O6dYq+vv+qbuab6nb7CftyYNiNaM+vDbN5qevgDgCEEYAFxIYWlVXU/f+g4PB1vo6Tu8d7imp0RrenKMRveJlJ8PPX0BoC0IwgDgRNU1tdqSXaBVaXXhd0tW8z19e4Y29PSN1tSB0YoK8e/aYgHAzRCEAaCLZeWVak39EsZrd+eqqLzpnr5+Pl6aUN/Td1pKtAbFhtLTFwA6EEEYADpZQ0/fhvC7p4Wevsk9Q+ramqXEaHzfHgr0o6cvAHQWgjAAdLCGnr6r03O0Ji1XGzPzWuzpO2VgtGbU3/XtFR7YxdUCgOciCANAB8gpqrA94PZtRq5yi5vu6evtZdHpfSLq5/rG6DR6+gKA0xCEAeAUVFTXaNO+fK1Kz9HqtFz9dMja7L6JPQLr5vkmx2jywCiF0dMXAFwCQRgAWqGup2+J7a5v6p48lVXVNLlvkJ+3Jg+I0vSUuvDbNyqIh9wAwAURhAGgGYWlVVq7u6Gnb64OFJQ1u+/w+DBNr5/ucDo9fQGgWyAIA0C9up6+hVqdlqPVDnr6xoT6a1pytGakxNDTFwC6KbcMwmvWrNGTTz6pzZs369ChQ4qOjtbgwYP1+9//XhdddJGzywPgQrLzT+jpm5Eraws9fcf37aFpydGanhKjwXH09AWA7s7tgvDChQv1u9/9TsOHD9fSpUs1atQoHThwQPfff79mz56t3/3ud3r66aedXSYAJymtrOvpuzotV6vTc7Qnp/mevgN7htRPd4jWhH5R9PQFADdjMcY084e/7qeqqkrR0dGyWq367rvvNG7cONtYWVmZEhISlJeXp59//lmDBg1q1WtarVaFh4ersLBQYWFhnVU6gE5SW2v002GrVtcvYbxxX74qa2qb3Dc80FdTB0Zrekq0piXHqHcEPX0BwJ251R3h/Px8Wa11LYyGDRtmNxYYGKgBAwYoLy9PmzdvbnUQBtD95BRV6NuMnPrwm6vc4oom9/P2smh0YkNP32iNSIigpy8AeBC3CsI9e/ZUYmKisrKytGPHDrs7wuXl5dq9e7ckqVevXs4qEUAnqKiu0abM/LrpDmk52tlCT9+EyMC6JYyTozVpQLTCA+npCwCeyq2CsCQtXbpUc+bM0a9//Wu9+OKLGjVqlLKzs3X//fcrLy9P48aN05QpU5xdJoB2MMZoT26J1qTlaHV6rlL3HFNpZfM9fSf1b+jpG61+0cE85AYAkORmc4Qb7N+/X3fddZfeeecd27aYmBhdccUVevTRRxUeHt7ssRUVFaqoOP5nVKvVqsTEROYIA05WWFaldRl1D7g56uk7rHdY/V3fGI1JoqcvAKBpbndH+KOPPtJ1112nvn376ttvv9XIkSOVnZ2tl156SYWFhSopKWkxCD/22GN6+OGHu7BiAE2pqTXakl2g1Wk5WpOeq81ZBapppqlvQ0/f6ckxmpocrWh6+gIAWsGt7ghnZmZq6NCh8vLyUkZGhmJjY21jxhiNGzdO+/bt0w8//KA+ffo0+RrcEQac50BBWf10hxx9m95CT19vL43rF6npyXVLGA/pRU9fAEDbudUd4TfeeEOlpaWaPXu2XQiWJIvFoquuukp33323HnroIf373/9u8jX8/f3l78/dJKArlFZWa8OevPrpDjna3UJP3wExwXXTHVJiNKFfDwX5udW3LwCAE7jVT5K9e/dKar4rRO/evSVJmzZt6rKaABxnjNFPh4q0Oj1Ha9Jz9P3e5nv6hgX4aGr9dIdpKTGKp6cvAKCDuVUQjo6OliQdPHiwyfGG7b6+tEsCukpucYW+rV/CeE1GrnKKmu7p62WRRveJtC1hPJKevgCATuZWQXj27Nn6v//7P61YsUJHjx5Vz549bWPGGL3xxhuSpFmzZjmrRMDtVVbXamNmntbUh98dB5vv6Rsfcbyn7+SB9PQFAHQtt3pYTpLuuusuPfHEExo7dqyefvppjRgxQtnZ2frf//1fLVu2TMOHD9eaNWsUERHRqtdjiWWgZcYY7c0tsXV3WN9CT99AX29NGhCl6cnRmpYSo/709AUAOJHbBWFJev/99/XCCy9o48aNys/PV2BgoAYNGqRLL71Uf/jDHxQcHNzq1yIIA41Zy+t6+q5Ky9Wa9Bxl57fc07dhCeMxSZHy9/HuwkoBAGieWwbhjkQQBurU1Bot//molqRmam1GbrM9faND/Ovv+EZr6sAYxYTShQUA4Jrcao4wgI6XU1ShtzZm6bUN+5tczc3P20tj+0baljAeEhcmLx5yAwB0AwRhAI0YY7QxM19L1mfqs+2HVFVjf/c3PiJQs4bGakZKjCb0p6cvAKB74qcXAJviimq9/+MBLU3N1M+Hi+zGLBbpjEE9NW9ikqanxNDaDADQ7RGEASjtSJGWpmbqvR8OqLjCflnjHsF+mjM2UddM6KPEHkFOqhAAgI5HEAY8VGV1rb7ceVhL1mdqw968RuOn94nQvElJOm94LwX40ukBAOB+CMKAhzlUWKbXN+zX699nNVrlLdDXW5eM7q1rJiRpeHy4kyoEAKBrEIQBD1Bba7Ru9zG9un6fvv7piE7ufNY/JljzJibp0tMTWN0NAOAxCMKAGyssrdLbm7K0bMN+7c0tsRvz9rLo7KGxmjcxSZMGRLHCGwDA4xCEATe0LbtQS1L36cMtB1VeVWs31jPUX1eN76OrxvdRXHiAkyoEAMD5CMKAmyivqtHHWw9pSWqmtmQVNBqf1D9K8yYladbQWPl6e3V9gQAAuBiCMNDNZR4r0bIN+/XWxiwVlFbZjYX6++iyMQmaO7GPBvYMdVKFAAC4JoIw0A3V1Bqt+PmolqRmalVaTqPxIb3CdO2kJM0e2VvB/nyZAwDQFH5CAt1IbnGF3vw+S69t2K8DBWV2Y37eXjr/tDjNm5Sk0/tE8vAbAAAOEIQBF2eM0abMfC1JzdSn2w6pqsa+91l8RKDmTkzSnLEJigrxd1KVAAB0PwRhwEWVVFTr/c0HtGR9pn4+XGQ3ZrFIM1NiNG9Skmak9JS3F3d/AQBoK4Iw4GLSjxRpaWqm3v3hgIorqu3GIoN8NWdcoq4Zn6Q+UUFOqhAAAPdAEAZcQFVNrb7ccURLUvcpdU9eo/HRfSI0b2KSzj+tlwJ8vZ1QIQAA7ocgDDjRocIyvf5dlt74br+OFlXYjQX4eumSUfGaOzFJw+PDnVQhAADuiyAMdDFjjNbtPqYl6zP11U9HVFNr//Bb/+hgzZ2YpMvGJCg80NdJVQIA4P4IwkAXKSyr0rubsrV0Q6b25JTYjXl7WTRrSKzmTUrS5AFRtD4DAKALEISBTrb9QKGWpmbq/c0HVF5VazcWE+qvq8b30VXjE9UrPNBJFQIA4JkIwkAnKK+q0afbDmlJaqZ+3F/QaHxi/x6aN7Gvzh4WK19vr64vEAAAEISBjrT/WKmWfZept77PUn5pld1YiL+PLju97uG35NhQJ1UIAAAaEISBdqqpNVqVdlRL1mdqZVqOjP2zbxocF6p5k5J0yah4BfvzJQcAgKvgpzJwio4VV+itjdlatiFT2flldmN+3l4677Q4zZuYpDFJkTz8BgCACyIIA21gjNEP+wu0NDVTn2w9pMoa+4ff4iMCdc3EPpozNlHRIf5OqhIAALQGQRhohdLKan2w+aCWrM/UzkNWuzGLRZqREqN5E5M0c1BPeXtx9xcAgO6AIAy0IONosZamZurdTdkqqqi2G4sI8tWVYxN19YQ+SooKdlKFAADgVBGEgZNU1dTq651HtCQ1U+t2H2s0PioxQvMmJumCEb0U4OvthAoBAEBHIAgD9Q4Xluv17/brje/364i1wm4swNdLF4+sa312WkK4kyoEAAAdiSAMj2aM0frdx7QkNVNf7jyimlr73mf9ooM1d2KSLj89QeFBvk6qEgAAdAaCMDyStbxK723K1pLUTO3OKbEb87JIs4bGat7Evpo8IEpePPwGAIBbIgjDo+w4WKilqZl6/8eDKquqsRuLDvHX1eMT9avxfdQ7ItBJFQIAgK5CEIbbK6+q0WfbD2nJ+kz9sL+g0fiEfj00b1KSzh4aJz8fr64vEAAAOAVBGG4rK69Uyzbs11sbs5RXUmk3FuLvo0tPr3v4LSU21EkVAgAAZyIIw63U1BqtTsvRktRMrdh1VMb+2TcNjgvV3IlJumR0vEL8ufwBAPBkJAG4hbySSr21MUvLNmQqK6/MbszX26LzhvfSvElJGpsUKYuFh98AAABBGN2YMUY/ZhVo6fpMfbztkCqra+3G4yMCdfWEPpozNlExof5OqhIAALgqgjC6ndLKan24+aCWpGZqx0Fro/EZKTGaNzFJZwzuKW9anwEAgGYQhNFt7M4p1tLUTL2zKVtF5dV2Y+GBvpozNkHXTEhS3+hgJ1UIAAC6E4IwXFp1Ta2+/umIlqRmam3GsUbjIxMjNHdCH100srcCfL2dUCEAAOiuCMJwSUet5Xr9uyy9/t1+HbaW2435+3jp4lG9NXdikkYkRDinQAAA0O0RhOEyjDFK3ZOnpamZ+mLHYVXX2vc+6xsVpLkTk3T5mARFBPk5qUoAAOAuCMJwOmt5lf7zwwEtSc1UxtFiuzEvi3TWkFjNm5SkKQOi5cXDbwAAoIMQhOE0Ow9atXRDpt7/8YBKK2vsxqJD/HXV+ERdNb6PekcEOqlCAADgzgjC6FJVNbX6dNshLVmfqY2Z+Y3Gx/froXkTk3TOsDj5+Xg5oUIAAOApCMLoMjlFFbpl6SZtOikAB/t569LTEzR3YpIGxYU6qToAAOBpCMLoEjsOFuqmVzbqYOHxDhCDYkM1d1KSfjk6XiH+XIoAAKBrkT7Q6T7ffkh3vrlFZVV184DjwgL0xJyRmjQgShYLD78BAADnIAij0xhj9MzyDD3xVZpt26jECL04b4x6hgU4sTIAAACCMDpJWWWN7nlniz7eesi27Zej4/XYpaexAhwAAHAJBGF0uMOF5brp1Y3adqBQkmSxSPeeM1i3zOjPVAgAAOAyCMLoUJuzCvSbVzfqaFGFpLqOEE/9arRmDY11cmUAAAD2CMLoMB9sPqB739mqiupaSVJCZKBenj9Wg+PCnFwZAABAYwRhtFttrdETX6Vp4YoM27bxfXvoubmnKyrE34mVAQAANI8gjHYpqajWnW9u1pc7j9i2/Wpcoh65eDgrwwEAAJdGEMYpy84v1Y2vbNTPh4skSV4W6YELhur6KX15KA4AALg8gjBOyff78nTLkk06VlIpSQoN8NHCq0/XjJQYJ1cGAADQOgRhtNlbG7P0P//ZpqoaI0nqGxWkl+eP08CeIU6uDAAAoPUIwmi1mlqjxz79SS9/u9e2berAaC28erQigvycWBkAAEDbEYTRKtbyKv3+9R+1cleObdt1k/vqgQuGyMebh+IAAED3QxCGQ/tyS3TDK99rd06JJMnHy6KHLx6mayYkObkyAACAU0cQRovWZeTqt8t+UGFZlSQpIshXz10zRpMGRDm5MgAAgPYhCKNZS1Iz9dCHO1RTW/dQXHLPEL08f6ySooKdXBkAAED7EYTRSFVNrR75aKeWpGbatp0xKEZPXzVaoQG+TqwMAACg4xCEYaegtFK3LvtB63Yfs237zfT+uu/cwfL2YpEMAADgPgjCsMk4WqQbX9mofcdKJUl+3l569NLTdPmYBCdXBgAA0PEIwpAkrdh1VL9/7UcVVVRLkqJD/PTCvDEak9TDyZUBAAB0DrdtAPvjjz9q3rx56tOnj/z9/RUTE6MJEybo7rvvVlFRkbPLcxnGGL28Zo9uWPy9LQQP6RWmD26fSggGAABuzS2D8HPPPadJkyZp0KBBWrt2raxWq7766it5e3vr73//u44dO+b4RTxARXWN7nt3q/78yU+qbwyhc4bF6p1bJik+ItC5xQEAAHQyt5sasXLlSt12221asGCB/vjHP9q2jxo1SosXL9bs2bPl60vng9ziCv126SZ9vy/ftu33Zw7UHWelyIuH4gAAgAewGGOMs4voSGPHjlVGRoaOHj0qPz+/dr+e1WpVeHi4CgsLFRYW1gEVOt9Ph6y68ZWNOlBQJkny9/HSgitG6qKRvZ1cGQAAQNdxqzvCO3fu1KZNm3Teeed1SAh2R1/sOKw739ys0soaSVJsmL9eunasRiREOLcwAACALuZWc4TXrl0rSUpKStKXX36pM844Q2FhYQoKCtKoUaP097//XdXV1U6u0jmMMXp2RYZuXrLJFoJHJoTrw9unEoIBAIBHcqsgnJ6eLkn64osvNGfOHN10003KzMxURkaGJk2apLvvvlsXX3yxampqmn2NiooKWa1Wu4/urryqRne8uVmPf7HLtm32yN568+ZJig0LcGJlAAAAzuNWQbiwsFCStHfvXj3++OO6+uqrFRkZqd69e+u5557T6NGj9emnn+pf//pXs6/x2GOPKTw83PaRmJjYVeV3iqPWcl35Yqo+2HzQtu2ecwbpH78apQBfbydWBgAA4FxuFYRPdM011zTaNn/+fEnSa6+91uxx999/vwoLC20fWVlZnVZjZ9uaXaDZC9dqS1aBJCnIz1svzBuj284YKIuFzhAAAMCzudXDcj161C0AER0draCgoEbj/fr1kySlpaU1+xr+/v7y9/fvnAK70EdbDurut7eoorpWkhQfEaiXrh2rob3do/MFAABAe7lVEB46dKgkqbKyssX93PluaG2t0VNfp+np5Rm2bWOTIvX8vDGKDun+AR8AAKCjuFUQPuuss+Tl5SWr1aq8vDzbHeIG+/btk3Q8MLub0spq/fHNLfp8x2HbtivGJOjPvxwufx/mAwMAAJzIreYI9+rVS7/61a8kSYsXL240/uqrr0qSbrzxxq4sq0scKCjT5c+tt4VgL4v0wAVD9LfLRxCCAQAAmuBWd4Ql6amnntKmTZv00EMPKT4+XrNnz1ZBQYH+7//+T5s2bdJNN92kK6+80tlldqhNmfm6eckm5RZXSJJC/X309NWjdcagnk6uDAAAwHW53RLLUl0btUcffVTvvfee9u/fr8DAQI0ePVo333yz7Y5xa7n6EsvvbMrWf7+3TZU1dQ/FJUUF6V/zx2pgz1AnVwYAAODa3DIIdyRXDcI1tUZ/+/xnvbB6j23b5AFRevbq0xUZzPLSAAAAjrjd1AhPUFRepT+8sVnLfz5q2zZ3Yh89eNEw+Xq71bRvAACATkMQ7mb2HyvVDa98r/SjxZIkby+LHrpoqOZN6uvcwgAAALoZgnA3sn73Md26bJPyS6skSeGBvvrnNadrysBoJ1cGAADQ/RCEu4nXNuzX//tgu6pr66Z0D4gJ1r/mj1Pf6GAnVwYAANA9EYRdXHVNrf78yU9avG6fbduMlBg9c/VohQX4Oq8wAACAbo4g7MIKS6t022s/6NuMXNu2G6f20/3nD5G3l/suEw0AANAVCMIuandOsW58ZaP25pZIkny9Lfq/S07TnHGJTq4MAADAPRCEXdDqtBzd9toPKiqvliRFBfvp+XljNK5vDydXBgAA4D4Iwi5myfp9evDDHap/Jk6D40L10rVjldgjyLmFAQAAuJkuWX1h9erVXfE23d6hwjL96YPjIXjW0Fi9+9vJhGAAAIBO0CV3hM844wzV1NR0xVt1awX1/YElaeagGL0wd4y8eCgOAACgU7Q7CFdWVuqrr77Srl27VFJSImNMR9Tl8XqFBxKCAQAAOlG7gvCOHTt0/vnnKzs7u8UAbLEQ6AAAAOBa2hWEb7/9dhUWFuqWW27R4MGDFRYW1mgfY4xuuOGG9rwNAAAA0OHaFYQ3bNigjz/+WGeeeWaL+/36179uz9sAAAAAHa5dXSMiIyM1YcIEh/vV1ta2520AAACADteuIDx//nytXLnS4X7cEQYAAICraVcQ/t///V+tWLFCCxcu1L59+1RRUdHkfq+88kp73gYAAADocO0Kwt7e3kpMTNRDDz2kAQMGKCgoSN7e3o0+AAAAAFfTroflFixYoHvvvVcWi0WRkZEKDQ1ttI8xRllZWe15GwAAAKDDtSsIP/fcc7r22mu1YMECRUdHN7ufl1eXrOQMAAAAtFq7gvDhw4f197//XVFRUS3u9+CDD7bnbQAAAIAO165btcOHD1dpaanD/a6//vr2vA0AAADQ4doVhB999FH96U9/Uk1NTYv79evXrz1vAwAAAHS4dk2NWLt2rQ4ePKiEhASdc8456tWrlwIDAzuqNgAAAKDTtCsIP/TQQ7b//+qrrza7n8Viac/bAAAAAB2uXUFYkhYtWtTiuDFGN9xwQ3vfBgAAAOhQ7Q7C8+fPd7gPSywDAADA1bTrYbn169e3ar+9e/e2520AAACADteuIDxhwoRW7ZeUlNSetwEAAAA6HEu+AQAAwCMRhAEAAOCRCMIAAADwSARhAAAAeCSCMAAAADwSQRgAAAAeqU1BeODAgbr33nu1bt26zqoHAAAA6BJtCsLx8fF68sknNW3aNPXu3Vu33nqrvv76a9XU1HRWfQAAAECnaFMQXrVqlQ4dOqQXXnhBo0eP1qJFi3TOOecoJiZG1157rd5//32VlZV1Vq0AAABAh2nzHOHo6GjdeOON+uSTT5STk6Nly5Zp1qxZev/993XppZcqJiZGl112mZYuXaqCgoJOKBkAAABov3Y9LBcSEqJf/epXevPNN5Wbm6sPP/xQc+bM0Zo1a3TttdcqNjZWZ599tp5//nkdPny4o2oGAAAA2q3Dukb4+fnpwgsv1L///W8dPnxYy5cv180336xdu3bp1ltvVUJCgqZMmdJRbwcAAAC0S6e0T/Py8tLMmTP19NNPKzMzU999953uvfde5efnd8bbAQAAAG3WJX2Ex44dq0cffVQ7d+7sircDAAAAHGJBDQAAAHgkgjAAAAA8EkEYAAAAHokgDAAAAI90ykF43bp1Ki4u7shaAAAAgC5zykF46tSpiomJ0UcffdSR9QAAAABdol1TIyoqKnTFFVe0GIa//vprffbZZ+15GwAAAKDDtSsIn3nmmerVq5euuOIKffzxx03us2zZMl144YXteRsAAACgw7UrCPfp00crV65UXFycLr/88mbDMAAAAOBq2t01IikpiTAMAACAbqdD2qf17dvXLgx/8sknHfGyAAAAQKfpsD7CJ4bhyy67jDDcTpXVtc4uAQAAwK116IIaJ98ZZppE2/j7HP/P8d6P2Xpx9W4ZY5xYEQAAgPs65SB82WWXqU+fPo22N4Th2NhYXXHFFUpNTW1XgZ6kX3SwLh7VW5JkjPTopz/r7re3qqK6xsmVAQAAuB+L6aRbjvv27dPMmTO1f/9+WSwW1dR0zzBntVoVHh6uwsJChYWFdfr71dYaPb08XU99nW7bNiYpUs/PHaOYUP9Of38AAABP0aFTI07Ut29frVq1StOmTVNQUFBnvY3b8fKy6I6zUvTs1acrwLfuP8+mzHxdvPBb7ThY6OTqAAAA3Een3RE+kTFGFouls9+mU3T1HeETbT9QqBtf2ajD1nJJUqCvt568cqTOHd6rS+sAAABwR512R/hE3TUEO9vw+HB9ePsUjUqMkCSVVdXolqU/6Jlv0nmIDgAAoJ26JAjj1PUMC9Abv5moX46Ot237+1dp+v0bm1Ve1T3nXQMAALgCgnA3EODrrSfmjNR95w5Ww831j7Yc1JwX1utwYblziwMAAOimCMLdhMVi0W9nDtCL88Yq2M9bkrQ1u1CzF36rzVkFzi0OAACgGyIIdzOzhsbq3VsnKyEyUJJ0tKhCV76wXh9sPuDkygAAALoXgnA3NDguTB/cNkXj+/aQJFVU1+oPb2zWgi92qbaWh+gAAABagyDcTUWF+GvpjRP0q3GJtm0LV2TolqWbVFJR7cTKAAAAugePCML79+9XWFiYLBaLVq5c6exyOoyfj5ceu/Q0/b8Lh8qr/iG6L3ce0WXPrVN2fqlziwMAAHBxHhGEb7zxRhUVFTm7jE5hsVj066n9tOj68QoN8JEk/Xy4SBcvXKvv9+U5uToAAADX5fZB+MUXX9TGjRs1ePBgZ5fSqWakxOg/t05R36i65ayPlVTq6pdS9dbGLCdXBgAA4JrcOgjv379f99xzj5588knFxsY6u5xON7BniN6/bYqmDIySJFXVGN37zlb9+eOdquEhOgAAADtuHYRvvPFGTZkyRfPnz3d2KV0mIshPi68fr/mTkmzbXv52r2545XtZy6ucWBkAAIBrcdsg/NJLL2nDhg164YUXnF1Kl/P19tLDFw/X//1yuHzqn6JbuStHl/5znfbllji5OgAAANfglkE4KytLd999txYsWKDExETHB5ygoqJCVqvV7qO7umZCkpbcMEERQb6SpIyjxbr42bVal5Hr5MoAAACczy2D8E033aRx48bppptuavOxjz32mMLDw20fbQ3SrmbSgCh9cNsUJfcMkSQVllVp3r+/05LUTCdXBgAA4FxuF4Rffvllffvtt3r55ZdP6fj7779fhYWFto+srO7fdSEpKljv3TpZZwyKkSTV1Br96f3t+tP721VVU+vk6gAAAJzDrYJwdna27rrrLv3lL39R3759T+k1/P39FRYWZvfhDkIDfPXy/HH6zfT+tm1LUjM1/9/fqaC00omVAQAAOIfFGOM2fbUWL16s66+/vtX7r1ixQjNnzmxxH6vVqvDwcBUWFrpNKH5nU7b++71tqqy/G9w3Kkgvzx+rgT1DnVwZAABA13GrINySmTNnatWqVa0KvydyxyAsSZsy83Tzkk3KLa67Gxzq76Onrx6tMwb1dHJlAAAAXcOtpkag9cYk9dAHt0/VkF514b6oolo3LP5eL6/ZIw/53QgAAHg4grAHi48I1Du3TNI5w+pW3as10p8/+Un3vbtVFdU1Tq4OAACgc7n11Ih9+/apX79+TY7NmDFDK1eudPga7jo14kS1tUZPfZ2mp5dn2LaN6xup5+aOUXSIvxMrAwAA6DxuHYQ7gicE4QYfbjmoe97eoorquofo4iMC9fL8sbbpEwAAAO6EqRGwmT2yt96+ZZJiw+ruAh8oKNNlz63TlzsOO7kyAACAjkcQhp0RCRH68PapGpkQLkkqrazRb5Zs0rMrMniIDgAAuBWCMBqJDQvQmzdP0uyRvW3bHv9il+54c7PKq3iIDgAAuAeCMJoU4Outf/xqlO45Z5Bt2webD+rKF1N11FruxMoAAAA6BkEYzbJYLLrtjIF6fu4YBfl5S5K2ZBVo9sK12pZd6OTqAAAA2ocgDIfOHR6nd26ZrPiIQEnSYWu5rnhhnT7actDJlQEAAJw6gjBaZWjvMH1w+xSNTYqUJJVX1ep3r/+oJ77cpdpaHqIDAADdD0EYrRYd4q9lN03QFWMSbNueXp6hW5f9oNLKaidWBgAA0HYEYbSJv4+3/nb5CD1wwRB5Weq2fb7jsC5/br0OFJQ5tzgAAIA2IAijzSwWi26c1l//mj9Oof4+kqSdh6y6eOFabcrMd3J1AAAArUMQxik7Y3BPvXfrZCVFBUmScosrdNWLqXpnU7aTKwMAAHCMIIx2SY4N1fu3TtHkAVGSpMqaWt399hY99ulPquEhOgAA4MIIwmi3yGA/vfLr8Zo7sY9t2wur9+imVzeqqLzKiZUBAAA0jyCMDuHr7aU/X3Ka/vfiYfKuf4pu+c9Hdek/12n/sVInVwcAANAYQRgdat6kvnr11+MVHugrSUo/WqyLn/1W63cfc3JlAAAA9gjC6HBTBkbr/dumaEBMsCQpv7RK8/61Qa9t2O/kygAAAI4jCKNT9IsO1n9um6IZKTGSpOpao//+zzY99OEOVdfUOrk6AAAAgjA6UViAr/593TjdOLWfbdvidft03aLvVVBa6cTKAAAAJIsxhh5XLbBarQoPD1dhYaHCwsKcXU639db3Wfqf97epqqbucgv289YvT4/X3IlJGhzHeQUAAF2PIOwAQbjjfL8vTzcv2aS8Evu7weP79tDcSUk6d1ic/Hz4IwUAAOgaBGEHCMId62BBmZ5ZnqEPNh9QaWWN3Vh0iJ9+Na6PrprQR/ERgU6qEAAAeAqCsAME4c5hLa/Sf344oCWpmco4Wmw35mWRfjEkVvMmJmnqwGh51fclBgAA6EgEYQcIwp3LGKPUPXlampqpL3YcVvVJyzL3jQrS3IlJunxMgiKC/JxUJQAAcEcEYQcIwl3niLVcb3yXpde+y9QRa4XdmL+Pl2aP7K15k5I0IiHCOQUCAAC3QhB2gCDc9apqavXNT0e0JDVTazMar0g3MiFccycm6aKRvRXg6+2ECgEAgDsgCDtAEHaujKPFWrYhU+9sylZRebXdWHigr+aMTdA1E5LUNzrYSRUCAIDuiiDsAEHYNZRWVuvDzQf16vpM7TxkbTQ+PSVG8yYm6czBPeXNw3UAAKAVCMIOEIRdizFGP2YVaOn6TH289ZAqT1quOT4iUFdP6KM5YxMVE+rvpCoBAEB3QBB2gCDsuvJKKvXWxiwt25CprLwyuzFfb4vOG95LcycmaVzfSFks3CUGAAD2CMIOEIRdX02t0eq0HC1JzdSKXUd18hU9OC5Ucycm6ZLR8Qrx93FOkQAAwOUQhB0gCHcvWXmlWrZhv97amNVoKecQfx9denq85k5MUkpsqJMqBAAAroIg7ABBuHuqqK7RZ9sO69X1+/TD/oJG4xP69dDciUk6Z1ic/Hy8ur5AAADgdARhBwjC3d+Og4Vamrpf7/94QGVVNXZj0SH+ump8oq4a30e9IwKdVCEAAHAGgrADBGH3YS2v0nubsrUkNVO7c0rsxrws0llDYnXtpL6aPCBKXrRgAwDA7RGEHSAIux9jjNbvOaalqZn6YscR1dTafwn0iw7WNRP66IoxiQoP8nVSlQAAoLMRhB0gCLu3w4XleuP7/Xr9u/06Yq2wGwvw9dLskb01b2JfnZYQ7qQKAQBAZyEIO0AQ9gxVNbX6eucRLUnN1LrdxxqNj0yM0LyJSbpwRC8F+Ho7oUIAANDRCMIOEIQ9T8bRIi1N3a93N2WrqKLabiwiyFdzxibqmgl9lBQV7KQKAQBARyAIO0AQ9lylldX6YPNBvbo+Uz8dsjYan5ESo3kTk3TG4J7y5uE6AAC6HYKwAwRhGGP0w/4CLU3N1CdbD6myptZuPD4iUFdP6KMrxyUqOsTfSVUCAIC2Igg7QBDGiY4VV+itjdlatiFT2flldmO+3hadf1ovzZuYpDFJkbJYuEsMAIArIwg7QBBGU2pqjValHdWS9ZlamZajk7+KBseFat6kJF0yKl7B/j7OKRIAALSIIOwAQRiO7D9WqmXfZeqt77OUX1plNxbi76PLTo/X3IlJSo4NdVKFAACgKQRhBwjCaK3yqhp9uu2QlqRm6sf9BY3GJ/bvoXkT++rsYbHy9fbq+gIBAIAdgrADBGGciu0HCrU0NVPvbz6g8ir7h+t6hvrrV+P76KrxieoVHuikCgEAAEHYAYIw2qOwrErvbsrW0tRM7cktsRvz9rJo1pBYzZuUpMkDoni4DgCALkYQdoAgjI5gjNG63ce0NDVTX+48oppa+y+7/jHBmjshSZeNSVB4oK+TqgQAwLMQhB0gCKOjHS4s1+vf7dfr3+3X0aIKu7EAXy9dMqru4brh8eFOqhAAAM9AEHaAIIzOUlVTq692HtGS9Zlav+dYo/HRfSI0b2KSzj+tlwJ8vZ1QIQAA7o0g7ABBGF0h/UiRlm3Yr3c3ZauootpuLDLIV3PGJeqa8UnqExXkpAoBAHA/BGEHCMLoSiUV1fpg80G9un6ffj5cZDdmsUgzUmJ07aQkzUjpKW8vHq4DAKA9CMIOEIThDMYY/bA/X0vWZ+rTbYdVWWPfgi0+IlDXTOyjK8cmKirE30lVAgDQvRGEHSAIw9lyiyv01sYsLUvdrwMFZXZjft5eOv+0OM2blKTT+0TSgg0AgDYgCDtAEIarqKk1WrnrqJakZmpVWo5O/sod0itM8yYm6eJRvRXs7+OcIgEA6EYIwg4QhOGKMo+V6LUN+/XmxiwVlFbZjYX6++iyMQmaO7GPBvYMdVKFAAC4PoKwAwRhuLLyqhp9svWQlqRmanNWQaPxSf2jNG9SkmYNjZWvt1fXFwgAgAsjCDtAEEZ3sS27UEtTM/XBlgMqr7J/uK5nqL+uGt9HV43vo7jwACdVCACAayEIO0AQRndTWFqld37I1tLUTO3NLbEb8/ay6OyhsZo3MUmTBkTxcB0AwKMRhB0gCKO7qq01Wrf7mJak7tNXO4+o9qSv9AExwZo7MUmXnp6g8EBf5xQJAIATEYQdIAjDHRwqLNPrG/br9e+zlFNUYTcW6OutS0b31tyJSRrWO9xJFQIA0PUIwg4QhOFOKqtr9eXOw1qyPlMb9uY1Gj+9T4TmTUrSecN7KcDX2wkVAgDQdQjCDhCE4a7SjhRpaWqm3vvhgIorqu3GegT7ac7YRF0zoY8SewQ5qUIAADoXQdgBgjDcXXFFtT7YfEBL1mfq58NFdmMWi3TGoJ6aNzFJ01Ni5O3Fw3UAAPdBEHaAIAxPYYzRxsx8LU3N1KfbDqmqxv5bQ2KPQF0zIUlzxiaqR7Cfk6oEAKDjEIQdIAjDE+UUVeitjVl6bcN+HSgosxvz8/HS+cPjdM6wOE0eGE3HCQBAt0UQdoAgDE9WU2u04uejWpKaqVVpOY3GvSzSqMQITU+J0fSUGI1MiGD6BACg2yAIO0AQBursyy3Ra9/t11sbs1RQWtXkPmEBPpqaHK3pyTGalhKj+IjALq4SAIDWIwg7QBAG7JVX1Wj9nmNanZajNem5yjha3Oy+A2KC6+4WJ8doQv8eCvLz6cJKAQBomdsF4eXLl2vZsmVas2aN9u/fLy8vL/Xt21cXXHCB7r33XsXExLTp9QjCQMsOFJRpTVqOVqfn6Nv0XFnLq5vcz8/bS+P6RWpacl0wHtIrlCWeAQBO5VZB+Nlnn9Xtt9+uESNG6Mknn9S4ceNktVq1bNkyPfDAA+rRo4dWr16tlJSUVr8mQRhovZpaoy3ZBVqTlqvV6Tn6cX9+o6WdG0SH+Gt6crSmp8RoanK0okP8u7ZYAIDHc6sgvGDBAv3P//yPMjIylJiYaDf2X//1X/rrX/+qc889V5999lmrX5MgDJy6wrIqrd+dq1VpuVqdltOoA8WJhvUO0/SUGE1LjtbYpB7y8/HqwkoBAJ7IrYLw0qVL9fXXX2vx4sWNxtatW6cpU6bIx8dH5eXl8vZu3fKxBGGgYxhjtDe3RKvTcrQ6PVfrdx9TWVVNk/sG+XlrUv8oTau/Y9wvOphpFACADudWQbglO3bs0PDhw+Xr66vS0lL5+LTuoR2CMNA5KqprtCkzX6vTcrUmPUc7Dlqb3TchMlDTkmM0IyVakwdGKyyA3sUAgPbzmCD87rvv6vLLL9cZZ5yh5cuXt/o4gjDQNXKKKrQ2I9d2xzi3uKLJ/by9LHW9i5NjND0lWiPoXQwAOEUeE4QvvvhiffTRR1q+fLlmzpzZ7H4VFRWqqDj+A9hqtSoxMZEgDHSh2lqjnw8XaXV6jlan5WjjvnxV1tQ2uW94oK+mDozW9JRoTUuOUW96FwMAWskjgvDnn3+u8847T3fddZcWLFjQ4r4PPfSQHn744UbbCcKA85RWVmvDnjytSsvRmvQc7c4paXbfgT1D6hf0iNbEflEK9Gvd8wAAAM/j9kE4PT1dU6dO1bRp0/Tmm286fEiOO8KA68vOL9Wa9Lq5xa3pXVw3jSJGg+PoXQwAOM6tg3BmZqZmzJihMWPG6PXXX5efn1+bX4M5woBrq66p1ZbsQq2pn0axOaug2d7FMaH+mpYcrRkpMZoykN7FAODp3DYIZ2Rk6Be/+IWmT5+uxYsXt7pd2skIwkD3UlhapXW7c+vnF+e22Lt4eHxY3TSK5BiNSYqkdzEAeBi3DMI7d+7UWWedpQsuuEAvvPCCvLyO/3BbsGCBrrzyykYLbjSHIAx0X8YY7anvXbzGQe/iYD9vTRoQVbcEdEqM+kYFMY0CANyc2wXhLVu2aNasWbryyiv19NNPN/pBZrFYtGLFihY7R5yIIAy4j4rqGm3al6/V6XVt2nYear53cWKPut7F05NjNHlgFL2LAcANuVUQ3rhxo84++2xVVFTooosuanKfN998kyAMQJJ0tKi8vndx3YN3ucWVTe7n7WXR6MQITU+pu1t8Wnw4vYsBwA24VRC+44479I9//MPhfgRhACerrTX66bBVq9Pq7hZvzMxTVU3T3x4jgnw1ZWC0ZtS3aesVTu9iAOiO3CoIdwaCMOCZSiurlbrnWF0wTs/RnhZ6Fyf3DKmfWxytCfQuBoBugyDsAEEYgHS8d/HqtBx9m5GrouZ6F/t4aUK/HpqWHK3pKTEaFEvvYgBwVQRhBwjCAE7W0Lt4dVqOVqfnaEsLvYt7hvrb7hZPHRitKHoXA4DLIAg7QBAG4EhhaZXW7q67W7w6LUcHC8ub3M9ikYb3Dtf0lGhNS47R6X3oXQwAzkQQdoAgDKAtjDHanVNiW+kudU+eg97F0ZqeEq3pyTHqGx3cxdUCgGcjCDtAEAbQHg29i1fVr3T3Uwu9i/v0CLLNLZ48IEqh9C4GgE5FEHaAIAygIx0tKte39Q/drUnP1bGS5nsXn94nQtPrV7obTu9iAOhwBGEHCMIAOkttrdHOQ1atTs/RmrRch72Lpw6MtgXjuPCALq4WANwPQdgBgjCArlJSUde7uKFN257c5nsXp8SGaHpyjKalxGhCvx4K8KV3MQC0FUHYAYIwAGfJyjveu3jtbse9ixvuFqfEhtC7GABagSDsAEEYgCuo611coFVpuVrjoHdxbFhD7+IYTR0YrR7Bfl1bLAB0EwRhBwjCAFxRQWml1mYcs7Vpa6l38Wnx4XXTKJKjdXpSpHy96V0MABJB2CGCMABXV9e7uFir03K1Oj1HqXuOqbyqtsl9Q/x9NGlAlKbXt2lLiqJ3MQDPRRB2gCAMoLspr6rRpsz8+iWgW+5dnBRV37s4OUaT6F0MwMMQhB0gCAPo7o5ay7UmvW5ucUu9i328LDq9T2TdSncpMRreO1xe9C4G4MYIwg4QhAG4kxN7F69Oy9GmzPxmexdHBvlqav3c4unJ9C4G4H4Iwg4QhAG4s4bexQ0r3bXUu3hQbKhtCejx9C4G4AYIwg4QhAF4kqy8Utvd4nUZx1RU0XTvYn8fL43v10MzUuratCX3pHcxgO6HIOwAQRiAp6quqdXmrAKtrl/UY2t2872L48ICbHeLpw6MViS9iwF0AwRhBwjCAFCnoXdxXTeKHB1qoXfxiPhwTU+J0bTkGI3uE0HvYgAuiSDsAEEYABpr6F28Kq3ubvGGvc33Lg6t7108LSVGM5Jj1CcqqIurBYCmEYQdIAgDgGPlVTXauC/fNr/458NFze6bFBWk6fVLQE8aEKUQf58urBQAjiMIO0AQBoC2O2ot1+oTehfntdS7OClSM1Lq2rTRuxhAVyIIO0AQBoD2qa012nHQvndxdTNP3fUI9tPUgdG2B+9iw+hdDKDzEIQdIAgDQMcqrqhW6u5jWl1/t3hvC72LB8cd7108ri+9iwF0LIKwAwRhAOhcWXmlWpWWozXpjnsXT+gfpenJ0ZqREqOB9C4G0E4EYQcIwgDQdarqexevScvRqvRcbc0uUHM/pXqFH+9dPGUAvYsBtB1B2AGCMAA4T35JpdburmvRtjotV4etLfQuTojQjORoTUuJ0ahEehcDcIwg7ABBGABcgzFGGUeL66dR5Cp1zzFVVDfduzjE30cjEsI1IiFCIxPCNSIxQr3DA5hKAcAOQdgBgjAAuKbyqhp9vy9Pa+qXgG6pd7EkRYf414fjcI1MiNCIhHBFhfh3UbUAXBFB2AGCMAB0D0es5bZQvGHvMR2xVjg8Jj4iUCMT6+4cj0gI12nx4QoN8O2CagG4AoKwAwRhAOiejljLtSWrQFuzC7Ulu+5/C8uqWjzGYpH6Rwfb7hiPSIzQ0F5htG0D3BRB2AGCMAC4B2OM9ueVakt2obbWB+TtBwtVWlnT4nE+XhYNigs9Pt84IUIpsSHy4WE8oNsjCDtAEAYA91VTW/cAXt0d47pw/NMhq6pqWv7RGODrpWG97ecb940KZnlooJshCDtAEAYAz1JRXaOfDxVpa3ZB3d3j7AKlHy1utp9xg9CAkzpVJESoF50qAJdGEHaAIAwAKKmo1vYDhXbzjffnlTo8LjrE3xaKRyTW3T3uwcIfgMsgCDtAEAYANCW/pFJbD9TNN264c3y0yHGnioTIwOMP4yVE6LSEcIX4+3RBxQBORhB2gCAMAGitw4XldvONW9upYkBMiN184yF0qgC6BEHYAYIwAOBUGWOUeazUNp1ia3aBth+wqqyq5U4Vvt6NO1Uk96RTBdDRCMIOEIQBAB2puqZWGTnF2pp1fL7xz4cdd6oI9PXWsN5hdeG4fhGQvlFBPIwHtANB2AGCMACgs1VU1+inhk4VWXV3jjNyHHeqCAvwsa2K1xCQ48LoVAG0FkHYAYIwAMAZim2dKo4/jJeVV+bwuJjQEzpV1M87jqRTBdAkgrADBGEAgKvIK6k84UG8uoCc04pOFYk9Au3mGw+Pp1MFIBGEHSIIAwBclTFGh63ltukUDQHZWl7d4nEWizQwJsRuvvGQXqHy96FTBTwLQdgBgjAAoDsxxmjfsVK7+cbbDxaqvKq2xeN8vS0aHBd2vI1bYrgGxtCpAu6NIOwAQRgA0N1V19Qq/Wix3Xzjnw8VqbqWThXwbARhBwjCAAB3VF5Vo58OWe2Wjd5Npwp4GIKwAwRhAICnKK6o1rb6O8ZbD9CpAu6PIOwAQRgA4MkaOlU0zDfekl2o3OK2d6o4LT5cwXSqgIshCDtAEAYA4DhjjA4VlttC8TY6VaAbIwg7QBAGAKBltbVGmXkd06kiuWeovL2Yb4yuQRB2gCAMAEDbtadTxfD4MLv5xkl0qkAnIQg7QBAGAKBjnGqnivBA3/ouFfWdKhIiFBce0DVFw60RhB0gCAMA0HmKyqu0/YDVtjLeluwCZec77lTRM9T/+MN4iREaER9Opwq0GUHYAYIwAABd61hxRV37tjZ2qujTI+j4fOOEcA2nUwUcIAg7QBAGAMC5Tu5U0XD3uMhBpwovizSwZ4hdG7fBdKrACQjCDhCEAQBwPbW1RvuOldjNN97Ryk4VQ3qF2c03HtgzhE4VHoog7ABBGACA7qG6plZpR+w7Vew67LhTRZCft4b3rn8YL7Hu7nGfHnSq8AQEYQcIwgAAdF/lVTXaeciqrVnHH8bbk1visFNFRJCvTos/Pt94ZGKEYsPoVOFuCMIOEIQBAHAvReVV2nagUFvr7xpvySrUgQLHnSpiw/zt5huPSAhXRBCdKrozgrADBGEAANxfbnGFtp0w33hrdoFyiysdHpcUFWQXjofHhynIj04V3QVB2AGCMAAAnscYo4OF5dqadXy+8bbsQhVVOO5Ukdwz1G4BEDpVuC6CsAMEYQAAINV1qth7rMQ2nWJrdoF2HLSqorrlThV+3l4a3CuUThUuiCDsAEEYAAA0p6qmVmlHiuzmG+86UqQaOlV0CwRhBwjCAACgLcqrarTjoP2y0XtyShweR6eKrkcQdoAgDAAA2staXqXt2YV2K+PRqcL5CMIOEIQBAEBnyC2usJtvvDW7UMdK6FTRlQjCDhCEAQBAVzDG6EBB2fFlo7MKte1AoYrb0qmifr7x4Lgw+fl4dVHl3RdB2AGCMAAAcJbaWqM9uSV28413HLSqshWdKob0CrVNpxiZGKEBMXSqOJlbBuHy8nL97W9/07Jly5SZmanIyEidddZZeuSRR9SvX782vRZBGAAAuJKqmlrtOnxCp4rsQqW1olNFsJ+3hsWH26ZUjEyIUGKPQI/uVOF2QbisrExnnXWWtm7dqn//+9+64IILlJ6ervnz52vfvn1atWqVRo4c2erXIwgDAABXV1ZZo52HCu3mG+/JddypIjLIV6edMN94ZEK4enpQpwq3C8L33Xef/va3v2nhwoW67bbbbNt3796tQYMGafDgwdqyZYu8vVu3wgtBGAAAdEeFZVXafuD4fOOt2QU6WFju8Li4sADbdIoRCeEaER+h8CDfLqi467lVEC4tLVXPnj1VVVWlY8eOKSQkxG783HPP1RdffKEPP/xQF110UatekyAMAADcRU5RhW06RcOd47xWdKroW9+poiEgD+vtHp0quv9ncILly5erpKRE48aNaxSCJWnq1KltDsIAAADuIibUX78YEqtfDImVVNepIju/7IT5xgXafsDaqFPFvmOl2nesVB9uOSiprlNFSqz9stGD4kK7XacKtwrCW7ZskaRmH4jr37+/3X4AAACezGKxKLFHkBJ7BOmCEb0kNXSqKLbNN96SXaidh+w7VdQa6efDRfr5cJHe2pgtqb5TRe8wu/nG/V28U4VbBeHDhw9Lknr06NHkeGRkpN1+TamoqFBFRYXt31artQMrBAAAcG1eXhYN7BmqgT1DddmYBElSZXWt0o4U2eYbb8kuUPrRYrtOFZU1tdqSVaAtWQWSMiVJf7n0NP1qfB8nfBat41ZBuLS0VJLk59f00oP+/v52+zXlscce08MPP9zxxQEAAHRTfj5eGh4fruHx4bpmQt22ssoa7Th4fNnobU10qjgtIdwJ1baeWwXhoKAgSVJlZdOTvhvu9Dbs15T7779ff/zjH23/tlqtSkxM7MAqAQAAur9AP2+N7dtDY/se/0v8iZ0qdhywKiU21IkVOuZWQTguLk6SlJeX1+R4fn6+JKlXr17Nvoa/v7/tzjEAAABaLzzQV1MGRmvKwGhnl9Iq3evRPgdGjRolSdq7d2+T43v27JEkjRgxoqtKAgAAgItyqyB8xhlnKDg4WFu3blVJSePVVNauXStJmj17dleXBgAAABfjVkE4KChIt99+uyoqKvTKK6/Yje3Zs0dff/21hg0bpvPPP99JFQIAAMBVuNXKcpJUVlamX/ziF9q+fbsWLVqk888/X+np6Zo/f7727Nmj1atXa+TIka1+PVaWAwAAcE9udUdYkgIDA7V8+XLddddduv/++xUZGalZs2Zp2LBh+vHHH9sUggEAAOC+3O6OcEfjjjAAAIB7crs7wgAAAEBrEIQBAADgkQjCAAAA8EgEYQAAAHgkgjAAAAA8EkEYAAAAHokgDAAAAI9EEAYAAIBHIggDAADAIxGEAQAA4JEIwgAAAPBIPs4uwNUZYyRJVqvVyZUAAACgJaGhobJYLK3enyDsQFFRkSQpMTHRyZUAAACgJYWFhQoLC2v1/hbTcMsTTaqtrdXBgweb/Q3DarUqMTFRWVlZbTrxnoLz0zLOT8s4Py3j/LSM8+MY56hlnJ+WueL54Y5wB/Py8lJCQoLD/cLCwlzmInBFnJ+WcX5axvlpGeenZZwfxzhHLeP8tKw7nx8elgMAAIBHIggDAADAIxGE28nf318PPvig/P39nV2KS+L8tIzz0zLOT8s4Py3j/DjGOWoZ56dl7nB+eFgOAAAAHok7wgAAAPBIBGEAAAB4JIIwAAAAPBJBGAAAAB7JY4NwamqqLrzwQkVFRSk4OFhjxozRiy++qLY+O1hWVqZ3331X119/vUaMGKGIiAiFhYVpyJAh+v3vf689e/Y0e2x5ebkeeeQRDRo0SAEBAerVq5fmzZunvXv3tvien332mc4880yFh4crLCxMU6dO1TvvvNOmuh3pqPNTXFys559/Xueee65iY2Pl6+urHj166IwzztArr7zS5OstXrxYFovF4cd1111nd9xDDz3Upv3bw5nnR5JmzpzZ4ue6ePHiJo+rra3V008/rZEjRyooKEjR0dH65S9/qc2bN7fxDDjmzHNUW1urDz74QFdffbX69esnf39/hYSEaPTo0XrkkUdUXFzc5Ht50jUk1S1Fetddd6lfv34KCAhQnz59dNtttyknJ6fF91y6dKkmTpyokJAQRUZG6uyzz9aKFSvaVLcjHXV+TrRr1y7NmDFDFotFK1eubHa/lStXtup70MyZM+2Oc/S96+T928OZ50eSrrvuuhY/14ceeqjZY939+pGk5cuX64YbblBKSooCAgIUFBSkoUOH6p577mn268uTrh/JhTKQ8UDvvvuu8fb2NhdffLHZt2+fsVqt5q9//auxWCzmuuuua9NrffbZZ0aSGT16tFm9erUpLi42R44cMf/85z+Nr6+vCQkJMZs3b250XGlpqZk8ebIJCQkxb731likpKTGbN282I0eONOHh4U0eY4wxTz75pJFkfvOb35jDhw+b3Nxcc+eddxpJ5sEHHzyV09FIR56fsWPHGknmlltuMWlpabbP89xzzzWSzGWXXWZqa2vtjlm0aJHx8fExgwYNavKjb9++RpL5y1/+Ynfcgw8+aMLCwpo97r/+67/afW6Mcf75McaYGTNmmLi4uGY/1/fee6/RMTU1Neayyy4zPj4+5p///KcpLi426enp5he/+IXx8/MzX3zxxSmfk5M5+xzdfffdRpKZMWOGSU1NNcXFxWbv3r3mnnvuMZLMoEGDzNGjRxu9lyddQ0ePHjXJyckmLi7OfP3116a0tNSsWbPG9OnTxyQkJJisrKwm3++OO+6wfb/Jz8832dnZ5qqrrjIWi8X8+9//PqXzcbKOPD/GGFNRUWEeeeQRExwcbEJDQ40ks2LFimb3X7Fihe06aeojOTnZds5PtGjRIhMYGNjscfPmzWtz7U1x9vkxxpj58+ebqKioZj/XZ555psnjPOH6WbhwoZFkRowYYb755htjtVpNdna2+etf/2p8fX1NbGys2bVrV6PjPOn6caUM5HFB+MiRIyY0NNQkJiaa0tJSu7H58+cbSeaNN95o9et99tlnxmKxmH379jUau/32240kc+211zYau/fee40ks3DhQrvtGRkZxtvb2wwbNsxUV1fbjW3bts14e3ubcePGNfrBNnPmTCPJrF+/vtW1N6Wjz8+wYcPMueee22h7WVmZLdCe/HqLFi0ySUlJzb7mk08+aQICAkxOTo7d9gcffNDMnz+/1bWdClc4P8bUBeFFixa1qfZ//vOfRpK5++677bbn5eWZyMhI07NnT2O1Wtv0mk1xhXN02223mdjYWFNYWNjouCuvvLLJEGOMZ11Dc+bMMZLMxx9/bLe9IQQ29ZqffvqpkWQuv/xyu+2VlZUmOTnZBAQENPm9sC06+vwYY8w555xjxowZY3744QczY8aMVgfh5vznP/8xksy2bdvsti9atMjMmDGjTbW1lSucn4b3amvw8JTr5/HHHzd+fn5m//79jcbuu+++Zr++POn6caUM5HFTI55//nkVFRVp3rx5CgwMtBu75ZZbJEl/+9vfWv1648aN08qVK5WUlNRobNCgQZKkgoICu+2lpaV69tln5efnp/nz59uNDRgwQGeddZZ27NihTz/91G7siSeeUE1NjW666SZZLBa7sZtvvlmS9Pjjj7e69qZ09PkZPny4brjhhkbbAwICNGvWLEnSF198YTfWr18/XXTRRU2+njFGzz77rK666ipFR0e3uo6O4grn51QtWLBAkvSb3/zGbntkZKTmzJmjo0ePNjuloi1c4RwlJSXp+uuvV1hYWKPjZs+e3eQxXcUVzk9mZqbefvttJSQk6Pzzz7cbmzlzpgYNGqTPP/9c27Ztsxtr+P5y8jXk6+urX//61yovL9fTTz/d6tqb0tHnp+G4DRs2aPTo0a3av2fPnrryyiubHX/66ad1xhlnaPjw4W2qoyO4wvk5VZ5y/cTFxemqq65SYmJio7GG7z9ff/21ampq2lRHR3CF8+NqGcjjgvCHH34oSZoyZUqjsbFjxyogIEA//PCDDhw40KrXi4qK0vTp05scW7dunSTp7LPPttu+fPlylZSUaOTIkQoJCWl03NSpU+1qbfDRRx81W3vDMZ9++qmqq6tbVXtTOvr8vPHGG7r88subHAsPD29y+4wZM/TMM880Ofbpp58qIyNDv/vd71r1/h3NFc7Pqdi+fbv27NmjmJgYJScnNxpv7po7Fa5wju655x499thjbTqmq7jC+fnoo49kjNGkSZMa/UCRmr4erFarVq1aJYvFokmTJrXqmFPR0edHki655BJ5e3u3ev+hQ4fqjTfeaHJs+/btWrFihdt8D5Lafn5OhSddP3Pnzm32pkLD16TFYmnXfNxT5Qrnx9UykEcF4ZqaGu3YsUNS3V3Hk/n4+Nh+g9uyZcspvUdZWZl++ukn/fGPf9Sbb76pP/zhD/rtb39rt0/DazdVgyT179+/UQ0HDhxQbm5us8fFx8fLz89P5eXlSktLO6Xau+L8nOjnn3+WpDY9APDMM89o6tSpzf7muXfvXl199dVKTEyUn5+fYmNjNXv2bH3zzTftrtfVzs/KlSt11llnqWfPnvL391e/fv104403Nvnf/1SuuVPhaufoVI7xhGvoVK6Hbdu2qba2VjExMU3+8Go4Zvfu3SopKTmlerv6/JyKhQsXKikpyXZn72S5ubm65ZZbbA9pRkVFadasWXr77bfb/d6udn42b96s2bNnq3fv3vL391d8fLx+9atfaePGjY325fqp0/A1OXXqVPn4+DQa94Trx9UykEcF4YKCApWXl0uSevTo0eQ+kZGRkqTDhw+3+fWXLl1qezL0zTff1KJFi7RgwQJ5edmf5obXbksNDf8/ICCg0Z8zpLrfLiMiIk65dqnzz8+Jjhw5oi+++EIpKSm66qqrWnXMrl279OWXX7Z4J+a7777T2LFj9d1336m4uFiffPKJ8vPzddZZZ+nPf/5zu2p2tfOzYsUK3XLLLcrIyFBeXp7+8Y9/6PPPP9fo0aP15Zdf2u3b2msuJydHtbW1p1y3q52jkxljtGTJEvn7++u+++5rch9PuIba8z3I0THGGB05cuSUau7K83MqCgoKtHTpUt16663N3gH76aefFB0drZUrV6qoqEhr1qxRSEiI5syZ02hKwKm8vyudnzVr1ujiiy/Wtm3bVFhYqKVLl2r79u2aNGmSXnnlFbt9uX7qvPrqq7JYLPp//+//NTnuCdePq2Wgxr+OuLhf/vKX+umnn9p0zKuvvqrx48ertLTUts3Pz6/Jff39/SXJbt/Wmjt3rq6++mrt27dPH330kX73u9/pmWee0euvv66BAwfa9mt47bbU4OiYE4+79tprm/yNuzlddX5O9Ic//EEWi0WvvfaafH19W3XMwoUL1bt3b1166aVNjl999dW6/PLL7ebtjR07Vh9//LH69eunP/3pT/ryyy919OjRNtXqiufnscceU9++fdWrVy/bttmzZysiIkIzZszQ1VdfrT179tjmyLb2mpOkiy++WOnp6W2q1xXPUVNeeOEFbdmyRc8884xtDv+JPOUa6ozvQSdeQ7/4xS/s/u2IM87PqfjXv/6l2tpa3XjjjU2On3XWWVq/fr3Gjx9v2zZ06FC9/fbbGjp0qF566SX9+OOPKioqatP7uuL5ufPOO/WnP/1JAwYMsG0744wz9MEHH2jw4MH6zW9+o+nTp9vu3nH9SJ9//rk+/PBD3XXXXU3+RcpTrp/OzkBtrb3bBeG9e/dq165dbTqm4aQEBQXZtlVWVja5b0VFRaN928LLy0v9+/fXH/7wB8XHx+uKK67Qeeedp507d9p+GDW8dltqcHTMice1ZW6P1LXnR5KeeuopvfPOO3rzzTc1ZsyYVh1jtVr1yiuv6N57723yz0mSlJKS0uT28PBwXXPNNVq4cKF27NihvLy8NtXriuenqTl2kjR9+nQNHz5c27dv14cffqi5c+fa1eOobknKyspy6a8x6dSuoXXr1umuu+7SHXfcodtvv73JfTzlGuqM70EnXkP79u1rU71dfX5ORW1trZ599lldffXVzd7JSkhIUEJCQqPtPj4+uvnmm3X33Xfr559/braPdXNc8fyMHDmyye0DBgzQueeeq48//ljLli3TAw88YFePp14/6enpmj9/vi677DL99a9/bXIfT7l+OjsDtbX2bheE29P4PyIiQgEBASovL1deXp5iY2Mb7ZOfny9JdnfaTtXll1+u6OhoZWRk6KuvvrI9nR0XFydJzf4wbaqGhmPKy8tVVlbW6E8Dxhhbd4qdO3dqyJAhba63K87PokWLdN9992nJkiW67LLLWn3c4sWLVVFRccp/Gmq4Iz948GCtXbv2lF7Dlc/PiQYOHKjt27fb/eWktddcz549Xf5r7FTO0Y8//qgLL7xQN9xwg5544olTel93uoba8z3I0TEWi0VFRUUKDg5uc+1d/T26LT7++GPt3bv3lB+Sa7h+wsLC2nxHr4Ern58TNXyup/I9yB2vn8zMTM2aNUtTp07Va6+9dkoPJrrT9dPZGaittXvUHGFvb2/bnzybWrmkurpaWVlZkqQRI0Z0yHs2/FnoxN9wR40a1WwNkmyr0Z1YQ3x8vK1dWFPHHThwQJWVlQoMDGyyK0BrdPb5ef7553Xrrbfq7bffbvWcTqnuAl+4cKGuvPJK9ezZs83vK8k277WpJ+Rby1XPz8ma+lxP5Zo7Fa54jjZs2KAzzzxTt9xyi55++ulTvgbc6Ro6levhtNNOk5eXl3Jycpp8mKnhmAEDBpxSiJGc8z26tZ5++mlNmzat2TuhjnSH66ejNPW5eur1k5GRoenTp2vatGl66623WvzTfkvc6fpxtQzkUUFYkq0/bUNrsxNt2rRJ5eXlOv300xUfH9+q1zv77LN1xx13NDt+8OBBSXVt1hqcccYZCg4O1tatW5v8htBwt+nkp5Jbqr3hmPPOO6/ZqQOt0dHnp8GTTz6pu+66S++//77d57Vjxw699NJLLR77+eefKz09vcU7MVlZWYqLi9PWrVubHG+Y8zp48OA21X0yVzg/b775Zotf6E19rsOGDVP//v2Vk5OjjIyMRsc0d82dClc4Rw1Wr16tWbNm6c4779Sjjz5q2261WrVgwQJZrVbbNk+6hi688EJZLBatX7++yRZOTV0PYWFhmjFjhowxWr9+fauOORWddX7aY+fOnfrmm28c3g2Oi4vTJ5980uSYq18/bbFu3TrFxcXZ7tydrKnP1ROvn507d2r69Ok6++yz9corr9jdCV6wYIEtdDbwlOvH5TJQm5fg6OYaVlXp06ePKSsrsxu7/vrrm11V5d577zVTpkwx27dvt9s+Y8YMk5KSYsrLyxsd8+WXXxpJJigoqNEqaA2ryzz77LN223fv3m18fHyaXFVl+/btxtvb20yYMKHRqipnnnlmh64s11Hnxxhj/vznP5uQkJAmV5pxtIqcMcace+65ZsKECS3us3fvXiPJPPbYY43G8vLyTHh4uJFk1qxZ0+LrOOIK52fRokXN/rf+5ptvjCTTo0ePRquqPffcc0aSueeee+y25+fnmx49enT4ynLOvoa+/PJLExQUZB5//PFGYw3Xy969extt84RryJjjK+x98skndttXrVrV7MpXDUvKX3HFFXbbKysrzaBBgzp0ZbCOPD8nau3KVye65ZZbTHx8vKmqqmpxP0nm5ptvbrS9srLS9O/f30gyS5YsafX7NsUVzk/Dynuvv/56o7Fdu3YZb29v4+vra/bs2WM35knXz+bNm01MTIy5/fbbm1zivKnX8JTrxxjXykAeF4SNMeadd94x3t7e5pe//KXJzMw0VqvVPP7448ZisTS5vGpOTo6RZCSZ22+/3W6sYVm/WbNmmXXr1hmr1WoOHTpkFi1aZKKjo423t7d55ZVXGr1maWmpmTRpkgkNDTXvvPOOKS0tNVu2bDGjRo0yYWFhDtfZvuWWW8yRI0fMsWPHzF133dWudbZP1pHn54EHHjCSzLBhw8yVV17Z6GP8+PEtBuG0tDRjsVjM0qVLW6x53759RpIJCAgwTz75pMnKyjIVFRXmu+++M5MnTzaSzEMPPXQqp6MRZ5+fxYsXG0mmT58+5r333jO5ubmmuLjYvP/++6Z3794mMDDQfPbZZ43qqKmpMZdeeqnx8fExzz//vCkuLjYZGRlm1qxZxtfX13zxxRcdcn5c4Rx98sknxt/f38TExDR5zIUXXtgoCHvSNWSMMUePHjXJycmmV69e5ptvvjFlZWXm22+/NUlJSSYhIaHJ5WGNMeYPf/iDkWQefvhhU1BQYA4cOGCuueYaY7FYzL/+9a+OOD0den5O1tYgXFBQYIKDg82f//xnh/tKMhaLxTzwwANm7969pqKiwuzYscPMnj3bSDLXX399k6GorZx9flauXGn7hXvRokXm8OHDpqyszCxfvtwMGTLEeHt7m5dffrnJYz3h+vn+++9NZGSkCQoKavJrsuGX0KaCsCdcP8a4VgbyyCBsjDHr1683559/vomIiDCBgYFm9OjR5vnnn2/yIqutrTUXXHCBiYqKavQf9+jRo+af//ynOfvss01cXJzx8fExAQEBZuDAgeb66683P/74Y7M1lJWVmYceesgkJycbf39/07NnT3PNNdeY3bt3t1j7J598YmbOnGlCQ0NNSEiImTx5snn77bdP5TQ0q6POT8NdtJY+WgrCv//9701sbKypqKhwWPPGjRvNvffea0aNGmVCQkKMj4+PiYuLM5dccon56quv2noKWuTM81NVVWU+++wzc/3119uuHX9/f9O/f39z0003mV27djVbd3V1tXnqqafMaaedZgICAkyPHj3MxRdfbH744YeOOC12nHmOLr74YofHnByEjfGca6hBfn6+ufPOO01SUpLx8/MzCQkJ5re//a05cuRIi7W/+uqrZvz48SYoKMiEhYWZs846y3zzzTeneiqa1FHnx5jjf0Vp7qMlTzzxhPH39zdHjx51WPPPP/9sHn74YTN+/HgTHh5uvL29TXR0tDn77LPNW2+91erPvTWceX5qa2vN6tWrzW233WaGDh1qAgMDja+vr0lMTDRXX321+f7771us3d2vn4aw7+jj5PfylOungatkIIsxTljjDwAAAHAyj3tYDgAAAJAIwgAAAPBQBGEAAAB4JIIwAAAAPBJBGAAAAB6JIAwAAACPRBAGAACARyIIAwAAwCMRhAEAAOCRCMIAAADwSARhAPAQM2bMUP/+/Z1dBgC4DIIwAHiAnJwcrV27Vr/85S+dXQoAuAyCMAB4gA8++EA1NTUEYQA4AUEYALrQ0qVLZbFYdP/99zc5ftlll8lisWjjxo0d+r7/+c9/FBsbq8mTJ3fo6wJAd0YQBoAu9OOPP0qSTj/99GbHvb29NXz48A57z6KiIn3zzTe6+OKL5eXVum/7H3/8sWbNmqWIiAgFBQUpJSVFl1xyiV555RWVl5dLklJTU2WxWGwfDz30kL7++muNHz9eQUFB6t+/vx5//HFJUnV1te6//3717t1bgYGBmjlzprZv395hnyMAnAqCMAB0oYYgPHr06EZjBQUF2rt3rwYNGqSAgIAOe89PPvlEFRUVuvTSS1u1/5/+9CdddNFFGjhwoLZv367c3FwtWrRI+/bt03XXXafPP/9ckjRx4kQZY7Ro0SJJ0tq1a/Xqq6/q7bffVmZmpqZNm6Z7771XTz31lO677z5NnjxZu3bt0kcffaQff/xR5513ni1UA4AzEIQBoAtt3rxZYWFhGjBgQKOxhpA8cuTIDn3P//znPwoPD9eZZ57pcN+vvvpKf/7znzVt2jQ999xzSkhIUFBQkKZMmaJ33nmnxWO3bt2qf/3rX0pKSlJMTIwWLFggi8WiRx55RElJSbrooosUGhqqs846S9dcc42ys7P11VdfddSnCQBtRhAGgC6SmZmp/Px8jRo1ShaLpdH4Dz/8IEkaNWpUh71nRUWFPv30U11wwQXy9fV1uP8zzzwjSbrmmmsajQ0cOFD/9V//pb59+zZ57DnnnGP3HjExMYqMjFR+fr4uvPBCu30HDx4sSdq1a1drPxUA6HA+zi4AADxFS9MiThzvyDvCX331lYqLi1s9LSI1NVWSNGTIkCbHH3vssWaP7dWrV6NtoaGhysvLazQWFhYmSSopKWlVXQDQGbgjDABdxFEQ/u677yTZB+FFixbZPZDW8DFo0KBWved7772nwMBAnXvuua3aPz8/X5IUEhLSqv1PFBgY2OYxY0yb3wcAOgp3hAGgi2zbtk2SmuwIsW3bNqWnpys2NlZxcXG27WeffbbWr19v+3dGRoZuuOEGzZkzx+H71dTU6KOPPtLZZ5+t4ODgVtUYGRmpnJwcFRcXt2p/AOjOuCMMAF0kOztbkhQdHd1o7IEHHpDUeFpEfHy8Jk6cqIkTJyoxMVF/+tOfdM011+iRRx5x+H5r1qxRbm5umxbRmDhxoiTpp59+anL8pZde0ooVK1r9egDgygjCANBFYmJiJNX16G1QXl6uO+64Qx999JGk5ucH5+bmatasWRozZoxeeumlJh+2O9l7770nHx8fXXTRRa2u8fe//72kuoU/TvbDDz/oN7/5jYqKilr9egDgygjCANBFGu7M3n777Zo8ebIuuOACJSYm6rPPPrN1aWiqY4TVatW5556rxMREvfbaa/L29m7V+73//vuaMWOGevTo0eoazzrrLD3wwAP69ttvdeutt+rAgQMqLS3V119/rTlz5uiSSy5p1AECALori+FJBQDoEsYY/fWvf9WLL76oQ4cOKSkpSVdeeaXuvfdeXXHFFfrss8+0fft2DRs2zHZMWVmZzj33XFVVVemrr75q9Vzf77//XuPHj9fChQt12223tbnWjz76SP/4xz+0ceNGVVVVqX///po/f75+97vfyd/fX5K0b98+9evXr9Gxe/fu1eLFi/Xwww/bbZ8/f74WL17c5N3sFStWaObMmW2uEwDagyAMAC6qqqpKl1xyibKzs7Vq1SpFRES0+tj//u//1l/+8hdlZ2erd+/enVckAHRjBGEAcFHXX3+93n33XS1ZskSxsbG27TExMU2uTHeiIUOGKDw83NYXGADQGEEYAFyQMUahoaFNLjhx77336q9//asTqgIA90IQBgAAgEeiawQAAAA8EkEYAAAAHokgDAAAAI9EEAYAAIBHIggDAADAIxGEAQAA4JEIwgAAAPBIBGEAAAB4JIIwAAAAPBJBGAAAAB6JIAwAAACP9P8BE75JQGzN69EAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(u*100,z)\n",
    "plt.xlabel(r'$u_z$ / cm')\n",
    "plt.ylabel(r'$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": 26,
   "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": 27,
   "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": 28,
   "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": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "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": 29,
     "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": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAIaCAYAAAApyCWzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZnUlEQVR4nO3deViVdf7/8ddhXxQQUVRAcDc1N9y13Mt9b8/Qlpmm/dc6NTVW00w12VRmy7SM2t43l3Kt3DPX3HcFBQQ3RBRk3+7fHyR6Aj0oB+7DOc/HdXFdeX8Ohzd3EM9uP9zHYhiGIQAAAMAFuJk9AAAAAFBdiF8AAAC4DOIXAAAALoP4BQAAgMsgfgEAAOAyiF8AAAC4DOIXAAAALoP4tcEwDGVkZIjbIQMAANR8xK8N586dU2BgoM6dO2f2KAAAAKgk4hcAAAAug/gFAACAyyB+AQAA4DKIXwAAALgM4hcAAAAug/gFAACAyyB+AQAA4DKIXwAAALgM4hcAAAAug/gFAACAyyB+AQAA4DKIXwAAALgM4hcAAAAuo0bEb25urv72t7/Jy8tLkyZNsvn4JUuWaMCAAQoMDFRAQID69Omj2bNnV/2gAAAAcGgOH78rV65U+/bt9cEHH6igoMDm499++20NGzZMLVq00MGDBxUfH69u3brppptu0osvvlj1AwMAAMBhWQzDMMwe4lK+/vpr/eUvf9G//vUv+fr66u6771ZMTIxmzpxZ7uN3796tjh07qnPnztq4caMsFkvpWv/+/bVq1SqtX79ePXr0qPAMGRkZCgwMVHp6ugICAir7KQEAAMBEDn3lt0mTJtq7d68eeOABq5C9lP/85z8qKirSfffdV+bxf/7znyVJb7zxRpXMCgAAAMfnYfYAl3MlV2glacGCBZKk3r17l1nr06ePJGnx4sUqLCyUh4dDf+oAAACoAg595fdKHD16VKmpqZJKrhj/UVhYmLy8vJSbm6uDBw9W93gVkplXqF9jU/VrbKriUs6ZPQ4AAIDTcZr4PXHihCTJx8dHvr6+ZdYtFouCgoKsHluevLw8ZWRkWL1Vl6S0bN356Ubd+elGffprQrV9XAAAAFfhNPGbnZ0tSfLy8rrkY7y9va0eW55XX31VgYGBpW8RERH2HRQAAACmcZr49fPzkyTl5+df8jF5eXlWjy3Ps88+q/T09NK3pKQk+w4KAAAA0zjNb301aNBAUskLYuTk5JTZ+mAYhs6ePStJatiw4SWfx9vbu/QKMQAAAJyL01z5DQsLU0hIiCQpPj6+zPrRo0eVn58vX19ftWjRorrHAwAAgANwmviVpJEjR0qS1q1bV2Zt7dq1kqShQ4dymzMAAAAX5VTx+8QTT8jd3V2ffPKJ/vjCdR999JEk6amnnjJjNAAAADgAp4rftm3baurUqdq4caMeeOABpaSkKC0tTU8++aRWrFihKVOmXPELZwAAAMB5OHz8WiwWWSwWTZ48WZI0a9as0mMzZ84s8/jHHntMixYt0v79+9W8eXNFRkZq/fr1+u677/Tiiy9W7/AAAABwKA6/+fWP2xcqYtiwYRo2bFgVTAMAAICazOGv/AIAAAD2QvwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfwCAADAZRC/AAAAcBnELwAAAFwG8QsAAACXQfw6qLyCIrNHAAAAcDrErwNpHOwnL4+SfyVr4lJVVGyYPBEAAIBzIX4diL+3h/q1rCdJOnUuT78lpJk8EQAAgHMhfh3M8PYNS/950c7jJk4CAADgfIhfBzPwmlB5/771Ycnu42x9AAAAsCPi18HU8vZQ/1b1JUmpmfnaGH/a5IkAAACcB/HrgNj6AAAAUDWIXwc0oHV9+XiW/Kv5ac8JFRYVmzwRAACAcyB+HZC/t4cGtL6w9WFTPHd9AAAAsAfi10ENv7ZR6T8v3MXWBwAAAHsgfh1U/9b15OvpLkn6cTdbHwAAAOyB+HVQfl4Xtj6kZeVrw2G2PgAAAFQW8evArO76sOuYiZMAAAA4B+LXgfVvVd9q60MBWx8AAAAqhfh1YL5e7hp4TcnWhzPZBVp/iBe8AAAAqAzi18GN4AUvAAAA7Ib4dXD9WtWXn9fvWx/2sPUBAACgMohfB+fj6a5B14RKktJzCrQ2LtXkiQAAAGou4rcGuPiuD4t5wQsAAICrRvzWAH1b1pP/71sfftpzUvmFbH0AAAC4GsRvDeDj6a5BbS7a+nCIrQ8AAABXg/itIYZfy10fAAAAKov4rSGub1lPtbw9JEk/7TnB1gcAAICrQPzWED6e7hr8+9aHc7mF+jXulMkTAQAA1DzEbw1y8daHhWx9AAAAuGLEbw1yXcsQ1f5968PSPSeVV1hk8kQAAAA1C/Fbg3h7uGtw29+3PuQVas1B7voAAABwJYjfGmbERS94sYgXvAAAALgixG8N06d5PdX2Kdn6sGzvSeUWsPUBAACgoojfGsbLw003tGkg6fetD7FsfQAAAKgo4rcGstr6sPOYiZMAAADULE4Zv2vWrNG4cePUtGlT+fr6KiIiQoMHD9aCBQvMHs0uejcPUcDvWx+WsvUBAACgwpwufqdPn67rr79esbGx+uKLL3T69GmtWLFCgYGBGjVqlB555BGzR6w0Lw833di2ZOtDVn6RVh/kBS8AAAAqwmIYhmH2EPZSUFCgkJAQZWRkaNOmTeratWvpWk5OjsLDw5WWlqb9+/erVatWFXrOjIwMBQYGKj09XQEBAVU1+hVbdSBFk2b8Jkka1aGRpt3WyeSJAAAAHJ9TXfk9c+aMMjIyJElt27a1WvP19VWzZs0kSdu3b6/u0eyud/MQBfp6SpKW7WPrAwAAQEU4VfzWr19fERERkqQ9e/ZYreXm5urQoUOSpIYNG5Z535rG091NQ37f+pCdX6RVB1JMnggAAMDxOVX8StIXX3yh0NBQ3X333Vq/fr1ycnIUGxurO++8U2lpaeratat69+5t9ph2Mfyiuz4s3MkLXgAAANjidPF7/fXXa9OmTWrdurV69eolPz8/tWzZUr/88oseeOABLV26VO7u7pd8/7y8PGVkZFi9OaqezeoqyK9k68PyfSnKyWfrAwAAwOU4XfwuWLBAnTp10uHDh/Xrr7/q3Llz2rdvnyZOnKi8vDxlZWVd9v1fffVVBQYGlr6d30bhiC7e+pBTwNYHAAAAW5zqbg+JiYlq06aN3NzcFBcXp9DQ0NI1wzDUtWtXJSQkaOvWrWrcuHG5z5GXl6e8vLzSP2dkZCgiIsLh7vZw3prYU5r46SZJJdsg3ru9s8kTAQAAOC6nuvL7zTffKDs7WwMGDLAKX0myWCy67bbbdPr0ab344ouXfA5vb28FBARYvTmynk3rqs7vWx9W7EtRdn6hyRMBAAA4LqeK3/j4eEmXvptDo0aNJElbtmyptpmqmoe7m4a0K/l8cwqKtHI/L3gBAABwKU4VvyEhIZKkY8eOlbt+/rinp2e1zVQdRlx014dFu8r/3AEAAOBk8Ttq1ChJ0sqVK5WSYv3LX4Zh6JtvvpEkDR48uNpnq0rdmwSrrr+XJGnF/hRl5bH1AQAAoDxOFb/dunXT448/rszMTA0fPlzr169XVlaWDhw4oIkTJ2rz5s1q166dnnnmGbNHtauSrQ8ld33ILSjWiv3c9QEAAKA8ThW/kvTmm29q3rx5CgkJ0ahRoxQYGKguXbpo//79+uc//6kNGzYoKCjI7DHtbvi1F2194AUvAAAAyuVUtzqrChkZGQoMDHTYW52dV1hUrB6vLldqZr68Pdy0+flBqu3jXHubAQAAKsvprvy6qou3PuQVFuvNnw+aPBEAAIDjIX6dyD19msrbo+Rf6cx1Cdpw+LTJEwEAADgW4teJNAnx19NDWpf++anZO7jzAwAAwEWIXyczuVeUukbVkSQlpeXotSX7TZ4IAADAcRC/TsbNzaI3JnSQj2fJv9rPNyRqbVyqyVMBAAA4BuLXCUWF+OuvF21/eHr2TmWy/QEAAID4dVZ39YxS9ybBkqSjZ3P0r8X7TJ4IAADAfMSvkzq//cHPy12S9NXGI/rl4CmTpwIAADAX8evEGtf107NDL2x/+OucncrILTBxIgAAAHMRv07uju6R6tWsriTpWHqu/rmQ7Q8AAMB1Eb9Ozs3NotfHt5f/79sfvt2cpFUHUkyeCgAAwBzErwuICPbT34a3Kf3zX+fsUnoO2x8AAIDrIX5dxG3dInRdixBJ0omMXP1j4V6TJwIAAKh+xK+LsFgsem18e9Xy9pAkzd6SrOX7Tpo8FQAAQPUifl1IWJCvXhhxTemfn527S+nZbH8AAACug/h1MTd3iVDflvUkSSnn8vTSgj0mTwQAAFB9iF8XU7L94VrV9inZ/jB321H9vOeEyVMBAABUD+LXBTUM9NWUkW1L//zcvN06k5Vv4kQAAADVg/h1UeM7h2lA6/qSpNTMPE2Zz/YHAADg/IhfF2WxWPTquGsV8Pv2h/k7junH3cdNngoAAKBqEb8uLDTARy+NvrD94W/zdut0Zp6JEwEAAFQt4tfFjekYpsFtQiVJp7Py9Xe2PwAAACdG/Lo4i8Wif45tpyA/T0nSop3HtWgn2x8AAIBzIn6h+rV99NKoC9sfXvhht1LZ/gAAAJwQ8QtJ0qgOjTSkbQNJUlpWvl74frcMwzB5KgAAAPsifiGpZPvDK2PbKdjfS5K0ZPcJLWD7AwAAcDLEL0qF1PLWP0a3K/3z33/YrZRzuSZOBAAAYF/EL6wMb99Qw9s3lCSdzS7Q3+ax/QEAADgP4hdlvDyqrer+vv1h6d6T+mH7MZMnAgAAsA/iF2XUreWtV8Zc2P4wZf4encxg+wMAAKj5iF+Ua+i1DTWqQyNJUnpOgZ6bu4vtDwAAoMYjfnFJL41qq5Ba3pKk5ftTNGfrUZMnAgAAqBziF5dUx99L/xp7YfvDSwv26EQ62x8AAEDNRfzism5o20BjO4VJks7lFuqvc3ey/QEAANRYxC9smjKyjerXLtn+sOrAKX23OdnkiQAAAK4O8Qubgvy89Oq4a0v//I+Fe3X0bI6JEwEAAFwd4hcVMvCaUE2IDpckncsr1F/nsP0BAADUPMQvKuyFEW3UIMBHkrQmNlXf/JZk8kQAAABXhvhFhQX6eurV8Re2P7yycK+Sz2SbOBEAAMCVIX5xRfq3qq9bukRIkrLyi/T07J0qLmb7AwAAqBmIX1yxv424Ro0CS7Y/rDt0Wl9uOmLyRAAAABVD/OKKBfh46vUJ7Uv//OrifUpKY/sDAABwfMQvrsp1Lerp9u6NJUnZ+UV6avYOtj8AAACHR/ziqj037BqFBflKkjYcTtPnGxJNnggAAODyiF9ctVreHvr3RdsfXluyX4mns0ycCAAA4PKIX1RK7+YhmtgjUpKUU1Ckp77j7g8AAMBxEb+otL8Oba2I4JLtD5sS0vT6T/t59TcAAOCQiF9Umr+3h96Y0KH0z/9dfViv/3iAAAYAAA6H+IVd9GhaV/8Y3bb0zx+uPqTXfuQKMAAAcCzEL+xmYs8o/WNMu9I//3f1Yb26hAAGAACOg/iFXU3sEal/jr0QwB/9clj/XLSPAAYAAA6B+IXd3dE9Uv8ae23pnz/5NV6vEMAAAMABEL+oErd3b6zXxl0I4E9/jdfLC/cSwAAAwFTEL6rMrd0a6/Xx18piKfnzjLUJemkBAQwAAMxD/KJK3dK1sV4f1740gGeuS9CL8/cQwAAAwBTEL6rczV0j9O/xFwJ41vpE/f0HAhgAAFQ/4hfV4qYuEZo6oUNpAH++IVEv/LCbl0IGAADVivhFtRkfHa43b7oQwF9sOEIAAwCAakX8olqN6xyu/9zcQW6/B/CXG4/ob98TwAAAoHoQv6h2YzuF661bOpYG8Nebjui5ebsIYAAAUOWIX5hidMcwqwD+5rckPTuXAAYAAFWL+IVpRncM0zu3dpL77wX87eYkPTNnJwEMAACqDPELU43s0Ejv3NqxNIC/25Ksp+fsVBEBDAAAqgDxC9ONaN9I0y66Ajx7S7Kenk0AAwAA+yN+4RCGt2+o6bd1ksfvATxna7Ke+m4HAQwAAOyK+IXDGHptQ02//UIAz912VE8SwAAAwI6IXziUIe0aavrtnUsDeN62o3r8/7arsKjY5MkAAIAzIH7hcIa0a6D37+gsT/eSAP5h+zE9/n87CGAAAFBpThu/27Zt08SJE9W4cWN5e3urXr166t69u5588kmdO3fO7PFgww1tG+j9O6JLA3j+jmP6fwQwAACoJKeM3w8++EA9e/ZUq1attHbtWmVkZGjp0qVyd3fXm2++qdOnT5s9IipgcJtQfXhntLzcS75MF+w4pke/ZQsEAAC4ehbDMJzqt4lWrVqlAQMGaOrUqXr88cet1g4ePKhRo0Zp+fLlCgsLq9DzZWRkKDAwUOnp6QoICKiKkWHDiv0ndf/nW5X/e/QOv7ah3r61ozzdnfL/3QAAQBVyuvjt0qWL4uLilJKSIi8vr0o/H/HrGFbuT9GfP99SGsBD2zXQtNs6EcAAAOCKOFU57N27V1u2bFGvXr3sEr5wHP1b19dHd0XLy6PkS3bJ7hN6+KttKmALBAAAuAJOFb9r166VJEVGRurnn39W//79FRAQID8/P3Xs2FFvvvmmCgsLTZ4SV6tfq/r6+K4upQH8454TeuirrcovJIABAEDFOFX8xsbGSpJ++ukn3XzzzbrvvvuUmJiouLg49ezZU08++aRGjx6toqKiSz5HXl6eMjIyrN7gOPq2rKdP7uoi798D+Kc9JwlgAABQYU4Vv+np6ZKk+Ph4vfHGG7r99ttVp04dNWrUSB988IE6deqkxYsX69NPP73kc7z66qsKDAwsfYuIiKiu8VFB17esp09iLgTwz3tP6oEvCWAAAGCbU8Xvxe64444yx2JiYiRJX3311SXf79lnn1V6enrpW1JSUpXNiKt3XYt6+jSma2kAL9t3Ug98uUV5hZe+qg8AAOBU8RscHCxJCgkJkZ+fX5n1Jk2aSCq55dmleHt7KyAgwOoNjqlPixD9b1JX+XieD+AUPfDFVgIYAABcklPFb5s2bSRJ+fn5l32cxWKpjnFQDXo3tw7g5ftTdP/nW5RbQAADAICynCp+Bw0aJDc3N2VkZCgtLa3MekJCgqQLkQzn0KtZiGZM6iZfT3dJ0soDp3T/FwQwAAAoy6nit2HDhrr11lslSTNnziyz/tlnn0mS7r333uocC9WgZ7O6mjG5a2kArzpwSn/mCjAAAPgDp4pfSXr77bfVqlUrvfjii/r222+Vk5Oj48eP66GHHtKWLVt033336ZZbbjF7TFSBHk3raubkrvLzKgng1QdP6b7PNhPAAACglNO9vLFUcsuzf/3rX5o7d66OHDkiX19fderUSX/+859LrwxXFC9vXPP8lpCmSf/bpKz8kui9rkWIPr6ri3x+vyoMAABcl1PGrz0RvzXT5oQ0xVwUwH2alwSwrxcBDACAK3O6bQ+AJHWJCtZn93RTLW8PSdKvcam697PflJPPFggAAFwZ8QunFR0ZrFl3XwjgtXGndc8sAhgAAFdG/MKpRUfW0Wf3dFPt3wN43aHTunvmb8rOLzR5MgAAYAbiF06vc2PrAF5/+LQmzyCAAQBwRcQvXEKnxnX0+b3dVdunJIA3xqdp0ozflJVHAAMA4EqIX7iMjhFB+uKeCwG8KT5NkwlgAABcCvELl9IhIkhf3ttdAecDOCFNk2ZsUiYBDACASyB+4XLahwfpq/t6KNDXU5L0W8IZTfofAQwAgCsgfuGS2oUF6st7u5cG8ObEM4r53yalZxeYPBkAAKhK1RK/v/zyS3V8GOCKnA/gIL+SAN6SeEbDpq3RtiNnTJ4MAABUlWp5eWN3d3cVFdXMFxbg5Y2d395jGbrz041Ky8qXJHm4WfTMkNa697omslgsJk8HAADsyaOyT5Cfn6+lS5fqwIEDysrKUjW0NGBXbRoFaOHDffTI19u0OfGMCosN/XPxPm04fFpTb+qgOv5eZo8IAADspFJXfvfs2aNhw4YpOTn5stFrsVi48guHV1BUrP8sPagPVh0qPdYo0EfTbuukLlHBJk4GAADspVLx279/f23btk133HGHWrduXW4cGoahe+65h/hFjbHqQIoe/78dpdsg3N0seuKGlrr/+mZyc2MbBAAANVml4tfPz08LFy7UgAEDLvs4Nzc3FRcXX+2HMRXx65pOpOfqkW+2aVN8Wumxvi3r6T83d1DdWt4mTgYAACqjUnd7qFOnjrp3727zcTU1fOG6GgT66Kt7u+uRAc11/nfeVh88pWHT1mjj4dPmDgcAAK5apeI3JiZGq1atsvm4u+++uzIfBjCFh7ubHr+hlT6/u7tCfr/aezIjT7d9vEHvLo9VUTG/3AkAQE1TqW0PRUVFeuaZZxQVFaURI0aoYcOG8vYu+1fC3OoMNV3KuVw99s12rTt04apv7+Z19dYtHVW/to+JkwEAgCtRqSu/7u7uioiI0IsvvqhmzZrJz89P7u7uZd6Amq5+bR99fk93PT64pc7/ztvauNMa9s6vWhuXau5wAACgwip15Xfq1Kl6+umnZbFYVKdOHdWuXbvMYwzDUFJSEld+4TTWHzqtR7/ZppRzeZIki0V6eEALPTqwhdy5GwQAAA6tUvHbrFkzXXfddZo6dapCQkIu+Tju9gBnk5qZp//37Xatib1w1bd7k2BNu62TQgPYBgEAgKOq1LaHEydO6M0337xs+ErSlClTKvNhAIcTUstbsyZ309NDWpVe7d0Yn6Zh76zR6oOnTJ4OAABcSqXit127dsrOzrb5uMmTJ1fmwwAOyc3Nogf6Ndc3f+qhhoElV3tPZ+Ur5n+b9PqP+1VYVDP/tgMAAGdWqW0Py5cv1+eff65PP/30sr/Yxt0e4OzSsvL15Hc7tGJ/SumxLpF1NO22TmoU5GviZAAA4GIelXnntWvX6tixYwoPD9eNN96ohg0byteXH/RwPcH+Xvrkri769Nf4kqu+xYY2J57RsGlr9J+bO2hA61CzRwQAAKrklV83t4rtmrBYLFz5hcvYeuSMHv5qm46ezSk99qfrm+qpG1vJ071SO40AAEAlVTp+Z8yYcdnHGIahe+65h/iFS0nPLtBTs3fo570nS491jAjS9Ns7KbyOn4mTAQDg2iodvxW5hRm3OoMrMgxDM9cl6F+L96mgqOTbLMDHQ2/c1EE3tm1g8nQAALimSv0d7Pr16yv0uPj4+Mp8GKBGslgsmty7ieb8pZcaB5dc7c3ILdSfP9+ilxbsUX5hzfwfQgAAarJKXfl1BVz5hT1k5Bbor3N2avGuE6XH2ocHavptndW4LtsgAACoLvz2DVANAnw89d7tnfWPMe3k5VHybbczOV3Dp63R4l3HTZ4OAADXwZVfG7jyC3vbcyxdD321TfGpWaXHJvaI1N+GXyMfz0vfLxsAAFQeV36Bata2UaAWPNxHozo0Kj32+YZEjXt/nVUQAwAA+yN+ARPU8vbQO7d21KvjrpX379sg9h7P0Ihpa/TD9qMmTwcAgPNi24MNbHtAVdt/IkMPfrlVh05duOp7W7cITRnZlm0QAADYGVd+AZO1bhCg+Q/10bjOYaXHvt6UpNHT1youJdPEyQAAcD5XFL/NmzfX008/rXXr1lXVPIBL8vf20H9u7qg3JrSX7+9Xew+cPKeR7/6qOVuSTZ4OAADncUXbHvr27at169apuLhYoaGhGjNmjMaNG6f+/fvL3d05/3qWbQ+obrEnz+nBr7bq4MkLV30nRIfr5dFt5eflYeJkAADUfFe85zc1NVXff/+95s2bpxUrVig/P1+BgYEaMWKExo0bpxtvvFG+vr5VNW+1I35hhpz8Ir04f4++3ZxUeqxF/Vp6747Oahla28TJAACo2Sr1C2+ZmZlauHCh5s2bpyVLligzM1N+fn668cYbNXbsWI0YMUJBQUF2HLf6Eb8w0/fbjuq5ebuUnV8kSfLxdNNLo9rq5i4RslgsJk8HAEDNY7e7PeTn5+vnn3/W3LlztXDhQqWmpsrT01N9+/bVuHHjNGbMGDVo0MAeH6paEb8w2+FTmXrwq23adzyj9NiYjo30ythrVcubbRAAAFyJKrnVWXFxsX755RfNnTtXP/zwg5KSkuTm5qbu3btr7dq19v5wVYr4hSPILSjSPxbu1Zcbj5Qeaxrir+m3d1abRnxdAgBQUdVyn9/Nmzdr7ty5+v7777V3796q/nB2RfzCkSzceUx/nbNLmXmFkiQvDzdNGdlGt3drzDYIAAAqgBe5sIH4haNJSM3SQ19v1e6jF7ZBDG/fUK+Nu1a1fTxNnAwAAMfHi1wANUxUiL/m/KWXJvWKKj22aOdxjXj3V+0+mm7eYAAA1ABc+bWBK79wZD/uPq6nZu/Uudzft0G4u+lvw6/RXT0j2QYBAEA5iF8biF84uqS0bD301VbtSL5w1XdI2wZ6fUJ7BfqyDQIAgIux7QGo4SKC/fTd/b10b58mpcd+3HNCw6et0faks+YNBgCAA7rq+F23bp0yMzNtPxBAlfPycNPzI9rok7u6lF7tTT6To5s+XKdP1hwWf8EDAECJq47fPn36qF69elqwYIE95wFQCYPahGrxo9epc+MgSVJBkaFXFu3TfZ9t1tnsfHOHAwDAAVRq20NeXp5uuummywbwsmXLtGTJksp8GABXICzIV9/+uaf+3Ldp6bFl+1I07J012pKYZuJkAACYr1LxO2DAADVs2FA33XSTFi5cWO5jvvzyS40YMaIyHwbAFfJ0d9OzQ6/RjMldFezvJUk6lp6rm/+7QR+uPqTiYrZBAABcU6Xit3Hjxlq1apUaNGigCRMmXDKAAZijf6v6WvzIdeoWFSxJKio29NqS/bp71m86nZln8nQAAFS/St/tITIykgAGHFiDQB99dV93PTyguc7f+nfVgVMaNm2NNh4+be5wAABUM7vc6iwqKsoqgBctWmSPpwVgJx7ubnrihlb67O5uCqlVsg3iZEaebvt4g6aviGUbBADAZdjtPr8XB/D48eMJYMABXdeinhY/cp16Nq0rSSo2pKk/H9Qdn2xUXAq3LgQAOD+7vsjFH68AswUCcDz1A3z0xb3d9f8GtZTb79sg1h8+rSFv/6JXF+9TZl6huQMCAFCFrjp+x48fr8aNG5c5fj6AQ0NDddNNN2nDhg2VGhCA/bm7WfTooBb64t7uCgvylSQVFhv67y+HNfDNVZq/4xgvjAEAcEoWo4p+wiUkJKhfv346cuSILBaLioqKquLDVLmMjAwFBgYqPT1dAQEBZo8D2F1OfpE+WH1IH64+pPzC4tLjPZoG66VR7dSqQW0TpwMAwL6qLH4lKTExUXfddZe2bt2qc+fOVdWHqVLEL1xF4uksvbxgr5bvTyk95u5m0aReUXp0UAsF+HiaOB0AAPZRpfF7nmEYspy/x1INQ/zC1Szfd1IvLdirI2nZpcdCannruWGtNbZTWI39XgYAQKqm+K3JiF+4otyCIn30y2G9tzJOeRdthegSWUcvj26nNo34XgAA1EzErw3EL1xZUlq2Xlm0Vz/tOVl6zM0iTewRqccHt1KgH1shAAA1C/FrA/ELSKsPntJL8/focGpW6bG6/l56ZkhrTYgOl5sbWyEAADUD8WsD8QuUyCss0qe/xuvd5XHKKbhw95aOEUF6eXRbtQ8PMm84AAAqiPi1gfgFrB07m6N/Lt6nRTuPlx6zWKTbujXWUze0Uh1/LxOnAwDg8ohfG4hfoHxr41I1Zf4eq5dFDvLz1FM3ttKtXRvLna0QAAAHRPzaQPwCl1ZQVKyZaxP09rKDysq/sBXi2rBAvTS6rTo3rmPidAAAlHXVL29ckxw5ckQBAQGyWCxatWqV2eMATsPT3U33Xd9UK5/spzEdG5Ue33U0XePeX6enZ+9QamaeiRMCAGDNJeL33nvvrbGvMAfUBPUDfPT2rZ307Z96qPVFL4f8f5uTNWDqKs1al6DCouLLPAMAANXD6eP3o48+0ubNm9W6dWuzRwGcXvemdbXw4T6aMrKNant7SJIycgs1Zf4ejZy+Vr8lpJk8IQDA1Tl1/B45ckRPPfWU3nrrLYWGhpo9DuASPNzdNLl3E614sp8mRIeXHt93PEM3fbhe/+/b7UrJyDVxQgCAK3Pq+L333nvVu3dvxcTEmD0K4HLq1fbW1Js6aM5feqrtRS+HPG/bUQ14c7U+WXNYBWyFAABUM6eN348//lgbN27Uf//7X7NHAVxadGSw5j/UR/8Y006BviUvh5yZV6hXFu3T8GlrtP7QaZMnBAC4EqeM36SkJD355JOaOnWqIiIiruh98/LylJGRYfUGoHLc3Sya2CNSK5/sp9u6Rcjy+y2AD57M1G0fb9DDX2/TiXS2QgAAqp5Txu99992nrl276r777rvi93311VcVGBhY+nal8Qzg0oL9vfTquPb6/oHe6hAeWHp8wY5jGvDmKn24+pDyC9kKAQCoOk4Xv5988ol+/fVXffLJJ1f1/s8++6zS09NL35KSkuw8IYAOEUGa90BvvTbuWtXxK9kKkZ1fpNeW7NeQd37RmthTJk8IAHBWThW/ycnJeuKJJ/Taa68pKirqqp7D29tbAQEBVm8A7M/NzaJbuzXWyif7aWKPSJ1/NeTDp7I08dNN+ssXW3T0bI65QwIAnI5TvbzxzJkzNXny5Ao/fuXKlerXr99lH8PLGwPVY/fRdP39h93aeuRs6TEfTzc91L+57ru+qbw93M0bDgDgNJwqfi+nX79+Wr16dYWC92LEL1B9iosNzd12VK8t2afUzPzS41F1/TRlZFv1b13fxOkAAM7AqbY9AKjZ3NwsmhAdruVP9NPk3lFy/30vRMLpbE2e+ZvunfWbjpzONnlKAEBNRvwCcDiBvp6aMrKtFj7cR92aBJceX7YvRYPeWq23lh5UbkGRiRMCAGoqp972kJCQoCZNmpS71rdvX61atcrmc7DtATCXYRiav+OY/rlon1LO5ZUeD6/jq7+PaKPBbUJlOX/jYAAAbHDq+LUH4hdwDOdyCzRteaxmrE1QYfGF/2z1a1VPU0a2VZMQfxOnAwDUFMSvDcQv4FhiT57TlPl7tO6il0X2cnfTfdc30YP9m8vPy8PE6QAAjo74tYH4BRyPYRhavOuEXlm0V8cvelnkRoE+en5EGw1t14CtEACAchG/NhC/gOPKyivUeyvj9PGawyoouvCfsj7NQ/TiqDZqXr+2idMBABwR8WsD8Qs4vsOnMvXigr365eCFl0X2cLPonj5N9PDAFqrlzVYIAEAJ4tcG4heoGQzD0M97T+rlBXutXhY5NMBbzw27RqM6NGIrBACA+LWF+AVqlpz8In2w+pA+XH1I+YXFpce7NwnWy6PbqVUDtkIAgCsjfm0gfoGaKfF0ll5esFfL96eUHnN3syimZ5QeG9xCAT6eJk4HADAL8WsD8QvUbMv3ndRLC/bqSNqFl0UOqeWtZ4e21rjOYWyFAAAXQ/zaQPwCNV9uQZE++uWw3lsZp7yLtkJ0iayjl0a3VdtGgSZOBwCoTsSvDcQv4DyS0rL1yqK9+mnPydJjbhbpzh6RemJwKwX6sRUCAJwd8WsD8Qs4n9UHT+ml+Xt0ODWr9Fiwv5eeGdJKN0VHyM2NrRAA4KyIXxuIX8A55RUW6X+/JujdFbHKzi8qPd4hIkj/GN1W7cODzBsOAFBliF8biF/AuR1Pz9Eri/Zp0c7jpccsFunWro319I2tVMffy8TpAAD2RvzaQPwCrmFtXKqmzN+juJTM0mNBfp568oZWuq1bY7mzFQIAnALxawPxC7iOgqJizVyboLeXHVTWRVsh2oUF6OXR7dS5cR0TpwMA2APxawPxC7ielIxc/WvxPn2//ZjV8Zuiw/XM0NYKqeVt0mQAgMoifm0gfgHXtfHwaU2Zv0f7T5wrPVbbx0NPDG6pO3tEysPdzcTpAABXg/i1gfgFXFthUbE+35Co//x8UOfyCkuPt25QWy+PbqduTYJNnA4AcKWIXxuIXwCSdOpcnl7/cb9mb0m2Ot63ZT09OqgF+4EBoIYgfm0gfgFcbEtimv7+wx7tOZZhdfz6lvX06MAWio4kggHAkRG/NhC/AP6oqNjQd5uT9O6KOB09m2O1dl2LED06sIW6RLEdAgAcEfFrA/EL4FLyC4s1d2uypq+MU/IZ6wju3byuHh3Ykj3BAOBgiF8biF8AthQUXYjgpDTrCO7VrK4eHdhC3ZvWNWk6AMDFiF8biF8AFVVQVKx5247qvZVxSjydbbXWo2mwHh3YUj2bEcEAYCbi1wbiF8CVKiwq1vfbj+ndFbFlIrhbk2A9NqiFejatK4uFl0wGgOpG/NpA/AK4WoVFxfph+zFNXxmn+NQsq7VuUcF6dFAL9WpGBANAdSJ+bSB+AVRWYVGxFuw8pneXx+nwHyK4S2QdPTqohfo0DyGCAaAaEL82EL8A7KWo2NDCncf0zvJYHT5lHcGdGwfpsUEtdV0LIhgAqhLxawPxC8DezkfwuyviFJeSabXWqXGQHh3YQn1b1iOCAaAKEL82EL8AqkpRsaHFu45r2vJYxf4hgjtGlERwv1ZEMADYE/FrA/ELoKoVFxtasvuE3ll+UAdPWkdwh/BAPTqohfq3qk8EA4AdEL82EL8AqktxsaEf95zQtOWx2n/inNVa+/BAPTKghQZeQwQDQGUQvzYQvwCqW3GxoZ/3ntA7y+O073iG1Vq7sAA9MqCFBrcJJYIB4CoQvzYQvwDMUlxsaOm+k3pnWaz2/iGC2zQM0KODWugGIhgArgjxawPxC8BshmFo6d6Temd5rPYcs47gaxoG6NGBzXVDmwZycyOCAcAW4tcG4heAozAMQ8v3peid5bHadTTdaq11g9p6ZGALDWlLBAPA5RC/NhC/AByNYRhaeSBFby+L1c5k6whuFVoSwUPbEcEAUB7i1wbiF4CjMgxDqw6c0tvLY7Uj6azVWsvQWnp4QAsNu7ah3IlgAChF/NpA/AJwdIZhaPXBU3pneay2HTlrtda8fi09PKC5RrRvRAQDgIhfm4hfADWFYRhaE5uqt5cd1NY/RHCzev56ZGALIhiAyyN+bSB+AdQ0hmHo17hUvbMsVpsTz1itNa3nr4cHNNfI9o3k4e5m0oQAYB7i1wbiF0BNZRiG1h06rXeWxWpTQprVWpMQfz3Uv7lGdySCAbgW4tcG4hdATWcYhtYfPq23l8VqU7x1BEfV9dNDA1poDBEMwEUQvzYQvwCcyfpDp/XO8oPacNg6giPr+unB/s01tlOYPIlgAE6M+LWB+AXgjDYePq13lsdq3aHTVscbB/vpwf7NNK5zOBEMwCkRvzYQvwCc2ab4NL2z/KDWxllHcHgdXz3Uv7nGdQ6XlwcRDMB5EL82EL8AXMHmhDS9szxWa2JTrY6HBfnqwf7NNSGaCAbgHIhfG4hfAK5kS+IZvbM8Vr8cPGV1PCzIV3/p10w3dQmXt4e7SdMBQOURvzYQvwBc0dYjZ/TOslit/kMENwr00V/6N9fNRDCAGor4tYH4BeDKth05o2nLY7XygHUENwz00V/6NdPNXSLk40kEA6g5iF8biF8AkHYkndW05bFavj/F6niDAB/d37epbu3WmAgGUCMQvzYQvwBwwc7kkghets86gkMDvHV/32a6jQgG4OCIXxuIXwAoa/fRdL2zPFZL9560Ol6vdkkE39GdCAbgmIhfG4hfALi0PcfSNW15rH7aYx3BIbW8dX/fprqje6R8vYhgAI6D+LWB+AUA2/Yey9C05bH6cc8Jq+Mhtbz05+ub6Y4ejeXn5WHSdABwAfFrA/ELABW373iG3l0Rq8W7rCO4rr+X/nR9U03sGUkEAzAV8WsD8QsAV+7AiXOatiJWi3cd18U/Zer6e+m+65tqYo9I+XsTwQCqH/FrA/ELAFfv4MlzmrY8Vov+EMHB/l6697omuqtnlGoRwQCqEfFrA/ELAJUXe/Kc3l0RpwU7j1lFcJCfp+67rqnu6hmp2j6e5g0IwGUQvzYQvwBgP3EpmZq+IlbzdxxT8R8i+N4+TRTTK4oIBlCliF8biF8AsL9DpzI1fUWcfth+1CqCA309dU+fJprUO0oBRDCAKkD82kD8AkDVOXwqU9NXxun7bdYRHODjobv7NNHk3k0U6EsEA7Af4tcG4hcAql5Capamr4zTvG1HVXRRBdf28dDdvZvo7j5EMAD7IH5tIH4BoPokns7S9BVxmvvHCPb20OTeUbq7TxMF+XmZOCGAmo74tYH4BYDqd+R0tt5bGac5W5NVeFEE1/L2UEyvSN3VM0qhAT4mTgigpiJ+bSB+AcA8SWnZen9VnL7bbB3BHm4WDWnXQJN6RSk6so4sFouJUwKoSYhfG4hfADBfSQQf0uwtSSoosv6x1bZRgGJ6RWlUh0by8XQ3aUIANQXxawPxCwCO40R6rr7cmKivNh7R6ax8q7U6fp66tVtj3dkjUmFBviZNCMDROV38rlixQl9++aXWrFmjI0eOyM3NTVFRURo+fLiefvpp1atX74qej/gFAMeTW1CkRTuPa9b6BO1MTrdac7NIN7RpoJheUerRNJgtEQCsOFX8vvfee3rooYfUvn17vfXWW+ratasyMjL05Zdf6vnnn1dwcLB++eUXtWzZssLPSfwCgOMyDEPbks5q1roELd51vMyWiFahtRXTK0pjOjWSn5eHSVMCcCROFb9Tp07V3/72N8XFxSkiIsJq7a9//atef/11DRkyREuWLKnwcxK/AFAzpJzL1dcbk/TFxkSdOpdntRbg46FbukZoYo8oNa7rZ9KEAByBU8XvF198oWXLlmnmzJll1tatW6fevXvLw8NDubm5cnev2C9FEL8AULPkFxZrye7jmrUuQVuPnLVas1ikga3rK6ZXlPo0D2FLBOCCnCp+L2fPnj1q166dPD09lZ2dLQ+Piv31F/ELADXXruR0zVyXoAU7jym/sNhqrVk9f8X0itK4zuGq5c2WCMBVuEz8zpkzRxMmTFD//v21YsWKCr8f8QsANd/pzDx981uSvtiQqOPpuVZrtbw9NCE6XHf1jFTTerVMmhBAdXGZ+B09erQWLFigFStWqF+/fpd8XF5envLyLuwVy8jIUEREBPELAE6gsKhYP+89qZnrErQpPq3Met+W9TSpV5T6tqwnNze2RADOyCXi98cff9TQoUP1xBNPaOrUqZd97IsvvqiXXnqpzHHiFwCcy77jGZq1LkHfbz+q3ALrLRFRdf00sWeUbuoSrgAfT5MmBFAVnD5+Y2Nj1adPH1133XX69ttvbf6iG1d+AcC1nM3O17e/JenzDYlKPpNjtebn5a5xncMU0zNKLUJrmzQhAHty6vhNTExU3759FR0dra+//lpeXl5X/Bzs+QUA11BUbGj5vpOatT5Ba+NOl1nv3byuYnpGaeA1oXJnSwRQYzlt/MbFxWngwIG6/vrrNXPmzArf2uyPiF8AcD2xJ89p1voEzd16VNn5RVZr4XV8NbFHpG7pGqEgvyu/qALAXE4Zv3v37tWgQYM0fPhw/fe//5Wbm1vp2tSpU3XLLbeUeRGMSyF+AcB1ZeQW6LvNyfp8fYISTmdbrfl4umlMxzDF9IrSNQ35+QDUFE4Xvzt27NDgwYN1yy23aNq0aWVuYG6xWLRy5crL3vHhYsQvAKC42NDqg6c0c12CVh88VWa9W5NgTeoVpRvahMrD3a2cZwDgKJzqrt6bN2/WDTfcoLy8PJ06dUq33Xab2SMBAJyAm5tF/VvXV//W9XX4VKY+W5+o2VuSlZlXKEnaFJ+mTfFpahjoozt7ROrWrhGqW8vb5KkBlMeprvw+9thjeuedd2w+jiu/AIDKyswr1NytyZq1LkGHTmVZrXl5uGlk+0aa1CtK14YHmjQhgPI4VfxWBeIXAHA5hmHo17hUzVqXoOX7U/THn6qdGwcppleUhrZrKC8PtkQAZiN+bSB+AQAVdeR0tj7fkKBvf0tSRm6h1Vq92t66o3tj3d69serX9jFpQgDErw3ELwDgSmXnF+r7bcc0a12CDpw8Z7Xm6W7RsGsbKqZXlDpFBJX5xWwAVYv4tYH4BQBcLcMwtOFwmmatS9DPe0+o+A8/cduHByqmZ5RGdGgob4+rux89gCtD/NpA/AIA7OHo2Rx9sSFR32w6ojPZBVZrdf29dFu3xrqjR2M1DPQ1aULANRC/NhC/AAB7yi0o0vwdJVsi9hzLsFpzd7NoSNsGiukVpa5RddgSAVQB4tcG4hcAUBUMw9CWxDOauS5BP+4+ocI/7Im4pmGAJvWK1OiOYfLxZEsEYC/Erw3ELwCgqp3MyNWXGxL11aYjSs3Mt1oL8vPULV0jNLFHpMLr+Jk0IeA8iF8biF8AQHXJKyzS4l3HNXNdonYknbVac7NIg64J1aReUerZrC5bIoCrRPzaQPwCAMywPemsZq1L0MKdx1RQZP2jumVoLd3VM0rjOofJz8vDpAmBmon4tYH4BQCY6dS5PH296Yi+3Jiokxl5Vmu1fTx0c5cI3dUzUpF1/U2aEKhZiF8biF8AgCMoKCrWj7tPaNa6BG1OPGO1ZrFI/VvVV0yvKF3XPERubmyJAC6F+LWB+AUAOJrdR9M1a12CfthxTPmFxVZrTUP8dVfPSI2PDldtH0+TJgQcF/FrA/ELAHBUaVn5+ua3I/pifaKOpedarfl7uWtCdLju6hWlZvVqmTQh4HiIXxuIXwCAoyssKtayfSc1c12CNhxOK7N+XYsQTeoVpX6t6sudLRFwccSvDcQvAKAm2X8iQ7PWJWretmTlFlhviWgc7Ke7ekbqpi4RCvRlSwRcE/FrA/ELAKiJ0rML9H+bk/TZhgQlpeVYrfl6umts5zDF9IxSqwa1TZoQMAfxawPxCwCoyYqKDa3cn6JZ6xO0Jja1zHrPpnUV0ytKg66pLw93NxMmBKoX8WsD8QsAcBZxKZn6bH2C5mxJVlZ+kdVaWJCv7uwRqVu7RqiOv5dJEwJVj/i1gfgFADibc7kFmr0lWZ+tT1R8apbVmreHm0Z3bKSYXlFq2yjQpAmBqkP82kD8AgCcVXGxoV9iT2nWugStPHCqzHrXqDqK6RWlG9s2kCdbIuAkiF8biF8AgCtISM3SZ+sT9d3mJJ3LK7RaaxDgozu6N9Zt3RsrpJa3SRMC9kH82kD8AgBcSVZeoeZuO6pZ6xIUl5Jptebl7qYR7RsqpleUOkQEmTMgUEnErw3ELwDAFRmGoXWHTmvmugQt33dSxX+ohY4RQZrUK0rDrm0oLw+2RKDmIH5tIH4BAK4uKS1bX2xI1De/JSk9p8BqLaSWt27v3lh3dG+s0AAfkyYEKo74tYH4BQCgRE5+kX7YflQz1yVo/4lzVmsebhYNvbahJvWKVOfGdWSx8DLKcEzErw3ELwAA1gzD0Kb4NM1an6Cf9pxU0R/2RLQLC1BMzyiN7NBIPp7uJk0JlI/4tYH4BQDg0o6dzdGXGxP19aYkpWXlW60F+3vp1q4RurNHpBoF+Zo0IWCN+LWB+AUAwLbcgiIt3Hlcs9YlaNfRdKs1dzeL+reqrwnRYRrQOpRfkIOpiF8biF8AACrOMAxtPXJWs9YlaPGu4yr8w5aIID9Pje7QSOOjw3VtWCB7g1HtiF8biF8AAK5OSkauvtx4RN/8dkQnM/LKrLcMraXxncM1tlOY6nOnCFQT4tcG4hcAgMopKjb0a1yq5mxJ1k97TiivsNhq3c0iXd+ynsZ3DtfgNqH8khyqFPFrA/ELAID9ZOQWaNHO45qzJVmbE8+UWQ/w8dCIDo00vnO4OjcOYlsE7I74tYH4BQCgasSnZmnu1mTN3XpUR8/mlFlvGuKv8dEl2yK4WwTshfi1gfgFAKBqFRcb2nD4tGZvTdaSXSeUU1BktW6xSL2bhWh8dJiGtG0oXy+2ReDqEb82EL8AAFSfzLxCLdl1XHO2JmvD4bQy67W8PTTs2gYa3zlc3ZoEsy0CV4z4tYH4BQDAHElp2Zq79ajmbE3WkbTsMusRwb4a3zlc4zuHKyLYz4QJURMRvzYQvwAAmMswDP2WcEZztiRr0a7jyswrLPOY7k2CNT46XMOubaha3h4mTImagvi1gfgFAMBx5OQX6ac9JzRna7J+jUvVHyvG19NdQ9s10PjocPVsWldubmyLgDXi1wbiFwAAx3TsbI7mbSvZFnH4VFaZ9bAgX43tFKbx0eFqEuJvwoRwRMSvDcQvAACOzTAMbUs6qzlbkrVgxzFl5JbdFhEdWUfjO4drePuGCvT1NGFKOAri1wbiFwCAmiO3oEjL9p3UnC3JWn3wlIr/UDneHm66oW0Dje8cputa1JM72yJcDvFrA/ELAEDNlJKRq++3H9XsLck6eDKzzHpogLfGdArThM7hahFa24QJYQbi1wbiFwCAms0wDO0+mqE5W5P1w/ajOpNdUOYxHcIDNT46XKM6NFKQn5cJU6K6EL82EL8AADiP/MJirdifojlbk7Vyf4oK/7AvwsvdTQOvqa8J0eG6vmU9ebq7mTQpqgrxawPxCwCAc0rNzNP87cc0e0uy9h7PKLMeUstLozuGaUJ0uK5pSAM4C+LXBuIXAADnt+94huZsSdb3248qNTO/zHqbhgEaHx2u0R0bKaSWtwkTwl6IXxuIXwAAXEdBUbF+OXhKc7Yma9neFOUXFVute7hZ1K9VfU2IDtOA1qHy8mBbRE1D/NpA/AIA4JrOZudrwY5jmr31qHYknS2zXsfPU6M6NNKE6Ai1CwuQxcJt02oC4tcG4hcAAMSlnNPsLUc1b1uyTmbklVlvGVpL4zuHa2ynMNUP8DFhQlQU8WsD8QsAAM4rKjb0a1yq5mxJ1k97Tiiv0HpbhJtFur5lPY3vHK7BbULl4+lu0qS4FOLXBuIXAACUJyO3QIt2HtfsLcnaknimzHqAj4dGdGikCdHh6hQRxLYIB0H82kD8AgAAW+JTszR3a7LmbEnWsfTcMutNQ/w1PrpkW0SjIF8TJsR5xK8NxC8AAKio4mJDGw6f1uwtyVqy+4RyCoqs1i0WqXezEE2IDteNbRvI14ttEdWN+LWB+AUAAFcjM69QS3aVbIvYGJ9WZr2Wt4eGXdtAE6Ij1DWqDtsiqgnxawPxCwAAKispLVtztx7VnK3JOpKWXWa9cbCfxncO17jOYYoI9jNhQtdB/NpA/AIAAHsxDEO/JZzR7C1JWrzrhDLzCss8pnuTYE2IDtewaxvK39vDhCmdG/FrA/ELAACqQk5+kX7ac0KztyRr7aFU/bHIfD3dNbRdA02IDlePpnXl5sa2CHsgfm0gfgEAQFU7djZH87Yd1ZwtyTqcmlVmPSzIV2M7hWl8dLiahPibMKHzIH5tIH4BAEB1MQxD25LOas6WZM3fcUzncstui4iOrKMJ0eEa3r6hAnw8TZiyZiN+bSB+AQCAGXILirRs30nN3pKsXw6eUvEfis3bw003tC3ZFtGneYjc2RZRIcSvDcQvAAAwW0pGrr7fflSztyTr4MnMMuuhAd4a2ylcE6LD1Lx+bRMmrDmIXxuIXwAA4CgMw9DuoxmavSVJP+w4prPZBWUe0yE8UBOiwzWyQyMF+XmZMKVjI35tIH4BAIAjyi8s1or9KZq9JVmrDqSo8A/7Irzc3TSoTX2N7xyu61vWk6e7m0mTOhbi1wbiFwAAOLrUzDz9sP2Y5mxJ1t7jGWXWQ2p5aUzHkrtFXNPQtXuG+LWB+AUAADXJ3mMZmrM1WT9sP6rUzPwy620bBWh853CN7thIdWt5mzChuYhfG4hfAABQExUUFWv1gVOaszVZy/elKL+o2Grdw82i/q1LtkUMaF1fXh6usS2C+LWB+AUAADXdmax8LdhZsi1iR3J6mfU6fp4a3TFM4zuHq11YgCwW571tGvFrA/ELAACcSezJc5q9NVnzth5Vyrm8MuutQmtrfHSYxnQMU/0AHxMmrFpOGb+5ubn697//rS+//FKJiYmqU6eOBg0apJdffllNmjS5oucifgEAgDMqLCrWr3GpmrP1qH7ec0J5hdbbItws0vUt62lCdLgGXRMqH093kya1L6eL35ycHA0aNEg7d+7U//73Pw0fPlyxsbGKiYlRQkKCVq9erQ4dOlT4+YhfAADg7NJzCrRo53HN2ZqsLYlnyqwH+HhoZIdGGh8drk4RQTV6W4TTxe8zzzyjf//735o+fboefPDB0uOHDh1Sq1at1Lp1a+3YsUPu7hX7vxfiFwAAuJLDpzI1d+tRzd2arGPpuWXWm9bz1/jO4RrXOUwNA31NmLBynCp+s7OzVb9+fRUUFOj06dOqVauW1fqQIUP0008/af78+Ro5cmSFnpP4BQAArqi42ND6w6c1Z0uyluw+oZyCIqt1i0Xq0zxE4zuH68a2DeTrVTO2RTjVPS1WrFihrKwsdejQoUz4SlKfPn0kSfPnz6/u0QAAAGoUNzeLejcP0X9u6ajfnh+kf09or+5NgkvXDUNaE5uqx77drq7/XKZnZu/Upvg0Ofp1VQ+zB7CnHTt2SNIlf6mtadOmVo8DAACAbbW8PXRzlwjd3CVCSWnZmrM1WXO2JispLUeSlJlXqG83J+nbzUmKrOuncZ1KtkVEBPuZPHlZTnXl98SJE5Kk4ODgctfr1Klj9bjy5OXlKSMjw+oNAAAAJSKC/fTYoJZa/WR/ffunHrq5S7j8L9rykHg6W28tO6jr/r1SM9fGmzhp+ZwqfrOzsyVJXl5e5a57e3tbPa48r776qgIDA0vfIiIi7D8oAABADefmZlH3pnX17wkd9Nvzg/TWLR3Up3mILr4RRNcm5V+QNJNTbXvw8yu5tJ6fX/Z1rKWSq7oXP648zz77rB5//PHSP2dkZBDAAAAAl+Hn5aGxncI1tlO4jp3N0bxtR7UrOV1tGwWaPVoZThW/DRo0kCSlpaWVu37mTMl96xo2bHjJ5/D29i69QgwAAIAr0yjIVw/2b272GJfkVNseOnbsKEmKjy9/f8nhw4clSe3bt6+ukQAAAOBAnCp++/fvL39/f+3cuVNZWVll1teuXStJGjVqVHWPBgAAAAfgVPHr5+enhx56SHl5eZo1a5bV2uHDh7Vs2TK1bdtWw4YNM2lCAAAAmMmpXuFNknJycjRw4EDt3r1bM2bM0LBhwxQbG6uYmBgdPnxYv/zyizp06FDh5+MV3gAAAJyHU135lSRfX1+tWLFCTzzxhJ599lnVqVNHgwcPVtu2bbVt27YrCl8AAAA4F6e78mtvXPkFAABwHk535RcAAAC4FOIXAAAALoP4BQAAgMsgfgEAAOAyiF8AAAC4DOIXAAAALoP4BQAAgMsgfgEAAOAyiF8AAAC4DOIXAAAALoP4BQAAgMvwMHsAR2cYhiQpIyPD5EkAAABwObVr15bFYrnsY4hfG86dOydJioiIMHkSAAAAXE56eroCAgIu+xiLcf7SJspVXFysY8eOVej/JOwpIyNDERERSkpKsvkvEeXjHFYO569yOH+Vw/mrHM5f5XEOK8es88eVXztwc3NTeHi4aR8/ICCAb7pK4hxWDuevcjh/lcP5qxzOX+VxDivHEc8fv/AGAAAAl0H8AgAAwGUQvw7K29tbU6ZMkbe3t9mj1Ficw8rh/FUO569yOH+Vw/mrPM5h5Tjy+eMX3gAAAOAyuPILAAAAl0H8AgAAwGUQvwAAAHAZxC8AAABcBvHrIDZs2KARI0aobt268vf3V3R0tD766CNd6e8jTpo0SRaLxebbzJkzq+YTMYm9zt95xcXF+vTTT9WvXz/VrVtXPj4+aty4sUaMGKHPP//cztM7BnueQ1tffwkJCfb/BExm76/Bi7311luyWCyKioqq/KAOyl7nLycnR3PmzNHkyZPVvn17BQUFKSAgQNdcc40eeeQRHT58uIo+A3PZ6/xlZmbqww8/1JAhQxQaGipPT08FBwerf//+mjVrll2+nh1RVXz/HjhwQH379pXFYtGqVavsN6wJcnNz9fLLL6tVq1by8fFRw4YNNXHiRMXHx1/xcxUXF2vatGnq0KGD/Pz8FBISorFjx2r79u32H/xSDJhuzpw5hru7uzF69GgjISHByMjIMF5//XXDYrEYkyZNuqLniomJMQICAoxWrVqV+1a/fn1DkrFhw4Yq+myqnz3Pn2EYRlZWljFw4ECjWbNmxg8//GBkZGQY6enpxldffWX4+fkZffv2tf8nYTJ7n0NJRtOmTS/5dZicnFwFn4V57H3+Lnbw4EHD19fXkGRERkbaZ2AHY8/zt2TJEkOS0alTJ+OXX34xMjMzjZMnTxrvv/++4enpadSqVcvYvn17FX0m5rDn+evSpYshybj//vuNgwcPGllZWcb27duNIUOGGJKM8ePHG8XFxVX0mZjD3t+/eXl5xssvv2z4+/sbtWvXNiQZK1eutP/g1SQ7O9vo1auXUatWLeP//u//Sr8mOnToYAQGBl7R91NRUZExfvx4w8PDw3j//feNzMxMIzY21hg4cKDh5eVl/PTTT1X4mVxA/Jrs5MmTRu3atY2IiAgjOzvbai0mJsaQZHzzzTcVfr6YmBgjJibmkuujR482unTpcrXjOhx7nz/DMIxJkyYZ/v7+RmJiYpm11157zZg4cWKlZnY0VXEOJRnx8fF2nNJxVcX5O6+oqMjo06eP0bNnT6eNX3ufvyVLlhgWi8VISEgos/bQQw8Zkoy77rqr0nM7Cnufv7Zt2xpDhgwpczwnJ8eIioqq1NezI6qK798bb7zRiI6ONrZu3Wr07du3xsfv008/bUgypk+fbnU8Li7OcHd3N9q2bWsUFhZW6Lnef/99Q5Lx5JNPWh1PS0sz6tSpY9SvX9/IyMiw2+yXQvya7KWXXjIkGc8991yZtfXr1xuSjM6dO1f4+d5//33j/fffL3ctPj7ecHNzM2bNmnXV8zoae5+/nTt3ll71cBX2PoeG4VrxWxXn77y33nrLqFevnrFp0yanjV97n7/U1FRj9erV5a69++67hiRj1KhRVz2vo7H3+bvllluM7777rty1++67z5BkTJ48+arndTRV8f07b9680his6fGblZVl+Pv7G15eXsa5c+fKrN94442GJGP+/PkVer6mTZsakoyDBw+WWfvzn/9sSDKmTZtW6bltIX5NFh0dbUgyFi1aVGatoKDA8PHxMSTZ5a+Jn3zySaN+/fpGbm5upZ/LUdj7/D311FOGJOPbb7+196gOqyq+Bl0pfqvqezg2Ntbw9fU1vvnmGyM+Pt5p47c6/xt42223lXsFqyar7p8hzha/VX3+anr8LliwwJBkdO3atdz1f/zjH4Yk495777X5XLt27TIkGfXq1St3/fPPPzckGYMGDarUzBXBL7yZqKioSHv27JEkNWnSpMy6h4eHIiIiJEk7duyo1MfKzs7Wp59+qj/96U8O+VKDV6Mqzt/atWslSaGhoXrllVfUunVreXt7q27duho+fLjWrFljp+kdQ1V+DX700Ufq3LmzAgMD5e/vr2uvvVYvvPCCzp49W+m5HUVVnT/DMHT33XdryJAhuuWWW+wzrAOqjv8G5uTkaN++fXr88cf17bff6tFHH9Vf/vKXqx/agVTnzxBJ2r9/vySpX79+lX4uR1Dd568mOv95l3d+JKlp06ZWj6uu56os4tdEZ8+eVW5uriQpODi43MfUqVNHknTixIlKfawvvvhC586d0/3331+p53EkVXH+YmNjJUl33323lixZoq+++kpnz57V4sWLFR8fr/79++vLL7+0w/SOoSq/Bn/99Ve99957OnnypBISEnTffffp3//+tzp37qzk5OTKDe4gqur8TZs2TXv27NH7779f+SEdWFX/N/CLL76Qn5+f2rRpo2+//VYzZszQ1KlT5ebmHD/6qvNnyMmTJ/XTTz+pZcuWuu222yr1XI6iOs9fTXX+87bH+anoc506dUrFxcVXPOuV8KjSZ3cBY8eO1b59+67ofT777DN169ZN2dnZpce8vLzKfez5q7QXP/ZqTJ8+XePGjVNYWFilnsfeHO38paenS5JSUlK0efPm0m/G7t276+uvv1bHjh11//33a+jQoZf8Bq5ujnYOJWnhwoUaNGhQ6fv6+PjokUce0dmzZzVlyhTdfffd+vnnn69o5qriaOfv0KFDeu655/Thhx+qQYMGVzSXGRzt/F3szjvv1O23366EhAQtWLBADz/8sN599119/fXXat68+RU/X1Vw5PN3sUcffVQWi0VfffWVPD09K/Vc9lRTzl9Ndf7ztsf5qehznX9srVq1rmjWK0H8VlJ8fLwOHDhwRe9z/gvAz8+v9Fh+fn65j83Lyyvz2Cu1cuVK7dq1yyGvIjnq+Rs2bFhp+J7XoUMHdezYUdu3b9f8+fM1adKkK3rOquKI53D48OHlHn/wwQc1ZcoULV26VMnJyQoPD6/wc1YVRzp/57c79O/fXxMnTryimcziSOevPG5ubmratKkeffRRhYWF6aabbtLQoUO1d+9eh4g4Rz9/kvT2229r9uzZ+vbbbxUdHX3Vz1MVasL5q8nOf972OD8Vfa6KPl9lEL+VVJmbMgcFBcnHx0e5ublKS0tTaGhomcecOXNGktSwYcOr/jjvvvuuOnbsqD59+lz1c1QVRzt/wcHBOnHihBo3blzuepMmTbR9+3YdPHjwque2N0c7h5dTt25dBQUF6ezZs9q3b59DxK8jnb/p06drx44dpfsQawJHOn+2TJgwQSEhIYqLi9PSpUs1bNiwSj9nZTn6+ZsxY4aeeeYZff755xo/fvxVz1pVHP381XTn//YpLS2t3PUrOT8Vfa769etX+dYk59j4VEO5u7urXbt2klTuq6QUFhYqKSlJktS+ffur+hiJiYmaP3++Hn744asf1EFVxflr06aNpEv/n+l5FovlSkZ1WNXxNfhH5/dyOcM5tPf5mzNnjtLT0xUeHm71injnf0EkMTGx9JgzvNqbGV9/58+lM7zKYFWfvw8//FAPPPCAvvvuO6fZ53sxM77+apqOHTtKKv/8SCp9xcSKnB97PldlEb8mGzlypCRp3bp1Zda2bNmi3Nxcde7c+ar36r7//vsKCgrS7bffXqk5HZW9z9+NN94o6dLfnOd/YJ6PZGdg73P46KOP6q677ip3LSUlRRkZGZKk1q1bX+XEjsWe52/VqlUySm5BafV2/usxMjKy9JgzxJtk/6+/G264QY899tgl148dOyap5G8hnEFV/Qx566239MQTT+j777/XqFGjSo/v2bNHH3/8ceWGdiBV/TO4puvfv7/8/f21c+dOZWVllVk/f4eki79GLqVt27Zq2rSpTp06pbi4uEo9V6VV+c3UcFnnX12mcePGRk5OjtXa5MmTL/nqMk8//bTRu3dvY/fu3Zd87uzsbCM4ONh45pln7D63o7D3+UtJSTECAwMNLy8v49ixY1ZrO3bsKL1HYXk3+66p7H0OY2JijPr16xvp6ell3ueFF14wJBkDBw607ydhoqr8Hj7Pme/za+/z17dvX6Nly5bl3s/8559/NiQZfn5+xqlTp+z7iZikKr7+XnnlFaNWrVrl3pt2xowZTvV1WNXfvzX9Pr+GYRjPPPOMIcl47733rI4fOnTI8PDwKPMKbzt37jR69eplPPvss2We64MPPjAkGU899ZTV8TNnzhjBwcG8wpsrmT17tuHu7m6MHTvWSExMNDIyMow33njDsFgs5b5U8alTpwxJhiTjoYceuuTzfvzxx4a7u3u5L9PrTOx9/s6/zvv5/7Dl5eUZGzduNNq2bWv4+fkZS5curYbPqnrZ8xxOmjTJkGRcf/31xvr1643MzEzj1KlTxltvvWV4enoakZGRTvc1WVXfw+c5c/wahn3PX79+/QxJxuDBg41169YZGRkZxvHjx40ZM2YYISEhhru7u1O9yqVh2Pf8Pf/884Yko23btsYtt9xS5q1bt25O93VYld+/zhC/2dnZRs+ePY3atWsbs2fPNrKzs40dO3YYHTt2NAICAozt27dbPf7BBx8sPT+pqalWa0VFRca4ceMMDw8P48MPPzQyMzONuLg4Y/DgwYanp6fx008/VcvnRPw6iPXr1xvDhg0zgoKCDF9fX6NTp07Ghx9+aBQXF5d5bHFxsTF8+HCjbt26l/2Gat++vTF27NgqnNpx2Pv8/fbbb8aYMWOMkJAQw9PT0wgLCzNiYmKM/fv3V/FnYh57ncP09HRj1qxZxqhRo4ywsDDD09PT8Pf3Nzp06GC88MILxpkzZ6rnE6pmVfE9PGXKlNIfIn98mzFjRtV9Miaw1/lLSUkx3n//feOGG24wGjRoYHh4eBg+Pj5G8+bNjcmTJxvbtm2rnk+omtnr/AUGBl7ya+78m7PFr2HY9/t3xowZlz1/NVFOTo7x4osvGi1atDC8vb2N+vXrG3fccYdx6NChMo9dunSpUadOHWPMmDHlnr/CwkLj7bffNq699lrDx8fHCA4ONkaPHm1s3bq1Oj4VwzAMw2IYhmG/TRQAAACA4+IX3gAAAOAyiF8AAAC4DOIXAAAALoP4BQAAgMsgfgEAAOAyiF8AAAC4DOIXAAAALoP4BQAAgMsgfgEAAOAyiF8AAAC4DOIXAFxE37591bRpU7PHAABTEb8A4AJOnTqltWvXauzYsWaPAgCmIn4BwAX88MMPKioqIn4BuDziFwCq0RdffCGLxaJnn3223PXx48fLYrFo8+bNdv248+bNU2hoqHr16mXX5wWAmob4BYBqtG3bNklS586dL7nu7u6udu3a2e1jnjt3TsuXL9fo0aPl5lax/+wvXLhQgwcPVlBQkPz8/NSyZUuNGTNGs2bNUm5uriRpw4YNslgspW8vvviili1bpm7dusnPz09NmzbVG2+8IUkqLCzUs88+q0aNGsnX11f9+vXT7t277fY5AkBFEb8AUI3Ox2+nTp3KrJ09e1bx8fFq1aqVfHx87PYxFy1apLy8PI0bN65Cj3/hhRc0cuRINW/eXLt371ZqaqpmzJihhIQETZo0ST/++KMkqUePHjIMQzNmzJAkrV27Vp999pm+++47JSYm6rrrrtPTTz+tt99+W88884x69eqlAwcOaMGCBdq2bZuGDh1aGtIAUF2IXwCoRtu3b1dAQICaNWtWZu18GHfo0MGuH3PevHkKDAzUgAEDbD526dKleuWVV3Tdddfpgw8+UHh4uPz8/NS7d2/Nnj37su+7c+dOffrpp4qMjFS9evU0depUWSwWvfzyy4qMjNTIkSNVu3ZtDRo0SHfccYeSk5O1dOlSe32aAFAhxC8AVJPExESdOXNGHTt2lMViKbO+detWSVLHjh3t9jHz8vK0ePFiDR8+XJ6enjYf/+6770qS7rjjjjJrzZs311//+ldFRUWV+7433nij1ceoV6+e6tSpozNnzmjEiBFWj23durUk6cCBAxX9VADALjzMHgAAXMXltjxcvG7PK79Lly5VZmZmhbc8bNiwQZJ0zTXXlLv+6quvXvJ9GzZsWOZY7dq1lZaWVmYtICBAkpSVlVWhuQDAXrjyCwDVxFb8btq0SZJ1/M6YMcPql8rOv7Vq1apCH3Pu3Lny9fXVkCFDKvT4M2fOSJJq1apVocdfzNfX94rXDMO44o8DAJXBlV8AqCa7du2SpHLv5LBr1y7FxsYqNDRUDRo0KD1+ww03aP369aV/jouL0z333KObb77Z5scrKirSggULdMMNN8jf379CM9apU0enTp1SZmZmhR4PADUNV34BoJokJydLkkJCQsqsPf/885LKbnkICwtTjx491KNHD0VEROiFF17QHXfcoZdfftnmx1uzZo1SU1Ov6IUtevToIUnat29fuesff/yxVq5cWeHnAwBHQ/wCQDWpV6+epJJ76J6Xm5urxx57TAsWLJB06f2+qampGjx4sKKjo/Xxxx+X+wtzfzR37lx5eHho5MiRFZ7xkUcekVTyYhx/tHXrVv3pT3/SuXPnKvx8AOBoiF8AqCbnr8A+9NBD6tWrl4YPH66IiAgtWbKk9O4K5d3pISMjQ0OGDFFERIS++uorubu7V+jjff/99+rbt6+Cg4MrPOOgQYP0/PPP69dff9UDDzygo0ePKjs7W8uWLdPNN9+sMWPGlLlzAwDUKAYAoFoUFxcbr776qtGkSRPDx8fHaNWqlfH3v//dyMzMNIYOHWpIMnbv3m31PtnZ2cb1119v9OzZ08jMzKzwx9q0aZMhyZg+ffpVzTp//nxj4MCBRmBgoOHn52e0a9fOeOONN4zc3NzSx8THxxuSyrzFx8cbU6ZMKXM8JibGMAyj3PdZuXLlVc0JAFfKYhj8qi0AOKKCggKNGTNGycnJWr16tYKCgir8vs8995xee+01JScnq1GjRlU3JADUMMQvADioyZMna86cOfr8888VGhpaerxevXrlvkLcxa655hoFBgaW3rcXAFCC+AUAB2QYhmrXrl3ui0A8/fTTev31102YCgBqPuIXAAAALoO7PQAAAMBlEL8AAABwGcQvAAAAXAbxCwAAAJdB/AIAAMBlEL8AAABwGcQvAAAAXAbxCwAAAJdB/AIAAMBlEL8AAABwGcQvAAAAXMb/B7IQJhrz+EN4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "u = np.linalg.solve(A,RHS)\n",
    "plt.plot(u*100,z)\n",
    "plt.xlabel(r'$u_z$ / cm')\n",
    "plt.ylabel(r'$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": 31,
   "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": 32,
   "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": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c535ec4fa69848dba2202d4890fe93d8",
       "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=r'nodes'))\n",
    "def plot(num_nodes=10):\n",
    "    fig,ax = plt.subplots()\n",
    "    ax.set_ylabel(r'$z$ / m')\n",
    "    ax.set_xlabel(r'$u_z$ / cm')\n",
    "    solution, z = run_sim(num_nodes)\n",
    "    ax.plot(reference_solution*1e2, reference_z, '--', color='k', label=r'reference solution')\n",
    "    ax.plot(solution*1e2,z,label=r'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": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwcAAAI/CAYAAAAiBqU7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABx2klEQVR4nO3de3zO9f/H8ee1i5lhG7FY5FRNG+aQQwck0jmxrzk0h/wSIjoplRxLSpLDMpWchkyodZC+osM3Jec5hFAIWWt2aGO2uX5/XHZlDNu167o+13Xtcb/drtt8Dtfn81yHj+t1vU8mi8ViEQAAAIBSz8foAAAAAADcA8UBAAAAAEkUBwAAAADOoTgAAAAAIIniAAAAAMA5FAcAAAAAJFEcAAAAADiH4gAAAACAJIoDOIDFYlF6erpYTw8AAMCzURygxDIyMhQYGKiMjAyjowAAAKAEKA4AAAAASKI4AAAAAHAOxQEAAAAASRQHAAAAAM6hOAAAAAAgieIAAAAAwDkUBwAAAAAkURwAAAAAOIfiAAAAAIAkigMAAAAA51AcAAAAAJBEcQAAAADgHIoDAAAAAJIoDgAAAACcQ3EAAAAAQBLFAQAAAIBzKA7gceLjpRo1pGXLjE4CAADgXUwWi8VidAh4tvT0dAUGBiotLU0BAQFOvVdSkhQaKqWmSkFB0t69UnCwU28JAABQatByAI9hsUiDBkkZGdbtjAxp8GBjMwEAAHgTigN4jPh4aeVKKS/Pup2XJ61YYd0PAACAkqNbEUrMFd2K8rsTpaVZWxDymUxSYCDdiwAAAByBlgNIkvLy8jR58mSVL19e33zzjdFxCji/O9GFpazFQvciAAAAR6E4gP7880+1b99eBw4c0OnTp42Oc5Fduwp2J7pQfveiXbtcmwsAAMDbUBxA//zzj2bOnKmRI0caHaVQ4eFSly6S2Vz4cbNZ6trVeh4AAADsR3FgkNOnT+ull16Sr6+v+vXrV6Tzx48fr9DQUPn5+alGjRrq3bu3fvvttxJnue6669S4ceMSX8dZTCYpNlaqVMn65wuPBQRIs2YZkw0AAMCbUBwYYN26dWrcuLFmzZqlnJycK55/6tQpdejQQZMnT9Yrr7yilJQUffnll9qxY4eaNm2q7du3uyC1sYKDrQVCYWMOYmMZjAwAAOAIZYwOUNosWbJEgwcP1sSJE1W+fHn179//iu8ZO3as1q9fr5kzZ6pbt26SpIiICC1fvlyhoaF6+OGHtX37dpnP9btJTU3Vn3/+edlr+vv769prry35L+RCUVHS0qVSQoJ1nIHZLHXubN0PAACAkqM4cLG6detq9+7dCgkJ0bx58654flZWlmJiYuTr66u+ffsWOFa/fn117NhRq1ev1hdffKEHHnhAkvThhx9q8BWm72nXrp3bzUp0Jfndi9ats66QTHciAAAAx6JbkYu1bt1aISEhRT5/7dq1yszMVEREhCpWrHjR8dtuu02SlJCQYNs3aNAgWSyWy748rTDIFxwszZ4tVa9u/Ul3IgAAAMeh5cDN5Y8nqFu3bqHH69WrV+C80iAqiq5EAAAAzkBx4Obyxw5UqVKl0OOVK1cucJ4rZGdnKzs727adnp7usnsDAADAeehW5OaysrIkSb6+voUeL1euXIHz7HXPPfeoR48ekqQnn3xSt99+u/76669Cz33ttdcUGBhoe9WqVatE9wYAAIB7oOXAzfn7+0uSzpw5U+jx/G/w88+z16pVq4p87gsvvKCnn37atp2enk6BAAAA4AVoOXBz1atXlySlpKQUevzkyZOSpBo1argsU7ly5RQQEFDg5S3i46UaNaRly4xOAgAA4HoUB26uSZMmknTJlZAPHjwoSW69wrGnSEqSBg6U/vxTeuwx6zYAAEBpQnHg5tq3b68KFSooMTFRmZmZFx3/4YcfJEkPPvigq6N5FYtFGjRIysiwbmdkSFdYKgIAAMDrUBy4OX9/fw0dOlTZ2dmaP39+gWMHDx7UmjVrFB4ernvvvdeghN4hPl5audK68rJk/blihXU/AABAaUFx4AHGjBmjm2++WSNHjtTy5ct16tQpJSYmKjIyUv7+/lq0aJHMZrPRMT1WUpK11cBkKrjfZLJ2M6J7EQAAKC1MFovFYnSI0sZ04afQ88ydO1f9+vW7aP/p06f1+uuva9GiRTp8+LACAwN15513avz48baF0IySnp6uwMBApaWledzgZItFioyUEhL+bTU4n9ksde4sLV/u+mwAAACuRnGAEvPk4mDnTqlRo6KdFx7u/DwAAABGolsR7BYTE6OwsDC1aNHC6Ch2Cw+XunSxthAUxmyWunalMAAAAKUDLQcoMU9uOZCsYwpCQ6W0NGs3o3wmkxQUJO3ZIwUHGxYPAADAZWg5QKkXHCzFxhYsDCTrdmwshQEAACg9KA4ASVFRBbsX5XcniooyNhcAAIArURwAsnYhio2VKlWybgcESLNmGZsJAADA1SgOgHOCg6XZs6Xq1a0/6U4EAABKGwYko8Q8fUAyAAAArGg5AAAAACCJ4gAAAADAORQHAAAAACRRHKAEvGGFZAAAAPyLAckoMQYkAwAAeAdaDgAAAABIojgAAAAAcA7FAQAAAABJFAcAAAAAzqE4AAAAACCJ4gAAAADAORQHAAAAACRRHABuKz5eqlFDWrbM6CQAAKC0YBE0lBiLoDleUpIUGiqlpkpBQdLevVJwsNGpAACAt6PlAHaLiYlRWFiYWrRoYXQUr2KxSIMGSRkZ1u2MDGnwYGMzAQCA0oGWA5QYLQeOtXSp1KNH4fujolyfBwAAlB4UBygxigPHye9OlJZmbUHIZzJJgYF0LwIAAM5FtyLATZzfnejCkt1ioXsRAABwPooDwE3s2iWtXCnl5RV+PC9PWrHCeh4AAIAzUBwAbiI8XOrSRTKbCz9uNktdu1rPAwAAcAaKA8BNmExSbKxUqZL1zxceCwiQZs0yJhsAACgdKA4ANxIcbC0QChtzEBvLYGQAAOBcFAeAm4mKKti9KL87EdOYAgAAZ6M4ANzM+d2LJLoTAQAA16E4ANxQcLA0e7ZUvbr1J92JAACAK7AIGkqMRdAAAAC8Ay0HAAAAACRRHKAEYmJiFBYWphYtWhgdBQAAAA5AtyKUGN2KAAAAvAMtBwAAAAAkURwAAAAAOIfiAAAAAIAkigMAkuLjpRo1pGXLjE4CAACMxIBklBgDkj1bUpIUGiqlpkpBQdLevSy6BgBAaUXLAVCKWSzSoEFSRoZ1OyNDGjzY2EwAAMA4FAdAKRYfL61cKeXlWbfz8qQVK6z7AQBA6UO3IpQY3Yo8U353orQ0awtCPpNJCgykexEAAKVRGaMDFEeDBg3UrFmzAq+goCCjYwEe5/zuRBd+PWCx/Nu9aPlyY/IBAABjeFTLgY+Pj0wmU4F9tWvXthUKzZs3V7NmzVStWjWDEpZOtBx4np07pUaNinZeeLjz8wAAAPfgUcXBm2++qS1btmjLli3av3+/zp49azt2ftEQEhJyUQvDNddcY0TkUoHiwPNYLFJkpJSQ8O94g/OZzVLnzrQcAABQ2nhUcXC+zMxMbdu2TVu2bNHWrVs1b968i845v2CoVq2amjVrpi+++MKFKUsHigPPdLkxB0FB0p49jDkAAKC08dji4EI+Pj6qUqWKPvvsM1vrwpYtW7Rr1y7l5ORIshYLeYV9TYoSoTjwXEuXSj16FL4/Ksr1eQAAgLE8akDylfj4+Kh169Zq3bq1bV9OTo527NhhKxYA/CsqyloI5Hcvyu9ORGEAAEDp5FUtB1WrVlVSUpLRUUqNmJgYxcTEKC8vT/v27aPlwEOdv0Jy5cp0JwIAoDSjOECJ0a3I88XHS8OHS9OnS926GZ0GAAAYheIAJUZxAAAA4B18jA4AAAAAwD14VHHwyy+/yEsaOgAAAAC341GzFTVs2FD+/v5q3LixmjVrpqZNm6pZs2Zq2LCh0dEAAAAAj+dRYw58fP5t6Dh/gbOyZcvqzJkz8vf315tvvqmmTZsqIiJCfn5+RsQsdRhzAAAA4B08qjjYv3+/Nm/ebFuzYOvWrUpJSSlwTn7RYDabFRoaqmbNmtleTZs2VcWKFY2I7tUoDgAAALyDRxUHhTl06FCBFZG3bNmiEydO2I6f38Lg4+NjWy0ZjkNxAAAA4B08vjgozLFjxy4qGP744w+ZTCbl5eUZHc/rUBwAAAB4B48akFxUISEhCgkJ0f3332/bl5ycrC1bthiYCgAAAHBvXtlyANei5QAAAMA7eNQ6BwAAAACch+IAgOHi46UaNaRly4xOAgBA6Ua3IpQY3YpQEklJUmiolJoqBQVJe/dKwcFGpwIAoHSi5QCAYSwWadAgKSPDup2RIQ0ebGwmAABKM4oDAIaJj5dWrpTyZxjOy5NWrLDuBwAArucVxUFSUpIOHz58xX0A3EdSkrXV4Lx1CiVZtwcOtB4HAACu5RXFQWRkpOrVq3fFfQDcw/ndiS4c9WSx0L0IAACjeEVxIEmFjatmrLVzxcTEKCwsTC1atDA6CjzMrl0FuxNdKL970a5drs0FAEBp5zXFAVxvyJAh2r17tzZu3Gh0FHiY8HCpSxfJbC78uNksde1qPQ8AALgOxQEAlzOZpNhYqVKlwsccBARIs2YZkw0AgNKM4gCAIYKDrQVCYWMOYmNZ6wAAACNQHAAwTFRUwe5F+d2JoqKMzQUAQGlFcQDAMOd3L5LoTgQAgNEoDgAYKjhYmj1bql7d+pPuRAAAGKeM0QEAICqKrkQAALgDWg4AAAAASKI4AAAAAHAOxQEAAAAASRQHAAAAAM6hOAAAAAAgyUtmK+rfv786dux4xX0AAAAALs1ksVgsRoeAZ0tPT1dgYKDS0tIUEBBgdBwAAADYiW5FAEql+HipRg1p2TKjkwAA4D5oOUCJ0XIAT5OUJIWGSqmpUlCQtHcvKzMDACDRcgCglLFYpEGDpIwM63ZGhjR4sLGZAABwFxQHAEqV+Hhp5UopL8+6nZcnrVhh3Q8AQGlHtyKUGN2K4CnyuxOlpVlbEPKZTFJgIN2LAADw6JaDIUOG6IMPPlBOTo7RUQC4ufO7E134lYjFQvciAAAkD2858PHxkclk0gMPPKDly5fLbDYbHalUouUAnmDnTqlRo6KdFx7u/DwAALgjj245kCSLxaJPP/1UvXr10uXqnBkzZui9995zYTIA7iQ8XOrSRbrUdwhms9S1K4UBAKB08/iWg6ZNm8rf31/r169XdHS05s+fX+i5bdq00fr165WXPwoRDkPLATzF5cYcBAVJe/Yw5gAAULp5fMtBhQoV9MUXX6h58+aKi4vTY489ZnSkUiMmJkZhYWFq0aKF0VGAIgkOlmJjCx9zEBtLYQAAgMcXB5JUqVIlrV69Wo0aNdKcOXM0bNgwoyOVCkOGDNHu3bu1ceNGo6MARRYVVbB7UX53oqgoY3MBAOAOvKI4kKTKlStrzZo1atCggWJiYvTcc88ZHQmAGzKZrK0ElSpZtwMCpFmzjM0EAIC78JriQJKqVq2qtWvXqn79+poyZYrGjBljdCQAbig4WJo9W6pe3fqT7kQAAFh5VXEgSVdffbXWrl2ra6+9Vq+88opee+01oyMBcENRUdLx41K3bkYnAQDAfZQxOoAz1KxZU2vXrlXbtm01atQo+fn5GR0JAAAAcHseXRzMnTtXp06dKvRY3bp1bQXCs88+qwoVKrg4HQAAAOBZPHqdg6LYvXu32rVrp7///lsmk4l1DpyAdQ4AAAC8g9eNObhQWFiY1qxZo5tvvln+/v5GxwEAAADclte3HJzPYrHIZDIZHcPr0HIAAADgHby+5eB8FAYAAADApZWq4gAAAADApVEcAIALxcdLNWpIy5YZnQQAgIuVqjEHcA7GHABFk5QkhYZKqalSUJC0dy+rMwMA3AstBwDgAhaLNGiQlJFh3c7IkAYPNjYTAAAXojgAABeIj5dWrpTyl1rJy5NWrLDuBwDAXdCtCCVGtyLg8vK7E6WlWVsQ8plMUmAg3YsAAO6jTHHfYDabnZHDo5hMJuXm5hodA4AHOL870YVfxVgs/3YvWr7cmHwAAJyv2C0HPj70RDKZTMrL7xsAWg6Ay9i5U2rUqGjnhYc7Pw8AAJdT7JYDqXQvJkYvLADFER4udekiJST8O97gfGaz1LkzhQEAwD3QDAAATmQySbGxUqVK1j9feCwgQJo1y5hsAABcyK6WA749B4CiCw62Fgg9ehTcb7FY9zMYGQDgLopdHIwZM8YZOQDAq0VFSUuX/tu9KL87UVSU0ckAAPgXU5mixBiQDBTN+SskV64s7dlDqwEAwL0w5gAAXCQ4WJo9W6pe3fqTwgAA4G5oOUCJ0XIAAADgHWg5AAAAACCJ4gAAAADAORQHAAAAACRRHAAAAAA4h+IAAAAAgCSKAwAAAADnUBwAAAAAkERxAAAAAOAcjywOfHx8VKZMGaNjAIDHiI+XatSQli0zOgkAwJ15ZHEgSSzsbLyYmBiFhYWpRYsWRkcBcBlJSdLAgdKff0qPPWbdBgCgMCaLB37K9vHxkclkUl5entFRICk9PV2BgYFKS0tTQECA0XEAnMdikSIjpYQEKS9PMpulzp2l5cuNTgYAcEce23IAALiy+Hhp5UprYSBZf65YYd0PAMCFitVykP+NfZEubDIpNzfX7mBFyUHLgXug5QBwT0lJUmiolJZmbUHIZzJJgYHS3r1ScLBx+QAA7qdYo3pHjx5d5OIAAGAci0UaNEjKyChYGOQfy8iQBg+mexEAoCDGHKDEaDkA3M/OnVKjRkU7Lzzc+XkAAJ6BMQcA4IXCw6UuXawDkAtjNktdu1IYAAAKKlHLQXJyshYvXqy9e/fq1KlTBS9sMmnOnDklDlgYWg7cCy0HgHu63JiDoCBpzx7GHAAACrJ7JbEDBw6oVatWys7OVlZWlqpWraqUlBTl5eWpcuXKCgwMdGROAEAxBQdLsbFSjx4F91ss1v0UBgCAC9ndrWjkyJGKiIjQiRMnZLFYtGrVKmVmZmrGjBny8/PTp59+6sicAAA7REUV7F6U350oKsrYXAAA92R3cbBhwwY9/vjj8vPzk2RdsdjX11dDhgzRgAED9MwzzzgsJADAPiaTtZWgUiXrdkCANGuWsZkAAO7L7uIgOTlZV199tXx8fFSxYkWlpqbajrVt21br1693RD4AQAkFB0uzZ0vVq1t/0p0IAHApdhcHISEhSklJkSTVqVNH33//ve3Yjh07VL58+ZKnAwA4RFSUdPy41K2b0UkAAO7M7gHJbdq00fr16/Xggw+qX79+GjlypI4ePaoyZcpo/vz56nHhCDgAAAAAbs3u4mDUqFE6duyYJGnYsGE6ePCglixZIkn6z3/+o6lTpzomIQAAAACXYIVklBjrHAAAAHgHVkgGAAAAIKkExcFTTz2l6OjoQo9FR0drxIgRdocCAAAA4Hp2FwcJCQnq1KlTocfuuusuffzxx/ZeGgAAAIAB7C4Ojh07pmuvvbbQY7Vq1dLRo0ftDgUAAADA9ewuDipVqqRDhw4Veuz333+Xv7+/3aEAAN4hPl6qUUNatszoJACAorB7tqJu3bpp586d2rhxoypWrGjb/88//6hly5YKCwvTRx995LCg52O2IvfCbEUACpOUJIWGSqmpUlCQtHcvqzMDgLuzuzj45Zdf1LJlS1WoUEG9evXSNddco6NHj2rJkiXKzMzUTz/9pLCwMEfnlURx4G4oDgBcyGKRIiOlhAQpL08ym6XOnaXly41OBgC4nBKtc7B161Y9/fTT+t///qe8vDyZzWa1bdtWU6ZMUZMmTRwYsyCKA/dCcQDgQkuXSj16FL4/Ksr1eQAAReOQRdBOnTqlkydPqkqVKvLz83NErsuiOHAvFAcAzpffnSgtzdqCkM9kkgID6V4EAO7MIYuglS9fXiEhIS4pDAAA7stikQYNkjIyChYG+ccyMqTBg43JBgC4smIVB2azWT///LP1jT4+MpvNl3yVKVPGKYEBAO5r1y5p5UrrOIPC5OVJK1ZYzwMAuJ9ifYIfPXq0atasafuzyWRySigAgGcKD5e6dPl3IPKF8gcmh4e7PhsA4MocMubA1Rhz4F4YcwDgfJcbcxAUJO3Zw5gDAHBXdo85OHz4sHJycgo9lpubq8OHD9sdCgDguYKDpdjYwsccxMZSGACAO7O7OKhbt662bt1a6LHt27erbt26docCAHi2qChr9yKz2bptNktduzKNKQC4O7uLg8v1RsrJyZGPj0MmQgIAeCCTydpKUKmSdTsgQJo1y9hMAIArK9aA5OTkZCUlJdm2f/vtN1WsWLHAOadOndL777+v6tWrOyYhAMAjBQdLs2dLw4dL06fTnQgAPEGxioOYmBiNGzdOJpNJJpNJvXr1KvQ8i8WicePGOSQgAMBzRUXRlQgAPEmxioOHHnpIderUkcViUf/+/TVq1CjVr1+/wDnlypVTWFiYGjdu7NCgAAAAAJyrWMVBRESEIiIiJEn79u3ToEGDFBIS4pRgAAAAAFzL7mWM33//fZlMJg0cOFDXXnutIzMBAAAAMIDdUwo98MADmjZtmurXr6/OnTtr9erVjswFAAAAwMXsLg7mzJmjo0ePavLkydq3b5/uueceXX/99ZoyZYpOnjzpyIwAAAAAXMBkudyCBcXw9ddf65133tGnn36qMmXKqHv37hoyZIhuuukmR1y+AB8fH5lMJuXl5Tn82ii+9PR0BQYGKi0tTQEBAUbHAQAAgJ0ctlJZhw4d9Pjjj6t169Y6ffq0Fi5cqFatWqljx446dOiQo24DAAAAwElKXBykpaVp2rRpuvHGG9WpUydlZWVp7ty5Sk1NVVxcnPbv368+ffo4IisAAEUSHy/VqCEtW2Z0EgDwLHZ3K9qyZYveeecdffjhh8rJyVFkZKSeeOIJ3XzzzQXO+/zzzxUZGanTp087JLBEtyJ3Q7ciAO4kKUkKDZVSU6WgIGnvXlZnBoCisnsq05tuuknVq1fXiBEjNHDgQFWvXr3Q8+rVq3dRwQAAgDNYLNKgQVJGhnU7I0MaPFhavtzYXADgKexuOViyZIm6deumMmXsri/sRsuBe6HlAIC7WLpU6tGj8P1RUa7PAwCexmGzFbkSxYF7oTgA4A7yuxOlpVlbEPKZTFJgIN2LAKAoHDZbEQAARjm/O9GFX3lZLP92LwIAXF6xWg7yv7Ev0oVNJuXm5todrCg5aDlwD7QcADDazp1So0ZFOy883Pl5AMBTFWvAwOjRo4tcHAAA4Crh4VKXLlJCglTY90Zms9S5M4UBAFwJYw5KuS+++EKzZ89Wenq6cnJylJOTo6efflrdu3cv8jVoOQDgDi435iAoSNqzhzEHAHAljDko5fr06aNOnTpp3bp1+t///qdRo0apR48eWrlypdHRAKBYgoOl2NjCxxzExlIYAEBRlKjlIDk5WYsXL9bevXt16tSpghc2mTRnzpwSBywMLQeOc9999ykhIUFms9m2Lzw8XA0aNNDyIk4MTssBAHdhsUiRkf92L8rvTsQ6BwBQNHYvUnDgwAG1atVK2dnZysrKUtWqVZWSkqK8vDxVrlxZgYGBjszpVU6fPq0JEyZo8uTJ6tWrl+bNm3fF89944w0tWrRIhw4dUuXKldWxY0eNHz9edevWLVGWzz///KJ9fn5+OnPmTImuCwBGMJmsrQTr1llXSA4IkGbNMjoVAHgOu7sVjRw5UhERETpx4oQsFotWrVqlzMxMzZgxQ35+fvr0008dmdNrrFu3To0bN9asWbOUk5NzxfNPnTqlDh06aPLkyXrllVeUkpKiL7/8Ujt27FDTpk21fft2h+ZLSUnRrl271KtXL4deFwBcJThYmj1bql7d+pPuRABQdHa3HGzYsEFTp06Vn5+fJMliscjX11dDhgzRX3/9pWeeeUZffvmlw4J6gyVLlmjw4MGaOHGiypcvr/79+1/xPWPHjtX69es1c+ZMdevWTZIUERGh5cuXKzQ0VA8//LC2b99u6xaUmpqqP//887LX9Pf317XXXlvosfHjx6t9+/bqUdgSowDgIaKiWBEZAOxhd3GQnJysq6++Wj4+PqpYsaJSU1Ntx9q2bau33nrLEfm8St26dbV7926FhIRcsSuRJGVlZSkmJka+vr7q27dvgWP169dXx44dtXr1an3xxRd64IEHJEkffvihBl9hpZ927drpm2++uWh/XFycvv32W33zzTdMWQsAAFAK2d2tKCQkRCkpKZKkOnXq6Pvvv7cd27Fjh8qXL1/ydF6mdevWCgkJKfL5a9euVWZmpiIiIlSxYsWLjt92222SpISEBNu+QYMGyWKxXPZVWGGwfPlyvfXWW/rqq68YLwIAAFBK2V0ctGnTRuvXr5ck9evXTxMnTtSAAQM0ePBgvfDCC7rvvvscFrK0yh9PcKlBx/Xq1Stwnr1WrFihV199VatXr1a1atUkSa+//nqJrgkAAADPY3e3olGjRunYsWOSpGHDhungwYNasmSJJOk///mPpk6d6piEpVj+2IEqVaoUerxy5coFzrPHihUrNGTIEH3wwQc6dOiQDh06JEmaNWuWnn/++ULfk52drezsbNt2enq63fcHAACA+7C7OKhfv77q169vvUiZMpo5c6ZmzpzpsGCwjjmQJF9f30KPlytXrsB59ujRo4dycnJ07733Fthfu3btS77ntdde07hx4+y+JwAAANwTKyS7MX9/f0m65JoD+d/e559njzNnzhQ6LuH333+/5HteeOEFpaWl2V5Hjhyx+/4AAABwH3YXB2fOnNG4cePUoEED+fv7y2w2F3iVKWN3owTOqV69uiTZBn5f6OTJk5KkGjVquCyTZG2xCAgIKPACAACA57P7E/zzzz+vGTNm6J577lHXrl1tXVzgOE2aNJEk/fbbb4UeP3jwoCSpcePGrooEALiM+Hhp+HBp+nTp3NI0AOBR7C4Oli1bprFjx2rUqFGOzIPztG/fXhUqVFBiYqIyMzNVoUKFAsd/+OEHSdKDDz5oRDwAwHmSkqSBA6XUVOmxx6R27VidGYDnsbtbUXp6um699VZHZsEF/P39NXToUGVnZ2v+/PkFjh08eFBr1qxReHj4RYOJAQCuZbFIgwZJGRnW7YwM6QrrUQKAW7K7OLjzzju1YcMGR2ZBIcaMGaObb75ZI0eO1PLly3Xq1CklJiYqMjJS/v7+WrRokcxms9ExAaBUi4+XVq6U8vKs23l50ooV1v0A4ElMFovFYs8b9+7dq8jISA0dOlT3339/oXPxl2QWncvx8fGRyWRSXv5T2IOYTKZLHps7d6769et30f7Tp0/r9ddf16JFi3T48GEFBgbqzjvv1Pjx420LoRkpPT1dgYGBSktLY3AygFInKUkKDZXS0qwtCPlMJikwUNq7l+5FADyH3cWBj4+10eFyH3ad9eG9sOIgf19RmEwm5ebmOiVbaURxAKC0slikyEgpIeHfVoPzmc1S587S8uWuzwYA9rB7QPLo0aOL/GHcFdwtT2kQExOjmJgYj2zBAQBH2LXL2p3oUvK7F+3aJYWHuy4XANjL7pYDI3lytyJvRMsBgNKKlgMA3oYVkgEAsJPJJMXGSpUqWf984bGAAGnWLGOyAYA9SrSMcXJyshYvXqy9e/fq1KlTBY6ZTCbNmTOnROFKwp2zAQC8R3CwtUDo0aPgfovFup/ByAA8id3FwYEDB9SqVStlZ2crKytLVatWVUpKivLy8lS5cmUFBgY6MqfXZAMAeJ+oKGnp0n+7F+V3J4qKMjoZABSP3d2KRo4cqYiICJ04cUIWi0WrVq1SZmamZsyYIT8/P3366aeXfX/r1q317rvv2nt7p2YDAKA4zu9eJNGdCIDnsrs42LBhgx5//HH5+flJkiwWi3x9fTVkyBANGDBAzzzzzGXf//PPP2vw4MGKjY296NiJEyeUkJBgb7QSZwMAoLiCg6XZs6Xq1a0/6U4EwBPZXRwkJyfr6quvlo+PjypWrKjU1FTbsbZt22r9+vVXvMaQIUM0bNgwLViwoMD+33//XV26dLE3mkOyAQBQXFFR0vHjUrduRicBAPvYXRyEhIQoJSVFklSnTh19//33tmM7duxQ+fLlr3iN6OhoxcTE6NFHH9VHH31kbxSnZAMAAABKG7sHJLdp00br16/Xgw8+qH79+mnkyJE6evSoypQpo/nz56vHhdM2XMKAAQOUlZWl6OholS9fXvfdd5+9kRyeDQAAAChN7C4ORo0apWPHjkmShg0bpoMHD2rJkiWSpP/85z+aOnVqka81fPhwZWZmqlu3bvrss89UoUIFe2M5PBsujRWSAQAAvIthKyT7+Pjop59+UsuWLW37XnzxRc2YMUNjx47Vc889d8kPnayQ7F5YIRkAAMA7lGgRtJLYsmWL6tatW2DfxIkTlZmZqeeee86gVAAAAEDpZfeA5JJq0qRJoYuRTZs2TUOHDlWl/Mmi7fDUU08pOjq60GPR0dEaMWKE3dcGAAAAvJXDi4Pu3btr0qRJWr16tZKSkuy6xrRp0wpMP1pcCQkJ6tSpU6HH7rrrLn388cd2XxsAAADwVg7vVrRs2bIC05JWr15dTZs2LfC6sDuRox07dkzXXnttocdq1aqlo0ePOvX+AAC4Sny8NHy4NH066ysAKDmHFwdvv/22tmzZos2bN2vPnj06fvy4jh8/rlWrVtnOCQgIUJMmTdS0aVM1adJEffr0cWiGSpUq6dChQ4Ue+/333+Xv7+/Q+wEAYISkJGngQCk1VXrsMaldO1ZmBlAyTp2tKCsrS9u2bdPmzZu1adMmrV271vatvclkksVisWvWoSvNVtStWzft3LlTGzduVMWKFW37//nnH7Vs2VJhYWEOXXSttGO2IgBwPYtFioyUEhKkvDzJbJY6d5aWLzc6GQBP5tKpTC0Wi9atW6eBAwfqxIkT6tOnj/bv368vv/yyWNe5UnHwyy+/qGXLlqpQoYJ69eqla665RkePHtWSJUuUmZmpn376SWFhYY74lSCKAwAwwtKlUmFrei5dKkVFuT4PAO9gyDoH6enp6tixo4KDg/XZZ58V+/1FWedg69atevrpp/W///1PeXl5MpvNatu2raZMmaImTZqUID0uRHEAAK6VlCSFhkppadYWhHwmkxQYKO3dS/ciAPYxbBG0zZs3q0WLFvrwww8VVcyvOIqzCNqpU6d08uRJValSRX5+fvbGxWVQHACA61zYnehCdC8CUBKGFQeSVLNmTTVo0EBr1qwp1vtYIdk9xMTEKCYmRnl5edq3bx/FAQC4wM6dUqNGRTsvPNz5eQB4F4cXB6mpqQoKCirSuQ0aNFBycrKSk5OLdY/CigOz2awff/xRLVu2tB2/FJPJpNzc3GLdE5dGywEAuA4tBwCcyeFTmV511VWqW7eumjVrpmbNmql58+Zq3ry5qlSpUuC8w4cP67ffflOFChUcct/Ro0erZs2atj9frjgAAMBTmUxSbKy0bl3hYw4CAqRZs4zLB8CzObzlwGw2K/+S539Ar1Wrlpo1a6aGDRvKZDJp4cKFOnTokNq2bat169YV6x50K3IvtBwAgOsxWxEAZ3B4cZCVlaWtW7dqy5Ytttcvv/xi68aTXzBYLBb5+flpzZo1uuWWW4p1jysVB4cPH1aNGjVUtmzZi47l5uZedgVlFB/FAQC4HuscAHAGlwxIzs7OVmJiorZs2aLdu3crNTVVderUUe/evXXdddcV+3pXKg7OH39woc2bN6tly5a0OjgQxQEAGCN/StPUVKlyZWnPHqYwBVAyDh9zUJhy5cqpRYsWatGihStup8vVOzk5OfLx8XFJDgAAnCk4WJo9Wxo+XJo+ncIAQMm5pDhwheTkZCUlJdm2f/vtN1WsWLHAOadOndL777+v6tWruzoeAABOERXFGAMAjuM1xUFMTIzGjRsnk8kkk8mkXr16FXqexWLRuHHjXJwOAAAAcH9eUxw89NBDqlOnjiwWi/r3769Ro0apfv36Bc4pV66cwsLC1LhxY4NSAgAAAO7La4qDiIgIRURESJL27dunQYMGKSQkxOBUAAAAgOfwmuLgfO+//75MJpMGDhzIlKUAAABAEXnltD0PPPCApk2bpvr166tz585avXq10ZEAAAAAt+eS4uDw4cM6evSoK24lSZozZ46OHj2qyZMna9++fbrnnnt0/fXXa8qUKTp58qTLcni7mJgYhYWFuWyKWgAAADiXSxZB8/HxUbVq1XTixAmHXe9yi6Bd6Ouvv9Y777yjTz/9VGXKlFH37t01ZMgQ3XTTTQ7JU9qxCBoAAIB3cFm3IhfUIJfUoUMHPf7442rdurVOnz6thQsXqlWrVurYsaMOHTpkWC4AADxVfLxUo4a0bJnRSQA4kleOOciXlpamadOm6cYbb1SnTp2UlZWluXPnKjU1VXFxcdq/f7/69OljdEwAADxKUpI0cKD055/SY49ZtwF4B6+crWjLli1655139OGHHyonJ0eRkZH64IMPdPPNN9vO6dmzpwICAhQZGWlgUgAAPIvFIg0aJGVkWLczMqTBg6Xly43NBcAxvLI4uOmmm1S9enWNGDFCAwcOVPXq1Qs9r169egUKBgAAcHnx8dLKlf9u5+VJK1ZY90dFGZcLgGO4bEBy1apVleSgdscrDUhesmSJunXrpjJlvLL2cTsMSAaA0iEpSQoNldLSrC0I+UwmKTBQ2rtXCg42Lh+AkvPKMQc9e/akMAAAwIHO70504deKFsu/3YsAeDavLA4AAIBj7dpl7U50qVnE87sX7drl2lwAHMtrigMfHx+ZzeYivWhVAACgeMLDpS5dJLO58ONms9S1q/U8AJ7Laz4ljx49WiaTyegYAAB4JZNJio2V1q0rfMxBQIA0a5Zx+QA4hlcOSIZrMSAZAEqPpUulHj0K389sRYDn85puRQAAwPmiogp2L8rvTkRhAHgHr+lWdKHk5GQtXrxYe/fu1alTpwocM5lMmjNnjkHJAADwXOd3L0pNpTsR4G28sjg4cOCAWrVqpezsbGVlZalq1apKSUlRXl6eKleurMDAQKMjAgDc1O+//662bdtqxIgReuKJJ4yO45aCg6XZs6Xhw6Xp01nbAPAmLulWNGbMGI0YMcIVt5IkjRw5UhERETpx4oQsFotWrVqlzMxMzZgxQ35+fvr0009dlgUA4DmSk5PVqFEjHTlyRMOGDdOECRN09uxZo2O5pago6fhxqVs3o5MAcCSXDEh2tCsNSL722ms1depUdenSRWXKlNHGjRvVvHlzSdLYsWP1008/6csvv3RlZK8UExOjmJgY5eXlad++fQxIBuDxkpOTVa9ePWVkZNj23XnnnZo7d66uueYaA5MBgGt45YDk5ORkXX311fLx8VHFihWVmppqO9a2bVutX7/euHBeZMiQIdq9e7c2btxodBQAcIiqVavql19+UYUKFWz7/vvf/6pRo0ZatmyZgckAwDW8sjgICQlRSkqKJKlOnTr6/vvvbcd27Nih8uXLGxUNAODmrrnmGn366acF1s45efKkoqKi1Lt3b6WlpRmYDgCcyyuLgzZt2thaB/r166eJEydqwIABGjx4sF544QXdd999BicEALiz9u3b6/nnn79of1xcnBo3bqxvv/3WgFQA4HxeOebgwIEDOnbsmNq0aaPc3Fw9+eSTWrJkiSTpvvvu04wZM5ixyIFYBA2ANzpz5oxuvfVWbdq0SZJUtmxZ5eTkSLJOif3ss89qwoQJKleunJExAcChvLI4gGtRHADwVr/++quaNm2qzMxMSVJYWJh2795tO964cWPFxcWpUaNGRkUEAIfyym5FAAA4wvXXX68ZM2bYto8fP66XXnpJvr6+kqTExETddNNNeuutt5jyFIBX8Mri4MyZMxo3bpwaNGggf39/mc3mAq8yZbxy7TcAgBP069dP3c5N5n/y5En973//008//aSGDRtKsv6d88wzz6hjx446cuSIkVFLhfh4qUYNicmjAOfwym5FTz31lGbMmKF77rlHjRo1KrQ/6JgxY5wds9SgWxEAb3fy5ElFRETYPvxPnDhRTz31lF566SW99dZbtvMCAwP1zjvvqFevXkZF9WpJSVJoqJSaKgUFSXv3sjoz4GjFKg7yP5QX6cImk3Jzc+0OVpQclyoOatasqUGDBmnUqFFOuT8KojgAUBp89913at++vc6ePasyZcrohx9+UMuWLbV27Vr169evQKtBjx499M4776hy5coGJvYuFosUGSklJEh5eZLZLHXuLC1fbnQywLsUqzgYO3ZskYsDyXnfzl+pOAgICNAnn3yi9u3bO+X+KIjiAEBpMWrUKL366quSpPr162vr1q2qVKmSUlNTNXToUC1atMh27jXXXKP58+erQ4cORsX1KkuXSj16FL4/Ksr1eQBv5ZXdiiIjI9WiRQuNHDnSxclKJ4oDAKVFTk6O2rRpow0bNkiyjkeYO3eu7fjSpUs1aNAgpaam2vY9+eSTmjhxIgtwlkB+d6K0NGsLQj6TSQoMpHsR4EheOSB54sSJiouLU2xsrP744w9lZWVd9AIAoLjKli2rRYsWqWLFipKkefPmaenSpbbj3bt3144dOwq0Frz99tu66aabtG3bNlfH9QoWizRokJSRUbAwyD+WkSENHmxMNsAblajlIDk5WYsXL9bevXt16tSpghc2mTRnzpwSByzMlVoOfHx8bBkuhTUSHIeWAwClzYIFC9S3b19J1kHI27dvV+3atW3Hz549qxkzZuj5559Xdna2JGthMWHCBD377LMym82G5PZEO3dKRVlGYudOKTzc+XkAb2d3cXDgwAG1atVK2dnZysrKUtWqVZWSkqK8vDxVrlxZgYGBOnjwoKPzSrpycVCUsRHMVuQ4FAcAShuLxaKHH35YS5YskSS1adNG69atu+hD/65duxQdHV2g1aBNmzZasGCB6tSp48LEnuvCgcgXYmAy4Fh2FwfdunVTSkqKPv30U1WsWFGbNm1Sw4YN9d5772nixIn66quvFO6kEp4Vkt0LxQGA0igtLU0RERE6dOiQJGnChAmFzpKXnZ2tMWPG6I033lD+X7mVKlXSjBkz1KdPn2JN9FFaXW7MQVCQtGcPYw4AR7F7zMGGDRv0+OOPy8/PT5L1WxRfX18NGTJEAwYM0DPPPOOwkAAAuJvAwEAtWrTI1pV17Nix+vHHHy86r1y5cpo0aZK+/fZbW9ejjIwM2+Jqf//9t0tze6LgYCk2tvAxB7GxFAaAI9ndcuDv76+vvvpKt912mwICArRy5UrbAKy1a9fqoYceUnp6ukPD5itKy4FR4yFKI1oOAJRmY8eO1bhx4yRJdevW1bZt2y75LExPT9ewYcM0f/58274aNWrogw8+0N133+2SvJ6KdQ4A17C75SAkJEQpKSmSpDp16uj777+3HduxY4ehU7YdOHBADRo00EsvvaTY2Fh9/vnnWrhwoebNm6dPPvlE69atu+z7W7durXfffddFaQEAnmzUqFG65ZZbJEm//fabhg4deslzAwICNG/ePH300UeqUqWKJOn48eO65557NHToUGbTuwyTydpKUKmSdTsgQJo1y9hMgDeyuzho06aN1q9fL8k6z/PEiRM1YMAADR48WC+88ILuu+8+h4UsrpEjRyoiIkInTpyQxWLRqlWrlJmZqRkzZsjPz0+ffvrpZd//888/a/DgwYqNjb3o2IkTJ5SQkOCs6B4lJiZGYWFhatGihdFRAMAwZcqUUVxcnK21YOHChVq8ePFl3xMZGamdO3cWaC2IiYlR8+bNtXnzZqfm9WTBwdLs2VL16tafdCcCnMBip/3791u+++47i8ViseTk5FiGDBliqVKliqVKlSqW3r17W1JTU+299BWZTCaLj4/PJY/XqlXL8tFHH1ny8vIsJpPJsmnTJtuxMWPGWO66664rXv+JJ56wlC1b1jJ//vwCx3766afL3rs0SktLs0iypKWlGR0FAAyzaNEiiySLJEtAQIDl4MGDV3zP2bNnLTNnzrT4+fnZ3lumTBnLhAkTLDk5OS5IDQAF2d1yUL9+fbVp00aS9VuTmTNn6u+//9bff/+tBQsWKDAw0BG1i12Sk5N19dVXy8fHRxUrViywUmXbtm1tLR6XEx0drZiYGD366KP66KOPnJgWAOANevXqpejoaEnWsQXR0dHKzc297HtMJpOGDBmirVu3qnnz5pKk3Nxcvfzyy2rbtq0OHDjg9NwAcD6vXCHZUeMhBgwYoMmTJys6Olqff/65U7ICALxHTEyM6tatK0lav369XnnllSK9r0GDBvrxxx81atQo2+xHP/74o5o0aaI5c+bYpkAFAGezuzh46qmnbN+QXCg6OlojRoywO1RJOXI8xPDhwzV69Gh169ZNa9eudVZkAIAXCAgI0KJFi2yLoU2YMEE//PBDkd6bv4Ly999/r3r16kmS/vnnHz366KPq0qWLkpKSnJYbAPLZXRwkJCSoU6dOhR6766679PHHH9t76RIbNWqUrQAYNmyYHnvsMa1YsULx8fH6z3/+o6lTpxbrei+++KKefPJJde7cWf/73/+cERkA4CVuvvlmjRkzRpJ09uxZPfzww0pLSyvy+2+55RZt27ZNjz76qG3fJ598okaNGtGKDcDp7F7noHz58lq1apVuv/32i4598803uvfee502JZuzV0jetm2b6tate9G4ieHDh2vmzJmSxOrM52GdAwAoKC8vT7fffrvtC6WePXtq0aJFxV4N+ZNPPtGAAQP0119/2fYNHDhQU6ZMUYUKFRyaGQCkErQcVKpUybZk/IV+//13+fv72x3KaE2aNCl0QPW0adM0dOhQVcqfZBkAgEKYzWbFxcXZ/i5ZsmSJ4uLiin2dzp07a8eOHbr//vtt+2bPnq0mTZpow4YNDssLAPnsLg7atWunSZMm6Z9//imw/59//tEbb7xRaIuCN5g2bVqB2Y8AAChM7dq1NXv2bNv2kCFD7Jp96Oqrr1ZCQoJmz55t++Jt//79uvXWWzV27Fjl5OQ4LDMuLT5eqlFDWrbM6CSAc9ndreiXX35Ry5YtVaFCBfXq1UvXXHONjh49qiVLligzM1M//fSTwsLCHJ1Xkv3dirp3766mTZvaXsGsnuIQdCsCgEvr16+f5s+fL0lq1aqVvv/+e5UtW9aua/3666/q3bt3gVaDFi1aKC4uTjfccIND8uJiSUlSaKiUmioFBUl797IAG7yX3cWBJG3dulVPP/20/ve//ykvL09ms1lt27bVlClT1KRJEwfGLMje4iD/ffmqV69eoFho2rSpbQo6FB3FAQBcWkZGhpo2bWprNRg1apQmTJhg9/Vyc3M1ceJEjR8/3vb3oL+/v6ZMmaKBAwcWe1wDLs9ikSIjpYQEKS9PMpulzp2l5cuNTgY4R4mKg3ynTp3SyZMnVaVKFfn5+Tki12XZWxxMnz5dW7Zs0ebNm7Vnzx7b+89/kAYEBKhJkyZq2rSpmjRpoj59+jg0uzeiOACAy/v555916623Kjc3Vz4+Plq3bp3atm1b4mtGR0fr119/te279957NWfOHFWvXr2kkXHO0qVSjx6F74+Kcn0ewNkcUhy4miNmK8rKytK2bdu0efNmbdq0SWvXrtXRo0clWYsFi8Xi1BmRvAnFAQBc2cSJE/XSSy9JkmrVqqXt27ercuXKJbpmZmamRowYoVmzZtn2Va1aVe+9954eeuihEl0b/3YnSkuztiDkM5mkwEC6F8E7ldri4EIWi0Xr1q3TwIEDdeLECfXp00f79+/Xl19+6bB7eCuKAwC4sry8PHXo0EHffvutJCkqKkoffvihQ7oBffHFF+rfv79OnDhh29e/f3+9/fbbzLBnpwu7E12I7kXwVsUqDsxms3788Ue1bNnyov77F13YZFJubq5DQl7ImescpKenq2PHjgoODtZnn33m8Ot7I4oDACiaI0eOKCIiQidPnpQkzZ07V/369XPItf/66y899thjBRYhrVu3rhYuXKhbb73VIfcoTXbulBo1Ktp54eHOzwO4SrGKg3HjxmnAgAEKCQnRuHHjrnh+/gqRjubsRdA2b96sFi1a6MMPP1QUHQqviOIAAIruo48+Urdu3SRJFSpU0NatW3X99dc75NoWi0Xz5s3TsGHDbFON+/j4aOTIkRozZox8fX0dcp/SgJYDlFZ2dSvKzs7WLbfcookTJ+quu+5yRq7LcnZxIEk1a9ZUgwYNtGbNGqfdw1tQHABA8Tz66KOaM2eOJOtUpD/88IPd05sW5uDBg+rTp49++OEH275mzZopLi5ON954o8Pu4+0uN+YgKEjas4cxB/A+di2CVq5cOR04cMDjvoEozuJlFStW1LZt25yWBQBQer399tu21oKNGzc6vKW9Xr16+vbbb/Xqq6+qTJkykqQtW7aoWbNmmjFjhs6ePevQ+3mr4GApNrZgYSBZt2NjKQzgnUq0QnL+oCpPcdVVV+m6665TVFSUJk2apP/+979KSUm56LzDhw/rt99+4+EJAHCKihUravHixbYP7pMmTdK6desceg+z2awXX3xRGzZssLUWnD59WsOGDdM999yjY8eOOfR+3ioqSurSxdqNSLL+7NqVaUzhveyerWjHjh166KGHNHDgQHXt2lU1atS4aIBy/jLvjmZvtyKz2az8X/f8rLVq1VKzZs3UsGFDmUwmLVy4UIcOHVLbtm0d/rD2RnQrAgD7vPHGG3r++eclSddcc40SExNVpUoVh9/n1KlTev755zVjxgzbvipVqmj27Nn6z3/+4/D7eZvzV0iuXJnuRPBudhcHPj7WRofLzVjkrDEB9hYHWVlZ2rp1q7Zs2WJ7/fLLL7ZZlfJ/F4vFIj8/P61Zs0a33HKLw/N7G4oDALDP2bNndeedd2rt2rWSpK5du+qjjz5y2irHX331lfr166fjx4/b9vXu3VszZsxQYGCgU+7pLeLjpeHDpenTpXPjyQGvZHdxMHbs2Cs+vDxhtqLs7GwlJiZqy5Yt2r17t1JTU1WnTh317t1b1113nQPSej+KAwCw39GjR9W4cWNbN9f33ntPjz76qNPu9/fff2vw4MFatmyZbd+1116rhQsXlnjVZgCej0XQUGIUBwBQMitXrlTXrl0lWbvkbtmyRaGhoU67n8ViUVxcnIYOHar09HRJ1tbzZ599VhMmTFC5cuWcdm8A7s3uAckAAMAxunTposcee0yStQtsr169dObMGafdz2QyqXfv3kpMTFS7du0kWQuGyZMnq1WrVtq5c6fT7g3AvTmkODh16pSOHj3qtBWRAQDwdm+99ZattWDLli0aNWqU0+9Zu3Ztff3113rjjTds6yxs375dzZs311tvvcWsfUApVKLi4KuvvlKrVq1UqVIl1a5dW4mJiZKkxx9/vEBfRninmJgYhYWFqUWLFkZHAQCPV6FCBS1ZssT2IX3y5MkuWYjTbDZrxIgR2rhxoxo2bChJOnPmjJ555hndeeedOnLkiNMzAHAfdhcHq1at0r333isfHx+NGjWqwLcLNWrU0Ny5cx0SEO5ryJAh2r17tzZu3Gh0FADwCk2bNtVrr71m2+7Tp4+Sk5Ndcu+IiAht3LhRTz/9tG3f2rVr1bhxYy1ZssQlGQAYz+7iYMyYMerWrZt+/PHHi5o+GzVqZGtFAAAARffUU0/pzjvvlCQdP35cjz76qFw1d4ifn5+mTJmir7/+WjVr1pQkpaamqlevXurZs6dOnjzpkhwAjGN3cbBjxw498sgjki5e66BKlSou+6YDAABv4uPjo/nz56tq1aqSpE8++UTvvvuuSzPccccdSkxMVK9evWz7PvzwQzVu3Fhff/21S7MAcC27iwN/f3/b9GcX+uOPP5yywiMAAKVBjRo1NGfOHNv2U089pV9++cWlGSpXrqxFixZpyZIlCgoKkmT9+71jx456+umndfr0aZfmKY3i46UaNSSGccKV7C4O2rZtqzfffFPZ2dm2fSaTSRaLRe+99546dOjgkIAAAJRGDz74oB5//HFJ1lkBe/bsWeDvXFfp0aOHEhMTdccdd9j2TZ06VTfddJO2bdvm8jylRVKSNHCg9Oef0mOPWbcBV7C7OJgwYYJ27dqlhg0b6qWXXpLJZNJ7772n2267TVu2bNHo0aMdmRMAgFLnzTffVFhYmCTrFKMvvviiITlq1aql//73v5o6daptgbRdu3apZcuWeuONN1iU1MEsFmnQICkjw7qdkSENHmxsJpQeJVohedu2bXr22Wf13XffKTc3Vz4+Pmrbtq2mTp2qiIgIR+YsgBWS3QsrJAOA8yQmJqpFixa2RdFWr16tTp06GZZn586dio6O1vbt22372rZtq/nz56tOnTqG5fImS5dKPXoUvj8qyvV5ULqUqDjId/r0aaWkpCgoKEj+/v6OyHVZFAfuheIAAJxr2rRpevLJJyVJ1atXV2JioqpVq2ZYnuzsbI0ePVqTJ0+2zaRUqVIlzZw5U717975oohIUXVKSFBoqpaVZWxDymUxSYKC0d68UHGxcPni/EhUHf/31lz7++GMdOXLkooFJJpNJr7/+eokDFobiwL1QHACAc1ksFt1777368ssvJUn333+/EhISDP8Q/t1336lPnz46dOiQbV9kZKRmz56tq666ysBknslikSIjpYQEqbCPOGaz1LmztHy567Oh9LC7OPjyyy8VGRmpU6dOFX5hJ354pzhwLxQHAOB8J06cUKNGjfTXX39JkmbOnKkhQ4YYnEpKS0vTsGHDtGDBAtu+/MVQ77rrLgOTeZ6dO6VGjYp2Xni48/OgdLJ7QPKIESPUpEkTbdu2TdnZ2Tp79myBFx/cAQBwnKuvvlrz5s2zbT/77LPatWuXcYHOCQwM1Pz587Vs2TLbNObHjx/X3XffrSeeeEJZWVkGJ/Qc4eFSly7WFoLCmM1S164UBnAuu1sOKlSooJUrVxoyKIqWA/dCywEAuM6wYcM0Y8YMSVKjRo30888/y8/Pz+BUVseOHVP//v21evVq274GDRooLi5OzZs3NzCZ57jcmIOgIGnPHsYcwLnsbjlo1KiR/v77b0dmAQAAV/DGG2+oYcOGkqQdO3bo+eefNzjRv0JCQrRq1SrNmDHDVrDs2bNHrVu31quvvqrc3FyDE7q/4GApNrZgYSBZt2NjKQzgfHa3HHz77bcaMmSIVqxYoRtuuMHRuS6LlgP3QssBALjWzp07ddNNN9kWRfviiy90zz33GJyqoD179ig6OlqbN2+27bvlllu0YMEC1a9f38Bk7u/CgckMRIYrlWi2opEjR2rKlCmqU6eOKleuXPDCJpM2bNhQ4oCFoThwLxQHAOB6M2fO1BNPPCFJCg4OVmJioq6++mqDUxV05swZjR8/Xq+99prOnj0rSapYsaLefvtt9e/f3/DZltxZfvei1FSpcmW6E8F17O5W9Oqrr+qNN95QlSpVFBwcrAoVKhR4uWK9AwAASqshQ4bovvvukyQlJSXpkUcekQOWLnIoX19fvfLKK/r+++9Vr149SdI///yjRx99VF27drXNvISLBQdLs2dL1atbf1IYwFXsbjmoXr26HnjgAcXGxsp8qWH1TkLLgXuh5QAAjJGUlKTGjRvrxIkTkqyLpQ0bNszgVIXLyMjQU089pTlz5tj2BQcH64MPPrAVOQCMZ3fLQWZmpnr27OnywgAAAFgFBwdr/vz5tu3nnntOiYmJBia6tEqVKun999/XypUrVbVqVUnW4ub+++/XoEGDlJmZaXBCAFIJioMOHTpo69atjswCAACK6a677tKTTz4pScrOzlbPnj0vuUCpO3jooYe0c+fOAq0Fs2fPVpMmTZw2VhFA0dldHLz22muaO3eu3nvvPf3xxx/Kysq66AUAAJxv0qRJioiIkCTt3r1bI0aMMDjR5V199dX69NNPFRsbaxujuH//ft16660aO3ascnJyDE4IlF52jznw8bHWFZebacBZYwIYc+BeGHMAAMbbvXu3mjdvrtOnT0uSPv30U91///0Gp7qyffv2qXfv3vr5559t+1q2bKmFCxe6fKp0ACUoDsaOHXvFKcjGjBljV6groThwLxQHAOAeZs2apccff1ySVLVqVSUmJqpGjRoGp7qy3Nxcvfrqq5owYYLt73Z/f39NmTJFAwcOZMpTwIVKtM6BUSgO3AvFAQC4B4vFooceekgJCQmSpE6dOmnVqlW21n53t2HDBvXu3Vu//vqrbd99992n999/X9WrVzcwGVB6eMbTAgAAXJHJZNKcOXNsrQVfffWVpk2bZnCqomvVqpW2bt2qQYMG2fZ9/vnnatSokT7++GPjggGlCMUBAABepGrVqgWmNx05cqS2bdtmXKBiqlChgmbNmqXPPvvMtuJzcnKyunTpov/7v/9TRkaGwQm9X3y8VKOGtGyZ0UlgBLu7FdWtW/eyfQBNJpMOHDhgd7DLoVuRe6FbEQC4nxEjRujNN9+UJDVo0ECbN2+2zQzkKf766y8NGDBAn3zyiW1fvXr1tHDhQt1yyy0GJvNeSUlSaKiUmioFBUl797I6c2ljd8tBu3btLno1bNhQaWlpslgsateunSNzAgCAYnj11VfVtGlTSdKePXv0zDPPGJyo+KpVq6aVK1dqzpw5qlixoiTp4MGDatOmjV566SWdOXPG4ITexWKRBg2S8htnMjKkwYONzQTXc/iA5JSUFN1111169tln1b17d0de2oaWA/dCywEAuKc9e/aoWbNmtkXRVq5cqYceesjYUHY6ePCgevfurfXr19v2NWvWTHFxcbrxxhsNTOY9li6VevQofH9UlOvzwBgOH3NQpUoVjRgxQuPGjXP0pQEAQDE0aNCgwIDkRx99VMeOHTMwkf3q1aun7777Tq+++qrKlCkjSdqyZYuaNWummTNnygMnX3QrSUnWVoMLe4ybTNLAgdbjKB2cMiA5MDBQhw4dcsalAQBAMTz66KPq0qWLJOnvv/9Wnz59dPbsWYNT2cdsNuvFF1/UTz/9pAYNGkiSTp8+rSeeeEJ33323xxY+Rju/O9GFNZbFQvei0sbhxUFycrImT56s66+/3tGXBgAAxWQymfTee+/pmmuukSR9/fXXmjJlisGpSqZ58+bavHmznnjiCdu+r776So0aNdJHH31kYDLPtGuXtHKldKne2nl50ooV1vPg/Rw6W1F2draSkpJkMpn08ccf695773VIyAsx5sC9MOYAANzf2rVr1bFjR1ksFpUtW1Y//vijmjdvbnSsElu9erUeeeQRHT9+3LavT58+mj59ugIDAw1M5jksFikyUkpIKLxAMJulzp2l5ctdnw2uZ3dx0K9fv4uKAz8/P9WuXVtRUVGqV6+eQwIWhuLAvVAcAIBnGDlypF5//XVJ0g033KAtW7aoQoUKBqcqub///luDBg0q0GpQu3ZtLViwQG3btjUwmefIn8I0La1g1yKTyTql6Z49TGlaWjh8tiJXoDhwDzExMYqJiVFeXp727dtHcQAAbu7MmTO65ZZbtHnzZknW8Qjvvfeewakcw2KxKC4uTkOHDlV6eroka5eqESNGaPz48SpXrpzBCd0fsxVBojiAA9ByAACeY9++fWrWrJkyMzMlSR999JEiIyMNTuU4hw4dUp8+ffTdd9/Z9kVERCguLk4NGzZ0+v1Pnz6tCRMmaPLkyerVq5fmzZvn9Hs6yoXdi+hOVDrZPSD5zJkzGjdunBo0aCB/f3+ZzeYCr/xpxgAAgPu44YYbNH36dNv2gAEDdOTIEQMTOVbt2rW1du1avf766ypbtqwkafv27brppps0depUp87UtG7dOjVu3FizZs1STk6O0+7jLCaTFBsrVapk3Q4IkGbNMjYTXM/uloOnnnpKM2bM0D333KNGjRoV2lw3ZsyYS76/devW6t+/vx577LFi35uWA/dCywEAeBaLxaKoqChbH/3bb79da9askdlsNjiZY23fvl0PP/ywdp03zc4dd9yhefPmqVatWg6915IlSzR48GBNnDhR5cuXV//+/dW3b1+PajnIFx8vDR8uTZ8udetmdBq4mt3FQc2aNTVo0CCNGjXKrhvnf8CPiYnRoEGDChw7ceKENmzYoAcffPCy76U4cA8UBwDgeU6ePKmIiAhbq8Frr72mkSNHGpzK8U6fPq0XX3xRU6dOte0LCgrSO++8o549ezrsPj/99JOuvfZahYSEaN68eXrkkUc8tjhA6WZ3t6L09HTdeuutJbr5kCFDNGzYMC1YsKDA/t9//922YAsAAHC8ypUra+HChbaZB19++WVt3LjR4FSO5+fnp7feektr1qxRzZo1JUmpqanq1auXevXqpZMnTzrkPq1bt1ZISIhDrgUYye7i4M4779SGDRtKdPPo6GjFxMTo0UcfZdESAABcrF27dnrhhRckSbm5uerVq5f++ecfg1M5R4cOHZSYmFigtWDJkiVq3Lixvv76awOTAe7F7uJg4sSJiouLU2xsrP744w9lZWVd9CqKAQMGaPLkyYqOjtbnn39ubxwAAGCHsWPHqmXLlpKk/fv3a9iwYQYncp7KlStr8eLFWrx4sW2BtD/++EMdO3bU008/rdOnTxucEDCe3cXBjTfeqN27d2vIkCGqXbu2KlWqdNGrqIYPH67Ro0erW7duWrt2rb2RAABAMZUtW1aLFy9WxYoVJUlz585VfHy8wamcq2fPntqxY4fuuOMO276pU6fqpptu0vbt2w1MBhjP7vlGR48efdEKySXx4osv6p9//lHnzp01duxYh10XAABcXv369TVz5kz169dPkvTYY4+pdevWuvbaa40N5kS1atXSf//7X02bNk0vvPCCsrOztWvXLrVo0UKvvPKKnnnmGa+bvQkoCsMWQdu2bZvq1q1ra9bLN3z4cM2cOVOSLjkbEbMVuRdmKwIAz2exWNSrVy99+OGHkqQ2bdpo3bp1peID8s6dO/Xwww8rMTHRtq9t27ZasGCBateuXezrMVsRPJnd3YpKqkmTJhcVBpI0bdo0DR06tFjdkgAAQMmYTCbNmjXL9mH4+++/12uvvWZwKtdo2LChfv75Zz333HO2XhHfffedGjdurAULFsig71EBQxhWHFzOtGnTlJqaanQMAABKlaCgIMXFxcnHx/rxYOzYsfrpp58MTuUa5cqV0+uvv65169bZulOlp6erb9++ioqK0t9//21wQu8XHy/VqCEtW2Z0ktLN4cVB9+7dNWnSJK1evVpJSUmOvjwAAHCi2267zbbAaV5enh5++GGlp6cbnMp12rVrp8TERPXu3du276OPPlKjRo20evVqA5N5t6QkaeBA6c8/pcces27DGA4vDpYtW6aXXnpJ9957r2rUqKFrrrlG999/v15++WWtWLFCv/32m6NvCQAAHOjll1/WzTffLEk6ePCghg4danAi1woMDNSCBQsUHx+vKlWqSJKOHz+uu+++W0888USRp2tH0Vgs0qBBUkaGdTsjQxo82NhMpZnDByRPnz5dW7Zs0ebNm7Vnzx7boOHzZzYKCAhQkyZN1LRpUzVp0kR9+vQp1j0YkOxeGJAMAN7nt99+U0REhDLOfWJbvHhxgQXESotjx47pkUce0VdffWXb9+yzz2ry5MkXnXu5WRznzp1rmw0KBS1dKvXoUfj+qCjX5yntnDpbUVZWlrZt26bNmzdr06ZNWrt2rY4ePWq9sckki8Vi14d8igP3QnEAAN5p0aJFio6OlmT9Ym/79u2qU6eOsaEMYLFYFBMToxEjRqhq1apKTExU5cqVjY7lFZKSpNBQKS3N2oKQz2SSAgOlvXul4GDj8pVGLp3K1GKxaN26dRo4cKBOnDihPn36aP/+/fryyy+LdR2KA/dCcQAA3is6OlqLFi2SJN1yyy369ttvVaaM3cskebRffvlFJ0+e1C233GJ0FK9gsUiRkVJCglTYRzqzWercWVq+3PXZSjND1jlIT09Xx44dFRwcrM8++6zY76c4cC8UBwDgvdLS0tSkSRP9/vvvkqwzGI0ZM8bYUPAKO3dKjRoV7bzwcOfngZUhU5kGBARo1qxZ+uKLL7x+iXYAADxZYGCgFi9ebFsMbfz48frhhx8MTgVvEB4udelibSEojNksde1KYeBqhq1z0Lx5c4WEhOjdd981KgIAACiCm2++WaNHj5YknT17Vg8//LDS0tIMTgVPZzJJsbFSpUrWP194LCBAmjXLmGylmcOLg+IsXlaxYkVt27bN0REAAICDvfjii7r11lslSYcOHdLjjz9ucCJ4g+Bga4FwYSd3i8W6n8HIrufwMQdms1l169ZVs2bN1KxZMzVv3lzNmze3zROc7/Dhw7r++utVoUIFpaSkFOsejDlwL4w5AIDS4dChQ4qIiLC1GixcuNA2mxFgrwsHJjMQ2VhOKQ7yL3n+fL+1atVSs2bN1LBhQ5lMJi1cuFCHDh1S27ZttW7dumLdg+LAvVAcAEDp8eGHH9rWO6hUqZK2bdumevXqGZwKni5/StPUVKlyZWnPHloNjOLw4iArK0tbt27Vli1bbK9ffvlFubm51hueKxgsFov8/Py0Zs2aYk8JRnHgXigOAKB06du3rxYsWCBJat26tb777juVLVvW4FTwdPHx0vDh0vTpUrduRqcpvVwylWl2drYSExO1ZcsW7d69W6mpqapTp4569+6t6667rtjXozhwLxQHAFC6ZGRkqEmTJjp48KAk6eWXX9b48eMNTgXAEQxZ56CkKA7cC8UBAJQ+GzZs0K233qq8vDz5+Pjom2++UZs2bYyOBaCEDJvKtCTOnj1LYQAAgIFatWqlcePGSbL+vRwdHV2sGQsBuCePLA4AAIDxRo4cqbZt20qyzkI4cOBAeWCHBADnoTgAAAB2MZvNiouLU1BQkCQpPj5e8+fPNzYUgBKhOAAAAHarVauW3n33Xdv20KFDtX//fgMTASgJigMAAFAi3bp1U//+/SVJmZmZ6tWrl3JycgxOBcAeFAcAAKDEpk2bpuuvv16StHHjRo0ZM8bgRMCVxcdLNWpIy5YZncR9eORUpnAvTGUKAJCsRcEtt9yi3NxcmUwmrV27VrfffrvRsYBCnb8qc1CQtHcvqzJLtBwAAAAHadGihV555RVJksViUXR0tFJSUgxOBVzMYpEGDZIyMqzbGRnS4MHGZnIXxW45MJvNzsriMUwmk3Jzc42O4TZoOQAA5Dt79qw6duyodevWSZIiIyO1bNkymUwmg5MB/1q6VOrRo/D9UVGuz+NOil0c+PjQ2MDqzAVRHAAAzvfHH3+ocePGOnnypCTp/fff1//93/8ZnAqwyu9OlJZmbUHIZzJJgYF0L7Lrk77JZCq1LwAAcHk1a9bU+++/b9seNmyY9u7da2AiwOr87kQXfj1usdC9SGLMAQAAcIKuXbtqwIABkqSsrCw9/PDDOnPmjMGpUNrt2iWtXCldqgNIXp60YoX1vNLKruLAYrGU2hcAACiaqVOnKjQ0VJK0efNmvfzyywYnQmkXHi516SJdagit2Sx17Wo9r7Qq9piDcePGOSuLR2H+5n8x5gAAcClbtmxR69atbYuirVmzRh06dDA4FUqzy405CAqS9uwp3WMOWOcAJUZxAAC4nDfffFMjRoyQJIWEhCgxMVFXXXWVwalQmjFb0aUx5gAAADjV008/rY4dO0qSjh07pkcffZSuujBUVFTB7kX53YlKe2EgURwAAAAn8/Hx0fz5822tBR9//LHeffddg1OhNDOZpNhYqVIl63ZAgDRrlrGZ3AXFAQAAcLqQkBB98MEHtu2nnnpKv/zyi4GJUNoFB0uzZ0vVq1t/luZxBudjzAFKjDEHAICiGjx4sGJjYyVJTZo00U8//aRy5coZnApAPloOAACAy0yZMkU33nijJGnbtm168cUXDU4E4HwUBwAAwGX8/f21ZMkS+fr6SpLeeustffXVVwanApCP4gAAALhURESEJk2aZNvu27ev/vrrLwMTAchHcVCKbdiwQb1791bbtm11++23Kzw8XAMHDlRqaqrR0QAAXm748OG66667JEl//vmn+vfvz/SmgBugOCjFli1bJl9fX33zzTf65ptv9MMPP+i///2vBg0aZHQ0AICX8/Hx0bx581StWjVJ0meffaZZzCUJGI7ZikqxvXv3qnLlygo+b+6uyMhIHTx4UFu3bi3ydZitCABgr88++0wPPPCAJMnPz0+bNm1SeHi4wamA0ouWAxc7ffq0XnrpJfn6+qpfv35FOn/8+PEKDQ2Vn5+fatSood69e+u3334rcZbQ0NAChcEPP/yg7777Ts8880yJrw0AQFHcf//9Gjp0qCTr33k9e/bU6dOnDU4FlF4UBy60bt06NW7cWLNmzVJOTs4Vzz916pQ6dOigyZMn65VXXlFKSoq+/PJL7dixQ02bNtX27dsdkmvx4sWqX7++OnfurOnTpys6Otoh1wUAoCjeeOMNNWzYUJK0Y8cOjRw50uBEgHPFx0s1akjLlhmd5GJ0K3KRJUuWaPDgwZo4caLKly+v/v37q2/fvpo3b94l3/P888/rjTfe0MyZMzVkyBDb/gMHDig0NFQNGjTQ9u3bZTabJUmpqan6888/L5vD399f1157baHHfvzxR3Xu3FkjR47U008/XeTfjW5FAICS2rFjh1q0aKHs7GxJ0hdffKF77rnH4FSA4yUlSaGhUmqqFBQk7d3rXqsz03LgInXr1tXu3bv1+OOPy2QyXfH8rKwsxcTEyNfXV3379i1wrH79+urYsaN27dqlL774wrb/ww8/1I033njZV58+fS55z5tvvlnDhw/Xiy++yIxFAACXatSokSZPnmzb7tevn06cOGFgIsDxLBZp0CApI8O6nZEhDR5sbKYLURy4SOvWrRUSElLk89euXavMzExFRESoYsWKFx2/7bbbJEkJCQm2fYMGDZLFYrns65tvvrGdn//tzPnCwsKUnZ2tffv2FeO3AwCg5IYOHap7771XkpSUlMT0pvA68fHSypVSXp51Oy9PWrHCut9dUBy4qfzxBHXr1i30eL169QqcZ4/Q0FAlJycX2Hf06FFJUpUqVey+LgAA9jCZTJo7d66uvvpqSdauRTNnzjQ4FeAYSUnWVoMLO5CYTNLAgdbj7oDiwE3ljx241If0ypUrFzjPXpMmTdLZs2clWQuDt99+W3fccYeuu+66S74nOztb6enpBV4AADhCcHBwgfF4I0aM0I4dO4wLBDjA+d2JLmwMs1jcq3sRxYGbysrKkiT5+voWerxcuXIFzrPHpEmTtGPHDrVs2VJt27bVfffdp27dumnlypWXfd9rr72mwMBA26tWrVp2ZwAA4EJ33323hg8fLsn6hVTPnj116tQpg1MB9tu1q2B3ogvldy/atcu1uQpDceCm/P39JUlnzpwp9Hj+eIH88+zRo0cPrV69Wps2bdJ3332nbdu26bXXXrvijEMvvPCC0tLSbK8jR47YnQEAgMJMmjRJjRs3liTt2rVLzz33nMGJAPuFh0tdukjnJpi8iNksde1qPc9oFAduqnr16pKklJSUQo+fPHlSklSjRg2XZcpXrlw5BQQEFHgBAOBIfn5+WrJkifz8/CRJM2fO1GeffWZwKsA+JpMUGytVqlT4mIOAAGnWLGOyXYjiwE01adJEki65EvLBgwclyfatCgAA3iYsLExTpkyxbT/yyCMlHmsHGCU42FogFDbmIDbWfdY6oDhwU+3bt1eFChWUmJiozMzMi47/8MMPkqQHH3zQ1dEAAHCZwYMH64EHHpAkJScnq2/fvraJNABPExVVsHtRfneiqChjc52P4sBN+fv7a+jQocrOztb8+fMLHDt48KDWrFmj8PBw23zQAAB4I5PJpDlz5ti623711VeaNm2awakA+5zfvUhyr+5E+SgO3NiYMWN08803a+TIkVq+fLlOnTqlxMRERUZGyt/fX4sWLZL5UiNbAADwEtWqVdOCBQts2yNHjtS2bduMCwSUQHCwNHu2VL269ae7dCfKZ7Kw9KDLmC4cgXKeuXPnql+/fhftP336tF5//XUtWrRIhw8fVmBgoO68806NHz/ethCa0dLT0xUYGKi0tDQGJwMAnObZZ5+1jUG48cYbtWnTphLN2gfgYhQHKDGKAwCAK2RnZ6t169a2VoNBgwZplrv1yQA8HN2KYLeYmBiFhYWpRYsWRkcBAJQC5cqV05IlS1S+fHlJUmxsrD755BODUwHehZYDlBgtBwAAV3r33Xc1cOBASdJVV12lxMREhYSEGJwK8A60HAAAAI8yYMAAdenSRZL0999/q0+fPkxvCjgIxQEAAPAoJpNJ7733nq214Ouvv9Zbb71lcCrAO1AcAAAAj3PVVVdpwYIFtpkAX3zxRW3ZssXgVIDnozgAAAAeqUOHDhoxYoQkKScnRz179lRmZqbBqQDPRnEAAAA81oQJE9S8eXNJ0r59+/TUU08ZnAjwbBQHAADAY/n6+mrx4sW2xdDee+89LV++3OBUgOeiOAAAAB7thhtu0PTp023bAwYM0B9//GFgIsBzURwAAACP179/f0VGRkqSTp48qT59+igvL8/gVIDnoTiA3VghGQDgLkwmk959913VrFlTkrRu3TpNnjzZ4FSA52GFZJQYKyQDANzFt99+q/bt28tisahMmTJav349X2IBxUDLAQAA8Brt2rXTyJEjJUm5ubnq1auX/vnnH4NTAZ6D4gAAAHiVcePG2VoL9u/fr2HDhhmcCPAcFAcAAMCrlC1bVosXL1bFihUlSXPnztWyZcsMTgV4BooDAADgda677jrNmDHDtv3YY4/p8OHDBiYCPAPFAQAA8Ep9+/ZV9+7dJUmpqamKjo5melPgCigOAACAVzKZTIqNjdW1114rSfr+++81adIkg1MB7o3iAAAAeK2goCDFxcXJx8f6kWfMmDHasGGDwakA90VxAAAAvFqbNm300ksvSZLy8vLUq1cvpaenG5wKcE8UBwAAwOuNHj1arVu3liQdPHhQTzzxhMGJAPdEcQAAALxemTJltGjRIlWqVEmStGDBAi1ZssTgVID7oTiA3WJiYhQWFsay9AAAj1CvXj298847tu1Bgwbp999/Ny4Q4IZMFovFYnQIeLb09HQFBgYqLS1NAQEBRscBAOCyHn74YS1evFiSdOutt+qbb75RmTJlDE4FuAdaDgAAQKnyzjvvqE6dOpKkH374QRMnTjQ2EOBGKA4AAECpEhgYqEWLFtmmNx03bpzWr19vcCrAPVAcAACAUueWW27R6NGjJUlnz57Vww8/rLS0NLuvt3btWv3f//2fbrjhBvn5+cnf319hYWEaMWKE/vrrL0fFBpyO4gAAAJRKL730km699VZJ0u+//64hQ4bYdZ2YmBh16NBBmzZtUmxsrP766y/9+uuv6tevn6ZNm6ZGjRpp3759jowOOA3FAQAAKJXKlCmjuLg422QaixYtUlxcXLGvc+rUKfn6+uqzzz7THXfcoUqVKumaa67Rc889p6efflonTpzQ8OHDHR0fcAqKAwAAUGrVqVNHsbGxtu3HH39cBw8eLNY1qlevrp49e6pWrVoXHXvwwQclSWvWrFFeXl7JwgIuwFSmKDGmMgUAeLo+ffpo4cKFkqTWrVvr+++/d8j0prt27VLDhg1VtmxZZWVlMWUq3B4tBwAAoNSbOXOm6tWrJ0n66aefNH78eIdcd8+ePZKk2267jcIAHoHiAAAAlHoBAQFavHixzGazJOnVV1/V999/X+LrLliwQCaTyTYzEuDuKA4AAAAktWrVSmPHjpVknd40Ojpaqampdl/vyy+/VEJCgp5++mndfvvtDskIOBtjDlBijDkAAHiLvLw8tW/f3tZq0L17dy1ZskQmk6lY1/n111912223qU2bNlq6dKmtRQJwd7QcwG4xMTEKCwtTixYtjI4CAIBDmM1mxcXFKSgoSJK0dOlSLViwoFjXOHTokO68807ddtttBboqAZ6AlgOUGC0HAABvEx8fr+7du0uSKlasqK1bt+q666674vv279+vDh06qG3btpo3bx6FATwOLQcAAAAXiIqK0iOPPCJJ+ueff9SrVy/l5ORc9j27d+9W27Zt1alTJ82fP79AYfDmm2/qyJEjTs0MOALFAQAAQCGmT59uay3YuHGjbbByYbZv367bb79dkZGRevfdd+XjU/Aj1ogRI3TgwAFnxgUcggl3AQAAClGxYkUtXrxYt9xyi3Jzc/Xaa6+pU6dOateuXYHzNm3apE6dOik7O1t//fWXevbsaVBioOQoDgAAAC6hRYsWmjBhgl544QVZLBZFR0dr+/btqlKliu2cuLg4nTx5UpJ1ADPgyRiQjBJjQDIAwJvl5eWpY8eO+uabbyRJzz77rCZPnmxsKMBJGHMAAABwGWazWQsXLlSVKlX0xBNPaPz48UZHApyGbkUAAABXULNmTe3Zs0fVqlUzOgrgVLQcAAAAFAGFAUoDigMAAAAAkigOAAAAAJxDcQAAAABAEsUBAAAAgHMoDgAAAABIojhACcTExCgsLEwtWrQwOgoAAAAcgBWSUWKskAwAAOAdaDkAAAAAIIniAAAAAMA5FAcAAAAAJFEcAAAAADiH4gAAAACAJIoDAAAAAOdQHAAAAACQRHEAAAAA4ByKAwAAAACSKA4AAAAAnENxAAAAAECSVMboAPB8FotFkpSenm5wEgAAAFxKpUqVZDKZLnsOxQFKLCMjQ5JUq1Ytg5MAAADgUtLS0hQQEHDZc0yW/K99ATudPXtWx44d0x133KFNmzbZdY0WLVpo48aNRTo3PT1dtWrV0pEjR674HziKrjj/DtyZO/0erszirHs58rolvVZJ3l/c9/KccTx3+n+zJNzp9+AZ4/hr2XsNT3nG0HIAl/Dx8VHNmjVVpkwZu/8DN5vNxX5vQEAAf2k7kD3/DtyRO/0erszirHs58rolvVZJ3m/ve3nOOI47/b9ZEu70e/CMcfy17L2GNz1jGJAMhxkyZIgh74VjeMu/A3f6PVyZxVn3cuR1S3otnjGezVv+HbjT78EzxvHXsvca7vTfRUnRrQgeJz09XYGBgUXqNwcA9uA5A8CZ3PkZQ8sBPE65cuU0ZswYlStXzugoALwUzxkAzuTOzxhaDgAAAABIouUAAAAAwDkUBwAAAAAkURyglMjLy9PkyZNVvnx5ffPNN0bHAeAlvvjiC3Xu3Fnt27fXbbfdplatWmnp0qVGxwLgJTZs2KDevXurbdu2uv322xUeHq6BAwcqNTXVafdknQN4vT///FNRUVEKCwvT6dOnjY4DwIv06dNH48aNs01j+Omnn+rBBx+Ur6+vunTpYnA6AJ5u2bJl8vX11TfffCMfHx+lpqaqWbNmSktL04cffuiUezIgGV5v//79ysrKUkBAgOrWrat169bp9ttvNzoWAC9w3333KSEhQWaz2bYvPDxcDRo00PLlyw1MBsAb7N27V5UrV1ZwcLBtX2RkpA4ePKitW7c65Z50K4JbOH36tF566SX5+vqqX79+RTp//PjxCg0NlZ+fn2rUqKHevXvrt99+u+jc6667To0bN3ZCagCewJnPl88//7xAYSBJfn5+OnPmjKPiA3BzznzGhIaGFigMfvjhB3333Xd65plnHPkrFEBxAMOtW7dOjRs31qxZs5STk3PF80+dOqUOHTpo8uTJeuWVV5SSkqIvv/xSO3bsUNOmTbV9+3YXpAbgCVz9fElJSdGuXbvUq1cvR/0KANyYq54xixcvVv369dW5c2dNnz5d0dHRjv5V/mUBDLR48WJLYGCgJSYmxvLBBx9YJFn69u172fc899xzFkmWmTNnFti/f/9+i9lstoSHh1tyc3Mvet9vv/1mkWRZt26dA38DAO7Klc+XfMOHD7fcfffdlrNnzzriVwDgxox4xqxfv95SrVo1y5QpUxzxKxSKlgMYqm7dutq9e7cef/xxmUymK56flZWlmJgY+fr6qm/fvgWO1a9fXx07dtSuXbv0xRdfOCsyAA/h6udLXFycvv32W3344YdFuh8Az2bEZ5ibb75Zw4cP14svvui0GYsoDmCo1q1bKyQkpMjnr127VpmZmYqIiFDFihUvOn7bbbdJkhISEhyWEYBncuXzZfny5Xrrrbf01VdfKTAw0P7QADyGK54x2dnZF50XFham7Oxs7du3z47UV0ZxAI+S3xevbt26hR6vV69egfMAoKjsfb6sWLFCr776qlavXq1q1apJkl5//XUnJgXgiex5xoSGhio5ObnAeUePHpUkValSxRkxKQ7gWf78809Jl/4fonLlygXOA4Cisuf5smLFCg0ZMkSvvvqqDh06pE2bNmnTpk2aNWuW8wMD8Cj2foaZNGmSzp49K8laGLz99tu64447dN111zklJ4ugwaNkZWVJknx9fQs9Xq5cuQLn5bvnnnt08uRJSdKTTz6poKAgLVu2zPYtHwDY83zp0aOHcnJydO+99xY4t3bt2k5KCcBT2fOMmTRpkubOnauWLVvK399f6enp6tatm1544QWn5aQ4gEfx9/eXpEvOIZ7fNy//vHyrVq1ybjAAHs+e5wvrGQAoKnueMT169FCPHj2cH+48dCuCR6levbok61zihclvHahRo4bLMgHwDjxfADiTpzxjKA7gUZo0aSJJha4iKEkHDx6UJFZEBlBsPF8AOJOnPGMoDuBR2rdvrwoVKigxMVGZmZkXHf/hhx8kSQ8++KCrowHwcDxfADiTpzxjKA7gUfz9/TV06FBlZ2dr/vz5BY4dPHhQa9asUXh4+EWDAwHgSni+AHAmT3nGUBzA44wZM0Y333yzRo4cqeXLl+vUqVNKTExUZGSk/P39tWjRIpnNZqNjAvBAPF8AOJMnPGNMFovFYmgClHqXW3J87ty56tev30X7T58+rddff12LFi3S4cOHFRgYqDvvvFPjx4+3LSICADxfADiTNz5jKA4AAAAASKJbEQAAAIBzKA4AAAAASKI4AAAAAHAOxQEAAAAASRQHAAAAAM6hOAAAAAAgieIAAAAAwDkUBwAAAAAkURwAAAAAOIfiAAAAAIAkigMAAAAA51AcAAAAAJBEcQAA8CD79++XyWTSnXfeqczMTI0bN0433nij/P39Vb9+fb3xxhuyWCxGxwQAj1XG6AAAABTV9u3bJUlVqlRRRESEsrOzdfPNN6tWrVr6+uuv9fzzz6tq1arq37+/wUkBwDPRcgAA8Bj5xUF8fLy6d++u3377TfHx8frqq680Z84cSdKKFSuMjAgAHo3iAADgMbZt2yZJioyM1KuvvqoyZf5tAL/77rslSUlJSS7JMmrUKNWpU8cl9wIAV6E4AAB4jPyWg9GjR190LDk5WZIUEhLi0kwA4E0oDgAAHiE1NVWHDx9W3bp11bhx44uOJyYmSpIaNmzo6mgA4DUoDgAAHiG/S1HLli0LPb5161ZJUpMmTWz7+vXrp9tuu02ffPKJQkNDValSJd177706fvx4gfdOmTJF9erVU7ly5RQWFqbFixcXOH727Fk999xzqly5sq666iq98MILOnv27EUZPvroIzVt2lR+fn6qWbOmJkyYUGD2pB07dqhTp04KCgpSpUqV1Lhx44vuBQBGojgAAHiE/C5FTZs2LfR4YcWBZJ3+dMKECXrttdf0/vvva/PmzRo6dKjt+MyZMzVixAj16tVLCQkJatOmjR5++GGtWrXKds6bb76pqVOn6rnnnlNcXJx27NihefPmFbjPkiVL1L17d91+++1KSEjQs88+q0mTJumtt96ynfPggw/K19dXS5Ys0cqVK/V///d/Sk1NLcE/FQBwLJOFCaEBAB6gf//+mjt3rlavXq1OnTpddPyqq65STk6O0tLSZDKZJFlbDhYvXqwDBw6oVq1akqytBCNHjlR2drYsFotq1qypyMhIzZw503ate+65RykpKdqwYYNyc3N1zTXXKDo6WlOmTJEknTlzRrVr11a5cuX0+++/y2KxqHbt2nrwwQcLXOfNN9/U66+/rmPHjiktLU3VqlVTYmKiGjVq5Mx/VABgN1oOAAAeIb9bUWEtB4cPH1ZKSooaN25sKwzy3XDDDbbCQJJuvPFG5ebmKikpSX/88Yf+/PNPRUZGFnhPt27dtHXrVuXl5enIkSNKSkrSgw8+aDvu6+uru+66y7a9b98+HTlyRP/5z3+Um5tre7Vv317Jyck6fPiwqlSpolq1amnw4MFatmyZ/v77b0f8YwEAh6I4AAC4vdzcXO3evVvXXHONqlWrdtHxS3UpkqSgoKAC276+vpKk06dP28YeBAcHFzjn6quvVk5OjpKTk3XixAlJuui+52/nz5TUvn17lS1b1va66aabJElHjhyRj4+PVq9erWrVqqlv3766+uqr1alTJ+3bt6+o/xgAwOlYIRkA4Pb27Nmj7OzsQj/8S5cvDi6nRo0akqxrI4SHh9v2nzhxQmXLllXVqlWVlZUlSfrrr78KvPf87SpVqkiS5s+fr7CwsIvuExoaKsnaarFy5UplZ2dr3bp1euaZZ9SzZ09t3ry5WLkBwFkoDgAAbq9hw4a63BC5sWPHauzYscW+bs2aNVW9enWtWLFC7du3t+3/6KOP1KxZM5nNZtWqVUvBwcFKSEhQu3btJFnHHKxevVrlypWTZP3wHxISoiNHjqhPnz5XvG+5cuV0991369dff9Xzzz9f7NwA4CwUBwCAUstsNuull17S8OHDVblyZd12221auXKlVq1aZZutqEyZMnrqqaf08ssvq2rVqmrSpIliYmIKjG3w8fHR5MmT9cgjj+jkyZPq2LGjfHx8tG/fPv33v//VJ598osTERI0YMULdu3dXvXr1lJSUpBkzZuiOO+4w6tcHgItQHAAASrWhQ4fq9OnTiomJ0aRJk1S/fn3FxcXp7rvvtp0zYsQIJScn6/XXX5ePj48effTRi9Yo6NWrlwICAjRx4kS988478vX11fXXX6+oqChJ1nEMVatW1YQJE3T8+HFVqVJF9957r15//XWX/84AcClMZQoAAABAErMVAQAAADiH4gAAAACAJIoDAAAAAOdQHAAAAACQRHEAAAAA4ByKAwAAAACSKA4AAAAAnENxAAAAAEASxQEAAACAcygOAAAAAEiiOAAAAABwDsUBAAAAAEnS/wMIQPLeDhtT1wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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(r'$n_\\mathrm{nodes}$')\n",
    "ax.set_ylabel(r'$\\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.12.7"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}