extra-FEA_material.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra Material: Direct stiffness methof for Finite Element analysis (FEA)\n",
    "\n",
    "The direct stiffness approach looks at the stiffness of each individual element and adds them to the overall stiffness of the FEA structure. Starting with a single element, show below we have the following relation between forces applied and deformation\n",
    "\n",
    "![Single link element with forces applied to both ends](../images/link_elem.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathbf{F}_{el} = \\mathbf{K}_{el}\\mathbf{u}$\n",
    "\n",
    "where \n",
    "\n",
    "$\\mathbf{F}_{el}=\\left[\\begin{array}{c}\n",
    "        F_{1x}\\\\\n",
    "        F_{1y}\\\\\n",
    "        F_{2x}\\\\\n",
    "        F_{2y}\\end{array}\\right]$\n",
    "\n",
    "$\\mathbf{u}=\\left[\\begin{array}{c}\n",
    "        u_{1x}\\\\\n",
    "        u_{1y}\\\\\n",
    "        u_{2x}\\\\\n",
    "        u_{2y}\\end{array}\\right]$, \n",
    "        \n",
    "and \n",
    "\n",
    "$\\mathbf{K}_{el}=\\frac{EA}{l}\\left[\\begin{array}{cccc}\n",
    "        \\cos^2\\theta & \\cos\\theta\\sin\\theta & -\\cos^2\\theta & -\\cos\\theta\\sin\\theta\\\\\n",
    "        \\cos\\theta\\sin\\theta & \\sin^2\\theta & -\\cos\\theta\\sin\\theta & -\\sin^2\\theta &\\\\\n",
    "        -\\cos^2\\theta & -\\cos\\theta\\sin\\theta & \\cos^2\\theta & \\cos\\theta\\sin\\theta\\\\\n",
    "        -\\cos\\theta\\sin\\theta & -\\sin^2\\theta & \\cos\\theta\\sin\\theta & \\sin^2\\theta &\\\\\n",
    "        \\end{array}\\right]$. \n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a linear-elastic, small deformation FEA, we assume that the angle $\\theta$ remains constant before/after deformation. So our stiffness matrix, $\\mathbf{K}_{el}$ is constant once the geometry of the simulation is initialized. So, for a single link element that is deformed by $\\mathbf{u}$, the force that is applied to create that deformation is given by our first equation. \n",
    "\n",
    "Consider a single aluminum bar $(E=70~GPa)$, with cross-section of $1~mm^2$, and length of 1 meter. If it is initially at an angle of $45^{o}$ pinned at one end and stretched 2 mm to along the x-direction, we have the following deformation and stiffness, \n",
    "\n",
    "$\\mathbf{u}=\\left[\\begin{array}{c}\n",
    "        0\\\\\n",
    "        0\\\\\n",
    "        2\\\\\n",
    "        0\\end{array}\\right]~mm$, \n",
    "        \n",
    "$\\mathbf{K}_{el}=70\\frac{N}{mm}\\left[\\begin{array}{cccc}\n",
    "        1/2 & 1/2 & -1/2 & -1/2\\\\\n",
    "        1/2 & 1/2 & -1/2 & -1/2 &\\\\\n",
    "        -1/2 & -1/2 & 1/2 & 1/2\\\\\n",
    "        -1/2 & -1/2 & 1/2 & 1/2 &\\\\\n",
    "        \\end{array}\\right]$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "rcParams['font.family'] = 'sans'\n",
    "rcParams['font.size'] = 16\n",
    "rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "u=np.array([0,0,2,0])\n",
    "K=70*np.array([[1/2,1/2,-1/2,-1/2],[1/2,1/2,-1/2,-1/2],[-1/2,-1/2,1/2,1/2],[-1/2,-1/2,1/2,1/2]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[F1x, F1y, F2x, F2y] =  [-70. -70.  70.  70.] N\n"
     ]
    }
   ],
   "source": [
    "F=K@u\n",
    "print('[F1x, F1y, F2x, F2y] = ',F,'N')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So to stretch the 1-meter by $1~mm^2$ link from its initial configuration to the right by 2 mm, it takes a $\\mathbf{F}=70\\hat{i}+70\\hat{j}$ force and a negative $\\mathbf{F}$ reaction force. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'forces and deformation of single Al-bar\\nscale factor = 10')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEXCAYAAABWNASkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xd4FPXWwPHvyaYBCb1IC0URUMEWAcUGoiKKeEWviHr12q5e9eqrgIoKiHgBFbkWLIjYFQVREUFs2BtYKEFAQErovZeU8/4xQ3Z22ZAFkp1Ncj7Pk+fZnfnN7NnJ5uTslDOiqhhjjClbEvwOwBhjTPGz5G6MMWWQJXdjjCmDLLkbY0wZZMndGGPKIEvuxhhTBllyP0Qi0lxEfhORrSLyH7/jKWkislhEOkU59hoR+TbKsRVE5EMR2SwiYw8tyuInIhkisk1EAjF+3Toi8rX7+Rp2AMudJiLziimGl0VkUHGsK8K6vxSR6/czX0XkiJJ47bIu0e8AyoA+wJeqerzfgZRylwB1gBqqmut3MCKyGLheVT8DUNWlQJoPodwIrAMq6wFclKKq3wDNSyyqAyAiAiwEdqnqUX7HU15Y5X7oGgFZB7OgiNg/16BGwPyDSexlfDs2AuYcSGKPQ6cDtYGmInKSHwGU8c9IRJbcD4GIfAF0AJ52v7IfKSJVRORVEVkrIktE5H4RSXDHXyMi34nIcBHZAAxwp98gIn+4X73niMgJ7vR6IvKuu66/vLt9RKSNiEwXkS0islpEHi8kxmoiMtFdx0b3cQPP/C9F5CE3rq0i8omI1PTMv8p9H+tF5L4itkcNEZngxvQzcHjY/BYi8qmIbBCReSLyd3f6g0A/4DJ3O14nIgnutlsiImvcbVrFHd/Y/bp+nYgsBb7wTPuniCxz3+tNInKSiMwUkU0i8rQnlsNF5Av3fa0TkTdEpKo77zUgA/jQjaePZ/2Jnt/NBPe9LBCRGzzrHiAi77gxbxWRLBHJ3M92O0VEpomzS2qaiJziTn8ZuBro48axz+4wEenifma2ishyEenlTj9TRLI94xaLSC93W2wWkbdFJNUzv4+IrBSRFSJyvexnd4iIXCAiv7vb9HsRaV3Ye3NdDXwATHIfH6guIrLI/T09KsG/p0J/h573fLeIzAS2S3lL8KpqP4fwA3yJ8/V97/NXcT7I6UBjYD5wnTvvGiAXuA1nl1gF4FJgOXASIMARONVaAvALTtJLBpoCi4Bz3XX9AFzlPk4D2hUSXw2gO1DRjWks8H5Y/AuBI914vgSGuPOOArbhVF4pwONu/J0Kea0xwDtAJeAY9319686rBCwD/um+9xNwdjcc7c4fALzuWde1wAL3facB44HX3HmNAXW3dSU37r3TngNSgXOAXcD7OFVjfWANcIa7jiOAs933VQv4Gvif5/UXe9+nZ/2J7vOvgGfc1zoOWAuc5Xkvu4AuQAAYDPxYyDarDmwErnK3y+Xu8xru/JeBQfv5/K0ETnMfVwNOcB+fCWSHvZ+fgXrua/4B3OTO6wysAo7G+Zy85r7XI8JjcH9va4C27nu72l13SiHxVQS2uNuiu/s7Ty7s7yfC8gpMdWPOwPl7uv4Afoe/Aw2BCn7nipjnJr8DKO0/3g+n+2HfDRzlmf8vnH3y4CT3pWHLTwFuj7DethHG3gu85D7+GngQqHmA8R4HbAyL/37P838DH7uP+wFjPPMqAXuIkNzd954DtPBM+y/B5H4Z8E3YMs8D/d3HAwhN7p8D//Y8b+6uP5Fgom3qmb93Wn3PtPXAZZ7n7wJ3FLJdLgJ+8zxfTCHJ3U0WeUC6Z/5g4GXPe/nMM+8oYGchr3sV8HPYtB+Aa9zHL7P/5L7U/YxVDpt+Jvsm9ys9zx8BnnMfjwYGe+YdQeHJ/VngobDXmof7TzNCfFfi/ONLxEnCm4C/Rfr7KWR5BTqHfT4/P4Df4bUH8vdRln5st0zxqolTZS/xTFuCUzXutSxsmYY4lXO4RkA996vvJhHZBPTFOegIcB1OtT3X/Sp/QaSARKSiiDzv7t7YgvNPoaqEnvWxyvN4B8EDh/W88arqdpyEGUktnD9g7/vzbodGQNuw93MFcFgh66vHvtsxkeD7h323JcBqz+OdEZ6nAYhIbREZ4+7K2AK8jvP7i0Y9YIOqbg2Lz/t7Dt+mqYXsFgh/n5HWtT/dcariJSLylYicvJ+xUf2eibxd92oE3BX2e2zoriOSq4F3VDVXVXfjfAOLuGvG3X21zf05rZB4lux9rSh/h/t7L2WaJffitQ6numzkmZaBs3tir/ADY8sI2zftmf6Xqlb1/KSrahcAVf1TVS/H2eUwFBgnIpUirOcunKq3rapWxtnFAs4uoKKsxPnDdRYQqYizmyeStTi7bBp6pmWEvZ+vwt5PmqreXMj6VrDvdswlNFkfykHGwe7yrd3tciWh22R/614BVBeR9LD4lhcyfn/C3+cBrUtVp6lqN5zPwfs4u8UO1Eqgged5w8IG4vweHw77PVZU1bfCB4pzbKcjcKWIrBKRVThnRXURz3Edz3s52v1MpKlztk+keDJwthkU/TuEQ/uMlGqW3IuRqubh/HE9LCLpItIIuBOnoijMKKCXiJwojiPc5X4GtrgHhCqISEBEjhH3bAMRuVJEaqlqPs5XXXB2FYRLx6lYN4lIdaD/AbylccAFInKqiCQDAynkM+O+9/HAAPfbwlGEVmgTgSPFOUCb5P6cJCItC3ntt4D/E5EmIpKGs4vnbS2+0yTTcY4nbBKR+kDvsPmrcfb370NVlwHfA4NFJNU9oHgd8MZBxDEJZ7v0FJFEEbkMZzfOxKIWFJFkEblCRKqoag7Ovu1In4GivAP8U0Rauv/A++1n7AvATSLS1v28VhKR88P+0e11Fc4+8uY4uwOPw/m2mY1zbCFavcU5MaAhcDvwtju9qN9huWbJvfjdBmzHOfj5LfAmzj7NiFR1LPCwO24rTvVV3U2WXXH+IP7C+VYwCqjiLtoZyBKRbcATQA9V3RXhJf6Hc8BxHfAj8HG0b0RVs4Bb3NhW4hzoy97PIrfifNVfhbOf9iXPurbiHOTsgVN5rcL5xpFSyLpG4xzY+xrn/e/C2bbF5UGcg4ObgY9w/jF5DQbud3c99Iqw/OU4++FXAO/hHDv49ECDUNX1wAU437DW41w3cYGqrotyFVcBi93dEjfhVK8HGsNk4EmcA5cLcPb5g3P8KHzsdOAG4Gmcz8MCnGNJkVwNPKOqq7w/OAe9D+SsmQ9wTi74Hed39aI7vajfYbkm7oEHY4wBwP02NRvnDBjfLygzB8cqd2MMIvI3dzdPNZxvVB9aYi/dLLkbY8A5nXItzplbeUBhB7pNKWG7ZYwxpgyyyt0YY8ogS+4mbonTo2V/p5Hub9ly1YrZmHCW3E1ZtbcVc7qqPnmwK5ES7GV+oETkVnGaxe0Wp6lY+PyzRGSuiOwQkanu9RKmnLLkbsqqg27FXJykeG/usQIYRITrJtwrPscDD+A02ZpO8GIfUw5ZcjfFxr2adrm7K2SeiJzlTg+ISF8RWejO+8W92hAReUKcFr1b3Omn7Wf97cRpMbtJRGaIyJmFjIvUivl8dzfNFvf1BoQtc6pn3cvEac98I07/m70tdz90x7YUp1XyJnH6oVzoWc/LIvKsiEwSke1uHMVCVcer6vtE7u9zMZClqmPdi9kGAMeKSIvien1TulhyN8VCRJrjXKF6kqqmA+fidOUDpwXD5TgNrirjtPPd4c6bhnMVbnWcK2HHiqfPuGf99XGuQhzkju0FvCsitcLHqmpH4BvgVrdPyXycq4b/AVQFzgduFpGL3HVnAJOBp3AaoB0H/K6qI3FaCjzirqeriCQBHwKf4PRzuQ14w33/e/XEueo4Hecq5fD38ox4Gm+F/cwsfCvv19HADM822I5zWuPRB7k+U8pZcjfFJQ+nlcBRIpKkqotVdW+3y+tx2grPU8cM97J7VPV1VV3vdg0c5q4j0u3hrgQmqeokVc13L/WfjvMPo0iq+qWqznKXnYnTu+YMd/YVOC1631LVHDee3wtZVTucFgtDVHWPqn6B0wfG2yvlA1X9zn2tfVpCqOq/wxpveX+KuvFFYdJwLsP32ozzD8aUQ5bcTbFQ1QXAHTi7A9aI04p1bxvYwtoaIyJ3iXMXqs3itI+tQuTWu42ASyW01eypQN1o4nMbXU0V545Um3H6sOx9nULji6AesMxt2LZXUW2dY2Ebzrcir8o4/YpMOWTJ3RQbVX1TVU/FScSKcxk7FNLW2N2/fjfwd6CaqlbFqTYjtSNehnMnJm+VW0lVh0QZ3pvABKChqlbBaV6193UKa7sM+7aMXQE0FPdWb66i2jqHEJHnJNi3PPznYA8CZwHHel6jEs578v2gsvGHJXdTLNzzyjuKSApOB8edBNvPjgIeEpFmTpdYaS0iNXB2GeTi3qlHRPqxb/W51+tAVxE51z1AmyrOfUIbFDI+XDrODTZ2iUgbnP3ie70BdBKRv4vTdreGiBznzgtv/fsTzv77PuK0LT4Tp3vnmCjjQFVv8vQtD/8pdB+5G1sqzl2v9m6DvTcAeQ84RkS6u2P6ATNVdW60cZmyxZK7KS4pwBCc1sKrcA429nXnPY7TM/wTnJ7jL+K0IZ6CcyBzPs6ujV0UskvD7aHezV3nWndcb6L/DP8bGCgiW3ESX8FNLVR1Kc6++7uADTitZfdWwS/iHEfYJCLvq+oe4ELgPPe9PgP8I0ZJ9H6cf5r34ByD2OlOQ1XX4tyV6WGcVrxtcdorm3LKessYY0wZZJW7McaUQZbcjTGmDLLkbowxZZAld2OMKYMSix5SMmrWrKmNGzf26+WNMSau5OXlsWrVKtasWUN+vnONXCAQoFWrVgQCwf5zv/zyyzpV3aftRjjfknvjxo2ZPn26Xy9vjDFxYffu3Tz77LMMGjSI9etDe8INGzaM22+/PWSaiCyJZr2+JXdjjCnP8vLyePPNN3nggQdYsmTffN24cWNuuummg16/7XM3xhgf/PTTT7z55pssX7484vxBgwaRkpJy0Ou35G6MMT445ZRTmDx5MoMHD95n3rHHHsvll18eYanoWXI3xhifjB49mj59+uwzfejQoSQkHFp6tuRujDE+GD16NNdffz2qSlpaGm+88QYAHTp04Jxzzjnk9dsBVWOMKQbv/7acR6fMY8WmndSrWoHe5zbnouM9bf4fbQbb1xQ8vRa4tl86q7crCy6aTPv27XnyyScZOnQoIpG6Xh8YS+7GGHOI3v9tOfeOn8XOHKfL9fJNO7l3/CyAYIL3JHavOpWEOu3bA/DCCy/QqlWrYokpqt0yItLZveHxAhG5J8L84SLyu/sz371LjjHGlAuPTplbkNj32pmTx6NT5h3QeoorsUMUlbuIBIARwNlANjBNRCao6py9Y1T1/zzjbwOOL7YIjTEmjq3esovlm/a5VS4AKzbthNw98M2wGEcV3W6ZNsACVV0EICJjcG6aMKeQ8ZcD/YsnPGOMiU+qyru/Lmfgh4XfyfCsysvg+dNh7R8xjMwRTXKvT+jdcbJx7vKyDxFpBDQBvihk/o3AjQAZGRkHFKgxxsSL1Vt20Xf8LD6fG3k/eiq7uTt5HNfsmQxr8yOOKWnRJPdIh20Lu31TD2CcquZFmqmqI4GRAJmZmXYLKGNMqaKqvPfbcgZMyGLLrtyC6RnVK3LhsXV577cVZGyZzmMpo6ivq4MLJlUCEdizbd+VVqpdIrFGk9yzgYae5w1w7gAfSQ/glkMNyhhj4k1h1frVJzfi7vNaUDF/O712PwO/vhJa/jbtAF2fgGqNYhpvNMl9GtBMRJoAy3ESeM/wQSLSHKgG/FCsERpjjI8Kq9YbVq/AI92P5eTDa8DcSfDRnbB1ZXDB1Cpw7mA4rqdTtcdYkcldVXNF5FacO9UHgNGqmiUiA4HpqjrBHXo5MEbtjtvGmDJizZZd9H1vNp/9sTpk+j9ObsTdnVtQKWcjjP0nZI0PXbDFBXD+MEg/LIbRhorqIiZVnQRMCpvWL+z5gOILyxhj/KOqvP/7cgZMmMPmnTkF0xtWr8DQ7q05pWkNmDUWJt8NOzcEF6xUG85/DI7q5kPUoewKVWOM8VizdRd9x++nWt+1Ct68DP6cErrgsT3h3IehYvUYRls4S+7GGINTrX/w+wr6T8gKqdYbVKvAI5e05pQm1eGXl+DT/rBna3DBKg2h6//giE4+RF04S+7GmHJvzdZd3PfebD6dE1qtX9WuEfec14JK25bAK1fDkm89cwXa3ABn9YOU9NgGHAVL7saYcktVmTDDqdY37Qir1ru35pQmVeHHETD1v5DraTFQoxlc+BQ0OtmHqKNjyd0YUy6t2bqL+9+bzSdh1fqV7TK457yWpG38A0Z1h5W/B2dKANrfDmfcDUmpMY74wFhyN8aUK4VV6/WrOvvW2zdOh6+HwLfDIT94XjuHtYJuI6DusT5EfeAsuRtjyo21W3dz//uzmJIVWq1f0TaDe7u0JG3Nr/DcrbDO06o3kAJn3gOn3AaBpBhHfPAsuRtjyjxV5cOZK+n/wWw2RqrWMyrA5/fDT88R0jsg42Rn33rNZrEP+hBZcjfGlGlrt+7mgfdn83HWqpDpPdtm0LdLS9Kyv4Fn/gOblgZnJqdBpwGQeR0c4o2q/WLJ3RhTJqkqE2eupF+Ean1o99ac2iAAH98Ov70euuDhZznnrVct3W3JLbkbY8qcdducan3y7NBq/fI2GfTt0oL0vz6GEXfBNs++99Sq0HkIHNvDl0Zfxc2SuzGmTJk4cwUPvB9arderksrQS1pzWl2FCdfDnPdDFzqqG3R5DNJKpre6Hyy5G2PKhHXbdtPvg9lMmhVerTek73ktSJ/3Ljx9D+zaFJyZVsfp3tiya4yjLXmW3I0xpd5HM1fywAez2bB9T8G0elVSGdK9NafX3gXvXg4LPgtd6Pgr4ZxBUKFajKONDUvuxphSa/223TwQoVrvcVJD+nZpTuVZr8K4AaG3t6uaAV2fhMM7xDbYGLPkbowplSJV63Xdav2M6pvgrW6w1HtjOIG2N0HH+yElLfYBx5gld2NMqbJ+2276Tcjio5krQ6ZfltmQ+847gsq/PgdvD4G83cGZNZtDt6ehYZsYR+sfS+7GmFJj0qyVPPD+bNaHVeuDL27FmZVXwmvnwqqZwQUSEuHU/4PTe0Niig8R+8eSuzEm7m3YvocHPpgduVrv3ITKPw2Hb/8HmhecWfc4p1o/rFWMo40PUSV3EekMPIFzg+xRqjokwpi/AwNwGjPMUNWexRinMaacmjxrJfeHVeuHVU5lcPdWdKiwCF7qAOv/DC6QmAod+kK7WyBQfuvXIt+5iASAEcDZQDYwTUQmqOocz5hmwL1Ae1XdKCJl50oAY4wvNmzfQ78PZjMxrFr/e2YD7js7gyrfDYafRxLS6KtRe+dMmJpHxDbYOBTNv7U2wAJVXQQgImOAbsAcz5gbgBGquhFAVdcUd6DGmPLj49lOtb5uW4RqPTATRl8Jm72NvtLh7AfhxH+W2kZfxS2a5F4fWOZ5ng20DRtzJICIfIez62aAqn4cviIRuRG4ESAjo3Q35THGFL8N2/fQf0IWH85YETL90hMb8MBZdan8VX+Y8WboQkec7TT6qtIghpHGv2iSe6QOOhr2PBFoBpwJNAC+EZFjVHVTyEKqI4GRAJmZmeHrMMaUYx/PXsX9788KqdbrVE5hyMWt6ZD/A7x4CWz37BSoUB3OGwqtLi0Tjb6KWzTJPRto6HneAFgRYcyPqpoD/CUi83CS/bRiidIYU2ZtdKv1CWHV+iUnNqDfmdWp/MWd8MeHoQsd0x06D4W0WjGMtHSJJrlPA5qJSBNgOdADCD8T5n3gcuBlEamJs5tmUXEGaowpe6ZkreK+92azblvwgqM6lVMY/Ldj6Ljrc3jxXti1ObhAel04/3Fo0cWHaEuXIpO7quaKyK3AFJz96aNVNUtEBgLTVXWCO+8cEZkD5AG9VXV9SQZujCm9Nm7fw4APs/jg99BqvfsJDeh/WhqVP70JFk0NXeiEq+HsgVChagwjLb1E1Z9d35mZmTp9+nRfXtsY459PslbRN6xar52ewuCLjuKsbR/CZw9CzvbgAtUaO6c3Nj0j9sHGIRH5RVUzixpXfs/wN8bE1KYdexgwIYv3w6r1i0+oz4MnJ5H+yTWw7KfgDEmAdv92LkhKrhTbYMsAS+7GmBL36ZzV9H1vFmu3hlXr3Vpw1oa34OVHIC94lgy1WkC3EdCgyALVFMKSuzGmxGzasYcHP5zDe78tD5l+8fH1efCkHNKn9IDVs4IzEpLgtLvgtDvLXaOv4mbJ3RhTIj6bs5p7w6r1WukpDO3ajI6rR8NrT4U2+qp3gtPoq87RPkRb9lhyN8YUq807cnjwwyzGh1Xrfzu+PgOP20z6JxfD+gXBGYkVoON9zv71hECMoy27LLkbY4rNZ+6+9TXh1foFTeiY/Qy8NSp0gcanQdcnoMbhMY607LPkbow5ZJt35PDgxCzG/xparV90XD0GHbOStE8uhC3ZwRkplZ1z1k+42hp9lRBL7saYQ/LF3NXcO34Wq7cEq/WaaSk82qU+HRb/D8aNCV3gyM7OVaZV6sc40vLFkrsx5qAUVq13O7YuDx+5gLTP/gU71gVnVKwB5z3i9IWxRl8lzpK7MeaAFVatP9a5FmcueAQ+nBi6QKtLnUZflWrEONLyy5K7MSZqm3fm8NDEOYz7JTtk+oWt6zK4ye9U+vRa2O1t9FUPLhgOzTvHOFJjyd0YE5Wpc9dwz/iZYdV6MsM6VeWMef1hytehC5z4T+fuSKlVYhypAUvuxpgiFFatd2tdh8H1v6fi54MhZ0dwRrUmcOFT0OS0GEdqvCy5G2MKNXXeGu59dxartuwqmFYzLZn/dUjh1Dl3w1RPZ1dJgJNvgTP7QnJFH6I1XpbcjTH72Lwzh0ET5zA2rFq/qHUt/lvrcyp+Pgzyc4Izah8N3Z6C+ifGOFJTGEvuxpgQX85bwz1h1XqNSsk8dYZyyuw7YH5WcHBCEpzRB9rfAYnJPkRrCmPJ3RgDwJZdTrX+zvTQav1vx1Tjv9UmUmHqs6D5wRn1M51GX7VbxjhSEw1L7sYYvpq/lnvencnKzaHV+oj222k3+xZY4LklclJF6PgAtP2XNfqKY5bcjSnHtuzK4eGJf/D29GUh07sfnc7DaWNJ/frV0AWanOE0+qreJIZRmoMRVXIXkc7AEzg3yB6lqkPC5l8DPArsvQ75aVUNa/9mjIknkar16pWSeb7tOk6a1QsWem6Hl1IFzh0Ex19lrQNKiSKTu4gEgBHA2UA2ME1EJqjqnLChb6vqrSUQozGmGG3dlcPDH/3BmGmh1fplLSswMPU1Ur4fH7pA8/Ph/GFQuW4MozSHKprKvQ2wQFUXAYjIGKAbEJ7cjTFx7mu3Wl/hrdYrJjHqxKWckPVf2LE+OLhSLejyKBx1kVXrpVA0yb0+4P0Xnw20jTCuu4icDswH/k9Vl4UPEJEbgRsBMjIyDjxaY8xB2borh/9O+oO3fg79s+zZIkD/wChSpk0JXaB1D+g8GCpWj2GUpjhFk9wj/cvWsOcfAm+p6m4RuQl4Bei4z0KqI4GRAJmZmeHrMMaUgG/+XMvd40Kr9RoVE3np2D9oPWcY7N4SHFy5AXT9HzQ724dITXGKJrlnAw09zxsAK7wDVNXzXY4XgKGHHpox5lA41fpc3vp5acj0fxyZx/06guTfvgtd4KTr4az+kFo5hlGakhJNcp8GNBORJjhnw/QAenoHiEhdVV3pPr0Q+KNYozTGHJBv/1zH3e/OZPmmnQXTalRI4NWjf+WouU8hucHpVD/cafTVuL0PkZqSUmRyV9VcEbkVmIJzKuRoVc0SkYHAdFWdAPxHRC4EcoENwDUlGLMxphDbdufy8Ed/7FOtX3vETu7NfZqk2b8FJ0oATrkNzrwHkirEOFJT0kTVn13fmZmZOn369KIHGmOiEqlar1UBXm/+HUfOH4l4G33VaeU0+qp3vA+RmkMhIr+oamZR4+wKVWNKuW27cxk86Q/e+Cm0Wv/X4RvovetpEufODU4MJMMZd0P72yGQFONITSxZcjemFPtuwTr6jAut1g+rkMfrTT/j8IWvIt4T2xq0cRp91WruQ6Qm1iy5G1MKFVat/6fJCm7f8RSBhUuCE5MqQaf+ztkw1uir3LDkbkwp8/2CdfR5dybZG4PVesMKe3gtYyKNl4wLHdy0g9Poq1qjGEdp/GbJ3ZhSYvvuXAZP/oPXfwyt1ns1WsjN20cQWLIqODG1Cpw7GI7raa0DyilL7saUAt8vdPate6v1Jqk7eK3+uzRYPjl0cMuu0GUYpNeJcZQmnlhyNyaObd+dy5DJc3ntR88+dJT7G87m2q3Pk7B8Q3Bypdpw/mNwVLeYx2nijyV3Y+LUDwvX0+fdGSzbEKzWj0zdzKt13uKw1V+HDj62J5z7sDX6MgUsuRsTZ7bvzmXox3N59YdgtS7k81D9afTc+iIJq7cFB1dp6DT6OqKTD5GaeGbJ3Zg4Eqlab5W6jpdqvErN9d4rugXa3ABn9YOU9NgHauKeJXdj4sCOPbkMnTyXVzzVeoA8htb9mu5bXkPWB9v1UqOZczFSRjsfIjWlhSV3Y3z246L19Bk3k6UbdhRMy0zN5oUqr1BtY1ZwoATg1Dvg9D6QlOpDpKY0seRujE927MnlkY/n8fL3iwumJZPDsDqfcMHWt5HNucHBh7WGbiOgbuvYB2pKJUvuxvjgp0Xr6R1WrZ+auohn0kZTefOi4MBAitOS95TbrNGXOSCW3I2JoUjVekV28UStiXTa+h6yzdPoK+Nk5yYaNZvFPlBT6llyNyZGfv5rA73HzWDJ+mC1fk7qHIZXeIlKW5cHByanQacBkHkdJCTEPE5TNlhyN6aE7a3WX/lhMXvvjVOZbTxdYzynb/8YPHe844hOcMFwqJrhR6imDLHkbkwJilStd0v9lSEpL1Nh+7rgwArVoPMQaH2ZNfoyxcKSuzElYOeePB6dMo+Xvv+roFqvyWaeqf4WbXZ8Dbs9g4+6CLo8Cmm1fYnVlE1RJXcR6Qw8gXOD7FGqOqSQcZcAY4HuEPtTAAAamklEQVSTVNVukGrKpWmLN9B77AwWF1TrSs/U7+mf9BopO7YEB6bVgfOHOV0cjSlmRSZ3EQkAI4CzgWxgmohMUNU5YePSgf8AP5VEoMbEu0jVen3W8lzV12i1azp47k/N8VfCOYOc3THGlIBoKvc2wAJVXQQgImOAbsCcsHEPAY8AvYo1QmNKgemLN9B73Ez+WrcdcBp9XZ/yBX0Sx5C0K7i/naoZ0PVJOLyDT5Ga8iKa5F4fWOZ5ng209Q4QkeOBhqo6UUQKTe4iciNwI0BGhp0NYEq/nXvyeOyTeYz+LlitN5UVPFv5ZZrvng15e0cKtLsZOtwHKWl+hWvKkWiSe6RD9wVXWohIAjAcuKaoFanqSGAkQGZmphYx3Ji49suSDfQaG6zWE8nl1pSPuS0wjsDuPcGBNZs7jb4atvEpUlMeRZPcs4GGnucNgBWe5+nAMcCX4pzCdRgwQUQutIOqpizalZPHY1Pm8aKnWj9aFjMibTSNcxZAvjswIRFOvRNO7wWJKb7Fa8qnaJL7NKCZiDQBlgM9gJ57Z6rqZqDm3uci8iXQyxK7KYt+WbKB3mNnssit1lPYQ6+U97ku4UMScvKCA+se51Trh7XyKVJT3hWZ3FU1V0RuBabgnAo5WlWzRGQgMF1VJ5R0kMb4bVdOHsM+mceob4PV+okyj6cqvUi93OzgjsrEVOjQF9rdAgG7jMT4J6pPn6pOAiaFTetXyNgzDz0sY+LHL0s20nvcDBatdar1Suykb8pYesoUJNdz6KhRe+dMmJpH+BSpMUFWWhhTiF05eTz+6XxGfbOIfDeHn54wg8crvETNvDXBgcnpcPaDcOI/rdGXiRuW3I2J4NelG+k1NlitV2EbD6a8wUXylef0RqDZOU6jryoN/AnUmEJYcjfGY1dOHsM/nc8Lnmq9c8LPDEl9har5G4MDK1SH84ZCq0ut0ZeJS5bcjXH9unQjvcfOYKFbrddiI/9NeYWz5efg6Y0Ax3SHzkMhrZY/gRoTBUvuptzblZPH8M/m88LXe6t15dLAVwxIfpNKui04ML0unP84tOjiV6jGRM2SuynXflu6kd7jZrJgjZPEG8hahia/SHuZ6bkOGzjhajh7IFSo6k+gxhwgS+6mXAqv1hPI5x+BT7gn+R1SdVdwYLXGzumNTc/wLVZjDoYld1Pu/L5sE73Gziio1g+X5TyW/ALHy/xgtS4J0O7fTqOv5Ir+BWvMQbLkbsqNXTl5PPH5nzz/1ULy1Wn09a/ARO5Ieo8kb7P1Wi2d1gENMv0L1phDZMndlAsz3Gr9T7daP0YW8VjySFrI0uCghCQ47S7nJzHZp0iNKR6W3E2Ztjs3j/99FqzWU9jDHYnvcmPiRwS85zfWO8Gp1usc7V+wxhQjS+6mzAqv1tvIHzySPIrGsjI4KLECdLzP2b+eEPApUmOKnyV3U+bszs3jic/+5PmvF5GXr6Sxgz6Jb/OPxE9DBzY+Dbo+ATUO9ydQY0qQJXdTpszMdqr1+audav3MhN8ZnPQidWV9cFBKZTjnIefcdWsdYMooS+6mTNidm8eTn//Jc1851Xo1tvBA0utcHPg2dOCR58EFj0Plev4EakyMWHI3pd6s7M30GjuDeau3Asr5CT8xMOllasiW4KCKNeC8R5y+MFatm3LAkrsptXbn5vHU5wt49quF5OUrtdnIoKTRnBP4JXRgq0udRl+VavgTqDE+sORuSqXwav3vgS+5P/ENKsuO4KD0ek6v9eadfYvTGL9Ycjelyp7cfJ764k+e+dKp1hvKaoYkjqJ9ICt0YOa10GkApFbxI0xjfBdVcheRzsATODfIHqWqQ8Lm3wTcgnOPmm3Ajao6p5hjNeXc7OVOtT531VYSyOfawBR6Jb5DRdkdHFS9KVz4FDQ+1b9AjYkDRSZ3EQkAI4CzgWxgmohMCEveb6rqc+74C4HHAfsubIrFntx8nv7iT0a41XozyeaRpJEcn7AgOEgS4ORb4cx7rdGXMURXubcBFqjqIgARGQN0AwqSu6p6TkugEqGdsI05aN5qPYlcbglM4NbE90gWz41Max/ttA6of4J/gRoTZ6JJ7vWBZZ7n2UDb8EEicgtwJ5AMdIy0IhG5EbgRICMj40BjNeXIntx8np66gGemLiA3X2ktCxmaNJKWCZ6PYkISnNEH2t9hjb6MCRNNco90UvA+lbmqjgBGiEhP4H7g6ghjRgIjATIzM626NxFlrdhMr7Ez+WPlFlLZTe/EcVwfmERAPB+Z+plOtV67pX+BGhPHoknu2UBDz/MGwIr9jB8DPHsoQZnyKbxab5cwh8GJL9AkYXVwUFJF6PgAtP2XNfoyZj+iSe7TgGYi0gRYDvQAenoHiEgzVf3TfXo+8CfGHABvtZ7ODu5JfIsrEj8PHdTkDKfRV/Um/gRpTClSZHJX1VwRuRWYgnMq5GhVzRKRgcB0VZ0A3CoinYAcYCMRdskYE8me3Hye+XIBT3/hVOsdE37l4aTR1JUNwUEpVeDch+H4K611gDFRiuo8d1WdBEwKm9bP8/j2Yo7LlANzVmzhrrEz+GPlFqqzhf5Jr9It8H3ooObnw/nDoHJdf4I0ppSyK1RNzOXk5TNi6t5qPZ8LE36gf9Ir1JCtwUGVakGXR+Goi6xaN+YgWHI3MTVnxRZ6jZ3BnJVbOIz1DEoaTafAb6GDWveAzoOhYnV/gjSmDLDkbmIiJy+fZ6Yu5Kkv/iQvP4/LA1O5N/FNKsvO4KDKDaDr/6DZ2f4FakwZYcndlLg/VjrVetaKLTSSVQxJGsXJgbDWQyfdAJ36Q0q6P0EaU8ZYcjclJicvn2e/dKr1/Lxcrg98zF2JY6kge4KDqh/uXIzU6BT/AjWmDLLkbkrE3FVOtT57+Raay1KGJo/kuIRFwQESgPb/gTPuhqQK/gVqTBllyd0Uq5y8fJ77ciFPfvEnkreH/0v8gH8HPiDJ2+irTivo9hTUO96/QI0p4yy5m2LjrdaPkwUMTR5J84Ts4IBAslOpt78dAkn+BWpMOWDJ3Ryy3Lx8nvtqIU98/ieJeTu5P3Es1wY+JsHb6KthW7jwaah1pH+BGlOOWHI3h2Teqq30GjuDWcs3c3JCFkOSX6BRwprggKRKzlkwJ90ACQn+BWpMOWPJ3RyU3Lx8nv96EU989iepeVsZnPgmlydODR3UtIPT6KtaI3+CNKYcs+RuDtj81U61PjN7M2cnTGdQymjqyKbggNQqcO5gOK6ntQ4wxieW3E3UvNV6et5Gnkp6ha6BH0MHtewKXYZBeh1/gjTGAJbcTZTmr95K77EzmJG9iYsSvqN/yqtUk23BAZVqw/mPwVHd/AvSGFPAkrvZL2+1XiNvLaOTXqRj4PfQQcddAecMskZfxsQRS+6mUH8W7FvfyBWBz7k7ZQzp3kZfVTKcRl9HnOVfkMaYiCy5m33k5uXzwjd/MfzT+dTPX86Y5BdomzDXM0KgzY1wVj9ISfMtTmNM4Sy5mxB/rt5Kr3Ezmb1sPdcHJvF/yeNIlZzggBrNnEZfGe38C9IYUyRL7gbwVOufzefwvL94P/l5WiUsDg6QAJx6B5zeB5JSfYvTGBOdqJK7iHQGnsC5QfYoVR0SNv9O4HogF1gLXKuqS4o5VlNCFqzZSq+xM5mzbC23Jr7Hzckfhjb6Oqw1dBsBdVv7F6Qx5oAUmdxFJACMAM4GsoFpIjJBVb13W/gNyFTVHSJyM/AIcFlJBGyKT16+8sI3i3j80/kckzeXj5JfoFnC8uCAQAp0uBdOvg0C9iXPmNIkmr/YNsACVV0EICJjgG5AQXJXVe915z8CVxZnkKb4LVizjV5jZzB/2SruSXyHa5KnhDb6yjgZLnwKajbzL0hjzEGLJrnXB5Z5nmcDbfcz/jpgcqQZInIjcCNARkZGlCGa4pSXr4z6ZhHDPp1Pm/wZTEkeRcOEtcEByWnQaQBkXmeNvowpxaJJ7pGag2iEaYjIlUAmcEak+ao6EhgJkJmZGXEdpuQsWLON3uNmsHBpNoMS3+DvyV+FDjiiE1wwHKraP15jSrtokns20NDzvAGwInyQiHQC7gPOUNXdxROeKQ55+cqL3y7isU/m0yH/J55PeYna3kZfFapB5yHQ+jJr9GVMGRFNcp8GNBORJsByoAfQ0ztARI4Hngc6q+qafVdh/LJw7TZ6j53BsqWLGZ70Mucn/xw64Oi/wXmPQFptX+IzxpSMIpO7quaKyK3AFJxTIUerapaIDASmq+oE4FEgDRgrTuW3VFUvLMG4TRHy8pXR3/7FY5/M5fz8rxid8hpVZXtwQFodOP9xaHmBf0EaY0pMVOe3qeokYFLYtH6ex52KOS5zCBat3UbvcTNZtWQ+I5Ne5IzkmaEDjr8KznnI2R1jjCmT7OTlMiQvX3npu794bMofXKqf8ErKGNJkV3BA1Qzo+iQc3sG/II0xMWHJvYxYtHYbfcbNZMPSLF5LGslJCfML5imCtLsZOt4PyZV8jNIYEyuW3Eu5vdX68ClZXK0fcnvyeFK8jb5qNke6PQ0N2/gXpDEm5iy5l2J/rdtO77Ez2Ln0N95OGskxnkZfmpCInHonnN4LElP8C9IY4wtL7qVQfr7y0veLeWLKTP6l4/hX8kQSJT84oO5xSLcRcNgx/gVpjPGVJfdSZvG67fQeN4P8JT/yXtJIDg+sLJinialIh77Q7hZr9GVMOWcZoJTIz1de/n4xI6b8xm36Fv9I/jS00Vej9siFT0GNw/0L0hgTNyy5lwKL122nz7iZpC6dygdJL9JA1hXM0+R05OwH4cR/WqMvY0wBS+5xbG+1/vyU6fThFbonfxM6oNk5yAXDoUoDfwI0xsQtS+5xasn67fQeN5MaSyYzMeklasmWgnlaoTpy3iPQ6hJr9GWMiciSe5zJz1de/WExoz/+kXt5kfOSp4UOOKa7k9gr1fQlPmNM6WDJPY4sWe+ct95o2Xt8mPg6VWRHwTxNr4uc/zi06OJjhMaY0sKSexzYW62//vE39GMkpyfNCh1wwtXIOQ9BahVf4jPGlD6W3H22dP0O+oz9lRbL3mZC4ttUlOB9TrRqY+TCJ6FpxBtbGWNMoSy5+yQ/X3ntxyW8M/kzBspznJj0Z8E8lQSk3b+RDvdBckUfozTGlFaW3H2wdP0O7h33C8ctfZXxieNJkdyCefm1WpDQbQQ0yPQxQmNMaWfJPYby85XXf1rChMmTeIhnaZm0NDgvIYmE03uRcOqdkJjsY5TGmLLAknuMLNuwg75jf+aUZaMYE/gopNFXfr0TnGq9zlE+RmiMKUssuZew/HzljZ+W8Mnk8TzI8zRNXBWcF0gl4awHSGh3MyQEfIzSGFPWRNWMREQ6i8g8EVkgIvdEmH+6iPwqIrkicknxh1k6Lduwg+tHfgEf3cVrCQ/SNMGT2BudRsItP8Apt1piN8YUuyIrdxEJACOAs4FsYJqITFDVOZ5hS4FrgF4lEWRpk5+vvPHzUr6b9CYPyQvUT1xfMC8vKZ1A50EknHC1tQ4wxpSYaHbLtAEWqOoiABEZA3QDCpK7qi525+VHWkF5smzDDga98w3nLn+S5wLfhszLa9aZQNfhULmeT9EZY8qLaJJ7fWCZ53k20PZgXkxEbgRuBMjIyDiYVcQtVeWNH5fw6+SXeFhGUzMQbPSVk1qDpAseJXD0xVatG2NiIprkHikbaYRpRVLVkcBIgMzMzINaRzzK3riDwe9MpVv2MB4P/BIyL++YS0k6byhUquFTdMaY8iia5J4NNPQ8bwCsKJlwShdV5c2fljBv0ggGy+tUDgQbfe2pVJfkbk8QOPJcHyM0xpRX0ST3aUAzEWkCLAd6AD1LNKpSIHvjDoa9/QmXLH+EKwJZIfNyT/gnyecMhNTKPkVnjCnvikzuqporIrcCU4AAMFpVs0RkIDBdVSeIyEnAe0A1oKuIPKiqR5do5D5RVd766S+WTRrOf2UMFQJ7CubtqtyY1ItHkNj4VB8jNMaYKC9iUtVJwKSwaf08j6fh7K4p05Zv2slTYyZw2YpH6JmwoGB6Pgnkt7uF1I59rdGXMSYu2BWqUVBV3vlpEWsnD2Eg75KckFcwb2e1FlS45FkS6p/gY4TGGBPKknsRlm/ayfNvjuXyVY/QMiF4RmiuJKGn96bCaf9njb6MMXHHknshVJWxP85n+8cP0Z+JBBKCZ25uq3U8aZc+B7Vb+BihMcYUzpJ7BMs37eTlN17jitWP0ThhdcH0PQmpyFn9SDv5JusHY4yJa5bcPVSV8d/PIe+Tftwnn4W0VdtStz2VLx0B1Zv4F6AxxkTJkrtrxaadjHnjBS5fM5y6sqFg+q5AGoHO/6Vy5j+sdYAxptQo98ldVfnguxkkf9qXO+W7kGYLGzPOodolT0Lluv4FaIwxB6FcJ/eVm3bw/mtP8vd1T1NDthZM35ZYjaSuw6jW2hp9GWNKp3KZ3FWVid9MJ/3zPtwsv4ZU6+sOv5ia3YdBxer+BWiMMYeo3CX3lZu2M+XVoXRfP5J02VkwfXNyHVL/9iQ1W3b2MTpjjCke5Sa5qyqTv/6emlN7cw1ZIdX66hZXUedvgyEl3b8AjTGmGJWL5L5q4za+euVBum18iVTJKZi+LiWD9Eufpc4R1ujLGFO2lOnkrqp8MvUL6n3dm8tYWFCt55HA6mP+Rb1uAyAp1dcYjTGmJJTZ5L5q/WZ+erUvXTa9RZIEG32trNCMapePpF6GNfoyxpRdZS65qypTP59Mw2/70I1lBdX6HhJZffwdNLzgHggk+RukMcaUsDKV3FevX8/vr/Tm7M3jSZBgo6+llVpR64qRNKx3lI/RGWNM7JSJ5K6qfDPlXZr+eC/nsqagWt9BKmtOupvG590BCQn7X4kxxpQhpT65r1mzmrmv3s7p2yaHTF+Q3ob6Vz1P49pNfYrMGGP8U2qTu6ry46TXOHxaP05nY8H0LaSxpn1/juh0g7UOMMaUW1EldxHpDDyBc4PsUao6JGx+CvAqcCKwHrhMVRcXa6SPNoPta4KvCZwcNiSr6pk0+cczHFG9frG+tDHGlDZF7ogWkQAwAjgPOAq4XETCj0xeB2xU1SOA4cDQ4g7Um9jDracq804fwdF3fEBFS+zGGBNV5d4GWKCqiwBEZAzQDZjjGdMNGOA+Hgc8LSKiqkoMpN4xneZVa8XipYwxplSI5hSS+sAyz/Nsd1rEMaqaC2wGaoSvSERuFJHpIjJ97dq1BxdxBJUssRtjTIhoknuko5LhFXk0Y1DVkaqaqaqZtWrtPyGvWLGi4PGcOXP2M9IYY0y4aJJ7NtDQ87wBsKKwMSKSCFQBNnCQZsyYwX/+8x/ASewdO3Y82FUZY0y5FE1ynwY0E5EmIpIM9AAmhI2ZAFztPr4E+OJQ9rffe++9fPTRR/z888907NiR1atXs2pbfuTBlWof7MsYY0yZVeQBVVXNFZFbgSk4p0KOVtUsERkITFfVCcCLwGsisgCnYu9xsAF9+eWXTJ7sXJB06qmnkpPjtOh9oUofHnjggYNdrTHGlCtRneeuqpOASWHT+nke7wIuPdRgVJW777674PnexD5w4EBL7MYYcwDi6grV8ePH8/PPP+8z/cUXX2TLli1ce+21tGzZ0ofIjDGmdImbblo5OTn07ds34rz169eTnp5Ow4YNI843xhgTKm4q99GjRzN//vyQaUlJSdx8883cd9991K5tB06NMSZacZHct2/fzoABA0KmXXHFFQwcOJCmTa2rozHGHCiJUYeAfV9YZC2w5CAXrwmsK8ZwSorFWbwszuJVGuIsDTFCbONspKpFXpbvW3I/FCIyXVUz/Y6jKBZn8bI4i1dpiLM0xAjxGWfcHFA1xhhTfCy5G2NMGVRak/tIvwOIksVZvCzO4lUa4iwNMUIcxlkq97kbY4zZv9JauRtjjNkPS+7GGFMGxXVyF5HOIjJPRBaIyD0R5qeIyNvu/J9EpHHso4wqztNF5FcRyRWRS/yI0Y2jqDjvFJE5IjJTRD4XkUZxGudNIjJLRH4XkW8j3NPX9xg94y4RERURX06Ti2JbXiMia91t+buIXB+Pcbpj/u5+PrNE5M1Yx+jGUNT2HO7ZlvNFZJMfcQJOJ8Z4/MFpL7wQaAokAzOAo8LG/Bt4zn3cA3g7TuNsDLQGXgUuiePt2QGo6D6+OY63Z2XP4wuBj+MtRndcOvA18COQGafb8hrgaT8+kwcYZzPgN6Ca+7x2PMYZNv42nBbpvmzXeK7cC27Mrap7gL035vbqBrziPh4HnCUikW75V5KKjFNVF6vqTKCQO47ERDRxTlXVHe7TH3HuuhVr0cS5xfO0EhFu6VjCovlsAjwEPALsimVwHtHG6bdo4rwBGKGqGwFUdU2MY4QD356XA2/FJLII4jm5F9uNuUtYNHHGgwON8zpgcolGFFlUcYrILSKyECd5/idGse1VZIwicjzQUFUnxjKwMNH+zru7u+LGiYgfrVejifNI4EgR+U5EfhSRzjGLLijqvyF3l2YT4IsYxBVRPCf3YrsxdwmLhxiiEXWcInIlkAk8WqIRRRbtzdZHqOrhwN3A/SUeVaj9xigiCcBw4K6YRRRZNNvyQ6CxqrYGPiP4TTiWookzEWfXzJk4FfEoEalawnGFO5C/9R7AOFXNK8F49iuek3vMb8x9kKKJMx5EFaeIdALuAy5U1d0xis3rQLfnGOCiEo1oX0XFmA4cA3wpIouBdsAEHw6qFrktVXW95/f8AnBijGLzivZv/QNVzVHVv4B5OMk+lg7ks9kDH3fJAHF9QDURWITz1WbvwYujw8bcQugB1XfiMU7P2Jfx74BqNNvzeJwDRs3i/PfezPO4K869fOMqxrDxX+LPAdVotmVdz+O/AT/GaZydgVfcxzVxdo/UiLc43XHNgcW4F4n69ePbC0e5MbsA892Ec587bSBOVQmQCowFFgA/A03jNM6TcP7rbwfWA1lxGudnwGrgd/dnQpzG+QSQ5cY4dX+J1a8Yw8b6ktyj3JaD3W05w92WLeI0TgEeB+YAs4Ae8Rin+3wAMMSP+Lw/1n7AGGPKoHje526MMeYgWXI3xpgyyJK7McaUQZbcjTGmDLLkbowxZZAld2OMKYMsuRtjTBn0/yEA277cDQp9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# scale factor\n",
    "s=10\n",
    "x0=np.array([0,1/2**0.5])\n",
    "y0=np.array([0,1/2**0.5])\n",
    "x=x0+u[0::2]*1e-3*s\n",
    "y=y0+u[1::2]*1e-3*s\n",
    "\n",
    "plt.plot(x0,y0,'o-',label='initial position')\n",
    "plt.plot(x,y,'s-',label='final position')\n",
    "plt.quiver(x,y,F[0::2]*s,F[1::2]*s)\n",
    "plt.title('forces and deformation of single Al-bar\\nscale factor = {}'.format(s))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combining multiple elements\n",
    "\n",
    "This method can be extended to multiple, connected bars. We just have to do some extra book-keeping. Let's consider the 3-bar element that we used in [02_Gauss_elimination](../notebooks/02_Gauss_elimination.ipynb).\n",
    "\n",
    "![Triangular truss](../images/truss.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have 3 elements and 3 nodes, so there are 6 degrees of freedom, with 3 constraints. Let's use 1-m aluminum bars with $A=1~mm^2$, but set $\\alpha=\\beta=60^{o}$. \n",
    "\n",
    "Each element still has the stiffness defined above,\n",
    "\n",
    "$\\mathbf{K}_{el}=\\frac{EA}{l}\\left[\\begin{array}{cccc}\n",
    "        \\cos^2\\theta & \\cos\\theta\\sin\\theta & -\\cos^2\\theta & -\\cos\\theta\\sin\\theta\\\\\n",
    "        \\cos\\theta\\sin\\theta & \\sin^2\\theta & -\\cos\\theta\\sin\\theta & -\\sin^2\\theta &\\\\\n",
    "        -\\cos^2\\theta & -\\cos\\theta\\sin\\theta & \\cos^2\\theta & \\cos\\theta\\sin\\theta\\\\\n",
    "        -\\cos\\theta\\sin\\theta & -\\sin^2\\theta & \\cos\\theta\\sin\\theta & \\sin^2\\theta &\\\\\n",
    "        \\end{array}\\right]$, \n",
    "\n",
    "but there are two extra things to consider. The angle $\\theta$ is different for each bar, and the connections are in different locations for each bar. \n",
    "\n",
    "so if the initial geometry, $\\mathbf{r}$, is \n",
    "\n",
    "$\\mathbf{r}=\\left[\\begin{array}{c}\n",
    "        x_1\\\\\n",
    "        y_{1}\\\\\n",
    "        x_{2}\\\\\n",
    "        y_{2}\\\\\n",
    "        x_3 \\\\\n",
    "        y_3\\end{array}\\right]=\n",
    "        \\left[\\begin{array}{c}\n",
    "        0\\\\\n",
    "        0\\\\\n",
    "        0.5 \\\\\n",
    "        \\sqrt{3}/2\\\\\n",
    "        1 \\\\\n",
    "        0\\end{array}\\right]~m$\n",
    "\n",
    "then, the displacement, \n",
    "\n",
    "$\\mathbf{u}=\\left[\\begin{array}{c}\n",
    "        u_{1x}\\\\\n",
    "        u_{1y}\\\\\n",
    "        u_{2x}\\\\\n",
    "        u_{2y}\\\\\n",
    "        u_{3x}\\\\\n",
    "        u_{3y}\\end{array}\\right]=\n",
    "        \\left[\\begin{array}{c}\n",
    "        0\\\\\n",
    "        0\\\\\n",
    "        u_{2x}\\\\\n",
    "        u_{2y}\\\\\n",
    "        u_{3x}\\\\\n",
    "        0\\end{array}\\right]$, \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The forces applied at each joint must equal the sum of the terms \n",
    "\n",
    "$\\mathbf{F} = \\mathbf{K_1 u}+\\mathbf{K_2 u}+\\mathbf{K_3 u}$\n",
    "\n",
    "where, $\\mathbf{K_i}$ is the stiffness of the structure due to element $\\mathbf{i}$, as such\n",
    "\n",
    "$\\mathbf{K}_{1}=\\frac{EA}{l}\\left[\\begin{array}{cccc}\n",
    "        \\cos^2\\theta_1 & \\cos\\theta_1\\sin\\theta_1 & 0 & 0 & -\\cos^2\\theta_1 & -\\cos\\theta_1\\sin\\theta_1 \\\\\n",
    "        \\cos\\theta_1\\sin\\theta_1 & \\sin^2\\theta_1 & 0 & 0 & -\\cos\\theta_1\\sin\\theta_1 & -\\sin^2\\theta_1 \\\\\n",
    "        0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "        0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "        -\\cos^2\\theta_1 & -\\cos\\theta_1\\sin\\theta_1 & 0 & 0 & \\cos^2\\theta_1 & \\cos\\theta_1\\sin\\theta_1 \\\\\n",
    "        -\\cos\\theta_1\\sin\\theta_1 & -\\sin^2\\theta_1 & 0 & 0 & \\cos\\theta_1\\sin\\theta_1 & \\sin^2\\theta_1 \\\\\n",
    "        \\end{array}\\right]$\n",
    "\n",
    "$\\mathbf{K}_{2}=\\frac{EA}{l}\\left[\\begin{array}{cccc}\n",
    "        \\cos^2\\theta_2 & \\cos\\theta_2\\sin\\theta_2 & -\\cos^2\\theta_2 & -\\cos\\theta_2\\sin\\theta_2 & 0 & 0\\\\\n",
    "        \\cos\\theta_2\\sin\\theta_2 & \\sin^2\\theta_2 & -\\cos\\theta_2\\sin\\theta_2 & -\\sin^2\\theta_2 & 0 & 0\\\\\n",
    "        -\\cos^2\\theta_2 & -\\cos\\theta_2\\sin\\theta_2 & \\cos^2\\theta_2 & \\cos\\theta_2\\sin\\theta_2 & 0 & 0\\\\\n",
    "        -\\cos\\theta_2\\sin\\theta_2 & -\\sin^2\\theta_2 & \\cos\\theta_2\\sin\\theta_2 & \\sin^2\\theta_2  & 0 & 0\\\\\n",
    "        0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "        0 & 0 & 0 & 0 & 0 & 0 \n",
    "        \\end{array}\\right]$\n",
    "\n",
    "$\\mathbf{K}_{3}=\\frac{EA}{l}\\left[\\begin{array}{cccc}\n",
    "        0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "        0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
    "        0 & 0 &\\cos^2\\theta_3 & \\cos\\theta_3\\sin\\theta_3 & -\\cos^2\\theta_3 & -\\cos\\theta_3\\sin\\theta_3 \\\\\n",
    "        0 & 0 &\\cos\\theta_3\\sin\\theta_3 & \\sin^2\\theta_3 & -\\cos\\theta_3\\sin\\theta_3 & -\\sin^2\\theta_3 \\\\\n",
    "        0 & 0 &-\\cos^2\\theta_3 & -\\cos\\theta_3\\sin\\theta_3 & \\cos^2\\theta_3 & \\cos\\theta_3\\sin\\theta_3 \\\\\n",
    "        0 & 0 &-\\cos\\theta_3\\sin\\theta_3 & -\\sin^2\\theta_3 & \\cos\\theta_3\\sin\\theta_3 & \\sin^2\\theta_3 \\\\\n",
    "        \\end{array}\\right]$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.  0.  0. -1. -0.]\n",
      " [ 0.  0.  0.  0. -0. -0.]\n",
      " [ 0.  0.  0.  0.  0.  0.]\n",
      " [ 0.  0.  0.  0.  0.  0.]\n",
      " [-1. -0.  0.  0.  1.  0.]\n",
      " [-0. -0.  0.  0.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "Ke1 = lambda a: np.array([[np.cos(a)**2,np.cos(a)*np.sin(a)],[np.cos(a)*np.sin(a),np.sin(a)**2]])\n",
    "Ke2 = lambda a: np.array([[-np.cos(a)**2,-np.cos(a)*np.sin(a)],[-np.cos(a)*np.sin(a),-np.sin(a)**2]])\n",
    "\n",
    "K1 = np.zeros((6,6))\n",
    "K2 = np.zeros((6,6))\n",
    "K3 = np.zeros((6,6))\n",
    "\n",
    "K1[0:2,0:2]=Ke1(0)\n",
    "K1[4:6,4:6]=Ke1(0)\n",
    "K1[0:2,4:6]=Ke2(0)\n",
    "K1[4:6,0:2]=Ke2(0)\n",
    "K2[0:4,0:4]=np.block([[Ke1(np.pi/3),Ke2(np.pi/3)],[Ke2(np.pi/3),Ke1(np.pi/3)]])\n",
    "K3[2:6,2:6]=np.block([[Ke1(2*np.pi/3),Ke2(2*np.pi/3)],[Ke2(2*np.pi/3),Ke1(2*np.pi/3)]])\n",
    "print(K1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 87.5         30.31088913 -17.5        -30.31088913 -70.\n",
      "   0.        ]\n",
      "[ 30.31088913  52.5        -30.31088913 -52.5          0.\n",
      "   0.        ]\n",
      "[-1.75000000e+01 -3.03108891e+01  3.50000000e+01  1.55431223e-14\n",
      " -1.75000000e+01  3.03108891e+01]\n",
      "[-3.03108891e+01 -5.25000000e+01  1.55431223e-14  1.05000000e+02\n",
      "  3.03108891e+01 -5.25000000e+01]\n",
      "[-70.           0.         -17.5         30.31088913  87.5\n",
      " -30.31088913]\n",
      "[  0.           0.          30.31088913 -52.5        -30.31088913\n",
      "  52.5       ]\n"
     ]
    }
   ],
   "source": [
    "E=70e3\n",
    "A=1\n",
    "l=1e3\n",
    "K = E*A/l*(K1+K2+K3)\n",
    "for k in K:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ufree=np.linalg.solve(K[2:5,2:5],np.ones(3)*70)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.529210992451761"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(K[2:5,2:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAEXCAYAAABriqeHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzs3XdcVtUfwPHPeZgyHLgVt2wUFbM0LVMrd5mpZetXmqVZ2dL2tD0sM0vNtMxRqbmyoaW2NRUHiHviRBEEARnP+f1xLjyogAP0YXzfr9fzknvPufd+H0C+z73n3vNVWmuEEEIIUTrYnB2AEEIIIc6fJG4hhBCiFJHELYQQQpQikriFEEKIUkQStxBCCFGKSOIWQgghShFJ3OWIUuozpdQLxd33HPtpqJTSSinXou7rPI71slLq6wvor5VSTc+z71Cl1GGlVIpSqurFR1n8lFL/U0r96ew4hBCXxyX/YypKDq31g5eib1mnlHIDPgCu0lqvd3Y8xUUptRz4Wmv9+SXa/25gsNZ66aXYvxDllZxxlxNKKRdnx1CK1QQ8gZgL3VAZpfL/2eW4SlKSjy9ESVUq/6AIQykVopRarpRKVErFKKV652mbqpT6VCm1WCl1ErjOWjc6T5+RSqmDSqkDSqnBeS8d5+2rlOqolIpTSj2hlDpibXNvnv30UEpFKaVOKKX2KaVevoD38LRSaodSKlkptUkp1SdP2/+UUn8qpd5TSh1XSu1SSnXL095IKbXC2nYJUO0cx3oqz/u974w2D+s4e61L4p8ppSoopQKBLVa3RKXUb1b/dkqp/5RSSda/7fLsa7lS6nWl1F9AKtDYWjdaKfW3dbl9oVKqqlJquvV9+08p1TDPPoKVUkuUUglKqS1Kqf552qoqpRZY260CmhTynj2VUl8rpY5Zvyf/KaVqKqVeBzoA46x4xln9tVLqIaXUNmBbfkMd1nsZnGf5fqVUbJ6fYSul1DSgPrDQ2v/InN+jM+LbrZTqYn39slJqthXvCeB/Silbnt+RY0qpb5VSfoX9nIUo87TW8iqFL8AN2A48C7gDnYBkIMhqnwokAVdjPqB5WutGW+1dgUNAGOAFTAM00DTP9jl9OwJZwKvWcbtjElKVPO3NrOM0Bw4DN1ttDa39uhbwPvoBdaxtBwAngdpW2/+ATOB+wAUYChwAlNX+D+YStgdwjfX+vy7gOF2tuMIBb2DGGe/3Q2AB4Af4AguBN/N7D1af48BdmOGm263lqlb7cmCv9b11tb5ny62fVxOgErAJ2Ap0sfp8BUyxtvcG9gH3Wm2tgKNAmNU+C/jW6hcO7Af+LOB9P2C9Fy/rexgJVMwT5+Az+mtgifUeK+T388u7nfXz2w9cASigKdDAatsNdMmzXUcg7ozj5fYBXrZ+3jdjfh8qACOAfwF/6+c8AZjp7P9/8pKXM19yxl16XQX4AG9prTO01r8BizBJJMd8rfVfWmu71jr9jO37YxJFjNY6FXjlHMfLBF7VWmdqrRcDKUAQgNZ6udZ6o3WcDcBM4NrzeRNa6++01gesbb8BtgFt8nTZo7WepLXOBr4EagM1lVL1McniBa31Ka3175gEVZCc9xuttT6JSRKAuZyN+XDwmNY6QWudDLwB3FbAvnoA27TW07TWWVrrmcBmoFeePlOt722W1jrTWjdFa71Da50E/Ajs0Fov1VpnAd8BLa1+PYHdWusp1vZrgTnArcoMefQFXtRan9RaR1vfl4JkAlUxH1CytdZrtNYnCukP5gNLgtY67Rz9AAYD72it/9PGdq31nvPYriD/aK3nWb8PaZgPHs9preO01qcwP7dblVxGF+WY/PKXXnWAfVpre551e4C6eZb3nWP71efZF+CYlWBypGI+OKCUuhJ4C3P25445M/ruHPvD2vZu4HHMmR3WPvNe8j6U84XWOtXk2Nw+x60knGMPUK+AQ9UB1pzRN0d1zBnpGmv/YM4eC7ovoM4Z2+fs71zf+8N5vk7LZ9nH+roBcKVSKjFPuyvmqkh16+u8+y8sUU7DfE9mKaUqA19jEmFmIduc63chr3rAjgvofy5nHrsB8L1SKu/veTbmvoP9xXhcIUoNOeMuvQ4A9dTpNz7V5/Q/ZoWVfjuIufyYo6CEdz5mYC4z19NaVwI+wyS+QimlGgCTgOGYy8yVgejz2RYTfxWllHeedfXP0T/ve8zb9ygmcYZprStbr0paax/ydwCTUPK6kO/9uewDVuSJpbLW2kdrPRSIxwxbFPReTmNdIXlFax0KtMOczd99jhjzrs/5YOSVZ12tM2ItaIz9zP2fzLsf6+pB9XNssw/odsb3wlNrLUlblFuSuEuvlZg/hCOVUm5KqY6YS7WzznP7b4F7lbnBzQt4sQix+AIJWut0pVQbYOB5bueN+UMdD2Dd8BZ+Phtal2NXA68opdyVUu05/VL1mb7F3OwUar3fl/Lsy475ADFGKVXDiqWuUurGAva1GAhUSg1USrkqpQYAoZihiuKwyNr/XdbP1k0pdYVSKsQaMpgLvKyU8lJKhQL3FLQjpdR1SqlmVpI8gbl0nm01HwYaFxaI1joe84HkTqWUi3VTX95E/TnwpFIqUhlNrQ9k+e1/K+CpzM2MbsDzmKszhfkMeD1nn0qp6kqpm86xjRBlmiTuUkprnQH0BrphzhjHA3drrTef5/Y/AmOBZZibpv6xmk5dRDjDgFeVUsmYDwDfnmcMm4D3rWMfxtzg9tcFHHcgcCWQgEnEXxVyrB8xN6D9hnm/v53RZZS1/l/rjualWGP4+ezrGObM9QngGDAS6Km1PnoBsRfIGmO/ATPGfgAzXPA2jiQ3HHNZ/RDmJsIpheyuFjAbk7RjgRWYy+UAH2HGi48rpcYWso/7gacw7zUM+DtPrN8Br2OuuiQD8zA3tgG8CTxv3c3+pDW2PwyT7PdjPniedpd5Pj7CXM35xfr9+hfzMxei3Mq5O1eUc0qpEMxlao8zxrKFEEKUIHLGXY4ppfpYl5mrYM7oFkrSFkKIkk0Sd/n2AGZ8eQdm3HOoc8MRQghxLnKpXAghhChF5IxbCCGEKEUkcYtyTV1gKVAhhHA2SdxCXGZKqdeUUhuVUlkqn4Is1vPhe5RSJ5VS86SohhAiL0ncQlx+2zHPfv9wZoNSKgxTSOMuzLSeqZhn9IUQApDELUoZpdQopdR+q4TkFqVUZ2u9i1LqWeUoEbpGKVXPavtImXKjJ6z1HQrZ/1XKlN5MVEqtt2akK1Za6y+tCWGS82m+A/NY3u9a6xTgBeAWpZSvUqqJMmU+W1mx1lFKHb0UMQohSi5J3KLUUEoFYWYNu0Jr7QvciCkLCaZQye2YkqMVgfswZ6sA/wEtMDN6zQC+U0p55rP/upiz4NFW3yeBOUqpM+fTzum/yErw+b0udvrTMGB9zoLWegeQAQRaX48CplvTtk7BVCFbfpHHEkKUQlIdTJQm2ZhpP0OVUvFa69152gYDI7XWW6zlvMkv781n7yulnsdMZ7qe090JLLbKlgIsUUqtxnwYOKt0pta6Z1HeTAF8MHXU80rCzAeP1nqSUqoXZq56jZn2VghRjsgZtyg1tNbbgRGYmsxHlFKzlFJ1rOYCy0sqpZ5QSsUqpZKsUpmVOL10aI4GQL+8Z85Ae0wN8MslBXPFIK+KnH5ZfRKmGMvHVo1qIUQ5IolblCpa6xla6/aYJKsxU7VCAeUlrfHsUUB/oIpVOjSJ/EuH7gOmnVFC0ltr/VZ+sSilflRKpRTw+vEi32IMEJHnGI0xVxm2Wss+mGIpkzEVwuSOcyHKGUncotRQSgUppToppTyAdEwN7ZwSlZ8DrymlAqzyks2VUlUxl5izMFO7uiqlXuTsM9ocXwO9lFI3Wje7eSqlOiql/PPrrLXuZtXJzu/VrZD34WaNsdusmDytspsA060YOihTa/xVYK5VMQxMtaw1WuvBmPH4z87neyeEKDskcYvSxAN4C1PG9BBQA3jWavsAU070F0wJy8lABeBn4EfMGeseTMLfl9/Otdb7gJusfcZb/Z6i+P+fTMJ86LgdeM76+i4rhhjgQUwCP4L54DEMwKpD3dVqB3NDXiul1B3FHJ8QogSTucqFEEKIUkTOuIUQQohS5JyJWyn1hVLqiFIqOs86P6XUEqXUNuvfKtZ6pZQaq5TarpTakDNRhNV2j9V/m1LqnkvzdoQQQoiy7XzOuKdixtXyehr4VWsdAPxqLQN0AwKs1xDgUzCJHngJuBJoA7yUk+yFEEIIcf7Ombi11r8DCWesvgnHhBRfAjfnWf+VNv4FKiulamNmuFqitU7QWh8HlnD2hwEhhBBCnMPFzpxWU2t9EEBrfVApVcNaX5fT79iNs9YVtP4sSqkhmLN1vL29I4ODgy8yRCGEKJ/WrFlzVGud71S9ovQr7ilP85vUQhey/uyVWk8EJgK0bt1ar169uviiE0KIckAptcfZMYhL52LvKj9sXQLH+veItT4OM/VkDn/gQCHrhRBCCHEBLjZxLwBy7gy/B5ifZ/3d1t3lVwFJ1iX1n4EblFJVrJvSbrDWCSGEEOICnPNSuVJqJtARqKaUisPcHf4W8K1SahCwF+hndV+MqaS0HVNS8V4ArXWCUuo1THlFgFe11mfe8CaEEEKIcyjRM6fJGLcQQlw4pdQarXVrZ8chLg2ZOU0IIYQoRSRxCyGEEKWIJG4hhBCiFJHELYQQQpQikriFEEKIUkQStxBCCFGKSOIWQgghShFJ3EIIIUQpIolbiAu1di0kyMR/QgjnKO7qYEKUXTt2wPPPw4kT8MMPzo5GCFFOSeIW4lyOHIHXXoPPPoPsbIiKcnZEQohyTBK3EAVJTob33zevlBSz7q67ICLCuXEJIco1SdxCnCkjAyZOhFdfhfh4x3p3d7NOCCGcSBK3EDnsdvj2W3juOdi58+z2YcOgYcOL3z1yN6gQougkcQsB5iz7/vvhq6/yb/f1NQm9CH4CRgDhZ7wCALci7VkIUZ7ICYAQYC6Df/klxMTkf1Y9ahRUq1akQ3QH2gDfA68BA4AwwBtoDgwE3gAWADsxZ+hCCHEmOeMWIofWMG4c7N59+vpatWDEiAveXRKwx3rttv49kk+/TGCj9cpRA3gSeBRwv+AjCyHKMkncQoBJ2g89BJ9+apabNIEID6i3FyqnwsQrofOL0Ly/6Q7E40jMeZNzzivpIsKIwFxOvw3wLNIbEkKUVZK4RbnVevQSjqZkOFZU7AGjelAt/QSrbk5CrXgepa3RpKR9nFr4CGOAqc37sxdIu8DjuQL1gTggz1FRQG9Mwr7WWhZCiIJI4hbl1mlJO+96z4rERY2mvs48bb1HZhq3/foqz1hn3WeqADQAGlr/nvmqDRwD6lj9fYBBwMNAk6K9FSFEOSKJW4h81EuKy3d9/aQ4epN/cq7Guc+Wv8GcdT8C3AtUKp5whRDliCRuIc7Q1/Z7gQnYVsmf+UXYdydgGOBShH0IIco3SdyiXFq2Jb/7uyFM7eJ1t8kFb3jVg0U6bliRthZCCHmOW5RDBxLTePybdWetr0wyE9zH4KmssW3f2uBb5/ROG+dAVv5j40IIcTlI4hblSma2nYdnRnE89fQbz2zYGes2Dn91FIBTLt5wz0J4IhbuXwY2a26zA2th2ejLHbYQQuSSxC3Klfd+2cKaPccBcLEpvn2gLbvf6sHO66O4xsUxBcrD6Q+w8oSfWajbCrq85NjJXx/B9l8vZ9hCCJFLErcoN36NPcyEFY7iIU/cEEibRn4QuxD+eD93/dism/kluzUPz4ziaMops/Kqh6BJZ8fOvn8QUvIfJxdCiEtJErcoF/YnpvHEd+tzlzsGVefBa5pA/Fb4fmju+vSG1/GV++0AHEk+xWPfrMNu12CzQZ/PwLuG6XjyCMwbaiqKCSHEZSSJW5R5mdl2Hp6xlkRrXLtWRU8+6N8CW2YKfHMHZCSbjpUb4Nn/C94b0Cp32z+2HeWTZdvNgk8Nk7xzbF8K/35yud6GEEIAkrhFOfDuz1tYuzcRMOPa4wa2xM/LzZwxH91qOrl6woCvwcuPjkE1GNbRMZfZmKVb+WfHMbPQtDNc/ahj50tfgf1rL9dbEUIISdyibFu66TATf3eMaz91YxCtG/rBXx+ase0cvcZC7ea5i49fb41/A3YNj8yKIj7ZGu++7nmoY52V2zNhziA4lXzJ34sQQkARE7dS6jGlVIxSKlopNVMp5amUaqSUWqmU2qaU+kYp5W719bCWt1vtDYvjDQhRkLjjqaeNa3cKrsGQDo1hxzL49VVHxzYPQMSA07Z1dbHx8e0tqeptimrGJ59ixDdRZNs1uLrDrZPB3dd0TtgJPzx5yd+PEEJAERK3UqouZsrl1lrrcMwsjrcBbwNjtNYBwHFMHQWsf49rrZsCY6x+QlwSGVl2hs+IIinNjGvXqeTJ+/0isJ3YB7PvA23dVFa/LdyQ/3PZNSt6MmZAC5Q1/+lf248x7jdrvNuvMfQc4+i8YRasn3Wp3o4QQuQq6qVyV6CCUsoV8AIOYqZjnm21fwncbH19k7WM1d5ZKSUVDMUl8fZPm1m3z4xru9oUHw9sRRX3bPjmTkhLMJ18akG/qeYMugDXBFZn+HVNc5c//HUrf283k7TQvB+0uMPR+Ycn4NiO4n4rQghxmotO3Frr/cB7wF5Mwk4C1gCJWussq1scUNf6ui6wz9o2y+pf9cz9KqWGKKVWK6VWx8fHX2x4ohz7JeYQk//clbs8smsQkfUrm8R60Lp0bnOD/l+Bb61z7u/RzgFcaY13aw2PzFrHkeR009jtHahqJfaMFHM2L1OiCiEuoaJcKq+COYtuhCkx7A10y6erztmkkDbHCq0naq1ba61bV69e/WLDE+XUvoRUnswzrt0lpAb3d2gMq7+AddMdHbu+CfWvPK995ox3V/MxZ+ZHU07x6Mx1Zrzbwwdu/QJcrLP2g+vg11eK7f0IIcSZinKpvAuwS2sdr7XOBOYC7YDK1qVzAH/ggPV1HFAPwGqvBCQU4fhCnMaMa6/lRLq54FO3cgXe6xeBivsPfhzl6BhxO1wx+IL2XaOiJx8OaJk73v3PzmOM/XWbWagdAV3yJOt/xsG2JUV5K0IIUaCiJO69wFVKKS9rrLozsAlYBtxq9bkHcssXL7CWsdp/01qfdcYtxMV688dY1sclATnj2i2pnH0cvr3bPLYFUKuZuansIm6vaB9QjYc7BeQuj/1tG39us8a7rxoKATc6On//ICQfvuj3IoQQBSnKGPdKzE1ma4GN1r4mAqOAx5VS2zFj2DnFjScDVa31jwNPFyFuIU7zU/Qhpvy1O3f56W7BtKrrA7PvheSDZmWFKmaSFbcKF32cRzsH0LaxuTVDaxjxTRRHTqSbDwI3jzc3vAGkHoXvh8iUqEKIYleku8q11i9prYO11uFa67u01qe01ju11m201k211v201qesvunWclOrfee59i/E+dh7LJWnZjvGta8Prcmg9o1gyYuw5y9rrYK+k6FKwyIdy8Wm+Oj2FlTz8QDgaEoGD8+MIivbDt7V4JYJ5N7OsXM5/D22SMcTQogzycxpolQ7lZXNQzPWkpx3XPvWCNTG2fDveEfHzi+Y6UqLQQ1fT8be5ni+e+WuBD7KGe9u3BHaP+bo/NtrELe6WI4rhBAgiVuUcm8u3szG/WZc281F8ckdrah0YgsseNjRKbgntH+8WI/brmk1Hu3sGO8et2w7v2+1Hl+87lnwv8J8bc8yj4ilJxXr8YUQ5ZckblFqLd54kKl/785dfqZbCC2qaVPxKyvNrKwaADd/elE3o53Lw50CuLqpY7z7sW/WcfhEOri4Qd/PwaOi6Zi4BxY9bjoJIUQRSeIWpdKeYycZNXtD7vKNYTW5t119mHM/HN9tVrr7mJvRPCtekhhcbIoPB7Skuq8Z7z52Ms94d5WG0OtDR+fo2bBuxiWJQwhRvkjiFqVOeqY1rn3KjGv7V6nAO7dGoFa8DdvzPD9983ioEXxJY6nu68HY21pis07oV+1KYMxSq1RoeF9oeZej8+In4ei2SxqPEKLsk8QtSp03FscSvf8EYI1rD2xFpb2/woo8dWuuHgGhN12WeNo2qcpjXQJzlz9ZtoPlW46YhW5vQzWrLTPVPJ6WdeqyxCWEKJskcYtSZdGGA3z1z57c5ee6hxDhdQzmDnF0atwROr1wWeMadl1TOgRUy11+/Nv1HExKA3dva0pUczmdQxthyUuXNTYhRNkiiVuUGruOnuTpORtzl7uF1+Ke1tVMxa9T1l3blepB3y/AxbWAvVwaLjbFmAEtqGGNdyeczOCRnPHuWs1OLx268lPY8tNljU8IUXZI4halQnpmNg9NX0uKNa5d38+Lt/s2Qy14BI5sMp1cPGDANPA+q+jcZVHNx4OxtzvGu//bfZz3l1jj3W3uh6Dujs7zh8GJg5c/SCFEqSeJW5QKry3axKaDZlzb3cXGJwNbUTFqIsTMdXTqOQbqtHRShMZVjavyxA1BucufLt/Bss1HzONoN30CvnVMQ+oxmHs/2LOdFKkQorSSxC1KvAXrDzB95d7c5ed7htAsc72Z0jRH6/ug5R1OiO5sQ69twjWBjpK0j3+7jgOJaeDlB7dMJHdK1N1/wJ9jnBOkEKLUksQtSrSd8Sk8M8fxvHaPZrW5K8QFvrsXtHW26n8FdH3LSRGezWZTjOkfQc2KZrz7eGomD8+MIjPbDo06wDVPOTovewP2rXJSpEKI0kgStyixzPPaUZzMMAm6QVUv3rwpEPXdPab6FoB3dej/Fbh6ODHSs1X18eDj21vhYg14r9lznPd+2WIarx0F9a4yX+tsmD0I0hKdFKkQorSRxC1KrFcWbiI2Z1zb1RrXXvYc7F9jOigX6PclVKzjxCgL1qaRH0/c4Hi+e8KKnfwae9jc8d53EnhWMg1Je2HhozIlqhDivEjiFiXS/HX7mbnKMa79Ys9Qwg/NgzVTHZ1ufB0aXn35g7sAD17ThI5BjvHuJ75bz/7ENKhcH3p/7Oi4aR6s/coJEQohShtJ3KLE2RGfwrNzHc9r92xemzv8482UoTma9YMrH3RCdBfGZlN80L8FtSt5ApCYmsnwGWvNeHfoTRB5r6Pzj6MgfouTIhVClBaSuEWJkpZhntfOGdduVM2bt26shfr2bsjOMJ1qhkOvjy5Jxa9Lwc/bnY9vb5k73h21N5F3ftpsGm98A6qHmK+z0kwJ0Mx0J0UqhCgNJHGLEuWVhTFsPpQMmHHtcbc1w2fhEDix33TwrGQmWXH3dmKUF651Qz+eutHxfPekP3axZNNhcPcyU6K6mjNyDkfDkss7XasQonSRxC1KjO+j4pj1377c5Zd7hRG26UPzvDMACm75HPwaOyfAIhrSoTGdgmvkLj/53XrijqdCzVAzXp9j1UTY/IMTIhRClAaSuEWJsP1IMs/Ojc5d7h1Rh9u918DfYx2dOj4DgTc4IbriYbMp3u8XQR1rvDspLZPhM6LIyLJD60EQ3NPRef5DkLTfSZEKIUoySdzC6cy4dhRpmWZcu3E1b97q4Iaa/5CjU2DX0ycuKaWqeLvz8cBWuFrj3ev2JfL2T5vNeH3vj6Giv+mYdtxUPJMpUYUQZ5DELZzupQXRbDlsxrU9XG18emtTvObeDZknTQe/xtBnAtjKxq9rZIMqjOoanLs8+c9d/BxzyEyJ2ncSKOt97vkT/njfSVEKIUqqsvGXUJRac9bE8e3quNzlV3qFEPTPU5Cww6xw84IB06FCZSdFeGkM7tCILiGO8e6nvlvPvoRUaNDOzKyWY/mbsOcfJ0QohCipJHELp9l2OJnn5znGtfu0rMuA9G9hy2JHp5vGmZu3yhilFO/1i6Bu5QoAnEjPYviMtWa8+5qnoIE1sYy2w5zBkJrgxGiFECWJJG7hFKkZWQybvjZ3XLtJdW/ebHYItewNR6e2wyG8r5MivPQqe7nz8cCWuePd6+OSePPHWLC5mCpintZVhhNxsPARmRJVCAFI4hZO8sK8GLYdSQHA083GpF7V8FzwAGAlp4YdoMsrzgvwMmlVvwpPd3OMd0/5azc/RR+ESv6mfneO2IWw+gsnRCiEKGkkcYvL7rvV+5iz1jGu/XqPJjT+9UFItypk+daBW6eYYhzlwKD2jbg+tGbu8lOzN7D3WCqE9IQrBjs6/vwsHN7khAiFECWJJG5xWW05lMwL8x3j2re0rMMt+9+Fw9bc5C7uZmY0n+oF7KHsUUrx3q0R+Fcx493J6Vk8NGMtp7Ky4YbRUCPMdMxKN1OiZqQ6MVohhLNJ4haXzclTWQybvob0TDsATWv48Jb/P6iN3zo6dX8X/Fs7KULnqeTlxriBrXBzMePdG/cn8cYPseBWwZoS1SR14mPhl+ecGKkQwtkkcYvLQmvNC/Oi2RFvns2u4ObClOsycf81z7zcre6GyP85J8ASoEW9yjzTLSR3+ct/9rB440GoEQzd3nJ0XP0FbJrvhAiFECVBkRK3UqqyUmq2UmqzUipWKdVWKeWnlFqilNpm/VvF6quUUmOVUtuVUhuUUq2K5y2I0uC71XHMjXJM4flu1xrUWzoU7FlmRZ1W0O1dJ0VXctx7dUNuDHOMd4+avYE9x05Cq3sg9GZHxwUPQ+K+fPYghCjrinrG/RHwk9Y6GIgAYoGngV+11gHAr9YyQDcgwHoNAT4t4rFFKbH50InTxrUHtKxJz9hRcPKIWeFV1Yxru3k6KcKSQynFO7dGUM/PGu8+ZR6bS8+ym1KmleqbjulJ5vnu7CwnRiuEcAalL/LZUKVURWA90Fjn2YlSagvQUWt9UClVG1iutQ5SSk2wvp55Zr+CjtG6dWu9evXqi4pPlAwpp7LoPe5PdlqXyANr+rC46QJc13xuOigb3DUPGl/rxChLng1xifT99G8ys81/rbuuasBrN4fD3pUwpRtoaw7za0fBdc86MVJREiml1mitT7tZZM2aNTVcXV0/B8KRYdKSzg5EZ2VlDY6MjDxyZmNRnrdpDMQDU5RSEcAa4FGgZk4ytpJ3zryOdYG81/birHWnJW6l1BDMGTn169cvQnjC2bTWPPf9xtykXcHNha9a7cB12eeOTl1ekaSdj+b+lXmuewgvLzSPf037dw9XNvajZ/OwG33aAAAgAElEQVQr4bpn4LfRpuPv70Kja6BheydGK0oDV1fXz2vVqhVSvXr14zabTWbzKcHsdruKj48PPXTo0OdA7zPbi/KpyxVoBXyqtW4JnMRxWTw/Kp91Z/3yaK0naq1ba61bV69efh4JKotm/beP+esO5C6P62ij1h/PODqE3gztHnZCZKXDPe0a0r1Zrdzlp+dsZNfRk9D+cTNBDVhTot4vU6KK8xFevXr1E5K0Sz6bzaarV6+ehLk6cnZ7EfYdB8RprVday7MxifywdYkc698jefrXy7O9P3AAUSZtOnCClxbE5C7/r4UvnTc8YZ5FBqgebGYGU/l9nhNgxrvf6tuc+n5egBl2eGj6WtKzMVOiVvAzHZMPmPrdMiWqKJxNknbpYf2s8s3RF524tdaHgH1KqSBrVWdgE7AAuMdadw+Q89zKAuBu6+7yq4Ckwsa3RemVcipPwQwgpIYXL5x6H5L2mg4eFU3FLw8fJ0ZZOlT0dGP8Ha1wdzH/VTcdPMFrizZBxTpwc577O7cshv8+L2AvQpQu1157bdOjR4+6FOc+t2zZ4h4QEBBWnPt0lqLOKfkwMF0p5Q7sBO7FfBj4Vik1CNgL9LP6Lga6A9uBVKuvKGO01jw7dyM7j5pxbS93F6Y3XYrL2uWOTn0mQLWmzgmwFAqvW4kXeobwwnxzBWP6yr20aeTHTS26wpUPwsrPTMefn4P6V0GtZk6MVpQFLV79JSIxNfOs/FDZyy1r3Ys3rL9Ux7Xb7WitWbFixfZLdYyyoEh3Fmqt11nj0c211jdrrY9rrY9prTtrrQOsfxOsvlpr/ZDWuonWupnWWm4XL4NmrNrLgvWOEZApVx7Gb+04R4drRkJwdydEVrrdeVUDejSrnbv87NyN7IxPgetfdSTq7FPWlKgnnRSlKCvyS9qFrb8QL7/8cs2AgICwgICAsFdffbXGli1b3Bs3bhx255131g8LCwvdsWOHe926dZsdPHjQFeCpp56q3ahRo7B27doF9OrVq9GLL75YE6BNmzZBQ4cOrdusWbOQhg0bhv/0008+YM6sIyMjg0JDQ0NCQ0NDlixZ4l3UmEua8lHFQVwW0fuTeGWhowjGI83tXLk+z/ScTbtAx8LuXxQFMePdzYg+kMSeY6mczMhm2PS1zHvoajxvnQITroHMVDi6FX56Gnp/7OyQRQnW8OkfIi/Vtrvf6rGmoLY//vjDa8aMGVXXrFkTq7UmMjIypHPnzsm7d+/2nDRp0u6vv/56b97+v//+u9fChQurbNy4cVNmZqZq0aJFaMuWLXMn68/KylIbN26M/eabbyq9+uqrdbp27bq1Tp06WX/88cdWLy8vvXHjRo/bb7+9cXR0dOzFvt+SSJ7lE8UiOT3ztHHtljVdGHHsFchINh0qN4BbJpla0+Ki+Hq68cnAVri7mv+2mw8lmw9K1QLMHO851n4F0XOdFKUQBVu+fLlP9+7dEytWrGivVKmSvUePHseXLVvmW7t27YzOnTufdalo+fLlPt26dUv08fHRVapUsV9//fWJedv79et3HKBdu3Yn4+Li3AEyMjLUwIEDGwYGBob269evyY4dO8rczE6SuEWRaa15eu5Gdh8zH4S93W1Mq/oltmPbTAfXCnDbdPDyc2KUZUN43Uq82DM0d3nmqr3MX7cfWtwB4X0dHReOgON7nBChEAUraMIvLy8v+4X0z+Hp6akBXF1dyc7OVgCvv/56zRo1amTGxsZuss7Uy1yek0vlosi+XrmXHzY4HhD4NnwVPpsWOzr0+khumCpGd1xZn393HmOR9T1/Zu5GwupUomnPMRC3GhL3wKkkmDMI7v0RXNycHLEoaQq7nA2FXw4/17aF6dSpU8p9993X8LXXXjuktWbx4sVVpk6duvPLL7/Md9KOjh07pgwdOrRBamrqwczMTLV06dLKd999d3xhx0hKSnLx9/fPcHFxYdy4cVWzs7MvNtwSq8x9EhGXV/T+JF7LM679YuhhwmI/dHS48kGIGOCEyMoupRRv3tKMRtXMPTepGdk8NH0taTYfUwLUZn0ej/sPlr/pxEhFaVXZyy3fSfALWn++2rdvnzpw4MBjrVq1ComMjAy566674qtVq1ZgZr322mtTu3btmhQaGhrWvXv3Js2bNz9ZqVKlQjPxiBEjjsycObNqRERE8NatWz0rVKiQ79l8aXbRc5VfDjJXecl2Ij2TXh//yR7rEvm1NdKYmvkUKs2axat+O7hngZzxXSKbDpzg5vF/5d5XMKB1Pd6+tTn8OQaWvmz1UnD3fJlWtpzJb67y9evX746IiDjqrJguVlJSkq1SpUr25ORkW9u2bYM+++yzPe3bt08995al3/r166tFREQ0PHO9nHGLi6K15uk5G3KTdlUPOxM9PnQkbZ9a0G+qJO1LKLRORV7u5ZhP4pvV+5i7Ng7aPQqNO1prNcwdAidL3d9rIQC48847GwQHB4c2b948pFevXsfLS9IujIxxi4sy7d89LN54yFrSfN9gDh57N5pFmxv0/wp8axa4vSget7epx8pdx3LnhH/u+2ia+1eiaZ8J8OnVkHoUUg7BvGEw8BuZYlaUOgsXLtzl7BhKGjnjFhdsY1wSoxc5Hov8qGkU9fd+7+jQ9U2of6UTIit/lFK83qcZja3x7rRM83x3mkd16POZo+O2nx0zrAkhSjVJ3OKCJKVlMmzGGjKyzbhq3+r76X3gI0eHiIFwxWAnRVc++Xi48skdrfCwnu/eejiFF+dHQ8D10Ha4o+OSF+HgJZutUghxmUjiFudNa82o2RvYl5AGQAOPFN62v4+yZ5oOtZpDzw/kcqwThNSuyCu9HePd362JY/aaOOj8ItSOMCuzM8yUqKdSnBSlEKI4SOIW523q37v5KcaMa7uSxdxqE3E9aY1zV6gCA74GtwpOjLB8G3BFPfq0rJu7/MK8aLYdy4Bbp4CbNV3zse3w4ygnRSiEKA6SuMV5Wb8vkTcWO8a1p9VbSNVj1qN6ymaeH67SwEnRCTDj3aNvDqdJ9dPHu1N9G0CP9x0d130NG2c7KUohHB5//PE6OUVD8nPgwAHX5s2bB4eEhITmFBG53M4VozNI4hbnlJSayUMz1pKZbZ75H14tirbx3zk6dHoBmnRyUnQiL28PV8bfEYmnm/mvve1ICi/Mi4EWt0PzPBPhLBwBCXKzrijEf5P9eC+wGS9XjuS9wGb8N/myz1m8aNEi36ZNm6bHxsZu6tq163mN8WRlFWmOmFJBErcolNaap2avJ+64GdeO9NjP4+l5Kk8F94T2jzkpOpGfoFq+vHpTeO7ynLVxfLt6nznrrtLIrMxINlOiZmc6KUpRov032Y+fn2lAymF30JBy2J2fn2lQHMl71KhRtRo2bBjerl27wG3btnkAxMTEeHTo0CEgLCwsJDIyMigqKsrz77//rvDSSy/5L1u2rFJwcHBoSkqKmjBhgl9gYGBoQEBA2NChQ3PHhby8vFqOGDGiTvPmzYN//fVXn7p16zYbPnx43RYtWgSHh4eH/Pnnn17t27cPqFevXvg777yTO73qCy+8UDM8PDwkMDAw9LHHHqtTWIwliTzHLQr1xV+7+WXTYQAqksI0n7HYTqabxqoBcPOncjNaCdQv0p9/dx5j7tr9ALw4P5oI//YE3foFTL4e7Fmwfw38Nhquf8XJ0YrL7uVKF17WM+uUjR8eb8QPjzcqfN9JhZb1/P777/2s4h/klOkcPHhwg4kTJ+5p1qzZqd9++8176NCh9f/999+tzzzzzIHVq1d7f/XVV3t3797t9vLLL9dds2ZNbPXq1bM6dOgQOG3atMp33XVXYlpami08PDztww8/PJBzrHr16mWsW7du86BBg+rdd999DVeuXLnZ6hc2cuTI+Llz51bcvn2754YNG2K11nTp0qXpjz/+6OPj42PPL8YL/n5dQpK4RYGi9h7nTWtcW2FnTo0peJ3YZxrdfUzFL8+KToxQFCRnvHtDXBLbj6SQnmln2PQ1LBjeHu/OL8GSF0zHvz4006HKUIe4DJYtW+bTvXv3RF9fXzvADTfckJienm6Liory6devX5OcfhkZGWedDfz555/eV111VXKdOnWyAAYMGJCwYsUKn7vuuivRxcWF//3vf8fz9u/fv38iQLNmzVJPnjxpq1Klir1KlSp2Dw8P+9GjR11++umnir///nvF0NDQUIDU1FTb5s2bPZOTk21nxnjpviMXRy6Vi3wlpmYwfEYUWXYzrv1mlR8IOPGPo8PNn0L1ICdFJ86Hl7sr4+9oRQU3UwN9R/xJnp8XjW77EDTp7Og49wFIKbTgkhDFRp1xhc5ut+Pr65u1efPmTTmvnTt3xpy5XWF1Ndzd3e2urqefh+aU/LTZbLi7u+dubLPZyMzMVFprRowYcTDnmHv37o1+7LHHjuYXY0kjiVucRWvNk99tYH+iGdfu6bmO29JmOjq0fwxCezspOnEhAmv68trNjvHu76P28+2a/WZWNe8aZuXJIzDvQbCXuSJKoiAvJ60p9NXjg124epz+C+HqYafHB7vOuW0hOnXqlPLDDz9UTklJUcePH7ctWbKkspeXl93f3z/jiy++qAImkf/zzz9nPVd6zTXXnFy5cqXvwYMHXbOysvjuu+/8OnbseNGTEnTr1u3EtGnTqiUlJdkAdu3a5bZ//37X/GK82GNcKpK4xVkm/7mLpbFmXLuhOsgY1/GOxsYdzV3kotS4NdKffpH+ucsvzo8hNtnz9ClRty+Ff8fns7Uol64YlMCNb+7Bp2YGKPCpmcGNb+7hikEJRdlt+/btU/v06ZMQHh4e1rNnzyZt2rRJAZg5c+bOKVOmVAsKCgoNCAgImzNnzlnJskGDBpkvvvji/muvvTYwJCQkrHnz5ql33nnnRV/GvuWWW07069cv4YorrggODAwM7dOnT5PExESXgmIsSaSspzjN2r3H6f/ZP2TZNV6ks7zyaGqk7zSNlerDkOXgXdWZIYqLkJaRzU2f/MnWw+ZvUONq3ix4uD0+K16Bv8eaTjY3GLwE6rR0YqSiOJSlsp7lmZT1FOeUmJrB8OlrrXFtzYRKUx1J28UDBkyTpF1KVXB3OW28e+fRkzz3/UZ0p+ehTivTyZ5pTYma7MRIhRDnIolbAGC3a574dj0HksyjXg95/kyHU787OvQcA3VaOCk6URya1vDl9T6O8e756w4wa+1huHUyuPualQk74YcnnRShEOJ8SOIWAHz+505+3XwEgLa2GJ5Q0x2NrQdByzucFJkoTre08mdA63q5yy8tiGFTejXzwSzHhlmwfpYTohNCnA9J3II1exJ4+6ctANTmGJ97jcems02j/xXQ9S0nRieK28u9wwiqac6wM7LsPDRjLSlBfUxJ1hw/PAHHdjgpQiFEYSRxl3MJJ83z2tl2jQcZfOXzMd5Z1jwG3jWg/1fg6u7cIEWxquDuwid3tMLL3Yx37zp6kmfmbkR3fwf8rDkwMlLMeHdWhhMjFULkRxJ3OWbGtddx0BrXft1zGgFZW02jzRX6fwkV6xSyB1FaNa3hwxt9muUuL1x/gOlRCabKm83NrDy4Dn6V6VCFKGkkcZdjE37fybItZsasAS7LuJVfHY03vA4N2jkpMnE53NyyLre3cYx3v7poE9G6EVz/qqPTP+Ng21InRCfKg5ySmSNGjKgzb9483wvdftGiRb7XXXdd00sRW3GbNm1a5TVr1ngWx75krvJy6r/dCbz3ixnXjlDbed19KuQ80t+sP1z5gNNiE5fPS73CiNqbyOZDyWRk2Rk+Yy0Lhw/Gd+cy2PaL6TTvQXjwL/AtUSWJRREpuPBCI4XQUOisaYXJWxykrJo3b17lrKyspMjIyPSi7kvOuMuhYymneNga165KEpMrjMVVW+Uda4ZDr4+k4lc54elmnu/2tsa7dx9L5envo9E3jQcfK1GfjIfvH5ApUUWxyK9kZt++fRtOmTKlCsCwYcPqNmnSJCwwMDB0yJAh/jntAwcOrB8ZGRnUsGHD8JkzZ1Y6c7/Lli3zatmyZXBISEhoy5Ytg9evX+8Bpj73kCFD/AMDA0MDAwNDX3/99RpgKpVdccUVQWFhYSHt27cP2LNnjxtAmzZtggYNGlSvdevWQY0bNw5bsWKF1w033NCkQYMG4Y888kju2OH48eP9mjVrFhIcHBw6cODABjl1wL28vFo+/PDDdYOCgkIjIiKC9+3b57pkyRLvpUuXVn7++ef9g4ODQ2NiYjxGjx5dI+d99uzZs/GFfA8lcZczdrvm8W/Xc+hEOi5k86nnOKrZrcmUPCuZSVbcvZwbpLisGlf34Y1bHOPdP2w4yNcbT8ItEwHrA9zOZY4Z1oS4SHnLei5atGj7+vXrvfO2Hz582GXx4sVVtm3bFrN169ZNb7zxxsGctn379nmsWrVqy8KFC7eNGDGiQWpq6mlnFxEREemrVq3aHBsbu+mll17aP3LkSH+A999/v/qePXs8YmJiNm3dunXT4MGDj506dUo98sgj9efPn78jJiYm9p577jn65JNP5tb3dnd3t69evXrLvffeG9+vX7+mkyZN2rt58+aYb775ptqhQ4dc1q5d6zl79my/1atXb968efMmm82mP/vss6oAaWlptrZt26Zs2bJlU9u2bVM+/vjj6tdff/3JLl26JI4ePTpu8+bNm8LCwk6NHTu2VnR09KatW7dumjp16p4L+T4W+VK5UsoFWA3s11r3VEo1AmYBfsBa4C6tdYZSygP4CnN55hgwQGu9u6jHFxfm0xU7WLHVjGuPdJ1FG3KK8CjoOxn8LuiDnygjbmpRl5W7Epixci8Ary2KpeWwdoS3fwz+/MB0+u01aNgB/Iv1CqsoR/Ir65m33c/PL9vDw8N+2223NejRo0fSgAEDknLa+vbtm+Di4kKzZs1O1atX79S6detOGy9OSEhwGTBgQKPdu3d7KqV0ZmamAvjtt98qPvjgg/Fubuamy5o1a2b/999/ntu2bavQqVOnQDCFTapXr56Zs68+ffokAkRERKQ1bdo0rUGDBpkA9erVO7Vz50735cuX+0RHR3tFRESEAKSnp9tq1KiRBeDm5qZvu+22JIDIyMiTS5cuzbf2cVBQUFqfPn0a9e7dO/GOO+64oDnXi+OM+1EgNs/y28AYrXUAcBwYZK0fBBzXWjcFxlj9xGW0cucx3rfGtXvY/uUB1x8cjdc9CwHXOykyURK82DOUkNrmb0xGtp1h09dyou1T5ll+AHsWzLkP0k84MUpR2hVWMtPNzY1169bF9u3bN3HevHmVO3bsGFDQdmcujxo1qu61116bvG3btpiFCxduz8jIsIGpdqiUOq0oh9ZaNW3aNC2npOfWrVs3/fXXX9ty2vOWBPXw8DitJGhWVpbSWqt+/fody9l+9+7d0R988MEBAFdXV22zmdTq6upKVlZWvm942bJl2x566KH4NWvWeEdERIRmZmbm1y1fRUrcSil/oAfwubWsgE7AbKvLl8DN1tc3WctY7Z1VSS96WoYcTTnFI7OisGsIVPt4z2OiozGwG3SQaS7Lu5zxbh8PcyFub0Iqo76PRd8yCTysk4bju2HRY1CCixOJkutcJTOTkpJs1plz0meffbYvNjY2d9xu7ty5VbKzs4mJifHYt2+fR0RExGk3eZ04ccLF398/A2DChAnVctZ36dLlxGeffVY9JzEePnzYpXnz5ukJCQmuS5cu9QY4deqUWr169Xnf8d21a9cTixYtqrJ//37XnH1u3bq10AkvfHx8sk+cOGEDyM7OZseOHe69evVKHj9+fFxycrJLUlKSy/kev6hn3B8CI4Gcu1aqAola6yxrOQ7IGTeoC+wDsNqTrP6nUUoNUUqtVkqtjo+PL2J4Asy49mPfrOPwiVP4ksokjw+poK3feb8mcMsEsMntDgIaVfPmrb6O8e4fow/x1WYFvT50dIqeDetmOCE6Udqdq2RmYmKiS9euXQMCAwNDO3ToEDR69Oh9OW1NmzY91aZNm6AePXoEfPjhh3u8vLxO+/Q4atSoQy+//LJ/q1atgrOzs3PXP/bYY/H+/v4ZwcHBYUFBQaGTJ0/28/T01LNmzdrx9NNP+wcFBYWGhYWFrlixwud830dkZGT6888/v79z586BgYGBoZ06dQrct2+fW2Hb3HHHHQljx46tFRISEhodHe0xcODARoGBgaHh4eGhDzzwwOFq1aplF7Z9Xhdd1lMp1RPorrUeppTqCDwJ3Av8Y10ORylVD1istW6mlIoBbtRax1ltO4A2WutjBR1DynoWj3G/beO9X7aisDPRbQzXu1hPbbh5w+ClUDPUuQGKEueFedFM+9fcL+PmopgztB3N1zwPUdNMBzcveOB3qBZQyF6Es5S1sp59+/Zt2LNnz6R77733uLNjuZwuRVnPq4HeSqndmJvROmHOwCsrpXJuevMHcp7PiwPqAVjtlYAiFWUX5/bPjmN8sMTMhvaQy3xH0ga4aZwkbZGv53qEEFbHXB7PzNY8NGMtSR1HQ7VA0yEzFWbfC1mnnBilEOXTRSdurfUzWmt/rXVD4DbgN631HcAy4Far2z3AfOvrBdYyVvtv+mJP98V5iU92jGt3tK3jcbfZjsa2wyH8FucFJ0q0nPFuX2u8e19CGiMXbEf3nQwu1lDeoY2w5CUnRinKizlz5uwub2fbhbkUA5ujgMeVUtsxY9iTrfWTgarW+seBpy/BsYUl2xrXjk8+RX11mLHun2DLmRqtYQfoInNQi8I1qOrN27c2z13+OeYwU3f6wg2jHZ1WfgpbfnJCdEKUX8WSuLXWy7XWPa2vd2qt22itm2qt+2mtT1nr063lplb7zuI4tsjfJ8u28+f2o3hyigluY6jISdNQsS7cOgVcZLZbcW7dm9XmnrYNcpffWBzLutr9zZMIOeYPgxMH89laCHEpyK3EZdDfO47y4dKtgOYtt0mE2MykGri4Q/9p4FPdqfGJ0uXZHiE0q2tmmMzM1gyfGcWJGz4E39qmQ+oxmHs/2M/7plghRBFI4i5jjiSn88jMddg1/M/lZ252+dvR2P09mfVKXDAPVxc+GdgKX09zlSbueBpPLI5D550Sdfcf8OcY5wUpRDkiibsMybZrRsxax9GUU1yhNvO823RHY6t7IPKegjcWohD1q3rxbp7x7iWbDjM5zh+uyTNxz7I3YN8qJ0QnypuxY8dWvfvuu+sDvPPOO9XHjRt31pwghfHy8mqZ3/rRo0fXaNy4cVjv3r0bFUecl4oMdJYhY3/dxt87jlGD44x3/whXrEuXdVpB93edG5wo9bqG1+Z/7Roy9e/dALz142YihzxAy12/w76VoLNh9iB48A+oULnwnQnnU6p4L79pfdFlPYti5MiRxTZT1+TJk6v/+OOP24KDgzPOp39mZiY5c6BfTnLGXUb8tf0oY3/bhhtZfOr+IdWVNTe/VzVT8cvVw7kBijLh2e4hRPib8e4su2b4rI0kdf/UVJYDSNoLCx+VKVFFvrp06dIkLCwspGnTpmHvvfde7rSkXl5eLe+//37/0NDQkLZt2wYeOHDAFUyJzfvuu69ey5YtgwMCAsKWLVt2VunCxx9/vM6LL75YEyAmJsajQ4cOAWFhYSGRkZFBUVFRngCbN292b9GiRXB4eHjIo48+WufMfQAMHDiwflxcnEfv3r2bvvLKKzUOHz7s0qVLlyaBgYGhERERwStXrqyQc7zbb7+9wdVXXx1wyy23NLrQsqFFKeeZQxJ3GXDkRDqPzopCa3jBdRqRNmuufGWDflOgkr9zAxRlhrurjXEDW1HRGu/en5jGE78koHvlKfm5aR6s/cpJEYqSbPr06btjYmJi161bt2nChAk1Dx065AKmFGarVq1SN23aFHv11VcnP/3007nJNTU11RYVFbV57Nixe4YMGVLoJezBgwc3GD9+/N6YmJjYd999N27o0KH1AYYNG1Z/8ODB8dHR0bG1atXKt5rHjBkz9taoUSNzxYoVW1966aUjI0eOrBMREZG6devWTa+99tr+e+65J/fYGzZs8Pr555+3L1y4cNeFlg0tSjnPHJK4S7lsu+aRWVEcTcngVpcV3O26xNF4/avQ6BrnBSfKpHp+XrzbLyJ3eWnsET4/1hwi/+fo9OMoiN9y+YMTJdrbb79dMygoKDQyMjLk0KFDbjExMZ5gqm4NHjw4AeC+++47tmrVqtx5wwcOHJgA0K1bt5SUlBTb0aNH8y3GkZSUZIuKivLp169fk+Dg4NBhw4Y1OHLkiBvA2rVrfe6///4EgAceeKDAabbzWrVqle+gQYOOAfTu3Ts5MTHR9dixYy4AXbt2TfTx8dGQf9nQDRs2eOSUDQ0ODg599913ax84cMANHOU8x48f7+fm5nZRl6ZkjLuU+2jpVv7dmUCY2sXrrl84GsL6mNnRhLgEbgyrxX1XN+KLv3YB8PZPm4kcNIpWe/+F+M2QlQaz74PBv4LbeRddEmXYokWLfFesWOG7evXqzb6+vvY2bdoEpaWl5XvymLdw5LnKeebIzs7G19c3a/PmzZvya7fZbBeUJPOb2DOnPKi3t7c9b7+CyoauW7du85n7WLZs2bYff/zRd968eZXfeeedOtu2bYu+0HFyOeMuxX7fGs/Hy7ZTmWQmuI/BQ1lXgKqHQO9xIFVTxSX0dLdgIuqZm9Cy7Jrh38aS1GMCuFqJ+nA0LHnBiRGKkiQxMdGlUqVK2b6+vvaoqCjP9evXe+e02e12pkyZUgVg6tSpVdu0aZOc0zZz5swqAD///LOPr69vdtWqVfOdMMDPz8/u7++f8cUXX1TJ2ec///xTAaBVq1YpkyZN8gOYNGnSed2BftVVVyVPmTKlKpgPHVWqVMny8/Ozn9nvQsqGFrWcZw5J3KXU4RPpPPbNOpS2M9ZtHP7KKvrjURFumw4e512hToiL4u5q45OBLXPHuw8kpfPY8kzsN7zu6LRqImz+wUkRipKkb54wvUgAABpqSURBVN++SVlZWSowMDD02WefrRMREXEyp61ChQr2mJiYCmFhYSG///6775tvvpk7FV+VKlWyW7ZsGTx8+PAGEyZM2F3YMWbOnLlzypQp1YKCgkIDAgLC5syZUxlg/PjxeydOnFgjPDw85HwT5dtvv31g7dq1XoGBgaHPPfdc3alTp+7Kr9+FlA3NyspSRSnnmeOiy3peDlLWM39Z2XYGfr6SVbsSeMp1Fg+5LnA03j4LgroVvLEQxWzJpsPc/5Xj/+kzXYN44NBLsHmRWVGhCjz4F1Sq66QIy5/SVtbTy8urZWpqatSZ69u0aRP03nvv7bvmmmtSnRGXs12Ksp7CST5cuo1VuxK40bbq9KR9zUhJ2uKyuz60Jvd3cNzs+84vW4lq+ZqZFx8g7TjMHSJTogpRTCRxlzIrtsbzyfLtNFH7ed/tM0dD0+uhoxRcE84xsmswLeub8e5su2bY3F2c6PGpeSQRYM+f8Mf7ToxQlGT5nW0DrFq1akt5PdsujCTuUuRQkhnX9tapTHAbg49KNw1VGsItE8F2wfc4CFEs3FzM892VKpi7Yw8mpfPI3xWwXzPS0Wn5m7DnHydFKETZIYm7lMjKtvPIzCgSTp7iXbcJNLUdMA2uFWDA1+Dl59wARblXt3IFPujveL57+ZZ4JupboH47s0LbYc5gc+lcOIPdbrfLoyalhPWzOusudpDEXWp8sGQrq3Yn8KDLQrq5/Odo6D0WajVzXmBC5NE5pCYPXOOYxfHdpTtYf+W74GnNXX4iDhY8LFOiOkd0fHx8JUneJZ/dblfx8fGVgOj82mUCllJg2ZYjjF++g/a2jTzl+o2j4cqh0Ly/8wITIh9P3hjE6j3HWbPnONl2zZD5h1jS7UMqzv+f6RC7ENZMgdb3OTXO8iYrK2vwoUOHPj906FA4ctJW0tmB6KysrMH5NcrjYCXcgcQ0eoz9A++0Ayxwfw4/lWIa6reDexaAy+WvTCPEueT83h5PNZNSXBNYnS9rfINa/bnp4OoJ9y+DmqFOjLLsyu9xMFF2yKeuEiwz287DM6NITT3Jp25jHEnbtzb0mypJW5RYdSpX4IP+LXKXf98azwTPe6FGmFmRlW6mRM1Mc1KEQpRekrhLsPd+2cKaPQm87vYFzWy7zUqbG/T/CnxrOjU2Ic7luuAaPHhtk9zld37dw4arPjA3VALEx8LPz/6/vXuPs6ne/zj++uy958JgDGPKNZeGU6dyTXLphFxydNSJhKLkuBQRXeh6eqTfkS4cinJSxyU5jkspurjUccpByF0xIQbDuBvM/fv7Y62ZvTGEmdlrr5nP8/HYj9nru75jv7+W9qe9vmt/l0PplHIvLdwhaulPB3jvPzt4wLuYzt5l/h13joKqjZ0LptRleLJtbRpdEwNAtoE+X6RwstVIf4fVH8CW+Rf4baVUXrRwh6C9x84wdNZ6Gsg2XvQF3Ne4Xg9o9IhzwZS6TD6vh/Hd61MuKhyAgyfTeHTLDZjr7/Z3mj8Qju1xKKFS7qOFO8RkZGUzaMZawk4nMyH874SLvUxkxbrwxzf1jl/KdSpGn/397v8mHGZS9GCIrmY1pB6HuX+BrEyHEirlLlq4Q8zrX/3Mht2HeDt8HFeLvVBFiXLWIithJZwNp9QVur1OHI/e7p/vfu3b/Wxs8gaIvdrf7v/BstEOpVPKXbRwh5DFWw4wadkOnvXN4BaPff918UDnyVC2mrPhlMqnoW1q07iGtcJftoHeS72cahqwJOqy12HXdw6lU8o9tHCHiMSjpxn27/V08nxHb9+X/h2tXoBarZwLplQB8Xk9jO9Wn/L2fHfyyTT67WqBuaa51cFkw5y/wOkjDqZUKvRp4Q4B6ZnZDJzxI5VSExgV9r5/x+86QvMnnAumVAG7qkwkY7rWy71U47tfjjE57llrOgjg5D74dKAuiarURWjhDgGjv/yJnXsSeS/sLUpIutUYWxvunqgXo6ki57baFRjY8trc7Ve/O8bWW0b5O/y8AH54P4/fVEqBFm7Hfb05icnf/cLYsHeo5km2GsNLQ9ePILKMs+GUKiSDW8dziz3fbQw8+F15TtcPWJb5q+cgKc/7KyhV7GnhdtCeI6d58t/rGeKbS0vvev+OeyZChdrOBVOqkOXMd8eWsua7D6Wk0T+pE+aqG6wOWWnWkqjppxxMqVRouuLCLSJVReQbEdkqIptFZLDdXk5EFonIdvtnjN0uIjJORBJEZIOINCioQbhRemY2Az/+kZvTVzLYN9e/o/lQuO4u54IpFSRxZSIZ27V+7mzQsp0nmVLpJQgraTUc+hm+HOFcQKVCVH4+cWcCw4wx1wFNgMdE5HpgOLDEGBMPLLG3Ae4E4u1HX2BiPl7b9UZ98RPHE7cyJmyCv7FmS2j1vHOhlAqy5vGxDGoVn7v98ooMtjd80d9h7RTYPM+BZEqFrisu3MaY/caYtfbzk8BWoDLQCZhid5sC5Kxt2AmYaiwrgLIiUvGKk7vYl5uSmPn9Vt4LG0MZse+OFF0NOn8AHq+z4ZQKssGt47m1ZnnAmu/u9kMtUusELok6GI7+6lA6pUJPgcxxi0h1oD6wErjKGLMfrOIOxNndKgOBCxIn2m3n/ll9RWS1iKxOTk4uiHghZffh0zw1ex2jwyZRx5MIgPFFQtdpULKcw+mUCj6vR/h7t3rElooA4NCpDPofexBT9hqrQ9pxmNMHsjIcTKlU6Mh34RaRUsAcYIgx5sTFuubRdt6XNY0xk4wxjYwxjSpUqJDfeCElLTOLgR+vpWvGfDp6V+S2S8cxUKneRX5TqaItrnQk4+73f7/721/TmFH1JfD4rIbEVfDtqAv/AUoVI/kq3CIShlW0PzLG5FxhdSDnFLj986DdnghUDfj1KsC+/Ly+2/xt4U9E7VvOCN8Mf+PNfaBed+dCKRUiml4by+DW/vnu51dHsvPGIf4O/30TdvzHgWRKhZb8XFUuwGRgqzHmrYBd84Fe9vNewKcB7T3tq8ubAMdzTqkXB19s3M9Xy9cwPmw8XrFPNFRpDO3+5mwwpULIoFbxNLvWP9/dZWNj0qrdZu81MLcvnDrsXEClQkB+PnE3Ax4EWonIOvvRARgFtBGR7UAbextgIbADSAD+ATyaj9d2lV8Pn+L52auZGD6GWLFmE0xUHNw3FXzhDqdTKnR4PcLYrvWpUNqe7z6dyaDUfpiSsVaHlCT49FFdElUVa/m5qvw7Y4wYY24yxtSzHwuNMYeNMa2NMfH2zyN2f2OMecwYU8sYc6MxZnXBDSN0pWVm8diMtQzLmkw9zw4AjMeH3DcFyhTLi+qVuqgKpSMYd399PPZ899e7hTnVnvV32PYlrHzPmXBKhQBdOa2QvbpgKzckfUJ33ze5bdLu/+Capg6mUiq03VqrPE/c4V898Ml1V7OnzsP+DotegP3r8/hNpYo+LdyFaMGG/axfsYSXff/0N954HzTu61gmpdzi0ZbX0iI+Nne78/a2ZMTdZG1kpVtLoqalOJROKedo4S4kuw6d4rU5/2Vi+FgiJBPAWof5rr/rHb+UugRejzCmaz3i7PnuA6cNw7IHYcKirA6HE+CLZxxMqJQztHAXgtSMLAZO/4HXssdQSY4AYCKika7TIbykw+mUco/YUhGM6+af756fGMWCak/6O6ybDhtnOxNOKYdo4S4EIxds4a5Dk7jVuwUAgyCdJ0O5Gg4nU8p9mtQsz7C2dXK3B26uTdI1f/J3+GwIHNnpQDKlnKGFu4B9tn4fR1fNop9vQW6btHwO4ts4mEopdxvwh1rcVjtnJUWh8+57yYyubm2mn4Q5j+iSqKrY0MJdgHYeOsX7cz5ndJj/qyqmzp3QYpiDqZRyP49HGHNfXa4qY813J54J41nvEEzOkqh718DSkQ4mVCp4tHAXkNSMLJ6atowxvEmUpAGQFVMTuec98Ohfs1L5Vb5UBOO7NcBrT3jP2hfH0soD/B2+Hwu/LHUonVLBoxWlgLzy2Sb6HRlNTU8SAFm+kni7zYDIaIeTKVV0NK5RjmFt/d/v7rP9Fg5f3cLfYV5/SCl6dxVUKpAW7gLw6bq9lFszjjbetblt3rvfgbjrHEylVNHU/7Za3F7Hmu82eOhyoCdZuUuiHoBP+kN2toMJlSpcWrjz6ZfkFBbOncoTvjm5bebWQXDDnx1MpVTR5fEIb91Xj4rRkQDsOBPFK2EBdxFLWAwrJjiUTqnCp4U7H1Izshg5dQGjZRwe+45fmdWaI3f81dFcShV15aLCGd+tfu589z8P1OR/V/fwd1j8V9j3ozPhlCpkWrjz4dVP1vD08ZFEy2kA0qMq4es6Bbw+h5MpVfQ1ql6Op9r5v9/dc1c7jpe70drIzrCXRD3pUDqlCo8W7is0b+0eGm54ies8uwHIkjDCu0+HqNjf+E2lVEHp26ImrX4XB0AGProf6Ut2eClr55EdsODJi/y2Uu6khfsKJBxMYesnb3C3d3lum6fjm1C5oYOplCp+PB7hzS51qWTPd29OLc/YiICviG2YCetnOpROqcKhhfsynUnPYuKUqTwl03LbMur1RBr2cjCVUsVXTFQ447s3wGfPd49Lrs/68h38HRYMg8O/OJROqYKnhfsyvTnnG4anjCJMsgA4E1ePsI5vOJxKqeKt4TUxPNP+d7nb3fZ25lSp6tZGeoo1352Z7kw4pQqYFu7LMPeHHXTY+gwV5DgAqeHlKNHjI/BFOJxMKdWnRQ3uuM6a7z5NJL1T+mM8YdbO/etgycsOplOq4GjhvkTbD5wk9bNnaOBJACALDxH3/xOiqzgbTCkFgIjwRpe6VC5bAoCVqdX4oETAFNb/3obtix1Kp1TB0cJ9CU6nZ/LJh6/T3fN1bltm65eRmn9wMJVS6lxlS4Yzvnv93PnuVw63ZFuZW/0dPukPJw84lE6pgqGF+xK8N3Muj5/xr8R0otZdRDQf5GAipdSFNKgWw/A7c+a7hW4He5Iaad8S9FQyzOunS6IqV9PC/Rs+Xb6BLr+MIEKse/0eKx1Pma7vgYjDyZRSF/JI8xq0uf4qAA4TzcDU/hjs/2Z3fAPLxzmYTqn80cJ9Edv2HyP2q0epIocAOOOJIrrXTAiPcjiZUupiRIQ3OtelSow137049TpmR3b2d1j6CiSucSidUvmjhfsCTqVlsubDoTSTjbltnnsnIbHXOphKKXWpokuG8Xb3BoR5rU/aI47dxZ6S11s7szNhTm9IPeFgQqWujBbuPBhjmDV1At3S/Xf8OtxwCBG/7+hgKqXU5apXtSwj7rRur5uJj25H+5Lhs5dEPboLPn8CjHEuoFJXQAt3Hr78dhldEl/N3d4f14Lyf3zJwURKqSv1cLPqtPu9Nd+daOIYkdHHv3PTbFg3w6FkSl0ZLdzn2LZ7L7W/7U8pSQXgUFglKj48DTz6V6WUG4kIozvXpWo5a757dlpjvg5v4++w8Ck4tN2hdEpdPq1GAU6lZnBgam9qyT4AUgknqtdMKBHjcDKlVH5ElwjjnYD57sEnupEcUc3amXHKXhI1zcGESl06Ldw2YwxLJ4+gReaK3LZjrd+kRJW6DqZSShWUm6qU5bkO1nz3GSLpdaI/WTlLoiZtgMV/dS6cUpfB53SAUPGfhf+iw8H3yfmq5/YaDxLfoqezoZRSBapX0+qs2nWEhRuT2GKq81rWAzwrH1o7V0yAmrdD7XaQkgK//gr79sH+/XDgACQnW4+jR6FSJXj7bZ1CU47Qwg0kbNtM3VXD8Ip1demOknWJf2CMw6mUUgVNRBh1701s2nuC3UdOMyntDlqV3kyTjFVWh390gXdTIOUiV5pXrapFWzkq6P/yRKS9iPwsIgkiMjzYr5+j0chFPP7sCPa+WJNaHzUlRlIASDIxVOwzE7xhTkVTShWiMpFhTOjRgHCvBxAGnOxNkrGvY4kSsoeWxrxYBp4vDe39d/4zwMru3dm+axdU0ZsLKecE9RO3iHiBd4A2QCLwg4jMN8ZsCWYOgKanlzIq7H1Kytn36N1vYtj52TvBjqOUCrInIo9xIjUTgG3ZVbjaexQAT85qxl6BxuEApC/J5sFp05jVtSsAZYCmwACgI3qxkAquYJ8qbwwkGGN2AIjITKATEPTC/bRv1nlFG6C+ZwfsfDvYcZRSQXYrwG+dWBPB3BzBNZO3klSpUm7zCeBL++EDbgC6YxXyUoWSVim/YP+PYmVgT8B2ot2WS0T6ishqEVmdnJxcOCkSE6lkrz+ulFIX5eGson2uLOAkcBBICVYmVawF+xN3XrfUOusqEGPMJGASQKNGjQp2LUJjYPp0GDSIfY/H5t48JNBJU4JNlboU6MsqpULTuj1Hc5/3836OR85/y8kS73ltEUAD4CH7EV5I+ZTKS7ALdyJQNWC7CrAvKK988CD07w/z5gEwOnPAeXPcp004z2U8zLh+fwtKJKWUs7oNX5D7vCRn6OldfNYdew3wbsOHACgP3A4MAZoHL6JS5wl24f4BiBeRGsBe4H6sqaHCNW8e9OtnfQcToEIFlnubMTzDmuuuJIfZZ8ozOvM+lpdsVehxlFKhQfCf8nspszcAPbxL8ZINHi9zGj7Erx3fYi9w4ZPlSgWXmCDfGUdEOgBjAS/wgTHm1Qv1bdSokVm9evWVv9ixY/D44zBtmr/tnnvg3XchLu7K/1yllAphIrLGGNPI6RyqcAR9ARZjzEJgYaG/0KJF0Ls3JCZa29HR1qIJPXpw1rkwpZRSykXc//XD7Oyzt0+dgsceg7Zt/UW7TRvYtAkeeECLtlJKKVdzd+HeuBEmTvRvL18O9erBhAnWdsmS1vOvvtKVjpRSShUJ7i3chw9Dp06wfj2kpcHw4dCiBSQkWPubNbP2DRign7KVUkoVGe68yUhmJnTtCjt3wvffw803W5++AcLDYeRIGDoUvOd//1IppZRyM3cW7qefhiVLrOdbAlZLrV8fpk6FG25wJpdSSilVyNx3qnzqVBiTxy03H3kEVqzQoq2UUqpIc1fhXrUK+vbNe9/06fDCC3D8eHAzKaWUUkHknsKdlAR//rN1IVpeGjSA2Fgt3EoppYo0d8xxp6XBvffC3r3+Np8PWra0inmnTlCxonP5lFJKqSAJ/cJtDAwaZH1Hu0QJaN/eWra0Y0eIiXE6nVJKKRVUoV+4Z82C1FSYOxfatbMWVVFKKaWKqdAv3F26WN/ZVkoppZQLLk7zhH5EpZRSKli0KiqllFIuEvT7cV8OEUkGfg3CS8UCh4LwOsFSlMZTlMYCRWs8RWksULTGU8cYU9rpEKpwhPQctzGmQjBeR0RWF6Wbzhel8RSlsUDRGk9RGgsUrfGIyGqnM6jCo6fKlVJKKRfRwq2UUkq5iBZuyySnAxSwojSeojQWKFrjKUpjgaI1nqI0FnWOkL44TSmllFJn00/cSimllIto4VZKKaVcpNgXbhFpLyI/i0iCiAx3Os9vEZGqIvKNiGwVkc0iMthuLycii0Rku/0zxm4XERlnj2+DiDRwdgTnExGviPwoIp/b2zVEZKU9ln+JSLjdHmFvJ9j7qzuZOy8iUlZEZovIT/YxutWtx0ZEnrD/jW0SkY9FJNJNx0ZEPhCRgyKyKaDtso+FiPSy+28XkV5OjMXOkdd4Xrf/rW0QkXkiUjZg3wh7PD+LSLuAdle956k8GGOK7QPwAr8ANYFwYD1wvdO5fiNzRaCB/bw0sA24HhgNDLfbhwOv2c87AF8AAjQBVjo9hjzGNBSYAXxub88C7refvwsMsJ8/CrxrP78f+JfT2fMYyxSgj/08HCjrxmMDVAZ2AiUCjslDbjo2wG1AA2BTQNtlHQugHLDD/hljP48JofG0BXz289cCxnO9/X4WAdSw3+e8bnzP08f5j+L+ibsxkGCM2WGMSQdmAp0cznRRxpj9xpi19vOTwFasN9lOWEUD++fd9vNOwFRjWQGUFZGQuXm5iFQB/gi8b28L0AqYbXc5dyw5Y5wNtLb7hwQRKYP15joZwBiTbow5hkuPDdYCTSVExAeUBPbjomNjjFkGHDmn+XKPRTtgkTHmiDHmKLAIaF/46c+X13iMMV8bYzLtzRVAFft5J2CmMSbNGLMTSMB6v3Pde546X3Ev3JWBPQHbiXabK9inI+sDK4GrjDH7wSruQJzdLdTHOBZ4Gsi2t8sDxwLejALz5o7F3n/c7h8qagLJwIf2qf/3RSQKFx4bY8xe4A1gN1bBPg6swb3HJsflHouQPUZ56I111gCKxnjUBRT3wp3XJwJXfD9OREoBc4AhxpgTF+uaR1tIjFFEOgIHjTFrApvz6GouYV8o8GGdypxojKkPnMI6HXshITsee+63E9Zp1kpAFHBnHl3dcmx+y4Xyu2JcIvIckAl8lNOURzfXjEddXHEv3IlA1YDtKsA+h7JcMhEJwyraHxlj5trNB3JOs9o/D9rtoTzGZsCfRGQX1im7VlifwMvap2fh7Ly5Y7H3R3P+qVAnJQKJxpiV9vZsrELuxmNzB7DTGJNsjMkA5gJNce+xyXG5xyKUjxFgXTwHdAR6GGNyirBrx6N+W3Ev3D8A8faVsuFYF9XMdzjTRdnzhpOBrcaYtwJ2zQdyrnjtBXwa0N7Tvmq2CXA851Sh04wxI4wxVYwx1bH+7pcaY3oA3wCd7W7njiVnjJ3t/iHzacEYkwTsEZE6dlNrYAsuPDZYp8ibiEhJ+99czlhceWwCXO6x+ApoKyIx9lmItnZbSBCR9sAzwJ+MMacDds0H7rev9q8BxAOrcOF7nsqD01fHOf3Aupp0G9aVls85necS8jbHOrW1AVhnPzpgzScuAbbbP8vZ/QV4xx7fRqCR02O4wLhux39VeU2sN5kE4N9AhN0eaW8n2PtrOp07j3HUA1bbx+cTrCuRXXlsgJeBn4BNwDSsK5Rdc2yAj7Hm5zOwPmk+ciXHAmvuOMF+PBxi40nAmrPOeS94N6D/c/Z4fgbuDGh31XuePs5/6JKnSimllIsU91PlSimllKto4VZKKaVcRAu3Ukop5SJauJVSSikX0cKtlFJKuYgWbqWUUspFtHArpZRSLvL/9RQQOUJzrGcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "scale = 10\n",
    "r0=np.array([0,0,0.5,3**0.5/2,1,0])*1000\n",
    "u=np.zeros((6))\n",
    "u[2:5]=ufree\n",
    "ix=range(-2,6,2)\n",
    "iy=range(-1,6,2)\n",
    "F=K@u\n",
    "r=r0+u*scale\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(r0[ix],r0[iy],'s-',label='orignal')\n",
    "plt.quiver(r0[ix],r0[iy],u[ix],u[iy],color=(0,1,1,1),label='displacements')\n",
    "plt.plot(r[ix],r[iy],'o-',label='deformed')\n",
    "plt.quiver(r0[ix],r0[iy],F[ix],F[iy],color=(1,0,0,1),label='applied forces')\n",
    "plt.axis([-100,1200,-100,1000])\n",
    "plt.legend(loc='center left', bbox_to_anchor=(1,0.5));\n",
    "plt.title('original and deformed structure\\nscale = {}x'.format(scale));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project Matrix creation\n",
    "\n",
    "In the project, we have 7 joints _(nodes)_ and 11 bars _(elements)_. That means we have 7 locations where $\\sum\\mathbf{F} = \\mathbf{0}$ and 11 stiffness matrices $(\\mathbf{K_1}...\\mathbf{K_{11}})$. The result is 14 equations with 14 unknowns, but we have some constraints on the system. $u_{1x}=u_{1y}=0$ and $u_{7y}=0.$ Which leaves us with 11 unknowns. There are 14 applied forces $(F_{ix},~F_{iy})$ that make up the vector $\\mathbf{F}.$\n",
    "\n",
    "$\\mathbf{u}=\\left[\\begin{array}{c}\n",
    "        u_{1x}\\\\\n",
    "        u_{1y}\\\\\n",
    "        u_{2x}\\\\\n",
    "        u_{2y}\\\\\n",
    "        u_{3x}\\\\\n",
    "        u_{3y}\\\\\n",
    "        u_{4x}\\\\\n",
    "        u_{4y}\\\\\n",
    "        u_{5x}\\\\\n",
    "        u_{5y}\\\\\n",
    "        u_{6x}\\\\\n",
    "        u_{6y}\\\\\n",
    "        u_{7x}\\\\\n",
    "        u_{7y}\\\\\\end{array}\\right]=\n",
    "        \\left[\\begin{array}{c}\n",
    "        0\\\\\n",
    "        0\\\\\n",
    "        u_{2x}\\\\\n",
    "        u_{2y}\\\\\n",
    "        u_{3x}\\\\\n",
    "        u_{3y}\\\\\n",
    "        u_{4x}\\\\\n",
    "        u_{4y}\\\\\n",
    "        u_{5x}\\\\\n",
    "        u_{5y}\\\\\n",
    "        u_{6x}\\\\\n",
    "        u_{6y}\\\\\n",
    "        u_{7x}\\\\\n",
    "        0\\\\\\end{array}\\right]$, and $\\mathbf{F} = \\mathbf{K_1 u}+...+\\mathbf{K_{11} u}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Node and element arrays\n",
    "\n",
    "In FEA, we typically store the node locations and their connections in two arrays called the node array: `nodes` and the element array: `elems`. The array `nodes` stores the node number and its x-y location. The array `elems` stores the node numbers that define each element. Our `nodes` and `elems` are printed below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "node array\n",
      "---------------\n",
      "[[  1.           0.           0.        ]\n",
      " [  2.         150.         259.80762114]\n",
      " [  3.         300.           0.        ]\n",
      " [  4.         450.         259.80762114]\n",
      " [  5.         600.           0.        ]\n",
      " [  6.         750.         259.80762114]\n",
      " [  7.         900.           0.        ]]\n",
      "element array\n",
      "---------------\n",
      "[[ 1  1  2]\n",
      " [ 2  2  3]\n",
      " [ 3  1  3]\n",
      " [ 4  2  4]\n",
      " [ 5  3  4]\n",
      " [ 6  3  5]\n",
      " [ 7  4  5]\n",
      " [ 8  4  6]\n",
      " [ 9  5  6]\n",
      " [10  5  7]\n",
      " [11  6  7]]\n"
     ]
    }
   ],
   "source": [
    "l=300 # mm\n",
    "nodes = np.array([[1,0,0],[2,0.5,3**0.5/2],[3,1,0],[4,1.5,3**0.5/2],[5,2,0],[6,2.5,3**0.5/2],[7,3,0]])\n",
    "nodes[:,1:3]*=l\n",
    "elems = np.array([[1,1,2],[2,2,3],[3,1,3],[4,2,4],[5,3,4],[6,3,5],[7,4,5],[8,4,6],[9,5,6],[10,5,7],[11,6,7]])\n",
    "print('node array\\n---------------')\n",
    "print(nodes)\n",
    "print('element array\\n---------------')\n",
    "print(elems)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the mesh in different ways, here I have chosen to loop through each triangle in the mesh, $node~1\\rightarrow node~2\\rightarrow node~3 \\rightarrow node~1 \\rightarrow node~2 \\rightarrow node~3 \\rightarrow ...$\n",
    "\n",
    "I create an array where each column is a loop around each triangle in the support structure for the x- and y-coords, `ix` and `iy`, respectively. \n",
    "\n",
    "The vector `r` is the coordinates for the structure, so $\\mathbf{r} = \\mathbf{r}_0+\\mathbf{u}$, where $\\mathbf{r}_0$ is the location of each node without forces applied and $\\mathbf{u}$ is the distance each node moves after a force is applied. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0.        ,   0.        , 150.        , 259.80762114,\n",
       "       300.        ,   0.        , 450.        , 259.80762114,\n",
       "       600.        ,   0.        , 750.        , 259.80762114,\n",
       "       900.        ,   0.        ])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ix = 2*np.block([[np.arange(0,5)],[np.arange(1,6)],[np.arange(2,7)],[np.arange(0,5)]])\n",
    "iy = ix+1\n",
    "\n",
    "r = np.block([n[1:3] for n in nodes])\n",
    "r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAE0CAYAAAACMmU0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJztnXmYFNXVuN8zM2yiLKIJCAIuqIAiQUD84q4RxS1GjFuMO8Y1xhU0RuPyqSS/mE9FEzcgLiwuUTRGVFwiRmWJiCAiiBAQUESQRQVm5vz+uHd6unu6Z6pruru6e877PPfpqntv1T11q7pO3eWcK6qKYRiGYWRKWdQCGIZhGMWJKRDDMAwjFKZADMMwjFCYAjEMwzBCYQrEMAzDCIUpEMMwDCMUpkAMwzCMUJgCMQoWETlLRD4UkW9FZKWI3C8i7bJchorIrtk8ZzYRkTEicmsOz/+GiJyXq/MbpY0pEKMgEZErgTuBq4G2wCCgG/CKiDQPcb6KkHKEOi5fRC1f1OUbEaOqFiwUVADaABuAnyfFbw18CZzj98cAt8alHwwsi9tfDFwLzAY2ARVJ5/sXoMBGX97JNefwx60EHgXOAqYmHavArn57CPARsB74HLjKx28HvACsBb4G3gLKUlyvAHf5a/vGy7snMAzYAmz28j2f7rri5UlTN8cDs4B1wKfAkcBtQBXwvT//vUB3f66KuGPfAM7z22cBb3t5v64pAzgHmAesASYD3aJ+jizkPtjXg1GI/A/QEngmPlJVN4jIP4GfAI8EPNepwNHAV6pamXS+A0VEgb1VdSGAiBwMdAS2xbV4ynCKpT4exim7t0SkPbCTj78Sp4y29/uDcC/nZI4ADgR2wymQPYC1qvqAiPwPTin+tr7rEpG0wonIQOBvwFBgCtAJ2EZVXxKRHwOPqepDPm/3Bq4VYF9gPPADoJmI/BS4DjgWWAAMB8bh7qNRwlgXllGIbEeKF75nhU8Pyt2qulRVv8vgmGrgRlXdFPC4LUAvEWmjqmtU9T9x8Z1wX+NbVPUtVU2lQLYA2+AUh6jqPFVd0UCZmVzXucAjqvqKqlar6ueq+nGA49KxXFXvUdVKX/4FwO1e7krgf4G+ItKtEWUYRYApEKMQ+QrYLk3/eiefHpSlIcpfparfZ5D/RFw31hIReVNE9vPxfwAWAi+LyCIRGZ7qYFV9Ddd9NAr4QkQeEJE2DZSZyXXtiOu2yhbJZXcD/k9E1opITXedAJ2zWKZRgJgCMQqRd3B9+z+LjxSR1sBRuG4YcGMXW8Vl6ZjiXGHcTScfk1COiCSUo6rTVfV4XJfOs8BEH79eVa9U1Z1x3TtXiMhhKQtUvVtV9wF647qyrm5A/uT4b0lfF0uBXQKeZ6P/ra9ek49ZClygqu3iQitV/XeaMo0SwRSIUXCo6jfA74F7RORIEWnm++afxI0pPOqzzgKGiMi2/qV+eYjivgB2biDPB0BvEekrIi2Bm2oSRKS5iJwuIm1VdQtukLrKpx0jIruKG6Coia9KPrmIDBCRfUWkGe4F/n1cviDygauL00SkXESOBA6KS3sYOFtEDhORMhHpLCJ7pDq/qq7CTQT4hT/XOaRXPjX8BRghIr399bQVkZMCyGwUOaZAjIJEVUfiBmb/iHv5vof70j1MVTf5bI/iXu6LgZeBCSGKugkY67tffp5Glk+Am4FXcYPEU5OynAEsFpF1wK+AX/j4Hv6YDbhW1X2q+kaKItoAD+JmMC0BVuOuG9zLv5eX79l6ruPXuFbOWuB0XEuoRv5pwNm4mVPfAG/iup0A/g8YKiJrRORuH3c+rgW0Gtciqrcloap/x025Hu/rYA6upWiUOJJ6TM8wDMMw6sdaIIZhGEYoTIEYhmEYoTAFYhiGYYTCFIhhGIYRClMgRgwRuUlEHsvSuXLqRTZXiMh1IvJQ1HIERUQOFpFlUctRHyKyWEQOT5MW2XNSn1xGMEyBGE0GEenu3ben9QGnqv+rqhm5NxeRFiLysIgsEZH1IvK+iNg01gKisYrKLy2QPH27yWMKxDAaTwXORuUgnOv5G4CJAR0TGkbRYgqkSBGRHUTkaRFZJSKfichlcWkDReQdb3y2QkTujV9DQ0R6i8grIvK1iHwhItfFnbq5iPzNf0nPFZH+9ciwR9x55qczxPN5jxGRWV6mf4tIn7i0xSJytYjMFpGN/mv+hyLyTy/Hq+K83NbkH+TPsVZEPhDnQbcm7Q0RuUVE3vbHviwiNc4X/+V/14rIBqn1WRUvZ6wbT0RaishjIrLalzVdRH6YfIyqblTVm1R1sXdW+ALwGbBPmro4S0SmisgfvQHfZ/EtFn9vJ/l6XSgi58eltfJf02tE5CNgQNK5G3ouZojIOn/f/5RGvvYi8oI/xxq/3SVgHSMiZ/jW2GoRuT5VGekI8Jxc5Z+Tb0RkgjjPADXp1/jnfbmInCd+sTARGYYzrrzG3/fn44rsm+58RgCi9idvIfOAU/wzgd8BzXGuKBYBg336PjjX4RW49R3mAZf7tG1wHm2vxLlM3wbY16fdhHOjMQQoB24H3k0jQ2vcV/fZvpx+OCeHvX36GGrXiuiHW+tiX3/eM3HW4y18+mLgXeCHOAd8XwL/AX4EtABew3nHxaev9jKW4Vy7rwa29+lv4BwH7ga08vt3+LTuJK11keK6bsK5NwfnZfZ5nF+ocl+vbQLcnx/6etwjTfpZOA+85/vzXggsp9aw903gPn9/+gKrcBb4AHfg1hXZFuckcQ5+DZQAz8U7wBl+e2tgUBr5OuAcRG7ln48ngWfj0uur4144y/sD/b37E1AJHJ6mrEyfk2nADv765wG/8mlH4tZv6e3lfpTENVti5cSVnfZ8ae7Z1FRpTTlYC6Q4GYB7Yd6sqptVdRHOFcYpAKo6U1XfVeduezHwV2p9Ix0DrFTV/6eq36tz+Pde3LmnquqLqlqF+xPunUaGY4DFqjral/Mf4GncmhPJnA/8VVXfU9UqVR2Lc5Y4KC7PPar6hap+jntBvqeq76tzW/J3nDIB5ybkRS9jtaq+AszAKZQaRqvqJ+pcjU/EvYTDsAX3Mt3Vyz1TVdfVd4A4f1aPA2O1fpfpS1T1QV/PY3Fehn8oIjsC+wPX+vszC3gI5y4F4OfAbar6taouBe6OO2e9z4W/nl1FZDtV3aCq76YSTFVXq+rTqvqtqq7HLTx1UFK2dHU8FHhBVf/l790NOPf4QQjynNytqstV9Wuccq8p9+deprmq+i3Ol1oQ0p3PCIApkOKkG7CDb+bXuNC+Dvfli4js5rsdVorzTfS/1K6h0ZBr75Vx298CLSX1oHM3YN8kGU4ntUfcbsCVSXl3xH351fBF3PZ3Kfa3jjvXSUnn2h/3Ak53DVsTjkdxq+uN990iI72CSImIlPljNgOXNHDumIz+hYeXcwfga//irmEJta7RdyDRnfqSuO16nwvcuiC7AR/77rhj0lzHViLyV98NtQ7X9ddORMpTyU9iHSfIp6obcS3EIAR5TgKVS3B399l6VpoktiJhcbIU+ExVe6RJvx94HzhVVdeLyOXUtgyW4lazy4YMb6rqTwLmvU1Vb8tSuY+q6vkN5qxLRo7f1HnX/T3we3ED4i8C83EODhMQEfHxPwSG+GPDsBzYVkS2iVMiXXEecsF1P+4IzI1Lq6He50JVFwCnekX3M+ApEengX/LxXAnsjuvaXCkifXHPU/plD2tZAfSs2RGRrXCtuCA05jlZAXSJ298xKd2c/uUAa4EUJ9OAdSJyrR9ULReRPUWkZkB1G5wH2w3i3HZfGHfsC0BHEblc3PTTbURk3xAyvADs5gdMm/kwQER6psj7IPArcS7LRURai8jRIrJNiHIfA44VkcH+uluKs4Xo0uCRbiyhmmDu0RGRQ0RkL//lvQ7XBVTHHbvnftyL81jNbPXDBHy31L+B2/219cG1HB73WSbiXKe399d8adzh9T4XIvILEdleVatxXntJcz3b4Fp9a0VkW+DGDC7hKeAYEdlf3MSNmwn+nmnMczIR57K+p1dav0tKD+oWvz7E35NYaOT5ih5TIEWI7zc/Ftdf+xlu8Poh3BRSgKuA04D1uD/lhLhj1+MGno/FNd8XAIeEkGE9bi3vU3BfzStxLr1bpMg7A9e/fS/OZflC3KBkxvgX7PG4rplVuK/WqwnwLPuuotuAt30XyaAGDumIeyGuww2wvolTYAmIW7r1Atz9WOln+mwQkdMDX1gip+IG/Jfjxn9u9GM94FpES3D3/WVq10YJ8lwcCcwVkQ04N+6naOqVF/+MGxz/Cje54aWggqvqXOBi4Alcq2ANbg2XIMeGfk5U9Z+48aDX/XHv+KQa1/9B3eLXx//gFGsspOnebTKYO3fDMEoO3xKeg5vBVRm1PKWKtUAMwygJROQEcStEtse1hp835ZFbTIEYhlEqXIDr1vwUN7ZzYf3ZjcZiXViGYRhGKKwFYhiGYYTCFIiRFhHp6mcTlTecu9FlqYjsmutyMkUCePA18oM432hnRi2HUYspkCaKBFgLQVX/q6pb++mhhhGabHwgqOpR3r1JkPLeEJGM3PIbmWMKxEiJfXFntw68YVxO/m+lcK9K4RqaIqZAmiAi8ijOBcbzvovqmriumnNF5L/Aa8ndNyJytojME+fCe5GIXBB3zoNFZJmIXCkiX4pzq312XHoHEXlenCvx6SJyq6RZoMdbyP9RRP4rzu34X0SkVZq8Z4lzK36XNxJbJCL/4+OXelnOjMt/tLgFn9b59Jvi0urUQYryTvSttz39fkOu5W8TkbdxfpZ29nIt8nX4WTpjQ2nYJb+KyMUisgBnDJqpe/1tRWS0OB9fa+KN60TkfHFu5L8W51Z+h7g0FZFficgCf9woERGftquIvCnONfpXIjLBx9e40f/AP28nxz0v14rISmC0BHMjf17cfU/pEl9EbgMOAO715d0rjrv88/CNOBfue6arHyMguXDxa6HwA86V9eFx+91x/oL+hnPV3ook9+fA0cAuOJ9IB+Feiv182sE4t903A81w3nG/Bdr79PE+bIVz+b2UOPfYJLre/jMwCediexucl9Tb01zHWb7cs3EuwG8F/guMwlnFH4GzyN86Ts69cB9PfXAuLn4apA58GQvj5AziWv6/OBfjFTiL8HXA7j69E979fYrrSuuSP66+XvF11IoG3OunOP8/cB4K2vv7dZCPP9Qf18/X3z3Av5LKfQFoh/sIWQUc6dPGAdf7umgJ7J/q/iY9LzXeC1oRzI38eXH3vT6X+LG8fn8wztV9O9zz2xPoFPX/sNhD5AJYiOjGp1cgO6eIS7l+BvAs8Gu/fTDOvUNFXPqX/iVY7v/su8el3UoKBeL/3BuBXeLS9sM5CUwlw1nAgrj9vfy5fhgXtxrom+b4PwN3BaiDq4CPgC5xadfiHDvGn28ycKbffgO4OS6tNc4H1YlAqwzv1+XA35Pq69C4/ZOBt5KO+St+HZWk+E44n2DtU6Q9DIyM29/a37vuceXGK4aJwHC//Tfggfg6Sr6/cfsH47wWt6znmvsCa+L23yBRgSyMS9vKl9ExOa/fPxT4xD+PZVH970otWBeWkUxaN9gicpSIvOu7Ntbivry3i8uyWhMtf2vcY29P7bKvDZWzPe5lMFNqXXq/5OPTkez6HVVN6Q5enKO+1303yTfAr5KuIZ1sVwOjVDXer1MQ1/LJrs1P9mWuEJF/iHN2WQep3yV/Kjkzca+/I85l/JoUaTsQ5yJeVTfgFHDnuDzpXKBfg/sAmCZuNctzUl1bHKs0zheXBHMjH086l/h1UNXXcD62RgFfiMgDItKmAfmMBjAF0nRJZ0GaMl5EWuAWjPoj7uu+Hc69eRAX36tw3RX1uduu4SvcC7+3qrbzoa2qZmudhidw3WM7qmpb4C/UvYZUdXAE8FsROTEursa1fLu40FpV70h3LlWdrM4FfifgY5yzy1Tc79N7qGobnPPI+uSsca8fL8vWqprKGnspzmV8uxRpy3HKCAARaY3rWvo8Rd5EYVRXqur5qroDzir8Pql/5lVyPce7kW+DW9UQgj1jDZ0bVb1bVffBdSnuhvsoMBqBKZCmS6burZvj+qpXAZV+wPKIIAeqmwb8DHCT/8rcA/hlmrzVuJfqXSLyAwAR6SwigzOQtT62wX19fy8iA3Fei4MwF+fNdpSIHOfjMnItL26d9+P8S3kTbunXdFOk63PJn4rA7vVVdQXwT9wLvr3PW/OyfgLnFr2v/2j4X9zqkIsbKB8ROSnu2tfgXuI11xfkeWuMG/lkEsrzdbGvuAXBNuKWHLbp6Y3EFEjT5XbcF/VaEbmqoczq3LdfhuvzXoN78U7KoLxLcIPIK3EuyMdR62o7mWtxg9Xv+q6MV3FfptngIuBmEVmPWzNiYtADVfUD3FK+D4rIUZq5a/ky3Ff2cuBr3ESEi9LkTeuSP41sgd3re87AjW18jBurutyfZwpuGdqnce7Yd6F2SdyGGAC8J85d/CTc+NhnPu0mYKx/3tLNDgvtRj4F/wcM9TO07gba4OpxDa6LbjWuNW00AvOFZUSCiNyJG/A0y2LDKFKsBWLkBXE2Cn38fPyBuFX2/h61XIZhhMesP418sQ2u22oHXJfJ/wOei1QiwzAahXVhGYZhGKGwLizDMAwjFCXdhbXddttp9+7doxbDMAyjqJg5c+ZXqlqf8S5Q4gqke/fuzJgxI2oxDMMwigoRWdJwLuvCMgzDMEJiCsQwDMMIhSkQwzAMIxSmQAzDMIxQmAIxDMMwQmEKxDAMwwiFKRDDMAwjFKZADMMwjFCYAjEMwzBCYQrEMAzDCEXkCsQvBfq+iLzg93cSkfdEZIGITBCR5j6+hd9f6NO7Rym3YRhGUydyBQL8GpgXt38ncJeq9sAtP3mujz8XWKOquwJ3+XyGYRhGRESqQESkC3A08JDfF+BQ4CmfZSzwU799vN/Hpx/m8xuGYRgREHUL5M/ANUC13+8ArFXVSr+/DOjstzsDSwF8+jc+fwIiMkxEZojIjFWrVuVSdsMwjCZNZApERI4BvlTVmfHRKbJqgLTaCNUHVLW/qvbffvsG3dkbhmEYIYlyPZAfA8eJyBCgJdAG1yJpJyIVvpXRBVju8y8DdgSWiUgF0Bb4Ov9iG4ZhGBBhC0RVR6hqF1XtDpwCvKaqpwOvA0N9tjOB5/z2JL+PT39NbUF3wzCMyIh6DCQV1wJXiMhC3BjHwz7+YaCDj78CGB6RfIZhGAYFsqStqr4BvOG3FwEDU+T5Hjgpr4IZhmEYaSnEFohhGIZRBJgCMQzDMEJhCsQwDMMIhSkQwzAMIxSmQAzDMIxQmAIxDMMwQmEKxDAMwwiFKRDDMAwjFKZADMMwjFCYAjEMwzBCYQrEMAzDCIUpEMMwDCMUpkAMwzCMUJgCMQzDMEJhCsQwDMMIhSkQwzAMIxSmQAzDMIxQmAIxDMMwQmEKxDAMwwiFKRDDMAwjFJEpEBFpKSLTROQDEZkrIr/38TuJyHsiskBEJohIcx/fwu8v9Ondo5LdMAzDiLYFsgk4VFX3BvoCR4rIIOBO4C5V7QGsAc71+c8F1qjqrsBdPp9hGIYREZEpEHVs8LvNfFDgUOApHz8W+KnfPt7v49MPExHJk7iGYRhGEpGOgYhIuYjMAr4EXgE+BdaqaqXPsgzo7Lc7A0sBfPo3QIcU5xwmIjNEZMaqVatyfQmGYRhNlkgViKpWqWpfoAswEOiZKpv/TdXa0DoRqg+oan9V7b/99ttnT1jDMAwjgYKYhaWqa4E3gEFAOxGp8EldgOV+exmwI4BPbwt8nV9JDcMwjBqinIW1vYi089utgMOBecDrwFCf7UzgOb89ye/j019T1TotEMMwDCM/RNkC6QS8LiKzgenAK6r6AnAtcIWILMSNcTzs8z8MdPDxVwDDI5C5KPnTn6BXL+jTBw47DJYsiVqi0uCpp0AEZsyIWpLiZ+JE94z27g2nnRa1NEZQKhrOkhtUdTbwoxTxi3DjIcnx3wMn5UG0kuNHP3Ivua22gvvvh2uugQkTopaquFm/Hu6+G/bdN2pJip8FC+D22+Htt6F9e/jyy6glMoJSEGMgRjgWL4aePeH8892X2xFHwHff1c13yCFOeQAMGgTLluVVzKIhaH0C3HCDU8QtW+ZVxKIiaH0++CBcfLFTHgA/+EFexTQagSmQImfBAvfnmzsX2rWDp5+uP//DD8NRR+VHtmIkSH2+/z4sXQrHHJN/+YqNIPX5yScu/PjH7gPnpZfyL6cRjsi6sIzssNNO0Lev295nH/fVl47HHnNdWW++mRfRipKG6rO6Gn7zGxgzJt+SFSdBns/KSqdo3njDtY4POADmzHEKxyhsrAVS5LRoUbtdXu7+jKl49VW47TaYNCnxmCAcfjiIaCwcfnh4eQudhupz/Xr3cjv4YOjeHd59F447LrOB9N694+vT7ZcqQZ7PLl3g+OOhWTOncHbf3SmUoFx0EVRUuAkNFRVu38gTqlqyYZ999tFS5rPPVHv3rt3/wx9Ub7yxbr7//Ed1551VP/kk8zIOO0wVqhU0LlQrvKQ4Q84SCt0UPozbv1LhxgaOeV1hnwzKmJWmPmcVwPVHVZ+DFcb47Q4K/1XYNmAZ96SoT9ULLwz7rzJUVYEZGuAday2QJsDVV8OGDXDSSa474bjjgh87ZQrUdQIgwBHZE7BJ0YfU9dknAlkKhcnAamAuzgzsaoLbCP+KVE4qHnggW7IZ9SFawrZ4/fv31xk2Sb9RiCjpvcjY90fmVGP1mU3S1adrixjhEJGZqtq/oXz2xBqhCdLEzWro1g1dtSp9+rHHor17R951mi4srm+GA1I49fnEE+iee6J77YUOHlx/nUcYbr31VqAqZW2Wl+fiiTeSMQVipOWtt94CXsZ9HcejbLPNexFIVA/PPANbbx21FPXy85//HJhNqvrs2PGrCCRKQWUl/PrX8PrrMHu2c19w771RS1WHTZs2ccsttwB/IVV9nnPOlgikanqYAjFSUl1dzZVXXgkcSa0SqQkvs379fkydOjU3hT/2GAwc6AZsLrgAqlJ/ZcbYsMH5a/ntb3MjTxaYOnUq06ZNw62dVqNEasJs1q/vTlVD1xmWTOpT/Tj0xo3ud9062GGH3MjVCIYNG8amTZuAS4FRlJXV1GUlMIquXW29uXxgCsRIycSJE5k+fToALVocz+LF/2XLlipat94Gp1TgjDPOyH7B8+Y5Pytvvw2zZrm+iMcfr/+YG26AK6+sNbcvQH7xi1/Etlu3/jFVVcpbb72N+wv2ZePGjVxxxRXZLzjT+mzWzPm72Wsvpzg++gjOPTd9/ghYuXIljz32WGz/0EP/TlWV8OCDD+PWpbuUO+64g5UrV0YmY1PBFIhRh++//57hw2t9VV5++eV069aNiooK7rjjjlj84sWLGTt2bKpThGfKFJg5EwYMcF/MU6bAokXp88+aBQsXwgknZFeOLDJmzBiWxHmwHDlyJGVlZey///4MHFjr9m3UqFF888032S080/rcssUpkPffh+XLXRfW7bdnV6ZGcvLJJ1NdXQ1AWVkZ48ePB+Dss89mzz33BGDjxo3cdNNNUYnYdIh6ICyXodTtQHLFyJEjY/0r2223na5duzYhvVOnTrH0tm3balVVVfYKv/tu1eHDU6d166a6alVi3H33qXbq5NI6d1Zt1kz1oIOyJ08jqaqq0jZt2sTqa4cddkhIX7x4sYqb6qaAHn300dkVINP6nDZN9dBDa/fffFP1qKOyK1MjmDFjRoIdyLnnnpuQ/tJLtfZJZWVlOmfOnIgkLW4IaAcS+Us+l8EUSOasWrVK27ZtG/sT3nvvvXXyxP9JAR0xYkT2BJg7V3XXXVW/+MLtr16tunix2071wosn2bKyABg+fHhCXb388st18gwdOjQhzydhLD7TkWl9fv65aseOql9+6fZ/+1vVK67InjyNpEePHrF6atmypW7atKlOniOOOCKWZ8iQIRFIWfyYAjEFEorLLrss9ufbbbfddPPmzSnz9enTJ5avWbNmunHjxuwJMX686t57q+61l2q/fqrvvOPii0yBbNy4UZs1axarp759+6bMt2HDhkD5QpNpfd5/v+oee7j8xxyj+tVX2ZUnJE8++WSCor3zzjtT5vvggw8SWnWvvPJKniUtfkyBmALJmE8++UQrKipif7xnn302bd6PP/444c980kkn5VHS4iC5ZTF//vy0eYO0VJoy1dXV2qFDh1j9bL/99vXmP+ecc2J59957b62srMyTpKWBKRBTIBlzwgknxP50Bx54oFZXV9ebf8iQIbH8IqJLlizJk6SFT6ZjGw2NlTR1brnllgQFW9/Hjarq559/rltttVUs/+jRo/MjaIlgCsQUSEb861//SviDTp8+vcFj1qxZo+Xl5bFj9t133zxIWhwMHDgwVi8VFRW6Zs2aBo8ZPXp0wj0YNWpUHiQtfL777jtt0aJFrF722GOPQMfdeOONCQo5q92sJY4pEFMggamqqtIBAwbE/mynn3564GMvueSShJfeW2+9lUNJi4O33noroU4uu+yywMd269Ytdlzr1q2t60VVf/nLXybU5wcffBDouPXr12vHjh1jx91yyy05lrR0MAViCiQw48aNi/3JWrRooYtrZukEYMuWLdq6devY8d27d8+hpMVBshLIZJpzY5RPKbJixQotqzUz10PjpxgH4MEHH0y4FytWrMiRpKVFwSsQYEec7+Z5OD/Ov/bx2wKvAAv8b3sfL8DdwEKcL4h+DZVhCqRhvvvuu4QX3rXXXpvxOe65556El96YMWNyIGlx8MgjjzS6Gyq++6u8vDxQ91epcuCBB8bqoqysTFfVNwsvBZWVlbrnnnvGzjFs2LAcSVpaFIMC6VSjBIBtgE+AXsBIYLiPHw7c6beHAP/0imQQ8F5DZZgCaZiGjAaDklPjwiIhWwPhOTcuLBKmT5+eoIyTjQaDYsaFmVPwCqSOIPAc8BNgPtBJa5XMfL/9V+DUuPyxfOmCKZD6CWI0GJRk48Lh6ayfS5hsTsXNZApwqRLEaDAoZlyYGUWlQIDuwH+BNsDapLQ1/vcFYP+4+ClA/xTnGgbMAGZ07do1u7VaYgQ1GgxKTo0LC5xsGwNUM6fWAAAgAElEQVQGNUIsVSZOnJigQNMZDQbFjAszo2gUCLA1MBP4md9Pp0D+kUKB7FPfua0Fkp5MjAaDkmxcOHTo0CxIWhzkwh1JUzUurK6u1m23rV0TvSGjwaCYcWFwikKB4HwvTwauiIuzLqw8EG80eNBBBzVoNBiUpmhcmKsxi6ZqXJip0WBQzLgwOAWvQHCD4X8D/pwU/wcSB9FH+u2jSRxEn9ZQGaZAUhPGaDAoTdG4MIzRYFDGjBnT6FldxURYo8GgmHFhMIpBgezvb+RsYJYPQ4AOvntqgf/dVmsVzijgU+DDVOMfycEUSF0aYzQYlKZkXJgPu42mZFwY1mgwKGZcGIyCVyD5CKZA6tIYo8GgVFZWJhgXduvWLetlFAqNMRoMSrKSuvTSS7NeRiHQWKPBoJhxYcOYAjEFUodsGA0GpSkYF2bDaDAo++67b6ycUjUubKzRYFDMuLBhTIGYAqlDtowGg1LKxoX5HuAudePCbBkNBsWMC+vHFIgpkASyaTQYlFI2LkyeYpsPu4JSNi7MptFgUMy4MD2mQEyBJJBto8Gg7L333rFyS8W4MOcrCKYh2bhw7733zku5uSbbRoNBMePC9JgCMQUSY/78+Vk3GgxKKRoXnnjiiQnXlNU1zBug1IwLq6qqcmI0GJRzzz03VnafPn1KeoZbJpgCMQUSI1dGg0EpJePCqMciSs24MFdGg0Ex48LUmAIxBaKquTUaDEopGRfm0mgwKMnGhfkYz8oFuTYaDIoZF9bFFIgpkLwYDQalFIwLC2mxp2T7ky1btkQmS1jOOOOMhPrMttFgUJKNC2+++eZI5CgkTIGYAtEnnngi9qfIldFgUErBuDAfRoNBKXbjwmSjwcMOOyxSecy4MBFTIE1cgSQbDRbCFNpiNi7Mp9FgUIrZuDBfRoNBMePCREyBNHEFkm+jwaAUo3FhoQ5cL1mypCiNC/NtNBgUMy6sxRRIE1YgURgNBqUYjQuvvfbaBJkLyV4g2bjw448/jlqkBonCaDAo8caFRx11VNTiRIYpkCasQKIyGgxKMRkXRmU0GJRiMy5MNhocOXJk1CIlkGxcWOx2NmExBdJEFUiURoNBKSbjwiiNBoNSLMaFURsNBsWMC7OsQID2QG9gZ6AsyDGFEJqiAonaaDAoxWBcGLXRYFCqqqoSuiwLZYwmmaiNBoNixoVZUCBAW+A63OJN84GpwAxgKfAkcEiQAqIMTU2BFILRYFCSjQsHDhwYtUh1SDYaLJSJCKlINi685557ohYpgWSjwZ49e0YtUr0kGxdu2LAhapHySjYUyCvAGUC7FGn7AH8Gzg1SSFShJBVIt26qKaY8Vo8YoSubNdP1RG80qKqp5dy4UXXIENXdd1ft1Utf7tcv4aUXiXFhmvpc27evfgz6vg/XnXde/mWLJ5Wc69ap7r13LKwuK9O74mwZIjEuTFOf9x5wgH4AOgf0TqIzGoyRRk697jrVLl20unXrBOPCMeeco/qjH6mWl6s++WTexc03NgbSxBTI5N//XjuCrid6o0FVTa9AXnvNbW/apNX7768/jfsqjcS4ME19vtOihe4T9zKOfLpxuhdeHOt3200PiFPIkRgXppBz5dy5ugR0Oy/X5I4dVV99Nf+yxZOuPt95R3X5ctXWrfWhhx6qbTG1aqWrpkxRPeMMUyBxoYwAiEgfETlORH5WE4IcZzSCxx6DgQOhb1+44AKoqkqb9fvvv2fYI4+w0u//5je/oVu3bgUnJ1ttBYcc4rabN0f69ePSE06IJS9ZsoSxY8dGLufo0aP5ftOm2P7IkSMpKwv0V2k8mdRnPAsWsPW337J54MBY1H333cfatWsjl/N3v/gFnwBfAWVlZex3ww3w9NO5kasRcgIwaBB06gTAWWedxZ577gnAvO++4/oJEyBfz0GR0GBtiMgjwCPAicCxPhyTY7maNvPmwYQJ8PbbMGsWlJfD44+nzX7PPfewZMkSAAQYPnx4QcqZwNq18PzzHHrbbXTyf1iAyy67jOrq6sjkrK6u5vLLLwdgNDCnooKLVq8G13WbWxpTn+PGwcknM/HJJxERAKqqqjj99NMjlXPGjBlMfP999gC6AeeeeSbbTJkCS5dmX65GyJmK8vJy/vjHP8b2H3roIdasWZMLSYuWIOp0kKr2V9UzVfVsH87JRuEi8oiIfCkic+LithWRV0Rkgf9t7+NFRO4WkYUiMltE+mVDhoJkyhSYORMGDHBfTlOmwKJFKbN+9dVX3HbbbbH95s2b07Zt24KTM4HKSjj1VLjsMth5Z8aMGRNLWrduHddff31kcl533XWsW7eO04E+wKpnnoG33oJHH82uTI2Usw7jx8Opp9K1a1dOPPHEWPSLL77I/PnzI5PztNNOYy1wIfCkCH/56CPo3h0qKrIrUyPlTMfgwYM54ogjAPdxMWPmzFxIWrw01McFPAz0CtIflmkADgT6AXPi4kYCw/32cOBOvz0E+CfuI3sQ8F5D5y/aMZC771ZNZ6Gd1Hd76aWXxvppd999d61u3To/MqpmJGcCZ5+tmtQ/n1PjwoByrl+/PrXR4OjRqhdfnD15GilnHWbNUu3RI7abc+PCgHJOmDAhtdHgX/+qevXV2ZWpEXKmJO5/NHv27Jjjx9Ggs3772+zKWYCQrUF0/5L/BjeVdzZuWu/sICcPJAB0T1Ig84FOfrsTMN9v/xU4NVW+dKFoFcjcuaq77qr6xRduf/Vq1ZpB8bgHP9lo8Lnnnkt48AtFzgSuv171Zz9TTRqUzqlxYUA5a4wGy0E71BgNbt6seuKJqvffnz15GilnHa69VvV3v0uIGjFiRO6MCwPImWw02LNDB5f+9ddu1lg+1nMPW5+qdf5HNcaFo0Gv6Nq15I0Ls6lAFgLHATvhujG7Ad2CnDyQAHUVyNqk9DX+9wVg/7j4KUD/FOcbhrNXmdG1a9fs12y+GD/e/dH22ku1Xz83O0Q14cGPNxp8oksXre7cWVVEtXNn1RtvLBg5Yyxd6h65PfaonX764IOx5GTjwqzOJGtAznijwa1AF7Rp4/L26qV62WWq+XphZFKfNey0k+q8eQlRycaFnTp1yqucyUaDSw84QLVnTxfGjcuuLI2Qsw5XX+3+P0n/oy//8Q9dJqIbQL8C/bpAjTWzRTYVyGtBThQ2ZKBA/pFCgexT37mLtgUSgGIyGgxKlMaFyUaD33zzTd7KzhVRGRcWm9FgUJqScWFQBRJkEP1jEXlCRE7N0zTeL0SkE4D//dLHLwN2jMvXBVieQzkKlurqaq688srY/umnn07//v0jlCg7tGvXjgsvvDC2P23aNKZOnZrzcqdOncq0adNi+xdddBFt2rTJebm55swzz0yYzj18+HAqKytzXu6wYcPYFDcNevz48TkvMx9cddVVdOzYEYDly5fzpz/9KWKJCoCGNAxuNmNyeCSIdgoSqNsC+QOJg+gj/fbRJA6iT2vo3KXaAimklQazTVVVVd5XLiyklQazzdSpUxNaIbk2Liy0lQazTbxxYSmvXEgxWKID44AVwBZcC+NcoAOue2qB/93W5xVgFPApbiC/zvhHcihFBVKIKw1mm+SVC3PpzC55pcH77rsvZ2VFRT5XLiy0lQazTfLKheeff37UIuWErCkQ3OD5n4BngEk1IcjJow6lqEAKdaXBbBO/cmGbNm1y0ipIXmmwc+fOWS+jEEheuXDIkCE5KSd5pcHzovYfliOSVy788MMPoxYp62RTgXwAXAYcAhxUE4KcPOpQagqkkFcazDaTJ09OeBnloqV1zTXXJJTxatT+mXJIPlYujF9psFWrVgW10mC2KfWVC7OpQBo02CvUUGoKJNlosNBWGsw2uTQuTGs0WKLk2rgwrdFgiRJvXEi27WwKgGwqkNOAG4H9cFbj/YB+QU4edSglBZLSaLDEyaVxYfxKgyJSkCsNZptcGRcWy0qD2aaUVy7MpgK53Q9wvwm87kNObUOyFUpJgRTLSoPZJhfGhckrDR5zzDFZkLTwyZVxYbLRYFP4uFGtu3LhI488ErVIWSObCuRjoHmQkxVaKBUFUopGg0HJhXFhKRoNBmXs2LEJz1JjjQtL1WgwKKVqXJhNBTIB+EGQkxVaKAUFUlVVpQMGDIg9pJGvNBgB8WM/0LiVC996662Ec1122WVZlLQ46N69e+z6G7ty4RlnnJFQn5GvNJhnNmzYkDBj8Oabb45apKyQTQXyBvA1MBmbxpt3ko0GlyxZErVIeSebxoWlbDQYlGwZF5a60WBQStG4MJsK5KBUIcjJow7FrkCagtFgUO69996El16Y/uamYDQYlGwYF5a60WBQStG4sNEKBJAGDw6QJ8pQ7AqkqRgNBmWHHXaI1UemxoVNxWgwKI01LmwqRoNBKTXjwmwokDeAS4GuSfHNgUOBscBZQQqJKhSzAmlKRoNBefnllxNeWpm0yJqS0WBQTjrppIQ6ycS4sCkZDQZl8ODBsTopduPCbCiQlsBFwNs4r7cfAYuAJcCDQN8gBUQZilmBNDWjwaCEMS5sakaDQQlrXNjUjAaDUkrGhVkbA3HnohludcB2QfIXSihWBdIUjQaDEsa4sCkaDQYl2bhw8uTJ9eZvqkaDQSkV48KsKpBiDcWqQJqq0WBQjj766ASFUJ9xYVM1GgxKpsaFN998c4LCsY+bRErFuNAUSAEokL59++qFF16oy5cvT59p+XLVAw9U9VP/3nzzzYQ/aFMyGgxKJsaF8TY0Tc1oMChBjQubutFgUNIaFyb91wsZUyAFoEAAbd68ubZs2TK9IrnwQtWyMtWLLjKjwQwIYlxoRoPBCWJcmGw0OHv27AgkLXzSGhfG/dcLnawpEOASoH2QkxVaKAQFUhNSKpLly1VbtnS3oVUrfea++2L5m6rRYFCCGBea0WBwko0LL7nkkoR0MxrMjGTjwpXvv5/wXy/0Vkg2FcitwEJgInAkBW77ER8KSYGkVCS//KVq8+aqoNXNm+vYrbcONUW1qVKfcaEZDWZOfcaFBxxwQCytKRsNBiXZuPD1nj1j/3Vt3rzgWyHZnoUlwGBgvFcm/wvsEuTYKEMhKpAERQJ6IehyUAXdCPpDzGgwE1IZF5rRYDjSGRea0WA4aowLO4J+6//jsVDgrZCgCqSMAPgTrvShEmgPPCUiI4Mcb9Rl8+bNfA88DOyMM7j5ArgBuOmmm2jbtm2U4hUNY8aMiW2vW7eO66+/nhEjRrBu3bpY/NixYyOQrPjo2rUrQ4cOje2/+OKLzJ8/n9NOOy0W16pVK0aNGhWFeEXH4MGDGTx4MDfgvsATqKqCW26JQKos05CGwS1nOxPnTPEkoJmPLwM+DaKlogqF3AJJDs1BW4KeD7rEZl5lRLxxYUVFhRkNNoJk48L4wWAwo8FM+WjKlLqtjyJohZDFFsh2wM9UdbCqPqmqW7ziqQaOaaT+yhgROVJE5ovIQhEZnu/yc8Vm4Hucf5jdBw3ioosuYsWKFRFLVRxMmDDBb82isnIzW7ZsAqqBD5g4cWKEkhUfW221FVdddZXfm8WKFZ/j6rIakTlcffXVEUpXfPR86ikqytK8ZkugFSJO2RQHIlIOfAL8BLdK4nTgVFX9KFX+/v3764wZM/IoYSIidRqugSkTARG6devGHnvsQatWrbIoWenxzDM3A71I7CxQtt56CUcccWVEUhUn1dXVPPvsTUAfrD7D0+6777hv8mRaVFenz9SqFSxaBB075k+wAIjITFXt32C+IlMg+wE3qepgvz8CQFVvT5W/mBWIkSnVpOhpxvW8BBrqMxKw+mwso4BzcE4F09K8OZx3HhTYuFJQBVJsT0JnYGnc/jIfF0NEhonIDBGZsWrVqrwKZxiGAdAROJsGlAfA5s0wejSsXJl7oXJARdQCZEi6T6LaHdUHgAfAtUDyIVQuqaiooKysjEMOOYShQ4fSvn37rJdx5IUX8tqdd7K5TZuE+B/feist16yhrKqKr3r25P3zzoPy8qyXnw3iJg/V4amnnsqfIKSvT9myhR89/DDbzZ0LIsw57TSWDxqUV9mCUgz12fuJJ+j65ps037iR5x57LBZftmUL/e+5h/aLFrF5661574or+PYHP8irzH0feIDmr70GlZUNZ64ZCymwVkgggoy0F0oA9gMmx+2PAEaky19Ms7CSQ3PQVqAXnXlm7pfI7NZNNZVhWI3fqOpq1Z/9THXcuNzK0Qh69Uo90aVXrwiESVefv/ud6vXXu+2qqtR5CoSiqM933nHeHFq3TowfNUr1ggvc9rhxqj//ec5FTCDew0TQUGAzssimHUgBMR3oISI7iUhz4BTcGu0lQ3OgFXAesKhZM0a1bk3HbA2wPfYYDBwIffvCBRe4L5/6qPniq6x0Te0CHtOZOxd69UqM69XLxeeMTOvzkUdgxAi3XVYG222XQ+EaR1HU56BB0KlT3fjnnoMzz3TbQ4fClCnuNZ0vbrkF6hs4T0WRzsgqKgWiqpU431yTgXnARFXN5SOdNxIUB24AruOWLdnrH503DyZMgLffhlmzXFfU4483fNzgwfCDH8A229Tfr1EAzJ2b+FmX05ddpvW5dq37veEG6NcPTjoJvvgihwI2noKuz/r4/HPYcUe3XVEBbdvC6tXZk7U+Vqxw/9nNmzM7rkjHQopKgQCo6ouqupuq7qKqt0UtT2NJqTjiM2Try2TKFJg5EwYMcF94U6a46YMNMXmy+1Ns2gSvvdZ4OUqFTOuzshKWLYMf/xj+8x/Ybz+I2VsYoZ/PVKRqbeSr9Rym9VFDEbZCim0QvWRoDpTjZmrcQJLSiKfmy+SGGxo3V1zVNetvTznjuX5atoTjjnNdAz/5SXgZSolM67NDB9hqKzjhBLd/0knw8MO5k6/YaMzzmUyXLrB0qfutrIRvvoFtt238eRsibOujhmz91/NI0bVAip0GWxypyMaXyWGHwVNPwZdfuv2vv4YlS9Ln37DB/SHA/QlffBH22KNxMpQSmdanCBx7LLzxhtufMqXuIENTJtP6rI/jjoMa/2dPPQWHHpqfFkhjWh81FFsrJMhIe7GGQpqFFZtVBboik9kZ2ZylMX686t57q+61l2q/fm4Wi2rqWS4rV6r27+/y9uqlesklqikWGWrSZFKfqqqLF6secIDLf+ihqrbeSyKZ1ufVV6t27qwq4n5vvNHFf/ed6tChqrvsojpggOqnn+Ze9jAzrwp4RhYBZ2EVlSV6phSCJXrgrqqGKFCLVcMwgIsucl2SYbuv4imA/3qpWqIXFP/6l5tQU1HhWsrJ9G3ZMrOuqvrYvBn+/e/GnKHg+ctfYK+93Bjq/vvDRyk9nBlBGTMGtt/e1WffvvDQQ1FLVLz85je19bjbbtCuXVKGd97JjvKAovqv2yB6I+ja1f1J//jH1Onvf/ddXuUpdk47DX71K7c9aRJccQW89FK0MhU7J58M994btRTFz1131W7fcw+8/35ShjoRTQNrgaRg8WLo2RPOPx9694YjjoBUuqB7d+jTx9mEGekJWp/xnio2bixou8VICVqfRsOEqctx4+DUU/MiXsFjr740LFgAF1/sjKfatYOnn45aouImaH2OGgW77ALXXAN3351fGYuJoPX59NPuI2foUDez1ahLJv/1JUvgs8/cxC7DFEhadtrJ9XcC7LOP+1IxwhO0Pi++GD79FO68E269NW/iFR1B6vPYY1387Nlw+OG13j2MRDL5r48f75RxgfoUzTumQNLQokXtdnl5MKeaRnoyrc9TToFnn82tTMVMkPrs0KE23/nnO0Nvoy6ZPJvjx1v3VTymQIyCYcGC2u1//AN69IhOllIgfkXkSZNcX78RnvnzYc0a54XGcNgsrEYwfbrzTLFmDTz/PNx4Y44dzpU4994Lr74KzZpB+/a1xsRGOO6+2ymOigrnyWPMmKglKm7GjXMtY5vcUYsZEhqGYRgJmCGhYRiGkVNMgRiGYRihMAViGIZhhMIUiGEYhhEKUyCGYRhGKEyBGIZhGKEwBWIYhmGEIhIFIiInichcEakWkf5JaSNEZKGIzBeRwXHxR/q4hSIyPP9SG4ZhGPFE1QKZA/wM+Fd8pIj0Ak4BegNHAveJSLmIlOPWZDoK6AWc6vMahmEYERGJKxNVnQduydckjgfGq+om4DMRWQgM9GkLVXWRP268z2tr1hmGYUREoY2BdAbiVy1Y5uPSxddBRIaJyAwRmbFq1aqcCWoYhtHUyVkLREReJfUy4Ner6nPpDksRp6RWdCmdeKnqA8AD4HxhBRDVMAzDCEHOFIiqHh7isGXAjnH7XYDlfjtdvGEYhhEBhdaFNQk4RURaiMhOQA9gGjAd6CEiO4lIc9xA+6QI5TQMw2jyRDKILiInAPcA2wP/EJFZqjpYVeeKyETc4HglcLGqVvljLgEmA+XAI6pqK28YhmFEiK0HYhiGYSRg64EYhmEYOcUUiGEYhhEKUyCGYRhGKEyBGIZhGKEwBWIYhmGEwhSIYRiGEQpTIIZhGEYoTIEYhmEYoTAFYhiGYYTCFIhhGIYRClMghmEYRihMgRiGYRihMAViGIZhhMIUiGEYhhEKUyCGYRhGKEyBGIZhGKEwBWIYhmGEwhSIYRiGEQpTIIZhGEYoTIEYhmEYoYhEgYjIH0TkYxGZLSJ/F5F2cWkjRGShiMwXkcFx8Uf6uIUiMjwKuQ3DMIxaomqBvALsqap9gE+AEQAi0gs4BegNHAncJyLlIlIOjAKOAnoBp/q8hmEYRkREokBU9WVVrfS77wJd/PbxwHhV3aSqnwELgYE+LFTVRaq6GRjv8xqGYRgRUQhjIOcA//TbnYGlcWnLfFy6+DqIyDARmSEiM1atWpUDcQ3DMAyAilydWEReBTqmSLpeVZ/zea4HKoHHaw5LkV9Jreg0Vbmq+gDwAED//v1T5jEMwzAaT84UiKoeXl+6iJwJHAMcpqo1L/plwI5x2boAy/12unjDMAwjAqKahXUkcC1wnKp+G5c0CThFRFqIyE5AD2AaMB3oISI7iUhz3ED7pHzLbRiGYdSSsxZIA9wLtABeERGAd1X1V6o6V0QmAh/hurYuVtUqABG5BJgMlAOPqOrcaEQ3DMMwAKS296j06N+/v86YMSNqMQzDMIoKEZmpqv0bylcIs7AMwzCMIsQUiGEYhhEKUyCGYRhGKEyBGIZhGKEwBWIYhmGEwhSIYRiGEQpTIIZhGEYoTIEYhmEYoTAFYhiGYYTCFIhhGIYRClMghmEYRihMgRiGYRihMAViGIZhhMIUiGEYhhEKUyCGYRhGKEyBGIZhGKEwBWIYhmGEwhSIYRiGEQpTIIZhGEYoTIEYhmEYoYhEgYjILSIyW0RmicjLIrKDjxcRuVtEFvr0fnHHnCkiC3w4Mwq5DcMwjFqiaoH8QVX7qGpf4AXgdz7+KKCHD8OA+wFEZFvgRmBfYCBwo4i0z7vUhmEYRoxIFIiqrovbbQ2o3z4e+Js63gXaiUgnYDDwiqp+raprgFeAI/MqtGEYhpFARVQFi8htwC+Bb4BDfHRnYGlctmU+Ll18qvMOw7VeADaIyPwsip0rtgO+ilqIHGLXV9zY9RUvYa+tW5BMOVMgIvIq0DFF0vWq+pyqXg9cLyIjgEtwXVSSIr/WE183UvUB4IFwUkeDiMxQ1f5Ry5Er7PqKG7u+4iXX15YzBaKqhwfM+gTwD5wCWQbsGJfWBVju4w9Oin+j0UIahmEYoYlqFlaPuN3jgI/99iTgl3421iDgG1VdAUwGjhCR9n7w/AgfZxiGYUREVGMgd4jI7kA1sAT4lY9/ERgCLAS+Bc4GUNWvReQWYLrPd7Oqfp1fkXNKUXW5hcCur7ix6ytecnptoppyKMEwDMMw6sUs0Q3DMIxQmAIxDMMwQmEKJGJE5EgRme/dtwyPWp5MEZEdReR1EZknInNF5Nc+flsRecW7nnmlxnNAfe5qChkRKReR90XkBb+/k4i8569vgog09/Et/P5Cn949SrmDICLtROQpEfnY38f9Sun+ichv/LM5R0TGiUjLYr5/IvKIiHwpInPi4jK+X9lwD2UKJEJEpBwYhXPh0gs4VUR6RStVxlQCV6pqT2AQcLG/huHAFFXtAUzx+5DGXU0R8GtgXtz+ncBd/vrWAOf6+HOBNaq6K3CXz1fo/B/wkqruAeyNu86SuH8i0hm4DOivqnsC5cApFPf9G0NdTxwZ3a+suYdSVQsRBWA/YHLc/ghgRNRyNfKangN+AswHOvm4TsB8v/1X4NS4/LF8hRpwdkdTgENxvtsEZ91bkXwfcdPL9/PbFT6fRH0N9VxbG+CzZBlL5f5R68ViW38/XsC5Rirq+wd0B+aEvV/AqcBf4+IT8gUN1gKJlsAuWooB39z/EfAe8EN1Njz43x/4bMV4zX8GrsFNOwfoAKxV1Uq/H38Nsevz6d/4/IXKzsAqYLTvontIRFpTIvdPVT8H/gj8F1iBux8zKZ37V0Om9ysr99EUSLQEdtFS6IjI1sDTwOWa6CyzTtYUcQV7zSJyDPClqs6Mj06RVQOkFSIVQD/gflX9EbCR2u6PVBTV9flumeOBnYAdcM5bj0qRtVjvX0M02j1UfZgCiZZ0rluKChFphlMej6vqMz76C+9JGf/7pY8vtmv+MXCciCwGxuO6sf6M8xRdY4gbfw2x6/PpbYFCNnpdBixT1ff8/lM4hVIq9+9w4DNVXaWqW4BngP+hdO5fDZner6zcR1Mg0TId6OFnhDTHDe5NilimjBARAR4G5qnqn+KSJgE1MzvOxI2N1MSncldTkKjqCFXtoqrdcffnNVU9HXgdGOqzJV9fzXUP9fkL9gtWVVcCS71nCIDDgI8okfuH67oaJCJb+We15vpK4v7Fken9yo57qKgHg5p6wLlu+QT4FOepOHKZMpR/f1zTdzYwy4chuH7jKcAC/7utzy+4mWefAh/iZsdEfh0Br/Vg4AW/vTMwDed250mghY9v6fcX+vbYhLMAAAJESURBVPSdo5Y7wHX1BWb4e/gs0L6U7h/we5y/vTnAo0CLYr5/wDjceM4WXEvi3DD3CzjHX+dC4OwwspgrE8MwDCMU1oVlGIZhhMIUiGEYhhEKUyCGYRhGKEyBGIZhGKEwBWIYhmGEwhSIYeQIEWklIm96p5nZPveroZzfGUYWMQViGLnjHOAZVa3KwbkfBS7KwXkNIzCmQAwjQ0RkgF9boaWItPZrTeyZIuvpeItgETnYt0YmisgnInKHiJwuItNE5EMR2cXnGyMi94tbY2WRiBzk13+YJyJj4s49CedR1TAio6LhLIZhxKOq00VkEnAr0Ap4TFXnxOfxrml2VtXFcdF7Az1xvpUWAQ+p6kBxi3BdClzu87XH+dw6Dnge54/rPGC6iPRV1VmqusYvftRBVVfn7GINox6sBWIY4bgZt+5Jf2BkivTtgLVJcdNVdYWqbsK5lnjZx3+IW9+hhufVuYj4EPhCVT9U1WpgblK+L3EeZg0jEkyBGEY4tgW2BrbB+U9K5rsU8Zvitqvj9qtJ7A3YlCJPqnwtfTmGEQmmQAwjHA8ANwCPk2LZU1VdA5SLSCrl0mi8Z9mOwOJcnN8wgmAKxDAyRER+CVSq6hPAHcAAETk0RdaXcd6Kc8E+wLtau6qeYeQd88ZrGDlCRH4EXKGqZ+Tg3P8HTFLVKdk+t2EExVoghpEjVPV94PVcGBICc0x5GFFjLRDDMAwjFNYCMQzDMEJhCsQwDMMIhSkQwzAMIxSmQAzDMIxQmAIxDMMwQvH/AWsrHIlD+0o5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(r[ix],r[iy],'-',color='k')\n",
    "plt.plot(r[ix],r[iy],'o',color='b')\n",
    "plt.plot(r[0],r[1],'^',color='r',markersize=20)\n",
    "plt.plot(r[0],r[1],'>',color='k',markersize=20)\n",
    "plt.plot(r[-2],r[-1],'^',color='r',markersize=20)\n",
    "# label the nodes\n",
    "for n in nodes:\n",
    "    if n[2]>0.8*l: offset=0.1\n",
    "    else: offset=-l/5\n",
    "    plt.text(n[1]-l/3,n[2]+offset,'n {}'.format(int(n[0])),color='b')\n",
    "# label the elements\n",
    "for e in elems:\n",
    "    n1=nodes[e[1]-1]\n",
    "    n2=nodes[e[2]-1]\n",
    "    x=np.mean([n2[1],n1[1]])\n",
    "    y=np.mean([n2[2],n1[2]])\n",
    "    # ----------------->need elem labels<-----------------\n",
    "    plt.text(x-l/5,y-l/10,'el {}'.format(int(e[0])),color='r')\n",
    "plt.title('Our truss structure\\neach element is 2 nodes and length L\\ntriangle markers are constraints')\n",
    "plt.xlabel('x (mm)')\n",
    "plt.ylabel('y (mm)')\n",
    "plt.axis(l*np.array([-0.5,3.5,-1,1.5]));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assemble the stiffness matrix\n",
    "\n",
    "Creating the matrix to represent $\\mathbf{K} = \\mathbf{K}_1+...\\mathbf{K}_{11}$ is called __assembling the stiffness matrix__. We loop through each element and add the individual stiffness matrices, the same way we did for the 3-element truss above. First, we define our $\\mathbf{K_{el}}$ function, `Kel(node1,node2)`. The nodes define the initial length and the angle of the beam element. We can multipy the stiffness matrix by $EA$ to get the actual N/m stiffness term. In our case, each beam is the same material and cross section, so we can factor out $EA$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Kel(node1,node2):\n",
    "    '''Kel(node1,node2) returns the diagonal and off-diagonal element stiffness matrices based upon\n",
    "    initial angle of a beam element and its length the full element stiffness is\n",
    "    K_el = np.block([[Ke1,Ke2],[Ke2,Ke1]])\n",
    "    \n",
    "    Out: [Ke1 Ke2]\n",
    "         [Ke2 Ke1]   \n",
    "    arguments:\n",
    "    ----------\n",
    "    node1: is the 1st node number and coordinates from the nodes array\n",
    "    node2: is the 2nd node number and coordinates from the nodes array\n",
    "    outputs:\n",
    "    ----------\n",
    "    Ke1 : the diagonal matrix of the element stiffness\n",
    "    Ke2 : the off-diagonal matrix of the element stiffness\n",
    "    '''\n",
    "    a = np.arctan2(node2[2]-node1[2],node2[1]-node1[1])\n",
    "    l = np.sqrt((node2[2]-node1[2])**2+(node2[1]-node1[1])**2)\n",
    "    Ke1 = 1/l*np.array([[np.cos(a)**2,np.cos(a)*np.sin(a)],[np.cos(a)*np.sin(a),np.sin(a)**2]])\n",
    "    Ke2 = 1/l*np.array([[-np.cos(a)**2,-np.cos(a)*np.sin(a)],[-np.cos(a)*np.sin(a),-np.sin(a)**2]])\n",
    "    return Ke1,Ke2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we loop through each element (1-11) and add each element stiffness term to the _global_ stiffness matrix, $\\mathbf{K}$. The result is a $14\\times 14$ array that satisfies our $\\sum \\mathbf{F} = \\mathbf{0}$ at all 7 nodes. \n",
    "\n",
    "$\\mathbf{F} = \\mathbf{Ku}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pretty_print as pp\n",
    "\n",
    "K=np.zeros((len(nodes)*2,len(nodes)*2))\n",
    "for e in elems:\n",
    "    ni = nodes[e[1]-1]\n",
    "    nj = nodes[e[2]-1]\n",
    "    \n",
    "    Ke1,Ke2 = Kel(ni,nj)\n",
    "    #--> assemble K <--\n",
    "    i1=int(ni[0])*2-2\n",
    "    i2=int(ni[0])*2\n",
    "    j1=int(nj[0])*2-2\n",
    "    j2=int(nj[0])*2\n",
    "    \n",
    "    K[i1:i2,i1:i2]+=Ke1\n",
    "    K[j1:j2,j1:j2]+=Ke1\n",
    "    K[i1:i2,j1:j2]+=Ke2\n",
    "    K[j1:j2,i1:i2]+=Ke2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \\left[ \\begin{array}{cccccccccccccc}\n",
      " 4.2 & 1.4 & -0.8 & -1.4 & -3.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " 1.4 & 2.5 & -1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " -0.8 & -1.4 & 5.0 & 0.0 & -0.8 & 1.4 & -3.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " -1.4 & -2.5 & 0.0 & 5.0 & 1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " -3.3 & 0.0 & -0.8 & 1.4 & 8.3 & 0.0 & -0.8 & -1.4 & -3.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " 0.0 & 0.0 & 1.4 & -2.5 & 0.0 & 5.0 & -1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & -1.4 & 8.3 & 0.0 & -0.8 & 1.4 & -3.3 & 0.0 & 0.0 & 0.0 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & -1.4 & -2.5 & 0.0 & 5.0 & 1.4 & -2.5 & 0.0 & 0.0 & 0.0 & 0.0 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & 1.4 & 8.3 & 0.0 & -0.8 & -1.4 & -3.3 & 0.0 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.4 & -2.5 & 0.0 & 5.0 & -1.4 & -2.5 & 0.0 & 0.0 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & -1.4 & 5.0 & 0.0 & -0.8 & 1.4 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -1.4 & -2.5 & 0.0 & 5.0 & 1.4 & -2.5 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -3.3 & 0.0 & -0.8 & 1.4 & 4.2 & -1.4 \\\\\n",
      " 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.4 & -2.5 & -1.4 & 2.5 \\\\\n",
      "\\end{array}\\right]\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pp.array_latex([K*1000],['[K]'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the inputs for the project\n",
    "\n",
    "Here, we save the arrays needed to work on the project as a Linear Algebra problem. We need the `nodes`, `elems`, and `K` arrays. So, we save them in the file `fea_arrays.npz`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.savez('fea_arrays',nodes=nodes,elems=elems,K=K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solution for our FEA problem\n",
    "\n",
    "In FEA, you have applied forces, in $\\mathbf{F}$, and displacement in $\\mathbf{u}.$ Typically, you have a mix of known and unknown forces and displacements. In our example, we have 14 degrees of freedom - 3 constraints for our displacements. So, we solve the problem in two steps\n",
    "\n",
    "1. solve for the free nodes `uf` $\\rightarrow \\mathbf{K}[2:13,2:13] \\mathbf{u}[2:13] = \\mathbf{F} [2:13]$\n",
    "\n",
    "2. plug in `u` to solve for `F` at constrained degrees of freedom $\\rightarrow \\mathbf{F}=\\mathbf{K}\\mathbf{u}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   0.    0.    0.    0.    0. -100.    0.    0.    0.    0.    0.]\n",
      "[ 1.42108547e-14  0.00000000e+00 -2.84217094e-14  0.00000000e+00\n",
      "  7.10542736e-15 -1.00000000e+02  2.84217094e-14  1.13686838e-13\n",
      " -7.10542736e-15  0.00000000e+00  0.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "E=200e3\n",
    "A=0.1\n",
    "Ff=np.zeros(2*len(nodes)-3)\n",
    "Ff[5]=-100\n",
    "print(Ff)\n",
    "# step 1 solve for uf (the joints without constraints)\n",
    "uf = np.linalg.solve(E*A*K[2:13,2:13],Ff)\n",
    "u=np.zeros(2*len(nodes))\n",
    "u[2:13]=uf\n",
    "\n",
    "# step 2 solve for F (the solution should include reactions and applied forces)\n",
    "F=E*A*K@u\n",
    "print(F[2:13])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "displacements:\n",
      "----------------\n",
      "u_1x:0.00 mm\n",
      "u_1y:0.00 mm\n",
      "u_2x:1.95 mm\n",
      "u_2y:-2.12 mm\n",
      "u_3x:0.43 mm\n",
      "u_3y:-4.00 mm\n",
      "u_4x:1.08 mm\n",
      "u_4y:-5.37 mm\n",
      "u_5x:1.73 mm\n",
      "u_5y:-4.00 mm\n",
      "u_6x:0.22 mm\n",
      "u_6y:-2.12 mm\n",
      "u_7x:2.17 mm\n",
      "u_7y:0.00 mm\n",
      "\n",
      "forces:\n",
      "----------------\n",
      "F_1x:-0.00 N\n",
      "F_1y:50.00 N\n",
      "F_2x:0.00 N\n",
      "F_2y:0.00 N\n",
      "F_3x:-0.00 N\n",
      "F_3y:0.00 N\n",
      "F_4x:0.00 N\n",
      "F_4y:-100.00 N\n",
      "F_5x:0.00 N\n",
      "F_5y:0.00 N\n",
      "F_6x:-0.00 N\n",
      "F_6y:0.00 N\n",
      "F_7x:0.00 N\n",
      "F_7y:50.00 N\n"
     ]
    }
   ],
   "source": [
    "xy={0:'x',1:'y'}\n",
    "print('displacements:\\n----------------')\n",
    "for i in range(len(u)):\n",
    "    print('u_{}{}:{:.2f} mm'.format(int(i/2)+1,xy[i%2],u[i]))\n",
    "print('\\nforces:\\n----------------')\n",
    "for i in range(len(F)):\n",
    "    print('F_{}{}:{:.2f} N'.format(int(i/2)+1,xy[i%2],F[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Postprocess FEA results\n",
    "\n",
    "The last step in FEA (and any applied linear algebra problem), is to post-process the results. I have chosen to show the deflections of each joint and the applied forces at each node. I imported the `interact` widget to animate the deformation of the structure. What you should see is that the sum of external forces $\\sum F_x=0$ and $\\sum F_y=0$. This is what we see, since our the sum of the reaction forces is equal to the applied force. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6f9b885a160e4e9a9bc89ccd96605ec4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=5, description='s', max=10), Output()), _dom_classes=('widget-interact',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "\n",
    "def f(s):\n",
    "    plt.plot(r[ix],r[iy],'-',color=(0,0,0,1))\n",
    "    plt.plot(r[ix]+u[ix]*s,r[iy]+u[iy]*s,'-',color=(1,0,0,1))\n",
    "    #plt.quiver(r[ix],r[iy],u[ix],u[iy],color=(0,0,1,1),label='displacements')\n",
    "    plt.quiver(r[ix],r[iy],F[ix],F[iy],color=(1,0,0,1),label='applied forces')\n",
    "    plt.quiver(r[ix],r[iy],u[ix],u[iy],color=(0,0,1,1),label='displacements')\n",
    "    plt.axis(l*np.array([-0.5,3.5,-0.5,2]))\n",
    "    plt.xlabel('x (mm)')\n",
    "    plt.ylabel('y (mm)')\n",
    "    plt.title('Deformation scale = {:.1f}x'.format(s))\n",
    "    plt.legend(bbox_to_anchor=(1,0.5))\n",
    "interact(f,s=(0,10,1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "steel=np.loadtxt('../data/steel_price.csv',skiprows=1,delimiter=',')\n",
    "al = np.loadtxt('../data/al_price.csv',skiprows=1,delimiter=',')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc8fd7b9c50>]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3XmcXFWZ8PHfU0uvWTo7IQthCTsKIbIJyGZY5DUwgMKAREVxFHx1XHGWFzfGZUQdR8VBQVERREBAUSBGURgMJGGHEBICJJ21k046Sa+1PO8f51bqVtWtrbu6u7rr+X4+/al7z13q3HTnPvcs9xxRVYwxxtSe0HBnwBhjzPCwAGCMMTXKAoAxxtQoCwDGGFOjLAAYY0yNsgBgjDE1ygKAMcbUKAsAxhhToywAGGNMjYoMdwYKmTx5ss6ZM2e4s2GMMSPKihUrtqnqlGL7VXUAmDNnDsuXLx/ubBhjzIgiIm+Wsp9VARljTI2yAGCMMTXKAoAxxtQoCwDGGFOjLAAYY0yNsgBgjDE1ygKAMbUg1gPta4c7F6bKFA0AInKIiDzr+9klIp8UkYkislhEVnufE7z9RUS+JyJrROR5EZnnO9cib//VIrJoMC/MGOOJ98H33wbfOwaW/mi4c2OqSNEAoKqrVPVoVT0aOBboAn4LXAcsUdW5wBJvHeBcYK73czVwE4CITASuB44HjgOuTwUNY8wgeu4O6Fjnlh/6/PDmxVSVcquAzgReU9U3gYXAbV76bcAF3vJC4OfqLAVaRGQ6cDawWFXbVXUHsBg4Z8BXYIwprKdjuHNgqlS5AeBS4A5veZqqbgLwPqd66TOA9b5jWr20fOnGmMEUquoRX8wwKjkAiEgd8G7gN8V2DUjTAunZ33O1iCwXkeVtbW2lZs8Yk08oPNw5MFWqnBLAucDTqrrFW9/iVe3gfW710luBWb7jZgIbC6RnUNWbVXW+qs6fMqXoYHbGmGIsAJg8ygkAl5Gu/gF4AEj15FkE3O9Lv9LrDXQC0OFVET0MLBCRCV7j7wIvzRgzmMQCgAlWUuWgiDQB7wQ+4kv+OnCXiFwFrAMu8dL/AJwHrMH1GPoAgKq2i8hXgGXefl9W1fYBX4ExpjBrAzB5lPSXoapdwKSstO24XkHZ+ypwTZ7z3ArcWn42jTH9ZgHA5GFvAhsz2lkAMHlYADBmtLNGYJOHBQBjak0yMdw5MFXCAoAxo50mM9cTseHJh6k6FgCMGe00633LRN/w5MNUHQsAxox2mlXlYwHAeCwAGDPa5VQBWQAwjgUAY0a77EZfCwDGYwHAmNHOGoFNHhYAjBntrA3A5GEBwJjRztoATB4WAIwZ7XK6gVoVkHEsABgz2lkjsMnDAoAxo51VAZk8LAAYM9rlNAJbFZBxLAAYM9pZCcDkYQHAmNEuuw0g3js8+TBVxwKAMaOdvQhm8rAAYMxoZ6OBmjwsABgz2mU3Aj93B9x8Ojzzy+HJj6kaJQUAEWkRkbtF5BURWSkiJ4rIRBFZLCKrvc8J3r4iIt8TkTUi8ryIzPOdZ5G3/2oRWTRYF2WM8cmuAlr3d9j4NNx/jVUH1bhSSwD/BTykqocCbwVWAtcBS1R1LrDEWwc4F5jr/VwN3AQgIhOB64HjgeOA61NBwxgziApNAdm1fejyYapO0QAgIuOAU4FbAFS1T1V3AguB27zdbgMu8JYXAj9XZynQIiLTgbOBxararqo7gMXAORW9GmNMruwSgF/ntqHLh6k6pZQADgDagJ+KyDMi8hMRaQamqeomAO9zqrf/DGC97/hWLy1fegYRuVpElovI8ra2trIvyBiTJbsNwK9rG/z1m3DjobDsJ0OXJ1MVSgkAEWAecJOqHgN0kq7uCSIBaVogPTNB9WZVna+q86dMmVJC9owxBfV15d+2/TX4yw2wexM8+OncHkNmVCslALQCrar6pLd+Ny4gbPGqdvA+t/r2n+U7fiawsUC6MWYwbX05/7ZNz2au9+4a3LyYqlI0AKjqZmC9iBziJZ0JvAw8AKR68iwC7veWHwCu9HoDnQB0eFVEDwMLRGSC1/i7wEszxgwWVdj8Yv7tG7MCwK5Ng5sfU1UiJe73ceB2EakD1gIfwAWPu0TkKmAdcIm37x+A84A1QJe3L6raLiJfAZZ5+31ZVdsrchXGmGA710FvR/7tm5/PXN+9EaYeOrh5MlWjpACgqs8C8wM2nRmwrwLX5DnPrcCt5WTQGDMAWwo8/QfZvXlw8mGqkr0JbMxotn1Nefvvsma5WmIBwJjRrNxxf6wEUFMsABgzmiULvAQWZLc1AtcSCwDGjGaF3gIOYiWAmmIBwJjRrNBbwEG6rWNeLbEAYMxoVm4JwAaHqykWAIypdm88Dnd/EFYvLv/YcgNAT4cNEV1DLAAYU+1+9i548R64/WKIl9mrJ99Q0ONnBacDdO8o7zvMiGUBwJhqFuvOXO8p8FZvkHwlgFnHw6V3QLQJJJy5zaqBaoYFAGOqWfZ4/TvXwS0L4EenwM71wcf45QsAzZPh0PPgE8/B516D2Semt1kAqBkWAIypZp1Zc2L88XOw/kk3hs9jNxY/Pl8AaJrsPsdMhcYJ0DQpvc0CQM2wAGBMNcsuAWxYnl5e8dPix+drA2ielLneNDG9bAGgZlgAMKaadQ1wysa8JYDsAOBb77QAUCssABhTzbKrgLJ1FXlxK9+LYKkqoL3rVgVUiywAGFPNik3avnVl4e2FGoH9LADUJAsAxlSzYgHg9b/CnZfDvR/J7TIK+QeDa5yYuW4BoCaVOiOYMWY4FKsC+us30suT58Kpn8ncHlQCaJ5SuA3AAkDNsABgTDULagSWUPCNfekPAwKArw3glE9DKAqHvxtCWYX/jF5A/RwQLplweRPp3/FmyFkVkDHVLKgK6Ix/C943Ec9N8weKKYfC6V+AaUfk7udvFO5sc5PJl2Pdk3DjofDj04OrokxVsgBgTLVSDa4Cmv7W4P2DZv/yBwAp8N+9fixEGt1yvBt6d5eeT4BfXAidW2HjM/D4d8o71gybkgKAiLwhIi+IyLMistxLmygii0Vktfc5wUsXEfmeiKwRkedFZJ7vPIu8/VeLyKLBuSRjRom+Toj35KZHm3O7cQIkA0bx9L8IVigAiMDYaen1PVtKzydArDO9vPHZ8o41w6acEsDpqnq0qs731q8DlqjqXGCJtw5wLjDX+7kauAlcwACuB44HjgOuTwUNY2qSqhvO4ff/DHsCnvTzNQBLCMbtm5uejOdWA5VaAgAYM4AAkPGdZU5CY4bNQKqAFgK3ecu3ARf40n+uzlKgRUSmA2cDi1W1XVV3AIuBcwbw/caMbGuWwJIvw/Jb4W/fzN2erwtoKAzjZwZvy75x9zcADGRqyHzDT5iqU2oAUOAREVkhIld7adNUdROA9znVS58B+IcpbPXS8qVnEJGrRWS5iCxvayvSBc6YkWL7a+5G73/S/9/vppefujn3mHzDQIjAuJz/Ok5Ha+a6PwCEsoZ9zlZqCSCZgJfvdxPUBDUWlzsJjRk2pXYDfbuqbhSRqcBiEXmlwL5BfcC0QHpmgurNwM0A8+fPL7MrgjFVKBGDn57rbqqv/AGuuNulF7sh56sCUg2uAgI3XPTs49PrpbYBQOltAC/eA/d+2C0v+h3sf2pW/iwAjBQllQBUdaP3uRX4La4Of4tXtYP3udXbvRXwTzc0E9hYIN2Y0W3nuvQNdc3i9JSLW17O3K93T+Z6viqg3t1uGOcgO97IXM+oAipWAtgnvby7QABI3fwB7r82d7sFgBGjaAAQkWYRGZtaBhYALwIPAKmePIuA+73lB4Arvd5AJwAdXhXRw8ACEZngNf4u8NKMGd1CWQXtLS/CXYtct0m/9rWZ6/kCQF8n1DUHb9vxeua6llECyGgD2FR435R4b26atQGMGKVUAU0Dfivu7b4I8CtVfUhElgF3ichVwDrgEm//PwDnAWuALuADAKraLiJfAZZ5+31ZVfv5yqExI0h2r5ibTwver30tTH9Lej1fFdCMY10QCVKoBJD99m+28b52hV0DKJxbCWDEKBoAVHUtkPPmiapuB84MSFfgmjznuhW4tfxsGjOClfpWbftrmev+RuAFX4W2VXDAaTBueu6T/t5zZKWX0wbg71nU0ery3Z9hHawb6IhhYwEZM1jivfD6YxBtKG3/7dlVQL4SwH5vh5M+nl7PVwW0e6MbiiHqvdXrDz7FAkDDeKgfB7273NvAXe25M4eVwkoAI4YFAGMGy8P/Cst+XPr+hdoAssfvrxuT/zxbXoaZx7rlchqBwZUCtnqN0x3riweAoBKCBYARw8YCMmYwqBa/+Z95PZz6ufR6+2vw2p/hvmtg/bLMAJA99EN2CWD/d6SXNz7ty0cZVUCQWw3UH/nmIDBVxwKAMYMhXx19SrgOTvgYvOPzbohmcF1Ff3EhPPtL+OVF6bF96sZAXVPm8dkBYO4708sb/AGgjBfBIDMAPPr18kcFzf5OU9UsABgzGPw34SD7znNtA+EITNgvd3tvR3o5e/IWcKN3HniGWz78ApgxP73NXwIopxEYMgPAlhfg+V9nbn/s25nrQV0+rRF4xLA2AGMGw8ZnCm+f9bb08sQDYfua/Ps2TwlO/8ffuPr6aUe4UUNTE8W0rXIvi9WPzeynH6kvnu/xszLXf/sReOulblkVlnwpc3vnVjckhF8yYF4CU5WsBGDMYFj/VOHt045KL088oPC++QJAOOLeGwiFXZXQlMO8DQqbnnOLfb63i/P1HPILGmSuc5sLJJtfCD7m9osz14NeDjNVyQKAMZXWuxs2rMi/vX4cHLwgvT7pwMLnK7Ur5ox56eV1S91nn2+c/kI9h1KCAsCL98APT4D/OaW0fATNYWCqklUBGVNpbz4RXA9+3rdg1waYuwAafVNhTDm08PnylQCyzT4RnvmFW173d/eZEQBKKAGMnZ6b9sfP5aYVYiWAEcMCgDGV9vrfctOmHQnHfTg3HWC/k1w1UPZ7AClBs38FnufE9PL6p9ygc3Hf/LypKR8LCUdL+65CbE7gEcOqgIyptLV/TS/Xj4dDzoMr7sm/fygMJ/3f/NtLLQFM2D89omfvLviar0E32lx8LKBKScZsQLgRwgKAMZXUud11nwT35u2nXoLL7oCx+xQ+7q2XwfjZwduy3wLORySzFOB/+i+l+mfveUp4X6AYawcYESwAGFNJb/iqf2Yc67piliLaAIsecMdkKzUAAMw+KTi9nABw1SOlVzvlY+0AI4IFAGMqyV//f8A78u8XZOL+cMa/5aaXWgUEmSUAv1J6AKXMnA+feTU4GJXK2gFGBAsAxlSSPwBkT5VYiqC3fst5Gp96eHB6OSUAcO0S4bryjvGzKqARwQKAMZXSsSH9Rm+kAWYeV/45sm/29eMhUsaNON94P+UGABhYjyALACOCBQBjKsX/9D/r+NLnAfDLLgH0Zzz+oHaA7MHkShGyADDaWQAwplK2rUovzz6hf+fIDhr9qYY59v25aeW0Aez97gIBoNjLa9YIPCLYi2DGVMruzenloCEVhkrQoG+VrgJa+EM36mjvLljy5dzt1gg8IpRcAhCRsIg8IyK/99b3F5EnRWS1iPxaROq89HpvfY23fY7vHF/w0leJyNmVvhhjhtXuTenloCEV+qUfc/JWKgAUqgJqbHFvNp/y6eDtVgIYEcqpAvoEsNK3/g3gO6o6F9gBXOWlXwXsUNWDgO94+yEihwOXAkcA5wA/FKnEGyfGVAl/CaDYi1+l6k9detCTe6WrgPxjGQXNMxC3EsBIUFIAEJGZwLuAn3jrApwB3O3tchtwgbe80FvH236mt/9C4E5V7VXV14E1QD+6SRhTpQajBNDdXv4x4SGoAqofl14O6rq65eXyv88MuVJLAN8FPgek5nqbBOxU1dTMD63ADG95BrAewNve4e2/Nz3gGGNGtr4u6PFm8QpFoXFiZc7b01F8n2xBDcf9CQDZ4/kceKY798n/7OYiSJl1fO6xT/0YuneW/51mSBUNACJyPrBVVf0DnAdVTGqRbYWO8X/f1SKyXESWt7W1FcueMdVhT1b1z0AGXjvtX3zLXyj/+EpVAaWGlE654h74wgY464uZ6e+6ERpaIOQLCr0d8OT/lP+dZkiV8lf6duDdIvIGcCeu6ue7QIuIpH7jM4GN3nIrMAvA2z4eaPenBxyzl6rerKrzVXX+lCllvAJvzHCqZP3/yZ+EE6+F4z7iJo4vl796JiXaj/cAjr4ivfzWy9xgc0EvpY3dBz69Cj79Klx4czp96Q/6V4IxQ6ZoAFDVL6jqTFWdg2vE/bOqXg78BUjNBbcIuN9bfsBbx9v+Z1VVL/1Sr5fQ/sBcoMi8ecaMEBn1/wMMAJF6OPsGOO+b0BBwMy9m4gG54wf1pwromCtg7tlu0vl33Vh432iDe2ntyIvSU1z2dMBf/gPWL3PzCZuqM5AXwT4PfEpE1uDq+G/x0m8BJnnpnwKuA1DVl4C7gJeBh4BrVIOmTTJmBMooAVSqC2g/hUJw0FmZaf2pAho3HS6/C95zW+kBJByBUz+bXn/yR3DLWfDE98r/fjPoygoAqvqoqp7vLa9V1eNU9SBVvURVe730Hm/9IG/7Wt/xN6jqgap6iKr+sbKXYswwqmQJoBLmvjNzvT8lgP466j0wYU5m2uL/N3Tfb0pmQ0EYUwnVVAIAOOD0zPX+jAXUX+EIzLty6L7P9JsFAGMqYTBeAhuIponuSRzcENHjZxXev9IGMpeAGTIWAIxJiffC649B7x633rkNfnEh3H5J8d4sgzIMxAAt/D588GH40BLXg2coTT86Ny3eC288nhkszbCyAGDM7i3Q1Q73XAW3nQ8/PQeSSfjf78Jrf4bVj8Bfv1nkHFVWAgDXm2j2CUNb/ZPS2JKbdt9H4Wfvgh+dDLs25W43Q84CgKkdK26DR/7dPdmnvPowfPtQuPEQWPk7l7b5BdjxOjzx3+n9/v79/Oft3Q19Xqkh0uheijK5Q1K8eI/77GyDB/MMImeGlAUAUxs2rIDf/V/XHTH1NK8Kf/oiaBISfZn7t6+F7LEKU33ZVWHl72H5rRDryXyaHbvP0Fe3VKtCcyKvehA2v1je+RJxeOZ2eO5OV0IzA2bzAZjasOn59PIGb1ST1/8KW/MMWrbtVTfKpf9VlV0bYfwMNw7+ry93ad07YMb89D7VUv9fDd71bfjukfm3L/53eN9vSz/fS7+F+703o6ONcPjCgeXPWAnA1Ah/Hf2ON9zn0h/l33/to5CMZaZtfNp9PuGrDlry5eqs/68GLbPg0PPzb3/tz7DuydLPd99H08v3Xt3/fJm9LACY2uAfrK1rmysRvPpQ/v1XP5KbtsELANnDH7f6RjSxEkCmfd5SePuGFYW3+/kDss05XBEWAExtyO56+PC/EDAYbWGty9xny+zM9GU/SS9bCSDT4QtzB6fb95j0csd6ShayGutKswBgasPurG6HbzxW/jlal7m+7FqgAdJKAJmmHgof+zsccp5bH7tv5qT1O9eVfq6giW7MgFhINaPbttVw5+WwbVXw9jHTYM+W0s4V73FBoH1t/n2sBJBr/Ey47A7YuR6aJ2f2/vEHgGW3uB5a8z8Ap13nels9dJ0bUC7Iqw/DwTa1+EBYCcCMbvdfm//mD3Dyp4qfY7+T08uvP2YBoL9aZrneOy2+YSlSVUCq8OCnXFvNo19z71a0Lst/8wf41XsGN781wAKAGd3WLy28fVaRaaklDPPel15/o0AAaJ6S2z5gcjVPTU9b2b3D3ey7tmfu09dV+Oafsm115fNXQywAmNrWML7w9glzMkfWXP8U7NqQu9/cs12f9ojVUxcVCmUOTrdzPex4M3OfnW/Cy/dTVOrtbdMvFgBMbQuaPtFv0kEwdlq650r2uwEp7/0F7HNUZfM2mmVXA+3MCgBLb4JkvPh5Xvl9ZfNVY6wR2Ixese7i+xSbcnHyXPd5zBWw8ZngfY55nz35l8tfVfar98CkuZnbX76vtPNsWAEdra6h2ZTNSgBm9Nq1sfg+xW7ckw5yn0deDJGGzG1HvQfe/6Ab8sCUZ3xWW8n2rLr8VFfbxonFz/XKg5XJUw2yAGBGr1ICAEDd2Pzbxkxzn40tuWPP7HsMzDkZInX9y18taylxgpoTP1Z8n5W/c79rm3i+bBYAzOgU7ys9AGSPXT/rBPcZbXLj6acc877M/SYe0P/81bpSekuF62De+10Du9/CH8I/v5Ref+Mx+PZhcPvFbl3V/f5NUUUDgIg0iMhTIvKciLwkIl/y0vcXkSdFZLWI/FpE6rz0em99jbd9ju9cX/DSV4mIvcFhKk8VfnkRfH02PP6d3O3+YQhSTvA9Zc67Ev7Pf8HbPgSX3OamVkyZczJMPNBbEZh2eEWzXlNKmaLyyItgzBQ49bOZ6WOmujeuJev2teZPsG4p/PexLiAUGmeor9MNBvjan8vP+yhSSiNwL3CGqu4RkSjwuIj8EfgU8B1VvVNEfgRcBdzkfe5Q1YNE5FLgG8B7ReRw4FLgCGBf4E8icrCqf7xdYwZozRJ3IwBoW5m57eyvuSfP1FDOqYbH4652DYnd7XDGv7sbzLtuzD23CFz0Y3j0G3DQWdbnfyBKGTLj+I+4z+wSWvMUCIXdoHydbZnb7roy/Wb3vR+Bjy/PPW8iBre/B9583K1fuwImH1Re/keJogFAVRXwpjsi6v0ocAbwj176bcAXcQFgobcMcDfwfRERL/1OVe0FXheRNcBxwN8rcSHGAG4mryAX/xSO/AdXQnjbh1yPnvP+020LR+Cc/yjt/DOOhcvvqkxea1m4yK1n5nHp0lrjhMxtzVO8z6m5AcA/rEd2w3LKki+lb/4Am5+3AFCIiISBFcBBwA+A14CdqprqqNsKzPCWZwDrAVQ1LiIdwCQv3f9apv8YYyojX0PguH3dp0jw070ZeuNmwq7W4G0nfTy9nD3FZioAjN0Htr5EWeJ98NRPMtOy30KuISU1AqtqQlWPBmbintoPC9rN+wyaD08LpGcQkatFZLmILG9raws4xJgCejuC01MBwFSP8QHPf6EoXHQLHP7udFo4Asd/FBA4/p/Sva7mf9DNwZxPdhsBQLzb/fh1tZed9ZKt/pObdnRnGcNeD6GyegGp6k7gUeAEoEVEUiWImUCqy0UrMAvA2z4eaPenBxzj/46bVXW+qs6fMmVKOdkzBvYEPTQIjLFB2qpO0DAcZ/8HHHVxbvq5X4cvtMK530inHXY+fHY1nHht8Pmzq44AkgFNjoNVAtjTBrdf5Doj3P3BwfmOASqlF9AUEWnxlhuBs4CVwF+A1G9qEZAauOMBbx1v+5+9doQHgEu9XkL7A3MB31RKxlRA59bctMlzra9+NQp6CW9MgYe++jEBaWNh9onB+2cHgK0rYdUfcvfLDgCJOKx6aGADzXW1w+8+kV5vfQrW/rXq3lUopQ1gOnCb1w4QAu5S1d+LyMvAnSLyVeAZ4BZv/1uAX3iNvO24nj+o6ksichfwMhAHrrEeQKbiOrflph14xtDnwxQXNMHLhP3LP0++Ibj9AWD7a3Dz6bnVP5AbAB77lhuSOtoE1y4Prqoq5lfvzZwqFODn73bVW0ElnGFSSi+g54GcztOquhbXHpCd3gNckudcNwA3lJ9NY0q0J6AE4B/N01SP7KE1jnkf7Ht0+edJva2dLRRNLz/zy+CbP+QGgEe/5j5jXS4YnB/wPkm8F3r3QPOk3G27Nube/FP+9q2qCgD2JrAZXYKqgOacnJtmhl92FdDC7/fvPPkCwLon3DDTqvDSvfmPL9QIvHVlbtruzfDtw+HGg11VUbbU3NFBSn07fYhYADCjRyLmJhjJFlR3bIZfdgmg3+cp0L7z4KfcOx873si/T6FG4I6AbqrP/gq6trnhqp/+ee72QgGgvsC4U8PAhoM2o0dQ/f+FNw99PkxphqJhfs2fYGpQr3WfeLerOrznKujMCgYd6yGZdJPYpLz+t/Tytldzz7e+QADYvdE9qISj+fcZQlYCMKOH/63QUAQ+tATeYvPGVq1KlQDANdjm82Seh4Boc3p5yZfdjT3oxTL/mELxXjfeUEr7Wvdy2Qt3wy1nwzO3w6Zn8+dFk8GlimFiAcCMHv6i/KwTYOZ89+avqU4T5lTuXJfdmb8tINEbnO6fD/qFu/Ofe6VvasrWZZmNyZpwQeCPn3PzT9//MYj3uG1j87x8mD372TCyAGBGnkSeqQK7fY15TSVMJGKG11GXuLGVwnVw4f8M7FwHvAM+vQo++ULpx0z0dTnN10MIYIuvVLD2r7nbt60KbkeY83b4YgdcvxPe8t50evb8x8PI2gDMyLBnq5skXBX+8lWYchi8//eZdaldFgBGlFDYVdP17i4+NWcpRArPINYyG3auc8snXpv//YFs/p5A/vr/vdtfCT5u5tvS+fKPHFtFJQALAGZkuPfDsPbR9Pr6pa5v9/wPpNMyAkBA/2xTfUQqc/NPqR/jJoy5P2AmsWM/4N7uTcbhtOuCn+aD7N7k/rbCdbAhYHjp7GHHU/xzT7Tsl15OBaEqYAHAjAz+m3/Kpucy1/1VQKXMJWtGp2Mud5P13HxaZnr9WLjwpvR6OTO6tb3iJpFJetWP4TpIeLOObXk5+Jg6XyPzBF8AqKIqIGsDMCNY1rgqVgVkUvbO3OYTCmeul9MIveWlzIeQo96THm1026o8B/k6IGSUACwAGDNwK34GT/04PcBW7+70tvoKViuYkadhHEw+JDMt+4Zf15Q5M5m/WyhAvW+00q0rM+v/D/s/6Tr+IKFIZhvDuBkgXgDaswUe+kL+toMhZAHAjGx/+Ay88Bu3rMl0evbTnqk9F97k5hU+6J1w1pdg/9Ny9/EPPnfMFZnb9jspvbz+Kdjs9TCSsNt2cJ5pzY++Av7h5sxSaDgC42em15f+EH6+sKzLGQwWAMzI99Jv3ac/AIgFgJo341i4+Fa44m44+ZOZb/OmHHKO+4w2wQkfzdzmHwV0ywvsrXIcO92VMOYGBICmyXDBD1zgydmW1TFhz+bgt9eHkAUAM/KtfRRiPe6lnJSg2aCMyXbix2HR7+CfHs98LwDyzzOQCiTTjnDTWvpNf0v+7wqa/yBosLkhZP9G2KuFAAAYWUlEQVRLTPV74r8Lb491wZv/m1UFZH/apgShEOx/KkzyGo0vv9v1IDvgdDjiQvdEny31cCECBy/I3JYvaIDrOZStbXjbAex/ialuf/8BPPJvxfdb/Ujm25jZDXrGlGLuO+Gzr8GV97l2pKDJYPyly4PPydw2+4T85w4qAVgAMCaP3Zth8fWl7bvqj5lT+E05eHDyZEY/f+kxu4oHMksFc05JP2xEGmHG/PznDSoBbC5j6IpBYAHAVK+tL0MyVtq+O99MD8I1Zp/gCcGNKde4gAHdTvp4ermuyfU2mnOKm9CmrsCopEElgNZlbvL4YWJvApvq1d8eElMPrWw+TO3KDgBv/yQc/u7MtMMXup9igoa/1iS88juY/8H+53EArARgqpd/fP/mqcH7zDklN21KkQlAjCmVv07/wDPhzP/X/3MFVQGBG+RwmBQNACIyS0T+IiIrReQlEfmElz5RRBaLyGrvc4KXLiLyPRFZIyLPi8g837kWefuvFpFFg3dZpuolk8W3+wOAfywVgAU3wLn/CRf+KPdYKwGYStnvJPiHH8OCr8J7fzmwFwyDqoAAXn8sdyayIVJKFVAc+LSqPi0iY4EVIrIYeD+wRFW/LiLXAdcBnwfOBeZ6P8cDNwHHi8hE4HpgPu6NihUi8oCqBkziakatWA/88iL3Ys1Ft7heF9me+D78+auZY7Rnz6V60rXp5X2OymxMsxKAqaRKzSoXKnC73bA8/5vFg6hoCUBVN6nq097ybmAlMANYCNzm7XYbcIG3vBD4uTpLgRYRmQ6cDSxW1Xbvpr8YyOpDZUYVVfcKvX/Mk5fvhzcfh54O+Nu3co/p3gmP/GvuBB2zCnSvy34jc8ohwfsZM5y2vJi5/t5fwru/D59dMyw3fyizDUBE5gDHAE8C01R1E7ggAaQqaWcA632HtXpp+dLNaPXiPXDLO+GHJ8Cm513a6ofT2zc96ybI9nv5vuBzHXi6m8Rj8sFw+T2Z2/z/ecbNgMaWgefdmErzDx43YX83oNy89w3ryLUlBwARGQPcA3xSVXcV2jUgTQukZ3/P1SKyXESWt7UNX/cok0dPB/zhs/C3/yxej39falIOhQc/5aZyXPOn9PZ4T+5T0bN3BJ+reTKcfQNcuwzmnpW5bebb3Ngr0SY49TNlXY4xQ2belW5E0nEz4Yp7iu4+FErqBioiUdzN/3ZVvddL3iIi01V1k1fFs9VLbwVm+Q6fCWz00k/LSn80+7tU9WbgZoD58+fnBAgzzB75d3jaq/mbfnRwHb4qPHdn5mTcW16C1qdcAPFrXZ6eOWn7a26mryDNU/LnScQN+pVM2Cigpnq1zIaPP+P+XiXoeXjoldILSIBbgJWq+m3fpgeAVE+eRcD9vvQrvd5AJwAdXhXRw8ACEZng9Rha4KWZkSKZTN/8IX/3tXV/h/v+KTMt1gWvBvy6W31T7D3/6+DzRRqgbkzx/NnN31S7UKhqbv5QWgng7cD7gBdE5Fkv7V+ArwN3ichVwDrgEm/bH4DzgDVAF/ABAFVtF5GvAMu8/b6sqr4pnEzV2/hM5nrD+OD9Vv0xOD0oAKTmWE0m4bl81T9Tquo/jTGjRdEAoKqPE1x/D3BmwP4KXJPnXLcCt5aTQVNFVj2YuV7um7qpybPDda6aKBmD7WvcVI471+WfLLs5YERGY8yA2ZvApnTZT/adW4P3KxYY5pzi+u6nbFjhqo1SJmcN5Fao/t8Y028WAEzptr2auZ5vEKt8gSHlwDNg36N9512dGQCyZ1OyAGDMoLAAYEqXjGeuZ9/o433w/G8yu3oGmXxw5mTcnVthna/3z6Hvyty/zsb2N2YwWAAwpdGAHrmd2zLfBVj6A7j3Q5n7vOvG3OMmzMl8qm9bBXu2uOVoM0w9InP/7MBjjKkICwCmNMlEbpomoNs3lNOfvpi7z+yTctNaZsMY3+ieHa3p5foxudM5Bn23MWbALACY0miet35T1UDx3txtjRNg8tzMtLqxEG3IHN459fQPIF5f/uN97xEcn/VOgTGmImxCGFMazfMUvmcrTD0seGq7vk4IRzPTUtU5Y3xVQBkBwHsmecfnoaHFBZBph/c/38aYvCwAmNLkq4ZJjdnfuix3W6IvIM0rKeTr2ZMKAE0T4fQvlJdHY0xZrArIZNrTBmuW5I7S2ZWnb/8erwooKACkpsk7yjee+jxv9JC65vRk2n7Z9f/GmEFjJQCT1rsHfnKmm2D9yIvh4lvS2/K9pZtqA/CP6dM0GZomwSneyJwLvuLe+BWB0/81vd+YKbCjM/N8YgHAmKFiAcCkLf2hu/kDvHh3VgBYH3zMnjZXCkgdF2mAT62EiG/+07H7wNV/yT22ZTbseCMzTWxAN2OGij1uGScRgxU/y0zr8z2d+0sAU3xz7nZuzXz6n/7WzJt/IQecnpuWb4A5Y0zFWQAwzov3wK4NmWm7NqWXO3wlgBnHppc72zLr//2zHhUzd0Fu2hEX5KYZYwaFBQDjPPk/uWn+gOAvAcyYl17ekx0A5pf+ndOOyE176z+WfrwxZkAsABg3zEP21IwAuzamlzMCgO8m37nVzfa1d5uvdFCMCBx9eXp97gJonlT68caYAbEAYNxwDkF99lMlgEQ8szTgH6450Qfdvnl9GieU991v/4TrNdQ8Fc7+WnnHGmMGxHoBmcw3cf1SJYCO9b43ePeBuqb85yq3G+eUQ+DTr7jeP/YOgDFDygKAgd2bg9NTAaB9bTpt4gGFz9WfbpzZw0UYY4aEBYBa8Mbj7im/p8NNtpLd1dJfAmjZL92nP1Xt4w8Ak4oFAHuKN2aksABQCx78NLS94pZnnZAbANpfTy/POs4XADZCrAde/1t6e6oEcNaX4E/XZ55nwhx7mjdmBCkaAETkVuB8YKuqHumlTQR+DcwB3gDeo6o7RESA/wLOA7qA96vq094xi4B/8077VVW9rbKXYvLy3/B7OtxnrAceuxFCEdj6cnr7nJPdOwGadOP/fPfI9IBvkA4Ax/+TG7BtqxdYGsbBW97jevYYY0aEUkoAPwO+D/zcl3YdsERVvy4i13nrnwfOBeZ6P8cDNwHHewHjemA+oMAKEXlAVX2ziZhBExQAVvwM/vbN3H2nHg7jZ6a7fXZmzfubCgDRBph3ZcWzaowZOkUrbFX1b0B7VvJCIPUEfxtwgS/95+osBVpEZDpwNrBYVdu9m/5i4JxKXIApQVAAWPKl4H0nz4U5p+Q/14T9K5cvY8yw6m+L3TRV3QTgfaamd5oB+EcNa/XS8qWboRAUAIL6/TdPdf34DzyjwLnGVTZvxphhU+kuG0EVwFogPfcEIleLyHIRWd7W1ha0iylXdgDo3hE80fqUQ9znAacNRa6MMcOsvwFgi1e1g/fpDQpPKzDLt99MYGOB9ByqerOqzlfV+VOm5Jk1ypQnIwDsDJ6+EdLz9zZPdqN6Ziv2DoAxZkTpbwB4APCmdmIRcL8v/UpxTgA6vCqih4EFIjJBRCYAC7w0MxSafOPrdLTCpueC9/MP8RBUDXT+dyubL2PMsCoaAETkDuDvwCEi0ioiVwFfB94pIquBd3rrAH8A1gJrgB8DHwNQ1XbgK8Ay7+fLXpoZClMPSy9vfgE2PR+8X8vs9HJ2ADjvW3DAOyqfN2PMsCnaDVRVL8uz6cyAfRW4Js95bgVuLSt3pjKmHu6GaNCEe6s33pveNvFAaH8NGlrcOwAps47PPMeYqRhjRhd7b78WRBt91TsKu1rdooThgw/BRbfAh/+c2VYQqYejLnHLoah7g9gYM6rYUBC1Yp+joG1lZtq+x7gn+6MuDj7m3G/C9KPdBDBjpw1+Ho0xQ8pKALViztsD0k7OTfNrmggnXQv7nTQ4eTLGDCsLALXi0PNzh2ouFgCMMaOaBYBa0TwZ9jkyMy27odcYU1MsANSSIy9KL4/d14Z1MKbGWQCoJfM/6LqEhuvhHJt/15haZ72Aakn9WPjoExDrLjyvrzGmJlgJoNaI2M3fGANYADDGmJplAcAYY2qUBQBjjKlRFgCMMaZGWS8gk0NV6Y4lqI+ECYeCJnMrXTKp7O6J01gXpi5S2vNGd1+CN9s7CYsQCgmqSldfgq6+BN3eZ1KVw6aPY8qYeuqjIeojIUSC85q6nu17+tjVE2NnV4zNHT3s7I7xxJptTBlbz2HTx7H/5Gaa68Osa+9i3fZuAD5w8hzqIyGioVBGXnZ2x+joitEdS3DQlDE01YdJJJWGaDgwD6l/i53dMTq6YzRGw9RHQuzqiREJh2jb3Ut7Zy/xhLJhZzdvmzORcEioi7hrq4uESCahO5agqy9OTyz979EdS9DZm2D7nl427eph6drtnHXYNJ5Zt4Pj9p9IU12Err443X1JeuMJmurCjKmP0lwfJqmKKjTVR2iuC9PSFGVcQ5TZk5qYOrah/F94mbr7Emzb08ue3twZ6hJJZXNHD/u2NLLfpCZCIoRCEBIhLIIIe3/nnb1x+uJJ4kmlLhJi255ewiKEQ8LE5jqa6yt/q0sm3d9VOCRs3dWLCLQ0RWmMhomE03/rffEk0bDk/fscTuJGcK5O8+fP1+XLlw/pdyaTyp6+OPGEsq69i9VbdnPEvm6UzPpoiMZomMlj6ku+mSWSSjyZpD4SRlXpjSfp6kvQE0sQCQt14RCRcIhoWPbeZPJRVXpiSfb0xunqi3ufCTp743R0x1i3vYvdvXEaomE27OimszfOxDF1TB5Tz+QxdSSTSlfMu2n0JeiKJeiNJRnb4P5zvLxxF607utje2UdvPLn3e5vqwoRF6OyLk1S3vm9LI4fsM5YDJzfTl1B297ib6nOtO+mNJ2mqC9MXT9LRHSOpEA4J4xoixBNKo3ej6Ysn6fTyEkskSaoSS/T/7zF1s6yPhL3PEF19CXZ0ZV5PfzREQzTXRejojhFP5s/jzAmNzJrQxK6eGN19CfoSSfriSfoSSfb0xAseW21EQBVmT2yiJ+YCTW/M+3cUN89r6p4meDdk3E1ZvH1SHyK523vjycAbfzlC4gJCoX9XEZjYVLc3ONRHQkRCQiyhxBLud5NMKuObokRCIUICkVCIpCqrt+7hlLmTiSWSbN/TR10kRG88SXtnHzu6+sh3+6yPhGisC9Pdl9j7txcNC4dNH0dTXZjmugj1URfUE6qoKomkklD2LidV+fApB3DmYeUPxCgiK1R1ftH9RmMA+MK9z9MX171/HKGQ+6PriSWIJ1z6zu4Y2/f0oSg7OmPs6o6hQGdfPO8v1W9GSyNzJjcRCbmnuGgoRHcsQYf3hNcXT9JcH2FXj7dcF6Y7lqDY//9ISIh6AaEuEiIadn+Inb2JkvNmjBkdbrjwSC4/fr+yjys1AIzKKqD7ntlIdywxqN+xYWc3G3Z2F9zHn4fOvtLyE08q8WSC7tiAsldVxtRH+h289p/cjIgrdTRGwzTWRWiKhumNJ1i5aTedfXF64+4puxT7TWpi2rgGpo1rYEJTlJ1dMX73/EZOnTuFzt443bEEMyc08r9rtjO+MRr4O26IhmhprGN8YxQReK1tz97AnigS4cc1RBjfFKUnlqS7L8G4hggJVbbs6mVMfWTvE/Hh08eRSOreEkRvPEFIhMa9/w7hvf8mDVG3PLYhyl9WbeW0g6eyYt0OZrY0MnVcvdvfO6bOKxXt7onR2Zsg4pU4u2KuJNm2u5flb+4o+d9zoKJhYfKYesY1RAmqIQmJ7K0iSiRddVVCdW/VVZBISJgxoRFV6Ikl2Lq7N3jHCmiMhulLJIv+3vsrOcglxlEZAJIVfExuiIbYb2IzivLqlj0DPl9dOERTfZiGSJiEKn3xJLFE6qd4vusiIcbUR7x6XPfZXB9hTH2EmRMaaWmqo7M3zr4tjYxvjNLe2ce2Pb1s7+wjEkrfQJrq3M20PuxKML3xJAdPG8vB08YweYy7aWzr7GVnV4yWxij1kTA7u/vo6I4xpt5VhazavJvWHd001oUZ2xBhfGOUI/YdR19c2d7Zy0FT3bmi4RA93g0mEg6xpzdOR1eM+qi7lsa6MHXhkKvbDQkhX91uqZLezbLXu1n2xtxnQzTMhKY6murCec/5vcuOyXveVAl5Y0cP0ZAwrjGaU8+vqogIffEkr27ZzfbOPloaozTXR1z7QdjV4TfVhQu2EVTCv59/eEXO0xNL8NLGDhqiqb+XiGu3iKarPlP/zRR3M1bcv4X6t2vudneMu1G7INq/unFVJanu/3vYd45UVU/Kjs4+Yskk0VCImPc3Eksk9/5eouEQqro3yLhqW1f9M7Y+QigkRELCuIYo3bEE4xojTGyuY0JT3d6/7Y7uGBOa6lCUunCInliSrr44TXURImFh6+5eunrjdPYl9n6mgrr/7z4Ucu0bqZqLAyaP6de/TalGZRXQb59pJZ5wf3RJVe+JAaLeDTCWUFoao0waUwfAhKY6JjTXIQLNdZGCDZ+7emKsb+/a2xgZSyQZ3xglllAaoiHGN0YZ3xilPhqmsze+9wYQTypNdWGi4fxtB+rVgacCQp8XFARo9m72hY43xhio8SqgC4+ZOWjnHtcQ3dsoXMyYMnseiAh1ESm5gdkYYwZiyO80InKOiKwSkTUict1Qf78xxhhnSAOAiISBHwDnAocDl4lIZSotjTHGlGWoSwDHAWtUda2q9gF3AguHOA/GGGMY+gAwA1jvW2/10owxxgyxoQ4AQd1rMrohicjVIrJcRJa3tbUNUbaMMab2DHUAaAVm+dZnAhv9O6jqzao6X1XnT5kyZUgzZ4wxtWSoA8AyYK6I7C8idcClwANDnAdjjDEMw4tgInIe8F0gDNyqqjcU2LcNeHOo8gZMBrYN4fcNtdF8fXZtI9NovjYYvuvbT1WLVqFU9ZvAQ01Elpfy9txINZqvz65tZBrN1wbVf332yqkxxtQoCwDGGFOjLABkunm4MzDIRvP12bWNTKP52qDKr8/aAIwxpkZZCcAYY2rUqA8AIjJLRP4iIitF5CUR+YSXPlFEFovIau9zgpd+qIj8XUR6ReQzWed6Q0ReEJFnRWRoJysOUOFraxGRu0XkFe98Jw7HNfnyU5FrE5FDvN9X6meXiHxyuK7Ly1Mlf2//7J3jRRG5Q0QGfyb3Aip8bZ/wruul4f6dpfTj+i4Xkee9nydE5K2+cw3/yMjqTUg8Wn+A6cA8b3ks8CpuJNJvAtd56dcB3/CWpwJvA24APpN1rjeAycN9TYN0bbcBH/KW64CW0XJtvnOGgc24PtIj/tpw42i9DjR663cB7x8l13Yk8CLQhJu35E/A3OG8tn5e30nABG/5XOBJ39/ia8AB3v+354DDh/p6Rn0JQFU3qerT3vJuYCXuP85C3E0P7/MCb5+tqroMqPpZeSt1bSIyDjgVuMXbr09Vdw7JReQxSL+3M4HXVHUoXy7MUeFriwCNIhLB3Sw3BuwzZCp4bYcBS1W1S1XjwF+BC4fgEgrqx/U9oao7vPSluOFvoEpGRh71AcBPROYAxwBPAtNUdRO4XyruSaQYBR4RkRUicvVg5bM/BnhtBwBtwE9F5BkR+YmINA9idstSgd9byqXAHZXO30AM5NpUdQPwLWAdsAnoUNVHBjO/5Rjg7+1F4FQRmSQiTcB5ZI4jNuz6cX1XAX/0lqtiZOSaCQAiMga4B/ikqu7q52nerqrzcEW5a0Tk1IplcAAqcG0RYB5wk6oeA3TiirHDrkK/N8SNPfVu4DeVyttADfTavHrmhcD+wL5As4hcUdlc9s9Ar01VVwLfABYDD+GqSOIVzeQAlHt9InI6LgB8PpUUsNuQd8msiQAgIlHcL+t2Vb3XS94iItO97dOBrcXOo6obvc+twG9xxbhhVaFrawVaVfVJb/1uXEAYVpX6vXnOBZ5W1S2Vz2n5KnRtZwGvq2qbqsaAe3F1zsOqgv/fblHVeap6KtAOrB6sPJej3OsTkbcAPwEWqup2L7noyMhDYdQHABERXN32SlX9tm/TA8Aib3kRcH+R8zSLyNjUMrAAV0wdNpW6NlXdDKwXkUO8pDOBlyuc3bJU6tp8LqNKqn8qeG3rgBNEpMk755m4OulhU8nfm4hM9T5nA/9AFfz+yr0+L+/3Au9T1Vd9+1fHyMhD3eo81D/Aybii1fPAs97PecAkYAnuqWIJMNHbfx9cdN4F7PSWx+HqyZ/zfl4C/nW0XJu37WhguXeu+/B6LoySa2sCtgPjh/t3NgjX9iXgFdzDyC+A+lF0bY/hHkSeA84c7t9bP6/vJ8AO377Lfec6D9eL6LXhup/Ym8DGGFOjRn0VkDHGmGAWAIwxpkZZADDGmBplAcAYY2qUBQBjjKlRFgCMMaZGWQAwxpgaZQHAGGNq1P8HuqP5BbSCbIkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(steel[:,0],steel[:,1])\n",
    "plt.plot(al[:,0],al[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}