IF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The integrate and fire neuron\n",
    "\n",
    "In this notebook we will use Python to simulate the integrate and fire (I&F) neuron model.  We'll investigate, in particular, how the spiking activity varies as we adjust the input current $I$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Background information about the I&F model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a video that describes a slightly more complicated model, the *leaky* integrate and fire model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"400\"\n",
       "            height=\"300\"\n",
       "            src=\"https://player.vimeo.com/video/140084447\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.VimeoVideo at 0x11c92a710>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.lib.display import VimeoVideo\n",
    "VimeoVideo('140084447')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's some additional intereting videos and references:\n",
    "\n",
    "- [Lecture by Prof. Gerstner](http://klewel.com/conferences/epfl-neural-networks/index.php?talkID=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Preliminaries\n",
    "Before beginning, let's load in the Python packages we'll need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize']=(12,3)                # Change the default figure size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Part 1:  Numerical solutions - Introduction\n",
    "\n",
    " \n",
    "How do we compute a numerical solution to the integrate and fire model?\n",
    "  The basic idea is to rearrange the differential equation to get $V(t+1)$ on\n",
    "  the left hand side, and $V(t)$ on the right hand side.  Then, if we know\n",
    "  what's happening at time $t$, we can solve for what's happening at time $t+1$.\n",
    "\n",
    "For example, consider the differential equation:\n",
    "\n",
    "$$\n",
    "  \\dfrac{dV}{dt} = \\dfrac{I}{C}\n",
    "$$\n",
    "\n",
    "In words, we can think of:\n",
    "\n",
    "$dV$ as the \"change in voltage V\",\n",
    "\n",
    "$dt$ as the \"change in time t\".\n",
    "\n",
    "Let's consider the case that we record the voltage $V$ in discrete time steps. So\n",
    "we observe:\n",
    "\n",
    "$V[0], V[1], V[2], \\ldots$\n",
    "\n",
    "at times:\n",
    "\n",
    "$dt, \\, 2*dt, \\, 3*dt, \\ldots$\n",
    "\n",
    "where $dt$ is the time between our samples of $V$.\n",
    "\n",
    "We can now write the \"change in voltage V\" as:\n",
    "\n",
    "$$\n",
    "  dV = V(t+1) - V(t)\n",
    "$$\n",
    "\n",
    "Notice that the change in voltage is the difference in V between two\n",
    "sequential time samples.  Now, let's rewrite $\\dfrac{dV}{dt}$ as,\n",
    "\n",
    "$$\n",
    "  \\dfrac{dV}{dt} = \\dfrac{ V(t+1) - V(t) }{ dt } \n",
    "$$\n",
    "\n",
    "where we've replaced $dV$.  Now, let's substitute this expression into the equation at the top of this file:\n",
    "\n",
    "$$\n",
    "  \\dfrac{ V(t+1) - V(t) }{ dt } = \\dfrac{I}{C}.  \n",
    "$$\n",
    "\n",
    "Solving this equation for $V(t+1)$ you'll find  that:\n",
    "\n",
    "$$\n",
    "  V(t+1) = V(t) + dt*(I/C)\n",
    "$$\n",
    "\n",
    "Notice that, in this expression, we use our current value of the voltage V(t) and the model (I/C) to determine the next value of the voltage V(t+1).\n",
    "\n",
    "Now, let's program this equation in Python.  First, let's set the values\n",
    "for the parameters $I$ and $C$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "C=1.0\n",
    "I=1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  We also need to set the value for $dt$.  This defines the time step for our\n",
    "  model.  We must choose it small enough so that we don't miss anything\n",
    "  interesting.  We'll choose:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt=0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's assume the units of time are seconds.  So, we step forward in time by $0.01$ s.\n",
    "\n",
    "The right hand side of our equation is nearly defined, but we're still missing one thing, $V(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  What value do we assign to $V(t)$?\n",
    "\n",
    "**A:**  We don't know --- that's why we're running the simulation in the first place!\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So here's an easier question:  what *initial* value do we assign to $V(t)$?\n",
    " \n",
    " To start, we'll create an array of zeros to hold our results for $V$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = zeros([1000,1])\n",
    "V.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This array `V` consists of 1000 rows and 1 column. We can think of each row entry as corresponding to a discrete step in time. Our goal is to fill-in the values of `V` (i.e., step forward in time), in a way consistent with our model.\n",
    "\n",
    "Let's choose an initial value for `V` of 0.2, which in our simple model we'll assume represents the rest state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "V[0]=0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  Given the initial state `V[0]=0.2`, calculate `V[1]`.  Then calcualte `V[2]`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  After the two calculations above, we've moved forward two time steps into\n",
    "  the future, from $t=0$ s to $t=0.01$ s, and then from $t=0.01$ s to $t=0.02$ s.  But what\n",
    "  if we want to know $V$ at $t=10$ s?  Then, this iteration-by-hand procedure becomes\n",
    "  much too boring and error-prone.  So, what do we do?  Let's make the\n",
    "  computer do it ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2:  Numerical solutions - implementation\n",
    "\n",
    "Let's computerize this iteration-by-hand procedure to find `V[999]`. To do so, we'll use a [for-loop](https://wiki.python.org/moin/ForLoop).  Here's what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(1,999):\n",
    "    V[k+1] = V[k] + dt*(I/C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  Does this loop make sense?  Describe what's happening here.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  Why does the `range` command end at `999`?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute this for-loop and examine the results in vector `V`.  To do so, let's plot `V`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAADCCAYAAABAOqrYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3daWyd2X3f8e8jkdS+jESNJEoaraQkkt5mxmN7xuPZF8+MPXbrnNhpXMdxMG+apWmKNAlaGGgT1AWCJC4QBBg4aezWcHLqGIjbGnECt2lQtHHsscewSEm8FLVLM9r3hdvpi/NQlGRpRuJ2yXu/H4Cg7j2XvH/o6JH+OjzP+RUpJSRJkqR6MKvaBUiSJElTxeZXkiRJdcPmV5IkSXXD5leSJEl1w+ZXkiRJdcPmV5IkSXWjYYrfz3PVJEmSNBWKWz051c0vR44cmeq3pLm5mRMnTkz5+2pqOc/1wXmuD85zfXCe60M15rmlpeW2Y257kCRJUt1425XfEMKfAC8Bx2KMneVzy4A/BzYA+4AQYzw9eWVKkiRJ43cnK79/Cjx/03O/AXwnxtgKfKd8LEmSJE1rb9v8xhj/Djh109MvA18uf/1l4GMTXJckSZI04ca653dljPEoQPn53okrSZIkSTNVunKJtOM1hr/xFYa//Y1ql/MTJv20hxDCK8ArADFGmpubJ/stf0JDQ0NV3ldTy3muD85zfXCe64PzXBuGz5+jf+ePGOh+nf6uHzLU1wPDwzB7NnM+8MS0m+exNr9vhhBWxxiPhhBWA8du98IY46vAq+XDVI0jTTxKpT44z/XBea4PznN9cJ5npnTuNFS6Sbt3kCpdcHg/pAQNjbCpjeKFn6Jo7YDN2xicM5fBwcFpddTZWJvfbwKfAb5Qfv7LMX4fSZIkTWPp1AlSzw6odOXPbxzOA01zYMt2igceoWjrgI1tFI1N1S32DtzJUWdfAx4HmkMIh4DPk5veGEL4HHAA+KnJLFKSJEmTL6UEx9/IK7o9ZbN74s08OG8+bGmneORpirZOuG8zRcOU56WN29tWHGP81G2GnprgWiRJkjSFUkrwxiFSTxf07Mifz5zMgwsXQWsHxVMfySu7azdQzJpd3YInwMxr1yVJkjQmaXgYDu8njTS6lS44fzYPLlmWm9y2DorWTli9lmJW7YUB2/xKkiTVqDQ0BAf2kHq68laGShdcupgHl99L0Xl/Xt3d2gkrVlMURXULngI2v5IkSTUiDQzAvsroyu6eXXD1ch5cuYbigUfKld0OiuX1GdNg8ytJkjRDpatXoW/X6Mpu324Y6M+Da9ZTPPwEtHZStLZTLF1W3WKnCZtfSZKkGSJdvgS9O0mVcmV3Xy8MDUIxC9ZtpHjsw3nfbms7xcLF1S53WrL5lSRJmqbShXPQ203aXa7sHuiDlNPT2NBK8czL+dixzdso5i+odrkzgs2vJEnSNJHOni5PYShXdg/vzwMNjbBpK8WLIa/sbtpKMWdudYudoWx+JUmSqiSdOp6DJHq6crP7ZpmeNmcubN5O8d5Hc1TwxjaKxsbqFlsjbH4lSZKmQE5PO3pjoMTJY3lw3oK8T/fRZ/PK7rpNMzI9bSbwd1WSJGkSpJTg6MHRld1KF5w5lQcXLs5Hjj3zcl7ZXbu+JtLTZgKbX0mSpAmQhofgUJmeVumCni64cC4PLl2Wb0wbCZRYtbYuAiWmI5tfSZKkMUiDgzk9rdJF2r0DenfC5TI9rXklxTsehK2deWV3xSqb3WnC5leSJOkOpIF+2FsZXdndswuuXsmDq9ZQPPgItHVStHVQLFtR3WJ1Wza/kiRJt5CuXoE9u/LKbk+ZnjY4kAfXrKd4+Kl8c1pbB8Xie6pbrO6Yza8kSRKQLl2EPTtzVHDPDtjfC0NDOT3tvk0UT7xQ7tttp1iwqNrlaoxsfiVJUl1KF85BpTtvY+jpgoN7y/S0BtiwheLZj5Xpadsp5s2vdrmaIDa/kiSpLqQzp66dwpAq16WnNTblxLSXQr45bdM2ijlzqlusJo3NryRJqknp5LEbAyWOHckDc+bBlm05PW1rJ6xvNT2tjtj8SpKkGS+lBMeO3hgoMZKeNn9BPl/3Q8/lbQz3baKYbaBEvbL5lSRJM04aHoajh3KzO3Iaw9kyPW3RktzsPvOxfBrDmvUUs2ZVt2BNGza/kiRp2kvDQ3BwH6myg7S7C3q74ML5PLh0OcXWd+Qjx9o685m7BkroNsbV/IYQfhX4BSABPwY+G2O8MhGFSZKk+pUGB2F/7+gZu73dcPlSHlyxiuKdD402u80rbXZ1x8bc/IYQ1gC/DLTHGC+HECLwSeBPJ6g2SZJUJ3J6Ws/osWN7dkH/1Ty4ai3Fez+Um93WdtPTNC7j3fbQAMwLIQwA84Ej4y9JkiTVunT1Cldf/weGv/9/877dvT0wOAhFAWs2UHzwmbxft7WDYvHSaperGjLm5jfGeDiE8LvAAeAy8Ncxxr+esMokSVLNSJcuQO916WkH9nBmaAhmzYL7NlM8+VLewrBlu+lpmlRFSmlMXxhCuAf4C+CngTPAfwW+HmP8Lze97hXgFYAY4wP9/f3jKngsGhoaGBwcnPL31dRynuuD81wfnOeZb/jcGfq7f8RA1w/p736dwb0VSAkaGmjc0k5Tx7uZ+84HmNXazqx5C6pdriZRNa7npqYmgFtuBB/Ptoengb0xxuMAIYRvAA8DNzS/McZXgVfLh+nEiRPjeMuxaW5uphrvq6nlPNcH57k+OM8zTzpzMu/VrXSRdu+AowfzQFNTTkx76ZN5G8OmrQw3zeEKsHBkni9ermrtmlzVuJ5bWlpuOzae5vcA8P4QwnzytoengO+P4/tJkqQZIp14czQ9rdIFx47mgTnzoHU7xfsfz9sYNmyhaDA9TdPHePb8fjeE8HXgB8Ag8ENGV3glSVKNSCnBm0duDJQ4dTwPzl8Ire0Uj304r+yuMz1N09u4TnuIMX4e+PwE1SJJkqaBNDwMRw7kFd2RqOCzp/PgoiV5Rfe5j+dmt8X0NM0sJrxJklTn0tAQHNo7ehJDpRsululp9zRTbHvnaKDEStPTNLPZ/EqSVGfS4ADs3zMaKNHbDVfKm85WrKJ490PQ1knR2mF6mmqOza8kSTUu9V8t09PKLQx7dsLI0aOr11G877EcJtHWSXHP8uoWK00ym19JkmpMunIZ9uwa3caw77r0tLUbKB59Lq/qtrabnqa6Y/MrSdIMly5dgMrOvI2h0gX7e2F4OKenrd9C8eRHrktPW1jtcqWqsvmVJGmGSefPXjuFIfXsgEP7rqWnsaGN4vlP5JMYNm+lmDu/2uVK04rNryRJ01w6ffLGM3avT0/bvJ3iI5/Kze7GNoqmOdUtVprmbH4lSZpGUkpw4s3yjN3yNIbjb+TBufNgSzvFB57I2xjWbzY9TbpLNr+SJFVRTk87nFd2R05jOHUiDy5YlG9Ke+LFvLK7dqPpadI42fxKkjSFcnra/ryiO9LsnjuTBxcvzSu6z5eBEqvXmZ4mTTCbX0mSJlEaGoKDfaOBEpVuuHQhDy5rpmh/92igxMoWAyWkSWbzK0nSBEqDA7Cvd/TYsd6do+lp966muP8DZaBEB0XzyuoWK9Uhm19JksYh9V+Fvt2jgRJ7d9+Ynvb+x0eb3aWmp0nVZvMrSdJdSFcuQe+u0ZXdvRUYKtPT1m3M6WltnflGtUVLql2upJvY/EqS9BbSxQv5fN2RM3YP7LkxPe3pj+aTGLZsp5hvepo03dn8SpJ0nXTuzLUwidSzAw7vH01P29hG8eEyPW3TNoq586pdrqS7ZPMrSapr6dSJGwMl3jiUB5rmwOZtFB/9VN7GsLGNorGpusVKGjebX0lS3biWnnZ9oMRIetq8+Tk97ZGn8rFjpqdJNcnmV5JUs1JK8MZ16Wk9O+DMyTy4cBFs6aB48kWK1k5Yt4FilulpUq2z+ZUk1Yw0PAyH94/u1610wfmzeXDJPeUpDCPpaWtNT5PqkM2vJGnGSkNDcKBv9NixShdcupgHl99L0XE/tJXN7r2rTU+TNL7mN4SwFPgS0Akk4OdjjP9vIgqTJOlmaWAA9lWuS0/bBVdH0tNaKO5/OEcFt3VQLL+3usVKmpbGu/L7ReCvYoyfCCE0AfMnoCZJkgBIV69C367RM3b7dsNAmZ7Wch/FB57IK7utHRRLl1W3WEkzwpib3xDCYuBDwM8BxBj7gf6JKUuSVI/S5UuwZyfnD/Yx9KPvwb7eMj1tVk5Pe+z5vIVhSzvFosXVLlfSDFSklMb0hSGEdwOvAt3Au4DXgF+JMV686XWvAK8AxBgf6O+f+v64oaGBwcHBKX9fTS3nuT44z7Vl+Pw5+rtfZ6D7dfq7Xmdwb09OT5s9m8bN22jseA9NHe+mcds7mbXA9LRa4/VcH6oxz01NTQC33OQ/nub3QeDvgUdijN8NIXwROBdj/Ddv8WXpyJEjY3q/8WhububEiRNT/r6aWs5zfXCeZ7Z07nR55Nh16WkADY2waWveq9vWSfN7H+bkhYtv/c0043k914dqzHNLSwvcpvkdz57fQ8ChGON3y8dfB35jHN9PklRj0qnjea/uyA1qbxzOA3Pm5vS0Bz9Ypqe13pCeVsydBza/kibBmJvfGOMbIYSDIYStMcbdwFPkLRCSpDqUUoLjb9wYKHHyWB6ctwC2bKf44DM5Pe2+zRQNnrYpaeqN92+eXwK+Wp700Ad8dvwlSZJmgpQSHD2YV3YrI+lpp/LgwsX5FIanP5pXdteuNz1N0rQwruY3xvg68OAE1SJJmsbS8BAc2l8eO5ZXd7lwLg8uWUbR1jEaKLF6nYESkqYlf+YkSbqlNDgIB8v0tJ4u6O2+MT3tHQ9cC5RghelpkmYGm19JElCmp+3tGV3Z3bMLrl7JgyvXUDzwSBko0UmxfEV1i5WkMbL5laQ6la5egb7do8eO9e2GwYE8uGY9xcNPQmsZFbzknuoWK0kTxOZXkupEunQR9uzMzW6lC/ZVYGgop6fdt4ni8Rcotnbk9LSFpqdJqk02v5JUo9KFc1DpHl3ZPbgXUk5PY0MrxbMfo2jtzEeQzZtf7XIlaUrY/EpSjUhnT98YKDGSntbYlNPTXgz55rRN2yjmzKlusZJUJTa/kjRDpZPH84puJccF8+b16WnbKd77aD52bEMrRWNjdYuVpGnC5leSZoCUEhw7emOze316Wms7xaPP5pXddZtMT5Ok2/BvR0mahlJKcOQgqTISFdwFZ8v0tEVLoLWD4pmX88rumvtMT5OkO2TzK0nTQE5P2zcaKFHpHk1PW7qMYmvnaKDEqrUGSkjSGNn8SlIVpMFB2N9bBkp0Qe9OuFympzWvpHjne8tAiQ5YscpmV5ImiM2vJE2BNNCf09N6rktP67+aB1etpXjvB/NWhrYOimWmp0nSZLH5laRJkK5egT278jaGShf09Yymp63dQPHI03krQ2s7xWLT0yRpqtj8StIESJcuQm/3aHra/t4b09OefDFvYWhtp1iwqNrlSlLdsvmVpDFI58/lI8cqI+lp+8r0tAbY2Erx7MfzzWmbTU+TpOnE5leS7kA6cyqv6I6cxnDkQB4YSU97KeRjxzZuNT1NkqYxm19JuoV08hhp90igxA44djQPzJkHW7ZRvO+xvLK73vQ0SZpJbH4l1b2UErx55MaV3VPH8+D8hXmf7mPPU7R25v27sw2UkKSZyuZXUt1Jw8Nw9GBe0R25Qe3s6Ty4aEk+X/fZj1Ns7YCW9RSzZlW3YEnShLH5lVTz0vAQHLwuPa23Cy6cz4NLl1NsfSds7cgru6vWGCghSTVs3M1vCGE28H3gcIzxpfGXJEnjcy097VqgxE64fCkPrlhF8a6HoLWMCm5eabMrSXVkIlZ+fwXYCSyegO8lSXctDfRDX89ooMRPpKd96FpUcLGsubrFSpKqalzNbwhhLfAi8DvAv5iQiiTpbaQrl6FvF2l3F6myA/b2wOAgFAWs2UDxwWfysWOt7RSLl1a7XEnSNDLeld8/AH4dMK5I0qRJly5AZSepsoOTfbsZ3rMLhodh1iy4bzPFkx/JWxi2tFMsWFjtciVJ09iYm98QwkvAsRjjayGEx9/ida8ArwDEGGlunvofOTY0NFTlfTW1nOfaMXz2NP3dP6K/64cMdL/O0L5eSAkaGpnV1s6Cj/8sjZ3voXFrJ7PmLah2uZoEXs/1wXmuD9NtnouU0pi+MITw74FPA4PAXPKe32/EGH/2Lb4sHTlyZEzvNx7Nzc2cOHFiyt9XU8t5nrnSmZP5FIaR0xiOHswDTU2waRtFW3lz2sY2VrSscZ7rgNdzfXCe60M15rmlpQXglnczj3nlN8b4m8BvApQrv//ybRpfScqBEiePjZ6x27MDjr+RB+fOgy3bKT7wBEVrB2zYQtFgepokaeJ4zq+kSZXT0w6PruxWuuBUuQIwf2E+heHxFyi2dsLajaanSZIm1YQ0vzHGvwX+diK+l6SZLQ0Pw5EDN6annTuTBxcvzSu6z//j/LnlPtPTJElTypVfSeOShobgYF8OlKh0QaUbLpbpafc0U2x/F4zs2V1pepokqbpsfiXdlTQ4APt6SZVyv27vTrhyOQ+uWEXx7vflrQxtnbD8XptdSdK0YvMr6S2l/quwt4e0u9yv27cL+vvz4Op1FO97LK/stnZQ3LO8usVKkvQ2bH4l3SBduQS9u0ZXdvdWYKhMT1u7geLR5/IWhtYOikVLql2uJEl3xeZXqnPp4gXo7c57dnt2wIE9o+lp67dQPP0RitZOaN1OMd/0NEnSzGbzK9WZdO4MVLpJI4ESh/eV6WkNsLGN4vlPUGztyOESc+dVu1xJkiaUza9U49LpkzcGSrxxKA80NcHm7RQf+VS+OW1jK0XTnOoWK0nSJLP5lWpISglOvJlXdCvlyu4N6WntFA8/lffsrt9sepokqe7Y/EozWEoJ3jhMquyA3eU5u6fL9LQFi/JNaU+8mFd2122gmGV6miSpvtn8SjNITk/bT9rdlRveni44fzYPLl6am9yRQInV60xPkyTpJja/0jSWhobgQB9pZAtDpQsuXcyDy1ZQdNyfAyVaO2Bli4ESkiS9DZtfaRpJAwOwvzIaKNG7C66W6Wn3tlDc/3DeyrC1k2L5vdUtVpKkGcjmV6qidPUq9I0ESnRB324YKNPTWu6j+MDjZXpaO8VS09MkSRovm19pCuX0tJ2jgRL7ekfT09ZtpHjs+byFobWDYtHiapcrSVLNsfmVJlG6eP7GQIkDfZCGYfbsMj3to/nmtC2mp0mSNBVsfqUJlM6dLsMkymPHDu8v09MaYVMbxQufyCcybNpqepokSVVg8yuNQzp1/NopDDk97XAeaJoDm7dRfPRn8sruxjaKxqbqFitJkmx+pTuVUoLjb+QV3ZHTGE68mQfnzc/paY88nffsrt9C0eDlJUnSdOO/ztJt5PS0Q6TdO8qV3S44czIPLizT0556KW9jWGt6miRJM4HNr1RKw8NwaN/osWOV69LTltxTpqd1ULR2wuq1pqdJkjQD2fyqbuX0tD2jx471do+mpy2/l6Lz/ry629YJ9642PU2SpBow5uY3hLAO+AqwChgGXo0xfnGiCpMmWhoYgH2V0WPH9uyEq1fy4Mo1FA88Uja7HaanSZJUo8az8jsI/FqM8QchhEXAayGEv4kxdk9QbdK4XEtPG1nZ7dsNgwN5cM16ig88eV162rLqFitJkqbEmJvfGONR4Gj56/MhhJ3AGsDmV1UxfOki6cev5ZXdShfsq8DQEBSzcnra4y+UgRLtpqdJklSnJmTPbwhhA/Ae4LsT8f2kO5EunIPebtLuHChx/GAfDJfpaRtaKZ75WG52N2+nmL+g2uVKkqRpoEgpjesbhBAWAv8b+J0Y4zduMf4K8ApAjPGB/v7+cb3fWDQ0NDA4ODjl76uJNXT6JAPdr9Pf9ToDXT9k8EBfHmhsorGtg7md9zO7/V00tXWYnlbDvJ7rg/NcH5zn+lCNeW5qagK45Z3q42p+QwiNwH8Hvh1j/L07+JJ05MiRMb/fWDU3N3PixIkpf1+NTzp5nFTZcS0umDfL9LQ5c/NqbltHDpTY2EbR2Og81wnnuT44z/XBea4P1ZjnlpYWuE3zO57THgrgj4Gdd9j4SreV09OO5iZ35DSGk8fy4LwF0NpO8egz+dixdZtMT5MkSWMyng7iEeDTwI9DCK+Xz/1WjPFb4y9LtS6lBEcP5lMYevKeXc6cyoMLF+cwiWdeziu7a9ebniZJkibEeE57+D/cZjlZulkaHsrpaSPHjlW64cK5PLh0WV7RLc/YZfU6AyUkSdKk8GfHmhRpcDCnp1W6SLt3QO9OuHxdeto7Hsyru22dsGKVza4kSZoSNr+aEGmgH/ZWRs/Y3bNrND1t1RqKBx8pAyU6KJavqG6xkiSpbtn8akzS1SuwZ1de2e3ZAX09N6anPfxU3sLQ2kGx5J7qFitJklSy+dUdSZcuwp6dpN3lyu7+3tH0tPs2UTzxwmizu2BRtcuVJEm6JZtf3VI6fw4qXaMruwf3QRqG2Q2wYQvFsx/L+3U3b6eYN7/a5UqSJN0Rm18BkM6cyiu6I6cxHDmQBxqbYNNWipdCPnZs0zaKOXOqW6wkSdIY2fzWqXTy2I2BEsfK5L2R9LSHPpRXdje0UjQ2VrdYSZKkCWLzWwdSSnDs6GigRM8OOHU8D85fkPfpfui53Ozet4litoESkiSpNtn81qA0PFymp3Xlfbs9XXC2TE9btCQ3u89+PN+gtuY+09MkSVLdsPmtAWl4CA7uI1V2kHZ3QW8XXDifB5cup9j6jjJQogNWrTVQQpIk1S2b3xkoDQ7C/t4cFVzpgt5uuHwpD65YRfHOh0bT05pX2uxKkiSVbH5ngDTQD309eWW3p0xP67+aB1etpXjvo2V6WjvFMtPTJEmSbsfmdxpKVy5D3668stuzA/b2wOAgFEVOT/vgM2WgRDvFYtPTJEmS7pTN7zSQLl2A3p2kkWPHDuzJ6WmzZsF9mymefClvYdiy3fQ0SZKkcbD5rYKcnrZjdM/uwb2QUk5P29hansTQCVu2Ucw1PU2SJGmi2PxOgXTm5OixY7t3wNGDeeBaeton8zaGTVspmkxPkyRJmiw2vxMspQQnj40GSlS64NjRPDhnHrRup3j/42V62haKBtPTJEmSporN7zillODNI7nZrYykp53Ig/MX5pvSHns+N7vrTE+TJEmqJpvfu5SGh+HIgbyiu3tH/nzuTB5ctCSfr/vcP8rbGFrWU8yaVd2CJUmSdI3N79tIQ0NwaC9ppNGtdMPFMj3tnmaK7e8aDZRYucZACUmSpGms5pvf9OYRLu/8AWy//85ePzgA+3pJla58k1pvN1y5nAdXrKJ490PQ2plXdk1PkyRJmlHG1fyGEJ4HvgjMBr4UY/zChFQ1gdI//B3n/tvXmPUfv3bLY8NS/1XY2zMaKNG3C/r78+DqdRQPPZZXdls7KJY1T3H1kiRJmkhjbn5DCLOBPwSeAQ4B3wshfDPG2D1RxU2EYkNrvilt/x7Y+o6cnrZn12igxL7r09M2UDz6HEXrSHra0mqXL0mSpAk0npXfh4DeGGMfQAjhz4CXgWnV/LJhCwDDf/nV3OTu74Xh4Zyetn4LxZMfyVsYtrRTLFhY5WIlSZI0mcbT/K4BDl73+BDwvptfFEJ4BXgFIMZIc/MUbx1obubk+s0M7q3Q2NpO08d/lsbO99C4tZNZ8xZMbS2aVA0NDVP/50tTznmuD85zfXCe68N0m+fxNL+3utMr3fxEjPFV4NWR8RMnTozjLcdm+X/4EidOnmC4aQ5XgCsAFy/nD9WM5uZmqvHnS1PLea4PznN9cJ7rQzXmuaWl5bZj4zmE9hCw7rrHa4Ej4/h+k6aYM8fYYEmSJI1r5fd7QGsIYSNwGPgk8DMTUpUkSZI0Cca88htjHAR+Efg2sDM/FbsmqjBJkiRpoo3rnN8Y47eAb01QLZIkSdKkGs+eX0mSJGlGsfmVJElS3ShS+onTySbTlL6ZJEmS6tatjuWd8pXfohofIYTXqvXefjjPfjjPfjjPfjjP9fxRxXm+Jbc9SJIkqW7Y/EqSJKlu1Evz++rbv0Q1wHmuD85zfXCe64PzXB+m1TxP9Q1vkiRJUtXUy8qvJEmSNL6Et5kghPA88EVgNvClGOMXqlySxiCEsA74CrAKGAZejTF+MYSwDPhzYAOwDwgxxtMhhII87y8Al4CfizH+oBq16+6FEGYD3wcOxxhfCiFsBP4MWAb8APh0jLE/hDCH/OfiAeAk8NMxxn1VKlt3IYSwFPgS0Ek+BvPngd14PdeUEMKvAr9AnuMfA58FVuP1PKOFEP4EeAk4FmPsLJ+763+PQwifAf51+W1/O8b45amov6ZXfst/QP8Q+DDQDnwqhNBe3ao0RoPAr8UYtwPvB/5ZOZe/AXwnxtgKfKd8DHnOW8uPV4A/mvqSNQ6/Auy87vF/AH6/nOfTwOfK5z8HnI4xbgF+v3ydZoYvAn8VY9wGvIs8317PNSSEsAb4ZeDBskGaDXwSr+da8KfA8zc9d1fXb9ksfx54H/AQ8PkQwj2TXjk13vySfzN7Y4x9McZ+8v80X65yTRqDGOPRkf8pxhjPk/+hXEOez5H/KX4Z+Fj565eBr8QYU4zx74GlIYTVU1y2xiCEsBZ4kbwqSLlq8CTw9fIlN8/zyPx/HXiqfL2msRDCYuBDwB8DxBj7Y4xn8HquRQ3AvBBCAzAfOIrX84wXY/w74NRNT9/t9fsc8DcxxlMxxtPA3/CTDfWkqPXmdw1w8LrHh8rnNIOFEDYA7wG+C6yMMR6F3CAD95Yvc+5nrj8Afp28vQVgOXAmxjhYPr5+Lq/Nczl+tny9prdNwHHgP4UQfhhC+FIIYQFezzUlxngY+F3gALnpPQu8htdzrbrb67dq13WtN7+3+h+jx1vMYCGEhcBfAP88xnjuLV7q3M9AIYSRPWSvXff0W82l8zwzNQD3A38UY3wPcJHRH5HeivM8A5U/wn4Z2Ai0AAvIPwK/mddzbbvdvFZtvmu9+T0ErLvu8VrgSJVq0TiFEBrJje9XY4zfKJ9+c+THn+XnY+Xzzv3M9Ajw0SHnG9QAAAG/SURBVBDCPvI2pSfJK8FLyx+bwo1zeW2ey/El/OSP4jT9HAIOxRi/Wz7+OrkZ9nquLU8De2OMx2OMA8A3gIfxeq5Vd3v9Vu26rvXm93tAawhhYwihibzR/ptVrkljUO77+mNgZ4zx964b+ibwmfLXnwH+8rrn/2kIoQghvB84O/LjGE1fMcbfjDGujTFuIF+v/zPG+E+A/wV8onzZzfM8Mv+fKF/vStE0F2N8AzgYQthaPvUU0I3Xc605ALw/hDC//Dt8ZJ69nmvT3V6/3waeDSHcU/6U4NnyuUlX00edxRgHQwi/SP7NnA38SYyxq8plaWweAT4N/DiE8Hr53G8BXwBiCOFz5L9of6oc+xb5WJVe8tEqn53acjXB/hXwZyGE3wZ+SHmjVPn5P4cQeskrRJ+sUn26e78EfLVcmOgjX6Oz8HquGTHG74YQvk4+zmyQfO2+CvwPvJ5ntBDC14DHgeYQwiHyqQ139e9xjPFUCOHfkRcqAf5tjHFKVvpNeJMkSVLdqPVtD5IkSdI1Nr+SJEmqGza/kiRJqhs2v5IkSaobNr+SJEmqGza/kiRJqhs2v5IkSaobNr+SJEmqG/8fWTsntTwo6WUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure()\n",
    "plot(V);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  What happens to the voltage after 1000 steps?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  This plot is informative, but not great.  Really, we'd like to plot the\n",
    "  voltage as a function of *time*, not steps or indices.  To do so, we\n",
    "  need to define a time axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = arange(0,len(V))*dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  What's happening in the command above?  Does it make sense?  (If not, trying printing or plotting `t`.)\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, with *time* defined, let's redo the plot of the voltage with the axes labeled appropriately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAADTCAYAAACoXCQbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXCd13nf8e9LAuC+iAJFEtwXgCQA27I229p3y5Zsyal9bGfspnZcpdMmcdO0mdiZTjLNJHFmMmnUTqaNxkljJ66TE8et3Y4nduNM6qZJvciWxwRI4oLgTkrcRFFcsZ3+cS4IUiZFQgLuC+B+PzMYLAfUfeQjkI8fnvf8ipQSkiRJkq7PjLILkCRJkqYSG2hJkiRpDGygJUmSpDGwgZYkSZLGwAZakiRJGgMbaEmSJGkMGsou4HXw3j1JkiTVQnGlL07FBppDhw7V/DWbm5s5duxYzV9XteU+1wf3uT64z/XBfa4PZexzS0vLVdc8wiFJkiSNQU0m0CGEPwKeAI7EGDurX1sC/DmwDtgDhBjjS7WoR5IkSXq9ajWB/mPgsVd97ZeBb8YYW4FvVj+XJEmSJrWaNNAxxm8BJ1715SeBz1U//hzwVC1qkSRJkt6IMs9AL4sxHgaovr+pxFokSZI0SaTzZ0nbnmP4y59n+OtfLrucHzMlbuEIITwNPA0QY6S5ubnmNTQ0NJTyuqot97k+uM/1wX2uD+7z9DD8yin6t/+Qge7n6e/6AUN9PTA8DDNnMusdD0y6fS6zgX4xhLAixng4hLACOHK1b4wxPgs8W/00lXFdjdfk1Af3uT64z/XBfa4P7vPUlE69BJVu0s5tpEoXHNwLKUFDI2xoo3j3ByhaO2DjFgZnzWZwcHBSXWNXZgP9VeCngM9U33+lxFokSZI0QdKJY6SebVDpyu9fOJgXmmbBpq0Ut95F0dYB69soGpvKLfY61Ooauy8C9wPNIYQDwK+SG+cYQvhpYB/wgVrUIkmSpImTUoKjL+TJck+1YT72Yl6cMxc2tVPc9TBFWyes2UjRMCVOFF+mJhXHGD98laWHavH6kiRJmhgpJXjhAKmnC3q25fcnj+fF+QugtYPioffkCfOqdRQzZpZb8DiYei2/JEmSSpOGh+HgXtJIs1zpgldezouLluRGua2DorUTVqyimDH9gq9toCVJknRVaWgI9u0i9XTlYxmVLjh7Ji/eeBNF5y15yry5E5auoCiKcguuARtoSZIkXZQGBmBPZXTCvGsHXDiXF5etpLj1ruqEuYPixvqM8bCBliRJqmPpwgXo2zE6Ye7bCQP9eXHlWoo7H4DWTorWdorFS8otdpKwgZYkSaoj6dxZ6N1OqlQnzHt6YWgQihmwej3Ffe/K55hb2ynmLyy73EnJBlqSJGkaS6dPQW83aWd1wryvD1JO+WNdK8UjT+Yr5TZuoZg7r+xypwQbaEmSpGkkvfxS9XaM6oT54N680NAIGzZTPB7yhHnDZopZs8stdoqygZYkSZrC0omjOaykpys3zC9WU/5mzYaNWyluvyfHYq9vo2hsLLfYacIGWpIkaYrIKX+HLw8tOX4kL86Zl88t3/NonjCv3jAlU/6mAv9XlSRJmqRSSnB4/+iEudIFJ0/kxfkL83VyjzyZJ8yr1k6LlL+pwAZakiRpkkjDQ3CgmvJX6YKeLjh9Ki8uXpIf9hsJLVm+qi5CSyYjG2hJkqSSpMHBnPJX6SLt3Aa92+FcNeWveRnFm26DzZ15wrx0uQ3zJGEDLUmSVCNpoB92V0YnzLt2wIXzeXH5Sorb7oK2Toq2DoolS8stVldlAy1JkjRB0oXzsGtHnjD3VFP+Bgfy4sq1FHc+lB/4a+ugWHhDucXqutlAS5IkjZN09gzs2p5jsXu2wd5eGBrKKX9rNlA88O7qOeZ2inkLyi5Xr5MNtCRJ0uuUTp+CSnc+ktHTBft3V1P+GmDdJopHn6qm/G2lmDO37HI1TmygJUmSrlM6eeLi7RipcknKX2NTTvZ7IuQH/jZsoZg1q9xiNWFsoCVJkq4iHT9yeWjJkUN5YdYc2LQlp/xt7oS1rab81REbaEmSJKqhJUcOXx5aMpLyN3devn/53nfmIxlrNlDMNLSkXtlAS5KkupSGh+Hwgdwwj9yS8XI15W/BotwwP/JUviVj5VqKGTPKLViThg20JEmqC2l4CPbvIVW2kXZ2QW8XnH4lLy6+kWLzm/J1cm2d+U5mQ0t0FaU30CGEXwA+ASTgR8DHYozny61KkiRNdWlwEPb2jt7B3NsN587mxaXLKd58x2jD3LzMhlnXrdQGOoSwEvh5oD3GeC6EEIEPAX9cZl2SJGnqySl/PaNXyu3aAf0X8uLyVRS335sb5tZ2U/70hpQ+gSbXMCeEMADMBQ6VXI8kSZoC0oXzXHj+Owx/7+/zOebdPTA4CEUBK9dR3P1IPr/c2kGxcHHZ5WoaKbWBjjEeDCH8DrAPOAd8I8b4jTJrkiRJk1M6exp6L0n527eLk0NDMGMGrNlI8eAT+TjGpq2m/GlCFSml0l48hHAD8JfAB4GTwF8AX4ox/umrvu9p4GmAGOOt/f39tS6VhoYGBgcHa/66qi33uT64z/XBfZ76hk+dpL/7hwx0/YD+7ucZ3F2BlKChgcZN7TR13MzsN9/KjNZ2ZsyZV3a5mkBl/Dw3NTUBXPFgfNlHOB4GdscYjwKEEL4M3Alc1kDHGJ8Fnq1+mo4dO1bTIgGam5sp43VVW+5zfXCf64P7PPWkk8fz2eVKF2nnNji8Py80NeVkvyc+lI9kbNjMcNMszgPzR/b5zLlSa9fEKuPnuaWl5aprZTfQ+4C3hxDmko9wPAR8r9ySJElSLaRjL46m/FW64MjhvDBrDrRupXj7/flIxrpNFA2m/GnyKPsM9LdDCF8Cvg8MAj9gdNIsSZKmiZQSvHjo8tCSE0fz4tz50NpOcd+78oR5tSl/mtzKnkATY/xV4FfLrkOSJI2fNDwMh/blyfJILPbLL+XFBYvyZPmd78sNc4spf5paSm+gJUnS1JeGhuDA7tEbMirdcKaa8ndDM8WWN4+Gliwz5U9Tmw20JEkaszQ4AHt3jYaW9HbD+eqDfEuXU9x8B7R1UrR2mPKnaccGWpIkXVPqv1BN+asex9i1HUaulV2xmuJt9+XAkrZOihtuLLdYaYLZQEuSpB+Tzp+DXTtGj2TsuSTlb9U6invemafLre2m/Knu2EBLkqSc8lfZno9kVLpgby8MD+eUv7WbKB58zyUpf/PLLlcqlQ20JEl1KL3y8sXbMVLPNjiw52LKH+vaKB57f74hY+Nmitlzyy5XmlRsoCVJqgPppeOX38F8acrfxq0U7/lwbpjXt1E0zSq3WGmSs4GWJGmaSSnBsRerdzBXb8k4+kJenD0HNrVTvOOBfCRj7UZT/qQxsoGWJGmKyyl/B/OEeeSWjBPH8uK8BflBvwcezxPmVetN+ZPeIBtoSZKmmJzytzdPlkca5lMn8+LCxXmy/Fg1tGTFalP+pHFmAy1J0iSXhoZgf99oaEmlG86ezotLminabx4NLVnWYmiJNMFsoCVJmmTS4ADs6R29Uq53+2jK300rKG55RzW0pIOieVm5xUp1yAZakqSSpf4L0LdzNLRk987LU/7efv9ow7zYlD+pbDbQkiTVWDp/Fnp3jE6Yd1dgqJryt3p9Tvlr68wP/y1YVHa5kl7FBlqSpAmWzpzO9y+P3MG8b9flKX8PvzffkLFpK8VcU/6kyc4GWpKkcZZOnbwYWJJ6tsHBvaMpf+vbKN5VTfnbsIVi9pyyy5U0RjbQkiS9QenEsctDS144kBeaZsHGLRTv/XA+krG+jaKxqdxiJb1hNtCSJI3BxZS/S0NLRlL+5szNKX93PZSvlDPlT5qWbKAlSXoNKSV44ZKUv55tcPJ4Xpy/ADZ1UDz4OEVrJ6xeRzHDlD9purOBliTpEml4GA7uHT2/XOmCV17Oi4tuqN6OMZLyt8qUP6kO2UBLkupaGhqCfX2jV8pVuuDsmbx4400UHbdAW7VhvmmFKX+Sym+gQwiLgc8CnUACPh5j/Idyq5IkTVdpYAD2VC5J+dsBF0ZS/loobrkzx2K3dVDceFO5xUqalEpvoIFngL+KMb4/hNAEzC27IEnS9JEuXIC+HaN3MPfthIFqyl/LGop3PJAnzK0dFIuXlFuspCmh1AY6hLAQuBf4JwAxxn6gv8yaJElTWzp3FnZt55X9fQz98Luwp7ea8jcjp/zd91g+jrGpnWLBwrLLlTQFFSml0l48hHAz8CzQDbwFeA74ZIzxzKu+72ngaYAY4639/bXvsRsaGhgcHKz566q23Of64D5PL8OvnKK/+3kGup+nv+t5Bnf35JS/mTNp3LiFxo630tRxM41b3syMeab8TTf+PNeHMva5qakJ4IoPPZTdQN8G/D/grhjjt0MIzwCnYoz/9jV+WTp06FBtCrxEc3Mzx44dq/nrqrbc5/rgPk9t6dRL1evkLkn5A2hohA2b89nltk6ab7+T46fPvPY/TFOeP8/1oYx9bmlpgas00GWfgT4AHIgxfrv6+ZeAXy6xHknSJJNOHM1nl0ce+nvhYF6YNTun/N12dzXlr/WylL9i9hywgZY0AUptoGOML4QQ9ocQNscYdwIPkY9zSJLqUEoJjr5weWjJ8SN5cc482LSV4u5Hcsrfmo0UDWXPgSTVo8nwO8/PAV+o3sDRB3ys5HokSTWSUoLD+/OEuTKS8nciL85fmG/HePi9ecK8aq0pf5ImhdIb6Bjj88BtZdchSZp4aXgIDuytXimXp8ycPpUXFy2haOsYDS1ZsdrQEkmTUukNtCRp+kqDg7C/mvLX0wW93Zen/L3p1ouhJSw15U/S1GADLUkaN2lgAHb3jE6Yd+2AC+fz4rKVFLfeVQ0t6aS4cWm5xUrS62QDLUl63dKF89C3c/RKub6dMDiQF1eupbjzQWitxmIvuqHcYiVpnNhAS5KuWzp7BnZtzw1zpQv2VGBoKKf8rdlAcf+7KTZ35JS/+ab8SZqebKAlSVeVTp+CSvfohHn/bkg55Y91rRSPPkXR2pmvl5szt+xyJakmbKAlSRell1+6PLRkJOWvsSmn/D0e8gN/G7ZQzJpVbrGSVJLXbKBDCDNijMO1KkaSVFvp+NE8Wa7kaGxevDTlbyvF7ffkK+XWtVI0NpZbrCRNEteaQB8MIfwJ8PkY47ZaFCRJmhgpJThy+PKG+dKUv9Z2insezRPm1RtM+ZOkq7jW747/DPgI8N0Qwnbgc8B/jTEenfDKJElvSEoJDu0nVUZisbvg5WrK34JF0NpB8ciTecK8co0pf5J0nV6zgY4xfgX4SghhMfBB4KPAb4cQvkFupr8aYxyY+DIlSdeSU/72jIaWVLpHU/4WL6HY3DkaWrJ8laElkvQ6Xdffz8UYTwJ/APxBCGE9uZH+99WvNU9ceZKkq0mDg7C3txpa0gW92+FcNeWveRnFm2+vhpZ0wNLlNsySNE7GdMAthDALuB14G7AM+PuJKEqS9OPSQH9O+eu5JOWv/0JeXL6K4va787GMtg6KJab8SdJEua4GOoRwN/CPgQAcAf4E+Ocxxr0TWJsk1bV04Tzs2pGPZFS6oK9nNOVv1TqKux7OxzJa2ykWmvInSbVyrWvsfo18XGMJ8BfA4zHG/1uDuiSp7qSzZ6C3ezTlb2/v5Sl/Dz6ej2O0tlPMW1B2uZJUt641gX478CvAf48xnq9BPZJUN9Irp/J1cpWRlL891ZS/BljfSvHo+/IDfxtN+ZOkyeRat3A8VqtCJGm6SydP5MnyyC0Zh/blhZGUvydCvlJu/WZT/iRpEvOWfEmaIOn4EdLOkdCSbXDkcF6YNQc2baF42315wrzWlD9JmkpsoCVpHKSU4MVDl0+YT1Qzp+bOz+eW73uMorUzn2eeaWiJJE1VNtCS9Dqk4WE4vD9Plkce+nv5pby4YFG+f/nR91Fs7oCWtRQzZpRbsCRp3NhAS9J1SMNDsP+SlL/eLjj9Sl5cfCPF5jfD5o48YV6+0tASSZrGJkUDHUKYCXwPOBhjfKLseiTpYsrfxdCS7XDubF5cupziLXdAazUWu3mZDbMk1ZFJ0UADnwS2AwvLLkRSfUoD/dDXMxpa8mMpf/dejMUuljSXW6wkqVSlN9AhhFXA48BvAP+q5HIk1Yl0/hz07SDt7CJVtsHuHhgchKKAleso7n4kXynX2k6xcHHZ5UqSJpHSG2jg94BfAozVkjRh0tnTUNlOqmzjeN9OhnftgOFhmDED1mykePA9+TjGpnaKefPLLleSNImV2kCHEJ4AjsQYnwsh3P8a3/c08DRAjJHm5tr/9WlDQ0Mpr6vacp+nj+GXX6K/+4f0d/2Age7nGdrTCylBQyMz2tqZ976P0Nj5Vho3dzJjzryyy9UE8Oe5PrjP9WGy7XORUirtxUMIvwV8FBgEZpPPQH85xviR1/hl6dChQ7Uo7zLNzc0cO3as5q+r2nKfp6508ni+HWPklozD+/NCUxNs2ELRVn3gb30bS1tWus91wJ/n+uA+14cy9rmlpQXgik+IlzqBjjF+CvgUQHUC/a+v0TxLUg4tOX5k9A7mnm1w9IW8OHsObNpK8Y4HKFo7YN0migZT/iRJ42cynIGWpNeUU/4Ojk6YK11wojqJmDs/345x/7spNnfCqvWm/EmSJtSkaaBjjH8L/G3JZUiaBNLwMBzad3nK36mTeXHh4jxZfuwf5fcta0z5kyTV1KRpoCXVrzQ0BPv7cmhJpQsq3XCmmvJ3QzPF1rfAyBnmZab8SZLKZQMtqebS4ADs6SVVqueXe7fD+XN5celyipvflo9ltHXCjTfZMEuSJhUbaEkTLvVfgN09pJ3V88t9O6C/Py+uWE3xtvvyhLm1g+KGG8stVpKka7CBljTu0vmz0LtjdMK8uwJD1ZS/Veso7nlnPo7R2kGxYFHZ5UqSNCY20JLesHTmNPR25zPMPdtg367RlL+1mygefg9Faye0bqWYa8qfJGlqs4GWNGbp1EmodJNGQksO7qmm/DXA+jaKx95PsbkjB5jMnlN2uZIkjSsbaEnXlF46fnloyQsH8kJTE2zcSvGeD+cH/ta3UjTNKrdYSZImmA20pMuklODYi3myXKlOmC9L+WunuPOhfIZ57UZT/iRJdccGWqpzKSV44SCpsg12Vu9hfqma8jdvQX7Q74HH84R59TqKGab8SZLqmw20VGdyyt9e0s6u3DT3dMErL+fFhYtzozwSWrJitSl/kiS9ig20NM2loSHY10caOY5R6YKzZ/LikqUUHbfk0JLWDljWYmiJJEnXYAMtTTNpYAD2VkZDS3p3wIVqyt9NLRS33JmPZWzupLjxpnKLlSRpCrKBlqa4dOEC9I2ElnRB304YqKb8tayheMf91ZS/dorFpvxJkvRG2UBLU0xO+ds+Glqyp3c05W/1eor7HsvHMVo7KBYsLLtcSZKmHRtoaZJLZ165PLRkXx+kYZg5s5ry9978wN8mU/4kSaoFG2hpkkmnXqoGllSvlDu4t5ry1wgb2ije/f58U8aGzab8SZJUAhtoqWTpxNGLt2PklL+DeaFpFmzcQvHen8wT5vVtFI1N5RYrSZJsoKVaSinB0RfyZHnkloxjL+bFOXNzyt9dD+czzGs3UTT4IypJ0mTjn87SBMopfwdIO7dVJ8xdcPJ4XpxfTfl76Il8JGOVKX+SJE0FNtDSOErDw3Bgz+iVcpVLUv4W3VBN+eugaO2EFatM+ZMkaQqygZbegJzyt2v0Srne7tGUvxtvoui8JU+Z2zrhphWm/EmSNA2U2kCHEFYDnweWA8PAszHGZ8qsSXotaWAA9lRGr5TbtR0unM+Ly1ZS3HpXtWHuMOVPkqRpquwJ9CDwizHG74cQFgDPhRD+V4yxu+S6JOCSlL+RCXPfThgcyIsr11K848FLUv6WlFusJEmqiVIb6BjjYeBw9eNXQgjbgZWADbRKMXz2DOlHz+UJc6UL9lRgaAiKGTnl7/53V0NL2k35kySpTpU9gb4ohLAOeCvw7ZJLUR1Jp09BbzdpZw4tObq/D4arKX/rWikeeSo3zBu3UsydV3a5kiRpEihSSmXXQAhhPvC/gd+IMX75CutPA08DxBhv7e/vr3GF0NDQwODgYM1fV+Nr6KXjDHQ/T3/X8wx0/YDBfX15obGJxrYOZnfewsz2t9DU1mHK3zTmz3N9cJ/rg/tcH8rY56amJoArPv1fegMdQmgE/ifw9Rjj717HL0mHDh2a4Kp+XHNzM8eOHav56+qNScePkirbLkZj82I15W/W7DxVbuvIoSXr2ygaG93nOuE+1wf3uT64z/WhjH1uaWmBqzTQZd/CUQB/CGy/zuZZuqqc8nc4N8ojt2QcP5IX58yD1naKex7JV8qt3mDKnyRJel3K7iDuAj4K/CiE8Hz1a5+OMX6txJo0RaSU4PD+fDtGTz7DzMkTeXH+whxY8siTecK8aq0pf5IkaVyUfQvH33GV0bj0aml4KKf8jVwpV+mG06fy4uIlebJcvYOZFasNLZEkSROi7Am0dFVpcDCn/FW6SDu3Qe92OHdJyt+bbstT5rZOWLrchlmSJNWEDbQmjTTQD7sro3cw79oxmvK3fCXFbXdVQ0s6KG5cWm6xkiSpbtlAqzTpwnnYtSNPmHu2QV/P5Sl/dz6Uj2O0dlAsuqHcYiVJkqpsoFUz6ewZ2LWdtLM6Yd7bO5ryt2YDxQPvHm2Y5y0ou1xJkqQrsoHWhEmvnIJK1+iEef8eSMMwswHWbaJ49Kl8fnnjVoo5c8suV5Ik6brYQGvcpJMn8mR55JaMQ/vyQmMTbNhM8UTIV8pt2EIxa1a5xUqSJL1ONtB63dLxI5eHlhypJkSOpPzdcW+eMK9rpWhsLLdYSZKkcWIDreuSUoIjh0dDS3q2wYmjeXHuvHxu+d535oZ5zQaKmYaWSJKk6ckGWleUhoerKX9d+RxzTxe8XE35W7AoN8yPvi8/9LdyjSl/kiSpbthAC6im/O3fQ6psI+3sgt4uOP1KXlx8I8XmN1VDSzpg+SpDSyRJUt2yga5TaXAQ9vbmWOxKF/R2w7mzeXHpcoo33zGa8te8zIZZkiSpyga6TqSBfujryRPmnmrKX/+FvLh8FcXt91RT/toplpjyJ0mSdDU20NNUOn8O+nbkCXPPNtjdA4ODUBQ55e/uR6qhJe0UC035kyRJul420NNEOnsaereTRq6U27crp/zNmAFrNlI8+EQ+jrFpqyl/kiRJb4AN9BSVU/62jZ5h3r8bUsopf+tbqzdkdMKmLRSzTfmTJEkaLzbQU0Q6eXz0Srmd2+Dw/rxwMeXvQ/lIxobNFE2m/EmSJE0UG+hJKKUEx4+MhpZUuuDI4bw4aw60bqV4+/3VlL9NFA2m/EmSJNWKDfQkkFKCFw/lhrkykvJ3LC/OnZ8f9LvvsdwwrzblT5IkqUw20CVIw8NwaF+eLO/clt+fOpkXFyzK9y+/8yfykYyWtRQzZpRbsCRJki6yga6BNDQEB3aTRprlSjecqab83dBMsfUto6Ely1YaWiJJkjSJ2UBfh/TiIc5t/z5sveX6vn9wAPb0kipd+cG/3m44fy4vLl1OcfMd0NqZJ8ym/EmSJE0ppTfQIYTHgGeAmcBnY4yfKbmkH5O+8y1O/Y8vMuM/fPGKV8Kl/guwu2c0tKRvB/T358UVqynuuC9PmFs7KJY017h6SZIkjadSG+gQwkzg94FHgAPAd0MIX40xdpdZ16sV61rzg357d8HmN+WUv107RkNL9lya8reO4p53UrSOpPwtLrt8SZIkjaOyJ9B3AL0xxj6AEMKfAU8Ck6qBZt0mAIa/8oXcKO/theHhnPK3dhPFg+/JxzE2tVPMm19ysZIkSZpIZTfQK4H9l3x+AHjbq78phPA08DRAjJHm5hofg2hu5vjajQzurtDY2k7T+z5CY+dbadzcyYw582pbiyZUQ0ND7f/7Us25z/XBfa4P7nN9mGz7XHYDfaWn59KrvxBjfBZ4dmT92LFjE1rUldz425/l2PFjDDfN4jxwHuDMufymaaO5uZky/vtSbbnP9cF9rg/uc30oY59bWlquulb2BcMHgNWXfL4KOFRSLa+pmDXLiGxJkiSVPoH+LtAaQlgPHAQ+BPxkuSVJkiRJV1fqBDrGOAj8LPB1YHv+UuwqsyZJkiTptZQ9gSbG+DXga2XXIUmSJF2Pss9AS5IkSVOKDbQkSZI0BkVKP3Zr3GQ35QqWJEnSlHSlK5en5AS6KOMthPBcWa/tm/vsm/vsm/vsm/tcz28l7vMVTcUGWpIkSSqNDbQkSZI0BjbQ1+/Za3+LpgH3uT64z/XBfa4P7nN9mFT7PBUfIpQkSZJK4wRakiRJGoPSkwinghDCY8AzwEzgszHGz5RcksZZCGE18HlgOTAMPBtjfKbcqjQRQggzge8BB2OMT5Rdj8ZfCGEx8Fmgk3z16cdjjP9QblUabyGEXwA+Qd7jHwEfizGeL7cqvVEhhD8CngCOxBg7q19bAvw5sA7YA4QY40tl1QhOoK+p+oft7wPvAtqBD4cQ2sutShNgEPjFGONW4O3Av3Cfp61PAtvLLkIT6hngr2KMW4C34H5POyGElcDPA7dVm6yZwIfKrUrj5I+Bx171tV8GvhljbAW+Wf28VDbQ13YH0Btj7Isx9gN/BjxZck0aZzHGwzHG71c/foX8B+7KcqvSeAshrAIeJ08nNQ2FEBYC9wJ/CBBj7I8xniy3Kk2QBmBOCKEBmAscKrkejYMY47eAE6/68pPA56offw54qqZFXYEN9LWtBPZf8vkBbKymtRDCOuCtwLdLLkXj7/eAXyIf09H0tAE4CvyXEMIPQgifDSHMK7soja8Y40Hgd4B9wGHg5RjjN8qtShNoWYzxMOSBF3BTyfXYQF+HK6XQeHXJNBVCmA/8JfAvY4ynyq5H4yeEMHKm7rmya9GEagBuAf5TjPGtwBkmwV/3anyFEG4gTyXXAy3AvBDCR8qtSvXEBvraDgCrL/l8Ff410bQUQmgkN89fiDF+uex6NO7uAt4bQthDPor1YAjhT8stSRPgAHAgxjjyN0hfIjfUml4eBnbHGI/GGAeALwN3lmSROy8AAANRSURBVFyTJs6LIYQVANX3R0quxwb6OnwXaA0hrA8hNJEfUvhqyTVpnIUQCvKZye0xxt8tux6Nvxjjp2KMq2KM68g/x38TY3RiNc3EGF8A9ocQNle/9BDQXWJJmhj7gLeHEOZWf/9+CB8Wnc6+CvxU9eOfAr5SYi2A19hdU4xxMITws8DXyU/5/lGMsavksjT+7gI+CvwohPB89WufjjF+rcSaJL0+Pwd8oTr06AM+VnI9Gmcxxm+HEL4EfJ98i9IPmGRJdXp9QghfBO4HmkMIB4BfBT4DxBDCT5P/z9MHyqswM4lQkiRJGgOPcEiSJEljYAMtSZIkjYENtCRJkjQGNtCSJEnSGNhAS5IkSWNgAy1JkiSNgfdAS9IkF0L4NLAhxviJCX6d+4G/Ac4C748x/tU1vn8XOak1GkojqZ7YQEtSyUIIpy/5dC5wARiqfv4zMcbfrGE5h2KMq67nG2OMG0MIvwZsmtiSJGlysYGWpJLFGOePfBxC2AN8Isb41+VVJEl6LTbQkjTJjUx5Y4wfCSGsA3YDHwf+HTAf+BTwHPCHwBrgT2OMP3vJr/848G+A5cB3gKdjjHuv87U3Vf+5NwMDwDdjjB8cn38zSZqafIhQkqamtwGtwAeB3wN+BXgY6ABCCOE+8gdPAZ8GfgJYCvwf4ItjeJ1fB74B3ACsAv7jONUvSVOWDbQkTU2/HmM8H2P8BnAG+GKM8UiM8SC5SX5r9ft+BvitGOP2GOMg8JvAzSGEtdf5OgPAWqCl+np/N87/HpI05dhAS9LU9OIlH5+7wucj56rXAs+EEE6GEE4CJ4ACWHmdr/NL1e//Tgihq3ocRJLqmmegJWl62w/8RozxC6/nF8cYXwD+KUAI4W7gr0MI34ox9o5jjZI0pTiBlqTp7T8DnwohdACEEBaFED5wvb84hPCBEMLItXYvAYnRK/YkqS7ZQEvSNBZj/G/AbwN/FkI4BWwD3jWGf8TtwLerd1V/FfhkjHH3+FcqSVNHkVIquwZJ0iQQQrgX+Do5yOWDMcavX+P7d5LPUscYo2ejJdUNG2hJkiRpDDzCIUmSJI2BDbQkSZI0BjbQkiRJ0hjYQEuSJEljYAMtSZIkjYENtCRJkjQGNtCSJEnSGPx/noXpvmt9DakAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure()\n",
    "plot(t,V)\n",
    "xlabel('Time [s]');\n",
    "ylabel('V');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's put it all together . . ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: I&F CODE (version 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  In Parts 1 and 2, we constructed parts of the I&F model in bits-and-pieces.\n",
    "  Let's now collect all of this code, compute a numerical solution to\n",
    "  the I&F model, and plot the results (with appropriate axes).\n",
    "\n",
    "First, let's clear all the variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAADTCAYAAACoXCQbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXTdZ33n8fdPluR9iS3Htrwvkm1JQMhK4uw72WnhC7ShlGXSng6FduhQoGcGBgYGerqQ6XRmSAM0LA39srQwLYfQpqfQNUAgHCTZ1pXlRd4tL/Fubc/88fzkKztermzd+9PV/bzO0ZHufa7u75s8kfP1o+f3fJIQAiIiIiIiUpiqrAsQERERESknaqBFREREREZADbSIiIiIyAiogRYRERERGQE10CIiIiIiI6AGWkRERERkBKqzLuAS6Nw9ERERESmF5FxPlmMDzc6dO0t+zbq6Onp6ekp+XSktzXNl0DxXBs1zZdA8V4Ys5rm+vv68Y9rCISIiIiIyAiVZgTazzwMPAXvdvSV9bjbwl8AyYAtg7n6wFPWIiIiIiFyqUq1A/zlw/1nPfRB43t0bgOfTxyIiIiIiY1pJGmh3/wFw4KynHwWeSb9+BnisFLWIiIiIiFyOLPdAz3P3XQDp5yszrEVERERExohw8jih9UUGv/lFBp/7ZtblvEJZnMJhZk8ATwC4O3V1dSWvobq6OpPrSmlpniuD5rkyaJ4rg+Z5fBg8cpje9T+jr/0lett+ykBXBwwOwoQJTLzxjjE3z1k20HvMbIG77zKzBcDe873Q3Z8CnkofhiyOq9ExOZVB81wZNM+VQfNcGTTP5SkcPgi5dsLGVkKuDXZshRCgugZWNJI88CaShmZYuYb+iZPo7+8fU8fYZdlAfxt4O/Cp9PO3MqxFRERERIokHOghdLRCri1+3r0jDtROhFVrSa5ZR9LYDMsbSWpqsy22AKU6xu5Z4Hagzsy2Ax8hNs5uZu8CtgFvKkUtIiIiIlI8IQTYtzuuLHekDXPPnjg4eQqsaiJZdzdJYwssWUlSXRY7is9Qkord/a3nGbqrFNcXERERkeIIIcDu7YSONuhojZ8P7Y+D06ZDQzPJXQ/HFeZFy0iqJmRb8Cgov5ZfRERERDITBgdhx1bCULOca4MjL8fBmbNjo9zYTNLQAgsWkVSNv+BrNdAiIiIicl5hYAC2bSJ0tMVtGbk2OH4sDs65kqTl6rjKvLoF5i4gSZJsCy4BNdAiIiIiclro64MtufwK86YNcOpEHJy3kOSadekKczPJnMqM8VADLSIiIlLBwqlT0LUhv8LctRH6euPgwqUkN90BDS0kDU0ks2ZnW+wYoQZaREREpIKEE8ehcz0hl64wb+mEgX5IqmDxcpLbXh/3MTc0kUybkXW5Y5IaaBEREZFxLBw9DJ3thI3pCvO2Lggx5Y9lDST3PBqPlFu5hmTK1KzLLQtqoEVERETGkfDywfR0jHSFecfWOFBdAytWkzxocYV5xWqSiZOyLbZMqYEWERERKWPhwL4YVtLRFhvmPWnK38RJsHItyXW3xFjs5Y0kNTXZFjtOqIEWERERKRMx5W/XmaEl+/fGwclT477lW+6NK8yLV5Rlyl850L9VERERkTEqhAC7uvMrzLk2OHQgDk6bEY+Tu+fRuMK8aOm4SPkrB2qgRURERMaIMDgA29OUv1wbdLTB0cNxcNbseLPfUGjJ/EUVEVoyFqmBFhEREclI6O+PKX+5NsLGVuhcDyfSlL+6eSSvuhZWt8QV5rnz1TCPEWqgRUREREok9PXC5lx+hXnTBjh1Mg7OX0hy7TpobCFpbCaZPTfbYuW81ECLiIiIFEk4dRI2bYgrzB1pyl9/XxxcuJTkprviDX+NzSQzrsi2WCmYGmgRERGRURKOH4NN62MsdkcrbO2EgYGY8rdkBckdD6T7mJtIpk7Puly5RGqgRURERC5ROHoYcu1xS0ZHG3RvTlP+qmHZKpJ7H0tT/taSTJ6SdbkyStRAi4iIiBQoHDpw+nSMkBuW8ldTG5P9HrJ4w9+KNSQTJ2ZbrBSNGmgRERGR8wj7954ZWrJ3ZxyYOBlWrYkpf6tbYGmDUv4qiBpoEREREdLQkr27zgwtGUr5mzI1nr98631xS8aSFSQTFFpSqdRAi4iISEUKg4Owa3tsmIdOyXg5TfmbPjM2zPc8Fk/JWLiUpKoq24JlzFADLSIiIhUhDA5A9xZCrpWwsQ062+DokTg4aw7J6lfF4+QaW+KZzAotkfPIvIE2s98G3g0E4OfAO9z9ZLZViYiISLkL/f2wtTN/BnNnO5w4Hgfnzid59fX5hrlunhpmKVimDbSZLQTeCzS5+wkzc+AtwJ9nWZeIiIiUn5jy15E/Um7TBug9FQfnLyK57tbYMDc0KeVPLkvmK9DEGiabWR8wBdiZcT0iIiJSBsKpk5x66YcM/vhf4z7mzR3Q3w9JAguXkdx8T9y/3NBMMmNW1uXKOJJpA+3uO8zsD4BtwAnge+7+vSxrEhERkbEpHD8KncNS/rZt4tDAAFRVwZKVJHc+FLdjrFqrlD8pqiSEkNnFzewK4BvAm4FDwNeAr7v7l8963RPAEwDufk1vb2+pS6W6upr+/v6SX1dKS/NcGTTPlUHzXP4GDx+it/1n9LX9lN72l+jfnIMQoLqamlVN1DZfxaRXX0NVQxNVk6dmXa4UURY/z7W1tQDn3Bif9RaOu4HN7r4PwMy+CdwEnNFAu/tTwFPpw9DT01PSIgHq6urI4rpSWprnyqB5rgya5/ITDu2Pe5dzbYSNrbCrOw7U1sZkv4feErdkrFjNYO1ETgLThub52IlMa5fiyuLnub6+/rxjWTfQ24DXmdkU4haOu4AfZ1uSiIiIlELo2ZNP+cu1wd5dcWDiZGhYS/K62+OWjGWrSKqV8idjR9Z7oF8ws68DPwH6gZ+SX2kWERGRcSKEAHt2nhlacmBfHJwyDRqaSG57fVxhXqyUPxnbsl6Bxt0/Anwk6zpERERk9ITBQdi5La4sD8Viv3wwDk6fGVeW73tDbJjrlfIn5SXzBlpERETKXxgYgO2b8ydk5NrhWJryd0UdyZpX50NL5inlT8rbeRtoM/tYge/R5+4fH6V6REREpAyE/j7YuikfWtLZDifTG/nmzie56npobCFpaFbKn4w7F1qB/iDwlQLe442AGmgREZFxLPSeSlP+0u0Ym9bD0LGyCxaT3HBbDCxpbCG5Yk62xYoU2YUa6FPu/o6LvYGZPTaK9YiIiMgYEE6egE0b8lsytgxL+Vu0jOSW++LqckOTUv6k4lyogS70r4/zRqMQERERyU44fhRy6+OWjFwbbO2EwcGY8rd0FcmdDw9L+ZuWdbkimTpvA+3uvWb2G8Cz7n7wQq8rSmUiIiJSNOHIy6dPxwgdrbB9y+mUP5Y1ktz/xnhCxsrVJJOmZF2uyJhysVM43g38oZn9LfAM8B13Hyh+WSIiIjKawsH9Z57BPDzlb+VakoffGhvm5Y0ktROzLVZkjLtgA+3uV5tZC/A24H8DE83sL4AvuvtPSlGgiIiIjEwIAXr2pGcwp6dk7NsdBydNhlVNJDfeEbdkLF2plD+REbroOdDu3gr8rpl9ELgbeBz4vpltAZ5x9z8obokiIiJyITHlb0dcYR46JeNATxycOj3e6HfHg3GFedFypfyJXKaCg1TcPQB/B/ydmX0B+ALwaUANtIiISAnFlL+tcWV5qGE+fCgOzpgVV5bvT0NLFixWyp/IKCu4gTazRcTV518BFgLfIO6LFhERkSIKAwPQ3ZUPLcm1w/GjcXB2HUnTVfnQknn1Ci0RKbILNtBmNhX4RWLTfAvwA+CTwDfd/XjxyxMREak8ob8PtnTmj5TrXJ9P+btyAcnVN6ahJc0kdTpNVqTULrYCvQfoBr4I/Kq7by9+SSIiIpUl9J6Cro350JLNG89M+Xvd7fmGeZZS/kSydrEG+m53//eSVCIiIlIhwsnj0Lkhv8K8OQcDacrf4uUx5a+xJd78N31m1uWKyFkudozd6ebZzG4BXgtMO+s1nyxOaSIiIuNDOHY0nr88dAbztk1npvzd/Ug8IWPVWpIpSvkTGesKuonQzP4EMOCfgBPDhkIxihIRESln4fCh04EloaMVdmzNp/wtbyR5fZryt2INyaTJWZcrIiNU6Ckcvwy0uPvOYhYjIiJSjsKBnjNDS3antwzVToSVa0geeWvckrG8kaSmNttiReSyFdpAdwOnilmIiIhIOTid8jc8tGQo5W/ylJjyt+6ueKScUv5ExqVCG+h3AX9mZs8ST+Y4zd1/MOpViYiIjBEhBNg9LOWvoxUO7Y+D06bDqmaSOx8kaWiBxctIqpTyJzLeFdpAXwO8HriVV+6BXjLaRYmIiGQlDA7Cjq35/cu5NjjychyceUV6OsZQyt8ipfyJVKBCG+hPAg+7+98XsxgREZFSCwMDsK0rf6Rcrg2OH4uDc64kab4aGtOG+coFSvkTkYIb6GPEFMJRZ2azgKeBFuKK9jvd/d+KcS0REZHQ1wdbcsNS/jbAqaGUv3qSq2+KsdiNzSRzrsy2WBEZkwptoP8r8Bkz+xiwd/iAuw9eZg1PAt919zeaWS0w5TLfT0RE5LRw6hR0bcifwdy1EfrSlL/6JSQ33hFXmBuaSWbNzrZYESkLhTbQn08//9qw5xLiivEl3y1hZjOI+6p/FcDde4HeS30/ERGRcOI4bFrPke4uBn72I9jSmab8VcWUv9vuj9sxVjWRTJ+RdbkiUoYKbaCXF+n6K4B9wBfM7DXAi8D73P3Y8BeZ2RPAEwDuTl1dXZHKOb/q6upMriulpXmuDJrn8WXwyGF621+ir/0letteYmBzBwwOcnzCBGpWrqHmkbdQ23wVNWteTdVUpfyNN/p5rgxjbZ6TELILEzSza4F/B9a5+wtm9iRw2N3/ywW+LezcWfo8l7q6Onp6ekp+XSktzXNl0DyXt3D4YHqc3LCUP4DqGlixOu5dbmyh7rqb2H/02IXfTMqefp4rQxbzXF9fD3HHxSucdwXazD5+kUZ26HX/zd0/com1bQe2u/sL6eOvAx+8xPcSEZFxKBzYF/cuD930t3tHHJg4Kab8XXtzmvLXcEbKXzJpMqiBFpEiuNAWjt8ys89zns57mPcCl9RAu/tuM+s2s9XuvhG4C2i/lPcSEZHyF0KAfbvPDC3Zn967PnkqrFpLcvM9MeVvyUqS6kJ3IoqIjJ4L/ckzFejk4g30ycus4TeBr6QncHQB77jM9xMRkTIRQoBd3XGFOTeU8ncgDk6bEU/HuPuRuMK8aKlS/kRkTDhvA+3uJYlWcveXgGtLcS0REclWGByA7VvTI+XiKjNHD8fBmbNJGpvzoSULFiu0RETGJP3uS0REiib090N3mvLX0Qad7Wem/L3qmtOhJcxVyp+IlAc10CIiMmpCXx9s7sivMG/aAKfSnX7zFpJcsy4NLWkhmTM322JFRC6RGmgREblk4dRJ6NqYP1KuayP098XBhUtJbroTGtJY7JlXZFusiMgoUQMtIiIFC8ePwab1sWHOtcGWHAwMxJS/JStIbn+AZHVzTPmbppQ/ERmfCm6gzWwO8ACwwN1/38zqgSp331606kREJFPh6GHItedXmLs3QxiECRNgWQPJvY+RNLTE4+UmT8m6XBGRkiiogTaz24BvAD8G1gG/DzQAvwM8XLTqRESkpMLLB88MLRlK+aupjSl/D1q84W/FGpKJE7MtVkQkI4WuQH8GeLO7P29mB9PnXgCuL05ZIiJSCmH/vriynIvR2OwZnvK3luS6W+KRcssaSGpqsi1WRGSMKLSBXubuz6dfh/Rz7wi+X0REMhZCgL27zmyYh6f8NTSR3HJvXGFevEIpfyIi51Hon47tZnafuz837Lm7gZ8XoSYRERkFIQTY2U3IDcVit8HLacrf9JnQ0Exyz6NxhXnhEqX8iYgUqNAG+v3A35jZ3wKTzeyzxL3PjxatMhERGZGY8rclH1qSa8+n/M2aTbK6JR9aMn+RQktERC5RQQ20u/+7mb0G+GXg80A3cL1O4BARyU7o74etnWloSRt0rocTacpf3TySV1+XhpY0w9z5aphFREZJwRvc3H0H8fQNERHJQOjrjSl/HcNS/npPxcH5i0iuuzluy2hsJpmtlD8RkWIp9Bi7L5G/eXC4U8B24K/d/WejWZiISKULp07Cpg1xS0auDbo68il/i5aRrLs7bstoaCKZoZQ/EZFSKXQF+mXgbcC3ids3FhP3QH8VWAv8rpn9urt/sShViohUgHD8GHS251P+tnaemfJ354NxO0ZDE8nU6VmXKyJSsQptoBuBB9z9X4aeMLMbgY+5+z1mdj/xrGg10CIiBQpHDsfj5HJDKX9b0pS/aljeQHLvG+INfyuV8iciMpYU2kDfQAxOGe7H5INUngMWjVZRIiLjUTh0IK4sD52SsXNbHBhK+XvI4pFyy1cr5U9EZAwrtIF+CfiEmX3E3U+a2STgo8DQvuflwIEi1CciUrbC/r2EjUOhJa2wd1ccmDgZVq0hueG2uMK8VCl/IiLlpNAG+u3AXwCHzewAMJu4Av3L6fhs4DdGvzwRkfIQQoA9O89cYT6wLw5OmRb3Ld92P0lDS9zPPEGhJSIi5arQc6C3ADeZ2RJgAbDL3bcNG/9xccoTERmbwuAg7OqOK8tDN/29fDAOTp8Zz1++9w0kq5uhfilJVVW2BYuIyKgp+BxoAHffZmbdQGJmVelzg0WpTERkDAmDA9A9LOWvsw2OHomDs+aQrH41rG6OK8zzFyq0RERkHCv0HOh64E+BW4FZZw1f9u8hzWwCcUvIDnd/6HLfT0Tkcp1O+TsdWrIeThyPg3Pnk7zmemhIY7Hr5qlhFhGpIIWuQH8WOA7cBXyf2Eh/FPjOKNXxPmA9MGOU3k9EZERCXy90deRDS16R8nfr6VjsZHZdtsWKiEimCm2gbwKWuPsxMwvu/jMzexfwr8CfXU4BZrYIeBD4BPCfLue9REQKFU6egK4NhI1thFwrbO6A/n5IEli4jOTme+KRcg1NJDPO/sWbiIhUskIb6AGgP/36kJnNBQ4DC0ehhs8AHwAUqyUiRROOH4XcekKulf1dGxnctAEGB6GqCpasJLnz4bgdY1UTydRpWZcrIiJjWKEN9AvAA8BfEUNT/hI4Qdy3fMnM7CFgr7u/aGa3X+B1TwBPALg7dXWl//VpdXV1JteV0tI8jx+DLx+kt/1n9Lb9lL72lxjY0gkhQHUNVY1NTH3D49S0vJaa1S1UTZ6adblSBPp5rgya58ow1ua50Ab6bcDQGUy/BbyfuGL8x5d5/XXAI2b2ADAJmGFmX3b3x4e/yN2fAp5KH4aenp7LvOzI1dXVkcV1pbQ0z+UrHNofT8cYOiVjV3ccqK2FFWtIHn5rXGFe3sgV9Qvp6enhJMCxE/FDxh39PFcGzXNlyGKe6+vrzztWaAN9j7t/DcDdTwD/HcDM3gh8/VILc/cPAR9K3+t24HfObp5FRM4WQoD9e/NnMHe0wr7dcXDSZFi1luTGO0gammHZKpJqpfyJiMjoKbSB/hzwtXM8/xSX0UCLiBQipvztyK8w59rgQLoSMWVaPB3j9gdIVrfAouVK+RMRkaK6YANtZivSL6vMbDkw/KDTFRB/Azoa3P0fgX8crfcTkfIVBgdh57YzU/4OH4qDM2bFleX7fzF+rl+ilD8RESmpi61AdwKB2DhvOmtsN/EsaBGRyxIGBqC7K4aW5Nog1w7H0pS/K+pI1r4GGtPQknlK+RMRkWxdsIF29yoAM/u+u99WmpJEZLwL/X2wpZOQS/cvd66Hk+mNfHPnk1x1Q9yW0dgCc65UwywiImNKQXug1TyLyOUIvadgcwdhY7p/uWsD9PbGwQWLSW64La4wNzSTXDEn22JFREQu4rwNtJn9E3H7xgW5+62jWpGIlL1w8jh0bsivMG/OwUCa8rdoGckt98XtGA3NJNNnZl2uiIjIiFxoBfrpklUhImUtHDsKne1xD3NHK2zblE/5W7qK5O6HSRpaoGEtyRSl/ImISHk7bwPt7s+UshARKR/h8CHItROGQkt2bElT/qpheSPJ/W8kWd0cA0wmTc66XBERkVFV6DnQmNk7iImEC4EdwJfc/QvFKkxExo5wcP+ZoSW7t8eB2lpYuTZN+WuB5Q0ktROzLVZERKTICmqgzez3gF8B/hDYCiwFPmBm9e7+iSLWJyIlFkKAnj1xZTmXrjCfkfLXRHLTXXEP89KVSvkTEZGKU+gK9LuB291969ATZvYc8ANADbRIGQshwO4dhFwrbEzPYT6YpvxNnR5v9LvjwbjCvHgZSZVS/kREpLIV2kBPBfad9dx+QJsbRcpMTPnbStjYFpvmjjY48nIcnDErNspDoSULFivlT0RE5CyFNtDfBb5iZh8EthG3cHwCeK5YhYnI6AgDA7CtizC0HSPXBsePxcHZc0mar46hJQ3NMK9eoSUiIiIXUWgD/R7gfwE/A2qAPsCB9xapLhG5RKGvD7bm8qElnRvgVJryd2U9ydU3xW0Zq1tI5lyZbbEiIiJl6IINtJlVufugux8GfsXMfhWoA3rcfbAUBYrIhYVTp6BrKLSkDbo2Ql+a8le/hOTG29OUvyaSWUr5ExERuVwXW4HeYWZfIh5Z9/O0ad5bgrpE5Dxiyt/6fGjJls58yt/i5SS33R+3YzQ0k0yfkXW5IiIi487FGuhfBx4Hfmhm64FngL9w97NvKBSRIgnHjpwZWrKtC8IgTJiQpvw9Em/4W6WUPxERkVK4YAPt7t8CvmVms4A3E4NUPm1m3yM20992977ilylSOcLhg2lgSXqk3I6tacpfDaxoJHngjfGkjBWrlfInIiKSgYJuInT3Q8Bngc+a2XJiI/3H6XN1xStPZPwLB/adPh0jpvztiAO1E2HlGpJHfimuMC9vJKmpzbZYERERKTzKG8DMJgLXATcA84B/LUZRIuNVCAH27Y4ry0OnZPTsiYOTp8SUv3V3xz3MS1eRVI/oR1RERERKoNAo75uJUd5GvInwS8BvDE8mFJFXiil/2wkbW9MV5jY4tD8OTktT/u56KG7JWKSUPxERkXJwsWPsPkrcrjEb+BrwoLv/SwnqEilLYXAQtm/JHymXG5byN/OKNOWvmaShBRYsUsqfiIhIGbrYCvTrgN8D/trdT5agHpGyElP+NuWPlOtsz6f8zbmSpOXquMrc2AJXLlDKn4iIyDhwsVM47i/mxc1sMfBFYD4wCDzl7k8W85oilyP09cGWXP5IuU3r4VT6d8t5C0muWZc2zM1K+RMRERmnsr5DqR94v7v/xMymAy+a2d+5e3vGdYkAw1L+hlaYuzZCf3py48KlJDfeOSzlb3a2xYqIiEhJZNpAu/suYFf69ZE0rGUhoAZaMjF4/Bjh5y/GFeZcG2zJwcAAJFUx5e/2B9LQkial/ImIiFSorFegTzOzZcBrgRcyLkUqSDh6GDrbCRtjaMm+7i4YTFP+ljWQ3PNYbJhXriWZMjXrckVERGQMSEIIWdeAmU0Dvg98wt2/eY7xJ4AnANz9mt7e3hJXCNXV1fT395f8ujK6Bg7up6/9JXrbXqKv7af0b+uKAzW11DQ2M6nlaiY0vYbaxmal/I1j+nmuDJrnyqB5rgxZzHNtbS3AOe/+z7yBNrMa4G+A59z9jwr4lrBz584iV/VKdXV19PT0lPy6cnnC/n2EXOvpaGz2pCl/EyfFVeXG5hhasryRpKZG81whNM+VQfNcGTTPlSGLea6vr4fzNNCZbuEwswT4HLC+wOZZ5Lxiyt+u2CgPnZKxf28cnDwVGppIbrknHim3eIVS/kREROSSZN1BrCMGtfzczF5Kn/uwu38nw5qkTIQQYFd3PB2jI+5h5tCBODhtRgwsuefRuMK8aKlS/kRERGRUZH0Kxz9znqVxkbOFwYGY8jd0pFyuHY4ejoOzZseV5fQMZhYsVmiJiIiIFEXWK9Ai5xX6+2PKX66NsLEVOtfDiWEpf6+6Nq4yN7bA3PlqmEVERKQk1EDLmBH6emFzLn8G86YN+ZS/+QtJrl2XhpY0k8yZm22xIiIiUrHUQEtmwqmTsGlDXGHuaIWujjNT/m66K27HaGgmmXlFtsWKiIiIpNRAS8mE48dg03rCxnSFeWtnPuVvyQqSOx7IN8xTp2ddroiIiMg5qYGWoglHDkOuLb/C3L0FwiBMqIZlq0jufSzuX165lmTylKzLFRERESmIGmgZNeHQgbiyPHRKxs5tcaCmFlasJnnI4pFyK9aQTJyYbbEiIiIil0gNtFyysH/vmaEle9OEyKGUv+tvjSvMyxpIamqyLVZERERklKiBloKEEGDvrnxoSUcrHNgXB6dMjfuWb70vNsxLVpBMUGiJiIiIjE9qoOWcwuBgmvLXFvcxd7TBy2nK3/SZsWG+9w3xpr+FS5TyJyIiIhVDDbQAacpf9xZCrpWwsQ062+DokTg4aw7J6leloSXNMH+RQktERESkYqmBrlChvx+2dsZY7FwbdLbDieNxcO58kldfn0/5q5unhllEREQkpQa6QoS+XujqiCvMHWnKX++pODh/Ecl1t6Qpf00ks5XyJyIiInI+aqDHqXDyBHRtiCvMHa2wuQP6+yFJYsrfzfekoSVNJDOU8iciIiJSKDXQ40Q4fhQ61xOGjpTbtimm/FVVwZKVJHc+FLdjrFqrlD8RERGRy6AGukzFlL/W/B7m7s0QQkz5W96QnpDRAqvWkExSyp+IiIjIaFEDXSbCof35I+U2tsKu7jhwOuXvLXFLxorVJLVK+RMREREpFjXQY1AIAfbvzYeW5Npg7644OHEyNKwled3tacrfKpJqpfyJiIiIlIoa6DEghAB7dsaGOTeU8tcTB6dMizf63XZ/bJgXK+VPREREJEtqoDMQBgdh57a4sryxNX4+fCgOTp8Zz1++7xfiloz6pSRVVdkWLCIiIiKnqYEugTAwANs3E4aa5Vw7HEtT/q6oI1n7mnxoybyFCi0RERERGcPUQBcg7NnJifU/gbVXF/b6/j7Y0knItcUb/zrb4eSJODh3PslV10NDS1xhVsqfiIiISFnJvIE2s/uBJ4EJwNPu/qmMS3qF8OFlDmIAAAeGSURBVMMfcPj/PUvV/3z2nEfChd5TsLkjH1rStQF6e+PggsUk198WV5gbmklm15W4ehEREREZTZk20GY2AfhT4B5gO/AjM/u2u7dnWdfZkmUN8Ua/rZtg9atiyt+mDfnQki3DU/6WkdxyH0nDUMrfrKzLFxEREZFRlPUK9PVAp7t3AZjZV4FHgTHVQLNsFQCD3/pKbJS3dsLgYEz5W7qK5M6H43aMVU0kU6dlXKyIiIiIFFPWDfRCoHvY4+3ADWe/yMyeAJ4AcHfq6kq8DaKujv1LV9K/OUdNQxO1b3icmpbXUrO6harJU0tbixRVdXV16f/7kpLTPFcGzXNl0DxXhrE2z1k30Oe6ey6c/YS7PwU8NTTe09NT1KLOZc6nn6Znfw+DtRM5CZwEOHYifsi4UVdXRxb/fUlpaZ4rg+a5MmieK0MW81xfX3/esawPGN4OLB72eBGwM6NaLiiZOFER2SIiIiKS+Qr0j4AGM1sO7ADeAvxStiWJiIiIiJxfpivQ7t4PvAd4Dlgfn/K2LGsSEREREbmQrFegcffvAN/Jug4RERERkUJkvQdaRERERKSsqIEWERERERmBJIRXnBo31pVdwSIiIiJSls515HJZrkAnWXyY2YtZXVsfmmd9aJ71oXnWh+a5kj8ynOdzKscGWkREREQkM2qgRURERERGQA104Z66+EtkHNA8VwbNc2XQPFcGzXNlGFPzXI43EYqIiIiIZEYr0CIiIiIiI5B5EmE5MLP7gSeBCcDT7v6pjEuSUWZmi4EvAvOBQeApd38y26qkGMxsAvBjYIe7P5R1PTL6zGwW8DTQQjz69J3u/m/ZViWjzcx+G3g3cY5/DrzD3U9mW5VcLjP7PPAQsNfdW9LnZgN/CSwDtgDm7gezqhG0An1R6f9s/xR4PdAEvNXMmrKtSoqgH3i/u68FXgf8R83zuPU+YH3WRUhRPQl8193XAK9B8z3umNlC4L3AtWmTNQF4S7ZVySj5c+D+s577IPC8uzcAz6ePM6UG+uKuBzrdvcvde4GvAo9mXJOMMnff5e4/Sb8+Qvwf7sJsq5LRZmaLgAeJq5MyDpnZDOBW4HMA7t7r7oeyrUqKpBqYbGbVwBRgZ8b1yChw9x8AB856+lHgmfTrZ4DHSlrUOaiBvriFQPewx9tRYzWumdky4LXACxmXIqPvM8AHiNt0ZHxaAewDvmBmPzWzp81satZFyehy9x3AHwDbgF3Ay+7+vWyrkiKa5+67IC54AVdmXI8a6AKcK4VGR5eMU2Y2DfgG8FvufjjremT0mNnQnroXs65FiqoauBr4P+7+WuAYY+DXvTK6zOwK4qrkcqAemGpmj2dblVQSNdAXtx1YPOzxIvRronHJzGqIzfNX3P2bWdcjo24d8IiZbSFuxbrTzL6cbUlSBNuB7e4+9BukrxMbahlf7gY2u/s+d+8DvgnclHFNUjx7zGwBQPp5b8b1qIEuwI+ABjNbbma1xJsUvp1xTTLKzCwh7plc7+5/lHU9Mvrc/UPuvsjdlxF/jv/B3bViNc64+26g28xWp0/dBbRnWJIUxzbgdWY2Jf3z+y50s+h49m3g7enXbwe+lWEtgI6xuyh37zez9wDPEe/y/by7t2Vcloy+dcDbgJ+b2Uvpcx929+9kWJOIXJrfBL6SLnp0Ae/IuB4ZZe7+gpl9HfgJ8RSlnzLGkurk0pjZs8DtQJ2ZbQc+AnwKcDN7F/EvT2/KrsJISYQiIiIiIiOgLRwiIiIiIiOgBlpEREREZATUQIuIiIiIjIAaaBERERGREVADLSIiIiIyAmqgRURERERGQOdAi4iMcWb2YWCFu7+7yNe5HfgH4DjwRnf/7kVev4mY1OoKpRGRSqIGWkQkY2Z2dNjDKcApYCB9/Gvu/skSlrPT3RcV8kJ3X2lmHwVWFbckEZGxRQ20iEjG3H3a0NdmtgV4t7v/fXYViYjIhaiBFhEZ44ZWed39cTNbBmwG3gl8DJgGfAh4EfgcsAT4sru/Z9j3vxP4z8B84IfAE+6+tcBrr0rf9yqgD3je3d88Ov9kIiLlSTcRioiUpxuABuDNwGeA3wPuBpoBM7PbiF88BnwY+AVgLvBPwLMjuM7Hge8BVwCLgD8ZpfpFRMqWGmgRkfL0cXc/6e7fA44Bz7r7XnffQWySX5u+7teA/+Hu6929H/gkcJWZLS3wOn3AUqA+vd4/j/I/h4hI2VEDLSJSnvYM+/rEOR4P7ateCjxpZofM7BBwAEiAhQVe5wPp639oZm3pdhARkYqmPdAiIuNbN/AJd//KpXyzu+8G/gOAmd0M/L2Z/cDdO0exRhGRsqIVaBGR8e3/Ah8ys2YAM5tpZm8q9JvN7E1mNnSs3UEgkD9iT0SkIqmBFhEZx9z9r4BPA181s8NAK/D6EbzFdcAL6VnV3wbe5+6bR79SEZHykYQQsq5BRETGADO7FXiOGOTyZnd/7iKv30jcS+3urr3RIlIx1ECLiIiIiIyAtnCIiIiIiIyAGmgRERERkRFQAy0iIiIiMgJqoEVERERERkANtIiIiIjICKiBFhEREREZATXQIiIiIiIj8P8Bl3g6ff8VEZUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pylab import *\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize']=(12,3)# Change the default figure size\n",
    "\n",
    "I=1                              #Set the parameter I.\n",
    "C=1                              #Set the parameter C.\n",
    "dt=0.01                          #Set the timestep.\n",
    "V = zeros([1000,1])              #Initialize V.\n",
    "V[0]=0.2;                        #Set the initial value of V.\n",
    "\n",
    "for k in range(1,999):           #March forward in time,\n",
    "    V[k+1] = V[k] + dt*(I/C)     #... updating V along the way.\n",
    "\n",
    "t = arange(0,len(V))*dt          #Define the time axis.\n",
    "\n",
    "figure()                         #Plot the results.\n",
    "plot(t,V)\n",
    "xlabel('Time [s]')\n",
    "ylabel('Voltage [mV]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:**  Adjust the parameter `I`.  What happens to `V` if `I=0`?  Can you set `I` so that `V` > 20 within 10 s?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4:  Voltage threshold\n",
    "\n",
    "  Notice, our model is missing something important:  **the reset**.\n",
    "  Without\n",
    "  the reset, the voltage increases forever (if $I>0$). Now, let's update\n",
    "  our model to include the reset.  To do so, we'll need to add two things\n",
    "  to our code.\n",
    "  \n",
    "  - First, we'll define the voltage threshold `Vth`, and\n",
    "  reset voltage `Vreset`.\n",
    "  - Second, we'll check to see if `V` exceeds\n",
    "  `Vth` using an [if-statement](https://docs.python.org/3/tutorial/controlflow.html);  if it does, then we'll set `V` equal to\n",
    "  `Vreset`.\n",
    "  \n",
    "Here's what we'll add to the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Vth = 1;        #Define the voltage threshold.\n",
    "Vreset = 0;     #Define the reset voltage.\n",
    "\n",
    "for k in range(1,999):            #March forward in time,\n",
    "    V[k+1] = V[k] + dt*(I/C)      #Update the voltage,\n",
    "    if V[k+1] > Vth:              #... and check if the voltage exceeds the threshold.\n",
    "        V[k+1] = Vreset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 5: I&F CODE (version 2)\n",
    "\n",
    "Now, let's put it all together to make a complete I&F model (with a thershold and reset), simulate it, and plot the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAADTCAYAAABDa58YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9d5Qd132g+d1GzqmZABAgSAAkQQAEEYgcmIMkUhKlIkVREmXJtKyRwxl7Zjwzux6v03q8O2csz9oecz2yZa/X3rI8Y8tj9gMBkP0a6G7kSGQQBJEJdM7d6K67f1S96tf1XndXv1fhVr36ztEh+oWq+1PVu/d3Q91PSClJSEhISEhISEhISBg5ZWEXICEhISEhISEhISGqJMl0QkJCQkJCQkJCQoEkyXRCQkJCQkJCQkJCgSTJdEJCQkJCQkJCQkKBJMl0QkJCQkJCQkJCQoEkyXRCQkJCQkJCQkJCgYwOuwBFkuzrl5CQkJCQkJCQEBTC+ULUk2muX78eynnLy8upq6sL5dwJwZBc49Iguc6lQXKdS4PkOpcGYV3n2bNn5309kGRa07QfAZ8Hbum6vjTP+wL4IfAy0AG8rev64SDKlpCQkJCQkJCQkFAoQa2Z/gvgxSHefwlYZP3vHeBPAihTQkJCQkJCQkJCQlEEkkzrul4FNAzxkVeBv9R1Xeq6vheYrmnafUGULSEhISEhISEhIaFQVNnNYw5wJevvq9ZrCQkJCQkJCQkJCcqiygOIOU9GMshOHZqmvYO5FARd1ykvL/ezXIMyevTo0M4dNEZrCx3/rDPpK28jRqtyy/hLX0MdHT/5MbO+/A1EmSp9Tn/p++w6XTUfMOlLb4VdlMDovfIJXRVpyl96LeyiBMadC6e588l5Jj73SthFCYyeU0fpOXOU8k3Phl2UwOg+shfZ0834tVvDLkpgdO1N03t9IuXL14RdlMDorHqfUeV3M3bJirCLEiiq5WCqZEZXgfuz/p4L5N2mQ9f1d4F3rT9lWE/tltITw8b/+H+Q7+l0Ll6GmPdQ2MUJBOOv/yuy8j3al65GlN8TdnECwfjRHyL3V9GxajNi4qSwixMIfe/+Jzh5hLY1pZNw9P3xf4RbN+h4YkPYRQkEKSXGH/4OowTwSGkkHLK3F+OHvw3TZ9L20GNhFycQZHcXxh/+FmMXP0bf7AVhFycQZGsLxn/5bcSKdZTdPTfs4gSCbG0Gw+CuhxYptZuHKkNuPwW+qWma0DRtHdCs6/qNsAuVYFbKcs/71h+lsa237OpE7v3Q+qNEYm5pQh6qsf4okZhv34STR0omXgB5+WP45BxII+yiBMfZE/DZtZK6zhzfD80NJRWz3F8FnR2lFXPNLujtRZbQ71nu/CnGr30Xo60l7KIMIKit8f4G2AaUa5p2FfgPwBgAXdf/K/Ae5rZ4FzC3xvt2EOVKcMGxfdDSZP67RCopub8KujqtP0ok5ppd0Neb+SvUsgSF3L3dvr5SSoTIt9osXsj0dusf4ZYjSGQ6Zf5XyrzrCeOIYcVcShc6c51Lps42DGRVJuZwyxIUsvcOcs8OeOwJyiZPhS51VgcEkkzruv61Yd6XwL8IoiwJI8OorOj/owQqKSklMl1iMRsGsmp7/wtGCcR85w5yz86sFwwQo8IrUADIzg7kvkrrj9IYyZLNjcgjtdYfJRLzZ9fh1FHzjxL4LQPIS+fh0wvmv0ugzgbgzDG4ZU3gl8i9zVFzcK9s60thlyQHVZZ5JCiIvHkVzhyHhY9aL5RAJXXpPFy+WFoxnzoKt2/2x1wCwxzySC20Nmdd53DLEwRyXyV0d5kxl0C8gDmK1dcHDz1SOjFXbYeyMliwmFIJWqZTMHYczF1QGnU21uzD5Klwz5xSucxmzLPuhsfUe/YhSaYTBkVWbYdRoxCbnrdeiP8vVqYrYNx4xLqnMq+EWp4gMNIpmDINsdJ6IK0krnMK7roXseQJ64V4xyylRFZWwLwHESWSZEmjD7n7fXhkOeLeObG/xgDyTg+yZiesWAfTZpZGzB1tyP1ViLVbYeKk0oi5qR6O7kNsfAbGjKEkfs/W4J7Y8gKiTL1ZxCSZTsiL7OlG1nyAeGI9Ytr0sIsTCLK9DXlgN2LtNhg/wXox1CL5jmyog+P7EZuehdFjMq+GWia/kTeuwLmPEFteMEfwzFdDLZPvXDwL1z5FbH0JhCiJhIOTR6D+FmXbXsLcfTX+MctDNdDWStm2l/JvOBtD5N5K6OlGbH3RvLdL4Trv2QGGgdjyIlAav2eZ3g6jRpttlYIkyXRCXuShGmhvNSuoTK0c8x+s3Psh9PQgtr5gVcoQ94pZ7tkBUiI2v9Df+MY7ZHNUetRoxMasSjn2MVfA+AmIJ7dQKomlkU7BtBnw+NqS6UDIdArung0PLyuJmM1nXFLwwCLE/IXWi+GWyW9kXx+y6n1Y8gTi7vtKotNkDu7tQqxcj5g6I+zi5CVJphPyItMV5lqsTKVsvhpqmfzErpQXLDb30hbx70DIPmsa/LEnEHfdSyl0mmR3N7L2A8SqDYgp00rj3m5vRR7Yg1j3FGL8hJLIpWX9bTh+ELHxOVM0JUT8Y772KVw4hdj6giWain8yzYXTcP2yNehDSXQg+OgQNNZRVkIxy4PV0NHWf50VJEmmE3KQVz+Bj88gtr5obhcW/3wDzp2EG1fMaXCgP7EMrUT+c/wANNUPrJSBOActD+6GjvaBjS/EujGS1bug94454wKUQpIld28HJGLL89mvhlWcQJDpChg9BrH+GfOFUhixrKyACZMQazabL5TAMg+jsgKmzwTb8lgCv+eqFNw7BxYvDbsog5Ik0wk5yHTKrJQ3PG29UgIJR7oCJk5CrN4EZOWVMa6YjXQFTJ8Fy6xKOf65tHlv33c/LMpY4eLdabJnXB56BDHXssLFPOEwRVM7YOkqxKy7zRdjPnonuzqRtR8iVm9ETJkKgIh5kiVbm5GHqxHrn0KMG5/1Roxjvn0TTh5GbHrenHGB7MYqlsgrjsE9RUmS6YQByK4OZG0lYs0mxKQp5osxH7GULU3Iw7WI9U8jxo0zX4z5iGXG/ic2P48YlXkyOuYxf2ra/wZUynHvNJ05DreuZ824EP9lHsf2Q3OjYy/amCeWB3ZDV+fAafC4d5os+5/5EJ5F3DtNu98HBGLzc443YhxzugLGjEWsf3r4D4dIkkwnDEDur4LuzoGNr/1m8OUJAllt2v+Es/E13w2jSL6T2YtWbM6aBo97p6kqBWPHItY/1f9i3DtN6RRMmoJYvTHr1XgnWUa6AmbeBctW9r8oiLXMQ6ZTMGc+PPRo/4sxXicuDcOMefFjiDnzst5Rd+SyWGz73+NrEDPv6n8jxp0m2dWB3JtGrM4a3FOUJJlOsLGnhOcugAcf7n8jxkmWrWR9eBnivrn9b8Q3ZNP+V70THn8SMWNWng8EXya/Me1/acSaLYiJk7Peie+Fls2NyKN7ERufQYwZ2/9GjEfv5GfX4fSx3L1oFZ4eLpaM/U9sfSl3Gjym15kzx+D2zYGj0hDrTpM8sg9am/ufcckQ59/zPmtwb5t6xkMnSTKd0I9l/8tZmxTn0btTR6HuszxPCcc35oz9L0fJGudOk2X/y5lxiW/Itv1PbH7B8U6MRywzoqmNzr1oY5xwpFOWaGrbwDdiPGKZsf/ZoqkMcU4s0xWm/S8jmhrwZvDl8RtbNHX/AsvmqTZJMp1gIysrYNwExLqtjnfim1ga6QrT/vfEuoFvxDmxTFfAXffCo48PfCOmnaZ++99D8MBCx7sxjdnoMxPLRx837X/ZCIjlfW3b/9Yips8c+KYgdtcYsOx/acSTWxATJjrejWdiKRsz9r9nEWPGON6NZwdC3rgKZ0+YA11ljrQtrp2mT87B1U/yz7goSJJMJwBk2f+2IsY7KmX17+OCkA11cOwAYtOziNGOSjmmubS8fhnOnbSmwZ0//3gmlnx8xrL/5XkaPK6dpo8OQ8Pt3ClhiO3onTxUbdr/8j3vEdckq7bSEk3lu86BFycQ+u1/zhkXYrtOXFZlRFPP5L4Z199zOmUO7q3dEnZRXJEk0wkAyNoP4E5P1l602cQzyZJ73gdknmlwiGs2bU6Dj84zDU58G990Ksv+5yCel3mg/S+HeF7oAfY/JzFMssxnXCoG2v+yiWGSZYumMva/vB+KWcwRsP95jSma2o1Yl2dwT1GSZDoh1/7nJIajd7n2PwcxXPIgu7uRNVn2vxxiGHNbC/Jglv0vhxje2/W34ESW/c+JfZljFPPVS3DhdJb9z0EME0vOn7JEU4NZ4WIY84mDA+1/TmK45EEe3DNQNOUkhvd2/+Ce+g8eZkiS6QTT/nfz6uBPzMYv37Dsfw2DTAlDHGUe8uBu6GwfvIKKY6ep5gPL/jdU40usGiN7L9p80+BALDtNVRnRVJ5p8P5PBVaeIJDplGX/G2QaPIYTEEY6Zdr/Hn8y/wfimFimU3Dv3CHsf/GKeYBo6v4FYRfHNUkynZBj/8slfo2vka6AGeWwbHX+D8RQ5tFv/1uS/wMxy6XtSnnho4i5DwzyqXh1mmz737JViFl35f9QzDpN/fa/TYjJU/N/KGZJlm3/25AlmsohZjFn7H8DRFP5PhhcmfxGXr4IF88Obf+LW6fp3Edw81rutoeKkyTTJY5saTTtfxueQYwdpFKOW+N760Ye+58DYf00YtIY9dv/hngyOmYx99v/hqiUy+J1b3Nsn2X/GyJmEbMOxP6qXPufk5itmZbVOy3732CzD1j3dnyClru3AwKx6fnBPyTKiFXM6dTw9j9RFp86GyvmiZMdoin1SZLpEse2/w3VC4zZVLht/9v03OAfsnOsmMScrrDsf9tcfDgeMRvpCtP+t2qoSjlzbxuBlMlvjHTKtP8tXTn4h0R8YrYfwpszHx56ZPAPCoGMQbyAJZrabtr/Zs8b4pMCjHj8lk37307L/lc++AcFYMTkOndlRFObEZMmD/5BIWLxWwaXg3uKkiTTJYytZHXa/5zEKJke1v5nE6OYO9qR+6vy2P8cxGgGQjY1WHvROux/TuITMvLmtfz2v0G/4H+ZfOfSBUs0NdxetDFa8nB6EPufkxg9jCeP7M1v/3MiRGwerJV706b9b7iYIS6X2eww9fUOsquY2iTJdClz6gjU33LxY43PtLA8XANtLUM8eGgRp8RyXzq//c9JnDpN1TtN+9+w6+5iFPNuy/431IwLxOveTlfkt/85iZGoxhZNOe1/TmK0tEWmU1B+T3773wDisYDYft7Djf0vJp0mWzT18DLEvUMM7ilKnn2T/EHTtBeBHwKjgD/Tdf33HO/PA34MTLc+82u6rr8XVPlKESOdym//cxKjh/FkVSq//c9JTPINexp8/kLEgkXDfDoenaYB9r97Zg/94ZgklvJOD7J6F2LFOsS0YfaijUmnSXa0IQ9Umdse5tj/nMQjsZSN9XBsP+L5L+Wx/+X9hu9l8hvb/vflb+Xf9jCbuJguL5417X/f+P7w9r+4PFx78ijU36LsK2+HXZKCCGRkWtO0UcAfAS8BS4CvaZrm3FLgfwF0XdefAN4A/jiIspUq/fa/53Ltf07i0vhm7H/5lKw5xCPmbPvfcIi4dJps+5+LPUrjkUub9r/2VndTwnHpNA1l/3MSk4RjSPufk7jEPJT9z4GIyXKeIUVTMcVIV8DU6bAin2hKfYJa5vEkcEHX9Yu6rvcAfwu86viMBDL7Gk0DrgdUtpLE3ItWIjYP8WS0TTwyDplOwejRw+xFaxGXEct0BUyYiFizefgPx6TTZFRWWPa/QfaiHUA8YpaVFXDPHHhk+fAfjkGnyZ5xGUw05SQGyzxknzXjMphoKofoJ5amaCpj/5s+/BdisORBtrdaoqlt7ux/Meg0yYbbcPygu8E9RQkqmZ4DXMn6+6r1Wja/AbyladpV4D3gF4IpWukhe3st+99Kd5VyDPJK2d1l7kW7cuMg9r/BvuhfmfxGtrYgD1ZblXI++5+T6F9oWX8LPjpkVcouVrHFoNMkr34CH58xHzwcbkoY4tFpOn9yGPufk+gnHJw4AE317mZcIBbLh/vtf25jjv5yHllj2f9c77Mc/Xtb7t6B+8E9NQlqzXS+n7Xz6n8N+Atd1/+Tpmnrgb/SNG2prusD9nzRNO0d4B0AXdcpLx9imxwfGT16dGjnLpauvWmamxuY9v1/w3gXMdxpbaABmDJliqvPq0jnzv9JS2c70199nbEuYuiZPp1GYNq0qa4+ryLt1e/T1nuHGa9+jTEuYuiaOpVmYPq06a4+ryJt2/+ediGY9eobjHIRQ+eUKbQAM2fMcPV5FWn5+7+gc8xYyj//VcqmDt9RbJ80mTZg1qyZlE2Y5H8BfaD5Lz+ke+Jk7nrxi4hx44f9fNvEibRLIltnAzTWfEDvzHLKn34RMWr4prtlwgS6hIh0zA01OzHmzGfWhq2uOopN48bRi4xszFJK6vfsoOzhpcx8Yo2r7zSOG4txp5tZUY25t5e66p2MWbmOGY885vp7quVgQSXTV4H7s/6eS+4yju8ALwLoul6radp4oBy4lf0hXdffBd61/pR1dXW+FHg4ysvLCevcxdL3TzrMKKd1/mLaXMQgm5oBaG1pdvV5Fen757+D2fNovmsOwk3MLS0ANDc1ufq8akjDwHjvv8PCR2meNA3cxNzaCkBTYyNiSgRj7u3FeP8fYekqGsVoVzEbbW0ANDQ0IERgz2N7huzqxKisQKzaSEPPHXcxd7QDUF9Xj5jQ6XcRPUe2NGHUfIjY9hL1rW3Q2jbsd4zOLkBGts6Wt29iHN2H+Pzr1Dc2ufqO0dWNNPqiG/PljzHOnUS8/l3q6+tdfcfo6aFMRvg6nzmOcf0yxrd/2XUMfXfuwJ070Y35cA1GYx296783ohjCysFmz87/UHtQrccBYJGmaQuAa5gPGL7p+Mxl4BngLzRNexQYD9wOqHwlg7x1A04dQbz65tBK1mwiPhMuP70Al84jvvaOu2lwIPJrac+eMO1/X3jD/XeivuTh2D5oaaJsm8spYYj8kgfb/jeSmCN+b8saSzTleokHkd/lQe7eDkIgNo9g/92ILx+W6e2WaGoI+18O0V7yINMpUzQ1IvtftGM2RVPlsGxV2EUpikDWTOu63gv8ANgOnDZf0k9qmvabmqa9Yn3sV4Cf1TTtGPA3wNu6rkf3DlEUV/a/HKKdZMl0CsaOQ6x7yv2XIr7e0EhXwOQpiFXD7EWbTbQvs1kpz7obHhtuL9psoptY2g/hzX0AHnzY/Rcj3Gnqt/8tRdx3//BfyBDhh7Rs+9/y4URTTiIcs1v7n5MId5pkSyPySC1iw9Mjs/+5HiBSD3nrOpw6itjsUjSlMIHNa1p7Rr/neO3Xs/59CoiWjD1i2Pa/FWsR00dQKUd49E52tJuV8pNbEBNHsj40wjFn7H/PvDK0/S+HCCdZGfvfF98aWaUc3XYILp037X9f/94IZlyI8mWGU0dN+98X3xrhF6N7oeXhWnf2PycR3tlC7q207H8jmXEh2onl7h2WaKoA+18E2ykodHBPTRIDYgnRb/8bYaVsNURR/L3KfZXQ0z2yKWGI9ujdngIr5Sh3mqpS7ux/OUQ45oz9b+22EX4zuvd2v2hq/ci+aF/m6MXcb/9bMbIvRnRnCymludXjvAfhgeFEU06iqROXRp+5w9Yjy0du/4top6l/cG8dYvrMsItTNEkyXULIdIVp/3tkGPufk4juS2srWecvRIy0Uo5ovmFWyi7tfzlEU+Yhe7qRNR+4s/85iWinSba3IQ/sRqzd6sL+5yCinSbZUAfH9yM2PuvS/pdNRGO+cQXOfYTY4kY0lfcInpfJdy6etUVTI5pxgegu8zh5xLT/jXigi8guYZKHqqGttbCYFWTQZR6apv2my2Pc0XX9tzwqT4JPyGuX4fwpxFfeHnmlHNHGl49Pm5XyN39QwJcjGvOJw9BQR5n23ZF/N6qdpkM1I7D/OYhmLo3c+6F7+18OEe007dkBUhY2DR7VTlPavf0vh6gmWekKy/63tYBvR3OU1kinirP/RS9k896++z53oqkIMFRW9WuY29kN979/7XMZEzxAVmXsf88W8O1oJpYynTLtf4UoWSPa+BrpCpg206X9z0FEO00yPQL7Xw7Ri9mecXFr/3MSwU6T7LOmwV3b/xxE8Ocsu7uRtR8gVm1wZ//LIXrLPGR7K/LAHsS6p1yKphxEcGmLrM/Y/54vzP4XwWUe8tqncOGUOftQ0IyLegz1AGK3ruvfHu4AmqZ90cPyJPiAbf9btRExZerwX3ASwWc6bPvf5uddSR0GP5B3ZfIb2/73Oc2d/S+HCCaWGfuf9p2RTwlDNDtNGfvf279U2Pej2GnK2P++/nMFHiB6Mffb/wqcBo+gQl3WfAC9dxBbC5h9gEgu85B73gckYkuh9r/ozUDIdApGj0GsL2DGRVGG6hK43e7hHi8KkuAfcn8VdI5AyZpDBBuiml1WpVxoQxS9JEtWbQcEYlOBlXIUO02ZSnnDSPaizXcgb8oTBLKyAiZOQqzeVOARondvG5UVMH0WLHNnhcshir/ndAXcdz8scm+FG0DElnnYMy4PPYKYu6DAo0RrlFb29pq7eCxdhZh1d2EHiVi9Lbu7kHs/RKwucHBPUQZNpnVd79E07fuapg35RI+u6z3eFyvBS2Q6BbPnwcJHCztAWbQaInMv2hQsXIKYM7+wgwjrpxGRxsjci3YHLF+NmHVXYQeJWsxdHcjaSsSaTYhJUwo7iD3FGJGYW5qQh2sR659GjBvBXrTZ2CPT3pXLT+Ttm3DyiDnL5FY05SRqMX/6MXxyrrCH8DJEbfr/zHH47FoRgz6Yv+eI1F8AHNsPzQ2UFROziFbM5uBeR+EDXYoy3GKV7wLXNU37iaZpX9A0Ldq7apcg8tJ5+PRCcZVyputrROQHe+Y43LpR3I/V7j9EJOajlv2vhGKW+6sK24t2ANG6t2V1AfY/J3ZiaXhTKJ+x96LdXOg0OBGMOWXa/0YimspBROa+BvrtfyMRTeU9UHRiNtIVlv1vZeEHidjSFplOwZz58FCBg3uKMmQyrev6SmAN8DHwx8ANTdP+QNO0Iq58QpDIqu0jt/85idgUqZFOjdz+5yRi60oLs/85iU7M9pTwSO1/Dvr7lxGIOTPjMlL7n5MI/ZztvWhHbP9zEp2RadlZoP3PSYRGpmVzI/Lo3pHb/5xEaGmL/Oy6KZoq0v4nIrRm2pvBPTUZ9jFKXdc/0nX93wDzgK8DM4C0pmknNE37Vb8LmFA4tv1v7dYR2v8GO2Dxh/Ab2VQPR/dae9GOxP7nJDoZh7x5Fc4cR2wpUskapU5Txv639aXiKuUodZpOHYW6zxDbihmJh0h1mo5Y9r9iY45Sp2lfJXR3Iba+XNyBopNLZ4mmip36j1BiWbW9QNGUgyh1mtIpUzRV1IyLmrh+5F/XdQnsAHZomvbnwJ8D/xH4P30qW0KRFGz/cxKhJEvu2QmGUdhetNlEJ+SsSrmQbQ/zHdCbw/iJaf+bgFhXyF602UTnQvfb/9YVd6Ao/Z7TKVM09egIRVNOItJp6rf/PQQPLCzyaNFIsmz736OPI+6dU9zBhPn/oepjnvJOD7JmJ6xYW7z9LyLbAcqONuT+qsJEUxHAdTKtadpc4C3gm8Ac4O+BH/tUroQisSvl+QsR84uslKPSEGXsf0tWIO4eqf3PSURi7ulGVu9CPLEeMXWE9j8nEUmybPvfuqcR44uslKMRsmn/O7Yf8eKXCtuLNu9BvTmMX8jrl03732vf8mAv2ohc6I/PmKKpb/yL4qfBo7Lk4aPDpv3vq8PuxDs8EVk60G//K3aWKXNA9a+z3FvpzeCeogyZTGuaNgl4DTOB3gxUAb8L/Hdd1zv8L15CwVw4DdcvF2j/cxKNxJITh0z73+sF2P+cRCWxPFgNHW3eVFBR6TTVfmDZ/4qcfQCicm/be9Fu9iDmqNzbVdst+58HMy7RCNkciR8/oTDRVA7RGLE00imYNgMeL9D+N4BodCBM+99seHhZ8QeLwDIP+xmXBxYVP7inKMONTH8GXAH+Enhb1/Wr/hcpwQtkuqJw+5+TaHT2rUp5JiwvwP7nJCqNb1XKtP95USlHILEs2v7nJAKJZdH2PycR6DTJ7m5kjWX/mzLNgyNGIOa2FuTBPYhNzxVm/3MSAWmLrL8FJw4iXvpqgaIpBxHY2UJevQQXTiO++m2P7H8R6ECcP2UO7n3rF8IuiW8Md/c+q+v63kBKkuAZsrUFeajafEq4GPufTQQaorrPirT/OYlAknWlSPufkyh0ms6dhJtXC7f/OVH/MsPxA9DUQNnXv+fRASPwez642xJNeTQlHIVOk23/8zBmha8xYHYSoQj7n5MIjNJWZURTHtn/IlBvy3QKJkxCrNkcdlF8Y8isIzuR1jRtM/AEMNnxmd/1p2gJhWLa/3qL3H83iyg0RLvfB0Rxe9FmE4XRu6oUjBlbvP3PJgIxpy3735pC7X9O1L+3jXQFzCgv3P7nJCqNbzH2PyeKX2Z7xmXho4i5D3h0VLUvtOztNXfxWLa6cPufE8UfxpNdncjaDxGrNyEme2T/U7zTJFubkYerEVte9GhwT01cDeFpmvZfAA3YDXRmvaXuFSxR7L1oFy1BzJnnzUFVb4iy7X8zC7T/5WAGLaWaTZJt/1tdhP3PieKdJtv+99TLxe1Fm43inSbb/vfKm4Xb/3JQPOaM/e+Ndzzci1bte9sUTV1HfOF1745pX2ap5p6+x/ZDc2NxoiknqieWB3ZDV6fHD+EpHrM9uBfPBw8zuJ0P/zqwVNf1634WJsEDzhwz7X+vvOnhQRVvfI9k7H8ejcSD8vvSyn0Z+5+XDVHm4N4d0ktk9U7T/lf0XrTZqC3zkOmUaf8rdi/abFTvNKUrTPvf+m3eHVTxTpORrrDsfxs9PGpWzAom06b97y5Y6rUDTs1rbO+wNWc+PPSIdwdW79LaSMMw67DFjyFmezS4pyhuV79fAbr9LEiCN5j2v6mIlUUqWbOJQuM76254bIV3B1W48bUr5bkLirL/5aJwzJlK+eFliPvmendghTtNtv3v8WLtfw4U/jnLjnZzL9o1WxATi7D/5aBup3j9+9QAACAASURBVEk2NcDRfYiNzxQpmnKgcL0tb14z7X/FiqacqDwyfekCXP64eNFUDgrHfOYY3L7p8QCImrgdmf4O8H9rmvY3mDt82Oi6XuV5qRIKwrT/7UM89ypijEd70YLaje+Nq3D2BOJL3/C2UlY4seSTc3D1E8TXf97bSlnhxpdTR6D+FuK1b3l7XJU7TYdroK3F2xkXQOV7W+5LW/Y/j2NWudNUvdO0/3mx7WE2Cv+c5W5LNOXFtocDUHfNtCmaGo9Yt83bAyu8NZ5RWeH94J6iuE2mVwEvAVvIXTMd77H7CCH37PDG/peDurWy5/a/DCpPnaVTHtn/BjuBP4ctBs/sfzkonFhWeWT/c6Jop8l8CM8r+58DRTtN0ugz6zAv7H85KBrznR5k9S5v7H9OFN0OUHa0IQ9UIdZu897+p2guLRvrTdHUc1/0dnBPUdwm078LfEHX9Z2FnkjTtBeBHwKjgD/Tdf338nxGA34D89Y4puu6lwt/Y43s60NWve+R/c+Bqg1RTzeyxiP7Xw6Kxpyx/633wP7nRNUkq6EOjh3w1v6XQdFOk2n/O+mR/W+wk/hz2ILx0v6Xg5r3Nh8dhobblGk/4/2xVf09H6qGdg/tf9kIgVSszgaQtZWWaMqHmBVd5uHf4J6auK2l2zHthwWhadoo4I8wR7eXAF/TNG2J4zOLgH8LbNR1/THglws9X0ny0SForPOnglJ0vaFt/9vmT6VsncX7YxeBrP0A7vT4G7NiFbOn9r8cFI25ajuM9sj+50TVezud8k405UTNkD22/w2CYjHLtCWaemS5D0dXL7G0Z1wWLEbM90A05UTBZR4DRFN33xd2cQLB7cj0rwN/oGnabwK3st/Qdd1w8f0ngQu6rl8E0DTtb4FXgVNZn/lZ4I90XW+0jnsr5ygJg2JUVsD0mbDco71os1F0vaFMV8C9c2DxUu8PrmDja+9F++DDiPsX+HAG9TpNsrfXqpRXemP/c6JgYim7u0z738qNHtn/HCjYabLtf5s9sv/loGDMGfvfyx7Z/5yoeG/b9r+f8We7PhWXeZw/CTeuIN7+RX+Or+Le2icOmoN7X3sn7JIEhttf8I+s//5c1muZ7pCbp77mYO4IkuEq4OyKLwbQNK3aOuZv6Lqeclm+ksbci/Yw4nOv+1spq9QQXb4IF88iXvfI/peDejFz7iPT/vdtj+x/TlTsNNn2v5/35/jq5RvmXrRe2v9yULDTVLPLtP/59dS/iollVUY05dM0uIr1djpj//NKNOVEvcTStv+t9tP+p1bQRjrl3+CeorjNvIodBsuX7Tiv/mhgEbANmAvs1jRtqa7rTdkf0jTtHeAdAF3XKS8vL7JohTF69OjQzu2kNfUTOoRg1quvM8qHMhnjx3IbmDRpEpMUibnlJ39O59ixlH/uq5RN8cgklUVvbzf1wJTJk5mgSMxNP/6AnklTuOuFLyLGeSQtyeJOcx0NwJQpUxivSMyNtbvonXU35U+9gBjlfUexe9p0moBp06YyVpGY66t3Iu9fwKz1W3zpKHZNnUozMGP6NEYrELM0DOr37KTskeXMXLHal3N0TJlMKzBzxkxGzVIg5t5e6mp2MmbVemY8/Kgv52ifNIk2oHzWLCXMc0ZnB3X7Khm/6RmmPfCgL+donTiRDimVaZuNpgZuH65hwgtfYuocrx8wNWmZMIFuIZSJue+z69SdPMwk7dtMvuce386jUg4GLpNpXdc/LfI8V4H7s/6eCzgFMFeBvbqu3wE+0TTtLGZyfcBRlneBd60/ZV1dXZFFK4zy8nLCOnc2svcOxvv/CMtW08go8KFMsqMdgPa2NjpViLmrA6MyhVi1iYbuHuj2IeYmsw/X2tpKuwoxtzRi1KYRT71MfWsrtLZ6f47mZgBaW5ppUyHmWzcwju5HvPIm9Y1Nw3+hkHO0tADQ3NSEUCHmTy9gXDiNeOMd6uvr/TlHWxsAjY2NiEkKxHz6GMaNKxgvf8W3OtVo7wCgob4eIcN/6lQeqsZoaqB3/dP+xdxpbrxVV1enRjJdtR3Z2UHP2qf8i7mrC5BKtM0ARurvobeX7ie3+hdzdzeyr0+dmP/xbwBB58pNdPlYprBysNmz82/wMOgDiJqm/ZabA2ua9r+5+NgBYJGmaQs0TRsLvAH81PGZfwCeso5Zjrns46KbMpQy8sg+aG326cFDC3u60M3yeP/xxf7nRLWYq3f5YP9zYsVsqDFlKKu2e2//c6LYVLhMp2DsOG/tf07sFQ9qxGykK2Cy1/a/QVAm5pQpmvLc/peNOnWY/RCe1/Y/J0KoU38ZhlmH+W3/U0hUI3vvIPfshMfXIGaqM2ocBEONTP+ypmk/YvjNo34R+A9DfUDX9V5N034AbMdcD/0jXddPWg80HtR1/afWe89rmnYK6AP+la7r/gzNxIh++98T/p1EoeWGtv3vfq/tf07USbJ8s/85UWhdqW3/W7HWW/ufE4WSadv+96TX9j8nCsWcsf8984q39j8nKt3bGfvfF9/yWDTlQJ2QLfvfRcTXv+fTMy4ZFNrZ4rRl//viWz6fSJ2Y5ZG9/g/uKcpQyfQk4ALDJ9Ndbk6k6/p7wHuO1349698S+JfW/xJcYNv/vvxN//aiBZSqlTP2v7e+72+lHP5McD+2/e9tf8+j0GXut//5raFVKLG07X8+x6xSYpmx//m9F61KnaaM/c/PGRdAqXs7Y/9bu83fEwmUiBesGZcp0xBPrPf3ROrk0uagT/k9sGRF2EUJnEGTaV3X/czOEjxAVqVg1GjExmf8PZFKDVFlhWn/W+vDXrQDUCdmozJTKfu4Fy2gUjYt0xWm/e8Rj+1/ThTpNNnT4PMXIh5Y5O/JFLnMA+x/93gsmspBjd+z7OlGVu9CrFiHmOa1aMqBIp0mUzRVhVj3lPf2PyeKLHmw7X/PfykA+58iMduDez6KphSm9CKOCbb9b6Uf9j8namylJdtbzb1o12313v7nRJWGqOE2HD+I2PSc9/Y/J4p0muS1y3D+FGLriwFUymrEzMenTfuf7yPxoEw2fcKy/wUxJaxKp+lQDbS3BnOdFbnMcu+Hlv0vyHs7XOTu94Oz/ykibZHpimAG9xQlSaYjijy4Bzr83Is2C0X2H7btf0E2vmE3RLt3ABKx+fkAzmYGHXZeKatSpv1vQwCVsiqdpoz9b42fe9FaKNJpMtIVlv3vyQDOpkbMMl3ho/3PSfj3ti2aWrAYMc8H+58T+zKHGHO2/c8P0ZQTBaQtsrsbWfuBNbg3PdzChESSTEcUmU7BvXP9sf85UaDx9d/+50SBmP22/zlRoNMku7uQtR/6Z/9zEn6+gWxtQR6sRqzb5pP9z0n4M02y/hZ8dMiacfFBNOVEgU6TvPoJfHwGseUFnx/Cs1Cg3rbtf4GMSoMK9TYnDkBTfcAP4YXcSbQH90rvwcMMSTIdQWz739YXg6mUVaigzn0EN68F92NVoPHl+AFobqBsW8Axh9mByNj/gopZgXtb1lr2v8DubfvMwZwvD77b/5wo8HOW6e2m/S+waXAFOk3pFEz02/6XhQL1tmn/mxWc/U+BdeIyXQH33Q+LHwu1HGHiekhA07RZwMvAfbqu/76mabOBMl3Xr/pWuoS8yHQKxoxFrPdLyepAgWVoZqU8GbE6gL1oQYnG10inYGY5LFsV0BkVSCzTKZg9Dxb6Y4XLIeTG19z2cDssXIKYMz+Yk4bcaZK9vcg978Py1YhZdwV01pBj7upE7v0QsWYTYtKUYE4acqdJtjQhD9Ugtr3ki7E1LyHX2/L2TTh5BPH5NxCjfNz2cADhLvOQlz+GT84h3vjZgAb31MTVyLSmaVuBs8DXgf/VenkR8Cc+lSthEGRXB3JfGrFmM2KSn3vRZhNyQ9TSiDxci9jwDGJsQJVyyLWyvHUdTh1BbH7e371oswm5HpSfXoBL5wOccck+ebCnszlzHG5dD3AaHML+PXN0L7Q0BbDtYRZhd5r2p6GrM9hp8LA7TRnRVAnd27JqOwgR0DMuFgJCnWVKb4exYxHrngqtDCrgdpnHHwCv67r+ItBrvbYPCOLJkYQs5N60//Y/J2E3RHt2WpVyQFPCEH5DFIT9L4eQY87Y/4KslEO+t410yrL/bQjupCF3mmz7n5+iqcEI4TLbz3vMfcBn0ZST8O5t0/6XgsVLEffdH9yJQ/w926Kp5U/6K5pyEuIyD9nZgdxXGfDgnpq4TaYf0HV9l/XvzFXrYQTLRBKKx66U718ACxYHd+IQ8w17L9qHlyHu9dH+l0N46w0H2P+mB1wpmyUI7pyZM3a0mzMuT25BTJwU3IlD7DTJpno4uhex4Vl/7X85hBjzzatw5niwMy4Qbqfp0nnT/hf0jEuY/cRTR6Hus4BHpbMIo94+UmvZ/4KOObyt8eS+Sks0VboPHmZwm0yf0jTNOSz4LHDC4/IkDMXFs6b9b+tLAU+DhzhiefKoaf8L+sca4npD0/7XGnylHGanaV8l9HSH0PiG2GnaszO4vWizCbPTVBWU/c9BmJ2moOx/OYQXs5FOBWP/cxLmvZ1OmaKpoO1/IeXSUkpTojbvQfBbNBUB3CbTvwL8taZpPwYmaJr2p8BfAP/Kr4Il5CLTKRgfhP3PQYgVlJGugKnTA7D/OQi78b37Pv/tfzmEE7M94xKE/c9JSJ0mafSZWuklKwKw/zkI6edsiqY+QDyx3n/7Xw7hdJpM+99uxNpt/tv/nIRUb8uGOji+H7Hp2QDsfw5CqrfljStw7iNz28PA7X8hLfO4eNYSTQU9uKcmrq66rut7gceBk8CPgE+AJ3VdP+Bj2RKy6Lf/bfPf/ufA/qEE3RAFaf/LIaSGKGP/2xKE/c9BWJ2mQO1/DsLqNJ04DA11IUwJQ2idpiDtfw5EWJ0m2/4X8OwDhNdp2rMDpAxu28MBhNRWpVOW/e/ZYE8MoRkQZbrCHNx7MuDBPUVxveZZ1/VrwO/7WJaEIZA1lv1vS0hr0MxSBHu23e9DYPY/B2E1ROmK4Ox/gxYi4NNVVpj2v1Aq5XAutGn/mwnLQ3iGO6wRy4z97+FlgZ4XCKXTFLj9L4cQYg7a/uckhE6Tbf9btSEY0ZSTEHJp2d6KPLDHHOgKRDSlPq6SaU3T/or8l6sbuAr8g67rx7wsWEI/dqX80CMB2f/yEPATw6b9bwcsXYUovyew8/YTQkPU3WXuRbtqI2LK1MDOaxNCkiVbW5CHqhGbX0CMGx/YeW1CyCtl3Wem/e9zWjD2v0ELEuCprlj2P+07IU0Jh3Chz1n2v7d/KbhzZhNGp+m4Zf/7+s8Fd85swug0Hdxt2f/CGugKPpuW1RnRVBizD2ridh65GXgV86pdtf77CtAHPArUapr2TV9KmABnT8Bn10IelQ54Y/jj+037X1gVVBiJ5f4q6OwI78noMBqiml3QG/RetNmEEPNuy/63KYQZFwjn3q5Kmfa/DQGJppyE0WlKV1j2v03BnTRvQYI7lZGuMO1/ywKy/+UQwr2dTpn2v0Uh2f+CHujKHtybG9LgnoK4TaYXAy/ruv4NXdf/na7r3wBeAh7Sdf0N4MvAv/OrkKVO4Pa/fJQF2/sN3v7nIIR14oHb/5wEnEzbe9EGaf9zUhZs4yt775hrSgO1/zkQVrUf1HXu6kDWVgZr/3MSdMwtTaZoav3Twdn/nNjPXAQUc8b+t/n5AO1/DgKut+Wnlv0vDNFUBhHwQJctmkq2w8vGbTK9FlPSks1B+qUt24EgNwEuGWRLI/JILWJjkPa/fAiQRiBnMu1/R82p/yD3os3GrpQDivnSefj0gqneDe3J6IBHac8ch1s3ENvCrJQDjvnoPsv+F2LMdv8hoCRrf5UlmiqhmG37nwL3thFQHZYRTYXxjEuGoOvtqpRp/1sfov1PBNc2gzXoM2lKuIN7CuI2mT4K/I6maeMBrP/+FpBZJ70AaPC+eAly9w7o6wt+L1onAQ5Mh2P/cxLwCEfVdtP+F/hetFkEPP1v2v+mIlYGaP9zEvBofL/9L+C9aAcQXMz99r8FAdv/nAR3b9szLg8vQ9wX4hhTgH1yWzT1eMD2PycBVmGm/S+NWLMFMTFM+19wF1o2NyKP7jUH9wIVTamP22T6W8BmoEXTtJtAC7DFeh1gJvB974tX2ph70b4PjywP2P6Xj2DWZck7PZb9bx1i+kzfzzcoAT4Vbtv/1m4N1v7nJMiGKGP/2/hM8HvRDsAMOohc2rb/bQlxxgWC7TR9ci4c+5+TIDtNp46Ea/+zCbDTdLjGsv+FPfUfYKdpb6Ua9j/7MgcQ8x5rcC+UbQ/VxtWj5LquXwI2aJo2D7gPuKHr+uWs9w/6U7wS5+QRqL9F2VfeDrskge1lKQ+FZP9zEmDja+5FG4b9z0mAMe/ZEY79z0mQnaZ0xv4Xwl602QTZaUqnYNwExLqt/p9sSIKbaeq3/63z/2RDEWCnSVZZ9r9HgxZNOQio3jZnXCpg3kPwwEJfzzU8WTH72GGVRp85g/ro44h75/h2nqgyIiuElUDvB65qmlamaVrQqp+SwkinYOp0WBGw/S8fAS3zkOmUZf9b7v/JhiTIStmy/80PuVIOqPGVfX3IqvdN+9/dAdv/nATV+PZ0I2t2mfa/qUHb/5wEFLNt/9sauGgqh4A6TbKhDo4dMO1/gYumHASUS8vrl+HcyZDsf04C6jR9fMYWTYVu/wuq0/TRYWi4Hf5Al6K43Wd6NvBHmEs7pjvedjVfqWnai8APrc//ma7rvzfI574C/B2wppRHvGW9Zf976bXwK2UgiGxaXvsULpxCfOXb4VfKQdWPF07D9cuIb/4goBMOQVADWR8dgsY6yt74rs8nckNAieXBauhoU2D2geA6TbWWaEqFvWiD6jTtyYimFIg5qHu7ant49j8nQXWa0il17H8B1dtGOgXTZsDjCgzuKYjbjOVPgR7gGaANWAn8FPiemy9rmjYKMxl/CVgCfE3TtCV5PjcF+EVydw4pOfor5RCfjM4mgL0sZTqzF22I9j+bgBqidAVMmKRGpRxQrWxWyiHZ/5wE1GmSVSm4NyT732D4eJnDt/858f/3HLr9z0kAnSbZ3Y2sCdH+5ySATpNsa0Ee3INY95Qi9r8AYq6/BScOIjY+F65oSmHcJtMbgJ/Rdf0oIC3b4XeAX3H5/SeBC7quX9R1vQf4W0wJjJPfwlSWd7k8biwJ3/6XD3/3srTtf6tDsv85CaIhytj/1j8Vjv3PSRANUcb+t/l5RSrlAGLO2P9UmBKGYEamz52Em1dD3vYwiyD+bz9+AJoaFHgIz4Gf9fbB3dDZHv5DeDYB1Ns1H1j2PwVmmSCYtiojmgr7GReFcdua9QG91r+bNE27C3NHD7er0OcAV7L+voq5d7WNpmlPAPfruv4/NU371cEOpGnaO8A7ALquU15e7rII3jJ69Gjfzt1V+yHNzQ1M/8K/ZVxI8Tm5VSaYMGE8U3wqT8eOn9La2cH0V95grAIxy74+bgETJ0xksk/lad+9nbbeXma++gajFYjZGDua28DkSZOY6FN5WlM/oUMIZr36OqMUiLm3p5N6YMrkyUzwqTwtP/lzOseOpfxzX6VMgY7inYZbNABTp0zxrX5p+vEueiZN4a4XvhietCSL7mnTaAKmTZvmW/3SWLOT3ll3U/7UC4hR4XcUu6ZOpRmYMX26b/VL/Z4dyPsXMGv9FiU6ih1TJtMKzJwxk1E+bNEnDYP6PTsoe2Q5M1es9vz4hdA+aRJtQPmsWb5sVyd7e6mr2cWYVeuZ8XBIQrE8+JmDFYLbX/w+4GXgf2AKWv4/oBNT3OKGfL8yuxtlPcj4n4G3hzuQruvvAu9mjlFXV+eyCN5SXl6OX+fu+ycdZpbTMn8hIqT4nEigs6ODbr9i/uefwJz5NJffp0TM0hIddLS30+VDeaRhYFT8PSxaQtPEqaBCzG0tALS1tdHhR8y9dzDe/0dYtppGRqkRc1MjAK2trbT7EXNXB0ZlCrFqEw3dPdCtQMzNzQC0NLf48luTLY0YtWnEUy9T39oKra2en2PkZTLL0NzU6E/Mt25gHN2PeOVN6hubPD9+IRhtbQA0NjYgJni/D7L89ALGhdOIN96hvr7e8+MXgtHWDkBDfT2iz/uRWnn6GMaNKxgvf8W39n+kGB0dANTV1fmTTB+qxmisp3f908rEDP7mYEMxe3b+h+bdLvP4BpC2/v3LwAfAR8DXXH7/KnB/1t9zgetZf08BlgKVmqZdAtYBP9U0TY2uX4DIz67D6WPh2v/y4uOWOxn7nyrT4OD/1NmZY3D7pkLTo/i+zEMe2afIXrTZ+Bzzvoz9T5EpYfD9IS3b/rdFwZh9+jmrIZpy4u/OFjKdsf9t8+cEheBzvW2kK0z73yqF7H8+19tGOgUz74KlK305flxwOzL9nK7rfweg63on8Ntg77zxExffPwAs0jRtAXANeAN4M/OmruvNgD1er2laJfCrpbibh5qVMtZmHj41vulU+PY/B3ZS71NDpIT9Lwd/GyKZrlDA/ufAx8ZXSomsrFDA/ufAx8ZXGob5ew7b/peDjzGrYv9zIOzHXHyIuaMdub9KAfufAx+rMNnUAEf3IZ75gmL2P//aKnnzmjm498W3FBvcUw+3I9P/bZDX3x3k9QHout4L/ABzichp8yX9pKZpv6lp2isuyxB75J0eZI0C9r+8+LObh+xoMyvlsO1/g+JDzE31ZqUcuv3PgZ8N0Y2rcPZE+Pa/wfCj//DJObj6iVozLoCvF/qUKZpSaiQe/O00Ha6BthbFZlzwt9O0L62G/S8HH2Ou3mna/1SacQFfZ5rk7oxoSrHBPQUZcmRa07QHrX+WWaPK2S3Cg4xg1w1d198D3nO89uuDfHab2+PGCXmo2rT/qfIEfDY+GRDl3kpF7H958Gk7wH77n2ox+9gQZfaiDdv+58TPJEsZ+58DH3NpZex/g+FHR1EV+18O/lxo2/43fyFiwSJPj100Pv2eB9j/7glZNOXEp06TvNODrN6FWLEOMS1s0ZT6DLfM4wLmXSmAjx3v3QR+w4cylSym/W+2WnvRZvAhl7anwR9YFL79Lx/986Se0W//ewJx933eHrxY/GqIMva/lSrY/xz41RC1t5r2vw1Ph2//y8GnmBtum/a/F7+siGgqC7/u7WuW/e8rb4cvmnLiV6fp49Om/U8F0dRgeB3zCcv+p33H4wN7gU/39qFqaG9Vc6BLQYZMpnVdLwPQNC2t67piwyvxQl69BBdOI76qgP0vLz5k0+dPwY0riG/9grfH9QwfYj5x0LL//ay3x/UCnxpfeXCPOva/HHxKLG37n6qzTOB547t7B6gkmsrGr05TVQpGj1ZENOXEp5jTKZgwEbFms6fH9QSf7m0jXWHZ/xQQTTnxq96urIB75sAjy709cExxlbUlibT/mJXyGMR6FStlfFnyYFbKk9SslMGXhy6NdAqmz4Tlazw9rjf4lGSlLfvf4qWeHtcTfFjKbNv/HnwYcf8C709QLD5cZlM09T48tlIN+18O3ieWsrsLWfshYuVGNex/TnxILGVrC/JgNWLdNkXsfw586DTJ+lumaGqTqvY/H2K+aommtryg2PMe6jLonaFp2m5c/Ap1XVfBgxxpZFenWSmrYv8bDC8b39Zm5OFqxJYX1bD/5cPjdeLy9k04eRjxudfVrJT9aIguX4SLZxGvf0fRStmH0btzH8HNa4hv/5J3x/QUH7Lp4weguYGyt37eu2N6iR+dpgMZ+5+KMy7402mq3WXa/5R73iODD3VYlWX/26yo/c+PTlN6uzm4t+Fpz44Zd4Zq0f8ssFKUOPLAbuhSbC9aJ2VlePpjrdkFvYrtRetElHnbEGWUrJsVfTLah+0AZVUKxoxFrFe0Ui7zoyFKwcTJiNWbPDumpwhrQtLDhMNIp2BGOSxTVA3gQ8wynYLZ82DREs+O6SmZmD26t81tD7fDwkcRcx/w5Jie43GHXfb2Ive8D8tWIWbd5emxPcPjelt2dSL3fohYvQkxWeHBPcUYNJnWdf3HQRakVLEfwpszHx5SR9WZi3fLPOy9aBc/hpgzz5Nj+oIApOHJoWTvHXMXj8fXIGaqXil7FHNXB3JvGrFmM2LSFE+O6T1WzIZH93ZLI/JwLeKpzyHGhq/Rzos9eCc9GbCVt27AqSOIV99EjFJw20PIGqX16Dp/egEunUd87R1FZ1zA63ubsyfg1nXEF97w5nh+4HEdxrF90NKEkjtsZfC63t5fZQ7uqRyzgriea9Y07duYJsQ5mOKVv9J1/c/9KljJcOkCXP4Y8eb3FK6U8Xb98GnL/vfq1705nm94t5uHPLLXsv8pPBLv8RSp3JtWz/7nxOMpUrlnp2n/26rolDB4vpxHVqXUFE0NwOOYM6KpdU95cjxf8Hj/YaOyAiZPQaxSSTTlwOOJJqMyI5p6wpsD+oJ3I9P2todzH1BLNBUBXD2AqGnavwd+Dfhb4Bet//5r6/WEIpDpChg3HrFuW9hFGQbvRqaNdIWC9r88eLhmWqZTZqW8ROFK2cO+nP0Q3v0LYMFi7w7sNR42vvZetA8vQ9yrkv3PiXeJpW3/W7EWMV0d+18OHnaaZEc7cl8a8eQWRUVTFh52mkzR1F7EhmcVs/858TDmm1fhzHHE5ufVFE1l8LLTdOk8XL6ooGhKfdyOTH8X2Kbr+qeZFzRN2w5UAb/jR8FKAdnRhjxQhVi7DTFBtb1oHXiUWMrGeji2H/HcF9Wy/+VBeLSDiW3/+/I3Fd32MIOHo3cXz5r2v7e+r3il7GE2ffKoaf977e3ij+UnXnaaDteYoimVZx+y8aLTtK9SXdHUADxMLPfstERTCs+4gLedpqqI2P+87DRlBvfWbiv6WKWG25Z9EnDb8Vo9oODeONFB1lZCT08EKmU823K53/6neKUMeLXMQ1alTPvfRkW3PczgZUOUsf+tVXyzHw8bIiNdYdn/1hZ9LH/xuPG9//vWAgAAIABJREFU6154RDX7nwOP7m17xmX+QsQDitn/nHjUaZJGn6mVVtH+NxhF3tqmaOqDiNj/PLq329tM0dTareoP7imI25HpFPDXmqb9GnAZmI85Ir3dr4LFHXttkqr2vxyKH6WVfX3mjhYq2v/yIaDoCkpl+58Tj3Jp2d6KPLhHUfufE2/WG8qG23D8oJr2PydeJZbXLsP5U2ra/5x41WmKgv3PxqOYTxyGhjrKtO8WXyS/8erePlQTHfufV/X23g+jM7inIG5rwB8ArcAxoA04CrQDqqrr1Of8SdP+F5UnZr1Y5pGx/0Up5mI7EAf3QEe7mia8HLxpfJW2/znxaL2hbf+LwoyLV42vbf97tugi+Y9HnaaM/e9JxWdcwLPE0rT/zVTT/ufEo06TTEfJ/ld8zPaMy4LFiHkPeVSu0mLIZFrTtDIAXddbdF3/JjARuA+YqOv6N3VdbwqgjLHEtv+tVtT+l48iGyIjXaGw/S8fxc+TysoKuHcuLH7Mg/L4jAeNr10pP/SImvY/Jx40vrb9b+kqRPk9HhXMTzyIubsLWfsBYpXioqkMHnSaTPvfHsS6p9QVTWXjQS4t6z4z7X+bVbX/OSm+02Tb/6LyEJ4XnabM4F4UBkAUZbhfxzVN0/4Kcxu8E7quG8CtAMoVa2RLE/JQDWLbS4hxiu5F66TIkWnT/ncE8fnX1d2L1kmRI9Py8sfwyTnE69+NRKVsl7GYTtPZE5b975c9KZP/eNAQHd9v2v+2ft+TEvmOF52m/VXQ2RGdxteLTlNGNBWZaXAPYs6IpjY9702RfEbYj7kUEXM6FS37nxedpsoKmDhJXdFUBBgumf4e8BawX9O008CPgf9X13Xnw4gJI0DW7DL3oo3ClHCGYhPL3duJUqUMFL23tkxvV9v+NyhFNkQTJyNWb/SwPD7iQUNkpFMwsxyWrfKkSIFRTOObsf8tVFk0lU1xF1oahrmsZeESxJz53hXLT4rsNNmiqeWr1bX/OSmy0yS7OpC1lYg1mxQWTTkpMuaWJlM0FaXBPQUZcpmHruv/qOv6VzGXdvwp8FXgiqZpP9U07TVN0xR/0kY9zEp5u2n/m62w/S8fBTa+ZqW807L/lXtbJl8pfDRednWYe9Gu2YyYNNnbYvlJEZ0m2dKIPFKL2PCMuva/HIpsiG5dh1NH1d+LNptik6xL5+HTC9GZBofiO01njsOtGxEalc6i0JiPWva/SMVc5L29v8oSTUVkxgWK/z1XW4N7kbrO6uFqEZS1NvpPgT/VNG0BpgnxP1uvRSk7Cp/I2P8cFLHMIxr2vzz0zxmOmEjY//JSRMx7dkJfn9r2PyfFNkRV2yNg/3NQ7Ohd1Xb17X85FBezkU6pb/9zUuS9bWREU0rb/xwUEbL9vEdU7X+FxJyZcVm8FHHf/d6XqYQY0X5GmqaNA9YAa4F7gBN+FCrOGJXWXrSq2/9yKGLEMp2C8nvUtv/lo8AOhF0pz3tQbftfPgrcDtC2/z2yXHH7n4NiGl/b/rdObftfDoUnlrb9b+1Wte1/TopILG3730bV7X8Oiug02fa/LS9EZ8YFKKrTZNv/XorOjAsU12k6dRTqPovOrmIK42pkWtO0TcA3AQ3zAcS/Ar6fbURMGB7ZUAfH9yOe/5Ly9r8cCqxc5I0rlv3vW+rvRZuPQirljP3vG6rb//JQ6DKPk0eg/hZlX3nb8yL5SxFJ1qHqaNn/MhRxS5p70UbB/uegmE5TpERT2RSRTKcz9r8obHuYRTGdpsoKUzS1bqu3ZfKbIjpNRjpliabWeVyo0mPIZFrTtN/AXNIxE/g74HO6rlcHUK5YEt1KmYIfxpPpiNj/8iCEKCivlOkKGD8hGnvR5lDYaLyRTsHU6bBCdfufg2JG79IpuPu+iOxFm01hMQ+w/0VCNJVNYUmW7OtDVr0PS1Yg7o6I/S9DgZ0mWzT1RAREU04KzKVt+9/6KIimnBT4e26og2P7ES9+SX3RVAQYbmR6HfDvgX/Qdb0rgPLEFtv+99gTiLvuDbs4BTDyEUvZ3W3uRbtyPWLqdJ/K5SMFLPOQ7a3IA3sQm56NYKVMQbm0rLfsfy+9FsFKubDtAOW1T+HCKcRXvh29GZdCR+8unIbrlyNi/3NQaKfpo0OmaOqNCNj/cigwyTpYDR1t0Zt9AAruNNmiqYgOdBWA3PM+IBGbIxizggyZTOu67tmvSdO0F4EfAqOAP9N1/fcc7/9L4LtAL3Ab+JlYLSM5cQCa6il78+fCLklhCIEcaQUVKftfHgpY8iBrPoDeO4gtUWyIoKBOk10pR2jbwwwFyjzMvWhHIzZEb8al4NG7dEV07H85mEFLObLcw0inTPvf8gjY/5wU2GmSVSnT/vfwMu/L5DcFdJqib/8rIObID+6pRyBDKpqmjQL+CHgJWAJ8TdO0JY6PHQFW67q+HPgJ8PtBlC0ojHQKps+KkP0vD4U0vvfdHw37X15G1uWPnP0vHyMcjTftfzsiZP9zUEjj292F3PthdOx/ORQQc2sL8lB1dOx/TgroNEXP/ueggFxaXomY/S+HAmaazp2Em1ejPegDjCjo4wegqSF6z3soTFDzk08CF3Rdv6jreg/wt8Cr2R/Qdf1DXdc7rD/3AkpuCSBbW5CHazGaG91/J2P/2/x8dOx/TsrKGFFDlLH/RbZSBspGOEp79gR8di26lTKMfDtA2/4X0ZgLSaYz9r+oPgFfNvLGt9/+F9GYhdXUjeQ6Z+x/UZxxgcJirkqZoqmo2P+cFHJvpy3735qI2v8KMNca6QqYUQ7LIjy4pxhBdbfnAFey/r6Kub3eYHwHqMj3hqZp7wDvAOi6Tnl5sNtc93x2hcY/+d8x5vxflD+20tV3Wt/T6RBlzHr1DUbNiua23PWjRzNqzBimu/z/u+Xv/hudY8dR/rnXKJscxdE7qBOCsWPHMs1lzE21H9AzeSp3Pf9KZE1St8rKmDB+HFNcxtxYs4ve8nso3xbNjqLs7eUWMHHiBCa7jLm+eifMe5CZazdHsqNojC7jNjB54iQmuohZGgb1e96nbMnjzHzcXZ2nGr3d7dQDUyZPZoKbmO/coa56J2NXb2D64qhYHgfSUz+dRmDqlCmMcxGz0dlO3d404zc9y7T50ZxZ6546jSZg2tSpjHURc19TA3VHapn40mtMmT3H/wL6QOfUqbQAM6ZPY7SLmHtvXqP+5BEmvfFdJt8TwdlEi9GjRwee/w1FUMl0vhYnbz9K07S3gNVA3v1pdF1/F3g3c4y6ujpPCugW2dYOQF93D00uzi3v3MHY+U+wfA2NUkDA5fWKvr4+eru7cfP/t+zswEhvR6zZRENXD3RFM2aJoLury13MzY0Y+9KIpz9PfWsrtLYGUELvkRI6OzrpdhPzZ9cxjh1AvPp16hvdz9SohOzrA6CjvZ0uNzFfOo/x8RnEmz9HfX2938XzBdnaDEBbWysdbmI+dQTj5jWMz7/h6regIrKxCYDWlhbaXcRgHNiDbG7kzvpnohtzs3mdW1qaEW5iTqeQXR30rN0W3ZhbWwBobmpyF3PFT6C3l641W13VeSpiWG1NY0MjYtzwe78b//A3UFZG58qNruo8VSkvLw/lPp09O/+uPkEl01eBbL3OXOC680Oapj2LuXvIVl3XuwMq28iw1s7J3juuPi6P1EbT/peD+yUPcl8ldHdFd0o4wwi2A5R7dpj2vyhue5jNCKQtkbT/ORnhrLBMp0z739ptfpUoAEY2LWza/6ZGUDSVxQjXlcp0hWX/W+FfmfxmBEuYpJTmPstzF0TT/mczgpgNw/w9P7wMcZ+Sq0rdMYJ72xZNPf4kYkaURFPqE1QyfQBYZKnIrwFvAG9mf0DTtCcw9eQv6rp+K6ByjZzMVHZfr6uP99v/Ilwpg+tn8exKed6D8MAif8vkO+4expOG9WR01Ox/+XC5g4m804Osydj/ZgZQML8YQePb0YbcXxU9+5+TETyMZ9r/9iGeezV6oqlsRpBLyxtXTdHUl74RMfufkxF0mj45Z4qmvv7zkVy6ZDOSTtMpUzQlXvuWr0XynZF0mg7XQFtLdJ9xUZhAHkDUdb0X+AGwHThtvqSf1DTtNzVNe8X62P8BTAb+TtO0o5qm/TSIso2YUVb/o3f4ZFreuALnPkJseTF6e9E6cbtN3MWzcO3TaD94mMHtw3gZ+1/kZx/Aba9JHqqJpv3PyUhGdfZWRtP+52QkjW+URVMDGEHMVRG1/zkZSacpnYqm/c/JCHLp+Nj/3HeaZFUK7roXHn3c3yKVIIHt96Pr+nvAe47Xfj3r39GouayRaeliZDrK9r9cXI5Y2va/iFfKWAZEFzVUZO1/+XC5tCW69r+BCJdPwkfb/ufEXcYRafufE5edpkjb/5y47DRF2/7nxGXMDXVw7EAs7H/9Yz7DxHz9Mpw7iXjtW9Ef3FOQ5P/RkTLK+uH19g35Mdndjaz5ALFqQzTtf05cjDLLthbT/rfuKcT4CQEUKgCGq5Qz9r9Nz0e+UjYZvtMkr14y7X9bYzDjYjNMNn3+lGn/i+p2eNm4Hb07cdC0/8VpSni4TtPBPab9Lw7X2W2nqXaXaf+LQ8xuO027Y2T/c9tpqtpuiqY2RmPcMmrEpSUMjszI9DAPIMqDu6GzPfpTwtkM92PN2P+iqGTNh4tlHnL3dkAitkR0L1onLqQtsioFo8cg1sdhxgVXS5hkOgUTJiHWbA6oUH7irvE1RVMzoy2ayuA2yUqn4N45sHip/2XyGxch2zMuDz4cXdFUPoaKubfXsv+tjIn9b/gLLbu7zMG9lRsRU6YFU6wSI0mmR0pmzfQwyzxkOmXa/xZF1f7nYJgka4D9b25MKuXhYu7tNdeULl2FmHV3cOXyk2FyadnViaz9ELE6qva/PAzTaZKtzcjD1Yj1EbX/OXGRWJqiqcPWjEsE7X9OXIzeycsX4eLZeDzvAbjqNJ37CG5ei8+gj5tO0/EDlmgqLjFb/x2qDjsQw8E9xUiS6ZFij0wPnkzLT2Ng/3My3OjdmeNw63r0t8MbwDAxH9sPzY3xmgYfJpuWB3ZDV2fMKuVhYs7Y/7bEJGY3ja9t/4vwtocDcJFMZ+x/6yNq/3PiptOUTpn2v9URtf85cdFpMtIpy/63OqBC+Y2Le9se3FsSUJlKjySZHikuRqZlVQrGjkWseyqgQgXEUI1vOgWTpiBWRXgvWifDdISMdAXMLIdl0bTC5WWYTpNMp2DOfHgomla4vAzx0KW9F+2iJYg584Itl28MnWTJ3jvmjMvy1YiZdwVXLD8ZZkxDdnUg96YRqzchJk0Jpkx+M0wuLVsakYdrERueQYyNprE1l6ETS3nrBpw6gtgcTWNrXobpNMlPL8Cl84itL8VncE9BkmR6pNjSlvzJtOzsQO5LI9ZsRkyaHGTJ/GWIJQ+yuRF5dC9iw9MxqpQZOsn67DqcPobY/ELE96LNw2CN76Xz8OmFeM24wNDLec4cg9s34zXjMszonTyyzxJNxSjm4ZKsfVXQHccZFxi03q7eBX0xmnGB4TtNcRBNORmu05QRTa3fFlSJSpIkmR4pmd0LBkum91Za9r+XgytTEAwxYtlv/4tRpQwMtcxDVqWsvWhjVCnD0J2mygoYNz5+My4MvmbaqKyIvv0vh6G3A+y3/z0RXJH8ZojRO1s0dX/U7X8Ohug0SaMvHva/HIaIOWP/W7E2Zva/Ie7tjnZTNPXkFsTEGA3uKUiSTI8QIQSMGp13n2nzIbwKmPcQPBD1vWjdIY0+s7f/6OOIe+eEXRxvGWSfaXmnxxzVWbE24va/PAzSaZIdbcgDVqU8Iep70ToYLObGeji2H7Hx2Wjb/5wMIfOw7X+x2vaQoUfvLp417X+xmwY3Y8k7HnDyqGn/i9XsA0N3mmz7X8wGfYbqNO2rtAb3YhazgsSotgyQUaPyr5n++Ex87H9OBhuZ/ugwNNyOXwUFg+7yIA9VQ3trzKbBM+S/b2VtJfT0xLNSFpC38Y2N/c/BUI1vVZxEU9kMEXPG/rd2S8Bl8pkhOk1GusKy/8VANJXNEJ0mma4w7X+PxM3+l3+maYBo6oFFgZeq1EiS6UIYNTrvmmmZTln2v5hVysBgSx6MdAqmzYDHY1Yp2wzS+N49Gx5eFkJ5fCbPOnF7xuWBRTGw/+Ujt9Mk+/rMHS2WPIG4+75wiuUb+RNL2/63Mgb2PyeDjFjK9lbkwT2IdVtjYP9zMEinSTZkRFPPxUQ0lc0gMV+7DOfjJpqyGKzT9PFpe3AvwX9idlcFxOjROWumZVuLVSnHyP6XTZ4BS1l/C04cRGx8Lh570ToQeUbjTfvfacTWF+JXKQN5O03nT8GNK/GtlPOtE7ftfzGMeZBJM9P+F9O9aAcZsZS1H5j2vzjPMjnrsN07AInYHBPRVDaDdZqqUqb9b0PcZlwYvNOUTsGEiTERTalPHLMB/xk1KmfNtLkX7Z14NkSQd5mHuRct8bH/OcmzzMO2/8WxUoa8iWW87H+55Os02fa/x58MqVR+MkTje+/ceNj/cshNsmJr/8uQbwAkdva/Qci6tWV3lymaiq39L/f3LFtbkAerEeu2xXNwT0GSZLoQRg0cmTb3ot0OCx9FzH0gvHL5iShjQEOUsf8tWx0f+58TZ8y2/W8TYnJM7H9OHB0I2/634el42P/y4Yw5Y/+L01602ZTlSSxjZ/9zIKymLrsDcfaEZf+L46g0+WM+vt+0/22Lacz2bGHWvb2/yrT/xTXmPL9XWZsZ3ItpzAqSJNOF4HwA8ewJy/4X01HpDEZWpXxsn2X/i3HMjvXDcn9VDO1/DoQAadh/yuqdlv0vZg/hZeOMefd2QCA2xXTGJTOSlfV7lumY2f+c2P0HR8wTJyNWbwynTH6TJ2YjnbJEU6vCKZPvDHJvz54HC2MkmsrGXuZhWP+xRFMLlyDmzA+xYKVFkkwXguMBRCNdYdn/YlopQ870v1kp3wVLY2T/c5I1/W8/hDdnPjz0SMgF8xnrMkvDMLc9XPwYYnZc7H/56B/ZMe1/O+HxNYiZ5SGWyUcc60plV0xFU9k41pXKlkbkkbjZ/5w4Yr51HU4dNWdc4iaayuB4GC+2oqlsnBNNZ47DrRvxHvRRkCSZLoRRo6H3DgCyqQGO7kNsfAYxZmzIBfOR7MTy5jXT/rclhva/AWStpb10AS5fjOFetA6yO02nLftf7GQ8DrLv7SN7LftffGO2799Mp2lvOob2PyeODsSenaZoamvMZ1wAO+Y42v+cODtNVdtN+1/sRFPZDIzZSKdg8hTEqjiJptQnSaYLYdQoZF8fQL/9b3OMK2UHcvd20/638dmwi+IvWUmzTGfsf9vCK08QZCWWRjqO9r88ZC3nkemUaf9bEiP736DI/ofw7l8ACxaHXSD/yMorbdHUw8sQ98bJ/jcIkoH2v+lxsv856b/QsqPdnHF5cgvi/2/v7mOkqs44jn/vzirK2/KygC4goGApWCpVwagBq6bVaMQaedAGY1Sqf1RrX1Kj9o82mlpNmlbSmDYGbTUa7RNrU/4oYsSKr0UFbRBpUxQqCwgsL8JKBZad/nHvwHhZhZ3OzN295/dJJswdzu59bs7O3GfOPfc8fftlGlVNlX1pKu7cBu/8neici/I9uNcDKZmuRGMjHOiIP5Rfzmn1v7Qkycp19b+0JMnKdfW/rhTJb/W/LsWj8bmt/teV0pem3Fb/Sysbvctr9b+08iRrxWvQvjvXV1yAz35pWvYi7Nub8ysuZYrJFZc8FprqBXJ+xqiRQiGeM71yBWxvy2klvLQkmc519b+0JMkqVf/L693g5ZJpHger/4VwIkpW8zhU/S/nV1yAUqGag4Wm8lb9L60ssexcuggGDspf9b+0sikPxaWLYPiJOaz+l1Z+zIFU/yv1c2lwb9LpRCNaso0pQEqmK1GIR6Y7ly5Kqv/lcS3alNJn1IuLYMRImDgl23jqoTQav3QRjDuV6KRTso6oDiI4kFwGnzw132vRlkQRfKb636CsI6q9CPikVGjq/PxV/0sr5dLb23Jc/S8tOehS9b8ZgVxxAYpr3gun+l/pmFe+lQzuBXDMPVDO31k1UmjkwOZN8O7y5EM5f9X/DhNFsHkjvP/P+MbDXF8SjkVRBGv/ne/qf2lRRHH1P2DntkCuPsSK7yxLqv8FcsxRRPHNV+Lqf3m/wRQoJZbF1/9Gbqv/pZUGQF5bkt/qf2mlY351SVz9b1rOr7gAB/+2X10CTUNgSgCDez1Q3bJAM7sYmA8UgAXufl/q//sAjwFnANuAOe6+rl7xdUuhQOeONogaArrxMII97Un1v5yuRduVPe1x9b8z81n97zAR8TEPGgJTzso6mrqIoojinvak+t/krMOpk+T9fMrEfFb/Syt9+d/THheaah6RbTx1ceiYo+kziQbktNDUZ5Qd89cvzW+hqXKlca097USXzQljcK8HqsvItJkVgAeBS4BJwDVmNinV7EZgh7uPB34N3F+P2CpSSP5Yv3IG0dBh2cZSL8nJKNfV/9JKx3zOBUR98roWbVpyzHmt/teVUj/neS3atOQwwxiVhvK1xIO5DF72txzSFZeDT0Pp54Nv5oYcF5rq+eo1zWMasMbdP3D3fcBTwKxUm1nAo8nzp4ELzaxHntlKSUYwH8plwvmA4lCSFdKd0VGUrEUb0IdyFMGxOa7+16Uo39X/0kpnklxX/0spHXOeq/8dJjnokKr/lb5ATDkznMG9Hqhe1wNGAuvLtluB9K3UB9u4e4eZfQwMBdrKG5nZTcBNSTuam+tfpWx3y2j2jRrLkJnhjN7tGnEiHXsnM3j6ecGM3n08ZBjR1LMZPCXHVR5TdjQPp2H8RJpOzXmVxzLbBzfTZ/pMBo4Zm3UodbN1SDPHz/gG/VtyvqRnonigg60DBtLv8qvpNzyEKR7Q2a8vW/scx8BvzeX4YWEkWQcaG2hrbKTpyrkcl0FukIWO/SezraHAoCvn0ieQYwZobGzMJP/7PFGxWDxyq/+Tmc0Gvunu85Lta4Fp7n5rWZtVSZvWZPv9pM22L/jVxY0bN9Yw8s/ZaWcnzU1NbNu9u+77zkqx80BcnCagheCHDh5E25atAayzfEixowMgqHl3Q5sG0rZ9RzBfjAGK+/bG9z/kfXWHMkMH9Kdt1+5gBgMAinv3BjRFLTZ0QH+27W7POoy6CrGfm5ubaWtrO3LDKmtpaYHyeWOJen2StgKjy7ZHAeks+GAbM2sEmoDtdYmum6KGhuD+cKOGQlCJNEBUaAwqkYY4iQ4pkQaIjjk2qEQaIDq2T1CJNEDU57igEmkguPMUEMZNhykh9nNPU6+z5pvABDMbB2wArga+nWqzELgOeB24CnjB3Ws/bC4iIiIiUqG6DE24ewdwC7AYWB2/5KvM7G4zuzxp9jAw1MzWAD8E7qhHbCIiIiIilarLnOkaymTONGQ3X0fqR30cBvVzGNTPYVA/hyHUOdMiIiIiIrmjZFpEREREpEK9fppH1gGIiIiISDByN80jyuphZsuz3L8e6mM91M96qJ/1UD+H+Mi4nw/T25NpEREREZHMKJkWEREREamQkunKPZR1AFJz6uMwqJ/DoH4Og/o5DD2qn3v7DYgiIiIiIpnRyLSIiIiISIUasw6gtzGzi4H5QAFY4O73ZRySVJmZjQYeA04AOoGH3H1+tlFJrZhZAXgL2ODul2Udj1SfmQ0CFgCnES+peoO7v55tVFJNZvYDYB5x/64Ernf3T7ONSqrBzB4BLgO2uPtpyWtDgD8CY4F1gLn7jqxi1Mh0NyQn3QeBS4BJwDVmNinbqKQGOoAfufuXgbOB76qfc+02YHXWQUhNzQeedfeJwFdRf+eKmY0EvgecmSRbBeDqbKOSKvoDcHHqtTuAJe4+AViSbGdGyXT3TAPWuPsH7r4PeAqYlXFMUmXuvsndVyTPdxOfeEdmG5XUgpmNAi4lHrWUHDKzgcAM4GEAd9/n7juzjUpqoBE43swagb7AxozjkSpx95eA7amXZwGPJs8fBa6oa1ApSqa7ZySwvmy7FSVZuWZmY4GpwLKMQ5HaeAC4nXg6j+TTycBW4Pdm9raZLTCzflkHJdXj7huAXwIfApuAj939uWyjkhob4e6bIB4AA4ZnGYyS6e7pqvKNlkPJKTPrD/wJ+L6778o6HqkuMyvNwVuedSxSU43A14DfuvtU4BMyviQs1WVmg4lHKscBLUA/M5ubbVQSEiXT3dMKjC7bHoUuJeWSmR1DnEg/4e7PZB2P1MS5wOVmto54ytYFZvZ4tiFJDbQCre5eurr0NHFyLflxEbDW3be6+37gGeCcjGOS2tpsZicCJP9uyTIYJdPd8yYwwczGmdmxxDc4LMw4JqkyM4uI51eudvdfZR2P1Ia73+nuo9x9LPF7+QV312hWzrj7R8B6M/tS8tKFwHsZhiTV9yFwtpn1TT6/L0Q3mebdQuC65Pl1wF8yjEVL43WHu3eY2S3AYuK7hR9x91UZhyXVdy5wLbDSzN5JXrvL3f+aYUwiUrlbgSeSQZAPgOszjkeqyN2XmdnTwAri1ZjepodVyJPKmdmTwPlAs5m1Aj8F7gPczG4k/jI1O7sIVQFRRERERKRimuYhIiIiIlIhJdMiIiIiIhVSMi0iIiIiUiEl0yIiIiIiFVIyLSIiIiJSISXTIiIiIiIV0jrTIiK9iJndBZzs7vNqvJ/zgReAPcBV7v7sEdq/T1wh1lX8RkRComRaRKQHMbP2ss2+wF7gQLJ9s7vfW8dwNrr7qKNp6O6nmNnPgPG1DUlEpGdRMi0i0oO4e//SczNbB8xz9+ezi0hERL6IkmkRkV6kNPrr7nPNbCywFrgBuBvoD9wJLAceBk4CHnf3W8p+/gbgx8AJwBvATe7+n6Pc9/jk954O7AeWuPuc6hyZiEjvpBsQRUR6v+nABGAO8ADwE+AiYDJgZjaT+MkVwF3AlcAw4GXgyW7s5x7gOWAwMAr4TZV+sDHdAAABUElEQVTiFxHptZRMi4j0fve4+6fu/hzwCfCku29x9w3ECfPUpN3NwC/cfbW7dwD3Aqeb2Zij3M9+YAzQkuzvlSofh4hIr6NkWkSk99tc9vy/XWyX5mGPAeab2U4z2wlsByJg5FHu5/ak/RtmtiqZMiIiEjTNmRYRCcd64Ofu/kQlP+zuHwHfATCz84Dnzewld19TxRhFRHoVjUyLiITjd8CdZjYZwMyazGz20f6wmc02s9JSeTuAIoeW7RMRCZKSaRGRQLj7n4H7gafMbBfwLnBJN37FWcCyZC3shcBt7r62+pGKiPQeUbFYzDoGERHpYcxsBrCYuGjMHHdffIT2/yKee+3urrnUIhIMJdMiIiIiIhXSNA8RERERkQopmRYRERERqZCSaRERERGRCimZFhERERGpkJJpEREREZEKKZkWEREREamQkmkRERERkQr9D2rKMbGAeNOOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from pylab import *\n",
    "%matplotlib inline\n",
    "rcParams['figure.figsize']=(12,3) # Change the default figure size\n",
    "\n",
    "I=1                          #Set the parameter I.\n",
    "C=1                          #Set the parameter C.\n",
    "Vth = 1;                     #Define the voltage threshold.\n",
    "Vreset = 0;                  #Define the reset voltage.\n",
    "dt=0.01                      #Set the timestep.\n",
    "V = zeros([1000,1])          #Initialize V.\n",
    "V[0]=0.2;                    #Set the initial condition.\n",
    "\n",
    "for k in range(1,999):       #March forward in time,\n",
    "    V[k+1] = V[k] + dt*(I/C) #Update the voltage,\n",
    "    if V[k+1] > Vth:         #... and check if the voltage exceeds the threshold.\n",
    "        V[k+1] = Vreset\n",
    "        \n",
    "t = arange(0,len(V))*dt      #Define the time axis.\n",
    "\n",
    "figure()                     #Plot the results.\n",
    "plot(t,V)\n",
    "xlabel('Time [s]')\n",
    "ylabel('Voltage [mV]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** Adjust the parameter `I`.  What happens to `V` if `I=10`?  If `I=100`?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** Adjust the parameter `C`.  What happens to `V` if `C=0.1`?  If `C=10`?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"question\">\n",
    "\n",
    "**Q:** What is \"spiking\" in this I&F model?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"donate\"></a>\n",
    "## Donate\n",
    "If you enjoy Case-Studies-Python, and would like to share your enjoyment with us, sponsor our coffee consuption <a href=\"https://www.paypal.com/donate/?hosted_button_id=DL8P5ZGS9962U\">here</a>."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md:myst"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}