04c_Erddruck_Newton.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "![Logo TUBAF](https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg)\n",
    "\n",
    "Exercise material of the MSc-level course **Numerical Methods in Geotechnical Engineering**.\n",
    "Held at Technische Universität Bergakademie Freiberg.\n",
    "\n",
    "Comments to:\n",
    "\n",
    "*Prof. Dr. Thomas Nagel  \n",
    "Chair of Soil Mechanics and Foundation Engineering  \n",
    "Geotechnical Institute  \n",
    "Technische Universität Bergakademie Freiberg.*\n",
    "\n",
    "https://tu-freiberg.de/en/soilmechanics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise on Earth pressure and the Newton-Raphson scheme\n",
    "\n",
    "![Erddruck.png](https://raw.githubusercontent.com/nagelt/Teaching_Scripts/master/Images/Erddruck.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometry \n",
    "Let the retaining wall have a vertical height $H$ and the origin be placed at the bottom of the slip surface which has length $L$. Then, the problem is characterized by the following points in the $x-z$ plane:\n",
    "\n",
    "* Origin: $O = [0,0]$\n",
    "* Top left: $A = [-H \\tan \\alpha, H]$\n",
    "* Top right: $B = [L\\cos \\vartheta, L \\sin \\vartheta]$\n",
    "\n",
    "The length of the slip surface can be found via the following trigonometric relation:\n",
    "\n",
    "\\begin{align}\n",
    "    L \\sin \\vartheta &= (L \\cos \\vartheta + H \\tan\\alpha) \\tan \\beta\n",
    "\\end{align}\n",
    "\n",
    "Thus we arrive at\n",
    "\n",
    "\\begin{align}\n",
    "    L &= \\frac{H(1 + \\tan \\alpha \\tan \\beta)}{\\sin \\vartheta - \\cos \\vartheta \\tan \\beta}\n",
    "\\end{align}\n",
    "\n",
    "The area of the triangle (i.e. the volume of the slip body per meter) can also be calculated. We use Gauß' equation for that:\n",
    "\n",
    "$$\n",
    "    A = \\frac{1}{2} \\left| (x_2 y_3 - x_3 y_2) - (x_1 y_3 - x_3 y_1) + (x_1 y_2 - x_2 y_1) \\right|\n",
    "$$\n",
    "\n",
    "Let Python do the math for us ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np #numerical methods\n",
    "#import sympy as sp #symbolic operations\n",
    "import matplotlib.pyplot as plt #plotting\n",
    "import sympy as sp\n",
    "sp.init_printing(use_latex='mathjax') #makes sympy output look nice\n",
    "\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "#Some plot settings\n",
    "import plot_functions.plot_settings\n",
    "%run plot_functions/earth_pressure_plots.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define our geometrical properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "H,alpha,beta,theta = sp.symbols('H alpha beta vartheta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{H \\left(\\tan{\\left(\\alpha \\right)} \\tan{\\left(\\beta \\right)} + 1\\right)}{\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}}$"
      ],
      "text/plain": [
       "       H⋅(tan(α)⋅tan(β) + 1)        \n",
       "────────────────────────────────────\n",
       "sin(vartheta) - cos(vartheta)⋅tan(β)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = H*(1+sp.tan(alpha)*sp.tan(beta))/(sp.sin(theta) - sp.cos(theta) * sp.tan(beta))\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the function that calculates the area of the triangle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def triangle_area_Gauss(A,B,C):\n",
    "    term1 = (B[0]*C[1] - C[0] * B[1])\n",
    "    term2 = (A[0]*C[1] - C[0] * A[1])\n",
    "    term3 = (A[0]*B[1] - B[0] * A[1])\n",
    "    return -(term1 - term2 + term3)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The points we defined above can now be inserted into this equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{H^{2} \\left(\\sin{\\left(\\vartheta \\right)} \\tan{\\left(\\alpha \\right)} + \\cos{\\left(\\vartheta \\right)}\\right) \\left(\\tan{\\left(\\alpha \\right)} \\tan{\\left(\\beta \\right)} + 1\\right)}{2 \\left(\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}\\right)}$"
      ],
      "text/plain": [
       " 2                                                           \n",
       "H ⋅(sin(vartheta)⋅tan(α) + cos(vartheta))⋅(tan(α)⋅tan(β) + 1)\n",
       "─────────────────────────────────────────────────────────────\n",
       "           2⋅(sin(vartheta) - cos(vartheta)⋅tan(β))          "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PO = [0,0]\n",
    "PA = [- H * sp.tan(alpha), H]\n",
    "PB = [L * sp.cos(theta), L*sp.sin(theta)]\n",
    "A = triangle_area_Gauss(PO,PA,PB).simplify()\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".Thus, the area of our triangle is given as:\n",
    "\n",
    "$$\n",
    "    A = \\frac{H^2 (\\sin \\vartheta \\tan \\alpha + \\cos \\vartheta) (\\tan \\alpha \\tan \\beta + 1)}{2(\\sin \\vartheta - \\cos \\vartheta \\tan \\beta)}\n",
    "$$\n",
    "\n",
    "Now we can write the force balance. For that purpose, we introduce the material's specific weight $\\gamma$, friction angle $\\varphi$ and the wall friction angle $\\delta_\\text{a}$ as well as the cohesion $c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "c, delta, phi, gamma = sp.symbols('c delta_a varphi gamma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the downward force due to gravity as well as the cohesive force along the slip surface. The active earth pressure force $E_\\text{a}$ and the reaction force $Q$ in the slip line are introduced as unknowns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ea, Q, G, C = sp.symbols('E_a Q G C')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write the vertical and horizontal equilibrium conditions:\n",
    "\n",
    "\\begin{align}\n",
    "    \\text{horizontal:}\\quad  0 &= C \\cos\\vartheta - Q \\sin (\\vartheta - \\varphi) + E_\\text{a} \\cos (\\alpha + \\delta_\\text{a})\n",
    "    \\\\\n",
    "    \\text{vertikal:}\\quad  0 &= C \\sin\\vartheta + Q \\cos (\\vartheta - \\varphi) + E_\\text{a} \\sin (\\alpha + \\delta_\\text{a}) - G\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fh = sp.Eq(C * sp.cos(theta) - Q * sp.sin(theta - phi) + Ea * sp.cos(alpha + delta),0)\n",
    "Fv = sp.Eq(C * sp.sin(theta) + Q * sp.cos(theta - phi) + Ea * sp.sin(alpha + delta) - G,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C \\cos{\\left(\\vartheta \\right)} + E_{a} \\cos{\\left(\\alpha + \\delta_{a} \\right)} + Q \\sin{\\left(\\varphi - \\vartheta \\right)} = 0$"
      ],
      "text/plain": [
       "C⋅cos(vartheta) + Eₐ⋅cos(α + δₐ) + Q⋅sin(varphi - vartheta) = 0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C \\sin{\\left(\\vartheta \\right)} + E_{a} \\sin{\\left(\\alpha + \\delta_{a} \\right)} - G + Q \\cos{\\left(\\varphi - \\vartheta \\right)} = 0$"
      ],
      "text/plain": [
       "C⋅sin(vartheta) + Eₐ⋅sin(α + δₐ) - G + Q⋅cos(varphi - vartheta) = 0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = sp.solve([Fh,Fv],(Ea,Q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{C \\sin{\\left(\\alpha + \\delta_{a} - \\vartheta \\right)} + G \\cos{\\left(\\alpha + \\delta_{a} \\right)}}{\\cos{\\left(\\alpha + \\delta_{a} + \\varphi - \\vartheta \\right)}}$"
      ],
      "text/plain": [
       "C⋅sin(α + δₐ - vartheta) + G⋅cos(α + δₐ)\n",
       "────────────────────────────────────────\n",
       "    cos(α + δₐ + varphi - vartheta)     "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get(Q).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".... for $Q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{C \\cos{\\left(\\varphi \\right)} + G \\sin{\\left(\\varphi - \\vartheta \\right)}}{\\cos{\\left(\\alpha + \\delta_{a} + \\varphi - \\vartheta \\right)}}$"
      ],
      "text/plain": [
       "-(C⋅cos(varphi) + G⋅sin(varphi - vartheta)) \n",
       "────────────────────────────────────────────\n",
       "      cos(α + δₐ + varphi - vartheta)       "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get(Ea).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... for $E_\\text{a}$. If we substitute the result for $G$, we start to see the issue ...."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{C \\left(\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}\\right) \\cos{\\left(\\varphi \\right)} + \\frac{H^{2} \\gamma \\left(\\sin{\\left(\\vartheta \\right)} \\tan{\\left(\\alpha \\right)} + \\cos{\\left(\\vartheta \\right)}\\right) \\left(\\tan{\\left(\\alpha \\right)} \\tan{\\left(\\beta \\right)} + 1\\right) \\sin{\\left(\\varphi - \\vartheta \\right)}}{2}}{\\left(\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}\\right) \\cos{\\left(\\alpha + \\delta_{a} + \\varphi - \\vartheta \\right)}}$"
      ],
      "text/plain": [
       " ⎛                                                        2                   \n",
       " ⎜                                                       H ⋅γ⋅(sin(vartheta)⋅t\n",
       "-⎜C⋅(sin(vartheta) - cos(vartheta)⋅tan(β))⋅cos(varphi) + ─────────────────────\n",
       " ⎝                                                                            \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                      (sin(vartheta) - cos(vartheta)⋅tan(β))⋅c\n",
       "\n",
       "                                                                 ⎞ \n",
       "an(α) + cos(vartheta))⋅(tan(α)⋅tan(β) + 1)⋅sin(varphi - vartheta)⎟ \n",
       "─────────────────────────────────────────────────────────────────⎟ \n",
       "                     2                                           ⎠ \n",
       "───────────────────────────────────────────────────────────────────\n",
       "os(α + δₐ + varphi - vartheta)                                     "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get(Ea).subs(G,A*gamma).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we would have to find the maximum with respect to $\\vartheta$ in order to identify the active earth pressure and the critical slip angle $\\vartheta_\\text{a}$. The equations get very complex and we decide to go for a numerical solution.\n",
    "\n",
    "First, we define the basic geometry of our problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 10. * np.pi/180.\n",
    "beta = 20. * np.pi/180.\n",
    "H = 4. #m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slip_length(theta):\n",
    "    L = H * (1 + np.tan(alpha) * np.tan(beta))/(np.sin(theta) - np.cos(theta)*np.tan(beta))\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def area(theta):\n",
    "    term1 = (np.tan(alpha)*np.tan(beta) + 1) * (np.sin(theta)*np.tan(alpha) + np.cos(theta))\n",
    "    term2 = np.sin(theta) - np.cos(theta) * np.tan(beta)\n",
    "    return H*H * term1 / term2 / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define the material properties and the interaction with the wall:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 19.e3 #N/m³\n",
    "phi = np.deg2rad(30.)\n",
    "c = 10.e3 #N/m²\n",
    "delta = np.deg2rad(10.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cohesive_force(theta):\n",
    "    return c * slip_length(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def weight(theta):\n",
    "    return gamma*area(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have all the ingredients to calculate the active earth pressure given by the equation above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ea(theta):\n",
    "    term1 = weight(theta) * np.sin(theta-phi) - cohesive_force(theta) * np.cos(phi)\n",
    "    term2 = np.cos(alpha+delta+phi-theta)\n",
    "    return term1/term2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvgAAAHvCAYAAAAhGatdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABIK0lEQVR4nO3deXzcV33v/9fRvtuS5X2Pl8RO4sSJE7KSQBIg0FAoW9iXUri9vW1ZCr9Ce7mlpS2UtrS0tCW0kLImZQskkBBCSMhuTBY7Tux43x3bkmXt65zfHzNyZMWLJEv6amZez8djHl/p+z0z+shf2X7rzFlCjBFJkiRJuaEg6QIkSZIkjR4DviRJkpRDDPiSJElSDjHgS5IkSTnEgC9JkiTlEAO+JEmSlEOKki4gl7zqVa+Kd911V9JlSJIkKfeFE12wB38UHTp0KOkSJEmSlOcM+JIkSVIOMeBLkiRJOcSAL0mSJOUQA74kSZKUQwz4kiRJUg4x4EuSJEk5xIAvSZIk5RADviRJkpRDDPiSJElSDjHgS5IkSTnEgC9JkiTlEAO+JEmSlEMM+JIkSVIOMeBLkiRJOcSAL0mSJOUQA74kSZKUQ4qSLkCSpHwUY6SnN0V3b4re3hS9felHz4CP+1KRVCq+6BhjJAJEjn4cY/p1Q4AAhBDSH2eOBSFQWBgoCIGCgkBhQaCwoIDCwkBRYQHFRQUUFhRQVBQoLiygqKiA4qJCCgtCcn9IkkbEgC9J0imkUpHO7l5aO3poG/Do6OqlvauXjs7MMfNxZ3cvnd19dPf0Hf24K/N5d2+KnswxGxQVBoqLCikuKqCkqICS4kJKSwopPXosorSkkLKSQspKiygrKaS8tIiykiLKSwspLy2mvKyIirIiKkqLqCgrpqIsfb3AXx6kMWHAlyTllVQq0trRw5HWLppau2hu7U4f27ppae+mpa2b5syxpb2b1vYe2jt7SMXRr6WosICS4nTveVFh/yPdo16Y+biwoOBoj3t/73t/zzz099inz0G6Jz8S08f+3v1UJBUjqRT0pVJH3w3oS0X6+lL09sWj7xr09qXfUTj67kJfpLevl46u0f3eQ4CK0iIqK0qoKiumsryYyvIiKsuLqa4ooaoicywvpqqihOqKYmoqS6mpLKGspPDo9y/pxQz4kqSc0JeKNLV0cqipg0NHOjnc3Mnhli4ON3fSOODjI23dpEaQ1stLC6k8GkSLj/ZEl2d6pctLi44+yvp7tEuKKC0d2NtdSElxf2944YTvwY4x0tsX6ento7snlXkHIv1uRFfPC8f+dyg6u3vpzLyT0dndl35H4+g7HD20d/bS3tlLR1cPHV19tHX20tbZy4Fh1lVcVEBNZcnRx6TKUiZVlzKpqoTJVaXUVJYyuaqUydWl1FaXUlZq3FF+8SdekjThxRhpae/h4OF2DhzueOHY1E5DUyeHjnRwuLlzyL3sleXFTKosYVImBNZkPq7O9BRXV5ZQU1FCdWW6B7myvJiiwvxblyKEQHFRoLiogIqy0X3tvlSko7PnhWFPnelja3sPLe09tHZ009Lek3kXpZuWth6a27tpbu2iuzdFw5FOGo50DulrlZUUUltdlg78NaXUVpdRV1NGXU0pdTXl1NaUUldTRk1lie8MKCcY8CVJE0JPb4rnG9vY39DO8w1t7GtoZ39DG883tvN8YxsdXX2nfI3JVaVMmVzGlJpy6iaVUVddSm1NGbWZY11NGZOqSikuyr+wPtEUFgSqKkqoqigZ9nM7u3vTQ6raujnSlg79Ta3dHGntyjy6aWrtpKm1m6bmTjq7+9jX0Ma+hraTvm5RYQF1k8qon1TGlEnlTMkc6yeXUT+5nKmTy5lcXebEY014BnxJ0rhJpSIHDrez92Abew62svdQK3sPtrH3UCsHGttP2gNfXlrEtNpyptZWMLW2nGm1FUydXE595lFXU0pxUeH4fTNKTFlJepLutNqKU7aNMdLe2UtTa9cLw7ZaOmk8khm61dxFQ2YYV1tHDwca2znQ2H7C1yssCNRNKjv6szettuLoz+W0zM+lQ4KUNH8CJUmjrq8vxd5Dbex6viXzaGXXgRZ2H2ilu+f4PfEFAabVVTBzSgUzplRmHhXMqKtk+pQKqsqLHT6hYQshHJ03MXtq1Unbdnb30tjcSUNTJw1H0nM5Go50pOd1NHVwqKmTptYuDh7u4ODhjhO+TnVFCdOnVDC9toLpdRXpj+teePiLqMaaAV+SNGIxRhqbO9m+r5kd+5ozxxZ2HWih5wTLQE6uLmXOtCpmT61iVn0Vs6dWMmtqFTOmGHyUrLKSImbVp38uT6S7p4+GI+nJ3Aeb0nNBDjS2c/BwBwcyc0Na2tMrMG3e1fSi54cAU2rKmFFfycwp6V9eZ06pZGZ9JTPrq6gqLx7D71D5woAvSRqSvlRk78FWtu45wra9R9iSOR5p7T5u+2m15cybUcOcaVXMnV7N3GnVzJ1eNaIx19JEUVJcmAnjlce9nkpFmlq7eL4hPXfk+cPtmY/b2d/YfnSVp0NHOnl6S8OLnl9TWcLM+kpmZQJ//y/As+orqSgz/GtoDPiSpBdJpSL7GtrYtKuJTbsOs2lnE1v3HqGr+8XDa6rKi1kwq4YFM2qYP7OG+TNqmD+z2jCivFRQEDIr9JSxbGHdi6739qU4eLiD/Q1tmUd7egLwofQk4Oa2bprbutm44/CLnltXU8qsqel3v2ZPrWLOtCrmTKtmWl2FE391DAO+JIkjrV1s3HmYDdsb2bjjMFt2N9HW2fuidlNryzlj1iTOmJ15zJrE1Npyx8ZLQ1RUWHDCdwD6h7ztO9TG3kNt7D3YesyxsbmLxuauF/X8FxUWMHtqJXOmVb/wjtn0amZPq6K02GFv+SjEOAZb8+WpVatWxTVr1iRdhiSdVCoV2bG/mWe2NbJhRzrQ7zv04uUD62rKWDJ3cuZRy6I5k5hUVZpAxZL6UpFDTR3sOdjKngOt7D3Yyu4Drew+0MKhE+wHEALMqKvMBP4q5s2oYd6MdPg3+OeEE/asGPBHkQFf0kTU05ti864m1m9rYP3WBp7d1vCi3vnSkkKWzJ3MmfNqOXN+HUvnTWbKpPKEKpY0HO2dPew92MbuAy3sOtB6dPWqvYfajrtrc0GA6VMqmT+jOj2kLjOsbvbUKgrzcEO3LGbAHw8GfEkTQU9vik27DrNu8yHWbj7Ehu2NdA9a0WZabTnLF07hrAV1nDm/loUza/yPXcoxPb0p9h7KBP79Lex4voWd+5vZc/D4wb+osIC506tYMLMm85jEwlk11NaM8jbGGi0G/PFgwJeUhFQqsnXPEZ7adJC1mw/xzLYGOgdNhp07vZqzz5jC2QvrWH7GlCFtECQpN/X09rHnYBs79zezY38LO/Y1s21f8wk3+JpcVcqCWTUsnDWJM2bVsHD2JObY2z8RGPDHgwFf0nhpONLBExsP8MTGgzzx3EFa2o9dqnLu9CrOXVTPiiVTOeeMKY6dl3RK7Z097NzfwvbMnhbb9h5h+75m2o8z4b6kqID5M2s4Y/YkFmUm3S+YNcmx/ePLgD8eDPiSxkpPb4pntjWw5tnneXzjAXbubznm+tTacs5fMpXzlkzl3MX11PmWuqRREGPk+cZ2tu1tZvveI2zb18yWPUeO29tfUBCYO62KRXMms2jOJBbPmczCWZMoL3XRxjFiwB8PBnxJo+lwSye/efZ5fv3s8zyx8SAdXS/0opWVFHLu4npWLp3GBWdNY1Z9pUtVSho3rR09bNuT3vBu654mtu45wq4DrS8a2x8CzJlWxZK5tSyek16Va+Fse/pHiQF/PBjwJZ2OGCO7D7TyyLp9PLZ+H8/tbDrm+tzp1Vy0bDqrlk3nrAV1FBc5/lXSxNHV08f2vUfYvPsIW3Y3sWX3EXbsb6ZvUOgvKAjMn1HNkrm1LJ03maXzapk3vdox/cNnwB8PBnxJw5VKRZ7bdZhH1+3j0af3sefgC+vRFxcVsGJxfTrUL5/B9DonxkrKLt09fWzf18ymXU1s3tXE5t1N7NzfzOBFfEqKC1k8ZxJL59WydF4tZ86rdRO9UzPgjwcDvqShSKUiz25v5MEn9/Dwur00NncdvVZdUcxFy2dw6bkzOX/JVMocuyopx3R29bJ17xGe29nEpl2HeW7nYfY3vHhMf211KWfOT+/Nceb8WpbMmey/iccy4I8HA76kE0mlIht2NPLgU3t56Km9NDa/sPPk1NpyLjlnJpeeM5PlC+t8m1pS3jnS2sWmXU1s2nmYDTsP89yOw7R29BzTpqAgsHBWDcvm13HmgjqWLahjWn738hvwx4MBX9JAMUa27W3ml7/ZxYNP7jlmO/lpdRVcsWIWl583iyVzJ+fzf1CS9CIxRvYeamPjjkY27DjMxh2H2b6v+UWTeOtqSjlrQR3LFkxh+cI6zpg9iaL86SQx4I8HA74kgAON7dz/xG7ue3z3MctZTq0t54rzZnOFoV6Shq2zq5dNu5p4dnsjz25vZOOORlraj+3lLyku5Mx5tSxfWJfZrbuWirLihCoecwb88WDAl/JXe2cPDz61l3vX7GL91oaj56srSrjy/FlcfcFczlpQa6iXpFESY2TPwVae3ZYO/M9sa2TPwdZj2hQEWDh7EmcvnMLZZ0xh+cIpTK7OmY3/DPjjwYAv5ZcYI09vbeCe1Tt5aO1eurr7gHQP0iVnz+CqC+dwwZnT8untYklKVFNLVybsN/DstkY272560TKdc6ZVcfYZUzhnUT3nLprClEnlCVV72gz448GAL+WHQ00d/OLXO/nFr3exr+GFZS3PPmMK1140j8tWzMzlt4QlKWt0dvWycedh1m9tYP3WBjbsOEx3T98xbWbWV3LO0cBfz9TarAn8BvzxYMCXclcqFXniuQPc+fB2fv3M/qNrOE+ZVMY1F83jmovmMqu+KtkiJUkn1dObYsueJtZvaWDdlkM8s63xmF3CAWZOqeTcxfWsyDxqa8oSqvaUDPjjwYAv5Z7DLZ3cs3onP3t0B883ptdpLiwIXHLOTF7xkvmct3QqhQWOq5ekbNTXl2LLniM8fTTwN9DeeWzgnzOtKh32l0xlxeJ6qitKEqr2RQz448GAL+WGGCMbdxzm9ge28vC6vfT2pf+dnFZbzqsuXcC1F82byD06kqQR6g/8azcfYt3mQ6zf1nB0fhVACLBo9iTOWzKVFUumsnxhHWUliW2+ZcAfDwZ8Kbv19KZ46Kk9/PiBrWza1QSkV2C4aPkMXnXpAlaeOc3eeknKIz29KTbtOszazYd4atNBNmxvPNrpA1BUWMDvvGwx77x+WRLlnfA/JPf7lZT3mlq6uPOR7dz58DYOt3QBUF1RzCsvWcD1ly1gWm1FwhVKkpJQXFTA8oXp5TVvvO5MOrt6eWZbI09tOshTmw+ydc8RJldNvGU3DfiS8tbuAy3cdv8W7l2zi57eFADzZ1Rzw5WLuPrCOZQWFyZcoSRpIikrLeKCs6ZxwVnTAGhu66ZgAr6za8CXlHc2bG/kB/dt5tGn9xFjekzlS86ewQ1XnsGKxfVuRiVJGpKaygkz4fYYBnxJeSGViqx59nm+/8tNPLOtEUiPnbzmorm87qpFzJlWnXCFkiSNDgO+pJyWSkUeXreXW3/+HNv3NQNQWV7Mqy9bwA1XnOFqOJKknGPAHySEMA34ArAqc2od8KEY4+7kqpI0XH2pyANP7uF/7tnIrudbAairKeP1Vy/iFS+Z706zkqScZcAfIIRQAvwceA44G4jAV4FfhhBWxhhbk6xP0qn19aW47/Hd/M89z7H3UBsAU2vLeePLl3DtRfMoceKsJCnHGfCP9W5gBfD6GGMvQAjh/wP2AL8PfD7B2iSdRCoVefCpPXzrrg1Hg/2MKRW86ZqlvOzCuRQXFSRcoSRJ48OAf6w3ADtjjFv7T8QY94cQnslcM+BLE0yMkdXr9/PNuzYcHWM/s76SG69bylUr51BYaLCXJOUXA/6xVpAenjPYNuCaca5F0ik89dxBvnHns2zceRiA+snl3HjdmVxz0VyKDPaSpDxlwD9WPfCb45xvBipCCOUxxo6BF0IIHwA+ADBv3ryxr1ASW/cc4Wu3r+fJTQcBmFxVypuuXcKrLlngGHtJUt4z4A/NCXe9iTHeBNwEsGrVqjhuFUl56FBTB9+481l++ZtdxAiVZUW84eVLuOGKMygr9Z8zSZLAgD/YIeB4u91UA+2De+8ljY/2zh6+d+8mfnT/Frp7UxQVBl5z+Rm8+dqlE3YXQUmSkmLAP9Za4KzjnF9Iej18SeOoLxW5+9HtfOtnGzjS2g3AFefN4l2vXs7M+sqEq5MkaWIy4B/rB8CXQwgLYozbAUII04FlwCeSLEzKN+u3NnDTD9exde8RAJYtqON9N5zNWQvqEq5MkqSJzYB/rJuB/wN8LoTwdiAFfJb0Kjr/nmBdUt5oONLBzXc8w32PpzePnlpbzvtuOJvLV8wihBNOh5EkSRkG/AFijN0hhOuALwDPkN7J9mng5e5iK42tnt4+brt/C/9zz3N0dvdRXFTAG162hDe8fDFlJf5TJUnSUPm/5iAxxueBtyVdh5RPntp0kH///lPsOZjegfbSc2fyvhvOZsYUx9lLkjRcBnxJiTnS2sVXb1/PvWt2ATBnWhUfeN25rDxzWsKVSZKUvQz4ksZdjJFf/HonX719PS3tPRQXFfCW65byO1cvobjIHWglSTodBnxJ42rX8y186XtPsX5rAwDnL5nK779xBbPqqxKuTJKk3GDAlzQu+vpS/OC+zXz7Zxvp7UsxqaqE97/2HK66YI6r40iSNIoM+JLG3I59zfzTrU+weVcTANddPI/33nA21RXuQitJ0mgz4EsaM319Kb7/y8185+50r3395HL+8M3nc4GTaCVJGjMGfEljYvu+Zv75lsfZvDu9E+0rL5nP+244m4qy4oQrkyQptxnwJY2qVCpy2/1b+Madz9DbF5laW84fvfl8zl9qr70kSePBgC9p1Bw83ME/3fI4azcfAuBVly7gvb+13F57SZLGkQFf0qh44Ik9fOn7T9HW0cPkqlL+6C3nc9HyGUmXJUlS3jHgSzotbR09fPmHa/nlb3YDcNHy6fzRm1cyubo04cokScpPBnxJI7ZhRyOf/+ZvONDYTklxIe//7XN41SXzXddekqQEGfAlDVuMkR/9ags33/EMfanI4jmT+OjbL2TOtOqkS5MkKe8Z8CUNS0t7N/98yxM8tn4/AL/90kW8+zXLKS4qSLgySZIEBnxJw7BxRyN/9401HDjcQWV5MR+6cSWXnDMz6bIkSdIABnxJp5QekrOVm+9YT18qsnTeZD7+zouYXleRdGmSJGkQA76kk+ro6uWfb3mCh9buBRySI0nSRGfAl3RCew+28tc3r2bn/hYqyor40I0XcOm5DsmRJGkiM+BLOq41zz7P339zDW2dvcyZVsWfvfdiV8mRJCkLGPAlHSOVinz33uf41l0biBEuOWcGH37rBVSUFSddmiRJGgIDvqSj2jt7+KdbnuCRdfsIAd7xqrN40zVLKShw4ypJkrKFAV8SAAca2/nL/3qUHftbqCwr4qNvv5CLls9IuixJkjRMBnxJbNzRyGe+tpqmli5mT63iU7/7EmZNrUq6LEmSNAIGfCnPPfjUHr7w7cfp7k2xYnE9n3j3RVRVlCRdliRJGiEDvpSnYox89xeb+MadzwLwipfM5/ffsIKiQte3lyQpmxnwpTzU09vHv373Ke5ds4sQ4D2vOZvXX72IEJxMK0lStjPgS3mmraOHv/7aatZtOURpSSEffduFbl4lSVIOMeBLeaSxuZO/+MojbNvbTG11KZ/63UtYPHdy0mVJkqRRZMCX8sSeg6186qZHONDYzuyplXz6A5cxva4i6bIkSdIoM+BLeeC5nYf59H8+SnNbN0vnTeZTv3sJk6pKky5LkiSNAQO+lOMe33CAv/3v1XR293HBmdP403dfRHmpf/UlScpV/i8v5bD7Ht/NP33ncfpSkasvnMMfv2Wly2BKkpTjDPhSjrrrke382/efIkZ4/dWLec9rllNQ4DKYkiTlOgO+lIN+/MAWvnLb0wC8+zXLeePLlyRckSRJGi8GfCnHfO/eTfz3T54B4AOvO5cbrjwj4YokSdJ4MuBLOSLGyLd/tpFbfr6REOAP3nger7xkQdJlSZKkcWbAl3JAjJGb73iGH9y3mYIAf3zjBbx81dyky5IkSQkw4EtZLsbITbet444Ht1FYEPiTd1zIFefNTrosSZKUEAO+lMVijHzlR09zx4PbKCos4BPvvoiLz56RdFmSJClBBnwpS/UPy7n9ga0UFRbwZ++9mFXLpiddliRJSpg73khZ6lt3beAH922msCDwp+9aZbiXJEmAAV/KSrf+fCO33vMcBQWBj71zFS85Z2bSJUmSpAnCgC9lme/fu4lv3rWBggAfeesFXL5iVtIlSZKkCcSAL2WRH/1qCzf/5BlCgD96y0quumBO0iVJkqQJxoAvZYmfP7aD//zR0wD8wRvP55qL5iVckSRJmogM+FIWeOzpffzrd58E4AOvO5dXXjI/2YIkSdKEZcCXJrintxzi776xhlSEG687kxuuPCPpkiRJ0gRmwJcmsG17j/CZrz5Gd2+K6y9dwNteeWbSJUmSpAnOgC9NUPsb2vh/Nz1CW2cvl6+YxQd/ZwUhhKTLkiRJE5wBX5qADrd08qkvP8Lhli5WLK7no2+/gMICw70kSTo1A740wbR39vAXX3mUfQ1tLJoziT9778UUFxUmXZYkScoSBnxpAunrS/G5b6xh654jzKyv5C/efykVZcVJlyVJkrKIAV+aIGKMfPm2dTy+4QA1lSV8+vcuZXJ1adJlSZKkLGPAlyaIH/1qK3c+vJ3iogL+/L0vYWZ9ZdIlSZKkLGTAlyaAR5/ex1dvT+9S+6EbV7JsYV3CFUmSpGxlwJcStnlXE3//rd8QI7zjVWfx0pVzki5JkiRlMQO+lKCDhzv4q68+Sld3Hy9fNZc3X7s06ZIkSVKWM+BLCWnv7OGvvvoojc1dnLNoCv/nTee7kZUkSTptBnwpAalU5AvfeZxte5uZPbWST77nYoqL/OsoSZJOn4lCSsCt9zzHo0/vp7K8mE/97iVUV5QkXZIkScoRBnxpnK1+Zj/f/tkGQoCPveNCZk2tSrokSZKUQwz40jjafaCFf/jWbwB45/XLuPCs6QlXJEmSco0BXxon7Z09/PXXVtPe2cvlK2bxxpcvSbokSZKUgwz40jhIpSL/+O3H2X2glfkzqvnjG1e6Yo4kSRoTBnxpHNx6z3M8tj49qfaT772Y8tKipEuSJEk5yoAvjbEXTaqtd1KtJEkaOwZ8aQw939jOP377ccBJtZIkaXwY8KUx0tOb4vPfWENbRw8XL5/hpFpJkjQuDPjSGPn6T59h487DTK0t50NvdVKtJEkaHwZ8aQw89vQ+brt/C4UFgY+/Y5U71UqSpHGTF0t5hBDuA6YB3YMu/WOM8euD2r4d+DgQSP8C9M8xxq+MR53KDQca2/mnW54A4F2vXs5ZC+oSrkiSJOWTvAj4Ga+OMW4/WYMQwo3AV4ErY4yrQwgrgEdCCAUxxi+PR5HKbr19Kf7um2to7ehh1bLpvO6qRUmXJEmS8oxDdDJCCAXA54FbY4yrAWKMa4Gbgb8JIZQmWJ6yxNd/+iwbdxymflIZH37rBRQUOO5ekiSNLwP+Cy4G5gC/HHT+XqAOeNm4V6SssvqZ/fzwvs0UFAQ+/s6LqKl03L0kSRp/+RTwPxJCWB1C2BBC+FUI4b2Drq/IHLcNOr9t0HXpRRqbO/mn72TG3V+/jGULHXcvSZKSkS8BvwnYTLoX/mzgi8C/hxD+fkCb+syxZdBzmzPHKWNZoLJXjJEv3voELe3dnL90Kq+/enHSJUmSpDyWdQE/hHBtCCEO4XFf/3NijK+LMX4xxtgWY+yLMX4P+C/gwyGEeaf6kqeo5wMhhDUhhDUHDx48/W9QWefOR7bzmw0HqCov5kM3rnTcvSRJSlQ2rqLzMLBsCO3aT3H9MeB/AxcBO4FDmfPVg9r1f95wvBeJMd4E3ASwatWqOIS6lEN2H2jhv368HoA/eNN5TJlUnnBFkiQp32VdwI8xtgMbhto+hFAClMcYjwy61Jc5FmaOazPHBYPaLRx0XQLSS2L+w7cfp7unj6svnMMV581OuiRJkqTsG6IzApcB/3Oc8xdmjk9kjquB3cDVg9q9DGgE7huD2pTFbvn5RjbvamJqbTn/6/XOwZYkSRNDPgR8gGtCCK/p/ySEcDXwQeAbMcZNADHGFPAx4MYQwqpMu3OB9wCfjDF2jnPNmsA2bG/ku/c8Rwjw4bdeQGV5cdIlSZIkAVk4RGcEHgc+DnwyhPA3QCXQDXyG9MZWR8UYbwkhFAJfDSFA+hegj2TG2UsAdHT18o/ffpxUhN+5ejHnLqo/9ZMkSZLGSc4H/BhjM/CPmcdQ2n8L+NaYFqWs9l8/fpp9DW0snFXDO64/K+lyJEmSjpEvQ3SkUfH4xgP87NEdFBUW8NG3XUhxUeGpnyRJkjSODPjSEHV09fKl7z4JwNteeSbzZ9YkW5AkSdJxGPClIfr6T5/hwOEOzpg9yd1qJUnShGXAl4bgmW0N/OShbRQWBP74LSspKvSvjiRJmphMKdIpdPf08cVbnyRGeMPLl3DG7ElJlyRJknRCBnzpFL5z90b2HGxl7vQqbrxuadLlSJIknZQBXzqJzbub+MF9mwkB/ujNK101R5IkTXgGfOkEevtSfPHWJ0ilIjdceQZnLahLuiRJkqRTMuBLJ/D9X25i295mptdV8M5XLUu6HEmSpCEx4EvHsftAC7fc/RwAf/im8ykrzflNnyVJUo4w4EuDxBj58g/W0duX4tqL5nHe0qlJlyRJkjRkBnxpkAef2suTmw5SXVHMe35redLlSJIkDYsBXxqgvbOH//zR0wC869XLmVRVmnBFkiRJw2PAlwb4zt0baWzuZOm8ybziJfOTLkeSJGnYDPhSxvZ9zfz4ga0UBPj9N5xHQUFIuiRJkqRhM+BLpCfW/vv3nyKVilx/2UIWz5mcdEmSJEkjYsCXgHvX7OKZbY1MrirlHde75r0kScpeBnzlvdb2br52x3oA3nvD2VSVFydckSRJ0sgZ8JX3vn7nsxxp7ebsM6bwsgvnJF2OJEnSaTHgK69t3tXEXY9sp7Ag8PtvWEEITqyVJEnZzYCvvBVj5D9//DQxwg1XnsH8GTVJlyRJknTaTivghxCmjVYh0nh7ZN0+1m9toKayhLdcd2bS5UiSJI2K0+3BXzMqVUjjrKe37+jE2re98iwn1kqSpJxRdKoGIYTXnuRy2SjWIo2bOx7cxv6GduZOr+JVl7hjrSRJyh2nDPjAD4H7gePNPqwe3XKksXektYtbf74RgPfdcA6FhU5FkSRJuWMoAX8T8Lsxxm2DL4QQdo1+SdLY+s7dG2nr7GXl0qlceJbTSCRJUm4ZStflfwP1J7j2H6NYizTmdu5v5s5HtlMQ4Hdfe47LYkqSpJxzyh78GOPfnuTaX49uOdLY+todz5BKRa6/dAHzZ7ospiRJyj3DHnwcQpg9FoVIY+3xjQdY8+zzVJQV8bZXnpV0OZIkSWNiJLMLnxj1KqQx1teX4r9+/DQAb75mKZOrSxOuSJIkaWyMJOA7aFlZ555f72Tn/ham1VVww5VnJF2OJEnSmBlJwI+jXoU0hrp7+vjO3ellMd/96mWUFBcmXJEkSdLYcQFw5byfPrydhiOdLJxVwxXnOYVEkiTlNgO+clp7Zw/f/cVzALzz+mUUFDjCTJIk5TYDvnLaj361lea2bpYtqGPVsulJlyNJkjTmRhLw3b1WWaG5rZsf3rcZgHe+epmbWkmSpLww7IAfY7xwLAqRRtv3791ER1cvF5w5jXMXnWgzZkmSpNziEB3lpIYjHdzx4FYgPfZekiQpX4xJwA8hfGMsXlcaqlvveY7u3hSXrZjJ4rmTky5HkiRp3BSN9IkhhMXAVcAMYPDC4tecTlHS6dh3qI27H91BQYC3v/KspMuRJEkaVyMK+CGE/w38Cyfe1dbNsJSYb9+9gb5U5OWr5jJvRk3S5UiSJI2rkQ7R+TjwB8DUGGPB4AewdvRKlIZux75m7n98N0WFgbfZey9JkvLQSIfoHIwx/sdJrr9zhK8rnZbv3L2RGOGVlyxgel1F0uVIkiSNu5H24P8whHDuSa5/cISvK43Yzv3NPLxuL0WFBbzpmiVJlyNJkpSIEfXgxxj/JoTwuRDCImAz0D6oyRuAPzzd4qTh+O4vNhEjXPeSeUyZVJ50OZIkSYk4nUm2HztJEyfZalztPdjKr57YTWFB4I0vs/dekiTlr5EO0fkk6WE49U6y1UTwvXs3kYrw8lVzmebYe0mSlMdGOsn2QIzxKye57iRbjZvnG9u5d80uCgK80bH3kiQpz420B//nTrLVRPH9X26iLxV56co5zKqvSrocSZKkRI20B78LuC2E8AROslWCGo508PPHdhICrpwjSZLEyAP+n2eOC09w3Um2Ghc/uG8zvX0pLl8xy11rJUmSGPkQnaeON7nWSbYaT00tXdz1yA4A3nLd0oSrkSRJmhhGGvA/dYrrDs/RmLvt/s109/Rx8fIZLJw1KelyJEmSJoQRBfwY4+2nuP7gyMqRhqa5rZufPrwNsPdekiRpoJH24EuJuv2BrXR09bFy6VSWzqtNuhxJkqQJw4CvrNPZ1ctPHtoKwFuuOzPhaiRJkiYWA76yzi/W7KKlvYcz59WyfGFd0uVIkiRNKAZ8ZZW+VORHv9oCwOuuXkQIIeGKJEmSJpYhB/wQwsvHshBpKFav38++Q21Mq6vg0nNmJl2OJEnShDOcHvy/H7MqpCH64X2bAfjtK8+gsNA3oCRJkgYbzk62S0II9w61cYzRHn+Nqo07Gnl2eyOVZUVce/G8pMuRJEmakIYT8PcDnz7J9ZnAF4DpwFdOpyjpeH54f3rs/asuXUBFWXHC1UiSJE1Mwwn4j8UY7z/ehRDC9cA/A+XA22OM3xmN4qR++xvaeGTtXgoLAjdceUbS5UiSJE1YQx7EHGN8x+BzIYSiEMI/AncAe4ELDfcaC7c/sJVUhJeunM2USeVJlyNJkjRhDacH/xghhMXALcAFwL8DH44xdo9WYVK/1vZu7n5sBwCvv3pxwtVIkiRNbCMK+CGEdwH/CqSAN8cYvzeqVUkD3PXoDjq7+zhvST0LZ01KuhxJkqQJbVjrDIYQqkII3wBuBjYAKw33Gks9vSluf2ArAK+7yt57SZKkUxnORlergCeAt5GeUHtZjHHbCdoedzKuNFwPPLmHxuZO5k6v5sKzpiVdjiRJ0oQ3nCE6D2Xa30Ym6IcQjtcuAPNPuzLlvRgjP/pVemnM11+1iBP8vEmSJGmA4QT8VuDHmY9fdoq2VSMrR3rBhu2H2brnCDWVJVx1wZyky5EkScoKwwn4O2OM7x1KwxDCEyOsZ9hCCAXAnwB/BXwwxnjzeH1tja2fPpweAfaKl8ynpLgw4WokSZKyw3Am2b5ljNqOWAhhHnAv8Fag5BRt3x5CeCqEsDaE8HQI4fdO0O4VIYTVIYR1IYQNIYRPZH6J0Dhqauniwaf2EgJcf+mCpMuRJEnKGkPuwY8xPtf/cQjhrafY0OpfgVecTmFD9FHgq8BO4JcnahRCuDHT7soY4+oQwgrgkRBCQYzxywPaXUF60643xxhvCyHMBR4DaoBPjOH3oUF+vnoHvX0pLl4+g2l1FUmXI0mSlDVG2jP9sRNdCCFcTXrzq/Hw0Rjj10/WINP7/nng1hjjaoAY41rSS33+TQihdEDzzwGPxhhvy7TbBXwB+GgIYdbol6/j6UtF7nxkOwCvvnxBorVIkiRlm5EG/HMzveJHhRCqQwj/AfwCqD3tyoYgxtg7hGYXA3N4cQ//vUAdmQnDIYSZwGUnaFcMvPa0itWQrXlmPwcPdzBzSiUrl7o0piRJ0nCMNOBvBd4fQngnQAjh1cB64B2ke/fXj055o2JF5jh4zf5tg66fO8R2GmM/fXg7ANdftoCCApfGlCRJGo7hrKIz0AdJj03/UQjh3aR7wR8AXhZj3BJCuGO0ChwF9Zljy6DzzZnjlGG2O0YI4QPABwDmzZs38ioFwN5DrTy+8QAlRQVce7F/npIkScM1oh78GON9McYO0sNWeoE7YoxXxxi3ZJr87nBfM4RwbQghDuFx30hqPt6XHI0XiTHeFGNcFWNcNXXq1NF4ybx2Z6b3/qUr51BdcdKFkSRJknQcQ+rBDyG86ySXfwB8JoTwGaB/pZ23AP/fMGt5GFg2hHbtw3zdQ5lj9aDz/Z83nKJdzaB2GiOd3b3cs3on4ORaSZKkkRrqEJ2bh9DmkwM+jsMtJMbYDmwY7vOGYG3muGDQ+YWDrq8bYjuNkQef3ENrRw9L5k5mydxxmactSZKUc4Ya8J8FXj3EtgH4ycjKGROrgd3A1Rz7i8rLgEbgPoAY474QwiOZdgxq1wPcPrZl5rcYIz95KD2f+TWXLzxFa0mSJJ3IUAP+F2OMO4b6oiGET4+wnlEXY0yFED4G3BxC+NcY45oQwrnAe0ivo985oPnHgXtDCK+NMf44hDAH+DDwDzHGPeNfff7YtKuJzbuPUF1RzBXnz066HEmSpKw1pIA/cLfXIbb/n5GVMzwhhCuBfwGqMqf+MoTwIeAzMcbvDajnlhBCIfDVEAKkJxd/JMZ406C6Hwwh3MALcwpKgC8Bnx3zbybP9ffeX3vxfEqLCxOuRpIkKXuNdJnMCSHG+ABw/hDbfgv41hDa/Qz42elVpuFobe/mgSf3EAJcf+mCpMuRJEnKaiPd6EoaNfc/vpue3hTnLZ7KzPrKpMuRJEnKagZ8Je7nv04vjXndS9zYSpIk6XQNOeCHEC4PmQHs0mjZuucIW3YfobK8mEvOmZl0OZIkSVlvOD34PwcOhhC+GUK4MYQweYxqUh65J9N7f/UFcyhxcq0kSdJpG07AnwK8D2gDPg8cCCHcH0L4WAhhKDvQSsfo6e3jvt/sAuDaix2eI0mSNBqGHPBjjB0xxh/HGD8YY5wLvIR0r/4bgHUhhK0hhC+GEF4RQigZq4KVOx5bv5+W9h4Wzqph0exJSZcjSZKUE0Y8yTbG+ESM8TMxxkuAmcBfZo7/AzSEEH4YQnjfKNWpHHTP6szk2ovn4/QOSZKk0TEqq+jEGA/GGG+OMb4JqAdeB2wnvTOs9CKHmjp4YuMBigoLuOqCOUmXI0mSlDNGfZnMGGNvjPEXMcYPxxjPGu3XV264d80uUhEuOWcGNZWO6JIkSRotroOvcRdjPGZ4jiRJkkaPAV/jbv3WBvY1tFE/qYzzlk5NuhxJkqScYsDXuPt5pvf+movmUVjg5FpJkqTRZMDXuGrv7OGhtXuBdMCXJEnS6DLga1w98OReurr7OHdRPTPrK5MuR5IkKecUjeaLhRCmAZ8GzgPK+s/HGC8Yza+j7HXP6h2AO9dKkiSNldHuwf8v0uvf1wP/D9gL/GSUv4ay1O4DLWzYcZjy0iIuWzEz6XIkSZJy0mgH/Lkxxs8BnTHG24HfAS4b5a+hLHX/43sAuHzFLMpKRvXNI0mSJGWMdsDv6j+GEGqBXsBtSkWMkfuf2A3A1e5cK0mSNGZGuxv1uRBCHfBNYDXQCjwxyl9DWWjTrib2HWqjtrqUcxbXJ12OJElSzhpywA8hhBhjPFmbGOM7Mx/+cwhhDVAL3Hka9SlH3P94uvf+ypWzXftekiRpDA2nB785hHBZjHHdUBrHGB8aYU3KMX2pyK+eTI+/v2qlw3MkSZLG0nDG4FcC5f2fhBAKQgi/CSEsHNgohFAWQqgZrQKV/dZtPkhTSxez6itZMndy0uVIkiTltNOZZBuAlaSH4Qx0HtB4Gq+rHNO/es5VF8whBIfnSJIkjaWx2snWFCcAunv6eHjdXgBeunJ2wtVIkiTlvrEK+BIAa559nvbOXhbPmcScadVJlyNJkpTzhhvwT7qKjjRY/9r3V7n2vSRJ0rgY7jr4vwghrAeeAp4hHfiLR70q5YS2jh5+/czzhABXnu/wHEmSpPEwnID/AeD8zONtQFXm/EMhhK3AOmAt6d1rJR5Zt4+e3hQrFtczZVL5qZ8gSZKk0zbkgB9j/M+Bn4cQlpAO+yszx0uB1/c3H53ylM36N7d6qWvfS5IkjZvhDtE5Ksa4CdgEfLf/XAhhGnAB6aUylccamztZu/kgRYWBy1fMTLocSZKkvDHigH88McYDwF2Zh/LYg0/uIRXhorOmU1VRknQ5kiRJecNlMjUmXD1HkiQpGQZ8jbq9h1p5bmcT5aWFXLR8etLlSJIk5RUDvkbdQ0+ld659yTkzKSsZ1VFgkiRJOgUDvkbdw+v2AXD5ilkJVyJJkpR/DPgaVQca29m8q4mykkJWnjkt6XIkSZLyjgFfo+qRp9O99xcum05pcWHC1UiSJOUfA75G1cNr0+PvLz/X4TmSJElJMOBr1Bxu7uTZ7Y0UFxVw4TKH50iSJCXBgK9R8+jT+4gRVi6dRkVZcdLlSJIk5SUDvkbNw2vT4+8vPXdmwpVIkiTlLwO+RkVLezdrtxyioCDwknNmJF2OJElS3jLga1Q89vR+UqnIikX1VFeUJF2OJElS3jLga1Q8ktnc6rIVDs+RJElKkgFfp629s4cnnjtACHDJOQZ8SZKkJBnwddp+8+wBenpTLFtQR21NWdLlSJIk5TUDvk7bQ+vSm1tdtsLNrSRJkpJmwNdp6erp4zfPPg/ApQ7PkSRJSpwBX6fliY0H6OzuY/HcyUyrq0i6HEmSpLxnwNdpeXhtZniOm1tJkiRNCAZ8jVhPb4rVz6SH5zj+XpIkaWIw4GvE1m0+RFtHD/NmVDN7alXS5UiSJAkDvk7Do+vTm1td6vAcSZKkCcOArxGJMbIms3rOxctnJFyNJEmS+hnwNSI79rdw8HAHk6tLWTxnctLlSJIkKcOArxH59TP7AVh11nQKCkLC1UiSJKmfAV8j0j88Z9Xy6QlXIkmSpIEM+Bq25rZuNmxvpKgwsHLp1KTLkSRJ0gAGfA3b4xsPkIpw9hlTqCgrTrocSZIkDWDA17D1j7+/yNVzJEmSJhwDvoalry/F4xsOAHDRMsffS5IkTTQGfA3Lhh2Hae3oYfbUSma5e60kSdKEY8DXsBxdHnOZw3MkSZImIgO+huXXmeUxL3J5TEmSpAnJgK8he76xnZ37WygvLWL5wilJlyNJkqTjMOBryNZkhudccOY0iov80ZEkSZqITGkasv7hOatcPUeSJGnCMuBrSDq7elm7+RAhwIXLpiVdjiRJkk7AgK8hWbv5ED29KZbMnUxtdVnS5UiSJOkEDPgaktXuXitJkpQVsj7ghxAKQggfDyF0hRDek3Q9uSjGyBrH30uSJGWFoqQLOB0hhHnA14FJQMlJ2t0MXAG0Drp0S4zxs4PavgL4DFAOFAP/DXwuxpgavcqzy7a9zTQc6aSupoxFsyclXY4kSZJOIqsDPvBR4KvATuCXp2j7/hjjfSdrEEK4ArgDeHOM8bYQwlzgMaAG+MTpl5udfv1s/+610wkhJFyNJEmSTibbh+h8NMb49VF8vc8Bj8YYbwOIMe4CvgB8NIQwaxS/TlZ5YuNBAFa5eo4kSdKEl9UBP8bYO1qvFUKYCVzGi98JuJf0UJ3XjtbXyibtnT1s2N5IQYBzF09NuhxJkiSdQlYH/GF6Zwjh4RDCsyGER0MIHwkhDByidG7muG3Q8/o/XzH2JU48T29toC8VWTKvlqry4qTLkSRJ0ilk+xj8oWrJPD4KHAFeBnyXdI/9GzNt6ge0Hag5c5xyvBcOIXwA+ADAvHnzRq/iCeLJ59LDc85fau+9JElSNpgwPfghhGtDCHEIj/uG+9oxxj+MMX4yxtgU0+4FPgu8IYRw+enUHWO8Kca4Ksa4aurU3AvB/QF/5VLH30uSJGWDidSD/zCwbAjt2kfp6z2WOV4CPAQcynxePahdTebYMEpfN2s0HOlg1/MtlJUUsnRebdLlSJIkaQgmTMCPMbYDG0b7dUMIhcCkGGPjoEt9mWNh5rguc1wwqN3CzHHtaNc20fX33p+zqJ7iognzZo8kSZJOIh9S21zgkeOcvzBzfBwgxrgv0+7qQe1eBvQAt49RfRPWC8Nzcm/okSRJUq7Kh4APsDQzGRaAEMIK0htX/Qr4xYB2HwcuCyG8NtNuDvBh4B9ijHvGsd7ExRh5cpMTbCVJkrJNVgf8EMKVIYQngf/MnPrLEMKTIYQ3Dmi2F/gg8JYQwtoQwkbgB8DXgVfHGGN/wxjjg8ANwP8NIawF7gG+BPzZ2H83E8uO/S00tXRRV1PG3OmDpyVIkiRpopowY/BHIsb4AHD+Kdp0AzdlHkN5zZ8BPzvt4rLck88dANK99yGEhKuRJEnSUGV1D77GzhOufy9JkpSVDPh6kZ7ePtZvTa8Kev4SA74kSVI2MeDrRTZsP0xXdx/zZ1RTW1OWdDmSJEkaBgO+XuSJo+Pv3b1WkiQp2xjw9SJPOv5ekiQpaxnwdYzW9m42726iqDBwzhlTki5HkiRJw2TA1zGe2nyIGOGsBXWUlWb1KqqSJEl5yYCvYzg8R5IkKbsZ8HWM/g2uVjrBVpIkKSsZ8HXU/oY29je0U1lezKI5k5MuR5IkSSNgwNdR/cNzViyup7AgJFyNJEmSRsKAr6Oe3JQO+Csdfy9JkpS1DPgCIMbI01sOAXDeEgO+JElStjLgC4DdB1o50tpNXU0pM+srky5HkiRJI2TAF8DR3vuzz6gnBMffS5IkZSsDvgB4emsDAOcscvdaSZKkbGbAV2b8fTrgn32GAV+SJCmbGfDF/oZ2Gps7qaksYd706qTLkSRJ0mkw4GvA+Pspjr+XJEnKcgZ8vTD+3uE5kiRJWc+Ar6MB3/H3kiRJ2c+An+cOHG7nQGM7lWVFLJg1KelyJEmSdJoM+HnumUzv/bKFUygscPy9JElStjPg5znH30uSJOUWA36e61//3g2uJEmScoMBP48dbu5kz8FWykoKWTRnctLlSJIkaRQY8PPY+m3p3vuzFtRRVOiPgiRJUi4w1eWxo8NzHH8vSZKUMwz4eWx9/wTbRfUJVyJJkqTRYsDPU81t3Wzf10xxUQFL5k5OuhxJkiSNEgN+nnomM/7+zPm1lBQXJlyNJEmSRosBP0/1j78/2/H3kiRJOcWAn6fWbz0EOMFWkiQp1xjw81BbRw9b9xyhsCBw1vy6pMuRJEnSKDLg56FntzeSirB47mTKSouSLkeSJEmjyICfh57e4vAcSZKkXGXAz0Oufy9JkpS7DPh5prunj827mwgBli1w/L0kSVKuMeDnmS27j9DbF5k3vZrK8uKky5EkSdIoM+DnmY07GwE409VzJEmScpIBP89s2H4YgLPm1yZciSRJksaCAT/PbNzR34NvwJckScpFBvw8cqipg0NHOqksK2LOtOqky5EkSdIYMODnkY070sNzls6rpaAgJFyNJEmSxoIBP49syAzPOcvlMSVJknKWAT+P9PfgO/5ekiQpdxnw80RPb4rNu5sAOHOeAV+SJClXGfDzxLa9R+jpTTFnWhVVFSVJlyNJkqQxYsDPExu2Z8bfu8GVJElSTjPg5wnH30uSJOUHA36ecAUdSZKk/GDAzwONzZ0cONxBeWkRc6e7wZUkSVIuM+DngRc2uJpMoRtcSZIk5TQDfh7YmBmec6YTbCVJknKeAT8PbMj04J/lBFtJkqScZ8DPcb19KTbtagJgqRtcSZIk5TwDfo7bvq+Z7p4+ZtVXMqmqNOlyJEmSNMYM+Dlu4/b+8ff23kuSJOUDA36O27AzM/7e9e8lSZLyggE/x23cntnB1vH3kiRJecGAn8OOtHaxr6GN0pJCFsysSbocSZIkjQMDfg7r3+BqydzJFBZ6qyVJkvKBqS+Hbejf4MrhOZIkSXnDgJ/D+nvwnWArSZKUPwz4OaovFXkus4KOS2RKkiTlDwN+jtp9oIXO7j6m1ZZTW12WdDmSJEkaJwb8HLVldxMAi+ZMTrQOSZIkjS8Dfo7avPsIAIsN+JIkSXnFgJ+jNu9qAgz4kiRJ+caAn4P6UpGte9M9+IvmTEq4GkmSJI2nrA34IYTaEMKfhBAeDiE8HkJ4NoRwdwjhihO0f3sI4akQwtoQwtMhhN87QbtXhBBWhxDWhRA2hBA+EULIqj+nPQda6OruY2ptOZOqSpMuR5IkSeOoKOkCTsMHgT8Bro4xPh1CKAT+Abg/hHB9jPHu/oYhhBuBrwJXxhhXhxBWAI+EEApijF8e0O4K4A7gzTHG20IIc4HHgBrgE+P3rZ0ex99LkiTlr6zqmT6O/4gxPg0QY+wD/hToA/6wv0Gm9/3zwK0xxtWZtmuBm4G/CSEM7OL+HPBojPG2TLtdwBeAj4YQZo35dzNKXlhBx+E5kiRJ+SabA/7ngf838ESMsRM4DAzc2eliYA7wy0HPvxeoA14GEEKYCVx2gnbFwGtHq/CxtmVPZvz97MnJFiJJkqRxl7UBP8bYl+m1PyqEUAvUA/cNOL0ic9w26CW2Dbp+7hDbTWipVGTrnibAITqSJEn5KGsD/gm8HzhAelhNv/rMsWVQ2+bMccow2x0jhPCBEMKaEMKagwcPDr/iUbbnYCsdXX3UTypjcrUTbCVJkvLNhAn4IYRrQwhxCI/7TvD8s0mPwX9rjLFhKF9yNOqOMd4UY1wVY1w1derU0XjJ0+IOtpIkSfltIq2i8zCwbAjt2gefCCHMA24Hfi/GeN+gy4cyx+pB5/s/bzhFu5pB7Sa0oyvozJ2cbCGSJElKxIQJ+DHGdmDDcJ+XWcrybuDjMcYfHKfJ2sxxwaDzCwddXzfEdhPa5kwPvuPvJUmS8tOEGaIzEiGEOcDPgT+PMX5vwPkfD2i2GtgNXD3o6S8DGslMyI0x7gMeOUG7HtLvEExo6Qm27mArSZKUz7I24GfC/X3A00BZCOEd/Q8GrHgTY0wBHwNuDCGsyjz3XOA9wCczS2v2+zhwWQjhtQO+xoeBf4gx7hn77+r07D3USkdXL3U1ZdRWlyVdjiRJkhIwYYbojMBHgEWZxxsGXdsx8JMY4y2ZnW6/GkKA9C82H4kx3jSo3YMhhBuAz4QQPgOUAF8CPjs238Lo2uIOtpIkSXkvawN+jPEjpEP+UNt/C/jWENr9DPjZaZSWmBfG3zs8R5IkKV9l7RAdvVh/D/4iV9CRJEnKWwb8HJFKRba4g60kSVLeM+DniP0NbbR39lJXU0pdjRNsJUmS8pUBP0dsdgdbSZIkYcDPGf072C6aPTnZQiRJkpQoA36O2OIKOpIkScKAnxNijC8EfFfQkSRJymsG/Bywv6Gdts5eJlc7wVaSJCnfGfBzwAsbXE0ms1OvJEmS8pQBPwdsObqCjuPvJUmS8p0BPwcM7MGXJElSfjPgZ7n0BFuXyJQkSVKaAT/LPd/YTmtHD5OqSqif7ARbSZKkfGfAz3L7DrVRVFjAIifYSpIkCQgxxqRryBmrVq2Ka9asGfev29OborW9m1qXyJQkScoXJ+zZtQc/BxQXFRjuJUmSBBjwJUmSpJxiwJckSZJyiAFfkiRJyiEGfEmSJCmHGPAlSZKkHGLAlyRJknKIAV+SJEnKIQZ8SZIkKYcY8CVJkqQcYsCXJEmScogBX5IkScohBnxJkiQphxjwJUmSpBxiwJckSZJyiAFfkiRJyiEGfEmSJCmHGPAlSZKkHBJijEnXkDNCCAeBHQl9+XrgUEJfW2PLe5vbvL+5y3ubu7y3uS1b7u+hGOOrjnfBgJ8jQghrYoyrkq5Do897m9u8v7nLe5u7vLe5LRfur0N0JEmSpBxiwJckSZJyiAE/d9yUdAEaM97b3Ob9zV3e29zlvc1tWX9/HYMvSZIk5RB78CVJkqQcYsCXJpAQwgMhhBhCWJB0LZIkKTsZ8LNACGFRCOHvQwi/yTyeywTB1wxqVxVC+NcQwsYQwjMhhLtDCGcnVbeGJ4TwBuCKE1zz3mahEMKCEEJrCOHJ4zwmD2jn/c1SIYQ3hBB+lfm3eWsIYU0I4Z0Drntvs0wI4eYQwubj/J3dEULoCiGUZ9p5b7NUCGFVCOHOEMKzIYR1IYTVIYQ3DWqT1ffXgJ8drgduBN4SY7wQOAt4BPhxCOGqAe2+C6wEVsYYlwOPAfeFEGaPd8EanhBCCfC3wE9P0MR7m73WxBjPP86jaUAb728WCiF8GPgz4G2Zf5vPBJ4DrhnQzHubnd4/+O8scBdwW4yxI9PGe5uFMu+Q/4L0RlbnxhjPBb4K/E8I4YYBTbP7/sYYfUzwB/B60v/YDDw3GYjAP2Y+vy7z+csHtCkBGoEvJf09+DjlPf4w8G3gLzL3ccGAa97bLH0AC4D7TtHG+5uFj8y97QIuGnR+FrDKe5u9D+BSYPqgc5XAEeAa7212P4D/nbl3KwedPwJ8J1furz34WSDG+MMY438OOl2TOR7MHN8A9AAPDnheN/BQ5pomqBBCHfAx4JMnaOK9zW3e3+z0TqApxvjrgSdjjHtjjGsyn3pvs1CM8ZEY4/ODTr8ZOADcm/nce5u9ejPHov4TIYRAelRLYeZU1t9fA34Wyrw99CXg8cwRYAWwN/MDONA2YHoIYdo4lqjh+RTwzRjj9hNc995mt+khhG+GEJ7IzJ/5dgjh3AHXvb/Z6TJge2YM/gMhhA0hhIdDCO8b0MZ7mzt+D/hKzHTl4r3NZrcAG4A/z4yzLyDdwVYK/EemTdbfXwN+FslMtt0M7Cb9W+brYozNmcv1QMtxntZ/fco4lKhhCiEsJt0z9Ncnaea9zV59pHuL/gW4EFhFulfosRDCRZk23t/sNBc4G/gT4E3AcuALwE0hhD/LtPHe5oAQwjLSf39vHnDae5ulMrnpGqCM9Dj8A8D7gOtijP3v0GT9/TXgZ5EY45YY42JgEumJXE+FEI676soAYewr02n4O+CzMcYjI3iu93aCizHuijGeG2N8LMaYyvzH8r+ANuBvTvF07+/EVkZ6XPbHYoz7M/f3u8CPgE+GECpO8lzvbXb5PeBHMcYDQ2jrvZ3gQghnAquBHUAdMI30ZPkfhBCuP9XTx7i8UWPAz0KZkPBh4Hng3zKnDwHVx2nef65hHErTMIQQrgTOAf79FE29tzkkplfgWAdckjnl/c1O/b17Tw46/wRQQbpH33ub5TIrnL0TuGnQJe9t9vor0guV/HGMsT3zy/ktwK+A/w4hFJED99eAnwVCCOWZCSBHZcYBrgPOCSGUAmuBWZl/jAZaCDw/xJ4Hja/rSA+1+nX/Osuke3cBfpo592q8t1krhDDpOPcN0kN3+idzeX+z04bMcfD/o30Dzntvs9/rSK+u8otB57232etcYHd8YbnTfs8BU0nfw6y/vwb87HAnL/T2DbSA9HiwbuAHQDHpiV/A0Z6Hy4Dvj32JGq4Y46dijIvisess90/weXXm3E/x3mazf2bQiguZe3cu6Uny4P3NVrdnjisGnT8H6ADW473NBe/n2Mm1/by32esAMDPTUz/QfNJLYx4mB+6vAT97fDqEMAXSyzmFEP4QuAj4Yky7G/gZ8FcDxn7+GZDi1GN9NYF5b7Pex0IIMwFCCIXA50n3En0avL9Z7Fbg18BnQghVcHTY3RuBv44xtnlvs1sIYT7wUuBrg695b7Pav5Beavwv+0dHhBBeBvwOcGuM8VAu3N/w4l9KNdGEEC4n3YtwEekVOcpIj//6N+Db/T0LIYRq4LOkh370kV5t50MxxvVJ1K2hywzF+RtgBjAdeBbozvTqe2+zVGY5zA8CV2ZO1ZO+t38dY/zlgHbe3yyU2cPic8ArgE7SG1/9S4zxKwPaeG+zVAjhL4FlMcY3neC69zZLhRBeCfwp6f9z+0gH92+Q7jTtyrTJ6vtrwJckSZJyiEN0JEmSpBxiwJckSZJyiAFfkiRJyiEGfEmSJCmHGPAlSZKkHGLAlyRJknKIAV+SJEnKIQZ8SZIkKYcY8CVJkqQcYsCXJI2qEML7QghPhxDaQwiPhBAWDeE5FSGEjhDClado964Qwq9DCKtDCG8fvaolKXcY8CVJoyaE8FHgS8DXgNcBZcC3h/DUa4EO4OGTvPalwN8CN2Re+/MhhJecXsWSlHuKki5AkpQbQgglwJ8Dn4kx/kPmXCXwgxDC/BjjjpM8/beAu2KMfSdpcxnwYIxxf+a1HwQuBx4blW9AknKEPfiSpNEyF5gMDAzyjZnjjFM899XAHadoswa4PIQwPYQwA8O9JB2XPfiSpNMWQigEngdagY+HEDYCG4EzM01aTvLclaR/AbjrZF8jxnh/COGvgJ8CKeCTMcaHRqF8ScopIcaYdA2SpCwWQqgGjgAfJx3q/xWYN6BJL1AVY+w6wfP/L/CKGONJJ9hKkobGITqSpNO1HAjAYzHG22OM84GZwHnAZuA3Jwr3Gb/FqYfnSJKGyIAvSTpd/ePrj469z0yE3QmcAfzkRE8MIUwDVmHAl6RRY8CXJJ2u/v9LwqDz7yI9Vv7rJ3nuq4FdMcb1Y1GYJOUjA74k6XRtzRwv7D8RQpgK/F/gy0NYHtPee0kaRU6ylSSdthDCU6Q3tfo/QDnwWdITb6+JMbaf4DnFwCHgLTHGk66gI0kaOnvwJUmj4XeAvcCPgX/LHK89UbjPeCnp5ZrvG/PqJCmP2IMvSUpECOELwBkxxt9OuhZJyiX24EuSkvIaHH8vSaPOHnxJkiQph9iDL0mSJOUQA74kSZKUQwz4kiRJUg4x4EuSJEk5xIAvSZIk5RADviRJkpRDDPiSJElSDjHgS5IkSTnk/wcBBXNd2gTKiwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "thetas = np.linspace(np.pi/6,np.pi/2*0.9,100)\n",
    "plt.plot(thetas*180/np.pi,Ea(thetas)/1e3)\n",
    "plt.xlabel('$\\\\vartheta$ / °')\n",
    "plt.ylabel('$E_\\\\mathrm{a}$ / kN m$^{-1}$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The active earth pressure is found at the peak of this curve. In other words we need to identify\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial E_\\text{a}}{\\partial \\vartheta} = 0\n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "    \\vartheta_\\text{a} = \\underset{\\vartheta \\in [0,\\pi/2]}{\\text{arg max }}E_\\text{a}\n",
    "$$\n",
    "\n",
    "As we saw above, this is cumbersome analytically. However, we can find the solution numerically in an iterative manner.\n",
    "\n",
    "As we have\n",
    "\n",
    "$$\n",
    "    \\left. \\frac{\\partial E_\\text{a}}{\\partial \\vartheta} \\right|_{\\vartheta_\\text{a}} = 0\n",
    "$$\n",
    "\n",
    "we find for other angles the following residual:\n",
    "\n",
    "$$\n",
    "    r(\\vartheta) = \\frac{\\partial E_\\text{a}}{\\partial \\vartheta} \\neq 0\n",
    "$$\n",
    "\n",
    "to which we seek the solution\n",
    "\n",
    "$$\n",
    "    r(\\vartheta_\\text{a}) = 0\n",
    "$$\n",
    "\n",
    "Since we don't want to go through the trouble to determine the first derivative of $E_\\text{a}$ with respect to $\\vartheta$ analytically, we do so numerically by numerical differentiation using a central finite difference scheme:\n",
    "\n",
    "$$\n",
    "    \\frac{\\partial E_\\text{a}}{\\partial \\vartheta} \\approx \\frac{E_\\text{a}(\\vartheta + \\epsilon) - E_\\text{a}(\\vartheta - \\epsilon)}{2\\epsilon}\n",
    "    \\quad\n",
    "    \\text{with}\n",
    "    \\quad \\epsilon \\ll 1,\\ \\epsilon > 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dEa_dtheta(theta,eps=1.e-4):\n",
    "    return (Ea(theta+eps) - Ea(theta-eps))/(2.*eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The former function defines the slope and equals our residual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residual(theta):\n",
    "    return dEa_dtheta(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we use Newton's method to find the solution to the non-linear equation $r(\\vartheta) = 0$. We depart from a Taylor series expansion to first order around an initial guess $\\vartheta_0$:\n",
    "\n",
    "$$\n",
    "    r(\\vartheta) \\approx r(\\vartheta_0) + \\left. \\frac{\\partial r}{\\partial \\vartheta} \\right|_{\\vartheta_0} \\Delta \\vartheta \\overset{!}{=} 0\n",
    "    \\quad\n",
    "    \\text{with}\n",
    "    \\quad\n",
    "    \\Delta \\vartheta = \\vartheta - \\vartheta_0\n",
    "$$\n",
    "\n",
    "We obtain the next guess for $\\vartheta$ by rearranging:\n",
    "\n",
    "$$\n",
    "    \\vartheta = \\vartheta_0 - \\left( \\left. \\frac{\\partial r}{\\partial \\vartheta} \\right|_{\\vartheta_0} \\right)^{-1} r(\\vartheta_0)\n",
    "$$\n",
    "\n",
    "As we linearized the non-linear function, the new guess is likely still not the solution we seek. Hence, we iterate:\n",
    "\n",
    "$$\n",
    "    \\vartheta_{i+1} = \\vartheta_i - \\left( \\left. \\frac{\\partial r}{\\partial \\vartheta} \\right|_{\\vartheta_i} \\right)^{-1} r(\\vartheta_i)\n",
    "$$\n",
    "\n",
    "We terminate the solution once $\\vartheta$ no longer changes, e.g. once\n",
    "\n",
    "$$\n",
    "    | \\vartheta_{i+1} - \\vartheta_i | < \\epsilon_\\text{tol}\n",
    "$$\n",
    "\n",
    "Another criterion, which is more direct, is to determine when the residual is close enough to zero:\n",
    "\n",
    "$$\n",
    "    | r(\\vartheta_{i+1}) | < \\epsilon_\\text{tol}\n",
    "$$\n",
    "\n",
    "In order to code this algorithm, we need the first derivative of the residual with respect to $\\vartheta$ (this is actually the second derivative of $E_a$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dr_dtheta(theta,eps=1.e-4):\n",
    "    return (residual(theta+eps) - residual(theta-eps))/(2.*eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now code the algorithm and give it a starting guess $\\vartheta_0$ and a termination tolerance $\\epsilon_\\text{tol}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_theta(theta_0=np.pi/6, tol=1.e-6):\n",
    "    theta = theta_0\n",
    "    iteration_counter = 0\n",
    "    res_vec = np.array([])\n",
    "    theta_vec = np.array([theta])\n",
    "    while (np.abs(residual(theta)) > tol and iteration_counter < 100):\n",
    "        res_vec = np.append(res_vec,residual(theta))\n",
    "        theta -= residual(theta) / dr_dtheta(theta)\n",
    "        theta_vec = np.append(theta_vec,theta)\n",
    "        iteration_counter += 1\n",
    "    return theta, iteration_counter, res_vec, theta_vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running this algorithm will give us the solution we seek."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_a, it, res, sols = find_theta()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slip angle of 62.43 ° was found after 8 iterations. The corresponding earth pressure is 0.02 MN/m.\n"
     ]
    }
   ],
   "source": [
    "print(\"The slip angle of %.2f ° was found after %i iterations. The corresponding earth pressure is %.2f MN/m.\" \n",
    "      %(np.rad2deg(theta_a),it,Ea(theta_a)/1.e6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result was found after several iterations. We may lower this number by improving our starting guess. For example, we could use the classical solution\n",
    "\n",
    "$$\n",
    "    \\vartheta_0 = \\frac{\\pi}{4} + \\frac{\\varphi}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slip angle of 62.43 ° was found after 4 iterations. The corresponding earth pressure is 0.02 MN/m.\n"
     ]
    }
   ],
   "source": [
    "theta_a, it, res, sols = find_theta(np.pi/4 + phi/2)\n",
    "print(\"The slip angle of %.2f ° was found after %i iterations. The corresponding earth pressure is %.2f MN/m.\" \n",
    "      %(np.rad2deg(theta_a),it,Ea(theta_a)/1.e6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We arrive at the same result, albeit in fewer iterations. We can check this result graphically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9ebc11089b8b4307a878e2e01b53265c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedFloatText(value=45.0, description='$\\\\vartheta_0 / °$', max=90.0, min=10.0), Outp…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_iterative_interactive()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the Newton algorithm reaches quadratic convergence once we are in the proximity of the solution.\n",
    "\n",
    "If we choose an initial guess that's outside of the convergence region for the Newton scheme, we may not converge at all or to the wrong solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare this result to the one from DIN4085, for example:\n",
    "\n",
    "$$\n",
    "    \\vartheta_\\text{ag} = \\varphi + \\text{arctan} \\left[ \\frac{\\cos(\\varphi - \\alpha)}{\\sin(\\varphi - \\alpha) + \\sqrt{ \\displaystyle \\frac{\\sin(\\varphi + \\delta_\\text{a}) \\cos (\\alpha - \\beta)}{\\sin(\\varphi - \\beta) \\cos (\\alpha + \\delta_\\text{a})} } }  \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def theta_ag_norm(phi,alpha,delta,beta):\n",
    "    t1 = np.cos(phi - alpha)\n",
    "    t2 = np.sin(phi - alpha) + np.sqrt((np.sin(phi + delta)*np.cos(alpha-beta))/(np.sin(phi-beta)*np.cos(alpha+delta)))\n",
    "    return (phi + np.arctan(t1/t2))*180/np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For the non-cohesive case the norm gives us 52.12 °.\n"
     ]
    }
   ],
   "source": [
    "print(\"For the non-cohesive case the norm gives us %.2f °.\" %theta_ag_norm(phi,alpha,delta,beta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tasks:\n",
    "\n",
    "* Check whether you get the same result by setting $c=0$.\n",
    "* Determine why some starting guesses don't work.\n",
    "* Play with the parameters and learn something."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bb4ad657f9b645d185aa2f793334b478",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='$c$ / kN m$^{-2}$:', max=50, step=5), IntSlider(value=3…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_earth_pressure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Erddruck.png](https://raw.githubusercontent.com/nagelt/Teaching_Scripts/master/Images/Erddruck.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}