Tutorial_14/Scipy-III-solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Scientific Computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we have seen how to get the required data, refine it so that it can be used for scientific analysis, and get scientific conclusions. The latter could be done by observing trends in the data and then building a suitable model, or using an existing model as predicted by theory or previous observations and trying to compare the fresh data with the existing model to present scientific conclusions. Thus, the requirement of a model is ubiquitous."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A model is mathematically presented in the form of differential equations, which could be both ordinary or partial. [Ordinary differential equations(ODEs)](https://en.wikipedia.org/wiki/Ordinary_differential_equation) involve derivatives with respect to only a single variable, whereas [Partial Differential Equations(PDEs)](https://en.wikipedia.org/wiki/Partial_differential_equation) can involve multiple variables participating in taking derivatives of a given function. Examples:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ordinary Differential Equation: $$\\frac{d^2y}{dt^2}  + a\\frac{dy}{dt} + bt = 0$$, here a and b are the parameters of the model. For a given model, they are constants."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Partial Differential Equation: $$\\frac{\\partial^2{y}}{\\partial{t}^2} + c\\frac{\\partial^2{y}}{\\partial{x}^2} = 0$$, c being the model parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, to understand the behaviour of a model, it is important to get solutions to such diffrential equations. Most of such differential equations are not feasable to be solved analytically. Thus numerical techniques must be employed to get their solutions for a set of given parameters. Here comes Python to our rescue!!!. In today's tutorial, we will try to understand how to use Python for solving ODEs. Numerical solutions to PDEs is a highly mathematically involved subject and thus beyond the scope of this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the required packages\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_ivp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First Order ODE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The package solve_ivp belonging to the module scipy.integrate will be used to numerically solve the differential equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will try to solve a first order differential equation: $$\\frac{dy}{dt} = ay + bt^2 + c$$, here a and b are the parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to define the derivative function, as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_derivative(t, y, a, b, c): \n",
    "    #t is the independent variable,y is the dependent variable, a and b are parameters. The order matters.\n",
    "    return a*y + b*(t**2) + c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An ODE of [order](https://en.wikipedia.org/wiki/Differential_equation#Equation_order) n requires n independent initial conditions so that it can be completely solved (Thus it is also known as an Initial Value Problem (IVP)). Our ODE is of order 1, thus we only need one initial condition, i.e. the value of y at x = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_span = (0, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "t_span denotes a 2-tuple of floats (t0, tf). t0 denotes the start point and it integrates till it reaches tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,5,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "t denotes the time points at which we would like the solver to compute and store the solution. If t is not specified, the solver will by default store values at those time points as specified by the algorithm for the numerical integration (NOT RECOMMENDED). Note that t must be within the interval specified by t_span."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = np.array([1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y0 is the array of initial conditions at t0. Since ours is a first order differential equation, this would be a one dimensional array. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we set the parameter values for our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1\n",
    "b = -2\n",
    "c = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our ODE is now ready to be put into the solver. Lets get going!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve_ivp(my_derivative,t_span, y0, t_eval = t, args = (a,b,c))\n",
    "#note that t_eval and args are keyword arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#uncomment below to see what object the solver returns\n",
    "#sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested in getting the solution at the time points as set by us, for that we retrieve the respective variables of the sol object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 100)\n"
     ]
    }
   ],
   "source": [
    "#uncomment below to see the solution array\n",
    "#print(sol.y)\n",
    "print(sol.y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus 'sol.y' gives a (n, n_points) shaped array, where n denotes the order of the differential equation, n_points denotes the length of the array of time points desired by us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = sol.y[0] #Why?\n",
    "time_points = sol.t\n",
    "#print(solution)\n",
    "#print(time_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with any scientific model, we would like to know how the solution looks like. So plot it!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGDCAYAAACBTdwmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deZzddX3v8fdn9n1fM5NJMslkISEkMJIAkS1RQKiIBYtYFe0Vtbi211av3nrtvVoX2trWlSq1VpBGEUVENmWHAAkkISHbZJ1JZk9m38/53j/OSTKEhCzMzPd3znk9H4/zOOuceWdOlne+v+/v+zXnnAAAAOBfku8AAAAAiKCYAQAABATFDAAAICAoZgAAAAFBMQMAAAgIihkAAEBAUMwAJCQzm2lmzsxSzvDr32dmD090LgCJjWIGIOaZ2Qoze9bMus3soJk9Y2ZvmcD3f12Jc87d6Zx7+0R9DwCQpDP6nyIABIWZ5Um6X9LHJa2WlCbprZKGfeYCgDPBiBmAWDdXkpxzP3fOhZxzg865h51zG80sycy+ZGZ7zazNzH5qZvnHexMz22Nmq8bd/z9m9rPo3Sej111m1mdmF5jZzWb29LjXX2hmL0ZH7V40swvHPfe4mf3f6Eher5k9bGYlE/+jABDrKGYAYt12SSEz+08zu8rMCsc9d3P0cpmkWkk5kr5zBt/j4uh1gXMuxzn33PgnzaxI0u8k/aukYkn/JOl3ZlY87mU3SfqQpDJFRvX+5xnkABDnKGYAYppzrkfSCklO0r9Lajez+8ysXNL7JP2Tc26Xc65P0hck3XimE/7fwNWSdjjn/ss5N+ac+7mkrZL+ZNxr/sM5t905N6jIIdclE5wBQBygmAGIec65Lc65m51z1ZIWSZom6dvR673jXrpXkbm15RMc4djvc/h7VY273zLu9oAio3cA8BoUMwBxxTm3VdJPFCloByTNGPd0jaQxSa3H+dJ+SVnj7leMf9uTfNtjv8/h77X/5IkB4CiKGYCYZmbzzeyvzaw6en+6pPdKWiPp55I+a2azzCxH0tck/bdzbuw4b7VekcOcqWZWL+n6cc+1SworMk/teB6QNNfMbjKzFDP7M0lnKXK2KACcMpbLABDreiUtk/RXZlYgqUuRQvQ5SX2KHGZ8UlKGpIckffIE7/O/FSlyhyQ9IekuSUWS5JwbMLOvSnrGzFIlXTn+C51znWZ2jaR/kfR9SQ2SrnHOdUzgrxNAAjDnTjZCDwAAgKnAoUwAAICAoJgBAAAEBMUMAAAgIChmAAAAAUExAwAACIi4WC6jpKTEzZw503cMAACAk1q3bl2Hc670eM/FRTGbOXOm1q5d6zsGAADASZnZsVu4HcGhTAAAgICgmAEAAAQExQwAACAgKGYAAAABQTEDAAAICIoZAABAQFDMAAAAAoJiBgAAEBAUMwAAgICgmAEAAAQExQwAACAg4mKvTACYSs45DY+FNTQa0vBYWMOjYQ2PRW6PhMIKh53Gwk6h6HXYOUmSRb/ezJRkUkpSktJSTKnJSdHbScpMS1ZmarKy0pKVnpIkMztxEABxh2IGICENjIyps29EhwZGdGhgVF0DI+oaGNWhgRH1DI6pd2hUPUOj6h0aU+/QmPqHx9Q/MqaB4ZAGRkMKhd2kZzSTMlOTlZeRqrzMlOh1qvIyUlSUna6S3DSVZKerOCdNxTnpKs9LV1luhpKTKHNArKKYAYgbI2NhtfUOqbVnWO29Q2rvHVZb77DaeobV0Tesjv4RdfYNq7NvRIOjoRO+T056inIzIkUoNyNFJTlpmlmSrazUZGWlJys7LeXIyFZ6apLSUyKjW+kpSZHRr2RTcpIpJSlJyUmR0bHDNS46eKawcxoNhTUWilyPhpxGQmENjYQ0MDKmgdGQhkZC6h8JRUri4Jh6hkbV1juk7a2jOtQ/ov6R1/8akpNMZbnpqszPUGV+pqoLM1VTnKUZRdmaUZylyvwMpSQziwUIKooZgJgwFgqrpWdIzd1DOtA1qANdkevm7kG19AyppTtSvo6VZFJJTnrkkpuu2SXZR0aYirLTVJiVpsKsVBVEr/MzU2OmuAyOhNTZHyma7b3Dau0dUkv3kA50DamlZ1Bbmnv0yJZWjYyFj3xNarKppihLdWW5qivP0ZyyHNWV5Wp2WbbSU5I9/moASBQzAAERDju19Axp38EBNR4cUOOhQTUdGtD+Q4NqOhQpX8cePszPTFVlfoYq8jO0aFq+KvIzVJGXofL8DJXlpqs0N13F2elxe2gvMy1Z1WlZqi7MOuFrDv9c93YOaN/Bfu3tHNCu9n5tb+vVI1taj/xMU5JMdeW5WjQtTwun5WlRVb7OmpanrDT+mQCmkjk3+fMkJlt9fb1bu3at7xgATmJ4LKTGg4Pa2xkpCHs7+7WnM1LEmg4NaiR0dGTHTKrMy1BVYaaqC7NUXZipaQWRS1VB5DBddjql4c0YHgtpT8eAtrf2amtLjzbt79Gm/d3q7B+RFBltXFCZp/NmFOq8GYU6t6ZQ1YWZnJAAvElmts45V3/c5yhmACZSKOzUdGhAuzr6tbu9X3s6+7W7I3I50DWo8YNeuekpqinO0szibFUXZaqmKEs1RVmaXpilaQWZSkuJjUOK8cQ5p9aeYW3a360NTV16ad8hrd/XdWQ+W1luui6YXayLZpfogtnFml504tE6AMdHMQMw4XqGRrWrvV872/q0s71Pu9r7taujT3s6B14zpyk3PUWzSrM1szhbM0uyNaskSzOKI/cLs1IZfYkBY6GwtrX26qV9XXph90E9t7PzyHy+mqIsXTi7WJfOK9OKuhLlMIoJnBTFDMAZcc6pvW9YDW19r7u09R6daJ+SZKopzlJtSY5ml2artjRbs0pyVFuareLsNMpXnHHOaXtrn57d2aFnGjr1/K5O9Q6PKS05Sctqi3T5/DJdPr9MM4qzfUcFAoliBuANOefU1jus7a292tHapx1th6/71D04euR1Oekpml2WozmlOZpdlh29zlFNUZZSY+RMRky80VBYL+45qMe2tukPW9u0q71fkjS/IldXn12pdyyu1OzSHM8pgeCgmAE4orNvWNtae7W9pVfb2/oi16296hkaO/KagqxUzS3L1ZzyHNWVHV1SoTwvndEvnNTezn498mqrfr+pRev2HpJ0tKS9c8k0RtKQ8ChmQALqGx7T9mgB2xotX9tbe9XRN3LkNfmZqZpXHlnPam70uq4sVyU5HH7ExGjuHtTvX2nRA680a220pL1lZqGuP69a7zi7UrkZqZ4TAlOPYgbEsbFQWLs7+rWlpVfbWnq0LVrEmg4NHnlNVlqy6spzNb88V3MrcjW3PEfzynNVmssIGKbOga5B3fvyft3zUpN2tfcrIzVJVyys0J/VT9cFs4v5vYiEQTED4sDhifhbmyNrTkWue9XQ1ndk/a/kJFNtSbbmVeRqfkWu5lXkaX5FrqoKMpUUp4usIvY457S+sUu/XNek3244oJ6hMc0uzdb7l8/Qu8+rVh6jaIhzFDMgxgyNhtTQ1qctzT3a2nK0iB1e+FOSyvPSNT9avOZX5mpeeR7b6iDmDI2GdP/GZv3Xmr3a0NilzNRkvWtplT544QzNr8jzHQ+YFBQzIKCci2yXs7W5V68eLmHNPdrV0X9kq5yM1CTNK8+NlLDK3CNlrDA7zXN6YGJtbOrSz9bs1W/WH9DwWFiXzivVLRfX6oJaDnMivlDMgAAYGg1pR2tkFGxLS8+R0bCugaPLUVQVZGpBZa4WVOYdKWIzi7Pjdq9H4Hi6Bkb0szV79ZNn96ijb0SLq/P10Ytn68pFFfxZQFygmAFT6PCaYK82R8tXc6+2HDMKlpmarLkVuTprXAmbV5Gr/Ezm1gCHDY2GdM9LTfrRU7u1u6Nfs0qy9amVc/TOc6ooaIhpFDNgkgyPHR0F29oSKWBbmnt06ASjYIcvNUVZ/MMCnKJQ2OnhzS361z82aEtzj2pLs/WZVXN19dmV/DlCTKKYAW+Sc07t0VGw8QVsZ/vx54IdORxZmccoGDBBwmGnh19t0T8/skPbWntVV5ajz75trq5aVMEcNMQUihlwGk7ljMhp+RmaX5n3mvlgs0qYCwZMhXDY6YFNzfr2ozvU0NanpTUF+tLVZ+m8GYW+owGnhGIGHIdzTs3dQ9ra0qMt0TXB3uiMyAWVuZEyVpGn/CxGwQDfQmGne9Y16baHt6mtd1hXn12pv71yvmqKs3xHA94QxQwJr394TNtae48uzhotYeP3h4zMBYssRRE5DMkZkUAs6B8e078/tUs/fGKXQmGnD144Q59aWcd2TwgsihkSxlgorD2d/dra0qttLb3a0tyrba09ajx4dHuinPSUIyvjRxZnjZwRyWrjQGxr7RnSPz68Tb9Y16TSnHT972vO0jWLK5l/hsChmCHuHF6YdVu0gB3eH7KhvU8jY0e3J5p1eHui8lzNi46EVRdm8hc1EMc2NHbpS7/epFf2d2vFnBL9/bULVVua4zsWcATFDDGts29Y21v7tL21V9tae7W9JXLdO+4wZEVehuZVRMrXvGgJm1OWo4xUticCElEo7HTn83v1rYe2aXg0rI9eUqtbL5vD3wkIBIoZYkL34Kh2RMvXjmgR297aq46+o2dD5mWkaH5FnuZW5EQLWJ7mlueoIIvtiQC8XlvvkP7hga269+X9mlOWo29dv1hLazh7E35RzBAo3QOj2tHWqx1tfdrR2qcdbZEC1tozfOQ1WWnJqivP1dyyHM2ryNXc8silPC+dw5AATtsT29v1+Xs2qrVnSB+5uFafXTWX0TN4QzHDlHPOqbN/RA1tfdrR1qeG1sj8rx2tfWrrPVrAMlOTNbss+0jxmlueo7qyXFUVZCqJsyEBTKCeoVF97XdbdPeLjZpdmq3bbjiH0TN4QTHDpAmFnQ50DaqhrU8NbX3a2R65bmjve83m3NlpyZpTnqs5pTmR8kUBA+DJk9HRs5aeId162Rx9emWdUpKTfMdCAqGY4U3rGx7T7vZ+7ero0872fu1s79POtj7t7ujXcPQsSEkqzk7T7LIczSnL0ezSHNWVRUpYRV4GhyABBEbv0Ki+8ttX9ct1TTpvRqG+/WdLNL2IhWkxNShmOCWjobCaDg1qd0efdrX3a3dHv3ZFy9j4+V9JJtUUZWl2aY5ml+Vodmm2aktzNKc0R4XZTMIHEDt+s36/vnTvJsmkf3j32bpm8TTfkZAA3qiYpUx1GPgVDjsd6B7Uno4B7e7s156OSAHb09GvfQcHNBY+WtTzMlJUW5qjFXNKVVuafaSAzSjOUnoKk2YBxL5rl1Tp3JpCferul/WJu17WU9s79OV3nqWsNP55hB/8zotDo6GwDnQNam/ngPYeHNDejn7t6ezXns4B7Ts4cGQBVimyF+TM4sgirFcsqlBtSbZqS7M1qyRHhVmpHH4EEPemF2Vp9Ucv0Lcf3a7vPb5T6xu79MP3n6eZJdm+oyEBcSgzRvUMjarx4ID2RcvW4cvezgHt7xo8sgm3dLR8zSjOil5na2ZJlmaVZKs8N4PJ9wAQ9eT2dn3q7pcVCjv9y41LdPn8ct+REIeYYxaDBkdC2t81oMZDg2o6NKimaPFqPDSgxoOD6h4cfc3rC7NSVVOUpZribM0oylJNcZZmFGVpRnG2ynLTKV8AcIoaDw7oYz9bp80HevSZVXX61OV1/B2KCcUcs4BxzqlncEz7uwYjl0MD425Hrsevdi9JaSlJqi7M1PTCLC2ZXqDphVmaUZyl6UWRCxtwA8DEmF6UpXs+fqG+eO8mffvRHdrY1K1/fs8S5Wfx9ywmH8VsEvQPj6m5e0jN3YOR667I7QPdQzrQNajmrkH1j4Re8zXpKUmqKshUVWGmzpqWp+rCLFUXZqq6MFNVBVmMegHAFMpITdZtNyzWkpoC/f1vN+u67z2jO25+C/POMOk4lHkaxkJhdfSNqLVnKHLpHVZr95BaovdborfHb659WElOmqYVZGpafqYqCzJUVZCpaQWZR8pYcXYaE+0BIIBe2H1QH/2vyL8xt3+gXm+ZWeQ5EWIdc8zepFcP9OiD//GCOvqGdeyPKznJVJqTrvL8DFXkpasiL0OVBZmqzM+I3M7PVHl+OstLAEAM29PRrw//5EU1HRrUN69frHctrfIdCTEsJueYmdmVkv5FUrKkHznnvu4rS0lumlbOL1NZXobK89JVnpuh8rwMleWlqyQnXckcYgSAuDazJFu/+ssL9bGfrdNn/nu99nT269Mr6zjSgQkXyBEzM0uWtF3S2yQ1SXpR0nudc68e7/XxeFYmACB4RsbC+sKvXtE9LzXp3Uur9I3rFyuVfTZxmmJxxOx8SQ3OuV2SZGZ3S7pW0nGLGQAAUyEtJUm33bBYM4uz9I+PbFf34Ki++75zlZHKdBVMjKDW/CpJjePuN0UfO8LMbjGztWa2tr29fUrDAQASl5npkyvr9NXrFumP29r0gR+/8Lq1JYEzFdRidlLOududc/XOufrS0lLfcQAACeZ9y2bo3967VC83HtKNt69Re++w70iIA0EtZvslTR93vzr6GAAAgXHN4mn60Qffoj0d/brhB8+q8eCA70iIcUEtZi9KqjOzWWaWJulGSfd5zgQAwOtcMrdUd35kmQ4NjOqGHzynPR39viMhhgWymDnnxiR9QtJDkrZIWu2c2+w3FQAAx3duTaH++6PLNRIK68bb11DOcMYCWcwkyTn3gHNurnNutnPuq77zAADwRuZX5Omujyw7Us52U85wBgJbzAAAiDXjy9l7KWc4AxQzAAAmEOUMbwbFDACACXZsOeNsTZwqihkAAJNgfkWe7vwfyzQwMqb3//h51jnDKaGYAQAwSRZU5uk/PnS+WnuG9YE72CEAJ0cxAwBgEp03o1A/fP95amjr1V/85EUNjoR8R0KAUcwAAJhkF88t1bf/bKnW7Tukj9+5TiNjYd+REFAUMwAApsDViyv1tevO1uPb2vXXv9igcNj5joQASvEdAACARPHe82vUNTCqbzy4VdPyM/SFdyzwHQkBQzEDAGAKfeySWh3oGtQPn9ylmuIsvW/ZDN+RECAUMwAAppCZ6ct/cpaaDg3o736zWVUFmbp0XpnvWAgI5pgBADDFUpKT9G83nau55bn6xF0va0tzj+9ICAiKGQAAHuSkp+iOm+uVnZ6sD//kRbX2DPmOhACgmAEA4EllfqbuuPkt6h4c1Yd/8qIGRsZ8R4JnFDMAADxaOC1f373pXG1p7tHnfrlRzrGMRiKjmAEA4Nll88v0uSvm63cbm/WDJ3b5jgOPKGYAAATAxy6p1dWLK/XNh7bqie3tvuPAE4oZAAABYGb61vWLNa88V5+86yXt7ez3HQkeUMwAAAiIrLQU3f7+epmZbvnpOvUPczJAoqGYAQAQIDXFWfrOTUu1o61Xn/vlBk4GSDAUMwAAAuatdaX62yvn64FXWvSjp3b7joMpRDEDACCAbrm4VlcurNA3Htyql/cd8h0HU4RiBgBAAJmZvnH9YlXkZ+gTd72s7oFR35EwBShmAAAEVH5mqv7tvUvV2jPEfLMEQTEDACDAltYU6vNXzdfDr7bqP5/d4zsOJhnFDACAgPuLFbO0akGZvvbAVr3S1O07DiYRxQwAgIAzM912wzkqyUnTrXe9pJ4h5pvFK4oZAAAxoCArTf9201Lt7xrU3/16k+84mCQUMwAAYsR5M4r0qcvr9Ov1B3T/xgO+42ASUMwAAIght142W+dML9AX792klu4h33EwwShmAADEkJTkJP3ze87RyFiYJTTiEMUMAIAYU1uaoy9evUBP7ejQf63Z6zsOJhDFDACAGPS+ZTW6dF6pvvbAFjW09fmOgwlCMQMAIAaZmb75p4uVkZqsv1q9XqOhsO9ImAAUMwAAYlRZXob+4bqztbGpW9/5Y4PvOJgAFDMAAGLYVWdX6tol0/S9xxu0taXHdxy8SRQzAABi3Jf/ZKFyM1L1t/e8olCYszRjGcUMAIAYV5Sdpi//yVna0Nil/3hmt+84eBMoZgAAxIF3njNNK+eX6baHt2lf54DvODhDFDMAAOKAmen/XbdIKUlJ+vyvNrLwbIyimAEAECcq8zP1+avm69mdnVq9ttF3HJwBihkAAHHkpvNrdP6sIv2/321Raw97acYaihkAAHEkKcn09XefrZGxsL78m82+4+A0UcwAAIgztaU5+tTKOj24uUWPb2vzHQengWIGAEAc+h9vnaXakmx9+b7NGhoN+Y6DU0QxAwAgDqWnJOsr1y7U3s4B/fCJXb7j4BRRzAAAiFNvrSvV1WdX6nuPN7C2WYygmAEAEMe+dM0CJSeZvvJbTgSIBRQzAADiWGV+pj6zqk5/2NqmR15t9R0HJ0ExAwAgzn3oolmqK8vR/7lvswZHOBEgyChmAADEudTkJP39tYu0v2tQ332swXccvAGKGQAACeCC2cV615Jpuv2pXWo8yIkAQUUxAwAgQfztVfOVZNLXH9zqOwpOgGIGAECCqMzP1C0Xz9bvNjZr3d6DvuPgOChmAAAkkI9dUquy3HT9/f1bFA4733FwDIoZAAAJJCstRZ+7Yp42NHbptxsP+I6DY1DMAABIMH96brUWVeXpG7/fyvIZAeOlmJnZt8xsq5ltNLN7zaxg3HNfMLMGM9tmZlf4yAcAQDxLSjJ96eqzdKB7SD96in00g8TXiNkjkhY55xZL2i7pC5JkZmdJulHSQklXSvqemSV7yggAQNxaXlusKxaW6/tP7FRbz5DvOIjyUsyccw8758aid9dIqo7evlbS3c65YefcbkkNks73kREAgHj3hasWaDQU1m0Pb/MdBVFBmGP2YUm/j96uktQ47rmm6GMAAGCCzSzJ1gcvmKlfrGvStpZe33GgSSxmZvaomW06zuXaca/5oqQxSXeewfvfYmZrzWxte3v7REYHACBh3HrZHOWkpTBqFhApk/XGzrlVb/S8md0s6RpJK51zhxdS2S9p+riXVUcfO9773y7pdkmqr69nIRYAAM5AYXaabrm4Vv/4yHa9tO+Qzq0p9B0pofk6K/NKSX8j6Z3OufEbdt0n6UYzSzezWZLqJL3gIyMAAIniwytmqSQnTd98cKuOjpXAB19zzL4jKVfSI2a23sx+IEnOuc2SVkt6VdKDkm51zrHACgAAkyg7PUWfuGyO1uw6qKd2dPiOk9AsHppxfX29W7t2re8YAADErOGxkC6/7QkVZafpN7depKQk8x0pbpnZOudc/fGeC8JZmQAAwLP0lGR99m1z9cr+bv1+U4vvOAmLYgYAACRJ1y2tUl1Zjv7xkW0aC4V9x0lIFDMAACBJSk4y/c8r5mlXe7/ueanJd5yERDEDAABHvP2sci2ZXqBvP7pDQ6OcfzfVKGYAAOAIM9Pnrpin5u4hrV7bePIvwISimAEAgNe4cHax6mcU6vuP79TwGKNmU4liBgAAXsPM9OlVdWruHtIv1jLXbCpRzAAAwOusmFOic2sK9P3Hd2pkjDM0pwrFDAAAvE5k1Gyu9ncN6pfrGDWbKhQzAABwXBfXleic6QX67mMNjJpNEYoZAAA4LjPTZ1bWaX/XoO59mVGzqUAxAwAAJ3TpvFItrs7Xdx5r0Ci7AUw6ihkAADghM9OnV9ap8eCg7n15v+84cY9iBgAA3tDl88u0qCpP332sgT00JxnFDAAAvCEz06cur9PezgHdv7HZd5y4RjEDAAAntWpBuerKcvSDJ3bKOec7TtyimAEAgJNKSjJ99JLZ2trSq8e3t/uOE7coZgAA4JS885xpqszP0A+f2Ok7StyimAEAgFOSlpKkv1gxS2t2HdT6xi7fceISxQwAAJyyG8+vUV5Gin7wOKNmk4FiBgAATllOeoo+cMFMPfRqi3a29/mOE3coZgAA4LTcfNFMpSUn6d+f3OU7StyhmAEAgNNSkpOuG+qr9auX9qutZ8h3nLhCMQMAAKftI2+t1Vg4rB8/s9t3lLhCMQMAAKdtRnG2rjq7Unet2aeeoVHfceIGxQwAAJyRj18yW73DY7r7hX2+o8QNihkAADgji6rytWxWkf7z2b1sbj5BKGYAAOCMfeiiWdrfNahHt7T6jhIXKGYAAOCMve2sclUXZuqOp/f4jhIXKGYAAOCMJSeZbr5wpl7Yc1Cb9nf7jhPzKGYAAOBNec9bpis7LVl3sHTGm0YxAwAAb0peRqquP69a929oVnvvsO84MY1iBgAA3rSbL5qlkVBYdz6/13eUmHbSYmZmnzSzwqkIAwAAYtOskmxdPr9MP1uzV8NjId9xYtapjJiVS3rRzFab2ZVmZpMdCgAAxJ4PXzRLHX0jun9Ds+8oMeukxcw59yVJdZJ+LOlmSTvM7GtmNnuSswEAgBhy0ZxizS3P0R3P7JZzznecmHRKc8xc5KfbEr2MSSqU9Esz++YkZgMAADHEzPShi2Zp84EevbD7oO84MelU5ph92szWSfqmpGckne2c+7ik8yT96STnAwAAMeS6pVUqyErVT9dwEsCZSDmF1xRJerdz7jU/Yedc2MyumZxYAAAgFmWkJuv6c6v1k2f3qL13WKW56b4jxZRTmWP25WNL2bjntkx8JAAAEMtuWlajsbDT6rWNvqPEHNYxAwAAE6q2NEcXzSnWXc/vUyjMSQCng2IGAAAm3J8vm6H9XYN6fFub7ygxhWIGAAAm3KqzylWWm66fcRLAaaGYAQCACZeanKQb3zJdj29vV+PBAd9xYgbFDAAATIobz6+RSfr5C/t8R4kZFDMAADApphVkauWCcq1e26iRsbDvODGBYgYAACbN+5bVqKNvRA9ubvEdJSZQzAAAwKS5uK5UNUVZupOTAE4JxQwAAEyapCTTTctq9Pzug9rR2us7TuBRzAAAwKS64bxqpSUn6c7nOQngZChmAABgUhXnpOvtC8v16/X7NTwW8h0n0ChmAABg0r2nfrq6Bkb1yKutvqMEGsUMAABMuovmlKiqIFOr1zb5jhJoFDMAADDpkpNMf3petZ7a0a79XYO+4wQWxQwAAEyJG86rlnPSPesYNTsRihkAAJgS04uydNGcYq1e26hw2PmOE0gUMwAAMGXeUz9dTYcGtWZXp+8ogUQxAwAAU+aKhRXKy0jR6rWNvtuelOEAABBzSURBVKMEktdiZmZ/bWbOzEqi983M/tXMGsxso5md6zMfAACYWBmpybp2SZV+v6lF3YOjvuMEjrdiZmbTJb1d0vhlgK+SVBe93CLp+x6iAQCASfSe+ukaHgvrvg0HfEcJHJ8jZv8s6W8kjZ/9d62kn7qINZIKzKzSSzoAADApFlXlaUFlnla/yOHMY3kpZmZ2raT9zrkNxzxVJWn8p9QUfQwAAMQJM9N76qv1yv5uvXqgx3ecQJm0YmZmj5rZpuNcrpX0vyT93Zt8/1vMbK2ZrW1vb5+Y0AAAYEq8a0mV0pKTOAngGJNWzJxzq5xzi469SNolaZakDWa2R1K1pJfMrELSfknTx71NdfSx473/7c65eudcfWlp6WT9MgAAwCQozE7TqrPK9NsNBzQaCvuOExhTfijTOfeKc67MOTfTOTdTkcOV5zrnWiTdJ+kD0bMzl0vqds41T3VGAAAw+a5bWq3O/hE9tYMjX4cFbR2zBxQZUWuQ9O+S/tJvHAAAMFkumVuqwqxU3fsyZ2celuI7QHTU7PBtJ+lWf2kAAMBUSUtJ0jWLp2n12kb1Do0qNyPVdyTvgjZiBgAAEsi7llZpeCysBze1+I4SCBQzAADgzbk1BZpRnKV7Xz7uuX4Jh2IGAAC8MTO9a0mVntvVqebuQd9xvKOYAQAAr65bWiXnpN+s5yQAihkAAPBqZkm2ltYU6NcczqSYAQAA/969tEpbW3oTfosmihkAAPDu6sXTlJJk+vX6xB41o5gBAADvirLTdOm8Uv1m/X6Fws53HG8oZgAAIBCuW1qt1p5hPbez03cUbyhmAAAgEFYuKFNueop+9XKT7yjeUMwAAEAgZKQm66qzK/TQphYNjYZ8x/GCYgYAAALjmsXT1D8S0uPb2nxH8YJiBgAAAuPC2cUqyk7Tbzc2+47iBcUMAAAERkpykq5cVKE/bmnTwMiY7zhTjmIGAAAC5ZrFlRocDemPWxPvcCbFDAAABMqyWcUqzU3X/RsS73AmxQwAAARKcpLpHYsq9Ni2NvUNJ9bhTIoZAAAInGvOmabhsbAefbXVd5QpRTEDAACBc15NoSrzM3T/xgO+o0wpihkAAAicpCTT1WdX6ont7eoeHPUdZ8pQzAAAQCBdc840jYacHt7c4jvKlKGYAQCAQDqnOl/TizJ1fwItNksxAwAAgWRmuvrsaXqmoUOH+kd8x5kSFDMAABBY1yyu1FjY6cEEOZxJMQMAAIG1cFqeZpVkJ8zZmRQzAAAQWJHDmZV6bmenOvuGfceZdBQzAAAQaFcuqlDYSY9uif/FZilmAAAg0BZOy1N1YaYe2kwxAwAA8MrMdOXCCj29o0O9Q/G92CzFDAAABN6Viyo0EgrrsW3tvqNMKooZAAAIvHNrClWam66HNsX3shkUMwAAEHhJSaa3nVWux7a1aWg05DvOpKGYAQCAmHDlwgoNjIT01I4O31EmDcUMAADEhOW1xcrLSNGDcXw4k2IGAABiQlpKklYtKNejW1o1Ggr7jjMpKGYAACBmXLGoQt2Do3ph90HfUSYFxQwAAMSMi+tKlZmaHLeHMylmAAAgZmSmJevSeaV6aHOLwmHnO86Eo5gBAICYcsXCCrX1Duvlxi7fUSYcxQwAAMSUy+aXKTXZ9NDm+DucSTEDAAAxJT8zVRfOLtGDm1rkXHwdzqSYAQCAmHPFwgrtOzigba29vqNMKIoZAACIOasWlEmS/rClzXOSiUUxAwAAMacsL0OLq/P16JZW31EmFMUMAADEpJXzy7W+sUsdfcO+o0wYihkAAIhJKxeUyTnpj1vj53AmxQwAAMSkhdPyVJmfoT/E0eFMihkAAIhJZqbL55fpqR0dGhoN+Y4zIShmAAAgZq1aUK6BkZDW7Or0HWVCUMwAAEDMumB2sTJTk+Nm2QyKGQAAiFkZqclaUVeiP2xpjYtdAChmAAAgpq1aUKYD3UPa0hz7uwBQzAAAQEy7bP7hXQBi/+xMihkAAIhpZbkZOmd6gR6Ng/XMKGYAACDmrZpfpg2NXWrrHfId5U2hmAEAgJi3ckG5JOmxGB81o5gBAICYt6AyV1UFmXo0xpfNoJgBAICYZ2ZauaBMT8f4LgDeipmZfdLMtprZZjP75rjHv2BmDWa2zcyu8JUPAADElpULyjU4GtJzMbwLQIqPb2pml0m6VtI5zrlhMyuLPn6WpBslLZQ0TdKjZjbXORe71RcAAEyJZbOKlJGapCe2teuyeWW+45wRXyNmH5f0defcsCQ55w4fEL5W0t3OuWHn3G5JDZLO95QRAADEkIzUZC2vLdaT29t9RzljvorZXElvNbPnzewJM3tL9PEqSY3jXtcUfex1zOwWM1trZmvb22P3AwAAABPn0rml2tXRr32dA76jnJFJK2Zm9qiZbTrO5VpFDqEWSVou6XOSVpuZnc77O+dud87VO+fqS0tLJ+FXAAAAYs0l0UOYj2+PzbMzJ22OmXNu1YmeM7OPS/qVi+w2+oKZhSWVSNovafq4l1ZHHwMAADipWSXZmlGcpSe2tesDF8z0Hee0+TqU+WtJl0mSmc2VlCapQ9J9km40s3QzmyWpTtILnjICAIAYdOncUj27szMml83wVczukFRrZpsk3S3pgy5is6TVkl6V9KCkWzkjEwAAnI5L5pVqcDSktXsO+Y5y2rwsl+GcG5H05yd47quSvjq1iQAAQLxYXlustJQkPb6tTSvqSnzHOS2s/A8AAOJKVlqKls0q0uMxuGwGxQwAAMSdS+aWqqGtT02HYmvZDIoZAACIO5fOiyyl9USMjZpRzAAAQNyZXZqjqoJMPbGNYgYAAOCVmemSeaV6pqFDI2Nh33FOGcUMAADEpUvnlqp/JKR1e2Nn2QyKGQAAiEsXzilRarLF1PZMFDMAABCXctJTVD+jKKbmmVHMAABA3Lp0Xqm2tvSqpXvId5RTQjEDAABx65LoshlPxsiyGRQzAAAQt+aV56o0N11PN3T4jnJKKGYAACBumZlWzCnRMw0dCoed7zgnRTEDAABx7aI5JersH9HWll7fUU6KYgYAAOLaijklkqRnYuBwJsUMAADEtYr8DM0py4mJeWYUMwAAEPdWzCnR87s7NTwW8h3lDVHMAABA3Fsxp0RDo2G9tLfLd5Q3RDEDAABxb1ltkZKTLPDzzChmAAAg7uVmpGrJ9AI9RTEDAADwb8WcEr3S1KXugVHfUU6IYgYAABLCiroShZ303K5O31FOiGIGAAASwpLpBcpOS9bTDcHdN5NiBgAAEkJqcpKW1xbrmQZGzAAAALy7aE6Jdnf0q+nQgO8ox0UxAwAACWNFXbC3Z6KYAQCAhFFXlqOy3HQ9HdDDmRQzAACQMMxMK+aU6JmGDoXDznec16GYAQCAhLKirkQH+0e0paXHd5TXoZgBAICEctGc4M4zo5gBAICEUp6Xodml2Vqz66DvKK9DMQMAAAlneW2xXth9UGOhsO8or0ExAwAACWd5bbH6hse0+UCw5plRzAAAQMJZXlssKXj7ZlLMAABAwinNTdecshytoZgBAAD4d0FtsV7cfVCjAZpnRjEDAAAJaXltsfpHQtq0v9t3lCMoZgAAICEtqy2SFKx5ZhQzAACQkEpy0jW3PCdQ65lRzAAAQMK6oLZYa/cEZ54ZxQwAACSs5bXFGhgJaWNTMOaZUcwAAEDCWhZdzywoy2ZQzAAAQMIqyk7T/IpcihkAAEAQLK8t1to9hzQy5n+eGcUMAAAktOW1xRocDWljU5fvKBQzAACQ2JbNKpJZMOaZUcwAAEBCK8xO0/yKvEAsNEsxAwAACW95bZHW7T2k4bGQ1xwUMwAAkPCW1xZraDSsDY1+1zOjmAEAgIQXlHlmFDMAAJDwCrLS9I5FlcrLSPGaw+93BwAACIjvvu9c3xEYMQMAAAgKihkAAEBAUMwAAAACgmIGAAAQEBQzAACAgKCYAQAABATFDAAAICC8FDMzW2Jma8xsvZmtNbPzo4+bmf2rmTWY2UYz87+gCAAAwBTxNWL2TUlfcc4tkfR30fuSdJWkuujlFknf9xMPAABg6vkqZk5SXvR2vqQD0dvXSvqpi1gjqcDMKn0EBAAAmGq+tmT6jKSHzOw2RcrhhdHHqyQ1jntdU/Sx5mPfwMxuUWRUTTU1NZMaFgAAYCpMWjEzs0clVRznqS9KWinps865e8zsPZJ+LGnV6by/c+52SbdLUn19vXuTcQEAALybtGLmnDth0TKzn0r6dPTuLyT9KHp7v6Tp415aHX0MAAAg7vk6lHlA0iWSHpd0uaQd0cfvk/QJM7tb0jJJ3c651x3GPNa6des6zGzvJGUdr0RSxxR8H5w6PpPg4TMJJj6X4OEzCaap+FxmnOgJX8XsI5L+xcxSJA0pOldM0gOS3iGpQdKApA+dyps550onI+SxzGytc65+Kr4XTg2fSfDwmQQTn0vw8JkEk+/PxUsxc849Lem84zzuJN069YkAAAD8Y+V/AACAgKCYnZ7bfQfA6/CZBA+fSTDxuQQPn0kwef1cLHL0EAAAAL4xYgYAABAQFLNTYGZXmtm26Obqn/edB5KZ3WFmbWa2yXcWRJjZdDN7zMxeNbPNZvbpk38VJpuZZZjZC2a2Ifq5fMV3JkSYWbKZvWxm9/vOAsnM9pjZK2a23szWesvBocw3ZmbJkrZLepsiW0S9KOm9zrlXvQZLcGZ2saQ+RfZWXeQ7D6TovraVzrmXzCxX0jpJ7+LPil9mZpKynXN9ZpYq6WlJn47uRwyPzOyvJNVLynPOXeM7T6Izsz2S6p1zXteWY8Ts5M6X1OCc2+WcG5F0tyKbrcMj59yTkg76zoGjnHPNzrmXord7JW1RZK9beOQi+qJ3U6MX/kfumZlVS7paR3e+ASRRzE7FiTZWB3ACZjZT0lJJz/tNAunIIbP1ktokPeKc43Px79uS/kZS2HcQHOEkPWxm68zslpO+epJQzABMKDPLkXSPpM8453p854HknAs555Yosv/w+WbG4X+PzOwaSW3OuXW+s+A1VjjnzpV0laRbo1NmphzF7OTYWB04RdE5TPdIutM59yvfefBazrkuSY9JutJ3lgR3kaR3Ruc03S3pcjP7md9IcM7tj163SbpXkalMU45idnIvSqozs1lmlibpRkU2WwcwTnSS+Y8lbXHO/ZPvPIgws1IzK4jezlTkRKatflMlNufcF5xz1c65mYr8m/JH59yfe46V0MwsO3rSkswsW9LbJXk5659idhLOuTFJn5D0kCKTmVc75zb7TQUz+7mk5yTNM7MmM/sL35mgiyS9X5H//a+PXt7hOxRUKekxM9uoyH80H3HOsTwD8Frlkp42sw2SXpD0O+fcgz6CsFwGAABAQDBiBgAAEBAUMwAAgICgmAEAAAQExQwAACAgKGYAAAABQTEDgGOYWYGZ/aXvHAASD8UMAF6vQBLFDMCUo5gBwOt9XdLs6CK53/IdBkDiYIFZADiGmc2UdL9zjs2+AUwpRswAAAACgmIGAAAQEBQzAHi9Xkm5vkMASDwUMwA4hnOuU9IzZraJyf8AphKT/wEAAAKCETMAAICAoJgBAAAEBMUMAAAgIChmAAAAAUExAwAACAiKGQAAQEBQzAAAAAKCYgYAABAQ/x+PW5UKj1RH8QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10,6))\n",
    "plt.plot(time_points, solution)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n",
    "plt.title('Solution')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation we discussed can be solved analytically. It is encouraged you try it out !!!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher Order ODEs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solver cannot solve a higher order ODE as it is, since it is more computationally expensive and error prone, given the algorithms implemented by the solver. Instead, we are supposed to manually break up the higher order ODE into several first order ODEs and then solve them together. In principle, an $n^{th}$ order ODE can be broken up into n first order ODEs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will try to solve the ODE representing a [Damped Simple Harmonic Oscillator(DSHO)](http://spiff.rit.edu/classes/phys312/workshops/w5b/damped_theory.html):\n",
    "$$m\\frac{d^2x}{dt^2} + b\\frac{dx}{dt} + kx = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This ODE can be solved analytically and has been extensively studied. We will demonstrate how to numerically solve this ODE and general things to be kept in mind when solving for any model. Here m denotes the mass of the particle under oscillation, k denotes the stiffness of the spring, b denotes the damping coefficient of the external medium. x denotes the displacement from the mean position and t the time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many parameters are there in this equation. The answer is 2 and not 3, since you can divide by m and get $\\frac{b}{m}$ and $\\frac{k}{m}$ as new independent parameters. We will call them $2\\gamma$ and $\\epsilon$ respectively. Hence our ODE becomes\n",
    "$$\\frac{d^2x}{dt^2} + 2\\gamma\\frac{dx}{dt} + \\epsilon x = 0$$, the factor of 2 in $\\gamma$ is for the purpose of convenience. It is always recommended to simplify the ODE as much as possible. Theory of [Damped Simple Harmonic Oscillator(DSHO)](http://spiff.rit.edu/classes/phys312/workshops/w5b/damped_theory.html) predicts that different kinds of behaviour can be seen for different parameter values:\n",
    "* Underdamped for $\\epsilon > \\gamma^2$, the particle exhibits decaying oscillatory motion around mean position.\n",
    "* Overdamped for $\\epsilon < \\gamma^2$, If particle is released from the extreme position, it will slowly come towards the mean position without overshooting it.\n",
    "* Critically damped for $\\epsilon = \\gamma^2$, the boundary between underdamped and overdamped."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned earlier, we need to break this equation into 2 first order ODEs. This can be done as follows:\n",
    "Let $y = \\frac{dx}{dt}$\n",
    "Then the ODE becomes\n",
    "$$\\frac{dy}{dt} + 2\\gamma y + \\epsilon x = 0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we get the second equation? Note that we already have it when we defined $y = \\frac{dx}{dt}$. Note that y also represents the velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus finally, the system of first order ODEs that we need to simultaneously solve is:\n",
    "$$\\frac{dx}{dt} = y$$ and $$\\frac{dy}{dt} = -2\\gamma y -\\epsilon x$$\n",
    "Now our model is ready to be solved. We follow almost the same procedure as we followed previously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We refer to the pair X = [x,y] at any time as the state of the system. So the differential equation is finally of the form $$\\frac{dX}{dt} = f(X)$$, where X is a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shm_derivative(t, state, gamma, eps):\n",
    "    #time is the independent variable, state is the dependent variable\n",
    "    x = state[0]\n",
    "    y = state[1]\n",
    "    return [y, -2*gamma*y - eps*x] #We return the list of derivatives of each individual ODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_span = (0,10) #t0 and tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0,10,1000) #evaluation time points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_state = [1,0] #The initial state of the system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have 2 ODEs, we require 2 initial conditions. The first element of the list is the initial codition for our first ODE, i.e. x(t0) = 1, and the second element of the list is the initial condition for our second ODE, i.e. y(t0) = 0. What does this mean physically in reference to our model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Setting the parameter values\n",
    "#Note that eps > gamma^2, thus this deontes underdamped oscillations\n",
    "gamma = 0.5\n",
    "eps = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets call the solver and integrate!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve_ivp(shm_derivative, t_span, init_state, t_eval = t, args = (gamma, eps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 1000)\n"
     ]
    }
   ],
   "source": [
    "#obtain the solution from sol object\n",
    "state = sol.y\n",
    "print(state.shape)\n",
    "#Uncomment below to see the solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus state is a (2,1000) array. The first row denotes the value of the first state variable (x) at the required time points and the second row denotes the value of the second state variable(y) at the required time points. Note that for our system, x denotes the displacement of the particle from mean position and y simply denotes the velocity of the particle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#access x\n",
    "x_underdamped = state[0]\n",
    "#access y (the velocity)\n",
    "y_underdamped = state[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAGDCAYAAABwRoerAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxcdb3/8ddnJluzN1uXtGmbdF/oFtqmVCibsggVRAWUTRTxqqhX71WvXkV+Xhf0XnEBhSubiIAiS71UkK2FQoGmG93bdE/bNEnTJk3TrPP9/TFTCCG0aZuZMzN5Px+PPDLLyZx3Mm3z7jnf7/eYcw4RERERiS4+rwOIiIiIyPuppImIiIhEIZU0ERERkSikkiYiIiIShVTSRERERKKQSpqIiIhIFFJJE5FTYmbXm9niXny9B8zsR731euFiZgvN7HMR3ueHzGxjp/vbzey80O1bzexPp/Daa81sbi/EFJFeopIm0seYmTOzkV0eO6Vf8PJ+ZjbezOabWb2ZHTKzl81s9qm8pnPuVefcmF7I9r4i7Jyb4JxbeKqvLSK9RyVNRDxjZn6vM4SDmZUArwGrgRHAYOBJ4J9mVuZlNhGJHSppIvIeZjbXzCrN7BtmVm1me83shk7P54aOEDWY2VtASZevH2tmz5tZnZltNLNPdnruATP7nZktMLPDwNlmNtXMloeONj0GpHTavr+Z/Z+Z1ZjZgdDtIZ2eX2hmPzKz182s0cz+Hsr3cCjfUjMb3ml7Z2a3mNlWM6s1s5+bma/T8581s/WhfT1nZsM6PXe+mW0IHRn7LWDH+DHeCixxzn3XOVfnnDvknPs18BDws9DrpZjZn8xsv5kdDGUdEHoux8zuN7M9oSxPdX5vevg+/tXMqkJ5XzGzCaHHbwI+Dfz70Z9Z6PHOp06TzeyO0P73hG4nd85wjD8fF5nZutD7udvMvtmTvCLyfippItKdgUAWUAjcCNxpZv1Dz90JNAODgM+GPgAwszTgeeDPQAFwJXCXmY3v9NpXA/8FZABvAU8RLC85wF+Bj3fa1gfcDwwDioAjwG+7ZL0SuCaUtQRYEvqaHGA98IMu218GlALTgHlH85vZPOA/gMuBfOBV4JHQc3nAE8D3gDxgC3DGB/zsAM4PfS9d/QU4w8z6AdcR/BkPBXKBm0PfH6GfRyowgeDP8ZfH2NcH+QcwKvT1y4GHAZxz94Ru3+6cS3fOXdLN134XmAVMASYDMwh+70cd68/HvcAXnHMZwETgpZPILiKopIlI99qA25xzbc65BUAjMCZ0evLjwPedc4edc2uABzt93UeB7c65+51z7c65FcDfgE902uZp59xrzrkAwRKQCNwR2tfjwNKjGzrn9jvn/uaca3LOHSJY7s7qkvV+59wW51w9wWKyxTn3gnOunWBRmtpl+5+Fjm7tBO4Argo9fjPwE+fc+tDX/hiYEjqadhGw1jn3uHOuLfR1Vcf4+eUBe7t5fC/Bf3dzCP6Mc4GRzrkO59wy51yDmQ0CLgRuds4dCP1cFh1jX91yzt0XOoLXQvDI3mQzy+rhl3+a4Ptf7ZyrAX5IsAgf1e2fj07PjTezzFD+5SeaXUSCVNJE+p4OgsWos0SCv1yP2h8qKkc1AekEjzAlALs6Pbej0+1hwMzQ6buDZnaQ4C/8gZ226fy1g4HdzjnX3euZWaqZ3W1mO8ysAXgFyO4ylm1fp9tHurmf3uV77Zp9cKfsv+qUu47gKc3C0DbvfF0ob+fX6aqW4JHGrgYBAeAAwaNlzwGPhk4p3m5miQSPrNU55w4c4/WPycz8ZvZTM9sS+rltDz2V18OXGMx739fOPyf44D8fECzxFwE7zGyRaQyeyElTSRPpe3YCw7s8NoL3/lL+IDVAO8EicVRRp9u7gEXOuexOH+nOuS922qZzIdsLFJpZ5/FdnV/vGwSP0Mx0zmUCZ4YeP9Z4sOPpmn1Pp+xf6JK9n3Pu9VDOd74ulLfz63T1Au89enjUJwmOVWsKHYX6oXNuPDCb4FHIa0M5csws+2S/QYKnlOcB5xE8LTn8aPTQZ9fN13S2h2BpParzz+mYnHNLnXPzCJ5mfYrgKV4ROQkqaSJ9z2PA98xsiJn5QoPFLwEeP94XOuc6CI7NujV0lGs8wbFVR/0fMNrMrjGzxNDH6WY27gNecgnB0ndLaNvLCY5/OiqD4NGwg2aWw/vHl52MfwtNSBgKfJXgzwPg98B3Og2wzzKzo0XrGWCCmV1uZgnALbz36GBXPwRmm9l/hSYBZJjZVwiWsG+FXv9sM5sUOirYQPBIZsA5t5fgadu7QjkTzezMD9rRB8gAWoD9BMe2/bjL8/uA4mN8/SME/4zkh8bjfR847hItZpZkZp82s6zQaeEGgkcOReQkqKSJ9D23Aa8Diwmedrsd+HRofFlPfJngqa0q4AGCg/QBCI0b+zDBwfx7Qtv8DEju7oWcc60EB+pfT/D04qcIlsCj7gD6ETx9+AbwbA8zHsvTwDJgJcHydW8oy5OhrI+GThGuITg2DOdcLcEjYz8lWHxGEVxio1vOuc3AHIKD7rcTPBL3ceAjzrmjXzeQYDFuIDjBYRHBU6AQHP/VBmwAqoGvneD3+EeCR0Z3A+sI/uw6u5fguLGDR2eOdvEjoBx4m+AyIstDj/XENcD20M/wZoKnu0XkJNh7h4KIiMQvM3PAKOdchddZRESOR0fSRERERKKQSpqIiIhIFNLpThEREZEopCNpIiIiIlFIJU1EREQkCiV4HaC35eXlueHDh3sdQ0REROS4li1bVuucy+/uubgracOHD6e8vNzrGCIiIiLHZWYfeLUXne4UERERiUIqaSIiIiJRSCVNREREJAqppImIiIhEIZU0ERERkSikkiYiIiIShVTSRERERKKQSpqIiIhIFFJJExEREYlCnpY0M7vPzKrNbM0HPG9m9mszqzCzt81sWqQzioiIiHjB6yNpDwAXHOP5C4FRoY+bgN9FIJOIiIiI5zwtac65V4C6Y2wyD/ijC3oDyDazQZFJ98EWbaqhfHsdDc1tXkcRERGROBXtF1gvBHZ1ul8Zemxv543M7CaCR9ooKioKe6gvPbycxpZ2fAZnjc7nlnNHMbWof9j3KyIiIn2H16c7e4Vz7h7nXKlzrjQ/Pz/s+3vsC7O4//rTufmsElbvbuDy373OD/++ltb2QNj3LSIiIn1DtB9J2w0M7XR/SOgxT00YnMWEwXD22AL+5eyR3P7sBu5/bTub9h3i95+ZTkZKotcRRUREJMZF+5G0+cC1oVmes4B659ze431RJKUnJ3DbvIn84hOTeXNrHdfd9xZHWju8jiUiIiIxzuslOB4BlgBjzKzSzG40s5vN7ObQJguArUAF8L/Av3gU9biumD6E31w1lRW7DnLLoyvoCDivI4mIiEgM8/R0p3PuquM874AvRSjOKbtw0iC+/9Hx/PDv67jz5QpuOXeU15FEREQkRkX76c6Yc/3s4XxsymDueGETb27d73UcERERiVEqab3MzPjRZZMoyknlG39dRVNru9eRREREJAappIVBenICt18xmcoDR/jl85u8jiMiIiIxSCUtTGaMyOGqGUXcu3gb6/Y0eB1HREREYoxKWhh9+4KxZPZL5McL1hOcAyEiIiLSMyppYZSVmshXzhnF4opaFm2q8TqOiIiIxBCVtDC7ZtYwinJS+cmCDVo7TURERHpMJS3MkhJ8/PsFY9i47xBPr/T8ilYiIiISI1TSIuCiiYMYMyCDuxZuIaCjaSIiItIDKmkR4PMZXzpnJBXVjTy3tsrrOCIiIhIDVNIi5OJJgxiRl8ZvXqrQTE8RERE5LpW0CPH7jC/OLWHd3gYWaqaniIiIHIdKWgR9bEoh+RnJ3P/adq+jiIiISJRTSYugpAQfn5k5jFc21VBR3eh1HBEREYliKmkRdvXMIpL8Ph58fbvXUURERCSKqaRFWH5GMpdMHszflldSf6TN6zgiIiISpVTSPHDDGcNpau3gieWVXkcRERGRKKWS5oGJhVlMKszisaW7tByHiIiIdEslzSOfOn0oG6oO8XZlvddRREREJAqppHnk0imDSUn08ejSXV5HERERkSikkuaRzJRELp40mL+v2kNTa7vXcURERCTKqKR56MoZQ2lsaeeZt/d6HUVERESijEqah0qH9ac4P43HdMpTREREulBJ85CZccX0IZTvOMCuuiav44iIiEgUUUnz2KWTBwPw9MrdHicRERGRaKKS5rEh/VM5fXh/nlq5R2umiYiIyDtU0qLAvCmFVFQ3sm5vg9dRREREJEqopEWBiyYNIsFnPL1yj9dRREREJEqopEWBnLQkzhqdz/yVe+gI6JSniIiIqKRFjXlTC6lqaObNbfu9jiIiIiJRQCUtSpw3roB+iX4WrNbCtiIiIqKSFjVSkxKYOyaf59buI6BTniIiIn2eSloUuWDiQGoOtbBs5wGvo4iIiIjHVNKiyDljC0hK8PGP1VVeRxERERGPqaRFkYyURM4clceza/ZqYVsREZE+TiUtylwwcRB76ptZVVnvdRQRERHxkKclzcwuMLONZlZhZt/u5vkiM3vZzFaY2dtmdpEXOSPp/HEDSPAZ/1ijWZ4iIiJ9mWclzcz8wJ3AhcB44CozG99ls+8Bf3HOTQWuBO6KbMrIy0pNZPbIPJ5dU6VTniIiIn2Yl0fSZgAVzrmtzrlW4FFgXpdtHJAZup0F9InrJl04cSA79jexfu8hr6OIiIiIR7wsaYXArk73K0OPdXYr8BkzqwQWAF+JTDRvfXj8AHwGz67VLE8REZG+KtonDlwFPOCcGwJcBDxkZu/LbGY3mVm5mZXX1NREPGRvy01PZlpRf15cv8/rKCIiIuIRL0vabmBop/tDQo91diPwFwDn3BIgBcjr+kLOuXucc6XOudL8/PwwxY2s88YPYO2eBvbWH/E6ioiIiHjAy5K2FBhlZiPMLIngxID5XbbZCZwLYGbjCJa02D9U1gPnjSsA4MX11R4nERERES94VtKcc+3Al4HngPUEZ3GuNbPbzOzS0GbfAD5vZquAR4DrXR+Z8liSn86w3FSd8hQREemjErzcuXNuAcEJAZ0f+36n2+uAMyKdKxqYGeeOHcCf3txBU2s7qUmevlUiIiISYdE+caBPO29cAa3tAV7dXOt1FBEREYkwlbQodvqIHDJSEnTKU0REpA9SSYtiiX4fc8cU8NKGGgKBPjEUT0REREJU0qLceeMKqG1sYVXlQa+jiIiISASppEW5uaML8PtMS3GIiIj0MSppUS4rNZFpRdks2tQnlocTERGREJW0GHDW6HxW766ntrHF6ygiIiISISppMWDumODVB17R0TQREZE+QyUtBowflEleepJOeYqIiPQhKmkxwOczzhydzyubaujQUhwiIiJ9gkpajDhrdD4HmtpYvbve6ygiIiISASppMeJDo/Ixg0UbdcpTRESkL1BJixE5aUlMHpLNwk1aL01ERKQvUEmLIWeNzmfVroMcONzqdRQREREJM5W0GHLWmHwCDl6tqPU6ioiIiISZSloMmTwkm+zURI1LExER6QNU0mKI32d8aFQ+izbVENBSHCIiInFNJS3GnDU6n9rGFjbuO+R1FBEREQkjlbQYM2dkHgCLN2tcmoiISDxTSYsxA7NSGFmQzmJNHhAREYlrKmkxaM7IPN7ctp+W9g6vo4iIiEiYqKTFoDkj82huC7B8x0Gvo4iIiEiYqKTFoJnFOfh9xuIKLcUhIiISr1TSYlBGSiJThmazuGK/11FEREQkTFTSYtSckXmsrjxIfVOb11FEREQkDFTSYtScUXkEHCzZqlmeIiIi8UglLUZNGZpNWpJfS3GIiIjEKZW0GJXo9zGrOFeL2oqIiMQplbQYNmdUHtv3N7GrrsnrKCIiItLLVNJi2NFLRL2mU54iIiJxRyUtho0sSGdAZrLGpYmIiMQhlbQYZmacMTKP17fsJxBwXscRERGRXqSSFuPmjMyj7nAr66savI4iIiIivUglLcaVleQCsGSLrj4gIiIST1TSYtygrH6MyEtTSRMREYkzKmlxoKwklze31dHeEfA6ioiIiPQSlbQ4UFacS2NLO2v2aFyaiIhIvPC0pJnZBWa20cwqzOzbH7DNJ81snZmtNbM/RzpjLJhVHByX9voWLcUhIiISLzwraWbmB+4ELgTGA1eZ2fgu24wCvgOc4ZybAHwt4kFjQH5GMmMGZGhcmoiISBzx8kjaDKDCObfVOdcKPArM67LN54E7nXMHAJxz1RHOGDPKSnIp336A1naNSxMREYkHXpa0QmBXp/uVocc6Gw2MNrPXzOwNM7uguxcys5vMrNzMymtqasIUN7qVleRypK2DVZUHvY4iIiIivSDaJw4kAKOAucBVwP+aWXbXjZxz9zjnSp1zpfn5+RGOGB1mjcjFDF6v0ClPERGReOBlSdsNDO10f0josc4qgfnOuTbn3DZgE8HSJl1kpSYyYXAmS7Zq8oCIiEg88LKkLQVGmdkIM0sCrgTmd9nmKYJH0TCzPIKnP7dGMmQsKSvOZfmOgzS3dXgdRURERE6RZyXNOdcOfBl4DlgP/MU5t9bMbjOzS0ObPQfsN7N1wMvAvznndD7vA8wuyaO1I8DyHQe8jiIiIiKnKMHLnTvnFgALujz2/U63HfCvoQ85jtNH5OD3Ga9v2c/skXlexxEREZFTEO0TB+QEpCcncNqQLC1qKyIiEgdU0uJMWXEub1fW09jS7nUUEREROQUqaXFmdkke7QHH0u11XkcRERGRU6CSFmemD+tPot94Q5eIEhERiWkqaXGmX5KfqUX9eV0lTUREJKappMWhsuJc1u6pp76pzesoIiIicpJU0uLQ7JJcAg7e3KajaSIiIrFKJS0OTSnKJiXRp1OeIiIiMUwlLQ4lJ/gpHZbDG1tV0kRERGKVSlqcKivJZUPVIeoOt3odRURERE6CSlqcmlWcA8CbOpomIiISk1TS4tRpQ7JJTfKzRCVNREQkJqmkxalEv4/S4Tks0eQBERGRmKSSFsfKinPZXN1IbWOL11FERETkBKmkxbGj49I0y1NERCT2qKTFsUmFWaQnJ+iUp4iISAxSSYtjCX4fpw/vr8kDIiIiMUglLc6VleSyteYw+xqavY4iIiIiJ0AlLc6VFecBGpcmIiISa1TS4tz4wZlkpCSopImIiMSYHpU0M/tqTx6T6OP3GTNHaL00ERGRWNPTI2nXdfPY9b2YQ8JoVnEu2/c3sbf+iNdRREREpIcSjvWkmV0FXA2MMLP5nZ7KAOrCGUx6T1lJLhAcl3bZ1CEepxEREZGeOGZJA14H9gJ5wH93evwQ8Ha4QknvGjcwk6x+iSzZopImIiISK45Z0pxzO4AdQFlk4kg4+I6OS9PkARERkZjR04kDl5vZZjOrN7MGMztkZg3hDie9p6wkl111R6g80OR1FBEREemBnk4cuB241DmX5ZzLdM5lOOcywxlMete749I0lFBERCQW9LSk7XPOrQ9rEgmr0QUZ5KQlaSkOERGRGHG8iQNHlZvZY8BTQMvRB51zT4QllfS6o+PS3ti6H+ccZuZ1JBERETmGnh5JywSagA8Dl4Q+PhquUBIeZSW57D54hF11Wi9NREQk2vXoSJpz7oZwB5HwKysOjktbsrWWotwij9OIiIjIsfR0dudoM3vRzNaE7p9mZt8LbzTpbSML0slLT9LkARERkRjQ09Od/wt8B2gDcM69DVwZrlASHmbGzOJclmwJjksTERGR6NXTkpbqnHury2PtvR1Gwq+sOJeqhma279d6aSIiItGspyWt1sxKAAdgZlcQvFyUxJij66VpKQ4REZHo1tOS9iXgbmCsme0GvgZ8MWypJGyK89IoyEjmDV0iSkREJKr1qKQ557Y6584D8oGxzrk5zrntp7pzM7vAzDaaWYWZffsY233czJyZlZ7qPvs6M6OsJJclWzUuTUREJJr1aAkOM8sGrgWGAwlHF0J1zt1ysjs2Mz9wJ3A+UAksNbP5zrl1XbbLAL4KvHmy+5L3mlWcy9Mr97Cl5jAjC9K9jiMiIiLd6OnpzgUEC9pqYFmnj1MxA6gIHaVrBR4F5nWz3f8DfgY0n+L+JOTd9dJ0ylNERCRa9fSyUCnOuX/t5X0XArs63a8EZnbewMymAUOdc8+Y2b/18v77rGG5qQzKSuGNrfu5ZtYwr+OIiIhIN3p6JO0hM/u8mQ0ys5yjH+EMZmY+4H+Ab/Rg25vMrNzMymtqasIZKy6YGWXFubypcWkiIiJRq6clrRX4ObCEd091lp/ivncDQzvdHxJ67KgMYCKw0My2A7OA+d1NHnDO3eOcK3XOlebn559irL5hVnEutY2tbK5u9DqKiIiIdKOnpzu/AYx0ztX24r6XAqPMbATBcnYlcPXRJ51z9UDe0ftmthD4pnPuVMuh8N710kYPyPA4jYiIiHTV0yNpFUCvLlHvnGsHvgw8B6wH/uKcW2tmt5nZpb25L3m/oTmpFGb306K2IiIiUaqnR9IOAyvN7GWg5eiDp7IER+jrFxCcOdr5se9/wLZzT2Vf8n5lJbm8uH4fgYDD5zOv44iIiEgnPS1pT4U+JI6UFefy+LJKNu47xLhBmV7HERERkU56VNKccw+aWT+gyDm3McyZJEJmdRqXppImIiISXXo0Js3MLgFWAs+G7k8xs/nhDCbhV5jdj6KcVC1qKyIiEoV6OnHgVoJXCDgI4JxbCRSHKZNEUFlxLm9tqyMQ0HppIiIi0aSnJa0ttCRGZ4HeDiORV1aSS/2RNtbtbfA6ioiIiHTS05K21syuBvxmNsrMfgO8HsZcEiGzQtfxfEOnPEVERKJKT0vaV4AJBJff+DNQD3wtXKEkcgZmpTAiL03rpYmIiESZns7ubAK+G/qQODOrOJf/W7WHjoDDr/XSREREokJPZ3c+b2bZne73N7PnwhdLIqmsJJdDLe2s3dN12KGIiIh4paenO/OccweP3nHOHQAKwhNJIm1WcQ6ATnmKiIhEkZ6WtICZFR29Y2bDAK3ZECcKMlIoyU/TemkiIiJRpKeXhfousNjMFgEGfAi4KWypJOLKSnJ5cvlu2joCJPp72t1FREQkXHr029g59ywwDXgMeBSY7pzTmLQ4Ulacx+HWDtbs1rg0ERGRaHAih0w6gGqgARhvZmeGJ5J4YebRcWk65SkiIhIVejq783PAK8BzwA9Dn28NXyyJtLz0ZEYPSNfkARERkSjR0yNpXwVOB3Y4584GphK6jqfEj7LiXMq3H6C1XVf8EhER8VpPS1qzc64ZwMySnXMbgDHhiyVeKCvJ5UhbB6t3q3+LiIh4raclrTK0mO1TwPNm9jSwI3yxxAszRwSv46lTniIiIt7r6ezOy5xzB51ztwL/CdwLfCycwSTy+qclMXZghiYPiIiIRIFjrpNmZjndPLw69DkdqOv1ROKpspJc/vzmTlraO0hO8HsdR0REpM863mK2ywheWaC7q247oLjXE4mnyopzuf+17azaVc+MEd11dBEREYmEY5Y059yISAWR6DBzRC5mwXFpKmkiIiLe6fFitmZ2uZn9j5n9t5lpPFqcykpNZMLgTJZsrfU6ioiISJ/W08Vs7wJuJjgebQ1ws5ndGc5g4p1ZI3JZvvMgzW0dXkcRERHps3p6JO0c4CPOufudc/cDF4UekzhUVpJLa3uA5TsPeB1FRESkz+ppSasAijrdHxp6TOLQ6SNy8Bm8sVWTd0VERLzS05KWAaw3s4VmthBYB2Sa2Xwzmx+2dOKJzJREJhVm8YYWtRUREfHM8ZbgOOr7YU0hUWdWcS73vbaNI60d9EvSemkiIiKR1tMjaeXAq865RcBeIAt43Tm3KPSYxJlZJbm0dTiW7dC4NBERES/0tKS9AqSYWSHwT+Aa4IFwhRLvnT48B7/PeEOXiBIREfFET0uaOeeagMuBu5xznwAmhi+WeC09OYHThmTpOp4iIiIe6XFJM7My4NPAMyf4tRKjyopzWbXrII0t7V5HERER6XN6WrS+BnwHeNI5t9bMioGXwxdLosGckXm0B5xmeYqIiHigR7M7Q5MDFnW6vxW4JVyhJDpMH96flEQfiytqOW/8AK/jiIiI9CnHLGlmdodz7mtm9nfAdX3eOXdp2JKJ55IT/MwYkcurm2u8jiIiItLnHO9I2kOhz78IdxCJTmeOyuNHz6xnb/0RBmX18zqOiIhIn3HMkuacWxb6vMjM8kO3dVilD5kzKg+AVzfX8snSoR6nERER6TuOO3HAzG41s1pgI7DJzGrMrFeuQGBmF5jZRjOrMLNvd/P8v5rZOjN728xeNLNhvbFf6bkxAzLIz0hm8eZar6OIiIj0KccsaWb2r8AZwOnOuRznXH9gJnCGmX39VHZsZn7gTuBCYDxwlZmN77LZCqDUOXca8Dhw+6nsU06cmTFnZB6vVdQSCLxvWKKIiIiEyfGOpF0DXOWc23b0gdDMzs8A157ivmcAFc65rc65VuBRYF7nDZxzL4cW0QV4AxhyivuUkzBnZB77D7eyvqrB6ygiIiJ9xvFKWqJz7n3nuULj0hJPcd+FwK5O9ytDj32QG4F/dPeEmd1kZuVmVl5ToyFzve3ouDSd8hQREYmc45W01pN8rleZ2WeAUuDn3T3vnLvHOVfqnCvNz8+PVKw+Y0BmCqMHpLO4QiVNREQkUo63BMdkM+vuHJcBKae4791A5+mCQ0KPvXdHZucB3wXOcs61nOI+5STNGZnPw2/uoLmtg5REv9dxRERE4t4xj6Q55/zOucxuPjKcc6d6unMpMMrMRphZEnAlML/zBmY2FbgbuNQ5V32K+5NT8KFRebS0ByjffsDrKCIiIn2CZxdJd861A18GngPWA38JXRf0NjM7eiWDnwPpwF/NbKWZzf+Al5Mwm1mcQ6LfeLVCY/5EREQioUfX7gwX59wCYEGXx77f6fZ5EQ8l3UpNSmBaUf/g5IELvU4jIiIS/zw7kiax50Oj8li7p4H9jRoaKCIiEm4qadJjc0YFZ85qlqeIiEj4qaRJj00qzCKrX6LWSxMREYkAlTTpMb/PmF2Sy+KKWpzTJaJERETCSSVNTsiHRuWzt76ZiupGr6OIiIjENZU0OSFnjVAAfgoAACAASURBVAmOS1u0SUtxiIiIhJNKmpyQwux+jCpIZ+FGlTQREZFwUkmTEzZ3TD5vbavjcEu711FERETilkqanLC5Ywpo7QiwZMt+r6OIiIjELZU0OWGlw/uTmuRn4SZdTlVERCRcVNLkhCUn+JldksvCjTVaikNERCRMVNLkpJw1poDKA0fYWnvY6ygiIiJxSSVNTsrc0cGlODTLU0REJDxU0uSkDM1JpSQ/jYUbNS5NREQkHFTS5KTNHVPAm9vqONLa4XUUERGRuKOSJidt7ph8WtsDvLFVS3GIiIj0NpU0OWmnD8+hX6JfpzxFRETCIMHrABK7UhL9lJXkslDX8ZQY0Noe4EhbB+0dAToCjg7nSE7w0y/RT3KCD5/PvI4oIvIeKmlySuaOyeelDdVsrWmkOD/d6zjSh7W2B9hYdYjN1YfYXnuYbfubqDzQxP7GVg4cbuXQcS5jlpGSQH5GMgUZyeRnpDA4K4UReWkU56dTnJ9GbloSZipyIhI5KmlySs4ZW8D3n17Li+urVdIkohpb2nm9opbXt+xnVeVB1u5poLU9AIDPoLB/P4pyUikqyqZ/ahI5aUmkJvlJ9PtI8Bs+M1raOjjSFjzCVt/USk1jCzWHWlhdeZDn1jTT2hF4Z3/9UxOZWJjFaUOymFSYzWlDshiUlaLiJiJho5Imp2RI/1TGDszghfX7+PyZxV7HkTi35+ARnnl7Ly9tqKZ8Rx1tHY6URB+nFWZzXdkwJg/NZuzATIbm9CM5wX9K++oIOPYcPMKWmka21hxmY9UhVu+u5+5FW2kPBK+0MTgrhVnFucwqyaWsOJehOam98W2KiAAqadILzhs3gN8t2sLBplayU5O8jiNx5lBzG8+8vZenVu7mzW11OAdjB2bw2TkjOGt0PqXDckhK6P05UH6fMTQnlaE5qcwd8+7jzW0drN/bwNuV9by1rY5Fm2p4YsVuAAqz+/GhUXmcM7aAOaPySE3SP7EicvIs3q69WFpa6srLy72O0aes3HWQj935Gnd8agofm1rodRyJE1trGnnw9e08vqySw60dFOelMW9KIfOmDGZ4XprX8d7hnGNzdSNLtuzn9S21vFaxn8aWdpISfJQV53LuuALOGVvAkP46yiYi72dmy5xzpd0+p5ImpyoQcMz8yYvMHJHDb6+e5nUciXErdh7gNy9V8NKGahL9xiWnDeaasmFMGZodE+O/WtsDLN1ex4vrq3lpwz62728C4LQhWVw0aRAXTxqk06Ii8g6VNAm7bz3+NgvW7GXZ984Py6kniX+rK+v55QubeGlDNTlpSVxbNoyrZxZRkJHidbRTsrWmkX+u28eC1Xt5u7IegMmhwnaRCptIn6eSJmH3z7VV3PTQMh7+3EzOGJnndRyJIZUHmvjxgvUsWF1FdmoiN51ZzHVlw0lLjr/xXDv3N7Fgzd73FLYpQ7O5fFohHz1tMDlpGtMp0teopEnYNbW2M/W257l6ZhE/uGSC13EkBhxp7eD3i7bw+0Vb8JnxhbOKuXHOCDJSEr2OFhE79zfxzOq9PL1yNxuqDpHgM+aOKeDyaYWcM7aAlMRTm50qIrFBJU0i4rMPLKWiupFF/zY3JsYOiXde3ljN955cw+6DR7hk8mC+c+FYBmf38zqWZ9bvbeDJFbt5asVuqg+1kJGSwMWTBnHZ1EJOH56jqyGIxLFjlbT4O58gnjl3XAEvbaimorqRUQMyvI4jUai+qY3/98w6Hl9WyaiCdB67aRYzi3O9juW5cYMyGTcok29dMJbXt9Ty5PLdzF+1h0eX7qIwux+XTS3ksmmFlGjBaJE+RUfSpNdU1Tcz6ycv8u8XjOFf5o70Oo5EmZc3VPPtJ96mtrGVm88q5pZzR53ygrPxrKm1nefWVvHkij0s3lxDwAUnHFw2tZBLJg8mNz3Z64gi0gt0ulMi5pLfLCYpwcffvjjb6ygSJVrbA/z0Hxu477VtjBmQwS8+MZlJQ7K8jhVTqhuamb9qD39bvpv1extC49fyuWzqEM4dp/FrIrFMpzslYs4dV8CvXtxMbWMLefqffp+3vfYwX3lkBat313Nd2TC+c9E4FYqTUJCZwuc+VMznPlTMhqoGnly+m6dW7uaF9dUavyYSx3QkTXrVuj0NXPTrV/np5ZO4ckaR13HEQ8+uqeKbf12Fz+D2KyZzwcSBXkeKKx0Bx5It+3lieSXPrq2iqbWDIf1D49emFlKs8WsiMUGnOyVinHOc9fOFjMhL48HPzvA6jnggEHD8+qXN3PHCZiYPyeLOT0/TJZHC7Oj4tSeW7+a1itrg+LWh2VweGr+m9ddEopdKmkTUTxas577XtlH+vfPJ6tc31rySoMMt7XzjL6t4dm0Vl08r5MeXTdLpzQjb19DM/JV7+Nvyyk7rrwXHr50ztoB+SXo/RKKJxqRJRH1k4kDufmUrL23Yx2VTh3gdRyJkV10Tn/9jOZv2HeJ7F4/jxjkjtF6eBwZkpvD5M4v5/JnF71l/7YX11fRL9HPOuAIumjiIs8fmk5qkXwEi0czTv6FmdgHwK8AP/ME599MuzycDfwSmA/uBTznntkc6p5yYKUOyGZiZwj9WV6mk9RFr99Rz/f1LaWnr4IEbZnDm6HyvIwnvXX/tza37eWb1Xp5bW8Uzb++lX6Kfs8fmc9GkQZw9piAuL8MlEus8+1tpZn7gTuB8oBJYambznXPrOm12I3DAOTfSzK4EfgZ8KvJp5UT4fMYFEwfyyFs7aWpt1//W49xrFbV84aFlZKQk8PAXZzNaCxlHHb/PmD0yj9kj87ht3kTe2lbHgtV7+ceaKhasriIl0ceZo/I5d1wBZ48poCAzti9qL+HlnKOptYP6I200NLdxuKWdptYODrd0cKStPfi5tYPDre0cae0IPtfaTnNbB63tAVo7HG3tAdo6gh+tHY7W9g7aOtw7j7UHHIGAwwE4cKH9Bhw4HC70GKH7ZobfDL/P8Bkk+H34zPD7CD7uDz7v8727XaLfR6I/+DkpwUeS3xd8LHQ7KcEYmNmPr543yrOftZe/PWcAFc65rQBm9igwD+hc0uYBt4ZuPw781szMxdtAujj0kQkDeeD17SzcWMNFkwZ5HUfC5OmVu/nmX1dRnJfOA589nUFZfffSTrHC7zPKSnIpK8nl1ksnsHR7sLC9sG4f/1y3DwgumnvO2AGcO66ACYMzddo6jrV3BKhramV/Y+jjcAu1ja0cONxK/ZG293w0HP3c3EZbR89+DScn+EhN8pOalEC/JH+wCCX4SPIbyYk+0lMSgiUpVJiSEoJFKVi2gn/uzMAIli8zMDMMwMAXuh1w0BEI0BGAgHN0BNw7Ra/DBT+3d70dCJbClvYAjS3twcLYHqCtw4U+ByjKSe2zJa0Q2NXpfiUw84O2cc61m1k9kAvURiShnLQZI3LITUviH2uqVNLi1P2vbeOHf1/HzBE53HNtqSaJxCC/z5hVnMus4lx+eOkENlQd4qUN1by4fh93vLiJX76wiYKMZOaMyuOMkjzOGJnHwCwdZYt2zjkamtvZ19DM3vpmquqPsK+hhf2NwQJW29jC/sOt7G9s4UBTW7ev4TPI7JdIVqePwv793rmdmfLu42nJftKSE+iXGPwcLGV++iX6SfD7Ivzdx5e4OA9lZjcBNwEUFWltrmjg9xkfnjCA+Sv30NzWoRl+ceauhRXc/uxGPjJhAL+6cqre3zhgZu+MYfvS2SPZ39jCwo01vLSxmpc3VPPE8t0AFOenMbsklzNK8jh9RI4WrY6wQMBR19RKVX0zVfXN7G0IlrC99c2dSlkzTa0d7/vazJQE8tKTyUtPZlRBOrOKc8hNSyYvPYnc9GRy04Kf89KTyExJ1MLIUcDLkrYbGNrp/pDQY91tU2lmCUAWwQkE7+Gcuwe4B4JLcIQlrZywj0wYyCNv7eK1ilrOHTfA6zjSC5xz3PHCZn714mYunTyY//nkZP1POU7lpifz8elD+Pj0IQQCjvVVDbxesZ/XttTyxPLd/OmNnQAMy01lWlF/phVlM7WoP2MHZujPxElq7whQ09jyTtGqqm+mqtPRsKqGZvbVt9DaEXjP1/l9RkFGMgOzUhg7MIO5owsYmJXMwKx+DMpKYWBmCgWZybpWbgzysqQtBUaZ2QiCZexK4Oou28wHrgOWAFcAL2k8WuyYXZJHRkoCC1ZXqaTFAeccP312A3cv2sonpg/hpx8/Db/+p90n+HzGhMFZTBicxefPLKa1PcDblQdZtuMAy3ceYHFFLU+uCP4fOzXJz/jQEbnxg4OfxwzI6PPrsx1p7aCqIVi89jU0v3O789GwmkMtBLr8hktK8L1TtKYV9WdgVgqDMlMYmNUveDsrhbz0ZP1djFOelbTQGLMvA88RXILjPufcWjO7DSh3zs0H7gUeMrMKoI5gkZMYkZTg48PjB/LPdVW0tE/U/+JimHOO/3pmPX9YvI3PzCritksn6lRIH5aU4KN0eA6lw3OA4J+PygNHWL7zACt2HmTtnnqeXLGbh97YAQTHNw3PS2NUQTrD89IYkZsW/JyXRkFGckxPTGjrCFB3uJWaQy1UH2qmqr7lnaNeVQ0t7AsdDas/8v6xX+nJCe8UrdEF+QzKSmFA6P7AzOBRsOzUxJj++cip0RUHJKwWbarhuvve4u5rpvORCbp2Y6z6+XMbuPPlLVw/ezg/uGS8fmnIcQUCweK2bm8D60MfW2oa2VnX9J6ZgalJfob078eAzODRooFZKe/czs9Ifneger/EsB8tam0P0NAcnMV4qLk9dLudA6HZj8EB950G3ze2dlu+zCA/Pfk938vArC7fX1YK6VqbTtAVB8RDZ5TkkpOWxPxVe1TSYtRvXtzMnS9v4aoZQ1XQpMd8PqMoN5Wi3FQumPju3/2OgGPPwSNsqz3M9v2H2VZ7mN0HjrCvoZlN+w51e8rvqIyUBLL6JZKRkkhKoo/kBB/JCX6SE3ykJPpJDI2FC62uRedPR5daaGkP0NzWEbwd+ny4JVjImtsC3e84JDs1kdy0JPLSkxk3MJPc9ODt3PQkctOSKchMZlBWCvnpyRqXJ71CJU3CKsHv46JJA3l8WSWHW9q1qnmM+d9XtvLfz2/i8qmF/NfHJqmgySnz+4yhOakMzUnlTN5/ZYr2jgC1ja1UNTRTe6il27W6GprbaWkPFqyDR9reKVut7e+WrKN/VI9+TvT5SE70h4qdj+x+iSRnJJOSGFwuIrNfIpkpCaHPiWT2SyAzJVgIs1MTyUlLeqcEikSKfmNK2F06uZA/vbGT59ft42NTC72OIz300JLt/NeC9Vw8aRC3X3GaxqBJRCT4fcFTg1qPTQT9t0DCrnRYfwZlpTB/1R6vo0gPPbViN//59FrOGzeAO66colM3IiIe0L+8EnY+n3HJ5MG8sqmGA4dbvY4jx/HKphq++ddVzCrO4bdXT9UpHhERj+hfX4mISycPpj3g+MeaKq+jyDGs2nWQm/+0jFEDMrjn2lJdSUBExEMqaRIREwZnUpyXxtMru15UQqLFttrD3PDAUnLSknjwhtPJTNG1OEVEvKSSJhFhFjzl+db2Oqrqm72OI11UNzRz7X1vAvDHz86gIFODtkVEvKaSJhEzb8pgnIOndDQtqhxqbuO6+5eyv7GV+68/neL8dK8jiYgIKmkSQcX56UwryubxZZXE25UuYlVHwPHVR1eyad8hfveZ6Uwemu11JBERCVFJk4j6ROlQKqobWVVZ73UUAX727AZe2lDNrZeM56zR719YVEREvKOSJhF18WmDSE7w8fiyXV5H6fP+Wr6Le17ZyjWzhnFN2XCv44iISBcqaRJRmSmJXDBxIPNX7qG5rcPrOH3W0u11/MeTqzljZC7fv2S813FERKQbKmkScVdMH0JDczsvrN/ndZQ+aVddE194aBlD+qdy19XTtVitiEiU0r/OEnGzS/IYlJXC48sqvY7S5zS2tPO5B8tp7wjwh+tKyUrVWmgiItFKJU0izu8zPj5tCK9sqtGaaRHUEXB87dEVVNQ0ctenp1OipTZERKKaSpp44orpQwg4NIEggm5/bgMvrK/mB5eMZ86oPK/jiIjIcaikiSeG56VxxshcHnlrFx0BrZkWbo8vq+TuRVv5zKwirtVMThGRmKCSJp759Mxh7D54hIUbq72OEtfKt9fxH08EZ3L+4JIJXscREZEeUkkTz5w/fgD5Gck8/OZOr6PEraMzOQdnp3Dn1dM0k1NEJIboX2zxTKLfx5WnD+XljdVUHmjyOk7caWxp5/N/LKe1I8Afrjud7NQkryOJiMgJUEkTT105owgDHnlLR9N6UyDg+NqjK9lc3cidV09jZIFmcoqIxBqVNPFUYXY/zh5TwGNLK2ltD3gdJ27c/txGXli/j/+8eBxn6pqcIiIxSSVNPPeZWcOobWzh2bVVXkeJC39bVsnvF23h6plFXDd7uNdxRETkJKmkiefOHJ3P8NxU7l28Dee0HMepWLajju88sZrZJbn88NIJmJnXkURE5CSppInn/D7js3NGsGrXQZbvPOB1nJhVeeDdmZx3fVozOUVEYp3+FZeocMX0IWT1S+QPr27zOkpMOhy6JmdLu2ZyiojEC5U0iQqpSQlcPbOI59ZWsXO/luM4EYGA4+uPrWTTvkOaySkiEkdU0iRqXD97OH6fcf/rOpp2In7xz438c90+/vOj4zWTU0QkjqikSdQYkJnCJacN5i9Ld1F/pM3rODHhieWV3LVwC1fNKOJ6zeQUEYkrKmkSVT47ZwSHWzt4+M0dXkeJest2HODbf1vNrOIcbpunmZwiIvFGJU2iysTCLM4anc8fXt1GU2u713Gi1u6DR/jCQ+UMyk7hd5+erpmcIiJxSP+yS9S55dyR1B1u5c+68Hq3GlvaufGBpbS0Bbj3ulL6p2kmp4hIPFJJk6gzfVgOs0tyufuVrTS3dXgdJ6p0BBxffWQFm6sb+e2npzGyIMPrSCIiEiYqaRKVvnLOKGoOtfDY0l1eR4kqP16wnhc3VHPrJeM5SzM5RUTimkqaRKVZxTmcPrw/v1u4RUfTQh5+cwf3Lt7G9bOHc03ZcK/jiIhImKmkSVQyM75+/miqGpp5aIlmei7eXMv3n17L3DH5fO/icV7HERGRCPCkpJlZjpk9b2abQ5/7d7PNFDNbYmZrzextM/uUF1nFO7NL8jhzdD53LqygobnvrptWUd3IFx9exsj8dH5z1VQSNJNTRKRP8Opf+28DLzrnRgEvhu531QRc65ybAFwA3GFm2RHMKFHg3z8yhoNNbdy9aIvXUTxRd7iVGx9cSnKCjz9cV0pGSqLXkUREJEK8KmnzgAdDtx8EPtZ1A+fcJufc5tDtPUA1oJHSfczEwiwumTyYexdvo7qh2es4EdXc1sHn/1jO3vpm7r6mlKE5qV5HEhGRCPKqpA1wzu0N3a4CBhxrYzObASQB3R5OMbObzKzczMpramp6N6l47psfHk17h+MX/9zodZSI6Qg4bnlkBct3HuCOT01h+rD3jQgQEZE4F7aSZmYvmNmabj7mdd7OOecAd4zXGQQ8BNzgnAt0t41z7h7nXKlzrjQ/Xwfb4s2w3DQ+O2cEfymvZOWug17HCTvnHD+Yv4Z/rtvHDz46nosmDfI6koiIeCBsJc05d55zbmI3H08D+0Ll62gJq+7uNcwsE3gG+K5z7o1wZZXo95VzRpKfkcwPnl5DIPCBnT4u3LVwC396YydfOKuY688Y4XUcERHxiFenO+cD14VuXwc83XUDM0sCngT+6Jx7PILZJAplpCTyHxeNZVVlPX9dFr8L3P61fBc/f24jH5symG99ZKzXcURExENelbSfAueb2WbgvNB9zKzUzP4Q2uaTwJnA9Wa2MvQxxZu4Eg0+NqWQ0mH9+dmzG6k73Op1nF73/Lp9fOeJ1cwZmcftV0zG5zOvI4mIiIcsOCQsfpSWlrry8nKvY0iYbKhq4JLfLOaiSYP41ZVTvY7TaxZvruWzDyxl3OBM/nTjDC21ISLSR5jZMudcaXfPaVVMiSljB2bypbNH8vTKPTy/bp/XcXrF0u11fP6P5RTnp/HgDaeroImICKCSJjHoX+aOZOzADL775Grqm2L7SgRvVx7khvuXMig7hYdunEl2apLXkUREJEqopEnMSUrw8YtPTGb/4Va++9RqYvWU/Zrd9Vx731tkpyby8Odmkp+R7HUkERGJIippEpMmFmbxr+eP5v/e3ssjb8XebM8VOw9w9f++QVpSAn/+3CwGZfXzOpKIiEQZlTSJWV88q4QPjcrj1r+vZd2eBq/j9Nhb2+r4zB/epH9aEo99YRZFubrck4iIvJ9KmsQsn8/45aemkN0vkS//eTkNzdE/Pm3x5lquu+8tBmal8NhNZQzpr4ImIiLdU0mTmJaXnsyvr5rKzromvvTwcto7ur1yWFR4ckUlNzzwFsNyU3n0pjIGZqV4HUlERKKYSprEvFnFufzoYxN5dXMtt/59bdRNJHDOcefLFXz9sVVMH9afx75QpkkCIiJyXAleBxDpDVfOKGJb7WHufmUreenJfO280V5HAqClvYMfPL2WR5fuYt6Uwdx+xWkkJ/i9jiUiIjFAJU3ixrcuGMv+w63c8cJmUhL93HxWiad59tYf4Yt/Ws7KXQf50tklfOP8MbrUk4iI9JhKmsQNn8/42cdPo7mtg5/+YwMtbQFuOXckZpEvRq9vqeWWR1ZwpLWD339mGhdMHBTxDCIiEttU0iSu+EMzPpMT/PzyhU1UH2rmh5dOIMEfmeGXLe0d/M/zm7jnla0U56Xx6E2zGFmQEZF9i4hIfFFJk7iT6Pfxi0+cRl5GEncv2kpFdSO/uXoqBRnhnU25atdBvvW3t9lQdYirZhTxvYvHkZasv2IiInJy9BtE4pKZ8Z0LxzG6IIPvPrWai3+9mJ9cNonzxg/o9X3VHW7l589t4NGlu8hLT+a+60s5Z2zv70dERPoWlTSJax+fPoQJhZl89ZGVfO6P5Vw8aRDfvnAsQ3NOfRHZg02t3Ld4G/e/tp2mtg5uPGMEXz1vFBkpib2QXERE+jqLtjWlTlVpaakrLy/3OoZEmdb2AHcv2sJvXq7AOccnS4dywxnDT2q82Jrd9fy1fBd/W76bxpZ2Lpw4kK+fP5rRAzT2TEREToyZLXPOlXb7nEqa9CV7649w58sVPLZ0F20djmlF2Zw7bgBnjMxj7MAMUhLfv4ZZzaEW1u9t4LUttSzcUMPGfYdISvBx4cSBfHFuCWMHZnrwnYiISDxQSRPpouZQC0+uqOTJFXtYv/fdi7MPykohMyWRpAQfh1va2X+4lfojwWuCJvqN6cP6c9GkQcybXEhWqk5riojIqVFJEzmGmkMtvLWtjorqRnbUHeZwSzst7QHSkhPon5rIiLx0xg7MYMrQbM3WFBGRXnWskqbfONLn5Wckc/FpWmxWRESiiy6wLiIiIhKFVNJEREREopBKmoiIiEgUUkkTERERiUIqaSIiIiJRSCVNREREJAqppImIiIhEIZU0ERERkSikkiYiIiIShVTSRERERKKQSpqIiIhIFFJJExEREYlCKmkiIiIiUcicc15n6FVmVgPsiMCu8oDaCOxHek7vSXTS+xJ99J5EJ70v0ScS78kw51x+d0/EXUmLFDMrd86Vep1D3qX3JDrpfYk+ek+ik96X6OP1e6LTnSIiIiJRSCVNREREJAqppJ28e7wOIO+j9yQ66X2JPnpPopPel+jj6XuiMWkiIiIiUUhH0kRERESikEraCTKzC8xso5lVmNm3vc4jYGZDzexlM1tnZmvN7KteZ5IgM/Ob2Qoz+z+vs0iQmWWb2eNmtsHM1ptZmdeZ+joz+3ro3641ZvaImaV4nakvMrP7zKzazNZ0eizHzJ43s82hz/0jmUkl7QSYmR+4E7gQGA9cZWbjvU0lQDvwDefceGAW8CW9L1Hjq8B6r0PIe/wKeNY5NxaYjN4fT5lZIXALUOqcmwj4gSu9TdVnPQBc0OWxbwMvOudGAS+G7keMStqJmQFUOOe2OudagUeBeR5n6vOcc3udc8tDtw8R/KVT6G0qMbMhwMXAH7zOIkFmlgWcCdwL4Jxrdc4d9DaVAAlAPzNLAFKBPR7n6ZOcc68AdV0engc8GLr94P9v7w5CrKriOI5/f6Qtsl2FVAYjES6qhblJ24i2KyxCWhUSbQtaBbXJTVAhEgS1CSNQgjAhidKCIlpJWBODRZvKnGEmXRVG1MT8W9w7OUIMDAPv3N77fuDx7rub+9u8x++dc+49wMOjzGRJW5tbgQsrPs9iGRiUJFPAduBM2yQCXgWeBZZaB9G/tgKXgLf6aeg3k2xqHWqSVdUccAj4GZgHfq2qj9um0gqbq2q+P14ANo/y4pY0jY0k1wPvAc9U1W+t80yyJA8CF6vqbOssusoG4B7gjaraDvzOiKdvdLV+jdNDdAX6FmBTksfaptJ/qe5xGCN9JIYlbW3mgNtWfN7Sn1NjSTbSFbRjVXWidR5xH7AvyU90ywL2JDnaNpLoRv9nq2p5pPk4XWlTO/cDP1bVpapaBE4Auxpn0hW/JLkZoH+/OMqLW9LW5kvgjiRbk1xLt7jzZONMEy9J6NbYfFdVh1vnEVTVc1W1paqm6L4nn1aVowONVdUCcCHJtv7UXuDbhpHUTXPem+S6/rdsL97MMSQngQP98QHg/VFefMMoL/Z/V1V/J3kKOE13B86RqjrXOJa6UZvHgZkk0/2556vqw4aZpKF6GjjW/9H8AXiicZ6JVlVnkhwHvqK7U/1r3HmgiSTvALuBG5PMAi8ALwHvJnkSOA88OtJM7jggSZI0PE53SpIkDZAlTZIkaYAsaZIkSQNkSZMkSRogS5okSdIAWdIkTaQkNySZ7l8LSeb648tJXm+dT5J8BIekiZfkIHC5qg61ziJJyxxJk6QVkuxO8kF/fDDJ20m+SHI+ySNJXkkyk+RUvx0ZSXYkjcbulwAAAOBJREFU+TzJ2SSnl7eRkaT1sKRJ0upuB/YA+4CjwGdVdTfwB/BAX9ReA/ZX1Q7gCPBiq7CSxofbQknS6j6qqsUkM3TbwZ3qz88AU8A24C7gk27rRa4B5hvklDRmLGmStLo/AapqKcliXVnIu0T3GxrgXFXtbBVQ0nhyulOS1ud74KYkOwGSbExyZ+NMksaAJU2S1qGq/gL2Ay8n+QaYBna1TSVpHPgIDkmSpAFyJE2SJGmALGmSJEkDZEmTJEkaIEuaJEnSAFnSJEmSBsiSJkmSNECWNEmSpAGypEmSJA3QP8oJ707OmSmjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Lets plot the solution\n",
    "plt.figure(figsize = (10,6))\n",
    "plt.plot(t, x_underdamped)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Displacement')\n",
    "plt.title('Underdamped Oscillations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, one can get the velocity vs time plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAGDCAYAAABwRoerAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZycVZ32/8+3eqve9z37vkASSCeEsGtwYZAgDiCboCIu46jP6OPo6OOooyPqOKP+1BkRFXBDRxRQQfY1CyQkgYQkJOnOnt73fak6vz+qGpvQSTrdXXVXdV/v16tfXcvdVVd3An3l3Oec25xziIiIiEhs8XkdQERERETeTCVNREREJAappImIiIjEIJU0ERERkRikkiYiIiISg1TSRERERGKQSpqIjImZ3WJmz4/j691lZl8br9eLFDN72sxujfJ7XmBmrw25f8DM1oRvf9nMfjmG137VzC4eh5giMk5U0kQmGTNzZjbnuMfG9Ate3szMFpnZg2bWambtZvaUma0ey2s6555zzs0fh2xvKsLOucXOuafH+toiMn5U0kTEM2aW4HWGSDCz2cA6YDswEygD/gg8ambneplNROKHSpqIvIGZXWxmR8zs02ZWZ2bVZvb+Ic/nh0eI2szsRWD2cV+/wMweM7MmM3vNzK4Z8txdZvbfZvaQmXUCl5jZWWa2JTza9FvAP+T4XDP7s5nVm1lz+PaUIc8/bWZfM7P1ZtZhZn8K5/tVON8mM5sx5HhnZp8wsyozazCzb5uZb8jzHzCzXeH3esTMpg957lIz2x0eGfsBYCf5MX4Z2OCc+4Jzrsk51+6c+z7wC+Cb4dfzm9kvzazRzFrCWYvDz+WZ2c/N7Fg4y/1D/2xG+Of4v2ZWE877rJktDj9+G3AD8NnBn1n48aGnTlPM7Lvh9z8Wvp0yNMNJ/n5cZmY7w3+eR83sMyPJKyJvppImIsMpAbKBcuCDwA/NLDf83A+BHqAU+ED4AwAzSwceA34NFAHvBX5kZouGvPb1wNeBTOBF4H5C5SUP+F/gPUOO9QE/B6YD04Bu4AfHZX0vcFM462xgQ/hr8oBdwL8ed/y7gQrgbGDtYH4zWwv8C3AVUAg8B/wm/FwB8Afgi0ABUAmcd4KfHcCl4e/leL8DzjOzVOBmQj/jqUA+8JHw90f455EGLCb0c/yvk7zXiTwMzA1//RbgVwDOuTvCt7/lnMtwzr1rmK/9ArAKWAYsBVYS+t4Hnezvx0+BDzvnMoEzgCdHkV1EUEkTkeH1A191zvU75x4COoD54dOT7wG+5JzrdM7tAO4e8nWXAweccz93zg0457YC9wFXDznmAefcOudckFAJSAK+G36v3wObBg90zjU65+5zznU559oJlbuLjsv6c+dcpXOulVAxqXTOPe6cGyBUlM467vhvhke3DgHfBa4LP/4R4BvOuV3hr/13YFl4NO0y4FXn3O+dc/3hr6s5yc+vAKge5vFqQv/fzSP0M84H5jjnAs65l5xzbWZWCrwT+Ihzrjn8c3nmJO81LOfcz8IjeL2ERvaWmln2CL/8BkJ//nXOuXrgK4SK8KBh/34MeW6RmWWF82853ewiEqKSJjL5BAgVo6GSCP1yHdQYLiqDuoAMQiNMicDhIc8dHHJ7OnBO+PRdi5m1EPqFXzLkmKFfWwYcdc654V7PzNLM7MdmdtDM2oBngZzj5rLVDrndPcz9jOO+1+Ozlw3J/r0huZsIndIsDx/z+teF8w59neM1EBppPF4pEASaCY2WPQLcGz6l+C0zSyI0stbknGs+yeuflJklmNntZlYZ/rkdCD9VMMKXKOONf65Df05w4r8fECrxlwEHzewZ0xw8kVFTSROZfA4BM457bCZv/KV8IvXAAKEiMWjakNuHgWecczlDPjKccx8dcszQQlYNlJvZ0PldQ1/v04RGaM5xzmUBF4YfP9l8sFM5PvuxIdk/fFz2VOfc+nDO178unHfo6xzvcd44ejjoGkJz1brCo1Bfcc4tAlYTGoV8XzhHnpnljPYbJHRKeS2whtBpyRmD0cOf3TBfM9QxQqV10NCf00k55zY559YSOs16P6FTvCIyCippIpPPb4EvmtkUM/OFJ4u/C/j9qb7QORcgNDfry+FRrkWE5lYN+jMwz8xuMrOk8McKM1t4gpfcQKj0fSJ87FWE5j8NyiQ0GtZiZnm8eX7ZaPzf8IKEqcAnCf08AP4H+PyQCfbZZjZYtP4CLDazq8wsEfgEbxwdPN5XgNVm9vXwIoBMM/tHQiXsn8Ovf4mZnRkeFWwjNJIZdM5VEzpt+6NwziQzu/BEb3QCmUAv0Ehobtu/H/d8LTDrJF//G0J/RwrD8/G+BJxyixYzSzazG8wsO3xauI3QyKGIjIJKmsjk81VgPfA8odNu3wJuCM8vG4mPEzq1VQPcRWiSPgDheWNvIzSZ/1j4mG8CKcO9kHOuj9BE/VsInV68llAJHPRdIJXQ6cONwF9HmPFkHgBeArYRKl8/DWf5YzjrveFThDsIzQ3DOddAaGTsdkLFZy6hLTaG5ZzbC5xPaNL9AUIjce8B3u6cG/y6EkLFuI3QAodnCJ0ChdD8r35gN1AHfOo0v8d7CI2MHgV2EvrZDfVTQvPGWgZXjh7na8Bm4BVC24hsCT82EjcBB8I/w48QOt0tIqNgb5wKIiIycZmZA+Y65/Z5nUVE5FQ0kiYiIiISg1TSRERERGKQTneKiIiIxCCNpImIiIjEIJU0ERERkRiU6HWA8VZQUOBmzJjhdQwRERGRU3rppZcanHOFwz034UrajBkz2Lx5s9cxRERERE7JzE54tRed7hQRERGJQSppIiIiIjFIJU1EREQkBqmkiYiIiMQglTQRERGRGKSSJiIiIhKDVNJEREREYpBKmoiIiEgMUkkTERERiUEqaSIiIiIxSCVNREREJAZNuGt3ikw2fQNBjjR30dzVR2dvgJREH+kpiUzLTyPLn+R1PBERGSWVNJE4Ewg6XjrYzGM7a1i3r5G9de30B9ywxxZkpLB8eg6rZuXz9sUllOWkRjmtiIiMljk3/P/c41VFRYXbvHmz1zFExl17Tz/3vniYezYe4HBTN0kJxooZeZw5JZt5RZnkZySTkZJI70CQtu5+DjZ1saemnU0Hmzjc1A1AxfRcbjp3OpedWUpSgmY7iIh4zcxecs5VDPecRtJEYlx/IMivXzjE957YS1NnHytn5PGZt83nLQuKyBzh6cwDDZ38ZXs19710hE/eu41vPLSbD54/k5vOnY4/KSHC34GIiIyGRtJEYtjumjb+6bcvs7O6jVWz8vj8OxeydGrOqF8vGHQ8s6eenzxXxfrKRspzUvnM2+exdmk5Pp+NY3IRERmJk42kqaSJxCDnHD9fd4BvPLyL7NQkvnblmbx9cTFm41ek1u1r4BsP72LH0TZWzszj9qvOZFZhxri9voiInJpKmkgc6ekP8IU/7uC+LUdYs7CYb77nTPIzUiLyXsGg439fOszX/rKL3oEgn750Hh+6YJZG1UREokRz0kTiRHtPPx+8azMvHmjiU2vm8om3zI1oYfL5jGtXTOOS+UV88f4dfOPh3ayrbOQ/r1lKQYSKoYiIjIyWd4nEiObOPm648wW2HGrme+9dxqfWzIvaiFZRlp8f37Scr7/7DDZWNXLZ955jQ2VjVN5bRESGp5ImEgPae/q56WcvsLumnf+5cTlrl5VHPYOZccM507n/Y+eRkZLIjT99gV9uPBj1HCIiEqKSJuKxnv4At969md3V7fz4xuWsWVTsaZ5FZVk88PHzuHBuAV+8fwdf+dOrBIITa+6qiEg8UEkT8ZBzjk//7mVe2N/Ed65ZyiULiryOBECmP4k7b17BB86byc/XHeBD92ymuy/gdSwRkUlFJU3EQz96upK/bK/m8+9c4MkpzpNJ8BlfetcivnblGTz9Wh03/fQFWrv7vY4lIjJpqKSJeOTJ3bX8x6OvsXZZGbddOMvrOCd046rp/OD6s3n5SAvX3bGRho5eryOJiEwKKmkiHjjU2MUnf7ONRaVZ3H7VknHdpDYSLjuzlDtvXkFVQwfX/M8GjrZ0ex1JRGTCU0kTibKBQJBP/XYrGPzPjctJTY6Pa2deNK+QX37wHOo7ernujo1Ut6qoiYhEkqclzczeYWavmdk+M/vcMM//k5ntNLNXzOwJM5vuRU6R8fSjpyvZcqiFr115BlPz0ryOc1oqZuTxiw+eQ3NnH9fdsZHath6vI4mITFielTQzSwB+CLwTWARcZ2aLjjtsK1DhnFsC/B74VnRTioyvbYdb+N4Te1m7rCzmFgqM1LKpOdz1gZXUt/dy3U82UteuoiYiEglejqStBPY556qcc33AvcDaoQc4555yznWF724EpkQ5o8i46R0I8OnfbaM4M4Wvrj3D6zhjsnx6Lnd9YCU1rT3c8JMXaOrs8zqSiMiE42VJKwcOD7l/JPzYiXwQeDiiiUQi6MfPVFFZ38nXrzqT7NQkr+OM2YoZefzslhUcauri/XdtorN3wOtIIiITSlwsHDCzG4EK4NsneP42M9tsZpvr6+ujG05kBKrqO/jBU/u4fEkpl8yPjQ1rx8OqWfn84Pqz2XG0lY/88iX6BoJeRxIRmTC8LGlHgalD7k8JP/YGZrYG+AJwhXNu2A2anHN3OOcqnHMVhYWFEQkrMlrOOb54/w5SEn186fLjp13Gv0sXFfONq87kub0N/NPvthHUJaRERMZFoofvvQmYa2YzCZWz9wLXDz3AzM4Cfgy8wzlXF/2IImP3p1eqWV/ZyNeuPIOiLL/XcSLimoqpNHX2cfvDu8lLT+YrVyyO+b3fRERinWclzTk3YGYfBx4BEoCfOedeNbOvApudcw8SOr2ZAfxv+H/4h5xzV3iVWeR09fQH+ObDu1lclsX1K6d5HSeiPnLRbJo6+7jj2SrKclL5yEWzvY4kIhLXvBxJwzn3EPDQcY99acjtNVEPJTKOfrZuP0dbuvn21Uvw+Sb+yNLn3rGAYy3d3P7wbqbmpvF3S0q9jiQiErc8LWkiE1l9ey8/eqqSNQuLWT27wOs4UeHzGf9x9VKqW3v4P7/bRkm2n+XTc72OJSISl+JidadIPPru43vo6Q/w+csWeB0lqvxJCfzkfRWUZvv50D2bOdjY6XUkEZG4pJImEgEHGzu5d9Nhrj9nGrMLM7yOE3V56cn8/JYVBJ3j/XdtoqVLm92KiJwulTSRCPj+E/tI9Bkfv2SO11E8M6swgztuquBIUze3/eIlegcCXkcSEYkrKmki46yyvoM/bj3CTaumT9gtN0Zq5cw8vn31El7c38QX/7gD57SHmojISGnhgMg4+97je0lJTOAjF2sLCoC1y8qprOvg+0/uY35JJrdeMMvrSCIicUEjaSLjaG9tO3965Rg3r55BQUaK13FixqfWzOPti4v594d28dRr2pdaRGQkVNJExtF/P11JalICt12o0aKhfD7jP69ZxvySLD7x663sq+vwOpKISMxTSRMZJ0eau3jg5WNct3IaeenJXseJOekpifzkfctJSfJx691a8SkicioqaSLj5M7n9uMzuPWCmV5HiVlTctP48U3LOdbSwz/8egv9gaDXkUREYpZKmsg4aOzo5d5Nh7hyWTml2alex4lpy6fn8fV3n8G6fY382593eh1HRCRmaXWnyDi4e/0BegeCfPgizUUbiasrprK3roM7nq1ibnEmN62a7nUkEZGYo5E0kTHq7gtwz8aDXLqwmDlFmV7HiRv//I4FXDK/kC8/+Crr9zV4HUdEJOaopImM0QPbjtLS1a/9v05Tgs/4/nVnMbMgnY/9egsHGnSNTxGRoVTSRMbAOcdd6w+wsDSLFTNyvY4TdzL9Sfz05goAbr1nM209/R4nEhGJHSppImOwsaqJ3TXtvH/1DMzM6zhxaXp+Oj+64WwONHTyid9sJRDUpaNEREAlTWRM7lq/n9y0JK5YVuZ1lLi2enYBX1m7mKdfq+f2h3d5HUdEJCZodafIKB1p7uKxnbV85KLZ+JMSvI4T9244Zzp7atr5yXP7mVucyTUVU72OJCLiKY2kiYzSLzYexMy4UdtHjJv/d/kizp9TwBf+uJ1NB5q8jiMi4imVNJFR6OkP8NtNh3nbomLKcrR57XhJTPDxw+vPZkpuGh/5xUscburyOpKIiGdU0kRG4a87amjp6tcoWgRkpyVx580V9AWCfOiezXT2DngdSUTEEyppIqNw76ZDTMtL49xZ+V5HmZBmF2bww+vPZk9tO5/67TaCWvEpIpOQSprIadrf0MnGqiauXTEVn0/bbkTKhfMK+X+XL+KxnbX8x6OveR1HRCTqtLpT5DTdu+kQCT7j6uVTvI4y4d2yegZ7ajv40dOVzCvO5Mqzyr2OJCISNRpJEzkNfQNB7nvpCG9dUERRlt/rOBOemfGVKxZzzsw8PnvfK2w91Ox1JBGRqFFJEzkNT+yqpaGjj+tWTvM6yqSRnOjjv29cTnFWCrf94iWOtXR7HUlEJCpU0kROw282HaY028+F8wq9jjKp5KUn89ObV9DVO6AVnyIyaaikiYxQdWs3z+2t5+rlU0jQgoGom1ecyQ+uP5td1W3842+2MhAIeh1JRCSiVNJERuj+rcdwDt6jBQOeuWRBEV9dewZP7q7jXx98Fee0NYeITFxa3SkyAs45/rDlCBXTc5men+51nEntxlXTOdLczf88U8mU3DQ+evFsryOJiESERtJERuDVY23srevg3WdrC4hY8Nm3z+ddS8v45l938+DLx7yOIyISERpJExmB+7YcITnBx+VnlnkdRQCfz/iPq5dQ29bDZ373MsWZKZyjqz+IyASjkTSRU+gPBHlw2zHWLCoiOy3J6zgSlpKYwB03LWdqXiofumcz++ravY4kIjKuVNJETuG5vfU0dvbx7rO0YCDW5KQlc9f7V5KcmMAtP99EbVuP15FERMaNSprIKdy35Sh56clcpL3RYtLUvDR+fssKmjv7uOmnL9DS1ed1JBGRcaGSJnIS7T39PLazlnctKSU5Uf+5xKozp2Tzk5srONDYxS0/36TNbkVkQtBvHZGTeGxnLX0DQa5YplWdsW717AL+v+vO4pUjLXz4Fy/ROxDwOpKIyJiopImcxJ9fqaY8J5Wzp+V4HUVG4O2LS/jme5bw/L4GPvmbbboqgYjENW3BIXICrV39PLe3nvefNxMzXQYqXlxdMZW2ngH+7c87+Zc/bueb71miPz/xhHOO9t4BWrv6ae7qo6N3gN7+ID39AXoGAvSEbw8EHA4X/hoYeh2NpAQfKYmhj+REHymJCaQkhe6nJyeS6U8k059Epj8Rf1KCN9+oRIxKmsgJPPJqDf0Bx+VLSr2OIqfpg+fPpLWrj+8/uY/kRB//tvYMFTUZN8456tp7OdrSTU1rD9WtPdS0dlPd2kNtWw+NnX20dPXT2t1PIBi9S5clJ/jCpS2RrNQk8tKTyU9PIT8jmfz0ZPLSkynISAk9nhG6rWIX21TSRE7gT68cY1peGmeWZ3sdRUbh/1w6j96BID9+tgpARU1OW3dfgD217VQ1dLC/vpOqhk72hz+6+t445zEl0UdZTirFWSksKs0iJy2JnNRkctKSyE4NfWSER7tSkxLwJyXgT/LhT0wgMcEwMwb/dpqBYTgc/QOO3kCA3v4gvQNBegcC9A0E6ekP0tk7QHtvP+09A7T3DNDW009H+HZrdz+NHX3sre2goaOX3oHhT/1n+RMpzvJTku2nOMtPcVYKJVmDt0OPF2SkkODTfzteUEkTGUZjRy/rKxv58IWz9Is9TpkZn3vnAgB+/GwVhvHVtYv15ynDaurs49Vjrew81sbO6jZePdZGVX0HgwNhPgtt9zKzIJ2VM/OYWZDOlNxUSrJSKcvxk52aFJm/W8kAY9tE2zlHV1+Axo4+Gjt7aezoo6mzj/qOXuraeqhp66G2rZd9dQ3Utfe+afTPZ1CY+bfyVprtpzg7/DnLT2l2KiVZflKTNSo33jwtaWb2DuB7QAJwp3Pu9uOeTwHuAZYDjcC1zrkD0c4pk89fX60hEHRcvkSXgYpng0XNAXc8W4UZfOUKFbXJLhh07KvvYPOBZjYfbGLLwWYONHa9/nxZtp9FZVlcdmYpi0ozmVOUybS8tLjdhsfMSE9JJD0lkWn5aSc9NhB0NHb2UtvaS224wA0WuZq2Xg42drGxqpG2njdvc5OdmjSkuL250JVkRbDMTlCelTQzSwB+CFwKHAE2mdmDzrmdQw77INDsnJtjZu8FvglcG/20Mtn8+eVqZhWms7A00+soMkZmxufDI2p3PFtF0Dm+esUZ+HT6ZtIIBh27atpYt6+B9ZWNbDnY/HrJyE9PZvn0XN67chpnlmezqDSL3PRkjxN7J8FnFGX6Kcr0cyYnnurR1TdATWu4vIXn5dW2Dc7P62FndRsNHb2446bk+ZN8lGaHTguHPr+x0On06ht5OZK2EtjnnKsCMLN7gbXA0JK2Fvhy+PbvgR+YmTl3/B97dPX0B+gPBMn06zqOE1Fdew8v7G/k45fM0b/4JojBomYGP36mirbuAb5zzVKSEuJzZERO7UhzF8/vbeD5cDFr6gxdiWJuUQZ/t6SMium5LJ+ey/T8NP13PgppyYnMKsxgVmHGCY/pDwSpa+99fVFFzeBHuNhtOtBEbVsP/YE3/koPFcUUSsKjb4OfCzIGF0GEPuelJ0/4hQ9elrRy4PCQ+0eAc050jHNuwMxagXygISoJh+GcY9lXH+V9587gXy5b6FUMiaCHt9cQdHD5Up3qnEhCRW0hOanJfPOvu2nv6edHNyzXPJoJYiAQZPPBZp7YVcsTu+qoaugEoDgrhYvnF3L+nALOm1NAcZbf46STR1KCj/KcVMpzUk94TDDoaOzs+9soXNsbV8q+VtvOM3vq37RQY1B6cgL54RWrBeHilpeeQm5aElmpSWT5k8hKTSTLH1q8kZUa2q4kXv6BNiEWDpjZbcBtANOmTYv0e1Gc5aemVRdynqge3lHN3KIM5hXrVOdE9NGLZ5OdmsQX7t/O+372AnfevILsVI2Kx6O2nn6e3VPP4ztreeq1elq7+0lO8HHu7HxuXDWdC+YWMKcoQyNlMcznMwozUyjMTOGME6ykH9xvLrTgoZeG8MKHps4+Gjp6X799tKWH7Udbaerse9Po3PHSkhPICu8vl56SSFpyAmnJiaSnJLx+Oy05gZJsPzecMz0S3/qIeFnSjgJTh9yfEn5suGOOmFkikE1oAcEbOOfuAO4AqKioiPip0OIsPzVtKmkTUVNnHy/ub+IfLpnjdRSJoOvPmUZ2ahKf+u1Wrv3xBn56y4qT/mtfYkd3X4BHd9bwwLZjPLe3nv6AIy89mTULi7l0URHnzy0kI2VCjD9ImJmFRsT8ScwsSD/l8YOrWdt6+mnrDm1N0trVH77fT1vPAG3doX3s2nr66eoL0NUXoLmrm66+gdD93gG6+gPMK8qctCVtEzDXzGYSKmPvBa4/7pgHgZuBDcDfA096PR8NQiVt+5EWr2NIBDy+q5agC11eSCa2v1tSSnZqEh/91Uus/cE67ry5gmVTdfmvWDQQCLK+spH7tx3lkR01dPYFKMv284HzZnLpomLOmpariebyuqGrWUvHsM2lc+6E+8tFi2clLTzH7OPAI4S24PiZc+5VM/sqsNk59yDwU+AXZrYPaCJU5DxXkpXCY209OOc0jD7BPLKjhvKcVBaXZXkdRaLg/LkF/OGjq/nA3Zu49scb+K9rl3HZmbrCRCxwzrHjaBv3bzvKgy8fo769l0x/Iu9aWsaVZ5WzckaeVuhKRJmZ5wsTPB0Tds49BDx03GNfGnK7B7g62rlOpTjLT09/kLaeAc1lmUA6egd4bl8DN5wzTeV7EplbnMn9HzuPD//iJT72qy18+tJ5/MMlc1QAPHK4qYsHth3lj1uPUlnfSXKCj0sWFPLus8q5eH6R5780RaJJJ+5HYXB1UG1bj0raBPLMa/X0DQR1qnMSys9I4Ze3nsPn/7Cd7zy2h5ePtPCdq5eRnab/vqOhubOPP2+v5oGtR9l8sBmAlTPzuPWCWVx2Rqn+HGTSUkkbhaElTSsAJ45HXq0hLz2ZFTPyvI4iHvAnJfCf1yxl6ZRsvv7QLi7/wXP89w3LT7jiTMampz/A47tquX/rMZ7ZU0d/wDGvOIPPvmM+VywtY0ruyXfGF5kMVNJGoSRc0rQNx8TRNxDkqd11vPPMEk1AnsTMjFvOm8mSqTn8w6+2cNV/r+cLly3kfedO1ynwcRAIOjZWNXL/1qM8vKOGjt4BirNSeP95M7lyWTkLSzP1cxYZQiVtFIqyUgCoa+/1OImMl/WVDbT3DuhUpwBw9rRc/vyP5/Pp/32Zf33wVZ7YXce3/36JNkIdBeccO6vbeGDbMR7YdpTatl4yUhJ55xklvPuscs6Zla9/GImcgEraKPiTEshJS9JI2gTyyKu1pCcncN6cAq+jSIzIz0jh57es4JcvHOLrf9nJ27/7LF+78gz+7sxSjfaMwNGWbh7YdpT7tx5lT20HiT7j4vlF/L/Ly1izsFgLAERGQCVtlIoz/dRqQ9sJIRB0PLazVivH5E3MjJtWTWf17Hz+6bfb+Pivt/KHBUf5yhWLmZqnOVPHa+3q56Ed1fxx61Fe3N8EQMX0XP4tXG7zJvGFy0VGQyVtlIqzVdImim2HW2jo6OVti4u9jiIxanZhBvd9dDV3rT/Afz62h0v/6xk++dZ53HrBzLi5BmCk9PQHeHJ3HQ9sO8pTu+vpCwSZVZjOpy+dx9pl5UzLV5kVGS2VtFEqzkxhT0271zFkHDy5u5YEn3HxvCKvo0gMS0zwhbaEOLOULz/4Kt/8625+t/kw//ft83nnGSWT6hRofyDI8/sa+NO2Yzy6s5aO3gEKMlK4cdV03n1WOWeUZ02qn4dIpKikjVJJtp/6jl4CQadJr3Huyd31VEzP1V5MMiJlOanc8b4Kntpdxzce3sXHfrWFZVNz+Nw7F7BqVr7X8SImGHRsOtDEgy8f4+EdNTR19pHpT+SyM0tYu6ycVVoAIDLuVNJGqSjLTyDoaOzopUgrvuLWsZZudlW38S+XLfA6isSZSxYUceG8Qu576UU9FGsAACAASURBVAj/+dge3nvHRlbMyOWjF8/mkvlFE2IkqT8Q5MX9TTz6ag2PvFpLTVsP/iQfaxYWc8XSMi6aX0hKouZxikSKStoovb5XWluPSloce3J3HQBvWaD5aHL6EnzGNSum8q6lZdy76RA/ebaKD9y1mQUlmXzgvJlcvrSUtOT4+t9sV98Az7xWz6M7a3liVy1tPQP4k3xcMLeQzy9ZwJqFxaSnxNf3JBKv9F/aKBWH90qrae1hyRSPw8ioPbm7jun5acwuTPc6isSx1OQE3n/eTG5cNZ0Htx3jx89W8tn7XuGrf97JFcvKuLZiKkumZMfk6Jpzjl3V7Ty3t57n9jbw4oEm+gaC5KQlcemiEt62uJgL5xaSmqwRM5FoU0kbpcGRtFptaBu3uvsCrNvXwPW6oLqMk6QEH+9ZPoWrzi7npYPN/ObFw/xhyxF+/cIhpuSm8o7FJbz9jBLOnpbr2fytYNCxt66Dlw42s+lAE8/tbaChI/T/sfnFmdy0ajprFhazYkYuiZN85aqI11TSRik/I4UEn1GrDW3j1vrKBnoHgrxVpzplnJkZFTPyqJiRx79esYiHt1fz1x013LPhIHc+v5/s1CRWzMhj1aw8zpmZz7ySjIjM7QoGHYeauthd08bO6nZePtzClkPNtPcMAFCQkcx5cwq4YG4h588poCRbUzdEYolK2igl+IzCjBTtlRbHnthdR3pyAitn6oLqEjlZ/iSuXTGNa1dMo72nn6deq2fd3gY27m/k8V21ACT6jJkF6SwozWJWQTplOX5KslMpzfaT5U8iLSWB9OTE10ffgkHHQNDR0TtAS1cfzV39tHT1cay1hyNNXRxp7uZwcxf76jro6gsA4DOYU5TB5UvKqJiey/LpuUzPT9MoskgMU0kbg+KsFGpU0uKSc44nd9Vx4bxCkhN1SkeiI9OfxBVLy7hiaRkQmtO66UATu2va2F3dzpaDzfzp5WMn/PpEnxFwDudO/B7JCT7Kc1OZkpvKtSumsrAkiwWlmcwtytS8MpE4o5I2BsVZfg42dnkdQ0ZhZ3UbNW09vGWBNrAV75Rk+3nX0jLeFS5tAH0DQWrbeqhp66G6tYeOngG6+gbo7A3QOxAg0Wf4fEaCGRn+RHLTkslOSyInNYmSbD/FmX582q9MZEJQSRuD4iw/L4SvTyfx5clddZjBxfNV0iS2JCf6mJqXpmuDigg6zzMGJdl+Wrv76ekPeB1FTtMTu+tYOiWHwswUr6OIiIgMSyVtDIrCv+C1eCC+1Lf38vKRFt6qU50iIhLDVNLGYHC5em2b9kqLJ8/uqce50GV9REREYpVK2hgUD7k0lMSPZ/fWU5CRzKLSLK+jiIiInJBK2hgMlrQ6lbS4EQw6ntvbwAVzC7UCTkREYppK2hhk+RPxJ/mo0VUH4saOY600dfZx0bxCr6OIiIiclEraGJgZJVl+Xb8zjjzzWj0A588t8DiJiIjIyamkjVFRll/X74wjz+6t58zybAoytPWGiIjENpW0MSrJ8mvhQJxo6+lny6EWLpynUTQREYl9KmljVJwVusi6O9nF9CQmrN/XQCDouGiett4QEZHYp5I2RsVZfnoHgrR293sdRU7hmT0NZKQkcta0HK+jiIiInJJK2hgNbmirU56xzTnHs3vqWT07n6QE/bUXEZHYp99WY1SanQpAdYtKWiyrrO/kaEs3F83X1hsiIhIfVNLGqDQ8klatFZ4x7Zk9oa03LpyrkiYiIvFBJW2MijJT8BlUt3Z7HUVO4tk99cwqTGdqXprXUUREREZEJW2MEhN8FGX6OabTnTGrpz/AxqpGjaKJiEhcUUkbB6U5fmraNJIWq17c30TvQFDz0UREJK6opI2DsuxULRyIYc/sqSc50ceqmfleRxERERkxlbRxUJLtp7pVG9rGquf21rNyRh6pyQleRxERERkxlbRxUJrtp7s/oA1tY1BdWw97ajt0QXUREYk7Kmnj4PW90rQNR8xZV9kAwPlzVNJERCS+qKSNg9Kcwb3StHgg1qzb10hOWhKLSrO8jiIiInJaVNLGweCGttqGI7Y451i3r4HVs/Px+czrOCIiIqdFJW0cFGX6SfAZNTrdGVOqGjqpbu3hPJ3qFBGROORJSTOzPDN7zMz2hj/nDnPMMjPbYGavmtkrZnatF1lHIsFnFGWmcEynO2PKun2ajyYiIvHLq5G0zwFPOOfmAk+E7x+vC3ifc24x8A7gu2aWE8WMp6U026+RtBizbl8DU3JTmaZLQYmISBzyqqStBe4O374buPL4A5xze5xze8O3jwF1QMxuGV+ak6rVnTEkEHSsr2zk/DkFmGk+moiIxB+vSlqxc646fLsGKD7ZwWa2EkgGKk/w/G1mttnMNtfX149v0hEqzfJzrKVbG9rGiO1HW2nvGdB8NBERiVuJkXphM3scKBnmqS8MveOcc2Z2wmZjZqXAL4CbnXPB4Y5xzt0B3AFQUVHhSUsqzUmldyBIS1c/uenJXkSQIQbno62erUtBiYhIfIpYSXPOrTnRc2ZWa2alzrnqcAmrO8FxWcBfgC845zZGKOq4eH0bjtZulbQYsG5fAwtLs8jPSPE6ioiIyKh4dbrzQeDm8O2bgQeOP8DMkoE/Avc4534fxWyjMljStHjAe919ATYfaOb8ORpFExGR+OVVSbsduNTM9gJrwvcxswozuzN8zDXAhcAtZrYt/LHMm7inNnhpqGMqaZ7bfLCJvkBQ89FERCSuRex058k45xqBtw7z+Gbg1vDtXwK/jHK0USvMTAlvaKu90rz2/L4GkhKMlTPzvI4iIiIyarriwDhJ8BnFmSlU69JQnlu/r5Gzp+WSluzJv0FERETGhUraOCrNSdVVBzzW3NnHjmOtOtUpIiJxTyVtHOmqA97bUNWIc6ikiYhI3FNJG0el2X6qW3u0oa2Hnt/XQEZKIkunZHsdRUREZExU0sZRaXZoQ9vmrn6vo0xa6/c1sGpWHokJ+qstIiLxTb/JxlFZTnhD2xbNS/NCdWs3Bxq7OHe2TnWKiEj8O2VJM7OXzOwfzCw3GoHiWVlOaK+0oyppnthQ2QjAubO0ia2IiMS/kYykXQuUAZvM7F4ze7uZWYRzxaXywZLWrJLmhQ2VjeSmJbGgJNPrKCIiImN2ypLmnNvnnPsCMA/4NfAz4KCZfcXMtFvoEHnpyfiTfBpJ88iGqkbOmZmPz6d/Q4iISPwb0Zw0M1sCfAf4NnAfcDXQBjwZuWjxx8woz0nVSJoHDjd1caS5m1Wz9O8GERGZGE65JbuZvQS0AD8FPuec6w0/9YKZnRfJcPGoPDdNI2ke2FgVno+mRQMiIjJBjOS6OVc756qGPmBmM51z+51zV0UoV9wqz0llx9FWr2NMOhuqGslPT2ZecYbXUURERMbFSE53/n6EjwlQnuOnqbOPrr4Br6NMGs45NlY2smpWPlrTIiIiE8UJR9LMbAGwGMg2s6EjZlmAP9LB4lV5bmiF57GWHuYUaVQnGg41dXGstYePztbWGyIiMnGc7HTnfOByIAd415DH24EPRTJUPCvPSQNCe6WppEXH3/ZH06IBERGZOE5Y0pxzDwAPmNm5zrkNUcwU1wZH0rTCM3o2VDVSmJnC7EKVYhERmThOdrrzs865bwHXm9l1xz/vnPtERJPFqeLMFBJ8xtGWLq+jTArOOTZWaT6aiIhMPCc73bkr/HlzNIJMFIkJPkqy/BpJi5L9DZ3UtvXqUlAiIjLhnOx055/Cn++OXpyJoTw3VXulRcmG1/dHU0kTEZGJZSQXWH/MzHKG3M81s0ciGyu+TdFVB6JmQ2UjxVkpzMhP8zqKiIjIuBrJPmmFzrmWwTvOuWagKHKR4l95bio1bT0MBIJeR5nQQvPRmjhX89FERGQCGklJC5jZtME7ZjYdcJGLFP/KclIJOqhp6/E6yoRWWd9BQ0evTnWKiMiENJLLQn0BeN7MngEMuAC4LaKp4lx5zt+24ZiSq9NwkfK3/dF0vU4REZl4TlnSnHN/NbOzgVXhhz7lnGuIbKz49vpeaVo8EFEbqhopz0llal6q11FERETG3UhG0gBWAxcOuf/nCGSZMIaOpElkBIOh+WgXzy/UfDQREZmQRrK683bgk8DO8McnzezfIx0snvmTEijISNZIWgTtqWunqbNP+6OJiMiENZKRtMuAZc65IICZ3Q1sBf4lksHiXXmO9kqLpNfno2nRgIiITFAjWd0JoYusD8qORJCJRhvaRtbGqkam5qVqYYaIiExYIxlJ+waw1cyeIrS680LgcxFNNQGU56Ty5O46nHOaMzXOgkHHC/ubeNuiYq+jiIiIRMxIVnf+xsyeBlaEH/pn51xNRFNNAOU5qfT0B2ns7KMgI8XrOBPKrpo2Wrr6dapTREQmtBOWtPC2G0MdCX8uM7My59yWyMWKf+Xh03BHmrtV0sbZ4Hy0VVo0ICIiE9jJRtK+c5LnHPCWcc4yoQzu3XWkuYtlU3NOcbScjo1VjczIT6M0W/ujiYjIxHXCkuacuySaQSaaqeGRtENNXR4nmVgC4floly8p9TqKiIhIRI1kn7Q0M/uimd0Rvj/XzC6PfLT4lp6SSH56ModV0sbVzmNttPcM6FSniIhMeCPZguPnQB+hqw4AHAW+FrFEE8jUvDQON2kbjvG0oSp0RTJtYisiIhPdSErabOfct4B+AOdcF6GtOOQUpual6XTnONtQ2ciswnSKsvxeRxEREYmokZS0PjNLJbRYADObDfRGNNUEMS0vlWMt3QwEgl5HmRAGAkE2HWjWKJqIiEwKJyxpZvZDMzsf+DLwV2Cqmf0KeAL4bHTixbepuWkMBB3VrT1eR5kQth9tpaN3QPujiYjIpHCyLTj2AN8GSoHHgMeBLcAnnXMNUcgW96blhVZ4Hm7qYmqeLl80VhurmgDtjyYiIpPDCUfSnHPfc86dC1wE7AOuIrR32sfMbF6U8sW1wWJ2uFnz0sbDxqpG5hVnaHNgERGZFE45J805d9A5903n3FnAdcC7gV0RTzYBlGb7SfCZFg+Mg/5AkE0HmjSKJiIik8ZI9klLNLN3heejPQy8RmhUTU4hMcFHeU6qtuEYB9uPttLVF1BJExGRSeNkCwcuNbOfEbpm54eAvxDajuO9zrkHxvKmZpZnZo+Z2d7w59yTHJtlZkfM7AdjeU+vTM1L1UjaONhYFbpe5zkz8zxOIiIiEh0nG0n7PLAeWOicu8I592vnXOc4ve/ngCecc3MJrRb93EmO/Tfg2XF636iblpfGEc1JG7MNlY3ML84kX/PRRERkkjjZwoG3OOfudM41R+B91wJ3h2/fDVw53EFmthwoBh6NQIaomJqXRkNHH529A15HiVv9gSCbDzSzapZG0UREZPIYyWa2kVDsnKsO364hVMTewMx8hFaTfuZUL2Zmt5nZZjPbXF9fP75Jx2jwQuta4Tl6rxxppbtf89FERGRyiVhJM7PHzWzHMB9rhx7nnHOEr2ZwnI8BDznnjpzqvZxzdzjnKpxzFYWFheP0HYyPv+2VpsUDo/X6fDSVNBERmUROtpntmDjn1pzoOTOrNbNS51y1mZUCdcMcdi5wgZl9DMgAks2swzl3svlrMWdwrzQtHhi9jVWNLCjJJC892esoIiIiUePV6c4HgZvDt28G3rRa1Dl3g3NumnNuBqFTnvfEW0EDyE1LIiMlkcMqaaPSNzA4H02jaCIiMrl4VdJuBy41s73AmvB9zKzCzO70KFNEmBlT89JU0kZp+9GW8Hw0LRoQEZHJJWKnO0/GOdcIvHWYxzcDtw7z+F3AXREPFiFTc1PZ3zBeu5dMLoPX61w5UyNpIiIyuXg1kjapTMtL43BzF6E1EnI6NlRqPpqIiExOKmlRMDUvjZ7+IPUdvV5HiSt9A0E2H9T1OkVEZHJSSYuCv23DoXlpp+OVIy309AdV0kREZFJSSYuCafmhknagQSXtdGyobMQMLRoQEZFJSSUtCqbmpuEzONioxQOnY+P+RhaUZJGTpvloIiIy+aikRUFyoo8puWnsb9RI2kj1DgR46aCu1ykiIpOXSlqUTM9P44C24RixV460aj6aiIhMaippUTKzIJ0DjZ3ahmOEBuejnTNTI2kiIjI5qaRFyYz8dNp7Bmjq7PM6SlzYWNXIQs1HExGRSUwlLUpmFIRXeGrxwCn9bT6aTnWKiMjkpZIWJTPy0wHYr204Tunlw630DgS1aEBERCY1lbQomZKbRoLPtA3HCPxtPppG0kREZPJSSYuS5EQf5Tm60PpIbKxqZFFpFtlpSV5HERER8YxKWhTNCK/wlBPr6Q+w5ZDmo4mIiKikRdHM/DQONnRpG46T2Ha4hd6BIOeqpImIyCSnkhZF0/PTae8doFHbcJzQxqrQfLQV2h9NREQmOZW0KJpZEFrhqSsPnNjGqkYWl2WRnar5aCIiMrmppEXRjILBbThU0oYTmo/Wwiqt6hQREVFJi6YpuanhbTi0V9pwth5qoW8gyLmzVdJERERU0qIoKcHHlNxU9muF57A2VjXiM6iYofloIiIiKmlRNiM/nf31KmnDCc1Hy9Z8NBEREVTSom52YQb7GzoJBrUNx1A9/QG2Hm7RpaBERETCVNKibHZROt39AarberyOElO2HGzWfDQREZEhVNKibHZhBgCVdR0eJ4kt6yobSPAZK7WyU0REBFBJi7rXS1q9StpQ6ysbWTolm4yURK+jiIiIxASVtCgryEgmy5+okjZEe08/rxxpZfXsAq+jiIiIxAyVtCgzM2YXZVBZpxWeg17c30Qg6Fg9R6c6RUREBqmkeWB2YYZG0oZYt6+RlEQfZ0/L9TqKiIhIzFBJ88Ccogzq2ntp6+n3OkpMWF/ZQMWMXPxJCV5HERERiRkqaR4YXDxQpU1taejoZXdNu+ajiYiIHEclzQOzC0MXWtc2HKGrDACs1v5oIiIib6CS5oGpeWkkJZjmpRGaj5aRksiZ5dleRxEREYkpKmkeSErwMT0/XSUN2FDZwDkz80hM0F9FERGRofSb0SOzC9OpnORz0o62dHOgsYvVczQfTURE5HgqaR6ZXZjBwcZO+gNBr6N4Zv2+BkDz0URERIajkuaR2YUZ9Acch5u6vI7imQ2VjeSnJzO/ONPrKCIiIjFHJc0js4tC23Dsm6QrPJ1zrKtsYNXsfHw+8zqOiIhIzFFJ88iccEnbO0lLWlVDJ7VtvZyn/dFERESGpZLmkYyURKbkpvJaTbvXUTyh+WgiIiInp5LmofnFmeypnaQlrbKRsmw/0/PTvI4iIiISk1TSPDSvJJPK+o5Jt8IzGHRsqGpk9ZwCzDQfTUREZDgqaR5aUJJJf8Cxv2Fy7Ze2s7qNlq5+neoUERE5CU9KmpnlmdljZrY3/Dn3BMdNM7NHzWyXme00sxnRTRpZ88JbT0y2eWnP7Q3NRztfm9iKiIickFcjaZ8DnnDOzQWeCN8fzj3At51zC4GVQF2U8kXFrMJ0Enw26ealPbe3ngUlmRRl+b2OIiIiErO8KmlrgbvDt+8Grjz+ADNbBCQ65x4DcM51OOcm1M6vKYkJzCxIn1Qjad19ATYfaNYomoiIyCl4VdKKnXPV4ds1QPEwx8wDWszsD2a21cy+bWYJw72Ymd1mZpvNbHN9fX2kMkfE/OJMXptEI2kv7G+kLxDkgnmFXkcRERGJaREraWb2uJntGOZj7dDjnHMOcMO8RCJwAfAZYAUwC7hluPdyzt3hnKtwzlUUFsbXL/95xZkcauqiq2/A6yhR8dzeBpITfaycked1FBERkZiWGKkXds6tOdFzZlZrZqXOuWozK2X4uWZHgG3Ouarw19wPrAJ+GpHAHplfkoFzoctDLZmS43WciHtubz0rZ+SRmjzsoKiIiIiEeXW680Hg5vDtm4EHhjlmE5BjZoNDY28BdkYhW1RNphWetW097Knt4Py5mo8mIiJyKl6VtNuBS81sL7AmfB8zqzCzOwGccwFCpzqfMLPtgAE/8ShvxEzPTycl0TcpVngObr1xgUqaiIjIKUXsdOfJOOcagbcO8/hm4NYh9x8DlkQxWtQl+Iy5xRnsngQjac/tracgI5mFJVleRxEREYl5uuJADJhXnDnhS1ow6Fi3r4Hz5hTg8+lSUCIiIqeikhYDFpVmUd/eS317r9dRImZXTRsNHX1cMDe+Vt+KiIh4RSUtBiwuywbg1WOtHieJHM1HExEROT0qaTFgUVlojtarx9o8ThI5z+9tYH5xJsW6FJSIiMiIqKTFgOzUJKblpbFzgpa07r4ALx5o0tYbIiIip0ElLUYsLstixwQ93fnC/kb6BoIqaSIiIqdBJS1GLC7L4mBjF209/V5HGXdPv1ZPSqKPc2flex1FREQkbqikxYjBxQO7JuApz6dfq2P17Hz8SboUlIiIyEippMWIxRN08cD+hk4ONHZxyYIir6OIiIjEFZW0GFGU5acgI2XClbSndtcBcPE8lTQREZHToZIWQxaXZU24vdKeeq2O2YXpTMtP8zqKiIhIXFFJiyFnlGexr66D3oGA11HGRVffAC9UNXHJfI2iiYiInC6VtBiyuCybgaBjT02H11HGxfp9jfQFgpqPJiIiMgoqaTFkcPHA9qMT45TnU6/VkZ6cQMWMXK+jiIiIxB2VtBgyLS+NnLQkXjnS4nWUMXPO8fRr9Zw3p4CURG29ISIicrpU0mKImbF0Sg7bDsd/Sdtb18HRlm6d6hQRERkllbQYs3RqDntq2+nsHfA6ypi8vvXG/EKPk4iIiMQnlbQYc9bUHIIu/uelPfVaHQtKMinNTvU6ioiISFxSSYsxS6fmAMT1Kc+Wrj42HWjmLTrVKSIiMmoqaTEmLz2ZaXlpvBzHJe3J3XUEgo63LS7xOoqIiEjcUkmLQcumxvfigcd21lKUmcKS8myvo4iIiMQtlbQYtGxqDtWtPdS29Xgd5bT19Ad4Zk89axYV4/OZ13FERETilkpaDFo2LX7npW2obKSrL8Cli4q9jiIiIhLXVNJi0KLSLJITfGw52Ox1lNP26M4a0pMTWD073+soIiIicU0lLQb5kxJYMiWbTQeavI5yWoJBx+O76rh4fpGuMiAiIjJGKmkxasXMPLYfbaWnP+B1lBHbdqSF+vZeneoUEREZByppMWrFjFz6A46th+JnXtpjO2tJ8BmXzNf+aCIiImOlkhajlk/Pw4y4OeXpnOPRV2s4Z2Ye2WlJXscRERGJeyppMSo7NYn5xZlxU9L21HZQWd/JO8/QBrYiIiLjQSUthq2YkceWg80MBIJeRzmlv2yvxgzerpImIiIyLlTSYtiKmXl09gXYVd3udZSTcs7x0PZqzpmZR1Gm3+s4IiIiE4JKWgxbOSMPgBf2N3qc5OT21Hawr66Dvzuz1OsoIiIiE4ZKWgwryfYzsyCdDZWxXdJ0qlNERGT8qaTFuPPm5LOxqpH+GJ6X9tD2albO0KlOERGR8aSSFuPOn1NAZ1+Al2P0Op57atvZV9fB5Ut0qlNERGQ8qaTFuHNnFWAGz+9r8DrKsB7cdgyfTnWKiIiMO5W0GJedlsSS8mzWxWBJCwYd9287ynlzCnSqU0REZJyppMWB8+YUsPVQCx29A15HeYPNB5s50tzNVWeXex1FRERkwlFJiwPnzylgIOh4oSq2Vnn+cesRUpMSeNsineoUEREZbyppceDs6bn4k3w8u6fe6yiv6+kP8OdXqnnHGSWkpyR6HUdERGTCUUmLA/6kBM6fU8Dju+pwznkdB4CndtfR3jPAu8/SqU4REZFI8KSkmVmemT1mZnvDn3NPcNy3zOxVM9tlZt83M4t21ljx1oXFHG3pZndNbFwi6g9bj1KYmcLq2fleRxEREZmQvBpJ+xzwhHNuLvBE+P4bmNlq4DxgCXAGsAK4KJohY8lbFxQB8MSuWo+TQH17L0+/VsfapWUkJmgwVkREJBK8+g27Frg7fPtu4MphjnGAH0gGUoAkwPuG4pGiLD9Lp2Tz+K46r6Pw+5eO0B9wvHflNK+jiIiITFhelbRi51x1+HYNUHz8Ac65DcBTQHX44xHn3K7oRYw9axYWs+1wC3XtPZ5lCAYdv3nxECtn5jGnKMOzHCIiIhNdxEqamT1uZjuG+Vg79DgXmgn/ptnwZjYHWAhMAcqBt5jZBSd4r9vMbLOZba6vj50VkOPtrQtDXfap3d6Npq2rbOBQUxc3nKNRNBERkUiKWElzzq1xzp0xzMcDQK2ZlQKEPw/XOt4NbHTOdTjnOoCHgXNP8F53OOcqnHMVhYWFkfqWPLewNJOpean8ZXuNZxl+/cIhctOSeIcuAyUiIhJRXp3ufBC4OXz7ZuCBYY45BFxkZolmlkRo0cCkPt1pZrxrSRnr9jXQ0NEb9feva+/hsZ21/P3yKaQkJkT9/UVERCYTr0ra7cClZrYXWBO+j5lVmNmd4WN+D1QC24GXgZedc3/yImwsWbusnEDQ8dD26lMfPM5+++JhBoJaMCAiIhINnmwV75xrBN46zOObgVvDtwPAh6McLebNL8lkfnEmD247xvvOnRG19+0dCHDPxoNcOK+Q2YVaMCAiIhJp2uQqDl2xrCx8cfOuqL3nn16upr69l1vPnxm19xQREZnMVNLi0BVLywB4YNuxqLyfc447n6tifnEmF8wtiMp7ioiITHYqaXFoal4aK2fk8bvNhwkGI38tz3X7Gtld084HL5jJJL4yl4iISFSppMWpG1ZN42BjF+srGyP+Xj94ai+FmSmsXVYW8fcSERGREJW0OPWOM0rIS0/mVy8cjOj7bKhsZGNVEx+9aLa23RAREYkilbQ4lZKYwN8vn8JjO2upbu2OyHs45/ivx/dQlJnC9brCgIiISFSppMWxm1ZNxwE/e35/RF5/Q2UjL+5v4mMXz8afpFE0ERGR/7+9u4+tsrzDOH5dtBBs5wvvUYrAnIAVdUK3IUvEIMlgoiwTdVvmzLI/tmw6Z8yE/c7h/QAAB49JREFULdn0D90Lc0T34gxRlE2mMcVkBhF0aNREMSuIIiCCWASsgDCpVrGV/vbHOSA4h1S6c9+c8/0kTZ/zHPo8V3Knpxf381ZKlLSj2JC+NZp65on6+7Ovafe7Hd267b2doV89vFYnHt+bm9cCAJAAJe0o9/1zT1Fb+1799Znmbt1u4/LNenFrq2ZOGcUsGgAACVDSjnL1Jx2nSacN1JwnN2pXW3u3bLN1T4d+t2SdGob22X9PNgAAUFqUtDJw3eRRamv/QH96bEO3bG/2Iy9rZ1u7rr/wdO6LBgBAIpS0MjBi0LG6ZOwQ/W1Zs17Z8c4RbWvZxp26++lmfWfcUJ1Rd3w3JQQAAF1FSSsT135lhGp6VWtG4wuf+ikErXs69NPG5zWsX41mTBnVzQkBAEBXUNLKxMBje+sXU+vVtOnfuuvp5i7/fGdn6Jr7VqrlrT36/aVnqaZXdfeHBAAAh42SVkYuHjNYk04bqF8vWqtlGw//cVERoVlL1mnpS9v1ywvrNXZo3/9jSgAAcDgoaWXEtmZf9nmd3K9GP5y/Qi9ve/sTfyYidOvS9br9iVf0rS+drMvHDS1BUgAA8EkoaWXmuN49decVX1B1D+uS259RU/Ou//lv93Ts1YwFL+iWf67XxWPqdOO00VzNCQBAJihpZWh4/1o1/mC8TqjpqcvmLNNND63RttY9+9/v2NupRata9NU/PKX7m7boqomf06zpZ6pHDwoaAAC5cMSnuxIwVw0NDdHU1JQ6RhZ2v9ehGxeu0YIVW9QZ0mf716p3zyo172zTu+17Nbx/rW646HRNGDEgdVQAACqS7eUR0fCx71HSyt+rb7bpoRde15qWVr3f0am6PsdowsgBmjBioKqYPQMAIJlDlTTus1ABhvev1ZUTT00dAwAAdAHnpAEAAGSIkgYAAJAhShoAAECGKGkAAAAZoqQBAABkiJIGAACQIUoaAABAhihpAAAAGaKkAQAAZIiSBgAAkCFKGgAAQIYoaQAAABmipAEAAGTIEZE6Q7eyvUPSphLsqr+kN0uwHxw+xiRPjEt+GJM8MS75KcWYDI2IAR/3RtmVtFKx3RQRDalz4EOMSZ4Yl/wwJnliXPKTekw43AkAAJAhShoAAECGKGmf3pzUAfBfGJM8MS75YUzyxLjkJ+mYcE4aAABAhphJAwAAyBAlrYtsT7a9zvYG2zNT54Fke4jtx22vsb3a9tWpM6HAdpXt52wvTJ0FBbZPsN1o+yXba22fkzpTpbN9TfGz60Xb99runTpTJbI91/Z22y8esK6v7Udtry9+71PKTJS0LrBdJenPkqZIqpf0Tdv1aVNB0geSro2IeknjJP2IccnG1ZLWpg6Bg9wqaXFEjJJ0lhifpGwPlvRjSQ0RMVpSlaRvpE1Vse6WNPkj62ZKWhoRp0paWnxdMpS0rvmipA0RsTEi2iXdJ2la4kwVLyJaImJFcfltFf7oDE6bCrbrJF0g6Y7UWVBg+3hJ50q6U5Iioj0i3kqbCpKqJR1ju1pSjaTXE+epSBHxpKRdH1k9TdK84vI8SV8rZSZKWtcMlrT5gNdbRBnIiu1hks6W9GzaJJB0i6TrJHWmDoL9hkvaIemu4mHoO2zXpg5VySJiq6SbJb0mqUXS7oh4JG0qHGBQRLQUl9+QNKiUO6ekoWzY/oykBZJ+EhGtqfNUMttTJW2PiOWps+Ag1ZLGSPpLRJwtqU0lPnyDgxXPcZqmQoE+SVKt7W+nTYWPE4XbYZT0lhiUtK7ZKmnIAa/riuuQmO2eKhS0+RHxQOo80JclXWS7WYXTAibavidtJKgw+78lIvbNNDeqUNqQziRJr0bEjojokPSApPGJM+FD22yfKEnF79tLuXNKWtf8S9Kptofb7qXCyZ0PJs5U8WxbhXNs1kbE7NR5IEXEzyKiLiKGqfB78lhEMDuQWES8IWmz7ZHFVedLWpMwEgqHOcfZril+lp0vLubIyYOSriguXyHpH6XceXUpd3a0i4gPbF8paYkKV+DMjYjViWOhMGtzuaRVtlcW1/08IhYlzATk6ipJ84v/0dwo6buJ81S0iHjWdqOkFSpcqf6cePJAErbvlXSepP62t0i6XtJvJN1v+3uSNkm6tKSZeOIAAABAfjjcCQAAkCFKGgAAQIYoaQAAABmipAEAAGSIkgYAAJAhShqAimS7n+2Vxa83bG8tLr9j+7bU+QCAW3AAqHi2b5D0TkTcnDoLAOzDTBoAHMD2ebYXFpdvsD3P9lO2N9n+uu1ZtlfZXlx8HJlsj7X9hO3ltpfse4wMABwJShoAHNopkiZKukjSPZIej4gzJL0n6YJiUfujpOkRMVbSXEk3pQoLoHzwWCgAOLSHI6LD9ioVHge3uLh+laRhkkZKGi3p0cKjF1UlqSVBTgBlhpIGAIf2viRFRKftjvjwRN5OFT5DLWl1RJyTKiCA8sThTgA4MuskDbB9jiTZ7mn79MSZAJQBShoAHIGIaJc0XdJvbT8vaaWk8WlTASgH3IIDAAAgQ8ykAQAAZIiSBgAAkCFKGgAAQIYoaQAAABmipAEAAGSIkgYAAJAhShoAAECGKGkAAAAZ+g/fsick4b0lrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10,6))\n",
    "plt.plot(t, y_underdamped)\n",
    "plt.xlabel('Time')\n",
    "plt.ylabel('Velocity')\n",
    "plt.title('Underdamped Oscillations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, you can get the plots for overdamped as well as critically damped oscillations, by tweaking the parameters. It is encouraged you try it out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus so far, we have seen how to solve ODEs using solve_ivp. solve_ivp has many other functionalities that you can check out [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your Assignment..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. should you chose to accept it"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stars are a result of complex hydrodynamic, thermodynamic, chemical and nuclear processes happening inside them. For long astrophysicists have desired more and more accurate mathematical models in order to explain structure as well as processes going on inside stars. Most of these models are very complex, both in theory and also in solving them numerically. Refer [here](https://en.wikipedia.org/wiki/Stellar_structure) to know more. But there have been instructive yet simplistic models proposed which can be used to get a good overview about the physics of these celestial objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One such model is the [polytropic](https://www.astro.princeton.edu/~gk/A403/polytrop.pdf) model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The polytropic model can be obtained by combining hydrostatic equilibrium with the equation of conservation of mass and then assuming an equation of state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final equation, also known as the Lane-Emden equation is given by: $$\\frac{1}{\\xi^2}\\frac{d}{d\\xi}\\left(\\xi^2\\frac{d\\theta}{d\\xi}\\right) = -\\theta^n$$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here $\\xi$ is called the scaled radius and $\\theta$ the polytropic temperature. The solution of the equation for each value of n is called a polytrope. The following relations are obeyed for this model:\n",
    "* $P = K\\rho^{1+\\frac{1}{n}}$ . P, $\\rho$ being the pressure and density at each point respectively, n is called the polytropic index and K being the polytropic constant.\n",
    "* $\\rho = \\rho_c\\theta^n$. Here $\\rho_c$ is the density at the center of the star.\n",
    "* $P = P_c\\theta^{n+1}$. Here $P_c$ is the pressure at the center of the star. (You would have to use this to obtain the initial conditions).\n",
    "* $r = \\alpha\\xi$. r is the radial coordinate within the star and $\\alpha$ is called the length constant, given by $$\\alpha^2 = \\frac{K(n+1)\\rho_c^{1-\\frac{1}{n}}}{4 \\pi G}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your task is to numerically solve the Lane-Emden equation for n = 0,1,2,3,4,5 (Hint: You need to perform a task repeatedly. What do we do when we need to perform similar tasks repeatedly?) and make the plots of $\\theta$ vs $\\xi$ in each case. Before deciding to jump aboard, you should answer the following questions:\n",
    "* What are the dependent and the independent variables you need to consider?\n",
    "* What are the model parameters?\n",
    "* How many initial conditions will you need?\n",
    "* What are the initial conditions? Note that you can't supply arbitrary initial conditions, since the fundamental scientific constraint you have is that the star you have created must be physical !!! (Hint: What happens at the center of a star?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that since you don't know the radius of the star beforehand, assume a reasonable span of $\\xi$ on which you would want to integrate the model. You could try different such spans and then evaluate in each case. Also, since $\\xi$ appears in the denominator of the model, you might get ZeroDivisionError while solving, so instead of starting from $\\xi = 0$, choose a small value of $\\xi$ close to 0 to start with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After you have done this much, try to answer the following question:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you are given the pressure and density values at the center and the polytropic index of a star obeying the polytropic equation, how would you calculate its radius?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that its a second order ODE and thus needs to be converted to two first order ODEs before solving. Let $\\frac{d\\theta}{d\\xi} = \\phi$. Then the two simultaneous first order equations are:\n",
    "$$\\frac{d\\theta}{d\\xi} = \\phi$$ and $$ \\frac{d\\phi}{d\\xi} = \\frac{-2}{\\xi}\\phi - \\theta^n$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state would be represented by [$\\theta$, $\\phi$] and $\\xi$ is the independent variable.\n",
    "Since at the center of the star, pressure is maximum and is equal to the central value, we have the initial conditions as:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At $\\xi$ = 0,  $\\theta = 1$ and $\\frac{d\\theta}{d\\xi}  = \\phi = 0$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def polytropic(xi, state, n):\n",
    "    theta = state[0]\n",
    "    phi = state[1]\n",
    "    return [phi, (-2/xi)*phi-(theta**n)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_polytrope(index):\n",
    "    xi_span = (0.001,10)\n",
    "    xi_evaluate = np.linspace(0.001,5,1000)\n",
    "    init_state = [1, 0]\n",
    "    sol = solve_ivp(polytropic, xi_span, init_state, t_eval = xi_evaluate, args = [index])\n",
    "    theta = sol.y[0]\n",
    "    plt.plot(xi_evaluate, theta, label = n)\n",
    "    plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAGHCAYAAAATNz5KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3Rc933n/fed3gcDYNAbAXaxiLIkW7KaLcuW5JY42Wxix5uc3cSO/SROtp5s8uxmT7Il2d3HafskfnadrNOcOIlTNrYkS7JsSrZsVVaRBEAQIEESwKBM73fu7/njdzEzYKdIECD5fZ1zz+0zF6Atfvj9lWsopRBCCCGEEDeeY60fQAghhBDidiVBTAghhBBijUgQE0IIIYRYIxLEhBBCCCHWiAQxIYQQQog1IkFMCCGEEGKNSBATQtwSDMOYMgzjfWv9HEIIcTUkiAkh1h07VBUNw8gZhjFnGMaXDMMIXafP/knDML5zPT5LCCGulQQxIcR69WGlVAi4C7gb+L9v1BcbhuG8Ud8lhLi9SRATQqxrSqkzwNPADsMwPmIYxluGYaQMw/i2YRjbzr3eMIwuwzAKhmG0NR27yzCMecMwdgJfAO6zq20p+/yXDMP4fcMwnjIMIw+8xzCMbfZ3pOzv/EjT533JMIwvGIbxnGEYWcMw9hqGMdh0fqt9bskwjFHDMH6k6dyThmEcse87YxjGv1qVX5wQ4qYgQUwIsa4ZhtEPPAlkgT8HfgGIA08B/2AYhqf5eqXULPBt4EeaDn8S+Aul1CHgZ4DvKaVCSqmWpms+DvwnIAy8AvwD8CzQAfwc8GeGYWxpuv4TwK8B7cB+4M/s5w0CzwFftu/9UeD3DMPYbt/3B8CnlVJhYAfwwtv6xQghbgkSxIQQ69Xf2RWr7wB7gSPA15VSzymlqsB/B/zA/Re494+AH4d6M+OPAX9yme/7e6XUd5VSFnAnEAJ+XSlVUUq9AHzN/pxlX1dKvaiUKgO/jK6y9QMfAqaUUv9bKWUqpfYBXwX+kX1fFdhuGEZEKZVUSr15db8WIcStRIKYEGK9+gGlVItSalAp9VmgBzi5fNIOTNNA7wXu/Xt02NkAPAaklVKvXub7ppu2e4Bp+zuWnTznu+rXK6VywJJ93yDwTrtJM2WHyU8AXfblP4Su8J20mzTvu8xzCSFuYa61fgAhhLhCZ4GdyzuGYRhAP3Dm3AuVUiXDMP4SXRXbyspqmLrI5zcfPwv0G4bhaApjA8BY0zX9Tc8SAlrt+6aBvUqpxy74JUq9BnzUMAw38LPAXzZ/lhDi9iIVMSHEzeIvgQ8ahvGoHWL+JVAGXr7I9X8M/CTwEVYGsTmg79y+Zed4BSgA/8YwDLdhGI8AHwb+oumaJw3DeMD+nF8Dvq+UmkY3YW42DOOT9r1uwzDusTv/ewzD+IRhGFG7eTUDWAghblsSxIQQNwWl1Ci6wvW7wAI6GH1YKVW5yPXfRYecN5VSJ5tOvQC8BcwahrFwkXsr9uc/YX/X7wH/RCl1rOmyLwO/gm6SfIf9bCilssD70Z30zwKzwG8AXvu+TwJThmFk0AMHPnHlvwUhxK3GUOpiVXohhLi5GYbxAvBlpdQXr/Pnfgk4rZS6YXObCSFuTdJHTAhxSzIM4x70ZLAfXetnEUKIi5GmSSHELccwjD8Cngd+wW4qFEKIdUmaJoUQQggh1ohUxIQQQggh1ogEMSGEEEKINXJTdtZvb29XQ0NDa/0YQgghhBCX9cYbbywopeIXOndTBrGhoSFef/31tX4MIYQQQojLMgzj5MXOSdOkEEIIIcQakSAmhBBCCLFGJIgJIYQQQqwRCWJCCCGEEGtEgpgQQgghxBqRICaEEEIIsUYkiAkhhBBCrBEJYkIIIYQQa0SCmBBCCCHEGlkXQcwwjD80DCNhGMbhtX4WIYQQQogbZV0EMeBLwONr/RBCCCGEEDfSunjXpFLqRcMwhtb6OZZ968//mJNv7sfhdGA4HBgOA8PpxFFfXDhcTpwuF56AH384QqitlUh7B21d3UTbO3F7PGv9YwghhBBinVsXQexKGIbxKeBTAAMDA6v6XceefYlCYeYaP8UFODEMJwZOwIGBE8NwYGCHOwMMp4HD5cDpceLyefBFwsR6uundspWR3XcRjMWuw08khBBCiPXIUEqt9TMAYFfEvqaU2nG5a++++271+uuvr9qzvPbs1zjx2utYtRrKsvS6ZtnbFspSKMveNxWWaaJMhWVZqJoChb4GQCmUAqUPYlFDUQOqKFUFapd5GjeG4cGBW4c6w4HD4cDhcuDyewjEIrQPDrDprnsY2LULt1cqcUIIIcR6YhjGG0qpuy94ToLY2rFqNRZmpjk9NsbC6VNkEwlyi0kqmTxm0cQyLaya0qFOWShqWFRRqgxUL/KpPhyGF4fh0oHN6cIT8BLuaGPgjp1suuedtA5243TeNMVQIYQQ4qYmQewWU61WGT+yj8nXX2Vp6hTFpSy1QhWrqrBquvpmKROLMkoVAeucT3BgGDqwOR0u3F4vwViEruERhva8g+7tGwlGoxiGsRY/nhBCCHFLWfdBzDCMPwceAdqBOeBXlFJ/cLHrb/cgdqWKpTJjx/Zx4tXvkDpxinIyj1VUKNOBpRQ1VbPDWp7zw5oLh+HH6fLiD4SIdXfTt/0Oendtp2OwF28guBY/khBCCHHTWfdB7GpJELs+qmaNyVPHGdu/l4W3jlKaTWHlLKi6sCwnFoqaUcayckBlxb0GXpxOPz5/iEi8m67hEbq2b6ZrpJ9oRwcOp3NtfighhBBinZEgJq5axbQ4PXOW8dHvMHvkFYonE1gpBWUv1NxYykENk5pRRFlZVlbUDJzOEP5gC61dfcQ3jdCzZYTODf2E2+M4HBLShBBC3D4kiInrRinFXLrE1NQYUxPfZXHyVdTMPI6UG0cpDLUglnJhOixMCigrzcqBBQ7cngihaJzW/gHiw0P0bByira+XcFs7hmO9zDEshBBCXB8SxMSqU0oxmy5y/MRxTk28wtyZV2FuHG/KwFVohUoMLB8WTipOkxp5VC0FmPXPMAwn3kCMSLyHtoF+ukeGiA/209rbjz8ckcEDQgghbkoSxMSaSRerHD2T5PTEYZKnXiGVfhNHZhJftopKd+AqdOAxIyh8lNwOKo4KysqgrBTNc6w5XT5Cbd209vbRuWGQjg2DtHb30dLVhdPlXrsfUAghhLgMCWJiXSlVaxyZyXDs+AmSk2+Qm38NamP4ywmMTA0r04k/342/3I5yhCl4nZSdNSyVxqolQeUbH2Y4CETaaensIT44QMeGAVp7+qSKJoQQYt2QICbWvWS+wv7TKY5NTpM58RquhX2EHMcw3LPkixVyGS+uVDudmS7C5Q5MZ4Scz0vRDTWVQVlJVC1JcxXN5Q0QjfcQH+ynfWCA1t4+qaIJIYS44SSIiZuOUoqTiwVeP5nk6PEJClOv0Z45TNx3HMM3wyw1khkPaslF72KUzmwngUoX+UArWV+AgsdBjRzKSmLVllZU0QzDQai1g7a+Ptr6+3UFraeX1p4+/BGZyFYIIcT1JUFM3BIWc2Vem0ry2uQiJ0+MEpx/k12M0usdp+qbZ9TtYrbgoZz20J1Q9C1Eiee7cTi6yQY7SAdCFLwuXUGrJVHWku6LphoDBjz+IK29fbT19hHr1uEs1tNLS1cPLrdU0YQQQlw9CWLilpQvm+yfTvHq5BJvTJzBOv0mu9UYe5xjdPpOMOmpctjrYUL5qaQc9M3BQAIG03FCxU6K/i5ywW6yoRbybieWlbabOJeAFJaZrX+XYRhEOjrr1bNYt11F6+0nEG2RKpoQQoiLkiAmbgvFSo3XppZ4eWKRl4/Pkzk7xl3GGPe6JrjHO07BOMthr4dDPj/HfCHMBZOhOcXQnGLLvI+2TIyyu5NcsId8uJd8pJucYdnNm0mwkjgcKczKEspqzI3m8Qd0OOvpo7W7l9bePmI9fcS6enB5PGv4GxFCCLEeSBATt6V0ocr3Jxd5+fgCL08skkjMco9jlIe94zzsHaWlfJyjXheHfX4OR+O85TQwFopsmFVsTDjYvuine9aiqmLkg93kgr0UO0bI+rooWiZWTVfPDFI4nWlq5hLVUqrxAIZBpL3DDmm9tPb017dDsTapogkhxG1CgpgQQCJT4qXxBfaOzfPi+DzVQoa7HWN8qGWK+13H6M4dYcGwOOj1cqC1lwOBAG9VM0TSJhvmFDsXg2xf9NN1poAjXSUf6CYf7KbYtYVCbIisM0axquz+Z7p65nKlsawklfwCNbNcfxa3z0+su8du6uyzg1ofse4e3F7fGv6WhBBCXG8SxIQ4R81SHDyd4tuj8+wdm+fA6RQeVeFB/xQ/1DbFOzhKPHUA0yxxzOvjQMcwB8IxDlg5ZspJonnFpoSLezJtbFnwEJ/O4Z5ZoOIOkg90U+zcTKl7C7lgF1kzSKkMqBxWLYnTmcbjzQApKsUFSrklaPr/Ybgtrps3u3vt/mh6sEAkLu/pFEKIm5EEMSEuYylf4aXxefaO6mrZQq6C32Hyo10zfCQyzrbyfnyz+0DVSHj8HOzdzoFoBwco81b2FBWrgr+keEemlXszcTYmHLRNpzGmzqAsi6o7TKFjI+UNuynEhsi5WkkX3FTKFkpVUVYKtyeD15fFIE21skAxnaBaLtaf0eF0Ee3sItbVTUtXD7GuHlq6uol198jL1IUQYh2TICbEVbAsxcEzaZ4/MsfzR+c4NqtHT+7ucPDJ7jM85DpCfOEVjLnDAFS9EY717+FAaw8HnBYHslPM5GcACFteHq5s4J50K0OzFpGTi5hjx1HVKgqoRLuobr6HUs8W8sEeMlaIVMqiWqqhlAJVwOPL4Q/lcboyWGaKcmGB/FKCarlUf+bmkBbr7qGls4eW7h5iXd0S0oQQYo1JEBPiGpxaLPD80TmeOzLHq1NL1CxFPOzlo5vcfDR6gu3l/TinXoSlE/qGYAeJDe9mf3yIfS7Ynxzl6NJRaqqGgcGm8DAPmRvYvRShb6aM6/g0paPHUIWCvt/tRm29i/KG3RTbh8m62kgX3STnipgVC9AT3gYjJoFIAbcvh6HSmOUl8ukE6cQMZrnRH83pchHt6KpXz1q6evV2Vw/h9nYJaUIIscokiAlxnaQLVb41muC5o3PsHZ0nVzYJe108uq2Dj40o7jMO457aCye+BYVFfVPXLgrDD3M4PsQ+KuxbOMiB+QPkqjkA4v44d7bv5l3WIHcs+Gk9laJ6bJTSW0eopdP6MxwOPCPD1LbcRal3O/lQL1kVJjlfJjlboFa16s8YinmJtlt4AzkczixWNUkpP092fpbk3IVDWrSzS4e1c9ZunwwcEEKIayVBTIhVUDZrvDyxyDOHZvnGkVlShSpBj5NHt3Xy5I4O3hOZwXvy2zDxAky/ApYJ7iBseJDa8Hs43rmJ/eV53kzsY39iP2fzZwHwu/zsbN/JnfHd3M0QIwkDY2yS0pGjlI4cwUwk6s/gHhjAu3Ub5shuivFhcp44qbRi6Wye1FyBmtkIaP6wm5bOAKFYDY83C6R0FS2VIJNIkJqboVIsrPgZA9EWop1dtDSHtY4uol1dhFpaMRyOG/GrFkKIm5oEMSFWWbVm8f0Tizx1aIZvvDXHUr5CwOPkPVs7+ODObt4z5MN/5mUdyo5/E5KT+saWARh5FDY+ymznVvanJ9iX2Me+xD5Gk6NYytLNmbFN7OnYw56OPex2DhI9uUT5yFFKR3U4q05P15/F1dmJb9s2PNu2UxvcRqFlgEzFS2q2QHI2T3K2QLnQeK2T2+sk1hUg2uknHAOPNw9kqJSWyM7PkU7MkpqbI7swj1KNYOd0u4nGO8+rpi0HN6mmCSGEJkFMiBvIrFm8Mrlkh7JZPQLT7eR92zv5gTt7eHBTHE9mSgeyiW/B5ItQyYLhhIF3wabHYNP7yccGObhwiP2J/exL7OPA/AEKpq5YdQQ66sFsT8ceRpxdmKPHKR05QunoEUpHjlA5MQmWDk7OaBTfHdvxbd+OZ+s21NAW8q4YqbkiS7MFUnZAyyUbzZYOh0G0w0+sK6iDWtyD21cEK00+NU9qbpZ0Ypb03NxFq2mReAeReCeR9jjReKe9rxePz3/D/kyEEGItSRATYo3ULMWrk0t87eBZnjo0Q7JQpSXg5oM7u/nonb3cPRjDoUyYfhUmvgnjz8LsIX1zpK8eytjwEKbbx3hyvF4x25fYx1xhDoCAK8DO+E4dzOJ72BXfRaDmpDw6SvGIDmblI0cpjY9DVb+eyREI4N22Dd/27fi2bcN3x3aMngFSi1VSs3mWZgskZ3RAS88XUVbjvxWhmJeWzgAtHQFaOgNE4j4CYYtaNUVmYY70ckibT5BdSJCZT1AzzRW/G1840hTQ4nZgawQ1XzB0Y/6QhBBilUkQE2IdqNYsXhqf5+/2neW5I3MUqzV6W/x8eHcPH72zh23dEX1h5iwcf16Hsolv62qZ0wOD74bNH9DBrG0EgJncDPsS+3gz8Sb7E/sZS46hUDgMB5tjm7kzfme9atYd6kZVKpQnJnTl7K0jumnz2DFUUc9XZng8eDdvxrdtG96tW/Bt3Yp382bwB0nPF3XT5kyB5Fye1FyR1FyBSrERsBxOg2jcT9QOaC0dfmJdASLtPpQq1EPZckBLz+v9zHwCs1Je8fvyBoIrKmiR9o56VS3c1o4/EpXXRAkhbgoSxIRYZ/Jlk+eOzPH3+8/w4vgCNUuxpTPMR+7s4Qf29NLbYjfbmRU49T0dysafg4VRfbx1WAeyTe/XAc2t+2PlKjkOzh9k37yumB2cP0jR1CGrM9DJno493Nmhw9nm2GZcDheqVqNy8qQOZkfscHb0KNbyiE3A3duLd+tWfFu26IC2ZQvu/n4wDEq5Kqm5AqlEoR7OUokC6URxxWABt89Zr6C1dPhp6bIrah0B3D4nxWyGTGKOzIqANkdmYZ7M/ByVYmNyW9B91MKt7YTa2gi3xQm3tTet2wm3x/EFQxLWhBBrToKYEOvYYq7MU4dm+Lv9Z3njZBLDgHePtPOP7u7j/du78Hua5vlKTulANv6s7ltmlsAdgA0PN5oxW/rrl5uWyWhytN7PbF9iH4mCHnUZcAXYFd9VD2e747sJuoOAnqfMnJ2lNDpK+dgopdFjlI+NUjl5st7vzBEI4N28uR7MvFt09cwZ0p9hWYrcUqkR0BIFHdLmCmSXStD0n55AxFMPaJG4n2g8oCtrcT8evwulFOV8nvS8DmrZhQWyi/PklhbJLs6TXVwgt7SIVaut+N26vN6V4WzFtt73BgKr8KcqhBANEsSEuEmcWizw1TdP89dvnOZMqkjY6+JDu3v44Xf0cddAy8rqTrUIky/Z1bJvQOqUPt6xXTdhbn4C+u6GpglblVLM5GfqoexKmjObWcUi5ePHKR07Rnl0jPKxY5RGR7Gy2fo17oEBfFs2492yFd/WLXi3bsXd27vi2c1qjfR8kXRzQLPXxWx1xXf6w24i7X6iHX6i7f5602ek3Y8/7K5/rmXVKKRSZBcX6uEsuzhvh7YFsksL5JPJFSM/ATz+QL2CFoq1EWptIxRrJdTaWt8PRKIyVYcQ4m2TICbETcayFN+fXOSv3zjN04dmKVZrDMeD/PA7+vjYnj66oudMDaEULIzrUDb2jG7OtEwItMGmD8CWx2HkveANn/ddl2rO7Ap2sSfeaM7cFNuEy+E656sV5tmzlEbHKI8eo3RslPKxY1ROnaq/zNwRDOLdskUHsy1b8W7ehHfTJpyh8zvkV0om6fkimfmiDmv1xR7V2fSfLLfPqYOZHdR0YNPVtFCLF8OxslmyZprkU0v1ipoOa43tXHKJQip1XlgzHA6CsVYd0GJtK0Ja/XhrG95AUJpChRDnkSAmxE0sVzZ56uAMf/XGNK9NJXEY8OCmOD9ydz+Pbe/E47pApaaY0h3+x57RTZmllO7wP/SArpRteVzPYXYBpmUylhxbMTrzSpozz2UVCpTHx3UwGz1mB7VRrFyufo2rpxvvpk34Nulg5t20Cc/wMI6LzEFWq1pkFu1gliiSXmgEtsxCEavW+O+Zw2UQbddNnZE2P5F2H5E2P+E2H5F2H96A+8LPXauRTyfJLy2RTS6SX1oil1wit7RILrlI3t4u5XPn3evyeOuhrDmghWKtBGOtBFtiBFtiePwBCWxC3EYkiAlxi5hayPPVN0/z1TdOczZdoj3k4R/d3c+P3tPPYNuFAxE1E6a/D6NP62C2eFwf77hDB7LNT0DvO+AiTW+Xa87cEttSr5jt6dhDV7Dros+vlKJ65gzl8XHK48cpj41RHh+ncuIEyp5WA4cDz8BAPZgtV888g4MYLtdFP3u5T1p6QYe0zHyxsb1YpFpa2X/MG3DpUNbmJ2yHtEibr77t9l76HZzVStkOaYtNQU2v80n7+NLSeaNBAVxuD4GWGMGWFgJRvV4OaYGWGMHo8nYLbo/3ks8hhFj/JIgJcYupWYoXx+b58qun+ObROSwFD2xs5+PvHOB92y5SJVu2cBzGnoZRuwlT1SAYbzRhDr8HvJeew+tqmjM3xzbjvMyLxVW1SuXUKR3QxsbtoDaumzftwQGG241neLgR0OyQ5u7puWz/LaUU5YJJZqFIZqFEdrFEZrGo1wt6bVZXNkf6w27CF6ikRdr8hFt9ON2X7zOmvzdPbmmRQjpFPpWsL4VUknzTsWI2U2/KbeYNBOuhrRHQYvXgtrz4wxEcTnmBuxDrkQQxIW5hs+kSX3ltmq+8durKq2TLCkt6hv+xp2H8eSinwemFDQ/C5sdhyxMQ7bvsMzQ3Z+5P7OfNxJvnNWfe2XEnu9p3sSu+i6g3ekU/m1UqUTlxoh7MSvbaPDtTv8YIBPBu3Ih3eBjPyDDekY14R4Zx9/VhXGEwUUpRzFZ1UKsHtEZIyy6VVjR7YkAw4iHU6iPc6rPXXkIxvR9u9eENuq6q+bFmmhQzaR3M0sth7Zzwlk6ST6XOe4vBMl8oTCASJRBtwR+JEIi0EIhG8Uei9e1ARO/7Q2EZgCDEDSJBTIjbQM1S7B1L8OVXpnnhWKNK9ol3DvDY9k5czsv8pVur6grZ6DM6mC2d0Mc7dzaaMHv2XLQJs9mFmjPHU+NYdif4ocgQu+K72B3fze74bkZaRs4bBHDJR81mKR8/3mjiHB+nMjGBOT9fv8bwePAMDeHdOIJneATvyDCe4RE8G4ZweDxX/F2gmz3zqXK9kpZZKJFb0gEtlyyTXSpRO6ei5vI4GiEt5iXctryt16GYF+elKpeXUC2XVlbYkkkKmTTFbJpCOk0hk6KQTlPMpCnmshestBmGww5rUR3WwjrALQc1Hdpa6tsyEEGIt0+CmBC3mZl0kb987XS9StYd9fHj7xrkH9/TT3voCvocLY/CXG7CnP4+KAtCnXqusi1PwvAj4LnyObgK1QJvLb7FgfkDHJg/wMH5gyyVlgDwu/zsaN/B7vjuetWszd921T93LZOhPDGhq2gTJ6hMTFCemKB65kwjjDgcePr78Yw0wtnyenkOtKullKKUq5JdDmdLZXut97PJMsVMZeVNhp4/Ldzqsytp3npAC7Z4CbX4CETcOC4XoC/DqtUoZjMUM2kKmTSFdIpCJkPRDmsF+/jyfrmQv+DnOJwu/JEI/lAYfziCL6zXy4svFLbPR+rnJbwJoUkQE+I2VbMU3zw6xx99b4rvHl/E43Twod3d/MR9Q+zub7nyDyos6dGXY0/rpsxyBlw+PZHslid0M2ak+/Kf00QpxencaQ7OH+Tg/EEOzB9gdGkUU+lXJvWF+tgV16HszvidbI5txu288EjHy7GKRSpTU5QnTlCeOE5l4gTlExNUTp6qv3sTwNXVZTdxjuAd3oBnaAjPhg24OjuvOVCY1Rq5ZNkOZ2VyyVJTWNPHz+2nZhgQiNrBrB7QmtYxvXZ5rl/fMLNarVfWVoa3NMVMhlIuo4NdNksxm6GUy6Is64Kf5XA6dUBbDmrhCP6mAOerB7lwfd8XCEqTqbjlSBATQjA+l+WPv3eSr755mkKlxp39LfzE/YM8ubMbr+sq/iI3K3Dyu3oE5ujTkDqpj3ffqStlWx6Hrl06RVylklniyOIRHc4WDnIgcYBEUfc18zq9bG/brqtm8V3sat9FZ7Dzqr+jmapWqUyfXhnOJk5QPnGi/v5NAMPvxzM4iGfDkG7u3GCHtKEhnJHINT1D/VnsqlouVSafKpNL2utz9pvf7bnMG3TZAc1HqMVjBzffigDnDVxdn7Urfm7LolwoUMxl7KCWtYOaXkp2YDv3/LlvQVhmGA58oRC+UBhfMIQvFMIbtPdDIXxBe20f94fC+nwwhOsqm5yFuFEkiAkh6jKlKl994zR//L2TTC7kaQ95+LF7B/jEOwfPnyj2cpSCxNFGE+bp1wAFkd5GZ/+hB+vvwrxaSinmCnMrmjOPLB6haukqVmegk13xXexo38GOth1sb9tOyHPpEZ9X9L2WhZlIUJmaojI5qatpk5NUpk5SPX26PpITwNnWZoeywUZA27ABd3//VfdFuxKVknnBgNYc3IrZyoqJbwGcbgeBiIdg1Esw6iEQ9RKIeurbwaiHQMSLP+Q+byLc600pRaVYtIOaHdpy2RUVt1IuRymfs9dZyrkcpUL+gv3dlrk8XnzBIL7lcNYU5nzBEN7Q+cf0tUEclxnZK8S1kCAmhDiPZSleOr7AH788xQujCZyGwYd2dfNTDw6zo/fKRjWeJzevX7c0+jRMvADVAriDMPIeHco2fQBC8Wt67kqtwrGlY/XmzEMLhziTOwOAgcGG6AYdzNp3sLN9J5tjm/E4r18gUpUKlelpHdLqAW2KytRJagsLjQsdDty9vfUqmmdwEM/AIJ6Bfj3lhvvtNbNeiZppkU+XyacqdkgrkU9XKKTL9XUhU6FcOL+65nAY+COelWEt0hTW7LU/4sF5jf3XrpZl1agUipRyWTukZethrZzPUcxlKZ9zvJTPUc7lqJZLl/xsjz+ANxjEG1heAniDocZ+0D4WCK3ctq9zreKfp7j5SRATQlzSqcUCX3p5iq+8dop8pca7hlv56QeHec+WDhxvtzpSLcHUSzqUjT4N2bOAAX33NEZhdmx7W02Y50qWknubvJQAACAASURBVBxeOMzhxcN6vXC4PhDA5XCxNbaVO9rvqIezocjQZec2eztqmQyVkydXVtLskKYKTVNOOJ24e3rw9PfjHhzA0z+AZ3AAd38/nv5+HH7/dX+2CzErtZUBLdMU1NKV+rFz3/8JgAH+kJtAxA5mYR3OAmEPgYi7sR/x4A9d+6CDa1UzqzqY1Stt54e2ciFPuVCgXMjpdT5PxT527muvzuV0u+3AZoezFQHuMtuBAG6fT6pytzAJYkKIK5IuVvnKa6f439+dYiZdYjge5J89sIGP7enDfy0dwpWC2YO6+XL0KZjZr4+3DDY6+w++G1zXp3KllGI2P8uhhUMcXjzMWwtv8dbiW+SrekRgwBXQwaxtR7161h3sXrURfkopzPl5qqdOUTk1TeXUSaqnpqmcOkVlehornV5xvaujA/dAf6OC1t/YdkbfZrXyGtRqFsVM5byqWj5TobAc3LIVipkqNfPCgcUXdNvBTIe0QHNQq++7CYQ913XwwfWgLItKqWQHNXvJX8F2oUA5n6NSKGBWK5f9HrfPj9fvx+0P4PX78fgDulIXCNS3PfZxr9+Pp+m4t+mc0+2W0arrjAQxIcRVqdYsnjo0wxdfmuTQmTSxgJsff9cgn7xvkI7w2+vvtULmrN3Z/xmY3AtmCbwR2PiorpRtegwCrdf+PU0sZTGVntLhzK6ajSZH6/3NWn2t7GjfwR1td7CtdRvb2rbRGbj20ZJXopZK6ebOU6eawprebp4bDcAZjeIeGNDVtL4+3H29ugm0rw93dzfGGnZYV0pRKdUoZip2MKtQzFYoZCoUslWK9rGCfbxSunCHfbfPqYNZ2A5qEQ/+sBt/SFfXfGE3/pDe94Xcb3s+thvJrFQuWHErF3JUikUqxQKVYoFyobG9fLy8vF8oXrYyB3qqEU+gOcz5Vwa2gD7m9Qdw22uPz4/b58ft8+Hx+fH49bbL45VQdx1IEBNCvC1KKV6dXOKL35nk+aNzuB0OPnpnD59+eISNHdfeKR6ASh5OfNt+F+Y3IJ8AwwkD77I7/D8J7Ruvz3ed+9W1CuPJ8RXhbDIzWZ94ttXXyra2bWxv3c72tu1sa9tGT7Dnhv7FZBUKVKZPU50+ReXkKSrTjbBWnZ0Fs6mfl2Hg6uzE3deLp7cXd28f7t5e3H19ePp69TQcl3hf541mVmo6sGWrjfBmB7diU3grZCqU8tXzBiAs8/ic+EK60uYL6ZDmWw5t9r4/7NFVubAbj391RpCuNqUUZrlsB7Om8FYsUCmcG96aAl1BX1MtFSkX9DUXegfqhRiGww5nvkZQ8/t1cPPqY3rftzLI2dvnnvf4/LdlxU6CmBDimk0u5PnD70zyV29MUzYt3r+9k88+svHq5iO7HMuCs282XlA+d1gfb9vYCGX97wTn6oWJQrXAWHKMI4tHOLp0lKOLR5lITdTnN4t6o/WK2fZWHc76w/04jBtflVGmiTk3R+XMGaqnz1A9c4bq6dNUzpymeuYs5uzsylGGLhfuri5dSevt0VU0O6i5e/twxdvX7RxelqUo56sUs1VKeTu85aqUchWKueXjOriVcvrcuW87WOZwGPiWQ1rYjS+oK27Loc0b0Nu+gBtfyIU34Mbrd636aNIbqWaaVErFRoArlaiWilRLJX3c3q6WijrcLZ8vl+z9xrXVkj5/qRGtzQyHox7UlitxOqj58Ph1fzm3d3nxNu17L3HOh8vjWbf/+5UgJoS4bhZzZb708hR/9PIUmZLJ/SNtfPaRjbx7Y9v1/1du6lTjlUuTL4FVBX8MNj6m+5ZtfBR8q99nqlwrM54c58jikXpAG0+O15s1Q+4QW1u36nDWtp3trdsZjAyuyoCAq6EqFaqzs3Y4s8Pa6dNUz5yhcuY0tfmFlTe43bg7O3F3d+Pu6cbV3Y27uwd3T7c+1t2NI/j23j5woymlMCuWDmbLAW45tOWqOqytOFelVLh41c0w0IEs6MIXdOMLnrNthzZfwI3XPuYLum7a6tvV0r/vMpVic0BbGdSWQ1213HS+2Hy+RKVU0AGwXKJaKl9RU2wzl9d78RDn9dX39XX62PBd99La07tKvxlNgpgQ4rrLlU2+/MpJvvjSJIlsmZ29UT7zyAgfuKML52pUDkoZPSXG6NMw/iwUl8Dh0p38lyeSjQ1d/++9iGqtyvHUcY4uHa2Hs9GlUco13eTjd/nZFNvEltgWvbRuYVNsE0H3+gkyVqlE9ezZeiWtenaG6szychZzLgHnTLzqiEbroawe2Lq662HNFY+vq+bPq2HVLMoFk1K+qtd2OCvn9bHlpZyvUrKPlfPVi/Z1AzAcBt6Aqx7MvEG70tYU5Dx+F96AXXkLuOqLy72+Bi3caEopatWqDmV2MKtvn7tfKlEt633znP2L3WvVdJX7w//8F9n8rgdW9WeRICaEWDVls8bfvnmGL+ydYGqxwHB7kE8/PMwP7unDs1qdqK0aTL+qR2COPQMLY/p4fJuulG15AnrfATe4ImVaJpPpyXo4G10aZTQ5SraSrV/TH+6vB7Pl9WqO2LwWyjT1aM+ZGTukncWcmW0KazPnjfjE6cTV2aErad3duLs6cXV26WOdnbi6unC1td20Ye1CajWLct6kXFgZ0BrBrRHkmgNe9RIBDsDpcqwIZh7/yqDmbdr3BHQ1bjnUefyutz/1zG2iZppUyyVcbs+qv5VBgpgQYtXVLMUzh2f5vW8f562zGboiPj798DA/du8AvtX+l/3iROOVSydfBlWDQLvdr+xxGH4PeK/T4IKrtDyVxmhytB7MRpdGmc5Oo+x2sLAnzObY5hUBbaRlBJ/rOoxQXWW1XB5zduacsGZvnz1LNZFY8T5PABwOXO3tuLq6cHd24OroxNXVqYNaR6cd3jpv2Hxqa2U5wFWKJqVClUrBpFwwKRd1qCsv7xf0vr7O1NcVTZR1ib+/DfD4XHj9LrxBex1w4wnobY9/ee3E49P7eu3U234XLrdjXf4D4WYkQUwIccMopXhpfIH/8a3jvDq5RDzs5WceHuHj9w5c21xkV6qYhPHndb+y8eehnAanFzY81JhINrq6/UGuxPKggLHkWD2gjSXHKJr6HZdOw8lQZIjNrZvZHNvMxpaNbGzZSE+oZ00GBrxdyrKoJZOYc3NUZ+cwE3NU5+YwZ+f0scQc5lwCK5s9715HJKKDWnNFrXk7HsfZ2rpuO2ivJqUU1XKtHtQqRV2NqxQbwa0R6hr7y2HOLF+6Ggd6UIPH3xTOlgNbU3jz+l14fI3w1gh1jWNSmZMgJoRYI98/schvPz/O904s0h7y8umHhvnEuwYIeG5Qs1StqitkY/ZEsskpfbxrV2Mi2e47YZ38RW4pi+nsdCOYLY0xmhxlJj9Tv8bv8tdD2caWjWyMbWRTyyba/e03dfXCyuepziV0UJudxZxLNIKaHdrMhYXzR+Y5nbja2nSFLR7H1RHXAc3ed8ft/Xh8Vd79ebOyahaVUo1K0WysiyaVkh3mzjteqx+vluzzxdqlq3I2l9eJtzms+V14vE7cPiduX2Pb43Pptddln2tse+xrb9ZQJ0FMCLGmXp1c4ndfGOel8QVagx5++sFh/sl9gwS9N7CfkFIwP9p4Qfn0K4CCYFyPwtz0Phh5rx6Vuc5kK1kmUhMcTx1vLMnjLJYW69dEvdGVAa1lI5tim4h6b/xM/KtFVauYCws6oM0lMBfmMeeblwXMhXlqi0srXsy+zBGN4orbgW15aY+v3O+I4wgGb+pQe6MopTCrViOsFWv1MFc+J8CtDHk1qiWTarlGpaS3rzSKuNyORoDzOXF7mwOcPn6hAHdu2HN79XKj/pwliAkh1oU3Tib5nW+Os3dsnljAzU/ZgSzsW4MXJucX4PjzMP6cXpdSYDj0PGUb3web3g9dO6/LuzBXy1JpiYnUBOPJ8RUBLVttNPPF/XE2tmxkpGWEkZYRNkQ3MBwdJuZbf4HzelGmibm0hDk/T21h4ZywZgc2e1tVzn/1kOH342prw9nWiqu1DVd7G87WtsaxtuXtNpwtLbdl0+j1tBzoqqUalZJJtVSjWjbtkGYfawptlbI+Xi3Z15Qb91XKtStqdgXAALfXyaM/sY2RPR2r+jNKEBNCrCv7p1P8zjfHeeFYgqjfzaceGuYn7x+6sRWyZjUTzryhp8U4/hzMHNDHw912KHsMhh+5IXOWXSulFIlCoh7MlkPaRGqCUq1Uvy7mjbEhuqEezIZbhtkQ3UB3sPum6oN2LZRSWJmMDmXLgS1hh7WlRWoLi5hLS9QW9frcqTwAcDhwtrbiam09J7BdOLg5vN4b/4PeZixLYS4Ht3qga4S2cwPc1vu6ifeHV/WZJIgJIdalQ6fT/PY3x3j+aIK2oIfPPDLCj79rcPVHWV5OdrZRLZv4lu7w73DBwH2NalnHtnVdLTuXpSxm8jNMpic5kTrBifQJJtOTTKYnSZaT9ev8Lj9DkaGVIS06zEBkAI/z9u1jpSyLWjqtQ9niErXFBczFJczFBWqLSzqwLSzoStziIqpQuODnOEKheqXNGYvhjLXgisXs7dZz9mM4wmFpJr0FSBATQqxr+04l+fxzY7w0vkBXxMfPvncjP3J3/+rNQ3Y1alU9Z9nx53QwW37tUqRP9yvb9H7Y8PCaTY9xPSRLyXowO5HWIW0qPcWZ3Jn6NU7DSV+4jw2RDQxGBhmMDjIYHmQwMkhHoEPCwjmsQqFRTWsObIuL9QpbLZmsL+rcKT6WuVw6nLU0wpmzNabDWss5+8vhzbf+pz253UgQE0LcFL5/YpH//o1RXj+ZpL/Vz88/upkfuLMHl3MdBLJl6TN2texZ/bLySg6cHv2S8pH36qVz57oZiXktimaRqfTUipA2mZ5kOjtdf4MA6CraQHiAgcgAQ5GhFeuYNyYh7TKUUlj5ArVUkpod0MxkkloyVQ9qZnJpxX4tlbroux0Nv9+urLXibGnBGY3ibIniiEb1drQFZ0vzOoozEsFwr0FfzduEBDEhxE1DKcXesXn+n2fHOHQmzXA8yD9/32Y+uLN7/Q1dNytw6nu6Wnb8BUi8pY8H2mHkPTDyqF6Hu9b2Oa8zS1nM5ec4mT3JyfRJvc7o5Uz2TP0F6aAnqx0MNypozSEt7Fndfjm3MlWrUctkVlTVzGSS2lJyZXhLp7FSad2smslccDTpMkcwiDMaxdHSFNii0UaYswPd8rbDPifTglyeBDEhxE1HKcWzR+b4/LNjjM5l2doV5l9/YAvv3bqOm8EyM7pKNvECnPgW5Of18Y477GD2Xhi8H9y37ozxVavK2dzZejBrXmbzs/W3CYAeMNAX7tNLqI/+cD99Yb3uCHTcNoMGbhRlWVjZrA5l6TS15YCWTunAlk5TS6WajjeWCw5UsBl+fyOohcM6oIXDOCJhnOEIzkgYRziCMxrBEQ7jjETs8xE9VcgtUD2+HAliQoibVs1SfO3gWX7r+XEmF/Lcu6GVX3xiK3cNrPPpFyxL9yebeEEvp74HtYqe5X/w/qZmzDtuqk7/16JklpjOTnMqc4qpzBSnc6eZzk5zOnua2fwsNdX4y97tcNMb6l0RzpbDWm+4F7/r1g2z641uOs3rgJZKUUundGg7N9ClUjroZbPUMhmsTAYrn7/0hzscOpxdKLhFIhc4pgPccpgzAoH1+w+zJhLEhBA3vWrN4iuvTfNbz4+zkCvzxI4u/tUHtjASv0k6yVcKepb/5WA2f1QfD3Xqd2EOPwIbHoRo31o+5ZqpWlVmc7NM53QwO521Q5od1vLVlX+ht/vb6+GsL9xHT6iH3lAv3cFuOoOduB3S32k9UKaJlcutCGe1TBYrq9e1bAYr03Qu23wue9HRp3UuF85gEEc4rANdKKS3Q0GcIftYOIQjFMYRDunAFwrp0av2PY5AYNWrchLEhBC3jHzZ5IsvTfI/X5ygZFr86D39/Pyjm+iI3GQjxTJn9dQYy82YBXuW/NYR/V7MDQ/C0EMQiq/tc64DSilS5dSKcNa8PZefW9Hk6TAcdAQ66An20BPqoTvYrUNaqLse1m7nqThuJqpa1eHMDmm1TEZX3ZpDXS5LLZvTx3NZrFze3tbHLtUvDqDnv/03oh/+0Kr+HBLEhBC3nIVcmd/95jh/9sop3E4HP/3gBn76oeG1maX/WlkWJI7A5It6OfldKGf0uY7tdjB7CAbfDf6WtX3WdahSqzCbn+Vs/ixnc3qZyc9wJneGmdwMc4W5Fc2eoN840B3qpje4MqD1hnrpCnYRcAfW6KcR15NSClUoUMvlsXJZu+k0p8NbLoeVzRF65BG8wxtW9TkkiAkhbllTC3n++7OjfO3gDG1BDz/33o18/J2D62MOsrerZurZ/Sf36mB26vtgFvUrmLp3N4LZwH3gCa710657pmWSKCR0MGsKaMvBbSY/g2mZK+4Je8J0BjrpCnbV183bnYFOCWviikkQE0Lc8g5Mp/j1p4/xvROLDLQG+MUntvLEjq6boiPvZZllOP16o2J2+jWwquBwQ+9duvP/4Luh/96b4jVM642lLOYL8/WQNpufZa4wx2x+tr69VFo6776IJ3LRoLa87XPdZE3mYlVIEBNC3BaW5yD79aePcWw2yz1DMf7dh7azq+8Wa86r5HWVbHKvHgBwdh9Ypq6Yde7QoWzwfr0E29f6aW8J5VqZRD7BbGF2RVCby8/Vt5tfFbWsxdtCR6CDeCBOh/+cdaCDuD9Om78Nl2ON3rMqboh1H8QMw3gc+G3ACXxRKfXrl7pegpgQ4lJqluIrr03z+edGWchV+NieXv7141vojt6iUx5U8rpidvJl3b/s9Gtg2i/4bt/cqJgN3n/bjsq8EUpmiUQhoStpBR3SZvOzJAoJEsUEC4UFFkoLWGpl53EDgzZ/G3F//PzQFuig3d9OR6CDmDeG07HG72EVb8u6DmKGYTiBMeAx4DTwGvBjSqkjF7tHgpgQ4kpkS1V+79sT/MF3JnEY8OmHRvj0w8MEPLd49cGswMx+HcpOvgynXtEvLgeIDsDgfdB3j27K7LgDnLf472MdqVk1lkpLJIoJ5gvzJAoJ5ovzK7YThcQFm0KdhpM2f1s9pLX522jztdHubz9vO+C6OebXul2s9yB2H/AflFIfsPf/LYBS6r9c7B4JYkKIqzG9VODXnznG1w/O0Bnx8q8/sJWP7eldf69MWi1WDebealTMpl+B3Jw+5w7qfmb990LfvTqgBdvW9nkFVavKYnFRB7QLhbZigsXiIslScsXUHcv8Lj+tvlba/G20+3Q4a/e30+ZrO29bBh2svvUexH4YeFwp9VP2/ieBdyqlfvZi90gQE0K8Ha9PLfFrXzvCgdNpdvZG+Xcf2s69G1rX+rFuPKUgdUo3YU6/AtOvwuwhWJ7ioXVEB7PlcNaxDaRJbF0yLZNUOcVicZHF4iILpQW9Li6wWFqsby+Vli4Z2poDWquvlZgvptfeGK1+e+1rpcXXIpPlvg23RBAzDONTwKcABgYG3nHy5Mkb/qxCiJufZSn+/sAZ/uszo8ykSzy5s4tfenIbfbHbvCpQKehO/6df1cFs+lUoLOhznjD03Ak9e3T1rOcuaBm4bV7NdKswLZNkKbkioC2W7HVxsX58qbREqpw6ry/bsrAnXA9p9cDmi523v7z2Or03+Cddf9Z7EJOmSSHEDVes1PifL57g9/ceB+AzD2/k0w8P43NL5QfQVbPkJEzbVbOzb8LsYT1tBkCgTQeznrsaAS3ctbbPLK6bmlUjU8mQLCV1Na2crIe4ZClZX5bKutKWKqUwlXnBzwq4AucFthZvCy2+FqLeKFFPlBav3l4+fquFt/UexFzozvqPAmfQnfU/rpR662L3SBATQlwvZ1JF/vPXj/L1QzP0xfz8+w9t57HtndLR+ULMsu5rdvZNXT07s0+/M3O5chLusUPZHujeA107Idy5ts8sbgilVD24JcvJelPouUFueT9dTlOqlS76eX6Xvx7Mop5oY7sprLV4W4h4Inrf20LYE163o0rXdRADMAzjSeC30NNX/KFS6j9d6noJYkKI6+3l4wv8h394i7G5HA9tjvMrH95+87xQfC1V8rp/2Rk7nJ19ExaPN84HO3Qg69oBXbv0dttG6XMmKJklUuUU6XKadDlNqpyq7y9vZ8qZFcfTlfRFm0wNDCLeSD2wNVfaIt5IPdBFvVEinkj9mhsR4NZ9ELtaEsSEEKuhWrP4k++d5DefG6Nk1vin797Azz26iZBXpne4KsUUzB3WAW32MMwehMTRRrOmyw+d2+2AtlMHtI7t4JXgKy7NUhbZSvaS4a353HJ4y1fzF/3M33jwN3hy+MlVfW4JYkIIcRXms2X+6zPH+Ks3TtMR9vJvn9zKD9zZK82V18KswMKYHc4O6XA2ewhKqcY10QHo2KpHaca36e32LeC5zQdSiGtWtaorAlymkqmvH+h9gMHI4Kp+vwQxIYR4G/adSvIf/s9bHDid5u7BGL/60R1s74ms9WPdOpSCzBkdyOYOQ+KYrpwtjkOtYl9kQGywEcyW122bJKCJm4YEMSGEeJssS/FXb0zzX58ZJVWs8pP3D/HPH9sszZWrqWbC0gk9ECBxrLFeHNfv1FwW6YP2jbrPWX0Z0ZU1eVuAWEckiAkhxDVKF6r8xjeO8eevnqIz7ONXPrydx3d0SXPljWRWYGnCrppN6GC2eBwWjjde4QTgcEPrBl01axvRS2wIWgb1uzadMiGpuLEkiAkhxHXy5qkkv/y3hzk6k+GRLXF+9SM7GGiTJrI1pRQUFnUoW14WxnVYWzoBtXLjWsMJ0V4dymKD0DJkr+39UKdMVCuuOwliQghxHZk1iz/63kk+/+wopqX4ufdu5KcfGsbrkikZ1h2rBunTkDoJyZPnr3OzK693+SDar8NapBciPfbS21j7YxLWxFWRICaEEKtgJl3k1752hKcOzTIcD/Iff2AH94+0r/VjiatRLer3btbD2RSkpyFzVi/ZmcaEtctcvkYoC3dDqEMvwQ4Ixe11BwTapa+aACSICSHEqvrWaIJ///eHmV4q8oN7evmlJ7cRD99ar2i5bdVMyCfsYHbmnLW9nZsHs3iBmw0ItOrmzmBchzN/q66oLS+Bc/Z9UZns9nqqmVDJQjkLpYxel7NQSuupU0op2PYRiG9Z1ceQICaEEKusVK3x/37rOF/YO4Hf7eSXntzGP76nXzrz3w6UgkoOcgnIz9vrhA5o+cTK48WkDgFc7O9eQ4cxfwy8Yb14QvZ2qLG94lgYPEFw+/RkuS4vuP26cufy6vV6D3dWDcySrlBWC+evKxc4Vi3qNzuUs1DO2Ms5oesSE7nW/dAfwM4fXtUfT4KYEELcIMcTOX75bw/xyuQS79zQyn/52E6G5VVJoplV02GsmGwshaWV+8UlO1TkGhWdck4Hvmrh6r/T4bbDmR3MnB4dzhwuPYDB4bzAvqvpmEMHTmUBSm+fu24+p5R+k0Jteak07VdWHq9VuHgwvQTDAe6gDqS+iB1MI40A64s2ts87HgF/i952rX71WoKYEELcQJal+Mrr0/znp45SNi1+/tFNfOqhYdxOx1o/mrgV1EwdyCo5Hc7KWb1tlnUTqVm2q0slva4vZV1FMst6JKlVA1XTa6um52hTTdsrjll6gIJhABdaO1YeMxw6wDk99tK0fbHjTrcOVm6/vQT0pL3uQGO/+ZzTc9MMmpAgJoQQayCRKfEr/+ctnj48y9auML/xQ7vY3d+y1o8lhLjBLhXE5J9nQgixSjoiPn7/x9/B//fJd5AsVPjB3/suv/a1IxQq5uVvFkLcFiSICSHEKvvAHV089y8e5sfuHeAPvjPJ+3/zRfaOza/1Ywkh1gEJYkIIcQNEfG7+0w/u5C8/fR8el4Of+MNX+Rdf2c9SvnL5m4UQtywJYkIIcQPdu6GVpz73IJ9770b+z4GzvO/ze/n6wZm1fiwhxBqRICaEEDeYz+3kX7x/C1/73AP0xfz8X19+k8/+2Rss5MqXv1kIcUuRICaEEGtka1eEv/nM/fybx7fw/JEEj31+L/9w4Cw342h2IcTbI0FMCCHWkMvp4LOPbOTrn3uAgbYgP/fn+/jMn77JfFaqY0LcDiSICSHEOrCpM8xXf+Y+fvGJrbwwmuCx39zL3+8/I9UxIW5xEsSEEGKdcDkd/MzDIzz1uQcYagvy83+xn0//yRsksqW1fjQhxCqRICaEEOvMxo4wX/3M/fzSk1vZOzbPY59/kb/bJ9UxIW5FEsSEEGIdcjoMPvXQCE/9/IOMxIP8wlf286k/kZGVQtxqJIgJIcQ6NhIP8Vc/cz+//OQ29o7N84HffJFn35pd68cSQlwnEsSEEGKdczoMfvqhYf7hZx+gK+rjU3/yBv/qrw6QKVXX+tGEENdIgpgQQtwktnSF+dvPvpufe+9G/ubN0zzxWy/x8sTCWj+WEOIaSBATQoibiMfl4F++fwt//Zn78bgcfPx/vcKvfe0IpWptrR9NCPE2SBATQoib0F0DMb7+uQf4J/cN8gffmeTDv/sdDp1Or/VjCSGukgQxIYS4SQU8Ln71ozv44396L9mSyQ/+3nf5nW+OY9astX40IcQVkiAmhBA3uYc2x/nGLzzEB3d18/nnxvihL3yPyYX8Wj+WEOIKSBATQohbQDTg5rd/dA//4+N7mFrI88HfeYm/fG1aJoEVYp2TICaEELeQD+3q4ZlfeJDdfS38m68e5LN/9iapQmWtH0sIcRESxIQQ4hbTHfXzpz/1Tn7xia08d2SOx2WaCyHWLQliQghxC3I6DH7m4RH+9rPvJuBx8okvvsKvP32Miikd+YVYTySICSHELWxnX5Svfe4BfvSefr6wd4If+v2XOTGfW+vHEkLYJIgJIcQtLuBx8V8+tosv/Pg7mE4W+ODvfIe/ePWUdOQXYh2QICaEELeJx3d08czPP8SegRZ+8W8O8Zk/fZNkXjryC7GWJIgJIcRtc8PqJgAAIABJREFUpCvq40//2Tv5t09s5ZvH5njit1/ilROLa/1YQty2JIgJIcRtxuEw+LTdkd/ndvBj/+v7/M43x6lZ0lQp/v/27jy6qvJ8+/h1Z4YkECFEIGHGQiCEEMKo0mqrdQJEFGUQARm0OFW7+rO1rW8HW6WtrQrKIIOKErXVKg5UWnFCQgiTMouAEkSIIGMIZHjeP5JaraABTvKcnPP9rJVlDmdzcq21l4sr9977eVDbKGIAEKYyUhvqpVvOVf+uzXX/wk26duZS7T5Q4jsWEFYoYgAQxhJio/TXq7M06cpMrfx4ny5+4G29sXG371hA2KCIAUCYMzMNyWmh+TefreSEWI2avUx/eHW9Stk8HKhxFDEAgCSpfUqiXrjpbA3r1VLT3tyiIdOWaPveYt+xgJBGEQMAfCEuOlK/H9RFU4Zla/OuQ7r0wbe1YM1O37GAkEURAwB8zaWZzfTyLeeqTXK8bpi7Qr/8xxqVlJb7jgWEHIoYAOC4Wjaur2dv6Kvx/drqibyPNOjhd7X1s8O+YwEhhSIGADihmKgI/fySdM0e1UM79x9R/4fe0avvc6kSCBSKGADgW53XMUUv33Ku2qck6MYnV+g389fpWBlPVQKniyIGAKiW1KR6emZCH43q21qzFm/VNdOX6JN9R3zHAuo0ihgAoNpioiL0/wZ01pRh2dr46UFd+uDbenNTke9YQJ1FEQMAnLRLM5vpxZvPUUpinEbNztf9CzexVyVwCihiAIBT0q5Jgv4x8WwNzk7Tg//+QNfNytdnh476jgXUKRQxAMApqxcTqT9d1VWTBmdq2ba9uvTBt1Wwba/vWECdQREDAJy2IT1a6Lkf9VVcdKSunp6nGW9tkXNcqgS+DUUMABAQnZs31Pybz9EF6WfqnlfWa+JTK3ToaJnvWEBQo4gBAAKmQVy0HhmRrZ9d3FEL1nyqQVMWa0vRId+xgKBFEQMABJSZacJ322nu9b205/AxDZy8WAvX7fIdCwhKFDEAQI3o2z5Z828+R22axGvc4wW6/7WNLHEB/A+KGACgxvxnNf6ruqfpwdc36/rHlml/canvWEDQoIgBAGpUXHSkJl2ZqXsGZWjx5s/Uf/I7Wr/zgO9YQFCgiAEAapyZaXivVnp6Qh8dLSvXoIcX64VVO3zHAryjiAEAak12yzM0/+ZzlJmapFtzV+k389eptLzCdyzAG69FzMyuMrO1ZlZhZjk+swAAakdKYpyeHNdLo89urVmLt2r4o0tVdJCtkRCefE/E1ki6QtJbnnMAAGpRdGSE7u7fWX+9OkvvFe7TwMnvaM2O/b5jAbXOaxFzzq13zm30mQEA4M/l3VL1txv6SpKunPqu5q/+xHMioHb5nogBAMJcRmpDvXDTOeqS2lA3z1upSQs2qIL1xhAmaryImdm/zGzNcb4GnuTnjDezAjMrKCoqqqm4AAAPmiTG6smxvTW0Zws9/MaHGvd4gQ6WsN4YQp855/+3DjN7Q9JPnHMF1Tk+JyfHFRRU61AAQB3inNMTeR/p1/PXqU1yvGaMzFGb5HjfsYDTYmbLnXPHfSiRS5MAgKBhZhrZp3XlPpWHjmrg5Hf01iaugiB0+V6+YpCZFUrqI+llM/unzzwAgODQp11jvXjTOWqeVE+jZufr0be3KBiu4ACB9q2XJs2staSJktpJ2itplaT5zrmPajrciXBpEgDCw+GjZbrjmdVasPZTXZGdqt8P6qK46EjfsYCTcrqXJl+QtEHSFEkXSOoq6S0zm2JmsYGLCQDAV8XHRunh4dn68Q++o+dW7NDV0/O060CJ71hAwFSniEU652Y65/4taa9zbpwqp2PbJE2vyXAAAEREmG79wVmaOqK7Pth1UAMnL2bxV4SM6hSxf5nZTVXfO0lyzpU55/6oynu7AACocRdlNNXfb+yryAjTVVOXaMGanb4jAaetOkXsdkkNzaxAUvOq9bxGmNkUSXtqNh4AAP+V3qyBnp/YVx2aJuqGuSv08BubuYkfddq3FjHnXIVz7h5J/SSNl9RUUndV7hN5cc3GAwDgq1IS45Q7vrcGdG2uSQs26o5nV+toWbnvWMApiarugc65YkkvVn0BAOBNXHSkHrgmS+2aJOgv/9qk7XuLNXVEdzVO4Bky1C0s6AoAqJPMKm/if2hoN71XuF+XP7xYH+w66DsWcFIoYgCAOq1/1+bKHd9bR45V6IqH39WbrMSPOoQiBgCo87q1PEMv3HS20hrV1+jZ+Xrs3W2+IwHVQhEDAISE1KR6+tsNfXR+xzN194tr9ct/rFFZeYXvWMA3oogBAEJGfGyUpl3bXRP6tdUTeR9p9JxlOlhS6jsWcEIUMQBASImMMP3sknTdN7iLlny4R1dNXaKd+4/4jgUcF0UMABCSru7RUrNH99COz4/o8imLtfYTtkVC8KGIAQBC1rlnNdGzN/ZRpJmGTF2iRRt3+44EfAVFDAAQ0jo2baDnJ56t1snxGvtYgZ5c+pHvSMAXKGIAgJB3ZoM4PTOhj777nSa66/k1+sOr61VRwR6V8I8iBgAIC/GxUZp+bXeN6N1S097copvnrVRJKXtUwq9q7zUJAEBdFxUZod8OzFCrRvG655X1+vRAiWaMzFGj+Bjf0RCmmIgBAMKKmWlcv7Z6eHi21uzYryseXqytnx32HQthiiIGAAhLl3RppqfG9daBkjINenixlm3b6zsSwhBFDAAQtrq3OkPP/6ivGtWP0fBHl+rV93f6joQwQxEDAIS1Vo3j9fcb+yqjeQP96KkVbBiOWkURAwCEvTPiY/Tk2N76ftWG4fct2CDnWN4CNY8iBgCApHoxkZo6IlvDerXUI298qDueWa3S8grfsRDiWL4CAIAqUZERuufyDDVrEKc/L9ykokNH9ciI7kqI5Z9L1AwmYgAAfImZ6ebvn6VJgzP17od7dM30Jdp9sMR3LIQoihgAAMcxpEcLPToyRx/uPqzBj7yrLUWHfEdCCKKIAQBwAud1TNG88b11+Gi5rpy6RCs//tx3JIQYihgAAN8gq0WS/n5jXyXERmnojDz9e/0u35EQQihiAAB8izbJlWuNnZWSqPFPLFdu/se+IyFEUMQAAKiGJomxyh3fW2e3T9adz72vKYs2s9YYThtFDACAaoqPjdLM63J0eVZz/fGfG/Xbl9arooIyhlPHwigAAJyE6MgI3T8kS0n1YzRr8VZ9XnxMk67MVHQksw2cPIoYAAAnKSLCdHf/TmocH6M/L9yk/UdKNWVYturFRPqOhjqG+g4AwCn4z8Kvv7s8Q4s27ta1M5dqf3Gp71ioYyhiAACchhG9W2ny0GytLtynq6cv0a4DrMKP6qOIAQBwmi7NbKbZo3rq473FGvzIu9r22WHfkVBHUMQAAAiAc85K1rxxvVV8rFxXTn1Xa3bs9x0JdQBFDACAAOnaIknPTOijmMgIDZ2ep7wte3xHQpCjiAEAEEDtUxL0txv76syGcRo5K1+vrf3UdyQEMYoYAAAB1jypnp6d0EfpzRrohrnL9bflhb4jIUhRxAAAqAFnxMfoqbG91Lddsn7y7Go9vmSb70gIQhQxAABqSHxslB69LkcXdDpTv3phraYs2uw7EoIMRQwAgBoUFx2ph4dna2DV/pT3LdjAZuH4AlscAQBQw6IjI/SXIVmKj43SI298qEMlZfr1gM6KiDDf0eAZRQwAgFoQEWG65/IMJcZGadpbW3T4aJkmXZmpKDYLD2sUMQAAaomZ6c6LOyoxLkp/em2TDh8r04NDuyk2is3CwxU1HACAWmRmuun8s3R3/07659pdGvtYgYqPlfmOBU8oYgAAeDD67DaadGWmFm/+TNfNyteBklLfkeABRQwAAE+G5LTQQ0OztWr7Pg2bkae9h4/5joRaRhEDAMCjSzObafrIHH2w65CunrZEuw6U+I6EWkQRAwDAs/M6pOixMT21c3+Jrpq6RIWfF/uOhFpCEQMAIAj0bttYT47tpX3Fx3T1tDx9tOew70ioBRQxAACCRNcWSZo3vreKj5VpyLQl+rDokO9IqGEUMQAAgkjn5g2VO76Pyiucrp6Wp027DvqOhBpEEQMAIMh0aJqo3PF9FGHSNdPztO6TA74joYZQxAAACELtUxL0zIQ+iouK0NAZeXqvcJ/vSKgBFDEAAIJU6+R4PT2hjxrUi9LwGUu1/KPPfUdCgFHEAAAIYi0a1dfT4/soOTFWI2cu1dIte3xHQgBRxAAACHLNk+rp6fG91bRhnK6bna/Fmz/zHQkBQhEDAKAOSGkQp6cn9FHrxvEaPWeZFm3c7TsSAoAiBgBAHZGcEKt543rrrJQETXh8uRau2+U7Ek4TRQwAgDrkjPgYPTW2t9KbN9CNc5frlfd3+o6E0+C1iJnZH81sg5m9Z2bPm1mSzzwAANQFDetHa+71PZXVIkk3z1tJGavDfE/EFkrKcM5lStok6Wee8wAAUCckxkVrzpie6lZVxl5+jzJWF3ktYs6515xzZVUv8ySl+cwDAEBdkhAbpTljeiq7ZZJuyaWM1UVRvgN8yRhJT5/oTTMbL2m8JLVs2fJr75eWlqqwsFAlJSU1FvB0xcXFKS0tTdHR0b6jAABCREJslGaP7qnRs/N1S+5KSdKlmc08p0J1mXOuZn+A2b8kNT3OW3c5516oOuYuSTmSrnDVCJSTk+MKCgq+8mdbt25VYmKiGjduLDMLQPLAcs5pz549OnjwoNq0aeM7DgAgxBw6WqbRs/O14uN9euCaLF2W2dx3JFQxs+XOuZzjvVfjEzHn3A++6X0zGyXpMknfr04JO5GSkhK1bt06KEuYJJmZGjdurKKiIt9RAAAh6MuTsVtzV0kSZawO8P3U5EWSfippgHOuOACfd/qhalCw5wMA1G0JsVGaM7rynrFbc1dp/upPfEfCt/D91ORkSYmSFprZKjOb6jnPaVuwYIE6dOig9u3b69577/UdBwAQZuKrylj3lmfotqcpY8HO6836zrn2Pn9+oJWXl2vixIlauHCh0tLS1KNHDw0YMECdOnXyHQ0AEEbiY6M0e3QPjZ69TLc9XXmZsn9XLlMGI98TsZCSn5+v9u3bq23btoqJidE111yjF154wXcsAEAY+k8Z697qDN2au5LJWJAKpuUrAubX89dq3ScHAvqZnZo30N39O3/jMTt27FCLFi2+eJ2WlqalS5cGNAcAANUVHxul2aN6aPScZbq1amkLJmPBhYkYAAAh7D9lLKd1I92au1IvvcdkLJiE5ETs2yZXNSU1NVXbt2//4nVhYaFSU1O9ZAEA4D++mIzNXqZbc1cpKiJCF2Ucb4lP1DYmYgHUo0cPffDBB9q6dauOHTum3NxcDRgwwHcsAAAUHxulWaN7qGtaQ908b4X+vX6X70gQRSygoqKiNHnyZP3whz9Uenq6hgwZos6d/UznAAD4X//ZmzK9WQPdOHeF3tzEIuO+UcQC7JJLLtGmTZv04Ycf6q677vIdBwCAr2gQF63Hx/RU+5QEjX+8QO9u/sx3pLBGEQMAIMwk1Y/R3LG91LpxvK5/rEBLt+zxHSlsUcQAAAhDjeIry1jzpDiNmbNMyz/63HeksEQRAwAgTDVJjNVT43qrSWKsRs3K1+rt+3xHCjsUMQAAwtiZDeL01LjeSoqP1rUzl2rNjv2+I4UVihgAAGGueVI9PTW2txLjKsvYhk8DuzsNTowiBgAA1KJRfT01rpdioyI1fMZSbd590HeksEARC6AxY8YoJSVFGRkZvqMAAHDSWjWO15PjeikiwjR0xlJtKTrkO1LIo4gF0KhRo7RgwQLfMQAAOGXtmiToqbG9VFHhNGzGUm3fW+w7UkijiAVQv3791KhRI98xAAA4LWedmai5Y3vpSGm5hj+6VLsOlPiOFLJCctNvvXqn9On7gf3Mpl2ki+8N7GcCABCk0ps10GNjemr4jDwNf3Spnh7fW40TYn3HCjlMxAAAwHFltUjSzFE9tH1vsUbOytf+I6W+I4Wc0JyIMbkCACAgerdtrGnXdte4xws0Zs4yPT6mp+JjQ7M++MBEDAAAfKPvdUjRQ0O7adX2fRr3eIFKSst9RwoZFLEAGjp0qPr06aONGzcqLS1NM2fO9B0JAICAuCijmf54Zabe/XCPbnpqhUrLK3xHCgnMFgNo3rx5viMAAFBjrshO0+Fj5frlP9box0+v0gPXdFNkhPmOVadRxAAAQLVd27uVio+W6Q+vblD9mEjde0WmIihjp4wiBgAATsqE77bT4WPlevDfH6h+TJTu7t9JZpSxU0ERAwAAJ+3HPzhLh4+WaeY7W5UQG6Wf/LCD70h1EkUMAACcNDPTLy5NV/GxMk1etFnxsVG68XvtfMeqcyhiAADglJiZfnd5FxUfK9d9CzYoPjZSI/u09h2rTqGIAQCAUxYZYfrTVV1VfKxcv3phrRrWi9bArFTfseoM1hELoO3bt+u8885Tp06d1LlzZz3wwAO+IwEAUOOiIyP00NBu6t22ke54ZrUWbdjtO1KdQRELoKioKP35z3/WunXrlJeXpylTpmjdunW+YwEAUOPioiM1Y2SO0ps10A1zl2vZtr2+I9UJFLEAatasmbKzsyVJiYmJSk9P144dOzynAgCgdiTGRWvO6B5KPaOexsxZprWf7PcdKeiF5D1i9+Xfpw17NwT0Mzs26qj/6/l/1T5+27ZtWrlypXr16hXQHAAABLPGCbGae30vXfnIu7puVr6evaGv2iTH+44VtJiI1YBDhw5p8ODB+utf/6oGDRr4jgMAQK1qnlRPT4ztpQonjXh0qT7dX+I7UtAKyYnYyUyuAq20tFSDBw/W8OHDdcUVV3jLAQCAT+2aJOix0T01dEaerp25VM9M6KMz4mN8xwo6TMQCyDmn66+/Xunp6br99tt9xwEAwKsuaQ01Y2SOPtpbrNFzlunw0TLfkYIORSyAFi9erCeeeEKvv/66srKylJWVpVdeecV3LAAAvOnTrrEmD+2m93fs14QnlutoWbnvSEElJC9N+nLOOefIOec7BgAAQeXCzk113+BM/eTZ1botd5UmD8tWZASbhEtMxAAAQC24snuafnFpul5d86l+/tz7DC6qMBEDAAC1Yuy5bbWvuFSTF21WUny0fnZxuu9I3lHEAABArbnjwu9o35FjmvbmFjWOj9H4fu18R/KKIgYAAGqNmek3AzL0eXGpfv/KBiUnxOqK7DTfsbyhiAEAgFoVEWG6f0hXfX74mH76t/fUKD5G3+uQ4juWF9ysDwAAal1sVKSmXdtd3zkzUTfOXaFV2/f5juQFRSyASkpK1LNnT3Xt2lWdO3fW3Xff7TsSAABBKzEuWnPG9FByYozGzFmmD4sO+Y5U6yhiARQbG6vXX39dq1ev1qpVq7RgwQLl5eX5jgUAQNBKSYzTE2N6ySSNnJmvXQfCa19KilgAmZkSEhIkVe45WVpaKjMWrAMA4Ju0To7X7NE99HnxMV03K18HSkp9R6o1IXmz/qe//72Ort8Q0M+MTe+opj//+bceV15eru7du2vz5s2aOHGievXqFdAcAACEosy0JE27trvGzFmmcY8V6LExPRUXHek7Vo1jIhZgkZGRWrVqlQoLC5Wfn681a9b4jgQAQJ1w7llN9Kerumrp1r26LXeVyitCf/X9kJyIVWdyVdOSkpJ03nnnacGCBcrIyPAdBwCAOmFgVqqKDh7V715er7tfXKPfDswI6dt8mIgFUFFRkfbtq3z89siRI1q4cKE6duzoORUAAHXL2HPbakK/tpqb97Eeen2z7zg1KiQnYr7s3LlT1113ncrLy1VRUaEhQ4bosssu8x0LAIA65/8u6qiiQ0d1/8JNSk6I1bBeLX1HqhEUsQDKzMzUypUrfccAAKDOi4gw3Tc4U3sPH9Mv/vG+GifE6Iedm/qOFXBcmgQAAEEpOjJCDw/PVpe0JN0yb6WWbdvrO1LAUcQAAEDQqh8Tpdmjeig1qZ7GPlagzbsP+o4UUBQxAAAQ1BrFx+ixMT0VHWm6btYy7T4YOqvvU8QAAEDQa9GovmaN6qG9h49pzJxlOny0zHekgKCIAQCAOiEzLUlThnfTuk8O6EdPrlBpeYXvSKeNIgYAAOqM8zueqd9d3kVvbirSL55fI+fq9ur7FLEAKy8vV7du3Vg/DACAGjKsV0vdfH57PV2wvc4v+Mo6YgH2wAMPKD09XQcOHPAdBQCAkHX7Bd/Rjn1HdP/CTWrWME5X5bTwHemUMBELoMLCQr388ssaO3as7ygAAIQ0M9O9V2TqnPbJ+tlz7+utTUW+I52SkJyIvf3MJn22/VBAPzO5RYLOHfKdbzzmtttu06RJk3TwYGitcQIAQDCKiYrQIyOyddXUJbpx7nI9c0MfdW7e0Hesk+J1ImZmvzWz98xslZm9ZmbNfeY5HS+99JJSUlLUvXt331EAAAgbiXHRmjO6pxrUi9bo2cu0Y98R35FOiu+J2B+dc7+UJDO7RdKvJN1wuh/6bZOrmrB48WK9+OKLeuWVV1RSUqIDBw5oxIgRmjt3bq1nAQAgnDRtGKc5o3vqyqnvatSsfP3thr5qWD/ad6xq8ToRc859+Y72eEl19hnUP/zhDyosLNS2bduUm5ur888/nxIGAEAt6dA0UdOu7a5tew5r/BMFOlpW7jtStXi/Wd/M7jGz7ZKGq3IidqLjxptZgZkVFBXVzRvyAABAzenbLll/uqqrlm7dq588+54qKoJ/vlPjRczM/mVma47zNVCSnHN3OedaSHpS0k0n+hzn3HTnXI5zLqdJkyY1Hfu0fO9739NLL73kOwYAAGFnYFaq7ry4o+av/kT3LdjgO863qvF7xJxzP6jmoU9KekXS3TUYBwAAhLgJ/dpqx+dHNO2tLWrZuL6G92rlO9IJ+X5q8qwvvRwoKfirKwAACGpmprv7d9J5HZroVy+s1Rsbd/uOdEK+7xG7t+oy5XuSLpR0q+c8AAAgBERFRmjysGx1ODNRE59coXWfBOeON76fmhzsnMtwzmU65/o753b4zAMAAEJHfGyUZo3qocS4aI2Zs0yf7i/xHelrfE/EAAAAakzThnGaNaqHDpaUasycZTp0tMx3pK+giAEAgJDWqXkDTR6erY27Durmp1aorLzCd6QvUMQCrHXr1urSpYuysrKUk5PjOw4AAJB0XocU/XpAZy3aWKRfz18n54JjjTHfWxyFpEWLFik5Odl3DAAA8CUjerfSx3uLNf2tLWrVuL7GntvWdySKGAAACB93XtRR2/cW655X1ivtjPq6KKOp1zwhWcQWzZmu3R9tCehnprRqq/NGjf/W48xMF154ocxMEyZM0Pjx3/53AABA7YiIMP3l6iztnJ6n255eqdyGfZTVIslbnpAsYj698847Sk1N1e7du3XBBReoY8eO6tevn+9YAACgSlx0pB69LkfjHi/wHSU0i1h1Jlc1JTU1VZKUkpKiQYMGKT8/nyIGAECQSU6I1XM39pWZec3BU5MBdPjwYR08ePCL71977TVlZGR4TgUAAI7HdwmTQnQi5suuXbs0aNAgSVJZWZmGDRumiy66yHMqAAAQrChiAdS2bVutXr3adwwAAFBHcGkSAADAE4oYAACAJxQxAAAAT0KqiAXLvlEnEuz5AABA7QqZIhYXF6c9e/YEbdlxzmnPnj2Ki4vzHQUAAASJkHlqMi0tTYWFhSoqKvId5YTi4uKUlpbmOwYAAAgSIVPEoqOj1aZNG98xAAAAqi1kLk0CAADUNRQxAAAATyhiAAAAnliwPmX4TcysSNJHNfxjkiV9VsM/AyeP8xJ8OCfBifMSfDgnwak2zksr51yT471RJ4tYbTCzAudcju8c+CrOS/DhnAQnzkvw4ZwEJ9/nhUuTAAAAnlDEAAAAPKGIndh03wFwXJyX4MM5CU6cl+DDOQlOXs8L94gBAAB4wkQMAADAE4rYcZjZRWa20cw2m9mdvvNAMrNZZrbbzNb4zoJKZtbCzBaZ2TozW2tmt/rOFO7MLM7M8s1sddU5+bXvTKhkZpFmttLMXvKdBZXMbJuZvW9mq8yswFsOLk1+lZlFStok6QJJhZKWSRrqnFvnNViYM7N+kg5Jetw5l+E7DyQzayapmXNuhZklSlou6XL+X/HHzExSvHPukJlFS3pH0q3OuTzP0cKemd0uKUdSA+fcZb7zoLKIScpxznld242J2Nf1lLTZObfFOXdMUq6kgZ4zhT3n3FuS9vrOgf9yzu10zq2o+v6gpPWSUv2mCm+u0qGql9FVX/y27ZmZpUm6VNKjvrMg+FDEvi5V0vYvvS4U/7gA38jMWkvqJmmp3ySougS2StJuSQudc5wT//4q6aeSKnwHwVc4Sa+Z2XIzG+8rBEUMwGkxswRJf5d0m3PugO884c45V+6cy5KUJqmnmXEp3yMzu0zSbufcct9Z8DXnOOeyJV0saWLVLTC1jiL2dTsktfjS67SqPwPwP6ruQ/q7pCedc8/5zoP/cs7tk7RI0kW+s4S5syUNqLofKVfS+WY2128kSJJzbkfVf3dLel6VtybVOorY1y2TdJaZtTGzGEnXSHrRcyYg6FTdGD5T0nrn3P2+80AysyZmllT1fT1VPnS0wW+q8Oac+5lzLs0511qV/5687pwb4TlW2DOz+KqHjGRm8ZIulOTlqXyK2P9wzpVJuknSP1V58/Ezzrm1flPBzOZJWiKpg5kVmtn1vjNBZ0u6VpW/4a+q+rrEd6gw10zSIjN7T5W/VC50zrFcAvB1Z0p6x8xWS8qX9LJzboGPICxfAQAA4AkTMQAAAE8oYgAAAJ5QxAAAADyhiAEAAHhCEQMAAPCEIgYAAOAJRQwAAMATihiAsGdm3zWz3WZWbmZbzewO35kAhIco3wEAIAg0lfSspLuq9mgEgFrBRAwApJGS/iVpv+8gAMILRQwApIckPSFpn5kN8x0GQPjg0iSAsGZmHSVNktRf0huODXgB1CI2/QYQ1swz9RL/AAAAl0lEQVTsL5JWO+fm+M4CIPxwaRJAuItT5c36AFDrmIgBCGtm1lbSDElnStopaaRzbqffVADCBRMxAGHNObfFOfd951yGpE8kZfvOBCB8UMQAQJKZXSYpXpXLWABAreDSJAAAgCdMxAAAADyhiAEAAHhCEQMAAPCEIgYAAOAJRQwAAMATihgAAIAnFDEAAABPKGIAAACe/H8t69p5FcrNAQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10,6))\n",
    "plt.xlabel(r'$\\xi$')\n",
    "plt.ylabel(r'$\\theta$')\n",
    "plt.title(\"Polytropes\")\n",
    "for n in [0,1,2,3,4,5]:\n",
    "    solve_polytrope(n)"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}