04d_Newton_arclength.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": [
    "# Softening behavior and its consequences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example problem\n",
    "\n",
    "We consider a minimal system to describe buckling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![minimal_buckling_model.png](https://raw.githubusercontent.com/nagelt/Numerical_Methods_Introduction/master/Images/minimal_buckling_model.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use our intuition for two limit cases. For a tiny mass and a stiff spring we expect the bar to stand upright, whereas for a huge mass and a soft spring, we expect it hanging down (either to the left or right). \n",
    "From these limits we can conclude that a stability loss of the upright position occurs. This stability loss, in terms of mechanics, corresponds to passing a state of zero stiffness. Thinking scalarly of F=ku we anticipate that solving for u=F/k fails when $k=0$. Indeed applying Newton's method the iterations do not converge (compare to exercise 04a starting Newton iterations on a point with a horizonal tangent).\n",
    "However non-converging iterations are neither a reliable nor an efficient indicator. So we go for an advanced method. In addition we obtain a glance into the exciting world of nonlinear systems.\n",
    "It goes without saying that you won't find nonlinear effects in a _linear world_, so make sure to enable nonlinear geometry in your FEM-code, if you want to analyse a structure for buckling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np #numerical methods\n",
    "from scipy.optimize import root # Newton's method outsourced\n",
    "import matplotlib.pyplot as plt # plotting\n",
    "from mpl_toolkits.mplot3d import Axes3D\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/Newton_arclength_plots.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principle of Minimum Potential Energy\n",
    "\n",
    "Since nature is lazy, all stable equilibria minimize potential energy (statics, conservative system). Note that maxima describe instable equilibria.\n",
    "For our system we have the potential energy of the rotational spring and of a mass in a gravitational field\n",
    "\\begin{align}\n",
    "\\Pi_\\mathrm{spring} &= \\frac{1}{2}c(\\varphi-\\varphi_0)^2, \\\\\n",
    "\\Pi_\\mathrm{grav} &= mgl\\sin \\varphi.\n",
    "\\end{align}\n",
    "The necessary condition for an equilibrium\n",
    "\\begin{equation}\n",
    "\\frac{\\partial \\Pi}{\\partial \\varphi}=0,\n",
    "\\end{equation}\n",
    "corresponds to an equilibrium of forces.\n",
    "From the second derivative we may identify the equilibrium as minimum $\\left(\\frac{\\partial^2 \\Pi}{\\partial \\varphi^2}>0\\right)$ or maximum of potential energy $\\left(\\frac{\\partial^2 \\Pi}{\\partial \\varphi^2}<0\\right)$, or we have to check higher derivatives $\\left(\\frac{\\partial^2 \\Pi}{\\partial \\varphi^2}=0\\right)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the spring relaxed at standing upright ($\\varphi_0=\\pi/2$) we define functions for potential energy and its derivatives. Note that we make potential energy nondimensional, thus leaving only one parameter\n",
    "\\begin{equation}\n",
    "\\frac{\\Pi}{\\Pi_\\mathrm{ref}} = \\frac{1}{2}(\\varphi-\\varphi_0)^2 + \\mu \\sin \\varphi,\n",
    "\\end{equation}\n",
    "with $\\mu=\\frac{mgl}{c}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi0 = 0.5*np.pi   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Epot = lambda phi, mu: 0.5*(phi-phi0)**2 + mu*np.sin(phi) # energy\n",
    "dEpot_dphi = lambda phi, mu: phi - phi0 + mu*np.cos(phi)  # force\n",
    "d2Epot_dphi2 = lambda phi, mu: 1 - mu*np.sin(phi)         # stability\n",
    "d2Epot_dphi_dmu = lambda phi, mu: np.cos(phi)             # for tangent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "146b900b2dd045bba9412894af5d6d3e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='$\\\\mu$', max=2.0), Output()), _dom_classes=('widget-…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_energy_landscape()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe how the energy landscape changes from single-well ($\\mu<1$) to a double well ($\\mu>1$). Basically this figure already explains everything, for a small load there is one equilibrium (upright) and for a huge load the upright position gets instable (maximum) and there occur two stable equilibria (minima), one to the left and the other to the right.\n",
    "\n",
    "The following graph shows the same, but plotting the potential energy and its derivative in $(\\phi,\\, \\mu)$ space. The gray plane indicates $\\text{d}E_\\text{pot}/\\text{d}\\varphi=0$. The trace of this intersection is the blue pitchfork, indicating the potential equilibrium positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "997922878145463685af118593b9afb5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='$\\\\mu$', max=2.0), Output()), _dom_classes=('widget-…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#%run plot_functions/Newton_arclength_plots.ipynb\n",
    "plot_energy_and_pitchfork()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution with the pseudo arc-length method\n",
    "\n",
    "The equation (equilibrium) to solve is transcendental\n",
    "\\begin{equation}\n",
    "f(\\varphi) = \\varphi-\\varphi_0 + \\mu \\cos \\varphi = 0,\n",
    "\\end{equation}\n",
    "so we resort to numerical methods.\n",
    "Newtons method, however, runs into trouble for $\\mu \\to 1$, because then the tangent \n",
    "\\begin{equation}\n",
    "\\frac{\\mathrm{d}f}{\\mathrm{d}\\varphi} = 1 - \\mu \\sin \\varphi\n",
    "\\end{equation}\n",
    "at the equilibrium ($\\varphi=\\varphi_0=\\frac{\\pi}{2}$) approaches zero.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it is time for the clever move, instead of solving for $\\varphi$ with the fixed parameter $\\mu$, we think of a solution curve in $\\varphi$-$\\mu$ space. Starting from a known point of this curve, e.g. the trivial unloaded case $\\mu=0$, we progress pointwise by imposing a fixed distance (called pseudo arc-length because coordinate $\\varphi$ and parameter $\\mu$ do not measure a spatial distance) from one point to another (index _p_ denotes previous point)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_squared = lambda phi,mu, phi_p,mu_p: 0.1*(phi-phi_p)**2 + (mu-mu_p)**2   # pseudo arc-length squared\n",
    "d_distance_squared_dphi = lambda phi,mu, phi_p,mu_p: 2*0.1*(phi-phi_p)  # for tangent\n",
    "d_distance_squared_dmu = lambda phi,mu, phi_p,mu_p: 2*(mu-mu_p)   # for tangent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the weigthing for the arc-length required some divine inspiration, we are not teaching here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The extended system (two equations for two unknowns), we have to solve, reads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sys(z, zp, ds):\n",
    "    phi = z[0]\n",
    "    mu = z[1]\n",
    "    phi_p = zp[0]\n",
    "    mu_p = zp[1]\n",
    "    \n",
    "    r0 = dEpot_dphi(phi, mu)   # equilibrium\n",
    "    r1 = distance_squared(phi, mu, phi_p, mu_p) - ds**2 # distance constraint\n",
    "    \n",
    "    res = np.array([r0, r1])\n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And its tangent reads "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsys(z, zp, ds):\n",
    "    phi = z[0]\n",
    "    mu = z[1]\n",
    "    phi_p = zp[0]\n",
    "    mu_p = zp[1]\n",
    "    \n",
    "    J00 = d2Epot_dphi2(phi, mu)\n",
    "    J01 = d2Epot_dphi_dmu(phi, mu)\n",
    "    J10 = d_distance_squared_dphi(phi, mu, phi_p, mu_p)\n",
    "    J11 = d_distance_squared_dmu(phi, mu, phi_p, mu_p)\n",
    "    \n",
    "    J = np.matrix([[J00, J01], [J10, J11]])\n",
    "    return J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from the unloaded initial state we increase the load (direction guided by start value of the iterations). In addition we observe the stability of the calculated states (from nonlinear systems theory it is known that a change of stability typically indicates branching of solutions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAImCAYAAACLl+YVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQXklEQVR4nO3de3xU1b3///dkMFeSgQBKIEEiEVQSrhK5GIFK/VJUwBiLWqqB/qoo4AWPcECx3mpsK+USPXK0CAKKIkYQaqtYAgaTCiohhIsCJ2jAiERghlsCTOb3x24iYyYwk8ueSfJ6Ph7zGLLW2nt/hvLYvruyZm2Ly+VyCQAAAECDC/J3AQAAAEBzQfgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4NpHL5ZLD4RDPNQIAAGieCN8mOnbsmGw2m44dO+bvUgAAAOAHhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJC38XQAAoHFyOqWcHKmkRIqJkVJSJKvV31UBQGBj5hsA4LOsLCkhQRo6VLrzTuM9IcFoBwDUjPANAPBJVpaUliYlJUl5edKxY8Z7UpLRTgAHgJpZXC6Xy99FNBcOh0M2m012u11RUVH+LgcAfOZ0GjPcSUnSypVS0DlTOBUV0ujRUmGhtHs3S1AAwBNmvgEAXsvJkfbtk2bMcA/ekvHz9OlSUZExDgBQHeEbAOC1khLjPTHRc39le+U4AIA7wjcAwGsxMcZ7YaHn/sr2ynEAAHes+TYRa74BNHas+QaAumHmGwDgNatVmjVLWrPGCNrn7nYyerTR/sILBG8AqAkP2QEA+CQ1VVqxQnrkEWngwJ/a4+ON9tRU/9UGAIGOZScmYtkJgKaEJ1wCgO+Y+QYA1IrVKg0Z4u8qAKBxYc03AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgkoAP32VlZXrssccUHBys9PR0n49fv369LBbLBV9DhgxxO27RokXnHT969Oh6+XwA0Fg5ndL69dKyZca70+nvigAg8LXwdwHnk52drXvvvVelpaU6c+ZMnc7VrVs3j+0VFRXavXu3rrzyymp9YWFh6tSpk8fjYmNj61QPADRmWVnSI49I+/b91Na5szRrlpSa6q+qACDwBWz4XrZsme677z798Y9/VHh4uMaPH1+n8+3atctj+8qVK3XLLbdo4sSJ1fqSk5O1fv36Ol0XAJqarCwpLU266SZj1jsxUSoslJ57zmhfsYIADgA1CdhlJ/Hx8dqxY4cmTpwoi8VS6/O0a9dOt956a4398+bN09ChQ5WYmFjrawBAc+F0GjPeN90krVwp9e8vtWxpvK9cabT/13+xBAUAahKwM9/9+/evl/N0795dK1as8Ni3fft2ZWdnKysrq16uBQBNXU6OsdRk2TIp6GfTN0FB0vTp0sCBxriffZUGAKAAnvk2Q2Zmpi699FKNHDnSY39paakmTJig+Ph4hYaGqm3btho+fLjef/99kysFgMBQUmK81/TLwsr2ynEAAHfNNnwfPXpUS5cu1f333y+r1epxzM6dO9W2bVtt2LBBDodD2dnZkqRRo0ZpypQpF7xGeXm5HA6H2wsAGrOYGOO9sNBzf2V75TgAgLtmG74XLFigiooK/e53v/PYP2zYMOXl5enZZ59Vp06dFBwcrKSkJK1cuVIdOnTQ7NmztWbNmvNeIyMjQzabreoVFxfXEB8FAEyTkmLsavLcc1JFhXtfRYWUkSHFxxvjAADVNcvwXVFRoZdeekl33nmn2rRp43FMbGyskpOTq7WHhoZW7byycOHC815n+vTpstvtVa/i4uK6Fw8AfmS1GtsJrlkjjR4t5eVJx44Z76NHG+0vvGCMAwBUF7BfuGxIa9asUVFRkSZNmlSr4xMSEiTVvH1hpZCQEIWEhNTqGgAQqFJTje0EH3nE+HJlpfh4thkEgAtpluF73rx5SklJUa9evWp1fMXPf9cKAM1Maqo0apSxq0lJibHGOyWFGW8AuJBmF7537typf/3rX1q+fPl5x7Vv316vvfaaRowYUa1v7969kqQrrriiQWoEgMbAamU7QQDwVZNZ852dna3k5GRlZmaed9y8efPUsWNH3XLLLecdd/DgQa1evbpae1lZmV577TVJUnp6eq3rBQAAQPPTZML3rFmztHnzZs2YMaPGMXa7XUuWLNF9992nFi0uPOn/yiuv6Mknn9Q333yj06dPq6CgQKNGjVJJSYkmT56sm2++uT4/AgAAAJq4gA7fFotFFotF48aNkyS9/vrrVW2LFi1yGztmzBhFRkZW7UTiyWuvvaazZ8/qnnvuueC1d+7cqSeeeEIfffSRBgwYoJYtW2rYsGGyWq1atWqV5s2bV6fPBgAAgObH4nK5XP4uorlwOByy2Wyy2+2KiorydzkAAAAwWUDPfAMAAABNCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMMmFn7EOAIAHTqeUkyOVlEgxMVJKimS1+rsqAAhszHwDAHyWlSUlJEhDh0p33mm8JyQY7QCAmhG+AQA+ycqS0tKkpCQpL086dsx4T0oy2gngAFAzi8vlcvm7iObC4XDIZrPJbrcrKirK3+UAgM+cTmOGOylJWrlSCjpnCqeiQho9WioslHbvZgkKAHjCzDcAwGs5OdK+fdKMGe7BWzJ+nj5dKioyxgEAqiN8AwC8VlJivCcmeu6vbK8cBwBwR/gGAHgtJsZ4Lyz03F/ZXjkOAOCONd8mYs03gMaONd8AUDfMfAMAvGa1SrNmSWvWGEH73N1ORo822l94geANADXhITsAAJ+kpkorVkiPPCINHPhTe3y80Z6a6r/aACDQsezERCw7AdCU8IRLAPAdM98AgFqxWqUhQ/xdBQA0Lqz5BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAEzSwt8FAAAaJ6dTysmRSkqkmBgpJUWyWv1dFQAENma+AQA+y8qSEhKkoUOlO+803hMSjHYAQM0I3wAAn2RlSWlpUlKSlJcnHTtmvCclGe0EcAComcXlcrn8XURz4XA4ZLPZZLfbFRUV5e9yAMBnTqcxw52UJK1cKQWdM4VTUSGNHi0VFkq7d7MEBQA8YeYbAOC1nBxp3z5pxgz34C0ZP0+fLhUVGeMAANURvgEAXispMd4TEz33V7ZXjgMAuAv48F1WVqbHHntMwcHBSk9Pr9U50tPTZbFYanzNmTOnxmOXLFmia665RhEREYqOjtaIESOUl5dXuw8DAI1cTIzxXljoub+yvXIcAMBdQIfv7Oxs9ejRQy+//LLOnDlTp3O1adNG3bp18/iKjo72eMyUKVN011136aabbtJ3332n/Px8hYaGKiUlRe+9916d6gGAxiglRercWXruOWON97kqKqSMDCk+3hgHAKguYPf5XrZsme677z798Y9/VHh4uMaPH1+n802aNElPPvmk1+PXrl2r2bNnKy0tTTNnzpQk2Ww2vfHGG+rWrZvGjRunlJQUtW3btk51AUBjYrVKs2YZu5qMHm2s8U5MNGa8MzKkNWukFSv4siUA1CRgZ77j4+O1Y8cOTZw4URaLxfTrP//885Kke+65x609LCxMY8eOld1u1/z5802vCwD8LTXVCNjbtkkDB0pRUcZ7YaHRnprq7woBIHAFbPju37+/OnTo4JdrOxwObdiwQRaLRQMGDKjWP2jQIEnSqlWrzC4NAAJCaqq0Z4+UnS29+abxvns3wRsALiRgw3d9y8/P18iRIxUTE6OQkBDFxcVp7NixKigoqDa2sLBQTqdT7dq1U8uWLav1x8fHV42r+PmiRwBoJqxWacgQ6Y47jHeWmgDAhTWb8J2Tk6PRo0dr+/btOnr0qBYsWKB///vf6tevn9555x23sSX/2SOrpi9itm7dWpKxE8vRo0cbtG4AAAA0Hc0ifD/88MPatGmTxo8fr+joaIWFhemGG27Qu+++qzNnzmjcuHE6ePBg1fiTJ09KkoKDgz2eLyQkpNpYT8rLy+VwONxeAAAAaL6aRfju2bOnunTp4rF90KBBOnHihJYvX17VHh4eLkk6ffq0x/OVl5dXG+tJRkaGbDZb1SsuLq62HwEAAABNQLMI3+eTkJAgSdq1a1dVW/v27SVJhw8f9njMkSNHJEmhoaFq1apVjeeePn267HZ71au4uLieqgYAAEBjFLD7fJvF0xcmk5KSFBQUpEOHDunEiROKiIhw6y8qKpIkJSYmKiio5v//EhIS4rZEBQAAAM1bk5/5zs3NVfv27Wv8YuTevXslSVdccUVVW1RUlAYPHiyXy+XxUfK5ubmSpJEjR9Z/wQAAAGiymkz4zs7OVnJysjIzM93aT58+rYMHD+rDDz+sdkx+fr5yc3MVHh6u2267za1v6tSpkqRXX33Vrb2srExLly5VVFSUJkyYUM+fAgAAAE1Zkwnfs2bN0ubNmzVjxgy39sqnY06cOFFLlizRDz/8oLKyMn300UdKS0tTixYttGDBgqp13pWGDx+uBx54QMuXL9ezzz4rh8Oh4uJijR07Vvv379eCBQvUrl070z4fAAAAGj+Ly+Vy+buImpzvsfILFy5Uenp61c9LlizRxIkTNW7cOM2dO7eq3eVyaePGjXrrrbeUm5urb7/9VsePH1dMTIwGDx6sKVOmqGfPnjVeZ/HixcrMzNSOHTt00UUXacCAAZo5c6YGDhzo8+dxOByy2Wyy2+2Kiory+XgAAAA0bgEdvpsawjcAAEDz1mSWnQAAAACBjvANAAAAmKTZ7/MNAKgdp1PKyZFKSqSYGCklRbJa/V0VAAQ2Zr4BAD7LypISEqShQ6U77zTeExKMdgBAzQjfAACfZGVJaWlSUpKUlycdO2a8JyUZ7QRwAKgZu52YiN1OADR2Tqcxw52UJK1cKQWdM4VTUSGNHi0VFkq7d7MEBQA8YeYbAOC1nBxp3z5pxgz34C0ZP0+fLhUVGeMAANURvgEAXispMd4TEz33V7ZXjgMAuCN8AwC8FhNjvBcWeu6vbK8cBwBwx5pvE7HmG0Bjx5pvAKgbZr4BAF6zWqVZs6Q1a4ygfe5uJ6NHG+0vvEDwBoCa8JAdAIBPUlOlFSukRx6RBg78qT0+3mhPTfVfbQAQ6Fh2YiKWnQBoSnjCJQD4jplvAECtWK3SkCH+rgIAGhfWfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACZp4e8CAACNk9Mp5eRIJSVSTIyUkiJZrf6uCgACGzPfAACfZWVJCQnS0KHSnXca7wkJRjsAoGaEbwCAT7KypLQ0KSlJysuTjh0z3pOSjHYCOADUzOJyuVz+LqK5cDgcstlsstvtioqK8nc5AOAzp9OY4U5KklaulILOmcKpqJBGj5YKC6Xdu1mCAgCeMPMNAPBaTo60b580Y4Z78JaMn6dPl4qKjHEAgOoI3wAAr5WUGO+JiZ77K9srxwEA3BG+AQBei4kx3gsLPfdXtleOAwC4Y823iVjzDaCxY803ANQNM98AAK9ZrdKsWdKaNUbQPne3k9GjjfYXXiB4A0BNeMgOAMAnqanSihXSI49IAwf+1B4fb7SnpvqvNgAIdCw7MRHLTgA0JTzhEgB8x8w3AKBWrFZpyBB/VwEAjQtrvgEAAACTEL4BAAAAkwR8+C4rK9Njjz2m4OBgpaen1+oc69at0+9+9zt17dpVoaGhCg8PV/fu3TV16lSVlpZ6PGbRokWyWCw1vkaPHl37DwUAAIBmKaDXfGdnZ+vee+9VaWmpzpw5U6tzvPTSS5o0aZJ69Oih+fPnq1+/fnI4HHrjjTf0+OOPa+nSpfrkk0+UkJBQ7diwsDB16tTJ43ljY2NrVQ8AAACar4AN38uWLdN9992nP/7xjwoPD9f48eNrdZ5Tp04pODhYa9asUVxcnCQpMjJSU6dO1eHDh/WnP/1JDzzwgD744INqxyYnJ2v9+vV1+RgAAABAlYBddhIfH68dO3Zo4sSJslgstT5P+/btdccdd1QF73ONHDlSkrR27Vo5nc5aXwMAAADwRsDOfPfv379ezjN27FiNHTvWY5/NZpMkWSwWsd05AAAAGlrAznybYdeuXZKkQYMGqUWL6v8/pLS0VBMmTFB8fLxCQ0PVtm1bDR8+XO+//77ZpQIAAKAJaNbhe/HixbJYLJo5c6bH/p07d6pt27basGGDHA6HsrOzJUmjRo3SlClTLnj+8vJyORwOtxcAAACar2Ybvv/5z3/q/fff10MPPaRf/OIX1fqHDRumvLw8Pfvss+rUqZOCg4OVlJSklStXqkOHDpo9e7bWrFlz3mtkZGTIZrNVvTytOwcAAEDz0SzD91dffaW77rpLqamp+stf/uJxTGxsrJKTk6u1h4aGVu28snDhwvNeZ/r06bLb7VWv4uLiuhcPAACARitgv3DZUIqKivTLX/5SAwcO1LJly2S1Wn0+R+We4JVrxmsSEhKikJCQWtUJAACApqdZzXzv3r1b1113na699lqtWLFCwcHBtTpPRUVFPVcGAACA5qDZhO8dO3bouuuu0w033KClS5e67W7ywgsvVFsS0r59e48P3pGkvXv3SpKuuOKKhisYAAKc0ymtXy8tW2a887gEALiwJhO+s7OzlZycrMzMzGp9+fn5Gjx4sNLS0vS3v/1NQUHuH/vRRx+tCtSVDh48qNWrV1c7V1lZmV577TVJUnp6ev19AABoRLKypIQEaehQ6c47jfeEBKMdAFCzJrPme9asWdq8ebN27typyZMnV7V//vnnuuGGG1ReXq5Dhw7pjjvu8Pqcr7zyii655BKNGzdOMTEx2rVrlx599FGVlJRo8uTJuvnmmxviowBAQMvKktLSpJtuMma9ExOlwkLpueeM9hUrpNRUf1cJAIHJ4grgRzue77HyCxcudJt5XrJkiSZOnKhx48Zp7ty5Ve0PPfSQ2881yc7O1pAhQ6p+3rVrl95++219+OGH2rdvn0pLS9WqVStdffXVmjBhQtWj6X3hcDhks9lkt9sVFRXl8/EA4G9OpzHDnZQkrVwpnfuLxIoKafRoI4jv3i3V4vvsANDkBXT4bmoI3wAau/XrjSUmeXlS//7V+/PypIEDpexs6Zz5DADAfzSZNd8AgIZXUmK8JyZ67q9srxwHAHBH+AYAeC0mxngvLPTcX9leOQ4A4I5lJyZi2QmAxo413wBQN8x8AwC8ZrVKs2ZJa9YYQTsvTzp2zHgfPdpof+EFgjcA1KTJbDUIADBHaqqxneAjjxhfrqwUH882gwBwISw7MRHLTgA0JU6nlJNjfLkyJkZKSWHGGwAuhJlvAECtWK1sJwgAvmLNNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYJIW/i4AANA4OZ1STo5UUiLFxEgpKZLV6u+qACCwMfMNAPBZVpaUkCANHSrdeafxnpBgtAMAakb4BgD4JCtLSkuTkpKkvDzp2DHjPSnJaCeAA0DNLC6Xy+XvIpoLh8Mhm80mu92uqKgof5cDAD5zOo0Z7qQkaeVKKeicKZyKCmn0aKmwUNq9myUoAOAJM98AAK/l5Ej79kkzZrgHb8n4efp0qajIGAcAqI7wDQDwWkmJ8Z6Y6Lm/sr1yHADAHeEbAOC1mBjjvbDQc39le+U4AIA71nybiDXfABo71nwDQN0w8w0A8JrVKs2aJa1ZYwTtc3c7GT3aaH/hBYI3ANSEh+wAAHySmiqtWCE98og0cOBP7fHxRntqqv9qA4BAx7ITE7HsBEBTwhMuAcB3zHwDAGrFapWGDPF3FQDQuLDmGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwScCH77KyMj322GMKDg5Wenp6nc7z1FNPqVu3bgoNDVWHDh00fvx4HThw4LzHLVmyRNdcc40iIiIUHR2tESNGKC8vr9Z1AAAAoPkK6PCdnZ2tHj166OWXX9aZM2dqfZ6ysjJdf/31mjVrljIyMnTkyBGtXr1aeXl56t27t/bs2ePxuClTpuiuu+7STTfdpO+++075+fkKDQ1VSkqK3nvvvVrXAwAAgOYpYMP3smXLdMstt+jBBx/UrFmz6nSup59+Wrm5ucrIyFBqaqrCwsLUt29fLV++XKWlpR5n1NeuXavZs2crLS1NM2fOlM1mU6dOnfTGG2+oQ4cOGjdunEpLS+tUFwAAAJqXgA3f8fHx2rFjhyZOnCiLxVLr85w6dUqZmZkKDg7W3Xff7daXlJSk/v3769NPP9XGjRvd+p5//nlJ0j333OPWHhYWprFjx8put2v+/Pm1rgsAAADNT8CG7/79+6tDhw51Ps+6det0/Phx9ezZUy1btqzWP2jQIEnSqlWrqtocDoc2bNggi8WiAQMGeHUMADQ3Tqe0fr20bJnx7nT6uyIACHwt/F1AQ9u6daskYybdk8r2goKCqrbCwkI5nU5dfPHFHgN75TGFhYWqqKhQUFDA/n8YAGgQWVnSgw9K+/f/1BYbK82dK6Wm+q8uAAh0TT41lpSUSJKio6M99rdu3dptnC/HlJWV6ejRozVeu7y8XA6Hw+0FAI1dVpaUluYevCXpwAGjPSvLP3UBQGPQ5MP3yZMnJUnBwcEe+0NCQtzG+XLMz4/7uYyMDNlstqpXXFycb8UDQIBxOo0Zb5erel9l20MPsQQFAGrS5MN3eHi4JOn06dMe+8vLy93G+XLMz4/7uenTp8tut1e9iouLfSseAAJMTk71Ge9zuVxScbExDgBQXZNf892+fXtJ0uHDhz32HzlyRJIUExPj8zGhoaFq1apVjdcOCQlxmyUHgMbunBV69TIOAJqbJj/z3atXL0lSUVGRx/7K9h49elS1JSUlKSgoSIcOHdKJEydqPCYxMZEvWwJoVs6Zp6iXcQDQ3DT55Dh06FBFRESooKDAY5DOzc2VJI0cObKqLSoqSoMHD5bL5fL4KHlPxwBAc5CSYuxqUtPjFywWKS7OGAcAqK7JhO/s7GwlJycrMzPTrT08PFyTJk1SeXm5Fi9e7Na3fft25eXlacCAAUr52X8ppk6dKkl69dVX3drLysq0dOlSRUVFacKECQ3wSQAgcFmtxnaCUvUAXvnznDnGOABAdU0mfM+aNUubN2/WjBkzqvU98cQT6t+/v6ZNm6b33ntPp06d0pdffqlf//rXatOmjRYtWlTtmOHDh+uBBx7Q8uXL9eyzz8rhcKi4uFhjx47V/v37tWDBArVr186ETwYAgSU1VVqxQurY0b09NtZoZ59vAKiZxeXytGFUYDjfY+UXLlyo9PT0qp+XLFmiiRMnaty4cZpbOS1zjrKyMmVkZOjNN9/Ut99+q+joaA0fPlzPPPOMYmNja7zO4sWLlZmZqR07duiiiy7SgAEDNHPmTA0cONDnz+NwOGSz2WS32xUVFeXz8QAQSJxOY1eTkhJjjXdKCjPeAHAhAR2+mxrCNwAAQPPWZJadAAAAAIGO8A0AAACYxKfwnZeX53G7PgAAAAAX5tMTLgcNGqSgoCB16dJFvXr1cnud+4TIe+65R3379tW9995b7wUDAAAAjZVPX7j8zW9+o61bt+qrr76S0+k0TvCfHUnatm2rXr16qUuXLnrnnXfUokULlfB8YTd84RIAAKB5q9VuJ+Xl5dq2bZvy8/OrXps3b9aZM2dksVjkcrkUFxenb775piFqbrQI3wAAAM1bvW01aLfbtXDhQj322GNKSkrS4sWL1bVr1/o4dZNB+AYAAGje6m23E5vNpoceekhvvfWWNm/erG3bttXXqQEAAIAmoUEesnPVVVcpLCxMX3zxRX2fulFj5hsAAKB5a5B9vi+77DLt3LmzIU4NAAAANFo+bTX44IMPqnfv3urdu7e6d++uFi08H7537161b9++XgoEAAAAmgqfwndmZmbV1oIXXXSRunfvXhXGExMTFRISoqVLl+rrr7/Wyy+/3CAFAwAAAI2VT+E7KytLW7ZsqfaqDOSVEhIS9OOPP+qf//yn+vTpo4svvrheiwYAAAAaozp94bK0tLRaEN+zZ48qKircAnlMTIz69u2rVatW1UvRjRVfuATQlDidUk6OVFIixcRIKSmS1ervqgAgsNX7bicnTpzQ1q1b3QL59u3bdebMmaqnYjZXhG8ATUVWlvTgg9L+/T+1xcZKc+dKqan+qwsAAl2DbDX4c2fPntX27dvVs2fPhr5UQCN8A2gKsrKktDTp5//1qPyF54oVBHAAqIkp4RsGwjeAxs7plDp3dp/xPpfFYsyAFxWxBAUAPGmQfb4BAE1TTk7NwVsyZsOLi41xAIDqCN8AAK+VlNTvOABobgjfAACvxcTU7zgAaG4I3wAAr6WkGGu6f/Z4hyoWixQXZ4wDAFRH+AYAeM1qNbYTlKoH8Mqf58zhy5YAUBPCNwDAJ6mpxnaCHTu6t8fGss0gAFwIWw2aiK0GATQlPOESAHzXwt8FAAAaJ6tVGjLE31UAQOPCshMAAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkLfxdAACgcXI6pZwcqaREiomRUlIkq9XfVQFAYCN8AwB8lpUlPfigtH//T22xsdLcuVJqqv/qAoBAF9DLTsrKyvTUU0+pW7duCg0NVYcOHTR+/HgdOHDA63OsX79eFovlgq8hQ4a4Hbdo0aLzjh89enT9flgAaCSysqS0NPfgLUkHDhjtWVn+qQsAGoOAnfkuKyvT9ddfr23btmnRokX61a9+pR07dmjs2LHq3bu3cnNzlZCQ4PX5unXr5rG9oqJCu3fv1pVXXlmtLywsTJ06dfJ4XGxsrNfXBoCmwuk0Zrxdrup9LpdksUgPPSSNGsUSFADwJGDD99NPP63c3Fy9+OKLSv3P7zD79u2r5cuXq2fPnkpPT9fGjRu9Pt+uXbs8tq9cuVK33HKLJk6cWK0vOTlZ69evr1X9ANAU5eRUn/E+l8slFRcb4372C0UAgAJ02cmpU6eUmZmp4OBg3X333W59SUlJ6t+/vz799FOvwne7du1066231tg/b948DR06VImJiXWuGwCaupKS+h0HAM1NQIbvdevW6fjx4+rZs6datmxZrX/QoEGSpFWrVl3wXN27d9eKFSs89m3fvl3Z2dmaPHly3QoGgGYiJqZ+xwFAcxOQ4Xvr1q2SpPj4eI/9le0FBQV1uk5mZqYuvfRSjRw50mN/aWmpJkyYoPj4eIWGhqpt27YaPny43n///TpdFwAaq5QUY1cTi8Vzv8UixcUZ4wAA1QVk+C75z+8ro6OjPfa3bt3abVxtHD16VEuXLtX9998vaw3fCtq5c6fatm2rDRs2yOFwKDs7W5I0atQoTZky5YLXKC8vl8PhcHsBQGNmtRrbCUrVA3jlz3Pm8GVLAKhJQIbvkydPSpKCg4M99oeEhLiNq40FCxaooqJCv/vd7zz2Dxs2THl5eXr22WfVqVMnBQcHKykpSStXrlSHDh00e/ZsrVmz5rzXyMjIkM1mq3rFxcXVul4ACBSpqdKKFVLHju7tsbFGO/t8A0DNAjJ8h4eHS5JOnz7tsb+8vNxtnK8qKir00ksv6c4771SbNm08jomNjVVycnK19tDQUI0fP16StHDhwvNeZ/r06bLb7VWv4uLiWtULAIEmNVXat0/KzpbefNN4LyoieAPAhQTkVoPt27eXJB0+fNhj/5EjRyRJMbX8Rs+aNWtUVFSkSZMm1er4yv3Fa9q+sFJISEjVLD0ANDVWK9sJAoCvAnLmu1evXpKkoqIij/2V7T169KjV+efNm6eUlJSq6/iqoqKiVscBAACgeQvI8D106FBFRESooKBAJ06cqNafm5srSTXuUnI+O3fu1L/+9a8Lbi/Yvn17ffDBBx779u7dK0m64oorfL4+AAAAmq+ADN/h4eGaNGmSysvLtXjxYre+7du3Ky8vTwMGDFDKOXtZZWdnKzk5WZmZmec997x589SxY0fdcsst5x138OBBrV69ulp7WVmZXnvtNUlSenq6l58IAAAACNDwLUlPPPGE+vfvr2nTpum9997TqVOn9OWXX+rXv/612rRpo0WLFrmNnzVrljZv3qwZM2bUeE673a4lS5bovvvuU4sWF17u/sorr+jJJ5/UN998o9OnT6ugoECjRo1SSUmJJk+erJtvvrmuHxMAAADNSMCG7/DwcGVnZ+vhhx/W1KlT1apVK914441KTk7Wli1b1LVrV7fxY8aMUWRkZNVOJJ689tprOnv2rO65554LXn/nzp164okn9NFHH2nAgAFq2bKlhg0bJqvVqlWrVmnevHl1/owAAABoXiwul8vl7yKaC4fDIZvNJrvdrqioKH+XAwAAAJMF7Mw3AAAA0NQQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAk1z4STMAAHjgdEo5OVJJiRQTI6WkSFarv6sCgMBG+AYA+CwrS3rwQWn//p/aYmOluXOl1FT/1QUAgY5lJwAAn2RlSWlp7sFbkg4cMNqzsvxTFwA0Bjzh0kQ84RJAY+d0Sp07Vw/elSwWYwa8qIglKADgCTPfAACv5eTUHLwlyeWSiouNcQCA6gjfAACvlZTU7zgAaG4I3wAAr8XE1O84AGhuCN8AAK+lpBhrui0Wz/0WixQXZ4wDAFRH+AYAeM1qNbYTlKoH8Mqf58zhy5YAUBPCNwDAJ6mp0ooVUseO7u2xsUY7+3wDQM3YatBEbDUIoCnhCZcA4DuecAkAqBWrVRoyxN9VAEDjwrITAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCQt/F0AAKBxcjqlnByppESKiZFSUiSr1d9VAUBgI3wDAHyWlSU9+KC0f/9PbbGx0ty5Umqq/+oCgEDHshMAgE+ysqS0NPfgLUkHDhjtWVn+qQsAGgOLy+Vy+buI5sLhcMhms8lutysqKsrf5QCAz5xOqXPn6sG7ksVizIAXFbEEBQA8YeYbAOC1nJyag7ckuVxScbExDgBQHeEbAOC1kpL6HQcAzQ3hGwDgtZiY+h0HAM0N4RsA4LWUFGNNt8Xiud9ikeLijHEAgOoCOnyXlZXpqaeeUrdu3RQaGqoOHTpo/PjxOnDggE/nSU9Pl8ViqfE1Z86cGo9dsmSJrrnmGkVERCg6OlojRoxQXl5eHT8ZADROVquxnaBUPYBX/jxnDl+2BICaBGz4Lisr0/XXX69Zs2YpIyNDR44c0erVq5WXl6fevXtrz549Pp2vTZs26tatm8dXdHS0x2OmTJmiu+66SzfddJO+++475efnKzQ0VCkpKXrvvffq42MCQKOTmiqtWCF17OjeHhtrtLPPNwDULGC3GpwxY4YyMjL04osvauLEiVXt27ZtU8+ePTVw4EBt3LjRq3Olp6erc+fOevLJJ72+/tq1a3XDDTcoLS1N77zzTlX7qVOn1K1bNzkcDu3Zs0dt27b1+pxsNQigKeEJlwDgu4Cc+T516pQyMzMVHBysu+++260vKSlJ/fv316effup1+K6N559/XpJ0zz33uLWHhYVp7Nixstvtmj9/foNdHwACndUqDRki3XGH8U7wBoALC8jwvW7dOh0/flw9e/ZUy5Ytq/UPGjRIkrRq1aoGub7D4dCGDRtksVg0YMAA068PAACApikgw/fWrVslSfHx8R77K9sLCgq8Pmd+fr5GjhypmJgYhYSEKC4uTmPHjvV4jsLCQjmdTrVr185j+K+8fmFhoSoqKryuAQAAAM1bQIbvkv88naGmL0K2bt3abZw3cnJyNHr0aG3fvl1Hjx7VggUL9O9//1v9+vVzW9Pty/XLysp09OhRr2sAAABA8xaQ4fvkyZOSpODgYI/9ISEhbuMu5OGHH9amTZs0fvx4RUdHKywsTDfccIPeffddnTlzRuPGjdPBgwd9vv6FaigvL5fD4XB7AQAAoPkKyPAdHh4uSTp9+rTH/vLycrdxF9KzZ0916dLFY/ugQYN04sQJLV++3OfrX6iGjIwM2Wy2qldcXJxX9QIAAKBpCsjw3b59e0nS4cOHPfYfOXJEkhRTD88vTkhIkCTt2rXL5+uHhoaqVatWNZ57+vTpstvtVa/i4uI61wsAAIDGq4W/C/CkV69ekqSioiKP/ZXtPXr0qPO1PH1hMikpSUFBQTp06JBOnDihiIgIj9dPTExUUFDN//8lJCTEbYkKAAAAmreAnPkeOnSoIiIiVFBQoBMnTlTrz83NlSSNHDnygufKzc1V+/bta/xi5N69eyVJV1xxRVVbVFSUBg8eLJfL5fFR8r5cHwAAAKgUkOE7PDxckyZNUnl5uRYvXuzWt337duXl5WnAgAFKSUmpas/OzlZycrIyMzPdxp8+fVoHDx7Uhx9+WO06+fn5ys3NVXh4uG677Ta3vqlTp0qSXn31Vbf2srIyLV26VFFRUZowYUKdPicAAACal4AM35L0xBNPqH///po2bZree+89nTp1Sl9++aV+/etfq02bNlq0aJHb+FmzZmnz5s2aMWOGW7vFYpEkTZw4UUuWLNEPP/ygsrIyffTRR0pLS1OLFi20YMGCqnXelYYPH64HHnhAy5cv17PPPiuHw6Hi4mKNHTtW+/fv14IFC9SuXbsG/TsAAABA0xKw4Ts8PFzZ2dl6+OGHNXXqVLVq1Uo33nijkpOTtWXLFnXt2tVt/JgxYxQZGanx48e7tV933XX65JNPNGbMGP31r3/VlVdeKZvNpnvuuUeDBg3S5s2bdfvtt3usYe7cuXr99de1atUqxcTEKCkpSSdOnFBOTo7S0tIa7LMDAACgabK4XC6Xv4toLhwOh2w2m+x2u6KiovxdDgAAAEwWkLudAAACn9Mp5eRIJSVSTIyUkiJZrf6uCgACG+EbAOCzrCzpwQel/ft/aouNlebOlVJT/VcXAAS6gF3zDQAITFlZUlqae/CWpAMHjPasLP/UBQCNAWu+TcSabwCNndMpde5cPXhXsliMGfCiIpagAIAnzHwDALyWk1Nz8JYkl0sqLjbGAQCqI3wDALxWUlK/4wCguSF8AwC8FhNTv+MAoLkhfAMAvJaSYqzp/s/Dg6uxWKS4OGMcAKA6wjcAwGtWq7GdoFQ9gFf+PGcOX7YEgJoQvgEAPklNlVaskDp2dG+PjTXa2ecbAGrGVoMmYqtBAE0JT7gEAN/xhEsAQK1YrdKQIf6uAgAaF5adAAAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJmnh7wIAAI2T0ynl5EglJVJMjJSSIlmt/q4KAAIb4RsA4LOsLOnBB6X9+39qi42V5s6VUlP9VxcABDqWnQAAfJKVJaWluQdvSTpwwGjPyvJPXQDQGFhcLpfL30U0Fw6HQzabTXa7XVFRUf4uBwB85nRKnTtXD96VLBZjBryoiCUoAOAJM98AAK/l5NQcvCXJ5ZKKi41xAIDqCN8AAK+VlNTvOABobgjfAACvxcTU7zgAaG4I3wAAr6WkGGu6LRbP/RaLFBdnjAMAVEf4BgB4zWo1thOUqgfwyp/nzOHLlgBQE8I3AMAnqanSihVSx47u7bGxRjv7fANAzdhq0ERsNQigKeEJlwDgO55wCQCoFatVGjLE31UAQOPCshMAAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJAEdvsvKyvTUU0+pW7duCg0NVYcOHTR+/HgdOHDAp/OsW7dOv/vd79S1a1eFhoYqPDxc3bt319SpU1VaWurxmEWLFslisdT4Gj16dD18QgAAADQnAbvbSVlZma6//npt27ZNixYt0q9+9Svt2LFDY8eOVe/evZWbm6uEhIQLnuell17SpEmT1KNHD82fP1/9+vWTw+HQG2+8occff1xLly7VJ5984vFcYWFh6tSpk8fzxsbG1vkzAgAAoHkJ2PD99NNPKzc3Vy+++KJS//PEhr59+2r58uXq2bOn0tPTtXHjxgue59SpUwoODtaaNWsUFxcnSYqMjNTUqVN1+PBh/elPf9IDDzygDz74oNqxycnJWr9+fb1+LgAAADRfAbns5NSpU8rMzFRwcLDuvvtut76kpCT1799fn376qVfhu3379rrjjjuqgve5Ro4cKUlau3atnE5n/RQPAAAA1CAgw/e6det0/Phx9ezZUy1btqzWP2jQIEnSqlWrLniusWPHatGiRR77bDabJMlisYgHfQIAAKChBWT43rp1qyQpPj7eY39le0FBQZ2us2vXLklGmG/RovoKnNLSUk2YMEHx8fEKDQ1V27ZtNXz4cL3//vt1ui4AAACap4AM3yUlJZKk6Ohoj/2tW7d2G1dbixcvlsVi0cyZMz3279y5U23bttWGDRvkcDiUnZ0tSRo1apSmTJlywfOXl5fL4XC4vQAAANB8BWT4PnnypCQpODjYY39ISIjbuNr45z//qffff18PPfSQfvGLX1TrHzZsmPLy8vTss8+qU6dOCg4OVlJSklauXKkOHTpo9uzZWrNmzXmvkZGRIZvNVvXytO4cAAAAzUdAhu/w8HBJ0unTpz32l5eXu43z1VdffaW77rpLqamp+stf/uJxTGxsrJKTk6u1h4aGavz48ZKkhQsXnvc606dPl91ur3oVFxfXql4ACEROp7R+vbRsmfHO99YB4MICcqvB9u3bS5IOHz7ssf/IkSOSpJiYGJ/PXVRUpF/+8pcaOHCgli1bJqvV6vM5KvcEr1wzXpOQkJCqWXoAaEqysqQHH5T27/+pLTZWmjtX+s/usAAADwJy5rtXr16SjKDsSWV7jx49fDrv7t27dd111+naa6/VihUralzWciEVFRW1Og4AmoKsLCktzT14S9KBA0Z7VpZ/6gKAxiAgw/fQoUMVERGhgoICnThxolp/bm6upJ/26fbGjh07dN111+mGG27Q0qVL3XY3eeGFF6otCWnfvr3HB+9I0t69eyVJV1xxhdfXB4CmwOk0Zrw97c5a2fbQQyxBAYCaBGT4Dg8P16RJk1ReXq7Fixe79W3fvl15eXkaMGCAUlJSqtqzs7OVnJyszMzMaufLz8/X4MGDlZaWpr/97W8KCnL/2I8++mhVoK508OBBrV69utq5ysrK9Nprr0mS0tPTa/sRAaBRysmpPuN9LpdLKi42xgEAqgvINd+S9MQTT2jDhg2aNm2a2rdvr+HDh2vnzp367W9/qzZt2lR7cM6sWbO0efNm7dy5U5MnT65q//zzz3XDDTeovLxchw4d0h133OF1Da+88oouueQSjRs3TjExMdq1a5ceffRRlZSUaPLkybr55pvr6+MCQKPg7Q6vddwJFgCarIAN3+Hh4crOzlZGRoamTp2q22+/XdHR0Ro+fLieeeYZxcbGuo0fM2aMPvnkE40bN86tfenSpVVf0Hz77be9vv7OnTv19ttv68MPP9Qrr7yi0tJStWrVSldffbVWrVrl05IXAGgqvP2eey2+Dw8AzYLFxXPVTeNwOGSz2WS32xUVFeXvcgDAZ06n1Lmz8eVKT//1sFiMXU+KiqRabCYFAE1eQK75BgAEJqvV2E5QMoL2uSp/njOH4A0ANSF8AwB8kpoqrVghdezo3h4ba7SzzzcA1IxlJyZi2QmApsTpNHY1KSkx1ninpDDjDQAXErBfuAQABDarVRoyxN9VAEDjwrITAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCQt/F0AAKBxcjqlnByppESKiZFSUiSr1d9VAUBgI3wDAHyWlSU9+KC0f/9PbbGx0ty5Umqq/+oCgEDHshMAgE+ysqS0NPfgLUkHDhjtWVn+qQsAGgNmvgEA51V2slzf7tyvCmeFKpzStPukSJeHgS7JYpH++37pihgpqm1LdUxoL4vFYnrNABCoLC6Xy9MtFA3A4XDIZrPJbrcrKirK3+UAgEfOs07t2rRHX64t0Jf/KtCuz3br7Blnrc4V3b6Vev0iUX2G9VCfYT3ULrZNPVcLAI0L4dtEhG8AgerQ/h+16YMvtekfW5SfXaiTjlMNcp24Kzrq6ht66pob+yjpuqsUHHJRg1wHAAIV4dtEhG8AgcLlcmnPliJtzPpM//77F/q/rd/UOLbj5TFKHHSFwiJDdeCA9O67Fz5/aqpkKftehTk7dep4mccxYS1D1WdYkgaM7KcBI69WVHRkbT8OADQahG8TEb4B+Nvh74/o3b+u0Scr8vT9vkMex9jaRqr3sB7qO6yH+gxL0sWd2lX1OZ1S587Glys9/dfDYjF2PSkqMrYdPHvmrHZ9tltfrC3Qlx8XaNemPapwVlQ7ztrCqp5Du+v/pQ/VkDEDFRTEfgAAmibCt4kI3wD85eyZs3pv7gda+swKnTxWfUlJ16u76JoRfXTNjX10ed/Lzht+K3c7kdwDeOX3KlesqHm7weNHT+iLj7bqsw++1KYPvpS99Fi1MVcN6KqJ88ara98uXn8+AGgsCN8mInwD8Icd//5as++Zr32FxVVtQdYg9fpFoq695RoNHNVPbWJa+3ROT/t8x8VJc+Z4v8+30+nUzn/v1qfvbdLGrH+7zcQHBVk0atKvNO7Z2xXWMsyn2gAgkBG+TUT4BmCm0+Vn9PoTb+udF95X5a3eYrFoxO+H6e6nx6j1xbY6nb8+n3Dpcrn0+UdbNX/KIn2780BVe/v4izXt9UlKvPbKOtUKAIGC8G0iwjcAs3y393s9O+av2v1lUVXb5X3i9fArE3R5n8v8WNn5nT1zVu/O/rsWP/m2TpedkWTMgo994jbd+ViqrDy/HkAjR/g2EeEbgBk2f5ivP94+WyfsJyVJLS6y6q4nx+jXj46UtUXjCK8H9pToz+kvaUfuV1Vt/X7VW4+9+aAibBF+rAwA6obwbSLCN4CG9vdX1mre/a+qosK4tXe8PEaPv/WwEnrH+7ky3zmdTi177j0teWp51efp3D1Oz/3jMR7WA6DRInybiPANoCG988L7emXqkqqfB47qp6mvT1JEVLgfq6q7Leu26Zlf/1XHDh+XJF1yaTv9Zd0fFBN/iZ8rAwDfsZEqADQBq176p1vwTptys/7w7n81+uAtSb1/kaTMfz+nmMuMsH3wm0OaOuxplX532M+VAYDvCN8A0Mjlrf5cLz3wWtXP6U/frntfuKtJPaimY0KMZuc8o7grOkqSvi/6QTNvfl6nTnh+eiYABKqmc2cGgGao5P8O6vnfzqvaSvD2aaP1m8dv9XNVDaNNTGv9+eMndMmlxhM392wpUuakv/m5KgDwDeEbABqpiooK/Tn9RZ10GE+svO62ARr/3J1+rqphte0QrWfXTFdoRIgkae3rG5ST9ZmfqwIA7xG+AaCR+tfSHBVu3CVJirnsEk15dYIslc94b8I6d4/TA//z+6qfX35oocpPlfuxIgDwHuEbABoh51mnFj+1vOrnh+bf0yS+XOmtYWOvU79f9ZYkHdr/o/7+vx/7uSIA8A7hGwAaoZx3/63vi36QJPW9oaf6DOvh54rMZbFY9Pvnf1P184rZq+V0Ov1YEQB4J6DDd1lZmZ566il169ZNoaGh6tChg8aPH68DBw6Ydq4lS5bommuuUUREhKKjozVixAjl5eXV9iMBQL344G//qvrzmKmj/FiJ/8QnXarkEf+Z/S7+UZ9/uNXPFQHAhQVs+C4rK9P111+vWbNmKSMjQ0eOHNHq1auVl5en3r17a8+ePQ1+rilTpuiuu+7STTfdpO+++075+fkKDQ1VSkqK3nvvvfr6qADgkyM/2LU1u1CS1CGhvXoNTfRzRf5z4+9/WfXnDctz/VgJAHgnYMP3008/rdzcXGVkZCg1NVVhYWHq27evli9frtLSUqWnpzfoudauXavZs2crLS1NM2fOlM1mU6dOnfTGG2+oQ4cOGjdunEpLS+vvAwOAl75cW1D1uPXrbu3vty9ZOp3S+vXSsmXGuz9WfVw9vJfCWoZKkj7/MF88tBlAoAvI8H3q1CllZmYqODhYd999t1tfUlKS+vfvr08//VQbN25ssHM9//zzkqR77rnHrT0sLExjx46V3W7X/Pnza/PxAKBOduR9VfXn3tcn+aWGrCypc2dp6FDpzjuN986djXYzBYdcpMSUKyVJRw7a9f2+H8wtAAB8FJDhe926dTp+/Lh69uypli1bVusfNGiQJGnVqlUNci6Hw6ENGzbIYrFowIABdbo+ANS3fduLq/58ed/LTL9+VpaUlibt3+/efuCA0W52AO/a56e/g6Jt35p7cQDwUUCG761bjS/NxMfHe+yvbC8oKGiQcxUWFsrpdKpdu3YeA3vlMYWFhaqoqLhgDQBQnw4V/yhJimwdocjW1e9RDcnplB58UPK0uqOy7aGHzF2C0v6yS6r+XPl3AwCBKiDDd0lJiSQpOjraY3/r1q3dxtX3ubw9pqysTEePHq3x2uXl5XI4HG4vAKirE/aTkqSWJgdvScrJqT7jfS6XSyouNsaZJSr6p7+Hyr8bAAhUARm+T540bp7BwcEe+0NCQtzG1fe5vD3mQjVkZGTIZrNVveLi4i5YLwBcSOWXCv3xPUsv5jx8Glcvmv5DPQE0IQEZvsPDjae0nT592mN/eXm527j6Ppe3x1yohunTp8tut1e9iouLaxwLAN4KjwyT5J9Z3piY+h1XH879e6jc+QQAAlVAhu/27dtLkg4fPuyx/8iRI5KkGC/u7rU5l7fHhIaGqlWrVjVeOyQkRFFRUW4vAKirNh2NJXH20mM6dfyUqddOSZFiY2uedbdYpLg4Y5xZDu47VPXntrFtzLswANRCQIbvXr16SZKKioo89le29+hx4ccp1+ZcSUlJCgoK0qFDh3TixIkaj0lMTFRQUED+FQJowjpd0bHqz3vz95l6batVmjvX+PPPA3jlz3PmGOPMsif/p/t7pys7nmckAPhfQCbHoUOHKiIiQgUFBR7Db26u8RSzkSNHNsi5oqKiNHjwYLlcLo+Pkvfl+gBQ367s37Xqz/nZ202/fmqqtGKF1PFnOTc21mhPTTWvFqfTqW2f7JQktWwVodiuJq53AYBaCMjwHR4erkmTJqm8vFyLFy9269u+fbvy8vI0YMAApZzze83s7GwlJycrMzOzzueSpKlTp0qSXn31Vbf2srIyLV26VFFRUZowYUKdPysA+KrvL3/6Td3G9z7zSw2pqdK+fVJ2tvTmm8Z7UZG5wVuSCjbskOPHY5Kk3tcnymrmlDsA1IYrQJ04ccLVv39/V2RkpCsrK8t18uRJ1xdffOG66qqrXG3btnV99dVXbuNvvPFGlyRXy5Yt63yuSg888IBLkuuZZ55x2e1217fffuu69dZbXUFBQa533nnH589kt9tdklx2u93nYwHgXBOTp7mGWdJcwyxprq8+3+Pvcvzmud/Mqfp7WLdso7/LAYALCsiZb8mYsc7OztbDDz+sqVOnqlWrVrrxxhuVnJysLVu2qGvXrm7jx4wZo8jISI0fP77O56o0d+5cvf7661q1apViYmKUlJSkEydOKCcnR2lpaQ3yuQHAG8PHX1/153dmrfZjJf5z8JtD2rDcWBoYGd1Sg0b383NFAHBhFpfL03PK0BAcDodsNpvsdjs7nwCok7KT5fpt/P06esghi8Wi//niT0ro5flJvk3Vn8e9qLWvb5AkjZ2ZprufGuPnigDgwgJ25hsAULPQ8BCNmTZakvHQnb/+fr7Onjnr36JM9OXHBVXBu2WrCKU+dKOfKwIA7xC+AaCRGjlxeNXWeru/+D/9bdpSP1dkjh9Ljuj5386r+nncs3cosnXL8xwBAIGD8A0AjVRwyEWa+vpktbjI2OHj3Tl/199fWevnqhrWqeOnNHPk8zpy0C5J6ntDT9004Zd+rgoAvEf4BoBGrNvVXXT/3J++aD73vle17s0cP1bUcE6dKNPMkX/S7i/+T5J0yaXtNG3xZB52BqBR4Y4FAI3czRNuUNqUmyUZ67+f/22mVs//yM9V1S/H4WP67xue0db1xkOFIltH6Nk109X6YpufKwMA3xC+AaAJuOcvv9VN9xrLL1wul+bd/6rmT1kk51mnnyuru293HdDk/jO0I+9rSVKELVwZ/3xcnbvH+bkyAPAd4RsAmgCLxaIH/uf3+vV/jaxqe3fO3zX1l0/rx5IjfqysbtYt26hJyf+t7/Z8L0lqdbFNL2Q/qW79EvxcGQDUDvt8m4h9vgGYYc3/rtWLkxdUzXpHtYnUgy//XtelDfBzZd47fvSEXnxggf619Kf165f1uFRPr5qmSy5t58fKAKBuCN8mInwDMEvhp7v0x9tnq/TA4aq2lFuv0f1zx6tth2g/VnZ+LpdLG9/bpBcnL9Dhc2bsf3nXYE1+6f9TWESoH6sDgLojfJuI8A3ATPZSh+be94py3v2sqi08MkxjZ6Zp9AO/0kXBF/mxuuq+3XVA86cs0uZ/5le1RdjCNSnzd7r+NymyWCz+Kw4A6gnh20SEbwBmc7lcyn7rU7380EIdPeSoao9u30qDRifr2tRr1GPwVWpxUQu/1Hf4+yP6dOVmbcz6t/Kzt6vCWVHVd82NffTA//xeF8e19UttANAQCN8mInwD8BfHj8f02ow39cHf/qWf3/YjW0eo7//rpWtG9FG/4b1ka9tw9yeXy6U9W4r02d+/1Gd//0Jfbd5brZ52cW00YdbdSrm1P7PdAJocwreJCN8A/G3PliItfeYdbfpHvs6Un6nWHxRkUdd+Ceo7rIf6/LKHrux/eZ2Xp/xYckRfflxgvNYW6PD3Rz2Oa9+5nf7fuF8o7ZGbFRoeUqdrAkCgInybiPANIFCcOn5Kmz7Yoo3vfabN/8zXCftJj+NCI0KUdN1VShx0hcIjw7w+v8vl0nd7vlf++kLtKyyucVznxDgNHNlPKbf2V5denZnpBtDkEb5NRPgGEIjOnjmr7blfadPfv9Smf2zRvu01h+W6Co0IUc8h3XXNiD665sY+urgT2wYCaF4I3yYifANoDEq/O6wtH2/TFx9v1ZaPt9W4TMQbQUEWdb26i/r8ZxnLVQO6BtwuKwBgJsK3iQjfABobl8ul/V9/p735+9x2IvFGy9Yt1X1gV0XYIhqoOgBofAjfJrLb7WrVqpWKi4sJ3wAAAAEsMjKyQb6H4p+NXZupH3/8UZIUFxfn50oAAABwPj/88IPatav/76UQvk0UHW080vnbb7+VzWbzczVoyhwOh+Li4vgtCxoc/9ZgFv6twSyV/9aCg4Mb5PyEbxMFBQVJkmw2GzcOmCIqKop/azAF/9ZgFv6twSwNtfVpUIOcFQAAAEA1hG8AAADAJIRvE4WEhOgPf/iDQkJ4bDIaFv/WYBb+rcEs/FuDWRr63xpbDQIAAAAmYeYbAAAAMAnhGwAAADAJ4RsAAAAwCeEbaELKysr0hz/8QV27dlVoaKg6dOig8ePHa//+/f4uDQAAiC9cAk1GWVmZrr/+euXm5iomJkYpKSnat2+fNm3apHbt2ikvL09dunTxd5kAADRrzHwDTcRzzz2n3NxcDRgwQF9//bXefvttffbZZ5o1a5YOHTqk8ePH+7tEAACaPcK3CT755BONGjVKl156qSwWi5588kl/l4Qm5syZM8rMzJQkvfTSS2rZsmVV35QpU9SjRw998skn+uKLL/xVIhoR7lnwt02bNslisXj1uvnmm/1dLhqBFStWKDk5WdHR0QoNDVVCQoIef/xxnT592vRaWph+xWbo+PHjuuqqq3TnnXfqoYce8nc5aII2btyoo0ePqkuXLurdu3e1/rS0NBUUFGj16tXq27evHypEY8I9C/526tQp/eY3v9Ebb7yhli1b6tZbb3Xrf+utt3T69GndddddhG94JTo6WlOnTtWVV16piIgIbdmyRffee68cDofmzZtnai2EbxOMGDFCI0aMkCRNmzbNz9WgKdq6daskqU+fPh77K9srxwHnwz0L/jZ48GC1adNGb7zxhq6++motWrSoqs/hcGjx4sXq0qWLWztwPr/4xS/cfu7cubM++eQTffzxx6bXwrIToAn49ttvJUmxsbEe+yvbK8cBQKArKCiQJPXo0cOtPT8/Xy6XS0lJSf4oC03Ezp079Y9//ENDhw41/dqEb6AJOH78uCQpPDzcY39ERITbOAAIdNu2bZOkaiE7Pz9fUvVQDnijZcuWCgkJ0VVXXaUhQ4Zo9uzZptfAshOgCajcMdRisZy3HwAai/PNfHtqB7yRn5+vU6dO6fPPP9d///d/65JLLtFTTz1lag2Eb6AJiIyMlCSdOHHCY//JkyclyW0XFAAIZNu2bVNQUJASExPd2ivDN8tOUBsJCQmSjH8/QUFBGj9+vKZOnVr1G2IzsOwkALHFEnzVqVMnSarxSZaV7ZXjgPrC/QoN4ejRoyouLlaXLl3cltOdPXtWO3bsUHh4OA8NQ71wuVw6c+aMW1tD39eY+TbB8ePHtWfPHknS6dOn9f333ys/P1/BwcG66qqrqo1niyX4qmfPnpKkL7/80mN/ZTu/poU3fLlncb9CQ6hc7/3ze9auXbtUXl6ufv36KSiI+UN475lnntE111yjyy67TC6XS5s2bdK0adM0atQotWrVym1sg9/XXGhw2dnZLknVXpdeemmNx2zbts0lyTVkyBC3drvd7rJYLK6EhIQGrhqNSXl5uctms7kkub788stq/T169HBJcm3atMkP1aGx8fWexf0K9e3FF190SXI99dRTbu3vvPOOS5Jr7NixfqoMjdX06dNdXbt2dYWFhbkiIyNdSUlJrueff9518uRJj+Mb8r7G/200wZAhQ+Ryuaq99u3bV+MxbLEEXwQHB2vSpEmSpEmTJrmt/f7rX/+qgoICXXvtterXr5+/SkQj4us9i/sV6ltNM9+VOzZt3bpVf//73/3ydEI0Ts8995y++uornTx5Ug6HQwUFBZo2bZrCwsI8jm/I+xrLTryQlpamwsJCn45ZvHixkpOTa31NtliCrx5//HF9/PHHys3N1eWXX66UlBR98803+uyzz9SmTRstXLjQ3yWiieJ+hfpWU/i+/vrr1aNHD+3cuVPTpk2rehgUUN8a8r5G+PbCvn379NVXX/l0TOXuErXFFkvwVWhoqLKzs5WRkaE333xTK1euVOvWrXX33XfrmWeeUVxcnL9LRBPF/Qr17dNPP/XYHhcXx5N6YYqGvK9ZXC42AA5EnTp10oEDB3Ts2DG3b3r36dNHW7Zs0ddff63LL7/cjxUCgIH7FYCmpiHva4TvBlbTQ088qfyf4ujRo2rdurUuv/xyff3111X9Z8+eVcuWLWW1WnXs2DG+6Q2g3vl6z+J+BSDQBdp9jbthA/P0paWaXpUutMVS9+7d+Q8ZgAbh6z2L+xWAQBdo9zXuiAGopnVGu3btkiR169bN9JoAwBPuVwCamoa+rxG+TbBw4UINHTpU7dq1U2RkpPr27as33nijxvFssQTAn3y5Z3G/AtAYBNJ9jfBtgn/9618aOXKkPvjgA23ZskW33367fvvb3+rtt9/2OP5CWyzt2rVL06ZN00UXXdTgtQNofny5Z3G/AtAYBNJ9jS9c+smIESMUFhamd99919+lAMAFcc8C0NT4677GzLefHD16VG3btvV3GQDgFe5ZAJoaf93XCN9+8Prrr+vzzz/Xvffe6+9SAOCCuGcBaGr8eV8jfJts1apVuvfee/W///u/6tOnj7/LAYDz4p4FoKnx932N8G2it956S2PGjNH8+fM1btw4f5cDAOfFPQtAUxMI9zXCt0leffVVpaena9GiRUpPT/d3OQBwXtyzADQ1gXJfY7cTE/z1r3/V1KlT9dJLL2nUqFFV7VarVe3atfNjZQBQHfcsAE1NIN3XCN8m6Ny5s7755ptq7Zdeeqn27dtnfkEAcB7cswA0NYF0XyN8AwAAACZhzTcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAwGffffed7r//fk2cOFHl5eWaPHmyLr74Yl188cV66qmnxMOTAcAzwjcAwCdr165VYmKiXn75ZV155ZV64okn9OKLL6pbt246efKknnzySf35z3/2d5kAEJAsLqYnAABe2rt3r3r27Kny8nK9+eabuuWWW9SuXTtFRkZq3759euedd3T77berbdu2+uGHH2SxWPxdMgAEFGa+AQBemzZtmk6cOKGpU6fqtttu0xdffKGjR49q8ODBCgoK0m233aaWLVuqtLRU33//vb/LBYCAQ/gGAHjlhx9+0KpVqxQZGamZM2dKkrZs2SJJ6tWrlyQpKChIbdq0kSQ5nU6/1AkAgYzwDQDwyrp163T27Fn96le/UmhoqCQpPz9fktSzZ8+qcaWlpbJYLLrkkkv8USYABDTCNwDAK1u3bpXkHrQr2ypnvnfv3q0TJ06oe/fuuuiii0yvEQACHeEbAOCV0tJSSVKrVq0kSRUVFdq2bZs6duyotm3bSpL+8Y9/SJJGjBjhlxoBINARvgEAXomMjJQk7dmzR9JPs9yVM+GnTp3SnDlzFBwcrHvvvddvdQJAICN8AwC8cu2110qSsrKy5HA4qtZ79+rVS0ePHtWtt96qoqIizZw5U5dddpkfKwWAwMU+3wAArzidTg0YMECbN29Wp06d1Lp1a23dulV9+vTRN998ox9//FGTJk1SZmamv0sFgIBF+AYAeO3HH3/Uo48+qlWrVunw4cOSpIiICF199dWaMmWKRo4c6ecKASCwEb4BALVSuZ/3jz/+6OdKAKDxYM03AMBnRUVFOnz4sPr27evvUgCgUSF8AwB89sUXX0gS4RsAfET4BgD4rDJ89+nTx8+VAEDjwppvAAAAwCTMfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJvn/AXAghfysTaavAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = 0.11   # arc-length between two successive points (resolution of solution points)\n",
    "\n",
    "z0 = np.array([0.5*np.pi, 0]) # unloaded state\n",
    "zp = z0\n",
    "dz= np.array([0, ds])  # assumed step to the next point (determines direction)\n",
    "\n",
    "stable_phi = []\n",
    "stable_mu =[]\n",
    "nonstable_phi = []\n",
    "nonstable_mu = []\n",
    "stable_prev = True   # we known by intution that the initial state is stable\n",
    "\n",
    "for k in range(17):\n",
    "    ZS = root(sys, zp+dz, args=(zp,ds), jac=dsys, method='hybr') \n",
    "    zs =ZS.x   # extract value from solution structure\n",
    "    stable = d2Epot_dphi2(zs[0], zs[1]) > 0\n",
    "    if stable:\n",
    "        stable_phi.append( zs[0] )\n",
    "        stable_mu.append( zs[1] )\n",
    "    else:        \n",
    "        nonstable_phi.append( zs[0] )\n",
    "        nonstable_mu.append( zs[1] )\n",
    "        \n",
    "    if not stable==stable_prev:\n",
    "        branchpoint = 0.5*(zp+zs) # stability change appeared somewhere in between\n",
    "    \n",
    "    # move on\n",
    "    dz = zs - zp \n",
    "    zp = zs\n",
    "    stable_prev = stable\n",
    "  \n",
    "# plot points of mu-phi curve\n",
    "plt.plot(stable_phi, stable_mu, 'bo')\n",
    "plt.plot(nonstable_phi, nonstable_mu, 'bo', fillstyle='none')\n",
    "plt.xlabel('$\\\\varphi$')\n",
    "plt.ylabel('$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = ['$-\\\\frac{1}{2}\\pi$', '$0$', '$\\\\frac{1}{2}\\pi$', '$\\pi$', '$\\\\frac{3}{2}\\pi$']\n",
    "plt.xticks(tick_pos, labels)\n",
    "\n",
    "# plot distance constraint for initial point\n",
    "plt.plot(z0[0], z0[1], 'bo')   # initial point\n",
    "phi_range = np.linspace(-np.pi, np.pi)\n",
    "mu_range =  np.linspace(-1.1*ds, 1.1*ds)\n",
    "PHI, MU = np.meshgrid(phi_range, mu_range)\n",
    "ZERO = np.zeros(PHI.shape)\n",
    "DS = np.sqrt(distance_squared(PHI, MU, ZERO, ZERO)) # refer to initial point (delta_phi=delta_mu=0)\n",
    "plt.contour(PHI+phi0, MU, DS, [ds])  # distance constraint\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A filled circle denotes a stable equilibrium and an empty circle an instable equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ellipse around the initial point shows where the next point, in this case the first point, has to be located, as dictated by the distance constraint. We will find intersections of the constraint-ellipse with the $\\phi$-$\\mu$-curve for all directions, even when the $\\phi$-$\\mu$-curve goes horizontal and the formulation for fixed $\\mu$ fails. For a quantitative analysis and to receive some inspiration for setting the pseudo arc-length metric it is advisable to calculate the determinant of the Jacobian (in function ``dsys``)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that the upright position is an equilibrium and at $\\mu=1$ it becomes unstable.\n",
    "Now we take a closer look at this point (stability loss) and try to follow the stable equilibrium branching off into the buckled state. Therefore we start the iterations from a state without load increase but rotated out of the vertical position ($\\varphi\\ne\\pi/2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAImCAYAAACLl+YVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABIaUlEQVR4nO3de3RU1d3/8c9kMAkBEq7KwASTkoJWwkU0gnQKqGWhVYgx1stDNeBaXgoI4k94QLGi1tinUoHo0mqRCChF0xGU2qqtAaJJK1ZugaDIChrSFI1IhlsCTub3x+kExkxgJiRnbu/XWrOm2XvPOd9pWaef7Oyzj8Xj8XgEAAAAoN3FhboAAAAAIFYQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL5N5PF45HK5xHONAAAAYhPh20SHDh1SSkqKDh06FOpSAAAAEAKEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCQdQl0AACAyud1SSYlUUyPZbJLDIVmtoa4KAMIbM98AgKA5nVJGhjR2rHTrrcZ7RobRDgBoGeEbABAUp1PKzZUyM6WyMunQIeM9M9NoJ4ADQMssHo/HE+oiYoXL5VJKSorq6uqUnJwc6nIAIGhutzHDnZkprVkjxZ0yhdPYKGVnS+Xl0u7dLEEBAH+Y+QYABKykRNq7V5o3zzd4S8bPc+dKlZXGOCASuN3S+vXSqlXGu9sd6ooQ7bjhEgAQsJoa433QIP/93nbvOCCcOZ3SjBnSvn0n2+x2afFiKScndHUhujHzDQAImM1mvJeX++/3tnvHAeHKe+/CqcFbkqqruXcB7Ys13yZizTeASMeab0QDt1tKS2sevL0sFmMGvLKSf8doe8x8AwACZrVKCxdK69YZQfvU3U6ys432p54isCC8lZS0HLwlyeORqqq4dwHtgzXfAICg5ORIRUXS/fdLl19+sj093WhnrSzCXaD3JHDvAtoD4RsAELScHGniRJ5wicgU6D0J3LuA9sCabxOx5hsAgNDzrvmurjaWmHwfa77RnljzDQAAYorVamwnKBlB+1TenxctInijfRC+AQBAzPHeu9C3r2+73c69C2hfLDsxEctOAAAIL2439y7AXNxwCQAAYpbVKo0ZE+oqEEtYdgIAAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYJOzDd319vR588EHFx8crLy8v6M+vX79eFovljK8x33u2bGFh4WnHZ2dnt8n3A4BI5XZL69dLq1YZ7253qCsCgPDXIdQFnE5xcbHuuusu1dbW6sSJE2d1rIEDB/ptb2xs1O7du3XhhRc26+vYsaP69evn93N2u/2s6gGASOZ0SvffL+3de7ItLU1auFDKyQlVVQAQ/sI2fK9atUr33HOPfv3rXyspKUlTpkw5q+Pt2rXLb/uaNWt0/fXXa+rUqc36srKytH79+rM6LwBEG6dTys2Vrr3WmPUeNEgqL5eeeMJoLyoigANAS8I2fKenp2vnzp3q06ePCgsLW32cXr166YYbbmixf8mSJRo7dqwGDRrU6nMAQKxwu40Z72uvldaskeL+u3hxxAjj5+xs6f/9P2niRMlqDWGhwH+53VJJiVRTI9lsksPBv02EVtiG7xEjRrTJcS666CIVFRX57duxY4eKi4vldDrb5FwAEO1KSoylJqtWnQzeXnFx0ty50uWXG+O+dysNYDqnU5oxQ9q372Sb3S4tXsxfZxA6YX/DZXsqKCjQ+eefrwkTJvjtr62t1d1336309HQlJiaqZ8+eGj9+vN58802TKwWA8FBTY7y39MdCb7t3HBAq3uVRpwZvSaquNtqZd0OoxGz4PnjwoFauXKlf/vKXsrbw96eKigr17NlTGzZskMvlUnFxsSRp4sSJmjVr1hnP0dDQIJfL5fMCgEhmsxnv5eX++73t3nFAKLjdxoy3x9O8z9s2cyY79CA0YjZ8L126VI2Njbrjjjv89l911VUqKyvT448/rn79+ik+Pl6ZmZlas2aN+vTpo6efflrr1q077Tny8/OVkpLS9EpNTW2PrwIApnE4jF1NnnhCamz07WtslPLzpfR0YxwQKiUlzWe8T+XxSFVVxjjAbDEZvhsbG/Xss8/q1ltvVY8ePfyOsdvtysrKataemJjYtPPKsmXLTnueuXPnqq6urulVVVV19sUDQAhZrcZ2guvWGTdXlpVJhw4Z79nZRvtTT3FDG0Ir0GVPLI9CKITtDZftad26daqsrNS0adNa9fmMjAxJLW9f6JWQkKCEhIRWnQMAwlVOjrGd4P33GzdXeqWns80gwkOgy55YHoVQiMnwvWTJEjkcDg0dOrRVn2/8/t9aASDG5OQY2wmyhRvCkcNh7GpSXe1/3bfFYvSzPAqhEHPLTioqKvT3v/9d06dPP+243r176+233/bbt2fPHknSBRdc0Ob1AUCksFqN7QRvucV4J3gjXFitxnaCkhG0T+X9edEi/s0iNKImfBcXFysrK0sFBQWnHbdkyRL17dtX119//WnH7d+/X2+99Vaz9vr6er300kuSpLy8vFbXCwAA2o93eVTfvr7tdjvLoxBaUbPsZOHChdq0aZMqKipanNWuq6vTihUrNHfuXHXocOav/sILL+i8887T5MmTZbPZtGvXLj3wwAOqqanR9OnTdd1117X11wAAAG2E5VEIR2E9822xWGSxWDR58mRJ0ssvv9zU9v1Hzt90003q0qVL004k/rz00kv67rvvdOedd57x3BUVFXr44Yf17rvvauTIkercubOuuuoqWa1WrV27VkuWLDmr7wYAANofy6MQbiwej79bEdAeXC6XUlJSVFdXp+Tk5FCXAwAAAJOF9cw3AAAAEE0I3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEk6hLoAAEBkcrulkhKppkay2SSHQ7JaQ10VAIQ3Zr4BAEFzOqWMDGnsWOnWW433jAyjHQDQMsI3ACAoTqeUmytlZkplZdKhQ8Z7ZqbRTgCHWdxuaf16adUq493tDnVFwJlZPB6PJ9RFxAqXy6WUlBTV1dUpOTk51OUAQNDcbmOGOzNTWrNGijtlCqexUcrOlsrLpd27WYKC9uV0SjNmSPv2nWyz26XFi6WcnNDVBZwJM98AgICVlEh790rz5vkGb8n4ee5cqbLSGAe0F+9fX04N3pJUXc1fXxD+CN8AgIDV1Bjvgwb57/e2e8cBbc3tNma8/f3d3ts2cyZLUBC+CN8AgIDZbMZ7ebn/fm+7dxzQ1kpKms94n8rjkaqq+OsLwhfhGwAQMIdDSkuTnnjCWON9qsZGKT9fSk83xgHtIdC/qvDXF4QrwjcAIGBWq7RwobRunXFz5am7nWRnG+1PPcXNlmg/gf5Vhb++IFyx24mJ2O0EQLRwOqX77zduvvRKTzeCNztNoD253cZfX6qr/a/7tliMXU8qK/klEOGJ8G0iwjeAaMITLhEq3t1OJN8AbrEY70VF/BKI8EX4NhHhGwCAtuFvn+/UVGnRIoI3whvh20SEbwAA2g5/fUEk6hDqAgAAAFrDapXGjAl1FUBw2O0EAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMEmHUBcAAIhMbrdUUiLV1Eg2m+RwSFZrqKsCgPBG+AYABM3plO6/X9q792RbWpq0cKGUkxOqqhAt+MUO0YxlJwCAoDidUm6ulJkplZVJhw4Z75mZRrvTGeoKEcmcTuMXubFjpVtvNd7T0vh3hehh8Xg8nlAXEStcLpdSUlJUV1en5OTkUJcDAEFzu6WMDCNor1kjxZ0yhdPYKGVnS+Xl0u7dzFQieN5f7L6fTCwW472oiL+sIPIx8w0ACFhJibHUZN483+AtGT/PnStVVhrjgGC43dKMGc2Dt3SybeZMYxwQyQjfAICA1dQY74MG+e/3tnvHAYEqKZH27Wu53+ORqqr4xQ6Rj/ANAAiYzWa8l5f77/e2e8cBgQr0FzZ+sUOkC/vwXV9frwcffFDx8fHKy8tr1THy8vJksVhafC1atKjFz65YsUKXXXaZOnXqpO7du+uaa65RWVlZ674MAEQ4h8O4+e2JJ4w13qdqbJTy86X0dGMcEIxAf2HjFztEurAO38XFxRo8eLCee+45nThx4qyO1aNHDw0cONDvq3v37n4/M2vWLN1222269tpr9e9//1tbtmxRYmKiHA6H3njjjbOqBwAikdVqbCe4bp1xc+Wpu51kZxvtTz3FzZYInsMh2e0nb678PotFSk3lFztEvrDd7WTVqlW655579Otf/1pJSUmaMmWKbr/9dhUWFgZ9rLy8PKWlpemRRx4J+DPvvfeexo0bp9zcXL3++utN7ceOHdPAgQPlcrn0+eefq2fPngEfk91OAEQLf/t8p6cbwZvdKNBa3t1OJN8bL9ntBNEkbGe+09PTtXPnTk2dOlWWln4NbkdPPvmkJOnOO+/0ae/YsaMmTZqkuro6Pf/886bXBQDhICdH+vxzqbhYevVV4333boIRzk5OjhGw+/b1bbfbCd6IHmH7hMsRI0aE7Nwul0sbNmyQxWLRyJEjm/WPGjVKkrR27Vo99NBDZpcHAGHBapXGjAl1FYg2OTnSxIk84RLRK2xnvtvali1bNGHCBNlsNiUkJCg1NVWTJk3Stm3bmo0tLy+X2+1Wr1691Llz52b96enpTeMav3/HEQAAOCveX+xuucV4J3gjmsRM+C4pKVF2drZ27NihgwcPaunSpfrHP/6hSy+91GdNtyTV/Hcfo5ZuxOzWrZskYyeWgwcPtnjOhoYGuVwunxcAAABiV0yE7/vuu08fffSRpkyZou7du6tjx44aN26c/vSnP+nEiROaPHmy9u/f3zT+6NGjkqT4+Hi/x0tISGg21p/8/HylpKQ0vVJTU9voGwEAACASxUT4HjJkiPr37++3fdSoUTpy5Ihee+21pvakpCRJ0vHjx/0er6GhodlYf+bOnau6urqmV1VVVWu/AgAAAKJATITv08nIyJAk7dq1q6mtd+/ekqQDBw74/cy3334rSUpMTFTXrl1bPHZCQoKSk5N9XgAAAIhdMR++/d0wmZmZqbi4OH399dc6cuRIs/7KykpJ0qBBgxQXF/P/FQIAACBAUZ8cS0tL1bt37xZvjNyzZ48k6YILLmhqS05O1ujRo+XxePw+Sr60tFSSNGHChLYvGAAAAFErasJ3cXGxsrKyVFBQ4NN+/Phx7d+/X++8806zz2zZskWlpaVKSkrSjTfe6NM3e/ZsSdKLL77o015fX6+VK1cqOTlZd999dxt/CwAAAESzqAnfCxcu1KZNmzRv3jyfdu/TMadOnaoVK1boq6++Un19vd59913l5uaqQ4cOWrp0adM6b6/x48fr3nvv1WuvvabHH39cLpdLVVVVmjRpkvbt26elS5eqV69epn0/AAAARD6Lx+PxhLqIlpzusfLLli1TXl5e088rVqzQ1KlTNXnyZC1evLip3ePx6IMPPtAf//hHlZaW6ssvv9Thw4dls9k0evRozZo1S0OGDGnxPMuXL1dBQYF27typc845RyNHjtT8+fN1+eWXB/19XC6XUlJSVFdXx82XAAAAMSisw3e0IXwDAADEtqhZdgIAAACEO8I3AAAAYJIOoS4AAABEN7dbKimRamokm01yOCSrNdRVAaFB+AYAtAqBCoFwOqUZM6R9+0622e3S4sVSTk7o6gJChWUnAICgOZ1SRoY0dqx0663Ge0aG0Q54OZ1Sbq5v8Jak6mqjnX8viEWEbwBAULyBKjNTKiuTDh0y3jMzCVQ4ye02Zrz97anmbZs50xgHxBK2GjQRWw0CiHRutzHDnZkprVkjxZ0yhdPYKGVnS+Xl0u7dLEGJdevXG38ROZPiYmnMmPauBggfzHwDAAJWUiLt3SvNm+cbvCXj57lzpcpKYxxiW01N244DogXhGwAQMG9QGjTIf7+3nUAFm61txwHRgvANAAiYNyiVl/vv97YTqOBwGLuaWCz++y0WKTXVGAfEEsI3ACBgDoeUliY98YSxxvtUjY1Sfr6Unk6ggrHmf/Fi4z9/P4B7f160iHsDEHsI3wCAgFmt0sKF0rp1xs2Vp+52kp1ttD/1FIEKhpwcqahI6tvXt91uN9rZ5xuxiN1OTMRuJwCihdMp3X+/cfOlV3q6EbwJVPg+HsgEnET4NhHhG0A0IVABQPB4vDwAoFWsVvZnBoBgseYbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADBJh1AXAAAAIoPbLZWUSDU1ks0mORyS1RrqqoDIQvgGAABn5HRKM2ZI+/adbLPbpcWLpZyc0NUFRBqWnQAAgNNyOqXcXN/gLUnV1Ua70xmauoBIZPF4PJ5QFxErXC6XUlJSVFdXp+Tk5FCXAwBnhSUIscHtltLSmgdvL4vFmAGvrOR/fyAQzHwDAILmdEoZGdLYsdKttxrvGRnMgEajkpKWg7ckeTxSVZUxDsCZEb4BAEHxLkHIzJTKyqRDh4z3zEyWIESjmpq2HQfEOpadmIhlJwAindttzHBnZkpr1khxp0zhNDZK2dlSebm0ezdLEKLF+vXGXzbOpLhYGjOmvasBIh8z3wCAgJWUSHv3SvPm+QZvyfh57lxj7S9LEKKHw2Gs6bZY/PdbLFJqqjEOwJkRvgEAAfMuLRg0yH+/t50lCNHDajW2E5SaB3Dvz4sW8ZcOIFCEbwBAwGw247283H+/t907DtEhJ0cqKpL69vVtt9uNdvb5BgLHmm8TseYbQKRjzXdsY3tJ4OzxhEsAQMCsVmnhQmNXk+xsY433oEFG4M7Pl9atM2ZCCWTRyWrlpkrgbBG+AQBB8S5BuP9+6fLLT7anp7MEAQDOhGUnJmLZCYBowhIEAAgeM98AgFZhCQIABI/dTgAAAACTEL4BAAAAk4R9+K6vr9eDDz6o+Ph45eXlteoY77//vu644w4NGDBAiYmJSkpK0kUXXaTZs2ertrbW72cKCwtlsVhafGVnZ7f+SwEAACAmhfWa7+LiYt11112qra3ViRMnWnWMZ599VtOmTdPgwYP1/PPP69JLL5XL5dIrr7yihx56SCtXrtTGjRuVkZHR7LMdO3ZUv379/B7Xbre3qh4AAADErrAN36tWrdI999yjX//610pKStKUKVNadZxjx44pPj5e69atU2pqqiSpS5cumj17tg4cOKDf/OY3uvfee/X22283+2xWVpbWr19/Nl8DAAAAaBK2y07S09O1c+dOTZ06VRaLpdXH6d27t2655Zam4H2qCRMmSJLee+89ud3uVp8DAAAACETYznyPGDGiTY4zadIkTZo0yW9fSkqKJMlisYjtzgEAANDewnbm2wy7du2SJI0aNUodOjT/PaS2tlZ333230tPTlZiYqJ49e2r8+PF68803zS4VAAAAUSCmw/fy5ctlsVg0f/58v/0VFRXq2bOnNmzYIJfLpeLiYknSxIkTNWvWLDNLBQDgrLnd0vr10qpVxjsrLgHzhe2yk/b217/+VW+++abuu+8+XXHFFc36r7rqKpWVlSkrK6upLTMzU2vWrFH//v319NNP64orrtC1117b4jkaGhrU0NDQ9LPL5WrbLwEAQICcTmnGDGnfvpNtdru0eLGUkxO6uoBYE5Mz359++qluu+025eTk6Le//a3fMXa73Sd4eyUmJjbtvLJs2bLTnic/P18pKSlNL383fQIA0N6cTik31zd4S1J1tdHudIamLiAWxVz4rqys1E9/+lNdfvnlWrVqlaxWa9DH8O4J7l0z3pK5c+eqrq6u6VVVVdWqmgEAaC2325jx9revgLdt5kyWoABmianwvXv3bv3kJz/Rj3/8YxUVFSk+Pr5Vx2lsbAxoXEJCgpKTk31eAACYqaSk+Yz3qTweqarKGAeg/cVM+N65c6d+8pOfaNy4cVq5cqXP7iZPPfVUs1np3r17+33wjiTt2bNHknTBBRe0X8EAALSBmpq2HQfg7ERN+C4uLlZWVpYKCgqa9W3ZskWjR49Wbm6u/vCHPyguzvdrP/DAA02B2mv//v166623mh2rvr5eL730kiQpLy+v7b4AAEQYds6IDDZb244DcHaiZreThQsXatOmTaqoqND06dOb2j/++GONGzdODQ0N+vrrr3XLLbcEfMwXXnhB5513niZPniybzaZdu3bpgQceUE1NjaZPn67rrruuPb4KAIQ9p1O6/35p796TbWlp0sKF7JwRbhwOY1eT6mr/674tFqPf4TC/NiAWhXX4/v5j5V9++WW9/PLLkoydRk6deb7pppu0ceNGTZ482eczK1eu1LfffitJWr16dcDnrqio0OrVq/XOO+/ohRdeUG1trbp27apLLrlEa9eubXo0PQDEGu/OGddea8x6DxoklZdLTzxhtBcVEcDDidVqbCeYm2sE7VMDuPf/ZhctMsYBaH8WD89VN43L5VJKSorq6uq4+RJARHK7pYwMKTNTWrNGOnUVX2OjlJ1tBPHduwlz4cbfPt+pqUbw5pclwDyEbxMRvgFEuvXrpbFjpbIyacSI5v1lZdLll0vFxdKYMWZXhzNxu41dTWpqjDXeDge/JAFmC+tlJwCA8OLdEWPQIP/93nZ2zghPViu/FAGhFjW7nQAA2p93R4zycv/93nZ2zgAA/1h2YiKWnQCIdKz5BoCzw8w3ACBgVquxneC6dUbQLiuTDh0y3rOzjfanniJ4A0BLWPMNAAhKTo6xneD99xs3V3qlp7PNIACcCctOTMSyEwDRhJ0zACB4zHwDAFqFnTMAIHis+QYAAABMQvgGAAAATEL4BgAAAEzCmm8AAMIQN7QC0YnwDQBAmHE6pRkzpH37TrbZ7dLixWzlCEQ6lp0AABBGnE4pN9c3eEtSdbXR7nSGpi4AbYN9vk3EPt8AgNNxu6W0tObB28tiMWbAKytZggJEKma+AQAIEyUlLQdvSfJ4pKoqYxyAyET4BgAgTNTUtO04AOGH8A0AQJiw2dp2HIDwQ/gGACBMOBzGmm6LxX+/xSKlphrjAEQmwjcAAGHCajW2E5SaB3Dvz4sWcbMlEMkI3wAAhJGcHKmoSOrb17fdbjfa2ecbiGxsNWgithoEAASKJ1wC0YknXAIAEIasVmnMmFBXAaCtEb4BAK3CzCwABI813wCAoDmdUkaGNHasdOutxntGBo8+B4AzIXwDAILidEq5uVJmplRWJh06ZLxnZhrtBHAAaBk3XJqIGy4BRDq325jhzsyU1qyR4k6ZwmlslLKzpfJyafdulqAAgD/MfAMAAlZSIu3dK82b5xu8JePnuXOlykpjHACgOcI3ACBgNTXG+6BB/vu97d5xAABfhG8AQMBsNuO9vNx/v7fdOy5aud3S+vXSqlXGu9sd6ooARArWfJuINd8AIh1rvo0bSmfMkPbtO9lmtxuPhefpkwDOhJlvAEDArFZp4UJp3TojaJ+620l2ttH+1FPRHbxzc32DtyRVV7PTC4DAMPNtIma+AUQLp1O6/37j5kuv9HQjeEfr7K/bLaWlNQ/eXhaLMQNeWRm9v3wAOHuEbxMRvgFEk1h7wuX69cbDhM6kuJjHwgNoGY+XBwC0itUaWyEz0B1c2OkFwOmw5hsAgAAEuoNLtO/0AuDsEL4BAAiAw2Gs6bZY/PdbLFJqqjEOAFpC+AYAIABWq7GdoNQ8gHt/XrQoute9Azh7hG8AAAKUkyMVFUl9+/q22+1Ge7Tu9AKg7bDbiYnY7QQAokOs7fQCoO2w2wkAAEGKtZ1eALQdlp0AAAAAJmHmGwAQsVj+ASDSEL4BABHJ6ZRmzPB93LvdbuxIwo2PAMJV2C87qa+v14MPPqj4+Hjl5eWd1XEWLFiggQMHKjExUX369NGUKVNUXV192s+tWLFCl112mTp16qTu3bvrmmuuUVlZWavrAACcPadTys31Dd6SVF1ttDudoakLAM4krMN3cXGxBg8erOeee04nTpxo9XHq6+t15ZVXauHChcrPz9e3336rt956S2VlZRo2bJg+//xzv5+bNWuWbrvtNl177bX697//rS1btigxMVEOh0NvvPFGq+sBALSe223MePvbq8vbNnOmMQ4Awk3Yhu9Vq1bp+uuv14wZM7Rw4cKzOtajjz6q0tJS5efnKycnRx07dtTw4cP12muvqba21u+M+nvvvaenn35aubm5mj9/vlJSUtSvXz+98sor6tOnjyZPnqza2tqzqgsAELySkuYz3qfyeKSqKmMcAISbsA3f6enp2rlzp6ZOnSpLS8/yDcCxY8dUUFCg+Ph43X777T59mZmZGjFihD788EN98MEHPn1PPvmkJOnOO+/0ae/YsaMmTZqkuro6Pf/8862uCwDQOjU1bTsOAMwUVPguKyvTkSNH2qsWHyNGjFCfPn3O+jjvv/++Dh8+rCFDhqhz587N+keNGiVJWrt2bVOby+XShg0bZLFYNHLkyIA+AwAwh83WtuMAwExB7XYyatQoxcXFqX///ho6dKjPy3bKVe7OO+/U8OHDddddd7V5wcHaunWrJGMm3R9v+7Zt25raysvL5Xa7de655/oN7N7PlJeXq7GxUXFxYfsHBABoN6Ha5s/hMHY1qa72v+7bYjH6HY72rwUAghVU+L7lllu0detWffrpp9q9e7def/31piUhPXv21NChQ9W/f3+98cYbeuutt8IifNf89++O3bt399vfrVs3n3HBfKa+vl4HDx5scVxDQ4MaGhqafna5XEFWDwDhKZTb/Fmtxnlyc42gfWoA965SXLSI/b4BhKegpmxfeeUVlZeX6/Dhw/roo4/0wgsv6J577tHIkSN18OBBvffee/r973+vb775RvHx8e1Vc1COHj0qSS3Wk5CQ4DMumM98/3Pfl5+fr5SUlKZXampqcMUDQBgKh23+cnKkoiKpb1/fdrvdaGefbwDhqlUP2UlISNAll1yiSy65pKmtrq5Oy5Yt04MPPqjMzEwtX768zYo8G0lJSZKk48eP++33zkx7xwXzme9/7vvmzp2rWbNmNf3scrkI4AAi2pm2+bNYjG3+Jk5s/5nnnBzjPDzhEkAkabMnXKakpGjmzJnq37+/srOztX37dg0YMKCtDt9qvXv3liQdOHDAb/+3334rST5r1gP9TGJiorp27driuRMSEnxmyQEg0gWzzd+YMe1fj9VqznkAoK20+Z2C1113nQYOHKgnnniirQ/dKkOHDpUkVVZW+u33tg8ePLipLTMzU3Fxcfr666/97u7i/cygQYO42RJATGGbPwA4O+2SHH/wgx+ooqKiPQ4dtLFjx6pTp07atm2b3yBdWloqSZowYUJTW3JyskaPHi2Px+P3UfL+PgMAsYBt/gDg7AQVvmfMmKHCwkJt3bpV3333XYvj9uzZ07R0wyzFxcXKyspSQUGBT3tSUpKmTZumhoaGZuvQd+zYobKyMo0cOVKO7+1JNXv2bEnSiy++6NNeX1+vlStXKjk5WXfffXc7fBMACF/ebf5aevaZxSKlprLNHwC0JKg13wUFBU1bC55zzjm66KKLNGzYMA0bNkyDBg1SQkKCVq5cqc8++0zPPfdcuxTckoULF2rTpk2qqKjQ9OnTffoefvhhbdiwQXPmzFHv3r01fvx4VVRU6Be/+IV69OihwsLCZscbP3687r33Xi1ZskSZmZm69957VVdXp/vuu0/79u3T6tWr1atXL5O+HQCEB7b5A4CzY/F4/N2z7t+aNWu0efPmpld1dbVxkO9NgWRkZCgvL0/Dhg3TxRdfrHPPPbd1xZ3msfLLli1TXl5e088rVqzQ1KlTNXnyZC1evLjZ+Pr6euXn5+vVV1/Vl19+qe7du2v8+PF67LHHZLfbWzzP8uXLVVBQoJ07d+qcc87RyJEjNX/+fF1++eVBfx+Xy6WUlBTV1dUpOTk56M8DQLjwt893aqoRvNnmDwBaFlT4/r7a2lqfML5582Z9/vnnamxs9AnONptNw4cPj/nHsRO+AUSTUD3hEgAi2VmFb3+OHDmirVu3+gTyHTt26MSJE3K73W15qohD+AYAAIhtbR6+/fnuu++0Y8cODRkypL1PFdYI3wAAALHNlE2qO3ToEPPBGwAAAOAJMQAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEk6hLoAAEBkcrulkhKppkay2SSHQ7JaQ10VAIQ3wjcAIGhOpzRjhrRv38k2u11avFjKyQldXQAQ7lh2AgAIitMp5eb6Bm9Jqq422p3O0NQFAJHA4vF4PKEuIla4XC6lpKSorq5OycnJoS4HAILmdktpac2Dt5fFYsyAV1ayBAUA/GHmGwAQsJKSloO3JHk8UlWVMQ4A0BzhGwAQsJqath0HALGG8A0ACJjN1rbjACDWEL4BAAFzOIw13RaL/36LRUpNNcYBAJojfAMAAma1GtsJSs0DuPfnRYu42RIAWkL4BgAEJSdHKiqS+vb1bbfbjXb2+QaAlrHVoInYahBANOEJlwAQPJ5wCQBoFatVGjMm1FUAQGRh2QkAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYJIOoS4AABCZ3G6ppESqqZFsNsnhkKzWUFcFAOGN8A0ACJrTKc2YIe3bd7LNbpcWL5ZyckJXFwCEO5adAACC4nRKubm+wVuSqquNdqczNHUBQCQI6/BdX1+vBQsWaODAgUpMTFSfPn00ZcoUVVdXB3yM9evXy2KxnPE1ZswYn88VFhaednx2dnbbflkAiAButzHj7fE07/O2zZxpjAMANBe2y07q6+t15ZVXavv27SosLNTVV1+tnTt3atKkSRo2bJhKS0uVkZER8PEGDhzot72xsVG7d+/WhRde2KyvY8eO6tevn9/P2e32gM8NANGipKT5jPepPB6pqsoY9705DQCAwjh8P/rooyotLdUzzzyjnP8uIBw+fLhee+01DRkyRHl5efrggw8CPt6uXbv8tq9Zs0bXX3+9pk6d2qwvKytL69evb1X9ABCNamradhwAxJqwXHZy7NgxFRQUKD4+XrfffrtPX2ZmpkaMGKEPP/wwoPDdq1cv3XDDDS32L1myRGPHjtWgQYPOum4AiHY2W9uOA4BYE5bh+/3339fhw4c1ZMgQde7cuVn/qFGjJElr164947EuuugiFRUV+e3bsWOHiouLNX369LMrGABihMNh7Gpisfjvt1ik1FRjHACgubAM31u3bpUkpaen++33tm/btu2szlNQUKDzzz9fEyZM8NtfW1uru+++W+np6UpMTFTPnj01fvx4vfnmm2d1XgCIVFarsZ2g1DyAe39etIj9vgGgJWEZvmv+u1iwe/fufvu7devmM641Dh48qJUrV+qXv/ylrC38v0RFRYV69uypDRs2yOVyqbi4WJI0ceJEzZo164znaGhokMvl8nkBQKTLyZGKiqS+fX3b7XajnX2+AaBlYXnD5dGjRyVJ8fHxfvsTEhJ8xrXG0qVL1djYqDvuuMNv/1VXXaWysjJlZWU1tWVmZmrNmjXq37+/nn76aV1xxRW69tprWzxHfn6+FixY0OoaASBc5eRIEyfyhEsACFZYznwnJSVJko4fP+63v6GhwWdcsBobG/Xss8/q1ltvVY8ePfyOsdvtPsHbKzExUVOmTJEkLVu27LTnmTt3rurq6ppeVVVVraoXAMKR1WpsJ3jLLcY7wRsAziwsZ7579+4tSTpw4IDf/m+//VaSZGvl7fTr1q1TZWWlpk2b1qrPe/cXb2n7Qq+EhISmWXoAAAAgLGe+hw4dKkmqrKz02+9tHzx4cKuOv2TJEjkcjqbzBKuxsbFVnwMAAEBsC8vwPXbsWHXq1Enbtm3TkSNHmvWXlpZKUou7lJxORUWF/v73v59xe8HevXvr7bff9tu3Z88eSdIFF1wQ9PkBAAAQu8IyfCclJWnatGlqaGjQ8uXLffp27NihsrIyjRw5Uo5TNpItLi5WVlaWCgoKTnvsJUuWqG/fvrr++utPO27//v166623mrXX19frpZdekiTl5eUF+I0AAACAMA3fkvTwww9rxIgRmjNnjt544w0dO3ZMn3zyiX7+85+rR48eKiws9Bm/cOFCbdq0SfPmzWvxmHV1dVqxYoXuuecedehw5uXuL7zwgh555BF98cUXOn78uLZt26aJEyeqpqZG06dP13XXXXe2XxMAAAAxJGzDd1JSkoqLi3Xfffdp9uzZ6tq1q372s58pKytLmzdv1oABA3zG33TTTerSpUvTTiT+vPTSS/ruu+905513nvH8FRUVevjhh/Xuu+9q5MiR6ty5s6666ipZrVatXbtWS5YsOevvCAAAgNhi8Xg8nlAXEStcLpdSUlJUV1en5OTkUJcDAAAAk4XtzDcAAAAQbQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEnO/Ix1AAD8cLulkhKppkay2SSHQ7JaQ10VAIQ3wjcAIGhOpzRjhrRv38k2u11avFjKyQldXQAQ7lh2AgAIitMp5eb6Bm9Jqq422p3O0NQFAJHA4vF4PKEuIla4XC6lpKSorq5OycnJoS4HAILmdktpac2Dt5fFYsyAV1ayBAUA/GHmGwAQsJKSloO3JHk8UlWVMQ4A0BzhGwAQsJqath0HALGG8A0ACJjN1rbjACDWEL4BAAFzOIw13RaL/36LRUpNNcYBAJojfAMAAma1GtsJSs0DuPfnRYu42RIAWkL4BgAEJSdHKiqS+vb1bbfbjXb2+QaAlrHVoInYahBANOEJlwAQPJ5wCQBoFatVGjMm1FUAQGRh2QkAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYJIOoS4AABCZ3G6ppESqqZFsNsnhkKzWUFcFAOGN8A0ACJrTKc2YIe3bd7LNbpcWL5ZyckJXFwCEO5adAACC4nRKubm+wVuSqquNdqczNHUBQCSweDweT6iLiBUul0spKSmqq6tTcnJyqMsBgKC53VJaWvPg7WWxGDPglZUsQQEAf5j5BgAErKSk5eAtSR6PVFVljAMANEf4BgAErKambccBQKwhfAMAAmazte04AIg1hG8AQMAcDmNNt8Xiv99ikVJTjXEAgOYI3wCAgFmtxnaCUvMA7v150SJutgSAloR1+K6vr9eCBQs0cOBAJSYmqk+fPpoyZYqqq6uDOk5eXp4sFkuLr0WLFrX42RUrVuiyyy5Tp06d1L17d11zzTUqKys7y28GAJErJ0cqKpL69vVtt9uNdvb5BoCWhe1Ddurr63XllVdq+/btKiws1NVXX62dO3dq0qRJGjZsmEpLS5WRkRHw8Xr06KGePXv67evevbvf9lmzZunpp5/Wo48+qnfffVd1dXWaOXOmHA6HXn/9dV1//fWt+m4AEOlycqSJE3nCJQAEK2z3+Z43b57y8/P1zDPPaOrUqU3t27dv15AhQ3T55Zfrgw8+COhYeXl5SktL0yOPPBLw+d977z2NGzdOubm5ev3115vajx07poEDB8rlcunzzz9vMdD7wz7fAAAAsS0sl50cO3ZMBQUFio+P1+233+7Tl5mZqREjRujDDz8MOHy3xpNPPilJuvPOO33aO3bsqEmTJqmurk7PP/98u50fAAAA0Scsw/f777+vw4cPa8iQIercuXOz/lGjRkmS1q5d2y7nd7lc2rBhgywWi0aOHGn6+QEAABCdwjJ8b926VZKUnp7ut9/bvm3btoCPuWXLFk2YMEE2m00JCQlKTU3VpEmT/B6jvLxcbrdbvXr18hv+vecvLy9XY2NjwDUAAAAgtoVl+K7576PRWroRslu3bj7jAlFSUqLs7Gzt2LFDBw8e1NKlS/WPf/xDl156qc+a7mDOX19fr4MHDwZcAwAAAGJbWIbvo0ePSpLi4+P99ickJPiMO5P77rtPH330kaZMmaLu3burY8eOGjdunP70pz/pxIkTmjx5svbv3x/0+c9UQ0NDg1wul88LAAAAsSssw3dSUpIk6fjx4377GxoafMadyZAhQ9S/f3+/7aNGjdKRI0f02muvBX3+M9WQn5+vlJSUpldqampA9QIAACA6hWX47t27tyTpwIEDfvu//fZbSZLNZjvrc3n3Ct+1a1fQ509MTFTXrl1bPPbcuXNVV1fX9KqqqjrregEAABC5wvIhO0OHDpUkVVZW+u33tg8ePPisz+XvhsnMzEzFxcXp66+/1pEjR9SpUye/5x80aJDi4lr+/SUhIcFniQoAAABiW1jOfI8dO1adOnXStm3bdOTIkWb9paWlkqQJEyac8VilpaXq3bt3izdG7tmzR5J0wQUXNLUlJydr9OjR8ng8fh8lH8z5AQAAAK+wDN9JSUmaNm2aGhoatHz5cp++HTt2qKysTCNHjpTD4WhqLy4uVlZWlgoKCnzGHz9+XPv379c777zT7DxbtmxRaWmpkpKSdOONN/r0zZ49W5L04osv+rTX19dr5cqVSk5O1t13331W3xMAAACxJSzDtyQ9/PDDGjFihObMmaM33nhDx44d0yeffKKf//zn6tGjhwoLC33GL1y4UJs2bdK8efN82i0WiyRp6tSpWrFihb766ivV19fr3XffVW5urjp06KClS5c2rfP2Gj9+vO6991699tprevzxx+VyuVRVVaVJkyZp3759Wrp0qXr16tWu/x0AAAAguoRt+E5KSlJxcbHuu+8+zZ49W127dtXPfvYzZWVlafPmzRowYIDP+JtuukldunTRlClTfNp/8pOfaOPGjbrpppv0u9/9ThdeeKFSUlJ05513atSoUdq0aZNuvvlmvzUsXrxYL7/8stauXSubzabMzEwdOXJEJSUlys3NbbfvDgAAgOhk8Xg8nlAXEStcLpdSUlJUV1en5OTkUJcDAAAAk4XtzDcAAAAQbcJyq0EAQPhzu6WSEqmmRrLZJIdDslpDXRUAhDfCNwAgaE6nNGOGtG/fyTa7XVq8WMrJCV1dABDuWHYCAAiK0ynl5voGb0mqrjbanc7Q1AUAkYAbLk3EDZcAIp3bLaWlNQ/eXhaLMQNeWckSFADwh5lvAEDASkpaDt6S5PFIVVXGOABAc4RvAEDAamradhwAxBrCNwAgYDZb244DgFhD+AYABMzhMNZ0Wyz++y0WKTXVGAcAaI7wDQAImNVqbCcoNQ/g3p8XLeJmSwBoCeEbABCUnBypqEjq29e33W432tnnGwBaxlaDJmKrQQDRhCdcAkDweMIlAKBVrFZpzJhQVwEAkYVlJwAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASTqEugAAQGRyu6WSEqmmRrLZJIdDslpDXRUAhDfCNwAgaE6nNGOGtG/fyTa7XVq8WMrJCV1dABDuWHYCAAiK0ynl5voGb0mqrjbanc7Q1AUAkcDi8Xg8oS4iVrhcLqWkpKiurk7JycmhLgcAguZ2S2lpzYO3l8VizIBXVrIEBQD8YeYbABCwkpKWg7ckeTxSVZUxDgDQHOEbABCwmpq2HQcAsYbwDQAImM3WtuMAINYQvgEAAXM4jDXdFov/fotFSk01xgEAmiN8AwACZrUa2wlKzQO49+dFi7jZEgBaQvgGAAQlJ0cqKpL69vVtt9uNdvb5BoCWsdWgidhqEEA04QmXABA8nnAJAGgVq1UaMybUVQBAZGHZCQAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgkrAO3/X19VqwYIEGDhyoxMRE9enTR1OmTFF1dXVQx3n//fd1xx13aMCAAUpMTFRSUpIuuugizZ49W7W1tX4/U1hYKIvF0uIrOzu7Db4hAAAAYknYbjVYX1+vK6+8Utu3b1dhYaGuvvpq7dy5U5MmTdKwYcNUWlqqjIyMMx7n2Wef1bRp0zR48GA9//zzuvTSS+VyufTKK6/ooYce0sqVK7Vx40a/x+rYsaP69evn97h2u/2svyMAAABiS9iG70cffVSlpaV65plnlPPfx6UNHz5cr732moYMGaK8vDx98MEHZzzOsWPHFB8fr3Xr1ik1NVWS1KVLF82ePVsHDhzQb37zG9177716++23m302KytL69evb9PvBQAAgNgVlstOjh07poKCAsXHx+v222/36cvMzNSIESP04YcfBhS+e/furVtuuaUpeJ9qwoQJkqT33ntPbre7bYoHAAAAWhCW4fv999/X4cOHNWTIEHXu3LlZ/6hRoyRJa9euPeOxJk2apMLCQr99KSkpkiSLxSKPx9P6ggEAAIAAhGX43rp1qyQpPT3db7+3fdu2bWd1nl27dkkywnyHDs1X4NTW1uruu+9Wenq6EhMT1bNnT40fP15vvvnmWZ0XAAAAsSksw3dNTY0kqXv37n77u3Xr5jOutZYvXy6LxaL58+f77a+oqFDPnj21YcMGuVwuFRcXS5ImTpyoWbNmnfH4DQ0NcrlcPi8AAADErrAM30ePHpUkxcfH++1PSEjwGdcaf/3rX/Xmm29q5syZuuKKK5r1X3XVVSorK9Pjjz+ufv36KT4+XpmZmVqzZo369Omjp59+WuvWrTvtOfLz85WSktL08rfuHAAAALEjLMN3UlKSJOn48eN++xsaGnzGBevTTz/VbbfdppycHP32t7/1O8ZutysrK6tZe2JioqZMmSJJWrZs2WnPM3fuXNXV1TW9qqqqWlUvAAAAokNYbjXYu3dvSdKBAwf89n/77beSJJvNFvSxKysr9dOf/lSXX365Vq1aJavVGvQxvHuCe9eMtyQhIaFplh4AAAAIy/A9dOhQSUZQ9sfbPnjw4KCOu3v3bl1xxRVyOBxavny535ssA9HY2NiqzwFANHG7pZISqaZGstkkh0NqxXwGAMSUsFx2MnbsWHXq1Enbtm3TkSNHmvWXlpZKOrlPdyB27typn/zkJxo3bpxWrlzpE7yfeuqpZktCevfu7ffBO5K0Z88eSdIFF1wQ8PkBIJo4nVJamjR2rHTrrcZ7WprRDgBoWViG76SkJE2bNk0NDQ1avny5T9+OHTtUVlamkSNHyuFwNLUXFxcrKytLBQUFzY63ZcsWjR49Wrm5ufrDH/6guDjfr/3AAw80BWqv/fv366233mp2rPr6er300kuSpLy8vNZ+RQCIWE6nlJsr7dvn215dbbQTwAGgZRZPmD5d5ujRo7ryyiu1Y8cOvfzyyxo/frwqKir0i1/8Ql999ZU+/PBDDRgwoGn8tddeqz//+c/q3LmzDh061NT+8ccfa9y4cWpoaNB1113n91yrV69WcXGxxowZ09RmsVgUFxen+fPna/LkybLZbNq1a5ceeOABvfvuu5o+fbqWLFkS1HdyuVxKSUlRXV2dkpOTg/svBADCgNttzHB/P3h7WSyS3S5VVrIEBQD8Ccs135Ix+11cXKz8/HzNnj1bN998s7p3767x48frsccek91u9xl/0003aePGjZo8ebJP+8qVK5tu0Fy9enXA56+oqNDq1av1zjvv6IUXXlBtba26du2qSy65RGvXrg1qyQsARIuSkpaDtyR5PFJVlTHulPkMAMB/he3MdzRi5htApFu1yljjfSavvirdckv71wMAkSYs13wDAMJToDu8tmInWACICYRvAEDAHA5jTbfF4r/fYpFSU41xAIDmCN8AgIBZrdLixcZ//n4A9/68aBE3WwJASwjfAICg5ORIRUVS376+7Xa70Z6TE5q6ACAScMOlibjhEkA04QmXABC8sN1qEAAQ3qxWthMEgGCx7AQAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMEmHUBcAAIhMbrdUUiLV1Eg2m+RwSFZrqKsCgPBG+AYABM3plGbMkPbtO9lmt0uLF0s5OaGrCwDCHctOAABBcTql3Fzf4C1J1dVGu9MZmroAIBJYPB6PJ9RFxAqXy6WUlBTV1dUpOTk51OUAQNDcbiktrXnw9rJYjBnwykqWoACAP8x8AwACVlLScvCWJI9HqqoyxgEAmiN8AwACVlPTtuMAINYQvgEAAbPZ2nYcAMQawjcAIGAOh7Gm22Lx32+xSKmpxjgAQHOEbwBAwKxWYztBqXkA9/68aBE3WwJASwjfAICg5ORIRUVS376+7Xa70c4+3wDQMrYaNBFbDQKIJjzhEgCCxxMuAQCtYrVKY8aEugoAiCwsOwEAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABM0iHUBcQSj8cjSXK5XCGuBAAAAKfTpUsXWSyWNj8u4dtE33zzjSQpNTU1xJUAAADgdL766iv16tWrzY9L+DZR9+7dJUlffvmlUlJSQlwNopnL5VJqaqqqqqqUnJwc6nIQxfi3BrPwbw1m8f5bi4+Pb5fjE75NFBdnLLFPSUnhwgFTJCcn828NpuDfGszCvzWYpT2WnEjccAkAAACYhvANAAAAmITwbaKEhAT96le/UkJCQqhLQZTj3xrMwr81mIV/azBLe/9bs3i8+98BAAAAaFfMfAMAAAAmIXwDAAAAJiF8AwAAACYhfANRpL6+Xr/61a80YMAAJSYmqk+fPpoyZYr27dsX6tIAAIC44RKIGvX19bryyitVWloqm80mh8OhvXv36qOPPlKvXr1UVlam/v37h7pMAABiGjPfQJR44oknVFpaqpEjR+qzzz7T6tWr9c9//lMLFy7U119/rSlTpoS6RAAAYh7h2wQbN27UxIkTdf7558tiseiRRx4JdUmIMidOnFBBQYEk6dlnn1Xnzp2b+mbNmqXBgwdr48aN+te//hWqEhFBuGYh1D766CNZLJaAXtddd12oy0UEKCoqUlZWlrp3767ExERlZGTooYce0vHjx02vpYPpZ4xBhw8f1o9+9CPdeuutmjlzZqjLQRT64IMPdPDgQfXv31/Dhg1r1p+bm6tt27bprbfe0vDhw0NQISIJ1yyE2rFjx/Q///M/euWVV9S5c2fdcMMNPv1//OMfdfz4cd12222EbwSke/fumj17ti688EJ16tRJmzdv1l133SWXy6UlS5aYWgvh2wTXXHONrrnmGknSnDlzQlwNotHWrVslSRdffLHffm+7dxxwOlyzEGqjR49Wjx499Morr+iSSy5RYWFhU5/L5dLy5cvVv39/n3bgdK644gqfn9PS0rRx40b97W9/M70Wlp0AUeDLL7+UJNntdr/93nbvOAAId9u2bZMkDR482Kd9y5Yt8ng8yszMDEVZiBIVFRX6y1/+orFjx5p+bsI3EAUOHz4sSUpKSvLb36lTJ59xABDutm/fLknNQvaWLVskNQ/lQCA6d+6shIQE/ehHP9KYMWP09NNPm14Dy06AKODdMdRisZy2HwAixelmvv21A4HYsmWLjh07po8//lj/+7//q/POO08LFiwwtQbCNxAFunTpIkk6cuSI3/6jR49Kks8uKAAQzrZv3664uDgNGjTIp90bvll2gtbIyMiQZPz7iYuL05QpUzR79uymvxCbgWUnYYgtlhCsfv36SVKLT7L0tnvHAW2F6xXaw8GDB1VVVaX+/fv7LKf77rvvtHPnTiUlJfHQMLQJj8ejEydO+LS193WNmW8THD58WJ9//rkk6fjx4/rPf/6jLVu2KD4+Xj/60Y+ajWeLJQRryJAhkqRPPvnEb7+3nT/TIhDBXLO4XqE9eNd7f/+atWvXLjU0NOjSSy9VXBzzhwjcY489pssuu0w/+MEP5PF49NFHH2nOnDmaOHGiunbt6jO23a9rHrS74uJij6Rmr/PPP7/Fz2zfvt0jyTNmzBif9rq6Oo/FYvFkZGS0c9WIJA0NDZ6UlBSPJM8nn3zSrH/w4MEeSZ6PPvooBNUh0gR7zeJ6hbb2zDPPeCR5FixY4NP++uuveyR5Jk2aFKLKEKnmzp3rGTBggKdjx46eLl26eDIzMz1PPvmk5+jRo37Ht+d1jV8bTTBmzBh5PJ5mr71797b4GbZYQjDi4+M1bdo0SdK0adN81n7/7ne/07Zt2/TjH/9Yl156aahKRAQJ9prF9QptraWZb++OTVu3btWf//znkDydEJHpiSee0KeffqqjR4/K5XJp27ZtmjNnjjp27Oh3fHte11h2EoDc3FyVl5cH9Znly5crKyur1edkiyUE66GHHtLf/vY3lZaW6oc//KEcDoe++OIL/fOf/1SPHj20bNmyUJeIKMX1Cm2tpfB95ZVXavDgwaqoqNCcOXOaHgYFtLX2vK4RvgOwd+9effrpp0F9xru7RGuxxRKClZiYqOLiYuXn5+vVV1/VmjVr1K1bN91+++167LHHlJqaGuoSEaW4XqGtffjhh37bU1NTeVIvTNGe1zWLx8MGwOGoX79+qq6u1qFDh3zu9L744ou1efNmffbZZ/rhD38YwgoBwMD1CkC0ac/rGuG7nbX00BN/vP9THDx4UN26ddMPf/hDffbZZ0393333nTp37iyr1apDhw5xpzeANhfsNYvrFYBwF27XNa6G7czfTUstvbzOtMXSRRddxP+RAWgXwV6zuF4BCHfhdl3jihiGWlpntGvXLknSwIEDTa8JAPzhegUg2rT3dY3wbYJly5Zp7Nix6tWrl7p06aLhw4frlVdeaXE8WywBCKVgrllcrwBEgnC6rhG+TfD3v/9dEyZM0Ntvv63Nmzfr5ptv1i9+8QutXr3a7/gzbbG0a9cuzZkzR+ecc0671w4g9gRzzeJ6BSAShNN1jRsuQ+Saa65Rx44d9ac//SnUpQDAGXHNAhBtQnVdY+Y7RA4ePKiePXuGugwACAjXLADRJlTXNcJ3CLz88sv6+OOPddddd4W6FAA4I65ZAKJNKK9rhG+TrV27VnfddZd+//vf6+KLLw51OQBwWlyzAESbUF/XCN8m+uMf/6ibbrpJzz//vCZPnhzqcgDgtLhmAYg24XBdI3yb5MUXX1ReXp4KCwuVl5cX6nIA4LS4ZgGINuFyXWO3ExP87ne/0+zZs/Xss89q4sSJTe1Wq1W9evUKYWUA0BzXLADRJpyua4RvE6SlpemLL75o1n7++edr79695hcEAKfBNQtAtAmn6xrhGwAAADAJa74BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAEH797//rV/+8peaOnWqGhoaNH36dJ177rk699xztWDBAvHwZADwj/ANAAjKe++9p0GDBum5557ThRdeqIcffljPPPOMBg4cqKNHj+qRRx7R//3f/4W6TAAISxYP0xMAgADt2bNHQ4YMUUNDg1599VVdf/316tWrl7p06aK9e/fq9ddf180336yePXvqq6++ksViCXXJABBWmPkGAARszpw5OnLkiGbPnq0bb7xR//rXv3Tw4EGNHj1acXFxuvHGG9W5c2fV1tbqP//5T6jLBYCwQ/gGAATkq6++0tq1a9WlSxfNnz9fkrR582ZJ0tChQyVJcXFx6tGjhyTJ7XaHpE4ACGeEbwBAQN5//3199913uvrqq5WYmChJ2rJliyRpyJAhTeNqa2tlsVh03nnnhaJMAAhrhG8AQEC2bt0qyTdoe9u8M9+7d+/WkSNHdNFFF+mcc84xvUYACHeEbwBAQGprayVJXbt2lSQ1NjZq+/bt6tu3r3r27ClJ+stf/iJJuuaaa0JSIwCEO8I3ACAgXbp0kSR9/vnnkk7Ocntnwo8dO6ZFixYpPj5ed911V8jqBIBwRvgGAATkxz/+sSTJ6XTK5XI1rfceOnSoDh48qBtuuEGVlZWaP3++fvCDH4SwUgAIX+zzDQAIiNvt1siRI7Vp0yb169dP3bp109atW3XxxRfriy++0DfffKNp06apoKAg1KUCQNgifAMAAvbNN9/ogQce0Nq1a3XgwAFJUqdOnXTJJZdo1qxZmjBhQogrBIDwRvgGALSKdz/vb775JsSVAEDkYM03ACBolZWVOnDggIYPHx7qUgAgohC+AQBB+9e//iVJhG8ACBLhGwAQNG/4vvjii0NcCQBEFtZ8AwAAACZh5hsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMMn/BzQT+MFlgOraAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "zp = branchpoint\n",
    "dz = np.array([0.5, 0.0])   # guide direction to find other branch\n",
    "for k in range(10):\n",
    "    ZS = root(sys, zp+dz, args=(zp,ds), jac=dsys, method='hybr') \n",
    "    zs =ZS.x\n",
    "    if d2Epot_dphi2(zs[0], zs[1]) > 0:\n",
    "        stable_phi.append( zs[0] )\n",
    "        stable_mu.append( zs[1] )\n",
    "    else:\n",
    "        nonstable_phi.append( zs[0] )\n",
    "        nonstable_mu.append( zs[1] )\n",
    "    dz = zs -zp\n",
    "    zp = zs\n",
    "\n",
    "plt.plot(stable_phi, stable_mu, 'bo')\n",
    "plt.plot(nonstable_phi, nonstable_mu, 'bo', fillstyle='none')\n",
    "plt.xlabel('$\\\\varphi$')\n",
    "plt.ylabel('$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = ['$-\\\\frac{1}{2}\\pi$', '$0$', '$\\\\frac{1}{2}\\pi$', '$\\pi$', '$\\\\frac{3}{2}\\pi$']\n",
    "plt.xticks(tick_pos, labels)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In accordance with our initial intuitive expectation we recognize the buckled state branching off when the critical load is exceeded.\n",
    "By symmetry there is also a mirrored branch (buckling to the other side). It can be found by slightly modifying the starting guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAImCAYAAACLl+YVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABJy0lEQVR4nO3de3zT1eH/8XcabEuBlqtSSJGODnRSLqKVyzJAHQ90CLXWefkyLfh4eBkgiF/hh4oTddZ9J+NSfeh0CALKUBZBmZu6WaDabuKUS6Eo8ihaug6tSMOtBdP8/vgshdAUktJ+8knyej4eecScc/LJSUlP356cz/nYvF6vVwAAAABaXVy4OwAAAADECsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCt4m8Xq/cbre4rhEAAEBsInyb6NChQ0pJSdGhQ4fC3RUAAACEAeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwSZtwdwAAEJk8HqmoSKqqklJTJadTstvD3SsAsDZmvgEAIXO5pIwMafRo6dZbjfuMDKMcANA0wjcAICQul5SbK2VmSiUl0qFDxn1mplFOAAeAptm8Xq833J2IFW63WykpKaqpqVFycnK4uwMAIfN4jBnuzExp7Vop7pQpnPp6KTtbKi2Vdu9mCQoABMKabyAMWCuLSFVUJO3dK61a5R+8JePxnDnS8OFGu1GjwtFDIDSMxzAb4RswmcslTZ8u7dt3sszhkBYtknJywtcvIBhVVcZ9//6B633lvnaAlTEeIxxY8w2YyLdW9tSBXpIqK1kri8iQmmrcl5YGrveV+9oBVsV4jHBhzbeJWPMd2zweqXfvxgO9j81mzLiUl/OVJ6yLNd+IBozHCCdmvgGTFBU1PdBLktcrVVQY7QCrstul+fOl9euNoH3qbifZ2Ub5008TWGBtjMcIJ9Z8AyYJdg0sa2VhdTk50po10v33GydX+qSnG+WslYXVMR4jnAjfgEmCXQPLWllEgpwcacIEdolAZGI8Rjix5ttErPmObb41hpWVxleap2ONIQCYg/EY4cSab8AkdruxfZVkDOyn8j1euJCBHgBaG+MxwonwDZjIt1a2Z0//coeDtbIAYCbGY4QLy05MxLIT+HBFNQCwBsZjmI3wbSLCNwAAQGxj2QkAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsuH79raWj300EOKj49XXl5eyM/fsGGDbDbbWW+jRo3ye96yZcvO2D47O7tF3h8ARCqPR9qwQVq1yrj3eMLdIwCwvjbh7sCZFBYW6q677lJ1dbVOnDhxTsfq169fwPL6+nrt3r1bF198caO6tm3bqlevXgGf53A4zqk/ABDJXC7p/vulvXtPlvXuLc2fL+XkhKtXAGB9lg3fq1at0j333KNf//rXSkpK0uTJk8/peLt27QpYvnbtWl1//fWaMmVKo7qsrCxt2LDhnF4XAKKNyyXl5krjxhmz3v37S6Wl0pNPGuVr1hDAAaAplg3f6enp2rlzp3r06KFly5Y1+zjdunXTDTfc0GT94sWLNXr0aPXv37/ZrwG0JI9HKiqSqqqk1FTJ6ZTs9nD3CjB4PMaM97hx0tq1Utx/Fy8OHWo8zs6W/vd/pQkT+NzCGhhTYTWWDd9Dhw5tkeNccsklWrNmTcC6HTt2qLCwUC6Xq0VeCzhXLpc0fbq0b9/JModDWrSImURYQ1GRsdRk1aqTwdsnLk6aM0caPtxod9qpNIDpGFNhRZY/4bI1FRQU6MILL9T48eMD1ldXV+vuu+9Wenq6EhMT1bVrV40dO1ZvvvmmyT1FLPB9lX/qHwlJqqw0yvl/RFhBVZVx39SXhb5yXzsgXBhTYVUxG74PHjyolStX6pe//KXsTXz/VFZWpq5du2rjxo1yu90qLCyUJE2YMEEzZ84862vU1dXJ7Xb73YBAPB5jdsbrbVznK5sxg90kEH6pqcZ9aWngel+5rx0QDoypsLKYDd9LlixRfX297rjjjoD1V199tUpKSvTEE0+oV69eio+PV2ZmptauXasePXpowYIFWr9+/RlfIz8/XykpKQ23tLS01ngriAJFRY1nZ07l9UoVFUY7IJycTmNXkyeflOrr/evq66X8fCk93WgHhAtjKqwsJsN3fX29nn32Wd16663q0qVLwDYOh0NZWVmNyhMTExt2Xlm6dOkZX2fOnDmqqalpuFVUVJx75xGVgv2Knq/yEW52u7Gd4Pr1xsmVJSXSoUPGfXa2Uf7005zQhvBiTIWVWfaEy9a0fv16lZeXa+rUqc16fkZGhqSmty/0SUhIUEJCQrNeA7El2K/o+SofVpCTY2wneP/9xsmVPunpbDMIa2BMhZXFZPhevHixnE6nBg0a1Kzn15/+XStwjpxO4wz8ysrAaxRtNqOer/JhFTk5xnaCbOEGK2JMhZXF3LKTsrIy/f3vf9e0adPO2K579+56++23A9bt2bNHknTRRRe1eP8Qm+x2Y+sryfijcCrf44ULCTawFrvd2E7wlluMez6fsArGVFhZ1ITvwsJCZWVlqaCg4IztFi9erJ49e+r6668/Y7v9+/frrbfealReW1url156SZKUl5fX7P4Cp/N9ld+zp3+5w8FX+QAQKsZUWFXULDuZP3++Nm/erLKysiZntWtqarRixQrNmTNHbdqc/a2/8MILuuCCCzRp0iSlpqZq165deuCBB1RVVaVp06bpuuuua+m3gRjHV/kA0HIYU2FFNq830Gooa7Cd/l3RKZYuXeo387xixQpNmTJFkyZN0iLfd02nWbBggebMmaOKigp169btjK+9a9curV69Wu+884727t2r6upqdezYUZdddpnuvvvuJi/McyZut1spKSmqqalRcnJyyM8HAABAZLN0+I42hG8AAIDYFjVrvgEAAACrI3wDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmaRPuDgAAIpPHIxUVSVVVUmqq5HRKdnu4ewUA1sbMNwAgZC6XlJEhjR4t3XqrcZ+RYZQDAJpG+AYikMcjbdggrVpl3Hs84e4RYonLJeXmSpmZUkmJdOiQcZ+ZaZQTwGEWxkJEIpvX6/WGuxOxwu12KyUlRTU1NUpOTg53dxChXC5p+nRp376TZQ6HtGiRlJMTvn4hNng8xgx3Zqa0dq0Ud8oUTn29lJ0tlZZKu3ezBAWti7EQkYqZbyCC+GYcT/1jI0mVlcw4whxFRdLevdKDD/oHb8l4PGeOVF5utANaC2MhIhnhG4gQHo8xyxPouypf2YwZfO2K1lVVZdz37x+43lfuawe0NMZCRDrCNxAhiooaz/KcyuuVKiqYcUTrSk017ktLA9f7yn3tgJbGWIhIR/gGIkSwM4nMOKI1OZ1S797Sk08aa7xPVV8v5edL6elGO6A1MBYi0hG+gQgR7EwiM45oTXa7NH++tH69cXLlqbudZGcb5U8/zcmWaD2MhYh07HZiInY7wbnweIwZx8rKwGsdbTbjTP/ycoIPWp/LJd1/v3HypU96uhG82WkCrYmxEJGOmW8gQtjtxhZakvHH5VS+xwsX8scG5sjJkb74QioslF591bjfvZvgjdbHWIhIx8y3iZj5RksItLdtWprxx4bgAyBWMBYiUhG+TUT4RkvxeIwz+auqjHWNTiezPABiD2MhIhHh20SEbwAAgNjGmm8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJG3C3QEAQGTyeKSiIqmqSkpNlZxOyW4Pd68AwNoI30AUIxyhtbhc0v33S3v3nizr3VuaP1/KyQlXrxAtGLsQzVh2AkQpl8sIQ6NHS7featz37m2UA+fC5ZJyc6XMTKmkRDp0yLjPzDTK+YzhXDB2IdrZvF6vN9ydiBVut1spKSmqqalRcnJyuLuDKOYLR6f/dttsxv2aNcxOonk8Hikjwwjaa9dKcadM4dTXS9nZUmmptHs3M5UIHWMXYgHh20SEb5jB4zFmifbtC1xvs0kOh1ReTjhC6DZsMGYiS0qkoUMb15eUSMOHS4WF0qhRZvcOkYyxC7GCZSdAlCkqavqPl2TMKFVUGO2AUFVVGff9+weu95X72gHBYuxCrCB8A1Em2NBDOEJzpKYa96Wlget95b52QLAYuxArLB++a2tr9dBDDyk+Pl55eXnNOkZeXp5sNluTt4ULFzb53BUrVuiKK65Qu3bt1LlzZ1177bUqKSlp3psBTBBs6CEcoTmcTmNpwJNPGmu8T1VfL+XnS+npRjsgFIxdiBWWDt+FhYUaMGCAnnvuOZ04ceKcjtWlSxf169cv4K1z584BnzNz5kzddtttGjdunP79739ry5YtSkxMlNPp1BtvvHFO/QFai9NprIv0naB0OptNSksjHKF57HZjO8H1642TK0/d7SQ72yh/+mnW5CJ0jF2IFZY94XLVqlW655579Otf/1pJSUmaPHmybr/9di1btizkY+Xl5al379569NFHg37Oe++9pzFjxig3N1evv/56Q/mxY8fUr18/ud1uffHFF+ratWvQx+SES5jFt2OA5L9rADsGoKUE2uc7Pd0I3ny20FyMXYgFlp35Tk9P186dOzVlyhTZmvrf4Fb01FNPSZLuvPNOv/K2bdtq4sSJqqmp0fPPP296v4Bg5OQYf6R69vQvdzj444WWkZMjffGFsavJq68a97t389nCuWHsQiyw7BUuhwbaw8okbrdbGzdulM1m07BhwxrVjxgxQpK0bt06Pfzww2Z3DwhKTo40YQJXiUPrsdvZThAtj7EL0c6y4bulbdmyRePHj9fmzZt14MABnX/++Ro5cqRmzZqlAQMG+LUtLS2Vx+PR+eefr/bt2zc6Vnp6ekO7+vp6xcVZ9gsExDjCEYBIxNiFaBYzqbGoqEjZ2dnasWOHDh48qCVLlugf//iHLr/8cr813ZJU9d99jJo6EbNTp06SjJ1YDh482ORr1tXVye12+90AAAAQu2IifN9333366KOPNHnyZHXu3Flt27bVmDFj9Kc//UknTpzQpEmTtH///ob2R48elSTFx8cHPF5CQkKjtoHk5+crJSWl4ZaWltZC7wgAAACRKCbC98CBA9WnT5+A5SNGjNCRI0f02muvNZQnJSVJko4fPx7weHV1dY3aBjJnzhzV1NQ03CoqKpr7FgAAABAFYiJ8n0lGRoYkadeuXQ1l3bt3lyQdOHAg4HO+++47SVJiYqI6duzY5LETEhKUnJzsdwMAAEDsivnwXX/6JdokZWZmKi4uTt98842OHDnSqL68vFyS1L9/f062BAAAQNCiPjkWFxere/fuTZ4YuWfPHknSRRdd1FCWnJyskSNHyuv1BryUfHFxsSRp/PjxLd9hAAAARK2oCd+FhYXKyspSQUGBX/nx48e1f/9+vfPOO42es2XLFhUXFyspKUk33nijX92sWbMkSS+++KJfeW1trVauXKnk5GTdfffdLfwuAAAAEM2iJnzPnz9fmzdv1oMPPuhX7rs65pQpU7RixQp9/fXXqq2t1bvvvqvc3Fy1adNGS5YsaVjn7TN27Fjde++9eu211/TEE0/I7XaroqJCEydO1L59+7RkyRJ169bNtPcHAACAyGfzer3ecHeiKWe6rPzSpUuVl5fX8HjFihWaMmWKJk2apEWLFjWUe71effDBB/rjH/+o4uJiffXVVzp8+LBSU1M1cuRIzZw5UwMHDmzydZYvX66CggLt3LlT5513noYNG6a5c+dq+PDhIb8ft9utlJQU1dTUcPIlAABADLJ0+I42hG8AAIDYFjXLTgAAAACrI3wDAAAAJmkT7g4AsA6PRyoqkqqqpNRUyemU7PZw9wpApGNsAU4ifAOQJLlc0vTp0r59J8scDmnRIiknJ3z9gnURqBAMxhbAHydcmogTLmFVLpeUmyudPhr4Nhxas4Y/kvDnckn33y/t3XuyrHdvaf58Pis4ibEFaIw130CM83iMWalA/xvuK5sxw2gHSCcDVWamVFIiHTpk3GdmGuUuV7h7CCtgbAECY+bbRMx8w4o2bJBGjz57u8JCadSo1u4NrM7jkTIyjKC9dq0Ud8oUTn29lJ0tlZZKu3ezBCXWMbYAgTHzDcS4qqqWbYfoVlRkLDV58EH/4C0Zj+fMkcrLjXaIbYwtQGCEbyDGpaa2bDtEN19Q6t8/cL2vnEAFxhYgMMI3EOOcTmPnAd8JUKez2aS0NKMd4AtKpaWB633lBCowtgCBEb6BGGe3G1t+SY3/SPoeL1zI+l0YnE5jV5MnnzTWeJ+qvl7Kz5fS0wlUYGwBmkL4BqCcHGPLr549/csdDrYCgz+73dhOcP164+TKU3c7yc42yp9+mkAFA2ML0Bi7nZiI3U5gdVw0BcEKtM93eroRvAlUOB1jC3AS4dtEhG8A0YRABQCh4/LyAIBmsdvZnxkAQsWabwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJG3C3QEAkc/jkYqKpKoqKTVVcjoluz3cvQLQ0vhdB84d4RvAOXG5pOnTpX37TpY5HNKiRVJOTvj6BaBl8bsOtAyb1+v1hrsTscLtdislJUU1NTVKTk4Od3eAc+ZySbm50umjiM1m3K9Zwx9lIBrwuw60HMK3iQjfiCYej9S7t/8s2KlsNmNWrLycr6WjFUsQYgO/60DL4oRLAM1SVNT0H2PJmCGrqDDaIfq4XFJGhjR6tHTrrcZ9RoZRjujC7zrQsgjfAJqlqqpl2yFy+JYgZGZKJSXSoUPGfWamUU4Ajy78rgMti/ANoFlSU1u2HSKDxyPdf780bpy0dq00dKjUvr1xv3atUf6//2u0Q3Tgdx1oWYRvAM3idBrrPH0nXJ3OZpPS0ox2iB5FRdLevdKDD0pxp/0FiYuT5swx1v6yBCF68LsOtCzCN4BmsduNLcakxn+UfY8XLuQErGjjW1rQv3/gel85SxCiB7/rQMsifANotpwcY4uxnj39yx0Oth6LVr6lBaWlget95SxBiC78rgMth60GTcRWg4hWbDkXOzweY1eTzExjjfepS0/q66XsbCOA797NZyAa8bsOnDvCt4kI3wCigW+3k3HjjDXe/fsbgTs/X1q/nplQADgTLi8PAAiJbwnC/fdLw4efLE9PJ3gDwNkw820iZr4BRBOWIABA6Jj5BgA0i90ujRoV7l4AQGRhtxMAAADAJIRvAAAAwCSWD9+1tbV66KGHFB8fr7y8vGYd4/3339cdd9yhvn37KjExUUlJSbrkkks0a9YsVVdXB3zOsmXLZLPZmrxlZ2c3/00BAAAgJll6zXdhYaHuuusuVVdX68SJE806xrPPPqupU6dqwIABev7553X55ZfL7XbrlVde0cMPP6yVK1dq06ZNysjIaPTctm3bqlevXgGP63A4mtUfAAAAxC7Lhu9Vq1bpnnvu0a9//WslJSVp8uTJzTrOsWPHFB8fr/Xr1ystLU2S1KFDB82aNUsHDhzQb37zG9177716++23Gz03KytLGzZsOJe3AQAAADSw7LKT9PR07dy5U1OmTJHNZmv2cbp3765bbrmlIXifavz48ZKk9957Tx6Pp9mvAQAAAATDsjPfQ4cObZHjTJw4URMnTgxYl5KSIkmy2Wxiu3MAAAC0NsvOfJth165dkqQRI0aoTZvG/x9SXV2tu+++W+np6UpMTFTXrl01duxYvfnmm2Z3FQAAAFEgpsP38uXLZbPZNHfu3ID1ZWVl6tq1qzZu3Ci3263CwkJJ0oQJEzRz5kwzuwpEPY9H2rBBWrXKuGclGNDy+D0Dws+yy05a21//+le9+eabuu+++3TllVc2qr/66qtVUlKirKyshrLMzEytXbtWffr00YIFC3TllVdq3LhxTb5GXV2d6urqGh673e6WfRNAlHC5pOnTpX37TpY5HNKiRVJOTvj6BUQTfs8Aa4jJme/PPvtMt912m3JycvTb3/42YBuHw+EXvH0SExMbdl5ZunTpGV8nPz9fKSkpDbdAJ30Csc7lknJz/QOBJFVWGuUuV3j6BUQTfs8A64i58F1eXq6f/vSnGj58uFatWiW73R7yMXx7gvvWjDdlzpw5qqmpabhVVFQ0q89AtPJ4jJm4QOc7+8pmzOCrceBc8HsGWEtMhe/du3frJz/5iX784x9rzZo1io+Pb9Zx6uvrg2qXkJCg5ORkvxuAk4qKGs/EncrrlSoqjHYAmoffM8BaYiZ879y5Uz/5yU80ZswYrVy50m93k6effrrRrHT37t0DXnhHkvbs2SNJuuiii1qvw0AMqKpq2XYAGuP3DLCWqAnfhYWFysrKUkFBQaO6LVu2aOTIkcrNzdUf/vAHxcX5v+0HHnigIVD77N+/X2+99VajY9XW1uqll16SJOXl5bXcGwBiUGpqy7aDudg5IzLwewZYS9TsdjJ//nxt3rxZZWVlmjZtWkP5xx9/rDFjxqiurk7ffPONbrnllqCP+cILL+iCCy7QpEmTlJqaql27dumBBx5QVVWVpk2bpuuuu6413goQM5xOY7eFysrA61FtNqPe6TS/bzgzl0u6/35p796TZb17S/Pns3OG1fB7BliLpcP36ZeVf/nll/Xyyy9LMnYaOXXm+aabbtKmTZs0adIkv+esXLlS3333nSRp9erVQb92WVmZVq9erXfeeUcvvPCCqqur1bFjR1122WVat25dw6XpATSf3W5sc5abawSAU4OB79d/4UKjHazDt3PGuHHGrHf//lJpqfTkk0b5mjUEcCvh9wywFpuX66qbxu12KyUlRTU1NZx8CZwi0P7DaWlGICDEWYvHI2VkSJmZ0tq10qmr+OrrpexsI4jv3k2Ysxp+zwBrIHybiPANNM3jMXZbqKoy1p46nYQ3K9qwQRo9WiopkYYObVxfUiINHy4VFkqjRpndO5wNv2dA+Fl62QmA2GG3E9YigW9HjP79A9f7ytk5w5r4PQPCL2p2OwEAtD7fjhilpYHrfeXsnAEAgbHsxEQsOwEQ6VjzDQDnhplvAEDQ7HZjO8H1642gXVIiHTpk3GdnG+VPP03wBoCmsOYbABCSnBxjO8H77zdOrvRJT2ebQQA4G5admIhlJwCiCTtnAEDomPkGADQLO2cAQOhY8w0AAACYhPANAAAAmITwDQAAAJiENd8AohInAyLS8RkGohPhG0DUcbmk6dOlfftOljkc0qJFbIOHyMBnGIhebDVoIrYaBFqfyyXl5kqnj2w2m3HPPtSwOj7DQHQjfJuI8A20Lo9H6t3bf7bwVDabMXtYXs7X97AmPsNA9OOESwBRo6io6dAiGTOJFRVGO8CK+AwD0Y/wDSBqVFW1bDvAbHyGgehH+AYQNVJTW7YdYDY+w0D0Y823iVjzDbQu33rZysrGJ6tJrJeF9fEZBqIfM98AoobdbmzFJp3cGcLH93jhQkILrIvPMBD9CN8AokpOjrEVW8+e/uUOB1u0ITLwGQaiG8tOTMSyE8A8XB0QkY7PMBCdCN8mInwDAADENi4vDwBoFmZmASB0rPkGAITM5ZIyMqTRo6VbbzXuMzKMcgBA0wjfAICQuFxSbq6UmSmVlEiHDhn3mZlGOQEcAJrGmm8TseYbQKTzeIwZ7sxMae1aKe6UKZz6eik7WyotlXbvZgkKAATCzDcAIGhFRdLevdKDD/oHb8l4PGeOcQGYoqKwdA8ALI/wDQAIWlWVcd+/f+B6X7mvHQDAH7udAEAzxeJuH6mpxn1pqTR0aOP60lL/dtEqFv/tAbQM1nybiDXfQPRwuaTp06V9+06WORzGpcGj+QqErPmO3X97AC2DZScAECLfbh+nhi9JqqyM/t0+7HZp/nxp/XojaJ+620l2tlH+9NPRHbxj9d8eQMtg5ttEzHwDkc/jkXr3bhy+fGw2Yxa0vDx6A6hkhMz77zdOvvRJTzeCd7TO/vJvD6AlEL5NRPgGIt+GDcYFZc6msFAaNaq1exNesbbumX97AC2BEy4BIATB7uIRC7t92O2xFTL5twfQEljzDQAhCHYXj2jf7SMW8W8PoCWw7MRELDsBIp9v3W9lpRRo9GTdb/Ti3x5AS2DmGwBCYLcbW8pJRtg6le/xwoWEr2jEvz2AlkD4BoAQ5eRIa9ZIPXv6lzscRnm07vYB/u0BnDuWnZiIZSdAdIm13T5wEv/2AJqL8G0iwjcAAEBsY9kJAAAAYBL2+QaACMOSh5P4WQCINIRvAIggLpc0fbr/Jc4dDmMXjlg72Y+fBYBIZPllJ7W1tXrooYcUHx+vvLy8czrOvHnz1K9fPyUmJqpHjx6aPHmyKisrz/i8FStW6IorrlC7du3UuXNnXXvttSopKWl2PwCguVwuKTfXP2xKxr7TublGfazgZwEgUlk6fBcWFmrAgAF67rnndOLEiWYfp7a2VldddZXmz5+v/Px8fffdd3rrrbdUUlKiwYMH64svvgj4vJkzZ+q2227TuHHj9O9//1tbtmxRYmKinE6n3njjjWb3BwBC5fEYs7yBTpH3lc2YYbSLdvwsAEQyy4bvVatW6frrr9f06dM1f/78czrWY489puLiYuXn5ysnJ0dt27bVkCFD9Nprr6m6ujrgjPp7772nBQsWKDc3V3PnzlVKSop69eqlV155RT169NCkSZNUXV19Tv0CgGAVFTWe5T2V1ytVVBjtoh0/CwCRzLLhOz09XTt37tSUKVNkO/1SYiE4duyYCgoKFB8fr9tvv92vLjMzU0OHDtWHH36oDz74wK/uqaeekiTdeeedfuVt27bVxIkTVVNTo+eff77Z/QKAUFRVtWy7SMbPAkAkCyl8l5SU6MiRI63VFz9Dhw5Vjx49zvk477//vg4fPqyBAweqffv2jepHjBghSVq3bl1Dmdvt1saNG2Wz2TRs2LCgngMArSk1tWXbRTJ+FgAiWUi7nYwYMUJxcXHq06ePBg0a5HdLPWWUu/POOzVkyBDdddddLd7hUG3dulWSMZMeiK9827ZtDWWlpaXyeDw6//zzAwZ233NKS0tVX1+vuDjLfoEAIEo4ncZOHpWVgdc622xGvdNpXp/Ctc2fFX8WABCskML3Lbfcoq1bt+qzzz7T7t279frrrzcsCenatasGDRqkPn366I033tBbb71lifBd9d/vHTt37hywvlOnTn7tQnlObW2tDh482GS7uro61dXVNTx2u90h9h4ADHa7sYVebq4RLk8Nnb6VeQsXmrfHdTi3+bPazwIAQhHSlO0rr7yi0tJSHT58WB999JFeeOEF3XPPPRo2bJgOHjyo9957T7///e/17bffKj4+vrX6HJKjR49KUpP9SUhI8GsXynNOf97p8vPzlZKS0nBLS0sLrfMAcIqcHGnNGqlnT/9yh8MoN2tvayts82eVnwUAhKpZF9lJSEjQZZddpssuu6yhrKamRkuXLtVDDz2kzMxMLV++vMU6eS6SkpIkScePHw9Y75uZ9rUL5TmnP+90c+bM0cyZMxseu91uAjiAc5KTI02YEL6rOp5tmz+bzdjmb8KE1u9TuH8WANAcLXaFy5SUFM2YMUN9+vRRdna2tm/frr59+7bU4Zute/fukqQDBw4ErP/uu+8kyW/NerDPSUxMVMeOHZt87YSEBL9ZcgBoCXa7NGpUeF47lG3+zOhjOH8WANAcLX6m4HXXXad+/frpySefbOlDN8ugQYMkSeXl5QHrfeUDBgxoKMvMzFRcXJy++eabgLu7+J7Tv39/TrYEEFPY5g8Azk2rJMcf/OAHKisra41Dh2z06NFq166dtm3bFjBIFxcXS5LGjx/fUJacnKyRI0fK6/UGvJR8oOcAQCxgmz8AODchhe/p06dr2bJl2rp1q77//vsm2+3Zs6dh6YZZCgsLlZWVpYKCAr/ypKQkTZ06VXV1dY3Woe/YsUMlJSUaNmyYnKftSTVr1ixJ0osvvuhXXltbq5UrVyo5OVl33313K7wTALAu3zZ/TV37zGaT0tLY5g8AmhLSmu+CgoKGrQXPO+88XXLJJRo8eLAGDx6s/v37KyEhQStXrtTnn3+u5557rlU63JT58+dr8+bNKisr07Rp0/zqHnnkEW3cuFGzZ89W9+7dNXbsWJWVlekXv/iFunTpomXLljU63tixY3Xvvfdq8eLFyszM1L333quamhrdd9992rdvn1avXq1u3bqZ9O4AwBrY5g8Azo3N6w10znpga9eu1aefftpwq6ysNA5y2hRIRkaG8vLyNHjwYF166aU6//zzm9e5M1xWfunSpcrLy2t4vGLFCk2ZMkWTJk3SokWLGrWvra1Vfn6+Xn31VX311Vfq3Lmzxo4dq8cff1wOh6PJ11m+fLkKCgq0c+dOnXfeeRo2bJjmzp2r4cOHh/x+3G63UlJSVFNTo+Tk5JCfDwBWEWif77Q0I3izzR8ANC2k8H266upqvzD+6aef6osvvlB9fb1fcE5NTdWQIUNi/nLshG8A0SRcV7gEgEh2TuE7kCNHjmjr1q1+gXzHjh06ceKEPB5PS75UxCF8AwAAxLYWD9+BfP/999qxY4cGDhzY2i9laYRvAACA2GbKJtVt2rSJ+eANAAAAcIUYAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJG3C3QEAQGTyeKSiIqmqSkpNlZxOyW4Pd68AwNoI3wCAkLlc0vTp0r59J8scDmnRIiknJ3z9AgCrY9kJACAkLpeUm+sfvCWpstIod7nC0y8AiAQ2r9frDXcnYoXb7VZKSopqamqUnJwc7u4AQMg8Hql378bB28dmM2bAy8tZggIAgTDzDQAIWlFR08FbkrxeqaLCaAcAaIzwDQAIWlVVy7YDgFhD+AYABC01tWXbAUCsIXwDAILmdBprum22wPU2m5SWZrQDADRG+AYABM1uN7YTlBoHcN/jhQs52RIAmkL4BgCEJCdHWrNG6tnTv9zhMMrZ5xsAmsZWgyZiq0EA0YQrXAJA6LjCJQCgWex2adSocPcCACILy04AAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJO0CXcHAACRyeORioqkqiopNVVyOiW7Pdy9AgBrI3wDAELmcknTp0v79p0sczikRYuknJzw9QsArI5lJwCAkLhcUm6uf/CWpMpKo9zlCk+/ACASWDp819bWat68eerXr58SExPVo0cPTZ48WZWVlUEfY8OGDbLZbGe9jRo1yu95y5YtO2P77Ozsln2zABABPB5jxtvrbVznK5sxw2gHAGjMsstOamtrddVVV2n79u1atmyZrrnmGu3cuVMTJ07U4MGDVVxcrIyMjKCP169fv4Dl9fX12r17ty6++OJGdW3btlWvXr0CPs/hcAT92gAQLYqKGs94n8rrlSoqjHanzWkAAGTh8P3YY4+puLhYzzzzjHL+u4BwyJAheu211zRw4EDl5eXpgw8+CPp4u3btCli+du1aXX/99ZoyZUqjuqysLG3YsKFZ/QeAaFRV1bLtACDWWHLZybFjx1RQUKD4+HjdfvvtfnWZmZkaOnSoPvzww6DCd7du3XTDDTc0Wb948WKNHj1a/fv3P+d+A0C0S01t2XYAEGssGb7ff/99HT58WAMHDlT79u0b1Y8YMUKStG7durMe65JLLtGaNWsC1u3YsUOFhYWaNm3auXUYAGKE02nsamKzBa632aS0NKMdAKAxS4bvrVu3SpLS09MD1vvKt23bdk6vU1BQoAsvvFDjx48PWF9dXa27775b6enpSkxMVNeuXTV27Fi9+eab5/S6ABCp7HZjO0GpcQD3PV64kP2+AaAplgzfVf9dLNi5c+eA9Z06dfJr1xwHDx7UypUr9ctf/lL2Jv5KlJWVqWvXrtq4caPcbrcKCwslSRMmTNDMmTPP+hp1dXVyu91+NwCIdDk50po1Us+e/uUOh1HOPt8A0DRLnnB59OhRSVJ8fHzA+oSEBL92zbFkyRLV19frjjvuCFh/9dVXq6SkRFlZWQ1lmZmZWrt2rfr06aMFCxboyiuv1Lhx45p8jfz8fM2bN6/ZfQQAq8rJkSZM4AqXABAqS858JyUlSZKOHz8esL6urs6vXajq6+v17LPP6tZbb1WXLl0CtnE4HH7B2ycxMVGTJ0+WJC1duvSMrzNnzhzV1NQ03CoqKprVXwCwIrvd2E7wlluMe4I3AJydJWe+u3fvLkk6cOBAwPrvvvtOkpTazNPp169fr/Lyck2dOrVZz/ftL97U9oU+CQkJDbP0AAAAgCVnvgcNGiRJKi8vD1jvKx8wYECzjr948WI5nc6G1wlVfX19s54HAACA2GbJ8D169Gi1a9dO27Zt05EjRxrVFxcXS1KTu5ScSVlZmf7+97+fdXvB7t276+233w5Yt2fPHknSRRddFPLrAwAAIHZZMnwnJSVp6tSpqqur0/Lly/3qduzYoZKSEg0bNkzOUzaSLSwsVFZWlgoKCs547MWLF6tnz566/vrrz9hu//79euuttxqV19bW6qWXXpIk5eXlBfmOAAAAAIuGb0l65JFHNHToUM2ePVtvvPGGjh07pk8++UQ///nP1aVLFy1btsyv/fz587V582Y9+OCDTR6zpqZGK1as0D333KM2bc6+3P2FF17Qo48+qi+//FLHjx/Xtm3bNGHCBFVVVWnatGm67rrrzvVtAgAAIIZYNnwnJSWpsLBQ9913n2bNmqWOHTvqZz/7mbKysvTpp5+qb9++fu1vuukmdejQoWEnkkBeeuklff/997rzzjvP+vplZWV65JFH9O6772rYsGFq3769rr76atntdq1bt06LFy8+5/cIAACA2GLzer3ecHciVrjdbqWkpKimpkbJycnh7g4AAABMZtmZbwAAACDaEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAk5z9GusAAATg8UhFRVJVlZSaKjmdkt0e7l4BgLURvgEAIXO5pOnTpX37TpY5HNKiRVJOTvj6BQBWx7ITAEBIXC4pN9c/eEtSZaVR7nKFp18AEAlsXq/XG+5OxAq3262UlBTV1NQoOTk53N0BgJB5PFLv3o2Dt4/NZsyAl5ezBAUAAmHmGwAQtKKipoO3JHm9UkWF0Q4A0BjhGwAQtKqqlm0HALGG8A0ACFpqasu2A4BYQ/gGAATN6TTWdNtsgettNiktzWgHAGiM8A0ACJrdbmwnKDUO4L7HCxdysiUANIXwDQAISU6OtGaN1LOnf7nDYZSzzzcANI2tBk3EVoMAoglXuASA0HGFSwBAs9jt0qhR4e4FAEQWlp0AAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACZpE+4OAAAik8cjFRVJVVVSaqrkdEp2e7h7BQDWRvgGAITM5ZKmT5f27TtZ5nBIixZJOTnh6xcAWB3LTgAAIXG5pNxc/+AtSZWVRrnLFZ5+AUAksHm9Xm+4OxEr3G63UlJSVFNTo+Tk5HB3BwBC5vFIvXs3Dt4+NpsxA15ezhIUAAiEmW8AQNCKipoO3pLk9UoVFUY7AEBjhG8AQNCqqlq2HQDEGsI3ACBoqakt2w4AYg3hGwAQNKfTWNNtswWut9mktDSjHQCgMcI3ACBodruxnaDUOID7Hi9cyMmWANAUS4fv2tpazZs3T/369VNiYqJ69OihyZMnq7KyMqTj5OXlyWazNXlbuHBhk89dsWKFrrjiCrVr106dO3fWtddeq5KSknN8ZwAQuXJypDVrpJ49/csdDqOcfb4BoGmWvchObW2trrrqKm3fvl3Lli3TNddco507d2rixIkaPHiwiouLlZGREfTxunTpoq5duwas69y5c8DymTNnasGCBXrsscf07rvvqqamRjNmzJDT6dTrr7+u66+/vlnvDQAiXU6ONGECV7gEgFBZdp/vBx98UPn5+XrmmWc0ZcqUhvLt27dr4MCBGj58uD744IOgjpWXl6fevXvr0UcfDfr133vvPY0ZM0a5ubl6/fXXG8qPHTumfv36ye1264svvmgy0AfCPt8AAACxzZLLTo4dO6aCggLFx8fr9ttv96vLzMzU0KFD9eGHHwYdvpvjqaeekiTdeeedfuVt27bVxIkTVVNTo+eff77VXh8AAADRx5Lh+/3339fhw4c1cOBAtW/fvlH9iBEjJEnr1q1rldd3u93auHGjbDabhg0bZvrrAwAAIDpZMnxv3bpVkpSenh6w3le+bdu2oI+5ZcsWjR8/XqmpqUpISFBaWpomTpwY8BilpaXyeDzq1q1bwPDve/3S0lLV19cH3QcAAADENkuG76r/XhqtqRMhO3Xq5NcuGEVFRcrOztaOHTt08OBBLVmyRP/4xz90+eWX+63pDuX1a2trdfDgwaD7AAAAgNhmyfB99OhRSVJ8fHzA+oSEBL92Z3Pffffpo48+0uTJk9W5c2e1bdtWY8aM0Z/+9CedOHFCkyZN0v79+0N+/bP1oa6uTm632+8GAACA2GXJ8J2UlCRJOn78eMD6uro6v3ZnM3DgQPXp0ydg+YgRI3TkyBG99tprIb/+2fqQn5+vlJSUhltaWlpQ/QUAAEB0smT47t69uyTpwIEDAeu/++47SVJqauo5v5Zvr/Bdu3aF/PqJiYnq2LFjk8eeM2eOampqGm4VFRXn3F8AAABELkteZGfQoEGSpPLy8oD1vvIBAwac82sFOmEyMzNTcXFx+uabb3TkyBG1a9cu4Ov3799fcXFN//9LQkKC3xIVAAAAxDZLznyPHj1a7dq107Zt23TkyJFG9cXFxZKk8ePHn/VYxcXF6t69e5MnRu7Zs0eSdNFFFzWUJScna+TIkfJ6vQEvJR/K6wMAAAA+lgzfSUlJmjp1qurq6rR8+XK/uh07dqikpETDhg2T0+lsKC8sLFRWVpYKCgr82h8/flz79+/XO++80+h1tmzZouLiYiUlJenGG2/0q5s1a5Yk6cUXX/Qrr62t1cqVK5WcnKy77777nN4nAAAAYoslw7ckPfLIIxo6dKhmz56tN954Q8eOHdMnn3yin//85+rSpYuWLVvm137+/PnavHmzHnzwQb9ym80mSZoyZYpWrFihr7/+WrW1tXr33XeVm5urNm3aaMmSJQ3rvH3Gjh2re++9V6+99pqeeOIJud1uVVRUaOLEidq3b5+WLFmibt26terPAAAAANHFsuE7KSlJhYWFuu+++zRr1ix17NhRP/vZz5SVlaVPP/1Uffv29Wt/0003qUOHDpo8ebJf+U9+8hNt2rRJN910k373u9/p4osvVkpKiu68806NGDFCmzdv1s033xywD4sWLdLLL7+sdevWKTU1VZmZmTpy5IiKioqUm5vbau8dAAAA0cnm9Xq94e5ErHC73UpJSVFNTY2Sk5PD3R0AAACYzLIz3wAAAEC0seRWgwAA6/N4pKIiqapKSk2VnE7Jbg93rwDA2gjfAICQuVzS9OnSvn0nyxwOadEiKScnfP0CAKtj2QkAICQul5Sb6x+8Jamy0ih3ucLTLwCIBJxwaSJOuAQQ6TweqXfvxsHbx2YzZsDLy1mCAgCBMPMNAAhaUVHTwVuSvF6posJoBwBojPANAAhaVVXLtgOAWEP4BgAELTW1ZdsBQKwhfAMAguZ0Gmu6bbbA9TablJZmtAMANEb4BgAEzW43thOUGgdw3+OFCznZEgCaQvgGAIQkJ0das0bq2dO/3OEwytnnGwCaxlaDJmKrQQDRhCtcAkDouMIlAKBZ7HZp1Khw9wIAIgvLTgAAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAk7QJdwcAAJHJ45GKiqSqKik1VXI6Jbs93L0CAGsjfAMAQuZySdOnS/v2nSxzOKRFi6ScnPD1CwCsjmUnAICQuFxSbq5/8Jakykqj3OUKT78AIBLYvF6vN9ydiBVut1spKSmqqalRcnJyuLsDACHzeKTevRsHbx+bzZgBLy9nCQoABMLMNwAgaEVFTQdvSfJ6pYoKox0AoDHCNwAgaFVVLdsOAGIN4RsAELTU1JZtBwCxhvANAAia02ms6bbZAtfbbFJamtEOANAY4RsAEDS73dhOUGocwH2PFy7kZEsAaArhGwAQkpwcac0aqWdP/3KHwyhnn28AaBpbDZqIrQYBRBOucAkAoeMKlwCAZrHbpVGjwt0LAIgsLDsBAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExi6fBdW1urefPmqV+/fkpMTFSPHj00efJkVVZWhnSc999/X3fccYf69u2rxMREJSUl6ZJLLtGsWbNUXV0d8DnLli2TzWZr8padnd0C7xAAAACxxLJbDdbW1uqqq67S9u3btWzZMl1zzTXauXOnJk6cqMGDB6u4uFgZGRlnPc6zzz6rqVOnasCAAXr++ed1+eWXy+1265VXXtHDDz+slStXatOmTQGP1bZtW/Xq1SvgcR0Oxzm/RwAAAMQWy4bvxx57TMXFxXrmmWeU89/LpQ0ZMkSvvfaaBg4cqLy8PH3wwQdnPc6xY8cUHx+v9evXKy0tTZLUoUMHzZo1SwcOHNBvfvMb3XvvvXr77bcbPTcrK0sbNmxo0fcFAACA2GXJZSfHjh1TQUGB4uPjdfvtt/vVZWZmaujQofrwww+DCt/du3fXLbfc0hC8TzV+/HhJ0nvvvSePx9MynQcAAACaYMnw/f777+vw4cMaOHCg2rdv36h+xIgRkqR169ad9VgTJ07UsmXLAtalpKRIkmw2m7xeb/M7DAAAAATBkuF769atkqT09PSA9b7ybdu2ndPr7Nq1S5IR5tu0abwCp7q6WnfffbfS09OVmJiorl27auzYsXrzzTfP6XUBAAAQmywZvquqqiRJnTt3DljfqVMnv3bNtXz5ctlsNs2dOzdgfVlZmbp27aqNGzfK7XarsLBQkjRhwgTNnDnzrMevq6uT2+32uwEAACB2WTJ8Hz16VJIUHx8fsD4hIcGvXXP89a9/1ZtvvqkZM2boyiuvbFR/9dVXq6SkRE888YR69eql+Ph4ZWZmau3aterRo4cWLFig9evXn/E18vPzlZKS0nALtO4cAAAAscOS4TspKUmSdPz48YD1dXV1fu1C9dlnn+m2225TTk6Ofvvb3wZs43A4lJWV1ag8MTFRkydPliQtXbr0jK8zZ84c1dTUNNwqKiqa1V8AAABEB0tuNdi9e3dJ0oEDBwLWf/fdd5Kk1NTUkI9dXl6un/70pxo+fLhWrVolu90e8jF8e4L71ow3JSEhoWGWHgAAALBk+B40aJAkIygH4isfMGBASMfdvXu3rrzySjmdTi1fvjzgSZbBqK+vb9bzACCaeDxSUZFUVSWlpkpOp9SM+QwAiCmWXHYyevRotWvXTtu2bdORI0ca1RcXF0s6uU93MHbu3Kmf/OQnGjNmjFauXOkXvJ9++ulGS0K6d+8e8MI7krRnzx5J0kUXXRT06wNANHG5pN69pdGjpVtvNe579zbKAQBNs2T4TkpK0tSpU1VXV6fly5f71e3YsUMlJSUaNmyYnE5nQ3lhYaGysrJUUFDQ6HhbtmzRyJEjlZubqz/84Q+Ki/N/2w888EBDoPbZv3+/3nrrrUbHqq2t1UsvvSRJysvLa+5bBICI5XJJubnSvn3+5ZWVRjkBHACaZvNa9OoyR48e1VVXXaUdO3bo5Zdf1tixY1VWVqZf/OIX+vrrr/Xhhx+qb9++De3HjRunP//5z2rfvr0OHTrUUP7xxx9rzJgxqqur03XXXRfwtVavXq3CwkKNGjWqocxmsykuLk5z587VpEmTlJqaql27dumBBx7Qu+++q2nTpmnx4sUhvSe3262UlBTV1NQoOTk5tB8IAFiAx2PMcJ8evH1sNsnhkMrLWYICAIFYcs23ZMx+FxYWKj8/X7NmzdLNN9+szp07a+zYsXr88cflcDj82t90003atGmTJk2a5Fe+cuXKhhM0V69eHfTrl5WVafXq1XrnnXf0wgsvqLq6Wh07dtRll12mdevWhbTkBQCiRVFR08FbkrxeqaLCaHfKfAYA4L8sO/MdjZj5BhDpVq0y1nifzauvSrfc0vr9AYBIY8k13wAAawp2h9dm7AQLADGB8A0ACJrTaazpttkC19tsUlqa0Q4A0BjhGwAQNLtdWrTI+O/TA7jv8cKFnGwJAE0hfAMAQpKTI61ZI/Xs6V/ucBjlOTnh6RcARAJOuDQRJ1wCiCZc4RIAQmfZrQYBANZmt7OdIACEimUnAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJ2oS7AwCAyOTxSEVFUlWVlJoqOZ2S3R7uXgGAtRG+AQAhc7mk6dOlfftOljkc0qJFUk5O+PoFAFbHshMAQEhcLik31z94S1JlpVHucoWnXwAQCWxer9cb7k7ECrfbrZSUFNXU1Cg5OTnc3QGAkHk8Uu/ejYO3j81mzICXl7MEBQACYeYbABC0oqKmg7ckeb1SRYXRDgDQGOEbABC0qqqWbQcAsYbwDQAIWmpqy7YDgFhD+AYABM3pNNZ022yB6202KS3NaAcAaIzwDQAImt1ubCcoNQ7gvscLF3KyJQA0hfANAAhJTo60Zo3Us6d/ucNhlLPPNwA0ja0GTcRWgwCiCVe4BIDQcYVLAECz2O3SqFHh7gUARBaWnQAAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmaRPuDsQSr9crSXK73WHuCQAAAM6kQ4cOstlsLX5cwreJvv32W0lSWlpamHsCAACAM/n666/VrVu3Fj8u4dtEnTt3liR99dVXSklJCXNvEM3cbrfS0tJUUVGh5OTkcHcHUYzPGszCZw1m8X3W4uPjW+X4hG8TxcUZS+xTUlIYOGCK5ORkPmswBZ81mIXPGszSGktOJE64BAAAAExD+AYAAABMQvg2UUJCgn71q18pISEh3F1BlOOzBrPwWYNZ+KzBLK39WbN5ffvfAQAAAGhVzHwDAAAAJiF8AwAAACYhfAMAAAAmIXwDUaS2tla/+tWv1LdvXyUmJqpHjx6aPHmy9u3bF+6uAQAAccIlEDVqa2t11VVXqbi4WKmpqXI6ndq7d68++ugjdevWTSUlJerTp0+4uwkAQExj5huIEk8++aSKi4s1bNgwff7551q9erX++c9/av78+frmm280efLkcHcRAICYR/g2waZNmzRhwgRdeOGFstlsevTRR8PdJUSZEydOqKCgQJL07LPPqn379g11M2fO1IABA7Rp0yb961//ClcXEUEYsxBuH330kWw2W1C36667LtzdRQRYs2aNsrKy1LlzZyUmJiojI0MPP/ywjh8/bnpf2pj+ijHo8OHD+tGPfqRbb71VM2bMCHd3EIU++OADHTx4UH369NHgwYMb1efm5mrbtm166623NGTIkDD0EJGEMQvhduzYMf3P//yPXnnlFbVv31433HCDX/0f//hHHT9+XLfddhvhG0Hp3LmzZs2apYsvvljt2rXTp59+qrvuuktut1uLFy82tS+EbxNce+21uvbaayVJs2fPDnNvEI22bt0qSbr00ksD1vvKfe2AM2HMQriNHDlSXbp00SuvvKLLLrtMy5Yta6hzu91avny5+vTp41cOnMmVV17p97h3797atGmT/va3v5neF5adAFHgq6++kiQ5HI6A9b5yXzsAsLpt27ZJkgYMGOBXvmXLFnm9XmVmZoajW4gSZWVl+stf/qLRo0eb/tqEbyAKHD58WJKUlJQUsL5du3Z+7QDA6rZv3y5JjUL2li1bJDUO5UAw2rdvr4SEBP3oRz/SqFGjtGDBAtP7wLITIAr4dgy12WxnrAeASHGmme9A5UAwtmzZomPHjunjjz/W//t//08XXHCB5s2bZ2ofCN9AFOjQoYMk6ciRIwHrjx49Kkl+u6AAgJVt375dcXFx6t+/v1+5L3yz7ATNkZGRIcn4/MTFxWny5MmaNWtWwzfEZmDZiQWxxRJC1atXL0lq8kqWvnJfO6ClMF6hNRw8eFAVFRXq06eP33K677//Xjt37lRSUhIXDUOL8Hq9OnHihF9Za49rzHyb4PDhw/riiy8kScePH9d//vMfbdmyRfHx8frRj37UqD1bLCFUAwcOlCR98sknAet95XxNi2CEMmYxXqE1+NZ7nz5m7dq1S3V1dbr88ssVF8f8IYL3+OOP64orrtAPfvADeb1effTRR5o9e7YmTJigjh07+rVt9XHNi1ZXWFjoldToduGFFzb5nO3bt3sleUeNGuVXXlNT47XZbN6MjIxW7jUiSV1dnTclJcUryfvJJ580qh8wYIBXkvejjz4KQ+8QaUIdsxiv0NKeeeYZryTvvHnz/Mpff/11ryTvxIkTw9QzRKo5c+Z4+/bt623btq23Q4cO3szMTO9TTz3lPXr0aMD2rTmu8b+NJhg1apS8Xm+j2969e5t8DlssIRTx8fGaOnWqJGnq1Kl+a79/97vfadu2bfrxj3+syy+/PFxdRAQJdcxivEJLa2rm27dj09atW/XnP/85LFcnRGR68skn9dlnn+no0aNyu93atm2bZs+erbZt2wZs35rjGstOgpCbm6vS0tKQnrN8+XJlZWU1+zXZYgmhevjhh/W3v/1NxcXF+uEPfyin06kvv/xS//znP9WlSxctXbo03F1ElGK8QktrKnxfddVVGjBggMrKyjR79uyGi0EBLa01xzXCdxD27t2rzz77LKTn+HaXaC62WEKoEhMTVVhYqPz8fL366qtau3atOnXqpNtvv12PP/640tLSwt1FRCnGK7S0Dz/8MGB5WloaV+qFKVpzXLN5vWwAbEW9evVSZWWlDh065Hem96WXXqpPP/1Un3/+uX74wx+GsYcAYGC8AhBtWnNcI3y3sqYuehKI75/i4MGD6tSpk374wx/q888/b6j//vvv1b59e9ntdh06dIgzvQG0uFDHLMYrAFZntXGN0bCVBTppqambz9m2WLrkkkv4QwagVYQ6ZjFeAbA6q41rjIgW1NQ6o127dkmS+vXrZ3qfACAQxisA0aa1xzXCtwmWLl2q0aNHq1u3burQoYOGDBmiV155pcn2bLEEIJxCGbMYrwBEAiuNa4RvE/z973/X+PHj9fbbb+vTTz/VzTffrF/84hdavXp1wPZn22Jp165dmj17ts4777xW7zuA2BPKmMV4BSASWGlc44TLMLn22mvVtm1b/elPfwp3VwDgrBizAESbcI1rzHyHycGDB9W1a9dwdwMAgsKYBSDahGtcI3yHwcsvv6yPP/5Yd911V7i7AgBnxZgFINqEc1wjfJts3bp1uuuuu/T73/9el156abi7AwBnxJgFINqEe1wjfJvoj3/8o2666SY9//zzmjRpUri7AwBnxJgFINpYYVwjfJvkxRdfVF5enpYtW6a8vLxwdwcAzogxC0C0scq4xm4nJvjd736nWbNm6dlnn9WECRMayu12u7p16xbGngFAY4xZAKKNlcY1wrcJevfurS+//LJR+YUXXqi9e/ea3yEAOAPGLADRxkrjGuEbAAAAMAlrvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAIfv3v/+tX/7yl5oyZYrq6uo0bdo0nX/++Tr//PM1b948cfFkAAiM8A0ACMl7772n/v3767nnntPFF1+sRx55RM8884z69euno0eP6tFHH9X//d//hbubAGBJNi/TEwCAIO3Zs0cDBw5UXV2dXn31VV1//fXq1q2bOnTooL179+r111/XzTffrK5du+rrr7+WzWYLd5cBwFKY+QYABG327Nk6cuSIZs2apRtvvFH/+te/dPDgQY0cOVJxcXG68cYb1b59e1VXV+s///lPuLsLAJZD+AYABOXrr7/WunXr1KFDB82dO1eS9Omnn0qSBg0aJEmKi4tTly5dJEkejycs/QQAKyN8AwCC8v777+v777/XNddco8TEREnSli1bJEkDBw5saFddXS2bzaYLLrggHN0EAEsjfAMAgrJ161ZJ/kHbV+ab+d69e7eOHDmiSy65ROedd57pfQQAqyN8AwCCUl1dLUnq2LGjJKm+vl7bt29Xz5491bVrV0nSX/7yF0nStddeG5Y+AoDVEb4BAEHp0KGDJOmLL76QdHKW2zcTfuzYMS1cuFDx8fG66667wtZPALAywjcAICg//vGPJUkul0tut7thvfegQYN08OBB3XDDDSovL9fcuXP1gx/8IIw9BQDrYp9vAEBQPB6Phg0bps2bN6tXr17q1KmTtm7dqksvvVRffvmlvv32W02dOlUFBQXh7ioAWBbhGwAQtG+//VYPPPCA1q1bpwMHDkiS2rVrp8suu0wzZ87U+PHjw9xDALA2wjcAoFl8+3l/++23Ye4JAEQO1nwDAEJWXl6uAwcOaMiQIeHuCgBEFMI3ACBk//rXvySJ8A0AISJ8AwBC5gvfl156aZh7AgCRhTXfAAAAgEmY+QYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATPL/AevLoM6xqjsXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "zp = branchpoint\n",
    "dz = np.array([-0.5, 0.0])   # guide direction to find other branch\n",
    "for k in range(10):\n",
    "    ZS = root(sys, zp+dz, args=(zp,ds), jac=dsys, method='hybr') \n",
    "    zs =ZS.x\n",
    "    if d2Epot_dphi2(zs[0], zs[1]) > 0:\n",
    "        stable_phi.append( zs[0] )\n",
    "        stable_mu.append( zs[1] )\n",
    "    else:\n",
    "        nonstable_phi.append( zs[0] )\n",
    "        nonstable_mu.append( zs[1] )\n",
    "    dz = zs -zp\n",
    "    zp = zs\n",
    "\n",
    "plt.plot(stable_phi, stable_mu, 'bo')\n",
    "plt.plot(nonstable_phi, nonstable_mu, 'bo', fillstyle='none')\n",
    "plt.xlabel('$\\\\varphi$')\n",
    "plt.ylabel('$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = ['$-\\\\frac{1}{2}\\pi$', '$0$', '$\\\\frac{1}{2}\\pi$', '$\\pi$', '$\\\\frac{3}{2}\\pi$']\n",
    "plt.xticks(tick_pos, labels)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is much more to discover. In terms of physics it may be interesting to study the influence of an imperfection ($\\varphi_0\\ne \\pi/2$), whereas in terms of numerics it is desireable to automate the pseudo arc-length metric, to predict branched curves and to make sure all branches have been found. For further reading we recommend the textbooks _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_ by J. Guckenheimer and P. Holmes and _Elements of Applied Bifurcation Theory_ by Y. Kuznetsov."
   ]
  }
 ],
 "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.18"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}