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": 12,
   "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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi0 = 0.5*np.pi   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b2b7b2e459540be90bce9b60de2f920",
       "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": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "16d43d4068314b71a4e0917692a680f6",
       "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": 17,
   "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": 18,
   "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": 19,
   "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": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAInCAYAAABAyzWwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABS00lEQVR4nO3de3RU1f3+8WeYkCtkuAoxpCYQQA0h4aKiECEItVKgaUBtsQpFK1d/WFEq9AJaa9qCFiiowX4FqoJaiBQQq0VAU6KIcgkBwSBBQhokcpkEyAUm8/tjmsiQCcyE5MxM8n6tNWvI3nvO+Qx1HZ7u7LOPyW632wUAAACgwTXzdgEAAABAU0H4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4NpDdbldxcbHYWh0AAKBpInwbqKSkRBaLRSUlJd4uBQAAAF5A+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMEuDtAgAA/slmkzIzpcJCKSJCSkqSzGZvVwUAvo2ZbwCAxzIypNhYKTlZGjPG8R4b62gHANSO8A0A8EhGhjR6tBQfL338sVRS4niPj3e0E8ABoHYmu91u93YRTUVxcbEsFousVqvCw8O9XQ4AeMxmc8xwx8dLa9ZIzS6awqmslFJSpJwcKTeXJSgA4Aoz3wAAt2VmSocPS7NmOQdvyfHzzJlSXp5jHACgJsI3AMBthYWO9x49XPdXtVeNAwA4I3wDANwWEeF4z8lx3V/VXjUOAOCMNd8GYs03AH/Hmm8AuDrMfAMA3GY2S889J61f7wjaF+92kpLiaJ83j+ANALXhITsAAI+kpkqrVknTp0u33fZde0yMoz011Xu1AYCvY9mJgVh2AqAx4QmXAOA5Zr4BAHViNkuDBnm7CgDwL6z5BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAzi8+G7rKxMv/71rxUYGKhx48Z5/PnDhw/LZDJd8RUdHe30uS1btng0HgCaGptN2rJFWrnS8W6zebsiAPB9Ad4u4HI2b96sCRMm6Ntvv9X58+ev6lidO3dW8+bNXfbl5ubqhhtuqNEeEBCgLl26uPxMZGTkVdUDAP4sI0OaPl06fPi7tuho6bnnpNRUb1UFAL7PZ8P3ypUrNWnSJD377LMKCQnR+PHjr+p4H3zwgcvZ6l27dqlXr16aOnVqjb7IyEjt37//qs4LAI1NRoY0erQ0fLhj1rtHDyknR3r2WUf7qlUEcACojc8uO4mJidG+ffs0efJkmUymOh8nLCxM9957r8LCwlz2L1y4UF26dNFdd91V53MAQFNhszlmvIcPl9askfr1k1q0cLyvWeNof/xxlqAAQG18Nnz369dP11577VUfp3379nrjjTfUvn37Gn0nTpzQypUrNWXKFDVr5rN/FQDgMzIzHUtNZs2SLr1sNmsmzZwp5eU5xgEAamrSifPll1+W2WyudUlLaWmpZsyYoW7duikkJEStWrXSgAEDtGTJElVWVhpcLQB4X2Gh471HD9f9Ve1V4wAAzpps+LbZbHrhhRf0wAMPyGKxuBxz/PhxWa1WrVu3TlarVTt37lRcXJwmTJigESNGXPEm0PLychUXFzu9AMCfRUQ43nNyXPdXtVeNAwA4a7Lhe82aNcrPz3d5o6UkxcfHa8uWLUpPT1f37t0VGBiomJgYpaenKzk5WRs2bNDcuXMve460tDRZLJbqV1RUVEN8FQAwTFKSY1eTZ5+VLv0FYGWllJYmxcQ4xgEAamqy4XvhwoW64447dOONN7rsb9u2rQYOHOiyb8qUKZKkV1555bLnmDlzpqxWa/UrPz//6ooGAC8zmx3bCa5fL6WkSB9/LJWUON5TUhzt8+Y5xgEAavLZrQYbUnZ2tj766COtWbOmTp+PjY2VJH311VeqqKhQYGCgy3FBQUEKCgqqa5kA4JNSUx3bCU6fLt1223ftMTFsMwgAV9Ikw/fChQsVHR2tESNG1OnzF99seTXbIAKAv0pNlX70I8euJoWFjjXeSUnMeAPAlTS5ZScnT57UihUrNHny5MtuL3jTTTfppZdectmXm5srSerSpUutT80EgMbObJYGDZJ++lPHO8EbAK6s0YTvPXv2qH///po1a9Zlx7388suSpAcffPCy44qKimpdlrJo0SJJ0rhx4zyuEwAAAE1Xownf6enpysrKUlpamk6cOOFyTNX2gvfdd5/atGlz2eOZTCa99957mjBhgg4cOKDy8nIdOnRIv/jFL5SZmak777xTM2bMaIivAgAAgEbKp8O3yWSSyWTSz3/+c0nS8uXLq9uWLVvmNDYlJUWtW7dWSkpKrcF67dq1OnLkiB555JErnvujjz7SvHnztH//fvXv319hYWHq06ePvvzySy1ZskQbNmyo9UZLAAAAwBWT3W63e7uIpqK4uFgWi0VWq1Xh4eHeLgcAAAAG8+mZbwAAAKAxIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABgnwdgEAAP9ks0mZmVJhoRQRISUlSWazt6sCAN/GzDcAwGMZGVJsrJScLI0Z43iPjXW0AwBqR/gGAHgkI0MaPVqKj5c+/lgqKXG8x8c72gngAFA7k91ut3u7iKaiuLhYFotFVqtV4eHh3i4HADxmszlmuOPjpTVrpGYXTeFUVkopKVJOjpSbyxIUAHCFmW8AgNsyM6XDh6VZs5yDt+T4eeZMKS/PMQ4AUBPhGwDgtsJCx3uPHq77q9qrxgEAnBG+AQBui4hwvOfkuO6vaq8aBwBwxppvA7HmG4C/Y803AFwdZr4BAG4zm6XnnpPWr3cE7Yt3O0lJcbTPm0fwBoDa8JAdAIBHUlOlVauk6dOl2277rj0mxtGemuq92gDA17HsxEAsOwHQmPCESwDwHDPfAIA6MZulQYO8XQUA+BfWfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYJ8HYBAAD/ZLNJmZlSYaEUESElJUlms7erAgDfxsw3AMBjGRlSbKyUnCyNGeN4j411tAMAakf4BgB4JCNDGj1aio+XPv5YKilxvMfHO9oJ4ABQO5Pdbrd7u4imori4WBaLRVarVeHh4d4uBwA8ZrM5Zrjj46U1a6RmF03hVFZKKSlSTo6Um8sSFABwhZlvAIDbMjOlw4elWbOcg7fk+HnmTCkvzzEOAFAT4RsA4LbCQsd7jx6u+6vaq8YBAJwRvgEAbouIcLzn5Ljur2qvGgcAcObz4busrEy//vWvFRgYqHHjxtXpGHPmzJHJZKr1dbnjvvvuuxo8eLAsFovCw8M1YMAArVq1qm5fBgD8XFKSFB0tPfusY433xSorpbQ0KSbGMQ4AUJNPh+/NmzerZ8+eevHFF3X+/PmrOlZ4eLi6d+/u8hVRyxTN/PnzNWzYMHXt2lVffvml8vLydPPNN+vuu+/WnDlzrqoeAPBHZrP03HPS+vWOmysv3u0kJcXRPm8eN1sCQG18dreTlStXatKkSXr22WcVEhKi8ePHa+zYsVq2bJnHx5ozZ44OHz7s0WdzcnKUmJio3r17a9u2bTKZTNV9ycnJ2rJliz7++GP169fP7WOy2wmAxiIjQ5o+3XHzZZWYGEfwTk31WlkA4PN8duY7JiZG+/bt0+TJk52Cr1Gef/552Ww2/eIXv6hx/gkTJkiS5s6da3hdAOALUlOlgwelzZulFSsc77m5BG8AuBKffby8JzPKDWHdunWSpP79+9foGzBggCRpw4YNunDhggICfPavEQAajNksDRrk7SoAwL/47Mx3fcvLy9OYMWMUFRWlwMBAdejQQSNHjtQHH3xQY2xBQYG+/fZbSY4Z+EtFRkYqMDBQZWVl+vLLLxu8dgAAADQOTSZ8f/rpp+rbt68+/fRTnTlzRu+8845OnTqlIUOG6JlnnnEae+zYMUlScHCwQkJCahzLZDKpVatWTmNdKS8vV3FxsdMLAAAATVeTCN9jxozR9u3b9dhjjykiIkKBgYHq27ev1q9fr9atW+u3v/2ttm7dWj3+3LlzkqTAwMBajxkUFOQ01pW0tDRZLJbqV1RUVD19IwAAAPijJhG+u3Xrph4uHsdmsVh03333SZKWLl1a3R4aGipJqqioqPWY5eXlTmNdmTlzpqxWa/UrPz+/TvUDAACgcWjydwrGxsZKkr744ovqto4dO0pyPOCntLS0xtITu92u06dPS1Kte4RLjtnxqhlyAAAAoEnMfF9O5f8e0XbxdoKRkZFq166dJMeNmpcqKChQRUWFQkJC1LVrV2MKBQAAgN9r9OE7Pz9fHTt2VHZ2tsv+3NxcSdL111/v1D5ixAhJUlZWVo3PVK0Pv+uuu9hmEAAAAG5rNOF7z5496t+/v2bNmuXUbrPZ9M0332jDhg01PnPq1CmtWLFCkjRu3DinvunTp8tsNutvf/ubLn0I6JIlSyRJTzzxRD1+AwAAADR2jSZ8p6enKysrS2lpaTpx4kR1e9Vykqeeekrz58/X0aNHVVFRoe3bt2v48OGyWq2aM2dO9YNzqsTFxWnevHnatm2bJk+erOPHj+vkyZN6/PHHtWnTJs2ePdvrDwICAACAfzHZL53W9SGXe6z80qVLnWarN27cqHvuuUcDBw5URkaG02c///xzvfXWW3r//fd18OBBlZWVqV27durXr5+mTJmiIUOG1HqeDRs2aO7cufr8889lt9vVs2dP/fKXv9To0aM9/j7FxcWyWCyyWq0KDw/3+PMAAADwbz4dvhsbwjcAAEDT1miWnQAAAAC+jq06AAB1YrNJmZlSYaEUESElJUlms7erAgDfxsw3AMBjGRlSbKyUnCyNGeN4j411tAMAakf4BgB4JCNDGj1aio+XPv5YKilxvMfHO9oJ4ABQO264NBA3XALwdzabY4Y7Pl5as0ZqdtEUTmWllJIi5eRIubksQQEAV5j5BgC4LTNTOnxYmjXLOXhLjp9nzpTy8hzjAAA1Eb4BAG4rLHS89+jhur+qvWocAMAZ4RsA4LaICMd7To7r/qr2qnEAAGes+TYQa74B+DvWfAPA1WHmGwDgNrNZeu45af16R9C+eLeTlBRH+7x5BG8AqA0P2QEAeCQ1VVq1Spo+Xbrttu/aY2Ic7amp3qsNAHwdy04MxLITAI0JT7gEAM8x8w0AqBOzWRo0yNtVAIB/Yc03AAAAYBDCNwAAAGAQwjcAAABgEMI3AAAAYBDCNwAAAGAQwjcAAABgEMI3AAAAYBDCNwAAAGAQwjcAAABgEMI3AAAAYBDCNwAAAGAQwjcAAABgEMI3AAAAYJAAbxcAAPBPNpuUmSkVFkoREVJSkmQ2e7sqAPBtzHwDADyWkSHFxkrJydKYMY732FhHOwCgdoRvAIBHMjKk0aOl+Hjp44+lkhLHe3y8o50ADgC1M9ntdru3i2gqiouLZbFYZLVaFR4e7u1yAMBjNptjhjs+XlqzRmp20RROZaWUkiLl5Ei5uSxBAQBXmPkGALgtM1M6fFiaNcs5eEuOn2fOlPLyHOMAADURvgEAbissdLz36OG6v6q9ahwAwBnhGwDgtogIx3tOjuv+qvaqcQAAZ6z5NhBrvgH4O9Z8A8DVYeYbAOA2s1l67jlp/XpH0L54t5OUFEf7vHkEbwCoDQ/ZAQB4JDVVWrVKmj5duu2279pjYhztqaneqw0AfB3LTgzEshMAjQlPuAQAzzHzDQCoE7NZGjTI21UAgH9hzTcAAABgEMI3AAAAYBCfD99lZWX69a9/rcDAQI0bN65Ox9i+fbseeeQR9ejRQ6GhoQoODlZsbKwmTZqkvLw8l5/ZsmWLTCZTra/o6Oi6fykAAAA0ST695nvz5s2aMGGCvv32W50/f75Ox3jnnXc0fPhwXXfddVqwYIFuv/12nT9/XmvXrtVjjz2m1157Te+9955uu/iW/f8JCAhQly5dXB43MjKyTvUAAACg6fLZ8L1y5UpNmjRJzz77rEJCQjR+/Pg6Hae0tFSS9Oabb+qWW26pbn/ooYd04cIFTZo0SQ8++KC++OKLGp+NjIzU/v376/YFAAAAgEv47LKTmJgY7du3T5MnT5bJZKrzcdq0aaOhQ4c6Be8qI0eOlCTt379fR44cqfM5AAAAAHf47Mx3v3796uU4gwcP1uDBg132WSyWejkHAAAA4A6fnfk2QtWSks6dO+t73/tejf7S0lLNmDFD3bp1U0hIiFq1aqUBAwZoyZIlqqysvOLxy8vLVVxc7PQCAABA09Wkw/ff//53SdLvfvc7l/3Hjx+X1WrVunXrZLVatXPnTsXFxWnChAkaMWLEFW8CTUtLk8ViqX5FRUXV+3cAAACA/2iy4Ts7O1svvviiRo0apbFjx9boj4+P15YtW5Senq7u3bsrMDBQMTExSk9PV3JysjZs2KC5c+de9hwzZ86U1WqtfuXn5zfU1wEAAIAfaJLhu6ioSHfffbduvvnm6tnvS7Vt21YDBw502TdlyhRJ0iuvvHLZ8wQFBSk8PNzpBQAAgKaryYXvU6dO6Qc/+IHat2+vd955R6GhoR4fIzY2VpL01VdfqaKior5LBAAAQCPls7udNISioiINHTpU7dq109q1a+sUvCU53Wx5NdsgAoA/s9mkzEypsFCKiJCSkiSz2dtVAYBvazIz3//97381cOBARUZGav369U7B++WXX9bevXudxt9000166aWXXB4rNzdXktSlSxc1b9684YoGAB+VkSHFxkrJydKYMY732FhHOwCgdo0mfO/Zs0f9+/fXrFmzavQdOXJEt99+u2644Qa9/fbbCg4Odur/wx/+oO3btzu1FRUVac2aNS7PtWjRIknSuHHj6qV2APAnGRnS6NFSfLz08cdSSYnjPT7e0U4AB4DaNZplJ+np6crKylJWVpamT5+utm3bSpIOHTqk5ORkFRQUKDExUQ888ECNzxYVFdVoM5lMeu+99zRhwgQ99thjio6OVkFBgdLS0pSZmak777xTM2bMaPDvBQC+xGaTpk+Xhg+X1qyRmv1vCqdfP8fPKSnS449LP/oRS1AAwBWfDt+Xrqdevny5li9fLklaunSp08xzSkqKVqxYoYEDB6pNmzbV7WvXrq1+dPzq1avdPvdHH32kt956S2vXrlX//v11+vRptWzZUj179tSSJUv04IMPqlmzRvOLAwBwS2amdPiwtHLld8G7SrNm0syZ0m23OcYNGuSNCgHAt5nsdrvd20U0FcXFxbJYLLJarWw7CMAvrVzpWONdUiK1aFGzv6RECg+XVqyQfvpT4+sDAF/H1C0AwG0REY73nBzX/VXtVeMAAM6Y+TYQM98A/J3N5tjVJD7eec23JFVWOtZ85+RIubms+QYAV5j5BgC4zWyWnntOWr/eEbQv3u0kJcXRPm8ewRsAauPTN1wCAHxPaqq0apVj15PbbvuuPSbG0Z6a6r3aAMDXsezEQCw7AdCY8IRLAPAcM98AgDoxm9lOEAA8xZpvAAAAwCCEbwAAAMAghG8AAADAIIRvAAAAwCCEbwAAAMAghG8AAADAIIRvAAAAwCCEbwAAAMAghG8AAADAIIRvAAAAwCCEbwAAAMAghG8AAADAIAHeLgAA4J9sNikzUyoslCIipKQkyWz2dlUA4NuY+QYAeCwjQ4qNlZKTpTFjHO+xsY52AEDtCN8AAI9kZEijR0vx8dLHH0slJY73+HhHOwEcAGpnstvtdm8X0VQUFxfLYrHIarUqPDzc2+UAgMdsNscMd3y8tGaN1OyiKZzKSiklRcrJkXJzWYICAK4w8w0AcFtmpnT4sDRrlnPwlhw/z5wp5eU5xgEAaiJ8AwDcVljoeO/Rw3V/VXvVOACAM8I3AMBtERGO95wc1/1V7VXjAADOWPNtINZ8A/B3rPkGgKvDzDcAwG1ms/Tcc9L69Y6gffFuJykpjvZ58wjeAFAbHrIDAPBIaqq0apU0fbp0223ftcfEONpTU71XGwD4OpadGIhlJwAaE55wCQCeY+YbAFAnZrM0aJC3qwAA/8KabwAAAMAghG8AAADAIIRvAAAAwCCEbwAAAMAghG8AAADAIIRvAAAAwCCEbwAAAMAghG8AAADAIIRvAAAAwCA+H77Lysr061//WoGBgRo3blydj1NZWamFCxcqISFBoaGhateunX784x9r165dl/3cJ598ouHDh6tt27YKCwtTnz59tGTJEtnt9jrXAgAAgKbJp8P35s2b1bNnT7344os6f/58nY9TWVmpe+65R9OnT9fEiRNVVFSkTz75RCUlJbrlllv0/vvvu/xcRkaGBgwYoICAAO3YsUPHjh3Tvffeq4kTJ2r8+PF1rgcAAABNk8nuo1O4K1eu1KRJk/Tss88qJCRE48eP19ixY7Vs2TKPj/Xiiy9q8uTJevzxxzV37tzq9lOnTqlLly5q3ry5Dh48qJYtW1b3HT9+XLGxsWrVqpUOHDigkJCQ6r5x48Zp+fLleuONN3Tvvfe6XUdxcbEsFousVqvCw8M9/h4AAADwbz478x0TE6N9+/Zp8uTJMplMV3WsefPmSZIefvhhp/bWrVvrnnvu0fHjx2uE+pdeekklJSW6//77nYK3JE2cOFGS9Oc///mq6gIAAEDT4rPhu1+/frr22muv+jg5OTk6dOiQ2rdvr65du9boHzBggCRp7dq1Tu1VP/fv37/GZ/r27avg4GDt2LFDBQUFV10jAPgjm03askVaudLxbrN5uyIA8H0+G77ry+7duyU5ZtJd6dy5s9M4SbLZbNq7d2+tnwsICFBUVFSNz12qvLxcxcXFTi8AaAwyMqToaCk5WRozxvEeHe1oBwDUrtGH72PHjkmS2rRp47K/devWkqSioiJVVlZKkk6fPq2ysjK3Pld1fFfS0tJksViqX1WBHQD8WUaGNHq0dPSoc3tBgaOdAA4AtWv04fvcuXOSpMDAQJf9QUFBNcZWvbvzuYvHXmrmzJmyWq3Vr/z8fM+KBwAfY7NJ06ZJrm7Vr2p79FGWoABAbQK8XUBDCw0NlSRVVFS47C8vL68xturdnc9dPPZSQUFBTuEeAPxdZmbNGe+L2e1Sfr5j3KBBhpUFAH6j0c98d+zYUZJ08uRJl/2nTp2SJF1zzTVq1szx19GqVSsFBwe79bmIiIh6rRcAfFlhYf2OA4CmptGH78TERElSXl6ey/5Dhw5Jknr27FndZjab1aNHj1o/d+HCheolJBd/DgAaO3fnG5iXAADXGn34jouLU+fOnVVUVKSDBw/W6N+6daskaeTIkU7tI0aMkCRlZWXV+Mznn3+usrIy9e7dW5GRkQ1QNQD4pqQkqVMnqbbHL5hMUlSUYxwAoKZGE7737Nmj/v37a9asWTX6nnjiCUnSkiVLnNpPnz6tt956S9dcc43GjRvn1Ddx4kS1bNlSr776avXOJ1XS09MlSTNmzKjHbwAAvs9slhYscPz50gBe9fP8+Y5xAICaGk34Tk9PV1ZWltLS0nTixAmnvocfflipqan6y1/+ovT0dJ09e1ZfffWV7rnnHpWUlOjVV191erS85FgDvnTpUhUUFGjMmDE6cuSISkpKNG/ePC1btkxjx4716NHyANBYpKZKq1ZJl/7ir1MnR3tqqnfqAgB/YLLbXW0Y5Rsu91j5pUuXOs1Wb9y4Uffcc48GDhyojIyMGp+12WxatGiR/u///k+5ubkKDQ1VUlKSZs+erV69etV6nk8++US///3vlZWVpfLycl1//fWaMGGCHn74YY8fe19cXCyLxSKr1arw8HCPPgsAvsZmc+xqUljoWOOdlMSMNwBciU+H78aG8A0AANC0NZplJwAAAICvI3wDAAAABvEofG/evLn64TIAAAAAPOPRmu9mzZrJZDKpU6dOSkxMdHrFxMRUj3vwwQfVp08fTZ48uUGK9les+QYAAGjaPArfDz30kHbt2qWcnBxVVFQ4DvC/HT/Cw8PVs2dPdevWTW+//bYCAgJ07NixhqnaTxG+AQAAmrY67XZy4cIF7du3T7t27dLOnTu1c+dOffLJJ6qoqJDJZJLdbtf3vvc9HT58uAFK9l+EbwAAgKat3rYaPHPmjJYtW6Ynn3xSPXr00GuvvabY2Nj6OHSjQfgGAABo2uptt5MWLVpo6tSpeuutt7R9+3bt2LGjvg4NAAAANAoN8pCduLg4BQYGaufOnfV9aL/GzDcAAEDT1iD7fMfExOjAgQMNcWgAAADAbwV4MnjatGlKTExUr169FBcXp+bNm7scd/DgQXXs2LFeCgQAAAAaC4/C91//+tfqrQWbN2+u66+/Xr169VJiYqLi4+MVGhqq1157Tbm5uXrppZcapGAAAADAX3m05nv9+vXVWwvu3LlTX3/9teMg/wvkVTp37qyxY8eqd+/e6tOnjzp06FC/Vfsp1nwDaExsNikzUyoslCIipKQkyWz2dlUA4Nuu6obLU6dOOYXxHTt26Msvv1RlZaXj4P8L5R07dlSfPn20du3a+qnaTxG+ATQWGRnStGnS0aPftXXqJC1YIKWmeq8uAPB19b7byblz55SdnV0dxnfu3KmcnBydP39eNputPk/ldwjfABqDjAxp9Gjp0n89qn4JumoVARwAatMgWw1equqJmD179mzoU/k0wjcAf2ezSdHRzjPeFzOZHDPgeXksQQEAVxpkq8FLBQQENPngDQCNQWZm7cFbcsyG5+c7xgEAajIkfAMAGofCwvodBwBNDeEbAOC2iIj6HQcATQ3hGwDgtqQkx5ruS3aYrWYySVFRjnEAgJoI3wAAt5nNju0EpZoBvOrn+fO52RIAakP4BgB4JDXVsZ1gZKRze6dObDMIAFdiyFaDcGCrQQCNCU+4BADPBXi7AACAfzKbpUGDvF0FAPgXlp0AAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGCfB2AQAA/2SzSZmZUmGhFBEhJSVJZrO3qwIA30b4BgB4LCNDmjZNOnr0u7ZOnaQFC6TUVO/VBQC+jmUnAACPZGRIo0c7B29JKihwtGdkeKcuAPAHPh2+y8rK9PTTT6t79+4KDg5WRESE7r//fuXl5bl9jMOHD8tkMl3xFR0d7fS5LVu2eDQeAJoCm80x42231+yranv0Ucc4AEBNPrvspLS0VEOGDFF2drZeeeUV/fCHP1Rubq7Gjh2rXr166cMPP1RCQoLbx+vcubOaN2/usi83N1c33HBDjfaAgAB16dLF5WciIyPdPjcANBaZmTVnvC9mt0v5+Y5xgwYZVhYA+A2fDd9z5sxRVlaWFi1apLvvvluSlJCQoNWrV6t79+667777tHv3bpndvLvngw8+cDlbvWvXLvXq1UtTp06t0RcZGan9+/df1fcAgMaksLB+xwFAU+OTy07OnTunxYsXKzAwUGPHjnXq69Kli4YMGaK9e/dqw4YNVzxWWFiY7r33XoWFhbnsX7hwobp06aK77rqrXmoHgMYsIqJ+xwFAU+OT4XvTpk06e/asEhIS1KJFixr9AwYMkCStXbv2isdq37693njjDbVv375G34kTJ7Ry5UpNmTJFzZr55F8FAPiUpCTHriYmk+t+k0mKinKMAwDU5JOJc/fu3ZKkmJgYl/2dO3d2GldXL7/8ssxms8aPH++yv7S0VDNmzFC3bt0UEhKiVq1aacCAAVqyZIkqKyuv6twA4I/MZsd2glLNAF718/z57PcNALXxyfB97NgxSVKbNm1c9rdu3dppXF3YbDa98MILeuCBB2SxWFyOOX78uKxWq9atWyer1aqdO3cqLi5OEyZM0IgRI3T+/PnLnqO8vFzFxcVOLwDwd6mp0qpV0qX3nXfq5Ghnn28AqJ1P3nB57tw5SVJgYKDL/qCgIKdxdbFmzRrl5+e7vNFSkuLj47VlyxYNHDiwui0mJkbp6enKzc3Vhg0bNHfuXM2aNavWc6Slpempp56qc40A4KtSU6Uf/YgnXAKAp3xy5js0NFSSVFFR4bK/vLzcaVxdLFy4UHfccYduvPFGl/1t27Z1Ct4XmzJliiTplVdeuew5Zs6cKavVWv3Kz8+vc70A4GvMZsd2gj/9qeOd4A0AV+aTM98dO3aUJJ08edJl/6lTpyRJEXW8nT47O1sfffSR1qxZU6fPx8bGSpK++uorVVRUXHaGvmqWHgAAAPDJme/ExERJqvVJlocOHZIk9ezZs07HX7hwoaKjozVixIg6ff7imy1Ntd3yDwAAAFzCJ8N3cnKywsLClJ2drbNnz9bo37p1qyRp5MiRHh/75MmTWrFihSZPnnzZ7QVvuukmvfTSSy77cnNzJTn2HK/tqZkAAADApXwyfIeGhmrq1KkqLy/X8uXLnfoOHTqkjRs3Ki4uTsOGDatu37Nnj/r373/ZGyAlx/aCkvTggw9edlxRUVGty1IWLVokSRo3btwVvgkAAADwHZ8M35I0e/Zs3XrrrXryySe1evVqlZaWKjs7W6NGjVJoaKhef/11p0fLp6enKysrS2lpaTpx4oTLY1ZtL3jffffVuo1hFZPJpPfee08TJkzQgQMHVF5erkOHDukXv/iFMjMzdeedd2rGjBn1+p0BAADQuPls+A4JCdGmTZs0ffp0zZw5U61bt9bQoUMVFxennTt3KiEhwWl8SkqKWrdurZSUlFqD9dq1a3XkyBE98sgjVzz/Rx99pHnz5mn//v3q37+/wsLC1KdPH3355ZdasmSJNmzYUOuNlgAAAIArJrvdbvd2EU1FcXGxLBaLrFarwsPDvV0OAAAADOazM98AAABAY0P4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMEuDtAgAA/slmkzIzpcJCKSJCSkqSLnr2GQDABcI3AMBjGRnStGnS0aPftXXqJC1YIKWmeq8uAPB1LDsBAHgkI0MaPdo5eEtSQYGjPSPDO3UBgD/gCZcG4gmXAPydzSZFR9cM3lVMJscMeF4eS1AAwBVmvgEAbsvMrD14S5LdLuXnO8YBAGoifAMA3FZYWL/jAKCpIXwDANwWEVG/4wCgqSF8AwDclpTkWNNtMrnuN5mkqCjHOABATYRvAIDbzGbHdoJSzQBe9fP8+dxsCQC1IXwDADySmiqtWiVFRjq3d+rkaGefbwCoHVsNGoitBgE0JjzhEgA8xxMuAQB1YjZLgwZ5uwoA8C8sOwEAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAxC+AYAAAAMQvgGAAAADEL4BgAAAAwS4O0CAAD+yWaTMjOlwkIpIkJKSpLMZm9XBQC+jfANAPBYRoY0bZp09Oh3bZ06SQsWSKmp3qsLAHwdy04AAB7JyJBGj3YO3pJUUOBoz8jwTl0A4A9Mdrvd7u0imori4mJZLBZZrVaFh4d7uxwA8JjNJkVH1wzeVUwmxwx4Xh5LUADAFWa+AQBuy8ysPXhLkt0u5ec7xgEAaiJ8AwDcVlhYv+MAoKkhfAMA3BYRUb/jAKCpIXwDANyWlORY020yue43maSoKMc4AEBNhG8AgNvMZsd2glLNAF718/z53GwJALXx6fBdVlamp59+Wt27d1dwcLAiIiJ0//33Ky8vz6PjzJkzRyaTqdbXuHHjav3su+++q8GDB8tisSg8PFwDBgzQqlWrrvKbAYD/Sk2VVq2SIiOd2zt1crSzzzcA1M5nw3dpaanuuOMOzZ07V88884xOnjypf/3rX9qzZ4969eql3bt3e3S88PBwde/e3eUropbFifPnz9ewYcPUtWtXffnll8rLy9PNN9+su+++W3PmzKmHbwkA/ik1VTp8WNq8WVqxwvGel0fwBoAr8dl9vn/1q1/pz3/+sxYtWqQpU6ZUt3/11Vfq3r27rr/+eu3evVtmN363OWfOHB0+fFjLli1z+/w5OTlKTExU7969tW3bNpku+v1qcnKytmzZoo8//lj9+vVz+5js8w0AANC0+eTM97lz57R48WIFBgZq7NixTn1dunTRkCFDtHfvXm3YsKHBanj++edls9n0i1/8wil4S9KECRMkSXPnzm2w8wMAAKDx8cnwvWnTJp09e1YJCQlq0aJFjf4BAwZIktauXdtgNaxbt06S1L9//1rPv2HDBl24cKHBagAAAEDj4pPhu2o9d0xMjMv+zp07O41zR15ensaMGaOoqCgFBgaqQ4cOGjlypD744IMaYwsKCvTtt9/WWkNkZKQCAwNVVlamL7/80u0aAAAA0LT5ZPg+duyYJKlNmzYu+1u3bu00zh2ffvqp+vbtq08//VRnzpzRO++8o1OnTmnIkCF65plnXJ4/ODhYISEhNY5lMpnUqlWrK9ZQXl6u4uJipxcAAACaLp8M3+fOnZMkBQYGuuwPCgpyGnclY8aM0fbt2/XYY48pIiJCgYGB6tu3r9avX6/WrVvrt7/9rbZu3er2+d2tIS0tTRaLpfoVFRXlVr0AAABonHwyfIeGhkqSKioqXPaXl5c7jbuSbt26qUePHjXaLRaL7rvvPknS0qVL3T6/uzXMnDlTVqu1+pWfn+9WvQAAAGicArxdgCsdO3aUJJ08edJl/6lTpySp1v25PREbGytJ+uKLL2qcv6ysTKWlpTWWntjtdp0+ffqKNQQFBVXPkAMAAAA+OfOdmJgoSbU+yfLQoUOSpJ49e171uSorKyXJaTvByMhItWvXrtYaCgoKVFFRoZCQEHXt2vWqawAAAEDT4JPhOzk5WWFhYcrOztbZs2dr9Fetzx45cuQVj5Wfn6+OHTsqOzvbZX9ubq4k6frrr3dqHzFihCQpKyur1vPfddddCgjwyV8eAAAAwAf5ZPgODQ3V1KlTVV5eruXLlzv1HTp0SBs3blRcXJyGDRtW3b5nzx71799fs2bNchpvs9n0zTffuHwgz6lTp7RixQpJ0rhx45z6pk+fLrPZrL/97W+69CGgS5YskSQ98cQTdf6OAAAAaHp8MnxL0uzZs3XrrbfqySef1OrVq1VaWqrs7GyNGjVKoaGhev31150eLZ+enq6srCylpaXpxIkT1e1Vy0meeuopzZ8/X0ePHlVFRYW2b9+u4cOHy2q1as6cOdUPzqkSFxenefPmadu2bZo8ebKOHz+ukydP6vHHH9emTZs0e/Zsjx4tDwAAAJjsl07r+pCysjL96U9/0uuvv64jR47IYrFo6NChevrpp6sftFNl48aNuueeezRw4EBlZGQ4reH+/PPP9dZbb+n999/XwYMHVVZWpnbt2qlfv36aMmWKhgwZUmsNGzZs0Ny5c/X555/LbrerZ8+e+uUvf6nRo0d7/H2Ki4tlsVhktVoVHh7u8ecBAADg33w6fDc2hG8AjYnNJmVmSoWFUkSElJQkXfQLSQCAC9wtCADwWEaGNG2adPTod22dOkkLFkipqd6rCwB8nc+u+QYA+KaMDGn0aOfgLUkFBY72jAzv1AUA/oBlJwZi2QkAf2ezSdHRNYN3FZPJMQOel8cSFABwhZlvAIDbMjNrD96SZLdL+fmOcQCAmgjfAAC3FRbW7zgAaGoI3wAAt0VE1O84AGhqCN8AALclJTnWdF/0KAUnJpMUFeUYBwCoifANAHCb2ezYTlCqGcCrfp4/n5stAaA2hG8AgEdSU6VVq6TISOf2Tp0c7ezzDQC1Y6tBA7HVIIDGhCdcAoDneMIlAKBOzGZp0CBvVwEA/oVlJwAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBArxdAADAP9lsUmamVFgoRURISUmS2eztqgDAtxG+AQAey8iQpk2Tjh79rq1TJ2nBAik11Xt1AYCvY9kJAMAjGRnS6NHOwVuSCgoc7RkZ3qkLAPyByW63271dRFNRXFwsi8Uiq9Wq8PBwb5cDAB6z2aTo6JrBu4rJ5JgBz8tjCQoAuMLMNwDAbZmZtQdvSbLbpfx8xzgAQE2EbwCA2woL63ccADQ1hG8AgNsiIup3HAA0NYRvAIDbkpIca7pNJtf9JpMUFeUYBwCoifANAHCb2ezYTlCqGcCrfp4/n5stAaA2hG8AgEdSU6VVq6TISOf2Tp0c7ezzDQC1Y6tBA7HVIIDGhCdcAoDneMIlAKBOzGZp0CBvVwEA/oVlJwAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBfDp8l5WV6emnn1b37t0VHBysiIgI3X///crLy/PoONu3b9cjjzyiHj16KDQ0VMHBwYqNjdWkSZNqPdaWLVtkMplqfUVHR9fDNwQAAEBT4rO7nZSWlmrIkCHKzs7WK6+8oh/+8IfKzc3V2LFj1atXL3344YdKSEi44nHeeecdDR8+XNddd50WLFig22+/XefPn9fatWv12GOP6bXXXtN7772n2267rcZnAwIC1KVLF5fHjbx0g1sAAADgCnw2fM+ZM0dZWVlatGiR7r77bklSQkKCVq9ere7du+u+++7T7t27Zb7CprKlpaWSpDfffFO33HJLdftDDz2kCxcuaNKkSXrwwQf1xRdf1PhsZGSk9u/fX4/fCgAAAE2ZTy47OXfunBYvXqzAwECNHTvWqa9Lly4aMmSI9u7dqw0bNlzxWG3atNHQoUOdgneVkSNHSpL279+vI0eO1E/xAAAAQC18Mnxv2rRJZ8+eVUJCglq0aFGjf8CAAZKktWvXXvFYgwcP1vvvv++yz2KxXF2hAAAAgAd8Mnzv3r1bkhQTE+Oyv3Pnzk7j6qpqSUnnzp31ve99r0Z/aWmpZsyYoW7duikkJEStWrXSgAEDtGTJElVWVl7x+OXl5SouLnZ6AQAAoOnyyfB97NgxSY4lI660bt3aaVxd/f3vf5ck/e53v3PZf/z4cVmtVq1bt05Wq1U7d+5UXFycJkyYoBEjRuj8+fOXPX5aWposFkv1Kyoq6qrqBQAAgH/zyfB97tw5SVJgYKDL/qCgIKdxdZGdna0XX3xRo0aNqrGuXJLi4+O1ZcsWpaenq3v37goMDFRMTIzS09OVnJysDRs2aO7cuZc9x8yZM2W1Wqtf+fn5da4XAAAA/s8nw3doaKgkqaKiwmV/eXm50zhPFRUV6e6779bNN99cPft9qbZt22rgwIEu+6ZMmSJJeuWVVy57nqCgIIWHhzu9AKCxsNmkLVuklSsd7zabtysCAN/nk+G7Y8eOkqSTJ0+67D916pQkKSIiwuNjnzp1Sj/4wQ/Uvn17vfPOO3UK8LGxsZKkr776qtb/gwAAjVlGhhQdLSUnS2PGON6jox3tAIDa+WT4TkxMlKRanz556NAhSVLPnj09Om5RUZGSk5PVunVrvf/++3Xe7eTimy1NJlOdjgEA/iojQxo9Wjp61Lm9oMDRTgAHgNr5ZPhOTk5WWFiYsrOzdfbs2Rr9W7dulfTdPt3u+O9//6uBAwcqMjJS69evd5rxfvnll7V3716n8TfddJNeeukll8fKzc2V5NhzvHnz5m7XAAD+zmaTpk2T7PaafVVtjz7KEhQAqI1Phu/Q0FBNnTpV5eXlWr58uVPfoUOHtHHjRsXFxWnYsGHV7Xv27FH//v01a9asGsc7cuSIbr/9dt1www16++23FRwc7NT/hz/8Qdu3b3dqKyoq0po1a1zWt2jRIknSuHHj6vDtAMB/ZWbWnPG+mN0u5ec7xgEAavLZx8vPnj1bH330kZ588kl16NBBw4YNU25ursaOHavQ0FC9/vrrTo+WT09PV1ZWlrKysjR9+nS1bdtWkiOsJycnq6CgQImJiXrggQdqnKuoqKhGm8lk0nvvvacJEyboscceU3R0tAoKCpSWlqbMzEzdeeedmjFjRsP9BQCADyosrN9xANDU+Gz4DgkJ0aZNm/SnP/1JM2fO1H333SeLxaKhQ4dq9erV1Q/aqZKSkqIVK1Zo4MCBTvuDr127tvrR8atXr3b7/B999JHeeustrV27Vv3799fp06fVsmVL9ezZU0uWLNGDDz6oZs188hcHANBg3L3PvQ73wwNAk2Cy212t3ENDKC4ulsVikdVqZdtBAH7JZnPsalJQ4Hrdt8kkdeok5eVJF/1yEgDwP0zdAgDcZjZLCxY4/nzpZk9VP8+fT/AGgNoQvgEAHklNlVatkiIjnds7dXK0p6Z6py4A8AcsOzEQy04ANCY2m2NXk8JCxxrvpCRmvAHgSnz2hksAgG8zm6VBg7xdBQD4F5adAAAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABgnwdgEAAP9ks0mZmVJhoRQRISUlSWazt6sCAN9G+AYAeCwjQ5o2TTp69Lu2Tp2kBQuk1FTv1QUAvo5lJwAAj2RkSKNHOwdvSSoocLRnZHinLgDwBya73W73dhFNRXFxsSwWi6xWq8LDw71dDgB4zGaToqNrBu8qJpNjBjwvjyUoAOAKy04AAC5VlFUo5z/7tWtzjnK27te54lKdOSNFHpUia/uQXVK+NPYGydLKrK69YpSQ3EO97uihVu0tBlYPAL6JmW8DMfMNwJfZ7XYdyv5an7+/W59vzFZO5heqKDtfb8fvkhit3nfEq8/3ExSfdIMCgwPr7dgA4C8I3wYifAPwNWXnyrVjY7Y+fWeHPn13p4qOnqh1bPOg5qqslCoqrnzcwECp8sIF1fZPTHBYkHoP6albhvXWzT/srXbXtqnrVwAAv0L4NhDhG4AvOFdSqk/WfabMjG3a/u5OlZe6TtPto9qqz5Ce6jWkpxIGxaltROvqNd8FBZKrfz0uXvNdUVamfVkHtPODPdqxMVu5O/JqrenGW7tpQGo/JY26RR2jr6mnbwoAvofwbSDCNwBvsl2wad1L7+vvs99UyamzNfqbBzVXYnKcbr6rt/p8v6c6dbtWJpOpxriq3U4k5wBeNXTVKtfbDVq/LdbOD/bo03/t1PYNO3W6qLjGGJPJpGEP3aFxz/yENeIAGiXCt4EI3wC85YttuZo/IV2Hsr92am/VPly3jrxJt47sq8TBPRQSFuzW8Vzt8x0VJc2f794+35WVlTqw/SttW/+5tq75VIf35jv1t2gVpof+eJ/ueugONWvGrrgAGg/Ct4EI3wCMVl5arldmrdTbCzc4rb++474k3fXQHeox4HqZ67gnYH0+4TL/QIG2vJGlVc+v07mS0ur2+Ntv0BOvTFFE5w51OzAA+BjCt4EI3wCMdCj7a/3hp3/RkS8Kqttie8XokcUP6cZ+3bxYWe1OHjull3/1mja++lF1W0iLYD2y6CENfWCgFysDgPpB+DYQ4RuAUd5fvkULJi2p3iqweVBz/fz3P1Hqoz+UOcD3n36z44M9+svDL+lY3vHqth+MH6xHFj+kwKDmXqwMAK4O4dtAhG8ADa2yslIvz3hNq55fV90W2ytGT776iK67McqLlXnuXEmpXnx0qf61dHN12423dtOct2eo9TXcjAnAPxG+DUT4BtCQLpy/oLk/X6xNK/5T3TZ8wlBN+ss4v36gzcbXPtL8CenVWyJ26hahP73/W13zvfZergwAPEf4NhDhG0BDsdls+uP9f9WWN7ZKkpqZm2nqXx/UiInf93Jl9SN3xyH9duQfdeK/pyRJEZ076PkPn1K7yLZergwAPMP+TQDg5+x2u16YtrQ6eDcPaq7Zqx9vNMFbkrr27qwFW/+gyK4RkqTCQ99o5l1/0Nnic16uDAA8Q/gGAD/37t8+0NoX3pMkmQPMmr1qum4beZOXq6p/Ha5rr3mbZldvO3g4J19/euCvtT7CHgB8EeEbAPzY1/vytXjaK9U//3LJBN3ywz5erKhhtYtsq2ff/bVatmkhSfp47Wda89d3vVwVALiP8A0Afsput+svE9KrtxMcOflO3Tku2ctVNbxOXSM087X/V/3z/818Xcfzv/ViRQDgPsI3APiprH9u196tByRJkV0j9PDc+71ckXFu+kEvjZh0pySpvLRCr855y8sVAYB7CN8A4IfsdrtWpmVU//zw3PsVFBLkxYqMN+7396pFqzBJ0r9f/YjZbwB+gfANAH7owPaDOrD9K0lSl8Ro3Tqir5crMl54m5b60dQfSJJsF2zasGSjlysCgCvz6fBdVlamp59+Wt27d1dwcLAiIiJ0//33Ky8vz+NjVVZWauHChUpISFBoaKjatWunH//4x9q1a9dlP/fJJ59o+PDhatu2rcLCwtSnTx8tWbKEu+sBeNUHr2VW//lHU34gk8nkxWq8Z/jE76uZ2fFP2Qevf8S1GYDP89nwXVpaqjvuuENz587VM888o5MnT+pf//qX9uzZo169emn37t1uH6uyslL33HOPpk+frokTJ6qoqEiffPKJSkpKdMstt+j99993+bmMjAwNGDBAAQEB2rFjh44dO6Z7771XEydO1Pjx4+vrqwKAR+x2uz5e95kkqXlggG4f3c8rddhs0pYt0sqVjnebzfga2l3bRomDe0iSjh0uUt6eI8YXAQAe8NnwPWfOHGVlZemPf/yj7r77boWGhiohIUGrV6/WmTNndN9998nm5pU+PT1dq1ev1qOPPqpJkyYpLCxMsbGx+sc//qGwsDDdf//9KikpcfrM8ePHNW7cOF177bVauXKlrrvuOrVs2VIzZszQAw88oGXLlunNN99siK8OAJd1/Mi3+ubrIklS3IDrFWYJM7yGjAwpOlpKTpbGjHG8R0c72o12y7De1X/O/nCf8QUAgAd8MnyfO3dOixcvVmBgoMaOHevU16VLFw0ZMkR79+7Vhg0b3DrevHnzJEkPP/ywU3vr1q11zz336Pjx41q2bJlT30svvaSSkhLdf//9CgkJceqbOHGiJOnPf/6zJ18LAOpF7o5D1X++sV83w8+fkSGNHi0dPercXlDgaDc6gN9463d/Bwd3er4sEQCM5JPhe9OmTTp79qwSEhLUokWLGv0DBgyQJK1du/aKx8rJydGhQ4fUvn17de3a1e1jVf3cv3//Gp/p27evgoODtWPHDhUUFFz5CwFAPSo8dLz6z9fd2MnQc9ts0rRpkqul1VVtjz5q7BKUi/8O/nvomHEnBoA68MnwXbWeOyYmxmV/586dncbV97FsNpv27t1b6+cCAgIUFRV1xRrKy8tVXFzs9AKAq1X87XfXktYdWxl67szMmjPeF7Pbpfx8xzijhLQIUXCYY5vF4m9LrjAaALzLJ8P3sWOOmYs2bdq47G/durXTuPo4VlFRkSorKyVJp0+fVllZ2VXXkJaWJovFUv2qCuwAcDXOV1yo/nNgcKCh5y4srN9x9aV5UHNJzn83AOCLfDJ8nzt3TpIUGOj6H5WgoCCncfVxrIvHXnzcq6lh5syZslqt1a/8/Pwr1gsAV1IVNCXH0x2NFBFRv+Pqy/my85KkwIv+bgDAFwV4uwBXQkNDJUkVFa7/USkvL3caVx/Hunjsxce9mhqCgoKcwj0A1AdLu5bVfz5ZeMrQcyclSZ06OW6udLXu22Ry9CclGVfT2eJzKjvnuCZb2ocbd2IAqAOfnPnu2LGjJOnkyZMu+0+dcvxjE+HG1Iq7x7rmmmvUrJnjr6NVq1YKDg6utxoAoD5dG9ux+s9f77vMAuwGYDZLCxY4/nzpc32qfp4/3zHOKBf/HUR07mDciQGgDnwyfCcmJkpSrU+yPHTIsc1Wz549G+RYZrNZPXr0qPVzFy5cqF5C4k4NAFCfYnt9dyP43qz9hp8/NVVatUqKjHRu79TJ0Z6aamw9e7ceqP5z196djT05AHjIJ8N3cnKywsLClJ2drbNnz9bo37p1qyRp5MiRVzxWXFycOnfurKKiIh08eNDtY40YMUKSlJWVVeMzn3/+ucrKytS7d29FXvqvDwA0sGui2unaLo4Z3n1ZX6rk1BnDa0hNlQ4fljZvllascLzn5RkfvCXp0w2fV/85YdCNxhcAAB7wyfAdGhqqqVOnqry8XMuXL3fqO3TokDZu3Ki4uDgNGzasun3Pnj3q37+/Zs2aVeN4TzzxhCRpyZIlTu2nT5/WW2+9pWuuuUbjxo1z6ps4caJatmypV199tXrnkyrp6emSpBkzZtT5OwLA1bh15E2SJNsFm7a8sdUrNZjN0qBB0k9/6ng3cqlJleNHirR7i+OpltfGdtT3bjB233MA8JRPhm9Jmj17tm699VY9+eSTWr16tUpLS5Wdna1Ro0YpNDRUr7/+uswXXenT09OVlZWltLQ0nThxwulYDz/8sFJTU/WXv/xF6enpOnv2rL766ivdc889Kikp0auvvqqWLVs6feaaa67R0qVLVVBQoDFjxujIkSMqKSnRvHnztGzZMo0dO1b33nuvIX8XAHCpIfffXv3nfy7+V/VWqU3NPxe/J/v/7vwcev9AmS5diA4APsZnw3dISIg2bdqk6dOna+bMmWrdurWGDh2quLg47dy5UwkJCU7jU1JS1Lp1a6WkpNTYm7tZs2Z66623NG/ePC1evFjt2rXTzTffrNDQUG3btk3f//73XdYwatQo/ec//1F5ebkSEhLUoUMHrVixQi+++KKWLl3aYN8dAK4kNjFGcf27S3LccJi56hMvV2S8U8etWvfie5Kk5oEBuuuhO7xcEQBcmclud7VZFBpCcXGxLBaLrFarwsPZDgvA1dn+r52aNexZSVKH69rr5T3PKaRFiJerMs688S/ovWWbJUkjJt2p/7f4IS9XBABX5rMz3wCAy+t7Z6J63REvSfrm6yIt+n+veLki42Su/qQ6eIe2DNHPfjvKyxUBgHsI3wDgp0wmkx596WEFhzke5vX+si1a89d3vVxVwzuU/bXmjX+h+ucJz41Vm46tvVgRALiP8A0AfuzaLh31y/QJ1T+/+Mul2vKmd3Y/McJ/vzqmWcP+oHMlpZKkQT/pr7seHOzlqgDAfYRvAPBzg8ck6SdP/liSVFlpV9rPFuqD1zO9XFX9yz9QoMcHz9GJ/zqeMNz9pi6a/rdJ7HACwK8QvgGgEfj5Mz/RXQ86dvuotFXqj/cv1Jt//qcayz31Of/5Qr9M+q2K8h1byV53Yyc9s36mgkODvFwZAHiG8A0AjUCzZs007aVfaPiEodVtf3vyNaX9bIFKz5R6sbKrY7fbtfaF9/TEHU/J+m2JJKlzwnWa+8FstWpv8XJ1AOA5tho0EFsNAmhodrtdrz+zWstnv1nd1qlbhJ589f+p+02xXqzMc9Zvi7Vg0hJlrt5W3dZ7aE/99s3H1KJVmBcrA4C6I3wbiPANwCiZGds07+eLq29MbGZupp/8KkX3/WaUAoMDvVzdlf3n7W3665S/6eSx09Vtox79oX7x5/tlDvDCc+wBoJ4Qvg1E+AZgpKO5hXp2zHzlfn6ouq1TtwhNWfig+n4/4TKf9J7CvG/00mPLlfXP7dVtLdu00GMvT9SAH9/ixcoAoH4Qvg1E+AZgtAvnL2jFHzK04tkM2S7Yqtv7De+jux68Q32+31NBId69adFut+vA9oPa8sZWrX3xfZ0vP1/dd8sPe+uXSyaqbQT7eANoHAjfBiJ8A/CWvJwjmj9xifZlHXBqDw4L0s3Deqvf8D66+a5esrQz5tpUUX5e2R/u07b1n2vrmk9VdPSEU3+bjq00ef7Pdfvdt7KVIIBGhfBtIMI3AG+qrKzUB69n6m9Pvq6Thadq9JtMJl1/S6xuGdZHvYf2VLe+nWU219/66sJD3+jzf2dr+792asfGbJWdLa8xJqC5WanTfqgxvxmlsPDQejs3APgKwreBCN8AfEFF+Xnt/GCP/rP6E2Wt/UzFJ0pcjmvRKkyJg3uo95CeSkyOU8eYazw6z5nT55T94T7t/GCPdmzMVuGhb1yOC2huVq8hPZWUeotuHdmXLQQBNGqEbwMRvgH4GtsFm/ZmHdC2d3bo0w07dHhvviHnbdU+XDcN66VbhvVR3+/3VJiFrQMBNA2EbwMRvgH4umOHj+vz93drxwd7tPODPSo5eaZejhvQ3Ky4/ter1x3x6jO0p7r17aJmzXjOG4Cmh/BtIMI3AH9is9l0cOdh7dqUo5z/fKGzxec8+rw5wKxuvTsrIbmH4m+/QSFhwQ1UKQD4D8K3gaxWq1q1aqX8/HzCNwAAgA9r2bJlg+y2FFDvR0StTpxwbKUVFRXl5UoAAABwOcePH1f79u3r/biEbwO1adNGknTkyBFZLNzNj4ZTXFysqKgofsuCBsd/azAK/63BKFX/rQUGBjbI8QnfBqq6uchisXDhgCHCw8P5bw2G4L81GIX/1mCUhnrAF7eaAwAAAAYhfAMAAAAGIXwbKCgoSLNnz1ZQUJC3S0Ejx39rMAr/rcEo/LcGozT0f2tsNQgAAAAYhJlvAAAAwCCEbwAAAMAghG8AAADAIIRvAAAAwCCEb6ARKS0t1e9+9zt169ZNwcHBuvbaazV+/HgVFBR4uzQAACDCtyE+++wzPfDAA4qNjZXJZNJvfvMbb5eERqisrEyDBw/W73//e505c0Y/+tGPFBUVpaVLl6pXr146dOiQt0uEH+G6BaCxWb58ufr27atWrVopLCxMvXv31htvvGF4HTxe3gBbt27VJ598ogEDBujbb7/1djlopJ555hl98sknuvXWW/X++++rRYsWkqTnn39e06dP1/jx47VlyxbvFgm/wXULQGNz6tQppaSkKDExUcHBwVqzZo1++tOfKjg4WCkpKYbVwT7fBqisrFSzZo5fMkRHR+tnP/uZnnnmGS9XhcakoqJC11xzjaxWq3bs2KFevXo59SckJCg7O1ufffaZ+vTp46Uq4U+4bsHbPvvsM910001ujR0yZIj+/e9/N3BFaIwGDBigiIgI/eMf/zDsnMx8G6DqHzCgoWzdulVWq1VdunSpEbwlafTo0crOzta6desI33AL1y14W2lpqX72s5/ptddeU4sWLTRq1Cin/jfffFNlZWV64IEHdNddd3mpSvi7tm3b6vz584aek/ANNAK7d++WJPXu3dtlf1V7dna2YTUBwNVISkpSu3bt9Nprr6lv375atmxZdd+ZM2f06quvqlOnTlq+fLn3ioRfunDhgs6dO6d3331X//73v7V69WpDz0/4BhqBI0eOSJI6derksr+q/euvvzasJgC4Wnv27JEkxcfHO7VnZ2ersrJSPXv29EZZ8GPHjh1TRESEJMlsNuuFF14w/DcnhG+gEThz5owkKTQ01GV/WFiYJKmkpMSwmgDgalX9tu7SkL1r1y6X7cCVtGvXTtu3b1dJSYn+9a9/aerUqWrbtm2NZU0NifANAAB8UlX4vnTmm/CNugoICFDfvn0lScnJyTp58qRmzpxpaPjmjhqgEajaVvDcuXMu+8+ePStJatmypWE1AcDVys7OlslkUo8ePZzaq8J3QkKCF6pCY5KYmGj4czAI3z7os88+k8lkcus1dOhQb5cLH/C9731PknT06FGX/VXt1113nWE1oengmoWGUFxcrK+//lqdO3euXjonSTabTTk5OQoKClL37t29WCEag6ysLEVHR9dob8jrGstODFBUVKQPP/xQkmNmcv/+/Vq1apXCwsJcLvJneyV4qmr2Z8eOHS77q9r5FS3c5cl1i2sWGkJt671zc3NVWlqqXr16yWw2e6M0+Knk5GSNGjVK119/vcrKyvTPf/5TK1as0JIlS2qMbdDrmh0NbvPmzXZJNV7XXXddrZ/Zt2+fXZJ90KBBTu0lJSX2Zs2a2Tt16tTAVcOflJeX2y0Wi12SfefOnTX6e/bsaZdk/+yzz4wvDn7J0+sW1yzUt8WLF9sl2X/3u985tb/99tt2Sfaf/OQnXqoM/mratGn2bt262UNCQuzt2rWz33777fZ169bVOr6hrmssOzHAoEGDZLfba7wOHz5c62fYXgmeCAwM1NSpUyVJU6ZMqV7jLTkeL5+dna2BAwfygB24zdPrFtcs1LfaZr6rdnfKzs7WO++8o9LSUsNrg3+aP3++Dhw4oHPnzlX/dm/48OG1jm+o6xrLTtzw4x//WF988YVHn/n73/+um2++uc7nZHsleOo3v/mNNm7cqKysLHXt2lVJSUn6+uuvtW3bNrVv316vvPKKt0tEI8Y1C/WttuAzePBgJSYmau/evXrssce0b98+b5SHJqChrmuEbzfk5eXpwIEDHn2mtl0n3MX2SvBUcHCwNm/erLS0NK1YsUJr1qxRmzZtNG7cOP3+97+v9QE8QH3gmoX6tnXrVpft1157rXbu3GlwNWiKGuq6ZrLb7farqgwNIjo6WkeOHFFJSYnTXd4333yztm/frr179+rGG2/0YoUA8B2uWQAam4a6rhG+G5jJZHJ7bNX/FMXFxbJYLOrSpYsOHjxY3W+z2dSyZUtVVlbq7Nmz3OUNoEF4et3imgXA1/nSdY0bLhuY3W7Xm2++qWHDhqljx44KDw9XUlKSMjMza9zIVOVK2yvdeOON/CMGoMF4et3imgXA1/nSdY3wbYD58+erXbt2Wrx4sf7xj38oMjJSd9xxh3bv3u1yfG1rjPbv3y9JPFQAQIPz5LrFNQuAP/CV6xo3XBpg3bp1atu2bfXPQ4YMUXx8vBYvXuxyY3d3t1caPHiwQkJCGrByAE2VJ9ctrlkA/IGvXNeY+TbAxf9DS1KzZs3Uo0cP5eXluRx/pe2VcnNz9dhjjykwMLBhCgbQ5Hly3eKaBcAf+Mp1jRsuvcBms6lr166688479eKLL3q7HAC4Iq5bABobb13XmPn2gkWLFunIkSOaPHmyt0sBALdw3QLQ2Hjrukb4Nti2bdv05JNP6je/+U2NX2UAgC/iugWgsfHmdY1lJwY6fPiw+vXrp9tvv11vvvmmR3tOAoA3cN0C0Nh4+7pG+DbI6dOn1b9/f7Vs2VKbN2/mjn8APo/rFoDGxheua2w1aICKigqlpqbq3Llz2rRpE/+AAfB5XLcANDa+cl0jfBtg8uTJ+vDDD/Xyyy8rLy+vekuboKAg9erVy8vVAUBNXLcANDa+cl1j2YkBoqOj9fXXX9dov+6663T48GHjCwKAK+C6BaCx8ZXrGuEbAAAAMAhbDQIAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwAAAAYhfAMAAAAGIXwDAAAABiF8AwA8VlhYqClTpmjChAmqqKjQtGnT1KFDB7Vr106//e1vVVlZ6e0SAcAnEb4BAB7597//rbi4OL3wwgu68cYb9dRTT2nhwoXq2rWrysvL9cwzz+iZZ57xdpkA4JNMdrvd7u0iAAD+4eDBg0pMTFRZWZlef/11jR49Wh06dFBQUJCOHDmiNWvWaPTo0bJYLDp58qSaNWOOBwAuxlURAOC2X/3qVzp79qwef/xx3Xvvvdq1a5dOnDih22+/XWazWampqWrVqpWsVqvy8/O9XS4A+BzCNwDALd98843++c9/qkWLFpo9e7YkaceOHZKkxMRESZLJZFL79u0lSTabzSt1AoAvI3wDANyyefNm2Ww23XnnnQoJCZEk7dq1S5KUkJBQPe7EiROSpIiICMNrBABfR/gGALhl9+7dkqRbbrmluq0qfFfNfOfn5+vkyZO6/vrrqwM6AOA7hG8AgFu+/fZbSVKrVq0kSXa7XdnZ2erQoYM6duwoSXrnnXckScOGDfNKjQDg6wjfAAC3tGzZUpJ09OhRSVJubq7OnDlTveSkrKxMzz//vJo3b66JEyd6rU4A8GWEbwCAWwYMGCBJ+vvf/y6r1eq05MRqteruu+9Wbm6uZs2apa5du3qxUgDwXezzDQBwi81m06233qrt27crKipKbdq00e7du9WrVy/l5+fr22+/1aRJk7R48WKZTCZvlwsAPonwDQBw24kTJ/TEE0/on//8p06ePClJCgsLU58+ffTLX/5SKSkp3i0QAHwc4RsAUCft2rXT+fPndfr0aWa6AcBNrPkGAHjs66+/1okTJ9S7d2+CNwB4gPANAPBY1ZMt+/Tp4+VKAMC/EL4BAB4jfANA3bDmGwAAADAIM98AAACAQQjfAAAAgEEI3wAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBCN8AAACAQQjfAAAAgEEI3wAAAIBBCN8AAACAQf4/CZIZkLbR+GkAAAAASUVORK5CYII=",
      "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(r'$\\varphi$')\n",
    "plt.ylabel(r'$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = [r'$-\\frac{1}{2}\\pi$', r'$0$', r'$\\frac{1}{2}\\pi$', r'$\\pi$', r'$\\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": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAInCAYAAABAyzWwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJeUlEQVR4nO3de3RTVd7/8U8ItrRAw1WpbbWFclFoAa8oVKjicmRAO7XqM6jAknGQyywcUEZwRtBx7JoBfYoPXorziDoKowMRuc2oCNU81gsqUIqiVYqUTh2QS4rQFkzz++P8UghJISntSdK8X2tlxey9c/KtZh0/3d1nH4vb7XYLAAAAQItrE+oCAAAAgGhB+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+DaR2+1WdXW12FodAAAgOhG+TXT48GHZbDYdPnw41KUAAAAgBAjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASdqGugAAQGRyuSSHQ6qqkhITpawsyWoNdVUAEN6Y+QYABM1ul9LTpexsaexY4zk93WgHADSO8A0ACIrdLuXlSRkZ0ocfSocPG88ZGUY7ARwAGmdxu93uUBcRLaqrq2Wz2eR0OpWQkBDqcgAgaC6XMcOdkSGtXCm1OWkKp75eysmRSkulsjKWoACAP8x8AwAC5nBIu3ZJc+Z4B2/JeD17tlRebowDIoHLJRUVScuWGc8uV6grQmvHBZcAgIBVVRnPAwb47/e0e8YB4cxul6ZPl/bsOdGWnCwtXCjl5oauLrRuzHwDAAKWmGg8l5b67/e0e8YB4cpz7cLJwVuSKiu5dgEtizXfJmLNN4BIx5pvtAYul5Sa6hu8PSwWYwa8vJzvMZofM98AgIBZrdITT0hr1hhB++TdTnJyjPYFCwgsCG8OR+PBW5LcbqmigmsX0DJY8w0ACEpurrR8uTRzpnT11Sfa09KMdtbKItwFek0C1y6gJRC+AQBBy82Vbr6ZO1wiMgV6TQLXLqAlsObbRKz5BgAg9DxrvisrjSUmp2LNN1oSa74BAEBUsVqN7QQlI2ifzPO6oIDgjZZB+AYAAFHHc+1CUpJ3e3Iy1y6gZbHsxEQsOwEAILy4XFy7AHNxwSUAAIhaVqs0YkSoq0A0YdkJAAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBJuLw8AaBKXS3I4pKoqKTFRysoybtUNAGhc2M9819bW6qGHHlJMTIwmTJgQ9Pt37doli8VyxkdqaqrX+4qKioIaDwDRxG6X0tOl7Gxp7FjjOT3daAcANC6sZ743btyoSZMm6YcfftDx48fP6lg9e/bUOeec47evrKxMF110kU9727Zt1atXL7/vSUpKOqt6ACBS2e1SXp40erS0bJk0YIBUWio9/rjRvny5lJsb6ioBIDyFbfhetmyZJk+erMcff1xxcXG6++67z+p47777rt/Z6i1btmjw4MGaNm2aT19SUpJ27NhxVp8LAK2JyyXNnGkE75UrpTb//++nQ4YYr3NypPvvl26+mSUoAOBP2C47SUtL0xdffKEpU6bIYrE0+Tjt27fX7bffrvbt2/vtf+qpp9SrVy/deOONTf4MAIgWDoe0a5c0Z86J4O3Rpo00e7ZUXm6MA8KByyUVFRl/pSkqMl4DoRS2M99DhgxpluN0795df//73/327d+/X8uWLdPjjz+uNqf+XwQA4KOqyngeMMB/v6fdMw4IJbtdmj5d2rPnRFtysrRwIUujEDpRnTiff/55Wa3WRpe01NTUaNasWerTp4/i4uLUqVMnDRs2TIsXL1Z9fb3J1QJA6CUmGs+lpf77Pe2ecUCoeK5NODl4S1JlpdHOxcEIlagN3y6XS88884zGjRsnm83md8zevXvldDq1evVqOZ1Obd68Wf3799ekSZM0ZsyYM14EWldXp+rqaq8HAESyrCwpNdW4uPLUOYj6eik/X0pLM8YBoeJyGTPebrdvn6ftvvtYgoLQiNrwvXLlSlVUVPi90FKSMjIyVFRUpMLCQvXt21cxMTFKS0tTYWGhsrOztW7dOs2fP/+0n5Gfny+bzdbwSElJaYkfBQBMY7VKTzwhrVljXFz54YfS4cPGc06O0b5gARdbIrQcDt8Z75O53VJFBdcmIDSiNnw/9dRTuu6663TxxRf77e/atauGDx/ut2/q1KmSpBdeeOG0nzF79mw5nc6GR0VFxdkVDQBhIDfX2E5w2zbp6qulhATjubSUbQYRHgK95oBrExAKYXvBZUsqKSnR+++/r5UrVzbp/enp6ZKkb7/9VseOHVNMTIzfcbGxsYqNjW1qmQAQtnJzje0EucMlwlGg1xxwbQJCISrD91NPPaXU1FSNGTOmSe8/+WLLs9kGEQAimdUqjRgR6ioAX1lZxq4mlZX+131bLEY/1yYgFKJu2cmBAwe0dOlSTZky5bTbC15++eV67rnn/PaVlZVJknr16tXoXTMBAEBoWK3GdoKSEbRP5nldUMBfahAarSZ8b9u2TUOHDtWcOXNOO+7555+XJE2cOPG04/bt29fospRFixZJkiZMmBB0nQAAoOV5rk1ISvJuT07m2gSEVqsJ34WFhSouLlZ+fr7279/vd4xne8E77rhDXbp0Oe3xLBaL3nrrLU2aNElfffWV6urqtHPnTt1zzz1yOBy64YYbNGvWrJb4UQAAQDPIzTXuyLpxo7R0qfFcXk7wRmhZ3G5/q6HCw+nWUy9ZssRr5nn9+vW67bbbNHz4cNntdr/vfeONN5Sbm6utW7cqMzPztJ9dUVGh119/XatWrdL27dt16NAhdezYUZmZmbrzzjs1ceLEoO+KWV1dLZvNJqfTqYSEhKDeCwAAgMgX1uG7tSF8AwAARLdWs+wEAAAACHeEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCRtQ10AACAyuVySwyFVVUmJiVJWlmS1hroqAAhvzHwDAIJmt0vp6VJ2tjR2rPGcnm60AwAaR/gGAATFbpfy8qSMDOnDD6XDh43njAyjnQAOM7lcUlGRtGyZ8exyhboi4PQsbrfbHeoiokV1dbVsNpucTqcSEhJCXQ4ABM3lMma4MzKklSulNidN4dTXSzk5UmmpVFbGEhS0PLtdmj5d2rPnRFtysrRwoZSbG7q6gNNh5hsAEDCHQ9q1S5ozxzt4S8br2bOl8nJjHNCSPH+BOTl4S1JlJX+BQXgjfAMAAlZVZTwPGOC/39PuGQe0BJfLmPH297d7T9t997EEBeGJ8A0ACFhiovFcWuq/39PuGQe0BIfDd8b7ZG63VFHBX2AQngjfAICAZWVJqanS448ba7xPVl8v5edLaWnGOKClBPqXFf4Cg3BE+AYABMxqlZ54Qlqzxri48uTdTnJyjPYFC7jYEi0r0L+s8BcYhCN2OzERu50AaC3sdmnmTOPiS4+0NCN4s8sEWprLZfwFprLS/7pvi8XY9aS8nF8EEX4I3yYifANoTbjDJULJs9uJ5B3ALRbjeflyfhFEeCJ8m4jwDQBA8/G3z3dKilRQQPBG+CJ8m4jwDQBA8+IvMIg0bUNdAAAAQFNZrdKIEaGuAggcu50AAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmaRvqAgAAkcnlkhwOqapKSkyUsrIkqzXUVQFAeCN8AwCCZrdLM2dKu3adaEtNlZ54QsrNDVVVaE345Q6tFctOAABBsdulvDwpI0P68EPp8GHjOSPDaLfbQ10hIp3dbvwyl50tjR1rPKem8t1C62Bxu93uUBcRLaqrq2Wz2eR0OpWQkBDqcgAgaC6XlJ5uBO2VK6U2J03h1NdLOTlSaalUVsYsJZrG88vdqenEYjGely/nryuIbMx8AwAC5nAYS03mzPEO3pLxevZsqbzcGAcEy+WSpk/3Dd7Sibb77jPGAZGK8A0ACFhVlfE8YID/fk+7ZxwQDIdD2rOn8X63W6qo4Jc7RDbCNwAgYImJxnNpqf9+T7tnHBCMQH9p45c7RDLCNwAgYFlZxoVvjz9urPE+WX29lJ8vpaUZ44BgBfpLG7/cIZKFffiura3VQw89pJiYGE2YMKFJx5g3b54sFkujj9Md95///KeuvfZa2Ww2JSQkaNiwYVq+fHnTfhgAiHBWq7Gd4Jo1xsWVJ+92kpNjtC9YwMWWaJqsLCk5+cTFlaeyWKSUFH65Q2QL6/C9ceNGZWZm6tlnn9Xx48fP6lgJCQnq27ev30diI79CFxQUaNSoUerdu7e+/vprlZeX64orrtCtt96qefPmnVU9ABCpcnONHSe2bZOuvlpKSDCeS0vZiQJnx2qVFi40/vnUAO55XVDAL3eIbGG71eCyZcs0efJkPf7444qLi9Pdd9+t8ePH68UXXwz6WPPmzdOuXbuCem9paakGDRqkSy65RB9//LEsJ50FsrOzVVRUpA8//FBDhgwJ+JhsNQigNeEmKGgpdrux68nJF1+mpBjBm1/uEOnC9g6XaWlp+uKLL3T++ec3KXCfrSeffFIul0v33HOPV/CWpEmTJqmoqEjz58/XihUrTK8NAMKB1SqNGBHqKtAa5eZKN9/ML3doncI2fAczo9wSVq9eLUkaOnSoT9+wYcMkSevWrdNPP/2ktm3D9l8jAAARiV/u0FqF9Zrv5lReXq6xY8cqJSVFMTExOu+883TTTTfp3Xff9RlbWVmpH374QZIxA3+qpKQkxcTEqLa2Vl9//XWL1w4AAIDWIWrC9yeffKLLLrtMn3zyiX788UetXbtWBw8e1MiRI/XYY495jf3+++8lSe3atVNcXJzPsSwWizp16uQ11p+6ujpVV1d7PQAAABC9oiJ8jx07Vps2bdKMGTOUmJiomJgYXXbZZVqzZo06d+6sP/zhD/rggw8axh89elSSFBMT0+gxY2Njvcb6k5+fL5vN1vBISUlppp8IAAAAkSgqwnefPn00wM+9kG02m+644w5J0pIlSxra4+PjJUnHjh1r9Jh1dXVeY/2ZPXu2nE5nw6OioqJJ9QMAAKB1iPorBdPT0yVJX375ZUNbjx49JBk3+KmpqfFZeuJ2u3Xo0CFJanSPcMmYHffMkAMAAABRMfN9OvX///7IJ28nmJSUpG7dukkyLtQ8VWVlpY4dO6a4uDj17t3bnEIBAAAQ8Vp9+K6oqFCPHj1UUlLit7+srEyS1K9fP6/2MWPGSJKKi4t93uNZH37jjTeyzSAAAAAC1mrC97Zt2zR06FDNmTPHq93lcuk///mP1q1b5/OegwcPaunSpZKkCRMmePXNnDlTVqtVf/3rX3XqTUAXL14sSXrggQea8ScAAABAa9dqwndhYaGKi4uVn5+v/fv3N7R7lpM88sgjKigo0J49e3Ts2DFt2rRJo0ePltPp1Lx58xpunOPRv39/LViwQB9//LGmTJmivXv36sCBA7r//vu1YcMGzZ07N+Q3AgIAAEBksbhPndYNI6fe1v1kS5Ys8ZqtXr9+vW677TYNHz5cdrvd672fffaZXn/9db399tv65ptvVFtbq27dumnIkCGaOnWqRo4c2ejnrFu3TvPnz9dnn30mt9utzMxM/fa3v1VeXl7QP091dbVsNpucTqcSEhKCfj8AAAAiW1iH79aG8A0AABDdWs2yEwAAACDcEb4BAAAAk7BPHgCgSVwuyeGQqqqkxEQpK0uyWkNdFcIV3xfAQPgGAATNbpdmzpR27TrRlpoqPfGElJsbqqoQrux2afp0ac+eE23JydLChXxfEH1YdgIACIrdLuXlSRkZ0ocfSocPG88ZGUa73R7qChFOPN+Xk4O3JFVW8n1BdGK3ExOx2wmASOdySenpRtBeuVJqc9IUTn29lJMjlZZKZWUsKYDxfUlN9Q3eHhaLMQNeXs73BdGDmW8AQMAcDmOpyZw53sFbMl7Pnm0EKYcjJOUhzDgcjQdvSXK7pYoKvi+ILoRvAEDAqqqM5wED/Pd72j3jEN0C/R7wfUE0IXwDAAKWmGg8l5b67/e0e8YhugX6PeD7gmjCmm8TseYbQKRjzTeC4VnzXVlpLDE5FWu+EY2Y+QYABMxqNbYTXLPGCNon73aSk2O0L1hAkILBajW2E5SMoH0yz+uCAr4viC6EbwBAUHJzpeXLpW3bpKuvlhISjOfSUqOdfZtxMs/3JSnJuz05me8LohPLTkzEshMArQl3LEQw+L4ABsK3iQjfAAAA0Y1lJwAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEnahroAAAAQOVwuyeGQqqqkxEQpK0uyWkNdFRA5CN8AACAgdrs0fbq0Z8+JtuRkaeFCKTc3dHUBkcTidrvdoS4iWlRXV8tms8npdCohISHU5QDAWWEGNLrY7VJennRqarBYjOflywngQCAI3yYifANoLex2aeZMadeuE22pqdITTxDAWiOXy/jve/KM98ksFmMGvLycX8CAM+GCSwBAUDwzoBkZ0ocfSocPG88ZGUa73R7qCtHcHI7Gg7dkzIZXVBjjAJwe4RsAEDCXy5jxHj1aWrlSGjJE6tDBeF650mi//35jHFqPqqrmHQdEM8I3ACBgDoex1GTOHKnNKf8HadNGmj3bWHrADGjrkpjYvOOAaEb4BgAEzDOzOWCA/35POzOgrUtWlrGm23Nx5aksFiklxRgH4PQI3wCAgHlmNktL/fd72pkBbV2sVmM7Qck3gHteFxRwsSUQCMI3ACBgWVnGrhePPy7V13v31ddL+flSWhozoK1Rbq6xnWBSknd7cjLbDALBYKtBE7HVIIDWwLPbyejRxhrvAQOMGe/8fGnNGoJYa8f+7sDZIXybiPANoLXwt893Wpq0YAHBGwBOh/BtIsI3gNaEGVAACF7bUBcAAIhMVqs0YkSoqwCAyMIFlwAAAIBJCN8AAACASQjfAAAAgEnCPnzX1tbqoYceUkxMjCZMmNCkY2zatEm/+c1vNGDAAMXHx6tdu3ZKT0/X5MmTVV5e7vc9RUVFslgsjT5SU1Ob/kMBAAAgKoX1BZcbN27UpEmT9MMPP+j48eNNOsbatWs1evRoXXjhhVq4cKGuueYaHT9+XKtWrdKMGTP0yiuv6K233tLVV1/t8962bduqV69efo+bdOpdBgAAAIAzCNvwvWzZMk2ePFmPP/644uLidPfddzfpODU1NZKk1157TVdeeWVD+69+9Sv99NNPmjx5siZOnKgvv/zS571JSUnasWNH034AAAAA4BRhu+wkLS1NX3zxhaZMmSKLxdLk43Tp0kXXX3+9V/D2uOmmmyRJO3bs0O7du5v8GQAAAEAgwnbme8iQIc1ynGuvvVbXXnut3z6bzdYsnwEAAAAEImxnvs3gWVLSs2dPXXDBBT79NTU1mjVrlvr06aO4uDh16tRJw4YN0+LFi1VfX292uQAAnBWXSyoqkpYtM55drlBXBESfqA7fL7/8siTp4Ycf9tu/d+9eOZ1OrV69Wk6nU5s3b1b//v01adIkjRkz5owXgdbV1am6utrrAQBAKNjtUmqqlJ0tjR1rPKemGu0AzBO14bukpETPPvusbrnlFo0fP96nPyMjQ0VFRSosLFTfvn0VExOjtLQ0FRYWKjs7W+vWrdP8+fNP+xn5+fmy2WwNj5SUlJb6cQAAaJTdLuXlSXv2eLdXVhrtBHDAPFEZvvft26dbb71VV1xxRcPs96m6du2q4cOH++2bOnWqJOmFF1447efMnj1bTqez4VFRUXF2hQMAECSXS5o+XXK7ffs8bffdxxIUwCxhe8FlSzl48KB+9rOfqXv37lq7dq3i4+ODPkZ6erok6dtvv9WxY8cUExPjd1xsbKxiY2PPql4AAM6Gw+E7430yt1uqqDDGjRhhWllA1Iqqme99+/YpOztbnTt31ttvv93k3U5OvtjybLZBBACgpVVVNe84AGcnasL3v//9bw0fPlxJSUlas2aN14z3888/r+3bt3uNv/zyy/Xcc8/5PVZZWZkkqVevXjrnnHNarmgAAM5SYmLzjgNwdlpN+N62bZuGDh2qOXPm+PTt3r1b11xzjS666CK98cYbateunVf/n/70J23atMmrbd++fVq5cqXfz1q0aJEkacKECc1SOwBEIratiwxZWVJystTYH2otFiklxRgHoOW1mjXfhYWFKi4uVnFxsWbOnKmuXbtKknbu3Kns7GxVVlZq0KBBGjdunM979+3b59NmsVj01ltvadKkSZoxY4ZSU1NVWVmp/Px8ORwO3XDDDZo1a1aL/1wAEI7sdmnmTGnXrhNtqanSE09Iubmhqgr+WK3SwoXGriYWi/eFl55AXlBgjAPQ8sI6fJ+6nvqll17SSy+9JElasmSJ18xzTk6Oli5dquHDh6tLly4N7atWrWq4dfyKFSsC/uz3339fr7/+ulatWqWhQ4fq0KFD6tixozIzM7V48WJNnDhRbdq0mj8cAEDAPNvWjR5tzHoPGCCVlkqPP260L19OAA83ubnGf5fp070vvkxONoI3/70A81jcbn+bD6ElVFdXy2azyel0KiEhIdTlAEDQXC4pPV3KyJBWrpROnoOor5dycowgXlbGTGo4crmMXU2qqow13llZ/HcCzBbWM98AgPDicBhLTZYt8w7ekvF69mzp6qvZti5cWa38dwFCjXUTAICAebajGzDAf7+nnW3rAMA/wjcAIGCe7ehKS/33e9rZtg4A/GPNt4lY8w0g0rHmGwDODjPfAICAWa3GdoJr1hhB+8MPpcOHjeecHKN9wQKCNwA0hgsuAQBB8WxbN3OmcXGlR1oa2wwCwJmw7MRELDsB0JqwbR0ABI+ZbwBAk7BtHQAEjzXfAAAAgEkI3wAAAIBJCN8AAACASVjzDQBAmOKiVqD1IXwDABCG7HZp+nRpz54TbcnJ0sKFbOcIRDKWnQAAEGbsdikvzzt4S1JlpdFut4emLgBnj32+TcQ+3wCAM3G5pNRU3+DtYbEYM+Dl5SxBASIRM98AAIQRh6Px4C1JbrdUUWGMAxB5CN8AAISRqqrmHQcgvBC+AQAII4mJzTsOQHghfAMAEEaysow13RaL/36LRUpJMcYBiDyEbwAAwojVamwnKPkGcM/rggIutgQiFeEbAIAwk5srLV8uJSV5tycnG+3s8w1ELrYaNBFbDQIAgsEdLoHWhztcAgCahGDY8qxWacSIUFcBoDmx7AQAEDS7XUpPl7KzpbFjjef0dO68CABnQvgGAATFc+vzjAzpww+lw4eN54wMbn0OAGfCmm8TseYbQKRzuYwZ7owMaeVKqc1JUzj19VJOjlRaKpWVsQQFAPxh5hsAEDCHQ9q1S5ozxzt4S8br2bOl8nJufQ4AjSF8AwAC5rml+YAB/vs97dz6HAD8I3wDAALmuaV5aan/fk87tz4HAP9Y820i1nwDiHSs+T6BrRYBNAUz3wCAgFmt0hNPSGvWGEH75N1OcnKM9gULWn8Itdul1FTvrRZTU9npBcCZMfNtIma+AbQWdrs0c6Zx8aVHWpoRvFv7rc89Wy2e+n9Pi8V45vbvAE6H8G0iwjeA1iQal124XMYM9549/vstFik52djxpbX/uwDQNNxeHgDQJNF463OHo/HgLRmz4RUVxrho+3cDIDCs+QYAIECBbqHIVosAGkP4BgAgQIFuochWiwAaQ/gGACBAWVnGmm7PxZWnsliklBRjHAD4Q/gGACBAVqu0cKHxz6cGcM/rggIutgTQOMI3AABByM01thNMSvJuT05mm0EAZ8ZWgyZiq0EAaD2icatFAGePrQYBAGiCaNxqEcDZI3wDACIaM9AAIgnhGwAQsex2afp07xvfJCcbF0Wy9hpAOOKCSwBARLLbpbw83ztOVlYa7XZ7aOoCgNMJ+/BdW1urhx56SDExMZowYUKTj1NfX6+nnnpKAwcOVHx8vLp166Zf/OIX2rJly2nf99FHH2n06NHq2rWr2rdvr0svvVSLFy8W16kCQOi4XMaMt79TsaftvvuMcQAQTsI6fG/cuFGZmZl69tlndfz48SYfp76+Xrfddptmzpype++9V/v27dNHH32kw4cP68orr9Tbb7/t9312u13Dhg1T27Zt9fnnn+v777/X7bffrnvvvVd33313k+sBAJwdh8N3xvtkbrdUUWGMA4BwErbhe9myZfrFL36h++67T0888cRZHauwsFArVqzQfffdp8mTJ6t9+/ZKT0/XP/7xD7Vv31533XWXDh8+7PWevXv3asKECTr//PO1bNkyXXjhherYsaNmzZqlcePG6cUXX9Rrr712VnUBAJqmqqp5xwGAWcI2fKelpemLL77QlClTZGnsPr4BWrBggSTp17/+tVd7586dddttt2nv3r168cUXvfqee+45HT58WHfddZfi4uK8+u69915J0l/+8pezqgsA0DSJic07DgDMElT43rhxow4ePNhStXgZMmSIzj///LM+TmlpqXbu3Knu3burd+/ePv3Dhg2TJK1atcqr3fN66NChPu+57LLL1K5dO33++eeqrKw86xoBIBK5XFJRkbRsmfFs5vrqrCxjV5PG5mYsFiklxRgHAOEkqPB93XXXqVu3brrwwgt18803a+7cuXrjjTdUXl7uNW7ixIl65plnmrXQptq6daskYybdn549e3qNkySXy6Xt27c3+r62bdsqJSXF530AEC3sdik1VcrOlsaONZ5TU83bYcRqNbYTlHwDuOd1QQH7fQMIP0Ht83333Xdry5YtKi0t1erVq7V69eqGJSEJCQnKzMxUnz599Oabb2rt2rWaMmVKixQdjO+//16S1KVLF7/9nTt3liTt27dP9fX1atOmjQ4dOqTa2tqA3uc5vj91dXWqq6treF1dXR38DwAAYcazxd+pO414tvhbvtycPbZzc43P8rfPd0EB+3wDCE9Bhe+//vWvkqSffvpJX3zxhbZs2aLNmzdr8+bN+uijj+RwOPR///d/crvduuCCC1qk4GAdPXpUkhQTE+O3PzY21mtshw4dGt4TyPtOHnuq/Px8PfLII0HXDADh6kxb/FksxhZ/N99szqxzbq7xWdzhEkCkaNIdLtu2bavMzExlZmZq3LhxkqQff/xRL774oh588EENGDBAr7zySrMW2lTx8fGSpGPHjvntP3lm2jPW8xzI+04ee6rZs2drxowZDa+rq6sblqsAQCQKZou/ESPMqclqNe+zAOBsNdvt5Tt06KBp06apZ8+eGjNmjD7//HOlp6c31+GbrEePHpKkAwcO+O33XEB67rnnqk0bYwl8p06d1K5dO9XW1urAgQM677zzGn1f4mkupY+NjfWaWQeASMcWfwBwdpp9q8FRo0apX79+ys/Pb+5DN8mgQYMkyeeiUI+dO3dKkjIzMxvarFarBgwY0Oj7fvrpJ1VUVPi8DwBaO7b4A4Cz0yL7fKelpemrr75qiUMHrX///urZs6f27dunb775xqf/gw8+kCTddNNNXu1jxoyRJBUXF/u857PPPlNtba0uueQSJSUltUDVABCe2OIPAM5OUOF7+vTpWrJkibZs2XLa271/8803Dcs9zLJt2zYNHTpUc+bM8el74IEHJEmLFy/2aj906JBef/11nXvuuZowYYJX37333quOHTvqb3/7W8POJx6FhYWSpFmzZjXjTwAA4Y8t/gDg7AS15vt//ud/GrYWPOecc9SvXz8NHjxYgwYNUkZGhuLj4/XKK6+orKxMzz33XIsU3JjCwkIVFxeruLhYM2fOVNeuXRv6fv3rX+udd97Rf//3f6tXr16688479f3332vy5Mk6fPiw1qxZo44dO3od79xzz9WSJUt0++23a+zYsSooKFDnzp1VWFioF198UePHj9ftt99u6s8IAOGALf4AoOksbre/DaP8W7NmTcPWgps3b9Z3331nHOSU6Y+ePXtq/PjxuuSSS3TppZf6vWAxoOJOc1v5JUuWeM1Wr1+/XrfddpuGDx8uu93u816Xy6VFixbpf//3f1VWVqb4+HhlZWVp7ty5Gjx4cKOf89FHH+mPf/yjiouLVVdXp379+mnSpEn69a9/HfRt76urq2Wz2eR0OpWQkBDUewEg3LhcbPEHAMEKKnyf6uDBg15h/PPPP9fXX3+t+vp64+D/P5z26NFDl156qc8t3KMN4RsAACC6nVX49ufo0aMqKSlpCOObN29WaWmpjh8/LpfL1ZwfFXEI3wAAANGt2cO3P547Ykb7tnyEbwAAgOjWIlsNnspzR0wAAAAgmpkSvgEAAAAQvgEAAADTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTtA11AQCAyORySQ6HVFUlJSZKWVmS1RrqqgAgvBG+AQBBs9ul6dOlPXtOtCUnSwsXSrm5oasLAMIdy04AAEGx26W8PO/gLUmVlUa73R6augAgEljcbrc71EVEi+rqatlsNjmdTiUkJIS6HAAImsslpab6Bm8Pi8WYAS8vZwkKAPjDzDcAIGAOR+PBW5LcbqmiwhgHAPBF+AYABKyqqnnHAUC0IXwDAAKWmNi84wAg2hC+AQABy8oy1nRbLP77LRYpJcUYBwDwRfgGAATMajW2E5R8A7jndUEBF1sCQGMI3wCAoOTmSsuXS0lJ3u3JyUY7+3wDQOPYatBEbDUIoDXhDpcAEDzucAkAaBKrVRoxItRVAEBkYdkJAAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgkrahLgAAEJlcLsnhkKqqpMREKStLslpDXRUAhDfCNwAgaHa7NH26tGfPibbkZGnhQik3N3R1AUC4Y9kJACAodruUl+cdvCWpstJot9tDUxcARIKwDt+1tbV69NFH1bdvX7Vr106JiYm66667VF5eHvAxdu3aJYvFcsZHamqq1/uKioqCGg8A0cDlMma83W7fPk/bffcZ4wAAvsJ22UlNTY1GjhypkpISvfDCC/r5z3+usrIyjR8/XoMHD9Z7772ngQMHBny8nj176pxzzvHbV1ZWposuusinvW3bturVq5ff9yQlJQX82QDQWjgcvjPeJ3O7pYoKY9yIEaaVBQARI2zD97x581RcXKxFixbp1ltvlSQNHDhQK1asUN++fXXHHXdo69atsgZ4dc+7777rd7Z6y5YtGjx4sKZNm+bTl5SUpB07dpzVzwEArUlVVfOOA4BoE5bLTo4ePaqnn35aMTExGj9+vFdfr169NHLkSG3fvl3r1q0747Hat2+v22+/Xe3bt/fb/9RTT6lXr1668cYbm6V2AGjNEhObdxwARJuwDN8bNmzQkSNHNHDgQHXo0MGnf9iwYZKkVatWnfFY3bt319///nd1797dp2///v1atmyZpk6dqjZtwvJfBQCElawsY1cTi8V/v8UipaQY4wAAvsIycW7dulWSlJaW5re/Z8+eXuOa6vnnn5fVatXdd9/tt7+mpkazZs1Snz59FBcXp06dOmnYsGFavHix6uvrz3j8uro6VVdXez0AIJJZrcZ2gpJvAPe8Lihgv28AaExYhu/vv/9ektSlSxe//Z07d/Ya1xQul0vPPPOMxo0bJ5vN5nfM3r175XQ6tXr1ajmdTm3evFn9+/fXpEmTNGbMGB0/fvy0n5Gfny+bzdbwSElJaXK9ABAucnOl5culU687T0422tnnGwAaF5bh++jRo5KkmJgYv/2xsbFe45pi5cqVqqio8HuhpSRlZGSoqKhIhYWF6tu3r2JiYpSWlqbCwkJlZ2dr3bp1mj9//mk/Y/bs2XI6nQ2PioqKJtcLAOEkN1fatUvauFFautR4Li8neAPAmYRl+I6Pj5ckHTt2zG9/XV2d17imeOqpp3Tdddfp4osv9tvftWtXDR8+3G/f1KlTJUkvvPDCaT8jNjZWCQkJXg8AaC2sVmM7wV/+0nhmqQkAnFlYbjXYo0cPSdKBAwf89h88eFCSlNjEy+lLSkr0/vvva+XKlU16f3p6uiTp22+/1bFjxxqdoQcAAABOFpYz34MGDZKkRu9kuXPnTklSZmZmk47/1FNPKTU1VWPGjGnS+0++2NLS2CX/AAAAwCnCMnxnZ2erffv2Kikp0ZEjR3z6P/jgA0nSTTfdFPSxDxw4oKVLl2rKlCmn3V7w8ssv13PPPee3r6ysTJKx53hjd80EAAAAThWW4Ts+Pl7Tpk1TXV2dXnrpJa++nTt3av369erfv79GjRrV0L5t2zYNHTpUc+bMOe2xn3/+eUnSxIkTTztu3759jS5LWbRokSRpwoQJZ/hJAAAAgBPCMnxL0ty5c3XVVVfpwQcf1IoVK1RTU6OSkhLdcsstio+P16uvvup1a/nCwkIVFxcrPz9f+/fv93tMz/aCd9xxR6PbGHpYLBa99dZbmjRpkr766ivV1dVp586duueee+RwOHTDDTdo1qxZzfozAwAAoHUL2/AdFxenDRs2aObMmZo9e7Y6d+6s66+/Xv3799fmzZs1cOBAr/E5OTnq3LmzcnJyGg3Wq1at0u7du/Wb3/zmjJ///vvva8GCBdqxY4eGDh2q9u3b69JLL9XXX3+txYsXa926dVxoCQAAgKBY3G63O9RFRIvq6mrZbDY5nU62HQQAAIhCYTvzDQAAALQ2hG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkbUNdAAAgMrlcksMhVVVJiYlSVpZ00o2HAQB+EL4BAEGz26Xp06U9e060JSdLCxdKubmhqwsAwh3LTgAAQbHbpbw87+AtSZWVRrvdHpq6ACAScHt5E3F7eQCRzuWSUlN9g7eHxWLMgJeXswQFAPxh5hsAEDCHo/HgLUlut1RRYYwDAPgifAMAAlZV1bzjACDaEL4BAAFLTGzecQAQbQjfAICAZWUZa7otFv/9FouUkmKMAwD4InwDAAJmtRrbCUq+AdzzuqCAiy0BoDGEbwBAUHJzpeXLpaQk7/bkZKOdfb4BoHFsNWgithoE0Jpwh0sACB53uAQANInVKo0YEeoqACCysOwEAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADBJ21AXAACITC6X5HBIVVVSYqKUlSVZraGuCgDCG+EbABA0u12aPl3as+dEW3KytHChlJsburoAINyx7AQAEBS7XcrL8w7eklRZabTb7aGpCwAigcXtdrtDXUS0qK6uls1mk9PpVEJCQqjLAYCguVxSaqpv8PawWIwZ8PJylqAAgD/MfAMAAuZwNB68JcntlioqjHEAAF+EbwBAwKqqmnccAEQbwjcAIGCJic07DgCiDeEbABCwrCxjTbfF4r/fYpFSUoxxAABfhG8AQMCsVmM7Qck3gHteFxRwsSUANIbwDQAISm6utHy5lJTk3Z6cbLSzzzcANC6sw3dtba0effRR9e3bV+3atVNiYqLuuusulZeXB3WcefPmyWKxNPqYMGFCo+/95z//qWuvvVY2m00JCQkaNmyYli9ffpY/GQBEttxcadcuaeNGaelS47m8nOANAGcStuG7pqZG1113nebPn6/HHntMBw4c0L/+9S9t27ZNgwcP1tatW4M6XkJCgvr27ev3kdjIlUEFBQUaNWqUevfura+//lrl5eW64oordOutt2revHnN8FMCQOSyWqURI6Rf/tJ4ZqkJAJxZ2N5k53e/+53+8pe/aNGiRZo6dWpD+7fffqu+ffuqX79+2rp1q6wBnO3nzZunXbt26cUXXwz480tLSzVo0CBdcskl+vjjj2U5aXFjdna2ioqK9OGHH2rIkCEBH5Ob7AAAAES3sJz5Pnr0qJ5++mnFxMRo/PjxXn29evXSyJEjtX37dq1bt67FanjyySflcrl0zz33eAVvSZo0aZIkaf78+S32+QAAAGh9wjJ8b9iwQUeOHNHAgQPVoUMHn/5hw4ZJklatWtViNaxevVqSNHTo0EY/f926dfrpp59arAYAAAC0LmEZvj3rudPS0vz29+zZ02tcIMrLyzV27FilpKQoJiZG5513nm666Sa9++67PmMrKyv1ww8/NFpDUlKSYmJiVFtbq6+//jrgGgAAABDdwjJ8f//995KkLl26+O3v3Lmz17hAfPLJJ7rsssv0ySef6Mcff9TatWt18OBBjRw5Uo899pjfz2/Xrp3i4uJ8jmWxWNSpU6cz1lBXV6fq6mqvBwAAAKJXWIbvo0ePSpJiYmL89sfGxnqNO5OxY8dq06ZNmjFjhhITExUTE6PLLrtMa9asUefOnfWHP/xBH3zwQcCfH2gN+fn5stlsDY+UlJSA6gUAAEDrFJbhOz4+XpJ07Ngxv/11dXVe486kT58+GjBggE+7zWbTHXfcIUlasmRJwJ8faA2zZ8+W0+lseFRUVARULwAAAFqntqEuwJ8ePXpIkg4cOOC3/+DBg5LU6P7cwUhPT5ckffnllz6fX1tbq5qaGp+lJ263W4cOHTpjDbGxsQ0z5AAAAEBYznwPGjRIkhq9k+XOnTslSZmZmWf9WfX19ZLktZ1gUlKSunXr1mgNlZWVOnbsmOLi4tS7d++zrgEAAADRISzDd3Z2ttq3b6+SkhIdOXLEp9+zPvumm24647EqKirUo0cPlZSU+O0vKyuTJPXr18+rfcyYMZKk4uLiRj//xhtvVNu2YfnHAwAAAIShsAzf8fHxmjZtmurq6vTSSy959e3cuVPr169X//79NWrUqIb2bdu2aejQoZozZ47XeJfLpf/85z9+b8hz8OBBLV26VJI0YcIEr76ZM2fKarXqr3/9q069CejixYslSQ888ECTf0YAAABEn7AM35I0d+5cXXXVVXrwwQe1YsUK1dTUqKSkRLfccovi4+P16quvet1avrCwUMXFxcrPz9f+/fsb2j3LSR555BEVFBRoz549OnbsmDZt2qTRo0fL6XRq3rx5DTfO8ejfv78WLFigjz/+WFOmTNHevXt14MAB3X///dqwYYPmzp0b1K3lAQAAAIv71GndMFJbW6s///nPevXVV7V7927ZbDZdf/31evTRRxtutOOxfv163XbbbRo+fLjsdrvXGu7PPvtMr7/+ut5++2198803qq2tVbdu3TRkyBBNnTpVI0eObLSGdevWaf78+frss8/kdruVmZmp3/72t8rLywv656murpbNZpPT6VRCQkLQ7wcAAEBkC+vw3doQvgEAAKIbVwsCAJrE5ZIcDqmqSkpMlLKypJNWAwIA/CB8AwCCZrdL06dLe/acaEtOlhYulHJzQ1cXAIS7sL3gEgAQnux2KS/PO3hLUmWl0W63h6YuAIgErPk2EWu+AUQ6l0tKTfUN3h4WizEDXl7OEhQA8IeZbwBAwByOxoO3JLndUkWFMQ4A4IvwDQAIWFVV844DgGhD+AYABCwxsXnHAUC0IXwDAAKWlWWs6T7pPmZeLBYpJcUYBwDwRfgGAATMajW2E5R8A7jndUEBF1sCQGMI3wCAoOTmSsuXS0lJ3u3JyUY7+3wDQOPYatBEbDUIoDXhDpcAEDzucAkAaBKrVRoxItRVAEBkYdkJAAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgEsI3AAAAYBLCNwAAAGASwjcAAABgkrahLgAAEJlcLsnhkKqqpMREKStLslpDXRUAhDfCNwAgaHa7NH26tGfPibbkZGnhQik3N3R1AUC4Y9kJACAodruUl+cdvCWpstJot9tDUxcARAKL2+12h7qIaFFdXS2bzSan06mEhIRQlwMAQXO5pNRU3+DtYbEYM+Dl5SxBAQB/mPkGAATM4Wg8eEuS2y1VVBjjAAC+CN8AgIBVVTXvOACINoRvAEDAEhObdxwARBvCNwAgYFlZxppui8V/v8UipaQY4wAAvgjfAICAWa3GdoKSbwD3vC4o4GJLAGgM4RsAEJTcXGn5cikpybs9OdloZ59vAGgcWw2aiK0GAbQm3OESAILHHS4BAE1itUojRoS6CgCILCw7AQAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMEtbhu7a2Vo8++qj69u2rdu3aKTExUXfddZfKy8uDOs6mTZv0m9/8RgMGDFB8fLzatWun9PR0TZ48udFjFRUVyWKxNPpITU1thp8QAAAA0SRstxqsqanRyJEjVVJSohdeeEE///nPVVZWpvHjx2vw4MF67733NHDgwDMeZ+3atRo9erQuvPBCLVy4UNdcc42OHz+uVatWacaMGXrllVf01ltv6eqrr/Z5b9u2bdWrVy+/x0069e4SAAAAwBmEbfieN2+eiouLtWjRIt16662SpIEDB2rFihXq27ev7rjjDm3dulXWM9zRoaamRpL02muv6corr2xo/9WvfqWffvpJkydP1sSJE/Xll1/6vDcpKUk7duxoxp8KAAAA0Swsl50cPXpUTz/9tGJiYjR+/Hivvl69emnkyJHavn271q1bd8ZjdenSRddff71X8Pa46aabJEk7duzQ7t27m6d4AAAAoBFhGb43bNigI0eOaODAgerQoYNP/7BhwyRJq1atOuOxrr32Wr399tt++2w229kVCgAAAAQhLMP31q1bJUlpaWl++3v27Ok1rqk8S0p69uypCy64wKe/pqZGs2bNUp8+fRQXF6dOnTpp2LBhWrx4serr68/qswEAABB9wjJ8f//995KMJSP+dO7c2WtcU7388suSpIcffthv/969e+V0OrV69Wo5nU5t3rxZ/fv316RJkzRmzBgdP378tMevq6tTdXW11wMAAADRKyzD99GjRyVJMTExfvtjY2O9xjVFSUmJnn32Wd1yyy0+68olKSMjQ0VFRSosLFTfvn0VExOjtLQ0FRYWKjs7W+vWrdP8+fNP+xn5+fmy2WwNj5SUlCbXCwAAgMgXluE7Pj5eknTs2DG//XV1dV7jgrVv3z7deuutuuKKKxpmv0/VtWtXDR8+3G/f1KlTJUkvvPDCaT9n9uzZcjqdDY+Kioom1QsAAIDWISy3GuzRo4ck6cCBA377Dx48KElKTEwM+tgHDx7Uz372M3Xv3l1r165tUoBPT0+XJH377bc6duzYaWfoPbP0ANDauFySwyFVVUmJiVJWlnSG3V8BIOqF5cz3oEGDJKnRu0/u3LlTkpSZmRnUcfft26fs7Gx17txZb7/9dpN3Ozn5YkuLxdKkYwBAJLPbpdRUKTtbGjvWeE5NNdoBAI0Ly/CdnZ2t9u3bq6SkREeOHPHp/+CDDySd2Kc7EP/+9781fPhwJSUlac2aNV4z3s8//7y2b9/uNf7yyy/Xc8895/dYZWVlkow9x88555yAawCA1sBul/LypD17vNsrK412AjgANC4sw3d8fLymTZumuro6vfTSS159O3fu1Pr169W/f3+NGjWqoX3btm0aOnSo5syZ43O83bt365prrtFFF12kN954Q+3atfPq/9Of/qRNmzZ5te3bt08rV670W9+iRYskSRMmTGjCTwcAkcvlkqZPl9xu3z5P2333GeMAAL7Ccs23JM2dO1fvv/++HnzwQZ133nkaNWqUysrKNH78eMXHx+vVV1/1urV8YWGhiouLVVxcrJkzZ6pr166SjLCenZ2tyspKDRo0SOPGjfP5rH379vm0WSwWvfXWW5o0aZJmzJih1NRUVVZWKj8/Xw6HQzfccINmzZrVcv8CACAMORy+M94nc7uligpj3IgRppUFABEjbMN3XFycNmzYoD//+c+aPXu27rjjDtlsNl1//fVasWJFw412PHJycrR06VINHz7ca3/wVatWNdw6fsWKFQF//vvvv6/XX39dq1at0tChQ3Xo0CF17NhRmZmZWrx4sSZOnKg2bcLyDwcA0GKqqpp3HABEG4vb7e+Ph2gJ1dXVstlscjqdSkhICHU5ABC0oiLj4soz2biRmW8A8IepWwBAwLKypORkqbGNniwWKSXFGAcA8EX4BgAEzGqVFi40/vnUAO55XVDAft8A0BjCNwAgKLm50vLlUlKSd3tystGemxuaugAgErDm20Ss+QbQmnCHSwAIXtjudgIACG9WKxdVAkCwWHYCAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJikbagLAABEJpdLcjikqiopMVHKypKs1lBXBQDhjfANAAia3S5Nny7t2XOiLTlZWrhQys0NXV0AEO5YdgIACIrdLuXleQdvSaqsNNrt9tDUBQCRwOJ2u92hLiJaVFdXy2azyel0KiEhIdTlAEDQXC4pNdU3eHtYLMYMeHk5S1AAwB9mvgEAAXM4Gg/ekuR2SxUVxjgAgC/CNwAgYFVVzTsOAKIN4RsAELDExOYdBwDRhvANAAhYVpaxptti8d9vsUgpKcY4AIAvwjcAIGBWq7GdoOQbwD2vCwq42BIAGkP4BgAEJTdXWr5cSkrybk9ONtrZ5xsAGsdWgyZiq0EArQl3uASA4HGHSwBAk1it0ogRoa4CACILy04AAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAk7QNdQHRxO12S5Kqq6tDXAkAAABOp2PHjrJYLM1+XMK3ifbv3y9JSklJCXElAAAAOJ29e/eqe/fuzX5cwreJunTpIknavXu3bDZbiKtBa1ZdXa2UlBRVVFQoISEh1OWgFeO7BrPwXYNZPN+1mJiYFjk+4dtEbdoYS+xtNhsnDpgiISGB7xpMwXcNZuG7BrO0xJITiQsuAQAAANMQvgEAAACTEL5NFBsbq7lz5yo2NjbUpaCV47sGs/Bdg1n4rsEsLf1ds7g9+98BAAAAaFHMfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDrUhNTY0efvhh9enTR+3atdP555+vu+++W5WVlaEuDQAAiPBtik8//VTjxo1Tenq6LBaLfv/734e6JLRCtbW1uvbaa/XHP/5RP/74o26++WalpKRoyZIlGjx4sHbu3BnqEhFBOG8BaG1eeuklXXbZZerUqZPat2+vSy65RH//+99Nr4Pby5vggw8+0EcffaRhw4bphx9+CHU5aKUee+wxffTRR7rqqqv09ttvq0OHDpKkJ598UjNnztTdd9+toqKi0BaJiMF5C0Brc/DgQeXk5GjQoEFq166dVq5cqV/+8pdq166dcnJyTKuDfb5NUF9frzZtjD8ypKam6s4779Rjjz0W4qrQmhw7dkznnnuunE6nPv/8cw0ePNirf+DAgSopKdGnn36qSy+9NERVIpJw3kKoffrpp7r88ssDGjty5Ei98847LVwRWqNhw4YpMTFR//jHP0z7TGa+TeD5HxjQUj744AM5nU716tXLJ3hLUl5enkpKSrR69WrCNwLCeQuhVlNTozvvvFOvvPKKOnTooFtuucWr/7XXXlNtba3GjRunG2+8MURVItJ17dpVx48fN/UzCd9AK7B161ZJ0iWXXOK339NeUlJiWk0AcDaysrLUrVs3vfLKK7rsssv04osvNvT9+OOP+tvf/qbk5GS99NJLoSsSEemnn37S0aNH9c9//lPvvPOOVqxYYernE76BVmD37t2SpOTkZL/9nvbvvvvOtJoA4Gxt27ZNkpSRkeHVXlJSovr6emVmZoaiLESw77//XomJiZIkq9WqZ555xvS/nBC+gVbgxx9/lCTFx8f77W/fvr0k6fDhw6bVBABny/PXulND9pYtW/y2A2fSrVs3bdq0SYcPH9a//vUvTZs2TV27dvVZ1tSSCN8AACAsecL3qTPfhG80Vdu2bXXZZZdJkrKzs3XgwAHNnj3b1PDNFTVAK+DZVvDo0aN++48cOSJJ6tixo2k1AcDZKikpkcVi0YABA7zaPeF74MCBIagKrcmgQYNMvw8G4TsMffrpp7JYLAE9rr/++lCXizBwwQUXSJL27Nnjt9/TfuGFF5pWE6IH5yy0hOrqan333Xfq2bNnw9I5SXK5XCotLVVsbKz69u0bwgrRGhQXFys1NdWnvSXPayw7McG+ffv03nvvSTJmJnfs2KHly5erffv2fhf5s70SguWZ/fn888/99nva+RMtAhXMeYtzFlpCY+u9y8rKVFNTo8GDB8tqtYaiNESo7Oxs3XLLLerXr59qa2v15ptvaunSpVq8eLHP2BY9r7nR4jZu3OiW5PO48MILG33PF1984ZbkHjFihFf74cOH3W3atHEnJye3cNWIJHV1dW6bzeaW5N68ebNPf2ZmpluS+9NPPzW/OESkYM9bnLPQ3J5++mm3JPfDDz/s1f7GG2+4Jbn/67/+K0SVIVJNnz7d3adPH3dcXJy7W7du7muuuca9evXqRse31HmNZScmGDFihNxut89j165djb6H7ZUQjJiYGE2bNk2SNHXq1IY13pJxe/mSkhINHz6cG+wgYMGetzhnobk1NvPt2d2ppKREa9euVU1Njem1ITIVFBToq6++0tGjRxv+ujd69OhGx7fUeY1lJwH4xS9+oS+//DKo97z88su64oormvyZbK+EYP3+97/X+vXrVVxcrN69eysrK0vfffedPv74Y3Xv3l0vvPBCqEtEK8Y5C82tseBz7bXXatCgQdq+fbtmzJihL774IhTlIQq01HmN8B2A8vJyffXVV0G9p7FdJwLF9koIVrt27bRx40bl5+dr6dKlWrlypbp06aIJEyboj3/8Y6M34AGaA+csNLcPPvjAb/v555+vzZs3m1wNolFLndcsbrfbfVaVoUWkpqZq9+7dOnz4sNdV3ldccYU2bdqk7du36+KLLw5hhQBwAucsAK1NS53XCN8tzGKxBDzW85+iurpaNptNvXr10jfffNPQ73K51LFjR9XX1+vIkSNc5Q2gRQR73uKcBSDchdN5jQsuW5jb7dZrr72mUaNGqUePHkpISFBWVpYcDofPhUweZ9pe6eKLL+Z/YgBaTLDnLc5ZAMJdOJ3XCN8mKCgoULdu3fT000/rH//4h5KSknTddddp69atfsc3tsZox44dksRNBQC0uGDOW5yzAESCcDmvccGlCVavXq2uXbs2vB45cqQyMjL09NNP+93YPdDtla699lrFxcW1YOUAolUw5y3OWQAiQbic15j5NsHJ/6ElqU2bNhowYIDKy8v9jj/T9kplZWWaMWOGYmJiWqZgAFEvmPMW5ywAkSBczmtccBkCLpdLvXv31g033KBnn3021OUAwBlx3gLQ2oTqvMbMdwgsWrRIu3fv1pQpU0JdCgAEhPMWgNYmVOc1wrfJPv74Yz344IP6/e9/7/OnDAAIR5y3ALQ2oTyvsezERLt27dKQIUN0zTXX6LXXXgtqz0kACAXOWwBam1Cf1wjfJjl06JCGDh2qjh07auPGjVzxDyDscd4C0NqEw3mNrQZNcOzYMeXm5uro0aPasGED/wMDEPY4bwFobcLlvEb4NsGUKVP03nvv6fnnn1d5eXnDljaxsbEaPHhwiKsDAF+ctwC0NuFyXmPZiQlSU1P13Xff+bRfeOGF2rVrl/kFAcAZcN4C0NqEy3mN8A0AAACYhK0GAQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAEGrqqrS1KlTNWnSJB07dkzTp0/Xeeedp27duukPf/iD6uvrQ10iAIQlwjcAICjvvPOO+vfvr2eeeUYXX3yxHnnkET311FPq3bu36urq9Nhjj+mxxx4LdZkAEJYsbrfbHeoiAACR4ZtvvtGgQYNUW1urV199VXl5eTrvvPMUGxur3bt3a+XKlcrLy5PNZtOBAwfUpg1zPABwMs6KAICA/e53v9ORI0d0//336/bbb9eWLVu0f/9+XXPNNbJarcrNzVWnTp3kdDpVUVER6nIBIOwQvgEAAfnPf/6jN998Ux06dNDcuXMlSZ9//rkkadCgQZIki8Wi7t27S5JcLldI6gSAcEb4BgAEZOPGjXK5XLrhhhsUFxcnSdqyZYskaeDAgQ3j9u/fL0lKTEw0vUYACHeEbwBAQLZu3SpJuvLKKxvaPOHbM/NdUVGhAwcOqF+/fg0BHQBwAuEbABCQH374QZLUqVMnSZLb7VZJSYnOO+889ejRQ5K0du1aSdKoUaNCUiMAhDvCNwAgIB07dpQk7dmzR5JUVlamH3/8sWHJSW1trZ588kmdc845uvfee0NWJwCEM8I3ACAgw4YNkyS9/PLLcjqdXktOnE6nbr31VpWVlWnOnDnq3bt3CCsFgPDFPt8AgIC4XC5dddVV2rRpk1JSUtSlSxdt3bpVgwcPVkVFhX744QdNnjxZTz/9tCwWS6jLBYCwRPgGAARs//79euCBB/Tmm2/qwIEDkqT27dvr0ksv1W9/+1vl5OSEtkAACHOEbwBAk3Tr1k3Hjx/XoUOHmOkGgACx5hsAELTvvvtO+/fv1yWXXELwBoAgEL4BAEHz3Nny0ksvDXElABBZCN8AgKARvgGgaVjzDQAAAJiEmW8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJIRvAAAAwCSEbwAAAMAkhG8AAADAJP8P/9B+6AliO34AAAAASUVORK5CYII=",
      "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(r'$\\varphi$')\n",
    "plt.ylabel(r'$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = [r'$-\\frac{1}{2}\\pi$', r'$0$', r'$\\frac{1}{2}\\pi$', r'$\\pi$', r'$\\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": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAInCAYAAABAyzWwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKyUlEQVR4nO3de3STVaL+8ScECy3QcFVqW22hCArloqgoFCjicmRAO7XqGVRg4TiIOAsHlCM4I3gZu2ZAp3DwUpwj6iiMDkTkNqMioD0UFRUoRdEqRUqnCAKmCG3BNL8/3l8KISkkvby5fT9rZcXsvfNmp6Q7j7v73a/F5XK5BAAAAKDZtQh2BwAAAIBoQfgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvg2kcvlUmVlpdhaHQAAIDoRvk109OhR2Ww2HT16NNhdAQAAQBAQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJO0DHYHAADhyemUCgqkigopIUHKyJCs1mD3CgBCGzPfAICA2e1SWpqUmSmNHWvcp6UZ5QCA+hG+AQABsdulnBwpPV3avFk6etS4T083ygngAFA/i8vlcgW7E9GisrJSNptNDodD8fHxwe4OAATM6TRmuNPTpRUrpBanTeHU1kpZWVJxsVRSwhIUAPCFNd9AELBWFuGqoEDas0dautQzeEvG45kzpWuvNdoNHx6MHgKBYTyG2QjfgMnsdmnqVGnfvlNlSUnS/PlSdnbw+gX4o6LCuO/Tx3e9u9zdDghljMcIBtZ8AyZyr5U9faCXpPJy1soiPCQkGPfFxb7r3eXudkCoYjxGsLDm20Ss+Y5uTqeUkuI90LtZLMaMS2kpf/JE6GLNNyIB4zGCiZlvwCQFBfUP9JLkckllZUY7IFRZrdLTT0urVxtB+/TdTrKyjPJ58wgsCG2Mxwgm1nwDJvF3DSxrZRHqsrOlZcuk6dONkyvdUlONctbKItQxHiOYCN+ASfxdA8taWYSD7Gzp5pvZJQLhifEYwcSabxOx5ju6udcYlpcbf9I8E2sMAcAcjMcIJtZ8AyaxWo3tqyRjYD+d+3FeHgM9ADQ3xmMEE+EbMJF7rWxiomd5UhJrZQHATIzHCBaWnZiIZSdw44pqABAaGI9hNsK3iQjfAAAA0Y1lJwAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJCN8AAACASQjfAAAAgEkI3wAAAIBJWga7AwCA8OR0SgUFUkWFlJAgZWRIVmuwewUAoS3kZ76rq6v1yCOPKCYmRhMmTAj4+Xv27JHFYjnnLSUlxeN5GzduDKg9AEQTu11KS5MyM6WxY437tDSjHABQv5Ce+d6wYYMmTZqkH374QSdPnmzUsbp166bzzjvPZ11JSYkuvfRSr/KWLVuqe/fuPp+TmJjYqP4AQLiy26WcHGn0aGnpUqlPH6m4WHrqKaN82TIpOzvYvQSA0BSy4Xvp0qWaPHmynnrqKcXGxmrixImNOt7777/vc7Z627ZtGjBggO6//36vusTERO3atatRrwsAkcTplKZPN4L3ihVSi///99NBg4zHWVnSgw9KN9/MEhQA8CVkw3dqaqq++OILXXjhhXr55ZcbfJw2bdro9ttvV5s2bXzWL1iwQN27d9eNN97Y4NcAmhpraRGqCgqkPXuMGe8WZyxcbNFCmjlTuvZao93w4cHoIeCJ8RShJmTD96BBg5rkOF26dNE//vEPn3WHDh3S0qVL9dRTT6nFmd8iQJDY7dLUqdK+fafKkpKk+fP5Uz6Cr6LCuO/Tx3e9u9zdDggmxlOEoqhOnC+++KKsVmu9S1qqqqo0Y8YMXXLJJYqNjVX79u01ZMgQLVq0SLW1tSb3FtHAvZb29C8KSSovN8o5mQ3BlpBg3BcX+653l7vbAcHCeIpQFbXh2+l06rnnntO4ceNks9l8tjlw4IAcDodWrVolh8OhrVu3qnfv3po0aZLGjBlzzpNAa2pqVFlZ6XED6uN0GjM0Lpd3nbvsgQeMdkCwZGRIKSnGyZVnzkHU1kq5uVJqqtEOCBbGU4SyqA3fK1asUFlZmc8TLSUpPT1dGzduVH5+vnr27KmYmBilpqYqPz9fmZmZWrt2rebOnXvW18jNzZXNZqu7JScnN8dbQYQoKPCeoTmdyyWVlRntgGCxWqWnn5ZWrzZOrty8WTp61LjPyjLK581jTS2Ci/EUoSxqw/eCBQt03XXX6bLLLvNZ36lTJw0bNsxn3ZQpUyRJL7300llfY+bMmXI4HHW3srKyxnUaEc3fNbKspUWwZWcb2wnu2GGcXBkfb9wXF7PNIEID4ylCWciecNmcioqK9OGHH2rFihUNen5aWpok6dtvv9WJEycUExPjs12rVq3UqlWrhnYTUcbfNbKspUUoyM42thNkFwmEIsZThLKoDN8LFixQSkqKxowZ06Dnn36ypcViaapuIcplZBhn4ZeX+16naLEY9aylRaiwWtlOEKGJ8RShLOqWnRw+fFhLlizRfffdd9btBa+88kq98MILPutKSkokSd27d6/3qplAoKxWY/sryfhiOJ37cV4eM4sAcC6MpwhlERO+d+zYocGDB2vWrFlnbffiiy9Kku6+++6ztjt48GC9y1IWLlwoSZowYULA/QTOxr2WNjHRszwpibW0ABAIxlOEqogJ3/n5+SosLFRubq4OHTrks417e8E77rhDHTt2POvxLBaL3nnnHU2aNElfffWVampqtHv3bt1zzz0qKCjQDTfcoBkzZjTHW0GUy842riC4YYO0ZIlxX1rKFwUABIrxFKHI4nL5Wg0VGs62nnrx4sUeM8/r1q3TbbfdpmHDhslut/t87ltvvaXs7Gxt375dffv2Petrl5WV6c0339TKlSu1c+dO/fjjj2rXrp369u2rO++8U3fffXfAV8WsrKyUzWaTw+FQfHx8QM8FAABA+Avp8B1pCN8AAADRLWKWnQAAAAChjvANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYpGWwOwAACE9Op1RQIFVUSAkJUkaGZLUGu1cAENqY+QYABMxul9LSpMxMaexY4z4tzSgHANSP8A2EKadT2rhRWrrUuHc6g90jRAu7XcrJkdLTpc2bpaNHjfv0dKOcAA4zMRYi3FhcLpcr2J2IFpWVlbLZbHI4HIqPjw92dxDG7HZp6lRp375TZUlJ0vz5UnZ28PqFyOd0GjPc6enSihVSi9OmcGprpawsqbhYKilhCQqaH2MhwhEz30CYcc86nv5lI0nl5cw6ovkVFEh79kizZnkGb8l4PHOmVFpqtAOaE2MhwhXhGwgjTqcxy+Pr71Xusgce4M+uaD4VFcZ9nz6+693l7nZAc2AsRDgjfANhpKDAe5bndC6XVFbGrCOaT0KCcV9c7LveXe5uBzQHxkKEM8I3EEb8nU1k1hHNJSNDSkmRnnrKWON9utpaKTdXSk012gHNhbEQ4YzwDYQRf2cTmXVEc7FapaefllavNk6uPH23k6wso3zePE62RPNiLEQ4Y7cTE7HbCRrL6TRmHcvLfa91tFiMM/1LSwk/aF52uzR9unHypVtqqhG82WUCzY2xEOGMmW8gjFitxhZakvHlcjr347w8vmzQ/LKzpW++kTZskJYsMe5LSgjeMAdjIcIZM98mYuYbTcXX3rbJycaXDeEHQLRgLEQ4InybiPCNpuR0GmfyV1QY6xozMpjlARB9GAsRbgjfJiJ8AwAARDfWfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACZpGewOAADCk9MpFRRIFRVSQoKUkSFZrcHuFQCENsI3EOEISGgOdrs0fbq0Z8+pspQU6emnpezsYPUKkYSxC5GKZSdABLPbjUCUmSmNHWvcp6QY5UBD2e1STo6Uni5t3iwdPWrcp6cb5Xy+0FiMXYhkFpfL5Qp2J6JFZWWlbDabHA6H4uPjg90dRDh3QDrzN9xiMe6XLWOGEoFzOqW0NCNor1ghtThtCqe2VsrKkoqLpZISZinRMIxdiHSEbxMRvmEWp9OYJdq3z3e9xSIlJUmlpQQkBGbjRmMWcvNmadAg7/rNm6Vrr5U2bJCGDze7dwh3jF2IBiw7ASJQQUH9X16SMaNUVma0AwJRUWHc9+nju95d7m4HBIKxC9GA8A1EIH+DDwEJgUpIMO6Li33Xu8vd7YBAMHYhGhC+gQjkb/AhICFQGRnGsoCnnjLWeJ+utlbKzZVSU412QKAYuxANQj58V1dX65FHHlFMTIwmTJjQoGPMmTNHFoul3tvZjvuvf/1LI0aMkM1mU3x8vIYMGaJly5Y17M0AJsnIMNZFuk9QOpPFIiUnE5AQOKvV2E5w9Wrj5MrTdzvJyjLK581jPS4ahrEL0SCkw/eGDRvUt29fPf/88zp58mSjjhUfH6+ePXv6vCXU87/QeXl5GjVqlHr06KGvv/5apaWluuqqq3Trrbdqzpw5jeoP0JysVmn+fOO/z/wScz/OyyMgoWGys40dJ3bsME6ujI837ouL2YkCjcPYhWgQsrudLF26VJMnT9ZTTz2l2NhYTZw4UePHj9fLL78c8LHmzJmjPXv2BPTc4uJi9e/fX5dffrk+/vhjWU4bBTIzM7Vx40Zt3rxZg3yd7l8PdjuB2ex2aepUzxOYkpONLy8CEhqLi6CguTB2IZKF7BUuU1NT9cUXX+jCCy9sUOBurGeeeUZOp1P33HOPR/CWpEmTJmnjxo2aO3euli9fbnrfAH9lZ0s330xAQvOwWtlOEM2DsQuRLGTDdyAzys1h1apVkqTBgwd71Q0ZMkSStHbtWv38889q2TJkf4wAAQlAWGLsQqQK6TXfTam0tFRjx45VcnKyYmJidMEFF+imm27S+++/79W2vLxcP/zwgyRjBv5MiYmJiomJUXV1tb7++utm7zsAAAAiQ9SE708++UQDBw7UJ598op9++klr1qzRkSNHNHLkSD355JMebffv3y9Jat26tWJjY72OZbFY1L59e4+2vtTU1KiystLjBgAAgOgVFeF77Nix2rJli6ZNm6aEhATFxMRo4MCBWr16tTp06KA//vGP2rRpU13748ePS5JiYmLqPWarVq082vqSm5srm81Wd0tOTm6idwQAAIBwFBXh+5JLLlEfH9dCttlsuuOOOyRJixcvriuPi4uTJJ04caLeY9bU1Hi09WXmzJlyOBx1t7Kysgb1HwAAAJEh6s8UTEtLkyR9+eWXdWVdu3aVZFzgp6qqymvpicvl0o8//ihJ9e4RLhmz4+4ZcgAAACAqZr7Ppvb/Xx/59O0EExMT1blzZ0nGiZpnKi8v14kTJxQbG6sePXqY01EAAACEvYgP32VlZeratauKiop81peUlEiSevXq5VE+ZswYSVJhYaHXc9zrw2+88Ua2GQQAAIDfIiZ879ixQ4MHD9asWbM8yp1Op77//nutXbvW6zlHjhzRkiVLJEkTJkzwqJs+fbqsVqv+9re/6cyLgC5atEiS9NBDDzXhOwAAAECki5jwnZ+fr8LCQuXm5urQoUN15e7lJI899pjy8vK0b98+nThxQlu2bNHo0aPlcDg0Z86cugvnuPXu3Vvz5s3Txx9/rPvuu08HDhzQ4cOH9eCDD2r9+vWaPXt20C8EBAAAgPBicZ05rRtCzrys++kWL17sMVu9bt063XbbbRo2bJjsdrvHcz/77DO9+eabevfdd/XNN9+ourpanTt31qBBgzRlyhSNHDmy3tdZu3at5s6dq88++0wul0t9+/bV73//e+Xk5AT8fiorK2Wz2eRwOBQfHx/w8wEAABDeQjp8RxrCNwAAQHSLmGUnAAAAQKgjfAMAAAAmYZ88AB6cTqmgQKqokBISpIwMyWoNdq8QivisIBB8XgAD4RtAHbtdmjpV2rfvVFlSkjR/vpSdHbx+IfTY7dL06dKePafKUlKkp5/mswJvjC3AKSw7ASDJ+HLMyfH8cpSk8nKj3G4PTr8QetyflfR0afNm6ehR4z49nc8KvDG2AJ7Y7cRE7HaCUOV0GrOWZ345ulksxixVaSl/Jo52TqeUlmYE7RUrpBanTeHU1kpZWVJxsVRSwmcFjC2AL8x8A1BBQf1fjpLkckllZUY7RLeCAmOpyaxZnsFbMh7PnGkEKT4rkBhbAF8I3wBUUdG07RC53J+BPn1817vL+axAYmwBfCF8A1BCQtO2Q+RyfwaKi33Xu8v5rEBibAF8Yc23iVjzjVDlXpdZXm78GfhMrMuEG2u+EQjGFsAbM98AZLUaW35Jxpfh6dyP8/L4coTxGXj6aWn1aiNon77bSVaWUT5vHp8VGBhbAG+EbwCSjL12ly2TEhM9y5OSjHL24oWb+7OyY4d07bVSfLxxX1zMZwXeGFsATyw7MRHLThAOuAod/MVnBYHg8wIYCN8mInwDAABEN5adAAAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJmkZ7A4ACH9Op1RQIFVUSAkJUkaGZLUGu1cAmgO/70DjEL4BNIrdLk2dKu3bd6osKUmaP1/Kzg5evwA0PX7fgcazuFwuV7A7ES0qKytls9nkcDgUHx8f7O4AjWa3Szk50pmjiMVi3C9bxhdyJGMGNLrw+w40DcK3iQjfiCROp5SS4jkDdjqLxZgRKy0lkEUiu12aPl3as+dUWUqK9PTTBLBIxO870HQ44RJAgxQU1P9FLBmzY2VlRjtEFvcMaHq6tHmzdPSocZ+ebpTb7cHuIZoav+9A0yF8A2iQioqmbYfw4HQaM96jR0srVkiDBklt2xr3K1YY5Q8+aLRD5OD3HWg6hG8ADZKQ0LTtEB4KCoylJrNmSS3O+AZp0UKaOdNYesAMaGTh9x1oOoRvAA2SkWGs8XSfbHUmi0VKTjbaIXK4Zzb79PFd7y5nBjSy8PsONB3CN4AGsVqN7cUk7y9k9+O8PE6+ijTumc3iYt/17nJmQCMLv+9A0yF8A2iw7Gxje7HERM/ypCS2HYtUGRnGrhdPPSXV1nrW1dZKublSaiozoJGI33egabDVoInYahCRiv2eo4t7t5PRo4013n36GDPeubnS6tUEsUjH7zvQOIRvExG+AUQKX/t8p6ZK8+YRvAHgbAjfJiJ8A4gkzIACQOBaBrsDAIDwZLVKw4cHuxcAEF444RIAAAAwCeEbAAAAMAnhGwAAADBJyIfv6upqPfLII4qJidGECRMadIwtW7bod7/7nfr06aO4uDi1bt1aaWlpmjx5skpLS30+Z+PGjbJYLPXeUlJSGv6mAAAAEJVC+oTLDRs2aNKkSfrhhx908uTJBh1jzZo1Gj16tC6++GLNnz9fQ4cO1cmTJ7Vy5UpNmzZNr732mt555x1de+21Xs9t2bKlunfv7vO4iWdeZQAAAAA4h5AN30uXLtXkyZP11FNPKTY2VhMnTmzQcaqqqiRJb7zxhq6++uq68t/85jf6+eefNXnyZN1999368ssvvZ6bmJioXbt2NewNAAAAAGcI2WUnqamp+uKLL3TffffJYrE0+DgdO3bU9ddf7xG83W666SZJ0q5du7R3794GvwYAAADgj5Cd+R40aFCTHGfEiBEaMWKEzzqbzdYkrwEAAAD4I2Rnvs3gXlLSrVs3XXTRRV71VVVVmjFjhi655BLFxsaqffv2GjJkiBYtWqTa2lqzuwtENKdT2rhRWrrUuHc6g90jIPLwewYEX1SH71dffVWS9Oijj/qsP3DggBwOh1atWiWHw6GtW7eqd+/emjRpksaMGXPOk0BrampUWVnpcQPgzW6XUlKkzExp7FjjPiXFKAfQNPg9A0JD1IbvoqIiPf/887rllls0fvx4r/r09HRt3LhR+fn56tmzp2JiYpSamqr8/HxlZmZq7dq1mjt37llfIzc3Vzabre6WnJzcXG8HCFt2u5STI+3b51leXm6UEwyAxuP3DAgdURm+Dx48qFtvvVVXXXVV3ez3mTp16qRhw4b5rJsyZYok6aWXXjrr68ycOVMOh6PuVlZW1riOAxHG6ZSmTpVcLu86d9kDD/CncaAx+D0DQkvUhe8jR47oF7/4hbp06aI1a9YoLi4u4GOkpaVJkr799ludOHGi3natWrVSfHy8xw3AKQUF3jNxp3O5pLIyox2AhuH3DAgtURW+Dx48qMzMTHXo0EHvvvtug3c7Of1ky8ZsgwhEu4qKpm0HwBu/Z0BoiZrw/Z///EfDhg1TYmKiVq9e7THj/eKLL2rnzp0e7a+88kq98MILPo9VUlIiSerevbvOO++85us0EOESEpq2HQBv/J4BoSViwveOHTs0ePBgzZo1y6tu7969Gjp0qC699FK99dZbat26tUf9n/70J23ZssWj7ODBg1qxYoXP11q4cKEkacKECU3SdyBaZWRISUlSfX9Aslik5GSjHUIP29aFB37PgNASshfZCVR+fr4KCwtVWFio6dOnq1OnTpKk3bt3KzMzU+Xl5erfv7/GjRvn9dyDBw96lVksFr3zzjuaNGmSpk2bppSUFJWXlys3N1cFBQW64YYbNGPGjGZ/X0Aks1ql+fON3RYsFs8TwtxBIS/PaIfQYrdL06dLe/acKktJkZ5+WsrODlav4Au/Z0CIcYUwSfXeFi9e7NH2vffec3Xo0MGVlZXlqq2trSv/61//etbj1He8vXv3uubNm+caOnSoq1OnTi6r1epq3769a+jQoa5Fixa5nE5nwO/H4XC4JLkcDkdDfhxAxFq+3OVKSnK5jFhg3JKTjXKEnuXLXS6LxeUaM8bl2rzZ5Tp61LgfM8Yo598tNPF7BoQGi8vla/MhNIfKykrZbDY5HA52PgHO4HQauy1UVBhrTzMymIkLRU6nlJYmpadLK1ZILU5bvFhbK2VlScXFUkkJ/36hiN8zIPgI3yYifAMIdxs3GldG3LxZGjTIu37zZunaa6UNG6Thw83uHQCEvog54RIA0Pzc29H16eO73l3OtnUA4BvhGwDgN/d2dMXFvuvd5WxbBwC+sezERCw7ARDuWPMNAI3DzDcAwG9Wq7Gd4OrVRtDevFk6etS4z8oyyufNI3gDQH0iZp9vAIA5srOlZcuMfb6vvfZUeWqqUc4+3wBQP5admIhlJwAiCdvWAUDgmPkGADSI1cp2ggAQKNZ8AwAAACYhfAMAAAAmIXwDAAAAJmHNN4CIxQmBCHd8hoHIQ/gGEJHsdmnqVGnfvlNlSUnS/PlshYfwwGcYiExsNWgithoEzGG3Szk50pmjm8Vi3LMXNUIdn2EgchG+TUT4Bpqf0ymlpHjOFp7OYjFmD0tL+fM9QhOfYSCyccIlgIhSUFB/aJGMmcSyMqMdEIr4DAORjfANIKJUVDRtO8BsfIaByEb4BhBREhKath1gNj7DQGRjzbeJWPMNND/3etnycu+T1STWyyL08RkGIhsz3wAiitVqbMUmndoZws39OC+P0ILQxWcYiGyEbwARJzvb2IotMdGzPCmJLdoQHvgMA5GLZScmYtkJYC6uDohwx2cYiDyEbxMRvgFEEoIhAASOZScAgIDZ7VJampSZKY0da9ynpRnlAID6Eb4BAAFxX/o8PV3avFk6etS4T083ygngAFA/lp2YiGUnAMKd02nMcKenSytWSC1Om8KprZWysqTiYqmkhCUoAOALM98AAL8VFEh79kizZnkGb8l4PHOmsf80lz4HAN8I3wAAv7kvad6nj+96dzmXPgcA3wjfAAC/uS9pXlzsu95dzqXPAcA31nybiDXfQOSJtu32WPN9SrT92wNoGsx8A0AD2e1SSorndnspKZG924fVKj39tLR6tRG0T9/tJCvLKJ83L/JDaDT+2wNoGsx8m4iZbyByuLfbO3MEtViM+0i/BLjdLk2fbpx86ZaaagTvSH7fEv/2ABqH8G0iwjcQGZxOY5Zz3z7f9RaLlJRk7PoRyTPA0bjsgn97AI3VMtgdAIBwU1BQf/iSjBnRsjKj3fDhpnXLdFZrZL8/X/i3B9BYrPkGgAD5u40e2+1FHv7tATQW4RsAAuTvNnpstxd5+LcH0Fis+TYRa76ByOBe91te7n3SncS630jGvz2AxmLmGwACZLVK8+cb/+3e4cLN/Tgvj/AVifi3B9BYhG8AaIDsbGNLucREz/KkJLaai3T82wNoDJadmIhlJ0Dkicbt9mDg3x5AQxC+TUT4BgAAiG7s8w0AYYhZ11P4WQAIJ4RvAAgzdrs0darnxV6SkowTAaNtvTE/CwDhhmUnJmLZCYDGstulnBzvbe7cO21E0wl//CwAhKOQ3+2kurpajzzyiGJiYjRhwoQGH6e2tlYLFixQv379FBcXp86dO+tXv/qVtm3bdtbnffTRRxo9erQ6deqkNm3a6IorrtCiRYvE/7MAMJvTaczy+hp+3GUPPGC0i3T8LACEq5AO3xs2bFDfvn31/PPP6+TJkw0+Tm1trW677TZNnz5d9957rw4ePKiPPvpIR48e1dVXX613333X5/PsdruGDBmili1b6vPPP9f+/ft1++23695779XEiRMb3B8AaIiCAs/lFWdyuaSyMqNdpONnASBchWz4Xrp0qX71q1/pgQce0NNPP92oY+Xn52v58uV64IEHNHnyZLVp00ZpaWn65z//qTZt2uiuu+7S0aNHPZ5z4MABTZgwQRdeeKGWLl2qiy++WO3atdOMGTM0btw4vfzyy3rjjTca1S8ACERFRdO2C2f8LACEq5AN36mpqfriiy903333yXLmZcQCNG/ePEnSb3/7W4/yDh066LbbbtOBAwf08ssve9S98MILOnr0qO666y7FxsZ61N17772SpL/85S+N6hcABCIhoWnbhTN+FgDCVUDhe8OGDTpy5Ehz9cXDoEGDdOGFFzb6OMXFxdq9e7e6dOmiHj16eNUPGTJEkrRy5UqPcvfjwYMHez1n4MCBat26tT7//HOVl5c3uo8A4I+MDGMnj/rmIywWKTnZaGcGp1PauFFautS4N3N9daj9LADAXwGF7+uuu06dO3fWxRdfrJtvvlmzZ8/WW2+9pdLSUo92d999t5577rkm7WhDbd++XZIxk+5Lt27dPNpJktPp1M6dO+t9XsuWLZWcnOz1PABoTlarsYWe5B063Y/z8szZ49pul1JSpMxMaexY4z4lxSg3Qyj9LAAgEAGF74kTJ2rAgAH6/vvvtWrVKj3xxBPKyclRWlqaOnTooGHDhumee+7R22+/rccff7y5+hyQ/fv3S5I6duzos75Dhw6SpIMHD6q2tlaS9OOPP6q6utqv57mP70tNTY0qKys9bgDQGNnZxhZ6iYme5UlJ5m2t597i78wTHsvLjXKzAngo/CwAIFABXWTnb3/7myTp559/1hdffKFt27Zp69at2rp1qz766CMVFBTo//7v/+RyuXTRRRc1S4cDdfz4cUlSTEyMz/pWrVp5tG3btm3dc/x53ultz5Sbm6vHHnss4D4DwNlkZ0s33xycqzqea4s/i8XY4u/mm83pTzB/FgDQEA26wmXLli3Vt29f9e3bV+PGjZMk/fTTT3r55Zf18MMPq0+fPnrttdeatKMNFRcXJ0k6ceKEz/qamhqvtu57f553etszzZw5U9OmTat7XFlZWbdcBQAaw2qVhg83/3UD2eLPrP4F62cBAA3RZJeXb9u2re6//35169ZNY8aM0eeff660tLSmOnyDde3aVZJ0+PBhn/XuE0jPP/98tWhhrMJp3769Wrdurerqah0+fFgXXHBBvc9LOMup9K1atfKYWQeAcMcWfwDQOE2+1eCoUaPUq1cv5ebmNvWhG6R///6S5HVSqNvu3bslSX379q0rs1qt6tOnT73P+/nnn1VWVub1PACIdGzxBwCN0yz7fKempuqrr75qjkMHrHfv3urWrZsOHjyob775xqt+06ZNkqSbbrrJo3zMmDGSpMLCQq/nfPbZZ6qurtbll1+uxDPP9AGACMYWfwDQOAGF76lTp2rx4sXatm3bWS/3/s0339Qt9zDLjh07NHjwYM2aNcur7qGHHpIkLVq0yKP8xx9/1Jtvvqnzzz9fEyZM8Ki799571a5dO/3973+v2/nELT8/X5I0Y8aMJnwHABD62OIPABonoDXf//M//1N3tcnzzjtPvXr10oABA9S/f3+lp6crLi5Or732mkpKSvTCCy80S4frk5+fr8LCQhUWFmr69Onq1KlTXd1vf/tbvffee/rrX/+q7t27684779T+/fs1efJkHT16VKtXr1a7du08jnf++edr8eLFuv322zV27Fjl5eWpQ4cOys/P18svv6zx48fr9ttvN/U9AkAocG/xN3Wq58mXSUlG8GaLPwCon8Xl8rVhlG+rV6+u21pw69at+u6774yDnDH90a1bN40fP16XX365rrjiCp8nLPrVubNcVn7x4sUes9Xr1q3TbbfdpmHDhslut3s91+l0auHChfrf//1flZSUKC4uThkZGZo9e7YGDBhQ7+t89NFHeuKJJ1RYWKiamhr16tVLkyZN0m9/+9uAL3tfWVkpm80mh8Oh+Pj4gJ4LAKHG6WSLPwAIVEDh+0xHjhzxCOOff/65vv7667qL1bjDadeuXXXFFVd4XcI92hC+AQAAolujwrcvx48fV1FRUV0Y37p1q4qLi3Xy5Ek5nc6mfKmwQ/gGAACIbk0evn1xXxEz2rflI3wDAABEt2bZavBM7itiAgAAANHMlPANAAAAgPANAAAAmIbwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmKRlsDsAAAhPTqdUUCBVVEgJCVJGhmS1BrtXABDaCN8AgIDZ7dLUqdK+fafKkpKk+fOl7Ozg9QsAQh3LTgAAAbHbpZwcz+AtSeXlRrndHpx+AUA4sLhcLlewOxEtKisrZbPZ5HA4FB8fH+zuAEDAnE4pJcU7eLtZLMYMeGkpS1AAwBdmvgEAfisoqD94S5LLJZWVGe0AAN4I3wAAv1VUNG07AIg2hG8AgN8SEpq2HQBEG8I3AMBvGRnGmm6LxXe9xSIlJxvtAADeCN8AAL9ZrcZ2gpJ3AHc/zsvjZEsAqA/hGwAQkOxsadkyKTHRszwpyShnn28AqB9bDZqIrQYBRBKucAkAgeMKlwCABrFapeHDg90LAAgvLDsBAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABM0jLYHQAAhCenUyookCoqpIQEKSNDslqD3SsACG2EbwBAwOx2aepUad++U2VJSdL8+VJ2dvD6BQChjmUnAICA2O1STo5n8Jak8nKj3G4PTr8AIByEdPiurq7W448/rp49e6p169ZKSEjQXXfdpdLSUr+PsWfPHlkslnPeUlJSPJ63cePGgNoDQDRwOo0Zb5fLu85d9sADRjsAgLeQXXZSVVWlkSNHqqioSC+99JJ++ctfqqSkROPHj9eAAQP0wQcfqF+/fn4fr1u3bjrvvPN81pWUlOjSSy/1Km/ZsqW6d+/u8zmJiYl+vzYARIqCAu8Z79O5XFJZmdFu+HDTugUAYSNkw/ecOXNUWFiohQsX6tZbb5Uk9evXT8uXL1fPnj11xx13aPv27bL6eXbP+++/73O2etu2bRowYIDuv/9+r7rExETt2rWrUe8DACJJRUXTtgOAaBOSy06OHz+uZ599VjExMRo/frxHXffu3TVy5Ejt3LlTa9euPeex2rRpo9tvv11t2rTxWb9gwQJ1795dN954Y5P0HQAiWUJC07YDgGgTkuF7/fr1OnbsmPr166e2bdt61Q8ZMkSStHLlynMeq0uXLvrHP/6hLl26eNUdOnRIS5cu1ZQpU9SiRUj+KAAgpGRkGLuaWCy+6y0WKTnZaAcA8BaSiXP79u2SpNTUVJ/13bp182jXUC+++KKsVqsmTpzos76qqkozZszQJZdcotjYWLVv315DhgzRokWLVFtbe87j19TUqLKy0uMGAOHMajW2E5S8A7j7cV4e+30DQH1CMnzv379fktSxY0ef9R06dPBo1xBOp1PPPfecxo0bJ5vN5rPNgQMH5HA4tGrVKjkcDm3dulW9e/fWpEmTNGbMGJ08efKsr5GbmyubzVZ3S05ObnB/ASBUZGdLy5ZJZ553npRklLPPNwDULyTD9/HjxyVJMTExPutbtWrl0a4hVqxYobKyMp8nWkpSenq6Nm7cqPz8fPXs2VMxMTFKTU1Vfn6+MjMztXbtWs2dO/esrzFz5kw5HI66W1lZWYP7CwChJDtb2rNH2rBBWrLEuC8tJXgDwLmEZPiOi4uTJJ04ccJnfU1NjUe7hliwYIGuu+46XXbZZT7rO3XqpGHDhvmsmzJliiTppZdeOutrtGrVSvHx8R43AIgUVquxneCvf23cs9QEAM4tJLca7Nq1qyTp8OHDPuuPHDkiSUpo4On0RUVF+vDDD7VixYoGPT8tLU2S9O233+rEiRP1ztADAAAApwvJme/+/ftLUr1Xsty9e7ckqW/fvg06/oIFC5SSkqIxY8Y06Pmnn2xpqe+UfwAAAOAMIRm+MzMz1aZNGxUVFenYsWNe9Zs2bZIk3XTTTQEf+/Dhw1qyZInuu+++s24veOWVV+qFF17wWVdSUiLJ2HO8vqtmAgAAAGcKyfAdFxen+++/XzU1NXrllVc86nbv3q1169apd+/eGjVqVF35jh07NHjwYM2aNeusx37xxRclSXffffdZ2x08eLDeZSkLFy6UJE2YMOEc7wQAAAA4JSTDtyTNnj1b11xzjR5++GEtX75cVVVVKioq0i233KK4uDi9/vrrHpeWz8/PV2FhoXJzc3Xo0CGfx3RvL3jHHXfUu42hm8Vi0TvvvKNJkybpq6++Uk1NjXbv3q177rlHBQUFuuGGGzRjxowmfc8AAACIbCEbvmNjY7V+/XpNnz5dM2fOVIcOHXT99derd+/e2rp1q/r16+fRPisrSx06dFBWVla9wXrlypXau3evfve7353z9T/88EPNmzdPu3bt0uDBg9WmTRtdccUV+vrrr7Vo0SKtXbuWEy0BAAAQEIvL5XIFuxPRorKyUjabTQ6Hg20HAQAAolDIznwDAAAAkYbwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmKRlsDsAAAhPTqdUUCBVVEgJCVJGhnTahYcBAD4QvgEAAbPbpalTpX37TpUlJUnz50vZ2cHrFwCEOpadAAACYrdLOTmewVuSysuNcrs9OP0CgHDA5eVNxOXlAYQ7p1NKSfEO3m4WizEDXlrKEhQA8IWZbwCA3woK6g/ekuRySWVlRjsAgDfCNwDAbxUVTdsOAKIN4RsA4LeEhKZtBwDRhvANAPBbRoaxptti8V1vsUjJyUY7AIA3wjcAwG9Wq7GdoOQdwN2P8/I42RIA6kP4BgAEJDtbWrZMSkz0LE9KMsrZ5xsA6sdWgyZiq0EAkYQrXAJA4LjCJQCgQaxWafjwYPcCAMILy04AAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJMQvgEAAACTEL4BAAAAkxC+AQAAAJO0DHYHAADhyemUCgqkigopIUHKyJCs1mD3CgBCG+EbABAwu12aOlXat+9UWVKSNH++lJ0dvH4BQKhj2QkAICB2u5ST4xm8Jam83Ci324PTLwAIBxaXy+UKdieiRWVlpWw2mxwOh+Lj44PdHQAImNMppaR4B283i8WYAS8tZQkKAPjCzDcAwG8FBfUHb0lyuaSyMqMdAMAb4RsA4LeKiqZtBwDRhvANAPBbQkLTtgOAaEP4BgD4LSPDWNNtsfiut1ik5GSjHQDAG+EbAOA3q9XYTlDyDuDux3l5nGwJAPUhfAMAApKdLS1bJiUmepYnJRnl7PMNAPUL6fBdXV2txx9/XD179lTr1q2VkJCgu+66S6WlpQEdZ86cObJYLPXeJkyYUO9z//Wvf2nEiBGy2WyKj4/XkCFDtGzZska+MwAIb9nZ0p490oYN0pIlxn1pKcEbAM4lZMN3VVWVrrvuOs2dO1dPPvmkDh8+rH//+9/asWOHBgwYoO3btwd0vPj4ePXs2dPnLaGeM4Py8vI0atQo9ejRQ19//bVKS0t11VVX6dZbb9WcOXOa4F0CQPiyWqXhw6Vf/9q4Z6kJAJxbyF5k57//+7/1l7/8RQsXLtSUKVPqyr/99lv17NlTvXr10vbt22X1Y7SfM2eO9uzZo5dfftnv1y8uLlb//v11+eWX6+OPP5bltMWNmZmZ2rhxozZv3qxBgwb5fUwusgMAABDdQnLm+/jx43r22WcVExOj8ePHe9R1795dI0eO1M6dO7V27dpm68Mzzzwjp9Ope+65xyN4S9KkSZMkSXPnzm221wcAAEDkCcnwvX79eh07dkz9+vVT27ZtveqHDBkiSVq5cmWz9WHVqlWSpMGDB9f7+mvXrtXPP//cbH0AAABAZAnJ8O1ez52amuqzvlu3bh7t/FFaWqqxY8cqOTlZMTExuuCCC3TTTTfp/fff92pbXl6uH374od4+JCYmKiYmRtXV1fr666/97gMAAACiW0iG7/3790uSOnbs6LO+Q4cOHu388cknn2jgwIH65JNP9NNPP2nNmjU6cuSIRo4cqSeffNLn67du3VqxsbFex7JYLGrfvv05+1BTU6PKykqPGwAAAKJXSIbv48ePS5JiYmJ81rdq1cqj3bmMHTtWW7Zs0bRp05SQkKCYmBgNHDhQq1evVocOHfTHP/5RmzZt8vv1/e1Dbm6ubDZb3S05Odmv/gIAACAyhWT4jouLkySdOHHCZ31NTY1Hu3O55JJL1KdPH69ym82mO+64Q5K0ePFiv1/f3z7MnDlTDoej7lZWVuZXfwEAABCZWga7A7507dpVknT48GGf9UeOHJGkevfnDkRaWpok6csvv/R6/erqalVVVXktPXG5XPrxxx/P2YdWrVrVzZADAAAAITnz3b9/f0mq90qWu3fvliT17du30a9VW1srSR7bCSYmJqpz58719qG8vFwnTpxQbGysevTo0eg+AAAAIDqEZPjOzMxUmzZtVFRUpGPHjnnVu9dn33TTTec8VllZmbp27aqioiKf9SUlJZKkXr16eZSPGTNGklRYWFjv6994441q2TIk/3gAAACAEBSS4TsuLk7333+/ampq9Morr3jU7d69W+vWrVPv3r01atSouvIdO3Zo8ODBmjVrlkd7p9Op77//3ucFeY4cOaIlS5ZIkiZMmOBRN336dFmtVv3tb3/TmRcBXbRokSTpoYceavB7BAAAQPQJyfAtSbNnz9Y111yjhx9+WMuXL1dVVZWKiop0yy23KC4uTq+//rrHpeXz8/NVWFio3NxcHTp0qK7cvZzkscceU15envbt26cTJ05oy5YtGj16tBwOh+bMmVN34Ry33r17a968efr4449133336cCBAzp8+LAefPBBrV+/XrNnzw7o0vIAAACAxXXmtG4Iqa6u1p///Ge9/vrr2rt3r2w2m66//no9/vjjdRfacVu3bp1uu+02DRs2THa73WMN92effaY333xT7777rr755htVV1erc+fOGjRokKZMmaKRI0fW24e1a9dq7ty5+uyzz+RyudS3b1/9/ve/V05OTsDvp7KyUjabTQ6HQ/Hx8QE/HwAAAOEtpMN3pCF8AwAARDfOFgQANIjTKRUUSBUVUkKClJEhnbYaEADgA+EbABAwu12aOlXat+9UWVKSNH++lJ0dvH4BQKgL2RMuAQChyW6XcnI8g7cklZcb5XZ7cPoFAOGANd8mYs03gHDndEopKd7B281iMWbAS0tZggIAvjDzDQDwW0FB/cFbklwuqazMaAcA8Eb4BgD4raKiadsBQLQhfAMA/JaQ0LTtACDaEL4BAH7LyDDWdJ92HTMPFouUnGy0AwB4I3wDAPxmtRrbCUreAdz9OC+Pky0BoD6EbwBAQLKzpWXLpMREz/KkJKOcfb4BoH5sNWgithoEEEm4wiUABI4rXAIAGsRqlYYPD3YvACC8sOwEAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwCeEbAAAAMAnhGwAAADAJ4RsAAAAwSctgdwAAEJ6cTqmgQKqokBISpIwMyWoNdq8AILQRvgEAAbPbpalTpX37TpUlJUnz50vZ2cHrFwCEOpadAAACYrdLOTmewVuSysuNcrs9OP0CgHBgcblcrmB3IlpUVlbKZrPJ4XAoPj4+2N0BgIA5nVJKinfwdrNYjBnw0lKWoACAL8x8AwD8VlBQf/CWJJdLKisz2gEAvBG+AQB+q6ho2nYAEG0I3wAAvyUkNG07AIg2hG8AgN8yMow13RaL73qLRUpONtoBALwRvgEAfrNaje0EJe8A7n6cl8fJlgBQH8I3ACAg2dnSsmVSYqJneVKSUc4+3wBQP7YaNBFbDQKIJFzhEgACxxUuAQANYrVKw4cHuxcAEF5YdgIAAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmCSkw3d1dbUef/xx9ezZU61bt1ZCQoLuuusulZaWBnScLVu26He/+5369OmjuLg4tW7dWmlpaZo8eXK9x9q4caMsFku9t5SUlCZ4hwAAAIgmIbvVYFVVlUaOHKmioiK99NJL+uUvf6mSkhKNHz9eAwYM0AcffKB+/fqd8zhr1qzR6NGjdfHFF2v+/PkaOnSoTp48qZUrV2ratGl67bXX9M477+jaa6/1em7Lli3VvXt3n8dNPPPqEgAAAMA5hGz4njNnjgoLC7Vw4ULdeuutkqR+/fpp+fLl6tmzp+644w5t375d1nNc0aGqqkqS9MYbb+jqq6+uK//Nb36jn3/+WZMnT9bdd9+tL7/80uu5iYmJ2rVrVxO+KwAAAESzkFx2cvz4cT377LOKiYnR+PHjPeq6d++ukSNHaufOnVq7du05j9WxY0ddf/31HsHb7aabbpIk7dq1S3v37m2azgMAAAD1CMnwvX79eh07dkz9+vVT27ZtveqHDBkiSVq5cuU5jzVixAi9++67PutsNlvjOgoAAAAEICTD9/bt2yVJqampPuu7devm0a6h3EtKunXrposuusirvqqqSjNmzNAll1yi2NhYtW/fXkOGDNGiRYtUW1vbqNcGAABA9AnJ8L1//35JxpIRXzp06ODRrqFeffVVSdKjjz7qs/7AgQNyOBxatWqVHA6Htm7dqt69e2vSpEkaM2aMTp48edbj19TUqLKy0uMGAACA6BWS4fv48eOSpJiYGJ/1rVq18mjXEEVFRXr++ed1yy23eK0rl6T09HRt3LhR+fn56tmzp2JiYpSamqr8/HxlZmZq7dq1mjt37llfIzc3Vzabre6WnJzc4P4CAAAg/IVk+I6Li5MknThxwmd9TU2NR7tAHTx4ULfeequuuuqqutnvM3Xq1EnDhg3zWTdlyhRJ0ksvvXTW15k5c6YcDkfdraysrEH9BQAAQGQIya0Gu3btKkk6fPiwz/ojR45IkhISEgI+9pEjR/SLX/xCXbp00Zo1axoU4NPS0iRJ3377rU6cOHHWGXr3LD0ARBqnUyookCoqpIQEKSNDOsfurwAQ9UJy5rt///6SVO/VJ3fv3i1J6tu3b0DHPXjwoDIzM9WhQwe9++67Dd7t5PSTLS0WS4OOAQDhzG6XUlKkzExp7FjjPiXFKAcA1C8kw3dmZqbatGmjoqIiHTt2zKt+06ZNkk7t0+2P//znPxo2bJgSExO1evVqjxnvF198UTt37vRof+WVV+qFF17weaySkhJJxp7j5513nt99AIBIYLdLOTnSvn2e5eXlRjkBHADqF5LhOy4uTvfff79qamr0yiuveNTt3r1b69atU+/evTVq1Ki68h07dmjw4MGaNWuW1/H27t2roUOH6tJLL9Vbb72l1q1be9T/6U9/0pYtWzzKDh48qBUrVvjs38KFCyVJEyZMaMC7A4Dw5XRKU6dKLpd3nbvsgQeMdgAAbyG55luSZs+erQ8//FAPP/ywLrjgAo0aNUolJSUaP3684uLi9Prrr3tcWj4/P1+FhYUqLCzU9OnT1alTJ0lGWM/MzFR5ebn69++vcePGeb3WwYMHvcosFoveeecdTZo0SdOmTVNKSorKy8uVm5urgoIC3XDDDZoxY0bz/QAAIAQVFHjPeJ/O5ZLKyox2w4eb1i0ACBshG75jY2O1fv16/fnPf9bMmTN1xx13yGaz6frrr9fy5cvrLrTjlpWVpSVLlmjYsGEe+4OvXLmy7tLxy5cv9/v1P/zwQ7355ptauXKlBg8erB9//FHt2rVT3759tWjRIt19991q0SIk/3AAAM2moqJp2wFAtLG4XL7+eIjmUFlZKZvNJofDofj4+GB3BwACtnGjcXLluWzYwMw3APjC1C0AwG8ZGVJSklTfRk8Wi5ScbLQDAHgjfAMA/Ga1SvPnG/99ZgB3P87LY79vAKgP4RsAEJDsbGnZMikx0bM8Kckoz84OTr8AIByw5ttErPkGEEm4wiUABC5kdzsBAIQ2q5WTKgEgUCw7AQAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABM0jLYHQAAhCenUyookCoqpIQEKSNDslqD3SsACG2EbwBAwOx2aepUad++U2VJSdL8+VJ2dvD6BQChjmUnAICA2O1STo5n8Jak8nKj3G4PTr8AIBxYXC6XK9idiBaVlZWy2WxyOByKj48PdncAIGBOp5SS4h283SwWYwa8tJQlKADgCzPfAAC/FRTUH7wlyeWSysqMdgAAb4RvAIDfKiqath0ARBvCNwDAbwkJTdsOAKIN4RsA4LeMDGNNt8Xiu95ikZKTjXYAAG+EbwCA36xWYztByTuAux/n5XGyJQDUh/ANAAhIdra0bJmUmOhZnpRklLPPNwDUj60GTcRWgwAiCVe4BIDAcYVLAECDWK3S8OHB7gUAhBeWnQAAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmIXwDAAAAJiF8AwAAACYhfAMAAAAmaRnsDkQTl8slSaqsrAxyTwAAAHA27dq1k8ViafLjEr5NdOjQIUlScnJykHsCAACAszlw4IC6dOnS5MclfJuoY8eOkqS9e/fKZrMFuTeIZJWVlUpOTlZZWZni4+OD3R1EMD5rMAufNZjF/VmLiYlpluMTvk3UooWxxN5mszFwwBTx8fF81mAKPmswC581mKU5lpxInHAJAAAAmIbwDQAAAJiE8G2iVq1aafbs2WrVqlWwu4IIx2cNZuGzBrPwWYNZmvuzZnG5978DAAAA0KyY+QYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGIkhVVZUeffRRXXLJJWrdurUuvPBCTZw4UeXl5cHuGgAAEOHbFJ9++qnGjRuntLQ0WSwW/eEPfwh2lxCBqqurNWLECD3xxBP66aefdPPNNys5OVmLFy/WgAEDtHv37mB3EWGEcQtApHnllVc0cOBAtW/fXm3atNHll1+uf/zjH6b3g8vLm2DTpk366KOPNGTIEP3www/B7g4i1JNPPqmPPvpI11xzjd599121bdtWkvTMM89o+vTpmjhxojZu3BjcTiJsMG4BiDRHjhxRVlaW+vfvr9atW2vFihX69a9/rdatWysrK8u0frDPtwlqa2vVooXxR4aUlBTdeeedevLJJ4PcK0SSEydO6Pzzz5fD4dDnn3+uAQMGeNT369dPRUVF+vTTT3XFFVcEqZcIJ4xbCLZPP/1UV155pV9tR44cqffee6+Ze4RINGTIECUkJOif//ynaa/JzLcJ3F9gQHPZtGmTHA6Hunfv7hW8JSknJ0dFRUVatWoV4Rt+YdxCsFVVVenOO+/Ua6+9prZt2+qWW27xqH/jjTdUXV2tcePG6cYbbwxSLxHuOnXqpJMnT5r6moRvIAJs375dknT55Zf7rHeXFxUVmdYnAGiMjIwMde7cWa+99poGDhyol19+ua7up59+0t///nclJSXplVdeCV4nEZZ+/vlnHT9+XP/617/03nvvafny5aa+PuEbiAB79+6VJCUlJfmsd5d/9913pvUJABprx44dkqT09HSP8qKiItXW1qpv377B6BbC2P79+5WQkCBJslqteu6550z/ywnhG4gAP/30kyQpLi7OZ32bNm0kSUePHjWtTwDQWO6/1p0Zsrdt2+azHDiXzp07a8uWLTp69Kj+/e9/6/7771enTp28ljU1J8I3AAAISe7wfebMN+EbDdWyZUsNHDhQkpSZmanDhw9r5syZpoZvzqgBIoB7W8Hjx4/7rD927JgkqV27dqb1CQAaq6ioSBaLRX369PEod4fvfv36BaFXiCT9+/c3/ToYhO8Q9Omnn8pisfh1u/7664PdXYSAiy66SJK0b98+n/Xu8osvvti0PiF6MGahOVRWVuq7775Tt27d6pbOSZLT6VRxcbFatWqlnj17BrGHiASFhYVKSUnxKm/OcY1lJyY4ePCgPvjgA0nGzOSuXbu0bNkytWnTxucif7ZXQqDcsz+ff/65z3p3OX+ihb8CGbcYs9Ac6lvvXVJSoqqqKg0YMEBWqzUYXUOYyszM1C233KJevXqpurpab7/9tpYsWaJFixZ5tW3Wcc2FZrdhwwaXJK/bxRdfXO9zvvjiC5ck1/Dhwz3Kjx496mrRooUrKSmpmXuNcFJTU+Oy2WwuSa6tW7d61fft29clyfXpp5+a3zmEpUDHLcYsNLVnn33WJcn16KOPepS/9dZbLkmu//qv/wpSzxCupk6d6rrkkktcsbGxrs6dO7uGDh3qWrVqVb3tm2tcY9mJCYYPHy6Xy+V127NnT73PYXslBCImJkb333+/JGnKlCl1a7wl4/LyRUVFGjZsGBfYgd8CHbcYs9DU6pv5du/uVFRUpDVr1qiqqsr0viE85eXl6auvvtLx48fr/ro3evToets317jGshM//OpXv9KXX34Z0HNeffVVXXXVVQ1+TbZXQqD+8Ic/aN26dSosLFSPHj2UkZGh7777Th9//LG6dOmil156KdhdRARjzEJTqy/4jBgxQv3799fOnTs1bdo0ffHFF8HoHqJAc41rhG8/lJaW6quvvgroOfXtOuEvtldCoFq3bq0NGzYoNzdXS5Ys0YoVK9SxY0dNmDBBTzzxRL0X4AGaAmMWmtqmTZt8ll944YXaunWryb1BNGqucc3icrlcjeoZmkVKSor27t2ro0ePepzlfdVVV2nLli3auXOnLrvssiD2EABOYcwCEGmaa1wjfDczi8Xid1v3P0VlZaVsNpu6d++ub775pq7e6XSqXbt2qq2t1bFjxzjLG0CzCHTcYswCEOpCaVzjhMtm5nK59MYbb2jUqFHq2rWr4uPjlZGRoYKCAq8TmdzOtb3SZZddxpcYgGYT6LjFmAUg1IXSuEb4NkFeXp46d+6sZ599Vv/85z+VmJio6667Ttu3b/fZvr41Rrt27ZIkLioAoNkFMm4xZgEIB6EyrnHCpQlWrVqlTp061T0eOXKk0tPT9eyzz/rc2N3f7ZVGjBih2NjYZuw5gGgVyLjFmAUgHITKuMbMtwlO/4eWpBYtWqhPnz4qLS312f5c2yuVlJRo2rRpiomJaZ4OA4h6gYxbjFkAwkGojGuccBkETqdTPXr00A033KDnn38+2N0BgHNi3AIQaYI1rjHzHQQLFy7U3r17dd999wW7KwDgF8YtAJEmWOMa4dtkH3/8sR5++GH94Q9/8PpTBgCEIsYtAJEmmOMay05MtGfPHg0aNEhDhw7VG2+8EdCekwAQDIxbACJNsMc1wrdJfvzxRw0ePFjt2rXThg0bOOMfQMhj3AIQaUJhXGOrQROcOHFC2dnZOn78uNavX88XGICQx7gFINKEyrhG+DbBfffdpw8++EAvvviiSktL67a0adWqlQYMGBDk3gGAN8YtAJEmVMY1lp2YICUlRd99951X+cUXX6w9e/aY3yEAOAfGLQCRJlTGNcI3AAAAYBK2GgQAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAAAExC+AYAAABMQvgGAAAATEL4BgAErKKiQlOmTNGkSZN04sQJTZ06VRdccIE6d+6sP/7xj6qtrQ12FwEgJBG+AQABee+999S7d28999xzuuyyy/TYY49pwYIF6tGjh2pqavTkk0/qySefDHY3ASAkWVwulyvYnQAAhIdvvvlG/fv3V3V1tV5//XXl5OToggsuUKtWrbR3716tWLFCOTk5stlsOnz4sFq0YI4HAE7HqAgA8Nt///d/69ixY3rwwQd1++23a9u2bTp06JCGDh0qq9Wq7OxstW/fXg6HQ2VlZcHuLgCEHMI3AMAv33//vd5++221bdtWs2fPliR9/vnnkqT+/ftLkiwWi7p06SJJcjqdQeknAIQywjcAwC8bNmyQ0+nUDTfcoNjYWEnStm3bJEn9+vWra3fo0CFJUkJCgul9BIBQR/gGAPhl+/btkqSrr766rswdvt0z32VlZTp8+LB69epVF9ABAKcQvgEAfvnhhx8kSe3bt5ckuVwuFRUV6YILLlDXrl0lSWvWrJEkjRo1Kih9BIBQR/gGAPilXbt2kqR9+/ZJkkpKSvTTTz/VLTmprq7WM888o/POO0/33ntv0PoJAKGM8A0A8MuQIUMkSa+++qocDofHkhOHw6Fbb71VJSUlmjVrlnr06BHEngJA6GKfbwCAX5xOp6655hpt2bJFycnJ6tixo7Zv364BAwaorKxMP/zwgyZPnqxnn31WFosl2N0FgJBE+AYA+O3QoUN66KGH9Pbbb+vw4cOSpDZt2uiKK67Q73//e2VlZQW3gwAQ4gjfAIAG6dy5s06ePKkff/yRmW4A8BNrvgEAAfvuu+906NAhXX755QRvAAgA4RsAEDD3lS2vuOKKIPcEAMIL4RsAEDDCNwA0DGu+AQAAAJMw8w0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACYhPANAAAAmITwDQAAAJiE8A0AAACY5P8B+anrerSvMjwAAAAASUVORK5CYII=",
      "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(r'$\\varphi$')\n",
    "plt.ylabel(r'$\\mu$')\n",
    "tick_pos= [ -np.pi/2, 0, np.pi/2, np.pi, 3*np.pi/2]\n",
    "labels = [r'$-\\frac{1}{2}\\pi$', r'$0$', r'$\\frac{1}{2}\\pi$', r'$\\pi$', r'$\\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.12.7"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}