{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Universidade Federal do Rio Grande do Sul (UFRGS)   \n",
    "Programa de Pós-Graduação em Engenharia Civil (PPGEC)   \n",
    "\n",
    "# PEC00025: Introduction to Vibration Theory\n",
    "\n",
    "\n",
    "### Class 14 - Vibration of beams\n",
    "\n",
    "[1.   The vibrating beam equation](#section_1)   \n",
    "[2.   Free vibration solution](#section_2)   \n",
    "[3.   Vibration modes and frequencies](#section_3)  \n",
    "[4.   Solution by approximation](#section_4)   \n",
    "[5.   Assignment](#section_5)  \n",
    "\n",
    "---\n",
    "_Prof. Marcelo M. Rocha, Dr.techn._ [(ORCID)](https://orcid.org/0000-0001-5640-1020)  \n",
    "_Porto Alegre, RS, Brazil_ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing Python modules required for this notebook\n",
    "# (this cell must be executed with \"shift+enter\" before any other Python cell)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. The vibrating beam equation <a name=\"section_1\"></a> \n",
    "\n",
    "The static analysis of a beam under Bernoulli's hypothesis (plane sections remain plane after beam\n",
    "deformation) leads to the well known differential equations:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{dQ}{dx}   &= -q(x) \\\\\n",
    "\\frac{dM}{dx}   &=  Q(x) \\\\\n",
    " EI \\psi^\\prime &=  M(x)\n",
    "\\end{align*}\n",
    "\n",
    "where $q(x)$ is the distributed transversal loading, $Q(x)$ is the shear, $M(x)$ is the bending\n",
    "moment, $\\psi(x)$ is the section rotation, and $EI$ is the flexural stiffness (regarded \n",
    "hereinafter as constant along the beam length).\n",
    "\n",
    "<img src=\"images/dynamic_beam.png\" alt=\"Dynamic beam equilibrium\" width=\"360px\"/>\n",
    "\n",
    "Disregarding shear strains, $\\gamma(x) = 0$, implies that the section rotation is approximated as:\n",
    "\n",
    "$$  \\psi \\approx -w^\\prime(x) $$\n",
    "\n",
    "what implies that:\n",
    "\n",
    "\\begin{align*}\n",
    " EI w^{\\prime\\prime} &\\approx -M(x) \\\\\n",
    " EI w^{\\prime\\prime\\prime\\prime} &\\approx  q(x)\n",
    "\\end{align*}\n",
    "\n",
    "where $w(x)$ is the beam transversal displacement, also called _elastic line_. The solution of this\n",
    "last differential equation is straighforward once the load $q(x)$ and the boundary conditions (two\n",
    "for each beam extremity) are specified. \n",
    "\n",
    "We shall now include the inertial forces in these equations, as well as regard now the section mean shear strain,\n",
    "\n",
    "$$ \\gamma(x) = \\psi(x) + w^\\prime(x) = \\frac{Q(x)}{GA_{\\rm s}}$$\n",
    "\n",
    "as relevant, where $GA_{\\rm s}$ is the shear stiffness (also regarded hereinafter as constant along\n",
    "the beam length). Although $\\gamma(x)$ is indeed negligible for actual slender beams, \n",
    "the following analysis may also be applied to other slender structures, like\n",
    "tall buildings, trusses, etc., for which we may define _equivalent stiffnesses_, \n",
    "$\\left(EI\\right)_{\\rm eq}$ and $\\left(GA_{\\rm s}\\right)_{\\rm eq}$.\n",
    "The dynamic equilibrium equations now become:\n",
    "\n",
    "\\begin{align*}\n",
    " Q^\\prime &= -q +   \\mu \\ddot{ w}   \\\\\n",
    " M^\\prime &=  Q - i_\\mu \\ddot{\\psi}\n",
    "\\end{align*}\n",
    "\n",
    "where $\\mu$ is the beam mass per unit length and $i_\\mu$ is the cross section rotational inertia\n",
    "per unit length. Derivating the equation for $\\gamma$ and replacing the solution \n",
    "$EI \\psi^{\\prime} = M(x)$ gives the elastic line equation accounting for shear strains:\n",
    "\n",
    "$$ w^{\\prime\\prime} = -\\frac{M}{EI} + \\frac{Q^\\prime}{GA_{\\rm s}} $$\n",
    "\n",
    "Now we replace the shear (with inertial loads):\n",
    "\n",
    "$$ w^{\\prime\\prime} = -\\frac{M}{EI} + \\frac{-q +   \\mu \\ddot{ w}}{GA_{\\rm s}} $$\n",
    "\n",
    "and derivate the whole equation replacing for $M^\\prime$:\n",
    "\n",
    "$$ w^{\\prime\\prime\\prime} = \\frac{i_\\mu \\ddot{\\psi} - Q}{EI} + \\frac{\\mu \\ddot{w}^\\prime - q^\\prime}{GA_{\\rm s}} $$\n",
    "\n",
    "The angular acceleration, $\\ddot{\\psi}$, may be safely disregarded, for the rotations are usually very small.\n",
    "Derivating the equation a last time and replacing for $Q^\\prime$ finally gives:\n",
    "\n",
    "$$ EI w^{\\prime\\prime\\prime\\prime} = q - \\mu \\ddot{ w} + \n",
    "\\frac{EI}{GA_{\\rm s}} \\, \\left(\\mu \\ddot{ w}^{\\prime\\prime} - q^{\\prime\\prime} \\right) $$\n",
    "\n",
    "which is the dynamic elastic line equation for a constant section beam under forced vibration due to \n",
    "dynamic load $q(x,t)$, with shear deformation accounted for (although plane section hypothesis still kept along). \n",
    "The last term may be disregarded whenever the shear stiffness is much larger that the bending stiffness.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Free vibration solution <a name=\"section_2\"></a> \n",
    "\n",
    "In this section, we take the vibrating beam equation derived above, disregard the shear deformation and look for\n",
    "free vibration solution, which implies that $q(x, t) = 0$. The beam equation becomes simply:\n",
    "\n",
    "$$ EI w^{\\prime\\prime\\prime\\prime} = - \\mu \\ddot{w} $$\n",
    "\n",
    "To solve this equation we separate the time and space independent variables through the hypothesis:\n",
    "\n",
    "$$ w(x,t) = w_0 \\sin \\left(\\omega t + \\theta \\right)\\, \\varphi(x)\\,   $$\n",
    "\n",
    "which is pretty much alike we have previously done for multiple degrees of freedom systems.\n",
    "The free vibration equilibrium equation then become:\n",
    "\n",
    "$$ \\varphi^{\\prime\\prime\\prime\\prime} - p^4 \\varphi = 0 $$\n",
    "\n",
    "where we have defined:\n",
    "\n",
    "$$ p^4 = \\left(\\frac{\\mu}{EI}\\right) \\omega^2 $$\n",
    "\n",
    "It can be shown that, in the general case, the space dependent function $\\varphi(x)$ has the form:\n",
    "\n",
    "$$ \\varphi(x) = C_1 \\left(\\cos px + \\cosh px \\right) + \n",
    "                C_2 \\left(\\cos px - \\cosh px \\right) + \n",
    "                C_3 \\left(\\sin px + \\sinh px \\right) + \n",
    "                C_4 \\left(\\sin px - \\sinh px \\right) $$  \n",
    "\n",
    "The corresponding space derivatives will be required to apply the boundary conditions:\n",
    "\n",
    "\\begin{align*}\n",
    "\\varphi^\\prime(x) = p^1&\\left[C_1 \\left(-\\sin px + \\sinh px \\right) + \n",
    "                              C_2 \\left(-\\sin px - \\sinh px \\right) + \n",
    "                              C_3 \\left( \\cos px + \\cosh px \\right) + \n",
    "                              C_4 \\left( \\cos px - \\cosh px \\right)\\right] \\\\\n",
    "\\varphi^{\\prime\\prime}(x) = p^2&\\left[C_1 \\left(-\\cos px + \\cosh px \\right) + \n",
    "                                      C_2 \\left(-\\cos px - \\cosh px \\right) + \n",
    "                                      C_3 \\left(-\\sin px + \\sinh px \\right) + \n",
    "                                      C_4 \\left(-\\sin px - \\sinh px \\right)\\right] \\\\\n",
    "\\varphi^{\\prime\\prime\\prime}(x) = p^3&\\left[C_1 \\left( \\sin px + \\sinh px \\right) + \n",
    "                                            C_2 \\left( \\sin px - \\sinh px \\right) + \n",
    "                                            C_3 \\left(-\\cos px + \\cosh px \\right) + \n",
    "                                            C_4 \\left(-\\cos px - \\cosh px \\right)\\right] \\\\\n",
    "\\varphi^{\\prime\\prime\\prime\\prime}(x) = p^4&\\left[C_1 \\left( \\cos px + \\cosh px \\right) + \n",
    "                                                  C_2 \\left( \\cos px - \\cosh px \\right) + \n",
    "                                                  C_3 \\left( \\sin px + \\sinh px \\right) + \n",
    "                                                  C_4 \\left( \\sin px - \\sinh px \\right)\\right]\n",
    "\\end{align*}\n",
    "\n",
    "The last equation above proves that the assumed general solution is correct.\n",
    "\n",
    "Now, to have a particular solution for the vibrating beam the kinematic boundary conditions must be applied.\n",
    "Let us assume a cantilever beam, fixed at the left end ($x = 0$) and free at the right end ($x = L$).\n",
    "\n",
    "<img src=\"images/cantilever_beam.png\" alt=\"Cantilever beam\" width=\"360px\"/>\n",
    "\n",
    "The corresponding boundary conditions are:\n",
    "\n",
    "\\begin{align*}\n",
    " \\varphi(0)                      &= 0 \\\\\n",
    " \\varphi^ \\prime(0)              &= 0 \\\\\n",
    " \\varphi^{\\prime\\prime}(L)       &= 0 \\\\\n",
    " \\varphi^{\\prime\\prime\\prime}(L) &= 0\n",
    "\\end{align*}\n",
    "\n",
    "The two last conditions implies that bending moment and shear force are zero at the right end, respectively.\n",
    "Applying these conditions at the corresponding derivatives:\n",
    "\n",
    "\\begin{align*}\n",
    "\\varphi(0) &=  C_1 \\left( 1 + 1 \\right) + \n",
    "               C_2 \\left( 1 - 1 \\right) + \n",
    "               C_3 \\left( 0 + 0 \\right) + \n",
    "               C_4 \\left( 0 - 0 \\right) = 0\\\\\n",
    "\\varphi^\\prime(0) &= p \\left[C_1 \\left(-0 + 0 \\right) + \n",
    "                             C_2 \\left(-0 - 0 \\right) + \n",
    "                             C_3 \\left( 1 + 1 \\right) + \n",
    "                             C_4 \\left( 1 - 1 \\right)\\right] = 0\n",
    "\\end{align*}\n",
    "\n",
    "which implies that $C_1 = 0$ and $C_3 = 0$. The other two conditions become:\n",
    "                              \n",
    "\\begin{align*}                              \n",
    "\\varphi^{\\prime\\prime}(L) &= p^2 \\left[C_2 \\left(-\\cos pL - \\cosh pL \\right) + \n",
    "                                       C_4 \\left(-\\sin pL - \\sinh pL \\right)\\right] = 0\\\\\n",
    "\\varphi^{\\prime\\prime\\prime}(L) &= p^3 \\left[C_2 \\left( \\sin pL - \\sinh pL \\right) + \n",
    "                                             C_4 \\left(-\\cos pL - \\cosh pL \\right)\\right] = 0\n",
    "\\end{align*}\n",
    "\n",
    "These two equations can be put into matrix form as:\n",
    "\n",
    "$$ \\left[ \\begin{array}{cc}\n",
    "                 \\left(-\\cos pL - \\cosh pL \\right) & \\left(-\\sin pL - \\sinh pL \\right) \\\\\n",
    "                 \\left( \\sin pL - \\sinh pL \\right) & \\left(-\\cos pL - \\cosh pL \\right) \n",
    "          \\end{array} \\right] \n",
    " \\left[ \\begin{array}{c}\n",
    "                  C_2 \\\\\n",
    "                  C_4 \n",
    "        \\end{array} \\right]  =                \n",
    " \\left[ \\begin{array}{c}\n",
    "                  0 \\\\\n",
    "                  0 \n",
    "        \\end{array} \\right] $$\n",
    "        \n",
    "In order to obtain a non trivial (non zero) solution for the unknown coefficients $C_2$ and $C_4$, \n",
    "the determinant of the coefficients matrix must be zero. This condition yields a nonlinear equation \n",
    "to be solved for $pL$. We can use the HP Prime for this purpose, as shown in the figure below.\n",
    "\n",
    "<table>\n",
    " <tr>\n",
    "   <td><img src=\"images/det_cantilever_1.jpg\" alt=\"HP Prime determinant\" width=\"320px\"/></td>\n",
    "   <td><img src=\"images/det_cantilever_2.jpg\" alt=\"HP Prime cantilever\"  width=\"320px\"/></td>\n",
    " </tr>\n",
    "</table> \n",
    "\n",
    "There will be infinite solutions $\\alpha_k = \\left( pL \\right)_k$, $k = 1, 2, \\dots, \\infty$, \n",
    "each one associated to a vibration frequency and a modal shape,\n",
    "$\\left[ \\omega_k, \\varphi_k(x) \\right]$. The natural vibration frequencies are obtained\n",
    "by recalling the definition of $p$, what finally gives:\n",
    "\n",
    "$$ \\omega_k = \\left( \\frac{\\alpha_k}{L} \\right) ^2 \\sqrt{\\frac{EI}{\\mu}}$$\n",
    "\n",
    "For instance, the fundamental vibration frequency, $f_{\\rm n}$, is given by:\n",
    "\n",
    "$$ f_{\\rm n} \\approx \\frac{1}{2\\pi} \\left( \\frac{1.8751}{L} \\right) ^2 \\sqrt{\\frac{EI}{\\mu}}$$\n",
    "\n",
    "which is a very useful formula for estimating the fundamental frequency of slender piles and towers\n",
    "with constant cross section.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 10, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Vibration modes and frequencies <a name=\"section_3\"></a> \n",
    "\n",
    "The table below was taken from a _Sonderausdruck_ (special edition) of the german _Betonkalender_ \n",
    "(concrete almanac), 1988. It summarizes the solutions for some other support conditions of slender beams. \n",
    "If more accuracy is desiredfor the $\\alpha_k$ constants, one can solve the so-called \n",
    "_characteristic equation_ with the help of a calculator like the HP Prime, ou HP 50G.\n",
    "\n",
    "<img src=\"images/beams.png\" alt=\"Beam solutions\" width=\"640px\"/>\n",
    "\n",
    "The characteristic equations are the determinants of respective coefficients matrix, \n",
    "which can also be solved with the ``fsolve()`` method from ``scipy``, as shown below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cantilever beam frequency parameter: 1.8751040687119394\n"
     ]
    }
   ],
   "source": [
    "def char_eq(x):\n",
    "    \n",
    "    x = x[0]\n",
    "    A = np.array([[-np.cos(x)-np.cosh(x), -np.sin(x)-np.sinh(x)], \n",
    "                  [ np.sin(x)-np.sinh(x), -np.cos(x)-np.cosh(x)]])\n",
    " \n",
    "    return np.linalg.det(A)              # from coefficientes matrix\n",
    "\n",
    "#   return np.cos(x)*np.cosh(x) + 1      # from characteristic equation\n",
    "\n",
    "#-----------------------------------------------------------------------\n",
    "\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "ak = fsolve(char_eq, 1)\n",
    "\n",
    "print('Cantilever beam frequency parameter: {0}'.format(ak[0]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that the result is exactly the same (within the required precision) as previously obtained with the \n",
    "HP Prime. One can use directly the characteristic equation, or one can program the determinant calculation\n",
    "using ``np.linalg.det()`` in the user function ``char_eq`` above.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Solution by approximation <a name=\"section_4\"></a> \n",
    "\n",
    "The solutions for the beam vibration frequencies presented above may also be calculated by means of\n",
    "Rayleigh quotient, as long as an educated guess for the $\\varphi(x)$ function is assumed.\n",
    "As an example, let us take a simply supporte beam, for which we assume:\n",
    "\n",
    "$$ \\varphi(x) = 4x \\left( \\frac{L - x}{L^2} \\right) $$\n",
    "\n",
    "with $\\varphi(L/2) = 1$, which obviously _is not_ the correct modal shape. \n",
    "\n",
    "<img src=\"images/simply_supported.png\" alt=\"simply supported beam\" width=\"400px\"/>\n",
    "\n",
    "We also assume that the beam is subjected to a constant distributed load, $q$ corresponding \n",
    "to its self weight. The maximum displacement at the beam center is known to be:\n",
    "\n",
    "$$ w_{\\rm max} = \\frac{5 q L^4}{384EI}$$\n",
    "\n",
    "In the following we shall estimate both the fundamental vibration frequency and the\n",
    "maximum displacement with the assumed modal shape. The reference kinetic energy is given by:\n",
    "\n",
    "$$ T_{\\rm ref} = \\frac{1}{2} \\int_0^L {\\mu \\varphi ^2(x) \\, dx} $$\n",
    "\n",
    "while for the elastic potential energy we need the curvature function, $\\varphi^{\\prime\\prime}(x)$:\n",
    "\n",
    "$$ V = \\frac{1}{2} \\int_0^L {EI \\left[ \\varphi^{\\prime\\prime}(x) \\right] ^2 \\, dx}\n",
    "     = \\frac{1}{2} \\int_0^L {q  w(x) \\, dx}$$\n",
    "\n",
    "On the other hand, the modal properties are evaluated with a continuous version of the same formula\n",
    "presented on [Class 11](https://nbviewer.jupyter.org/github/mmaiarocha/PEC00025/blob/master/Class_11_FreeVibrationMDOF.ipynb?flushcache=true). In particular, the modal mass and the modal load are:\n",
    "\n",
    "\\begin{align*}\n",
    " \\vec\\phi_k^{\\intercal}{\\mathbf M} \\vec\\phi_k &\\implies M_k = \\int_0^L{\\mu \\, \\varphi^2(x) \\, dx} \\\\\n",
    " \\vec\\phi_k^{\\intercal} \\vec{F}               &\\implies F_k = \\int_0^L{ q  \\, \\varphi(x)   \\, dx} \n",
    "\\end{align*}\n",
    "\n",
    "The static modal response can be calculated as:\n",
    "\n",
    "$$ u_k = F_k/K_k$$\n",
    "\n",
    "where $K_k = \\omega_k^2 M_k$. \n",
    "\n",
    "Let us apply this approach to the beam example above. \n",
    "For the assume modal shape we have the curvature:\n",
    "\n",
    "$$ \\varphi^{\\prime\\prime}(x) = -\\frac{8}{L^2}$$\n",
    "\n",
    "Hence:\n",
    "\n",
    "\\begin{align*}\n",
    " T_{\\rm ref} &= \\frac{1}{2} \\int_0^L {\\mu \\left[  \\frac{4x(L - x)}{L^2} \\right]^2 \\, dx} =  \\frac{ 4}{15 }\\mu L \\\\\n",
    " V           &= \\frac{1}{2} \\int_0^L { EI \\left[ -\\frac{8}        {L^2} \\right]^2 \\, dx} =  \\frac{32}{L^3} EI\n",
    "\\end{align*}\n",
    "\n",
    "The Rayleigh quotient results:\n",
    "\n",
    "$$ \\omega_k^2 = \\frac{V}{T_{\\rm ref}} = \\frac{32EI}{L^3} \\frac{15}{4\\mu L}  \n",
    "                                      = \\frac{ 120}{L^4} \\left( \\frac{EI}{\\mu} \\right)$$\n",
    "\n",
    "which gives a fundamental frequency comparing to the exact solution as:\n",
    "\n",
    "$$ \\omega_k \\approx \\left( \\frac{3.31}{L} \\right)^2 \\sqrt{\\frac{EI}{\\mu}} \n",
    "            \\approx \\left( \\frac{ \\pi}{L} \\right)^2 \\sqrt{\\frac{EI}{\\mu}} $$\n",
    "\n",
    "with an error of approximatelly 11%. The modal shape approximation may also be used to estimate\n",
    "the displacement at beam center, for which we calculate the modal mass and the modal load as:\n",
    "\n",
    "\\begin{align*}\n",
    " M_k &= \\int_0^L{\\mu \\, \\left[  \\frac{4x(L - x)}{L^2} \\right]^2 \\, dx} =  \\frac{8}{15}\\mu L\\\\\n",
    " F_k &= \\int_0^L{ q  \\, \\left[  \\frac{4x(L - x)}{L^2} \\right]   \\, dx} =  \\frac{2}{ 3}  q L\n",
    "\\end{align*}\n",
    "\n",
    "The modal stiffness is then:\n",
    "\n",
    "$$ K_k = \\frac{ 120}{L^4} \\left( \\frac{EI}{\\mu} \\right) \\cdot \\frac{8}{15}\\mu L = \\frac{64EI}{L^3}$$\n",
    "\n",
    "and the modal displacement is:\n",
    "\n",
    "$$ u_k = \\frac{2}{3} qL \\cdot \\frac{L^3}{64EI} = \\frac{4 q L^4}{384EI} \\approx \\frac{5 q L^4}{384EI} $$\n",
    "\n",
    "The modal displacement is already the displacement at the beam center, for $\\varphi(L/2) = 1$.\n",
    "The implied error is hence 20%, not bad for an arbitrary elastic line shape.\n",
    "The following scripts show how to numerically accomplish these calculations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAngAAAEWCAYAAAD1geqNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABK20lEQVR4nO3dd3xUVfrH8c+ThBA60pHeAgJKpKmoCFawgA2VtWFZFlfs/efq2tfesLIuoLtYsCOiqBRRkS4dAggICIKA0qQlOb8/ngFiCBBCkkkm3/frdV+ZmXvuzLlcZuaZU55jIQREREREJHbERbsCIiIiIpK3FOCJiIiIxBgFeCIiIiIxRgGeiIiISIxRgCciIiISYxTgiYiIiMQYBXgiUuyZWX0zC2aWEIXXHmRmD+Ww7Bgzu/oAnz9q5yYi0aMAT0RyxcyWmNkWM9tkZqvMbKCZlY12vURERAGeiBycs0IIZYHWQDvgH1kLqOVIRKTgKcATkYMWQvgZ+AxoCRDpErzWzBYACyKP/dXMFprZOjMbamaH7jw+Uv56M1tkZmvM7Akzi4vsizOzf5jZT2a22szeMLMKkX1JZvY/M1trZr+b2SQzqx7ZV8HM/mNmK83sZzN7yMziI/vizezJyGstAs7Y1/lFWitvM7MZZrY58rzVzewzM9toZl+Z2SGZynczs9mROo0xs8My7TvSzKZGjnsHSMq07xAzG2Zmv5rZb5HbtXNyDcysvZlNNrMNkRbVp7MUudjMlkbO+e4sx30fqetKM3vBzBJzcm0i+680s7mR+o4ws3o5qa+I5C8FeCJy0MysDnA68EOmh88GjgKam9mJwL+AC4CawE/A21me5hygLd4a2B24MvJ4r8jWGWgIlAVeiOy7HKgA1AEqA32ALZF9rwNpQGPgSOBUYOf4tb8CZ0Yebwucn4PTPA84BUgGzsID2v8DquCfpddH/i2SgbeAG4GqwHDgEzNLjAROHwH/BSoB70aed6c4YCBQD6gbOZcXyJnngOdCCOWBRsCQLPuPA5oCJwH3Zgo604GbIudxTGT/37Mcm+21MbOzI/8G50bO9ZvIuYtItIUQtGnTpu2AN2AJsAn4HQ/YXgJKRfYF4MRMZf8DPJ7pfllgB1A/U/kumfb/HRgZuT0S+HumfU0jxybggcY44IgsdasObNtZn8hjPYHRkdujgD6Z9p0aqUPCPs714kz33wdeznT/OuCjyO17gCGZ9sUBPwOdgI7ACsAy7R8HPLSX100Bfst0fwxw9V7KjgXuB6pkebx+5NxqZ3psInDRXp7nRuDDTPf3dW0+A67Kcq5/APWi/f9Tm7bivqkFT0QOxtkhhIohhHohhL+HELZk2rcs0+1D8SAQgBDCJmAtUGsv5X+KHLPHsZHbCXgQ919gBPC2ma0ws8fNrATeAlYCWBnpevwdeBWoluk5s77e/qzKdHtLNvd3TjDJeq4ZkdeqFdn3cwghZPfaZlbazF6NdEdvwIO2iju7lvfjKrx1cV6kq/rMLPt/yXT7j531NbPkSFfwL5HXfARvzctsb9emHvBcpn/jdYDx5+sqIlGgAE9E8kvmIGYFHgwAYGZl8C7VnzOVqZPpdt3IMXscG9mXBqwKIewIIdwfQmgOdMC7XS/DA5JteGtWxchWPoTQIvIcK7N5vbyS9Vwt8lo/R163VuSx7F77FryF8qjgXa0ddz7N/l40hLAghNATD2IfA96L/Dvvz8vAPKBJ5DX/L5vX29u1WQb8LdO/ccUQQqkQwrgcvK6I5CMFeCJSEN4ErjCzFDMribcSTQghLMlU5rbIJIM6wA3AO5HH3wJuMrMG5mlYHgHeCSGkmVlnMzs80sK1Ae+6TQ8hrAS+AJ4ys/KRiRqNzOyEyHMOAa43s9qRyRF35uG5DgHOMLOTIq2Jt+DB5jjgezw4vd7MEszsXKB9pmPL4a2Bv5tZJeCfOX1RM7vEzKpGWgx/jzycnoNDy+H/dpvMrBlwTTZl9nZtXgHuMrMWkTpUMLMeOa2ziOQfBXgiku9CCCPxsWnv461YjYCLshT7GJgCTAM+xcftAQzAu2LHAouBrfiYN4AawHt4gDIX+Br4X2TfZUAiMAf4LVKuZmTfv/Gu3enAVOCDvDhPgBBCKnAJ0A9Yg0/IOCuEsD2EsB2fkNArUqcLs7z2s0CpyHHjgc8P4KW7ALPNbBM+4eKiEMLWHBx3K/AXYCP+7/JONmWyvTYhhA/x1sK3I927s4CuB1BnEckn9uehICIiBc/MAt5FuDDadZE/07URKZrUgiciIiISY6Ia4JnZAPPEpbP2sr+Tma03s2mR7d5M+7qYWap54tS8HD8jIiIiUqRFtYvWzDriebTeCCG0zGZ/J+DWEMKZWR6PB+bjSUeXA5OAniGEOfldZxEREZHCLqoteCGEsXjepAPVHlgYQlgUGbT8Np5dXURERKTYKwqLgB9jZtPxvEu3hhBm40k0MyfeXI4vibQHM+sN9AZISkpqU7duXqa7koKUkZFBXJyGjRZFunZFm65f0aVrV7TNnz9/TQiham6OLewB3lR8yZtNZnY6voZjE7JP+pltX3MIoT/QH6Bp06YhNTU1n6oq+W3MmDF06tQp2tWQXNC1K9p0/YouXbuizcxysspOtgp1WB9C2BBZ0ogQwnCghJlVwVvsMmdWr83uzOoiIiIixVqhDvDMrMbOJX3MrD1e37X4pIomkcz2iXjC1KHRq6mIiIhI4RHVLlozewvoBFQxs+X4sjwlAEIIrwDnA9eYWRq+fM9FkUW608ysL56JPh4YEBmbJyIiIlLsRTXAiyyMva/9LwAv7GXfcGB4ftRLRERECtaOHTtYvnw5W7fmZIW92JKUlETt2rUpUaJEnj1nYZ9kISIiIsXA8uXLKVeuHPXr1ycyOqtYCCGwdu1ali9fToMGDfLseQv1GDwREREpHrZu3UrlypWLVXAHYGZUrlw5z1suFeCJiIhIoVDcgrud8uO8FeCJiIiIxBgFeCIiIlLsrV27lpSUFFJSUqhRowa1atXadX/79u0ADB06lEcffXSfzzNo0CD69u1bEFXeJ02yEBERkWKvcuXKTJs2DYD77ruPsmXLcuutt+7an5aWRrdu3ejWrVuUanhgFOCJiIiIZKNXr15UqlSJH374gdatW3P44YczefJkXnjhBT755BMeeughtm/fTuXKlRk8eDDVq1ePdpV3UYAnIiIihcqNN0KkMS3PpKTAs88e+HHz58/nq6++Ij4+nkGDBu16/LjjjmP8+PGYGa+99hqPP/44Tz31VF5V96ApwBMRERHZix49ehAfH7/H48uXL+fCCy9k5cqVbN++PU9z2OUFBXgiIiJSqOSmpS2/lClTJtvHr7vuOm6++Wa6devGmDFjuO+++wq2YvuhWbQiIiIiB2j9+vXUqlULgNdffz3KtdmTAjwRERGRA3TffffRo0cPjj/+eKpUqRLt6uxBXbQiIiIimeytu7VXr1706tULgO7du9O9e/d9lokmteCJiIiIxBgFeCIiIiIxRgGeiIiISIxRgCciIiISYxTgiYiIiMSYqAZ4ZjbAzFab2ay97L/YzGZEtnFm1irTviVmNtPMppnZ5IKrtYiIiEjhFu0WvEFAl33sXwycEEI4AngQ6J9lf+cQQkoIoW0+1U9ERESKifj4eFJSUnZtjz76aJ4997Rp0xg+fHiePd/+RDUPXghhrJnV38f+cZnujgdq53ulREREpFgqVaoU06ZNy5fnnjZtGpMnT+b000/Pl+fPKtoteAfiKuCzTPcD8IWZTTGz3lGqk4iIiMSw9evX07RpU1JTUwHo2bMn//73vwG45ppraNu2LS1atOCf//znrmMmTZpEhw4daNWqFe3bt2f9+vXce++9vPPOO6SkpPDOO+/ke70thJDvL7LPCngL3rAQQst9lOkMvAQcF0JYG3ns0BDCCjOrBnwJXBdCGJvNsb2B3gBVq1ZtM2TIkHw4CykImzZtomzZstGuhuSCrl3RputXdBWla1ehQgUaN24MQMk77iBu5sw8ff6Mww9n22OP7bNMxYoVadGixa77N998M+eddx6jRo3i4Ycf5pprrmHw4MF8+OGHAKxbt45KlSqRnp7OWWedxeOPP05ycjJt27Zl4MCBtGnThg0bNlC6dGneeecdpk6dylNPPZXtay9cuJD169f/6bHOnTtPye0wtEK/VJmZHQG8BnTdGdwBhBBWRP6uNrMPgfbAHgFeCKE/kbF7TZs2DZ06dSqIaks+GDNmDLp+RZOuXdGm61d0FaVrN3fuXMqVK+d3EhMhPj5vXyAxkcSdz78XpUqVYsaMGXs83r17dz799FNuvfVWpk+fvquegwcPpn///qSlpbFy5Up++uknypYty6GHHrrr331n2aSkJBITE3efYxZJSUkceeSRB3GCf1aoAzwzqwt8AFwaQpif6fEyQFwIYWPk9qnAA1GqpoiIiOSlZ5+Ndg3+JCMjg7lz51KqVCnWrVtH7dq1Wbx4MU8++SSTJk3ikEMOoVevXmzdupUQAmYW7SpHPU3KW8D3QFMzW25mV5lZHzPrEylyL1AZeClLOpTqwLdmNh2YCHwaQvi8wE9AREREYt4zzzzDYYcdxltvvcWVV17Jjh072LBhA2XKlKFChQqsWrWKzz7zaQLNmjVjxYoVTJo0CYCNGzeSlpZGuXLl2LhxY4HVOdqzaHvuZ//VwNXZPL4IaLXnESIiIiK5s2XLFlJSUnbd79KlC1deeSWvvfYaEydOpFy5cnTs2JGHHnqI+++/nyOPPJIWLVrQsGFDjj32WAASExN55513uO6669iyZQulSpXiq6++onPnzjz66KOkpKRw1113ceGFF+bruRTqLloRERGRgpKenp7t43Pnzt11++mnn951e9CgQdmWb9euHePHj9/j8Z2tegWhKKVJEREREZEcUIAnIiIiEmMU4ImIiEihEO3cvNGSH+etAE9ERESiLikpibVr1xa7IC+EwNq1a0lKSsrT59UkCxEREYm62rVrs3z5cn799ddoV6XAJSUlUbt27Tx9TgV4IiIiEnUlSpSgQYMG0a5GzFAXrYiIiEiMUYAnIiIiEmMU4ImIiIjEGAV4IiIiIjFGAZ6IiIhIjFGAJyIiIhJjFOCJiIiIxBgFeCIiIiIxRgGeiIiISIxRgCciIiISYxTgiYiIiMQYBXgiIiIiMSaqAZ6ZDTCz1WY2ay/7zcyeN7OFZjbDzFpn2tfFzFIj++4suFqLiIiIFG7RbsEbBHTZx/6uQJPI1ht4GcDM4oEXI/ubAz3NrHm+1lREBAghZ5uISDQlRPPFQwhjzaz+Pop0B94IIQRgvJlVNLOaQH1gYQhhEYCZvR0pOyefqywiRcSWLbBuHaxbG1i/dD1bflrNtmWrSf91HeH39bBhPXEb1hO/aT0Jf2wgcct6Sm5bT9L2jZTI2EqJjG0kZmwlMWMbiWEriWEbSWwliW37fe0MjK2UZCtJbLeSbLcktseVZHtcEjviSrIjPoltJcqxLakC20tVIK1MBdLLViCjXAUoXx6rWIGE6pUpWacapetXo0LtclSqbFSqBCVLFsA/nogUeVEN8HKgFrAs0/3lkceye/yo7J7AzHrjrX9UrVqVMWPG5EtFJf9t2rRJ16+IyotrFwJs3JjAunWJvq0twfafNxG/Yg1Jq1dRet0qKm74hfJb11Jx+xqqhtVUYzVNWU0iO/b6vFtIYnNcOTYllGdzQjm2JpRlU0JF0uMTSYsvQVqJkqTHJ5JeIpH0hBKkl0gEs33W1TIyiE/bQdyO7SSkbSc+siWk+5aYtpWK21dQZvM8yqRvoFzGBkqyfa/Pt5WSrKYaM6jOmriq/F6iChtKVWF9+Rr8Ubk626pXI6NWFUrWLE2lyjuoVGk7lSptp0yZtP1VNUf03iu6dO2Kr8Ie4GX30RT28fieD4bQH+gP0LRp09CpU6c8q5wUrDFjxqDrVzTl5Npt3QrLlsHSpbu3Xxf8Tpi/gKRlC6iwZiF10hZTl6Ucy1Jqs5xSbP3Tc2yPK8mG0jX445BqbK9wKJsrp7CoajXiD61OYu1qJNWpSum6lSlVvQIJlb21rFTJkpQCquTf6efM1q2wfj071qxnyy/r+WP5OrYuXc325avIWLka+3U1h6xbTY31qym9aRbl1/9Cid93wNLdT7GZ0iyjDkupyzLqsLxEQ9ZXa8K2uk2IS25M9cblqFuXXVvt2pCYuP+q6b1XdOnaFV+FPcBbDtTJdL82sAJI3MvjIlKIbd4MP/4ICxfu3pbM307GvPlUWTWbJsynCQtoxgLOZCFVWbPr2AyMzRVrsbVaXUKtNvzR4BxoWoek5LpY3TpQpw6JVapQJS+arKIhKQmSkihRvTolWkD5/ZXPyIDVq2HZMsJPS/lj/jK2zl9GpcVLqfLzMjqsHk7Zjb/Az/j2PfxCdRbQhAU0YQRNWGjJrDu0JQlNG9GoaQKNG7Nra9jQqyQiRVNhD/CGAn0jY+yOAtaHEFaa2a9AEzNrgH90XQT8JYr1FJFMfv0VZs+GOXN8GzeuFatXplPql0W0YDYtmUVLZnFW/CyaZKSSENJ2Hbu1Wh1CoyYkNj8XmjWBJr7FNWxIuaQkykXxvAqVuDioUQNq1MDataMMUCZrmZ0R9YIFsGABVeYtpPycBbRf9Dkl1w70fo+fYduKkswf04zpGS0ZR0v605LZtCSjdl2aHhZH+fKNWbAAmjf37ZBDCv50ReTARDXAM7O3gE5AFTNbDvwTKAEQQngFGA6cDiwE/gCuiOxLM7O+wAggHhgQQphd4CcgUsxt3AgzZsC0aTBr1u6A7vc1O2jBbNowhWNKTKF3/ESa7phDSbbsOjajXgPijmgJLbtBy5bQogUkJ5NUqlT0TijWlCkDRxzhG/6Bv+tDf9MmSE2F2bMpOWsWh8+aRfMZY4n/efCuw7f8UpbU31rx7dY2jH2/Lc/QhlSaUq1G/K5gr2VLSEmBww+H0qUL+gRFZG8sFKP5/E2bNg2pqanRrobkksaSRE8IsHIl/PCDB3M7t4ULIY50WjKLjqUmcWL5KaRkTKHObzNISPPZpqFcOdY3bEjFzp09CmjZ0iODsmWjeUqyN+vXe/PrrFkwcyZMnUr61KnEb/XxjtsTy7CkYgo/WBtG/t6Gb7a1I5WmWFwcTZt6sJeSAkce6X+rVo3myYg+N4s2M5sSQmibm2MLexetiETBmjUwaRJMnOjbpEne7QpQnvWcW3M8D1ccR9tG46izYgIltmyELUCJ8tC6NbTpC23aQJs2WOPGTBs7Vl8yRUWFCtChg28R34wcSacaNWDKFBInTyZ5yhSSp73GhdueB2B72UNYUuMYJoYODBvZgfvfas8fkQ7jQw+Fdu2gfXvf2rXzlxCR/KUAT6SY27IFpk7dHcxNnAiLFvk+Mzi58RIeafo1RzcZR4OV4yi9ZDa2MsCqOO/6u+JSDwbat4dGjXxsmMSW+HjvQm/RAi67zB9LT4d582DiRBLHjSN53DiS5wznEiDEx7OxYSt+rNaBcXTg7ZkncPfHh+56umbNdgd87dtDq1Y5m80rIjmnAE+kmFm7Fr77Dr791rfJk2FHJE1cnTrQteUyurUdTesNY6g+dzRxC5bAArzZ5Zhj4IoLdgd05TTlodjKHPRdcYU/9ttvMH48Nm4c5ceN48gJAzly8wtcC6Q3TmZl085MLNuZj37rxIgR1XnjDT8sKQmOOgqOO863Dh2g/H6nEYvIvijAE4lxS5bAN994MPfNNzB3rj9eooR3l93z11/oWnIUzVeNpvSE0fDZj16gUiXo1AluvQVOOMG/yNU6J/tyyCHQtatvAGlpMH06jBlD/OjR1B77JrU3vsq5QGjenE2ndmJm1c58urkzX0ypzKOPesNgXKRxeGfAd/zx3tUrIjmnAE8kxqxYAaNHw6hRvi1Z4o9XqADHHguX99xO1/LfcdjSzykxcgS8NN0LVKzogVzfvrBzQoQCOjkYCQm7xmJyyy0e8E2dCmPGYKNHU+6D1+mw+SU6mPFwu3Zsv60LM2qexvA17fnm+wQGDoQXXvCnSk6Gk06CE0/03x1Vop6ZWqRwU4AnUsStWwdjxuwO6Ha20FWs6HHazTfDKQ1/JHnR58R9OQIeG+X50UqU8IjvX/+CU07xKY/x8VE8E4l5CQm7B97dfruPDZg8Gb74AkaMIPHxh2ib8QBtK1aEk08m7anTmF3rNL5KrcOoUfDf/8LLL/tTpaR4sHfiid7Cpy5dkT9TgCdSxKSl+USIzz/3bfJkT2NSpox/0V15JZx4/A5abfqO+E+HwvOfeD4T8OUJLr8cTjvNoz+NoZNoKlHCx3Uecwz885/+a2XkSP+PPWIECe+9RyugVYsW3HLWWaTd2Y1JcUcxakwcI0fCiy/C00/775JjjoEuXXw78kg1PosowBMpAn7+GUaM8O+9L7+E33/3L7Cjj/bvxZNPhnbJ60kc9TkMHQoPD/dCJUt6E8cNN/g3X+PG0T4Vkb2rVAl69PAtBM/HN2IEfPopPPEECY8+yjHVqnHMmWdy9w3d2PLuKXw/vTRffeXF/vEP36pV898wXbp447Ry8UlxpABPpBBKT4cJEzxW++wzXy0CfKD5uef6F9fJJ8Mhm5bBRx/BfUO9nzYtzQcnnX02dOvm325KKCxFkZknxW7Z0sfv/fab/8IZOhTeew8GDKBUUhInnnwyJ3brxiOfdWdVqMYXX3ix4cO9S9cM2rb1eR/dunmaxqK6XLHIgVCAJ1JIbN7srXNDh8KwYZ5YOCHBZxE+9pgHdYcfDvbTEnj/fej6rkeB4InFbr7Zv8GOPlpj6ST2HHII9Ozp2/btPiV86NDdb5g+fajesSOX9ujBpU+eQ/obNZk61YO9zz6Dhx6CBx6AWrX8bdKtm49SKFky2icmkj8U4IlE0cqV/t308cfw1VewbZvPdj39dP8C6to1kvV/0SJvtbj6PV9WArwp4pFH4LzzfIqhSHGRmOhTak86CZ591pdUe/99ePdduPZa6NuX+OOOo12PHrS78lzuuacWv/7qPb1Dh8Lrr/tkjbJl/YdTt27+nqtcOdonJpJ3FOCJFLClS/276L33YNw4f6x+ffjb3/yLpmNHH3vOTz/By2/5l9bUqV6wXTtvzjv/fJ8wIVLcmXnSvCOOgPvvhzlz/M317rtw/fW+dehA1QsuoNeFF9KrVw22bvW5HEOHwiefePH4eE+/cv75cM45UL16tE9M5OBonpFIAVi8GJ580ntP69Xz3tTNm73LaMYMb6B77jk4qdUaSrz2svfL1q8Pd93l0d6TT/qTTJzo6SUU3Ilkr3lzuPdeb9WbO9f7Zjdvhhtv9P7ZU08l6e1BnHH8Bl59FZYv3/22WroUrrnGx7p27uyzdFesiPYJieSOAjyRfLJokTe2tW3r8dhtt3nar0cegfnzYdo0uOceOLzhZuytN+HMM6FmTfj7330G7COP+JOMH++DzOvXj/IZiRQxzZrB3Xf7m23OHPi//4Mff/Sl1apVgx49iBv6Ee2O2MYjj0Bqqv/g+sc/YNUqz/ldu7anH3ruOZ/NLlJUKMATyUOrVkG/ft5S16gR3HmnpzN5/HH/XpkyxRvlmjRM9+Sul1zifUEXX+zfLDff7Es7zZzpBRs0iPYpicSGww6DBx/0nJDffw+9e8PYsd4fW6MG/PWv2HffcnjLsKund/ZsuO8+WL/eGwDr1PGsQ6+95pN6RQozBXgiB2nDBnjjDc+7deihPuRn61ZvvdvZq3rbbZFe1fnzvRWhXj0/YPhwD/K+/trXFHvsMR9LpDwOIvnDzH+BPf+8N8l9/jmcdRa8/bY31SUnw8MPw9Klu3p7Z8yAefM85+Ty5fDXv3pMeM45PtRvy5Zon5TInhTgieTC9u0+8/WCC7wB7vLLPXa7806YNct7hG6/PdKrumGD/+Q/9lho2tSDuJQU/2ZYuRJeecVnVij1vkjBSkjwH1pvvAG//OLTa2vX9j7a+vU9j+Sbb8Iff9C0qQd4qak+kf3aaz1L0c7PgF69vFE+PT3aJyXiNItW5ABMmwaDBsHgwbBmjWfIv/pq+MtfvFFgV8NbRoYnHh4wAD74wH/iH3aY99VecomPtRORwqNMGbjsMt8WL/agb9AgHz5RvjxceCFcdRXWvj1t2xpt28ITT3jj+5tv+kzc11/3eRyXXeY/+po2jfZJSXEW1SYDM+tiZqlmttDM7sxm/21mNi2yzTKzdDOrFNm3xMxmRvZNLvjaS3Hx668+wDolxde4fPlln2H36ac+w65fP18H0wxYvdqDuKZNPUfXsGH+ST9hgg/oue02BXcihV2DBt5c9+OPMHq0rwwzeLD/iktJ8em1v/9OfPzuMXmrVnmjfEqKfwQ0awYdOsC//+1j+EQKWtQCPDOLB14EugLNgZ5m1jxzmRDCEyGElBBCCnAX8HUIYV2mIp0j+9sWVL2leNixw3NknXOOj6u78UbPVvLii96rOmSIJ0ZNSMBb60aO9F/4tWvDHXd4EPff/3rhl1+G9u01rk6kqImL8+R4r7/uXbivvupv+r59/YPhyit9lnsIlCzpOfSGDYNly7x1b/16n8tRo4Y3BH71lX9ciBSEaLbgtQcWhhAWhRC2A28D3fdRvifwVoHUTIqtRYt8DkSdOtC9u0+2u/FGn9Q6aZJnMKlUKVJ41SofT5ec7AvDfvWVf/DPmeOz8y65BEqViubpiEheKVfOo7UpU2DyZLj0Um+yO+YYaNUKXnjB0xvhv+9uvdXH406c6FlZhg/3IX0NG/ocDuXXk/xmIYTovLDZ+UCXEMLVkfuXAkeFEPpmU7Y0sBxovLMFz8wWA78BAXg1hNB/L6/TG+gNULVq1TZDhgzJj9ORArBp0ybKli2b58+blmZ8911lhg07lMmTKxEXFzj66LWcccZK2rdfR0JCpvdICJSfM4daH35I1a+/Ji4tjd+POIIVZ53Fmo4dyUhMzPP6xYL8unZSMHT9shf/xx9UGzWKmsOGUT41lfSSJVl18sn8fM45bG7U6E9lt2+P49tvK/Ppp4cydeohxMUFOnRYw5lnrqRt23X5tny0rl3R1rlz5ym57qUMIURlA3oAr2W6fynQby9lLwQ+yfLYoZG/1YDpQMf9vWZycnKQomv06NF5+nwLF4Zw550hVK8eAoRQp04I998fwrJl2RTesiWEgQNDaNPGC5cvH8INN4Qwd26e1ilW5fW1k4Kl65cDU6aEcPXVIZQq5Z8RHTuG8O67IWzfvkfRBQtCuP32EKpW9aL16oXw4IMh/Pxz3ldL165oAyaHXMZZ0eyiXQ7UyXS/NrC3RuuLyNI9G0JYEfm7GvgQ7/IV2af0dB9bd+qp0LixD4Y+6iifMLF4see8ql070wFLl3rC4dq1vZ9lyxZ46SXPn/Xssz6SWkSkdWufUbF8uQ/AW7oUevTwCRsPP+wTsCIaN/bRHcuX+3jeJk18VZu6dX3c78iREKXONYkh0QzwJgFNzKyBmSXiQdzQrIXMrAJwAvBxpsfKmFm5nbeBU4FZBVJrKZLWrvVgrlEjH1s3d66vS750qeezO/10dneRhOAz5847zz+cH3/cE6COHOmDaq65BtTlISLZqVTJB+AtXOi/Jps397x6dep4/pRJk3YVTUz0GPDLL734rbfCt9/6kN4WLfy35MaNUTwXKdKiFuCFENKAvsAIYC4wJIQw28z6mFmfTEXPAb4IIWzO9Fh14Fszmw5MBD4NIXxeUHWXouOHH+Cqq3ZPbm3QwPNV7Wytq1UrU+E//vBZckcc4bkPxozxtCaLFsGHH/pjmgkrIjkRH+8rZHzxhf+i7N3bP0fat/dug//9zzOmRzRqBI8+6jNwX3/d0/Jde61/Rl1/vSdYFjkQUc2DF0IYHkJIDiE0CiE8HHnslRDCK5nKDAohXJTluEUhhFaRrcXOY0XAU5y8/bYvHNG6td++/HKfCbuzYS4hc4rvVau8f6ROHejTx/OhDBjg/SePPurLiomI5FazZp4w8+effYm033/3WbgNGnhfbaaFbZOSdjf0jR/vPQ6vvOJPcdpp8MknWi1DckZrI0nM+O03701t0AB69vQhL88845+pr7wCLVtmOWDOHF+Gom5dHyPTsaOnN5kyxcfbKcWJiOSl8uXhuuu8Re+zz7z79s47/cflDTd410ImRx3l6TSXLYMHH/Rc6d26ebD30kuwefNeXkcEBXgSA3780bsw6tTxbthmzTzZaGqq57CrWDFT4RBg1Cg44wwf5PLmm96Hm5rq3SfHH69uWBHJX3Fx0KWLD76bNg3OPdcjtsaNfXHbiRP/VLx6dR/Gt3ix90hUquTdt3XqeN5O5dST7CjAkyIpBPjuO+9ubdLEW+jOO88/K7/6yuO3uMz/u3fs8KWG2rTxJcQmT4YHHvBZFi+95E8iIlLQWrXydW8XL/ZZFl984U13xx/vM8AyLX1RooQvmDN+vH/+de7so0jq1/dhKNOnR+80pPBRgCdFSnr67uTxxx3nY+ruuguWLPGBya1aZTlgwwZ48klPH3/JJbB1qy8c+dNPPu6uSpVonIaIyJ/Vru3j8ZYt87Ely5b5GrjNmvkv2C1bdhU183Vu338fFizwif3vv+/r4J58svf+Ks2KKMCTImHbNhg2rCbNmnkPxtq1vi7ssmU+fO7QQ7Mc8Ouv3qdRt67PhG3SxPttZ83yLtmkpKich4jIPpUr52NLFi6Ed97xMSbXXOPNdI8+6gvcZtKoETz3nH8WPvYYzJvnaZ+OPNK7c9PTNeSkuFKAJ4Xahg2eM7RBA3jqqaZUqOBpTubN83Vhy5TJcsCyZf7hWK8ePPKIL/44efLucXdx+i8vIkVAQoL/mp0wwVM2HXmkd1fUq+c/Xn/99U/FDzkEbr/dszoNGuQZWHr2hMsua88rr3jnhRQv+raTQmn1av8Mq1fPP7SaN4cnn5zGpEk+1m6PdRvnz/eWuUaNvGnvwgt9luy77/q4OxGRosgMTjgBPv/cf6yecor/eK1Xz3/MLlv2p+KJiT4eb9YsnzdWvvyOXQ2Ajz22RwOgxDAFeFKoLF0KffvuboA76SSfUPbVV9Cmze97TnCdNs2DuWbNfEZsnz4+rXbgQC0jJiKxpU0b/9E6Z45/7r34ov+ovfpqH4yXSVycD+F76aWpjB7t45PvvNNHrdx11x4NgBKDFOBJobBkCfztb54loH9/+MtfPFXUe+9Bu3bZHPDdd97leuSR/sv2zjv9SZ5/3j/BRERiVbNm/iN24UL/4Bw82B+76KI9ptKaQadOMGKEp/js0sVb8urX90m7q1ZF5QykAOQowDOzODM70szOMLMTzax6fldMiodFi/zHZ5MmPm7kr3/1z6z//AeaNs1SeOcasZ06+RTaiRN9hsVPP3lzX3X9txSRYqRePV8hY8kSH8syfLhPpT3rLO/OzaJ1a5+3MWeOD3V55hkP9G66CVauLOjKS37bZ4BnZo3MrD+wEHgU6An8HfjSzMab2RVmplZAOWALFvhiEcnJviTjNdd4z+qLL2bTABcJ7FJuvNHXg12wAJ591gO7//u/LJmMRUSKmerV4V//8jEuDz4I48Z518dZZ1Eum0VsmzXz1Hvz5nmjX79+PpHt+ut9hUaJDfsLzh4C/gc0CiGcFkK4JIRwfgjhCKAbUAG4NL8rKbFj/nxfZ7FZM5/Cf911nt/z+ec9DdSfZG6xO/FESq1Y4Z9EP/7oy/qULh2NUxARKZwqVty95MXDD8O4cbTp0wfOPNMXt82iSRPv6Z0/39OEvvyyD+n7+9/3mLshRdA+A7wQQs8QwtgQ9kyZGEJYHUJ4NoTwev5VT2LFkiVw5ZVw2GGekPPmm/0z6JlnoGbNLIWzBHYsXAj9+jFh8GCfgaEcdiIie1e+vPduLF7Moquugu+/h/bt9xroNWzo+d939qy89pqPh77hBvjllyjUX/JETsfgPWhmCZnulzezgflXLYkVP//svwaTk32S6403emD3xBNQo0aWwnsJ7PjxR+jbl4zExCicgYhIEVW+PEsvucR/YT/yyH4Dvfr1fdGMhQu9p2XnJN077/Tk8lK05HT8XAIwwcyOMLNTgUnAlPyrlhR1q1fDLbf4r8B//9snUvz4Izz1FFSrls0B+wjs1GInInIQypXbvaZj5kDvjDOyDfTq1vXP7Xnz4Jxz4PHHfYzeffcpj15RkqMAL4RwF3AHMAEYBJwRQnghH+slRdRvv8Hdd3uT/7PP+gDe+fPhpZegVq1sDvj+ew/qFNiJiOSvrIHe+PEe6J19NsycuUfxxo19EtzMmXDqqXD//f7Z/thjsHlzgddeDlBOu2g7As8BDwBjgBfMLOvqn1KMbdnib/oGDfxz46yzfCr+wIH+2B6mT/dCHTrA7NkeDSqwExHJf5kDvQcfZFcm5Isv9h/aWbRo4TlJJ0+Go4/2LttGjXxSxo4dBV99yZmcdtE+CfQIIfwrhPAXoD8wKv+qJUVFWprnrGvSxN/0xx/vi0u89VY2eewAUlO9WS8lBb791qPBnbNiFdiJiBSccuV2z7q94w5f26xZM09Ims002jZt4NNP/aM7OdnHV+8M/vaciinRltMA75gQwpydd0IIHwDHHuyLm1kXM0s1s4Vmdmc2+zuZ2XozmxbZ7s3psZK/QoChQ/1H39VXe4qTMWPgk0/8sT389JNPo23eHIYN837cxYv9V2TZsgVdfRER2alSJc+jt2iRR22vv+79szfemO1SF8ceC19/7Z/3iYnQo4e37H39dcFXXfZuf4mOLzGzuBBCetZ9IYS1kUTIx+Xmhc0sHngR6Ao0B3qaWfNsin4TQkiJbA8c4LGSD8aNg44doXt3b8F7/30fSnfCCdkU/uUXT3bXpIkvp3P99f4h8tBDSlAsIlKY1KjhSUkXLPDEeP36+aC7u+/2AdaZmPlk3OnTYcAAWLHC58mdcUa2w/kkCvbXglcZ+MHMBpjZtWZ2gZldZmYPmNnXwONAbleyaw8sDCEsCiFsB94GuhfAsZJLqak+o+rYY32YxiuvwKxZcO65/mb/k99+8z7bhg19oEavXn7QM8/sZRqtiIgUCvXq+dibOXN8rPQjj/hn+cMP7zG7Ij7ec+fNn+/jsMeN816cXr2ULDnaLJscxn8u4K1lJ+JdsjWBLcBc4LMQwtJcv7DZ+UCXEMLVkfuXAkeFEPpmKtMJeB9YDqwAbg0hzM7JsZmeozfQG6Bq1apthgwZktsqF1vr1yfw+uv1GTr0UEqWzOCii5Zx/vnLKFUqY4+ycdu2UevDD6k7eDAJmzez+sQTWdKrF1v2WKbiwG3atImy6s4tknTtijZdv6IrL65dmYULaTBwIFXGjWNbpUosufxyfjn9dEJCwh5lN2xI4M036/LBB7WJiwtccMEyevZcRqlSe3QESg507tx5Sgihba4ODiFEZQN6AK9lun8p0C9LmfJA2cjt04EFOT02uy05OTlIzm3bFsLTT4dQsWIIcXEh9OkTwqpVeymclhbCwIEh1KkTAoTQtWsI06blaX1Gjx6dp88nBUfXrmjT9Su68vTaffttCB06+Gd8cnII778fQkZGtkWXLAnhoou8aI0aIfznP/41IQcGmBxyGWftbwze7ZG//czs+axbriLK3ZYDdTLdr4230u0SQtgQQtgUuT0cKGFmVXJyrOReCPDRRz476uabPU3S9One07pH72oIPmmiVStvp69RA0aNguHD9zLbQkREiqRjj/UptB9/7H2z550HxxwDY8fuUbRePc+mMG6cr5Bx1VXQtq1nZJGCsb8xeHMjfyfjK1dk3Q7GJKCJmTUws0TgImBo5gJmVsPMR3eZWftIfdfm5FjJnR9+8JzD55wDJUp4nPb559CyZTaFx4/3UbVnnQXbt8O778KECdC5c0FXW0RECoIZdOsGM2b4OL3ly32G3Vln+aDsLI45xoO8t96Cdev8+6V7dx+zJ/lrnwFeCOGTyN/Xs9sO5oVDCGlAX2AEHkgOCT6+ro+Z9YkUOx+YZWbTgeeBiyKtltkeezD1Ke5WrfIsJm3a+Hv0xRf9/du1azYTKFJTd/9yS031ZSpmz4bzz8+msIiIxJyEBP/SmD8fHn0UvvkGjjjCe3KyzK4w8/Sn8+Z5NpbRo72H6Kab4Pffo1P94iCnK1kkm1l/M/vCzEbt3A72xUMIw0MIySGERiGEhyOPvRJCeCVy+4UQQosQQqsQwtEhhHH7OlYO3I4d8PTTnrTyf//z9WMXLPBUSHuMn/3lF+jTx9+ZX3wBDzzgM2Ovucab+0REpHgpXdqTJC9a5GN63nzT02LdfvseqVVKlfLkCgsWeBz43HOeEH/AAMjYc86eHKScJjp+F/gB+AdwW6ZNirAvv/QfXLfc4kMrZs2CJ57IJj3dH3/4cjaNG/s78e9/99Un7rlHSYpFRMSTJT/5pLfoXXSR327c2KO47dv/VLR6dejf35c+a9zYx+cdfbSP8JG8k9MALy2E8HIIYWIIYcrOLV9rJvlm8WIfY3fqqd6C98knvvxMcnKWghkZntE8ORnuvRe6dPG8SM8/r1x2IiKyp3r1YNAgH9DdurWvhtGypc/cy5KWrXVrn7Pxv//5UL6jj/aWvV9+iUbFY8/+ZtFWMrNKwCeRRMc1dz4WeVyKkD/+8DjtsMO8h/WRR3zo3JlnZjN0bvRon/LUqxcceqiPr3jvPf+5JSIisi+tWvkXzaef+nifc87xSXlT/tw2ZAYXX+zDue+4wxc8Sk6Gp57ao+FPDtD+WvCm4DNoLwduBb6L3N+5SREQAnzwgQd2Dz7o8yNSU30Z2JIlsxROTfUZUieeCGvX+rtt/Hg4Llcr0omISHFlBqef7jP2Xn4Z5s71hoPLLttjIka5cj5XY/ZsOP54uPVWjxFHjoxS3WPA/mbRNgghNMTXe30BmA5MA/oBLfK9dnLQFi3yFrrzzvOxdWPHesy2x8ISa9b4mrEtW8KYMT7Vad48+MtfIC6nPfkiIiJZJCT4BL2FC32WxZAh3kz3j3/Axo1/KtqkiTf6DRvmLXgnn+xfQytXRqnuRVhOv7lfBw7DU5X0i9w+qDQpkr+2bfPWuhYtPKh75hlvGT/++CwFt271mRWNG/svrL/+dfebsFSpqNRdRERiUPny3niQmuqLmD/8sEd0/ftDWtqfip5xhk/8u/deeP99aNYM+vWDdK14lmM5DfCahhCuDiGMjmy9gab5WTHJva++gsMP9zdGt27eEHfjjVnSnoQA77zj/ba33+7TaGfM8Jx2mkAhIiL5pV4970qaMMEbF/72N0hJ8az6mZQqBfff74HeUUfB9ddDu3YwcWJ0ql3U5DTA+8HMjt55x8yOwsfjSSGyYoXPTj/lFI/fRozwGK5WrSwFJ0yADh28cPnyni/l00+hefOo1FtERIqh9u13T+DbutUz6592mg/Ey6RJk93fZ7/84rNt+/TZI82eZJHTAO8oYJyZLTGzJcD3wAlmNtPMZuRb7SRH0tM9c0mzZj4T/f77YeZMT4PyJytW+ODWo4+GJUt8mZmpU32Qg4iISEEz80His2d71v2JE312xXXX+dpmmYpdcIH3SN1wA/z7354k+b//3SP7ikTkNMDrAjQATohsDYDTgTOBs/KnapITM2d6Y9wNN/jfnWMWkpIyFdq61XOiJCf7T6C77vJklFde6QtGi4iIRFPJkr522YIF3mX70kveffvCC38an1e+/O4x5Y0aeZvFqad67n35sxwFeCGEn/a15XclZU9bt/oEpNatfabs4MHw2WdZ0tRlzo9y993+Lpg714O9cuWiVncREZFsVanii6FPmwZHHukteSkpPrg8k5QU+O47Lzphgo87f/xxT94vTvkviqCxY70F++GHffr43Ln+90/JimfO9K7X887z5cS++sqDvYYNo1ZvERGRHDn8cP/e+vBD2LLFB5d37+5ZHiLi4nzlzDlzvP3ijjt8EsakSVGsdyGiAK8I+f13b7k+4QT/lfLFF76SWJUqmQqtXQvXXus/b6ZN8583P/wAJ50UnUqLiIjkhhmcfbaPz/vXv2DUKJ8MeMcdsGHDrmK1a/v48w8+gNWrfZj5TTfBpk1Rq3mhoACviPjgA/9//dprnuF75kz/QbPLjh2eJKhJE3j1VQ/yFizwnzd/yo8iIiJShCQleW7W+fN9XbPHH/cx5QMG+JrpEeec4z1af/sbPPus54H99NPoVTvaFOAVcj//7P9pzzsPqlf3CUZPPAFlymQq9OWX3mJ3/fXQpg1Mn+7TaitpuWAREYkRNWvCwIH+RdiwIVx1lada+W531rYKFXx+xrff+vfkmWf6RIxME3KLDQV4hVRGhjfENW/uuR8fe8z/T7dpk6nQwoU+JuHUU33pio8/9n7bFlpFTkREYlS7dh7UDR7sifGOOw569vzT+rbHHuujk+65B956y78WP/ooelWOBgV4hdBPP3nM1qePr8s8c6YvNlGiRKTAxo0+BqF5cx+TsHOF5m7dssy0EBERiUFmPrswNdVzg330kSeDfeQRb/DAM6888IA3jlSv7r1hf/mLL71eHCjAK0RC8Fa7li192verr/okol2pTzIy4I03fOzB44/DJZf4mIQ77vD/ySIiIsVJmTKe3X/uXOjSxVOCtWwJw4fvKnLkkR7k3X8/vPuut+a9/34U61xAohrgmVkXM0s1s4Vmdmc2+y82sxmRbZyZtcq0b0lkJY1pZja5YGue95Yu9RVa+vTxNfdmzYLevTM1yP3wAxx/PFx+ua/jN3GiDzCtWTOq9RYREYm6+vU9ahsxwhP4n3GG92pFMiAnJnpD35QpPuv2/PN9ZYzVq6Nb7fwUtQDPzOKBF4GuQHOgp5llXQx1MXBCCOEI4EGgf5b9nUMIKSGEtvle4XwSgs+MbdkSxo2Dl1/2ORP16kUKrFvnM2HbtvVZsQMHesF27aJabxERkULn1FNhxgyfjTh6tDfX3XMP/PEHAEccAePHex7Zjz7y3UOGxOZyZ9FswWsPLAwhLAohbAfeBrpnLhBCGBdC2Lmc8HigdgHXMV8tX+5rK//1r7vH2vXpE2m1S0/3xfaSk72vtm9f747t1cuzO4qIiMieEhM9n1hqqjfVPfSQr+j03nsQAiVKwP/9ny/FXr8+XHihF1u1KtoVz1sWohS2mtn5QJcQwtWR+5cCR4UQ+u6l/K1As0zlFwO/AQF4NYSQtXVv53G9gd4AVatWbTNkyJA8P5cDFQJ89lkNXnqpMenpRu/eP9K9+4pdcVu5uXNp8txzlE9N5fcjjmDB9dezuVGj6Fa6ENi0aRNly5aNdjUkF3TtijZdv6JL1w4qzJhBk+efp+yPP/Jb69YsuO46/qhfH4D0dGPIkNoMHNiAUqXSuemm+XTq9Gt0K5xJ586dp+S6lzKEEJUN6AG8lun+pUC/vZTtDMwFKmd67NDI32rAdKDj/l4zOTk5RNvy5SF07RoChNCxYwgLF2bauXp1CFdd5Ttr1gxh8OAQMjKiVtfCZvTo0dGuguSSrl3RputXdOnaRezYEcKLL4ZQsWIICQkh3HxzCOvX79o9Z04I7dr51+/FF4fw22/Rq2pmwOSQyzgrmn19y4E6me7XBlZkLWRmRwCvAd1DCGt3Ph5CWBH5uxr4EO/yLdTefNP7+8eM8TzEo0dDo0ZAWpovKZac7GuP3XabNy3vscCsiIiIHLCEBB/PPn8+XHEFPPOMf+e+8QZkZHDYYZ5a77774O23fSncUaOiXemDE80AbxLQxMwamFkicBEwNHMBM6sLfABcGkKYn+nxMmZWbudt4FRgVoHV/AD99pvHahdf7KnrZsyA666LDKX79lsfgNe3r2cxnjHDU6CUKxftaouIiMSWqlWhf3/PRFG/vmemOP54+OEHSpSAf/4Tvv/es6+cdJKvabtlS7QrnTtRC/BCCGlAX2AE3v06JIQw28z6mFmfSLF7gcrAS1nSoVQHvjWz6cBE4NMQwucFfAo5MmqUz9p591148EEYOzaS127lSrj0Uv+PtW6dD/788ksfCCoiIiL5p21bz0gxYIBnqGjTBq65Btato107n4DRt6+vadumjd8vaqI6HTOEMDyEkBxCaBRCeDjy2CshhFcit68OIRwSPBXKrnQowWfetopsLXYeW5hs2+aTeE46CUqX9v9H//gHJIQd8PTT0LSpz82++25P0HjeeeqOFRERKShxcd5dO3++r+X+73/7d/OAAZROyqBfP0+rt36956d9+GEfUVVUKN9GPpg509PUPfWU/yCYOjWStm7MGEhJgVtu8Za72bN9+naZMlGusYiISDFVsaI31U2d6sudXXWVf0dPm8app/p3+vnneyNNx46+DHxRoAAvD2VkeONc27aeT2fYMHjpJSiz8RdfVqxzZ0+2OHSo79y1BpmIiIhE1RFH+DiqQYN2d9vecAOV4tfz1ls+UXLuXG+nefXVwp8cWQFeHlm2DE45xRvnunb1pcbO6JIOL7zgTb7vvuvZtOfMgbPOUnesiIhIYWPmEy9SU33lgX79vFVv8GB6XhSYORM6dPBd3bvDmjXRrvDeKcDLA2+/7YH/hAnehf/hh1B10QRo396ny7Zv7228DzwApUpFu7oiIiKyL4cc4unLJk2CunW9F+7EE6m9YQ6ff+49uiNGQKtWhTedigK8g7Bpk68c1rOnB/jTpsHV567D+vwNjjkGfvkF3nkHvvjC8+2IiIhI0dGmjedNefVVmD4dWrUi7s7bueGqTUyY4BnNTj7Zlz7bsSPalf0zBXi59MMPft3/+1+491745usMGn8z0Ltj//MfT54zbx5ccIG6Y0VERIqquDjo3du7bS+7DJ54Ag47jJSF7zFlcuCqq+Bf//J5GYsWRbuyuynAO0Ah+CoURx/tLXgjR8L9580g4cSOcOWV3lI3dapPoVWyYhERkdhQtao34Hz3HVSuDD16UOa8Lvz79gUMGeJtOikp8NZb0a6oU4B3ANas8UGVN9wAp54KM77bSKehN0Pr1h7ZDxgA33zjA/JEREQk9nToAJMnw3PPwfjx0LIlPWbcw4zxf3D44b5yVa9esHFjdKupAC+Hvv7aI/MRI+C5ZwNDLxlC5WOb+UjLq6/2AO+KKyLrj4mIiEjMSkjw5Mjz5kGPHvDQQ9Tt2oKxt33Cvff68K3WrWHKlOhVUdHIfqSl+eLDJ57oK1JMfXs+1w87FbvoQqhe3QdfvvIKVKoU7aqKiIhIQapZE/73Pxg9GkqXJv6cbtz/Qze+f3MxW7f6fMsnn/Q8uQVNAd4+LF/uS43dfz9ccdEWZp59Dy0uOtwXKe7Xz6dPH3VUtKspIiIi0dSpk6fSePxxGDWK9r2ak3rpQ5zddRu33QZnnFHwOfMU4O3F0KGe32bKFBh50zBe+74FJZ94yGfFpqb6KsTx8dGupoiIiBQGJUrAbbf5chdnnknpf93DO3MP55PrvmD0aB/m9d13BVcdBXhZbN8ON97okymOrvkTvxx9Nic+cxYkJXkT7H//CzVqRLuaIiIiUhjVqeOrV33+OQac2e80Vna8gHoJP3PCCd5lWxDLnCnAy2TpUjjhBHjpue18fMyjDFt0GGW//xIee8ybXjt1inYVRUREpCg47TRfxerBBznkm0/4dm0zXmvxDHfelkb37rBuXf6+vAK8iM8/9xkvlaePYk2tVnT7/i6sSxdvar39dkhMjHYVRUREpCgpWRL+8Q+YPRvr2JFeM27ml0Pb8PvwcbRu7UP680uxD/DS0+Gee+DKrit5Pe1ihm05ifKJ22DYMPjgA1+DTkRERCS3Gjb0uOL996kSt46x6cfy6Nqr6XbsWp5/Pn+6bIt1gLd6NXQ9JY31Dz3PwhLNOH3Le77u2OzZPuVFREREJC+Ywbnnes/grbdy4ZZBzI9ryvQb/sMF52ewfn3evlyxDfC++QaubD6ex8e043luoHTno7FZszwnSqlS0a6eiIiIxKKyZeGJJ7AffqBc+8P4D1dz0wfH85eWM5g6Ne9eJqoBnpl1MbNUM1toZndms9/M7PnI/hlm1jqnx+5NCNDvvrXMO6E3w9YeQ/Oqv8KQIT4Ir0mTvDo1ERERkb07/HDs669h4EDaVZzPx8tbM7bdLfzn2Y150mUbtQDPzOKBF4GuQHOgp5k1z1KsK9AksvUGXj6AY/eQkW682Po/9Ly/KVcygG19byFx4VxfZsQsz85NREREZL/i4qBXL0osnEfapVdyY8bTnHbTYbzQ6T3+2HxwUV40W/DaAwtDCItCCNuBt4HuWcp0B94IbjxQ0cxq5vDYPSQsWk7faVezvUEz4qb9QMl+T0K5cnl7ViIiIiIHonJlkt7oT8Z335NQowrXje3BDzW6HtRTJuRR1XKjFrAs0/3lQNZ1v7IrUyuHxwJgZr3x1j9aEc/IS/5J/BUdmb9uLYwZc1AnIAVr06ZNjNE1K5J07Yo2Xb+iS9eu6LH/PcXs50bQ7pNXDup5ohngZdcnmrU9cm9lcnKsPxhCf6A/QJOGTcJJ/73vAKoohcmYMWPopGTTRZKuXdGm61d06doVUSedxPIJN8DRtXP9FNEM8JYDdTLdrw2syGGZxBwcu4e4EsV20rCIiIgUIbWPqnVQx0cz4pkENDGzBmaWCFwEDM1SZihwWWQ27dHA+hDCyhweKyIiIlIsRa0FL4SQZmZ9gRFAPDAghDDbzPpE9r8CDAdOBxYCfwBX7OvYKJyGiIiISKETzS5aQgjD8SAu82OvZLodgGtzeqyIiIiIFOOVLERERERilQI8ERERkRijAE9EREQkxijAExEREYkxCvBEREREYowCPBEREZEYowBPREREJMYowBMRERGJMQrwRERERGKMAjwRERGRGKMAT0RERCTGKMATERERiTEK8ERERERijAI8ERERkRijAE9EREQkxijAExEREYkxCvBEREREYowCPBEREZEYE5UAz8wqmdmXZrYg8veQbMrUMbPRZjbXzGab2Q2Z9t1nZj+b2bTIdnrBnoGIiIhI4RWtFrw7gZEhhCbAyMj9rNKAW0IIhwFHA9eaWfNM+58JIaREtuH5X2URERGRoiFaAV534PXI7deBs7MWCCGsDCFMjdzeCMwFahVUBUVERESKqmgFeNVDCCvBAzmg2r4Km1l94EhgQqaH+5rZDDMbkF0Xr4iIiEhxZSGE/Hlis6+AGtnsuht4PYRQMVPZ30II2QZpZlYW+Bp4OITwQeSx6sAaIAAPAjVDCFfu5fjeQG+AqlWrthkyZEiuz0mia9OmTZQtWzba1ZBc0LUr2nT9ii5du6Ktc+fOU0IIbXNzbL4FePt8UbNUoFMIYaWZ1QTGhBCaZlOuBDAMGBFCeHovz1UfGBZCaLm/123atGlITU09uMpL1IwZM4ZOnTpFuxqSC7p2RZuuX9Gla1e0mVmuA7xoddEOBS6P3L4c+DhrATMz4D/A3KzBXSQo3OkcYFY+1VNERESkyIlWgPcocIqZLQBOidzHzA41s50zYo8FLgVOzCYdyuNmNtPMZgCdgZsKuP4iIiIihVZCNF40hLAWOCmbx1cAp0dufwvYXo6/NF8rKCIiIlKEaSULERERkRijAE9EREQkxijAExEREYkxCvBEREREYowCPBEREZEYowBPREREJMYowBMRERGJMQrwRERERGKMAjwRERGRGKMAT0RERCTGKMATERERiTEK8ERERERijAI8ERERkRijAE9EREQkxijAExEREYkxCvBEREREYowCPBEREZEYowBPREREJMYowBMRERGJMVEJ8Myskpl9aWYLIn8P2Uu5JWY208ymmdnkAz1eREREpDiKVgvencDIEEITYGTk/t50DiGkhBDa5vJ4ERERkWIlWgFed+D1yO3XgbML+HgRERGRmGUhhIJ/UbPfQwgVM93/LYSwRzermS0GfgMC8GoIof+BHB/Z1xvoHbnbEpiVV+chBa4KsCbalZBc0bUr2nT9ii5du6KtaQihXG4OTMjrmuxkZl8BNbLZdfcBPM2xIYQVZlYN+NLM5oUQxh5IPSJB4c7AcHKWrl4pQnT9ii5du6JN16/o0rUr2jLPPzhQ+RbghRBO3ts+M1tlZjVDCCvNrCawei/PsSLyd7WZfQi0B8YCOTpeREREpDiK1hi8ocDlkduXAx9nLWBmZcys3M7bwKns7l7d7/EiIiIixVW0ArxHgVPMbAFwSuQ+ZnaomQ2PlKkOfGtm04GJwKchhM/3dXwO9M+rE5Co0PUrunTtijZdv6JL165oy/X1i8okCxERERHJP1rJQkRERCTGKMATERERiTExGeCZWRczSzWzhWa2xyoX5p6P7J9hZq2jUU/ZUw6u3cWRazbDzMaZWato1FOyt7/rl6lcOzNLN7PzC7J+snc5uXZm1imydORsM/u6oOsoe5eDz84KZvaJmU2PXL8rolFP2ZOZDTCz1WaWbZ7e3MYsMRfgmVk88CLQFWgO9DSz5lmKdQWaRLbewMsFWknJVg6v3WLghBDCEcCDaABxoZHD67ez3GPAiIKtoexNTq6dmVUEXgK6hRBaAD0Kup6SvRy+964F5oQQWgGdgKfMLLFAKyp7Mwjoso/9uYpZYi7Aw3PlLQwhLAohbAfexpc2y6w78EZw44GKkXx6El37vXYhhHEhhN8id8cDtQu4jrJ3OXnvAVwHvI/yVxYmObl2fwE+CCEsBc9PWsB1lL3LyfULQDkzM6AssA5IK9hqSnYiCzis20eRXMUssRjg1QKWZbq/PPLYgZaRgneg1+Uq4LN8rZEciP1ePzOrBZwDvFKA9ZL9y8l7Lxk4xMzGmNkUM7uswGon+5OT6/cCcBiwApgJ3BBCyCiY6slBylXMkm8rWUSRZfNY1lwwOSkjBS/H18XMOuMB3nH5WiM5EDm5fs8Cd4QQ0r0hQQqJnFy7BKANcBJQCvjezMaHEObnd+Vkv3Jy/U4DpgEnAo3w5T+/CSFsyOe6ycHLVcwSiwHecqBOpvu18V8sB1pGCl6OrouZHQG8BnQNIawtoLrJ/uXk+rUF3o4Ed1WA080sLYTwUYHUUPYmp5+ba0IIm4HNZjYWaAUowIu+nFy/K4BHgye/XWhmi4Fm+EICUrjlKmaJxS7aSUATM2sQGUB6Eb60WWZDgcsiM1OOBtaHEFYWdEVlD/u9dmZWF/gAuFQtB4XOfq9fCKFBCKF+CKE+8B7wdwV3hUJOPjc/Bo43swQzKw0cBcwt4HpK9nJy/Zbira+YWXWgKbCoQGspuZWrmCXmWvBCCGlm1hefoRcPDAghzDazPpH9rwDDgdOBhcAf+C8bibIcXrt7gcrAS5FWoLQQQtto1Vl2y+H1k0IoJ9cuhDDXzD4HZgAZwGshhGzTOkjByuF770FgkJnNxLv87gghrIlapWUXM3sLn9lcxcyWA/8ESsDBxSxaqkxEREQkxsRiF62IiIhIsaYAT0RERCTGKMATERERiTEK8ERERERijAI8ERERkRijAE9EREQkxijAExEREYkxCvBERHLBzNqZ2QwzSzKzMmY228xaRrteIiKgRMciIrlmZg8BSUApYHkI4V9RrpKICKAAT0Qk1yLrfk4CtgIdQgjpUa6SiAigLloRkYNRCSgLlMNb8kRECgW14ImI5JKZDQXeBhoANUMIfaNcJRERABKiXQERkaLIzC4D0kIIb5pZPDDOzE4MIYyKdt1ERNSCJyIiIhJjNAZPREREJMYowBMRERGJMQrwRERERGKMAjwRERGRGKMAT0RERCTGKMATERERiTEK8ERERERizP8DeQwXKJY5QdUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L  =  1         # bar length (m)\n",
    "EI =  2.6       # bending stiffness (Nm2)\n",
    "mu =  0.260     # mass per unity length (kg/m)\n",
    "q  =  mu*9.81   # distributed load is self weight (N/m)\n",
    "\n",
    "# Proposed modal shape for first mode (second order polynomial)\n",
    "x  =  np.linspace(0, L, 200)\n",
    "qk =  4*x*(L - x)/L/L          # guessed modal shape\n",
    "q0 =  np.sin(np.pi*x/L)        # exact modal shape!!!\n",
    "\n",
    "plt.figure(1, figsize=(10,4))\n",
    "plt.plot(x, qk, 'b', x, q0, 'r')\n",
    "plt.legend(('Trial', 'Exact'))\n",
    "\n",
    "plt.xlim( 0.0, L  );  plt.xlabel('x');\n",
    "plt.ylim(-0.5, 1.5); plt.ylabel('phi(x)');\n",
    "plt.title('Proposed modal shape')\n",
    "plt.grid(True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The same calculation could be carried out with the correct modal frequency, with much more accurate results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Available fundamental vibration frequency:  4.97 Hz\n",
      "Modal mass (integrated over bar length):    27.6 kg\n",
      "Modal stiffness (from mass and frequency): 26880 N/m\n"
     ]
    }
   ],
   "source": [
    "wk =  ((np.pi/L)**2)*np.sqrt(EI/mu)     # exact fundamental frequency\n",
    "fk =    wk/(2*np.pi)\n",
    "Mk =    np.sum(mu*qk*qk)                # modal mass from guessed modal shape\n",
    "Kk =    wk*wk*Mk                        # improved modal stiffness\n",
    "\n",
    "print('Available fundamental vibration frequency: {0:5.2f} Hz'.format(fk))\n",
    "print('Modal mass (integrated over bar length):   {0:5.1f} kg'.format(Mk))\n",
    "print('Modal stiffness (from mass and frequency): {0:5.0f} N/m'.format(Kk))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum displacement approximation:  12.59mm\n",
      "Theoretical maximum displacement:    12.77mm\n"
     ]
    }
   ],
   "source": [
    "Fk =  np.sum(q*qk)         # modal force\n",
    "uk =  Fk/Kk                # modal displacement\n",
    "\n",
    "wp =  np.max(uk*qk)        # approximated elastic line\n",
    "w0 = (5*q*L**4)/(384*EI)   # exact elastic line\n",
    "\n",
    "print('Maximum displacement approximation: {0:6.2f}mm'.format(1000*wp))\n",
    "print('Theoretical maximum displacement:   {0:6.2f}mm'.format(1000*w0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fundamental frequency from Rayleigh quotient:  4.97 Hz\n"
     ]
    }
   ],
   "source": [
    "w    =  uk*qk\n",
    "V    =  np.sum( q*w  )/2    # potential energy calculated with external work\n",
    "Tref =  np.sum(mu*w*w)/2\n",
    "\n",
    "wk   =  np.sqrt(V/Tref)\n",
    "fk   =  wk/(2*np.pi)\n",
    "\n",
    "print('Fundamental frequency from Rayleigh quotient: {0:5.2f} Hz'.format(fk))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}