04c_Erddruck_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 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)⋅tan(α) + cos(vartheta))⋅(tan(α)⋅tan(β) + 1)⋅sin(varph\n",
       "-⎜C⋅(sin(vartheta) - cos(vartheta)⋅tan(β))⋅cos(varphi) + ─────────────────────────────────────────────────────────────────────────\n",
       " ⎝                                                                                                 2                              \n",
       "──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────\n",
       "                                      (sin(vartheta) - cos(vartheta)⋅tan(β))⋅cos(α + δₐ + varphi - vartheta)                      \n",
       "\n",
       "             ⎞ \n",
       "i - vartheta)⎟ \n",
       "─────────────⎟ \n",
       "             ⎠ \n",
       "───────────────\n",
       "               "
      ]
     },
     "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": "iVBORw0KGgoAAAANSUhEUgAAAt0AAAIXCAYAAABAROIgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABaCUlEQVR4nO3deXxU5d3///dk3zeyr+z7joCoBRSXoqJVqUpdEG5trXut1tr259aK9q7Wr1a9bd2rFa1CraKgSBGlLkBAQCBsISRkIQlZJnsmM+f3x5CRmAAhzOTMZF7Px2OamXOuM+czPZG8c+U612UxDMMQAAAAAI8JMLsAAAAAoK8jdAMAAAAeRugGAAAAPIzQDQAAAHgYoRsAAADwMEI3AAAA4GGEbgAAAMDDCN0AAACAhxG6vZhhGLJarWL9IgAAAN9G6PZidXV1io2NVV1dndmlAAAA4CQQugEAAAAPI3QDAAAAHkboBgAAADyM0A0AAAB4GKEbAAAA8DBCNwAAAOBhhG4AAADAwwjdAAAAgIcRugEAAAAPI3QDAAAAHkboBgAAADyM0A0AAAB4GKEbAAAA8DBCNwAAAOBhhG4AAADAwwjdAAAAgIcRugEAAAAPI3QDAAAAHhZkdgEAAPgawzDUZneosblNza12Nbe0qbm1TS02u1ptDrW02tVicz5sbXa1tTlksztka3O4njvshuwOQw7DkN1uyO5wyOFwvreh7746/0eyWCSLxSKLDn+1SAEBFgUGWBQYGOD8evh5UKBFwUGBrq/BgRYFBwcqJDhQocEBCnE9D1RoSKDCQ4IUGhKosNAghQQFyGKxmPd/LtBHEboBAH7JMAw1t9pVW98ia0Or6hpbVdfQqrpGm/N5Y6sammxqbG5TfZNNjc021+umljbZHYbZH8EjAixSaEiQIsIOP0KDFR4WpPDQIEWGBSsy3PmI+t7XmMgQRUUEKyo8RMFB/CEd+D5CNwCgT7G12XWotlnV1hZV1TWr2tqsKqvzdXVds2rrW1Tb0Krauha1tjnMLtfrOAypqcX5i8Wh2p69R3hooKIjQhQTGaKYqFDFRIYoNjJUsVEhiokMVXx0qOKiQxUb5fwaGhzo3g8BeCFCNwDAZxiGoZr6FpVXNaq8qkkHqxtVWdOkypomHaptUmVNs2rqWzxy7qDAAEWFf9fr2/4ICwlUeKhzeIZzuEaQQoMDFHp4CEdwcKCCgwIUHBigoKAA1/Ouhoa0DxtpH0Yii3T4mQwZMozDw04Of7U7DDkczq92h8P51W50MZzFOeyl9fCQl/bnrTa7c3hMa1uHr02He/ObWmxqarGf8P9XTS12NbU0qby6qVvtw0ODFB8dqviYMCXEhCk+JlQJ0WFKiHW+TowLV7+YMIWFElvgu/juBQB4FVubXWWHGlV2qEGllQ0qPfz1YFWjyqsaT7p3OsAixUSGKiYqRHGHe2FjIkMUHRmimIgQRR3uoXUOlQh2DakI8dPeWIfDUHNrmxqa2tTY4hxic+Sj/vDDOTzH+bW+qdU5ZKehVd0ZhdPes15S2XDMdpFhQUqIDVdirDOIJ8WFKyk+XIlxzkdSfAS95vBaFsMw+uagtD7AarUqNjZWtbW1iomJMbscAHAbwzBUXdeiA+V1Ki6v14GKeh0or1dxeb3KqxvVk59MARYpISZM/eLClRgbroTYMMVHhx7uOT3cgxodquiIEAUEcKNgb3A4DNU32WRtaFFtfavra019i2rqDj/qW1RT16wqa4uaWtpO+pxxUaFKTnAG8JT4CCUnRCjliIe//vIE8xG6vRihG0BfUFvfooJSqwrL6lR4sE6FZc7n9U22E3qf0JBAJcdHKDk+3Bmk4iOUHB/h6umMjw5VYCA38Pmy5pY2Vde1OMfg1zWrqrZZhw4/KmubDr9uOqm/diTEhCklIUKp/SKUlhiltH4RSkuMVFpilKIjgpm5BR5D6PZihG4AvsTuMFRcXqd9JVbtK6l1fa2u6/4Y6/DQIGUkRyk9MVJp/SKVlhip1H6RSk+MVFx0KIEIMgxD1oZW11j+ypomVbQ/qptUXt2oKmtzj/5aEhkWpLTESKUnRSkjyfl9mJ4UpfSkKEWFB7v/w8CvELq9GKEbgLdqD9h7DtRoz4Fa7SmqUX5JrVpau3fTXWJsmLJTY5SZEqXMpChlJkcrIzlK8QRruIGtzaHKGmcAP1h1+HGoUWVVDTp4qLFHN9vGRYcqM/nw92pS1OHnUUqOj2C4ErqF0O3FCN0AvEWVtVk791dr5/4q7Sqs0e6iajV3I2BHRwRrQHqs+qfFKDs1Rjmp0cpKiVYkvYYwUXNLm8qqGp036lY2uG7aLTnUoIoTvKcgJDhQmclRyk5xfm9npUQrOzVaqf0iFUgYxxEI3V6M0A3ADHaHocIyq7blH9L2fVXK21+lim5M/ZaSEKHBmXEakBGjAemxGpgeq36xYfRcw6e0z55TXFGvkop6lVQ2OG/yrahXzQkMlQoJClDm4QCec/gXzpy0GCXFhfPfhJ8idHsxQjeA3mBrc2hXYbW25R/Stn2HlFdQpcbmY88ikRQfriFZcRqc6XwMyoxTTGRIL1UMmKO+yabi8jodKHfOtlN0sE4HyutUWtnQrakRJee48f6H//rTPy1G/dNjlJMao3DmIO/zCN1ejNANwBPsdod2H6jR1j2V2rKnUjsKqo45Fjs0JFBDsuI0LDtew3ISNCwnXgkxYb1YMeDdWm12lVQ2qKisTvsPWlV0sE77S+tUWlnfrTBusUjpiZHOvxBlxLq+8t9Z30Lo9mKEbgDuYBiGSisbtGlnuTbtqtCWPZXHnA85LipUIwcmaOSAfho1oJ8GpMcwFR/QA602uw6U16uwzKqCUqv2l9WpoKRWlbXN3To+PjpUgzLjNCgjVoMyYzUoI05J8QxP8VWEbi9G6AbQU43NNm3eXaGNOyu0cWe5yqsaj9q2X2yYxgxO1JhBiRo9sJ/SEiP5oQ54UH1jqwpKnUE8v7hW+0qt2l9qla0b849HR4RoSFacBmXGHh7iFa/EOO6d8AWEbi9G6AZwIkoq6rV+x0Gt316mbfmH1Gbv+p/32KgQjR2cpLGDEzV2cCIhG/ACdrtDByrqta+4VnuLa7WvpFZ7D9R2axGpuKhQDcmO09DseA3NiteQ7DhFR3CPhbchdHsxQjeAY7E7DOUVVOmrb0u1bluZSiobumwXFBigkQMSNGFYsiYMTdKA9FjmFQZ8gGEYKq9uUn5xjfYeqD08L36Nautbj3tsWmKkhmbFa1iO8zEgPVbBQQwTMxOh24sRugF8X6vNri17KvXlVmfQPtoiH8kJEZo8IkWnjEjR6IH9FMbMCECfYBiGKmuatedAtXYX1WhPUY12F9Uct0c8OChAgzPjNCwnXsNzEjS8f7z6xYb3UtWQCN1ejdANQHIG7Y07y7X2mxKt216qppbOM40EBFg0on+Cpox0Bu2slGiGjAB+wjAMlR5qcC5cVVitXYXVyi+uVetxxognx4dreP8EjeifoOH9EzQgjZumPYnQ7cUI3YD/srXZtWlXhdZ+U6yvt5V1OW92SHCgJg1P1qmj0zR5ZApjOAG4tNkdKii1amdBlfIKq7Vzf7VKjzIErV14aKCG5ThnLho5IEHDsuP5K5kbEbq9GKEb8C8Oh6EdBVVanVuktZtL1NDFn4ujwoM1ZVSqpo1J0/ihSQoL4QcigO6prW/RzsJq5RVUKa+gWjsLq9VqO/oc/QEBFg3KiNWogf00emA/jRrYT1H8ct9jhG4vRugG/EPRwTqtzi3Smo0HVN7FcuuRYUGaOjpNPxifoXFDkrgZCoBbtNkdKiixakdBlbbvO6Tt+6pUZT36HOIWi5STGqPRA/tp9KBEjR7UT7FRob1YsW8jdHsxQjfQd9U32fT5pgNaua5Qu4tqOu0PCwnUqWOcQXvC0CQFBwX2fpEA/Er7bCnb8g8dDuGHVHSw/pjHZKdGa8ygRI0Z7JznnxB+dIRuL0boBvoWh8PQt/mVWvl1ob7YUtLpJqeAAIsmDE3SzElZOnVUKmMpAZiutr5F2/dVaVv+IX2bX6l9xbXHXNq+f1qMxg5O1LghSRo1sJ8iw4N7r1gvR+j2YoRuoG+oqWvRynX79fHX+1V2qPPKkAMzYjXrlCz9YEKG4qPDTKgQALqnocmmHQVV+nZvpbburdSeopqjhvCAAIuGZMZp7BBnCB85IMGv/2pH6PZihG7AdxmGoW35h7T8iwJ9sbWk0+qQUeHBmjkpU+dMydHAjFiTqgSAk9PYbNP2fVXauqdSW/ZWKv/A0UN4SHCgRg1I0PihSRo/NFn902L8aqEuQrcXI3QDvqex2abVG4r04ZcFKiyr67DPYpHGDUnSuVNyNHV0qkKC/bfHB0DfVN/Yqm/zD2nz7gpt2VPZ6d/BI8VGhWj8kGRNHO4M4QkxffsvfYRuL0boBnzHwapGLVubr4+/3t9pTu2YyBCdMyVbP5zWX6n9Ik2qEAB6X7W1WZv3VGrzrgp9s6tclbVHnx2lf1qMJg5L1sRhyRo5sO8NRSF0ezFCN+DdDMPQ9n1Veu/zvfpqa2mnP6mO6J+g80/rr9PHpfe5Hx4AcKIMw1BxRb0276rQpl3OnvCmls4Lf0nOGZzGDk7SpBHOEN4XOiwI3V6M0A14J7vD0JdbS7Rk9R7t+d50f8FBAZo5MVNzfjBQA9IZqw0AR9Nmd2jn/mpt2lmujTvLtedAjY6WSjOSonTKiBRNHpGikQP7+eR6BYRuL0boBrxLq82u/2wo0tJP93RaTjk+OlQXnD5AP5zWn3lqAaAHautbtHl3hXLzyrUxr1w19S1dtgsPDdL4oUmaPCJFp4xIUbyPjAUndHsxQjfgHRqabPrwi3167/N81dR1/CEwMCNWP5oxSGeMy/DJnhcA8EYOh6H8klrl5h1U7o5y7dxfddRZUYZkxWnKqFRNGZmqAekxsli8c0YUQrcXI3QD5qpvsun9z/bq35/nq6HJ1mHfuCGJuuzMIRo/NMlr/4EHgL6irrFVm3aWa/0OZwiva2ztsl1iXLgmj0zRlJGpGjs40atmiSJ0ezFCN2CO+sZW/fuzfL3/+V41HDETicUinTYmXZedNVhDsuJNrBAA/JfdYWh3YbXW7zio9dvLtK/E2mW7Wy8fr3On5vRydUdH6PZihG6gd9U3turdNXv1/tr8DtP+BQRYdNakLM2dNUQZSVEmVggA+L7y6kat335Q67aXacvuSrXZHbJYpFfvP8+rVvkldHsxQjfQO5pb2/T+5/lasnpPh2EkgQEWzZqcrR/PGtInpqsCgL6usdmmb3ZVaH9ZneadO8zscjogdHsxQjfgWbY2hz7+er/eWrlT1UfcIBkU2B62hyolIcLECgEAfUWQ2QUAQG9zOAx9tumAXl+Rp4NVja7tARbprFOydeW5wwjbAAC3InR7QF5enn73u99p7dq1cjgcysnJ0R133KGrrrrK7NIAv/ft3kq9+N632nOgtsP2aWPSdM3sEcpKiTapMgBAX0bodrNNmzZpxowZmjVrlvLy8hQbG6tXX31V1157rXbv3q0HHnjA7BIBv1RSWa9Xlm3Xl1tLO2wfPyRJ15w/QkOzmY0EAOA5jOl2I8MwNGHCBOXn56uoqEixsd8tAX3ttdfqH//4h3JzczV+/PhuvR9juoGTV9/YqjdX7tIH/81Xm/27f+76p8Vo4ZxRmjAs2cTqAAD+guXT3Ojzzz/X5s2bdf7553cI3JI0b948ORwOPfPMMyZVB/gXh8PQR18V6KePrNK/P9vrCtzx0aG69fLx+n93ziRwAwB6DcNL3GjVqlWSpEmTJnXa175t5cqVvVoT4I92FVbruaVbtLuoxrUtJDhQl8wcpMvOHKLwUP7pAwD0Ln7yuNGOHTskSRkZGZ32JSUlKTg4WPv371dTU5PCw8M7tWlpaVFLy3fTllmtXa+wBKBrtfUtem35Dn389X4dOXBu+vgMXXfhKCXFd/7vDgCA3kDodqOamhpJUmRk50U0LBaLIiIiVFtbq5qami5D9yOPPKIHH3zQ02UCfY7DYWjlukK9+sE21TV+t7hNdmq0brxkrMYMTjSxOgAACN1e5d5779Wdd97pem21WpWVlWViRYD3KzpYp2fe2axt+Ydc28JDg/ST84brwjMGKCiQW1cAAOYjdLtRXFycJKmhoaHTPsMw1NjY2KHd94WGhio0NNRT5QF9iq3NoSWrd+utlbvUZne4ts+cmKkFc0YpISbMxOoAAOiI0O1GI0aMkCQVFxd32ldRUSGbzaacnJwuh5YA6L7t+w7p6bc3q+hgnWtbar8I3Tx3nMYPZUYSAID3IXS70axZs/TQQw8pNze30772beecc05vlwX0Gc2tbXpt+Q69/3m+60bJgACLLpkxSFeeO0xhIfyTBgDwTvyEcqMf/OAHGjdunD788EPV1tZ2mKt78eLFCggI0M0332xihYDvyiuo0v97c6OKK74bvjU4K063/ni8BmbEHuNIAADMxx1GbmSxWPTyyy/LMAwtWLBAtbW1MgxDr7zyiv7xj3/od7/7XbdXowTg1Gqz65Vl23TP05+7AndwUIAWzhmlx26bTuAGAPgEerrdbMKECVq/fr1++9vfaujQoXI4HMrJydGrr76qq6++2uzyAJ+yu6haTyze1GHs9tDsON1x5URlpUSbWBkAACfGYhhHLiEBb2K1WhUbG6va2lrFxMSYXQ7Qa+wOQ0tX79Y/VuTJ7nD+ExUUGKCfnDdMl84crECmAQQA+Bh6ugF4lYrqJv15ca6+3fvdvNuDMmP1iysnKieNXz4BAL6J0A3Aa/x3S4me/uc3qm9yrippsUg/njVU884dxiI3AACfRugGYLrmljY9/+9v9fHX+13bEmPDdOdVkzRmEEu4AwB8H6EbgKkKy6x69O/rVXSw3rXt9HHpumXuOEVFhJhYGQAA7kPoBmCa1blFeuadzWpptUuSwkIC9bNLxmjW5GxZLBaTqwMAwH0I3QB6XavNrr+9u1UfffXdcJL+aTG659pTlJnMVIAAgL6H0A2gV5VWNujRv69XfnGta9vZk7P1s0vHsIw7AKDP4iccgF6zfnuZHv9Hrhqa2yRJIUEB+vllY3X2lByTKwMAwLMI3QA8zjAMvb1qt15fsUPty3GlJ0bq1/Mna0A6y7gDAPo+QjcAj2puadP/e2uT/ru5xLVt2pg03XHlBEWEBZtYGQAAvYfQDcBjyg416OGX16mg1OradvUPh+vys4cyOwkAwK8QugF4xJY9FXr01Q2qa2yVJIWHBumuqyZpyqhUkysDAKD3EboBuN3HX+/Xs+9slt3hHMCdnhip3y2cqqwUpgMEAPgnQjcAt3E4DL2+YofeXrXbtW3i8GTdffUpigpn/DYAwH8RugG4RavNrv/35iZ9/k2xa9tF0wdq4ZzRCgxg/DYAwL8RugGctNr6Fj388jrtKKiSJAVYpBt+NEYXnjHQ5MoAAPAOhG4AJ+VAeZ0efOErlR1qlCSFhQTq7mtO0ZSR3DAJAEA7QjeAHttVWK0Hnv/KNUNJQkyY7vufqRqUGWduYQAAeBlCN4Ae2bSzXIteWafmVrskqX9ajO6//lQlxoWbXBkAAN6H0A3ghH2+qVh/XpyrNrtzSsDRg/rpdwumKpIZSgAA6BKhG8AJ+eC/+/TXf22R4czbOnV0qu6++hSFBAeaWxgAAF6M0A2gWwzD0OKPd2rxxztd286Zkq2b545TYGCAiZUBAOD9CN0AjsswDL3w72/13uf5rm0/njVE18weIYuFObgBADgeQjeAY3I4DD33ry1a/kWBa9v/XDRaP5oxyLyiAADwMYRuAEflcBh6+u1vtHJdoSTJYpFuu3y8zp6SY3JlAAD4FkI3gC7ZHYaeemuT/rOhSJJzlck75k3UmZOyTK4MAADfQ+gG0Ind7tCf39ioz74pliQFBFh011WT9IPxGSZXBgCAbyJ0A+igze7QY6/n6r9bSiRJQYEW/eqaUzRtTLrJlQEA4LsI3QBc7A5Df35j4xGBO0D3XjdZU0ammlwZAAC+jdANQJLzpsmn3tqkzw8PKQkOCtDvFkzVxOHJJlcGAIDvY0ULADIMQ88t3eK6aTIo0KLfXDeFwA0AgJsQugE/ZxiGXnjvWy3/skCS86bJX11zik4ZkWJuYQAA9CGEbsDPvbZ8h977zLnSpMUi3TlvIjdNAgDgZoRuwI+99clOvb1qt+v1bZeP14yJmSZWBABA30ToBvzUii8L9PryPNfrGy8dy0qTAAB4CKEb8ENfbi3R/y3Z7Hq94MJRuuD0ASZWBABA30boBvzMt3sr9afXc+UwnK8vnTlYl5452NyiAADo4wjdgB8pKLXqDy99LVubQ5J01ilZmn/BSJOrAgCg7yN0A37iYFWj7v/bF2pobpMkTRqerFsvH6+AAIvJlQEA0PcRugE/UFvfovv/9oWqrC2SpKHZcfr1tZMVFMg/AQAA9AZ+4gJ9XKvNrodfXqfiigZJUkZSlO77n1MVFhpkcmUAAPgPQjfQhxmGoafe+kY7CqokSQkxoXrop9MUGxVqcmUAAPgXQjfQh725cpfWbDogSQoNCdT/9z+nKjkhwuSqAADwP4RuoI9as/GA3vjIufiNxSL98ieTNDgzztyiAADwU4RuoA/asa9KT761yfX6ugtGadqYNBMrAgDAvxG6gT6m7FCDHn7lu7m4z52ao0tmDjK5KgAA/BuhG+hDGptteujFr1Vb3ypJGjs4UT+/bKwsFubiBgDATIRuoI9wOAz9+Y2NKjpYJ8k5NeC985mLGwAAb8BPY6CP+OeqXfp6W5kkKSo8WPddP1VRESEmVwUAACRCN9AnrNte1mGmkruvPkXpiVEmVwUAANoRugEfV1JRrz//I1eG4Xx9zewRmjg82dyiAABAB4RuwIc1Ntv0h5fXqaG5TZJ02tg0zT1riMlVAQCA7yN0Az7KMAw9+dYm142TWSnRuv2KCcxUAgCAFyJ0Az7qnf/s1hdbSiVJEWFB+u2CKYoICza5KgAA0BVCN+CDtu6t1OvLd7he//KqScpI4sZJAAC8FaEb8DE1dS167PUNchy+cXLeucM0ZWSquUUBAIBjInQDPsThMPTE4o2qsrZIcq44ecU5w0yuCgAAHA+hG/AhS1bv1sad5ZKkuOhQ3XXVJAUGcOMkAADejtAN+Iht+Yf0+orvFsC56yeTFB8TZnJVAACgOwjdgA+wNrQ6x3EfHsh9+dlDNW5okslVAQCA7iJ0A16ufRx3ZW2zJGn0oH6ad+5wk6sCAAAngtANeLn3Pt+rDTsOSpJiIkMYxw0AgA8idANebF9JrV794Lv5uO/8yUT1iw03sSIAANAThG7AS7Xa7PrzGxvVZndIki6ZOViThqeYXBUAAOgJQjfgpV5bvkMFpVZJUv+0GF0zm3HcAAD4KkI34IW27KnQvz/bK0kKCgzQL6+apOCgQJOrAgAAPUXoBrxMfZNNTyzeJOPwMu/Xnj9C/dNizC0KAACcFEI34GX++q8tqqxpkiSNGZSoi6cPMrkiAABwsgjdgBf5/JtifZp7QJIUERakO+ZNUADTAwIA4PMI3YCXOFTbpGff2ex6feOlY5UcH2FiRQAAwF0I3YAXMAxDz7yzWfVNNknSGePSNXNipslVAQAAdyF0A17gs03FWr/duepkXHSofn7ZOFksDCsBAKCvIHQDJqutb9Hf3t3qev3zS8cqJjLExIoAAIC7EboBk/3tX1tlbWiVJJ0+Nl2njU03uSIAAOBuhO6jeOWVVxQeHq7U1NQuHwUFBV0e9+GHH2r69OlKTk5WSkqKZs+erXXr1vVu8fAZX39bqs++KZYkRUcE62eXjjG5IgAA4AmE7mO44oorVFZW1uWjf//+ndq/8MILuuCCC3TeeeeppKREBQUFGjBggM444wytXLmy9z8AvFp9k03PLvlutpLrLx6j+OgwEysCAACeQuh2k+LiYt12222aPn26fvvb3yooKEjh4eF66qmnlJWVpQULFqipqcnsMuFFXnrvW1VZWyRJk4Yn68xJzFYCAEBfReh2kxdffFFNTU2aN29eh+1BQUGaO3euiouLtXTpUpOqg7f5Zle5Vq4rlCSFhwbp5rnjma0EAIA+jNDtJqtWrZIkTZo0qdO+9m0MMYEkNbe26S9vfzesZMGFI5UUH25iRQAAwNMI3cewZ88ezZs3T0OGDFFSUpJGjx6tO+64Q8XFxZ3a7tixQ5KUkZHRaV/7try8vGOer6WlRVartcMDfc8/P9ml8qpGSdLoQf103qn9zS0IAAB4HKH7GPbs2aO5c+dq27ZtKiws1GOPPaYlS5Zo7Nix2rx5c4e2NTU1kqTIyMhO79O+rbq6+pjne+SRRxQbG+t6ZGVlueeDwGsUHazTvz7dI0kKCgzQzXPHKSCAYSUAAPR1FsMwDLOL8JRFixaptbW12+2vu+4616wkjY2Nstvtio6O7tBm1apVOvvsszVx4kTl5ua6toeEhMhms6mmpkaxsbEdjvnmm280YcIEDR06VDt37jzq+VtaWtTS0uJ6bbValZWVpdraWsXExHT7c8A7GYah3z33hbbsqZQkXXH2UF09e4TJVQEAgN4QZHYBnrRo0SI1NDR0u/3MmTNdoTsiIqLLNrNmzVJ8fLw2btyoffv2acCAAZKkuLg4VVRUqKGhoVPobq8hPj7+mOcPDQ1VaGhot+uFb/lsU7ErcCcnRGjurCEmVwQAAHpLnx5eUl9fL8Mwuv2YOXNmt943Pd25YmBpaalr24gRzh7LrsZ7t28bPnz4SX4i+KqGJptefO9b1+ufXTJGYSF9+ndeAABwhD4duk/GAw88oIqKii73lZSUSJKSk5Nd22bNmiVJHYactGvfds4557i7TPiIf3yUp+o659ChqaNSNWVkqskVAQCA3kToPooHH3xQq1ev7rR9zZo1qq6u1ogRIzR48GDX9oULFyo8PFyLFy/u0L6trU3vvPOOMjIydMkll3i8bnifvQdq9MHafElSSHCgfvojlnoHAMDfELqP4Ze//KU+/fRTORwO2e12ffbZZ1qwYIEiIyP14osvdmibmZmpJ598Up999pkWLVoku92u5uZm3XbbbSosLNRLL7101HHi6LscDkP/t3SLHIdvV77ynKFKTuD7AAAAf0PoPopPP/1UF198sW677Talp6crISFB8+fP16xZs7R582ZNmzat0zE33HCDPvjgAy1fvlypqanKzs5Wfn6+1q5dq3PPPdeETwGzrVxXqJ37nVNFZiZH6UczBh/nCAAA0Bf16SkDfZ3ValVsbCxTBvqo+iabfvbIJ7I2OKet/MONp2nckCSTqwIAAGagpxvwkLdW7nQF7h+MzyBwAwDgxwjdgAeUVNZrWfvNk0EBuu6CkSZXBAAAzEToBjzglWXb1WZ3jtz60czB3DwJAICfI3QDbrZlT4W+3OpcOCk+OlRzz2LlSQAA/B2hG3Aju8PQi//e5np97fkjFB7KypMAAPg7QjfgRqvWFyq/pFaSNDAjVmedkm1yRQAAwBsQugE3aWy26bXlO1yvb7h4tAICLCZWBAAAvAWhG3CTd/6zWzV1LZKk08amafSgRJMrAgAA3oLQDbjBwapGvbtmryQpKDBACy4cZXJFAADAmxC6ATd4fcUO2dockqSLpw9Uar9IkysCAADehNANnKSCUqvWbDwgSYqOCNGPZw01uSIAAOBtCN3ASXrtwx0ynOvg6PKzhygyPNjcggAAgNchdAMnYfu+Q1q3vUySlBgbpvNPG2ByRQAAwBsRuoEeMgxDf//wuykCrzx3uEKCA02sCAAAeCu3he7q6mr9/e9/d9fbAV4vN69c2/IPSZIykiJ19uQskysCAADeym2hu7CwUAsWLHDX2wFezeEw9NoRvdxXzx6hwED+cAQAALoW1N2GhYWFx9xfUlJy0sUAvmLt5mLXcu+DMmN12ph0kysCAADerNuhu3///rJYjr6ktWEYx9wP9BVtdodeX5Hnen3t+SNZ7h0AABxTt0N3fHy8fv/732vGjBld7s/Ly9Pll1/utsIAb7VyXaFKKxskSWMGJWrC0CSTKwIAAN6u26F70qRJqq6u1qhRXS9v3dbWJqN9smKgj2qx2fXmxztdr6+9YAR/4QEAAMfV7dD985//XA0NDUfdn52drZdfftktRQHe6qMvC1RlbZYkTR2VquE5CSZXBAAAfIHFoHvaa1mtVsXGxqq2tlYxMTFml+P3Wm123bDoE1fofuqXMzUgPdbkqgAAgC9gjjOgm1auK3QF7lNHpxK4AQBAtxG6gW6wtdn1zqpdrtdXnDPMxGoAAICvOanQHRgYqF27dh2/IeDjVq0vUmWts5d78sgUDc6MM7cgAADgU04qdDMcHP6gze7Q20f0cl9JLzcAADhBDC8BjmP1hiKVVzdJkiYOT9bQ7HiTKwIAAL6G0A0cg93u0D+P6OWeRy83AADoAUI3cAxrNh1Q2aFGSdL4IUka3p95uQEAwIkjdANHYXcYemvlEWO5z6WXGwAA9AyhGziKz78pVkmlcxXWMYMSNWpgP5MrAgAAvorQDXTB4TD0z0+O7OUeamI1AADA151U6L7nnnvUrx+9f+h7NuQdVNHBOknSiP4JGjMo0eSKAACALws6mYMfeeQRd9UBeJWlq/e4nl925mBZLBYTqwEAAL6u14aXDBkypLdOBZyUXYXV2pZ/SJKUmRylySNTTa4IAAD4upPq6f6+9evXq7CwUA0NDZ32VVZWuvNUgMcc2cv9oxmDFRBALzcAADg5bgndu3bt0o9+9CPt3Lmzy/2GYfDnefiE0soGfbm1RJIUFx2qMydlmlwRAADoC9wSum+++WY1NjbqN7/5jbKyshQaGtphv2EY+sUvfuGOUwEe9e6aPXIYzudzzhiokOBAcwsCAAB9gltC94YNG7R582ZlZ2cftc0TTzzhjlMBHlNb36JP1hdJksJCAnX+af3NLQgAAPQZbrmRMiUl5ZiBW5I2b97sjlMBHvPhf/ep1WaXJJ07NUdRESEmVwQAAPoKt4TuBQsW6OOPPz5mmzlz5rjjVIBHNLe2adl/90mSAgIsunj6IJMrAgAAfYlbhpfcc889euKJJ/Txxx9r6tSpSkpKUkBAxzz/5ZdfuuNUgEf8Z0ORrA2tkqQzxqUrOSHC5IoAAEBf4pbQXVJSojfffFMbNmxwx9sBvcruMPTup3tdry+dOdjEagAAQF/kltB9yy23KC8vTz/5yU+UnZ2tkJCOY2ENw9Cf//xnd5wKcLuvvi1V6SHn3PLjhiRqUGacuQUBAIA+xy2he/Xq1Vq7dq3GjBlz1DZ//etf3XEqwO3e/zzf9fzSmaycCgAA3M8tN1LGxMQcM3BL0v79+91xKsCt9pXUupZ8z0qJ0oRhSSZXBAAA+iK3hO7zzz9fubm5x2xzxx13uONUgFstW7vP9fyC0weycioAAPAIi2EYxsm+ybZt23T33Xfrhz/8oaZNm9Zp9hLDMDRlyhQdPHjwZE/lV6xWq2JjY1VbW6uYmBizy+lz6hpbdd1DH6vVZldEWJBeue88hYe6ZcQVAABAB25JGGPGjJHFYtFHH33kjrcDesXKrwtdi+HMmpxN4AYAAB7jlpQRFhamK6644qj7DcPQ22+/7Y5TAW5hdxj64Isjh5YMMLEaAADQ17kldMfGxurll18+Zht6weFNcnccVHlVoyRp4rBkZSRFmVwRAADoy9xyI+Xrr79+3DaEbniT99d+N03ghWfQyw0AADzLLaF71qxZx20zduxYd5wKOGlFB+v0za4KSVJqvwhNHJ5ickUAAKCvc0voBnzJh//tOJY7MIBpAgEAgGcRuuFXGpttWrWhUJIUGhKosydnm1wRAADwB4Ru+JX/bChSU4tzmsCZEzMVFRFickUAAMAfELrhNwzD+N4KlNxACQAAegehG35jy+5KFVfUS5JGDeynAemxJlcEAAD8RY9C9+7du91dB+BxH3293/WcXm4AANCbehS6Z8yY4e46AI+qrW/Rl1tLJUkxkSE6dXSqyRUBAAB/0qMVKVtaWvT555/LMIwTOm769Ok9OR1w0lbnHlCb3SFJOuuULAUHBZpcEQAA8Cc9Ct21tbWaP3/+cUO3xWJRYWGhDMNQYGCgWltbe1QkcDIMw9DHRwwtOXdqjonVAAAAf9Sj0J2UlKT8/PxjtmlsbNQNN9yggoICJSYmavHixT0qEDhZO/dXq+hgnSRp5IAEZaVEm1wRAADwNz0K3ZMmTTrm/l27dumyyy7T9u3bNWnSJC1ZskTZ2SxCAnN89BW93AAAwFw9upFy2bJlR9339ttva/Lkydq2bZvmz5+vtWvXErhhmsZmmz7fXCxJiggL0ulj002uCAAA+CO3zdNtt9t1xx136Morr1RLS4ueffZZvfTSSwoNDXXXKYAT9tmmYrW0OlegnDExU2GhPfrjDgAAwElxSwIpKSnR5Zdfri+//FLp6el65513NHXqVHe8NXBSPuIGSgAA4AVOuqf7P//5jyZOnKgvvvhCP/jBD5Sbm0vghlfIL67VnqIaSdKgzFgNzowztR4AAOC/Tip0L1q0SOedd57Ky8t1++23a9WqVUpOTnZXbcBJWUkvNwAA8BI9Gl5SU1Oja665Rh988IEiIiL0/PPPa968ecc8ZsyYMdq6dWuPigROVIvNrtUbD0iSQoIDNWNCpskVAQAAf9aj0D1hwgQVFhZq8ODBWrp0qUaPHn3cYyorK3tyKqBHvthSooYmmyTpjHHpigwPNrkiAADgz3oUuvfvd/7Zfs+ePRo3btxx2xuGIYvF0pNTAT3CCpQAAMCb9HhFynXr1nW7vWEY3FyJXlN2qEHf7j0kScpIitLIAQkmVwQAAPxdj26kDAgIUE5OTrcf/fv3V0CA26YEPymNjY26++67FRQUpAceeOC47T/88ENNnz5dycnJSklJ0ezZs4/5C0dJSYkWLlyo9PR0JScna/To0XrqqafkcDjc+ClwLJ8eHsstSbMmZ/FXFgAAYLoeJeH24SWePsbdVqxYoVGjRmnJkiWy2+3Hbf/CCy/oggsu0HnnnaeSkhIVFBRowIABOuOMM7Ry5cpO7QsLCzV58mTl5eUpNzdX5eXl+tOf/qRf//rXWrhwoSc+Er7HMAyt3lDkej1jIjdQAgAA8/UodIeEhHR4/eabbx73mFGjRvXkVG7z5Zdf6sYbb9TTTz+t++6777jti4uLddttt2n69On67W9/q6CgIIWHh+upp55SVlaWFixYoKampg7H3HLLLSotLdVrr72mtLQ0SdLs2bN111136dVXX9V7773nkc+G7+wuqlFJZYMkacygRCXHR5hcEQAAgJuWgb/zzjuPuX///v0qLy93x6l6bMiQIdq2bZsuuOCCbrV/8cUX1dTU1GkqxKCgIM2dO1fFxcVaunSpa3thYaHef/99TZkyRYMGDepwTPt7/OUvfznJT4HjObKX+8xJ9HIDAADv4JbQXVtbq3379nW5780339T48eNVX1/vjlP1WGJioiIjI7vdftWqVZKkSZMmddrXvu3IISbHaj98+HBFRERozZo1stlsJ1Q3uq/N7tBn3xRLkkKCAnTa2HSTKwIAAHByS+hua2vTeeedp7KyMte2hoYGzZ8/X1dddZWCgoIUFhbmjlP1mh07dkiSMjIyOu1r35aXl9et9haLRenp6bLZbNq7d68nyoWkjTvLZW1olSRNGZXK3NwAAMBruCV0JyQk6De/+Y3OOeccVVdXa/369Ro/frxee+01zZo1S1u2bFFMTIw7TtVrampqJKnL3vH2bdXV1d1qf7Rjvq+lpUVWq7XDA93XYWjJKVkmVgIAANBRj+bp/r5Vq1Zp5MiRqq6u1pQpU1RYWChJevTRR3X33XfLYrG4ZfaSRYsWqbW1tdvtr7vuOvXv3/+kz9tbHnnkET344INml+GTGpps+nqb8y8tMZEhmjgs2eSKAAAAvuOW0D1y5EhJ0i9+8QtVV1dr0aJFWrNmjU4//XRXm4ceekh/+MMfTuo8ixYtUkNDQ7fbz5w5s8ehOy4uThUVFWpoaFBsbGyHfe01xMfHd2h/5L7v6+qY77v33ns73JRqtVqVlUWPbXd8saVEtjbnXOjTx2coKNA75oUHAACQ3DS85EgPPfSQbrrpJj322GMyDMO1/aWXXjrp966vr5dhGN1+zJw5s8fnGjFihCTn1IHf175t+PDh3WpvGIZKSkoUHBzcaWaTI4WGhiomJqbDA92zOve7BXEYWgIAALxNj3q6AwMDu/fmQW7pSDfFrFmz9Nlnnyk3N1eTJ0/usC83N1eSdM4557i2nXXWWR32HSkvL0+NjY06++yzFRzMzX3uVl7dqK17KyVJ6YmRGpIVZ25BAAAA39OjVBwUFKTTTjut2+0Nw9BXX33Vk1OZZuHChXr00Ue1ePFi3Xjjja7tbW1teuedd5SRkaFLLrnEtT0nJ0dz5szRsmXLtHfv3g492osXL5Yk3Xrrrb33AfzImo0de7lZ9h0AAHibHoXuhIQErV69+oSOaV+h0VdkZmbqySef1E9/+lMtWrRI99xzj2w2m+68804VFhbqgw8+UEREx9UOn376aW3YsEHXXnutlixZotTUVK1YsUKPP/645s+fr4suusikT9N3GYbRYWjJTJZ9BwAAXqhHY7p///vfn/AxL7zwQk9O5VbTpk1Tamqqbr/9dknSY489ptTUVGVmdh3UbrjhBn3wwQdavny5UlNTlZ2drfz8fK1du1bnnntup/bZ2dlav369hg4dqokTJyo5OVm//OUvtWjRIreMaUdn+cW1KjpYJ0ka0T9Bqf26vwASAABAb7EYR97tCK9itVoVGxur2tpabqo8ihff+1bvrnEuOHTT3HGaPa2/uQUBAAB0gXnV4LMcDkOfH172PSjQojPGsew7AADwToRu+Kyd+6t1qLZZkjR+aLKiI0JMrggAAKBrhG74rLVbvpsTnV5uAADgzQjd8EkOh6H/bi6R5BxaMnW0b82OAwAA/AuhGz7p+0NLosJZdAgAAHivHoXu5557TuXl5e6uBeg2hpYAAABf0qPQfdNNNykjI0MzZszQ008/reLi4uMfBLgJQ0sAAICv6VHo3r9/v/73f/9Xdrtdt99+u3JycjRt2jT9+c9/VkFBgZtLBDpiaAkAAPA1PQrdWVlZ+sUvfqG1a9fqwIEDevLJJxUeHq577rlHgwYN0qRJk/Too49q165d7q4X0NrNDC0BAAC+xa0rUlZWVupf//qXlixZotWrV6utrU0jR47Uj3/8Y1166aUaPXq0u07lF1iRsjOHw9DCP3ysQ7XNCgq06LUHZ9PTDQAAvJ5bZy9JTEzUDTfcoBUrVujgwYN66aWXNGDAAD3yyCMaN26chg8frmXLlrnzlPAzDC0BAAC+yGNTBsbFxWn+/Pl67733VFFRoddff12jR4/Wzp07PXVK+IEjh5b8YDxDSwAAgG8I6o2TREVFad68eZo3b15vnA59lMNh6L9bvpu1ZMooZi0BAAC+gcVx4DPy9le5hpZMGMbQEgAA4DsI3fAZ7XNzS8xaAgAAfAuhGz6BoSUAAMCXEbrhExhaAgAAfBmhGz7hq2/LXM9PH8vQEgAA4FsI3fAJ67aVSpICAiyaPDLV5GoAAABOTK9MGShJpaWlstlsHbZlZ2f31unhw4oO1qm4okGSNHJAgmIiQ0yuCAAA4MR4NHRXVVXplltu0dKlSzsFbkmy2+2ePD36iHXbvhtaMnUUvdwAAMD3eHR4ye233659+/bpo48+UlhYmJYuXapnn31WAwcO1D//+U9Pnhp9yNcdQjezlgAAAN/j0Z7ulStX6r333tOUKVMUGBioESNG6OKLL1ZaWpr+93//V5dddpknT48+oKauRXn7qyRJWSnRSkuMNLkiAACAE+fRnu6mpiYlJydLkhISElRW5uyxHD16tDZu3OjJU6OPWL+9TIbhfH7qaIaWAAAA3+TRnu7Ro0fr22+/Vf/+/XXKKafoscceU1hYmJ5//nnl5OR48tToI44cWjKF8dwAAMBHeTR0//a3v1Vzs3NBk4cfflgXXXSRTj31VCUkJGjx4sWePDX6gObWNm3aVSFJio8O1dCseJMrAgAA6Jkehe5t27Zp5MiRslgsx2x3/vnnu54PGzZMO3fu1KFDh5SQkHDcY4EtuyvVanPOcDNlVKoCAvieAQAAvqlHY7rHjh2rhx9+uEcn7NevH4Eb3fLVt6Wu5wwtAQAAvqxHodswDDkcjg7b/va3v+nqq692S1GAw2Fo/faDkqTQkECNG5JkckUAAAA957bZS0pLS486TvuRRx7R1KlT3XUq+IFdhdWqqW+RJE0YmqTQ4ECTKwIAAOg5j04Z2K61tVUbNmzojVOhj2BBHAAA0Jf0SugGTtTX25zjuQMs0uSRKSZXAwAAcHII3fA6JRX1KjpYL0ka3j9BsVGhJlcEAABwcgjd8Dodh5YwawkAAPB9PV4c55lnntE333yjU045RaeccopqamrcWBb8WYfQPZrx3AAAwPf1KHRPmDBB27Zt07vvvqt33323w7zbV1xxhcaOHauxY8dq3Lhxys7Odlux6PvqGlu1Y98hSVJGUpQykqJMrggAAODk9Sh05+bmymazaevWrdq4caNyc3O1ceNGbdmyRW+//bbefvttVxCPiYlRRESEW4tG3/XNzgo5DOdzbqAEAAB9RY+HlwQHB2vixImaOHGirr/+ekmS3W7Xtm3blJub6wrimzdvVm1tLatQoltydx50PZ80PNnESgAAANynx6G7K4GBga6hJQsWLJAkORwO7dixQ7m5ue48FfogwzC0Ma9cknMVypED+plcEQAAgHu4NXR3JSAgQKNGjdKoUaM8fSr4uIJSq6rrnKtQjhmUqBBWoQQAAH0EUwbCa+Qe7uWWGFoCAAD6FkI3vMbGI0L3REI3AADoQwjd8AqNzTZtPzxVYFq/SKUnMlUgAADoOwjd8Apb9lTKfniuQHq5AQBAX0PohldgaAkAAOjLCN0wnWEYyt3pDN1BgQEaMyjR5IoAAADci9AN0xVX1Ku8qlGSNGpggsJDPT6TJQAAQK8idMN0HYaWDGPpdwAA0PcQumG69qElEvNzAwCAvonQDVO12Oz6dk+lJKlfbJiyU6NNrggAAMD9CN0w1ba9h9Ta5pAkTRyWLIvFYnJFAAAA7kfohqlydx50PZ80nPHcAACgbyJ0w1TtN1EGBFg0bmiSydUAAAB4BqEbpjlY1agD5fWSpGHZ8YoKDza5IgAAAM8gdMM0G5m1BAAA+AlCN0yz6YjQPWEYoRsAAPRdhG6Ywu4wtOXwVIHREcEanBlnbkEAAAAeROiGKfYV16qhySZJGjM4UQEBTBUIAAD6LkI3TLFlT4Xr+djBzFoCAAD6NkI3TLH58NASSRo7ONHESgAAADyP0I1e12Z3aHv+IUlSQkyoMpOjTK4IAADAswjd6HW7C2vU3GqXJI0ZlMTS7wAAoM8jdKPXdRjPPYShJQAAoO8jdKPXbWE8NwAA8DOEbvSqVptdOwqqJEnJCRFK7RdpckUAAACeR+hGr8rbXyVbm0OSNI5ebgAA4CcI3ehVW3Z/N7RkDKEbAAD4CUI3ehXjuQEAgD8idKPXNDbbtKuwWpKUkRSlfrHhJlcEAADQOwjd6DXb91XJ7jAkMVUgAADwL4Ru9JqtRwwtGTc4ycRKAAAAehehG73myEVxRg/qZ2IlAAAAvYvQjV5R39iqvcW1kqT+aTGKjQo1uSIAAIDeQ+hGr9i695AM53BuxnMDAAC/43ehu7GxUXfffbeCgoL0wAMPHLOtxWJRampql48//vGPXR5TUlKihQsXKj09XcnJyRo9erSeeuopORwOD3wa33Hk0BLGcwMAAH8TZHYBvWnFihX6+c9/LovFIrvd3q1jysrKuv3+hYWFmjZtmnJycpSbm6u0tDQtX75cl112mTZu3KhXXnmlh5X7vvabKAMs0qiBjOcGAAD+xW96ur/88kvdeOONevrpp3Xfffd55By33HKLSktL9dprryktLU2SNHv2bN1111169dVX9d5773nkvN6uuq5Z+8vqJEmDMuMUGR5sckUAAAC9y29C95AhQ7Rt2zZdcMEFHnn/wsJCvf/++5oyZYoGDRrUYd+8efMkSX/5y188cm5v9+2eQ67nrEIJAAD8kd8ML0lM9GzYW7VqlSRp0qRJnfYNHz5cERERWrNmjWw2m4KD/aun99v87+bnHkPoBgAAfshverp76v7779e4ceOUmpqq7Oxs/ehHP9LKlSs7tduxY4ckKSMjo9M+i8Wi9PR02Ww27d271+M1e5vt+6okOcdzj+ifYHI1AAAAvY/QfRyNjY366KOPVFpaqv/85z8KCAjQueeeq9///vcd2tXU1EiSIiMju3yf9u3V1dVHPVdLS4usVmuHh6+rb7Jpf5nzc/RPj1VEmH/18gMAAEg+Nrxk0aJFam1t7Xb76667Tv379+/x+UpLS5Wamup6PXjwYL311lsaPny47r//fl100UUaN25cj9//+x555BE9+OCDbns/b5BXUOWan3vkAHq5AQCAf/K50N3Q0NDt9jNnzjyp0H1k4G4XHBysiy++WE888YSWLl3qCt1xcXGSdNT62rfHx8cf9Xz33nuv7rzzTtdrq9WqrKysnpbvFXYUVLmejxzAVIEAAMA/+VTorq+vN7sESVJ6erokZ094uxEjRkiSiouLO7U3DEMlJSUKDg7uNLPJkUJDQxUa2reWR9++77uZS+jpBgAA/oox3Ufx7rvvdnnDpORcdVKSkpOTXdvOOussSVJubm6n9nl5eWpsbNSMGTP8auYSW5tDu/Y7x7CnJESoX2y4yRUBAACYg9B9FO+++66eeOKJTtvb2tq0bNkySdJFF13k2p6Tk6M5c+Zo3bp1nWYoWbx4sSTp1ltv9WDF3mdvcY1a2xyS6OUGAAD+jdB9DMuXL+8wjry0tFTXXHONdu/erTvvvFNTpkzp0P7pp59Wamqqrr32Wtfy8StWrNDjjz+u+fPndwjp/mB7PuO5AQAAJMliGO1zS/R906ZN0759+9TU1CSr1arIyEhFRUUpKChIBw4c6NC2qKhIb731lt5//33l5+erqalJhmFo0qRJuvHGG3XppZd2eY7i4mL97ne/00cffaS2tjYlJSXppz/9qW699VYFBJzY7zhWq1WxsbGqra1VTExMjz+3Wf7w0tf6epvzl4+n7z5TOam+9xkAAADcwa9Ct6/x5dBtGIauum+F6hpbFRUerH88NFsBARazywIAADAFw0vgEQfK61XX6JxTfcSABAI3AADwa4RueET70u8S47kBAAAI3fAI5ucGAAD4DqEbHrHjcE93cFCAhmTFmVsMAACAyQjdcLsqa7NKDzmnWRySFafgoECTKwIAADAXoRtut4Px3AAAAB0QuuF2jOcGAADoiNANtzsydI/oT+gGAAAgdMOtmlralF9ilSTlpEYrKiLE5IoAAADMR+iGW+3aXy2Hw7nIKeO5AQAAnAjdcCvGcwMAAHRG6IZbsRIlAABAZ4RuuI3d7lDefmfoTowNU1J8uMkVAQAAeAdCN9ymoNSq5la7JGnEgH6yWCwmVwQAAOAdCN1wm12F1a7nw3LiTawEAADAuxC64Ta7Cmtcz4dlE7oBAADaEbrhNjsP93QHBlg0ICPW5GoAAAC8B6EbbtHYbNOB8jpJ0oD0GIUGB5pcEQAAgPcgdMMtdhfVyHCuiaMhDC0BAADogNANt+hwEyWhGwAAoANCN9xid1GN6/lQQjcAAEAHhG64xc79zp7uiLAgZSRFmVwNAACAdyF046Qdqm1SlbVZkjQkK04BASyKAwAAcCRCN05aey+3xNASAACArhC6cdKOvImS0A0AANAZoRsnjZsoAQAAjo3QjZNidxjaXeTs6U6MC1dCTJjJFQEAAHgfQjdOyoHyOjW12CVJQ7PjzC0GAADASxG6cVJ2HXkTZRZDSwAAALpC6MZJ2XXkeO4cQjcAAEBXCN04Ke093QEWaXBmnLnFAAAAeClCN3qsubVNBWVWSVJ2aozCQ4NMrggAAMA7EbrRY/nFtXI4DEnOlSgBAADQNUI3eoxFcQAAALqH0I0e21VY43o+jJsoAQAAjorQjR7bebinOzQkUNkp0SZXAwAA4L0I3eiRmroWlVc1SnLOWhIYyLcSAADA0ZCU0CO7ir4bz81NlAAAAMdG6EaPHHkTJeO5AQAAjo3QjR7ZfcRNlCz/DgAAcGyEbpwwwzC0+/Dy77FRIUqKDze3IAAAAC9H6MYJq6xpVl1jqyRpUGacLBaLyRUBAAB4N0I3Ttje4hrX80EZseYVAgAA4CMI3Thh+cW1rucDCd0AAADHRejGCSN0AwAAnBhCN07Y3sOhOzw0SKkJkSZXAwAA4P0I3TghtfUtqqxpkuTs5Q4I4CZKAACA4yF044QcObSEmygBAAC6h9CNE8J4bgAAgBNH6MYJIXQDAACcOEI3Tkj7TZTBQQHKSok2uRoAAADfQOhGtzW1tKmksl6SlJMWo6BAvn0AAAC6g9SEbttXUivDcD7nJkoAAIDuI3Sj2xjPDQAA0DOEbnQboRsAAKBnCN3otvabKAMsUv+0GJOrAQAA8B2EbnSLrc2hwjKrJCkjOVphIUEmVwQAAOA7CN3olsIyq9rszrsouYkSAADgxBC60S2M5wYAAOg5Qje65cjQPSiT0A0AAHAiCN3olr1H9nSnE7oBAABOBKEbx+VwGNpX4gzdyQkRiooIMbkiAAAA30LoxnGVVNarudUuiZsoAQAAeoLQjePiJkoAAICTQ+jGcXW4iZLQDQAAcMII3TiuvfR0AwAAnBRCN47JMAxXT3dcVKgSYsJMrggAAMD3ELpxTJU1zbI2tEpy9nJbLBaTKwIAAPA9hG4cU35xjes5Q0sAAAB6htCNY8ovsbqesxIlAABAzxC6cUwFpd/dRNk/LcbESgAAAHwXoRvHtL+0TpIUEhSgtMQok6sBAADwTYRuHFWrza7SynpJUlZqtAIDuIkSAACgJwjdOKqig3VyGM7nOakMLQEAAOgpvwndmzZt0i233KKhQ4eqX79+iomJ0YQJE/THP/5RLS0tRz3uww8/1PTp05WcnKyUlBTNnj1b69atO2r7kpISLVy4UOnp6UpOTtbo0aP11FNPyeFweOJjedT+sjrX85zUaBMrAQAA8G1+Ebo3btyoiRMnav369Xr33Xd16NAhlZeX65prrtGvf/1rnXXWWWptbe103AsvvKALLrhA5513nkpKSlRQUKABAwbojDPO0MqVKzu1Lyws1OTJk5WXl6fc3FyVl5frT3/6k379619r4cKFvfFR3aqw7LuZS3K4iRIAAKDH/CJ0t/cyv/HGGxo5cqQkKSwsTHfeeafmzp2rL774Qi+99FKHY4qLi3Xbbbdp+vTp+u1vf6ugoCCFh4frqaeeUlZWlhYsWKCmpqYOx9xyyy0qLS3Va6+9prS0NEnS7Nmzddddd+nVV1/Ve++91wuf1n0KSo8I3QwvAQAA6DG/CN3p6elatGiRBg0a1GnfD37wA0nSF1980WH7iy++qKamJs2bN6/D9qCgIM2dO1fFxcVaunSpa3thYaHef/99TZkypdN52t/jL3/5i1s+T29pH14SGRakfrEs/w4AANBTfhO677333i73tQ8r6devX4ftq1atkiRNmjSp0zHt244cYnKs9sOHD1dERITWrFkjm83Wg0/Q+xqabKqscfbk56TFsPw7AADASfCL0H0s7T3c3+/R3rFjhyQpIyOj0zHt2/Ly8rrV3mKxKD09XTabTXv37nVP4R62v4yhJQAAAO4SZHYBZtq5c6eWLVuma6+9VlOmTOmwr6amRpIUGRnZ6bj2bdXV1d1qf7Rjvq+lpaXDTCpWq/WobT2NmUsAAADcx6dC96JFi7qcZeRorrvuOvXv37/Lfc3Nzbr22ms1evRoPfPMM26q8OQ88sgjevDBB80uQ5JUWMrMJQAAAO7ic6G7oaGh2+1nzpzZZehua2vT5ZdfLqvVqs8++0xRUZ2XN4+Li1NFRYUaGhoUGxvbYV97DfHx8R3aH7nv+7o65vvuvfde3Xnnna7XVqtVWVlZR23vSQVHDC/JZngJAADASfGp0F1fX3/S79Ha2qorrrhC+/bt06effqqkpKQu240YMUIVFRUqLi5Wenp6h33FxcWSnDdIHtn+yH1HMgxDJSUlCg4O7nIGlXahoaEKDQ094c/kboZhaH+pc3hJQkyoYiJDTK4IAADAt/nVjZQtLS269NJLVVhYqE8//VQpKSmSpLq6Om3durVD21mzZkmScnNzO71P+7ZzzjnHte2ss846avu8vDw1NjZqxowZCg4Ods+H8aCauhbVNTqH8XATJQAAwMnzm9Dd1NSkOXPmqKqqSv/5z386TBGYm5urOXPmdGi/cOFChYeHa/HixR22t7W16Z133lFGRoYuueQS1/acnBzNmTNH69at6zRDSft73Hrrre7+WB5RwHhuAAAAt/Kp4SU9VV9frwsvvFBffPGFbrvtNj3xxBMd9hcUFHQ6JjMzU08++aR++tOfatGiRbrnnntks9l05513qrCwUB988IEiIiI6HPP0009rw4YNuvbaa7VkyRKlpqZqxYoVevzxxzV//nxddNFFnvyYbsPMJQAAAO7lF6H7k08+0Zo1ayRJjz/+eJdtcnJyOm274YYblJGRoUceeURPPPGELBaLJk6cqLVr12rq1Kmd2mdnZ2v9+vX63e9+p4kTJ6qtrU1JSUlatGiRz/RyS1JhGT3dAAAA7mQxDMMwuwh0zWq1KjY2VrW1tYqJ6b3we+f/W6PdRTWyWKR/PnyBwkL94nczAAAAj/GbMd3oHofDUOFB5/CS1IRIAjcAAIAbELrRQXl1o1pa7ZKknDTGcwMAALgDoRsddJi5hOkCAQAA3ILQjQ72lxG6AQAA3I3QjQ7aV6KUpGyGlwAAALgFoRsdtPd0BwValJEUZXI1AAAAfQOhGy62NoeKy+slSZnJ0QoK5NsDAADAHUhVcCmuqJfd4Zy2PZuVKAEAANyG0A2X/UfMXNKflSgBAADchtANF2YuAQAA8AxCN1w6zFzC8BIAAAC3IXTDpb2nOzw0UMnxESZXAwAA0HcQuiFJampp08GqRklSdkqMAgIsJlcEAADQdxC6IUkqPGI8N0NLAAAA3IvQDUlSeXWT2ju3mbkEAADAvSyGYRhmF4GuWa1WxcbGqra2VjExng/CrTa7ig7WKS46VP1iwz1+PgAAAH8RZHYB8B4hwYEalBlndhkAAAB9DsNLAAAAAA8jdAMAAAAeRugGAAAAPIzQDQAAAHgYoRsAAADwMEI3AAAA4GGEbgAAAMDDCN0AAACAhxG6AQAAAA8jdAMAAAAeRugGAAAAPIzQDQAAAHgYoRsAAADwMEI3AAAA4GGEbgAAAMDDCN0AAACAhxG6AQAAAA8LMrsAHJ1hGJIkq9VqciUAAAA4lujoaFkslqPuJ3R7sbq6OklSVlaWyZUAAADgWGpraxUTE3PU/RajvTsVXsfhcKikpOS4vzm5i9VqVVZWloqKio75TQPvw7XzTVw338W1801cN9/lC9eOnm4fFhAQoMzMzF4/b0xMjNd+Q+PYuHa+ievmu7h2vonr5rt8+dpxIyUAAADgYYRuAAAAwMMI3XAJDQ3V/fffr9DQULNLwQni2vkmrpvv4tr5Jq6b7+oL144bKQEAAAAPo6cbAAAA8DBCNwAAAOBhhG4AAADAwwjdfqCtrU3//Oc/tWDBAo0aNUrp6elKSkrS+PHj9cc//lH19fWdjnE4HHryySc1evRoJScnKz09XQsXLlRpaakJnwBHuuWWW2SxWHTdddd1ub+5uVn333+/hgwZouTkZOXk5OiOO+5QbW1t7xbq52bOnKmEhASlpqZ2elx88cWd2nPdvIvdbtezzz6radOmKTs7W3FxcRo0aJCuvPJKbd68uUNbrp35XnnlFYWEhHT531tqaqrCwsIUGRnZ4Rium/ew2+169dVXNW3aNKWnpyslJUVjx47Vww8/3GVG8dlrZ6DPKy0tNSQZEyZMMLZu3WoYhmHYbDbjxRdfNAICAoypU6caNputwzFXX321ERERYXz44YeGYRhGSUmJMW3aNCMjI8MoKirq9c8Ap5UrVxoWi8WQZMyfP7/T/tbWVuPMM880kpOTjXXr1hmGYRi7du0yhgwZYowZM8aora3t5Yr914wZM4zVq1d3qy3Xzbs0NTUZ55xzjjFt2jTj22+/NQzD+W/mY489ZkgyXn75ZVdbrp13ePnll40ZM2Z0uc9msxnp6enG//zP/7i2cd28y/XXX29IMh566CGjpaXFcDgcxrJly4zw8HBjwoQJRktLi6utL187QrcfaA/dGzdu7LRv7ty5hiTjk08+cW1bsmSJIcm47777OrTdvXu3YbFYjEsuucTjNaOz6upqIzMz07jmmmuOGroff/xxQ5Lx0ksvddi+atUqQ5Lxi1/8opeqxYmEbq6bd7n99tuNfv36GdXV1Z32XXrppcb777/ves218w4ff/yxcffdd3e5r/1nWm5urmsb1817HDhwwJBkjB8/vtO+O+64w5BkLF682LXNl68dodsPtLa2GmvXrjUcDkenfb/85S8NScaSJUtc284880xDkrFjx45O7adMmWJYLBZ6u01w1VVXGRdeeKGxevXqo4bugQMHGoGBgUZdXV2H7Q6Hw0hJSTGio6ONpqamXqrYv51I6Oa6eY+DBw8aQUFBxu23396t9lw773fOOecYU6dO7bCN6+Y9vvrqK0OScfnll3fa98wzzxiSjEcffdS1zZevHWO6/UBwcLBOP/10WSyWTvvWrVunsLAwTZ06VZLU2tqqtWvXKjIyUsOGDevUftKkSTIMQ5988onH68Z3lixZoo8++kjPP//8Udvs27dP+fn5GjZsmKKiojrss1gsmjhxourq6vTVV195ulycAK6bd3n33XfV1tamU0899bhtuXbeb+/evfrkk0900003ubZx3bzL0KFDFRYWph07dnTa175tzJgxknz/2hG6/ZDD4dC+fft08803a+PGjXr55ZeVkZEhSdqzZ49sNpvS09O7DOnt7fLy8nq1Zn928OBB3Xjjjfq///s/paamHrVd+z9O7dfo+7h2ve+tt97SGWecoaysLKWmpmrGjBn661//Krvd7mrDdfMuubm5kqTY2Fjdd999GjVqlJKSkjRkyBDdeOONKiwsdLXl2nm/v/71r0pISNDll1/u2sZ18y7x8fH6y1/+ory8PP3mN7+R1WpVa2ur3n77bT3//POaN2+ezj//fEm+f+0I3X5m2bJliouL08CBA7VixQq98cYbuvLKK137a2pqJKnTXd7t2rdXV1d7vFY4XX/99Tr33HM1d+7cY7bj2nmfgoIC/fWvf1VRUZE2b96sM888UzfddJMuvPBC2Ww2SVw3b1NSUiJJmj9/vgoLC/XJJ5+opKREzzzzjN59911NnDhRO3fulMS183atra165ZVXtHDhQoWFhbm2c928z/XXX6933nlHb775pmJjYxUVFaUbbrhBixYt0htvvOFq5+vXjtDtZy688EJZrVaVl5frjjvu0Lx58zRnzpwup+SB+V588UVt3LhRTz/9tNml4AS9/fbbWrZsmUaNGiVJSklJ0QMPPKCf/OQnWrFihf7v//7P5ArRlaamJklSVFSUXnjhBaWlpSk4OFjnnnuu/vznP+vQoUO64447zC0S3bJkyRIdOnRIN954o9ml4BgMw9DPfvYzXXLJJbrxxhtVXV2turo6vfHGG/rTn/6k888/3xW2fR2h208lJSXp1ltv1e9//3stW7ZM9957ryQpLi5OktTQ0NDlce3b4+Pje6VOf1ZQUKA777xTL774Yrf+/+baeZekpCQFBgZ22v7jH/9YkjMQSFw3b9PeUzZr1iwFBQV12HfRRRdJkj755BM1Nzdz7bzcc889p/POO08DBw7ssJ3r5l1eeeUV/e1vf9PVV1+tX/3qV4qLi1NoaKjOP/98Pfnkk1q+fLnrF11fv3aEbj934YUXSnLePCRJgwcPVnBwsEpKSmQYRqf2xcXFkqThw4f3Wo3+6sMPP5TD4dB1113XYZGHSy+9VJJzvHD7tn//+98aMWKEpO+u0fdx7bxDenq6JLkWmuK6eZf+/ftLkhITEzvti4qKUmRkpNra2lRVVcW182I7duzQZ5991uEGynZcN++yYsUKSc5fdL+vfdu//vUvSb5/7QjdfuDTTz/V22+/3eW+iIgISdKhQ4ckSSEhITrjjDPU0NDgGrd4pNzcXFksFp199tmeKxiSpJtuukl1dXUqKyvr8Fi6dKkk6YorrnBtu/jiizVgwAANHDhQu3bt6jRcyDAMbdq0SdHR0d2alQEn55tvvjnqTDPtY4aTk5MlievmZaZPny5JKisr67SvsbFRDQ0NCgoKUkJCAtfOiz333HPq37+/6wa8I3HdvEv7Nehq8ob2bfX19bLb7T5/7QjdfuDTTz/V/fffL4fD0WnfypUrJck1ZaDkXGZckhYvXtyh7Z49e7R+/Xr96Ec/UmZmpgcrRk/dfPPNamtr6/RL1urVq1VWVqbrr7++ww1F8IxvvvlG99xzj5qbmzvta++xaR+qIHHdvMkFF1ygzMxMrVy5Uq2trR32LV++XJI0e/Zs1/Xg2nmfpqYm/f3vf9fPfvYzBQR0HXO4bt6jPSCvWbOm077PPvtMknTKKae4huv59LUzb4pw9Jb777/fkGQsXLjQKC4uNgzDMFpaWoy3337biI+PN2JiYowNGzZ0OOaqq64yIiIijOXLlxuG4VzVkmXgvcOxFsdpbW01Zs6c2Wl53KFDh3r98rh9ycsvv2xIMubOnWscOHDAMAzDqKurMx599FEjICDAOP300zss3sB18y4ff/yxERISYlx99dXGoUOHDIfDYXzxxRdGZmamkZaWZuTn57vacu28z0svvWSEhoYa5eXlR23DdfMe1dXVxrBhw4zAwEDjueeeM5qbmw2Hw2F8/vnnRnZ2thEaGmqsWbPG1d6Xrx2h2w9UVFQYzz77rHHOOecY/fv3N5KSkoyoqChj2LBhxs0332zs27ev0zF2u9144oknjJEjRxpJSUlGWlqacd1117lCO3rfH/7wByMlJcWIj483JBlhYWFGSkqKMWbMmA7tmpqajPvuu88YOHCgkZSUZGRlZRm33367UVNTY1Ll/sdqtRqvvPKKMWfOHNd1iI6ONk455RTjscceM1paWjodw3XzLhs2bDAuuugiIyEhwYiLizMGDBhg3HrrrUZZWVmntlw77zJ16lTjqquuOm47rpv3qKmpMX7zm98Yo0aNMqKjo43Y2FgjMzPTuOqqq4zNmzd3au+r185iGF3cLQcAAADAbRjTDQAAAHgYoRsAAADwMEI3AAAA4GGEbgAAAMDDCN0AAACAhxG6AQAAAA8jdAMAAAAeRugGAAAAPIzQDQAAAHgYoRsAAADwMEI3AKBHVqxYoenTpysiIkJJSUn61a9+JbvdbnZZAOCVCN0AgBP2xhtvaPbs2Tpw4IBuvPFGDR8+XH/605903333ndD7FBYWymKx6Nlnnz1u28bGRt1xxx0aOHCgRo8e3a1jAMBbWAzDMMwuAgDgO2pqapSenq5hw4bp888/V1RUlNra2jRmzBiVlJSooqJCISEh3XqvJ554QnfddZeKi4uVmpp6zLZXXXWVqqqq9Oabb6q0tFRnn322HnjgAV1//fXu+FgA4FH0dAMATsjSpUvV1NSke+65R1FRUZKkoKAgnX/++bJardq5c2e33+udd97RaaeddtzA3djYqDfffFMPPvigYmNjNXz4cN1yyy166aWXTuqzAEBvIXQDAE5Ifn6+JCkuLq7D9vbebavV2q33KSkp0Zdffqm5c+cet237H2UtFotrW0BAgPhjLQBfQegGAJyQ5ORkSdILL7zQIfQWFRVJkmJjY7v1Pv/6178kSZdeeulx20ZGRuriiy/Wgw8+qLq6Ou3atUtPP/20rr322hMtHwBMwZhuAEC3vPTSSyosLNSCBQs0efJkVVRUaMyYMTrttNOUmJioV155RaWlpaqrq1NERMRx3+/MM89UY2Ojvv76626dv7a2Vr/61a/00UcfKSIiQj/72c902223dej9BgBvRegGAHRLWlqaAgICVFxcrH379mnRokX673//q6KiIjU2NsrhcGjo0KHdGtNdUVGhtLQ0PfLII7r77rt7oXoAMBfDSwAAx1VdXa2ysjJNmjRJkjRgwAA9//zz2r59u+rq6vTJJ59IcvZed8e7774ru92uyy67zGM1A4A3IXQDAI6rvLxckpSTk9Pl/lWrVkmSLrzwwm693zvvvKPx48dr4MCB7ikQALwcoRsAcFzNzc2SpMDAwE777Ha7lixZopSUFP3whz887ntVV1dr9erV3Zq1BAD6iiCzCwAAeL+kpCRJUkFBQad9zz//vPLy8vT4448rKOj4P1bee+892Ww2hpYA8CvcSAkA6JbBgwfr4MGD2rt3r2vawG+++UZnnnmmBg4cqK+++krBwcHHfZ85c+YoPz9f27Zt83TJAOA1GF4CAOiW3/zmN6qvr9cZZ5yhhx9+WDfffLNOP/10JScn67333utW4K6rq9PKlSsZWgLA7xC6AQDdsnDhQj333HOyWCz6/e9/r+XLl+u2227Thg0blJGR0a33WLZsmVpaWhhaAsDvMLwEANBrLrvsMm3dulW7du0yuxQA6FXcSAkA6DWnnnqq5s+fb3YZANDr6OkGAAAAPIwx3QAAAICHEboBAAAADyN0AwAAAB5G6AYAAAA8jNANAAAAeBihGwAAAPAwQjcAAADgYYRuAAAAwMMI3QAAAICHEboBAAAAD/v/AaUMPVJXm20qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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": "5ea07b9a6b1a4925809c9a17e1ecb2b0",
       "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": "c7376fe85d154494b409643fff15705b",
       "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.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}