03_Generalized_Midpoint_Picard.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": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3 - Generalized midpoint rule, Picard iterations\n",
    "\n",
    "The ODE $\\dot{y} = k (y_\\text{eq} - y)$ is to be integrated numerically.\n",
    "\n",
    "First, we import some libraries for numerical basics and for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Some plot settings\n",
    "%run plot_functions/plot_settings.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first code the analytical solution of the ODE in order to have something to compare against.\n",
    "\n",
    "$\\displaystyle \\frac{\\text{d}y}{y_\\text{eq}-y} = k\\,\\text{d}t \\quad \\rightarrow \\quad \\int \\limits_0^t k\\,\\text{d}\\bar{t} = \\int \\limits_{y_0}^y \\displaystyle \\frac{\\text{d}\\bar{y}}{y_\\text{eq}-\\bar{y}}$\n",
    "\n",
    "$kt = \\displaystyle \\left. - \\ln (y_\\text{eq}-\\bar{y}) \\right|^y_{y_0}$\n",
    "\n",
    "$\\exp (-kt) = \\displaystyle \\frac{y_\\text{eq}-y}{y_\\text{eq}-y_0}$\n",
    "\n",
    "$y = y_\\text{eq} - (y_\\text{eq} - y_0) \\exp (-kt)$\n",
    "\n",
    "Let's define the equilibrium value and the rate constant as global variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_eq = 1.\n",
    "k = 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def y_analytical(t,y_0):\n",
    "    return y_eq - (y_eq - y_0)*np.exp(-k*t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the ODE definition $\\dot{y}=f[y(t),t]$ we define the rate function itself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dy_dt(y):\n",
    "    return k*(y_eq - y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the function and its rate for some exemplary values ($y_0 = 0$) in a time interval $t \\in [0,4]$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = np.linspace(0,4,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$y$, $\\\\dot{y}$ / s$^{-1}$')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAIZCAYAAABd4AK4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABoP0lEQVR4nO3dd3xUZf7+/2vSGxkTkhBKRDqIS0ABkY4CKwqI8hUVVkHXjm1VVNxV0I+/ZVVcFXERVyliRXFXBGwQFJcoRaSIqEAoQVogIb3P+f0xaUMSkpCZnCmv52PnMWdOuec9x2H24nCf+7YYhmEIAAAAQK38zC4AAAAAcHeEZgAAAKAOhGYAAACgDoRmAAAAoA6EZgAAAKAOhGYAAACgDoRmAAAAoA6EZgAAAKAOhGYXMgxDWVlZYv4YAAAAz0ZodqHs7GxZrVZlZ2ebXQoAAAAagdAMAAAA1IHQDAAAANSB0AwAAADUgdAMAAAA1IHQDAAAANSB0AwAAADUgdAMAAAA1IHQDAAAANQhwOwCAAAAfJVhGCouLpbNZjO7FK/g5+enwMBAWSwWp7dNaAYAAGhiRUVFOn78uPLy8lRaWmp2OV7F399fYWFhiouLU1BQkNPaJTQDAAA0oby8PKWmpsrf319RUVEKDQ2Vv7+/S66O+hLDMFRaWqr8/HxlZmZq//79atOmjcLCwpzSPqEZAACgCZ04cUKBgYFq27at/P39zS7H60RERCg6OloHDhzQiRMndO655zqlXa+7EfDgwYOKjIyUxWLR119/3eDjCwoK9PTTT6tLly4KCQlRy5YtdeONN2rfvn3OLxYAAPiUkpIS5ebmKjo6msDsQv7+/oqOjlZubq5KSkqc0qbXheZbb71V2dnZZ3Vsfn6+LrvsMj3//PN65plnlJ6ers8//1w7duxQr169tG3bNidXCwAAfEl5gAsODja5Eu9Xfo4JzTV4/fXXtXnzZnXt2vWsjp85c6aSk5P1j3/8Q9dee63CwsKUmJioZcuWKScnR5MmTaKzPgAAaDT6L7ues8+x14TmgwcPatq0aXrxxRfVokWLBh+fl5enV199VUFBQZo8ebLDtg4dOmj48OHauXOnVq1a5aySAQAA4CG8JjTfeuutGjBgQLXAW19JSUnKzc1VYmKiIiIiqm0fOHCgJGn58uWNqhMAAACexytGz/j3v/+tDRs26KeffjrrNsr7K7dr167G7e3bt3fYDwAAAL7D40NzamqqHn74Yc2ePVsJCQln3c7Ro0clSdHR0TVuj4qKctjP7Rz9Sfp+npR7XPrDBKnHtWZXBAAA4DU8PjTfdttt6tOnj2677bZGtZOXlydJtc4cU34HZvl+NSksLFRhYWHF66ysrEbV1CB5J6Wtb9uX47pJIjQDAAA4i0eH5jfeeEP/+9//GtUto1z5bDFFRUU1bi8Pw2eaVWbWrFl66qmnGl3LWYmIq1zOSTOnBgAAAC/lsTcCHjp0SA899JD+8Y9/6Lzzzmt0e/Hx8ZKk9PT0GrdnZGRIklq2bFlrG9OnT1dmZmbFIzU1tdF11VtElRFDco833fsCAAA00ptvvimLxaK77rqr1n2uvfZaWSwWffrpp01YWSWPDc2rV69WVlaW7r33XlksFofHN998I0kaNmxYxbq6Zgfs2bOnJNU6819KSookqUePHrW2ERwcrMjISIdHkwk5R/Ir+4eDHEIzAADwHP3795ckbdy4scbtSUlJ+uijj3T55ZdrzJgxTVlaBY8NzVOmTJFhGDU+hgwZIklau3ZtxbqhQ4eesb1hw4YpPDxc27dvV25ubrXt69evlySNHTvW6Z/FKfz8pPBY+zKhGQAAeJCuXbsqOjpa27dvV35+vsO2kpIS3XfffQoMDNRLL71kToHy8D7NZ2Pt2rV69NFHdeONN+ree++tWB8WFqZ77rlHzz77rBYvXqy77767YltKSopWr16t7t2764orrjCj7PoJj5Wyj0i5aZLNZg/SAADAo4x55X9Kyy6se0c3ENssWJ/eO7DR7VgsFvXr10+rVq3Sli1bNGDAgIptc+fO1c6dO/Xwww+rS5cujX6vs+VzofmFF17Qpk2btGvXLofQLEkzZszQunXr9Nhjj6lFixa64oortHv3bk2ePFlhYWF655135O/vb1Ll9VDer9kolfIzpPDm5tYDAAAaLC27UEezCswuo8n1799fq1at0oYNGypC8/HjxzVz5kzFx8frySefNLU+rwnN+/fvrzYxybBhwyRJQ4YMqejTfN1112ndunW6+eabq7URGhqqpKQkPfvss5o+fbomTZokq9WqESNGaNmyZRUTnLitqiNo5B4nNAMA4IFimwWbXUK9ObPW8n7NGzZsqFj32GOPKTMzUy+//LKaNWvmtPc6GxbDMAxTK/BiWVlZslqtyszMbJqbAr+aIa1/yb580ydS+6Guf08AAFBvBQUF2rdvn9q1a6eQkBCzy3Erubm5Ouecc9S6dWvt379fGzduVL9+/XTxxRcrOTlZFoulQe05+1zT6dWbMFYzAADwUOHh4erRo4cOHDigI0eO6J577pEkzZkzp8GB2RUIzd6EsZoBAIAHK++icdddd2nTpk26+eab1adPH5OrsiM0e5PyIeckhp0DAAAepzw0f/LJJ7JarZo1a5bJFVUiNHsTh+4ZhGYAAOBZykOzZB/VLC4u7gx7Ny1CszcJP230DAAAAA+Sl5cnSerWrVu1oYHNRmj2JqFRTKUNAAA81rRp0yRJL7/8sgIC3GtkZEKzN6k6lXYuo2cAAADPsWDBAq1cuVLXXXedRowYYXY51bhXhEfjlU+lnXOcqbQBAIBb27Vrl1544QUdPHhQq1ev1rnnnqu5c+eaXVaNSFTepvxmwPKptAEAANxUUlKS3nzzTW3cuFFjxoxRUlKSYmJizC6rRlxp9jan3wzIVNoAAMBNTZ06VVOnTjW7jHrhSrO3Ydg5AAAApyM0extCMwAAgNMRmr0NYzUDAAA4HaHZ20QwlTYAAICzEZq9TUSLymXGagYAAHAKQrO3qdo9I+eYeXUAAAB4EUKztwmNkiz+9mW6ZwAAADgFodnbMJU2AACA0xGavVH5sHO5afaptAEAANAohGZvVB6abSVMpQ0AAOAEhGZvxFjNAAAATkVo9kaM1QwAAOBUhGZvxFjNAADAwxUXF+uXX37RwYMHzS5FEqHZOzFWMwAA8HC///67unXrpptuusnsUiQRmr0T3TMAAACcKsDsAuACDjcC0j0DAAB4nvPOO0+GYZhdRgWuNHujqn2audIMAADQaIRmb+QwlTZ9mgEAABqL0OyNmEobAADAqQjN3iqiSmhmKm0AAODG3nzzTVksFt1111217nPttdfKYrHo008/bcLKKhGavVV5v2ZbiVRwytRSAAAAzqR///6SpI0bN9a4PSkpSR999JEuv/xyjRkzpilLq0Bo9laM1QwAADxE165dFR0dre3btys/P99hW0lJie677z4FBgbqpZdeMqdAEZq9F2M1AwAAD2GxWNSvXz+VlJRoy5YtDtvmzp2rnTt36v7771eXLl1MqpBxmr0XYzUDAOCZ5g/xnAteEXHSHd84pan+/ftr1apV2rBhgwYMGCBJOn78uGbOnKn4+Hg9+eSTTnmfs0Vo9laM1QwAgGfKOS5lHza7iiZX3q95w4YNFesee+wxZWZm6uWXX1azZs3MKk0Sodl7OXTPoE8zAAAeIyKu7n3chRNr7du3rwICAipC88aNG7Vo0SL169dPN910k9Pe52wRmr0V3TMAAPBMTuru4GnCw8PVo0cPbdmyRUeOHNE999wjSZozZ44sFovJ1XnBjYCbN2/Wk08+qYEDB6p169YKCQlR27ZtdeWVV2rFihUNamvmzJmyWCy1PqZMmeKaD+EKVf/mR/cMAADgAcq7aNx1113atGmTbr75ZvXp08fkquw8PjSPHj1azz33nG688Ub9+OOPOnnypN5++22lpqZqzJgxDe40HhkZqS5dutT4aNmypYs+hQuERldOpZ1LaAYAAO6vPDR/8sknslqtmjVrlskVVfKK7hmPP/647rjjjorXgwYN0sKFC9W7d28999xzmj59ukJDQ+vV1tVXX61Fixa5qNImVD6Vds5RrjQDAACPUB6aJWnGjBmKi3Of/t0eH5o//PBDnX/++dXWl4/jV1hYqIKCgnqHZq8SURaay6fS9vP4f1gAAABeLC8vT5LUrVs33XvvvSZX48jjQ/OgQYNqXJ+cnCxJ6t27t6KiopqyJPdRfjNg+VTaYdGmlgMAAHAm06ZNkyS9/PLLCghwr5jqVZceS0tLdfjwYS1ZskRTpkzRRRddpPfff79Bbezbt08TJ05UQkKCgoKC1KJFC40dO1Zr1qxxUdUuxM2AAADAQyxYsEArV67UddddpxEjRphdTjVeFZqDg4PVunVr3XzzzRo3bpxWrVqlDh06NKiNjRs3qnfv3tq4caNycnK0cuVKZWRkaPjw4XrmmWfOeGxhYaGysrIcHqZyCM2M1QwAANzLrl27dOutt2rkyJG69dZbde6552ru3Llml1UjrwrNJSUlSktL02effaYNGzaoY8eOWrJkSb2PnzhxojZt2qQHH3xQLVu2VFBQkHr37q0VK1YoKipKTzzxhNavX1/r8bNmzZLVaq14JCQkOONjnT3GagYAAG4sKSlJb775pjZu3KgxY8YoKSlJMTExZpdVI68KzZIUExOjESNGaO3atYqIiNDkyZO1du3aeh3buXNnXXDBBdXWW61WTZo0SZK0cOHCWo+fPn26MjMzKx6pqaln9yGche4ZAADAjU2dOlWGYejUqVP65JNPGtxDoCl5XWguFxkZqUmTJskwDM2bN6/R7XXs2FGS/Z8RahMcHKzIyEiHh6nCq0ylzVjNAAAAZ81rQ7MktWvXTpK0f//+Rrdls9kkyS2mcay3iBaVy1xpBgAAOGseHZqXLFmimJgYGYZR4/bDhw9Lkpo3b15nW6mpqYqPj9f27dtr3L57925JUteuXc+yWhPQPQMAAMApPDo0l5aW6uTJk/r666+rbcvPz68Ybm7s2LEV63fs2KEBAwbo8ccfr9bWsWPHtGrVqmptZWRk6N1335UkTZkyxXkfwNWYShsAAMApPDo0l3eVmDhxot566y0dOXJEOTk5+u677zRq1Cjt3btXo0aN0m233VZxzPz585WcnKxZs2bp5MmT1dp66qmn9NJLL+nQoUMqKirSpk2bNHr0aGVmZmrmzJkaOHBg037IxvDzk8LL7kDNYfQMAACAs+XRoXnSpEn67LPPNGrUKD3zzDPq0KGDoqKidNVVVykoKEiLFi3SihUrHGaUGTdunKKiojRu3DhFR1fOkNe2bVtt3rxZ9913nxYvXqxu3bopPDxcY8eOVVxcnL766ivNmDHDjI/ZOOVdNHKP26fSBgAAQINZjNo6BKPRsrKyZLValZmZad5IGkuukfaWzWb4yD6m0gYAwEQFBQXat2+fzjvvPIWGhppdjlfLz8/X/v371a5dO4WEhDS6PY++0ox64GZAAADcRvm/fhcWFppcifcrP8dVexw0BqHZ2zmE5qPm1QEAABQQEKDw8HClp6ertLTU7HK8VmlpqdLT0xUeHu600OycVuC+mrWqXM46Yl4dAABAkn324tTUVO3bt09Wq1WhoaHy9/f3rLkg3JBhGCotLVV+fr4yMzNls9nUsmVLp7VPaPZ2kVW+LNmHzasDAABIksLCwtSuXTsdP35cGRkZOnHihNkleRV/f3+FhYUpLi5OQUFBTmuX0OztuNIMAIDbCQoKUps2bWQYhoqLiytmHkbj+Pn5KTAw0CVX7QnN3i6yamjmSjMAAO7EYrE49WooXIcbAb1dRAvJUvafme4ZAAAAZ4XQ7O38A6TwshE06J4BAABwVgjNvqC8i0bOMam02NxaAAAAPBCh2RdU9Gs27MEZAAAADUJo9gXNqgw7RxcNAACABiM0+wKHETR+N68OAAAAD0Vo9gVVQ3M2V5oBAAAaitDsCxy6ZzDsHAAAQEMRmn1BZOvKZUIzAABAgxGafUFklSvNdM8AAABoMEKzLwgKl4Kt9mWuNAMAADQYodlXlN8MmHVYMgxzawEAAPAwhGZfUd5Fo7RQys8wtxYAAAAPQ2j2Fc2qjtVMFw0AAICGIDT7ikhCMwAAwNkiNPsKhxE0CM0AAAANQWj2FQ5jNTPsHAAAQEMQmn2Fw6yAv5tXBwAAgAciNPuKqn2ameAEAACgQQjNviKsueQfZF+mewYAAECDEJp9hcVS2UWD7hkAAAANQmj2JeVdNApOScX5ppYCAADgSQjNvoSxmgEAAM4KodmXOIygQWgGAACoL0KzL2EEDQAAgLNCaPYldM8AAAA4K4RmX9KM0AwAAHA2CM2+JLJKn+ZsQjMAAEB9EZp9SUS8JIt9mQlOAAAA6o3Q7EsCgqTwWPsy3TMAAADqzaND8+bNm/Xkk09q4MCBat26tUJCQtS2bVtdeeWVWrFiRYPbs9lsmjNnjhITExUWFqaYmBhdffXV2rp1q/OLN0t5F42cY5Kt1NxaAAAAPIRHh+bRo0frueee04033qgff/xRJ0+e1Ntvv63U1FSNGTNGTz75ZL3bstlsmjBhgh566CHdeeedSktL0/fff6/s7GxdfPHF+vLLL134SZpQZGv7s1Eq5Rw3txYAAAAPYTEMwzC7iLMVHx+vu+++u1o4/uGHH9S7d28FBwcrIyNDoaGhdbY1b9483X333Xr44Yf1/PPPV6zPyMhQhw4dFBgYqD179qhZs2b1ri8rK0tWq1WZmZmKjIys/wdzpRUPSpvftC/fmiS1ucjcegAAADyAR19p/vDDDzV16tRq67t06SJJKiwsVEFBQb3amj17tiTp9ttvd1gfFRWlCRMm6Pjx41q0aFHjCnYHjKABAADQYB4dmgcNGqTmzZtXW5+cnCxJ6t27t6Kioups56efflJKSopiY2PVqVOnatsHDhwoSVq+fHkjK3YD5d0zJEbQAAAAqKcAswtwltLSUh07dkxr1qzRo48+qosuukjvv/9+vY7dtm2bJKldu3Y1bm/fvr3Dfh6tWZUrzVm/m1cHAACAB/Ga0BwcHKzS0lL5+/vr9ttv18yZMxUXF1evY48ePSpJio6OrnF7+dXqtLQ02Ww2+fl58AX6qleas7nSDAAAUB8enP4clZSUKC0tTZ999pk2bNigjh07asmSJfU6Ni8vT5IUFBRU4/bg4OBq+9aksLBQWVlZDg+3U7VPM2M1AwAA1IvXhGZJiomJ0YgRI7R27VpFRERo8uTJWrt2bZ3HhYWFSZKKiopq3F5YWFht35rMmjVLVqu14pGQkNDAT9AEgptJQWUjgBCaAQAA6sWrQnO5yMhITZo0SYZhaN68eXXuHx8fL0lKT0+vcXtGRoYkKS4u7oxdM6ZPn67MzMyKR2pq6llU3wQiW9mfs49InjviIAAAQJPxmj7Npyu/qW///v117tuzZ09J0r59+2rcnpKSIknq0aPHGdsJDg526MrhtiJbSid+lYrzpIJMKfQcsysCAABwax57pXnJkiWKiYlRbXOzHD5s73pQ05B0p+vevbvat2+vtLQ07dmzp9r29evXS5LGjh3biIrdSLNWlct00QAAAKiTx4bm0tJSnTx5Ul9//XW1bfn5+RXDzVUNujt27NCAAQP0+OOPVztm2rRpkqTXX3/dYf2pU6e0dOlSxcXFacqUKc77AGaKrBKameAEAACgTh4bmi0WiyRp4sSJeuutt3TkyBHl5OTou+++06hRo7R3716NGjVKt912W8Ux8+fPV3JysmbNmqWTJ086tHf77bfrmmuu0Ysvvqj58+crNzdXe/fu1YQJE5Sdna0lS5Y0aAptt+YwggbDzgEAANTFY0PzpEmT9Nlnn2nUqFF65pln1KFDB0VFRemqq65SUFCQFi1apBUrViggoLLb9rhx4xQVFaVx48ZVG5PZz89PS5cu1ezZs/Xqq68qJiZGffv2VVhYmDZs2KCRI0c29Ud0HbpnAAAANIjFqK1TMBotKytLVqtVmZmZioyMNLucSoe3Sq8PsS9fNEUa87KZ1QAAALg9j73SjEao2qeZ7hkAAAB1IjT7orAYyS/Qvkz3DAAAgDoRmn2Rn5/UrOxmwKzfza0FAADAAxCafZW1jf05P10qzDG3FgAAADdHaPZV55xbuZzpptN9AwAAuAlCs6+Kalu5fOqgeXUAAAB4AEKzr6p6pZnQDAAAcEaEZl/lEJoPmFcHAACAByA0+yquNAMAANQbodlXRbaWLGX/+TO40gwAAHAmhGZf5R9oD84SV5oBAADqQGj2ZeVdNPLTpcJsc2sBAABwY4RmX+bQr5mxmgEAAGpDaPZl3AwIAABQL4RmX3YOE5wAAADUB6HZlzFWMwAAQL0Qmn0Z3TMAAADqhdDsyyJbSxZ/+zJXmgEAAGpFaPZl/gGM1QwAAFAPhGZfVzFWc4ZUkGVuLQAAAG6K0OzrqvZrzmSsZgAAgJoQmn0dNwMCAADUidDs66IYqxkAAKAuhGZfx5VmAACAOhGafV3V0Jyx37QyAAAA3Bmh2dc1a1VlrGauNAMAANSE0Ozr/AMkK2M1AwAAnAmhGdI5ZTcDFpySCjJNLQUAAMAdEZpx2s2AjNUMAABwOkIzGEEDAACgDoRmEJoBAADqQGhGZZ9midAMAABQA0IzTrvSfMC8OgAAANwUoRlSs5aSX4B9mdAMAABQDaEZ9rGaIxmrGQAAoDaEZtiVd9EoyJTyT5laCgAAgLshNMOu6s2AmYzVDAAAUJVHh+akpCT9+c9/VufOnRUSEqKwsDCdf/75mjZtmtLS0hrU1pQpU2SxWGp9zJw50zUfwl0w7BwAAECtAswu4Gy9+uqruueee9SjRw+99tpr6tOnj7KysvTOO+/ob3/7m5YsWaJ169apc+fO9W6zefPmiomJqXFbbeu9BqEZAACgVh4bmvPz8xUUFKQVK1YoISFBktSsWTM98sgjSk9P17PPPqv7779fn332Wb3bvOeee7z/inJtoqp0z8hgBA0AAICqPLZ7Rnx8vG644YaKwFzV2LFjJUmrV69WaWlpU5fmmbjSDAAAUCuPvdL8pz/9SX/6059q3Ga1WiVJFotFhmE0ZVmeq3ysZlsJoRkAAOA0Hnul+Ux++eUXSdLAgQMVEFD/vxds3bpVY8eOVatWrRQcHKzWrVvr+uuv1+bNm11Vqvvw85esbezLhGYAAAAHXhma33rrLVksFj355JMNOu7bb7/VVVddpR07digzM1Nvv/22fvrpJ11yySVavHhxnccXFhYqKyvL4eFRyrtoFDJWMwAAQFUWw8v6L3z++ecaNWqUHnroIc2ePbvex23btk0RERHq0KGDw/q9e/eqa9eu8vPz0y+//KJ27drV2sbMmTP11FNPVVufmZmpyMjI+n8Is3wyVfrxbfvyHd9KLXuYWw8AAICb8Korzbt379bkyZM1fvx4Pfvssw06NjExsVpglqQOHTro8ssvV1FRkd55550ztjF9+nRlZmZWPFJTPWySkKoTnNBFAwAAoILH3gh4ugMHDmjEiBEaOHCg3n33Xfn7+zut7Y4dO0qSdu3adcb9goODFRwc7LT3bXIOI2gw7BwAAEA5r7jSvGfPHg0ePFiDBg3S0qVLFRQU5NT2bTabJPtoHF6tamhmrGYAAIAKHh+af/75Zw0ePFgjR47U4sWLHa4wz549u15dJJKTkxUfH6+MjIwat+/evVuS1LVrV+cU7a6izqtczthvVhUAAABux6ND87Zt2zR06FCNHz9er7/+uvz8HD/OtGnTtHfv3orXa9euVd++ffXKK6847FdUVKRjx47piy++qPYev/32m7788ksFBgZq0qRJrvkg7qJZSykwzL58co+5tQAAALgRj+3TvHnzZo0cOVKFhYVKS0vTDTfcUOcxL7zwgjZt2qRdu3bp3nvvrVhf3u1i6tSpKigo0KhRo2S1WvXdd99p6tSpkqR58+adceQMr2CxSNEdpGM77FeaS4sl/0CzqwIAADCdx4bmt99+u6I7xQcffFCvY6677jqtW7dON998s8P6wYMHa926dfrggw/0/PPP6+6771ZJSYni4+M1aNAgvfXWW+rdu7fTP4Nbal4Wmo1S+wgazauPKAIAAOBrvG6cZneSlZUlq9XqOeM0S9Ka/5O+LRvfeuJSqfMfza0HAADADXh0n2a4QPOOlcv0awYAAJDkgtD8+eefq3379s5uFk2lancMQjMAAIAkF4Tm3NxcHTjAGL8ey+FK897a9wMAAPAh9b4R8Mknn6zXfnXNmgc3FxYthUZJ+RmEZgAAgDL1Ds3PPPOMzjnnHFmt1jPul5eX1+iiYLLmHaVDm6SsQ1JRnhQUZnZFAAAApqp3aO7QoYMGDRqkBQsWnHG/jz76SNddd12jC4OJojvYQ7MkpadI8ReYWw8AAIDJ6t2n+ZJLLlFycnKd+1ksFjGKnYdjBA0AAAAH9b7SfP/992v9+vV17jdkyBCtXbu2UUXBZFVH0EinXzMAAEC9Q/NFF12kiy66qM79YmJiNGTIkEYVBZMxggYAAIADJjdBddFVxtmmewYAAEDjQ3NqaqouvfRSZ9QCdxEcITVraV8mNAMAADQ+NOfl5embb75xRi1wJ+VdNPJO2sdsBgAA8GF0z0DNHKbTTjGvDgAAADdAaEbNGHYOAACgAqEZNYuueqWZ0AwAAHwboRk1q3qlmbGaAQCAjyM0o2ZR50mWsq8HV5oBAICPIzSjZgFB0jlt7csn90pMjQ4AAHyYU0KzQaDyTuUjaBTlSDnHzK0FAADARI0OzfHx8Zo3b54zaoG7YTptAAAASU4IzVarVXfccYczaoG7Ydg5AAAASfRpxpk0Z9g5AAAASQpoijf5/vvvtXTpUoWHh+vWW29V27Ztm+Jt0Vh0zwAAAJDk5ND8+eef68orr5Qk7d69W+3bt9c333yj4cOHy2azyTAMvfbaa/rhhx907rnnOvOt4QqRbST/YKm0kLGaAQCAT3Nq94xly5apU6dO+uKLLyquJk+fPl2lpaX629/+pmXLlqlVq1Z64YUXnPm2cBU/Pym6vX05PUWylZpbDwAAgEmceqX5+++/17///W8NGjRIkv1q8/fff68xY8boqaeekmQfbePmm2925tvClZp3kNJ2SaVFUmaqfdITAAAAH+PUK8379+9Xnz59Kl5/+umnslgsDiG5V69eOnTokDPfFq7ECBoAAADODc3h4eHKycmpeP3hhx8qNDRUl19+ecW6/Px8hYaGOvNt4UoOI2ikmFcHAACAiZwamrt376433nhDkv2mwA0bNmjs2LEKCQmp2Gft2rVq06aNM98WrsSVZgAAAOeG5vvvv1+PP/64YmNjNXr0aPn7++svf/mLJOnEiROaN2+e7rrrLg0ePNiZbwtXIjQDAAA4NzSPHTtWr732mtq2batevXrprbfequjj/OOPP+q5555TWFiYJk6c6My3hSuFx0rBkfZlQjMAAPBRFsMwDLOL8FZZWVmyWq3KzMxUZGSk2eWcvdeHSod/lCx+0l+PSgHBZlcEAADQpJhGG3WLLrsZ0LBJGftNLQUAAMAMhGbUrWq/5hO7zasDAADAJIRm1C22S+Xy8V3m1QEAAGASQjPqFnd+5fLxn82rAwAAwCQeHZqTkpL05z//WZ07d1ZISIjCwsJ0/vnna9q0aUpLS2twewUFBXr66afVpUsXhYSEqGXLlrrxxhu1b98+F1TvQZp3kPwC7ctcaQYAAD7IY0Pzq6++qssuu0ybN2/Wa6+9prS0NO3evVtTpkzRyy+/rD/84Q/67bff6t1efn6+LrvsMj3//PN65plnlJ6ers8//1w7duxQr169tG3bNhd+GjfnH1jZRePkbqmkyNx6AAAAmlijQnPHjh31yCOPKDk52Vn11Ft+fr6CgoK0YsUKXXrppWrWrJlat26tRx55RA8++KCOHTum+++/v97tzZw5U8nJyfrHP/6ha6+9VmFhYUpMTNSyZcuUk5OjSZMmqbS01IWfyM3FdbM/20rswRkAAMCHNCo0t27dWi+++KIGDRqkVq1a6e6779bq1aubJFzGx8frhhtuUEJCQrVtY8eOlaR615KXl6dXX31VQUFBmjx5ssO2Dh06aPjw4dq5c6dWrVrlnOI9UdV+zcfo1wwAAHxLo0LzN998oyNHjmj+/Pnq1auXFi5cqD/+8Y+KjY3VTTfdpP/+97/Kz893Vq0O/vSnP2nRokU1brNarZIki8Wi+szdkpSUpNzcXCUmJioiIqLa9oEDB0qSli9ffvYFezpuBgQAAD6s0X2aY2JidOutt2rlypVKS0vTO++8oxEjRui///2vrrnmGsXGxmr8+PF6++23derUKSeUXLdffvlFkj3sBgQE1Ll/eX/ldu3a1bi9ffv2Dvv5pPLuGRI3AwIAAJ/j1BsBIyIidP311+uDDz7QiRMntHz5ck2YMEHffvutbrrpJrVo0UIjR47Ua6+9pqNHjzrzrR289dZbslgsevLJJ+u1f3kt0dHRNW6Piopy2M8nnXOuFFR2Ff74TnNrAQAAaGIuGz0jKChIo0eP1oIFC3T06FElJSXpjjvu0K+//qq7775bbdq00YABA5z+vp9//rmWL1+uBx98UEOHDq3XMXl5eRU11yQ4ONhhv9oUFhYqKyvL4eE1LJbKq82nDkqF2ebWAwAA0ISaZMg5Pz8/DR06VHPmzNGBAwe0ceNGPfLII8rIyHDq++zevVuTJ0/W+PHj9eyzz9b7uLCwMElSUVHNQ6kVFhY67FebWbNmyWq1VjxquknRozl00fjFvDoAAACamCnjNPfu3Vt///vf9fPPzruh7MCBAxoxYoQGDhyod999V/7+/vU+Nj4+XpKUnp5e4/bycN+yZcsztjN9+nRlZmZWPFJTU+tdg0eI6165zM2AAADAh9R9l5wH2LNnjy677DINHjxYixYtalBglqSePXtKUq0z/6WkpEiSevToccZ2goODK7pyeCVuBgQAAD7KY2cELPfzzz9r8ODBGjlypBYvXuwQmGfPnl2vq73Dhg1TeHi4tm/frtzc3Grb169fL6ly/Gef1aLqlWZuBgQAAL7Do0Pztm3bNHToUI0fP16vv/66/PwcP860adO0d+/eitdr165V37599corrzjsFxYWpnvuuUeFhYVavHixw7aUlBStXr1a3bt31xVXXOG6D+MJwmOk8Fj7MleaAQCAD/HY7hmbN2/WyJEjVVhYqLS0NN1www11HvPCCy9o06ZN2rVrl+69916HbTNmzNC6dev02GOPqUWLFrriiisqbiwMCwvTO++80+BuH14prpu0L03KTZNy0qSIWLMrAgAAcDmPDc1vv/12xQ16H3zwQb2Oue6667Ru3TrdfPPN1baFhoYqKSlJzz77rKZPn65JkybJarVqxIgRWrZsWcUEJz4vrru0b519+fjPUsQQc+sBAABoAhajPvNM18PKlSt15ZVXOqMpr5GVlSWr1arMzExFRkaaXY5z/LBY+vQ++/Llz0r97jS3HgAAgCbgtD7NY8aM0ejRo7Vnzx5nNQl3xM2AAADABzktNI8fP16rVq3SBRdcoMcee6zGUSjgBWK7VC5zMyAAAPARTgvNH374oVavXq2OHTvqueeeU5cuXfT22287q3m4i+Bm0jnn2peP75Kc07sHAADArdU7NKempurSSy/VZZddVus+l156qbZt26Z//vOfys3N1eTJkzVw4EBt2bLFKcXCTZTPDFiUI506aG4tAAAATaDeoTkvL09ff/21vv766zPu5+/vrwceeEC//fabJk+erO+++059+/bV7bffrhMnTjS2XrgDZgYEAAA+pt6huUOHDtq3b1/FlNJ1iY2N1YIFC/Tdd9/pwgsv1BtvvKHOnTtrzpw5stlsZ10w3EDc+ZXL3AwIAAB8QL1Dc0BAgNq2bau2bds26A369u2r5cuX66qrrtKpU6f0l7/8RYmJiUpKSmpwsXATLaqGZq40AwAA7+f0yU0KCwu1ZcsWbdiwoeJx4MCBiu2GYWjnzp0aMWKErr/+er3yyiuKjo52dhlwpeadJL8AyVZCaAYAAD7BaaF56tSp2rhxo7Zv366SkhJJ9oAcEBCgXr16qX///hWPvLw8Pfzww3rvvff03Xffac2aNWrXrp2zSoGrBQTZg3PaLintV6m0WPIPNLsqAAAAl3HajIB+fvaeHtHR0erXr19FQO7bt6/CwsJqPGbBggW67bbbdNVVV+njjz92RhluxStnBCz34c3SzrL/ZndvkOK6mlsPAACACzntSvMbb7yh/v37q2vX+oenW265RR988EGdI3LADcWdXxmaj+8kNAMAAK/mtMlNbrnllgYF5nIJCQnKzMx0VhloKtwMCAAAfIjTbwRsqIcfflgXXnih2WWgoRirGQAA+BDTQ3PXrl3P6go1THbOeVJgmFScJx1jrGYAAODdnNY9Az7Gz0+KLfvLTsZ+qTDb1HIAAABcidCMs9cysWzBkI5sN7UUAAAAVyI04+y16lW5fPhH8+oAAABwMUIzzh6hGQAA+AhCM85eXDcpIMS+TGgGAABejNCMs+cfKMX/wb6cvlfKP2VqOQAAAK5CaEbjVO2icWSraWUAAAC4UpOH5gcffFCvvfaa0tPTm/qt4Qr0awYAAD6gyUPzSy+9pKlTp6pNmzZ6/fXXm/rt4WyEZgAA4AOafEbAf//73/r999/16aef6q677pJhGLrjjjuaugw4S0xnKTBcKs4lNAMAAK9lMQzDMOvN58yZo3nz5mnXrl1mleBSWVlZslqtyszMVGRkpNnluM6CUdLBZPvytBQpvLm59QAAADiZ07pnrFy5ssHH3Hfffdq0aZOzSoBZ6KIBAAC8nNNC85gxYzR69Gjt2bOnQcdFREQ4qwSYhdAMAAC8nNNC8/jx47Vq1SpdcMEFeuyxx5Sbm+uspuHuCM0AAMDLOS00f/jhh1q9erU6duyo5557Tl26dNHbb7/trObhzqLbS8FlfbYJzQAAwAs5dci5Sy+9VNu2bdM///lP5ebmavLkyRo4cKC2bNnizLeBu/Hzk1r1tC9nH5ayj5paDgAAgLM5fZxmf39/PfDAA/rtt980efJkfffdd+rbt69uv/12nThxwtlvB3fh0EVjq2llAAAAuILLJjeJjY3VggUL9N133+nCCy/UG2+8oc6dO2vOnDmy2WyueluYxSE08y8LAADAu7h8RsC+fftq+fLluuqqq3Tq1Cn95S9/UWJiopKSklz91mhK3AwIAAC8mNNnBCwsLNSWLVu0YcOGiseBAwcqthuGoZ07d2rEiBG6/vrr9corryg6OtrZZaCpndNWCo2S8jPsodkwJIvF7KoAAACcwmmheerUqdq4caO2b9+ukpISSfaAHBAQoF69eql///4Vj7y8PD388MN677339N1332nNmjVq166ds0qBGSwW+9XmvUlSbpqU9btkbWN2VQAAAE7htGm0/fzsPT2io6PVr1+/ioDct29fhYWF1XjMggULdNttt+mqq67Sxx9/7Iwy3IrPTKNdbs3/Sd/Oti9f97bUbYy59QAAADiJ0640v/HGG+rfv7+6du1a72NuueUWffDBB/r666+dVQbMdHq/ZkIzAADwEk67EfCWW25pUGAul5CQoMzMzEa9d0FBgf76178qKChIU6ZMOas2Zs6cKYvFUuvjbNv1KVVD8++MoAEAALyH028EbKiHH35YF1544Vkfv3btWt1xxx06ceKEiouLG1VLZGSkWrZsWeO22tajishWUniclHucmwEBAIBXMT00d+3a9ayuUEvSe++9p7vuukt///vfFRoaqltuuaVRtVx99dVatGhRo9rwaeU3A+7+Qio4JWXsl6K5wRMAAHg+l4/T7Ert2rXTzz//rLvvvlsWrmi6h9ZV/tWA8ZoBAICX8OjQ3K9fP7Vq1crsMlAVk5wAAAAv5NGh2dn27duniRMnKiEhQUFBQWrRooXGjh2rNWvWmF2a56gamg9tNq8OAAAAJyI0V7Fx40b17t1bGzduVE5OjlauXKmMjAwNHz5czzzzTJ3HFxYWKisry+HhcyLipKiyfsy//yCVFJpbDwAAgBMQmstMnDhRmzZt0oMPPqiWLVsqKChIvXv31ooVKxQVFaUnnnhC69evP2Mbs2bNktVqrXgkJCQ0UfVupm1/+3NpIUPPAQAAr0BoLtO5c2ddcMEF1dZbrVZNmjRJkrRw4cIztjF9+nRlZmZWPFJTU11Sq9s795LK5QNn/osGAACAJzB9yDlP0LFjR0nSrl27zrhfcHCwgoODm6Ik91Z+pVmSDn5nXh0AAABOwpXmerDZbJLEsHb1Fd1eimhhXz64QbKVmlsPAABAIxGaJaWmpio+Pl7bt2+vcfvu3bsl6awnYfE5FktlF42ibOnoDnPrAQAAaCSfCs07duzQgAED9PjjjzusLy0t1bFjx7Rq1apqx2RkZOjdd9+VJE2ZMqUpyvQOdNEAAABexKdC8/z585WcnKxZs2bp5MmTFevLu1089dRTeumll3To0CEVFRVp06ZNGj16tDIzMzVz5kwNHDjQrNI9T9XQfCDZvDoAAACcwGIYhmF2EY1xpn7GCxcudLg6vHr1ak2YMEFDhgzRxx9/7HDsDz/8oKVLl+rLL7/Unj17VFBQoJiYGPXr109Tp07V8OHDG1xbVlaWrFarMjMzFRkZ2eDjPZqtVHq2nVSYKYXHSg/vtnfbAAAA8EAeH5rdmU+HZkl651pp95f25Xs2SzGdzK0HAADgLPlU9ww0MbpoAAAAL0Fohuucy82AAADAOxCa4TqtekkBIfZlZgYEAAAejNAM1wkIklr3ti+fOihl/m5uPQAAAGeJ0AzXYrxmAADgBQjNcK22l1Qu00UDAAB4KEIzXKtNX8nib18+wJVmAADgmQjNcK3gCKllon05bZeUl25uPQAAAGeB0AzXo18zAADwcIRmuN65Vfs1M8kJAADwPIRmuF7V0MyVZgAA4IEIzXC98OZSbFf78uGtUmGOqeUAAAA0FKEZTaP8arNRKh3aZG4tAAAADURoRtNoO6ByOeVr08oAAAA4G4RmNI0OwyRZ7Mt71phaCgAAQEMRmtE0wmOkVj3ty8d2SFlHTC0HAACgIQjNaDodh1cu700yrw4AAIAGIjSj6VQNzXtWm1cHAABAAxGa0XRa95aCrfblvUmSrdTcegAAAOqJ0Iym4x8gdRhqXy44Jf2+xcxqAAAA6o3QjKbVcUTlMl00AACAhyA0o2l1vKxymdAMAAA8BKEZTSuylRTX3b78+w9S7klz6wEAAKgHQjOaXsXVZkNKWWtqKQAAAPVBaEbTY+g5AADgYQjNaHrn9pMCw+3Le9ZINpu59QAAANSB0IymFxAstRtsX849bp9WGwAAwI0RmmEORtEAAAAehNAMczj0a15jXh0AAAD1QGiGOaLbSc072pdTN0gFmebWAwAAcAaEZpin/GqzrUTat87cWgAAAM6A0AzzVO2isfsr8+oAAACoA6EZ5mk7QPIPti/vWSMZhrn1AAAA1ILQDPMEhUnnDbAvZx2SjjL0HAAAcE+EZpir65WVyz//17QyAAAAzoTQDHN1GytZyr6GO/9LFw0AAOCWvCI0FxQU6K9//auCgoI0ZcqUs27HZrNpzpw5SkxMVFhYmGJiYnT11Vdr69atTqsVp4mIs/dtlqT0vdKxnebWAwAAUAOPD81r165Vjx49NG/ePBUXF591OzabTRMmTNBDDz2kO++8U2lpafr++++VnZ2tiy++WF9++aUTq4aD86+qXKaLBgAAcEMeHZrfe+89XX311XrggQf0wgsvNKqt+fPna9myZXrggQd01113KTw8XB07dtSHH36o8PBw3XjjjcrOznZS5XDQbawki32ZLhoAAMANeXRobteunX7++WfdfffdslgsjWpr9uzZkqTbb7/dYX1UVJQmTJig48ePa9GiRY16D9SiWQupbX/78snd0vFd5tYDAABwGo8Ozf369VOrVq0a3c5PP/2klJQUxcbGqlOnTtW2Dxw4UJK0fPnyRr8XakEXDQAA4MY8OjQ7y7Zt2yTZr1zXpH379g77wQVO76IBAADgRgLMLsAdHD16VJIUHR1d4/aoqChJUlpammw2m/z8+LuG00W2lM7tJx38Tjrxq3T8Fymuq9lVwSSGYajEZqik1FCxzabSUvvrUpuhEput7Lnsdakhm2FUrCtfttkMlZYvG4ZKbZLNsK+3GVKpYcgwjLJ1VV+X7WfY6zBqeW1IFcuSZLM5rjNUub+h8mdVvC77X8XnPX171a79Rll7qnq8Krv/V91WdU3F9or9Kht1aN/h3FddX8tONR5nnGnXam3XpfY2Gn/Pg0vvmuCWDAcGJwSNMH1UNyVEh5ldRgVCs6S8vDxJUlBQUI3bg4ODHfaNiIiocb/CwkIVFhZWvM7KynJilT7g/HH20CzZu2jEPWZmNT6tuNSmvMJS5RaVKL+4VPlFpcovLlVekX25oLjKo8SmguJSFZbYVFhsU2GJfbmopHK5uNT+uqjUUHGJTUWl9nUlpYaKSm0qKbWpuNSwrysLwwAA33bPsOpdZs1EaJYUFmb/W0xRUVGN26sG4fJ9azJr1iw99dRTzi3Ol3QbI33+qH1553+loYTms2GzGcrML9ap/GKdyitSZn6xMvOLlVX2nF1QoqyCEuUUlii7oFg5BSXKLihRblGJcgtLlFtUqqISm9kfAwAAt0JolhQfHy9JSk9Pr3F7RkaGJCkuLu6MXTOmT5+uBx98sOJ1VlaWEhISnFipl7O2lhIullI3SGm7pLRfpdguZlflFopKbDqWVaDj2YU6kVOotOyyR06hTmQXKiOvSOm5RcrIswdld75QGxTgpyB/PwX6WxTgb18O8Lco0N9PAX5lz/4WBfr5yd/PogB/iwL87PsG+Fns6/ws8vfzk7+fKp4D/PzkZ7HI30/y87PI32Lf16/suXzZzyL5+1lksVjkb5EsFov8/MrWW+z7qGzZYpH8anqWZCk71r5sP77GdbIfZ5Gk016Xj/pT+dq+U9Xt5eMCWcqOVcV+clg+fVut61X92NO3OLTvsH/1UYpqG7fIsY0zj27UyMGPzqqNxo64dMa2XdYyGsOF/8nhIs3Dg+veqQkRmiX17NlTkrRv374at6ekpEiSevToccZ2goODHbpy4CycP84emiXp50+kIY+YWk5TMAxDJ3OLlJqep0MZ+UrNyNORUwU6klmgo1n5OppZoBM5Nf8riLOEBPopIjhAYUEBCg8OUHiQv8KDAxQW5K/QIH+FBfkrLChAIYH+Cg30V2ign0IC/csefgoO9FdIgL+CA/0UHFD+8K8IyEEB9kdAWVgFAMDTEJolde/eXe3bt1dKSor27Nmjjh07Omxfv369JGns2LFmlOdbzh8rfTHdvrzzv14TmsuDcUparvadyFFKWq5STuRq/4lcHcrIV35xaaPfIzzIX1HhQYoOD1JUWJCiwgJlDbU/IkMdl5uFBKhZsP05IiRAgf7c3AoAwJn4VGjesWOH7rzzTg0ZMkR///vfHbZNmzZNd911l15//XU999xzFetPnTqlpUuXKi4uTlOmTGniin2QtY3Upo90aJN0fKd0YrcU4143AtQlM79Yvx7N1q9Hs/TL0Wz9ejRbvx3LVlZBSYPb8vezqEWzYMVbQxRvDVFcsxDFNgtWbESwYpoFKTYiRDHN7CE5JNDfBZ8GAABIPhaa58+fr+TkZCUnJ+uhhx5S8+bNK7bdfvvt+uqrr/Tiiy+qQ4cO+tOf/qSjR4/qrrvuUnZ2tlasWKFmzZqZWL0POX+cPTRL9lE0Bk8zs5ozyioo1k+HMrXtUKa2Hzql7Ycy9fup/HofHxTgpzZRoUqICrM/R4cpISpMrc4JUatzQhUTESx/P7ozAABgNovhjEEvTXSm/pELFy50uDq8evVqTZgwQUOGDNHHH39c7djS0lLNnTtXb775pnbv3q2wsDANGjRIM2bMUK9evRpcW1ZWlqxWqzIzMxUZGdng433WqVTppQvsy7Fdpbu/d4s7OAzD0KGMfH2fclIb96Vry8EM7U3Lrdexrawh6tiimdrHhKtdTLjax9qfW1lD5UcoBgDA7Xl8aHZnhOZGWHB55ZjNf14tJfQxpYzU9Dx9u/uENuyzB+UjmQVn3D8syF/nt4xU15bN1CU+Ul3jm6lzi2ayhgY2UcUAAMAVfKp7BjzIhTdVhuYti5ssNBeWlGrTvgyt/fW4vv71+BmvJAf6W9StZaR6tLGqR5tzlNjmHHWMi6A7BQAAXogrzS7EleZGKMqVZneRirKlwHDp4V+lYNf0Kc8uKNbqXce0asdRrd9zQnlFNY9kERror4vaRqlvu2j1bRetngnncPMdAAA+givNcE9B4dIf/p/0w0KpOFfa+R/71WcnySks0Zpdx7Ry+xF9/VtajTPg+VmkC8+N0tAuserfMUYXtLIqKICh2QAA8EWEZrivC2+0h2ZJ2rKk0aHZZjO0fu8Jvb8pVat/PqbCGoJyTESQhnSO09AusRrUKUbnhAU16j0BAIB3IDTDfbW6UGpxgXTsJ+nQRun4LimuW4ObOZ5VoA9/OKT3Nx1Uanr14eBimwXrigviNTqxlS46N4rRLAAAQDWEZrgvi0XqdaP0+aP211uWSJf//czHlDEMQ8l7T2pR8n4l/XJcpTbHrvvNw4N0xR9a6soeLdXnvGhu3gMAAGfEjYAuxI2ATpCXLr3QVSotlMKaSw/ukgKCa9291Gbo85+O6rVv9mrH75kO2ywWaVCnWN3QJ0GXdWtB/2QAAFBvXGmGewuLlrqNkX76SMo7Kf26Sup+dbXdCopLtWzLIf17XYr2n8xz2BbXLFjX9UnQhN4JSogOa6rKAQCAFyE0w/1deKM9NEv2LhpVQnNxqU3vbzyoOUl7lJZd6HDY+S0jdefQDrrigngF+HNVGQAAnD1CM9zfeYOlc9pKpw5Ie5OkUwdlWBP0xc6jevbzX7XvhOMEJP07NNedQzpoUKeYM06zDgAAUF+EZrg/Pz/7DYFrn5Fk6Pev39S9h0dqy8FTDrtd3j1edw/roB5tzjGjSgAA4MUIzfAMPSfK+Prvshg2GT++ra2FvSXZu1z0bRetx6/opp4J55haIgAA8F6EZrg9m83Qu7tKdK4tUYMtP6qN5YQG+e3Q7zEDNH1UV13aNY5uGAAAwKUIzXBre9NyNH3ZDm3cn64/+g3V4KAfJUn/aPU/xd71GDf4AQCAJkFohlsqLrXp9XUpennNbhWVTXe92nah0gNbKrr4iFqeWC+l/SzFX2BypQAAwBdwmQ5u5+DJPF3zr2Q9/8WvFYH53OgwLbm1v6KH/6Vyx+RXTKoQAAD4GkIz3MrnPx3Vla98WzGbn59FumNwe33xwGD17xgj9fqTFBpl3/mnj6TMQyZWCwAAfAWhGW6hqMSmpz/9WXe+/YOyC0okSe1iwvXJ1IGafkU3hQb523cMCpf63GpftpVI388zqWIAAOBLCM0w3aGMPF07/zstWL+vYt3oHi21/J4B+kMba/UD+t4u+Qfbl39YLBVkNlGlAADAVxGaYapvd6fpyjn/07bUU5KkIH8//d9V3fXKDb3ULCSw5oMi4qSeN9iXi7KlzQubplgAAOCzCM0wzbIfDunmhZuUmV8sSUqIDtWyu/rrxkvOq3vc5UvulVS2z4bXpJIi1xYLAAB8GqEZTc4wDL26do8e+nCbSmyGJGl4tzituHdQzd0xahLTUep6pX05+4i040MXVQsAAEBoRhMrtRn6239/0vNf/Fqx7sZ+bTX/xt6yhtbSHaM2A+6vXE6eI9lsTqoSAADAEaEZTSa/qFR3vv2D3tlwsGLdo5d31dNXdZe/31lMg53QV0roZ19O+0Xa85WTKgUAAHBEaEaTyMwv1qQ3vtdXPx+TJAX4WfTidYm6a2iHuvsvn8mA+yqX189pZJUAAAA1IzTD5XIKSzRl4UZtOXhKkhQRHKCFN/fR1b3aNL7xzqOk5p3sywf+Jx1IbnybAAAApyE0w6Xyi0p16+JN+rEsMDcPD9IHd/TToE6xznkDPz9p4AOVr7+aIRmGc9oGAAAoQ2iGyxSWlOqOt3/Q9ynpkiRraKCW/PlidW9VzxEy6qvH9VJsV/vyoY3SLyud2z4AAPB5hGa4RHGpTfe++6PW/ZYmyd4l461b+ur8VpHOfzP/AOmyJytfr3laKi1x/vsAAACfRWiG05XaDD20dJu+LLvpLyTQTwum9FFiwjmue9MuV0gJF9uXT/wqbXvXde8FAAB8DqEZTmUYhv723x1avu2wJCkowE9v3NRHfdtFu/aNLRZp+FOVr9fOkorzXfueAADAZxCa4VRv/m+f3tuYKsk+rNy8SRdqYKeYpnnztpfYR9OQpOzD0ob5TfO+AADA6xGa4TTf/Jamv6/aVfH6n9f11GXdWjRtEZc9Kals3Of//VPKz2ja9wcAAF6J0Ayn2HciV/e+u0W2stHe7ru0o8Ymtmr6QlqcL/WcaF8uyJT+92LT1wAAALwOoRmNll1QrNve2qysAvuIFSPOb6EHhnc2r6Ch0yX/YPvyhvlS5u/m1QIAALwCoRmNYrMZeuD9rdpzPEeS1LlFhF68rqf8/BoxNXZjnZMg9b3NvlxSIK39u3m1AAAAr0BoRqO88NWvWvPLcUn2yUv+fVNvRQQHmFyVpEEPScFlY0JvfZvptQEAQKN4fGguKCjQ008/rS5duigkJEQtW7bUjTfeqH379jWonZkzZ8pisdT6mDJlims+gAdbuf2IXl27V5LkZ5FenXih2jYPN7mqMmHR0rC/Vr7+9AGppNC0cgAAgGdzg0uCZy8/P1/Dhw/X9u3btWDBAl155ZXavXu3Jk+erF69eumbb75RYmJivduLjIxUy5Yta9xW23pf9fupfD328faK13+98vymG1quvvreJm3/QDq8xT7hyf9ekoY+anZVAADAA3l0aJ45c6aSk5M1d+5cXXvttZKkxMRELVu2TF26dNGkSZO0bds2+fv716u9q6++WosWLXJhxd7BZjP08NJtyi678W9sYivdMuA8c4uqiZ+/NHaONH+IZJRK386Wul8txZp4kyIAAPBIHts9Iy8vT6+++qqCgoI0efJkh20dOnTQ8OHDtXPnTq1atcqkCr3XwuT9+i7lpCSppTVE/zfuAlksJt74dybxf5D632tfLi2SVjwg2WymlgQAADyPx4bmpKQk5ebmKjExUREREdW2Dxw4UJK0fPnypi7Nq/12LFvPfv5LxesXrk2UNTTQxIrqYcij0jlt7csH1ttvDAQAAGgAjw3N27ZtkyS1a9euxu3t27d32K8+9u3bp4kTJyohIUFBQUFq0aKFxo4dqzVr1jS+YC9QVGLTXz7YqqIS+5XaWwa0U/+ObtaPuSZBYdLoKpOcfPk3Kee4efUAAACP47Gh+ejRo5Kk6OjoGrdHRUU57FcfGzduVO/evbVx40bl5ORo5cqVysjI0PDhw/XMM8/UeXxhYaGysrIcHt7k5TW/aedh+2fqFBehRy7vYnJFDdDxMqnHdfblgkzp88fMrQcAAHgUjw3NeXl5kqSgoKAatwcHBzvsV5eJEydq06ZNevDBB9WyZUsFBQWpd+/eWrFihaKiovTEE09o/fr1Z2xj1qxZslqtFY+EhIQGfCL39sOBdM372j68XKC/RS9e11MhgfW7wdJt/PHvUqj9L1P6aZn062fm1gMAADyGx4bmsLAwSVJRUVGN2wsLCx32q0vnzp11wQUXVFtvtVo1adIkSdLChQvP2Mb06dOVmZlZ8UhNTa3Xe7u73MIS/eWDbbIZ9tcPDO+sC1pbzS3qbITH2INzuU+mSlmHzasHAAB4DI8NzfHx8ZKk9PT0GrdnZGRIcs74yh07dpQk7dq164z7BQcHKzIy0uHhDV5a/ZsOptuv2F/UNkp3DG5vckWNkHiD1HmUfTnvpLTsVqm0xNyaAACA2/PY0NyzZ09JqnXmv5SUFElSjx49Gv1etrIhytx2WDUX2nM8WwvX75ckBQf46YVrExXg77FfG8likcb9S4psbX99YL207jlzawIAAG7PY9PPsGHDFB4eru3btys3N7fa9vL+x2PHjq2zrdTUVMXHx2v79u01bt+9e7ckqWvXro2o2PMYhqGnPv1ZJWX9Mu4c0kHnxbjJNNmNERYtjX9TspT1yf7mOSnlG3NrAgAAbs1jQ3NYWJjuueceFRYWavHixQ7bUlJStHr1anXv3l1XXHFFxfodO3ZowIABevzxxx32Ly0t1bFjx2qcCCUjI0PvvvuuJGnKlCnO/yBu7Iudx/Tt7hOSpNbnhOquoR1MrsiJ2l4iXfrXsheG9PFtDEMHAABq5bGhWZJmzJihSy65RI899piWLVum/Px8bd++XePHj1dYWJjeeecdhym058+fr+TkZM2aNUsnT56sWF/e7eKpp57SSy+9pEOHDqmoqEibNm3S6NGjlZmZqZkzZ1ZMmOILCopL9X8rfq54/cTobp43WkZdBvxFaj/MvpxzTPr4dmYLBAAANfLo0BwaGqqkpCQ99NBDmj59uqKiojRixAh1795dP/74oxITEx32HzdunKKiojRu3DiH8Z3btm2rzZs367777tPixYvVrVs3hYeHa+zYsYqLi9NXX32lGTNmNPXHM9W8r/fq91P5kqSBHWP0x+7xJlfkAn5+0jWvSxEt7K9T1kr/+6e5NQEAALdkMQzDMLsIb5WVlSWr1arMzEyPGkkjNT1Pw//5jQpLbArws+jzBwapY1wzs8tynX3rpMVjJRmSxU+auFTqNMLsqgAAgBvx6CvNcI1nVv6swrKpsqf0P8+7A7MktRssDXnUvmzYpA+nSEfqP/06AADwfoRmOPh2d5q+2HlMkhQTEaz7h3cyuaImMuRRqVvZSCtFOdI7E6RT3jE5DQAAaDxCMyqUlNo0c/nOitePjeqqZiGBJlbUhMr7N7fpa3+dc1R651op/5SpZQEAAPdAaEaFj3/8XXvT7GNe9zr3HF3Tq7XJFTWxwFDphvel6LIZD9N2SUtvlEpqnqodAAD4DkIzJEnFpTbNTdpT8frxK7rJz8/3ZkBUeHNp0kdSaNnoKvvWSZ/eJ3G/LAAAPo3QDEnSf7b8roPpeZLsQ8z1OS+6jiO8WPMO0sQPpIAQ++tt70lJz5hbEwAAMBWhGSoutemVtbsrXj/gKzf/nUlCX3sfZ5Vdbf92trTm/7jiDACAjyI0Qx9vOaTUdPtEJoM6xai3L19lrur8q6RRz1a+/na29NWTBGcAAHwQodnHFZfa9EqVvsxcZT7NxXdIV8yufJ08R/p8OsEZAAAfQ2j2cct+OKRDGZVXmS9qy1XmavreJo1+SRVdNTbMk1Y+KNlsZlYFAACaEKHZhxWV2DR3bdWrzJ1NrMbN9b5ZuupVVQTnzQukT++VbKWmlgUAAJoGodmHLdtSeZV5cOdYXdQ2yuSK3FyvSdI1/5YsZX9sfnxbWnqTVJhjbl0AAMDlCM0+qqjEcVxm+jLXU49rpf+3QPILsL/+ZYW04I/SqYPm1gUAAFyK0OyjPvrhkH4/Zb/KPKRzrC48l6vM9db9avs4zsFW++tjP0n/vlQ6uMHcugAAgMsQmn1Qqc3Qa9/srXjNVeaz0HG4dOvqyim3c9OkxaOlH98xty4AAOAShGYftPaX4w6z//XiKvPZie0s3bpGajfY/rq0SPrkbumLv0qlxebWBgAAnIrQ7IMWf7e/YvnmAeeZVodXCIuW/vSx1OfWynXfzbX3cz65t/bjAACARyE0+5g9x7P17e4TkqRzo8M0tEucyRV5Af9A6coX7JOgWPzt637/QXptkPTDYiZCAQDACxCafczi5AMVyzdd0lb+fhYTq/EyfW+T/vxVZT/n4lzp0/ukD/4k5Z40tzYAANAohGYfklVQrGVbDkmSQgP9dW3vBJMr8kJtLpLu+Fa6cHLlul9WSPMukX77wry6AABAoxCafchHmw8pr8g+g901F7aWNTTQ5Iq8VHCENHaOdN07UmjZtOQ5x6R3J0jvT5IyDpz5eAAA4HYIzT7CZjP0VpUbAKf0P8+0WnxGt9HS3d/Zh6cr98sK6dW+0jfPS8UF5tUGAAAahNDsI775LU37T9qHmRvQsbk6tWhmckU+olm8NOkj+/TbES3s60oKpLXP2Lts7P7K3PoAAEC9EJp9xKLk/RXLky85z7Q6fJLFIvWYIN2zSep3d+UIG+kp0jv/T3rrKil1k7k1AgCAMyI0+4CUtBx981uaJKlNVKgu69bC5Ip8VIhVunyWdOe30rn9K9enfC29OVx651rp8I+mlQcAAGpHaPYBb33HMHNupUV36eZV9i4b57StXL/7S+n1ofabBY/+ZFp5AACgOkKzl8spLNFHP1QOM3dd73NNrgiSKrts3PuDNOZlKbJN5bZfVkivDZDeGmcfps5mM61MAABgR2j2cv/Zckg5hSWSpHG9WssaxjBzbsU/ULpoinTfFvuMgs1aVm5LWWsfpm5ub2nDfKkw27QyAQDwdRbDYI5fV8nKypLValVmZqYiIyNNqeGquf/TtkOZkqTP7h+kbi3NqQP1VJwvbXlL+n6elLHPcVtwpJR4vZR4g9Sql/1qNQAAaBKEZhcyOzTvTcvRZS98I0nq3ipSK+8b1OQ14CzZSu19nL+fJ+37pvr22K72AN3jOimyVdPXBwCAj6F7hhf7z5bfK5av7tXaxErQYH7+UpdR0uTl0l3J0oU3SQEhldvTfpFWz5T+eb697/MPi6Wc42ZVCwCA1+NKswuZeaXZZjM06Lm1+v1Uvvws0vePX6a4ZiF1Hwj3VZAp7fyvtO196WByDTtYpHP7SV2vlLqOlqLbNXWFAAB4LUKzC5kZmjeknNR1r38vSRrSOVaLb+nbpO8PF0tPkbZ9IG17Tzp1oOZ9YrtJHS6VOgyT2vaXgsKbtkYAALwIodmFzAzNjy3brvc3pUqSXr6+p67qSfcMr2QY0uEt0q4V0i8rpRO/1ryfX6CUcLHUYah03mCpZaIUyL88AABQX4RmFzIrNBcUl6rP/7da2QUlCg/y1+a/jVBokH+TvT9MdGK3fZznX1ZKhzZLquWPt3+QfQSOhL72MJ1wsRQR16SlAgDgSQjNLmRWaF65/YimvrtFkjT+wjZ6YUJik7033EheurT/W2nvWvtU3acPYXe6yNb2K9DxPezPLXvY1zG0HQAACjC7gMYqKCjQc889p3feeUcHDhxQVFSUhg8frqefflrt2jXsRiibzaa5c+fqzTff1O7duxUWFqZBgwZpxowZ6tmzp2s+gAv858dDFcvXXEi3DJ8VFi2df5X9IUkZ+6WUb6TUDdLB76X0vY77Z/1uf/y6qnJdaLR9eLu4rvbn2C7254gWhGkAgE/x6CvN+fn5Gj58uLZv364FCxboyiuv1O7duzV58mTt379f33zzjRIT63eV1WazacKECfrkk080Z84c3XTTTTpy5IjuvPNOffvtt/r00081cuTIBtVnxpXmkzmFuvjva1RiMxQfGaL1j10qfz/CDWqQe0JK3WgP0Yc2SUd3SIVZ9Ts2qJl9dI7o9o6Pc861z2ro7/F/HwcAwIFHh+ZHH31Uzz33nObOnaupU6dWrN+7d6+6dOmirl27atu2bfL3r7s/77x583T33Xfr4Ycf1vPPP1+xPiMjQx06dFBgYKD27NmjZs2a1bs+M0LzovX7NPPTnyVJdwxpr+mjujXJ+8IL2GzSqf3SkW3Ske325+M/S9lHGtaOxd8+4Yo1QbK2sT8iW9mvTjdrKTWLty8HBLnkYwAA4AoeG5rz8vIUFxen4uJinTx5UhEREQ7bL7/8cn3xxRdavny5xowZU2d7HTp0UEpKin777Td16tTJYdudd96p+fPna86cObr33nvrXaMZobnqtNlfPDBYXeLrH/KBGuWfkk78Zp9QJe1X6fgu+5B3pw5KRunZtxsabb/5MDxWCo8pe46VQqPsXUtCo8ueo+zLQeF0CQEAmMZj/w01KSlJubm56tOnT7XALEkDBw6sd2j+6aeflJKSotjY2GqBubyt+fPna/ny5Q0KzU1tb1pORWA+v2UkgRnOEXpO2Sgbp431XVpsD87p++whOj1FykwtexyS8k6eud38dPsj7Zf61eEXIIVYHR/BzaTgSPtzUETZ6wj7clC4/RFY9hwUJgWGSYGh9mf/wLM6HQAA3+SxoXnbtm2SVOvNfu3bt3fYr6naMlPVabO5ARAu5x8oNe9gf9SkKFfK/F3KOiRlH7V38yh/zjoi5abZ+1UX59bv/Wwl9iBeVxivL78Ae3gOCLGPWR0QIgUESwGhZc9lr/2Dqj/7B0r+wWXP5a8D7eNhVzwH2J/9AsqWyx9l6/z8KtdZ/O1Tp/v5Vy5b/Cv3s/idtt6/bB1X3gGgqXhsaD569KgkKTo6usbtUVFRDvs5o620tDTZbDb5+fk1uF5Xs9kM/edHe2j2s0hjE1uZXBF8XlC4FNvZ/jiTolx7eM5Nsz/yyq5A52c4LhdkOj4MW+Pqs5XYb3ys782PbslSFp7LHhVhujxQn7bdUva6Yr2lyn5Vt5227PDsJ1lUw7aa1pWtr2v59PBftb0zvq5jXbX1qsc+NexfnzbqtX8jNbjdJv5LFX+Jg7ON+D8pqq3ZVVTw2NCcl5cnSQoKqvlmouDgYIf9nNFW+b41dQeRpMLCQhUWFla8zspquv8z3rQ/Xb+fypckDewUq7hIZnuDhyjvRtGQH0bDkAqzKx9FOWUBOMf+ujjPvq4oVyoqWy7Ok4rza3jOl0oKpZIC+6O0yHWf1ekMe7/y8r7ljehiDgBuZ/A0sytw4LGhOSwsTJJUVFTz/8GVh9fy/ZzRVl3tzZo1S0899VSd7+cKhzPzFRkSoKyCEl3Ti64Z8HIWixQSaX84m81WFp4LpZKiyiBdUmhfV1psf11aVLlcUiTZiu2vbcVSaYn92VZStnzaa6O0bN1pz0apfdkotddhK7FfUa9Yb6vcbtgcH7ay/QzDcb0M+zEyyrYZjtvKX1dbluN6h2fbaevkuB0AvJDHhub4+HhJUnp6eo3bMzIyJEktW7Z0WltxcXFn7Joxffp0PfjggxWvs7KylJCQUOf7O8PVvdpo1AUttfaX4xrSJbZJ3hPwSn5+9psGVfdfuHEG5WG6apCuablq6K56XJ2v61hXbX0ddZ5+bJ1t1LK+QQNSNWDfBg901cR/efHMgbjg7sLdK894bGgun6Fv376apwZOSUmRJPXo0aPJ2goODnboytHUQgL9NeoPdf8lAQBczlJTX2MA8Fzud0dbPQ0bNkzh4eHavn27cnOr332/fv16SdLYsWPrbKt79+5q37690tLStGfPnka1BQAAAO/jsaE5LCxM99xzjwoLC7V48WKHbSkpKVq9erW6d++uK664omL9jh07NGDAAD3++OPV2ps2zd7Z/PXXX3dYf+rUKS1dulRxcXGaMmWK8z8IAAAA3J7HhmZJmjFjhi655BI99thjWrZsmfLz87V9+3aNHz9eYWFheueddxym0J4/f76Sk5M1a9YsnTzpONbr7bffrmuuuUYvvvii5s+fr9zcXO3du1cTJkxQdna2lixZ0qAptAEAAOA9PDo0h4aGKikpSQ899JCmT5+uqKgojRgxQt27d9ePP/6oxMREh/3HjRunqKgojRs3rtqYzH5+flq6dKlmz56tV199VTExMerbt6/CwsK0YcMGjRw5sik/GgAAANyIxTC45dVVsrKyZLValZmZqchIFwyNBQAAgCbh0VeaAQAAgKZAaAYAAADqQGgGAAAA6kBoBgAAAOpAaAYAAADqQGgGAAAA6kBoBgAAAOpAaAYAAADqQGgGAAAA6kBoBgAAAOpAaAYAAADqQGgGAAAA6hBgdgHezDAMSVJWVpbJlQAAAOBMmjVrJovFUut2QrMLZWdnS5ISEhJMrgQAAABnkpmZqcjIyFq3W4zyy6FwOpvNpsOHD9f5NxdnycrKUkJCglJTU8/4H90XcW5qxnmpHeemZpyXmnFease5qRnnpXZmnRuuNJvIz89Pbdq0afL3jYyM5A9gLTg3NeO81I5zUzPOS804L7Xj3NSM81I7dzs33AgIAAAA1IHQDAAAANSB0OxFgoODNWPGDAUHB5tditvh3NSM81I7zk3NOC8147zUjnNTM85L7dz13HAjIAAAAFAHrjQDAAAAdSA0AwAAAHUgNAMAAAB1IDS7sYKCAj399NPq0qWLQkJC1LJlS914443at29fg9uy2WyaM2eOEhMTFRYWppiYGF199dXaunWr8wt3MWedl5kzZ8pisdT6mDJlims+gIsVFBTor3/9q4KCghr1GbzpOyM557x423cmKSlJf/7zn9W5c2eFhIQoLCxM559/vqZNm6a0tLQGt+fM3ywzOfO8TJky5YzfmZkzZ7rmQ7jI5s2b9eSTT2rgwIFq3bq1QkJC1LZtW1155ZVasWJFg9vzlt8ZZ54Xb/udOd3BgwcVGRkpi8Wir7/+usHHm/o7Y8At5eXlGf379zciIiKMpUuXGrm5ucbWrVuNxMREw2q1Glu3bq13W6Wlpcb48eONgIAA41//+peRk5Nj7N6927jsssuMoKAg44svvnDhJ3EuZ56XGTNmGJGRkUaXLl1qfDz22GMu/CSukZSUZHTq1MmIiooyJBmTJ08+q3a86TtjGM47L970nZk7d64hyejRo4exZs0aIysryzh06JDx7LPPGoGBgUaLFi2MX3/9td7tOfPPppmcfV4mT55sNG/evNbvzCuvvOLCT+N8LVq0MIKDg43XXnvNOHbsmJGTk2OsW7fO+MMf/mBIMp544ol6t+VNvzPOPC/e9DtTkxEjRhiSDEnG2rVrG3Ss2b8zhGY39cgjjxiSjLlz5zqs37Nnj+Hv7290797dKCkpqVdb//rXvwxJxsMPP+ywPj093YiKijLi4uKMrKwsp9XuSs48LzNmzDjr8OSO3n33XcNqtRqvvvqqsWDBgkaFQ2/6zjjzvHjTd+b55583goKCjIMHD1bb9uijjxqSjMsvv7ze7Tnzz6aZnH1eJk+ebMyYMcOJFZqrRYsWxlNPPVVt/ebNmw1JRnBwsJGXl1evtrzpd8aZ58WbfmdON3/+fCMqKsro2rXrWYVms39n6J7hhvLy8vTqq68qKChIkydPdtjWoUMHDR8+XDt37tSqVavq1d7s2bMlSbfffrvD+qioKE2YMEHHjx/XokWLnFK7Kzn7vHibdu3a6eeff9bdd98ti8XSqLa85TsjOfe8eJP4+HjdcMMNSkhIqLZt7NixkqTVq1ertLS0zra86c+mM8+LN/rwww81derUauu7dOkiSSosLFRBQUG92vKm3xlnnhdvdfDgQU2bNk0vvviiWrRo0eDj3eF3htDshpKSkpSbm6vExERFRERU2z5w4EBJ0vLly+ts66efflJKSopiY2PVqVOnRrVlNmeeF2/Ur18/tWrVqtHteNN3RnLeefE2f/rTn2oNJFarVZJksVhk1GMof2/6s+nM8+KNBg0apObNm1dbn5ycLEnq3bu3oqKi6mzH235nnHVevNmtt96qAQMGVAu89eUOvzOEZje0bds2SfYrZDVp3769w35N1ZbZXPFZ9u3bp4kTJyohIUFBQUFq0aKFxo4dqzVr1jS+YA/lTd8ZV/CF78wvv/wiyf5/QgEBAXXu7yvfmYael3Jbt27V2LFj1apVKwUHB6t169a6/vrrtXnzZleV2iRKS0t1+PBhLVmyRFOmTNFFF12k999/v17HevN3pjHnpZy3/c78+9//1oYNGzR//vyzbsMdvjOEZjd09OhRSVJ0dHSN28v/tlq+nzPaSktLk81ma3CtTcmZ56Xcxo0b1bt3b23cuFE5OTlauXKlMjIyNHz4cD3zzDONL9oDedN3xhV84Tvz1ltvyWKx6Mknn6zX/q74s+mOGnpeyn377be66qqrtGPHDmVmZurtt9/WTz/9pEsuuUSLFy92UbWuV/4XgJtvvlnjxo3TqlWr1KFDh3od682/M405L+W86XcmNTVVDz/8sGbPnl1jt6f6coffGUKzG8rLy5MkBQUF1bi9fC728v2c0VZ92zOTM8+LJE2cOFGbNm3Sgw8+qJYtWyooKEi9e/fWihUrFBUVpSeeeELr1693TvEexJu+M87mC9+Zzz//XMuXL9eDDz6ooUOH1usYZ//ZdEdnc14k6S9/+Ys2btyoP//5z2revLlCQkI0bNgwffLJJ5Ls/Xk9bUi+ciUlJUpLS9Nnn32mDRs2qGPHjlqyZEm9jvXm35nGnBfJ+35nbrvtNvXp00e33XZbo9pxh98ZQrMbCgsLkyQVFRXVuL2wsNBhP2e0Vd/2zOTM8yJJnTt31gUXXFBtvdVq1aRJkyRJCxcuPJtSPZo3fWeczdu/M7t379bkyZM1fvx4Pfvss/U+ztl/Nt3N2Z4XSUpMTKzxKmOHDh10+eWXq6ioSO+8846zSm1yMTExGjFihNauXauIiAhNnjxZa9eurfM4b/+dOdvzInnX78wbb7yh//3vf3rjjTca3ZY7/M4Qmt1QfHy8JCk9Pb3G7RkZGZKkli1bOq2tuLg4+fm599fBmeelLh07dpQk7dq1q9FteRpv+s40JU//zhw4cEAjRozQwIED9e6778rf37/exzbln82m1pjzUhdP/85UFRkZqUmTJskwDM2bN6/O/X3ld6ah56UunvSdOXTokB566CH94x//0Hnnndfo9tzhd8Zzv4lerGfPnpJU6z/ZpaSkSJJ69OjRpG2ZrSk/S3kfOl8cosybvjNNyZO/M3v27NHgwYM1aNAgLV26tNZ//qyNt35nGnte6uLJ35malN+gtX///jr39dbvTE0acl7q4knfmdWrVysrK0v33ntvtVkNv/nmG0nSsGHDKtbVNTugO3xnCM1uaNiwYQoPD9f27duVm5tbbXt5X6byMUPPpHv37mrfvr3S0tK0Z8+eRrVlNmeel9TUVMXHx2v79u01bt+9e7ckqWvXro2o2DN503fGmbz1O/Pzzz9r8ODBGjlypBYvXuxwJXX27NlKTU2tsw1n/tl0F844L8nJyYqPj6+4AnY6T/vOLFmyRDExMbUOt3f48GFJqnHotdN50++MM8+LN/3OTJkyRYZ9Er1qjyFDhkiS1q5dW7GurnsF3OJ3xmXTpqBRymeeevXVVx3W79271wgICKg268327duN/v37G9OnT6/W1rx58wxJxrRp0xzWZ2RkGNHR0R4165Kzzsu+ffsMScasWbOqvUd6erphtVoNSca3337rmg/iYgsXLqxz5jtf+c5U1Zjz4o3fma1btxqxsbHGPffcY9hstmrbddqMXUlJSUafPn2MOXPmVNu3oX823ZmzzsvatWsNScZ7771XrY1ff/3V8Pf3NwIDA42UlBSnfwZXKP/zk5SUVG1bXl6e0aFDB0OS8a9//ativS/8zjjzvHjj70xNhgwZUuuMgO78O0NodlN5eXnGJZdcYjRr1sz46KOPjLy8PGPbtm1Gz549jcjIyGrzq0+dOrViLvcTJ044bCstLTWuueYaIyAgwHjttdeMnJwcY8+ePcaIESOMwMBA44svvmjKj9Yozjov+/fvNyQZISEhxosvvmikpqYahYWFxsaNG43+/fsbkoyZM2c29cdzmvqEQ1/5zlTVmPPibd+ZTZs2GVFRUUZYWJhx3XXX1fg4/f/UrrzySkOSERERUa29hv7ZdFfOPC9ff/21IcmIjo42Fi5caBw9etTIz883kpKSjG7duhn+/v7GG2+80cSf8OwtWrTIkGTEx8cbixcvNg4fPmxkZ2cbycnJFSFo1KhRRnFxccUxvvA748zz4m2/M7U5U2h2598ZQrMby8/PN2bOnGl06tTJCA4ONuLi4oxJkyYZe/furbbvV199ZURFRRnjxo2r8cpISUmJ8dJLLxl/+MMfjJCQECM6Otq46qqrjC1btjTFR3EqZ52XzZs3G4888ojRs2dPIyIiwggICDDi4+ONcePGGV999VVTfRynKv8RrumxcOFCh3196TvjrPPiTd+Z+++//4znpfxR9f/U3nrrLaNZs2bGfffdV2ObDfmz6a6ceV5sNpuxbt06Y+rUqcb5559vhIaGGoGBgUZCQoIxceJEY9OmTU386RqnuLjY+Oyzz4ybb77Z6NSpkxEaGmoEBAQYsbGxxogRI4xFixYZpaWlDsf4wu+Ms8+LN/3OVFV+Fb2mx5AhQyr2c+ffGYth+OhcoAAAAEA9cSMgAAAAUAdCMwAAAFAHQjMAAABQB0IzAAAAUAdCMwAAAFAHQjMAAABQB0IzAAAAUAdCMwAAAFAHQjMAAABQB0IzAAAAUAdCMwD4mC+++EIWi0V33HHHWR0/ZMgQtW/f3slVAYB7IzQDgI/ZsGGDJKlv374NPjYtLU3r16/X1Vdf7eyyAMCtEZoBwMeUh+Y+ffo0+NhPPvlEpaWlhGYAPofQDAA+4rnnnpPFYtGqVaskSYmJibJYLLJYLLrsssvq1cZ//vMftWjRQv3793dlqQDgdgLMLgAA0DTi4uJ07bXX6sMPP1RcXJxGjRpVsW3YsGF1Hp+dna01a9Zo8uTJ8vOr3zWX5cuX66WXXtLOnTuVk5Oj8847T5dcconGjx/v8P4A4O4IzQDgI6ZMmaL27dvrww8/1B//+EctWrSoQcevXLlShYWFuuaaa+q1/+uvv6477rhDzzzzjN5//31FRERo06ZNuv322/Xmm2/KMIyz+BQAYA66ZwCAD9m6daskqWfPng0+9j//+Y+sVqsuvfTSeu0/f/58RUdH669//avi4uIUFhamIUOGaO7cuQ1+bwAwG6EZAHzI2YbmwsJCrVq1SldeeaUCAwPrdYzFYlFGRoY+++wzh/VDhw7V7t27G/T+AGA2QjMA+JBt27ZJst8E2BBfffWVcnJy6t01Q5IefPBBWSwWXXHFFerXr59eeOEF/fbbbwoMDFTHjh0b9P4AYDaLQacyAPAJJSUlioiIUGxsrFJTUxt07C233KL3339faWlpCg8Pr/dxGzZs0Isvvqjly5crPz9fktSvXz/NmjVLQ4cObVANAGAmrjQDgI/YtWuXCgsLG3yVubS0VJ9++qlGjhzZoMAsSRdffLHef/99nThxQsuWLdPo0aP1/fff67LLLlNycnKD2gIAMxGaAcBH7Nq1S5J0wQUXNOi4b7/9VidOnGjwhCZffvmljh8/LkkKCwvTNddco08//VRPPPGEbDabli5d2qD2AMBMhGYA8BHFxcWSpLy8vAYd9/HHHysgIEBjxoxp0HG33367Pvzww2rry29CDAsLa1B7AGAm+jQDgI/Ys2ePzj//fNlsNo0cOVJxcXEKCQnRa6+9dsbjzj33XHXu3FmrV69u0Pudd955ys3N1WuvvaZhw4YpNDRUW7du1R133KFDhw7phx9+ULt27RrzkQCgyXClGQB8RMeOHfXBBx8oMTFR69at0+LFi7V58+YzHrNp0yalpqY2uGuGJL333nu6/vrr9fTTT6tTp05q3ry5Jk+erMGDB2vLli0EZgAehSvNAIBaPf744/rHP/6hQ4cOqVWrVmaXAwCmITQDAGrVrVs3Wa1Wff/992aXAgCmIjQDAAAAdaBPMwAAAFAHQjMAAABQB0IzAAAAUAdCMwAAAFAHQjMAAABQB0IzAAAAUAdCMwAAAFAHQjMAAABQB0IzAAAAUAdCMwAAAFAHQjMAAABQh/8fmbh/kNGMJoYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time,y_analytical(time,0),label=r'$y$')\n",
    "plt.plot(time,dy_dt(y_analytical(time,0)),label=r'$\\dot{y}$')\n",
    "plt.legend()\n",
    "plt.xlabel(r'$t$ / s')\n",
    "plt.ylabel(r'$y$, $\\dot{y}$ / s$^{-1}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to get the same result with a generalized midpoint rule. It can be considered as a weighted average between an explicit forward Euler and an implicit backward Euler scheme. For $\\alpha > 0$ the update proceeds implicitly with the *recursive* update formula\n",
    "\n",
    "$y_{n+1}^{i+1} = y_n + [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1}^i)] (t_{n+1} - t_n)$\n",
    "\n",
    "We define the corresponding Euler update as\n",
    "\n",
    "$\\Delta y_\\text{MP}^{i} = y_{n+1}^{i+1} - y_n = [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1}^i)] \\Delta t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Mid_Point_Update(y_i,y_n,al,dt):\n",
    "    return ((1.-al)*dy_dt(y_n) + al*dy_dt(y_i))*dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the first step is a forward Euler step. After that, a mixture of BE and FE schemes occurs.\n",
    "\n",
    "We now need to define the time stepping itself. So we choose a time-step size and the same starting value (initial condition) of $y_0 = y(t=0) = 0$ to initialize our solution field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, for every time step we append to this array the value of the next time step by applying the generalized Midpoint update rule. Each new time step value then serves as the initial condition for the following time step.\n",
    "\n",
    "We will *recursively* re-fresh the update with a new estimate for the rate until the estimated solution does not change anymore, i.e.\n",
    "\n",
    "$\\left| y_{n+1}^{i+1} - y_{n+1}^{i} \\right| < \\epsilon_\\text{abs} \\cup \\left| \\frac{y_{n+1}^{i+1} - y_{n+1}^{i}}{y_{n+1}^{i}} \\right| < \\epsilon_\\text{rel}$\n",
    "\n",
    "with suitably chosen absolute and/or relative tolerances.\n",
    "\n",
    "For simplicity, let's set only an absolute tolerance of $10^{-6}$ here.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate(dt,alpha):\n",
    "    #Startwerte\n",
    "    t_end = 4.\n",
    "    absolute_tolerance = 1.e-6\n",
    "    max_iter = 1000\n",
    "    iteration_counter = np.array([0])\n",
    "    y = np.array([0.])\n",
    "    times = np.array([0.])\n",
    "    #\n",
    "    while times[-1]+dt < t_end: #repeat the loop as long as the final time step is below the end point\n",
    "        times = np.append(times,times[-1]+dt) #here define the next time point as the previous time point plus the time increment dt\n",
    "        y_old = y[-1] #Starting value for recursive update\n",
    "        i = 0\n",
    "        #\n",
    "        while True:\n",
    "            y_new = y[-1]+Mid_Point_Update(y_old,y[-1],alpha,dt) #one backward Euler step\n",
    "            i += 1\n",
    "            #if change is below tolerance, above maximum iterations, or scheme is fully explicit, stop iterations\n",
    "            if (np.abs(y_new - y_old) < absolute_tolerance or i > max_iter or alpha == 0.0): \n",
    "                break\n",
    "            y_old = y_new #preparation of next recursion\n",
    "\n",
    "        y = np.append(y,y_new) #append the new found solution to the solution vector\n",
    "        iteration_counter = np.append(iteration_counter,i) #store how much iterations this time step took to converge\n",
    "    return times, y, iteration_counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c4070b1de7cb4a999b65b2ef82f9f38a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=0.1, description='$\\\\Delta t$ / s', max=1.0, min=0.001, step=0.1)…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "@interact(dt=widgets.BoundedFloatText(value=0.1,min=1e-3,max=1,step=0.1,description=r'$\\Delta t$ / s'),\n",
    "          alpha=widgets.RadioButtons(options=[0,0.5,1], value=1, description=r'$\\alpha$'))\n",
    "\n",
    "def plot_CN(dt=0.1,alpha=1):\n",
    "    time_n, y_n, iter = integrate(dt,alpha)\n",
    "    fig, ax = plt.subplots(ncols=2,nrows=2,figsize=(18,12))\n",
    "    #solution\n",
    "    ax[0][0].plot(time,y_analytical(time,0),label=r'analytical')\n",
    "    ax[0][0].plot(time_n,y_n,label=r'numerical',marker='d',ls=':')\n",
    "    ax[0][0].legend()\n",
    "    ax[0][0].set_xlabel(r'$t$ / s')\n",
    "    ax[0][0].set_ylabel(r'$y$')\n",
    "    #error\n",
    "    ax[0][1].plot(time_n,y_n-y_analytical(time_n,0),ls=':',marker='d',label=r'absolute')\n",
    "    ax[0][1].plot(time_n[1:],(y_n[1:]-y_analytical(time_n[1:],0))/y_analytical(time_n[1:],0),\n",
    "               ls=':',marker='d',label=r'relative')\n",
    "    ax[0][1].legend()\n",
    "    ax[0][1].set_xlabel(r'$t$ / s')\n",
    "    ax[0][1].set_ylabel(r'error')\n",
    "    #iterations\n",
    "    ax[1][0].plot(iter)\n",
    "    ax[1][0].set_xlabel(r'time step')\n",
    "    ax[1][0].set_ylabel(r'number of iterations')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the error:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "Tasks:\n",
    "\n",
    "- How does the peak relative error change when you divide the time step size by 10?\n",
    "- How does the solution change when you increase the rate constant by a factor of 2 and why?\n",
    "- How does the solution change again for $k=4$ and time step sices of $\\Delta t= 0.2$ and $\\Delta t = 0.3$?\n",
    "- Can you recover the forward and Backward Euler solutions?\n"
   ]
  }
 ],
 "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
}