{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Series Solutions With Python\n", "\n", "We are going to see how we can use Python to graph solutions (or approximations to solutions) of ODEs. The \n", "example we are using is: \n", "$$\n", "(1-t^2)x''(t) - 2tx'(t) + \\alpha(\\alpha + 1) = 0.\n", "$$\n", "\n", "\n", "Solutions can be expanded in a power series centered at zero on the interval $(-1,1)$. We can write this as: \n", "$$\n", "x(t) \n", "= a_0\\sum_{n=0}^{\\infty}a_{2n}t^{2n} + a_1\\sum_{n=0}^{\\infty}a_{2n+1}t^{2n+1}\n", "=: a_0x_{E}(t) + a_1x_{O}(t).\n", "$$\n", "Recall that the relation that the coefficients satisfies is: \n", "$$\n", "a_{k+2}\n", "= \\frac{k(k+1) - \\alpha(\\alpha + 1)}{(k+2)(k+1)}a_{k}.\n", "$$\n", "\n", "Writting this interms of $k$ instead of $k+2$ gives: \n", "$$\n", "a_{k} \n", "= \\frac{(k-2)(k-1) - \\alpha(\\alpha+1)}{k(k-1)}a_{k-2}\n", "$$\n", "\n", "This is not a relation that is easy to solve for and we will use python to compute the coefficients. The first thing we \n", "will do is to write a function that computes the coefficients in terms of $\\alpha, x(0), x'(0)$. As you can see, \n", "this is kind of diffucult to reat, but it's doable. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "zoo*a*x0*(a + 1)*(-a*(a + 1)/8 + 3/4)*(-a*(a + 1)/24 + 5/6)*(-a*(a + 1)/48 + 7/8)*(-a*(a + 1)/80 + 9/10)*(-a*(a + 1)/120 + 11/12)*(-a*(a + 1)/168 + 13/14)*(-a*(a + 1)/224 + 15/16)*(-a*(a + 1)/288 + 17/18)*(-a*(a + 1)/360 + 19/20)\n" ] } ], "source": [ "import sympy as sp \n", "\n", "def A(k):\n", " a = sp.symbols('a')\n", " x0=sp.symbols('x0')\n", " x0p=sp.symbols('x0p')\n", " if k==0:\n", " return x0\n", " if k==1:\n", " return x0p\n", " return (((k-2)*(k-1) - a*(a+1))/((k-2)*(k)))*A(k-2)\n", "\n", "print(A(20))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, maybe the symbolic function has some uses, but it seems to be complicated here. Actually, it's really just the \n", "fact that $\\alpha$ (```a``` in the code) is symbolic that is causing problems. In the code below, we keep ```x0``` and \n", "```x0p``` as symbolic variables, but define ```a``` as a literal number. \n", "\n", "In this case, we can also get python to write out a polynomial. Of course, the code below doesn't print the polynomial in\n", "the nicest format -- and you can always do \"easy\" things to make this polynomial look better. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(x0p)t^1+(-1.66666666666667*x0p)t^3+(0)t^5+(0)t^7+(0)t^9+(0)t^11+(0)t^13+(0)t^15+(0)t^17+(0)t^19+\n" ] } ], "source": [ "def A(k):\n", " a = 3\n", " x0=sp.symbols('x0')\n", " x0p=sp.symbols('x0p')\n", " if k==0:\n", " return x0\n", " if k==1:\n", " return x0p\n", " \n", " return (((k-2)*(k-1) - a*(a+1))/((k-1)*(k)))*A(k-2)\n", "\n", "poly = \"\"\n", "for k in range(0,10):\n", " poly = poly + f\"({A(2*k+1)})t^{2*k+1}+\"\n", "print(poly)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is good, but if we want to have Python plot stuff, we need to have literal values for ```x0``` and ```x0p```. We do \n", "that below: " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1)t^1+(-1.6666666666666667)t^3+(-0.0)t^5+(-0.0)t^7+(-0.0)t^9+(-0.0)t^11+(-0.0)t^13+(-0.0)t^15+(-0.0)t^17+(-0.0)t^19+\n" ] } ], "source": [ "def A(k):\n", " a = 3\n", " x0=1\n", " x0p=1\n", " if k==0:\n", " return x0\n", " if k==1:\n", " return x0p\n", " \n", " return (((k-2)*(k-1) - a*(a+1))/((k-1)*(k)))*A(k-2)\n", "\n", "#And odd polynomial\n", "poly = \"\"\n", "for k in range(0,10):\n", " poly = poly + f\"({A(2*k+1)})t^{2*k+1}+\"\n", "print(poly)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As we saw in class, if $\\alpha=2n$ is even, $x_{E}(t)$ will be an $2n$ degree polynomial and if $\\alpha=2n+1$ is odd, then \n", "$x_{O}(t)$ will be a polynomial of degree $2n+1$. \n", "\n", "We need to fix the function a little and make it so that we can pass in the arguments for ```a```, ```x0```, ```x0p```. \n", "This is just basic good programming practice. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1)t^1+(-1.6666666666666667)t^3+(-0.0)t^5+(-0.0)t^7+(-0.0)t^9+(-0.0)t^11+(-0.0)t^13+(-0.0)t^15+(-0.0)t^17+(-0.0)t^19+\n" ] } ], "source": [ "def A(k, a, x0,x0p):\n", " if k==0:\n", " return x0\n", " if k==1:\n", " return x0p\n", " \n", " return (((k-2)*(k-1) - a*(a+1))/((k-1)*(k)))*A(k-2, a, x0, x0p)\n", "\n", "poly = \"\"\n", "for k in range(0,10):\n", " poly = poly + f\"({A(2*k+1, 3, 1, 1)})t^{2*k+1}+\"\n", "print(poly)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, we plot the solution for the parameters $a=3$, $x(0)=1$ and $x'(0)=1$. We plot degree $3,4,5,6$ polynomials. Each of\n", "these is an approximation to the actual solution. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3yN1x/A8c/J3okMsULE3iGx96y91aqqamlVKYru0v50lxZdlJoVqvbeqyhB7JkgQ0IEGbIk9/z+eC4NEjJucjPO+/W6L7n3Oc95vvfG/d6Tc88QUkoURVGUgsvE2AEoiqIoOaMSuaIoSgGnErmiKEoBpxK5oihKAacSuaIoSgGnErmiKEoBpxK5kmeEEL8KIT7Oo2s1F0JczItrGZsQwlMIIYUQZtk8f7AQYpuh41LyjlDjyAsPIcQeoA5QQkqZZORw8pQQQgKVpJRXjBjDHmCJlPL3bJ7fDPgGqAGkAueBd6SUR59znidwFTCXUqYYqqxScKgWeSGhf4M2ByTQPRevk61Wn/JsQggHYAMwC3AGSgNTgSL1gaxkj0rkhcfLwGFgATA07QEhxAJ9t8Z2IUSsEGKvEKJcmuNSCDFGCBEkhLgthPhWCGGiP/aKEOIfIcQMIcQdYIoQwlEIsUgIESmEuC6E+EgIYSKEcBZChAohuunPtRNCXBFCvJwmjv/pf26lLztJCHFLCBEuhOgphOgshLgkhLgjhPggTYwNhBCHhBD39GVnCyEs9Mf26YudFELECSH6P6w/zfnVhBB79OefFUJ0T3NsgRDiJyHERv3r868QooL+mNA/91tCiGghxCkhRM0nX3whxDS0D9LZ+hhm6x9vIoQ4qj/3qBCiSQa/v8oAUsplUspUKWWClHKblPKUvh4T/et8XR/LIiGEY3oVCSGuCSHapbk/RQixRH/34Wt1Tx9nY/3v+ECa8hnGrH8NP9f/n4gVQmwTQrjqj1kJIZYIIaL0r/NRIYR7Bs9XMSQppboVghtwBRgF+AAPAPc0xxYAsUALwBL4ETiQ5rgEdqO1BMsCl4DX9MdeAVKAtwEzwBpYBKwF7AFPffnh+vIdgAigODAXWPlEHP/T/9xKX+8ngDnwOhAJ/KmvtwaQCHjpy/sAjfQxePJft0Pa51Axzf1WQKj+Z3P96/MBYAG00b8eVdLEdQdooK9/KeCnP/YCcAxwAgRQDSiZwe9gz8PXTX/fGbgLDNHXO1B/3yWdcx2AKGAh0Ako9sTxV/XPwQuwA1YBi/XHPPXP30x//xrQLs25U9C6fJ4qm+Z3fCAzMeufYyDaB4+1/v5X+mMjgfWADWCq/505GPu9URRuqkVeCOj7VssBK6SUx9DeaIOeKLZRSrlPan3nHwKNhRAeaY5/LaW8I6UMBn5AewM/dENKOUtqfarJQH/gfSllrJTyGvA92hsfKeU24C9gJ9AF7c2dkQfANCnlA8APcAV+1Nd7FjgL1NbXe0xKeVhKmaK/5m9Ay0y+RI3Qkt9XUspkKeUutG6MtM9xlZTyiP45LgW808RoD1RF+07pvJQyPJPX7QJcllIu1se9DLgAdHuyoJQyBmiGlmTnApFCiHVpWrSDgelSyiApZRzwPjAgF7q6MhPzH1LKS1LKBGAFj79WLmgfqKn631mMgeNT0qESeeEwFNgmpbytv/8nT3SvACEPf9AngjtAqfSOA9efccwVrVV7/YnypdPcnwPURHvDRz0j7igpZar+5wT9vzfTHE9AS8AIISoLITYIISKEEDHAF/pYMqMUECKl1D0j5og0P8c/vK4+6c8GfgJuCiHm6PuzM3vd60889uR1H9F/SLwipSyD9vqVQvtQTa+u62gtZkN3XWQm5nRfK2AxsBXwE0LcEEJ8I4QwN3B8SjpUIi/ghBDWwItAS32SiwDGAXWEEHXSFPVIc44d2p/QN9I7jta9kvZY2qFNt9FaXuWeKB+mr9sUrbW8CHhTCFExm0/tSb+gtQwrSSkd0LpJRCbPvQF4POz3fzLm55FSzpRS+qB191QGJmZUNJ3rlnvisUxdV0p5Aa3L52F//JN1lUXrmrrJ0+6jdW88VOIZMT4pJzE/kFJOlVJWB5oAXdG+u1FymUrkBV9PtKFq1dH+xPVG68fdz+Nvos5CiGb6Lwg/B/6VUqZtaU8UQhTTd7eMBZandzF9C3oFME0IYS+0L03HAw+/THv4BeWrwHfAIn1yzyl7IAaIE0JUBd584vhNtP7j9PyLltwmCSHMhRCt0LoK/J53USFEfSFEQ33L8j5av31qBsWfjGETUFkIMUgIYSaE6I/2e9qQznWqCiEmCCHK6O97oHX9HNYXWQaME0KU138QfwEsl+kPIQxA63YxF0L4An3THIsEdGT8WmU65nSeQ2shRC397zsG7QM/o9dKMSCVyAu+oWhdGMFSyoiHN7TugMFp+lD/BD5F61LxQetzTWst2pd6AcBGYN4zrvk2WlILAg7o654vhPBBS+ov6xP+12gtwPdy/jR5F63fPxatD/nJD5opwEL9aIkX0x6QUiajDcnshPYXxc/6GC9k4roO+uvdRetiiEL7gErPj0BfIcRdIcRMfbdSV2CC/rxJQNc0XWBpxQINgX+FEPfREvgZ/bkA89G6LvahjQNPRPs9pOdjoII+5qlovx8ApJTxwDTgH/1r1SjtiVmM+UklgJVoSfw8sJf/PuCVXKQmBBUBQogFaCM4PsrguNEn0yiKkn2qRa4oilLAqUSuKIpSwKmuFUVRlAJOtcgVRVEKOKMsgOTq6io9PT2NcWlFUZQC69ixY7ellG5PPm6URO7p6Ym/v78xLq0oilJgCSGenHULqK4VRVGUAk8lckVRlAJOJXJFUZQCTiVyRVGUAk4lckVRlAJOJXJFUZQCziCJXAgxX7+P4BlD1KcoiqJknqHGkS9AWzZ1kYHqU5R8KSUllds3wwi/doY7YYHE3Qom6c4tUqLvQkI8Ju6lcKnZjBotelPS1QUhMrv3haJkn0ESuZRynxDC0xB1KYoxpKakcGjjfML2b0DGxsD9BMwSkjFPTMEiMRWrRIlVItgkgkUqWPH4Xnj/CQOOckfM4JIT3CtmRpyrPSklPbCuUp/yDXtQrYIn9lZqBzTFcPJsZqcQYgQwAqBs2bJ5dVlFydDdWyHsXfgFD/yPUPJqPC4x2s7BAMmmEG8NCVaQZCWIKWZGlLUZqdYWSFsbTOzsMXdyxtrFHQd3D1xKV8LZrTSX/bdx68y/PAgOweJ2DE5RKVS4dhcz3V3gFDCPU/YQVcyUaGdb4t1LI7y8ca/bkRY+3hSztTDiK6IUVAZb/VDfIt8gpaz5nKL4+vpKNUVfMYZT+1ZzfvVcbM4HUzY0FYsUSDSH4LKmJFbzokKnQVSo3Zxibunuj5wtydE3CTq8husnDnD/6lVMb0VjG5WC612wfPBfuZDigvN1qlN+0Md0rF8LCzM1FkF5nBDimJTS96nHVSJXCrPEuGj2LP6SmIO7cQuMocQd7fFIJ4jwssOyQRNavfwhjs7F8zy21NhbhAds4fqJvUSdv4DdmTuUjIQkMzhV0Y7QZv3o8uJr1PYopvraFUAlcqUIuXp6P8f+/AGzM5fwuJ6CTTI8MIXg0ibEVi6FZ5ehNOo4ON8lR11CDGeWf8K1jdspc1GHdTKEOwsOV6uCQ7ex9GndkJKO1sYOUzGiXE3kQohlQCvAFW0n8U+llBlu3qsSuZIbdi2YQvyyv6hwXQfAPTsI9bQCb2+aDPmA0uUqGTnCTJKSuFNrObzoG+TxO5QJF6SYwInydpyq25HGPYfTqXYZbCyMsnipYkS53iLPCpXIFUPR6XTsmfcRScvX4hmqI8YGAr1dcG3bi1b93sbCooB/eXjnKpdXfcbFHQdxv2iCXQJEOgh2V/YivsUwerZtRcPyzpiY5K+/LpTcoRK5UqjodDp2/zqJlL82UTZccs8OAhuUptV7v1GmbAVjh2d4SXEkHV3I4VVzSDiTSLlgE3QCTpa1ZX+VJtTsPpLXWlXF0szU2JEquUglcqVQSE1NYc/s8ci/d1D6liTKAYIaeNDqvbmULVPO2OHlPinh6j6CN3/P6X/O4nzJAqc4iHAyYVmTLgwYMYE2Vd2NHaWSS1QiVwq0lOQk9sx+B5NVeyl5WxLpBFcbetJy4lw8y5QxdnjGcfcaDw79ypGtK0k5akbxO4JDXo4c7fgOkwd1w9PV1tgRKgamErlSICUn3mfvD2OwWHuQ4nchwgWuNqhA83G/Uams4cZ6F2iJ0SStHceujXsoedwSKWCFdzUcX/yYMS/UUl+KFiIqkSsFSlJ8LHumj8JmvT+u0XDDDYIaVKbZmF+pVq6kscPLn06vJOSvdzntb0X5IEGIsxmLG/Wg39BRdKtdMt8Nt1SyTiVypUBITU1h51fDsVtzhGKxEOoOgfWr0fStn6lVvoSxw8v/okNh9Rsc8g8g9agDLtGSvZWcOdR6LO8P6kjVEg7GjlDJAZXIlXzvypHNBH0wEY/QVIJLwWWfmjQa8SO+ldJfnkrJgE4Hh2Zzf9tn7Alyw+O4NiFqab062HefzPiONXG0UYt2FUQqkSv5VnJCPDs+GUyZTRdIMofjjYpT7e3FNK/uoboDciLiNPz9OkHXg7h4piSeQckEuZnzR4M+9O0/jBd9PdT48wJGJXIlXzq9y4/IKZ9T8paOS5XgbOvhvDlyjFoF0FAeJMCOqegO/8KBWE/M96fgFKtjW9XiHGg6hg/6t6Fu2WLGjlLJJJXIlXwlLjqKPe+/RPnd14i2hVMNimPT/xeGtKimWuG54cpOWDOK2Ji77A2vhueBW9y3Eiys64N9x3FM6VkHaws1mSi/yyiRq3UylTx3ZPUvBHRoToVd17hcXbKyx0s0/WAdL7esrpJ4bqnYFt48iH319nQtHYD5MHfuuVvz9kF/Gv4ynKHfLyXsXoKxo1SySbXIlTxzJ+I6Bye9QoUjEUQWk1yq786lBtP4uG9DtWNOXpESApbC5snohCl7ZAuc/I6RZCaY2aoP494eR4PyzsaOUsmAapErRqPT6dj3xzQud+lEOf8ILtXT8Web/tgPWczXg5uqJJ6XhIC6L8Eb+zFxq0yb5PU4vFGTZGtTPtiyEr+po1h88Kqxo1SySLXIlVx1I+g0ARNHUv7sXW64S274FGdN6Yl8+VIrKrvbGzu8oi01BfZ9A3u/JrpcBw5sDMfrTBS7K7sQOvhbpvRpqHYpymcyapGrubtKrkhNTWH3rPdwnr+RUjoIbJTKphJ9cGwwiCVdq6sv1vIDUzNo/QHYuOC4eRIde3dmSxUvWv99lMCf3mBE8AS+e7M/rnaWxo5UeQ7VIlcMLuxKAGfGvEbZoPtc99CRWM+VH6zGMKF3c7rVUZN78qV/f4PNk6BqVw7ii+WXc0gxEfzcuicTJ0ymZmlHY0eooFrkSh459NcszL74BbdUSUiLB2x16UVgqe4sGuRDWRcbY4enZKThSO3fzZNoUhWCFv7C9bfHMGnzahbdvkTjd2bQs56HcWNUMqQ6wBSDSE6KZ+P4F3H6+Gdi7XTEd7HkU+dPcGk6jJVvNlVJvCBoOBI6fQMXNuB1YS6+qzcRVLs4rxw5y61P+vPl6sOk6vL+L3jl+VSLXMmxsCsBnHvrNbyu3yewZgrWNcszVjeKrwc0pl11tclBgZKmZW4PdFm8lc1fvE1TvwOEfD+Ct6+O58s3B+NorUYa5ScGaZELIToKIS4KIa4IId4zRJ1KwXBwxUzC+w7CJfw+ka3vc7dOM8bJccwZ3kIl8YIqTcvcZNVwunz0E3FfjsUlLoVhC79m4kcfceVWnLGjVNLIcSIXQpgCPwGdgOrAQCFE9ZzWq+RvyUnxbBzXj2Kf/EK0g6RY53vsKfUS38uBLH29CT7l1KSSAi1NMmflMBp2e5XiyxYQU8ySsRvXsWjCEHacDTd2lIqeIVrkDYArUsogKWUy4Af0MEC9Sj4VevkEe7u3wGvzGQJrSRp2iGeazXusE23wG9GIWmXUCIdC4YlkXq6SNw3X7uSKT0kGHr1A8MRezNx4BGOMfFMeZ4hEXhoISXM/VP/YY4QQI4QQ/kII/8jISANcVjGGg8t/JKLfYK0rpU0C7ZpYMyDlM86Y1WT5yEZq44LCpuFI6PTto2RuZ2tPt8U7uPpyaxpejabqtOF8uXiNSuZGZohEnt4qR0/9VqWUc6SUvlJKXzc3NwNcVslLyQnxbBjXl2Kf/kp0MROcO0XhXacabe59RLSNB3+90ZgKbnbGDlPJDQ1HPJbMTXQpdP7gZ+K/nYhzfArNZ33E9wv9jB1lkWaIRB4KpB1gWga4YYB6lXwi9PIJ9vZoToXNZwn0taFZsxCsavWkaehbWNm7sGJkYzyc1fDCQu2JZE5KMvW7vortL9MxT5U0m/k5P877w9hRFlmGSORHgUpCiPJCCAtgALDOAPUq+cA/fjOI6DcY55vx3O5mQ9eKgQT7vEeLCz0p7eLA8pGNKelobewwlbyQTjKv1rgTrgt+44GZoOnMb/jpt5+MHWWRlONELqVMAUYDW4HzwAop5dmc1qsYV3JCPBve6YvzlDlEu1jg2huaO4VzrPFPdD5ah6olHFj2eiPc7NU6HEVKOsm8Yu3mlF2ymFgbU5rMns2vs78zdpRFjkHGkUspN0kpK0spK0gppxmiTsV4bl4/z57eLamw5SyBzUvRrOUtytkLdjVdSL89TtQu48SS1xqq7diKqrTJ/K9XIPUB5SrXo4rf30Q6mdPkl3n8Pn2KsaMsUtQUfeUxATv9uNqnL8VD44gcWpuupY9hWaIyq3wWM3xrMo28XFj0agMc1BriRdvDZH5xI+ycCkCpclWovXwdwe5WNJq7nD++nGjkIIsOlcgVQNv8Yfv08Zi+PZVkSxNs3mlCi6QtUK0bC6v8zPjNEbSq7Mb8V+pja6lWdlDQkrnvcDg4Cy5sBKBEKU8a+G3iYllbGi3cwMJPRxk5yKJBJXKF+Lh7bBreiTJzNhNS3YXaE9pTOWQl+L7KT24f8+nmq3SsUYLfhvhiZa7WEVfS6PgllPSG1W/CHW1noeLFS9LKbxsnKjjSYPlulkx82chBFn4qkRdxwReOcqhHa8ofCuZqv4a88GpznE7NQ/q+xndmI/l222V6epdi9qC6arcY5WlmlvDiQm02yV+vwINEAFyKOdPxz60cquaKz/qjLBvdF51OZ9RQCzP1zizCjqydw82BQ3GMSiL2f6Pp3LQUpkd+hvqv85vdm8zeE8iA+h58/6I3Zqbqv4qSgWKe0PNXCA+ArR88etjZ0ZGei7axq3YpvHecZeXrXUlNTTFenIWYencWQTqdjs2fvY7t5BnEOVrgumQejexC4fBP0GAkW8pO4OutF+lWpxRf9KqFqUl6k3cVJY2qnaHJGPCfB6f+evSws701L/6+gU0+5an1z1XWvvwCyUnxRgy0cFKJvIiJuRPBlkGt8fzzAFcblKb+mh14hqyBwz9Dwzc5XesD3lkRgLeHE9/2rY2JSuJKZrX9BMo2hvVjIfLio4fdHKwZ8usqVjauTrVjN9g8qAMJ8TFGDLTwUYm8CLlyYg8nenSg7MlbBA9rR+c/tmJ38Gv491doNIrwxp8wfJE/LraWzFFfbCpZZWoOfeeDuTWseBmS7z86VNzeitd+WMKiFj5UPhvFrhfbExd924jBFi4qkRcRB5Z+R8zQN7GKTyFpxge8MGkmJlsmw5E50Hg091t9xvCFx4hPTmX+K/XVjE0lexxKQZ/ftRb5hvGQZlXEko7WjPpqDr+1bopnYAz/9HuBu5Ehz6hMySyVyAu5lAfJbJw4EJfP5xFVwgaPv5bj02EwbJwAR+dC49Gktvucd1ac5EJEDLMH1aVKCXtjh60UZBVaQ6v34JQfHF/42CEPZxve+fxHZrZtR8mweI4O6kHMnQgjBVp4qEReiEWFX2VH35Z4rQ8gsFVFWqzdSynPGrBpgvalVJMx0OF/fL31ItvP3eTTbjVoVaW4scNWCoMWE8GrNWyaBOEnHzvk6WrLxE++5IfWHSgZlsChwd2Ii44yUqCFg0rkhVTATj8u9uhGySv3CB/Tm66/rsfS0gY2jgP/+dD0HWj/GX5HQ5izL4ihjcsxtImnscNWCgsTU62LxcYFVgyFxOjHDlcsbsdr73zM981aUOZaHAde6kJ83D0jBVvwqUReyKSmprDlfyMxGz2VVDOB+dxvaTNqGuh0sOEdOLYAmo2HdlP4JzCKj9acoWVlNz7uqrZZVQzM1hX6/QH3gmHtW4/1lwM0qehKsyGTmdG0CR6Xo9n7cmeSEtSmztmhEnkhcifiOtv6t6Lckn1cr1uC2uu3Ub1JVy2Jrx+j9Vc2nwBtPyHw9n3eXHIMLzdbZg2qqyb8KLmjbCNoPxXOr4fDvzx1eHiz8ti0Hc1PTerjee4uO4d2VuPMs0G9ewuJkztXcKF7F0qfjyJ0RGc6LdmJo0tJLYmvextOLNb6Ldt8zN34B7y64CgWZibMG1pfrWSo5K7Go6FqV9j+MYQceeyQEIKv+9Tmap3XmdekLuVPRbJ9WBceJCcaKdiCSSXyAi41NYUt097AdPSnpJoJzOZ+S/vx32NiYgK6VFg3GgKWQMvJ0PpDklJ1jFxyjPDoRH4b4qu2aFNynxDQ4ydwLKOtx3L/8S82rS1MmTPEhx2er7C8eR28jkewdUQ3Uh4kGyfeAkgl8gLsUVfK4r2Pd6WA1h+5cQIELIWW70HrD5DAB6vOcOTqHb7rVwefcsWMGr9ShFg7Qb+FcP82rHpd+0sxDQ9nG2YNrMsi1yFsaVmLCodD2fJWL7U2SyapRF5Andy1gvM9tK6UkNc7/teV8tC+7+DYH9rolNbvA/DznkD+Ph7KO+0q0b1OKSNFrhRZpbyh01cQuBP2f//U4eaV3JjUsSo/FhvKwTY1qLAviM1j1aqJmaESeQGj0+nY+uUoTEZ/ijQRmP72FR0mzNC6Uh46sQR2/w9qD4B2UwDYfDqcb7depId3Kca2rWSU2BUFn2FQ60XY8wUE7X3q8MgWXnSpXZJpDkM52bYaFXZcZPOE/iqZP0eOErkQop8Q4qwQQieE8DVUUEr67t4KZkv/lpRduJvgOu7UXL+FGs16PF7o8nZYN0abjNF9FgjBqdB7jFsRgE+5YnzdpzZCqIWwFCMRArrOAJdK8PdwiAl/4rDg2761qVzCkS+Lj+BCq0p4bT7D5g/V5hTPktMW+RmgN7DPALEoz3Bqz0rOdeuMx9nbWlfK0l04uZZ+vFDYMW2xIvca0H8xmFlw414Cwxf642pnyW9DfNRCWIrxWdpp/z+T4mDzpKcO21iY8dsQH6QU/FblHa40L4/X6mNsnjLcCMEWDDlK5FLK81LKi88vqWTXw64UMepjpInAZE46XSkAd4Jg6YvaJIzBK8HSnvtJKQxf6E+ifiEsVzu1EJaST7hVgRbvwvl1cHnHU4fLudgyc2BdLt5KYHvLqQQ28sDT7yBbv1R7gKZH9ZHnY8EXjrK1TzPKLtxNSO0MulIA4iJhcW+QqfDSKrB3R0rJ+6tOczEihtmD61HZXS2EpeQzTd4Gl4qw6d1HW8Sl1apKcd7tUIV1Z25xY+hMguqVoOzC3WyfPt4IweZvz03kQogdQogz6dzSySjPrGeEEMJfCOEfGRmZ/YiLgOSkeDZ/9jp3+r6Me+BdQkd0ouOf6XSlgLbm858vQmw4DFoBrtoXmWsCwlh38gbj2lWmZWW3PH4GipIJZpbQ+Tu4exX++SHdIqNaVaBTzRJ8tS0Qh48Wc7W2G2XmbGbn7PfzONj8Tcgn1j/IViVC7AHelVL6Z6a8r6+v9PfPVNEiJ2CnH/c++xL3m8lc9XbH+4tZlPKqlX7h1BTwGwhXdkD/pdp2W0BwVDydZ+6nWkl7/EY0Vlu1Kfnbylfh/AZ46zA4ez11OC4phV4//cPtuCRWvubNpdE9KXvuLrcmDKD1658aIWDjEUIck1I+NbCkwHettGrViq5du+b5dRcsWICDgwNWVlZYWVnRsGHDHNV391YwG0Z0wfKtqZgnpnLv81F09tuTcRKXEjaMhcvboMv3j5J4SqqOd5afQAiY0d9bJXEl/+swDUwtYNPEpxbWArCzNGPOy76k6CRjVl6k0bz1hFZywG26H//8Od0IAec/OR1+2EsIEQo0BjYKIbYaJqz8LyIiAktLS8LDw0lMTGT9+vXZqken07Fn7hQudeqI5/4ggjrXos62vTTu9/azT9zzpTZevMUk8H310cOzdl3hePA9pvWqRZliavq9UgA4lITWH2h/WZ5P/31U3tWWHwd4cy48hs+2hdBk8QbCy9piP20uR9bOyeOA85+cjlpZLaUsI6W0lFK6SylfMFRgz9K+fXssLCxwdnbm6tWrjx7ftWsXrq6u2NjY4OjoyKZNmx49bmdnh62tLS1atDDIOOp9+/ZRq1YtihXTprkXL571DRmCTh9gW48muH+/nBhXG8wWzKDL9BXYObo8+0T/P2Dv1+D9kvYGePjwtTvM2nWZ3vVKq5mbSsHSYAS414Qt72vDEtPRpqo749pVZvWJMFaejcN3yWpuF7fE4sMZBGxflscB5y8FrmtlyZIl7N+/nxs3bnDq1Clu3Ljx6Fjfvn1ZtmwZ8fHxTJ8+ncGDBwMwcOBAhg0bxv379ylTpkyGdTs6OmJtbf3U7ZtvvnmqbEREBEeOHMHW1hYXFxcWL16c6eeQEB/Dpg+HEjfgddyCYwh7oyttNh2masOOzz/54mbYOB4qtoduP2gTLICYxAe8szyAMsVsmNq9RqZjUZR8wdRM6yKMCYV9T7/fHhrduiLtq7vzv43nuRhrQ83FK4hxMif13c84d3BDHgacz0gp8/zm4+Mjs6tnz56yefPmj+7Xq1dPdunSRYaHh0tAWllZPbpZWFhIKaUUQsiEhAQppZQhISFSe9o50759e9m9e3eZnJwsP//8c2lpaSnv3r373POObpgv9zatJc9VqSrXD2krbwZfyPxFg49I+bm7lL+1lDIx9rFDY5cdl17vb5T+1+5k7YkoSn6yZpSUU52lvHkuwyIxCcmy9be7ZaMvdsjohGQZduWkPNCghjziXU1ePr4rD4PNe4C/TCenFqgWeQokUWcAACAASURBVHJSPFLKdLtGUlJSMDExISEh4dEtKSkpS/VnpUXu5eVFr169MDc35/333yclJYXLly9nWHdk2BU2vNIe2wnfIHSS2K/H0XXRDop7VMlccLevaMMM7UvAoL+02XF6a06EsSbgBmPaVFIrGioFW7upYGEHG99N94tPAHsrc6b39+ZWbBKfrT9HqQq1KT1/LqmmgsgRo7l+/ki65xVmBSqR75w2mmqXAzh86CAR4aHcuHGD06dPA1CmTBlsbGwYP16bLKDT6VixYgUArq6uTJ48GYD33nsvw/qjo6Mf+yB4eJs06elpxC1atGD37t2A1t1jamqKl9fTQ6di791i+/TxXO/SnXJHQgnqWQ+fbftp0GNE5p947E1Y0huECbz0N9j9Ny485E48H685g2+5YrzVukLm61SU/MjWVVvo7foBOLUiw2LeHk682bICK4+FsuPcTTxrNMb51x8xS5GEvDqc8Ktn8izk/KBAJXLrcp50cnLAx9wCj9IeVK3oRVmP//q8N27cyKJFix61pH/99VcAli5dyrx587CzsyMsLOzp6e3ZcOrUKTZs2ICVlRVvvfUWX3zxBS4u2peUOp2OkztXsOG1Tlxp3pIyczYTVdoO66W/0OWrpdjYOWX+Qkmx8Gc/uB+pTfhx+S9Za0MNAwBtqKHark0pFOoNhdI+sO0jSMh4Q+YxbStRtYQ97606zd37yVT2bYf1zC+wuZ/CxaGDuX0jMA+DNi6DTAjKqpxMCEpNTeHflT8Ru2gpZQNjSbCAG21rUPetjyhd0Tvdc27fvo2zszMmJiaMGTOGv/76i/Dw8HTL5kRk2BWOLfgeq60Hcb+VTKI5hDX0xGPgMGq17pv1DxBdKiwbAFd2wsBlUPnxQUE/7rjMjB2X+HGANz2805n1qSgF1Y0TMLcN1H8NOn+bYbFzN2Lo8dMBOtYsyayBdQE4vnUJphOmcdvdCt8Vmx5fp7+AKzQTgkxNzWjSfywvbDyC7vevCK9TCs8tZ7nbbSAbh7Tl7IG1T53j5+eHra0t1tbWLFy4kCVLlhgsngfJiRxe+RObBrQkon03yi3ewwMbc2683YuK+/fS9ffN1Gn7Yvb+Ctj5mTbhp/O3TyXxY9fvMnPXZXp6l1JJXCl8StUF3+Fw9He4EZBhseqlHBjTphLrT95g4ymtcVbvhZdI/HwsxcMT+fel7sRFR2V4fmFR4Frk6Qm9fIKAn6dRaudZrJMhuII99i8PpmHftzA1NTPYddK6fv4IZxb+SLFdARSL0RFjK7jdogZVXh5Fxbqtc36B0yu19Zp9hmnDDNOITXxA55n7kRI2jW2uNk9WCqeEezDbF5zKwvAdkEFjKCVVR59fDhJ8J55t41riZq+t8rl/8Tc4f/EHIZWdaO23HUtru3TPL0gyapEXikT+UHRUOIfnfI7dmn04R6cS6WpOct8XaDz8fWztnXNcf3zcPY4sn0nSmo2UvRyDTkBwdWfse/eiQe9RWFgbaCbljQCY31HbGuvldWBm8djh8csDWBMQxoqRjfH1zPnzUpR86+RyWD0Cuv4AvsMyLHblViydZx6gZWU35gzxeTSybdevH1Pyh5VcrVOc9ou3Ym5hlVeR54oikcgfSk6K59DS6aT8uYpSoQnEWQsiO/ng/dokLK3tSEqIIyk+lqT4WB7Ex/Eg4T7JCXGkJMSTkhhPakICqYkJ6BIT0CUmIZOSkHfuUeLYdWwTJVHFzIh9oQF1Xn4n47VQsisuEua00n4esRvsHp8xujYgjLF+AYxtW4lx7Ssb9tqKkt9ICQu6ws0z8PYxbVRLBubuC2LapvNMf7EOvev9Nwhi23fv4PH7VgIblaHTvM259ld6XihSifwhnU7HiW1LuTlvDuVP3852PSkmkGQhiKhdihL9B1O345Dc+c+QkgyLesCN4/DqFq2fMI2QO/F0/nE/ldztWDGysRqlohQNt87Dr82gzgDo8VOGxVJ1kv6/HeLizVi2jWtBSUfrR8c2T3kNT79/CGxdkc4/rTXIyDVjKJKJPK3AU/u5smk5wswMUytrTCytMLe2wczaBnNrO8xtbLGwtsPCxg4rGwcsrG2xsnXEysYBM3OL51/AEDaMB/950Pt3qN3vsUMpqToGzj3M+fBYNo9tjoezWhBLKUK2fQwHZ8KrW6FsowyLXbt9n04/7qd+eWcWDqv/qItFp9OxefIgvNafJKhLbTp9u6xAJvMin8jzPf/5sGEcNBkDHT5/6vCsnZf5fvslZvSvQ6+6Ga8XoyiFUlIc/NQArIvBiL3a2iwZWHToGp+sPcuXvWsxsEHZR4/rdDo2je5JhV2XufZiYzp9Nj8PAjesQjP8sFC6fkhbi7liO21W2xOOB9/lh52X6V6nFD3VUEOlKLK0g45faX3lR569bO1LDcvRpIIL/9twjpA78Y8eNzExodOsVdr+nysOsfmz13M76jyjErmxRYfCiiHgVA76/A4mj+9yH5eUwjt+AZRwsOLznjUNsgSvohRI1bppjZ3dX0BMxhP6TEwE3/StjRCCiStPotP91+tgampGx7kbCGxQGs8/DxSaZK4SuTE9SAC/QdrGswOXaX82PuHz9ecIvRvPDwO8cbRW48WVIkwI6PQNpCbDtg+fWbRMMRs+7lqNw0F3WHjo2mPHzMwt6DhvU6FK5iqRG4uUsG4MhJ+CPnPB7elVEP8NimK5fwivt/CivhovrijaWkPNxsGZvyFw9zOLvujrQesqbny95QJBkY9vVlHYkrlK5MZycBacXgGtP4QqnZ46nJyi46M1ZyjtZM3YtpWMEKCi5FPNxmldkds/znCpWwAhBF/1qY2lmSnv/nWSVN3jZQtTMleJ3Biu7IAdn0L1HtDi3XSL/H4giMu34visRw1sLAruBAZFMThzK2j1HkSchgsbn1nU3cGKz3rU4HjwPebuD3rqeGFJ5jndfPlbIcQFIcQpIcRqIUQW1mctoqICYeWrULw69Pj50VZtaYXciWfmzst0qO5O22ruRghSUfK5Wi+Csxfs/eqZrXKA7nVK0bFGCaZvu8Slm7FPHS8MyTynLfLtQE0pZW3gEvB+zkMqxBJjYNlAEKYwYOlju/w8JKXk03VnMRGCT9Xem4qSPlMzaDEpU61yIQT/61UTeyszxq8I4EGq7qkyBT2Z5yiRSym3SSlT9HcPA2qmSkZ0Olg9EqKuQL8FUMwz3WJbz95k14VbjGtXmdJO1umWURQFqNUv061yVztLpvWqyZmwGH7enf6GEwU5mRuyj/xVYLMB6ytc9nwJFzdBxy/Bq2W6Re4npTB1/VmqlrDnlaaeeRufohQ0j7XKNzy3eMeaJenhXYpZuy5zJiw63TIFNZk/N5ELIXYIIc6kc+uRpsyHQAqw9Bn1jBBC+Ash/CMjIw0TfUFxbi3s+wa8X4IGGe/V+cOOS4RHJzKtV03M1YJYivJ8tfqBcwXY87X2V+9zTO1eA2dbC9796yQp6XSxQMFM5s/NFlLKdlLKmunc1gIIIYYCXYHB8hkLt0gp50gpfaWUvm5ubhkVK3wiL8HqN6G0L3Sdnu6XmwDnw2OY/881BjbwwKecGjOuKJliagYtJ8HN03Dx2X3lAE42FkztXoMLEbEs/Tc4w3IFLZnndNRKR2Ay0F1KGf+88kVOcjz8NVQbLtV/MZhZpltMp5N8uPo0jtbmTO5YNY+DVJQCrmbfLLXKO9YsQdOKLny/7SJRcUkZlitIyTynf7/PBuyB7UKIACHErwaIqfDYPFFbS7n3HHAolWGx5f4hHA++xwedq+Fkk0dL5ipKYZG2VZ6JvnIhBFO61SA+OZXvtl18ZtmCksxzOmqlopTSQ0rprb+9YajACryAP+HEEmg+QVvoJwO345L4avMFGpZ3pk89tbKhomRLzb7gUhH2Zq5VXsndnqFNPPE7GsKp0HvPLFsQkrn6Ri033DqvbRJRrhm0evbQ+i82ned+Ugr/UysbKkr2PRzBcvNMplrlAGPbVcLF1pIp684+tkJiep5M5htG9+RBcqIhIjcIlcgNLSkOVgzVJvv0nffMBfAPBUax6ngYI1p4UcndPg+DVJRCqGafLLXKHazMmdyxCseD77H6RNhzyz9K5p1qUGHHRXYMaEvMnQhDRJ5jKpEbkpSwcQLcvqStLW5fIsOi2qJYpylTzJq326hFsRQlxx5rla/P1Cl96pXB28OJLzdfIDbxwXPLm5lb0HXGSsJGdcfj/B2O9e5I6OUTOY08x1QiN6QTi+GUn7agj1erZxaduz+IwMj7fNajBtYWps8sqyhKJtXS95VncgSLiYlgavcaRN1PYubOy5m+TLsxX5PwzQQc7iYTOvAlTu/5OydR55hK5IYScUbbrs2rFbSY+MyiwVHaolgda5SgTVW1KJaiGIyJKbScDLfOZrpVXsfDiRd9PPjjn2tcuRX3/BP06nd7DccFP5NiboJu9EfsW/RldqPOMZXIDSEpVhsvbuUIvec+tV1bWlJKPll3BjMTwafdq+dhkIpSRNTsAy6VMt0qB5jYsQrWFqZMXX+WrGxIX7FuK6r+vYabHna4fbGIzVOGo8vkNQ1JJfKckhLWj4U7QdB3PtgVf2bxLWci2HMxknHtK1PSUS2KpSgGl7ZVfn5dpk5xtbNkfPvK7L98m23nbmbpcq6lKtD8750ENiyDp99BNg3vRFJC5lv2hqASeU75z9e2nWr9AXg2e2bRuKQUpq4/R7WSDrzSxDNv4lOUoqhmb61VnskRLAAvNSpHZXc7Pt9wjsQHqVm6nLWNA53/2MrVFxtR4VAwe/u0ISr8anYizxaVyHMi/CRseR8qtIVmE55bfMb2S9yM1RbFMlOLYilK7nnUKj+X6Va5uakJU7rXIPRuAr/tfXo3oede0sSEzp/9QeTkIbhfj+Vcn+4Entqf5XqyQ2WT7EqM1saL27jo+8Wf/VKevRHNH/9cZUD9stQrWyyPglSUIqxmb3CtnKVWeZMKrnSpVZKf91wh9G72lo9qMewDmDUVi6RU7r48Ev+Nf2SrnqxQiTw7pIR1b8O9YOj3B9i6PLO4tijWGYrZWDC5Y5U8ClJRirjHWuVrM33aB12qIYQ26zq76rR5kZJLFxLnYI7lxG/YOTt3N09TiTw7jszV1hhv+wmUbfTc4suOBhMQco8Pu6hFsRQlT9XopbXKszCCpbSTNW+1qsim0xH8c+V2ti9dtmp96q7eQljlYpSavYaNE/qTmpry/BOzQSXyrAo7Bls/gEovQJMxzy0eHf+Ab7ZcpJGXM73qqkWxFCVPPWyVR57PUqv89RZeeDhbM2Xd2XT3+MwsR5eStF2xi8A2lfDaeIqtg9oSFx2V7foyohJ5ViTchb9eATt36PXrc/vFAX7ec4WYxAd80rWGWhRLUYyhRi9wrZKlVrmVuSkfd6nO5VtxLDp0PUeXN7ewovPsNYQM70C5U7c4uWlhjupLj0rkmSUlrHkLYm5omyfbPH8Xn7B7Cfxx8Bq9vEtTvZRD7seoKMrTTEy19cojz8O5NZk+rX11d1pUduOH7Ze4/YwNKDIVgokJHSb+iPWKuTQdOD5HdaVbv8FrLKwO/6xtJdX+M/Con6lTZmy/BBLGd6icy8EpivJMD1vlWRjBIoTg027VSUxJ5ZstFwwShletZ881yS6VyDMj7Bhs/wSqdoVGozJ1yoWIGP4+HsrQJuUoU8wmlwNUFOWZHrXKL2SpVV7BzY5Xm5ZnhX8oASHP3oDCmFQif56kOPj7dbArAT1mZ7h58pO+3nwBO0sz3mpdMZcDVBQlU2r0Areq+lZ55mdujm5TETd7Sz5de+a5G1AYi0rkz7NlsraOSu/fwDpzE3kOBUax+2Iko1pVVMMNFSW/yGar3N7KnPc7VeVkaDQrj4fmYoDZpxL5s5xdo993c/xz11F5SErJV1suUNLRimFNPXM3PkVRsqZ6T32r/FttAEMm9apbGp9yxfhmywViMrEBRV7LUSIXQnwuhDglhAgQQmwTQmS8VXxBEx0K68dAaZ/n7ruZ1qbTEZwMuce49pWxMlcbRihKvmJiqs3/iDwP1zK/DooQDzegSOaH7ZnfgCKv5LRF/q2UsraU0hvYAHxigJiMT5cKq0Zq//aeC6bmmTrtQaqOb7deoLK7HX3qlcnlIBVFyZaavcHaGY7MydpppR0ZUL8sCw9d4+rt+7kTWzblKJFLKWPS3LUF8uc3AVn1zw9w/QB0+gZcKmT6NL8jwVyLimdyx6qYmqjJP4qSL5lbQ72X4cJGuBeSpVPHt6+MhamJNrQ4H8lxH7kQYpoQIgQYzDNa5EKIEUIIfyGEf2RkZE4vm3vCjsHuL7RvuL0HZfq0uKQUftx5mQblnWlT9dmbSyiKYmT1h2v/+s/P0mlu9pYMa+rJ+lM3OB8e8/wT8shzE7kQYocQ4kw6tx4AUsoPpZQewFJgdEb1SCnnSCl9pZS+bm5uhnsGhpQUB3+/pg017Doj00MNAebuC+J2XDLvd6qqpuIrSn7nVBaqdIbjC+FBYpZOHdmiAnaWZny/Lf+0yp+byKWU7aSUNdO5PbkCzZ9An9wJM49sngx3r0HvOZkeaghwKzaRufuD6FyrBHXVWuOKUjA0eB3io+Ds6iyd5mhjzsgWXuw4f5MTwXdzKbisyemolUpp7nYHDDOP1RjOroaAJdBsPHg2zdKps3ZeISlFx7sd1FrjilJglG+pLXGbxS89AYY1LY+LrQXfbbuYC4FlXU77yL/Sd7OcAjoAYw0QU967F6JtoFzaF1q9l6VTr96+z7IjwQxs4IGXm10uBagoisEJAQ1GwI3jEOqfpVNtLc0Y1boi/1yJ4mBg9tcsN5Scjlrpo+9mqS2l7CalDDNUYHlGlwqr9UMN+2R+qOFD3269gIWZCWPbqoWxFKXAqTMALOyz1Sof3LAsJR2t+G7rRWQWJhflBjWz88AMuP4PdP4WnL2ydOqJ4LtsOh3B6829cLO3zKUAFUXJNZb22ui0s6sh7laWTrUyN+XtNpU4HnyPXReydq6hFe1EHnoM9nwJNXpDnYFZOlVKyVebL+BqZ8HrLbL2AaAoSj5S/zVITdZGsGRRP98ylHOx4bttl4y6oFbRTeRJcbDqNbAvmeWhhgC7L97i36t3GNO2EnaWZrkUpKIouc6tMni1hqPzITVr66iYm5owrl1lzofHsPF0eC4F+HxFN5E/NtTQKUunpuokX2++iKeLDQMblM2d+BRFyTsNR0LsDW22ZxZ1q1OKKu72zNh+iZQc7O+ZE0UzkT8cath8ApRrkuXTVx0P5eLNWCa+UBVz06L5EipKoVKpgzZJ6MjcLJ9qaiIY36EyQbfvs+qEccZ7FL0slHaoYcvJWT498UEq07dfok4ZRzrXKpELASqKkudMTLW+8usH4ObZLJ/eobo7dco48uOOyySlZH7TCkMpWok8h0MNARYcvEZ4dCLvdaqmpuIrSmFSdwiYWWVrKKIQgndfqELYvQT8jmRtIS5DKFqJ/NFQw++yPNQQ4F58Mj/vvkLrKm40ruCSCwEqimI0Ns5Qqx+cWgEJWZ9636yiK428nJm16wrxySm5EGDGik4iD9MPNazZR5sEkA0/7wkkNimFyZ2qGjg4RVHyhQavw4N4CPgzy6cKIZj4QhVuxyWx8OD1XAguY0UjkSfHaxtF2JWALtOzPNQQIOxeAgsOXqN33TJULeGQC0EqimJ0JeuARyPtS09d1keg+JRzpnUVN37dG0h0Qt5tCVc0EvnOqRB1GXr+nOWhhg99r18cZ3wHNRVfUQq1hiPg7lW4siNbp0/oUIXohAfM2x9k4MAyVvgTedAe+PdXaPgmeLXMVhWXb8ay+kQYw5p4UtrJ2rDxKYqSv1Ttpv31no0vPUHbEq5LrZLMO3CVqLgkAweXvsKdyBPuwZpR2lKV7T7NdjU/7b6CtbkpI1tmfts3RVEKKDML8B0GV7ZDVGC2qhjXvjIJD1L5ZU/2zs+qwp3It7wHsRHQ61dtn75suHb7PutO3uClRuVwtrUwcICKouRLPq+AiRkcnZet0ysWt6N3vTIsOnyd8OgEw8aWjsKbyM+tg5PLoMVEKO2T7Wp+2ROImakJrzUvb8DgFEXJ1+xLQPWecGKJti5TNoxtWwkpJbN2XTFwcE8rnIk87hZseAdKekOLd7NdTdi9BP4+HsrA+h4Ut7cyYICKouR7DUZAUjScXpGt0z2ctbWYVhwN4XrUfQMH97jCl8ilhHVjtE/R3nOyNXvzod/2BiIEjFB944pS9Hg0gBK1taGI2dw4YnTripiZCn7YcdnAwT2u8CXygKVwaTO0mwJu2d9D81ZMIn5HQ+hTr4waqaIoRdHDreBundNmhGdDcQcrhjbxZE1AGJduxho4wP8UrkR+9zpsfg88m0PDN3JU1dz9QaSk6nizlWqNK0qRVasvWBeDf3/LdhVvtKiAnYUZ07ddMmBgjzNIIhdCvCuEkEIIV0PUly06nTbUELSJPybZf2p37iez9N9geniXppyLrYECVBSlwDG3hnova+uUR4dmq4pitha81tyLLWcjOBV6z8ABanKcyIUQHkB7IDjn4eTAv79oS1B2+lpbVzgH/vjnKgkPUhmlWuOKovgOB6kD/z+yXcWrzTwpZmPOd7nUKjdEi3wGMAkw3oZ1ty7AjqlQpYu2kWoORCc8YME/1+hYowSV3O0NFKCiKAVWsXJQpRMcWwAp2ZupaW9lzqhWFdl3KZJ/g6IMGx85TORCiO5AmJTyZCbKjhBC+Ash/CMjI3Ny2celPoDVI7TdsLv9mK0FsdJafOgasUkpvNW6omHiUxSl4GswAuJva7uLZdOQxuV4pYknZZxtDBiY5rm7BgshdgDpbYXzIfAB0CEzF5JSzgHmAPj6+hqu9b7vWwg/Cf2XgJ1bjqq6n5TCvANXaVO1ODVLOxooQEVRCjyvVuBSSVt/JZvLYFuZmzKlew2DhvXQc1vkUsp2UsqaT96AIKA8cFIIcQ0oAxwXQuTd/mehx5gy9XO+u1odqnXLcXV//hvM3fgHWWqNr1y5EiEE/v7+Ob6+oij51MOhiGHHIPSYsaN5Sra7VqSUp6WUxaWUnlJKTyAUqCeljDBYdM+SHK9t22ZpB1W75Li6xAepzNkfRNOKLviUK5apc2JjY5k5cyYNGzbM8fUVRcnn6gwACzs4mvUNmnNbgRxHPm3aNKp4laHdjAAuWtQGc236fGBgIB07dsTHx4fmzZtz4cKFR483atSI+vXr88knn2BnZ/dUnSv8Q4iMTWJ060qZjuPjjz9m0qRJWFmp6fuKUuhZOWiDKc78DfF3jB3NYwyWyPUt89uGqi8jx44dw2/RPE4MTWHVt29z9Px/WyqNGDGCWbNmcezYMb777jtGjdLGlY8dO5axY8dy9OhRSpUq9VSdySk6ft0TSNzKD3izT1u8vb0fu+3Y8fQC8ydOnCAkJISuXbvm3pNVFCV/qTsEUpPh7CpjR/KY537Zmd/s37WNXuVisSlZBbp/Qfd/tFEqcXFxHDx4kH79+j0qm5SkDRU6dOgQa9asAWDQoEG8++7jC2mtPhHKjehEVm7cTqsqxZ8bg06nY9y4cSxYsMBAz0pRlAKhRC0oXh1OLof6rxk7mkcKXCLn/HpEctxTa4zrdDqcnJwICAjIUnUpqTp+2RNIrdKOfPxaH2Jjn14P4bvvvqNdu3aP7sfGxnLmzBlatWoFQEREBN27d2fdunX4+vpm73kpipL/CQG1+8OOT7VNJ1zyx6TBgtVHfm4dLcxPsTrEiQTn6sTGxrJ+/XoAHBwcKF++PH/99RcAUkpOntSGtzdq1Ii///4bAD8/v8eq3Hg6nGtR8bzVuiL79+8nICDgqVvaJA7g6OjI7du3uXbtGteuXaNRo0YqiStKUVGrHyDg1HJjR/JIwUrkN89Qz8eX/q+8ibe3N3369KF58+aPDi9dupR58+ZRp04datSowdq1awH44YcfmD59Og0aNCA8PBxHR22MuE4nmb3rCpXd7ehQ3d0oT0lRlALGsTSUb6El8mwub2toQhohEF9fX5ntcdcPEh+NUsms+Ph4rK2tEULg5+fHsmXLWLt2LVvOhPPGkuP8OMCbHt6lsxePoihFT8CfsOZNeHUrlG2UZ5cVQhyTUj71p3/B6yPPYhIHbaTL6NGjkVLi5OTE/PnzH23BVN7Vlq61nx7JoiiKkqFq3WDDeDjpl6eJPCMFL5FnQ/PmzR/1lz+0+8Itzt6I4Zu+tTE1ydn6LIqiFDGW9lCtq7b2SqevwczSqOEUrD5yA5FSMnPXZUo7WdOrrupSURQlG2oPgMR7cGmrsSMpmon8UGAUJ4Lv8UarCpibFsmXQFGUnPJqBbbF88XolSKZxWbtukJxe0v6+ZQxdiiKohRUpmbaUMRLW40+Zb/IJXL/a3c4FBTFiBZeWJmbGjscRVEKsjr9QffA6FP2i1win737Cs62FgxqmLPt4BRFUShRG9yqaVP2jahIJfIzYdHsuRjJ8GblsbEoEgN2FEXJTUJorfLQI9qUfSMpUol89q4r2FuZMaRxOWOHoihKYVHrRbQp+yuMFkKRSeSXbsay5WwEw5p44mBlbuxwFEUpLBxLQ/nmRp2yX2QS+e/7g7AyN2FY0/LGDkVRlMKm9gC4exVCjhjl8kUikUfFJbEm4AZ96pWhmK2FscNRFKWwqd4dzKzhlN/zy+aCIpHI//w3mOQUHcOaeho7FEVRCiNLe23v4DOrICUpzy9f6BN5coqORYev07KyGxWL2xs7HEVRCqs6+in7l7fl+aVzlMiFEFOEEGFCiAD9rbOhAjOUjadvEBmbxKvNVN+4oii5yKu1NmX/ZN53rxiiRT5DSumtv20yQH0GI6Vk3oGrVCxuR4tKrsYOR1GUwszUDGr1NcqU/ULdteJ//S5nwmJ4pYknQqilahVFyWW1H07ZX52nlzVEIh8thDglhJgvhCiWUSEhxAgh+gNZJgAAC7JJREFUhL8Qwj8yMtIAl32++Qeu4mhtTu96aqlaRVHyQMk64FY1z1dEfG4iF0LsEEKcSefWA/gFqAB4A+HA9xnVI6WcI6X0lVL6urm5GewJZCTkTjxbz0YwsEFZNR1fUZS8IYTWKg/5F+4E5dlln5vIpZTtpJQ107mtlVLelFKmSil1wFygQe6HnDmLD19HCMHLajq+oih5qXbeT9nP6aiVkmnu9gLO5Cwcw7iflMKyI8F0qlmCUk7Wxg5HUZSixLEMeDbTRq/k0ZT9nPaRfyOEOC2EOAW0BsYZIKYc+/t4KLGJKWrIoaIoxlFHP2U/9GieXC5HiVxKOURKWUtKWVtK2V1KGW6owLJLp5P88c81vD2cqFc2w+9eFUVRck+17mBmlWdjygvd8MM9l25x9fZ91RpXFMV4rBy0KftnV0FKcq5frtAl8vkHrlHCwYpONUsYOxRFUYqy2gMg4W6eTNkvVIn8YkQsB67c5uX/t3fvsVnVdxzH3x+5TAWhclGpOkFTiXZI5bZOeeZlZEEyxwzZpDplw4SAYaljG2FZQkzMkrmwxi2yNU66oeGyGPCKiG6dd8UVbAG5WFGm1a50Smi5rB32tz/OaT3r7bmd0+c55PtKnvT0nPP7ne/zOz++nP6ec37PNZcwZNBp9daMMXFz2Y0wbOyAzIh4WmW7P732AWcOOYOy6fZ9nMaYHBs0GL7iP7J/8kikhzptEvlnx9t5/O2PueVqm3PcGJMnJt8Kn7dH/sj+aZPIN7z1IW2nOlhoc44bY/LFuBIYMxHqon1k/7RI5O2nOnjkjUMkisZQdL7NOW6MyROSd1X+0Zvw2QeRHea0SORb9zTS1GJzjhtj8tCk73k/I3xkP/aJvHPO8UvHDuO6ougn4zLGmLQUXAzjE97dKxE9sh/7RL7zwyPsajjKD6+dwBln2Jzjxpg8dNWt3myIDTWRVB/7RF716iFGnDmYeTbnuDEmX10513tkP6J7ymOdyBuOnGDrnkbKvmpzjhtj8tiZI2DiHNizKZJH9mOdyB99o3PO8fG5DsUYY/o32X9k/70XQq86ton8RLs35/js4gu40OYcN8bku8tuhFseggnXhV51bMcjNu38mJb/nGLhzPG5DsUYY5IbNMS7pzwCsbwi9+Yc/4DJF43kqarfsmrVqgGPobKykkmTJlFSUsLMmTPZu3fvgMdgjDEQ00T+Un0z7zd7c45Lubnl8LbbbmP37t3U1tayfPlyli1blpM4jDEmlon8xytW0rRmMb/9yQIOHDjQtf7gwYPMnj2bqVOnkkgk2L9/f9f60tJSpk+fzsqVKxk+fHjWMYwYMaJr+fjx4zn7D8UYY2I3Rr75+ZfZ++pz3P/Is9x17SVMmTKFqVOnArBo0SIqKyspKipi+/bt3H333VRXV1NeXk55eTllZWVUVlb2WXcikaC1tbXH+lWrVjFr1qwe61evXk1FRQXt7e1UV1eH9yaNMSYNclk+MirpR8BS4BSwxTm3PFmZadOmuZqazJ5wmrVwOW/XN1D/3J8ZNWwoy5Yto7CwkMWLFzN27FgmTpzYtW9bWxv79u1j9OjRNDU1MXjwYFpaWigsLOTYsWMZHb8369evZ9u2baxduza0Oo0xpjtJO5xz07qvz+qKXNINwFzgKudcm6TzsqkvmSPH29nVcJQrxo1gVLc5xzs6OigoKKC2tjbj+tO9Iu80f/58lixZkvFxjTEmG9mOkS8BfuWcawNwzh3OPqS+rX/rQwYVXknznlc5efIkra2tPP3004A3Zj1hwgQee+wx/Fioq6sDoLS0lE2bNgGwcWPfj8i+8sor1NbW9nj1lsTr6+u7lrds2UJRUVFo79MYY9KRbSK/HEhI2i7pJUnT+9pR0iJJNZJqmpubMzrY2HO+xB3fup47by+jpKSEefPmkUgkuravW7eONWvWMHnyZIqLi3nyyScBeOCBB6ioqGDGjBk0NjYycuTIjI4f9OCDD1JcXExJSQkVFRU2rGKMyZmkY+SS/gr09pX0vwB+CVQD5cB04C/ApS5JpdmMkWfixIkTnHXWWUhi48aNbNiwoSvJG2NMXGQ8Ru6c63NwWNISYLOfuN+S1AGMATK75I7Ijh07WLp0Kc45CgoKqKqqynVIxhgTmmxvP3wCuBF4UdLlwFDg31lHFbJEItE1Xm6MMaebbBN5FVAlaQ/QDixINqxijDEmXFklcudcO/D9kGIxxhiTgVg+om+MMeYLlsiNMSbmLJEbY0zMWSI3xpiYy3rSrIwOKjUD/8yw+Bjy8BZHLK50WVzpsbjSk69xQXaxXeKcG9t9ZU4SeTYk1fT2ZFOuWVzpsbjSY3GlJ1/jgmhis6EVY4yJOUvkxhgTc3FM5A/lOoA+WFzpsbjSY3GlJ1/jgghii90YuTHGmP8XxytyY4wxAZbIjTEm5vIykUv6rqR3JHVI6vM2HUmzJR2Q9J6kFYH1oyS9IKne/3luSHElrVfSREm1gVeLpHv8bfdK+jiwbc5AxeXvd0jSbv/YNemWjyIuSRdL+rukff45Lw9sC7W9+uovge2S9Dt/+y5JU1ItG3Fct/vx7JL0uqTJgW29ntMBiut6SUcD52dlqmUjjutngZj2SPpc0ih/WyTtJalK0mF/Jtjetkfbt5xzefcCrgAmAi8C0/rYZxBwELgUbx70OuBKf9uvgRX+8grg/pDiSqteP8Z/4d3ED3Av8NMI2iuluIBDwJhs31eYcQHjgCn+8jnAu4HzGFp79ddfAvvMAbYCAkqB7amWjTiua4Bz/eWbOuPq75wOUFzXA89kUjbKuLrtfzNQPQDt9XVgCrCnj+2R9q28vCJ3zu1zzh1IstsM4D3n3PvOm053IzDX3zYX6PwSzbXAd0IKLd16vwEcdM5l+hRrqrJ9vzlrL+dco3Nup7/cCuwDLgzp+EH99ZdgvI84z5tAgaRxKZaNLC7n3OvOuSP+r28CF4V07Kziiqhs2HWXARtCOnafnHMvA5/1s0ukfSsvE3mKLgQ+CvzewBcJ4HznXCN4iQI4L6RjplvvfHp2oqX+n1ZVYQ1hpBGXA56XtEPSogzKRxUXAJLGA1cD2wOrw2qv/vpLsn1SKRtlXEF34V3ZderrnA5UXF+TVCdpq6TiNMtGGReSzgZmA5sCq6Nqr2Qi7VvZfkNQxtTPlzo751L5ZmT1si7reyn7iyvNeoYC3wZ+Hlj9B+A+vDjvA34DLBzAuK51zn0i6TzgBUn7/SuJjIXYXsPx/sHd45xr8Vdn3F69HaKXdd37S1/7RNLXkhyz547SDXiJfGZgdejnNI24duINGx7zP794AihKsWyUcXW6GXjNORe8Uo6qvZKJtG/lLJG7fr7UOUUNwMWB3y8CPvGXmySNc841+n++HA4jLknp1HsTsNM51xSou2tZ0h+BZwYyLufcJ/7Pw5Iex/uz7mVy3F6ShuAl8XXOuc2BujNur17011+S7TM0hbJRxoWkq4CHgZucc592ru/nnEYeV+A/XJxzz0r6vaQxqZSNMq6AHn8RR9heyUTat+I8tPIPoEjSBP/qdz7wlL/tKWCBv7wASOUKPxXp1NtjbM5PZp1uAXr9hDuKuCQNk3RO5zLwzcDxc9ZekgSsAfY55yq6bQuzvfrrL8F47/TvMCgFjvpDQqmUjSwuSV8GNgN3OOfeDazv75wORFwX+OcPSTPw8smnqZSNMi4/npHAdQT6XMTtlUy0fSvsT2/DeOH9o20A2oAmYJu/vhB4NrDfHLy7HA7iDcl0rh8N/A2o93+OCimuXuvtJa6z8Tr0yG7lHwV2A7v8kzVuoOLC+1S8zn+9ky/thTdM4Pw2qfVfc6Jor976C7AYWOwvC1jtb99N4I6pvvpaSO2ULK6HgSOB9qlJdk4HKK6l/nHr8D6EvSYf2sv//QfAxm7lImsvvIu2RuC/eLnrroHsW/aIvjHGxFych1aMMcZgidwYY2LPErkxxsScJXJjjIk5S+TGGBNzlsiNMSbmLJEbY0zM/Q9+IVu5ZcjGlQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "a = 3\n", "x0 = 1\n", "x0p = 1\n", "\n", "h=.1\n", "T = np.arange(-1,1+h,h)\n", "plt.figure()\n", "plt.title('Approximations to Solutions')\n", "\n", "for deg in range(3,7):\n", " X = np.zeros(T.size)\n", " for k in range(0,deg):\n", " X = X + A(k,a,x0,x0p)*(T**k)\n", " plt.plot(T,X)\n", " plt.annotate(f\"deg = {deg}\", xy=(-1,X[0]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The *Legendre Polynomial of degree $n$* is the polynomial solution to the equation with $\\alpha = n$ and with $P(1)=1$. \n", "We will denote these as $P_{n}(t)$ for now. Graph $P_0,\\ldots,P_5$ below. To normalize, we will divide ```X``` by ```X[-1]``` since ```x[-1]``` is the last entry of ```X``` and this corresponds to $x(1)$. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3gUVRfG35uEhBp67wgqQVCKBQQVBUEUFETF3hEVEf1UiqJgx45gQ0XsiiIqgoCiVAVEepEWQgmdEEJ6snu+P94ddpNs35ndDbm/59knm507M3dnZ+6591QlItBoNBpN2SUm0h3QaDQaTWTRgkCj0WjKOFoQaDQaTRlHCwKNRqMp42hBoNFoNGUcLQg0Go2mjKMFgaZMo5SaqpR6PsJ9uEQptTeSfSiOUqqJUipTKRXrR9uo678mMLQg0ISEUipFKdUj0v2IBhwDot0xgJ5QSm1RSt0Z6X4Fg4jsFpHKImKLdF801qMFgUbjAaVUXBC77RORygASAYwA8KFSKsncnmk05qIFgcYSlFIxSqmRSqkdSqmjSqlpSqkaLttvU0rtcmwb47qyUEqNdbT/zDGz3qiU6uSybwOl1HSl1GGl1E6l1DCXbb72ba+UWuXY9i2A8i7bLlFK7VVKjVBKHQDwia/v4QkhPwI4BiBJKZWglHpLKbXP8XpLKZXg5ro9rpSaXuyziUqptxzvFyilnlNKLXV8h3lKqVoubfs5vnO6o21rl20pjuOvU0plKaU+VkrVVUr96jjW70qp6o62zZRSYghDpdSdSqnNjnbJSqn7PH13x/VLdVkVXebremkiixYEGqsYBuAaABcDaAAOiO8AgGOG/C6AmwHUB1AVQMNi+/cD8A2AagB+BjDJsW8MgJkA1jr2uQzAcKVULz/2jQfwI4DPAdQA8B2Aa4udt55jW1MAg719D284BEh/Rx/WA3gSwAUAzgFwNoDzADzlZtcvAPRWSlVzHCcOwA2OPhvcBOBOAHUAxAN4zNH2dABfAxgOoDaA2QBmOr63wbUAegI4HUBfAL8CGA2gFjgeDIN7DgG4Clzp3AngTaVUBzff+wwAQwGcKyJVAPQCkOLhmJpoQUT0S7+CfoEPeQ83n28GcJnL//UBFACIA/A0gK9dtlUEkG8cB8BYAL+7bE8CkON4fz6A3cXONQrAJ37sexGAfQCUy/a/ADzveH+Jox/l/fkebr7zJQDsANIBpAFYA2CQY9sOAH1c2vYCkOKy316Xbb8CuNfx/ioAm1y2LQDwlMv/DwCY43g/BsA0l20xAFIBXOLyW93ssn06gPdc/n8IwI+O980AiLvv6dj+I4CHi/cfQEtQaPQAUC7S96d++fcKRgeq0fhDUwAzlFJ2l89sAOqCM+s9xocikq2UOlps/wMu77MBlHfMjpsCaKCUSnfZHgtgsR/7NgCQKo4Ry8GuYuc9LCK5fn6PVJRkn4g0cvN5g2Ln2uX4zB2fArgfwIcAbkHR1QBQ8vtVdncOEbErpfag6GrroMv7HDf/V4YblFJXAHgGXEnEgMJ7ffF2IrJdKTUcFMhtlFJzATwqIvvcHVcTHWjVkMYq9gC4QkSqubzKi0gqgP0ATg6WSqkKAGoGcNydxY5bRUT6+LHvfgANlVLK5bMmxdoUT8fr7XsEwj5QqLie19Pg+COAdkqps8AVwZfBnMPxPRvDvcDyG4ctYzqA1wDUFZFqoNpJuWsvIl+JSFdHXwTA+FDOr7EeLQg0ZlBOKVXe5RUH4H0ALyilmgKAUqq2UupqR/vvAfRVSnVx6K/HwcOg4oYVADIcBskKSqlYpdRZSqlz/dj3bwCFAIYppeKUUgNAXb03vH2PQPgawFOO/WuB6rEv3DV0rEi+B/AVgBUistvPc0wDcKVS6jKlVDkA/wOQB6q/QiEeQAKAwwAKHauDy901VEqdoZS61CE8csFVhnZBjXK0INCYwWzwgTdeYwFMAA2185RSJwAsA/X7EJGNoD76G3CWfgLUK+f5OpHQr70vaHTdCeAIgI9Ag7OvffMBDABwB2j0vQHADz528/g9AuR5ACsBrANVKqscn3niUwBtUVIt5BER2QKqkiaC16UvgL6O7x00InICNCJPA6/bTeA1cUcCgJcd5z8AGrRHh3J+jfWooupSjSb8KKUqgwbWViKyM9L9iQaUUk0A/AegnohkRLo/mlMbvSLQRASlVF+lVEWlVCVQ97we2s0QwEkX2UcBfKOFgCYcaK8hTaS4GlR7KFBlMkj08hQOwXgQ9P7pHeHuaMoIWjWk0Wg0ZRytGtJoNJoyTqlUDdWqVUuaNWsW6W5oNBpNqeLff/89IiK1i39eKgVBs2bNsHLlykh3Q6PRaEoVSqnikfQAtGpIo9FoyjxaEGg0Gk0ZRwsCjUajKeNoQaDRaDRlHC0INBqNpoxjiiBQSk1RSh1SSm3wsF0ppd5WSm13lMnr4LKtt6Oc3Xal1Egz+qPRaDQa/zFrRTAV3sPhrwDQyvEaDOA9AFBKxYJl/64AK0ndqHShb41GowkrpsQRiMgipVQzL02uBvCZI5fMMqVUNaVUfbAc3nYRSQYApdQ3jrabzOhXCS68ENi8GWjXzpLDazQajVX8c6gpXtl8Fe44fzOuXPa0qccOl42gIVxKEwLY6/jM0+clUEoNVkqtVEqtPHz4cHC92L4dOHYMyMwMbn+NRqMJMyLAd4cvxvubu+F7XId/U/wt5uc/4Yosdld9Srx8XvJDkckAJgNAp06dgsuUN3Ik8OijQO3awK+/BnUIjUajCRf79wMPPADMX5SB2/ApAKD/t4NMP0+4VgR7wdqpBo3A+qqePreG8xxVCefMAXSKCo1GE6WIAFOmAK1bc7ia3fMtpKEGAKB5hxqmny9cguBnALc5vIcuAHBcRPYD+AdAK6VUc0ft2kHwXAIvdFq35t+KFYGnnrLsNBqNRhMsyclAz57A3XcDZ58NbFh4FKf9/QXykIBEdQKVq/hb3tt/zHIf/RosDH6GUmqvUupupdQQpdQQR5PZAJIBbAfwIYAHAEBECgEMBTAXwGYA0xz1bK2hRg2gTh3grLOAuXOBhQstO5VGo9EEgs0GvPUW0LYtsGIF8N57wJ9/AqdNfwV/Zp6LPJRH/conLDm3WV5DN/rYLgAe9LBtNigowkNSEpCTAzRoADz5JLB4MaDMl7AajUbjL5s2cQWwbBnQpw/w/vtA48agkWDiRCxv9jWOpNRC4/oFlpy/7EUWt24N/PcfMGYMsHSpNhprNJqIkZ8PPPcc0L49sG0b8MUXwC+/OIQAALzwAiS/AKknErEXjdC4aawl/SibguD4cYrdFi2A0aMBuz3SvdJoNGWMlSuBc88Fnn4aGDCAq4Kbb3ZRUKSkAJMnY9vAUah2dDsOoB4anlnZkr6UTUEAMKbg2WeBtWuB776LbJ80Gk2ZITsbeOIJ4PzzgSNHgJ9+Ar7+mubLIowbB8TE4LezHkFj7IENcWh4ZqIlfSp7giDJkcFi82Zg0CAajseMAQoLI9svjUZzyrNwIT2BXn2VNoFNm4B+/dw03LAB+Owz4MEH8dvK6qiXcAwA0LCxNUN22RME9esDiYn8BWJjgeefp3Lu008j3TONRnOKkpEB3H8/cMkl1ETPnw9MngxUrephh9GjgSpVUPD4aPzxB1ABuQCAhm7zLoRO2RMESlE9tHkz/+/Xj4FmY8cCubkR7ZpGozn1mDULaNOGA/+jjwLr1wOXXuplh8WLgZkzgZEjsWJHTciJE8jISwAANGpkTR/LniAAigoCpYCXXwb27gUmToxsvzQazSnDkSPALbcAV13Fmf9ffwGvv854Vo+I0IDQoAEwbBh++w1ojf+QioYoF2dHrVrW9LVsCoKkJODAASagA4Du3YErrgBefBFIS4ts3zQaTalGBPjmG843p02jsmHVKhqHffLjjwwmGDcOqFgRv/0G9Gm+GaloiAZ1bIixaMQum4LA8BwyVgUAMH483UpffDEyfdJoNKWe1FTg6quBG2+kd/qqVcAzzwDx8X7sXFgIjBrF8emOO3D8OLB8OdC9/mbsVY3RsKl1OUK1IDBo2xa44w6qh1JSItErjUZTShEBPvyQyobff6cK6K+/6JToN1OmAFu2AC+9BMTFYcECpp1oozZjT7kWaNzEugwIZVMQNGsGJCQUFQQA4wpiYnRCOo1G4zc7dgCXXQYMHgx07Ehj8KOP0inRb7KyqEPq0uWkP+lvvwGVKgE1Dm7G3sJ6zmhjCyibgiA2FjjjjJKCoFEjYPhw4MsvgdWrI9M3jUZTKrDZgDfeoDLh33/pFTR/PnDaaUEcbMIE5hV65ZWTocVz5wKXdc3D0R3pyLUnoEkTc/vvStkUBADXcJvcVMQcORKoWRN4/HGu9zQajaYYGzZw8v6//wE9enAouffeIPNXHjlCG+XVV7OcLpj4YPt24Pr227Bb6DOqVwRWkJREW0DxspVVqzLSeP58YN68iHRNo9FEJ/n51OB06ADs3EnvoJ9+CjHQ64UXOA65OKrMmcO/l9XbiD2O2l1aEFhB27b8u9FN+YP776fJf8QIrv80Gk2ZZ8UKCoBx44Drr+cq4IYbQsxiv3078M47wJ13OtPfgEmRW7YE6h1ejz2qKQAtCKyhXTv+Xbeu5Lb4eErptWtpL9BoNGWW7GyqgDp3pof5L78wXbQpwV0jRnC8ee65kx/l5rIgzRVXAFi3DrtrnI2EBJZatwqzKpT1VkptUUptV0qNdLP9caXUGsdrg1LKppSq4diWopRa79gWvkLCzZrRJL9+vfvt118PdOpED6KcnLB1S6PRRA9//knlwRtv0Cto40bgyitNOvjChcAPP9AuWb9+kY9zchyCYP167Kl0Jho3trZ+VsiCQCkVC+AdAFcASAJwo1IqybWNiLwqIueIyDkARgFYKCKuIbzdHds7hdofv4mJoZOvJ0EQEwO89hqwZw/vAo1GU2Y4fpwD/6WXcihYsIClIxPNygJtt9PHtHFjLjdcmDOH3u0Xt88AUlKwB40tVQsB5qwIzgOwXUSSRSQfwDcArvbS/kYAX5tw3tBp146CwJN30MUXs2LESy8B+/aFt28ajSYi/Pwz1fUff0znwbVrORSYyuefM+z4pZeAChWKbPr1V2YprZi8AQCwO7tWqRAEDQHscfl/r+OzEiilKgLoDWC6y8cCYJ5S6l+l1GBPJ1FKDVZKrVRKrTx8+LAJ3QbXfEeP0n/XE6+8AhQUsL6xRqM5ZTl0iCVKrr6aHuTLl/Px95okLhiysphm+txzmYvChZ07GVxsqIUKEYt9aeVLhSBwp7ny5IDfF8DSYmqhC0WkA6haelApdZG7HUVksoh0EpFOtc2ymhieQ57UQwCjQ4YPB6ZOZdSIRqM5pRChT0hSEjBjBu22K1fSRGgJr71GDcObb6J4FjmjhLohCPZXPh12uyoVgmAvANduNgLgSY8yCMXUQiKyz/H3EIAZoKopPPgjCACuBurUoUDQQWYazSnDnj1A375MF92qFRMKPPWUn0nigiE1lcuM6647GTzmypw5QPPm7AvWrcOeZt0AWOs6CpgjCP4B0Eop1VwpFQ8O9j8Xb6SUqgrgYgA/uXxWSSlVxXgP4HIAG0zok3/UrMm83+5cSF1JTGQlsyVLgO+/D0/fNBqNZdjtwPvvs2DMn38Cb73Fxzspyfe+IfHkk8wyOn58iU05OUxY16cPoCDA+vXYVfdcAEDTptZ2K2RBICKFAIYCmAtgM4BpIrJRKTVEKTXEpWl/APNEJMvls7oAliil1gJYAWCWiMwJtU8B0bat7xUBANx1F43Ljz+uK5lpNKWYbdvoDXT//awRsGED8PDDASaJC4Z//2VJ3OHDOe0vxh9/UBj07QuuHNLTkVKJ6UujXhAAgIjMFpHTReQ0EXnB8dn7IvK+S5upIjKo2H7JInK249XG2DestGvH5HO+itfHxnLasGsXdXsajaZUUVjIovHt2gFr1tAraN48t2Oy+YjQXbR2bRqK3TBzJlC5Mj2GDC3FLjRFzZr83ErKbmSxQdu2QF4epwm+6N4duOYa5gTx5mmk0WiiirVrgQsuYBXI3r2ZHuKuu6wN0irCtGnAokVMde+mYr0II5Yvv5wxBIaWIiWzFpo1s757WhAYBmNfdgKDV1+l4NA1CzSaqCcvjzkkO3WiYXjaNAbzNmgQxk5kZQGPPQa0b88UpW5Ys4baoL59HR+sXw80aoRd+8pZrhYCtCBgtbLYWP/sBAAzQT38MPDJJ/Qx02g0Ucnff3Psff554KabuAq47rowrgIMXnwR2LuX1Q89GCJmzmS/+vRxfLBuHaRtO+zaBb0iCAsJCSxS468gALgaqFMHGDqU7gcajSZqyMqiPfbCC5ndefZs2mhr1oxAZ7ZvZ9zALbe4dRc1mDmTqqs6dcAA1v/+w+HTLkBOjvWGYkALAtK2rf+qIYA6vldeYejh1KmWdUuj0QTG778zhdiECcADDzBJ3BVXRLBDjzzCoIRXXvHYZN8+KhdOqoW2bAEKCpBSm66jekUQLtq1Y5Ga48f93+fWWynhR4wAjh2zrGsajcY3x44Bd98N9OzJcXfRImDSJKBKlQh2avZsWoCffrpIdtHizJrFvycFwdq1AIBdlVoD0IIgfHTowL9r1vi/j1IsKJGWxh9ao9FEhBkzGAj26afM6Lx2LdCtW4Q7lZdHW+IZZ/CvF2bOpPqnTRvHB6tWAeXLI6WQJSq1aihctG/Pv6tWBbbf2Wdz/fnuu4EJEY1GEzIHD7JsyIABQL16rCD20ktA+fKR7hkYa7R9O3VUXvJVZGUBv/0G9OvnYsRetQpo1w679sSiWjW33qamowUBANSty6KjgQoCgH7BNWvScKzzEGk0liMCfPYZHf5++onFBI0yklFBaipdla6+GujVy2vTuXOZqKB/f8cHdjvHoQ4dkJISntUAoAWBkw4dghME1asDL78MLF3KHOMajcYydu+mi+Xtt1MQrF3LQN1y5SLdMxcef5xhzH4UtJoxg/PIk6qsnTuBjIyTgiAc9gFACwInHToA//3HtVqg3HEHk5Y88URgBmeNRuMXdjtNcm3aAIsX0yV/8WLgzDMj3bNi/Pkn8PXXHAtatPDaND+f9oF+/YC4OMeHjsmotO+A5GSfhzANLQgMOnTg3RaIG6lBTAzv0kOHgLFjTe+aRlOW2bKFFcKGDgW6dGGSuKFDS6Tyjzx5ecCQIRy9R43y2XzBAs4bT6qFAAqCuDgcrH0WcnK0IAg/hoIxGPUQAHTsCNx3H6cqwQgTjUZThIICal3PPpvxAFOnMl9/uNQlAfPKK8DWrZwUFis/6Y4ffgAqVaLL60lWrQLOOgvJqQkAWBcrHGhBYNCwITMDBisIAFqtatSgQNARxxpN0KxeTW3rqFHAVVcxPcTtt0cgPYS/bN/O5//665nVzgd2Ow3dffq4eDmJnDQUJyfzI70iCDdKBW8wNqhRgwaiZcuAyZPN65tGU0bIzWXtlnPPZcTt99/zVa9epHvmBRHgwQeZrsbPFPXLlgEHDhRTC+3dCxw5clIQKKW9hiJDhw5UQOblBX+Mm28GLruMkS06VbVG4zdLlwLnnMMcbbfeylXAtddGuld+8O23LGzwwgt+pzX94Qd6Op1MMgc4J6EOQdCwYfhiIkwRBEqp3kqpLUqp7UqpkW62X6KUOq6UWuN4Pe3vvmGlQwe6fW0IoVqmUsB773FqM3y4eX3TaE5RTpwAHnqILpS5ufSt/+QTLrCjnvR05hPq1Iklz/xAhG6jPXoUCxZbtYoW8HbtwuoxBJggCJRSsQDeAXAFgCQANyql3FX+XCwi5zhezwa4b3gI1WBs0KoVM5ROm8Z8IxqNxi1z5zJJ3DvvUBhs2MDiLKWGp56it+D77/td63LVKiA52c1qZ9Uq+sNWqlT6BAGA8wBsd5SdzAfwDYCrw7Cv+TRvDlSrxtqiofLEE4x4eeCB4GITNJpTmLQ0ht/07g1UrMiYgAkTrC/JaCr//MP0MkOH0mvQT6ZNY9xAEfsAcNJQnJvL4OTSJggaAtjj8v9ex2fF6ayUWquU+lUpZaRX8ndfKKUGK6VWKqVWHj582IRuuz1J6AZjg/h44IMPWOP42WdDP55Gc4owfTqTxH3xBQ3Dq1d7TdUfnRQW0juwfn3guef83k2EgqBnz2KqrwMHaB3v2BEpKfyotAkCdw5dxZPurALQVETOBjARwI8B7MsPRSaLSCcR6VS7du2gO+uTDh0YB1BQEPqxunUD7rkHeP11HVugKfPs3091yMCBtKmuXMmUPFGRJC5Q3nqLEmzCBCAx0e/dVq5kxvvrry+2YfVq/nVxHW3e3JSe+oUZgmAvgMYu/zcCsM+1gYhkiEim4/1sAOWUUrX82TfsdOpEr6FAKpZ5Y/x4iv7BgwGbzZxjajSlCBEGgyUlMff+yy8zSdw550S6Z0GybRsLIV9zTcBuTdOm0Vvo6uIK8BUrqJFo3x7btvGjVq3M6a4/mCEI/gHQSinVXCkVD2AQgJ9dGyil6inFUBCl1HmO8x71Z9+wc/75/Lt8uTnHq1GDvsXLl1OfqNGUIVJSmIDzzjudhQBHjHDJrVPasNtZgD4hgc9zABFuhlro8suZq7IIy5YxkVKVKti2jd5EtWqZ23VvhPxziEihUmoogLkAYgFMEZGNSqkhju3vAxgI4H6lVCGAHACDREQAuN031D6FRNOmTEu9bJnf7mA+uekm4MsvnWGS4VzzhUhhIR/ejRsZPb9tG3D4MBMkZmTw5i5fnq+aNYHGjfk64wxq2Vq0iMKcMBrLsdnoCTR6tLOG05Ahp8C98OGHwMKFwMcfe6065o7ly5k9tYRJQYQbHauLbduAli3DG0Vtilx2qHtmF/vsfZf3kwBM8nffiKIUVwVmrQiMY37wASX+PfewsGrUxsozE+706cD8+QzyOXGCn8fEMM9L/frMxnHaafwaublATg696P79l4LCIDGRl7NnT86E2rY9BQYDjVc2b+Zt/tdf9Ar64AOgSZNI98oE9uxhiunLLuMSJ0CmTaMPSb9+xTZs28Zamw5txLZtTsVE2BCRUvfq2LGjWMqLL4oAImlp5h73gw943A8+kLw8kd27RfbsESkoMPc0wZCeLjJhgsj557OLgEhSksiQISJffSWyebNIXp5/x8rOFlm1SuSjj0Tuv1+kTRvnMRs2FHnkEZHly0Xsdmu/kya85OeLPP+8SHy8SI0aIp99dgr9xna7yJVXilSsKJKcHPDuNptIo0Yiffu62fjZZ3w41q+XvDyRmBiRMWNC77I7AKwUN2NqxAf1YF6WC4L583lp5swx9bDZWXb5qPWr0jlmmcTH208OjuXLi1x0kcj774tkZJh6Sp/s3MmBuUoV9qV9e5Hx4/m5mezdK/LJJyL9+omUK8dztWol8tprIkePmnsuTfj591+Rs8/m73r99SIHDkS6Rybz5Zf8cm+9FdTuCxZw9y+/dLPxgQdEKlcWKSyUzZvZ7rPPQuuuJ7QgCISMDBGlRMaONe2QCxaInHYar3gbtVEebz5NJn9gl/ff50CclMRttWuLvP22SGGhaad2y6FDIkOHisTF8XXzzSIrV1p7ToO0NJGPPxa58EKnILzrLq46NKWL7GyRESNEYmNF6tcXmTEj0j2ygIMHRWrWFLnggqAfzHvuEalUSSQz083Gjh1FLr1URER+/pnPxN9/h9BfL2hBEChnnSVyxRWmHGriRMqVFi1E5s4Vsb89kZd+6tSTbex2kaVLeT8AIl26iOzYYcrpi5CXJ/LSSyKJiXx4hwyheipSrF3LPlSowGt07bWcXWqin4ULuaoDRO6+W+TYsUj3yCIGDaK+a+PGoHbPyRGpWlXkllvcbMzO5kxs1CgREXn9dV7PI0dC6K8XtCAIlHvuoaIzRCWnYW7o189lNmCziXTrJlKtmkhqapH2drvI559zU/XqIr/9FtLpi7B6tXP53revyKZN5h07VA4dEnnyST4whnph69ZI90rjjuPHqc0ARJo3F/n990j3yEJ++IFf9Lnngj7E9OniWdO8ZAk3/vSTiHBSVL160KfyiRYEgfLhh7w8IYxGX3zBQ9x8sxuD8Nat1In06+dW2GzfTiNrbCyPEwoFBSJPP82JR716J++5qCQ9nYaySpXY3wceiEIbQl6eyJo1ItOm0To6eLDIjTeKXHWVSM+eItdcI3LrrSIPPUSd8qxZXN6dApbT2bNFGjfm6m34cA+qjlOFgwepq+3QgZbwIBkwQKROHQ9OIa+9xkHCYVS57DKR884L+lQ+0YIgUNavl1CsNqtWcTV58cVevG2Mm8DDSJ+RIdK9Ox+6Dz8Mqhuyfz8N0QDHpqgbVD1w4ACFQGysSK1atCnYbBHqjN0usmIFbUaXXko9lmHpNww7LVvS0n7++SJt24o0a0b9m2u7unVFrrtOZNIkuoyVIo4c4f1jeJNZpcOOGux2CvSEBJENG4I+TFoax4GHH/bQ4LrreK84aNKEE0er0IIgUAoLacl/8MGAd83JoYmhfn0fur7CQpHOnakH8jAw5OSI9O5Nl7Lp0wPrx5Il7EOFCqGvKiLF2rUiXbvyTu3cmeqtsLFxo8gTT/BBBSiR27fnU/3NN1wVeJsS2+3UeS1dStfhW27hdNoQDF26cMUQxdLZbhf59lvKurg4rtZycyPdqzDw6af8jV59NaTDGIqFFSs8NGjcWOSGG0SEt1KIWiifaEEQDJdeSot+gIwcySs7a5Yfjbdtox7ksss8TnkzM+mwkJBA7yN/mDqVD27LlhxMSzN2O79P7doUiI89RgFpCYWFdN3o0YM/YlwcJfEnn5g3YG/ZIvLCC06DTYUKVC+tX2/O8YuRmUn99Msvi9x0E2/pM87gGFS3rsjpp4tcconIffcx9sNwUkhNFbn6anaxU6fSfx/5ze7dXM116xay+97FF/P6utUKpqby4r75pojQSQIQ+f77kE7pFS0IgmHUKA4E2dl+77JtG/3kb789gPMYgWYTJnhscuSISOvWvD/XrPF8KLudXkEAZcup5MmRlsbx0lBPmOruarfzCTT8eBs25GB96JCJJ3HD2rV0TChfnue95hpTBEJmJp0O+vcvqslq0kSkVy9OQnKed4cAACAASURBVO+4Q+Tee2mYNxamru3Kl+fk49VXoyPoMSzYbHxwKlUK2W1v505ey3HjPDQwDNF//SUizlCFEDRRPtGCIBhmzuQl8ncaLnyOK1US2bcvgPMYUYvly3t15dm9m+NT06buVU42G+2TAG2X/kYClzbmzOF1iIuj2j4EOx6ZP5/TZIBT5a++MuGgAXLkCL9MYiJVUDfdJJKSEvBhDhwQeeopOrwBIg0aMF5k3jzfgfJ2O73UWrZ0CoTYWE5q/vsvuK9V6pjocO3+4IOQDzV2LH9Kjz/jo49S0jqWt2PGcMVrpepNC4JgSEvjL/nss341NzzBnn8+iHPt30+raIcOXkfw5ctpfOrVq+iq1WajLzfAALWIGVbDRFoaVe4Ax/CgnLtSU+kjDlC6fvJJ5Ke+R44wQqtCBb5eeskviX7smHNcUYoTkoUL/b8PCgupoahYkVHm779Pz7Xhw50xHrfeegpGDLuyZQu/7BVXhOzhZbPRtHTZZV4adepETw4H111HIWwlWhAEy9lnU1/sB717U48dtEudsVR86imvzQxNkpGPxG6n/zFAX/xTwEvRb6ZP5+y3ShXab/3CZqPnTpUqHDnHjrXQ6BAkKSkczQGRM8+kwdkNhYW8H2rX5mB9550czwJhwwZnjqkrrywZYHjggMjjj3MCUr26yJQpp+A9VlBAQ1z16iVie4LByFLjNqWECF0CY2KKPOtt29ID2Uq0IAiWhx7iNMmHqmDlSl7NF18M8Xx33MEbxKE3dIfdzpQMAO2aQ4fy/RNPnIIPqB/s2kUdN0CB6HVM37uXvv6AyOWXc9obzcyaxallTIzI6NFFVgfbt9PxCKBdM9CI7Lw86q/LlWMGhS+/9H7/bNrk9ODq2ZNu9qcM48bxi331lSmHu+UWBkd6NC/OmcPzzZsnIhToCQl0hLASLQiC5bvveJmWLfPabOBA/vDp6SGe7/hxqilathQ5ccJjs+xsejIahsBHHimbQsAgP5+CEOAizu2s+NtvaRGtWJHT6NJywTIynHq/9u3FvmmzTJ5MW1TVqsFl+VyxgjNQw57kr03cZhN57z2asxo2PEXiCZYsoaA1yYE/PZ3P5ZAhXhqNHk0DjOMZ37GDv8VHH5nSBY9YKggA9AawBcB2ACPdbL8ZwDrH6y8AZ7tsSwGwHsAaT50s/gqrIDhwgJfplVc8Ntmzh7/p44+bdM4FC7jOv+cer82ef15OOric6jYBf5k1i7PbxEQX9928PKcV/fzz6dpVGvnxRzlRq5kMjP3hpFdYoHFpWVmcdcbE8L75+efgurJ6NdNLlCtHe0Kp5dgxuki1aMFJmAlMnsxbbflyL426di0SQjxrFvdZssSULnjEMkEAVhbbAaAFgHgAawEkFWvTBUB1x/srACx32ZYCoFYg5wyrIBChJ4kX5d3TT3PcNjVJ3KhR/Hk8KL7nzaPXzJlnstnbb5t47lJOSgpXS0qJTByxR+wXXOBcNoXbG8hEtm4VaXNGvsSgUF7BY2J76OGAXMP+/NPpETR4cOir16NHaVc1zFqlZYF1ErudvrNxcT5X/IFwwQVMD+PxemRn0+Diogd69VVeR6tjC60UBJ0BzHX5fxSAUV7aVweQ6vJ/9AuCe+/lGtxNcEl+PqN3+/Qx+Zz5+byjEhNLSJgNG2jnbNuWD3OfPlyqW+l/XNrIyhIZ0/NvOYA6khVXRXI+/y7SXQqJOXN4C9asKfL7nAK68wCMWPKRqjI9ncFiAFOh//GHef0qKHBqre6/v5StTD/+WMwx7DnZsIGHfO01L42M4gQuy7Hbb+c4YjVWCoKBAD5y+f9WAJO8tH+sWPudAFYB+BfAYC/7DQawEsDKJk2aWHqxSmBUEHITyTVjBjfNnGnBeXfu5NN/3nknZ7Lp6Uz9W6+e07vjwAFnbqxIez9GDd98I/aEBDlWo4UkqU3Srl1QbvlRwSefUPV49tnFCgZ98QVnlq1aeXQVmjmTKqCYGJH//Y8C0mzsdqpFDXtDqbgH//uPtqLu3U0t/vHQQ/xJDh/20ujZZ7lcdQns6NiRBnirsVIQXOdGEEz00LY7gM0Aarp81sDxt45DrXSRr3OGfUWQkiKe9C8DBjBM37Kb3zBWP/GE2GwM+Y+LE1m82H2zl1+2qB+lBbvdaTzp2lXk8OGTs+l69cJXfMcM7HYGNwP0YHarwl66lLOA6tWLBD4eOsRBGWDeK6/6apMwItpvvz3KVwa5udQd1qxJLzKTyMzkfebT5tyjh0i7dif/tdloXH7kEdO64pGIq4YAtHPYEk73cqyxAB7zdc6wCwIRGpQGDizykZFZcPhwi8/tCBJ48Y4tAniuljdgAF3QykwUaHFsNiaEM3J/u4RobtxIZ6yKFaM7DbeB3S4ybBi/yi23+DAFJCcz/0hCgth/mCFffcXYxHLl6BUZzghzwwtz6NAothk8+ig7+eOPph7W0DQVn6QVIT+fN+HQoSc/2r5dwuIxJGKtIIgDkAyguYuxuE2xNk0cHkVdin1eCUAVl/d/Aejt65wREQS33MKpv8vdbXgHWD7LzM6W+c3ukhgUyo39sz0+YPv20Tuya9con5FZQUEBYzAACgM3F2D/fgZzKuU1rVPEsdmcOn2/o8SPHJHc9udLIWLkTnws551nWQ47r9jtVEEZwY1Rx+zZ7NwDD5h+6HPP9WEkFuHSDGAtCwc//ij+eKibgtXuo30AbHXM+J90fDYEwBDH+48AHHO4iJ50E3V4Gq11vDYa+/p6RUQQGKO+Sy6giy+mQ5HVM58jR0Qa1MmXM9Vmyex+ldeR4ZNP2M1Jk6ztU1SRk8PsakaGLy8/SGamM6PmsGHW14YOlMJCRgcDdBzz596y2RgWUb/KCZkXc7kIILZXvFkrrcVup3+FSSl7zGPXLoaht2sXUCJJfzACSidO9NFw/Hg23L//5EeGJjMjw9QuucVSQRDuV0QEQXKyuNoJDh7kzPLpp609rd3OOr7lyomsGuNIQfHCC17bX345SymUVuNoQJw4QYd6H9lbXSksdDrdFCkhGmFsNuei5pln/BMC27YxhTRAu+f2TXlMWuMj9sVqCgroWhoXFyWlLPPy6HRRpYolNVDvvZcaH58uuT160Gjjwo03Um0ZDrQgMIOWLU/GExgFJ6wulDJlCs8zfrxwZBg0iC4gXp6unTsZdXr55VGspzWDrCwuy2JjWUgkQCZO5KXs3DnytWHsdmfM2zPP+G5fUEAXxQoV6GH84Ycuv3VBAfNMR1gYHD9OVUm1alFgtzIubqDVnfwgPZ3P2913+2iYnU0jXjGrcNu2zPEUDrQgMIP77+cvnpcnffowBYyVA+327ZzZX3KJiwrjxAnmzK9Vy2tYqZFNN4jxsXSQnc3ZVUxMSPlhpk+nwT8pqWSytXAyZoz4nSpk3Trqo40VjVvHl4ICZ2bV8eMt6bM/7NzprOTpKw22ZXzzDa/Do49acvjXX+fhfeZ6mjuXDX/99eRHublcNY0caUnXSqAFgRk4ggaOz14i8fHWunvZbJzsVq3qZrzfvJkS4vzzPSYvt9k4061Z04dPc2kkN5epXpUyRdL98Qc1Bk2aRGbmapSuvvtu70IgN5eqyLg4Dq7ffONDaBQUOH1II+hXvHQp+9y3bwScGIxnpUsXS6LKCws5IezWzY/G//sfZx0uwRyrV/Pn+fZb07vmFi0IzCA9XSQ2VqYN+EoAkUWLrDuVkWraY9F6I3DAi/fD+vV8AAOqlhbt5OVxRPF6cQLn3385uNaqJfLPP6Yd1ieGivG667wbrpcto5rFcCf1W7i7CoMI5iGZMIFdeOmlMJ40M5MXrXZtU+MFXDECSv0qL9muHcvfumA4dwSaOjxYtCAwiy5d5K5aP0n16tYFke3dS71v9+4+ZnyGn97nn3tsYqQsmj/f/H6GncJCWs4BkXfeMf3wW7dydle5Mit1Wc2MGVzU9O7t2dc/M5MrT6VEGjUS+eWXIE5UUOCsbeDlXrESI61PTIy5KS68nvCWW3jhLPwxL7mEK0mfY8G+fW4l4cMP08gcLu81LQhMwv7MWGmAvXL9NdbUk7Pb6d5YvrwfSTILCljhqEIFj5XFs7OZX6ZlS9M95sKLq0/iG29YdprUVBrvypWztoj433/zNz7/fM9pH+bPZ1JMI49PSMkxc3I4G42NDT7laIhkZDBJYp06RbwnreHdd3nh/KwuGAyGWscve7yRpqaYIeHii5lSLFxoQWASaz9dLYDIlPtXWHJ8Q+Pjt7PH/v3MVtWypcdK9b/9JtEb4OMvhjV19GjLT3XsGFXKMTF8fs1m+3aqoFq0cF/c5dgxZiAH+LMGUDLbOxkZjKgrX551LCPAhg08fa9eFtoLFi2iTrRPH0uNEnfeydm8X0bwW2/lj+7SH7udNkCvdQtMRgsCkxj/UqEAIqk3mV9KKDOTy/9zzglQ7bR4sU9r3K23skmpzFA6aRJv1bvuCps/bGYmwxOUMjff/uHDzBFXo4Z7vfCPP1Kux8Sw0I7pq7jDhzktT0wMvKSZSbzzjnhNlRISu3dzyXH66R4nRmZw8CA9Qe+/34/GdjsTXQ0aVOTjnTt5HcJZz0ELApPo3l2kXeJORoCYPCiNHs1fJKjiFIa/qIcZ86FD9CDq0qWUpZ+YNo2jcd++YU9rmZND/26ALoKhkp3N65+QUPI3PnjQ6frfrp3FBus9e6jYrlu3WDrT8GC38+eMj/eo0QyO7Gym4K1SpUgGACt48kneln55ma1bxx92ypQiHxupJcJZ5U0LAhPIyeFD/OhlDuWgib6G27bxwbjlliAP4KpD9+BXb3golJqKUn/+yYty4YXW5E/2gzyXQF0f2Su8YrM5j+OSZkbsdtpva9TgV33uuTDVztm0iZFeZ54ZEQf/Q4c4SU5KMmnVY7cz0aBSlttAMjJ46QYM8HMHo+pMsUCVMWO48gtnZLsWBCbwxx+8YjM/PihmR2327UtvlX37QjhIXh4zzpUv7zYTnt3OFU3VqiGeJxxs28bUyq1bRzzst6CALriObOBBCYNnnpESsV27d1ONDdBguHGjWT32kwULKH0uushjPIqVGPXbTSnxagRjPPecCQfz71R+p/bu1q1I2mmD3r3pmBBOtCAwgaefpgRPTxcq8rt2NeW4RkJEU+TKwYNc9jds6NY1Y+tWrmquu86Ec1nFsWOcqdasaXL9z+Cx2agPBkQefDAw9dq333K/22+nELHZ6NRSpQqNjW+9FcHkd19+yc4NGhQRneHgwXymQlKPzJrFg1x7reU2pNxckQYNOKHyi8OH2bcxY4p8bLfz9r7rLvP76A0tCEygWzc6XYiIUyocOhTSMQsKuDxu1crEvPFr1nCEueAC6rOK8dxz/OWD8km3moICupTExZnoLmMOdjvLzAJMDufP4P3PP1ygXXghB5EtW3gfGcVmkpOt77dPjIoy4cpz4MLx45y3nHFGkCqitWu5lG7fPiw6FqPmwJw5fu4wdSp3KLZCNwzF775rehe9ogVBiGRlcRV9chn777/izgAUKB99JP5HJgbC99/zwDfeWGKml5dH4dOkCVMXRRVGYZnJkyPdE7fY7SJjx7KL11/vXZ+fmsrZY9OmfD9+PIVCtWq8baImIaDd7iyAYGK0tr/Mmycn1W4BsX+/SOPGvMgWRQ67YrNRYJ1zTgC/Xf/+XJ0X22HaNH7ncEaxi2hBEDK//86rNWuW4wO7nTfh1VcHfcysLN7DF1xg0aBgzPTcBBAsXiwnk5xFDUbNh4cfjnRPfGLY//r2dbvokuxsJoarVIkPfYcObN+/f5TaZwoKmK42Qiuxe+/lAtvv4ixZWUwrXbFi2NxgjcH766/93CE7m/1z42P6xBMMWgy3aUYLghAZO5YOCUXyjT/4IKN6g/RoMerRWpazyG53RiZ9/HGJzffdx4cvKur4LljAQahXr1JS/dwZvHr55UVvASNbuFK0xcTF0bX9u++iaBXgjmPHOOWtWZNRb2Hk+HHOq848071gLYLNxrKxSplebtLbKdu0Yf/8tufMnCme9Ejdu7uomcOI1RXKegPY4ihHOdLNdgXgbcf2dQA6+Luvu1ckBMHll7sx/Btr2iCK4B4+zJiefv3M6Z9H8vOdM71iNQyOHaMLX4cOER57d+zg4HPmmX5U9ogupkzheHTxxc4KU4YNpnZt/r3tNlaZKxW4emuF+bcwvIhGjPDRcORINnz11bD0S8Rp8A8o4/k999AjoNi0v7CQH4czotjAyprFsY4SlS1cahYnFWvTB8CvDoFwAYDl/u7r7hVuQVBYyEG7xA+Xl8cNPitSlOSRRzgbD4vLYHo6qyJVrVrihMYNbmH6Hu8cP06DRfXqllSOCgdffcUUPp07O7PGApzhuqSeLz388QcnDr17h32GcNddvJZr1nhoYESZ33df2JZXhYW8RVu3DmA1YLMxYO/660tsMuLLrEhf4gsrBUFnAHNd/h8FYFSxNh8AuNHl/y0A6vuzr7tXuAWB1x/uhhu47g/A/2/fPhoNw5oeetcuTv+bNCliWLPb6ctesWIESlsWFvLkcXGlPj3q9OkcwAwhMGRIeGrQWoYh0YYPD+tpjx7lSuq889w8UtOnc/nVr19YBZRR1+abbwLY6e+/udOXX5bY9P773BRm7ZuIWCsIBgL4yOX/WwFMKtbmFwBdXf6fD6CTP/u6bBsMYCWAlU2aNLH4chXF6w/31VcSqKL/4Yc5aIT9Rli1imvSs84qEk2akkJBcOWVYdZhG2m033svjCc1n7Q0Z3oIgLL2wIFI98oEDA+uMFeg/+ILnnbSJJcPFy9mAEznzmGNMi8s5EogKSnAMIsnnuAEx02+o9tu49wxEvYiKwXBdW4G84nF2sxyIwg6+rOvu1e4VwRef7iMDBqMvRSIcSU1lffznXea20e/mT/fmbbBxXHbKLfnmv7AUoxizA8+GKYTWsMPP1ADANALZOJECtUzzgiLR6O1FBRQPRQXx3QfYcJuF+nZk3OW1FShOrN6dV7UMBtbjOzRAVUQs9k4G/BQiLhlS5aHiARaNRQCrVr58BK94QammPUjScywYVwNRDRg9rvvSiRyKyig0bhePUuTNpLFizlq9uhRajyEirN/Px1XAI5RgLMu+qJFHMRatIiAus1s0tM5Ja5RI6w37fbtVJ9e2yeLg2q9emFPkJeTw1N37BjgamDRIvGkFjroyE4TRjt3EawUBHEAkgE0dzH4tinW5spixuIV/u7r7hVOQZCWxqvktcTeTz+x0ezZXo+1dy9XA+EOK3eL4fvoktp55UoasC31Zti5k0KzVasIVjMPHrudwaLVq/O37NWLl/H554u2W76cgWONG/tRYCjaMTyJkpJCrI4TGC+OPiGAyM8VrqdaM8wYsSIBm6+GDOGy0E20ppFxdOlSc/oYKJYJAh4bfQBsdXgAPen4bAiAIY73CsA7ju3rAXTytq+vVzgFgeEh6rXaXV4eHxQfqUOHDeMqOyrSCogwTQZAXb1DGAwfzo+CSoXti4wM2ieqVYtMlfgQSUlxDvxdulCHrRSDt92pDVevpsyrX9/yrMjWM38+l7JXXhmexEjHjkleu07SRm2QJnVzwh4Bn5bGR7p37wB3zMvj6unGG91ufvRRTiB8xkpYhKWCINyvcAqCF1/kVfI5eb33XoaRejBkHT5MU8Idd5jfx6Cx26mjBygUhJOYJk04+TMt95EIB4++fTmYzJtn4oGtx2aj7r9SJb4mTuQEtVIlRg97y5GzYQO1GrVre3GJLC0YFWUCzgURIJmZlLTlysnSN5YJwAE0nDz+OIV8wPUSfvmF12jmTLebO3RgneNIoQVBkPTvT+OOT/78U7z5mBmT76ibGdpsVA8BIi+/LCLOgMixY008z4gRPOjEiSYe1Hr++492dYCrgZSUogle/UkXsWULK89Vry6ywpoKp+HDSMH66afWHD83l5bimJiTCbiMCPhwFVRLSeGs/bbbgtj5ppu4InAzi0pLo3AZNy70PgaLFgRB0qiRx1VeUWw2jgx9+5bYdOIE740Q0hJZS2EhvyQg8vbbIsL7OTbWpKRYn34q4Q4CCpX8fK4GExI4gE+dyq7n5lIwVKgQWGqO5GSR5s0Zf2iJ2i1c5OczP0J8vPmltfLy6E4DsIqSg7Q0emZ17BgerdTAgfx9d+0KcMfMTNoG7rvP7WbDPhChctEiogVBUOzfL4FF3T72GI0AxVzc3nyTxwlnSbqAyc93PoQffihpaZRrZ54ZYgWppUs5aHTvHqbSW6GzahWzGgMcFIyyDna7s0BNMG62e/awlG6lSgzeLbUcOSJy2mkcnXfvNueY+flcfntYNRpBXRMmmHM6TxjJJZ99NoidjdoOHmKKhg+nJ1QEagCdRAuCIDBUJH7Hiq1eXWRWLcJJTuPGLAIV9eTm0jqmlMhHH8lvv/HrDBsW5PF27WIAxmmnlYpkOzk5IqNGcSVUt67THdTAyCEUisps/34mLytfvpSmnzDYuJE+smbUAXAVAi7Pjit2O1VzVapYF5+Rn0/bWPPmQRpze/SgztCDr+nZZ4tcemlofQwVLQiCYNw4jokBeSycdx7vJocKxKhL4cOzNHrIzqYwcESUDhsmvr2m3HHiBO/8xMQoNIyUZPFiztYBBvsVdw4wtFu33Ra6duvwYea0L1cubMkzrWHWLD4gAwcGX90sP5/Ff/2Y7u/Y4YgtuDa4U/nCWLkH9Zts3Spu/YgdGPEDHjaHDS0IgqB/fw4OAWFEzC5aJHY7a5K2bVtqVOMkJ4dugoBkv/m+tG5Nzxe/0ybYbFQzxcRE/bQ3I8PpONWsmcjcuSXbzJ9Pjd+ll5rnSZWWJnL++TxuQDlsog3D2T6YZVJ+Pkd1gPU6/cDw4vPglBM0Bw5wztKrV5DPqqEW9uA98Pnn7HekU75rQRAELVq4TR7onawsZvm88UaZP59X2E0pgOgnN5fJvQBZN+ILKV+e2az9mviNHs0v/uablnczFObM4UpeKaq/3K381q/nz9mmjfkR18ePs2xlTEwR22jpIljDSU4OHSsCvE/y8vhbNGlibmXKG26gKSuo8JacHKZRHzjQY5Obb6YLcQTKQhdBC4IASU/n1XnhhSB2HjZMpFw5ubXXQalVK3LBIyGTl3dSd/vBgDm+I6xFnAaze+6J2mXQ0aNU8QA0hnuK8kxNpX2nfv0gPEj8JDOTqmWA+Z5KJbm5TAZXoYJ/EcAnTnB5BTA2IUCWLOGujz0WRF/d8PPPEryBWMQ53S9W78PAZmNgoY9407CgBUGAGOlCTpamDISNG0UAeQLj5amnTO9aeMnPF7nhBrEDcv2ZayU21u7Z/XHZMvpbXnSRydFo5mC3M81SnTpcxT/5pGchffQog6ArV7Y+u0FurjNv0YgRUSs/vWPUD27UyOlm5Y60NNZmjY0NKSH/vff6qFvgJ+np9I5r2zaEW/bCC5kyxcN0/59/+Nt+8UXw/TQLLQgC5O23eXVSU4Pbf1vDi2Q7WkjqngivBc2gsFDkvvskHYlyWuJBqV/fXlIVunMnR9gWLUQOHYpEL72yb5/TMaVDBzp4eSIjg/r7+HiPkzzTKSxkihoj/VOpzMW3ahX96C+4wL2EPXCADgTx8UzbGgJG3YLzzw8ttmDIEKrmgg70M4qVvPaaxyaG08nBg0Gew0S0IAiQu+7ijRbM7Oz4cZE7yn8dwpIiCrHbRZ58UtbhLKkYmyOdL7A5/aHT06m4rVpVZPPmiHazOHY77ffVqtHjZPx474NsTg61FrGx4ffosdtFnnmGt02/fiHGb0SK778Xt+5VO3Zw1lyhgnuLfBAYdQvefTe4/Y1kACGlrxgyhKtgL+7R7dszY0Y0oAVBgHTsyEj3YJgwQaQc8iS3TiNaA08lXn9dvsO1NAPcnif2/AKPNZEjTXKyU//erRtTPXgjP/+kfVw+/zw8fXTHxImcQXbrFoaU4FYwbhwv4iuv8P9//uFqsXp1U8Oq7Xb+vomJ/qX6cCUtjZqsVq1CMDofOsTZhZd0wsnJEtG008XRgiAACgoo5P/3v8D3tdtZP+OCC4QSIcDqZaWCqVNltHpRAJFJ507ld/zoo0j36iSFhbz0FStSx//uu769NQoL6dkRpP3SdL7+mnEG7doxIrlUYbfT3U4pLnEqVRJp2tSS1eLWrXxWb7ghsO4NGsS5S0i5n4zlm5c4mTfeYJNIlKV0hxYEAfDff7wyU6cGvu8ff3Dfzz4TupLWri1yxRWm9zHSFM6eK31jZkoMCuXHSy2O+w+ATZvowALwsvvj7VNYKHLrrdznxRet76O/zJvHSNoGDbzbNKKSrCwO/oZrVqBT9gAwIr791cIaTj4hBXdlZdEVyE1uMVe6daMhOlrQgiAApk/nlQkm4dr11zPB3En9rhEBE4HCGpby00+SiUpynlohFZAlf70R2URK+fkcEOLjef0//9w/+05BAWeHAPePNtatowqjUqVSZG6y2USeeooXNT6eXyBYrws/yM1lMH+jRjRXeSM5maqkrl1DTGBn1Hb1ourat4+G6GeeCeE8JqMFQQAYVv5AdYcHDnC5+cgjLh+mp/POu+46U/sYUf79l3qXTp3k4LJkaRm/S2risPz3zFcR6c7KlVShABTE/npn5OU53TbHj7e2j6GQmkqDY0xMdKitvJKZ6UwZcdddLNVWuTJ/IF+jdAisWEEDv7fqfzk59BirWjXEqpdZWbR5XHaZ12aGrPBlmwonlggCADUA/AZgm+NvdTdtGgP4E8BmABsBPOyybSyAVABrHK8+/pzXakFw3XXMkxYoxuS/RHTi6NGULFHmURMUW7dS3dWkycnl/vY1J6R2/DFpiD2ybdBTYYshyM5mjZSYGAZ9zZjh/77Hj/M5LgUB0CLCGCxH1g/53/8iH6Hqll276B4aE8OLaizJ5s3jDKl7d0tTbxoB7e5WTna7s+zGzz+HeCJjhPdh+2vfYyHKwAAAIABJREFUnqnHogmrBMErAEY63o8EMN5Nm/oAOjjeV3GUpUwSpyB4LNDzWi0IkpLoPRIIhYXMVdO9u5uNhw5xVmRVtqxwkZrKL1mrVglpt3ZVodSqcEIaYo9s7TgogMREwbFwIT0+jCDmQLxr9u3jeBUXZ119FSsoKHDmRerfP8BkiFazdClnyYmJ7vNLffYZOz5okGVSLDeXQYANGpS8HyZP5ulDDvBMT2c6CR+rgQ0bxFsy1YhhlSDYAqC+OAf8LX7s8xOAnhKlgiA3lwPE6NGB7Td7Nq/mt996aGC41EV1UQIvHDtGq1flyh6NJ+vWidSqkiMNkCpb6nYzqapNUY4fdxbJat48cI/Vdetow6xUibmGSht2OyfbMTH8OXbsiIIOvfMOXZxatvS+6n35ZQndcd87K1dSReRaEnbZMpoqevUyobCNsezwkT3uiSfYj2gIInPFKkGQXuz/Yz7aNwOwG0CiOAVBCoB1AKa4Uy257DsYwEoAK5s0aWLZhTICBb8KUN199dWcEHnUipw4wST3XbuWvhwC2dnsd7lyPvNRr1snUqtagdSKOSJL4y8R+fBD077vrFk0CCpFO0ygNpxvv6Vpo359S2RUWJkzh0Fy1asHkSLcLDIzWcoOoN7q6FHv7e12kaFDxWrL/JNP8hS//ELjcJ06nDSEXBJjzx4GxPkoWZiTw0XzgAEhns8CghYEAH4HsMHN6+pABAGAygD+BTDA5bO6AGIBxAB4AcAUX/0Ri1cEXzsCgtet83+f/fsp/R9/3EfDDz+UqEk64i/5+XSRU8rv7JJbtoi0bF4gCTF58jVuoIN+CHqMw4edPv5JSYEvqgoKOEMDGOFpoSdjWNm2jQHdMTFUW4d1frF5M3+MmBhmZvRX3WOzOX11LYqyMlREdevS1le9uknmuRtuYABZcrLXZoZ7asQEtBciqhoCUA7AXACPejlWMwAb/DmvlYJgzBje24HYtIyU7D5vNptN5Nxzmdz/+PGQ+hkWXPPFv/deQLsePizSratdAJGnMU4KTk8SWbs2oGPY7czVX7s21XVPPx24rXHrVgb3AVQpRWEuvJDIyHDmUBo0KEy31ZdfUkVYu3Zw0eQFBRxULVSiL1/OuYtSJtUINgKE/Ki70Lkz65hEo0HfKkHwajFj8Stu2igAnwF4y822+i7vHwHwjT/ntVIQDBxII6S/2O2cGHXu7OcOxh364INB9S9sFBTQfSqAoiHFyc2lrhYQubDcctkZf7rI++/7NXVNTXWme+jUKWAZInY7ZVfFipwRluriLz6w2Tgpj4nhvWtZ8Fl6OnMpA8y4GUrNSNca2R98YF4fpWjlS1PSe2dnc2Rv3txnAqi//5ao9kSzShDUBDDf4T46H0ANx+cNAMx2vO8KQBx2gCJuogA+B7Dese1nV8Hg7WWlIDjrLJ/BgkVYtoxXcfLkAE5i1H9csCDg/oWF/HxnlJUJSfK/+EIksYpNEuMyZQruENuVfT2mKrbbqUGrWpWr8NdeCzwT54oVzujiyy+3rsZttLFwIT1mEhL8lrf+s2QJPcZiY5m434z0qLm5In36cGI0ZUrox5Oi85cJEygQypULMZXE44+Lv7qeq65iQGNUeXS5YIkgiNTLKkFQWEjvAp+6fhfuu4/2o4CW5JmZTNfcokX0qYhycpxTcROjrHbuFOnqUBV1VCtlQWLfEtXhd+xw1iu5+GLqwANhxw5nsay6dTm2lDa7fKgcOkTvGIDal+K1lwMmP586uZgY3q9me73l5FBam6AmKix02q6N+UtaGkNemjf3bct2y5Il/O733uuz6Zo1PHfQBW7CgBYEfrB9O6+Iv6Uls7LoNn3rrUGcbPFi3mA33RQ9o1VGBgMhLMq8ZrNxddCoXr4AIn3wi/x66SuSf+iYvPEGBWqVKtQUBKJfXbaMs8CYGAryESOiT76GE5uNwY2xsSy6ErSbrGvI9m23WXdRc3Odupxx44J6HnJzGVUO0EvVlb//5n3RvTvlmt8cPUop4ueErV8/3r8hC18L0YLAD375hVfEU+nC4hj50P/8M8gTGtmyAtIrWcT+/TRkx8Za7tWUlSXy/LhCqVPphAAizbBTkrBeunS2+xX6b7NxjHrySZHWrXkJq1YVGTnS0pQ2pY6VK2m/Arhy9VtdkZ1NaRobS1/bn36ytJ8iQp2OsZwbNiwgh//0dOf8xZMjkhHPNniwn3LGZuPI7meK0t9/5/F9lnKNMFoQ+MFrr/GK+Otv3KMHl5xBewcUFnIdX64cVwiRYu1aJgarWNGE+Hv/yMtjlGeLmGTpiJUnDXtVKtukZ0+6mz//PNP4vvAC2w4ezIpUFSuybWwsVUnvvMPFjKYkOTms7asUJ7Y+Jy0LF9IwCojcfXd4CyLYbAwQMeIS/JiFp6Y6o8R91ZAYNUpO2g58YgQj+NG4oIALp2bNor8+uRYEfnDPPfSI84e9e/lwPf10iCc9dowPXq1akclF9PPPDLNt2JDJ5MLAihU0ygPUjB3aVyCHxr0r0+JvlvviPpIOjQ5IYqL9pHAAqPapWVPkkks4Yfz0U7qoavxj8WIKAoAxGSVs9ampzmCNpk2ZHyhSvPcepfxZZ3nNDrd4MRcslSr5V/TMZqOjUkyMj5CYKVOcgtCP5cPzz7P599/77kOk0YLAD7p29b+gmLF62LrVhBNv3crwx0aNQkyLGAD5+dSlKEX/zDDoVLKymDAtJoZyp8TiY+dOkd69eWFbt5bcn+bIsWPU/0aLGaU0k5XFlVV8PG1bEyaIFGTlUZ9SuTI3PPUUG0aa335j6HSNGiVUU3Y7PZrj4pjVIpDgz8xMPudxcR4Wv99/zxv08sv9CjpZvZoL+kAK40QSLQj8oFYtv5wDRETknHNMziy4Zg1v/MaN3aQvNZmdO51RVvfeG5YH/88/GeVp6Ks9ZiS221ks2Gjcp4/116OMsWWLSM+eXHGdnbBJfselVMUE6qZlNVu3Mm80QF1hTo4cOOC0K19zTXCZrdPTaQ6Li2MmgZN8/jlXIp07+5W/5PBhrrLq1TMhfUWY0ILAB0ePit9R70ZmQb90jYGwejVXBjVrhmCB9kJhocikSZwOJiZ6yZBnHunp1O0DHNv/+MPPHXNzWfO2ShU+sQ88oC3BZrFkidi7XCjTMFAax6UKwPrcYdIMBkZursjw4SKAHGvURvpUWSQJCVyRh7JKPH5c5KKLuCB+8QW72J9/gTdp9+5+GZyMlUX58qUrj6QWBD5YvpxX48cffbcdNYoTB0syLW/dSleY2FhaSc0I3BGhj6Uxu+rRw2e+FDP4+WcGOMXE0GAZ1MLjwAEuIeLi+NQ99hid5TWBs2GDM0akXj2R996TnIx8ef11amAAbo62gW3DBpExnWZLCppQIFx9mykPX3a2yI39cwQQ6Y/pcnDAEL+svUeOcEHt09YQhWhB4APDFXTjRu/tbDa6FvfubXoXnBw/7nSK7tTJf39WdyxZwlBpgJa1b76xXOF+6BATNAJMlRxSVKfB9u30ZY+JoT77ySe1QPCXf/5x6lMSEznBKKb6SE9nGh1DIHTvzkmRWfOQYNi3j6vJmBi6B094MVNsI0ZRKV+pEj2Mdu8O7uAFBSIffij2mrXk9djHJT62QGrVssvEiZ5NA3a7MwNuQoLIDz8E/90ihRYEPnjmGS4TfU0IFi3iVQtLAtFvv+XgbeRKmDnTv6xpaWn0fLjwQu5bsyafcoujrOx25iOrWZPP6rhxFiR527TJmUOgfHmRIUOiT7cdDdjtTGFiRO1Wq0YXNx/K7BMnGJXbsCF3a9SIkbIpKWHqt9Cb+Y47aLuOi6OXWBEPsf/+YxRnbCwb3Hwz3Yb8iRY7coTqUcNFtksXkY0bZcMGqooAeg4+9BDnTEuWiMyfz2vSqRO3JyWZNLmJAFoQ+OCmm+g154vBg+nHHrZcIpmZDBM1BELVqly/jx3LwX7aNBZPeOst6tHPO49TKEMp//bbgSfuD4I9e5hnBaCv/4YNFp9w82b6+8bHU4IPGEApXdbdi7KyRD76iHUSjVwb48cHPAkoKOCMt0cPOenCe8EFvM2sKIZz8KDIu+/Saw/gM/bgg1wIeiQlhSN21arcqXJlLtVHjqQL6rff8vXuu3RX69qVwgPgczJ9epH7xW6n1+zAgZzxu7ovGwJg0iRLq21ajhYEPjj3XN703sjP59L5pptMP71v8vK4Irj7bpEzz+TgV/xOTUzkkzRmDI0eYRgUbTYmOKtShSki3njDhCpQgbB/P6tGVasmhtupvPlmkIllSjHbtnGwq16d1+Gss/jD+MiW6Q87dnAuYmSbAOgtc999jOdYvz5wFVJGBqNxx4zhTNwYn1u3ZnRuQD9fTg71WPffT11kuXIln42EBM5QRo/2yyqelyeyahUXGvPne8yRWOrwJAgUt5UuOnXqJCtXrjTteCJAjRrAoEHAe+95bjdnDnDFFcDPPwN9+5p2+uDIywP27QOys4GYGKB2baBmTUCpsHVh2zbg3nuBhQuBSy8FPvwQaNEibKcvSlYW8O23wOTJwPLlQEICMHAgcPPNQI8eQLlyEeqYhRw9CkybBnzxBfDXX0BcHDBgAPDgg0C3bpbcC1u3AvPm8bVgAXDiBD8vXx5o3hxo3JivatX4WUICb9WMDOD4cWD3bh5j3z7uFxMDtG8P9OrF5++ss0zodkEBcOgQkJ7O/2vUAOrUAWJjQzxw6Ucp9a+IdCrxuRYEwJEjHEdffx149FHP7e66C5g+nfdYQoJppy91FBYCb70FjBkDxMfzut19d1hlkHfWrqVA+PJLjj41agDXXgvccANwySWle0BISwN+/RX47jtg9mwOem3aALfeyleDBmHris3GQX3VKmD1amDnTg70e/ZQQOTmAnY7B/vERL4aNQJOPx1o1YoC4MIL+bkmPHgSBBFX8wTzMls1ZBST8JZmJy+P2ofbbjP11KWOtWudRrN+/aI8139uLqNSb7qJ+mOAur1Bg0SmTi0d6327nfaQ114rqkOpX5+qoNWro9ouUlAQ1d0rc8CDaiguROlSA8C3YJnJFADXi8gxN+1SAJwAYANQKA6J5O/+VrNtG/+2auW5ze+/c6V5/fXh6VO0kZcHvPgiX9WrUwtz3XVRtApwR0IC0K8fXzk5nEH/8gt1fN98wzbnnEM1SufOQJcuQJMmkf1SdjuwaROwaBF1bosWAQcOcFu7dsDIkdRLnnsup9pRTlxII4wmXISkGlJKvQIgTUReVkqNBFBdREa4aZcCoJOIHAlm/+KYrRp6+mnghReobvek8rnjDuDHH6kWio837dSlgmXLqPrZtAm45RaqhWrWjHSvQkCE6qM5c6jsXr6cPz4A1K/PQbZNG76SkoAzzwQqVDC/D0ePAjt2sC9r1wJr1gDr1gGZmWzTsCFw8cXARRcBvXsDTZua2wdNmcMSG4FSaguAS0Rkv1KqPoAFInKGm3YpcC8I/Nq/OGYLgltuAZYsAVJS3G/PywPq1gX69wc++cS000Y9WVm0A7z1FsekDz4A+vSJdK8soLCQA/Dff9Pouno1l4mFhdyuFG+ARo14IRo1AurVA6pU4atyZf6NieGM3mbjq6AAOHaMr7Q0/t23j8r0lBReYIPERODss/nq2JGDf/PmUb7k0pQ2rBIE6SJSzeX/YyJS3U27nQCOARAAH4jI5ED2d2wbDGAwADRp0qTjrl27gu53cbp0oYfDH3+43z5zJrULs2fTa6gsMH8+PYJ27gTuvx94+eUyZtTLz6cw2LSJr927gdRUvvbudXqk+Eu5ctSp1avHAb55c6BZM/5t25bv9aCvsRhPgsCnBk8p9TuAem42PRnA+S8UkX1KqToAflNK/SciiwLYHw7hMRngiiCQfX2RnAxcdZXn7dOm8Rnu0cPMs0Yn6enAY48BH39Mm8nChZycljni453qIXfk51OFc+KE86/dTo+kmBj+jYvjjVO9OlCpkh7oNVGLT0EgIh6HP6XUQaVUfRfVziEPx9jn+HtIKTUDwHkAFgHwa38rycoCDh707P+el8e4gYEDT01XdFd++omz/0OHgBEjgGeeMV81fsoQH0+31Bo1It0TjSZkQnU7+BnA7Y73twP4qXgDpVQlpVQV4z2AywFs8Hd/q9m5k389CYI//mAwzIAB4etTuDl4kC7211zDuJvly6kK0kJAoykbhCoIXgbQUym1DUBPx/9QSjVQSs12tKkLYIlSai2AFQBmicgcb/uHk+Rk/vUkCGbMoC3wssvC16dwIQJ8/jkdY378EXj+eeCff2ir1Gg0ZYeQvHxF5CiAEkOkQxXUx/E+GcDZgewfTrwJApuN6pIrr6Qx+VRi925gyBAGqXbuTJtA69aR7pVGo4kE0R+RYjHJyfT8c+cX/9df1Jf37x/+flmF3Q68+y5toAsXAhMmAIsXayGg0ZRlynzcX3IyVwPuHDpmzKBN8FRxGd26FbjnHg78PXowHU/z5pHulUajiTR6RZAMnHZayc9FKAh69Cj9/vOFhcD48cxQsH49MGUKA2q1ENBoNEAZFwR2O72G3NkH1qxh8Gdp9xZauxY4/3ymqOnTh7FRd96pXdo1Go2TMi0IDh5kqtxmzUpumzGDcUH9+oW9W6aQmws89RTQqRODYb//HvjhB6bS0Wg0GlfKtI3AyFLhSRB07co6BaWNv/5ikrj//gNuvx144w0d96TRaDxTplcERpK54kkdd+4ENmxggFVpIjMTGDaMAiw7m8k1p07VQkCj0XinTAsCY0VQXBDMnMm/ES9HGQDz5rHM36RJrFS4YQPL/2k0Go0vyrQgSEnhbLlKlaKfz5zJFPQtW0akWwFx7BiNv716Meht0SJg4sSS30mj0Wg8UaYFwa5dJe0DGRkMtCoNq4EffmB6iM8/B0aNoqdT166R7pVGoyltlGljcUoKZ/6uzJ3LeiLRLAgOHACGDgWmT2elxdmzWQhco9FogqHMrghE3K8IfvmF6qLOnSPSLa+I0PiblMR+vvgisGKFFgIajSY0yuyK4OhReta4GoptNs6u+/SJvqLbKSnAfffRKHzhhcBHH5VczWg0Gk0wlNkVgeE66roiWLYMOHIkutRCdjuNv2edxfiASZNoENZCQKPR/L+9c42xqjrD8PNqivGKooijgNqKRow3OoLaqlikARKlNvFCqtLUSNSQaE1raUwqxj9qxEYTJYJisFaxxtsgooUplQCKDqjIVUBBBkZQNGDVcP36Y+3jbIdznX32mRn29yQne+112eudtdec76x7tehkv3trR76po9Onh5ZAZ5l2uXJl2CRu/vyg6fHH953q6jiOkxRvEZzU6jd9OlxyCXTv3hGKWtm1K/T/n3122Bto6tRwboAbAcdx0iCRIZDUQ9IsSauj61F54pwm6YPYZ7uk26Ow8ZI2xsJGJNFTCZ99FubaH3lkuF+/PnzpFjvEvhYsXgwDB8Jdd4V9jlasgBtu8E3iHMdJj6QtgnFAo5n1Axqj+x9hZqvM7BwzOwf4OfAd8HIsyt9z4Wb2etv0adHcDH36tN7PnBmuHXX2wPffh7UAAweG6aEvvggvvAC9enWMHsdxskNSQzASmBq5pwKlducZAqw1s/UJ803Mhg3Qu3fr/cyZoZvo1FNrr2XevLAe4L77wq//5cu7/vbXjuN0HZIagl5m1gIQXY8tEf9a4Lk2fmMlLZE0JV/XUg5JYyQ1SWr64osvkqnmxy2CHTugsTG0BmrZBfPNN2Fh2EUXwc6dYWrolClwVMFScBzHqT4lDYGk2ZKW5vmMrCQjSd2AK4AXYt4TgZ8B5wAtwIRC6c1skpnVm1l9z4R7Q+/aFbpfci2CefPg229r2y30xhthSuhjj8Ftt4WTw4YOrV3+juM4OUpOHzWzywqFSdosqc7MWiTVAVuKPGo4sNjMNsee/YNb0mTgtfJkJ6OlJazSzRmCN94IZxP/6lfp5711K9xxBzz9dDgwfv78zrmK2XGc7JC0a6gBGB25RwOvFok7ijbdQpHxyHElsDShnrJobg7XnCGYORMuvhgOPTS9PM3C4G///vDss+H0sPffdyPgOE7Hk9QQ3AcMlbQaGBrdI+l4ST/MAJJ0SBT+Upv0D0j6SNIS4FLgjwn1lMWGDeHau3dwL1sGw4all19LSxj8vfrqMC7R1AT33gsHHZReno7jOOWSaGWxmW0lzARq678JGBG7/w44Ok+865Pk315yLYI+feD554M7jfEBM3jqqdAVtGMH3H9/cHe2fYwcx8k2mfxKam6Gww6DI44I3UJ9+4b++mry6acwZgzMnh26nSZP7pipqY7jOKXI5BYTzc2hW2j37upPG92zBx5+OMwIWrgQJk6EOXPcCDiO03nJZIsgt5hs4cIwl79am8wtXw433hh2MR0+PGwSF1+97DiO0xnJdItg9mw44AAYPDjZ83buDIO/554Lq1fDM8/AjBluBBzH6RpkrkWwe3eYxdOnTzAE9fXJVvI2NYVWwJIlcM018MgjcGyp9dWO4zidiMy1CD7/PBz2cvTRoQvnsoLL5Yrz/fdw550waFA4zOaVV2DaNDcCjuN0PTLXIshNHf3qqzCw2x5D8NZb4cCYNWvgppvggQdat7N2HMfpamSuRZBbTLZmDRx8MFx4Yflpt2+HW24JYwp794YZR5MmuRFwHKdrkzlDkGsRLFoU5veXu7p3xgw444zwxX/HHWFMoBZ7EzmO46RNJg3BwQfDqlXldQt9+SVcd104ueyII8IB8hMmpLsvkeM4Ti3JpCHInUlczBCYhcHf008P21DcfXc4RnLQoNrodBzHqRWZGyzOjRH07AlnnZU/zsaNcOut0NAA550HTz4JZ55ZO42O4zi1JJMtgm3bYMiQsJgsjlnYE6h/f5g1Cx58EN5+242A4zj7N5lqEezZA5s25Z82unZtmAo6Z06YFTR5MpxySofIdBzHqSmZahFs3hyMALQagj174KGHwq/+RYvC/kCNjW4EHMfJDplqEeTGB+rq4MQTYenSsD3Eu++GWUETJ7aeWuY4jpMVErUIJF0laZmkvZLqi8QbJmmVpDWSxsX8e0iaJWl1dE2w609p1q0L1wsugHvugQED4JNPwtGRDQ1uBBzHySZJu4aWAr8F5haKIOlA4FHC4fX9gVGS+kfB44BGM+sHNEb3qbFgQbi+8w6MHw9XXRW2jh41qnrnETiO43Q1kh5VuQJAxb9FBwJrzOyTKO40YCSwPLoOjuJNBf4L/CWJpmK8/HK47t0L06eH7iDHcZysU4sxghOADbH7ZiC3LKuXmbUAmFmLpIJ7d0oaA4wB6Nu3b7uEnHRS2C9o5crWRWWO4zhZp6QhkDQbOC5P0F1m9moZeeRrLlgZ6X6cwGwSMAmgvr6+4vQAcwt2YDmO42SXkobAzNq5Y/8PNAPxs7p6A5si92ZJdVFroA7YkjAvx3Ecp0JqsY7gPaCfpJMldQOuBRqisAZgdOQeDZTTwnAcx3GqSNLpo1dKagYuAGZIejPyP17S6wBmthsYC7wJrAD+ZWbLokfcBwyVtBoYGt07juM4NURm7epu71Dq6+utqampo2U4juN0KSQtMrN91nxlaosJx3EcZ1/cEDiO42QcNwSO4zgZxw2B4zhOxumSg8WSvgDWtzP5McCXVZRTLVxXZbiuynBdldFZdUEybSeaWc+2nl3SECRBUlO+UfOOxnVVhuuqDNdVGZ1VF6SjzbuGHMdxMo4bAsdxnIyTRUMwqaMFFMB1VYbrqgzXVRmdVRekoC1zYwSO4zjOj8lii8BxHMeJ4YbAcRwn4+yXhkDSVZKWSdorqeA0K0nDJK2StEbSuJh/D0mzJK2OrkdVSVfJ50o6TdIHsc92SbdHYeMlbYyFjaiVrijeOkkfRXk3VZo+DV2S+kiaI2lF9M5vi4VVtbwK1ZdYuCQ9EoUvkTSg3LQp6/pdpGeJpAWSzo6F5X2nNdI1WNK22Pv5W7lpU9b155impZL2SOoRhaVSXpKmSNoiaWmB8HTrlpntdx/gdOA0whnI9QXiHAisBX4KdAM+BPpHYQ8A4yL3OOD+Kumq6LmRxs8Ji0AAxgN/SqG8ytIFrAOOSfp3VVMXUAcMiNyHAx/H3mPVyqtYfYnFGQHMJJzKdz6wsNy0Keu6EDgqcg/P6Sr2TmukazDwWnvSpqmrTfzLgf/UoLwuBgYASwuEp1q39ssWgZmtMLNVJaINBNaY2SdmthOYBoyMwkYCUyP3VOA3VZJW6XOHAGvNrL2rqMsl6d/bYeVlZi1mtjhyf0M48+KEKuUfp1h9iet92gLvAEcqnLxXTtrUdJnZAjP7Orp9h3BKYNok+Zs7tLzaMAp4rkp5F8TM5gJfFYmSat3aLw1BmZwAbIjdN9P6BdLLzFogfNEAx1Ypz0qfey37VsKxUdNwSrW6YCrQZcC/JS2SNKYd6dPSBYCkk4BzgYUx72qVV7H6UipOOWnT1BXnRsIvyxyF3mmtdF0g6UNJMyWdUWHaNHUh6RBgGPBizDut8ipFqnWr5JnFnRVJs4Hj8gTdZWblHHmpPH6J59IW01Xhc7oBVwB/jXlPBO4l6LwXmAD8oYa6fmFmmyQdC8yStDL6JdNuqlhehxH+YW83s+2Rd7vLK18Wefza1pdCcVKpayXy3DeidCnBEPwy5l31d1qBrsWEbs//ReM3rwD9ykybpq4clwPzzSz+Sz2t8ipFqnWryxoCM7ss4SOagT6x+97Apsi9WVKdmbVEza8t1dAlqZLnDgcWm9nm2LN/cEuaDLxWS11mtim6bpH0MqFZOpcOLi9JPyEYgX+a2UuxZ7e7vPJQrL6UitOtjLRp6kLSWcATwHAz25rzL/JOU9cVM9iY2euSHpN0TDlp09QVY58WeYrlVYpU61aWu4beA/pJOjn69X0t0BCFNQCjI/dooJwWRjlU8tx9+iajL8McVwJ5ZxikoUvSoZIOz7mBX8fy77DykiRaklpiAAABKklEQVTgSWCFmT3UJqya5VWsvsT13hDN8Dgf2BZ1aZWTNjVdkvoCLwHXm9nHMf9i77QWuo6L3h+SBhK+j7aWkzZNXZGe7sAlxOpcyuVVinTrVrVHvzvDh/BP3wzsADYDb0b+xwOvx+KNIMwyWUvoUsr5Hw00Aquja48q6cr73Dy6DiH8Q3Rvk/4fwEfAkuhl19VKF2FWwofRZ1lnKS9CN4dFZfJB9BmRRnnlqy/AzcDNkVvAo1H4R8RmrBWqa1Uqp1K6ngC+jpVPU6l3WiNdY6N8PyQMYl/YGcoruv89MK1NutTKi/CjrwXYRfjuurGWdcu3mHAcx8k4We4achzHcXBD4DiOk3ncEDiO42QcNwSO4zgZxw2B4zhOxnFD4DiOk3HcEDiO42Sc/wOnTf86kkyX9QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "h=.01\n", "T = np.arange(-1,1+h,h)\n", "plt.figure()\n", "plt.title('Legnendre Polynomials')\n", "\n", "\n", "\n", "\n", "for n in range(0,3):\n", " Xe = np.zeros(T.size)\n", " Xo = np.zeros(T.size)\n", " ae = 2*n\n", " ao = 2*n + 1\n", " for k in range(0,ao+1):\n", " Xe = Xe + A(k, ae, 1, 0)*(T**k)\n", " Xo = Xo + A(k, ao, 0, 1)*(T**k)\n", " Xe = Xe/Xe[-1]\n", " Xo = Xo/Xo[-1]\n", " plt.plot(T,Xe, 'r-')\n", " plt.plot(T,Xo, 'b-')" ] }, { "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }