{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Blocksim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`tbcontrol.blocksim` is a simple library for simulating the kinds of block diagrams you would encounter in a typical undergraduate control textbook. Let's start with the most basic example of feedback control."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](../assets/simple_feedback.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tbcontrol\n",
    "tbcontrol.expectversion(\"0.1.1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tbcontrol import blocksim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our first job is to define objects representing each of the blocks. A common one is the LTI block"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gp = blocksim.LTI('Gp', 'u', 'y', 10, [100, 1], 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LTI: u →[ Gp ]→ y"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use a PI controller"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gc = blocksim.PI('Gc', 'e', 'u', 0.1, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PI: e →[ Gc ]→ u"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Gc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have the blocks, we can create a Diagram.\n",
    "\n",
    "Sums are specified as a dictionary with the keys being the output signal and the values being a tuple containing the input signals. The leading + is compulsory.\n",
    "\n",
    "The inputs come next and are specified as functions of time. `Blocksim.step()` can be used to build a step function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagram = blocksim.Diagram([Gp, Gc], \n",
    "                           sums={'e': ('+ysp', '-y')}, \n",
    "                           inputs={'ysp': blocksim.step()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LTI: u →[ Gp ]→ y\n",
       "PI: e →[ Gc ]→ u"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diagram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Blocksim is primarily focused on being able to simulate a diagram. The next step is to create a time vector and do the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time vector also specifies the step size for integration. Since blocksim uses Euler integration internally you should choose a time step which is at least 10 times smaller than the smallest time constant of all the blocks. The timespan is of course dependent on what you are investigating."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = numpy.arange(start=0, stop=1000, step=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "611961eabbb94bfd81876f7a6fb9fe1e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "simulation_results = diagram.simulate(ts, progress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of `simulate()` is a dictionary containing the simulation results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c19f872b0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3zU9f3A8df7LnfZg2zIYO8NAUQEUQRRcdbW0SpWHK2taG1tVVoH1lZ/Vquto+LWDkXUSl2oKIoIMsPeO4QMAmTPu8/vj+8BISRAkkvuLnk/H4/v4+6+830Z3/f3+/l+hhhjUEop1X7ZfB2AUkop39JEoJRS7ZwmAqWUauc0ESilVDuniUAppdq5IF8H0BTx8fGmS5cuvg5DKaUCyooVKw4YYxLqzg/IRNClSxeWL1/u6zCUUiqgiMju+uZr0ZBSSrVzmgiUUqqd00SglFLtnCYCpZRq5zQRKKVUO6eJQCml2jmvJAIReUVE8kRkXQPLx4tIoYhkeqb7ay2bLCKbRWSbiNzjjXiUUkqdPm+1I3gNeAZ44yTrLDTGTKk9Q0TswLPARCALWCYic40xG7wUlwJcbheZ+ZlsLNhIaXUpHUI60D++P/1i+yEivg5PKeVjXkkExphvRKRLEzYdCWwzxuwAEJG3gEsBTQReUOOuYc6WOby87mVySnNOWJ4WmcaNA27ksh6XEWQLyLaFSikvaM3//tEishrIBn5jjFkPpAB7a62TBYyqb2MRuQW4BSA9Pb2FQw18uaW5/PrrX7M6fzXDk4bz64xfk5GUQbQzmvzyfJbmLGX25tk8tPgh3t3yLo+MfYRu0d18HbZSygda62HxSqCzMWYw8Hfgv5759ZVL1DtkmjFmljEmwxiTkZBwQlcZqpY9RXu49qNr2XJoC4+OfZRXz3+VyV0mEx8aj8PuoFNEJy7rcRn/uvBfPD7ucbJKsrjmw2v4JusbX4eulPKBVkkExpgiY0yJ5/3HgENE4rHuANJqrZqKdcegmii7JJsb591IlbuKNy94k4u6XdTgcwARYXLXycy5eA6dozpz+5e387/t/2vliJVSvtYqiUBEksVzNhKRkZ7jFgDLgJ4i0lVEnMDVwNzWiKktqqip4I6v7qCsuowXJ71I79jep7VdUngSr01+jRHJI/j9ot8zb9e8Fo5UKeVPvPKMQET+A4wH4kUkC3gAcAAYY/4BXAn8XERqgHLgamOMAWpE5JfAPMAOvOJ5dqCa4E/f/4lNBzfx7IRn6RPbp1HbhjnC+Ns5f+NnX/yMe765h5jgGEZ1rPdxjVKqjRHrfBxYMjIyjHZDfbwFexdw+5e3c/PAm5k+bHqT91NcVcx1H1/HgYoD/Oei/5AWmXbqjZRSAUFEVhhjMurO15bFbUBxVTEPL3mYHjE9+PngnzdrX5HOSP527t8wxjD9y+mU15R7KUqllL/SRNAGPL/6efLL8pl55kwcdkez95celc7j4x5n2+Ft/GXZX7wQoVLKn2kiCHB7i/byn03/4YqeVzAwYaDX9ntmypn8tP9Pmb1lNvP3zPfafpVS/kcTQYB7etXTOGwOfjHkF17f9+1Db6dvbF8e+O4B8svyvb5/pZR/0EQQwLYc2sK8XfO4vt/1JIR5v5Gdw+7gsXGPUV5dzp+X/tnr+1dK+QdNBAHs5bUvExYUxnX9rmuxY3SN7srPh/ycz3d/zvzdWkSkVFukiSBA7S3ey6e7PuVHvX9EdHB0ix5rav+p9O7Qm0e+f4SiqqIWPZZSqvVpIghQr657FbvYub7f9S1+LIfNwUNnPkRBRQF/W/m3Fj+eUqp1aSIIQIWVhczdPpdLul/SIs8G6tM/vj/X9LmGd7a8w6aDm1rlmEqp1qGJIAB9sO0DKl2VXNPnmlY97s8H/5xoZzR//v7PBGKLdKVU/TQRBBi3cTN7y2yGJAw57U7lvCU6OJrpw6azMm8ln+76tFWPrZRqOZoIAsyS/UvYXbSbq/pc5ZPjX97jcvrG9uUvy/9CWXWZT2JQSnmXJoIAM2fLHDoEd2BS50k+Ob7dZufeUfeSV5bHS2tf8kkMSinv0kQQQAorC1mwdwEXdbsIp93psziGJg7lgq4X8OaGN8ktzfVZHEop79BEEEA+2/0Z1e5qpnSf4utQmD50Oi7j4rnVz/k6FKVUM3klEYjIKyKSJyLrGlj+YxFZ45m+E5HBtZbtEpG1IpIpIjrIwEl8uP1DukV3o19sP1+HQmpkKlf3uZr/bvsv2w5t83U4Sqlm8NYdwWvA5JMs3wmcbYwZBDwMzKqz/BxjzJD6BkxQlqziLFbmreTi7hc3OAZxa7tl4C2EB4Xz1MqnfB2KUqoZvJIIjDHfAAdPsvw7Y8whz8clWIPUq0b4aMdHAFzU9SIfR3JMTEgMNw26ia+zvmZZzjJfh6OUaiJfPCOYBnxS67MBPhORFSJyS0MbicgtIrJcRJbn57e/LpE/3vkxGUkZdIzo6OtQjnNtn2tJCkviieVP4DZuX4ejlGqCVk0EInIOViL4Xa3ZY4wxw4ALgF+IyLj6tjXGzDLGZBhjMhISWqdbBX+x4/AOdhTuYFIX31QZPZmQoBBuH3o76wvW89muz3wdjlKqCVotEYjIIOAl4FJjTMGR+caYbM9rHvA+MLK1YgoUn+/+HIAJ6RN8HEn9pnSbQs8OPXl65dNUu6p9HY5SqpFaJRGISDrwHnCdMWZLrfnhIhJ55D0wCai35lF7Nn/PfIYkDCExLNHXodTLbrNz1/C7yCrJYvaW2b4ORynVSN6qPvofYDHQW0SyRGSaiPxMRH7mWeV+IA54rk410STgWxFZDSwFPjLGaCc2tewt3svGgxs5r/N5vg7lpMZ0GsPI5JG8sPoFSqpKfB2OUqoRgryxE2PMSbvBNMbcBNxUz/wdwOATt1BHHBkVzF+LhY4QEe4afhdXf3Q1r65/lduH3u7rkJRSp0lbFvu5L/Z8Qd/YvqRG+n+N2/7x/ZncZTJvbnhTB7tXKoBoIvBjBeUFrMlfw7np5/o6lNM2feh0ql3VPL/6eV+HopQ6TZoI/Ni3+77FYBiXWm+NWr+UFpXGD3v/kPe2vsfOwp2+DkcpdRo0EfixhfsWkhCaQN/Yvr4OpVFuHXQrwfZgnl75tK9DUUqdBk0EfqraXc13+77jrJSz/KZvodMVFxrHTwf8lPl75pOZl+nrcJRSp6CJwE9l5mVSXF0cUMVCtV3f73riQuJ4csWTOr6xUn5OE4GfWrhvIUG2IM7oeIavQ2mSMEcYtw25jVV5q1iwd4Gvw1FKnYQmAj+1MGshwxOHE+GM8HUoTXZ5z8vpEtWFp1Y+RY27xtfhKKUaoInAD2WXZLPt8DbGpo71dSjN4rA5uGPYHewo3MEH2z7wdThKqQZoIvBD3+77FiDgEwFYLaIHJQziucznKK8p93U4Sql6eKWLCeVdS/YvITk8ma5RXX0dSrMd6Xrihk9v4J8b/snNg272dUgBK78sny2HtrDl0BayirPIK8sjrzyPsuoyqt3VVLmqCLIFERoUSmhQKLEhsSSHJ5McnkxqRCo9OvSga1RXHHaHr7+K8jOaCPyMy+1iac5Szkk7J+CqjTZkeNJwxqeO55V1r3BlryvpENLB1yEFhKKqIhZmLeT7/d+zNGcp+0r2HV0WExxDYlgiCWEJpEWmEWwPxmFzUO2uprymnLKaMvLL81l3YB2HKg8d3S7IFkSXqC4MjB/I0MShDE0cSueozm3mb001jSYCP7Pp0CYKKwsDtrZQQ+4cfidXzL2C5zKfY8YZM3wdjt+qdFXyxe4v+GTnJyzKXkSNu4YoZxQjk0fyk74/oXdsb3p16EV0cPRp77OipoK9xXvZemjr0TuKL/d+yfvb3gcgNiSW4UnDGdNpDGNSxpAcntxSX0/5KU0EfmZJ9hIARnUc5eNIvKt7THd+2OuHzN4ymyt7XUnv2N6+DsmvZJdkM3vzbN7b+h6HKg+RHJ7Mj/v8mEldJtE/rj92m73J+w4JCqFnh5707NCTC7kQALdxs6twF6vyVrEybyVL9i85OgBS9+jujEkZw5hOYxiePJxge7BXvqPyXxKIjX0yMjLM8uXLT71iALr5s5spqCjgvUve83UoXldYWciU96fQLbobr01+TYsjgL1Fe5m1dhb/2/4/DIbxqeO5pu81jEweiU1ary6HMYbth7ezKHsR3+77lhW5K6h2VxMaFMqo5FGMTR3L2JSxfjdmtmocEVlhjMmoO1/vCPxIRU0FK3NXclWfq3wdSouIDo5m+rDpzFw8k493fsxF3S7ydUg+k1eWxzOrnmHu9rkE2YK4ps81XN/vep+daEWEHh160KNDD6b2n0pZdRnLc5fzTdY3LMxayIKsBQD0iOnBuNRxjE0Zy+DEwThs+uC5LfBKIhCRV4ApQJ4xZkA9ywV4GrgQKANuMMas9CybCvzes+ofjTGveyOmQJSZn0mVu6rNPR+o7YoeVzBnyxyeXP4k56SdQ5gjzNchtapKVyVvrH+DF9e+SI27hqv7XM20AdNICEvwdWjHCXOEMS51HONSx2GMYWfhTisp7FvIG+vf4JV1rxDpiOTMlDMZmzKWMSljiA+N93XYActt3FS7q6l2VVuvnqm+hpidwjt5veaXt+4IXgOeAd5oYPkFQE/PNAp4HhglIrHAA0AGYIAVIjLXGHOogf20aUuylxAkQWQknXDn1mbYbXbuG3UfP/n4J/xjzT+4a/hdvg6p1SzMWsgj3z/CvpJ9nJt2Lr/J+A1pUWm+DuuURIRuMd3oFtONGwbcQElVCYv3L2Zh1kIW7lvIvF3zABgQN+BoEVL/+P6tWrTlL4wxFFcXk1uaS25ZLrmlueSV5ZFfnk9xVTEl1SWUVpdar1WlFFcXU1Zdhsu4TvsYH1z2Ad2iu3k1bq89IxCRLsCHDdwRvAAsMMb8x/N5MzD+yGSMubW+9RrS1GcEOTedT+WOrEZvd5yoFIhsmVoVGwo2YBMbfWL7tMj+/cnOop0UlBfQL64fYUFt+66gxtSwt2gvByoOEGoPJT0qnShnlK/D8goDlNWUUVh5mMLKQkqqSwGDw+Yg2hlNpDOSCGcEwfYQAv2JkAFq3NVUuaupdlVR5ao69t597LP7hJO64LAFYRe7Ndnsx723iR2bCIIgIgi24z7XVtmtEwMefpJIZ2STvoOvnxGkAHtrfc7yzGto/glE5BbgFoD09PSmRRESBc7wpm0LUFkMRftaJBHUmBpKa8pICe/k9X37o7TINAorC9lVuIu+cX2RgD9N1O9Q5SF2F+2mxl1Dp/BOdIzohK0NfVcBwoPCCA8Ko1N4J6rdNRRVFVJYWcjhysMcqDgAQJAEEeGMIMIRQZgjjNCgMBw2h9/8JAzmaKO8Klc11e76TvTVGNzHbScIDpsTp91BmCOMGJsTh92J0+bAaXfisDlx2B1e+50HR3dtchI4mdZKBPX9FMxJ5p8405hZwCyw7giaEkTyM+80ZbNjvp8Fn9wNt8+FuO7N21cdX+35ioe+ms5rkx+jc9Jwr+7bX23d9Sl3f303vx4+gRsG3ODrcLyqrLqMR75/hLnbV9E3ti8zx8xsF3d6tbmNm52FO8nMy2RFfiaZeZnsKtp1dHlMcAw9O/SkW3Q3UiNS6RTRiZTIFDqFdyImOMYrtcpcbpeVkMoPkF+eT35ZPvnl+eSV5XGg/AB5ZXnkluZyoOIAbnP8ST7YHkxSWBJJ4UkkhiVa7z2fk8OSSQpPIjYktk0UgbVWIsgCaheGpgLZnvnj68xf0EoxNV6vSVYi2PIpjP6FV3e9PHc5TpuTAfEnlKy1Wed3Pp+P0j7i2cxnmZA+ISDKy0/H5oOb+c3Xv2F30W5+Nvhn3DLolnZZu8YmNrrHdKd7THd+0OsHgFWF+Eijtq2HtrL10FY+3vExxdXFx21rFzvRwdFEB0cTExxDlDMKp91JkC0Ih81x9OdZ9wFrSXUJxVXFFFcVU1RVRGl1ab2xRTmjSAhNICEsgTNTzjx6gj96sg9LIjo4ut1UcW6tRDAX+KWIvIX1sLjQGLNfROYBfxKRI30OTALubaWYGq9DF0jo2yKJYEXuCgYlDGpXjXdEhN+P+j2XfXAZMxbN4JXzXyHIFrg1mo0xvLPlHR5b+hjRwdG8NOklRnYc6euw/Ep0cDQjkkcwInnEcfOLqorILslmX8k+skuyOVRxiMOVhznsefaQV5Z3rDaNq5oqdxU2bDjsjqPJIcgWRLgjnJSIFCKdkUQ5o6wp2DrpJ4YlEh8aT3xoPCFBIT76Cfgnb1Uf/Q/WlX28iGRh1QRyABhj/gF8jFV1dBtW9dGfepYdFJGHgWWeXc00xhz0Rkwtptf5sPgZqCiEkNNv5n8yJVUlbDy4kZsHtr8O2ZLCk/j9Gb/nnoX38OLaF/n54J/7OqQmKaoq4sHvHuTz3Z8zJmUMj4x5hLjQOF+HFTCinFFExUa1u+Izf+GVRGCMueYUyw1Q7yW0MeYV4BVvxNEqel8Ai56CbfNhwBVe2WVmfiZu4yYjue1WGz2Zi7pdxLf7vuWF1S8wuuNohiQO8XVIjbI2fy13f3M3uaW53DX8Lqb2n9omyo1V+6F/rY2VOgJCO8CWeV7b5fKc5QRJEIPiB3ltn4FmxqgZJIcnc8/CeyiuKj71Bn7Abdy8vv51rv/keowxvHbBa/x0wE81CaiAo3+xjWWzQ89JsPUzcJ9+I5CTWZG7gv7x/dtdK9vaIpwRPDr2UXJLc7l34b0n1ODwN4cqDvHL+b/kL8v/wvi08cy+eDaDEwb7OiylmkQTQVP0Oh/KD0LWslOvewrlNeWsK1jH8HZSZfRkhiQO4bcjf8vXWV/zzKpnfB1Og5blLOPKuVfy/f7vmTFqBk+Of7JR3UIr5W80ETRF9wkgdtj2RbN3tSZ/DTXumjbdrURjXN37aq7oeQUvrn3xaNcF/qLGXcOzmc9y02c3EeYI418X/Yur+1zdbqoYqrZLE0FThMZAagZs/7LZu1qeuxyb2BiaONQLgQU+EWHGqBkMSRjCfQvvY1lO8++6vCGnNIdp86bxj9X/YEq3Kbw95W2t4aLaDE0ETdX9XNi3EsqaV9t1Re4K+sT2IcIZ4aXAAp/T7uTv5/6d1MhUbv/ydjYUbPBpPPN3z+cHc3/ApoOb+PPYP/PIWY+06+c5qu3RRNBU3ScABnYsaPIuqlxVrMlfo88H6hETEsMLE18gyhnFrZ/fyvoD61s9hqKqIu5fdD93LriT1MhU3rn4HaZ0m9LqcSjV0jQRNFWnoVaDsu3zm7yLdQfWUemq1OcDDUgOT+blSS8T7ghn2mfTWrWY6Jusb7j8g8v5YPsHTBswjX9e8E/So5rY2aFSfk4TQVPZg6DbeNj+FTSxK+8VuSsAGJY4zHtxtTFpUWm8Pvl1ksKSuPXzW3lva8sO4Zlfls+9C+/lF/N/QZQzin9d+C/uHH6n1wcCUcqfaCJoju7nWt1SH9jSpM0z8zPpHt2dmJAYLwfWtiSFJ/H65NfJSMrgge8e4MHvHqSsusyrx6hyVfHquleZ8v4U5u2ax80Db+btKW+3q04AVfuliaA5up9rvW5rfPGQ27hZnb864LpT8JWYkBieP+95pg2Yxrtb3+WKuVewaN+iZu+3ylXF25ve5sL3LuTJFU8yInkE/730v0wfNh2n3emFyJXyf4Hb1aM/iEmHuJ5WNdLRtzVq012FuyisLNTWqI1gt9m5c/idnJVyFg8ufpCfffEzRnUcxa2DbiUjKaNR9fmzS7KZs2UO7297nwPlBxiSMISZY2ZyZqczW/AbKOWfNBE0V/dzYeUbUF0BjtPv2jYzPxNA7wiaICM5g3cveZfZm2fz0tqXuHHejXSJ6sL5Xc5nVMdR9Intc8IoTkf6wV+Vt4qv937N2gNrARibOpYf9/0xozuO1oZh7VR1dTVZWVlUVFT4OhSvCQkJITU1FYfj9J5taSJorh4TYOkLsGcxdD/ntDdblbeKmOAYukR1abnY2rBgezDX9buOH/T8AZ/v/pz3tr7Hi2tf5IU1LwDQIbjD0br+RZVFxw18MiBuALcNuY1Lul9Cp4j2MTSoalhWVhaRkZF06dKlTVwMGGMoKCggKyuLrl27ntY2mgiaq/MYsAXBzq8blQgy8zIZkjCkTfzh+VKYI4xLe1zKpT0upaiqiJW5K9lRuIN9xfuocFXgMi5igmNIDEukT4c+9InrQ2xIrK/DVn6koqKizSQBsFrnx8XFkZ+ff9rbeGtgmsnA04AdeMkY82id5X8Fjpwlw4BEY0yMZ5kLWOtZtscYc4k3Ymo1wRGQkgE7F572JocqDrGraBeX9ri0BQNrf6KcUYxPG8/4tPG+DkUFmLaSBI5o7PdpdiIQETvwLDARawziZSIy1xhztF8AY8yvaq1/O1C7Y51yY0xgF5R3HQsLn4SKIgiJOuXqq/NXA2j/Qkopv+CN6qMjgW3GmB3GmCrgLeBkl7rXAP/xwnH9R5exYFzWc4LTkJmXSZAE0T+ufwsHppRSp+aNRJAC7K31Ocsz7wQi0hnoCtTutjNERJaLyBIRuayhg4jILZ71ljem7KtVpI0EezDs/Oa0Vs/Mz6RvXF8dQFspBcCuXbsYMOBY48W//OUvPPjgg612fG88I6ivMKqhPheuBuYYY2oP7ZVujMkWkW7AlyKy1hiz/YQdGjMLmAWQkZHRtD4dWooj1EoGp5EIql3VrDuwjh/1/lErBKaUaoyH/reeDdlFXt1nv05RPHCxf9/9e+OOIAtIq/U5FchuYN2rqVMsZIzJ9rzuABZw/PODwNFlLOSsPWW31JsObqLSVcmQhMB+LKKUaju8cUewDOgpIl2BfVgn+2vrriQivYEOwOJa8zoAZcaYShGJB8YA/+eFmFpf13Gw4E+w+zvo23BXxavyVgHakEwpf+SrK/egoCDc7mPjdLd247Zm3xEYY2qAXwLzgI3AbGPMehGZKSK1q4JeA7xlzHFddfYFlovIauAr4NHatY0CSspwcISdsngoMz+TTuGdSAxLbKXAlFL+Likpiby8PAoKCqisrOTDDz9s1eN7pR2BMeZj4OM68+6v8/nBerb7DhjojRh8LsgJaaNgV8PtCYwxrM5bTUayjj+glDrG4XBw//33M2rUKLp27UqfPq07DKq2LPamrmNh/kwoyYeIhBMWZ5dmk1eep8VCSqkTTJ8+nenTp/vk2NoNtTd1Pdt6beCuIDPP09GcPihWSvkRTQTe1HEIOCNPmghCg0Lp2aFnKwemlFIN00TgTfYg6Dwadn1b7+LV+asZFD+IIJuWyCml/IcmAm9LH20NXVl64LjZZdVlbDm0hUEJg3wUmFJK1U8Tgbd19oxwVaffofUF63EZlz4oVkr5HU0E3tZpKASFwO7jE8GRB8U6NKVSyt9oIvC2oGBrfII93x03e3X+arpGdyU6ONpHgSmlVP00EbSEzqNh/xqoLAE8DcnyV+vdgFLKL2kiaAnpo63xCbKWArC7aDeHKw9rIlBKNeif//wnI0eOZMiQIdx66624XK5Tb+QlWo+xJaSNBLFZHdB1P/foiGTakEwpP/fJPVYvwt6UPBAuePSkq2zcuJG3336bRYsW4XA4uO222/jXv/7F9ddf791YGqCJoCUER0LyoKMPjDPzM4l0RNItppuPA1NK+aP58+ezYsUKRowYAUB5eTmJia3XMaUmgpbS+UxY/grUVFoNyRIGYRMtiVPKr53iyr2lGGOYOnUqf/7zn31yfD0ztZT00VBTQfGe79h2aJs+H1BKNWjChAnMmTOHvLw8AA4ePMju3btb7fiaCFpK+mgA1m77CINhcKImAqVU/fr168cf//hHJk2axKBBg5g4cSL79+9vteNr0VBLiUiAuJ6szl2BIAyK164llFINu+qqq7jqqqt8cmyv3BGIyGQR2Swi20TknnqW3yAi+SKS6ZluqrVsqohs9UxTvRGP3+g8mtVl++kR04MIZ4Svo1FKqXo1+45AROzAs8BErIHsl4nI3HqGnHzbGPPLOtvGAg8AGYABVni2PdTcuPyBO200aw58yfkR6b4ORSmlGuSNO4KRwDZjzA5jTBXwFnDpaW57PvC5Meag5+T/OTDZCzH5hR2xKRTbbQxx6aMYpZT/8sYZKgXYW+tzlmdeXT8QkTUiMkdE0hq5LSJyi4gsF5Hl+fn5Xgi75WVWWDUABhfm+TgSpZRqmDcSgdQzz9T5/D+gizFmEPAF8HojtrVmGjPLGJNhjMlISDhxPGB/tDJvFbEE0Xnfal+HopRSDfJGIsgC0mp9TgWya69gjCkwxlR6Pr4IDD/dbQPZyryVDAtPRQ7vgeJcX4ejlFL18kYiWAb0FJGuIuIErgbm1l5BRDrW+ngJsNHzfh4wSUQ6iEgHYJJnXsDLKc1hX8k+hnU6w5rh6YBOKaX8TbMTgTGmBvgl1gl8IzDbGLNeRGaKyCWe1aaLyHoRWQ1MB27wbHsQeBgrmSwDZnrmBbyVuSsBGN5jCtidsFcTgVLKP3mlQZkx5mPg4zrz7q/1/l7g3ga2fQV4xRtx+JOVeSsJd4TTO34AdBwMWct8HZJSyk/94Q9/ID4+njvuuAOAGTNmkJSUxPTp01vl+NqyuIWsyF3BkIQh2G12SBsFS1+EmioIcvo6NKVUAx5b+hibDm7y6j77xPbhdyN/d9J1pk2bxhVXXMEdd9yB2+3mrbfeYunS1itF0AruLaCwspBth7cxLGmYNSN1BLgqvd/PuVKqTejSpQtxcXGsWrWKzz77jKFDhxIXF9dqx9c7ghawKm8VAMMSPYkgbZT1uvd7SB3ewFZKKV871ZV7S7rpppt47bXXyMnJ4cYbb2zVY+sdQQtYkbsCh83BwISB1oyojhCdpjWHlFINuvzyy/n0009ZtmwZ559/fqseW+8IWsDK3JUMiB9AsD342MzUEdYdgVJK1cPpdHLOOecQExOD3W5v1WPrHYGXlVWXsaFgw7FioSPSRkHRPijM8k1gSim/5na7WbJkCdOmTWv1Y2si8LK1B9ZSY2qOPSg+Is0ai1TbEyil6tqwYQM9evRgwoQJ9OzZs9WPr0VDXrYydyWCMCRxyPELkgdBUKjVnmDAFb4JTinll/r168eOHTt8dgRyO5UAACAASURBVHy9I/CyFXkr6NWhF1HOqOMX2B3Qaag+J1BK+R1NBF5U7a5mTf6aE4uFjkgbCfvXQHV56wamlFInoYnAi9YfWE95TTkZSRn1r5A2EtzVkJ3ZuoEppdRJaCLwoqU51oPgEckj6l8hdaT1qsVDSik/oonAi5bmLKVXh150COlQ/woRCRDbTWsOKaX8iiYCL6lyVZGZl8nI5JEnXzFtlHVHYOodiE0ppVqdJgIvWZ2/mkpX5eklgrIDcNB3VcWUUv7lD3/4A08//fTRzzNmzOChhx5i3LhxDBkyhAEDBrBw4UIAIiIi+PWvf82wYcOYMGEC3hjDXdsReMmynGXYxMbw5FN0Kne0A7qlENe95QNTSp22nD/9icqN3u2GOrhvH5Lvu++k69TXDfVtt93G+eefz4wZM3C5XJSVlQFQWlrKsGHDeOKJJ5g5cyYPPfQQzzzzTLNi9ModgYhMFpHNIrJNRO6pZ/ldIrJBRNaIyHwR6VxrmUtEMj3T3LrbBoqlOUvpE9vnxPYDdSX0geBo2LukdQJTSvm9+rqhHjFiBK+++ioPPvgga9euJTIyEgCbzcZVV10FwE9+8hO+/fbbZh+/2XcEImIHngUmYg1Gv0xE5hpjNtRabRWQYYwpE5GfA/8HXOVZVm6MqdMMN7CU15SzOn811/W97tQr22xWdxP6wFgpv3OqK/eWVLcb6nHjxvHNN9/w0Ucfcd1113H33Xdz/fXXn7CdiDT72N64IxgJbDPG7DDGVAFvAZfWXsEY85UxpszzcQmQ6oXj+o3MvExq3DUNVxutK20U5G2E8sMtG5hSKmDU7YZ69+7dJCYmcvPNNzNt2jRWrrTGQXe73cyZMweAf//735x11lnNPrY3nhGkAHtrfc4CRp1k/WnAJ7U+h4jIcqAGeNQY89/6NhKRW4BbANLT05sVsLcty1mGXewNtyiuK20UYCBrOfQ8r0VjU0oFhrrdUC9YsIDHH38ch8NBREQEb7zxBgDh4eGsX7+e4cOHEx0dzdtvv93sY3sjEdR3X1Jv3UgR+QmQAZxda3a6MSZbRLoBX4rIWmPM9hN2aMwsYBZARkaGX9W9/D7ne/rHDaCk3M7+Q8XUuA3GQERwEDFhDiKCg46/fUsZDmKzqpFqIlBKcawb6nfeeQeAqVOnMnXq1HrXffjhh3n44Ye9dmxvJIIsIK3W51Qgu+5KInIeMAM42xhTeWS+MSbb87pDRBYAQ4ETEoE/KSyvZvH2A2TuLWTd/lzW2NZRVXA2oxbOr3f9iOAguieE0yMxkowuHRjZNZZuSQMQbWGslMLqhnrKlClcfvnlAdsN9TKgp4h0BfYBVwPX1l5BRIYCLwCTjTF5teZ3AMqMMZUiEg+MwXqQ7HeMMXyQmc0bi3eRufcwbgMOu5DSaQ9EuLm411gGjxlAdKgDh826+i+urOFwWRX7DpWzPb+UBZvzeHelNTDNY2EpXMYCFm/IZnSvJIKDWndEojartABy11rtNCoKoaYKgiMhIhES+0J8L6snWKX8SGO6oS4pKfH68ZudCIwxNSLyS2AeYAdeMcasF5GZwHJjzFzgcSACeMdTRLLHGHMJ0Bd4QUTcWA+uH61T28gv7D1Yxoz/ruObLfn0TorkF+f0YGzPBAanRfP0qid4e5OTP15wMSFBISfdjzGGHQdKWbrzIKUrhxOc8wmPv/k+uxzdOb9/MldmpHJG1zhstubXAmhXCvdB5r9g41zIWXvyde3B0GUM9JoM/S6DyKTWiVEpPyYmALs6yMjIMMuXL2/x47jchlcX7eSJz7ZgE/jt5D785IzO2GudqC//4HISQhOYNWlW43Z+eC88NYAtwx/g5crz+Hjtfoora0iLDeUHw1K5akQaHaNDvfyN2piC7bDgUVg3B4wb0s+EHuda40PH9YCwOLA5oLIIinMgbwPsWwFbP4OCbWALgj4XwYiboMtY8EI1PBV4Nm7cSJ8+fbxSDdNfGGPYtGkTffv2PW6+iKwwxpzQPbImggZsyC7i3vfWsDqrkHP7JPLHywbQKeb4E3NOaQ4T50zk18N/zQ0DbmjcAYyBJ/tB5zPhypcpr3Ixb30O76zYy6JtBdhtwuQBydw4pgvD0ju0qT/SZqsuh68fg0V/A7sTRt4EGdMgtuvp7yN/M6x8AzL/DeUHIe0MGH8PdBuvCaGd2blzJ5GRkcTFxbWJ/zNjDAUFBRQXF9O16/H/Ew0lAu1ioo6Kahd/m7+VWd/sICbMwd+vGcqUQR3r/QNZnL0YgDNTzmz8gUSs8Qk8DctCnXYuG5rCZUNT2HuwjDcW7+KtZXv5aM1+BqZE89MxXbhoUEd9lpC7AWZfDwVbYciPYcIDTSveSegN5z8C5/4BVr0JC5+ENy+D7hNg8qOQ0Mv7sSu/lJqaSlZWllf67PEXISEhpKaefnMtvSOoZfH2Au57fy07D5Tyw+GpzLioLzFhzgbXv/vru1mRu4L5P5zftCuJJc/Dp/fAXRshqtMJi0sra3hv1T5eW7ST7fmlJEQGc90ZnfnxqHTiIoIbf7xAt3YOzL3devh7+T+g+7ne23dNJSx7ySpqqi6DUT+Dc+4DZ7j3jqGUjzV0R6C9jwKFZdXc8+4arnlxCS634Z/TRvH4DwefNAm43C4W71/M6E6jm347mebpqXRP/f0OhQcHcd0Znfn8V2fz+o0j6dcxiic/38KZj37JPe+uYUtucdOOG4gWPgnvToOOg+HWb7ybBACCgmH0L+D2lTD4alj8DDw/BnYt8u5xlPJD7bpoyBjDp+tyuH/ueg6WVnHruG7ceV4vQp2nLn7ZeHAjhZWFnNmpCcVCRyQPAkc47FkMA65ocDWbTTi7VwJn90pgW14xryzaxXsrs3hr2V7G9oznxrO6cnbPhLZZ28gY+OJBWPQUDPwhXPocBDWcoJstIgEufRYGXwMf/AJeuxBG3grnPaB3B6rNareJIKewgvs/WMdnG3Lp3ymKV28YwYCU6NPe/rvs7wA4o+MZTQ/C7oD0UY266uyRGMmfLh/I3ZN68++le3j9u1389NVldE8I58azunLF0NTTSmQB40gSyLgRLnzC6rSvhVTWuCipqKG4ooZS50DcF39C4rLHSFr6AmWbPmfDmKcoielLkM1GkF0IsglBdhtBNsFhtxHmtBPmtBMeHERwkK1NPHhU7UO7e0bgdhv+vXQPj32yiSqXm7sm9mLaWV0JsjfuBHPDpzdQVl3G7ItnNymOo755HL78I/x2J4TFNnrzqho3H63N5uVvd7JuXxEdwhxcOyqd60d3ISnq5O0a/N6ip+Hz+60kcNGTTarNU+1yk324nD0Hy8gprCCvuJJ8z5RXXMGBkiqKyqsprqyhqsZd7z5G29bzlONZYijhkZof84ZrEvX3rHKM3SaEe5JCmNNOZIiDmDAHMaEOYsKcRId6Poc5iAl1El1rWVRIUKP/HpU6HVprCNiWV8K9761h2a5DjOkRx58uH0jnuMbf7pdUlbA6bzVT+9ffD0ijdB5jve5ZbNVpbyRnkI3Lh6Zy2ZAUlu48yCuLdvLcgu288PUOJvRN5IfD0zi7dwKOQDuxrH7LSgL9r4AL/3LSJGCMIetQOVtyi9maV8LugjL2HCxlz8Eysg9X4HIff7ETGRxEQlQwCRHB9O8URXSog4iQIKJCrH6hIoKDCA+247DbPFf8o8iquozg737HzKzX+VW3fWwd/SgVQTHUuN3UuAzVLkNZVQ1lVS5KKmsoq6qhtNJFaaU1r6iimkOlVezIL+VwWRVFFTUn/fqRnniiQ49NUaFBdT5bU3So4+i6EZ67kUAoJnS5DdUuN1UuN1U1buv90VdDlevYvCqXm+ojr0fnmWPzPNtVuty4XAa3AYPV55fbGM9k/a243Rz/2fMewCZWt84iYBNBsF6tG1HBdmR+rVfxzD8yz2YT7LVfxSretXs+i1gXCnabYJMjrxx9f2S+9Zla61j7zOjcgfBg756621UiuO+9tWzJLeHxKwdx5fDUJt+6L8tZRo2pYUzKmOYHlTLcau26+7smJYIjRIRR3eIY1S2OPQVlvLlkF++v2se89bnERwRz+dBOXDEslT7Jkf5fZJG1AuZOtxp5Xf4C2I4VdR0srWJzTjGbc4rYnFvM5pxituSWUFJ57MQaG+4kLTaMoWkduHRwGOlxYaTHhtEpOpSEyOAmFp3FQ9/3YMnzdPj8fkZ+cgn84CWrlXITuNyGovJqDpdXc7isisPl1RSWHXt/uKyaovJqCj3TjgMlR99XVNd/51JbqMNOqNN+9DXMaSfEYb067VbRlt1mFWvZbVLn1YZNwO05iRrPibLuydMY63tUnuokXedEb80zJyRob3AGWd+p9sn62Mn72PsjJ3yb7dgJH6zeMo98N2OOvT+WLEytn0vtdWv/jKzv1gJfD4Av7jqbHokRXt1nuyoa2nWglPDgIBIim1f1cubimXy04yMWXr0Qp90LDy5fvQiqSuDWr5u/r1qqXW6+3pzPOyv2Mn9jHjVuQ5e4MM4fkMzk/skMTo3xvyvH4lyYNR63LYiNF3/A+sMOz8m+mE05xeQXH+2vkJgwB72TIumdbE19kiPpkRhJdGgL9yWUnQlzfgqHdsH4+2DsXcclq5ZWWeOiqLzmaGIoKq+mqMJ6X1rporzaRXlVDeXVLsqqXFR4XsurrGVVNW5cbutkVXP01X3ss8vgMubo1avt6BUqR0+mR66CbTbBGWTDabcdfXXYbTiOzpN65tlw2AWn3Y7Ds7zhbe047HLctk7POg67EFxrH3ab+NVFjjmSLI8mhmMJwu22fsZuz+cj712e+cYYXG7qbGdN/TtFE+Jo2t+btiz2EmMM5805j8EJg3ly/JPe2elXf7KeFfxuN4ScYqjLJiooqWTe+lw+XZ/Dd9sOUOM2xEcEc2b3OM7sHseYHvGkdght9X+kqho3uwpK2ZxTzNacQ1y86hZSK7byg6oH2eC2RjQNDrLRMymC3klR9EmOpJfnpJ8YGey7f/zKYvjwV7D2Heg6Dq54ESKTfROLUqdJnxF4yaaDm8gry2Nc6jjv7bTzmVZfOXu/h54TvbffWuIigrl2VDrXjkqnsLyaLzfl8vXmfBZtL2DuaqvX8PiIYAamRDEgJZr+naLoEh9OemwYYc7m/ZlU1rjIKaxgd0EZuw+WsaeglN0FZewqKGXngVKqXdbFyK8c79LTvo5XkmYwscd53O650u8cF35c/05+ITjSOvl3Gw8f/Qb+cZZVjNVjgq8jU6rRNBE00tdZXyMIY1PGem+nqSOtDtB2L2qxRFBbdKiDy4emcvnQVIwxbM8v4bvtBazeW8i6fYV8vSX/uPLN+IhgOkaHEBPmIDbcSUyog2CH/ejDL5tNqKxxUVFlFUGUVbsoKLFq5hwoqaKwvPq44wcH2UiPDaNzXDjn9U2id3Ikg90b6fy/92HQNdx4+W9b/GfgFSIw9CeQkmEVFf3zCjjrV3DODO3qWgUUTQSN9E3WNwxMGEhcaJz3duoMg07DfNKKVUTokWiVrTPamlde5WJrXrGn9k0ZewrKyC2u4FBZNbsLyjhcVkWV61g5s9uA0247+lAy1GknLtxJ7+RIxkRYtXOSokJIjwujc1wYSZEhxz+bKD8M//gVxKTDhY+3+s+g2RL7wE3zre5Cvv2r9eD/By9DTNqpt1XKD2giaIQD5QdYe2Attw+93fs77zIGvvs7VJVZicGHQp12BqXGMCg15rTWN8Y0r6z+k99B0T6Y9plV5BKInGFwyd+s5wX/u9MqKrrsuWbVBGt33C4oyobDe6AkB8oOWlP5QetioaYCXNXgqrIm47buvOxO647a7gRHGITGQEi0Z4qxPocnQGRH69Wup726vPITEZHJwNNYA9O8ZIx5tM7yYOANYDhQAFxljNnlWXYv1oD2LmC6MWaeN2JqCQuzFgJwdurZp1izCTqPsa4ms5Za5c4BpFlJYOsXsOYtGPdbSD3hGVbgGXgldBoKc26Et661uqeYOBMcAd64z5sqSyB/E+Sus3qTzd8Ih3ZbFwPuetpXOCMhtAM4Qq2T/ZGTv9g8yaHqWIKoKrNGpqsurf/YYrOSQUSSlRgikyEqBaI6Wh0/RqVYr8FR7ao78mYnAhGxA88CE7HGL14mInPrjDQ2DThkjOkhIlcDjwFXiUg/rKEt+wOdgC9EpJcxxtXcuFrC11lfkxyeTK8OLdBFcdoo649016KASwRNVlkMH94J8b1h3G98HY33xHW37m6+eBCWPAc7Flj9F6WN8HVkrc/thgNbrIoQR6aCbceWOyOsLsFTR0CHK63iwZh06yQdFmdd0TelbylXtZUQKgqh/BCU5ELxfqt6ckmONVBR8X7IXgml9XQ/7YywEkJkx2PJ4Wii8MwLi2szycIbdwQjgW3GmB0AIvIWcClQOxFcCjzoeT8HeEasy8hLgbc8g9nvFJFtnv0t9kJcJ/hox0fsKdrT5O2/y/6OS7pf0jJVFkOirJ41d33r/X37q/kzoTDLOmkGtbFutYOCYfKfrVpEc++AlydavZueM8PnRX8trjALtn5uTbu/tU7GAKGx1gXPwB9BUn9riuncMv1H2R0QHm9Np1JTZSWFomwozrZei7KtO5SibNj5jbW87vWpPfhYUjghaXheIxJbtY1JU3kjEaQAe2t9zgJGNbSOZ4zjQiDOM39JnW1T6juIiNwC3AKQnp7epEA/2vERC/ctbNK2AA6bgwu6XtDk7U+p69mw+Fnr1jnYuy0H/c7epbD0RRh5y7HuuNuiHufBbYutu4PFz8Cmj6wuM3qe5+vIvMdVYxVpbplnnfzz1lvzo9Oh36XW6G9po6w7JX+8gg5yQofO1tQQtwtK8o4liOL9xxJFUTZkLbdeXZXHbyd2T/FTJ6t4KzjKeg4WHHn8e0eI51mH41jR15FiMKmTKBP6eL2o0RuJoL7fbN1Wag2tczrbWjONmQXMAqtBWWMCPOLZCc82ZbPjtGgDpm7jrZ429yxulWqkPuN2wUd3WVdQE/7g62haXkgUTHnS6mp87nT41w+g5/lw/p8gvoevo2uaknzY9gVsnQfbv7Su+m1BkD4aJj4MPSdZRT7+eOJvCpvdc/XfEetRZz2MsR5uH0kQde8uSvOtcbYri62pprxpsfximddH0PNGIsgCateTSwWyG1gnS0SCgGjg4Glu6zX+1Py8XulnWLebOxa07USw4jXIWQtXvhK4tYSaostZcNsS+P4f8PX/wXNnwIhpVtsDf2+V7HbD/lXWFf+WeZC9CjAQngh9LoZek6wLmZDT78q9zRGB8Dhr6jjo1Ou7qj1JocgqnjpSG6p2zShXlZVgaovq6PXQvZEIlgE9RaQrsA/r4e+1ddaZC0zFKvu/EvjSGGNEZC7wbxF5EuthcU9gqRdiCkyOUCsZ7Fjg60haTtlB+PJhq0O5/g0PxtNmBTlhzHQYdJX1c1j6opUYM26EMXf4V0IoP2xd7W/9zEoAZQcAsWp3nXOfddWfPKhFx4ho0+wOq+v5JnQ/723NTgSeMv9fAvOwqo++YoxZLyIzgeXGmLnAy8CbnofBB7GSBZ71ZmM9WK4BfuGvNYZaTbfxMP8hqzwyItHX0Xjflw9DRRFc8FjbKTZoisgkuPQZ625g4RPw/QtWUuh3CYy42bogaO2fz5Gr/m3zrWKfrOXWA9KQGOtZR89J1mu4FxtTKr+gnc75m30r4cVzrJapA6/0dTTelZ0Js8bDqFutRKCOObgDlr4Emf+0ytsT+lh3TP0v93p58FFHqnbuWWzVVtvxFZQVWMs6DYXuE6wiypQMbYTVRmjvo4HC7YL/6wZ9p1h1z9sKY+DVC60Tz+0rrNae6kRVpbB2Dqx52+qqAgMdulotlruOs07QHbo0vkqi2w2HdkLeBqsRV/ZK2LMEKg5by8MTofu51hV/93NOr9qlCjja+2igsNmtf/jtC6yTZ1spPtn8Cez5Di56QpPAyTjDYfhUayraDxv/Z12pr38fVr5ureMIg/ieVl31yI5WGbM92NMWw1ita6tLoewQFGVZtVYKs6C6zHMQgbge0Pdiq5ZP+hkQ263t/K2pRtNE4I+6jYeNc63igrjuvo6m+Vw1Vj36uB4wzAvDe7YXUR1h1C3W5KrxdMmwDnLXW3dWh3ZbxTrlhzmh1nVQqFWDJzrFKmbqMdHqHC+xv/XqbPwQrart0kTgj7qNt163f9k2EsHqf8OBzfCjN7R75qayB0GnIdZUlzFWHz01FYBYdwxak0c1gv61+KPYblY58NbPfR1J81WVWSOwpWRA30t8HU3bJGIl2OBIq0W6JgHVSPoX449ErJanO7+B6ia2PvQX3z9vNcefOFPLoJXyU5oI/FWvSVYT9EDuhK60AL59CnpNtsZbUEr5JU0E/qrzWVZZ7xa/HZ7h1BY+AVUlcN6Dvo5EKXUSmgj8lSPE6o1067wT+xoJBId2w7IXYci1kNjX19EopU5CE4E/6zXJGrbvwBZfR9J4X/7R6j53/H2+jkQpdQqaCPxZD08PpIFWPLR/NaydDWf83KrHrpTya5oI/FlMmtUAaOtnvo6kcb540BqEY8ydvo5EKXUaNBH4u96TrT5nyg76OpLTs/0rqyHc2N9oVxJKBQhNBP6u78VWV8CbP/Z1JKfmdsPn91tDFI682dfRKKVOkyYCf9dxiHVi3TDX15Gc2vr3IGcNnPv7tjcYvVJtmCYCfydi3RXs+Moa0MVf1VTC/JmQNBAG/tDX0SilGqFZiUBEYkXkcxHZ6nntUM86Q0RksYisF5E1InJVrWWvichOEcn0TPX0qKXod4k1dqk/PzRe/goc3g0TH9S+bpQKMM39j70HmG+M6QnM93yuqwy43hjTH5gMPCUitZ8i3m2MGeKZMpsZT9uUOhIikqyuqf1RRaE1GHu38daoVkqpgNLcRHAp4Bktg9eBy+quYIzZYozZ6nmfDeQBCc08bvtis0GfKVZvpFWlvo7mRN8+BeUHtWM5pQJUcxNBkjFmP4Dn9aSjrYvISMAJbK81+xFPkdFfRaTBJ4wicouILBeR5fn5+c0MOwANvNIaYWrjh76O5HiF+2DJczDoKug42NfRKKWa4JSJQES+EJF19UyXNuZAItIReBP4qTHG7Zl9L9AHGAHEAr9raHtjzCxjTIYxJiMhoR3eUKSdYdUeWvO2ryM53ld/AuOGc2b4OhKlVBOdcoQyY8x5DS0TkVwR6WiM2e850ec1sF4U8BHwe2PMklr73u95WykirwK/aVT07YnNBoN+BN8+CcU5EJns64isIRNX/xvOuA06dPZ1NEqpJmpu0dBc4MggtFOBD+quICJO4H3gDWPMO3WWdfS8CtbzhXXNjKdtG3y1dfW9do6vI7F88aA1KtbYX/s6EqVUMzQ3ETwKTBSRrcBEz2dEJENEXvKs8yNgHHBDPdVE/yUia4G1QDzwx2bG07bF94ROw2DNW76OxBo9betnVhIIi/V1NEqpZmjW4PXGmALghPqCxpjlwE2e9/8E/tnA9uc25/jt0uCr4ZPfQs5aSB7omxhcNfDpvRCdBiNv9U0MSimv0ZY/gWbgDyEoBJa97LsYVr4Guetg0h+tAXSUUgFNE0GgCYuFAVfCmtlWQ67WVnbQGnSmy1jo16iKY0opP6WJIBCNmAbVpbDaB1VJv3rE6vPogse08ZhSbYQmgkCUMsx6aLzspdYdzzhnrdWn0IibIKl/6x1XKdWiNBEEqlE/gwObW28YS7cLPrwLQmLgnHtb55hKqVahiSBQDbjCamm88InWuStY9jJkLYXJj1rDUCql2gxNBIHK7oAx062T8+5FLXusw3th/kNWz6KDftSyx1JKtTpNBIFs6E8gPAG+ebzljmEMfHSX1aJ5yl/1AbFSbZAmgkDmCIWzfgU7FsC2+S1zjBWvWi2IJ9yv/Qkp1UZpIgh0I26CmM7WoPFul3f3nbcJPr0Pup+rLYiVasM0EQS6oGA47wGrpe+K17y33+oKePcmcIbBZc/r8JNKtWH6390W9L8Cuo6zegMtym7+/oyB/90BuWvh0uf8o8trpVSL0UTQFojAxU+Dqxr+d2fzq5N+93erh9Px90Hvyd6JUSnltzQRtBWx3WDiQ7B1Hix6qun7Wf229byh7yUw7m7vxaeU8luaCNqSkbdYxUTzZzatxfH6/8J/fwZdx8IVs/S5gFLtRLP+00UkVkQ+F5Gtntd6m5yKiKvWoDRza83vKiLfe7Z/2zOamWoqEbjk79Y4BW9f17gqpUv+AXN+Cqkj4Or/WFVTlVLtQnMv+e4B5htjegLzPZ/rU26MGeKZLqk1/zHgr57tDwHTmhmPCo6A6/4L8b3g31fB0hdP/syg7CDMuRE+/R30ugCue9/ah1Kq3WhuIrgUeN3z/nWscYdPi2ec4nOBIwPwNmp7dRJhsXDD/6DbePj4N/DqhbDlM+th8hGH98DXj8Pfh8GGD+CcGXDVm+AM91XUSikfadZQlUCSMWY/gDFmv4gkNrBeiIgsB2qAR40x/wXigMPGmBrPOllASkMHEpFbgFsA0tPTmxl2OxDaAa6dDavetAaS+fcPwe6EyI5QVQplB6z1ek6yWg37athLpZTPnTIRiMgXQH0VyWc04jjpxphsEekGfOkZsL6onvUaLMMwxswCZgFkZGS0Yif8Acxmg+FTYfA1sO1z2LMYSvKsRmgJfa2qobHdfB2lUsrHTpkIjDHnNbRMRHJFpKPnbqAjkNfAPrI9rztEZAEwFHgXiBGRIM9dQSrghdZQ6gRBTuhzkTUppVQdzX1GMBeY6nk/Ffig7goi0kFEgj3v44ExwAZjjAG+Aq482fZKKaVaVnMTwaPARBHZCkz0fEZEMkTkJc86fYHlIrIa68T/qDFmg2fZ74C7RGQb1jODl5sZj1JKqUYS05pj3npJRkaGWb58ua/DUEqpgCIiK4wxGXXna9NRpZRq5zQRKKVUJcTk1gAABKlJREFUO6eJQCml2jlNBEop1c5pIlBKqXYuIGsNiUg+sLuJm8cDB7wYTiDQ79w+6HduH5rznTsbYxLqzgzIRNAcIrK8vupTbZl+5/ZBv3P70BLfWYuGlFKqndNEoJRS7Vx7TASzfB2AD+h3bh/0O7cPXv/O7e4ZgVJKqeO1xzsCpZRStWgiUEqpdq5dJQIRmSwim0Vkm4jc4+t4vEFE0kTkKxHZKCLrReQOz/xYEflcRLZ6Xjt45ouI/M3zM1gjIsN8+w2aTkTsIrJKRD70fO4qIt97vvPbIuL0zA/2fN7mWd7Fl3E3lYjEiMgcEdnk+X2Pbuu/ZxH5lefvep2I/EdEQtra71lEXhGRPBFZV2teo3+vIjLVs/5WEZla37Ea0m4SgYjYgWeBC4B+wDUi0s+3UXlFDfBrY0xf4AzgF57vdQ8w3xjTE5jv+QzW9+/pmW4Bnm/9kL3mDmBjrc+PAX/1fOdDwDTP/GnAIWNMD+CvnvUC0dPAp8aYPsBgrO/eZn/PIpICTAcyjDEDADtwNW3v9/waMLnOvEb9XkUkFngAGAWMBB44kjxOizGmXUzAaGBerc/3Avf6Oq4W+J4fYA0StBno6JnXEdjsef8CcE2t9Y+uF0gT1tCm84FzgQ8BwWptGVT39w3MA0Z73v9/e+fzElUUxfHPgSkjo9KgsAxsINpmK6kWQWEgURs3ERTVP9AqiFbtI9xFkLSIKKgkxE2Lam0lRElFjRg5ZSlEBq2ETot7no7DlM4PfMx95wOPmXvuWdzv+z447/5gJmd5kraGKvVuBCbLxx2zz8AOYApoN99GgKMx+gx0AeO1+gqcBG6UxJfkLXdlZkbA4kOVULRYNNhUuBsYBbap6jSAfW61tFjuwwBwEfhj7S3ATw3/fw1LdS1otv45y28m8sAscMuWw26KSCsR+6yqX4CrwGdgmuDbGHH7nFCtr3X5naVCIBVi0ZydFZENwEPggqr++l9qhVhT3QcROQbMqOpYabhCqq6gr1nIAfuA66raDfxmcbmgEk2v2ZY2TgC7gO1AK2FppJyYfF6Of2msS3uWCkER2FnS7gS+pjSWhiIiawhF4I6qDln4u4h0WH8HMGPxGO7DAeC4iHwC7hGWhwaAzSKSs5xSXQuarX8T8GM1B9wAikBRVUet/YBQGGL2+QgwqaqzqjoPDAH7idvnhGp9rcvvLBWCF8BuO3GwlrDpNJzymOpGRAQYBN6p6rWSrmEgOTlwhrB3kMRP2+mDHmAumYI2C6p6SVU7VbWL4ONTVT0FPAP6La1cc3Iv+i2/qd4UVfUbMCUieyx0GHhLxD4TloR6RGS9PeeJ5mh9LqFaXx8DvSLSZjOpXoutjLQ3SVZ5Q6YP+ABMAJfTHk+DNB0kTAFfA6/s6iOsjT4BPtpnu+UL4fTUBPCGcCIjdR116D8EjNj3PPAcKAD3gRaLr7N2wfrzaY+7Rq17gZfm9SOgLXafgSvAe2AcuA20xOYzcJewBzJPeLM/X4uvwDnTXgDOVjMG/4kJx3GcjJOlpSHHcRynAl4IHMdxMo4XAsdxnIzjhcBxHCfjeCFwHMfJOF4IHMdxMo4XAsdxnIzzFxtipqempKxbAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for signal, value in simulation_results.items():\n",
    "    plt.plot(ts, value, label=signal)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Re-using parts of a diagram\n",
    "\n",
    "Let's compare the output of a PI and a PID controller on this system. We've already got the PI response, which we should store."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pi = simulation_results['y']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's swap out the PI controller for a PID."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gc_pid = blocksim.PID('Gc', 'e', 'u', 0.1, 50, 25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagram.blocks = [Gp, Gc_pid]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cab96b610b634938a7cedc4c1f9b4231",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=1000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "simulation_results = diagram.simulate(ts, progress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c1a113d68>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3wc1bn/8c+zq967bEtWsS33gm25ADbNBkwJkMQQm5KEAIYUyC8JySU3CXBzk3tzU0hCAgGHGgIYMAQMMdgJvdm4914lq/deVnt+f8xayLJsraVdbdHzfr302t3Z2Z1ntPZXs2fOnCPGGJRSSgU+m68LUEop5Rka6EopFSQ00JVSKkhooCulVJDQQFdKqSAR4qsNp6SkmJycHF9tXimlAtKGDRsqjDGpPT3ns0DPyclh/fr1vtq8UkoFJBE5cqrntMlFKaWChAa6UkoFCQ10pZQKEj5rQ1dKKU9qb2+nsLCQlpYWX5fiEREREWRmZhIaGur2azTQlVJBobCwkNjYWHJychARX5fTL8YYKisrKSwsJDc31+3XaZOLUiootLS0kJycHPBhDiAiJCcnn/G3DQ10pVTQCIYwP64v+6JNLgGkvqWd9/eWc6SyiagwO5Mz45k6PBGbLXj+ESul+k4DPQC0Ojp46N0D/PWDgzS3d5zw3IiUaO65bCyXTBjio+qUUsfZ7XYmTZqEw+Fg3LhxPP3000RFRRETE0NDQ4PXt69NLn6uqrGNrzy6hgff3se8cWksv+Nsdv18Aet+Mp8HrpuCzSYseWYD//mPbbQ5nL4uV6lBLTIyks2bN7N9+3bCwsJ45JFHBnT7eoTux+pa2lm8dA2HKhv5yw3TuGzS0M7nIsPsfGlaJldNGcZvV+/lkfcPUFLbwsM3TCMi1O7DqpVSAHPnzmXr1q0Duk0NdD/V4TR857lNHChv4KmbZzInL6XH9ULsNu65bCxZSVH85z+28b0XNvPQ9dO0XV0Nav/1+g52FtV59D3HD4vjvi9McGtdh8PBm2++yYIFCzxaQ2+0ycVPPfHRIT7YW87Pr554yjDv6vpZWfz0inG8ub2EB/61dwAqVEp119zczFlnnUV+fj5ZWVnccsstA7p9PUL3QwfKG/jt6j1cPD6dxTOHu/26W+bksr+sgT+/u5/p2YlcODbNi1Uq5b/cPZL2tONt6L6iR+h+xhjDz17dTkSonV9eM/GM+qKKCPdfNYGxQ2L54fItVDS0erFSpZS/0UD3Mx/sq+CTA5V8b34eaXERZ/z6iFA7f1w0lboWB/e+tt0LFSql/JUGuh9xOg2/enM3w5MiuX5Wdp/fZ8yQWO66aBQrt5Xw3p4yD1aolDqdU/U1H4g+6KCB7ldW7yxlV3EdP7h4DGEh/ftobjtvBCNSo7lvxQ5aul2MpJQKThrofuSxDw8yPCmSKycP7X3lXoSH2PnF1RM5UtnEYx8e9EB1Sil/12ugi8gTIlImIqdtkBWRGSLSISILPVfe4LHhSDXrj1Rzy7m5hNg983f2nFEpXDw+nUffP0hVY5tH3lMp5b/cSY6ngNP2jhcRO/B/wCoP1DQoPfHRIeIjQ7k23/1uiu740aVjaGxz8NC7+z36vkop/9NroBtjPgCqelntTuBlQM/A9UFFQyurd5awcHom0eGevTQgLz2WhdMzeebTIxRWN3n0vZVS/qXf3+1FJAP4ItDrKDQiskRE1ovI+vLy8v5uOmj8Y+Mx2jsMi2Z49uj8uO9dPBoRePDtfV55f6WUf/BEY+0fgP8wxvTalcIYs9QYk2+MyU9NTfXApgOfMYZl644yLSuBvPRYr2xjaHwki2YM55WNxzhW0+yVbSilrOFzzzrrLCZOnMi1115LU5P1rTgmJgaAw4cPExkZydSpUxk3bhwzZ87k6aef9tj2PRHo+cAyETkMLAQeFpFrPPC+g8KmghoOlDfyFS8dnR+35PyRACx9/4BXt6PUYObO8LkjR45k06ZN7Nq1i2XLlvH73/+eJ5980iPb73egG2NyjTE5xpgcYDnwLWPMq/2ubJB4fUsRYSE2Lp/U/66Kp5OREMmXpmWwbF0B5fU6JIBS3jZ37lz27z99Z4QRI0bwwAMP8OCDD3pkm72egROR54ELgBQRKQTuA0IBjDEDO3p7kHE6DSu3FXPB6FRiI0K9vr1vXjCK5RsKeeyjg/z4snFe355SPvPmPVCyzbPvOWQSXPYrt1Y9k+Fzp02bxu7du/tbHeBGoBtjFrv7ZsaYr/ermkFm3eEqSutaucIDFxK5IzclmismD+Pvnx7hWxeMIj7S+39ElBpMjg+fC9YRujvD5xpjPLZ9HT7Xh/65rZiIUBvzx6UP2DZvP28Er28p4qX1Bdw6d8SAbVepAeXmkbSn9WX43E2bNjFunGe+Meul/z7S4TSs3FbCRWPTPN73/HQmZsQzMzeJJz8+jKND5yBVypcOHz7M3XffzZ133umR99NA95HNBdVUNLRy2cSBaW7p6pY5uRyraWb1ztIB37ZSg92BAwc6uy1ed9113Hnnndx8880eeW9tcvGRt3eVEWITzhs98P3x549LJyspiic+OuT13jVKDSa9DZ+bk5NDc7P3rgXRI3QfeWd3GTNyknxyYtJuE75+Tg7rj1SzpaBmwLevlPIODXQfKKxuYndJPRf5cM7P62YMJzY8hCc+PuSzGpRSnqWB7gPv7rbGMLtonO8CPSY8hIX5mazcVqxzj6qg4ckugL7Wl33RQPeBt3eXkZMcxYiU6P6/mbMDyvfA/n/D9pdh1xtw5BNoru71pTfMyqa9w/DCuoL+16GUj0VERFBZWRkUoW6MobKykoiIM5tXWE+KDrCW9g4+PVDJ4plZiEjf3qS2EDY8DQffs66Gc5ziJEvKaBh7JUy+DtJO7uc6Ki2Gs0ck89zao9xx/kjstj7Wo5QfyMzMpLCwkGAZyTUiIoLMzMwzeo0G+gDbeLSaVoeTOaNSzvzFxsD6J2DVf0JHG2TkQ/43rEuSE3MgMsFa3lAOJVvg8Efw8R/howcg7xI470cwfMYJb3nj7Gy+/dxG3t9bxkVjB+4CJ6U8LTQ0lNzcXF+X4VMa6APsk/2V2G3CrBFJZ/ZCRyusvBs2/g1GzYcrfw8JWadeP28+zP2BFe4bn4Y1D8Pj82HK9XDxzyHG6i55yYR0UmPD+fuaoxroSgU4bUMfYB8fqGByZvyZDcbV2gDPLrTCfO7dcP2Lpw/zrmJS4by74btbYc73YNtL8PBs2PMWAKF2G4tnDOfdPWUUVOmMRkoFMg30AVTf0s7WwlrOHXkGzS3N1fDMNXD4Y/jiozDvZ2Czn/nGw2Ng/v1wx4cQOxSe/wq88X1oa2LRzCwEeO6zo2f+vkopv6GBPoDWHqyiw2k4Z1Syey9oa4S/L4SizXDd0zBlUf+LSBsHt70N59wJ6x+Hx+YxzFnMvHHpvLiugFZHrxNPKaX8lAb6APpofwXhITamZSX2vnJHO7z0dSjaCNc+BeO+4LlCQsLhkl/Aja9AXREsvZC7so9S2djGW9tLPLcdpdSA0kAfQJ8cqGBGThIRoW40mfzzB7BvNVzxOxh3pXcKGjUPlrwHcRlMfO8Wvh33EX9fc8Q721JKeZ0G+gCpbGhlb2kDZ490o7ll4zNWz5Q537O6JXpTUi7cshoZOY8ftj3MjIIn2VtS591tKqW8otdAF5EnRKRMRLaf4vkbRGSr6+cTEZni+TID3/oj1pWbs3J76a5Yss3qnph7Hlz0swGoDOuE6eLnaR2/kB+FvkjlK3dbfd6VUgHFnSP0p4DTTYx3CDjfGDMZ+G9gqQfqCjrrD1cRFmJjUmb8qVdqa7LazSMT4ctP9K03S1/ZQwlf+FfeS1zI2WUv0P7GDzTUlQowvQa6MeYDoOo0z39ijDk+cMga4MyuVR0k1h2uZkpmPOEhpwnpd38Jlfut7okxAz9OOjYbUV/4NY84riR0w+Ow8oca6koFEE+3od8CvHmqJ0VkiYisF5H1wTLegjua2zrYfqyW/JzTNLcUfAafPmS1mY84f+CK62ZGbhIvJ97GK5FfhnV/1VBXKoB4LNBF5EKsQP+PU61jjFlqjMk3xuSnpvrgCNRHNhfU4HAaZp4q0B1t8Nq3IT7Tuizfh0SEG2Zn8/3qL1E+6XYr1N+6R0NdqQDgkUAXkcnAY8DVxphKT7xnMFl/uAoRTt3//LOlULEXrngAwmMHtrgefHFaJhGhdh7gBpj9bVj7CLz7P74uSynVi34HuohkAa8ANxlj9va/pODz2eEqxqTHEh/Vw/gtjZXw/q9h1MUw+pKBL64H8ZGhXDVlGK9tKabu/Pth6k3wwa/hkz/5ujSl1Gm4023xeeBTYIyIFIrILSJyh4jc4VrlXiAZeFhENovIei/WG3AcHU42HqkmP+cUR+fv/Q+0NcClvxzYwnpxw6xsmto6eG1zEXzhjzD+Glj9U2scdqWUX+p1+FxjzOJenr8VuNVjFQWZPaX1NLZ1kJ/dQ/t5TQFseAqmfx1Sxwx0aac1OTOeiRlxPLv2KDfOzka+9FfrD8/r37X6rU/8sq9LVEp1o1eKetmWgloApmYlnPzk8SaMOd8bwIrcIyLcMCub3SX1bDxaDSFhcN0zkHU2vLIE9q72dYlKqW400L1sc0E1iVGhZCVFnfhEU5V1ef/kRZAw3DfF9eKqKcOICQ/h2TWuYXXDouD6ZZA+AV68yRrSVynlNzTQvWxLQS1ThiecPH/o9pfB0QKzlvimMDdEh4fwxakZvLGtmOrGNmthRLw1SmNCFjz3FSja5NsilVKdNNC9qKHVwd6yeqZk9tDcsuV5SJsAQyYPfGFn4PpZWbQ5nLy8sfDzhdEpcNOr1hAFz3wJynb7rkClVCcNdC/aVliLMXBW9/bz8r1wbAOctRi6H7n7mXFD45ienciza49iul5cFJ8BX30VbCHWjErVh31Wo1LKooHuRZsLagBOPkLf8hyIHSZd54OqztwNs7I4VNHIpwe6XTOWPNIK9fZm+Ns1UK+TYyjlSxroXrS5oJrs5CiSosM+X9jeApuftyaXiE33XXFn4PJJQ0mICuXZtT3MOZo+AW5YDg1l8MwXrZO9Simf0ED3oi0FtScfnW94EhpKrDk9A0REqJ2F0zJZtaOEsvqWk1cYPgMWP2eNFPnsQmitH/gilVIa6N5SUttCSV0LZw3vEuhtTfDhA5Az15rAIoAsnpWFw2l4aX1hzyuMuAAWPmlNaL3seuubiFJqQPV6pajqm872866Bvu6v0FgGX3nGR1X13cjUGM4Zmcxza49yx/kjsdt6OJk77kq45mH4x+2w/Btw3dNg72H8mkGusdXBpqM1bCms4VBFI4crGqlqbKOuxUGro4Mwu41Qu434yFDS4yMYGhfBqLQYxg+LY9zQuBOb8JTqQgPdSzYX1BBiEyYMi7MWtDbAx3+EkfMga7Zvi+ujG2Zl8+3nNvLB3nIuHJvW80pTFkFLHbz5Q2tI4GseAZt+ESyoauKt7SW8ub2YLYW1dDitHkPpceFkJ0czblgccRGhRITaaO9w0uZwUtPUTkldCzuLanlhfUHne41IiebskcmcMzKFc0Ymk6gBr1w00L1k+7FaxgyJJSLUNUPRhqegqRIuuMendfXHxePTSYkJ59m1R04d6GBdLNVaC+/8AkLC4co/DspQ73Aa3t5Vyt8+PcJH+ysAmDAsjm+eP5IZuUlMy0ogNsK9bzCVDa3sKq5ne1Etnx2q4tVNx3h27VHsNmFWbhILJg7h0glDSI+L8OYuKT+nge4Fxhh2FNVy6YQh1gJHqzVuS85cGD7Tt8X1Q1iIja/MyOQv7x2gsLqJzMSoU688926rO+OHvwNnB1z1p4GdI9WHnE7D61uL+MO/93GoopGh8RH84OLRXH1WBlnJp/mdnUZyTDhz8sKZk5fCHeePpL3DydbCWt7dXcZbO0q497Ud3PvaDmaPSGLh9OFcNnEI0eH633uw0U/cC4pqW6huamdChmtC6M3PWT1bvviIbwvzgBtmZfPI+wd5+pPD/OSK8adeUQQu+hnYw+C9/wWnA65+GOzB/U/uo30V/OKfO9ldUs/YIbE8dP00Lp2QTojds99QQu02pmcnMj07kbsvHcP+snpWbivhH5uOcfdLW7j3te1cPmkoX56WyazcJGw9nfNQQSe4/3f5yPZj1giLE4fFgdMJnzwIw6ZZPUEC3LCESK6YNJRlnxVw17y80zcZiFhNTLYQeOe/rVD/4qNBeaK0rL6FX7yxixVbishOjuLBxVO5ctLQAQvSUWmx3DUvljsvGsXGo9Us31DIG1uKWb6hkKykKL4yYzgLp2dqk0yQ00D3gh3HarHbhHFD4+DgO1B1EL78uN9f5u+uW+fmsmJLES+sK+DWuSN6f8F5d1sh/q97reanLz8GoZHeL3SArNhSxE//sY2WdiffnZfHNy8Y+fm5kwEmIkzPTmJ6dhL3XjmBVTtKeGFdAb9ZtYffrd7DhWPSuG7GcC4am0aoh781KN/TQPeC7UV1jEqNsf5Tr3scolNh3Bd8XZbHTM5MYGZOEk9+fJivn5PjXnPCud+FkEh480fWFaWLn7cG9wpgDa0O7l+xg+UbCpmWlcBvrp3CyNQYX5fVKTLMzjVTM7hmagZHKht5cX0BL60v5O3dZaRGh7LorCQWTk4mO85mne9wNFt/cI3pcvAh1n2xWX+EQyIgNMq6HxppNakFyYFKMOg10EXkCeBKoMwYM7GH5wX4I3A50AR83Riz0dOFBpLtx2qZk5cCNUdh71vWBBYh4b4uy6NumZvL7c9sYNWOUq6YPNS9F81aYo3U+I/b4YkFcOPLEJ/p3UK9ZFdxHd/8+waOVjVx10WjuGtensfbyfvEGGsYhpqjUHMEagug5ijZDWX8sKmKu+MrabeVEdJai22DEzb0c3tig7BYiIyHiATrj3Rkguv+8ceJ1kFNdKr1+UenQliM/iHwAneO0J8C/gz87RTPXwbkuX5mAX9x3Q5KZXUtlNW3MnFYvNVVEawp5oLM/HHp5CRH8dhHB90PdICJX7L+Qy+7Hh67GG54EYZM8l6hXvDW9hK+/+JmYsJDeP622cwakTzwRRhjhXXZLijbaQ1hXL7LGsnT0XziupGJEDsUopKRtLGE5cyBqGTqiWJTcSufHGmgsAEIjWT6yCGcPzqN3JRoBAO4Rth0Oq33bW+G9ibXret+awO01EBzDTRXQ32xdb+lBjraeq7fHn5iwJ90PxWik63bqBQI9W3bf31LO+X1rVQ1tlHR0EZVYxtVja1UNrZR3+Kgqc1BU1sHTa0dNLVb91vbnRhjcBpwum6tx4abz83lrnl5Hq/TnTlFPxCRnNOscjXwN2ONrbpGRBJEZKgxpthDNQaUHUV1AEwcGgOvPg+jLrYmgwgydpvwjTm53PvaDjYcqWJ6T3OmnkruXPjGW/D3hVaoX/1nmLTQe8V6iDGGP72znwf+tZcpwxNYetP0gTvJ2FxtDblcuAGOrbfuN3UZ/TJ2GKSNhfybIWmE9W8uIQvih1tzwPYgFjgPmGsMnx2q4oV1BfxqWzH/tdPJuKGhLJoxnGvOyiA+qo8nsY2xQr+5ChorXD/lXX4qoMm1rHy39c2io7Xn9wqNhqhkiEqEyCSISjrFrev58Dhrv938ZtzY6rCu2q1s5Fh1M8dqmimqaabQdb++xdHj62LCQ4iLCCEqPISoMDtRYXbSYiOIDLMTbrdhswl2EWw26/yGTUAQxg6J7dvvtBeeaEPPAAq6PC50LRuUgX68h8uk9i1QXwSX/tLHFXnPwumZPPCvvfzlvQM89rUzCHSwRmlc8h689DV4+RYo3gLz7/fbvuptDic/XL6F1zYX8aWpGfzPlyZ598RnUxUc/ggOfQCHP7QCDwCxJhQffRlkTIP0idbjyB4mUXGTiDBrRDKzRiRz31UTXCe8j3Lfih38cuUuLhidysXj05k3Lv3Mhh0QsaYtDItyr2nNGGsi8sZyaKw8Mfybqqw/DMdva45at801dH6L6Ikt1Ar2sFhMeAzt9igaiaTOGU5VeyjVrUJlC1S1Cm2E0mZCaCOUqNBwJkdGcU5UFPFZEcRHhRMfFUZMRCgxEWFEh4cQExFKWEiIq+nI1XxkOqzrLozTddth9e7qvixmKjDE/d+lmzwR6D01hPX4GxaRJcASgKys4DtqBdheVMuIlGgidy2D8HgYc5mvS/KaqLAQvnFuLg/8ay87imqZMCz+zN4gNh2+ugLeusfq2lm82erWGDfMOwX3UUOrg2/+fQMf7qvgh5eO4VsXjDx5SsH+am+xAvzgu1aIl2wDjHVkmn02TLoWMmfAsKkQEefZbXcRHxnKTbOzuWl2NtuP1bJ8QyFvbS9h9c5SbAIzcpK4aGwa54xMYfywuJ7H9OkrEQiPtX6S3Og9BVZAttR2CfxK2hoqqa6qpKamioa6GpoaamlvqqOjpo4IZzOxUk00LaRLKzk2B+HiIDSsnRBnO4Lz8/dudv1Unmrj/XDu/4Psczz+tnLCLDSnWslqcnnjFCdFHwXeM8Y873q8B7igtyaX/Px8s379+r7U7NfO/dU7zB4ezu8OX2s1I1z1oK9L8qra5nbm/Ood5uSl8Jcbp/f9jTY+Y/WACQm3rir1k15BlQ2t3PzUOnYU1fG/X5rEdfkenNC7vhT2rYK9q+DAu9DeaLUtD59pjcaZe551/UKIb8dqMcaw7Vgt/9pZyuodpewptYZHjosIYdaIZGbkJDIxI54Jw+KJjxyYawyMMVQ0tHGwvIED5Y0cKG/ovF9Y3YSzS6xlJEQyIjWakakxjEyLYWRqNKNSY0iNDT/5D3OHw2r2cbRa7f+OVuuI2hjrB9etcZ58H6yJa2x2163NugbjhGV260RySESfzwuIyAZjTH5Pz3niCH0F8B0RWYZ1MrR2sLafVze2caymmStyt1n/Oacs8nVJXhcfGcrXz83hT+/sZ29pPaPT+9g2OO0ma9Cyl2+FF26Es26ES/7bahv1kWM1zdz42FqKapp59MbpzB/fzwlJjLGOvPe+Zf0cc3Uxicu0/q2MuQxy5vhdH30RYXJmApMzE/jBJWMorWthzcFKPtlfyacHK/nXztLOdbOTo8hLiyU7OYrs5CiykqIYEh9BYlQYCVGhhIf03kzV5nBS3XT8xGMblY1tFNU0c6y6mcLqJo652rab2jo6XxMRaiM3JYbJmfF8cWpGZ4CPSI0mKuwMYs4eYv2ERZ/R78hf9HqELiLPAxcAKUApcB8QCmCMecTVbfHPwAKsbos3G2N6PfQOxiP0j/ZVcOPja9mU+xcSm4/CXZsHxaBU1Y1tzPm/d5g3Lp0HF0/t35s52qyhAj7+o9UufMkvrbAb4C5uBVVNLP7rGmqb23ny6zPIz+njH5b2Zjj0Iex90zoSrzsGCGRMhzELYPQCqx08gLvwVTa0sr2oju3Hatl+rJYD5Q0crWqipd150rrRYXYiw0IItQshdiHEZqPDaWhp76DV4ey87Ul8ZCgZCZFkJkaSkRjJ8MSoziPuYfGRg2Z4g34doRtjFvfyvAG+3cfagsqu4jriaCCh5BNrRqJBEOYAidFh3HR2Dks/OMBd8/IYldaPi2tCwmD+fVb3xje+B6/eAZuegfn/Zc2MNAAKqppYtHQN9S3t/P2WWSeOae+OumLYt9o6Cj/4ntW1LzQaRl4IF/4n5F0CMacZrTLAJMeEc/7oVM4fndq5zBhDeX0rR6qaOrv71TS1UdXYTnN7B44OJw6nob3DSYhNCA+xExFqIyLUTlRYCEkxYSRHh5EUbd0OiY9we2TKwUyvFPWgXSV1XBO9A+nogLFX+rqcAXXb3Fz+vuYIv1u9p39t6ccNmQTfWA0bn4J3/wcen2/9Ti/4MQw56VSOxxyttI7MG1odPHvrbCZlunGi1+m0TujuXWWFePFma3n8cDjrBusoPGeOz/tSDyQRIS0ugjQdO2ZAaaB70O7ieu4P2wj2dOsr9SCSHBPOkvNG8MC/9rLxaDXTsjxwWb/NBvnfgEnXwZqH4eMHYfcb1iBnZ3/HmizEg9+CjlQ2snjpGhrbOnj21llMzDhNmDdWwuEPYP/b1tF4Q6l1sitzJsy71wrxtPEB3ZSiAo8Guoe0dzg5WlbFlPANMPW6QdPc0tUtc3L526dH+NXK3bxw+2zPde0Lj4HzfwQzbrWuvv1sqTUZdVyG1Z1v8nX9Ds/DFY0s/usamts7eO62WSd3wWyuhoLPrC6FB9+H0m2u2uJh1DwrwEfNt65uVMpHNNA95GB5I/lmG+HOpkHX3HJcdHgI352fx89e3c67e8q4aGw/e4V0F5UEc79vHZ3vfh22LLMmDvn4D1ZPkVHzXF39plr9mN0M+EMVjSxa+iltDifP3Tqb8fHtcORTKNlq9UQ5tgEq91srH+9WeNFPIfd8a1tBOBywCkxu9UP3hmDr5fLa5mM0Lv8OiyLXYvuPQ0E3GJe72jucXPL7DxCBN787161uav3SUG41w+z/t3X03GoNvUBEPCTnQcJwqy07OsUaJTAsGhCrj3FHG1WVZfx73TYSTA1z0tqIqjtoHY0fF5MOGfnWVZmZM6ww97NuhWpw8XY/dAXsKqrlG/aN1tgtgzTMwZpJ594vjOfmJ9fx2IeH+PaFo7y7wZhUa/yS/Juho90aqKpoExRttsahL94Cu/95ykGikoDLiCIsYSjhkcMg8xpIybP+GKSPt5p1tB1cBQgNdA9pP/oZaVID467wdSk+d+GYNC6dkM6f3tnH1WcNO/3co55kD4WhU6yfrueknU7XCIFN0NYIGA7VOLj9uW00mAieWnJe3y+IUsqPDL4zd16SVfYuHdgh72Jfl+IX7v3CBAB+/vpOH1eCdYI6PMbq+52Uy35HGtc9f5Qq4nhaw1wFEQ10D6hqbONcx2cUJ04P+Fl4PCUjIZK75uWxemcp/9zqPyNB7C+rZ9HSNRgDz982mzwNcxVENNA94PCezYyyFdEyYoGvS/Ert80dweTMeH7y6jbK6lp8XQ77SutZtHQtAMuWzNIwV0FHA90DOna+AUDitGt8XIl/CbXbeOC6s2hu6+BHL2/FVz2qADYX1HDto58iAsuWzGZUmoa5Cj4a6B6Qeqy2vFgAABWTSURBVOzf7GIEyRkjfV2K3xmVFsOPLxvLe3vKeezDQz6p4eP9FVz/1zXERYTy8h3n9G+sGaX8mAZ6f9WXktW8kx1xc3xdid/62jk5XDZxCP/75i4+2lcxoNt+a3sxNz+5juGJUSy/42yykgeox41SPqCB3k8du1diw1CVqb1bTkVE+M21UxiZGsOdz2/kaGWT17dpjOGR9w/wzWc3MiEjjhdun60DRamgp4HeTy3bX6fAmUrSiH6OAx7kYsJDWPrVfJwGbnx8rVdPkrY5nPxo+VZ+9eZuLp80lOdvm01ClG9n/VFqIGig90drAxEFH7Lamc/Yod6b5zFY5KZE89TNM6hoaOXGx9dSXn+KGd77oaCqiWsf/ZSXNhRy17w8/rRoqncnc1bKj2ig98eBt7E723jb5JOXrifa3DE1K5HHvprP0aomFj7yCYcrGj323iu3FXP5gx9ysKyBh2+YxvcvHj1oZrFRCjTQ+2fX6zTYYqlMmub9QaiCyDmjUnj+ttnUNbfzxYc/5t9d5qTsi5LaFm5/Zj3fenYjI1JjWPnduVw+aaiHqlUqcLgV6CKyQET2iMh+Ebmnh+ezRORdEdkkIltF5HLPl+pn2ppgz5u8I7MZNfQMpyhTTM1K5JVvncvQ+Ehu/dt6fvzKNioazqwJpra5nd+t3sO8373He3vK+dGCMSy/42yGJ2lPFjU49To4l4jYgYeAi4FCYJ2IrDDGdB2k46fAi8aYv4jIeGAlkOOFev3HvtXQ1sCythnM0otU+iQ3JZp/fPscfrtqD098fJjXtxTx1bOzWTQj65TdC40x7CiqY/mGQl7eWEh9i4PLJw3hPxaMJTs5MGdqV8pT3BltcSaw3xhzEEBElgFXA10D3QDHzwrGA0WeLNKvHP7IGo51z5u0RQ9lbcs4vjpE28/7KjzEzk+uGM9XZmTx21V7eOT9Azz83gFGpcUwJTOBjIQIQu02apvbOVrVxMaj1VQ0tBFmt3HpxCHccf6Ik2cXUmqQcifQM4CCLo8LgVnd1rkfWC0idwLRwPye3khElgBLALKyss60Vv/wwW+smdyBrZPvp6PSrmOCeMCotBgeuWk6RTXNvL6liE8OVPLx/gpK61swBsJDbGQmRnLuqBTOHZXCxePSSYzWrohKdeVOoPfUTaD7oByLgaeMMb8TkbOBZ0RkojHGecKLjFkKLAVrxqK+FOxTxkDJNmsGm7nfZ/WBXMLsR8jWNluPGZYQye3nj+T2861hFDqcBmMMdpt4bo5SpYKUOydFC4HhXR5ncnKTyi3AiwDGmE+BCCDFEwX6lboiaKqEyV+BsVewt6yBEanRhNi1s5C32G1CiN2mYa6UG9xJonVAnojkikgYsAhY0W2do8A8ABEZhxXo5Z4s1C+UuGZ6HzoZgH2lDTo5glLKb/Qa6MYYB/AdYBWwC6s3yw4R+bmIXOVa7QfAbSKyBXge+Lrx5Vip3lKyFRBIn0B9SzvHapoZM0QDXSnlH9yaU9QYsxKrK2LXZfd2ub8TONezpfmh4i2QNALCY9l31JoZPk+HYlVK+Qlt/D0TJdu6NLfUA2iTi1LKb2igu6u5BmqOwJBJAOwtbSA8xKZXJSql/IYGuruOnxAdMgWAvaX15KXHYNfBn5RSfkID3V3derjsLa1ntF7yr5TyIxro7irZCjHpEJNGbXM7pXWteoWoUsqvaKC7q3grDOl+QlR7uCil/IcGujvaW6BiT5fmlgZAe7gopfyLBro7yneB09Glh0s9UWF2MhIifVyYUkp9TgPdHcVbrdshn58QzUuL0enNlFJ+RQPdHSXbICwWEnMBq8lFT4gqpfyNBro7irdYzS02G9WNbVQ0tOoJUaWU39FA701Hu9VlcdhUwGpuAT0hqpTyPxrovSnfDY4WyJgGaKArpfyXBnpvijZZt51H6A3EhocwND7Ch0UppdTJNNB7U7QJwuO7nBCtZ1R6jM6go5TyOxrovSnaBMOmgM2GMYZ9ZQ06hotSyi9poJ+OoxVKtnc2t1Q0tFHV2MZonaVIKeWH3Ap0EVkgIntEZL+I3HOKda4TkZ0iskNEnvNsmT5SugOc7TDMOiGqY7gopfxZr1PQiYgdeAi4GCgE1onICte0c8fXyQN+DJxrjKkWkTRvFTygup0Q3eMK9DHaw0Up5YfcOUKfCew3xhw0xrQBy4Cru61zG/CQMaYawBhT5tkyfaRoE0QmQUIWYPVwiY8MJTU23MeFKaXUydwJ9AygoMvjQteyrkYDo0XkYxFZIyILenojEVkiIutFZH15eXnfKh5IRZuto3NXj5Z9pfWMSY/VHi5KKb/kTqD3lF6m2+MQIA+4AFgMPCYiCSe9yJilxph8Y0x+amrqmdY6sNqboWxn5wVFxhj2uKadU0opf+ROoBcCw7s8zgSKeljnNWNMuzHmELAHK+ADV8l2MB2d7eelda3UtzgYoz1clFJ+yp1AXwfkiUiuiIQBi4AV3dZ5FbgQQERSsJpgDnqy0AFXtNG67XZCNE/7oCul/FSvgW6McQDfAVYBu4AXjTE7ROTnInKVa7VVQKWI7ATeBX5ojKn0VtEDomiTNYdo7FBAuywqpfxfr90WAYwxK4GV3Zbd2+W+Ab7v+gkOxzZY/c9dJ0D3lNSTEhNGcoz2cFFK+Se9UrQnTVVQsReGz+hctLesQUdYVEr5NQ30nhxztZ9nWoHudBr2ldZroCul/JoGek8KPwOxdV7yf6ymmaa2Du2yqJTyaxroPSlcB2kTINwK8L16yb9SKgBooHfndELhBsjM71y0t7QBQCeGVkr5NQ307ir2QmstDJ/ZuWhfaT1D4iKIjwz1YWFKKXV6GujdFX5m3WZ+3sNFL/lXSgUCDfTuCtdBRAIkjwKgw2nYX9ag7edKKb+ngd5dwTrr6Nx1QdHRqiZaHU7tsqiU8nsa6F211EL57hOaW473cNFp55RS/k4DvatjGwFz4hWiJccH5dI2dKWUf9NA76rgM0AgY3rnoj2l9WQkRBId7tawN0op5TMa6F0d+QjSJ0JEfOei3SX1jBsa58OilFLKPRroxznarBOiOed2Lmpp7+BgeQPjhmr7uVLK/2mgH1e0CRzNkP15oO8va8BpYOwQPUJXSvk/DfTjjnxk3XYJ9F3FdQCM1SN0pVQA0EA/7vDHkDoOopM7F+0uqSc8xEZOcrQPC1NKKfe4FegiskBE9ojIfhG55zTrLRQRIyL5p1rHL3U4oGAtZJ9zwuLdJXWMGRKL3SY+KkwppdzXa6CLiB14CLgMGA8sFpHxPawXC9wFrPV0kV5XsgXaGk44IWqMYVdxPeO0/VwpFSDcOUKfCew3xhw0xrQBy4Cre1jvv4FfAy0erG9gHP7Yus2e07movKGVqsY2bT9XSgUMdwI9Ayjo8rjQtayTiEwFhhtj3jjdG4nIEhFZLyLry8vLz7hYrznysTUYV2x656LdxdYVotrDRSkVKNwJ9J4akE3nkyI24PfAD3p7I2PMUmNMvjEmPzU11f0qvcnZAUc+7bH9HGCsjuGilAoQ7gR6ITC8y+NMoKjL41hgIvCeiBwGZgMrAubEaOkOa0KLLs0tYB2hD4mLIDE6zEeFKaXUmXEn0NcBeSKSKyJhwCJgxfEnjTG1xpgUY0yOMSYHWANcZYxZ75WKPe2Iq/28ywlRgF0l9dp+rpQKKL0GujHGAXwHWAXsAl40xuwQkZ+LyFXeLtDrDn8ECdkQn9m5qL3Dyf6yem0/V0oFFLeGEDTGrARWdlt27ynWvaD/ZQ0QpxOOfAJjLjth8cHyRto7jI7hopQKKIP7StHy3dBcdZoTonqErpQKHIM70I+3n2d3az8vrifMbmNEql7yr5QKHBrocRmQmHPC4t0ldYxKiyHUPrh/PUqpwDJ4E8sYOLoGsmZ3Tgh93O5i7eGilAo8gzfQawuhvhiGzz5hcWVDKyV1LTqGi1Iq4AzeQC9wjSE2fOYJi3cUWSdEJ2RooCulAssgDvTPIDTKmkO0i85AHxrf06uUUspvDd5AL90OQyaB/cSu+DuKaslMjCQ+KtRHhSmlVN8M3kAv3wOpY05avLOojgnDtLlFKRV4BmegN1VBUwWknBjoDa0ODlU2MmGYNrcopQLP4Az0ir3WbcroExbvKq7DGPQIXSkVkAZnoNcds24Thp+weMexWgA9QldKBaTBGegNrtmSotNOWLyjqI7k6DDS48J9UJRSSvXP4Az0xjIQO0QmnrB4R1Ed44fFIdLTJE1KKeXfBmegN5RCTBrYPt/9NoeTfWX12tyilApYgzTQyyH6xDlN95bW095h9ISoUipgDc5AbyyzjtC72Hn8ClENdKVUgHIr0EVkgYjsEZH9InJPD89/X0R2ishWEXlbRLI9X6oHNVaedIS+o6iW6DA7Ock6BrpSKjD1GugiYgceAi4DxgOLRWR8t9U2AfnGmMnAcuDXni7Uo9oaICzmhEU7iuoYNzQOm01PiCqlApM7R+gzgf3GmIPGmDZgGXB11xWMMe8aY5pcD9cAmfiz9iYIjex82OE07CzWS/6VUoHNnUDPAAq6PC50LTuVW4A3+1OUVzmd4GiBsM+bVg5VNNDU1sHEDO3hopQKXCG9r0JPbRCmxxVFbgTygfNP8fwSYAlAVlaWmyV6WLvri0SXI/QtBdYVolOGJ/iiIqWU8gh3jtALga7XyGcCRd1XEpH5wE+Aq4wxrT29kTFmqTEm3xiTn5qa2tMq3tfebN2GRnUu2lpYQ1SYnZGpMad4kVJK+T93An0dkCciuSISBiwCVnRdQUSmAo9ihXmZ58v0oM4j9M8DfUthLRMz4rHrCVGlVADrNdCNMQ7gO8AqYBfwojFmh4j8XESucq32GyAGeElENovIilO8ne91a3Jp73Cys7iOydp+rpQKcO60oWOMWQms7Lbs3i7353u4Lu85Huiuk6J7SuppcziZrO3nSqkAN/iuFG078Qh9a6HrhGimHqErpQLb4Av0bidFtx2rIT4ylKykqNO8SCml/N8gDPRG69YV6FsKapmcGa9D5iqlAt4gDPTjR+iRtLR3sKe0nsna3KKUCgKDMNA/Pym6s7iODqdhcqaeEFVKBb7BF+hdTopuLagBYIoGulIqCAy+QHe0WLchkWwtrCU1NlznEFVKBYXBF+gdbSA2sIew9VgtU/SEqFIqSAzOQLeFUt/SzoHyBm0/V0oFjUEY6A6wh7KloBZjYGqWBrpSKjgMwkBvA3soG45UIwJn6SX/SqkgMfgC3dkOtlA2Hq1mdFossRGhvq5IKaU8YvAFeocDYw9j09FqpmXr0blSKngMwkBvw4GduhYHU7MSfV2NUkp5zOALdGc7zR3Wbk/TQFdKBZHBF+gd7TR2CPGRoYxIie59faWUChCDMtAb2oRpWQnYdMo5pVQQGXSB3t7eRr1DtLlFKRV03Ap0EVkgIntEZL+I3NPD8+Ei8oLr+bUikuPpQj2lsamJdkKYlq2BrpQKLr0GuojYgYeAy4DxwGIRGd9ttVuAamPMKOD3wP95ulBPaWxuoQM7U/SCIqVUkHFnkuiZwH5jzEEAEVkGXA3s7LLO1cD9rvvLgT+LiBhjjAdrBWDrey8T98F9fX79EGcJxWGTiAl3a35spZQKGO6kWgZQ0OVxITDrVOsYYxwiUgskAxVdVxKRJcASgKysrD4VHBYdT1VUbp9eC1BFLjLxuj6/Ximl/JU7gd5TV5DuR97urIMxZimwFCA/P79PR+9jZ8yHGfP78lKllApq7pwULQSGd3mcCRSdah0RCQHigSpPFKiUUso97gT6OiBPRHJFJAxYBKzots4K4Guu+wuBd7zRfq6UUurUem1ycbWJfwdYBdiBJ4wxO0Tk58B6Y8wK4HHgGRHZj3VkvsibRSullDqZW109jDErgZXdlt3b5X4LcK1nS1NKKXUmBt2VokopFaw00JVSKkhooCulVJDQQFdKqSAhvupdKCLlwJE+vjyFblehDgK6z4OD7vPg0J99zjbGpPb0hM8CvT9EZL0xJt/XdQwk3efBQfd5cPDWPmuTi1JKBQkNdKWUChKBGuhLfV2AD+g+Dw66z4ODV/Y5INvQlVJKnSxQj9CVUkp1o4GulFJBIuACvbcJqwOViAwXkXdFZJeI7BCR77qWJ4nIv0Rkn+s20bVcRORB1+9hq4hM8+0e9I2I2EVkk4i84Xqc65pofJ9r4vEw1/KAmYi8NyKSICLLRWS36/M+O5g/ZxH5nuvf9HYReV5EIoLxcxaRJ0SkTES2d1l2xp+riHzNtf4+EflaT9s6lYAKdDcnrA5UDuAHxphxwGzg2659uwd42xiTB7ztegzW7yDP9bME+MvAl+wR3wV2dXn8f8DvXftbjTUBOQTQRORu+CPwljFmLDAFa/+D8nMWkQzgLiDfGDMRawjuRQTn5/wUsKDbsjP6XEUkCbgPa5rPmcB9x/8IuMUYEzA/wNnAqi6Pfwz82Nd1eWlfXwMuBvYAQ13LhgJ7XPcfBRZ3Wb9zvUD5wZr96m3gIuANrKkMK4CQ7p831nj8Z7vuh7jWE1/vQx/2OQ441L32YP2c+Xy+4STX5/YGcGmwfs5ADrC9r58rsBh4tMvyE9br7SegjtDpecLqDB/V4jWur5lTgbVAujGmGMB1m+ZaLRh+F38AfgQ4XY+TgRpjjMP1uOs+nTAROXB8IvJAMwIoB550NTU9JiLRBOnnbIw5BvwWOAoUY31uGwj+z/m4M/1c+/V5B1qguzUZdSATkRjgZeD/GWPqTrdqD8sC5nchIlcCZcaYDV0X97CqceO5QBICTAP+YoyZCjTy+dfwngT0fruaC64GcoFhQDRWc0N3wfY59+ZU+9mv/Q+0QHdnwuqAJSKhWGH+rDHmFdfiUhEZ6np+KFDmWh7ov4tzgatE5DCwDKvZ5Q9AgmuicThxn4JlIvJCoNAYs9b1eDlWwAfr5zwfOGSMKTfGtAOvAOcQ/J/zcWf6ufbr8w60QHdnwuqAJCKCNTfrLmPMA12e6joB99ew2taPL/+q62z5bKD2+Fe7QGCM+bExJtMYk4P1Ob5jjLkBeBdronE4eX8DfiJyY0wJUCAiY1yL5gE7CdLPGaupZbaIRLn+jR/f36D+nLs40891FXCJiCS6vt1c4lrmHl+fROjDSYfLgb3AAeAnvq7Hg/s1B+ur1VZgs+vncqz2w7eBfa7bJNf6gtXj5wCwDasXgc/3o4/7fgHwhuv+COAzYD/wEhDuWh7herzf9fwIX9fdj/09C1jv+qxfBRKD+XMG/gvYDWwHngHCg/FzBp7HOk/QjnWkfUtfPlfgG6793w/cfCY16KX/SikVJAKtyUUppdQpaKArpVSQ0EBXSqkgoYGulFJBQgNdKaWChAa6UkoFCQ10pZQKEv8fdvr7HDuTB3gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ts, y_pi, label='PI')\n",
    "plt.plot(ts, simulation_results['y'], label='PID')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that adding the derivative action has improved control."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Disturbances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can simulate a more complicated block diagram with a disturbance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](../assets/standard_feedback.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "Km = blocksim.LTI('Km', 'ysp', 'ytildesp', 1, 1)\n",
    "Gc = blocksim.PI('Gp', 'e', 'p', Kc=8, tau_i=10)\n",
    "Gv = blocksim.LTI('Gv', 'p', 'u', 1, 1)\n",
    "Gp = blocksim.LTI('Gp', 'u', 'yu', [1], [10, 1])\n",
    "Gd = blocksim.LTI('Gd', 'd', 'yd', [1], [10, 1])\n",
    "Gm = blocksim.LTI('Gm', 'y', 'ym', [1], [1, 1])\n",
    "\n",
    "blocks = [Km, Gc, Gv, Gp, Gd, Gm]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "sums = {'e': ('+ytildesp', '-ym'),\n",
    "        'y': ('+yd', '+yu')}\n",
    "\n",
    "inputs = {'ysp': blocksim.step(),\n",
    "          'd': blocksim.step(starttime=50)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagram = blocksim.Diagram(blocks, sums, inputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = numpy.arange(start=0, stop=100, step=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e825d893e847437e93607261224bcbec",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=2000), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "results = diagram.simulate(ts, progress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1c19feccf8>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5RcdZnu8e9b1bd0Ot1JX0jMje6MISRECKQFFYmJcAQUwbhcIgLiIRBdjspB1zk6JwdF0RkZvC8ZNMNgwDMQCYoDHgzMQjCMEkggEHMDYhJIE5J0mlw76a7be/6o6qT6lq50qrqyd57PWrW6a+9dVW+tXXny63f/am9zd0REJPgixS5ARETyQ4EuIhISCnQRkZBQoIuIhIQCXUQkJEqK9cL19fXe2NhYrJcXEQmkF154YZe7N/S1rmiB3tjYyMqVK4v18iIigWRmr/e3Ti0XEZGQUKCLiISEAl1EJCQU6CIiIaFAFxEJCQW6iEhIKNBFREIilIG+ee9mHnr1ITqTncUuRURkyBTti0WFEk/F+fx/fp5t7dvYun8rN8+8udgliYgMidCN0F/a+RLb2rdhGL957TfEk/FilyQiMiRCF+jPvfUcUYty2/m3sbdzL6t3rS52SSIiQyJ0gf7q7ldprG5k9oTZGMbz258vdkkiIkMidIG+ee9mJo2cRE15DVNqp7Bqx6pilyQiMiRCFejxZJyt+7fSWN0IwNTaqWx4ewO6ELaInAxCFehbD2wl6UmaapoAOL32dHZ37mbHwR1FrkxEpPBCFejb27cDMLZqLABT66YCsL5tfdFqEhEZKgMGupndY2Y7zWxNP+uvNrPVmdtfzOys/JeZm9aDrQCcMuwUAKaMmgKkD5SKiIRdLiP0RcAlR1m/GfiAu58J3AYszENdg7Lz4E4AGirTV2eqLK1kzPAxbN63uVgliYgMmQG/Keruy8ys8Sjr/5J1dzkw/vjLGpydB3dSXVZNRUnF4WVN1U1s3qtAF5Hwy3cPfR7wh/5Wmtl8M1tpZitbW1vz/NLpQD+l8pRuyyaNnMTmvZs100VEQi9vgW5mc0gH+tf628bdF7p7s7s3NzT0edHq47Lr0C4ahnV/3qbqJg4lDmmmi4iEXl4C3czOBO4GrnD3tnw852Ds7tzNqIpR3ZZNGjkJgE17NxWjJBGRIXPcgW5mE4HfAte6e1Gnk+zp3ENNeU23ZV1z0tVHF5GwG/CgqJk9AMwG6s2sBfgmUArg7j8HvgHUAf9iZgAJd28uVMH9SaaS7I/tZ2T5yG7L6yrqGFE6QoEuIqGXyyyXqwZYfwNwQ94qGqR9sX0AvUboZkZTTRNb9m4pQlUiIkMnNN8U3du5F+gd6ACNNY2aiy4ioReaQN/TuQegV8sF0n30nQd30h5vH+qyRESGTGgC/fAIvayPEXrm7Itb9m0ZwopERIZWeAI9lg70/kbooJkucnz2xfZx91/v5gcrf8C6tnXFLkekl9BcJHpPR7rlUlPRe4Q+YcQEIhbRgVEZtLZDbVz7h2vZun8rJZESFq1dxFdnfpXPTv9ssUsTOSxUI/SIRagqreq1rixaxviq8Wq5yKB97/nvsfPgThZdsohlVy7j4saL+cELP+CBDQ8UuzSRw0IzQt8f209VaRUR6/v/qMaaRrVcZFDW7FrD0i1L+cJZX2Dm6JkA3H7B7XQmO7n9+ds5tfpU3jf2fUWuUiREI/T2eHufo/MujdWNvLHvDVKeGsKqJAzuX38/w0uH85kzPnN4WTQS5fYLbqeppomvLfsa2w5sK2KFImmhCvTK0sp+1zfVNNGR7Dh8VSORXLTH23l8y+NcNukyhpcO77ausrSSH8/5MYlUgpufvpnOZGeRqhRJC1WgDzRCB810kWPzzJvPEEvFuLTp0j7Xn1p9Kv90wT+xrm0d313+XZ2mWYoqVIHecwSVrbGmEdBcdDk2T77+JLUVtcxomNHvNrMnzGb+mfN5eOPDPPTaQ0NYnUh3oQr0o7Vc6irqGFGmk3RJ7pKpJH9+88/MnjCbaCR61G2/cNYXOH/c+fzjc//I6tbVQ1ShSHehCfQD8QNHbbmYGU3VOkmX5G7D2xvYH9/PeWPOG3DbroOkoytHc/PTN+tYjRRFaAL9YPzgUVsuoJN0ybFZsX0FAO8e8+6ctq8pr+HHc35Me7ydax67hlfefqWQ5Yn0EopAd/cBe+igk3TJsVmxYwWN1Y00VOZ+ucTTa0/n3kvuxd256v9dxS9e/gUH4wcLWKXIEaEI9EOJQzg+8Ai96yRdarvIAFKe4sUdL+Y8Os82pXYKD370QWZPmM3PXvoZH/rNh7jt2dt4euvTtB7M/8XRRbqE4puiXSPugQJ98qjJALyy+xXOqD+j4HVJcG3Zt4UD8QOc2XDmoB5fN6yOH87+IS+3vsyv1v2KRzc9yoOvPghAbUUtTTVN/F3N3zFp5CSaapqYVDOJ0ZWjyVz1S2RQQhHoB+IHgIEDfcKICQwvHc76tvUweSgqk6Bau2stANPrph/X85zVcBZnfeAsOhIdrG1by/q29by25zU27dnE0i1LD19pC9Kf36bqJiaNnMTkkZM5o/4MptVNG/BzLdIlFIHe1aMc6IMfsQhTRk1hw9sbhqIsCbA1u9YwrGTY4VMvH6+Kkgpmjp55+FwwkD7209bRxqY9m9i098jt2W3P8sjfHgHAMCbVTOKM+jOYXj+d6XXTmVI7hbJoWV7qknAJRaDn2nIBmFo3ld++9luSqeSAc4vl5LWmbQ1Ta6cW9DNiZtQPq6d+WD3nvuPcbuvaDrWxtm0ta3etZU3bGv7rzf86HPIlkRJOG3UaZ9SdwdS6qTRWNzJxxEROqTxFLZuTXCgCPdeWC8DU2qkcShzi9f2vM6lmUqFLkwCKp+K88vYrXDnlyqLVUDesjlnjZzFr/CwgPZrfcXAHa3atSd/a1rB081KWvLrk8GOGlQxj4oiJTKyeyLiqcYwZPoYxlWMYU5X+WVtRq8APuVAEetcI/WhfLOpyeu3pAKxvW69Alz79bc/f6Ex2Mr3++Prn+WRm6YAePoaLTr0ISM/E2d6+ndf3vd7t9uruV/nT1j8RS8W6PUdZpIzRw0czZvgY6ivqqRtWl75VdP9ZW1Grlk5ADRjoZnYPcBmw0917fcIt/V/+T4APAweBz7r7i/ku9Gi6Av1oX/3vMmnkJIaVDGN162o+MukjhS5NAmjNrjXA8R8QLbSIRRhbNZaxVWN579j3dlvn7uzu3M1b7W+xvX0729u3s6N9B2+1v8WOgztY27aWto62fr+TUVVaRXVZNdXl1emfPX4fUTaC6rJqqsqqqCyppLK0stvPYSXD1NIsglxG6IuAnwH39bP+UtJzRiYD5wF3ZX4OmYOJ3A6KApRGSpnRMIMVO1Yc8+ts3b+VJ7Y8wbq2dezp3EPUojRUNjBhxAQaaxppqm6isaaR8mj5MT+3nDjW7FpDdVk140eML3Ypg2Zm1FbUUltRyxl1/U/R7Uh00NbRRtuhzC3z++7O3ezr3Me+WPq2ee/mw7/neprg8mj54ZAfVjKMypJKykvKKYuWUR4ppzxaTmm0lPJo+veyaFl6XbScskhZt2WlkVJKIiXpm6V/RiPRbvdLIiVELXpku573M9uFue00YKC7+zIzazzKJlcA93n6vKHLzWykmb3D3d/KU43drG7Zw78vf6PbsldjWwD45sOv5rSz2mLjeS3+LDcv+TNlNmLA7WN+gA2x+2lJLAOcShtNuY3ESfCib6DTd2dtbVRaA8MjY6m0Uyi1SqJUELVSHMc9iZMiRRLPuqUOL09klqXwAl+M49g/2IX8h3Bsz22FqsVhb/xZJiTLWPHTqwvzGieo4ZnbxD7XVmZuY4iTot2SHLQkHZaiw1J0kqIzc7/z8P3MOmunk/10WIpDlv6Ex82JkyJunrmfIo6TsKE5/XDE058hI/3tSsPSP7st777egIinf3Zfn1nmWc+T9Rh63nd4b+37+OrHf5j395WPHvo4YGvW/ZbMsl6BbmbzgfkAEyf2/bEZSOv+Tv70avdv28Wq98DwMpa9tiun50iWvQMa4Jmtz1HScdZRt02VbKOj7l/x6D5KDsymtH0WlhxFV3eyBIhajFS0FS/dQapkB50lO+ko3UFrdANYZ3pP98cj4FHSH5soeAQjmlnW9dHJVSH/MRzrcxf6H2Zhnt9I0cBB5u5po/GATrA11FJAwiBmEMOIR6DTIIGRsCPrk1nLkqTvJ81I0LWu97IkkDAjBXjmuTzzmik78rsDqa7tstYlM487/JjMdt0f1/35yCw/cjPcILGvMJ+tfAR6X4nT5782d18ILARobm4e1L/IC6eO5sKpo7st+87yZ3liy3CW/e8Lc3qOeGoWs3+9iNnntvHd9/f/mE17N/Hfl36bEZFSfjLn/w7qIJm7cyhxiFgyRjQSJWpRopEoEYtQYuH+8y+Q2tvgjklw6R1w3vxiVyNyTPIR6C3AhKz744EhvcBiR6KDipKKnLcvjZQyZ8Ic/vjGH4kn45RGS3tt07K/hRufuBHDuOfiezi1+tRB1WZm6QNFORywFRE5Hvk4OdcjwGcs7T3A3kL1z/vTkexgWMmwY3rMhxo/xP74fpa9uazXup0Hd3LjEzfSkehg4YcWDjrMRUSG0oCBbmYPAM8CU8ysxczmmdnnzezzmU0eAzYBG4F/Bb5QsGr7cShx6JhG6ADvG/s+xg4fy6I1i7pdB7LtUBvzn5jP2x1v8/OLfs5po07Ld7lyQtM1QSW4cpnlctUA6x34+7xVNAgdiQ4qoscW6CWREq6ffj3fee47/G7j75g7eS4t+1v40h+/xLb2bdx54Z28q+FdBapYTng6tiEBFIpvinYkOhhRNvD0w54+cdoneOL1J7j12Vt5bPNjrG5dTcQi3HnhnYM6D7aISDGF4wIXyWNvuUD6OpA//eBP+eRpn6Sto42LTr2IJR9dojAXkUAKzQh9MIEO6W+XLnjPgjxXJIHl6qFLcIVihN6ROPZZLiIiYROKQD+UOHTMB0VFRMImFIGuEbqISAgCPZ6Kk/DEoHvoIt2phy7BFfhA70h0AKjlIvmleegSQKEJ9GGlarmIyMkt8IF+KHEI0AhdRCQ0ga6DopIXmocuARb4QO9IZnroOigqeaUeugRP8ANdB0VFRIAQBboOiorIyS7wgX4omemhRxXokg/qoUtwBT7Qu0bo5SXlRa5EQkXz0CWAAh/osWQMgPKoAl1ETm6BD/TOZCegQJc80bRFCTAFuohISIQm0EsjpUWuRMJFPXQJnlAEenm0HNNBLBE5yQU+0GPJGGXRsmKXIaGhHroEV06BbmaXmNkrZrbRzL7ex/qJZvaUma0ys9Vm9uH8l9q3rhG6iMjJbsBAN7MocCdwKTANuMrMpvXY7P8AD7r72cCngH/Jd6H9iSVjCnTJP7XwJIByGaGfC2x0903uHgMWA1f02MaB6szvNcC2/JV4dJ3JTrVcRETILdDHAVuz7rdklmW7FbjGzFqAx4Av9fVEZjbfzFaa2crW1tZBlNubWi6SV5qHLgGWS6D39bdnz0/9VcAidx8PfBj4lZn1em53X+juze7e3NDQcOzV9kEHRUVE0nIJ9BZgQtb98fRuqcwDHgRw92eBCqA+HwUORCN0KQz10CV4cgn0FcBkM2syszLSBz0f6bHNG8CFAGY2lXSg56enMgCN0EVE0gYMdHdPAF8EHgfWk57NstbMvm1ml2c2+ypwo5m9DDwAfNZ9aJqRnclOyiMaoUu+qIcuwVWSy0bu/hjpg53Zy76R9fs64Pz8lpYbTVsUEUkL/DdFNW1RCkLz0CWAQhHoGqGLiIQg0HVQVPJK89AlwAIf6Bqhi4ikBTrQU54inoor0KUA1EOX4Al0oHddT1QtFxGRgAe6Lj8n+aceugRXOAK9RIEuIhKOQNcIXfJN89AlgAId6Oqhi4gcEehAPzxC17lcJF80D10CLNCB3jVCV8tFRCTggd41QlfLRfJPPXQJnlAEukboIiIBD3QdFJX8Uw9dgivQga4RuojIEYEOdB0UlYLRPHQJoEAHug6KiogcEYpA1whd8kbz0CXAAh3oarmIiBwR6EDvTHZiGCWRnK51LXIM1EOX4Al0oMeSMcqj5ZgOYImI5BboZnaJmb1iZhvN7Ov9bPNJM1tnZmvN7P78ltm3zmSnDohKnqmHLsE1YK/CzKLAncB/A1qAFWb2iLuvy9pmMvAPwPnuvtvMTilUwdkU6CIiR+QyQj8X2Ojum9w9BiwGruixzY3Ane6+G8Ddd+a3zL7FU3HKIgp0KQC18SSAcgn0ccDWrPstmWXZTgNOM7M/m9lyM7ukrycys/lmttLMVra2tg6u4izxZJzSaOlxP4+ISBjkEuh9DVV6NhpLgMnAbOAq4G4zG9nrQe4L3b3Z3ZsbGhqOtdZeYqkYpREFuuSR5qFLgOUS6C3AhKz744FtfWzzH+4ed/fNwCukA76g4qm4Al1EJCOXQF8BTDazJjMrAz4FPNJjm98BcwDMrJ50C2ZTPgvtSywZ00FRKRD10CV4Bgx0d08AXwQeB9YDD7r7WjP7tpldntnscaDNzNYBTwH/093bClV0F43QRUSOyOkrlu7+GPBYj2XfyPrdga9kbkMmnoxTVVY1lC8pInLCCvQ3RTVCFxE5ItCBrh66FIzmoUsABTrQ46m4TswlIpIR6ECPpWL6pqjkl+ahS4AFOtD1TVERkSMCHegaoUvhqIcuwRPoQE+kEprlIiKSEehA1ywXyT/10CW4AhvoyVSSpCfVQxcRyQhsoMdTcQC1XKQwNA9dAijwga6DoiIiaYEN9FgyBqCWi+SX5qFLgAU20DVCFxHpLriBnsz00DVCFxEBghzoGqGLiHQT2ECPpTI9dM1ykbxSD12CK7CBrpaLiEh3gQ10jdCloDQPXQIosIF+uIeur/6LiAABDvTD89A1Qpd80jx0CbDABnpXD10jdBGRtOAGus7lIgWlHroET06BbmaXmNkrZrbRzL5+lO0+YWZuZs35K7FvXQdFNQ9dRCRtwEA3syhwJ3ApMA24ysym9bHdCODLwHP5LrIvmrYohaEeugRXLiP0c4GN7r7J3WPAYuCKPra7DfhnoCOP9fVL0xZFRLrLJdDHAVuz7rdklh1mZmcDE9z993ms7ag0QpeC0jx0CaBcAr2vT/bhv0vNLAL8CPjqgE9kNt/MVprZytbW1tyr7IN66CIi3eUS6C3AhKz744FtWfdHANOBp81sC/Ae4JG+Doy6+0J3b3b35oaGhsFXjWa5SIFoHroEWC6BvgKYbGZNZlYGfAp4pGulu+9193p3b3T3RmA5cLm7ryxIxRnxZJyoRYlGooV8GRGRwBgw0N09AXwReBxYDzzo7mvN7NtmdnmhC+xPPBXX6FwKSD10CZ6SXDZy98eAx3os+0Y/284+/rIGFkvGdEBURCRLoL8pqhG65J966BJcgQ30WDKm87iIiGQJbKBrhC4FpXnoEkCBDnTNQRcROSK4gZ6M66Co5J/moUuABTbQY6mYRugiIlkCG+jxlEboUkjqoUvwBDbQY8mYDoqKiGQJbKDHU3FNW5QCUA9dgiuwga4RuohId4EN9EQqoYOiUjiahy4BFNhA17lcRES6C2yg65uiUhCahy4BFthAj6XUQxcRyZbT6XNPRPGkZrlIIamHHjTxeJyWlhY6OobkOvUFV1FRwfjx4yktzX3gGthA1whdRLK1tLQwYsQIGhsbsYAf1HZ32traaGlpoampKefHBbLl4u7pWS4aoUveqYceVB0dHdTV1QU+zAHMjLq6umP+ayOQga4LRItIX8IQ5l0G814CGeixZAxAI3QpnBAFg5w8AhnoXSP0kkhgDwGIiORdIANdI3QpGM1DlwAL5BBXPXQROZpvPbqWddv25fU5p42t5psfPaPf9bfccgv19fXcdNNNACxYsICysjKefPJJ9u3bRyKR4K677uKCCy6gqqqKz33uczz11FOMGjWKxYsX09DQcNw1BnOEnsqM0HUuFykY9dDl2MybN497770XgFQqxeLFi6mqquLiiy/mpZde4uWXX2bGjBkAtLe3c8455/Diiy/ygQ98gG9961t5qSGnEbqZXQL8BIgCd7v793qs/wpwA5AAWoHr3f31vFTYh3gyM0LXuVxEpA9HG0kXSmNjI3V1daxatYodO3Zw9tln8+53v5vrr7+eeDzOxz72scOBHolEuPLKKwG45ppr+PjHP56XGgYcoZtZFLgTuBSYBlxlZtN6bLYKaHb3M4GHgH/OS3X96Gq5aIQu+aceugzeDTfcwKJFi/jlL3/J9ddfz6xZs1i2bBnjxo3j2muv5b777uvzcfmabplLy+VcYKO7b3L3GLAYuCJ7A3d/yt0PZu4uB8bnpbp+qIcuIieiuXPnsnTpUlasWMHFF1/M66+/zimnnMKNN97IvHnzePHFF4F0S+ahhx4C4P777+f9739/Xl4/l5bLOGBr1v0W4LyjbD8P+ENfK8xsPjAfYOLEiTmW2FvXLBe1XKRgNA9dBqGsrIw5c+YwcuRIotEoTz/9NHfccQelpaVUVVUdHqEPHz6ctWvXMnPmTGpqavj1r3+dl9fPJdD7+mT3+XepmV0DNAMf6Gu9uy8EFgI0NzcP+m9bjdBF5ESUSqVYvnw5S5YsAeC6667juuuu63Pb2267jdtuuy2vr59Ly6UFmJB1fzywredGZnYRsAC43N0781Ne3zQPXQpG89BlkNatW8c73/lOLrzwQiZPnlyUGnIZoa8AJptZE/Am8Cng09kbmNnZwC+AS9x9Z96r7EEjdBE50UybNo1NmzbltO2BAwcKUsOAI3R3TwBfBB4H1gMPuvtaM/u2mV2e2ewOoApYYmYvmdkjBak2QyN0KTz10CV4cpqH7u6PAY/1WPaNrN8vynNdR5VIJQCN0EVEsgXzm6IaoUvBqIcuwRXIQFcPXUSkt0AGete5XBToUjCahy7H6dZbb+X73//+kL5mIAO961wuarmIiBwRyNPnxlIxDCNq0WKXImGjFno4/OHrsP2v+X3OMe+CS7931E2++93vct999zFhwgQaGhqYOXNmfmsYQCADPZ6KUxYtC9X1A0Uk2F544QUWL17MqlWrSCQSnHPOOQr0XMSTcfXPpcA0WAi0AUbShfDMM88wd+5cKisrAbj88ssHeET+BbOHnhmhi4icSIrdNQhkoMeSMV0gWgpETXQZnFmzZvHwww9z6NAh9u/fz6OPPjrkNQQyFeOpuC5uISInlHPOOYcrr7ySGTNmcOqpp3LBBRcMeQ2BDPRYMqZzoUth6YC7DMKCBQtYsGBB0V4/kC0XjdBFRHoLZKDHUjHNcpHC0PnQJcACGeiJZEKzXKTA1HKR4AlkoGuELiLSWzADXQdFRUR6CWSgx1P6pqgUinroElyBDPRYMqYeuhSWWugSQIEMdI3QRUR6C+QXi+JJnctFRPp3+/O3s+HtDXl9ztNrT+dr536t3/W33HIL9fX13HTTTUD6S0ajR4/my1/+cl7rOBqN0EWyaR66DNK8efO49957AUilUixevJirr756SGsI5Ahd0xal8NRED7KjjaQLpbGxkbq6OlatWsWOHTs4++yzqaurG9Iacgp0M7sE+AkQBe529+/1WF8O3AfMBNqAK919S35LPSKejGvaooiccG644QYWLVrE9u3buf7664f89QdsuZhZFLgTuBSYBlxlZtN6bDYP2O3u7wR+BNye70K7uDuxVEznchGRE87cuXNZunQpK1as4OKLLx7y189lhH4usNHdNwGY2WLgCmBd1jZXALdmfn8I+JmZmXv+G5LPrPgpAKUr7oG//CrfTy8nu/jBYlcgAVZWVsacOXMYOXIk0ejQX/M4l0AfB2zNut8CnNffNu6eMLO9QB2wK3sjM5sPzAeYOHHioAqurmzg0uhILqweA5Fhg3oOkaNqvADGnVPsKiSAUqkUy5cvZ8mSJUV5/VwCva+jQz1H3rlsg7svBBYCNDc3D2r0PmP6p5kx/dODeaiISMGsW7eOyy67jLlz5zJ58uSi1JBLoLcAE7Lujwe29bNNi5mVADXA23mpUEQkAKZNm8amTZuKWkMu89BXAJPNrMnMyoBPAY/02OYR4LrM758A/liI/rmIyNGEKXYG814GDHR3TwBfBB4H1gMPuvtaM/u2mV2e2ezfgDoz2wh8Bfj6MVciInIcKioqaGtrC0WouzttbW1UVFQc0+OsWG++ubnZV65cWZTXFpHwicfjtLS00NHRUexS8qKiooLx48dTWtr9Ozdm9oK7N/f1mEB+U1REpKfS0lKampqKXUZRBfJcLiIi0psCXUQkJBToIiIhUbSDombWCrw+yIfX0+NbqCcBveeTg97zyeF43vOp7t7Q14qiBfrxMLOV/R3lDSu955OD3vPJoVDvWS0XEZGQUKCLiIREUAN9YbELKAK955OD3vPJoSDvOZA9dBER6S2oI3QREelBgS4iEhKBC3Qzu8TMXjGzjWYWyrM6mtkEM3vKzNab2VozuymzvNbM/tPMXsv8HFXsWvPJzKJmtsrMfp+532Rmz2Xe768zp28OFTMbaWYPmdmGzP5+b5j3s5ndnPlMrzGzB8ysIoz72czuMbOdZrYma1mf+9XSfprJtNVmNujLZQUq0HO8YHUYJICvuvtU4D3A32fe59eBJ919MvAk4TtN8U2kT9Hc5XbgR5n3u5v0xcjD5ifAUnc/HTiL9PsP5X42s3HAl4Fmd58ORElfXyGM+3kRcEmPZf3t10uByZnbfOCuwb5ooAKdrAtWu3sM6Lpgdai4+1vu/mLm9/2k/5GPI/1e781sdi/wseJUmH9mNh74CHB35r4BHyR90XEI2fsFMLNqYBbp6wng7jF330OI9zPpM7wOy1zZrBJ4ixDuZ3dfRu+rtvW3X68A7vO05cBIM3vHYF43aIHe1wWrxxWpliFhZo3A2cBzwGh3fwvSoQ+cUrzK8u7HwP8CUpn7dcCezAVWIJz7ehLQCvwy02q628yGE9L97O5vAt8H3iAd5HuBFwj/fu7S337NW64FLdBzuhh1WJhZFfAb4H+4+75i11MoZnYZsNPdX8he3MemYdvXJcA5wMKaxPkAAAGHSURBVF3ufjbQTkjaK33J9IyvAJqAscBw0u2GnsK2nweSt8960AI9lwtWh4KZlZIO8393999mFu/o+lMs83NnserLs/OBy81sC+k22gdJj9hHZv40h3Du6xagxd2fy9x/iHTAh3U/XwRsdvdWd48DvwXeR/j3c5f+9mveci1ogZ7LBasDL9M//jdgvbv/MGtV9sW4rwP+Y6hrKwR3/wd3H+/ujaT36R/d/WrgKdIXHYcQvd8u7r4d2GpmUzKLLgTWEdL9TLrV8h4zq8x8xrveb6j3c5b+9usjwGcys13eA+ztas0cM3cP1A34MPAq8DdgQbHrKdB7fD/pP7lWAy9lbh8m3Vd+Engt87O22LUW4L3PBn6f+X0S8DywEVgClBe7vgK83xnAysy+/h0wKsz7GfgWsAFYA/wKKA/jfgYeIH2cIE56BD6vv/1KuuVyZybT/kp6FtCgXldf/RcRCYmgtVxERKQfCnQRkZBQoIuIhIQCXUQkJBToIiIhoUAXEQkJBbqISEj8f2A67/q43XRoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for name in ('ysp', 'd', 'y'):\n",
    "    plt.plot(ts, results[name], label=name)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Algebraic equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes it is useful to be able to handle non-linear calculations in block diagrams. This deviates from the strict interpretation of block diagrams but can be useful for instance in calculating the response of a controller with output limits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def limit(t, u):\n",
    "    return numpy.clip(u, 0, 0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "Gp = blocksim.LTI('Gp', 'ulimited', 'y', 10, [100, 1], 50)\n",
    "Gc = blocksim.PI('Gc', 'e', 'u', 0.1, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "limiter = blocksim.AlgebraicEquation('Limiter', 'u', 'ulimited', limit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagram = blocksim.Diagram([Gp, Gc, limiter], \n",
    "                           sums={'e': ('+ysp', '-y')}, \n",
    "                           inputs={'ysp': blocksim.step()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = numpy.arange(start=0, stop=1000, step=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_results = diagram.simulate(ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LTI: ulimited →[ Gp ]→ y\n",
       "PI: e →[ Gc ]→ u\n",
       "AlgebraicEquation: u →[ Limiter ]→ ulimited"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diagram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1c1a16ad68>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXzcVb34/9c7k33f0zR723RfUpouQCnIZkGliCBwVRa54vWKcEW9inrhXtT7latewZ9chQuCKFCQK1CxWNkKIm1pum+0SdNsTdLsW7PPnN8fM4EhJM1MMpNPMvN+Ph7zyMz5nM/nvD+dZt6Zz/mcc8QYg1JKqeATYnUASimlrKEJQCmlgpQmAKWUClKaAJRSKkhpAlBKqSAVanUA3khNTTX5+flWh6GUUtPKrl27mowxacPLp1UCyM/Pp6SkxOowlFJqWhGRypHKPboEJCLrReSoiJSJyHdG2H6niBwWkf0i8pqI5LnKi0Rkm4gccm271m2fx0XkhIjsdT2KxntySimlvDdmAhARG/AgcBmwELheRBYOq7YHKDbGLAWeA/7LVd4N3GCMWQSsB+4XkUS3/b5ljClyPfZO8FyUUkp5wZNvAKuAMmNMuTGmH9gIbHCvYIx5wxjT7Xq5Hch2lR8zxpS6ntcCDcBHrkMppZSafJ4kgCyg2u11jatsNLcALw8vFJFVQDhw3K34R65LQz8XkYiRDiYit4pIiYiUNDY2ehCuUkopT3iSAGSEshEnEBKRzwPFwE+GlWcCvwNuNsY4XMV3AfOBlUAy8O2RjmmMedgYU2yMKU5L0y8PSinlK54kgBogx+11NlA7vJKIXAx8D7jCGNPnVh4P/Bn4vjFm+1C5MabOOPUBj+G81KSUUmqSeJIAdgKFIlIgIuHAdcAm9woishx4COeHf4NbeTjwPPCEMeYPw/bJdP0U4Erg4ERORCmllHfGHAdgjBkUkduALYAN+I0x5pCI3AuUGGM24bzkEwv8wfl5TpUx5grgs8A6IEVEbnId8ibXHT9PikgazktMe4F/8u2pTa5dlS1sL28hOSacyxdnkhAdZnVISil1RjKd1gMoLi42U20gWO+AnW/+YR8v7a97vyw1NoJfXF/EObNTLYxMKaWcRGSXMaZ4eLnOBTQBdofhq0/u5s8H6rjzkrkc+PdLeeGr55IYHcYXH9/JnqpWq0NUSqlRaQKYgIfeOs5r7zXwH1cs4vaLComLDKMoJ5Gnv7SG1NgI7ti4l66+QavDVEqpEWkCGKeq5m4eeLWUyxbP4Atr8j60LS0ugp9fW0RVSze/fL3MogiVUurMNAGM0/2vHkME7vnUIlwd3x+yMj+ZTy/P4jd/P0FtW48FESql1JlpAhiHE02neWHvSb6wJo8ZCZGj1rvzkrnYHYbH36mYvOCUUspDmgDG4eG3jhNmC+FL62adsV5OcjSXLZ7B0zuqtC9AKTXlaALwUmfvAC/urWVD0UzS40b/63/ILWsL6Owb5LmS6jHrKqXUZNIE4KVN+2rp7rdz/apcj+ovz01iSVYCG3dWM53GXCilAp8mAC89u7Oa+TPiKMpJHLuyy7Urc3ivvpMDJ9v9GJlSSnlHE4AXqlu62VfTzqeXZ414589oriiaSWRYCBt36mUgpdTUoQnAC385WA/AZYszvdovPjKMyxdn8qe9tfT02/0RmlJKeU0TgBc2H6xjcVY8uSnRXu/72ZU5dPYNsuVQvR8iU0op72kC8FBtWw97qtq8/ut/yKr8ZGYmRPKnfR9ZSkEppSyhCcBDrx05BcD6xTPGtX9IiPDJZTN5q7SRtu5+X4amlFLjognAQ2+VNpGTHMWs1JhxH+NTS2cyYDd6GUgpNSVoAvDAgN3BtuPNnFeY5tXdP8MtzoonPyWaP+2rG7uyUkr5mSYAD+ypaqOrb5B1hRNblF5E+NSymbxzvInGzr6xd1BKKT/SBOCBv5U2YgsRzpmTMuFjfWrZTBwGXj6o3wKUUtbSBOCBv5U2sTwnkfjIia/zOzcjjllpMbxy+JQPIlNKqfHTBDAGYwzHTnWyJDvBZ8e8ZGEG2443094z4LNjKqWUtzxKACKyXkSOikiZiHxnhO13ishhEdkvIq+JSJ7bthtFpNT1uNGtfIWIHHAd8xcykd5VP2o+3U93v53cZO8Hf43m0oUZDDoMW482+OyYSinlrTETgIjYgAeBy4CFwPUisnBYtT1AsTFmKfAc8F+ufZOBe4DVwCrgHhFJcu3zK+BWoND1WD/hs/GD6pZuAHKSfJcAinKSSI0N18tASilLhXpQZxVQZowpBxCRjcAG4PBQBWPMG271twOfdz3/OPCKMabFte8rwHoR2QrEG2O2ucqfAK4EXp7Q2fhBdUs3uXKKWSH10NTpk2PagGsL+thy9Ah9g8uICLX55LhKKeUNTxJAFuA+jWUNzr/oR3MLH3yQj7RvlutRM0L5lJNw8Le8FXEfbPTtcb8FfEtg7/ZEitZ+wrcHV0opD3iSAEa6Nj/iyiYi8nmgGDh/jH29OeatOC8VkZvr2SIsvhTVcpgOYoi/6gGfHrffbsfxwm10730eNAEopSzgSQKoAXLcXmcDH5nRTEQuBr4HnG+M6XPb94Jh+251lWePdUwAY8zDwMMAxcXFk76kVtTpGupDs4hfeo1PjxsO7P/ro+Q2v+3T4yqllKc8uQtoJ1AoIgUiEg5cB2xyryAiy4GHgCuMMe63tmwBLhWRJFfn76XAFmNMHdApImtcd//cALzog/PxuaT+OjojZ/rl2D35F5Ft6qgpO+CX4yul1JmMmQCMMYPAbTg/zI8AzxpjDonIvSJyhavaT4BY4A8isldENrn2bQF+gDOJ7ATuHeoQBr4CPAKUAceZgh3A9sFB0h2NDMT759JT1qorAajdOSVzn1IqwHlyCQhjzGZg87Cyu92eX3yGfX8D/GaE8hJgsceRWqCx9gQzxE5IUt7Ylcchu2A+JySHmMrXgO/7pQ2llBqNjgQ+g+aaUgCiMmb5rY261HOY03OA3u4uv7WhlFIj0QRwBqdPHQcgaWah39qIXnAxETJAWcmrfmtDKaVGogngDOwtFTiMkJ7jvwRQuPLj9BsbnYf/6rc2lFJqJJoAziC0o4oGSSE8ItJvbcTEJVAWsYjUhm1+a0MppUaiCeAMYrtraAkf3yLw3ujKWkuho5yG+pqxKyullI9oAjiDlIF6Tkf7f4aK5CXOefCqSqbcnbBKqQCmCWAUvT2nSacFu5/GALgrWHouHcTgOP7G2JWVUspHNAGMoqHaeQtoaEq+39uyhYZSGn0WuW3vgpn02S6UUkFKE8Ao2k6WARAzY/aktNeXu44ZppG68oOT0p5SSmkCGEVPYzkAKVlzJ6W9GUUfB+Dk3lcmpT2llNIEMArTUkG/CSU10z/TQAxXMHcJTSQilX+flPaUUkoTwCjCu6o5FZJOiG1yVuuSkBAq4paT07Eb43BMSptKqeCmCWAUcT21tEX4fwyAO0fOOaTTQvXxQ5ParlIqOGkCGEWavZ6emOyxK/pQ5rJLAKjdp/MCKaX8TxPACDrbW0ikC0fC5Fz/H5JduIwWErQfQCk1KTQBjKCh6hgAYan+mwZ6JBISQmVcEbmdezA6HkAp5WeaAEbQUeccAxCfOTljANwNZJ9DJk2crDg66W0rpYKLJoAR9DU61wFIy5mcMQDu0pZcBECtjgdQSvmZJoARSGslXSaKhOT0SW87b95ZtBJHiPYDKKX8TBPACCJP19AQmoGETP4/T4jNRnn0Mma27570tpVSwUUTwAgS+mrpiJhpWfu9M1cz05yipb7SshiUUoHPowQgIutF5KiIlInId0bYvk5EdovIoIhc7Vb+MRHZ6/boFZErXdseF5ETbtuKfHda42ccDtLtDfTG5lgWQ/L8tQBU7n3dshiUUoFvzAQgIjbgQeAyYCFwvYgsHFatCrgJeMq90BjzhjGmyBhTBFwIdAPui99+a2i7MWbv+E/Dd5obThItfUjS5I4BcDdryTn0mHD6Tmy3LAalVODz5BvAKqDMGFNujOkHNgIb3CsYYyqMMfuBM01iczXwsjGme9zRToKmGucYgMi0yR0D4C4iIpITEfNIbtZ+AKWU/3iSALKAarfXNa4yb10HPD2s7Ecisl9Efi4iESPtJCK3ikiJiJQ0NjaOo1nvdNU7bwFNmDnH722dSUfqcgoGjnO6q9PSOJRSgcuTBCAjlHk1TFVEMoElwBa34ruA+cBKIBn49kj7GmMeNsYUG2OK09LSvGl2XAabKgBIzy30e1tnEjPnXMLETvn+ty2NQykVuDxJADWAe49oNlDrZTufBZ43xgwMFRhj6oxTH/AYzktNlgtpr6SFeKJjEy2NI7/oAgDaj/7N0jiUUoHLkwSwEygUkQIRCcd5KWeTl+1cz7DLP65vBYiIAFcCU2ItxOjuGppCZ1gdBnHJM6gOySKmYZfVoSilAtSYCcAYMwjchvPyzRHgWWPMIRG5V0SuABCRlSJSA1wDPCQi709oLyL5OL9BvDns0E+KyAHgAJAK/HDipzNxSf11dEZN7jTQo6lPKCK/5xAOuy4Qo5TyvVBPKhljNgObh5Xd7fZ8J85LQyPtW8EIncbGmAu9CXQy2AcHSXc0URVn3RiAD8leRVLrn6ks20/evCkxTEIpFUB0JLCbhpPHCRM7tuR8q0MBIG3ROgAaDr9lcSRKqUCkCcBNs2sMQFT65E8DPZLcwmW0EwPV71odilIqAGkCcNNzqhyAlGxrbwEdEmKzUR65mIy2KTFIWikVYDQBuBlsqcBuhLTsqfENAOB0+lnkOqrpamuyOhSlVIDRBOAmrKOKRkklLHzEQcmWiJ29GoDKAzogTCnlW5oA3MT21NISbv0YAHd5S84DoPP4DosjUUoFGk0AblIG6jgdPTXGAAxJSk6lUrKJaNB+AKWUb2kCcOnt7iKNVuzxuVaH8hEN8YvI6T6MceiAMKWU72gCcDlVXQZAaGq+tYGMwMxcQSpt1FaVWR2KUiqAaAJwaastBSA2Y+rcATQkZd7ZAJw8pB3BSinf0QTg0tvgGgOQM8/iSD4qb8FK+k0oA5UlVoeilAogmgBcTGslfSaMlIwpMg+Qm9CIKCrCZ5PQut/qUJRSAUQTgEt4ZzWnbOmE2GxWhzKi1sSlFPQfY3BgYOzKSinlAU0ALvG9J2kLn2l1GKMKzV1BjPRReXSP1aEopQKEJgCXdHs9PTHjWep4cqTPPxeA5qPvWByJUipQaAIA2lubiOc0JjHP6lBGlTVrER3EwEldIUwp5RuaAIDGqqMARKQWWBzJ6EJsNioj5pHSPiVWzlRKBQBNAEBnvXOAVVzmHIsjObOOlGXkDVbQ19NpdShKqQCgCQDoazwBQPoUHAPgLjJvJaHioOrQdqtDUUoFAE0AgLRV0kkU8clpVodyRpmLnB3BbaU6M6hSauI8SgAisl5EjopImYh8Z4Tt60Rkt4gMisjVw7bZRWSv67HJrbxARHaISKmIPCMi4RM/nfGJPF1Dg20GiFgVgkcys/JoIJmQ+n1Wh6KUCgBjJgARsQEPApcBC4HrRWThsGpVwE3AUyMcoscYU+R6XOFWfh/wc2NMIdAK3DKO+H0isa+OjsipOwZgiIhwMmouqZ1HrA5FKRUAPPkGsAooM8aUG2P6gY3ABvcKxpgKY8x+wKP5ikVEgAuB51xFvwWu9DhqHzIOB+n2U/TFTr0pIEbSm7qEHHsN3V3tVoeilJrmPEkAWUC12+saV5mnIkWkRES2i8jQh3wK0GaMGRznMX2muaGGKOlHkvKtaN5r0XkrCBFD5SHtB1BKTYwnCWCkC+PGizZyjTHFwD8A94vIbG+OKSK3uhJISWNjoxfNeqap+hgAkelTdwyAu6xFzqmh246/a3EkSqnpzpMEUAO4Xx/JBmo9bcAYU+v6WQ5sBZYDTUCiiISOdUxjzMPGmGJjTHFamu/v0umqPw5AYmahz4/tD6mZ+TSRSEi9zgyqlJoYTxLATqDQdddOOHAdsGmMfQAQkSQRiXA9TwXOBQ4bYwzwBjB0x9CNwIveBu8LA82uMQC5c61oflzqouaR1vWe1WEopaa5MROA6zr9bcAW4AjwrDHmkIjcKyJXAIjIShGpAa4BHhKRQ67dFwAlIrIP5wf+j40xh13bvg3cKSJlOPsEHvXliXnK1l5FE4lExcRZ0fy49KQtIc9eRUendgQrpcYvdOwqYIzZDGweVna32/OdOC/jDN/vHWDJKMcsx3mHkaViTtfQHDqDVKsD8UJ0/gpsVY9QcWgHS9dcanU4SqlpKuhHAicN1NEZNXWngR5J1sJzAOgo1yUilVLjF9QJYHCgn3RHE4Px02MMwJCkjDxaicemHcFKqQkI6gTQUFNOqDiwJedbHYp3RDgZNZ+0Lh0RrJQav6BOAC0nnWMAotNnWRyJ93rTlpBvr6K9o8PqUJRS01RQJ4DuU+UAJGVPn1tAh0TlrSBUHFQe3ml1KEqpaSqoE4C9pQK7EdKzpt83gOwFawBoL9cEoJQan6BOAGGd1ZwKSSM0PMLqULyWkDmLNuKw6dTQSqlxCuoEENt9ktawTKvDGB8RaqPmaUewUmrcgjoBpA7WcTp6eo0BcNeTOtQR3GV1KEqpaShoE0BvdxeptGFPyLU6lHGLyl9BmNipOKIzgyqlvBe0CeBUlfMW0LCU6TEN9EiyFjinhm4/rh3BSinvBW0CaKstAyB2xmyLIxm/hMzZtBOrHcFKqXEJ2gTQ2+BcByB1Go4BeJ8ItdHzSNM1gpVS4xC0CcC0VtJrwkiZMb3mARquJ3UJ+Y5K7QhWSnktaBNARFc1p2wZSMj0/ieIyismXOycOKwdwUop70zvT78JiO+tpT18mo4BcJOlI4KVUuMUtAkgzV5PT+z0vvwDEJ85hw7tCFZKjUNQJoD2lkbi6cYk5lkdysSJcDJ6PunaEayU8lJQJoDG6qMARKTmWxuIjwx1BLd1dFodilJqGgnKBNBR5xwDEJdZaHEkvhGdt8LVEaz9AEopzwVlAuhvOgFAWu48iyPxjZlDI4K1I1gp5QWPEoCIrBeRoyJSJiLfGWH7OhHZLSKDInK1W3mRiGwTkUMisl9ErnXb9riInBCRva5HkW9OyYPzaauigxgSklInq0m/is+crR3BSimvhY5VQURswIPAJUANsFNENhljDrtVqwJuAr45bPdu4AZjTKmIzAR2icgWY0yba/u3jDHPTfQkvBXZVU2DbQbxk92wv4hwMnqedgQrpbziyTeAVUCZMabcGNMPbAQ2uFcwxlQYY/YDjmHlx4wxpa7ntUADkOaTyCcgsb+OjsiZVofhU72pSyhwVNLarh3BSinPeJIAsoBqt9c1rjKviMgqIBw47lb8I9eloZ+LyIjLconIrSJSIiIljY2N3jb7EcbhIMN+iv7Y7AkfayqJfn9EsPYDKKU840kCkBHKjDeNiEgm8DvgZmPM0LeEu4D5wEogGfj2SPsaYx42xhQbY4rT0ib+5aG5vppIGUCS8yd8rKkkc6GzI7hDO4KVUh7yJAHUAO5DZrOBWk8bEJF44M/A940x24fKjTF1xqkPeAznpSa/a3KNAYhMm34LwZ9J/AztCFZKeceTBLATKBSRAhEJB64DNnlycFf954EnjDF/GLYt0/VTgCuBg94EPl5dp5xXoBJnzpmM5ibP0NTQukawUspDYyYAY8wgcBuwBTgCPGuMOSQi94rIFQAislJEaoBrgIdE5JBr988C64CbRrjd80kROQAcAFKBH/r0zEYx0FwBQEbuNF4HYBQ9qUuY5aikRTuClVIeGPM2UABjzGZg87Cyu92e78R5aWj4fr8Hfj/KMS/0KlIfsbVX0UQiqdGxVjTvV9H5KwivepwTR0pIXvMxq8NRSk1xHiWAQBLTXUNTWCaBMQTsw2YuPBvego7jO0ETgM+c7myj8uA2Oir2YLoaCOnvxBEWjUSnEJW1iKwFa0id5gsLqeAUdAkgub+Ok/HLrA7DL+IynFNDh9bvtTqUac8+OMi+V39PyP6NLDy9k4Uy6Cw3wmmJIsr0ESZ2KAPehDLbbBpnfoy8i25lZn5gTDGiAl9QJYCB/j7STRMV8blWh+If2hE8YcbhYPdfHiN95085y9RyihR2Z3yGqPkXkzV/NckZ2cTbbBiHg/a2Zk4e20X70b+RULOVVVWPIo89yr7olUR87F+Zv+oSq09HqTMKqgTQUFNOlhhCA2wMgLuetKUsqniC5rZ2UhITrA5nWqmvKqXhqS+zoncX5SH57Fp5P0WXfoGM0I/+mkhICAnJaSSsWQ9r1r+//4lXfs286mdJ3nw1+15fSewnfsjsJWsm+1SU8khQzQbaevIYAFEZgTUGwF10/irCxU7loe1jV1bv2/f6RiJ/cwGzew6xfd63yfvuLlZcfjO2ET78RzMjt5Czb/kZkd88yLZZt5PXe4T859az/Vdfpquj1Y/RKzU+QZUAuhvKAUjOCrxbQIfMXLwWgK7yHRZHMn1sf+oHLHvryzTb0mj9wmusuf67Xn3wDxcdm8DZN/wAuX0vJalXsObURrr/+yz2b/0/H0at1MQFVQKwt1QwaEJIzyqwOhS/iUvLpUFSiKjfY3UoU55xONj28O2sOfZTdsesI+sbb5M9Z7HPjp+QnMbqrz3Be5/4P7pDYli69Yts/58v0dtz2mdtKDURQZUAwjqqaQhJIzQs3OpQ/Ko2ZhFZpw+PXTGIGYeDHQ/fxtm1v2VHygaWff15Iv00NmT+youZ8c3t7Ei7mjUNz1L3k7OpPLLLL20p5Y2gSgBxPSdpCc+0Ogy/68tYTjb1NDd4PGVT0Nn++3tYU/8kO1KvYtVXH5/QJR9PREbHsvqrj7Lv/EeId7SRtvEydr/8mF/bVGosQZUAUgbr6Y72eibraSd2tvOuk5pDb1scydS084Vfcnb5LyiJu4iVX3kECZm8X4NlH7uGwX98g+qwfM7a8S9se+irDA70T1r7SrkLmgTQ3dVOKm3YEwJ0DICb3MVnYzdC74l3rQ5lyjm2eyvL9tzDwYgilt72FCE226THkJE9m/xvbmVHygbOrvs9R356Ka2NdZMeh1JBkwAaqksBCEsJ3A7gIXHxSVTYcolu1BHB7prqq0jcdDNNIclkf+kZwiMiLYslIjKa1V97gp1L72Vu70F6/ud8ThzSO7fU5AqaBNBWWwZA7IzZFkcyORrjF5PTcwTjcIxdOQg47HbqH/sCcaaLnqt+R2LqDKtDAmDlVXdQ8alnCTP9ZDz7KXb/5XGrQ1JBJGgSQK9rDEBqTuCOAXBnslaQSBf1FTotBMC7T/+AxX17ObD0e1NuZO684gvh1q3OfoHtd7D9kTtx2O1Wh6WCQNAkANoq6THhpKQH1lrAo0mddy4AtYf+ZnEk1ju+/x3OKv0Fe2LWsvLTt1sdzojSZuaT9403eDfxctbUPMq+n32SzvYWq8NSAS5oEkBEZzWnbBmTeseHlQoWrKDbRDBQVWJ1KJbq7TmN7YVbaZd48m+a3Dt+vBUZFcPK259k+7xvs+T0dpofWEdN2aQslKeC1NT9bfCx+L5a2iNmWh3GpAkNC6MiYi5JrfutDsVSe566m3xHNbXn/4SktKk/BkRCQlhz/Xd575LfkuhoJf73l3LgzT9aHZYKUEGRAIzDQfpgPb0xwXH5Z0hnyjIKBo7T29tjdSiWqDyyixVVj1ESfzHLPnaN1eF4ZfHaKzh9w6s0h6Sx8PUvsv3Je7VDX/lcUCSAjtZG4qQHk5RndSiTKqpgNeEyyPH926wOZdI57HZ6/ngb3RJFwecesDqcccmatYD0r7/Jvti1rCn9GSUPXKfzCCmfCor1ABqqj5EAhKcG7jTQI8leej68A23H3oZVlizBbJmS5x9g1cBh3l32Q1ZlTN9vfjFxiRTd+SLbnvguZ1f+mmM/u4DEm58N6AkNfam/r5f6yvdor6+kp7kKe1sNIV31hPZ3EDrQSbj9NJH2bkKNczS2YAADCP0SQX9IFAO2SAZtUQyGxWKPTMHEpGGLSyciIYPopBkkZeSRnJFtyaDCifIoAYjIeuABwAY8Yoz58bDt64D7gaXAdcaY59y23Qh83/Xyh8aY37rKVwCPA1E4F5y/wxhjJnQ2o+isc44BSMgMjjEAQ5Jn5FEn6UTW7bQ6lEnV0dZM4cGfcyRsESs3fNXqcCYsxGbj7JvvY89flzD379+g538v4L1PPMr8lRdbHdqU0lhbQc2BN+mr3kd4aykpPeVk2WvJlQ9fOmslji6Jozckmr7QWNrCknDYwnF+/AuIgDHY7D2E2nsIt/cQO9BKZPdpEtvaiZKPTt3Rb0JpCEmlLXwGPVGZDMZlY0vOJTotn+SsuaRnz5qSk1COmQBExAY8CFwC1AA7RWSTMcZ9uskq4Cbgm8P2TQbuAYpxptVdrn1bgV8BtwLbcSaA9cDLEz2hkfQ3nQAgLTf41mo9GbeMvI4SjMMxpe+A8aVDG7/PatNJ8yfuC6hzXn7p5zmRPY/wP3yOWS9dy87qu1l51R1Wh2UJ43BQcWQnp/ZuIbxuJ1ldh8igmTSc6zbXhmTSFFVAbdLFhKbPIzotn8QZeaRk5pMUFUPSBNru7mqnrbGOzuZaelrr6WupxtFWRXhXLTE9deS17SC19S+EVH/w9+ygCaE2JJXWsExOR2dhT8glLKWA2BmzSc2eS8qMHEv+r3ryDWAVUGaMKQcQkY3ABuD9BGCMqXBtG95L9XHgFWNMi2v7K8B6EdkKxBtjtrnKnwCuxE8JQNoqaSeGhMQUfxx+SrNnrSSt4xVqK48xs2C+1eH4XVXpflbUPcPu5MsoLjrP6nB8rmDhStpv+xtHH76OlfvvZkfdfpb/44OWTmsxWZpqKznx7p+Q8jfI7yihgDYKgFrJoDquiBOZZ5E09xzyFq0mJyqGHD/FER2bQHRsApzh96m/r5fGkydorS2l+1Q59pYKwjqrie0+yay2d0ht2wyVH9TvNWGcsmXQFp5Jb2wOJjGX8NQC4jPnkJYzj4TkNL+ciycJIAuodntdA6z28Pgj7ZvletSMUO4XUadraLTNIBhXyE1ZsA6O/Cd1B98IigTQ/Mdvk0IYBZ+9z+pQ/CYhJfS1YJsAABQfSURBVIMF39zC9kduZ82ppyn7rz3YPvMQBQtXWh2az9VXl1Hxt6dJOPEy8/oPkyqGZhI4Eb+S8oLzyS2+nJk5c5hqN3iHR0SSNWsBWbMWjLi953QnDVXHaKsrpbfhBKa1koiuauJ7ayloOkR8UzeUfVC/g2harv0T+QuKfRqnJwlARijz9Fr9aPt6fEwRuRXnpSJyc8c3k2fClT+ht6ttXPtOd/kLiuk0UTgqdwBfsTocvzr49iaW97zDjtm3sTozsGd9DQ0LZ81Xfs2ev64l7527iH1mPdtn/zMr/+Eev69t4G+1J96j6u9Pk1TxMvMGjzIDKA/JZ0f+l0kv/jQFC1eSMg07XN1FxcSRt2AFeQtWjLi9vaWRxuqjdNSV0d90AmmrYv6MfJ/H4cn/lBr40LepbMDTlUZqgAuG7bvVVZ49rHzEYxpjHgYeBiguLh5XJ3He/LPGs1tACA0LozxyAamtgb1EpHE4iNz6H9SRRtFnv2t1OJNm+aWfp6XoYxz87ZdZU/4Ljv34L3D5T5l71vlWh+aV6tJ91LzzDGlVf2GO/TgzgVLbHLYV3Eb2udcxa84SgukevoTkNOdln2Vr/dqOJwlgJ1AoIgXASeA64B88PP4W4D9FZKjP5VLgLmNMi4h0isgaYAdwA/D/eRe68tTp9GKWVP0vXe3NxCYEZj/I3ld+x/LBMnYs+xGZkTFWhzOpktOzSPrGJko2P0J+yY9IfnED7779SQqvn9qjnyuP7KJ22zNk1GxhlqOCHOBo6Dy2z/kXcs+9nsKC+RRaHWSAE0/uvBSRy3He5mkDfmOM+ZGI3AuUGGM2ichK4HkgCegF6o0xi1z7fhEY+pPsR8aYx1zlxXxwG+jLwNfGug20uLjYlJQE99w243Hgby+y5LUb2H/Boyy94Gqrw/E5++AgNf9ZBBiy7tpLaFiY1SFZprO9hUNP3UVx/bP0Ec7+nH9g4VXf9VsnojeMw8HxA9to3PkcM2tfIc9RjcMIR8MX0l5wOfnnXceMnDlWhxmQRGSXMeYjHQgeJYCpQhPA+HR3tRH+kwJ2Zt/I2V+63+pwfO7d53/Jqn3fY/eaBzhr/U1WhzMlVB7ZRfOf/52zut6ig2gOZV/H7Mtun/QBZA67nWMlr9G2+4/kNrzGTNOA3QjvRSyha/YnmX3edaTODK4R+lYYLQFM794i5ZHo2ESOh80irnGX1aH4XF9vNzn77qc0dA7LL73B6nCmDGcH4584fmA7HX/5AaurH8Px8OPsjl2LbcXnWHDulX67dbSpvooTO15Cjr9OQce7zKedfhPKkegVnCy8jdlrr2ZReuCvzT0daAIIEq0py1lU/yK9vT1ERkZZHY7P7HnhAdbQyIHzAmvQl6/MXrIGlvzZeWfNll8wv/5FEt96i463vs6+hLVQcAHZyy8mM298gyQddjvVpfs4deTvmJoS0lv3UuCoIBVoIZ7y+FWcmLeeeWs/w7KEZN+enJowvQQUJPb/9QmWvvM1Dq7/A4vXXGp1OD7R19NF532LaAjPYcF33tIE4IH+vl6O/H0T/fv/SGHb30ikC4AGkjkVkU93/CxMYi62mFTC41MIsTmnL3DYBxnoamawqxlzupHw9goSe6rIHDxJtPQB0GmiqIicT1fWeaQtW8+sxWum5fw4gUgvAQW5/LMuhXeg48gbECAJYN+mB1lFG3XrHtQPfw+FR0Sy7MLPwoWfxWG3U36khIYDr2Gr30Pi6RPManiJmMbeMx5j0IRQH5JBc2QO+1OKCclcQsaCteQULmOJfuBPK5oAgkR86gxO2PKJqw+MqaEH+nvJPfIwR8IWsvicy60OZ1oKsdmYtXg1sxZ/MLDfOBx0dLTS2VJPV8spjMO1NnFICDEJqcQlZRCXmEp2aCjTd45VNUQTQBBpSC5macMmBvp7CQuf3nPH7HvpVxTTRP05P9a//n1IQkKIT0whPjEFZi2yOhzlZ/qbE0RCZ68jSvo5sfctq0OZEPvgAJkHfs0xWyHLzv+M1eEoNW1pAggieWc5r/23HXnD4kgmZu/Lj5Bl6ula/S/6179SE6C/PUEkNT2T4yH5RNdO334Ah91O2p5fUh6SR9FF11sdjlLTmiaAINOYsoo5vQfp7em2OpRx2ffKE+Q6amg+63a9xVCpCdIEEGSi5p5PpAxQuudNq0PxmnE4SNx5P1WSxfKP32h1OEpNe5oAgszslR/HYYT2Q69aHYrXDmx9hgJ7BXXL/jmoJ3xTylc0AQSZ2MQ0jofPJfXU21aH4hXjcBD1zn9zUjI46/IvWR2OUgFBE0AQap6xlsKBo7Q1N1gdiscOv/0ihYPHqFr4ZcLCI6wOR6mAoAkgCCUuWY9NDOXv/tnqUDwW8vZPOUUKyz8Z2MtaKjWZNAEEodnLL6DDRGM/Nj36Ad7bvpkF/Qcpn/uPREZFWx2OUgFDE0AQCgsLpzRmBblt2zEOh9XhjGlw609oIpFlG75mdShKBRRNAEFqIP8CMkwTVcf2WR3KGZXtfp3Fvbs5NutGomPirA5HqYCiCSBI5a+5AoCTJX+yOJIzO/3qfbQRy5Irv251KEoFHE0AQWpG7lwqQ3KJq3rN6lBGVX7gHZZ1b+dQ7ueJi0+yOhylAo4mgCBWl3khC/r2094yNW8HbdvyYzpNFIuv/JbVoSgVkDxKACKyXkSOikiZiHxnhO0RIvKMa/sOEcl3lX9ORPa6PRwiUuTattV1zKFt6b48MTW2pOUbCBUHpW//0epQPqLyvV0Udb7FgezrSEhOtTocpQLSmAlARGzAg8BlwELgehFZOKzaLUCrMWYO8HPgPgBjzJPGmCJjTBHwBaDCGLPXbb/PDW03xkzNP0MDWOHy82kkiZBjm60O5SMaN/8/egln/oZ/tToUpQKWJ98AVgFlxphyY0w/sBHYMKzOBuC3rufPAReJiAyrcz3w9ESCVb4VYrNxPOk85nXtYKCvx+pw3ldz/BDL219l/4yrSE6faXU4SgUsTxJAFlDt9rrGVTZiHWPMINAOpAyrcy0fTQCPuS7//NsICQMAEblVREpEpKSxsdGDcJU3whd/ihh6eW/b1BkVfPKl/2SQUOZceZfVoSgV0DxJACN9MBtv6ojIaqDbGHPQbfvnjDFLgPNcjy+M1Lgx5mFjTLExpjgtLc2DcJU3Fp7zSU6bSE7ve8HqUACorSxlecvL7Eu/gtTMPKvDUSqgeZIAaoAct9fZQO1odUQkFEgAWty2X8ewv/6NMSddPzuBp3BealKTLDIqmqMJ5zCvdSv9fX1Wh0P1ph8iQN4V+te/Uv7mSQLYCRSKSIGIhOP8MN80rM4mYGiFjquB140xBkBEQoBrcPYd4CoLFZFU1/Mw4JPAQZQlbEuvIYlOjrwz/G2dXPWVR1ne9Cd2pX6KjJxCS2NRKhiMmQBc1/RvA7YAR4BnjTGHROReEbnCVe1RIEVEyoA7AfdbRdcBNcaYcreyCGCLiOwH9gIngf+d8NmocZl/3qdpJ4bBvc9aGkf1C/+BIYS8T99jaRxKBYtQTyoZYzYDm4eV3e32vBfnX/kj7bsVWDOs7DSwwstYlZ9ERESxJ/FjLGl9ld7uTiKjJ3/OnerSfSxveZmSjGtYkz1r0ttXKhjpSGAFQGzxdcRILwdef8aS9us33Us/Ycy56t8saV+pYKQJQAGwcM1lnJIUwg9M/lCN4wffZUXHaxzIupbUGTlj76CU8glNAAqAkNBQTmR/miW9u6itODqpbXduvptuiWTBZ74/qe0qFew0Aaj35V38ZQCqX/v1pLV55J2XKOrexoGCW4hPyZi0dpVSmgCUm8y8ueyPWsms6ucZGOj3e3t2u52w1+6mnlSKrvmu39tTSn2YJgD1IbLiJtJoZe8rT/m9rZJNv2KO/Tg1K/6VqOgYv7enlPowTQDqQ5Z87LPUSgZxu38FZviMH77T2d5Cwb6fURo6lxWf+Ee/taOUGp0mAPUhIaFhVM37IvMH3+Pozr/6rZ3DT91FqmnFXHYfEmLzWztKqdFpAlAfseST/0wLcfRt/Zlfjl+27+8U1z/DjpQrmLviQr+0oZQamyYA9RExsfEcy/scS7t3cGzv2z499uDAAIN/+jptEs/Cz/snwSilPKMJQI1o8VX/Shux9P/FtyNzd228l/mDR6le+X0SknV6b6WspAlAjSg2IYXDc77M4t7dHHjTN2sGl+1/h+VlD7IrZh3LLtOOX6WspglAjWrF1d+kVjKI3/pv9PWentCxurvasL3wZdolnlk3/y8Sov/1lLKa/haqUUVERtN8wf8jz9Sw9/ffG/dxjMPB0V/fQK69moaLHyApdYYPo1RKjZcmAHVGS87/DDsS1rOi+rcc3vHquI7x98e/y/KuN9kx+3YWrd3g4wiVUuOlCUCNaeHNv+RUSBrpL99C08kTXu377rP/xdqqX7E74RLO/vy/+ydApdS4aAJQY4pLTKP7qt8TbXpof/TTtDfVjbmPcTjY8eS9rDr8I/ZGrWHxPz+p1/2VmmL0N1J5pHDJKo5f9Guy7DW0/s8l1JTuG7Vud2cr7/7i86wu/Rm7Ys5j4R3PEx4RMYnRKqU8oQlAeWzJuqso//jjxDvaSP79JWx79Bs0ul0S6mg+xY5n7qPjZ8WsbN3MjuybKbrzRcIjoy2MWik1GjF+nPDL14qLi01JSYnVYQS9hpMnqH76DlZ0vQlAI0nYsZFmmrGJ4b3QBXDJvcxffanFkSqlAERklzGmeHi5R4vCi8h64AHABjxijPnxsO0RwBM4F3pvBq41xlSISD5wBBhaYmq7MeafXPusAB4HonAuOH+HmU7ZKIilZxWQ/s1N1JTuo2bbc0hLGWIcVCbmkVz0CeYtO0+v9ys1DYyZAETEBjwIXALUADtFZJMx5rBbtVuAVmPMHBG5DrgPuNa17bgxpmiEQ/8KuBXYjjMBrAdeHveZqEmXXbiM7MJlVoehlBonT/5MWwWUGWPKjTH9wEZg+M3cG4Dfup4/B1wkIjLaAUUkE4g3xmxz/dX/BHCl19ErpZQaN08SQBZQ7fa6xlU2Yh1jzCDQDqS4thWIyB4ReVNEznOrXzPGMQEQkVtFpEREShobGz0IVymllCc8SQAj/SU//Fr9aHXqgFxjzHLgTuApEYn38JjOQmMeNsYUG2OK09J09killPIVTxJADZDj9jobqB2tjoiEAglAizGmzxjTDGCM2QUcB+a66mePcUyllFJ+5EkC2AkUikiBiIQD1wGbhtXZBNzoen418LoxxohImqsTGRGZBRQC5caYOqBTRNa4+gpuAF70wfkopZTy0Jh3ARljBkXkNmALzttAf2OMOSQi9wIlxphNwKPA70SkDGjBmSQA1gH3isggYAf+yRjT4tr2FT64DfRl9A4gpZSaVDoQTCmlAtxoA8F0tI5SSgWpafUNQEQagcpx7p4KNPkwnOlAzzk46DkHh4mcc54x5iO3UU6rBDARIlIy0legQKbnHBz0nIODP85ZLwEppVSQ0gSglFJBKpgSwMNWB2ABPefgoOccHHx+zkHTB6CUUurDgukbgFJKKTeaAJRSKkgFRQIQkfUiclREykTkO1bH4wsikiMib4jIERE5JCJ3uMqTReQVESl1/UxylYuI/ML1b7BfRM6y9gzGT0RsrinGX3K9LhCRHa5zfsY1ZxUiEuF6Xebanm9l3OMlIoki8pyIvOd6v88O9PdZRL7u+n99UESeFpHIQHufReQ3ItIgIgfdyrx+X0XkRlf9UhG5caS2RhPwCcBtRbPLgIXA9SKy0NqofGIQ+IYxZgGwBviq67y+A7xmjCkEXnO9Buf5F7oet+JckW26ugPnUqND7gN+7jrnVpwr1IHbSnXAz131pqMHgL8YY+YDy3Cee8C+zyKSBdwOFBtjFuOcg2xopcFAep8fx7kSojuv3lcRSQbuAVbjXLzrnqGk4RFjTEA/gLOBLW6v7wLusjouP5zniziX7TwKZLrKMoGjrucPAde71X+/3nR64Jw6/DXgQuAlnGtLNAGhw99vnBMYnu16HuqqJ1afg5fnGw+cGB53IL/PfLDAVLLrfXsJ+Hggvs9APnBwvO8rcD3wkFv5h+qN9Qj4bwB4tqLZtOb6yrsc2AFkGOd027h+pruqBcq/w/3AvwIO1+sUoM04V6KDD5/XmVaqmy5mAY3AY67LXo+ISAwB/D4bY04CPwWqcC4q1Q7sIrDf5yHevq8Ter+DIQF4vPrYdCQiscD/Af9ijOk4U9URyqbVv4OIfBJoMM7Fhd4vHqGq8WDbdBEKnAX8yjhX1jvNB5cFRjLtz9l1CWMDUADMBGJwXgIZLpDe57GMdo4TOvdgSACerGg2LYlIGM4P/yeNMX90FZ8SkUzX9kygwVUeCP8O5wJXiEgFsBHnZaD7gURxrkQHHz6vEVeqm8yAfaAGqDHG7HC9fg5nQgjk9/li4IQxptEYMwD8ETiHwH6fh3j7vk7o/Q6GBODJimbTjogIzoV4jhhj/tttk/vqbDfywUprm4AbXHcTrAHah75qThfGmLuMMdnGmHyc7+PrxpjPAW/gXIkOPnrOH1mpbhJDnjBjTD1QLSLzXEUXAYcJ4PcZ56WfNSIS7fp/PnTOAfs+u/H2fd0CXCoiSa5vTpe6yjxjdSfIJHW0XA4cw7km8fesjsdH57QW51e9/cBe1+NynNc+XwNKXT+TXfUF591Qx4EDOO+wsPw8JnD+FwAvuZ7PAt4FyoA/ABGu8kjX6zLX9llWxz3Ocy0CSlzv9QtAUqC/z8B/AO8BB4HfARGB9j4DT+Ps4xjA+Zf8LeN5X4Evus69DLjZmxh0KgillApSwXAJSCml1Ag0ASilVJDSBKCUUkFKE4BSSgUpTQBKKRWkNAEopVSQ0gSglFJB6v8HN7ZuWOkcotsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ts, simulation_results['u'])\n",
    "plt.plot(ts, simulation_results['ulimited'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the effect of the limiter clearly in the above figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}