04_Generalized_Midpoint_Newton.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 4 - Generalized midpoint rule, Newton 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": 37,
      "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": 38,
      "metadata": {},
      "outputs": [],
      "source": [
        "y_eq = 1.\n",
        "k = 2."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "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": 40,
      "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": 41,
      "metadata": {},
      "outputs": [],
      "source": [
        "time = np.linspace(0,4,100)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 42,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0, 0.5, '$y$, $\\dot{y}$ / s$^{-1}$')"
            ]
          },
          "execution_count": 42,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtQAAAIWCAYAAAB3DRpPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACBO0lEQVR4nO3dd3hUVf4G8PdOySSTMukZEpJQQ+8GJFRFEBRBV1wFBQG3CYr82HVdhBUUlMW2C4qCqAQVWUQUyyKI2ECkkyAl9HTSSDKTOplyf39MMmRIAimT3Cnv53nmmZlbzv3OdQxvbs49RxBFUQQRERERETWLTOoCiIiIiIhcGQM1EREREVELMFATEREREbUAAzURERERUQswUBMRERERtQADNRERERFRCzBQExERERG1AAM1EREREVELMFBLQBRF6PV6cE4dIiIiItfHQC2BkpISaDQalJSUSF0KEREREbUQAzURERERUQswUBMRERERtQADNRERERFRCzBQExERERG1AAM1EREREVELMFATEREREbUAAzURERERUQsopC6AiIiIiBrPaDTCbDZLXYZLksvlUCqVDm+XgZqIiIjIBej1ehQUFMBgMEhdiktTqVQIDQ1FQECAw9pkoCYiIiJycnq9HllZWfDz80NoaCiUSiUEQZC6LJciiiKMRiN0Oh2ysrIAwGGhmoGaiIiIyMkVFBTAz88P7du3Z5BuAR8fH/j7+yMzMxMFBQUOC9S8KZGIiIjIiRmNRhgMBmg0GoZpBxAEARqNBgaDAUaj0SFtMlATERERObGaGxBb42Y6T1VzLh11c6dLBWqDwYAtW7Zg4sSJ0Gq1CAkJQVhYGO6++2589913TW7vww8/RHx8PMLDw6HVavHAAw/g3LlzDW6fkpKCKVOmQKvVIjw8HPHx8di0aVNLPhIRERFRo/DqtOM4+ly6VKB++umn8dBDD6Fv3744f/48rl69iqNHj6Kqqgpjx47Fm2++2ei2Fi9ejJkzZ2Lu3LnIzc1FSkoKTCYT4uPjceLEiTrbHz9+HIMHD4bZbEZKSgpyc3Mxd+5czJgxA0uXLnXgpyQiIiIiVyKIoihKXURjPfHEEzhx4gR+/vlnu+UFBQWIiYmBxWJBTk4OAgMDb9jO0aNHER8fj+nTp2Pjxo225TqdDtHR0YiLi8Phw4dtv72IoogBAwbg0qVLyMjIgEajse0zY8YMbNq0CUePHkX//v0b9Tn0ej00Gg10Op1Dh2whIiIi91NZWYnLly+jY8eO8Pb2lroct+Doc+pSV6jHjx+PRYsW1VkeGhqK7t27w2Aw4NixYzdtZ82aNRBFEVOnTrVbrtFoMGHCBBw9ehT79++3Ld+7dy+Sk5Nx11132YVpAJg6dSosFgvWrFnTzE9FRERERK7MpQL1xIkTceedd9a7rqqqCgAQEhJy03b27NkDABg0aFCddTXLdu/e3eztiYiIiMhzuMU41AUFBTh//jx69uyJvn373nDb8vJypKenw8vLC2FhYXXWR0VFAbDegFjjzJkzdutqCwsLg1KpRFpaGioqKuDj49OSj9IqvjjzLS4WpsFgrsLCkXOlLoeIiIjIrbhFoF69ejVMJhNWr15907s2i4uLAQBqtbre9b6+vgCAoqKiOvvUrKtNEASo1WrodDoUFxfXG6gNBoPdNKF6vf6GNTra/vQjuFycAUEQUGWqgpfCq02PT0REROTOXKrLR30OHDiAf/3rX3jhhRcwZswYqcup14oVK6DRaGyP6OjoNj1+tCYSgPXmykx9Tpsem4iIiMgRXnzxRQiC0OCQxRs3boQgCFixYkUbV+bigfr06dOYOHEi5s2bV+/NivWpGQGkvLy83vVlZWUAgKCgoDr71KyrTRRFW1sNjS6ycOFC6HQ62yMjI6NRtTpKTGCk7XWGLrtNj01ERETkCEOHDgUAHDp0qM66kpIS/OMf/0Dnzp2xYMGCti7Ndbt8nDx5EnfccQdmz56Nl19+udH7qdVqxMTEID09Hfn5+XX6UWdlZQEAunfvblvWo0cPu3W15efnw2g0IjY2tsH+0yqVCiqVqtE1OlqM5lrf73Rd3c9AREREruv//v0jikoMN99QQkH+Kvz7/0a3qI3BgwdDLpfj4MGDddY9//zzyMnJwRdffCFJ5nLJQH38+HGMGzcOc+fOtZtUJTU1FV5eXoiMjGx4ZwBjxozBhg0bcPToUYwfP95u3dGjRwEAY8eOtdv+hRdesK272fbOpqbLBwCk8wo1ERGRWykqMeCqrlLqMlqdn58f+vTpg6SkJFRVVcHLy3pP2NmzZ7F69WqMGzcOkyZNkqQ2lwvUhw8fxp133ol//OMf+Pvf/263bunSpejQoYMtZIuiiMzMzDp9lufOnYvExERs3rzZLlDrdDp88803GDhwIBISEmzLR4wYgX79+mHHjh3Q6XR2Y1Fv3rwZMpkMc+c67+gZwT6B8FX6oMxYwS4fREREbibIX7q/gjeWo2pMSEhAUlISkpOTER8fDwCYP38+RFHEf/7zH4ccozlcKlDv378fEyZMQGRkJMrLy+tM+Z2UlIQOHTrY3s+bNw9vvvkm5s2bh1WrVtmWDxo0CM8++yxWrFiB22+/HTNmzIBer8fMmTMBAO+//77daCGCIGDDhg0YOXIkZs2ahQ0bNiAgIAAbN27Epk2bsHjx4kbPkigFQRAQExiFM/kXUFhRjNKqMvh51R2xhIiIiFxPS7tSuJKEhAS89dZbOHjwIOLj4/Hll19i586dmD9/vq2LrhRcKlC//PLL0Ov10Ov1eP755+vd5t5777W9jo6OhlqtrndUjeXLlyMuLg6rVq3C008/DUEQMGLECBw+fBjdunWrs/2AAQNw+PBhLFq0CHFxcbBYLIiNjcXGjRvxyCOPOOwztpZoTSTO5F8AAKQXZ6NneFeJKyIiIiJqmto3JhoMBixYsABhYWFYsmSJpHUJoiiKklbggfR6PTQaDXQ6HQICAtrkmN9e+AnvHv0vAOCxgQ/hzq6j2uS4RERE1DKVlZW4fPkyOnbsCG9vb6nLkZxWq4VGo8Gjjz6KRYsW4Z133sEf//jHJrXh6HPq0sPmUeNxpA8iIiJyB0OHDsX58+fx4osvYsCAAXjsscekLomB2lNwpA8iIiJyBwkJCbZ5QFavXg2ZTPo4K30F1CZ8vdQI8bFOVpOhywZ7+hAREZEr6tixIwBg6tSpGD58uMTVWDFQe5CaGRPLjRW4WlEkcTVERERETSOKIv79739DrVY3aWK/1sZA7UGia/ejLma3DyIiInItb7/9Nvbv348lS5agffv2Updj41LD5lHLxNj1o87CwMjeElZDREREdHOZmZn46KOPcO7cOXz44YcYPnw4/vrXv0pdlh0Gag9SO1BzxkQiIiJyBd9++y0WLlyI4OBgPPTQQ1i1ahXkcrnUZdlhoPYgUQFayAQZLKKFI30QERGRS5g9ezZmz54tdRk3xD7UHkQpV6KdfzgAIEufA5PFLHFFRERERK6PgdrD1IxHbbKYkFOaJ3E1RERERK6PgdrDxHCkDyIiIiKHYqD2MNeP9EFERERELcNA7WE40gcRERGRYzFQe5hwv1Co5F4AwJE+iIiIiByAgdrDyAQZ2mvaAQDySgtQaTJIXBERERGRa2Og9kA1NyaKEJGpuyJxNURERESujYHaA0Xb3ZjIbh9ERERELcFA7YE40gcRERGR4zBQe6CYwGtjUWcwUBMRERG1CAO1B9Ko/OGv8gMApLMPNREREbmo9PR0pKSkSF0GA7UnEgTB1u1DV6mHvrJE4oqIiIiImm7GjBno0aOH1GUwUHsquynI2e2DiIiIqNkYqD0UR/ogIiIiV/fjjz9CFEWpy2Cg9lQxDNREREREDsFA7aFqX6FOK86UsBIiIiIi18ZA7aF8lN7Q+oUBsF6hNlvMEldERERE5JoYqD1Yh6BoAIDRbMSVkjyJqyEiIiJq2IsvvghBELBp06Z612/cuBGCIGDFihVtXBkDtUfrENje9vpyUYaElRARERHd2NChQwEAhw4dqrOupKQE//jHP9C5c2csWLCgrUtjoPZkHauvUANAajEDNRERETmvwYMHQy6X4+DBg3XWPf/888jJycHrr78OlUrV5rUxUHuwDoEM1EREROQa/Pz80KdPHyQlJaGqqsq2/OzZs1i9ejXGjRuHSZMmSVKbQpKjklMI8tFA4x0AXaUeqUWZEEURgiBIXRYRERE1wT++XYHiSr3UZdxQoHcA/jVuYYvbSUhIQFJSEpKTkxEfHw8AmD9/PkRRxH/+858Wt99cDNQerkNgeyTnnEZJVRmuVhQhVB0sdUlERETUBMWVehRWFEtdRptISEjAW2+9hYMHDyI+Ph5ffvkldu7cifnz50s6BTkDtYfrGBSN5JzTAIDUokwGaiIiIhcT6B0gdQk35agaa9+YaDAYsGDBAoSFhWHJkiUOab+5GKg9XO2RPlKLM3FLVF8JqyEiIqKmckRXClfRqVMnRERE4ODBg3jttddw8eJFvPPOOwgMDJS0Lt6U6OHsAjWHziMiIiInN3ToUJw/fx4vvvgiBgwYgMcee0zqkhioPZ3WLxwqhXV4GY70QURERM4uISEBoiiivLwcq1evhkwmfZyVvgKSlEwmQ6wmCgCQV3YVZVXlEldERERE1LCOHTsCAKZOnYrhw4dLXI2VywbqkydPIiEhAYIgIDU1tUn7jh49GhqNBlqtts4jIiICgiBg7ty5dfYJDg6ud5/Jkyc78JO1vQ5B17p9pBVnSlgJERERUcNEUcS///1vqNVqvPzyy1KXY+NyNyVWVlZi2bJlWLt2bYsu8a9atQozZ86ss3zHjh24++67MXXq1DrrPvvsM4wePbrZx3RWtSd4uVyUgZ7hcRJWQ0RERFS/t99+G/v378fKlSvRvn37m+/QRlwuUM+fPx85OTlITk7GI488gp9++qnJbQwePBhRUVH1rlu7di369u3rNH9CaAv2U5DzCjURERE5j8zMTHz00Uc4d+4cPvzwQwwfPhx//etfpS7LjssF6oULFyI2NrZFbTT0J4KMjAzs2LEDb775ZovadzXRmkjIBBksooUjfRAREZFT+fbbb7Fw4UIEBwfjoYcewqpVqyCXy6Uuy47LBeqWhukbWb9+PXx9ffHII4+02jGckZdciagALTJ02cjUX4HRbIRSrpS6LCIiIiLMnj0bs2fPlrqMG3LZmxIdzWQy4b333sP06dPh5+dX7zZbtmzB8OHDER0dDa1Wi1GjRmHdunUwm81tXK3j1YxHbRYtyNBdkbgaIiIiItfBQF3tq6++QnZ2NubMmdPgNqmpqVi3bh0yMjKQnJyM2267DXPmzMHEiRNhNBob3M9gMECv19s9nA37URMRERE1DwN1tbVr12LUqFHo2bNnveu3bt2Kr7/+Gr169QIAREREYOnSpZg2bRp27tyJt99+u8G2V6xYAY1GY3tER0c3uK1UOGMiERERUfMwUAO4dOkSdu/efcOr02FhYfV2gH/ggQcAANu2bWtw34ULF0Kn09keGRnOF1hrD53HGROJiIiIGs/lbkpsDevWrYNWq8V9993X5H0jIyMBAFeuNNzvWKVSQaVSNbu+tuCn8kWoOhgF5YVILc6ERbRAJvD3LSIiIqKb8fjEVFVVhQ0bNuAPf/gDlMr6R7ZISkrC+vXr612XnZ0NAAgPD2+1GttKh+p+1JUmA/JKCySuhoiIiGoTRVHqEtyGo8+lWwdqURRv2r1i27ZtKCwsxJ///OcGt0lKSsIzzzyDysrKOus+//xzAMCkSZNaVqwTqN2P+jK7fRARETmFmi6nNxoAgZqm5lw6ajxrtw7U8+bNQ0xMDJ566qkGt1m7di0mTZrU4MyJNYqKijB9+nRkZWUBAEpLS7Fy5Up88MEHGDZsGObNm+fQ2qVgN9JHEUf6ICIicgZKpRIqlQo6nY5XqR1AFEXodDqoVKoGeyc0lcv1oT506JDtanBhYSEAID4+HnK5HA8//DBee+0127bR0dFQq9UNjqqRkpKCn3/+Gbt3777hMe+//34IgoBt27Zh5MiRKCkpQWVlJbp164aXX34ZTz75JLy8vBz0CaVjN9IHr1ATERE5jdDQUGRlZSEzMxMajQZKpRKCIEhdlksRRRFGoxE6nQ6lpaU3vZjaFILIX3XanF6vh0ajgU6nQ0BAgNTl2IiiiNnb/4ayqnIEeWuwbvK/pC6JiIiIqun1ehQUFMBgMEhdiktTqVQIDQ11aAZzuSvU1HoEQUCHwPY4lXcORZU6FFfqEejtPIGfiIjIkwUEBCAgIABGo9EtZmmWglwud1g3j9oYqMlOx8BonMo7B8Daj7p/u/onuiEiIiJpKJXKVgmF1HxufVMiNV2HWjcmXipKk7ASIiIiItfAQE12OgXH2F5fKkyXsBIiIiIi18BATXYi/SPgrbDO6niRV6iJiIiIboqBmuzIBBk6BVmvUl8tL0JxpV7iioiIiIicGwM11dE5ONb2+lIhr1ITERER3QgDNdVRO1BfYKAmIiIiuiEGaqqDV6iJiIiIGo+BmuoI9w2Fn5cvAOBiYRo4mSYRERFRwxioqQ5BENC5evg8naEEVyuKJK6IiIiIyHkxUFO9OgVd6/Zxkd0+iIiIiBrEQE31qt2PmoGaiIiIqGEM1FQvBmoiIiKixmGgpnoF+wQi0DsAgHWkD96YSERERFQ/Bmqql/XGROtV6jJjBXJL8yWuiIiIiMg5MVBTgzjBCxEREdHNMVBTgzjBCxEREdHNMVBTgzrXHjqviIGaiIiIqD4M1NSgAG9/hKmDAQCXijJgsVgkroiIiIjI+TBQ0w11qu72YTAZkFWSI3E1RERERM6HgZpuiONRExEREd0YAzXdEAM1ERER0Y0xUNMNdQqKsb3mSB9EREREdTFQ0w35eqnRzj8cAJBanAmT2SRxRURERETOhYGabqpm+DyjxYQM/RWJqyEiIiJyLgzUdFP2/ahTpSuEiIiIyAkxUNNN2QfqdAkrISIiInI+DNR0Ux2CoiEIAgBeoSYiIiK6HgM13ZS3QoX2Ae0AAOm6bFSZqiSuiIiIiMh5MFBTo3QJ7gAAsIgWXCrKkLYYIiIiIifCQE2N0jWko+31uauXJKyEiIiIyLkwUFOjxDFQExEREdWLgZoapb2mHXyU3gCA8wWXIYqixBUREREROQcGamoUmSBD12DrVeqiSh0KygslroiIiIjIOTBQU6PFhbLbBxEREdH1GKip0eJCOtlenyu4LGElRERERM6DgZoarUtIB9trXqEmIiIisnLZQH3y5EkkJCRAEASkpqY2ad/U1FTI5XJotdp6Hx9//HG9+6WkpGDKlCnQarUIDw9HfHw8Nm3a5IBP4xr8vHwR5a8FAKQWZXCCFyIiIiIACqkLaKrKykosW7YMa9euhUzW/N8HoqOjmxTEjx8/jlGjRmHMmDFISUmBRqPBxo0bMWPGDJw/fx5Lly5tdi2upGtoR2SV5MAsWnCpKB3dw7pIXRIRERGRpFzuCvX8+fNx6tQpJCcno1evXm1yTFEUMWvWLABAYmIiAgMDIQgCZs6ciYcffhjLli1DUlJSm9QiNbt+1Oz2QUREROR6gXrhwoXYvn072rdv32bH3Lt3L5KTk3HXXXdBo9HYrZs6dSosFgvWrFnTZvVIyW6CF96YSEREROR6gTo2NrbNj7lnzx4AwKBBg+qsq1m2e/fuNq1JKrUneDl39RIneCEiIiKP53KB2lHKy8uxYMEC9OrVCxEREejUqROmTZuGw4cP19n2zJkzAICoqKg668LCwqBUKpGWloaKiopWr1tqtSd4Ka7UI58TvBAREZGH89hAXVRUBK1Wi3379iE7Oxvbt29HRkYGhg4disTERLtti4uLAQC+vr512hEEAWq12m676xkMBuj1eruHK7Ob4KWA/aiJiIjIs3lkoI6OjsaVK1fw97//HUFBQZDL5ejbty+++OIL+Pn5Yc6cOcjNzXXY8VasWAGNRmN7REdHO6xtKfDGRCIiIqJrPDJQy+VyhIaG1lkeHByM22+/HRUVFdixY4dteWBgIACgrKyszj6iKKK8vNxuu+stXLgQOp3O9sjIyGj5h5BQ7Qlezl/ljYlERETk2VxuHOrWFhkZCQC4cuWKbVmPHj0AAFlZWXW2z8/Ph9FoRGxsLHx8fOptU6VSQaVStUK10vDz8kVUgBZZ+hzbBC9eCi+pyyIiIiKShEdeoU5MTKz35kMAyM7OBgCEh4fblo0ZMwYAcPTo0Trb1ywbO3aso8t0ajXdPmomeCEiIiLyVG4dqEVRrLd7RWJiIjZs2FBneXFxMX788Ud4eXlh/PjxtuUjRoxAv379sGPHDuh0Ort9Nm/eDJlMhrlz5zr+Azgxu/Go2Y+aiIiIPJhbB+p58+YhJiYGTz31VJ1169evx7vvvouqqioAwMWLFzFlyhQUFxfjlVdesZs4RhAEbNiwwTZjok6ngyiKSExMxKZNm7B48WL079+/rT6WU4gLrXVjIid4ISIiIg/mcoH60KFD0Gq10Gq12L9/PwAgPj4eWq0Wf/3rX+22jY6OhlqtrjOqxrp167B48WK8++676NixI4KDg5GQkAA/Pz/s2bMH8+bNq3PcAQMG4PDhwxAEAXFxcQgPD8ebb76JjRs34vnnn2+9D+ykogK0UCutfcY5wQsRERF5MkFkEmpzer0eGo0GOp0OAQEBUpfTbC/+tBrJOdZJb96cuBzhviESV0RERETU9lzuCjU5j64hnOCFiIiIiIGamo0TvBARERExUFML1J7ghVeoiYiIyFMxUFOz+Xn5IlpjnQjncnEGKoyVEldERERE1PYYqKlFeoR1AWAd8/ssr1ITERGRB2KgphapCdQAkFJwXsJKiIiIiKTBQE0t0iO0q+31mfwLElZCREREJA0GamqRYHUgIvzCAADnr6aiymyUuCIiIiKitsVATS1W0+3DZDHhwtVUaYshIiIiamMM1NRiPUKv9aM+k89+1ERERORZGKipxXqEX+tHnVLAftRERETkWRioqcUifEMR5KMBAKQUXILZYpa4IiIiIqK2w0BNLSYIAnqEWa9SG0wGXC7KkLgiIiIiorbDQE0O0TOsdj9qdvsgIiIiz8FATQ7RvfaNiexHTURERB6EgZocor2mHfy8fAEAKfkXYBEtEldERERE1DYYqMkhZIIM3au7fZRWlSFTd0XiioiIiIjaBgM1OQz7URMREZEnYqAmh2E/aiIiIvJEDNTkMB2DoqFSqABYZ0wURVHiioiIiIhaHwM1OYxcJkf30E4AgKIKHXLLCiSuiIiIiKj1MVCTQ9VM8AIAZ/LOS1gJERERUdtgoCaHYj9qIiIi8jQM1ORQXUI6QCFTAOBIH0REROQZGKjJobzkSnQN6QAAyC3NR2F5saT1EBEREbU2BmpyuB61x6MuYD9qIiIicm8M1ORwtW9MPJV7TsJKiIiIiFofAzU5XLfQzpDL5ACAk3lnJa6GiIiIqHUxUJPDeStUiAvpCADIKc1HQVmhxBURERERtR4GamoVvcO72V7/lpsiYSVERERErYuBmlpFn4jutte/sdsHERERuTEGamoVXYI7QKVQAQBO5qZAFEWJKyIiIiJqHQzU1CoUcgV6Vg+fV1ypR5Y+R+KKiIiIiFoHAzW1mt7htbp9sB81ERERuSkGamo1fSJq3ZjIftRERETkphioqdXEBEbBX+UHADiddw5mi1niioiIiIgcj4GaWo1MkKFXeBwAoNxYgctFGRJXREREROR4DNTUqvqwHzURERG5OZcN1CdPnkRCQgIEQUBqamqT9s3MzMTy5ctxyy23ICQkBIGBgejSpQueeOIJXLlypd59Ro8ejeDgYGi12jqPyZMnO+ATuafa/ahP5jFQExERkftRSF1AU1VWVmLZsmVYu3YtZLLm/T7Qs2dPeHt746OPPsIdd9wBANi1axceeeQRbN26Ffv370fnzp3r7PfZZ59h9OjRLSnf40T4hSFUHYyC8kKkFFxCldkIL7lS6rKIiIiIHMblrlDPnz8fp06dQnJyMnr16tWsNiwWC1asWIFx48ZBJpNBJpNhwoQJeOGFF5CXl4dFixY5uGrPJQgCeldfpTaajThXcFHiioiIiIgcy+UC9cKFC7F9+3a0b9++2W0888wzuOeee+osHzFiBABg//79zW6b6rLvR83h84iIiMi9uFygjo2NbXEb//znPxEeHl5neVVVFQAgJCSkxcega3rX7kfNGxOJiIjIzbhcoG5NNVemp06dWu/6LVu2YPjw4YiOjoZWq8WoUaOwbt06mM03Hl/ZYDBAr9fbPTxJkI8G7QPaAQAuFKWhvKpC4oqIiIiIHIeBulplZSXWrFmD7t2744knnqh3m9TUVKxbtw4ZGRlITk7Gbbfdhjlz5mDixIkwGo0Ntr1ixQpoNBrbIzo6urU+htOquUotiiJO55+TuBoiIiIix2Ggrvb3v/8d+fn52Lp1K9RqdZ31W7duxddff227ETIiIgJLly7FtGnTsHPnTrz99tsNtr1w4ULodDrbIyPD8yY46RPBftRERETknhioAaxcuRIbNmzAjh070Lt373q3CQsLg1wur7P8gQceAABs27atwfZVKhUCAgLsHp6mZ1hXCIIAgP2oiYiIyL14fKB+5ZVXsGLFCuzatQu33nprk/ePjIwEgAYnhCErXy81OgdZbyjN0F9BcYVO4oqIiIiIHMOjA/Xy5cuxcuVK7NmzBwkJCbblR44csY34AQBJSUlYv359vW1kZ2cDQL2jhpC92t0+knPOSFgJERERkeO4daAWRbHB/sqLFy/GG2+8gR9++AGDBg2yWxcfH28LyoA1UD/zzDOorKys087nn38OAJg0aZIDK3dP/dv1tL0+nnNKwkqIiIiIHMflph5vinnz5uHNN9/EvHnzsGrVKtvyp59+Gq+++ioeeughbNu27Yb9n2sUFRVh+vTp+M9//oOoqCiUlpZizZo1+OCDDzBs2DDMmzevNT+KW+ga0gk+Sm9UGCtxIucMLBZLs6ePJyIiInIWLheoDx06ZLsaXFhYCMB6RVkul+Phhx/Ga6+9Zts2OjoaarXabpi64uJivPrqqwCA//73v4065v333w9BELBt2zaMHDkSJSUlqKysRLdu3fDyyy/jySefhJeXl6M+ottSyOToG9EDBzOPo7SqDBcKUxEX2knqsoiIiIhaRBBFUZS6CE+j1+uh0Wig0+k8bsSP7y/9grWHPwIATOl1N37fe6LEFRERERG1DP/eTm2qn/ZaP+qkK+xHTURERK6PgZraVIg6CDGaKADAxcI06A2lEldERERE1DIM1NTmakb7ECHiRM5piashIiIiahmHBuqioiJ88MEHjmyS3FB/bS/b6+Ps9kFEREQuzqGBOj09HbNmzXJkk+SGuod2hrdCBQBIzjkNi2iRuCIiIiKi5mvSsHnp6ek3XF97MhSihijkCvSO6I4jWcnQG0pxuSgDnYNjpS6LiIiIqFmaFKg7dOgAQRAaXC+K4g3XE9UYoO2FI1nJAKyjfTBQExERkatqUqAOCgrCsmXLMGrUqHrXp6Sk4Pe//71DCiP3Vnsa8qQrp3B/r7skrIaIiIio+ZoUqAcNGoSioiL06tWr3vUmkwmcJ4YaI8w3BFEBWmTpc3Cu8DJKq8rg5+UrdVlERERETdakmxIff/xxxMY2/Kf5mJgYbNiwocVFkWeoGe1DFEWcyEmRuBoiIiKi5uHU4xLw5KnHazuRcwbLf1oNABjdYSjmDJkhcUVERERETceJXUgy3cO6QCX3AgAk5ZxidyEiIiJySS0O1HK5HOfOnXNELeRhvORK9AqPAwAUV+qRVpwpcUVERERETdfiQM2ritQS/dtdu8E1idOQExERkQtilw+SVO1AzWnIiYiIyBUxUJOktH5haOcXDgA4W3ARpYYyiSsiIiIiahoGapLcoKi+AACLaMGxKyclroaIiIioaRioSXLx1YEaAI5knZCwEiIiIqKmY6AmycWFdIJ/9SyJSTmnYDQbJa6IiIiIqPEYqElycpkcAyP7AAAqTQaczDsrcUVEREREjdfiQP3MM88gJCTEEbWQB4uP6md7fZjdPoiIiMiFtDhQr1ixgoGaWqyvtgeUMgUA4Gj2CVhEi8QVERERETUOu3yQU/BWqNAnojsAoKhCh0uF6RJXRERERNQ4DNTkNG6p1e3jSHayhJUQERERNZ7DAvX58+fxwgsv4K233rIt+/rrr9GpUycEBgZi/vz5nKacbuiW6hsTAeBI1m8SVkJERETUeA4L1O+//z5ef/11XL16FQCQm5uLqVOnIjs7G507d8batWuxfv16Rx2O3FCgjwZdgzsAANJ1WcgrLZC2ICIiIqJGcFig3r17Nz799FP885//BAB88MEHKCsrw9atW3H06FFs2bIF7733nqMOR27qFrvRPtjtg4iIiJyfwwJ1bm4u7rjjDtv77du3Iy4uDvfccw8AYPLkycjMzHTU4chNxdv1o+bweUREROT8HBaoa/ePzs3NxcGDBzFlyhS7bQRBcNThyE1FBWih9QsDAJzJv4BSQ5nEFRERERHdmMMCdXh4OPbt2wcAePvttyGKIu6//37b+vT0dPj5+TnqcOSmBEHALZF9AQAW0YJjV05KXBERERHRjTksUM+ePRsTJkzALbfcgmXLlmHIkCEYMGAAAODQoUOYPXs2Bg8e7KjDkRuzGz6PsyYSERGRk3NYoJ47dy6eeuopVFRUYPz48di0aRMAwGKx4MEHH8TFixcxbdo0Rx2O3Fi30E7w9/IFACTlnILRbJS4IiIiIqKGCSIHh25zer0eGo0GOp0OAQEBUpfjlN48mIifUw8CAJ4d+QT6t+slcUVERERE9eNMieSUao/2cSgzSbpCiIiIiG6CgZqcUj9tT3jJlQCAg1lJMFvMEldEREREVD8GanJK3goVBrazTkVeYijF6fzzEldEREREVD8GanJat0YPtL3+NeOYhJUQERERNcxlA/XJkyeRkJAAQRCQmprarDY+/PBDxMfHIzw8HFqtFg888ADOnTvX4PYpKSmYMmUKtFotwsPDER8fbxvNhBxvYGTva90+Mo+z2wcRERE5pWYH6rVr1yIvL8+RtTRKZWUlFi1ahFGjRuH8+eZ3A1i8eDFmzpyJuXPnIjc3FykpKTCZTIiPj8eJE3XHPj5+/DgGDx4Ms9mMlJQU5ObmYu7cuZgxYwaWLl3agk9EDWG3DyIiInIFzQ7Uc+bMQVRUFEaNGoU333wTWVlZjqyrQfPnz8epU6eQnJyMXr2aN5Ta0aNH8dJLL+GRRx7BzJkzIQgCAgMDkZiYCFEUMXv2bLup1EVRxKxZswAAiYmJCAwMhCAImDlzJh5++GEsW7YMSUlJjvh4dB12+yAiIiJn1+xAnZaWhpdffhlmsxlPPfUUYmNjMXToULz++uvN7oLRGAsXLsT27dvRvn37ZrexZs0aiKKIqVOn2i3XaDSYMGECjh49iv3799uW7927F8nJybjrrrug0Wjs9pk6dSosFgvWrFnT7HqoYez2QURERM6u2YE6Ojoa//d//4d9+/YhMzMTq1atgo+PD5555hl07twZgwYNwr/+9a8b9klujtjY2Ba3sWfPHgDAoEGD6qyrWbZ79+5mb0+Ow24fRERE5OwUjmikXbt2mDt3LubOnYuCggJ8/vnn2LZtG5YsWYJFixahZ8+eeOCBB/C73/0OvXv3dsQhm628vBzp6enw8vJCWFhYnfVRUVEArDcg1jhz5ozdutrCwsKgVCqRlpaGiooK+Pj4tFLlnuvW6IE4kGnt7vFrxjH0iegucUVERGSxiBBFERZRhEUEREut16JYvR62bUQR1vW1lou12hFFQMR1++LafqLl2vtr219bbhEBiPVsU10TYH8M8frtay+3tXH9fvbvUf0ZRNH+tYiaDerZT6xZimvtVB8bttcNtXHtOKi1TKy9IRpu8/pl9m3Uv77udmKd49TZr/6X9bYtNvDm+i1r79u5fSAmDO1Qb51ScUigri00NBR//OMf8cc//hHFxcX44osvsG3bNqxYsQLPP/88unbtildffRUTJ0509KEbpbi4GACgVqvrXe/r6wsAKCoqqrNPzbraBEGAWq2GTqdDcXFxvYHaYDDAYDDY3uv1+uaW75Fqun1UmY04mHkcjw18EHKZXOqyiMiDWSwizBYLjCYLTGYRJrPl2sNkgdkiwmiywGwWYbJYYDZbtzObLTBZRLv3Zsu1Z5PZ2q51mfW19Vii3XY1y66ts9RZVvvZItZ6XRN8G3otwu69WGtZTdi11J+5iNpEQt8q9w/UtQUGBuLRRx/Fo48+itLSUnz11VfYtm0bzp49K1mglkLNLxPUPDXdPg5kHrN1++BVaiLPZTZbYDCaUWW0oMpohqH6Yax5b7K+NhjNMJqs2xlNZlSZrOuNtZ6NZguMRmswrtnGZK55bw3HRnOt51qBmYioRqsG6tr8/PwwderUOjcCtrXAwEAA1q4f9SkrKwMABAUF1dmnZl1toija2qrZ7noLFy7EggULbO/1ej2io6ObWrpHY7cPItciiiIMRjMqDCZUGEyoNNR6XWVCpcGECoMZlVXWZYYqMyqrrO8NtZ4NRuvymteGKjNMZovUH88lyQRAJhMgEwTr8w1fAzJBgFC9TC4TIFTvLwgC5MJ172WC9Ri12hGq26i9XKjVriCD3TIAtuMJAiCger1MgFC9ThCqt7Fublsmq95BwHXbXLfMvh0AuHYs23Jrw/avUd/2td5ft85+H2sNqN6u5rV1+bV9q3erWQXB9plqvbe9rm7n2q61tqt1zNrrb7jPtbZrt2HXTu3lDdTVmDZQz+a1F13/Oevj662sf4WE2ixQOwu1Wo2YmBikp6cjPz+/Tj/qmuH/une/Fth69Ohht662/Px8GI1GxMbGNth/WqVSQaVSOeojeCR2+yBqO6IoorLKjLIKI8oqjCit9VxeaURZpRFlFSbr6wojyg0mVFRa35cbTCivNKGi0ujS3QIUcgFKhQwKuRxKhQCFQg6lXIBCLoNCIYNSLoNcbn1WKGRQ1Kyrfsjt3guQy2VQyKqfa72XXb9cZt1XJhOgkMkglwmQya2B1fq4tr7m/bXX15bXBFi5vHbAbSCdEFGLeVygBoAxY8Zgw4YNOHr0KMaPH2+37ujRowCAsWPH2m3/wgsv2NbdbHtyPHb7IGqeKqMZutIq6MsMKCmvQkmZEfryKujLqqzvy6tQWm5EaXkVSsqNKK2wvneWLg1eChlUXnKolPLqZwW8lNZlXkrrQ1X9sL6XWZ8V1tdKhRwqpQzK6m2UChm8FNZnpUIGpVIGpbzWe4U1BMtkDJ9E1HhuHahFUURmZmad7hVz585FYmIiNm/ebBeodTodvvnmGwwcOBAJCQm25SNGjEC/fv2wY8cO6HQ6u7GoN2/eDJlMhrlz57b+B/Jw7PZBZL0xrKS8CkUlBhSXVKK4tArFttcG6EqroCs1QF9mDdEVhrYbu91LKYfaWwG1SgEfbwV8VAqoVUr41Hrv4yWHt0oB79qvveTw9lJAVf3sXb1cpZQz2BKRS3DrQD1v3jy8+eabmDdvHlatWmVbPmjQIDz77LNYsWIFbr/9dsyYMQN6vR4zZ84EALz//vvX9eERsGHDBowcORKzZs3Chg0bEBAQgI0bN2LTpk1YvHgx+vfv38afzvOw2we5M1EUoS+rQqG+Eld11kehrgJX9ZUo0htQWFKJYn0likoMrXb1WO2tgJ/aC34+Svj5KOFb+1mthK+39bWvtxJqbwV8fZRQeyvhWx2W5fJmT21AROTSXC5QHzp0CJMmTQIAFBYWAgDi4+Mhl8vx8MMP47XXXrNtGx0dDbVaXe8NgMuXL0dcXBxWrVqFp59+GoIgYMSIETh8+DC6detWZ/sBAwbg8OHDWLRoEeLi4mCxWBAbG4uNGzfikUceaaVPS7Vd3+3jVN459NX2kLosokapNJiQV1SOvKIK5BeVI7+4AgXFFSgorrQ+6ypgNDnmZjtBAPx8vKDx84LGT4UAXy/bw19d/ezrhQC19bkmNCsYiImImkUQGxrBuwlWrlyJBQsWQKl0vrsunZFer4dGo4FOp0NAQIDU5biUAxnH8Pr+9QCAkR2G4IkhM6UtiKiayWxBXlE5cq6WI/dqGXKuliOnsAx5hdYQrS+ralH7ggBo/FQI9vdGYIAKQf4qBPl7I9BfBY2fCkF+Kttrf7WSV4uJiNqQQwK1TCZD165d8frrr+Puu+92RF1ujYG6+arMRvzpi2dQbqyAt0KFdyavhLeCI6hQ2zCZLci5WobsgjJcKShDdn4psgus7wuKyps9qoXaW4HQQB+EanwQovFGsMYbIQHeCNH4IDjA+l7j68WQTETkpBzS5WPDhg1YuHAhJk2ahDvvvBP/+c9/EBcX54imiex4yZW4NXogvr/0CypNBhzJSsbw2MFSl0VuprzSiPTcEmTmliIzrwSZeaXIzCtFztWyJvdflglAsMYH4UE+CA9SI6zWc2igD8ICfaB2wjFViYio8RxyhRoASktLsXTpUrzxxhsAgCeffBLPPfccr8DWg1eoW+Z03jks/eHfAIAB7Xph4cgnJK6IXJXBaEZ6jh5pV0qQnluCtBw90nNKUFBc0aR2fL0V0Ib6QhviC22w2vocokZEsC/CgnzYN5mIyM01OlCfPXsWPXr0gEwmg8lkuuF28+fPx65duxAeHo6XXnoJs2fPdljB7oCBumUsogVPfP1PFJQXQibIsHbSCgR68zzSjRXpK3ExS4fL2TqkZutx+YoOWXmlje6m4aWQISrcD1Fh1ke7UF9EhvohMswXAb5enDSDiMiDNTpQZ2RkYPr06RAEAT/88MNNt//yyy+xYMECXLp0CYMGDcLq1asxdOjQFhfsDhioW27ziS/w+ZmdAICZAx7AXXG3S1wROZOrugpczNThQmYxLmQW42JmMQr1hkbt6+utQIw2ADFaf0RH+KN9uB/ah/sjLNCHYyITEVG9HNbloz5VVVV47bXX8OKLL6KiogLTpk3DypUrERkZ2VqHdAkM1C2Xqb+CBd+8AADoHBSLFeP+IXFFJJVKgwkXMotxLr0IZ9OLcDatCFd1lTfdTyGXISbCHx0iA9AxMgAx2gDEav0RHODNq81ERNQkDh+H2mQyITk5GQcPHsShQ4dw8OBBVFRUQBRFbNq0Cdu3b8fSpUuxYMEC/qNFzdY+oB06BcXgUlE6LhalIUufg6gArdRlURsoKK7AmcuFOH35Kk5fLkRqjh6Wm/Tb8PVWoHP7QHSK0qBTlAYdIzVoH+7Hvs1EROQQDgnUmzdvtoXnpKQkGAwG1Fz4joiIwOTJk5GQkICOHTvi9ddfx9NPP42vvvoK27dvR2BgoCNKIA80InYwLhWlAwD2ph3EQ30mS1wROZooisguKMNvFwpw6tJVnL58FXlFN75h0EelQNfoQHSNDkTn9oHo0j4Q2hA1f4EnIqJW47BxqAFALpejT58+SEhIQEJCAoYOHYqOHTvW2f6tt97C/PnzMW3aNCQmJrb08C6HXT4co7hSj798uRAW0YIwdTDemLgMMoFXHF2ZKIrILSzHiQsF+O1CAU5cKEChvuHuG4IAxGoD0C02CHExQegWE4T2Ef6Qs68zERG1IYdcoV6+fDmGDh2KwYMHw9fX96bbz5kzB0eOHMFXX33liMOThwr0DkDfiO5IyjmN/PJCnC24iB5hXaUui5qorMKIExfycexsPo6fzUNuYXmD23op5egWE4SeHYPRs2MIusUGwdeHYzgTEZG0HBKon3322SbvExcXh8LCQkccnjzYyA5DkJRzGgDwc+ohBmoXIIoiLmbpcORMLo6l5OFselGDfaBVXnL07BCMPl1C0adLKDpHBUKp4F8hiIjIuTj8psTGmj59OiIiIqQ6PLmJW6L6QaVQwWAy4NeMo5g18PfwkvOKpbOprDIh+Vw+Dp/JxeHTuQ1241DIBfToEIJ+Xa0Bumt0EAM0ERE5PckCdVRUFGbNmiXV4clNeCtUGNK+P35OPYhyYwWOZf+GW6MHSl0WAdCXVeHgySvY/9sVnDifjyqTpd7tosL8MLB7OAbEhaFP51B4qyT7sURERNQs/JeLXN7I2CH4OfUgAGBv2iEGagkV6ivx629X8Otv2fjt4tV6u3J4KWTo2zUMg3tGYFD3CIQHqyWolIiIyHEYqMnl9Q7vhiBvDYoqdTiW/Rt0lXpoOBV5m9GXVeGXE9n46VgmTl++ivrGDQoOUCG+pxaDe2rRt2sovL34o4eIiNwH/1UjlyeTyTCywxB8kfItzKIFP6UexKTuY6Uuy61VGEw4eCoHPx3LxPGzeTDXcyU6IliNhL6RSOjbDnHRQZy2m4iI3BYDNbmF2zsNwxcp3wIAvr/0C+7pdgcn8nAwi0XEyUsF+O5QOvb/dgWGKnOdbaLC/DC8XyQS+kaiY2QA/xsQEZFHYKAmt9DOPxy9wuNwKu8csktycbbgIrqHdZG6LLeQV1iOPYfTsedIRr1jRIdqvDFiQHuMGhCFTlEahmgiIvI4bRaoT5w4gS5dukCt5g1I1Dpu7zgMp/LOAQD2XPqFgboFTGYLDp7Kwc79qUi+kF+nX7SvjxLD+0Vi9MD26NkxhN05iIjIozlk6vHGkMlk8PHxwaxZs/Dqq6/C29u7LQ7rlDj1eOuoMhvx5y//gbKqcnjJlVg36V/w9eIvcE1xVVeBXQfSsOtAWp2xogUB6N81DHcMjsGtvdvBSymXqEoiIiLn0mZXqKdOnYr09HSsW7cOx44dw08//QSlkhNwkON4yZUYETsYO8//iCqzEb+kH8a4LqOkLsvpiaKI3y4W4Ot9l3HwVE6doe60IWrcER+D226JRngQf0EhIiK6XpsF6k2bNgEALly4gClTpmDt2rV48skn2+rw5CHGdBqGned/BGDt9sFA3TCjyYK9SVn44qeLuJSts1snE4D4nlrcNawj+ncNY5cOIiKiG2jzmxK7dOmC3bt3s6sDtYrYwPboHByLi4VpuFyUgUuF6egUHCN1WU6lpLwKO39Nxdf7Ltfp1hHkr8K4W2Nx55AOCAvykahCIiIi1+KQQL1y5UosWLCg0V04wsLCHHFYonqN6TQcFwvTAFiH0GOgtrqqq8BnP17ArgNpdYa86xIdiHtHdsawfpFQyGUSVUhEROSaHHJTokwmQ9euXfH666/j7rvvdkRdbo03JbauCmMl/vTlP2AwGeCj9MY7k1ZCpfCSuizJ5BaWY9v357H7UDpMZottuSAAt/Zuh8kjO6Nnx2AOd0dERNRMDrkUtWHDBpSUlGDSpEm46667cO7cOUc0S9QsPkpvJEQPAmAN179mHJW4Imlk5ZfiP/89hj+v+A7f/JpqC9NeSjnuHtYRa/8xBs/OHIxenUIYpomIiFrAYcPmlZaWYunSpXjjjTcAAE8++SSee+45XoGtB69Qt75zBZeweM8rAIDuoZ3xwpi/SVxR28ktLMfmb1Pww5EM1B6ww0clx10JHTF5VGcE+XvusJVERESO5vBxqM+ePYv58+dj165dCA8Px0svvYTZs2c78hAuj4G69YmiiL/tXIYM/RUAwL8nLEFUgFbiqlpXkb4Sn+w5h52/psJkvva/ta+PEpNGdMI9IzrBX+25XV+IiIhai8PvPurWrRu++eYbbN++HX5+fvjDH/6A+Ph4/Prrr44+FFGDBEHA7Z2G2d7vubhPwmpaV2l5FT7YcRp/XPEdvt532RamfX2UmD6hB95fPBbT7uzOME1ERNRKWu12/kmTJuH06dN48cUXcebMGQwfPhzTp09HdnZ2ax2SyM7IDkOglFkHsvkh9VcYTFUSV+RYJrMFX+69iD++9B227jlvG7lD5SXH7++Iw7uLxuL3d8RB7c0JlIiIiFqTw7t8mEwmJCcn4+DBgzh06BAOHjyIc+fOoeYwvr6+WLp0KRYsWOCxN0Kxy0fbWXNwI35KPQAA+NMtD+OOzsMlrqjlRFHEkTO5eO/LU8jKL7UtV8hlmJDQAQ+M6co+0kRERG3IIeNQb9682Raek5KSYDAYbAE6IiICkydPRkJCAjp27IjXX38dTz/9NL766its374dgYGBjiiBqF4Tuo62Beqd53/EmE7DXPoXubQcPd774iSOn8u3Wz56UHtMH98D4cGcGpyIiKitOWwcagCQy+Xo06cPEhISkJCQgKFDh6Jjx451tn/rrbcwf/58TJs2DYmJiS09vMvhFeq2tfi7V3Du6iUAwNLb/g89w+MkrqjpyiqM+GjnGez45bLdyB09OgTjD5N7Iy4mSLriiIiIPJxDrlAvX74cQ4cOxeDBg+Hr63vT7efMmYMjR47gq6++csThiW5ofNdRtkC98/xPLhWoRVHEvuRsvPvFbyjUG2zLw4N8MPPuXhjeP9Klr7gTERG5A4cE6meffbbJ+8TFxaGwsNARhye6oVvbD8QH3ttQXKnHoawkFJQXIlQdLHVZN3WloAxrPzuBY2fzbMu8lHL8/o6uuHdUF6iUcgmrIyIiohoOCdTNMX36dEREREh1ePIgCrkCd3QegU9P/Q8W0YLdF/Ziat/JUpfVIKPJgs9+OI8t352D0XRtqvD4nhH48319EcF+0kRERE7F4aN80M2xD3XbK6rQYc5Xz8IsWuCv8sPb97wEL7nzDSd3IbMYq/57HKlX9LZloRpv/Om+vri1t5bdO4iIiJxQq41D3ZpSUlIwZcoUaLVahIeHIz4+Hps2bWr0/qNHj4ZGo4FWq63ziIiIgCAImDt3bp19goOD691n8mTnvdpJVkE+GgyJHggAKDGU4tf0oxJXZM9osmDTzhT8bdXPtjAtkwm4d1RnvPXMGAzt045hmoiIyElJ1uWjuY4fP45Ro0ZhzJgxSElJgUajwcaNGzFjxgycP38eS5cubVQ7q1atwsyZM+ss37FjB+6++25MnTq1zrrPPvsMo0ePbtkHIMlM6Doa+9OPALAOoTeywxCnCKmXsnT4z3+P4XL2tavSHdoFYP5DA9C5faB0hREREVGjuFSgFkURs2bNAgAkJiZCo9EAAGbOnInvv/8ey5Ytw7333ov+/fvfsJ3BgwcjKiqq3nVr165F3759MXy4608AQvbiQjqhY1A0Lhdl4GJRGs5fvYy40E6S1WM2W/DJd+ew5btzMFePhSeXCXhgTBx+f0cclAqX/AMSERGRx3Gpf7H37t2L5ORk3HXXXbYwXWPq1KmwWCxYs2bNTdt5+eWXMXbs2DrLMzIysGPHDjz++OMOq5mchyAIGN9ltO39zvM/SlZLXmE5Fr71Cz7+9qwtTHdoF4DXnhqJh8d3Z5gmIiJyIS71r/aePXsAAIMGDaqzrmbZ7t27m93++vXr4evri0ceeaTZbZBzGxZzC/y9rGOl/5p5DMUVujavYV9yFua99gPOpFqHjZTJBDw4Ng6vzx/FLh5EREQuyKUC9ZkzZwCg3u4aYWFhUCqVSEtLQ0VFRZPbNplMeO+99zB9+nT4+fnVu82WLVswfPhwREdHQ6vVYtSoUVi3bh3MZnOTj0fS8FJ4YUxna3ces8WMnRd+arNjVxpMeOOTJKz84AjKKk0AgPBgNVbOHY5HxvfgVWkiIiIX5VL/ghcXFwNAvbMxCoIAtVptt11TfPXVV8jOzsacOXMa3CY1NRXr1q1DRkYGkpOTcdttt2HOnDmYOHEijEZjg/sZDAbo9Xq7B0nnzi6jIBesX/1dF35CpbGy1Y95OVuH//vPT/j2YJpt2cj+UVi9YDS6d3D+SWaIiIioYS4VqFvT2rVrMWrUKPTs2bPe9Vu3bsXXX3+NXr16AQAiIiKwdOlSTJs2DTt37sTbb7/dYNsrVqyARqOxPaKjo1vlM1DjhKiDMDx2MACgrKoc31/e36rH+/5IOv626mdk5pUCALy95HjqwQH42yOD4OvjfGNhExERUdO4VKAODAwEAJSVldVZJ4oiysvL7bZrrEuXLmH37t03vDodFhYGubzuVM8PPPAAAGDbtm0N7rtw4ULodDrbIyMjo0n1kePd0+0O2+uvz+6ByeL4bjtGkwVvb0vGvzcfR1X1jIedojT49/+Nwh2DY5xiyD4iIiJqOZcaNq9Hjx4AgKysrDrr8vPzYTQaERsbCx8fnya1u27dOmi1Wtx3331NrikyMhIAcOXKlQa3UalUUKlUTW6bWk9MYBQGtuuNY1dOoqC8EL+mH8WIDoMd1v5VXQX+tfEwUtKKbMvuvDUWf76vD5SKur+YERERketyqSvUY8aMAQAcPVp3lruaZfUNh3cjVVVV2LBhA/7whz9Aqaz/z+9JSUlYv359veuys7MBAOHh4U06LklvUvdxttdfpnwLURQd0u7JiwWY/++fbGFaqZDhyd/3xxMP9GeYJiIickMuFahHjBiBfv36YceOHdDp7Ic727x5M2Qymd2U4aIo3rR7xbZt21BYWIg///nPDW6TlJSEZ555BpWVdW9e+/zzzwEAkyZNaspHISfQI6wLugZ3AACk6bKQnHOmxW3u2H8Zi9buR3GJAQAQFuSDlU8Mx7ghsS1um4iIiJyTSwVqQRCwYcMG24yJOp0OoigiMTERmzZtwuLFi+1mSZw3bx5iYmLw1FNPNdjm2rVrMWnSpAZnTqxRVFSE6dOn27qblJaWYuXKlfjggw8wbNgwzJs3zyGfkdqOIAiY1MP+KnVzmS0i1m//DW9vOwFL9UQt/buG4d/zR6FrdFCLayUiIiLn5VJ9qAFgwIABOHz4MBYtWoS4uDhYLBbExsZi48aNdSZkiY6OhlqtbnBUjZSUFPz88883nQzm/vvvhyAI2LZtG0aOHImSkhJUVlaiW7duePnll/Hkk0/Cy8vLYZ+R2k58ZD+08wvHldI8nMw7i4uFaegc3LSryRUGE1756AgOn861LbtvdBc8endPyGW88ZCIiMjdCaKjOo5So+n1emg0Guh0OgQEBEhdjsf77uJevHPkYwDA0OhB+L+EPzR634LiCix77yAuZVu7IMllAh6/vx/uvJVdPIiIiDyFS3X5IGoNIzvcCo3KHwBwIPMYckvzG7XfhYxi/HXVz7Yw7eutwPN/HMowTURE5GEYqMnjecmVmBB3GwDrjaxfn91z032OpeThH2/tQ6HeeqOqNkSNV+aNRL+4sFatlYiIiJwPAzURgHFdRkKlsI4V/v3l/Siu0DW47d7jWVj2/gEYqqyTwfToEIxX541EdIR/m9RKREREzoWBmgiAn5cvxnYaDgAwmo34MqX+G1V37L+MVzYdgclsvfUgoW87LP9LAjR+nLiHiIjIUzFQE1Wb1H0slHLr5D7fXvwZxZV62zpRFLFl91m8ve0Eam7jHTckFn+fHg8vJSdrISIi8mQM1ETVAn00GNd5JACgqtZVaotFxLtfnsRHO1Ns295/Wxc88UA/DotHREREDNREtU2ufZX6wk+4Wq7Dqi3H8eXPl2zbzJrYEzMn9oIgMEwTERERAzWRnUAfDcZ2HgHAepX6hS824fsj1unrZQIw7/f98bvbukpZIhERETkZBmqi60zuPg5KmfUqdbZ4ClAYIJcJeGZGPMYO4RjTREREZI+Bmug6AaoAhJi6AQAEuQVekal4ZkY8EvpGSlwZEREROSMGaqJazBYRq7ccx+XkUIgW6/8eqnYZ6NnVV+LKiIiIyFkxUBNVqwnT3x/JAIzesORHAwBMoqnBcamJiIiIGKiJYB1nes3WJNsNiHKZgL8Mv8824seuCz9BV2tcaiIiIqIaDNREADb+7zR2H0oHgOobEG/B2IHdbLMnVpmN+IJXqYmIiKgeDNTk8T774Ty2/XABACAIwF8fHoShfaw3IE7ucee1q9Tnf0RBeaFkdRIREZFzYqAmj7b7YBo2fH3a9v7x3/XFiP5RtvdBPhpM6DoaAGC0mPDJb1+3dYlERETk5BioyWP9+tsVvLk1yfb+kQndMSGhY53t7u1+J3yVPgCAn1IPIL04q61KJCIiIhfAQE0e6bcLBXjloyOwiNb3k0Z0wu/HxNW7rZ/KF/f1HA8AECHi4xPb26hKIiIicgUM1ORxLmfrsOz9gzCaLACA0YPa47FJvSEIQoP7jO96G0LUQQCAY1dO4nTeuTaplYiIiJwfAzV5lEJ9JV549wAqDCYAwC09IvDUgwMgkzUcpgHAS67Eg73vsb3/KPlziKLYqrUSERGRa2CgJo9RaTBh2XsHUKCrBAB0iwnCMzNugULeuP8NRsYOQbTGOvrHhcJUHMw83mq1EhERketgoCaPYLGIeH3zMVzI1AEAwoN8sGj2YHh7KRrdhkwmw8N977W933ziC5gsZkeXSkRERC6GgZo8wgc7TuPX364AAHxUCjz32K0I8vducjsD2vVGz7CuAIArpXn4/tI+h9ZJREREroeBmtze7oNptolbZALwzIxbENsuoFltCYKAh/vdZ3u/9dQOVBorHVInERERuSYGanJrJy7kY82nybb3f7q3DwZ1j2hRm11DOuLW9gMBALpKPban7GpRe0REROTaGKjJbeVcLcOKxMMwVw82fc+ITrh7eCeHtD2172TIZXIAwFcp3yG3NN8h7RIREZHrYaAmt2QwmrEi8TBKK4wArMPjPTapt8Pab+cfjrvjxgCwTkm+MWmbw9omIiIi18JATW5HFEW89WkyLmVbR/SIDPXF3x4eBPlNxppuqvt7TkCgt7Uv9pGsZCRdOe3Q9omIiMg1MFCT29n5ayq+P5IBAFB5yfHszMHw9VE6/Dg+Sm880u93tveJxz+ByWxy+HGIiIjIuTFQk1tJSSvEO9t/s71/6vcDmj2iR2MMj41HXIi1X3Z2SS52Xvix1Y5FREREzomBmtxGcYkB/9p4GCaz9SbEySM7Y8SAqFY9pkyQYfbA30OAtTvJ1pP/Q3GFrlWPSURERM6FgZrcgtlswcsfHsHV6mnFe3UKwcyJPdvk2J2CY3F7p2EAgApTJT4+8UWbHJeIiIicAwM1uYUPvzmD3y4WAACCA7zxzPRboJC33dd7ap9JUCt9AAA/pv6K81cvt9mxiYiISFoM1OTyjp/Ns82EKJcJ+MeMeAQFNH1a8ZYI8PbH73tPtL1//9gWWCyWNq2BiIiIpMFATS6tuMSAf28+Zns/c2JP9OgYLEkt47qMQnRAOwDAxcI0fHvxZ0nqICIiorbFQE0uSxRFrNpyHEUlBgDAwO7hmDSis2T1KGRyPDZoqu39xye2o6C8ULJ6iIiIqG0wUJPL+mrvJRw5kwsACPRXYf5DAyBz8OQtTdUzvKvtBsVKkwHvHf0vRFGUtCYiIiJqXQzU5JIuZemw4etrMxP+30MDEeTftv2mG/JIv/ugqZ5B8Wj2bziYeVziioiIiKg1uWSgTklJwZQpU6DVahEeHo74+Hhs2rSp0funpqZCLpdDq9XW+/j4449b5bjkGJUGE1756AhMZutNf/eO6oyB3cMlruoaPy9fzB74e9v7949tQWlVmYQVERERUWtyuUB9/PhxDB48GGazGSkpKcjNzcXcuXMxY8YMLF26tNHtREdHIycnp97HtGnTWu241HLvfnkSmXmlAIBOURrMuKuHxBXVdWv7gRgU2QcAUFypx6bk7dIWRERERK3GpQK1KIqYNWsWACAxMRGBgYEQBAEzZ87Eww8/jGXLliEpKcltjkt1/fpbNnYdSAMAqLzkePqRQVAq5BJXVZcgCHhs0EPwVqgAAHsu7cPpvHMSV0VEREStwaUC9d69e5GcnIy77roLGo3Gbt3UqVNhsViwZs0atzku2dOVGvDWpyds7/98bx+0D/eXsKIbC1UHY1rfe23v1x3ZhCqzUbqCiIiIqFW4VKDes2cPAGDQoEF11tUs2717t9scl+y98/lvKC61DpE3pJcWdwyOkbiimxvXeSS6hnQEAFwpycNnp3dIXBERERE5mksF6jNnzgAAoqKi6qwLCwuDUqlEWloaKioqbtpWeXk5FixYgF69eiEiIgKdOnXCtGnTcPjw4VY9LjXP/hPZ+DkpCwDg56PEnCn9IAjSDpHXGDKZDH++5WHIZdZuKdvPfIsLV1OlLYqIiIgcyqUCdXFxMQDA19e3zjpBEKBWq+22u5GioiJotVrs27cP2dnZ2L59OzIyMjB06FAkJiY69LgGgwF6vd7uQY2nKzXg7W21unrc1wfBbTy1eEvEBEbh/p4TAAAW0YI3Dm6AwVQlcVVERETkKC4VqB0lOjoaV65cwd///ncEBQVBLpejb9+++OKLL+Dn54c5c+YgNzfXYcdbsWIFNBqN7REdHe2wtj3B9V09Rg1sL3FFTXdvj/HoHBwLwNr1Y9OJzyWuiIiIiBzFpQJ1YGAgAKCsrO6YvqIoory83G67hsjlcoSGhtZZHhwcjNtvvx0VFRXYseNaX9eWHnfhwoXQ6XS2R0ZGxg3ro2tctavH9RQyOZ4cMhNKuRIAsPP8jziRc0biqoiIiMgRXCpQ9+hhHW84Kyurzrr8/HwYjUbExsbCx8en2ceIjIwEAFy5csVhx1WpVAgICLB70M25eleP60UGaPFI3/ts79869AEnfCEiInIDLhWox4wZAwA4evRonXU1y8aOHXvTdhITE+u9+RAAsrOzAQDh4ddm3nPUcalp3KGrx/Xu7DoKfSK6AwAKK4rx/rFPJK6IiIiIWsqlAvWIESPQr18/7NixAzqdzm7d5s2bIZPJMHfuXNsyURTr7V6RmJiIDRs21FleXFyMH3/8EV5eXhg/fnyzj0std/h0jlt09bieTJDh8cHToVZa/5qxL+0Qfs2o+4saERERuQ6XCtSCIGDDhg22mQt1Oh1EUURiYiI2bdqExYsXo3///rbt582bh5iYGDz11FN12lq/fj3effddVFVZR1u4ePEipkyZguLiYrzyyito3/7a1dCmHpdaprLKhLWf/2Z7/8d7Xburx/VC1cF4bOBDtvfrj2xGUYXuBnsQERGRM3OpQA0AAwYMwOHDhyEIAuLi4hAeHo4333wTGzduxPPPP2+3bXR0NNRqdZ1RNdatW4fFixfj3XffRceOHREcHIyEhAT4+flhz549mDdvXouOSy3zyXfnkFdovdGzb5dQ3DbI9bt6XG94bDxujR4IACitKsMbBzbAYrFIXBURERE1hyCKoih1EZ5Gr9dDo9FAp9PxBsXrZOSWYN5rP8BkFqGQy/DG30Y79fTiLVFiKMXfdi23XZ1+oNfdeKD3RImrIiIioqZyuSvU5L5EUcRb25JhMlt/x7v/ti5uG6YBwF/lh6dufczWN/zTUztwMjdF4qqIiIioqRioyWn8cDQDJy9eBQBoQ9R44I44iStqfT3Du+LB3vcAAESIWH1gA4orOZMmERGRK2GgJqdQUl6F9786ZXv/l9/1hUopl7CitnNvjzvRN8I61nlxpZ79qYmIiFwMAzU5hY3/Ow1dqXXElWF9IzGoe4TEFbUdmSDDk7fORJC3BgDwW24KPjuzU+KqiIiIqLEYqElyKWmF2HUgDQDgo5Ljj/f2lriitqfxDsC8obNt/am3nvoap/LOSVwVERERNQYDNUnKbBHtphd/eHwPhGiaP3W8K+sVHocHellH+RBFEat+fY/jUxMREbkABmqS1PeH03EpyxoaO0YGYOKwjhJXJK3f9RiPPhHdAFj7U7/6yzoYzUaJqyIiIqIbYaAmyZRXGvHhN2ds7/94bx/I5Z79lZTJZJh362yEqIMAAOevXsZ7R/8LDhdPRETkvDw7vZCkPv3+PIpKDACAoX3aoU/nUIkrcg4a7wA8PewvUMqVAIDvL+/Hrgs/SVwVERERNYSBmiSRV1iO7T9dBAAo5AJmTuwpcUXOpVNwDB6Pn257n3h8K29SJCIiclIM1CSJjf87DaPJOtbyPSM6IzLUT+KKnM/w2HhM6j4OAGARLXj9l3eQV3ZV4qqIiIjoegzU1OZSUgvxc1IWACDA1wsPesCMiM01rc9k9Ndar96XVJXhlX1rUWkySFwVERER1cZATW3KYhHx7hcnbe8fHt8dvj5KCStybjKZDE8NfQzt/MIBAGnFmVhzcCMsImdSJCIichYM1NSmfj6eibPpRQCAGK0/7hwSK3FFzs/XS42nR/wFPgpvAMDBzOP4KOkziasiIiKiGgzU1GYqq0zY+L/TtveP3dPb44fJa6z2Ae0wP+ExyATr+fr63B7sOPe9xFURERERwEBNbeiLny+iQFcJABjUPRwDu4dLXJFrGdCuN/44aKrt/cbjn+Jg5nEJKyIiIiKAgZraSEl5FT7/4QIAQCYT8Nik3hJX5JrGdB6O3/WcAAAQIWL1gQ04W3BR4qqIiIg8GwM1tYnPfriAskoTAGDMLdGIjvCXuCLX9WDvezCywxAAgNFsxMt730Z2Sa7EVREREXkuBmpqdUX6Sny59xIAQCGX4aFx3SSuyLUJgoC/3PII+kR0B2AdTm/FT29CV6mXuDIiIiLPxEBNre6T786hymgGANyV0AHhQWqJK3J9CrkCf034E2I0UQCA3LICvPjTGyitKpO4MiIiIs/DQE2tKrewHDsPpAIAvL3kmDKmq7QFuRG1lw8WjpyLEHUQACC1OBMrfl6DCmOlxJURERF5FgZqalX//fYsTGYRADBpZGcE+XtLXJF7CVEH4Z+jn4LGOwAAcP7qZby8721UmaokroyIiMhzMFBTq8nILcH3R9IBAL4+Stw3uovEFbmnSP8I/HPUPPh6WbvSnMo7h9f2vwOj2ShxZURERJ6BgZpazaZdKbBYL07j/tu6wI9TjLeamMAoLBr5pG02xeNXTmHVgfdhtpglroyIiMj9MVBTq7iQWYxfkrMBAIH+KtwzvJPEFbm/LiEd8MyIOfCSW39xOZSZhLcOfQCLxSJxZURERO6NgZpaxUffnLG9/v2YOHirFBJW4zl6hnfF34b9BQqZ9XzvTTuENw9t5JVqIiKiVsRATQ535nIhjqbkAQDCgnwwfmisxBV5lv7temL+0McgF6z/e+9LO4RVv74PE0M1ERFRq2CgJof7ZM852+uHxnaDUiGXsBrPNLh9f/x12J9sV6oPZB7D67/wRkUiIqLWwEBNDnUhsxhHzlinwQ4L8sHtt0RLXJHnuiWqH/4+/C9QVvepPpJ9Ai/vW8sh9YiIiByMgZoc6pPvrl2dnnJ7Vyjk/IpJqX+7Xlg4Yg5Uci8AQHLOaazYuwaVnPyFiIjIYZh2yGHScvT49bcrAIDgABXuiI+RuCICgN4R3bFo1LUh9U7lncPyn95AiaFU4sqIiIjcAwM1Ocyne87bXt83ugu8lOw77Sy6h3XBP0c/BV+lDwDg3NVLeG7Pa8gvuypxZURERK6PgZoc4kpBGX4+ngkA8Fd7YfytHaQtiOroEtIBS277PwRWT1OeVZKDxd+9gtSiDIkrIyIicm0M1OQQn35/3jYr4uRRnTjutJPqEBSN5WOeRjv/cABAUaUOS75/HSdzUySujIiIyHUxUFOL5RdV4Psj6QAAX28FJg7jrIjOLNwvFMvGPI2uIR0BABWmSrz485vYl3ZY4sqIiIhcEwM1tdhnP56HyWy9PH338E7w9VFKXBHdTIDKD8+Nno9bIvsCAMwWM1YfeB+fn94JURQlro6IiMi1MFBTixSVVOLbA2kAAJWXHJNG8Oq0q1ApvPDXYX/CHZ2G25Zt/u0LvHFgA8eqJiIiagIGamqRL366iCqTBQAwYWgHaPxUEldETSGXyfHHW6bhoT6TbMv2pR/Gkh9eR2F5sXSFERERuRCXDNQpKSmYMmUKtFotwsPDER8fj02bNjV6/8zMTCxfvhy33HILQkJCEBgYiC5duuCJJ57AlStX6t1n9OjRCA4OhlarrfOYPHmyoz6aSymrMGLH/ssAAKVChntHdZa4ImoOQRDwu54T8Ldhf4ZKYf2F6GJhGhbu/hcuXE2VtjgiIiIX4HKB+vjx4xg8eDDMZjNSUlKQm5uLuXPnYsaMGVi6dGmj2ujZsydWr16Nl156Cfn5+SgsLMQbb7yBzZs3o3///rh48WK9+3322WfIycmp8/jiiy8c+Aldx7cH01BhMAMAxsTHIETjI3FF1BKD2/fH8jF/Q5g6GEDNCCCv4efUgxJXRkRE5NxcKlCLoohZs2YBABITExEYGAhBEDBz5kw8/PDDWLZsGZKSkm7ajsViwYoVKzBu3DjIZDLIZDJMmDABL7zwAvLy8rBo0aJW/iSuz2y24Mu9l2zvJ49k32l3EBvYHivG/gM9wroAAIwWE948mIj3j26B0WyUuDoiIiLn5FKBeu/evUhOTsZdd90FjUZjt27q1KmwWCxYs2bNTdt55plncM8999RZPmLECADA/v37HVOwG9t/4goKiisAAPE9I9A+3F/iishRArz98c9RT2FMrZsVd174Ec99/xryOLMiERFRHS4VqPfs2QMAGDRoUJ11Nct2795903b++c9/Ijw8vM7yqirryAYhISEtKdPtiaKIz3+6YHs/eST7TrsbhVyBP90yDX8cNA1KmXWSnouFaXjm25dwJOuExNURERE5F5cK1GfOnAEAREVF1VkXFhYGpVKJtLQ0VFRUNKv9mivTU6dOrXf9li1bMHz4cERHR0Or1WLUqFFYt24dzGbzDds1GAzQ6/V2D1d2JrUQ5zOKAQAdIwPQt0uotAVRqxAEAWO7jMCyMU8jwtf637isqhwv73sbHyV/DpPlxt97IiIiT+FSgbq4uBgA4OvrW2edIAhQq9V22zVFZWUl1qxZg+7du+OJJ56od5vU1FSsW7cOGRkZSE5Oxm233YY5c+Zg4sSJMBob7l+6YsUKaDQa2yM6OrrJ9TmT7T9du2nz3lGdIQiChNVQa+sUHIOV457F4Kj+tmVfpnyL579/HXmlBdIVRkRE5CRcKlC3pr///e/Iz8/H1q1bbcG8tq1bt+Lrr79Gr169AAARERFYunQppk2bhp07d+Ltt99usO2FCxdCp9PZHhkZGa32OVpbztUyHDxpHVowOECFEf3bS1wRtQW1lw/+OuxPeLT/FMgF64+Ns1cv4W+7luPHy79ydkUiIvJoLhWoAwMDAQBlZWV11omiiPLycrvtGmvlypXYsGEDduzYgd69e9e7TVhYGORyeZ3lDzzwAABg27ZtDbavUqkQEBBg93BVX+29BEt1drp7WCcoFS71FaIWEAQBd3cbg+dv/yvCfa33GVSaDHjr0Ad4bf87KDGUSlwhERGRNFwqDfXo0QMAkJWVVWddfn4+jEYjYmNj4ePT+PGQX3nlFaxYsQK7du3Crbfe2uSaIiMjAaDBCWHcSWmFEbsPWacZ91LKMX5oB2kLIknEhXbCy3cuwuiOQ23LDmUm4W87lyPpymkJKyMiIpKGSwXqMWPGAACOHj1aZ13NsrFjxza6veXLl2PlypXYs2cPEhISbMuPHDliG/EDAJKSkrB+/fp628jOzgaAekcNcTffHqg9kUs0Any9JK6IpKJW+mDO4Bn467A/wc/Lek9DUaUOL/38BtYf+RjlxubdGExEROSKXCpQjxgxAv369cOOHTug0+ns1m3evBkymQxz5861LRNFscH+yosXL8Ybb7yBH374oc4wfPHx8bagDFgD9TPPPIPKyso67Xz++ecAgEmTJjX7c7kCk9mCr/bVnsiFQ+URMKT9ALw6fjH6aXvalu2+uBcLvnmBw+sREZHHcKlALQgCNmzYYJsxUafTQRRFJCYmYtOmTVi8eDH69+9v237evHmIiYnBU089ZdfO008/jRdffBG33347tm3bhqVLl9o96lNUVITp06fbupuUlpZi5cqV+OCDDzBs2DDMmzevtT62U9h/Its2kcvgnlpEhflJXBE5i2CfQDw78gnMHvggVAoVAKCwohgv73sb//n1PegqXXuYSCIioptRSF1AUw0YMACHDx/GokWLEBcXB4vFgtjYWGzcuBGPPPKI3bbR0dFQq9V2w9QVFxfj1VdfBQD897//bdQx77//fgiCgG3btmHkyJEoKSlBZWUlunXrhpdffhlPPvkkvLzcu/vDV7WmGb93FK9Okz1BEDC+62gMiuyDd458jOQca1/q/elHcCLnDGYOeAAjYgdziEUiInJLgsjxrtqcXq+HRqOBTqdziRE/LmfrMO+1HwEAsVp/vPG32xiMqEGiKGJv2iEkHt+K0qprI/L0DOuK2QMfRExg3YmZiIiIXJlLdfkgaezYn2p7PSGhI8M03ZAgCBjZYQhen/AcEqKv3Z9wOv88/v7tS0g89gnKq3jTIhERuQ8Garqh8kojfjxqvbHT20uO2wZxIhdqnEDvAMxP+AP+MWKObepyi2jBjvM/4KlvluKnywc4IQwREbkFBmq6oR+OZqKyyjpU3uhB0VB7KyWuiFzNwMg+eG3Cc/h973uglFu/P7pKPdYc2oh/7nkVZwsu3qQFIiIi58Y+1BJwlT7UoijiyVd/QFpOCQBg9V9Ho2OkRuKqyJXllV3FB8c/xaGsJLvlQ9oPwLS+96Kdv/uP505ERO7H5Ub5oLZz+nKhLUx3jw1imKYWC/cNwd+G/xlJV05j4/GtyCrJAQAczDyOI1nJGNt5JKb0ugsB3v4SV0pERNR47PJBDfrmupsRiRylf7ueeHX8Yvxx0DRovK1/pTGLFuy88COe/N9z+PTU/3jjIhERuQx2+ZCAK3T50JUaMPOFb2EyW+Cv9kLic+PgpZRLXRa5oUpjJb48+x2+StkNg7nKttzXS417ut2BCV1vg4/SW8IKiYiIboxXqKleuw+lw2S2AADuGBzDME2txlvpjd/3nojVd7+AOzoNh0yw/lgqqyrHf3/7Ek98vRjbz+xCpbFS4kqJiIjqxyvUEnD2K9QWi4g/rfgOuYXlAIB3Ft6BdqG+EldFniKnNB/bTu3Az2kH7YbVC1D54a642zGuy0j4efH7SEREzoOBWgLOHqiPnMnF8+8eAAAM7BaO5/80VOKKyBNll+Ri26kd2Jd2GCKu/ZjyVqgwtvMI3N1tDIJ9AqUrkIiIqBoDtQScPVC/8N4BHD6dCwBYNGswbu3dTuKKyJNl6q9g26kd2J9x1O6KtUKmwMjYwZjUfSwiA7QSVkhERJ6OgVoCzhyo8wrL8YeXdkMUgVCNN95dNBZyObvak/RySvPxVcpu/Hj5VxgtJrt1A9r1woSut6GvtoetDzYREVFbYaCWgDMH6o93pWDzt2cBANPu7I6p47pJXBGRveIKHXac/wG7LvyEiutuVIz0j8D4rqMxusOt8ObIIERE1EYYqCXgrIHaYhHxx5d2I6+oAjIBeG/xOIQG+khdFlG9yqsq8N2lfdh1/kfklxfarfNRemNUh1txR6fhiAmMkqhCIiLyFAzUEnDWQH3iQj4Wvb0fAG9GJNdhsVhwJPsEvjn/A07lnauzPi6kE+7oPBxDowdBpfCSoEIiInJ3nHqcbL47lG57fUd8jISVEDWeTCbD4Pb9Mbh9f6QVZ+Kbcz9gb/phGM1GAMC5q5dw7uolJB7fiuGx8RjdYSg6B8dCEASJKyciInfBK9QScMYr1OWVRkxfugtVRjN8fZT4YMmdnMyFXFZpVRn2pR3Gdxf3IV2XVWd9lL8WIzsMwYgOgxGqDpagQiIicicM1BJwxkC960Aa3tyaBAC4K6EDHr+/n7QFETmAKIo4f/Uy9lz6BfvTj9hNbQ4AAgT0Co/DyA5DEB/VD75eaokqJSIiV8YuHwQA2HO4VnePwezuQe5BEATEhXZCXGgnPNp/Cg5kHsNPqQdxJv88AECEiJN5Z3Ey7yzkMjn6aXsiIXoQbonqC7WSN+QSEVHj8Aq1BJztCnVmXgkeX/k9ACBW6483/nYb+5eSW8srLcDPaYfwc+oB5JTm11mvlCnQv10v3Np+IAZE9uJU50REdEMM1BJwtkC98X+n8en31it2j03qhXtHdZG4IqK2UdMlZH/GURzIOIbCiuI628gFGXqGxyE+qh9uierLPtdERFQHA7UEnClQmy0iZi/7FoX6SshlAhKfuxOB/ipJayKSgkW04FzBJexPP4oDmcdQXKmvd7tOQTEYGNkHA9r1QuegWMhknJmRiMjTMVBLwJkC9ZEzuXj+3QMAgCG9tFg8e4ik9RA5A4vFgpSCiziSlYxDWUnIK7ta73b+Kj/01/bEgHa90U/bA/4qvzaulIiInAEDtQScKVD/64PD+CU5GwCwaNZg3Nq7naT1EDkbURSRocvGoaxkHM5MwuXijHq3EyCgY1A0+kR0R5+I7uge2hlenEiGiMgjMFBLwFkCdUl5FWYs3QWT2YJAPxU2PDcOCjn/fE10I4XlxTh+5SSO55zCbzkpqDBV1rudUqZAt9DO6B3RDT3CuqBLcAco5co2rpaIiNoCh83zYD8dy4TJbAEAjB7UnmGaqBGC1YEY03k4xnQeDpPZhJSCizh+5SR+y01BanGmbTujxWQbkg+wBuyuIR3RM7wreoZ1RZfgDvBWekv1MYiIyIEYqD3Yd4c51ThRSyjkCvSO6IbeEd0AALpKPU7lncOJ3BT8lnMG+eWFtm2NFhNO55/H6eoxsGWCDLGBUegW0tk2VnaYOphDVhIRuSB2+ZCAM3T5SM/RY+4rPwAAukQH4t/zR0lSB5G7EkURuWUFOJN33hak8xu4ubFGoHcAugR3QOfgWHQJ6YDOQbHwU3EMbCIiZ8cr1B7q5+NZtte3DWwvYSVE7kkQBGj9wqD1C8NtnRIAAAVlhTidfx5n8i/g3NVLyNRdgYhr1zSKK/U4kn0CR7JP2JZF+IWhc1AMOgRFo2NQNDoGRiPA27/NPw8RETWMgdoDiaJoC9QyARjeP0riiog8Q6hvMEb6DsHIDtbhKcurKnC+8DLOFlzCuYJLuFCYinJjhd0+uaX5yC3Nx/6Mo7ZlIT5B6BDUHrGBUYjRtEdMYCTa+YVDLpO36echIiIrBmoPdD6jGFeulgEA+nQJRXAAb4wikoLaywf9tD3RT9sTgHVymZzSfFy8moaLham4UJiGy8UZMJqNdvtdrSjC1YoiHM3+zbZMKVOgfUA7RAdGon1AO+tD0w7h6hBOPkNE1MoYqD3QT8evjUQwcgC7exA5C5kgQ6R/BCL9IzCiw2AAgMliRrY+B5eLMnC5OAOXizKQWpyBCqP9cH1Gi8m6/rpxspVyJaL8IxAZoK1uOxyR/hFo5x8BH44yQkTkEAzUHsZsEbEvydrdQyGXIaEPJ3IhcmYKmRwxgVGICYzCKNwKwHolO7/sKtJ12UgrzkK6LgvpxVm4UpqH6+8zN5qNSC3OtBvSr0aQtwZa/zBEVPf1rnlE+IXB10vdJp+PiMgdMFB7mJMXC1CoNwAABnUPh5+aM7kRuRqZIENEdfCNj+pnW15lqkJ2SS4y9VeqHznI0uXgSmkeLKKlTjtFlToUVepwJv9CnXW+Sh+E+4YizC8E4eoQhPuFIsw3BGHqYISqg6H28mnVz0hE5EoYqD1M7dE9RnF0DyK34qXwQoegaHQIirZbbjKbkFOWjyslecjW5yK7JBdXSqzPekNpvW2VGSvq7UJSw0fpjTB1CELVQQipfgT7BFpf+wQi2CeQE9cQkcdgoPYgRpMZv5zIBgD4qOSI7xkhcUVE1BYUcoXtRkVcN6hPubECuaUFyC3NR071I7c0H3llV3G1vKjeK9sAUGGstHY10WXVux4AfBTeCPLRWB/eGtvrQO+A6of1ta+XmhPaEJFLY6D2IMdS8lBWYR0tYEivdvD24n9+Ik+nVvpYx7e+7qo2AJgtZlwtL0JeWQHyyq6ioLwQ+WWFKCgvREF5Ea6WF8FkMTXYdoWpEhUllcguyb1hDXKZHBqVPzQqfwR4Vz+r/BDgXf2s8oN/9SPAyw9qLx/IBI5cQkTOwyUTVUpKChYvXox9+/bBYrEgNjYW8+fPx8MPP9ykdj788EOsXr0aaWlpkMlkGDFiBF588UXExcW16nGlwu4eRNQUcpkc4X6hCPcLrXe9RbRAV1mCq+VFKKwoRmFFMa6WF+FqRTEKy4usfbQrdKg0GW54HLPFbNu/MQRBgJ+XL/y9fOHn5Qs/1bXXvl5q+Hmp4atUX3vtpYav0gdqLzW85MqmngYioptyuUB9/PhxjBo1CmPGjEFKSgo0Gg02btyIGTNm4Pz581i6dGmj2lm8eDFWrFiB9957D48++ih0Oh1mzZqF+Ph47N27F3379m2V40qlwmDCgVM5AAB/tRf6x4VJXBERuTqZILN147iRCmMliip1KCwvRnGlHrpKPYqve+grS6A3lMDcQBeT2kRRRImhFCUN9P++EaVMAbWXGmqlN9QKH6i9vOGj9LG+Vlpf+yhV8FH4wEfpbX0oVPBWWF97K1TwUajgpfDiVXIishHE68dYcmKiKGLAgAG4dOkSMjIyoNFc+yE+Y8YMbNq0CUePHkX//v1v2M7Ro0cRHx+P6dOnY+PGjbblOp0O0dHRiIuLw+HDh219+hx13Bp6vR4ajQY6nQ4BAQGNPwEt8OOxTLy2yTrT2oShHTBnSr+b7EFE1LYsogXlVRXQGUqgqw7YJYYylFSVWgN3VRlKqpeVVpWhpKqsznjcbUWAAC+FF7wVKnjLrc8qhar62QsquRdUChVUcqX1WeEFL3nNci94yZW2ZV5yZfXD/rVSrmBoJ3IRLnWFeu/evUhOTsaDDz5oF2oBYOrUqfjwww+xZs0arF+//obtrFmzBqIoYurUqXbLNRoNJkyYgE8++QT79+/HsGHDHHpcKf10rPZkLpxqnIicj0yQwU9l7cIRFaBt1D4mixmlVWUoNZShtKocZcZylFWVW5dVWV+XGytQZqxAeVW53XOFsRIimndNSYQIg8kAg8kAXbNaaByFTAGlXGEN2zIFlHIllDXP1a8VcgW8ZEoo5ArrOtt6BRQyORQyhbWdmme5AnJBDoW8Zt21Z7lQe5n1tVwmh1wmh0KQ273mDJxE17hUoN6zZw8AYNCgQXXW1SzbvXt3i9v55JNPsHv3blugdtRxpaIvq8Lxs3kAgFCNN3p2DJG4IiIix1DI5LZRQ5rKIlpgMFWhvDpclxsrrDdSGqsfpmvPlUaD9dlkQKXtvTVQV5oMqDRX1Zki3hFMFhNMFpNkV+JvRIAAmUxmF7Tlguzas3BtmUxW/V6QQVazTJBBLqt+rg7oNcuty669lslqvRZkkAlCI99fWy7A+lq4bp0gCLZ1MkGwrbcuEyDYvRbs2rl+uVCzP2r2Q/UyGWQQAAGQQVb9LADV2157xrV2aj/XtFNrWe1tSXouFajPnDkDAIiKqnuFNSwsDEqlEmlpaaioqICPT/2TDpSXlyM9PR1eXl4IC6vbj7im7ZSUFIceV0r7T2TDbLFehRkxoD1kMv7PR0QkE2S2ftKOYLaYYTBVwWCusl69NlfVeV9lMlqfbeuMqDJXocpsrH5Yg3mV2Qij2WR7XWWxPpvMRlRZTK0S3ptKhAizxQwzzIBZ6mpIsAV2a0CvCeLW1zXra7/GddvcZHvrQVBrS/s26tm3epNr9dXap3bN195fa6f2Pte33zuiO2b0v98Rp81hXCpQFxcXAwB8fX3rrBMEAWq1GjqdDsXFxQ0G25o21Or6p9WtabuoqMhhxzUYDDAYrt3lrtfr6z12a/npOLt7EBG1NrlMDrWXD9Ro/QsromgNs1UWI4xmI4wWE0xmE4zVYdv6bILJYobJYoTJYrYGdIsJZou5+sq32XYFvGa5dZ0ZJvHaNjXLzaLZ7r1JNMNiscAk1t7GArN47XXNeotogdnC1N2aRIiACJirn91ZhJ/zDazgUoHaVa1YsQLPP/+8JMc2mswwGq13zUeF+aFz1I3vxiciIucnCAIUcmv/aSid7y+jDbFcF7ItFos1bNcK4hZRtK4TLTBXr6//ITbyvQix5hnXL7e+FlH9XHsdrOut62q3Ub1cFGFBzT72219rF4AownLdPhABC0Sgnv2sWVi01Wdbct2xRdSEaLG6Let70e41aq2vvU+tba2Hs9ve/nXNJtf2q//9dftWv7etq1VvzbKaNmzram17bTv79gHrL6/OxqUCdWBgIACgrKyszjpRFFFeXm633Y3aqNn2ejVtBwUFOey4CxcuxIIFC2zv9Xo9oqPrTqLQGpQKOV59aiSy80txVV/JvlZERCQZmSCDTC6DwrXiB9FNudQ3ukePHgCArKy6U93m5+fDaDQiNjb2hv2Y1Wo1YmJikJ6ejvz8/Dr9qGva7t69u8OOq1KpoFKpbvLpWldkmB8iw/wkrYGIiIjIHbnUmDdjxowBYB1H+no1y8aOHevwdhx1XCIiIiJyPy4VqEeMGIF+/fphx44d0OnsR/7cvHkzZDIZ5s6da1smiiIyMjLqtDN37lwIgoDNmzfbLdfpdPjmm28wcOBAJCQkNPu4REREROQ5XCpQC4KADRs2QBRFzJo1CzqdDqIoIjExEZs2bcLixYvtZiucN28eYmJi8NRTT9m1M2jQIDz77LP46KOPsHHjRoiiCJ1Oh5kzZwIA3n//fbu+xk09LhERERF5DpcK1AAwYMAA27TgcXFxCA8Px5tvvomNGzfWGUkjOjoaarW63hsAly9fjg0bNmD16tWIiIhAXFwc5HI5Dh8+jH796k7L3ZTjEhEREZHnEERRdPPRCp2PXq+HRqOBTqdDQEDTZ/ciIiIiIufhcleoiYiIiIicCQM1EREREVELMFATEREREbUAAzURERERUQswUBMRERERtQADNRERERFRCzBQExERERG1AAM1EREREVELMFATEREREbUAAzURERERUQswUBMRERERtQADNRERERFRCyikLsATiaIIANDr9RJXQkREREQ34u/vD0EQbrgNA7UESkpKAADR0dESV0JEREREN6LT6RAQEHDDbQSx5nIptRmLxYLs7OxG/cbjKHq9HtHR0cjIyLjpl8KT8Lw0jOemfjwv9eN5aRjPTf14XhrGc1M/qc4Lr1A7KZlMhvbt20ty7ICAAP7PWQ+el4bx3NSP56V+PC8N47mpH89Lw3hu6ueM54U3JRIRERERtQADNRERERFRCzBQewiVSoUlS5ZApVJJXYpT4XlpGM9N/Xhe6sfz0jCem/rxvDSM56Z+znxeeFMiEREREVEL8Ao1EREREVELMFATEREREbUAAzURERERUQswULuolJQUTJkyBVqtFuHh4YiPj8emTZua3M6HH36I+Ph4hIeHQ6vV4oEHHsC5c+daoeK209Jzk5qaCrlcDq1WW+/j448/bsXqW9fJkyeRkJAAQRCQmprarDbc8TsDtOzcuNt3xmAwYMuWLZg4cSK0Wi1CQkIQFhaGu+++G999912T23On74yjzo27fWdMJhM++eQTzJo1C7169UJkZCTCwsLQv39/rFy5EqWlpU1qb8eOHRg5ciTCw8MRERGBCRMm4NChQ61Ufety5LkRBKHB78zKlStb8VO0jSeeeAKCIGDmzJlN3lfynzMiuZxjx46J/v7+4r333isWFRWJFotF3LBhgyiTycQlS5Y0up1FixaJMplM3LBhg2ixWMSioiLx3nvvFQMCAsTk5OTW+wCtyBHn5vLly2JsbGyr1tnWKioqxGeffVYMDg4WQ0NDRQDi5cuXm9yOO35nHHFu3O078+STT4oAxIULF4p6vV4URVFMS0sT77jjDhGA+MYbbzS6LXf7zjjq3Ljbd+bKlSsiAHHAgAHib7/9JoqiKBqNRvG9994TZTKZOGTIENFoNDaqrfXr14sAxOXLl4tGo1EsLy8XH3/8cVGpVIrffvtta36MVuHIc+POsW337t2iIAgiAPHRRx9t0r7O8HPGff/LuCmLxSL269dP9Pf3F4uLi+3WTZ8+XZTJZOLx48dv2s6RI0dEQRDEGTNm2C0vLi4W/f39xUGDBokWi8WRpbc6R50bd/uHThRF8c9//rM4efJkMSMjQxw1alSzQqM7fmdE0THnxt2+M3PnzhVHjBhRZ3l+fr7o4+MjqlQqsaio6KbtuON3xlHnxt2+MzWh8dixY3XWTZkyRQQgfvfddzdtJzMzU/Tx8RFHjhxpt9xoNIqdOnUSo6KixPLycofV3RYcdW5E0X0DdVFRkdi+fXtx+vTpTQ7UzvJzhl0+XMzevXuRnJyMu+66CxqNxm7d1KlTYbFYsGbNmpu2s2bNGoiiiKlTp9ot12g0mDBhAo4ePYr9+/c7tPbW5qhz444WLlyI7du3t2jKe3f8zgCOOTfuZvz48Vi0aFGd5aGhoejevTsMBgOOHTt203bc8TvjqHPjbkJCQrBv3z7079+/zrrY2FgAgE6nu2k77733HioqKup8ZxQKBaZMmYKsrCx89tlnDqm5rTjq3LizJ554Av3798fs2bObvK+z/JxhoHYxe/bsAQAMGjSozrqaZbt3726zdpyJO34mR6n5od0S7np+HXFu3M3EiRNx55131ruuqqoKgDUk3Iw7fmccdW7cjVKpxLBhwyAIQp11hw4dgre3N4YMGXLTdtzxO+Ooc+Outm3bhl27dmH9+vXN2t9ZvjMM1C7mzJkzAICoqKg668LCwqBUKpGWloaKiooG2ygvL0d6ejq8vLwQFhZWZ31N2ykpKQ6qum044tzUKC8vx4IFC9CrVy9ERESgU6dOmDZtGg4fPuzwul2Bu35nHMkTvjMFBQU4f/48evbsib59+95wW0/7zjTl3NRw5++MxWLB5cuXMXfuXBw7dgwbNmyo92fz9W70c9xdvjPNPTc1lixZgn79+kGr1SImJgb33nuvy/2SUSM3Nxd/+ctf8Pbbb0Or1TZ5f2f6OcNA7WKKi4sBAL6+vnXWCYIAtVptt92N2qjZ9no1bRcVFTW/UAk44tzUKCoqglarxb59+5CdnY3t27cjIyMDQ4cORWJiogOrdg3u+p1xJE/4zqxevRomkwmrV6+u92pbbZ72nWnKuanhrt+Zr7/+GoGBgejUqRN27tyJjz/+GA899FCj9r3Rz3F3+M605NzUKC8vx65du3DlyhV8//33kMlkGDduHJYtW9ZKVbeeP/zhDxg3bhymTJnSrP2d6eeMotWPQORioqOjceXKFYSGhtqW9e3bF1988QU6deqEOXPmYMKECYiIiJCwSnImnvCdOXDgAP71r3/hhRdewJgxY6Qux6k059y483dm4sSJ0Ov1yM/Px3//+19MnToVt99+OzZv3gw/Pz+py5NUS8/NlStX7K7kdunSBVu2bEH37t2xZMkSTJo0Cf369WvNj+Aw7733Ho4dO4aTJ09KXYpD8Aq1iwkMDAQAlJWV1VkniiLKy8vttrtRGzXbXq+m7aCgoOYXKgFHnBsAkMvldv/I1QgODsbtt9+OiooK7Nixo8X1uhJ3/c44irt/Z06fPo2JEydi3rx59d6QVx9P+c4059wA7v+dAaxd7Z588kksW7YMX3/9NRYuXHjTfW70c9xdvjNA884NgHq7RSiVSkyePBmiKLrMDZupqalYsGAB3nvvvRb993SmnzMM1C6mR48eAICsrKw66/Lz82E0GhEbGwsfH58G21Cr1YiJiUFVVRXy8/PrrK9pu3v37g6qum044tzcTGRkJADrVQJP4q7fmbbg6t+ZkydP4vbbb8fs2bPx6quvNno/T/jONPfc3Iyrf2euN3HiRADA9u3bb7rtjX6Ou8N35npNOTc34mrfmR07dsBisWDmzJl2k9P87ne/AwBs2bLFtuyLL75osB1n+jnDQO1iav6cePTo0TrrapaNHTu2zdpxJo76TImJiQ3eFJSdnQ0ACA8Pb26ZLssdvzOO4q7fmePHj+O2227DX/7yF7z88su25ampqbbPdSPu/J1p6blxt+/Mjz/+iK1bt9a7rqZ/69WrV2/ajjt+Zxx1brZv397gzYeu9p2ZM2cOSkpKkJOTY/eoucL+4IMP2pZNnjz5hm05zXem1Ue6Jodq6uQlFotFTE9Pr9POzQZCHzhwoMtNuOCoczNq1Cjx8ccfr7O8qKhIDAoKEr28vMSMjAyH199WbjZ5iSd9Z67X3HPjjt+ZQ4cOiUFBQeLKlSvrrHv00UftZh71tO+MI86Nu31nlixZIvbo0UM0m8111r3//vsiAHH06NF2yzMyMkSTyVRn2c0mdikrK3P8B2hFjjo3jz76qDhhwoQ6bRiNRrFr164iAPHgwYOOLb6N/fDDDw1O7OLsP2cYqF3QsWPHRD8/P/G+++4Ti4uL7abXfu655+y2feKJJ0QA4rx58+q0UzNVZ2JiomixWMTi4mLx3nvvFf39/cWkpKS2+jgO5YhzM2rUKFGhUIjr168XDQaDKIqieOHCBXHMmDGiIAjiqlWr2uzztIabhUZP+87U1txz427fmV9++UUMCAgQu3fvLi5ZsqTOo1+/fnah0ZO+M446N+72nVmyZIkIQJw9e7aYlZUliqIoGgwGcevWrWJQUJAYEBAgHjlyxLb9q6++KgIQJ02aVKetd955RwQgvvjii6LJZBIrKirExx9/XFQoFOKuXbva7DM5iqPOzaOPPmo7L6WlpaIoimJ2drb40EMPiQDEBQsWtN2HaiU3CtTO/nOGgdpFnTlzRvzd734nhoeHi6GhoeKgQYPEDz/8sM52K1euFNVqtfjKK6/U287GjRvFgQMHimFhYWJ4eLh4//33iykpKa1dfqtq6blJSUkRly5dKg4ZMkSMjIwUg4KCxPDwcHHy5Mni999/31Yfw6EOHjwoRkREiBEREaJSqRQBiKGhoWJERESdH8Ke9p1xxLlxt+/M5MmTRQA3fNQOjZ70nXHUuXG370x+fr741ltviWPHjhU7dOgghoWFiX5+fmK3bt3EuXPn1vklddOmTaKvr6/45JNP1tve//73P3H48OFiaGioGBYWJt55553igQMH2uCTOJ6jzk16err4yiuviCNHjhTbt28vhoSEiMHBweLYsWPFbdu2teEncrzly5eLERERYlBQkAhA9Pb2FiMiIsQ+ffrYtnH2nzOCKIpia3UnISIiIiJyd7wpkYiIiIioBRioiYiIiIhagIGaiIiIiKgFGKiJiIiIiFqAgZqIiIiIqAUYqImIiIiIWoCBmoiIiIioBRioiYiIiIhagIGaiIiIiKgFGKiJiIiIiFqAgZqIiOwsWrQIgiDgl19+adJ+6enpEAQBb731VitVRkTknBioiYjIztGjRyGTydC/f/8m7bdt2zbIZDL87ne/a53CiIiclCCKoih1EURE5DzCw8MRHByMlJSUJu03bNgwyGQy7N27t5UqIyJyTrxCTUREAID58+dDEATk5+fj7NmzEATB9jhz5swN983Ozsavv/6KKVOmtFG1RETOQyF1AURE5BwGDx6MBx98EFu2bMH48eMxZMgQAIAgCIiLi7vhvp9//jkANKq7h8ViwVtvvYX33nsPubm5EAQBWq0WY8eOxezZs296LCIiZ8NATUREAIBp06YhKysLW7ZswRNPPIG777670ft++umniI+PR3R09E23fe655/Dvf/8bO3bswKhRowAAu3fvxpQpU+Dt7Y2lS5c29yMQEUmCXT6IiMjm2LFjAIABAwY0ep/8/Hzs3bu30d09PvvsM3Tr1s0WpgFg7NixWLBgAUJDQ5tWMBGRE2CgJiIim+PHjyM8PByRkZGN3mf79u0wm824//77G7V9hw4dcPz4cSxZsgQ5OTm25UuWLMETTzzR5JqJiKTGUT6IiAgAUFZWhoCAANxxxx3YtWtXo/e78847kZeXh+PHjzdq+0uXLmHatGk4ePAgZDIZhgwZgnvvvRczZsyAVqttbvlERJLhFWoiIgIAJCUlwWKxNKm7R1FREX744Ycmje7RqVMnHDhwAEeOHMEzzzyDgoICPPPMM+jSpQu+/PLL5pRORCQpBmoiIgIAnDhxAgCaNKHLl19+CaPR2OjuHgBgNpsBAIMGDcJLL72Ec+fOYfv27TAajZg3b16TaiYicgYM1EREBAC4evUqACAoKKjR+3z66afo2bMnunfv3uh9OnfujAMHDtgtmzx5Mnr16oXCwsJGt0NE5CwYqImICMC1kT3mzZuHhQsXYunSpfjpp58a3L6kpMQ23F1TPfvss0hLSwMAiKKIL774AqdOncKjjz7avOKJiCTEcaiJiAgAcPfdd2PFihV455138Oqrr8JkMqFbt24Nbv/111/DYDA0qbsHALz11lv4+OOPMXbsWJSWlkIURbRv3x7/+c9/8Kc//amlH4OIqM1xlA8iImqW+++/H7/99hvOnTsndSlERJLiFWoiImqWW2+9lV00iIjAK9RERERERC3CmxKJiIiIiFqAgZqIiIiIqAUYqImIiIiIWoCBmoiIiIioBRioiYiIiIhagIGaiIiIiKgFGKiJiIiIiFqAgZqIiIiIqAUYqImIiIiIWoCBmoiIiIioBf4fj84NQUH/qAYAAAAASUVORK5CYII=\n",
            "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 the generalized $\\alpha$-scheme we find the approximation of the time-discrete rate as\n",
        "\n",
        "$$\\displaystyle \\frac{y_{n+1}- y_n}{t_{n+1} - t_n} = [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1})]$$\n",
        "\n",
        "This can be re-cast in a residual form, i.e based on an approximation error:\n",
        "\n",
        "$$r(y_{n+1}) = \\displaystyle \\frac{y_{n+1}- y_n}{t_{n+1} - t_n} - [(1-\\alpha) f(y_n) + \\alpha f(y_{n+1})]$$\n",
        "\n",
        "where we would like to find the solution $y_{n+1}$ such that $r(y_{n+1}) = 0$.\n",
        "\n",
        "Due to the non-linear dependence of $r$ on the sought solution $y_{n+1}$, the equation cannot be solved directly. We thus use the concepts of a Taylor series expansion truncated after the first member to arrive at a recursive expression:\n",
        "\n",
        "$$r(y_{n+1}) \\approx r(y_n) + \\left. \\frac{\\partial r}{\\partial y} \\right|_n \\underbrace{(y_{n+1} - y_n)}_{\\displaystyle \\Delta y_{n+1}} \\overset{!}{=} 0$$\n",
        "\n",
        "Thus, with the help of the linearization we proceed from the known solution to the new one in a linear step\n",
        "\n",
        "$$\n",
        "    y_{n+1} = y_n - \\left( \\left. \\frac{\\partial r}{\\partial y} \\right|_n \\right)^{-1} r(y_n)\n",
        "$$\n",
        "\n",
        "Because this is only a linear step, we make an error (truncation error of the Taylor series expansion). Thus, we take the new guess for $y_{n+1}$ to repeat the above step recursively until we arrive at a satisfactory solution (linearization always around the newly found state):\n",
        "\n",
        "$$\n",
        "    y_{n+1}^{i+1} = y_{n+1}^{i} - \\left( \\left. \\frac{\\partial r}{\\partial y} \\right|_{i} \\right)^{-1} r(y_{n+1}^i)\n",
        "$$\n",
        "\n",
        "To perform this update, we need two incredients: a residual evaluation and a Jacobian evaluation (aka \"linear slope\")."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 43,
      "metadata": {},
      "outputs": [],
      "source": [
        "def residual(y_i,y_n,alpha,dt):\n",
        "    return (y_i - y_n)/dt - ((1.-alpha) * dy_dt(y_n) + alpha * dy_dt(y_i))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 44,
      "metadata": {},
      "outputs": [],
      "source": [
        "def df_dy():\n",
        "    return -k"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 45,
      "metadata": {},
      "outputs": [],
      "source": [
        "def jacobian(alpha,dt):\n",
        "    return 1./dt - alpha * df_dy()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note, that the Jacobian is constant for a given time step size (usually it is not!!). Thus, in the present case a distinction between Newton-Raphson, modified Newton, or initial stiffness methods is immaterial. In other words, an initial stiffness approach can be chosen here without any loss of convergence quality. We will demonstrate the full Newton-Raphson scheme below.\n",
        "\n",
        "The Newton update function is then as follows (returns $\\Delta y_{n+1}^{i+1} = y_{n+1}^{i+1} - y_{n+1}^i$):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 46,
      "metadata": {},
      "outputs": [],
      "source": [
        "def Newton_Update(res,Jac):\n",
        "    return -res/Jac"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, for every time step we advance the system by applying the generalized midpoint Newton-Raphson update rule. Each new time step value then serves as the initial condition for the following linearization step.\n",
        "\n",
        "We will *recursively* apply the Newton update until the residual is below a specified tolerance:\n",
        "\n",
        "$\\left| r_{n+1}^{i+1}\\right| < \\epsilon_\\text{abs}$\n",
        "\n",
        "with a suitably chosen absolute tolerance set here to $10^{-6}$.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 47,
      "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",
        "            #evaluate residual\n",
        "            res = residual(y_old,y[-1],alpha,dt)\n",
        "            #if residual is below tolerance, above maximum iterations, stop iterations\n",
        "            if (np.abs(res) < absolute_tolerance or i > max_iter): \n",
        "                break\n",
        "            #evaluate Jacobian\n",
        "            Jac = jacobian(alpha,dt)\n",
        "            #perform linear step\n",
        "            y_new = y_old + Newton_Update(res,Jac)\n",
        "            #update counter\n",
        "            i += 1\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\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 48,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "4a2f024f37574ff580f6377489e82ebb",
              "version_major": 2,
              "version_minor": 0
            },
            "text/plain": [
              "interactive(children=(BoundedFloatText(value=0.2, description=r'$\\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.2,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": {},
      "source": [
        "Tasks:\n",
        "\n",
        "- Compare the number of iterations with the Picard implementation (Exercise 3)\n",
        "- Investigate the influence of the time-step size\n",
        "- Repeat the analysis for $\\alpha = 1$ (Backward Euler)\n",
        "- How does the solution change for the modified NR method and the initial stiffness method?\n",
        "- Repeat the analysis for a function with a non-constant Jacobian (one example is in lecture notes, but you may choose any other)\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.11.4"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}