{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Solving Ordinary Differential Equations\n",
    "====\n",
    "\n",
    "## Unit 13, Lecture 1\n",
    "\n",
    "*Numerical Methods and Statistics*\n",
    "\n",
    "----\n",
    "\n",
    "#### Prof. Andrew White, April 24 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import sqrt, pi, erf\n",
    "import seaborn\n",
    "seaborn.set_context(\"talk\")\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Goals:\n",
    "---\n",
    "\n",
    "0. Know the definition of ordinary differential equation\n",
    "1. Be able to determine the order of an ordinary differential equation\n",
    "2. Convert ordinary differential equations into standard form\n",
    "3. Implement ordinary differential equations into python\n",
    "4. Interpret the output from the ordinary differential equation solver in scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Classification of Ordinary Differential Equations\n",
    "===="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The word ordinary means that there are no partial derivative terms nor are there boundary conditions (only initial conditions)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "1$^\\textrm{st}$ Order vs. 2$^\\textrm{nd}$ Order\n",
    "---\n",
    "\n",
    "A 1$^\\textrm{st}$ order only has first derivatives and a second order has second derivatives and often first derivatives."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "*First order Example:*\n",
    "\n",
    "  * $\\frac{dC}{dt} = -rC$ The rate of chage of concentration of a chemical species as a reaction proceeds\n",
    "  \n",
    "  * $\\frac{dV}{dt} = kt$ Filling a tank at a constant rate ($k$)\n",
    "  \n",
    "  * $\\frac{dV}{dt} = kt - A_{out}\\sqrt{2g\\frac{V}{A}}$ Filling a tank with a hole in it\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Second order Examples:*\n",
    "  \n",
    "  * $\\frac{d^2x}{dt^2} + k x = 0$ An ideal spring-mass system\n",
    "  \n",
    "  * $\\frac{d^2x}{dt^2} + b \\frac{dx}{dt} + k x = 0$ A spring-mass system with dampening\n",
    "  \n",
    "  * $\\frac{d^2T}{dx^2} = x$ Fourier's law of heat conduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Solving ODEs in Python\n",
    "===\n",
    "\n",
    "Python can hanlde any order of system of ordinary differential equations. It cannot handle coupled (PDEs) or boundary problems. For those and all other things not in `scipy`, see this list of [numerical methods in python](http://www.scipy.org/topical-software.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To solve an ODE in python and all other programming languages/packages:\n",
    "\n",
    "1. Convert to standard form\n",
    "2. Implement the standardized equation as a Python function\n",
    "3. Create a grid of points where you want to evaluate the ODE\n",
    "4. Call `odeint` with the function, initial value, and grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Converting to Standard Form\n",
    "===="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The standard form is:\n",
    "\n",
    "$$ \\frac{dy}{dt} = f(y,t)$$\n",
    "\n",
    "All first order, second order and $n$th order ODEs can be conevrted into this form. Let's see some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Converting 1st Order to Standard Form\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The chemical reaction example:\n",
    "  $$\\frac{dC(t)}{dt} = -rC(t)$$\n",
    "is already in standard form. So are all the other first-order examples given above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Converting 2nd Order to Standard Form\n",
    "---\n",
    "\n",
    "The trick for the 2nd order form is to create a second dimension, the dimensions of the derivative, which turns the problem into 2D 1st order ODE. Let's see the spring mass system first:\n",
    "\n",
    "$$\\frac{d^2x(t)}{dt^2} + k\\, x(t)= 0$$\n",
    "\n",
    "Now we'll use a $x_1(t) = x(t)$ and $x_2(t) = \\frac{dx(t)}{dt}$. I'll drop the function notation for simplicity now\n",
    "\n",
    "$$\\frac{dx_2}{dt} + kx_1 = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\frac{dx_2}{dt} = -k x_1$$\n",
    "\n",
    "$$\\frac{dx_1}{dt} = x_2 $$\n",
    "\n",
    "Now we have two dimensions and their ODEs are both in standard form. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's see the Fourier equation as another example:\n",
    "\n",
    "$$\\frac{d^2T}{dx^2} = x$$\n",
    "\n",
    "$$T_1 = T$$\n",
    "\n",
    "$$T_2 = \\frac{dT}{dx}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Original diff eq: $$\\frac{dT_2}{dx} = x$$\n",
    "\n",
    "Result of our definition: $$\\frac{dT_1}{dx} = T_2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Implementing the functions in Python\n",
    "==="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We need to implement the right-hand side of this equation:\n",
    "\n",
    "$$\\frac{dC}{dt} = -rC$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def rxn_d(C, t, r=1):\n",
    "    #I added a default argument so that I can easily solve\n",
    "    return -r * C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's see the spring-mass system:\n",
    "\n",
    "\n",
    "$$\\frac{dx_1}{dt} = x_2 $$\n",
    "\n",
    "\n",
    "$$\\frac{dx_2}{dt} = -k x_1$$\n",
    "\n",
    "Our function should return the RHS of the governing ODEs for each equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def ideal_spring(x,t,k=1):\n",
    "    #This will get x as a vector!\n",
    "    return np.array([x[1], x[0] * -k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We'll do Fouier's law now:\n",
    "\n",
    "$$\\frac{dT_1}{dx} = T_2$$\n",
    "\n",
    "$$\\frac{dT_2}{dx} = x$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def fourier(T, x, k=1):\n",
    "    return np.array([T[1], k * x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Using ODEInt In Python\n",
    "===\n",
    "\n",
    "Let's solve the reaction system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEOCAYAAACaQSCZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXGWd7/HPqaW36n1JOntnfbJvkEDCHhQEBEZREUcRRJyZO6Ozeee6vdxGZeZeR0cYdXTGGUEQETCKBibsIpAQEgjZn+7sC0l6S+9bVVfdP04ldppOunqpOt1d3/fr1a/qeuqkz++XQL55zvIcJxaLISIi0h+f1wWIiMjooMAQEZGEKDBERCQhCgwREUmIAkNERBIS8LqAZKmpaR705V+O41BSEqKurpV0uYos3XpOt35BPavnxJSV5Tnn+kwzjD74fO5vui+NfnfSred06xfUc7pIZs9p9NsoIiJDocAQEZGEKDBERCQhCgwREUmIAkNERBLi2WW1xphPAd8FvmKt/fY5trkd+DxunXXAp621r6euShEROc2TGYYx5vvAu4A959lmMXAv8EFr7WzgB8CvjDEZqalSRER68mqG8bC19mVjzIvn2eajwDpr7Q4Aa+0DxphvAVcCT/e3g8Feh3zg7Sbue3wbN10+izXLJgz8B4xSPp9z1utYl279gnpOF8ns2ZPAsNa+nMBmc4EtvcaqgAUkEBglJSEcZ+C/Ya9X1lHX1Mmjz1Vy0+UzCPjT6zRPYWHI6xJSKt36BfWcLpLR80heGiQEtPcaa4+P96uurnVQM4yKcTkAtLSHeW3bMeZNKxr4DxmFfD6HwsIQDQ2tRKNjfwmFdOsX1LN6Tkxxce45PxvJgdECZPcaC8XH+xWLxejuHvhOi/OymDoul8PVLWzeU82cyYUD/yGjWDQao7s7Pf7HgvTrF9RzukhGzyP5eMtOwJx+Y4xxcA9TbUv2jpebMgDesDVps2CZiEh/RnJgPAhcb4xZFH//SdzZxUvJ3vEF8cCob+7k0MnmZO9ORGRUSPkhKWOMH3f2ADAVmG+M+SSwNj7Waq39hrV2lzHmL4BfxC+lPQ7cbK2NJLvGKeNyGVecQ3V9G29U1lBRnp/sXYqIjHgpDwxrbTfuoaVEtn0YeDi5Fb2T4zhcvLCcJ17az5uVtbz/8pmpLkFEZMQZyYekPHXxQvcejGO1rZysb/O4GhER7ykwzmF+RTF5OUEA3qiq8bgaERHvKTDOwe/3sXRWKQBvVCowREQUGOdx+mqpfceaONXc6XE1IiLeUmCcx4IZxWRl+AHYbKs9rkZExFsKjPPICPhZOts9LLVljwJDRNKbAqMfF5pxAFQdbaShRYelRCR9KTD6sXB6MZkZfmLAFquT3yKSvhQY/cgI+lkyswSAzTosJSJpTIGRgBVz3cNSlUcaaGzt8rgaERFvKDASsHBGCRlBHzHgDV0tJSJpSoGRgMygnyUz3aulXtdhKRFJUwqMBF0YPyxljzTQqKulRCQNKTAStHhmiXu1VEyzDBFJTwqMBGUG/SyL38T32u6THlcjIpJ6CowBuGjeeMBdW6q2od3jakREUkuBMQALphcTynKfObVJh6VEJM0oMAYg4PedOfn92i4dlhKR9KLAGKCV8cNSR6pbeLu21eNqRERSR4ExQGZKIQW5GYBmGSKSXhQYA+TzOayc684yXtt9klgs5nFFIiKpocAYhIsXuIFRfaqd/cebPK5GRCQ1FBiDUFGex/jiHAA27DjhcTUiIqmhwBgEx3FYvbAccM9jRLqjHlckIpJ8CoxBWjXfPSzV2hFh+746j6sREUk+BcYglRZmY6YUAvCqDkuJSBpQYAzBqvhhqbf21dLSHva4GhGR5FJgDMGFZhzBgI9Id0wr2IrImKfAGIKcrMCZFWx1tZSIjHUBL3ZqjFkB3AeUAmHgHmvtA31s92fAZ3CDrQn4orX22VTW2p/VC8vZtLuavccaOVHfRnn8clsRkbEm5TMMY0wmsBb4obV2FvAB4F5jzKJe260C/gm4zlo7D/gC8GtjTEmqaz6fBdOLKQi5S4W8sv24x9WIiCSPF4ekrgaw1t4ff90JrANu67XdEmCPtfZwfLvngExgeupK7Z/f52P1Ivfk98vbj9Md1T0ZIjI2eXFIai5Q1WusEljea+x54JvGmEXW2u3GmJuBE8CORHbiOA6+Qcahz+ec9dqfK5ZO5KmNh2ls6WLXwVMsjZ/XGE0G2vNol279gnpOF8ns2YvACAG9H1fXHh8/w1pbaYz5EvCmMeYU7uziVmttRyI7KSkJ4ThD+w0rLAz1vxFQXJzL/OnF7DpQz8bd1ay5qGJI+/VSoj2PFenWL6jndJGMnr0IjBYgu9dYKD5+hjHmeuBzwBxr7f74OY4XjDHXWGvf6G8ndXWtQ5phFBaGaGhoJRpNbDXaVQvGs+tAPZt2nuDgkXry4+c1RovB9DyapVu/oJ7Vc2KKi3PP+ZkXgbET+GyvsXnAtl5j1wPPWmv3A8QPS70FrAH6DYxYLEZ399AKjUZjdHcn9ht+wZwyHsyopLOrm1e2HeealVOHtnOPDKTnsSDd+gX1nC6S0bMXJ71fACLGmDsBjDFLgGuAh3pttwO43BhTGt9uKrAU2JrCWhOWlRFgZfzxrX/YdlzPyRCRMSflgWGtDQM3A3cbY6qAB4G7rLXWGHNP/LwFwI+BR4FXjTF7gPXA10fafRg9XbZ4IgDHalvZ/7aekyEiY4snN+5Za7cCq/sY/3yP76O49158IYWlDcnMSflMLA3xdm0rL249xsxJBV6XJCIybLQ0yDByHIcrlrqzjE27q2nt0IKEIjJ2KDCG2eqF5QQDPsKRqJY9F5ExRYExzEJZQVbOc09+/37r2zr5LSJjhgIjCa5cOgmAt2tbqTra6HE1IiLDQ4GRBDMm5jNlnHvzy4tbj3lcjYjI8FBgJIHjOFwZP/m9eU81zW1dHlckIjJ0CowkuXhBOZlBP5HuGH/YpmXPRWT0U2AkSXZmgNXxZ36/8MZRLXsuIqOeAiOJ1lwwGYC6pk7e2lvncTUiIkOjwEiiSaUh5k0rAuC5LUc9rkZEZGgUGEn2rvgsY/ehUxyraelnaxGRkUuBkWRLZpVSkp8FwPNv6BJbERm9FBhJ5vM5rFnu3sj36o4TtGl9KREZpRQYKXDZkokEAz46w9289JYusRWR0UmBkQK52cEzl9g+u+WILrEVkVFJgZEi16yYAkB9Uyeb99R4XI2IyMApMFJkQkmIxTNLAFi/6bBWsRWRUUeBkULXxmcZB080axVbERl1FBgpNHda0ZlVbNdvOuxxNSIiA6PASCHHcbh2pTvL2FpVy8lTbR5XJCKSOAVGiq2cN57C3AxiwPpNR7wuR0QkYQqMFAv4fVyzYioAL287TmNLp8cViYgkRoHhgSuWTiQnM0CkO8rTmzXLEJHRQYHhgezMAGsucJcLefHNY7R1RDyuSESkfwoMj7zrwilkBHy0d3bzwpta+lxERj4FhkfyczK4bLH73O9nNh+lK9ztcUUiIuenwPDQtRdNwec4NLV26bnfIjLiKTA8VFqQzcULxgPw5MZDhCNalFBERi4Fhsfeu7oCx4FTzZ28sl2zDBEZuRQYHisvzuGi+e4sY92Gg0S6NcsQkZEp4MVOjTErgPuAUiAM3GOtfaCP7ZYC/w5MBNqBf7DW/iaVtabCe1dV8NrOk9Q1dfLqjhNcvmSi1yWJiLxDymcYxphMYC3wQ2vtLOADwL3GmEW9tgsBTwLfsdZOBe4G/toY40nIJdPE0hAr5o0D4HevapYhIiOTF4ekrgaw1t4ff90JrANu67XdTUCNtfaX8e1estausdaOybvcblxdAUBtYwcbdpzwthgRkT548a/1uUBVr7FKYHmvsWXAAWPMT4DLgJPAF621LyWyE8dx8A0yDn0+56zXVJhanseKeeN4fXc1v331IJcsnkAwkLo896JnL6Vbv6Ce00Uye/YiMEK45yN6ao+P91QErAGuBT4J3A48YYyZZa2t7W8nJSUhHGdov2GFhb1LSq47b1zIlj3PU9vYwebKWm64dEZK9w+p79lr6dYvqOd0kYyevQiMFiC711goPt5TA/C6tXZD/P39xph7gNXAE/3tpK6udUgzjMLCEA0NrUSjqXuUaijosGphOa9sP8HDT1uWzy4hM+hPyb696tkr6dYvqGf1nJji4txzfuZFYOwEPttrbB6wrdfYXtwZRk8xIKFzGLFYjO4hrrYRjcbo7k7tf2Q3XjKdjTtP0tjaxTOvH+G6i6aldP9e9OyldOsX1HO6SEbPXpz0fgGIGGPuBDDGLAGuAR7qtd0jwGxjzHXx7W7GnZlsYAwbV5jNZfHLap/ccIj2zjF5jl9ERqGUB4a1NgzcDNxtjKkCHgTustZaY8w9xpgvxbdrwD1/8W1jzD7gS8DN1tpTqa451W5cXUEw4KO1I6Jnf4vIiOHJPQ3W2q245yJ6j3++1/sNwIJU1TVSFOVlsmb5JNZvOsL6TUe4atkkCnIzvS5LRNKclgYZoW5YVUFOZoDOcDe/eeWg1+WIiCgwRqrc7CA3rHZPeL+09W2O17V6XJGIpLsBH5IyxswGKnBPQNcA26y1+tssCd51wWSe33KUuqZOHntxH5++ZbHXJYlIGksoMIwxFcBfAh8ByoGed8RFjDEvAz8CHrXWaiGkYRIM+Hnf5TP4z9/t5s2qWiqPNDBnSqHXZYlImur3kJQx5jvAdmAG8DlgPlAAZOCGx3uAl4BvAG8ZYy5MWrVp6OIF5UwZ595I88jzVURj6XUtuYiMHInMMDKA2dbavlbEqwaej3991RjzPmA2sHn4SkxvPsfh1jWz+PYvtnLgeDMbdpzgkkUTvC5LRNJQv4Fhrf2rRH+YtXbt0MqRvsyvKGbZ7FLerKrlsd/v4wJTRlbGmFvlXURGuAFdJWWMqTPGFCerGDm3D62ZRcDv0NjSxboNh7wuR0TS0EAvqy2ij1mJMabIGPPI8JQkfRlflMO7L5wCwPpNR6hp6L3gr4hIciUUGMaYV40x/4q7+N8CY0yw1ybZwC3DXZyc7b2rK8gPZRDpjvLL5/d6XY6IpJlED4Q/A6zEvZz2OaDLGLMLeBPYCswCjialQjkjOzPALZfP4L+f2sOWyhp27K9j4YwSr8sSkTSR0AzDWvsVa+11QCewCPfxqb8AcoA/x11ttveS5ZIElyyewMxJ+QA8+Ewl4cgQ13AXEUnQQC+1ybXWduM+0+LpJNQj/fA5Dh+7xvC1n75O9al2nnrtMDddMt3rskQkDQzopHc8LMRjU8fnsWb5ZADWbTikE+AikhKJ3Om9MNEfZowJxteakiR732UzyA9lEI5EeeiZSmK6A1xEkiyRGcZ6Y8wvjDFrjDFOXxsYY8qNMZ8FqoCrhrVC6VNOVoBb18wCYNu+Ol7fU+1xRSIy1iVyDsPgriH1CJBhjNkCHAM6gGLcBxzNAv4AfNxa+/sk1Sq9XDx/PK9uP87Og6f4+bNVLJheTCir9xXPIiLDo98ZhrW2xVr7JWAy8HHcy2hzgSlAG+4qtQuttVcpLFLLcRw+9p65ZAR8NLV28egLujdDRJIn0eXNrwK+ivtM7V8ntSIZkHGF2dx82XQefWEfL711nFULyjFTi7wuS0TGoESvkvpbYIu1tqH3B8aYbGPM7caY0uEtTRJ1zYopTI0vgf7Tp/bQGdbFbCIy/BINjGXAw319YK1tB24A/s9wFSUD4/f5uOP6ufgch5On2ln70n6vSxKRMSjRwCgF+noexmkP4t7tLR6pKM/n+lVTAXjm9SNUHnnHZFBEZEgSDYwTwPluJ94BTB16OTIUN66ezqSyEDHgv57crUNTIjKsEg2M54G/OM/nwQH8LEmSYMDHXTfMw+c4VJ9q5/Hf7/O6JBEZQxL9S/7/Ae83xnz5HJ9fCeiazhGgojyfG1ZNA+DZzUfZdbDe44pEZKxIdLXaPbj3YHzRGPOGMeZOY8wCY8xUY8ztwLeA+5NZqCTuxksqmDrevWrqJ+t209Ie9rgiERkLEj6MZK39BXAp7s16PwG2AQeAn+Le5f39JNQngxDw+/jUjQsIBnycau7kgfVWa02JyJANdLXa1621lwJzgNuATwAXWmvfp5VsR5aJpSE+dJW71tTmPdVs2Hm+i9xERPo30OdhAGCt3YvOWYx4a5ZPYtu+Orbvr+PBpyuZOamA8UU5XpclIqOUrmwawxzH4RPXzyUvJ0hHVzf//pudhCNRr8sSkVHKk8Awxqwwxmw0xuw1xuyOnzg/3/YXG2O6jTF3pKjEMaMgN5O73zsfgEMnmnn0RU0MRWRwUh4YxphMYC3wQ2vtLOADwL3GmEXn2D4L+E/cJdVlEBbOKOG6i937Kp/dfJQ3q2o8rkhERiMvZhhXA1hr74+/7gTW4Z5E78s3gN8BWiBpCN532QxmTsoH4L/W7aZWj3UVkQEa1EnvIZqL+2S+niqB5b03NMasAq4FLgTWD2QnjuPgG2Qc+nzOWa9jgd/v53+9byFf/s9NtHZE+P6vd/Clj19ARsAPjM2ezyfd+gX1nC6S2bMXgRECev/ztj0+foYxJhv3UNSd1tpOY8yAdlJSEsJxhvYbVlgY6n+jUaS4OJfPfvRCvv6TjRw60cwvX9jPZ25ddtY2Y63n/qRbv6Ce00UyevYiMFqA7F5jofh4T98AfmOt3TSYndTVtQ5phlFYGKKhoZVodGzd8DZjfIibL53Or/9wgGc2HWZyaQ5XLps0pnvuS7r1C+pZPSemuDj3nJ95ERg7gc/2GpuHe+d4T7cAPmPMR+Lvy4GFxpgl1tq/7W8nsViM7iHeShiNxujuHnv/kb13dQX7jjWxfX8dP1tfyaTSXGZPKQDGbs/nkm79gnpOF8no2YuT3i8AEWPMnQDGmCW4z9J4qOdG1toKa+3U+GsFsBH4bCJhIefncxzuvnE+pQVZRLqj3PerbZxq7vS6LBEZ4VIeGNbaMHAzcLcxpgr34Ut3WWutMeYeY8yXUl1TOsrNDvKZWxaTGfTT2NLFvY9to0vPzxCR83DG6qJ0NTXNg27M73coLs6lvr5lzE9jt9gavr92OwBXXTCZj187h2ga3AyeTn/Gp6ln9ZyIsrK8c14tpKVB0twFpow/ucx9mOILW47yu1cPeVyRiIxUCgzhxtUVrJw3DoDHXtzHa7tOelyRiIxECgzBiZ8EnzutCHAfulR1tMHjqkRkpFFgCAAZQT9f+sRFlBVmu1dOPb6dE/VtXpclIiOIAkPOKMjN5O9uXUIoK0BLe5jvPLKVhhZdbisiLgWGnGViaYhP37KYYMBHbWMH33nkLdo6Il6XJSIjgAJD3mHOlEL+/KYFOA4crWnhvse3EY7oHg2RdKfAkD4tm1PGx98zFwB7pIEfrN1BpDsNbtAQkXNSYMg5Xb5kIrdcMQOAt/bV8R+/3ZU2C7iJyDspMOS8blhVwfUXTwPg9T3V/PSpPUTH6OoAInJ+Cgzp1y1XzODq5ZMBeHn7cR58ulKhIZKGFBjSL8dxuO3ds7l08QQAXnzzGA+utwoNkTSjwJCE+ByHO66by6WL4qGx9W2FhkiaUWBIwnyOwx3Xnx0aP31qj06Ei6QJBYYMyOnQuCx+eOrlbcf58W936pJbkTSgwJAB8zkOH79u7pkT4Zt2V/ODtTt0c5/IGKfAkEHxOQ4fefdsbljlXnK7dW9tfBmRsMeViUiyKDBk0BzH4ZYrZp65uc8eaeCfHnpTzwcXGaMUGDJkN6yq4I7r5p5Ze+pbP9vC8bpWr8sSkWGmwJBhcfmSifzV+xcRDPioa+rgWz/bwp5Dp7wuS0SGkQJDhs2y2WX87w8vIzc7SGtHhH95ZCuvbD/udVkiMkwUGDKsZk0u4Iu3X8D4omy6ozF+sm43j/9+n27wExkDFBgy7MYX5fDF2y9kzpRCANZtOMS/Pb6d9k49iElkNFNgSFLkZgf5+1uXnrnBb+veWr75sy2cPKXnhIuMVgoMSZpgwMcd183lT989B5/j8HZtK1//6Wa2VtV6XZqIDIICQ5LKcRyuvmAyf3/rEnKzg7R3Rrj38W3ueQ2tQSUyqigwJCXmVRTzlTtWMH1CHuCe1/iXR7bS0KKb/ERGCwWGpExJQRaf+9MLuGrZJAB2HzrFV/5rEzv213lcmYgkQoEhKRUM+PjYtYZP3TSfzAw/zW1hvvPLt/jlC3sJR7TirchIpsAQT1w8v5yv3rmCaePdQ1T/89phvvnAZo7VtHhcmYicS8CLnRpjVgD3AaVAGLjHWvtAH9t9Bvgz3DrbgH+w1j6TylolecYX5fCFj13Ar17ax9ObjnC4uoWv/XQzH7xyJldfOBmf43hdooj0kPIZhjEmE1gL/NBaOwv4AHCvMWZRr+1uBD4HXGutNcA9wGPGmKxU1yzJEwz4uHXNbD572zKK8jKJdEd5+Lkq/u9Db+ieDZERxosZxtUA1tr74687jTHrgNuA7T222wd8yFp7NP5+PZAPTANsfztxHAffIOPQ53POek0HXve8cEYx3/zURfxsfSUbdpyg8mgjX/nJJj541UzedeGUYa/L6369oJ7TQzJ79iIw5gJVvcYqgeU9B6y1u3ptczNwDNifyE5KSkI4QzykUVgYGtKvH4287LkY+MKdF7Fxx3G+/9hbNDR38tAzVbxua/irDy5l+sSCYd+n/ozTg3oeHl4ERgho7zXWHh/vkzHmSuC7wIettQk90q2urnVIM4zCwhANDa1pc3PZSOp5zsQ8vnn3RTz0TCWvbj9B5eEG/va7v+c9F0/l5kunkxn0D3kfI6nfVFHP6jkRxcW55/zMi8BoAbJ7jYXi4+9gjLkd+DZwq7X22UR3EovF6B7iI6aj0Rjd3enxH9lpI6Xn7IwAn7xhPhfPH88D/2Opbexg3auH2LjjBB++eg7L55QOeQYJI6ffVFLP6SEZPXtxWe1OYE6vsXnAtt4bGmPuAr4GXDmQsJCxY+H0Ev7xkxdx3cVT8fsc6po6+f7a7Xz30bf0VD+RFPMiMF4AIsaYOwGMMUuAa4CHem5kjJkP/DNwdR/nMySNZAb9fPDKWXztEyuZO9VdMn3H/nq+/JNN/PzZSlo7EjpKKSJD5MQ8eLCNMWYp8AOgDOgAvmqtfdwYcw/Qaq39hjHmR8BHcE909/R31ton+9tHTU3zoBvz+x2Ki3Opr29Jm2nsaOk5FouxaXc1j764l/omdx2q3OwgN66u4MplkwgGEvs30GjpdzipZ/WciLKyvHMe6/UkMFJBgTEwo63nznA36zcd5smNh+gKu0uKlBZk8f4rZrBy3vh+b/obbf0OB/WsnhNxvsDQ0iAyKmUG/dx0yXTu+dQqLl08AceB2sYOfvzELr7236+ztaqWsfqPIRGvKDBkVCvKy+QT18/ja59YyeKZJQAcqW7h3se38c2fbWH7/joFh8gw8WQtKZHhNrksl7/54BKqjjaw9qX97DncwP63m/juL99i+oQ8blw9nSWzSoblUlyRdKXAkDFl9uRC/vdty9h16BS/efkAe482cuB4M/c+vo3JZblcf/FUVswbh98/9Jv/RNKNAkPGHMdxWFBRzPxpRew53MATLx/AHmngaE0LP/7tLn710n7ec9FUbrpytteliowqCgwZsxzHYd60IuZNK6LqaANPbTzM1r211DZ28ODTlax9aT9XLJvEVUsnUVKgRZBF+qPAkLQwe3Ihsz9QyNGaFtZvOsxru07S2hHhyQ2HeGrjIZbNLuPqCyYzd2qhznOInIMCQ9LK5LJc7rphPh9aM4tXd1Xz5CsHaG4L80ZlDW9U1jChJIcrlkxk9aIJ5GYHvS5XZETRjXt90M0+Y7/n0/2eqG5k446TPLflKAdPNJ/5POD3caEp49LFE5g7rWhMPP0v3f6MQT0P9417mmFIWssI+Llk0QRWLyznwPFmXtx6jE27T9IVjrJx10k27jpJSX4WlywqZ9XCcsYX5XhdsohnFBgiuCfIZ0zMZ8bEfD68Zjav7TrBH7Yd5+CJZuqaOnjilYM88cpBZk7KZ9WCci6cO478nAyvyxZJKQWGSC85WQGuWj6Zq5ZP5kh1Cy9vO85ru07Q1BZm37Em9h1r4ufPVDGvooiVc8exbE6ZzndIWtA5jD7ouOfY73mg/XZHo+w8cIpXdxxn697aMwseAvh9DnOnFXGBKWPZrFIKcjOTWfqgpdufMahnncMQ8YDf52PxzBIWzyyhs6ubt/bV8tquk2zfX0+kO8rOA/XsPFDPz7DMnFTAstmlLJlVyoSSHF2mK2OGAkNkgDIz/KycN56V88bT3hlh2746Nu+pZvuBOrrCUfYea2TvsUYefXEfpQVZLJlVyqIZJZiphcPyPHIRrygwRIYgOzPARfPHc9H88XSFu9l18BRvVNawbV8tTW1hahs7eG7LUZ7bcpRgwIeZUsj8imLmVxQxZVyuZh8yqigwRIZJRtDP0tmlLJ1dSjQW4+DxZrburWX7/joOnWgmHImy40A9Ow7UA5CfE2TutCLmTi1i7rQixhdlK0BkRFNgiCSBr8dluu+/fAaNrV3sPFDHzgP17Dp4isbWLprawmzaXc2m3dUAFORmYKYUusuYTC5gclkuPp8CREYOBYZIChSEMli9cAKrF04gFotxrLaVXQdPsefQKeyRBto7IzS2dJ0VINmZfmZMLGDWpAJmTsxn+sR8Qlm6fFe8o8AQSTHHcZhclsvkslyuWTGFaDTG4epmKg83YI80UHW0kZb2MO2d3WeuvjptfHEO0yfkMX1CPhXleUwdl0dmhk6kS2ooMEQ85vM5VJTnU1GezzUrpxKLxThR30bVUfdqq/1vN/F2bSsAJ+vbOFnfxsadJwFwHJhQEmLa+FymjMtj6vhcpozLJU93oUsSKDBERhjHcZhQEmJCSYjLl0wEoK0jzIHjzRw43sSB403sP95EY0sXsRi8XdvK27WtbIiHCLjnQ6aU5TJ5XC6TSkNMKgsxeVyuVy3JGKHAEBkFcrKCLJhezILpxWfGGlo6OXSimUMnmzl8soXDJ5upbewAoLGli8aWP16RBeAAZcU5TCjKprwkh/LiHCaUhCgvziEvJ6grtKRfCgyRUapRxj82AAAJ3UlEQVQwN5PCWZksmVV6ZqytI8zRmlaO1rRwpLqFYzWtHKttpb0zQgyorm+jur6Nt/bVnfWzsjP9jC/KYXxxDuMKsxlX5H6VFWZTEMpQmAigwBAZU3KygsyZUsicKYVnxmKxGKeaOzlxqo3GtghVh+t5u7aNE3WtNLWFAWjv7ObgieaznglyWkbAR0lBFmWF2ZQUZFFakEVJfpb7fX4WeaGMMfG8EOmfAkNkjHMch+L8LMqKsuOL0o0/syhda0eYE/VtVNe3c6K+jZOn2qhpaKf6VDutHREAuiJRjte1cbyurc+f7/c5FOVlUpyfRXFeJkV5mRTmZVKUG/8+N5OC3AwCfl/KepbkUGCIpLFQVpCZEwuYObHgHZ+1doSpbeigpqGdmoZ2ahs74l/t1DV1nFmxtzsaO/PZ+eRmBynIzaAgdPork/xQBvmhIPk5GeSHMsjLySAvJ6hwGaEUGCLSp1BWkFB5kGnlee/4LBaL0dIepr6pk/qmDuqaOqhv7qShuZP65k5ONXfQ0NJFOPLHZeBb2sO0tIc5VtPa776zMwPk5QTJyw6Se/orx30NZZ1+DRCKv8/JCpCV4de5liRTYIjIgDmOE58NZPQZKOCGSltnhFPNnTS2dNHQ0kljaxeNLV00tXXR2NJJU1uYptYuWtrDZ/3a9s4I7Z0Rqk+1J1yTz3HIyQq4X5nuaygrSGF+Fn4HsoJ+sjIDZGf6yc4IkJ0ZICvTT1ZGgOwMP1kZ7vdajuXcFBgikhSO47izlKwgk8vOv22kO0pLe5jmtjBNbV20tIVpbus6M3Z6dtLSHqa1w33t+RArgGh81tM7fAYqI+AjMx4gmcFA/NVHRtDvjgf9ZMS/MoM+Ms+895EZ8BMM+sgIuO+DAT8ZAR8Zgfj3QR9+nzNqZ0KeBIYxZgVwH1AKhIF7rLUP9LHd7cDnceusAz5trX09lbWKSPIF/D73MuEBPK0wHInS1hGmpSNCa3uYts4I7R0RWjvi38e/IjFobO6ktT1Me1c3HfHxrki0z5/bFYnSFYnS3Da04DkXx4FgwEfQ73Nf42ES9PsIBJz4q/t5IP4VDDg9vndD5/R7v9/9NX6/Q8DnIxj0sWyeQ0YSMinlgWGMyQTWAl+01t5vjFkAvGKMedNau73HdouBe4FLrbU74uHxK2PMTGttV6rrFpGRJRjwUZCbed5H4p7vcaXd0SgdXd20d0bo6OqOf0XoPPN9N11h97Uz3OMrPt4Zjp4ZC0ei8bFuIv08FjUWg65w9B0zpOGUuW439/71pQT9w7vOmBczjKsBrLX3x193GmPWAbcB23ts91FgnbV2R3y7B4wx3wKuBJ7ubyeO4+Ab5IUWp49hptOxzHTrOd36BfXcm9/vHkrKDw3vulvRaIxwJPrHIIl00xWOEu6OEg530xWJEo5/dUWiRCLuZ13hbiLdUcLdsTNjPV8j3TH3+9NfkRjdUXfcHXPfhyNRFswoISszAMP8GHMvAmMuUNVrrBJY3sd2W3qNVQELSCAwSkpCQz5OWFgYGtKvH43Sred06xfUswyeF4ERAnpf+tAeHx/Mdn2qq2sd0gyjsDBEQ0Mr0egwR/QIlW49p1u/oJ7Vc2KKi8+9SKUXgdECZPcaC8XHB7Ndn2KxGN3dg6rvjGg09o7jnmNduvWcbv2Cek4XyejZi9spdwJzeo3NA7b1sZ05/cYY4+Aepuq9nYiIpIAXgfECEDHG3AlgjFkCXAM81Gu7B4HrjTGL4u8/iTu7eClVhYqIyB+lPDCstWHgZuBuY0wVbjDcZa21xph7jDFfim+3C/gL4Bfx7T4G3GytjaS6ZhER8ejGPWvtVmB1H+Of7/X+YeDhVNUlIiLnpiUhRUQkIQoMERFJiBOLpdelZiIiMjiaYYiISEIUGCIikhAFhoiIJESBISIiCVFgiIhIQhQYIiKSEAWGiIgkRIEhIiIJUWCIiEhCPFl8cCQzxqwA7gNKgTBwj7X2AW+rSj5jzKeA7wJfsdZ+2+t6kskYczXwLaAA8AM/sNZ+19uqkssY8x7gG0Au7pOe/91a+z1vq0oNY0wh7vN1nrHW3uFxOUlljGkFjgM9V/W+0Vrb+7HYg6IZRg/GmExgLfBDa+0s4APAvT2eyTEmGWO+D7wL2ON1LclmjCkHfgN8wVo7F3gP8HVjzCpvK0ueeM+PAX8d7/kG4B+NMZd5W1nKfA/o9LqIZDPGBIEcYJW1dm6Pr2EJC1Bg9HY1gLX2/vjrTmAdcJuXRaXAw9baDwHNXheSAt3Ax6y1zwFYa/cBe4HFnlaVXDHgT621rwBYa/cDVcACT6tKAWPMe4FZuM/dGesK468NydqBDkmdbS7u/0g9VQLLPaglZay1L3tdQ6pYa2twZ5EAGGNm4j4i+BXPikoya+1J3FkVAMaYNcA04BnPikoBY0wR8K+4M6qx/o8+gCLcQ1H3G2MW486qfmSt/fFw7UAzjLOFgPZeY+3xcRljjDGTgd/inqfa4XU9yWaMud4YcwT4JfDn8dnVWPY93PNT1utCUqQTdyb1PWvtQtzHWv+TMeZPhmsHmmGcrQXI7jUWio/LGGKMWY77r+5/s9b+s9f1pIK19klgijFmHvBbY0ww/lTLMccYcyMwA7jD41JSxlp7CLizx/s3jTE/B/4E+PVw7EMzjLPtBOb0GpsHbPOgFkmSeFg8CfxNOoSFcd14+r21djfwBHCTd1Ul3a24gbHfGHMQ+BvgA8aYDV4WlUzGmAJjTO+/v3y4V3sOCwXG2V4AIsaYOwGMMUuAa4CHPK1Kho0xJgt4FPhLa+3jXteTIkXAz+PHtU9fZvou4A1Pq0oia+1HrbUTrbUV1toK3HMZj1lrx+zVcMBC4DVjjAEwxswCPkSPc3ZDpSfu9WKMWQr8ACgDOoCvjuW/WIwxftyZFcBU3MNv9cBaa+3nPSssSYwxt+Ee5+19ccMvrLVfTX1FqWGMuR34Iu59Jw7uDOMfrLXD9q/PkcwY81WgIg3uw7gL+HvcP+cu4F5r7X8M189XYIiISEJ0SEpERBKiwBARkYQoMEREJCEKDBERSYgCQ0REEqLAEBGRhCgwREQkIQoMkRQzxnzNGLPR6zpEBkqBIZJ6FwGbvC5CZKB0p7dIihhjArhLr2T2GG6w1hZ5VJLIgGiGIZI6UeDC+PfvBiYAs70rR2RgNMMQSSFjzDXA74A8a+2Yf860jC2aYYik1jJgl8JCRiMFhkhqLQPe9LoIkcFQYIik1mJgq9dFiAyGAkMktXzALGPMJGNMsdfFiAyEAkMktb6M+yztI8CPPK5FZEB0lZSIiCREMwwREUmIAkNERBKiwBARkYQoMEREJCEKDBERSYgCQ0REEqLAEBGRhCgwREQkIf8fBp9hhEMsww8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbad8886390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.integrate\n",
    "\n",
    "#my initial condition\n",
    "c_init = 1.0\n",
    "\n",
    "#where I want to evaluate my function\n",
    "t_points = np.linspace(0,5, 100)\n",
    "\n",
    "c_points = scipy.integrate.odeint(rxn_d, c_init, t_points)\n",
    "\n",
    "plt.plot(t_points, c_points)\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$C(t)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's solve the spring-mass system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEOCAYAAABM5Pr8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8XNd16PvfmRn03gFWgASwCRLsnao2RVlWcZGV6yKXJLZy7ev48/Lu9U3iPN8kzk3s58R5TtwSOXFsy3YkP6vYsmT1QkkkJbGBJEBiAyABgiBBovc+M/ePMwNCEEgMMOWcmVnfz4cfigcjztqcObNm77323obX60UIIYQIhsPqAIQQQkQ/SSZCCCGCJslECCFE0CSZCCGECJokEyGEEEFzWR2AVTo7BxdVxmYYBnl5aXR3DxMvlXDSZmlzrJI2L7zNBQUZxlzXpWeyQA6H+WI44uhfTtocH6TN8SFcbY6jf0IhhBDhIslECCFE0CSZCCGECJokEyGEEEGTZCKEECJoti0NVkr9EfBt4K+01t+6xmM+DXwFsx3dwJe01ocjF6UQQgiwac9EKfV94Dag/jqP2QB8B/g9rXUF8APgcaVUYmSiFEII4WfXnsnDWus3lFKvXucxnwSe1lrXAmitH1JKfR24FXh+vidYbJ31vz91mrfqrpCWksCygnSqV+Wyc20RORlJC//LooTDYbzj93ggbQaP18uZll6O6g4aL/TT0TfK5JSHlCQXxXmpVK3MYde6IpYVpFsZdlDkdQ4dWyYTrfUbATxsDXB01rVGYB0BJJO8vDQMY+H/mC2XB5mY8jAxOE7v4DinznXzq1eaeM/W5Xzy/VXkZiYv+O+MFtnZaVaHEHHx2OasrFQOnmzn58+eoa1j6F0/HxqdpKmtn6a2fn57oIVNlQV85s61lC/PtiDa0IjH1znUbbZlMglQGjA669qo7/q8uruHF9Uz+X8+vZULXaO0Xe6n8UI/NY1djIxP8cLbrbxec5H791Vy08aSRSUqu3I4DLKz0+jrG8bjiY8tJ+K1zR6Hk3/8xRFONnVPX19ZnMGG1XksLUgjKcHJ0Ogk5y8Pcryhi+6BMWoaOjnRsJ/37VzBfbeuJsFly9HzOcXr6xxMm3Nz5+6JRnMyGQJSZl1L812fl9frxe1e+JOmJLrYVV1Cz5IMbtm4lPEJN/tPXOK3B5oZHpviR0+fof58L5++Y01U3VSB8Hi8uN3xccP5xVObm9r6+e7jp+gbHAdgXWkO996ymrKSzHc99obqEj723gqON3by2P5zXO4Z4dm3WtGtvfzxvRuibtg3nl5nv1C3OZo/7eoA5f+DUsrAHPo6GckgkhKd3L59OX/3wC42rs4D4EDtZf7pVycYn1xEthLCAqdbevjmfx6jb3Cc5EQnD9y9lv/+0U1zJhI/h8Ngqyrka3+4gzt3rQSguX2Qr//sCO3dw5EKXdhENCeTnwN3KqXW+/78OcxeyWtWBJOZlsiX7tvAB24oBeDM+V6+99hJJqckoQh7O9HUxT/96iQTkx7ys5L5y9/fxu7q4oCHahNcDu67dTVfunc9iS4H3QPjfOuRGrr6Z49Ci1hm2G3bZaWUE7PXAbACM0H0AE/4rg1rrf/W99iPA18FEoF24L/5q7vms9gt6J1Og9zcdHp6hq7ZRXz2rVb+/1eaANiwOo8/vnc9Lmf05u1A2hxr4qXNTW39/P3Dx5hyeynITuEbX7yRBDyLbnNTWz//+MsaxifdFOWm8pVPbiEz1b7V+vHyOs8UbJuvtQW97ZJJpIQzmQA8dbCFx187B8Dercu4f1/l4gK1AbnhYrPNXf2j/O1PjzAwMklBdjJ/8amtlJfmB93muuYe/ulXJ3B7vJQWZ/Bn928hKcEZwshDJx5e59nClUyi9+uyzd29p5Tbty8H4KWjbRyqvWxxREJcNTYxxXcePcXAyCQpSU7+r/s2hqysfV1ZLg/csxYDs5T+4RcbQvL3CnuTZBJGv/ee1axZYdbe//TZelqvDFockRCmXzzfQFvnEIYB//UD1SzJD+2agx1VRXzwpjIAXjvRLl+m4oAkkzByOhx8/kPV5GYmMTHl4cEn62RCXljuSH0HB3wf7h+5ZTUbfFWIoXb37lLWluYA8NBzWiq8YpwkkzDLTE3kCx+sxjCgvXuEX7/RbHVIIo71DY3z0HMaALU8mzt2rgjbczkcBg/cs46stETGJ938+1Nn4mZhYDySZBIBq5dm8b4d5k377FutnLs0YHFEIh55vV5+8kw9Q6OTJCc6+ezdVTjCvFNDVloin727CoDm9gFeOtYW1ucT1pFkEiEfvqmMkrxUvF74j9+dYXLKY3VIIs4ca+jk5Flzm5T791WSnzV7A4nwqC7LY/e6YgAe33+O7v6xiDyviCxJJhGS4HLyh3dWYQCXuoZ56ah8QxORMz7h5pGXGgGoWpnDnuriiD7/x/aWk56SwPikm589r4nXJQmxTJJJBK1emsUtm5YA8OSBZvqHJyyOSMSLpw610D0wjtNhcP++yohvRJqRmsjHb6sA4OTZbk6d657n/xDRRpJJhH3o5lWkJLkYm3Dz+P6zVocj4sCV3hGee7sVgNu3Lw95GXCgdq0tonJZFgC/fLkJt0eGemOJJJMIy0xN5IM3mvX3b5xsp+WyTMaL8Hps/zmm3F5yMpK4x7d3nBUMw+Cje83eSXv3CK+daLcsFhF6kkws8N4tS83JeOBXr0jvRIRPy+UBjtR3APChG8tITrT21Imykkx2rS0C4Nevn2N0fMrSeEToSDKxgMvp4L5bVgPm7sL153stjkjEqsf2m/vDleSlsmd9ZCfdr+XeW1bhcjoYHJnk2bdarQ5HhIgkE4tsqsintDgDgCdePyfVLSLkzrT0UNfcA8C9N6/CuZijRcMgPyuF27YtA+DFoxcYGp20OCIRCvZ4d8UhwzD48M2rAGhs65++6YUIBa/XO71rdVlJBlsqCyyO6J3u2LmCxAQHo+NuXjh8wepwRAhIMrFQdVku5b7qlsdfk96JCJ36872c9e208OGbV0W8FHg+mamJ7N0ivZNYIsnEQoZh8OGbzN5Jy+VB6lqkdyJC46lD5wGzV7KuNNfiaOb2PumdxBRJJhZbsyKb1UvNc7afeVMmI0Xwzl7s54yvqOPu3aW265X4ze6dSGVXdJNkYjHDMLhz50rArOxqbpd1JyI4T/t6JUvz09hYkW9xNNf3vh0rcDnN3sn+mktWhyOCIMnEBjZW5FOSlwrAM2+etzgaEc3aOoeoaeoC4K7dK8O+K3CwMtMSucFXsvzCkQtMuWVVfLSSZGIDDsOYPlfiqO7kSs+IxRGJaOWfe8jPSmZ7VaHF0QTm9u3LMYDewXEO+xZYiugjycQmdq8rJicjCS/wvExGikUYGJngUN0VAPZtW26bdSXzKclLY5NvOO65t1qlqjFKRce7LQ64nA72bjUnIw/WXmZkTCYjxcLsP36RKbeH5EQnN24osTqcBfEfHtfaMTRdPCCiiyQTG7l54xISXA7GJ928cUo2wROBm3J7ePnYRQBu3FBCSpK1e3AtVMWyLMpKzKpGOesnOkkysZH0lITpTfBeOnpBzssWATt8poP+4QkM4DZfDzeaGIbB3q1LAahp6qKrf9TiiMRCSTKxGf9QV2ffGCflACERoBePmvNsmyryKcxJtTiaxdm+ppCM1AS8Xnj1uJQJRxtJJjazoiiDyuXZgHT3RWBaLg/Q3D4IwHujsFfil+BycvNG8yTS105cYnLKbXFEYiEkmdiQf5iirrmHjl4pExbX5/8WX5STQtXKHIujCc6tm5ZiGDA0OsnbZ6RMOJpIMrGhTRX5ZKUlArD/hHT3xbWNjk/x1mmzHPiWTUttv0hxPnlZyWwqN8uEXzl+0eJoxEJIMrEhl9MxXdp54GS7rAoW1/Rm3WXGJ924nI7pleTR7j1bzIn4c5cGaOscsjgaEShJJjblHzseGJmkprHL4miEHXm93ulv79vXFJCRmmhxRKGxtjSXvMwkAN44KSXy0UKSiU0VZKewrszcOnx/jXT3xbudax+grXMYgFs3L7U4mtBxGAY3rDd75gdrL0vPPEpIMrGxW3y9k7qWXjr6pO5evNOBU5cB83z38qVZFkcTWjeuL8HAnIiXnnl0kGRiY5sq8slMTQDMuRMh/Can3NMT7zesL7HtmSWLlZ+dQlWpWZn2urz3o4IkExtzOR3sWmdOqh6svYxHNsATPscauhgdn8IwzE1CY9FNG8yeeW1zNz0DYxZHI+YjycTm/GPH3QNjNLT2WRyNsIsDvr3bqsvyyMlIsjia8NhSmU9asguv92p7hX1JMrG55YXprChMB+BArdxQwjz3o66lByBmyoHnkuBysmut2b43TrVLz9zmJJlEgT3V5g11RHcyPiFbTMS7g7XteL2QmuRis82P5Q3WTRvNnnln3xhaeua2JskkCuxcV4zDMBifcHOsodPqcISFvF4vb/iquHauLSLB5bQ4ovBaUZTByqIMAA7KUJetSTKJAllpiaxfZa45kaGu+Hb20sD0sc7++bRYt9vXMz/a0MnEpPTM7cp2J+gopbYD3wXygUngG1rrh+Z4zAHg3IzLvVrr3RELNML2rC/hxNluzrT00jMwRm5mstUhCQv4J6JL8lIpK8mwOJrI2FlVyC9fbmRswk1NUxc7qoqsDknMwVY9E6VUEvAE8C9a63LgPuA7Sqn1sx6aDTRqrdfM+BWziQRgU3keqUkuvMChustWhyMsMDnlmd5J98YYXFtyLVnpSawtNXvmb/rOuBf2Y6tkAuwF0Fr/1Pd7HfA08PFZj8sG4mo2LsHlZIfvFMaDtZfxSmVL3Dl1rttcW4I5XxJP/CeQnjrXzdDopMXRiLnYbZhrDdA461oDsGXWtRygQCn1PLASuAj8ldb69UCfyDAMHItIpQ6H8Y7fI+mmDSW8evwi7d0jtHUOUeo7MzvcrGyzVezY5rfPmN/KK1dkU5CTEvK/345t9tteVcjPntNMTHk4qjtCdgiYndscLuFqs92SSRowexOqUd/1mc4DzwLfBC4BnwKeUUqt0VoHdDxhXl5aUMME2dmzQwq/7TlpFOam0tEzwolzvWxZtySiz29Fm61mlzaPjk9R02Qe47x3x0pyc9PD9lx2afNsO6tLeL3mIod1J/ftWxPSv9uubQ6nULfZbslkCJj9lSvNd32a1vo54LkZlx5SSv0pcBvwk0CeqLt7eNE9k+zsNPr6hvF4Ij/UtH1NAU8fPM/+Y23cs2dFRA5DsrrNVrBbmw/WXmZi0o3TYbBueRY9PaE/58NubZ5ta0Uer9dc5HRzDw3NneRnBd87s3ubwyHYNl/ri4zdkkkd8OVZ16qAkzMvKKWWAmitZ+7N7sCs/gqI1+vFHUSVocfjxe2O/Jtvx5oinj54nu6BMRpb+ylfFrndYq1qs5Xs0uZDtWbRxbqyXFKSXGGNyS5tnm1taS7pKQkMjU5y8NRl7tpdGrK/265tDqdQt9luE/CvAFNKqT8AUEptBG4HfjHrcb8PPKGUyvQ97kPAUuClyIVqjWUFaZTkpQJM7xorYtvQ6CR1zeb2KfE28T6Ty+lg+5pCAA7VXZEiFJuxVTLRWk8CHwQeUEo1Aj8HPqu11kqpbyilvup76Dcx15kcU0rVA/8TuEtrHfM1s4ZhsNNXZ3+4/gpujxwcFOuO1Hfg9nhJdDlifvuU+exaZ773L3UNc6FDjvS1E7sNc6G1rgH2zHH9KzP+ewr4v32/4s6OtUX8+o1mBkYm0a190zX4Ijb5e6Aby/NJTrTdLRtR5UuzyMtMontgnCO6gxVF8bFwMxrYqmciAlOcmzq9X5G/XFTEpp6BMRoumEuqdsXxEJefYRhs8w11HT7TIUNdNiLJJErtWGveUEd1p5yRHcMO13fgxdwhuHpVntXh2ML2NWZSvdI7KkNdNiLJJErt8N1Qw2NT1PomZ0XsOVJvbp+ypbKABJfcrgBlJRnk+famO+z79xHWk3dnlMrLSp4uC5ahrtjUMzDG2UsDAGxbU2BxNPZhGMZ0VdfhehnqsgtJJlFsh++GOtHUxeSUDHXFmqPaPLsmJclF1Uopsphpe5X53u/oHaX1igx12YEkkyi2VZk31Oi4m9MtMtQVa45ocwhnU3m+DHHNUlqcQX6WDHXZibxDo1hORhKrl5qbPfo/eERs6B0cp6mtH5AhrrnMHOo6IkNdtiDJJMptrTRvqJrGLqnqiiHHGjrxAkmJTqrLZIhrLv4S4Y4+GeqyA0kmUW6bMr+1Do9NoVvj6oiXmHb0HUNcsX3O+2LNHOp6u16KUKwmySTK5WensLLYXMAoQ12xoX94Au1bqOj/siDezTCM6Yl4GeqyniSTGOD/wDnW0Bk322jHsuMNnXi9kJjgkIWK89jmK0Lp7BuTBYwWk2QSA/xVXYMjk9Nbb4jo5e9hblidT1KCDHFdT2lxBrmZSYD5ZUpYR5JJDCjOTWVZgXlqmn9tgohOgyMT1J+XIa5AGYbBloqrPXNhHUkmMcLfOzna0IFHxo6j1vHGLjxeLwkuBxtWyxBXILb6km5b5zBXekYsjiZ+STKJEf5vsX1DE5y7OGBxNGKx/ENc61flxf1284GqWJZNRmoCIL0TK0kyiRFL8tMozjVPYJSqrug0PDbJmZZeQIa4FsLhMNhUbh4adlSSiWUkmcQIwzCmu/tHdaeUSUahmsYu3B4vLqfBxvL4PlFxofzv/XOXBugdHLc4mvgkySSG+MskuwfGaLk8aHE0YqH8QzRrS3NJSZIhroWoWplLcqJZ+SZDXdaQZBJDVhSlT68Ilqqu6DI+6abOdy7NlkoZ4lqoBJdjujcnycQakkxiiGEY072T441yQ0WT0y09TEx5MECGuBbJn4R1ax9Do5MWRxN/JJnEmM2V5gdRe/cI7d3DFkcjAnW8oQuA1UuzyEpLtDia6LR+VS4upwOP10tNY5fV4cQdSSYxZvWSrOkySbmhooPH46WmyXyt/F8GxMIlJ7qmd1iWoa7Ik2QSY2aWSR6Toa6o0HSxf3pYZnOFzJcEw1/VVdvcw+j4lMXRxBdJJjHI/4F07uIA/UNSJml3/vmtkrzU6bVCYnE2lufjMAym3B5Oneu2Opy4IskkBq0tzSExwYEXpodPhD15vd7p+RKp4gpeekoCakU2IENdkSbJJAYlJjhZX2bu63Rc5k1s7WLXMB19owBsqpD5klDwD3WdPNvN5JScPhopkkxilP+D6XSLjB3bmT/ZZ6UnUlaSaXE0scE/Zzg24aa+tdfiaOKHJJMYdXXs2Du9GE7Yz3HfUMzmigIchmFxNLEhNzOZshLz9NHjMtQVMZJMYlR6SgKVy7MAqeqyq54Z295sliGukPIXofi39BfhJ8kkhvlvqJNN3Uy5ZezYbvzFEcmJTtasyLE4mtjiL2boH57g3CU5kiESJJnEMP+33ZHxKTnO14b8QzDrV+WR4JJbMZRK8lIp8pVZy1BXZMg7OIblZ6ewojAduLpdh7CHkbFJ6lvNBC+r3kPPPM736saPciRD+EkyiXH+qq7jTXJD2cnJc924PV6cDoMNqySZhIN/qOtK7yiXuuU433Bb8KEJSqkKoBRIATqBk1pr2VHQprZUFvDkgRZ6BsZpvTLEyuIMq0MSXO0prlmRTWqynF0SDmVLMslKS6R/eILjDZ0szU+zOqSYFtC7WClVCnwR+ARQDMysYZxSSr0BPAj8SmstM702srwwnbzMZLoHxjjW0CnJxAYmp65u9bFZVr2HjcMw2FyRz6s1lzje2Mnde0qtDimmzTvMpZT6/4BTwCrgz4G1QBaQiJlY7gBeA/4WOKGU2ha2aMWCGb4bCmQ1vF3Ut/YyNuEGri6wE+HhT9bN7YP0DIxZHE1sC6RnkghUaK0vz/GzDuBl36+/Vkp9GKgAjoQuRBGszZUFvHi0jbbOITr6RinMTrE6pLjmry4qLc4gNzPZ4mhiW9XKHFKSnIyOuzne2MXercusDilmzdsz0Vr/sT+RKKW6lVK513nsE1rrh0MZoAhe5fIs0nzj8jVSJmkpj9fL8emzS2SIK9xcTgfrV/n3qZP3fjgttJorhzl6M0qpHKXUL0MTkgg1p8PBhtUy1GUHze0D9A9NALLqPVJmHuc7PCbH+YZLoBPwB4G3AS+wTinVq7We+aqkAB8JRUBKqe3Ad4F8YBL4htb6oTke92ngK5ht6Aa+pLU+HIoYYtGWynwO1V2moa2PwZEJMlLlaFgr+Ku4CrNTpLooQtavysPlNPepO9nUze7qYqtDikmB9kxeABRmFddLwKBS6phS6kdKqS8Bfwa0BRuMUioJeAL4F611OXAf8B2l1PpZj9sAfAf4Pa11BfAD4HGllHxCXsO6MvN8bK8XTjTJoUFW8Q+1bK7Mx5CNHSMiJclF1Urfcb4y1BU2ASUTrfVfaa3fD4wD64EPAI8AqcDngduBL4cgnr2+5/up7/c64Gng47Me90ngaa11re9xD2EmultDEENMSk50sa7U3P9Jxo6tcblnhHbf4jk5njey/LsM1J7rYWLSbXE0sWmhq6XStdZuoA54PgzxrAEaZ11rALbM8bijs641AusCjcswDByLWP/vcBjv+D2abF1TyImz3dQ19zDl8ZCU4Azo/4vmNi9WONp8wjfxnpFqngZot3/PWH6dt6oCfvasZnzSjb7QN70zRCy3+VrC1eYFJRNfIgmnNGB01rVR3/XFPO6a8vLSghpmyM6OvvHu92xfyY9/d4aJKQ/nO0fYVV2yoP8/GtscrFC2+eQ581yZXdUl5Ofbd/FoLL7OubnpqJU51J/vpball/fuLH3Hz2OxzfMJdZvnTSZKqWr/cFIAj00ASrXWs3sXgRrCnMyfKc13fTGPu6bu7uFF90yys9Po6xvG44m+va7Kl2bR2NbP/qMXqFwS2AdatLd5MULd5v6hcepbzGSydmU2PT0Bv1UjJtZf542r86g/38ubte18Ym85DocR822eS7Btzs1Nn/N6ID2T55RSrwM/BF7RWr/r2ZVSxZjzGH8MfJ13D1UFqo53z71UASfneJya8fwG5tDX7Mddk9frxR1EP8vj8eJ2R9+bb3NFAY1t/dQ0djEx6ca5gIwarW0ORqjafKyhCy+QmOCgakWOrf8dY/V13lSezy9fbmJwZBLd2kfl8uzpn8Vqm68n1G0O5JNEAU3AL4E+pdTLSqmfKaX+TSn1mFKqHrOS6y7gM1rrHwYRzyuYe339AYBSaiPm5P4vZj3u58CdM6q8PofZK3ktiOeOC/61DUOjkzS19VscTfw45lssuq40l8QA56pEaBXlprLEV459TBbvhlwgK+CHtNZfBZYBnwFqgHRgOTCCucFjtdb6PVrr/cEE41u78kHgAaVUI2bS+KzWWiulvqGU+qrvcaeBLwCP+B73KeCDWuupYJ4/HhTlpk6vb5AFjJExNjHF6ZZe4OoCOmGNzXLGSdgEPAGvtR5XSp3RWv86nAFprWuAPXNc/8qsPz8MyNYti7C5Mp+LXcMca+jko+8tl/UOYVZ7rocptwfDgA2r86wOJ65tqSzg6UPn6eofo61zmNIS+xZCRJuFTkG/4Vuh/i5KKTmUIUr41zj4bygRXv51PWp5tuw8YLHS4gxyMpIAOc431BaaTP5f4EWl1L6ZF5VSnwB0yKISYfWOG0oWMIbVlNszveOALFS03swjGWTeJLQWlEy01v+IeUjW40qpjyul9imljgM/An4TjgBF6BmGcfU4XzkbPqwaLvQxMm5O5cnGjvbg3625tWOIzr7Zy9XEYi14pYXW+ufA1zArrJ7G3ACyXGv930McmwijLb5vyeevDNLdL4cGhYs/WS8vTCdfzpGxBbU8m9Qkc1RehrpCZ0HJRJkewRzueg5zr67TWuuL4QhOhI9akU2K/4aSoa6w8Hq90xsLSq/EPlxOBxvLzUKIo5JMQmahPZNaYAlwk2/jx/cBf6mU+nrIIxNh5XI62Ljaf2iQDHWFw/krg/QOjgNSEmw3/vkr3dpH/9C4xdHEhoUmk49orW/WWh8C0FofBG4CPqmU+nHIoxNh5Z83kUODwuOYb4grLzOZ5YVzb0EhrFG96uqRDIdPX7E6nJiw0An4J+e4dhq4EdgZqqBEZPgPDfJ4zUODRGjJ2SX2lZzoorrMPOPkzdp2i6OJDYvY6vDdtNatmAlFRBE5NCh8OnpHuOhbw7NFSoJtafN0RWMn43LGSdBCkkwAtNY9ofq7ROTIoUHh4Z+HSkt2UbE8y+JoxFw2VuRjGDAx6ebUWemZBytkyUREp83l+RjA+KSb0+d7rQ4nZvhLTjeV5y9oZ2YROZmpiVQsM3cOlgWMwZN3eZzLSk9i1dJMAGpkqCskBoYnaLxo7si8Waq4bG2rMl+fmsYu3B6PxdFEN0kmYrpMsqaxK24OCAqnE01deL2Q6HKwzjfJK+zJX7I9PDZFQ2ufxdFEN0kmYnoicmBkkrOX5IyTYPnnS9aW5pIkZ5fYWmFOCqUlZs/8mKy3CookE0FJXholeamA7NUVrLGJKWqbzVoUf3GDsLdd1SWAWcotZ5wsniQTAVwd6jomN1RQ6pqvnl2ysVySSTTYvd5MJj0D47ReGbI4muglyUQAV79Fd/SOcql7xOJoopd/1XvFsmwy5eySqFC2JJP8rGRA9uoKhiQTAUBZSSZZ6eaHn+ykujhTbg8nz5rJZIts7Bg1DMOYnoiXTU8XT5KJAMBhGGz2DcvIDbU4jRf6GB4zzy7ZJCXBUWWLr0T4YucwV3qlZ74YkkzENP+aiOb2QXoG5IyThTri69EtK0inUM4uiSqVy7NIT0kApAhlsSSZiGlrVuSQnGiWssq29Avj8Xo5ps1ksm2N9EqijdNx9YwT2aducSSZiGkJLgcbfGecyPYSC9PU1k//8AQA21ShxdGIxfBvyHl2xmspAifJRLyD/4OwvrWXAbmhAnZEdwBQkpfKkvw0i6MRi7G2LJdElwMv5i4GYmEkmYh3WL86j8QE89Ag6e4HxuP1ctQ/xCW9kqiVlOCkepX0zBdLkol4h6QEJxt8N9TR+g6Lo4kOze0D08fzblsjySSa+bcWOt3Sw+j4lMXRRBdJJuJd/B+JAwaqAAAX/0lEQVSIZ873MTQqx/nOx98rKcxJYVmBDHFFs43l+TgMgym3d3pbHBEYSSbiXTasziPB5cDj9coCxnl4vV6O+Hpw21ShHM8b5dJTEqj0HWYm7/2FkWQi3iU50cV631DXYS1DXdfTemWIrn5zTY7/bAwR3fyr4U+c7WbKLWecBEqSiZjTNt8H45mWXoZlqOua/FVceZnJlBZnWByNCAX/pqej41PUt8rpo4GSZCLmtLE8H5fTwO3xygLGa/B6vRzxzZdsVQUyxBUj8rKSWVlkfjGQ1fCBk2Qi5pSS5KK6zDfUdUaGuuZysWuYKz3mPk5SxRVb/LtoH2volNNHAyTJRFyTf1uQ2uZuGeqag3/iPTs9kVVLMi2ORoSSf71Q//AEDRfkON9ASDIR17SpPB+nwyyTfPv0ZavDsR3/2RdbVSEOGeKKKUvy01hWkA7A27LeKiCSTMQ1pSYnsK4sF4ADJy5ZHI29tHcPc7FzGLharCBiy44qs3dyVHfg9khV13wkmYjr8pe7HtMdsiJ4hrdOXwEgKy2RimXZFkcjwsGfTAZHJqk/L0Nd85FkIq5rc0UBTofB5JSHGqnqAswqrrd9RQnb1xTicMgQVywqzEllpa/c+60zVyyOxv4kmYjrSk9JYG2pOdT15mm5ocBcqHjZV8W1c22RxdGIcPL3To7pTlnAOA9JJmJeu6vND8xTZ7tlry6ufkvNz0qWKq4Yt91X8j0yPkWd7NV1XZJMxLy2VBaQ6HLg9ninV3zHK4/Xy9u+ZLKjqkgWKsa4/KwUVi81vzC8LeutrkuSiZhXSpKLHeuKAXirLr6Hus5e7KdnwNxuXoa44sOONebrfLyxk8kpt8XR2JfL6gBmU0r9KfA5zETXCjygtT47x+O+D/wXoHvG5f/QWv99RAKNMzdvXsYbJy7RcKGPnoExcjOTrQ7JEv4qLnMdgmw3Hw+2rSnkkZcaGZtwc/Jsj2zoeQ226pkope4GvgS8V2tdDrwKPHyNh2cDf6+1XjPjlySSMNlWVUhqkgsv8dvdd3s806ved1bJdvPxIicjicrlZvn3W7J495rs1jP5NPCQ1rrN9+dvAf9LKVWptW6Y9dhsYNHF34Zh4FhEKvWXgcZTOajDYZDgcrK9qpD9NZd468wV7tqz0uqwwmqu1/n0+T4GRswChN3VxTidsfUeiNf39szfr2XP+mL0hT5ONHUzNjlFWnJCJMILi3C9znZLJmuAp/1/0FqPKKXagHXA7GSSA3xEKfVHQBawH/hzrXU3AcjLSwvqm2V2dvwNcezbVcr+mkucvzzI8KSX5UWxv+X6zNe55vlGAMqXZ1NVHrsbO8bje3u+Nt++ZxU/f76BySkPp1v7ed+u0sgEFkahfp0jnkyUUh8DvjfHj/p9v4/Ouj4KzNXql4ArwA+BZOBnwI+ADwUSR3f38KJ7JtnZafT1DcfNbqL+Ni/PSyE7PZG+oQmeO9jMvbessjq0sJn9Ok9Mujlw8iIA21UBPT1DFkcYevH83g6kzZvK8zlc38Hzb55nu29X4WgU7Oucm5s+5/WIJxOt9SPAI3P9TCl1AkiZdTkNeNedq7X+XzP+OKGU+t/AQaVUgtZ63sUQXq8XdxCFGR6PF7c7Pm64mbavKeKFIxc4VHeZD9xQGvPzBv7X+XB9B6PjbgzDXHsQy699PL63A2nzrnVFHK7voOFCH1e6R8jPnv1RFV1C/TrbagIeqAOU/w9KqQxgKXBq5oOUUg6l1Aal1MyBSwfg8f0SYbJrnVkm2dE7ytlLAxZHEzkHTpkTr+tX5ZGdnmRxNMIK61flkZ5ifuQcqpOJ+Nnslkx+AnxaKbXM9+c/Bw7MURrsBR4F/geAL6l8GfiN1loKwcOotDiDkrxUAN442W5xNJHROzjO6RZz9fOe6mKLoxFWcTkd7Kwyv0wdrLuC1xtfvbf52CqZaK2fB/4ReEkp1QhsBu73/1wpVa+UWqe19mLOjexTSjUAtZhDYV+wIOy4YhgGN24oAeBw/RXGJ2M/dx+qu4zXC6lJLjZXRO9YuQjebt+XiSs9IzS3D1ocjb3YrZoLrfW3gW9f42drZvz3aWBvpOISV+1eV8xjr55jdNzNsYZOdq+L3W/rXq+XA6fMHtiOqkISXE6LIxJWKivJoCg3lSs9IxyqvSx7s81gq56JiA7Z6UlUrzJ3Eo71oa5zlwZo7zZ3CL5hfYnF0QirGYbBHt+84VtnrshOwjNIMhGLcqPvg7X+fC9d/bOruWOHv1dSlJsq30IFwHRPfGh0khNNcsaPnyQTsSibKvJJT0nACxysjc3KlskpN2/6Nra8obo45sugRWDys1OoWpkDwH45znqaJBOxKC6ng12+XXMPnGrHE4OVLW/XXWF4bAoDqeIS73TLpiUA1J3roasvdnvmCyHJRCyafw6hs2+Mxguxd0b282+dB6CqNCdud0kWc9tcUTDdM389xucNAyXJRCzayuIMlheaWyvE2g3V2TfK8QZzh+CbNy6xOBphNwkux/S84esnL+H2yES8JBMRlKtrTjpi6kjf/ccv4fVCZloCWyrl/ArxbjdtNN/7fUMTnDorR/pKMhFBuaG6mESXg8kpDwdPxUbvZMrt4TXfxOpNG5bgcsptIt6tJC8N5TvnZH/NRYujsZ7cJSIoqckJ7PBNxL9Scykmtpioaeyif3gCuDrRKsRcbva9P06e66ZnYMziaKwlyUQE7T2blwLmFhP153stjiZ4/nLPTRUFFOWmWhyNsLNtqoC0ZBdeb+wv4J2PJBMRtLKSTEqLzYOyXjke3d399u5h6prN8e/37Y7t0yRF8BJczun9uvafiO+JeEkmIiRu9fVOjjd20Ts4bnE0i/fyUTMZ5mQksatatk8R8/P3zHsHxznWEL8r4iWZiJDYubaItGQXbo+XV463WR3OooyMTfFGrTlUsXfrMpl4FwEpyUujuszcq+6FIxcsjsY6creIkEhKcE5PRr56/BITUbg1/RsnLzE+4SbB5eDWzTLxLgJ327blADS19dNyOX4OjZtJkokImb1bluF0GAyNTnIwyk6i83i8vHTM7FHtWltERmqixRGJaFK9Kne6WOOFw9HZMw+WJBMRMrmZyWxbUwjAC4cvRFWZ8ImmLjr7zNJO/7dMIQLlMAxu22oeEPv2mSv0D0XvvOFiSTIRIXX7dvODuL17hNrm6FgV7PV6+Z1vH641K7Knt4gRYiFuWF9MSpI5b/ji0fjrnUgyESFVVpJJ+bIsAJ5587zF0QSm4UIfZy+a49x3SjmwWKTkRBfv3WJWdr18rI2RsSmLI4osSSYi5O7caX4g17f20dTWb3E08/vdm60ArChKZ11prsXRiGi2b9tyElwORsfdUVvVuFiSTETIbSzPY1mBOVT01KEWS2OZT+uVQU6d6wbgzl0r5QAsEZTMtERu8m1++sLhC1FZ1bhYkkxEyBmGwd17zN7JybPdnL88aHFE1/Y731BcYU4K21ShxdGIWHDHjhU4DIOBkUneiJHNTwMhyUSExTZVSLGvVPKpgy3WBnMNbZ1DHD5jnlny/p0rcDikVyKCl5+dws615heTpw+dZ3IqPnonkkxEWDgcBnf5JrOPNnTasnfyxGvn8AL5WcnTp0YKEQr33FCGwzDoHRznlePxcU68JBMRNjvXFlGUkwLAY/vPWhzNOzW3D3C80dxH6YM3lsnWKSKkinNTuWG9uQHk04daGB2P/couuYNE2LicDj5yy2oAapt7OGOj7ekff+0cACV5qexeV2xxNCIWfeCGMlxOg8GRSV6Mgz27JJmIsNqqCigrMbenf/TVJlusij/T0jO9zfyHb1olcyUiLPKykrl1k7nu5Nm3W21xrLXX66WzbxSPJ/T3oSQTEVaGYXDfreUANLcPcri+w9J43B4P//lSIwClxRlsVXK+uwifu/aUkpTgZHTcPd0bttKj+8/y5e8f5DevhX7YWZKJCLuqlTlUrzIXA/7y5SbGJqwbP371+CUudg4D8Il9lbKuRIRVVloi99xQCsD+4xctLURp6xziubfM4TanM/Tve0kmIiI+cVslLqdZ3fLkgRZLYhgcmeAJ37fDPdXFlC/NsiQOEV/2bVtOUU4KXuAXLzZYMtTr9Xr5zxca8Hi9FGQnc8eu0pA/hyQTERHFuancsXMFYK4MbuscingMj+0/y8j4FEmJTu67dXXEn1/EpwSXg4/fVgmY550crI388Qyvn2ynvrUPML/YJSY4Q/4ckkxExNy1u5T8rGTcHi8/e06HZRLwWmrPdfPaCXM18gduKCU7PSlizy3EhtV5bCrPB+CRlxrpi+AW9b2D4/zyZXOecJsqYEuY5gklmYiISUpw8ol95je0xrZ+nnu7NSLPOzI2yY+fqQfMXY392+QLEUmfvL2SlCQXw2NT/OSZ+ogMd3m95he30XE3acku7vfdf+EgyURE1Kby/OnFXI+/di4iE5IPv9hI7+A4LqeDz95VhdMhb3sRebmZyXzitgrA3LPu9ZPh37frjZPt1DSZi3M/treCrDD2yOWuEhH3idsqKcg2h7t++Ns6xsO4s+rB2nYO+MaoP3LLKpbkp4XtuYSYz57qYjZXmMNdD7/UyMUwzh22Xhnk5y80AOYw257q8C7OlWQiIi4lycUD96zDYRi0d4/wo6fP4AlDl7+5fYCfPKMBszx5nxzHKyxmGAafuWMNWWmJjE+4+e5jp8KymHFkbIofPFHL5JSH3MwkPnf32rCXwUsyEZYoX5rFvbesAuBIfcd0yW6o9A9P8L3HTzHl9pCXmcznP7hOVroLW8hMS+SL967H5TTo6Bvlwd/U4vZ4Qvb3T7k9PPhkHR19ozgdBl/4UDXpKQkh+/uvRZKJsMz7d66YsRneeV47EZrdVYdGJ/n2L2voHRwn0eXgSx9ZT0ZqYkj+biFCoXxpFp+6XQFQ19LLT56pD0l1o8fj5d+fOj194NvH9laweklk1lNJMhGW8Xf51fJsAH76TD0vHwvuqNOh0Um+9fBxWjuGMAz47N1rWVGUEYpwhQipmzYumR56PXDqsjncG0RC8Xi9PPSc5m3fGT137V7J3q3LQhJrICSZCEu5nA6+eO96ykoy8AI/f76B3x5oXlTZZFf/KH//n1cTyQN3r2X7Gjk9UdjXR/eW854t5maQh+ou8+CTiytImZh08+Bv6qZ79+/ZspR7b14V0ljn44roswVAKZUD/BC4DyjQWndd43EpwL8CNwJe4ADwea31aKRiFaGRnpLAlz+2me89fooz53t54vVmWjuG+MwdawIe6z3e2MmPf1fP0OjkdCLZJVvLC5tzGAaf3FeJ0zB48Wgbh+s7uNwzwgP3rGVZQXpAf8fFziEefLKONt+eczdvXML9Fuw7Z6tk4kskbwKPBPDw/w0UA1WYyeQp4GvAn4YtQBE2KUku/uT3NvBvT53hSH0HR3UnurWPe24o5eaNS0i6xvYPFzqG+PXr56YPukpJcvJfP7CODavzIxm+EItmGAYfv62CrPREHt9/jgsdQ/zNTw7z3i3LuGPnimvu1tA3NM6zb7Xy0tE23B4vBnDvLau4c9dKSzYwNexwvoSfL5ksAYaBZq7fM+kA7tdav+D78z7gp1rrJYE8V2fn4KIa7nQa5Oam09MzhNttn3+7cIpkm71eLy8fu8ijr56d7u6nJDnZsDqfspJMstISmZzycKV3hLrmHlpmLHqsWJbFZ++qojAnNeg45HWWNlvhTEsP//G7M3QPmNutOB0GVStzUCuyyc1MBqC7f4yGtj5ON/dOl9TnZSbzh3dVUbUyZ97nCLbNBQUZc2YqWyUTP6VUKddJJkqpXKAbKNVan/ddWwm0ALla63mP9OvqGvIuZiG0w2GQnZ1GX99wRPeWspIVbe4eGOPx/ec4VHsZ9zzPWZCdwodvLmN3dTGOEH0jk9dZ2myVsYkpnnmzleffvsDIPMf9pia5uH3Hcu7ctZKkxMA2bwy2zbm56XPeZBEf5lJKfQz43hw/6tdaB7qVq38Z88z5kdEZP5s3meTlpQXVFczOjr+V1JFsc25uOn9Wmk/PwBiv11zkVFMX5y8PMDA8gcvpoDA3lYrl2exaV8KmyoKwrSGR1zk+2K3Nn/1QNve/fy2Hats5cuYKzZcG6B0YAyAnM4nVS7PZuqaQXdUlJCct7mM81G2OeDLRWj9CYHMi1+PfgyBlxrW0WT+7ru7uYaRnEhir23xTdRE3VRdd8+d9fcMhf06r22wFabP92ryxLIeNZdceuhoZHmNkgW//EPRM5rxuqwn4QGmte5VS7YACzvsuVwEXtNZ9gfwdXq8XdxBbQnk8XluMsUaStDk+SJvjQ6jbHJXJxOcnwJeVUq8CBvBl4MdWBiSEEPHKVosWlVL3KaXqgZd8lw4ppeqVUjt8P69XSq3z/exvgFbgNFAL1AN/F+mYhRBC2KxnorV+FHj0Oj9fM+O/x4DPRSIuIYQQ12ernokQQojoJMlECCFE0CSZCCGECJokEyGEEEGz5XYqQgghoov0TIQQQgRNkokQQoigSTIRQggRNEkmQgghgibJRAghRNAkmQghhAiaJBMhhBBBk2QihBAiaJJMhBBCBM1WW9DbnVJqO/BdIB+YBL6htX7I2qjCSym1F/g6kAU4gR9orb9tbVSRoZTKBuqAF7TWv29xOGGllMoFHgR2YR4290Ot9d9YG1V4KaVuBv4B8709Bfyb1vqfrY0q9JRSfwR8G/grrfW3fNfygR8B1YAHeBL4n1prz2KfR3omAVJKJQFPAP+itS4H7gO+o5Rab21k4aOUKgZ+A/yF7yyZO4C/UUrttjayiPlnYNzqICLkx0AHsALYCuxTSlVaG1L4KKVSMd/bf+d7b98GfFUpdYe1kYWWUur7mG2rn/WjfwW6gXLM13sv8PlgnkuSSeD2Amitf+r7vQ54Gvi4lUGFmRv4lNb6JQCt9VmgCdhgaVQRoJS6G/NG+7nVsYSbUmoJcCfw11prr9b6itb6Jq11g9WxhdEKIBt4DkBrfRk4gflNPZY8rLX+L8Cg/4JSKgP4EObIildrPQB8H/hkME8kySRwa4DGWdcagHVzPDYmaK07tdZP+P+slFoNVAEHrIsq/JRSOcA/AX+IOQQQ6zZh9kr+QCl1Sil1Qin1BauDCrMmzPv3fgCl1CpgPVePDI8JWus35rhcgTmUeXbGtaA/yySZBC4NGJ11bdR3PeYppZYBv8X8NlNrdTxh9s+Yc0Pa6kAiJAcoBMa11uuBTwHfVErtszas8NFaTwG/D/yDUqoL84vi97TWxy0NLDLSMF/rmV+Ugv4sk2QSuCEgZda1NN/1mKaU2gIcAn6qtf6a1fGEk1LqHmAVZs8kXvQBXuB7AFrrk8BTwPutDCqclFIlmJPO92ut84Ei4B6l1JesjSwihoAkpdTMz/+gP8skmQSuDpg9IVkFnLQglojxJZLfAX+itf6m1fFEwEcxk8k5pVQL8CfAfUqpQ1YGFWZNQALv/GbqxaxwilU3AANa62cBtNZdmD3v91kaVWQ0YM6Hls+4FvRnmSSTwL0CTCml/gBAKbURuB34haVRhZFSKhn4FfBFrfVjVscTCVrrT2qtl2itS7XWpZg9lEe11jFbweYbzjsA/AWAUqoUc0L+aQvDCrfTwFJfub+/umsfUGNpVBGgtR4GHgW+opQyfCXw/w2zom/R5KTFBVBKbQJ+ABQAY5jVLzH7IauU+jhmNdPswoNHtNZ/HfmIIk8p9ddAaRysMykGfobZ+x4GvqO1/ldrowovpdT9wJ8BSZhfrF8E/ofWesTSwEJEKeXEHFEBs3ptCOjBXOLwD8C/YRZfuIFHMNehLDohSDIRQggRNBnmEkIIETRJJkIIIYImyUQIIUTQJJkIIYQImiQTIYQQQZNkIoQQImiSTIQQQgRNkokQNqKU+ppS6k2r4xBioSSZCGEvO4G3rQ5CiIWSFfBC2IBSyoVvN9cZl/u01jkWhSTEgkjPRAh78ADbfP+9DyjBPMRIiKggPRMhbEIpdTvmOSIZWut4OXtexAjpmQhhH5uB05JIRDSSZCKEfWwG4uHYWBGDJJkIYR8biIPDmURskmQihH04gHKl1FKlVK7VwQixEJJMhLCPvwQ+AFwAHrQ4FiEWRKq5hBBCBE16JkIIIYImyUQIIUTQJJkIIYQImiQTIYQQQZNkIoQQImiSTIQQQgRNkokQQoigSTIRQggRtP8D6FY9jxQYAH4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbad889ce80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#This is just pasted from above\n",
    "def ideal_spring(x,t,k=1):\n",
    "    #This will get x as a vector!\n",
    "    return np.array([x[1], x[0] * -k])\n",
    "\n",
    "x_init = [1.0, 0] #The position is at 1.0 and has no velocity\n",
    "t_points = np.linspace(0,10, 250)\n",
    "\n",
    "x_points = scipy.integrate.odeint(ideal_spring, x_init, t_points)\n",
    "\n",
    "plt.plot(t_points, x_points[:,0]) #Only plot x[0], the actual x-position\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$x(t)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "What about the dampened system?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\frac{dx_1}{dt} = x_2 $$\n",
    "\n",
    "\n",
    "$$\\frac{dx_2}{dt} = -k x_1 - b x_2$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEOCAYAAABM5Pr8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJztnXeYY2d97z+SRpreZ3Z3tvd3q9e7XttgjMEVDIQSCAnFvgETLinc5xZSIJWSkEJCCyQkIcEEB4MBh+JgMI5xB7Mu2/fd3mdnpxdpRjOSzv3jFB1pVM6RtLtTfp/n2WdHR2denVfS6KtfDxiGgSAIgiCUQ/BKX4AgCIIw+xExEQRBEMpGxEQQBEEoGxETQRAEoWxETARBEISyqbrSF3Cl6O0dLSmNLRAI0N5eT39/lPmSCSd7lj3PVWTP/vfc2dkYyHVcLBOfBIPmixGcR8+c7Hl+IHueH1yqPc+jp1AQBEG4VIiYCIIgCGUjYiIIgiCUjYiJIAiCUDYiJoIgCELZzNjUYKXU+4FPA3+qtf5UnnPuBj6MuY9+4INa619cvqsUBEEQYIZaJkqpLwC3AYcKnHMV8DngV7TW64AvAt9RSkUuz1UKgiAINjPVMvm61voppdRPC5zzbuAhrfU+AK31V5VSfwG8GvhxsQcoNc/6n79/gFM9Y3zkrh3UVc/Up6+yBIOBjP/nA7Ln+YHsuXLMyE9DrfVTHk7bADyfdewIsBkPYtLeXk8g4P/JfPFwL9GJBMcvjHHT9qW+f38209JSf6Uv4bIje54fyJ7LZ0aKiUfqgfGsY+PW8aL090dLskzam2uIToxx7MwgW1a0+F9gFhIMBmhpqWdoKEoqNT9aTsieZc9zlXL33NbWkPP4bBaTMaA261i9dbwohmGQTPp/0IVtdZzuGeNCf4xkcn68+WxSKUP2PA+QPc8PKr3nGRmA98h+QNk3lFIBTNfXnkv5oIta6wDoGYhdyocRBEGYVcxmMfka8Dql1Fbr9vswrZInLuWDLmwzjaELIiaCIAgOM87NpZQKYVodAMuBTUqp9wEPWseiWutPaK0PKKV+E7jfSgfuBt6ktU5cyutb2GZaJtGJBGPjUzTUhi/lwwmCIMwKZpyYaK2TmO4qL+d+Hfj6pb2iTBZZYgKmq6thSfPlfHhBEIQZyWx2c10RGuvC1NeYGiyuLkEQBBMRE58EAgG6Oszs477hiSt8NYIgCDMDEZMSaKgzO7aMxy9peEYQBGHWIGJSArVWG5WJSRETQRAEEDEpiboaW0xKqHoUBEGYg4iYlIBtmYzHRUwEQRBAxKQk6mrM2hJxcwmCIJiImJSAWCaCIAiZiJiUgATgBUEQMhExKQEJwAuCIGQiYlICYpkIgiBkImJSArZlkkgaTCVSV/hqBEEQrjwiJiVQ65r9LtaJIAiCiElJZIqJxE0EQRBETErArjMBERNBEAQQMSkJt2UizR4FQRBETEqiJhIiYP0slokgCIKISUkEAgFqqkOABOAFQRBAxKRkaiNSuCgIgmAjYlIijmUiMRNBEAQRk1KpsSyTcbFMBEEQRExKpRItVZKpFM/s66ZnMFapyxIEQbgiiJiUSE3EDsCXbpl89WHNv/zgIF/+wcFKXZYgCMIVQcSkRNIzTUqzTLr7ozy5pxuAo+eGMQyjYtcmCIJwuRExKZFyLZMf/+JMxu3ohATyBUGYvYiYlEh12BSTyanSxOTi4HjG7f7hibKvSRAE4UohYlIi4SrzqZsssQX9aGwq43b/iIiJIAizFxGTEok4lkmpYjKZcbtPLBNBEGYxIiYlErEsk6mEfzdXyjCmWyYiJoIgzGJETErEsUxKcHPFJhKkrOytJZ31gLi5BEGY3YiYlIhtmZQSgHe7uFYtagLEMhEEYXYjYlIiYcfN5d8yGYm6xKSrERDLRBCE2Y2ISYm43Vx+Cw7teEkkHKSr3XRzjY1PEZc+X4IgzFJETErEdnMBJJL+rJMRy83VVBehoS49Ajg6MZXvVwRBEGY0IiYlEnaJid8gvG2ZNNZFqHONAI5JO3tBEGYpIiYlYru5wH+tiW2ZNNaFZZ68IAhzAhGTEomE3ZaJv1jHaDTt5qqJhAhYA+Vj0p9LEIRZiohJiUSq0pbJlG/LxHJz1YcJBAKOq0ssE0EQZisiJiUSKStmkrZMIN3OXmImgiDMVkRMSiTscnP5balit5uvrzEzuepqxDIRBGF2I2JSIuFQECvU4dsysUf92jNRbDeXxEwEQZitiJiUSCAQSLeh9xEzSaUM5/yaalNMxM0lCMJsR8SkDNIzTby7ueKuXl41YVNEKhGAf+zFczz07ElSKRn/KwjC5aeq+CmXF6XUtcDngQ5gCvik1vqrOc55GjjuOjyotX75ZbtQzFqT6ETCV38u95hf281VW1Oem2tgZIJ//5EGYCQ6xTtuW1fSOoIgCKUyo8REKVUNPAj8odb6XqXUZuBppdSLWuu9rlNbgCNa681X5EItSukcbMdLAKqzYialWiZne6POz4/sOsO1GxewdklzSWsJgiCUwowSE+BWAK31vdb/+5VSDwHvALLFZKicBwoEAgRLcPIFgwHnf7sKPpEyCIUChX7NYcrVx6u+topQKEB9rZnVFYsnPK/j5sJANOP20bNDqOUtvtfJh3vP8wXZ8/xA9lw5ZpqYbACOZB07DOzIOtYKdCqlfgysAM4Bf6q1ftLrA7W31xMIlP5ktrTUU2el9lZVhWhra/D0e+cG0q3muxY2EwmH6GwzOwfHp5Ke13HTP5o5Ang0Xto6xWhpqa/4mjMd2fP8QPZcPjNNTOqB8axj49ZxN6eAh4G/As4DdwE/VEpt0Fqf9fJA/f3Rki2TlpZ6hoaiBAJmsHt4dIKBgTFPv3+xbxSAUDDA6EiMQCCAkTTdZGOxKc/ruDlxfjjj9vme0ZLWyYd7z/MlwC97lj3PVcrdc74vqjNNTMaA2qxj9dZxB631j4AfuQ59VSn1e8BtwFe8PJBhGCTLGB+SShmEQ6aba2IySTLp7UWx039rIiFSKQCDatdslPhkkqqQd5UzDIPzVsykq72O7v4YvcMTnq/HD6mUcUnWncnInucHsufymWmpwfuB9VnHNgJ73AeUUkuUUkuyzgtiZn9dNiLOtEU/AXjzXDv4DukKePBfazISnXR+Z+vqdsAcAex3YJcgCEI5zDQxeQxIKKXeA6CU2gbcAdyXdd6vAw8qpZqs894MLAEevXyXmm6p4qcC3p6mWBNJC0hGG3qf6cHn+9LB9y2r2szHmEo6LVsEQRAuBzPKzaW1nlJKvQn4olLqI8AEcI/WWiulPglEtdafwIyVdAAvKKUSQD/weq31hct5vXbnYD9dgx3LxDUPpZwBWQOjcQAaasMsXZD2ZfYPT9BQG873a4IgCBVlRokJgNb6JeCGHMc/7Po5Afwf698Vw6kz8VW0mNmXCzItE79iEh23pzaGaaqPUBUKkEga9A1PsGJRo6+1BEEQSmWmublmFbaby0/MJO3mSotJVSjotGaJT/octGWJSUNtmGAgQFtTDQD9w9lJcYIgCJcOEZMysN1c8RLcXG4xgbTby6+YjLnEBKDdEpO+kYm8vyMIglBpREzKoJxsLncAHqDasnLiPlqzgFmbAmkxaW2sBswsL0EQhMuFiEkZhMuImVRnWyaWuEyU6uaqM8XEFhXbYhEEQbgciJiUgd2bq5SuwdPdXP6bRoIrAF9rjgCur4CYdPdH+btvvMTe4/0lryEIwvxCxKQMSukabLuxpru5rGp6n2IymhUzabAKIKNliMm3fnqMfScG+PcfaVJS/CgIggdETMogXFVJy8QO5nsXE8Mw0jETy83lWCYlFi2OxCbZc8y0SPqGJzhypqzmzIIgzBNETMog7ATgU57bl+QVk4j/bK7xeMKxHBprM2Mm8ckkiaS/2fQAzx3oIelq/vb0vstaByoIwixFxKQMbDExIOMDOB+GYaQD8OHyLZNRlyurIUtMoLS4yQuHe4G02L2ge6XPlyAIRRExKQNbTMCbq2sykcL+XK6EZeIWi+xsruz7vdIzaBY7XrO+EzAr8iUzTBCEYoiYlEE45E9M4hnz33MH4P1YJna8JBgIOC1Z6mvSYuI3CD+VSDFk9frauLLVOX5xUKrpBUEojIhJGfi1TNxCEQlnPvU1ZVgmDbVVBK2pkZFw0JmH4teiGBiZwHZorepqcqyli0MiJoIgFEbEpAwyxMRDsNtd3JgdM4mUEjNxMrkizrFAIEBDrZUe7DOjq2843YKlvamGzmZzTlmvWCaCIBRBxKQMIn5jJhmWSaaY1JRQZ+JYJjWZLrNSq+D7rOaQzfURczZ9i9nnq1csE0EQiiBiUga+A/AuoXD/LqQD8H4KIMetdvV1NZlzS0oXE9My6bBEZEGraZlIzEQQhGKImJRBppgUFwG34EwTE9symUz6qFmxZqNUZ1o5dhC+ZDGx3FudLZaYlGmZ9A9PcPLCSFlrCIIws5lxw7FmE1V+s7msVvXhqqATMLexxcQwIJFMOdX1hRiP527NYlfB+83m6rNEo6PZtExsMRkcifvqjOzm4Z+f5jtPHCORNPiju3eyenFTSesIgjCzEcukDAKBgCMoXsTE/kCOVE1/2t1dhL12DrYtk9qsmpWy3VxZYmJQmqvr1IVRvvnYURJJ09J67mCP7zUEQZgdiJiUidNSxUc2V3bwHTKzu7xmdI3nac1SipgkkimGrRko9oCt1oZq5/6BEoZt6ay+Xi8clmp6QZiriJiUibs/VzHs4Houy8QtCF4nNzp9vqoz3Vx1VnbXuI958m7haao3U42rIyHnuoZG4p7Xsjl2bhiA9iZTlPqGJzjbG/W9jiAIMx8RkzIJ+3BzFbJM3Me8Fi5OWGKRbZnY1fB2TMUL7smMja66FVtYBsf8WybHz5ticsuOpc4EyD3H+nyvIwjCzEfEpEzsSvbKWiZeYybmebVZAfja6nQBZDLlzcrJ1TQSzJoTwGmz4pXB0Tj9ljWzZkkza5Y0A3CuTywTQZiLiJiUiWOZlBkzqQoFCQXNDC8vlkkqZaQHbWWlBtdVp8XAq3UyGjMtk9rqqoy0Zccy8enmOn7eTAUOBQOsXNRIV1sdAN39MV/rCIIwOxAxKZOwj2mL9jnZNSY2flqquDO+8lkmYHb99YLdmqWxLrMA0rFMxvyJybm+MQC62uuJhEN0tZticmEgVnIQ3jAMegZinq0tQRAuH1JnUiaVyuYC09U1Hk94FJO0SGTHTOpcAflxj/25bMukyRUvAZdlMuovZmKnGdtV9F3t9YBpdQ2NTToxFD9854njPPTsKRa21vLe129k3dIW32sIgnBpEMukTKpKyOaqLmaZeHBzjbstk+psy8QlJpWyTHzGTLILIBe21Tr3dff7j5uMxxP8ZNdZwJy58qXv7Zc0Y0GYQYiYlIkdM0l4EhOrAj6fZeKj2eNEPL9lEq5Kx1/KFxPTghgajTsjgr2QXQBZE6mizUoRLiVu8sy+CxkW28BInJMXRn2vIwjCpUHEpEz81JkUqoAHqA57j79MFBi0FQgEnFoT7zET083VmMfNlUwZxMa9rZVMpRiwAvZ2ny/ACcJfKEFMntx9HoAbtixyBOrFI72+1xEE4dIgYlImfmImcSdmkkdMLFHw0k7FjplEwkGCwcC0+21Xl1cxGXEsk0wxsd1cAENRb66uQZcVY3cgBljUZsZNLgz4c3NNTCY4c9EM6F+3cQE7rJHCLx6WmhVBmCn4DsArpdYBK4FaoBfYo7Wet8UDEashoyfLxHJzRfI0cbQtEy8BeDvlNzuTyyZduOhNTMYcyyTTzdVUn749Ep2kyxKEQvQNpYP1thUB0G793O8zzfh0z5gzAXLloiYiVSF+/IsznOuLMjQWp6XBfzAf4Nj5YSanUmxc0Vr8ZEEQCuJJTJRSK4HfBt4JLALcX4UTSqmngC8BD2it51Xepq92KrabK0/MpNrH6F6n/Xwk91p1PsQkkUw5Uxmzs7nCVSHqaqqITSQYHpvM9evTsOMlDbXhDBecLSYDoxMYhkEgMN2iysUpKzbS2lhNU32E6nCIAGYDytM9YyWJydN7u/nyQwcBeMet67j92mW+1xAEIU1RN5dS6u+AvcBq4A+ATUAzEMEUltcCTwCfAHYrpXZesqudgTh1Jj4C8PljJt7rTNJNHsu3TNx9ubItE0gLzEjMq5hkZnLZtFnpwJNTKV8jhe1A+8pFjYApuousupXTPf6D8Kd7RvnX/zro3P76o0c4eGrQ9zqCIKTxYplEgHVa6ws57rsI/Lf178+UUm8B1gG7KneJMxt/vblsyyRfzMS/ZVJbXdgyiXmogLczuWB6zMQ8FubCAIxGvXUhzs7ksmlrSt8eGJnIaNtSiFOWYKywxARg2YIGuvtjnLZiKX54Zt8FDMPsjlxTHeJcb5Qnd58Xd5cglEFRMdFa/479s1KqH1NYBvKc+2AFr21W4K9rcLGYiY8K+DyDsWwcy2SiuACMZfTlmr5eo0/LxK6Wb23MFJPm+gihYIBkyqB/ZILlCxtz/XoG8ckk3VY/r5UuMVm+sJHnDl7kjE/LxDAMdumLALx8yyKa6yPc98hhXjzSR3wqmTEKQBAE7/jN5molhwAppVqVUt+ozCXNLqo8ZnMZhlHUMqnx1U4l98heG9ti8WKZxCzBiYSDOSc82q4vtwVTCDu20tKQaeUEgwGn8n3AYxC+ZzDmBN+XdDQ4x5cvNH++ODjuq9X+ie5R57F3qk52blhAIGA+53uP9Xtex83x8yP89X+8wP/67JP84JmTJa0hCLMdT2KilHpGKfUZzJjnZqVUtn+iFnhrpS9uNmBbJokiY20TSQO75i/fSN5IxEfMpIhl4icAb8cv7Nnx2dgxk1Gflklzw3SXWZsjJt7as/RalfRVoUBGC5ZlC0wrxQDO9np3db101KxN6WypYdmCBprrI2xYbrq3Xjjsv25lPJ7gc9/ew6HTQ4yNT/GdJ46z69BF3+sIwmzHa2rwI8B1mFlcjwKTSqkDwIvAS8Ba4OwlucIZjteYiXuGelHLxE/MJE82l58AvD0rPp+YNNZ7n9w4lUgH1+3qeTdtzTVwdpgBj+1ZLjptWWoz6mma6yM01oUZjU3R3R/z3Kfr2Dmzm/HmVe1ONtmW1W0cPDXIkbNDhX41Jw89e4qR6KSTXQZw78OHuGpNe96sPUGYi3gSE631nwIopcaBncAS4GpgO/ABTAvnQ5foGmc09gdGMTeXe3pidQViJvlG9tr4mbY4NmGLSe63gx0o9xIzcQ/Zym2Z2LUmXi2TzIaRbrra6hiNDXuuqE+lDE50m2KyuqvJOW4LUf9InP7hCSeFuRgTkwl+susMAHe+bAU3b1/C7//js0QnEuzSF7lhS5enddyMxCa594eHGI8n2Lyqjddct5yqkNQWCzMfv0WLDVrrJLAf+PEluJ5Zh1fLZNJlmYSLZHMlkgaJZKrgh0i+kb02tmUymUgVXStqtUmpz5NdZbu5ouMJkqkUoWD+tYbdYlI/XUzsEb6DXsVk0BSKzubpYrKovY7DZ4e5MOBNTLr7o87ztnpxWkxWLmokXBVkKpHiyNkh2psXeVpvz7F+JhMpgoEAd1y3jKa6CNvWtvPikT4ef+m8bzEZjk7yV/e94Ozn0OkhhsYmedft632tIwhXAl9feSwhEVw4MZOkQSqVvxHilMsyKZbNBcX7c+Ub2Wvj7hxcrKVKrIhl4q49GSvSn2vYipeEgoGcqb+tVnrw4OhkwefLxrZMOnNYJnZ7lm6PYnLcskpqq9N1KmAOJrMtlSNnhz2tBfC8NmMsanmLI7ivunqxs45XkbN58Injzu8ssvqYPfr82ZJjMCnDYN+Jfr771AkOnhyQLsvCJcVL0eIWr4sppcJWu5V5g3vQVSFXV9xDzKTaJQzF+nPlG9lrU+ejDb0TgM9jmTS6LIxiQfghyzJpbojkrHBvt8QkZRhFB24lUynHHdbZMt31ZH/g9g2Nk/DQG+2ENf1x5aImglnXtm6ZOVbYa9xkcirJHiv7a+eGBc7xLavaHRF96Yj33mEXB2M8tacbgLffvJaPv+86Niw33W8PPnncV8dmMLsafOaB3fzdN3bz3adO8Df3v8Tnv73X0/MkCKXgxTL5kVLqfqXULUqpnP0vlFKLlFIfAo4AN1f0Cmc4GWJSwNXl1zIpFDcpNLLXprbGh5iMe7dMiqUH25ZJLhcX4LShh+Lpwf0jcZKW9bKgJbebC8yOxnbWVyFOWJX0bheXzZrF6Rn17sFj+Th8Zsh5DXas63COB4MBrlrTDsCeY97F5OGfnyZlGLQ0RLhlxxJCwSBvv2UtYLbsf0H7yzT7j0cOs++4WQ5mZ8G9dLSPr/34sK91bPqGx/mXHxzgY1/5BZ95YDeHpGOAkIWXmInCbKPyDSCilHoeOAdMAG3AZsxsrieB/6G1fvwSXeuMJBzyJiZ2zCQQMNNcc+G2TCan8q9VaGSvjZ9pi9GJwtlckaoQtdUhxuPJopaJHTPJlcllX1d1JER8MsnA6ARmZ57cuAWiI4eYdDTXOEWQF/pjzjTHXKRSBuet4sdlCxqm3b/KcnMZhtkLbNOqtrxrAegzpgWztLOe5qzeYFetaeeZfRc4cnaY2MQUdXmeV5tEMsVzB01X1u07lzlJHSsXNbFlVRv7Tgzw8HOnMyygQhw8NchPXzJb9t95/XJ+5ea1fPO/j/Lwc6d5Yvd5dqzv4Ko1HUVWSXPg5AB//529Ge+7Pcf6efONq3jjjas8r2MzMZlg16FejnePUF9TxfZ1nTkFXphdeKmAHwP+SCn1ceBO4CZgFdCJ2TX4S8APtdaHKnFBSqlrgc8DHcAU8Emt9VdznHc38GFrD/3AB7XWv6jENfjBq5vLXf2er8Gh2zIp9O240Mhem6pQ0AkqFytcLObmAmiqr2Y8HiuaHmwXLObK5AJz1kpbYzXd/bGiGV29g+POWrkq06tCQTpbarkwECsan+gZjDlin0tMmuojtDfV0D8ywYnu4mJy6LT5zVwtm96CZcuqNkfk9p0Y4LqNCwuute/4gBPXun5T5rl3XLeMfScGOH5+hLMXx1ia49rdpAyDb/73UcBMLHjrq9cA8Lab13Ds/DBHzg5z3yOH2bC81VPq8vm+KF940BSShtowr96+mL3HBjjVM8p/PnWCcDjIndevKLqOzeEzQ3zpe/sZdKWGP/TsKa5Rnbznzo1OFqJX+obGeWpvN2cujhEKBVm1qJFXbO1y5vCUQiKZIpk0Mr7cCcXx/MppreNKqYNa6/+8VBejlKoGHgT+UGt9r1JqM/C0UupFrfVe13lXAZ8DbtRa77OE5TtKqTVaa2+VdRXCq5urWPU7mA0g7XqFeAHLpNDIXje11VVMJSaJxfMLQCKZcr5x5nNzgfmB3jMQK+7mihZ2c4EZN+nujxV1c9mWSS4Xl82itjouDMSKBuHP9ppWSbgqmDPNGGDV4iZLTEYKrhWfTHKy23SZqeXT61vqasKsW9rModND7D7aX1RMfn6wB4D1S5sz+pcBbFrZ5ojcE7vP884imV0v6F6nl9mv3rLWiQ0FAwHefYfio//2C3qHJnj0+bPc+bLCIjCVSPEP393HeDxJY12YP7p7J50ttfzSDav44oN72X2sn2/99Biru5pQy4v3Ndt/YoDPPLCbZMogFAywflkLA6NxegZiPK97Od8X5f++/WpPqdlTiRQPPnmcHz93JiOetOvQRR588gRvunEld16/Iuesn1xEx6f4wTOneF5fpMf6EtPaWM2WVW28evsSx3L1Ss9gjIMnB+kZjJFMGjTWR1i1qJENK1pLSvUeG5+if3iCyUSS2uoqFrbW5i1+vlL4TQ1+Sin1ulwWgFKqSmvtva9Fbm4F0Frfa/2/Xyn1EPAOzM7FNu8GHtJa77PO+6pS6i+AV3OZU5bdYlJodG+xjsFgfmuPWC6gQjGTQiN73dRWVzESnXSq5XMRc7nA8rm5wLRMwEMA3mmlkr8tvB03KVYFbxcsdhYSk/Y6OEpRy8QerrW4vT5vavOqrkZ2HbpYVEyOnh92Yjnrl+Uulty2toNDp4fYe7yfVMrI+6EWn0w6EyOv2zRddIKBAK/c1sV/PnmCZ/df4FduXpP3Q8QwDH7481MAbF7VNu0DftmCBl519WIee/EcDz17ilduW1yw2eZDz57kXG+UAPBbb97ivA7hqiAfeNMWPv7VXZzvi/KP393Pn733uoJfIE5dGOXvH9xLMmWwoLWW3/nlrSztbCBlGDz2wjnuf/QI3f0x/vK+5/nQO7azsLUu71rDY3E+9+29zuvUVBdm65p2EkmDPcf6GI8n+fbjx9l9tJ973rCx4FqJZIrvP3mc+3500EmRtxkcjfPknm6e3NPNVWvaecsrV2c0G8211vO6l5/sOsOx87nfQ3XVVVyjOnnZ5kWo5S3TEkHca+nTQ7x0tI89x/qcrEabUDDAsgUNrF3SzMYVrajlrXmtukQyxZmLYxw7N8yx8yMMjsZ57xu30NVS2hygfPgVk78EfqKUepvW+hH7oFLqncDHgTVlXs8GzCC+m8PAjhznPZ917Ahm/MaTmAQCAQqUS+TF/lCw/3cHwJNGilCeeEgiZU9ZDOU9B8wq+PhkkqlEMu957pqV+tpw3g8q29KYmEzkXWvc5TJrrA/nPC8YDDhug7HxqbxrpQzDKVpsbazOe549yndgNF7wubAtk4VttXnPW9xhxkl6BmIF1zrfZ4rJsoUNec9bu8SM3/QNTzA2PkUb5Hxuj1jxksUd9bQ25f6D3L6+g2/891HGxqc42TOSt0J/z/E+JqfMWpXrNy3MeW2vunox333qBNGJBC8e6ePlW3LXwRw6NcQJy2J6ww0rcq715ptW8fS+bmLxBA///BS/ems6+dL93j7dM8pDz5rCdPt1y6a5/epqq/jgW7fyZ//6C4ajk/zz9/fzu+/YnvP5ujg4zmce2E18MklLQ4Tff9d25z0QwqzR6Wqv43Pf2kP/SJy//NoL/N47t+d06Z3uGeXT39ztWLV3vmw5v3zTasdlF5tIcP+jR3j8pfMcPTfMn/3rL3jn7eu46erFGR/chmGw+2i/I2Jgeg1uu2Ypmyw35Ykgez5eAAAgAElEQVTzIzy5p5vu/hh7jvWz51g/16hOXnv9ctYtbXbc1b1D4zz+0nme3N2dkaFYX1PFsoUNRKpC9A2Pc74vRiyecASqrama6zctZP3SFloaq0kkU3T3x9h/wnysQl8CkymDkxdGOXlhlJ88f5ZAwIz7rVjYSFN9hEQqxcDwBBeHxjl9YWyaC/4XBy7w5htX5l2/FHyJidb6b5VSPZgupfcDfcBfY364/0MFrqceyE7LGbeOl3JeXtrb6z0PZ8pFS4v5UEnXi1RTU01bW26fdqjKfKpra8J5zwHTRTIcnSQUrsp7XtVZ81tPdSRER0f+b0pNlnVgBIJ51+oZTlsayxa35A0W20HmialU3rWGx9LZV8uXtOQ9b1mX+aE9OBrPe45hGE4r+9VLW/Oet36lmTk1GpsiXBPJ2UIf4Fyf+YGhVrbnXWt7fQ3BAKQMuDgyyfIl6dfZzVHrW+e29Z1512ptraero57uvij67AjXX7U053kvHDHTi69e38nKZbnjNG1tDVyzYSG7DvbwzP4eXn/T2pznPfKdfQCsWdrMK7Yvy/n+bmtr4C2vXss3HjnMI7vO8rbbNkyr4WlsrOXer+wimTJY2FbHb7z5qpzFsW1tDfzO26/mb+97ngMnB3l411ne/dqNGecMj8X59Dd/xnB0ktrqKj76/htYvWR60sWr2hpob6vnY1/+OcPRSf7yvhf4k/e9jA0r0s/JT58/wxe/vZvxeJJIOMT/fecOXnHV4sxrAj5017XcdM0FPv/NlxgajfNv/3WIp/Ze4DUvW8niTvM1+clzpzl4Mt38/Jady7j7dRtpdxXH3rgD3vW6zTy1+xz/8SPNud4xnte9PK97aayLsLCtltHYFD1ZVvHOjQv5pRtXs21dByGXS+viYIynXjrH4y+c4/j5YQZG4vzwZ6f5IaenPR82qxc3c+3mhWxb18nSzgaqIyFGY1McPTvE4VOD7Dvex9EzQ6QMs+Ho8TwWEUBLYzUbVrSycWU7r7thZd6C51LxvZrW+mtKqUXAfUAC+DfgDVrrcxW4njHMppFu6q3jpZyXl/7+aMmWSUtLPUNDUafoLhgIkDIM+gejDAzkdskMj5jaFwrAwED+y7QzvQaHYnnP6+03j9eEQwXXClvfEgeH86/VfXHE2cN4dIKJ2PQ4RjAYcFwYA8Pjedc665otEkgm855nW+Mj0Uku9AznDASPxiYdF1xdOJh3rYZI+gPz4NFe1i6d/kE1Hk84f/DtDeGCz1lXRz3neqPsOXyRnRsXZrzOYNaXHLaC76sWNhRca+vqNrr7ovxsbzdveNnyafdHx6fYZcVLrlnfUXCtGzabYrLnaB/7j/RMy1w72zvmrHXHzqUMDuafpH3zti7+6+kTjMam+Of/3MMH3rQZSL+3v/HjQxy1rK//8VpFLDpBLM9yV61s4eYdS3jshXN885HDNERC3GQVbg6Nxfm7+3dzvi9KKBjgg2/dSktt/vfs4tYafu+d2/nU/S8yGpvi9z7/JNdvWsiitjr2nxhwCkqb6yP8n7dvY9XiprxrrV3UwMfvuY57f3iIXbqXI2eGOHLmpWnnqeUtfOCXt9HRGCaV5z27ZUULH3/ftTy7r4f/+tkpzvVGGY1NZrh8G+vCvPKqLl61fYlT/zQ8nCkyVcCrt3Xx6m1dnL04xjP7LrDv+ABnLo45cZ+6mipWL27i6rUd7FjfmRE/MhIJJhIJwsDGpU1sXNrEm16xgujEFPrUEIdOD3JhIEZ0PEEgYE42bW+qYfnCBtYsaaajucbyyASoqa6a9t72St4vuX4WUUop4KPA24AfATcCByokJGC2acnu8bUR2JPjPOW6rgCmdZR9Xl4MwyBZRj1/KmWQTJovRLgqSHwqSXwy6RzLxg5yR6qCec+BdEZXLJ5/LftDtiYSKriWHU+Jjifynmf/QdTVVGF64nKfZ2dnjcam8q7lDqg31IbzntfiShvuHZpw/vjcXOhPG57tTTV516qrDlNfU0V0IsG53mjOQOnpnvQHxOL2+oLP2apFTZzrjXL8nPnB5X6dAQ6fGSZh3V63pLngWletbufHz53hzMUxLg6MTwssP3fwIsmUQVUoyNVrOwqutWVVG62N1QyOxnn452e4+zUq4/6HnjFdUh3NNexY31lwrUhViDe+YhX3PXKYZ/ddYKfqZPu6TgDO9Izy7Z8eA+DGq7rYsLy14FoAv3bLWk52j3Kie4QvP3SQAycH6Wyp4fHd5xkeM5tg3vP6jZ7WWrGwkd9/5w6+8J299AyO87P9PRn3b17Zyj1v2ERLQ3XRteprwvzWW7Zy4OQAjz5/lv0nBphMpKiJhFDLTBG8el0H7e2NDAyMFVkvwMs3L+JlmxZytjfKsfPDjEYnqQ6HWLaggXXLWpzAerHrAuhqr+etr1rDW1+1hmTKbI4aDgWpiWRme3pZqyZcxba1HWxbWzjdO/vvO/u9XS5+LZN9wLPAK7XWzyqlbgC+r5RaqLX+SAWu5zHMmfLv0Vr/m1JqG3AH8MdZ530NeEYptdXK8nofplXyRAWuwTe2mBRKDbYzvYplYFRb2V6F2qkUG9lrU+tMW8yfF1GsL5eN7TIbG5/KO7/d9hc31IYLZqy4W8kPjOQWk4tD5re66nAo5yhhN4va6zh2biRvEN62mJrqI0VTRld1NfLU3m6Od4/kbD+iLatkUVvdtPqSbNYva6EmEmJiMsmeY33cvCPT1fXzA+YH5bY17QWz8sBMg75t51IeeOwYT+/t5s2vXOW0cDnbO8az+81BqK+5bnnB3mk2N29fwnMHezhydpgv/+Agv/3LVTQ3RPjMA3uYTKRoro/wq7fkdqdlE64K8f9+dRufeWAPR88NO9diX/f/fOMmrlHeamQAlnY28LF7rufxl86x51g/o7EpOltreeVVXWxZ1ebbPb1pZRubVrZhGIYzAM1ew+9agYAZ+M6VXl4qoWDQeS1nM37F5K1a6+/ZN7TWzyilXgk8rJTq0lq/p5yL0VpPKaXeBHxRKfURzMLIe7TWWin1SSCqtf6E1vqAUuo3gfuVUhGgG3hTBbLJSsLLtEVbHAqlBgNUR+ygeYFsriIje228dA62CxYbiuT3226uZMogFk/kzPwacQoWC/9hRCyBGI1N5a01sWtMOltqi/7BL2ozxaS7P7cv5ow172RZZ/GQ2iqreG40NsXFwXGyd2lXfm/IkRKcTVUoyOZVbTyve3npaH+GmAyPxZ1alezakny8atsSvv/0SSYmk/znE8e5+7UbMAyDBx47hmGYVslN2xYXXwjTpXXPGzbxiXt3MTY+xd98/UXnvkg4yG+/ZWvB7L5s6mrC/O47tvPknvM8taebqUSK1YubeP3LV7CgQDZVPsJVQW7buYzbdi7z/bv5CAQCRb+ACaXjNwD/vRzHDiilbgQersQFaa1fAm7IcfzDWbe/Dny9Eo9ZLl46B8cThUf22nixTIqN7LXxMtPES8EikPEtfDQ2lfODZqhIwaKbtqYaRmNTDOapNSnUej6bxVb8wK5wz+acZZks6Sz+bXJpZ4NT7HngRD/bVqXTaycmE07K56aVhYsabbav6+B53cv+EwMMjsYdq+zJPd0YhumKtNuvFKOuporXv3wF3378OD996TwbVrRyvi/K3uNmEP+Xb1qdkapejAUttXzkrmv47AO7ndqKBa21vPd1G3LGnooRrgpyy46l3LIjd7KBMLepyKAErfVpzPjJvMRuKV+4N1flLZN8fblsvExbtC2TYpXHTR6aPXopWLSxJy7ms0zSNSbFC9hskbg4OD6tPscwDM705m+jkk1VKOikCO8+ktkPS58eIpkyCAAbVhQv0gO4Ri2grrqKlGHw5B6zxUkimeLRF8xZcjdu7fI1ROs11y13Rhb/43f3872nTwJw3cYFni0cN4va6vjEb1zPn/z6Tn7/Xdv5pw/fxkaPQikIbio2dUdrPVD8rLmJY5kUaqeSSNeZFMK2TAoVLRYb2WuTjpnkX8uJmRRxadRWVzmZZvmq4NOtVIoXQ9ndg/NNXPRS/W6z1HJfGUy3TgZG4o6YLvVgmQBsWmkKxe4jfRlxkwMnTbfU8kWNBYv93FSHQ9yw1awLefyl88Qnkzyz74ITmL51p79v8VWhIL/zlq2O4IHZC+y9r9tYcqp7KBhk5aImNq1sy0hlFQQ/iAOxAqRjJvk/tJ12KkXcEF6mLRYb2Wtji0kimWIqkcwZ/I8WmWViEwgEaKyLMDgaz2uZ2O3nW7xYJraY5LBMphJJhiyRKVT9btPaWE1ddRWx+PSMLru1SCgYYHGHN9/9ppVtfPvx4/QNjdMzME5nSy2GNRvEvN+bVWJz8/YlPPr8WQZH43zoi087Xwa2r+8sWJ2dj46WWv7g3Ts4dWGU+pqqkmISglBp5GtIBfAWgPdomUSKi0mxkb02dRkDsnKv57Sf9/BN286qytfsccRyczV5ipnYLVXi07KmeocmnATGXEOxsgkEAiyxrJOzvZl1AnbLDTMW4s2dtGJho+P2s+MRp3vGnErpbT467kI6DRTMGFXKMGhuiHDXHaVPUAwGAqzqahIhEWYMIiYVwEsA3g6oFwuQ1tiWScGYSeGRvTZeZpo4AXgP3Vrt6vJcbq74VNL5xt2Sp/28G9syiU8lnWuwsV1cgUDaHVYM24V1LktM7IrgVT5anAeDAbauNoPiT+w+j2EYPLPPTHftaK5hXQnB6TuvX87NO5bQ0VzD2iXN/K+3XuXJHSgIswVxc1UAT5aJdV+uVupuIl7cXEVG9toUm7aYMoyis0zcNFmWSS4317CrJ5GXbC63SAyMTGTEIOzge3tTjecOq3bc5Exv1KmDSRkGJy9YYtKVv+1MLm69Zgk/P9DD6Z4x9hzr5+cHTDG5YcuikmITgUCAu+5QxU8UhFmKWCYVwBGTgvNMPFomEe+WSb7BWDbuOpRYjgFZE/EktofJm5srv2Xink9RqGOwTXN9hJDV7iW7FX2vh27B2ay04iQj0UkuWmmuPQMxx1ry20J8/bIWx5r57Lf2MBKbIhQMcMPWLl/rCMJ8QcSkAnixTKY815mY908mUjn75ngZ2WsTCgadGEyuKnjbKgGvbi7bMpkuJnaNSU0kVLSaG6w+UJbo9A5n9uzsddU8eGX5wgZHPO1iQNvFVR0OObUoXgkEArz7zo24jZB33r7eU3aZIMxHREwqgB3YzScmhmG43FzF6kwKz4H3MrLXjS0S0RxB80wx8eDmsrK0RnK4uWzLxItVYmPPcM9ug2IPuvIjJqFgkPVWm/dDp80mhS8dMWewr1va7HlIkpvrNi3iY/dcx8s2L+TXbl3HzduX+F5DEOYLIiYVoJhl4j5ezM3ljqnkqoL3MrLXjS0SbuGwcQe+vYxLtau3h8cmp1lNdl8ud9+tYnRZYtLtqg2JTyUdy2RJh7/+R3Yh4aFTg8Qnk04m1g7V6WsdN8sXNvL+X9rMHddWrq2HIMxFREwqQLFsrknXca+pwZDbMvE6stfGsUxyxExsa6UmEvIU6LatjpRhTLNObDHxY5nYbdTdI3cv9MectOClHnppudlgTRYcjk7yrcePMZlIEQjAjnWli4kgCN4QMakAxSwTt4XhtWgRcrdUKdUyiRWwTLw29HNbHe6JcoBTZNjS6L376WLLMhkeS88usetEaqtDvqwcMKcoLrdapjz6vNmuZP3SlqKdggVBKB8RkwpQLJvLl2WS4eaavt6Eq/jQSwfUOidmkt8y8RJ8B2ioCzsZWEOjmZbJoO3m8hUzSVse3QOmq+uc5fJa3OF/EmYwEOAu15yPQADe+IqVvtYQBKE0REwqQLF2KhmWiQ8xmZjKkc5rWSaRcNBTUNlO+c0VM4l57BhsEwwEaLFqSAZdlolhGE42lx83V1Nd2BGybmusrt1by2+8xGbNkmbeeds6Nq1s5Q/v2ilNCwXhMiFFixXAV8ykiJsrGAw4LdDjk9PXs+smvGRyQdrqyFVnMuaxL5ebloZq+kfijlsLTHeZvXc/rqlAIOAMtrJnkdgV7Et8xkvcVHoOhiAIxRHLpAL4iZl4mTeRbvaY3zLxEi8Bc2gR5LZMxmLe+3LZ2JaHO2Yy5LNg0c2yBWZl+uGzQ/QNjdNvFTAur+AkO0EQLj0iJhXAqTPJFzOxYh+hYMBT1lRaTHJYJh5H9toUyuYaHTddU8XG4rpJi0k6ZmILSwBvrVTcbF1tuqGOnxvhiT3dgNkGZs0S//2vBEG4coiYVAC3ZZJrbrjTfr5IwaJNoZYqXkf22tiZWhOTSRJZYmdXsjfWehcAO1vLbZn0Dptt5JsbIp57adlsWtlGuCqIAfzgmZOAOZ/D7zqCIFxZ5C+2AthiYhjmjPRsbMvEawv0Qs0evY7stamvdbehz7RObDdXKZaJuxdXj1UnUspsjupwiI1ZUwuvXuevxbsgCFceEZMK4I6D5IqbTHkcjGXjxTIp1pfLps5VQ+JuqZJIphxxsRs4eqHFCrCPjU85e7UbKy5sK61v1bUbFjg/11ZXsWWVt5nogiDMHCSbqwKEQ5liUpsVg7ZjH8Xaz9sUmrbodWSvjTtTy53R5R5w5ccy6WhOt46/ODTOko56egbtXlqlDWp6+ZZF1NVUEZtIsHpxk6fWLoIgzCzkr7YCeLVMvGRyQbqlSqEK+GIje21qq6sIYM5Hd2d0uTv/+rFMOptrqQoFSSRTdPdF6Wqrc1rGL/TRmNFNMBBgu7Q8EYRZjbi5KkCGmOTI6LLrTIoVLNrYnYVzNXr0OrLXJhgIpKvgXZaJe8BVg4/U4GAw4DRoPNcXZWB0gkTSjBOVEjMRBGFuIGJSAYpZJra7ymvMpDpsfvhPFGhBX2xkrxunc/D4dMuktjrk2WKyWdJhFhSe74vSM5ieReJlXrsgCHMTEZMKECkiJpN+YyYRc73c2Vz+ihaBgpaJn7Rgmy6XmFy0MrlaG6s9708QhLmHiEkFyLRMpguA3zoTJwCfM2bir50KuIZaRdOurdES0oJt7KmFFwZinOgeBUqPlwiCMDcQMakAVaEibi5LAMrN5vIzsteN3ZzRXWg4Om6LiX/LxO6blUwZPLXXrFrfvEoaKgrCfEbEpAIEAuk2KTndXH4D8HnqTPyO7LVprp/eAsV2czWUYJl0ttRkDPEC2KkW5DlbEIT5gIhJhSg008QJwHtup5I7AO93MJaNbZkMR12WSRlurlAwyOtetsK53VwfYWGbZHIJwnxGxKRCFOocbKf4enVz1bosE/esdb8je23sFijDY5NO77ByAvAAd16/nBULzY6/r7lueUlrCIIwd5CixQphZ3Tlqg2xs7kiHntz1boqwMcnE+lmjfHSLJNmS0ySKYOx8SkaasNO/MRvl1+bqlCQ333Hdo6cHWLraml/IgjzHbFMKoQdD5ksVGfi0c1VV527BYqd2hvAX51Ji0swhsYmicUTTluWzubSs7DqaqrYtrbD08RHQRDmNmKZVIjClok/N5e7OaNbTGJWO5S6miqCPuaj26nBAMNj8QzXWbur15YgCEKpiJhUiEiBgVZ+s7nqXGm/4/HplonfRohVoSCNdWFGY1MMjU06llJVKFiym0sQBMGNuLkqhG11ZFsmqZThBOW9urnCVSEn1dg9gyRtmfjPwLLTg4ejcXqHzGFW7c01viwcQRCEfIiYVAhbKCazKuDdt/20G7Gtj1wxk/oSWrSnCxcn6bcmI3aKi0sQhAohbq4KYWdqTWa5udy3vWZzgRmEH4lOZri5Yo6bqwTLxFUFb1tKHSImgiBUCBGTCmG3jc9ugeK+7dXNBek6klhGzMR0c5VimdhZW+f7oo5rq6NF+mkJglAZREwqRL7UYHcMpVw3V6zEADzAmqXNAHT3x5xjYpkIglApJGZSISJ5mjO6xcVrNheka01i8fQMknTMxL+ba83ipmnBdkkLFgShUoiYVIh80xHdzRqrfbi5bOvDLi6EtLCUYpnURKpYvrDBud1UF2b5goYCvyEIguAdEZMKkTcAb2VzBQKZreqL4cRMJipjmQAsc4nH7dcuI+wjIUAQBKEQIiYVIm9q8FS6YDHgo6ajLisAn0imHCunFMsEYOcGs018VSjIzduXlrSGIAhCLiQAXyGcAHyWZWLHUKp9zlnPDsC7s7pKyeYC2Lq6nd968xYWtdWVLEiCIAi5mHGfKEqp3wPeh2k1nQZ+Q2t9LMd5XwDeDvS7Dv+r1vqvL8uFZmG7uaYF4J0mj/5cSrZlYteZuLO6SqkzsbGtE0EQhEoyo8REKfUG4IPAy7XWZ5VSfwJ8Hbgux+ktwF9rrf/mcl5jPqoj6XkmKcNwMqfsXl1+0oIhs87EMAynxgRKt0wEQRAuFTPtU+lu4Kta67PW7U8Bf6yUWq+1Ppx1bgswVOoDBQIBgiVEjOx269lt193DqlIpg7AtLsm0ZRIKeY+Z2ON0DcOc3miP7A1Y913Onlr59jyXkT3PD2TPlWOmickG4CH7htY6ppQ6C2wGssWkFXirUur9QDPwOPAHWut+PNDeXu8rIJ5NS0t9xu2OmMsN1VDjDKQKhkyLpKEuQlub91Tcrsl07KW6tpqAtU5dbZiO9saSr7scsvc8H5A9zw9kz+Vz2cVEKfVrwN/nuGvY+n886/g4kGvXjwI9wD8BNcC/A18G3uzlOvr7oyVbJi0t9QwNRTPmgkzE0vPVe3pHSE6arUpGRs2mikEMBgbGPD9OwlWseLZ7iIt95u/WVYd8rVMJ8u15LiN7lj3PVcrdc74vxZddTLTW9wP357pPKbUbyG4YVQ9M+/TUWv+x6+akUurjwDNKqbDWeir7/GwMwyA5fY6VZ1Ipg2Qy/UJUuZQpNpEk2WDeZ89tD1cFM84vRr01ACtlGPQNTdA/YopSU33E1zqVJHvP8wHZ8/xA9lw+M83NtR9Q9g2lVCOwBNjrPkkpFQS2AAddwhEEUta/y061aya7uwre75RFm1AwSEtjhIGROAMjaTFpb5IWKIIgzDxmWtHiV4C7lVJ2Rd0fAE/nSA02gG8B/w9AKRUGPgR8V2tdhr1ROhFXHUmmmPibsuimzRKOgdE4A8MiJoIgzFxmlJhorX8M/C3wqFLqCLAdeJd9v1LqkFJqs9bawIyN3K6UOgzsw3SF/eYVuGzAdGPZuJs7xp06E/9PtS0c/W7LRJozCoIwA5lpbi601p8GPp3nvg2unw8At16u6ypGIBAgEg4yOZWqiJsLoK3RzAjrG5pgcHQSEMtEEISZyYyyTGY7uarg7QB8TcS/btturmPnhkkZZqBMxEQQhJmIiEkFSbehT7u57HYotdWlxExMy8SdbyFuLkEQZiIiJhUk3ewxbZlMTNpi4t8yybZC6qqrSlpHEAThUiNiUkGcaYtWAD5lGExYw61KEYG2LDERq0QQhJmKiEkFsdvM25ZJfDLpuKhqS4iZ1NdUZWSBSbxEEISZiohJBcmeaTLumkFSSswkEAiglrU6tzcsbynzCgVBEC4N4oCvII6by7JMMsWktKf6g2/dypmLY0TCIRa315V/kYIgCJcAEZMKUmu1VLGD7nZaMJSWGgzmiN1VXU3lX5wgCMIlRNxcFaTWGloVtaYiTrgsk5oS3FyCIAizBRGTClJvjdPNntteEwld1mFWgiAIlxsRkwpiz22PWSN27emIUhsiCMJcR8SkgtTVpOe2g7v6XcREEIS5jYhJBXHEZCKBYRhpMYlIvEQQhLmNiEkFsd1cyZTB5FSKcav6vUYsE0EQ5jgiJhXEDsCD6eoaL6MvlyAIwmxCxKSC2G4ugOjElLi5BEGYN4iYVBC3mMQmEk6diVgmgiDMdURMKkh1OF1PYrq5JDVYEIT5gYhJBQkEAq6MLnFzCYIwfxAxqTDu9GCpMxEEYb4gYlJh0lXw4uYSBGH+IGJSYeoty2Q4OkncEhN3YF4QBGEuImJSYWqtWpPTF0edYzJuVxCEuY6ISYWx3VynLowBEEDG7QqCMPcRMakwtpsrkTRH97Y2VVMVkqdZEIS5jXzKVZjs+EhHc+0VuhJBEITLh4hJhVnYmjmnvVPiJYIgzANETCrM+mUtGbc7WsQyEQRh7iNiUmGa6iMs7qh3bneIZSIIwjxAxOQSsHxhg/Nzp1gmgiDMA0RMLgErFzY6P4tlIgjCfEBKsy8BN129mF/oi7Q21tDaWH2lL0cQBOGSI2JyCaiJVPGHd+280pchCIJw2RA3lyAIglA2IiaCIAhC2YiYCIIgCGUjYiIIgiCUjYiJIAiCUDYiJoIgCELZiJgIgiAIZSNiIgiCIJRNwDCMK30NgiAIwixHLBNBEAShbERMBEEQhLIRMREEQRDKRsREEARBKBsRE0EQBKFsREwEQRCEshExEQRBEMpGxEQQBEEoGxETQRAEoWxkbK8PlFLXAp8HOoAp4JNa669e2auqLEqpW4G/AJqBEPBFrfWnlVIdwJeBLUAK+B7wu1rr1BW72AqjlGoB9gOPaK1/fS7vWSnVBnwJeBkQAP5Ja/2xOb7nm4C/wXxvJ4B/1lp/VilVC/wjcCNgAE8DH9Baj1+xiy0TpdT7gU8Df6q1/pR1LO9rq5QKYj43b7KW2A/co7Xu8/qYYpl4RClVDTwI/IPWei3wNuBzSqmtV/bKKodSahHwXeAjWusNwGuBjymlXo75x9YPrAWuAW4FPnClrvUS8Vkg7ro9l/f8b8BFYDnm3m5XSq1nju5ZKVWH+d7+c+u9fRvwR0qp1wIfBxYBG61/i4CPXqlrLRel1Bcw93co665Cr+1vWbevAdYBfcAX/TyuiIl3bgXQWt9r/b8feAh4x5W8qAqTBO7SWj8KoLU+BhwFrgPejGmJGVrrEeALwLuv2JVWGKXUGzD/yL5m3W5kju5ZKbUYeB3wZ9beerTWrwS6maN7xhTNFuBHAFrrC8BuzG/pdwOf0lpPaq2ngE8xu/f8da3124FR+4CH9/PdwBe01sNaawP4S+AtSiclfIQAAAQASURBVKl6rw8qYuKdDcCRrGOHgc1X4FouCVrrXq31g/ZtpdQazG9qL2K6Qo65Tp8ze1dKtQKfAd6Laf6D+e1sru75akyr5D1Kqb1Kqd1Kqd9kbu/5KOZe3gWglFoNbAUeAzqt+2wOA13W+2LWobV+KsfhYq/tBjKfg2OY+rDe6+OKmHinHsj2oY5bx+ccSqmlwPeBT2L6keNZfvO5tPfPYsaGtOtYPXN3z63AAsz9bQXuAv4KeD1zdM9a6wTw68DfKKX6ML8Y/j2mqELm37b986zft4ti7+eMzzfrvDg+ngMRE++MAbVZx+qt43MKpdQO4FngXq31RzH3WG0F6WzmxN6VUr8ErMa0TNzM2T0DQ5hfEP4eQGu9B/gBcAtzdM9KqS7MgPO7tNYdwELglzAFBjL/tu0P0Fm/bxfF3s8Zn29KqRBQjY/nQMTEO/uZbvJtBPZcgWu5ZFhC8l/A/9Za/5V1+DBmPGWt69S5svdfxRST40qpk8D/xkyu+Ffm7p6PAmEyv3UawC7m7p5fAYxorR8GsLKUvg9cjxkrUq5zNwJntNZDl/0qLx3F/ob3k/kcKMyMN7e1XhARE+88BiSUUu8BUEptA+4A7ruiV1VBlFI1wAPAb2utv20f11pHgW8BH1ZKBawU2t/CzAia1Wit3621Xqy1Xqm1XolpoXxLa72dubtnjZn++hEApdRKzID8D5ijewYOAEus9H47u+t24CXgK8CHlFIRK2vzQ8yNPTt4+Bv+CvDbSqlmpVQA+DBwv5/0aKkz8YjWekop9Sbgi0qpjwATmHnYnpV7FvAWYCXw50qpP3cdvx/4HeCfMb/VJq1jX7nM13e5mct7fhvw70qpU0AU+LDW+nGl1F7m4J611geUUvcAX7YEIwj8BLOmKoXp8juAaaE9Avx5vrVmMpZ7ar91czmwSSn1PsyyhkLv53/BtNB3YQbqdwH/089jy9heQRAEoWzEzSUIgiCUjYiJIAiCUDYiJoIgCELZiJgIgiAIZSNiIgiCIJSNiIkgCIJQNiImgiAIQtmImAjCDEIp9VGl1M+u9HUIgl9ETARhZnE98NyVvghB8ItUwAvCDEApVYXV2dV1eEhrPStnagjzD7FMBGFmkAJ2Wj/fDnRhDjQShFmBWCaCMENQSt2B2bm3UWsdL3a+IMwkxDIRhJnDduCACIkwGxExEYSZw3bgxSt9EYJQCiImgjBzuApzWJMgzDpETARh5hAE1iqlliil2q70xQiCH0RMBGHm8CfAG4EzwJeu8LUIgi8km0sQBEEoG7FMBEEQhLIRMREEQRDKRsREEARBKBsRE0EQBKFsREwEQRCEshExEQRBEMpGxEQQBEEoGxETQRAEoWz+PzC0L4oTCUe2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbad0f16c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#This is just pasted from above\n",
    "def damp_spring(x,t,k=1, b=0.1):\n",
    "    #This will get x as a vector!\n",
    "    return np.array([x[1], x[0] * -k - b * x[1]])\n",
    "\n",
    "x_init = [1.0, 0] #The position is at 1.0 and has no velocity\n",
    "t_points = np.linspace(0,100, 250)\n",
    "\n",
    "x_points = scipy.integrate.odeint(damp_spring, x_init, t_points)\n",
    "\n",
    "plt.plot(t_points, x_points[:,0]) #Only plot x[0], the actual x-position\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$x(t)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Since we're in python, we can do some programming to explore the parameter space. Let's do a for-loop to see the effect of different dampenings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAETCAYAAADDIPqYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlgVNXd8PHvvbMlmZnsCxDWBLjsm7IK4o5Wi6gVqwKCFltptY9VfNQ+D1WpKA/q01beujyvL/oI1boWRUAEhKLWhX1JuAnIEgLZt5lJZr/vHxPGBDIkk0ySiTmfdkzmzr3nnDm53N895957jqRpGoIgCILQFLmzCyAIgiBELxEkBEEQhJBEkBAEQRBCEkFCEARBCEkECUEQBCEkESQEQRCEkPSdXYBzKYpyKbACSAC8wP+oqvrnzi2VIAhC9xRVLQlFUeKAtcDTqqoOAa4C/kNRlGs7t2SCIAjdU1QFCaAvkAh8CqCqahGwDxjRmYUSBEHorqItSBwB8oA7ARRFyQJGAls6s1CCIAjdlRRtw3IoijIZWAdoQBLwhKqqS1uyraZpmiRJ7Vk8QRCEH6OQB86oChKKovQE9gNzVVXdqChKKrAeeFNV1Reb276szK7JrWgbybJEYqKZqioHfn/01Ee0EvUVPlFn4RH1FZ621ldysiVkkIi2u5suAWpUVd0IoKpqmaIoHwMzgGaDhKZp+Hytz9zv1/D5xA7ZUqK+wifqLDyivsLTHvUVbdckcoBMRVHGQ/Bup6uBvZ1aKkEQhG4qqoKEqqo5wD3Aa4qiqATubDoELOvUggmCIHRT0dbdhKqqa4A1nV0OQRAEIcpaEoIgCEJ0EUFCEARBCEkECUEQBCEkESQEQRCEkESQEARBEEISQUIQBEEISQQJQRAEISQRJARBEISQRJAQBCFqPfTQA7z55uudXYxuLeqeuBYEQTgrL09l1qxb2z0fn8/Hyy+vZMOGj3G53EyYMJHFi39PYmJiq9ZvSXqbN3/KBx+8y5Ej+bhcTrZv/6bdv2driJaEIAhRqaysjLKyMrKzB7Z7XqtXv86OHdt55ZXX+fDD9QAsXbqk1eu3JD2rNZ6bbvoZDzzwu0h/nYgSLQlB6Ka8Pj8VNc4OySs5Pga9Lrxz0pycHKxWK/v27eG1116hqqqS6dOv4JFHfo9eH9lD10cffciCBb8gM7M3AIsW/ZbbbpvFmTOn6dmzV9jrtyS9iRMnA7B7986IfpdIE0FCELohr8/P469+TVl1xwSJ1IQYlt07KaxAkZOTQ11dHYWFp1i1ag2VlZUsWnQP69atZdasW5rc5rnnnmXz5o0h07zzzvnMnTu/0TK73U5xcRGKMjS4LDOzN2azmaNH888LEs2tb7XGh5VetBNBQhCEqJSTk8OECZNYsGAhAHFxZqZMmcaRI3kAbN++FaPRyOTJU4PbPPzwozz88KNh5eNw2AEwmy2NllssVhwOR9jrh5tetBNBQhC6Ib1OZtm9k6K6uyk3N5d77vllo2WVlRUMGqQAMH36FREpW1ycGfjh4H+W3W7DbDaHvX646UU7ESQEoZvS62TSk+I6uxhNcjjsFBQUkJycHFxWUVHOrl07WbDgXgAWLpzHypWvYjLFBNdZsWIZmzZtCJnu3LkLmDfv7kbLrFYrGRk9UNXDwQBUWHgKh8NBdvag89Jobv1w04t2IkgIghB18vLy0Ol0bNy4npEjx1JRUc7SpUu46qprUJQheL1efD5/owABsHjx4yxe/HjY+c2ceRNr1rzBuHEXk5CQwEsvvciECZNDXj9obv2WpOfz+fB6vXi9XgBcLhcARqMRSZLC/g7tRQQJQRCiTl7eYSZNmkRqahozZ84gJiaG66+fyV133QPAsWNHycrKjlh+c+bMx2azsXDhPNxuD+PHT2TJkqXBz1esWEZRURHPP/+XFq3f3OcAn366nmXLngy+v/LKSwB4992PouritqRpWmeXIWJKS22t+jI6nURysoWKCjs+34+nPtqLqK/wiToLT3P1tW7dWmpra5k9+/ZOKF30aev+lZZmDdl0EQ/TCYLQ5eTnqyjKkM4uRrcgupsEQehyHnzwkc4uQrchWhKCIAhCSCJICIIgCCGJICEIgiCEJIKEIAiCEJIIEoIgCEJIIkgIgiAIIYkgIQiCIIQkgoQgCIIQkggSgiAIQkhR98S1oijJwCvAJEACXlVV9anOLZUgCJ3hoYceYOTIsefNJid0nKgLEsAq4BTQF0gH3lMU5W1VVfM6t1iCIHS0vDyVWbNubfd8fD4fL7+8kg0bPsblcjNhwkQWL/49iYmJTa7/17/+ha+++oKSkmJiY2OZMmUq9913P/HxCQA8/fQTbNq0AaPRGNzmvvse4Oabf/guzeU5Z85siovPNFjfj9vt4rXXVqMoQ87LQ5IkFi26P+L1FVVBQlGUXsBPgF6qqmpAMTCtc0slCEJnKCsro6ysjOzsge2e1+rVr7Njx3ZeeeV1EhISeeaZJ1m6dElwaPBz6XQ6lix5iqysgdhsNv74xz+wbNmTPPvsC8F1rrvuBh599D9bnefq1e80Wv+VV/4PO3ZsbzSw4dk82nOU4agKEsAYoARYoCjKXMAPvKyq6kst2ViSJORWXGWRZanRT+HCRH2FT9RZeGRZ4uDBHKxWKwcO7OGBB16hqqqSyy67gkcf/T16vSGi+X388YcsWLCQvn37APCb3/yWW2+dRUnJmSbndli06DfB302mZG69dTZ/+MPv0ekCf19JCrzOvm9rnl6vl/XrP2LevAVN5tGe+1e0BYkkAl1MLlVVRyqKMgr4QlGUI6qqftbcxikp5jbN6JSY2PXmn+1Mor7CF0115vV5Kaur7JC8UmOT0OvCO9zk5OTgdDopLy/ho4/WUlFRwR133MHWrZ/y85//vMltnnjiCdatWxcyzXvvvZd777230TKbzUZRURETJ15EcrIFgOTkoVgsFoqKChg+fHCzZT14cC9DhgwJbm8yGdiyZQv//Oc2kpKSuOKKK/jNb34TnOM63Dw3btyIw+HgjjtuIz6+uTwszZY3HNEWJKoADVgJoKrqfkVR1gHXAc0GifJyR6tbEomJZqqqHPj9YkKY5oj6Cl+01ZnX7+UPX/4X5c6OCRIpMUk8eckj6OWWHXJkWSInJ4cJEyZyxx3zcbvBYklm8uRL2Lt3P9dccwPbtm3FaDQyZcrU4HYPPPAwDzzw8AXTrqiwN3pfXFwMgM8nN/rMYrFQXFx+3vrn+vzzLbzzzjv89a+vBtedOfMWfvGL+0hMTOL48WP88Y9PUlBQyFNPLWtVnqtX/40rr7war1duMo+TJ4+zbNlT/Pu/PxbMIxxnA1VToi1IHAEMgBmoqV+mAd6WbKxpGj5f6zP3+zUxa1gYRH2FL1rqzNcJgcrn05DCmAkzNzeXe+75ZaP6qqioYNAgBZ9PY9q0y4PptoXJFAdATY2NjIwf0rLZbMTGxl0w/a1bN7NixTKeffYFBg4cElx30KDAdQNNg379srj//ge5//5f8thjf8BoNIaVZ2HhKXbu/JaXX17VaHnDPPr3z+Kxxx5j7ty5wTwiJaqChKqqqqIoXwKPA48qitKfwIXsWZ1aMEH4kdHLepZMWkyls7pD8kuKSWhxKwLA4bBTUFBAcnJycFlFRTm7du1kwYJAd9HChfNYufJVTKaY4DorVixj06YNIdOdO3cB8+bd3WiZ1WolI6MHqnqYQYMUIHBgdjgcZGcPCpnWJ598xMqVf2L58hcYNWrMBb+PXN/FcXa66HDyXLv2fQYOHMTw4SPCyiNSoipI1PsZ8KaiKCcAB/CYqqrbO7lMgvCjo5f1pMWldHYxmpSXl4dOp2PjxvWMHDmWiopyli5dwlVXXYOiDMHr9eLz+RsFCIDFix9n8eLHw85v5sybWLPmDcaNu5iEhAReeulFJkyY3ORFa4B3332bVav+hxdeeJGhQ4ef9/nmzZ8yceIUrFYrBQUnWbnyT0ydeikmkymsPD0eD+vXr+MXv/hVs3ksX76cadOmN8ojEqIuSKiqWgRc3dnlEASh8+TlHWbSpEmkpqYxc+YMYmJiuP76mdx11z0AHDt2lKys7IjlN2fOfGw2GwsXzsPt9jB+/ESWLFka/HzFimUUFRUFb0/985+fQ6fT8cADjQ/en322A4B//ON9nn9+OR6Pm6SkZC699DLuvrvxBfPm8gTYvn0rLpeLa6657rwyn5vHjBnXMGfOgojUR0NSpJsmnam01NaqL9Oe9xj/GIn6Cp+os/A0V1/r1q2ltraW2bNv74TSRZ+27l9padaQt4WKsZsEQehy8vPVRg+VCe0n6rqbBEEQmvPgg490dhG6DdGSEARBEEISQUIQBEEISQQJQRAEISQRJARBEISQRJAQBEEQQhJBQhAEQQhJBAlBEAQhJBEkBEEQhJBEkBAEQRBCEkFCEARBCEkECUEQotZDDz3Am2++3tnF6NbE2E2CIEStvDyVWbNubfd8fD4fL7+8kg0bPsblcjNhwkQWL/49iYmJTa7/9NNPsGnThkYzwN133wPcfPOtLU6vuXWay6OjiJaEIAhRqaysjLKyMrKzB7Z7XqtXv86OHdt55ZXX+fDD9QAsXbrkgttcd90NfPbZjuCr4cG7Jem1ZJ0L5dFRRJAQBCEq5eTkYLVa2bdvD7fdNosZM6azbNmTeL0tmvI+LB999CFz5swjM7M3FouFRYt+yzfffMWZM6fbLb1I59leRHeTIHRTmteLp6KiQ/IyJCcj6cM73OTk5FBXV0dh4SlWrVpDZWUlixbdw7p1a5k165Ymt3nuuWfZvHljyDTvvHM+c+fOb7TMbrdTXFyEogwNLsvM7I3ZbObo0fyQU5hu27aV7ds/JyEhkWnTprNgwULi4uJalF5L8wyVR0cSQUIQuiHN6+XYfzyKt6ysQ/LTp6Yy4I/PhhUocnJymDBhEgsWLAQgLs7MlCnTOHIkDwhM7Wk0Gpk8eWpwm4cffpSHH340rLI5HHYAzGZLo+UWixWHw9HkNrfcchv33Xc/iYlJHD9+jGeeeZLly4t58sllLUqvJetcKI+OJLqbBEGISrm5ucyY0Xhu58rKCpKSkgGYPv2KRgGiteLizMAPB+6z7HYbZrO5yW2GDBlKcnIKsiyTlZXN/ff/jm3btuB2u1uUXkvWuVAeHUm0JAShG5L0egb88dmo7W5yOOwUFBSQnJwcXFZRUc6uXTtZsOBeABYunMfKla9iMsUE11mxYhmbNm0Ime7cuQuYN+/uRsusVisZGT1Q1cMMGqQAUFh4CofDQXb2oBaVV5YD59uaprUovdbk2TCPjiSChCB0U5JejzE9vbOL0aS8vDx0Oh0bN65n5MixVFSUs3TpEq666hoUZQherxefz98oQAAsXvw4ixc/HnZ+M2fexJo1bzBu3MUkJCTw0ksvMmHC5JDXIzZv/pSJE6dgtVopKDjJypV/YurUSzGZTC1Or7l1msujo4ggIQhC1MnLO8ykSZNITU1j5swZxMTEcP31M7nrrnsAOHbsKFlZ2RHLb86c+dhsNhYunIfb7WH8+IksWbI0+PmKFcsoKiri+ef/AsA//vE+zz+/HI/HTVJSMpdeehl3331vi9NryTrN5dFRpI5uurSn0lJbq76MTieRnGyhosKOz/fjqY/2IuorfKLOwtNcfa1bt5ba2lpmz769E0oXfdq6f6WlWaVQn4kL14IgdDn5+SqKMqSzi9EtiO4mQRC6nAcffKSzi9BtiJaEIAiCEJIIEoIgCEJIIkgIgiAIIUXtNQlFURKBQ8BnqqrO7+TiCIIgdEvR3JL4M+Dq7EIIgiB0Z1EZJBRFuQEYCKzu7LIIgiB0Z1HX3aQoShLwJ+B6IKwnZSRJQm5F2JNlqdFP4cJEfYVP1Fl4RH2Fpz3rK+qCBIFupr+qqqoqihLWhikpZiSp9ZWUmNj0iI9C00R9hU/UWXhEfYWnPeorqoKEoig/BbKA+a3Zvrzc0eqWRGKimaoqB36/GDKhOaK+wifqLDyivsLT1vpKTraE/CyqggRwG4Eg8X19KyIR0CuKoqiqOrm5jTVNw+drfeZ+vybG1QmDqK/wiToLz4MP3s/IkWPPm01OaFp77F9RFSRUVZ3T8L2iKE8A/cUtsILQPeXlqcyadWu75+Pz+Xj55ZVs2PAxLpebCRMmsnjx70lMTGxy/TlzZlNcfKbB9n7cbhevvbYaRRnC008/waZNGzAajcF17rvvAW6++dYG21w4z7/+9S989dUXlJQUExsby5QpU7nvvvuJj08AOC8PSZJYtOj+iNdXVN7dJAiCUFZWRllZGdnZA9s9r9WrX2fHju288srrfPjhegCWLl1ygfXf4bPPdgRft912B/37ZzUadPC6625otE7DANGSPHU6HUuWPMX69Vt4/fW3KCkpYdmyJxulcTaPrVu/YM+ePdxyy+w218W5ojpIqKr6hGhFCEL3lJOTg9VqZd++Pdx22yxmzJjOsmVP4vV6I57XRx99yJw588jM7I3FYmHRot/yzTdfcebM6Wa39Xq9rF//ETfeeHNE8/zlL3/N4MFD0Ov1JCUlccsts9mzZ1ervl9bRFV3kyAIHcfn82Ov6ZjnVS3xJnS68M5Jc3JyqKuro7DwFKtWraGyspJFi+5h3bq1zJp1S5PbPPfcs2zevDFkmnfeOf+86xt2u53i4iIUZWhwWWZmb8xmM0eP5oecne6sHTu2Ybfbufba6xst37ZtK9u3f05CQiLTpk1nwYKFxMXFtTrPXbu+PW9q07N5JCYmcvXVV3HnnQswmWIvWN5wiSAhCN2Qz+fnrVe/xVbt7JD8rAkx3H7vhLACRU5ODhMmTGLBgoUAxMWZmTJlGkeO5AGwfftWjEYjkydPDW7z8MOP8vDDj4ZVNofDDoDZ3PgOH4vFisPhaHb7tWs/4Morr8FqtQaX3XLLbdx33/0kJiZx/PgxnnnmSZYvL+bJJ5e1Ks9t27bw0Uf/YOXKV5vMo6DgOM8+u5SCgkL+8IenW/jNWyaqu5sEQei+cnNzmTHjukbLKisrSEpKBmD69CsaBYjWiosLPFtw9sB9lt1uw2y+8HMHhYWn2LXrO268sXHLZsiQoSQnpyDLMllZ2dx//+/Ytm0Lbrc77Dy3bt3M8uVPs3z5C42ueZybx2OPPcbWrZuDeUSKaEkIQjek08ncfu+EqO1ucjjsFBQUkJycHFxWUVHOrl07WbAgMM/zwoXzWLnyVUymmOA6K1YsY9OmDSHTnTt3AfPm3d1omdVqJSOjB6p6mEGDAg/wFhaewuFwnNe9c661a99n4MBBDB8+4oLryfUPcJ2dLrqleX7yyUesXPknli9/gVGjxoSVR6SIICEI3ZROJ5OQFNn+60jJy8tDp9OxceN6Ro4cS0VFOUuXLuGqq65BUYbg9Xrx+fyNAgTA4sWPs3jx42HnN3PmTaxZ8wbjxl1MQkICL730IhMmTL7g9QiPx8P69ev4xS9+dd5nmzd/ysSJU7BarRQUnGTlyj8xdeqlmEymFuf57rtvs2rV//DCCy8ydOjwZvNYvnw506ZNb5RHJIggIQhC1MnLO8ykSZNITU1j5swZxMTEcP31M7nrrnsAOHbsKFlZ2RHLb86c+dhsNhYunIfb7WH8+IksWbI0+PmKFcsoKiri+ef/Ely2fftWXC4X11xz3Xnp/eMf7/P888vxeNwkJSVz6aWXcffd94aV55///Bw6nY4HHmgchD77bEeTecyYcQ1z5iyISH00JEW6adKZSkttrfoyOp1EcrKFigq7eBq2BUR9hU/UWXiaq69169ZSW1vL7NlhjQH6o9XW/SstzRpy0Dtx4VoQhC4nP19tdBFXaD+iu0kQhC7nwQcf6ewidBuiJSEIgiCEJIKEIAiCEJIIEoIgCEJIIkgIgiAIIYkgIQiCIIQU9t1NiqIMAvoDsUApsF9V1eZHwRIEQRC6nBYFCUVR+gO/Bu4AegANH7zwKoryBfAK8K6qqv5IF1IQBEHoHM12NymK8gJwgMDc048Cw4AEwEggYFwL/BP4I7BPUZSL2620giAIQodqSUvCCAxSVbWoic9KgK31rycURbkJGATsjFwRBUEQhM7SbJBQVfU3Z39XFKWcQMCoCLHuhxEsmyAIgtDJwr27KYkmAouiKEmKovw9MkUSBEEIeOihB3jzzdc7uxjdWksvXH8FfAtowHBFUSpVVfU0WCUWaHrSWUEQhFbKy1OZNevWds9n8+ZP+eCDdzlyJB+Xy8n27d9ccH2fz8fLL69kw4aPcbncTJgwkcWLf09iYmK7l7WjtbQl8RmgELiraQtgUxRlt6IorymKcj/w78CpdiqjIAjdUFlZGWVlZWRnD2z3vKzWeG666Wc88MDvWrT+6tWvs2PHdl555XU+/HA9AEuXLmnPInaaFrUkVFX9A4CiKHXAxUAmMAYYC/yKQLB5uJ3KKAhCN5STk4PVamXfvj289torVFVVMn36FTzyyO/R6yM7gPXEiZMB2L27ZffcfPTRhyxY8AsyM3sDsGjRb7nttlmcOXP6grPZdUXh1rRFVVUfcAjY1A7lEQShg2h+H15PdYfkpTckIMm6sLbJycmhrq6OwsJTrFq1hsrKShYtuod169Yya1bTvdvPPfcsmzdvDJnmnXfOZ+7c+WGV41x2u53i4iIUZWhwWWZmb8xmM0eP5nfvIFEfIARB6OI0v4/Tuf8Hn7uqQ/LTGRPpNfTXYQWKnJwcJkyYxIIFCwGIizMzZco0jhzJAwLThxqNRiZPnhrc5uGHH+Xhhx+NbOHP4XDYATCbLY2WWyxWHI4f3+ATLXmYbkRLE1MUxVA/bIcgCEKb5ObmMmNG4/mjKysrSEpKBmD69CsaBYiOEhdnBn4IFmfZ7TbMZnOHl6e9taQl8amiKDuAV4HPVVU9bwJVRVF6AHOA3wDLgPyIllIQhIiSZB29hv46arubHA47BQUFJCcnB5dVVJSza9dOFiy4F4CFC+excuWrmEwxwXVWrFjGpk0bQqY7d+4C5s27uxXf4AdWq5WMjB6o6mEGDVIAKCw8hcPhIDv7x3eO3JIgoRAYjuPvgFFRlF1AIeAEkoHhwEBgB3CXqqrb26msgiBEkCTrMJiSm1+xE+Tl5aHT6di4cT0jR46loqKcpUuXcNVV16AoQ/B6vfh8/kYBAmDx4sdZvPjxsPPz+Xx4vV68Xi8ALpcLAKPRiCRJ560/c+ZNrFnzBuPGXUxCQgIvvfQiEyZM/tFdj4CWPXFtB/5DUZSlwHXApcAAII3AKLCvABtUVT3cngUVBKH7yMs7zKRJk0hNTWPmzBnExMRw/fUzueuuewA4duwoWVnZEcvv00/Xs2zZk8H3V155CQDvvvsRPXv2YsWKZRQVFfH8838BYM6c+dhsNhYunIfb7WH8+IksWbI0YuWJJpKmndd7FJKiKIqqqmo7lqdNSkttLf8yDeh0EsnJFioq7Ph8rUqiWxH1FT5RZ+Fprr7WrVtLbW0ts2ff3gmliz5t3b/S0qznN5fqhTssxxeKooxv6gNFUSJy47KiKFcqivKNoiiHFUXJVxTlwUikKwjCj0d+voqiDOnsYnQL4QaJZ4HNiqJc3XChoih3AG1uYdRfAF8LPK6q6hACw5A/pSjK5LamLQjCj8eDDz7C6NFjO7sY3UK4z0k8ryhKMfCBoij3AmXAfwFDgJciUB4fMFdV1S31+R1VFOUIMAr4VwTSFwRBEMIQdheRqqqr68/41wBeYBVwg6qqhW0tjKqqpUBwuHFFUbKBocCXLdlekiTkVszaLctSo5/ChYn6Cp+os/CI+gpPe9ZX2BeugSeBnxEY9G8q8B+qqv450gVTFKU3gaE//q6q6pPNrQ+gaZrW1O1qgiAIwgWFPHCG25I4SKDbZ5qqqv9SFGUK8LGiKBmqqoZ/c3IIiqKMI3BtYqWqqstbul15uaPVLYnERDNVVQ78fnHnSXNEfYVP1Fl4RH2Fp631lZxsCflZuEHiFlVVPzr7RlXVrxRFmQZsVBSlp6qqC8Iu3TnqA8R64Neqqr4fzraapuFrw+hSfr8mbk8Mg6iv8Ik6C4+or/C0R32Fdd7dMEA0WJZDoNtpYlsLoyhKDPAurQgQgiAIQuRF5NkGVVVPKooSiZG2bgL6A08rivJ0g+Vvq6r6RATSFwRBEMIQsZk7VFWtiEAabwFvRaA4giAIQgS04jKvIAiC0F2IICEIgiCEJIKEIAiCEJIIEoIgCEJIIkgIgiAIIYkgIQiCIIQkgoQgCIIQkggSgiAIQkgiSAiCIAghiSAhCIIghCSChCAIghCSCBKCIAhCSBEb4K8rO7j7W0qPfI2mSfjlONKzRjB01MXIrZnBSGgRr8/P6TIHp8sd2Gs9SJKEOVZPz2QzvdPN6ETdtyuP109hmZ3yaidOtw+9TibJaqJ3moW4GHFYaG8uj4+i8lqq7C78fo1Yk570pFiSrCaibXZNsTcAtoId9M+y17+rAO0U+V9vobwklpT+Exg8enLU/eG6Ik3TyD1RyY79Z9h/tIw6V9MzRJmMOkZnp3DZmEyUvomi7iPE5/ezSy3lywNF5J6oxOvzN7le/x5WJgzN4NLRvUTAiKBap4cv9p/hO7WE42ds+JqYQS7BbGRk/b6f1Su+E0p5vrDmuI52paW2Vn2Zf37yNgmmAowmP3FmF+cek6orYzFaRjN4/NXigAXodBLJyRYqKuwtngUrr6CK97cfJf9UdaPlJqOOhDgjANW1blzuxoFjYGYCM6f2Z8SAlMgUvpO0ps4ixef3s33vaTZ+c5Kyamejz/Q6mTiTDo/Pf17QjjHquGxsJj+Z1A9LrKEji9yp9RVpNQ43n357ks/3FOI8Z/+WJNDJEt4mvmN2r3humNKf0QNTm82jrfWVlmYNeWATQYLGFVxVXsHeLzeDs5C0Hk7i4tzB9Ww1MaT1u5b0/qMiVuauKJwdstbp4c1NeXyTUxxc1r+HlUtG9mREVjLpibHBwKtpGqXVTg5+X86OfWc4UWwLbjNpeAZ3XDW4ww9WkdJZB72TxTZWbTjMiaIf6nJEVjKThmUwuE8iKfExwfq313n4/nQ1u/NK+dehYjzeQEsjPs7AndcoXKykddhJ0o8hSGilNiR+AAAgAElEQVSaxteHivnb5jwcTi8ARr3M+KHpjBmYyoCe8SRaTciShMvt43S5g8MnKvnyYBGnyxzBdC4anMYdVw8myWoKmZcIEi0UiSDRsIJLSsr4btM79M6oIyn1hz+aw55G1sU/xxSb1PZCd0Et3SHVk5X833U5lNe4AOibbuHm6dmMzEpu9mCjaRoHj1Ww9otjfH+6Bgg0xRf+dBjD+idH7st0kI4+6GmaxvqvT/DhP4/hr/83Pml4BtdP6kdmWuhJ78+qqXWzZecpPv3uJG5PIFhcrKSx4CdDiTW1fxdUVw8S9joP/++TXPYeKQMCrbKrLu7DVRf3Jr6+5RyKpmnkHK/koy+PBVveMUYd91w/lIuU9Ca3EUGihSIdJM7KO3mKnI3vMHCwn/iEWgC8XhlL6uWk9Z/S7bqgWrJDbttTyJubVDQt0KVx62XZXHlxb+Qw68qvaWzZeYr3tx/F7fUjSxJ3Xj2Iy8f1jsRX6TAdedBzuX2s2pDLt7klAKQnxnLXtQpDWxFcS6vqeGPjYXKOVwKQmWbm/ltGkZ4YG9Eyn6srB4nTZQ7+8t5+SqrqABiZlcJd1yokx8eElY5f0/hy/xne+fxIsCVyw5R+zJqahSw3/nckgkQLtVeQgEC/7pv/+IC+VafoO7wWozHwR/P6M+gz/GcYYrp2n3k4LlRffk3jg+3fs/7rEwD0SjXzq5nD6Z3e/NnrhRRX1LLygwMU1jfDr7yoN7dfNSjsoNNZOuqgV1Pr5r//vi/YVTdxWAbzrxuCyaBrdZqapvHZzlP8fWs+mgbmGD2/vXU0AzMTIlXs83TVIHHwWDkv/eMgdS4fep3EnVcP5tLRvdp0Illpc/HXDw9wtL5FPXZQKr+6cQQG/Q93ALZnkNA98cQT4Zc6StXWup9ozXayLBEba6Suzk2omClLEmOGDuOIFs+RHTkYTYlYrE5kyUFNyS40yUiMJbNbtCpC1Zdf03h9w2G27DoFwLD+Sfxu9hhSEsI7g2qKJdbA5OE9OFVqp7iyjmNnaqiocTF6YGqXqPOW7GNtVWlzseKtPZwqdSABt16ezW1XDESva9vtxJIkkZ2ZwMDMBPYfLcPh9PJtbgmDeieQmtA+LYqOqK9I25tfxsoPDuD2+ImPM/Dg7DGMG9z26zixJj2Th/fAVuvmRJGNoopajp2u5qLB6cG/bVvry2w2PRnqMxEkCK+Cs3umofUczPZvv0R3Oh1LmgeDwYvbfhRH1QniEgYg69p+UIxmTdWXX9P4340qO/afAeCSET341Y0jMBlbfwZ7LoNeZsLQDKrsbk4U2zhZYqes2smYLhAo2vugV1Hj5L/+tofiyjp0ssS9M4dz2ZjInrSkJ8UydlAqu/NKcTi9fJdbwoCe8aQnRT5QdLUgsfNwCS+tPYjPr5GRFMujc8a1ufXckE6WGJ2dgk6WOHyyitIqJ+rJKi4eko5BL4sg0VIdESQAMhItpPQaxSfff4Vlt4wnIRGrtRbNW01N6R4MpngMMelRf+BqrXPrS9M01mzKY/ve0wBcOroX838ypF0eiJMkiVEDU7DVejheZKOgxE6lzRX1gaI9D3o1te5ggNDrJH5908iQFzjbyhpnZOygVPbkl2Kv87JLLWFI36Sw+9ub05WCxL4jZfz1w4P4NY2eKXE8cse4iNcHBPZ9pW8SsSY9B49VUGFzcexMDROGpqPXyyJItERHBQmAlPhYeqeP4N2qg/TYW0SRczBJqTb0Og911Yfx1BUTYxmArLvwnQxd0bn19fGXx9n47UkApo7syV3XDWnXawWSJDEqOwV7nYdjZ2ycLLbj12Bov+i926y9Dnp1Li8v/H0vp0od6GSJ+28Z1aL76tvCHGvgosHp7M4rqQ8UpYwamEK8OXL7elcJEvmnqvjLe/vx+jUyU808csc4Eiyhb1WNhOzMBCyxBg58X05ZtZPCUgfjh6ZjjjOJINGcjgwSACkJMWQmDOIt+0kGHd9PwakRGBM0zGYnXlc5jvJ96E0pGGLa9x9tR2tYX//cd5q3txwBYPyQdH5xw7Dz7rxoSNM0yp2V5FcdZX/ZIfaWHuS74j3sKdnP3tIDHCo/zImaAoprS3B6XcToTRibCLSSJDEiK4XiyjoKSx3kFVRhiTVEzVOq52qPg57P7+fF9/eTf6oaCVj40+FcpKRFJvFmxMXoGZmVwre5JTicXnbnlzJeSY/YE9pdIUicKrXz/Nt7cXp8pMTH8Mgd40hs5wBxVlaveGQJDp+soqiilkqbmymjeuF0eiIeJMQz9200fEAKdzpu4g3pA2YUfEHh56MpGj6IYcpR9NRSduzvmJNHk9R7xo/uWsWhYxW8sUEFQOmTGDJA1HpqOVCWy6HywxyuzMfhqQ0rn5SYZAYnZaMkDWR4yhDiDIE+cFmSuOf6odhr3Rw6XsnfNueRnhTLyKzucafZe9uOcqj+1tQ7rxnMxGEZHZp/zxQzD84ezX+9tYdqu5sX39/PY3Muiuh1qGhlr/Pwl/f2U+vyYok18NDPx1zwYbf2cMOU/jicXjZ9V8A/951mYN+jTB/VI+L5iJYEbT9r6ZNuRa5L5xOtkl7+gyTn+TlQM4aEBAexsS48dcU4Kg5ijM1Ab4reLpGWkmWJmjovT636BrfHT69UMw/9fAwxxh/OOTRNI7cij7Xfb+Rvh99jT+kBzjiK8fg9gTQkmYy4NDItPRmQ0Jc+1kwyLT1Ji00hVh8IAk5f4CG8Om8dp+yn2Vt6kM8LdnDCdgo0jeSYZIx6A2MHpbH/aDk1Djf7jpQzbnAq1mYeWOpokT4z/jqniHe2HgXgqot7M/OSAW1PtBWSrCb6plv5JreYaoebM+W1XDyk7dfjorkl4fdrrPzwACeKbOh1Eg/9fAx9060dXg5Jkhg2IJkz5bWcLnNgr/Nw6eheEW9JiOckiMw92Zqm8X8/yWGnbSsTKlRG5cRwoMdlZA6qZPCg4+jkQLqWtAkk9roSWe6aw0tAoJvjmTV7+L6wGnOMnv+cPz74cJXP72Nn8V42n9zOaUdRcBu9rGdwYjbDUhQGJg6ghzkDg3zhhmytp5ZC+xmOVh9HrTjC99XH8Wo/jH1j1BmZ0GMc0zOnYPQm8NQbO7HXechIjuM/5l2EOSZ66jiS9/2fLLax7M1duL1+lD6JPPTzMW2+zbWtPv32JH/fGuh2vGFKP26+NLtN6UXzcxLvbjvChq8D1+AWXDeEaaN7dWp5PF4/Xx48w7ihPUiK04uH6S6kM4MEgMfr45k1uymM/ZIRNUeZ/p2bnIzpeNITGDNSJSE+MNKs3pRCcp/riLFmtTqvzqJpGq+tz+WrA0VIwL/NHs3IrJTAMBrluXx4ZD3FtSXB9QfE92Vyz/GMyxgVbCG0Vq2njv1lh9hVvI/Dlfn4tR9GMR2UmMXg2LF8+LEDnx+GD0jm324dFTVDjkdqH7PXeXjq9e8oq3aSZDXxh/njI3rBuLU0TWPVhsN8UX8L9L0/Hcak4a3v+ojWIPFtbjEvrz0EwOVjM5k7Q+nkEgW058N04ppEBBn0Ou6/eRRPvlHHoR4+nFNPct0Xn3HSNZYvHWMYlH2S7KwCvK5ySo6sJjZhCEmZ16A3JXZ20Vts6+5CvjoQaCHcPD2LkVkpFNgK+SB/HXlVR4PrjU4dzlX9LiMroV+j7TVNw1dTg6ekBE9FGT6HA7/djs/hQPO40QL31AIgGYzIMTHIMTHo4szoExMZk5TKxQN+Ru0QiW+Kd7Oj8F/1F8K/J7/qe1InJFOWn8mh437e/fwoP79yUMdVTjvz+f28vPYgZdVO9DqZ39w8skUBQtM0qt01FDtKKXdWUuO2YXPbsLntuHwuvH4fXs2Lz+/HIOsx6AwYZQMmvYl4ozX4SolNIj02LXhNqCFJkpg3Q6Gkopa8U9Ws2nCYHilx9O8RnTcStMapEjv/b30uAAN7J3D7VT+efetCoi5IKIoyHngRSAU8wDOqqv5v55aq5ZKsJu6/aTTL33JzNNvLB1fK3Lh9DwnOUg75plFcksKIYfkkJtipqz5MXU0+8emTsaZPRtfGM+32ln+qire35AMwcXgPrp7Ui/fyP2JbwZdoBA7sAxMHcPPAG+gX3wfN58NVUIDz2Pc4j3+P8/hx3MXFaC7nhbJpGZ2OgcnJDM3oQU1iXw7rK8jVV1AZX4Yxqxx/73y2nMiiz6FYLhnetcZ5CuWD7d8Hx1CaN0NhQM/zD8B+zU9JbRnHa05yvKaAk7ZTFDtKgtd3IsFqsJAel0ZGXCq9LD3pY82kt6UXMXoTi24eydLXd1Je4+T/fHCA/5w/vtkB7boCe52HFz/Yj9vjJ9Fi5NezRnR6F19HiaruJkVRTMBR4Peqqr6hKMpw4EtgmqqqB5rbvrO7mxrasf80qzYexKjsJEUrY9bnVRhdcRzKmEZ1bDq9exUzVDmO0RgYilySjVjTJkRtsKiyu3hy1XdUO9z0SI5j4dyevL73LcrqKgBIj01l1sCfMFRLozbnEI5DB6k7nIu/ri5kmnJsLDqLFSxW/HHxePUx+CQdfnT4Ab/Hi8/txe92o9XV4XfYkHxeJPxImoas+dFpHnR+D3q/F1nzIgFOo0R5gp7yBD1lVhP9sscxbdJPsCR23q3Ibd3HGnZzXDmuN3deMzj4mdPr4nBFHvvLcjhYnnvBu8esBgvxph9aBzF6EzpJh17Wo5NkvH4fbr8bj89DrddJjdsWeLlqcNffdNAUCYmMuDT6WHtjJZVN/6zBY7MypHcKD/18TNjdftHU3eT3a/zp3X0cPFaBXifx73eMI7sF41ZpmobT56TGZcPmceD0OnH6XMGfLp8bv+Zv9ALq/xaBv4leDvw0ykbi9DHEGWKJ1ccFfzfpTMiS3H0G+FMU5SfAq6qq9m6wbA1wQlXVx5vbPpqCBMD/fqqybf9xTEO+xSJXceO2alIrfZxIGsmxlLHodD6yswrIGnAGWQoMGCjJRiwpY7GkXhQ1z1d4fX7+6297OFJYjSlG46LLS9lTvhsI3KX007hxjCnSU7trF+7ThcHtNMCli8NpTcfToz/ehHRcRitODNS5JeqcXtxOb4M619Dp/Mhyw5eGLNX/lJueSQ1AkjTQQNa8yJoPyedD5/Mh+X3IXi86nxdJ0jDGxWKymDHFxxObkkRMWhomSxxGkxGjyYDBaMAUo0enkyP6BHdb9rGCEjtPv7kTt8fP4N4JPHz7WOxeGwfKctlfdoi8iiONLugDmHRG+pt7M8CYQU9DEik6C0mSGaMP/C4XmssVCL4+L/h8aH4/ms/X+HdNQ5JlkGVAwqV5sHtrsXtrqfE5sLkdlHtrsGl1eHUSPlnCp5Pw6sCnk/DIEl63mRRzT6YMGUbvpL70SshEZzA2W7fRFCTe336UT/4VGLBy/nVDuLT+QrXH76W8roKyunLK6iooc5ZT6ayi2mULBlfPBQJrJMiSjMVgJt5o4cZh1zAyYcSPPkj8DvipqqqXN1j2B2Ccqqo3Nrd9WZlda811SlmWSEw0U1XlwN/ElIKt5fH6eWb1Lo4WlxIz7FuMOhs//cJGnzMuqk2pHO5zJXZiMRrcDB50hr59CpHwBrePsfbHmnYRcQmDO/XJ7f/dqLJl1ykkczXpow9T460krs7P5DMmRp3w4j9TjFfSY4tJpcaUii0uHaclFY9Rh97oJsbkJsbkIibGjcnkxmDwotd50et9GPSBn2df0SBwWUQKvJBAkwM/kYCzwSNw4ESSA79LEhJy/XsJSZJ/+AwJSZLQ6WR8Pj9afUpn/9sw3x/yD/zu82sUldfi8/nQy37iYsDjdeP3+ZA0kDWQNAlJAx1S4H9aIJHG/7QlkKjvFGyQryShNVynwQ/tnPcX0jCV0Os08UZq/Emjt5LUxCrS2f9fIPEmkw6hcbmlc5b5/Ro+vz9YHEkiULf1HazS+Uk06ULFkBp+qp33S4iSnvuhRFG1zA13PtyqY1hysqXLXLg2A+f2T9TVL29WSoq5TWd/iYktyiYs/3nPJP7the1U5Y5HN/JbPpwucfk+JyNzyxh/5F0Ke0/ke90QDub0Q83vxcgRFfTqeRrNZ8NpO47TdhxJ1hOfMpjE9OEkpA7BYOq4i4FbvjvJll0F6Hsew5SZT8oJJ5ed9JDuNONKiONMHwv+oVnIZhmTyU1GjIt+ptOYTMfPmwa2qwgcDDR++GfZiuDVcPN6fh/BUNNkvk0s08vQr2Meoha6OL0prl2OYdHWkngQmHlOS+JJYLSqqrOa2z7aWhJnHT5ZyfLVe9CMdiwjd+KV6xhzws/0b6rA68UZm8ix4Tdyuqr+zEnSUIY4yRpQjOQr4Nyjjd6YiMnShxhLbwyx6RhMKegMljYFSE3z4/fW4fM68Hns+Dx2SstL+TYnjxRLGSloxBokDDEasi7cOpLRGS3oDVZ0BiuyPhZZZwq8ZFPwd0lnQpaNSLIeSdKBrEOSAr9Lsi6wrPEp6Dm/nv1FA82PpvnR/D7e2pzHt7lFyDo3Ey9x833tIZxeJ3FuP4kOH/3d8fTyWjE6PLirbfh8fnw6PX6dHr+sC/zU6dBkXf17Hcj1p5UyDX6XzjY2Ai0LCSRZqw+W59TZuSfp9d1lUvD8Xas/Oz/bhgk0L852qyEF7iiSJAlJDvxEDrRgJFlCkuVAV1Hw/dl161s90g8FkKRACyR4hi4RfB9oITXt/L1ACjaDNE0Dv4bP5+d0qT1w55TOT7xFh9fvweP14PV7wF///eqD6tnfJQ1kZGSk4EvSzp65n1MmreFPrcGvgbrSGi5DCrYCtOAn538X7bz9q77OztY5MrIkB8YoC9andM76Z7c/Z589W78NFgVbTcG/S1P7dxO0s6tL9B92Ccn9h0S8JRFtQeIaYJWqqpkNlr0DHFZVdUlz20fbNYmGPvuugLe25CPF2rCO3IkHF/2r9Nz4lR2qApOJ1I6YxtH4kZSU/HDhMSVNx7ARblISi/HWHUfTvE2mL8kG9MZEZH0csi42cCCWjdQfBZCQ8GteNL+n/uXG73Pj99bi8zrwe2tpUbu5AU0Dr88Ekhm9KZ5YcyKm2AR0Bis6gwWdIb4+KMR16gitHq+PZ1bv5niRjViTnkfnjCS/7gDbTn1FubMiuJ7VYOGi9FGMtyhkuIx4KyrwlJfhLS/DW12Nz2bDV1OD127D73DQ5keBJQk5Lg4tNgaXUaJG56VCduI0SLiMEk6jhNMkoxnjGdJ/KAN7DSMzPRu92Yykj7ZOgNCOF9Ww7M3deH1+hvZL4ne3jUYny/j8Pk47ijlRc5ITNac4ZS+kuLYUl8/dfKLtwCDrSYtNJTU2hZMn/ZQUy0huM4uum8iofpnIUvTezdSdLlwbgCPAE6qqrlIUZTSwHZioqqra3PbRHCQ0TeN/Ps7h65xiZHM18SP34PI7SfAamJ8bDwcC919LcXF4rpxNviOB0yerG6XRs7eZrGw/qal29HIp7tpC/F5HU9m1mdutx+ky4XQZcTmNuJwGdJoec1I6yZm9yVb6gz4OTYvefzgNlVc7efL177DXechMNfP7eRdhNMgcKj/M9lNfkVuR12j9jLg0RqQOZViyQnbigPOeDtf8fvxOJ5rbhd/lrv/pQvOfc3Fd05D0evQxRhJTE7HVean21nG8rohc+zFyKvKodFU12kRGxleTjKcijT4x2Tx+29Quf7vllwfO8NongX382gl9mX3FwCbXO/tMR6mzDJtWzYny01TWVVPtrqHKVUO1qxqPv+kTpQuJ0ZmI1ccSb7KSZEogsf6VZEogMSaRlJgkEkzxyJLMB/88yrqvAheq77pWYfqYzGZS73zdJkgAKIoyBvgrkAY4CQSM91uybTQHCQjMPfz0m7s4VWrHnFRH7NCdOLwOjJKeX7rGoP/k8+AtozFZ2UhXzuJYtYkjOSW4nI3/Yeh0EqkZVlJS9SQmeTGb69DratHr3GiaE/xO0Lz4/X78fh9+nx+fT8LtlnA6JZx1Gh6PjNttwO024HIbcbkNuF1GPC4dsa5q4p1lYLIz5JKLGDj9EnRGQ4fWV6Tlnqjkubf3oGmBEWt/dePwYAunpLaM74r38G3RbsrqyhttZ5QNZCcOoF98H/paM8mISyMlJhmD7sLDfnj9XsrqyimuLaXUWUaxq4TDJUepcFaet67FYGZossLA+EF8vMFBaYWXBIuRP8wf32Eji7a3NZ/lBWctvHfmMCYNC/1Edqh9TNM0PH4PLp+7/vXDraRnu39kKdBRZdQZA7eM6mLQyS0bdHDn4RL++o+DAEwf04u7rh3Shm/ccbpVkGiLaA8SACWVtTz1+k5qXV4yM8Gf9S+qXNXoJR1zMq8jc+sB7Du/C64fN3wECVdfR5kxnZPfV3Dy+wrsNZF7MEpCw+KuwuIsxeqqwOoqR0clan8jFVnjWXTTXec1s7tqkADY+M1J3vk8MMbQ7MsHcu3Evo0+1zSN4zUn2VNygNyKvEbjTzUkIWExmDEb4ogzxCJLMhISLp+bWm8dtZ5a6rzOBr3ejcmSTP/4PgxLHsKwlMH0sWbi98N/v7OP3BOV6HUSj9wxrl3nke5oXp+f597eS15BFUa9zONzL6JvRtMD43XGPnaiyMYza3bh9vjJ7hXPI3eMazSPdDQTQaKFukKQANh/tJw/v7sPDbh4pIWS5M8prT97vabf5Vzp7E35++/hOnkiuI2xVyYJUy/FPG4cTr2FkjM2SotsVFfWUVNZh8PhxlV3/ljykgQxsQZizUZiTDpi/bWYaorRFx4hproIs7samUAXSUG6gYMDY8lPTSbDOZXHbrmiyW6OrhwkNE3jlY8O8W1uCZIED982hqH9k0OuX+WqJrc8j2M1JzhZc4pCR1GjMaNaQpZk0mJT6JeUSWZcL/pb+9HXmnnePBlvblL5fHfgOZO7fzKUqaN6hv8Fo1yNw82Tr39Hpc1FakIMS+aPxxJ7fouso/exKruLpW/spNLmIslq4j/vurhLteBEkGihrhIkAD768hj/2HEMgFuv6sNheTNHqgLvh6UozFVuRTp4mIqNG3AdP9ZoW1OfvsQMHERsVhaGjJ4YkpPRWa1okoTX4w/0j9c68JaX4ykpwlN4CufRI7hOFcA5feau1Hh29dY43M+IzazDW9yH+OrRLJk3KeRw2105SECg2++Pb+6ksNSBJdbAkvkXk5rQsqfcfX4f5c5KSuvKqXHbcHgcgRaDpuHX/Bh1BuIMccTpYzEb4kiNTSE1JhmjQX/BOvt89yne3BS4LnKhPvsfg2NnanhmdeBC9vD+STw4e8x585B05D7m9vhY/rc9HDtTg9Eg89idF9GvR8cP/d0WIki0UFcKEn5NY+X7B9h7pAydLPGbnw3noGs7X50JdDVZjRbmDb2NocmDcX5/lOptn2Pfu/uCw1wAgabDhf6mkkTswEE4BmWyMfYk+caawF023jjqjg7FUNuD38+96IKTuHf1IAFQXFHLU2/spM7lpWdKHI/Pbd+hxS9UZznHK3jh7/vwaxqjslN44JZRF5zd78fgi/1ngoPlXTuxL7MvbxwUO2of0zSNVz/O4ZucYgB+fdOIdpsfvD21Z5AQkw7ROROcSJLEyKwU9uSXUlPrYW9+ObePn0pWWg8OV+ZT563ju+I9FNeWMnjAWFLHTyHpmmuJHTgIfWISkizjs9vBd+EHvXQJCcT0H4B1/ASSZ1yHd+bVrO1RyifaYSr0biRJJqFOofLACDSnlV/dOIIhzcwVHc0TwrSUJdZAvx4Wvs0toabWw9HCGiYOy0DXTgfnUHV2osjGf7+zD4/PT8+UOB68dQxGw49/Zre+GVZstW6OnbFxpLCauBg92b1+uP7SEfuYpmm8+/lRtu89DcDNl2Zx2djov5OpKW2tLzHpUDM688y4rLqOp9/cRbXdTYLZyONzL8JvtPPGobc5YSsAAuPwXNZ7KtN7X0KC6YdmsKZp+Ow2vBUV+Gtr0XxeNL8fXawZ2RwYWlsXF4emaeRXfc9nJ7eRU/7DncT94/tiLBnFvv2BO6duv3IQV4/v02yZfwwtibMantFerKTxyxuHt8scFE3VWXFlLc+8uYuaWg8JZiOPzb0oOHlTd+D1+fnze/s5dCzwrMrCnw5jcv0cFB2xj63/+gTvbQsMbz9lRA/uuX5opz7P0xaiJdFCXaklcVZcjIFh/ZL4OqcYh9PLnrwypg3vx5UDJpNgsvJ99QmcPhdHq4+x/dSXlNSVISORFJOEXtYhm0zoExMxpKVhTM/AmNEDQ0oKssXMaVcZXxT+i7+p77O1YEfw4nhabAq3KzdTd3wwX+8NPItx3cS+/LSFU2D+GFoSZ/XNsCIRmFD+dHktpZV1jB2UFvGDxbl1VlZdx3Nv7aXS7ibWpOPhn4+lV2rkh1SIZrIsMW5wKjnHK6myu9h3pIzeaRZ6ppjbfR/7fE9hcNj7MQNTWfjTYVEzQVVriJZEC3XFlsRZuccr+NN7+/F4/aQmxPDvd4wjJSEGu9vBloJ/sqPwX9R5f5iHwSAb6GXuQU9LBgnGeAyyHr/mp8pVQ4WzkhO2gkbrA/SxZnJ138sYnTqct7YcDd5Jc8mIHiy4fmhgiIEWiIb6iiRN03h7yxE+2xlouV0yogcLfjI0otcFGtbZmTIHK97aQ3mNC71O5qHbRqP07fpzn7eWrdbNs2t2c6a8Fp0sBVoUI3q02z626duTvF0/1ergPon8bvboLt/FJy5ct1BXDhIAh45V8Of39uP1+UmJj+HfZo8ms/7s0ul18q8zO9ldsp9j1SdC3n9/rnijlYsyRjM+Yyx9rb3x+TVWb1L5577ANJOThwea2eEcEKOlviJJ0zT+9lk+W3YHHvYaNziNe386LGIHj7N1dkAtYsVbe6m0uTDqZe7/2SiGX+AW3O6iyu7iubf3crrMgSTBPdcP5cbLB0d0H9M0jXVfHefD+rsKB/VO4N9uHU2sqesMcRKKCBIt1NWDBMDBY+W8+P4BPF4/sSN4LycAAAs4SURBVCY9v7lpxHn38Ve7bORXHqHQUcQZRzH/v717D46yOuM4/k02ISEkIURQQI3QBB4RjSD1gtYR7473qkWdqdPa2k6rtXWm7R9OO20dZ6ptdWxtsbVOa9VqnfGC98uM1mm9FFCrolEeFLmIQuRigIQEspf+8W5CGnnjbszyvsn+Pv9gDgQfHw/89px995zt3Z0ks0cV1FXUMrZiLPtVT6Kxbir7VO3aOtne1c0tD73Ve7vZYF8xx6lfQ6l/UDROruXKC5qH5Ga1RKKE91s7+PVdL9O5I0VFeYKrvtJc1CuI/rZu38mN977OBx8Hd8Gff3wTZx7VMCTbTTu7U9z5tPPSW8GHI2dOGcf3zmumYtTwXkH0UEjkaCSEBMB7H27h5vuX0t7ZTaK0hPOPa+SUI/bPeTtod9a0buPWR1pYtyk4PPC0Ixq4YF7joLZU4tavoZTJZHhy8ZreNzTH1VTwrTMP+swnvgaSTmd4askaHvjXCjIZqK0q58rzm3O63azYtHd2s+DBN/EPgvOsmhv34tLTZzA2h7u8w2xo6+SPD73FqvXbAJhjwSqxvGxkBAQoJHI2UkICguM7fnvfUtZvDv5SP7Chjm+cMSPnD331SKXTPLFoDY+8sJJUOkOitIRLTrXe27UGI479GmqLWtZz+5PL6E6mKQFOO6qBs46eQuWo/LYmPtzYwV8ff4eV64KTfhv2qebK85rZa2xlAaoeGZKpNPc++y7/zL5nNqayjItPmsbcmRPzeqAglU7zzCtrWfj8++zsDj5E+uVjp3LG0VM+1wuuOFJI5GgkhQTA9q4k9zyzvHeJXJYoYd6sfTlj7gGM/YwjA1LpNIvfbuXRl1bTmg2a8WMruezMg5i+f93nqiuu/RpqH25o59ZHWli7IThpt3bMKM4+ZgrHNk/+zDN9Wj/ZzpOL1vDim+tIZc/3P+nwBubP+8KwP9F1T0gkSnjl3U3c9tCbdO0MPgvUuG8tpx95AIdOGz/gX/LJVDD3n1qyhg97/t9VlfP102cwqykeVwIPNYVEjkZaSPR4ednH3PW0094Z3Jdblihl9rTxHNq0F1Mn1VJfU0lpaQntnd18tLGDd1Z/wuK3W9m0ddfTTcfNmsz845uG5E26uPdrKHUnUyx8fiXPvPIByex/a/XocubYBA6aUs/k8WOoGV1OKp1h87YuVq3bxhsrNtKycnPvXnp9bQWXnj6DeYcfUBQ9Gwo9c+zdVRv52xPLWLpi18m8E+oqmTN9bxr3rWVC3WgqyhN07kzSurmT5Wvb+K9vYEvHrjspvnTIJOaf0LTbM6JGCoVEjkZqSAB07kjyzKtreWrxGjp35H6e/uxp4zn7mKlDehbNcOjXUNvY1snC51eyqGV9zlcz1VWP4tQjGjhu1mTGjC4vup59Hn3nWDKZpmXVZp5ctIZ3Vn/6mPXdKSF4Qu3UIxtG1Em6YRQSORrJIdFje1eSJctaeXXZxyxfu4Xu5KdPJJ20VxWHNo7n6EMmst+E8DOYBms49WuobdrSxUst61m6YiOr12/rXV30qB5dju1fxxybwBzbu3dbqph7Nhhh/Vq9fhtL3mnl7VWfsHZDe+9WHkB5WSkNe1czc2o9cw+eyD7jqqIoPRIKiRwVQ0j0lUqn2dDWxdaOnaTSGaoqyphQN5qqysI+9z1c+zXUkqk0n2zbwfauJCUlUFddQU1V+W7fXFXP8pNLv1LpNFs7uulOpqgYVUbN6PIRfzBimEKGxPD/FEkRS5SWMrG+ion1xfOKKU7KEqVMKKKzluImUVrKuJrhc+fDcKXHLEREJJRCQkREQikkREQklEJCRERCKSRERCSUQkJEREIpJEREJJRCQkREQikkREQklEJCRERCKSRERCSUQkJEREIpJEREJFSsToE1s1LgGuACIAFsAL7v7q9GWpiISJGK20riCuAcYK67TwceBu6JtiQRkeIVt5BYBHzN3duyXz8NTDezyghrEhEpWrHabnL3l/sNnQO87O5duXx/SUkJpYOIvZ7brIr1Vqt8qV/5U8/yo37lp5D92uPXl5rZRcAfdvNTW9y9sc+vuxBYAJzg7ktz+b0zmUxmd1dHiojIgIbXHddmdjVwOXCWu7+e6/dt3NieGexKoq5uDG1tHaTT8etH3Khf+VPP8qN+5efz9qu+vnr43HFtZtcCZwJHuvtH+XxvJpMhlRr8vzudzuiS+jyoX/lTz/KjfuWnEP2KVUiY2cnAJcAcd98UdT0iIsUuViEB/AioBV40s77jF7r7G9GUJCJSvGIVEu5+atQ1iIjILnH7nISIiMSIQkJEREIpJEREJJRCQkREQikkREQklEJCRERCKSRERCSUQkJEREIpJEREJJRCQkREQikkREQklEJCRERCKSRERCSUQkJERELF8vpSERGJB60kREQklEJCRERCKSRERCSUQkJEREIpJEREJJRCQkREQikkREQklEJCRERClUVdQNTM7HDg98B4oBu4zt3vjLaq+DKzDmAdkOwzfJa7vxtRSbFkZt8GbgJ+7u43ZMfGA38BDgbSwCPAj909HVmhMRHSr8eB2cDWPr/0Wne/O4ISY8PMTgR+CYwFEsAt7n5ToeZXUYeEmVUAC4GfuPsdZjYTeNHMXnP3NyMuL3bMrByoAua6+4ao64krM1sATACW9fupPwGbgCagBvg38B3glj1aYMwM0K864Afuft+eryqezGwi8DBwjrs/a2aNwOtmtgj4IQWYX8W+3XQigLvfkf2xBXgcuDjKomKsLvtjW6RVxN8/3H0+sK1nwMxqgHMJVqoZd98KLAC+GlGNcfKpfmXVobnWXwq4xN2fBXD3FcB7wBEUaH4V9UoCOBDov02yHDgsglqGg3EE20x3mFkzsAO41d3/HG1Z8eLuL+xmeBpQAqzoM7YcmLlHioqxkH5BMN+uMLPrCVawjxFsR23fY8XFTHYFv7Dn6+xKYgbwGgWaX8W+khgDdPYb68yOy6ftAP4O/M7dDwYuA643s3OjLWtYGAPs6Lc/rLk2sIeB+4DDgXnAsQR78QKY2X7Ao8B1QIYCza9iD4l2YHS/sTHZcenH3Ve7+6Xuvjj79WvAPQTLXBlYO1BhZn3/zGmuDcDdv+vud7t72t1bgRvQXAPAzA4D/gPc4e7XUMD5Vewh0QJM7zc2A1gaQS2xZ2Zjzax/v0oJngqTgS0n2E9u6jOmuRbCzEaZ2ax+w5pr9AbEE8BV7v6r7HDB5lexh8RzQNLMLgUws0OBU4CifsRuAAcDi83MAMysCZhPnz1S2T137wDuB642sxIzqwMuB26PtrLYKgeeM7OLAMysGrgSeDDSqiJmZpUEW3BXuPsDPeOFnF9Ff+lQ9tXKLQSP4HUBv+jbfPl/ZvZNgkftEsBO4GZ3vy3aquLDzBIEK1SABoLl/maCIP0NcBswi+BV370Eb8QW7R/Cz+jXY8CNBG9gZwiePPypu/d/H7FomNnFBO8L9n/g5l7gZgowv4o+JEREJFyxbzeJiMgAFBIiIhJKISEiIqEUEiIiEkohISIioRQSIiISSiEhIiKhFBIie4CZXZM9819kWFFIiOwZRwJLoi5CJF/6xLVIAZlZGdkTOvsMt7n7uIhKEsmLVhIihZUGvpj955OBSQQXEIkMC1pJiBSYmZ1CcFhdjbvviLoekXxoJSFSeLOBtxUQMhwpJEQKbzbBHcQiw45CQqTwmoHXoy5CZDAUEiKFVwo0mdm+ZlYfdTEi+VBIiBTez4CzgQ+AWyOuRSQverpJRERCaSUhIiKhFBIiIhJKISEiIqEUEiIiEkohISIioRQSIiISSiEhIiKhFBIiIhJKISEiIqH+B+xH/ZeTHViwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbad0cbd630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_init = [1.0, 0] #The position is at 1.0 and has no velocity\n",
    "t_points = np.linspace(0,20, 250)\n",
    "\n",
    "for bi in np.linspace(0.001, 1.0, 5):\n",
    "    x_points = scipy.integrate.odeint(lambda x,t: damp_spring(x,t,k=1.0, b=bi), x_init, t_points)\n",
    "    plt.plot(t_points, x_points[:,0], label='$b_i = {}$'.format(bi)) #Only plot x[0], the actual x-position\n",
    "\n",
    "plt.xlabel('$t$')\n",
    "plt.ylabel('$x(t)$')\n",
    "plt.ylim(-3, 8) #this is to leave enough room for a legend\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Integration Method\n",
    "===\n",
    "\n",
    "We're using the [lsoda](http://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html) solver when we call `odeint`. You can choose between many more by using the more complex `ode` interface in `scipy`. We won't cover this in lecture, but this allows you to explore other types of solution methods including the famous runge-kutta method and methods for complex numbers"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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": 1
}