{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sympy ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will have a quick look at the Symbolic Python or sympy library, which is described in full [here](https://www.sympy.org/en/index.html). This library provides Python tools for solving equations, differentiating and integrating functions, expanding functions as series and other symbolic manipulations.\n",
    "\n",
    "## Sympy introduction - defining equations\n",
    "\n",
    "In order to use sympy, we have to import it in the same way we import numpy. Let's do that and then look at a first example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The quadratic equation is Eq(a*x**2 + b*x + c, 0)\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "#\n",
    "x, a, b, c = sp.symbols('x a b c')\n",
    "quadEq =  sp.Eq(a*x**2 + b*x + c, 0)\n",
    "print(\"The quadratic equation is\",quadEq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see above, the first thing we have done is defined the symbols we want to use. The statement:\n",
    "```Python\n",
    "x, a, b, c = sp.symbols('x a b c')\n",
    "```\n",
    "tells sympy that `x`, `a`, `b` and `c` are to be treated as algebraic symbols. (These don't need to be single letters, for example, we could define `theta = sp.symbols('theta')`.)\n",
    "\n",
    "Our next step was to define the quadratic equation $ax^2 + bx + c = 0$. To do this, we can't write `a*x**2 + b*x + c = 0`, because Python uses the `=` character to assign values to variables (e.g. `y = 18.0` or `y = np.cos(1.3)`). We also can't define the quadratic equation by writing `a*x**2 + b*x + c == 0`, because `==` tests whether what's to its left is identical to whatever's to its right and returns the value `True` or `False` accordingly. Hence the syntax:\n",
    "```Python\n",
    "sp.Eq(a*x**2 + b*x + c, 0)\n",
    "```\n",
    "is used to indicate that the left hand side (LHS) `a*x**2 + b*x + c` is equal to the right hand side (RHS) `0`.\n",
    "\n",
    "## Solving equations\n",
    "Having defined an equation, can we solve it?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of Eq(a*x**2 + b*x + c, 0) is x = [(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "quadSol = sp.solve(quadEq, x)\n",
    "print(\"Solution of\",quadEq,\"is x =\",quadSol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `sp.solve` statement solves the equation we have set up (`sp.solves`'s first argument) for the variable indicated (`sp.solve`'s second argument), and we see that it has returned the expected solutions for the quadratic equation:\n",
    "\n",
    "$$\n",
    "x = \\frac{-b + \\sqrt{b^2 - 4ac}}{2a}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "x = \\frac{-b - \\sqrt{b^2 - 4ac}}{2a}.\n",
    "$$\n",
    "\n",
    "We can find out from sympy how many solutions there are and get at the solutions separately as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution 0 is (-b + sqrt(-4*a*c + b**2))/(2*a)\n",
      "Solution 1 is -(b + sqrt(-4*a*c + b**2))/(2*a)\n"
     ]
    }
   ],
   "source": [
    "nQuadSols = len(quadSol)\n",
    "for n in range(0, nQuadSols):\n",
    "    print(\"Solution\",n,\"is\",quadSol[n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in order to use `sp.solve`, the RHS of the equation has to be zero. For example, if we were given the equation $ax^3 + bx^2 = cx$ to solve, we would first have to rewrite it as $ax^3 + bx^2 - cx = 0$, then we could determine the solution as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution 0 is 0\n",
      "Solution 1 is -b/(2*a) - sqrt(4*a*c + b**2)/(2*a)\n",
      "Solution 2 is -b/(2*a) + sqrt(4*a*c + b**2)/(2*a)\n"
     ]
    }
   ],
   "source": [
    "cubicSol = sp.solve(a*x**3 + b*x**2 - c*x, x)\n",
    "nCubicSols = len(cubicSol)\n",
    "for n in range(0, nCubicSols):\n",
    "    print(\"Solution\",n,\"is\",cubicSol[n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use sympy to solve the equation $4\\sin \\theta + 5 \\cos \\theta = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution 0 is 2*atan(2/3 - sqrt(10)/3)\n",
      "Solution 1 is 2*atan(2/3 + sqrt(10)/3)\n"
     ]
    }
   ],
   "source": [
    "theta = sp.symbols('theta')\n",
    "sinEq = sp.Eq(4*sp.sin(theta)  + 5*sp.cos(theta) - 1, 0)\n",
    "sinSol = sp.solve(sinEq, theta)\n",
    "nSinSols = len(sinSol)\n",
    "for n in range(0, nSinSols):\n",
    "    print(\"Solution\",n,\"is\",sinSol[n])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, as we have done above (e.g. `quadSol = sp.solve(quadEq, x)`) we can assign a symbolic expression to a variable. Here is another example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of a*x**2 + b*x + c = 0 is x = [(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]\n"
     ]
    }
   ],
   "source": [
    "y = a*x**2 + b*x + c\n",
    "print(\"Solution of\",y,\"= 0 is x =\",sp.solve(y, x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This does not define a function $y(x)$, it just assigns the symbolic expression `a*x**2 + b*x + c` to the variable `y`. The following expression will therefore produce an error message. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'Add' object is not callable",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-7-5eb93071928a>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"y(2) =\"\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: 'Add' object is not callable"
     ]
    }
   ],
   "source": [
    "print(\"y(2) =\",y(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to define a symbolic function, we can use the normal Python syntax for creating functions (having defined the relevant symbols):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of quad(x) = 0 is x = [(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]\n"
     ]
    }
   ],
   "source": [
    "x, a, b, c = sp.symbols('x a b c')\n",
    "#\n",
    "def quad(x):\n",
    "    f = a*x**2 + b*x + c\n",
    "    return f\n",
    "#\n",
    "print(\"Solution of quad(x) = 0 is x =\",sp.solve(quad(x), x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use this Python function to evaluate `quad(x)` for a specific value of `x`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad(2) = 4*a + 2*b + c\n"
     ]
    }
   ],
   "source": [
    "print(\"quad(2) =\",quad(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we get a symbolic answer, as `a`, `b` and `c` are symbols. We can check this, e.g. for `a`, by looking at the type of `a`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type of a is <class 'sympy.core.symbol.Symbol'>\n"
     ]
    }
   ],
   "source": [
    "print(\"Type of a is\",type(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we assign numerical values to `a`, `b` and `c`, `quad(2)` will return numbers for the solutions of the equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad(2) = 15\n",
      "Solution of quad(x) = 0 is x = [-1, -1/2]\n"
     ]
    }
   ],
   "source": [
    "a = 2\n",
    "b = 3\n",
    "c = 1\n",
    "print(\"quad(2) =\",quad(2))\n",
    "print(\"Solution of quad(x) = 0 is x =\",sp.solve(quad(x), x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The disadvantage of giving `a`, `b` and `c` specific values in this way is that we can no longer use them as symbolic objects. \n",
    "E.g. if we re-check the type of `a`, we will find that it is now an integer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type of a after assigning integer value is <class 'int'>\n"
     ]
    }
   ],
   "source": [
    "print(\"Type of a after assigning integer value is\",type(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an alternative method. Let's first turn `a`, `b` and `c` back into symbols and check that they are behaving as we want them to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type of a after setting as symbol is <class 'sympy.core.symbol.Symbol'>\n",
      "quad(2) = 4*a + 2*b + c\n"
     ]
    }
   ],
   "source": [
    "a, b, c = sp.symbols('a b c')\n",
    "print(\"Type of a after setting as symbol is\",type(a))\n",
    "print(\"quad(2) =\",quad(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now substitute the numerical value 3 for `a` using `subs(a, 3)`. We can check that `a` is given the value 3 when we do this, and that `a` is still a symbol, as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad(x).subs(a, 3) = b*x + c + 3*x**2\n",
      "quad(2).subs(a, 3) = 2*b + c + 12\n",
      "Type of a is <class 'sympy.core.symbol.Symbol'>\n"
     ]
    }
   ],
   "source": [
    "print(\"quad(x).subs(a, 3) =\",quad(x).subs(a, 3))\n",
    "print(\"quad(2).subs(a, 3) =\",quad(2).subs(a, 3))\n",
    "print(\"Type of a is\",type(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can substitute more than one numerical value into an expression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad(x).subs(a, 3).subs(b = -2.4).subs(c = 17/3) = 3*x**2 - 2.4*x + 5.66666666666667\n",
      "quad(2).subs(a, 3).subs(b = -2.4).subs(c = 17/3) = 12.8666666666667\n"
     ]
    }
   ],
   "source": [
    "print(\"quad(x).subs(a, 3).subs(b = -2.4).subs(c = 17/3) =\",quad(x).subs(a, 3).subs(b, -2.4).subs(c, 17/3))\n",
    "print(\"quad(2).subs(a, 3).subs(b = -2.4).subs(c = 17/3) =\",quad(2).subs(a, 3).subs(b, -2.4).subs(c, 17/3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also substitute another symbolic value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad(x).subs(a, x) = b*x + c + x**3\n",
      "quad(x).subs(a, y) = b*x + c + x**2*(a*x**2 + b*x + c)\n"
     ]
    }
   ],
   "source": [
    "print(\"quad(x).subs(a, x) =\",quad(x).subs(a, x))\n",
    "print(\"quad(x).subs(a, y) =\",quad(x).subs(a, y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solveset and the future of solve\n",
    "\n",
    "In time, the authors of sympy hope to replace `sp.solve` with a new routine, `sp.solveset`(see [here](https://docs.sympy.org/latest/modules/solvers/solveset.html) for a description), which is used in a similar but not quite identical way. For the moment `sp.solve` can cope with more types of equations than `sp.solveset`. An example using `sp.solveset` is shown below. Notice that it removes the restriction that the equation must be presented with the RHS being zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This quadratic equation is Eq(p*z**2 + q*z, -r)\n",
      "Its solution is {-q/(2*p) - sqrt(-4*p*r + q**2)/(2*p), -q/(2*p) + sqrt(-4*p*r + q**2)/(2*p)}\n"
     ]
    }
   ],
   "source": [
    "z, p, q, r = sp.symbols('z p q r')\n",
    "newQuadEq =  sp.Eq(p*z**2 + q*z, -r)\n",
    "print(\"This quadratic equation is\",newQuadEq)\n",
    "#\n",
    "newQuadSol = sp.solveset(newQuadEq, z)\n",
    "print(\"Its solution is\",newQuadSol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The curly brackets indicate that the solutions are not given as a list as with `sp.solve`, but as (for us) a new data type, a set. This allows `sp.solveset` to deal with situations where there is an infinite number of solutions, but means these have to be accessed in a slightly different way. If the number of solutions is finite, the following works.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution 0 is -q/(2*p) - sqrt(-4*p*r + q**2)/(2*p)\n",
      "Solution 1 is -q/(2*p) + sqrt(-4*p*r + q**2)/(2*p)\n"
     ]
    }
   ],
   "source": [
    "nSol = 0\n",
    "for sol in newQuadSol:\n",
    "    print(\"Solution\",nSol,\"is\",sol)\n",
    "    nSol += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fractions\n",
    "\n",
    "Fractions (rational numbers) can be entered as below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha = 1/137\n"
     ]
    }
   ],
   "source": [
    "alpha = sp.Rational(1, 137)\n",
    "print(\"alpha =\",alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They can be added, subtracted, multiplied and divided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "half + quarter = 3/4\n",
      "half - quarter = 1/4\n",
      "half*quarter = 1/8\n",
      "half/quarter = 2\n"
     ]
    }
   ],
   "source": [
    "half = sp.Rational(1, 2)\n",
    "quarter = sp.Rational(1, 4)\n",
    "print(\"half + quarter =\",half + quarter)\n",
    "print(\"half - quarter =\",half - quarter)\n",
    "print(\"half*quarter =\",half*quarter)\n",
    "print(\"half/quarter =\",half/quarter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, you don't have to assign variable names to fractions before manipulating them, and fractions can be turned into floats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans = -2290/1781\n",
      "ans as float = -1.2857944974733295\n"
     ]
    }
   ],
   "source": [
    "ans = sp.Rational(3, 137) - sp.Rational(17, 13)\n",
    "print(\"ans =\",ans)\n",
    "print(\"ans as float =\",float(ans))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiation\n",
    "\n",
    "Sympy provides a routine for differentiating. You must give it the expression you want to differentiate and the variable with respect to which the differentiation should be peformed, as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The differential of a*x**5 + b*x**2 + c/x**3 w.r.t. x is 5*a*x**4 + 2*b*x - 3*c/x**4\n",
      "The differential of quad(x) w.r.t. x is 2*a*x + b\n"
     ]
    }
   ],
   "source": [
    "print(\"The differential of a*x**5 + b*x**2 + c/x**3 w.r.t. x is\",sp.diff(a*x**5 + b*x**2 + c/x**3, x))\n",
    "print(\"The differential of quad(x) w.r.t. x is\",sp.diff(quad(x), x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second (and higher) differentials can also be calculated (in several ways)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Second differential is 20*a*x**3 + 2*b + 12*c/x**5\n",
      "Second differential is 2*(10*a*x**3 + b + 6*c/x**5)\n",
      "Second differential is 2*(10*a*x**3 + b + 6*c/x**5)\n"
     ]
    }
   ],
   "source": [
    "print(\"Second differential is\",sp.diff(sp.diff(a*x**5 + b*x**2 + c/x**3, x),x))\n",
    "print(\"Second differential is\",sp.diff(a*x**5 + b*x**2 + c/x**3, x, x))\n",
    "print(\"Second differential is\",sp.diff(a*x**5 + b*x**2 + c/x**3, x, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the output from the first of these methods is written differently to that from the latter two, but we can check they are the same by subtracting them and using `sp.simplify` to simplify the result!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(sp.diff(sp.diff(a*x**5 + b*x**2 + c/x**3, x),x) - sp.diff(a*x**5 + b*x**2 + c/x**3, x, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the first derivative w.r.t $x$ of the function:\n",
    "\n",
    "\\begin{align*}\n",
    "s(x) = A\\left( {\\tanh \\left[\\frac{x - x_0}{L} + \\frac{w}{2L}\\right] - \\tanh \\left[\\frac{x - x_0}{R} - \\frac{w}{2R}\\right]} \\right).\n",
    "\\end{align*}\n",
    "\n",
    "Using the values $A = 2.0$, $x_0 = 0.5$, $L = 1.1$, $R = 0.5$ and $w = 0.4$, plot the function and its derivative in the range $-2 < x < 3$ using numpy and matplotlib.\n",
    "\n",
    "*Hint, you can copy the output from the differentiation and paste it into your code for doing the plotting. It will need some editing, but not much! This can save you a lot of time and helps avoid typos.* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A*(-(1 - tanh(-w/(2*R) + (x - x0)/R)**2)/R + (1 - tanh(w/(2*L) + (x - x0)/L)**2)/L)\n"
     ]
    }
   ],
   "source": [
    "A, x, x0, L, R, w = sp.symbols('A x x0 L R w')\n",
    "def s(x):\n",
    "    s = A*(sp.tanh((x - x0)/L + w/(2*L)) - sp.tanh((x-x0)/R - w/(2*R)))\n",
    "    return s\n",
    "#\n",
    "s_ = sp.diff(s(x), x)\n",
    "print(s_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUddbA8e+hSJDepCOogNQEEgUVETs2FHYVFEUWFWHXurorrrvqrg3LvrqWVbHrKqjYEFBEAUHdVYqAQkBalFCkKYKCQHLeP85EAqbMTGbmziTn8zz3mUzm3jvnpsy5vy6qinPOOVecSkEH4JxzLrl5onDOOVciTxTOOedK5InCOedciTxROOecK5EnCueccyXyROFSgog8JiJ/CzqOkohIaxFREalSzOuLRKRPHN63j4jkluH4v4jIk7GMqdC53xGRi+Nxbpc44uMoXGlEJAdoDOQV+nY7VV0bp/cbClyqqr3icf54EZHWwCqgqqruKWXfW4HDVPXCGLxvH+A/qtqirOcqYxy3EqNrcsnFSxQuXGepas1CW1yShEus4ko/zhXmicJFragqDxHJEZGTQl/fKiKviMjzIrItVPWSVWjfliLyuohsFJHNIvKwiHQAHgOOEpHtIvJ9aN9nReT2QsdeJiLLRWSLiEwQkWaFXlMRGSEiy0TkOxF5RESkmGs4UkT+KyLfi8i6UAwHhHMuEaksIveJyCYRWQmcUcrPK0dEThKRvsBfgIGha1wQen2oiKwM/axWicjgYs5TPfTz+E5EFgNH7Pd6MxF5LfRzXSUiVxV67VYRGS8i/xGRH4Choe/9J/T6uyJyxX7nWyAiA0Jf/0tEVovIDyIyV0SODX2/uGuaISKXiki10M+4c6HzNhKRHSJyUOj5mSIyP7TfJyLStdC+N4jImtDPZqmInFjSz9rFlicKF2/9gHFAXWAC8DDYhywwEfgaaA00B8apajYwAvhvqORSd/8TisgJwF3AeUDT0DnG7bfbmdgHaHpov1OLiS8PuBZoCBwFnAj8PsxzXRZ6rRuQBfy2pB9EAVV9F7gTeDl0jekiUgN4EDhNVWsBRwPziznFLcChoe1U4Jc2ABGpBLwNLMB+picC14hI4es/GxiP/U5e3O/cLwHnFzpfR+BgYFLoW7OBDKB+aN9XRSStqGva75p/Bl4vfG7sZ/mhqm4Qke7A08DlQAPgcWBCKMG0B64Ajgj9bE4Fcor52bg48EThwvVm6E7vexF5M4LjPlLVyaqaB7yAfdgCHAk0A/6kqj+q6k5V/SjMcw4GnlbVeaEPoBuxEkjrQvuMVtXvVfUbYDr24fYrqjpXVf+nqntUNQf7gDpuv92KO9d5wAOqulpVt2DJqyzygc4iUl1V16nqomL2Ow+4Q1W3qOpqLMEUOAJopKr/UNVdqroSeAIYVGif/6rqm6qar6o79jv3G0CGiBwcej4YeD30c0ZV/6Oqm0M/r38C1YD2YV7fPkkIuCD0PbCk+7iqfqqqear6HPAz0BNL5tWAjiJSVVVzVHVFmO/pYsAThQvXOapaN7SdE8Fx6wt9/ROQFqoXbwl8XVqjbzGaYaUIAFR1O7AZu4Mu7n1rFnUiEWknIhNFZH2oKuZOrHRR0jUUnKsZsLrQa18TJVX9ERiIlabWicgkETm8mN1Let+DgWaFkvr3WJVQ40L7FD52/zi2YaWHgsQyiEKlDhG5TkSyRWRr6Nx1+PXPqzjTgOoi0iOUiDKwxFQQ93X7xd0SaKaqy4FrgFuBDSIyrnBVo4s/TxSuLH4EDix4EqpOahTmsauBVlJ0Y2ppXfHWYh8sBe9bA6uuWBPmexf2KLAEaKuqtbEP1SLbM4qwDvswK9Aqgvf91TWq6hRVPRmrTluClQQifd/VwKpCSb2uqtZS1dNLeu/9jAXOF5GjgOpYKYpQe8QNWImmXqhacCt7f14lnldV84FXsFLFBcDEUGIqiPuO/eI+UFXHho59KdQL7uDQ+9xdyjW4GPJE4criK6yEcIaIVAX+ilURhOMz7ANvtIjUEJE0ETkm9Nq3QIvCjcr7eQn4nYhkiEg1rBTwaajqKFK1gB+A7aE7+JERHPsKcJWItBCResCoCI79FmgdalNARBqLSL9Q0vsZ2M6+3ZH3f98bRaSeiLQAriz02mfAD6HG3+qhBvfOInJE0acq0mTsA/kfWJtDfuj7tYA9wEagiojcDNQu7pqK8RJWchrM3monsKQ4IlTakNDfxBkiUktE2ovICaHf9U5gB8X/bFwceKJwUVPVrVjD75PY3fyPQFgDv0JtFmcBhwHfhI4bGHp5GrAIWC8im4o49gPgb8BrWLI5lH3r4CNxPXZ3uw37sHo5gmOfAKZgDcfzsMbacL0aetwsIvOw/8XrsNLSFqydZP9G9QJ/x6qbVgHvYW0/wD4/14zQ65uw30+dcAMr1PB8Evt+mE8B3sFuEL7GPrQLV2Ptf01FnftT7O+kWehcBd+fg7VTPAx8BywHhoZergaMDl3LeuAgrOTnEsQH3DnnnCuRlyicc86VyBOFc865EnmicM45V6JAE4WIPC0iG0Tky2Je7xPqrz0/tN2c6Bidc66iC3pCsGexXg7Pl7DPLFU9M5KTNmzYUFu3bh1VQMu2LKNt/bZRHZuq/JrLv4p2veDXHKm5c+duUtUix0EFmihUdeZ+0y7EROvWrZkzZ05Ux2aNyWLO8OiOTVV+zeVfRbte8GuOlIgUO7NA4N1jQ4lioqp2LuK1Plhf+Vysf/n1xc1/IyLDgeEAaQ3SMjvd2SmqeLI3ZdOhYYeojk1Vfs3lX0W7XvBrjtTcy+fOVdWsIl9U1UA3bObQL4t5rTZQM/T16cCycM6ZmZmp0cp8PPpjU5Vfc/lX0a5X1a85UsAcLeYzNal7PanqD2oTvqGqk4GqIhLuBGTOOediIOjG7BKJSBPgW1VVETkS66W1OeCwnHMJtnv3bnJzc9m5c2fYx9yTcQ/Z2dlxjCr5hHPNaWlptGjRgqpVq4Z93kAThYiMBfoADcVWSrsFqAqgqo9hC8GMFJE92ERgg0JFJOdcBZKbm0utWrVo3bo1UvRihb+iG5UOjSpWG0Vp16yqbN68mdzcXNq0aRP2eYPu9XR+Ka8/TGhFNOdcxbVz586IkoQrmojQoEEDNm7cGNFxSd1G4ZxzBTxJxEY0P8ekbqNwLkjbtkF2NnzzjW27dkGzZtC8OXToYF87VxF4onCukG+/hTfegLfegg8+gN27i9/32GNh8GD47W+hQYPExeiC8eCDD/Loo4/SvXt3XnzxxdIPKEVOTg6ffPIJF1xwAQBz5szh+eef58EHHyzlyMTzqifngJUrYcQIaNUKRo6EZcvg6qstYcyfD1u2WAlj6VKYNg1uuw02brRjDj4Y/u//YE80q3+7lPHvf/+byZMnxyRJgCWKl17auy5UVlZWUiYJ8EThKrj162HoUGjbFp55Bn73O/jiC0sU994L/fpBejrUqwc1a0K7dnD88fDXv8LixTBvHvTpA9ddB1lZ8OmnQV+Ri4cRI0awcuVK+vXrR506dbjvvvt+ea1z587k5OSQk5NDhw4duOyyy+jUqROnnHIKO3bsAGD58uWcdNJJpKen0717d1asWMGoUaOYNWsWGRkZ3H///cyYMYMzz7Rp7bZs2cI555xD165d6dmzJwsXLgTg1ltvZdiwYfTp04dDDjkkYYnFq55chaT5lXj4YbjpJti5E669Fv74x8jaHUSgWzd4+21480246io45hgYMwaGDYtf7BXdNddYKa80P+0+mAPDHCqQkQEPPFD864899hjvvvsu06dP5+GHi++IuWzZMsaOHcsTTzzBeeedx2uvvcaFF17I4MGDGTVqFP3792fnzp3k5+czevRo7rvvPiZOnAjAjBkzfjnPLbfcQrdu3XjzzTeZNm0aQ4YMYX7oopcsWcL06dPZtm0b7du3Z+TIkRGNiYiGJwpX4SxdCktGP8OVX8PJJ8Mjj1iJIloi0L8/nHginHceXHIJrF1rScg76lQsbdq0ISMjA4DMzExycnLYtm0ba9asoX///oANeCvNRx99xGuvvQbACSecwObNm9m6dSsAZ5xxBtWqVaNatWocdNBBfPvtt7Ro0SJOV2Q8UbgKQxWeftru/H+W5owdCwMHxu7DvHZtK11ccgn87W+wZo0loUpewRtTJd35F7Z449d0bNQx5u9fpUoV8vPzf3leeLR4tWrVfvm6cuXK7Nixg2jGCBd1TEG31v3fY08CGsf8T9hVCNu2waBBcOml0LMndPzb+QwaFPs7/qpV4bnn4M9/hsces7YMV760bt2aefPmATBv3jxWrVpV4v61a9emRYsWvPnmmwD8/PPP/PTTT9SqVYtt27YVeUzv3r1/aTSfMWMGDRs2pHbt2jG8ish4onDl3rJllhxeew3uugveew8OqBfZyNRIiMDo0TB8uL3fM8/E7a1cAH7zm9+wZcsWMjIyePTRR2nXrl2px7zwwgs8+OCDdO3alaOPPpr169fTtWtXqlSpQnp6Ovfff/8++996663MmTOHrl27MmrUKJ577rl4XU5YvOrJlWvvvAPnnw9VqsCUKdaOkAgi8PDDsGqVJYyDD4YTTkjMe7v4yMnJ+eXr9957r8h9vvxy76rO119//S9ft23blmnTpv1q/w8++GCf53369AGgfv36vPXWW7/a/9Zbby32/eLJSxSu3HrwQTjjDGjdGubMSVySKFC1Krz6qnWp/c1vbKyGc6nIE4Urd/LyrMH66qvh7LPh448tWQShTh2YOBHy82HIEIvNuVTjicKVKz/+CAMGwEMP2diI8eOhRo1gY2rTxno/ffyxDeJzLtV4onCpIy8PNm2yEXJF2LwZTjrJ7uAfesim1ahcOcExFmPwYDj3XLj5Zvj886CjcS4ynihc8tqwAR5/3AY7dO1qRYNGjaB6dUhLsxbiyy+H997jmxW76dXLPoRffRWuuCLo4PclYt1lGzWCCy+E0MwOzqUE7/XkksuePfDKK/DUUzBjhlXuH3wwdOkCp54KLVpY/dLWrbBiBbz4IowZQy2pz+nVbuPsd0fQu09y3v/Ur29dZU89Fe6+G/brwOJc0gr0P0pEnhaRDSJSZB8vMQ+KyHIRWSgi3RMdo0uQnTvtlrt9e6un+eYb+MtfYMEC62P69ttWwX/11fb9u++G8eOZ/c4mLqg5gS+qduOfO/9A77/0gkWLgr6aYp1yinXXHT3a8pxLHQ8++CAdOnRg8ODBTJgwgdGjR4d97P4zxZZk6NChjB8/PqLYHnvsMZ5//vmIjolE0LdezwJ9S3j9NKBtaBsOPJqAmFwiqdpIuA4dbH7vhg1tQYilS20u765dix0+/cEHcPxpaXx60Fm0WDwVnn8evvrKZuobOzbBFxK+++6zrrPXXBN0JC4ShacZ79evH6NGjfrVPsVNpxFJoojUnj17GDFiBEOGDInL+SH4NbNnikjrEnY5G3hebeKT/4lIXRFpqqrrEhKgi68vv7R+rNOnW9XSe+9Za3QY82q8+aY1XbRrZ4c1bSpw6EVw2mnWanzhhZaEQovCJJNmzaza6frrraB01llBR+RKU3ia8WHDhlGvXj3mzJnDww8/zNChQ6lfvz6ff/453bt3p1+/flx99dWAzc80c+ZMRo0aRXZ2NhkZGVx88cVce+21v5xbVbnyyiuZNm0abdq02Weep7lz5/LHP/6R7du307BhQ5599lmaNm1Knz59OProo/n444/p168f27Zto2bNmrQ/qj1DzxjKZ599BliC6tev3y/TlEcr2dsomgOrCz3PDX3vV4lCRIZjpQ7SGqSRNSYrqjfM3pQd9bGpKtHXXDlP+d0767hk8np+TKvEoxe04s1eVclbdSM8cWOpx2/+3+nkPHczNQ7OpvKwqznr7R/2eT1tQB4PrD2QbhcO5pYPbuLdHr9efi7o37MeWJm0pi/x22HV6HTLQCod8HNc3y/o6y2rezLuQTfaB2jjv95F2pdLSz2mpebzo4RXabKzc3u+vb34v72rbruKtye/zWOvPka9BvV4Y9wbbNmxhcUbF/P9zu/J+TKHh8Y+ROXKlfn9hb/n+tuvp3uP7vy4/UdWbV/FiFEjePbfz/LvF/8NwOKNi38599SJU5n3xTxenvYymzdupl+vfpz825NZsHYBl464lIeef4j6Devzzpvv8Ifr/sDt/7qdn3b/xMp1K3l0vFWyPHLPI/zIj5xwyAn88NMPTJk9hZatW/LkU0/S54w++7wfwPpt67lozEVh/Wwg+RNFUbeWRU7FqKpjgDEAWVlZOmf4nKjeMGtMFtEem6oSes0LFthKQfPXwfnnU/fBB7mxYUNKTw/m4YfhymdsOoy33upCzZq/nhYBgGE/wplncvuzM7l9wMM2RLuQZPg9Tz/cruPcbR9zY7g/gCglw/WWRXZ2Nh0adbAn1RtA1QNLPebH3T9RI4z9AGpUb0CDUmaarVqpKu0btqdhw4Z8Vusz1lVfR8dGHambVpf+g/vTpUkXAPoe35eHbnuIwYMHM2DAAFo0bcGWuluoeUDNImezHbNgDJdefKkd3wROOvEkWtZuSeUtlVmxdAVXnG9d+PLy8mjatCkdG3XkwKoHMnLoyF/O16hGI2rWrElalTSGXDCEzz/4nFNHncqMSTN4+eWXadto33n0ZZP86u9BLi++JJ/siSIXaFnoeQtgbUCxuLLIz7f5oUeNsu4/b7wB55wT9uGqcOedNhvr2WfDuHHWQ7ZYNWrYgIpjjtm7bF3jxmW/jhg6/nirdrr7bpsPytfdDlOY84x/vXFxXKYZL0qNQqM6R40axRlnnMHkyZPp2bMn77//fqnHSxHVrapKp06d+O9//1vqexY2cOBAzj33XAYMGICI0LYsi62EBN2YXZoJwJBQ76eewFZvn0hB69db28F119md/aJFESeJP//ZksSFF9o4iTDWfrFk8dJLNsf4sGF2oiRz553www82y6wrH1asWEGXLl244YYbyMrKYsmSJaVOKT5u3Djy8vJYt24d06dPB6B9+/Zs3Ljxl0Sxe/duFoXRo+/QQw+lcuXK3HbbbQwcODAm1xR099ixwH+B9iKSKyKXiMgIERkR2mUysBJYDjwB/D6gUF20PvjAFp2eNcu6v77+ekS3znl5drd9333whz/YWg8RrfrYsaN1q508Gf7978jjj7POneHii61K7Ztvgo7GxcIDDzxA586dSU9Pp3r16px22mklTinev39/2rZtS5cuXRg5ciTHHXccAAcccADjx4/nhhtuID09nYyMDD755JOwYhg4cCD/+c9/OO+882JzUapa7rbMzEyNVubj0R+bquJyzXv2qN56q6qIaseOql9+GfEpdu5UPfdcVVD9619V8/OjjCU/X7VvX9W0NNXFi1U1uX7PX3+tWq2a6tCh8XuPZLreaCwO/d4isWjDojhEktzCveaifp7AHC3mMzXZq55cKtq40aqabr3V6oo++ww6dYroFNu3W/39q6/CP/9pQyqiXo1OxIZEp6XBn/4U5Unip1Urm3Lkueesx7BzycYThYutzz6DzEyYORPGjLFPvwinby2Y3G/aNPt8/+MfYxBXkybWkD5pEnz0UQxOGFs33gi1avm0Hi45eaJwsaFqbRDHHmtTtn78MVx2WcTFgNWr7RTz59uA7aFDYxjjlVdC06aWMJKsYbtBA5ud5LXXoIxjo8otTbLfWaqK5ufoicKV3Y4d1qto5Ejr8zlnjpUqIrRoERx1FKxZY8uWnn12jOM88ECb5/vjjznmyx9K3z/BrrkGateGf/wj6EiST1paGps3b/ZkUUaqyubNm0kLq9vgXsk+jsIlu5wcW+dz3jz429/glluiWgTio4+sTaJ6desg1bVr7EMF4JJL4L77+MOba+Bf+VApee6V6te3UsVtt1mpIm4/gxTUokULcnNz2bhxY9jHrN+2HtkUbcNWagrnmtPS0mjRokVkJy6ulTuVN+/1FJmor3nKFNX69VXr1FGdMCHq93/lFev1066d6qpVUZ8mfC+9ZF2pXnopAW8Wmc2bVWvXVv3Nb2J7Xv+7rhjKcs14rycXU/n5NlKsb1+b4W727KhmtlO1Hk3nnQdZWfDJJwla23rgQFY2TQt7hG8iFZQqXnvNZjtxLhl4onCR2brVqppuugkGDYL//Q+imCJgzx5rW77+epvs9f33EziFRaVKjD+ukfXQmpN88x9de621VdxxR9CROGc8UbjwLVxot/4TJ9rd+IsvRtz1FSzXnHUWPPKIzepR6rxNcTCpZwOLPQlHa9erZ6PQx4+H7Oygo3HOE4UL1wsvQM+etgzp9OlWPxLFCLicHJun7/33bZjFffcF0578Y/XKNhhw7FjYsiXxAZTi2msteUawiJpzceOJwpVsxw6bbGnIEOjRw3o39eoV1almzIAjjtjb/fWyy2IbasRGjrQlWJ99NuBAfq1RI7j8ciu0rVoVdDSuovNE4Yr31Vc2sOGJJ2yQ2tSpNsI5Qqo26d1JJ9lKp59+auswBC493Yo3jz5qDfRJ5vrrrafx3XcHHYmr6DxRuKKNHWuD5nJzbebVu+6CKpEPu9mxAy691BquTz/d2r7btYtDvNH6wx9g+XJLgkmmeXNbSuOZZ6wU5lxQPFG4fW3fbp9OF1xgd9yff24T/EVh5Uq7YX/6aVtL4s03oU6dGMdbVgMGWD3PU08FHUmRbrjBplq/776gI3EVmScKt9e8edar6bnnbJT1jBnQsmWphxVl0iQrkKxaBRMm2GjjJBoEvVe1avDb31pPrh9/DDqaX2nTBgYPhscft0l5nQtCMv7rugSrlK/WvaZnT1sN7oMPbMKhKKqadu2yuvUzz7QPublzoxqLl1gDB1od2aRJQUdSpBtvtDb3JBwf6CoITxQVXU4Oj/3zK/s0OvtsW1v6+OOjOtXy5VbV9M9/WtX/J5/AIYfEON546NXLGulfeSXoSIp0+OE2xvHhh+H774OOxlVEQS+F2ldElorIchEZVcTrQ0Vko4jMD22XBhFnuaRq9RldutAu9yerbnrlFZtDIopTPfMMdOtmyeL11+1DLdGD6KJWubJVP02aZG00Segvf7G1tR95JOhIXEUUWKIQkcrAI8BpQEfgfBHpWMSuL6tqRmh7MqFBllfffAOnngojRkCPHgy6uaONk4hiAN2GDdC/v80ynplp60j07x+HmOPtvPOsfuftt4OOpEjdulmvsfvvT8qmFFfOBVmiOBJYrqorVXUXMA6I9QoErrC8PHjwQejY0eqFHn0Upk5lfYNqUZ1u/Hjo0gXeecd65UybBgcfHOOYE+WYY2yCwyStfgKbXmvzZhvR7lwiBZkomgOrCz3PDX1vf78RkYUiMl5EouuC46zt4ZhjbOqN3r1tlaARI6IqRaxfb3Xm554LLVrYvHrXXZekvZrCVamSXdA771gdTxI6+mhrPrrnHmt7dy5RRANaMUpEzgVOVdVLQ88vAo5U1SsL7dMA2K6qP4vICOA8VS1yTK+IDAeGA6Q1SMvsdGenqOLK3pRNh4Ydojo2GdXYkcfwt9cycPoGth1YhfsGtmTKEfX2SRDhXrPmC5s/6Ufua1eRvyuNZmeNofHJ/0Eq58XzEuKiqGvuumI7T9+zlL8Oa827PRI1lW1kti3N5Kv/e5yWA+/loBNeDvu48vZ3HQ6/5sjMvXzuXFXNKvLF4haqiPcGHAVMKfT8RuDGEvavDGwN59y+cJGq5ufbwjxNm6qKqA4frrppU5G7hnPNCxaoHnWUrfdz7LGq2dmxDjixirzmvDzVFi1Uzzkn8QFFoHdv1WbNVHfsCP+YcvN3HQG/5siQpAsXzQbaikgbETkAGARMKLyDiDQt9LQf4JMuh2P2bKtmuuACaNrU5s14/PGoFnzYsgWuugq6d4dly2z+vA8/tC6b5U6lSnDGGTaOZNeuoKMp1s03w9q1NuLduUQILFGo6h7gCmAKlgBeUdVFIvIPEekX2u0qEVkkIguAq4ChwUSbInJy4KKL4Mgjbf6Mp56yxXmOPDLiU+3ebV0x27a1x+HDYelSuPjiqJo1Ukffvjbo8L//DTqSYp1wgrVX3HUX/Pxz0NG4iiDQ5kdVnayq7VT1UFW9I/S9m1V1QujrG1W1k6qmq+rxqrokyHiT1qZNtoBB+/bWFWnUKJv5ddgwGyMQAVVbhrNzZ7jiCsjIsC6v//53VEMsUs8JJ9iI9HffDTqSYonALbfYfI1JOEO6K4dSuZ+K27LF+ky2aWPdXocMsfqhu+6ytTQjoGo1LkcdZWPPqlSxSfzef9+6wFYYtWtbtV0SJwqAk0+2GVfuuMOGfzgXT54oUtHGjZYgWre2pHDGGfDll7ZuRIsWEZ1KFX7IPoLevW29iDVrrMZqwQKb0aNcVzMVp29fK0atWxd0JMUSgdtvh9WrrfnJuXjyRJFKvvnGWpYPPhjuvNNGVy9caItOd4isS1x+vs3q2qsXLHvgUVatsmk3li+3Gqso5gMsP/r2tccpU4KNoxQnnmg1ZXfembQzj7hywhNFKvj0Uxg0yGbYe/RRm+108WJ49VVrTIjAzp3WW6ZLFysxrF0LLQfdzfLlNpFftegGaZcv6ek2SWCSVz+BVT1t2GA1j87FiyeKZLVzJ7zwgjUa9OxpI4avuQZWrLAZ+CIsQaxZY0tMtGoFl1xiJYYXX7QmjYOOfzV1JvBLBBErrb33nk17ksR69oR+/Wy09nffBR2NK688USSb7Gz405+srWHIEJvc58EHrYvLfffZJ32Y8vIsv5xzjtVW3XGH5Z1p06wK/oILKngVU0lOO80+eWfPDjqSUt12m806cu+9QUfiyiv/mEgG331n1UjPPGOD46pUsXqhkSOtEjrCFuVly2zW8Oeft8bORo1sMaHhw1NkfYhkcNJJNgDv3Xfttj2Jde0K559vM8uOGBHRvYRzYfESRVB++snGPPTvb/Xhl19ut4X33Welh/HjrbUyzCSxfr0VPI46Ctq1s85QnTrByy/b6UaP9iQRkQYN4IgjrPopBdx1lz3eeGOwcbjyyUsUifTDD1YXNH48TJ5syaJJE/j9760eKCsrotJDbq6NdXj9dZtWIz/f7i5Hj7YB2if8gPkAACAASURBVM2axfFaKoI+fWy5vp9+ggMPDDqaErVqZaXG22+HK69M+kJQyfbsgc8/t1kFli+3drmvv7YR89u32++jShWoWtV6X9SrZ4m9YUNo3tzWeW/Vig45P8LWrVCnTtBXlPI8UcTbihWWHCZMgBkzbG6MJk1g6FCbq/u448IePZ2fb/8/kybBxIl7q887drRhFYMG2dcuRnr3hrvvturAE4qctDip3HCDjYG55hqbgSSlxsD8+KMVf19/HWbN2jvV+4EHwqGHWiNb3bpQs6Z9b88e+1/audOqbjdvtqnzp0z5pa/wCwB31bUEcvjh1gGkQwcranfpYv+HKfVDCo4nilj77ju7vX//favfXrHCvn/44fYffNZZNlFPmMlh/XqYOnXvtn69/W0feaQ1Tg8YUE4n6EsGxxxj7RQzZ6ZEoqhZ08ZU/O53MHasFVKT3tKl8MAD1gVv2zZLCuefb6W5Xr2shBDJh7mqlSK++YbrHhvAP9tcbo12S5ZYEtq8ee++DRpYwkhPt6J4166WRKpXj/llpjpPFGW1eTN89JHdBc2YAfPm2R/rgQfah8s119gArsMOC+t0ubnw8cd2qhkz7O8b7G/6pJNsOcy+feGgg+J1Qe4XderYZFczZwYdSdiGDIGHHrKOc2eeGfFMLomzeTP8/e82iVjVqrYU7WWXWXIuy12+iJU86tblw4y6MPxP+76+caPNYvDll7aY18KFNqPBTz/Z65Uq2ZxpXbvuTSDp6ZEnrHLGE0Ukdu2y4u1nn1l1xP/+t/eTvFo16NHDZms74QS75S9l9NqOHdZN9bPPbEzdxx/b4GuAWrXg2GPt7vCkk+zzKqVXkEtVvXvDY4/Z7/6AA4KOplSVKtln71FHWcP2I48EHdF+VK133/XX253/5ZdbwmjUKDHv36iRLRN4/PF7v5efbyX/hQttW7DA/iFfLrQwVN26ljS6dLFBrp07W+mjXr3ExB0wTxTF2b7d7jgWLLBP83nz7I+oYF7nRo2sxXDIEPtEz8qipFFr33+/93Tz5tm2aJFVtYLdsBx9NPzxj3ZTlZHhYxySQu/eVjUyZ479glJAjx7WoP3QQzB4cBKFvWOHTUn89NNWtfTQQxHPLBAXlSrZfPpt21q7YYGtW+1/vqDk8cUX1ud827a9+zRtagmjY0dr/yhoCznooHJVAvGPogK7d8Pf/sb/TV4Od7WxtR0K1KkD3brZf19Wlm2HHPKrPwRVm/F76VLbFi+2bdEiG89QoFEjyMy0aqQjj7TNeyglqV697HHmzCT6xC3d7bdbj7jLLrObksDl5NiH8Lx5NkXALbdEPAV+wtWpYzeBxx6793uq9s9cUH1V8E/+1FPWIF/42HbtbDvssL3boYda43qKJRFPFAWqVoX//Iem7IJje9o8FwX1k61a/fKL3b3b2hFyZsCqVbY+0IoVti1bZiWHAtWr2w1G7957S61du1pSSLG/k4qrUSO7W5w509b5SBG1atm0YGecYR23aBJgMMuX2z/BTz/B229b40mqErHPg1at7E6vgKp9MCxZYrMrfPWVbR99BC+9ZK8XqFXLlgZo08Z6c7VuvfecLVtaaSTJ6pk9UYSowg9ffEP/h87loWPHsm4drP0K1ky3eZJWr7Zt7dp9f+eVK9vvuqCzRvv2dhPRvr39/pPs9+2i0bu39crJy0v+u+BCTj/d/iZvvx0OvT6yucFiJifH2ux277YPzWSoaooHEfuQb9nSFgsp7Oef7Y5y+XJ7LLx98MGvp/6tWtXqogu2Zs2siqtpU+vS26QJNG5sPVwS9PfoiSJEFRo1rsTu3a9RuCNkjRr2u2rRwn7/LVtaYmjTxhJBy5b2e3XlWEGD9oIFtnh4Cnn4Yft8XvXknfwwKsG9oHJzbXaBbdtg+vTymyRKU63a3jEc+1O1Bci++cbuRL/5xn5ua9bY4/z5Nji3qHnkRSxZNGpkW8OG/PH71TA89pcQaKIQkb7Av4DKwJOqOnq/16sBzwOZwGZgoKrmxCOWSpVsCoy7Z/+VpwbfTtOmlshr1/ZqogqvoI565syUSxT169uYil69mzJihBWMEvL3/OOP1o9740a7a87ISMCbpqCCD/sGDawdtDjbttkgqm+/tQW1NmzYd9u0CZYsoVucFiYJLFGISGXgEeBkIBeYLSITVHVxod0uAb5T1cNEZBBwNzAwXjGNGAFPVnqXE064PV5v4VJRixbWeeHDD21cTIo55hhoduYTjB07kpNPti7XcaVq09IsXmxzZR1xRJzfsAKoVcu2tm1L3O2iMVnMicPbB1mDfiSwXFVXquouYBxw9n77nA08F/p6PHCiiN/fuwAceyx88sm+DVQppMlpz3D88bY4VdxnTn/mGetGevPNNgjIpbwgq56aA4U6jZIL9ChuH1XdIyJbgQbApv1PJiLDCdXOpTVII2tMVlRBZW/KjvrYVOXXXLrf7N7IjRs2cNZdXVnXMPWWAVyyJZvDTj+FvC+e5ugTq3P4qN9RreHamL/PoWt28Nxd2Sw4vBZXNnmb/DETY/4e4SqPf9f5u6uSt7MG+TtrkLezBnk7DyT/5wPJ/7k6eT9XZ9327mQR+2sWDegOSUTOBU5V1UtDzy8CjlTVKwvtsyi0T27o+YrQPpuLOmeBrKwsnTMnugJY1pgs5gyPR+Etefk1h2HePBv8Mm6cLUWbYgqud8kSG7XdpIkVkGI6sHjXLqtn37LFGmEbN47hySOXrH/Xe/ZYk0LBtnmzbVu22Pbdd7Z9//3ebetWmydx166Sz12l1hZ2/1A/qrhEZK6qFpllgixR5AItCz1vAex/i1OwT66IVAHqAFsSE55zhXTpYiPvP/00JRNFgcMPt4F4p5xiS6FMnhzDGdTvv9/aJSZNCjxJJFpB56U1a2xbt8629ev3tkF/+6217Ze0ZG21atYBoV49mzWkcWPral+njm21a9tW0GRRs6ZtBV+f8/q5wAcxv74gE8VsoK2ItAHWAIOA/ee7nABcDPwX+C0wTYMqArmKrWpV6/H06adBR1Jmxx0Hzz5r03v07WtT1pe52+zq1bYm69ln7zsQrZzYtct6rubk2Pb11/a8oFdrbu7e2X0Kq1Nn77CH9HQbS3fQQb/0ZqVhw72dnurXL/vEtVVqbi3bCYo7b1zOGoZQm8MVwBSse+zTqrpIRP4BzFHVCcBTwAsishwrSQwKKl7n6NHDhjvv3p3yg2fOP9+6hF94oY2He/dd+9CK2nXX2YDEBx6IWYyJtnOnzbDw1Vc2y8Ly5XvXTcrNtbkDC1SqtHeNpKwsK50VHiNXMDauvMxYHug4ClWdDEze73s3F/p6J3BuouNyrkg9elj1ysKF1l6R4gYOtOqK3/7WShlvvRX2bPj7ev99W/P9H/+wUahJbuvWvXOwZWfbtmSJlRQK11c0amQ/j2OPtd7Rhxxil9e6tSWDFL9XiIiPzHYuXD1CnfI+/bRcJAqwuaDefdcWwOre3ZZmiKgJZs8emyzz0ENtEYwkonmVf5n4tWDy1y++sNJBgbQ0awPo0cMmgi6Ygueww3wF1cI8UTgXroMPtgrmTz+1AWXlxHHH2RK7gwbZNn063HuvNZCWatw4ux1/7bUSp9mPt23brKPV/Pl2LfPnw/wvZtE1NI1/1ao2g8Zxx+1dSqJTJ5+PLVyeKJwLl4jdepaDBu39tWplA8//+le45x6rhrrzTrj44hI+SPPybD3ezp3hnHMSFut331kymDt379ouy5btrTZq1Mh66R50wjj+OWQI6elWUqhIVUWx5onCuUj06GFTZX//vfVfLEeqVrUpyQcMsJlKhg2zSQVHjbI88KsP2tdes9LEyy/H7bZ8wwZLCgUJYd48m3S1QMuWVgt44YWWHLp12zuNf9aYBxk8eEhc4qpoPFE4F4mCdorZs389nXQ50aOHDcYbO9ZKGOedZ714Lr0Uzj3XChCi+TZ/+eGH77sqXJR27bJSQUGbQsHCkmsLjaw65BBrR7n0UksO3buXsaeWC5snCuciccQRdrv66aflNlGAXeIFF1jD9rvv2jrct99uQyWaNYMbD3+LK774giU3vUCj7ytTv37ps9L+8IM1JOfmWg+jgu6nS5ZYkihYFrhKFWtPOP54SwbdutnksxVkeeqk5InCuUjUqWN30eWwnaIolStbz6gzzrARx1OmwLvvKL3euJ3lHErnOwaRd8feEcX169vXqrbt2LF3Gor9B6QdcIB1lmrXzsYhdOxoDcwdOtg5XPLwROFcpHr0sLkvVCvUYiXNm1u7xbDDZsH4eWwZPYaJ6VVYvNimp9iyxeYs2r3b9hexAWd169rWqJHN2N68uTWet2iRUgsGVmieKJyLVGamzYGxdq196lU0jz0GdepQ/8rB9D3QpgFx5Zv3IHYuUgWr3M2bF2wcQdiwAcaPt36zMZtN0CU7TxTORSo93epVKmKiePZZq1u6/PKgI3EJ5InCuUjVqGEN2nPnBh1JYuXnw+OP2/Dmjh2DjsYlkCcK56LRvXvFK1FMnWqj3UaMCDoSl2CeKJyLRvfu1l/022+DjiRxHnvMui717x90JC7BPFE4F42CBu3PPw82jkRZt86mLhk2zAc5VECeKJyLRrdu9lhRqp9eecUmARw6NOhIXAA8UTgXjTp1bNGCitKgPXaszaNx+OFBR+ICEEiiEJH6IjJVRJaFHoucxUVE8kRkfmibkOg4nStRRWnQXrnSpiw5//ygI3EBCapEMQr4QFXbAh+Enhdlh6pmhLZ+iQvPuTB0726z223ZEnQk8TVunD1GtPSdK0+CShRnA8+Fvn4OSNyqJ87FSkVp0B43Do4+2lb4cxWSaOHVxBP1piLfq2rdQs+/U9VfVT+JyB5gPrAHGK2qb5ZwzuHAcIC0BmmZne7sFFVs2Zuy6dCwQ1THpiq/5ujU2b6HD65bwL8GNOeFU5vEKLL4iPZ6D1m7g1f+vph7BrXkleMPikNk8eN/15GZe/ncuaqaVeSLqlriBlwB1CttvyKOex/4sojtbOD7/fb9rphzNAs9HgLkAIeG896ZmZkarczHoz82Vfk1l0GrVqqDBsXmXHEU9fXedJNqpUqq69fHNqAE8L/ryABztJjP1HBmj20CzBaRecDTwJTQSUukqicV95qIfCsiTVV1nYg0BTYUc461oceVIjID6AasCCNm5xKje/fy2/NJ1Xo7nXgiNG4cdDQuQKW2UajqX4G2wFPAUGCZiNwpIoeW4X0nABeHvr4YeGv/HUSknohUC33dEDgGWFyG93Qu9rp3t+XZtm0LOpLYK1igetCgoCNxAQurMTtUglgf2vYA9YDxInJPlO87GjhZRJYBJ4eeIyJZIvJkaJ8OwBwRWQBMx9ooPFG45JKebo9ffBFsHPEwYQJUqgT9vMNhRVdq1ZOIXIXd9W8CngT+pKq7RaQSsAz4c6RvqqqbgROL+P4c4NLQ158AXSI9t3MJVZAoFiywnkHlyYQJdk0NGwYdiQtYOG0UDYEBqvp14W+qar6InBmfsJxLEa1a2TqfCxYEHUlsrV4N8+fD3XcHHYlLAqUmClW9uYTXsmMbjnMpRsRKFfPnBx1JbE2caI9e7eTwuZ6cK7v0dGujyMsLOpLYmTDB5rJq3z7oSFwS8EThXFmlp8NPP8GKctJze/t2mDYNzjrLSkyuwvNE4VxZZWTYY3mpfpo6FXbt8mon9wtPFM6VVceOULly+WnQnjDBGuiPOSboSFyS8EThXFmlpdk6DeUhUeTlwaRJcPrpULVq0NG4JOGJwrlYyMgoH4li9mzYuNHaJ5wL8UThXCykp0NuLmzeHHQkZTNlijVgn3JK0JG4JOKJwrlYKDxCO5W99x5kZUH9+kFH4pKIJwrnYqE8JIqtW23JUy9NuP14onAuFho3hiZNUruL7IwZ1ph98slBR+KSjCcK52IlPT21SxRTp0KNGnDUUUFH4pKMJwrnYiUjAxYvtsFqqei996BPHzjggKAjcUnGE4VzsZKeDrt3w9KlQUcSua+/tgWYvNrJFcEThXOx0rWrPaZi9dPUqfboicIVwROFc7HSrp1V2yxcGHQkkZs6FZo1gw4dgo7EJaFAEoWInCsii0QkX0SyStivr4gsFZHlIjIqkTE6F7GqVaFTp9RLFHl58P771i3WZ4t1RQiqRPElMACYWdwOIlIZeAQ4DegInC8iHRMTnnNR6to19aqePv8ctmzxaidXrEAShapmq2ppLX5HAstVdaWq7gLGAWfHPzrnyqBrV1i/HjZsCDqS8E2bZo8n/moZe+eA8NbMDkpzYHWh57lAj+J2FpHhwHCAtAZpZI0ptkarRNmbsqM+NlX5NcfOEat/4FHg93f14rMOtWN+/miVdL3/emEZTZqmMfCtMxIcVXz533UMqWpcNuB9rIpp/+3sQvvMALKKOf5c4MlCzy8CHgrnvTMzMzVamY9Hf2yq8muOoQ0bVEH1n/+Mz/mjVOz17t6tWrOm6siRiQ0oAfzvOjLAHC3mMzVuJQpVPamMp8gFWhZ63gJYW8ZzOhdfjRrZVB6p0qA9b54tfdqnT9CRuCSWzN1jZwNtRaSNiBwADAImBByTc6VLpak8Zsywx+OOCzQMl9yC6h7bX0RygaOASSIyJfT9ZiIyGUBV9wBXAFOAbOAVVV0URLzORaRrV5vKY/fuoCMp3YwZNnaiceOgI3FJLJDGbFV9A3ijiO+vBU4v9HwyMDmBoTlXdl272nxPX31l4yqS1Z49MGsWXHRR0JG4JJfMVU/OpaZUWZvC2ydcmDxROBdr7dvbKO1kb9D29gkXJk8UzsXaAQdYvX8qJApvn3Bh8EThXDwke8+ngvYJr3ZyYfBE4Vw8dO0Ka9fCpk1BR1I0b59wEfBE4Vw8FDRoJ2v104cf2mPv3sHG4VKCJwrn4qEgUcyfH2wcxZk1y9bPaNIk6EhcCvBE4Vw8HHQQNG2anO0U+fnw0Udw7LFBR+JShCcK5+IlIyM5SxSLFsF333micGHzROFcvKSnQ3a2jdJOJrNm2aMnChcmTxTOxUtGhs33tHhx0JHsa9YsaN4c2rQJOhKXIjxROBcvyTiVh6olimOP9fWxXdg8UTgXL23bQvXqydVOkZMDa9Z4tZOLiCcK5+KlcmXo0iW5EsXMmfboicJFwBOFc/GUkWFVT7acb/BmzYJ69ZJ7+nOXdDxROBdP6enWFXX16qAjMbNmQa9eUMn/9V34/K/FuXjKyLDHZGjQ/vZbW0zJq51chIJaCvVcEVkkIvkiklXCfjki8oWIzBeROYmM0bmY6NLFHpOhneKjj+zRE4WLUCBLoQJfAgOAx8PY93hVTdIpOJ0rRa1acNhhyVGimDXLemF17x50JC7FBLVmdjaAeD9uVxFkZMDnnwcdhfV46tnTFlZyLgLJ3kahwHsiMldEhgcdjHNRSU+HFStg27bAQqixI89KNT6tuItC3EoUIvI+UNQcxjep6lthnuYYVV0rIgcBU0VkiarOLOb9hgPDAdIapJE1ptimjxJlb8qO+thU5dccX73Wf88DwCW3HcGCw2om5D331+LzbMjPZ+TWl5g9ZmIgMSSa/13HkKoGtgEzgKww970VuD6cfTMzMzVamY9Hf2yq8muOs7VrVUH1gQcS9577eeq0JqpVqqhu3x5YDInmf9eRAeZoMZ+pSVv1JCI1RKRWwdfAKVgjuHOppWlT2+YE13Gv27Lt1ohdo0ZgMbjUFVT32P4ikgscBUwSkSmh7zcTkcmh3RoDH4nIAuAzYJKqvhtEvM6VWVYWzJ0bzHvv3EmnnB+9W6yLWlC9nt4A3iji+2uB00NfrwTSExyac/GRmQkTJ1qDdq1aiX3v2bM5YI96Q7aLWtJWPTlXrmRl2XxPQQy8K1io6JhjEv/erlzwROFcImRm2mMQ7RQzZ7KiWRo0aJD493blgicK5xKhSRNo1izx7RR5efDJJ3weULdcVz54onAuUYJo0F6wALZt4/O2CW4XceWKJwrnEiUzE5YuTewI7VD7xHwvUbgy8EThXKIUNGgnct6nDz+E1q35tr7P7+Si54nCuURJdIN2fr4liuOPT8z7uXLLE4VzidK4MbRokbh2ioULYcsWTxSuzDxROJdImZmJK1FMn26Pffok5v1cueWJwrlEysy05Uh/+CH+7zVjBhx6KLRsGf/3cuWaJwrnEikrNAV0vKuf8vK8fcLFjCcK5xKpRw97/OST+L7P/PmwdasnChcTniicS6T69aFzZ/joo/i+j7dPuBjyROFcovXqZSWKvLz4vceMGdCunU0b4lwZeaJwLtF69bLG7C/jtA7Xnj0wc6ZXO7mY8UThXKL16mWP8ap+mjfPpgnxROFixBOFc4nWqpUNvItXoihonzjuuPic31U4niicSzQRK1XMmmVzP8Xa++9Dx442tblzMRDUmtn3isgSEVkoIm+ISN1i9usrIktFZLmIjEp0nM7FTa9esGYNfPNNbM+7fbu1T5x2WmzP6yq0oEoUU4HOqtoV+Aq4cf8dRKQy8AhwGtAROF9EOiY0SufiJV7tFNOmwa5dcPrpsT2vq9ACSRSq+p6q7gk9/R/QoojdjgSWq+pKVd0FjAPOTlSMzsVV585Qu/be9axjZfJkqFlzbyJyLgaqBB0AMAx4uYjvNwdWF3qeC/Qo7iQiMhwYDpDWII2sMVlRBZO9KTvqY1OVX3Mw/tUSmrz9HAO7x2iSQFUmvvIFiw+rwZ+fPXqfl5LhehPNrzmGVDUuG/A+8GUR29mF9rkJeAOQIo4/F3iy0POLgIfCee/MzEyNVubj0R+bqvyaA3L77aqgunlzbM73xRd2viee+NVLSXG9CebXHBlgjhbzmRq3EoWqnlTS6yJyMXAmcGIoyP3lAoWnvWwBrI1dhM4FrHA7Rb9+ZT/f5Mn26A3ZLsaC6vXUF7gB6KeqPxWz22ygrYi0EZEDgEHAhETF6Fzc9exp7QkFH/BlNXkypKdD8+axOZ9zIUH1enoYqAVMFZH5IvIYgIg0E5HJAGqN3VcAU4Bs4BVVXRRQvM7FXrVqcOqp8PbbtmxpWWzdaiUT7+3k4iCQxmxVPayY768FTi/0fDIQo9st55JQv37w2ms27UZWGRohp061SQY9Ubg48JHZzgXp9NOhUiUrVZTFpElQt65VZzkXY54onAtSw4Zw9NEwoQzNbzt2wBtvwFlnQZVk6PHuyhtPFM4FrV8/W5Eu2uk8JkywNoqLL45tXM6FeKJwLmhnnWWPEydGd/xzz0HLlj6tuIsbTxTOBa19e2jbNrrqp3XrYMoUuOgia+twLg78L8u5oIlYqWL6dFtwKBIvvmhda4cMiU9szuGJwrnk0K+fzfoayeA7Vat26tnTSiXOxYknCueSwTHHQJs28MAD4S9m9Pnntu62N2K7OPNE4VwyqFIF/vQn+N//4MMPwzvmuedsdPfAgfGNzVV4niicSxa/+x00bgx33VX6vrm58OST8JvfQL168Y/NVWieKJxLFmlpcO218N57MHduyftef701Yt9xR2JicxWaJwrnksnIkVCnDoweXfw+06fDyy/DjTdC69YJC81VXJ4onEsmtWvDH/5gEwUuXfrr13fvhiuvtIbvP/0p8fG5CskThXPJ5uqroUYN6NsXFi7c97V//QsWLYL774fq1YOJz1U4niicSzYHHQTTptm4iqOPtgn/pk+HE0+0UsTpp8dmRTznwuRTTTqXjI44AmbPhv79YcAA+16TJnDvvdaOIRJsfK5C8UThXLJq1szGVPzjH7a86bBhXt3kAhFIohCRe4GzgF3ACuB3qvp9EfvlANuAPGCPqpZhCTDnUlBaGtx5Z9BRuAouqDaKqUBnVe0KfAXcWMK+x6tqhicJ55wLRiCJQlXfU9U9oaf/A1oEEYdzzrnSiYY7AVm8AhB5G3hZVf9TxGurgO8ABR5X1TElnGc4MBwgrUFaZqc7O0UVT/ambDo07BDVsanKr7n8q2jXC37NkZp7+dy5xdbcqGpcNuB94MsitrML7XMT8AahhFXEOZqFHg8CFgC9w3nvzMxMjVbm49Efm6r8msu/ina9qn7NkQLmaDGfqXFrzFbVk0p6XUQuBs4ETgwFWdQ51oYeN4jIG8CRwMxYx+qcc654gbRRiEhf4Aagn6r+VMw+NUSkVsHXwClYicQ551wCBdXr6WGgFjBVROaLyGMAItJMRAqW+GoMfCQiC4DPgEmq+m4w4TrnXMUVyDgKVT2smO+vBU4Pfb0SSE9kXM45534t8F5P8SAiG4Gvozy8IbAphuGkAr/m8q+iXS/4NUfqYFVtVNQL5TJRlIWIzNEKNrjPr7n8q2jXC37NseSzxzrnnCuRJwrnnHMl8kTxa8WO/i7H/JrLv4p2veDXHDPeRuGcc65EXqJwzjlXIk8UzjnnSuSJoggicq+ILBGRhSLyhojUDTqmeBKRc0VkkYjki0i57k4oIn1FZKmILBeRUUHHE28i8rSIbBCRCjP9jYi0FJHpIpId+ru+OuiY4k1E0kTkMxFZELrmv8fy/J4oihbJwkrlwZfAAMr5hIsiUhl4BDgN6AicLyIdg40q7p4F+gYdRILtAa5T1Q5AT+APFeD3/DNwgqqmAxlAXxHpGauTe6IoglawhZVUNVtVlwYdRwIcCSxX1ZWqugsYB5wdcExxpaozgS1Bx5FIqrpOVeeFvt4GZAPNg40qvkIzhW8PPa0a2mLWU8kTRemGAe8EHYSLiebA6kLPcynnHyAVnYi0BroBnwYbSfyJSGURmQ9sAKaqasyuOZBJAZOBiLwPNCnipZtU9a3QPjdhxdgXExlbPIRzvRWAFPE97x9eTolITeA14BpV/SHoeOJNVfOAjFCb6hsi0llVY9I2VWETRSwWVkolpV1vBZELtCz0vAWwNqBYXByJSFUsSbyoqq8HHU8iqer3IjIDa5uKSaLwqqcihLOwkktJs4G2ItJGRA4ABgETAo7JxZiIFaLXNAAAAWhJREFUCPAUkK2q/xd0PIkgIo0KemeKSHXgJGBJrM7viaJoRS6sVF6JSH8RyQWOAiaJyJSgY4qHUAeFK4ApWAPnK6q6KNio4ktExgL/BdqLSK6IXBJ0TAlwDHARcELo/3e+iJwedFBx1hSYLiILsRuiqao6MVYn9yk8nHPOlchLFM4550rkicI551yJPFE455wrkScK55xzJfJE4ZxzrkSeKJxzzpXIE4VzzrkSeaJwLs5E5IjQ2iZpIlIjtF5A56Djci5cPuDOuQQQkduBNKA6kKuqdwUcknNh80ThXAKE5paaDewEjg7N9OlcSvCqJ+cSoz5QE5tDLC3gWJyLiJconEsAEZmArajXBmiqqlcEHJJzYauw61E4lygiMgTYo6ovhdbt/kRETlDVaUHH5lw4vEThnHOuRN5G4ZxzrkSeKJxzzpXIE4VzzrkSeaJwzjlXIk8UzjnnSuSJwjnnXIk8UTjnnCvR/wMDurMYZkmjXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#\n",
    "def func(x):\n",
    "    f = A*(np.tanh((x - x0)/L + w/(2*L)) - np.tanh((x-x0)/R - w/(2*R)))\n",
    "    return f\n",
    "#\n",
    "def dfunc(x):\n",
    "    f = A*(-(-np.tanh(-w/(2*R) + (x - x0)/R)**2 + 1)/R + (-np.tanh(w/(2*L) + (x - x0)/L)**2 + 1)/L)\n",
    "    return f\n",
    "#\n",
    "A = 2.0\n",
    "x0 = 0.5\n",
    "L = 1.1\n",
    "R = 0.5\n",
    "w = 0.4\n",
    "#\n",
    "xArr = np.linspace(-2, 3, 100)\n",
    "yArr = func(xArr)\n",
    "dyArr = dfunc(xArr)\n",
    "#\n",
    "plt.figure(figsize = (6, 4))\n",
    "plt.title(\"Function and its derivatives\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.plot(xArr, yArr, color = 'b', label = \"function\")\n",
    "plt.plot(xArr, dyArr, color = 'r', label = \"first deriv\")\n",
    "plt.grid(color = 'g')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration\n",
    "\n",
    "The syntax for performing indefinite integration (i.e. without limits) is similar to that for differentiation. Note that sympy does not add a constant of integration to indefinite integrals, you have to think about what you need to do about these yourself!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral of a*x**2 + b*x + c w.r.t. x is a*x**3/3 + b*x**2/2 + c*x\n"
     ]
    }
   ],
   "source": [
    "intQuad = sp.integrate(quad(x), x)\n",
    "print(\"Integral of\",quad(x),\"w.r.t. x is\",intQuad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to evaluate this integral with particular values of `a`, `b`, `c` and `x` we must `sp.subs` all of these, including `x`, as `intQuad` is not a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral of a*x**2 + b*x + c with a = 1, b = 2, c = -1 and x = 4 is 100/3\n"
     ]
    }
   ],
   "source": [
    "print(\"Integral of\",quad(x),\"with a = 1, b = 2, c = -1 and x = 4 is\",intQuad.subs(a, 1).subs(b, 2).subs(c, -1).subs(x, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to do a definite integral, give `sp.integrate` a tuple containing the variable and the range over which the integration should take place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral of a*x**2 + b*x + c w.r.t. x over range -1 to 2 is 3*a + 3*b/2 + 3*c\n"
     ]
    }
   ],
   "source": [
    "intQuadLim = sp.integrate(quad(x), (x, -1, 2))\n",
    "print(\"Integral of\",quad(x),\"w.r.t. x over range -1 to 2 is\",intQuadLim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerical values of `a`, `b`, and `c` can be provided as previously described. The substitution can be made before or after the integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral of a*x**2 + b*x + c w.r.t. x over range -2 to 2 is 3\n",
      "Integral of a*x**2 + b*x + c w.r.t. x over range -2 to 2 is 3\n"
     ]
    }
   ],
   "source": [
    "intQuadLim1 = sp.integrate(quad(x).subs(a, 1).subs(b, 2).subs(c, -1), (x, -1, 2))\n",
    "intQuadLim2 = sp.integrate(quad(x), (x, -1, 2)).subs(a, 1).subs(b, 2).subs(c, -1)\n",
    "print(\"Integral of\",quad(x),\"w.r.t. x over range -2 to 2 is\",intQuadLim1)\n",
    "print(\"Integral of\",quad(x),\"w.r.t. x over range -2 to 2 is\",intQuadLim2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting with sympy\n",
    "\n",
    "We can demonstrate the plotting capabilities of sympy by plotting `quad` and its integral. Sympy does this via an interface to matplotlib, but the syntax for making the plots is a little different to that we have seen so far. Here are the two plots as two  separate figures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWcAAAD3CAYAAADBqZV6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3zOdf/A8dfXZk45ZgubQ85rwzDFXbdSDh0UqTv8uEvjpoO6b0QqKmUhpKTSbodwKyFqOUV0pGgOZUiKYQ5ZsjCHMZ/fH28kNru26/D9Xtf1fj4ee5jtur7ft23e12efz/vz/ljGGJRSSjlLEbsDUEopdSlNzkop5UCanJVSyoE0OSullANpclZKKQfS5KyUUg6kyVkFHcuyJlqWNdTuOJS6HEvrnJWvWZaVBlwF5Fzw4brGmL1euFcPoJcx5gZPX1spbwq1OwAVtO40xnxqdxBKOZVOayhHsCzrJsuy0i/6WJplWa3Pvv+8ZVmzLcuablnWEcuyNlmWFX/BY6taljXPsqwMy7IOWpY1wbKsaGAi0MKyrKOWZWWefew7lmUNv+C5/7Is62fLsn63LCvZsqwqF3zOWJb1kGVZ2yzLOmRZ1huWZVne/noopclZ+ZO7gFlAOSAZmABgWVYIsADYCdQAIoFZxpgtwEPAN8aYK4wx5S6+oGVZNwMjgPuAymevMeuih7UHmgGNzj6unaf/YUpdTJOzssuHlmVlnn370MXnfG2MWWSMyQFmIMkS4FqgCjDQGJNljDlhjPnaxWt2A6YYY9YZY04CTyEj7RoXPGakMSbTGLML+AyIc/HaShWaJmdll47GmHJn3zq6+Jz9F7x/DChuWVYoUBXYaYw5XYg4qiCjZQCMMUeBg8joO6/7XlGI+yhVIJqclVNkASXP/eXsVEW4i8/dDVQ7m6gvll850l6g+gX3LQVcCexx8d5KeYUmZ+UUPyEj4TssyyoKDAGKufjcNcA+YKRlWaUsyypuWdb1Zz/3KxBlWVZYHs99F3jQsqw4y7KKAS8Bq40xaYX+lyjlAZqclSMYY/4AHgEmIaPWLCD9sk/687k5wJ1AbWDX2ed1PvvpFcAmYL9lWb/l8tzlwFDgAyTB1wK6uPNvUcoTdBOKCkqWZUUA1yNzzseBVCDFGHPG1sCUOkuTswoqlmW1AgYDFYD1wAGgOFAXGTXPBcYaYw7bFqRSaHJWQcayrNHA62fL4i7+XChS0xxijPnA58EpdQFNzkop5UC6IKiCkmVZMyzLKnvB32tYlrXczpiUupC7jY902K380sSJExk3blz3RYsWsWfPHurUqcPYsWNBf6aV97nUm8XdaQ39QVZ+6+uvv6ZVq1ZUrFiR9evXU6lSJbtDUsHBpeSs0xoqKM2YMYOEhASmT59Ojx49uP322/n+++/tDkup83TkrIJSx44dSUpKIiIiAoA1a9bQp08f1q9fb3NkKgjotIZSBZGdnU1YWF67vJXyGJ3WUOpiw4cP5/fff8/1c2FhYaxYsYIFCxb4OCqlLqXHVKmg0qBBA+68806KFy9OkyZNCA8P58SJE2zbto0NGzbQunVrnn76abvDVEpHziq4zJ07l5UrV9KuXTtiYmLIycmhTJkydO/enTVr1jBu3DjCwy/fqTQhIYGIiAhiY2Mv+dyYMWOwLIvffrukx5JSBaIjZxVU1q5dy86dO5k5cyafffbZXz53/PhxSpQoke81evToQd++fbn//vv/8vHdu3ezbNkyqlWr5tGYVXByKzkbA3rUpfInDz30ELfeeivbt28nPv78+bAYY7Asi+3bt+d7jZYtW5KWlnbJx/v168fLL79Mhw4dPBmyCiBnzkARF+cr3JrWeOEFd56tlO89/vjjbNmyhYSEBLZv337+bceOHS4l5rwkJycTGRlJo0aN8n1sUlIS8fHxxMfHExMTU+h7Kv8zcKDrj3WrlK5SJczOnVDM1fMqlAoQaWlptG/fntTUVI4dO0arVq1YunQpZcuWpUaNGqSkpFCxYsV8rxMfH09KSooPIlZ2O3oUIiPhjz98UEr3668wZ447V1DK//3yyy/s2LGDRo0aUaNGDdLT02nSpAn79+/P/8kqaEybBocL0CXcrTnnevXg9dehe3d3rqKUf2vQoAEHDhw4//eCjJxVcDhzRnJls2auP8etkfNjj8GaNbB6tTtXUcq/dO3alRYtWrB161aioqKYPHmy3SEph1u2DLZuhX//2/XnuDXnfOQIJioK2reHmTMLfRmlgpbOOQeH22+H9eth504IC/PBnHPp0pCQALNnw9697lxJKaUC008/weLF8PDDUJDWLW7vEHz0UcjJgbffdvdKSikVeCZMgKJFoU+fgj3P7eRcuzbccYck55Mn3b2aUkoFjj/+gKlToUsXuOqqgj3XI701Hn9cy+qUUupi77wj9c0FWQg8xyPJuXVraNcOFiyQLd1KKRXscnKkfO5vf4OmTQv+fI8kZ8uCjh3h/ffh6689cUWllPJvixfDL7/IzEJheKxl6P33Q4UKMG6cp66olFL+a/x42a7dqVPhnu+x5FyypKxGfvghuNE/Riml/N7mzTJqfuQRqdQoDI8223/0UQgJkVcMpZQKVq++Kns//vWvwl/Do8k5MhI6d4bJk6WERCmlgk1GBkyfDg88APkcqnNZHj+mql8/KR3RdgNKqWD05puy5+M//3HvOh5Pzk2bQsuWMrVx+rSnr66UUs514gS88YZszKtf371reeWA1379pMHHhx964+pKKeVMM2fKtEb//u5fy62udECuT87Jgbp1oVIlWLnSncsrFdi0K13gMAZiY6U6Y/36y56v6v2udHkJCZHtiqtWaa9npVRwWLpUSuj69/fMwddeSc4ADz4IZcrAjBneuoNSSjnHK69A5crS5MgTvJacS5eW0fPEiZDLKfJK+a2EhAQiIiKIjY09/7GBAwdSv359GjZsyN13301mZqaNESpfS02VkXPfvgXr2Xw5XkvOIAXYliUF2UoFih49erBkyZK/fKxNmzakpqbyww8/ULduXUaMGGFTdMoO48bJLumHHvLcNb2anKtWha5dYdIkOHTIm3dSyndatmxJhQoV/vKxtm3bEhoq5yU3b96c9PR0O0JTNti/X6o0evSQ/kKe4tXkDDBgAGRlyfSGUsFgypQp3HbbbXl+Pikpifj4eOLj48nIyPBhZMobJkyQgai7m04u5vXk3KgRtG0rm1L0pBQV6BITEwkNDaVbt255PqZ3796kpKSQkpJCuDv7e5XtjhyRTScNG0KdOp69tteTM8DAgX8O/ZUKVNOmTWPBggXMnDkTyxO1VMrxJk2CzEwYNMjz1/ZJcr7lFoiLgzFj4MwZX9xRKd9asmQJo0aNIjk5mZIlS9odjvKBU6ekfO7GG+G66zx/fZ8kZ8uCJ56ALVtg0SJf3FEp7+natSstWrRg69atREVFMXnyZPr27cuRI0do06YNcXFxPOTJZXvlSO+9B+np3hk1g5e2b+fm1CmoVQtq1oTPP3fnlkoFDt2+7Z+MgQYNZOD5ww8F3hFo3/bt3BQtKquZX3yhW7qVUv5t8WLYtElGzd5aXvDZyBlkZbNjR6kFnDPHndsqFRh05OyfbrwRduyQo6gKcQyVs0bOIFu6mzeHDz6ArVt9eWellPKMb7+FL7+U1siFPR/QFT5NziD9NooVg9GjfX1npZRy3+jRUL68e+cDusLnyTkiAnr2lDO2dIerUsqfbN0K8+fLqdpXXOHde/k8OYOU1Z05I81ClFLKX4wdK13nHnvM+/eyJTnXqCENkd5+G37/3Y4IlFKqYPbvh2nTpFf9VVd5/362JGeQEpSsLGkaopRSTjdtmrQFHTDAN/ezLTk3aADt20tDpKwsu6JQSqn8HToEiYlw++1Qu7Zv7mlbcgYYPBgOHpTmIUop5VQTJsg+DW9t1c6NTzeh5KZlSznG6uefPXe8i1L+QjehON/Ro1C9Ovztb/Dxxx65pPM2oeRm8GD5h7//vt2RKKXUpc4VLjzzjG/va/vI2Rho0kTmnbdsgZAQd6+olP/QkbOznTgBV18N11wDy5d77LL+MXK2LBgyBLZtg9mz7Y5GKaX+NHWqlND5etQMDhg5g2xIiY2FIkWk/V4R218ylPINHTk716lTcvRU5cqwapVHu8/5x8gZJBk/84y04PvoI7ujUUopePdd2LlTcpMdp445IjkDdO4s9YPDh8s8tFJOlZCQQEREBLGxsec/9vvvv9OmTRvq1KlDmzZtOHTokI0RKnfl5MCIEXJA9R132BODY5JzaCg89RSsWwdLltgdjVJ569GjB0su+iEdOXIkt9xyC9u2beOWW25h5MiRNkWnPGHePGly9PTT9oyawSFzzudkZ8scT2QkrFxp3xdFqfykpaXRvn17UlNTAahXrx6ff/45lStXZt++fdx0001sdaFpuc45O48xcN99sHGjTLV6oYLMf+aczwkLgyefhG++gc8+szsapVz366+/UrlyZQAqV67MgQMH8nxsUlIS8fHxxMfHk5GR4asQlYuSk2HuXHj2WXtLex01cgapK6xZE6KjPVpXqJRHXTxyLleuHJmZmec/X758eZfmnXXk7CzGQNOmslV7yxaZbvUC/xs5AxQvDsOGwY8/wldf2R2NUq656qqr2LdvHwD79u0jIiLC5ohUYXz8MaxfLxUaXkrMLnNccgbo1k1WS4cNszsSpVxz1113MW3aNACmTZtGhw4dbI5IFZQxknNq1YLu3e2OxqHJuWRJmXtevlxHz8p5unbtSosWLdi6dStRUVFMnjyZwYMHs2zZMurUqcOyZcsYPHiw3WGqAlqwQKrFnDBqBgfOOZ9z7JjMPcfGwqefeusuStlL55ydwRho1kwaHG3d6t1TtfHXOedzdPSslPKVRYtg7Vrp8+PlxOwyx46cAY4fl9GzhztCKeUYOnK2nzFw7bVy8IcPRs3g7yNngBIlZPS8YgV8+aXd0SilAtHixZCSInPNThk1g8NHzqCjZxXYdORsL2Pg7rvh5EnZfOKj5Oz/I2eQ0fOgQTJRr6NnpZQnJSdLJ8x//MNZo2bwg5EzyOi5Vi3pWvfFF9pzQwUOHTnb58wZiIuT/OLF3YC5CYyRM8joecgQqdpYtszuaJRSgWDuXGluNGyYM+qaL+YXI2eQjnV160J4OKxZo6NnFRh05GyP06dlD0VIiJy+5OMGR4EzcgbpWPfcc7KqmpxsdzRKKX/27rtSNvfCC849VNpvRs4gr3bXXCPNkTZs0LMGlf/TkbPvnToF9epBuXIy2LMhjwTWyBlkXmjYMJkn0pO6lVKFMXUq7NgBL77o7AGeX42cQVZYGzWSOehNm5w5ka+Uq3Tk7FsnTshpS1FRHj9RuyACb+QM8kr34ovw008wY4bd0Sil/Mn06ZCeLgdJO72owO9GzvDnXviyZWHhQihWzI4olHKfjpx95+hR2S/Rpo0M7GxMzoE5cgb5oo4aJdu5337b7miUUv7glVfgwAF47DHnj5rBT5MzQKtWcPPN8uvJkSN2R6OUGDduHDExMcTGxtK1a1dOnDhhd0gKyMiAMWOgUye47jq7o3GN3yZny4KRI+WLPnas3dEoBXv27GH8+PGkpKSQmppKTk4Os2bNsjssBbz0EmRlQWKi3ZG4zm+TM8jJBffeK8n5MifRK+Uzp0+f5vjx45w+fZpjx45RpUoVu0MKejt3wptvwoMPQv36dkfjOr9OziDTGseP+9crogpMkZGRPPHEE1SrVo3KlStTtmxZ2rZte8njkpKSiI+PJz4+noyMDBsiDS7PPSdVXs8/b3ckBeP3yblePUhIgLfeksJypexy6NAhPvroI3bs2MHevXvJysrif//73yWP6927NykpKaSkpBAeHm5DpMEjNVXK5x57TGqb/YnfJ2eQV8aQEPlTKbt8+umnXH311YSHh1O0aFE6derEqlWr7A4rqD39NJQpA/54GHpAJOfISBg4UOaWNmywOxoVrKpVq8a3337LsWPHMMawfPlyoqOj7Q4raK1aBQsWyFF3FSrYHU3B+eUmlNxkZkqBedOmsHSp3dGoYPXcc8/x/vvvExoaSuPGjZk0aRLFLrNLSjeheIcx0KKFnG7yySdQsqTdEf2FS1XWAZOcAV59Ffr1gyVLoF07u6NRKn+anL1j9mzo3BkmT5Y1KYcJvuScnQ3R0fIquWGDc/u0KnWOJmfPO3lS8kDp0rBunSPzQOBu385LWBiMGCErtNOm2R2NUsoOEyZI5daYMY5MzC4LqJEz/DnXtHu3dK4rVcruiJTKm46cPevgQTkIunlzWLzY7mjyFHwjZ5Bt3WPGwN69MG6c3dEopXzpxRfh8GEYPdruSNwXcMkZ4IYboGNH+O9/Yf9+u6NRSvnCtm3wxhvQs6cc3urvAjI5A7z8svRvffZZuyNRSvnC4MHS2/2FF+yOxDMCNjnXqQMPPACTJunGFKUC3ddfw7x5suGkUiW7o/GMgFsQvNChQ5KkGzSAFSv8o8G2Ci66IOg+Y2QBMD3db4oAgnNB8ELly8sCweefw0cf2R2NUsob3nsPiheXns1+kJhdFtAjZ4DTpyEuTk7d3bRJzxtUzqIjZ/ccPSqdKatUgdWrpTWoH9CRM0BoqJTU/fILjB9vdzRKKU8aMULKZseP95vE7LIA++fkrk0baN9epjh+/dXuaJRSnvDLL7KnoXt32XgWaIIiOYN8E48f19I6pQLFgAHSdW7UKLsj8Y6gSc716knj7TVrpBmKUsp/LV0qi/xDhsh8cyAK+AXBC2VmSpKuVUvqIgNtjkr5H10QLLhTp6BRI+lC6aeL/LogeLFy5eRXoG++gVyOdlNK+YE33oAtW2Sh3w8Ts8uCKjkD3H+/FKwPGgR//GF3NCrQZGZmcu+991K/fn2io6P55ptv7A4poBw4IKdot2sni/yBLOiSc5Ei8Prr8k0OlD34yjn+/e9/c+utt/Ljjz/y/fff6xmCHvbMM5CVJaceBfqO36BLzgDx8dCrF7z2msxZKeUJhw8f5ssvv6Rnz54AhIWFUa5cOZujChxr18oi4OOPQ/36dkfjfUGZnEG2epYpI99o99ZElRLbt28nPDycBx98kMaNG9OrVy+ysrIueVxSUhLx8fHEx8eTkZFhQ6T+JycHHnpI1o2CpRw2aJNzxYqyKSUrC+bOtTsaFQhOnz7NunXrePjhh1m/fj2lSpVi5MiRlzyud+/epKSkkJKSQnh4uA2R+p+334aUFBg2DMqWtTsa3wja5AzQp4/03nj8cV0cVO6LiooiKiqK6667DoB7772XdVpU77b9++Gpp6B1a+jSxe5ofCeok3NoKEycKFu6hwyxOxrl7ypVqkTVqlXZunUrAMuXL+eaa66xOSr/N2CANC57443AXwS8UKjdAdgtPh4efVS+8Q88IH9XqrBef/11unXrRnZ2NjVr1mTq1Kl2h+TXPv0U3n0XnnsO6ta1OxrfCqodgnn54w9Z/a1SRbZ3+/Nx6sq/6A7BvJ04AQ0byoL9xo3SszlA6A5BV5UtK3WT69bBm2/aHY1SCuQc0HOHtgZQYnaZjpzPMgZuvVW2dv/4Y+A2U1HOoiPn3G3bJsfLdewIs2bZHY3H6ci5ICxLXqGzs6FfP7ujUSp4GSPrQMWKwSuv2B2NfTQ5X6B2bdkeOnu2tCRUSvne7NmwbBkkJgb3b7A6rXGRkyeha1eZf05NhSuusDsiFch0WuOvMjOhWzfpffPttwG7OK/TGoVRrBj07w87d8LQoXZHo1RwGTAAPvkEkpICNjG7TJNzLm64AR55RBojrV5tdzRKBYdly2DKFBg4EBo3tjsa++m0Rh4OH4aYGGm0snYthIXZHZEKRDqtIY4ehdhYKZnbsCHgS+d0WsMdZcrAW2/JvHMuvWuUUh701FOwaxdMnhzwidllmpwvo317abQyfDhs3mx3NEoFpq+/hgkToG9fuP56u6NxDp3WyMeBA9CkCdx8M0ydqosUyrOCfVrj+HGIi5P9BRs3Bk11lE5reEJEBIweDTNmyIGSSinPGTYMfvpJqjOCJDG7TJOzC7p0gQ4dpK3ojz/aHY1SgWHtWhgzBhISoE0bu6NxHk3OLrAs6ftcqhQ8+KAcmaOUKrzsbEnKEREwdqzd0TiTJmcXVaokixbffhvc+/2V8oTERChfXiqi9Azc3GlyLoAuXeDuu2Xn4JYtdkejlH9avVqSc/XqMl2ocqfJuQAsS17po6JkIePUKbsjUk6Uk5ND48aNad++vd2hOE5WFvzznxAZCePH2x2Ns2lyLqCrrpJNKe+/L6/+Sl3stddeIzo62u4wHGngQPj5Z5g2LXhO0S4sTc6FcO+98uo/fLj23lB/lZ6ezsKFC+nVq5fdoTjOkiXym2e/fnDTTXZH43yanAvp9dflV7N//lN+VVMK4D//+Q8vv/wyRYrk/V8rKSmJ+Ph44uPjycjI8GF09jl4UKozYmL0N05XaXIupLJl5Vezn3+GJ56wOxrlBAsWLCAiIoKmTZte9nG9e/cmJSWFlJQUwsPDfRSdfYyBhx+G336TzVzaO8M1mpzdcNNN0n924kRYtMjuaJTdVq5cSXJyMjVq1KBLly6sWLGC7t272x2W7d59F+bMkUV0bQXqOu2t4aaTJ6FZMyhRAhYsgCAYCCkXfP7554wZM4YFCxZc9nGB3ltj1y5o2FCmM778UnvTnKW9NXyhWDEZGfz0E/ToIb/CKaVkJ22fPtI4bPp0TcwFpcnZA2JjpXJj0SI5PUWpm266Kd9Rc6BLTJQKjQcegFq17I7G/2hy9pBHHoG77oJBg+RwWKWC2VdfyRxz9+5w//12R+OfdM7Zgw4ehEaNpEHS2rXaAlHlLxDnnA8elB7NxYvLQKV0absjchydc/a1K6+EmTNh2zZ4/HG7o1HK94yBnj3h119h1ixNzO7Q5OxhN94ofZ/few9mz7Y7GqV864034KOPYNQoyKfcW+VDpzW84PRp6WC3ZAl89x1omwWVl0Ca1tiwAZo3h9at4eOPpVGYypVOa9glNFSqNkqWlD4cur1bBbrMTPlZv+UWOWtTE7P7NDl7SWSkTG1s2SK1nlr/rAKVMVLjv3MnPP20bsTyFE3OXnTLLVJONHOmHGCpVCAaPVrmmceMgeuvtzuawKFzzl525gzcfjt89hmsWqWLJOqv/H3O+fPPZRBy771SnaHTGS5x6aukydkHfvtNGr4ULQopKVChgt0RKafw5+S8d6/8XFeoAGvWaNlcAeiCoFNUrChdua6+WnZM6endyt+dOgX33SeL3R98oInZGzQ5+0jz5vLDvHix1EEr5c9eeglWroT//heuucbuaAJTqN0BBJM+fWQ768iRsr21c2e7I1Kq4KZMgeefhxdfhK5d7Y4mcOmcs49lZ8PNN0uS/uYb6cWhgpe/zTmvXAmtWslO2MWLpaZfFZjOOTtRWBjMnSuLKB06yGKhUv5g1y7o1AmqV5fT5zUxe5cmZxtUqgQffgj798s89KlTdkekPGH37t20atWK6OhoYmJieC2AmntnZUlL3BMnZGu2Vhx5nyZnm8THy2LKgQPw5JO6gzAQhIaGMnbsWLZs2cK3337LG2+8webNm+0Oy21nzkjD/I0bpZa5fn27IwoOmpxt9M9/Qvv2MG6cnqASCCpXrkyTJk0AKF26NNHR0ezZs8fmqNz3wgtSLvfyy3DbbXZHEzx0QdBmZ87AP/4B8+fLW4cOdkekPCEtLY2WLVuSmppKmTJl/vK5pKQkks7u58/IyGDnzp12hOiSefPgnnukd8aUKboD0EN0h6C/OHYMbroJNm2SE4p1i7d/O3r0KDfeeCPPPPMMnTp1uuxjnVytsWoVtGsnyfntt+UwY+URWq3hL0qWhORk2Ul4552we7fdEanCOnXqFPfccw/dunXLNzE72dat8rNYubI0NNLE7HuanB2iUiVYuFBWxdu3h8OH7Y5IFZQxhp49exIdHU3//v3tDqfQ9u+HW2+VUrnFi2XQoHxPk7ODxMZKD45y5WT3YHa23RGpgli5ciUzZsxgxYoVxMXFERcXx6JFi+wOq0COHIE77pAqogULoFYtuyMKXlpG7jBt20q3rwcflEWY//0PiuhLqF+44YYbcHMNx1bnmhl9/71MszVrZndEwU2TswP16CGnFw8eLKdKvPqqrpIr7zJGer8sWSL197ffbndESpOzQw0aJAl63DiZj37qKbsjUoFs9Gg5++/ZZ6FXL7ujUaDJ2bEsS1bJDxyQc9kiIqBnT7ujUoFoxAj5GRsyRLrNKWfQ5OxgRYpI4X9mpvx5xRXaZlR51rhxkpj/7/8kMev0mXPoUpPDhYVJB7CQEOjWTXYRKuUJb74J/fvLJpNp0+RnTDmHJmc/UKqU1EA3ayYjZz+rzlIONGUKPPqobDR5911t/+lEmpz9ROnSsiGgYUPpqfvpp3ZHpPzVzJmy6NeundTVh4XZHZHKjSZnP1KuHHzyCdSrJ711v/zS7oiUv5kzB+6/X3q5zJun27KdTJOzn7nySli2DGrUgIEDNUEr1yUny8JfixbyfsmSdkekLkeTsx+KiIDly6Xd6K23wtKldkeknG7RIknMTZrI+1dcYXdEKj+anP1U5cryn6xuXVnUSU62OyLlVO+9J33C77pLdgBe1F5aOZQmZz8WHg6ffQaNG8si4fvv2x2Rcpq33pISzL/9DSZOhPLl7Y5IuUqTs58rX17moK+/Xn5tnTrV7oiUExgjO/8eeUS6zOmI2f9ocg4A58rsWreW5jVvvWV3RMpOxsihwU8/LaPmefOgRAm7o1IFpck5QJw7TeXhh2W01L+/LBiq4JKTA//6lzQyevRRmD4diha1OypVGLovKIAUKwavvCLvjxsHO3dKP2gdNQWH48flhfmdd2DoUBg2THtl+DMdOQeYkBB47TVJzvPnwy23wG+/2R1V8FiyZAn16tWjdu3ajBw50mf33b8fWrWSuvfx4+GFFzQx+ztNzgHqP/+R3WDr18umg59/tjuiwJeTk8Ojjz7K4sWL2bx5M++99x6bN2/2+n2//x6uvRY2bpTpjMce8/otlQ9ocg5g99wjm1UOHZIEvWqV3REFtjVr1lC7dm1q1qxJWKe3/YcAAAuYSURBVFgYXbp04aOPPvLqPZOTpVLnzBn46ispqVSBwXLnzLOYmBhTwuETmhkZGYSHh9sdRr68GefJk7B9O5w+LZtX3DlNWb+eeTt06BCHDx+mevXqABw8eJCsrCyqVat2SWy/nZ1rOnnyJHFxcYW63/79sGePLAbXru3dhT/9vnvO2rVrNxljYvN9oDGm0G9NmzY1TucPMRrj/TgzMoxp08YYMKZXL2OOHy/cdfTrmbfZs2ebnj17nv/79OnTTd++fS/7nJIlSxb4PsePG5OQIN/Lf/zDmKysAl+iwPT77jlAinEhv+q0RpCoWFFqoZ9+GiZNgr//HXbtsjuqwBIVFcXu3bvP/z09PZ0qVap49B47d8KNN8LmzVKRMWuWNjAKVJqcg0hICCQmShXH1q3SBEf7QntOs2bN2LZtGzt27CA7O5tZs2Zx1113eez6H38sW/W3bIEnnpCKjCL6PzhgufWt7d27t6fi8Bp/iBF8G2fHjvDdd3DVVdJwfcQI1zes6Nczb6GhoUyYMIF27doRHR3NfffdR0xMzGWfU9GFBYBTp+Q09rvugurVYd06Wez1Jf2+e1SSKw9ya0EQcOvJyl5Hj0JCgiwslSolfTkqVbI7quASHx9PSkpKnp9PT4cuXWDlSnjoIalfL17chwEqb3CpAl1/KQpiV1whney6dYPPP4cGDbT1qJPMmQM9e0od87vvSs8UTczBQ5NzkLMsaZa0bh1UrSp9f/v0gawsuyMLXgcPQteucN99cOQIrF0rf1fBxa3k/PzzzxMZGUlcXBxxcXEscvix0GPGjMGyrPM1pk4zdOhQGjZsSFxcHG3btmXv3r0+u3d0NHz7rcxt/ve/0L173kdgDRw4kPr169OwYUPuvvtuMjMzfRZnQcyZM4eYmBiKFCly2akDO5zb5p2amvqXbd4ffwyxsfDBB/Dii/I9qFvXvjgTEhKIiIggNjb/sly77N69m1atWhEdHU1MTAyvvfaa3SHl6sSJE1x77bVYlvW9ZVmbLMsadtknuFJvl9fbc889Z0aPHu2z+kB37Nq1y7Rt29ZUq1bNZGRk2B1Orv7444/z77/22mumT58+tsTx1VfG1KwpdbQPPWTMBWEZY4z55JNPzKlTp4wxxgwaNMgMGjTIhijzt3nzZvPjjz+aG2+80Xz33Xd2h3Pe6dOnTc2aNc0vv/ximjRpYho2bGhWrdpiHnhAvuYNGxqzYYPdUYovvvjCrF271sTExNgdSp727t1r1q5da4wx5vDhw6ZOnTpm06ZNNkd1qTNnzpgjR44YI+t8RYHVQHMT7HXO/fr14+WXX8ZycDeYMhd0Q8/KyrIt1htugB9+kLajSUkQEwMLF0qfYIC2bdsSGioNDZs3b056erotceYnOjqaevXq2R3GJS7c5m1ZFvXqDadt20jWrYMhQ6SSplEju6MULVu2pEKFCnaHcVmVK1emSZMmAJQuXZro6Gj27Nljc1SXsiyLK/48vLHo2bc8iyrcTs4TJkygYcOGJCQkcOjQIXcv5xXJyclERkbSyCk/8ZfxzDPPULVqVWbOnMkLL7xgWxylSsHYsfDNN1Crlsx53nEH/PTTXx83ZcoUbrvtNnuC9FN79uyhatWqbNwo9eZz5txJqVK/MW2aTGWEhdkdof9KS0tj/fr1XHfddXaHkqucnBwsy9oAHACWGWNW5/XYfPs5W5b1KZBbgdUz+/fvZ+jQoViWxdChQxkwYABTpkwpfORuaN26Nfv377/k44mJibz00kssdcgR1ZeLs0OHDiQmJpKYmMiIESOYMGECw4ZdflrKWy6M05hQSpbsyiefPEJMTCn69y/CkCEwfnwioaGhdOvWzZYYL47zQue+nk6UlRXC6tX30rix/D0h4RtKlHiPxo3H2xuYnzt69Cj33HMPr7766l9+C3WSkJAQjDFxlmWVA+ZblhVrjEnN9cF5zXe4+Hbejh07HDkv9cMPP5jw8HBTvXp1U716dRMSEmKqVq1q9u3bZ3dol5WWlua4r+e+fcb06CHzotHRe8zVV480mZk+aOzgJqfMOR87Zszo0caUKZNtypdfb3r3NqZRo6bmpZdeMi+99JLd4eXJqf+3L5SdnW3atm1rxo4da3corsDIHOFzwBPGG3PO+/btO//+/PnzHbmi26BBAw4cOEBaWhppaWlERUWxbt06Kjlwt8W2bdvOv5+cnEz9+vVtjOZSlSrJRpVx474hPT2VHTuepEmTksyYIccjqdydPg2TJ0vVxcCB0KJFCMWLP8PgwTsICTEe3+YdbIwx9OzZk+joaPr37293OHnKyMg4X9lkWVYJoDXwY55PyCtru/LWvXt3Exsbaxo0aGDuvPNOs3fvXh+++BRO9erVHVut0alTJxMTE2MaNGhg2rdvb9LT0+0OKVe1atUykZFR5uqr+5rixbcYMOaaa4yZO9eYnBy7o/vTvHnzTGRkpAkLCzMRERGmbdu2Pr3/qVPGzJplTP368tvGddcZ89ln8rmFCxeaOnXqmLCwMDN8+HCfxlUQXbp0MZUqVTKhoaEmMjLSTJo0ye6QLvHVV18ZwDRo0MA0atTINGrUyCxcuNDusC7x/fffm7i4OAP8AKQCz5rL5Ffdvq3ccuYMzJ0Lzz4LFSpAZiYMGCC7DoN1N9uRIzJSfvVV6SJ3221y6GrHjpceHZXf9m0VkHT7tvK+IkVkJ1tqKvTrJ4fM9uolDXqGD5fdbsFizx4YPFh2WvbrJ39++CEsWAB3361n+qmC0ZGz8ihjYMUKKcNbvBhuvlmSVO/eclRWoCUoY+CLL6T3xTvvyNz7PffIbw+uVHPpyDkoufS/QJOz8ppNm2DaNGnYc/So/Hr/97/LlMdFJzf5nW3b4L335IXoiy+gbFk5VPf++6FmTdevo8k5KGlyVs5w9KjMS0+eDF9/DS1bQna2jDA7dSpYMrPTTz/JNMUnn0hStizpGteypfxbCnMiiSbnoKTJWTnP9u3w0UcwYwasXy+JLSNDRtW33ipTH3/ucLXX4cNyovXSpZCW9mc71W7d5ESSzp0hKsq9e2hyDkqanJWzbd8uo9D582HVKvnYyZMyV9u8ufT4aNBARtbenqs+c0amKtaulWOgFi+GX3+VZvclSsgLR6tW0lI1tymZgQMH8vHHHxMWFkatWrWYOnUq5cqVy/e+mpyDkiZn5T+ysuS0jy++gJ9/ltF1s2awd6+MrDt2hPBwKde7+mpZZKxSRT5WunT+ydsYGQlnZEhVxa+/SplbaqpMsfz8M5zLka1by9FQf/+7LGi2aJF/WeDSpUu5+eabCQ0N5cknnwRg1KhR+f67NTkHJU3Oyn9lZ8PGjTKKXb1aEuvBg9Lf+MgRScy1a0vSDQ39s1nQzp3yZ7VqklA3b4bff5cFu8hIefzKlTIq37ABypWTKZUyZeCaa6BpU/mzaNHCxz5//nzmzp3LzJkz832sJufg8d1339GzZ082btxYAggB1gCdTR69NTQ5K7+SlQW7dsHu3XDokGx6+eMPOH4cjh2DffvgxAmoWFFqrg8flmmJKlXgyislGV95pYy8q1WTz3nanXfeSefOnenevXuun09KSiIpSc74PH78OJs2bfJ8EMqRhgwZQmJi4ligBJBujBmR12M1OSvlIlc64CUmJpKSksK8efMc3Ttc2SM7O5tixYr9AJwA/maMybMrTb4tQ5VS4tNPP73s56dNm8aCBQtYvny5JmaVq99//x3gCqTRfnEgz9M6dfu2Uh6wZMkSRo0aRXJyMiULU/CsgkLv3r0BhgIzgcuuGOvIWSkP6Nu3LydPnqRNmzaAHN81ceJEm6NSTjJ9+nRCQ0MxxrxrWVYIsMqyrJuNMStye7zOOSullG9pVzqllPJXmpyVUsqBNDkrpZQDaXJWSikHcrdaQ4s5lVLKC3TkrJRSDqTJWSmlHEiTs1JKOZAmZ6WUciBNzkop5UCanJVSyoH+H5cQSGTK7ERqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWcAAAD3CAYAAADBqZV6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXiU1d3G8e+ByCKLrUAIJApFaIRASGkQLC2KslURpHRBoIpRUVHbqpT6qlGwIKhYFkElLLIIRWlBUJBVqWCtMRRS2TQVohACJOxEWZKc948TKUICCZnMM8v9ua65Mpk888wPGO6cOc9ZjLUWEREJLJW8LkBERM6lcBYRCUAKZxGRAKRwFhEJQApnEZEApHAWEQlACmeRUjDGDDPGvO51HRI+FM7iKWNMpjGmcymOW2OMudsfNYkEAoWzhD1jTITXNYicTeEsAcEYM9AYs84YM8YYc9AYs8MY8/Oin40EfgZMNMYcM8ZMLHr8amPMSmPMAWPMZ8aYX59xvjrGmLeNMUeMMZ8YY0YYY9ad8XNrjHnAGJMBZBQ9Nt4Ys7PoOeuNMT/z61+CyBkUzhJI2gGfAXWB54FpxhhjrX0CWAs8aK2taa190BhTA1gJzAUigduAl40xcUXnmgTkAVHAHUW3s91a9Jotir7/BEgALi8673xjTDXf/zFFLkzhLIHkS2vtFGttATATaADUL+HYHkCmtfY1a22+tfbfwN+BXxpjKgN9gKettV9ba7cUne9so6y1B6y13wBYa1+31u4vOt+LQFUg1sd/RpFSUV+bBJI9396x1n5tjAGoWcKxjYB2xphDZzwWAcwG6hXd33nGz868X+xjxphHgbuBhoAFauNa8SJ+p3CWYHH28ok7gX9Ya7ucfWBRyzkfiAE+L3r4ivOds6h/+U/AjcBma22hMeYgYHxQu0iZqVtDgsVeoMkZ378D/NAY81tjzCVFt7bGmOZF3SILgGHGmEuNMVcDt1/g/LVwgZ4DRBhjnsK1nEU8oXCWYDEe15980BgzwVp7FOgK9AV247pEnsP1EwM8CFxW9Phs4K/AifOcfznwLq6l/SVwnOK7QkT8wmixfQllxphIoANwP1ADeBhIs9YWelqYyAUonCUkGWNuB+4FqgC7gG7ACtynxauAvwEvWmuPeFakyHkonCUkGWNm4yauRAL7gMnAaGutLZoR2AOobK39u4dlipRI4SwiEoB0QVBCmjFmtjHmsjO+b2yMWe1lTSKlUd5xzmp2S0B79dVXGTt27IClS5eSlZVFs2bNePHFF0HvXfFOqcbOl7dbQ29wCXjr1q2jU6dO1K1blw0bNhAVFeV1SRLeShXO6taQkDZ79mySkpKYNWsWAwcO5KabbiI9Pd3rskQuSC1nCWm33norKSkpREZGApCamsq9997Lhg0bPK5Mwpi6NUSKc/LkSapUqeJ1GRK+1K0h4WvEiBEcOHCg2J9VqVKF9957j3feecfPVYmUnlalk5DUqlUrbrnlFqpVq0abNm2oV68ex48fJyMjg40bN9K5c2cef/xxr8sUKZFazhKS/va3v/Hhhx/SrVs34uLiKCgooHbt2gwYMIDU1FTGjh1LvXr1LurcSUlJREZG0rJly9OPDRs2jOjoaBISEkhISGDp0qW++qNImFKfs4SkFi1a8O6779KzZ0/ef//9c35++eWXX/S5P/jgA2rWrMntt9/Opk2bABfONWvWZMiQIRd9XgkbpepzVreGhKT77ruP7t27s337dhITE08/bq3FGMP27dsv+twdO3YkMzPTB1VK2LEWTOn2b1C3hoSk3/3ud2zdupWkpCS2b99++rZjx45yBfP5TJw4kfj4eJKSkjh48GCJx6WkpJCYmEhiYiJxcXElHichaOjQUh+qbg2Ri5CZmUmPHj1Od2vs3buXunXrYowhOTmZ7Oxspk+ffsHzJCYmkpaWVtHlSiD4+muIjoaDBzWUTsRf6tevT+XKlalUqRL33HMPqampXpckgeavf4VDhy58XBGFs4gPZGdnn76/cOHC74zkEMFamDQJyvC+0AVBkTK67bbbWLNmDbm5ucTExDB8+HDWrFnDxo0bMcbQuHFjJk+e7HWZEkj+9S/YsAFeeaXUT1Gfs4iH1OccJgYMgLffhqwsqFlTfc4iIp7buxfefBPuuANq1iz10xTOIiIVaepUOHUKBg8u09MUziIiFSU/H159FTp3hquvLtNTFc4iIhXl7bdh1y544IEyP1XhLCJSUSZNgiuugB49yvxUhbOISEXYuhVWr4b77oOIso9aVjiLiFSESZOgShW4++6LerrCWUTE1w4dgqVL4eGHoWj/yrJSOIuI+NqUKbBjB/zmNxd9CoWziIgvnToFEybADTfAj3500afR2hoiIr705ptu+Fw511dRy1lExFeshRdfhObNoXv3cp1KLWcREV9Zs8atPjdlClQqX9tXLWcREV8ZM8aNzhgwoNynUjiLiPjC1q1u+NwDD0C1auU+ncJZRMQX/vIXF8r33++T0ymcRUTKa+9emD3brdlcr55PTqlwFhEpr0mToH17NyPQRxTOIiLlcfAgjBsH9etDbKzPTqtwFhEpjwkT4OhReOIJn55W4SwicrGOHIHx46FXL4iP9+mpFc4iIhfr5Zddt8aTT/r81ApnEZGLkZfnpmp37w6JiT4/vcJZRORipKRAbi4kJ1fI6RXOIiJldfw4vPACdOoEP/lJhbyEFj4SESmr6dMhOxvmzKmwl1DLWUSkLE6ehNGjoUMHuP76CnsZtZxFRMpi1izYudP1ORtTYS+jlrNIGSUlJREZGUnLli1PP3bgwAG6dOlCs2bN6NKlCwcPHvSwQqkw+fkwapQbndGtW4W+lMJZpIwGDhzIsmXLvvPY6NGjufHGG8nIyODGG29k9OjRHlUnFWrePNi+3Y3QqMBWMyicRcqsY8eOXH755d95bNGiRdxxxx0A3HHHHbz11ltelCYV6eRJGD4c+vaFHj0q/OUUziI+sHfvXho0aABAgwYN2LdvX4nHpqSkkJiYSGJiIjk5Of4qUcpr6lT473/ht78t9xZUpaFwFvGzQYMGkZaWRlpaGvV8tPavVLC8PHjmGejYEX7+c7+8pEZriPhA/fr1yc7OpkGDBmRnZxMZGel1SeJL48e7BfUXLKjwvuZvqeUs4gM9e/Zk5syZAMycOZNevXp5XJH4zP798Nxz0LNnhc0GLI7CWaSMbrvtNq699lo+++wzYmJimDZtGo899hgrV66kWbNmrFy5kscee8zrMsVXRo926zWPHOnXlzXW2vI8v1xPFgl3iYmJpKWleV2GlGTXLmja1I3QmDHDV2ctVb+IWs4iIiUZPhyshWHD/P7SCmcRkeJs2+YWOLr/fmjc2O8vr3AWESnOpElQsyY8/rgnL6+hdCIiZ1uzBiZOhDFjwKNhkWo5i4icqaAAfv97aNQIBg/2rAy1nEVEzjR1KvznPzB/PlSv7lkZajmLiHzr4EF44gm47jro08fTUhTOIiLfGj7cBfT48X6bpl0ShbOICMCWLe4i4KBB0Lq119UonEVEsBYefhhq1XKrzwUAXRAUEVmyBFasgHHjIECWcVXLWUTC2zffwNixbp1mD4fOnU0tZxEJb3/+M7z3nrtdconX1ZymlrOIhK9PP4UXXoCBA6FTJ6+r+Q6Fs4iEp8JCuPde+N733DTtAKNuDREJT5Mnw0cfwaxZUKeO19WcQy1nEQk/u3fDY4/BjTfCgAFeV1MshbOIhJ/f/x5OnIBXXvF8JmBJFM4iEl6WLoXFiyE5GZo187qaEimcRSR85OZCUpLrzvjjH72u5rwUziISHqyF++6DAwfcjtpVqnhd0XlptIaIhIc5c+Dvf3fBHB/vdTUXpJaziIS+nTvhwQehQwcYMsTrakpF4Swioa2wEO68E/LzYeZMqFzZ64pKRd0aIhLaJk6E1ashJQWuusrrakpNLWcRCV3btsGf/gQ33wx33+11NWWicBaR0HTiBDz9NNSo4TZtDdDJJiVRt4aIhKZHHoE333STTqKivK6mzBTOIj7UuHFjatWqReXKlYmIiCAtLc3rksLT3Lnw8svw6KNuEf0gpHAW8bH333+funXrel1G+Nq8Ge65B376Uxg1yutqLpr6nEUkdBw9Cn36uI1a33gjoHY2KSuFs4gPGWPo2rUrP/7xj0lJSSn2mJSUFBITE0lMTCQnJ8fPFYYwa92IjIwMmDcPGjb0uqJyMdba8jy/XE8WCTW7d++mYcOG7Nu3jy5duvDSSy/RsWPHEo9PTExUv7SvTJjglgIdNcqt1Ry4SjVsRC1nER9qWNRai4yMpHfv3qSmpnpcUZj46CN38e+WW2DoUK+r8QmFs4iP5OXlcfTo0dP3V6xYQcuWLT2uKgxkZcGwYRAb66ZnVwqNWNNoDREf2bt3L7179wYgPz+ffv360b17d4+rCnHHjrnWckaGaz1///teV+QzCmcRH2nSpAnp6elelxE+CgqgXz9IT4d33oEQ+5SicBaR4DRkCLz9NkyaFLQTTc4nNDpnRCS8vPwyjBsHf/gDDB7sdTUVQuEsIsFl4UKYPt31NY8Z43U1FUbhLCLBY9Uq6NsXqlVz62cEycL5F0N9zsGuoACys2H3bti7F/bsgZwc+PprtwMEuCva33zjvr/8cjel1Vo35KhBA3eFu1YtiIx03zdsGNTTXiVEffQR9OrlhswtXgw1a3pdUYVSOAcLa+G//4UtW2DrVvjnP+HLL13YWutCOjcXPv/chXCTJhAdDZdd5tYb+OYbd0xUFJw86QJ961Zo3Bjq1IEjR9zC5Nde6447dcoFdtOm0KYNxMRAXBz84Ach3VqRAJWeDjfd5BoOK1a4932IUzgHqsJC2LjRtRaWLXNfIyNda2HXLtfabdcOmjVzb9ioKPfz+vXdG7datQu/hrWQlwf797vt4nfvdvdzc+Grr2D9ercW7urV7hdBhw6wYYMbstS+PbRqBW3bQosWamlLxdm2Dbp2de/9VauCcm3mi6G1NQLJ0aOwfLn7yLZ/vwvGKlVcALdrB9df74KxRQuoWtV/dR0+7P6DfP65C+zPP3ct98OH3c/j4lx3yM9+5m7t2sGll/qvviCmtTUuYNMmt8VUs2ZuyFxsrNcV+UKp1tZQOHvtxAkXyLNnw8cf/69VPGAAXHMN3HCDC75AU1gIX3wBqanwr3/B2rXwn/+41vgll0Biorua3ratW1e3NC35MKRwPo+0NOjWzb13Vq50jZLQoHAOaBkZ8MorMGuWaw1v3Qq/+pW7degAEUHY43ToEHz4oQvqtWtdq3/NGqhe3bX6u3d3/9l++MOg28+toiicS7BunetjrlPHdas1aeJ1Rb6kcA441ro33QsvuNZmXp57A959N9x4Y3AG8vnk5blwXr7c9ZtnZLjH27Z1t169XGhXqeJllZ5SOBdj5Ur33rjyStfHHBPjdUW+pnAOGAUFsGjR/0K5Th148EG4776wubgBwPbtLqg//RRmzHAjSGrXdlNve/VyX7/3Pa+r9CuF81kWLIDbboPmzd2ojMhIryuqCApnz1nr5v7PmeN2AW7SxK05O3CgLph9841rFS1a5P6O9u1z3Tm1asEvf+nCOgz24VM4F7EWxo6Fv/wFEhLcNZgQWmHuLFps31MbN0Lnzi5kPv0U/vpXN8ph8GAFM7h+6FtugalT3RC+f/7TtZw//9x180RFQZcu7j9pbq7X1UpFOnnSbcj66KNunP0bb4RyMJeaWs6+lp0NTz4Jr73mxhsPHw6DBmkccGlZ636x/e1vMH++a2FnZ7s++d/8Bnr3Dqn/uGHfcs7NdRuyfvABJCe7RfNDZLH881C3hl+dOOHGYT71lGsJPPSQC+kQChK/Kyx0nzreeMPdtm93FxKjoly/ZM+eUKOG11WWS1iH86ZN7pNlVpZbyKhfP68r8heFs9+kp8Ptt7uP6g0bwvPPu2nP4jvWugkwy5bBq6+6/9CXXgq33uqCumvXoBz1EbbhPHu224Q1Ksot/9mundcV+ZPCucLl58Nzz7muizp1ICXF9aNKxSosdEMS5851XR8HDrhJLz/+MfTv7y4sBslH47AL57w8N1Jpxgzo2NH9G0ZHe12VvymcK9TWrXDHHfDJJ24Jw4kTXUCLf5086UZ9zJkDb73lVuO78krXmu7XD+Ljva7wvMIqnNPT3f+Vzz5z/cvJyaE3tr90NFqjQljr+sc6dHB9oG++6UZiKJi9UaWKm8gzZ45bMvX1192MyzFj4K673OJMo0ZBZqbXlYavkyfdp8v27d3CXKtWue/DM5hLTS3nsjh+HO6/330ku+suGDnSvdkk8OTkuE0/p01zU8rB/ULt18+Now6QyQ0h33Jevx7uvNNd2O3fH8aPV0NGLWcf++ort+LajBluREZKioI5kNWr50Jh3TrYscO1ng8fdrM0GzZ063zMmOHWAxHfO3oUnn3WXejbv9+ttPj66wrmMlDLuTTefx9+/Ws3XG72bDf8R4LTpk3uItS8eS60q1Rxk146dHAXc2vV8ms5IddyLix0i3n93/+5gE5KgmeeCbtp+ReglnO5WQsvveRmqtWt6y7+KZiDW8uWrkX3xRdunZPBg92sxP79XWu7d293DeHoUa8rDT7r1rmW8p13QqNGrm95wgQF80VSy7kkhYVuOulHH7mPwTNn+r1VJX5SWOimj7/5ppuZmJ3tJg9df72b6NKjR4Wt8xESLee1a+HPf3YjZTIzYfRo17cfJMMZPaCW80XLz3cX/MaNc3P9589XMIeySpXchgATJrjNDtauhd/9zi32fued7tpCp07uYpZGfTjWupbxdde58crp6e5Tx2efuY0iFMzlppbz2Y4fd7/1Fy50fWVPPqmF4cOVtfDvf7vx0wsXwubN7vEbbnC71Nx0k/vlXY4hYUHXcj5xwi3r+fbbrvunYUMYOtQtXKQFvUpLk1DK7OhR99t/9WrXinroIa8rkkCSkQHvvuuCet069wnrssvc1PFevVzru1GjMp0yaMJ52zb3S2rsWLe86/XXu4Wo7rzTv/tZhgaFc5ns3+9aQuvXuxXlfvtbryuSQHb4sPtY/+67biPeqCjYsIGvo6JY8vXXfFStGk3uuosHn332vKcJ6HDeudP92aZNcxfDGzWC1q3dRdQuXdR1cfEUzqW2e7dr/fz3v+6iUM+eXlckwcRa2LSJwtWrWfXEE9wYEUHlI0cAyGvblhotWsBPfuJGMsTFfacbJKDCOT/fdeMsW+aGj65ZAzVruo1Vf/MbNyU+EDcbDj6lCmfNn9y92138+eor1wrq1MnriiTYGAOtWvHxsWO8+NOf0nXpUkhPZ2VyMlfu2EHs0qWwZIlbu7h6dTeZqVkz1wrNy4MjR9x2Xf529Chs2OBm7737rgvmffvc6JX+/d2oi9693Ya84nflajnHxcXZ6tWr+7Ac38vJyaFevXrF/zA/311dPnnSvQE9XBv4vHUGENVZsoMHD3LkyBEaFfU779+/n7y8PK688kp3IS0vD/LyOHHsGJd8/TWVgC1AC3CbMVx66f++Vq7sJshccom7XWwXQmEhnDrl3uPffj150m1icOzY/97z+fluRFLt2u7rWRc59e/uO+vXr99srW15oePK1XKuXr164HwkK0GJHxuPHXPbSFnrdvu94Qb/F3eGgPp4ex6qs2Tz589n+fLlTJ06FYDZs2eTmprKSy+9dO7BhYWQmUn9Fi1IGzbMTYTZts11rTVq5FY9zM2F2Fg35nrPHheoP/rR/yZ1nDjhgvzSS1245ue78x486Na7rlrVhTy493utWm6z4apV3bjtFi3caJPWrS+41oj+3X3HGHO8NMeFZ7fGiRNukfa0NPj73z0PZgkNMTEx7Ny58/T3u3btomHDhsUfXKkSNGnCscqV3aLzZzpyxF2M27nTrbR37Jjrbti507VojXFrghw65AL7xAnXyIiIgGrV3PC26GgXwPXqua8NGkBMjAt+rQYXFMLvXyk/341jXr3azfrTdGzxkbZt25KRkcGOHTuIjo5m3rx5zJ07t+wnql3bXTiMi/N9kRI0yhXOgwYN8lUdFeY7NVrrNltdsMDN/rv9du8KO0sw/F2C6jyfiIgIJk6cSLdu3SgoKCApKYm4CwRs3QqaFu5r+nf3qZTSHBReQ+mGD4f33nMD6IcP97oakaDoIxWf09oa3/Hmm27b9auugqef9roaEZHzCo9wTktz+/116ACvvKKZTSIS8MqVUsOGDSM6OpqEhAQSEhJYunSpr+rynawsd9Gvfn1e7twZU60aubm5XldVrOTkZOLj40lISKBr167s3r3b65KK9cc//pGrr76a+Ph4evfuzaEA3U1k/vz5xMXFUalSpYDrOli2bBmxsbFs2rSJ0aNHe11OiZKSkoiMjKRlywsOy/XMzp076dSpE82bNycuLo7x48d7XVKxjh8/zjXXXIMxJt0Ys9kYc/6+VWvtRd+efvpp+8ILL9iAlZdnbWKitTVq2Ozly23Xrl3tlVdeaXNycryurFiHDx8+fX/8+PH23nvv9bCaki1fvtyeOnXKWmvt0KFD7dChQz2uqHhbtmyx27Zts9ddd5395JNPvC7ntPz8fNukSRP7xRdf2DZt2tj4+Hi7efNmr8sq1j/+8Q+7fv16GxcX53UpJdq9e7ddv369tdbaI0eO2GbNmgXk32dhYaE9evSote463yXAx0B7W0K+hu7ne2vdFjnr18PcuTyYksLzzz+PCeDlP2ufMYU3Ly8vYGvt2rUrEUVjZdu3b8+uXbs8rqh4zZs3JzY21usyzpGamkrTpk1p0qQJxhj69u3LokWLvC6rWB07duTyyy/3uozzatCgAW3atAGgVq1aNG/enKysLI+rOpcxhpo1a3777SVFtxIHVZQ7nCdOnEh8fDxJSUkcPHiwvKfznZdecmsFjBrFYiA6OprWrVt7XdUFPfHEE1xxxRXMmTOHZ555xutyLmj69On8/Oc/97qMoJKVlcUVV1xx+vuYmJiADJNglJmZyYYNG2jXrp3XpRSroKAAY8xGYB+w0lr7cUnHXnCcszFmFRBVzI+e2LNnD8nJyRhjSE5O5tFHH2X69OkXX3k5dO7cmT179gDwo7w8XsvMZHnt2uTHxvLsyJGsWLHCk7rOdmadZxo5ciS9evVi5MiRjBw5klGjRjFx4kSGezTk70J1fns/IiKC/v37+7u800pTZ6CxxQxfDdRPScHk2LFj9OnTh3Hjxn3nU2ggqVy5MtbaBGPM94CFxpiW1tpNxR17wXC21nYuzYvec8899OjRo4yl+s6qVavcndxcSEiApk25ef16Pv3yS3bs2HG61bxr1y7atGlDamoqUVHF/c7xU50X0K9fP26++WbPwvlCdc6cOZN33nmH1atXexospf37DCRlmuYtpXLq1Cn69OlD//79+cUvfuF1ORdkrT1kjFkDdAeKDedydWtkZ2efvr9w4ULvr+gWFrpZfzk5blxz7dq0atWKffv2kZmZSWZmJjExMfz73//2JJgvJCMj4/T9xYsXc/XVV3tYTcmWLVvGc889x+LFi7lUWxOV2ZnTvK21zJs3j55aQ/yiWWu56667aN68OY888ojX5ZQoJyfn9MgmY0x1oDOwrcQnlHSlsDS3AQMG2JYtW9pWrVrZW265xe7evdsP1zzPY/Roa8Hal18u8ZBGjRoF7GiNX/ziFzYuLs62atXK9ujRw+7atcvrkop11VVX2ZiYGNu6dWvbunXrgB1VsmDBAhsdHW2rVKliIyMjbdeuXb0u6bQlS5bYZs2a2SpVqtgRI0Z4XU6J+vbta6OiomxERISNjo62U6dO9bqkc6xdu9YCtlWrVqffk0uWLPG6rHOkp6fbhIQEC/wH11p+yp4nX0Nn+va6dW5adp8+MG+eNmWVoKDp22EpjKZv5+ZC377QuDFMmaJgFpGgF/zhbK3bJbtpU5g/35vtfkREfCz4w/n11103Ro8ebpcIEZEQENzhvGuXazV36AAPP+x1NSIiPhO84Wwt3HWX27Ryxoz/7ZUmIhICgnebqsmTYcUKmDTJ9TeLiISQ4Gw5f/EFDBkCXbrA/fd7XY2IiM8FXzgXFMCdd7pujGnTNGxORILCJ598Qnx8PMaYasaYGkVrOpc4rTr4wnnGDPjwQ5gwAc5Y2UtEJJC1bdv222n6I4DngddtCYseQbBt8PrVV9C8OfzqV/Daa2o1S9DTDMHwcvLkSapWrfof4DjwE2ttQUnHBlfL+aGH3NfhwxXMIhJ0Dhw4AFATqAVUO9+xwRPOb70Fixe7HbQbNfK6GhGRMhs0aBBAMjAHeO58xwbHULqjR12ruVUr+MMfvK5GRKTMZs2aRUREBNbaucaYysA/jTE3WGvfK+744OhzfuQRGDfOXQi89lq/vKSIP6jPOSyFyKp0GzbA+PEwaJCCWUTCRmCHc0EB3Hsv1K0Lo0Z5XY2IiN8Edp/zlCnwySdu5bnvf9/rakRE/CZwW87798Ozz0JSEvTr53U1Iuc1bNgwoqOjSUhIICEhgaVLl3pdkgS5wG05Dx8OWVluKVCNaZYg8PDDDzNkyBCvy5AQEZgt561b4eWXXX+z1zt6i4h4IDDDecgQqFnTtZ5FgsTEiROJj48nKSmJgwcPlnhcSkoKiYmJJCYmkpOT48cKJZgE3jjn5cuhe3cYMwYefdTnpxe5WJ07d2bPnj3nPD5y5Ejat29P3bp1McaQnJxMdnY206dPv+A5Nc45LJWqnzawwjk/H1q3hhMnYPNmqFrVp6cX8YfMzEx69OjBpk0lLjh2msI5LAXhJJQpU2DLFnjhBQWzBJXs7OzT9xcuXEhLXSuRcgqc0RqHDkFyMlx/Pdx6q9fViJTJ0KFD2bhxI8YYGjduzOTJk70uSYJc4ITziBFw4ACMHauhcxJ0Zs+e7XUJEmICo1sjMxNSU+HBByEhwetqREQ8Fxgt5+HDXTjPnet1JSIiAcH7lvO2bTBrFjzwAMTEeF2NiEhA8D6cn3oKLr0UHnvM60pERAKGt+G8YQPMn+/Wz6hXz9NSREQCibfh/OSTbilQzQQUEfkO78L5ww9h6VL405/gsss8K0NEJBB5E87WwhNPQFSUGz4nIp4qGBYAAAQKSURBVCLf4c1QulWr4B//gJdegho1PClBRCSQ+b/lbC08/jg0agT33OP3lxcRCQb+D+elS92iRk89pcWNRERK4N9uDWvdbMDcXLj9dr++tIhIMPFvOK9c6XbTTkmBiMCYOS4iEoj8260xYoSboq1Ws4jIefmv+frBB7B2LUyYoL5mEZEL8F/LecQIiIyEu+/220uKiAQr/4Tzxx+7/uZHH4Xq1f3ykiIiwcw/4TxypFtD4/77/fJyIiLBruLDOT0d3n4b/vAHqFWrwl9ORCQUVHw4z54NsbHw0EMV/lIiIqGiYsM5M9Nt2HrLLa5bQ0RESqViw3n8eKhUCX7/+wp9GRGRUFNx4XzoEEydCn37am9AEZEyqrhwTkmBY8e0y4mIyEWomHA+edLNBLzhBkhIqJCXEBEJZRUzffuNNyArC6ZMqZDTi4iEOt+3nK2FF1+EFi2ge3efn15EJBz4vuX83ntu4sm0aWCMz08vIhIOfN9yHjMG6teH/v19fmoRr82fP5+4uDgqVapEWlrad342atQomjZtSmxsLMuXL/eoQgkVvg3nTZtg2TK3o7aWBZUQ1LJlSxYsWEDHjh2/8/iWLVuYN28emzdvZtmyZQwePJiCggKPqpRQ4Ntw/stf3KpzWuBIQlTz5s2JjY095/FFixbRt29fqlatyg9+8AOaNm1KamqqBxVKqPBdOOfmwpo1MGQI1Knjs9OKBIOsrCyuuOKK09/HxMSQlZXlYUUS7Hx3QXDaNNixw80IFAlinTt3Zs+ePec8PnLkSHr16lXsc6y15zxmSrggnpKSQkpKCgA5OTnlqFRCmW/CuaAAXnkFOnVyQ+hEgtiqVavK/JyYmBh27tx5+vtdu3bRsGHDYo8dNGgQgwYNAiAxMfHiipSQ55tujSVL4Msv4YEHfHI6kWDTs2dP5s2bx4kTJ9ixYwcZGRlcc801XpclQcw34TxpEkRHQwkf+URCxcKFC4mJieGjjz7i5ptvplu3bgDExcXx61//mhYtWtC9e3cmTZpE5cqVPa5Wgpkprq+sDCyff+4W03/mGUhO9lVdImEhMTHxnPHSEvJKNTuv/C3nV16BSy6Be+4p96lERMQpXzjn5cFrr0GfPhAV5aOSRESkfOE8dy4cPqwLgSIiPla+cJ40CeLjoUMHH5UjIiJQ3nHO6ekwebJWnxMR8bHytZwvu0yrz4mIVIDyhfPAgVCjhm8qERGR08oXzn/+s4/KEAlPdevW9boECVDln4QiIiJl4adJKCIi4nMKZxGRAKRwFhEJQApnEZEApHAWEQlACmcRkQCkcBYRCUAKZxGRAFTeDV614pGISAVQy1lEJAApnEVEApDCWUQkACmcRUQCkMJZRCQAKZxFRALQ/wOeZzK3eFMz5wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#\n",
    "sp.plot(quad(x).subs(a, 1).subs(b, 2).subs(c, -1),(x, -5, 3), line_color = 'blue', title = \"Function\")\n",
    "sp.plot(intQuad.subs(a, 1).subs(b, 2).subs(c, -1).subs(x, x),(x, -5, 3), line_color = 'red', title = \"Integral\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here they are plotted together. (Getting two different line colours is a little convoluted!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWcAAAD3CAYAAADBqZV6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVhVVfcH8O8GnBAVB1ABh8wZUlIcGrQcUnLMLFPTMkwzK8vyzSanEjX7WZY4oWnOmFP6mqk5R5moSZqW+SooAqIgDogy3fX7Y4GiMly4595z7mV9nuc+Mpy7z+IK6+6zz9p7KyKCEEIIY3HSOwAhhBD3k+QshBAGJMlZCCEMSJKzEEIYkCRnIYQwIEnOQghhQJKchVUppeYppcbpHUdBlFJ1lVKklHLJ5/vHlVJP2jisYlFKfaeUmqx3HMJykpwdmFIqWil1UymVkuvhZcXzDVFKhef+GhGNIKLPrHVOWyAiXyLaY86x2a95ZyuHJEoASc6OrycRueV6xOkdkCie/Hr2wjFJci6BlFJPKqXO3/O12z0+pdREpdT3SqmlSqnr2Zf1AbmOraWUWq+UuqSUSlJKhSilmgCYB+CR7B76lexj77rMVkoNU0r9Tyl1WSm1KXdPPntoYYRS6pRSKlkpNVsppfL5GVorpfYrpa4opeKzYyhtTltKKWel1P8ppRKVUmcAdC/k9TLrtVFKLQNQG8B/s1+D97O/3lYp9Vt2rH/mHiJRSj2glNqX3daO7DiXZ38vZ7hlqFLqHIBd2V9fo5S6oJS6mv1c34LiF/ZJkrPITy8AYQDcAWwCEAJwYgOwGcBZAHUBeAMII6K/AYwAsD+7h+5+b4NKqY4ApgLoB6Bmdhth9xzWA0ArAM2zj+uaT3xZAEYDqAbgEQCdAIw0s61h2d97GEAAgOcKeiHykOdrQ0SDAZzDnauV6UopbwA/ApgMoAqAMQDWKaU8sttaCSACQFUAEwEMzuN8TwBokiv+nwA0AOAJ4A8AK4oYv7ADkpwd3w/ZPbYrSqkfivC8cCLaQkRZAJaBExwAtAbgBeA/RHSDiG4RUXi+rdztRQCLiOgPIkoD8CG4p1031zHTiOgKEZ0DsBuAf14NEdFhIvqdiDKJKBrAfHASyy2/tvoBmElEMUR0GfyGURT5vTZ5GQRgS/bxJiL6GcAhAN2UUrXBbx7jiSg9+3XclEcbE7Nf65vZP/siIrqe/RpOBNBcKVWpiD+DMDhJzo7vGSJyz348U4TnXcj1cSqAstljnrUAnCWizGLE4gXuLQMAiCgFQBK4953fed3yakgp1VAptTn78v4agCngXnRBP0NOW14AYnJ97yyKJr/XJi91ADyf6w3yCoDHwVcOXgAuE1FqruNj8mjj9teyh2SmKaVOZ//c0dnfuvdnF3ZOknPJdAOAa84n2UMVHvkffpcYALXzSUaFLXEYB05WOectD76cjzXz3LnNBfAPgAZEVBHARwDyHJ/OQzz4TSZH7WKcPz/3vgYxAJbleoN0J6LyRDQtO44qSinXXMfXwv1ytzkQQG8AnQFUAg8tAeb/7MJOSHIumf4F9/a6K6VKAfgEQBkznxsBTirTlFLllVJllVKPZX8vAYBP7htz91gJ4BWllL9Sqgy4t3sge1iiqCoAuAYgRSnVGMDrRXju9wBGKaV8lFKVAXxQjPPnJwFAvVyfLwfQUynVNbvXWzb7hqwPEZ0FD3FMVEqVVko9AqBnIe1XAJAGvuJwBb+GwgFJci6BiOgq+ObZQnCv9QaA8wU+6c5zs8AJpD745td5AC9kf3sXgOMALiilEvN47k4A4wCsAyf4BwH0L+aPMQbci7wOYAGA1UV47gIA2wD8Cb6htr6YMeRlKoBPsocwxhBRDLin+xGAS+Ce9H9w52/vRfANzSTwTcPV4OSbn6XgYZhYACcA/K5h7MJAlCy2LxyZUsoTwGPg8d2bAP4CcIiITLoGlg+l1GoA/xDRBL1jEfqS5CwcklKqA3i4ogqAIwAuAigLoCG4x74WwAwiuqZbkACUUq0AXAYQBaALgB8APEJER/SMS+hPZhwJR9UNwLDsMrq7ZN/M7AHgKfAQi55qgIdVqoKHiF6XxCwA6TkLIYQhyQ1B4dCUUstyT9DInhK9U8+YhDCHpcMa0u0WhjZv3jx89dVXg7Zs2YLY2Fg0aNAAM2bMAOR3V+jHrJp0S4c15BdcGF54eDg6dOiAatWq4ciRI6hRo4beIYmSzazkLMMawqEtW7YMQUFBWLp0KYYMGYJu3brhzz//1DssIQolPWfh0J555hmEhobC09MTABAREYHXXnsNR45IQYTQjQxrCJGX9PR0lC6d3wxzIaxOhjVEyTV58mRcvnw5z++VLl0au3btwubNm20clRDmk0kowiE99NBD6NmzJ8qWLYsWLVrAw8MDt27dwqlTpxAZGYnOnTvjo48+0jtMIfIlPWfhkNauXYtff/0VXbt2ha+vL7KyslCxYkUMGjQIERER+Oqrr+DhYe4qqXcLCgqCp6cn/Pz8bn9t4sSJ8Pb2hr+/P/z9/bFlyxatfhRRQsmYs3BITZs2xU8//YRevXph9+7d932/SpUqxW573759cHNzw0svvYS//voLACdnNzc3jBkzptjtihLDrDFnGdYQDmnEiBEIDAzEmTNnEBBwe29aEBGUUjhz5kyx227fvj2io6M1iFKUOERA3nsW30eGNYRDGjVqFP7++28EBQXhzJkztx9RUVEWJeaChISEoFmzZggKCkJycnK+x4WGhiIgIAABAQHw9ZWNs0uU9983+1AZ1hCiGKKjo9GjR4/bwxoJCQmoVq0alFIYN24c4uPjsWjRokLbCQgIwKFDh6wdrjCC1FTA2xtITpZSOiFspXr16nB2doaTkxOGDRuGiIgIvUMSRrNqFXDlitmHS3IWQgPx8fG3P96wYcNdlRxCgAgICQGK8HshNwSFKKIBAwZgz549SExMhI+PDyZNmoQ9e/YgMjISSinUrVsX8+fP1ztMYST79wORkcC8eWY/RcachdCRjDmXEC++CGzeDMTGAm5uMuYshBC6S0gA1qwBhgwB3NzMfppFyXnlSkueLYQQJcCCBUBGBjByJL791vynWZSchw0DTp60pAUhhHBgmZk8zvzUU4i82QgjR5r/VIuSc9mywMCBQHq6Ja0IIYSD2rgRiI1F2qtvYMAAoGpV859qUXL+9lvgjz+ATz6xpBUhhHBQs2cDtWtj9M4eOHkSWLbM/KdalJyfeQZ47TXgiy+AHTssaUkIIRzMiRPA7t043v51zA11xpgxQKdO5j/d4lK61FSgZUvg6lXg6FGgWjVLmhOiZJFSOgf2xhughQvRyPU8KtTzwP79QPYGPLYppXN15VmJSUnA0KE8EUYIIUq0q1dB/57C5mpDEJvugVWrbidms2lS5+zvD3z+ObBpU5EmwAghhGMKDYXa8TMmxL2Gb74BGjYsehOaTUIZNQro1w9YvJiHN4QQokRKT0f6F1/jsHtHNB7QAkFBxWtGs+Ts5ATMmgXExAAvvACkpGjVshBC2I8bi8JQ+lIsvin9H4SEmL22/n00nb7t6cmzBk+eBN58U8uWhRDC+MhESBz7f/gLfnj9h66wYDc07dfW6NABGD8eWLKEH0IIUVJsHrUdda4dQ8wLY9D2kWJ2mbNZZVW6rCygc2cgIgI4fBho3NiSUwjhuKSUznFERgJJLZ5C89InUOVKFJzK5lueod+qdM7OwIoVXGbXrx9w86Y1ziKEEMZw/TowodcRdKIdKDf27YISs9mstmSolxdPVTx2DBg92lpnEUIIfRHxEhb9YmYgs5wbyo8erkm7Vl3POTAQGDsWCA8HwsKseSYhhNDHggXAD9+cwwCnMLi8Phxwd9ekXasvtv/ppzyl+9VXeaq5EEI4isOHgbfeAr6s8zWXzL39tmZtWz05ly7N48/lywPPPSf1z0IIx5CcDDz/PFC/2hX0SQyF6t8fqF1bs/Ztsk2Vtzevv3HyJDB8uKy/IYSwb0S861RMDLClTyicbqQAY8Zoeg6b7SHYsSPw2WecpOfMsdVZhRBCe//3f7yW0Mxpt1Dn+E+cqf39NT2HTTd4/eADoHt3rt44cMCWZxZCCG2EhwMffcRDGiNLLQD27AFeflnz81hlEkpBLl/m9Z+zsoAjR4q2bYsQjkYmodiXuDigRQugTRtg2YJbqPjwg0D9+sDevUVpRr9JKAWpUoV3CU9I4HefrCxbRyCEEEWXnn6nqGHKFKDiusWcrcePt8r5XKzSaiECAnj/wcGDuec8ZYoeUQghhPlGjwb27we+/x7wrZ8GPD0VePRRvqFmBbokZwAYNAj45Rdg6lS+THjuOb0iEUKIgi1dyoUM773HY80IXcKlGgsXFn9N0ELYfMw5t7Q0XsXu6FHg998BPz9LWhPC/siYs/H98Qfw2GNA27bAzz8DLqZ03tqkRg3uShc9ORtzzDm3MmWAtWuBChV4J+/kZD2jEUKIuyUm8t6ojz0GrF4NuLiAFw06exaYMMFqvWZA5+QM8AJJ69YB584BL74oNwiF8QUFBcHT0xN+uS71Ll++jKeeegoNGjTAU089hWTpadi9jAwewvj7b2DaNN5MBBkZQHAw3zgLDLTq+XVPzgCPqc+aBfz0E78ZCWFkQ4YMwdatW+/62rRp09CpUyecOnUKnTp1wrRp03SKTmjlvfe4hHnBAs7FAHgtiqgoq/eaAQBEZMlDMyYT0auvEgFE69Zp2bIQ2ouKiiJfX9/bnzds2JDi4uKIiCguLo4aNmxoVjstW7a0SnzCMt9+y7lo9OhcX8zIIHrkEaIWLThhFZ9Z+VW3ao17KQWEhACpqcDIkUC9eprPhhTCahISElCzZk0AQM2aNXHx4sV8jw0NDUVoaCgA4NKlSzaJT5jv99+B11/n3ZymT8/1jcWLgYMHgc2brd9rhs7VGnmJjwdateLdvA8eBKpX1/oMQlguOjoaPXr0wF9//QUAcHd3x5UrV25/v3LlymaNO0u1hrGcP8+Ls505A/z6a64ZzKmpQIMGQN26PH/bsuRs/GqNvNSsyQuKJCYCffpwuZ0QRle9enXEx8cDAOLj4+Hp6alzRKKobtwAevfm3Lt27T1LS4SE8GzAadNs0msGDJicAZ6UsmQJlxDKEqPCHvTq1QtLsrebX7JkCXr37q1zRKIoTCbgpZd4k9awsHvmXCQn82y5bt2Adu1sF5S5g9P5PKxq4kQelJ8+3dpnEsJ8/fv3pxo1apCLiwt5e3vTwoULKTExkTp27Ej169enjh07UlJSklltyQ1BY/j4Y841X36Zxzc//JBIKaLISK1OZ1Z+NdyY812NE9C/Py+UtHEj0LOnNc8mhO3JmLP+Vqzg5SRefRUIDb1n1CIujlede/ZZYPlyrU5pn2POuSnFN0hbtOCF+v/8U++IhBCOZP9+ngH4xBPA7Nl5DCd/9hlPPPn0U5vHZujkDACurtxrTkoCevTgNzIhhLBUzgxsHx+epVy69D0H/PsvsHIl8NprXNtrY4ZPzgDvQbhuHXDlCido2SRWCGGJK1f4/l5O2XKem3785z+AhwcwbpzN4wPsJDkDPCFl9Woe2hgwQNbgEEIUT3o60LcvcOoUsH490LhxHgdt3841vcOG6TbZwm6SM8DvdCEh/E43erTe0Qgh7A0R59tdu3jDjw4d8jgoI4MTzIMPAu+8Y/MYcxhm+ra5Xn8dOH0amDGDX7u339Y7IiGEvfj0U144f9Ik3okpT/PmASdO8M2uMmVsGl9uhi6ly4/JxDunXLrEb2x9++oRhRCWk1I621myBBgyhDfKXrw4n4l+iYk8TbtVK2DbNmvNBrT/Urr8ODlxyWFWFq8B/csvekckhDCyn3/mpT8DA/OoZc5t/Hjg+nXgq69sNk07P3aZnAEusfvvf3kdkl69gOPH9Y5ICGFEhw/zHJJr14BVq/Iomctx9Cgwfz4vi+nra9MY82K3yRng8pdt24By5fgdMSZG74iEEEZy+jQXElSpAmzdCri753MgEY+RVq7MA9IGYNfJGQDq1OEdVK5dA55+WvYhFEKwhASga1ce/ty2jbfEy9eGDcDu3TwjsHJlm8VYELtPzgDQvDnwww9ct9i7N3Drlt4RCSH0dP060L07zyjevDmfWuYcN2/ynlQPPcR1dgbhEMkZ4HrFpUt52Ojtt4HMTL0jEkLoIS2NRygiI4Hvvwfati3kCVOmALVqATNnZm+vbQzGiUQDL7wAXL3KU+HT0oBFi7iyQwhRMmRmcgXXhg1c0dWjRyFPOHaMF9AfMADo2NEmMZrLoZIzwIvzX7jAC5pUqsRvhjpXxAghbMBk4r//deu4Em7AgEKekJXFwxju7sCXX9okxqJwuOQM8DolycmcmCtXBiZO1DsiIYQ1EfGw8eLF3DEza9b17NnAgQO8oHO1alaPsagcMjkrxdO7r1zhqhh3d12nyAshrOyzz7gz9vbbnJwLde4c8NFHXOJVaBdbHw6ZnAEea16wgEvsRo8GKlYEgoL0jkoIobWvv+aEPGQIj04UOoxJxIv0AMDcuYYd93To22UuLrxW9jPP8LtqWJjeEQkhtLRyJXe+nn2WO2NmFQCEhQFbtgDBwTxRwqAcOjkDvKjU8uU89jxoEG95LoSwf8uX89/0sGGcpM2qgktK4rGP1q2BN9+0eoyWcNhhjdzKl+dC9MBAHl4qVYonqwgh7NPq1by6XIcOfFVs9sqe48dztcDOnYCzs1VjtJTD95xzVKjA07xbtACefx748Ue9IxJCFMe6dVzL/PjjvFlJuXJmPnHtWl7Y6IsveDagwZWY5AzwTcFt2/j/5dln+WMhhP3YuBHo3x9o04avhsuXN/OJMTE8/tGyJfDGG1aNUSslKjkDXFb3889AkyZ8o3DnTr0jEkKYY+tWngXcogVfBVeoYOYTs7J4cDozkwenS5WyapxaKXHJGeDlA3fs4AQ9YQJ/LIQwro0bed32557jK96KFYvw5M8/B/bt40knDz5otRi1ViKTM8ATgrZu5TroHj24skYIYTxr13JSbtGCN3jOd03mvBw4wDcB+/cvYNNAYyqxyRkAPD15CVdfXx7i2LhR74iEELmFhXFebd0a2L69iIn5+nVg4EDAx8fQk03yU6KTM8C7qezcCTz8ML87r1mjd0RCCABYtoyrMh57rBhDGQDw1ltAdDSvnVGkrG4MJaLOuTA5Nwm7deN36bQ0vn8gRFHVrVsXFSpUgLOzM1xcXGRn7WJasgR45RWuY960qQhVGTnCwriR8eM5u9shSc7ZKlbkMehevYCXXuIEPXSo3lEJe7R7925UM+AqZ/Zi5kxg6lTuKH37bRHqmHNERwMjRgCPPMJLVNqpEj+skZubG9dO9uwJzJplyCVehXBYRJxLR48G2rXj5T+LnJhv3uQud/PmPJxhoJ1NikqS8z1cXXlrmwYNeH3YDz7gXxohzKGUQpcuXdCyZUuEhobmeUxoaCgCAgIQEBCAS5cu2ThCYzKZeKmLyZP5inX16iJMyc5BxBNN9u7lP94HHrBKrDZDRJY8HFZmJtFrrxEBREFBRBkZekck7EFsbCwRESUkJFCzZs1o7969BR7fsmVLW4RlaOnpRAMH8t/af/5DZDIVs6HPP+dGJk/WND4rMCu/Ss85H87OXH0zfjzvRfjcc3zFJERBvLy8AACenp7o06cPIiIidI7I2FJTgT59eOLetGnA9OnFrHjbsoUvc194gRfRdwCSnAugFO+kMmsW3zEODOQNZIXIy40bN3D9+vXbH2/fvh1+fn46R2Vcycn8N7VlC69HNHZsMRv6+29ebtLfn3tSdlbPnB/7HS23oTff5BmFL70EPPEEV3XUqKF3VMJoEhIS0KdPHwBAZmYmBg4ciMDAQJ2jMqaYGN4hytMTWLWKO7zFkpzM6/+WLcuzyFxdNY1TT4osu9tVom6Vbd/Oq9n5+nKBfMOGekck7F1AQECJq4WOjAS6dwdu3AA2bOBa5mLJzOTJCXv28FRf+6lnNqtrL8MaRdClC7BrF3/8yCPAL7/oG48Q9mbbNi6Tc3YGwsMtSMwA8P77PHts7lx7Ssxmk+RcRK1b880LDw+gc2feKkcIUbhFi7jH/OCDwP79gEXD8XPn8kpzo0Y57GwxSc7FkPPL9dhjvNDVxIlSCy1Efoh485GhQ4FOnTinentb0OCqVbxgfq1awIwZmsVpNJKci6lyZb4x+MorXNExeDBP+RZC3JGaytOwJ03imX+bNxdjAaPcfvyR78y3b1+EXV3tk+P+ZDZQujTP/W/QgEsrz57lGxw2XVYhKwu4cAE4fx64eBFISOBHairfMHF25n+J+Guurjz1qnRp4NYtDrZ6dV5ZpkYNfnh58feFsMD581xIceQI1y+/956FVW579/KEA3//Im4eaJ8kOVtIKeDDD4F69fgN/dVXeQqq5uWtRMDp08DRo8ChQ0BcHCfhvXu5y/744/zXcO4cJ+M6dTjJVqvGvQuleFv48uUBJycgI4Nvpjg5AY8+yteaVavy8V5ewP/+x+sTeHgAAQH8NV9fHtNxkgsuUbADB3iN9Bs3gP/+l8eaLXL4MC9688ADvEeVRd1v+yCldBo6eJBXtbt+HVi6lMvuis1k4pqjHTt4la3Vq4HGjYGTJ3kmTPfuvFKTjw+/M9Spw4m0Rg0eczGnV2EycW86KQlITATi44FLl/jzo0c56f/yC4/tnT7NX3d15a1jqlThzTIDAvidyIEvL63JEUvpli/nToq3Nyfmpk0tbPCff7jEw82NSzwsGrA2BLOuHyQ5aywujpPygQPAJ5/wWJvZHc2UFK41+uEHnvV04gT3gtu25aLqJ5/k3muTJlx0bytXrwL//gscO8aPq1eB9ev537Zt+WutWnHNafPmfKe0yAvwlkyOlJyzsviqceJE/lVdu5Yvxixy/Djw7rv8t7B7N1C/vgaR6k6Ss17S0oCRI7l0qEcP7klUqlTAwT/9xDc3kpK4kLpKFSAoiLdn6dABqFnTpvGbxWS6M8yybx/w22/cq963j3c3btMG6NiRH23bFmOJsZLBUZJzYiLvCHXyJA9n/N//abDJ9YED/IZfpgz/XTRurEmsBiDJWU9EwJw5wDvvcKf3+++503vbuXNcq3n4MI/9enjw9iu9e3PP0x6HCVJSgF9/5R7Orl38s5lM3MsfPJj/uLp35xfEQdY/sJQjJOeICL5Pd/Eib3CtSdnxzp38t1C9Ov991KunQaOGYd4vv7nL1+XzEIXYu5eoWzciV1eipUtMRLt2EfXpQ+TkxI++fYn++19eN9HRJCcTbdpE9M47RI8+yss5AkT16hG99RbRtm1Et27pHaWu7HnJUJOJaO5colKliOrWJTp0SKOG168nKl2ayM+PKC5Oo0YNxaz8KsnZBmLPm+iTJmtpBQYQAWSqWpXogw+Izp7VOzTbiooimjOHqHt3onLliB55hMjNjahfP6JVq4iuXNE7Qpuz1+R84wbRSy9xBnn6aaKkJI0a/u477rS0batho4YjydkQfv6ZKCCACKDTtdrTECyiAN9U+ucfvQPTWWoq0datRMOGEXl68q9iqVJEgYFEy5YRXbyod4Q2YY/J+fhxfj9VimjiRKKsLI0a/uor/j3o3Jno+nWNGjUkSc66OnCAqGNHfolr1yZavJgoM5O2bCGqWpU7jCtX6h2kQWRmEv3yC9F77xE98ABRrVpEzs78RxoaSnTpkt4RWo09JWeTiWjePL7oqVaNaMcOjRpOTyd64w2idu2Inn22JAx1SXLWRVTUnf2tPDyIZs6875ft3Lk7Q7AffUSUkqJPqIaUlUV05AjRhx8S1a/PL5KzM9Hw4UQLFxJdvqx3hJqyl+SclMS3SgCip57ScCg4IYGofXtueMwYfqN2fJKcbSori+ibb4jKlydq04Zo0iSia9fyPTw9neizz4hcXIgaNiQ6eNCGsdoLk+lOos55NytViqhXL77scIB3NXtIznv2EPn48Ev/xRcaDmMcOsRXSWXLEq1YoVGjdkGSs82cOHEneQQGEkVHm/3UnTv5F9/FhSg4uKR0HIrBZCKKiCB6910ib29+rR98kKh/f652SUvTO8JiMXJyzsgg+uQTHltu0EDDagwivq9QtiwP+R0+rGHDdkGSs9Wlp3NGLV2aqEoVoqVLi7V18OXLfIMFIHr8cR4ZEQXIyuIaxXHj+HUH+N8RI4j27dOwa2d9Rk3OJ08SDRrEL+0rr2h4f+7aNaKhQ/nq8sknS8yN33tIcraqEyduV2FQv35EFy5Y1JzJxJ2JihX5PuKCBRZsEV+SpKVxz3nAAC4mr1WLe2NjxxL9+afe0RXKaMk5M5Poyy+5U1utGtHatRo2Hh7ON3ydnIjGj3fM2n7zSHK2muXL+Zb1008TbdigadNnzvCklZwbL0UYIRHXrxOtWcMvoLMzv4h+fnxT1qCXI0ZKzqdO8ZUbQNSjB1FsrEYNp6VxXb+TE09ACg/XqGG7JclZc2lpRG++yS9b+/ZE8fFWOU1WFs/VcHPjx5w5dnWlbgwXLxLNnk302GN3qgEefZQoJIQrBAzCCMk5M5MrFsuVI6pUiWjJEg2v2v74g8jfn1//V18t8CZ5CSLJWVPnz/OMNoDrcW1wSRYVxaW+ANHLL5NMXCmuM2eIpkzhXnRO3XlgIM9G03lWot7J+c8/iVq35mGMgQP511wTSUlcu+zszG+OGzdq1LBDkOSsmd27eRZb+fJE339v01ObTDz+/MADfN9x/HiimzdtGoJjOXqUaPp0fkEBflGfeYZo3TpdZqXplZxTU3mkwcWFy/FXrNCot5yZyZd6VarwMMabbzryNOzikuRsMZOJaP58fvdv3JhvAuokPp57NgDPzdi+XbdQHIPJRPT777woU82anKzLluWZFitX2uzyW4/k/NNPXCyRU4mRmKhRw3v2EDVrxg136MBvhCIvkpwtkpVFNGoUUeXKvMKLQcbKfv6Za04BLvHV7DK0JMvM5BK8t97iRA1won71VR76sGLPz5bJ+cwZot69+cfr3ZsXSNTEvn1EXbtyYq5dm0s8pNSoIJKciy0tjTMfwD0rg92Nu3mTF5xp0YJvGAYHy1CHZrKyeCU6I6kAABTMSURBVJ2Pd9/lRJMzfbxjR54BqnH5jC2Sc2oq0YQJRGXK8Mjc559rMGfHZCLasuVOeYeHB99svXFDi5AdnSTnYrl27c5duM8/N3QP4NSpO+sd1KnDw+EGDtf+mEw8r/6jj4iaNLkzjvrQQzylPDzc4imd1kzOWVlc9dmrF/+ODBigwZXWrVtEq1ffqcCoVYvftCQpF4Uk5yJLSOCJJc7OvIqcndi5885QX7t2sk6H1Zw8yYtLPPkk30kD+HL++ee5Fq0YvWprJeddu4hatuQQW7bkCZUWiYzkYZ8qVbjerkULokWL7HbavM4kORdJXBz/kZUrxzPO7ExmJt+7bNeO/1f79tX1/qXjS07mHuR77xF5eVHOLi8p3t4UWrkyjfL0pFkfflhoM1on5+PHeS+DnE7t0qUWjMrFxBDNmsWJOKey5YUX+G60wYb67IwkZ7MlJPBla/nyfHPDjl29yuOLbm58Bf7yy4adHOc4TCai48cp68svaaurK2VWqnQ7Wac2bsz/CXPm8MpB99THa5WcT58mGjKEt4uqWJFo2jQeay6SzExeh3ziRE7COduKdevGQxdSEqcVSc5mSUzkMcRy5TS49jOOS5e4U1emDC/1+OabGk7HFXn67bffqEuXLryc28GD9HO3bnSodes7O73kVIEMHkw0ciTR/PnUsnFjfkctpuho3kzGxYX/rydPLsLeBFevcg3/1Kk8X7tSJZ4wohQP3QQH81CO0JpZ+dWi3bd9fX2pXLlyxX6+LVy6dAkeHh55fzMrC/j3X+DmTaB+faBiRdsGl0uBcVogPR2IjweSk3nD68qVgRo1gNKli9eeteLUmh5xJicn49q1a6hTpw4AICkpCTdu3EDt2rX5P+LGDSAlBWkpKSiVmgonACcANAUAV1fA2RkoUwYoVw4oVeruh5PTXedKSwMSE4GEBP68WjWgZk0+9C4mE5CRwY/09Lsf16/z9wGgQgU+d8WK/PE9u7/L/7t2Dh8+fJyI/Ao7zqWwAwpSrlw5w2/rnu/W89euAV268C/t5s1At262Dy6XfOPUyJkzwJQpwNKlnKgHDQLGjgUaNy5aO9aOUyt6xLlmzRps27YNCxcuBAAsW7YMERERmDVr1v0Hm0xAdDSqN22KQxMmAElJwK+/Av/7H+DmxgkyIgLw8gKqV+fPMzOReiUNsRec8XtiA7i5AQ/WSka9NtXg5p6dlW/eBC5c4E6Hhwc/z2Tid2aAk3zp0kDt2kCdOkDr1kDLlpzdCyD/79pRSt0y5ziLkrPdunED6N4dOHwYWLNG98RsC/XqAQsXAhMmADNmAKGhwJIlwBtvAP37A48+eufvVxSPj48PYmJibn9+/vx5eHl55X2wkxNQrx5SnJ2BDz+8+3vXrgHnzgHnzwPx8aArV3Dmr5uIOnwZ8X8lwMspAX4N09Go7i24JsYBqQQ4uQFly3Lv28ODk3q9epx0q1ThbnWtWpyU7+teCyMqeck5MxMYMQI4exZYsQJ45hm9I7KpWrWAmTOBjz8GZs8G5s8HQkKAgADgnXeA558v/pBHSdeqVSucOnUKUVFR8Pb2RlhYGFauXFn0hipWBPz8cKu+H1atAmZ+Bxw9CjRrBjw7Eej2JlC1qtbRC8Mxd3A6r8f8+fNtNYBebHfFaDLd2Xx1wQL9gsqDXq9lSgoXEjRqxC/L008Tffxx/hUe9vB/TqRfnD/++CM1aNCA6tWrR5MnTy70+Nq1a9/3tVOneK+AwEC6vST1woXFqL7QkPy/awfAcJJqjXtMn84/shn1pyVNVhbRjz9yfbRS/AgM5L0ESu6GFdaXU0qXlsYzPDt1otszxl95hWjHDpn16YCsX60BwKIn29TatXzN/sILwMqV9939FnecOwd8+y2PUcfFAY89BrRqBQweDDz8sIxNa4UIaNo0AJ06HUJkJN8PrF0bGDYMCAriYWPhkMz6CyoZyfn334EOHYAWLYCdO/nGiShUZiawbRsn6s2bubClaVPgtdeAwECgYUO9I7RP//wDhIUBy5cDp08HoGzZQ+jdG3jlFaBzZ76nJxyaWcnZou7jxIkT4e3tDX9/f/j7+2PLli2WNGcdUVFAr16AlxdmP/UUVLlySExM1DuqPI0bNw7NmjWDv78/unTpgri4OF3jcXHhopb167k6a948rpP+7LM4NGoElC17Co0bhyE8/Bose4+3jjVr1sDX1xdOTk66llcRAYcO8U3Ypk25MiY42IT4+N/g4nIeY8d+ibAwoGtX4yXmoKAgeHp6ws+v0LJc3cTExKBDhw5o0qQJfH198fXXX+sdUp5u3bqF1q1bQyn1p1LquFJqUoFPMHf8I6/HhAkT6IsvvrDVOE3RpaQQPfEEUdu2FLd7N3Xp0oVq165Nl8yeQmVbV3PNFPv666/ptdde0zGa/C1btpdmzMjMXscjiwDeTOStt3gVST1vXOV24sQJ+ueff+iJJ56ggzZeDerqVd6Z6ZNPeI2LnHHkDh2IZs7Molq12tLp06epRYsW1KxZMzp+/LhN4zPX3r176fDhw+Tr66t3KPmKi4ujw4cPExHRtWvXqEGDBoZ8PU0mE13n3XYAoBSAAwDaUj751XFL6YiA4cOBffuAbdvwVkgIpk+fjt69e+sdWb4q5pqheOPGDSiDDu4OGtQeAPDuu8DixVsxd2483N2HYsECYNYsHjUaOBBo1Ah44gkeTdKjtLZJkyY2O1dGBveOf/6ZH7//zsNCNWoAbdsCn34K9OzJJXD79x9AkyYVUa9ePSil0L9/f2zcuBFNmza1Wbzmat++PaKjo/UOo0A1a9ZEzZo1AQAVKlRAkyZNEBsba7jXUykFNze3nE9LZT/yvea0ODmHhIRg6dKlCAgIwIwZM1C5cmVLm9RGSAjf+AsOxqabN+Ht7Y3mzZvrHVWhPv74YyxduhSVKlXC7t279Q6nUOvXz8WoUS9g0CBgzhx+L9y6ledPjB3Lx5QvDzz+OA/7BwTwo1IlfeO21OXLnIB/+40fx47xzEuTid+MxozhCaiPPsqzonOLjY1FrVq1bn/u4+ODAwcO2PgncEzR0dE4cuQI2rRpo3coecrKyoKLi0skgPoAZhNRvv/xhSZnpdQOADXy+NbHFy5cwLhx46CUwrhx4/Dee+9h0aJFxY/cAp07d8aFCxcAAA+npmJRVBR+qVAB15s0wZTgYGzfvl2XuO6VO87cgoOD0bt3bwQHByM4OBhTp05FSEgIJk0qeFjKWgqLM+djFxcXvPjiiwB4SYiuXfkB8LoP+/YBe/YAR47wRDgirv5ITATat+clTR56iB/e3kWvBDEnTkuYTFy9cvw4Pw4d4p/l+nX++ZydgebNeZZlp078JlTITGhQHgP0Rr1KsicpKSno27cvZs6ceddVqJE4OzuDiPyVUu4ANiil/Ijor7yOLTQ5E1Fnc046bNgw9OjRo4ihamfHjh38wYUL3HWpXx8dDx7EsZgYREVF3e41nz9/Hi1atEBERARq1MjrPcdGcRZi4MCB6N69u27JubA4lyxZgs2bN2Pnzp35Jpbq1bl68fnn+fMrV3i5iGPHgPBwIDYWWLAA8PXlxFe5MvDgg5ywmzfnIYA6dbi8zMuLl4m491Tmvp75IQKuXuVYzp/nRJyYyAn41CleH+jECT72oYeAlBQuJ3zySf68VSu+MiiKIk3zFmbJyMhA37598eKLL+LZZ5/VO5xCEdEVpdQeAIEAipecCxIfH397rGfDhg3639HNyODBzitX+Nra3R0Pubvj4sWLtw+pW7cuDh06hGqFdW90cOrUKTRo0AAAsGnTJjQu6qpENrJ161Z8/vnn2Lt3L1xdXc1+nrs7X+p36QK89x5/LTmZE/OxYzyjPjKSE/ihQ7wGUKNGwMmTfOzDD3PibNeOkyYRLx/h6soLDJpMvIyEkxMvTwHwVPSzZwfhiy98UK4crykUG8vnqlqV20hNBfbv515wVhYn3thYfpNo04aXpfD15Ye7u+WvX+5p3kRU/GneAgBfiQwdOhRNmjTBu+++q3c4+bp06RJKlSoFd3d3KKXKAegM4PP8jreoznnw4MEUGRkJpRTq1q2L+fPn307Wuhg/HvjhB+D993nZtTwYOTn37dsXJ0+ehJOTE+rUqYN58+bB29tb77DuU79+faSlpaFq9gIPbdu2xbx58zQ9R2bmnd7s2bM8GebaNf7Y1RW4eJHrhb29ueQvJoY/Dwjgz8PDAVfXDGRkHENm5g0ADVGmDNCzZ3UkJfEQTNOm3Ft3d+dHzrpAXl73rZipuS1btuCdd97B2bNnMX78eHz88cfWPWExDRgwAHv27EFiYiKqV6+OSZMmYejQoXqHdZfw8HC0a9cODz30EJyyJ5dNmTIF3Qy2oNnRo0fx8ssvIzIy8hi4jPl7Ivo0v+MdZxJKeDiXBgweDHz3nd7RCJ3krGqvFP9r9Img9rDEpdCcWTcYHKOU7upV7inXrcu1XKLEUurOuLTcYxP2zDGS8xtv8PVveDjfNRJCCDtn/8l5xQp+fPopV/sLIYQDMPiIXCGiooCRI7lw9t7dJIQQwo7Zb3LOyuL5wwAv72Xt2+tCCGFD9pvRvvkG2L4dWLyYbwQKIYQDsc+e85kzvP5ix453pp8JIYQDsb/kTMRbRbi4AHPnSr2UEMIuHDx4EM2aNYNSqqxSqnz2ms75Tqu2v+S8aBGwaxfwxReAj4/e0QghhFlatWqFXr16AcBkANMBLM9v0SPA3mYIxsXxnFt/f07QRp/+JUQhZIZgyZKeno4yZcocBXALwKNElJXfsfaT3Yi4bC4tjZcyk8QshLAzly9fBgA3ABUAFLiZqf1kuLVrgY0bebJJ9sptQghhT4YPHw4A4wCsQAEr0gH2Ukp3+TLwzjtAy5bA6NF6RyOEEEW2dOlSuLi4gIhWKqWcAfymlOpIRLvyOt4+xpxHjOAtNVau5PFmIRyEjDmXSGaVmBl/WOPwYSA0lPc+ksQshCghjJ2cTSbgzTcBDw9g4kS9oxFCCJsx9pjzkiW8xfF339n/ds1CCFEExu05X7kCjB3Le8sPHqx3NEIUaOLEifD29oa/vz/8/f2xZcsWvUMSds64PecJE3g3z23bpKZZ2IXRo0djzJgxeochHIQxs97Ro0BICFdpPPyw3tEIIYTNGS85E/FEk8qVgcmT9Y5GCLOFhISgWbNmCAoKQnJycr7HhYaGIiAgAAEBAbh06ZINIxT2xHh1zj/8ADz3HLBwITBkiObNC1FcnTt3xoULF+77enBwMNq2bYtq1apBKYVx48YhPj4eixYtKrRNqXMukcyqczZWcs7IAPz8eIz52DHZ3UTYpejoaPTo0QN//ZXvgmO3SXIukcxKzsbKft9+C/z7L6+hIYlZ2JH4+HjUrFkTALBhwwb4+eW7TK8QZjFOBkxJ4Ykm7doBPXvqHY0QRfL+++8jMjISSinUrVsX8+fP1zskYeeMk5xnzAASEnjMWXY3EXZm2bJleocgHIwxqjUuXOCdTZ57DmjbVu9ohBBCd8ZIzpMm8SL6U6boHYkQQhiC/sn55Ene2WTECFlEXwghsumfnD/8EChXDhg3Tu9IhBDCMPRNzr/+CmzYwAsceXrqGooQQhiJfsmZCJg1C+jYUbaeEkKIe+iXnHftAlavBvr0AcqX1y0MIYQwIn2SMxEvCertDbz6qi4hCCGEkekzCWXHDh5vnj0bKFtWlxCEEMLIbN9zzuk116oFDB1q89MLIYQ9sH3Peds2YP9+YO5coEwZm59eCCHsgW17zjm95tq1gaAgm55aCCHsiW17zj/9BEREAKGhQOnSNj21EELYE9v1nImAVauA+vVlhxMhhCiE7XrOO3cCy5dzr7lUKZudVggh7JHtes5TpgA1awIvvWSzUwohhL2yTc/599+B3bt5QX2p0BBCiELZpuc8dSpQpQowfLhNTieEEPbO+sn5r7+ATZuAUaMANzern04IIRyB9ZPztGm8sNFbb1n9VEII4Sism5zPnOHyuREjeFhDCCGEWaybnBcvBvz9gXfftepphBDC0VgvOSclcXWGvz/g5WW10wghhCOyXnKePx+4eVN2ORFCiGKwTnJOS+MtqLp0Afz8rHIKIYRwZNaZhLJ6NXDhAvDdd1ZpXgghHJ32PWci4MsvgaZNuecshBCiyLTvOe/ZA/z5J7BwIaCU5s0LIURJoH3P+csvAQ8P4MUXNW9aCL2tWbMGvr6+cHJywqFDh+763tSpU1G/fn00atQI27Zt0ylC4Si0Tc4nTwKbNwNvvCEbtwqH5Ofnh/Xr16N9+/Z3ff3EiRMICwvD8ePHsXXrVowcORJZWVk6RSkcgbbJeeZMXnXu9dc1bVYIo2jSpAkaNWp039c3btyI/v37o0yZMnjggQdQv359RERE6BChcBTaJefERGDJEmDwYMDTU7NmhbAHsbGxqFWr1u3PfXx8EBsbq2NEwt5pd0MwNJQnnbzzjmZNCqGHzp0748KFC/d9PTg4GL17987zOUR039dUPjfEQ0NDERoaCgC4dOmSBZEKR6ZNcs7I4HU0hg0DfH01aVIIvezYsaPIz/Hx8UFMTMztz8+fPw+vfJYtGD58OIZnr20eEBBQvCCFw9NmWGPTJuB//wN69tSkOSHsTa9evRAWFoa0tDRERUXh1KlTaN26td5hCTumTXKeMweoUwfo1k2T5oQwqg0bNsDHxwf79+9H9+7d0bVrVwCAr68v+vXrh6ZNmyIwMBCzZ8+Gs7OzztEKe6byGisrAsLff/NswKlTgQ8+0CouIUqEgICA++qlhcMza3ae5T3nuXOB0qWBoUMtbkoIIQSzLDmnpHD5XL9+PCtQCCGEJixLzitWANeu8YxAIYQQmrEsOc+ZAzz8MNCmjUbhCCGEACytcz56FFiwQFafE0IIjVnWc65UCRg4UKNQhBBC5LAsOb/yCuDqqlEoQgghcliWnD/9VKMwhCiZqlWrpncIwqAsn4QihBCiKGw0CUUIIYTmJDkLIYQBSXIWQggDkuQshBAGJMlZCCEMSJKzEEIYkCRnIYQwIEnOQghhQJZu8CorHgkhhBVIz1kIIQxIkrMQQhiQJGchhDAgSc5CCGFAkpyFEMKAJDkLIYQB/T8iqDnx5f5BLgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "thisPlot = sp.plot(quad(x).subs(a, 1).subs(b, 2).subs(c, -1),\n",
    "                   intQuad.subs(a, 1).subs(b, 2).subs(c, -1).subs(x, x),\n",
    "                   (x, -5, 3), show = False, title = \"Function and integral\")\n",
    "thisPlot[0].line_color = 'blue'\n",
    "thisPlot[1].line_color = 'red'\n",
    "thisPlot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving differential equations\n",
    "\n",
    "A range of ordinary differential equations can be solved with sympy using `dsolve`. We will look at the equation:\n",
    "\n",
    "$$\n",
    "\\frac {dg}{dt} = \\exp(-t) - g,\n",
    "$$ \n",
    "\n",
    "with the inital condition $g(0) = 1$. In order to set up this, or any, differential equation in sympy, we first have to define the function we want to solve for, in this case $g$. As we don't know what this will be until we have solved the equation, we can't write down an explicit form for it. Sympy allows us to define a general function using `sp.Function` as below. When defined like this, $g$ can be a function of anything and can have any functional form. We have to make sure we refer to it consistently; in this example, we must refer to it explicitly as a function of $t$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of differential equation is Eq(g(t), (C1 + t)*exp(-t))\n"
     ]
    }
   ],
   "source": [
    "t = sp.symbols(\"t\")\n",
    "g = sp.Function(\"g\")\n",
    "diffEq = sp.Eq(g(t).diff(t), sp.exp(-t) - g(t))\n",
    "solDiffEq = sp.dsolve(diffEq, g(t))\n",
    "print(\"Solution of differential equation is\",solDiffEq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how we have defined the derivative w.r.t. $t$ in the above: `g(t).diff(t)`.\n",
    "\n",
    "An alternative method of solving the equation requires that we rewrite it so that the RHS is zero:\n",
    "\n",
    "$$\n",
    "\\frac {dg}{dt} - \\exp(-t) + g = 0.\n",
    "$$ \n",
    "\n",
    "We can then solve it as below. Note also the alternative way of defining $\\frac{d}{dt} g(t)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of differential equation is Eq(g(t), (C1 + t)*exp(-t))\n"
     ]
    }
   ],
   "source": [
    "g_ = sp.Derivative(g(t), t)\n",
    "print(\"Solution of differential equation is\",sp.dsolve(g_ - sp.exp(-t) + g(t), g(t)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which ever method we choose, we now have to find the arbitrary constant `C1`. This can be done using symbolic methods as below. Note that the statement `Eq(1, C1)` is an equation that tells us $1 = C1$, as described above!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equation for C1 is Eq(1, C1)\n"
     ]
    }
   ],
   "source": [
    "equationForC1 = solDiffEq.subs([(t, 0), (g(0), 1)])\n",
    "print(\"Equation for C1 is\",equationForC1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve symbolically the equation that describes radioactive decay,\n",
    "\n",
    "$$\n",
    "\\frac {d N}{d t} = -\\frac{N}{\\tau},\n",
    "$$ \n",
    "\n",
    "using the fact that at time $t = 0$ the number of nuclei is $N_0$.\n",
    "\n",
    "(This equation describes the situation in which the number of decays per unit time in a sample is proportional to the number of remaining nuclei, $N$. Hence the rate of change in $N$ with time, $\\frac {d N}{d t}$ is proportional to $-N$.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of decay equation is Eq(N(t), C1*exp(-t/tau))\n",
      "Equation for C1 is Eq(N0, C1)\n"
     ]
    }
   ],
   "source": [
    "t, tau, N0 = sp.symbols('t tau N0')\n",
    "N = sp.Function(\"N\")\n",
    "decayEq = sp.Eq(N(t).diff(t), -N(t)/tau)\n",
    "solDecayEq = sp.dsolve(decayEq, N(t))\n",
    "print(\"Solution of decay equation is\",solDecayEq)\n",
    "equationForC1 = solDecayEq.subs([(t, 0), (N(0), N0)])\n",
    "print(\"Equation for C1 is\",equationForC1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence result is $N(t) = N_0 \\exp{\\left[-\\frac{t}{\\tau}\\right]}$, as expected!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Series\n",
    "\n",
    "Sympy knows about Taylor series. To use them, we must specify the function we want to express as a series and the variable in which we want the expansion to be done. We can also specify the point around which the expansion should be made ($\\pi/4$ in the example below) and the maximum power in the expansion (the 3 below means we get terms up to the power 2). Look at [the documentation](https://docs.sympy.org/latest/modules/series/index.html) for more information! The example is the determination of the Taylor expansion for the sine function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First terms in Taylor series for sine function about pi/4 are sqrt(2)/2 + sqrt(2)*(x - pi/4)/2 - sqrt(2)*(x - pi/4)**2/4 + O((x - pi/4)**3, (x, pi/4))\n"
     ]
    }
   ],
   "source": [
    "sinSeries = sp.sin(x).series(x, sp.pi/4, 3)\n",
    "print(\"First terms in Taylor series for sine function about pi/4 are\",sinSeries)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using numpy and matplotlib, plot a graph showing a cosine function in the $\\theta$ range $0...\\pi$. Superimpose on this a plot of the sum of the Taylor series for the cosine up to and including the $\\theta^3$ term, expanded about $\\pi/2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cosSeries\n",
      " pi/2 + (x - pi/2)**3/6 - x + O((x - pi/2)**5, (x, pi/2))\n"
     ]
    }
   ],
   "source": [
    "# <-- Demo -->\n",
    "#\n",
    "cosSeries = sp.cos(x).series(x, sp.pi/2, 5)\n",
    "print(\"cosSeries\\n\",cosSeries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVkAAAEWCAYAAADM/ORiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gU5drH8e+dAkFBWkCp0pUeSGiiAWlWQKrAkaIiNvRgfW0HFMEuKooCooDl0KyIoFJElCagFOlFEaSHIh1D7vePmXCWkLIh2Z1Ncn+ua67s9N8O5M7sszPziKpijDEmMMK8DmCMMbmZFVljjAkgK7LGGBNAVmSNMSaArMgaY0wAWZE1xpgAsiJrgk5ExonIkCDv8y4RmRXMffpDRKqJyMEAbLeriPwlIkdEpHp2b9/4z4psLiAiPURkqfsLtVNEZojIlV7nyiwRecJ9D0dE5ISInPYZX+11vkBQ1Q2qWiQAm34NuE1VC6rq2gBs3/jJimwOJyIPAq8DzwEXA+WBt4H2XuY6H6r6nFsUCgJ3AQuTx1W1pheZRCRMRALyeyIiEQHabiRQGjivP0wiEp69ifI2K7I5mIgUBgYD96rqZ6p6VFX/UdWvVPURd5n8IvK6iOxwh9dFJL87L1pEponIQRHZLyI/plVQROQNEdkmIn+LyDIRucpn3tMiMllEPhCRwyKyWkTifObXE5Ff3HmTgKgsvOd3RGS7m+NnEWnsTi8vIkdF5CKfZZu67/mc9yQizdxMh0RkkYg08Jm3SEQGi8hi4BhOwUq5/n/cTw1/i8ja5OMhIuHuvC0isk9EPhaRIu68y0UkUUTuEJFtwPTkaT7bLeYex13u8R6UnN9d9ic3814R+SCVXIWBA+7o+uRPACJS2/33PSgiK0XkOp91JorIcBH5TkSOAk1S2W60T64D7r9j8rx7RWSziCSIyGcicrHPsXjLzXpIRFaIyGWp/8vmYqpqQw4dgGuBRCAinWUGA4uAkkAJYAHwrDvveWAkEOkOVwGSxnZuAYoDEcBDwC4gyp33NHACuB4Id7e7yJ2XD9gKPODuozPwDzAkg/fWB/gplem9gKLutp4EtgGR7rw5wK0+y74DvOy+vguY5b4uCfwNdHXfTx9gL1DYnb8I2AJc5u4nIkWGuu78iwEBKgEV3XmPAT/iFOYoYBww1p13OaDAGOACoIA7LdFn2zOAN935pYBfgd7uvM+Bh919FgCapnHsotz9lPUZ3+r+u0UC1wBHfDJPBPYDjXBOvPKnss3ZwIdAEfffNN6dfr37f6GOu5/RwEx3XntgIXCRu92aQEmvf2+C/nvqdQAbsvCPB/8CdmWwzGbgep/xa4A/3NeDgS+BKuex7wNAXff108kFzB2vARx3X8cDO/Ap3jiF/ryKbIplBOdM8zJ3vDcw232dD0gA6rjjvkX2DmBeim39CnRzXy8CnkhnvzWBncDVnFuAf/ctfkBFN6P4FNnSPvPPFFngUuAo7h8Nd9qtwAz39WTgLaBUBsclZZFt7RZZ33+Dz4HH3NcTgdHpbK8icAoolMq8j4HBPuNFgCTgErcArwYaksYf77wwWHNBzpYARGfQtlca5xcs2Vb+9/H3ZWAT8J378faxtDYiIg+5H4sPifNteGEg2meRXT6vjwFRbq7SwF/q/gb6ZDgvIvK4iKwXkUM4hT7KJ8enQAMRKYPzC75dVVemspmUxyQ5Uxmf8W1pZVDV1ThnrEOBPW6TwMUiIkA5nGaAg+5x+hXnLK64u3qSqu5IY9OXuu9nr8/6b+CcMYPzaeAC4Ff3I/8taWVM5f3+mcq/gV/v131Pe1T1cBrbPnMsVfUgzqeEMjhn5e8Bo4DdIvK2iBT0M3OuYUU2Z1uI8zH9pnSW2YHzy5usvDsNVT2sqg+paiWgLfCgiLRMuQG3vfH/cD5eF1Xn2/BDOGdnGdkJlHELkG+GTBOR1sB9QAecM6ZiwPHkHKp6BOcMrQfQE+fjbWpSHpPkTH/5jKf7eDpVHa+qV+A0FUThnJmru40WqlrEZ4hS1X1+bHcbzsf4oj7rXqSq9d19/qWqt+E0I9wPvC8i/hzLHZx7zDPzfrcBJdMokGcdS7dN+CLcP6yqOkxV6+E0J9QF/u1H3lzFimwOpqqHgIHACBG5SUQuEJFIEblORF5yF5sAPCUiJUQk2l3+IwARuVFEqrgF8G/gtDukVAin7XcvECEiA3F+kfyx0F33fhGJEJGOOB8fz0chnPbcvTjNAYM590u0D4C+OO3VH6exnalAPRHp7GbqhVN0vvEnhIjUcL84y49T5I/zv+M2EnhBRMq5y5YUkbb+bFdVf8dpqnhJRAqJc2VDVXEvxxORm0WktFvMk6+tTUxrez5+BMJEZID7flsDbYApmcg1D3hLRAqLSD4RiXdnTwDuEJFaIhIFvAjMUdVdItJYROLcTzRHcZocUvv/latZkc3hVHUY8CDwFE7x2Qb0B75wFxkCLAVWAquAX9xpAFWBWThnTwuBt1V1biq7+Rbno98GnI+GJ0j/46VvvlNAR5w21gPAzcBn/r/Ds3yF88u+GeeLp30479nX9zhfCv2kqjvTyLQbaIfzxVkCzvG60f2o648CwKvu/ncCBXH+eAG8hHNM54jIYZz25/p+bhegO85Z+jqcL6Mm8b/mgibAMhE5glMg+6XT9HCGqp4AbsT50jEBGAbcrKqbM5krEtiI0zR0t7vtaThfdE7FOau9BOdTBO77GIfzB2ELzv+d4ZnYZ64gZzfTGJPzicgCnD8YH3mdxRg7kzW5iog0BarhfAlmjOcCcseJMV4QkYk4l6jdq6rHvc5jDFhzgTHGBJQ1FxhjTADlqeaC6OhorVChQqbW2bh/I1WLVQ1MoGxg+bIu1DNavqwLdMZly5btU9USqc70+pazYA6xsbGaWbGjMr9OMFm+rAv1jJYv6wKdEViqdlutMcYEnxVZY4wJICuyxhgTQHnqiy9jcqJCEYX4/fffOXHihNdRUvVSzEusXRvaPdxkV8aoqCjKli1LZGSk3+tYkTUmxPWp1IdChQpRoUIFzn6YWWjQvUr1EqHdV2N2ZFRVEhIS2L59OxUrVvR7PU+bC0TkfRHZIyK/pTFf3G4xNrnPz6zvM6+3iGx0h97BS21McJW9oCzFixcPyQKbl4gIxYsXz/QnCq/bZMfhPJIuLdfhPCmqKtAPpzsRRKQYMAinu4yGwCARKRrQpMZ4RBArsCHifP4dPG0uUNV5IlIhnUXaAx+416EtEpEiIlIKaI7Tj9B+ABGZiVOsJ2RXtq1b4b33YMfyvjyfAAUKQOHCUKQIlCwJpUpB6dIQdd5dAhpj8oJQb5Mtw9nPLd3uTktr+jlEpB/OWTBRxaOIGx2X2mLnOLyhPhteHQ3cxRNfpb1cZJE95I/eTlSpPyhQZhMFymzkgvLrCI8KzvNJ1u5b6/d78kKo54PQzzi0zlDW7F3jdYw0nUg8kW6+35b/xtTJU3niuSeCmOpsGWXMjF2Hd9FzdM+MF3SFepFN7dxc05l+7kTV0Tg9aBIXF6dL+y31e+f6MsSNasj8Pj9z7BgcOgQHD8Lu3bBzJ2zbBlu2lGTTppKsXlWfbT8664WFQc2acOWV0Lw5NGsGF1+c7q7OW9zoODLznoIt1PNB6Gec8/McapSo4XWMNK3ZuybdfDVa16Br665BTHSujDJmhuyTc/6/yJ1pNyOEepHdjtOJW7KyOE9f347TZOA7fW5271wEJCyJqCinWaBYsbSXVYUdO2D5cvj5Z1i8GD76CN55x5lfrx5cfz3ceCM0bOgUYmNyig8++IBXXnkFEaFOnToMGTKE2267jb1793JBkQuY/NFkypcvz5QpU3jmmWcIDw+ncOHCzJs3j7lz5/LKK68wbdo0nn76af7880+2bNnCn3/+yYABA7j//vsB+Oijjxg+fDinTp2iUaNGvP3224SHh3v8zrMu1IvsVKC/+5zQRsAhVd0pIt8Cz/l82dUGeNyrkOAU5DJlnOGGG5xpiYnw668wezZMnw4vvABDh0LZstCpE3Tv7hRc+07D+GvAAOcPeXaKiYHXX097/urVqxk6dCjz588nOjqa/fv307t3b3r16kXv3r0Z8sYQ7r//fr744gsGDx7Mt99+S5kyZTh4MPXefNatW8f333/P4cOHueyyy7j77rvZtGkTkyZNYv78+URGRnLPPffw8ccf06tXr+x9sx7wtMiKyAScM9JoEdmOc8VAJICqjgSm43TtvAmnm+lb3Xn7ReRZYIm7qcHJX4KFkogIaNDAGR57DA4cgGnT4NNPYeRIeOMNqFYNevWCW291vkgzJtTMmTOHzp07Ex3t9LxerFgxFi5cyGefOV21te3Sltefdap006ZN6dOnD127dqVjx46pbu+GG24gf/785M+fn5IlS7J7925mz57NsmXLaNCgAQDHjx+nZMmSQXh3gef11QXdM5ivwL1pzHsfeD8QudJ09Chs2QK1azvjL70EP/4I+/bB8eNw8iSUKwfffefMf+IJ5zKFiy+GihUpWrkyPa+qQc+eFTh0CD75BD78EJ56CgYNgrZtoX9/aNHCzm5N6tI74wwUVc3w0qXk+SNHjmTx4sV8/fXXxMTEsDyV0+78+fOfeR0eHk5iYiKqSu/evXn++eezN3wIsJbB9OzcyXWLEuD226FGDShUyPkmK9nmzbB9O1x4IVSo4BTfOnX+N/+vv2DRIue09f77nXaE7s7flcKF4fbTo5n76HQ2/3KIhx6C+fOhVStnE++959RsY7zWsmVLJk+eTEJCAgD79+/niiuuYOLEiQBM+3QaV155JQCbN2+mUaNGDB48mOjoaLZt86tTY1q2bMknn3zCnj17zuxj69atAXg3wRfqbbLeGjaMZ8f+AUU/h6ZNoWtX5xssVedUc9So9NcfP975qepckrB5s/ManAbbRx6Bv/+mUng4LzZtypD7rmNqvs4M/m8V+vaFgQPhoYegXz8oWDCg79SYNNWsWZMnn3ySZs2aER4eTr169Rg+fDi33XYbL7/88pkvvgAeeeQRNm7ciKrSsmVL6tatyw8//JDhPmrUqMGQIUNo06YNSUlJREZGMmLECC699NJAv73AS+tBs7lxyPRDuzdv1h5PVldNTMzcev46elR1zhzVJ55QjYlRBdXBgzUpSfW7aSe1a+OtCqolSqi+8oqzeEqh/sDkUM+nGvoZZy+e7XWEdK3es9rrCBnKzoxr1qw5Zxr20O7zVKkS68tfAIG6jOSCC+Dqq51LDn791bnw9q67EIHWJ75i0uIKHGjQhvtKfcL/PZxI5cpOy8M//wQmjjEm+1mRDSVly0IJt5ughg1h4ECK7F7Pf1Z24WjJijyT/zkG3H2CWrXgiy/+1/JgjAldVmRDVbly8PTTztUMX35J/jqXc0eBD5n8eT7Cw6Frh1O0bg3H/6rsdVJjTDqsyIa68HBo1w5mzkR+/pl2N4Wxcv5hEopWodNPD7Dv2WE88AAcPux1UGNMaqzI5iSFCgEQkXiCQu1bctc/w9kilSn6+iDqVzvCJ59YE4IxocaKbE5UogSMHYusXcv8ehcwkMEs3FeF/l120bGj8wwFY0xosCKbk1WrxhP9KsHChRR7oA8PvXQJ33wDzS/fxdixdlZrskdCQgIxMTHExMRwySWXUKZMmTPjp06d8ns7iYmJFClSJIBJ0/f555/z8ssvB32/djNCbtC4MWGNG/MI0CV2C6Va12Tcbb3oPulFXhtbhFKlvA5ocrLixYufuT326aefpmDBgjz88MMB329iYiIREdlTohITE+nQoUO2bCuz7Ew2l6nQoAT5BtzLHTKG176tzhPVpvD5516nMrlV27Zt6dKqCzVr1mTMmDEAjBo1ikceeeTMMu+88w6PPvroWeslJSXx4IMPUqtWLWrXrs0nn3wCwKxZs2jVqhXdunWjXr16Z62TmJhIz549qV27NrVq1WL48OEAbNy4kWuuuYbY2Fji4+PZsGEDALfccgsPPfQQV199Na8PfZ0xY8YwYMAAAHbv3k3Hjh2Ji4ujYcOGLFq0CHAehlO3bl1iYmKoX78+R48ezfIxsjPZ3KZQIeTVV5B/9aBor36MXd2VSR27csdtE3jtjTC7PTc38H1+RrKuXeGee+DYMefBxSn16eMM+/ZB585nz5s797yjjB8/nl2nd1HhwgrExcXRqVMnevToQUxMDM8//zwRERGMHTuWcePGnbXelClTWLNmDStWrGDv3r00aNCA+Ph4ABYtWsSaNWsoX778WessW7aMffv2sWrVKoAzj1Ls168fY8aMoXLlysyfP5/+/fvznfuQps2bNzN79mzWJaxjwZcLzmzr/vvv59FHH6Vx48b88ccf3Hjjjfz222+8/PLLjB49mkaNGnHkyBGisqF/KSuyuVX9+kQtX0Tii69wyVcHeG9sGD/Oh8mTz36GjTFZ8dprrzH5s8lERUSxfft2Nm/eTFxcHPHx8cyYMYNKlSoRHh5OjRo1SExMPLPeTz/9RI8ePQgPD+eSSy7hyiuvZOnSpeTLl48mTZqcU2ABqlSpwvr16/n3v//N9ddfT5s2bTh48CCLFi2iU6dOZ5bz3U+XLl0IS+UJ+bNmzWL9+vVnxg8cOMDx48dp2rQpAwYMoEePHnTq1ImC2XBWYkU2N4uIIOLJx2j2JMyaA693mc/SeuNYMuw1bru/oD1OMadK78zzggvSnx8dnaUzV1+zZs1i3rx5TJgxgfrl63PllVee6S67b9++DBs2jAoVKnDrrbees66m863shRdemOr04sWLs3LlSmbMmMHw4cP59NNPefHFF4mOjk71kYrpbUtV+fnnn8mXL99Z05966inatWvH119/TYMGDZg7dy5Vq1ZNM6s/rE02j2jRAj5+YCl9kt6j2YAYBt6wjGxobjJ52KFDhyhWrBhRBaJYvXo1S5YsOTOvadOmbN68mSlTpnDzzTefs258fDwTJ07k9OnT7N69m/nz5xMXl35nlnv37kVV6dKlC8888wy//PILRYsWpVSpUnzufvGQlJTEihUrMszeqlUrRowYcWY8uUhv3ryZOnXq8Pjjj1OvXr2zznbPl6dFVkSuFZH1IrJJRB5LZf5rIrLcHTaIyEGfead95k0NbvKcqdBT/4a5P1Ci8Cn+M6MJw6sMZ/06u87LnJ8bbriBY8eO0aF5BwYPHkyjRo3Omt+5c2fi4+MpXLjwOet27tyZyy+/nLp169KqVSuGDRuWYU8I27ZtIz4+npiYGO644w6ee+45ACZOnMjIkSOpW7cuNWvWZNq0aRlmHzFiBPPnz6dOnTrUqFGDd999F4BXXnmFWrVqUadOHYoUKUKbNm38PRxpS+vxXIEegHBgM1AJyAesAGqks/x9wPs+40cyu89MP+pQQ/8xeOeVb98+3dO4rSrorQUm6NSp2Z8rWagfP9XQz5hTH3V4zTXX6Ny5c4OcJnV59VGHDYFNqrpFVU8BE4H26SzfHZgQlGS5XfHilFjwJXvf+YSVl3elXTsYOugUSUleBzO5QUJCAtWqVaNo0aI0a9bM6zieE/XotiAR6Qxcq6p93fGeQCNV7Z/KspcCi4CyqnranZYILAcSgRdU9Ys09tMP6AcQVTwqtuZzNTOVc+2+tVSPrp6pdYIpq/mSTuXnn/f68dXyQQyp2INfHvyBsHzZ1+9NqB8/CP2MQ+sMpVzlcl7HSNOJxBNERWT9UqdAys6Mu37fxaPLz77ud9mdy5apauqNymmd4gZ6ALoAY3zGewJvprHs/6WcB5R2f1YC/gAqZ7RPay5IXdLOXbqt4pWqoGMvflS3/ZF9PUGE+vFTDf2MsxbP0qSkJK9jpCkv9YyQlJSUo5oLtgO+f57LAmk92qQbKZoKVHWH+3MLMBeod+5qxh9yycWUXTebrdfdRZ/dL7HhshtZ/sMhr2MZ1/Zj20lISEj3sicTeKpKQkJCpm9Q8PI62SVAVRGpCPyFU0h7pFxIRC4DigILfaYVBY6p6kkRiQaaAi8FJXVulS8fl05/h78GxnDVs/0Z3WoQ2z57nbZtvQ5mxm0ZR/yl8ezdu9frKKnadXgXsi+0L7rOroxRUVGULVs2U+t4VmRVNVFE+gPf4lxp8L6qrhaRwTin3smXZXUHJurZf8arA6NEJAnnMrQXVHVNMPPnVmUG38m++jWZPLguP7aHN99I4t777HJqLx1OPEzFihW9jpGmnqN7srTfUq9jpMvLjJ7e8aWq04HpKaYNTDH+dCrrLQBqBzRcHhZ905VMbw23dT1CzP3X8PHsAXT/rAup3J1ojMmA/dqYVF14IXw89hQXXyL868uuTIx7hVMnrU3QmMyyImvSFFGyGJW3zGJtrS70+PURpl32IEf+totpjckMK7ImXVIgiuorJrKq1QA6bn2dzy5/nH37vE5lTM5hRdZkLCyM2t8NY2W/txiacDfNmlk/Ysb4y4qs8Y8IdUbdy6hvK7BtaxITaz7LH8sSvE5lTMizImsypXlzWDD6N+45OJSTjePZ+IOd0hqTHiuyJtNq9ajDzve/oczpP4loEc+aGVu9jmRMyLIia85LxVubs3/SLIpqAoVvuJKVn2zwOpIxIcmKrDlv5bs04ui0uWh4OE/0/ovFi71OZEzosSJrsqTM9XVJWrOedaWupnVrWPztwYxXMiYPsSJrsqx81fzMnQt3XPAxla+tyvLxGfexZExeYUXWZIuyZeHRzxrzT0QU5W5tyYoPrNAaA1ZkTTa6+IrKyNy5nAovQNk+LVn50UqvIxnjOSuyJltd0rQyzPmeU2EFKN2rJfnWVPA6kjGesiJrsl2pq6qgs+fwZtFB/Pzue6ywlgOTh1mRNQFRullV+iztT1j+4zzYfBkbZ2zyOpIxnvC0yIrItSKyXkQ2ichjqczvIyJ7RWS5O/T1mddbRDa6Q+/gJjf+qFgRqt9/J2MO30yBti3Z+uOfXkcyJug8K7IiEg6MAK4DagDdRaRGKotOUtUYdxjjrlsMGAQ0AhoCg9x+v0yIiSz9F6f/O5lCSYdIatGSHct2eh3JmKDy8ky2IbBJVbeo6ilgItDez3WvAWaq6n5VPQDMBK4NUE6TRVW61mfnezMokbiTI1e0Zu86e3qXyTu87OOrDLDNZ3w7zplpSp1EJB7YADygqtvSWLdMajsRkX5AP4Co4lHEjY7LVMi1+9Zmep1gykn5at4cz6hJc3i74aN8+Owmwgsc9TidIycdw1AU6vnA44yq6skAdAHG+Iz3BN5MsUxxIL/7+i5gjvv6EeApn+X+AzyU0T5jY2M1s2JHZX6dYMpp+ea/tlgLhJ/Uq65SPXrUo1Ap5LRjGGpCPZ9q4DPi9LCdat3xsrlgO1DOZ7wscNbDSVU1QVVPuqPvArH+rmtC0xUDGjLuv/lY/+Mevq7zGP8cT/Q6kjEB5WWRXQJUFZGKIpIP6AZM9V1AREr5jLYD1rqvvwXaiEhR9wuvNu40kwN07QoT+nxDl80vsrDWHSSdtl5wTe7lWZusqiaKSH+c4hgOvK+qq0VkMM6p91TgfhFpByQC+4E+7rr7ReRZnEINMFhV9wf9TZjz1mJsL+Zt3UL8988wp2Fxrl76CiJepzIm+3n5xReqOh2YnmLaQJ/XjwOPp7Hu+8D7AQ1oAuqqWYP4sd4+WvzyKnNuuIQW0x/2OpIx2c7u+DKekTCh6dI3WFi+K2VmvMsHo457HcmYbOfpmawxYZHhxP72AT3aHuaLewtQrAzceKPXqYzJPnYmazyXr1B+xn4VTVydU+y5qR+/jV2S8UrG5BBWZE1IKFQIpn5wkNYyk4tvv4HfZ9oDZUzuYEXWhIyStUpy+utvCSOJsOuvYc9ve7yOZEyWWZE1IaVCm2rsencaJRJ3sqdxW47uCY1bb405X1ZkTcipeXtjfntiAqWPbuTJTutItJvCTA5mRdaEpIZD2/P5K1t446dY7rsP1G4KMzmUXcJlQtbtDxVhw24l8eVhzN0GV097yOtIxmSancmakPb889Cp3M9c/fXDLHhgitdxjMk0K7ImpIWFC/VXjmdloabUf70nK0cu8DqSMZliRdaEvKgiUZRd9iW7IstR+p72/DF7s9eRjPGbFVmTIxSrWhy+no4ivH3LfPbbM9dMDmFF1uQYFVpXZePXG3ljfy86doRTp7xOZEzGrMiaHOWK6wozdixE/fAN8+rehybZtV0mtHlaZEXkWhFZLyKbROSxVOY/KCJrRGSliMwWkUt95p0WkeXuMDXluib36tEDHmu+mFbr3uKH617wOo4x6fKsyIpIODACuA6oAXQXkRopFvsViFPVOsAnwEs+846raow7tAtKaBMyms0eyIJLu9P8uydY+KBd2mVCl5dnsg2BTaq6RVVPAROB9r4LqOr3qnrMHV2E02GiMUiYUH/5+6wsdAUxr/VizXh7PKIJTaIe3a8oIp2Ba1W1rzveE2ikqv3TWP4tYJeqDnHHE4HlOP1/vaCqX6SxXj+gH0BU8ajYms/VzFTOtfvWUj26eqbWCaa8nq/QjgsZ++zvfBrRjvHP/E6+YrszvY28fgyzKtTzQeAzLrtz2TJVjUt1Zlp9hQd6ALoAY3zGewJvprHsLThnsvl9ppV2f1YC/gAqZ7TP2NjM970e6n3KWz7VtfP2aKGCSRoTo3r4cObXt2OYNaGeTzXwGXE6f0217njZXLAdKOczXhbYkXIhEWkFPAm0U9WTydNVdYf7cwswF6gXyLAmdF1+VQkmTxFOrFjP3Lr/JikxyetIxpzhZZFdAlQVkYoikg/oBpx1lYCI1ANG4RTYPT7Ti4pIfvd1NNAUWBO05CbkXHstvN1pNjduGc68q570Oo4xZ/j1FC4RqYVzBUBU8jRV/SArO1bVRBHpD3wLhAPvq+pqERmMc+o9FXgZKAhMERGAP9W5kqA6MEpEknD+ULygqlZk87jmk+5mXu1VNF/0Aj/1q86Vo3t5HcmYjIusiAwCmuMU2ek4l1z9BGSpyAKo6nR3m77TBvq8bpXGeguA2lndv8ldJExosmQ4v5TdQIN372BVbBVq33mF17FMHudPc0FnoCXON/u3AnWB/AFNZcx5irwgkopLprAzsjzb73+RrVu9TmTyOn+K7HFVTQISRU+jVQwAABzKSURBVOQiYA/ON/rGhKSilYuR+M1s+kRNol07OHLE60QmL/OnyC4VkSLAu8Ay4Bfg54CmMiaLqrQoz4dToti26iBfNHzOrjgwnsmwyKrqPap6UFVHAq2B3m6zgTEhrU0bmHTzZ9yy9knmNR+Y8QrGBECGRVZEZie/VtU/VHWl7zRjQlmrj29l3mV9aT5/KAvum+B1HJMHpVlkRSRKRIoB0e51qcXcoQJQOlgBjckKCRMaLx3B8oviqffWbfaMAxN06Z3J3onTBns5TjvsMnf4EufpWcbkCPkK5qPsok9ICL+YE3f0Z8df9gxaEzxpFllVfUNVKwIPq2pFn6Guqr4VxIzGZFl09RIcnTKdLhFf0KGjcOKE14lMXuHP1QXvi8hTIjIaQESqisiNAc5lTLa7rEMNXvm4FMt+TmRs6/9arwomKPwqssApIPnWme3AkIAlMiaAOnSAKZ0mcfdP/+KHdq96HcfkAf4U2cqq+hLwD4CqHgckoKmMCaCbJvdgYdkuxH/9KEuemZ7xCsZkgT9F9pSIFAAUQEQqAyfTX8WY0CVhQp1lY1lfIIZqT3en1KpiXkcyuZg/RXYQ8A1QTkQ+BmYDjwY0lTEBdmHJC7lo9heckiiGjtzGgQS7I8wERoZP4VLVmSLyC9AYp5ng36q6L+DJjAmwMk3Ks2LE59zRXyjaI4yvv4YIvx7+aYz//H1odxRwAPgbqCEi8YGLZEzw1L37Cvb8axbffQdv9l7qdRyTC/nzPNkXgZuB1UDyZyoF5gUwlzFBE33llzyVUIu+/+3ITxeO48rRvb2OZHIRf85kbwIuU9UbVLWtO7TLjp2LyLUisl5ENonIY6nMzy8ik9z5i91bepPnPe5OXy8i12RHHpN39Z5yI78UbUHcu3fy23uLvY5jchF/iuwWIDK7dywi4Ti3516H0+tCdxGpkWKx24EDqloFeA140V23Bk6fYDWBa4G33e0Zc14iL4ikwuLJ7I4oQ4l+N7Fz6V9eRzK5RHoPiHlTRIYDx4DlIjJKRIYnD9mw74bAJlXdoqqngIlA+xTLtAfGu68/AVqK09lXe2Ciqp5U1d+BTe72jDlvxaoW59TkL7kg6Qj7m3fg+OFEryOZXCC9NtnkbwGWkaIXWdxrZrOoDLDNZ3w70CitZdyOFw8Bxd3pi1KsWya1nYhIP6AfQFTxKOJGx2Uq5Np9azO9TjBZvqxLmTHu2hYc/qYr37T+jgq3DkQ8vvUm1I9hqOcDjzOqaroDziVbGU7L7AB0Acb4jPcE3kyxzGqgrM/4ZpwiOwK4xWf6e0CnjPYZGxurmRU7KvPrBJPly7rUMg4erAqqbw3c7UGis4X6MQz1fKqBz4jTw3aqdcefNtnUvmrtc95V/X+2A+V8xssCO9JaRkQigMLAfj/XNea8PfUUDIr/nlsHV2DJ4BlexzE5WHptst1F5CugoohM9Rm+BxKyYd9LgKoiUlFE8uF8kZWyWWIq/yvynYE57l+NqUA39+qDikBVrN8xk41E4JEpDfkz6jKqDerOlhnrvY5kcqj02mQXADuBaMD3cUWHgZVZ3bE6baz9gW+BcOB9VV0tIoNxTr2n4jQDfCgim3DOYLu5664WkcnAGiARuFdVT2c1kzG+Lix5IQVnfcGpqxog7dtxaONiCl9axOtYJodJs8iq6lZgK9AkUDtX1enA9BTTBvq8PoHTdpvaukOBoYHKZgxA2aaXsmL4p9S4rwXL47pT/69phOezqwWN//y9rdaYPKtu/6tY+K8RzN5Xl8cf9zqNyWnscRjG+CH+o35MKgxvD4M6tf7hlluz/f4ck0tl6kzW7bW2TqDCGBPKXn8dbo9bwZW3VbNbb43fMiyyIjJXRC5yuwdfAYwVkWGBj2ZMaImMhBc/KkNYhBDdr4Pdemv84s+ZbGFV/RvoCIxV1VigVWBjGROail8WzcnJU7kw6TAHmt3E8f3HvY5kQpw/RTZCREoBXYFpAc5jTMir2qEWax7/iBrHlvJL/b7W661Jlz9FdjDOtaybVXWJiFQCNgY2ljGhrdFz7ZnbaggHth7ilSEnvI5jQpg/3c9MAab4jG8BOgUylDE5QbNvn+CWHklMeDqcy2KgXbY8ZdnkNv588VVWRD4XkT0isltEPhWRssEIZ0wokzBhzNhwrq39F0U6Xs2GT1d5HcmEIH+aC8biPCugNM7jBL9ypxmT5xUoAO+9B9V0PQW6tWPf2r1eRzIhxp8iW0JVx6pqojuMA0oEOJcxOUapuDIkjPmCEok72d6kM6eOnPI6kgkh/hTZfSJyi4iEu8MtZM9TuIzJNWre2pBl97xPzKF5LIq71644MGf4U2Rvw7l8axfOU7k6u9OMMT6ajujB3KZPUmL9T4x84aDXcUyI8Ofqgj8B+97UGD/Ezx3MLR3+j0n/KUSFenDddV4nMl7z5+qC8SJSxGe8qIi8H9hYxuRMYRFhjJ5QiAa1jrOrfT82TV3jdSTjMX+aC+qo6pnPPqp6AKgXuEjG5GwFC8KnoxO44fRUIju1JWH9Pq8jGQ/5U2TDRKRo8oj7oJgsPSJRRIqJyEwR2ej+LJrKMjEislBEVovIShG52WfeOBH5XUSWu0NMVvIYk93KNCrLntFfcnHiX2xv2NGuOMjD/CmyrwILRORZt2uYBcBLWdzvY8BsVa0KzHbHUzoG9FLVmsC1wOu+zRbAI6oa4w7Ls5jHmGxX6/ZG/NJ/LHX//pHF9e60Kw7yqAyLrKp+gHMb7W5gL9BRVT/M4n7bA+Pd1+OBm1LZ7wZV3ei+3gHswa7PNTnMFW92Z26zQdTa9AUjH9/qdRzjAXE6fw3yTkUOqqrvl2kHVPWcJgOf+Q1xinFNVU0SkXE4fY+dxD0TVtWTaazbD+gHEFU8KrbmczUzlXXtvrVUj66eqXWCyfJlXaAzapJycuQDrF5xC5XufISi9b/P1PqhfgxDPR8EPuOyO5ctU9W4VGeqakAGYBbwWypDe+BgimUPpLOdUsB6oHGKaQLkxym+A/3JFBsbq5kVOyrz6wST5cu6YGQ8dky1caMkfSLiRV3zwZJMrRvqxzDU86kGPiNOD9up1p2A9fGlqmk+2Nt90EwpVd3pPqt2TxrLXQR8DTylqot8tr3TfXlSRMYCD2djdGOyXYEC8OUHhzhV420i+rzOjmqLKd2onNexTBB41VvtVKC3+7o38GXKBUQkH/A58IE6j1v0nVfK/Sk47bm/BTStMdmgZLUiHJ8yjQuSjnC4+Y0c3nHY60gmCLwqsi8ArUVkI9DaHUdE4kRkjLtMVyAe6JPKpVofi8gqYBUQDQwJbnxjzk/VDrXYMGQKlU+sZl1MNxJPJHodyQSYJ12Cq2oC0DKV6UuBvu7rj4CP0li/RUADGhNAcU9ew7zVb3HFhP682WMBAz6NR8TrVCZQvDqTNSZPi//vXQy7bTUPfh7PG294ncYEkidnssYYePjdy1h8EGY/MI0muxNp9Pw5l4ubXMDOZI3xSFgYfDg+iaEFn6f2Cz1Y/f5iryOZALAia4yHLigYRqlFn7M3ohQX923Ln3O3eB3JZDMrssZ4rETNkpyeOp1wEjnd5jr2b7SOR3ITK7LGhIBK113Gn29OpdQ/W/mozQecOOF1IpNdrMgaEyLq3nslc19bzr//GEDPnpCU5HUikx2syBoTQq4dcDmvviqs/mQN3zYdbI9HzAXsEi5jQsyDD0KtyRNos2gI89rmdx6pZHIsO5M1JgS1+ukZfirfnfjpjxH/cRmv45gssCJrTAgKiwijwaqx/FKkBS/Om86yod94HcmcJyuyxoSo/Bflp/LKz1kdWY09g95iyc/WPpsTWZE1JoQVLncRAx4vyYAyn3D9DcKGDV4nMpllRdaYEHekzGG+mhlFET3Aunrd2PXzn15HMplgRdaYHKBaNfhixHbij33DsavacGDDXq8jGT95UmRFpJiIzBSRje7PVDtRFJHTPg/snuozvaKILHbXn+T2omBMrlbz5tpsef0rSp3ayq7613Fkx99eRzJ+8OpM9jFgtqpWxe1tNo3ljqtqjDu085n+IvCau/4B4PbAxjUmNNT/91Usf2IKVY6uYHPNdpw8eNzrSCYDXhXZ9ji9zOL+9PtBmm6/Xi2AT85nfWNyuiZDb2ThXR9Q6OCf3N9tD4nWg01IE6c32yDvVOSgqhbxGT+gquc0GYhIIrAcSAReUNUvRCQaWKSqVdxlygEzVLVWGvvqB/QDiCoeFVvzuZqZyhrqfcpbvqwL9Yxp5Tv4XSc2f/o4xRp9RaVez6Ae3b8Z6scPAp9x2Z3LlqlqXKoz0+orPKsDMAunF9mUQ3vgYIplD6SxjdLuz0rAH0BloASwyWeZcsAqfzLFxma+7/VQ71Pe8mVdqGdML9/gZ5J0FHfogsv6aFLi6SCm+p9QP36qgc8ILNU06k7AmgtUtZWq1kpl+BLY7dOtdylgTxrb2OH+3ALMBeoB+4AiIpL8d7sssCNQ78OYUPbUf4SKTcvQZP04FsbdZw+UCUFetclOBXq7r3sDX6ZcQESKikh+93U00BRY4/7V+B7onN76xuQFItBq3kBmxTzMFcvfZuEVD4EHTYAmbV4V2ReA1iKyEWjtjiMicSIyxl2mOrBURFbgFNUXVHWNO+//gAdFZBNQHHgvqOmNCSESJrRY+hKza9zHFYtf46cWA72OZHx40lSuqglAy1SmLwX6uq8XALXTWH8L0DCQGY3JScLCheYr3uDb2mG8NrcpLV+GRx7xOpUBu+PLmFwjPEJouep1itx8LY8+Ch8++Ks1HYQAK7LG5CIREfDRR/DU1fP512uxLLpmoBVaj1mRNSaXiYiAgTOa8H2F22g8cwiLWz9phdZDVmSNyYUi84cRv240Myv1o9Hs5/m52SNWaD1iRdaYXCoyfxhXrxvJt1XvpeGPr/Lf22d5HSlPsiJrTC4WESm0XP0mL7WYwb/Gtuapp+yENtisyBqTy0VECg/PvJa+feHrob+yqG4/9OQpr2PlGVZkjckDwsJg1Ch4otkCmqx6l9+q3sTpw8e8jpUnWJE1Jo8IC4PO39/L59eNpsa2b9lUsTWndh/wOlauZ0XWmDxEBDpMv4Ov/jWJiglL2FElnqN/WFc2gWRF1pg86KaPOjNzwHSWHKnONTcXYf9+rxPlXlZkjcmjbnitFeGfTmbpikjaN97N7i8Weh0pV7Iia0we1rEjfPMN3Pf7gxTucDV/vjzJ60i5jhVZY/K45s2h5uzhrMjXgPKPdmPzHS/YxbTZyIqsMYaa8cUp/dtMpl/UjcpjHmdjs9vhlF1Lmx2syBpjAChXNYor/vgvH1QcRNKP8xn6xFE7oc0GnhRZESkmIjNFZKP7M7Weaq8WkeU+wwkRucmdN05EfveZFxP8d2FM7lOkqNB9/dO80esXnnq1KD27nOD40tVex8rRvDqTfQyYrapVgdnu+FlU9XtVjVHVGKAFcAz4zmeRR5Lnq+ryoKQ2Jg+IjIQR4y7kpZcg5tP/IA0bsP/tiV7HyrG8KrLtgfHu6/HATRks3xmYoap2H6AxQSDidF9Te/zD/CKxFLu3OztvHgD//ON1tBxH1INGFxE5qKpFfMYPqOo5TQY+8+cAw1R1mjs+DmgCnMQ9E1bVk2ms2w/oBxBVPCq25nM1M5V17b61VI+unql1gsnyZV2oZ/Q6X+KfFej7amH6n3ifJSXK85+HC7OvSL6QyeePQGdcdueyZaoal+pMVQ3IAMwCfktlaA8cTLHsgXS2UwrYC0SmmCZAfpwz4YH+ZIqNjdXMih2V+XWCyfJlXahnDIV8CQmqQ2tP0L8opU923aDHj/9vXijky0igMwJLNY26E7DmAlVtpaq1Uhm+BHaLSCkA9+eedDbVFfhcVc98TlHVne57OwmMxXquNSagihWD//u1GyMf3szQyVVpeoWyd9iHdpmXH7xqk50K9HZf9wa+TGfZ7sAE3wk+BVpw2nN/C0BGY4yP8HAY/HIBvvwSLt7wIyUe6sXBGldw6a4TXkcLaV4V2ReA1iKyEWjtjiMicSIyJnkhEakAlAN+SLH+xyKyClgFRANDgpDZGAO0awdvrYzn4Uqfkbj5Dz58ZiP/vDECkpK8jhaSPCmyqpqgqi1Vtar7c787famq9vVZ7g9VLaOqSSnWb6Gqtd3mh1tU9Uiw34MxeVmlSjB0TQfevGMVPyS1IHJAfw616+l1rJAU4XUAY0zOlD8/PDO6FFUiLqfth534Y2Zp2rwDd93+DxIe5rQvGLut1hiTNUVifuTRDX051vx67rkHPqj5Iqdim8Byu0cIrMgaY7JBqVIwYwa8+SbM/PMyDq7cSlL9WLT/fXAgb3dxY0XWGJMtwsKgf394enUXejdcxwi9h6QRb3O6SjWYOtXreJ6xImuMyVZVqsC0+UVh+Js0jfqF+Qdr8fGcUiQmAkeP5rln1VqRNcZku/BwuO8+mLSuLi9e+z23vNGA2FjYffN90LCh07aQR4qtFVljTMBceilMmwaffQb798MjXzdnz9p9cP310KiRMyOXX19rRdYYE1Ai0KEDrF0L5Z/sReV/1tM/3ygSNu2HTp1gSO6+l8iKrDEmKAoWdOrpynX52NO+HyUPrKffRRP5iFucRyDMmwd33w0//5yrmhKsyBpjgqpiRZg8GRYsCmd9zM30HFSJqlVhwciV6PjxTjNCrVrw7LOwbp3XcbPMiqwxxhONGsHcuU6X5JdcAk0n9KdO9E7mdBvN6aLFYdAgaN36f2e1a9bAsQA/t3/r1mw/i7Yia4zxjAhccw0sWgRffw2Fyxem5cQ7KLVhHq8M2E7Cm/91FkpKghYtoEgRaNIEHnwQPvwQNm8+/50fOQLHjzuvv/vOeSBDhQqwaVO2vLdk9uwCY4znRJwLDq6/Hn78EV56CR59vTRPjihN167Q73blyjHvIT/96Czwzjtw4oRTbF991TnDbdMGSpaEEiWgUCGIioK2baFRI4of+gf69YO9e2H3btiyxfk5YQJ06+bcshYTAw88AEXT7KTlvFiRNcaElKuucoaNG53bdMePh48+Cqd69Rvo3fsG/jUJyl6SCBs2QIECzkoHDzo9QG7YAD/95Nz0cPy4UzAbNaLAydPOXWclSjjDjTc6d03UquWsX7u2czlZAFiRNcaEpKpVYfhweP5554uyd9+Fxx6Dxx+HZs0i6Ny5Bh06QGmA0qXh++/P3oDqmfbV7SWjYNeuoL8HsDZZY0yIu/BCuPVWWLDAObsdOND5pN+/P5Qt69xA9vTTTrvuWZ3pijgPVPCYJwlEpIuIrBaRJBFJvYdHZ7lrRWS9iGwSkcd8plcUkcUislFEJolIvrS2YYzJPapUcQrqmjWwejUMHgwREc7PJk2cvsiuu865Hve770LjAWBeNRf8BnQERqW1gIiEAyNwuqfZDiwRkamqugZ4EXhNVSeKyEjgduCdwMc2xoSKGjWc4amnICEB5sxxLglLviwsWdmy8HfhN3hgLVSu7FxEULq0811XdHTgny3uSZFV1bUATj+IaWoIbFLVLe6yE4H2IrIWaAH0cJcbDzyNFVlj8qzixaFLF2cAOHQIli51hlWr4LMfSjBq1P+u2PJ14YVQuDDky+d8d1amzLnNu1kRyl98lQG2+YxvBxoBxYGDqproM71MWhsRkX5AP4Co4lHEjU6zdSJVa/etzfQ6wWT5si7UM1q+LCgKxIPUWMvlxauT+HdxTu4rwz+HovnnUDSJRwpz+nhBTp64kBOnI9HTERxIOkLc6OezL4OqBmQAZuE0C6Qc2vssMxeIS2P9LsAYn/GewJtACZwz3OTp5YBV/mSKjY3VzIodlfl1gsnyZV2oZ7R8WRfojMBSTaPuBOxMVlVbZXET23EKaLKywA5gH1BERCLUOZtNnm6MMSHH++sb0rYEqOpeSZAP6AZMdf9qfA90dpfrDXzpUUZjjEmXV5dwdRCR7UAT4GsR+dadXlpEpgO4Z6n9gW+BtcBkVV3tbuL/gAdFZBNOG+17wX4PxhjjD6+uLvgc+DyV6TuA633GpwPTU1luC87VB8YYE9JCubnAGGNyPCuyxhgTQFZkjTEmgKzIGmNMAInmog7LMiIie4GtmVwtGufa3FBl+bIu1DNavqwLdMZLVbVEajPyVJE9HyKyVFVD9J5By5cdQj2j5cs6LzNac4ExxgSQFVljjAkgK7IZG+11gAxYvqwL9YyWL+s8y2htssYYE0B2JmuMMQFkRdYYYwLIiixpd9joMz+/22HjJrcDxwohmLGPiOwVkeXu0DeI2d4XkT0i8lsa80VEhrvZV4pI/WBly0TG5iJyyOf4DQxyvnIi8r2IrHU7Gf13Kst4dhz9zOf1MYwSkZ9FZIWb8ZlUlgn+73JaT/POKwMQDmwGKgH5gBVAjRTL3AOMdF93AyaFYMY+wFseHcN4oD7wWxrzrwdmAAI0BhaHYMbmwDQvjp+7/1JAffd1IWBDKv/Gnh1HP/N5fQwFKOi+jgQWA41TLBP032U7k/XpsFFVTwETgfYplmmP02EjwCdAS8mgF0gPMnpGVecB+9NZpD3wgToW4fRsUSo46Rx+ZPSUqu5U1V/c14dxnqGcsu86z46jn/k85R6XI+5opDuk/GY/6L/LVmRT77Ax5X+eM8uo8zDxQzgPCw8WfzICdHI/Rn4iIuVSme8Vf/N7rYn7UXOGiNT0KoT7EbYezpmYr5A4junkA4+PoYiEi8hyYA8wU1XTPIbB+l22Iut8xEgp5V8/f5YJJH/2/xVQQVXr4HRiOf7cVTzj9fHzxy8495/Xxemw8wsvQohIQeBTYICq/p1ydiqrBPU4ZpDP82OoqqdVNQan77+GIlIrxSJBP4ZWZNPusDHVZUQkAihMcD96ZphRVRNU9aQ7+i4QG6Rs/vDnGHtKVf9O/qipTo8ckSISHcwMIhKJU8A+VtXPUlnE0+OYUb5QOIY+WQ7i9IZ9bYpZQf9dtiKbRoeNKZaZitNhIzgdOM5Rt+U8VDKmaJtrh9NmFiqmAr3cb8cbA4dUdafXoXyJyCXJbXMi0hDndyMhiPsXnL7q1qrqsDQW8+w4+pMvBI5hCREp4r4uALQC1qVYLOi/y5708RVKVDVRRJI7bAwH3lfV1SIyGKcv9ak4/7k+FKfjxv04RS7UMt4vIu2ARDdjn2DlE5EJON8sR4vTQeYgnC8dUNWROP20XQ9sAo4BtwYrWyYydgbuFpFE4DjQLch/SJsCPYFVbpsiwBNAeZ+MXh5Hf/J5fQxLAeNFJBynwE9W1Wle/y7bbbXGGBNA1lxgjDEBZEXWGGMCyIqsMcYEkBVZY4wJICuyxhgTQFZkTa4nIkVE5B73dXMRmZbJ9fuISOnApDO5nRVZkxcUwXn60vnqA1iRNefFrpM1uZ6IJD+1bD3wD3AU2AfUApYBt6iqikgsMAwo6M7vg3MR/jjgL5wL7JsAjwBtgQLAAuDOIF90b3IQK7Im13OfGjVNVWuJSHPgS6Amzn3/83GK5mLgB6C9qu4VkZuBa1T1NhGZCzysqkvd7RVT1f3u6w9x7iz6KrjvyuQUef62WpMn/ayq2wHcW0QrAAdxzmxnurffhwNpPRfgahF5FLgAKAasxnkKmjHnsCJr8qKTPq9P4/weCLBaVZukt6KIRAFvA3Gquk1EngaiAhXU5Hz2xZfJCw7jdJmSnvVACRFpAs5j/XweOu27fnJB3ec+W7Vzdoc1uYudyZpcT1UTRGS+OJ0oHgd2p7LMKRHpDAwXkcI4vxuv4zQFjANGikjyF1/vAquAP3AeQ2lMmuyLL2OMCSBrLjDGmACyImuMMQFkRdYYYwLIiqwxxgSQFVljjAkgK7LGGBNAVmSNMSaA/h9RjEChrmUtXAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "#\n",
    "def cosTaylor(Q):\n",
    "    c = np.pi/2 + (Q - np.pi/2)**3/6 - Q\n",
    "    return c\n",
    "#\n",
    "thetaArr = np.linspace(0, np.pi, 100)\n",
    "plt.figure(figsize = (5, 4))\n",
    "plt.title(\"Cos and Taylor series for cos\")\n",
    "plt.xlabel(\"theta\")\n",
    "plt.ylabel(\"cos theta\")\n",
    "plt.plot(thetaArr, np.cos(thetaArr), color = 'b', label = 'cosine')\n",
    "plt.plot(thetaArr, cosTaylor(thetaArr), color = 'r', linestyle = '--', label = \"Taylor series\")\n",
    "plt.legend()\n",
    "plt.grid(color = 'g')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrices and linear equations\n",
    "\n",
    "This section is here for information only as you have yet to look at matrices in maths! \n",
    "\n",
    "Matrices can be defined and inverted as shown here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is the matrix A\n",
      " Matrix([[3, 7], [4, -2]])\n",
      "Its inverse is\n",
      " Matrix([[1/17, 7/34], [2/17, -3/34]])\n",
      "A multplied by its inverse is\n",
      " Matrix([[1, 0], [0, 1]])\n",
      "The inverse of A multplied y A is\n",
      " Matrix([[1, 0], [0, 1]])\n",
      "As expected, both multiplications give the unit matrix!\n"
     ]
    }
   ],
   "source": [
    "A = sp.Matrix(([3, 7], [4, -2]))\n",
    "invA = A.inv()\n",
    "AinvA = A*invA\n",
    "invAA = invA*A\n",
    "print(\"This is the matrix A\\n\",A)\n",
    "print(\"Its inverse is\\n\",invA)\n",
    "print(\"A multplied by its inverse is\\n\",AinvA)\n",
    "print(\"The inverse of A multplied y A is\\n\",invAA)\n",
    "print(\"As expected, both multiplications give the unit matrix!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Systems of linear equations can be defined and solved using matrices. For example, consider the equations:\n",
    "\n",
    "\\begin{align*}\n",
    "3x + 4y &= 7 \\\\\n",
    "x - 3y &= 2.\n",
    "\\end{align*}\n",
    "\n",
    "These can be written:\n",
    "$$\n",
    "\\left( {\\begin{array}{*{20}{c}}\n",
    "3&4\\\\\n",
    "1&{ - 3}\n",
    "\\end{array}} \\right)\\left( {\\begin{array}{*{20}{c}}\n",
    "x\\\\\n",
    "y\n",
    "\\end{array}} \\right) = \\left( {\\begin{array}{*{20}{c}}\n",
    "7\\\\\n",
    "2\n",
    "\\end{array}} \\right)\n",
    "$$\n",
    "\n",
    "Defining\n",
    "$M = \\left( {\\begin{array}{*{20}{c}}\n",
    "3&4\\\\\n",
    "1&{ - 3}\n",
    "\\end{array}} \\right)$,\n",
    "$\\vec x  = \\left( {\\begin{array}{*{20}{c}}\n",
    "x\\\\\n",
    "y\n",
    "\\end{array}} \\right)$\n",
    "and\n",
    "$\\vec c  = \\left( {\\begin{array}{*{20}{c}}\n",
    "7\\\\\n",
    "2\n",
    "\\end{array}} \\right)$,\n",
    "we have\n",
    "$M \\vec x  =  \n",
    "\\vec c$.\n",
    "The solution is obtained by multiplying through by $M^{-1}$, the inverse of $M$, and using the fact that $M^{-1}M = I$, the unit matrix:\n",
    "\n",
    "\\begin{align*}\n",
    "{M^{ - 1}}M  \\vec x  &= {M^{ - 1}} \\vec c \\\\\n",
    "  \\vec x  &= {M^{ - 1}} \n",
    "\\vec c.\n",
    "\\end{align*}\n",
    "\n",
    "Using sympy, this can be done as follows. (Note we start by redefining $x$ and $y$ so we do not use the assignments above!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equations to solve Matrix([[3*x + 4*y], [x - 3*y]]) = Matrix([[7], [2]])\n",
      "Solution is Matrix([[x], [y]]) = Matrix([[29/13], [1/13]])\n"
     ]
    }
   ],
   "source": [
    "x, y = sp.symbols('x y')\n",
    "M = sp.Matrix(([3, 4], [1, -3]))\n",
    "X = sp.Matrix((x, y))\n",
    "C = sp.Matrix((7, 2))\n",
    "invM = M.inv()\n",
    "solVect = invM*C\n",
    "print(\"Equations to solve\",M*X,'=',C)\n",
    "print(\"Solution is\",X,\"=\",solVect)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is therefore:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "x \\\\\n",
    "y\n",
    "\\end{array} \n",
    "\\right) =\n",
    "\\left( \\begin{array}{c}\n",
    "\\frac{29}{13} \\\\\n",
    "\\frac{1}{13}\n",
    "\\end{array} \\right),\n",
    "\\end{align*}\n",
    "\n",
    "or\n",
    "\n",
    "\\begin{align*}\n",
    "x &= \\frac{29}{13} \\\\\n",
    "y &= \\frac{1}{13}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This system of equations can of course easily be solved by hand, but sympy will cope with systems of equations that would cause even Bradley Cheal to despair!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Odds and ends\n",
    "\n",
    "You can print equations from sympy in LaTeX format, so they can be copied and pasted into Markdown cells or documents. An example is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sp.latex(sp.Integral(sp.sqrt(1/x), x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Student -->\n",
    "You can also get sympy to write its results to the screen in a more legible format by giving the command `sp.init_printing()`. This looks on your computer to find the most attractive printing mode available (e.g. LaTeX if you have it installed) and uses that. The disadvantage of using this mode is that you can't copy and paste the output into code cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()\n",
    "#\n",
    "sp.solve(a*x**2 + b*x + c, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to turn \"posh printing\" off, you can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-b + sqrt(-4*a*c + b**2))/(2*a), -(b + sqrt(-4*a*c + b**2))/(2*a)]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.init_printing(pretty_print = False)\n",
    "#\n",
    "sp.solve(a*x**2 + b*x + c, x)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}