01_phugoid/.ipynb_checkpoints/01_03_PhugoidFullModel-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, G.F. Forsyth, I. Hawke. Partly based on content by David Ketcheson, also under CC-BY."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Full phugoid model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the third IPython Notebook of the series on the _phugoid model of glider flight_, our first learning module of the course [\"Practical Numerical Methods with Python\"](http://openedx.seas.gwu.edu/courses/GW/MAE6286/2014_fall/about). In the first notebook, we described the physics of the trajectories known as phugoids obtained from an exchange of potential and kinetic energy in an idealized motion with no drag. We gave you a neat little code to play with and plot various phugoid curves.\n",
    "\n",
    "In the second notebook, we looked at the equation representing small perturbations on the straight-line phugoid, resulting in simple harmonic motion. This is a second-order ordinary differential equation, and we solved it numerically using **Euler's method**: the simplest numerical method of all. We learned about convergence and calculated the error of the numerical solution, comparing with an analytical solution. That is a good foundation!\n",
    "\n",
    "Now, let's go back to the dynamical model, and take away the idealization of no-drag. Let's remind ourselves of the forces affecting an aircraft, considering now that it may be accelerating, with an instantaneous upward trajectory. We use the designation $\\theta$ for the angle, and consider it positive upwards."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Image](./figures/glider_forces-lesson3.png)\n",
    "#### Figure 1. Forces with a positive trajectory angle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Figure 1, $L$ is the lift, $W$ is the weight, $D$ is the drag, and $\\theta$ the positive angle of the trajectory, instantaneously. \n",
    "\n",
    "In lesson 1, we wrote the force balance in the directions perpendicular and parallel to the trajectory for a glider in _equilibrium_. What if the forces are _not_ in balance? Well, there will be acceleration terms in the equations of motion, and we would have in that case:\n",
    "\n",
    "\\begin{align}\n",
    "m \\frac{dv}{dt} & = - W \\sin\\theta - D \\\\\n",
    "m v \\, \\frac{d\\theta}{dt} & = - W \\cos\\theta + L\n",
    "\\end{align}\n",
    "\n",
    "We can use a few little tricks to make these equations more pleasing. First, use primes to denote the time derivatives and divide through by the weight:\n",
    "\n",
    "\\begin{align}\n",
    " \\frac{v'}{g} & = - \\sin\\theta - D/W \\\\\n",
    "\\frac{v}{g} \\, \\theta' & = - \\cos\\theta + L/W\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Recall, from our first lesson, that the ratio of lift to weight is known from the trim conditions—$L/W=v^2/v_t^2$— and also from the definitions of lift and drag, \n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "L &=& C_L S \\times \\frac{1}{2} \\rho v^2 \\\\\n",
    "D &=& C_D S \\times \\frac{1}{2} \\rho v^2\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "we see that $L/D=C_L/C_D$. The system of equations can be re-written:\n",
    "\n",
    "\\begin{align}\n",
    " v' & = - g\\, \\sin\\theta - \\frac{C_D}{C_L} \\frac{g}{v_t^2} v^2 \\\\\n",
    " \\theta' & = - \\frac{g}{v}\\,\\cos\\theta + \\frac{g}{v_t^2}\\, v\n",
    "\\end{align}\n",
    "\n",
    "It is very interesting that the first equation has the factor $C_D/C_L$, which is the inverse of a measure of the aerodynamic efficiency of the aircraft. It turns out, this is the term that contributes damping to the phugoid model: if drag is zero, there is no damping. Drag is never zero in real life, but as engineers design more aerodynam-ically efficient aircraft, they make the phugoid mode more weakly damped. At altitude, this is nothing but a slight bother, but vertical oscillations are unsafe during final approach to land, so this is something to watch out for!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The initial value problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to visualize the flight trajectories predicted by this model, we are going to need to integrate the spatial coordinates, which depend on both the forward velocity (tangential to the trajectory) and the trajectory angle. The position of the glider on a vertical plane will be designated by coordinates $(x, y)$ with respect to an inertial frame of reference, and are obtained from:\n",
    "\n",
    "\\begin{align}\n",
    "x'(t) & = v \\cos(\\theta) \\\\\n",
    "y'(t) & = v \\sin(\\theta).\n",
    "\\end{align}\n",
    "\n",
    "Augmenting our original two differential equations by the two equations above, we have a system of four first-order differential equations to solve. We will use a time-stepping approach, like in the previous lesson. To do so, we do need *initial values* for every unknown:\n",
    "\n",
    "$$\n",
    "v(0) = v_0 \\quad \\text{and} \\quad \\theta(0) = \\theta_0\\\\\n",
    "x(0) = x_0 \\quad \\text{and} \\quad y(0) = y_0\n",
    "$$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve with Euler's method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know how to apply Euler's method from the previous lesson. We replace each of the time derivatives by an approximation of the form:\n",
    "\n",
    "$$v'(t) \\approx \\frac{v^{n+1} - v^n}{\\Delta t},$$\n",
    "\n",
    "where we are now using a superscript $n$ to indicate the $n$-th value in the time iterations. The first differential equation, for example, gives:\n",
    "\n",
    "$$\\frac{v^{n+1} - v^n}{\\Delta t} = - g\\, \\sin\\theta^n - \\frac{C_D}{C_L} \\frac{g}{v_t^2} (v^n)^2$$\n",
    "\n",
    "Alright, we know where this is going. At each time iteration $t^n$, we want to evaluate all the known data of our system to obtain the state at $t^{n+1}$—the next time step. We say that we are _stepping in time_ or _time marching_.\n",
    "\n",
    "The full system of equations discretized with Euler's method is:\n",
    "\n",
    "\\begin{align}\n",
    "v^{n+1} & = v^n + \\Delta t \\left(- g\\, \\sin\\theta^n - \\frac{C_D}{C_L} \\frac{g}{v_t^2} (v^n)^2 \\right) \\\\\n",
    "\\theta^{n+1} & = \\theta^n + \\Delta t \\left(- \\frac{g}{v^n}\\,\\cos\\theta^n + \\frac{g}{v_t^2}\\, v^n \\right) \\\\\n",
    "x^{n+1} & = x^n + \\Delta t \\, v^n \\cos\\theta^n \\\\\n",
    "y^{n+1} & = y^n + \\Delta t \\, v^n \\sin\\theta^n.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we've learned before, the system of differential equations can also be written as a vector equation:\n",
    "\n",
    "$$u'(t) = f(u)$$\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{align}\n",
    "u & = \\begin{pmatrix} v \\\\ \\theta \\\\ x \\\\ y \\end{pmatrix} & f(u) & = \\begin{pmatrix} - g\\, \\sin\\theta - \\frac{C_D}{C_L} \\frac{g}{v_t^2} v^2 \\\\ - \\frac{g}{v}\\,\\cos\\theta + \\frac{g}{v_t^2}\\, v \\\\ v\\cos\\theta \\\\ v\\sin\\theta \\end{pmatrix}.\n",
    "\\end{align}\n",
    "\n",
    "It's a bit tricky to code the solution using a NumPy array holding all your independent variables. But if you do, a function for the Euler step can be written that takes any number of simultaneous equations. It simply steps in time using the same line of code:\n",
    "\n",
    "```Python\n",
    "def euler_step(u, f, dt):\n",
    "    return u + dt * f(u)\n",
    "```\n",
    "\n",
    "This function can take a NumPy array `u` with any number of components. All we need to do is create an appropriate function `f(u)` describing our system of differential equations. Notice how we are passing a _function_ as part of the arguments list to `euler_step()`. Neat!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### And solve!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, we start by loading the modules and libraries that we need for this problem. We'll need a few transcendental functions, including the $\\log$ for a convergence study later on. And remember: the line `%matplotlib inline` is a magic function that tells Matplotlib to give us the plots in the notebook (the default behavior of Matplotlib is to open a pop-up window). In addition, we are importing the module `rcParams` to define notebook-wide plotting parameters: font family and size. Here we go!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from math import sin, cos, log, ceil\n",
    "import numpy\n",
    "from matplotlib import pyplot\n",
    "%matplotlib inline\n",
    "from matplotlib import rcParams\n",
    "rcParams['font.family'] = 'serif'\n",
    "rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to set things up to start our numerical solution: the parameter \n",
    "values and the _initial values_. You know what the acceleration of gravity is: 9.8 m/s$^2$, but what are good values for $C_D/C_L$, the inverse of the aerodynamic efficiency? Some possible values are given on a table in the Wikipedia entry for [lift-to-drag ratio](http://en.wikipedia.org/wiki/Lift-to-drag_ratio): a modern sailplane can have $L/D$ of 40 to 60, depending on span (and, in case you're interested, a flying squirrel has $L/D$ close to 2).\n",
    "\n",
    "For the _trim velocity_, the speed range for typical sailplanes is between 65 and 280 km/hr, according to Wikipedia (it must be right!). Let's convert that to meters per second: 18 to 78 m/s. We'll pick a value somewhere in the middle of this range.\n",
    "\n",
    "Here's a possible set of parameters for the simulation, but be sure to come back and change some of these, and see what happens!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# model parameters:\n",
    "g = 9.8      # gravity in m s^{-2}\n",
    "v_t = 30.0   # trim velocity in m s^{-1}   \n",
    "C_D = 1/40  # drag coefficient --- or D/L if C_L=1\n",
    "C_L = 1   # for convenience, use C_L = 1\n",
    "\n",
    "### set initial conditions ###\n",
    "v0 = v_t     # start at the trim velocity (or add a delta)\n",
    "theta0 = 0 # initial angle of trajectory\n",
    "x0 = 0     # horizotal position is arbitrary\n",
    "y0 = 1000  # initial altitude"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll define a function `f()` to match the right-hand side of Equation (15), the full differential system in vector form. This function assumes that we have available the parameters defined above. If you re-execute the cell above with different parameter values, you can just run the solution without re-executing the function definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(u):\n",
    "    \"\"\"Returns the right-hand side of the phugoid system of equations.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    u : array of float\n",
    "        array containing the solution at time n.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    dudt : array of float\n",
    "        array containing the RHS given u.\n",
    "    \"\"\"\n",
    "    \n",
    "    v = u[0]\n",
    "    theta = u[1]\n",
    "    x = u[2]\n",
    "    y = u[3]\n",
    "    return numpy.array([-g*sin(theta) - C_D/C_L*g/v_t**2*v**2,\n",
    "                      -g*cos(theta)/v + g/v_t**2*v,\n",
    "                      v*cos(theta),\n",
    "                      v*sin(theta)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the code defining function `f(u)` with the differential equations, and convince yourself that it's right!\n",
    "\n",
    "\\begin{align}\n",
    "u & = \\begin{pmatrix} v \\\\ \\theta \\\\ x \\\\ y \\end{pmatrix} & f(u) & = \\begin{pmatrix} - g\\, \\sin\\theta - \\frac{C_D}{C_L} \\frac{g}{v_t^2} v^2 \\\\ - \\frac{g}{v}\\,\\cos\\theta + \\frac{g}{v_t^2}\\, v \\\\ v\\cos\\theta \\\\ v\\sin\\theta \\end{pmatrix} \\nonumber\n",
    "\\end{align}\n",
    "\n",
    "Now, Euler's method is implemented in a simple function `euler_step()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def euler_step(u, f, dt):\n",
    "    \"\"\"Returns the solution at the next time-step using Euler's method.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    u : array of float\n",
    "        solution at the previous time-step.\n",
    "    f : function\n",
    "        function to compute the right hand-side of the system of equations.\n",
    "    dt : float\n",
    "        time-increment.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    u_n_plus_1 : array of float\n",
    "        approximate solution at the next time step.\n",
    "    \"\"\"\n",
    "    \n",
    "    return u + dt * f(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After defining a final time for the solution, and the time step $\\Delta t$, we can construct the grid in time using the NumPy function [`linspace()`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html). Make sure you study the decisions we made here to build the time grid: why do we add 1 to the definition of `N`, for example?\n",
    "\n",
    "Look at the code below, and make sure you understand the following aspects of it.\n",
    "\n",
    "* The NumPy array `u` contains the solution at every time-step,  consisting of the velocity, angle and location of the glider. \n",
    "* The first element of the array `u` is set to contain the initial conditions. \n",
    "* In the `for`-loop, the function `euler_step()` is called to get the solution at time-step $n+1$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T = 100                          # final time\n",
    "dt = 0.1                           # time increment\n",
    "N = int(T/dt) + 1                  # number of time-steps\n",
    "\n",
    "# initialize the array containing the solution for each time-step\n",
    "u = numpy.empty((N, 4))\n",
    "u[0] = numpy.array([v0, theta0, x0, y0])# fill 1st element with initial values\n",
    "\n",
    "# time loop - Euler method\n",
    "for n in range(N-1):\n",
    "    \n",
    "    u[n+1] = euler_step(u[n], f, dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the trajectory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to plot the path of the glider, we need the location (`x`, `y`) with respect to time. That information is already contained in our NumPy array containing the solution; we just need to pluck it out. \n",
    "\n",
    "Make sure you understand the indices to `u`, below, and the use of the colon notation. If any of it is confusing, read the Python documentation on [Indexing](http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get the glider's position with respect to the time\n",
    "x = u[:,2]\n",
    "y = u[:,3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to plot the path of the glider and get the distance travelled!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiAAAAGcCAYAAADgaRuPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4FFXWwOHfgYQlQRZZRNwiiriMgKAygBJQRBwUFdlk\nREE2EVDgU8CREUcdFXEQBwRllEVkWEWRZRDZXBBZBAQFxWEIssm+Y4CQ8/1R1dA2HVIJSaqSnPd5\n+kmq+lbV6ZNKcrrurduiqhhjjDHG5KQCfgdgjDHGmPzHChBjjDHG5DgrQIwxxhiT46wAMcYYY0yO\nswLEGGOMMTnOChBjjDHG5DgrQPIpEWkqIp+IyGYR2SEi+0VkhYi8LSJNRKSQ266yiGwXkd9EJDVi\nH4VFZJ2IjPJwvGkiskdEUkXkkex6XUEgIgtFZIHfcQSZiDQQkSUi8quIbBORY+7jZGTuRKSKe+48\nmsljxbnn8CH3/Ls0a15FhuPIl+eFiFwlIt9E/v1Io+1jIvKD+zfpvyLyFxGJ+n9KRBJEZLJ7Du0Q\nkf+ISJUMxhYjIs+6x9rhHrvzWdpXFZHZbttfRWSSX+dTXmAFSD4jIvEi8gkwGBgNVFTVC4BywIvA\nncDHQEcAVf1JVS8EJgCRk8bEAKWA0ukdV1XvBZqGFs/9lWSMiDyfg/98SuMhJ+cqh19TlhGR8jjn\n2ApVLQ/8AdgNVAJ+4czzIw44DyiZmeOp6lH3HH49g3GmeimuM7BNjpwXQSIi3YEvgctJ5/deRF4A\n/gF0d/8mPQD0AM7Ip4hcBCx2FysCFwMbgEUi8ocMhDgSeAJo6h7zSeANEflblGNWARYBPwEXAVfi\n/A38RkQqZOCYxmUFSP4zCqgH1FPVqap6AkBVj6vqR0Aj4Dhn/rEQ93GKqh4BLgXu9XhsSb9JnlAN\nuMHvIAKsJk5R8TGAqu4FqgBbiXKOqOo3QAlVHXSOx83M+ZeZYjmtbfLVeSEirYFmwB+Bdem0vQr4\nCzBEVecDqOp3OG+K2ohIvYhNXgGKAx3dAvME0As4BAz1GN9twENAf1Vd7R5zLjAc+IuIXBmxyVBg\nL9BLVVNU9TDQCTgfeNnLMc3vWQGSj4jI7Th/EEao6v+itVHVn4AZQIqXfbqFS26aTjfbiyBVPamq\nJ7P7OGFyW2FXyv2aHFqhqnvOljNV/S3bo4ouM7mNuo0P54Xf5qpqoqpu9NC2Hc7/o6kR6z90v3YI\nrRCR84CWwOequj+0XlWP4/ztqhuleIimI06xGO2YBYG2Yce8CrgFmBH+M1TV3cAXQEsRKebhmCaM\nFSD5S6gPffrZGqlqM1UdcbY2IlI9rbEh7vM1RORzETkoIv8TkTeBomfZX2sRWe729e8RkU9F5Jaw\n5383FkVELhKRqSKyxV1Os29dRH4A/s9dXObuZ5uIXCEiz7nLqSIySkTud+M4EBqvIiIFRORJEflC\nnDEz+0RktYg8FnGchHRycpPbTx16jd9GG9cgjsdFZJXb17xFRL4UkT4iUvwsr2m7iFSM+BnMCOsj\nXyYiD0Yca7acHpvTzu0PXyfOeIxU97nj7vfbxR2/E/bzTxGR3SKSmFb+o/wsBruLU919zDpL+25u\nmzPGhrjPtxaRte7P5CcR6SciL7jx7hSRd6Ls9nwR+UBEtrr7fkdE4sL22U5EtruLLcNyO+wscaa5\njYiUcX8GvxuDEuWcvlhEZonILhFJEpEn3Hat3HNhv4jME5HLoxy/kIj0F5Gf3Z/ZdhF5X3zsnlPV\nnRloXhenGFgdsY9tOFcd6oat/iMQG9nWFVrn5XysC+yKEufqsOdDEiOei2xfGOfKnskIVbVHPnkA\nG4GTQJlMbDsaSI2yfhRwMmJdZeAgMBfn3a7gXHlZBaQCD0e074tzxaWN27YoMASnK+iOKMdLxbl8\nf7W77u/A/HTi7+9ud2kaz6cCP+Bcfi0KxAM/Ag8DxdznewHitm8JnACe9piTejjv+IcCRXDeYbVx\nX/fzEW3fBo4Cf3KXC7mvMRW418trijheYXdd6Hh/iWib6O5nFdAe543J1W7+LwXecp+/MmK70sAe\nIC6D59Ij7v7qRnkuKdrP0j1350es+7O7nwE4ffFFcC6F73DP80sj2j/vtp8DVHXX3eq+zsFpnBMj\nM/ja0twmrZ9X2Dn9AVAe53fgDXfdYOAF93y5APgvsChi+wLALGA7cLO7rgLOO/PtQHkPcT/ltvXy\nmJKJvx8LifidiHh+F7AvjedWuz/PIu5yVzc3T0Zp29R97pV04ol3261M4/mDwPaw5YFE/P6FPdfL\nfa5zRvOS3x92BSR/Ke9+3ZeF+zxjbAjwEs4/7Z6quk8dU4BlZ2wskuC2H6+qY922vwE9cf6RRA4c\nDB3rA1X90f3+HdLv9/VyKb0c0ENVf1NnfEsPYCnOP+0ZqjpI3b84qjoRmILzhzvasU4dT5xR/COA\nbTh/NJPVuRw/Fudy7zMicoHb9hacfuUhqjrLPdZxVX0WWM/vxxdEfU1hxzuE8zM45u5nLPAf4G/h\nV0rC9rNFVd9T1VQ3t11wzpWR7vPtIw71MPCRqh6NFsdZnHO3hojE4PxT2Ak8o06ffLKq/gU4kM6+\nZqozvgBV/RJnMON9mYgpo9J63aH1/1bVX91zLHSVqAVOgXpSVXcA44Ba4gzkDWmNM3arn6ouhVNX\nDjrjFC3PpBeYqr6uqhd6fDTL8CtPXwmcojuao2Ftwr9Gax9al96A5bPtI7Q+fB9ZcUwTwQqQ/OmM\nP4QicnfYZeODIvJ9pnYsUhC4C+fdw5qIpxdG2eQBnPPws/CVqpqC887nehG5MMp2X4S1/UVVI/tx\nM2NZ6J+1u9/Zqvqj+4+tSZT2PwNlRaRsOvu9AWfE/AI9cwzAcpzLybe5yy3dr7Oj7KcFp0f+ezne\nZ+oOMg7zCc676QeibPdF+IJbjBxS1W+BNcDD8vtbIttyujjJaTVwCup5qhrZ3fU5Zy9yFkUsb8W5\nYuC3pWHfh7pzVkS8vq3u14vD1oXOmcjfoXXAbzjFiTGBE+N3ACZH/YpzSf18nHeOp6jqDOBCABFJ\nwrlLITPKutv+GOW5X6OsCw0We0NEBkQ8Vxg4jPPPYXv4E5qx/mWv0tynOKPwewHX4Vy+VZyrPMpZ\nxra4Qq+xhYj8KeK5GJzXeEFE222RO1F3pL4HoX1sj/JcaN0VUZ47W05HAYNwisuZInITTtfO1x5j\nymqhcRDRXmO08yzc7ojl4wTgb6E6AxpD3x8XEYgeK/z+9zP0817ibhMuhdODfoPsAGn/zYkLaxP+\nNVr7yLZnO15a+witD99HVhzTRLArIPnLIpx3htfnwLEyepm9XZRLveeranH3HXhOiDpRkog0AeYD\n+4GaqlpeT88rkZHX+a8or7Gs+xoHR7QtlKlXcG7ONlHUBzhjXkKDZh8F3sv2iNKXme6cdCfEChCv\nsSpQOcr5VUJVy6W3sYg8FXYFNL3HlHN7SVGtB4qLSJEoz1UAtqlq6K6pn9yv0a6Mhq5k/Xy2g7ld\nrNuj7UOcu1mKRexj/bke05zJCpD8JTShz/3ZeIxdOH2i5aM8F21d6I/JxZFPiEgJcWbMLJiF8WXG\nozh/4HuFv0vF+z+/NF8jgIjUF5Ey7mKaf+hEpLSIlIhcH0VoH9G6FTL1x9J93dOBu0XkMpxBxe9n\nZB9ZLHQbudfzLC9bj3MuRvsdShCRm9PbQQDGgIS6zapGxF8B54rtwrDV3+BcCfpdW1doJtSFUZ6L\ndsxyofFX6ezjc/drWsdMduMyGeBrASIiF4pzG2BuekeSa6nqPJxBj+1F5NrM7CK99e4Yh5nAhXLm\ntMj1omw7FecycdMoz3UGBkUZN5GZeUeOuF9jAESkrYg097jtMaIXG55ucVTVVThFyB3ht3u6cdwM\nfMrp38VJ7td7ItoVwLlLJfyPf1qvaRVOgXGHuFPqh2mCk+/MjJkZiTNeZQLOnRg7MrGPrPItzjvY\n3xWo4vRBhG7pzIho7Y9yOrflRGSoiMSns59o23jpzjyXuXQmuF+j/Q4Nw7lbKOhG4dzpEvnmKHS+\nn7raps4EYJOARBE51b3knuv3AF+o6obwnbg/i8jfhXfdr9GOmYJz51/omOuBr3AK8FPdde74r0Rg\nkntVxWSAbwWIiDTF6RJI4Cy/fCJSzP0l/lGcefo/jfbPU0RiReRFceYwWCMii0SkThr77OHu6ztx\n5mLwOpNnXtAWWADME5EHQ5c83fzVF2ea9ks5/a49XHqj+EP+ijOuYZCInC/OPBpNgTsi26vqLzij\n9G8TZ66NWHE0xpkZsXcG4jibte7X60SkMPA0Z75TTmu/oaJgYOifiTuWo3UG4umIM6blndBVDPc8\nHgW8FhrToqqLcG7DfVREGrntiuLclnkAGJ/ea3LvouiIM1ZlkIgUcXPaBmdAYn+NPjlUeq9jNs7Y\nlJpEDD4VZ66SVR7y4OV46a53i9KnccYcveKeN0Vw7qg6eZZ9pLX/aOvWAVe7Rc2dQEsP/2Qit2kV\n5S4hr8f3tF5Vx+Ncnfo/cSYbDM0L8hec2VcHphNzTknzZ6KqP+PMbtpNnBlKEZFqQD/gfVVdGLFJ\nX5wu0X+J8/EShXB+R+KBbr87qEhtnPN2ecQx5+N0LT4ferMkIg2Ax4CXI4sYnNt/S+L8TsW6XTUj\ncMbppHunkYlCfbr/l9PFx2iizC8R1u4/OKPzQ/eAv4AzWK5CRLu3cQY+lnaX2+O8Q6wa0a4vTjfB\n5e5yA5zLeY38yoVP+W+BM3fANpx3kruA79w83h7WrrL7/FGcP+zbcX7Zboi2Pmy76jiXMA8Cm3De\nbdyH06e9H1gTEU8znDs8duN8HshcoH7Y8+fjDC4MP95Z5/6I2L/gzPHxq/uax+MMHu3i7uuku+/t\nRLmfH2emxjXu6/keGItz+2+qez4+H9Y26jmN889gmpvrrcBKoFMa8T6OcxfQDjcfo4ALPbymImHP\n18CZGXKH+1gKPBixj7E4c3mcdH8u24Hrz5LHl91jFQhbVxjnKtF4Dz+HH9zjnHSPux3nylhb9/sT\n7r62ueu7RqzfDiSG7e9Bd5/7cAqyLpye7+PisHYbcG5LPun+vPrhjLOJPKfah21TC+dq0k6cwqKx\nh9d3xjZAGfc44ccfhjM4NPL4D7qP8HNyG86bxWkRuRsTdtwYnN/LH939b8T553qFj39jCuJ0Tfzm\nxnzS/f43wv7GRGzzmPtz3OH+zP4Sfq5FtL0MmBx2fv8n2rkLXOs+Py2NGJ91j7XDPXaa83ngdMHM\nDjvmJNKYW8ge6T9CkyrlOBERVVURGY0zMdUZV2NE5A6cy9O3qVsBi0gszi/teFXt5q6rjHPitFfV\n0WHbfw8kqerd7nJJnD/8A1X1+bB2M4AEVc3IhxgZE5V7Tj+kqr7fWZHVxPnAsMKq2idsXUuc4qeh\nOp+l4SsRGY4zl8p5mvE5SowxOcS3Lhj1Vvk8gHN14quw7U7gXD0Jn8fgfpx3g5FTNS8AGob1wTbC\nedcbrd21biFjTIaIMx34TWGrSgOb/YonK4nzceWhbqcCOB/eNTLs+Uo4l74H53TxIc70552iPPVH\nnPkzrPgwJsCCfhdMFZzbryI/GC0JuCDszoEqOJf3folotxHn0uS1Ye1C6yPbQc7cnmrynpuAnu5Y\nlytwBqVNSGeb3OIWTn920CPAT+p8YGHIHpyp3XvleGTO3Cn9w/rvi4pIf5y5Wv7iQzzGmAwI+iXi\nMjj9ppEOul9L44wZKAMcjXJVJbxdaH9E2WdkO2MyYjbO2KTQmJhhOJ/7kRfsAyqLyC6cfvKW4U+q\n6l7C7hbIYTOAq4AP3QGBRXHujmmoZw5aNMYETNALEGMCT1U/xRmrlOeo85kpUecv8Zuqfs/pidGM\nMblM0AuQ3USfVKi4+3VPWLv40MDWdNoBnMfvP5Atst0pIuLPKF1jjDHGJ6qamekOMiToY0BWAxeF\nT/ziuhz4VU/PSvkdzmu5JEq7E5yeL+E792tClHah453hyJEj7Nu3jx07drBx40YWLVrEBx98QM+e\nPalZsyYS9vkLV199NVOnTiU1NdX3W5xy8tG/f3/fY8gND8uT5cryZLkK+iOnBKUASesVf4gz8+Kp\nCcXcCWfquM+FfOTuo37E9vWBOXp6NPxsnD76aO1+UGe2uzPExcVRsmRJypUrR0JCArVr1+bPf/4z\ngwYN4ptvvmHLli0MHDiQSy65hB9//JGmTZvSsGFDNm3a5OGl5w1JSUl+h5ArWJ68s1x5Y3nyznIV\nLEEpQKJe6lHVz3D61l90Z4MEZ9KYEzgTIoXarceZke4ZESkNzi16OFc2ng1rdwB4EegqIpe77RoA\nDYGnMht8hQoVeOqpp9iwYQNDhw6lVKlSzJ07l6pVqzJx4sTM7tYYY4zJs/ycin2YiGzEmc9DRWSj\niPzPnWgsXDOcrpFVIrIW5x7/eqoa+THc3XFmxVskImuADjij4X/XraKqA4C/AzNE5DvgNaCZOgMJ\nz0lsbCxdu3Zl/fr13HfffRw4cIBWrVrx3HPP5ehlLT+0bdvW7xByBcuTd5YrbyxP3lmugsW3mVBz\nizPHtXqjqgwdOpQePXqQmprKgw8+yJgxY4iNjayvjDHGmOAQEdQGoeZeIkL37t2ZMWMGxYoVY/z4\n8bRq1YoTJ074HVq2WLhwod8h5AqWJ+8sV95YnryzXAWLFSDZ7K677mLBggWUKFGCqVOn0rp1a1JS\nIid2NcYYY/IX64JJR2a7YCItW7aMBg0acPDgQTp06MCIESN+d/uuMcYYEwTWBZPH3HTTTcyePZsi\nRYrw7rvv8sorr/gdkjHGGOMbK0ByUK1atRg3bhwiwrPPPpunbtG1vlVvLE/eWa68sTx5Z7kKFitA\ncljTpk0ZNGgQAO3atWP16qiTrxpjjDF5mo0BSUdWjQEJp6q0a9eOMWPGcMUVV7Bs2TJKlSqVpccw\nxhhjMiOnxoBYAZKO7ChAAH777Tfq1KnDypUr+dOf/sT06dMpUMAuSBljjPGXDULN44oWLcrUqVM5\n//zzmTVrFgMHDvQ7pHNifaveWJ68s1x5Y3nyznIVLFaA+CghIYGxY8cC0K9fP5YsWeJzRMYYY0zO\nsC6YdGRXF0y4Xr168cYbb5CQkMCqVasoUaJEth7PGGOMSYuNAQmInChAjh07Ru3atVmxYgUtW7Zk\n/PjxNkmZMcYYX9gYkHykcOHCTJgwgWLFijFx4kRGjhzpd0gZZn2r3lievLNceWN58s5yFSxWgARE\npUqVGDZsGADdu3dn3bp1PkdkjDHGZB/rgklHTnTBhHv44YcZO3YsVapUYcmSJRQpUiTHjm2MMcZY\nF0w+9dZbb3HllVeyevVqnn76ab/DMcYYY7KFFSABc9555zFhwgRiY2MZOnQon3zyid8heWJ9q95Y\nnryzXHljefLOchUsVoAEUI0aNXj11VcB5/NitmzZ4nNExhhjTNayMSDpyOkxICGpqancfffd/Oc/\n/yExMZF58+ZRsGDBHI/DGGNM/mJjQPK5AgUKMHr0aMqXL8/nn3/Oyy+/7HdIxhhjTJaxAiTAypUr\nx9ixYxERnn/+eb766iu/Q0qT9a16Y3nyznLljeXJO8tVsFgBEnANGjSgd+/epKam8uCDD7Jz506/\nQzLGGGPOmY0BSYdfY0DCnThxgsTERBYvXkxiYiKfffYZsbGxvsZkjDEmb7IxIOaU2NhYpkyZcmo8\niM0PYowxJrezAiSXqFChAh9++CGxsbG8+eabjBkzxu+Qfsf6Vr2xPHlnufLG8uSd5SpYrADJRWrX\nrs3QoUMB6NChA3PmzPE5ImOMMSZzbAxIOoIwBiRS7969GThwIPHx8SxcuJAbb7zR75CMMcbkETk1\nBsQKkHQEsQBJTU3lkUce4YMPPqBs2bJ8+eWXVK5c2e+wjDHG5AE2CNWkqUCBArz33ns0bNiQXbt2\nUbduXdasWeNrTNa36o3lyTvLlTeWJ+8sV8ES+AJEROqKyAIR+VlEkkRkpohcE6VdFRGZLiL/E5EN\nIvK5iNSO0i5WRF4UkXUiskZEFolInZx5NVmnUKFCTJ06lTvuuIOdO3dSr149li1b5ndYxhhjjCeB\n7oIRkURgLvCMqr7urnsLaAHcoKpb3HVXA8uAUar6hLuuN/A8cIuqrgjb59tAPaCOqu4RkfbAP4Ha\nqvpdlBgC1wUTLjk5mRYtWjB9+nSKFCnCu+++y5///Ge/wzLGGJNL2RgQQES+ABJU9dKwdecBO4Ax\nqtrFXfc+0Bwoq6qH3XUCJAFrVfUud11lYC3QXlVHh+3zeyBJVe+OEkOgCxBwJip7/PHHeffddwF4\n8skneeWVVyhatKjPkRljjMltbAyI40bgx/AVqnoI2AQ0iWi3OVR8uO0Up9hoICJF3NX3AwIsiDjO\nAqChiMRlbfg5IzY2lhEjRjBs2DBiYmJ48803ueGGG/jiiy9yLAbrW/XG8uSd5coby5N3lqtgCXoB\ncpjoMaYC5UWkmLt85CztCgJXustVgJPALxHtNgIxwLXnGrBfRIQuXbqwaNEirrnmGn766ScSExNp\n0qQJS5cu9Ts8Y4wx5neC3gXzKVAVuFhVU9x1JYDtQGHgElXdJiLvAG2BCqq6x21XEKewuBhnHMjX\nIjIHqKmqJSKO0wEYAdylqp9GPBf4LphIycnJDBgwgIEDB3LkyBEAatSoQfPmzbnnnnu45pprcHqo\njDHGmN+zMSCAiNwEfA28ATyLc5XiLaAZEA+UUdV9InIp8APwCdAROAH8FegBFANqqOrK/FKAhOzY\nsYN//OMfvPfee+zdu/fU+pIlS1KlShWuvPJKLrjgAsqWLUtcXByFCxemUKFCFCpUiMKFC3PeeedR\nokQJSpYsyUUXXUShQoV8fDXGGGNyghUgLvdOmGeAisB+4N/AVcBDqlo8rF1VoD9wPXAQmIFTiLwA\nlFPV3SLyb5w7aGLDqwoR6QW8Dtysqssjjq+PPPIICQkJgPPPu1q1atSrVw843acY5OVjx45x6NAh\npk+fzieffML+/fvJqIIFC5KQkMAll1xCtWrV6N69OxUrVvxdn2q9evUC8XqDvDx48OBcd/74tRx5\nbvkdT1CXV61aRY8ePQITT5CX7fcv+nLo+6SkJADGjBljBUhaRGQWEKOqDdNpNwxooKpXuct9gFdw\n7qz5JazdEKATUEpVj0bsI9deAYlGVdm+fTurV6/ml19+YceOHezatYvk5GSOHz9+6pGcnMyhQ4c4\ncOAAe/fuZcuWLUTm4ZZbbqFLly60aNGCr7766tRJbdK2cOFCy5NHlitvLE/eWa68sSsggNu1UlJV\nV4etKwFsxrmVdrK7rgxQSVUXh7UrCPwMDFfVge66q4B1wKOqOias7ffARlW9J0oMeaoAyazk5GT+\n+9//smjRIubPn8+MGTM4etSp1a688koGDBjA/fffb2NLjDEml7MCBBCRtkBvnEnD9olIPDAKKKiq\nD4S1qwdMBG5S1V9EJBYYCNwMJKrqibC2w4H6nJ6IrB0wFKgVXuiEtbcCJIrDhw8zfvx4BgwYwIYN\nGwBo3Lgx//rXv7jwwgt9js4YY0xm2TwgjjXAVmCliKwEPge+xRnHEW4TsBT4UkRWud//htP9ciKi\nbXdgMrBIRNYAHYCG0YoPk7ZixYrRsWNHfvzxR4YOHUpcXBwzZ86katWqzJkzx+/wAiu8z9WcneXK\nG8uTd5arYInxO4CzUdVvgTs8tNsInNF9kkbbFJw7ZP56btEZgJiYGLp27Ur58uUZPnw48+bN4667\n7mLIkCE8/vjjfodnjDEmoALdBRME1gXj3cmTJ3nhhRd44YUXAOjTpw+vvPKKjQsxxphcxMaABIQV\nIBk3atQoOnbsyMmTJ3nyySd54403rAgxxphcwsaAmFwlvG+1Xbt2TJs2jdjYWN5880369u17xi28\n+ZX1QXtnufLG8uSd5SpYrAAx2aJx48ZMnjyZmJgYXnvtNV577TW/QzLGGBMg1gWTDuuCOTeTJk2i\nZcuWAEyePJlmzZr5HJExxpizsS4Ykye0aNGCAQMGANCmTRv7ZF5jjDGAFSAmi5ytb/Xpp5+mQ4cO\nJCcn06RJE7Zu3ZpzgQWM9UF7Z7nyxvLkneUqWKwAMdlORBg2bBj169dnx44dtGzZkhMnIueHM8YY\nk5/YGJB02BiQrLNz506qV6/O1q1b6dmzJ4MGDfI7JGOMMRFsHpCAsAIka3399dckJiaSkpJig1KN\nMSaAbBCqyVW89q3Wrl2b119/HYBHH32U9evXZ2NUwWN90N5ZrryxPHlnuQoWK0BMjnviiSdo3rw5\nhw4dokWLFiQnJ/sdkjHGmBxmXTDpsC6Y7HHw4EGqV6/Ohg0b6NKlC8OGDfM7JGOMMdgYkMCwAiT7\nrFixglq1anH8+HEmTpxIixYt/A7JGGPyPRsDYnKVzPStVq9e/dSdMB06dOC///1vFkcVPNYH7Z3l\nyhvLk3eWq2CxAsT46vHHH6dZs2Y2HsQYY/IZ64JJh3XBZL8DBw5QvXp1/ve//9G1a1eGDh3qd0jG\nGJNv2RiQgLACJGcsX76c2rVrc+LECZsfxBhjfGRjQEyucq59qzfeeCP/+Mc/AGjfvj0bNmzIgqiC\nx/qgvbNceWN58s5yFSxWgJjA6NatG02bNuXgwYO0bNmSY8eO+R2SMcaYbGJdMOmwLpictX//fm64\n4QaSkpLo3r07//znP/0OyRhj8hUbAxIQVoDkvGXLllGnTh1OnDjBhx9+SNOmTf0OyRhj8g0bA2Jy\nlazsW73pppsYOHAg4HxezP/+978s27ffrA/aO8uVN5Yn7yxXwWIFiAmkJ554gvvuu48DBw7QqlUr\njh8/7ndIxhhjspB1waTDumD8s2/fPm644QY2bdpk40GMMSaHWBeMyfdKlSrFpEmTiImJYciQITZB\nmTHG5CFWgJgskV19qzfffDPvvvsu4HTLTJ06NVuOk1OsD9o7y5U3lifvLFfBYgWICbxHHnmEl156\nCVXlz3/+M19++aXfIRljjDlHgR8DIiJ1gb8BFwOxwA/AU6q6LqLdVcDfgT8AJwAF3lHVYRHtYoHn\ngGZACnBaDNYpAAAgAElEQVQQ6K2qi9I4vo0BCQBVpUuXLrzzzjsUK1aMWbNmceutt/odljHG5Dk2\nBgQQkURgHjBTVSupagKQBHwhIheHtSsBzAVKAjeoahWgDzBERHpE7HYI0By4RVWvB0YCc0Skana/\nHpN5IsLQoUN56KGHOHz4MI0aNWLBggV+h2WMMSaTAl2AAC8C21X19bB1fYF44NmwdXVwrpAMVdVk\nAFWdDawBWoUaiUhloCPwqqrucdu9B2zEuXpiMikn+lZjYmIYPXo0bdu25ejRozRq1IgxY8Zk+3Gz\nkvVBe2e58sby5J3lKliCXoDcCPwYvkJVDwGbgCZhq1Pcr7ER28fidMeE3A8IEPnWeQHQUETizjVg\nk70KFizIe++9x5NPPsnx48dp27Ytffr0ISUlJf2NjTHGBEagx4CIyE5gtao2iFj/A3A1UEJVD4tI\nQZyumsLAPaq6W0TaAO8ATd2rIYjIv3G6XwqFD+wQkV7A68DNqro84lg2BiSg3nnnHbp160ZKSgo3\n33wzY8eO5aqrrvI7LGOMydVsDIhjJfAHEYkJrXDHe1zuLhYHUNWTwN04XSnbRORX4FWgRaj4cJUB\njkapKA66X0tn/Usw2aVz587MnTuXSy65hKVLl3L99dfzzDPPcPDgwfQ3NsYY46ugFyD9cIqCl0Uk\nVkSKAm9wusvlNwARKQssAeKAsqpaHmgDjBORp3I+7PzHr77VxMREVq9eTdu2bTl+/Divvvoql112\nGc888wwbNmzwJaazsT5o7yxX3lievLNcBUugCxBVXQY0AKrg3H77ObAa+ADnSsY+t+nTwDVAV1U9\n4G47H3gPp3ip4LbbDcSLSOSlpeLu1z3Z9VpM9ilZsiSjRo3im2++4dZbb2X//v28+uqrXHnllVSv\nXp1+/foxffp0tmzZwsmTJ/0O1xhjDBCTfhN/qernOIXHKSIyC1gctup6IFlVt0Zs/jPOa7wO2IZT\nvLQCLgF+CWt3Oc5g1bXRYmjbti0JCQmA88+uWrVq1KtXDzhdUduy/8s1a9bkhRdeYO3atSxevJhp\n06axcuVKVq5cSUiBAgUoW7YsFSpUID4+nmPHjlGkSBEqVqxI8eLF2bdvH/Hx8dx6661ceeWV7N69\nm/j4+CyLN7QuCPkK+nK9evUCFU+Ql0OCEk9Ql0PrghJPUJZD3yclJZGTgj4I9VKgpKquDltXAtgM\ntFfVye660cDDwAWquius7atAb6Caqq52JytbBzyqqmPC2n0PbFTVe6LEYINQc6nk5GTmzp3LV199\nxZIlS1i7di07d+7M0D4KFChA1apVqVevHi1atKBmzZqceQHNGGPyjpwahBr0AqQtTgFRR1X3iUg8\nMAooqKoPhLWrCXwJ/BvoqKonROR6nCsn36lq/bC2w4H67j73iEg7YChQK7zQCWtvBYgH4e8qgiw5\nOZmtW7eyf/9+jhw5cupx6NAhDh48yMGDB9m/fz+bN2/m559/5ocffuDEidN3cleuXJmnnnqKNm3a\nULhw4QwfP7fkKQgsV95YnryzXHmTUwVI0Ltg1gBbgZUisg84CUzGuWX2FFVdIiJ1cKZYXy0iJ3DG\ntwwBXovYZ3egP7DIbXcQaBit+DB5T5EiRbjiiis8tz969CjffPMNM2fOZPz48fz000907NiRl19+\nmTfeeIMmTZrYFRFjjMmEQF8BCQK7AmJCUlJSmDRpEn//+99Zu9YZLvTAAw8wYsQIzj//fJ+jM8aY\nrGFdMAFhBYiJlJKSwttvv81f/vIXDh06xEUXXcRHH33ETTfd5HdoxhhzzmwiMpOrRI7Gz8tiYmLo\n1q0bq1atolatWmzdupXExEQ+/PDDdLfNT3k6V5YrbyxP3lmugsUKEGMyqWLFiixcuJD27dvz22+/\n0axZM4YOHep3WMYYkytYF0w6rAvGpEdVee211+jbty8A//znP+nevbvPURljTObYGJCAsALEePXW\nW2/RrVu3U98//vjjPkdkjDEZZ2NATK5ifavQtWtX3nrrLQC6devGxx9/fEYby5N3litvLE/eWa6C\nxQoQY7LQ448/zosvvoiq8uCDD/LNN9/4HZIxxgSSdcGkw7pgTEapKp06deLdd9+lTJkyLFmyhIoV\nK/odljHGeGJjQALCChCTGSdOnKBJkybMnj2bqlWr8vXXXxMXF+d3WMYYky4bA2JyFetb/b3Y2FjG\njx/PlVdeyXfffUenTp1QVctTBliuvLE8eWe5ChYrQIzJJiVLluSjjz4iLi6OcePG2RwhxhgTxrpg\n0mFdMOZcTZw4kVatWhETE8PChQupU6eO3yEZY0yarAvGmDyiZcuW9OzZk5SUFB588EH27t3rd0jG\nGOM7K0BMlrC+1bMbMGAANWvWZPPmzbRv3x67qpY+O6e8sTx5Z7kKFitAjMkBoUGp8fHxfPzxx6cm\nLDPGmPzKxoCkw8aAmKw0efJkWrRoQaFChViyZAnVqlXzOyRjjPkdGwNiTB7UvHlzOnfuzPHjx2nZ\nsiWHDx/2OyRjjPGFFSAmS1jfqjcLFy7kjTfe4A9/+APr16+na9eufocUWHZOeWN58s5yFSxWgBiT\nw4oWLcrEiRMpWrQo77//Pu+//77fIRljTI6zMSDpsDEgJruMHDmS9u3bExcXx7Jly7j22mv9DskY\nY2wMiDF5Xbt27XjooYc4evQozZs358iRI36HZIwxOcYKEJMlrG/Vm/A8iQjDhw/nmmuuYe3atXTp\n0sXmBwlj55Q3lifvLFfBYgWIMT4qVqwYU6ZMIS4ujrFjxzJy5Ei/QzLGmBxhY0DSYWNATE4YO3Ys\nDz/8MEWKFGHJkiVUqVLF75CMMfmUjQExJh9p06YN7du3Jzk5maZNm7Jnzx6/QzLGmGxlBYjJEta3\n6s3Z8jRkyBCqV6/Ohg0baN68OSdOnMi5wALIzilvLE/eWa6CxQoQYwKiaNGiTJs2jfLly7NgwQKe\nfPJJv0MyxphsY2NA0mFjQExOW7JkCYmJiRw7doxBgwbRs2dPv0MyxuQjNgbEmHyqZs2ap+6G6dWr\nl82UaozJkwJfgIhIXRFZICI/i0iSiMwUkWsi2jwvIptEZGXEY72IpIrI5WFtY0XkRRFZJyJrRGSR\niNTJ+VeWt1jfqjde89S6dWsGDRoEwKOPPsrHH3+cjVEFk51T3lievLNcBUugCxARSQTmATNVtZKq\nJgBJwBcicnFYUwX+qqo3hD+AUcASVd0Y1nYI0By4RVWvB0YCc0Skag68JGM869mzJ88++ywnT56k\nefPmTJw40e+QjDEmywR6DIiIfAEkqOqlYevOA3YAY1S1i7vueuBweKEhIgKsB15W1VHuusrAWqC9\nqo4Oa/s9kKSqd0eJwcaAGN+oKn369GHgwIGICCNGjKBDhw5+h2WMycNsDIjjRuDH8BWqegjYBDQJ\nW7cm4ioHQD2gLDAhbN39gAALItouABqKSFzWhG1M1hARBgwYwEsvvYSq0rFjR/7v//6PlJQUv0Mz\nxphzEvQC5DDRY0wFyotIsbNs2xH4QFV/C1tXBTgJ/BLRdiMQA9jHkWaS9a16k5k8iQjPPvssb7/9\nNjExMQwaNIg777yTbdu2ZX2AAWLnlDeWJ+8sV8ES9AJkJfAHEYkJrRCREkBoUGnxaBuJSCngPmBE\nxFNlgKNR+lQOul9Ln3PExmSTzp07s2DBAi644ALmz5/Pddddx+jRo0lNTfU7NGOMybCgjwG5Cfga\neAN4FucqxVtAMyAeKKOq+6Js9wTQWlX/GLF+DlBTVUtErO+AU6zcpaqfRjxnY0BMoGzbto2OHTsy\na9YsAKpXr85LL71Eo0aNcIY+GWNM5tkYEEBVlwENcLpOfgA+B1YDHwBHohUfrvacefUDYDcQL2f+\nlQ5dSbEP4DCBV6FCBWbMmMH7779PhQoVWLFiBX/605+49tprGTBgAOvXr/c7RGOMSVegr4CkRURm\nATGq2jDKczcDc4ALI8Z/ICJ9gFdw7qz5JWz9EKATUEpVj0Zso4888ggJCQkAlCxZkmrVqlGvXj3g\ndJ9ifl8OrQtKPEFdHjx4cJaeP7Nnz+ajjz5i5syZbN26lZCKFSty0003Ubx4ccqVK8edd95JhQoV\n+PHHH4mLi6N+/fqByMfZliPPLb/jCeryqlWr6NGjR2DiCfJyVv/+5ZXl0PdJSUkAjBkzJkeugAS6\nABGRS4GSqro6bF0JYDPOrbSTo2wzAjihql2jPHcVsA54VFXHhK3/HtioqvdE2ca6YDxYuHDhqZPa\npC278nTixAlmzpzJ1KlTmT59Ovv370+zbaFChShbtizlypWjbNmyv/u+XLlyXH/99VSvXp2YmJg0\n95ET7JzyxvLkneXKm5zqggl6AdIW6A3UUdV9IhKPM7lYQVV9IEr7eGAbUFdVv0tjn8OB+u4+94hI\nO2AoUCu80AlrbwWIyVVSUlJYu3Yty5cv5/vvv2fz5s1s3ryZHTt2sHPnTo4ePZruPooXL84dd9xB\n27ZtadSoke/FiDEm51gBAohIDeBVoBKwD+cW2snA66p6Mkr7R4HOqlrzLPuMAfrjzIZ6AucOmN6q\nuiiN9laAmDzl6NGj7Nq169Rj586dp77fvn07ixcv5r///e+p9hUqVKBv37507NiRIkWK+Bi5MSYn\nWAESEFaAeGOXNr3JLXnatGkTEyZMYNSoUfz0008AXHbZZQwdOpS77z5jwuBskVty5TfLk3eWK2/s\nLhhjjG8uu+wy+vTpw7p16/j444+pUqUKmzZt4p577uHhhx/m8OHDfodojMnl7ApIOuwKiDHOuJIh\nQ4bQr18/jh49yjXXXMOUKVO49lqbPNiYvMa6YALCChBjTlu7di3NmjVj3bp1FC9enE8++YTExES/\nwzLGZCHrgjG5Svj95CZtuT1P1157LUuXLuWBBx7g4MGD3HnnnXz00UfZcqzcnqucYnnyznIVLFaA\nGGMypFixYkycOJEuXbpw7NgxmjdvztSpU/0OyxiTy1gXTDqsC8aY6FSVfv368fLLLxMbG8u0adO4\n6667/A7LGHOObAxIQFgBYkzaVJVevXoxePBgihQpwvz586lVq5bfYRljzoGNATG5ivWtepPX8iQi\nDBo0iA4dOpCcnMy9997Lxo0bs2TfeS1X2cXy5J3lKlisADHGnBMRYfjw4TRs2JBdu3bRuHHjs34W\njTHGgHXBpMu6YIzx5sCBA9SuXZu1a9fSsGFDZs2aRcGCBf0OyxiTQdYFY4zJVUqUKMHMmTMpW7Ys\nc+bM4W9/+5vfIRljAswKEJMlrG/Vm7yep4SEBMaPH0+BAgV48cUXmTlzZqb3lddzlVUsT95ZroLF\nChBjTJa6/fbbeemllwBo06ZNlg1KNcbkLTYGJB02BsSYjEtNTeW+++5j+vTpVK9enUWLFlGkSBG/\nwzLGeGBjQIwxuVaBAgUYM2YMFStWZMWKFTz55JN+h2SMCRgrQEyWsL5Vb/JTnkqVKsWUKVMoXLgw\nI0aMYNy4cRnaPj/l6lxYnryzXAWLFSDGmGxzww038M9//hOAzp07s27dOp8jMsYEhY0BSYeNATHm\n3Kgqbdq0Ydy4cVx33XUsWbKE+Ph4v8MyxqTBxoAYY/IEEeHtt9/m6quv5ocffuDxxx/HinpjjBUg\nJktY36o3+TVPxYoVY8qUKRQtWpT333+fUaNGpbtNfs1VRlmevLNcBYsVIMaYHHHdddcxfPhwALp2\n7crq1at9jsgY4ycbA5IOGwNiTNZq3749I0eOpFKlSixfvpzixYv7HZIxJkxOjQGxAiQdVoAYk7WO\nHj3KH//4R9asWUPLli0ZP348Itn+t84Y45ENQjW5ivWtemN5gri4OCZPnkyxYsWYOHHiqW6ZSJYr\nbyxP3lmugsUKEGNMjqtcuTLvvvsuAD179mT58uU+R2SMyWnWBZMO64IxJvt07dqVYcOGkZCQwIoV\nKyhVqpTfIRmT79kYkICwAsSY7HPs2DHq1KnDt99+S/369Zk9ezaFChXyOyxj8jUbA2JyFetb9cby\n9HuFCxdm6tSplC9fngULFtCxY8dTk5RZrryxPHlnuQoWK0CMMb669NJLmTFjBnFxcbz//vv87W9/\n8zskY0wOCHwXjIjUBf4GXAzEAj8AT6nqGZ9qJSL1gGeBUkAJ4DdgtKoOCmsTCzwHNANSgINAb1Vd\nlMbxrQvGmBwwY8YM7r33XlJTUxk0aBA9e/b0OyRj8iXrggFEJBGYB8xU1UqqmgAkAV+IyMURbZsB\n44D/U9UbVbUSMAu4O2K3Q4DmwC2qej0wEpgjIlWz9cUYY87q7rvvZsSIEQD06tWLt956y+eIjDHZ\nyXMBIiIVszOQNLwIbFfV18PW9QXica50ACAi5wFvA8+pavj8zq8A/cPaVQY6Aq+q6h4AVX0P2Aj8\nPbteRH5gfaveWJ7Orn379gwbNgyAbt26MXDgQPvgunTYOeWd5SpYMnIF5OtsiyJtNwI/hq9Q1UPA\nJqBJ2Or7gPOBGRFtD6jql2Gr7gcEWBBxnAVAQxGJy6K4jTGZ1KVLF4YMGQJA7969efLJJzl58qTP\nURljslpGCpByIrJBRPpFdn9ko8NEjzEVKC8ixdzl2sAx4HIRmS0iP4jItyLytIiEb18FOAn8ErG/\njUAMcG3Whp9/1KtXz+8QcgXLkzfdunVjwoQJFCpUiCFDhtC4cWN27tzpd1iBZOeUd5arYMlIAbID\n50pDWeBb9x99c3dQZ3ZZCfxBRGJCK0SkBHC5uxj6FKtLcAqI0UAPVb0O6An0w+maCSkDHI0yqvSg\n+7V0lkZvjMm0li1bMmfOHEqXLs2nn35K1apV+eyzz/wOyxiTRTJSgDymqmtU9Umcf/jvAY8CW0Rk\nsIhUyYb4+uEUBS+LSKyIFAXewLl7BZy7XACKAAWBQar6I4CqfgH8C2gvIpdlQ2wmjPWtemN58m7h\nwoUkJiayatUq6taty6+//krDhg1p3bo1W7Zs8Tu8wLBzyjvLVbDEpN/EoarTwr4/DkwGJruDU8cB\n3UVkBc5dJf9W1QPnGpyqLhORBsAzOLff7gf+DSQDbVR1n9v0kPt1VcQuVuGM+aiBM25kNxAvZ95b\nG7qSsidaHG3btiUhIQGAkiVLUq1atVOX8kIndH5fDglKPEFdXrVqVaDiyS3L8+bNY+DAgfTv35/x\n48fz4Ycf8uijj3LzzTdz+eWX+x6fn8urVq0KVDxBXrbfv+jLoe+TkpLISZ7nARGR21R1ftjytUB7\noA1O10Yy8DHO1Ym7gClA36woRKLEMguIUdWG7vIrQB+gtqp+E9buQZziqIWqThGRvsDLQIKq/hLW\nbgjQCSilqkcjjmXzgBgTEElJSfTt25dJkyadujumevXqNG3alPr163PjjTfaVO7GnKPAfRaMiGwH\nrgJa43S93OQ+FX7VY7/b9jzgKZyC4I5MBydyKVAy/NZadwzIZqC9qk5219UCFuF0E40Ia/t3nNt2\nr1DVJBG5ClgHPKqqY8LafQ9sVNV7osRgBYgxAbNu3TrefPNNJkyYwIEDp9/jxMTEULFiRSpXrsxF\nF11EmTJl0nzEx8f7+AqMCa4gFiAncO40icPpqhgHjIyYdyNymx2qekGmgxNpC/QG6qjqPhGJB0YB\nBVX1gYi2HwJ/AOqp6nYRuQb4Apiqqp3D2g0H6rv73CMi7YChQK1or8UKEG8WLlx46rKeSZvlyTsv\nuUpOTmbWrFnMnTuXBQsW8NNPP3meN6R48eJcd911XH/99dSqVYtGjRpRvnz5LIg8Z9k55Z3lypuc\nKkA8jwHBGbD6Oc7Vjk9U9URaDUWkCDAMZ/zFuVgDbAVWisg+nFtoJwOvR2n7Z+AF4GsRSXbXDXQf\n4brjTE62yC2qDgINz1ZIGWOCqUiRIjRt2pSmTZsC8Ntvv/Hzzz+zfv16duzYwe7du9mzZw+7du1i\nz5497N69m127drFr1y4OHjzI4sWLWbx48akZWGvWrEmnTp1o1aoVcXE2LZAx2SkjV0B+VVVPbw9E\npBBOV80OVf3POcTnO7sCYkzeo6rs2rWLNWvW8N133zFv3jwWLFjAb785N9aVKlWKvn370r17d4oW\nLepztMbkrCB2wRRU1Xw3HaEVIMbkD0ePHmXSpEkMGzaMZcuWAc4n9Y4YMYI777zT5+iMyTmB+zC6\n/Fh8GO/Cb+cyabM8eZfTuYqLi6Nt27YsWbKE2bNnU7VqVX755RcaNWrEY489xrFjx3I0Hq/snPLO\nchUsGZmIzBhj8jwR4c4772T58uW8+uqrFC5cmHfeeYdbb72VzZs3+x2eMXmG5y6Y/Mq6YIzJ31as\nWEHTpk3ZtGkTF1xwAXPmzKFKleyY+NmYYAhcF4wxxuRH1atX59tvv+W2225jx44dJCYmsnjxYr/D\nMibXswLEZAnrW/XG8uRdkHJVunRpZs2axf3338/+/ftp0KABixYt8jssIFh5CjrLVbBYAWKMMR4U\nLlyYSZMm8dBDD3H06FEaN27MypUr/Q7LmFzLxoCkw8aAGGPCpaSk0KpVKz788EPKli3LokWLqFSp\nkt9hGZNlAjcPSH5lBYgxJtKxY8do0qQJc+bM4aqrruKbb76hVKlSfodlTJawQagmV7G+VW8sT94F\nOVeFCxdmypQpVKlShfXr19OsWTNOnEjz0ymyVZDzFDSWq2CxAsQYYzLhvPPOY/r06VxwwQXMnz+f\nbt26+R2SMbmKdcGkw7pgjDFns3TpUhITE0lOTubdd9+lffv2fodkzDmxMSABYQWIMSY977//Po88\n8giFCxfm66+/pnr16n6HZEym2RgQk6tY36o3lifvclOuHn74YTp37syxY8d44IEH2Lt3b44dOzfl\nyW+Wq2CxAsQYY7LA4MGDufHGG0lKSqJNmzakpqb6HZIxgWZdMOmwLhhjjFdJSUnUqFGDvXv38uKL\nL9KvXz+/QzImw2wMSEBYAWKMyYjZs2fzpz/9CRHhs88+47bbbvM7JGMyxMaAmFzF+la9sTx5l1tz\n1ahRI/r160dqaioPPvgg27Zty9bj5dY8+cFyFSxWgBhjTBbr378/t99+Ozt37qRVq1akpKT4HZIx\ngWNdMOmwLhhjTGbs2LGDG264ge3bt9O7d28GDBjgd0jGeGJjQALCChBjTGZ9+eWX1K9fn5MnTzJt\n2jSaNGnid0jGpMvGgJhcxfpWvbE8eZcXcnXrrbfyyiuvAPDII4+wcePGLD9GXshTTrFcBYsVIMYY\nk42eeuopmjRpwv79+2nevDnHjh3zOyRjAsG6YNJhXTDGmHO1b98+qlevTlJSEl26dGHYsGF+h2RM\nmmwMSEBYAWKMyQrffvsttWvX5vjx44wbN47WrVv7HZIxUdkYEJOrWN+qN5Yn7/JarmrUqMGbb74J\nQKdOnVi3bl2W7Dev5Sk7Wa6CxQoQY4zJIZ07d6Z169YcOXKEe+65h927d/sdkjG+sS6YdFgXjDEm\nKx05coS6deuyYsUK6tSpw9y5cylSpIjfYRlzinXBGGNMHhQfH8/06dO5+OKLWbRoEe3atbNPzjX5\nUuALEBGpKyILRORnEUkSkZkick2UdsdFZGWUx9UR7WJF5EURWScia0RkkYjUyblXlDdZ36o3lifv\n8nKuKlSowMyZMylWrBgTJkygR48eZPZKa17OU1azXAVLoAsQEUkE5gEzVbWSqiYAScAXInJxRPOt\nqnpDlMePEe2GAM2BW1T1emAkMEdEqmbvqzHGmNOqVKnC1KlTKVSoEEOGDKFv376ZLkKMyY0CPQZE\nRL4AElT10rB15wE7gDGq2iVs/UZVvTyd/VUG1gLtVXV02PrvgSRVvTvKNjYGxBiTbT755BMeeOAB\nUlJS6NOnD6+88goi2d79bkyabAyI40bgd1cwVPUQsAnIzIcq3A8IsCBi/QKgoYjEZSZIY4zJrCZN\nmjBu3DgKFizIgAED6NChg316rskXgl6AHCZ6jKlAeREpFrYuTkSGicjXIrJeRD4WkVsitqsCnAR+\niVi/EYgBrs2qwPMb61v1xvLkXX7KVYsWLZg2bRpFixZl5MiR3HPPPezdu9fTtvkpT+fKchUsQS9A\nVgJ/EJGY0AoRKQGEulqKh7U9AkxR1do4hcZaYKGIhF8pKQMcjdKnctD9WjorgzfGGK8aN27MvHnz\nKF26NLNnz6Z69eosXbrU77CMyTZBHwNyE/A18AbwLM5VireAZkA8UEZV951l+zVAYVW9yl2eA9RU\n1RIR7ToAI4C7VPXTiOdsDIgxJsds2rSJ5s2bs2zZMmJiYujVqxfPPfcc8fHxfodm8omcGgMSk34T\n/6jqMhFpADwD/ADsB/4NJAMPna34cC0D2opIKbftbiBezqwqQldS9kTbSdu2bUlISACgZMmSVKtW\njXr16gGnL+nZsi3bsi1n1fKXX35J3759GTx4MK+99hoTJ07kueee45JLLiE2Ntb3+Gw5by2Hvk9K\nSiInBfoKSFpEZBYQo6oN3eV44KSqJke0+xfQHiinqrtFpC/wMs6dNb+EtRsCdAJKqerRiH3YFRAP\nFi5ceOqkNmmzPHlnuYKlS5fSqVMnvvvuOwAuueQSOnXqROvWralYsSJgecoIy5U3dhcMICKXikiV\niHUlgFuAf4WtfhroEWUXNYAtqhr6wIWpgAL1I9rVB+ZEFh/GGOOnm2++meXLl/PBBx9w7bXXsnnz\nZv76179yxRVXULVqVZ544gnmzZvHihUrOHz4sN/hGpMhgb4CIiJtgd5AHVXd517pGAUUVNUHwtr1\nBx4FblPVDe66p4DXgEcj5vwYjlNw1FHVPSLSDhgK1FLV1VFisCsgxhjfpaamMmfOHMaNG8dHH33E\nkSNHzmhTsmRJypYte+pRrlw5LrzwQipUqECFChWoVKkSlSpVomDBgj68ApNb5NQVkKAXIDWAV4FK\nwD6cW2gnA6+r6smwdgnAY0BDnHk+SuPcavsPVf0oYp8xQH+c2VBP4NwB01tVF6URgxUgxphASU5O\nZivxajoAAB+gSURBVOnSpXz++eesXLmS9evX8/PPP3P8+PF0ty1atChVqlTh9ttv5+6776ZmzZoU\nKBDoi+Emh1kBEhBWgHhjfaveWJ68s1x5E8pTamoqe/fuZdeuXaceO3fuZPv27Wzbto2tW7eybt06\nNm3a9LvtL7/8ch577DHat29P6dJ5eyYCO6e8sbtgjDHGeFagQAHKlClDmTJluOaaMz6v85R9+/ax\nZMkSZs2axbRp09i4cSN9+vThpZdeonfv3vTs2dNu+TU5wq6ApMOugBhj8qqTJ08ye/Zs3nzzTT77\n7DMAKlasyKhRo6hbt67P0Rm/WBdMQFgBYozJDxYuXEiPHj1O3fLbq1cvBgwYQEyMXSjPb+w2XJOr\nhE9oY9JmefLOcuVNVuWpXr16LF26lP79+xMTE8OgQYO488472bMn6vyMuZKdU8FiBYgxxhgAChUq\nxPPPP8/ChQspV64c8+fPp2bNmmzcuNHv0EweZF0w6bAuGGNMfrR582buvfdeVq5cycUXX8zcuXOp\nXLmy32GZHGBjQALCChBjTH518OBBGjduzFdffUW5cuX4/PPPufrqq/0Oy2QzGwNichXrW/XG8uSd\n5cqb7MxT8eLFmT17Ng0aNGDnzp00bNiQLVu2ZNvxspudU8FiBYgxxpg0xcfHM23aNGrXrs3mzZtp\n2LBhnhqYavxjXTDpsC4YY4yBvXv3UrduXX744Qdq1arF/PnzKVKkiN9hmWxgXTDGGGMC4/zzz+fT\nTz/l0ksvZfHixTz22GPYmzNzLqwAMVnC+la9sTx5Z7nyJifzdNFFFzFt2jTi4uIYM2YMgwcPzrFj\nZwU7p4LFChBjjDGeVatWjdGjRwPw1FNP8emnn/obkMm1bAxIOmwMiDHGnKl///688MILlCxZkm+/\n/ZaKFSv6HZLJIjYPSEBYAWKMMWdKTU3l/vvv55NPPqFGjRp89dVXNig1j7BBqCZXsb5VbyxP3lmu\nvPErTwUKFGDMmDFcfvnlfPvtt/Tq1cuXODLCzqlgsQLEGGNMppQsWZLJkydTqFAhhg8fzr///W+/\nQzK5iHXBpMO6YIwx5uzefvttunTpQnx8PMuWLeOaa67xOyRzDmwMSEBYAWKMMWenqv/f3r2HSVFe\n+x7/LmBQGLmooKKCjEQEoiNPREVUkEAIKIgg3iK4UTT75CSbbaKHeAlgNBoT3ficQ7xGUYN4Q9Co\njIKXGRBwthcuY0BB5SJuwSsKgiAy6/xRNdg2PUyBPV0107/P8/TT1NtvV69evAxr6n2rigsuuICH\nH36YLl268Oqrr1JYWBh3WLKHtAZE6hTNrUajPEWnXEWThDyZGXfffTedOnVi6dKl/OpXv0rkRcqS\nkCv5jgoQERH5wfbZZx8ef/xxmjZtyuTJk7nrrrviDkkSTlMwNdAUjIhIdFOmTGH48OE0btyYl19+\nmeOPPz7ukGQ3aQpGRETqnAsuuIBf//rXfPPNNwwbNoxPP/007pAkoVSASFZobjUa5Sk65SqaJOZp\nwoQJdO/enTVr1vCLX/yC7du3xx0SkMxc5TMVICIiklWNGzdm6tSptGrViueff55rr7027pAkgbQG\npAZaAyIismdefPFF+vXrR2VlJU8//TQDBw6MOySJQGtARESkTuvTpw9/+tOfgGBtyJIlS2KOSJJE\nBYhkheZWo1GeolOuokl6nn7/+98zbNgwNmzYwOmnn85HH30UWyxJz1W+qZcFiJk9aGaVZtYz7lhE\nRPJZgwYN+Mc//sHxxx/P6tWrGTx4MF9//XXcYUkCJH4NSFhE/BE4FCgAlgBXuPtb1fTvBrwKONDb\n3eekvV4AjAOGAd8CG4Ax7j6vmv1pDYiIyA+0bt06TjjhBN5//30GDRrEtGnTKCgoiDssyUBrQAAz\n6wW8CMxw9yPcvT2wCphjZodW87b/AmYA1SVvInA2cLK7Hw1MAmaZ2THZjF1ERL5z0EEHUVJSwr77\n7svTTz/N8OHDE3N6rsQj0QUIcD2w1t1vSWm7EigErknvbGZDCL7T45l2ZmZHApcCN7n7ZwDufi+w\nErghu6HnF82tRqM8RadcRVOX8vTjH/+YWbNm0bx5cx577DFGjRqV0yKkLuUqHyS9AOkGvJ3a4O4b\ngdXAGant4dTKTcDlVH/0Y0j4WmlaeynQz8yaZiFmERGpRrdu3SgpKaFp06Y88MADnH/++WzdujXu\nsCQGiV4DYmYfAxXu3jetfQnQCWjh7l+FbZcB3d39PDMbSTC1cmrqGhAze4hg+qVx6sIOM/sdcAtw\nvLu/nvZZWgMiIpJls2fP5owzzmDDhg307t2bJ554ghYtWsQdlqA1IFUWAkeZWaOqBjNrARSFm83D\ntn2BMQTTM7vSCticoaLYED7v/4MjFhGRGvXq1Ys5c+Zw0EEHUVpaSrdu3Vi8eHHcYUkOJb0A+QNB\nUXCjmRWYWRPgVoKzVwCqzuUaCzzk7qtyH6KA5lajUp6iU66iqct5OuaYY5g/fz5du3bl3XffpXv3\n7tx2221UVlbWyufV5VzVR41q7hIfd3/NzPoCVxGcfvsF8BCwBRju7uvNrANwPtA5wy7SDyF9ChTa\nzvMqzcPnzzLFMXLkSNq3bw9Ay5Yt6dq1K6eeeirw3YDO9+0qSYknqduLFi1KVDzarvvbixYtSlQ8\ne7I9f/58Ro8ezT333MNvfvMbJk+ezJ133skXX3yR1c/Tv7/M21V/XrVqFbmU6DUg1TGzEqCRu/cz\ns0uBqwmKkyr7AW2B94CvgInuPsnMrgRuBNq7+/sp+5sI/BLY1903p32W1oCIiOTAtGnTGD16NB9+\n+CEAQ4cO5eqrr+bYY4+NObL8kqs1IIkuQMysHdDS3StS2loAa4BR7j61mvf9G3AfOy9C7Qi8BVzs\n7g+ktP8LWOnugzLsSwWIiEiObNiwgWuvvZbbb799x9kxXbt2Zfjw4Zx22ml06tQJs1r/vzGvaRFq\n4KfAI+EiU8ysEPg78Hx1xUfI0p4BcPflwN3AVWa2f7jPiwgWte50XRGJLvVQnlRPeYpOuYqmvuWp\nefPmTJgwgZUrV3L55ZfTsmVLFi1axBVXXEGXLl1o06YNAwcOZMyYMUyaNIlZs2ZRUVHBxx9/XOPa\nkfqWq7ou0WtAgDeB/wEWmtl6YDswleCU2Z2El2GfCuxDcCn2R8zsa4Krnn4YdvsPYDwwz8y2EZwB\n0y/1KIuIiMSrTZs23HLLLdxwww2UlJQwdepUSktLWbduHTNmzGDGjBk7vadhw4a0adOGtm3b0rZt\nW9q1a0dxcTHHHnssRx55ZAzfQnYl0VMwSaApGBGRZHB33nnnHSoqKnjrrbdYvnw5a9euZe3ataxb\nt47PP/+82vc2b96cAQMGMHjwYM4880yaNGmSw8jrFq0BSQgVICIidcPWrVv58MMPWbNmDWvWrGHl\nypUsXLiQ119/nfff33HeAfvttx+XXHIJl112GW3atIkx4mRSAZIQKkCiKSsr23Fql1RPeYpOuYpG\neYpmxYoVTJgwgfLyct544w0ACgsLGTNmDJdffjmFhYUxR5gcWoQqIiKSJYcffjjDhg3jtddeo7y8\nnMGDB7Np0ybGjx9PcXExr7zyStwh5h0dAamBjoCIiNRPs2fPZvTo0VRUVNCgQQPGjh3LuHHjaNAg\nv3831xRMQqgAERGpv7Zu3cq4ceO4+eabcXeGDBnC5MmT83pKRlMwUqfo/PpolKfolKtolKfoMuVq\nr7324i9/+QszZ86kRYsWPPHEE/Ts2ZOPP/449wHmGRUgIiKS9372s59RXl5Ohw4dWLBgAb1792bd\nunVxh1WvaQqmBpqCERHJHx999BF9+vRhyZIldOzYkbKysrw7VVdTMCIiIjl24IEHUlpaSnFxMcuX\nL2fAgAF8+eWXcYdVL6kAkazQPHQ0ylN0ylU0ylN0UXPVunVrXnzxRTp27MjixYs588wz2bJlS+0G\nl4dUgIiIiKRp1aoVM2fOpE2bNpSVlTFixIgab3Ynu0drQGqgNSAiIvmroqKCnj178uWXXzJ27Fiu\nu+66uEOqdboOSEKoABERyW+zZs1iwIABVFZW8thjj3H22WfHHVKt0iJUqVM0Dx2N8hSdchWN8hTd\nnuaqX79+3HzzzQCMHDmSRYsWZTGq/KUCREREpAa//e1vufDCC9m8eTNnnnkmn376adwh1XmagqmB\npmBERARgy5Yt9OrVi1dffZX+/fszY8aMennfGE3BiIiIJMjee+/N1KlT2X///Xnuuee44YYb4g6p\nTlMBIlmheeholKfolKtolKfospGrdu3a8eCDD2JmjB8/nhdeeOGHB5anVICIiIjshv79+zN27Fjc\nnfPPP58PPvgg7pDqJK0BqYHWgIiISLrt27czYMAAnn/+eXr06EFZWRkFBQVxh5UVWgMiIiKSUA0b\nNmTKlCkccsghzJ8/nyuvvDLukOocFSCSFZqHjkZ5ik65ikZ5ii7buWrdujWPPfYYjRo1YsKECUyf\nPj2r+6/vVICIiIjsoR49enzvImXLly+POaK6Q2tAaqA1ICIisivuzjnnnMPjjz9Oly5dKC8vp1mz\nZnGHtce0BkRERKQOMDMmTZpE586dWbp0KSNHjkS/uNZMBYhkheaho1GeolOuolGeoqvNXDVr1own\nn3ySFi1aMH36dP785z/X2mfVFypAREREsqBjx45MmTIFM+MPf/gDzz77bNwhJZrWgNRAa0BERGR3\nXH/99YwbN47mzZszb948jjrqqLhD2i25WgOiAqQGKkBERGR3VFZWcs455zBt2jTatm1LeXk5Bx98\ncNxhRaZFqCEz62lmpWb2jpmtMrMZZtY5rc8BZnadmZWb2etm9paZvWZmF2bYX4GZXR/2edPM5pnZ\nSbn7RvWT5qGjUZ6iU66iUZ6iy1WuGjRowOTJk+nRowdr1qzh9NNPZ+PGjTn57Lok0QWImfUCXgRm\nuPsR7t4eWAXMMbNDU7oeD/wn8Gt37+bunYF7gfvN7IK03U4EzgZOdvejgUnALDM7pna/jYiI5Ism\nTZrwz3/+kyOOOIJFixYxePBgNm/eHHdYiZLoKRgzmwO0d/d2KW3NgI+AB9z9V2Fbd+Cn7n5j2vvX\nAy+6+7Bw+0hgKTDK3e9P6fcvYJW7D8wQg6ZgRERkj7z33nucfPLJrFu3jr59+/LUU0/RpEmTuMPa\nJU3BBLoBb6c2uPtGYDVwRkpbeYbiowBoAnyS0jwEMKA07XNKgX5m1jR7oYuISL7r0KEDpaWlHHjg\ngbzwwgsMHjyYr7/+Ou6wEiHpBchXZI6xEjjIzPbJ9CYz2w+4A1gD3JDyUjGwHXg/7S0rgUZAlx8a\ncL7SPHQ0ylN0ylU0ylN0ceWqU6dOvPTSSxxwwAE8//zz9OnTh08++aTmN9ZzSS9AFgJHmVmjqgYz\nawEUhZvNUzubWRMzewv4mKCYGOzuH6R0aQVszjCnsiF83j+bwYuIiAB06dKF0tJS2rVrxyuvvEKP\nHj1YtmxZ3GHFKulrQI4D5gO3AtcQHKW4DRgGFAKt3H19hvc1JliUei0wwt2nh+2zgBPcvUVa/0uA\nu4EB7j4z7TWtARERkaxYu3YtAwcOZMGCBRQWFnLHHXcwYsSIuMP6nlytAWlUc5f4uPtrZtYXuApY\nAnwBPARsAYZnKj7C930D3GxmfYC/m9mT7l4JfAoU2s5VRdWRlM8y7W/kyJG0b98egJYtW9K1a1dO\nPfVU4LtDetrWtra1rW1tR9mePXs2l156KY888ggXXnghJSUlTJgwYccRkVzHU/XnVatWkVPuXuce\nQAkwK2W7MdAgQ7+/EqwX6RBuXxlut0vrNxHYCjTNsA+XmpWWlsYdQp2gPEWnXEWjPEWXpFxVVlb6\nvffe602bNnXAmzVr5n/961/9q6++ijs0D//fq/X/yxO9BsTM2plZcVpbC+Bk4O8pzXcD52TYRXuC\nguPzcHs64EDvtH69CQoanaQtIiK1zsy4+OKLefPNNxk0aBAbN25kzJgxtGvXjvHjx+f+aEQMkr4G\nZCQwBjjJ3debWSFwH9DQ3c9K6Xcf0BUY5OGiUzMbCjwKTHH3kSl97yAoOE5y98/M7CLgb8CJ7l6R\nIQZPco5ERKTue/bZZ7nuuusoLy/f0da9e3dOO+00TjnlFI477jgKCwtzEovuBQOY2bHATcARwHqC\nU2inAre4+/aUfkcB/w70JDjisRewieAqp3en9W0EjCe4Guo2gjNgxrj7vGpiUAEiIiK1zt2ZO3cu\nt99+O0899dROV05t164dnTt3pkOHDhQVFXH44YdTVFREUVERLVu2zFocKkASQgVINGVlZTsWNkn1\nlKfolKtolKfo6lKuNm3axMyZM5k9ezZz5sxhyZIlbNu2rdr+++67745ipKioiOLiYvr27UubNm12\n+7N1FoyIiEieKiwsZOjQoQwdOhSAb7/9lhUrVvD222+zYsUKVq5cueN55cqVrF+/nvXr17NgwYLv\n7adr165ceOGFDB8+nNatW8fxVaqlIyA10BEQERFJMnfnk08++V5hMn/+fMrKynZM4xQUFDBq1Ciu\nueYaDj300F3uT1MwCaECRERE6qKtW7dSUlLCvffeS0lJCe5O48aNueKKKxg7dix77713xvfpZnRS\np6Re0EaqpzxFp1xFozxFl2+52muvvRgyZAjPPPMMS5Ys4bzzzmPbtm3ceOONdO3alVdeeSXW+FSA\niIiI1HOdO3fm4YcfZu7cuXTq1Illy5bRs2dPJk6cSFxH+TUFUwNNwYiISH2yZcsWrr76am699VYA\nRowYwT333EPjxo0BrQFJDBUgIiJSHz366KOMGjWKTZs28fOf/5xp06ZRWFioNSBSt+Tb3OqeUp6i\nU66iUZ6iU66+79xzz2XOnDm0bt2amTNn0rdvX7788sucfb4KEBERkTz1k5/8hLlz53LYYYdRXl5O\n//79c/bZmoKpgaZgRESkvlu9ejW9evVi9erVAJqCERERkdp32GGHUVpayoknnpizz1QBIlmhudVo\nlKfolKtolKfolKtdKyoqYt68jPdlrRUqQERERAQITsHN2WdpfcOuaQ2IiIjkE52GKyIiIvWWChDJ\nCs2tRqM8RadcRaM8RadcJYsKEBEREck5rQGpgdaAiIhIPtEaEBEREam3VIBIVmhuNRrlKTrlKhrl\nKTrlKllUgIiIiEjOaQ1IDbQGRERE8onWgIiIiEi9pQJEskJzq9EoT9EpV9EoT9EpV8miAkRERERy\nTmtAaqA1ICIikk+0BkRERETqLRUgkhWaW41GeYpOuYpGeYpOuUoWFSAiIiKSc4lfA2JmPYE/AocC\nBcAS4Ap3fyulz4+A/w30AyqBRsCrwHh3X522vwJgHDAM+BbYAIxx93nVfL7WgIiISN7QGhDAzHoB\nLwIz3P0Id28PrALmmNmhKV3/BnQDTnH3YuBE4HDgNTM7MG23E4GzgZPd/WhgEjDLzI6p1S8jIiIi\nOyS6AAGuB9a6+y0pbVcChcA1KW0OXOvu6wHc/UvgWqAVcEFVJzM7ErgUuMndPwv73gusBG6ova9R\n/2luNRrlKTrlKhrlKTrlKlmSXoB0A95ObXD3jcBq4IyU5kHu/lLae9eGzy1T2oYABpSm9S0F+plZ\n0x8ccZ5atGhR3CHUCcpTdMpVNMpTdMpVsiS9APmKzDFWAgeZ2T4A7v5thj4dw+eylLZiYDvwflrf\nlQTrRrr8kGDz2RdffBF3CHWC8hSdchWN8hSdcpUsSS9AFgJHmVmjqgYzawEUhZvNd/HeXwKz0o6M\ntAI2Z1hVuiF83v8HxisiIiIRJL0A+QNBUXCjmRWYWRPgVoKzVwC+zvQmMxsJdAb+LRdBCqxatSru\nEOoE5Sk65Soa5Sk65SpZ6sJpuL2AqwjOavkCeIhgemW4u+90BMTM+gF3A33d/d201x4CzgEKUo+C\nmNnvgFuA49399bT3JDtBIiIiWZaL03Ab1dwlXu4+G5id2mZmJUB5el8z6wvcAfRLLz5Ci4HzgLZ8\nfx1IEbANWJrh82v9L0FERCTfJHoKxszamVlxWlsL4GTg72ntfYC7gAHuvjxsO9bMxqV0e4LglN3e\naR/Vm2C9yOYsfwURERHJINFTMOFajjHASe6+3swKgfuAhu5+Vkq/nwJPAxOAZSm76Awc7O4XpfS9\ng6DgOMndPzOziwguZHaiu1fU9ncSERGR5E/BvAn8D7DQzNYTnEI7lWC9Rqqbgb35/sXJIDja8UBa\n238A44F5ZraN4AyYfio+pCZm1oagAO7n7ok+ehgn5UlEInF3PdIewAHAFIKLoL1NUPQcEndcOfru\n7Qmuv7Iww6NFSr99CI4cvU1wf56ZQJcM+ysguKLtWwQF5TyCo0+xf9fdzMtQYEX4fbfvol/W8wJc\nFu5rMfAGMDjufGQhT99UM8Y65UmeuhJMIy8FKsK4/y/QSuNpj/KU1+MpjLcDwS/nr4ePZcAc4LSk\njqnYk5a0B9A4TOKjBGtkGgD3A8uBwrjjy8H3bw+URuj3bDi49w63rwM+JpjySu13ZzjQ9w+3RwGb\ngGPi/q67mZd5YW7uBypzlReCWw98AhSF233DH7b9487JD8zTyoj7q695qvrFpkm4fXD4g35Z1djR\neNqtPOX1eApj/A2wBjg83DbgRoLLVvRM4piKPWlJexDcK6YSaJ/SdmD4l3hF3PHl4Pu3p4YCBPhZ\nmKNTU9oKgM+Av6W0HUkwbTYy7f3/Ap6J+7vuZl6q1kvdTzX/sWY7LwS3EdhEcJ+j1H7PAP+KOyd7\nmqfw9ZUR9lWf87S06j+KlLaLw/EzVOMpep40nnbEdyZwcVpbizBXtyRxTGl+dmdnAavdfVVVg7t/\nRPAP4azq3pRnziKocOdWNbj7NoLfflNzVG/uvePhv6oaZDsv/YEm1fTrEt5cMVEi5imqepsnoNjd\nV6S1pd+/Ku/HE9HyFFV9zhPu/qS7T0prbhE+fxI+J2pMqQDZWTHBvWHSrQKOzm0osTnQzCab2X+b\n2TIzm2JmR6W8Xgx86Dvfg2dV+N5WKf3y6d472c5LcUp7ej+o2+OxqZndbmbzzWy5mT1pZien9am3\necowRiC4wKITHB4HjaeoeYI8H0+ZmNkhwG0EazJuC5sTNaZUgOysFbAxQ/sGgkG+V47jybXtBNNN\nE9z9BII7Em8D/tvMuoV9dpUj+O6eOvl2751s56Xqh0H6PutD/jYBj7t7D4IfYkuBMjNLvct13uTJ\nzBoSzLHf499dRFHjKU01eQKNpx3MrIOZvUuwHsSAIe7+VfhyosaUCpCdZfMQcp3j7mvcvdjdF4bb\nG4H/RfAP/MZYg5N6w90P9/BGke6+xd2vJlhcmH6Kfb4YC2wlOJtAqpcxTxpP33H399z9RwTTL+8A\ni83spJjDykgFyM4+BZplaG8ObHL3rTmOJ3buvoVg8dEJYdOucgTBgqaqfoVmln45+/R+9UW28/Jp\n+Jy+z/qav9eAH5nZvuF2XuQpvBjiMIKrOKfeYFPjKcUu8lSdvBxPVdx9o7v/FvgIuD1sTtSYUgGy\nswqCe8OkKyI4F7peM7PmZlaQ4aXtQMPwzxXAIWaWfiG7ImCdu1cNysUEY6xthn4Z771Tx2U7L4vD\n5/YZ+lV9Xp1jZoVmtneGl7aHz6njrF7nycxGAL8DfpoyPqpoPIV2lSeNp0A1OYDgl8cfm1ljEjam\nVIDsbDpwmJkdVtVgZgcCnYBpsUWVO/+PtLN9woF7NLAgbJpGcOrWSWl9TuL7Oaqv996pbpou23l5\nDthcTb8lHt7zKMGqy9P/IfNUw7HAByk/BKdTj/NkZsMJbjXRx90/DtsGmtmlYReNJyLlSeMp8JyZ\ndc/Q3h7Y4O7fEOQgOWOqts9NrmuP8C9nMfAIQeXcgOCy0suApnHHl4Pvfx/BqumDwu2GBFce3Ebw\nA6CqX9XFbKouEPRHgkN9bdL2dwffv5jNRQTrSYrj/q57mJ/7iXYhsqzkBfg9wUWC0i/y8/O4c7Gn\neSK4FcJqoENK2xUE1ycYmQ95Ai4AvgYuB4anPO4Cxms8Rc+TxtOOmEsJrmq6X7htwOgwD9cmcUzF\nnrQkPvjuUuzLyL9LsR8FTAQWhY8PgFlAr7R+hQSX811GcDhuJtA5w/4aEVzO923q9qXYbyc4tWwj\nwaHdlQSXHC+o7bwA/8l3lzleAJwRdz5+SJ4IfiO7KfwuCwlO9ZtLsFo/X/L0WZifyrTHdmCcxlP0\nPGk87Yi3BzAp/F4LCRbhvgycn9YvMWMq0XfDFRERkfpJa0BEREQk51SAiIiISM6pABEREZGcUwEi\nIiIiOacCRERERHJOBYiIiIjknAoQERERyTkVICIiIpJzKkBEREQk51SAiIiISM6pABEREZGcUwEi\nIiIiOacCRERiZ2YzzGyTmVWa2edmdp6ZNTOzdWa2zcw+M7Or4o5TRLJHd8MVkUQws1OAMuBxdz83\nbDsOeAD4ibtviTE8EckyHQERkURw95eBvwFnm9kwM9sLmAT8u4oPkfpHR0BEJDHMrClQATQDpgPb\n3H10vFGJSG1QASIiiWJmpwIvAV8Ah7r75ngjEpHaoCkYEUmaCmA90BI4LeZYRKSW6AiIiCSKmT0I\nvAuMAPYBurj7Z/FGJSLZpiMgIpIYZjYIOBz4I3AJ0Aq4LdagRKRW6AiIiCSCme0LvAGc5u5vh213\nAr8EznL3J+KMT0SyS0dARCR2ZvZfwHvAwcDTYdu5wDDAgYfM7K34IhSRbNMREBEREck5HQERERGR\nnFMBIiIiIjmnAkRERERyTgWIiIiI5JwKEBEREck5FSAiIiKScypAREREJOdUgIiIiEjOqQARERGR\nnFMBIiIiIjn3/wFNtOKMHH60ggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb22b260780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualization of the path\n",
    "pyplot.figure(figsize=(8,6))\n",
    "pyplot.grid(True)\n",
    "pyplot.xlabel(r'x', fontsize=18)\n",
    "pyplot.ylabel(r'y', fontsize=18)\n",
    "pyplot.title('Glider trajectory, flight time = %.2f' % T, fontsize=18)\n",
    "pyplot.plot(x,y, 'k-', lw=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grid convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's study the convergence of Euler's method for the phugoid model. In the previous lesson, when we studied the straight-line phugoid under a small perturbation, we looked at convergence by comparing the numerical solution with the exact solution. Unfortunately, most problems don't have an exact solution (that's why we compute in the first place!). But here's a neat thing: we can use numerical solutions computed on different grids to study the convergence of the method, even without an analytical solution.\n",
    "\n",
    "We need to be careful, though, and make sure that the fine-grid solution is resolving all of the features in the mathematical model. How can we know this? We'll have a look at that in a bit. Let's see how this works first.\n",
    "\n",
    "You need a sequence of numerical solutions of the same problem, each with a different number of time grid points.\n",
    "\n",
    "Let's create a NumPy array called `dt_values` that contains the  time-increment of each grid to be solved on. For each element `dt_values[i]`, we will compute the solution `u_values[i]` of the glider model using Euler's method. If we want to use five different values of $\\Delta t$, we'll have five solutions: we put them in an array ... but each one is also an array! We'll have an array of arrays. How meta is that?\n",
    "\n",
    "We have one more trick up our sleeve: `enumerate()`. To get all the numerical solutions—each with its value of $\\Delta t$—done in one fell swoop, we will loop over the elements of the array `dt_values`. Within the loop, we need to access both `dt_values[i]` and the index `i`. It turns out, `enumerate()` is a built-in Python function that will give us consecutive `index, value` pairs just like we need. \n",
    "\n",
    "Read the code below carefully, and remember: you can get a help pane on any function by entering a question mark followed by the function name. For example, add a new code cell below and type:  `?numpy.empty_like`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dt_values = numpy.array([0.1, 0.05, 0.01, 0.005, 0.001])\n",
    "\n",
    "u_values = numpy.empty_like(dt_values, dtype=numpy.ndarray)\n",
    "\n",
    "for i, dt in enumerate(dt_values):\n",
    "    \n",
    "    N = int(T/dt) + 1    # number of time-steps\n",
    "\n",
    "    # initialize the array containing the solution for each time-step\n",
    "    u = numpy.empty((N, 4))\n",
    "    u[0] = numpy.array([v0, theta0, x0, y0])\n",
    "\n",
    "    # time loop\n",
    "    for n in range(N-1):\n",
    "       \n",
    "        u[n+1] = euler_step(u[n], f, dt)   ### call euler_step() ###\n",
    "    \n",
    "    # store the value of u related to one grid\n",
    "    u_values[i] = u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In [Lesson 2](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_02_Phugoid_Oscillation.ipynb), we compared our numerical result to an analytical solution, but now we will instead compare numerical results from different grids. \n",
    "\n",
    "For each solution, we'll compute the difference relative to the finest grid. You will be tempted to call this an _\"error\"_, but be careful: the solution at the finest grid is _not the exact_ solution, it is just a reference value that we can use to estimate grid convergence.\n",
    "\n",
    "To calculate the difference between one solution `u_current` and the solution at the finest grid, `u_finest`, we'll use the $L_1$-norm, but any norm will do.\n",
    "\n",
    "There is a small problem with this, though.  The coarsest grid, where $\\Delta t = 0.1$, has 1001 grid points, while the finest grid, with $\\Delta t = 0.001$ has 100001 grid points.  How do we know which grid points correspond to the same location in two numerical solutions, in order to compare them?  \n",
    "\n",
    "If we had time grids of 10 and 100 steps, respectively, this would be relatively simple to calculate.  Each element in our 10-step grid would span ten elements in our 100-step grid. \n",
    "\n",
    "Calculating the _ratio_ of the two grid sizes will tell us how many elements in our fine-grid will span over one element in our coarser grid.\n",
    "\n",
    "Recall that we can _slice_ a NumPy array and grab a subset of values from it.  The syntax for that is\n",
    "\n",
    "```Python\n",
    "my_array[3:8]\n",
    "```\n",
    "\n",
    "An additional slicing trick that we can take advantage of is the \"slice step size.\"  We add an additional `:` to the slice range and then specify how many steps to take between elements. For example, this code\n",
    "\n",
    "```Python\n",
    "my_array[3:8:2]\n",
    "```\n",
    "\n",
    "will return the values of `my_array[3]`, `my_array[5]` and `my_array[7]`\n",
    "\n",
    "With that, we can write a function to obtain the differences between coarser and finest grids. Here we go ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_diffgrid(u_current, u_fine, dt):\n",
    "    \"\"\"Returns the difference between one grid and the fine one using L-1 norm.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    u_current : array of float\n",
    "        solution on the current grid.\n",
    "    u_finest : array of float\n",
    "        solution on the fine grid.\n",
    "    dt : float\n",
    "        time-increment on the current grid.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    diffgrid : float\n",
    "        difference computed in the L-1 norm.\n",
    "    \"\"\"\n",
    "    \n",
    "    N_current = len(u_current[:,0])\n",
    "    N_fine = len(u_fine[:,0])\n",
    "   \n",
    "    grid_size_ratio = ceil(N_fine/N_current)\n",
    "    \n",
    "    diffgrid = dt * numpy.sum( numpy.abs(\\\n",
    "            u_current[:,2]- u_fine[::grid_size_ratio,2])) \n",
    "    \n",
    "    return diffgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the function has been defined, let's compute the grid differences for each solution, relative to the fine-grid solution. Call the function `get_diffgrid()` with two solutions, one of which is always the one at the finest grid. Here's a neat Python trick: you can use negative indexing in Python!  If you have an array called `my_array` you access the _first_ element with\n",
    "\n",
    "`my_array[0]`\n",
    "\n",
    "But you can also access the _last_ element with \n",
    "\n",
    "`my_array[-1]`\n",
    "\n",
    "and the next to last element with\n",
    "\n",
    "`my_array[-2]`\n",
    "\n",
    "and so on.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dt = 0.1\n",
      "dt = 0.05\n",
      "dt = 0.01\n",
      "dt = 0.005\n",
      "dt = 0.001\n"
     ]
    }
   ],
   "source": [
    "# compute difference between one grid solution and the finest one\n",
    "diffgrid = numpy.empty_like(dt_values)\n",
    "\n",
    "for i, dt in enumerate(dt_values):\n",
    "    print('dt = {}'.format(dt))\n",
    "\n",
    "    ### call the function get_diffgrid() ###\n",
    "    diffgrid[i] = get_diffgrid(u_values[i], u_values[-1], dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time to create a plot of the results! We'll create a *log-log* plot with the Matplotlib function [`loglog()`](http://matplotlib.org/api/pyplot_api.html?highlight=loglog#matplotlib.pyplot.loglog). Remember to skip the difference of the finest-grid solution with itself, which is zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAGNCAYAAACrEY57AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYVNW19/HvamYERAQFoqJBHBMjjuCUhlwFghhocABF\nuBFRo4DBKBoVWzMYIsaJwVwVCIoyCNJIC4L3plEE8SqQiMMNoAyvgAwXRW3G7vX+UQUXmh6rq+qc\nqvp9nqcfqXN2nbPgKXvV3mufvc3dERERCaOsoAMQEREpi5KUiIiElpKUiIiElpKUiIiElpKUiIiE\nlpKUiIiEVs2gAwgbM9OcfBGRKnJ3S8R11ZMqhbun/M9DDz2UFveMxzVjuUZV3lPZthW1q+75VPkJ\n6u8Rxs9nqnw2K2qTSEpSaSo7Ozst7hmPa8Zyjaq8p7JtK2pX0fk1a9ZU6j5hF8RnM1H3re41U+Wz\nWdX7xpMlOgumGjNz/ZtIGPXv358JEyYEHYbIYcwM13Bf8uTm5lJQUBB0GCKH6N+/f9AhiByioKCA\n3NzchN5DPakS1JMSEaka9aRERL17yUhKUiIiEloa7itBw30iIlWj4T4REclISlIiKUI1KclESlIi\nIhJaqkmVoJqUiEjVqCYlIiKhk5+fT6dOnRJ6DyUpkRShmpSESX5+PkOGDGHevHkJvY+SlIiIVNnT\nTz/N6tWrE34fJSmRFBHUKtQipdm9e3dS7qMkJSIiVbZz586k3EdJSiRFqCYlYbFp0yZWrlyZlHsp\nSYmISKXt3buXa665hu3bt3PmmWdyxRVXJPR+ek6qBD0nJSJStqFDh/LEE0/QsmVLPvzwQ5o3b57Q\n56SUpEpQkhIRKd3kyZPp3bs3NWvWZMGCBVx00UWAHuYVEVSTkmB9/PHH3HTTTQA88cQTBxJUoilJ\niYhIub755htycnIoLCzkhhtu4Pbbb0/avTXcV4KG+0RE/k9xcTE9e/Zk5syZnHXWWSxevJj69esf\n0kbDfSIiEog///nPzJw5kyOPPJLp06cflqASLW5Jysx+aGZtzUyJTyQBVJOSZHvrrbe4//77AXjp\npZc4+eSTkx5DlROKmfUwsy/MbP5Bx54AVgEfAh+Z2TFxjFFERJJs3bp19O7dm+LiYh588EGuvPLK\nQOKock3KzKYDtYF73f1jM7sIWAi8B7wM9AaWufsd8Q42GVSTEpFMt2vXLi699FI++OADOnXqRH5+\nPjVq1CizfSJrUjVjeM85wLnu/r/R1wOAPUCOu28ys2nA4ngFKCIiyTV48GA++OADTjzxRF5++eVy\nE1SixVI/qrM/QZlZTeAqIM/dNwG4+1dA3fiFKCKgmpQkxwsvvMBzzz1H3bp1mT59Ok2aNAk0nliS\n1Hdm1jz65+5AE2Dy/pNm1gDYG4fYREQkiT788MMDz0CNHTuWc845J+CIYqtJPQp0Bt4C+gM7gR+6\n+z4zawqMAJq7e9c4x5oUqkmJSCbatm0b5557LmvXruWWW27h2WefrfR7w1aT+h3Qgkgt6ivg36MJ\nqiawOdqmX5ziC0Rubi7Z2dnaZE5EMkJRURF9+vRh7dq1XHDBBTz11FOVel9BQUHCh6G14kQJ6klJ\nWBUUFOiLkyTEgw8+yO9//3uaNm3K0qVLOf7446v0fq04ISIiCTFr1ix+//vfk5WVxZQpU6qcoBIt\npp6UmRmRSROdgabu3tPM2gA/Ama5e1F8w0we9aREJFOsXLmS8847jx07djBixAjuueeemK4Tqv2k\nzKwekA9kRw995+6NzKwtkckU/wSucvdv4xlosihJiUgm+P7772nXrh0rVqygR48eTJ8+nUj/o+rC\nNtz3IHAqcDtwIVAI4O7LgOOAbcBv4xWgiEToOSmJF3dn4MCBrFixglNPPZUJEybEnKASLZbZfVcT\nWV1iCUR6HvtPuPtOM7sVWATcF58QRUQknkaNGsXLL7/MEUccwYwZM2jUqFHQIZUpluG+be5+9EGv\nN7p7ixJtNrh7yzjFmFQa7hORdLZw4UI6dOjAvn37mDJlCtdcc021rxm24b49ZtairJNm1hoojj0k\nERFJhI0bN3L11Vezb98+hg4dGpcElWixJKk3gKlmdtjGImZ2CfAq8Hp1AxORQ6kmJdWxd+9err32\nWjZt2sRPf/pTRowYEXRIlRJLknoAOAH4l5mtBpqY2VIz+wp4GzgSeCiOMYqISDXdc889vPPOO7Rs\n2ZIpU6ZQs2YsUxKSL9bnpI4B/gj0AvZX3L4BpgL3u/vWuEWYZKpJiUi6mTx5Mr1796ZWrVosWLCA\n9u3bx/X6oXpO6pA3R7aKbxZ9ucXdU74WpSQlIulkxYoVXHjhhRQWFjJq1KgDq5zHU9gmThzg7sXu\n/lV0D6ka0W06RCQBVJOSqvrmm2/IycmhsLCQG264gV/96ldBh1RlVU5SZnaVmW03s/81sw4HnToW\n+MrMHo5feCIiEovi4mL69evHypUrOeuss/jrX/8a2gd2yxNLT+pGYBXQFSg46PiXRNbz62lmg6sf\nmogcTCugS1WMGDGCvLw8GjduzIwZM6hfv37QIcUklod5VwKXu/uaMs63IbKd/BnVDy/5VJMSkVQ3\nf/58OnfuTHFxMbNnz6Zr18TuQRu2mlTDshIUgLuvBI6KOSIRKZVqUlIZa9eupXfv3hQXFzN8+PCE\nJ6hEiyVJ7a1gxYkWQMpu1SEikqp27dpFr1692LZtG507d2b48OFBh1Rtsa44Md3MTi95wszOBaYR\n2cpDROJINSmpyODBg/nggw848cQTmTRpEjVq1Ag6pGqLpSbVHFgCHA+sBzYAtYEWQHNgHXBhdFp6\nylFNSkRS0QsvvMCAAQOoW7cuixYtom3btkm7d6hqUu6+CTgfGEdkCaQLgbZAfeAF4PxUTVAiYaaa\nlJTlgw8+OPCQ7rPPPpvUBJVoMS3e5O6bgQFmNpBSVpwws9Pd/dM4xSgiImXYunUrPXv2ZPfu3dx6\n663069cv6JDiqlrLIpV5Ue0nJSKScEVFRXTp0oX58+dz4YUXsmDBAurUqZP0OBI53BdTT8rMugO/\nAFoCtUqeBo4+7E0iIhJXDz30EPPnz6dZs2a8+uqrgSSoRItl4sSdwF+AfcBWYE8pzX7g7iWTV0pQ\nT0rCqqCgQDP85IC8vDy6d+9OVlYW8+fPp2PHjoHFErae1B3Ab4G/uHtpCQoz21itqEREpEwrV67k\nxhtvBODRRx8NNEElWiw9qS3AMeV1N8ws290LqhlbINSTEpEw+/7772nXrh0rVqwgJyeHV199NfCF\nY0M1BR34GGiSgOuKiEg53J2bb76ZFStWcNpppzF+/PjAE1SixZJM7gDGmNlZ5bSZFGM8IlIGPScl\nzzzzDK+88goNGjRgxowZNGrUqOI3pbhYalKvAw2Aq82skMjkiYN35DWgaRxiExGRqIULF3LXXXcB\nMG7cOE4//bCV6dJSLDWp3cAiIsmoLO3cvW51AguKalIiEjYbN27knHPOYdOmTdx1112MHDky6JAO\nkciaVCxJaqO7l7kKemXbhJWSlIiEyd69e+nYsSMLFy4kOzub+fPnU7NmTI+4JkzYJk70rUSbzjFc\nV0TKoZpUZrr77rtZuHAhP/jBD5g8eXLoElSixbLA7FuVaPOP2MIREZH9XnnlFZ566ilq1arFtGnT\nOPbYY4MOKeliWrvPInMeuxPpMTV1957RbeN/BMxy95Td9FDDfSISBitWrODCCy+ksLCQ0aNH86tf\n/SrokMoUqhUnzKwekU0Ns6OHvov+twHwPPBPM7vK3b+NS4QiIhnmm2++IScnh8LCQvr27cttt90W\ndEiBiaUm9SBwKnA7kb2kCgHcfRlwHLCNyLJJIhJHqkllhuLiYvr168fKlSv5yU9+wrPPPpv2D+yW\nJ5YK3NVAjrsvgcjw2P4T7r7TzG4lMkX9vviEKCKSOUaMGEFeXh6NGzdm+vTp1K9fP+iQAhXLFPRt\n7n70Qa8Pm26u/aRERKpu/vz5dO7cmeLiYmbPnk3Xrl2DDqlSwjYFfY+ZlfkMlJm15tAVKELBzGqZ\n2TAz+87MTgg6HhGRg61du5bevXtTXFzMQw89lDIJKtFiSVJvAFPN7OSSJ8zsEuBVIksnhc0tQAGQ\n2X1nSVmqSaWvXbt20bNnT7Zt20aXLl0YPnx40CGFRiw1qQeA94B/mdkXQBMzWwr8AGgGrAEeiluE\nceLuo4CMLkCKSDgNGjSIDz/8kJNOOomXXnqJrCxtJLFfLA/zbgTOB8YR2Sa+FnA2UBt4DrjA3TfH\nK0Aza2Fmc80sdEOIIsmkXXnT0/PPP8/zzz9P3bp1mT59Ok2aVLQTUmaJaX2NaBIaYGYDifSeALa4\ne1wTiZnlACOJbFFf3iaLxwBPAOdGD30E3OnuX8YzHhGReMjPz+fpp59m69atLFu2DIBnn32Wtm3b\nBhxZ+FS5J2Vm481snJn9m7sXu/tX0Z9E9HTuAjoSGV4sdZzOzGoD84kk3DOiP98DfzezIxIQk0gg\nVJNKD/n5+QwZMoR58+axdOlS3J1GjRrRtKl2OCpNLAOf/YgM8W2KcyylucTd11TQph/wY2BYNGkW\nA8OAHwJlPaatwpSIBOLpp59m9erVhxzbsWMHzzzzTEARhVssSWqzu/d19xVxj6aESj6w1BNYe3Ay\nc/evgE+i5wAws0vN7Bkiw4a5ZvaLOIcrklCqSaWHXbt2Vel4poslSa0ws3If1DWzv8UYTyzOAr4o\n5fgaIj0sANz9HXcf5O413P3f3T0vWQGKiOz35Zell8rr1k3JfWITLpYkNRgYa2bnldPmihjjiUVT\noLTFbHcA9c2sThJjEUkY1aRS34QJEw4b6gNo3bo1gwYNCiCi8Itldl8+kRXP3zezncAWDl1hwogk\njmTRGkYiEnrz5s3j5ptvBuC2227j888/Z9euXdStW5dBgwZphYkyxJKkWhJZQLa8yQfJ3Dp+K9Cw\nlOONgO/dfXdVL9i/f39OPPFEABo3bszZZ599oB6w/9usXut1sl9nZ2eHKh69rvzrxo0b07NnT/bt\n28d1113HmDFjQhVfVV/v//OaNWtItFgWmD1sQdlY2lTxnhOAG939sOFJM5sDnObuJ5U4/hHwrbtf\nVMV7aYFZEYmbdevW0a5dOzZu3Ejv3r3TckWJsC0w27cSbTrHcN2KlJU5ZgCtzKzV/gNmdixwGjA9\nAXGIBOLgb7GSGrZv306XLl3YuHEj2dnZjB8/Pu0SVKLFsizSWxbRw8z+ambTAcysTfRYDXf/R/xD\nLXN4cQKRFSZGmFkNM8sC/gR8DoxNQBwiIhXavXs3PXr04JNPPuGMM87gtddeo04dzeOqqliG+w7b\nPt7dG5lZW+At4J9AXLaPN7MxQBciEzHqA+uI9KhOdfe9B7XbvyzSedHzMS+LpOE+Eamu4uJibrjh\nBl555RVatGjBe++9xwknpO8OQYkc7oslSf2RyCoPvwc+AF539+bRc/WAF4GV7p6SO/MqSYlIdd17\n772MGDGCBg0a8M4773D22WcHHVJCha0mtX/7+LHu/t8cVCty953ArRy00oOIxIdqUqlh9OjRjBgx\ngpo1azJ9+vS0T1CJFkuSauLuS8o66e5biTxHlbJyc3P1C0FEqmzWrFkMHjwYgOeee44rrkjmugbJ\nV1BQQG5ubkLvEdMUdOCc6L5Sh003j24fv8Ddj4trpEmi4T4RicWSJUvo0KEDO3fu5OGHH86o3XXD\nNtyXqtvHi4gkxKpVq7jyyivZuXMnN910Ew8++GDQIaWNWHpSLYjs73Q8kYVdjwM+5tDt49vFc3fe\nZFJPSsKqoKDgwJP/Eh5btmzhoosuYtWqVXTq1InXX3+dWrVqBR1WUoWqJ5Xs7eNFRMKqsLCQq666\nilWrVtG2bVumTZuWcQkq0SrsSZnZA0R2x73a3beVOJdFArePD4J6UiJSGUVFRfTs2ZO8vDxatWrF\n4sWLadEimcuWhkfQPal/J9JD2h4Npt/+E0nYPl5EJHTcnTvvvJO8vDwaN27MnDlzMjZBJVplktQR\n7v7KQUnoTxW9wcyWVy8sESlJj0WEx+OPP86oUaOoXbs2eXl5nH766UGHlLYqk6QKzax9Fa97bCzB\niIiE3eTJk7n77rsBmDhxIpdddlnAEaW3yuwnNQl418y+AnYBzczs83LaJ3vTw7jLzc09sH+PSFjo\n8xi8BQsW0K9fpOLx2GOPce211wYcUbAKCgoS3sOvzMSJGsBtwM+Ao4D2wOIKrtvO3evGJcIk08QJ\nESnNJ598wsUXX8zXX3/NoEGDeOqppzBLyFyBlBO2BWaTvulhMilJSVjpOangbNiwgfbt27Nu3Tp6\n9OjBtGnTqFGjRtBhhUags/vM7AEz+y8zOzp66IZKXLcyGyOKiITet99+S9euXQ/ssDtp0iQlqCSq\nzHDfauABYIq7F5tZP3f/W1KiC4B6UiKy3969e+nWrRtvvvkmJ598MosWLaJZs2YVvzHDBP2clKag\ni0jGcXduvfVW3nzzTZo1a8bcuXOVoAKgKegiKULPSSXXI488wrhx46hXrx6zZ8+mdevWQYeUkTQF\nXUSkhPHjx5Obm0tWVhaTJ0/mggsuCDqkjFWZmlQW8CuqNgW9vbvXiUuESaaalEhmmzdvHl27dmXf\nvn2MGTOG2267LeiQQk9T0JNISUokcy1fvpxLL72U7777jmHDhvGnP1VYgheCnzhRUmWml98Tw3VD\nQ9vHSxjpM5lY69at4+c//znfffcdvXv35o9//GPQIYVeKLePr9RFzTa4e8u4XzgJ1JOSsNLDvImz\nfft2LrnkEj755BOys7OZO3cudeqkZMUiEIEO95lZLtDE3QdHX38BlPcmA45z95Tc+UtJSiSz7N69\nm06dOrFgwQLOPPNMFi5cSOPGjYMOK6UEnaQ2AI2Apu6+y8z2AO8SSUZl0dp9IhJ6xcXFXH/99Uye\nPJmWLVuyePFiTjjhhKDDSjmJTFKVmYJ+DlDH3XdFX29z9w7lvcHMNlY7MhE5hIb74u++++5j8uTJ\nNGjQgPz8fCWoEKowSbn7phKHbq3EdSvTRkQkMKNHj+bPf/4zNWvWZPr06Zx99tlBhySlSMjEiVSm\n4T6R9JeXl0dOTg7FxcWMHz+e/v37Bx1SSgvbFPQKmdlHibiuiEh1LVmyhN69e1NcXMzDDz+sBBVy\nFQ73mVk/yp/Nd9hbgONjjkhESqWaVPWtWrWKK6+8kp07d3LTTTfx4IMPBh2SVKAyEyfGJzwKEZEE\n27JlC126dGHr1q107tyZsWPHamfdFFCZKejbgBz+b8p5PeBR4L+BAmD/xIrmQDZwOXCnu8+Mf7iJ\np5qUSPopLCykY8eOLFmyhLZt27JgwQIaNmwYdFhpI+gp6AvdfcFBwTwF/N7dXy2l7SQzywGuBFIy\nSUFkWaTs7GwNrYikgaKiIvr06cOSJUto1aoV+fn5SlBxUlBQkPDlumJZYPZzoHVZ3Y3oqulfuHur\nOMSXdOpJSVipJlV17s6gQYMYPXo0Rx11FO+++y6nn3560GGlnbDN7msKlLdmSGMiW3qIiARq5MiR\njB49mtq1a5OXl6cElYJi6UktAPYCv3b3j0qcOwv4C1DT3bPjFWQyqSclkh4mT55M7969D/z52muv\nDTii9BW2/aTOBf4ONAA2Rn+MyMSJFsB3QLa7L41vqMmhJCWS+hYsWMAVV1zBnj17GDlyJHfddVfQ\nIaW1UA33ufuHwPnAdOBIImv7tY3+eRpwXqomKJEw035SlfPJJ5/QvXt39uzZw6BBgxg6dGjQIUk1\nVGZ232Hc/X+Aq82sBtAseniLuxfFLTIRkSrasGEDXbp04euvv6ZHjx488cQTehYqxWntvhI03CeS\nmr799lsuu+wyli9fTvv27fnP//xP6tWrF3RYGSFUw30iImGzd+9err76apYvX06bNm2YNWuWElSa\nUJISSRGqSZXO3bnlllt48803adasGXPmzKFp06ZBhyVxoiQlIintkUceYfz48dSrV4/Zs2fTunXr\noEOSOFJNqgTVpERSx/jx4/nlL39JVlYWM2fOpFu3bkGHlJFUkxIRKeHNN9/k5ptvBmDUqFFKUGlK\nmx6WIjc3V+P/Ejr6TP6fZcuW0atXL4qKihg2bBi33XZb0CFlpIKCAnJzcxN6j8ps1RHLpodPuXt5\n6/uFlob7JKy0wGzE2rVrad++PRs3bqRPnz68+OKLZGVpUChIgS6LZGbFsVzY3VPyU6MkJRJe27dv\n5+KLL+bTTz8lOzubuXPnUqdOnaDDynhB7ye1nRg2PYxnkCIiu3fvpkePHnz66aeceeaZvPbaa0pQ\nGaAyPak8d//FQa+fAt4pY9NDopse/tzdB8Q10iRRT0rCKpOH+4qLi7n++uuZPHkyLVu2ZPHixZxw\nwglBhyVRgc7uOzhBRXUjsrhsWWYS6U2JiMTFfffdx+TJk2nYsCFvvPGGElQGiWWrjh1AK3ffXsb5\nJsAad28Uh/iSTj0pkXAZPXo0d9xxBzVr1uSNN97g8sv1HThswvac1DJgmpn9uOSJ6KaHUwFt1SEi\n1ZaXl8fgwYMBeP7555WgMlAsW3UMJbLp4T/MrMxND+MVoIhEZFpN6r333qN3794UFxfzyCOP0K9f\nv6BDkgBo00MRCZ1Vq1bRrVs3du7cyYABA3jggQeCDkkCUq21+9Jx00PVpESCtWXLFtq3b8/q1avp\n3Lkzs2bNolatWkGHJeUIW03qAHcvcvdN0Z8DCcrMbqh+aCKSaQoLC+nWrRurV6/mnHPOYdq0aUpQ\nGS5Rq0I8lqDrimSsdF+7r6ioiD59+rBkyRJatWpFfn4+DRo0CDosCViFEyfMLBdo4u6Do6+/oPy1\n/AzQjmMiUmnuzpAhQ8jLy+Ooo45izpw5NG/ePOiwJAQqs+LEBqAR0NTdd5nZHmBRBddt5+514xRj\nUqkmJZJ8jz32GPfccw+1a9fmrbfe4tJLLw06JKmCoBeYbQ7Ucfe10dcb3b1FBe+psE1YKUmJJNcr\nr7xCnz59AJgyZQrXXHNNwBFJVQW9LNKm/Qkqqo6ZjTOzfyvnbbdWP7TgaD8pCaN0/EwuWLCA/v37\nAzBy5EglqBQTiv2kDntDZOuOScAId1+RkKgCpJ6UhFW6Pcz7ySefcPHFF/P1118zePBgnnzyScwS\n8mVcEizQ4b5Sgtnk7mlb0VSSEkm8DRs20L59e9atW0ePHj2YNm0aNWrUCDosiVHYnpNaYWYty2tg\nZn+LMR4RSXPffvstXbt2Zd26dbRv355JkyYpQUmZYklSg4GxZnZeOW2uiDEeESlDOtSk9u7dS69e\nvVi+fDlt2rRh1qxZ1KtXL+iwJMRiWWA2H2gAvG9mO4EtwMFbzOs5KRE5jLtzyy23MG/ePJo1a8ac\nOXNo2lS/KqR8sdSkdhN5Tqq88Uc9JyUih8jNzeXhhx+mXr16FBQUcMEFFwQdksRJ2CZO6DkpEamS\ncePGcdNNN5GVlcXMmTPp1q1b0CFJHIVt4kTfSrTpHMN1RaQcqVqTevPNNxk4cCAQ2WVXCUqqIpb9\npN6qRJt/xBaOiKSTZcuW0atXL4qKirj33nu59daUfs5fAhDLcN9H7n7Y1vHpQsN9IvGxdu1a2rVr\nx6ZNm+jTpw8vvvgiWVmJ2nhBghS2mtTXRKahl6cY2Aa85+7bY4wtEEpSItW3fft2Lr74Yj799FM6\ndOjAnDlzqFOnTtBhSYKELUkVV9zqgN3An939oSrdJEBKUhJWqbIs0u7du7niiit4++23OfPMM1m4\ncCGNGzcOOixJoEQmqViek+oFPAHMAwqATdHjzYFs4FzgAaAh0B4YamYb3P2v1Q1WRMKtuLiYfv36\n8fbbb9OyZUvmzJmjBCXVEktP6gVgjru/Wsb5HOAyd78z+ro7kOvuZ1c32GRQT0okdvfccw+PPfYY\nDRs25J133uEnP/lJ0CFJEoRtuG8NcFJZv8nNLAtY5e4/jL6uAWxz95T4OqUkJRKbUaNGMWjQIGrW\nrMkbb7zB5ZdfHnRIkiRhe06qCVBewmnMocsiFQOFMdxHRA4S5uekZs6cyeDBkflUzz//vBKUxE0s\nSWoZMM3MDpuGbmZnAVOBpQcdvhrYGlt4IhJ27733Hr1798bdeeSRR+jXr1/QIUkaiWW471zg70QW\nmd0Y/TEiEydaAN8B2e6+1MyeA/oDj7v7vXGMO2E03CdSeStXruSiiy5i69atDBgwgP/4j//QxoUZ\nKFQ1KQAzOxX4PdAFqB89XEhkhfQH3f1f0XY/Bo4A/idVnpdSkhKpnM2bN3PRRRexevVqunTpwqxZ\ns6hZM5YJw5LqQpekDrw5MimiWfTlFncviktUAVKSkrAK03NShYWFdOjQgffff59zzjmHBQsW0KBB\ng6DDkoCE7TmpA6JJaVOFDVNMbm4u2dnZofmFIBImRUVF9OnTh/fff58TTzyR/Px8JagMVVBQkPAJ\nPdXqSaUj9aREyubu3HHHHYwZM4ajjjqKRYsWcdpppwUdlgQsbFPQRSRDPfbYY4wZM4Y6deqQl5en\nBCUJpyQlkiKCfk7qlVdeYdiwYQBMnDiRSy+9NNB4JDMoSYlIhQoKCujfvz8Ajz/+ONdcc02wAUnG\nUE2qBNWkRA718ccfc/HFF/PNN98wePBgnnzyST0LJYcI3RR0i3xCuxPZJr6pu/c0szbAj4BZqTwV\nXUlK5P9s2LCBdu3asX79enr06MG0adOoUaNG0GFJyIRq4oSZ1QP+E5gO3AzsX6SrAfA88JaZNYxb\nhCICJL8mtWPHDn7+85+zfv162rdvz6RJk5SgJOliqUk9CJwK3A5cSHTxWHdfBhxHZEfe38YrQBFJ\nvr1799KrVy/+8Y9/0KZNG2bNmkW9evWCDksyUCxr960EbnD3JdHXG929xUHnmwKL3P2UuEaaJBru\nk0zn7vzyl79kwoQJNGvWjMWLF9O6deugw5IQC9VwH9Bkf4IqjbtvJTL0JyIp6OGHH2bChAnUr1+f\n/Px8JSgJVCxJao+ZtSjrpJm1JrKHlIjEUTJqUuPGjePhhx8mKyuLKVOmcP755yf8niLliSVJvQFM\nNbOTS577NYS7AAAar0lEQVQws0uAV4HXqxuYiCTX3LlzGThwIACjR4/myiuvDDgikdhqUi2A94Dj\ngS+ITJb4GPgBkRXR1wDt3H1zXCNNEtWkJBMtW7aMyy67jO+++457772XRx99NOiQJIWE8TmpY4A/\nAr2ARtHD3xDZlff+aF0qJSlJSaZZu3Yt7dq1Y9OmTVx//fVMnDiRrCwtRiOVF7okdeDNZlkcup9U\nyteilKQkrBKxn9T27du5+OKL+fTTT+nQoQNz586ldu3acb2HpL+wze47wN2L3f2r6M+BBGVmp1c/\nNBFJpN27d9O9e3c+/fRTzjzzTGbMmKEEJaGTkLX7zGyDu7eM+4WTQD0pyQTFxcX06dOHKVOm0LJl\nS9577z2OP/74oMOSFBW6nXnNrDvwC6AlUKvkaeDoasYlIgk0bNgwpkyZQsOGDZkzZ44SlIRWlZOU\nmd0J/AXYB2wF9pTSTFVXkTiLV01q1KhRjBw5kpo1azJ9+nTOOuus6gcnkiCx9KTuAO4FnnT30hIU\nZraxWlGJSELMnDmTwYMHA/DCCy9w+eWXV/AOkWDF8pzUFuCY8go3Zpbt7gXVjC0QqklJulq8eDEd\nO3Zk165d/O53v+OBBx4IOiRJE2Gb3fcx0CQB1xWRBFm5ciXdunVj165dDBgwgPvvvz/okEQqJZZk\ncgcwxsx+Uk6bSTHGIyJliHXtvs2bN9OlSxe2bdtGly5dGDt2rHbWlZRRYU3KzL4ASo5/NQSuNrNC\nIpMnDn6I14CmcYtQRGJWWFhIt27dWL16Neeccw5Tp06lZs2YJvWKBKLCmpSZ7QYWEUk+ldXO3etW\nJ7CgqCYl6aKoqIicnBxmzZrFiSeeyOLFi2nevHnQYUkaCvo5qf919w5Vuahm94kEy90ZPHgws2bN\n4qijjmLOnDlKUJKSKlOT6lvZi5lZ/egfO8UWjoiUpSo1qccee4wxY8ZQp04dZs2axWmnnZa4wEQS\nqMIk5e5vHfzazMaU0/zPZraVyDYeIhKAl19+mWHDhmFmvPjii1xyySVBhyQSs1iek9ro7qXuzBvd\nwqMT8Ft3T8lFZlWTklRWUFDAFVdcwd69e3n88ccZOnRo0CFJBgjbc1Jlim50+BIpvnZfbm5uUrbq\nFomnjz/+mO7du7N3716GDBnCr3/966BDkjRXUFBAbm5uQu9RqZ6UmY0nMg3dgGuByWU0rQmcBtR2\n97PjFWQyqSclYVXe2n0bNmygXbt2rF+/npycHKZOnUqNGjWSG6BkrKBn9wH0K/G6fxntCoFPgdti\nDUhEqmbHjh38/Oc/Z/369Vx00UW89NJLSlCSNuJak0oH6klJKtm7dy9du3Zl/vz5tGnThkWLFtG0\nqZ6ll+QKW03qt3GPQkSqzN0ZOHAg8+fP55hjjmHu3LlKUJJ2qpyk3H18IgIRkfKVnMyTm5vLhAkT\nqF+/PrNnz+aHP/xhMIGJJJBWKxdJQS+88AKPPPIIWVlZTJkyhfPPPz/okEQSoso1qXSnmpSE3dy5\nc7nyyispKiri2Wef5ZZbbgk6JMlwYatJiUhAli5dSq9evSgqKuK+++5TgpK0p55UCepJSdjk5+fz\n9NNPs379elavXs2ePXu4/vrrefHFF7UvlIRCoM9JmdkDQEfganffloggRKR0+fn5DBkyhNWrVx84\nVq9ePXr16qUEJRmhMsN9/w48B2wHMLOSD/aKSII8/fTThyQogJ07d/Lss88GFJFIclUmSR3h7q+4\n+/7dd/9U0RvMbHn1whIRgF27dlXpuEi6qUySKjSz9lW87rGxBCMih9q0aVOpx+vWTcmNr0WqrDJr\n900C3jWzr4BdQDMz+7yc9gbosXeRapo1axb/+te/DjveunVrBg0aFEBEIslXmST1EPAV8DPgKKAl\nsK6C96Tt2n4iybB8+XL69OkDQN++fdm8eTO7du2ibt26DBo0iK5duwYcoUhyJGSB2VRehFZT0CVo\nmzZt4oILLmD9+vXccMMNTJw4UTP5JNTC9jBv3zi1EZESdu7cSffu3Q9su/Hcc88dSFDaiFMyUWX3\nkzrA3d+yyP81NwLXAK2jp1YCU939RXd/K44ximQEd+emm25iyZIlnHDCCbz22muaICEZL5bhvtpA\nHtCpjCZvAle5+95qxhYIDfdJUH73u98xfPhwGjRowKJFi/jxj38cdEgilRK24b57gbbA3cCpQKPo\nz2nAPdFz98UrQJFMMG3aNIYPH46Z8corryhBiUTF0pP6F3Cduy8t4/y5wGR3bxOH+JJOPSlJtg8+\n+IDLLruMnTt38vjjjzN06NBS2xUUFJCdnZ3c4EQqIWw9qUZlJSgAd/8QaBh7SCKZ48svv+Sqq65i\n586d3HTTTfz6178OOiSRUImlJ7UeOM/dvyrjfHPgA3c/Lg7xJZ16UpIs33//PZdddhlLly7lpz/9\nKfPmzaN27dpBhyVSZWHrSc0FXjWzc0qeiA71vQrMqW5gIumsuLiYG2+8kaVLl9K6dWumT5+uBCVS\niliS1IPACcAHZrbBzD6M/mwE/hs4LtpGRMowfPhwZsyYwZFHHsns2bM5+uijK3yPnpOSTFTlJOXu\nm4DzgXFAfSKz+doC9YAXgPOjbUSkFC+99BJ/+MMfqFGjBlOnTuW0004LOiSR0KrWzrxmlgU0i77c\nctB2HilLNSlJpMWLF5Odnc2ePXsYNWoUt99+e9AhiVRbImtS2j6+BCUpSZS1a9dywQUXsHnzZm6/\n/XZGjRoVdEgicRG2iRMiUkXffvst3bp1Y/PmzVx++eU8+eSTVb6GalKSiZSkRBKsqKiIPn368NFH\nH3HqqacydepUatas8rKZIhlJw30laLhP4u3uu+9m5MiRNGnShCVLlnDyyScHHZJIXGm4TyRFvfDC\nC4wcOZKaNWsyffp0JSiRKlKSEkmQBQsWcOuttwIwduzYaq+7p5qUZKKEJCkz+30iriuSKlavXk1O\nTg779u1j6NChDBgwIOiQRFJSdZ+Takrkgd5DDgP/7e7HVCewoKgmJdX19ddf0759ez777DO6du1K\nXl4eNWrUCDoskYRJZE2qylOMoonpKSAHqFNGM/2Wl4y0b98+rr32Wj777DN+9KMf8fLLLytBiVRD\nLPNg/wP4KTAV+BLYU0qbu6oTlEiq+vWvf828efNo1qwZr7/+Oo0aNYrbtbWflGSiWJJUB+ACd19Z\nVgMzuyr2kERS05gxYxg1ahS1a9fmtdde48QTTww6JJGUF8t+Up+6++kJiichzOw44BlgE3A88FB0\nc8bS2qomJVU2f/58unTpQlFRERMnTqRv375BhySSNGF7TmqMmQ0sr4GZfRBjPIkyBpjm7rcBvwWm\nRhfHFam2zz77jKuvvpqioiLuu+8+JSiROIppdp+ZPQlcBLwPbAUOXv3cgN+4eyi2kDezJsAW4Ch3\n3xE9tgbo5+4LSmmvnpRU2rZt22jXrh2rVq0iJyeHadOmkZWVmO8/qklJWIVtdt91wCAiyei8MprF\n5be8mbUAxgNXuHus/+efCOzZn6CiNkWPH5akRCprz5499OrVi1WrVtG2bVsmTpyYsAQlkqlimTiR\nC8wARgEbKX123/vViAkAM8sBRkavX2bSM7NjgCeAc6OHPgLudPcvqxuDSFncnTvuuIOCggJatGjB\nrFmzOOKIIxJ6T/WiJBPFkqSOAXq7+76yGpjZyNhDOuAuoCORpHhKGfepDcwHPgPOiB4eB/zdzNq6\n+/fAWqC2mTU6qDfVHFgThxglQz355JM899xz1K1bl7y8PI477rigQxJJS7GMTXxEZKv48iyK4bol\nXeLuaypo0w/4MTDM3YujOwMPA34I3Abg7tuAfKAbgJn9BCgC3o5DjJKB8vPz+c1vfgPA3/72N84/\n//yk3Fdr90kmiiVJ3Q6MMrMzymkzNcZ4Dqjk7IWewNqDk5m7fwV8Ej2336+Aq81sLPBH4DrNjpBY\nrFixgt69e1NcXMzDDz/MNddcE3RIImktluekvgAaAEcD3wPbOHx23/HuHpdd3cxsAnBjaRMnzGwD\n8Jm7dyxxfBbQ0d0bxHA/5S8p1ebNm7nwwgtZs2YN1113HS+//DJmCZnQJJJSQjW7D2hJZDivvIBa\nxBZOlTUFvi3l+A6gvpnVcffdSYpF0tju3bvJyclhzZo1XHjhhYwbN04JSiQJYklS/+vuHcprYGYb\nY4ynqtTlkYRzdwYOHMi7777L8ccfz8yZM6lXr6KybPzpOSnJRLEkqRsq0aZzDNeNxVagtIeGGwHf\nx9qL6t+//4F11xo3bszZZ5994JfD/uK1XmfO65dffpmJEydyxBFHMHz4cD777DOaN28emvj0Wq+T\n/Xr/n9esWUOixVKTGk+kB/Oyu7+VkKgOvd8Eyq5JzQFOc/eTShz/CPjW3S+K4X6qSckBr732Gjk5\nOZgZM2bMoHv37kGHJBI6YVu7rx9Qi8iqDclSVtaYAbQys1b7D5jZscBpwPRkBCbpa9myZdxwQ2Tg\n4NFHH1WCEglALElqs7v3dfcVcY+mbGVl6AlEntsaYWY1oovG/gn4HBibpNgkDW3cuJGrrrqKwsJC\n+vXrxz333BN0SIcMtYhkiliS1Aoza1leAzP7W4zxHHyNMdHp7j0BN7MvzOxzM6u1v4277wUuJ/Jw\n7ifRnwZEpp8XVjcGyUw7d+7kF7/4Bf/v//0/LrnkEv76179qJp9IQGKpSZ0BPAr8zt1L3ZLDzDa6\ne7KmoceValKZzd3p3bs3U6ZM4aSTTmLJkiU0a9Ys6LBEQi1sz0nlE+mtvG9mO4lsg1HyYd6mcYhN\nJOkeeeQRpkyZQsOGDXn99deVoEQCFktPajcVP8zbzt3rViewoKgnlbmmTJnCddddR1ZWFrNnz6ZL\nly5Bh3SIAj0nJSEVtp5UmB7mTYjc3Fyys7P1CyGDvP/++/Tv3x+Av/zlL6FLUCJhVFBQkPAJPbH0\npP6touejzOwn7v6PakUWEPWkMs/69eu54IIL2LRpEwMHDuTZZ5/VRAmRKkhkTyqm7ePTmZJUZvnu\nu++49NJLWb58OR07dmTu3LnUqlWr4jeKyAFhe5gXi+hnZvlm9ln053Uz6xvvAEUSpbi4mL59+7J8\n+XLatGnDtGnTQp2g9JyUZKIq16Siu+HmAZ1KnDoF6GpmfYCros8wiYTW/fffz8yZM2ncuDGzZ8+m\nSZMmQYckIiXE0pO6F2gL3A2cSmQx10ZEliK6J3ruvngFKJIIEydO5E9/+hM1atTg1Vdf5ZRTTgk6\npAppIo9kolgmTvyLyM62S8s4fy4w2d3bxCG+pFNNKv0tXLiQn/3sZ+zZs4exY8dy6623Bh2SSEoL\nW02qUVkJCsDdP6T07TNEArdmzRp69OjBnj17GDx4cEolKNWkJBPFkqT2RlcaL5WZNQf2xR5S8HJz\nc/ULIQ3t2LGDK6+8kq1bt9KpUycef/zxoEMSSWkFBQXk5uYm9B6xDPc9R6T+NKRkjyo61PcU8Km7\n3xy3KJNIw33pqaioiKuuuoo33niD008/ncWLF3PkkUcGHZZIWgjVc1LRntIS4Hgie0rtX12iJXAs\nsI7IskjJ3G8qbpSk0tPQoUN54oknOProo1myZAmtW7cOOiSRtBGqmlQ0+ZwPjAPqE5nN1xaoB7wA\nnJ+qCUrSS35+Pp06deKUU07hiSeeoEaNGrz22mspm6A0BC2ZKJa1+3D3zcAAMxsI7F8meou7FwOY\n2enu/mmcYhSpsvz8fIYMGcLq1asPHDv66KPZsWNHgFGJSFUlZFkkM9vg7uVujBhWGu5LD506dWLe\nvHmlHp87d24AEYmkr7Ctgo6ZdQd+QaQOVXIdGQOOrmZcItWye/fuUo/v2rUryZGISHXEsizSncBf\niEwz3wrsKaVZTGsCisRLnTp1Sj1et25KbnMGaD8pyUyxJJM7gN8CDdy9pbufWPKHSPISCczgwYMP\nmyDRunVrBg0aFFBEIhKLWKagbwGOKa9wY2bZ7l5QzdgCoZpU+sjPz+eZZ55h165d1K1bl0GDBtG1\na9egwxJJO2GrSX0MNAG2ldMmpYf7tDNveujatauSkkgChXVn3h8BDwJ/cPd/ltFmo7u3iEN8Saee\nlISValISVmHrSb0ONACuNrNCIvWn4oPOG9A0DrGJiEiGi6UntRtYRCQZlaWdu6fkNCr1pEREqiZs\nPan/dfcO5TUws43lnRcREamMWCY49K1Em84xXFdEyqG1+yQTxbLA7FsHvzazdqW0+Ud1ghIREYE4\nrN2XyjP5SqOalIhI1YRqqw4REZFkUZISSRGqSUkmUpISEZHQikeS6hmHa4RKbm6uvrVK6Gi1CQmb\ngoICcnNzE3qPhGx6mMo0cUJEpGo0cUJE1LuXjKQkJSIioZWU4T4zO8Xd/5XwG8WBhvtERKomHYb7\nbkzSfUREJI3EpSdlZuOAk8pp8mN3T4ntO9STkrDSflISVmFbBb00fwAeACZQ+hYev4nTfUREJIPE\nJUm5+2oz+9DdF5R23szaxOM+IplMvSjJRHpOqgQN94mIVE3oJ06Y2eVm1tXMNKVdJEH0nJRkorgk\nFXefH73W62b2ZzM7PR7XFRGRzFbp4T4za+vuyypoUwuYCnRz93hNykgqDfeJiFRNWIb77q6ogbvv\nBQYAe2KOKAS0wKyISMVCtcCsmf2nu/+skm3fcvd/q1ZkAVFPSsJKz0lJWIWlJ9XBzP5uZveb2UVm\nVqOctl9VNzAREZGq9KS+Ab4ETose+g54F/g78F/Ah/u7IGb2orv3jX+4iaeelIhI1YRlxYk33f0a\nM2sOdDzop1P0/Ddm9jaRpNUqvmGKiEgmqkpP6lx3/7CU4634v4TVAWgJuLuXNxwYWupJSVipJiVh\nFYqaVGkJKnp8rbuPd/e+7n4ccAawMl4BiohI5krIskhmNtndr4v7hZNAPSkRkapJZE8qUUkqZTY5\nLElJSkSkakIx3FcVqZqgRMJMD5hLJtKCsCIiElraqqMEDfeJiFRNyg33iYiIxIOSlEiKUE1KMpGS\nlIiIhJZqUiWoJiUiUjWqSYmISEZSkiqFNj2UMNJnUsImVJseZgoN90lYaYFZCauUWxYplSlJiYhU\njWpSIiKSkZSkRFKEalKSiZSkREQktFSTKkE1KRGRqlFNSkREMpKSlEiKUE1KMpGSlIiIhJZqUiWo\nJiUiUjWqSYmISEZSkhJJEapJSSZSkhJJEcuXLw86BJGkU5ISSRFff/110CGIJJ2SVJoKYmgoEfeM\nxzVjuUZV3lPZthW1y5ThvKD+nmH8fKbKZ7Oq940nJak0pSRVvWuEMUmtWbOmUvcJOyWp6r0/05KU\npqCXYGb6BxERqSLtJyUiIhlHw30iIhJaNYMOIBWZ2ZNAQ2A7cDYw3t0nBRuVCJhZE2Ak8B3gwA+B\noe6+MtDARGKknlRs9rj7Te7+G+C3wPNmpn9LCYPjgUJ3H+zuQ4A3gXEBxyRygJnVMrNhZvadmZ1Q\nUXv9Yo2Bu99z0MvTgI/dvTioeET2c/d/uPsdBx36AvhBUPGIlOIWoACoX5nGaZ2kzKyFmc01s7gn\nEDM7y8ymAoOBq+N9fUl/ifx8HqQbMCqB1xepEncf5e5LKts+bZOUmeUA7wInEhmbL6vdMWY2ycw+\ni/5MM7MfHHR+sJmtjP503H/c3f/p7tcAtwNvR2sBIpWS6M9n9FxXoIG7/yVBfw1JY0n6ElWhtE1S\nwF1AR+A9oNT5+2ZWG5hPZALJGdGf74G/m9kRAO7+tLu3if78l5ll7T8XPb8E2AlcltC/jaSbhHw+\nD3pvVyK9qL4J/VtIWkrGl6jKSuckdYm7r6mgTT/gx8Awdy+O1pWGEZkRdVsZ7zkBeH7/CzM7GmgO\nrK52xJJJEvX5xMyuBi5391vd3c3sqXgFLRkjoV+iqiJtk1Qldy7sCaw9+JeFu38FfBI9V5ptgJnZ\neDN7HHgRGOLuH1UzZMkgifp8mtlZwMvANWa20cw2AgOqH7FkmIR9iSqhwlUqMv05qbOAz0o5vobI\nt4jDuPu3wHUJjElkv1g+n/8EaiUwJskA1fkSZWb7v0SNLO1NZnYpcA2RYcRcM5vp7nll3STTk1RT\n4NtSju8A6ptZHXffneSYRPbT51PCrMpfogDc/R3gHWBQZW6StsN9laSFCyXM9PmUMKvwS1Q8bpLp\nSWorkeWNSmoEfK9vqRIwfT4lzJLyJSrTk9Q/gZNKOX4SoIkQEjR9PiXMkvIlKlOSVFkZfwbQysxa\n7T9gZscSWepoejICE0GfT0lNSfkSlSlJqqxpjhOI/GOOMLMa0UVi/wR8DoxNUmwi+nxKmAX6JSpt\nk5SZjTGzL4hMhXQz+8LMPjezA9Nz3X0vcDlQROTZk0+ABkBHdy8MIm7JDPp8SgoJ9EuUduYVEZFD\nmNkYoAuRGXz1gXVEelSnRr887W93DPAEcF70/EfAne7+ZdxiUZISEZGwStvhPhERSX1KUiIiElpK\nUiIiElpKUiIiElpKUiIiElpKUiIiElpKUiIiElpKUiIiElpKUiIiElpKUiIiElpKUiIiElpKUiIh\nZGatzGyLmZ1Sibanm9nPkhGXSLIpSYmE0/1AE+APlWj7Z6BHYsMRCYaSlEjImNmPgM3ANKCnmZ1X\nTtss4GLgv5IUnkhSaasOkZAxs/HAEOAYIhsdvu3u/1ZG27bAB8Cx7r41eVGKJId6UiIhYmYXA/90\n9x3uvgp4DuhoZpeXaJdjZi8CfwO2AE+Y2ejkRyySWOpJiYSImU0C/t3d90RfHwusAv7H3Q8b9jOz\nmcB6dx+U3EhFkkM9KZGQMLNuwJv7ExSAu39FZHvuc8zsmhLtawA/BRYkNVCRJFJPSiQEohMgXnT3\n60s51xBYDWwHznD3oujxC4D3gObuvjmZ8Yoki3pSIuFwAzCptBPu/i2RqehtgJsOOtWRyDCgEpSk\nLSUpkYCZWW2go7u/UU6zMcAaYLiZ1Y0e60h0qM/MapvZgwkNVCQANYMOQES4DXAzu7OCdp8CXYhM\nTx8BNAdmRc/dCbyYsAhFAqKalEiAor2itUCzSr7Fga+BVkA3oDewFHjL3RcmJEiRAClJiYhIaKkm\nJSIioaUkJSIioaUkJSIioaUkJSIioaUkJSIioaUkJSIioaUkJSIioaUkJSIioaUkJSIioaUkJSIi\nofX/ARE7ZqFnKYAXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb22b25ff60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# log-log plot of the grid differences\n",
    "pyplot.figure(figsize=(6,6))\n",
    "pyplot.grid(True)\n",
    "pyplot.xlabel('$\\Delta t$', fontsize=18)\n",
    "pyplot.ylabel('$L_1$-norm of the grid differences', fontsize=18)\n",
    "pyplot.axis('equal')\n",
    "pyplot.loglog(dt_values[:-1], diffgrid[:-1], color='k', ls='-', lw=2, marker='o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Order of convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The order of convergence is the rate at which the numerical solution approaches the exact one as the mesh is refined. Considering that we're not comparing with an exact solution, we use 3 grid resolutions that are refined at a constant ratio $r$ to find the *observed order of convergence* ($p$), which is given by:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "p = \\frac{\\log \\left(\\frac{f_3-f_2}{f_2-f_1} \\right) }{\\log(r)}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $f_1$ is the finest mesh solution, and $f_3$ the coarsest.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The order of convergence is p = 1.014\n"
     ]
    }
   ],
   "source": [
    "r = 2\n",
    "h = 0.001\n",
    "\n",
    "dt_values2 = numpy.array([h, r*h, r**2*h])\n",
    "\n",
    "u_values2 = numpy.empty_like(dt_values2, dtype=numpy.ndarray)\n",
    "\n",
    "diffgrid2 = numpy.empty(2)\n",
    "\n",
    "for i, dt in enumerate(dt_values2):\n",
    "    \n",
    "    N = int(T/dt) + 1   # number of time-steps\n",
    "\n",
    "    # initialize the array containing the solution for each time-step\n",
    "    u = numpy.empty((N, 4))\n",
    "    u[0] = numpy.array([v0, theta0, x0, y0])\n",
    "\n",
    "    # time loop\n",
    "    for n in range(N-1):\n",
    "\n",
    "        u[n+1] = euler_step(u[n], f, dt)         ### call euler_step() ###\n",
    "    \n",
    "    # store the value of u related to one grid\n",
    "    u_values2[i] = u\n",
    "    \n",
    "\n",
    "#calculate f2 - f1\n",
    "diffgrid2[0] = get_diffgrid(u_values2[1], u_values2[0], dt_values2[1])\n",
    "\n",
    "#calculate f3 - f2\n",
    "diffgrid2[1] = get_diffgrid(u_values2[2], u_values2[1], dt_values2[2])\n",
    "\n",
    "# calculate the order of convergence\n",
    "p = (log(diffgrid2[1]) - log(diffgrid2[0])) / log(r)\n",
    "\n",
    "print('The order of convergence is p = {:.3f}'.format(p));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See how the observed order of convergence is close to 1? This means that the rate at which the grid differences decrease match the mesh-refinement ratio. We say that Euler's method is of *first order*, and this result is a consequence of that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Paper airplane challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you wanted to participate in a paper-airplane competition, and you want to use what you know about the phugoid model to improve your chances. For a given value of $L/D$ that you can obtain in your design, you want to know what is the best initial velocity and launch angle to fly the longest distance from a given height.\n",
    "\n",
    "Using the phugoid model, write a new code to analyze the flight of a paper airplane, with the following conditions:\n",
    "\n",
    "*  Assume $L/D$ of 5.0 (a value close to measurements in Feng et al. 2009)\n",
    "*  For the trim velocity, let's take an average value of 4.9 m/s.\n",
    "*  Find a combination of launch angle and velocity that gives the best distance.\n",
    "*  Think about how you will know when the flight needs to stop ... this will influence how you organize the code.\n",
    "*  How can you check if your answer is realistic?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Feng, N. B. et al. _\"On the aerodynamics of paper airplanes\"_, AIAA paper 2009-3958, 27th AIAA Applied Aerodynamics Conference, San Antonio, TX. [PDF](http://arc.aiaa.org/doi/abs/10.2514/6.2009-3958)\n",
    "\n",
    "* Simanca, S. R. and Sutherland, S. _\"Mathematical problem-solving with computers,\"_ 2002 course notes, Stony Brook University, chapter 3: [The Art of Phugoid](https://www.math.sunysb.edu/~scott/Book331/Art_Phugoid.html). (Note that there is an error in the figure: sine and cosine are switched.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "###### The cell below loads the style of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 750px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "    margin-top:0.8em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Alegreya Sans' sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 125%;\n",
       "    font-weight: 400;\n",
       "    width:600px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: 400;    \n",
       "    font-size: 45pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 30pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0,51,102);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.3em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    margin-top:16px;\n",
       "    font-size: 22pt;\n",
       "    font-weight: 600;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: rgb(102,102,0);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-size: 14pt;\n",
       "    text-align: center;\n",
       "    margin-top: 0em;\n",
       "    margin-bottom: 2em;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Nixie One', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 16pt;\n",
       "    color: rgb(163,0,0);\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.8em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "    font-family: \"PT Mono\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"],\n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "css_file = '../../styles/numericalmoocstyle.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}