{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "2fn71JmtwN_n", "pycharm": {} }, "source": [ "# Response of a First Order System to Step and Square Wave Inputs\n", "\n", "Keywords: Simulator, ipopt usage\n", "\n", "```{index} Simulator\n", "```\n", "```{index} ipopt\n", "```\n", "\n", "This notebook demonstrates simulation of a linear first-order system in Pyomo using the `Simulator` class from Pyomo. The notebook also demonstrates the construction and use of analytical approximations to step functions and square wave inputs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": {}, "colab_type": "code", "id": "I-G6nxedwUAF", "pycharm": {} }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "from math import pi\n", "\n", "import shutil\n", "import sys\n", "import os.path\n", "\n", "if not shutil.which(\"pyomo\"):\n", " !pip install -q pyomo\n", " assert(shutil.which(\"pyomo\"))\n", "\n", "if not (shutil.which(\"ipopt\") or os.path.isfile(\"ipopt\")):\n", " if \"google.colab\" in sys.modules:\n", " !wget -N -q \"https://ampl.com/dl/open/ipopt/ipopt-linux64.zip\"\n", " !unzip -o -q ipopt-linux64\n", " else:\n", " try:\n", " !conda install -c conda-forge ipopt \n", " except:\n", " pass\n", "\n", "assert(shutil.which(\"ipopt\") or os.path.isfile(\"ipopt\"))\n", " \n", "from pyomo.environ import *\n", "from pyomo.dae import *" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "DPakBPtBwN_q", "pycharm": {} }, "source": [ "## First-order differential equation with constant input\n", "\n", "The following cell simulates the response of a first-order linear model in the form\n", "\n", "$$\n", "\\begin{align}\n", "\\tau\\frac{dy}{dt} + y & = K u(t) \\\\\n", "\\end{align}\n", "$$\n", "\n", "where $\\tau$ and $K$ are model parameters, and $u(t)=1$ is an external process input. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "colab_type": "code", "executionInfo": { "elapsed": 5593, "status": "ok", "timestamp": 1557657354444, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg", "userId": "09038942003589296665" }, "user_tz": 240 }, "id": "tHGcynz1wN_r", "outputId": "40db578c-9045-4e62-ba7f-0060cc699a97", "pycharm": {} }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApmklEQVR4nO3deXxcdb3/8dcne5uke5uudGdpQSgtFASlZVG4IIuKoOAVAdGfckGvVwW8CN57vXLV65XrdUdERaiyKUsFEVtWy1KgdIWW0iVN26RLmqRt1vn8/jgnMA1JO0kzcyZn3s/HYx6ZOevnO5m8c+Y7Z77H3B0REYmfvKgLEBGR9FDAi4jElAJeRCSmFPAiIjGlgBcRiSkFvIhITCngJauY2UlmttrMGszs/IPc1s1mdmd4/5Bwm/m9UuhBMrP/MLNtZrYl22rrjJktNLMro65DukcBHyEzW2dme8M/7i1mdoeZlUVdV8T+Dfg/dy9z9z/21kbdfUO4zbbe2mZPmdk44MvANHcfeTC1mdllZvZM71eZXmZWbGbfNrMN4d/AajP7iplZ0jILzazRzOrNrM7MFpvZdWZWnLTMzWbWEv4Ntd9qI2lUFlLAR+9D7l4GHAPMAK6PtpzIjQeWR11EbzGzgk4mjwe2u3t1CuubmWXs77S399dF+wHuAU4D/gEoBz4JXAXc2mG5q929HBhF8E/xYmB+8j8C4PfhP8j226Deqr+vU8BnCXffAjxGEPQAmNkJZvacmdWa2RIzm5M07zIzWxse3bxlZpckTX/WzH5oZrvMbJWZnZa03mgze9DMdpjZGjP7TNK8m83sD2b2m3C7y81sVtL8r5nZpnDe6+3bNbO88MjqTTPbHm5jSFdtNbPPhPveEdYyOpz+JjAJeCg8EivuZN32/dSb2QozuyCV59fMJpiZtwdOeHT47+FzVW9mfzGzYSk+9582s5XhemvN7LNJ8+aYWWX4XG0BftWhjtOBx4HRYRvv6KK2b5nZs8AeYFJnv28zOwL4KXDigY5czey9ZvZi+Jp40czemzSvs/2dEb52dpnZ/wHWYXuXh8/BTjN7zMzGJ81zM/uCma0GVndSy2nAB4CPuPsyd29190XApcAXzGxKx3Xcfbe7LwTOBU4Ezu6qrZLE3XWL6AasA04P748FlgK3ho/HANsJjnDygDPCx8OBUqAOOCxcdhQwPbx/GdAKfAkoBC4CdgFDwvlPAj8GSgj+mdQAp4XzbgYaw33mA98GFoXzDgM2AqPDxxOAyeH9LwKLwjYUAz8D7u6izacC24Bjw2V/CDzV2XPSxfoXAqPD5+QiYDcwqotlbwbuTKrXgYLw8ULgTeBQoF/4+JYDPffh/LOByQShdwpBKB4bzpsTPv//FbavXyd1zQEqkx53VtsGYDpQAAw8wO/7mQO8zoYAOwmOkguAj4ePh3axv+Hh/j5K8Br6UtimK8PlzwfWAEeEy/8r8FzS/pzgn9iQLtp/C/BkF7WuBz6bVNeVnSzzFPBfHX/Hur37piP46P3RzOoJwrMauCmcfikw393nu3vC3R8HXiIIHYAEcKSZ9XP3ze6e3K1RDfzA3Vvc/ffA68DZFvT9ngx8zd0b3f1V4DaCP/x2z4T7bAN+CxwdTm8jCKxpZlbo7uvc/c1w3meBr7t7pbs3EfzRfbSLt+eXALe7+8vhstcTHIFOSOXJcvd73L0qfE5+T3CEeHwq63biV+7+hrvvBf7AO++e9vvcu/sj7v6mB54E/gK8L2m7CeAmd28Kt90Td7j7cndvJQjX/f2+D+RsYLW7/9aDo+W7gVXAh7rY31nACne/191bgB8AW5KW/SzwbXdfGS7/n8AxyUfx4fwdXbR/GLC5i1o3h/P3p4rgn0e7j4XvtNpvCw6wfs5QwEfvfA/6GOcAh/POi3s8cGHyC5cgnEe5+26Co9fPAZvN7BEzOzxpm5s8PLwJrSc46h0N7HD3+g7zxiQ9Tv5D3gOUmFmBu68hOFK/Gag2s3ntXSthrQ8k1bmS4B9CRSftHR3uEwB3byA4Oh7TybLvYmb/aGavJu3rSA4cCF3p2Nb2D7i7fO7DGs4ys0VhF1MtQfAn11Dj7o09rKndxvY7Kfy+32bvnJHTYGYN4eR9nvNQx9/7xqT7ozvs3zvMHw/cmvTc7CB4N9PV9jraRvhcdmJUOH9/xoT7bPcHdx+UdJt7gPVzhgI+S4RHgncA3wsnbQR+2+GFW+rut4TLP+buZxD8QawCfpG0uTFm+3wIdQjBUU8VMMTMyjvM25RijXe5+8kEf+BO0A3RXutZHWotcffOtlsVrg+AmZUCQ1OpITxC/AVwNUH3wiBgGR36h3tBl899+LnAfQS/p4qwhvkdauiNIVr32cZ+ft8dl2s/I6fMgw/vocNzHur4e0/ezmZgXPuD8LU0Lmn+RoJulOTnp5+7P9dV/R38FZgdvqN8m5kdH+7nb12tGK4zE3h6P9uXkAI+u/wAOMPMjgHuBD5kZh80s3wzKwk/wBtrZhVmdm4Yjk1AA8ERc7sRwDVmVmhmFxL0lc53943Ac8C3w+29B7gC+N2BCjOzw8zs1DDgGoG9Sfv8KfCt9rfoZjbczM7rYlN3AZ82s2PCbf0n8Ly7r0vh+SklCI6acD+fJjiC721dPvdAEUFXVQ3QamZnEXxgmDYH+H1vBcaaWdF+NjEfONTMPmFmBWZ2ETANeLiL5R8BppvZh8NutmuAkUnzfwpcb2bTw/oGhq+zlLj7X4EngPvMbHr4HJ9A8Dr8ibt39sFsfzM7BfgT8ELYJjkABXwWcfca4DfAjWEYnwfcQBAmG4GvEPzO8ghOGasieKt6CvD5pE09D0wleKv7LeCj7r49nPdxgg/1qoAHCPqKH0+hvGKCD8e2EXRtjAhrg+DUtgeBv4SfJywCZnfRxieAGwmOgjcTfFh5cQr7x91XAP8N/J0g2I4Cnk1l3e7Y33Mfdm9dQ9BnvxP4BEHb02l/v++/EZxWusXMOu3aCH/354Tb2A58FTjH3btafhvBh9m3hMtPJel5dvcHCN69zTOzOoJ3UWd1s00fARYAjxL8w7oT+CXwTx2W+7/wNbWV4ADoPuBMd08kLXOR7XsefIOZjehmPbFk+3bVSl9nZpcRnHlwctS1iEi0dAQvIhJTCngRkZhSF42ISEzpCF5EJKa6GggoEsOGDfMJEyb0aN3du3dTWlrauwVlObU5/nKtvaA2d9fixYu3ufvwzuZlVcBPmDCBl156qUfrLly4kDlz5vRuQVlObY6/XGsvqM3dZWYdv6X8NnXRiIjElAJeRCSmFPAiIjGVVX3wnWlpaaGyspLGxv0Pzjdw4EBWrlyZoar2VVJSwtixYyksLIxk/yIincn6gK+srKS8vJwJEyaw7wCJ+6qvr6e8vLzL+eni7mzfvp3KykomTpyY8f2LiHQlrQFvZuuAeoKR71rdfdb+13i3xsbGA4Z7lMyMoUOHUlNTE3UpIiL7yMQR/NyuRq1LVbaGe7tsr09EclPWd9GI5Cp3py3hNLW2kUhAayJBIgFt7vvcb2vz4GfinVvCndak++0/cUg4OB78dMdpvzYznU5zIPH2fA9re2da+3zapyXPT95O0rAovk872Wf6mnUtvPnMW3Q1jMq+y3sX07u3fMfnvTe22eFKLJ3vLLSlspl0nPqf1rFozOwtgjGzHfiZu/+8k2WuAq4CqKiomDlv3rx95g8cOJApU951kfV3aWtrIz8/vzfK7pE1a9awa9eujO6zoaGBsrKyAy8YI5lqc1vC2dsKTW1OUxu0JJzmNmhqg+Y2pzkR/gwfB8sEy7cmoPXtIA7ut3qwzfb7wXQPl3vnfrBcewhLnOzvfX55ofO/p/XsdT137tzFXXV/p/sI/iR3rwoH33/czFa5+1PJC4Sh/3OAWbNmecdvc61cuTKlD0+j+pC1XUlJCTNmzMjoPvWNv665O7ub29i5u5kdu5vZsaeZ2j3N7NjdQu2eZuobW2loaqW+sYWGplYaGlupb2ylPry/t6XtgPvoqLggj5LCfIoL8ijMz6OoII/CfKOwMHjcPz+PwgKjMD+cnx/Of3vZ4HFBfh75ZuTnGRvWr2PK5Enk5xn5ZuTlGQV5wc98S7qfB/l57euF9/Mgz4yCvDzyLOhK3PcnQIdpGBbOa7+fZ+0/g+XfngZv3wfIy9t3WnA/3F77/aTnK7lnM3nOM88+w8knn9zJMsnrWhfTO99mV72oqSzf7f32oMs2XX/LaQ14d68Kf1ab2QPA8cBT+18ru9x4440MGzaMa6+9FoCvf/3rVFRUcM0110RcWW5qTThVtXvZWtfI1rqm8Gdwv7q+kZr6JnbsbmbnnmZa2jo/BjaDsuICBpQUUlZcQFlJAYP6FzFuSH/KSwooKy6gPJzXvyiffkX59CvMp6Rw3/slhXn0C6eVFOSTl9f7n8UsXFjFnDkHfgcbJ6WFxsB+OuW4N6Qt4MPrR+a5e314/wPAvx3MNr/50HJWVNV1Oq+nXTTTRg/gpg9N73L+FVdcwYc//GGuvfZaEokE8+bN44UXXuj2fiQ17s7WuiY27NjDxh17gp87g/sbdwTB7n/Z95rM+XnGiPJiKgaUMHZwf44eO4jBpUUMKS1kcP+i4FZaxJDSIob0L6K8pCAtYSySbdJ5BF8BPBC+XSkA7nL3R9O4v7SYMGECQ4cO5ZVXXmHr1q3MmDGDoUOHRl1Wn9eWcDbu2MPq6gbWVDewurqeNeH9Pc3vdI+YwcgBJYwb3J+TpgyjdddWZr/ncEYOLGZEeQkVA0oYWlqkwBbpRNoC3t3XAkf35jb3d6Sdzj74K6+8kjvuuIMtW7Zw+eWXp2UfcdbalmBNTQNLK3exbNMullXVsaKqbp9+7pEDSpgyooyPzRrH5OGlHDK0lHGD+zFmcD+KC955Z7Zw4ULmzD4kimaI9Dk6TTIFF1xwAd/4xjdoaWnhrrvuirqcrFff2MLi9Tt5cd0OXnhrB0s37aKxJQFA/6J8po8ewMXHj+OIkQOYUlHGlBFlDChRn6tIb1PAp6CoqIi5c+cyaNCgSE/FzFZ7m9tY9NZ2nnqjhhfe2sHKzXUkPOgbP3LMQD5x/HiOGjuAo8YMZOKwMvLVnSKSEQr4FCQSCRYtWsQ999wTdSlZY21NAwter+HJN2p4fu12mloTFBfkcewhg7n61KnMnjiEGYcMon+RXmIiUdFf3wGsWLGCc845hwsuuICpU6dGXU6k1lTX88hrW5i/dDOvb60HYPLwUi6ZPZ45hw3n+IlDKCnUOxyRbKGAP4Bp06axdu3aqMuIzOZde7n/5U386dVNvLG1ATM4bvwQbvrQNE4/ooJxQ/pHXaKIdEEBL+/S1NrGX1dU84eXNvL06hoSDsdPGMI3z53OmUeOpGJASdQlikgKFPDytuq6Ru5ctJ7fPb+B7bubGTWwhC/MncJHZ45l/NDcusq9SBwo4IXlVbv45TNv8dCSKloTzmmHj+DSE8bzvqnDdcaLSB+mgM9hyzbt4tYnVvP4iq30L8rnktnj+dR7JzBxmI7WReJAAZ+Dllft4gd/DYJ9QEkBXzr9UC47aYIGeBKJGQV8Dtla18h3H3ud+16upLy4gH8+Iwh2fYtUJJ4U8Aewbt06zjnnHJYtWwbA9773PRoaGrj55pujLawbGlva+PlTa/npk2/S2uZc9b5JfH7uFB2xi8Rc3wr4P18HW5Z2OqtfWyvk96A5I4+Cs245yMKy17NrtnHDA0tZv30PZx05kuvOOlxnxIjkiL4V8JKynbub+db8ldy7uJIJQ/tz15Wzee+UYVGXJSIZ1LcCfj9H2nvTNFxwQUEBiUTi7ceNjY29vo/e9rdVW/nqva9Ru6eFz8+ZzDWnTdUQAiI5KC/qArJdRUUF1dXVbN++naamJh5++OGoS+pSY0sb3/jTMi6/4yWGlRXz4NUn89UzD1e4i+SovnUEH4HCwkK+8Y1vMHv2bCZOnMjhhx8edUmdWr21ni/c9TJvbG3gipMn8tUzD9vnQhkiknsU8Cm45pprsvoi239eupkv37OE/kX5/Pry4znl0OFRlyQiWUAB34cl3PnuY6v40YI3OWbcIH566UxGDtRAYCISUMD3UbubWrn15SaW1LzJxceN45vnTVeXjIjso08EvLtjlr2DXrl7Rve3raGJy+94kaU1bfz7edP55IkTMrp/Eekbsv4smpKSErZv357xEE2Vu7N9+3ZKSjLTNbJu224+8pPneGNrPdccW6xwF5EuZf0R/NixY6msrKSmpma/yzU2NmYsZDsqKSlh7Nixad/Piqo6PvnL50m4c9dnTqBu7ZK071NE+q6sD/jCwkImTpx4wOUWLlzIjBkzMlBRNJZX7eKS256nX2E+d145m8nDy1iYu1cSFJEUZH3ASzBu+yW3PU9pUT53X3WCxpIRkZQo4LPcys11XHLb85QVFzDvqhN0kWsRSVnWf8iayzZs38M/3v4C/YvyFe4i0m0K+CxVU9/EJ29/npa2BL+5/HiFu4h0mwI+C9U3tvCp21+guq6J2y87jqkVvT9KpojEn/rgs0xbwvnivFd5fWs9v/zULI49ZHDUJYlIH6Uj+CzzncdW8cSqam7+0DTmHDYi6nJEpA9Le8CbWb6ZvWJm2TuQepa4b3ElP3tyLZeecIi+oSoiBy0TR/DXAiszsJ8+7dWNtVx//1JOnDSUmz40PepyRCQG0hrwZjYWOBu4LZ376et27Wnh6rteZnh5MT++5FgK89VzJiIHz9I5iJeZ3Qt8GygH/sXdz+lkmauAqwAqKipmzps3r0f7amhooKys7CCqjYa788NXmlhS08YNs0uYPCj1IX/7apsPRq61OdfaC2pzd82dO3exu8/qdKa7p+UGnAP8OLw/B3j4QOvMnDnTe2rBggU9XjdKtz+z1sd/7WH/xVNvdnvdvtrmg5Frbc619rqrzd0FvORdZGo6+wJOAs41s3XAPOBUM7szjfvrc5Zt2sV/zl/J6UeM4IqTDzygmohId6Qt4N39encf6+4TgIuBv7n7penaX1/T2NLGl37/KkNKi/juR4/O6guaiEjfpC86ReT7j7/B6uoG7vj0cQwuLYq6HBGJoYwEvLsvBBZmYl99wYvrdvCLp9fyidmH6MtMIpI2Oh8vw/Y2t/Ev9yxh7OB+3PAPR0RdjojEmLpoMuzWJ1azfvse7v7MCZQV6+kXkfTREXwGrdpSx21Pr+XCmWM5cfLQqMsRkZhTwGdIIuHccP9SyksKuF5dMyKSAQr4DJn34kZe3lDL18+exhCdNSMiGaCAz4DaPc1857FVzJ44hI8cOybqckQkRyjgM+DWJ1ZTt7eFm8+dri80iUjGKODT7M2aBn779/VcdNwhHDFqQNTliEgOUcCn2bfnr6SkMJ8vf+DQqEsRkRyjgE+jZ1Zv468rq7n61CkMKyuOuhwRyTEK+DRxd255dCVjB/fj0ydNiLocEclBCvg0eWz5FpZtquNLpx9KcUHqF/EQEektCvg0aEs4//2XN5g8vJTzZ+i0SBGJhgI+DR5csonV1Q388xmHkZ+n0yJFJBoK+F7W0pbgfx5fzbRRAzjryJFRlyMiOUwB38seeHkTG3bs4V8+eCh5OnoXkQgp4HtRW8L5yZNvcuSYAczVhTxEJGIK+F706LItvLVtN5+fM0VDEohI5BTwvcTd+dGCNUwaXsoHp6vvXUSip4DvJU++UcOKzXV87pTJOnNGRLKCAr6X/HjBm4weWML5x+i8dxHJDgr4XvDqxlpeWLeDK983iaICPaUikh2URr3gV8++RXlxAR87blzUpYiIvE0Bf5C27Grkkdc2c+GscZQVF0RdjojI2xTwB+nORetpc+ey906IuhQRkX0o4A9CY0sbd72wgdMOr+CQof2jLkdEZB8K+IPw4KtV7NjdzOUa711EspACvofcnTueW8dhFeWcOHlo1OWIiLyLAr6HXqvcxYrNdVx64ngNSyAiWUkB30PzXtxAv8J8zjtmdNSliIh0Km0Bb2YlZvaCmS0xs+Vm9s107SvTGppaefDVKs55zygGlBRGXY6ISKfSeeJ2E3CquzeYWSHwjJn92d0XpXGfGfHQkip2N7dx8fGHRF2KiEiX0hbw7u5AQ/iwMLx5uvaXSfNe2MBhFeUce8igqEsREemSBTmcpo2b5QOLgSnAj9z9a50scxVwFUBFRcXMefPm9WhfDQ0NlJWVHUS1qVlf18ZNzzVyyeFFnDEh2u6ZTLU5m+Ram3OtvaA2d9fcuXMXu/usTme6+wFvQH/gRuAX4eOpwDmprBsuPwhYABy5v+VmzpzpPbVgwYIer9sd3/jjUp/69fm+c3dTRva3P5lqczbJtTbnWnvd1ebuAl7yLjI11Q9Zf0XQp35i+LgS+I9U/8O4ey2wEDgz1XWyUUtbgode28wZ0yoY1L8o6nJERPYr1YCf7O7fAVoA3H0vsN+Tv81suJkNCu/3A04HVvW81Og9+XoNO3Y38+EZGvNdRLJfqh+yNoch7QBmNpngiH5/RgG/Dvvh84A/uPvDPa40CzzwyiaGlBbx/kOHR12KiMgBpRrwNwGPAuPM7HfAScBl+1vB3V8DZhxUdVmkrrGFx1du5ePHjaMwX98PE5Hsl1LAu/vjZvYycAJB18y17r4trZVlmT8v3Uxza4ILjh0bdSkiIilJ6VDUzE4CGt39EYIzYm4ws/HpLCzb3P/yJiYNK+XosQOjLkVEJCWp9jX8BNhjZkcDXwHWA79JW1VZZlPtXp5/awfnzxijgcVEpM9INeBbw/MtzwP+191vBcrTV1Z2eeS1KgDOP0Znz4hI35Hqh6z1ZnY9cCnw/vDMmJwZZeuRpVs4asxAXbVJRPqUVI/gLyI4LfIKd98CjAG+m7aqskjlzj0s2VjLWUeNjLoUEZFuSfUsmi3A95MebyBH+uAfXbYFgLOPGhVxJSIi3ZPqWTQfNrPVZrbLzOrMrN7M6tJdXDZ4ZOlmpo8ewPihpVGXIiLSLal20XwHONfdB7r7AHcvd/cB6SwsG1TV7uWVDbX8g47eRaQPSjXgt7r7yrRWkoX+HHbPKOBFpC9K9Syal8zs98AfSRqDxt3vT0dR2WL+0s0cMWoAE4epe0ZE+p5UA34AsAf4QNI0B2Ib8NV1jSxev5Mvn3Fo1KWIiPRIqmfRfDrdhWSbJ1ZVA3DG9IqIKxER6ZlUz6IZa2YPmFm1mW01s/vMLNajbj2xcitjBvXjsIqc+cKuiMRMd67o9CAwmuBLTg+F02Jpb3MbT6/exhnTKjT2jIj0WakG/HB3/5W7t4a3O4DYXvXi2TXbaGpNcNoRI6IuRUSkx1IN+G1mdqmZ5Ye3S4Ht6SwsSn9duZWy4gJmTxwadSkiIj2WasBfDnwM2BLePhpOi51EwnliVTWnHDqcogJduUlE+q5Uz6LZAJyb5lqywmubdlFT38Tp09Q9IyJ9W6pn0Uwys4fMrCY8k+ZPZjYp3cVF4YmVW8kzmHOoAl5E+rZU+yDuAv4AjCI4k+Ye4O50FRWlJ1ZWM2v8EAaXFkVdiojIQUk14M3df5t0Fs2dBN9kjZWa+iZWbK7jlMNie4KQiOSQVIcqWGBm1wHzCIL9IuARMxsC4O470lRfRj29ugaA909VwItI35dqwF8U/vxsh+mXEwR+LPrjn169jaGlRUwfHfuRkEUkB6R6Fs3EdBcStUTCeXp1DSdPHUZenr69KiJ9X6pn0VxoZuXh/X81s/vNbEZ6S8usFZvr2NbQrO4ZEYmNVD9kvdHd683sZOCDwK+Bn6avrMx7Kux/f9/UYRFXIiLSO1IN+Lbw59nAT9z9T0CsziN8+o1tHD6ynBEDSqIuRUSkV6Qa8JvM7GcEwxXMN7Pibqyb9XY3tfLS+h2ccqi6Z0QkPlIN6Y8BjwFnunstMAT4SrqKyrRFa7fT0ua8XwEvIjGSUsC7+x6gGjg5nNQKrE5XUZn29OptlBTmMWvC4KhLERHpNameRXMT8DXg+nBSIXDnAdYZZ2YLzGylmS03s2sPrtT0WbR2O7PGD6G4ID/qUkREek2qXTQXEIwmuRvA3auAA13LrhX4srsfAZwAfMHMpvW00HTZsbuZVVvqOWHSkKhLERHpVakGfLO7O+H4M2ZWeqAV3H2zu78c3q8HVhJc7i+rPL82uG7JiZN1cQ8RiRcLcns/CwQXJb2RIJzPAL5NMETBXe7+w5R2YjYBeAo40t3rOsy7CrgKoKKiYua8efO62YRAQ0MDZWVl3V7vtyuaeHpTKz8+rT8FfewbrD1tc1+Wa23OtfaC2txdc+fOXezuszqd6e4HvAEvE4T7d4HvAWeksl64bhmwGPjwgZadOXOm99SCBQt6tN4Z31/ol962qMf7jVJP29yX5Vqbc6297mpzdwEveReZmupgY38Hat29W6dGmlkhcB/wO3e/vzvrZsK2hibe2NrA+TOyrudIROSgpRrwc4HPmtl6wg9aAdz9PV2tEHbt/BJY6e7fP6gq0+T5tcEoxydOUv+7iMRPqgF/Vg+2fRLwSWCpmb0aTrvB3ef3YFtp8fe12ygtyufIMQOjLkVEpNelOlzw+u5u2N2fAbL6U8u/v7md4yYOoTA/NqMuiIi8LWeTrbq+kTdrdqt7RkRiK2cD/sW3dgIwWwEvIjGVuwG/bgf9CvN1eT4Ria2cDfjF63dy9LiB6n8XkdjKyXTb3dTKis11zBqv8WdEJL5yMuCXbKylLeHM1PDAIhJjORnwL63fiRkce4gCXkTiK2cD/tAR5QzsVxh1KSIiaZNzAd+WcF5Zv1PdMyISezkX8G9srae+qZXjFPAiEnM5F/AvrQsGGNMZNCISd7kX8Ot3MqK8mLGD+0VdiohIWuVewK/byawJgwlGMxYRia+cCvjqukY21e7V6ZEikhNyKuCXVO4C4Jhxg6ItREQkA3Ir4DfWkp9nTB+tC3yISPzlVsBX1nJYRTn9ivKjLkVEJO1yJuDdnSUbazla3TMikiNyJuDXbd9DXWMrx4xT94yI5IacCfglG2sBdAQvIjkjZwL+1Y219C/KZ+qI8qhLERHJiJwJ+CWVtRw5ZiD5efqCk4jkhpwI+ObWBMur6nT+u4jklJwI+Ne31NPcmuDosYOiLkVEJGNyIuBfrawF4GidQSMiOSQnAn7JxlqGlRUxZpBGkBSR3JETAb+0chdHjRmoESRFJKfEPuAbW9pYU9PAkWPUPSMiuSX2Ab9qSz1tCWf66AFRlyIiklGxD/jlVcEQwRpBUkRyTdoC3sxuN7NqM1uWrn2kYnlVHQP7FeoSfSKSc9J5BH8HcGYat5+S5Zt2MW3UAH3AKiI5pyBdG3b3p8xsQrq2v48/X8cxq56GtwbtMzmBc0PNDkYOKIFflWaklEw6prb2XW2Ou1xrc661F3KzzVNaB8OcOb2+3bQFfKrM7CrgKoCKigoWLlzY7W1MqaykX1sbtbW1+0xvbHPcgdZmamtbDr7YLNPWSZvjLtfanGvthdxsc3NxWY+y74DcPW03YAKwLNXlZ86c6T21YMGCd02796WNPv5rD/sbW+p6vN1s1lmb4y7X2pxr7XVXm7sLeMm7yNRYn0WzvKqOksI8Jg0vi7oUEZGMi3XAL6vaxRGjBmiIYBHJSek8TfJu4O/AYWZWaWZXpGtfnUkknJVVdfqCk4jkrHSeRfPxdG07FRt37qG+qZUj9QUnEclRse2iWbapDtA3WEUkd8U24Fds3kV+njG1Qh+wikhuim3Av76lnknDSikpzI+6FBGRSMQ24FdtqeewkeVRlyEiEplYBnxDUyuVO/dyuAJeRHJYLAP+9S31ABw+UqdIikjuinXAq4tGRHJZTAO+jrLiAo0BLyI5LZYBv3JLPYdWlGkMeBHJabELeHfn9S31HKb+dxHJcbEL+K11Teza26IzaEQk58Uu4FdtCYYo0AesIpLrYhfw75wiqYAXkdwWy4AfOaCEQf2Loi5FRCRSsQt4DVEgIhKIVcC3tiVYU9Og7hkREWIW8Ou276a5NaEjeBERYhbwq7c2AHBohQJeRCReAV8dBPyk4aURVyIiEr1YBfya6gbGDu5H/6K0XWpWRKTPiF3ATxmhS/SJiECMAj7hztptDUwZroAXEYEYBfz2vU5jS0JH8CIiodgE/KaGBIACXkQkFJuA37zbAQW8iEi72AR8VUOCYWXFGoNGRCQUm4DfvDvBlBE6/11EpF0sAt7d2dSgD1hFRJLFIuBr6pvY24pOkRQRSRKLgF8TDlEwVWPQiIi8LR4BXxMEvLpoRETekdaAN7Mzzex1M1tjZtelaz+rtzbQrwBGlBenaxciIn1O2gLezPKBHwFnAdOAj5vZtHTsa011A6NK8zCzdGxeRKRPSucR/PHAGndf6+7NwDzgvHTsaE1NA6PLYtHbJCLSa9I5ru4YYGPS40pgdseFzOwq4CqAiooKFi5c2K2dtCWcQ8vbmFTa0u11+7qGhga1OeZyrb2gNvemdAZ8Z/0l/q4J7j8Hfg4wa9YsnzNnTrd3dNqpsHDhQnqybl+mNsdfrrUX1ObelM5+jUpgXNLjsUBVGvcnIiJJ0hnwLwJTzWyimRUBFwMPpnF/IiKSJG1dNO7eamZXA48B+cDt7r48XfsTEZF9pfXipe4+H5ifzn2IiEjndG6hiEhMKeBFRGJKAS8iElMKeBGRmDL3d333KDJmVgOs7+Hqw4BtvVhOX6A2x1+utRfU5u4a7+7DO5uRVQF/MMzsJXefFXUdmaQ2x1+utRfU5t6kLhoRkZhSwIuIxFScAv7nURcQAbU5/nKtvaA295rY9MGLiMi+4nQELyIiSRTwIiIx1ecDPlMX9s4WZjbOzBaY2UozW25m10ZdU6aYWb6ZvWJmD0ddSyaY2SAzu9fMVoW/7xOjrindzOxL4et6mZndbWYlUdfU28zsdjOrNrNlSdOGmNnjZrY6/Dm4N/bVpwM+kxf2ziKtwJfd/QjgBOALOdDmdtcCK6MuIoNuBR5198OBo4l5281sDHANMMvdjyQYZvziaKtKizuAMztMuw54wt2nAk+Ejw9anw54Mnhh72zh7pvd/eXwfj3BH/2YaKtKPzMbC5wN3BZ1LZlgZgOA9wO/BHD3ZnevjbSozCgA+plZAdCfGF4Fzt2fAnZ0mHwe8Ovw/q+B83tjX3094Du7sHfsw66dmU0AZgDPR1xKJvwA+CqQiLiOTJkE1AC/CrulbjOz0qiLSid33wR8D9gAbAZ2uftfoq0qYyrcfTMEB3HAiN7YaF8P+JQu7B1HZlYG3Ad80d3roq4nnczsHKDa3RdHXUsGFQDHAj9x9xnAbnrpbXu2CvudzwMmAqOBUjO7NNqq+ra+HvA5eWFvMyskCPffufv9UdeTAScB55rZOoJuuFPN7M5oS0q7SqDS3dvfnd1LEPhxdjrwlrvXuHsLcD/w3ohrypStZjYKIPxZ3Rsb7esBn3MX9jYzI+iXXenu34+6nkxw9+vdfay7TyD4Hf/N3WN9ZOfuW4CNZnZYOOk0YEWEJWXCBuAEM+sfvs5PI+YfLCd5EPhUeP9TwJ96Y6NpvSZruuXohb1PAj4JLDWzV8NpN4TXv5V4+Sfgd+HBy1rg0xHXk1bu/ryZ3Qu8THC22CvEcNgCM7sbmAMMM7NK4CbgFuAPZnYFwT+6C3tlXxqqQEQknvp6F42IiHRBAS8iElMKeBGRmFLAi4jElAJeRCSmFPDS54WjLn4+6fHo8HS7dOyr0Mxy6Ru10ocp4CUOBgFvB7y7V7n7R9O0r5OB59K0bZFepYCXOLgFmGxmr5rZd81sQvtY22Z2mZn90cweMrO3zOxqM/vncACvRWY2JFxuspk9amaLzexpMzu8i32dCfw5eUI4Tv0d4RjmS83sS/vbpplVmNkDZrYkvOXK1/Elw/r0N1lFQtcBR7r7MfD2KJvJjiQYdbMEWAN8zd1nmNn/AP9IMFLlz4HPuftqM5sN/Bg4tZN9zQW+2WHaMcCYcAxzzGxQOL2rbf4v8KS7XxBe06CsZ80W2T8FvOSCBeHY+fVmtgt4KJy+FHhPODLne4F7giFQACjuuBEzGw3scPc9HWatBSaZ2Q+BR4C/HGCbpxL8Y8Hd24BdB99EkXdTwEsuaEq6n0h6nCD4G8gDatvfAezHWQTjHu3D3Xea2dHAB4EvAB8DvpjiNkXSRn3wEgf1QHlPVw7H03/LzC6EYMTOMLA7elf/e7j8MCDP3e8DbgSOPcA2nwD+Xzg9P7x6k0ivU8BLn+fu24Fnww85v9vDzVwCXGFmS4DldLj0Y9hXPtXdV3Wy7hhgYTi65x3A9QfY5rXAXDNbCiwGpvewZpH90miSIikws5OBS939c1HXIpIqBbyISEypi0ZEJKYU8CIiMaWAFxGJKQW8iEhMKeBFRGJKAS8iElP/H212iEQbhZwgAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "tfinal = 10\n", "tau = 1\n", "K = 5\n", "\n", "# define u(t)\n", "u = lambda t: 1\n", "\n", "# create a model object\n", "model = ConcreteModel()\n", "\n", "# define the independent variable\n", "model.t = ContinuousSet(bounds=(0, tfinal))\n", "\n", "# define the dependent variables\n", "model.y = Var(model.t)\n", "model.dydt = DerivativeVar(model.y)\n", "\n", "# fix the initial value of y\n", "model.y[0].fix(0)\n", "\n", "# define the differential equation as a constraint\n", "@model.Constraint(model.t)\n", "def ode(m, t):\n", " return tau*model.dydt[t] + model.y[t] == K*u(t)\n", "\n", "# simulation using scipy integrators\n", "tsim, profiles = Simulator(model, package='scipy').simulate(numpoints=1000)\n", "\n", "fig, ax = plt.subplots(1, 1)\n", "ax.plot(tsim, profiles, label='y')\n", "ax.plot(tsim, [u(t) for t in tsim], label='u')\n", "ax.set_xlabel('time / sec')\n", "ax.set_ylabel('response')\n", "ax.set_title('Response of a linear first-order ODE')\n", "ax.legend()\n", "ax.grid(True)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "eXvU5FXLd1Wt", "pycharm": {} }, "source": [ "## Encapsulating into a function\n", "\n", "In following cells we would like to explore the response of a first order system to changes in parameters and input functions. To facilitate this study, the next cell encapsulates the simulation into a function that can be called with different parameter values and input functions." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "colab_type": "code", "executionInfo": { "elapsed": 5880, "status": "ok", "timestamp": 1557657354747, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg", "userId": "09038942003589296665" }, "user_tz": 240 }, "id": "UfJ6CgtMd2ME", "outputId": "acf26c63-b532-47e9-facd-c5e625976389", "pycharm": {} }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABvEElEQVR4nO2dZ3gbx7Ww3wFAEuxdrBKpLsvqci9xjVsc24ntFJf0ODc39ctNrlKu026KU256dxLHJW5xt+y4i5bc1HtvlNh7AwmQBDDfj8GSEAWSKLvAgtr3efBIxLYZ7O6cMmfOEVJKLCwsLCxOPWyJboCFhYWFRWKwBICFhYXFKYolACwsLCxOUSwBYGFhYXGKYgkACwsLi1MUSwBYWFhYnKJYAsAiqRBCnC+EOCiEcAkhbojxXN8VQjwQ+P+MwDntujQ0RoQQPxBCtAshms3WtlAIIWqEEJ9KdDssIsMSACZGCFErhHAHXv5mIcQ/hBBZiW5Xgvk+8DspZZaU8im9TiqlPB44p0+vc0aLEGI68F/AQillaSxtE0J8TAjxhv6tNBYhRJoQ4sdCiOOBd+CgEOJrQggRtE+NEMIjhOgTQvQKITYLIb4uhEgL2ue7QojhwDukfboT0ikTYgkA8/NeKWUWsAxYDnwjsc1JOFXA7kQ3Qi+EEI4QX1cBHVLK1jCOF0KIuL3Hel9vnP4D/Au4DLgGyAZuB+4Afj1mv89LKbOBMpTQ/BDwfLCgAB4JCFDtk6dX+5MdSwAkCVLKZuBFlCAAQAhxjhDiLSFEtxBiuxDi4qBtHxNCHAloR0eFELcGff+mEOK3QogeIcQ+IcRlQceVCyGeEUJ0CiEOCSE+HbTtu0KIR4UQ9wXOu1sIcUbQ9lVCiIbAtv3aeYUQtoBmdlgI0RE4R8F4fRVCfDpw7c5AW8oD3x8GZgHPBjS5tBDHatfpE0LsEUK8L5zfVwhRLYSQ2oAU0C7/N/Bb9QkhXhJCFIX5239cCLE3cNwRIcRngrZdLISoD/xWzcA9Y9pxOfAyUB7o4z/GadsPhRBvAgPArFD3WwhxGvAn4NzJNF8hxHlCiI2BZ2KjEOK8oG2hrvfuwLPTI4T4HSDGnO8Tgd+gSwjxohCiKmibFEJ8TghxEDgYoi2XAVcAN0opd0kpvVLKd4DbgM8JIeaMPUZK2S+lrAGuA84F3jNeXy2CkFJaH5N+gFrg8sD/K4GdwK8Df1cAHSgNyQa8O/B3MZAJ9ALzA/uWAacH/v8xwAv8PyAF+CDQAxQEtr8O/AFwooRNG3BZYNt3AU/gmnbgx8A7gW3zgTqgPPB3NTA78P8vA+8E+pAG/Bl4aJw+Xwq0AysC+/4WWBvqNxnn+JuB8sBv8kGgHygbZ9/vAg8EtVcCjsDfNcBhYB6QHvj7rsl++8D29wCzUYPiRahBc0Vg28WB3/8ngf6lh2jXxUB90N+h2nYcOB1wALmT3O83JnnOCoAulJbtAD4c+LtwnOsVB653E+oZ+n+BPn0qsP8NwCHgtMD+/wO8FXQ9iRJyBeP0/y7g9XHaegz4TFC7PhVin7XAT8beY+tz8seyAMzPU0KIPtTg2gp8J/D9bcDzUsrnpZR+KeXLwCbUoATgBxYJIdKllE1SymC3SSvwKynlsJTyEWA/8B6hfM8XAKuklB4p5Tbgr6iBQeONwDV9wP3A0sD3PtSAtlAIkSKlrJVSHg5s+wzwLSllvZRyEPVS3jSO+X8r8Hcp5ZbAvt9AabDV4fxYUsp/SSkbA7/JIygN86xwjg3BPVLKA1JKN/Aoo9bXhL+9lPI5KeVhqXgdeAm4MOi8fuA7UsrBwLmj4R9Syt1SSi9q8J3ofk/Ge4CDUsr7pdK2HwL2Ae8d53pXA3uklI9JKYeBXwHNQft+BvixlHJvYP8fAcuCrYDA9s5x+l8ENI3T1qbA9oloRAkXjQ8ELDXts2aS408ZLAFgfm6Qysd5MbCA0Ye/Crg5+MFGDd5lUsp+lPb7H0CTEOI5IcSCoHM2yIB6FOAYSmsuBzqllH1jtlUE/R38og8ATiGEQ0p5CKXpfxdoFUI8rLluAm19Mqide1ECoyREf8sD1wRASulCadcVIfY9CSHER4QQ24KutYjJB4zxGNtXbQJ+3N8+0IarhRDvBFxY3SjBENyGNimlJ8o2adRp/wnjfo8gRiOKXEIIV+DrE37zAGPve13Q/8vHXF+O2V4F/Drot+lEWUPjnW8s7QR+yxCUBbZPREXgmhqPSinzgj6XTHL8KYMlAJKEgCb5D+Dnga/qgPvHPNiZUsq7Avu/KKV8N+qF2QfcHXS6CiFOmCSbgdKaGoECIUT2mG0NYbbxQSnlBagBQKLcHFpbrx7TVqeUMtR5GwPHAyCEyAQKw2lDQMO8G/g8yn2RB+xijH9aB8b97QPzEo+j7lNJoA3Pj2mDHil4TzjHBPd77H5aRFGWVMEFMOY3DzD2vgefpwmYrv0ReJamB22vQ7lpgn+fdCnlW+O1fwyvAGcHLNIRhBBnBa7z2ngHBo5ZCayb4PwWASwBkFz8Cni3EGIZ8ADwXiHElUIIuxDCGZhgrBRClAghrgsMnoOAC6Vxa0wDviiESBFC3Izy1T4vpawD3gJ+HDjfEuCTwD8na5gQYr4Q4tLAAOgB3EHX/BPwQ80FIIQoFkJcP86pHgQ+LoRYFjjXj4D1UsraMH6fTNTA0ha4zsdRFoDejPvbA6koV1gb4BVCXI2a0DSMSe53C1AphEid4BTPA/OEELcIIRxCiA8CC4HV4+z/HHC6EOL9ATfeF4HSoO1/Ar4hhDg90L7cwHMWFlLKV4BXgceFEKcHfuNzUM/hH6WUoSaOM4QQFwFPAxsCfbKYBEsAJBFSyjbgPuDOwGB9PfBN1GBTB3wNdU9tqJC4RpQpfBHwn0GnWg/MRZnSPwRuklJ2BLZ9GDXp2Ag8ifJVvxxG89JQk3ftKNfJtEDbQIXuPQO8FJjPeAc4e5w+vgrcidKim1CTqR8K4/pIKfcA/we8jRr4FgNvhnNsJEz02wfcZ19EzRl0Abeg+m4kE93v11Bhs81CiJCuk8C9vzZwjg7gv4FrpZTj7d+Ommy/K7D/XIJ+Zynlkyjr72EhRC/KCrs6wj7dCKwBXkAJtAeAvwFfGLPf7wLPVAtKQXocuEpK6Q/a54PixHUALiHEtAjbMyURJ7qCLaY6QoiPoSInLkh0WywsLBKLZQFYWFhYnKJYAsDCwsLiFMVyAVlYWFicolgWgIWFhcUpyniJmExJUVGRrK6ujurY/v5+MjMz9W1QgrD6Yj6mSj/A6otZiaUvmzdvbpdSFo/9PuECQKgc55tQq1OvnWjf6upqNm3aFNV1ampquPjii6M61mxYfTEfU6UfYPXFrMTSFyHE2JXegDlcQF9CpQawsLCwsIgjCRUAgZWT70ElHLOwsLCwiCMJjQISQjyGSimcDXw1lAtICHEHqhAEJSUlKx9++OGoruVyucjKmhrFtKy+mI+p0g+w+mJWYunLJZdcsllKecZJG2Qcc08Hf1BLz/8Q+P/FwOrJjlm5cqWMljVr1kR9rNmw+mI+pko/pLT6YlZi6QuwSZqsHsD5wHVCiFrgYeBSESjQbWFhYWFhPAkTAFLKb0gpK6WU1ahkX69JKW9LVHssLCwsTjXMEAVkYWFhYZEATCEApJQ1cpI1AGZESsm6g238Ze1hjnX0J7o5FhYWFhGR8IVgycwvXz7Ab147FPj/QR741FmsrCqY5CgLCwsLc2AKCyAZeetQO7957RA3r6xkzVcvZlpOGl96eBvuId/kB1tYWFiYAEsARIGUkp++uJ/K/HT+94ZFzCzK5K73L6G+y80jG48nunkWFhYWYWEJgCjYWNvFtrpuPvOuWThT7ACcO7uQFTPy+NubR/H5rRTbFhYW5scSAFHw8Ibj5Dgd3LRy+gnff/z8mdR1ull/pGOcIy0sLCzMgyUAIsQ95OPF3c1cvaiM9FT7CdsuP62E9BQ7z+1sSlDrLCwsLMLHEgARUrO/lf4hH9ctKz9pW3qqnUtPm8YLu5pPGTfQs9sbuebX6/j4PRuo6xxIdHPiis8v+cXLB7jqV2u5f88gnuFTKwCgs3+I//fINq7+9Truf7s20c2JO7sbe7jl7ne48Y9v8fbh5LT6LQEQITX728h2Ojh7ZuhwzysWltDRP8Tuxp44tyz+rNnXyhce2orX72fTsS5u+es7uAa9iW5W3Pj5S/v5zasHSUux89pxL1/91/ZENylueH1+Pn3fJp7b0YQA7nx6N/9cHzLl/JSkodvNrX9dz4EWF809Hj52zwb2NPYmulkRYwmACNAWfl0wpwiHPfRPd97sIgDeONQez6bFnUGvj/95ahfzS7J55vMXcM/HzqS+y81vXzuY6KbFhYMtffxl7RE+cEYlT3/ufN43N4XVO5pYe6At0U2LC49sqmPzsS5+etMSnv3CBVwwp4gfP7+Pzv6hRDctLvz4+b14hn386z/O5ZnPn0+208GdT+/SEl0mDZYAiIDDbS4aezxcOPekymojFGensaA0mzenuAB4dnsTDd1uvvme03Cm2DmjuoDrlpZz/9vH6HEPJ7p5hvO3N46SYhd8/erTALhmZgpluU7+svZIgltmPH6/5O61R1g+I4/rl5VjtwnuvHYhrkEvD22Y+mHQdZ0DPLeziU+cP5OZRZkUZqXxpcvnsflYF+uPdia6eRFhCYAIWHtADeoXzi2acL/zZhexqbaLIa8/Hs1KCI9sPM6sokzeFfRbfPKCmQwM+Xh+ik+C9wwM89S2Bm5YVkFBZioADpvg1rNn8Mahdg63uRLcQmN541A7tR0DfOy8aoQQAMwvzeb8OYU88M4xvL6p+9wDPLThOAK4/dyqke9uXllJjtPBoxvrEtewKLAEQAS8faSDqsIMphdkTLjfiqo8Br1+9jUnn08wHI6297OxtosPnDl9ZAAAWFyRy5xpWTyxpT6BrTOeF/c04xn2c8vZM074/oNnzsAm4OmtDQlqWXx4Yks9eRkpXLWo9ITvbzmriqYeDxtruxLUMuPx+yWPb6nnkvnTKMtNH/nemWLnmsVlvLi7OamyAVgCIEyklGw93sXKqvxJ9102PQ+AbXXdxjYqQbyypwWA9y49MRJKCMENy8rZWNtFa68nEU2LCy/tbqYiL53FFbknfF+cncaZ1QW8FPh9piLDPj+v7WvlsgUlpDlODIO+eH4xqQ4bL+5uTlDrjGdXYw8tvYO8Z0nZSduuX1ZB/5CP1/a1JqBl0WEJgDA53jlAu2uIFTMmFwAVeekUZ6ex7Xi38Q1LAK/ua2FBaTYVeeknbbt0QQkAr0/RydCBIS9rD7bz7oUlJ1g/Gu9eWMK+5j6Od0zNkNgNRzvp9Xi54vSSk7Zlpjl419xiXtrdnHSToeHyyp4WbAIumT/tpG1nVueT43Tw+gFLAEw5thxXZm04FoAQgmXT89g6BS2AXs8wG2u7uOy0k18AgNPKsinOTqNmigqAjYG5nUsXhO7/5aepgXHtwanZ/7UH2kixi3HnwS5dMI3GHg9H2qdmevQ1+9tYWZVPfmDuJxiH3cYFc4tYe6A9aQSgJQDCZPOxLrLSHMwryQ5r/2XT8zja3k/3wNQKi9t8rAufX3L+7NADgBCCi+YV88bBdvxTcDHc+iMd2G1iXEWgqjCDslwnb0/RdCAbajtZWplHRmroTPLnzi4ESNqFURPR5xlmd2MP584qHHefd80tprnXw8HW5AgEsARAmGw93s3S6bnYbSeb/aFYWpkHkJSLQyZiU20ndptg2Yy8cfc5a2YBPe5hjrQnx0sQCRuOdrK4IpfMtNADoBCCc2YVsv5IR9JogeEyMORlZ30PZ42zCBKgujCD0pypKQC3HO/GL+HMCfqvCcANSRIOagmAMBjy+jnQ0sfiirywj1lQpiyFvc19BrUqMWys7WJRec64GiAwMk+y+djUigZxD/nYXt/N2bMmLvpzzqwC2l1DUy4cdOvxbrx+OaEAUAKwgPVHOqecANx4VCk/E80DzijIoDAzla1JMv+XMAEghHAKITYIIbYLIXYLIb6XqLZMxoGWPoZ9ktPLc8I+pigrjaKsNPY1TR0LYNDrY3tdN2dUTzwAzirKJC8jZcoJgK3Huxj2Sc6ZOb4LAEbniZJlEAiX9Uc7sYnJ58GWz8in3TVIU8/UigTbUNvJ6eU541p/oATg8hn5bD2eHM9+Ii2AQeBSKeVSYBlwlRDinAS2Z1w0N04kAgDUhOi+KWQB7GroZdDr58zqiQcAW0BLmmoCYFt9N8CkkWAzi7LITLWzs2Fq5YPaeryLBaU5ZDtTJtxvSaUKj90+hYIghn1+ttd1c+Ykyg/A8hl5HGnvpysJ0mIkTABIhWYjpwQ+prQZdzf2kJlqp7owM6LjFpRmc6Clb8qsjNwRGACXhxEKu7Iqn8NtU2sSfHdDLzMKMsjNmHgAtNsEiypy2V4/dQSAlJJdDT0jg/tEnFaWQ4pdTKn+H2p1Mej1h9V/TUHYWmd+BSihReGFEHZgMzAH+L2Ucn2Ife4A7gAoKSmhpqYmqmu5XK6oj31rr5vyDFi79vXIDuweZtDr59F/11CepZ+sjaUvsfDqzkFyUmHP5rfZGyIGPhjZoVZDPvjvdSwstI+7X6L6Eg0bDg1QlWML2d6x/ciXQ7zSMMwrr63BEWbggFkIdU863H66BoZJ7W8J635VZApe33mUc9ITuyhMr+drXb3Kb9Vfv5+a7okTHg56JQJ4at12bM0nh4tGiyHvipQy4R8gD1gDLJpov5UrV8poWbNmTVTH+Xx+ufDOf8tvP7Uz4mN3NXTLqlWr5TPbGqK69nhE25dYufpXa+Xtf1sf1r5tfR5ZtWq1/MvrhyfcL1F9iZTu/iFZtWq1/P2agyG3j+3HM9saZNWq1XJnfXccWqcvoe7JC7uaZNWq1XLLsc6wzvHNJ3bIRd9+Qfp8fp1bFxl6PV/feXqXPO3Of0tvmP259Odr5Kfu3ajLtTVi6QuwSYYYU00RBSSl7AZqgKsS25KTOdY5QP+Qj9PLJzf9xjK7OAshmBLRIENePwdb+1hYFt48SFFWGtOy09g7RSbBtfoOi8J8DrQw4B1TxA2yu6EHu01wWpj3f2llHn2DXo5NkSJBuxt7OK0sJ+ww8NPKcpLi2U9kFFCxECIv8P904HJgX6LaMx7aBPDCCCeAQSWIqshL53Bb8q+KjCYSamF5DnuS4CUIh12aAKgITwBML0gnK83B/imSEHBXYy9zirNwpozvzgtmfqkKg94/BYIg/H7JnsZeFkXw7J9WlkN9l5tej7lToyfSAigD1gghdgAbgZellKsT2J6QHGp1IQTMmZYV1fGzirM4MgUsgGgioRaW5QQmz5InO+J47GvqozTHOZL+eTKEEMwryWJ/S/IPgKA04Eju/bySbISYGgIgGi/AaYF1QPuazN3/REYB7ZBSLpdSLpFSLpJSfj9RbZmII+0uKvLSw9Z8xjK7OJMjbf1JnxZhb3MvGRFGQp1WloPXLznYkvwC8GCri7klkSkB80uz2d/cl/QLono9w7T0DjI3zDQooOpjVxVksL8l+S0gTYhpVk04aK4ys7uBTDEHYGYOt7mYVRyd9g/KAnAP+2hO8vTIh1pdzC7OwhZBRIv2EiT7Wgi/X3Ko1cXcaeEPAKC04K6BYdpcgwa1LD4cCuS1mRuhFTyvJHtKWADaHF4kXoDSHCd5GSmWAEhmpJQcaetndnFk8f/BzC5Sxx5J8nmAQ62uiN1gVYUZpNhF0k+CN3S7cQ/7IrcAAhrzgebk7v+hgAUXaf8XlGZT2zGAZzi5XYAHW/ooz3VOuAJ4LEII5peodUBmxhIAE9Dc62FgyMfsGC0AIKkTo7kGvTT1eCIWACl2G1WFmRxOksyI43GwVb3EkWrAIxOhJh8EJuNgax9pDhuV+RNXwhvLvNJsfAHrKZk51OZiTgTuL43Z07I43NZvahegJQAm4HCr0tpnxWABlOSkkZlqT2oLQBvAoxGEs4szk94C0OYwIhWAhVlpFGWlJn0k0MFW5QYNNwRSQ7OANAGajGjuvzlRPftZ9LiH6TRxSghLAEyAprVHc/M1hBDMKs5K6kFQ0+CiiYSaMy2LYx0DDCdxOoyDrS6Ks9PIy4h8Vees4iyOJnlxFDX/Efm9n1GYgU3A0SRWfhq63XiG/RG7v2BUcTRzGLglACbgcKuLrDQHxdlpMZ1nViASKFk51ObCYRNUFUbmAgClBXn9kuNJvCDoYJQDIKjMqMl87weGvNR3uaPqf5rDTmV+RlJXB4tJ+dHcvyZW/iwBMAGHAxPAoWq/RkJVYSaNPe6kjYc/1OqiuiiTFHvkj4vmNkrWeQApJYda+sKuBDeWmUWZdPQP0TNg7gVB46G5QaNfB5OZ1BbQiACIwgtQnpdOmsNmauvfEgATcKTNFdMEsEZVQQZSQkOXW4dWxZ/DUfpAYdQMPmTil2Aimns99A/5oo4EmxmIAjvakZyDoOa/j1YAzCxSAsDME6ETcajVRVFWasgawJNhtwlmFmVaLqBkpH/QS2OPJ6YJYI0ZAddJMrpBhrx+jnUORD0AZDtTKMlJG9Ekk43adnXPZhbFJgCPJmkUWG17PzahrNhomFWUycCQj5be5FwLEes6oNkmzwRgCYBx0MxWPSyAGQXJKwDqugbw+eWIJhsNc6Yl7yT4sYDmHs38B8D0guSeCD3WOUB5XjqpjuiGCk1wJmsY9LHOAaqjvPegouCOdw6Y1v1rCYBx0AasWKS/xrTsNNIcNo53JJ8A0Noc7QAIo26AZKS2Y4AUu6A8Lz2q45N9IrS2YyC2ez9iASVf//sHvbT1DUZt/YBaC+CXcMyk774lAMbhcJtm+kb/8GsIIZhRkJGUFoCmAc+I4XeoKsikxz2clBOhte39TC/IiDgGPphkFoDHO/pjGgDLcpw4U2xJGQmlva+xjAFa7qxak95/SwCMw5E2F9MLMqJOAjeWpBUAnQNkpNopzoo+FFYTHsc6zfkSTERtR3/EpUDHkqwToT3uYboGhqkqiH4AtNkE1YXJKQA1rb2qIPr7b3b3ryUAxuFwWz+zYvB7j2VGoRIAyTYIHO8YYEZBRkyhsJoGZVYzeDyklBzrGIhZAMwqTs6J0FH3X+z9T04BELv1m5eRQrbTYQmAZMLvlxxt1ycEVGNGQQYDQz46TLwsPBTHOmPzAYP5taDxaOsbxD3so7ootv5rA+ixJAsFrY1xAlxjekEGDV1ufEmWEv1Y5wB5GSnkpqdEfQ4h1AJKsyo/lgAIQWOPWv6txwSwRjIOgv7ACt5YNcCMVLWaOvkGQH004On5agK5LsnWgejhAweYnp/BkM9PS5KlRD/eEfuzD+rdrzPpe28JgBBoCzdiSQM9Fu0lMuuDEIqWPg9DXv+I8IqFqgLzakHjoWnAsYQBAlTkpyNEct17UBOXxdlpZKSGnwY5FNrzk3T97+iPaf5DY0ZB5kg4tdmwBEAIjugYAqqhpdJNplDQYzqEgGok4yR4bXs/DpugIsoQUI00h53SHCd1XcnV/1hj4DWmawIgiSygIa+fxm63Ls9+VWEGwz5JU4/5+p/IovDThRBrhBB7hRC7hRBfSlRbxnK4zUWO00FRVuTLv8fDmWKnMDOVRhM+BONxXIcoCI0ZhRk093qSqjjI8cAiKEcUOZDGMj0/g/rO5Ln3oOYsZuhw78vznAiRXO7Phm43fhm7+w/M7f5NpAXgBf5LSnkacA7wOSHEwgS2Z4TDrf3MnpYVcxK4sVTkp9PQnTx+0GOdSgMuz3PGfK6qQpUPqT6JtOD6LjfTC2LT/jUqC9KTygJwB6KW9LAA0hx2ynKc1JtwABwPvSbAIUgAmND6T2RR+CYp5ZbA//uAvUBFotoTzJF2F7OizP0yEeW56TQk0SBwrGOAinx9NGBNkzSjFjQeDd1uKvNiHwBAWQDNvR7TpgQYy8gEsE6h0JUFGUklAEet39jvf3leOg6b4JgJn/3YZnd0QghRDSwH1ofYdgdwB0BJSQk1NTVRXcPlcoV1rNsraekdRLhaor7WeEjXIHWdXtasWROTdRFuX2Jld62b7BShy7V6B9UE2Kvv7MDWPBpWF6++RMqQT9LWN8hwTzM1NZ2T7j9ZP/pbh5ESnnzxdUozzT315nK5eH6tehXbju6lputAzOdMGRxkd4cv7vc62ufrrb2DpNph9+a32aODJ6DQCZv31VLjbI76HIa8K1LKhH6ALGAz8P7J9l25cqWMljVr1oS137bjXbJq1Wr5751NUV9rPP667oisWrVadrgGYzpPuH2JlaXfe1F+84kdupzL7/fL+f/zvPzfZ3ef8H28+hIph1r7ZNWq1fLxzXVh7T9ZP9Yf6ZBVq1bLmv2tOrTOWNasWSP/8eZRWbVqtWzpdetyzl+9fEBWrVot3UNeXc4XLtE+X3fct1Fe9n81urXjtr++I9/723UxnSOWdwXYJEOMqQlVRYQQKcDjwD+llE8ksi0aI2Ugp+kXAqqhRZM0dpt/MrB/0Ev3wDAV+fr4wIVQCdUakqDvMFq7IdJC6OOhzSUkSyhkQ7ebVIeNoszYquFpaP1Pmvvf7Y45+iuYyvx0U773iYwCEsDfgL1Syl8kqh1jOdzaj90mdIl+GEtlYDCtT4JwOO1F1fclyEiaAaB+RADo0/+SbCepdlvS+MEbutQAaIshCV4w05NsLUBDl1s35QfUe9TuGjJdFFwiLYDzgduBS4UQ2wKfaxLYHkBZADMKMqLOfz4R5UlkAWgDtV4DIKiXIFmqotV3DeCwCUpyYo+AApUUrSI/PWlCQet11oCTaTHYwJCXroFhXfuvvftmU4ASNgkspXwD0DfOUgcOt/brugI4mPyMFNJT7KZ7CEKhDdTR5sEPRWV+Oh39QwwMeWNeXWo0Dd1uyvKcMaWBHktlfvKEgjZ0ublswTTdzleclUaqw5YUi8EadLb+YNSSbuhy65pjLFbMHY4QZ3x+ydGOfl1XAAej/ODOpLEAHDbBtGx9NGBIrjmQ+i79QkA1pps4J0wwQz5Ju2tQVxeIzSaUAEyC/tcb4P7UfkuzPfuWAAiiocvNkNdvmAUAUJEkfvCGLv014IokmgOp7xrQVQMEtRaga2CY/kGvrufVm06PCtnVcwDUzme2ATAUmgWgpwAszVHvktnefUsABKGVgTTSRKtIIgtA7wFAG1DN9hKMZdDro7VPXw0YGFlRbcacMMG0uwMCQOf+V+Qlx0r4RgOsX4fdRmmO03RzYJYACELPOsDjYdZogLE0drup0NkFMi3bicMmTPcSjKWp24OU+oWAaoz4gU0+CHa4/YD+FkB5XjrtrkHTr4Y2Yv4H1O9ZbzLlxxIAQRxu6yc/I4WCTP2SwI3FrL7AYIYDudsrdMgBFIzdJijLc5reAtA7BFSjLEnmQNo9EpuA0lx9778WUNDcY24BqIXA6k15nmUBmJrDbS5DtX9Q+YDA3G6Q5h4Pfqm/CwACWpDJXoKxNHSriUq9BUBJdho2YX4B0OGWlOY4SdEhB1Qw5QGBYuZnHzT3p77WH6j3qbnXY6q6AJYACOJIm3EhoBrJYAFoA7QRL0FlfobptKCx1He5sdsEpTqtAdDQ/MCNSeACMkL4j66DMW//R6xfQ5SfDHx+aarKaJYACNAzMEy7a9DwGN2SHCc2gakHQU04GWUBaJXGzEp9l5vSHKcuWVDHUp4EkTDtbmmIC0RzKZm5/5r1W2lA/ytMGARhCYAAh9uNnwAGSLHbKMlx0mhiP6j2gJbp7AMG9RJIaW4/sN5pAIIpy0s3dVEgr89P16A0pP/OFDtFWWmmjoKqNyAEVCN4MZhZsARAgCMG1AEej9Jcp+kHwKKsNJwpdt3PrWlW9d3mXRDU1Ose8VfrTXmek6ZuD34T+YGDaekbVPM/Brj/QIVBmzkKSlN+9FwBr1FhwnQQlgAIcLjNhcMmRpJWGUlZrpNmE/kBx9LQbawGDOa1APx+SUvPIKW5xvS/Ii+dIZ+fjv4hQ84fK1rVLqPuv9ldYJp2boT1m55qpyAz1RIAZuRwq4uqwgzdIx9CUZJjcgug222IDxQYmVhtMmn/O/qHGPL5dSmDGQotCsysg6DmnjJiDgBGBYBKUW8+GroHKM42xvoF8yVEtARAgCPt/XFL0lSW68Q16KXPMxyX60WClJKGbrdhA2B6qp38jBTT+oE1wax3BJBGWZ65J0I1wWyEBqydd2DIR4/bfM8+qP4b4f7RqDBZTQxLAKAmvo4ZmARuLJp7wYxWQGf/EENeP2UGuUBA9b/JpH5gTTAZ1f+RhHgmvPegnskMB2SmGZOt1Yx+8GCaejyGzf+AsoCaTGQBWQIAVQB72CfjMgEMo9qVGd0gWpuMsgBALQgyY99htP96r4LVyE1PISPVbmoLoMBpXJZ2Tbs2qwLQ3OMx7N6Devf7h3z0mSQhoCUAgEOtWhnIOFkAAfeCGS2AEReIgRZAWZ7TtC6gph4PqXYbhQalAxFCUJZr3oSAzT0eCpzGDQsji8FMeP/7PMO4Br2Gub9gVLEwy7tvCQDgkJYFNE4CQKsyZcZIoKZeY33A6tzpdA0MmzIhXnOPm5LcNN1KIYbCzJEwTT0e8g20AAozU0l12EzpAoqL8mMy698SACgLoCQnjRxnSlyul+qwUZSVapqHIJjmHpUGoShLn2LgoTDbSxBMY4+HshzjBgAI5MU3Yd8HvT7aXYOGuoBsNkF5rjnTYRg9AQ7BFoA5BGBCBYAQ4u9CiFYhxK5EtuNwqytu7h8NtRjMHA9BME09Hkqy03RPhRuM9hI0mVQLLDNw/gOUBdDWZ760yK29gwCGWgCgLEAzWkBGR4CBSokuhHmUn0RbAP8ArkpkA6SUHG7rZ06c63SW5qSb5iEIxuhJMBiNhTdb/6WUcel/mcn8wBra/TDSAgA1B2S2vsNo/0sMFADK+k8zTf8TWplbSrlWCFGdyDY093pwDXrHtwD626FuPfiGoWwJFMzS5bpluU42HevU5Vx60tzr4bTSnNEv2g9By05wOGH62ZBREPM1RiyAHjeFiVZBgtAWgZVpA4DfD41boKsWMotV/1NiHxzKgsKAqwrjE3kWDtrEfL42CTzshuNvw0AnFM6GsmUgYhcOZblOWnpVOgwj51oipblXpUBJdQT672qD+g2Bd38pFMzU5TplJoqCS6gACAchxB3AHQAlJSXU1NREdR6XyxXy2F3tygx3NR6mpqZ25Hu7d4DZh/9BWdPLCEYzV3YUrODg3M/gSS+Nqh0a7o4hugeGeenVNaTaI3sJxutLrEgpaegcYG7GIOuff4h5B/5AfveOke1+Yaex/CqOzPoIfntsA2F2Cmzee4TqqmFD+hINtT3qWeioP8yOx59gzqG7yXA3jWwfdmRzdOYtNJZffdJAGMk9aXSp5+m1d7biPm6eV/DNoyo9Reqwi0P3f5nq2kdx+PpHtvdnTOfg3Dvozl8S03V6W4bx+iXPvLyGvDRjNYBI7svuIx6ybJJ1rzzP7MN/p6zp1THv/hkcnHsHnvSSmNrkGPJwqNsf8XNvxHtvnqdvHKSUfwH+AnDGGWfIiy++OKrz1NTUEOrYo28eBfZw4xXnj9YA7T4O/7wZ2g/CWZ+CRTcpze/gyxS+8SsKt6+CW/8F08+KrlNAZ049jx3czrylZ1FdFJkWOF5fYqXHPczgiy/x3vJuzt76X2BzwLu/D7MvhUEXth0PU7n5XiqHa+H2JyE7+hdh+vZ1iEwnWVn9hvQlGl7e0wJvb+KW7E2UbfgxFC+Aq76nLL+uY6S883vmHfwz85xdcMMfwT76+kRyT1yDXr75xosUVMzk4otmG9SbyKnp3U1h2hEubvg9xe3vwNwr4aw7ILcCGjaTufbnLNvxHbjmZ3Dmp6K+zvCeFu7fs4lZp69gSWWefh0IQST35a5tazk3v4cL930fOg7D2XfAohvBkQYHXqLwzV9RuGMV3PoYVJ4RdZvW9Ozi4NaGiJ97I9570wsAoznU6iLH6aBYi3rpb4f73wf9bWqQm3XR6M5lS9UD8cD74b7r4RMvqsEhCoJz4kQqAIyiucfDSrGfq7ffBYWz4LbHILdydIeqc+G06+CR2+D+G+Djz0N6flTXKjdhVsimHjeftD9P2YYH4PT3w/v+pF5+gGmnwbwrYd3P4bUfgM2uhEAULpGsNAfZaQ7TuAE0Wrr7+XXK7yhqXw9X3QXnfHZ047TTYOEN8Ngn4Ln/gpQMWHZLVNcJXgezpHKSnePIYHczXx34Ltj64aPPQPUFoxvLlsLiG9XYcN/18MmXoOT0qK5TlpdOn8eLa9BLlkErrsPFRB7YxHC4TUUACSGUz/exj0NPPdzy6ImDv0bBTPj4v9XA99CHlcCIgpFwsF7zREN0NB/jj6m/ZjizDD62+sTBX2POZfDhh5R19NR/QpRL2ktzzbcYLO1YDd90/BO54L1w419HB38NIeBdX4OLvwnbH4J3/hD1tUoCfnAzcWnzX7nA+w6H5nzyxMFfIy0LPvRPmHkRPPMFqNsQ1XVGn33z9H/AM8iP/b8gx9uuNPzgwV+jYJZ691Oz4KEPqbmRKDBTEECiw0AfAt4G5gsh6oUQn4x3Gw61upg7LVv98c4f4OhauPqnMOOc8Q/KLoUPPQj9rfDcV6IaBEvNFgsvJTPf+G+ycNN93T2QWTT+vrMuhiv+F/Y/DxvujupyZbnpdA8MM+gzR04U3N1cdej7HLNVIt73J6Xhj8dF/w3z3wMvfxuatkd1OTNNBAJw7G1uHHiUDXnX0FD53vH3s6fAB+6DnHJ44g4YdEV8qcLMVFLswlT9d7/+K86x7WXH0jth+pnj75hTrt793kZ4/mtRXctMmQDCEgBCiAwhxJ1CiLsDf88VQlwb68WllB+WUpZJKVOklJVSyr/Fes5IaO3z0O4aYkFZtop2efV7sOBaWPGRyQ8uXwYXfx32PA17nor42hmpDnLTU0zxEACw8zHK2t/kZ74PUjBz2eT7n/0fMOdyePX76mWIEE0L6vKYRAC89C2yvF38ueBrStOdCCHght8rK/DZL4M/8nh+U6UE9w4in/5P6mQx6+eHMail58ENf1LRUa/9IOLL2WyCadlOWszS/7b95K//Gc/5zmLo9A9Pvn/lSrhoFex6DPaujvhyZSNh0Im3gMO1AO4BBoFzA3/XA5HfeZOxr6kPgAWlOfDSt8CeBu/5Rfh+3fO+BCWL4aU7VchchJSZpTKYpwde+DrH0xfyvPPa0TC4iRBCTQb6h+GFb0R8Se0l6DSDADi+HrY+wMMpNzBQFOacTno+XPkjFSa66e8RX7Is10mbaxCvzwS1kd/5I6LzCHcOf5zCggksv2Cqz4eVH4UNf4HWfRFf0jQWkJTwwjfw2pzcOfyJkYJFk3LB/4NpC9W44R2M6JLTcpRr0QzvfrgCYLaU8qfAMICU0g2YJ4A3SvY29QKw2LMJDrwA7/pqZJEtdgdc+UPoqYvKH1xqlspgb/4aBtq5O/fzlORFMCFdMAsu+IqygI6/E9EltWyjnZ4ED4BSwst3IrNK+Kn7usjSACy+GWa+C9b8CLs3shKXJTlOfH5JuyvBlcH6WmDtz+iefjlr/Usj6/+ldyp/+IuRKwAlZnn2D74Eh19lQ9Wn6SQn/EWA9hT17nfVwvo/R3RJZ4qdwszUkbxbiSRcATAkhEgHJIAQYjbKIkhq9jX3UZaTRta6H0BeVeiJr8mYdRHMvwbW/SLiSaHSHBNoQX3N8M4fYdFNrHdXRr4M/rzPQ+Y05QqIYC5EW22ZcAtg33NQt56B81bR402NbAAUAi7/Lrg7qax/JqLLlpllIvSNX8Kwm62nfRWIMA12ZpGaDzn8GhxdF9FlywIusITmxff7lQuzYDYvZ76X/IyUyCqBzb4U5l4Ba38O7u6ILm2WuuDhCoDvAC8A04UQ/wReBf7bsFbFib1NvXwoby8071AP8tioj3C59H9gyBWxJlCa66TdNciQN4Fa8LpfgG8ILvmmyoMTaRqE1Ey48L+gdh0cqQn7ME0LSqgAkBJq7oLCuRydfgMQRSKwipWw4Fqm1z0dkQIwkhE2kX7gvhbYfA8s/TCHfcryjbj/Z34Sskrg9Z9EdFhprhP3sI9edwLz4u9/Hlp2wUWraOjzRZcF9NI7YbAn4mAIs6QED0sASClfBt4PfAx4CDhDSlljXLOMZ8jr51BrHzcPPAS5M2DJB6M/WcnpavJ4/R/B0xv2YWW5TqRUk9EJob8dttwHSz5Ef1YVvR5vdC/BGR+HnApY938RHVaa60ysADj4skpzceFXaO7zBtoURf8v+RYO30BEg4ApQgHf+o0S/hd+heYeD+kpdnLTI8yIm5Ku/OG166D2zbAPS3goqJSw9qeQPxMW3UhTNMoPqHVA866Cd34fUUSUWdy/4UYBnQ94pJTPAXnAN4UQVUY2zGgOtbo4h52Uu3bDBV9WPr1YuPC/1GTqxr+GfUjCS0Nu+At43XD+F0cextLcKKwgRxqc859qEGjYEvZhZbnOxEYBvfELyKmExTeP+GOjKgdYspCOgjPU7xlmMEBBZiqpdlvi/MD9HWryevEHoHA2Tb1qABTR5PpZ+THlBnzzV2EfMpoSPEFa8KFXVAjvhV8BuyO2JIAXfhXcXcqaChMtDNo9lNiMsOG6gP4IDAghlgJfA44B9xnWqjiwt6mXT9ufw5sxDZbfFvsJK1bA7MuUP90b3sReQvPiD/WrAWv+NVA8PygVbpS58Fd8BNJy4K3fhn1Iaa6TrkRNAh9/RyU6O+8LYE+hqduNwyYojLIOQt30G2CgHbY9GNb+QghKctMSFwq55R8wPKC0d2LMApuSrlJDHHwJ2g6EdUhJomPhN/wFskphyYfwDPvo6B8aTQIYKdPPVIvj3v69ShwXBqUmKQoVrgDwSjVbcz3wGynlr4Fs45plPA2HtnORfQe2Mz8Zve9/LOf+p1ocFua6AO0lSMiK0G0PKq3l/C8DOhTDcOYoTXDPUyoyIgzKctNxDZOYymAb7gZn7siaj+YeDyU5zqjrIHTnLYLy5fD278JeF1CWqJTgvmHY+De1oG/aAkCHNOBnfEKFUYcZDaflxU/IANhTryyA5beBI3WkDkJM/T/3c9DXBPvCWxeQcAsoQLgCoE8I8Q3gNuA5IYQdiE/5LIOYdfRBvDiwnfkJHU96KRTMVtpFGOQ4HWSk2uOvBUkJm++F0iUw42xgdDIyppfgnM+CsIXtC0/YisiBTtj7DCz5EKRmAETvA9YQAs79PHQeUVExYZCwUMh9q6G3QS3mA3x+SUtvjP3PKoYlH4DtD4c1GZ7qsFGYmaC8+Fv/CdIPK24HRgfhslhKQc65HPKrYX1477623iDRkUDhCoAPosI+PymlbAYqgJ8Z1iqD8bi6uNj9CvsKL4esafqd2GZT2RPrN4blCxdCqFDQeA8CTdvU5GfQiuemHk/kYXBjySlXLqVtD8Lw5H1KmAtsx6Nq8vOE/rtjLwRz2nWqbkCYC8O0hYBxD4Vc/2c1WM29AoAO1yBev4y9Fu45n1VzSlvC8w4nZDGY3wdb71fWT341QND8Vwz332aHMz8Nx9+C5p2T7h6cDDKRhBsF1Cyl/IWUcl3g7+NSyqSdA2hd9w+yhIfepQakHlp2i1ocE6YVUJKTgCXxW+5TBV4W3zzylXIB6FAL94xPgDugYU9CQhLiSQlb7oXyFVC6KPCVjN0CAHCkwvLb1aLC7rpJdy/JcTLo9dM9EJ7fWBeatqu5j7PuGMl3NOL+i7USVsnpMONc9fuGIdRKE5EQ78gatXBzxUdHvtL6H7MCsPxWlSU1jHDw9FQ7eRkpyeECEkK8XwhxUAjRI4ToFUL0CSHCj3c0E1KSvet+tvlnMXvZu/Q/vzMHln4Idj0elikcdy1oaAB2PgYLr1c5XQLoMgCCmgwrmBWWFpyQhHgNW6B1zwnaf/fAMINevz4CcOXHAkJmcv0oIYvBNOEflMq5SQ/3n8bKjyk3WO3kC8MSshBy872QUQgL3jPyVXOPh2ynI/bUzOn5yg228zEVETgJpTlOmnsSu542XBfQT4HrpJS5UsocKWW2lDJn0qPMSMMW8vsP82LqFcbV/lzxUeVi2PmvSXctyXXS2qfK48WFPU/DYO9JCe9aenWqhWuzwcqPKy2zZc+Eu2akOshMibMfdMu9SktbdOPIVzFPgAeTXwVz360G2kkiQuIeCTPsVs/kae89oY6Drv1feD0482DT5CGRpblOetxxDIV0tarFX0s/fELgR1OPW5++gxKAXrcSApNQlutMeDr4cAVAi5Ryr6EtiRNyy324SaV39nXGXaRsiSogseX+SU3hslwnwz5JR3+ccsJsuU9p6FXnj3wVcxjcWJbdCvbUsOKi89PimBZ40KUss9Pfpyy1AJobQjeF4IxPgKtZDTYTEHcLYO9qpZkuv/2Er5t7PKTabRRkpsZ+jZR0NcDufXbSWhlxD4Xc/hD4vScpP7q5P0HVTS5ZHJYFaIZ0EOEKgE1CiEeEEB8OuIPeL4R4v6EtM4KhAfw7H+N539msnG/wOrblt6uJ1qZtE+4W11DQ9kNqkmr57SdkPNUlDC6YzEKlCW5/RLmcJqDAaYvfS7D7SZWyY8wAoKsGDGpyNadSuRsmoDg7DSHi6ALbej/kzYDqC0/4uikQAhrVIrBQrPyoyhI7yZqIuIZCam65GedC8fwTNjX1ePRTfoRQ4aVN2yadDC7NSafdNcSgN3GLwcIVADnAAHAF8N7AJ+Z6AHFnz9PYh1084r2Y8+eEmfY2WhbfrHytW+6fcLe4RsJsvQ+E/aRSfrqEwY1lxUdVjpRJJoPznXG0ALbeD0XzYPrZJ3zd3OtBCDUg64LNrgaBw69B17Fxd0ux2yjOitNisK5jcPR1WHabctMF0dzjiTwJ4ERMOw2mnzPpZLCmcMRF+Tn2JnQcOkn4D/v8tLkGKdFL+IOaB7Cnhv3uawpYIgg3CujjIT46BtDHia3302gvp6/kLOP8/xrpeSoscOdjE2rBcTODfQGNbN5VqqJZEDGlgRiP6guUq2kSU7jAKeKTEK91H9StVwPAGE23ucdNcVYaKXYdC+Qtv1X9u+2fE+5WmhunMOBtDwIiZB3fZr3mf4JZ+VE14B5/e9xd4hoEsOU+SMtVdY2DaO0bREodrT+AjAKVG2zHIxOGQ5uhKmC4UUCVQognhRCtQogWIcTjQggTlXOenPSBRjj2Jg94LuTapeXxueiK2yfVgguz0nDYhPFZIQ+8oArdr7j9pE2jYXA6WgBCqMH22JuqfvA45DvVYGy4Frj1frClqMVfY2juHdR/AMyboeonb31gwpXBKhLE4Hvv9ylBNPsSyJt+wiYpJc2xLgILxcLrITV7QgUgI9VBjtNhvAXk7lLBD0tuHln4p6HLAshQrLgdPN0Trgw2w2rgSCqCPQOUoxaBPRv4LiaEEFcJIfYLIQ4JIb4e6/kmorT5VfzYeMz3Lq5dUmbkpUapukAtNpnAFLTbBNOy04wPB9tyv8p9MufdJ21q7vGQnaZDGNxYlt6iXE4TDAIFAQFgqAXkHVQTgAuuUStWx9Dc49bXBaKx4iNqxe2hV8bdJS4TgUdfV7HvIXJedQ0MM+T1628Rp2bC4ptg91MThkSWxiMMesej4PWELPWqvXe63/+ZF6ssw1vHf/dH1sGY3QIAiqWU90gpvYHPP4CT36QICKST+D1wNbAQ+LAQYmEs5xwXn5eS5td4Uyxn/ty5VBVGUPUqFmw29dIdewM6Do+7W6nR4WA9DXDoZeWWsJ88yDf1uPX1gWpkl8D8q9XgO06CvHynLdAGA1+C/c/DQMe4tZ5jzoMzHvOuViuDJxCApblOej1eBoYMzIu/9QEV9rng5Gm70fkfgwTgJCGRpbnpxgp/Le1J2TIVmTcGw/pvs6n37UjNuPNA2c4UstIcCU0IF64AaBdC3CaEsAc+twEdMV77LOCQlPKIlHIIeBiVbE53jm18BudQJw8MvYsvXz7XiEuMz9JbVH6cCSIiDNcCtz2ocp+Mk/U0qkIw4bLiI8r1dOCFkJtHLAAjzeAt90HudJh1yUmbBoa8gToIBvTfkapCIvf/W1VeC4HhdQEGOlX45+IPhEx6OBICa0T/y5dDyaKJBWCOsfmAsvsOQutuNScRgpZeD2kOW+R1EMJh2S2ASOy7Pwnh2vyfAH4H/DLw95uB72KhAgheL18PnD12JyHEHcAdACUlJdTU1ER8odR37iFT5lAw+2z6ju6g5mh0DY6WxfnLyVp/D2/bzlMukTF4ewdp6PSG3TeXyxX+7yD9nL3+bjx5i9m+4zhw/KRdjrUNsLjIHtVvO/n1HZybWojrlV+ys/XktYM+Tz9Ou2Dj7kPM80+ePiFSnO4Wzj68htrqD3Fs7cmrU5v71eRzV8NRamrqo77OePck3buAs6WPI0/8gONVN518/Q41P/DC2vUsLIwhD9M4VNQ/x1zfIJt8C3CFaN/rx9Vitdo9W+k9ovTBiJ6vya6ffS5zD93Npmf/jit71knbB7uGaOsb5pXX1uCIMhPrRMw6/hw+WxpvdZfiC9Gn7Qc95KVKXn/9dd2vDbAkfykZ7/yNdzhbKYJjSPO52V83ENbvred9GUFKmZAPcDPw16C/bwd+O9ExK1eulNHQ2N4lX3v871Edqwu7n5LyOzlSHngp5OY/v35IVq1aLXvdQ2Gdbs2aNeFf+/Aade3tj4bcPOz1yZlfXy3/78V94Z8zUl79gZTfyZWy6/hJm9asWSMv+78a+R/3bzLm2q/9cNxrSynlmwfbZNWq1fLNQ20xXWbCe/L3q6X81VIpfb6TNh1pc8mqVavl45vrYrr+uPzxAin/dOG4m//vxX1y5tdXy2HvaNsier4mo79Dyu8XS/ncV0NufnD9MVm1arWs7xrQ75oanl7p/d40KZ/8z3F3uemPb8oP/vkt/a+tsfMx9f4dejXk5q8+uk2e/cNXwjpVLPcF2CRDjKnhRgHNEkI8K4RoC0QCPS2EOFmcR0Y9EBySUAk0xnjOkJQV5iEKZhpx6vCYd7XKPzLOhJChlcG23Kfy3p8WetlGm2sQv9Q5AmgsmutpnJBIw/Ih+X3K/z3nspOiXzRGQmCNDAte8VHoOqrmgsZgaFbIpu2q3vXykyO/Rnbp8VCcnYZDzxDYYDIKVOqJHY+ErJZm6EToriew+z3jun8gsAjOyHs//z0qNcbWB0JuLg2kgvH6ElMYKdy7/iDwKFCGigT6F6o2cCxsBOYKIWYKIVKBD6EijaYejlQVfrjv+ZDL4w1bCzDQqZbkL/mgWqIfgtEQUB3XAIwlv0ql3x0nJFKFQhowABx+TUXhjDP5CzpmgpyIhdcpIRxiZXB6qqrDa4zwv18VaVl8sutJQ60BMFD4g/r9PT1qLmIMhs6BbLmP/ozpUHlmyM3+QB0EQ/uf4lQLw/auDpkcsjTXiV8qRSwRhCsAhJTyfjkaBfQAEFP2MimlF/g88CKwF3hUSrk7lnOamuW3qeXxOx49aZNhq4FD5L0fS8ylIMNl5UdVKOLhNSdtKjNKC9pyL2QUKQtsHFp6PYHCPDqHwAaTkq6E8N5nQg8COQYUhhl2w85HlfAJSvw2FrUK2EDhDyr1RH61uh9jKAs8d7rHwrfshoZNNJVdcdLCP43OgSGGfdL4/i+/HXyDIaOhEloWlvAFwBohxNeFENVCiCohxH+jKoMVCCEKor24lPJ5KeU8KeVsKeUPoz1PUlCyUOWg3/rAScvjpwUeQF0XxGi5T8qWQenicXfTBIBhUUAa869RbrAQg0Bpbrr+WpCrVUXfLPuwssDGQaXBNlj4wWiG2B2PnLTJkEiQcRK/jaU5Hv3XwqFr16lU0UHkpDtITzGgKt7me8GeSnPpxePuMqL8GN3/siWq+l4IF7CmeCUqEiiSimCfAdYANcBnUVFAm4FNhrRsKrLidhWS1rj1hK/THHYKM1P1TQnQsEVdawLtH5QLINVhIy/D4AqfjrRASOTz4Go7YZMhWpCW+XH5xP1v6fUYEwI5ltJFSgHYfHJ+HEPmQLbcC3lVJyV+C8Y16KVv0Gt8WhRQGWKF7SRfuBBC9V/PZ3/YDTsehtPeizdl/Kz1zfFw/2ms+Iiaj2nafsLXhocBT0K4uYBmTvCJdTL41GHRjeBID6kJ6F4ZbPM9kJJ5QtWvUGiFYHTLBDkRy29Xg/L2E6ePdJ8I9PvVQDvjXCieN+GuumaCnIyVH4W2vVB/os5UmuvUNx+SVpBl+e0nJX4LJm7WH6hyoXPererx+k5c9Ka7BbTnGWX9rBh/8hcYETpx6f/im9R8zBgBmJeRQprDlrDFYOFGAd0shMgO/P9/hBBPCCGWG9u0KYgzV+VICZEgTlct0NOr8t4vvvGEvPehaOnxxEcDBJi2QGXi3HLfCVqw7hZA7VroPKwK00zAsM9Pu96ZICdi0Y1KKG/5xwlfa/3XLR/S1n8qbTtE4rdgdK+DMBkrPqLqJIxJjaG7ANj0dyiYPaH1A+rZt9sERVkGzwGAmoc57b1qXi4oQdyIBWRmCwC4U0rZJ4S4ALgSuBf4k3HNmsIsv01V5BqTJKpEz/qoOx+F4QFVnWgSmnp1rIYUDis+Ch0HT8gSmZuegjPFpt9q4E33qBdu4cQLyw3JBDkRadmw6P2w6wklpAOMhAHrcf/9PrXydPZlkFsx4a6610GYjHlXQua0k1YGlwWefZ8eVfFadkPdO6oozwTWD6j+T8tOw27AArSQLL9NJYjb/9wJXysBmJiEcOEKAC127z3AH6WUTwM6lA86Bak6PxARMeYlyHHS0a9DcQgpYdM/1MRv+YpJdpW09BiQCXMiTr/hpCyRSgtK10cL6mtRwnXZrSoEbwJGI6DiLACHB5SFFqBcTwvowIvQ1xgy6+tYWkbSgMep//YUNSl/4IUTUmOU5qbj9Us69AgC2HSPcrVMYv1AYP4nnvd+5kUqQdxJAlCnZz8KwhUADUKIPwMfAJ4XQqRFcKxFMDabGgRq10HraJVNzQ0Rc3GIxi2qEtnKj40b/qbR2T/EkM8f3wEwOEuku3vka93WAmy9X80zhGH9xHUSUKPyDJi28IRBYHQORActcMOfIbtcRV1NQlOPm7yMFJwp+qegGJflHwHpO2EeqEyvxXCDLtj+sCr5mTF5cKKutYDDwWZTbrAjNSekSC8NWEBxqwse3KQw9/sAKl7/KillN1AAfM2oRk15VnxUVQtb/+eRr3Tzg2/4a1iTv8HXiutLAKNZIneNxkXr4gf1+1T0S/WFUDR50r+4rAIeixDq/jduGSkZqGWFjLn/bfvV4HLmJ5S2PQnNPYPx7TtA0RxlBQfVy9atMMqux2GoT7l/wqCldzC+FgAoxcSeesK7X5oT57rgQYQbBTQAtAIXBL7yAuNX+bCYmMxCNUBvf1gVq0Cn1cB9LWpQXX6rmnCehFEXQBzi4IMpX65cVEErY3XRgvb/G7qPhz0ANPe4SYtHCOxYlnxAuSnG9D9mC2jD3WpwWfGxsHZv7nXH1/rRWH67mqQ/9iYQHAoZgwUkJWy8G6adDtPPmnT3Ps8wrkFv/JWfrGJYdJOapwnUSUhkXYBwo4C+A6wCvhH4KgUIndzCIjzO/ozSggPFYnRxA2z8qyr9ePZ/hLV7UyJ84DCqBTfvIKvvEABlecoP3N4fgwvs7d8pH+tp14W1u1YJLC4hsMFkFKgVujseHcmPE7MF5OlVbpVFN4YsehOKhFgAoCbn0/PhnT8CUJCZSqrdFttagCM1yqI657OTuj4hQe4/jbPvgOH+kZDQRFYGC9cF9D7gOqAfQErZCGQb1ahTgtLFqmLYhrvBN0y2M4XMVHv0lcGGPbDpb6rmb+HssA5pDoTB6VYMPRIW3wypWUyvexoY9QNHrQXVb1KRRed8NmTRm1AYVgksHFZ8RJUL3fkvQAnhmAaALffCkAvOuiOs3Ye8KgQ2IQNgagac+WnY9xy0HUAIEbsF9NZvIKtEWVdhkBD3n0b5cph+jnID+X0j98DwsqghCFcADAVSikoAIUScSmpNcc77PPQcHx0EYqkMtu2fqurVuf8Z9iHNvR6Ks+IYBhdMeh6c8QmmtapqaTH7gd/6rSr6HUb0i4YhxdDDpfpCVaHqjV+CzxvIhzTIcDT5kIbdqv8zL4KKiSO/NFr7EjgAgrKAHWlq4CbG0pDNO1XiP+2cYTA6/xVn96fGeZ+H7mOw63GKMlVd8EREAk0qAISyj1cHooDyhBCfBl4B7ja6cVOeeVcpS2Dtz8HnjV4L8g6pgaTyrEkXvwRjWCnEcDn3c0hhhzd/HduS+Lb9quj3GR9XsfZhkJAQ2GCEgHd9Ta3a3f0kpbnpSAltfVFYgFsfAFeLOl+YJNQFApBZpOLidzwCvY2UxWIBvPFLFfgQ5twPjObdmmZ0IrjxmP8eFQ229mfYUDWZTTkHEND8bwAeAx4H5gPfllL+1timnQIIARetUhNiu5+gNCc9uodg2z9Vps2LV4Xl/9SIexjcWLJLaSq7HLY9SIG3RfmBo+l/zY9VeOl5Xwz7kISEwI5FGwTW/ZyyXLWsJuL+e4fgjV8pl0L1BZPurtEc7zUAoTjvC6pU6dqfjSg/UkYYBNC8S0X/nP2ZCbOejqWp10NBZmp8Q2CDsdngov+G9gOw+8mErQYO1wX0NtAtpfyalPKrUsqXjWzUKcX896jIhZofU55to7VvMLIVkUMDyoKoOEOt/oyAhITBjeH4jBtB2BCv/TC6FZHNO2H3k2riO7Mw7MMSFgIbjM2mtPa2fcxvUatDI1YANv4Veuvhoq9FJPxH8gAZnQZ8IvKrVbqOzfcyz97CkM9PZ6ShkGt+qFx/54cv/CHOKVDG47TroXgBvP4TynMSUxw+XAFwCfC2EOKwEGKH9jGyYacMNhtc/l3oPMJF3U9GviLyrd+oAeCK/41oAEhYGNwYBp3FauJ2x8Ocm1EXmRYkJbzwDVVx6bzPR3TduOfBGY+FN0DFGZRu+AkZeCKbCB7ohNfvgtmXRiz8m3s8OFNs5KQbWAchHN71NXCkcf6xPwARWkBH16nssud9ISLtX7tOop99bDa47NvQfoBrPc/R1OOO3AKKtQlh7nc1MBu4FHhv0MdCD+ZdAXOvYOmRP1NET/iaQHedMv9Pfz9UnRfRJRPuAw7mwq9ARiGf6f8zrT394R+38zG1ovqyb0c1AEACJwE1bDa46i5s/S18KfWZyCyANT+EwT648kcRCX9QLpCy3PT4h8COJbsELvh/lDa8yMW2reH33zsEz31FpbyOUPhDAtJAjMf8a2DWJVzU9FcyhrvpcQ/H9fLhLgQ7FupjdONOKa78EXbfED9K+StN3WFogX4/PP05lfXx3d+P+HIJDYMbizMXrvwRs9y7uLLvyfC0oP52eOlbKqQujLQPY2np9WATUJRlgpRW08+EpbfwKdszpLVsCe+Yo2uV++esz8C00yK+pHKBJGgCdCznfwlvwTx+kHIPbZ0nV0wLyRu/VP7za34+brnT8Rj0+ujoH0q8BQBKcF/9E1J8g/wg5e/hvfs6YuXzMQtFcxm48FtcYd9M7p7QxeNP4O3fwdHX4aofjVvwfCJMowFrLPkgx4sv5v/ZHqH36CSDoN8PT31W5RK67rdgi3wir9noYuiRctWP6bQXcVtjQKufCFcbPPlZlfL4sm9Hdbm4VUILB0ca4rrfUEYHS7d++6SCOSdx7C3l+lp8s7KeI0TLt2UK5QegeD6NK/+La+wbENvCePd1JCFPf6C+wG4hhF8IcUYi2mBGMt71BV73L+XMvXfBoVfH3a+wfT28/G2VX3ySohfjkfAwuLEIwaGzf0gH2WQ8fgv0NITeT0ql+R98Ca784YTlLiciLsXQIyE9j0crv0WxrxkeuV25OEIxNACP3AYD7XDT39SiqgjRiqGbwgUSwF59Lnen3MrCjpdg7c/G37H9IDz6EcifCdf+MqprNZnJ/RnAdt4XWOdbxLyN31GrmuN13bhd6UR2Ae8H1ibo+qbEZrfzg/T/piW1Ch768MlFpKWErQ9w+u6fKNfH+/4Sse9XI+FhcCEoLKnkk0NfQwz2wd+vVCF+wXiH4N+r4J0/wNmfhTM/FfW14lIMPUIGys/hm95Pw5E18M+bTi4g39cC910P9Rvghj+qZyAKOvqH8PqlOVwgQbxccAvr0i9Tcxsvf1ulNQmmfjPcG5h6vOWRsNd8jEWbaDdT/4tzM/i890t0ps+ABz+oakbEgYSEAEgp9wKJn4AyIbl5BXyXu/hL2i/h8U+qGP8F16p46T1PQ+06evKWkP+Rp6LS/jTUAGieFwDUC7lXVvHSmX/jmp1fhr9cDEs/BDPOUQudtv5TrZk49/Pw7siinsbS3OPh/DlFurVdD8py0/m99yK+9Z6F5L7yNfjdGWp+o3AOtO5RKaSHPXDzvSqXUJRoE61msgAASvPS+V7/53hl5XR489dw4CWV1z89X8157HoMcivhtifCyvY6HiMRYCYSACl2G87sAn4/41d8p/+H8NjH1bs//xr17h9+DWfeDbpfN8ExYJMjhLgDuAOgpKSEmpqaqM7jcrmiPjae2AY9bO/18/oF/0Wl41kq654l7fBrAHjSiqmb8ykO5L6LzHe2TnKmiTnU4CbfKRL+mwTfF7+U2AX8+5id3CU/YebRBynZ/i/sgRrKvdlzqV38bTrTVsLa6I1Ht1fSN+ilv72Bmpq2yQ8IAz2er/ZWVSv38ZbpLFn2E2YevZ+Cdb9E4Edio73oTI7M+gju1hxojf5aWwPXaTy8m5r2fSdtT9S74u0d5Hi3lzVZN1B8ehnVtQ+T9fKdaps9k+by91Bb/UG8e1thb2tY5wzVl417B3HaYfPbb5hKCc0Uw2ys8/L6yq8x3fE0FcefIy1QPtOTVoisWqb/fZFSGvJBpYvYFeJzfdA+NcAZ4Z5z5cqVMlrWrFkT9bHx5H+f3S0X/M+/pd/vV1/4fFJ210nZXa/+L/Xpy/LvvyS/8cSOmM8TK2P7ct6PX5VfeWTb6BfDg1J2HJGyv0O3ax5o7pVVq1bLp7bW63ZOPe7JroZuWbVqtfz3zsbRL909UrYfknKwP+bza9z31lFZtWq1bOlxh9yeqHfl7rWHZdWq1bK7f2j0S1ebuv/eofEPnIBQffnsA5vkJT8/+ftE85n7NsnL/69m9AufT8qu41L2NEjp98d0X4BNMsSYapgFIKW83KhzT2VKc524h330ur3kZqSoOPHcSl2v4Rn20dk/ZDoXEIRIiOdIhYKZul6jMeACKc8z0SQwoxFZJyyGcuaoj440dHtIscepGHoEjPS/162efVA5gzL1ddU1dnsoN1MAQIDSXCdvHmof/cJmiyrCLxJMEgNnoTFSF8DAZeHNJh0AIcaskGHSGIi1Nlv/8zNSSHXYDE8K1titCsHYEpEFdgJ0qww2CY3dbsrzzKf8lOU66Rv00ueJ32KwRIWBvk8IUQ+cCzwnhHgxEe0wI6Uj9VGNWxAyOgCa8CXIiTIpWAQ0dbuxCShJRB2ECRBCxCUpWFOP25QacEwZYcNkyOunzTVonjUQQSSiLkBCBICU8kkpZaWUMk1KWSKlvDIR7TAj8XgIGgICoMJkGjCo/g8M+ej1eA27RkO3ioE3zSKwIErjkBa4sdtjOusHoDg7DZsw1gJo6fUgpUmf/ZGiSDFUxYsQ870BpzjTso03gxu7zbcQRkPTzIwcBBOeBnsCynKdNBpo/fn8kuZejymtvxS7jeLstNjKok5Cg0ndfxA8BxS/dBCWADAZqQ4bRVlphloAjd1uirLSSHOYZxGYhiaUjBwEG7vdlJlwAABVG7ml14M/kpTgEdDa58Hnl6Z0gQCU5qYbrPyY1/2prcqPZ2EYSwCYkNLcNGNfgh43FSZ8AcB4P7CUksYejyldAKD6P+yTdESaFz9MGk3s/oPROSCjMF0OrCCcKXYKM1NpmupzABYTE3VlsDBRURDmewHAeD9wR/8QQ16/aV1Ao35gY/qvuf/Mev9jLg4/CQ3dbgoyU0lPNZ/1C8b3fyyWADAhpblphoWBSilNOwkIxvuBzRoCqmG0H1jrf5mJLUAjQyHNGgKqEe/SkJYAMCFluel0DwzjGfbpfu7ugWHcwz7TDoBgrB94RAM2oQsAjI+Fb+rxkJ3mIMeZYsj5Y0WbmzGs/90mSoMdgqjKosaAJQBMiOYGaDSgOMRIFIRJXSBgrB/YzJOAAIWZqaQ6bIbce1D336zaP4zOTTR0GWcBmXX+A5Ty12WQ8hcKSwCYkMr8wEtgwCBgdhcIGOsHbepxk+awUZBpgkpgIbDZBBV56dQbJACaesw7/wOjz74R/e/1DNM36DWt8Adjlb9QWALAhFTkG6cFNZk4DYSGkX5gbf7DTFkgx1KZn069YRqwuV0gxVlppNptxjz7Jp8AB2OVv1BYAsCElOY4sduEYRZAqsNGoUk1YDB2NXSjiReBaVTkpRsyALqHVBJAs4YAg7KAyvKc1HcN6H7ukQlwEwvAygJV48MoF9hYLAFgQhx2G6U5TkO0wIZuN+UmTAQWTMismDph5hBYjYq8dNpdg7r7gbXIIrP3vzI/3RDlx8wpUDRKstOw24RhFuBYLAFgUiryjdECG7vdptaAYHQxmGay68Wwz09r36CpJ8ABKguMcQNoEVBmv/9GWUCN3W4cNkGxyZIABuOw2yjLNcYCCoUlAExKZZ4xWpCZ1wBolOY6sQn9JwJVllHza8AVeca4AbT0GmbWgEH1v7VvkEGv3haQSgJoN7H1CwEBaM0BnNpU5KfT1ONm2OfX7ZzDPj8tfR5T+4BBLQYry02nvlNfLcjMicCC0YIA9HYDaD7wklzzasAw2v9GnS3ABpOHgGpU5mdYLqBTncr8dPxS35QAWipcsw+AoAYBvV+CuoBAmR6YaDMrJdlpOGyChm59BWBdp5uSHHMmAQzGqLUADV3mXgWsUZmfTnOvhyGvfsrfeFgCwKSMuAF0NAXrOgMugHzzCwAVCqnvAFjf5UYI8y4C03DYbZTm6h8EUNc1wPR8cws/CA6F1O/+D/v8NPW4TS/8Qb2fUsYnLbQlAEyKEW6AusCAmgyDwPT8DJp01oLqugYozXGaXgMGYyZC6zsHkmIAHJkD0rH/Td0e/DI5nv1KA9cBjcUSACZF01L1fAjqOwewieRwAVUaoAXVd7pHXi6zU5mfoav1N+T109zrYXoS9D8lEAat57OvKT9ahJWZ0YRUPOYBElUT+GdCiH1CiB1CiCeFEHmJaIeZSXPYmZadpqsZfLxzgLLcdFId5pf7lYGXQHNb6UF9krhAQFmAevqBm3rc+OXoQiOzU5GvbzqMkfmfJLj/oxaQ8aGgiRoJXgYWSSmXAAeAbySoHaZG74nQui4305NAAwJG2qnXSzDk9dPU60maAbAyT1lAegUBaII0GQZACFhAOlsAdpsw/SpwGLWAjMoHFUyiisK/JKXUqn6/A1Qmoh1mR283QF1n8mjAWjoMvQRgY7cbKUkKFwgEJ0XTRwCOzP8kiQJQkacsIK9OYdB1nSoCyGE3v/UL8QsFdRh+hcn5BPDIeBuFEHcAdwCUlJRQU1MT1UVcLlfUxyYKX+8Q9Z3DvLZmDbag5GXR9GXIJ2ntG8Tb02Kq32GivuSnwaZ9R6lJa4r5Orvb1aKijuMHqHEdjvl8Y9H7+WodUAPfq29vZagu9tz9bx4Ywi7gwLb1HJokEZ4Z3hVX6zA+v+Spl2ooSo9+0Nb6srvWTZadhPcrXOxDHg50+k9oryH3RUppyAd4BdgV4nN90D7fAp4ERDjnXLlypYyWNWvWRH1sorjv7VpZtWq1bOp2n/B9NH052NIrq1atlk9uqdepdfowUV8+9Oe35fv/8KYu13lw/TFZtWq1rOvs1+V8Y9H7+Roc9snqr6+Wv3hpvy7n+/yDW+SFP3ktrH3N8K6sPdAqq1atlm8fbo/pPFpfVv7vy/K//7Vdh5bFh5+/uE/O+sZzctjrG/kulvsCbJIhxlTDLAAp5eUTbRdCfBS4Frgs0ECLMYy4AboGRjJkRsuIDzhJXACg+v/6gTZdzlXXOYDDJkyfB0cj1aH8wHU6rYau6xxIqns/IzBXc7xzgHNmFcZ0LveQj3bXYFL1vyIvHZ9f0tTjMTR0N1FRQFcBq4DrpJTxyXqUhFQFbvyxjth/omRaA6BRma9ywuiRFbOuS2UBNXsemGCqCjM4ppMAqO9yU5mXPPe+Ii8dh01wrKM/5nPVdyXHCvBgRqLgDI4EStSMyO+AbOBlIcQ2IcSfEtQOU1OZn4FNoMsgcLxjgDSHzdSZEMcyXcesmPVdyaUBA1QXZuoyACajBuyw26jMT6dWR+WnMomUn6rCgAWkQ/8nIiGTwFLKOYm4brKR6rBRnpeuyyBQ16VWgZq5EtZYpge5AWYXZ8V0rrpON5ctmKZHs+JGVWEm7a4h+jzDZMdQxD0ZNWCAGToJwGR0f5bnpZNiF7oIwIlIjpioU5iqwgx9tKBOd9KEQGpUF2YCUNse2yAwMORNOg0YoLpQHxdgMmrAoPp/rH2AWKcI6zoHcKbYKM5KHuvXbhNML8jQRQBOhCUATE5VYSbHY3wIpJTUdQ2MTKwlC0VZqWSlOWIWALXtagCcWRSbFRFvqgICMFYBoPVfcyskC1WFmfQNeukaiK029PHA+pdksn5BKUCWBXCKU12YQdfAMD3u6F+Czv4h+jxeZgQGlGRBCEF1UewWUG1AgFYXJdsAqNpbG6MCUNvRT7bTYeo60KGo1qn/R9v7mVmUXM8+BIIAOvpjtoAmwhIAJmdGgXpwY5kMOhrQoGcVJ+NLkKnLAACjLqVkITPNQXF2WsxugKPt/cwqykw6DXjUAoq+/34pOdY5kJQCoLowk4EhH22uQcOuYQkAk6NprbEMgkc0AZCEL8HMwkzqu2KrjHa0vZ9p2Wlkpplh4XtkVOswB3SkrZ/qJLz30wvSEWLUhRUNHW7JkNeflAKgSqc5oImwBIDJmTGyFiAGAdDWT4pdJEU5vLFUF2Xi88uY8qIkqwsAlBYcy733DPto7HEnZf/THHbKc9M5HkMYdEsgpUYy9l+vIIiJsASAyclIdTAtOy0mLeBou4sZBRlJkwgrmJmaBRTDS1CbxAKgujCDlt5BBoa8k+8cguOdA0iZnAMgaFFw0d/75n7lP5+ZhO7Piny1cNGyAE5xqgszR/zY0aA04OSKgNHQ/MDR9r/HPUxH/1BSukBgtP/RasFH2tTvlrwCIDOmAbC5309mqj2pQkA1UkYWwxlnASSfU3QMw8PD1NfX4/FMnDc9NzeXvXv3xqlVJ+J0OqmsrCQlJbrFPLOnZfLvXc1IKSOeyPP5JbUdA1wyP7kWQWkUZqaSneaI2g2iWQ7JOgCOugEGWFCaE/HxoxFQydn/qsIMOvuH6HEPk5se+fvTPCCZWZx8E+Aa1YWZI0LcCJJeANTX15OdnU11dfWEN7mvr4/s7Ow4tkwhpaSjo4P6+npmzpwZ1TnmTMume6COjv4hiiLUZBq73Uk7CQZaKGjmyER2pGgDYLL2XwsCONzmiur4o239FGWlkhPDSuJEoq0AP9zmYsWM/IiPb+n3c05Fclq/AHOnZfHOkQ58fmNCQZPeBeTxeCgsLDSthBdCUFhYOKmFMhFzpqkH+FBr5IPAkSTXgEH1P5q+g3KBCEHSLYLTyHamUJbrjLr/yTwBDjCvRD37B1v6Ij520Ouj3S2Tuv9zS7IY9PoNKxCf9AIAMO3grxFr++bGIACOBjTHZJwE05hbkkVTj4deT+SL4Q61upien4EzxW5Ay+LDnGlZHIhiAAQ40u5K6gGwMj+DNIeNgy2RP/vHOwaQjAYSJCNzpimvxcHW6O7/ZEwJATDVKct1kplqj0oAHGx1ke10JOUkmMY87SWIYhA40NI3okUmK/NKsjnU6orYDdDhGqTdNcS8kvi7PvXCbhPMLs7iYBTP/oHA8zJ3WvL2X7P+o+l/OFgCIAkQQkTtBjnQ0sf8kmzTW0kToQ1gkboBhrx+jrb3MzeJB0BQbpBBr38kq2e4aANgMgsAUP2PxgW0v6UPweggmozkpqdQkpMWlfITDpYASBJmT8uK2AyUUrKvuY/5pck9AFTmp5OeYh8Z0MKltqMfr18mvQWguQEi7b/mNkp2ATC3JJvGHg99EboADzT3MS1DJLX7D5QFc8ggF1DSRwEF871nd7OnsTfkNp/Ph90e+YOwsDyH77z39HG333nnnRQVFfGlL30JgG9961uUlJTwxS9+MeJrTcTcadk8saUhIj94c6+HPo836QWAzaYsoEgF4NQZADU3QB/vXlgS9nH7W/rIcTooyUle9x+MavCH2/pZNj0v7OMOtPRRmZ38Ou6caVk8uqkOKfW/j8n/6ySYT37yk9x7770A+P1+Hn74YW699VbdrzPiC4xAC9zXrAbA+Uk+AIIaBCOdCD3Q4sImiLmYTKLJCUQCReoGONiirL9kdv9BdC5Az7CP2o5+KrKSf4ibW5LFwJCPTo/+oaBTygKYSFM3ah1AdXU1hYWFbN26lZaWFpYvX05hYWxFrEOxIKDF72vupSLMY/ZrAiDJLQBQg8ATWxoiWhB0oLmPqsLMpHcBQOSRQFJK9jf3ce3ScgNbFR+m56eT6rBFNBF6qNWFX0LlVBAAARdggyv6hIjjkaii8P8rhNgRqAf8khAiqZ/ST33qU/zjH//gnnvu4ROf+IQh16jMTyc3PYVdDaFdXKE40NxHSU4aeRnJlQc+FJofP5JBcG9z74jgTHbmByKBvGFmRW3u9dDr8U4J689htzGvJGtc924otOdkKlgAC8tz+NNtK5mZq78ik6hf52dSyiVSymXAauDbCWqHLrzvfe/jhRdeYOPGjVx55ZWGXEMIwaKKHHY39oR9zO7G3qjSB5iR08tzAdjVEF7/e9zDHOsYYFFFrpHNihuLK3MZ9PrD1oJ31qvfacr0vyKXnQ09YRdH2dXQizPFRmlmcru/ALLSHFy1qJTsVP37khABIKUMFuWZgHElb+JAamoql1xyCR/4wAeimmgOl9PLc9nX1Ic3jHjwgSEvB1v7WFo5NQaAkhwn07LTRga2ydgdEBSLp9AACITd/50NPdhtgoVlU0MBWFSRS497OOy04Dsbujm9PBe7LfkFgJEkbA5ACPFD4CNAD3DJBPvdAdwBUFJSQk1NzQnbc3Nz6eub3C3g8/nC2i8a/H4/b731Fvfee++41/B4PCe1PVJs3V6GfH4Ot/XjmORc+zt9+CXYuuuoqWmK6bpG4nK5wv5dyp1e3jnQSE1N96T7Pn90CIDe2l3UNBo/CETSj2jwS4nTDi9s2MO0/sOT7v/6Dg/lmYL1b62L+FpG9yUahnp8ADz80lucWTrxsOXzS7bXDXBRpQOXa9h0fYkWQ+6LlNKQD/AKsCvE5/ox+30D+F4451y5cqUcy549e076LhS9vb1h7Rcpu3fvljNnzpRf+cpXJtwv3HZOxKHWPlm1arX8/v0vTbrv3WsPy6pVq2Vrryfm6xrJmjVrwt73Vy8fkNVfXy37PMOT7vu5f26W5/341RhaFhmR9CNaPvjnt+R1v3tj0v38fr9c8f2X5Fcf3RbVdeLRl0jxDHvlnG8+J+/6995J993X1CurVq2WT2ypM2VfoiWWvgCbZIgx1TALQEp5eZi7Pgg8B3zHqLYYycKFCzly5EhcrjWzMJPsNAeHeyafCNxe30N5rpPi7OSOAQ9myfRcpFTunbNnTRxptauhZ8q4fzQWV+Ry79vHGPb5SZmguE9jj4eO/iGWTBH3H6jqYPNKssNyge2o7wZgcUUeYXrMTlkSFQU0N+jP64B9iWhHsmGzCVZW53OgyzfpvtvrulkawaKZZEAb0LfVdU+4X4drkNqOAZZMnzoDIMCSyjyGvP5Jo2G2B36fxZV5xjcqjiybnse2uu5JI6F21PeQleZIyhrY8SZRUUB3CSF2CSF2AFcAX0pQO5KOs2YW0OiSdPYPjbtPS6+H450DLJ+RF7+GxYGirDRmFWWy/mjnhPttrO0C4OyZBfFoVtzQ+rP+aMeE+2042kl6ip3Ty6fGBLDG2bMKcQ162dM0sQDcWNvJ8hl52KwJ4ElJVBTQjVLKRVKFgr5XStmQiHYkI2dVq0FgY+34g+Dbh9UAce6sori0KZ6cM7uQDUc7J9QCNxztJM1hY3FFXvwaFgem5TiZVZw5cn/H450jHaysyp/QTZSMnKMJwCPjP/ud/UPsa+7jnElchBaKqfWEnAIsrswlxaYGufF4+3AHOU4HC6eYBghw3mylBe6cYD3AhtoOls/II9Ux9R7vc2cVsrG2a1wB2D0wxP6Wviln/YASgDOLMnnnyPgCcEPAOjpn1tTrvxFMvTdkipPmsDMnz8YbB9tDbpdS8ubhds6aWTglY6A1ze7tcQaBtr5BdjX0cv7sqWf9gOq/a9DL7nHmAd481IGUcO7sqakBnzOrgA21nePWRnjzUAfpKfYpZ/0ZhSUAkpBl0xzsb+njeMfJ+eH3NvVR3+Xm8tOSswj8ZBRlpbGgNJvX97eF3L5mfysAl07R/p87uxAhRvs5llf3tpCfkcLyKOrnJgMXzCmmz+NlUwgXqJSSV/a2cOHcoilp/RnBlEoGx7+/Ds07Q25K93nBHkV3SxfD1XdNuEttbS3XXnstu3btAuDnP/85LpeL7373u5FfLwxWTLPz0D54aU8zn7pw1gnbXtjdjE3A5RGkDU42rlhYwu/WHKKtb/CkMNdX97ZQmuOcMitgx1KUlcaZVQW8sKuZL18+74RtXp+f1/a3cumCaVPS+gO4eH4xqQ4bL+xuPikUeHdjL009Hr7y7nnjHG0xFktMJiHFGTZOL8/hiS0NJ+RG8fslq7c3ckZ1AUVJXAJyMq5ZUoZfKmEXTFf/EGv2tXHVotKkT4E8EVcvLmVfc99JFeLWHWyne2CYKxaWJqhlxpOZ5uBdc4t5YVfzSW6gZ3c0YrcJLl0wNa0/I5haFsAEmrrboHTQieLDZ83gf57axda6blYEzP11h9o50t7PFy6bk+DWGcv8kmzml2Tz4Prj3Hb2jJHB/qltDQz5/HzgjOkJbqGxvGdxGT98bi//XH/shBToD244TlFWGpdNUfeXxo0rKnhlbwuv7m3hitOVsBv0+vjXpnouP20ahVNY+dEbywLQAYfDgd8/GpXh8XgMv+YNyyvITnPw+9cOAcr/+Yc1hyjKSuOaxWWGXz+RCCH45AUz2dvUy7rAZPiQ1889b9aytDJ3SkY/BTMtx8l7l5bz6MY6ugfUepDDbS5e29fKzWdUTrnwz7G8e2EJFXnp/HXd0REL+MktDXT2D3HbOVUJbl1yMbWflDhRUlJCa2srHR0dDA4Osnr1asOvmZXm4HOXzuHVfa38seYwv3zlIOuPdvKly+eS5kj+AiiTcf3yciry0vn+6j30eob53WsHOd45wJcunzv5wVOAz1w0C4/Xz/ee3cOg18f/PLkLp8PGJy+YmeimGY7DbuMzF81iQ20nj26qo7Hbzc9f2s/yGXlcMGdqRn8ZxdRyASWIlJQUvv3tb3P22Wczc+ZMFixYEJfrfuqCmWyq7eQnL6hMGtcuKePWs2bE5dqJJs1h56c3LeGjf9/AGT94hSGvnxtXVHLJ/Knt/tBYUJrDFy+dyy9fOcBzO5sY8vr5+c1Lp/TcTzC3nl3FC7uaWfX4TlLtNlLsgrvev2RKz/0YgSUAdOKLX/yi7oXgJ8Nht/GX289g3aF2HDbBebMLT6kX4Pw5RTx0xzk8vrme+aXZfOTc6lOq/1+8bA4V+elsONrBFQtLp3Tk11jsNsFfP3oGf113lJZeDx89r3qkdrBF+FgCIMmx2QQXzStOdDMSxpnVBZxZfWqu+hRCcNPKSm5aWZnopiSEjFQHX7zs1HD5GYU1B2BhYWFxijIlBEBwLLwZMXv7LCwsTk2SXgA4nU46OjpMO8hKKeno6MDpdCa6KRYWFhYnkPRzAJWVldTX19PWFjo3jIbH40nYIOx0OqmsPDX9tBYWFuYl6QVASkoKM2dOHvtcU1PD8uXL49AiCwsLi+Qg6V1AFhYWFhbRYQkACwsLi1MUSwBYWFhYnKIIs0bPhEII0QYci/LwIiB0Ga3kw+qL+Zgq/QCrL2Yllr5USSlPWjGaVAIgFoQQm6SUZyS6HXpg9cV8TJV+gNUXs2JEXywXkIWFhcUpiiUALCwsLE5RTiUB8JdEN0BHrL6Yj6nSD7D6YlZ078spMwdgYWFhYXEip5IFYGFhYWERhCUALCwsLE5RTgkBIIS4SgixXwhxSAjx9US3JxaEELVCiJ1CiG1CiE2Jbk+4CCH+LoRoFULsCvquQAjxshDiYODf/ES2MVzG6ct3hRANgfuyTQhxTSLbGA5CiOlCiDVCiL1CiN1CiC8Fvk+6+zJBX5LxvjiFEBuEENsDffle4Hvd78uUnwMQQtiBA8C7gXpgI/BhKeWehDYsSoQQtcAZUsqkWtwihHgX4ALuk1IuCnz3U6BTSnlXQDDnSylXJbKd4TBOX74LuKSUP09k2yJBCFEGlEkptwghsoHNwA3Ax0iy+zJBXz5A8t0XAWRKKV1CiBTgDeBLwPvR+b6cChbAWcAhKeURKeUQ8DBwfYLbdMohpVwLdI75+nrg3sD/70W9sKZnnL4kHVLKJinllsD/+4C9QAVJeF8m6EvSIRWuwJ8pgY/EgPtyKgiACqAu6O96kvTBCCCBl4QQm4UQdyS6MTFSIqVsAvUCA9MS3J5Y+bwQYkfARWR6t0kwQohqYDmwniS/L2P6Akl4X4QQdiHENqAVeFlKach9ORUEgAjxXTL7vc6XUq4ArgY+F3BHWCSePwKzgWVAE/B/CW1NBAghsoDHgS9LKXsT3Z5YCNGXpLwvUkqflHIZUAmcJYRYZMR1TgUBUA9MD/q7EmhMUFtiRkrZGPi3FXgS5eJKVloCvlvNh9ua4PZEjZSyJfDS+oG7SZL7EvAxPw78U0r5RODrpLwvofqSrPdFQ0rZDdQAV2HAfTkVBMBGYK4QYqYQIhX4EPBMgtsUFUKIzMAEF0KITOAKYNfER5maZ4CPBv7/UeDpBLYlJrQXM8D7SIL7Ephs/BuwV0r5i6BNSXdfxutLkt6XYiFEXuD/6cDlwD4MuC9TPgoIIBD69SvADvxdSvnDxLYoOoQQs1BaP6hyng8mS1+EEA8BF6NS2rYA3wGeAh4FZgDHgZullKafXB2nLxej3AwSqAU+o/lrzYoQ4gJgHbAT8Ae+/ibKd55U92WCvnyY5LsvS1CTvHaUkv6olPL7QohCdL4vp4QAsLCwsLA4mVPBBWRhYWFhEQJLAFhYWFicolgCwMLCwuIUxRIAFhYWFqcolgCwsLCwOEWxBIDFKYEQIk8I8Z9Bf5cLIR4z6FopQojNRpzbwkJPLAFgcaqQB4wIACllo5TyJoOudQHwlkHntrDQDUsAWJwq3AXMDuSE/5kQolrL5y+E+JgQ4ikhxLNCiKNCiM8LIb4ihNgqhHhHCFEQ2G+2EOKFQCK+dUKIBeNc6yrg38FfBJJ7/UMIsUuoeg7/b6JzCiFKhBBPBnLCbxdCnGfYL2NxyuJIdAMsLOLE14FFgQRbWsbIYBahMkg6gUPAKinlciHEL4GPoFaS/wX4DynlQSHE2cAfgEtDXOsS4HtjvlsGVATVD8gLfD/eOX8DvC6lfF+gpkVWdN22sBgfSwBYWCjWBPLI9wkheoBnA9/vBJYEskyeB/xLpZ0BIG3sSYQQ5aiiHQNjNh0BZgkhfgs8h0rpPdE5L0UJHqSUPqAn9i5aWJyIJQAsLBSDQf/3B/3tR70nNqBbsyAm4GrgxbFfSim7hBBLgSuBz6EqVX05zHNaWBiCNQdgcarQB2RHe3Agt/xRIcTNoLJPBgb0sZzk/w/sXwTYpJSPA3cCKyY556vAZwPf24UQOdG23cJiPCwBYHFKIKXsAN4MTML+LMrT3Ap8UgixHdjNmNKiAV/9XCnlvhDHVgA1gSpP/wC+Mck5vwRcIoTYiapve3qUbbawGBcrG6iFhU4EUhLfJqX8j0S3xcIiHCwBYGFhYXGKYrmALCwsLE5RLAFgYWFhcYpiCQALCwuLUxRLAFhYWFicolgCwMLCwuIUxRIAFhYWFqco/x/aEkzownniVQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def first_order(K=1, tau=1, tfinal=1, u=lambda t: 1):\n", " model = ConcreteModel()\n", " model.t = ContinuousSet(bounds=(0, tfinal))\n", " model.y = Var(model.t)\n", " model.dydt = DerivativeVar(model.y)\n", " model.y[0].fix(0)\n", " @model.Constraint(model.t)\n", " def ode(m, t):\n", " return tau*model.dydt[t] + model.y[t] == K*u(t)\n", " \n", " tsim, profiles = Simulator(model, package='scipy').simulate(numpoints=1000)\n", "\n", " fig, ax = plt.subplots(1, 1)\n", " ax.plot(tsim, profiles, label='y')\n", " ax.plot(tsim, [u(t) for t in tsim], label='u')\n", " ax.set_xlabel('time / sec')\n", " ax.set_ylabel('response')\n", " ax.set_title('Response of a linear first-order ODE')\n", " ax.legend()\n", " ax.grid(True)\n", " \n", "first_order(5, 1, 30, sin)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "XEubS-8Kf24e", "pycharm": {} }, "source": [ "## Analytical approximation to a step input\n", "\n", "The math functions supported by Pyomo are limited to that are imported to the standard arithmetic operations of Python (\\*, /, \\**, \\*=, /=, \\**=) and a particular set of nonlinear functions from the Python `math` library. Simulating the response of a system to a discontinuous step input, for example, requires an approximate.\n", "\n", "An infinitely differentiable approximation to a step change is given by the *Butterworth function* $b_n(t)$\n", "\n", "$$ b_n(t) = \\frac{1}{1 + (\\frac{t}{c})^n} $$\n", "\n", "where $n$ is the order of a approximation, and $c$ is value of $t$ where the step change occurs." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "colab_type": "code", "executionInfo": { "elapsed": 6178, "status": "ok", "timestamp": 1557657355055, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg", "userId": "09038942003589296665" }, "user_tz": 240 }, "id": "XN4fHvsPisM6", "outputId": "d3479b3f-0680-49ba-bc8d-f9a102fc460a", "pycharm": {} }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxHElEQVR4nO3de3xcZZ348c83c8k9TZM2adP7hVJKBUord5cWQUAUV1dFV2QRFf2tLujP9cqKuK6ru7L+vK3XVVERisiiCAgitiBIgRaKXFKglN5o06RJ02SSTJKZ+f7+OGfKEHKZJHNyZs5836/XvDIz55znfJ85mW+ePOc5zxFVxRhjTPCU+B2AMcYYb1iCN8aYgLIEb4wxAWUJ3hhjAsoSvDHGBJQleGOMCShL8CaviMjpIvK8iMRE5G8nWdY1InK9+3y+W2YoJ4FOkoj8m4gcFJGWfIttOCKyUUQ+4HccZnwswftIRHaKSJ/75W4RketEpMrvuHz2r8B3VLVKVX+Tq0JVdbdbZjJXZU6UiMwDPgGsUNVZk4lNRC4VkQdyH6W3RKRURL4iIrvd78DzIvJJEZGMdTaKSFxEukWkS0S2iMhnRKQ0Y51rRGTQ/Q6lH52+VCoPWYL335tVtQo4AVgFfNbfcHy3AHja7yByRUTCw7y9AGhX1dYsthcRmbLvaa73N0L9AW4GXg+8EagG3gtcDnxzyHofVdVqYDbOH8V3AXdm/iEAbnL/QKYftbmKv9BZgs8TqtoC3I2T6AEQkVNE5C8i0ikiT4jI2oxll4rIDrd186KIvCfj/QdF5NsiclhEtonI6zO2axKR20SkQ0S2i8gHM5ZdIyK/EpGfu+U+LSJrMpZ/WkRecpc9my5XRErcltULItLullE3Ul1F5IPuvjvcWJrc918AFgO/c1tipcNsm95Pt4g8IyJvzebzFZGFIqLphOO2Dr/kflbdIvIHEZmR5Wf/PhFpdrfbISIfyli2VkT2up9VC/DTIXGcDdwDNLl1vG6E2L4sIg8CvcDi4Y63iBwDfB84dayWq4icJiKPur8Tj4rIaRnLhtvfOe7vzmER+Q4gQ8q7zP0MDonI3SKyIGOZishHROR54PlhYnk98Abg71T1KVVNqOom4GLgIyKydOg2qtqjqhuBC4FTgQtGqqvJoKr28OkB7ATOdp/PBZ4Evum+ngO047RwSoBz3NczgUqgCzjaXXc2cKz7/FIgAXwciAAXAYeBOnf5fcB3gTKcPyZtwOvdZdcAcXefIeArwCZ32dHAHqDJfb0QWOI+/xiwya1DKfAD4MYR6nwWcBA40V3328D9w30mI2z/DqDJ/UwuAnqA2SOsew1wfUa8CoTd1xuBF4BlQLn7+qtjffbu8guAJThJ70ycpHiiu2yt+/n/h1u/8mHiWgvszXg9XGy7gWOBMDBtjOP9wBi/Z3XAIZxWchh4t/u6foT9zXT393ac36GPu3X6gLv+3wLbgWPc9f8F+EvG/hTnj1jdCPX/KnDfCLHuAj6UEdcHhlnnfuA/hh5je7z6YS14//1GRLpxkmcr8AX3/YuBO1X1TlVNqeo9wGacpAOQAlaKSLmq7lfVzG6NVuAbqjqoqjcBzwIXiNP3ewbwaVWNq+pW4H9wvvhpD7j7TAK/AI5330/iJKwVIhJR1Z2q+oK77EPAVaq6V1X7cb50bx/h3/P3AD9R1cfcdT+L0wJdmM2Hpao3q+o+9zO5CaeFeFI22w7jp6r6nKr2Ab/i5f+eRv3sVfUOVX1BHfcBfwBel1FuCviCqva7ZU/Edar6tKomcJLraMd7LBcAz6vqL9RpLd8IbAPePML+zgeeUdVfq+og8A2gJWPdDwFfUdVmd/1/B07IbMW7yztGqP8MYP8Ise53l49mH84fj7R3uv9ppR8bxti+aFiC99/fqtPHuBZYzsu/3AuAd2T+4uIk59mq2oPTev0wsF9E7hCR5RllvqRu88a1C6fV2wR0qGr3kGVzMl5nfpF7gTIRCavqdpyW+jVAq4isT3etuLHemhFnM84fhMZh6tvk7hMAVY3htI7nDLPuq4jIJSKyNWNfKxk7IYxkaF3TJ7hH/OzdGM4XkU1uF1MnTuLPjKFNVeMTjCltT/pJFsf7CHl5RE5MRGLu26/4zF1Dj/uejOdNQ/avQ5YvAL6Z8dl04Pw3M1J5Qx3E/SyHMdtdPpo57j7TfqWqtRmPdWNsXzQswecJtyV4HXCt+9Ye4BdDfnErVfWr7vp3q+o5OF+IbcCPMoqbI/KKk1DzcVo9+4A6EakesuylLGO8QVXPwPmCK043RDrW84fEWqaqw5W7z90eABGpBOqzicFtIf4I+ChO90It8BRD+odzYMTP3j0vcAvOcWp0Y7hzSAy5mKL1FWWMcryHrpcekVOlzsl7GPKZu4Ye98xy9gPz0i/c36V5Gcv34HSjZH4+5ar6l5HiH+KPwMnuf5RHiMhJ7n7+NNKG7jargT+PUr5xWYLPL98AzhGRE4DrgTeLyLkiEhKRMvcE3lwRaRSRC93k2A/EcFrMaQ3AFSISEZF34PSV3qmqe4C/AF9xyzsOeD/wy7ECE5GjReQsN8HFgb6MfX4f+HL6X3QRmSkibxmhqBuA94nICW5Z/w48rKo7s/h8KnESR5u7n/fhtOBzbcTPHojidFW1AQkROR/nhKFnxjjeB4C5IhIdpYg7gWUi8vciEhaRi4AVwO0jrH8HcKyIvM3tZrsCmJWx/PvAZ0XkWDe+ae7vWVZU9Y/AvcAtInKs+xmfgvN7+D1VHe7EbIWInAn8FnjErZMZgyX4PKKqbcDPgc+7yfgtwOdwkske4JM4x6wEZ8jYPpx/Vc8E/jGjqIeBo3D+1f0y8HZVbXeXvRvnpN4+4FacvuJ7sgivFOfk2EGcro0GNzZwhrbdBvzBPZ+wCTh5hDreC3wepxW8H+dk5buy2D+q+gzwX8BDOIntNcCD2Ww7HqN99m731hU4ffaHgL/HqbuXRjvef8IZVtoiIsN2bbjH/k1uGe3Ap4A3qepI6x/EOZn9VXf9o8j4nFX1Vpz/3taLSBfOf1Hnj7NOfwdsAO7C+YN1PfBj4J+GrPcd93fqAE4D6BbgPFVNZaxzkbxyHHxMRBrGGU8gySu7ak2hE5FLcUYenOF3LMYYf1kL3hhjAsoSvDHGBJR10RhjTEBZC94YYwJqpImAfDFjxgxduHDhhLbt6emhsrIytwH5JCh1CUo9wOqSj4JSD5hcXbZs2XJQVWcOtyyvEvzChQvZvHnzhLbduHEja9euzW1APglKXYJSD7C65KOg1AMmVxcRGXqV8hHWRWOMMQFlCd4YYwLKErwxxgRUXvXBG2OMHwYHB9m7dy/x+GQnAZ2YadOm0dzcPOo6ZWVlzJ07l0gkknW5luCNMUVv7969VFdXs3DhQl45EevU6O7uprq6esTlqkp7ezt79+5l0aJFWZfraYIXkZ1AN87MdwlVXTP6FsYYM/Xi8bhvyT0bIkJ9fT1tbW3j2m4qWvDrRpq1zhhj8kW+Jve0icRnXTQFKpFM0TuYpLc/Se9Agt6BJP2JJINJJZFUBlMpEkklkUwxmFIGEykSqRSDSfd+jUAq5f7U9L15QVH3NaTcaSzS66Xfc+66qURCJYRDJURCwoyqUo6fV8uiGcG48MSYIPB0LhoReRFnzmwFfqCqPxxmncuBywEaGxtXr1+/fkL7isViVFVVjb1iHutLKAf7lL0dvcQo5VBc6R5QYoPOz+4BpWdQiSchkRq7PD+c2hTispWlREokEMckzeqSf3JZj2nTprF06dKclDURyWSSUCg05nrbt2/n8OHDr3hv3bp1W0bq/va6BX+6qu5zJ9+/R0S2qer9mSu4Sf+HAGvWrNGJXs1VaFe17T3Uy6M7O2je303z/i6a93dzMNbvLhVggEhIqKuMMr0iSuOMKMsrotRWRKgsDVMRCVNZGqI8GqIyGqY8GqI0XELUbVWHQ0KkxP0Zkpdb2yWCiFAivPwTQUqcvZaIIOL8JON15jIRQVVJppTBpDKQSHGgO85tW/fxnQ3bOW7pLD77xmMK7piMxuqSf3JZj+bm5lFPcnptrJOsaWVlZaxatSrrcj1N8Kq6z/3ZKiK3AicB94++VTDFB5P8+fmD3PNMCw/taGdPh3Oz+Wi4hGWNVaw7eiZLGqqYO72c1h3NXHj2GdRXRvO2X1BECIeEcAjKoyGmVUT453OPpqUrzk8f3MmHzlzid4jGFIwvfelLzJkzhyuvvBKAq666isbGRq644opJletZgnfvH1miqt3u8zcA/+rV/vLV1j2d/Pyhnfzh6QPE+hPUlIU5dUk9l52+iJMX1bOssYpw6JXXm23seI4ZVaU+RTw5H3jdIn69ZS+/3foS2Q/mMiZ/fPF3T/PMvq6clrmiqYYvvPnYEZdfcsklXHLJJVx55ZWkUinWr1/PI488Mun9etmCbwRudVugYeAGVb3Lw/3lDVXl3uZWvrNhO1v3dFJVGuaNr5nFG18zm9OXziASCu4FxMtn1bCgvoIHt7ezaIHf0RhTGBYsWEB9fT2PP/44Bw4cYNWqVdTX10+6XM8SvKruAI73qvx89de9nXzp9md4dOchFtRX8MULj+XvVs+lqrR4BiydvKiOu58+wN/Pz/6KO2PyxWgtbS994AMf4LrrrqOlpYXLLrssJ2UWT9bxWH8iyTf++Dw/uO8F6qtK+fJbV/LONfMC3VofyZoFdfxq815ae+3Xy5hsvfWtb+Xqq69mcHCQG264ISdl2jcwB/Z09HL5L7bQvL+Li9bM46o3HUNNWfG2Xo9qdIau7Yvl6VhOY/JQNBpl3bp11NbWZjVkMhuW4CfpkRc7+NAvNpNMKT+5dA1nLW/0OyTfLWlwE3yPJXhjspVKpdi0aRM333xzzsosvv6DHLr/uTYu+cnDTK+M8tuPnmHJ3VVTFmFWTRn7YnZDd2OysW3bNpYuXcrrX/96jjrqqJyVay34Cfrz82184GebWdJQxfXvP4n6Ah3W6JUF9RW0Her0OwxjCsLy5cvZsWNHzsu1FvwEPNvSzf+5/jEWz6zkxg+ebMl9GHNqy+mIWwveGD9Zgh+ng7F+LrvuUSpLQ/z0fa+ltiLqd0h5qam2nEP9zmRnxhh/WIIfB1Xln29+grZYPz/+h9cye1q53yHlrabaclIKrd39Y69sjPGEJfhxuO4vO9n4bBv/csExrJwzze9w8lpTbRkA+zr7fI7EmOJlCT5LO9pifOX32zj7mAbee4pdgz+W9H83LV3+3OPSGGMJPiuqytW/fZrScAn//rbX5O0Mj/lkRpVzbuKgddEY4xtL8Fm448n9PLD9IJ8892gaqsv8Dqcg1FZEEeBgbMDvUIwpCDt37mTlypVHXl977bVcc801kyrTxsGPYSCR4j/u2sYxs2t4z8nWNZOtUIlQHZWMm5gYUyB+/xloeTK3Zc56DZz/1dyWmQVrwY/hpkd3s6ejj0+ddzShEuuaGY9ppZbgjfGTteBHER9M8q0/beekhXWsXTbT73AKTk3UumhMAfKhpQ0QDodJpV6+biQen/wABWvBj+LXW/bS1t3Px89ZZidWJ6DGWvDGZK2xsZHW1lba29vp7+/n9ttvn3SZ1oIfQSql/PiBFzl+7jROWVzndzgFaVpUeLytH1W1P5DGjCESiXD11Vdz8skns2jRIpYvXz7pMi3Bj+CPzQd48WAP3373KktOE1QTFeKDKXoGkkV1RytjJuqKK66Y9I22M1kXzQh+9tBO5tSWc/7KWX6HUrBqSp0/jDYW3hh/WIIfxu72Xh7c3s5Fr51HuAhvuZcrNVEnwbf3WII3xg+WvYZx85Y9lAi8ffVcv0MpaFVugj/UM+hzJMaMTTW/p7eeSHyW4IdIJFPcvHkvZy6bSVOtzRY5GVURN8H32lBJk9/Kyspob2/P2ySvqrS3t1NWNr4r6e3M1xCbdnTQ0hXnC29e4XcoBS+d4Dt7rQVv8tvcuXPZu3cvbW1tvuw/Ho+PmbzLysqYO3d8vQqW4Ie448l9VEZDrFve4HcoBa88DOESsRa8yXuRSIRFixb5tv+NGzeyatWqnJdrXTQZBpMp7nqqhbNXNFIWCfkdTsETEWorIhyyFrwxvrAEn+GhF9o51DvIG18z2+9QAqO2IkqnteCN8YUl+Ax3PrmfymiIM23emZyZXhGxLhpjfGIJ3qWq3LutlbVHN1j3TA45LXjrojHGD5bgXU/v66Ktu5+1R1vrPZemV0To6LEWvDF+sATvuu85Z3jUmZbgc2q624LP1/HFxgSZJXjXhm2tvGbONLslX47VVkQZSKboHUj6HYoxRcfzBC8iIRF5XEQmP7mxRzp7B3hs9yHWWes956ZXRAC7mtUYP0xFC/5KoHkK9jNhm3a0k1L4Gxs9k3O1FVHArmY1xg+eJngRmQtcAPyPl/uZrE07OiiLlHDc3Fq/Qwkca8Eb4x/x8uSXiPwa+ApQDfyzqr5pmHUuBy4HaGxsXL1+/foJ7SsWi1FVVTWhbT//YB81Ufjka/NjcrHJ1CWfxGIxDmsFVz3Yx4ePL+WU2YU7M0ZQjgkEpy5BqQdMri7r1q3boqprhl2oqp48gDcB33WfrwVuH2ub1atX60Rt2LBhQtsd6unXhZ+5Xb/1x+cmvO9cm2hd8s2GDRu0tSuuCz59u/7sLy/6Hc6kBOWYqAanLkGph+rk6gJs1hFyqpddNKcDF4rITmA9cJaIXO/h/ibk4Rc7UIVTltT7HUog1aa7aGxOeGOmnGcJXlU/q6pzVXUh8C7gT6p6sVf7m6iHd3RQGi7huLnT/A4lkCKhEqpLw9YHb4wPin4c/OZdHZwwr5bSsE1P4JXayohNOGaMD6YkwavqRh3mBKvf4oNJmvd3sWr+dL9DCbTpFVGbMtgYHxR1C/6Z/V0MJpUT5tX6HUqg2ZTBxvijqBP81t2dAKyaX+trHEE33W76YYwvijvB7+lk9rQyGmts/hkvOV001oI3ZqoVfYI/3q5e9VxtRYTueILBZMrvUIwpKkWb4Ntj/ezu6OUE657xXF2lMx/N4T7rpjFmKhVtgn/ypcMA1oKfAukJxw7ZjT+MmVJFm+Cb93cDsGJ2jc+RBF+dm+Dtzk7GTK0iTvBdzKktZ5p7Kb3xzvTK9IyS1kVjzFQq6gR/zOxqv8MoCuk+eBtJY8zUKsoEHx9MsuNgD8dY98yUmG5dNMb4oigT/PMHYiRTagl+ipRFQpRHQnaS1ZgpVpQJvnl/F4Al+ClUVxmlw7pojJlSRZngn9nfRUU0xIK6Cr9DKRrTKyN2X1ZjplhRJvhtLV0cPauakhLxO5SiMb0ian3wxkyxokzw21tjLGuwETRTqa7S5qMxZqoVXYLv7B3gYGyAJQ2VfodSVKwFb8zUK7oE/0JbDIClDcG4G3uhmF4RtQnHjJlixZfgW3sAWDLTEvxUqjtyNau14o2ZKkWX4Le3xYiGS5g73UbQTKXp7tWsNpLGmKlTdAn+hdYYi2dUErIRNFPKJhwzZuoVXYLf3haz7hkf2JTBxky9okrw8cEkezp6WWInWKdcesIxu5rVmKlTVAl+Z3sPKYUlM22I5FSrdadltha8MVOnqBL8jjYbQeOXskiIymjI5oQ3ZgoVVYLf1d4LwIJ6G0Hjh+mVUWvBGzOFiirB7+7ooa4ySnWZ3cXJD9MrbEZJY6ZSkSX4XubbDJK+sRa8MVOrqBL8rvZe657xUV1FhHZL8MZMmaJJ8AOJFPs6+2wOeB/VV5XahU7GTKGiSfAvdfaRUphnCd43M6pK6R1I0juQ8DsUY4qCZwleRMpE5BEReUJEnhaRL3q1r2zsaneGSC6otzHwfqmvci52OthtrXhjpoKXLfh+4CxVPR44AThPRE7xcH+j2t1hQyT9NrOqFICDPf0+R2JMcQh7VbCqKhBzX0bch3q1v7Hsau+lLFJCQ3WpXyEUvRnpBN9tCd6YqSBOHvaocJEQsAVYCvy3qn56mHUuBy4HaGxsXL1+/foJ7SsWi1FVNfIVqt98LE5rb4ovn5H/Lfix6lIohtajvS/FJ+7r49Jjo6ydV1jXIgTlmEBw6hKUesDk6rJu3botqrpm2IWqOuYDqAA+D/zIfX0U8KZstnXXrwU2ACtHW2/16tU6URs2bBh1+Ru+fp++/7pHJlz+VBqrLoViaD3igwld8Onb9Vt/fM6fgCYhKMdENTh1CUo9VCdXF2CzjpBTs+2D/ylOn/qp7uu9wL9l+xdGVTuBjcB52W6Ta/s6++wmHz4rDYeoKQtzMGZdNMZMhWwT/BJV/U9gEEBV+4BR75ghIjNFpNZ9Xg6cDWybeKgT1xUfpLs/QVNtmR+7NxlmVJVy0MbCGzMlsj3JOuAmaQUQkSU4LfrRzAZ+5vbDlwC/UtXbJxzpJOzvjAPQVFvux+5NhhlVpXaS1Zgpkm2C/wJwFzBPRH4JnA5cOtoGqvpXYNWkosuRfZ19gCX4fDCjOsqzLd1+h2FMUcgqwavqPSLyGHAKTtfMlap60NPIcuglN8HPsQTvu/rKUtp72v0Ow5iikFUfvIicDsRV9Q6cETGfE5EFXgaWS/s6+wiXyJFx2MY/M6pK6ewdZDCZ8jsUYwIv25Os3wN6ReR44JPALuDnnkWVY/s6+5g1rYxQyajnhc0UmFHtTFfQHrMTrcZ4LdsEn3DHW74F+JaqfhOo9i6s3NrXGbf+9zxRX+lezWpDJY3xXLYJvltEPgtcDNzhjowpmEsRX+rss/73PDHTbcFbgjfGe9km+ItwhkW+X1VbgDnA1zyLKoeSKaWlK25j4PPEyy1466IxxmvZjqJpAb6e8Xo3BdIH39odJ5lS66LJEw01ToJv7Y77HIkxwZftKJq3icjzInJYRLpEpFtEurwOLhdsDHx+qYiGqS4L09plXTTGeC3bC53+E3izqjZ7GYwXXkpfxTrNEny+aKwp40CXteCN8Vq2ffAHCjG5A7S6iWRWjfXB54vGmlJL8MZMgWxb8JtF5CbgN2TMQaOq/+tFULnU2t1PabiEmnLP7m1ixqmxuoyHX+zwOwxjAi/brFcD9AJvyHhPgbxP8Ae64jTWlCFiFznli4aaMlq746RSSoldfGaMZ7IdRfM+rwPxipPgbYqCfDKrppTBpHKod4B6mz7CGM9kO4pmrojcKiKtInJARG4RkbleB5cLrV39NFRb/3s+aXTPhxywkTTGeGo8d3S6DWjCucjpd+57ea+1u//I2GuTHxrSCd7GwhvjqWwT/ExV/amqJtzHdcBMD+PKiVh/glh/4kiL0eSHdJdZq42kMcZT2Sb4gyJysYiE3MfFQN5P6p1OINYHn19mVjvHw7pojPFWtgn+MuCdQIv7eLv7Xl5LJ5BG64PPK6XhEHWVUVqsBW+Mp7IdRbMbuNDjWHIuPd9Jg3XR5J2G6lLrojHGY9mOolksIr8TkTZ3JM1vRWSx18FNVnq+EzvJmn9mTSuzLhpjPJZtF80NwK+A2TgjaW4GbvQqqFw50BWnPBKiutSuYs03jdVl1kVjjMeyTfCiqr/IGEVzPc6VrHntQHc/jTWldhVrHppdW0Zbdz/9iaTfoRgTWNkm+A0i8hkRWSgiC0TkUzh3dqoTkTovA5yMA11x63/PU+npm1sOWyveGK9k23dxkfvzQ0PevwynJZ+X/fFt3f0c21TjdxhmGHPdBP9SZx8L6it9jsaYYMp2FM0irwPxwsHu/iNjrk1+Sbfg93VaC94Yr2Q7iuYdIlLtPv8XEflfEVnlbWiTEx9M0t2fYIZNZpWXZk1zus5eOtTncyTGBFe2ffCfV9VuETkDOBf4GfB978KavI4e56bO9ZVRnyMxwymLhJhZXXrklorGmNzLNsGnhzpcAHxPVX8L5HXmPBhzxljbdLT5q6m2nH2HLcEb45VsE/xLIvIDnOkK7hSR0nFs64v2mNOCn1GV13+Hitrc2nLrojHGQ9km6XcCdwPnqWonUAd80qugciHdgrc++PzVVFvGS519qOb9JRXGFKSsEryq9gKtwBnuWwngea+CyoX2dB+8teDzVlNtOf2J1JFjZYzJrWxH0XwB+DTwWfetCHD9GNvME5ENItIsIk+LyJWTC3V82mP9lEdCVERtmoJ8NefIUEnrpjHGC9l20bwVZzbJHgBV3QdUj7FNAviEqh4DnAJ8RERWTDTQ8WqPDVjrPc81WYI3xlPZJvgBdTpKFUBExrz0UFX3q+pj7vNuoBnndn9Toi3Wb/3veW7udCfB77UTrcZ4QsY6wSXOTF2fx0nO5wBfwZmi4AZV/XZWOxFZCNwPrFTVriHLLgcuB2hsbFy9fv36cVbBEYvFqKqqOvL66gf7qCsTPra68OaiGVqXQjVWPVSVf7y3l9Oawrx3RX7/MQ7KMYHg1CUo9YDJ1WXdunVbVHXNsAtVdcwH8BhOcv8acC1wTjbbudtWAVuAt4217urVq3WiNmzY8IrXJ335Hv3UzU9MuDw/Da1LocqmHhd86359748f9j6YSQrKMVENTl2CUg/VydUF2Kwj5NRsz0A+BHSq6riGRopIBLgF+KWq/u94tp0MVbU++AKxoK6Sp/cd9jsMYwIp2z74dcBDIvKCiPw1/RhtA7dr58dAs6p+fbKBjkdXX4JESu0q1gIwv76CvYf6SCRTfodiTOBk24I/fwJlnw68F3hSRLa6731OVe+cQFnj0nbkIidrwee7hfUVJFLK/sNx5tVV+B2OMYGS7XTBu8ZbsKo+APhyK6V2u4q1YMyvcwZk7WrvtQRvTI7l9XwyE2VXsRaOBfVOUt/Z3uNzJMYET6ATfJ1NFZz3ZtWUEQ2XsLuj1+9QjAmcQCb4Q26Cry23BJ/vSkqE+XUV7LIWvDE5F8wE3ztAVWmYaDiQ1QucBXUV7Gq3FrwxuRbIDNjZO0htRcTvMEyW5tdXsLuj16YNNibHApngD/UOML3CumcKxZKZVfQOJNl/2G7AbUwuBTTBDzLdTrAWjKUNzhwc21tjPkdiTLAEMsF39g4w3bpoCoYleGO8EcgEf6jHumgKSX1llNqKCNvbLMEbk0uBS/CJZIqueMJOshYQEWHpzCprwRuTY4FL8J19gwDWgi8wS2ZW8YIleGNyKngJvte9yMla8AVlaUMV7T0DRy5SM8ZMXuAS/KFea8EXoiMnWq0f3picCV6Cd1uAluALi42kMSb3ApfgO90WvHXRFJY5teVUREM829LtdyjGBEbgEvwhtw/eLnQqLCUlwvJZ1TTv7xp7ZWNMVgKX4Dt6B4iEhMpoyO9QzDgdM7uGZ/Z32Zw0xuRI4BJ8Z88gtRVRnFvCmkKyoqmG7niCvYf6/A7FmEAIXII/ZNMUFKwVs2sArJvGmBwJXIJ3pgq2/vdCdPSsakTgGUvwxuRE4BK8teALV0U0zKIZldaCNyZHApjgB20MfAFLn2g1xkxeoBK8qtLVN8g0a8EXrGObatjT0XdkygljzMQFKsH3J1IMJFPUlFmCL1QnzK0F4Im9h/0NxJgACFSCP+zOJDmt3BJ8oTpuXi0i8PjuQ36HYkzBC1SC73ITfI0l+IJVVRpmWUM1W/d0+h2KMQUvUAk+3YKvKQv7HImZjBPm1bJ1T6dd0WrMJAUqwXfFrYsmCE6YX0tn7yA723v9DsWYghasBN+XAKyLptCtml8LwNY91g9vzGQEKsHbSdZgOKqhmqrSMJt3WoI3ZjI8S/Ai8hMRaRWRp7zax1BHTrLaMMmCFioR1iyczsMvdvgdijEFzcsW/HXAeR6W/yqH+wYpj4SIhgP1j0lROmVxPdtbY7R19/sdijEFy7PhJqp6v4gs9Kr8V3jiJmbtf4qjYpW8O9oFWzNbfkOmDX7FNMIy9vuvWjba+xMp79XbzGh7Cpq7c1be+OPLZpuMZaEINKyAijpy5dTF9QA8/GI7bzquKWflGlNMxMuhaG6Cv11VV46yzuXA5QCNjY2r169fP+79vO7+dxJKWUvPT0oJLbPW8vxRHyYVKgUgFotRVVU1ofKSKeUj9/ZyWlOYS44tzWWoEzKZuuSboNQlKPWAydVl3bp1W1R1zXDLfB8wrqo/BH4IsGbNGl27du34CzlhM5seeohfvlhKMpXiu+85MV340L1l7nikgLLcZsh6Iy2bQHmPbt7Ma1evzll5r36Zi/IyXg/2INv/xOxN32V2QwO87QcAbNy4kQkdT9epux5hV0fvpMrIlcnWJZ8EpS5BqQd4VxffE3xO1M4nXr6DHYkQjTVlMH2h3xFNSk9VO8w+zu8wxmfp2RAuhQe+Dmd8DBqOmXSRpy6pZ8OzbbQcjjNrWtnkYzSmyATqbOThvkEbIumnUz8KoVJ47Oc5Ke7MZQ0AbHy2NSflGVNsvBwmeSPwEHC0iOwVkfd7ta+0rr5Bm6bAT5X1sOA0eGFDTopb1ljFnNpy/rTNErwxE+FZglfVd6vqbFWNqOpcVf2xV/sCSKnS3Z+wFrzfFq+FtmaITT4piwjrls/kge0H6U8kJx+bMUUmMF00fQnnvJ9NU+Czua91fu7/a06KW3d0A70DSR6xi56MGbfAJPjeQWdUhyV4n6VPrrY+k5PiTlsyg9JwCfc2WzeNMeMVnASfcBO8TVPgr4o6qGqE1uacFFceDXHmspn8/qn9pFI2fbAx4xGYBN/jTENjffD5YOZyaNuWs+IuOG42B7r62bzLJh8zZjwCk+Bf7qKxUTS+q50Ph/fmrLizj2mkLFLC7X/dl7MyjSkGwUnwbheNteDzQO186GmlJDmQk+IqS8OctbyBO59sIWndNMZkLTgJ3u2iqbY+eP9NmwtAaX97zoq88PgmDsb6uf+5tpyVaUzQBSbB97kt+KpS66Lx3ZEEn7tkfNbyRmZURbnxkd05K9OYoAtUgq+MhgiVjDC1r5k60+YBUBbP3dDGaLiEvztxLvdua6W1K56zco0JsgAleKiyaQryQ80cQCjtP5jTYi967TySKeXmLbk7gWtMkAUowav1v+eLcBSqZ1EWz21/+eKZVZy2pJ5fPLSLgUQqp2UbE0SBSfDxhPW/55WaOTntg0+7/G8W09IV57YnbMikMWMJTIJ3WvCW4PNGTVNOR9GknblsJstnVfOj+3fg5d3IjAmC4CT4pFoLPp/UzPEkwYsIl//NYp490M3dT7fkvHxjgiQwCd66aPJMzWzCyT6Id+W86AuPb2JpQxVfu/tZEknrizdmJIFJ8H0JtVE0+aRmjvOze3/Oiw6HSvjnNxzNC2093PKYjagxZiSBSPCplBJP2FWseaWmyfnZ9ZInxZ97bCOr5tfytbuf5XD6MmZjzCsEIsH3DCRQoNq6aPJH9WznZ5c3o11EhC+9ZSUdPQN87Q+5m7nSmCAJRIKP9ScAu9AprxxJ8LnvoklbOWcal562iF8+vJtNO3J/QteYQheMBB93E7y14PNHpIyByDTPumjSPvGGZSysr+TjN23lUE9uZq80JigCkeC7rQWfl/pL6zzrokmrLA3z7Xev4mCsn0/d8le765MxGQKR4NMteOuDzy/9pTM8T/DgdNV87o3HcM8zB7j2D896vj9jCkUgMmK6D95G0eSX/tIZ0PHQlOzr0tMW8nxrjO9ufIGm2nIuPmXBlOzXmHwWjAQfty6afNRXPgvindDb4dyM20Miwr9eeCwth+P8y2+eQgTec7IleVPcAtFF0xV3xkHbSdb80lfujqRpf2FK9hcOlfDd95zIWcsbuOrWp/juxu02X40paoFI8EeGSVqCzyt95e7VrB1Tk+AByiIhvn/xat58fBP/edez/NONj9M7kJiy/RuTT4KR4OMJSkPY3ZzyTF95I0gJHHxuSvcbDZfwrXedwKfOO5o7ntzP+d/8s42TN0UpGAm+P0F52JJ7vtGSCDSsgH1bp3zfIsI/rl3KjR88BVV41w838cmbn2BfZ9+Ux2KMXwKR4Lv7E9j51Tw150R4aQv41Bd+yuJ67vrY6/jg6xbx2637WHvtRq657Wl2HuzxJR5jplIgEnwsnqDCWvD5ac4aZyRNa7NvIVREw1x1wQru/cSZXHh8E9dv2sXaazfyDz95hFsf33vkJL0xQROIdq/TReN3FGZYy84FBJpvg8YVvoYyr66Ca99xPJ8692h++fBubnp0Dx+/6QmioRJOW1rPaUvqOWVxPStm1xAOBaLtY4qcp2lRRM4DvgmEgP9R1a96sZ9YPEGVteDzU/UsWPQ62HIdnPZPEK30OyIaasr4+DnLuPL1R/H4nkPc8dcWNj7XysZnnXvIlkdCLJtVzTGzqlk+q5pFM6vYF0vRN5CkPBryOXpjsudZgheREPDfwDnAXuBREblNVZ/J9b6644PMrLIEn7fWXQU/ORdufh+c9xWoWwzi//EqKRFWL6hj9YI6rmYFrd1xHt7RwWO7D/Fsi3NLwPWP7jmy/uceuIv6yiizppVRVxmltiJKXUWE6ZVRassjVETDlEdDlEdCVERDlEdDVETDlEVKCIdKiJQI4VAJ4ZAQLhHCJSVEQoLkwWdhgsnLFvxJwHZV3QEgIuuBtwC5T/D9Ccprc12qyZn5p8AF/wW//wx8+0QIl0FFPZSEQELuzxLA30TXALzZfQBoHSRrlcFkip6+OCXhCINJJdmdInlYSaaUpCqpMe4aOOA+RiROzSXj9SufvvpzGW7d4dd8tQWq7Lxvgp91Hv0xmp9KsfP+YHSl1VMJax/OebleJvg5wJ6M13uBk4euJCKXA5cDNDY2snHjxnHvaOV0ZXZpYkLb5qNYLBaIuryyHkspPel71LdvpizeQmSwC9HUkQfk6RWnIeeRKEkQDocRnC9N5hdHVUkqJBVSGY+k6pHnijOQaPSfL38GQz8NzXgy3CelQ16M9mmq6oT+a5BxHiOvj6iKIl5e+zKFv5J9Uu7Nd15VPXkA78Dpd0+/fi/w7dG2Wb16tU7Uhg0bJrxtvglKXYJSD1WrSz4KSj1UJ1cXYLOOkFO9/P9mLzAv4/VcwPu5Y40xxgDejoN/FDhKRBaJSBR4F3Cbh/szxhiTwbM+eFVNiMhHgbtxejJ/oqpPe7U/Y4wxr+TpOHhVvRO408t9GGOMGV4wxhgZY4x5FUvwxhgTUJbgjTEmoCzBG2NMQInm0T0rRaQN2DXBzWcAB3MYjp+CUpeg1AOsLvkoKPWAydVlgarOHG5BXiX4yRCRzaq6xu84ciEodQlKPcDqko+CUg/wri7WRWOMMQFlCd4YYwIqSAn+h34HkENBqUtQ6gFWl3wUlHqAR3UJTB+8McaYVwpSC94YY0wGS/DGGBNQBZ/gReQ8EXlWRLaLyGf8jmcyRGSniDwpIltFZLPf8YyHiPxERFpF5KmM9+pE5B4Red79Od3PGLM1Ql2uEZGX3GOzVUTe6GeM2RCReSKyQUSaReRpEbnSfb/gjssodSnE41ImIo+IyBNuXb7ovp/z41LQffDujb2fI+PG3sC71YMbe08FEdkJrFHVgrt4Q0T+BogBP1fVle57/wl0qOpX3T++01X1037GmY0R6nINEFPVa/2MbTxEZDYwW1UfE5FqYAvwt8ClFNhxGaUu76TwjosAlaoaE5EI8ABwJfA2cnxcCr0Ff+TG3qo6AKRv7G2mmKreD3QMefstwM/c5z/D+ULmvRHqUnBUdb+qPuY+7waace6VXHDHZZS6FBz3Tnsx92XEfSgeHJdCT/DD3di7IA+6S4E/iMgW92bkha5RVfeD8wUFGnyOZ7I+KiJ/dbtw8r5bI5OILARWAQ9T4MdlSF2gAI+LiIREZCvQCtyjqp4cl0JP8MPdUr1w+5zgdFU9ETgf+IjbVWDyw/eAJcAJwH7gv3yNZhxEpAq4BfiYqnb5Hc9kDFOXgjwuqppU1RNw7lV9kois9GI/hZ7gA3Vjb1Xd5/5sBW7F6YIqZAfcvtN0H2qrz/FMmKoecL+UKeBHFMixcft4bwF+qar/675dkMdluLoU6nFJU9VOYCNwHh4cl0JP8IG5sbeIVLonjxCRSuANwFOjb5X3bgP+wX3+D8BvfYxlUtJfPNdbKYBj457M+zHQrKpfz1hUcMdlpLoU6HGZKSK17vNy4GxgGx4cl4IeRQPgDov6Bi/f2PvL/kY0MSKyGKfVDs69cm8opLqIyI3AWpxpTw8AXwB+A/wKmA/sBt6hqnl/8nKEuqzF6QZQYCfwoXR/ab4SkTOAPwNPAin37c/h9F0X1HEZpS7vpvCOy3E4J1FDOI3sX6nqv4pIPTk+LgWf4I0xxgyv0LtojDHGjMASvDHGBJQleGOMCShL8MYYE1CW4I0xJqAswZuCJyK1IvKPGa+bROTXHu0rIiJbvCjbmFyzBG+CoBY4kuBVdZ+qvt2jfZ0B/MWjso3JKUvwJgi+Cixx5wP/mogsTM/lLiKXishvROR3IvKiiHxURP6viDwuIptEpM5db4mI3OVO9PZnEVk+wr7OA36f+YY7cdR1IvKUOPP5f3y0MkWkUURudecDf0JETvPskzFFLex3AMbkwGeAle7kTenZBjOtxJl9sAzYDnxaVVeJyP8DLsG5EvqHwIdV9XkRORn4LnDWMPtaB3xxyHsnAHMy5o6vdd8fqcxvAfep6lvdexpUTazaxozOErwpBhvcOcS7ReQw8Dv3/SeB49wZCk8DbnamPAGgdGghItKEc0OG3iGLdgCLReTbwB04Uz6PVuZZOH9YUNUkcHjyVTTm1SzBm2LQn/E8lfE6hfMdKAE60/8BjOJ84O6hb6rqIRE5HjgX+AjOXYY+lmWZxnjG+uBNEHQD1RPd2J1X/EUReQc4Mxe6CXuoV/W/u+vPAEpU9Rbg88CJY5R5L/B/3PdDIlIz0diNGY0leFPwVLUdeNA9yfm1CRbzHuD9IvIE8DRDbv3o9pUfparbhtl2DrDRvUPPdcBnxyjzSmCdiDyJc2/RYycYszGjstkkjcmCO13txar6Yb9jMSZbluCNMSagrIvGGGMCyhK8McYElCV4Y4wJKEvwxhgTUJbgjTEmoCzBG2NMQP1/KjBM43LFAicAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "u = lambda t: 1/(1 + (t/10)**100)\n", "\n", "first_order(5, 1, 30, u)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "colab_type": "code", "executionInfo": { "elapsed": 6401, "status": "ok", "timestamp": 1557657355289, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg", "userId": "09038942003589296665" }, "user_tz": 240 }, "id": "eozE3YD2ivdl", "outputId": "994a72b4-9383-46d2-ad16-bb293cacd678", "pycharm": {} }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnaUlEQVR4nO3de3xcdZ3/8dcntyZNQts0JfRGWy4CpQilFVhhf7YiCoIiKCqru6vA4u7KD9zf/lxEFkF3FXbV/eG6XhZ3FW9YQUQFUURtZAELUihKL9I7vadN0iZpmut8fn+ckzJNc5kkc+bMnLyfj8c8MnOun8+czCfffM+Z7zF3R0REkqco7gBERCQaKvAiIgmlAi8iklAq8CIiCaUCLyKSUCrwIiIJpQIvecXMzjez9WbWZmbvGOO27jCz74TPjw+3WZyVQMfIzP7ZzPaZ2e58i20gZlZvZtfFHYeMjAp8jMxsi5kdCj/cu83sXjOrijuumH0K+A93r3L3H2Vro+7+SrjN3mxtc7TMbDbw98B8dz9uLLGZ2QfM7MnsRxktM5tgZnea2SvhZ2C9mX3UzCxtmXoz6zCzVjNrMbOVZvYxM5uQtswdZtYdfob6HvtjSSoPqcDH723uXgWcBSwEbok3nNjNAVbHHUS2mFnJAJPnAI3u3pDB+mZmOfucZnt/g+QP8ABwIfBWoBr4c+B64Av9lrvB3auB6QR/FN8LPJr+hwD4fvgHsu8xOVvxFzoV+Dzh7ruBxwgKPQBmdp6ZPW1m+83sRTNbkjbvA2a2KWzdbDaz96VNf8rMvmhmB8xsnZldmLbeDDP7iZk1mdkGM/urtHl3mNn9ZvatcLurzWxx2vybzWxHOO+Pfds1s6KwZbXRzBrDbdQMlquZ/VW476Ywlhnh9I3ACcDDYUtswgDr9u2n1czWmNkVmby/ZjbXzLyv4IStw38K36tWM/uFmdVm+N5/0MzWhuttMrMPpc1bYmbbw/dqN/CNfnG8CXgcmBHmeO8gsX3azJ4C2oETBjreZnYa8FXgT4ZruZrZ683sd+HvxO/M7PVp8wba30Xh784BM/sPwPpt75rwPWg2s8fMbE7aPDezD5vZemD9ALFcCLwZeKe7v+TuPe6+Ang/8GEzO6n/Ou5+0N3rgbcDfwJcOliuksbd9YjpAWwB3hQ+nwX8AfhC+Hom0EjQwikCLgpfTwMqgRbglHDZ6cDp4fMPAD3A3wGlwHuAA0BNOP83wJeBcoI/JnuBC8N5dwAd4T6LgTuBFeG8U4BtwIzw9VzgxPD5R4AVYQ4TgP8EvjdIzm8E9gFnh8t+EXhioPdkkPWvAmaE78l7gIPA9EGWvQP4Tlq8DpSEr+uBjcBrgIrw9V3Dvffh/EuBEwmK3hsIiuLZ4bwl4fv/L2F+FQPEtQTYnvZ6oNheAU4HSoBJwxzvJ4f5PasBmglaySXA1eHrqYPsb1q4v3cR/A79XZjTdeHy7wA2AKeFy/8j8HTa/pzgj1jNIPnfBfxmkFi3Ah9Ki+u6AZZ5AviX/sdYj6MfasHH70dm1kpQPBuA28Pp7wcedfdH3T3l7o8DzxEUHYAUsMDMKtx9l7und2s0AHe7e7e7fx/4I3CpBX2/FwA3u3uHu68C/ovgg9/nyXCfvcC3gTPD6b0EBWu+mZW6+xZ33xjO+xBwq7tvd/dOgg/duwb59/x9wNfd/flw2VsIWqBzM3mz3P0Bd98ZviffJ2ghnpPJugP4hru/7O6HgPt59b+nId97d/+pu2/0wG+AXwB/mrbdFHC7u3eG2x6Ne919tbv3EBTXoY73cC4F1rv7tz1oLX8PWAe8bZD9XQKscfcfuHs3cDewO23ZDwF3uvvacPnPAGelt+LD+U2D5F8L7Bok1l3h/KHsJPjj0efd4X9afY/lw6w/bqjAx+8dHvQxLgFO5dVf7jnAVem/uATFebq7HyRovf41sMvMfmpmp6Ztc4eHzZvQVoJW7wygyd1b+82bmfY6/YPcDpSbWYm7byBoqd8BNJjZsr6ulTDWh9LiXEvwB6FugHxnhPsEwN3bCFrHMwdY9ihm9hdmtiptXwsYviAMpn+ufSe4B33vwxguMbMVYRfTfoLCnx7DXnfvGGVMfbb1PcngeB9mr16R02ZmbeHkI97zUP/jvi3t+Yx++/d+8+cAX0h7b5oI/psZbHv97SN8LwcwPZw/lJnhPvvc7+6T0x5Lh1l/3FCBzxNhS/Be4HPhpG3At/v94la6+13h8o+5+0UEH4h1wNfSNjfT7IiTUMcTtHp2AjVmVt1v3o4MY7zP3S8g+IA7QTdEX6yX9Iu13N0H2u7OcH0AzKwSmJpJDGEL8WvADQTdC5OBl+jXP5wFg7734XmBBwmOU10Yw6P9YsjGEK1HbGOI491/ub4rcqo8OHkP/d7zUP/jnr6dXcDsvhfh79LstPnbCLpR0t+fCnd/erD4+/klcG74H+VhZnZOuJ9fD7ZiuM4i4H+G2L6EVODzy93ARWZ2FvAd4G1m9hYzKzaz8vAE3iwzqzOzt4fFsRNoI2gx9zkWuNHMSs3sKoK+0kfdfRvwNHBnuL3XAtcC3x0uMDM7xczeGBa4DuBQ2j6/Cny67190M5tmZpcPsqn7gA+a2Vnhtj4DPOPuWzJ4fyoJCsfecD8fJGjBZ9ug7z1QRtBVtRfoMbNLCE4YRmaY470HmGVmZUNs4lHgNWb2Z2ZWYmbvAeYDjwyy/E+B083syrCb7UbguLT5XwVuMbPTw/gmhb9nGXH3XwK/Ah40s9PD9/g8gt/Dr7j7QCdmJ5rZG4AfA8+GOckwVODziLvvBb4F3BYW48uBjxMUk23ARwmOWRHBJWM7Cf5VfQPwt2mbegY4meBf3U8D73L3xnDe1QQn9XYCDxH0FT+eQXgTCE6O7SPo2jg2jA2CS9t+AvwiPJ+wAjh3kBx/BdxG0AreRXCy8r0Z7B93XwN8HvgtQWE7A3gqk3VHYqj3PuzeupGgz74Z+DOC3KM01PH+NcFlpbvNbMCujfDYXxZuoxH4B+Aydx9s+X0EJ7PvCpc/mbT32d0fIvjvbZmZtRD8F3XJCHN6J7Ac+DnBH6zvAP8N/O9+y/1H+Du1h6AB9CBwsbun0pZ5jx15HXybmR07wngSyY7sqpVCZ2YfILjy4IK4YxGReKkFLyKSUCrwIiIJpS4aEZGEUgteRCShBhsIKBa1tbU+d+7cUa178OBBKisrsxtQTJKSS1LyAOWSj5KSB4wtl5UrV+5z92kDzcurAj937lyee+65Ua1bX1/PkiVLshtQTJKSS1LyAOWSj5KSB4wtFzPr/y3lw9RFIyKSUCrwIiIJpQIvIpJQedUHP5Du7m62b99OR8fQg/NNmjSJtWvX5iiqI5WXlzNr1ixKS0tj2b+IyEDyvsBv376d6upq5s6dy5EDJB6ptbWV6urqQedHxd1pbGxk+/btzJs3L+f7FxEZTKQF3sy2AK0EI9/1uPviodc4WkdHx7DFPU5mxtSpU9m7d2/coYiIHCEXLfilg41al6l8Le598j0+ERmf8r6LRqSQpFLOoe5e2rt6OdTVS3t3Dx3dKXp6U3T3Oj2pFD29Tk/Kg2kppzcVzut1Uu7BnTLCn6lU8NOd8KeHzz1tGqTCIUdenR9MH8iWLV282HPUkOtH8THet2Sso6AMt/rWLV083/1ydAHk0O7tXURxSX+kY9GY2WaCMbMd+E93v2eAZa4Hrgeoq6tbtGzZsiPmT5o0iZNOOuom60fp7e2luLg4G2GPyoYNGzhw4EBWttXW1kZVVdXwC+a5JOTh7rR2w7bGgxyycho7nOYOp63bae1yDoY/27udjl7oTg2/TcmdQvnfurrU+fcLR/dZWbp06crBur+jbsGf7+47w8H3Hzezde7+RPoCYdG/B2Dx4sXe/9tca9euzejkaVwnWfuUl5ezcOHCrGwrKd/QK7Q8Uiln3e5WXtjWzLpdrazb3cK63a20dvQQlIpOACaUFDG1sozJE8s4blIZp1WWMamihMqyEirKiqkoLWZiWTEVZSVMLCumvLSIkqIiSoqN0uIiSoos7XXwvLgomFdkwa4Mo8iC7j8DLJxmRYSvX51eFHYRHl7GgmWKBuk6rP9NPUvesCSj92SsvY9Rdl8W2u/XUKLKJdIC7+47w58NZvYQcA7wxNBr5ZfbbruN2tpabrrpJgBuvfVW6urquPHGG2OOTLKhsa2Tx1bvYfkfG3h2cxMHDnUDUD2hhFOnV/OOs2Yyr7aS5h0becsFr2PWlAomVZQW9HmXIjOKigo3fslcZAU+vH9kkbu3hs/fDHxqLNv85MOrWbOzZcB5o+2imT/jGG5/2+mDzr/22mu58soruemmm0ilUixbtoxnn312xPuR/NHdm+JnL+1m2bOvsGJTIymH2TUVXHz6cZx7Qg2vm1vDrCkVRxTx+vqtLJg5KcaoRUYuyhZ8HfBQ+CEpAe5z959HuL9IzJ07l6lTp/LCCy+wZ88eFi5cyNSpU+MOS0bhUFcv9z69hW88tZmG1k6Or5nI3y45ibeeMZ3TplcXdKtcZCCRFXh33wScmc1tDtXSjrIP/rrrruPee+9l9+7dXHPNNZHsQ6KTSjkPrNzG53/xMg2tnfzpybXc9c4zWPKaY9VVIYmmyyQzcMUVV/CJT3yC7u5u7rvvvrjDkRF4pbGdmx/8Pb/d1MiiOVP40vvO5nVza+IOSyQnVOAzUFZWxtKlS5k8eXKsl2LKyPz8pd38/f2rMDPuvPIM3vu62eqGkXFFBT4DqVSKFStW8MADD8QdimTA3fnCr9Zz9y/Xc+bsyXz5fWczc3JF3GGJ5JyGCx7GmjVrOOmkk7jwwgs5+eST4w5HhuHufPLhNdz9y/W88+xZfP/681TcZdxSC34Y8+fPZ9OmTXGHIRlwd+74yWq++dutXHP+PG677DR1yci4pgIviXHPE5v45m+3ct0F87j1UhV3EXXRSCL8/KXd3PmzdVz62ul8/K0q7iKgAi8JsL25nX/4wYucOWsSn7/qTF3bLhJSgZeC1tOb4qZlq3CHL159NuWluoxVpI/64KWg/deTm1m5tZkvvPcsjp86Me5wRPKKWvBSsLY1tXP3L1/mzfPruPysmXGHI5J3VOCHsWXLFhYsWHD49ec+9znuuOOO+AKSwz71yBqKzLj97YOPUSQynhVWF83PPga7/zDgrIreHigeRTrHnQGX3DXGwCTXfrelicfX7OGjbzlFX2QSGYRa8FJw3J1/+dk6jq2ewDXnz4s7HJG8VVgt+CFa2ociGi64pKSEVOrVG212dHRkfR8yMvUv7+W5rc18+ooFVJTpqhmRwagFP4y6ujoaGhpobGyks7OTRx55JO6Qxr2v1m9kxqRy3r14dtyhiOS1wmrBx6C0tJRPfOITnHvuucybN49TTz017pDGtVXb9vPM5ib+8dLTKC1W+0RkKCrwGbjxxht1k+08cc8TG6kuL+G95xwfdygieU9NICkYDS0dPLZ6D1efczxVE9Q2ERmOCrwUjAdWbqc35Vyt1rtIRgqiwLt73CEMKd/jS4JUyvn+77Zx3gk1zKutjDsckYKQ9wW+vLycxsbGvC2i7k5jYyPl5eVxh5JoKzY18kpTu1rvIiOQ9x2Zs2bNYvv27ezdu3fI5To6OmIrsuXl5cyaNSuWfY8XP1q1g6oJJbzl9OPiDkWkYOR9gS8tLWXevOG/rVhfX8/ChQtzEJHkWldPip+/tJs3z6/TcMAiI5D3XTQiT27YS0tHD5edOT3uUEQKigq85L1Hfr+LSRWlXHDStLhDESkoKvCS17p7Uzy+Zg9vnl9HWYl+XUVGQp8YyWvPbWmmtaOHN82vizsUkYKjAi95bfkfGygtNs4/qTbuUEQKjgq85LXl6xo4d95UDU0gMgoq8JK3tjW1s76hjaWnHht3KCIFKfICb2bFZvaCmWkgdRmR+peDL7ctPUVXz4iMRi5a8DcBa3OwH0mY327cx8zJFRp7RmSUIi3wZjYLuBT4ryj3I8nj7qzY1MS5J9RgZnGHI1KQLMpBvMzsB8CdQDXwf939sgGWuR64HqCurm7RsmXLRrWvtrY2qqqqxhBt/khKLmPJY0drilufOsS1C8r401mlWY5s5JJyTCA5uSQlDxhbLkuXLl3p7osHnOnukTyAy4Avh8+XAI8Mt86iRYt8tJYvXz7qdfNNUnIZSx73PrXZ59z8iL/SeDB7AY1BUo6Je3JySUoe7mPLBXjOB6mpUXbRnA+83cy2AMuAN5rZdyLcnyTIik2NzJxcweyaiXGHIlKwIivw7n6Lu89y97nAe4Ffu/v7o9qfJIe788zmJs47YWrcoYgUNF0HL3lnS2M7TQe7WDx3StyhiBS0nHw90N3rgfpc7EsK36ptzQAsPH5yvIGIFDi14CXvrHplP5VlxZx8bHXcoYgUNBV4yTsvbNvPGbMmUVyk699FxkIFXvJKR3cva3e1sPB49b+LjJUKvOSV1Ttb6O51zpo9Oe5QRAqeCrzklVXb9gOwUAVeZMxU4CWvrN55gLpjJnDsMeVxhyJS8FTgJa+s2dnCadOPiTsMkURQgZe80dWTYuPeNuarwItkhQq85I31Da1097pa8CJZogIveWPNzhYA5s9QgRfJBhV4yRtrd7VSUVrM3Km6g5NINqjAS95Ys+sApxxXrW+wimSJCrzkBXdn7a5Wdc+IZJEKvOSFhtZODhzq5tTjNMCYSLaowEte2NDQBsBJ05Jxj02RfKACL3nhcIE/VgVeJFtU4CUvbGhoo7q8hGnVE+IORSQxVOAlL2xoaOPEaVWY6QoakWxRgZe8sGFvm7pnRLJMBV5id+BQN3tbO1XgRbJMBV5ipytoRKKhAi+x26graEQioQIvsdu4t42y4iJm10yMOxSRRFGBl9htaTzI7JoKjUEjkmUq8BK7rY3tzNEIkiJZpwIvsXJ3XmlqZ85Udc+IZJsKvMRqX1sX7V29zFH/u0jWqcBLrLY2HgRQF41IBFTgJVZbG9sBOF5dNCJZpwIvsdra1E6RwawpFXGHIpI4kRV4Mys3s2fN7EUzW21mn4xqX1K4Xmk8yPRJFUwoKY47FJHEKYlw253AG929zcxKgSfN7GfuviLCfUqB2aoraEQiE1kL3gNt4cvS8OFR7U8KU3ANvAq8SBTMPbqaa2bFwErgJOBL7n7zAMtcD1wPUFdXt2jZsmWj2ldbWxtVVckYyyQpuQyXx6Ee529+2c5Vrynl0hPKchjZyCXlmEBycklKHjC2XJYuXbrS3RcPONPdh30AE4HbgK+Fr08GLstk3XD5ycByYMFQyy1atMhHa/ny5aNeN98kJZfh8li3q8Xn3PyI/2TVjtwENAZJOSbuycklKXm4jy0X4DkfpKZm2kXzDYI+9T8JX28H/jnTvzDuvh+oBy7OdB1Jvh37g0skZ0zWFTQiUci0wJ/o7v8KdAO4+yFgyJGhzGyamU0On1cAbwLWjT5USZod+zsAXSIpEpVMr6LpCou0A5jZiQQt+qFMB74Z9sMXAfe7+yOjjlQSZ+f+Q5QWG9OqdKNtkShkWuBvB34OzDaz7wLnAx8YagV3/z2wcEzRSaLtaD7EcZPKKdIwwSKRyKjAu/vjZvY8cB5B18xN7r4v0sgk8XbuP8RM9b+LRCajPngzOx/ocPefElwR83EzmxNlYJJ8O/Yf0glWkQhlepL1K0C7mZ0JfBTYCnwrsqgk8bp7U+xp6WCWCrxIZDIt8D3h9ZaXA//u7l8AqqMLS5JuT0sHKdclkiJRyvQka6uZ3QK8H/hf4ZUxpdGFJUm3o/kQoAIvEqVMW/DvIbgs8lp33w3MBD4bWVSSeDsPBAV+pq6BF4lMplfR7Ab+Le31K6gPXsbgcAt+kgq8SFQyvYrmSjNbb2YHzKzFzFrNrCXq4CS5dh7ooKayjIoyjQMvEpVM++D/FXibu6+NMhgZPxpaOqg7pjzuMEQSLdM++D0q7pJNu1s6qDtGQxSIRCnTFvxzZvZ94EekjUHj7j+MIihJvj0tnZw+fVLcYYgkWqYF/higHXhz2jQHVOBlxHp6U+xr66RukrpoRKKU6VU0H4w6EBk/9rV14Y66aEQilulVNLPM7CEzazCzPWb2oJnNijo4SaY9LcE48HXVasGLRGkkd3T6CTCD4EtOD4fTREbscIHXVTQikcq0wE9z92+4e0/4uBeYFmFckmCvFnh10YhEKdMCv8/M3m9mxeHj/UBjlIFJcu1p6aS4yJiqOzmJRCrTAn8N8G5gd/h4VzhNZMT2tHQwrWoCxbqTk0ikMr2K5hXg7RHHIuPEntZOdc+I5ECmV9GcYGYPm9ne8EqaH5vZCVEHJ8nU0NLBsTrBKhK5TLto7gPuB6YTXEnzAPC9qIKSZNujYQpEciLTAm/u/u20q2i+Q/BNVpER6ezppbm9W9fAi+RApkMVLDezjwHLCAr7e4CfmlkNgLs3RRSfJExDSzCUka6BF4lepgX+PeHPD/Wbfg1BwVd/vGRkX1tQ4Gury2KORCT5Mr2KZl7Ugcj40NjWBcDUSvXBi0Qt06torjKz6vD5P5rZD81sYbShSRK92oJXgReJWqYnWW9z91YzuwB4C/BN4KvRhSVJ1Vfgp1aqi0YkapkW+N7w56XAV9z9x4A+oTJi+9q6qJ5QQnmp7sUqErVMC/wOM/tPguEKHjWzCSNYV+SwfW2d6p4RyZFMi/S7gceAi919P1ADfDSqoCS59rV1Ululf/5EciGjAu/u7UADcEE4qQdYH1VQklyNbV26gkYkRzK9iuZ24GbglnBSKfCdYdaZbWbLzWytma02s5vGFqokQdBFoxa8SC5k+kWnK4CFwPMA7r6z77LJIfQAf+/uz4fLrjSzx919zejDlULW3Zuiub1bLXiRHMm0D77L3Z1w/BkzqxxuBXff5e59fxBagbUEt/uTcar5YPAlJ51kFckNC+r2EAuYGXAbQXG+CLiTYIiC+9z9ixntxGwu8ASwwN1b+s27HrgeoK6ubtGyZctGmEKgra2NqqqqUa2bb5KSS/88trb0cvvTHdxw1gQWH5fpP4/5ISnHBJKTS1LygLHlsnTp0pXuvnjAme4+7IOga+Yi4LPA54CLMlkvXLcKWAlcOdyyixYt8tFavnz5qNfNN0nJpX8e9X9s8Dk3P+LPbm6MJ6AxSMoxcU9OLknJw31suQDP+SA1NdNm1G+B/e4+oksjzawUeBD4rrv/cCTrSvI09g1ToHuxiuREpgV+KfAhM9sKHOyb6O6vHWyFsGvnv4G17v5vY4pSEuHwODS6Dl4kJzIt8JeMYtvnA38O/MHMVoXTPu7uj45iW5IA+9q6KCspompCYfW/ixSqTIcL3jrSDbv7k4CNOCJJrH1tndRWlhH8cyciUdN4MpIzzQe7qFH3jEjOqMBLzjS1dzNlogq8SK6owEvO7G/vokbjwIvkjAq85EzTwS614EVySAVecqK7N0VrR48KvEgOqcBLTjS3B+PQ1FSWxhyJyPihAi850XywG4Ap6oMXyRkVeMmJwy14ddGI5IwKvORE31DBasGL5I4KvOREU9iC10lWkdxRgZec6GvBT56ok6wiuaICLznRdLCbyrJiykuL4w5FZNxQgZec2N/epf53kRxTgZecaNIwBSI5pwIvOdGsYQpEck4FXnKiqb2LKTrBKpJTKvCSE80Hu9UHL5JjKvASua6eFG2dPfoWq0iOqcBL5Pa361usInFQgZfINR0eSVIFXiSXVOAlck36FqtILFTgJXJ9QwWrBS+SWyrwEjkNFSwSDxV4iVzfQGOT1EUjklMq8BK55vZgoLEJJRpoTCSXVOAlcvvbu5is7hmRnFOBl8g1t3cxRTfbFsk5FXiJXHN7twYaE4mBCrxErlldNCKxUIGXyAVDBauLRiTXIivwZvZ1M2sws5ei2ofkv57eFC0dPWrBi8Qgyhb8vcDFEW5fCsCBQ8G3WNWCF8m9kqg27O5PmNncqLYvhaG5va/Ahy34A9uhaTMcaoJUD6RS4L3gHmOUw6vbvRZW7Yo7jKxISi5JyQNgWsNGYEnWtxtZgc+UmV0PXA9QV1dHfX39qLbT1tY26nXzTVJyaWtrY/lTzwBw6KWf0fLrr3BM64aYoxqd0wDWxR1FdiQll6TkAXBiySTq6y/I+nZjL/Dufg9wD8DixYt9yZIlo9pOfX09o1033yQll/r6euYeexp1z/yCq7b/EyVlFfCWz0DdApg4FYpKoKgYrCh45LEVzzzDeeeeG3cYWZGUXJKSB8ALz/4uks987AVekq25vYuPlDxIcXc7/NWvoPakuEMalY6KrVAzL+4wsiIpuSQlD4DO8q2RbDe/m01S8FpbD/D24qfpPuPdBVvcRQpVlJdJfg/4LXCKmW03s2uj2pfkrym7nqLSOil97bviDkVk3InyKpqro9q2FI7apufpooSyOefHHYrIuKMuGonUcW1r2Fh8ApToi04iuaYCL9HxFMd3vsyWCafGHYnIuKQCL5GZ0NlIuXfQOPGEuEMRGZdU4CUyFYd2AtBeNSfmSETGJxV4iUxFe1DguyarBS8SBxV4iUzZwZ10eCklk2fGHYrIuKQCL5Ep7djLDq9lSuWEuEMRGZdU4CUyZZ1N7PIajQUvEhMVeIlMRVcju5mq+7GKxEQFXqKR6qWyu4ndPoWpVSrwInFQgZdoHNxLESl2ew216oMXiYUKvESjZQcAe62WYyo0KrVIHFTgJRotwTXwHRV1mFnMwYiMTyrwEo2W4F6ZvVXTYw5EZPxSgZdotOwIhgmuro07EpFxSwVeotGykwavYWp1RdyRiIxbKvASCW/dyS6fQm21rqARiYsKvEQitX87u7yG2ioVeJG4qMBL9qVSFLXsZIfXUqsvOYnERgVesu9gA5bqYofXcmx1edzRiIxbKvCSffu3AbDDa5k5WSdZReKiAi/ZdyAo8Lt8KnWT1AcvEhcVeMm+A9sBaCubxoSS4piDERm/NEiIZN+BbbTbRMomVMYdici4pha8ZN/+beyiltoKjUEjEicVeMm61L6XebmnjhlV+vUSiZM+gZJd3Yew5s287LNV4EVipk+gZNe+lzFP8cfULGZW6tdLJE76BEp27XoRgM3F86irVB+8SJxU4CW7Nv2GJptC7Zz5lBSpwIvESQVesqeni96Ny3mi5zTOO1HjwIvELdICb2YXm9kfzWyDmX0syn1JHlj9EMWHGnk4dT7vWDgz7mhExr3ICryZFQNfAi4B5gNXm9n8qPYnMUr10rV+OZ2P/ANrUnOYfMZbNQaNSB6I8pus5wAb3H0TgJktAy4H1mR7Rxv+6Wzm9Haw9YkiwA9Pt7TnAObOQGyIdY7c3sDrZLqce/99DbzcfJy99YPFN3jsg8Uw2uWOzvHV9dPXK6ObMnrY7rV8qfZWPnP5GYOsJyK5ZD5I0Rvzhs3eBVzs7teFr/8cONfdb+i33PXA9QB1dXWLli1bNuJ9lTz9WSzVjVlf2Xm1/BydXfo8G2hyv/I11PaOXM77lexMYsCOXi6VSlFUVJRRrJns14Zb7qhzoSPLo9eKaZgwl9a6c5lfV0lRmFNbWxtVVVVHrVWIlEv+SUoeMLZcli5dutLdFw80L8oW/ECXUBxVI9z9HuAegMWLF/uSJUtGvqclS6ivr2dU6+ahpOSSlDxAueSjpOQB0eUS5UnW7cDstNezgJ0R7k9ERNJEWeB/B5xsZvPMrAx4L/CTCPcnIiJpIuuicfceM7sBeAwoBr7u7quj2p+IiBwp0vHg3f1R4NEo9yEiIgPTN1lFRBJKBV5EJKFU4EVEEkoFXkQkoSL7JutomNleYOsoV68F9mUxnDglJZek5AHKJR8lJQ8YWy5z3H3aQDPyqsCPhZk9N9jXdQtNUnJJSh6gXPJRUvKA6HJRF42ISEKpwIuIJFSSCvw9cQeQRUnJJSl5gHLJR0nJAyLKJTF98CIicqQkteBFRCSNCryISEIVfIFP0o29zWyLmf3BzFaZ2XNxxzMSZvZ1M2sws5fSptWY2eNmtj78OSXOGDM1SC53mNmO8NisMrO3xhljJsxstpktN7O1ZrbazG4KpxfccRkil0I8LuVm9qyZvRjm8slwetaPS0H3wYc39n4ZuIjgBiO/A65296zf9zUXzGwLsNjdC+7LG2b2v4A24FvuviCc9q9Ak7vfFf7xneLuN8cZZyYGyeUOoM3dPxdnbCNhZtOB6e7+vJlVAyuBdwAfoMCOyxC5vJvCOy4GVLp7m5mVAk8CNwFXkuXjUugt+MM39nb3LqDvxt6SY+7+BNDUb/LlwDfD598k+EDmvUFyKTjuvsvdnw+ftwJrgZkU4HEZIpeC44G28GVp+HAiOC6FXuBnAtvSXm+nQA96yIFfmNnK8Gbkha7O3XdB8AEFjo05nrG6wcx+H3bh5H23RjozmwssBJ6hwI9Lv1ygAI+LmRWb2SqgAXjc3SM5LoVe4DO6sXcBOd/dzwYuAT4cdhVIfvgKcCJwFrAL+Hys0YyAmVUBDwIfcfeWuOMZiwFyKcjj4u697n4Wwb2qzzGzBVHsp9ALfKJu7O3uO8OfDcBDBF1QhWxP2Hfa14faEHM8o+bue8IPZQr4GgVybMI+3geB77r7D8PJBXlcBsqlUI9LH3ffD9QDFxPBcSn0Ap+YG3ubWWV48ggzqwTeDLw09Fp57yfAX4bP/xL4cYyxjEnfBy90BQVwbMKTef8NrHX3f0ubVXDHZbBcCvS4TDOzyeHzCuBNwDoiOC4FfRUNQHhZ1N28emPvT8cb0eiY2QkErXYI7pV7XyHlYmbfA5YQDHu6B7gd+BFwP3A88Apwlbvn/cnLQXJZQtAN4MAW4EN9/aX5yswuAP4H+AOQCid/nKDvuqCOyxC5XE3hHZfXEpxELSZoZN/v7p8ys6lk+bgUfIEXEZGBFXoXjYiIDEIFXkQkoVTgRUQSSgVeRCShVOBFRBJKBV4KnplNNrO/TXs9w8x+ENG+Ss1sZRTbFsk2FXhJgsnA4QLv7jvd/V0R7esC4OmIti2SVSrwkgR3ASeG44F/1szm9o3lbmYfMLMfmdnDZrbZzG4ws/9jZi+Y2QozqwmXO9HMfh4O9PY/ZnbqIPu6GPhZ+oRw4Kh7zewlC8bz/7uhtmlmdWb2UDge+Itm9vrI3hkZ10riDkAkCz4GLAgHb+obbTDdAoLRB8uBDcDN7r7QzP4f8BcE34S+B/hrd19vZucCXwbeOMC+lgKf7DftLGBm2tjxk8Ppg23z34HfuPsV4T0NqkaXtsjQVOBlPFgejiHeamYHgIfD6X8AXhuOUPh64IFgyBMAJvTfiJnNILghQ3u/WZuAE8zsi8BPCYZ8HmqbbyT4w4K79wIHxp6iyNFU4GU86Ex7nkp7nSL4DBQB+/v+AxjCJcBj/Se6e7OZnQm8BfgwwV2GPpLhNkUioz54SYJWoHq0K4fjim82s6sgGLkwLNj9HdX/Hi5fCxS5+4PAbcDZw2zzV8DfhNOLzeyY0cYuMhQVeCl47t4IPBWe5PzsKDfzPuBaM3sRWE2/Wz+GfeUnu/u6AdadCdSHd+i5F7hlmG3eBCw1sz8Q3Fv09FHGLDIkjSYpkoFwuNr3u/tfxx2LSKZU4EVEEkpdNCIiCaUCLyKSUCrwIiIJpQIvIpJQKvAiIgmlAi8iklD/H9bqKjViii9uAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "u = lambda t: 1 - 1/(1 + (t/10)**100)\n", "\n", "first_order(5, 1, 30, u)" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "VxHD0mXNjPGq", "pycharm": {} }, "source": [ "## Analytical approximation to a square wave input\n", "\n", "An analytical approximation to a square wave with frequency $f$ is given by\n", "\n", "$$\\frac{4}{\\pi} \\sum_{k=1, 3, 5,\\ldots}^N \\frac{sin(k\\pi/N)}{k\\pi/N}\\frac{sin(2\\pi ft)}{k}$$\n", "\n", "where the first term is the *Lanczos* sigma factor designed to suppress the Gibb's phenomenon associated with Fourier series approximations." ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "colab_type": "code", "executionInfo": { "elapsed": 7318, "status": "ok", "timestamp": 1557657356220, "user": { "displayName": "Jeffrey Kantor", "photoUrl": "https://lh5.googleusercontent.com/-8zK5aAW5RMQ/AAAAAAAAAAI/AAAAAAAAKB0/kssUQyz8DTQ/s64/photo.jpg", "userId": "09038942003589296665" }, "user_tz": 240 }, "id": "AtvueB7KZndR", "outputId": "a913c9a4-8c19-4728-911a-b8ff43bc51ec", "pycharm": {} }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABPOUlEQVR4nO2deXgcxZn/P++MZjS6b8mWZEs+AdscxuYmwVwBEs4cJCEESDYhu0l+kD0SINkQdrO5s0nY3MluQi4ghHAFSDiCZU5jsLHxfcu2JOu+j5FGM/X7o3tkWZqRRtKMpnu6Ps8zjzTdM91vTXXVt+qtt6pEKYVGo9FonIcr2QZoNBqNJjloAdBoNBqHogVAo9FoHIoWAI1Go3EoWgA0Go3GoWgB0Gg0GoeiBUBjK0TkPBHZKyK9InLtDK91j4j83vx/vnlNd1wMnSEi8l8i0ioijVazLRIiUiMin0i2HZqpoQXAwohIrYgMmIW/UUTuE5HsZNuVZP4T+JFSKlsp9Vi8LqqUOmxeMxiva04XEZkH/CuwTCk1Zya2icgtIvJy/K1MLCKSLiLfEJHDZhnYKyKfFxEZ9ZkaEfGLSI+IdIvIRhG5U0TSR33mHhEJmGUo/OpMSqIsiBYA63OVUiobOA1YCdyVXHOSThWwPdlGxAsRSYtwuApoU0o1x/B9EZFZK8fxvl+U9AP8CbgYeDeQA3wUuBW4d8znPquUygHmYojmh4CnRwsF8EdTQMOv/HjZb3e0ANgEpVQj8AyGEAAgImeLyKsi0ikiW0Rkzahzt4jIAbN1dFBEPjLq+Csi8kMR6RKRXSJy8ajvlYvIEyLSLiL7ROSTo87dIyIPichvzetuF5HVo87fISL15rnd4euKiMtsme0XkTbzGoXR0ioinzTv3W7aUm4e3w8sBP5ituTSI3w3fJ8eEdkhItfF8vuKSLWIqHCFZLYuv2r+Vj0i8qyIFMf4239MRHaa3zsgIp8adW6NiNSZv1Uj8OsxdlwCPAeUm2m8L4ptXxORV4B+YGGk/BaRk4CfAedM1vIVkXNF5A3zmXhDRM4ddS7S/S41n50uEfkRIGOu93HzN+gQkWdEpGrUOSUinxGRvcDeCLZcDLwLeJ9SaptSalgptR64EfiMiCwe+x2lVJ9Sqga4GjgHeE+0tGpGoZTSL4u+gFrgEvP/SmArcK/5vgJow2ghuYBLzfclQBbQDZxgfnYusNz8/xZgGPhnwAN8EOgCCs3z64CfAD4MsWkBLjbP3QP4zXu6gW8A681zJwBHgHLzfTWwyPz/c8B6Mw3pwM+BB6Kk+SKgFTjd/OwPgRcj/SZRvv8BoNz8TT4I9AFzo3z2HuD3o+xVQJr5vgbYDywFMsz335zstzfPvwdYhFEpXoBRaZ5unltj/v7fMtOXEcGuNUDdqPeRbDsMLAfSgLxJ8vvlSZ6zQqADo5WdBnzYfF8U5X4l5v3ej/EM/bOZpk+Yn78W2AecZH7+34FXR91PYYhcYZT0fxNYF8XWQ8CnRtn1iQifeRH41tg81q/xL90DsD6PiUgPRuXaDHzFPH4j8LRS6mmlVEgp9RzwJkalBBACVohIhlLqqFJqtNukGfiBUiqglPojsBt4jxi+5/OBO5RSfqXUZuB/MSqGMC+b9wwCvwNONY8HMSq0ZSLiUUrVKqX2m+c+BXxJKVWnlBrEKJTvj9L9/wjwK6XUJvOzd2G0YKtj+bGUUn9SSjWYv8kfMVqYZ8by3Qj8Wim1Ryk1ADzEsd7XhL+9UuoppdR+ZbAOeBZ4x6jrhoCvKKUGzWtPh/uUUtuVUsMYle9E+T0Z7wH2KqV+p4zW9gPALuCqKPe7AtihlHpYKRUAfgA0jvrsp4BvKKV2mp//OnDa6F6Aeb49SvqLgaNRbD1qnp+IBgxxCXO92VMLv9ZO8n3HoAXA+lyrDB/nGuBEjj38VcAHRj/YGJX3XKVUH0br9x+BoyLylIicOOqa9cpsHpkcwmg1lwPtSqmeMecqRr0fXdD7AZ+IpCml9mG09O8BmkXkwbDrxrT10VF27sQQjLII6S037wmAUqoXo3VdEeGz4xCRm0Rk86h7rWDyCiMaY9MaHoCP+tubNlwhIutNF1YnhjCMtqFFKeWfpk1hjoT/iSG/R5BjEUW9ItJrHj7uNzcZm+9HRv1fPub+asz5KuDeUb9NO0ZvKNr1xtKK+VtGYK55fiIqzHuGeUgplT/qdeEk33cMWgBsgtmSvA/4rnnoCPC7MQ92llLqm+bnn1FKXYpRYHYBvxx1uQqR4wbJ5mO0mhqAQhHJGXOuPkYb71dKnY9RASgMN0fY1ivG2OpTSkW6boP5fQBEJAsoisUGs4X5S+CzGO6LfGAbY/zTcSDqb2+OS/wZI5/KTBueHmNDPJbgPe4aE+T32M+FI4qylRFcAGN+c5Ox+T76OkeBeeE35rM0b9T5IxhumtG/T4ZS6tVo9o/heeAss0c6goicad7nhWhfNL+zCnhpgutrTLQA2IsfAJeKyGnA74GrROQyEXGLiM8cYKwUkTIRudqsPAeBXowWd5hS4DYR8YjIBzB8tU8rpY4ArwLfMK93CvAPwB8mM0xEThCRi8wK0A8MjLrnz4CvhV0AIlIiItdEudT9wMdE5DTzWl8HXldK1cbw+2RhVCwt5n0+htEDiDdRf3vAi+EKawGGReQKjAHNhDFJfjcBlSLineASTwNLReQGEUkTkQ8Cy4Ano3z+KWC5iLzXdOPdBswZdf5nwF0isty0L898zmJCKfU88HfgzyKy3PyNz8Z4Dn+qlIo0cJwpIhcAjwMbzDRpJkELgI1QSrUAvwW+bFbW1wBfxKhsjgCfx8hTF0ZIXANGV/gC4NOjLvU6sASjK/014P1KqTbz3IcxBh0bgEcxfNXPxWBeOsbgXSuG66TUtA2M0L0ngGfN8Yz1wFlR0vh34MsYreijGIOpH4rh/iildgD/DbyGUfGdDLwSy3enwkS/vek+uw1jzKADuAEj7Ylkovx+ASNstlFEIrpOzLy/0rxGG/AF4EqlVLTPt2IMtn/T/PwSRv3OSqlHMXp/D4pIN0Yv7Ioppul9wFrgbxiC9nvg/4D/N+ZzPzKfqSaMBtKfgcuVUqFRn/mgHD8PoFdESqdoT0oix7uCNamOiNyCETlxfrJt0Wg0yUX3ADQajcahaAHQaDQah6JdQBqNRuNQdA9Ao9FoHEq0hZgsSXFxsaqurp7Wd/v6+sjKyoqvQUlCp8V6pEo6QKfFqswkLRs3bmxVSpWMPW4rAaiurubNN9+c1ndrampYs2ZNfA1KEjot1iNV0gE6LVZlJmkRkbEzvQHtAtJoNBrHogVAo9FoHIoWAI1Go3EothoD0Gg0mmQQCASoq6vD75/pIq7TJy8vj507d074GZ/PR2VlJR6PJ6ZragHQaDSaSairqyMnJ4fq6mqOX0h39ujp6SEnJyfqeaUUbW1t1NXVsWDBgpiuqV1AGo1GMwl+v5+ioqKkVf6xICIUFRVNqZeiBUCj0WhiwMqVf5ip2qhdQBbGHwhytMtPS88gOxq6WFCSzQVLx83l0Gg0mmmhBcACDAwFebuukz1NPexr7mWv+WrpGRz32T//07lJsFAzmsYuP2/UtvNGbTsNnQP8y6UnJNskxzM4HOTtui7eqG3nrcOdnFqZx2cuXJxssyyPFoAkMDQcYsPBdtbububN2na2N3QzHDIW5ctOT2NxaTZrlpZQXZzFnFwfRdle5hVmct2PX+GBDYe5UncCZpX2viGe39nE+v1tvHGonSPtxj7mmV43SsGR9s3ceZpeVHE2CYYUGw91ULO7mTdq29lS18XQsLEHTGGWl+d2NHFyZX5yjbQBWgBmiVBIsf5gGw+/WcezO5roHRzGm+bitHn53PrOhayqKmBZeS5zcn1R/XjnLynmtf1tXFmih24STe/gME+93cAjm+p5o7adkIKiLC9nVBdyy7kLOKO6gGVzc/njm0f40qPbaOjNSLbJKY9Sii11XfzpzSP8bVsjbX1DpLmEFRV53HxOFaurC1ldVUCOz8OZX3+ev2xpSKnG0le/+lUqKiq4/fbbAfjSl75EWVkZt91227SvqQUgwQwNh3hkUx0/W7ef2rZ+ctLTeM/Jc7l0WRnnLi4i0xt7Fiybm8vTWxsZGM5MoMXOpqnbz8/W7eePbxyhfyjIopIsPnvhYt61fA7Ly3PHifP5i4sB2NsZjHQ5TRwIhhRPbz3KT2r2s/NoNxkeN5csK+Oy5WVcsLSEHN/4mPczqwvZeKgjIQLwH3/Zzo6G7rhec1l5Ll+5avmEn7npppu46aabuP322wmFQjz44INs2LBhRvfVApAglFI8vrmB7zyzm/rOAU6pzOMHHzyNy5bPIcPrntY1F5caMcBH+0KTfFIzVboGAvzg+T384fXDBEOKa04r5yNnVXH6/PwJIyvmFWTi87g42qvzJBE8v6OJb/x1J/tb+lhcms3XrzuZq06dG7HSH80Jc3J4fmcTgVDqNJaqqqooKirirbfeoqmpiZUrV1JUVDSja2oBSAAHWnr54qNbWX+gnVMq8/jadSu4YGnJjMPIFpdmA+jKJs78ZUsD//nkDtp6B7l+9Tw+c+Fi5hXGVnG4XMKikmwa+voSbKWzaOzyc/fj23h2RxOLS7P58Q2nc8WKObhcsZWhxaXZhBQ09cV/bGaylnoi+cQnPsF9991HY2MjH//4x2d8PS0AcebhjXXc/fg20lzC165bwYfPmB/zQzsZVUWZuAQa+/WAYzzoGxzm7se38+dNdZxSmcevbj6DkyvzpnydxaXZvLyrJwEWOpMXdjXxrw9tYSAQ5M4rTuQfzl+Axz21ca9wY6khxXrL1113HXfffTeBQID7779/xtfTAhAngiHFf/5lO7957RBnLSjk3g+tZE6eL6738LhdFGen0zWo/c0zpaFzgI/9+g32NPdw+8VLuO3iJbinKdRVhZk84VcMB0OkTbGi0hxDKcWPXtjHfz+3h2Vzc/nRDStZWJI9rWtVFRkbp7QOpJYAeL1eLrzwQvLz83G7p+dKHo0WgDgwOBzktgfe4pntTXzyHQu484qTpl2ZTEZZro9Of29Cru0U9jb1cNOvNtDrH+Z3Hz+L85cUz+h6pbk+FNDWN0RZbnxF3ykEQ4q7H9/GH14/zHUrK/jGe0/G55l+BZednkaW102XP7V6y6FQiPXr1/OnP/0pLtfTzZUZEgiG+Oz9RuV/95XL+NJ7liWs8gcozUmnYzC1HurZ5GBrHx/+5XqCIcUfP3XOjCt/YKTSb+pO3kqRdkYpxb8/ZlT+n7pgId+7/tQZVf5hynJ9KVVWdu3axeLFi7n44otZsmRJXK6pewAzQCnF5/+0hed2NPEfVy/n5nOrE37P0lwfGw6kVrd2tjjaNcCN//s6IQUP3nr2iJ94ppTlpgPQ1D1+5rZmcr71t908sOEwn16ziC9cfmLcrluam057x0DcrpdsTjzxRA4cOBDXa+oewAz46br9PLa5gX9719JZqfzBqGx6hoyehyZ2/IEgn/rdRroGAvz242fGrfIH3QOYCQ9vNObIfOSs+Xz+svguqVGW66MzhXoAiUALwDR5cU8L331mN1edWj6ra46U5hiVTWuvbm1OhXue2M7bdV187/pTWVEx9UifiSjK8iJoAZgq2+q7+NKjWzlnYRH/cfXyuK+2WZqTTqdfoZQWgWhoAZgGrb2D/PMfN7OkNIdvve/kWV0mtiDTmADT0ReYtXvancc31/PgG0f4zIWLeNfyOXG/fprbRbbXGATWxEb/0DCfuX8TBZlefnjDyoRETxVnpzMUgv4hHTUXDS0AU0QpxZcf20aPf5gf3rBySks5xIM8UwA6B3RlEwvNPX6+8sR2Vs7PT+iqndkeobNf50msfOeZ3Rxq6+cHHzqN4uz0hNyjINMLQIfOl6hoAZgif9vWyF+3NfK5S5ewtCz69myJIvxQd/brHkAs3PPEdvqHgnzn/acmNDoryyM6T2Jk46EO7nu1lpvOqeLshTNbymAiRhpLOl+iknQBEBG3iLwlIk8m25bJ8AeCfP2vOzlxTg63vmNhUmzQAhA7r+1v4+mtjfy/CxfHddA3EtkeoUPnyaSEzAmTZTk+7ohjxE8kdFmZnKQLAHA7MPFW9xbhvldrOdI+wL+/Z1nSZnzmh8cAdLd2QkIhxX89tYPyPB+ffGfixTrLI3TpPJmUJ7Y0sKWui89fdgJZ6Yl1n+anmLv00KFDrFixYuT9d7/7Xe65554ZXTOp8wBEpBJ4D/A14F+SactktPUO8uMX9nHxiaVxmTw0XXweN16XsXqlJjqPvFXP9oZu7v3QaXGZVDQZ2V7oaNF5MhH+QJBv/20XKypyuW5lRcLvd6yxFOd8+eud0Lg1vtecczJc8c34XjMGkj0R7AfAF4CoznQRuRW4FaCsrIyamppp3ai3t3fa3wV4eM8QvYPDXFzcM6PrxIPMNMXO/YepqWlKqh3xYKb5EolgSPGtlwaoznWR07GHmpq9cb1+JLwqwEBAePbva/G6rb95+EQkIk8Anj8UoKFriBuXwosvrov79ccS3mXvre27mec/OKNr5eXl0dNjLPiXHhjCFRyesX2jCQWGGOyZeEHBUChEKBQasWNwcJDBwcGR92H8fn/M+Zc0ARCRK4FmpdRGEVkT7XNKqV8AvwBYvXq1WrMm6kcnpKamhul+t6s/wGfXvsC7T5nLDVeePq1rxJOcV/6KL6+INWtWJ9uUGTOTfInGY2/V0zKwma+9/3QuSkDYZyReOPwcMMSpZ5xj+/WAEpEnQ8MhvvjaWs6oLuDT75u9fa19f3+KgrJK1qxZNqPr7Ny5k5wcs5169ffiYNl4vJOd9xqfCNuhlCI9Pf2YXSY+n4+VK1fGdM9kjgGcB1wtIrXAg8BFIvL7JNoTld+8Vkvv4DCfWWONTaazPOiQwyiEQoqf1OxjaVk2l5xUNmv3zfYarX49NhOZxzfX09Dl59OzvFF7lkdSJk9KS0tpbm6mra2NwcFBnnxy5nEzSRMApdRdSqlKpVQ18CHgBaXUjcmyJxpDwyF++1otF55QwrLy3GSbA4RjzrW/ORLr9rawp6mXf7xgUdz2YYiFbI9xL50v41FK8b8vHeSkubmsWTq7m/Rme4WuFMkTj8fD3XffzVlnncWVV17JiSfOPIoq2WMAlueZ7Y209g7N2lo/sZDlEWo7U+Ohjjd/WH+Y4mwvV55SPqv3zTJ3KNQ9s/G8UdvB7qaeWZ81D0a+pEoPAOC2226b0SbwY7FCGChKqRql1JXJtiMSv1t/iPmFmbxzyey2XCYi2yN0DQzpNU7GUN85wAu7mrh+9Ty8abP7aOseQHR+v/4QOb40rjp1dkUZzN6yjpiLiiUEwKrsbephw8F2PnJW/LZ1jAdZXggElV7jZAwPbjiMAj585vxZv3dYAPRksONp6x3kr9uO8r7TK2d92RTQM7QnQwvABDzyVj1ul/C+VZXJNuU4stKMykbPBThGKKR4ZFM971xSEvOG7vHE6wav26XzZAxPbGkgEFTccNbsizKYE/QGAnHpLduhxz1VG7UARCEUUjz+Vj3vXFKcsMWqpkumRwvAWN481EF958CsTDCKhIiQm5Gm82QMj29uYNnc3KSsmwWQ6THmhfTNsLfs8/loa2uztAgopWhra8Pniz0MWQ8CR+GN2nYauvzccUVi1yuZDllaAMbx+OZ6fB4Xly6bvdDPseRmeOjWeTJCbWsfm490clcSy9Do3nL2DJaeqKyspK6ujpaWlniZNmX8fv+klbvP56OyMnaPhRaAKDy2uZ5MrzupFUo0Ms1c0wJgMDQc4qmtR7l02ZyEry8zEXkZHp0no3hiSwMiJGXwN8xIb7k/QEV+xrSv4/F4WLBgQbzMmhY1NTUxT/CKFe0CikAwpHhmexOXnFSWlIGrydA9gON5dX8rnf0BrkliRQNaAMby5NsNnFFdSPkMKt6ZosvKxGgBiMCmwx209w3xruXWa/3DsVaNdjcYPL+ziQyPO6mL9IEWgNEcbutnT1Mvl83SUhzRCM/P0PkSGS0AEXh+RxMet/DOWZ61GCsZaSCiBQCMga+/72zmnUuLZ2XVz4nQAnCM53caCxVeclJpUu3ITNONpYnQAhCB53Y2cfbCInJ9nmSbEhGXCLk+XdkAbG/o5miXn4tncd2faORleOj2BwiFrBspMls8v7OJJaXZVBVlJdUO7QKaGC0AY9jf0suBlj5LDv6ORrc2DZ7f2YQIXHRicluaYOSJUtAzGN+lgu1G10CADQfbucQCZciXBi7RAhANLQBjWLfbCPO68ITkVygToWPODdbuamblvHxLzNXIzTB6jE53N7y8t5XhkEq6+wfM3rJuLEVFC8AYXtnXSnVRZlJmk04F3QMwQvu21nfxDous05RnCoDT8+Xlfa3kpKdxamV+sk0BdFmZCC0AowgEQ6w/0MZ5i5MbTRIL+qGG9QfbCCk4d1FRsk0BjgmA03sAr+5v5ayFhUnbN3ss4bEZzXiskUMW4e26TvqGgpxvGwFwtq/51X2tZHjcrJxfkGxTAN0DAKjr6OdQWz/nLrJOGdKNpehoARjFy3vbEIFzLNKinIjwsgNWXpsk0byyv40zFhTO+tLP0dACAK/ubwOwVC9ajwFExxolxyK8sq+VFeV55GdOtjtn8snL8DAUDOEPhJJtSlJo6vazr7mX8ywk1loAjF5ZcXY6S8uyk23KCHl6jaaoaAEwGRwOsrmuk7MWFCbblJgY8Tc71Lf5+sF2wFq9tUyvmzSXOFoA1h9o5+yFhbO+89dEhF1ATu4tR0MLgMn2hm6GhkOsrraGP3kynN7a3HSogwyPm2VzrbFPM4SXhHauu6Ghc4DGbj+rq6xVhnJ9HgJBxUBAb6A0Fi0AJpsOdQBwukUGFCfD8QJwuINT5+VZJtIkjJMHHDeaZWhVlbV60U4vKxNhrdKTRDYd7qCyIIPS3Ng3U0gm4WUquhy43d3AUJAdDd2sslhLE5w94LjpcAc+j4sT5yZn85doaAGIjhYAjAXFNh7qsGSFEg0nP9Rv13UyHFKW7K05ecBx0+FOTqnMx2PBXhk4s7E0GdbKqSRR3zlAU/egFgCbsPGw4WqwSvz/aJzqAvIHgmyv77JkGXJyWZkMLQAYLRewj/8fjq0748SHetOhThYWZ1GYZb1w3TyHrtG0tb7L0r0ycGZZmQwtAMC2+i68aS5OmGMt3+VEuF1CTrozK5u36zo5bV5+ss2IiLHswLDjQg63HOkEsGS+aAGIjhYADAE4aU6O5XyXk+HETcibu/009wyyvCIv2aZEJC/DQzCk6HXYktDbG7opy02nJCf5q7KOJceXpjdQioK9arwEoJRiW32XZSuUiXDiIlfbG7oBWFFunfj/0Ti1tbmtvosV5dYsQy4H95Ynw/ECcKR9gG7/sGUf3olw4oDj9oYuAJZpAbAMA0NB9rf0styieQKQl+m8shILjheAbWaFsqLCug9vNJy4Kcy2+m6qizLJseh2nU4cnN/Z2E1IYeletBMbS7GgBaC+izSXsLTMPgPAYZz4UG8/am13nRP3BBhxy1k8X5xWVmJBC0BDN0vLcvB53Mk2Zco47aHu6g9wpH3A2q6GEQFwziDw9vou8jM9lOdZdxZ9ODpLczyOF4AdDV2WrlAmIi/Dgz8QYnDYGYtcbT9quussPF7jxDGA7Q3drCjPs9QKoGNxWmMpVpImACIyT0TWishOEdkuIrfPtg3tfUO09g7ZKv5/NE6rbHY39gBYbq2Z0WSnp+ES5+RJMKTY09TDiRYvQ05eo2kiktkDGAb+VSl1EnA28BkRWTabBuxpMiqUxaXW2bxiKuQ6zN+8p6mX/EwPJdnWizUP47QloY+09zM4HLL8GFquz8PQcAi/XhL6OJImAEqpo0qpTeb/PcBOoGI2bdhrCoDVH95oOK0HsLeph6WlOZZ2NYCz3A17m3sBWGKhHcAi4bSyEitpyTYAQESqgZXA6xHO3QrcClBWVkZNTc207tHb2zvuuzU7BslIg91vrWePxSuV0YTTsr/TaM28vGETPQctkZVTJlK+REIpxc6Gfs6YkzbtZyCRjE6HK+DnYH2TJe2MhVjzBOCZ/UMANO7ZTM0B65WhcFrqjxoDwM+ve5WKHHsOfU4lX2JGKZXUF5ANbATeO9lnV61apabL2rVrxx374M9fVdf++OVpXzNZhNOyv7lHVd3xpHp0U11yDZoBkfIlEk1dA6rqjifVr18+kFiDpsnodNz4v+vVNT+y33MVJtY8UUqp2x/YpM75+vOJM2aGhNOybnezqrrjSbXhYFtyDZoBU8mXsQBvqgh1alKlUEQ8wJ+BPyilHpnt++9t6mWJTf3/4KxJR8dcDdZ31zlpjaa9zb22yBO9J0BkkhkFJMD/ATuVUt+b7fu39Q7S1jdkW/8/QL75UHf0DyXZksQTHrC3uq8ZjHzpdIAABEOKfc32aETlZxplxQn5MhWS2QM4D/gocJGIbDZf756tm+9psk+LMhppbhe5vjQ6+lJfAPY295KXYe0IoDCFWV46+4cIhVJ7Sei6DntEAAEUmHtHOKGsTIWkjRwqpV4GkjZqtK/FFAAbtF4moig7nXYHdGv3NvWwtCzb8hFAAAWZXkLKcM0VWHDTmnhxrBFl/TKUk56Gxy20O6C3PBXsORweBw629JHhcTPHJpvAR6Mg00N732CyzUg4+1v6WFRi/YoGoCjbqPTbUry1ud9sRC2yQSNKRCjI9NLem9p5MlUcKwC1bX1UFWXiclm/RTkRhVnptPeldg+gayBAe98Q1cVZyTYlJgoyTXdDirc2a1v7KMrykmvRlVnHUpjl1T2AMThWAA629rGwxB4VykQUZqV+D+BQWx8A1UX2yK/wXsVtKd7arG3rs40ogykAKd4rmyqOFIDhYIgj7f22qVAmojArnY6+QErvQXuw1RCABTapbMICkPo9AHuVocIsrx4EHoMjBaCuY4DhkLJV6yUahVkehoIh+oZSd42T2tZ+AOYXZibZktgIC0AqtzYHhoI0dvupLrJHnoCRL6k+LjNVHCkA4RblwhQQgLC/OZUHtw619TE3z0eG1x57Nvg8bjK97pQWgNqwW85GZagwy0vXQIDhYCjZplgGRwuAnR7eaIQjTlJ5cOtgW5+tXA2Q+v7m8LiMXdxyMNo1l9pBE1PBkQJQ29ZHji+NohSI0R7pAaTwQHBtq70GGyH1BeCg6ZarspkLCFJ/bGYqOFIADrb2saA4yxaTiiajKMuYGZuqoaBd/QE6+gO28jVD6gtAbWsfxdlecmwSAgpQmOmM6Kyp4GgBSAUKssz1gFK0srGjrxmMyiaVBcCWbrls3QMYi+MEYDgY4miX3zYRJZORnZ6G1+1K2eiGWpvNAQiT6j2AQ219VNktTzKdMUN7KjhOAI52+QmGFBX5Gck2JS6ICAVZnpTtAdR1DAAwr9Be+VWQ5WUgEGQgBcNz/YEgTd2DtmtE6QXhxuM4AajvNCqUygJ7PbwTUZCZuvHNdR0DFGV5yfTaa8ezcIBBKkZnHe3yA1BZYC9R9pir56Zyz2yq2KtUzZTBHua+8DnuTFNUB0th52sw0AluL6R5jb/udHC5QYWASLNrpzlwPN0B54wC2PkX6DgEV3wbsorGfaQo22tvv+bWh2H7o7DqFvDlQctuIw/cXubUH+KaLIH95hwAFSmGewaD+dPJl7QMGOiAt34Hp1wPy68b95HRrU1b9jZ7W+CvX4CSE+Ckq6B1Dwz2Qlo6/uYBLnEdYln/EOzfa5GyIpCRD1v+CIM98O5vg3e8iyrVXXNTxVkC8PrPqap7gn9MAx78S7KtmTq55fCur447XJDpZUdDdxIMmjnu4T54/DMw7IddT447f3v4n9/Nqlmxc/AlWHzpuMPhHoBte2YvfRe2m5v01XzjuFMnAf/rBf4+61bFTulJcO5nxx3WAnA8MQmAiGQC/wrMV0p9UkSWACcopcaXWCuz91kOZ5zEvwT+iYevSoeixZBdBsEh4zU8aPwNDYO4x7dEpr3ezjS/p0LQ3WA8zM/+O+x7PqIAFNl4intBx1aj8v/oo0bLzeUx0guo4BDX3LuW604u4mNnV5r5Icfny2znCRh2DvvBlQYPfAiOrGdsUbK9v3nPM7D0CrjoS9C0A0qWQmYRBAPc9+JuHn3zAH++dTVphKxRVkJBo6xUnA4P3QR7n40qAA2d/mnalnrE2gP4NcbG7eeY7+uAPwH2EQAVhMatbPe+i2DhYjj1vGRbNDUqVsOL3zFEagyFWel0DQQIBEN43PYa1snu3W9UIPPPAc/xrpLWnkHeHp7H++Yth6rq5Bg4EYPGNpUc3QKsOu5UuAfQ2mu/CXppgV7oOAinfxTmnGy8RrFlqIfWHB9pVWcnycJJqFwNOx6PeKowy8vW+q5ZNsi6xFpbLFJKfRsIACilBkjibl7TIX2wHQL9bA/MtadPtnAhoKDzyLhTJTnGZDA7TnDJ7K+HgqpxlT8YWw4C1s2v9BzIKjHGZ8aQl+HB4xZabZgnGQNHjX+KT4h4vr5jgAorDwAXLjTGaPzjK/ri7HTaelN/u85YiVUAhkQkA7N/JiKLAFs1bdIHWwHY0Z9r7Yc3GgXVxt+O2nGnis0JLi09tsoSwMyXvMqI50YitqwcAlpQHTFPRITi7HRb9gDCZYX8eRHP13X0WzsCaKSsjBfmkpx0hkOKLr05PBC7AHwF+BswT0T+gDH884WEWZUAfP4WAI4Ei+wZAjryUB8cdyrcA7BjZePzt0Je5Iqm3pwDYNkeAEQVADDyxbaiDJA7XpgDwRCN3X4qrZ4nEDFfwmWlxYZlJRHEJABKqeeA9wK3AA8Aq5VSNYkzK/6EH+qjqtDaD280sssgzQed41s1xdnmQ223yiYYwDvUDrkVEU/XdQyQl+Gx9noz+VXQVYeEhsedKs62pwD4/C1GqGtm4bhzjV1+Qsri82jyq4y/EXvLNi0rCSImARCR8wC/UuopIB/4oohUJdKweJM+2EogLYdeMu3pAnK5DBHoaRp3yratmp6jCAryIgtAfeeAtVv/ALlzQQXxBMaH4ZbY2QWUVxExHj88M9vSZSgjHzyZ0Bu9rNgxXxJBrC6gnwL9InIq8HngEPDbhFmVANIH2+j2lgIWdylMRFYx9LeOO+zzuMnxpdmvVdNtDjZG7QFY3NcMkFkMEFEAinOM8NygzQYc0wfbjDknEQgPzNsiX/rGlxXdAzieWAVgWBmbzl4D/I9S6l4gJ3FmxR9PoJsuyaEwy0tWuk3nv2WVQF9LxFMl2en26wH0txl/I7gaABo6/ZRbXayzSgDwDnWOO1WSnU4wpGw3S9sT6DZi/iMQXgZiTp5vNk2aOlnFEctKri8Nb5pLC4BJrALQIyJ3ATcCT4mIG7CwY3Y8nkAPraEc+7b+wXyox7dqAIpz0mm120M90G78jVDZ9A4O0zs4bI+Khsg9gJIcw3a7uRs8gZ4JBaAoy0t6msW358wqidhbFhF7NpYSRKwC8EGMsM9/UEo1AhXAdxJmVQLwBLppGs60ftd1IsLd2gizLG35UPebApARebARYK7lBcDoAXgCkWLObRieGwqSNtwbMU8Amrr9lOVaPE9g0saSrfIkgcQaBdSolPqeUuol8/1hpZR9xgCUIm24l/pBn70FIKsEQgHShvvGnbJlyGF/GyFJMyZUjaGp2xAAy1c2vnwQN96h8QIwMjhvp3zxdxkD81F6AI1dfuv3yuCYAERrLNkpTxJIrFFA7xWRvSLSJSLdItIjIvZZfWywG5cK0hLMtnb42mRM0trs8Q/jD9ho/fmBdgKenIjRJiO+ZqsLgMsFmUUR88SWESeTjMs0ddtFAEogOIg7ODDuVElOui1naCeCWF1A3wauVkrlKaVylVI5SqncRBoWV8yHulPl2LsHYLbKIvub7VjZtBPwRH6Mwj0AW1Q2mUUR8yQ7PY10uw04ht1yEQRgcDhIW9+Q9UUZJi4r2V7a+wZtF52VCGIVgCal1M6EWpJI+jsAaCfH2vHLk2EWyokFwEYtm/52htMiB5M1dvnJz/Tg81h8sBGiCoCI2M81NxB9XKa520iH7QUgJ52QgrY+G+VLgog1HvJNEfkj8Bij1gBSSj2SCKPijtkD6CLbdnvLHseIAPSMO2XL+OaBdgKegoinGrv99qhoADIL8LSMn6ENNnQ3TOACagyPy9ihV5YRvayMNJZ6hijNsUFaEkisPYBcoB94F3CV+bpypjcXkctFZLeI7BORO2d6vaiYrZrsglJ7tCijMdKqif5Q20oA+tuiuoAau2wSbQJRewBgw+Ug+qOH5jbaZVwGJuwtjzSW7OQuTRAx9QCUUh+L943NuQQ/Bi7F2F/gDRF5Qim1I973Cj/UpXMizzi1Dd5scHkiPtRFWTYTAKVgoINAURQXULefZXNtMsyUUWiETio1bkC7JCedjYc6kmTYNDAjs1ze7HGn7DUuM3kPwDZlJYHEGgVUKSKPikiziDSJyJ9FJPIavrFzJrBPKXVAKTUEPIgx0zjudLY1ElTCknmRp7fbBhHILCRtePxD7U1zUZTlpanHJrsdDXZDaDhiDyAQDNHaO2gPVwNAZhEuFTTSNIbSnHTa+4YYHLZJdNYEkVmNXX4yPG5yfTaYSZ+eB+KKWFbCbp+woFmegU7Y/Vc8EUKNZ8pUdgS7H/iA+f5G89j4zVBjpwIYvbtJHXDW2A+JyK3ArQBlZWXU1NRM+UaBt3ewimyKBo5QU1M/PWstwmrlwzXQEfF3yHINs/1APTU1bbNv2BTxDRzlbKA36BmXlraBEEpBd+MhamoakmLfVChrbOYkYP3ap/FnzD3uXFejse78E8+uoyTT+ru1LT+8m3R3Fq9FeL7e3ucn1xNi3bp1s2/YNDg3LQf62yKXFQ9s2rmfGqmbfcOmSG7XTk5/607SlnyBmpq8uF47VgEoUUr9etT7+0TkczO8d6QdxcbFZSmlfgH8AmD16tVqzZo1U75Rz7IFvPD8M7z/ioum/F3LcXAew53tRPodltS+QUOXnzVr3jH7dk2Vuo3wOrizS8alZeOhDlj3KheccSprTixNjn1TYc8g7LqXs09eCpXHbw0pe1r49bYNVC87jTOqI8fWW4oD36IzkB/x+frRzldZMEdYs+ac8d+zItvKyMQfMS3zNr+IKzuTNWtWz75dU2VXP7wFaTmlEdMyE2JtkrSKyI0i4jZfNwIzbWbWAaN3AqkEEtLcyymtIq98aSIuPftkFkT0a4IRnWGbbq0ZbRLJBWSbWcBhwiGT/eOLRHjANDyxzfKEXUARsFVkFkBGYdTB+bJc38igtuUxg1iiBUzMhFgF4OPA9UCj+Xq/eWwmvAEsEZEFIuIFPgQ8McNrpj4TRJzMzfXR3jdkj9nAIw919GUgbDHYCMdCJsMx9KMIp6HJLpVNlMgspRTN3TYalwGzrERuLM3N842EtVqe/uhlZabEGgV0GLg6njdWSg2LyGeBZwA38Cul1PZ43iMlmSDiJFw4m7sHmV9k8SUvJnioW3oGSXMJ+Rk2WXA2M3oPINeXRqbXbY8ewASRWd0DwwwFQ5SYIZS2ILNgwh5Aa+8gQ8MhvGkWH5vpbwOXh6A7/pNYY40CWigifxGRFjMS6HERWTjTmyulnlZKLVVKLVJKfW2m13MEE0SchFfOPNo1fv0Ty9HfBuJiOG38xLzW3kGKsr24XJGGiSxIeh4K17EY+lGICHNybeKamyAyKxwzHw6htAXhHkCEBeHm5vlQCprtEDU30G40MiJEZs2UWKXvfuAhYC5QDvwJY29gzWwzQWsz7J+1Rdd2oB0yCkDGP4ItPYMjk3Vsgctl9GQiuIDAcAPZRpSJ3isD7NUDyCjEpQIQ6B93asQ1Z4ey0t8edXXWmRKrAIhS6ndKqWHz9XsiROxoZoGRAcfxk4vCD7UtBrf626KuOd/aO2SvliZmpRlBlMHIl6ZuG0w6Mp+pSOszhRcZLLZTvkzUWMqz0eB8f3vUsjJTYhWAtSJyp4hUi0iViHwBY2ewQhGxQWxbChFuCUR4qHN8HrK8bnv0ACZo1bT22qwHgBmhEcEFBIy4gCy/+uQE0SYjAmCnfBkpK+PzZW6u4U+3RWMp7AJKALHOA/ig+fdTY45/HKMnMOPxAE2MTBBxAkbLxh4PdQfkzx93WCllSwEYTsuJKgBz83wMhxRtvYOUWjmMchIXkK0G5mHC8NzcjDR8Hpc9ykp/W3IFQCm1ICF310ydkW7tRP5mmzzU5aeNO9w1ECAQVCPbKdqFgCcH+g5HPDcnz2xtdvstLgAT9wBsNTAPx3oAA+PdpSLC3LwMjlq9t2xGZiXVBSQiHxCRHPP/fxeRR0RkZUIs0kzMSMRJFH9zbob1B7aUiurXbLVjtAlhF1BbxIgT20wGG2g3I7PGhxC39g7Zrlc20RgAQFluuvV7AGZkVrIHgb+slOoRkfOBy4DfAD9LiEWaiXG5CHiyo7qA5ub5aO6x+G5HQ70QHIz4ULf0GGvn2yraBLMHEByCofH7NdtmcL6v1RDlVIjMAmO/ZpjANZdhjzyBpAtAeGrpe4CfKqUeB+zVR08hJhpwLM/PIBhS1u4FdJsrfuSOX567xY7RJoxym0QQ5qIsL163iwarh4J210Nu5BVzW3sHbdcrw51GIG3ixlJTt5/hYGiWDZsC3ebilVHyZabEKgD1IvJzjOUgnhaR9Cl8VxNnhrz50HM04rnwnsd1HRaubLrMRWDzxgtAa48No00w8wSge3y+uFxCRUGGtfMEoKse8sav8m7XgXkw86U78hJjlQWZDIcUTVbeF6DLFIAI+RIPYq3Er8dYsuFypVQnUAh8PiEWaSbF75sDHZG3IDwmAOMnv1iGCR7q1l4bRptg5glAZ/R8sbwAdNdF7JXZdWAezHyZIE8A6totXFa6zeWqk9kDUEr1A83A+eahYWBvQizSTMpARhn0NsLQ+Ae3PN8GPYDuekAgZ+64U7aMNgH8PnPZ6o7aiOcrCzKot7IoD/aAvyuqKIP9BubBLCsdhyIOztujt1wHmcXgif86QBB7FNBXgDuAu8xDHuD3CbFIMynHWpvjww59HjelOekW7wHUQc4ccI9v5dtysBEIub2GoEUVgExae4cYGLLoSq1dZkszggDYdWAezLIy2B0xFNQWjaWuuoiu0ngRqwvoOozVQPsAlFINQPzXJtXExEBGmfFP+/6I5y3vb27eAcWR92ewZbhhmIJqaIucJ+HWZn2nRYW52dyKu3jJuFN2HZiHUWUlQr6EG0uWzROA5l1Ry0o8iHUm8JBSSomIAhCR8Us4amaNvqxqcHng8GuAGD2B0pNgeBAG2rnWdYhDrX7Y1jBmBcFR/yfieCgIKMM2l9t4PzwA6TkQChm+2MozoGk7nP1PEdPW2jvICXNs2rYoXwlv/gpa98Gevxl+26wS6GvmlLYOrnYdYHBzM1Tkj/pSgvNEKVBBEDe40oxzgQEj1NOXazw73mw49Bq406HkJNj96nHJarXjQnAmvdmLjH8OvQJdh43ouZITjXBdfyc3ZOzDXZ8G2/Ynoaxg9ILFDaEADPvBlwcBv+EmnXOKMQYw97QppXkqTCoAIiLAk2YUUL6IfBJjCYhfJswqzYSE3OmwcA28+kPjNYabw/88PJtWTZET3jPukFKKtt4himw42AjA0sth/U/gR6vGnVoA/I8XeHXcKetw4pWQNv63b+sbxO0S8mw2MA8w6Cs2KtLnvxLx/OfC/1i6rFyRsEtPKgBmy/9ajDGAbuAE4G6l1HMJs0ozOe/+Drz8fag6Dxa8A1p2GS25jEKeeLuee5/bzYOfOIOSHHPpgeMGwUb9H8/jLrfxf3AYVAhcLkjzGQOMYEwy2vIAZJfC/LPGJalvKMhQMERRlk0FYME74V1fM0J0V38c/J1G2rNKCbm8XPGDdbzv9Lnc+g5zZZXZyBNxGS1PFTLyBcDjM1qg/i6jl9K4FQ6ug3f8W8RktfcFKMi038D8CNf+BF7/OSy9zOiltew2WtoZBfzyxf08/MYhnr7tPNzh9M1WWVHKyAcVMt6npRt54vKANwve+j2ULYeiRdNP+yTE6gJ6DehUSunQT6tQuACu/p9j70eFieVV5rFf9VLrrqKk1GKLtV4SuSUG0N5rDDYWZNpUAETg3M9GPOUChgrq2DKYa7jrrETxEljx3qin2/sGKcyyX+t/hDknwzU/OvZ+1EB3VrmH3aEATb4FI4PCluHyryf8FrEOAl8IvCYi+0Xk7fArkYZppo8t5gJEoL3fEIBCu/YAJsEWcwEi0GH2AFIRW4SCJpBYewCJc0Jp4k5FOLyt3V4PdUef2QNIYQF4bkdTss2YMu39QywpzU62GQlhdGPpzAUW6y3PArEuBx15Kp3GkoTD2w5ZeYZjBNpNAbDtGMAkhOcC9A0Ok5Uea9sr+XT0DaVsr6yiIAMROGyzshIv9Ho+KUp1cRa1reNXprQyHf2p3QOoLjKip2vb7JMvoZCioz91BSA9zU15Xobtykq80AKQoiwoyrJVRQPQ1jdEmkvIsVHreCpUFxvr7Ne22qe12TUQIKRsPDAfAwuKszjYZp88iSdaAFKUBSVZtPYO0e0PJNuUmOnoG6Igy4uITcMNJyHcAzjY2ptkS2In1QfmwRDmgy29qAjrBaU6WgBSlBF3g426tu19QxSmcEszKz2Nstx0DtqoB5DqA/MAC4qz6fYP09Fvn8ZSvNACkKIsKA63Nu0jAKnsaw5TbTPXXKoPzAMsMF1zdior8UILQIpSVWQ/f3N7CkebhFlYYq/B+VQfmAd79pbjhRaAFMXncVOe57Nda7PAzjNOY6C6KIu2viG6BuzhbmgzewCp7JqbV5iJ2yW6B6BJLaqLs2zzUAdDis6BQEpXNGDkCdintdnRN4TP4yLD6062KQnD43YxryCDgzZqLMULLQApjJ0EoGsggFKpHW0CsLDYXnMB2vtSX5TBnvNm4oEWgBRmYXEWXQOBkYE8K9PugGgTMNwNInCgxR6VTUf/EIV2XZ57ClQXGQLgtFBQLQApzGJz/Za9TT1JtmRyOhwQbw7G2Mz8wkz2NdtjLkB731BKTwILs7g0m76hIA1d/mSbMqskRQBE5DsisstcVfRREclPhh2pztIyY2etPTaobNrsvhT0FFhSmsMeG4gyOCMyC0aVFZvkS7xIVg/gOWCFUuoUYA/HNpvXxJG5eT5y0tN0D8BiLC3L5mBrH0PDoWSbMikdDukBLC2zT285niRFAJRSzyqlzO2JWA9UTvR5zfQQERaXZduiVRMeA3CGAOQwHFKWHwgOBEP0DA47QgDyM72U5qSzu9H6veV4YoVVtz4O/DHaSRG5FbgVoKysjJqammndpLe3d9rftRpTSUtOaJC36oYtm/ZwWrbtGSLNBetfeSnZJk2LqeRJV7exIfijL7zOWXOtUASPJ5yW7kFjQLSlvpaamvokWzU9ppIvJd4Am/Y1UFPTkVijpklC6jClVEJewPPAtgiva0Z95kvAo4DEcs1Vq1ap6bJ27dppf9dqTCUtv3xxv6q640nV2uNPnEEzIJyWOx7eolb/13PJNWYGTCVPBoaG1YI7n1T//cyuxBk0A8Jp2dfco6rueFI9uqkuuQbNgKnkyz1PbFMn/vtfVTAYSpxBM2AmdRjwpopQpyas+aGUumSi8yJyM3AlcLFpoCYBHBvc6uWc7PQkWxOdroEAeRmpPQs4jM/jprooiz1N1nY3hGcrOyVflpblMBAIUtcxwHxzKZVUJ1lRQJcDdwBXK6Xss1iNDQkLwN5ma48DOEkAAJaUZbPHBnkCkOuQfHFiJFCyooB+BOQAz4nIZhH5WZLsSHnKctPJ9aWxq9HaD3W331kCsLQsh0Nt/fgDwWSbEpVuh/UAlpiRQLsdJABJGYFSSi1Oxn2diIiwrDyX7Q3dyTZlQroGAiwpzUm2GbPGSXNzCYYUe5p6OKUyP9nmRMRpApDr81BZkMGOo9YuK/FEzwR2ACvK89h5tJtA0Lpx5139AXJ91ouISRQnV+QBsLW+K8mWROeYC8hZ+bLNwnkSb7QAOIAVFXkMDYfY32LNQcdQSNEzOOyYliZAZUEGeRkettVbt7XZNRDA53GRnpa6K4GOZUVFHofa+m2zXPdM0QLgAFZU5AJYtrLp8Q+jlHMGG8Fwza2oyGV7g3Vbm04bmAdYXm6UlR0Wd5nGCy0ADmBBcTaZXrdlu7ZOCzcMs6I8j11Heyy7JIQTBWCF6ZqzalmJN1oAHIDbJSyba93WZrffoQJQkcdQMGTZEN3uAWe55QCKs9OZm+djm0XLSrzRAuAQVlTksb2hm2DIenPuHNsDMFub2y3qmnNiDwBgeblzBoK1ADiEFRV59A8FLTkQPCIAmc6qbKoKM8lJT2NLXWeyTYlI10CAXJ+z8gSMSKADrX30+FN/IFgLgEM4fX4+AJsOWW+hq5FwQ4dVNi6XcNr8fDYd7ky2KRHpHgg4amA+zMr5+SgFm490JtuUhKMFwCEsKM6iINPDRgsLgBPdDauqCtjd2G251mbQgaG5YVbOz0cES5aVeKMFwCGICKuqCth42HoPdddAgDSXkOl1Trx5mFVVBYQs2Np02izg0eT4PJxQlqMFQJNanF5VwIGWPsttEh8ebBSRZJsy65w2z5qtTadGZoVZVVXA5sOdlgyaiCdaABzEqvkFALxlsV5At0OjTcC6rU0nu+XAEICeweGUXxlUC4CDOKUynzSX8KYFKxsnDjaGWV1dwFsWa206NTIrzOqqQsB6PbN4owXAQWR43ZxcmcfrB9qSbcpxODXaJMwZ1YX0Dg5baqKeUyOzwswrzKA0J531Fisr8UYLgMM4b1ExW+q6Rny8VsCpE47CnLOoCIBX9lmnsnG6C0hEOG9xMa/tbyNkoZ5ZvNEC4DDOW1xMMKR4/UB7sk0ZwRAA5yw5PJbSHB8nlOXw6v7WZJsygtMFAIyy0tY3lNIbxGgBcBinV+Xj87h4ZZ81KhulFN1+Z8abj+a8xcVsONhumR3CugYCeN0ufB7nVhHnLQ73zKxRVhKBc3PXoaSnuTlzQREvW+Sh9geNSUdOF4DzlxQxOByyzEzt7oFhch0amhtmbl4Gi0qyLFNWEoEWAAdy3qIi9jX3crRrINmm0B8w/KtOF4AzFxSR5hJeskhl0+1wt1yY8xcX8/qBdgaHrdEzizdaABzImhNKAXhhV3OSLYE+UwCcGm0SJjs9jVVVBbywM/l5Ajo0N8w7l5YwEAiy3kJjZvFEC4ADWVqWzbzCDJ7f0ZRsU+gzg5Gc3gMAuHRZGbubejjc1p9sUxwfmRXmvMXFZHjcligriUALgAMRES49aQ6v7G+jb3A4qbb0D5s9AF3ZcOmyMgCe25n8ykYLgIHP4+adS4t5fmcTSqVeOKgWAIdy6bIyhoZDvLS3Jal29OkxgBGqirJYWpbNczsak22KFoBRXHJSGUe7/JbdU3smaAFwKGdUF5CX4eHZ7cltbfaHXUAOXXJgLJcuK+ON2g46krhgX0gpevxaAMJcfFIZLoFnLSDM8UYLgENJc7u4bHkZz+5oSmrsed+wwiWQ7dURJwBXrJhLMKR4etvRpNngH4aQ0r2yMIVZXs5eWMRftjSknBtIC4CDufa0CnoHh3k+iT7n/oAiN8ODy+XcePPRLC/PZXFpNo+9VZ80G0Yis7QAjHDtaRXUtvVbbt+GmaIFwMGctbCIObm+pFc2Tg8BHY2IcN3KCt6o7aCuIznRQCMD8zpfRrj85Dl401xJLSuJQAuAg3G7hKtPK6dmd0vSNonpC0CB9v8fx9WnlgPw+OaGpNw/HJqr8+UYuT4Pl5xUypNvHyUQDCXbnLihBcDhvPf0CoZDikc21SXl/n0BRV6mNyn3tirzCjM5c0Ehf3zjSFJWouw1XUD5Ol+O470rK2nrG0qpOQFaABzOiXNyOaO6gN+vP5SUyqYvoMjXvuZxfPTsKg6397MuCWG6fUNhAdD5MpoLTyylIj+D3752KNmmxA0tABpuPLuK2rb+pCx61RtQuqKJwGXL51Ccnc7vk1DZ6LkZkXG7hBvOms9rB9rY15waS0QnVQBE5N9ERIlIcTLtcDpXrJhLcbaX+16tndX7hkKK/gC6BxABb5qLD585jxd2N3OorW9W790bAJ/Hhc/jntX72oEPnjEPr9vFb15NjV5A0gRAROYBlwKHk2WDxsCb5uKmc6p5YVczO4/O3mzHHv8wCh1uGI0bz67C43Lxs3UHZvW+/cN6ee5oFGenc+3Kch568wjNPf5kmzNjktkD+D7wBSC1ZlbYlJvPqSY7PY0fr903a/cM7zqlBxsjU5br4/ozKnl445FZXbrbGJfReRKNf1qzmEAwxP+9dDDZpsyYpEy/FJGrgXql1JbJNpwQkVuBWwHKysqoqamZ1j17e3un/V2rkai0XFAuPPX2Uc7NfYHy7MS3DQ52GTOQj+zfRU3P7AlPIkhUnpzqDfFASPHlP7zIjcvS4379SHQNDON296VEeUlUvpw5x819rxxgRVojOd7ZmcSYkLQopRLyAp4HtkV4XQO8DuSZn6sFimO55qpVq9R0Wbt27bS/azUSlZbWHr9acfff1Md+vSEh1x/Lut3NquqOJ9UbB9tm5X6JJJHP151/flstuusptb+5J2H3GM15X31KffI3b8zKvRJNovJlb1O3WnjXU+rLj21NyPUjMZO0AG+qCHVqwpp5SqlLlFIrxr6AA8ACYIuI1AKVwCYRmZMoWzSxUZSdzmcvWswLu5pnZZXQzhEXkPY3T8S/XLoUn8fN15/eNSv36wvoPJmMxaU5fOSs+fzh9cPssfGm8bM+BqCU2qqUKlVKVSulqoE64HSlVOottWdDbjmvmvmFmfzXkzsTPuOxq9+YfZyn/c0TUpKTzmcuXMzzO5t4eW/iQ3X7AkqPy8TA5y5ZSqbXzVef3GHbReL0PADNcaSnufnSe05id1MPP1+3P6H36jTXgtYRJ5PzsfOqqSrK5IuPbk3oJj7+QJChkM6TWCjM8vLPlyzlpb2tPGrTNYKSLgBmT8AaO2FrAGMS0pWnzOXev+9NaFhoW98QPrcRhqqZGJ/HzXfefypHOvr51t8S5woKrwlVoHsAMXHzudWsqirgnie209hlv7BQXfI0EfnPa1aQl+Hls/dvoscfSMg9WnoHyU/Xy0DHypkLCvnYuQv47WuH+Nu2xHhMW3sHAcPtpJkct0v47gdOZSgY4rYH37LdQnFaADQRKczy8j8fPo3atn7+7U9bEuLjbO0ZJFcLwJS444oTOHVePv/60Gb2NffG/fotPVoApsqC4iy+8d6T2XCwnW/M0kB9vNACoInKuYuKueuKE3lmexPffXZ33K/f0jtInhaAKZGe5uZnN55OhtfNJ3/75kiFHS/C1yvO1i6gqXDdykpuObeaX71ykPtft8/iBloANBPyD+cv4ENnzOPHa/fzszgPCrf2DJI3S5NoUom5eRn8/KOraOzy89H/e53O/vjt5RB2ARVn6x7AVPnSe05izQklfOmxrTy+2R6DwloANBMiInztupO56tRyvvnXXfzw73vj4g7yB4J0+4d1D2CarKoq5Jc3reZASx83/PJ1mrvjMwDZ0jNIZhp6Ibhp4HG7+NmNqzhrQSH/8tAWHnrzSFyuOzgc5HBbP8EELNeuBUAzKW6X8L3rT+W6lRX893N7+OKjWxkantlgV9jVoMcAps/5S4r55c2rqW3r47qfvMruxplPSGrW4zIzwudx8783n8G5i4r4wsNv8/3n9sx4n41dR3t453fW8nZrME5WHkMLgCYmPG4X37v+VD69ZhEPbDjC+376KrWt01+m+Ei7sd9taYZ+BGfCBUtLeOhT5zAUDHHNj1/m/tcPz6iHdri9X+fJDMlOT+NXt5zB+06v5N6/7+XmX2+Y0cqhteZy4KWZ8c8XndOamBERvnD5ifzsxlUcbu/nintf4qc1+6fVGzg48lDr1uZMWVGRx1O3nc/qqkK++OhWbvrVBg5OQ5yVUhxq69d5Egc8bhff/cApfO26FWw42M67vv8iD244PK3eQG1rPyJQkhH/fNECoJkyl6+Yw19vfwfvWFLMt/62i8vvfZHHN9dPyUd5qK0fb5qLAp+ubOJBaY6P3378TO65ahmbD3dy2fdf5KtP7qBpCmMDbX1D9A4OU5aAlqYTERE+clYVT/6/81lamsOdj2zlup+8wtpdzVPqpR1q62Nurg+vWwuAxiKU52fwi5tW8383r8Ytwu0PbuaS763jvlcO0tE3eVTKwdY+qgozcU2yHLgmdlwu4ZbzFvD3f72Aq08r575Xa3nHt9dy1yNb2XKkc9JKJ+zSK83SeRJPlpTl8MdPnc33rj+V1t4hPnbfG1z1o5d56M0jMS3rcbCtj6qirITYlpT9ADSpw8UnlXHhCaU8u6ORn9Ts556/7ODrT+/inUuLuXRZGRedWBZxUtH2+i5WVhUAs7cDmVMozfXx3Q+cym0XLeGn6/bxyKY6HthwmMWl2Vy2vIxLl83hlIo8XK7jK/pt9V0AVM7CXhBOQ0R47+mVXHVqOY++Vc/P1+3nCw+/zT1PbOfCE0t517Iy1iwtJW/MKqyBYIhdR3v44BnzgPgvNaEFQDNjXC7h8hVzuXzFXHY0dPPwxjqe2d7I8zubga0sKM5i5fx8VlUVsKqqAJcIDV1+/nFBIQxqAUgU84sy+cZ7T+Gud5/Ek1uO8vjmen5as58fr91Pji+N0+blc/p8I09Om5/Puj0tVORnUKQHgROGx+3i+tXz+MCqSjYd7uDhjXU8t6OJp94+iggsLc3h9Kp8Vpr50tTtZyAQ5MwFhdAW/yXatQBo4sqy8lzuLl/Gl688iR1Hu3lxTyubDnewbncLj2w6NjnG53Fx+Yo57NhYmzxjHUKuz8MNZ83nhrPm09E3xLo9LWyobWfToQ7+54W9jPYM3XbRYuBo0mx1CiLCqqpCVlUV8l/XKjYf6eTlvUZZefLtozyw4dgcgoJMD2tOKGHDq/Gfja8FQJMQRITl5XksL88DjAiTw+39bDrcwbb6bi5YWkJpjo8dSbbTaRRkebl2ZQXXrqwAoMcfYMuRLjYd7iAQDPHpCxez/hUtALOJ2yUjvWOAUEixv6WXTYc72Hm0hytPmUumNzFVtRYAzawgIlQVZVFVlMV1K5NtjSZMjs/D+UuKOX9JcbJN0Zi4XMKSshyWlOUk/l4Jv4NGo9FoLIkWAI1Go3EoWgA0Go3GoWgB0Gg0GoeiBUCj0WgcihYAjUajcShaADQajcahaAHQaDQahyLx2N5vthCRFuDQNL9eDLTG0ZxkotNiPVIlHaDTYlVmkpYqpVTJ2IO2EoCZICJvKqVWJ9uOeKDTYj1SJR2g02JVEpEW7QLSaDQah6IFQKPRaByKkwTgF8k2II7otFiPVEkH6LRYlbinxTFjABqNRqM5Hif1ADQajUYzCi0AGo1G41AcIQAicrmI7BaRfSJyZ7LtmQkiUisiW0Vks4i8mWx7YkVEfiUizSKybdSxQhF5TkT2mn8LkmljrERJyz0iUm/my2YReXcybYwFEZknImtFZKeIbBeR283jtsuXCdJix3zxicgGEdlipuU/zONxz5eUHwMQETewB7gUqAPeAD6slLLlboQiUgusVkrZanKLiLwT6AV+q5RaYR77NtCulPqmKcwFSqk7kmlnLERJyz1Ar1Lqu8m0bSqIyFxgrlJqk4jkABuBa4FbsFm+TJCW67FfvgiQpZTqFREP8DJwO/Be4pwvTugBnAnsU0odUEoNAQ8C1yTZJsehlHoRaB9z+BrgN+b/v8EosJYnSlpsh1LqqFJqk/l/D7ATqMCG+TJBWmyHMug133rMlyIB+eIEAagAjox6X4dNHwwTBTwrIhtF5NZkGzNDypRSR8EowEBpku2ZKZ8VkbdNF5Hl3SajEZFqYCXwOjbPlzFpARvmi4i4RWQz0Aw8p5RKSL44QQAkwjE7+73OU0qdDlwBfMZ0R2iSz0+BRcBpwFHgv5NqzRQQkWzgz8DnlFLdybZnJkRIiy3zRSkVVEqdBlQCZ4rIikTcxwkCUAfMG/W+EmhIki0zRinVYP5tBh7FcHHZlSbTdxv24TYn2Z5po5RqMgttCPglNskX08f8Z+APSqlHzMO2zJdIabFrvoRRSnUCNcDlJCBfnCAAbwBLRGSBiHiBDwFPJNmmaSEiWeYAFyKSBbwL2DbxtyzNE8DN5v83A48n0ZYZES6YJtdhg3wxBxv/D9iplPreqFO2y5doabFpvpSISL75fwZwCbCLBORLykcBAZihXz8A3MCvlFJfS65F00NEFmK0+gHSgPvtkhYReQBYg7GkbRPwFeAx4CFgPnAY+IBSyvKDq1HSsgbDzaCAWuBTYX+tVRGR84GXgK1AyDz8RQzfua3yZYK0fBj75cspGIO8boxG+kNKqf8UkSLinC+OEACNRqPRjMcJLiCNRqPRREALgEaj0TgULQAajUbjULQAaDQajUPRAqDRaDQORQuAxhGISL6IfHrU+3IReThB9/KIyMZEXFujiSdaADROIR8YEQClVINS6v0Jutf5wKsJurZGEze0AGicwjeBReaa8N8Rkerwev4icouIPCYifxGRgyLyWRH5FxF5S0TWi0ih+blFIvI3cyG+l0TkxCj3uhz46+gD5uJe94nINjH2c/jnia4pImUi8qi5JvwWETk3Yb+MxrGkJdsAjWaWuBNYYS6wFV4xcjQrMFaQ9AH7gDuUUitF5PvATRgzyX8B/KNSaq+InAX8BLgowr0uBP5jzLHTgIpR+wfkm8ejXfN/gHVKqevMPS2yp5dsjSY6WgA0GoO15jryPSLSBfzFPL4VOMVcZfJc4E/GsjMApI+9iIiUY2za0T/m1AFgoYj8EHgKY0nvia55EYbwoJQKAl0zT6JGczxaADQag8FR/4dGvQ9hlBMX0BnuQUzAFcAzYw8qpTpE5FTgMuAzGDtVfS7Ga2o0CUGPAWicQg+QM90vm2vLHxSRD4Cx+qRZoY9lnP/f/Hwx4FJK/Rn4MnD6JNf8O/BP5nG3iORO13aNJhpaADSOQCnVBrxiDsJ+Z5qX+QjwDyKyBdjOmK1FTV/9EqXUrgjfrQBqzF2e7gPumuSatwMXishWjP1tl0/TZo0mKno1UI0mTphLEt+olPrHZNui0cSCFgCNRqNxKNoFpNFoNA5FC4BGo9E4FC0AGo1G41C0AGg0Go1D0QKg0Wg0DkULgEaj0TiU/w8Skpjgt55d4gAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def square(t, f=1, N=31):\n", " return (4/pi)*sum((N*sin(k*pi/N)/k/pi)*sin(2*k*f*pi*t)/k for k in range(1, N+1,2))\n", " \n", "u = lambda t: square(t, 0.1)\n", "\n", "first_order(5, 1, 30, u)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "colab": { "collapsed_sections": [], "name": "05.01-Response-of-a-First-Order-System-to-Step-and-Square-Wave-Inputs.ipynb", "provenance": [], "version": "0.3.2" }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 4 }