{
"cells": [
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdaUBU9cIG8DNnVoZt2HcQFFRcURTEXUER1zTL1NRMy8y08nZvq3VbLLWsXEozNXPX1HDFDQUREGVTQEFWZd9h9uWc836gl7xmqTBwZnl+n2bODHMe/87MM2fnMAxDAAAAmCuS7QAAAABsQhECAIBZQxECAIBZQxECAIBZQxECAIBZQxECAIBZQxECAIBZQxECAIBZQxECAIBZQxECAIBZY7MIKyoq1q1b9+TP12q1HRfGVGHQ2gCD1gYYtDbAoLVBRwwam0VYUlJy9OjRJ3++SqXquDCmCoPWBhi0NsCgtQEGrQ06YtCwahQAAMwaihAAAMwaihAAAMwaihAAAMwaihAAAMwaihAAAMwa70meJJfLY2NjMzIylErl6tWrW6fTNL1169aYmBhXV9eVK1cGBAS0TC8qKlq7dm1ZWVl4ePjrr7/O5XI7JDsAAEC7PdESYWZm5tq1a7Ozsx86/v2bb77ZtGnTq6++6ubmNnLkSKlUShCESqUaOXKkjY3Na6+99ssvv3z++ecdEhwAAEAfnqgIw8LCrly58sEHHzw4kaKoDRs2fPfdd1FRUZ988km3bt327dtHEMThw4cdHBzWrFkzYcKEzZs3b9q0SaPRdEh2AAAwJxqa2F3EbVDr+WXbvo2woqKitLR02LBhLXeHDRt2/fp1giBu3LjROjEkJKSpqamoqKj9QQEAwJyduEf3PawuTbvO0Sj1+8pPtI3wkaqqqiwsLCwsLFruOjg43Lp1q2V69+7dWyaSJGlvb19ZWdk65UG1tbW3b98eM2ZM65TXXnttwoQJfzdHuVzO4XAYhlm6dKlCoWhzcoPC4/E2btwoFos76PVbBq2DXtxUYdDaAIPWBhi0J1Qg5byfyvQqvHi69ijX3oGr8JdxqCf8W5FIxOM9punaXoSWlpYajYamaZIkCYJQKpVWVlYEQYjFYrX6zwVXhUJhaWn5yFewt7f38PB4//33W+6SJBkUFNTyIo/EMIyVlZVGozlw4EDLalgTsGTJEq1W+w//6nZqGbQOenFThUFrAwxaG2DQHktFEWtTVVUJZ7+t+03p3iPdcaZNmUsPLcdLr+PW9iJ0d3dnGKa0tNTb25sgiOLiYi8vL4IgvLy87ty50/Kcuro6mUzWMv2vSJK0sbEJDw9/2lmTJDlz5sw2Jzcob731FtsRAAAMUVyp9vixswvKD1W59r1lNdfqvoTbQ+o32dnL002/M2r7NkIbG5vIyMjt27cTBFFTU3Py5MmWcpo5c+bZs2fLy8sJgtixY8eIESNcXFz0FRcAAEyeVMNsOnCZ3Lx0qKKwiLdQWTvYprfliE/6TVsQ7uxor/fZPdESYW1tbUBAAEVROp3O3t7ezc0tOzubIIivvvoqKirq4sWLhYWFzz///ODBgwmC6N2796uvvhoUFBQQEFBYWHj8+HG9hwYAAFOVfCOrLnqnncBPzrxMUHKXCZLBQ/p16MbUJypCBweHgoKC1rstGwUJgujTp09+fv6tW7ecnJx8fHxan7BmzZply5ZVVlb26dNHJBLpNzEAAJgkZV1N7C/bm5olTurnhQ5NXWY69ugR2gnzfaIi5HA4dnZ2j3xIKBQGBwf/dbqXl9ffbRoEAAB4EEPp0o8dyM2SOyrG8r3rA6Z5+HgO7rS5t31nGQAAgParzrpx8Wiag6y/yq229xJvN9dBnRwARQgAAOxQNDXEbI+2quxZb+vS5RXvl7p13lLgg1CEHYhhmPj4+NTU1Orq6jfeeMPDw4PtRAAABoFhmIuHT2pS7WQWkmthNv+dHCJi7+oMKMIOpFarly9fHhwcvHv37pkzZ6IIAQAIgriTnZtzqIjUiWLcm4ZPnLKmG8sXBMT1CPVj9+7d1dXVrXf37t1bVVUlEokyMzO3b98uEAhYzAYAYCBkcsWxzSdLfm1otiz5rM+gZfOi5rDdggSWCPXl4sWLxcXFH330EUEQt27devPNN+/fv892KAAAA5IYn9Z4WqkRNlzvxr/l83LsKJ6tYSwjmEgRzo+jyhVMp83ug/7cUW7/c3Tna6+9NnPmzPfff5/L5f70008LFizAAZQAAC0am5ov7kgSVVnw7eO+c10Y0cs1egCXNJjzjZtIES7vRer9ClX/oL/Dw/+BISEhzs7OZ8+eHT169L59+5KSkjovDQCAAbuWlNFwXEmJyux6aJ/l/mvDMOGzvuyvDn2QiRThQEf2f1q8+uqr27Ztq66u7t+/f0BAANtxAABYplKrT/8Sb1FsKRRFK/tPeUEZciqCG/SXBQnWmUgRGoLZs2e/9957eXl5q1atYjsLAADL8vNLcn4tpgRSf9sTG3v96xrjfm0qz9WC7ViPYljLp0bN0tJy9uzZ1dXV06ZNa504atSorl27KhSKqVOndu3atbKyksWEAACd4+KZxMJtVRzL9KFut1/t8UWVpcelKANtQQJLhPqlUqkWLlwoFApbp0RHR9M03XrX1taWjVwAAJ1EpVaf+jleUC5ytjroGBgyQj5zthf534Fcg1sf+gAUoX4UFBQcP3788OHDWVlZD05H8wGA+SirqLqx9TYjUgeSu2tHLBlaOujrEO5cAzhS8J+hCPWjpqamtrY2JiYGp48BAPOUeiOr9jcZ4VkWVnfuatTHywu9D43ljXA15EXBP6AI9SM0NDQ0tDOumwUAYIDORscTSUKRT0qQqnDbmG92lUriJnEDbI2gBQkUYYcqLy//9ddfU1JSGIYZPXr0kiVLcK41ADAxOkp34udLvPtCH7dTrkLhv/y/LGzkJ07mORrPOUVQhB3o/PnzxcXFc+bMIUly1apVt2/f/vHHH9kOBQCgN1Kp/PzmJIJmelrut/Ia+Cx/vgeXvBDFZfFSEm2AItSPjRs3Pvvss25ubi13N2/ePG3atPnz58+fP79lilgsnjdvHooQAExGWVll2tY8xlEVUndAM3zW8PrI5705nwUb9A6ij2ToO/MYi8zMzK1bt7be/uyzz5ycnB58Qn5+vre3NxvRAAD0Lzs779bGYo5P3ZDqPWURy8Oqxn/Qn/zcCFuQMJklwrodn1EN1Y9/nl5wOLaTFgoD+j84bdmyZZMnT/7www95PN6WLVteeumlBzcHFhUVrVq16sCBA52UEACgIyUlpElP6AS9i/rln0yI/O8bJX4HxvAeuhSBETGRIrSdvJBWKTptdnxnz4em9O/f393d/fTp02PHjt2/f/+NGzdaHyorK4uIiPj4448jIiI6LSEAQAc5fyKBvsqT9LrlU3p926h1uyuc441nB9FHMpEi5Dmxf/Tea6+91nLS7ZCQkG7durVMrKysDA8Pf+WVV5YvX85uPACAdmII5vivsbxcoU+3BMeGyhV919YqrZKn8OyFj/9bQ4ZthHrz/PPPJycnr1u37pVXXmmZUlNTExERMWfOnH//+9/sZgMAaCcdpTu26TxTRAa4nLLhaKPcP7OXWMdEGn0LEihCPbKwsHjxxRebm5unTJnSMmXHjh3Z2dnr16+3/386nY7dkAAAbaBQKKO/vkQ0E/0Eh5WOXUOFKxf3Fv44lMs3iQ4xkVWjBqKxsfHll1/m8/ktd1esWNG6dNiCx8OAA4CRaWhoituYSlhqByn3pwU+95oucv9Y3kij3TXmr/C9rB85OTlHjhyJjo7Ozs5unSgSiUQi4zm5AgDAX5SVVaZvucs4N4dUH9zZ641oweDkKK63lem0IIEi1Jfm5maapmNjY11dXdnOAgCgH7l5hUU7qxjv0uCKU290+8TevVvCUCM7a8yTQBHqB066DQAmJvXGrdrfFKT3Lf/6jKguX78Z6vxyd5PYJPgXKMIOVFtb++mnn+bk5Oh0uuDg4P/85z8PnW4GAMAwXb6QojlPCF1iKanqha5r9o6z6mtvUqtDH4Qi7EAajcbPz2/GjBkcDmf9+vXPPPNMQkIC26EAAB7jzOF4Thrf0vrIdbr7zZCXrgzlWfHZztSRUIT68cMPP0ydOrX1qrwt5+B2d3d/8803W6Y4OTn17t1bp9Nhx1EAMFg0zUT/fJF3j+9M7vzJZkbExIi3u5rm6tAHmf6/sHPcvHmz9coS169fX716taOjY8vdxsbG/Pz8DRs2TJ48GS0IAAZLqVIdW3+BLtc5Ub/+3GvFxwvHzTKDFiRMZonw/bjPaxX1nTa7Rf3nDnYb8OCUpUuXjh8//uOPP+bz+T/99NOiRYtajyYcNWpUWVmZUCg8fvx4pyUEAHgqNXX1CZsz1EyjmErJnL72p0EOJrtJ8C9MpAiXDVwk1cg6Z15cDult8/BJt/v27evn53fixIlx48YdPnw4IyOj9aGW2/v27QsPDy8qKrK1te2cnAAAT+hufnHe9jKl8I7ChhO6cLWfxKQ3Cf6FiRShuxX7R++1nnR76NChXbp0eejR2bNnL168ODc3d/DgwWykAwB4tMvx6eqTGoXFRbvQ0S9EDjefBcFWJlKEhmDmzJkrV67My8tbv359y5SioiIHBwcbGxuCIA4ePMjhcLp3785qRgCAP9EM8dP28z55lmrbs6MXLnF0c2M7ETtQhHojFArnz5+/b9++iRMntkxJTk5eunSps7OzVCoVCAQHDhzAelEAMBBnilUle867KoR2fXJCZn/E4ZpvHZjvv7wjVFRULFq0qHXX0BdeeGHGjBllZWUikcjNXH9qAYChuVHLfHupcHpGkR1PNmqOu6TXOLYTsQxFqB85OTn79u2LiYnJycl5cLpAIPD19WUrFQDAgzLrmU9SaZ+s0/MqHFXu955ZMosUidkOxT4UoX6o1Wp7e/v4+HicRA0ADFBKDbM6gy6uqFlZds6xsav9mOZBkQvZDmUoUIT6ERQUFBQUxHYKAICHnStj1mZSRU30al00P9eKFtgPejvA0QU/2f9kFmcNYF1sbOxzzz1348YNtoMAgLnQ0MSuu3T/o7p/XaPeFGd9X/Kd3S1vfg9mxqqJaMGHYImwwzU3N7/55psVFRVz585lOwsAmL4qJbH1Dr3lNtXXnvNd11KPpD0Zmd1F1GDvl516dA9hO50hwhKhfnz55ZfFxcWtdz///PPS0tKW2//+97+XLVsmkUjYSQYAZuNGLTPvMtXzN225nIkNrdlX/700endhUQTfx27cR2E9uvuxHdBAoQj1o7y8fMuWLS23k5KStm7d2nK8xOXLl+/cubN48WJW0wGAKdPQxP4Ceshx3XMXqT72nLtja768/73i109j83sy1HDP+fbTFoYLhQK2YxouE1k1mvVjkape2znz4pAcv2mudj2tH5y4dOnSkSNH/ve//xUKhS0n3eZyuQqFYtmyZYcPH+ZwzPCkRQDQ4aqUxJbb9NY7VKCE859+5ASyUB57uKEgK8Nuko1imEWwYvyMYTwzPlL+CZnIAHWf60Vp6E6bnfAvZ6Tt2bNnz549o6Ojx48ff+TIkaysLIIgVq1aNXHiRFdX14aGBoqiZDKZUqm0sLDotJwAYKpu1jPfZdHRJfRMX/J8JLdrxXXpyWO1NeW3vCYQ2tEER9rnX36uzo5sxzQOJlKEfGse6ydLbznpdlVV1ejRo729vQmCuHPnTmJi4rZt2wiCaGpqevXVV/Py8j755BOWgwKAMTtXxnx9k8ppJJYFknlTNMKMc7IfTjSLrW57jGquctSV67rOkvTtF8p2TGNiIkVoCKZPn/7WW2/l5eX98MMPLVNOnjzZ+qi/v/8333wzZcoUltIBgHGjGeJIMf1lBq1jiJV9yJmiInXSKeXBK5weAwv6z6pJ45PZpMNYy6EjB7Kd1PigCPVGIBAsWLBg3759kZGRf33Uxsam9VK9AABPjmaIg4X0Z+m0rYD4rJ9uVO0V+anTTc11FqET8oa/0XyNIfIJu5EWw0cPIknsjtAWKEJ9KioqWrx4MZfL/etDqampnZ8HAIwaQxDHiulVqbREQGzpXtG34Ixi+wWlTw9y+IzUQoq5zNNYqJ1G24aNGIA98toDRagf6enpu3btunz5cut6UQCA9rhSybyTQlE66ie7FP/bp3TX75Eh46XPfph9rcr6oFDu2hjwvEvffj3YjmkKUIT6IRQKe/funZiYaG9vz3YWADBuRVLmnRQ6v7xxk+Bc19uneA5uvNDIm7Ui6Q216IqCG0j2WdnFFadJ0x8UoX4EBgYGBgaynQIAjJuKIr7KpE6nlazRHA8oTRD3HVo27vW8mzLr32ya7Zqch9kOGRGE4wL1DgMKAGAQzpcxG2NyltYcXiTPZQZNuOmwVJvD5d1keD3JwBVeHu6D2A5oslCEAAAsa1AT352+FZSx71umoiJwwrWyMEmCg8K93ivKMXhQH+wL2tFQhAAAbIq9kVd/clck01DpMiav1F1zW2MZRAYvDrC1sX78H4M+oAgBANihqKtN/GmTprFBJxoua3KlmUbvOY69egWwncvsGF8RkiRJEETXrl3ZDqIf1dXVAgHOCg9gXhiN+k70gexbJVxqKJ8jtO1Bh44LtLa2ZDuXmTK+IuTxePfv35fJZGwH0Q+xWOzkhN2gAcxIQ9qV8ydTrBShOrG7xzi70KFBOByeXcZXhARBODs7Ozs7s50CAODpyGoqzu06Ia7p02zZ02mm4/AAd2trbAhkn1EWIQCAcaEo6uyuo2SuR7OFS/IQ+88mDxFyCalUynYuIAgUIQBAR7t2Obn6nELNFR3xoiPHT17rT7KdCP4HihAAoKOUlVUm/5oubLaody3+1Hn2oQjhQEdsDjQ4KEIAAP3TUbpzh6/w0ywZ67uJfYdetpyXGMFzsWA7FjwKihAAQM/u5hdn7ymiaJmrd8Ia19c9nGwvDeMKH3F9NjAIKEIAAL2haeb0oThBupi0SgoM8h2vePflrtwPg7BR0KChCAEA9KOysjpx+y1Cp/MS71eMXTK8qPf3odxZXdGChg5FCACgB0lX06THtYTz/SH09YTwT5bnOP42ljfcFbvGGAEUIQBAu1AUdfLXS4I8sZVL7EAJvbPfuh/yBLETuT0laEHjgCIEAGi7hoamSz/eIBjG12qPa8CgVXYvXrlHXJ3McxOznQyeGIoQAKCNcvMK83+pIDzlg6sOWEfOX6wIr2ti4ibxbPhsJ4OngSIEAGiLxPhU2SlK0K8mKOeA6IV/Ty/p6yAkTo/n4TAJo9OuIiwrK1Or1a13RSKRu7s7QRDFxcU0TbdMtLKywgmyAcDEnP4tjkwVOg4u8bh1krPwi/BMnzAXzoYhXFxM3hi1qwiXLVt28+bNlttVVVXjxo07evQoQRD9+vWTSCQ8Ho8giKlTp65fv779QQEADAFFUdHbYvnlwm4Dc2zykpQvfz0uyWGeP/kRDhY0Wu0qwmPHjrXcoGm6S5cuL774YutDcXFxXbp0ac+LAwAYGqVKdXrjFY6W7N8jXXi/qOrFrydcsfywP7mkJ1rQiOlnG+HZs2fVavWkSZNapxQXFyuVyq5du+Ly6wBgGhqamuI2pHKsiTD7RF6z6vaML5+N4/0QRs7wRQsaN/0U4fbt2+fNm8fn/7GnlEAgWL58uVqtbmho2L59++TJkx/5V1qttr6+/tChQ61TQkNDPT09/24uNE23bnqEJ4RBawMMWhuY/KBVVlan/pjHeOuGqmJ5XJv4Uateu0LuH80d5Uq0+R9u8oPWEZ520Ejy8T9T9FCEdXV1J0+eTEtLa51SUFBgY2NDEMSuXbvmzp177949W1vbv/6hVCqtra09ePBg6xSxWOzo6Ph3M1Kr1a1dC08Ig9YGGLQ2MO1BKykpK9pdw/RUDqk6Sbr77u/12pfXyGOjdP0kWpWq7S9r2oPWQZ520AQCQcsOK/9AD0X466+/DhgwIDAwsHVKSwsSBDF//vwVK1ZkZ2eHhYX99Q/t7e0DAgKOHDnyhDOiKEosxkGqTweD1gYYtDYw4UG7cye/ZFcdf7BqYP7vop6DvnWbty+PuDKZ62fd3g4z4UHrOB0xaHoowl9++eWNN9545EMVFRUymczJyan9cwEA6HyZmbcr9zcLhyt7pu0WDYlaKZx+s4y5OpnnJGI7GehPe4swJSWloKDgueeea51y6dKlo0ePBgcHy+XyzZs3T5061d/fv51zAQDofCnXbjYeVVmPUfombReMnT1bMY6jYGKjeFZYnWla2luEWq12x44dretCCYLo0aOHi4vL5cuXxWLxe++998ILL7RzFgAAnS/xSprspE4yVuZ1ZRs98dXwimGDnDibw7g87CJqctpbhEOHDn1oipub24cfftjOlwUAYFHcxRTVecY5vMklbkfjlLcjCwcs6Um+1w8daJpwrlEAgP9x8UyiNo70GF3vGL87d/KHz+V13xhGPu+HFjRZKEIAgD+dPZ7AJHG9h5fZJR49Ff75J4U+J8ZxQ5xxClFThiIEAPjD6cNxnHRBl+C7kvTLa0LWXGlwTZ7C9bREC5o4FCEAAEEQxIl9sWSOoGuvdHFBztyANa4W9lfGccX4jjQD+E8GACCid14gCwX+XvF0tXS02+rXA63e6o2NguYCRQgAZo1hmGM/XeCW83s4nC5V2y10+XjXGOEIV6wONSMoQgAwXzpKF735EtlI9uDtP0uEnuw299povrMF27Ggc6EIAcBMqdTqUxviCQ3VU7X7G5dZ3UZHnu5D4hLzZghFCADmSCqTn/8+meaoesgO/rfbynem9h/giA40UyhCADA7tbV1iRuzlIIaD+ml4+Ff7RzlKeKynQnYgyIEAPNSXFKatbW4UXRHxFRZL/v2Aw9cCMncoQgBwIykpN+uP9AkF11x7Oo3btYrPC5WhwKKEADMA8UQW35L7nadoxGfGPbs8x69+7KdCAwFihAATBxDEIcLqeyj50JqbUTu8RGL3iGtJGyHAgOCIgQAk0UzxOEi+usbsiV3LgWphT1HVvpGvktwsDoU/geKEABMkEJH/HqXXp9FhymzPsit0QqV4a/1sfL0YzsXGCIUIQCYlBIZ8+NtekcuPdKJWtd4UHQ3QONXO+2V5wkS5w6FR0MRAoApoBji9H36pzt0cjUzz59M6paSdj5bIO/lPJUbNPQFttOBQUMRAoBxy2tifsmjf81nfKyIxd3JfQH5Bcf336oKIS29Qv7T297Olu2AYOhQhABglBo1xKFC+te7dEEzM6cbeW4C119Z0nRm94UykZV0gkWIatz04RzsFwNPAEUIAMZERxNny5jdd+mzpXSEB/luP26kJ4e+f0d65NDN0soSZipJCrq+5urn6812UjAaKEIAMA63G5kdufTeArqLFWeeP/njUL4dn1LeSmj4/XdlU2OqZLylYoJFsHL8jOG4hAQ8FRQhABg0hY44WEhvu0PfkxPzunEuTeR1t+VQjbXyy2crks7w7F1vOg1TV7rQYnnvN73d3VzYzgvGB0UIAAbqbhOz+Ta95y491JV8vz93gheH1GlUWVdrr5/XlOSKg0aUDn+19KqOW8dxf8Y6eFAY23nBWKEIAcDgXKpgvrlJXa9lFnUn06fzvMSEOj+z+eAl5c1EgZe/eHB4cb8Z9y/WCeQcyXDRyHGhWBcK7YEiBABDwRBEdAm9OoOWaYm3+5C/jeWQpbmKc5crMq6Q1nbi4DGO4+YkpRc3nJTztDKrMNHo8SE8Lr7EoL3wHgIAg3DiHr0qleZyiA/6k1H8ElX65caDlwm+UDxgpNOytfWM6PL5m7xzBSoLpeNw66GjgrEUCPqCIgQAliVWMe+kUDItsaZ7XVjFZeWB2AaNymLAKIdFn9AOHteu3qzdedu61lbro/N/0alHz25s5wVTgyIEANaUypl/XaNTK9Wb7K4F3T+nSy2g+4+we3456RmQej2n7EiB9f1aqW2zbZBlyKheVmJcSh46BIoQAFigpYnvsuh918tWU2e+uXdR6NnNcsgEImBAempe5YkK61K5zEou7MnrM9PH1dWZ7bBg4lCEANDZUmuZb09kvlBx7Lj8jk3oeGr8lzfvVNfHSm335TZJGkXd+YEzfTxwRCB0FhQhAHQeLU3sPHHN4/qBVbxmy+FTs6RjFdlam/iGRkepTW9x0MKu9na4djx0NhQhAHSSu+np94/90oOgGruNzqqQ2J6VSF3rnYJtBob1xPY/YBGKEAA6nLayJGvfzwWNjEoU7tborq2pdxhkHTKsl4VIxHY0ABQhAHQkWqUoif71+m2FSD2Ox6ccupMDxnazk+AagWBAUIQA0FEK485fv1RiLx/YYF83cJpzcP/ubCcCeAQUIQDo3/38gqSDKZKmrkV2ztUTfN4IGcR2IoC/hSIEAH1qaGy69Osl61J3pQP9eaDPD5GDhrniXGhg0FCEAKAfFEWd/z2BvCbSWTbeGOy6n3z+xDiujxVaEAwdihAA9CDvblFZdL2Wlnr7pfzk9Yqca3l1LM+az3YsgCeAIgSAdlGrNWf2XrHIteRYXhkyvM8U2VsDbDkHwrg8ku1kAE8GRQgAbXcnt/Du3lKdUOYl3m8x4+0R2QHz/MlVQehAMCYoQgBoo5ij8eQ1C55HVogutzBy1bhM94+DyMU90IJgZFCEAPDUpDL5ua2JvGZ+F7doL1vL22PXPnOZs2UYOb0LWhCMD4oQAJ5OYWFJ9o57hItuAP+Qjf/Q5L7z5sVRP4dopnQRsB0NoC1QhADwFJIT0ptPaIVBsuDbe2zGzznnHvlaPPV7BK+3hYrtaABthCIEgCd1+kgceV3oPLbZNX6X5PkVx8Qh/0qkYiJ5/R04Uinb4QDaCkUIAI9H00z0zxd494U9oqSiczscFn64Vxf4cQp9IYoXKMEh82DcUIQA8BhqtebkpjhSRQ6OVNJndzku+eJnqe9XmXRsFNffFi0IRg9FCAD/RCaTn9+QTIiIkSOU6vMHHV//6ocaz++z6csTub7WaEEwBShCAPhb9Q2NV77PYJzpiECl4tIRp2Vrvy132ZZLx03ielmiBcFEoAgB4NGqa+qSN2QzvroIP6k84YTTsnVf3XPck6/1bDAAACAASURBVE9fiuJ6oAXBhKAIAeARyiuq0jbf5fSkxro1yq+ecVq27r8FdkeL6UsTea4WbIcD0CsUIQA8rKysMv2HfLKvbrhjjTzprNOytR/kSWLuM5cm8pxEbIcD0DcUIQD8j7KyyszNBdwB9DD7WnnSWcdla9/OkVytYmIn8uyFbIcD6AAoQgD4U1l5ZebmAs4Aaph9tTzpjMPra5fekmQ3MBejeLY4gRqYKBQhAPyhorIqY3M+OYAeZlclT4qxW7r2pQxJuZw5G8mzwiV2wXShCAGAIAiiqqYubdNdTl/dMPtaeVKM7WtrZ92QaGnm1HieBb4nwKThmikAQNTU1ad8n00Haoc71cqTzoiXrJ2cbCvkEkcj0IJg+lCEAOauoaEp6ftbjL9ulGudPOkMd/HaiATbbjacvaO5AnxDgBnA2xzArDVLZfHfp9M+2tGedfKkM6oFa4bH2Y7z5GwZxuXioHkwD1jrAWC+ZHJF7LfXGTdduE+9/Nr5qjlrJsVL3u1Hvh6In8hgRlCEAGZKoVSe+y6JsKfH+lTJr8ffmr5m7lWbrcPIaT5oQTAvKEIAc6RWa2K+v0pYMaN9y1WZSccjvvwkzSY6ghfqjPWhYHZQhABmR6PVnvw+jhAQo7wK1Lez1g768sI968TJ3C64rBKYJRQhgHnRanUnvr/M4XCGu91Sldx/udvnIsr66mSuNQ6ZB3OFjQEAZkRH6aI3XOJQnDCHlKbqmjHOnw7xtjkWgRYEs4YlQgBzoaN00d9f4mg5IdaxBc3iWY4fbRspivTE6lAwdyhCALOgo3S/f3+JVBNB5PGLKv8dAa8mjeF74vq6AChCAHOg1eqOb7jEUTOB6gNbrCOsxzx/sT+J4+UBWqAIAUycWq05+X08rdV2b963yfelhc+OGeCIDgT4E4oQwJQpFIoz3yepdIruskOJY99bH94Lpw8FeAiKEMBkNTY1X/juhoyo9WcS7N5Yv8TTme1EAIYIRQhgmmJyqhV7b6v4RQM8mvos+IbDxwXmAR6tXUX4n//8Jy0treW2h4fHL7/80nI7ISHhww8/rKioGDNmzLp166ysrNqZEgCe3MVyZtuVwhfTKylR2rQJva1CIthOBGDQ2lWEaWlpoaGhI0eOJAhCLBa3TGxsbJw8efK33347atSo119//Z133vnxxx/1kBQA/pGaIg4W0t9n032qUl4qJHV2V5956QWeixfbuQAMXXtXjfbp0yc8PPzBKfv27evTp8+CBQsIglizZk1oaOjXX39taWnZzhkBwN/JaWS259J78ukBDpy3K485Fnpa9c4Jm/M2h4ttHwCP194dyNasWTN06NAlS5YUFRW1TMnOzh44cGDL7d69e1MUVVxc3M65AMBflcmZb7PoQb/rxp2hRFwiIaTspRsHbEoc/Wfxhs5biBYEeELt+qjMnz/f3d1dIBDs3r17yJAhWVlZjo6ONTU1gYGBrc+RSCTV1dW9evX6659XVVWlpqba2dm13OVwOOvWrXvuuef+bnYymaw9ac0TBq0NDHzQiuWcE6Xc46Xk3WZOlAf9US9qhE1z6dmjWbl+HLHNoOU9LK2spFJpJ6cy8EEzTBi0NnjaQROJRHz+Y86l264inDt3bsuNoUOHJiUlnThx4qWXXpJIJHK5vPU5zc3NrVX3EBcXl379+p07d651irW1NY/3T5Gsra3bE9g8YdDawNAGjWaIazXMiRL65H2mRslM8SE/CSbHuHP4WqUs/sSVxFtcxWRBP9m0FyI5HNYOlje0QTMKGLQ20Pug6WflCYfDcXBwaClqPz+/q1evtky/f/++RqPx8fH529nzeH9XkwAg0xJnS+kT95jT92k3MWeSN2fbMO4gJw7JIWh5s+zc8XuJ59JFUZbqKLfZkr79QtnOC2CU2l6Ecrk8IyMjLCyMw+EcOXIkOTl506ZNBEHMnj37888/z8jI6Nev39dffz1p0iRUHcBTqVYSv5fQx4rpxCpmqCtnsjf56UCet9Ufi3ra0vymhJPKmwl3fcbIta9ynRUjlve3ssL+aABt1PYiVKvVCxcuLCkp4XK5Hh4eu3fvbtkQ6O3tvWHDhpZdSf39/Q8ePKi3sAAmrUFN/FZEHyik0+uYSE/y5e7k4bGk1f9v3aBljYq0OHnKeUYhlQeNz3B8xaLEynGKePCQoaymBjB6bS9Ce3v73NxcjUZD07RIJHrwoYULFy5YsEChUOBQeoDHohniYjnzcy59rpQe50m+0YuM9CRF3D8epZrqVFnJiswE7f27ot4hwvHzrtyUia5Y8XpKRy0JEostWM0OYArau41QIHj0eZtIkkQLAvyzZi2xI5felEPbCohF3cmtw/iSls8TTWuKc1W3b6hyruvqKkSBg6yGTSL8+l46m0bupTXOVN/lHj6eHiynBzAVONIIgAXVSuK7LOqnO3SEJ7lnFDfUmUMQhK6mTJaXrs5NV+ff5No5i3oMtJ26WOgXWNckvXA6XbAvW+WoDZzv0aMndooB0CcUIUCnalATa25SP9+hZ3UlU5/heZFydV56Q2yqKi+NoChh9wEW/YZKZi7jWtsxDJORll246aKk1EHjpeuxyNW/Wxe24wOYIBQhQCfR0sQPOfTqTGqaD5k+vFySn6zakVJZXiTo2kfUfYD16Bmt5wXNybmbdy1NcNdCx6VEffgDXuxmZ2fLbngAE4YiBOgMlyuY169SgzllSZbx4itXGLWS6jPEetxsYbe+HB6fIIhmqSzz0o2anEZxqaWOpIiuhN+Lzj16dmM7OIDpQxECdKxGDfF+opxzM/6I8pytvMpiwCjx7LcF3t0JDqe+oSkt8Vb13XqilGcjs250aBT68QImuPv5/e05KABA71CEAB3o0p3qG9HH3qq/aNOtt234bFGPgUX3ym5k3Ws+fkFQLbJQiZscGnhepNdkp169/f9uH2wA6FAoQoAOoaypiNu33+t+cuSAcWTUqoyCBllMk/Wu6xRJyV0UVj7CLhFu/v6+XC738a8FAB0JRQigZ7SssfTEnqb0hAL38VLfly2zrDm3pAo3laSHZY/nfTw8XNkOCAD/A0UIoD80JY37/Xrc1QJOqA/xtnuznO9P+I9z8/f3ZTsZAPwtFCGAftTmZiYcieNK+2k4U0U96F7h3h7uWPgDMAIoQoD2qqmsit8Ta1PVrd7WN6aHw/pne3lasnZRQAB4WihCgLZrapLG7k+wKnSgJNQvA2xUTgN2j+KK8akCMCr4yAK0BUVRF04kEkkCjUWVy2jNSuWsse6cb0K4JBYFAYwNihDgqd25U5B74D7FKD0ksX7PLhuf5vFaT/KdviTbuQCgLVCEAE9BR+lO740TZVuTjmlD7aXlE/87KlbwRTC5IAAtCGCsUIQAT+peaXna9lyKT/nZ7PLuHXJn4JtTzlObw8gZvmhBACOGIgR4Ilfj02SntbwessGF+yVTFmV4jX7mvG7nCF6UF7YKAhg3FCHAY9A0c3J3LD9X7Dyizi1xr8P8965b951xQbdnFC/CAy0IYPRQhAD/RK5QxPxwlVTy+ofX8uMOOL62+hrXd/oF3f4xvLHuaEEAU4AiBPhbtbV1CZtvceyIkQObNFeOOb6x7gbtNuO8bu9otCCA6UARAjzavftlt34sJntQYz0a5YkxTm+sS6ecpp7X7RqJNaIAJgVFCPAId3ILi3dW88OYEeJa+bXzTsvWZlP2k8/qfh7OjfRECwKYFBQhwMNy7xTV/aa0GEcO4FTKk887LVubz9hNiKE2hnEne+NICQBTg081wP+4mXG7/jeV9ST+QF6ZPOmM0+tr7nHsxp2hvhpMPovjBQFMET7YAH+6mXG7cr9UFMn041bK4qOdXv+qkmcffpp6rx/5Yjd8WABMEz7bAH/IupVXuV9qM1XQW1DXfG6/49Iv60VOEWeoJT3JJT3xSQEwWfh4AxAEQdzJLSzf02A5idfPRqWK2e245Au5tduEGN2zvhycTRvAtOETDkAUFt0r3lktGkcGu3IbD35v+eJ7GkfvqLO64a6cTwdy2U4HAB0LRQjmrqyiKndbKXcEFRpgW7frS/sF72ucfaed1wVKOOtD0YIApg9FCGatoaEpbXMuHaQZHexTt+0T+xfeJH37zE8UOIk4W4dxccAggDnAcYRgvuQKRdzGNKIrHRXRr3rD27aTF/IDQ2fFUlwOsWskrjUPYC6wRAhmSqvVxWxIou10k2cG125533rENFFw+II4SqZldg7R8PHJADAb+LiDOWIY5sSPlzgcZvKiIXU/fyzqPUQ8cvriK1SFgjkazhNiyyCAOUERgjk6/ksst4E/bmlI069f8V27WEcteDWBKpAy0eN4FthcAGBm8KEHs3M2Op5XIAx9K1BzbBOHL7CdufzVq3RuE3N6PM8SHwgA84MlQjAvV+NSmSRBn1d9uZf2UdIm27nvLrrK3G1mTo/nWfHZDgcAbMAPYDAjN2/myk9TXnPtJTcvqEryJK99Nf8qt0rJnMKyIIAZw6cfzMW9e2UV+xosJ/K8a27Kb161ff3r564KdTRzcjxPhL1jAMwYVo2CWaipq7+1tYgTpgviVcuvnhK98uWkBLEFjzgWgRYEMHcoQjB9coUicVMmFaAZ5qZuPreffmn1mHjbnhLO3lFcHC8IAPgaABOno3RnNyZS9rqxfcimkzua5nwx9KrTs77kpjCcOwYACALbCMG0MQwTvfkSh+BEDBM0H9t6Z/oXL6S4fRtKvtAVPwEB4A8oQjBl0TtiuQ384RMY2bGtJ8d89kmO19FwbpgLlgQB4E8oQjBZpw5d5hULBo5XqE/uWjPg8+Qm72tTuZ6WaEEA+B8oQjBN509fJTMEASMrNeeOvtzti+7O3vEhXJxEFAD+CkUIJijuYgoVz/UaXKBNuPCc75fvDveY6YuNggDwaChCMDWJV9JU5xmRf0ZTWsZ3wWt+j3DywupQAPh7KEIwKdeSMmQndTLHOEVxXfPz6w71sUYHAsA/QxGC6UhKzGj+XdNgedKSZxPyny+crfD2BoDHwzcFmAItTWw7kdb1qo4WHeo3IKTn5JlsJwIAo4EiBONWryZ25NEZ8Slz73GE4oOhs+aIegxkOxQAGBMUIRgliiEuljO779In79ErFXFz7lnau5wZuPA9rsSJ7WgAYGRQhGBMlDriUgVzvIT+vYT2seLM6cqZmnVcVOrYZXBez2dWESSOEwSAp4YiBEOnpogbtUx8JXOpnE6uZgY6ciZ5k0lTeE7SkjPbr3HVtgNftHHrNZ/tmABgrFCEYIgKpUxKNZNSw1yrYTLrmJ4SznBXzhu9yCPhpDWfoBWywlP7z2d4EpYWke8PtxBbsp0XAIwYihAMQqOGSKpirtXQKTVMSjUj5nEGOXEGO3G+CCYHOXEs//99SitlzbHHU5PT1Yrp3F7SaS9Gcjg4UBAA2gVFCKxp1hKx5XRsOXO5gimRMsFOnCHOnNd6kjtGkK4WDz9ZV1MmSzjZfCP2ut1EK+UMh2fEg0KHsJEaAEwNihA6W7mCOVLERJfQKTXMEGdOuAe5fTgZ5MDhPepsoLRSpryZqLh+QVd1vzowvIi/hNDRQSv9nZ0cOj04AJgmFCF0EjVFHCmmd+bRabXMZG9yeS8y3IMU/80bUFt5T52bprp9XVN8RxjQnxc2KfGmxiLDVjxUHTF5KFaHAoAeoQihw1UpiU051LY7dH8Hzis9yMnepOgvhznQSpm2tEBz/66m+LamKIcjEAoDgizDJtot+PDy5QzNb4TOUd13pa+rsyMb/wIAMGUoQuhAVUriy0xq9116VlfyyiSev+0fS3KMTquruqetvKetKNZVlmgriilZE9/dT+Dtb9FvmGT6Eq7EiaKohEs3mg7foPg635ku/YOwRRAAOgSKEDqEUkesu0VvyKLm+ZM5z/KdOArtvTzpjTxtaYGmvIiqr+Q5uvNdfXhuXcSDI/juvjwHN+L/V3jW1NWnHLzMyRSoLVRu4+1ChwWx+28BANOGIgT9iyllXr9Khdmr0wJv2pRkqONvVtRVCjy7Cbz8Rb1CrCNm8Vy8ONyH33sNDU1pydlNNxU2tRKVt7b7bPfevQNYyQ8AZgVFCPok0xLvJsjo7Kun6Ks26Vl87+7c7gPsZr0p8Oz2yPOf1dbW594urrnbQJSS1lKbRhe5XX+rwcMCrSzFnR8eAMwTihD0JisrL+n4iTcbk6y79bYOHiMKfJcU/dlnGq22ory6qrymoUqqqFYx9RxxoxWP4jXbN3E9CM8Jjn369xDw+SzmBwDzhCIEPdAU38k9vFNVW94zeLLfhJdlDC+/sKzmUoasWkU1UKSUJ5JbWOhEcoFcbalmJJTQkW/fy8ari7unhyvb2QHA3KEIoX0Ypun0r/evXjzgPKVvoD2/hI75IkekEzVbNelsdXx70ra7pZ2zjYuro6urE5eLq0MAgMFBEULb1ddUX939m6zR1Um7IripmRBT9v2tfbt7enm645h3ADAWKEJoi4bGpthfkiVljo1i1zw/1xXT/B0lNmyHAgBoCxQhPLXKyuobm3J1wvIbLvcUo176LBjvIgAwYvgKg6dTVlaZsTlfaZGWaeXmEvHsf3o/6lTZAADGA0UIT6Gsoipzc4HaIiXJ2i9wwpRF3dGCAGD0UITwpKpr6jI23VVapFyx7hYUNWlBAFoQAExBu4pQoVBcvHixqKjI1dV18uTJFhZ/XE31999/12q1Lbe9vb1DQkLaGxPY1tjUnLwxSyPKSrbqghYEAFPSriIMCwtzcnLq3bv30aNH33vvveTkZCcnJ4Ig5s+fP2zYMEtLS4Ighg0bhiI0dgqF8vL31zX84mxLScCEaWhBADAl7SrCmJgYV1dXgiAYhhk8ePCePXveeuutloc2b97cpUuX9ucD1lEUdWbjFR1RXy7S2Ua88EoPtCAAmJR2FWFLCxIEweFwrK2tH3zo4sWLtra2wcHBqENjF731IqNUqYT3G0e/9d8+aEEAMDX62Vnm0qVLmZmZe/bsabnbpUuX2NhYhUKxYMGC1atXL1++/JF/pVAoysvLv/jii9Yp48eP79Onz9/NRa1WCwQCvQQ2H+0ctLNHErhlXNLiRnLYe9/1odRqSo/ZDBbeaW2AQWsDDFobPO2g8fl8knzML3g9FGFGRsYLL7ywe/dud3f3limZmZktNxISEsaMGTNr1ixnZ+e//iFN0xRFNTY2ttzlcDg0Tbc/D+hL4pV0/k0LnvjYzoEf/RyCZUEAME3tLcLs7OyoqKjNmzdHRUX99dFhw4ZZWVnl5eU9sgitrKy8vLzWrVv3hPPSaDRCobBdcc1PmwctOydPe47DEx/d2Os/h8bZiMzpdNl4p7UBBq0NMGht0BGD1q4izMvLi4yMXLdu3YwZM1onMgzTesLlnJycpqYmX1/fdmWETldVXXdvd61OfPbHrq/vneRkVi0IAOamXUU4bdo0hmGio6Ojo6MJghg3btyiRYtOnjy5YcOGAQMGKBSK/fv3v/POOx4eHnpKC51BpVYn/Zipscg64DxlyzNd7PCDFQBMWruKcO3atUqlsvVu165dCYIYNWqUTqcrKCgQiUQxMTHBwcHtzQid69SPl3V042WJ70czB3hb4WpKAGDi2lWEkyZN+utEa2vrZ555pj0vCyw6degSr4YssK6fOvO5IAe0IACYPuwKCH9KuZbJSxU1W1yzm7I40hMtCABmAUUIfygpLWs4qpRans4b8+biQJyNHQDMBb7vgCAIQq5QZG69q7ZIie3/6uYhVmzHAQDoPChCIBiGObM5TkfWnu8WuXmiJ4l1ogBgTrBqFIhTey9yGzmJTo5rnu2LQwYBwNygCM1d0tU03k1xjs295fMiHUVspwEA6HQoQrNWXFLafFxTYpM4ZsGibjZYJQoA5gjbCM2XVCa/9VN+vTjNbdbSIS74SQQAZgpFaKZomon5/rKGV8NMfHFKNzHbcQAAWIPlADN1dMtJjpKpCRuyINiJ7SwAAGxCEZqjc0cvWNyzLext8+b47mxnAQBgGYrQ7KQmpTHJ1nd8a/71/Ai2swAAsA9FaF6KCgprjqly3XKWL56BA+cBAAjsLGNWamtrb/18r9y+4JXlL/PxEwgAgCAIFKH5kCsUl76/LhXXzFu5EKePAQBohSI0C2q19ve1Zxmu7rl3XhDzsUoUAOBPWEFm+pQa6sCaIzyKN3lllJWIz3YcAADDgiVCE9ekps99t89W7TzmX6G21pZsxwEAMDgoQlNW3KC58MNRB7XXiLf629nZsh0HAMAQoQhN1tUSee4vB5w0vkPf6uvg6Mh2HAAAA4UiNE0/p1QKTp+213UJWhLgiBYEAPh7KEJT06ghvojOGpCVLmI8R70ziIvD5gEA/hGK0KRcKGNO/RYdVsuQIsfxbw0Xiy2kUinboQAADBqK0ETUqYlVV5qCrm0ZJRukc1dNXRrB4+I/FwDg8fBdafQ0NPFjDh0fd21+bYxIHkWHqqZOj2I7FACA0UARGjENTfx6l96WUvVO+c6Zaje+arzrHNu+/ULYzgUAYExQhEapVM5sz2UOZDW93Xh0dd1tlXaKykkdtqKPra0129EAAIwMitCYlMmZk/eZw4V0dUXVJ5pTUyqu5wsj1ZrpdhOEYSNwcUEAgLZAERo0iiHympjUWia5momrYOTN0mW86xtqLysba4vE4bXqV4T+zaNmBInFFmwnBQAwVihCg6CiiGolU6YgKhRMqZwokTFFUqKgmbnbzAQKpFFkwRTNnZW1N7k1xbluQ3MUQ22VTnxf6YBX/B3sJGxnBwAwbijCTqWjifxmJruByW8miqRMmUxb16SUS2UcjcqLr/TmKTxImSdH2oNqctLUSOTVwsZSWqurcetVwfdN1Y6ypZxlzVLbfsKhY/tgKRAAQC9QhB1OrqGvZt8ruJ0nq67gyRoljNqGkfem1UGUmsMQBI/P4QtIHo/D4+tIMc0XUBwLJSEsovxoKpBPiq11tspKpcpZIekl7hfi5+TkwPY/CADApKAIO9DlK9fLL5baqhystLZ+XFsNT0DxPGm+jiBJNZejFjxw8jOKICiC4dMETRACgiskBA48C4nIyc3Ox9fD1gb7ggIAdBQUYUdJzbitOqlVu2u6hLv6B/gK+LgiLgCAIUIRdojKmvrSQ+V17jULV8xiOwsAAPwTFKH+6XT0pU1xGgstWhAAwPCRbAcwQb9+c4BHC59dOYXtIAAA8HhYItSzvbt+d23w7rbYy1IsYjsLAAA8HpYI9SkhMc0x240zggro6sN2FgAAeCIoQr2prWuqO9FY6V00IWok21kAAOBJYdWo3pzefI5jwcxfhh1kAACMCYpQP3758ZCd0nHovwawHQQAAJ4OilAPLsWnuBf58CdxHB1wCmwAACODbYTtVVff3HxGVtmlZPSIwWxnAQCAp4YlwvY69UMMR8SZt/Q5toMAAEBboAjbZc/OaDu5y6C3erMdBAAA2girRtvuRkaO021n3kiOqzMujQQAYKxQhG2kUmmKf7tX7lo0IXIE21kAAKDtUIRttG/DES2pnofTagMAGDlsI2yLY0cvuNf6eL/sxuXilwQAgHHD9/hTK7pXLrpm0dyvLrC7L9tZAACgvbBE+HQYhkn6OUUjkS+YM4ftLAAAoAcowqezc+NhCWUzedk4toMAAIB+oAifwsXLyR6lPpbT+dbWYrazAACAfmAb4ZOqqW1UnFVXdi0ZFoozawMAmA4sET6pmM2xjIVm/qs4XgIAwKSgCJ/Iri1HJGrbIW8PZDsIAADoGYrw8S5dvuZW6MmPIpwd7djOAgAAeoZthI9RXlUjj1FXdi0ZPSqE7SwAAKB/WCL8JzTNxG1OVFvJFryKowYBAEwTlgj/ya71h3g0f/qKZ9gOAgAAHQVLhH/rwN6TLrUe3nOdbXDUIACA6cIS4aPFJaTaZTrpwlS9ewewnQUAADoQivAR8vLvy04oK/2KpkwJZzsLAAB0LKwafVhdQ3PWjjuNdhULl8xjOwsAAHQ4FOH/UKu157+9rBHJF/z7RbazAABAZ0AR/ommmcNrTnA5xPSV00iSw3YcAADoDNhG+AeaZnavPiLQ8Me+MdrK0oLtOAAA0EmwREgQBEHTzJ7VR8Qqi7DXB+M8agAAZgVFSKjV2sNrTog0wiGvD/Zwc2I7DgAAdCpzL8Lq2vrY7xO5JDHqzWFYFgQAMENmXYRpGXdKD1aoxU3PrpxhKRaxHQcAAFhgvkV4cO9pSaZDrVvpS2/O5XCwjygAgJnqwCLUarV8Pr/jXr/NqmvrY7bE2suclaFNC6fjeEEAALPWIYdP7Nmzx9nZ2cHBYdSoUeXl5R0xi7ahaWbPjujUdXcYWtN3Wfdp08exnQgAAFim/yIsKytbsmTJiRMnGhoaevbsuWLFCr3Pog1omjmw91T0RxctCoTKQc3zV8329nRhOxQAALBP/6tG9+7dO3LkyJCQEIIg3n33XX9//8bGRolEovcZPaHkG7duX85xrnXnc7lyv8ZZ86fxeOa7ZRQAAB6i/0rIz88PDAxsue3j4yMQCIqLi/v37//IJ1MU1dDQ0HrXxsaGy+X+3Ss3S+V1DfJ/nrtUKquvra+prJNWNfMaSXuFK8mQHCuKGqKePmU8dooBAICH6L8IGxsb3d3dW+9aW1s/WHUPqqqqSk9P9/Pza52yZs2a2bNn/90rX/g5zkPq93ePttCSWg2p0vG0HAGtsaXogdSwYf243H4EQcjljylRkySXy1H/TwuD1gYYtDbAoLXB0w6aSCR67FpA/Reho6Njc3Nz693GxkZnZ+dHPtPFxSU4ODgpKekJX3n6W1HW1tZ6iGhOGIaxsrJiO4WRwaC1AQatDTBobdARg6b/nWV69uyZnp7ecvv27dscDsfHx0fvcwEAANAL/Rfh7Nmz09LSDh48WFVV9cEHH8yZMwc/eQAAwGDpvwgdHBx+//33b7/9dtCgQRKJ50PLpAAACC9JREFU5JtvvtH7LAAAAPSlQw6oHz16dHJy8r1793bs2GFjY6Ovl925c2d1dbW+Xs1MYNDa4JdffqmqqmI7hZHBoLUBBq0NOmLQjOnCvPv377979y7bKYwMBq0NMGhtgEFrgwMHDuTl5bGdwsh0xKAZUxECAADoHYoQAADMGooQAADMGodhGLbmHRMTM336dDc3tyd8fkVFhb29vVAo7NBUJgaD1gYYtDbAoLUBBq0NnnbQZs+e/dlnn/3zc9gsQoIg7t69+w8nF32IWq3GO+ZpYdDaAIPWBhi0NsCgtcHTDpqbm5uFhcU/P4flIgQAAGAXthECAIBZQxECAIBZQxECAIBZQxECAIBZ0//1CDvImTNnzp8/7+rqunjxYjs7O7bjGKiGhobU1NTCwsLBgwf379+/dXptbe327durqqqioqLCw8NZTGiA8vLyTp8+ff/+fW9v73nz5rW+uzQazY4dO/Ly8vr16/fiiy+SJH41/unGjRsXLlyorq52cnKaPXt266XWNBrN9u3b7969279//7lz52LQHun8+fNlZWULFixoucswzMGDB1NSUnx9fRctWvTYXRzNSkJCQk5OTuvdV155peWGRqP5+eef8/Pzg4KC5syZ0853mnG8Tbdt2/bKK6/4+fllZGSMGDFCq9WynchAzZgx49133/3000/PnDnTOlGlUoWFhWVnZ/v6+s6bN2/Pnj0sJjRAkyZNys3N9fb2jo+P79u3b21tbcv0uXPn7t+/39/ff+PGjStWrGA3pKGJjo5ubm728/MrKSnp3bv37du3W6bPmjXr0KFD/v7+33333cqVK9kNaZhycnLmzJnz5ptvtk756KOPPvvsM39//1OnTk2bNo3FbAZo//79e/bsKfx/rYc5zJo16/Dhw/7+/t9++60e3mmMwaMoys/P7/jx4y23AwMDf/vtN7ZDGSiKohiGmTZt2urVq1sn7t69u3///jRNMwxz+PDhwMBA1vIZJKVS2XKDpulevXrt2LGDYZjc3FwLC4uGhgaGYUpKSkQiUXV1NZspDVhkZOSXX37JMExOTo5YLG5qamIYprCw0MLCora2lu10hoWiqOHDh69fv97W1rZlilQqtbGxyczMZBhGqVTa2dldv36d1YyGZenSpZ999tlDE7Ozs8VicXNzM8MwBQUF7X+nGcESYWlpaVFRUUREBEEQJEmOHTs2Li6O7VAG6pHrB+Lj4yMiIjgcDkEQEREROTk5uDDTg0QiUettlUrVch3pK1euBAcHSyQSgiC8vb19fX2Tk5NZi2jAampq8vLy+vTpQxBEfHx8SEhIy5XXfH19PT09U1JS2A5oWNavXx8cHDxkyJDWKenp6RYWFn379iUIQiQSDR8+HN9vD0lNTV27du3BgwfVanXLlJZ3mrW1NUEQfn5+Hh4e7XynGUERVlRUWFtbt35bubi4lJeXsxvJuFRUVDg5ObXctrW1FQqFFRUV7EYyTN99951AIJgyZQrxv4NGEISzszPedQ/ZsmWLt7e3h4fH888/P3HiRAKD9jhFRUU7d+789NNPH5z40KDh++0hHh4ebm5ujY2N33zzzYABA5qbmwmCqKys1O87zQiKkM/n63S61rtarRYnJXoqPB6vdQAZhqEoSiAQsBvJAB09enTdunW//fZby7uLx+NRFNX6qFarxaA9ZMGCBampqadPn96zZ8/BgwcJDNo/oml6wYIFa9asaVnl0Oqvg4bvtwe9//77P/zww+rVq5OSksRi8ZYtW4gOeKcZQRG6u7srFIqGhoaWu2VlZU9+nm4gCMLDw6P151JlZaVOp8MAPuT48eOvv/76mTNnAgMDW6Z4eHiUlZW1PqGsrMzd3Z2ldAZKJBI5OTmFh4cvXry4pQgxaP/g/v37SUlJn3zySXBw8EsvvSSTyYKDg/Pz8z08PCoqKmiabnkavt/+DpfLDQ0NLSwsJDrindaeDYydZtCgQVu3bmUYRiaTeXh4xMXFsZ3IoD20s0xMTEyXLl0UCgXDMBs2bBgxYgR70QxRTEyMi4tLSkrKgxOrqqosLCzy8/MZhklJSZFIJHK5nKWAhujB0ZgyZcqKFSsYhikvLxeJRC279iUmJjo4OLTuiARqtfrG/9u5c6eVldWNGzcUCoVWq3V3d4+JiWEYpqKiQiwWtwwgtGh9C8nl8p49e27YsIH5/3daUVERwzBXr15t/zvNOIrw/PnzDg4O8+bN69u37zPPPMN2HMP1xRdfDBw4UCKReHh4DBw4MDo6mmEYivq/9u5eRWEgCsOwvaWKCYqDkEYG/MEmjdjYiN6BpNArsLfwJgStvRIRwUKdXu3UQpjKQkizukXYILvl7pLIvE+XkyKHYcIHwyH5aLfbtVrN87xUKrVcLqNuM16SyWQ2m61/mU6nQX08HhcKhcFgYFnWZDKJtsm4sW270+l4niellFJer9egPhqNhBDBos1ms2ibjK31eh1OjT6fz/l8nslk+v2+4zjD4TDCxmIol8t1u91er5fP51utlu/7Qf1vd9rb/H3icrmsVivLsprNZjAAiZ9Op5PWOrwUQqTT6UQi8Xg8FouF1rrRaHBa9Y1S6vUtsG07XCKl1H6/L5fLUsqIuosprfVms7ndbkII13Vfx5W32+3xeKxUKuE5M7653+/BNwfCyuFw2O12xWLRdd0IG4uh8/mslPJ933Gcer3+eivYadVqtVQq/fIpbxOEAAD8hzcYlgEA4P8QhAAAoxGEAACjEYQAAKMRhAAAoxGEAACjEYQAAKMRhAAAoxGEAACjEYQAAKMRhAAAo30CieRsRpLUqW0AAAAASUVORK5CYII=",
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Two versions of the leapfrog method is used to solve a system consisting of\n",
"# a spring attached to a \"wall\" moving at velocity v. A mass m is attached to \n",
"# the other end of the spring. The mass experiences friction as it moves (it\n",
"# sits on a surface) and there is a rest friction as well that has to be overcome\n",
"# in order for the mass to move at all.\n",
"\n",
"# The function accel() adds up all forces and returns the acceleration\n",
"# a rest friction force fr0 has to be overcome before the mass can move\n",
"# a very small non-zero v is treated as v=0 and then the rest friction is applied\n",
"# fr1 is a coefficient regulating the v-dependent dissipation\n",
"# The spring has native (relaxed) length lspr and can be elongated without bounds\n",
"# and compressed to length 0. The final compression force is divergent to avoid the\n",
"# spring length becoming negative (in principle a max length should also be imposed)\n",
"\n",
" function accel(x,v,t,mass,kspr,lspr,xspr,kwall,fr0,fr1)\n",
" fspr=(xspr-x-lspr)*kspr-kwall/abs(xspr-x)^6 # force from spring\n",
" if (abs(v)<0.000001 && abs(fspr)