Exercise_02a_Scale-Selective_Damping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 02.a: Scale-selective damping\n",
    "\n",
    "**Due:** _25 November 2020_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Names:** David, Lukas, Mariana, Linus (for now:D)  (Up to 3 people can submit together.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the harmonic and biharmonic diffusion equations\n",
    "\n",
    "$$\n",
    "  \\frac{\\partial u}{\\partial t} = A_2 \\frac{\\partial^2 u}{\\partial x^2} \\quad\\text{(harmonic)}\n",
    "  \\qquad\\qquad\n",
    "  \\frac{\\partial u}{\\partial t} = -A_4 \\frac{\\partial^4 u}{\\partial x^4} \\quad\\text{(biharmonic)}\n",
    "$$\n",
    "\n",
    "where $A_2$ and $A_4$ are the positive definite harmonic and biharmonic diffusion coefficients,\n",
    "respectively.\n",
    "\n",
    "1. Discretize the harmonic equation in the x direction. Derive the damping time scale as a function of wave number.\n",
    "\n",
    "1. Discretize the biharmonic equation in the x direction. Derive the damping time scale as a function of wave number.\n",
    "\n",
    "1. Plot both damping time scales as a function of wave number. How do they differ for long and short waves and why?\n",
    "\n",
    "1. Derive the damping time scale as a function of wave number for the analytic solution of both equations and compare them to the time scales implied by the finite difference equations.\n",
    "\n",
    "1. Show that in the limit $\\Delta x \\rightarrow 0$, the damping time scales of the discretized and of the analytical solution are the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint:** Since the analytical solution of both diffusion equations has the form\n",
    "\n",
    "$$u(x, t) = \\hat{u}(t) {\\rm e}^{ikx}$$\n",
    "\n",
    "where $\\hat{u}$ is independent of $x$, $i=\\sqrt{-1}$ and $k$ is the wave number.\n",
    "After discretizing the equations in the x direction, one can substitute $u(j\\Delta x, t)$\n",
    "into the discrete equation to get the damping coefficient as a function of wave number $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.)\n",
    "To discretize the harmonic equation in space, we make use of the derivation of the discrete second order derivative in the script. (eq. 1.12)\n",
    "With this, the harmonic equation takes the form\n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t} u_{i} = A_{2}\\frac{\\partial^{2} u_{i}}{\\partial x^{2}} = A_{2}\\frac{u_{i+1}-2 u_{i}+u_{i-1}}{\\Delta x^2}\n",
    "$$\n",
    "\n",
    "Substituting the analytical solution for each gridpoint \n",
    "\n",
    "$$\n",
    "u_{j}(x,t) = \\hat{u}(t) e^{i k (j \\Delta x)}\n",
    "$$\n",
    "\n",
    "as given in the hint, yields\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial}{\\partial t} \\hat{u}(t) e^{i k j \\Delta x} &= \\hat{u}(t)\\frac{A_{2}}{\\Delta x^2}\\left( e^{i k (j+1) \\Delta x}-2  e^{i k j \\Delta x} +  e^{i k (j-1) \\Delta x} \\right) \\notag \\\\\n",
    "&= \\frac{A_{2}}{\\Delta x^2}\\left( e^{i k j \\Delta x} e^{i k \\Delta x}-2  e^{i k j \\Delta x} +  e^{i k j \\Delta x} e^{i k \\Delta x} \\right)\\hat{u}(t).\n",
    "\\end{align}\n",
    "\n",
    "Deviding by $e^{i k j \\Delta x}$ then leaves\n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t} \\hat{u}(t) = \\frac{A_{2}}{\\Delta x^2} \\left( e^{i k  \\Delta x} + e^{- i k \\Delta x} -2 \\right) \\hat{u}(t) \n",
    "$$\n",
    "\n",
    "Using $ e^{i k  \\Delta x} + e^{- i k \\Delta x} = 2 \\cos(k \\Delta x)$ simplifies above to \n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t} \\hat{u}(t) = \\frac{2 A_{2}}{\\Delta x^2} \\left( \\cos(k \\Delta x) -1 \\right) \\hat{u}(t).\n",
    "$$\n",
    "\n",
    "While this simple ODE is well understood, the solution for $\\hat{u}(t)$ is of form\n",
    "\n",
    "$$\n",
    "\\hat{u}(t) \\propto e^{\\frac{2 A_{2}}{\\Delta x^2} \\left( \\cos(k \\Delta x) -1 \\right) t},\n",
    "$$\n",
    "\n",
    "where the damping time scale can be identified as the negative inverse of the exponent (since $\\hat{u}(t)\\propto e^{-\\frac{t}{\\tau}}$):\n",
    "\n",
    "$$\n",
    "\\tau_{h} = \\frac{\\Delta x^2}{ 2 A_{2}\\left(1- \\cos(k \\Delta x) \\right)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.)\n",
    "\n",
    "To discretizing the biharmonic equation in space, we make use of the derivation of the discrete second order derivative as in 1.), but applying it twice. Namely\n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial t} u_{j} = - A_{4}\\frac{\\partial^{4} u_{j}}{\\partial x^{4}} = -A_{4}\\frac{\\partial^2 f_{j}}{\\partial x^2}\n",
    "$$\n",
    "\n",
    "with \n",
    "\n",
    "$$\n",
    "f_{j}=\\frac{\\partial^2 u_{j}}{\\partial x^2}\n",
    "$$\n",
    "\n",
    "Sparing algebraic details this yields\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial}{\\partial t} \\hat{u}(t) e^{i k j \\Delta x} &= - A_{4} \\left( \\frac{u_{j+2}-4 u_{j+1}+6 u_{j}-4 u_{j-1}=u_{j-2}}{\\Delta x^4} \\right) \\hat{u}(t) \\notag \\\\\n",
    "&= - \\frac{A_{4}}{\\Delta x^4} e^{ikj \\Delta x}\\left( e^{2ik \\Delta x}-4 e^{ik \\Delta x}+6-4 e^{-ik \\Delta x}+e^{-2ik \\Delta x} \\right) \\hat{u}(t) \\notag \\\\\n",
    "&= - \\frac{A_{4}}{\\Delta x^4} e^{ikj \\Delta x}\\left( 2 \\cos(2k \\Delta x) -8\\cos(2k \\Delta x)+6 \\right) \\hat{u}(t)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "The solution for $\\hat{u}(t)$ is of the same form as the previous case, therefore the damping time scale is given by:\n",
    "\n",
    "$$\n",
    "\\tau_{b} = \\frac{\\Delta x^4}{A_{4} 2\\left( \\cos(2k \\Delta x) -4\\cos(2k \\Delta x)+3 \\right) }\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N:       21\n",
      "L:       10\n",
      "A_2:     1.0\n",
      "A_4:     1.0\n",
      "dx:      0.5\n",
      "x:       [ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5\n",
      "  7.   7.5  8.   8.5  9.   9.5 10. ]\n"
     ]
    }
   ],
   "source": [
    "#Initializing diffusion object, grid and diffusion parameters are set to default values\n",
    "D = Diffusion()\n",
    "D.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcQAAAEICAYAAADBdYTQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABNZUlEQVR4nO3dd3QUVfvA8e9NAiGhdzL0jnTpIiBILwIiAoqCgmIHVsfuC1heFFldfX8qSpOiglhARDpKk957r2FCkdBLIMn9/bETTAIkm7K7SXg+5+zZ3btzZ5+R4z65c5vSWiOEEELc6QL8HYAQQgiREUhCFEIIIZCEKIQQQgCSEIUQQghAEqIQQggBQJC/A0gPAQEBOiQkxN9hCCFEpnL58mWttZaGkS1LJMSQkBAuXbrk7zCEECJTUUpd8XcMGYn8ZSCEEEIgCVEIIYQAJCEKIYQQgCREIYQQApCEKIQQQgA+HGVqGI5AYB1wzLJcnQzDUQD4ESgDHAJ6WJbrjH3sm0B/IAYYaFmueb6KUwghhJtTqfFAJ+CkqXX1eOUvAS8C0cAfptav2eUJfrtNrefZ5XWBCUAIMBsYZGqtnUoFA5OAusBpoKep9SG7Tl/gHfsrPzC1nujdq/VtC3EQsDPe+zeARZblqggsst9jGI6qQC+gGtAO+MpOpkIIIXxrAu7f4RucSrUAugA1Ta2rAU67/KbfbqdScb/do4ABQEX7EXfO/sAZU+sKgAsYYZ+rADAUaAg0AIY6lcrvnUv8l08SomE4SgAdgbHxirsAcRl/ItA1XvlUy3JFWZbrILAP93+QdHdk627+27grx/cd9sbphRAiUzO1XgpEJip+DvjI1DrKPuakXd4FmGpqHWVqfeO326lUGJDH1Hql6d5vcBIJf+/j8sDPQEunUgpoCywwtY40tT4DLCBRYvYGX7UQPwNeA2LjlRW1LFcEgP1cxC4vDhyNd1y4XZaAUmqAUmqdUmpddHR0qoI6ti+c4JW/Meb5d1NVXwghMrmguN9R+zHAgzqVgKZOpVY7lVriVKq+XX673+7i9uvE5QnqmFpHA+eAgkmcy6u8nhANw9EJOGlZrvUeVlG3KLtpF2Ot9WitdT2tdb2goNR1hd7zYEuii5Xn4uJZ7NxppeocQgiRiUXH/Y7aj9Ee1AkC8gONgFeBaXar7na/3Un9pqemjtf4ooV4L9DZMByHgKnA/Ybh+A44YRiOMAD7Oa7ZHQ6UjFe/BOC1bNVk8HMUvn6K4S/+D3drXgghRBLCgV9NrbWp9Rrcd/4Kcfvf7nD7deJy4tdxKhUE5MV9i9aneSCO1xOiZbnetCxXCctylcHd4fqnZbkeA2YCfe3D+gK/2a9nAr0MwxFsGI6yuDtg13grvrpPPQGBQUSvns/MmZu89TVCCJFVzADuB3AqVQnIDvyD/dvtVCrYqdSN325T6wjgglOpRnZLsg8Jf+/j8kB34E+7n3Ee0MapVH57ME0bu8yr/DkP8SOgtWE49gKt7fdYlms7MA3YAcwFXrAsV4y3gggpWJAKnR+gxrU9vP/ur1y6FOWtrxJCiEzFqdQUYCVQ2alUuFOp/sB4oJxTqW247/r1tVuLN/12m1rH/XY/h3tQ5T5gPzDHLh8HFHQqtQ94GXu2gal1JPA+sNZ+vGeXeZXKCrcJc+bMqdOy28X+339neufOzMjfmQ7m07z1Vqd0jE4IITImpdRlrXVOf8eRUchKNUCZdu0IKVyYNgUi+OabxezffzL5SkIIIbIUSYhAYLZsVH3sMXId2USeoOsMHTpDBtgIIcQdRhKirVrfvsRev84zDeDPP3cyf/52f4ckhBDCh6QPMZ5Jd98NgUF8m7s7Fy9GsWTJ64SGBqdDhEIIkfFIH2JC0kKMp1rfvpxcv44hz9Ti2LEzOJ2yprgQQtwpJCHGc1fv3gRkywar59G7dyPGjFnC9u3H/B2WEEIIH5CEGE9o4cJU7NaN7RMn8trLrcifP5TXX/+JmJjY5CsLIYTI1CQhJlLz6ae5euYMp/6cy7BhXdmw4TCTJ6/wd1hCCCG8TBJiIqVatCBf+fJsGT2aBx+sQ9Omlfjwwz84ceKcv0MTQgjhRZIQE1EBAdR4+mnCly7lzJ49fPhhd65di2bIkBn+Dk0IIYQXSUK8hepPPEFAUBBbxoyhXLnCDBzYmt9/38Sff+70d2hCCCG8ROYh3sZv3bsTvngxzxw7RgyBtG49kqtXr/Pnn6+RK1eOdP0uIYTwB5mHmJC0EG+j1oABXDl9mn3TpxMcHMQnn/Ti2LGzDB/+h79DE0II4QWSEG+jdKtW5ClThs2j3RtI169flv79mzJhwnJWr97v5+iEEEKkN0mIt6ECAqj51FMc/esvIvfsAeCNNzpQqlQBXnnlR65cuebnCIUQQqQnSYhJqNG/PwHZsrHpq68ACA0NZuTInhw4cAqnc66foxNCCJGeJCEmIWexYlTq3p1t337LtYsXAWjatBK9ezfim28Ws2nTET9HKIQQIr3IKNNkHFuxgin33kurUaOo/eyzAJw/f4UWLUaQN28oc+e+TPbsQV75biGE8KbkRpk6lRoPdAJOmlpXT/SZCYwECpta/2OXvQn0B2KAgabW8+zyusAEIASYDQwytdZOpYKBSUBd4DTQ09T6kF2nL/CO/XUfmFpPTJeLToK0EJNh3HMPRevUYeMXX9zYNDhPnhA++uhhdu2K4PPPF/g5QiGE8JoJQLvEhU6lSgKtgSPxyqoCvYBqdp2vnEoF2h+PAgYAFe1H3Dn7A2dMrSsALmCEfa4CwFCgIdAAGOpUKn86X9tNJCEmQylF7Rdf5PT27RxdvPhGeevW1Xjoobr83/8tZMuWo/4LUAghvMTUeikQeYuPXMBrQPxbjF2AqabWUabWB4F9QAOnUmFAHlPrlaa7VTEJ6BqvTlzL72egpVMpBbQFFphaR5panwEWcIvEnN4kIXqgSq9ehBQsyMYvvkhQ/v773ShUKDcvvfS9jDoVQtwRnEp1Bo6ZWm9O9FFxIH7rINwuK26/TlyeoI6pdTRwDiiYxLm8ShKiB7KFhFDjqafYN2MG54/8O5AmX75QPv20F3v3nmDEiNl+jFAIIVIlSCm1Lt5jQFIHO5UKBd4GhtziY3WLMp1EeWrreI0kRA/Veu45ADZ//XWC8ubNq9C3772MHr2Ev//e64/QhBAitaK11vXiPUYnc3x5oCyw2anUIaAEsMGpVDHcrbiS8Y4tAVh2eYlblBO/jlOpICAv7lu0tzuXV0lC9FDe0qUp37kzW8aMIfrq1QSf/ec/D1CuXGEGD57ChQtXb3MGIYTI3Eytt5paFzG1LmNqXQZ34qpjan0cmAn0cioV7FSqLO7BM2tMrSOAC06lGtn9g32A3+xTzgT62q+7A3/a/YzzgDZOpfLbg2na2GVe5fWEaBiOHIbhWGMYjs2G4dhuGI537fJhhuE4ZhiOTfajQ7w6bxqGY59hOHYbhqOtt2P0VJ2BA7nyzz/s/P77BOWhocF8/vmjREScZciQ6X6KTggh0pdTqSnASqCyU6lwp1L9b3esqfV2YBqwA5gLvGBqHWN//BwwFvdAm/3AHLt8HFDQqdQ+4GXgDftckcD7wFr78Z5d5lVen4doGA4F5LQs10XDcGQDlgODcI8YumhZLmei46sCU3APtTWAhUAly3LFcBvenIcYn9aayXXqEB0VxZPbtqECEv49MWLEbD7/fAHjx/ejXbsaXo9HCCHSQna7SMjrLUTLcmnLcl2032azH0ll4S7AVMtyRVmW68bQXS+H6RGlFPVeeYXInTs5OPfmpdscjjZUq1acV1+dxj//XPBDhEIIIVLLJ32IhuEINAzHJuAksMCyXKvtj140DMcWw3CMNwxH3KRLj4bbKqUGxI2Mio6O9mb4CVTu2ZNcxYuz7pNPbvose/Yg/u//enPx4lUcjilkhVWAhBDiTuGThGhZrhjLctXGPVKogWE4quNeuaA8UBuIAOIyjEfDbbXWo+NGRgUF+W7ptMBs2agzaBBH/vyTExs33vR5lSph/Oc/nVm0aCfjxi31WVxCCCHSxqejTC3LdRZYDLSzLNcJO1HGAmP497aoX4bbpkStAQPInjv3LVuJAE8+2YTWravxwQe/s23bMR9HJ4QQIjV8Mcq0sGE48tmvQ4BWwC7DcITFO+xBYJv9eibQyzAcwYbhuDF019txpkRw3rzUeOopdk2dyvmjNy/bppTi0097kT9/Tp5/fhKXL0f5IUohhBAp4YsWYhjwl2E4tuAePrvAslyzgI8Nw7HVLm8BOAAsy3XT0N2kRpj6S51BgwDY8L//3fLzggVz8b//9Wb//lMMGTLDh5EJIYRIDdn+KQ1mPfIIB2bP5pkjRwjOm/eWx/z3v7/z5Zd/Mnp0Xzp1qu3bAIUQIgky7SIhWakmDeq/+irXzp9n06hRtz3mtdc6ULt2KV59dRrh4Wd8GJ0QQoiUkISYBkXr1KFM27as+/RTrl++fMtjsmUL5KuvHicmJpYXXpjM9esZ7u6vEEIIJCGmWaO33+bKqVNsHTv2tseUKVOIESN6sHbtQT766A8fRieEEMJTkhDTqETTppRo2pS1I0cSc+32eyI++GAd+vRpzKhRfzFv3rbbHieEEMI/JCGmg4Zvv82F8HB2TJ6c5HHDhnWlRo0SDB78A0eOnPZRdEIIITwho0zTgdaa7+rXJ+rcOfrt3ElAEivnHD78D23bfkLZsoWZMWMgwcG+W2VHCCHik1GmCUkLMR0opWj41luc3beP3T/9lOSxpUsXwuV6lM2bjzJs2AzfBCiEECJZkhDTScWuXSlYtSqrhw9Hx8YmeWz79jV45pnmTJz4NzNmbPBRhEIIIZIiCTGdqIAAGr71Fv9s28be6clvEvzWW52oV68Mr746jb17T/ggQiGEEEmRPsR0FBsTw4Tq1VGBgTyxZctNGwgnZllnadv2E/LlC+WPPwaTJ0+IjyIVQgjpQ0xMWojpKCAwkMbDhnF6+3Z2T5uW7PGGkY/Ro/ty6NA/vPTS98Qmc6tVCCGE90gLMZ3p2Fgm1qpF7PXrPLFtW5IjTuOMHbuUIUOm88orbXnllXY+iFIIIZJvITqVGg90Ak6aWle3y0YCDwDXgP3Ak6bWZ+3P3gT6AzHAQFPreXZ5XWACEALMBgaZWmunUsHAJKAucBroaWp9yK7TF3jHDuUDU+uJ6XfltyYtxHSmAgJo/O67RO7ezc4pUzyq079/Ux5+uB6ffDJPJu0LITKSCUDiv9IXANVNrWsCe4A3AZxKVQV6AdXsOl85lQq064wCBuDezq9ivHP2B86YWlcAXMAI+1wFgKFAQ9x75Q51KpXfC9eXgCREL6jYtStFatdm5bvvEnP9erLHK6X46KOHqVmzBC+99J0MshFCZAim1kuByERl802to+23q3Bv4g7QBZhqah1lan0Q2Ac0cCoVBuQxtV5pum9JTgK6xqsT1/L7GWjpVEoBbYEFptaRptZncCdhr98+k4ToBSoggHvfe4+z+/cnu3pNnJCQ7Iwb14/g4Gz06zeeCxeuejlKIYQgSCm1Lt5jQArr9wPm2K+LA/F3TA+3y4rbrxOXJ6hjJ9lzQMEkzuVVkhC9pFynThSrX5+V771HdFSUR3WKF8/PN9/IIBshhM9Ea63rxXuM9rSiU6m3gWjge7tI3eIwnUR5aut4jSREL1FK0eSDDzh/+DCbv/7a43qNG1dg2LAuzJ+/jY8+mu3FCIUQInXsAS+dgN7mvyMzw4GS8Q4rAVh2eYlblCeo41QqCMiL+xbt7c7lVZIQvah069aUatmSVe+/T9S5cx7X69evKY89dg9ffLGIadPWeDFCIYRIGadS7YDXgc6m1vE3gp0J9HIqFexUqizuwTNrTK0jgAtOpRrZ/YN9gN/i1elrv+4O/Gkn2HlAG6dS+e3BNG3sMq+ShOhFSimajRjBldOnWTtyZIrq/fe/D9GkSUVefXUaq1cf8GKUQghxa06lpgArgcpOpcKdSvUHvgByAwucSm1yKvU1gKn1dmAasAOYC7xgah23I/pzwFjcA23282+/4zigoFOpfcDLwBv2uSKB94G19uM9u8yrZB6iD8x65BH2/fYbT+3bRy7D8Lje2bOX6dTpM86evcwffwymdOlCXoxSCHGnkZVqEpIWog80+e9/iY2OZsWwYSmqly9fKJMmPU1srKZPn7GcP3/FOwEKIYSQhOgL+cqVo/Zzz7F13DhO79yZorrlyhVmzJgnOHjwFM8+O4no6JjkKwkhhEgxSYg+0uidd8iWMyfL3nwzxXXvvbciH374MIsX72Lo0BlkhdvcQgiR0UhC9JHQwoVp8Prr7PvtN8KXL09x/d69G/HMM8359tvlfPPN4vQPUAgh7nBeH1RjGI4cwFIgGAgCfrYs11DDcBQAfgTKAIeAHpblOmPXSbBArGW5khxum9EH1cS5dukS31apQmiRIjy2dm2y20MlFhsby3PPTeb33zfx1VeP07VrHS9FKoS4E8igmoR80UKMAu63LFctoDbQzjAcjXAPr11kWa6KwCL7PYbhuGmBWMNwBN7qxJlN9pw5afbxx5zYsIFtEyakuH5AQACff/4oDRuWY/DgH1ixYl/6BymEEHcorydEy3Jpy3JdtN9msx+ahIu6TiThYq9TLcsVZVmuGwvEejtOX6nSqxfF772XZW++maLJ+nFy5MjG+PH9KF26EP36jWP37ggvRCmEEHcen/QhGoYj0DAcm4CTwALLcq0GilqWKwLAfi5iH+7Roq5KqQFxC9JGR0cn/jjDUkrR4vPPuXzqFCs/+CBV58ifPyfffz+AkJDsPProaCIizqZvkEIIcQfySUK0LFeMZblq416ProFhOKoncbhHi7pqrUfHLUgb5MEmvBlJsbp1qf7kk2z4/HMi9+xJ1TlKlCjA5MkDOH/+Co8/PkbmKAohRBr5dJSpZbnOAotx9w2eMAxHGID9fNI+zC+Luvpa0+HDCcqRg8WvvJLqc1SvXpwxY55gz57j9Os3nqtXk997UQghxK15PSEahqOwYTjy2a9DgFbALhIu6tqXhIu99jIMR7BhOG4sEOvtOH0tZ9Gi3DNkCAdmzeLg3LmpPk/z5lX47LNHWbFiH889JxP3hRB3LqdSE51K5Yv3Pr9TqfGe1vdFCzEM+MswHFtwL9K6wLJcs4CPgNaG4dgLtLbfY1mumxaItSxXlvyVrzNwIPkrVWLRSy8RfTX1GwJ361aXDz7oxrx52zDNH2UfRSHEnaqmqfXZuDem1meAuz2tLIt7+9nhhQv5qXVrGg8bRuOhQ9N0LpdrPiNHzuGpp5rx7rtdce+0IoQQt5bV5iE6ldoMNLcTIU6lCgBLTK1reFI/c41GyYJKt2pFlUceYfXw4dz16KPkr1gx1ecaPLg1Z89eYsyYpeTPnxOHo006RiqEEBneJ8BKp1I/2e8fBoZ7WllaiBnApePHGV+lCsXq16f7/PlpatnFxsby8stTmTZtLR980I1+/ZqmY6RCiKwkq7UQAZxKVQVa4J6xsMjU2uMdFTzuQzQMx8S4wTH2+/yG4fC4s1LcXs5ixWgyfDiHFy5k19SpaTpXQEAATmdP2ratzjvv/Movv6xLpyiFECJjciq13H6+AKzCPSblQ2CNU6nznp7H4xaiYTg2Wpbr7uTK/CGztxABYmNi+L5RIy4cPUq/XbvIkS9fms539ep1+vQZw4oV+xg1qg8PPFA7XeIUQmQdWbGFmBYp6UMMMAxH/ngLcBdIYX2RhIDAQFp//TXfN2jAsrfeovVXX6XpfDlyZOPbb/vTu/c3PP/8ZIKCAmjfvmY6RSuEuBPYUxY6ASdNravbZTdtzBBvEEuCjRlMrefZ5XWBCUAIMBsYZGqtnUoFA5OAusBpoKep9SG7Tl/gHTuUD0yt45b69JqUTLv4BFhhGI73DcPxHrAC+Ng7Yd2ZitWty90vvcTmUaMIX7YszefLmTOY774bQO3apXj22UnMn78tHaIUQtxBJuBeSCW+N3D3zSXYmMHuu0uwMYNTqbiNGUYBA3DPK68Y75z9gTOm1hUAFzDCPlcBYCjQEPda1kOdSuVPLlinUsFOpR51KvWWU6khcQ9PL9bjhGhZrklAd+AEcAroZlmuyZ7WF55p8sEH5ClThnn9+3P9StqXY8uVKwfffz+AatWKM2DABP76y+P+ZSHEHc7UeikQmag4yY0ZTK2jTK1vbMzgVCoMyGNqvdJ099FNSlQn7lw/Ay2d7lGFbYEFptaRdutzATcn5lv5zT5nNHAp3sMjKbrlaU+a356SOiJlsufKRduxY/mpVStWDB3KfR+nvRGeJ08IP/zwDD16fEW/fuOZMOEp7ruvcjpEK4TI5IKUUvFH3o3WWo9Opk5RU+sIAFPrCKdS8TdmWBXvuLiNGa7brxOXx9U5ap8r2qnUOaAgHm7ycAslTK09SZy3lGwL0TAcy+3nC4bhOB/vccEwHB6P3hGeK92yJTWffpp1n3xCxJr0WbUuX75Qpk59jvLli/Dkk+P4+++96XJeIUSmFh23SYL9SC4ZJuV2GzMktWFDauokZYVTKY8m4d9Ksi1Ey3I1sZ9zp/ZLRMrdN3IkB+fMYW6/fjy+fj1BwcFpPmeBAjn58cfn6N79Sx5/fAwTJvSnWTNpKQohUuSEU6kwu3XoycYM4fbrxOXx64Q7lQoC8uK+RRsONE9UZ7EHsTUBnnAqdRD35vQK0KbWHo0o9PiWqWE4goGHcI8sulHPslzveXoO4bngvHlp/c03/NqxI6v++1+avJc+/5kLFszFTz89T8+eo+jbdyzjxvXj/vvvSpdzCyHuCHEbM3zEzRsz/OBU6lPAwN6YwdQ6xqnUBadSjYDVQB/g/xKdayXuMSp/2qNP5wHD4w2kaQO86UFs7dNyYSkZZZqmzkqRcuU6dKBqnz6s+fBDjq9Lvwn2hQrl5qefXqBixaL06zeOefNk9KkQ4mZOpabgTlaVnUqFO5Xqj70xg1OpBBszmFrftDGDqXXcxgzPAWNxD7TZD8yxy8cBBZ1K7QNexh6xamodCbyPe0OItcB7dlmSTK0PA/mAB+xHPrvMIymZmL/NslxJbezrN1lhYv7tXD1zhgk1apA9d24eX7+ebKGh6Xbus2cv07v3N2zdGs5XXz1Op0610+3cQoiML6tNzHcqNQh4GvjVLnoQGG1q/X+3r/WvlLQQVxiGI9WdlSJ1cuTPT/uJE4nctYulr7+erueOG2hz992lefbZSfz66/p0Pb8QQvhYf6ChqfUQU+shQCPcCdIjKZl20QR4wjAcCTorLcsly594WemWLak7eDDrP/uMcp06UbZt23Q7d+7cOfjhh2fo23csL730PdeuRdOrV8N0O78QQviQwr1KTpwYbj1i9ZZSkhDT1Fkp0qbphx9yaMEC5j75JE9s3UpIwYLpdu6cOYOZNOlp+vUbx8svT+XSpSj692+WbucXQggf+RZY7VRquv2+K+5+So+kaPsnw3DUAuL2E1pmWa7NHlf2oqzchxjfyc2b+a5+fcp37kznn35K9w2Ar169zvPPT2bu3K28/HJbXnmlrWwyLEQWltX6EOHGuqn34m4ZLjW13uhp3ZQMqrllZ6VluTzqrPSmOyUhAqz5+GOWvv46bcePp8aTT6b7+aOjY3j11Wn8+OManniiCR988CABASnpahZCZBZZMSGmRUpumfYHGlqW6xKAYThG4B6O6/eEeCep98orHJo3j0UvvkhYw4YUqlo1Xc8fFBTIp5/2Il++UL75ZjFnz17ms88eIXt22dhECJGxOZWqB7wNlMad37wzMZ80dlaK9BEQGEiH775jUu3a/N6jB4+tWZOuUzEAlFIMGdKZggVzMXz4LM6du8yYMU8QGpr21XKEEMKLvgdeBbYCsSmtnJKE+C2w2jAcqeqsFOknV1gYHb77jp/btuXPgQNpO3Zsun+HUooXX2xJ/vyhvP76T/Tq9TWTJj1Nvnzpm3yFECIdnTK1npnayikdVJOgs9KyXB53VnrTndSHGN+yt99m9fDhdPjuO6r27u217/njj8288MJkypQpxHffDaBEiQJe+y4hhO9ktT5Ep1ItgUdw79MYFVduav3rbSvFk6KEmFHdqQkxNjqaH1u04OSmTTy+fj0FKlXy2nf9/fde+vcfT44c2Zg06Wlq1iyZfCUhRIaWBRPid0AV3NsUxt0y1abW/Typn5JRprfsrMwIE/Pv1IQIcCE8nEm1a5OrRAkeXbEi3fsT49u9O4LevUdz9uxlRo9+QhYFFyKTy4IJcaupdapXVEtJQtzNLTorLcuV5MKphuEoiXuH5GJ2vdGW5frcMBzDcE/jOGUf+pZluWbbdd7EPao1BhhoWa55SX3HnZwQAQ7MmcOvHTtStXdv2k+a5NW5g8ePn+Pxx8ewa1cEH330ML17N/LadwkhvCsLJsQxgMvUekdq6qckIS6P2xsxJQzDEQaEWZZrg2E4cgPrcQ/I6QFctCyXM9HxVYEpQAPcW4gsBCpZliuG27jTEyLAyvff5+8hQ7j/88+pM3CgV7/r4sWrPPPMRP76axeDB7fm1VfbywR+ITKhLJgQdwLlAe/uhwgMNQzHWBJ1VlqWK8nOSstyRQAR9usLhuHYCRRPokoXYKpluaKAg4bh2Ic7Oa5MQax3nEZvv83xdev46+WXKVy7NiWbeW/ptVy5cjBhwlO8+ebPfPbZAsLDz+B09pS5ikIIv3G6/yp/BvB4u6fEUvIL9iTuzspsxOus5N+Va5JlGI4ywN24N4m8F3jRMBx9gHXAK5blOoM7Wa6KVy2cWyRQpdQAYABA9uzZU3AZWZMKCKDDpEl816ABvz/8MI+vX0/uEiWSr5hK2bIFMnJkD0qWLMCIEbM5ejSSsWOfpGDBXF77TiGEuB17Y2GXqXXd1J4jJWty1bIsVz3LcvW1LNeT9sOjkTsAhuHIBfwCDLYs13lgFO6mbW3cLchP7ENvde/tpvu6WuvRWut6Wut6QUHSMgEIzpuXrjNmcP3yZWZ27050VFTyldJAKcWgQa356qvH2bz5KB06uNi1K8Kr3ymEEElY5VSqfmorpyQhrrL791LMMBzZcCfD7+NusVqW64RluWIsyxULjMF9WxTcLcL4Y/pLAFZqvvdOVPCuu2g/cSIRq1ez4Nln8cW0mq5d6/Drry9y7Vo0DzzwOQsWbPf6dwohxC20AFY6ldrvVGqLU6mtTqW2eFo5JQmxCbDJMBy7DcOxxTAcWw3DkewXGYZD4V7RZqdluT6NVx4W77AHgW3265lAL8NwBBuGoyxQEViTgjjveJW6deOeIUPYPmECaz7+2CffWbt2KWbPdlC+fGGeeGIcX331p0+SsRBCxNMe953H+4EHgE72s0c8GmVqJ7Wm3KKz0oNpF02AZSScrvEW7tUEauO+HXoIeMYegINhON4G+gHRuG+xzknqO2SU6c201sx65BF2//gjnX/5hUrduvnkey9fvobDMYXff99Ejx71GTGiB8HBcktbiIzIk1GmTqUcwFO4f6u34h5PEgr8CJTB/fvdw9T6jH18gmlzptbz7PK6wAQgBJgNDLL7/YJxT82rC5wGeppaH0rtNTmVyo+7IZUjrszUeqkndVMy7WK9ZblS3VnpTZIQb+36lStMa9GCU1u20GvZMorV9c0/n9Yal2s+Tudc6tUrw5gxT1C0aF6ffLcQwnPJJUSnUsWB5UBVU+srTqWm4U5mVYFIU+uPnEq9AeQ3tX7dqdQtp82ZWsc4lVoDDMI9aHI28D9T6zlOpZ4HappaP+tUqhfwoKl1z9Rcj1Opp+zvKAFsAhoBK02t7/ekfkr7EFPdWSl8L1tICF1nzCCkcGFmdO7MhWPHfPK9Silefrkto0f3Zft2i3btPmXt2oM++W4hRLoLAkKcSgXhbhlauKfHTbQ/n4h7bjl2+VRT6yhT64PAPqCBU6kwII+p9UrT3QqblKhO3Ll+Blo6Uz+xeRBQHzhsat0C96yGU0lX+VdKEmILYKVhOPanpA9R+FfOYsXo9vvvRJ0/z/QHHuDahQs+++5OnWrzxx+DCQ3NzkMPfcG33y6XfkUhMpYgpdS6eI8B8T80tT4GOIEjuGcDnDO1ng8UNbWOsI+JAIrYVYoDR+OdIm7aXHH7deLyBHVMraOBc0DBVF7PVVPrqwBOpYJNrXcBlT2tnJKEmKbOSuE/hWvW5IEff+TUli3M7N6dmGvXfPbdVaqEMXv2yzRvXoW33/6FwYN/4MoV332/ECJJ0XHT1+zH6Pgf2v1xXYCyuG+B5nQq9VgS57vdtLmkptN5NNXOQ+FOpfIBM4AFTqV+IwWzFDxOiPbgmfNAUdwLfMc9RCZQrkMH2nzzDYfmz2duv37o2BTvnZlqefOGMGFCf0yzHT//vJ4uXf7HkSOnffb9QohUawUcNLU+ZWp9HfdCLI2BE/ZtUOznk/bxt5s2F26/TlyeoI59WzYvEJmaYE2tHzS1PmtqPQz4D+4ZDl08re/x8D/DcNyysxJ3i1FkAjX69+fS8eMsf+cdcoaF0XzkSJ99d0BAAC+/3JZatUry4ovf0b79p3z55eM0b17FZzEIIVLsCNDIqVQocAVoiXtlsUtAX+Aj+/k3+/iZwA9OpT7F3aKsCKyxB9VccCrVCPdKZX2A/4tXpy/ufNId+NNMZd+KPWL1IdyjX+PyW23gPU/qp+SW6Y3OSstypbizUmQMDd96i9ovvMA6p5N1n36afIV01rJlVebMeZlixfLSu/doRo6cQ0yM71qrQgjPmVqvxj3QZQPuKRcBwGjcibC1U6m9QGv7PabW24FpwA5gLvCCqXXcxgzPAWNxD7TZD8RNpxsHFHQqtQ94GXgjDSH/hrtFGI07acc9PJKSaRdrLctV3zAcm4CGluWKMgzHJsty1U5xyOlMpl2kTGxMDL/37MneX36hw+TJVH0sqS4B77h8OYq33/6VH39cQ+PGFfjyy8dkaoYQPpYFd7vYZmpdPbX1U9JCDDcMRz7szkrDcKSos1JkHAGBgXT87jtKNm/OnCeeYO/06T6PITQ0GJfrET777BE2bjxC69ZOlizZ7fM4hBBZygqnUt7fIDg+w3Dch7vjc45lua6n9svTi7QQU+fahQv81Lo1JzZs4MGZMynbrp1f4tiz5zgDBkxg796TDBzYildeaUtQUKBfYhHiTpJVWohOpbbiHpkahLvf8gCp2A8xJbdMb9VZiWW5POqs9CZJiKl39exZpt1/P5E7d9Jt9mxKtWjhlzji30Jt1Kg8X331OMWKyS1UIbwpCyXEJGc8mFp7tEdiShLiXNwTJtfjXqMOAMtyfXLbSj4iCTFtLv/zDz/edx/nDx+m+/z5FG/c2G+x/PTTWt5442dCQrLx6ae9aNMm1d0BQohkZJWEmF5SkhC3WZYrQ/46SUJMu4sREfx4331cOnGCHosWUaxePb/FsnfvCZ57bhI7dlj07Xsv//lPZ0JDZRNoIdJbVkmITqWWm1o3cSp1gZsXAtCm1nk8OU9KEuJo4P8sy7U1xdF6mSTE9HH+6FGmNmtG1JkzdJ83j7CGDf0WS1RUNB999AfffLOYChWK8OWXj1OjRonkKwohPJZVEmJ6STYhGoYjyc5Ky3J51FnpTZIQ08/5I0f4sUULrpw6RbfZsynRpIlf41m6dDeDB0/h9OmLvP56B559tjkBASkZHC2EuJ2slhCdSuUAnse9f6/GvfXg13HrmybHk4SYZGdlcvsh+oIkxPR14dgxpt1/PxePHePBWbMo1by5X+OJjLzEa69NY/bsLTRpUpHPPnsUw8jn15iEyAqyYEKcBlwAvrOLHsG9NdXDntRP1bSLjEYSYvq7dPw401q25NzBgzw4cyalW7Xyazxaa6ZOXcN//vMr2bMHMXz4Q3Tpcjep3yVGCJEFE+JmU+tayZXdjtx7EreUs1gxevz1F/kqVODXTp3YP2uWX+NRSvHIIw2ZP9+kXLnCPP/8ZAYMmMDp0xf9GpcQIkPZaK+XCoBTqYbA355WlhaiSNKV06f5uW1bTm7aRLvx46nWp4+/QyI6Ooavv16M0zmH3LlD+Oij7nTs6NEfgEKIeLJKCzHexPxsuPc/PGK/Lw3s8HQ5N49biIbhWJGKOEUmF1KwID3/+su9zFvfvqz9xO/TTgkKCuTFF1syb94rFC+ej6efnsBzz00iMlL+KBLiDhW3P2873Hs33gc0t1939PQkKZl2scGyXHUSlTW1LNcyT7/MW6SF6H3RUVHMfuwx9vz8Mw1ef52mH36YIfrvrl+P4csvF+FyzSdfvlA+/rgHbdtmyOmyQmQ4WaWFmF5SkhAvAfOB7cA24AQw1rJc5b0XnmckIfpGbEwMi158kc1ff031fv1o8803BAR5vKWmV23ffoxBg35gxw6Lbt3q8u67XSlYMJe/wxIiQ5OEmFBKBtUcBIbj3seqLvAU8K43ghIZU0BgIK2++op7hg5l2/jx/PrAA0SdP+/vsACoVq04s2c7cDja8Pvvm2jW7CN+/nktWaGPXAjhG2m6ZZpRSAvR97aMHcuCZ5+lYNWqdJs1izylSvk7pBt2747ANH9k/frD3HdfZUaMeJhSpQr6OywhMpys1kJ0KrXC1DrVizGnpIV4X2q/RGQ9NZ96iu5z53L+8GG+b9iQ4+vW+TukGypXDmPGjIF88EE31q07RIsWH/P1138RHR2TfGUhRGaWI3GBU6mmnlb2+rQLw3CUBCYBxYBYYLRluT43DEcB4Efc20kdAnpYluuMXedNoD/uXTUGWpZrXlLfIS1E//ln+3Z+7diRyydP0vGHH6jYtau/Q0rg2LEzvPnmzyxcuIOaNUswcmRPWRNVCJsnLUSnUvmAsUB13FMZ+gG7SfT7bWp9xj4+we+3qfU8u7wuMAEIAWYDg0yttVOpYNw5oi5wGuhpan0oNdfjVOqWY11MrT0a6+KLifnRwCuW5boLaAS8YBiOqsAbwCLLclUEFtnvsT/rBVTDPYT2K8NwyG6xGVShatXovXo1hWrU4Ldu3Vg9YkSG6rcrXjw/Eyc+xddf98GyztKhg4thw2Zw4YJHSxsKIeBzYK6pdRWgFrAT+/fb1DrB77dTqZt+v51Kxf1+jwIG4F4Tu6L9ObiT5xlT6wqACxiRhljTNNbF6wnRslwRluXaYL++gPs/ZnGgCzDRPmwi0NV+3QWYalmuKMtyHQT2AQ28HadIvZxFi9Lzr7+o3KMHy954g1m9enEtA7XYlVJ07nw3S5a8Qa9eDRgzZinNmn3IjBkbMlTyFiKjcSqVB2gGjAMwtb5man2WZH6/Ta2jTK1v/H47lQoD8pharzTd/9NNSlQn7lw/Ay2dqZ/Tdc3Ueq2p9bem1q+aWvc2tZ7kaWWfLt1mGI4ywN3AaqCoZbkiwJ00gSL2YcWBo/GqhdtlCSilBiil1iml1kVHR3s1bpG8bKGhdJoyhWYjRrDn55/54Z57OHvggL/DSiB//pyMHNmT338fRJEieXj++cn07DmKvXtP+Ds0IfwlKO531H4MSPR5OeAU8K1TqY1OpcY6lcoJFDW1jgCwn5P7/S5uv05cnqCOqXU07o3oUzsKLk1jXXyWEA3DkQv4BRhsWa6kxurf6i+Dm/6M11qP1lrX01rXC8ogc+HudEopGrz2Gg/NmcPF8HC+q1ePQ/Pn+zusm9SpU5rZsx18+GF3tm4Np1WrkQwfPovLl6P8HZoQvhYd9ztqP0Yn+jwIqAOMMrW+G7iEfXv0Nm73+53U77pHv/meMLW+kJp6cXySEA3DkQ13Mvzesly/2sUnDMMRZn8eBpy0y8OBkvGqlwAsX8Qp0keZNm14bO1acpUowS/t27P6ww/RsbH+DiuBwMAA+va9l2XL3uTBB+vyxReLaNbsI2bP3iK3UYX4VzgQbmq92n7/M+4EecK+DYr9nNzvd7j9OnF5gjpOpYKAvEBkul+JB7yeEA3DoXDff95pWa5P4300E+hrv+4L/BavvJdhOIINw1EWd+frGm/HKdJXvvLleXTFCio9/DDL3nqLXzp25PKpU/4O6yaFCuXms88eYfr0l8iXL5SnnvqWnj1HsWOH/A0mhKn1ceCoU6nKdlFLYAfJ/H47lQp2KnXj99u+rXrBqVQju3+wT6I6cefqDvxp+umvUl9Mu2iCe9firbinXQC8hbsfcRpQCvfK5A9blivSrvM27qG90bhvsc5J6jtk2kXGpbVm8zff8NfgwYQUKkSnKVMo0dTjaUE+FR0dw6RJK/jkk7mcO3eF3r0b8eqr7SlUKLe/QxPCKzycdlEb97SL7MAB4EncjakEv9+m1pH28Ql+v02t59jl9fh32sUc4CV72kUOYDLu8SWRQC9Ta78MQJDtn4RPnNi4kd979ODcwYM0ef99Grz+OiogY27HeebMJVyu+UyYsJyQkOw4HG3o168p2bNLX7XIWrLaSjVpJQlR+EzU+fPMHzCA3T/+SJm2bWk/YQI5ixXzd1i3tXfvCd577zcWLdpJ2bKF+M9/OtO2bfUMscuHEOlBEmJCkhCFT8XdQl3scJAtVy7ajBmT4Va3Seyvv3YybNhv7N17giZNKjJkSBeqV79pJpAQmY4kxIQkIQq/+GfHDmY/9hgnN26ker9+3P/ZZ2TPnXH76q5fj2HyZHf/4tmzV+jWrQ6vvdaBkiUL+Ds0IVJNEmJCkhCF38Rcu8aKYcNYM2IEecqUocPkyRRvnOqF6n3i3LkrfPHFQsaNW0ZsbCx9+97LwIGtZe9FkSlJQkxIEqLwu/Dly5nTpw/nDx+m/muv0XjoUIJy3LRofYZiWWdxOucybdoacuYM5sUXW9K/fzNCQ7P7OzQhPCYJMSFJiCJDiDp/nsUvv8zWcePIX6kSbceNo0STJv4OK1m7d0fw4YezmT9/G8WK5eWVV9rRs2d9goJkPXqR8UlCTEgSoshQDi1YwPynn+b8kSPc/cILNP3wQ7Lnyvi3I1ev3s8HH8xi/fpDlCtXmFdeaUeXLrUJyKBTS4QASYiJSUIUGc61ixdZ9tZbbPziC/KUKkWb0aMp06aNv8NKltaauXO3MXLkHHbtiqBy5WKYZjvat68hiVFkSJIQE5KEKDKsY3//zbz+/YncvZuqjz3GfSNHZuh5i3FiY2OZOXMTn3wyj/37T1KtWnFefbUdrVtXkzmMIkORhJiQJESRoUVfvcqq//6XtR9/TFBICPd+8AG1n3uOgMCM30cXHR3D9Okb+PTTeRw+fJq77y7Fq6+25777KktiFBmCJMSEJCGKTCFyzx4WvfAChxcupGidOrT66ivCGjb0d1geuX49hmnT1vLZZ/M5duwMDRqUxeFoS7NmlSQxCr+ShJiQJESRaWit2f3TTyx2OLgYEUHNp56iyX//S2jhwv4OzSNRUdFMmbKK//u/hUREnKN27VIMGtSaNm3kVqrwD0mICUlCFJnOtQsX+HvYMDZ8/jnZcuak0TvvUGfgQIKCg/0dmkeioqL56ac1fPHFIo4ciaRqVYOBA1vRsWMtAgNl8I3wHUmICUlCFJnW6V27WGKaHPjjD/KWLUuzjz+m0kMPZZrWVnR0DDNmbOB//1vIvn0nKV++CAMHtqJr1zpky5bx+0hF5icJMSFJiCLTO7RgAYtffpl/tm2jeJMmtHC5KFavnr/D8lhMTCxz5mzhs88WsGOHRcmSBXjhhft5+OH6hITIyjfCeyQhJiQJUWQJsdHRbB03juX/+Q9XTp2ics+e3PveexSoVMnfoXlMa83ChTv4/PMFbNhwmIIFc9GvX1P69r2XAgXkN0ukP0mICUlCFFlK1PnzrBkxgvWffUZMVBTVn3ySe4YMIU/Jkv4OzWNaa1au3M+oUX+yaNFOQkKy06tXA555pjmlShX0d3giC5GEmJAkRJElXTpxgtXDh7P5668BqP388zR4801yFini58hSZvfuCEaNWsz06euJiYmlY8daPP98C2rVKuXv0EQW4GlCdCoVCKwDjplad3IqVQD4ESgDHAJ6mFqfsY99E+gPxAADTa3n2eV1gQlACDAbGGRqrZ1KBQOTgLrAaaCnqfWhdLxMj8mQNpEl5SxalPs//5x+e/ZwV+/ebPjf/xhbvjzL3n6by//84+/wPFa5chifffYIq1a9w7PPtmDx4l20b++ie/cvWbRoB7Gxsf4OUdwZBgE7471/A1hkal0RWGS/x6lUVaAXUA1oB3xlJ1OAUcAAoKL9aGeX9wfOmFpXAFzACO9eyu1JQhRZWt7SpWk3fjxPbN9O2fbtWf3hh4wpU4bFr77KpePH/R2ex8LC8vHOOw+wbt1QhgzpzMGD//D442No1uwjxo9fxsWLV/0dosiinEqVADoCY+MVdwEm2q8nAl3jlU81tY4ytT4I7AMaOJUKA/KYWq803bclJyWqE3eun4GWTj8NFZeEKO4IBatUofO0aTyxbRsVunZl/aefMqZsWRYNHMj5o0f9HZ7HcufOwbPPtmDlyrf54ovHyJcvlHfe+ZW6dd9l6NDpHDqUeVq/IkMIUkqti/cYcItjPgNeA+Lfjihqah0BYD/H9UUUB+L/DxVulxW3XycuT1DH1DoaOAf4pbNcEqK4oxSqWpWO331Hv127qPLoo2weNYqx5csz/5lnOLNvn7/D81j27EF061aXWbMGM2vWYFq1qsq33y7n3nuH06fPGJYu3U1WGB8gvC5aa10v3mN0/A+dSnUCTppar/fwfLdq2ekkypOq43OSEMUdKX/FirQbN47++/ZR46mn2D5hAuMqVWLGgw8Svnx5pkomdeqU5ssvH2fNmiEMHtyaTZuO0KvX17RoMYJJk/7m8uUof4coMq97gc5OpQ4BU4H7nUp9B5ywb4NiP5+0jw8H4g/pLgFYdnmJW5QnqONUKgjIC0R642KSI6NMhQAuHT/Oxi++YNOoUVyNjKRYgwbUN00qPvggAUFB/g4vRaKiopk5cyNjxy5l69ZwcufOwUMP1aNPn8ZUqRLm7/BEBpKSaRdOpZoDpj3KdCRw2tT6I6dSbwAFTK1fcypVDfgBaAAYuAfcVDS1jnEqtRZ4CViNe5Tp/5laz3Yq9QJQw9T6WadSvYBuptY90vtaPeH1hGgYjvFAJ+CkZbmq22XDgKeBU/Zhb1mWa7b9WYIhu5blmpfcd0hCFOnl2qVLbJ84kfUuF2f37SNPmTLUHTSI6v36EZwnj7/DSxGtNevWHWLy5BX8/vsmoqKiqV+/LH36NKZjx1rkyJHN3yEKP0tDQiwITANKAUeAh02tI+3j3gb6AdHAYFPrOXZ5Pf6ddjEHeMmedpEDmAzcjbtl2MvU+kC6XWQK+CIhNgMuApMSJcSLluVyJjq2KjCFf/+6WAhUsixXTFLfIQlRpLfYmBj2//476z75hGPLl5MtZ06qPv44tZ9/nsI1avg7vBSLjLzEtGlr+O67lRw4cIr8+UPp0aMBjz/emHLlMsduISL9ycT8hLzeh2hZrqV4fj+4CzDVslxRluW6MWTXa8EJcRsBgYFU7NqVR5Yt47G1a6ncowfbJ0xgYs2aTG3WjF1TpxJz7Zq/w/RYgQI5efbZFixd+gY//vgc995bkXHjltKkyXB69PiKWbM2ce1atL/DFMKv/Dmo5kXDcGwxDMd4w3Dkt8tuN2T3JkqpAXFDhaOj5X9k4T3F6tWj3fjxPBMezn1OJxcti1mPPMI3JUuy/J13OHvwoL9D9FhAQABNm1Zi9OgnWLduKG+80YFDh/5hwICJ1KkzjGHDZrBrV4S/wxTCL/yVEEcB5YHaQATwiV3u8fBbrfXouKHCQZls0IPInEIKFqT+K6/Qf88eHpozh7CGDVk1fDhjy5VjWqtW7PzhB65fueLvMD1WpEgeBg5szcqV7zB58tM0blyBb79dzv33f0yHDi4mTfqbc+cyz/UIkVZ+ySSW5ToR99owHGOAWfbb2w3ZFSLDUAEBlG3XjrLt2nH+yBG2T5zI1vHj+aN3b4Lz5qXKo49So18/itatmyn2ZgwMDKBly6q0bFmV06cvMn36eqZMWc0bb/zMsGG/0aFDTXr1akjjxuUJCJCZWiLr8sm0C8NwlAFmxRtUE2ZZrgj7tQNoaFmuXobhuOWQXRlUIzI6HRvL0SVL2DpuHHt/+YXoq1cpXLMm1fr2pUqvXuQyDH+HmCJaa7ZsOcrUqWuYPn09589fpWTJAvTs2YDu3evJrhtZhAyqScgXo0ynAM2BQsAJYKj9vjbu26GHgGfiJcgEQ3YtyzUnue+QhCgykqtnz7Jr6lS2jR/P8bVrQSlKtWjBXY8+SsWHHiJHvnz+DjFFrly5xty525g6dTXLlu0BoEGDsjz0UD06dapF/vzye5pZSUJMSCbmC+FFkbt3s3PKFHZ+/z1n9+0jMHt2ynXsSJVHH6Vcx45kCwnxd4gpEh4eyfTpG/j553Xs3XuC7NkDadmyKg89VI+WLasSHCz9+ZmJJMSEJCEK4QNaa46vW8euH35g19SpXDp+nOy5c1OhSxcqPvQQZdq2zVTJUWvN1q3h/PrreqZP38CpUxfImzeETp1q0b17PerXLyv9jZmAJMSEJCEK4WOxMTEcXbyYnd9/z74ZM7h65gxBoaGU69CBSg89RLmOHcmeO7e/w/RYdHQMy5fv5ddf1zN79hYuX75GiRL5efDBunTpcjd33RWWKQYX3YkkISYkCVEIP4q5fp3wJUvY88sv7J0+ncsnThAYHEyZNm2o2K0b5Tt3JqRAAX+H6bHLl6OYO3cbv/yyjqVL9xATE0uFCkV44IHadOlyN5UqFfN3iCIeSYgJSUIUIoOIjYnBWrGCvb/+yp5ffuHC0aOogACK33svZTt2pHynThSsWjXTtLZOn77I7NlbmDlzIytX7ic2VlO5cjE6d76bBx6oTYUKRZI/ifAqSYgJSUIUIgOK63PcP3MmB2bN4uSmTQDkKVOGcnZyLNm8OUE5cvg3UA+dPHmeP/5wJ8c1aw6itaZqVYPOnWvzwAO1KVtW1lP1B0mICUlCFCITuBAezoHZszkwaxaHFy4k+soVgkJDKd2qFWXatqVM69bkq1AhU7QeIyLO8scfm5k5cxPr1h0CoEaNEnToUJP27WtQsWLRTHEdWYEkxIQkIQqRyVy/coWjixdzYNYsDvzxB+cPHwYgT+nSlGnThtKtW1Pq/vsJKZjxJ8+Hh5/hjz82M2vWJtavd19HuXKFbyTH2rVLSXL0IkmICUlCFCIT01pzdt8+Ds2fz+EFCzjy119cO38elKJo3bqUad2a0q1bY9xzT4a/vRoRcZZ587YxZ85WVq7cR3R0LGFheWnXrgbt29ekUaNyBAUF+jvMLEUSYkKSEIXIQmKjo4lYs4bDCxZweMECrFWr0DExBAYHE9awISXuu4+S992Hcc89ZAsN9Xe4t3XmzCUWLtzBnDlbWbx4F1evXid//lBat65O+/Y1aNq0EqGh2f0dZqYnCTEhSYhCZGFR588TvmQJRxYvJnzJEk5u3IiOjSUgWzaK1a9Pyfvuo8R991G8ceMMO/fx8uUoFi/ezZw5W1iwYDvnz18lR45sNGlSkdatq9GyZVUMI5+/w8yUJCEmJAlRiDtI1LlzHPv7b44uWUL40qWcWLeO2OhoVGAgRevUwbjnHox77iHsnnvIUyrj9d9duxbNypX7WbhwOwsWbOfIEffe49WqFad166q0bl2NWrVKyio5HpKEmJAkRCHuYNcuXsRauZLwJUsIX76c42vXEn35MgA5w8IwGjUizE6SRevWzVDLy2mt2bv3BAsWbGfhwh2sXXuQ2FhNoUK5aNnSnRybNatErlwZu+/Un5JLiE6lSgKTgGJALDDa1Ppzp1IFgB+BMrg3aOhhan3GrvMm0B+IAQaaWs+zy+sCE4AQYDYwyNRaO5UKtr+jLnAa6GlqfSjdL9YDkhCFEDfEXL/OP1u3Yq1ceeNx7sABAAKyZaNI7dqENWpEsXr1KFq3LgWqVCEgMGMMdDlz5hJ//bWLBQu289dfOzl//irZswfSqFF57r//Lpo3ryJTOhLxICGGAWGm1hucSuUG1gNdgSeASFPrj5xKvQHkN7V+3alUVWAK/27htxCoZGod41RqDTAIWIU7If7P1HqOU6nngZqm1s86leoFPGhq3dNrF50ESYhCiCRdOnGCiFWr3Aly1aoErchsOXNSpHZtitoJsli9euSvVMnvSfL69RjWrTvIggU7WLhwO/v2nQTAMPLRvHkVmjevTNOmlcmbN+O0eP0hpbdMnUr9BnxhP5qbWkfYSXOxqXVlu3WIqfWH9vHzgGG4W5F/mVpXscsfses/E3eMqfVKp1JBwHGgsOmH5CR7tQghkpSzaFEqdOlChS5dAPcSc5G7dnFi/XqOr1vHifXr2TJ6NNFXrgDuJFm0Th2K1q1LkbvvpnCtWhS86y4Cs/tuVGi2bIHcc08F7rmnAkOGdCY8PJLFi3exePFufv99Ez/8sIqAAEWdOqXtBFmFWrVKEhh4x/U9Biml1sV7P1prPfpWBzqVKgPcDawGippaRwDYSTFuHb7iuFuAccLtsuv268TlcXWO2ueKdip1DigI/JOG60oVSYhCiBQJCAykULVqFKpWjWp9+gDu6R6nd+3ihJ0gj69bx+avvyb66lV3naAgCtx1F4Vr1nQ/atWicM2a5CxWzCe3MEuUKMBjjzXmsccaEx0dw4YNh28kyE8+mYfTOZf8+UNp2rQyzZtX5r77KhMWls/rcWUA0Vrreskd5FQqF/ALMNjU+rzz9v9mt/pAJ1GeVB2fk4QohEizgKAgClevTuHq1an+xBOAO0lG7tnDP1u2cGrLFk5u3kz4kiXs/P77G/VCChf+N0nWrEnBu+6iwF13EZwnj9diDQoKpEGDcjRoUI7XXuvA6dMXWbp0D0uW7GLx4l3MnLkRgPLli9CkSUWaNKlI48YVyJ//zhyM6VQqG+5k+L2p9a928QmnUmHxbpmetMvDgZLxqpcALLu8xC3K49cJt2+Z5gUivXIxyZA+RCGET12JjOSfrVs5tXkzp+xk+c+2bTduuQLkKl6cglWrUvCuuyhYtSoF7OfQQoW8GpvWmh07LJYt28Py5XtZtWo/ly9fQylF9erFadKkIk2bVqJBg7KEhgZ7NRZf8GBQjQIm4h5AMzhe+UjgdLxBNQVMrV9zKlUN+IF/B9UsAirag2rWAi/hvuU6G/g/U+vZTqVeAGrEG1TTzdS6h3euOGmSEIUQfhcbE8O5Awc4vXMnp3fsuPEcuXMn1+P9vx1SqNCNRFngrrvIX7Ei+StWJE+ZMgRmy5bucV2/HsPGjYdZvnwvy5fvZf36Q1y/HkO2bIHUqVPabkFW4u67S5E9e+a74eZBQmwCLAO24p52AfAW7qQ2DSgFHAEeNrWOtOu8DfQDonHfYp1jl9fj32kXc4CX7GkXOYDJuPsnI4FeptYH0vlSPSIJUQiRYenYWC6EhydIkKd37OD0jh1cPXPmxnEBQUHkKVPmRoLMZz/nr1iRPKVKERCUPsnq8uVrrFlzgOXL9/L333vZsiUcrTWhodmpX78sjRqVp1Gj8tSuXYrg4IyfIGVifkKSEIUQmY7WmsunTnF2717OxHvEvY/fqgzIlo185cqRr0IF8lWoQN6yZRM8sufKleo4zp69zMqV+1i2bC+rV+9n584IAIKDg6hTp/SNBFm3bukMeYtVEmJCkhCFEFmK1ppLx48nSJA3Eub+/TfmUMYJKVSIPGXKJEyU9vs8pUunaJeQyMhLrFlzgFWr9rNq1X62bTtGbKwmKCiAWrVK0ahRORo1Kk/9+mXJk8f/cyAlISYkCVEIcceIa1meP3iQcwcPcu7QIfez/Th/+DCx168nqJPLMMhTpgy5S5Z0P0qUuPE6T8mShBYpgrrN2qkXLlxl7dqDNxLk5s1HuX49hoAARbVqxWnUqBz16pWlXr0yfpnmIQkxIUmIQghhi42J4VJERIIkee7QIc4fOsSFo0e5EB5OTFRUgjqB2bOTq3jxfxNmoqSZu0QJQgoWRAUEcPnyNTZsOMzKlftYtWo/Gzce4epVdwIuXjw/9eu7k2O9emWoWtXw+v6PkhAT8npCNAzHeKATcNKyXNXtspsWhrUs1xn7swQLw1qWa15y3yEJUQjhC1prrvzzjzs5Hj3Kefv5wtGjXAwP57z9HBsdnaBeQLZs5CxWjJxhYeQKCyOnYZArLIzgwkWIvJ6dfSdj2HLoEqu3nSbixAUAQkOzc/fdpewkWZY6dUqTL1/67mEpCTEhXyTEZsBFYFK8hPgxEGlZro8Mw/EGkN+yXK8bhuOWC8Nalismqe+QhCiEyCh0bCyXTpz4N1EeO8bFiAguRURw0bK4ZL++cvr0TXVVYCA5ChUhNlc+LqpQTlwJ4th5xUUVwuWAUAqUKUGlOlWodW9N6jWpSoUKRdK01ZUkxIS8Pi7YslxLDcNRJlFxF6C5/XoisBh43S6falmuKOCgYTj24U6OK70dpxBCpAcVEEAuuyUY1qDBbY+Ljori0vHjNxLkxYgILlnWjeR5KSKCvOePUuLiKYhruGxyP6zx8CuBXAnKSY76zXlzxW++uLQsz18TZYpalisCwLJcEYbhSG5h2JsopQYAAwCy+3DRYCGESA9BwcHkLV2avKVLJ3lcbHQ0V/75h0snTnD5xAkuRkRwdOcBjmzfT/TBI+QuWTLJ+sJzGW3mqMeLvNorso8G9y1TbwYlhBD+EhAU5O5/LFbsRll1P8aTlflrr5MThuEIA7Cfk1sYVgghhPAqfyXEmUBf+3Vf4Ld45b0MwxFsGI6yQEVgjR/iE0IIcYfxxSjTKbgH0BQCTgBDgRkkWhjWslyR9vEJFoa1LNec5L5DRpkKIUTKySjThGRivhBC3KEkISbkr1umQgghRIYiCVEIIYRAEqIQQggBSEIUQgghgCwyqEYpFQtcScMpgnCPar1T3GnXC3LNdwq55pQJ0VpLw8iWJRJiWiml1mmt6/k7Dl+5064X5JrvFHLNIi3kLwMhhBACSYhCCCEEIAkxzmh/B+Bjd9r1glzznUKuWaSa9CEKIYQQSAtRCCGEACQhCiGEEEAWT4hKqXZKqd1KqX1KqTdu8blSSv3P/nyLUqpOvM8OKaW2KqU2KaXW+Tby1EvjNedTSv2slNqllNqplLrHt9GnTmqvWSlV2f73jXucV0oN9vkFpEIa/50dSqntSqltSqkpSqkcvo0+5dJ4vYPsa92eWf59waNrrqKUWqmUilJKmSmpK25Da50lH0AgsB8oB2QHNgNVEx3TAZgDKKARsDreZ4eAQv6+Dh9f80TgKft1diCfv6/J29ec6DzHgdL+viZvXjNQHDiIe0I2uLdhe8Lf1+TF660ObANCcU9gXwhU9Pc1pdM1FwHqA/8FzJTUlcetH1m5hdgA2Ke1PqC1vgZMBbokOqYLMEm7rQLyKaXCfB1oOkr1NSul8gDNgHEAWutrWuuzPow9tdLr37klsF9rfdj7IadZWq85CAhRSgXhThSWrwJPpbRc713AKq31Za11NLAEeNCXwadSstestT6ptV4LXE9pXXFrWTkhFgeOxnsfbpd5eowG5iul1iulBngtyvSVlmsuB5wCvlVKbVRKjVVKZYZ90tL67xynFzAl3aPzjlRfs9b6GODEvTF3BHBOaz3fi7Gmh7T8G28DmimlCiqlQnG3JEt6Mdb04sk1e6PuHS0rJ0R1i7LEc0ySOuZerXUdoD3wglKqWXoG5yVpueYgoA4wSmt9N3AJyAx9D2n9d0YplR3oDPyUjnF5U6qvWSmVH3droSxgADmVUo+lc3zpLdXXq7XeCYwAFgBzcd8+zAxrnXpyzd6oe0fLygkxnIR/CZbg5ltDtz1Gax33fBKYjvs2REaXlmsOB8K11qvt8p9xJ8iMLk3/zrb2wAat9QmvRJj+0nLNrYCDWutTWuvrwK9AYy/Gmh7S+v/yOK11Ha11MyAS2OvFWNOLJ9fsjbp3tKycENcCFZVSZe0WQC9gZqJjZgJ97BFqjXDfPopQSuVUSuUGsG8btsF96yWjS/U1a62PA0eVUpXt41oCO3wWeeql+prjff4Imed2KaTtmo8AjZRSoUophfvfeacvg0+FNP0bK6WK2M+lgG5kjn9rT67ZG3XvbP4e1ePNB+7+gj24R1y9bZc9Czxrv1bAl/bnW4F6dnk53LdWNgPb4+pmhkdqr9n+rDawDtgCzADy+/t6fHDNocBpIK+/r8OH1/wusAv3H3mTgWB/X4+Xr3cZ7j/uNgMt/X0t6XjNxXC3Bs8DZ+3XeW5XVx7JP2TpNiGEEIKsfctUCCGE8JgkRCGEEAJJiEIIIQQgCVEIIYQAJCEKIYQQgCREIYQQApCEKIQQQgCSEIXwGqVUK6XUZH/HIYTwjCREIbynFrDR30EIITwjCVEI76kFbFRKBSulJiilhtvrhwohMqAgfwcgRBZWCzgJzAPGaq2/83M8QogkyFqmQniBUiob8A9wGHhGa73SzyEJIZIht0yF8I6quLfhiQZi/ByLEMIDkhCF8I5awArce9F9q5Qq6ud4hBDJkIQohHfUArZprfcArwPT7NuoQogMSvoQhRBCCKSFKIQQQgCSEIUQQghAEqIQQggBSEIUQgghAEmIQgghBCAJUQghhAAkIQohhBAA/D+1m2j64rKZrAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEICAYAAADhmdstAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4C0lEQVR4nO3dd3xUdb7/8dc3mfReaCehE6oSUKS4WBHFBro2LLt2l93VdcedveveLf68966r66yzruV6sXcsa0FFWFEQFEFQQJrSSzgREpKQkD6Z7++PmcSQBDKTTM18no/HPGbOzDlnPqGcd77n+z3fo7TWCCGEEOEsJtQFCCGEEJ2RsBJCCBH2JKyEEEKEPQkrIYQQYU/CSgghRNizhLqAroqJidFJSUmhLkMIISJKTU2N1lpHXEMlYsMqKSmJ6urqUJchhBARRSlVG+oauiLi0lUIIUT0kbASQggR9iSshBBChD0JKyGEEGFPwkoIIUTYk7ASQggR9iSshBBChL2oC6sv53/C386ezf6tu0NdihBCCC9FXVjtWrWOmCWvsXfDtlCXIoQQwksBn8HCMKzPABcBB03TcUIHnyvgYeACoAa4wTQdXweqnvRemQCUHTgUqK8QQgjhZ8FoWT0HzDjO5+cDBZ7HbcD/BrKYrN45AFQcKAvk1wghhPCjgIeVaTqWAcdLhlnAC6bp0KbpWAlkGoa1X6DqyerrDqvKQ+WB+gohhBB+Fg59VnnAvlbLRZ73AiLbE1ZVhw4H6iuEEEL4WTjMuq46eE93uKJSt+E+VUh8fHyXviwxPQ2AmgoJKyGEaMuuVMs4A5vW7cYZ2JVqN87ApnXAxhk0C4eWVRHQv9VyPmB2tKLWeq7WeoLWeoLF0rWcjUtJAaC2orJL2wshRA/3HGE0zqBZOITVfOCnhmFVhmGdDBw2TUdxoL6sOazqKo8E6iuEECJi2bT2apyBTWtt03olkGlXKmDjDJoFY+j6q8CZQK5hWIuAe4A4ANN0PAEswN2c3I67SXljIOuJjY9Hq1gajkhYCSGikkUptabV8lyt9Vwftj/WOIOANTIgCGFlmo6rO/lcA78MdB2t6fgEGmvkLsNCiKjk1FpP6Mb2Xo8z8KdwOA0YdCohkabaGrQO+J+vEEL0NF6PM/CnqAyr2MRkYpsaqaysC3UpQggRaeYDP7UrpexKTQYO27QO6ClACI+h60FnSU4m7nAjZWVHyMhICnU5QggRNuxKtYwzsCt11DgDm9ZBH2fQLCrDKi41lTh9kLKyagYP7hXqcoQQImzYtD7uOAObu/8kqOMMIEpPAyampRKnGzl0SEYECiFEJIjOsMpII047KSuTEYFCCBEJojKsUjLTidONElZCCBEhojKsEtPdLSs5DSiEEJEhKsMqLjWVeKRlJYQQkSI6wyolBYurkUOlMpmtEEJEgqgNK4CyA3IDRiGEiARRHVblBw6FuBIhhBDeiOqwqiytkPkBhRAiAkR1WNFQx+HDtaEtRgghRKeiMqziU1MBiNONlJRUhbgaIYQQnYnKsGpuWcVpp4SVEEJEgCgPK2lZCSFEJJCwKpFrrYQQItxFdVglKiclJTLlkhBChLuoDKv4tDQAspIVBw9Ky0oIIcJdVIdVRiKUlkqflRBChLuoDKsYiwVLUhJp8ZqDByWshBAi3EVlWIG7dZUS2yQtKyGEiADRG1bp6STFNFFSUoXL5Qp1OUIIIY4jesMqLY14XY/T6aKiQqZcEkKIcBbVYRXnagCQa62EECLMRW9YpacT01gHwIEDElZCCBHOojes0tKgvgaA778/HOJqhBBCHE9Uh5WrVsJKCCEiQfSGVXo6jVVVZGYmU1wsYSWEEOEsesMqLQ1nbS39eqdKy0oIIcJcVIcVgJGbIGElhBBhLmrDKiE9HYA+GRYJKyGECHOWYHyJYVhnAA8DscBTpum4v83nGcBLwABPTXbTdDwbyJriPC2r3hkWSkqqcDqbsFhiA/mVQggREexKHXXMtml9f5vP2x2zbVoH9Jgd8JaVYVhjgceA84HRwNWGYR3dZrVfAptN01EInAn83TCs8YGsq7lllZMai8ul5Y7BQggB2JVqd8y2K9XhMdumdcsx265UQI/ZwTgNOBHYbpqOnabpaADmAbParKOBNMOwKiAVKAOcgSyquc8qM9G9LKcChRAC8ByzbVrvtGl93GO2XamgHbODEVZ5wL5Wy0We91p7FBgFmMAG4E7TdLSbXVYpdZtSao1Sao3T2b0/l+awSkvQgISVECJqWJqPo57HbW0+79Ix26Z1QGcED0ZYqQ7e022WzwPWAQYwDnjUMKzp7TbSeq7WeoLWeoLF0r3utnjPacDUOHcpcq2VECJKOJuPo57H3Dafd+mYbVeq3THbn4IRVkVA/1bL+bjTuLUbgbdM06FN07Ed2AWMDGRRzS2reFc9cXGx0rISQgg3r4/ZNq21TeugHLODMRpwNVBgGNbBwH5gNnBNm3X2AtOA5YZh7QOMAHYGsqjmsGo8coTevdMlrIQQwm01UGBXyqtjtl2poByzA96yMk2HE7gdWARsAV43Tccmw7DOMQzrHM9q/w2cahjWDcDHwO9M01EayLpi4+OJTUigoaqKvn0zKC6uCOTXCSFERLBp3e6YbdN6k12pOXaljjpm25VqOWbbtA7oMVtp3fZUZGRISUnR1dXV3drHY716Mfzyy/mXnsI33+xjxYo/+Kk6IYQIT0qpGq11Sqjr8FXUzmABkJCRQf3hw+TnZ2Ga5XJ7eyGECFPRHVaZmdRXVJCfn0VDQxMHD8qFwUIIEY4krDxhBVBUVB7iioQQQnQkqsMqMSuLuooK8vOzASgqKgtxRUIIIToS1WElLSshhIgMElYVFaSmJpKZmcz+/RJWQggRjqI6rBIzM3HW1uKsrycvL0tOAwohRJiK6rBKyMwEaBm+LqcBhRAiPElYAfUVFfTvn01RUTmRepG0EEL0ZBJW0DLIorq6noqKmtAWJYQQop2oDqtET1jVlZe3Gr4upwKFECLcRHVYtW1ZgYSVEEKEIwkr3GGVl+cOq/37ZUSgEEKEGwkroK6iguzsFJKT49m3T8JKCCHCTVSHVVxyMjEWC/UVFSilGDgwhz17JKyEECLcRHVYKaVaZrEAGDQol927S0JblBBCiHaiOqyAdmG1Z88hua+VEEKEGQmrzEzqPGE1eHAu9fVOTPNwaIsSQogexq7U83alMlstZ9mVesbb7aM+rBKPaln1AmD37tIQViSEED3SWJvWFc0LNq3LgfHebhz1YZWQlUVdmXtQxaBBuQDSbyWEEP4XY1cqq3nBrlQ2YPF2Y69X7KmScnKoO3QIAMPIICHBwq5d0rISQgg/+zvwhV2pNzzLVwD3ebtx1LesknJyqCsvR7tcxMTEMGBAjpwGFEIIP7Np/QLwY+B74ABwqec9r3gdVoZhfd4wrJmtlrMMw+p151i4SszJQbtcLYMs3MPXJayEEMIf7Ep95nmuAlYC9wN/Bb60K1Xp7X58OQ041jQdFc0LpukoNwyr151j4Sop191PVXfoEEnZ2QwenMvy5VvRWqOUCnF1QggR2WxaT/U8p3VnP76cBowxDGtL55hhWH3qHAtXSTk5ANR6+q0GDcqlrq6RAwe8DnwhhBAB5kvY/B1YYRjWNwENXAn8JSBVBVFHYQXu4et9+2aErC4hhOhJ7EolAJcBg2iVPTat/8ub7b1uWZmm4wXgctwdYyXAj03T8aIvxYajRE9YNY8IHDzYHVY7dsjwdSGE8KN3gVmAE6hu9fCKT6fxTNOxCdjkyzbhrm3LKj8/m4QEC9u3HwhlWUII0dPk27Se0dWNOw0rw7B+ZpqOqYZhrcJ9+q+ZArRpOtK7+uXhICEjAxUbS22pewRgbGwMQ4f2lrASQgj/WmFX6kSb1hu6snGnYWWajqme526N5AhXKiaGxOzsltOAAAUFffj66z0hrEoIIULHrtQM4GEgFnjKpvX9HaxzJvAPIA4otWl9Rie7nQrcYFdqF1CPp8Fj03qsNzV5fRrQMKwddo6ZpsOrzrFwlpST03IaEGD48L7Mn7+Ompp6kpMTQliZEEIEl12pWOAxYDpQBKy2KzXfpvXmVutkAo8DM2xa77Ur1duLXZ/fnbp8Gbrerc6xcNY2rAoKeqO1Zvv2gyGsSgghQmIisN2m9U6b1g3APNzH/tauAd6yab0XwKZ1pwdLm9Z7gEzgYs8j0/OeV3wZYJFvmo4udY4ZhvWoJqVpOto1KQ3DeiatmpSm6eisSek3iTk5VO7e3bI8fHhfALZvP8jYsf2DVYYQQgSDRSm1ptXyXK313FbLecC+VstFwKQ2+xgOxNmVWgqkAQ93NnWSXak7gVuBtzxvvWRXaq5N60e8KdqXltUKw7Ce6MP6ABiGtblJeT4wGrjaMKyj26yTibtJOdM0HWNwT3AYNEm5uUe1rAYNyiU2NoatW78PZhlCCBEMTq31hFaPuW0+72jqHt1m2QKcDFwInAf8ya7U8E6+92Zgkk3rP9u0/jMwGXd4ecWXltVU4AbDsB7VOWaajs46xyYC203TsRPAMKzNTcrNrda5BnjLNB17AUzTEdTzb61nXgeIj7cweHAu27bJiEAhRNQpAlqfUsoHzA7WKbVpXQ1U25VaBhQCW4+zXwU0tVpuouNg7JAvYdXVzjGvm5SGYV2Kp0npuQj5KEqp24DbAOLj47tYTntJubk46+poqK4mPiUFcI8I3LpVwkoIEXVWAwV2pQYD+4HZuBsUrb0LPGpXygLE4z6mOzrZ77PAKrtSb3uWLwGe9rYoX2awaNc55nmvM11qUhqGtV2TUms9t7nparH4b1rC5N7ugSw1B39o0A0f3pfdu0tpaHD67XuEECLc2bR2ArcDi4AtwOs2rTfZlZpjV2qOZ50twELgG+BL3MPbN3ay34eAm4AyoBy40ab1P7yty5eh6+06xwzDOtc0HZ11jnndpDRNRzVQbRhWb5qUftMSVgcOkDl4MOBuWTU1udi5s4SRI/sFowwhhAgLNq0XAAvavPdEm+UHgQd93O9XwFddqcmX5snNwCRPoGAY1geAL4DOwmo1UGAY1k6blIZh9aVJ6TfJffoA7rBq1hxQW7YUS1gJIUQ32ZWaAPwBGIg7e3y6KNiX0YBd6hwzTUe7JqVpOjYZhnWOYVjneNZp16Q0Tcdxm5T+lNIcVq1OAw4b1pv4+Fg2b94frDKEEKInexl3v9VluLuSLvI8e8WXltWzwCrDsPrcOWaajnZNStN0PNFm2ecmpb8k9eoFQHWrllV8vIWCgr5s2tT2jKUQQoguKLFpPb+rG3sdVqbpeMgwrJ8CP8LdorrRNB1ru/rF4cSSkEBCZuZRpwEBxozJY+nSLSGqSgghepR77Eo9BXyM+/InAGxav3XsTX7g6y1Cutw5Fu5S+vQ5qmUFMGaMweuvf0lJSRW9evXIeXyFECJYbgRG4p6lyOV5T/PDoL3j8mU0YIedY15cFBwRkvv06bBlBbBp037OPHNkKMoSQoieotCmtc+zIDXzZYBFtzrHwl1y795HDbAAGD3aANxhJYQQoltW2pUa3flqHfPlNGCJaTq63DkW7pL79KHm44+Pei8zM5m8vCwZZCGEEN03Fbg+4PezAu4xDGu7zjHTdHh1vjHcpfTpQ115OU0NDcS2msppzJg8Gb4uhBDdYFdKAT8DunxXW1/CqludY+EuudW1Vmn5+S3vjx5tsHjxJmpqGkhO9t98hEIIES1sWmu7Ug6b1id3dR++hFWhaTq63DkW7lrPD9g6rAoL++NyaTZt2s8ppwwOVXlCCBHpVtqVOsWm9equbOzLAIuVbe9D1ZOk9HXfcPFIcfFR748b557WcN26vUGvSQghepCzgC/sSu2wK/WNXakNdqW+8XZjX+9ndX0X7mcVEVLz3MPUq82jB1P06ZOBYWSydq2ElRBCdENXbzMFeBlWhmHtdudYuEvp2xeUomp/+8EUJ500UFpWQgjRDTat99iVygIKgMRWH3mVK16FlWk6tGFYHabp6HLnWLiLjYsjuXdvjnQQVuPGDeD999dz6NARcnJSQ1CdEEJENrtStwB34r5N1Drct7X/Ajjbm+197bM6xdcCI0laXl6HYTV+/AAA1q/f1+4zIYQQXrkTOAXYY9P6LGA8UOLtxr6E1VnAF4Zh3WEY1m8Mw7rBMKxed45FgtRjhNXYsf2JiVGsXdtjz4IKIUSg1dm0rgOwK5Vg0/pbYIS3G/sywKJbnWORINUwMFesaPd+SkoCw4f3lX4rIYTouiK7UpnAO8BHdqXKaX/X+GPy5RYhewzD2uXOsUiQmpdH7aFDOOvrsSQkHPXZ+PEDWLhwAy6Xi5gYXxqkQgghbFpf6nn5/+xKLQEygA+93d6XWde71TkWCZqHrx8xTTIHH30B8MSJQ3j11VVs23aAESPkNvdCCOELu1IJuCdCH8QP2TMO+C9vtvelidDSOWaaDp87xyJBS1h10G81adIQAFat2hnUmoQQood4F5gFOIHqVg+v+NJnVWeajjrDsGIY1gTTdHxrGFavO8ciwfHCauDAHPr2zWDlyp389Kc/CnZpQggR6fJtWs/o6sa+tKyKDMOaiadzzDCs7+JD51gkSDXc96/qKKyUUkycOIRVq3aitQ52aUIIEelW2JXq8vyyvgywaOkcMwyrz51jkSAxK4u4lBQq93Y86m/SpCHMn7+WoqJy+vfPDnJ1QggReexKbcB9hw4LcKNdqZ0E8n5WhmHtVudYJFBKkT5wIJV7Oh7g2NxvtXLlDgkrIYTwzkX+2IkvfVbvAoeBr2h188WeJn3QoGOG1ciRfcnISGLlyh1ccUWPnsxDCCH8wqa1Xy5v8iWs8k3T0eXOsUiRPnAgxV980eFnMTExTJ48lM8/3xbkqoQQIjLZlfrMpvVUu1JVuE8HqlYfa5vW6d7sx5ewWmEY1hNN07HBl0IjTcagQdSVl1NfWUlCevs/w9NOG86iRRvZvbuUQYNyQ1ChEEJEDpvWUz3Pad3ZT6dhZRjWozrHDMN6VOdYT7mfVbP0gQMBqNyzh14nth+4csYZ7tH6y5Z9J2ElhBBesiuVCPwC970RNbAceKJ5vsDOeNOy8kvnWKTIGDQIOHZYDRnSi7y8LJYt2yrXWwkheiS7UjOAh4FY4Cmb1vcfY71TgJXAVTat3+xkty8AVcAjnuWrgReBK7ypqdOwMk1Hj5n7zxvNLavDu3d3+LlSitNPH86CBd/Q1OQiNlbmCRRC9Bx2pWKBx4DpQBGw2q7UfJvWmztY7wFgkZe7HmHTurDV8hK7Uuu9rUuOtG0k9+mDJTHxmCMCAU4/fQSHD9fK/a2EEEH37LOf8c03AT32TAS227TeadO6AZiHe5qktu4A/gUc9HK/a+1KTW5esCs1Cfjc26IkrNpoudbqGC0rgKlTC1BK8emn3wWvMCFE1KuoqOFPf3qLhQsDOs4tD2idhkWe91rYlcoDLgWe6GxndqU22JX6BpiEexaL3XalduGeCP10b4vyOqwMw9r+Rk89VMaQIVRs337Mz3NyUhk/fgCLF28KYlVCiGi3dOm3uFyas88e3Z3dWJRSa1o9bmvzuepgm7ZzzP0D+J1N6yYvvu8i4GJgBjAYOAM40/P6Qq+L9nZFjr6HFQCGYT3NNB3LfdhHRMgaPpyiZcvQWqNUR39vMH36GB54YAEHDhymT5+MIFcohIhGixdvJisrhfHjB3RnN06t9YTjfF4E9G+1nE/7eWAnAPPs7uNjLnCBXSmnTet32u4sFBcFjzAM69vAJmAjcAB4Chja2YaGYT1qZIlpOjocWWIY1paRJabp6GxkScBkFRTQWF1N9fffk9qv43tXnXuuO6wWL97MtddOCXKFQoho09jYxMcfb+a8804I9MCu1UCBXanBwH5gNnBN6xVsWrfc8M+u1HPA+x0FlT/58hPvAu4DdgAnA7cA93a2kWFYm0eWnA+MBq42DGu7NqxnPV9GlgRMVkEBAOVbtx5znZEj+5Gfn8W//y2nAoUQgbdy5Q4OH67lvPO6PHG5V2xaO4HbcR+LtwCv27TeZFdqjl2pOQH98uPwpWXVYJqO1bhT1xcTge2m6dgJYBjW5pElm9us1zyyJOST7mUNHw5A+bZt9D/jjA7XUUpx7rkn8MorK6mpaSA5OT6YJQohoszChRtITIxrmZggkGxaLwAWtHmvw8EUNq1v8GafdqVW2LQ+tas1+dKy6vio3blOR5YYhtWrkSVKqduaOwWdTmcXy+lcWv/+xMbHU77t+HMATp8+hrq6RpYtk1GBQojA0VqzaNFGzjhjRCT/Ytxu3INdqdO83djrsDJNR5W367bh9cgS03Qcd2SJ1nqu1nqC1nqCxeJLo9A3MbGxZAwdetzTgABTpgwlMzOZ99/3+ro2IYTw2YYNRZhmBTNmBPYUYICNsCv1tl2p/7ErNduu1FnAc95uHIzrrLweWWIY1t3A5cDjhmG9JAi1HVP28OGdtqzi4y1ccMFYFi3aSG1tQ5AqE0JEm/ffX09sbAznnNOtIeuh1qVxD80C1zz5wWqgwDCsxxxZYpqOlpElhmF9DnjfNB3vBKG2Y8osKGDXwoVolwsVc+xMnzlzHK+8spJPPtnChRcWHnM9IYToCq017767ltNOG05OTmqoy+mOBpvWXRn3AAShZWWajnYjS0zTsckwrHMMwxqykSWdySoooKm+nsp9x5/W5NRTh5GTk8r8+WuDVJkQIpp89dUe9u0r49JLTwp1Kd3V1XEPQHBaVpimo93IEtN0dDiYwjQdNwSjps5kj3CPuCnbsoUMz+S2HbFYYrnookJee+1LqqvrSUlJCFaJQogo8M47X5OYGBfp/VXYtO7quAdA5gY8plzP7UFKN27sdN1Zs8ZTV9fIokWdryuEEN5yOpuYP38t06aNIi2t3WC6qCJhdQxJ2dmkGgYlGzqfMHLixMEYRiZvvNGlU7FCCNGh5cu3Ulp6hEsuifhTgN0mYXUcuSeeSKkXYRUTE8NVV01k2bKtFBWVB6EyIUQ0ePXVVWRnp3DOOWNCXUrISVgdR+6JJ3Jo82ZcXlyAPHv2JABef/3LQJclhIgCpaVVLFq0kcsvn0BCQlCGF4Q1Cavj6HXiiTTV11N+nNuFNOvfP5upUwt47bUvcblcQahOCNGTvfnmGhobm7j66kmhLiUsSFgdR/Mgi5JvvvFq/auvnsS+fWV89tnxLyYWQojj0VrzyisrOfnkQYwY0fGdH6KNhNVx5IwahYqN9arfCmDGjBPJykrmpZe+CHBlQoiebPXqXWzffpBrrpnc+cpRQsLqOCyJiWQVFHgdVomJccyePYkPP9wgAy2EEF329NPLSU9PZObMcaEuJWxIWHWi9/jxHPjqK6/Xv/FG9yTCzz33WaBKEkL0YEVF5SxY8A3XXDNFJhloRcKqE/0mTqSqqIgjZtu5dzuWn5/F+eef6LnPVX2AqxNC9DTPPrscgJtu8vruGVFBwqoTfSdOBOD71d5f8HvzzadTUVHDm2+uCVRZQogeqLq6npdf/oILLhhLfn5WqMsJKxJWneg9fjwqNtansJo4cTBjx+Yzd+6nNDXJMHYhhHfmzVtFZWUdt97arTlfeyQJq07EJSXRa+xYir/0/mJfpRR33HEOO3eWMH/+usAVJ4ToMerrnTz++BJOOWUwEyYMCnU5YUfCygv9Jk7k+9Wr0T5c7Hv++ScyYkRf/vnPj+QiYSFEp15//UuKiyu4667zQl1KWJKw8kLfiROpr6jo9M7BrcXExPCrX03nu+++58MPvRv6LoSITg0NTv75z8WcfPJATj99eKjLCUsSVl7oN9l9Yd7+FSt82m7mzHEMGdILh+Pf0roSQhzTG2+sYf/+cqzW81BKhbqcsCRh5YWcUaNI6tWLfUuX+rRdbGwMd911Hps3m7z99teBKU4IEdFqaxt46KFFjBs3gLPOGhnqcsKWhJUXlFL0P/NM9i1ditbap20vuWQ8J56Yz/33L6CurjFAFQohItVTTy2juLiCP/3pYmlVHYeElZf6n3kmVXv3cnj3bp+2i4mJ4Y9/vJj9+8tlVgshxFEOHTrCI48sZvr0MUyZMizU5YQ1CSsv9T/zTACfTwUCnHbacM46ayQPP/wRhw4d8W9hQoiI9dBDi6itbeSPf7w41KWEPQkrL3W136rZPffMorq6nv/5n/f8W5gQIiJt2FDE889/znXXTaGgoE+oywl7ElZeUkox4Kyz2PPRRz73WwEMH96Xn/3sTF577UtWrdoZgAqFEJHC5XLx+9+/SXZ2Cr/73QWhLiciSFj5YPAFF1BdXMzBtWu7tL3Vei55eVn8/vdv0tjY5OfqhBCR4uWXV/L113v4859nkpmZHOpyIoKElQ8Gn38+KMXODz7o0vbJyQn8939fyrffFvPEE0v8XJ0QIhIUF1dw333vM2XKUC67bEKoy+mQXakZdqW+syu13a7U3R18fq1dqW88jxV2pQoDXZOElQ9Sevem38SJ7Hj//S7v47zzTuCiiwqx2xeyebN3tx0RQvQMWmt+85vXaGho4sEHrwrLoep2pWKBx4DzgdHA1XalRrdZbRdwhk3rscB/A3MDXZeElY+GXHgh33/5JdUHDnRpe6UUf/3r5WRmJnPHHS9RX+/0c4VCiHD14osrWLr0W/70p4sZMqRXqMs5lonAdpvWO21aNwDzgFmtV7BpvcKmdfPt0FcC+YEuSsLKR0Mvdg8x3TF/fpf3kZOTit1+FVu2FPPggx/6qzQhRBjbtauE//qv+Zx++nCuv/5HoS7nePKAfa2WizzvHcvNQMAPZBJWPupVWEhWQQHfzpvXrf1Mnz6Ga6+dzP/+7xKWLNnip+qEEOGotraB2257nvh4C3//++xQn/6zKKXWtHrc1ubzjorrcAi0XamzcIfV7/xdZFsSVj5SSjFi9mz2LlnCkeLibu3r3nsvZdSoftx++0sUFZX5qUIhRLj54x/fYtOm/TzyyLXk5XXvDsClmzbhaurWaGKn1npCq0fb/qYioH+r5XygXQe7XamxwFPALJvWh7pTkDckrLpg1NVXg9ZsfeONbu0nOTmeJ5+8AafTxW23PS/9V0L0QPPmreLVV1dx553TmTat7TgF3zTW1vLKqaey5Ne/9k9xHVsNFNiVGmxXKh6YDRzV72FXagDwFvATm9ZbA1lMMwmrLsgZNYpehYVsfvnlbu9r8OBe/OMf17Bu3V5++9vXunTBsRAiPK1atZO7736DqVMLsNlmdHt/2995h4bKSgouu8wP1XXMprUTuB1YBGwBXrdpvcmu1By7UnM8q/0ZyAEetyu1zq7UmoAV5KEi9eCYkpKiq6urQ/b9axwOlt51F9evX0+vsWO7vb9//OPf/O1vH3LXXef55R+1ECK0du0q4aKLHiYrK5n58+8kOzul2/t849xzKd+2jVt37EDFdK2toZSq0Vp3v5ggswTjSwzDOgN4GIgFnjJNx/1tPr+WHzrojgA/N03H+mDU1lVjrr+e5b//Pev/7/8457HHur2/O++czu7dh3jooUUMGJDNlVdO9EOVQohQKCur5ic/eRLQvPDCrX4JqvJt29jz0Ueceu+9XQ6qSBbwn9gwrO0uMDMMa4cXmJmmI2gXmHVXUnY2I668ks0vvkjDke7PpK6U4m9/u8JzuuA1/v3vjX6oUggRbJWVtVxzzRPs31/OM8/c5LfrqdY9/jgxcXEU3tZ28F50CEY8TwS2m6Zjp2k6OrzAzDQdK0zTEdQLzPyhcM4cGqqq2PTCC37ZX3y8haeeupETTsjnttueY+nSb/2yXyFEcNTU1HPddXPZsqWYJ5+8gUmThvplvw1HjrDhmWcYccUVpPTt65d9RppghFVYXmDmD8aUKfSbNIk1djsup39G8qWnJ/HKKz+joKAvN930DMuXB2WgjRCim2pq6rnhhqf5+us9PPbYdZxzzhi/7Xvjc8/RUFnJ+Ntv99s+I00wwsrrC8wMw3rcC8yUUrc1X8jm9FM4dIdSiol3383hXbv4rpvD2FvLzExm3rw5DBqUy/XXP8WiRXJKUIhwVlFRw1VXPcGKFdt5+OFruOiicX7bt7O+ntUPPEDe1Kn0mzzZb/uNNMEIK68uMDMMa8sFZqbp6PACM6313OYL2SyWoIwN6dSwmTPJHjWKVffd190L9Y6Sk5PKG2/8gpEj+3HLLc/y2mtf+m3fQgj/OXiwkssue5QNG/bx5JM3+H0m9Y3PPktVURGn3nNPqGe+CKlghNVqoMAwrIMNw9rhBWaGYW25wMw0HRF13kvFxHDqn/9M6caNbH7pJb/uuzmwfvSjYVitr/Lww1278aMQIjA2bza58MJ/sGfPIV544VbOP7/7l7G01lhby6r77sOYMoUB06b5dd+RJuBhZZqOdheYmaZjk2FY5xiGtd0FZoZhXWcY1oBfYOZPI666ir6nnMLnf/wjjbW1ft13SkoCL7xwKz/+8ck88MACfv7zF6mpafDrdwghfPfhhxuYOfNhmppc/Otfv+T000f4/Tu+cjio2rePqX/5S1S3qkAuCvabfcuW8doZZzD1L39h8n/+p9/3r7Xm8cc/4b77PmDMGINnnrmJ/Pxsv3+PEOL4mppc/POfi3nwwQ8ZN24AzzxzE337Zvj9e44UF/N0QQEDp0/nkrff9tt+I/WiYAkrP3r3ssvYtWAB12/YQNawYQH5jsWLN/HLX75ETIzCbr+KCy8M+A06hRAeBw9WcvvtL/HZZ9u49NKTsNuvIikpPiDf9d5VV7Ht7be5cfNmvx5PIjWsou8y6ACa9sgjxCYksOiWW9AuV0C+45xzxrBw4V0MHtyLW299jv/4j9fltKAQQbB06becc86DrFmzmwcfvIpHH70uYEH13Ztv8t3rr3PqPfcE7BffSCMtKz/b8PTTLLrlFs5+5BFOCuA1EQ0NTv72tw95/PFPGDKkFw8+eCVTpsg/aiH87fDhWu69913mzVvFiBF9eeKJnzJiRL+AfV9NSQnPjhlDev/+XLNyJbFxcX7df6S2rCSs/ExrzVsXXcTexYu5+vPP6TvBv8NY2/rss2389revsWfPIa67bgp/+MPFZGQkBfQ7hYgWCxdu4Pe/f5PS0iP8/Odncddd55GY6N/waM3V1MS/ZsygaPlyrvvyS79Mkt2WhFWQhWtYAdQeOsQLJ52EUoqffPUVSTk5Af2+mpoG7PaFzJ27lNzcVO6++0KuvPIUYqJwsksh/GHbtgPce++7fPLJFkaPNnjoodmMHdu/8w27afkf/sCq++7jvKef5sSbbgrId0hYBVk4hxVA8erVzJs6lb6nnMLlH31EXFLgWzvr1+/lD394i6+/3sMJJ+Rx772XyKlBIXxQUVGDw7GIZ5/9jKSkeKzWc7npptOIjw/8JASbXniBD6+/nhNvuYXznnwyYN8jYRVk4R5W4O4kfe/KKxk2cyYz33yTmCDMuqG15p131vKXv7yHaVYwffoYfvOb84LyW6EQkaqqqo6nn17G//3fUior67j22sn8x3+cT25uWlC+f8cHH/DOrFn0P/NMfvzBB1gSEgL2XRJWQRYJYQWw9rHH+Pj22ym49FIufPXVgP4jbK22toEnn/yUJ55YSkVFjYSWEB2oqqrjuec+44knllBe7v5/8tvfns8JJxxvrm3/2j5/Pu9deSW5J5zAVUuWEJ8W2ICUsAqySAkrgK8ffZRP7riDgdOnM+utt4hPTQ3ad1dW1vLMM8uZO/dTKipqmDq1gFtvPYNp00ZJn5aIWkVF5Tz99DJeeWUlVVV1TJs2it/8Zgbjxg0Iah2bXnyRhTfeSJ+TTuLHCxaQnJsb8O+UsAqySAorgI3PP8+im24ie9QoLnnnnaBfO1FVVcfzz3/Os88up7j4MEOG9OLmm0/nsstOJj1dRg+Knk9rzapVO3nuuc/54AP3jcgvvngcP/vZGRQWBjekXE1NfP6nP7Hqr39lwNlnc8k77wS8RdVMwirIIi2sAPYsXsz7s2fjcjo596mnGHH55UGvobGxiQ8+WM+TT37K2rV7SUyM46KLCpk9exKTJw+R1pbocQ4erOT111czb94qdu4sIT09kWuvncJNN51GXl5W0Os5UlzMhz/9KXsWL2bsrbdy9iOPBK17ACSsgi4Swwrg8O7dzL/iCg6sWcOIK69k2qOPktzLP7e99tW6dXuZN28Vb7/9NVVVdQwcmMOPf3wyF188jhEj+kb9xJkiclVW1rJo0Ubee28dS5Z8S1OTi8mTh3L11ZO48MJCkpMDM/PE8Wit2fTCCyz59a9pqqtj2qOPcuLNNwe9DgmrIIvUsAJoamxk9d/+xop77yU+LY1T77mHwjlziI0P/n8gcF+ntXDhBubNW8WKFdtxuTQFBX246KJCCS4RMSora1m8eDPz569j6dItNDQ0kZeXxaxZ45k9exLDhvUOWW0Hvv6aT3/7W/Z+8gl5U6dy3tNPkz18eEhqkbAKskgOq2YlGzey5Ne/Zu/HH5M5bBiT//hHRl19dchCC6CkpIoFC77hvffW8cUXO9Bak5+fxdlnj+Lss0cxdWoBycnBO2UhxLForfn222I+/ngLn3yyhdWrd9HU5KJfv0xmzhzHxRePY/z4ASH9Rat00yZW/fWvbHn5ZZJycjj13nsZ9/Ofo0J4ul3CKsh6QliB+z/crg8/ZNndd1O6YQNp+fmc9OtfM+b664MyMuh4Dh6s5N//3sgnn3zLsmXfUVPTQHx8LJMmDWXq1AImTx5KYWH/oFwwKQTAvn1lfPHFdlau3Mmnn35HcXEFAGPG5DFt2iimTRvNyScPDGnfq3a52PvJJ6xxONi1YAGW5GRO+tWvmHT33SRk+P9WIr6SsAqynhJWzbTW7Fq4kC8feICiTz8lJi6OoTNncsINNzBw+vSgdsB2pL7eyZdf7uSTT7awdOm3fPfd9wAkJsYxYcIgJk8eysSJQxg7Nl9GFwq/aGpysWPHQdas2c3KlTv44osd7N9fDkBmZjKnnjqMadNGc9ZZIwNyPylflW3dyuaXXmLT889TtXcvSb16cdIddzDuF78I+JRrvpCwCrKeFlatlWzYwMZnn2XzSy9RW1JCXGoqg2fMYNisWQw891xSeofu3HuzQ4eOsGrVzpbfcjdvNmn+tzR0aG/Gjx9AYWF/xo0bwOjRRsBupSB6Bq01e/YcYt26vaxfv4/16/exYUMR1dX1AOTkpDJlylAmTRrClCnDGDmyb8hHrrqcTsyVK9kxfz473nuPsm+/BaUYdO65nHDDDQy75BIsiYkhrbEjElZB1pPDqllTQwN7Fi9m+7vvsv3dd6k5cACA7FGj6H/GGfQ/4wz6TpxIxuDBIR8AUVFRw7p1ez2Pfaxfv5cDByoBiIlRDBqUy8iR/TyPvowY0Y/Bg3OxWGJDWrcILq01paVH+O677/n222K2bv3hubKyDoCEBAtjxuRRWNifwsL+jB8/kGHDeof833hDVRXff/UV+5cvp2j5cswVK2isriYmLo7+Z5zB0JkzKbj0UtLy80NaZ2ckrIIsGsKqNe1y8f2aNexdsoSiTz9l/2ef0VBVBUBCRga9xo2jz/jx9CosJGv4cLKGDycpJyek/8GLiytYt24fGzcWtRycdu8uxeVy/5tLSLAwaFBuy2PIkF4MGpTL4MG5GEZmyH9zFl1XUVHDnj2l7NpVyp49h9i92/28dev3lJX98P82KyuZ4cP7MmJEX044IY/CwgGMGNE3pP2gzvp6Du/cSfm2bZR88w0l69dzcN06KrZvd6+gFLknnED+aafR/8wzGXTuuWHRF+UtCasgi7awasvldHJw/XoOfPUVB9eu5eDatZR88w3O2tqWdRIyM93BNWwYaf37k5qfT1qrR3Lv3kEflVRb28C2bQf47rvv2bKlmF27Sti1q4Q9ew5RX+/8ofYEC/n52RhGJoaRSb9+meTlZbYsG0YWaWnhd4olGjQ0ODlwoBLTrKC4uILi4sMUF1fw/feHKSoqZ/fuUioqao7apk+fdAYOzGHYsD6MHNmX4cP7MnJkP3r1SgvqL1Raa+oPH+aIaXJk/36OmCbVpknl3r2Ub9tGxfbtVO7dC62Oi5lDh9Jr3Dh6FxbS+6STyDv1VBKzgn8xsb9IWAVZtIdVR1xOJxWe3wgrtm2jbOtWKrZto3z7do7s34+rsfGo9WMsFpJyc4/9yMkhPi2N+PR09yMtjQTPa0tSkl8PMi6Xi+Liw+zeXeoJsFL27SvDNCswzQoOHKik7b/V1NQEevVKIycnlV690sjNTSUnJ63ltXs5lYyMJNLTk0hMjAv5qaRw5HK5OHKknrKyakpLj3DoUEePakpLqzhwoJKSkqp2+0hJSaBfvwwMI5PBg3sxcGAOAwe6W8kDBuT49SJcrTVNDQ00HjlCw5EjNHoedRUV1B06RG1ZGXVlZdQdOkRdWRm1zc+lpRwpLsZZU9Nun4nZ2WQVFJA5bBiZw4aR5XnOHTMmaNMgBYs3YWVXagbwMBALPGXT+v42nyvP5xcANcANNq2/DlDJgIRV1NAuFzUlJVQVFXGkqMj9vH8/taWlLY+a0lJqS0qoKytDu1zH3Z+KiWkJsLjkZGITE7EkJf3waL3c5rPY+Hhi4uKIiYsj1vPcdrn1OjFxcWhiKK+sp6S8lpJDNRw4eISDpUcoK6+hrKKWsrIaSg7VUF5RQ5NWaBRaeZ49j/gECxkZyaSnJ5GRkdQSYhkZSaSmJpKUFEdycgLJyfFHPZKSfnhOSoonPj6WuLjmh4W4uJignLJsanLhdLpoamqirs5JbW0DtbWNnucG6uoaWy03UlfXSE1NPVVVdVRW1lFVVcvhw7VUVdV53qulsrKWqqp69y8C+oc/rRhcxGgXyUlx5GQlkZudRHZGErk5SfTOTaVXTgq9cpLJzU4mJyuJpEQL2unE1dSEy+k86nVTfX27h7P5dV3dD687WMdZU3NUIDW/djmdnf55JWRkkJidTWJ2Nkk5OSRmZ5NiGKQaBml5eaQ0P/frR1xycsD//sJFZ2FlVyoW2ApMB4qA1cDVNq03t1rnAuAO3GE1CXjYpvWkQNYtF8hECRUTQ0qfPqT06QMnn3zcdV1NTdRXVFBXVkZDVZX7UVlJfWVly+sGz+v6ykqctbUtj6a6OmoOHmz3nrO2FmddnV9/pnTPY5CX62ulAIVWMaAULs+h2aXdr6uBIyh++PXth1aY5ugWWcuyarWuwtNyU7ifVMu6Cn7Yb0vrrvU+VEtgaE9ouJfxnJJyP5SmpcJjPnu2b11pqnL/WeWhUe5Saa6uOaQ41i8ou9v8OQIHPY9uUwpLQgKxHTwsCQnEpaaS3Ls38UOGEJeaSnxq6lHPrV8npKeTmJNDUk4OCZmZxMYF7vbzPdxEYLtN650AdqXmAbOAza3WmQW8YHO3dlbalcq0K9XPpnVxoIqSsBLtxMTGkuT5T+9PWmucdXW4GhpoamzE1erRermpoaHj91u91i6X+9HUhHa5cDU1gee59fu6qanD9zrapqnRibPRidPZhLOxCafT9cOy073c5GzC5XLhcmlcTS5c2vPs0u73m1zo5teesPkhRGgJHo07VDS09I/EKIiJiUF5HjExyvNatbwf41mOtcRiscQSa7FgiWt+HYvF09prfs8SZyEuLhYVGwtKucO09XNMTMvrGIsFFRtLjMVy1OuO3ouJjUV53mt53bx+69exse0CqPVyjMUip2bDTx6wr9VyEe7WU2fr5AESViLyKaWIS0qCJLloWIgQsiil1rRanqu1nttquaPfHtr2F3mzjl9JWAkhRHRxaq0nHOfzIqD1LcXzAbML6/iVhJUQQojWVgMFdqUGA/uB2cA1bdaZD9zu6c+aBBwOZH8VgFx1KYQQooVNaydwO7AI2AK8btN6k12pOXal5nhWWwDsBLYDTwK/CHRdMnRdCCGiSKReFCwtKyGEEGFPwkoIIUTYk7ASQggR9iSshBBChL2IHWChlHIBtZ2uGH4sQOcTm/Us8jP3fNH280Lk/sxJWuuIa6hEbFhFKqXUmk4uyOtx5Gfu+aLt54Xo/JlDKeLSVQghRPSRsBJCCBH2JKyCb27nq/Q48jP3fNH280J0/swhI31WQgghwp60rIQQQoQ9CSshhBBhT8IqCJRS/ZVSS5RSW5RSm5RSd4a6pmBRSsUqpdYqpd4PdS3BoJTKVEq9qZT61vP3PSXUNQWaUsrq+Xe9USn1qlIqMdQ1+ZtS6hml1EGl1MZW72UrpT5SSm3zPGeFssaeTsIqOJzAb7TWo4DJwC+VUqNDXFOw3In7NgPR4mFgodZ6JFBID//ZlVJ5wK+ACVrrE4BY3Pc/6mmeA2a0ee9u4GOtdQHwsWdZBIiEVRBorYu11l97XlfhPoDlhbaqwFNK5QMXAk+FupZgUEqlA6cDTwNorRu01hUhLSo4LECSUsoCJBPgO8aGgtZ6GVDW5u1ZwPOe188DlwSzpmgjYRVkSqlBwHhgVYhLCYZ/AP8BuEJcR7AMAUqAZz2nPp9SSkXcfYN8obXeD9iBvUAxcFhr/e/QVhU0fbTn7rie594hrqdHk7AKIqVUKvAv4Nda68pQ1xNISqmLgINa669CXUsQWYCTgP/VWo8Hqunhp4Y8/TSzgMGAAaQopa4LbVWiJ5KwChKlVBzuoHpZa/1WqOsJgh8BM5VSu4F5wNlKqZdCW1LAFQFFWuvmVvObuMOrJzsH2KW1LtFaNwJvAaeGuKZgOaCU6gfgeT4Y4np6NAmrIFBKKdz9GFu01g+Fup5g0Fr/Xmudr7UehLvD/ROtdY/+jVtr/T2wTyk1wvPWNGBzCEsKhr3AZKVUsuff+TR6+KCSVuYD13teXw+8G8JaejxLqAuIEj8CfgJsUEqt87z3n1rrBaErSQTIHcDLSql4YCdwY4jrCSit9Sql1JvA17hHva6lB05DpJR6FTgTyFVKFQH3APcDryulbsYd2leErsKeT6ZbEkIIEfbkNKAQQoiwJ2ElhBAi7ElYCSGECHsSVkIIIcKehJUQQoiwJ2ElhBAi7ElYCSGECHsSVkL4kVLqHKXUi6GuQ4ieRsJKCP8qxD2LgxDCjySshPCvQmCtUipBKfWcUuo+z5x5QohukLkBhfCvQtyzby8CntJa9/SZ5oUICpkbUAg/8dwGphTYA/xMa/1FiEsSoseQ04BC+M9oYDXu2cebQlyLED2KhJUQ/lMIrMB9/65nlVJ9QlyPED2GhJUQ/lMIbNRabwV+h/teR3EhrkmIHkH6rIQQQoQ9aVkJIYQIexJWQgghwp6ElRBCiLAnYSWEECLsSVgJIYQIexJWQgghwp6ElRBCiLD3/wGF1xaxoOKXoAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "D.k = np.linspace(0.05,0.1,1000)\n",
    "D.tau_plot(D.tau_h_dis(),D.tau_b_dis())\n",
    "\n",
    "D.k = np.linspace(1,11,1000)\n",
    "D.tau_plot(D.tau_h_dis(),D.tau_b_dis())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also note, that the damping time scales diverge periodically with $k=n\\frac{2 \\pi}{\\Delta x}$ with $n = 0,1,2... $. Within this periodicity, the damping timescales are decreasing with decreasing wavelength. Shorter waves are thus damped much faster than longer waves for both equations, with the biharmonic equation damping overall faster than the harmonic equation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.)\n",
    "To derive the damping time scales for the analytical solution, we substitute the solution as given in the hint into the harmonic equation.\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial u}{\\partial t} &= A_{2}\\frac{\\partial^2 u}{\\partial x^2}  \\\\\n",
    "e^{i k x} \\frac{\\partial \\hat{u}}{\\partial t} &= A_{2} \\hat{u} \\frac{\\partial^2 }{\\partial x^2} e^{i k x} \\\\\n",
    "e^{i k x} \\frac{\\partial \\hat{u}}{\\partial t} &= -A_{2} \\hat{u} k^2 e^{i k x} \\\\\n",
    "\\frac{\\partial \\hat{u}}{\\partial t} &= -A_{2} k^2 \\hat{u}\n",
    "\\end{align}\n",
    "\n",
    "As before, the harmonic damping time scale can be identified\n",
    "\n",
    "$$\n",
    "\\tau_{h} = \\frac{1}{k^2 A_2}\n",
    "$$\n",
    "\n",
    "The biharmonic damping coefficient can be derived analogously as \n",
    "\n",
    "$$\n",
    "\\tau_{b} = \\frac{1}{k^4 A_{4}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEICAYAAADocntXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzC0lEQVR4nO3deXxU5d3//9eVSUIWwpYActhlB9mEsql1V9xr626rtXorbW3t/Dqt3r+7m73dvu20o7a2FnFrrVq1VvharRsqtbgQBGRfZE1OCEsChCwkk1zfP2YSkxAgE5I5meT9fDzO48xZ5pyPEfLmnHNd1zHWWkRERBJJktcFiIiIxErhJSIiCUfhJSIiCUfhJSIiCUfhJSIiCSfZ6wKOV1JSkk1PT/e6DBGRhFJWVmattQl7AZPw4ZWenk5paanXZYiIJBRjTLnXNRyPhE1dERHpvBReIiKScBReIiKScBReIiKScBReIiKScOLW2tBx/E8AFwO7XDd0UhPbDfAQcCFQBnzTdUOfxqs+ERFJHPG88noKmH2U7RcAI6LTrcAf41CTiIgkoLiFl+uGFgFFR9nlMuDPrhuyrhv6COjhOP5+bVXPJwve5VdnXk3B59vb6hQiItJG2tMzr/7AjnrLedF1hzHG3GqMyTXG5IbD4RadbMuSFSS99wKf565q0fdFRMQ77WmEDdPEuibflGmtnQvMBcjMzGzR2zRzhgxkB7B7S15Lvi4i0qEFjalrpxCw9rB2CkFjDmunELA2bu0U2tOVVx4wsN7yAMBtq5OdMDxyquIdbXYKEZFE9hTtuJ1CewqvBcANjuM3juOfAex33VBBW52s/+gTATiQr/ASEWksYG2z2ikErLUBaz8CegSNabN2Co3Fs6n8c8AZQI7j+POAnwMpAK4behR4jcjl5yYil6A3tWU93ftkEzbJVO7a1ZanERHpqI7UTqHNLjrqi1t4uW7o2mNst8B341QOxhgOpWRSs3dPvE4pItKeJBtjcustz422J2iuZrdTaAvtqcFG3FVndKdm/9GuikVEOqywtXbqcXw/ru0UGuvU4WWyemCKdeUlItICC4Dbg8Y8D0wH9gesjcstQ+jk4ZXSI5vqgs1elyEi0u4EjalrpxA0pkE7hYC1cW+n0FinDq+0nN5Uh8uoqgqTktKpfxQiIg0ErD1qO4WAtXFtp9BYe2oqH3eZJ5yAjxoKPldHZRGRRNKpw6vHgEiXBHfDNo8rERGRWHTq8Oo1ODJ04u4tGpxXRCSRdOrw6nviIAD2btVtQxGRRNKpw6v/qKEA7M+PW+tOERFpBZ06vLIH98cCpYWFXpciIiIx6NTtw5N8Pg4lZ1KlIaJERBJKp77yAgind6OqeK/XZYiISAw6fXiZbr2gROMbiogkkk4fXinZfUgt3+d1GSIiEoNOH14ZJ5xAeriUkgNlXpciIiLN1OnDq/vA/vioYdvaLV6XIiIizdTpwytn6GAA8td+7nElIiLSXJ0+vPqOjHRULtyoKy8RkUTR6cNr8PgRABRpfEMRkYTR6cMrO3rbsCQ/bm+vFhGR49Tpwyu5SxcOpWRSsWun16WIiEgzdfrwAqjO6E51sYaIEhFJFAovIKlHDkajbIiIJAyFF5Ca04fUQweoqanxuhQREWkGhReQ2c8hs7qUXYX7vS5FRESaQeEF9Bg0gCQsO9Zu9roUERFpBoUX0HuYRtkQEUkkCi+g36gTAdj9+VZvCxERkWZReAEDxg4HoHjrDo8rERGR5lB4Ad0G9McCJfn5XpciIiLNoPACfCkpVHbpxqHCAq9LERGRZkiO58kcxz8beAjwAfNcN/RAo+3dgWeAQdHagq4bejIetdkevbHFu+JxKhEROU5xu/JyHL8PeAS4ABgLXOs4/rGNdvsusMZ1QxOBM4DfOI4/NR71pfbuR2pZMdXV6qgsItLexfO24TRgk+uGNrtuqBJ4Hris0T4WyHIcvwG6AkVAOB7FdR04kK7hEgoK9sXjdCIichziGV79gfrN+fKi6+r7PTAGcIGVwB2uGzrsUsgYc6sxJtcYkxsOt0629TpxCCmE2bpaHZVFRNq7eIaXaWKdbbR8PrAccIBJwO8dx9/tsC9ZO9daO9VaOzU5uXUe2/UbE3kp5Y6V61vleCIi0nbiGV55wMB6ywOIXGHVdxPwsuuGrOuGNgFbgNHxKG7wpMhpCtdtjMfpRETkOMSzteESYITj+IcC+cA1wHWN9tkOnA3823H8fYFRQFzu4/UZGe2ovGVbPE4nItLuBY1p0EI8YO0DjbYf1kI8YG1cWojH7crLdUNh4HbgDWAt8ILrhlY7jn+O4/jnRHf7X2CW4/hXAu8Ad7puKC5viUzPyaE6KZkyNy8epxMRadeCxhzWQjxoTJMtxAPW1rUQDxoTlxbice3n5bqh14DXGq17tN5nFzgvnjXVMsYQzsqmek+hF6cXEWlvpgGbAtZuBggaU9tCfE29fSyQFTQm7i3ENcJGPb7sE/CV7MHaxu1IREQ6nOTaVtvR6dZG21vUQjxgbVw6yyq86sno15/MqgMUFZV6XYqISFsL17bajk5zG21vUQvxoDGHtRBvCwqveroPGUzXmlJ2bNnpdSkiIl5rdgvxgLU2YG1cW4grvOrpO2oYAFtXbPC4EhERzy0BRgSNGRpthHENsKDRPrUtxAkaE9cW4gqvegZPHANA3ip1VBaRzi1g7WEtxAPWrg4aMydoTIMW4kFj6lqIB6yNSwtxk+iNEzIzM21paes8oyreuJHHR47kwIXf4Zf/fKRVjiki0h4ZY8qstZle19FSuvKqp9vgwVgMJVu3eF2KiIgchcKrHl9qKtVZ2YQL9UZlEZH2TOHViK+PQ0rJbg4diks/OxERaQGFVyNZQ4bSPbyfHTuKvC5FRESOQOHVSO9RI8msKWPzOg3QKyLSXim8Ghk0eRwAW3JXeVyJiIgcicKrkQETI4Mm71y9zuNKRETkSBRejfQcHnmv1/7NcekkLiIiLaDwaiStZ0+qUzOocHcce2cREfGEwqsJJsfBt6+Q6uq4jOwvIiIxUng1IWPgYLpV7SM/v9jrUkREpAkKrybkjBxOt+oDbNpQ4HUpIiLSBIVXEwaffBI+ati4RM3lRUTaI4VXE/pPjPT12r5M4SUi0h4pvJqQPSbyXq/iderrJSLSHim8mpDRty81qelU5muIKBGR9kjh1QRjDMnOYNIPFlJU1DovuhQRkdaj8DqC7sNH0itczMaNhV6XIiIijSi8jqD/pJPIqjnIhlVbvS5FREQaST7WDo7jLwFsE5sMYF031K3Vq2oHhkybxCpg6yfL4eazvC5HRETqOWZ4uW4oKx6FtDe1LQ4LV6/xuBIREWnsmOFVn+P4ewIjgLTada4bWtTaRbUHPYcPxxpD6ZbPvS5FREQaaXZ4OY7/FuAOYACwHJgBfAh0yHtqyWlpmF4nkFzkUlZ2iIyMLl6XJCIiUbE02LgD+BKwzXVDZwKTgd1tUlU70XXoMHqGi1i3bqfXpYiIdAhBY0zQmK8HjflZdHlQ0JhpsR4nlvCqcN1QBYDj+Lu4bmgdMCrWEyYSZ+JJ9AzvY83qPK9LERHpKP4AzASujS6XAI/EepBYnnnlOY6/B/AK8Jbj+IsBN5aTOY5/NvAQ4APmuW7ogSb2OQN4EEgB9rhu6PRYztGaBk+bzIbHw6z/cDl84xSvyhAR6UimB6w9OWjMMoCAtcVBY1JjPUizw8t1Q5dHP/7CcfzvAt2BfzX3+47j9xFJ13OBPGCJ4/gXuG5oTb19ehBJ5dmuG9ruOP4+zT1+W+gzcSIABZ8u97IMEZGOpCpojI9oF6ygMb2BmN/8G1Nrw1quG3q/BV+bBmxy3dBmAMfxPw9cBtRvi34d8LLrhrZHz7OrJfW1luxxkdHlD25aj7UWY4yX5YiIdAQPA/8A+gSNuRe4AvhJrAdp9jMvx/E/Hb0yql3u6Tj+J2I4V39gR73lvOi6+kYCPR3H/57j+Jc6jv+Gpg5kjLnVGJNrjMkNh8MxlBCb1K5dSertkFW6k/z8fW12HhGRTmQC8GPgfqAA+AowNdaDxNJgY4LrhvbVLrhuqJhIi8PmauqypfHIHcnAFOAi4Hzgp47jH3nYl6yda62daq2dmpzcoovHZus5Zhw5VXtYuzamx3siItK0cwPWrgtY+0jA2t8HrF0LXBDrQWL5zZ/kOP6e0dDCcfy9Yvx+HjCw3vIADm/wkUekkUYpUOo4/kXARGBDDOdpVYOnT2H3ordZvWIr5547zqsyREQSWtCYbwPfAU4MGvNZvU1ZwH9iPV4s4fMb4EPH8b8YXb4SuC+G7y8BRjiOfyiQD1xD5BlXffOB3zuOPxlIBaYDoRjO0eqcqZNJwrL5w6VELghFRKQFngVeJ3K78K5660sC1hbFejBjbVNj7jbNcfxjgTOJ3AJ8x3VDa2M5meP4LyTSDN4HPOG6oXsdxz8HwHVDj0b3+RFwE5HWJ/NcN/Tg0Y6ZmZlpS0vb7p1be9eu5cmxY1k64gqe2/Disb8gIpIAjDFl1trMo+0TNKZB96aAtYd1bwoacwb1ujcFrI1L96Zjhpfj+D9w3dCp9UaXr//syvNR5ds6vGrCYX6blkFu2gQeK/iArKy0Y39JRKSdO1Z4RZuzb6Be9ybg2oC1a+rt0wNYDMwOWLs9aEyfgLVHbSUejDTbvh44MWDtL4PGDAJOCFj7SSz1N2dU+VOj8045unxScjLpQ4aRk7eHzz7bwSmnjPC6JBGReJgGbApYuxkgaMwRuzcFrN0OcKzgivoDkTtrZwG/JDLCxt+JDD/YbLE0lf8/zVnXETlTJ5NTtZsVK3Yce2cRkcSQXNvlKDrd2mh7s7s3BY15L2jM0qAxTXZvamR6wNrvAhUQGWGDSBuHmMTSVP7cJtbF3LwxEQ2aOZ2uNWWsWrzC61JERFpLuLbLUXSa22h7i7o3BY05rHtTI/EZYcNx/LXNG4c5jv+4mzcmor5TpgBQsCTX40pEROKm2d2bAtaWAqVBY5rTvalVRthoTlP5IzZvdN1QzM0bE1GfSZPAGIy7maKiUnr1OmoDHRGRjmAJMCJozDG7NwWNaVb3pmhjjUXAUuBsIld3X4l2VI5Js5vKO47/e8Bf6o+y0R60dWvDWo8MGcaqnXD1v/7JGWeMbvPziYi0pWY2lW/QvSlg7b1BY+YABKx9NLpPg+5NAWsfPMYxlwasnXLc9ccQXvcQSd5PgSeAN1w31PxOYm0kXuH1ytXXsvyl/0vP3/6DO+5o6vGfiEjiaE54tYWgMY8ATwWsXXI8x2l2gw3XDf0EGAE8DnwT2Og4/vscxz/seApIFANmTKNrTSmf/eezY+8sIiJHcibwYdCYz4PGfBY0ZmWj4aKaJZbWhkSvtHZGpzDQE3jJcfy/ivXEiaa20Ubex0uIZVQSERFp4AJgGJF+XpcAF0fnMYmln9f3Hce/FPgVkVaG41039G0izSS/FuuJE03fyZPBGDKKtvP557u9LkdEJFH9EtgfsHZbwNptwAHg57EeJJaBeXOAr7puaFv9la4bqnEc/8WxnjjRpGZl0XXIME5wC1myZAvDh3v6kmcRkUQ1IWDtvtqFgLXFQWNieb0WENszr581Dq5622Ju5piIBp02i/7hApZ8stnrUkREElVS0JietQtBY2J9vRY05wv1BuStZfhigF7PB+aNp/6zZrHmz39mzb+XAtd6XY6ISCL6DbA4aMxL0eUrgXtjPUhzBubtlAPyNsWZOROAqk2r2Lv3INnZXT2uSEQksQSs/XPQmFwiDTYALm9JJ+WYWhs6jr+n4/inOY7/y7VTrCdMZNnjxuHLyKRfVQG5uVu9LkdEJOEEjbkS2BGw9vdEWqzfGzTm5FiPE0trw1uIDOvxBnB3dP6LWE+YyJJ8PpwZM+hfVcBHH33udTkiIonopwFrS4LGnEpkwPengT/GepBYrrzuIPK+lW2uGzoTmAx0ujbj/U+ZRU7VHj5etNrrUkREElF1dH4R8GjA2vm08StRKlw3VAHgOP4urhtaB4yK9YSJzpk5E4OlaMWnFBe3/bBUIiIdTH7QmD8BVwGvBY3pQoyPsIjxC3mO4+8BvAK85Tj++Rw+PH6H58yYAUC/Qy4ffqhbhyIiMbqKyGOn2dH+Xr2AH8V6kGYPzFuf4/hPB7oDr7tuqCrmA7SieA3MW99TEyawfNMBMn7wG+67r8MPLiIiHZBXA/O2lmZ3DHMcfxciw0ANqfe9SUSG+uhUBp15JrvW/JHX/t0p+maLiBy3oDEfBKw9NWhMk32HA9bG1Gc4ll7N84H9RF4idiiWk3Q0A884g08ffpiDa1dSWLifvn27e12SiEi7FrD21Oi8VfoOxxJeA1w3NLs1TproBnw50r1twKE8Fi3awJVXfsnjikREOpdYwmux4/jHu25oZZtVkyDSs7PJmTCBYZsLWLhwrcJLROQYgsb8f0fbHrD2t7EcL5bwOhX4puP4txC5bVg7tuGEWE7YUQw64wx2rXmUP7+3hnC4muRkn9cliYi0Z7W3C0cR6TO8ILp8CZEBMGISS3hdEOvBO7La515puzezdOk2pk8/0euSRETarYC1dwMEjXkTODlgbUl0+RfAi7Eer9nhdaTXoXRWA04/HZOUxNDKHSxcuFbhJSLSPIOAynrLlURascckpneoOI6/JzACSKtd57qhmC/3OoL0Xr3oN306Y9e5vPPOGv77vy/yuiQRkUTwF+CToDH/INJk/nIi4xvGRAPzHoehF1xA1r7tbFm5Cdfd53U5IiLtXsDae4GbgGJgH3BTwNr7Yz2OBuY9DkNmzwZrGVy5nXfeWeN1OSIiCSFg7acBax+KTstacgwNzHscTpgyhfScHMYl7+S11z7zuhwRkU4jlmdejQfmLSbGgXkdxz8beAjwAfNcN/TAEfb7EvARcLXrhl5qap/2wCQlMeS88yib/09e+fcGiopK6dUrYYcKExFJGLG0Nrw8+vEXjuN/l8jAvP9q7vcdx+8DHiHy8rE8YInj+Be4bmhNE/v9HyLP1Nq9IbNns/bZZ8lO38mbb67immume12SiEhCCBrzMnAQqAF2Baz9cXO/G1Nrw1quG3q/BV+bBmxy3dBmAMfxPw9cBjR+WPQ94O9Enq+1e0POPx+MYVKXnfzznysUXiIizbc6YO1PAYLGZMfyxZhfAHYc+gM76i3nRdfVcRx/fyLNJh892oGMMbcaY3KNMbnhcLjVC41FZp8+ODNnMrZ6K4sWbWD//nJP6xERSQRBYx4DTg0a86OgMaOBoli+H8/wMk2sa/wysQeBO103VN3Evl98ydq51tqp1tqpycktunhsVSO++lWSdm4hvbyIN99c5XU5IiLtXsDa/wIuBBYDsznGRUtjLfrN7zj+BvcpXTfUnPuUecDAessDOLzBx1TgecfxA+QAFzqOP+y6oVdaUme8jLj8ct4PBJjSxeUf//hUA/WKiBxB0JjFAWtnAQSsLQf+E51i0tLLltWuG/opgOP4m3ufcgkwwnH8Q4F84Brguvo7uG5oaO1nx/E/Bbza3oMLoMeJJ9J74kQm7dvBQ4vWs3Pnfk44Qe/4EhFpQlrjFUFjTgtY++9YDhLzbUPH8T8GnOo4/h85jr/Z9yldNxQGbifSinAt8ILrhlY7jn+O4/jnxFpHezPiq1/FbF9LWtVBXn55qdfliIi0V6OCxvwjaMw9QWOuCRpzJvBUrAcx1jZ+7HRsjuNPB04m0iJwjOuGbov5IK0kMzPTlpaWenX6OrtXruTpCRPYMP4KNmZPZeHCH2NMU4/5RES8Z4wps9YetWNq0JgGfXMD1jbZNzdoTF3f3IC1R+2bGzRmFZHhoU4CxgIO8EbA2j/HUn+zbxs6jn+x64ZmAbhuqMX3KTuqnJNOoufIkYwLb+bV9QNYuTKPCRMGHvuLIiLtUNCYw/rmBo1ZELB2TRP7xdI3tzJg7RIij5JaLJbbhofdp3Qc/2nHc/KOxBjDmOuvp2rdMrJ9Zbz44nH9fxER8do0YFPA2s0BayuB2r65jdX2zd3VzOOe3hrFxRJeoxzH/w/H8d/jOP5rHMffovuUHdnY668Ha7mo/z5eeimXsrLKY39JRMQbybX9ZaPTrY22H7NvbtCYZvXNra/2JZTHK5bw2gLcB3wOTAFuIfJqFInqMWwY/WbMYNCeFezfX878+S0aLFlEJB7Ctf1lo9PcRtub3Tc3YO1R++a2hViayle6bui471N2dGO//nXeuf12pk6dzVNPfcA110xTww0RSUTN7psbjPyOywEuDBoTDlj7SlsXF8uVV6vcp+zoRl11FUnJyZzfexcrV+axfPl2r0sSEWmJJcCIoDFDg8akEumbu6D+DgFrhwasHRKwdgjwEvCdeAQXxBBerhtqlfuUHV1G794MvfBCWLqQrIxknnrqA69LEhGJWcDaw/rmBqxdHTRmTtAYz/vmtqifV3vSXvp51ff5P//JPy6+mJLLf8jTn8JHH/1UI26ISLvSnH5e7Vk8B+btNIbOnk3WoEEM272UcLiGefMWeV2SiEiHovBqA0k+HxNuuYVdH7zHV07rx1/+spgDB/SqFBGR1qLwaiPjb74Z4/NxRvoWSkoqeOaZD70uSUSkw1B4tZGujsPwSy+l8P++yGkzB/PYY+9TXq5OyyIirUHh1YZOvuMOyvfs4ephJRQWHuAvf1nsdUkiIh2CWhu2IWstz3zpS1SWlPDGhO+xdl0hH330EzIzu3hdmoh0cmptKEdkjGHqD39I8YYN3Dwrjb17D/L442p5KCJyvBRebWzkFVeQNXAgxfOf4ZxzxvLHP77L/v1qeSgicjwUXm3Ml5LClB/8gLz33+e/zslh//5yfve7t70uS0QkoSm84mDCbbeR0acPBX9+hKuu+hKPPfY+mzfv9rosEZGEpfCKg9TMTKbdeSfb33mHb53Zk9TUZO6+e77XZYmIJCyFV5xMnDOHjL59WfvQr/nBD87lrbdW895767wuS0QkISm84iQlI4Ppd93FjnffZfYIGDo0h5/85GUqKqq8Lk1EJOEovOJowm230bV/fxb/z13ce89X2bx5Nw8++KbXZYmIJByFVxylpKdz2v33U5ibS++8XK6+ehqPPLKQlSvzvC5NRCShaISNOLM1NTwzfTqlBQVc8clyzrnwYXr37sZrr/lJSfF5XZ6IdBIaYUNiYpKSODMU4mB+Phvm/YH77ruC1avzCYV0+1BEpLkUXh4YcOqpjLzySj554AFmjcrkqqum8fDDb/HRR597XZqISELQbUOPlOTn8+SYMfSbMYMLXl7A7Nm/pby8irffDtCzZ8JeyYtIgtBtQ2mRrP79Oe3++9n21ltsn/93/vCHb7BnTwl+//PU1NR4XZ6ISLum8PLQxDlz6DdjBu/+4AcMd9L56U8v5c03V/Hb3+r5l4jI0Si8PJTk83He3Lkc2r+ft267jW9961SuuupL/Pa3b/D66595XZ6ISLul8PJY7/HjOfXee9n48susfvppHnjgSiZPHsT3v/8s69YVeF2eiEi7FNcGG47jnw08BPiAea4beqDR9uuBO6OLB4Fvu25oxdGOmagNNuqrqa7mxXPOYWduLjcuX055RjYXXhgiKSmJBQu+T//+Pb0uUUQ6GDXYaCbH8fuAR4ALgLHAtY7jH9toty3A6a4bmgD8LzA3XvV5Kcnn44KnnybJ5+PVa6+ld690nnnmVg4erOC66/5EcXFih7OISGuL523DacAm1w1tdt1QJfA8cFn9HVw3tNh1Q8XRxY+AAXGsz1PdBg1i9pNPsnPJEhZ+//uMG9efJ5+8mW3b9nDDDfMoKzvkdYkiIu1GPMOrP7Cj3nJedN2R3Ay83tQGY8ytxphcY0xuOBxuxRK9NeLyy5l21118NncuKx9/nFmzhvPII99g2bJtfP3rj1FaqgATEYH4hpdpYl2TD9wcx38mkfC6s6nt1tq51tqp1tqpycnJrVii90695x4Gn3MOb3/3uxR88gkXXTSR3/3u63zyyWa+/vW5HDxY4XWJIiKei2d45QED6y0PANzGOzmOfwIwD7jMdUN741Rbu5Hk83HRc8/R1XH4xyWXsG/zZi6//GQeeeQb5OZu1TMwERHiG15LgBGO4x/qOP5U4BpgQf0dHMc/CHgZ+IbrhjbEsbZ2JSMnh6+9/jo14TB/nz2bsj17uOyyyTz66A189tkOvvKV35GXV3zsA4mIdFBxCy/XDYWB24E3gLXAC64bWu04/jmO458T3e1nQDbwB8fxL3ccf2686mtveo0axeULFnBg+3ZeufRSqsrKuOiiiTz77BwKC/dz6aUPsXp1vtdlioh4QgPztnMb/v53Flx5JYPPPpuvLFhASno669YVcP31f6KkpIKHHrqeCy4Y73WZIpJgEr2fl8IrAax6+mn+ddNNDDnvPL7yyiskp6Xhuvu45ZYnWb58O3fccS6BwGx8Pg2YIiLN05zwChrTYGCJgLUPNNp+2MASAWuPOrBEa9FvuwRw0o03cv68eWx94w3mf/WrVJWV4Tg9ePnl27nmmuk89NBb3HDDY+zZU+J1qSLSQQSNOWxgiaAxTQ4sEbA27gNLKLwSxPhvfYvzHnuMLf/6Fy+ddx4VxcWkpaXwm99czQMPXMHixZs4++xf8+67a70uVUQ6hmnApoC1mwPWNjmwRMDaxQFrPRlYQuGVQCbccguXvPACO5cs4bnTTqMkPx9jDDfccAqvveYnO7sr118/l5/85GXKyyu9LldEElurDSzRFhReCWbUFVfwtddfp2T7dp6dOZPCZcsAGDPG4bXX/Nxyy5d54ol/c/bZv2bRovUeVysi7Vhy7UhF0enWRtubPbBE0JijDizRFhReCWjQWWdxzaJFYC3PnXIK6/72NwDS0lL45S8v54UXvkNSkuGaax7le997hr17D3pcsYi0Q+HakYqiU+PnVc0aWCJoTN3AEgFr4zawhFobJrDSwkIWXHEF+R98wLQ77+TUe+4hKTpcVkVFFQ8//BaPPLKQzMwu+P3nceONp5Ca2rGG0xKRljlWa8OgMcnABuBsIJ/IQBPXBaxdXW+fQcBC4IaAtYvbuOQGFF4JrrqykoXf/z4r/vQnnFmzuOivf6X7kCF129evL+BnP3uFf/97A0OH5vA//3MJF1wwHmOauiMgIomiqqwMX2pq3T9YY9XMpvIXAg8SaSr/RMDae4PGzAEIWPto0Jh5wNeAbdGvhAPWTm1RQTFSeHUQa599lre+/W2MMZw7dy6jr7qqbpu1loUL1/K//7uADRsKmTJlCD/84fmcfvoohZhIAsr/z394/ZvfZPzNNzP9rrtadIxE76SsZ14dxJjrruOG5cvpNWYMr159Na9edx1lu3cDYIzh7LPH8vbbP+JXv7qSgoJ9XHfdn7jkkod49921JPo/YEQ6i6ryct4LBHjutNOoCYfpN2OG1yV5RldeHUx1VRUf338/H91zD126dePMBx9kzPXXN7jCOnQozAsvfMLDD79Nfn4x48cP4L/+68tccslkunTRMzGR9mj7woW8/Z3vULR+PRPnzOH0X/+a1K5dW3y8RL/yUnh1UHtWr+aNW26h4KOPGHzOOZz54IPkjBvXYJ/KyjAvvpjL3LnvsXFjIX36ZHHjjafyjW/MJCcny6PKRaS+kvx83g8EWPf883QfOpRzH32UIeedd9zHVXh5TOF1ZDXV1Sz/4x/5z09/SmVJCRPnzOGUu+8mPTu7wX7WWt5/fz3z5i1i4cK1pKT4OP/8k7jmmumcfvoojZko4oGqsjI+/d3v+Oiee6ipqmLaXXcx7c47SUlPb5XjK7w8pvA6trI9e1j8i1+w4tFHSc3KYtqddzL59tubvOWwcWMhf/3rh7z0Ui5FRaX069eDq6/+El/72lSGDevjQfUinUtNOMzKJ57gw7vv5qDrcuJFF3HWQw/RY9iwVj2PwstjCq/m27N6Ne//+Mdsee010nNy+NKPf8yk73yH1MzD//xWVoZ5663VPPfcx7z33jpqaixjxzpcdtlkLr54IkOH9vbgv0Ck46oJh1n/wgssvvtuijdswJk5k9Puv5+Bp5/eJudTeHlM4RU796OPWPzzn7P1zTfJ6NOHk++4g4m33XbY7cRaO3fu59VXlzN//nKWLt0KwPjxA5g9ezxnnz2G8eMHqMm9SAtVlZWx6sknWRIMcmDrVrLHjeO0++5j2CWXtOnfK4WXxxReLZe/eDEf/vKXbH3jDZLT0xl3441M+cEP6DVq1BG/k5dXzD//uYJXX13Op59ux1rLCSd056yzxnDOOWM57bSRZGZ2ieN/hUhiOlhQwGePPcay3/+e8t27cWbOZNqdd0ZCK6ntnzMrvDym8Dp+u1et4tMHH2TNX/5CdWUlQ84/n/G33MLwSy/Fl5p6xO/t2VPCwoVrefvtNbz//npKSipITk5i0qRBnHLKCGbNGs7UqUNITz/yMUQ6E1tTw/aFC1n+6KN8Pn8+NeEwQy+8kOl33UX/U0+N6x0MhZfHFF6tp7SwkBWPPsrKefMoycsjPSeHsTfcwPibbyZnbON30DVUWRnmk0+2sGjRehYv3sSKFTuorq4hNdXHlClDmDFjGFOnDmHSpEH07Jmwf19EWmT/1q2sffZZVj/1FMUbN5Kenc24m25i4q230nPECE9qUnh5TOHV+mqqq9n65pusfPzxun8d9p44kdFXX82oq6+mx4knHvMYJSUVfPzx5/znP5tYvHgTq1bl143kMWxYH04+eRAnnzyEk08ezOjR/UhJ8bX1f5ZIXJXt2cP6F15g7V//irs4MmbtgNNOY8KttzLyiitITkvztD6Fl8cUXm2rdNcu1j33HOv/9jfcDz8EoO/UqYy68kqGXXIJvUaPbtatjoMHK1i+fAfLlm1j6dJtLF26te5VLampPkaN6sfYsQ7jxvXnpJP6M2aMQ/furdOfRSRe9m/bxucLFrBp/nzy3n+fmnCY7HHjGHv99Yy+9toGg2Z7TeHlMYVX/Ozfto0NL77Iur/9jcLcXAC6Dx3KiRdfzIkXXcTA009v9r8mrbXs2FHE0qXbWLkyjzVr8lm1Kp+ioi/+Xw4c2IsxYxxGjOjD8OF9GT48MleoSXtREw6zMzeXLa+/zqb589m9YgUAvcaMYfillzL62mvpPWFCu2yNq/DymMLLGwe2b2fza6+x+dVX2f7OO4QrKkjOyGDAqacy8MwzGXTmmfSdMiWm1zVYayksPMCaNS6rV+ezerXLunUuW7bsoaqqum6/Pn2yGD68LyNG9GXYsN4MHpzD4MHZDBjQi4wMNQ6RtmOtZe/atWx/5x22vf02O957j8oDBzBJSfQ/5RSGXXYZwy+91LPnWLFQeHlM4eW9qvJydrz7Lltef53t777L3tWRd9WlZmXR/7TTGHj66fSbMYO+U6Y02SH6WMLharZvL2LjxkI2bSpk06ZdbNpUyMaNhRw4UNFg3z59shg0KJvBg7MbzB2nB337dtfAwxKT6spKdi1bhvvhh+QvXkz+Bx9QWlAAQI9hwxh09tkMPuccBp111hH7SbZXCi+PKbzan9Jdu9jx3nvsePdddrz3HkXr1gFgfD5yTjqJftOn1029Ro8mydeyxhrWWvbuPci2bXvZvr2I7dv3Rj9HJtfdR01Nwz/fOTld6devR3TqjuNE5rXr+vbtpn5qnZStqaF40yZ2LVvGztxcCj78kJ25uVQfOgRAt8GDcWbNYtBZZzHo7LPpMXSoxxUfH4WXxxRe7V/Z7t0UfPxx3bTzk084tH8/AMlpaWSfdBK9J0ygz8SJ5EyYQO8JE0jv1eu4z1tZGSY/v5jt24soKNhHQcH+6PyLz8XFZYd9Lz09ld69u9K7dxY5OVkN5o3XZ2WltcvnGXJ0laWlFK1dy67ly9m1bBm7li9n94oVVEV/l/hSU+k7ZQr9Zs6k/6xZODNn0tVxPK66dSm8PKbwSjy2poaiDRso+Phjdq9Ywe7PPmP3ihWU79lTt0/X/v3JHjOGnqNG0WvkSHqOHEnPUaPoNmhQi6/UmlJWVsnOnZEgc9197Np1gD17DrJ7d0mDqaiotMmXdiYnJ9GjRwY9emTQs2cmPXtmRj/XTpn07JlBjx6188iUkZGq0Gtj1lrKdu1i79q1FK1bF5nWrmXvunWUbN9et19qVhZ9Jk2i96RJ9Jk8mb6TJ5M9duxRO+h3BAovjym8OgZrLWWFheyKhtmezz6jaP16itavp/LAgbr9fKmp9Bg+nJ4jRtBtyBC6DR5M9+i82+DBpPXq1SahEA5XU1RUWhdme/aUsGfPQYqLyyguLmXfvsg8slzGvn1llJdXHvF4Pl8SWVlpdO3ahW7d0snKSiMrK63uc/11keU0srLS6dYtjczMLmRmdiE9PZX09JROHYKVpaUc2LqV/Vu21E31l2uv8AGSMzLIHj2aXqNH02vMGHqNHk2fSZPoceKJcRmOqb1ReHlM4dWx1f7ruXjDBorWr6d4w4bItHEjB7Ztq7vNUyula9e6IOvqOGT260fXfv3I7NePzBNOqJsnd2n751oVFVUNQq328/795ZSUVHDgQO28goMHGy6XlJQTDtcc8xzGGDIzU8nISK0LtNpwi6xLJSOjS932L9Z3IS0tpRlTMkke/GK3NTWU7d7Nwfx8Drpu3VTqupTk59fNy3fvbvC95PR0ug8dSrchQ+g+dCi9Ro6sC6qs/v07ZUgdicLLYwqvzstaS0VREfu3buXAtm2Rqd7ngwUFlO3aBU38GU/r1YvMfv1Iz8mJTNnZpGdnkxad1/+clp1NWo8eMTX7b43/tvLyKkpKyikpOURJSXldqJWVVVJaWklZ2aG6eWTd4cv15xUVVS2qpUuX5Lowq//5iym17nN6egopKT5SU5NJTU2mS5dkkpMsKdWHSKmuwFdZRlL5QZLKD0DpAezBfdiD+6k+sI/w/n2E9xdTWbyXyuIibHV1w0KMIaNPH7o6Tt1UG1K1U0afPp36SjQWCi+PKbzkaGrCYcp27aJ0504OFhRQGp1qP5fv2UP53r1U7N1L+d69h//CrCc5PZ3Ubt3o0q0bqVlZpHbr9sVyvXUpGRkkp6d/MaWlNVxutM0kJ+NLSSEpOblNrwyqq2vqhdwhKiqqKC+toPxgOeWl5VSUllNRVkHFwXLKS0o5VHKQytLS6FROVVkZ4fIyqivKqS4vp+ZQBTWHKrCHKjCHykiqLMdXVU5yOBJUKTWH6FJTSQrhI9ZkgfKkdMpNWmRebzro60qpL5ODSZkc9HWlIjmT5C5d6NIlmZQUHykpPpKTa+dJ9T5Hlpv6fPi6ht9ruO7wbT5fw+/5fEn4fEkkJyc1+JyU1HCdz2fw+Xz4fIbkZB9JSaaujtrPPp+J61WuwisGjuOfDTwE+IB5rht6oNF2E91+IVAGfNN1Q58e7ZgKL2kt1loqDxygPBpkFfXmh/bv59CBA1SVlHDowAEqo1PjzzVVLbu6qWMMvpQUTHIySbVTNNhqp8Y1N1px2HZbXU1NOExNVRU1VVVUR+c14XCTV6UxlZuUREpmJl26d48EeffupHbvXhfoKd264cvoSlJGV5Iyo/OuPTBZ3TFdu1OTmkFVNVRVhTl0KExlZZiqqupjfq6sDBMO1xAOV1NVVR2df7FcXV1Tb3014XDD5Ybbv/ie14wxjQKufgAeHpTXXz+T2247o6XnSujwitt9EMfx+4BHgHOBPGCJ4/gXuG5oTb3dLgBGRKfpwB+jc5E2Z4yhS/fudOnevVmDDzclfOgQ4bIywuXlhMvLqYrOa6fqiooG66srKiLBUjtFQ6UmHMaGw1RXVWGjy7WfaXRb7LDbZI23+3x1IehLSYmEYW0gNrHOl5oauSqMXkHWXUlmZNR9rp0npXScBiPWWqqraxqF4pFDL7JfDTU1kXl19eFTZLslHK6mutpSXR2ZR5Zr6j433KepY9Q0eb7evbO8/rF5Jp7DDUwDNrluaDOA4/ifBy4D6ofXZcCfXTdkgY8cx9/Dcfz9XDdUEMc6RVosuUuXSGOQnj29LkViZEztrTwfkOJ1OXIM8Wx60x/YUW85L7ou1n0wxtxqjMk1xuSGw0e+ny4iIh1TPMOrqXsLjW+4N2cfrLVzrbVTrbVTk+PYAkxERNqHeIZXHjCw3vIAwG3BPiIi0snF87JlCTDCcfxDgXzgGuC6RvssAG6PPg+bDuzX8y4REWksblderhsKA7cDbwBrgRdcN7TacfxzHMc/J7rba8BmYBPwGPCdeNUnIiKJQ52URUQ6oUTv56WBvkREJOEovEREJOEk/G1DY0wNUO51HccpGY4yAFzno59HQ/p5fEE/i4aO5+eRbq1N2AuYhA+vjsAYk2utnep1He2Ffh4N6efxBf0sGurMP4+ETV0REem8FF4iIpJwFF7tw1yvC2hn9PNoSD+PL+hn0VCn/XnomZeIiCQcXXmJiEjCUXiJiEjCUXh5yBgz0BjzrjFmrTFmtTHmDq9r8poxxmeMWWaMedXrWrxmjOlhjHnJGLMu+mdkptc1eckY44/+PVlljHnOGJPmdU3xZIx5whizyxizqt66XsaYt4wxG6PzTvMWVIWXt8LAD621Y4AZwHeNMWM9rslrdxAZuFngIeBf1trRwEQ68c/FGNMf+D4w1Vp7EuAj8maKzuQpYHajdXcB71hrRwDvRJc7BYWXh6y1BdbaT6OfS4j8cjrszdGdhTFmAHARMM/rWrxmjOkGfBl4HMBaW2mt3edpUd5LBtKNMclABp3sXX/W2kVAUaPVlwFPRz8/DXwlnjV5SeHVThhjhgCTgY89LsVLDwI/Bmo8rqM9OBHYDTwZvY06zxiTsCOAHy9rbT4QBLYDBcB+a+2b3lbVLvS11hZA5B/DQB+P64kbhVc7YIzpCvwd+IG19oDX9XjBGHMxsMtau9TrWtqJZOBk4I/W2slAKZ3ollBj0Wc5lwFDAQfINMZ83duqxEsKL48ZY1KIBNdfrbUve12Ph04BLjXGbAWeB84yxjzjbUmeygPyrLW1V+IvEQmzzuocYIu1dre1tgp4GZjlcU3tQaExph9AdL7L43riRuHlIWOMIfJMY6219rde1+Mla+1/W2sHWGuHEHkQv9Ba22n/ZW2t3QnsMMaMiq46G1jjYUle2w7MMMZkRP/enE0nbsBSzwLgxujnG4H5HtYSV8leF9DJnQJ8A1hpjFkeXff/W2tf864kaUe+B/zVGJMKbAZu8rgez1hrPzbGvAR8SqSV7jI62dBIxpjngDOAHGNMHvBz4AHgBWPMzUQC/krvKowvDQ8lIiIJR7cNRUQk4Si8REQk4Si8REQk4Si8REQk4Si8REQk4Si8REQk4Si8REQk4Si8RNqAMeYcY8xfvK5DpKNSeIm0jYlERoEQkTag8BJpGxOBZcaYLsaYp4wx90XH5BORVqCxDUXaxkQiI3y/Acyz1nbmEfJFWp3GNhRpZdHX3OwBtgG3WWs/9LgkkQ5Htw1FWt9YYAmR0c+rPa5FpENSeIm0vonAYiLvJXvSGNPX43pEOhyFl0jrmwisstZuAO4k8r6lFI9rEulQ9MxLREQSjq68REQk4Si8REQk4Si8REQk4Si8REQk4Si8REQk4Si8REQk4Si8REQk4fw/0rf4PUhfgh8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "D.k = np.linspace(1,10,1000)\n",
    "D.tau_plot(D.tau_h_ana(),D.tau_h_dis(),ylabel=(r'$\\tau_{h}$ analytical',r'$\\tau_{h}$ discrete'))\n",
    "\n",
    "D.tau_plot(D.tau_b_ana(),D.tau_b_dis(),ylabel=(r'$\\tau_{b}$ analytical',r'$\\tau_{b}$ discrete'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The discrete solutions show smaller damping time scales for both equations, with a way stronger difference for the biharmonic case (Mind the scales!). However, the analytic solutions do not show any periodicity. Therefore the discrete solution can still have larger damping time scales for very short waves.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.)\n",
    "In the limit of $\\Delta x \\rightarrow 0$, the denominator and numerator get zero for both, harmonic and biharmonic damping time scale. Therefore, L'Hôpitals rule needs to be considered.\n",
    "\n",
    "$$\n",
    "\\lim_{\\Delta x \\rightarrow 0} \\frac{f(\\Delta x)}{g(\\Delta x)} = \\lim_{\\Delta x \\rightarrow 0} \\frac{f'(\\Delta x)}{g'(\\Delta x)}\n",
    "$$\n",
    "\n",
    "Applying this twice for the harmonic and four times for the biharmonic then yields:\n",
    "\n",
    "\\begin{align}\n",
    "\\lim_{\\Delta x \\rightarrow 0} \\frac{\\Delta x^2}{ 2 A_{2}\\left(1- \\cos(k \\Delta x) \\right)} &= \\lim_{\\Delta x \\rightarrow 0} \\frac{2}{ 2 A_{2} k^2 \\cos(k \\Delta x)} &&= \\frac{1}{A_{2} k^2} \\\\\n",
    "\\lim_{\\Delta x \\rightarrow 0} \\frac{\\Delta x^4}{A_{4} 2\\left( \\cos(2k \\Delta x) -4\\cos(2k \\Delta x)+3 \\right) } &= \\lim_{\\Delta x \\rightarrow 0} \\frac{24}{ 2 A_{4} k^4 (16\\cos(2k \\Delta x)-4\\cos(k \\Delta x))} &&= \\frac{1}{A_{4} k^4}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Diffusion:\n",
    "    def __init__(self, params=None):\n",
    "        if params is None:\n",
    "            params  = {'N': 21, 'L': 10, 'A_2': 1., 'A_4': 1., 'k_lim':(0.01,1.)}\n",
    "            \n",
    "        self.N     = params['N']\n",
    "        self.L     = params['L']\n",
    "        self.A_2   = params['A_2']\n",
    "        self.A_4   = params['A_4'] \n",
    "        self.k_lim = params['k_lim']\n",
    "        self.update_params()\n",
    "\n",
    "        \n",
    "    def A_harmonic(self):\n",
    "        #harmonic matrix:\n",
    "        A = np.zeros((self.N-2, self.N-2))\n",
    "        indices = np.arange(self.N -2)\n",
    "\n",
    "        # set diagonal\n",
    "        A[indices, indices] = -2 / self.dx ** 2\n",
    "\n",
    "        # set off-diagonal\n",
    "        A[indices[:-1], indices[:-1] + 1] = 1 / self.dx ** 2\n",
    "        A[indices[:-1] + 1, indices[:-1]] = 1 / self.dx ** 2\n",
    "        return A\n",
    "    \n",
    "    def A_biharmonic(self):\n",
    "        #harmonic matrix:\n",
    "        A = self.A_h * self.A_h\n",
    "        return A\n",
    "    \n",
    "    def display(self):\n",
    "        print('N:      ', self.N)\n",
    "        print('L:      ', self.L)\n",
    "        print('A_2:    ', self.A_2)\n",
    "        print('A_4:    ', self.A_4)\n",
    "        print('dx:     ', self.dx)\n",
    "        print('x:      ', self.x)\n",
    "        \n",
    "    def tau_h_dis(self):\n",
    "        return(self.dx**2 / (2*self.A_2*(1-np.cos(self.k*self.dx))))\n",
    "    \n",
    "    def tau_b_dis(self):\n",
    "        return((self.dx**4)/(2*self.A_4*(np.cos(2*self.k*self.dx)-4*np.cos(self.k*self.dx)+3)))\n",
    "    \n",
    "    def tau_h_ana(self):\n",
    "        return(1/(self.A_2*self.k**2))\n",
    "    \n",
    "    def tau_b_ana(self):\n",
    "        return(1/(self.A_4*self.k**4))\n",
    "    \n",
    "    def update_params(self):\n",
    "        self.dx    = self.L / (self.N -1)\n",
    "        self.x     = self.dx*np.arange(self.N)\n",
    "        self.A_h   = self.A_harmonic()\n",
    "        self.A_b   = self.A_biharmonic()\n",
    "        self.k   = np.linspace(self.k_lim[0],self.k_lim[1],1000)\n",
    "        \n",
    "    def tau_plot(self, tau_h, tau_b, ylabel=(r'$\\tau$ harmonic',r'$\\tau$ biharmonic')):\n",
    "        fig, ax1 = plt.subplots()\n",
    "        ax2 = ax1.twinx()\n",
    "        ax1.plot(self.k,tau_h,color='midnightblue',label='harmonic')\n",
    "        ax2.plot(self.k,tau_b,color='darkred',label='biharmonic')\n",
    "        ax1.set_ylabel(ylabel[0],color='midnightblue')\n",
    "        ax2.set_ylabel(ylabel[1],color='darkred')\n",
    "        ax1.set_xlabel(r'$k$')\n",
    "        ax1.tick_params(axis='y', labelcolor='midnightblue')\n",
    "        ax2.tick_params(axis='y', labelcolor='darkred')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ws20climnum]",
   "language": "python",
   "name": "conda-env-ws20climnum-py"
  },
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}