{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Föppl-von Kármán equations are given by: \n",
    "\\begin{align*}\n",
    "D \\nabla^4 w &= q + \\frac{\\partial^2 F}{\\partial y^2}\\frac{\\partial^2 w}{\\partial x^2} - 2 \\frac{\\partial^2 F}{\\partial x \\partial y} \\frac{\\partial^2 w}{\\partial x \\partial y} + \\frac{\\partial^2 F}{\\partial x^2}\\frac{\\partial^2 w}{\\partial y^2} \\\\\n",
    "\\nabla^4 F &= Eh \\left[ \\left( \\frac{\\partial^2 w}{\\partial x \\partial y}\\right)^2 - \\frac{\\partial^2 w}{\\partial x^2} \\frac{\\partial^2 w}{\\partial y^2} \\right]\n",
    "\\end{align*}\n",
    "\n",
    "Define the \"diamond\" operator as \n",
    "\\begin{gather*}\n",
    "\\diamondsuit (f,g) = \\frac{1}{2} \\left\\{ (\\nabla^2 f) (\\nabla^2 g) + \\nabla^2 \\left(f \\nabla^2 g + g \\nabla^2 f \\right) \\right\\} - \\frac{1}{4} \\left\\{ \\nabla^4(fg) + f \\nabla^4 g + g \\nabla^4 f \\right\\}\n",
    "\\end{gather*}\n",
    "\n",
    "(I learnt about the diamond operator from a 1962 paper: [E. H. Mansfield, \"Bending, buckling and curling of a heated thin plate\", Proc. R. Soc. Lond. A268, 316–327 (1962)](https://doi.org/10.1098/rspa.1962.0143). )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "x, y = sp.symbols('x, y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def laplacian(f):\n",
    "    return sp.diff(f,x,2)+sp.diff(f,y,2)\n",
    "\n",
    "def biharmonic(f):\n",
    "    return laplacian(laplacian(f))\n",
    "\n",
    "def diamond(f,g):\n",
    "    term1 = laplacian(f)*laplacian(g) + laplacian(f*laplacian(g)+g*laplacian(f))\n",
    "    term2 = biharmonic(f*g) + f*biharmonic(g) + g*biharmonic(f)\n",
    "    return (sp.Rational(1,2)*term1 - sp.Rational(1,4)*term2).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\partial^{2}}{\\partial x^{2}} f{\\left(x,y \\right)} + \\frac{\\partial^{2}}{\\partial y^{2}} f{\\left(x,y \\right)}$"
      ],
      "text/plain": [
       "Derivative(f(x, y), (x, 2)) + Derivative(f(x, y), (y, 2))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = sp.Function('f')(x,y)\n",
    "g = sp.Function('g')(x,y)\n",
    "laplacian(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\partial^{4}}{\\partial x^{4}} f{\\left(x,y \\right)} + \\frac{\\partial^{4}}{\\partial y^{4}} f{\\left(x,y \\right)} + 2 \\frac{\\partial^{4}}{\\partial y^{2}\\partial x^{2}} f{\\left(x,y \\right)}$"
      ],
      "text/plain": [
       "Derivative(f(x, y), (x, 4)) + Derivative(f(x, y), (y, 4)) + 2*Derivative(f(x, y), (x, 2), (y, 2))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "biharmonic(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2 \\frac{\\partial^{2}}{\\partial x^{2}} f{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial y^{2}} f{\\left(x,y \\right)} - 2 \\left(\\frac{\\partial^{2}}{\\partial y\\partial x} f{\\left(x,y \\right)}\\right)^{2}$"
      ],
      "text/plain": [
       "2*Derivative(f(x, y), (x, 2))*Derivative(f(x, y), (y, 2)) - 2*Derivative(f(x, y), x, y)**2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diamond(f,f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\partial^{2}}{\\partial x^{2}} f{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial y^{2}} g{\\left(x,y \\right)} + \\frac{\\partial^{2}}{\\partial y^{2}} f{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial x^{2}} g{\\left(x,y \\right)} - 2 \\frac{\\partial^{2}}{\\partial y\\partial x} f{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial y\\partial x} g{\\left(x,y \\right)}$"
      ],
      "text/plain": [
       "Derivative(f(x, y), (x, 2))*Derivative(g(x, y), (y, 2)) + Derivative(f(x, y), (y, 2))*Derivative(g(x, y), (x, 2)) - 2*Derivative(f(x, y), x, y)*Derivative(g(x, y), x, y)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diamond(f,g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above two results give us an indication that the RHS of the two Föppl-von Kármán equations can be written in terms of the diamond operator. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = sp.Function('w')(x,y)\n",
    "F = sp.Function('F')(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the expression: $\\displaystyle -\\frac{1}{2} \\diamondsuit(w,w)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle -\\frac{1}{2} \\diamondsuit (w,w) = - \\frac{\\partial^{2}}{\\partial x^{2}} w{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial y^{2}} w{\\left(x,y \\right)} + \\left(\\frac{\\partial^{2}}{\\partial y\\partial x} w{\\left(x,y \\right)}\\right)^{2}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import Math\n",
    "\n",
    "display(Math(r'-\\frac{{1}}{{2}} \\diamondsuit (w,w) = {}'.format(sp.latex(-sp.Rational(1,2)*diamond(w,w))))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, consider the expression: $\\diamondsuit(F,w)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\diamondsuit (F,w) = \\frac{\\partial^{2}}{\\partial x^{2}} F{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial y^{2}} w{\\left(x,y \\right)} + \\frac{\\partial^{2}}{\\partial y^{2}} F{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial x^{2}} w{\\left(x,y \\right)} - 2 \\frac{\\partial^{2}}{\\partial y\\partial x} F{\\left(x,y \\right)} \\frac{\\partial^{2}}{\\partial y\\partial x} w{\\left(x,y \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(r'\\diamondsuit (F,w) = {}'.format(sp.latex(diamond(F,w)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the Föppl-von Kármán equations can be rewritten as:\n",
    "\n",
    "\\begin{align*}\n",
    "D \\nabla^4 w &= q + \\diamondsuit(F,w) \\\\\n",
    "\\nabla^4 F &= -\\frac{1}{2}Eh \\; \\diamondsuit(w,w)\n",
    "\\end{align*}\n",
    "\n",
    "The primary motivation of rewriting the Föppl-von Kármán equations in this form is that it does not depend, representation-wise on any coordinate system. So we can apply this form to the polar coordinate system also whereby to proceed with the analysis we only to substitute the relevant forms for the Laplacians and biharmonics. In fact, this is exactly what we will do to study the buckling of a thin circular plate. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Buckling of a thin circular plate clamped at the periphery and loaded by an in-plane compressive load along the periphery"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The thickness of the plate is uniform ($h$). The origin of the cylindrical coordinate system is set at the centre of the mid-plane; the thickness varies from $z=-h/2$ to $z=h/2$. Because the analysis will be done throughout by carrying out integrations along the thickness, it effectively means we are using a polar coordinate system.\n",
    "\n",
    "The radius of the circle is $R$. The magnitude of the compressive load is $P$ per unit length of the circumference (so that the unit for P is N/m). The problem is axisymmetric.\n",
    "\n",
    "We first import the various definitions in the `polarUtilities.py` file we had created in our previous course on Applied Elasticity (ME60401) (Here is the GitHub repository [link](https://github.com/jeevanjyoti4/elasticity))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from polarUtilities import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $r$ and $\\theta$ symbols are already defined in `polarUtilies.py`, and so we do not have to be define them separately here. \n",
    "\n",
    "We set the symbol $D$ for the bending rigidity ($\\displaystyle \\frac{E h^3}{12(1-\\nu^2)}$), where $E$ is the Young's modulus and $\\nu$ is the Poisson's ratio. \n",
    "\n",
    "The deflection of the plate in the $z$-direction is represented by $\\zeta(r)$, where the dependence is solely on $r$ due to axisymmetry. \n",
    "\n",
    "Just like $F$ represented the (integrated) stress function in the rectangular Cartesian coordinate system, we will use $\\varphi(r)$ to represent the (integrated) stress function in the polar coordinate system. Again note that $\\varphi(r)$ is dependent solely on $r$ due to axisymmetry. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = sp.symbols('D')\n",
    "zeta = sp.Function('zeta')(r)\n",
    "varphi = sp.Function('varphi')(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We next define the diamond operator for the polar coordinate system using the Laplacian and biharmonic operator definitions present in the `polarUtilities.py` file that we have already imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polardiamond(f,g):\n",
    "    term1 = polarLaplacian(f)*polarLaplacian(g) + polarLaplacian(f*polarLaplacian(g)+g*polarLaplacian(f))\n",
    "    term2 = polarbiharmonic(f*g) + f*polarbiharmonic(g) + g*polarbiharmonic(f)\n",
    "    return (sp.Rational(1,2)*term1 - sp.Rational(1,4)*term2).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\frac{d}{d r} \\varphi{\\left(r \\right)} \\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)} + \\frac{d^{2}}{d r^{2}} \\varphi{\\left(r \\right)} \\frac{d}{d r} \\zeta{\\left(r \\right)}}{r}$"
      ],
      "text/plain": [
       "(Derivative(varphi(r), r)*Derivative(zeta(r), (r, 2)) + Derivative(varphi(r), (r, 2))*Derivative(zeta(r), r))/r"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polardiamond(varphi,zeta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{2 \\frac{d}{d r} \\zeta{\\left(r \\right)} \\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)}}{r}$"
      ],
      "text/plain": [
       "2*Derivative(zeta(r), r)*Derivative(zeta(r), (r, 2))/r"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polardiamond(zeta,zeta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the two Föppl-von Kármán equations in polar coordinate system under axisymmetric conditions become:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle D \\left(\\frac{d^{4}}{d r^{4}} \\zeta{\\left(r \\right)} + \\frac{2 \\frac{d^{3}}{d r^{3}} \\zeta{\\left(r \\right)}}{r} - \\frac{\\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)}}{r^{2}} + \\frac{\\frac{d}{d r} \\zeta{\\left(r \\right)}}{r^{3}}\\right) = q + \\frac{\\frac{d}{d r} \\varphi{\\left(r \\right)} \\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)} + \\frac{d^{2}}{d r^{2}} \\varphi{\\left(r \\right)} \\frac{d}{d r} \\zeta{\\left(r \\right)}}{r}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d^{4}}{d r^{4}} \\varphi{\\left(r \\right)} + \\frac{2 \\frac{d^{3}}{d r^{3}} \\varphi{\\left(r \\right)}}{r} - \\frac{\\frac{d^{2}}{d r^{2}} \\varphi{\\left(r \\right)}}{r^{2}} + \\frac{\\frac{d}{d r} \\varphi{\\left(r \\right)}}{r^{3}} = -\\frac{1}{2}Eh \\frac{2 \\frac{d}{d r} \\zeta{\\left(r \\right)} \\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)}}{r}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Math(r'{} = q + {}'.format(sp.latex(D*polarbiharmonic(zeta)),sp.latex(polardiamond(varphi,zeta)))))\n",
    "\n",
    "display(Math(r'{} = -\\frac{{1}}{{2}}Eh {}'.format(sp.latex(polarbiharmonic(varphi)),sp.latex(polardiamond(zeta,zeta)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will be easier to proceed from here if we define $\\displaystyle \\Psi(r):= \\frac{{\\rm d}\\varphi}{{\\rm d}r}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Psi{\\left(r \\right)} := \\frac{d}{d r} \\varphi{\\left(r \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Psi = sp.Function('Psi')(r)\n",
    "display(Math(r'{} := {}'.format(sp.latex(Psi),sp.latex(sp.diff(varphi,r)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To proceed with the solution of the two equations, we first take up the second, and note that in the situation just before buckling starts, $\\zeta = 0$. Then, the second of the Foppl-von Karman equations becomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d^{3}}{d r^{3}} \\Psi{\\left(r \\right)} + \\frac{2 \\frac{d^{2}}{d r^{2}} \\Psi{\\left(r \\right)}}{r} - \\frac{\\frac{d}{d r} \\Psi{\\left(r \\right)}}{r^{2}} + \\frac{\\Psi{\\left(r \\right)}}{r^{3}} = 0$"
      ],
      "text/plain": [
       "Eq(Derivative(Psi(r), (r, 3)) + 2*Derivative(Psi(r), (r, 2))/r - Derivative(Psi(r), r)/r**2 + Psi(r)/r**3, 0)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lhs = polarLaplacian(sp.diff(Psi,r)+Psi/r)\n",
    "eqn = sp.Eq(lhs,0)\n",
    "display(eqn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Psi{\\left(r \\right)} = \\frac{C_{1}}{r} + C_{2} r + C_{3} r \\log{\\left(r \\right)}$"
      ],
      "text/plain": [
       "Eq(Psi(r), C1/r + C2*r + C3*r*log(r))"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.dsolve(eqn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, it is important to note that since $\\displaystyle \\int_{-h/2}^{h/2}\\sigma_{rr}\\; {\\rm d}z := \\frac{1}{r}\\frac{{\\rm d} \\varphi}{{\\rm d} r}$ and $\\displaystyle \\int_{-h/2}^{h/2}\\sigma_{\\theta \\theta}\\; {\\rm d}z := \\frac{{\\rm d^2} \\varphi}{{\\rm d} r^2}$, therefore we have $\\displaystyle \\int_{-h/2}^{h/2}\\sigma_{rr}\\; {\\rm d}z := \\frac{\\Psi}{r}$ and $\\displaystyle \\int_{-h/2}^{h/2}\\sigma_{\\theta \\theta}\\; {\\rm d}z := \\frac{ {\\rm d} \\Psi}{{\\rm d} r}$.\n",
    "\n",
    "So, in order to have finite values of $\\displaystyle \\int_{-h/2}^{h/2}\\sigma_{rr}\\; {\\rm d} z$, it is necessary that in the expression for $\\Psi$, $C_1 = C_2 = 0$. Therefore, the expression for $\\Psi$ becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Psi(r) = C_{2} r$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "C2 = sp.symbols('C_2')\n",
    "Psi = C2*r\n",
    "display(Math(r'\\Psi(r) = {}'.format(sp.latex(Psi))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine $C_2$, we use the boundary condition that at $r=R$, $\\displaystyle \\int_{-h/2}^{h/2}\\sigma_{rr}\\; {\\rm d}z = -P$, where $P$ is the magnitude of the compressive force per unit length of the circular periphery. From this condition, we obtain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C_2 = - P$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Psi(r) = - P r$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "P = sp.symbols('P',positive=True)\n",
    "C2 = -P\n",
    "Psi = C2*r\n",
    "display(Math(r'C_2 = {}'.format(sp.latex(C2))))\n",
    "display(Math(r'\\Psi(r) = {}'.format(sp.latex(Psi))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we go to the first of the Föppl-von Kármán equations, set $q=0$, use $\\displaystyle \\frac{{\\rm d} \\varphi}{{\\rm d}r} = \\Psi(r)$, and use the expression of $\\Psi(r) = -Pr$ to obtain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle D \\left(\\frac{d^{4}}{d r^{4}} \\zeta{\\left(r \\right)} + \\frac{2 \\frac{d^{3}}{d r^{3}} \\zeta{\\left(r \\right)}}{r} - \\frac{\\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)}}{r^{2}} + \\frac{\\frac{d}{d r} \\zeta{\\left(r \\right)}}{r^{3}}\\right) = \\frac{- P r \\frac{d^{2}}{d r^{2}} \\zeta{\\left(r \\right)} - P \\frac{d}{d r} \\zeta{\\left(r \\right)}}{r}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rhs_FvK1 = 1/r*(Psi*sp.diff(zeta,r,2)+ sp.diff(Psi,r)*sp.diff(zeta,r))\n",
    "display(Math(r'{} = {}'.format(sp.latex(D*polarbiharmonic(zeta)),sp.latex(rhs_FvK1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set $\\displaystyle \\frac{P}{D} =: K^2$, and rewrite the equation as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\displaystyle \\frac{1}{r} \\frac{{\\rm d}}{{\\rm d} r} \\left[ r \\frac{{\\rm d}}{{\\rm d} r} \\left\\{ \\frac{1}{r} \\frac{{\\rm d}}{{\\rm d} r} \\left( r \\frac{{\\rm d}\\zeta}{{\\rm d}r} \\right) \\right\\} \\right] = -K^2 \\frac{1}{r} \\frac{{\\rm d}}{{\\rm d}r} \\left( r \\frac{{\\rm d \\zeta}}{{\\rm d}r} \\right) $\n",
    "\n",
    "We cancel the $\\displaystyle \\frac{1}{r}$ from both sides, integrate once, and divide by $r$ to obtain:\n",
    "\n",
    "$\\displaystyle \\frac{{\\rm d}}{{\\rm d} r} \\left\\{ \\frac{1}{r} \\frac{{\\rm d}}{{\\rm d} r} \\left( r \\frac{{\\rm d}\\zeta}{{\\rm d}r} \\right) \\right\\} = -K^2 \\frac{{\\rm d}\\zeta}{{\\rm d}r} + \\frac{A_1}{r}$\n",
    "\n",
    "\n",
    "We next define $\\displaystyle W := \\frac{{\\rm d}\\zeta}{{\\rm d}r}$, which we note is nothing but the slope of the plate as it deforms in the transverse direction. In terms of $W$, the equation becomes\n",
    "\n",
    "$\\displaystyle \\frac{{\\rm d}}{{\\rm d} r} \\left\\{ \\frac{1}{r} \\frac{{\\rm d}}{{\\rm d} r} \\left( r W \\right) \\right\\} = -K^2 W + \\frac{A_1}{r}$\n",
    "\n",
    "We immediately note that it in order to ensure that the solution for the slope comes out finite everywhere, it is necessary that $A_1 = 0$. We use this condition, and rewrite the above equation as\n",
    "\n",
    "$\\displaystyle \\frac{{\\rm d^2} W}{{\\rm d}r^2} + \\frac{1}{r} \\frac{{\\rm d} W}{{\\rm d}r} - \\frac{W}{r^2} = - K^2 W$ \n",
    "\n",
    "or, $\\displaystyle r^2\\frac{{\\rm d^2} W}{{\\rm d}r^2} + r\\frac{{\\rm d} W}{{\\rm d}r} + (K^2r^2 - 1)W = 0$\n",
    "\n",
    "We solve this equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle W{\\left(r \\right)} = C_{1} J_{1}\\left(K r\\right) + C_{2} Y_{1}\\left(K r\\right)$"
      ],
      "text/plain": [
       "Eq(W(r), C1*besselj(1, K*r) + C2*bessely(1, K*r))"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K = sp.symbols('K', positive=True)\n",
    "W = sp.Function('W')(r)\n",
    "\n",
    "sol = sp.dsolve(r**2*sp.diff(W,r,2)+r*sp.diff(W,r)+(K**2*r**2-1)*W)\n",
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above solution, $J_1$ and $Y_1$ are the Bessel functions of first order of the first and second kinds, respectively. We can see this if we extract the different parts making up the solution explicitly as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(W(r), C1*besselj(1, K*r) + C2*bessely(1, K*r))"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol.args"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, `besselj` and `bessely` are the Bessel functions of the first and second kinds, respectively. We express the arguments of $J_1$ and $Y_1$ by normalizing $r$ by the radius $R$ of the circle so that\n",
    "\n",
    "$W(r) = C_1 J_1(\\bar{K} \\bar{r}) + C_2 Y_1 (\\bar{K} \\bar{r})$,\n",
    "\n",
    "where $\\displaystyle \\bar{r} = \\frac{r}{R}$ and $\\displaystyle \\bar{K} = K R = \\sqrt{\\frac{P}{D}}R$\n",
    "\n",
    "Next, we plot the two Bessel functions for $\\bar{r} \\in [0,\\;1]$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$\\\\bar{r}$')"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuMAAAHeCAYAAAAxcQVKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZxkdXX//9eppffp6dlgWGZANkWCiI4b8tW4RUUlfmOMmhg1UdBEjUqMmpjEoCZG0WDcosS4/IwGjbti0Lh9FReQfRUBkWEGZpjpvbu6az2/P27d6uqeXqqmb9W9Xf1+PuzHdN9769YZ9XH7zKnzOR9zd0REREREpP1ScQcgIiIiIrJeKRkXEREREYmJknERERERkZgoGRcRERERiYmScRERERGRmCgZFxERERGJiZJxERFZc8ys28xuNrOnHsZrzzaz3Yscf5KZ3WhmOTP7oZmdFE20IiJLUzIuIiJripn1AP8FnHYYrz0d+CILfv+Z2Q7g68B/AruAfcDXzEy/J0WkpfSQERGRNcPMHgr8HDjxMF77SuCnwP5FTp8H3ODu73H3W4E/BXYAT15FuCIiK1IyLiIia8n/Ab4DPG7hCTM7zcy+b2YzZnanmf2lmVndJU8DXgJcvMh9Hwv8KPzB3XPAtYu9j4hIlDJxByAiItIod/9Y+H19nm1mvcDlwGeAVwInAZcABeCD1df+fvXaly1y66OA+xYc2w8cG1nwIiKLUGVcREQ6wR8CI+7+N+5+h7v/D/C3wOsbfH0fkF9wLA90RxijiMghVBkXEZFOcCpwmplN1R1LAd1m1uXuhRVeP8uhiXc3MBxhjCIih1AyLiIinSAD/BB41SLnSg28fi+wfcGx7cDNqwtLRGR5alMREZFOcDtwCvAbd7/T3e8EHg682d0rDbz+58DZ4Q9m1gecWT0uItIySsZFRKQT/CdBW8nHzexUM3sa8GFgpMHXfwJ4jJm9tTo+8T+A3cD3WhKtiEiVknEREVnz3H0SeAZwPMFIwk8DnwLe2uDrfwP8HvDHwNXAkcDvNlhVFxE5bObucccgIiIiIrIuqTIuIiIiIhITJeMiIiIiIjFZ16MNt27d6scff3zcYYiINO2aa6456O7b4o6jnfTMFpG1arln9rpOxo8//niuvvrquMMQEWmamd0Tdwztpme2iKxVyz2z1aYiIiIiIhITJeMiIiIiIjFJbDJuZt1mdomZjZrZPjN70zLXPsTMvm9mOTP7lZk9r52xioiIiIgcjiT3jF8EnAU8FTgW+IyZ7Xb3S+svMrMB4LvA94EzgGcC/2VmD3f3W9scs4iIiIhIwxKZjJtZP3Ae8Bx3vwa4xszeA7wGuHTB5S8BisDL3b0I3GFmvwM8DlAyLiIiIiKJlchknKDC3Q1cUXfsCuDvzCzj7qW6408Gvl5NxAFw92e3J0wRERERkcOX1J7xo4ARd5+tO7Yf6AIWzmg8EXjAzD5iZveb2bVmpmRcRERERBIvqcl4H5BfcCz8uXvB8Q3AXwFjwDnA54GvmtkjF7uxmZ1vZleb2dUHDhyIMGQRERERkeYkNRmf5dCkO/w5t+B4CbjJ3f/G3a9z93cDlwPnL3Zjd7/E3Xe5+65t29bV5nUiIiIikjBJTcb3ApvMrKvu2HaC6vjIgmvvA3654NjtwM7WhSciIiIisnpJTcavBwoEow1DZwPXLFi8CfAz4BELjj0U+E3LohMRERERiUAik3F3zwGfBj5iZo82s3OBNwIfADCz7WbWW738Y8ApZvZuMzvRzF5PMJv8kjhiFxGRlTWzsZuISCdLZDJedQHwC4LNfD4KvN3dP189dz/wAgB33w08DXgScAtBr/jz3P26tkcsIiKNqt/Y7ZXA35rZC+MNSUSk/ZI6Zzysjr+0+rXwnC34+efAo9sUmoiIrEKTG7uJiHS0JFfGRUSkMy21sdujzCyxRSIRkVbQQ09EEm/PaI53feuX/NP/PZ2Nfdm4w5HVW2ljt/tjiWoNyZfKjEwXGJkuMDVbIlcoM5UvkSuUmMqXmSmUKJSdUrlCsVyhWHZKlQqlss/7vlSp4L78e61wOrimkYsauFMj92ksngbeq6H7tCeWRjT2301E/x23871WvqShi9oZz7t+73ROOXJDA1c2Rsm4iCTeO795G5ffso8XPXonZ5+8Ne5wZPUa3tjNzM6num/Ezp3rZ2Ktu3P/+Cw37R3nnuFp7hnOsXskx57RGQ5O5pnMLxwstrhMysikjWw6RTadIpMKvs+kLTiXSmG28n2sgYsauE2D79XANQ28W2P3acAKN4rs793Qfdr4v0Mjd4rgEmskHmskHmvb/79SDf0fp3FKxkUk0X521zCX37IPgMnZYszRSEQa3tjN3S+hOh1r165d0ZQYE8jduX3/JD+8/QDX3jPK9feO8cDk3L9XNvZmOW5LHw89apBtp3SzdaCLzf3dbO7PsqEnS19XmoHuDH3dGQa6MvR2pcmmraHkTUTipWRcRBKrXHHeedmtbOjOMJkvMTnbWDVQEq+2sZu7F6rHltrYraPdMzzNF6/Zw1eu28ue0RkAjt/Sx1knbuHhO4Z42I4hTtw6oPYskQ6mZFxEEuvq34xwy30TvO05D+XCb9zKhCrjnaJ+Y7cfVo8ttbFbR7pxzxgf/sGdfPuW/ZjB/zl5G6950kk8+SFHcMRgT9zhiUgbKRkXkcTaNxGs7zv7pKBPXJXxzuDuOTMLN3Z7GUFV/I1Ue8M72QMTs7zrf37JV67by2BPhtc++ST+8DE7OWpj78ovFpGOpGRcRBLrQLVn9ogNPQx0Z1QZ7ywXAP9GsLHbBPM3dutIl914P2/50o3kSxVe/aQTedUTT2RDj9pPRNY7JeMiklgHpwpk08Zgb4YNPRlVxjvIchu7dZpiucI7v3krn/7ZPZy5c4j3Pf8MTtg2EHdYIpIQSsZFJLGGp/Js6e/GzKrJuCrjsrbkS2Ve87nr+N9b9/Pysx/Em5/xELoy2m9PROYoGReRxDo4lWfrhi4ABnuyqozLmpIvlXnFp6/mx3cc5B+e81Be9vgHxR2SiCSQknERSayDUwW29Afjpzf0ZDg4VVjhFSLJ4O789Zdu4sd3HOQ9z3sYf/CoHXGHJCIJpc/KRCSxhqfybB0Ik/GsFnDKmvFv/+8uvnzdXi542ilKxEVkWUrGRSSR3J2DU4Vam4oWcMpa8fNfD3PRt2/n3DOO5rVPPinucEQk4ZSMi0giTcyWKJQrbK22qQz2ZpmcLeLesTuiSwfIFUq8+Us3snNzH//8vNO1Hb2IrEjJuIgk0sGpYMZ4fWW8WHbypUqcYYks66Jv3849wzne/byH0delZVkisjIl4yKSSMPVxZr1PeMAEzPqG5dkuuW+cT7109/wkscdx2NP2BJ3OCKyRigZF5FEqlXGq8n4YE9QZZxQ37gk1EXfvp3Bnix/+TsPjjsUEVlDlIyLSCKFyfiWgbk2FUAb/0giXfnrYX54+wH+7LdPZGOvtrgXkcYpGReRRDo4VcAMNvfNbfoDaKKKJI67c9G3b+eIDd289HHHxx2OiKwxSsZFJJEOTuXZ3NdFJh08pjYoGZeEuuruEa6+Z5TXPvkkervScYcjImuMknERSaSDk/laiwrMtalo4x9Jmv/v5/ewsTfL7z9Sm/uISPOUjItIIg1PF2qLN0E945JMD0zM8u2b9/H8Rx6rqriIHBYl4yKSSAen8vOS8f6uDGZqU5Fk+dxVuym78+LHHhd3KCKyRikZF5FEWtimkkoZG7ozSsYlMYrlCp+7cjdPPGUbx2/tjzscEVmjlIyLSOLMFMpMF8rzKuMQLOLUpj+SFD+7a5gHJvO86NE74w5FRNYwJeMikjjhjPFthyTjGW36I4lx2Y33M9Cd4YmnbIs7FBFZw5SMi0jiLNzwJzTYk9UCTkmEYrnCt2/dx1NPPYKerBZuisjhUzIuIokzlgsS7k39C5LxXvWMSzL87K5hxnJFzjn9qLhDEZE1Tsm4iCTOeLUvfGjBtuIberJM5lUZl/h966agReUJalERkVVSMi4iiTOWKwCw8ZBkPMPEjCrjEq9SucLlt6hFRUSioWRcRBJnvJpwDy6SjE/lS7h7HGGJAHDDnjHGckV+57TtcYciIh1AybiIJM74TJH+rjTZ9PxH1IaeLOWKkyuUY4pMBK64YxgzOOvELXGHIiIdQMm4iCTO+EyRob6uQ44P9gSVci3ilDj95M6D/NbRGxf9/6iISLOUjItI4ozPFA5pUYGgTQVgQuMNJSbT+RLX3TvK40/aGncoItIhlIyLSOKMzxTZ2Js55PhANRlXZVzictVvRiiWnbOVjItIRJSMi0jijM8UGeo9tAWgu9pDXixX2h2SCAA/ueMgXZkUu47fFHcoItIhlIyLSOKM5YqHjDUEyGaUjEu8rrjzILuO26SRhiISGSXjIpI44zNFNvYtkoxXK+OlskYbSvuNTBf45b5J9YuLSKSUjItIoswWy+RLlcUr42kDoKDKuMTghj1jADxip1pURCQ6SsZFJFEmZoJJKYsn42pTkfjccO8YZnD6sRvjDkVEOoiScRFJlDEl45JQN+4Z5+QjBhjoPnTSj4jI4VIyLiKJMr5sMh60qRTVMy5t5u7ccO8YZxw7FHcoItJhlIyLSKKM55ZOxrtUGZeY7BmdYXi6wMN2KBkXkWgpGReRRAkr40PLTFMplpSMS3uFizcfrsq4iERMybiIJMpyPeMZtalITG7cM05XJsWDt2+IOxQR6TBKxkUkUcLK+IaepSvjGm0o7Xb9vWOcdvQgXRn92hSRaOmpIiKJMjFTZENPhnTKDjmnTX8kDuWKc/PecS3eFJGWUDIuIokyPlNctF8cIJ0y0inTAk5pq3uGp8kVypx29GDcoYhIB1IyLiKJMpYrLNovHsqmlYxLe/1q/xSA+sVFpCWUjItIoozPFJdPxlMp9YxLW92xfxKAE7cNxByJiHQiJeMikijjM0WGeruWPJ/NpFQZl7b61QNTHLupl37tvCkiLaBkXEQSZXymxOAKbSpawCntdMf+SU45Ui0qItIaSsZFJDHcnfGZlXrG1aYi7VMqV/j1gWlOPkItKiLSGkrGRSQxZoplimVfNhnvSqe06Y+0zT0jOQrlCierMi4iLaJkXEQSI9zwZ6nRhhDswlksqTIu7XFHdZLKKUeqMi4iraFkXEQSYywXJOMrtaloAae0iyapiEirKRkXkcQIK+MrJuMVtalIe2iSioi0mpJxEUmMRpLxrnRKbSrSNpqkIiKtpmRcRBKjocp4RjtwSntokoqItIOScRFJjPGwZ3y5BZwp9YxLe9w/PkuhXOGEbf1xhyIiHUzJuIgkxvhMkZTBQNfS/bnBnHH1jEvr7R7JAbBzs5JxEWkdJeMikhjjM0UGe7OkUrbkNV0Zo6TKuLRBLRnf0hdzJCLSyZSMi0hijM8UGVqmXxw02lDaZ/dIjmza2D7YE3coItLBlIyLSGKMzRSXXbwJYTKuNhVpvd0jOY7d1Ed6mU9qRERWS8m4iCRG2KaynGzaKKgyLm1w70iOHZvVoiIiraVkXEQSY2KmyFBf17LXqE1F2mX3SI4dm3rjDkNEOpyScRFJjPGZIht7l9/pMJtOUVKbirTY+EyRsVyRnaqMi0iLKRkXkURw92oyvnLPuNpUpNXurY01VDIuIq2V2GTczLrN7BIzGzWzfWb2pgZes9nM9pvZy9oQoohEaCpfolzxFZPxrnSwA6e7quNrnQW+Y2aviDuWhfaMBsm4esZFpNUSm4wDFwFnAU8FXgn8rZm9cIXXvB84otWBiUj0xmeC3TeHepfvGc+kU7hDuaJkfC0zsxTwAeBpcceyGM0YF5F2SWQybmb9wHnA6939Gnf/GvAe4DXLvOaZwKOBA+2JUkSiFCbjK09TCR5bGm+4dpnZMcD3gHOBsZjDWdTukRxDfVkGe5b//6OIyGolMhkHzgC6gSvqjl0BPMrMDlndZWYbgI8C5wOFtkQoIpEazwXJ+Mo948HM52JFfeNr2JnAXcAjgfGYY1nU7pEZ9YuLSFskNRk/Chhx99m6Y/uBLmDbIte/B7jc3X/UjuBEJHphZXzFnvFMtTJeUjK+Vrn7N939Fe5+MO5YlqIZ4yLSLsvPEItPH5BfcCz8ubv+oJk9EXgOcFojNzaz8wkq6OzcuXN1UYpIZGo9431qU1nrzKwHOHaJ0/vdfbKJe7X9mV2uOHtGczzjt7a35f1EZH1LamV8lgVJd93PufCAmfUCHwde6+4NfdTp7pe4+y5337Vt22JFdhGJw1iDlfFMdWtybfyTaLuAO5b4el4zN4rjmX1wKk+x7BwzpA1/RKT1kloZ3wtsMrMudw97wLcTVMdH6q57NHAS8BkzC4/1AR81s8e6+6vaFbCIrM74TJFMyujrSi97XdimolnjyeXuVwC24oUJtW886JDcPtgTcyQish4kNRm/nmAh5lnAD6vHzgaucfdS3XVXAScveO2PgYuBT7U2RBGJUrjhT90/rBcVtqloF05plX0T1WR8o5JxEWm9RCbj7p4zs08DH6lu4LMdeCPVvkEz2w6Mu/sMcGf9a82sDDzg7g+0N2oRWY3xmSIbV+gXh/qecVXGpTX2V5PxIwYXdkuKiEQvqT3jABcAvwC+TzC28O3u/vnqufuBF8QVmIhEbzxXXLFfHOZGG6pNRVpl3/gsmZSxtV/JuIi0XiIr4xBUx4GXVr8Wnlvyc2x3X2oFv4gk2PhMkS0Dy+++CXWVcY027AjufnzcMSy0b2KWIzZ0k0qt2bZ3EVlDklwZF5F1JOwZX4lGG0qr7Z+Y5Uj1i4tImygZF5FEGJ8pMtREm4p24JRW2Tc+q0kqItI2SsZFJHaVijMx22RlXG0q0iL7J/IcqWRcRNpEybiIxG5ytoQ7DKpNRWI2lS8xlS9prKGItI2ScRGJ3Xh1982hvkYWcGoHTmkdbfgjIu2mZFxEYhcm4820qWi0obRCOGNcbSoi0i5KxkUkdqO5AgBDDWz605XRDpzSOrXKuNpURKRNlIyLSOzCZHyTduCUmO2bUJuKiLSXknERid1YrvGe8Yx6xqWF9k/MMtiTobcrHXcoIrJOKBkXkdjV2lQa6BnvUs+4tND+iVm1qIhIWykZF5HYjeWKbOjJkEmv/EiamzOunnGJ3j7NGBeRNlMyLiKxG8sVGlq8CZBOGSmDknbglBbYr903RaTNlIyLSOxGc0U2NdAvHsqmU2pTkciVK86BqbzaVESkrZSMi0jsgsp4c8m42lQkaiPTBcoVZ9uG7rhDEZF1RMm4iMQuqIw31qYCwS6cmqYiURuZDhYSb+5v/B+GIiKrpWRcRGI3mis03aainnGJ2vB0HoAt/aqMi0j7KBkXkViVyhUmZ0sNL+CEas+42lQkYsNTQWV8y4Aq4yLSPkrGRSRWYzPBhj/NVMa7Mim1qUjk1KYiInFQMi4isRoLN/xpojKeSalnXKI3PJXHrLl/GIqIrJaScRGJ1Wiu+cp4Nq3KeJTMbHPcMSTB8HSwdiGdsrhDEZF1RMm4iMRqtNoa0FQynklRLKtnPEL/EncASTAyXWCLWlREpM2UjItIrMaqlfFm2lS6NNowauea2Y64g4jb8FRB/eIi0nYNJeNmNtTqQERkfRqt9oxvaiIJUptK5PYB/2Bm2+MOJE7D03lNUhGRtmu0Mv7ulkYhIuvWaK5INm30d6Ubfk0mnaKgNpUonQv8GfB6Mzs97mDiMjxd0IxxEWm7RpPxF5vZJ8zsxWZ2dEsjEpF1ZSxXYKivC7PGF811pY1iSZXxqLj7ne5ecPe3AOeY2dPjjqndSuUKY7mi2lREpO0aTcbvB7YD/wbca2a3mdmHzex5WoUvIqsR7L7ZeL84aAfOVnL3dwNHmNkrFp4zs0wMIbVFONVnq9pURKTNGk3Gv+Pu5wCbgCcCnwNOAz4LPGBmV5nZM1sUo4h0sNFckaEm5zoHPeNqU2mhzwIPMbOPmdnZZhb+a+miOINqpeHpPACb1aYiIm3WaDL+VQB3L7n7Fe7+Dnf/bWAIeCbwU+DTZvaI1oQpIp1q7DAr4wW1qUTGzN5Z/fNMM7sYuA+4APgj4EfAmJl9Dzgnvihba2RKu2+KSDwaSsbd/TtLHJ919/9199cDjyF4eIuINGw0V2Sot9nKuEYbRuxlZnYTcDVwHvB94NnARmBH9dhdBO2KHelgdd692lREpN0i6/9z97vN7Kqo7icinc/dgwWc/c1XxpWMR2o7cCtBG8qX3H267txegtbEzzWzyHatGZkK21SUjItIe0WSjJvZzcAngJ4o7ici68N0oUyx7E3tvgnVBZzqGY/Sl939Dxq47qstjyQmI9MFUkbT6xdERFYrqh04fwD8HnBlRPcTkXVgtNoa0HTPeMYoqDIepYsbucjdv9XqQOJycLrApr4u0qnOrf6LSDJFUhl399dGcR8RWV/GZ4Jxcs1WI7vUphIpd/9Z3DHEbWSqoBYVEYlFVJVxEZGmDVcr41uaTIIyqRQVh3JFrSoSjeHpPFu0eFNEYqBkXERiM1xdNLdloLnZztlM0Eqg6rhEZXi6wBbNGBeRGCgZF5HYDFdnOzdbkexKB48uJeMSlZHpgirjIhILJeMiEpuD03m60ik2dDe3fCVbS8bVpiKrVypXGMsV1TMuIrFQMi4isRmeCqqRzc6vzqoy3hJmttPMFh1tY2ZZM9vZ7pjaobaQuLe5qT4iIlFoKhlfrw9qEWmN4anDWzSXSQfJe6GkZDxidwNnLnHuEdXzHSdMxjc2OWJTRCQKzVbG1+WDWkRaY3i6wNYmF2+CesZb6O3AniXO7ame7zi1ZFyVcRGJQbPJ+Lp8UItIawxPHd4Ei7BNpaTRhpFy9wvd/b4lzu119wvbHVM7KBkXkTg1tWpquQexu+8FOvJBLSLRc3cOTuXZehhtKlm1qUTGzL7ezPXufm6rYomLknERiVMkO3CKiDRrulAmX6ocVs+4FnBGahCo/4jhLKAC/AzYB2wBHkvw++IbbY+uDSaqyfigknERicGKybiqJiLSCrUNfw6jTSWdCirj2oFz9dz9t8PvzexNwBDwTHffV3d8E0EivlSb4pqmyriIxKmRnvFBYEPd19OBpwEDwBTQDTwReAqQa02YItJpDh7mhj8AGSXjrfJG4O/rE3EAdx8F/hl4eSxRtdj4TJGebIruTDruUERkHVqxMq6qiYi0QlgZP5xpKqqMt0wa2LzEuR1AoY2xtM34TFFVcRGJTbPTVNZl1UREojc8vYrKeHUBp6apRO7LwEVm9vtmNgBgZoNm9jLgXcBn4gyuVZSMi0icml3AuS6rJiISvbAyfjhbkKdMlfEWeT1wFPAFwM2sCGQBA/4TeEuMsbWMknERiVOzyXhYNZkGLnf3KTMbBH6PoGry8agDFJHOdHCqwIaezGH16WZSwYd6Ssaj5e7TwLPN7GHA4wnaEoeBH7r7r2INroXGZ0ocM9QTdxgisk41m4yvy6qJiETvcHffhLmecbWptIa73wjcGHcc7TIxU+TUozbEHYaIrFPNbvqzLqsmIhK94ak8Ww6jRQXmesZVGY+WmaWA84BzgH4OXVfk7v6UtgfWYmpTEZE4HdamP+utaiIi0RueKnD81r7Dem3YM16qaNOfiL0PeB1wHcF0rI7/L7hUrjCVLykZF5HYNJWMr9eqiYhEb3g6zyOP33RYrw3njFdclfGIvRi40N0vjDuQdpmYLQHa8EdE4tNsZXzdVU1EJHrlijMyXWDrYbap1HrGy0rGI9YN/CjuINpJu2+KSNyaTcbXXdVERKI3litQcdhymAs41TPeMpcBzwV+EHcg7aJkXETi1mwyvu6qJiISvdVs+AOQNk1TaZHLgYvN7ATgSiC34Ly7+8XtD6t1lIyLSNyaTcbXXdVERKL3wESw4c9qRxuqZzxyn6z++azq10IOKBkXEYlQs8n4uquaiEj09k/MArB98PA2Wgk3/VHPeLTcfeGi/I6nZFxE4tZsMr7uqiYiEr191WT8yMNMxtPqGZeITFST8UEl4yISk2Y3/Vl3VRMRid4DE7MM9mTo7Uof1uvVMx4dM7sA+Ky7769+v5yO+/RzfKZITzZFT/bw/r8oIrJaKybj6/1BLSLR2zcxy/aNh1cVh7me8bI2/YnCe4ErgP3V75fTcZ9+jue0+6aIxKuRyvi6flCLSPT2T+QPu0UF5jb9KSsXX7X6Tzzb9emnmR1B8Lvidwh+b3wTuMDdx9rx/vXGZ5SMi0i8VnzwunvK3a+q+365L33OJyIr2j8xu6pkPJUyzJaujP9y3wRPeM8PGK2OUJTGmNnpy5wbMrNLInqrzwHHAk8j2NH5dOA/Irp3U5SMi0jcmqqCtPFBLSIdqlJxHpjMc+Tg4Y01DKXNluwZv/OBKXaP5Ng7NrOq91iHfmhmj1p40Mz+GPgVwcZvq2JmxwJPAc539+vd/WqCnZ2fa2Z9q71/s5SMi0jcmv1IsuUPahHpbAen85QrvqrKOAR940tNUwmP50vqY2nS94DvmtkTAMzsFDP7PvBp4CqCCvZqjRNM47qj7pgT/D5a3b/QDieYmaImqYhIrJpNxtvxoBaRDhZu+LPaZDyzTDIebgZUUDLerBcA/wn8T/WTzhuA44Dfdfdnu/tdq30Dd59092+5e/3/OK8DbnL30dXev1kTqoyLSMyanTP+AuBDBA/qzwJ/DNxH8KD+RtTBiUjn2b/KGeOhdGrpNpVwYWdBKzyb4u4OvNrMHgDeBnyb4PneVPO9mfUQ9IQvZr+7T9Zd+wbg+cDTl7jX+cD5ADt37mwmjBWVK85kvsRgj5JxEYlPs3PGI3lQi8j6tW+Vu2+GlmtTqVRUGW/UEiNrJ4FrgacCF5rZgerxRsfX7gJ+vMS5PwE+VX3vvwQuAl7r7t9d7GJ3vwS4BGDXrl2RDpafypcA2NDTbF1KRCQ6jc4ZX2i1D+oVmVk38EGCikke+Bd3f88S174A+DvgQcCdwN+qUi+STPsn8qQMtg50reo+6VRq6cq42lSasdLI2jfXfd/Q+Fp3vwKw5a4xswuBvwf+wt0/vNI9W0HJuIgkQaNzxpfT9IO6QRcBZxEk/McCnzGz3e5+af1FZvZ/gIq3WOMAACAASURBVM8ArwZ+QDAm68tm9mh3vy6iWEQkIvvHZ9k60E0mvbqR1pmU1SrgC4UV80K5vKr3WA/i2FnZzF5HUEB5ZbXyHYvpajLe361kXETis+ITKKYHdT9wHvAcd78GuMbM3gO8Brh0weUvBb7k7v9e/fkDZvZsgv52JeMiCbN/cnUzxkPL9YxrAWdymdlO4N3AvwFfN7PtdacPuHvb/gU1ORsk4wNKxkUkRm1PtBt0BsGIqyvqjl0BPMrMFj41Pwi8Y8ExB1b/215EIrdvPLpkfKlNf9QznmjnEjzf/xy4f8HXg9oZSFgZVzIuInFKajJ+FDDi7rN1x/YDXcC2+gvd/QZ3vzX82cxOI9hQ4kftCFREmhPFhj8QtKks3TMe/FkoR7reTyLg7h9yd1vi6852xhL2jA+oZ1xEYpTUZLyPYNFmvfDnJX+Lm9kRwFcIVvF/dYlrzjezq83s6gMHDix2iYi0SL5UZmS6sOpJKhBUxsN2lIVUGZdGhMl4f5eScRGJT1KT8VkOTbrDn3OLvaC6xfIPgTLw+ws2lKhx90vcfZe779q2bdtil4hIi0S14Q9Ue8aXqHxrmoo0YmpW01REJH5NPYGqC2/ud/fiIud6gIe7+88jiGsvsMnMuupmmG8nqI6PLPLeJxDsDpoDnuTuwxHEICIRC2eMH7lx9cl4Jr30nHFNU2mcmX29mevd/dxWxdJumqYiIknQbGX8buDhS5x7DPD91YVTcz1QIBhtGDobuMbdS/UXmtlm4H+BceCJ7r4/ohhEJGJ7RoMPtnZs6l31vdK2zDQVtak0YxDYUPf1dOBpwAAwRfCp5BMJ1uIs+snkWjWVL9GdSZFd5ZhNEZHVaGTTn48CR4c/Au8zs7FFLj0VOBhFUO6eM7NPAx8xs5cRVMXfSHVL5OoorHF3nwH+EdgKPA/I1I3JmnH38SjiEZFo3DsyA8DRQxEk48v0jKtNpXHu/tvh92b2JmAIeKa776s7vgn4BrCn7QG20FS+pEkqIhK7RsoB32KuYgLQz/wqygaCBZc3AC+MMLYLgF8QVNs/Crzd3T9fPXc/wRxxCHboHCSYKV4/IiuWHd1EZGl7RnMcsaGbnmx61ffKpFJL9ozXKuNlJeNNeiPw9/WJOIC7jwL/DLw8lqhaZCpf0iQVEYldI5v+fB34OoCZ/QD4c3e/rdWBuXuOYEOfly5yzuq+39rqWEQkGveOzLBjc18k9wrmjC9fGc+rMt6sNLB5iXM7CNoHO8Z0vqRJKiISu6Ya5dz9SQsTcTN7lJk9r9q7LSKypD1jOY6NoF8cwh04F0+2w4K42lSa9mXgIjP7fTMbADCzwWq74LuAz8QZXNQmZ1UZF5H4NZWMm9lOM/uJmb2t+vOrgJ8D/w3cYWZntiBGEekApXKF+8dmI03Gl9rTp6Ke8cP1euAq4AvAuJnNAqPAJwg+IX1LjLFFbrqgnnERiV+zT6H3AkcC3zezLuCfCBb1/AVBX/dFwFMjjVBEOsK+iVlKFWfHpmjaVDIpo7xkZVw944fD3aeBZ5vZw4DHEyzmHAZ+6O6/ijW4FpiaLXHCViXjIhKvZp9CTwFe4e4/NrNnEDyo/9Xdd5vZ+4EvRR6hiHSEPaPBJJVjI0rGl930p5qMF5WMHxZ3vxG4Me44Wm0qX9aMcRGJXbNPoSxzm+48C5gk2Ho+PNdRi3tEJDpzyXiEbSpLLOB0takcFjNLAecB5xBMzlrYyuju/pS2B9YiU/midt8Ukdg1+xS6Djiv2kf4QuBb7l4ysy3Am4Grow5QRDrDvSM5zKKZMQ5hz7jmjEfsfcDrCJ71e4CO/S+wVK4wW6xomoqIxK7Zp9CbgMuAPySokL+9evzW6p/PiCguEekwe0Zn2D7YQ1cmmt0OM8uNNqymkBpt2LQXAxe6+4VxB9Jq0/kygKapiEjsmnoKufuVZrYDOA34pbtPVU+9BLiqujGEiMgh9oxGN9YQIK1Nf1qhG/hR3EG0w1ShBMBA9+o3oBIRWY2mS1TuPuPuV9cl4rj7t5WIi8hy9ozORDZJBSCdYsVNf9Sm0rTLgOfGHUQ7TM2GyXg25khEZL3T53Mi0nLFcoX7x2cir4wv1TNeq4wrGW/W5cDFZnYCcCWQW3De3f3i9ocVval8kIz3qzIuIjFTMi4iLbdvfJaKRzfWEFboGXe1qRymT1b/fFb1ayEHOioZ1zQVEYnbik8hM5to4n7u7htXEY+IdKC7D04DsHNLlG0qRmmJZLusyvhhcfdoVteuAdO1yriScRGJVyNPoRcDnwGKwIcIKiMiIg2760CwxOTEbQOR3XO5OeMV9YwfFjPrdfeZJc6lgCF3H1ns/Foz1zOuZFxE4rXiU8jdv25mzwS+Dxxw9w+3PiwR6SR3HZhisCfD1oGuyO6ZSRmlJUcbBsdLFadScVIpi+x9O5GZvRF4I7DNzPYA73b3jyy47FHAT4GOaLIO21SUjItI3Br6SNLdfwr8PfB2MxtsbUgi0mnuemCaE48YwCy6pDidsloFfKH67hX1jS/PzF4NvAv4IvB64FfAh8zsUjPr2Ex1Sm0qIpIQzTyF3g/cTrBFcjN95CKyzt11YIonnLIt0nsuVxmvT9IL5Qo92Y4o5rbKnwPvcPdwE7cPmtkrgH8Dsmb2fHfvuH/RTOdLdGdSZNPrpk1eRBKq4aeQuxfc/Wvufn8rAxKRzjIxW+SByXyk/eIAqZThPjfGsN68ZFx94ys5Dvhx/QF3/zjwMoKZ4x+PIaaWm8yXNElFRBJBTyIRaalfHwgmqZx0RLTJeKbaB16qOF0LesLrF3YqGV/RbuAxwA/qD7r7Z83sCOB9ZjYKfCGO4FplOl9Si4qIJIKeRCLSUnc9EE5S6Y/0vulU8MHeYn3jqow35ePAO82sB/iKu98QnnD3i81sG/AW4ClxBdgKU7MlLd4UkUTQk0hEWuquA1Nk08aOzdHNGIf5lfGF5lXGtYBzJe8HBoE3AJuBv6g/6e5/Y2b7gX+OIbaWmVJlXEQSQitXRKSl7jowxXFb+iNfKBeOKyyXF6mM109TUWV8We5ecfd/AIYIpmYtds2/Ag8C/rSNobXUVL7EBiXjIpIAK/52NLNrzezM6vcvMbMtrQ9LRDrFXQemI29RgfrK+KHJdrmuTSWvZLwhHhhb5vw+d/90O2NqJfWMi0hSNFKqOhU4pvr9J4ETWheOiHSSYrnCPcPTkU9SgWDOOMxPvEPlipNNB+dVGZfFTOVLDGiaiogkQCNPouuAS83sbsCAz5nZotslExRXzogsOhFZ03aP5CiWnRNakIyHlfHyEqMNe7JpiuVSUz3jpXKFuw9Oc/KRGyKLU5JpKq8FnCKSDI08iV4AvI5gYc9DCTb+OdDKoESkM9x2f7A/2EO2R5/chj3jpUV6xssVpzebZnK2RLGJyvi3bt7HGz5/Pb9461PZ3N8VWaySLKVyhdlihf4uJeMiEr8Vn0Tufi/wRgAzexLw1vrRVyIiS7l57wTZtHFKCyrNy1XGyxWntyvYdbOZyvhYrkC54kzNltZNMm5m1wIvd/frzOwlwGXuPhx3XK2UK5YB6OvSzqwiEr+mxhu4+4Pc/QYLPNTMHmtmJ7cqOBFZ2265b5wHb99AVyb6wU3L9YxXPKiMw/ye8fvGZnjzF28kXyoves/w2nU2DnHdrQuaLQT/+/cqGReRBGj6N6SZ/SlwP3AT8FPgl2Z2v5m9MurgRGTtcndu3jvOaUdtbMn9M9VNf5aqjHcvkoz//NfDfP7qe7lj/9Si9wyT8OL6SsbDdUE3Mbcu6MYlvjriU9FcQZVxEUmOphrmzOxFBLu1XVr92g9sB14IfMTMxt390sijFJE15/7xWUZzRX7rmMGW3D+9TM+4O/Rmg2Q9X5dYh9cOTxcWvWexFJxfZ8n4ulsXFCbj4acnIiJxanb1yl8DH3X3P19w/GtmNgK8iSBJF5F17ua94wCcdkxrKuPp5XrGl2hTKVZnkg9P5Re9Z6EcJGnrKRlfj+uCZopqUxGR5Gi2TeVk4MtLnPsq8JDVhSMineLm+yZIGZy6vTWV8cwKc8ZrCzjrkvEwcR9ZqjJerZwXSofecz1YL+uCZmptKpqmIiLxazYZvwc4fYlzDwM6egW+iDTulr3jnLhtoGXVx7nK+KFV7ErF6cksUhmvJtsHpxZPxsNrF9vVc71YD+uCcoUSoJ5xEUmGZssCnwLebmaTwBfdfczMhoDnA/8AfCja8ERkrbrlvgked+KWlt0/s9yccXey6RTplNVaT4JrV2pTWZcLOGvWy7qgsE2lRz3jIpIAzSbj7wXOAC4BPmZmpeo9DPgS8PfRhicia9HBqTz7JmY57ejWtKjA3KY/i09TCc53pVPzKuOlykoLOKujDddpmwrrZF3QjKapiEiCNJWMu3sJeJGZ/SPwBGATMAJc4e43tSA+EVmDrrlnFIAzdgy17D2W6xmvuJNOQTZt85PxcJrKEpXx4jqvjBOsC7pgiXNfBf6kjbG0jEYbikiSHNbqFXe/Gbg54lhEpENcdfcI3ZkUDzu2NZNUoG604RJzxtNmdGXSFOraWMJe8CV7xpWMh+uCvrvIuY5ZF6Q2FRFJEi0lF5HIXXX3CGfuHKI707pkp7bpzyI945WKk0oZ3ZnUogs4h6fzuDtmNu91hfU5Z7zep1gH64JmCmVSBt0t2BlWRKRZehKJSKQmZ4vcct84j35Q6xZvAlRz8cUr4x5WxlPztrYPJ6/MFiu1VoV6YRJeWCTBXyfeC3yTYF3QsJnlCarhHwMup0PWBeUKZfq6Mof8Y0xEJA6qjItIpK65Z5SKw2MetLml7xNWxitLzBlP1xZwziXdxboke2S6QH/3/EdgWEUvltZnZXy9rAuaKZa04Y+IJIaScRGJ1JV3j5BJGY/Yuaml77Ncz3jFgzaVrszCaSpz3x+cyrNjc9+812kBZ6DT1wXlCuXaDq0iInGLrE3FzLJmtjOq+4nI2nTV3SM87NiNLa88Zpbb9MdZtE2lfib58CKLOJWMrw8zhbImqYhIYkTZM/4I4O4I7ycia8xMocyNe8Za3i8OdZXxxTb9qS7gXGzOePVlDE8fOt4wX1r3PePrwkyxrDYVEUmMKJPxPcDbI7yfiKwxP797mGLZeewJre0Xh7lkfGHPeKXatpIyDm1TKVfYtqEbWHy8oSrj60NOlXERSZDIknF33+vuF0Z1PxFZe7576376utI89oTWV8YzS/SMh5sAhW0q+frRhhWnvzvDQHdmiTaV6mjDdbqAc72YUc+4iCSIRhuKSCTcne/d9gBPOHlbWzZTSdd6xhck42FlPGxTqR9tWHYyKWPLQNeibSphFX2xRaHrXSetCwraVDS/QESSoamnkZmlgPOAc4B+Fknm3f3J0YQmImvJLfdNsG9ilqc+9Mi2vN9SPeNh20p6iWkqmVSK/v4MI9NLt6kU1KaymEcAPwXWfEk5VyjRp8q4iCREs6WB9wGvA64j6BHXbywRAeB/b91PyuBJD97WlvdbqTKetqAyXt//XSw72bSxpb+bvWMzh9wzTMLVprKojlkXlCtoAaeIJEezyfiLgQvVGy4iC333tv088rhNbBnobsv7hZv+lA9ZwBn8udSc8Uw6xdaBLm7cM3bIPWub/qgyfgh33wt0xLN/VtNURCRBmu0Z7wZ+1IpARGTtunckxy33TfDUU9vTogLLVMZrCzgXm6YS7My5ZaCLkelCbfJKaG6ainrGO1WxXKFYdrWpiEhiNFsZvwx4LvCDFsQiImvUl6/dixmcc/pRbXvPpXrGa20qqUU2/ak4PdkUW/q7KVWcidkiQ31dtdeFufl67RlfD+uCcoUygCrjIpIYzSbjlwMXm9kJwJVAbsF5d/eLI4lMRNaESsX572vu5fEnbj1ke/lWCjfvWbgDZ7iAM5ymUiw7leomQKVyhUx3hqG+LACjublkvL6Cvo7bVDp+XdBsMUjG+zRNRUQSotmn0Serfz6r+rWQA0rGRdaRn/96mD2jM/zV0x/c1vc1MzIpO6RnfN4CzkxQ2C2UK/Sk0rUFnOHoxfoEvFBWMs46WBc0VxnXZF8RSYYVn0Zmdq2ZnVn98U+AI9w9tcSXPvcTWWe+cPW9DPZkePpp29v+3umUHTITvL4y3l2XjEOQqKerFXNYuhpeLK3bnvGOXxeUK5QA6M2qMi4iydBIaeBU4Jjq958AHtS6cERkLRnPFfmfm/fxuw8/pi0b/SyUThnlhXPGqzn1vMp4OCWlOk0lPJ4vlWuvW6pKvs6E64JaysyONrMvm9m4me0zs3ebWVuy45lC2Kai2pGIJEMjD7/rgEvN7G7AgM+Z2aEDeqvc/WFRBSciyfbpn/2GfKnCHz4mno0ZF6uMl+s3/VlQAS+VnWx9xXypyvj6TcbbtS7oC8A48FhgG/BZYAL4xwjuvayZohZwikiyNJKMv4BgQc9m4KHA7cCBVgYlIsk3nS/xiZ/czVMecgSnHjUYSwyZlNXaUkJhz7gZZBck40GbSl1lfIkEfB0n4y1fF2RmG4B7gTe7+27gNjP7b+CJtCEZr/WMa7ShiCTEism4u98LvBHAzJ4EvNXdb2h1YCKSbJ+7cjdjuSKvfvJJscWQTqWW7BlPp4yUBSNXajtrlitk03PtK/niXNKdrybs3ZnUqueMf+eWffzjt27jgqedwrlnHI1V40giM7sWeLm7X0ewLugydz/Yqvdz90ngRXXvfxpwLvDvrXrPempTEZGkaWo5ubs/SIm4iMwWy1zy419z1olbeMTOTbHFkU5xSM/4otNUwjaVipNJG92Z6jSVedXw4HX93ZlVV8avv3eMe4ZzvO7S63nD56/HPdELQmNbF2RmPwFuBsaAD7XjPWc02lBEEkaznUSkaR/54V0cmMzzuqecHGscmUUq42EynkrNJePFusp4JpVatme8ryvdcDK+ZzTHx3/860OOT86W2Nib5YWP2sFXr7+P6UJ5kVcnRrgu6Cbm1gXduNRXIzc0sx4zO2mJrw11l74aeArQC/zXEvc638yuNrOrDxxYfYek2lREJGmUjItIU35zcJqP/r+7OPeMo3nMCVtijSW9SM94rU3FDl3AWa44mbokfbFpKv1dmYbbVL52/X2887LbGJ0uzDs+OVtkY2+29qnBwvMJ8wLgo8AvCHrCbweuWearEbuAO5b4el54kbtf7+7fB14OPMfMjl94I3e/xN13ufuubdu2HcZfb76ZcLSh2lREJCH0OZ2INMzdufAbt9CVTvHWZ50adzhkFpumUqmbplKrjAfHSmUnk168Mh62rPR1pymONVYZH58pAjAxW2RTf1ft+ORsiQ09czt9juWK7Njc9F+vLVqxLsjdryCosh/CzDab2Qvc/fN1h2+t/rkV+M1q3nsluUJ53j/IRETipqeRiDTs0l/cyw9uP8AbnnYKRw72xB1OMGe8Mj9xrt/0pzZNpRxUwIuVBQs469tU6irjjc4ZH89Vk/GZ0rzjE7NFNvRkagn62ExQGf/fW/fz3A//hCt/Pdz4X7KN2rQuaDNBW8wj6o49EigDv2rxezNTLKsqLiKJsmJl3MwuaOJ+Uc2gFZGEuXnvOG/7+i084ZRt/MlZx8cdDlCdM37IAs7quXltKk6l4rgvPn8c5irj/d2N94yHlfHwz9DkbIkdm/sY6g0q46PVpP0ndx7k+nvHeOG//5zXPvlkLnjaKc38dTuCu99pZpcDHzOz84CNwCXAB919otXvP1Moq19cRBKlkTaV9zZxv1XPoBWR5BmeyvPqz13L5r4uLv6DM0ilkjGqb7Ge8bkFnNCVmRttWKxW0LPpFJl0inTKFl3A2d+VoeLhTPLl/571bSr15tpUqpXxXFAZH5kucPTGHh6+c4gPfO8OXvq449gy0H1Yf/c17o+A9wPfAyrAZ4C3tOONc4WyxhqKSKKs2Kbi7qkmvvSEE+kw4zNFXvKJq9g3PsuH/+jMRCWPi/WM+7wFnMEjqViq1CromWqC3ZVOzVvAWSwF5/u6q69poDpeS8YXVMYnZosM9mTn9YwDDE/n2b6xh5c+7nggGIG4Hrn7iLu/xN23uPs2d7/A3duyyjVoU9FyKRFJjsPuGTez48zssWbWb2YDUQYlIskwOl3gZZ+8il/tn+Rjf/xIHnlcslYhBj3jCyrjdZv+ZOsq47VkvNqi0p1NzauM52ttKpnaa1ayWGW8UnGm8kFlPJtOMdCdYbRaGR+eKrBloJvTj91IOmVct3t9JuNxmlFlXEQSpunygJk9D/hn4ESCjxcfDbzNzCaBP3H34nKvF5G14c4HJnn5p6/m/rFZPviiR/DbDz4i7pAOsXjPeN2c8fTcnPFStU2lvjI+b9OfugWc9T8vZ2Lm0AWc04US7jDYE1TFh/qydZXxAmfuHKKvK8OpR23guntHm/wbR2+9rQvKFUra8EdEEqWpJ5KZ/QHBxgyfBP4a+EL11FeADwN3A38XZYAi0l7uzheuvpd3fPM2erIp/uv8xySuIh5Kp6xWCQ/VzxnP1o0wDNtZMulqMp5JkS8uvulP8PPys8bLFWcyHyTh9ZXxydng2Iae4PG6qa+LsVyBSsUZmS6wuTph5cwdm/jKdXsb6k1vsXW1LihXKCeq1UpEpNk2lb8H/tXdX0GQgAPg7p8C/pZgUY6IrFG33T/BH//HVbz5Szdx2tGDfO01Zyc2EYdgB85D2lSq+XWqfppKuVJLtrOpaptKJlVrTYG5ySoD1TaVlXrGJ+sS8Pqe8TAx31BXGR/NFRmfKVKuOFv6g0TwzJ1DTOVL3PHAZBN/4+itt3VBs0W1qYhIsjSbjJ8EfGuJc9cBR60unDlm1m1ml5jZqJntM7M3LXPtGWb2MzPLmdk1ZvaoqOIQ6XTuznW7R3nN567lnA/8mBv2jPGO3z2N/zrvsRwz1Bt3eMtKL7PpTyrF3JzxUqV2fK4ynl50mko4g3qlnvHxeQn4XJvKwsr4ULUyPlzdhXPLQLUyXt2dM6l94526LkjTVEQkaZptnNsNnA18d5FzjwbuXXVEcy4CzgKeChwLfMbMdrv7pfUXmVk/8D/A54E/BV4JXGZmJ7p7vCUnkYRyd+4ZznH5Lfv4xg33cct9E/R3pfmzJ57IK59wIhurU0CSbrlNf9Ipq30Vy5Va20nYEtKVSc3b9KdQdrrSqVo1fWEv+kLzkvGZ+jaVsDIetqlkGZspMjyVB6hVxo/f0semvizX7R7lRY/e2eTfvHU6fV3QTKFMj+aMi0iCNJuMfwh4r5kZQYXcgWOqO6m9FXhHFEFVE+zzgOe4+zXANWb2HuA1wKULLn8BUAT+0t0rZvYG4FnV4x+PIh6RtW46X+LXB6a5dvcov/jNCFf/ZpR9E7MAnHHsRt7x3N/i/555TK1FY60IkvH5x8IKeNrqFmqW5hZwhtXy7kyKQt1ow0KpQlcmVZu2slKbSpiMbx3oWqJnvNqm0ptlfKbIgTAZr1bGzYwzd25KVGV8PawLmlGbiogkTFO/ed39A2a2CXgzQY+4AV8jSIY/4O7NLARazhlAN3BF3bErgL8zs4y71+89/VjgJ+5eqcboZvYT4HEoGZcOVihVmM6XmKp+Tc6WGJnOc2Ay+HpgMs9vhqe5++A0+yfytddtH+zhUQ/azKOO38STHnwEOzb3xfi3WJ3MMpXxcGOibNoolv2QOePdmRRT+blHSbFcIZs2sum5cYjLCZPxYzf1sW98tnY8bFkZ7J1rU3GHuw9MA3PJOMCZO4a45p5RZouJqdaG64IuMLNaQO7+qeqz/7Ws4WQ8mKrj2oFTRBKl6TKYu19oZu8nSIK3AOPAle5+MMK4jgJG3H227th+oAvYBty/4NrbF7x+P/DwCOMBYPdwjo/88M6obxspX/6T9da+N/G9ebx/7+D9Haf6Hyru1WPV5LB63n3+OV/wfanitZaKUrlCoRz8XKoeCxci5grlef3OC6UMNvd3s3NzL2eftI0TtvXzoK39POzYjRwz1ItZMnbQXK3lesZrlfFMmvxi01TSqUN6xrN1bSorjTYMk/Edm/v41f65jriwTSUcbbipP/jzrgNTwc99c8n4eU84gVc/6aTE7GhKsC7o9Uuci3RdUBzCtqTujJJxEUmOw/pM2t3HgW8DVBdLPtHMfuDuIxHF1QfkFxwLf144k2qpaxedXWVm5wPnA+zc2Vyf5sRskR/c/kBTr4mDEd8v9jhzvDjTmTC5NQumeJgF8aSq31j1mpQF//uYBT9b9TVWPZ5JG9l0ip5simx105hs9VgmlaIrY2RSKfq602zoztDfnWGgO8OGnuD7zf1dHLGhh839XXGPy2uLRTf9qSzoDU9b7R80EExggUM3/QnbVMJxiCuNNqwl45t6yRXKtWR+YqZENm10V+8z1Bsk33cdmGaoL1trkwGSUg2v1851QW03Wwzaknqyh73fnYhI5JqdM76ToJ/wO9UK+asI+ggNGDWzp7r7dRHENcuhyXT4c67BaxdeB4C7XwJcArBr166m6qm/dcxGrvybpzbzEhFpocWS8YVtKl2ZIOkulg+tjM9fwFmhK52qJcuN9Ix3pVMcsSF4/EzMFNky0M3kbJENPdnaP9CG+uYq40dt7FnV37cN2rIuKC6qjItIEjVbHngvcCTwfTPrAv4J+AZwPHAlwQSUKOwFNlXfI7SdoOK9sPq+t3qOBdfej4h0tMyiyXjwZ9imkk2nKJbnRhuGyXaYpIfCynajPeMTM0UGe7O1yTNhr/jkbKk2SQWCnnFYG5vNuPsHgH8E/gr4KXPrgj4MXBLhuqBYhJXxblXGRSRBmn0iPQX4K3f/MfBkYIhgsc9u4P3AYyKK63qgQDDaMHQ2cM2CxZsAPwfOqlZyqP55VvW4iHSwdCq17JxxqKuMVxd6pmsLONPzEu5i2enK1PWM150bnspz6VW7573PxEyJjb2ZWm94ON4wqIzPJeOb6sZEbunvIunc/UKCMNSRcgAAIABJREFU3vBzgBcDzwGOcfe/ijWwCIQ7rqoyLiJJ0mwynmWuMv0sYBL4cd25QhRBuXsO+DTwETN7tJmdC7wR+ACAmW03s3A3ki8CA8AHzeyhwL8Agxw6AlFEOkw6xZJtKvWV8UK5UpumEu7A2ZVJkS/OH20Y9ufD/GT8v6/Zw1u+fNO8qSnjM0U29mYZ7A0r42EyXqol6BAs5AzXUtRPUkkydx9392+7++eABwjWBSV3K9YGzZbUMy4iydPsE+k64DwzewzwQuBb7l4ysy0E4w6vjjC2C4BfAN8HPgq83d0/Xz13P8Eccdx9guAfBmcB1wKPB87Rhj8inS+TStUWZoYOWcBZrYzXFnCm646X5/eMZ9N1CzhLc0n+ntFgCcpobq7eUEvGa5XxxdtUUiljYzVhDzf8SSoz22lmPzGzt1V/fhXBp4z/DdxhZmfGGuAqqTIuIknU7DSVNwGXAX9IUCF/e/X4rdU/nxFRXGF1/KXVr4XnbMHPvwAeEdV7i8jakE4ZCwrjdW0qcws1c4VSrZ0l7AnvzqQolp1KxUlVd+kc6M4s2jO+d3QGODQZP3Fbf22e+FxlvFjb8Ce0qa+LsVxxLVTGl1oX9BcERZGLCHZFXpNUGReRJGrqieTuVwI7CEZcHe/ut1VPvQR4SESTVEREGpJJWW1nzVBtmorNr4CXaj3jc20qMJd0F0rBNJXFesb3VJPxsdzcTpvj1QWcC3vGJxZUxmFuokrSK+O0b11QLFQZF5EkOpxNf2aoa0epzhkfIN4xzyKyDqUWnTMe/DnXM24USz432rCuYg7BuLuebLpumsr8ZNzd2Ts2vzJeqTgTs0GbSl9XmnTKmJgtUq44U/nSIZXxobBNJfmV8basC4pLXpVxEUmgpp5Ind5PKCJrS2aRHTjn5owHP3dVp6YsHG3YXd1wJ0zQwmkqYU95mLyP5orkCsE1YWV8Ml/CHTb2BvPEB3syTMyUmMoHfeODCyrj4a6ba2CaSjvXBbXd3GhDVcZFJDmSOmdcRGRF6ZThHlSqQ7UFnHWV8foFnLXRhtWkPJw1HkxTSdWmrYSV8XDxJsBYtTIetqSEk1QGe7NMzBaZrPaNDy6sjIfJeMLnjBOsC3oG8DOCTzvr1wU9GHhLTHFFItz0pyejyriIJEezbSpPAV7h7j82s2dQ109oZu8HvhR5hCIiSwhbTsrupKqdcgunqXRXe8bDSne2bpoK1CXj5QpdGSOVMjLVBZ0wt3gTgio5BP3iQG1KysbeLBMzRSarG/8s7Bk/9agNHDPUW2tXSSp3v9LMdgCnAb9096nqqZcAV7n7aHzRrZ4q4yKSRM0m4x3dTygia0s4MaVcccL8quKOGbXt6MMdOMMFnJmwTSUz1zMOQSU87CMPXhMk7+HizWOGemuV8YXJ+GBPlvF5yfj8pPv5u3bw/F07ovyrt0wnrwsKF3CqMi4iSZLkOeMiIssKK+OlBW0qYYsKBAs1C6VK7ZpMaonKeLVNBeZaWwD2js2woTvDcVv6lqyMD/ZmmJgt1dpXFlbG14pOXxc0WyqTTlntH2QiIknQ7BOpo/sJRWRtCccU1k9UqfhcxRwgm6lWxhdMUwnH24WjDYvlSm3Dn67qayDoGT9mU291Vvj8nvH6yvjETJHJ/NpOxunwdUH5YkVVcRFJnKZ+Y3R6P6HI/9/e/QdJdlWHHf+e/jEzu/q1q5UsIVYrRYbYFtgKeBV+CRfBQGwwrjiOSRy7YlMhcmwSCuNAiliCOJiYGMoVq2I7rA0BHLBNKlCOHKBciYEgkSLIICfYDsjEsCCBkNCvlXZnpn/c/PH6db/u7VntSN3z3uv3/VRt7UxP78y96tWZM2fPPVf10mlN94lD1qYyWxnvDdI4uW7PVMa3ekNSykYfTrep5Mn4KQ4f3McF+7vjaSrffCRLyvMpKefvy9pU/vLe7LDnbJtKjaz0uaDN/sB+cUmVs+sSQUrpVErpNuBQRDwzIs4BbjURl7TXWuM2lckFPYNhGifcMEm6T20P6LZj3Es+ufRnMK6O54/lPeMpJe68/xSHD+7n4P4uD5zqkVLiq/ef4qJz19i3liV252902OoPuem/38F5G51xxbyGVvpckJVxSVW0639LjYgfBt4CfCswJLuN840RcQJ4eUqpd6Y/L0mLMq8yPhgmCrn4uNp9sjeg05okYuuFnvHZSSvddrA9GPLQqT4ntvo88cA+IrLP/dBmP2tdObBv/Llees1lnNjs852HL+DZ33rROKmvofxc0CYreC5osz+0Mi6pcnZ76c/LgPcDHwdeVvjzHwR+CHjDQlcnSWfQ3qlNpdgzPkqwT271x8k7FNpU+kN6o8OaU20q/SFffSBrOzl8cN94VvgDJ7fH1fLcFYfO4fUv/g5+4Lsu48LqX+xzJit9LmirNxj/ECZJVbHbqPQGsv7BV5Al4ACklN4F3AD82OKWJklntlNlfLpNJauEntwejG/XhEnivdUfjttUZg9wjscaHpzMCL//ZI87H8j6yFdNSulTwOVk/+J5ZUrpz0cf+gfAt6eUPlva4hbAyrikKtptMv4k4EM7fOyzwBMe33Ik6ey154w2HKZEK06vjJ/qDaZG2q13C8n4qDLenZkznl/4c/jgfg6ekyXjd9x9gq3+cCWTcVjtc0FbvYE945IqZ7dR6Thw3Q4f++vAVx7fciTp7M1rU8l6xk9vRzm5PZhqU1lvj0Yb9ofjySnr4wOcWc/4V+8/xb5um4P7u+M2lT+96yEgq5avooj44Yi4A/hL4Bay9pT3RsR/jIjankwFK+OSqmm3yfi/A/5FRPwCcBRIwBMj4hXAzwPHFrw+SdrR/DYVpttURtXuR7b6020qhQOc4zaVmdGGdz5wksMH9xER4zGG/+fOBwGmesZXxaqfC7IyLqmKdjtn/KaIOEh2qv4GsgM+vw/0gJtSSm9b/BIlab780p/8Qh8YtakU8q3xaMPegG7hA5MDnAN6/XyayiQZ7w/SeMY4ZOMLAf4sr4wfWMnKeH4u6DURMS4hp5TeNYr9/xS4sbTVPU5b/SEbVsYlVcyuRxumlH5hdPnDM4FDwIPAp1JK9y56cZJ0Jp2d5oxP9YxP2lQO7p90WbRbQacVU5XxtUKbSm8w5OsPbfK0Iweyr9Vucf5Gdu39wf1dzlmv7S2bZ/Ik4NU7fKz254KcpiKpis46KkXE5aPbN0kpPQh8gSwhfyVwY0R823KWKEnz5W0nU20qKY0vA4LpS3/aremQt9ZpZcn4+ABnPme8xX2PbPPAyd5UO8rB0djCVWxRGVnpc0GbVsYlVdCjlnYi4gLgd4C/OXr/Q8BrgFuBg2S3tX0f8IqIeG5K6TPLW64kTeS94b1im8qOlfH+ONnOrXda2ZzxwfSc8bV2i2+c2AKm21EO7F/jy988ubKTVMjOBb0tsmtKP8TkXNDTyc4FvanMxT1eVsYlVdHZRKW3ANcALwf+NvAtZIn4ncAVKaVLgKvIKuW17SWUVD95ol1sU5m99CdPvoaJqWkqMKmM905rU5mExmLinc8aX9F+cVJKNwFvBl4LfJLJuaBfA47V/VyQlXFJVXQ2TY8vBW5MKb0HICK+CPxv4J+klO4CSCl9KSLehNNUJO2hzpw544MhM3PGW4Xnz2lTGcyZM96Z/PmpNpVRz/kKV8ZX9lxQbzBkMExWxiVVztlEpUuAOwrv529/eeZ5dwIHFrEoSTobnR2mqbTn9IwDU6MNAdY7bbb6g7mjDbOPt7jo3Mn19vms8VXsGV/1c0Fbox+4rIxLqpqzScbbwFbh/f7o996c58acxyRpKfLkuj+YnqZSPMBZ7BMv3sAJWW948QBnsWccsot9olBlz2eNr9KFPxFxwegs0JeAL0XEzRHxZOBTwE8D30020vAzo97xWtrsDYDJzauSVBVGJUm1lSfaxTaVYUoUC+DFynh3Ts/4Vn/IPQ9n9YYLR1XwvDI+2xt+9WXnc/F561xxaKUq4404FzSujHesjEuqlrMdlPtzEXH36O38u9lrI+KewnMuWdyyJOnRjS/9mZ0zPucGzuz586epHP/mSQ6ds8a5o9nheTI+247ywqsv4YVXr1yoa8S5ICvjkqrqbJLx42TzZYu+TNZLOO+5krQnOnNGGw6GaeoA51RlfLZNpdPixGaf4/ed5Eih2p0f4Fzlg5oFjTgXtNXLfmBbtzIuqWIeNRlPKV25B+uQpF3Lk+vBTJtKcWrK1DSVOXPG7+0PuefEFtdeeXDyZ1p5ZbwRyXgjzgVt9q2MS6omo5Kk2srbTmYPcBbbUTqtIC+Uz442XO+0eWSrz9cePMWRCwuV8XajKuONkFfG7RmXVDVn2zMuSZWTJ81TbSqJqWkqEUF3NDVl3qU/X73/JMMERw6dM3780LnrdNvBFYXHVtzKnwuyMi6pqkzGJdVWZ16bynB6mgrAep6Mz2lTyf9ocULK9z/1Uv7a5c/jonPXl7PwamnEuSAr45KqymRcUm2ND3CeYZoKQLfTgq35BzhzVxTaVDrtFpdfuFLjC3fUlHNBW1bGJVWUUUlSbXXGPePTBziL01RgMt5wNknPH9/otrj4vEZUwRtrXBn3Bk5JFWMyLqm25h3gHKZ5lfHs/dPaVEZV0iMX7p+6aVOrZ9wz3vHbnqRqMSpJqq2IoNOKqRs4B8M0dYATJhXw7sw0lbV2ViU9cmFjDmo2lpVxSVVlMi6p1jrt6WR8mKA9U+XOe8VnK+N5z/iKXW+vOcY3cFoZl1QxRiVJtdZtteidYc44TBKw2dGG6ybjjbHVH9JuxWmHeCWpbEYlSbXWbsfUaMPB8PQDnJPK+PxpKk2ZnNJkm72BVXFJlWRkklRrnVZr6tKf7ADn9HPWdqiMH9jfBeBJF5+73EWqdFv9of3ikirJZFxSrXXbMTVNZe6c8fwA50yW/qKrL+X3X/kcK+MNYGVcUlUZmSTVWrs13aYyTOm0MYV5Zfy0OeOdFtdcfmD5i1TpNq2MS6ook3FJtdZtt+jN9IzPTlMZjzZsO0u8qbasjEuqKCOTpFrrtM6mTWV06U/LkNdUm/0h61bGJVWQ35kk1Vq7dfqc8dlpKuMDnFbGG8vKuKSqMjJJqrVuuzVVGZ83TWU82tDKeGPZMy6pqvzOJKnWZm/gHAwTrTkHNfPnqloi4oaI+NKyv46VcUlVZWSSVGtZz/jMnHEPcNZCRHwHcMNefC3njEuqKpNxSbXWabXoD898gHMy2tCQVxUR0QLeAXx6L76elXFJVWVkklRrnXaMb+BMKc09wDm+9KdlZbxCXgWcBN69F18s6xn3W56k6jEySaq1TuHSn7x1fKfKeGf2ZKdKERFXAT8P/NRefc2sMm6biqTq6ZS9AEl6PDrtFr3RNJU8KT99zrgHOPdSRGwAh3f48N3AbwK/nFL6YkT8jUf5XNcD1wMcOXLkMa9psz+0TUVSJZmMS6q1bmGayjBlv+84Z9w2lb1yFPjEDh+7HjgA/MrZfKKU0jHgGMDRo0fTozx9rsEwMRgmK+OSKslkXFKttVutcUV8Uhmffs6aN3DuqZTSLcDcn3wi4qPAU4EHI/uhqQOsRcTDwNUppeOLXs92P/uXkzUr45IqyGRcUq11WzFpU3mUyrijDSvhx4F9hfd/BHgl8DzgrmV8QZNxSVVmMi6p1jrtyZzx4XB+Mv6sqy7iZUcPc8Whc/Z8fZqWUrqz+H5E3AP0U0p/sayvuT36YW3NH8YkVZDJuKRaa7da457xnQ5wXnrBBr/8d67Z87WpGsbJuJVxSRVkZJJUa9kBzizZykcbtjyoWRsppd9KKV25zK9hm4qkKjMySaq1Tqs1aVMZ9Yy3w2RcE+NkvO00FUnVYzIuqdY6hcr4TtNU1Gx5Mu4BXklV5LcsSbXWaU0OcA52OMCpZrNnXFKVGZkk1VqnnR3gTClN2lTsGVeBPeOSqszIJKnW8ls1+6NbFsFkXNPyyvi6ybikCjIySaq1zqgPeDCcVMZtU1GRBzglVZnJuKRa646uuO8NhowKoFbGNSW/obXb8e+FpOoxGZdUa3llvD9IHuDUXJPKuN/yJFWPkUlSrRV7xj3AqXk8wCmpyoxMkmqtM6p29odD54xrri1HG0qqMCOTpFobV8YHicGoMh62qaigZ5uKpAozMkmqtXHP+GjWOEDbZFwFXvojqcoqGZki8+aI+EZE3B8Rb4uIHWdSRcT3RsSnI+LhiPh8RPzDvVyvpPJ0RtNU+k5T0Q48wCmpyjplL2AHPwv8JPAjQADvBe4F3jL7xIh4MvAHwJuAHwWeAbwjIr6RUrp5rxYsqRzdUWW85zQV7WC7P6QVk/MFklQlVY1MrwbemFL6eErpY8A/B165w3P/LnB7Sulfp5T+IqX0XuA9wI/tzVIllak9qowPnKaiHWwPhnRNxCVVVOUq4xFxGXA58D8KD98CHI6Iy1NKX5n5I+8HPjzzWAI2lrdKSVWR94z3nKaiHWz3h/aLS6qsKkanJ4x+v6vw2N2j3w/PPjml9IWU0h/n70fEJcDfYzqZl7SiuuOe8ck0FdtUVLQ9GLJuMi6pokqpjEfEBnMS65H9o9+3Co/lb68/yuc9B/gAWSL/73d4zvXA9QBHjhw5yxVLqqr2+NKfIcOhbSo63XZ/6OFNSZVVVpvKUeATO3zsdaPf14Fe4W2Akzt9woi4gOwg51XAdSmluc9NKR0DjgEcPXo07W7ZkqomP8DZ9wCndmCbiqQqKyU6pZRuSSnFvF9kk1MALi38kfztr837fBFxEfBRskT8eSmlLy5t8ZIqpXgDZ3+UjHtYT0U9D3BKqrDKRaeU0l3AceC6wsPXAXfNObxJRKyRVcQvAr4npfT5PVmopEoo3sB5ansAwL7ujtcSqIGsjEuqsspNUxn5DeCXIuI4MAB+CfjV/IMRcTFwKqX0MNlM8u8Gvg94JCLyKvp2Sum+vV22pL1WvIHzVC9LxjfWTLw0sT0wGZdUXVVNxt8KXAz8Z7Jk/D8Abyt8/NPAu4B/SXYxUAf4bzOf41amq+uSVlB+A2dvMGRzlIzvX6tqaFMZtjzAKanCKvkdK6U0AH5u9Gvex68svH10j5YlqYLyNpXBcNKmsmEVVAXb/SHnbVTy250kVa9nXJJ2o1OYpnKyN2Ct3fLac03pDayMS6ouo5OkWsunZPSGQ05tD9joGtY0zQOckqrM6CSp1tqFNpXN3oB9a05S0TQPcEqqMqOTpFrrjg9wZtNUPLypWd7AKanKjE6Sam3SMz7k5PaADWeMa0ZvMKRrZVxSRRmdJNVa3qbSz9tU7BnXDEcbSqoyo5OkWssPcOY3cNozrlnb/SHrVsYlVZTRSVKttVtBBPSHQ071Buzr2jOuiZSSBzglVZrRSVLtdVqRHeC0Mq4Zg2Eipcm/oEhS1RidJNVep9ViMK6MG9Y0sT0YAlgZl1RZRidJtddpx3i04T6nqahguz9Kxq2MS6ooo5Ok2uu0IusZ3x6wYZuKCsbJuJVxSRVldJJUe512i14/sdUfst8DnCrYMhmXVHFGJ0m1120FD2/1Adi3ZljTRG9gm4qkajM6Saq9djt4aLMHYM+4pniAU1LVGZ0k1V631eLEZlYZ3zAZV4EHOCVVndFJUu112sGJvDLuAU4VeIBTUtUZnSTVXrtQGd9vMq6CvE3FS38kVZXRSVLtddthm4rmsjIuqeqMTpJqr9MKTvUGgAc4NS1PxtdNxiVVlNFJUu11WpNQZs+4ipymIqnqjE6Saq/TjvHbXvqjorwybs+4pKoyOkmqvU4h0drw0h8V9KyMS6o4o5Ok2uu0JpVxe8ZV5JxxSVVndJJUe8Vk3GkqKtpymoqkijM6Saq9vB+42w57gzUlP8DpNBVJVWV0klR77VFl3BYVzer1E+ABTknVZXSSVHv5NBXHGmrW9mBAuxXjH9gkqWpMxiXVXnc0Z9zKuGZt94ce3pRUaUYoSbXXHlXGPbypWdv9oYc3JVWaEUpS7XVHLQj7bVPRjO1Bsl9cUqUZoSTVXn7pjz3jmrXdHzpJRVKlGaEk1V7HaSrawfbANhVJ1WaEklR7HXvGtYPt/sADnJIqzQglqfY6TlPRDjzAKanqjFCSaq/b9gCn5usN0vjvhyRVkcm4pNprjyrjGybjtRARz46INPPr9mV8LSvjkqquU/YCJOnxyiuftqnUxtXAbcBLC4/1lvGFtgZDLljrLuNTS9JCmIxLqj2nqdTO1cCfppS+vuwv5A2ckqrOCCWp9trOGa+bq4HP78UX6g2GrHXsGZdUXSbjkmqva2W8bq4Gro2Iz0XE8Yh4e0RcsIwvZGVcUtXZpiKp9ryBs1oiYgM4vMOH7wYuBz4D/ARwCPgV4H3AS+Z8ruuB6wGOHDmy67V4gFNS1ZmMS6o9e8Yr5yjwiR0+9nLgIHAipTQAiIifAG6LiCMppePFJ6eUjgHHAI4ePZp2uxBv4JRUdSbjkmqv4zSVSkkp3QLsplH7z0e/PxE4fqYn7tZ2f0jXNhVJFWaEklR74xs4bVOpvIi4NiJORMRlhYefBgyAOxb99ayMS6o6I5Sk2nv6kQO85LuewF+95Lyyl6JH9yfAncA7IuIpEfE9wG8B70wp3bvoL/aUy87n8MH9i/60krQwtqlIqr1vOX+DX/v7Ty97GToLKaXtiHgx8G+BW4E+8F7gdcv4eh/8mecs49NK0sKYjEuS9lRK6f8BP1j2OiSpCmxTkSRJkkpiMi5JkiSVxGRckiRJKonJuCRJklQSk3FJkiSpJCbjkiRJUklMxiVJkqSSmIxLkiRJJTEZlyRJkkpiMi5JkiSVxGRckiRJKonJuCRJklQSk3FJkiSpJCbjkiRJUklMxiVJkqSSREqp7DWUJiLuAb5c9joW5CLg3rIXsYfc72pr0n4f616vSCldvOjFVNmKxWzw7/kqc7+ra+Exu9HJ+CqJiNtSSkfLXsdecb+rrUn7bdJeNa1Jr32T9grud5UtY6+2qUiSJEklMRmXJEmSSmIyvjqOlb2APeZ+V1uT9tukvWpak177Ju0V3O8qW/he7RmXJEmSSmJlXJIkSSqJybgkSZJUEpPxBoiI9Yh4X0T8z4j4w4j4K2Wvadki4g0RcWtE3BYRf6vs9SxTRKxFxO9FxCci4iMRcVHZa1q0iGhFxNsj4pMR8dGIuKrsNS1TE15TnVnT4rYxe7UYs3f3mpqMN8M/Ah5IKT0LeDVwU8nrWaqIeD7wnSml5wAvBp5U8pKW7UeBr6SUngv8LvD6ktezDD8EtFNKzwZuBN5a8nqWrQmvqc6sMXHbmL2S/38bs3ehs5QlqWquBj4MkFL6s4h4SsnrWbYXAl+IiJuBdeBVJa9nqVJK746I/P/lw8DdZa5nSZ4DfAQgpXRLRLyv5PUsVUNeU51Zk+K2MXv1GLN3wWS8GW4HXhoRfwA8k+wvyiq7GHgi8IPA04F3As8udUVLllLqR8R/Ba4FXlT2epbgfODBwvtR1kL2SgNeU51Zk+K2MXv1GLN3wTaVFRQRvxgRH4uID48eeifwCPBR4AeA20pb3BLM2e83gY+klHoppU8Bl5e4vIWbs18AUkovIatG/KdyVrZUDwHnFd4flLWQvbTir6kKmhS3jdmZFf//25i9C1bGV1BK6YaZh54BfCyl9LMR8QxWLNDN2e+twD8GfjUivo0V+yfA2f1GxM+MHv914GFWM+jdCnw/8IGIuA74bMnrWaqGvKYqaFLcNmY34v9vY/YuWBmvgdGp+s9FxAtmHjsWEfdHxNcj4nVn+BRfAF4VEZ8E3gic6bmlW8B+bwb+b0T8L+C3gZ9e9pofjwXs93eAF0fEx4H3kx38qrTHsOcPAsPR3+F/A7x2r9f8eDyG/dbuNdW0JsVtY7Yx25j9+F5TK+MVFxEbwPuA2cM7byXrqXsBWS/hb0fE8ZTS785+jpTSPcD3Lnuti7Cg/SbgNcte6yIsaL/3k/0zdi08lj2nlIbAT+3tShfjMe63Vq+ppjUpbhuzx4zZxuzH/pqmlPxV0V9kp+lvB/4ESMALRo+fA5zK3x89dgNwS9lrdr/u1z03d7/+atZr3qS9NnG/TdxzWfu1TaXangv8IfCsmcevIRv/dEvhsVuAa2MyWqeO3G9mVfcLzdtz0/arZr3mTdorNG+/0Lw9l7LfOv8HW3kppbfnb0dMTQV6AnBfSmmz8NjdwBrZiKiv7ckCF8z9jq3kfqF5e27aftWs17xJe4Xm7Reat+ey9mtlvJ72A1szj+Xvr+/xWvaC+13t/ULz9ty0/apZr3mT9grN2y80b89L3a/JeD1tcvqLn79/co/Xshfc72rvF5q356btV816zZu0V2jefqF5e17qfk3G6+lO4GBErBUeu5Tsp7T7ylnSUrnf1d4vNG/PTduvmvWaN2mv0Lz9QvP2vNT9mozX0+3ANtPXBV8H/HFKqV/OkpbK/a72fqF5e27aftWs17xJe4Xm7Reat+el7tcDnDWUUjoZEe8Gfj0ifpLsp7N/Blxf6sKWxP2u9n6heXtu2n7VrNe8SXuF5u0XmrfnZe/XZLy+XgP8BvBHwEPAv0op/V65S1oq97va+4Xm7blp+1WzXvMm7RWat19o3p6Xtt8YDS6XJEmStMfsGZckSZJKYjIuSZIklcRkXJIkSSqJybgkSZJUEpNxSZIkqSQm45IkSVJJTMYlSZKkkpiMS5IkSSXxBk5pgSLiBuDAzMMPpJR+sYz1SJJ2ZsxWFXgDpyRJklQS21QkSZKkktimIi1IRATweuCpwJuBpwGXAAdSSjeWuTZJ0jRjtqrCyri0OC8C/gtwHDgG3AzsB55f5qIkSXMZs1UJJuPS4hxKKX0OeCbwmymlB4H3AD9e7rIkSXMYs1UJHuCUFigi9gEPAE9OKR0vez2SpJ0Zs1UFVsalxXou8BWDuiTVgjFbpTMZlxbr+cAflb0ISdJZMWardCbj0mJ9O/CBshchSTorxmyVzp5My+LoAAAAU0lEQVRxSZIkqSRWxiVJkqSSmIxLkiRJJTEZlyRJkkpiMi5JkiSVxGRckiRJKonJuCRJklQSk3FJkiSpJCbjkiRJUklMxiVJkqSSmIxLkiRJJfn/YItwLdjXvPUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x533.333 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "\n",
    "rcParams['font.family']='sans'\n",
    "rcParams['font.sans-serif']='Arial'\n",
    "rcParams['font.size']=14\n",
    "rcParams['mathtext.fontset']='cm'\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "n_values = np.arange(-3.,11.,0.1)\n",
    "besselj_values = np.zeros(len(n_values))\n",
    "bessely_values = np.zeros(len(n_values))\n",
    "\n",
    "fig = plt.figure(figsize=(12,12/1.62))\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "#ax3 = fig.add_subplot(223)\n",
    "\n",
    "for idx, n in enumerate(n_values):\n",
    "    besselj_values[idx] = sp.N(sp.besselj(1,10**(-n)))\n",
    "    bessely_values[idx] = sp.N(sp.bessely(1,10**(-n)))\n",
    "\n",
    "                           \n",
    "ax1.semilogx(10**(-n_values),besselj_values)\n",
    "ax1.set_ylabel(r'Bessel fn. of 1st kind: $J$',fontsize = 16)\n",
    "ax1.set_xlabel(r'$\\bar{r}$',fontsize = 16)\n",
    "\n",
    "\n",
    "ax2.semilogx(10**(-n_values),bessely_values)\n",
    "ax2.set_ylabel(r'Bessel fn. of 2nd kind: $Y$',fontsize = 16)\n",
    "ax2.set_xlabel(r'$\\bar{r}$',fontsize = 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that as $\\bar{r} \\to 0$, the Bessel function of the second kind, $Y \\to -\\infty$.\n",
    "\n",
    "Now, we know that the solution for $W$ should be such that it has zero value at $\\bar{r}=0$ (from symmetry conditions). Therefore, the constant $C_2$ must be zero so that the solution reduces to:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle W(r) = C_{1} J_{1}\\left(\\bar{K} \\bar{r}\\right)$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "C1, Kbar, rbar = sp.symbols(r'C_1, \\bar{K}, \\bar{r}')\n",
    "W_soln = C1*sp.besselj(1,Kbar*rbar)\n",
    "display(Math(r'W(r) = {}'.format(sp.latex(W_soln))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we use the condition that at $\\bar{r}=1$, we must have $W = 0$, i.e. the slope is zero at the periphery of the circle, so we need to find values of $\\bar{K}$ which will satisfy $J_1 (\\bar{K})=0$. These will be the values which will lead to the critical values of $P$. \n",
    "\n",
    "Finding out the values of $K$ is a little tricky. We use `SciPy` to evaluate the numerical values of the Bessel function of the first kind as well as to find the desired values. Note that the root is very sensitive to the choice of the trial value. To obtain some idea of the trial value, we first plot the Bessel function of the first kind. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEUCAYAAAAbV1CxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1f3/8dcnOyRhD9nYdwiQCAiKoIC4oIJate5Vv661aq3VWlu1amtr1VZr677UrYprRakrAgICCpF9CTuEJQkhLNm3+fz+mOAvRpYsM3Pnznyej8c8IHfuzH1fI/OZc+6554iqYowxxgRKhNMBjDHGhBcrPMYYYwLKCo8xxpiAssJjjDEmoKzwGGOMCagopwO4QadOnbRHjx7Nem1paSnx8fG+DRTk7JzDg51zeGjJOWdnZxeqalLD7VZ4GqFHjx4sXry4Wa+dPXs248aN822gIGfnHB7snMNDS85ZRLYearvru9pEJFZEnhORvSKSJyK/OcK+A0RkpoiUicg6ETkvkFmNMcaERovnEWA0MBHoArwmIttUdWr9nUQkAZgBzAQygUnAmyKSpaqrA5zZGGPClqsLj4jEA9cCk1U1G8gWkYeBm4CpDXb/GVANXK2q1cB6ETkVOB6wwmOMMQHi6sKDt+USC8yrt20ecI+IRKlqTb3tE4AP64oOAKp6VmBiGmOMOUjcPFdb3TWaZ1W1U71tA/G2YNJUdVe97UuAd4F04FxgF3Cvqk4/zHtfB1wHkJycPHzq1IYNqMYpKSkhISGhWa91Kzvn8GDnHB5acs7jx4/PVtURDbe7vcXTGqhssO3gz7ENticCdwBPAWcApwIfiMioum66H1DV54DnAEaMGKHNHdVho2DCg51zeLBz9g23F54KflxgDv5c1mB7DbBCVX9X9/MSERmLt1Vzvf8iGmOMqc/thWcH0F5EYlS1qm5bCt5WT1GDfXcCGxtsywEG+TdiaCmtrGFp7j42FJRQXFFNRISQ3q4VGWlt6Z0Uj4g4HdEYE+TcXniWAlV4h1PPrts2BshuMLAAYAFweoNtg4AtfswXElSV+Rv38O+vtzBn/W6qajyH3C+9XSvOH96Fnp5DP2+MMeDywqOqZSLyCvCUiFyJt7VzO3WDAkQkBdivquXAs8AtIvJXvNduJuO992ekE9ndYn1+MfdMW8nCTUV0SojhslHdGdc/iQEpibRrHUONx8P2veVkb93LJyvzeGLmeqIjYJPkcOP4PsRFRzp9CsaYIOPqwlPnNuBpvDeGHgAeUNW36p7bBVwFvKyq20TkFOAJ4JfAJuA8VV3iQOagp6o8P3cTD3+aQ3xsFH88O4MLRnT9USGJIYJ+yYn0S07k4pHd2FxYyp2vz+WJmRv4ZGUej1+URUZaW4fOwhgTjFxfeFS1DLii7tHwOWnw80KshXNU5VW13Pb2Uj5ZmcekwSn86ZzBdExoOIbj0Hp2iufGrDh+kZbBHe8s49yn5vO3CzKZnJnm59TGGLdw/VxtxrcOVFTzs5e+4dNVefzujAE8demwRhed+k7ql8Snt55IZpe23PzmEp75quG4DmNMuLLCY75XUlnD5S98w5Jt+3jiomO47sTeLRql1iE+htevGcXkzDQe+mQtT87a4MO0xhi3cn1Xm/GNqhoPN7yWzcqdB3jmsuGcMijZJ+8bGxXJ4xdmESHwyGc5xEZFcM3YXj55b2OMO1nhMagqv31/OfM2FPLI+UN9VnQOiowQ/nZBJpXVHh78eA3dOrTm1IwUnx7DGOMe1tVmeH3hVt7/bge/mtiPC0Z09csxoiIjeOzCLIamt+WXU5eyaud+vxzHGBP8rPCEuSXb9vLA9NWcPKAzN0/o49djtYqJ5PmfjaBd62iuezWb/eXVR3+RMSbkWOEJYwcqqrnpjSUkt4nj7z/NIiLC/9PddG4Tx9OXDSf/QAW/fW85bp4d3RjTPFZ4wtgfP1rNrv3l/OuSYbRtHR2w42Z1bccdp/Xnk5V5vPHttoAd1xgTHKzwhKmZa/N5J3s7Px/Xm6yu7QJ+/GvH9uLEfkk88NFqthSWBvz4xhjnWOEJQwcqqvnteyvon5zILSf3dSRDRITwyPlDiYmK4M73luPxWJebMeHCCk8YevyL9ewuqeTh84cSG+XcJJ7JbeK458xBfLO5yLrcjAkjVnjCTE5eMa8s2MLFI7uR6UAXW0MXjOjC2L6deOiTteQfqHA6jjEmAKzwhBFV5Q8friQxLoo7Tu3vdBwARIQHzxlCVa2Hhz5Z63QcY0wAWOEJI9OX72LhpiJuP7U/7eNjnI7zvW4dW3Pd2F78d8kOFm9puHCsMSbUWOEJE1U1Hh75LIcBKd51c4LNjeN7k9o2jj98uIpaG2hgTEizwhMm3lqcy7aiMu48fQCRAbhRtKlax0TxuzMGsmrnAd5enOt0HGOMH1nhCQNlVTU88eV6RvbowLj+SU7HOayzhqYyont7HvtiHeVVtU7HMcb4iRWeMPDy/C3sLq7kN6f3b9H6Ov4mItw5aQAFxZW8PH+L03GMMX5ihSfEFVdU88zsjZw8oDMjenRwOs5RHdujAxMGdObp2RvYX2aTiBoTiqzwhLjXF27jQEUNt07s53SURrvjtP4UV9bwtC2XbUxIssITwiqqa3lx3iZO7JfEkC5tnY7TaANT23BOVjovz9/M7uJKp+MYY3zMCk8Ie2tRLoUlVfxiXG+nozTZzRP6UFXj4YV5m5yOYozxMSs8IaqqxsOzX23k2B7tGdWro9NxmqxXUgJnDU3jtQVb2Vta5XQcY4wPWeEJUdOW7mDn/gpuHO/fVUX96aYJfSirquWlrzc7HcUY40NWeEKQqvLcnE0MSm3DuH7Be9/O0fRLTmTS4BRe/nqLLZNtTAixwhOC5m0oZH1BCVeP6RnU9+00xk0T+lBcWcMrdl+PMSHDCk8I+vfXW+iUEMtZmalOR2mxjLS2TBjQmZfnb6Gi2mYzMCYUWOEJMZsLS5m5toBLR3VzdJE3X7p2bC+KSqt4/7sdTkcxxviAFZ4Q88r8LURHCpceF3wzUDfXcb06MDi9DS/M22RLZBsTAqzwhJADFdW8sziXyUPT6JwY53QcnxERrh3bi027S5mVU+B0HGNMC7m+8IhIrIg8JyJ7RSRPRH7TiNd0EJF8EbkyABED5t3F2ymtquWqE3o6HcXnzhiSSmrbOJ6fazeUGuN2ri88wCPAaGAicD1wt4hcdJTXPA509newQFJV3vx2G1ld27lqepzGio6M4KoTerBwUxErd+x3Oo4xpgVcXXhEJB64FrhVVbNVdRrwMHDTEV4zCRgJ7A5MysDI3rqX9QUlXDIqdK7tNHThsd2Ij4nkBWv1GONqri48QCYQC8yrt20ecKyIRDXcWUQSgWeA64CQmofljW+2kRgbxVlD3T+E+nDatormghFd+d+KXTZ5qDEu5vbCkwoUqWpFvW35QAxwqFv2HwY+VdU5gQgXKPvKqpi+YhfnHJNO65gf1duQcvnx3amuVVse2xgXc/unVGug4Vffgz/H1t8oIicBk4GMxryxiFyHt2VEcnIys2fPblbAkpKSZr+2sb7YUk1VjYe+EfnMnl3o12M1hr/PeVDHCF78ah0DNJfIiOCYmSEQv+dgY+ccHvxxzm4vPBU0KDD1fi47uEFEWgEvADeraqOuTKvqc8BzACNGjNBx48Y1K+Ds2bNp7msbQ1X585I5ZHZtzc+mnOC34zSFv8+5olMeN7yeTW3yQE7OSPHbcZrC3+ccjOycw4M/ztntXW07gPYiElNvWwreVk9RvW0jgT7AayJSIiIlQBrwjIg8E7C0fvDdtr2syy/hkpFdnY4SMBMHdia1bRyvLdzqdBRjTDO4vfAsxTtIYHS9bWOAbFWtqbftW6AvkFXvkQ/cW/dwrbcW5RIfE8lZQ9OcjhIwUZERXDKyG3PXF7Jpd4nTcYwxTeTqwqOqZcArwFMiMlJEpgC3A08AiEiKiLRS1XJV3VD/AdQCBarq2lvhy6tq+XhFHmcMSSU+1u29pk1z4ciuREcKry/c5nQUY0wTubrw1LkNWATMxDtU+gFVfavuuV3AhU4F87fPV+dRUlnDT4Z1cTpKwHVOjOP0wam8k51LeZXNWm2Mm7i+8KhqmapeoaoJqpqmqn+r95yo6suHeV2Xwz3nFu9/t4P0dq0Y1bOD01EccdmobhRX1PDxil1ORzHGNIHrC0+4yj9Qwdz1u/nJsHQigmRIcaCN7NmBnp3ieWuR3dNjjJtY4XGpaUt34FE495h0p6M4RkS48NiufLuliI02yMAY17DC40KqynvZOzimWzt6JSU4HcdRPxmWTlSE8La1eoxxDSs8LrR61wFy8ovDclBBQ50T45gwoDPvfbed6lqP03GMMY1ghceF3sveQUxkBJNDeELQprhoZFcKS6r4co1rR8YbE1as8LhMTa2HD5ftZMKAzrRrHXP0F4SBE/smkdImjrcW2T09xriBFR6X+WZzEYUllZydFT4zFRxNVGQEF4zowlfrdrNzX7nTcYwxR2GFx2U+WraT+JhIxg8IqQVUW+ynI7riUXg3e7vTUYwxR2GFx0Wqajx8uiqPUwYlExcd6XScoNK1Q2tO6NORd7Jz8XjU6TjGmCOwwuMiX28oZF9ZdVhNCNoU5w3rQm5ROYu37nU6ijHmCKzwuMhHy3fSJi6Ksf06OR0lKJ0+OIXWMZG8/511txkTzKzwuERFdS2fr8rn9MEpxEZZN9uhtI6JYtLgVP63fBcV1TZxqDHBygqPS8zO2U1JZY11sx3FecPSKa6s4YvV+U5HMcYchhUel5i+fCcd4mMY3buj01GC2nG9OpLWNo73rLvNmKBlhccFyqpq+HJNAZMGpxAVab+yI4mIEM4dls6cdbspKK5wOo4x5hDsU8wFZqwpoLy6lsmZ1s3WGOce0wWPwodLdzodxRhzCFZ4XOCTFbtISozl2B7hueBbU/XpnEBm13a8990Op6MYYw7BCk+QK6+qZXbObk7LSCYyTBd8a47zhqWzZtcBVu884HQUY0wDVniC3Jz1uymvruX0DJuJuikmD00jOlLsnh5jglCU0wHMkX26Mo92raMZ1cu62ZqifXwMEwZ05oOlO/ntpAE2KMMElKqyvqCEbzbtYX1BCblFZZRV1aIKbVtHk9o2joy0Ngzv3oHeSfGIhFdvhhWeIFZV42HGmnxOy0gh2j44m+zcY7rw2ap85m/cw4n9kpyOY8LAtj1lvJOdy3vZ29m53zuqMjE2im4dW5MQG0VEBOQWlbFg4x5eXbAVgB4dWzMlK51LRnYjpW2ck/EDxgpPEFuwaQ/FFTWcnpHidBRXGtc/icS4KKYt3WmFx/jV1j2l/OPL9XywxDugZWzfJH45sS+je3eiS/tWP2rReDzKlj2lfL1xD5+vyuOfM9fz5KwNTMlM47ZT+tG1Q2snTiNgrPAEsU9X5tE6JpIxfW1utuaIi45k0uAUPl6Rx4PVg21Gb+NzlTW1PDlzA09/tZEIEa4e05P/G9OT1Latjvi6iAihV1ICvZISuPy47mzbU8arC7bw2sKt/G/5Lq46oQe3TuxHq5jQ/H/W+m+CVK1H+WJ1HuMHdLYPzBY4OyudksoaZq61ZbGNb63auZ8zn5jHEzM3cNbQNOb+Zjy/P3PQUYvOoXTr2Jq7zxrE7DvGcXZWGs/O2cQZT8xl8ZYiPyR3nhWeIJW9dS+FJVXWzdZCx/XqSFJiLNOW2j09xnfeXpzLT56aT3FFNS9fdSyPXZhF5zYtvz6T2rYVj1yQyRvXjKK61sNPn13A07M3ohpaa0xZ4QlSn67MIyYqwlYabaHICGHy0DRmrd3N/vJqp+MYl/Oo8sBHq/nNu8sZ3r09/7tlLOP6+/7f6Og+nfj01hOZNCSVv366lhtez6a0ssbnx3GKFZ4gpKp8tiqPE/t2IiHWLsO11NlZaVTVevhsZZ7TUYyLVdV4eG55JS99vZkrR/fgtatH0Skh1m/HS4iN4l8XH8PdZw7ki9X5XPL8QvaUVPrteIFkhScIrdixnx37yjnNutl8YmiXtvTo2Jppy6y7zTRPRXUt17y6mIW7arnz9AH8YfKggMwkIiJcM7YXz14+grV5xZz/zAJyi8r8flx/s8IThD5dmUdkhDBxYLLTUUKCiDAlM435G/dQcMBmrDZNU13r4aY3vmPOut1clRHDz8f1DvgNn6cMSuaNa0dRVFrFRc8tdH3xscIThD5blceonh1oHx/jdJSQMSUrDVX4aPkup6MYF6n1KL96aykz1hTwx7MzOKlrtGNZhnfvwH+uGUVxRTWXvLCQnfvKHcvSUi0uPCLSR0Qmish5InKxiEwRkWEikuiLgOFmc2EpG3eXcuoga+34Up/OiWSkteHDZbZUgmm8P05fzfTlu7hr0gAuP76H03EYnN6W164exb7Sai55fiGFLr3m0+TCIyJtROQ6EflERIqBHOBz4B3gP8AHwGJgr4gsE5E/icgAn6b+YZ5YEXlORPaKSJ6I/OYI+14oIitFpLQu22R/5WquL9d4l2w+2brZfO7srDSW5e5jS2Gp01GMC7y2YAsvz9/C/53Qk+tP6u10nO9ldm3Hy/83krwDFVz9ymLKqtw32q3RhUdE4kTkfmATcDWwCrgMGAZ0BxKBWCANGAycCrwNjAC+EZGPRKSfb+MD8AgwGpgIXA/cLSIXHSL/WOA14B9AJvAi8L6IHOOHTM32xep8BqQkhvyUGU6YnJmGCNbqMUc1Z91u7vtoNScP6MzvzxzodJwfGd69Pf+46BiWb9/HLW8updbjrvt8GlV4RCQLWAAkASNVdZSq3q6q01R1marmqmqpqlarap6qrlbVmar6oKqeDqQDc4HPReRGX4UXkXjgWuBWVc1W1WnAw8BNh9j9CuA9VX1eVTeo6hPALOBCX+Vpqb2lVSzeutcGFfhJattWjOzRgQ+W7gi5G/KM7+QWlXHzm0vo2zmBJy4+JmjXwTotI4X7JmcwY00+f5y+2uk4TXLUwiMixwOPAWer6o2quqmpB1HVElV9GBgAZIrIn5se9ZAy8bay5tXbNg84VkQa3gDzT+CPDaMBQTMd7Ox1BdR6lIl2fcdvzs5KZ9PuUlbZAnHmEKpqvCPYPKo8e/lw4oP8ProrRvfg6jE9eXn+Ft7Nds/aU41p8ZwBnK6q21p6MFWtUNXrgY0i4ov2aypQpKr1x8jmAzF4W2f1j71MVb//WiAiGcDJwBwf5PCJGasLSEqMZWh6W6ejhKwzhqQQHSk2hY45pD9/vIZl2/fzyPmZdO8Y73ScRrlr0gCO79WR3/13BSt37Hc6TqOIm7scRORy4CFVTa+3rRewEeipqlsO87rOeFtGO4CTVdVziH2uA64DSE5OHj516tRmZSwpKSEhIeGo+9V4lJu+LGNUahRXDfbf3dCB0Nhzdspj2RXkFnt49KRWRPjofoxgP2d/CLVzXpRXw5NLKzmtexQXDzz0v8FgPecDVcp988sR4L7RrUiM8V33YEvOefz48dmqOuJHT6iqax/ABUBhg20D8XahdT7Ma7oAq4E1QMfGHGf48OHaXLNmzWrUfnPWFWj3O6frF6vymn2sYNHYc3bKB0u2a/c7p+s3m/b47D2D/Zz9IZTOOX9/uWbd/5lO/udcrayuPex+wXzOy3L3at/ff6yXvbBQa2s9PnvflpwzsFgP8Znaovt4RCRFRG4Tkb4teZ8W2AG0F5H6d1qmAJXAj+YTr2sNzcVbmMap6p6ApGyEGavziYuO4IQ+tvaOv00cmExcdAQf2hQ6Bu+X79++v4Kyqlr+/tMsYqLceV/90C7tuG9yBnPXF/Lc3CZfig+olv4Xfhh4FO+w6e+JyFUi8hcRadPC9z+apUAV3uHUB40BslX1B4PbRaQD8AWwHzhJVfP9nK3RVJUZawoY0ycpZBd+CibxsVFMHJjMxyvyqK79US+rCTNvL85l5toC7jx9AH06B183WlNcPLIrZwxJ4dHPcliau8/pOIfV0sKzD7gG74ix76nqv4GXgOdFpEcLj3FYqloGvAI8JSIjRWQKcDvwBHzfIju4KtODQCfgSiCq7rkUEXH8Sv7avGJ27CvnlEG2BEKgTM5Mo6i0iq83FDodxTgot6iMBz5azfG9OnLl6B5Ox2kxEeEv5w4luU0ct7y5hOKK4FwKpKWFJxaYpqovNXxCVdcDNwJ/aOExjuY2YBEwE3gGeEBV36p7bhf//z6dC4A2wJK67QcfT/o531HNWJ2PCEwYYMOoA2Vc/yQS46LsZtIwpqrc9f4KRIRHLhhKRJDer9NUbVtH8/hFWWzfW8a901Y5HeeQWjpI/Q7gRRFZAnwGfFd3QQkAVd0jIn7ty6hr9VxR92j4nNT7e9BePJmxJp/MLu1ISnT3aDY3iY2K5PSMFD5ZmUdFda0tLx6Gpi3dybwNhTxwdgZd2ofWTCHH9ujAzRP68o8v13NaRgqnDw6uJVZa2uI5BZgM/An4FiiqmxrndhEZVzcXWteWhgxl+QcqWLZ9P6fYTaMBNyUrjZLKGmatLXA6igmwfWVV/HH6arK6tuPSUd2djuMXN03oQ0ZaG+7+YAVFpVVOx/mBlhaea4GRQD/gEuCtur8/DHwJvAA80MJjhLSZdR96Nk1O4B3fqyOdEmKsuy0MPfTJWvaVV/Pnc4cE7ZQ4LRUdGcGjF2Syv7yae6etdDrOD7S08KxX1eXqnfvsLVW9QVX7471X5nrga8A98zg4YMbqfLq0b0W/ZHePpnGjqMgIzhySypdrC4L2IqzxvW83FzF1US7XjOnJoDR/D7x11sDUNtwyoS/Tl+/i4xXBsxZVSwuPikiXH21U3amqL+CdvfreFh4jZFVU1/L1xkJOHtA54CsaGq8pWWlU1Xj4fFXQjK43flRT6+HeaStJb9eKX0506vbDwLphXG+GpLflng9WsidI1u9paeG5D/iXiIxv+ISI3IN3aLO712j1o282F1FR7WHcABtG7ZRh3dqT3q6VdbeFiamLclmbV8zvzxxI65jgngDUVw52uRVX1HD/R8Exi3WLCo+qFgE/BYaIyKUNnj4Lb2Gy9ZsPY9baAuKiIzi+V0eno4QtEWFyZhrzNhQGzbdB4x/7y6r52+c5HNerA5OCbJSXv/VPSeTG8b35cNlOZuc4P5imxXNDqGqVqj6hqv9p8NSZeO+dua2lxwhFqsqsnAJG9+5kQ3kdNiUzjVqP8snKPKejGD96bMY674X2szLCsmv75+N60yspnnumraS8qtbRLI1Zj6eDiDR5kLuqFqrqe6pacpj37dbU9wwlmwtL2bqnjPH9k46+s/GrgamJ9OmcYN1tIWx9fjGvLdzKxSO7hfyAgsOJjYrkz+cOIbeonH98ud7RLI1p8UQBL9UtJeATInI+cJev3s+NZuXsBmBcf7u+4zQRYUpmGou2FLFrf7nTcYwfPDB9NfExkdx2Sj+nozjquF4d+emILjw/dxNrdjm3GOJRC4+qFgC/B94XkZ9JC9qoIpIuIk8DU4Cbm/s+oWB2TgF9OifQtUNo3THtVlMy01CF6cuCZ8ip8Y2563czd30ht5zcl44JNjvIXZMG0rZVNHe9vwKPx5n12Bp1jUdVN+JdifREIEdEficiWY0pQiKSKCKni8i/ge+AZar6s4azR4eT0soavtlUZN1sQaRHp3iGdmlr3W0hxuNRHvpkLV3at+Ly40NzhoKmah8fwz1nDWRp7j7+822LF5ZulkaPJ1TVA8A1IjIM7zDpu4EaEVmM9ybRfXiXHIgBOgDtgZ7AUGA33tmqB6vqbp+egQt9vaGQqloP420YdVCZPDSNBz9ew+bCUnp2cseyx+bIPlq+k1U7D/DYhZnERtkgnoPOyUrnncXbefSzHM4ckkqH+MAOPm7yqDZV/U5VLwGS8S4xkI230IwFLgLOBoYANcB7wHggXVXvtqLjNStnNwmxUYzo3sHpKKaeszJTEYEPl1qrJxRU1Xh49PMcBqa24ezMdKfjBBUR4f4pGZRW1vDIZ2sDfvxm30GlqsXA+3UP00iqyuycAsb06eTalQ5DVWrbVhzbowMfLtvBLSf3Ccsht6HkjW+2kltUzstXDQ6ZJQ98qW9yIleO7sGLX2/momO7kdm1XcCObZ98AZaTX8yu/RWMH2DXd4LRlMw0Nu4uZbWDI35MyxVXVPPEzA0c36sjJ/Wzf2uH88uJfemUEMu9H64K6EADKzwBdnA2ahtGHZzOGJJKVITYIAOXe37uZopKq/jtpAHWcj2CxLhofnfGAJbl7uPd7MDN52yFJ8Bmr91NRlobktvEOR3FHEKH+BjG9O3E9GW7HBtqalqmqLSKF+Zu4owhKQHtPnKrc7LSObZHex76dC37ywIzS3tAC4+InBnI4wWb/WXVZG/by3hr7QS1KZlp7NhXznfb9jodxTTDs3M2Ul5dy68mhvfNoo0lItw3JYN9ZVX8/YucgBwz0C2ecwN8vKAyd8Nuaj1q13eC3KkZKcRGRVh3mwsVllTy6vytTMlMo29yotNxXCMjrS2XHded1xZuZfVO/1/f9Mm84CLyJNCqEbueA1zji2O60ay1u2nXOpqsru2djmKOICE2ipMHdubjFbu496xBREVaj7RbPDdnE5U1tdxycnisteNLvz6lPx8t28kD01fx5rXH+fXamK/+Rc0HBgNylEfYdpp7PMpX6wo4qV9SyC61G0qmZKZRWFLFgk17nI5iGqmguIJXF2zh7Kx0eifZir5N1bZ1NL8+tT8LNxXxqZ9navdJi0dV/yMix6jq7UfaT0Re9MXx3GjFjv0UllTZ9R2XGNe/M4mxUXy4dCdj+1rXqBs8+9UmqmvVWjstcNGxXXl94VYe/HgN4wd09tuSLY1u8YjIG0fZ5dtGvM1/G3u8UDMrpwARONHuKXCFuOhITs1I4dNVeVTWOLt2iTm6ggMVvL5wK+dkpdt0Ry0QFRnBvWcNYvvecl6ct9lvx2lKV9vJInLYCX1U9e2jvYGqTm/C8ULKrJzdZHVtF/A5kUzzTclKo7iihtk5NtNTsHtq9kZqPMotJ/dxOorrje7TidMzUnhy1gbyD1T45RhNKTxJwGYReUVErhCRrn5JFIIKSypZvn2fdbO5zFBFivwAABiXSURBVAm9O9IhPsZGtwW5vP0VvPHtNs4blk73jtba8YXfnTGQmlrlr5/6Zx63pg4u6Ahcjnem6S0isk5EnhGRn4pIJ9/HCw1f5exGFSbYbNSuEhUZwRlDUvhyTT6llWG7ikfQe3r2Bjwe5eYJdm3HV7p1bM01Y3vy0bKdFFV4fP7+TSk8uUAicBJwPzAP6AZcB7wJ5InIMhF5TEQmiYiNQa0zK6eApMRYBqWG55K7bjYlM52Kag9frM53Ooo5hIIDFby5KJfzhnWxRRV97MbxffjklyfSIc73H+VNecciVa1W1bmq+oCqnoR3zZ3TgIeBJcAg4JfAdGC1iAz0eWKXqfUoc9btZly/JJsh14VGdG9Pats4624LUi/M20xNrYefj+vtdJSQkxAbRZ/O/hmW3pTC82bDDaparqpfqOpdqnos3q64c4BngM7AZyIS1ndLbtzv4UBFjS365lIREcLkzDTmrNvN3tIqp+OYevaWVvH6wq1Mzkyjh41kc5VGFx5VfbgR+xxQ1Q9V9RdAP7wrkx7x3p5Qt6yglqgIYUxfuwTmVlMy06jxKJ/4+aY60zQvz99CWVUtN46zkWxu47frMKpaCJwKjPLXMdxgeWEtI3q0p01ctNNRTDNlpLWhV6d4Ply2w+kopk5JZQ0vz9/CKYOS6Z9ic7K5jV8HAKhqCbDRn8cQkVgReU5E9opInoj85gj7ZorIAhEpE5FsETnWn9l27S8nt9hjw6hdTsTb3fbN5iK/3ddgmub1hVvZX17NTeOtteNGPi88IjJRRPaJyP11m/w9adIjwGhgInA9cLeIXHSIXPHAJ8BCYDgwF/ifiPjt69LBGw/t+o77TclKQxWmL9/ldJSwV1FdywtzNzO2bydbb8el/NHiiQHaAHeKyJvAND8cA/i+mFwL3Kqq2ao6De8Iu5sOsfuFQDXwa1VdA/wK2F+33S9mri2gY5zQ108jQ0zg9E5KICOtjY1uCwJvLcqlsKSSX1hrx7V8XnhU9WPgdOD3wGONmUqnBTKBWLz3FB00DzhWRBpOgHoc8LWqeupyKvA1cLw/gqkq5VW1ZHWOtKV3Q8SUzDSW5e5j655Sp6OEraoaD89+tZHh3dszqmcHp+OYZhLv5687ich5wLOq2qnetoHAaiBNVXfV2/4RkFN/Bm0R+SuQpaqnHeK9r8N7cyzJycnDp06d2qyMB4pLaJMYXi2ekpISEhJC75z3lHv49Vfl/KRvNFN6/3DOvVA95yNx4pznbq/mxZVV/Gp4LJlJPplcv0ns99w048ePz1bVEQ23B/4351utgcoG2w7+HNvIfRvuB4CqPgc8BzBixAgdN25cswLOnj2b5r7WrUL5nKdunc/KA9X8fdxJP9geyud8OIE+51qPcv/fvyIjLY5bzh/jSE+C/Z59w+3T2lTw48Jx8OeyRu7bcD9jDmtKZhrr8ktYm+f/5YHND328YhebC0v5xfg+1n3tcm4vPDuA9g2Wa0jB25IpOsS+KQ22pQA2TMk02qQhqURGCB8utUEGgaSqPDlrA72S4jkto+E/Y+M2bi88S4EqvMOpDxoDZKtqw+mEFwKjpe6rUt2fo+u2G9MonRJiGd27Ix8t34mbr4+6zZdrClibV8yN4/rY0vEhwNWFR1XLgFeAp0RkpIhMwTtFzxMAIpIiIq3qdn8X7z1F/xSRQcDf8Q77bt6oARO2pmSmkVtUzpLcfU5HCQuqyr9mbaBL+1acnZXmdBzjA64uPHVuAxYBM/FOTvqAqr5V99wu6u7TUdUDwJl4WznfAScAZ6hqccATG1c7bXAKMVER1t0WIAs27mFp7j6uP6k30ZGh8JFl3D6q7WCr54q6R8PnpMHPi4BhAYpmQlSbuGjG90/ifyt2cc9Zg6zrx8+enL2BpMRYLhjexekoxkfs64MxzXB2Vjq7iyuZv7HQ6SghbWnuPr7esIdrxvQkLjrS6TjGR6zwGNMMEwZ0pk1cFO9mb3c6Skh7atYG2sRFcelx3Z2OYnzICo8xzRAXHcmUrDQ+XZnHgYpqp+OEpPX5xXy+Op8rR/cgIdb1VwVMPVZ4jGmmC4Z3pbLGw/9sxmq/eHr2RlpFR3LlCT2djmJ8zAqPMc00tEtb+nZO4J3FuU5HCTm5RWVMW7aTi0d2o0N8zNFfYFzFCo8xzSQinD+8C99t28euEo/TcULKc3M2ESFw7YnW2glFVniMaYFzj0knMkKYt6PhRBmmuQqKK3hrcS7nDetCattWR3+BcR0rPMa0QOc2cZzUL4mvd9ZQ67EpdHzhpXlbqKn1cP1JvZ2OYvzECo8xLXTB8C7sq1TmbbB7elpqf3k1ry/cyhlDUunZKd7pOMZPrPAY00ITBnYmPhobZOADry3YQkllDT8fZ62dUGaFx5gWio2K5LjUKD5fnc/+Mrunp7nKq2p56estjO+fREZaW6fjGD+ywmOMD4xNj6KqxsOHy3Y4HcW1pi7aRlFpFb8Y38fpKMbPrPAY4wPd20QwKLUNb36ba+v0NENVjYfn5mxiZI8OjOjRwek4xs+s8BjjAyLCxaO6sXrXAZZv3+90HNf5YMkOdu2v4Mbxdm0nHFjhMcZHzs5Ko1V0JG98s83pKK5S61Ge+WojGWltOKlfktNxTABY4THGR9rERTMlM40Pl+2k2CYObbRPVu5iU2EpPx/Xm7qV6U2Is8JjjA9dPKob5dW1TLPVSRvF41H++eUGeifFM2lwqtNxTIBY4THGhzK7tGVQahve+GabDTJohM9X55GTX8zNE/raSq5hxAqPMT5kgwwaT1X5x5cb6NUpnsmZaU7HMQFkhccYHzvHBhk0yher81mz6wC/GN/HWjthxgqPMT6WWG+Qga1OemiqyhMz19O9Y2vOzrLWTrixwmOMH1x6nHeQwbuLtzsdJSjNyilg5Q5vaycq0j6Gwo39xo3xg6Fd2jGsWzteXbAFjy2X8AOqyj9mrKdL+1ace0y603GMA6zwGOMnV57Qky17yvhq3W6nowSVr9btZtn2/fxifB+irbUTluy3boyfTBqcQnKbWP49f4vTUYKGdyTbetLbteK8YV2cjmMcYoXHGD+JjozgslHdmbNuNxsKSpyOExTmrC9kybZ93DCuNzFR9vETruw3b4wfXTyqGzGREby6YIvTURynqjz6WQ7p7Vpx4YiuTscxDrLCY4wfdUqIZXJmGu9mbw/7odWfrcpjxY793Dqxr7V2wpz99o3xsytH96Csqpa3F4Xv0ti1HuVvn6+jV1K8jWQzVniM8bchXdoysmcHXpq3mepaj9NxHDFt6Q7WF5Tw61P62307xt2FR7weFJECEdkrIo+KSOQR9j9ZRBaJSImI5IjI1YHMa8LXDSf1Yuf+Cj4Mw1mrq2o8PD5jPYNS2zBpcIrTcUwQcHXhAX4FXAlcAJwLXAzccagdRaQvMB34L5AFPAA8KSKTA5LUhLXx/TvTPzmRZ+dsDLtZq99enMu2ojJuP60fETYnm8H9hedW4A+q+pWqzgbuBH5xmH0vBJaq6p9VdYOq/gd4Fbg0MFFNOBMRrj+pF+vyS5iVU+B0nICpqK7lnzPXM7x7e8b37+x0HBMkXFt4RCQN6ArMqbd5HtBFRA41VvNt4KYG2xSI809CY35ocmYaaW3jeOarTU5HCZiXvt5M/oFK7jitv60uar7n2sIDHFyusH6neX7dnz+6JVpV16lq9sGfRSQZuIgfFi5j/CY6MoKrx/bi281FfLdtr9Nx/K6wpJKnZm1k4sDOHNero9NxTBCRYO5vFpE4DlFE6qTiLRoxqlpdt38EUAuMr+t6O9z7xgOfAx2A4apadoh9rgOuA0hOTh4+derUZp1DSUkJCQkJzXqtW9k5H15FjfLrr8ro0y6SXw13d2P7aOf86upKZufW8KcTWpGW4ObvuP+f/b/dNOPHj89W1RE/ekJVg/YBjMHbHXaoxx11fybU279V3baRR3jPtsBcYBfQuzE5hg8frs01a9asZr/Wreycj+yfX67T7ndO16Xb9vovUAAc6ZzX5xdrr7v+p3f/d0XgAgWA/b/dNMBiPcRnalB/DVHVeaoqh3oA/6nbrf74zIN/33Wo9xORTsAsoBcwTlU3+i28MYdxxegetGsdzeMz1jkdxW8e+mQNraIj+eXEvk5HMUEoqAvPkajqTmAb3lbRQWOAnar6o1vERSQG73DqTsCJqpoTkKDGNJAYF821Y3sxK2c3S0LwWs/8jYXMWFPAjeN70ykh1uk4Jgi5tvDUeRr4i4hMEJGTgL8A/zj4pIgkicjBzslfAcOBq4BSEUmpe3QIeGoT9q4Y3YP2raN5fMZ6p6P4VE2thwc+Wk16u1b83wk9nY5jgpTbC88jwBvAe3WPN4FH6z2/CLi97u8XAFHADLxdcQcfHwYqrDEHJcRGcd2Jvflq3W4WbylyOo7PvLpgK2vzirnnrIHERR92EhET5lxdeFS1VlV/rartVbWTqt6hqp56z/dQ1fvq/j7iMNeLxhz2AMb40RWju5OUGMufP14TErMZFBRX8NgX6zixXxKnZdjUOObwXF14jHGz1jFR/PqUfny3bR+frsxzOk6LPfTJWipqarlv8iC7WdQckRUeYxx0wYiu9EtO4KFP11JV496ZqxdvKeL973Zw7dhe9EoKr/tcTNNZ4THGQZERwl2TBrJ1Txn/+War03GapbKmlrveX0Fa2zhumtDH6TjGBazwGOOwcf2TOKFPR/7x5Xr2llY5HafJnpy1kfUFJTz4kyG0jolyOo5xASs8xjhMRLj7zEEUV9Tw8Gfuur1sza4DPDVrAz85Jt1mnzaNZoXHmCAwMLUNV43uwdRF21wzgWitR7nzveW0bRXNPWcNcjqOcRErPMYEiVtP6UdyYhz3fLCSGhcskf2/zdUs376f+8/OoH18jNNxjItY4TEmSCTERnHv5EGs2nmAl+dvcTrOEX23bS8fbKhmSmYaZw5JPfoLjKnHCo8xQWTS4BQmDuzMI5/lsKGgxOk4h1RcUc2tU5fSPlb407mD7Z4d02RWeIwJIiLCn38yhFYxkfz6nWVB2eX2hw9XsX1vGTdkxtImLtrpOMaFrPAYE2Q6J8bxp3MGsyx3H0/PDq6VO978dhvvf7eDmyb0pW97m4vNNI8VHmOC0FlD05icmcbjX65nUZBMIvrdtr38YdoqxvbtxC9PtnV2TPNZ4TEmSD147mC6tm/FL/7zHbuLKx3NUlBcwc9fzya5bSz/vPgYIiPsuo5pPis8xgSpNnHRPHXpcPaXV3PLm0uoduh6T0llDVe/vJj95dU8e9kI2rW2odOmZazwGBPEBqW14c/nDmHBpj3c/d+VAV8+oarGw89fz2b1rgM8eckwBqW1CejxTWiyiZWMCXLnDe/C1j2lPDFzA13at+LmAF1fqfUod7y7jLnrC3n4/KGcPDA5IMc1oc8KjzEu8KtT+rF9bzl/+2IdCXFRXOXnZaWraz3c9vYyPlq2k9+c3p+fjujq1+OZ8GKFxxgXEBEeOm8opVU13P/Ramo9yjVje/nlWBXVtfxy6hI+W5XPbycN4IaTevvlOCZ82TUeY1wiJiqCf10yjDOGpPCn/63hr5+uxePx7TWf/AMVXPjsAj5blc8fJg+yomP8wlo8xrhIdGQET1x0DG1breLp2RtZn1/CoxcM9clIs/kbC/nVW0spqajhucuHc2pGig8SG/Nj1uIxxmWiIiP487mDuW/yIGbnFHDa43OYnVPQ7Pc7UFHNfR+u4pLnvyE+Jor3bhxtRcf4lbV4jHEhEeHKE3oyrHt7bnt7GVf+exHj+ydx+2n9yUhr26j3KK2s4e3Fufxz5gaKSqu44vju/HbSQFrF2FQ4xr+s8BjjYkO7tGP6zWN4Zf4W/jVrA2c+MY/h3dtzdlYao3t3oleneCLqzTJQVFrFkm17mbGmgOnLdlJcWcOonh24+8xBDOnSuIJlTEtZ4THG5eKiI7n+pN5ceGxX3s3ezhvfbuPeaasAiI4UOifGAVBaVcO+smoA4mMiOXlgMleM7sGwbu1saQMTUFZ4jAkR7VrHcM3YXlw9pifbispYuGkPW/aUkb+/AgRaRUfSs1M8A1PbMKJHe2KjrEvNOMMKjzEhRkTo3jGe7h3jnY5izCHZqDZjjDEBZYXHGGNMQFnhMcYYE1BWeIwxxgSUFR5jjDEBZYXHGGNMQFnhMcYYE1BWeIwxxgSUBHoNdzcSkd3A1ma+vBNQ6MM4bmDnHB7snMNDS865u6omNdxohcfPRGSxqo5wOkcg2TmHBzvn8OCPc7auNmOMMQFlhccYY0xAWeHxv+ecDuAAO+fwYOccHnx+znaNxxhjTEBZi8cYY0xAWeExxhgTUFZ4/EBEYkXkORHZKyJ5IvIbpzP5m4j0FpGP6s55u4j8TUTinM4VKCLygojMdjqHv4lItIj8XUQKRWSPiDwtIrFO5/InEWkvIq+LSJGI7BCRh0QkJJdvrfvsWikiE+tt6yAi74jIARHZIiJXtPQ4tgKpfzwCjAYmAl2A10Rkm6pOdTaWf4hIDPARsBrveXcGXqp7+tdO5QoUETkZuBr4yuksAfAIcA5wNqDAG8Ae4G4nQ/nZU0AqcCKQxP8/50ecDOVrdV8U3wAyGjz1MpAAnAAcCzwrIutVdX6zj2WDC3xLROLx3uU7WVVn1G27GzhdVcc4Gs5PRGQMMBPooKolddsuAf6uqimOhvOzut/3cmAXUKOq45xN5D8i0g7IB85S1S/qtl0JXKiqk5zM5k8ish+4QlU/qPv5b8CgUDpnERmEt+gIMBQ4RVVniEhvYAPQV1U31O37AhCnqpc193jW1eZ7mUAsMK/etnnAsSISqi3MHOCMg0WnjuL97xDqHgRm1z1C3RigHJhxcIOqvhxKH8CHsQe4VERai0gacDqQ7XAmXxsLfA4c32D7KGDXwaJTZ94h9msSKzy+lwoUqWpFvW35QAzeZnrIUdXdB1t3ACISAdwEzHUulf+JyPHABcDtTmcJkN7AFuBiEVklIltF5NG6rtZQdiMwDigGdgB5wH0O5vE5VX1WVX+jqmUNnkoFdjbYlo/3EkKzWeHxvdZAZYNtB38OhxYAwN+BY4C7nA7iL3UX1F8EblXVvU7nCZBEoCdwM3A9cANwPvCwk6ECoA+wBO81njOAHsCjTgYKoMN9nsWIiDT3TUO168dJFfy4wBz8ueG3iZBS9z/i43i/IZ6vqqscjuRP9wLrVfUdp4MEUA3QBrhMVTcCiMjteAfP3KaqHkfT+UHdNY7HgR6qur1u2zXAFyLyF1XNdzSg/x3u86xcWzBAwAqP7+0A2otIjKpW1W1Lwfstoci5WP5V1732InAp3ovN0xyO5G+XAKkicvC6VgwQKSIlqprgYC5/2ol3AMXGettygDi83cih+CE8HCg+WHTqZAORQHdC85zr24H386u+FLyDaZrNutp8bylQhXdY8UFjgGxVrXEmUkD8De+H8U9U9X2nwwTAOGAwkFX3eB5YXPf3ULUAiBKRIfW2DcJ77WOPM5H8bifQTkS61ts2sO7PTQ7kCbSFQLqI9Ki3bUzd9mazFo+PqWqZiLwCPFU31DQF78Xn6xwN5kcichxwK95rOotF5PtvSKqa51gwP1LVHywMKCJ78XY/bDjMS1xPVdeLyDTg3yJyPd7+/4eA50P4S9VCvF8m/y0it+E952eB11Q15BeEU9VNIvIZ8KqI3IS3BXgpML4l72stHv+4DViE996WZ4AHVPUtZyP51fl1f/4FbxP8+0cIDyEPV5fjvW9pJvAB8F/gt44m8qO6gnom3m7ymcD7eG8Uvt7JXAH2M2Af8A3ea5vXqOqClryh3UBqjDEmoKzFY4wxJqCs8BhjjAkoKzzGGGMCygqPMcaYgLLCY4wxJqCs8BhjjAkoKzzGGGMCygqPMcaYgLLCY4wxJqCs8BhjjAkoKzzGGGMCygqPMcaYgLKZg41xERHpCNwHCNAX7zpAXwCP4F1ssB1wp6rudCqjMUdjhccYlxCRWOAl4GZV3SYimcC3wHTgBmAK8AKwDHjUsaDGHIV1tRnjHjcAT6jqtrqfy/Auub1UVXfXbVsOfOREOGMaywqPMe5RpKpf1vt5WN2fnwKo6ouqmqmqOYGPZkzj2UJwxriUiDwDXAx0UNVap/MY01jW4jHGvSYA86zoGLexwmOMC4lIOt5RbV812P5/ziQypvGs8BjjAiKSJCLfisgf6jZNqvtzcb19+gH9Ax7OmCaywmOMO5wEHAuIiMQDZwKFQCJ8f3/Pg8BfHEtoTCPZ4AJjXEBEEoHHgCqgNfAA0AW4F8jF+yXyflXd5FhIYxrJCo8xxpiAsq42Y4wxAWWFxxhjTEBZ4THGGBNQVniMMcYElBUeY4wxAWWFxxhjTEBZ4THGGBNQVniMMcYElBUeY4wxAfX/ABGzkXP98oJVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.special import j1\n",
    "\n",
    "x_values = np.arange(0.01,10,0.01)\n",
    "j1_values = np.zeros(len(x_values))\n",
    "for idx, x in enumerate(x_values):\n",
    "    j1_values[idx] = j1(x)\n",
    "    \n",
    "plt.plot(x_values,j1_values)\n",
    "plt.xlabel(r'$x$',fontsize=22)\n",
    "plt.ylabel(r'$J_1(x)$',fontsize=22)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We note that there is a root just before 4. And, so we use that as our trial value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.8317059702075125"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import root\n",
    "\n",
    "\n",
    "def myfun_besselj(Kvalue):\n",
    "    return j1(Kvalue)\n",
    "\n",
    "sol = root(myfun_besselj,4)\n",
    "Kvalue_soln = sol.x[0]\n",
    "Kvalue_soln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, we must have for the first critical value $\\displaystyle \\bar{K} = \\sqrt{\\frac{P}{D}}R = 3.832$, and so the first critical value for the externally applied force per unit length, $P$ is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle P_{\\rm crit} = \\frac{14.68 D}{R^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R = sp.symbols('R')\n",
    "P_crit = sp.N(3.832**2*D/R**2,4) # This ensures that the result is shown correct to 2 decimal places\n",
    "display(Math(r'P_{{\\rm crit}} = {}'.format(sp.latex(P_crit))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above result matches with that given on page 390 (Section 9.9) of Theory of Elastic Stability, 2nd edition by Timoshenko and Gere."
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}