{
 "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 [HyperPython](http://nbviewer.ipython.org/github/ketch/HyperPython/tree/master/) by D.I. Ketcheson, also under CC-BY."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Riding the wave"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the fourth and final lesson of Module 3, _Riding the wave: convection problems_, of the course **\"Practical Numerical Methods with Python\"** (a.k.a., [#numericalmooc](https://twitter.com/hashtag/numericalmooc)). We learned about conservation laws and the traffic-flow model in the [first lesson](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb), and then about better numerical schemes for convection in [lesson 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_02_convectionSchemes.ipynb). \n",
    "\n",
    "By then, you should have started to recognize that both mathematical models and numerical schemes work together to give us a good solution to a problem. To drive the point home, [lesson 3](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_03_aBetterModel.ipynb) deals only with an improved model—and showed you some impressive SymPy tricks!\n",
    "\n",
    "In this lesson, we'll learn about a new class of discretization schemes, known as finite-volume methods. They are the _most widely used_ methods in computational fluid dynamics, and for good reasons! Let's get started ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite-volume method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are you curious to find out why the finite-volume method (FVM) is the _most popular method_ in computational fluid dynamics? In fact, almost all of the commercial CFD software packages are based on the finite-volume discretization. Here are some reasons:\n",
    "\n",
    "* FVM discretizations are very general and have no requirement that the grid be structured, like in the finite-difference method. This makes FVM very _flexible_.\n",
    "\n",
    "* FVM gives a _conservative discretization_ automatically by using directly the conservation laws in integral form."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conservative discretization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go right back to the start of this module, where we explained conservation laws looking at a tiny control volume. To simplify the discussion, we just looked at flow in one dimension, with velocity $u$. Imagining a tiny cylindrical volume, like the one shown in Figure 1, there is flux on the left face and right face and we easily explained conservation of mass in that case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![1Dcontrolvolume](./figures/1Dcontrolvolume.png)\n",
    "#### Figure 1. Tiny control volume in the shape of a cylinder."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The law of conservation of mass says that the rate of change of mass in the control volume, plus the net rate of flow of mass across the control surfaces must be zero. The same idea works for other conserved quantities.\n",
    "\n",
    "Conservation means that any change in the quantity within a volume is due to the amount of that quantity that crosses the boundary. Sounds simple enough. (Remember that we are ignoring possible internal sources of the quantity.) The amount crossing the boundary is the flux. A general conservation law for a quantity $e$ is thus:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{\\partial}{\\partial t}\\int_{\\text{cv}}e \\, dV + \\oint_{\\text{cs}}\\vec{F}\\cdot d\\vec{A} =0\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $\\vec{F}$ is the flux, and $\\text{cv}$ denotes the control volume with control surface $\\text{cs}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Why not make the control volume itself our computational cell?**\n",
    "\n",
    "Imagine that the one-dimensional domain of interest is divided using grid points $x_i$. But instead of trying to compute local values at the grid points, like we did before, we now want to follow the time evolution of _average_ values within each one-dimensional cell of width $\\Delta x$ with center at $x_i$ (the idea is that as long as the cells are small enough, the average values will be a good representation of the quantities we are interested in). \n",
    "\n",
    "Define $e_i$ as the integral average across the little control volume on the cell with center at $x_i$ (see Figure 2).\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "e_i = \\frac{1}{\\Delta x} \\int_{x_i - \\Delta x / 2}^{x_i + \\Delta x / 2} e(x, t) \\, dx.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "If we know the flux terms at the boundaries of the control volume, which are at $x_{i-1/2}$ and $x_{i+1/2}$, the general conservation law for this small control volume gives:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\frac{\\partial}{\\partial t} e_i + \\frac{1}{\\Delta x} \\left[ F \\left( x_{i+1/2}, t \\right) - F \\left( x_{i - 1 / 2}, t \\right) \\right] = 0.\n",
    "\\end{equation}\n",
    "$$ \n",
    "\n",
    "This now just requires a time-stepping scheme, and is easy to solve *if* we can find $F$ on the control surfaces."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![finite volume](./figures/finite_volume.png)\n",
    "\n",
    "#### Figure 2. Discretizing a 1D domain into finite volumes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've seen with the traffic model that the flux can depend on the conserved quantity (in that case, the traffic density). That is generally the case, so we write that $F = F(e)$. We will need to compute, or approximate, the flux terms at the cell edges (control surfaces) from the integral averages, $e_i$.\n",
    "\n",
    "If we had a simple convection equation with $c>0$, then the flux going into the cell centered at $x_i$ from the left would be $F(e_{i-1})$ and the flux going out the cell on the right side would be $F(e_{i})$ (see Figure 2). Applying these fluxes in Equation (3) results in a scheme that is equivalent to our tried-and-tested backward-space (upwind) scheme! \n",
    "\n",
    "We know from previous lessons that the backward-space scheme is first order and the error introduces numerical diffusion. Also, remember what happened when we tried to use it with the non-linear [traffic model](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb) in the green-light problem? *It blew up!* That was because the problem contains both right-moving and left-moving waves (if you don't remember that discussion, go back and review it; it's important!).\n",
    "\n",
    "To skirt this difficulty in the green-light problem, we chose initial conditions that don't produce negative wave speeds. But that's cheating! A genuine solution would be to have a scheme that can deal with both positive and negative wave speeds. Here is where Godunov comes in."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Godunov's method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Godunov proposed a first-order method in 1959 that uses the integral form of the conservation laws, Equation (1), and a piecewise constant representation of the solution, as shown in Figure (2). Notice that representing the solution in this way is like having a myriad of little shocks at the cell boundaries (control surfaces).\n",
    "\n",
    "For each control surface, we have *two values* for the solution $e$ at a given time: the constant value to the left, $e_L$, and the constant value to the right, $e_R$. A situation where you have a conservation law with a constant initial condition, except for a single jump discontinuity is called a **Riemann problem**. \n",
    "\n",
    "The Riemann problem has an exact solution for the Euler equations (as well as for any scalar conservation law). The [shock-tube problem](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_05_Sods_Shock_Tube.ipynb), subject of your assignment for this course module, is a Riemann problem! And because it has an analytical solution, we can use it for testing the accuracy of numerical methods. \n",
    "\n",
    "But Godunov had a better idea. With the solution represented as piecewise constant (Figure 2), why not use the analytical solution of the Riemann problem at each cell boundary? Solving a Riemann problem gives all the information about the characteristic structure of the solution, including the sign of the wave speed. The full solution can then be reconstructed from the union of all the Riemann solutions at cell boundaries. *Neat idea, Godunov!*\n",
    "\n",
    "Figure 3 illustrates a Riemann problem for the Euler equations, associated to the shock tube. The space-time plot shows the characteristic lines for the left-traveling expansion wave, and the right-traveling contact discontinuity and shock."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Riemann-shocktube](./figures/Riemann-shocktube.png)\n",
    "\n",
    "#### Figure 3. The shock tube: a Riemann problem for Euler's equations. Physical space (top) and $x, t$ space (bottom)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to solve many Riemann problems from $t$ to $t + \\Delta t$, one on each cell boundary (illustrated in Figure 4). The numerical flux on $x_{i+1/2}$ is \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F_{i+1/2}= \\frac{1}{\\Delta t} \\int_{t^n}^{t^{n+1}} F\\left(e(x_{i+1/2},t) \\right)\\,dt\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "To be able to solve each Riemann problem independently, they should not interact, which imposes a limit on $\\Delta t$. Looking at Figure 4, you might conclude that we must require a CFL number of 1/2 to avoid interactions between the Riemann solutions, but the numerical flux above only depends on the state at $x_{i+1/2}$, so we're fine as long as the solution there is not affected by that at $x_{i-1/2}$—i.e., the CFL limit is really 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![many_Rieman_problems](./figures/many_Rieman_problems.png)\n",
    "\n",
    "#### Figure 4. Riemann problems on each cell boundary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Riemann solution, even though known analytically, can get quite hairy for non-linear systems of conservation laws (like the Euler equations). And we need as many Riemann solutions as there are finite-volume cell boundaries, and again at each time step! This gets really cumbersome. \n",
    "\n",
    "Godunov solved the Riemann problems exactly, but many after him proposed *approximate* Riemann solutions instead. We'll be calculating the full solution numerically, after all, so some controlled approximations can be made. You might imagine a simple approximation for the flux at a cell boundary that is just the average between the left and right values, for example: $\\frac{1}{2}\\left[F(e_L)+F(e_R)\\right]$. But that leads to a central scheme and, on its own, is unstable. Adding a term proportional to the difference between left and right states, $e_R-e_L$, supplies artificial dissipation and gives stability (see van Leer et al., 1987).\n",
    "\n",
    "One formula for the numerical flux at $x_{i+1/2}$ called the Rusanov flux, a.k.a. Lax-Friedrichs flux, is given by\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F_{i+1/2}= \\frac{1}{2} \\left[ F \\left( e_L \\right) + F \\left( e_R \\right)  \\right] - \\frac{1}{2}  \\max \\left|F'(e)\\right| \\left( e_R - e_L \\right)\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $F'(e)$ is the Jacobian of the flux function and $\\max\\left|F'(e)\\right|$ is the local propagation speed of the fastest traveling wave. The Riemann solutions at each cell boundary do not interact if $\\max|F'(e)|\\leq\\frac{\\Delta x}{\\Delta t}$, which leads to a flux formula we can now use:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F_{i+1/2}= \\frac{1}{2} \\left( F \\left( e_{i} \\right) + F \\left( e_{i+1} \\right) - \\frac{\\Delta x}{\\Delta t} \\left( e_{i+1} - e_{i} \\right) \\right)\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's try it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from matplotlib import pyplot, animation\n",
    "from IPython.display import HTML\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the font family and size to use for Matplotlib figures.\n",
    "pyplot.rcParams['font.family'] = 'serif'\n",
    "pyplot.rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply Godunov's method to the [LWR traffic model](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb). In [lesson 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_02_convectionSchemes.ipynb) we already wrote functions to set the initial conditions for a red-light problem and to compute the fluxes.  To save us from  writing this out again, we saved those functions into a Python file named `traffic.py` (found in the same directory of the course repository). Now, we can use those functions by importing them in the same way that we import NumPy or any other library. Like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from traffic import rho_red_light, flux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "You've probably noticed that we have the habit of writing detailed explanations of what a function does after defining it, in comments.  These comments are called *docstrings* and it is good practice to include them in all your functions.  It can be very useful when loading a function that you aren't familiar with (or don't remember!), because the `help` command will print them out for you.  Check it out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function rho_red_light in module traffic:\n",
      "\n",
      "rho_red_light(x, rho_max)\n",
      "    Computes the \"red light\" initial condition with shock.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    x : numpy.ndarray\n",
      "        Locations on the road as a 1D array of floats.\n",
      "    rho_max : float\n",
      "        The maximum traffic density allowed.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    rho : numpy.ndarray\n",
      "        The initial car density along the road as a 1D array of floats.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(rho_red_light)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can write some code to set up our notebook environment, and set the calculation parameters, with the functions imported above readily available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set parameters.\n",
    "nx = 100  # number of cells along the road\n",
    "L = 4.0  # length of the road\n",
    "dx = L / nx  # cell width\n",
    "nt = 30  # number of time steps to compute\n",
    "rho_max = 10.0  # maximum traffic density allowed\n",
    "u_max = 1.0  # speed limit\n",
    "\n",
    "# Get the grid-cell centers.\n",
    "# x_i is now the center of the i-th cell.\n",
    "x = numpy.linspace(0.0 + 0.5 * dx, L - 0.5 * dx, num=nx)\n",
    "\n",
    "# Compute the initial traffic density.\n",
    "rho0 = rho_red_light(x, rho_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAELCAYAAAAP/iu7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEoFJREFUeJzt3X2M5VV9x/H3d/ZB3V0UhCkgyoO77SoVRKXGiOgoSpXWWIrV1kqamnaT1tpKq9WiBVGpIEYbrQ9d0LZJrY/VPiTiQ1pGWiTCYn1AhQUXuiJSRagwC8vM7Hz7x70zs5nusjPL/c059+z7lUxu9jf3Jt98z+R+9vc753d+kZlIktSlkdIFSJLaZ9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOreydAFdO/jgg3PDhg2ly2jOjh07WLt2bekymmNfu2Ffu3HdddfdmZmji3lv82Fz+OGHs2XLltJlNGd8fJyxsbHSZTTHvnbDvnYjIv57se/1MpokqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzho0kqXOGjSSpc4aNJKlzRcMmIo6MiM9HRJasQ5LUrWJhExFnAlcD6/fxvlUR8baIuCEiro+Ir0TEs5anSknSIJQ8s3kj8ALgqn28733Ay4FTM/NJwEeAL0XESR3XJ0kakJJhc0pm3vRgb4iIjcAm4KLM/DFAZl4GbAMu7L5ESdIgFAubzJxexNvOBAK4YsHxfwdOj4h1Ay9MkjRwK0sXsA8nAjPA9gXHb6FX+/HANctdlKTh8tUfTvPWd40zPeNapEE46uBH8LFNz1jSZ2oPm8OA+zJz14Lj9/RfD93ThyJiE73Lb4yOjjI+Pt5ZgQeqiYkJ+9oB+9qNL2/fyba7o3QZzdi58/4l/53WHjZ786B/NZm5GdgMsHHjxhwbG1uOmg4o4+Pj2NfBs6/duPiay4EZ3v2yJ/O0Yw4pXc7QWzESPPaQNUv6TO1hcyewJiJWLDi7Oaj/+pMCNUkaMrv6V88ee8gajjl0bdliDlC17yDwTXo1Pm7B8eOAaeC7y16RpKEzPdN7XbXCS2ml1B42nwUSGFtw/LnAFzPz3mWvSNLQmQ+b2r/y2lV15zPzRnpzL38WEYcBRMSr6O068KaStUkaHtPZu462emXVX3lNKzZnExGX0NtB4Oj+v7/e/9XTM3Nyt7e+BjgfuCoipoB7gdMz8+tI0iLs8symuGJhk5mvX+T7poA3938kacmcsynPmJfUvOn+arTVntkUY+clNW925wAvo5Vj5yU1b27OxgUCxdh5Sc1zzqY8w0ZS0zJzbs5m1YhfeaXYeUlNm52vWTkSjIx4ZlOKYSOpaVP9CRsXB5Rl9yU1bWp6diWaZzUlGTaSmjbZP7Nxq5qy7L6kpnkZrQ52X1LTDJs62H1JTZsPG+dsSjJsJDVtctqtampg9yU1bcoFAlWw+5Ka5pxNHey+pKZNOmdTBcNGUtOmdjlnUwO7L6lpU/0tn31wWll2X1LTnLOpg92X1LS5ORtXoxVl9yU1bX7OxgUCJRk2kpo2d5+Nl9GKsvuSmuacTR3svqSmTU4bNjWw+5KaNjdns9I5m5IMG0lNc86mDnZfUtOcs6mD3ZfUtEnDpgp2X1LTpqa9z6YGho2kpvk8mzrYfUlNc86mDnZfUtO8z6YOdl9S03x4Wh0MG0lN8z6bOth9SU3zSZ11sPuSmjbl82yqYPclNW1+gYBzNiUZNpKa5pxNHey+pKY5Z1MHuy+pad7UWQe7L6lpk3Pb1ThnU5JhI6lp83M2KwpXcmAzbCQ1bW7XZ89sijJsJDXNOZs62H1JTfPhaXWw+5Ka5n02dai++xFxckRcHhHfjYhvRcQ1EfFrpeuSNBzm77NxzqakqsMmIo4F/g24EzghM08APgJ8MiJeXLA0SUNg10yyayYJYMWIYVNS1WEDnAE8Enh3Zk4DZOaHgHuAV5QsTFL9Zi+hrRiBCMOmpNrDZrr/unL2QPT+YkYAF81LelCzYeOq5/JqD5uPAzcAb46IdRExApwLPAz4UNHKJFVvdr7GpwuUt3LfbyknM++JiNOAv6E3bzMB/BR4QWZ+eW+fi4hNwCaA0dFRxsfHl6HaA8vExIR97YB9Hay7d/bObEYi7WthVYdNRGykt0Dgc8CjgZ3Ay4DPRMQrM/PyPX0uMzcDmwE2btyYY2Njy1PwAWR8fBz7Onj2dbC+f9d9MH4Fq0ZG7GthtZ9cvg04GPijzLwvM2cy8+PAlcDfRUTVYSmprLk5m9q/6Q4AtQ/BCcBtmXn/guNbgVHguOUvSdKwcM6mHrUPwY+AI/dwBnMMkMDdy1+SpGExf2bjcrTSag+b99G7z+at/SXPRMRzgV8FPpGZd5YsTlLdZvdFc/OA8qqe88jMT0fEC4E3At+JiF3ADPAm4L1Fi5NUvalp52xqUXXYAGTmF4AvlK5D0vBxzqYeDoGkZs3vIOB1tNIMG0nNmtxtbzSV5RBIapb32dRj0XM2/WfIvARYB9wCfDYzr+yqMEl6qNyIsx6LyvuIOB/4BPBi4PHA2cB4RHytv6WMJFVnanp2gYBpU9piTy5fDXwKODQzT8zMw4BT6W2MeU1EPLGrAiVpfzlnU4/FDsGjgA/PPsAMIDOvAp4DfA14Zwe1SdJD4mW0eiw2bG4DHrfwYGYmvbv8xwZYkyQNhAsE6rHYIfggcH5EHLWX3+8cUD2SNDDzN3V6alPaYlejvQc4Dbg+Iv6K3vNlbgPWA28HLu2mPEnaf5PT7o1Wi0Wd2WTmLnor0S6m9wTM/wRupfdgs7XArRHxFJ8vI6kmXkarx6KHIDOnM/Mi4AjgmcAfAx+ltyLtA8AW4N6IuKaLQiVpqXzEQD2WfCbSXxTw1f4PABGxBngKcDLw1IFVJ0kPwdycjVlT3EAue2XmfcBV/R9JqoL32dTDIZDULJ9nUw+HQFKznLOph2EjqVnO2dTDsJHULOds6uEQSGqW99nUwyGQ1CwfC10Pw0ZSs+afZ1O4EBk2ktrlnE09HAJJzfJ5NvUwbCQ1y/ts6mHYSGrW/PNsChciw0ZSu3yeTT0MG0nN8j6bejgEkprlnE09DBtJzXLOph4OgaRmOWdTD8NGUpMyc+6mTs9synMIJDVpeqZ3CW3FSDDi3mjFGTaSmjS7OGCV19CqYNhIatLsJpyr3BitCo6CpCbNztesNmyq4ChIatL8ZTS/5mrgKEhq0lzYuOVzFQwbSU3yzKYujoKkJk32Fwg4Z1MHR0FSk2bPbFZ7R2cVHAVJTfIyWl0cBUlNmvSmzqoYNpKaNLvjs2c2dXAUJDVpatqbOmviKEhqknM2dXEUJDVpbs7G1WhVGIpRiIizIuLKiLguIrZFxJaIOLt0XZLqNT9n4wKBGlQfNhFxDvAm4BWZ+TRgI7AVOK1oYZKqNuVGnFVZWbqABxMRxwIXAc/KzNsAMnMqIl4HPKZgaZIq55xNXaoOG+Bs4H8z89rdD2bm7cDtZUqSNAwmpw2bmtQeNs8Ebo2Is4DXAqPAXcBlmfmRvX0oIjYBmwBGR0cZHx9fhlIPLBMTE/a1A/Z1cG7cNgnAHbffxsRRk/a1sNrD5nHAscDrgDOBHwFnAR+LiCMz88I9fSgzNwObATZu3JhjY2PLUuyBZHx8HPs6ePZ1cL656ybYupX1xx3DutU/tK+F1X5++XBgLfD6zLwjM2cy81PAPwPnRsSasuVJqpVzNnWpfRTu7b9+fcHx/wLWAMcvbzmShsWkYVOV2kfhhv7rwjp37eW4JAEwNe19NjWp/cv6X/uvJy44/iTgfuDby1uOpGHh82zqUvsofAK4Fnh7RKwDiIhTgZcCF2bmjpLFSaqXczZ1qXo1WmbuiogXAhcD346IncADwB9k5qVlq5NUM+ds6lJ12ABk5l3A75auQ9JwcW+0uhj5kprk82zq4ihIapJzNnVxFCQ1yefZ1MVRkNSk+TMb52xqYNhIatLsAgHnbOrgKEhqknM2dXEUJDXJ59nUxVGQ1KT57Wqcs6mBYSOpSfM3dfo1VwNHQVKTnLOpi6MgqUmGTV0cBUlNmnS7mqo4CpKaNDdn4wKBKhg2kprkZbS6OAqSmjMzk0zP9M5sVo54ZlMDw0ZSc6Zm5udrIgybGhg2kprjg9PqY9hIas7sg9N8vEA9HAlJzXFxQH1Wli6gazumkvdfcXPpMpqzbdsk3077Omj2dTDuuX8K8B6bmjQfNhNTcMkXbixdRptusq+dsK8Dc9DDm/+KGxrNj8TaVfD7Y+tLl9Gc7du3c/TRR5cuozn2dXAi4PTjjyhdhvqaD5t1q4I/feETSpfRnPHxOxgbs6+DZl/VKi9oSpI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOmfYSJI6Z9hIkjpn2EiSOjd0YRMR/xERGRHHlq5FkrQ4QxU2EXEW8KzSdUiSlmZowiYiVgPvAD5XuhZJ0tIMTdgArwa2ANeWLkSStDRDETYR8Wjg9cC5pWuRJC3dUIQNcB7w95l5a+lCJElLt7J0AfsSERuAlwFPXMJnNgGb+v98ICKu76K2A9xhwJ2li2iQfe2Gfe3GxsW+sfqwAd4JXJSZP13sBzJzM7AZICK2ZObJXRV3oLKv3bCv3bCv3YiILYt9b9VhExGnAk8CXl66FknS/qs6bIAXACuAayNi9tgR/dfPRcQkcG5muhxakipWddhk5nn0FgfMiYi3AOcDZyxywcDmwVcm7GtX7Gs37Gs3Ft3XyMwuCxm43cLmOFenSdJwGJqwiYgzgL+gdxntcOC7wGRmnlS0MEnSPg1N2EiShtew3NSpCkTEkRHx+YjwfyiSlrQLf5NhExE/ExEfjYgb+z+fjojHlq5rmEXEmcDVwPrStbQkIk6KiEsj4rqI+EZEfCci3hsRo6VrG2YRsT4i3tXv63URsbX/xfhLpWtrxVJ34W8ubPq7Q38JWA38PHA8sAO4IiLWlaxtyL2R3lL0q0oX0piPA48Gnp2ZT6bX49OBqyLiEUUrG24vAn4deHlmPg14Ar3/LP1LRDynaGUN2J9d+JsLG+C3gBOBN2TmdGbuAt4APB74vaKVDbdTMvOm0kU06g2ZuQMgM38AXAL8LHBG0aqG2w+At2TmzQCZOUNvgdEI8JKShTViybvwtxg2ZwHbM3Pb7IHMvAP4Tv932g+ZOV26hkadOPuFuJvb+6+HLHcxrcjMz2bmZQsOP7L/+uPlrqcl+7sLf4thcyJwyx6O3wKcsMy1SA8qMyf3cPjngASuXOZymhURRwHvB77Wf9X+269d+FsMm8OAe/dw/B5gjdfBVbOIWAG8CvhwZm4tXc+w6y8UuBm4jd7WV7+SmfcULmto7bYL/4VL/WyLYbM3se+3SMX9OTANnFO6kBZk5vcycwPwKGAr8I2IWPQKKv0/S96Ff1aLYXMncNAejh8E3JeZ9y9zPdKiRMRv0/tf44syc6J0PS3pn82cA/wP8IHC5Qyl3Xbh/+D+fL7qjTj30zfpLXNc6DjgW8tci7QoEXE28CfA8zLzR6XrGXb9y+U7c7ctUjIzI+JbwEsj4mGZ+UC5CofSQ9qFv8Uzm88Ax+x+R2tEHE7vSZ//WKgmaa8i4pX0luc/v79ykoj45f4TZ7V/LgeesYfjx9Kbv93Twgw9iMw8LzPXZ+ZJsz/Ah/q/PqN/bK/33bQYNn9L7wzm4ohYGREjwEX0VqPt1+mf1JWI+E3gUnp/t8+PiFf2w+fFwGNK1taACyLiUIDoeQ3wC8B7dz/j0fJociPO/pnMe4CT6S0hvR54bWZ+v2hhQywiLqF3Gn00vfs/vtH/1dP3snxXixARd7H3+2kuyMy3LGM5zYiIU4DfoRcu08DDgZ/Qm6/5B8PmodmfXfibDBtJUl1avIwmSaqMYSNJ6pxhI0nqnGEjSeqcYSNJ6pxhI0nqnGEjSeqcYSNJ6pxhI0nqnGEjSeqcYSNJ6pxhIxUUERsiYioiLlhw/IMRcW9EnFyqNmmQDBupoMy8GbgMOCciDgOIiPOAVwFnZuaWkvVJg+Kuz1JhEXEE8D1629/fAGwGfiMzP1m0MGmAWnwstDRUMvOOiPhLeo+FXgn8oUGj1ngZTarDTcDDgKsz8/2li5EGzbCRCouI5wF/DVwNnBIRTy5ckjRwho1UUEQ8FfgneosExoDt9B63KzXFsJEKiYgNwOXAF4HXZOYkcAFwRkQ8u2hx0oC5Gk0qoL8C7Sv0zmR+MTMf6B9fAVwP3J2ZzyxYojRQho0kqXNeRpMkdc6wkSR1zrCRJHXOsJEkdc6wkSR1zrCRJHXOsJEkdc6wkSR1zrCRJHXu/wBXwGkcXFfAcQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the initial car density on the road.\n",
    "fig = pyplot.figure(figsize=(6.0, 4.0))\n",
    "pyplot.xlabel(r'$x$')\n",
    "pyplot.ylabel(r'$\\rho$')\n",
    "pyplot.grid()\n",
    "line = pyplot.plot(x, rho0,\n",
    "                   color='C0', linestyle='-', linewidth=2)[0]\n",
    "pyplot.xlim(0.0, L)\n",
    "pyplot.ylim(4.0, 11.0)\n",
    "pyplot.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is a new function for applying Godunov's method with Lax-Friedrichs fluxes. Study it carefully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def godunov(rho0, nt, dt, dx, bc_values, *args):\n",
    "    \"\"\"\n",
    "    Computes and returns the history of the traffic density\n",
    "    on the road using a Godunov scheme\n",
    "    with a Lax-Friedrichs flux.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    rho0 : numpy.ndarray\n",
    "        The initial traffic density along the road\n",
    "        as a 1D array of floats.\n",
    "    nt : integer\n",
    "        The number of time steps to compute.\n",
    "    dt : float\n",
    "        The time-step size to integrate.\n",
    "    dx : float\n",
    "        The distance between two consecutive locations.\n",
    "    bc_values : tuple or list\n",
    "        The value of the density at the first and last locations\n",
    "        as a tuple or list of two floats.\n",
    "    args : list\n",
    "        Positional arguments to be passed to the flux function.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    rho_hist : list of numpy.ndarray objects\n",
    "        The history of the car density along the road\n",
    "        as a list of 1D arrays of floats.\n",
    "    \"\"\"\n",
    "    rho_hist = [rho0.copy()]\n",
    "    rho = rho0.copy()\n",
    "    for n in range(nt):\n",
    "        rhoL = rho[:-1]  # i-th value at index i-1/2\n",
    "        rhoR = rho[1:]  # i+1-th value at index i-1/2\n",
    "        # Compute the flux at cell boundaries.\n",
    "        F = 0.5 * (flux(rhoL, *args) + flux(rhoR, *args) -\n",
    "                   dx / dt * (rhoR - rhoL))\n",
    "        # Advance in time.\n",
    "        rho[1:-1] = rho[1:-1] - dt / dx * (F[1:] - F[:-1])\n",
    "        # Apply boundary conditions.\n",
    "        rho[0], rho[-1] = bc_values\n",
    "        # Record the time-step solution.\n",
    "        rho_hist.append(rho.copy())\n",
    "    return rho_hist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can run using a CFL of $1$, to start with, but you should experiment with different values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set time-step size based on CFL limit.\n",
    "sigma = 1.0\n",
    "dt = sigma * dx / u_max  # time-step size\n",
    "\n",
    "# Compute the traffic density at all time steps.\n",
    "rho_hist = godunov(rho0, nt, dt, dx, (rho0[0], rho0[-1]),\n",
    "                   u_max, rho_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_plot(n, rho_hist):\n",
    "    \"\"\"\n",
    "    Update the line y-data of the Matplotlib figure.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n : integer\n",
    "        The time-step index.\n",
    "    rho_hist : list of numpy.ndarray objects\n",
    "        The history of the numerical solution.\n",
    "    \"\"\"\n",
    "    fig.suptitle('Time step {:0>2}'.format(n))\n",
    "    line.set_ydata(rho_hist[n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video width=\"432\" height=\"288\" controls autoplay loop>\n",
       "  <source type=\"video/mp4\" src=\"data:video/mp4;base64,AAAAHGZ0eXBNNFYgAAACAGlzb21pc28yYXZjMQAAAAhmcmVlAAA1OG1kYXQAAAKuBgX//6rcRem9\n",
       "5tlIt5Ys2CDZI+7veDI2NCAtIGNvcmUgMTQ4IHIyNjQzIDVjNjU3MDQgLSBILjI2NC9NUEVHLTQg\n",
       "QVZDIGNvZGVjIC0gQ29weWxlZnQgMjAwMy0yMDE1IC0gaHR0cDovL3d3dy52aWRlb2xhbi5vcmcv\n",
       "eDI2NC5odG1sIC0gb3B0aW9uczogY2FiYWM9MSByZWY9MyBkZWJsb2NrPTE6MDowIGFuYWx5c2U9\n",
       "MHgzOjB4MTEzIG1lPWhleCBzdWJtZT03IHBzeT0xIHBzeV9yZD0xLjAwOjAuMDAgbWl4ZWRfcmVm\n",
       "PTEgbWVfcmFuZ2U9MTYgY2hyb21hX21lPTEgdHJlbGxpcz0xIDh4OGRjdD0xIGNxbT0wIGRlYWR6\n",
       "b25lPTIxLDExIGZhc3RfcHNraXA9MSBjaHJvbWFfcXBfb2Zmc2V0PS0yIHRocmVhZHM9OSBsb29r\n",
       "YWhlYWRfdGhyZWFkcz0xIHNsaWNlZF90aHJlYWRzPTAgbnI9MCBkZWNpbWF0ZT0xIGludGVybGFj\n",
       "ZWQ9MCBibHVyYXlfY29tcGF0PTAgY29uc3RyYWluZWRfaW50cmE9MCBiZnJhbWVzPTMgYl9weXJh\n",
       "bWlkPTIgYl9hZGFwdD0xIGJfYmlhcz0wIGRpcmVjdD0xIHdlaWdodGI9MSBvcGVuX2dvcD0wIHdl\n",
       "aWdodHA9MiBrZXlpbnQ9MjUwIGtleWludF9taW49MTAgc2NlbmVjdXQ9NDAgaW50cmFfcmVmcmVz\n",
       "aD0wIHJjX2xvb2thaGVhZD00MCByYz1jcmYgbWJ0cmVlPTEgY3JmPTIzLjAgcWNvbXA9MC42MCBx\n",
       "cG1pbj0wIHFwbWF4PTY5IHFwc3RlcD00IGlwX3JhdGlvPTEuNDAgYXE9MToxLjAwAIAAAA06ZYiE\n",
       "ABH//veIHzLLafk613IR560urR9Q7mffS/DfcIgDmIK8EUV8JDFX7eGduuOz73/AIdcyIgb3io39\n",
       "xyM7hq1HkbIOtWp3FhaPByIRoG4HjU4LmQqJ5Z21ubyVwpfHX252xjtdLlZUbcH8lD/+iMTkgRlf\n",
       "1PScYBJIBRXse15wHlBxdu/CfrdIXZ7F1p4oelRTh8xgEGNbLn9W2vdyQvyQTbdy6v/wMm5hqbcV\n",
       "txyHWni05C2jVkBE1AuVX/jntNdezff0w8CKGf6mHv0GSAtJlMVdV7CMd073XzAuq0qN841Ctrib\n",
       "QHfVATh9hj7qqiTVfR4Rubg+HBwxCWNLvwg+nFCWTn4TzNzUz9bjEdYHcHHTn8nQbj1nRib3AvoQ\n",
       "ND2wLAGW3zC80k9YkgqE3yo3ZMrnLwI5qFqNiGxv9nULUh88CUfn3zMDGPHgvsEFgt+LjdKmc+9X\n",
       "xOErWdjAN/hpvEm9073xzLwuvVDn69nswDTQ4nvIqxrLlmWFLg1fwPpC40bFKhI1N4Exv4hUsEgY\n",
       "wxnmZPJkGJkyQbxDOCZGcZIpmo67JjTi/xz10r/9HiqQVYjFJ1NRVNPwUHP6jtSYIrytx4nMcXig\n",
       "YkukyfWwo539/H22N87LJMru7iJ6LJwN8AXBp3rOPPtwxm5iTb/6FnXV/+u+7QUTvVzktLJQ6ETY\n",
       "AI4vIncE1U+iYY7f8LTelU19cI68Ktg6TaGv8Ah/YxRxSwLalIz1MLb2NNUe6CigS4GcP8Ti2oLI\n",
       "f3j0zvb+5mERb3ju/Hudnrfcf5tzbw/+hTgCLL8YjC9bFmEJdc8zqrrPny+EhZYgMoOHMCggpw2P\n",
       "PXiu1kzR64rgwnWDS+eAum6Wf//hhZmMnWtCqCk4NmTpP04GfNZguaKFkeFepKkRhr74j8APmb4Z\n",
       "2fbVNK7mPaZLvgAV4kvlilymX49aDZMmgsTND6NyGHVykE5nFGuUXyvEssQ3Iexrl4WcTvH8D5DV\n",
       "HRuGuhNOLRvRJr5gNUJvmZX/BTYBWKr7+ZfGmfxsjzunxQvRDzQGQ0NvJwsd9Igc8jklVsOFx+qf\n",
       "nwrbkb0UuQHD9H/PxJGujymFwcWIUQgwC3V9/nbfwLLZFa6KP7inD6s9ZoLtwuMhOaS5MqBGGO1F\n",
       "SIn/F/UYw7zqCTYF/C4KghE2/piBvvwXM+gvYTCSLtuDkiufwrcu44ckrSjjSIe9MVf1icksiSc0\n",
       "WPxbVdc+oFjvl42bPcozuq2Gls11ynQZyJBpGrT/PIuA7c1JdnBFpbfMPXX1ba41KyZQjY67Muvx\n",
       "iH3QGZzPKLGrT1SQOAHzrjT8/KN0LJulWZjUAUuMpPMXz4eNLeB0YY5VM/8Bv+F3SOOGJEgg9/Gr\n",
       "H4Gir2efrChioSI/D4PzHS5KoGScDFjEQgoXNHO4kjlhvlNvcHA90ngBnZukbs1JoK9h3I5ipmRZ\n",
       "IeFZLA3S0k9mVUEE5Ff+tEiR699iSIjTQ/MXrOFMRPdKkxQyF23Wxm1eSYxKnWWWbKk4Cyz6CpgY\n",
       "4ETIQqeoI3oOiWvJSqTvog4kWcdd0ZndifytkePd+CVLXdg2GUrJjaS9zlYytTpEfis2952D+lQI\n",
       "3PL5oi8v6MrpWj0K4aQZMbVx24fageVtdlKNexStRD1CUC6pKb56LLqIwiq1b2a8Ex1cRyqZiCPc\n",
       "0hgyD76UgM8eowZBmmYMcD7BqCaqRjn9jBTt+OX2jFmlNh0xKxd6xDQLP+nynzfnQIY3Yfq3oA5B\n",
       "21vCufmpbTDZ3evodxMAlc11ZgTQ75iJgaHTCxQAYUM2zLdIociEgk4WO9PKhqx/Ko58TIxug9gA\n",
       "PfBEiFl7vuOhQXPkbq84b+femcrPwtbMKVUfE9C9fYB8tbJ6d4hMRIbczRXZu4ncfP/m8SuDKo4I\n",
       "Gw7AdaB/6VB+MYkiubgCP28sJQv2XElpNF5sVnYqZ+tNJ7Zj3w1KVgZWC8r/MKlVaupA8lcnPKjQ\n",
       "lvZi1l7bCmkyxZr5HOIeY6G40zyBuVKytWHLmhcC4kNMyn2/np1kMq5kFfKiHUk9Mq/ZfLw+ftKG\n",
       "Lv0xRpaOUKfz21eXHtO/2q0Mn3kzcimyOwF2W4v6D+rEOXH4IVMJXKZ+xfftT+MwMjsfdJdHEuLD\n",
       "heij1EFqnBAjvLjAuwpvlFVjHeCHPH0/AgAq65c6fH/3G+c74kPRHzKmh7MbtRvGkwPllVpUANoG\n",
       "g7px1iL5RwwHqA3ny4hO9dItmYObVfTNedVZ3veP+dY4AZ4e/Ij6ryh4C5+D/hBNeq8RPFushOIO\n",
       "eYVqnnoUu1B6h9bJCEfxb2ho/k8hQ9LDD26MkZ6AU384KplOEOlA1x1UhMotH/Xyg1y/m+22Nlh6\n",
       "0D44gqh/5DfzRn86BRS2O5WDoZg9EXdU4XNX4Zq3v55Rq9HYUCMYIeP7pXIURkl8aNB0So3FNPGB\n",
       "I6xvxLRZxrU/sKB+73WxlcqB5IdP9nVWSPoLnhgRINeUgQSB+4OVUXR+AQNs8VtypE3UK+qjGPWo\n",
       "TM4RmCZJ0/iMHr0nBb0oY1YC8zrsR0TJY4VrXOWhcBsbK04tZUjZMUsNRyU2RmxTrQ7HBho0GL+T\n",
       "+B588/IudHtxB6AJMp3pDYGSNFcIzsze4Xf7dVs1vhYx/OGmXeT7iwe/pPQoFZyDkswOyPFtE2t0\n",
       "NkbReAYALYwKMtvkz/t95DUksg4ibcgr6kIOaw2bzH88zwJIO6rDvhr4VRfegYCStyV3tCfECabJ\n",
       "6e1vWDBxI58uckAH/jSQLMCYrU3Crob8kU/W5sQuNkjVCNCeANyaz78t9mONnkuJQqyvUvZP4pxE\n",
       "1difOhbh/lpzc8uFGKUB34lpzXmprFhOX0bM4alR42n/y6l0VYHlHXaro4MyC6zVVF4Iod4ORQPb\n",
       "E/vCjQ+b6hkPDLyftFCHVC4jzn6bRJQi2frQNA0B+OkRHtMJT7gRLr6CdV6q1X36evZBLMjmU00Z\n",
       "i6mjwMHn1MmyHySSKQ+6vhXaNkxIb6fSS0ct3Vw7zOWz2Ejy3VbhUNwl/8//lphDr8pXd+HXiouh\n",
       "4OvvdlPg/PG5Cq06hls/EkYu82tvWM53ec2lw8bHuNCzAK9Q9xXhTAifhJ3/ozZeMYQbn5d/5fuq\n",
       "/wzK4YfaTKddv8U5hCnGGu1eedNjSaHX32jCUZZODtDgwSe5/Nn+68K7jT1qHUoQQwFNE8G87nIW\n",
       "JXPkKWk7gIXohmOe+9OyitHCT+Jj58W68KNLZ4KH1fX6HrK4uXh0QtdxnXb2fgjFGAZvVg2Vabrw\n",
       "Yu+KfqPb+RN/mBSjgkFywQldkZEOHpTzNgffBRqEkJ4Gmke6rLmmS4YONplLuwaZMzgXUhO2rwAA\n",
       "YPmz3eY3rM0EegSRa7vqp7JrOXI26ia33qOApXwadX8eg0aETMAzLV6qB7j3awroLzyssOIrBtLj\n",
       "bN2/BRhPO5D9br2qkyRt1d40vCH0DZYwhGXS3Elt1z9x4Ph8G+4JYsFwvqIIFF+Be/pXtWJD+b7g\n",
       "KbgireBllBtEMy17yXvPKdi1MTD6Ctq+Cp60zV9nTe2uSvOWJ4MQaZw6oRdsaQptwLfjVZeKU2id\n",
       "Chyfz9taiQbZH4ZAgsTb4ARaGVJ6gj7lWdLxxVlL6oz093+SleeuU18T7vEm8UIlsEy0kXusXMwW\n",
       "eoBITarxZegVBGKNToKtPz7mrE5T8nvkIuQW2Bod/DRpNCbCUmyHNvcYKYI5zRbCEdNMugRc9j0X\n",
       "+gUHgZq1BjK+/g5dTcKng1QsuR7JHSJLKQWNBQJKybE1/iuY1k/H4wWf4Wkbo8vPzjH4bJPe00r9\n",
       "joUL7NLCncSV4YOCZtsKlzEic4Z9DKyuI9JuWBkFF1SNbHB/rBfCSDlNHcL8G14TzjejdXnRX5Yh\n",
       "xWzV5YO7HUF4R3VVB6/OFadGrVEX5jPk8d/z6KVlvCVsBtUaNqdVS8qgGLP7TIXcGOQStDXSVO5G\n",
       "v34/URrdkJ3/mppSJ4D3wA2wvAfxnelT9T9SdgV7C2ajZjpwdsCjtHCyOW+G0EC9rRER9BLZubq0\n",
       "V3dprMk8xEID8twGKvi0t7AxQ1Z8C9MBZIC+xJLw1u7R5+N1oJt0+PsOb64MeGwZTcVTvNvkq2RP\n",
       "MlFSOQ1/tDzd7AYlwzwFiAd8FnWFiKLVZ+TRhtXo3EPRqDXk7E1yK//03L1LPPM3JH4OoPOeDW1A\n",
       "zZjrJrRAYZk3wucN2fYlKdJzfX6mibHRyZWUbf6zsidZGV7PjZGDuJSc10c6rRWH0w6jjXyBW8rk\n",
       "8/+s0oHuBDrLzAYcmG0nnnTIWNrOZx+X/Hwyur+mm+aBOpBFqtJg1GMy0jq6YL5PpFGKpv3M4D+6\n",
       "sXUq52UWmPPqHKoXDa1rV/NSxDDRPpWBhum8vODPl4AEX579rar5pzg0GDp5hIJrx779oYnngSZy\n",
       "YSjigkkYDy21KeZPVhcbzHBfa886ch9b3y/cPiXOmmjqklHiXKlv/Lcv/jxc3KPv11nZLc9Xh0/+\n",
       "fcQ9uNxBzuHjEFsuAeXapQAJ22EAAAKsQZokbEEf/rUrxcrIAHaXn0d9hZXVFiIAKQfzLe68I+BI\n",
       "7jnIXFsN94XnFqRqu7s2olipqifX+pXvE0Qvo04r658uMqFRZBvewAXDbNl83gSIoQkWPztTpfQh\n",
       "F95RTAsGj3FHhoFFh/w98wSfAEqX3O7IgFjOd5IcNTrHUBGD0WAfbDcZsJjqj2QaLRuSpVM3oWmJ\n",
       "6VxtL2x5PCKjTjSk9/e9EC/CEwUtRNVa+AcupsQWk0KI2Xl9GXOYmP4f5+VmfbBxUQTd2szw5TRS\n",
       "sRqVpPZrRmR9nkqoN9+JM9eJug7YRHcx6XUWtlONN7O7HNPa923farlPpHM5Rymz0z3SH3mO3fLH\n",
       "/juVCZE9BEwCB1umyQRcJeBTtAEvHWL9n76IEzX3DnN+kP/q2EFmJXILuSFlqE+PTDvoRPap03Qo\n",
       "cfVMCmgWBd7MjTdTgkNKQULMp5mKBH2+dep+ADwv+ccb7GnxJ0v+fYTpIoL8GIodbILLe2fNXOsU\n",
       "Df19KspnA5IlarchtW2rrivsYUMiDGIibBR1cdgP+xLXFDTFciN1oNuNov6nJBszvECbK4kZ27o2\n",
       "k9NOYFxXaLX4BX+9UocICGzymiF836yT6ZYoOflUo/3ITcKS51SAtf+3hoXNxm5UfStiKF8IHq4F\n",
       "KXi6rm/S7FAd/FEmn80Dy8IKDKVQFHSl0oxtAeVRq+V5SeD2EErMSrJuIVFRF/wzyQM7FFkugxeh\n",
       "YdD/RW5eEzCfO3J9VIedtn9noXFRFwsrtMClPWBK0xzH2GzOkBiA425EYlorBTIbNWpnowYZHQWw\n",
       "7eTZfBst2fK91yAxDtgLZjnpMce5gLK9WNBoWe4rRg853I5lviNotomhUrk3rwKnjsMUSGpCUD6h\n",
       "JLkE/Oq9IbYuz/LpZIJDtgLRTyFumn2fAAABeEGeQniHfwWXItAA/pC1Xc5ZB3dYijBs72VldVvi\n",
       "dMOjV6qIwksFFOfcuddASkJNZpm8Vd2n+OfGFqgz4KBaNVXMuxkSloA8bof5ya+ZGoj3u66tWARG\n",
       "ilojysprSx3mBiyOE7lpwC3j1catZS0EjjGw576tSr6lgK1Txeusn1NidHfH5P7AZCo2SrY7E5wZ\n",
       "cFr5LC9xb2gP4KWIEM2uY0+0H5Hd2qmZvZx1l2yDsuK9CKpwJdeOq4A/eaEn+0PUKixMPWjKv2Tl\n",
       "4v5EuKFuE8jYVl3x4ibeFkB3HvvcFz13rTlASiqYqee+zZZODSdowZSxZlamv6lEGFs09XZ3Nrnc\n",
       "AI+LTfGLM5LlBKoF6f+yIwrAMLvnT+pNGZD+ALf8wp4W+KFQZczaGvr71nL9Tuw3yuQP4fpEGn02\n",
       "LH1hOu10iDi81KR4z2x1yxBNWGkGkaAsh0n8UWuCfU5O48HW+AyTG+hy+zXWOqz0xItHbA5jaH5U\n",
       "EQGjrW0AAAEYAZ5hdEN/BqzgjjhxABsoLMCjw4vvN+0qtx0H0xghlxwG9hOto21NDBerA017Lnpc\n",
       "uFAMxfoxwGHvGeyJhESk8rXEfZgYehhfAqLkvLAFP0/Wgsj1MWU6Oiy8m8yuTobAoHNa6r4Egaqe\n",
       "3zj6PgAra5FT1FVFxfbAn0AzaQH1TOqoqQVUkp+wKEgP07gCUCLE2LNTROZ1aMHysvS/o4iOIug+\n",
       "odiGk0K1uFSA6FrTq37lMvZZkz2ASqX4dLefMYak/oVy8XGIVkNz4X6DV4Fsx005lpPqP/J4gB86\n",
       "wKugGAbkEHfvJGjhpKx68ZAtMs4i17+2yPNenBcMIIyybwT/nLXU9s26bw+E0MGLJBYgNWuY26GX\n",
       "jIiqgAAAAW4BnmNqQ38HR8hAAIg4llU+OP0giDR+Jb5IfOkiQW4qw6PX6uMK7k7Usp92z5NiVCew\n",
       "JLtJJUx7mPa6HUhcvQ8ZjJiOYmLf62rsrFD2A/HGqOoaPXUGfoABRy0OBrLR4H3OxV4fGthjPWIj\n",
       "+5WMT7OkokT2ZdEhWweFH1/IMGplQ8u1gzrMzlvf9aBFPDXd88MIBZJmUg50MSUXGte/q9VoKxbG\n",
       "l5Jyv0OGG/ndd/8P/2IgLpVXx68MNgb8/wGjTVROQc7U+7q3lPdSUURQR3vGPGSFhlVbHZpOUYJ7\n",
       "0WCJRQAKpGUs0HTEgRV5bzGyx1LYFcHhNKd2FgVQ8Z04C+wUPSVOWKVd5PuGH2bUagSGrKCsk3NI\n",
       "kJ8EXnx9vvMks9v2lI56xhRL9wdKDDmxBIeZc7yzg1Dmtjr+sA/8H/oXU5ZogICPmyKpgnp5GliK\n",
       "UBTZ4b58HpRB4IMNNCCOwGyi44B7KeyecUHdD2i4ZYEAAAHFQZpoSahBaJlMCCP//rUrwhDgAAtU\n",
       "fVm711LUsFysGu+GJzGVQ8hl8f2W9YhHTPbDazkABJxHPd2jjRM8Lf8bAMhuvruZDpAoj0MSjJfI\n",
       "AGgZnc29pPCiv6HG0AxoEHYJiqPaLzO0u9WOUsXUJkmekQ5mYtlBewjRsDWfvWpj+AkJ92NeUDAZ\n",
       "t48MHhuLb9W5rhKzmxZSkVmym0olRaGttJydq9HKv0VdVJM7xMOqDAVLFBq+O5CJ72BB7hQ5nrD9\n",
       "b9msWBURMBg/Dhe9qAo6WIxbFUn+JxLrbFioB9Nkc8Zp0f+HblQ3h4T8vt5RytMKr4fIsPGgVBMR\n",
       "eN8OZDUPuFhN88vN7E1HzKa+tL4M4qvvb0qEkA9jnUVG+t8m6KSBFm32Up9bLDsPvoOsEUTzHgND\n",
       "7Mey/eJIhZWoL3CBmXmWMlsQKS+fE2DoTfEPcr00etssoqV9BlINwy8c8nBPrh9yPlbGH95o2nAg\n",
       "LVt8ffc5GUEy+J7tb1xyb1hVvwHkfhz4BG+UmBoaF0LI0dL01xNV80WU0RZIyoG7uv83NHlC19CJ\n",
       "sQtopn4smuxw3o39MEAFQPpu6CEiz0GwqSrlFyVGvc4NAAABPEGehkURLDv/BZaHywAhTzxjFVcV\n",
       "uOEUYjAprq4ZX3rQUuZuPGL7gNd7813Cd+lVvkmZ7NeV5Jfx8/RlYFM02C/7dkV3LVt++PcUNCoT\n",
       "B0gp/s8xGh5o6VfsObhI58MI30KRVS32Qy8PyB84CcUBLbVsFE1ugLLsOkXXEyYqF76PuVIneKUs\n",
       "uhz4YEF66K+mA3H+h9H1gSSjNm5N4SwVEHb0wlz7J5fv0B+QRBfJ+kdNrHsiF/NpD+lYldg1uf4t\n",
       "EP8GZJcLGESHP6l4GmZo3RS9X+uwTJtN7nuiqTep/OSS4yZy/MWy7HuQWo2osbj2NCNiKra7XSg5\n",
       "hBM2M/8W7MRRuvaPkyJIDLIKo1AVzyzviX8Haui/QE0dQaRw102f+vnlNnQoEBlO08rfX+435uo2\n",
       "JE57s9GdZRa3HSEAAAEgAZ6ldEN/B0S3ZABtVZ+qxBnCV0mvrQHwY6yyetLbmvwjsHckeR+/Whlq\n",
       "VGLpneW8sAZlMurNRQd/nmxWgtSBNjvUiMCZNTYc7hVzpnM17/m2neWhW+wOCIKofw87aHv4KDPt\n",
       "fchy6sCBoBLPko8bMpIGQIjNQxAcIrEmhdBu8MuFxBQ1lG09BHtldSy7MEtRzIJepQv68NviXwIu\n",
       "vu6q/iFIGCQ7e+tdk3avBzGrhIq07azN5nuNyBDKjRZxgd+jBlgxrQpfoYDHe21lr6wJNpY9FV6w\n",
       "nIWvzLU0lM6X0ksmKB+1OSH/nhKmwhvAKYGXwff7zNhTp3MscarzzHVBwvNnop9BECOyoVOYQCrP\n",
       "L2MGeHF1OHUI+L3PdkRTagn9AAABbwGep2pDfwcaqABaIAQp7zOV8zCreoZBxxdMwVIsHiPdLFpW\n",
       "k6+CDN05AxJLElmDc9i3X/P/nvEqDlao/rZ5rg7ZnsZxYbMQL1XjWv5GQ8FcHC9b3w2Cyr1wzkYO\n",
       "PBXaK4SDTMAPka9MNkwUemOVdgJHfsN2O0CKeW2DsEKz5rkOlD6FtNkwXj257cSaXM0ROi/CHw/n\n",
       "LaNfHl/mcrMBPQudpJC5YL8jyvZa56nFr9ziXEysKZQ8A1pIR5NBFuy4PnpobdDYAvlLGA0BBIBT\n",
       "bIkhExVCyCb5Cghd/sNXuJrnZvD9pGiyd+6JJF+U+fKOyAyVIaD0LIwC8aSTci7WS5GSxVgP/UVT\n",
       "Y3m1xlz/i1T6Z3TXWTD8OK+0KaaDDGDqsOIKVRwxIBZigZFIlMf+DAPrICAeepV95/kRRyMXNqfm\n",
       "cNxOadSI5IFTdpiy3+xYqH+ilrRbzZT4v/98dCLwf37OTZcKPJHXz0Ok+ot4R44AAAHVQZqsSahB\n",
       "bJlMCCP//rUrv3vYxlrW98JdOpHACdbgpiQ2Harp9ChCXVquNLZG5oh+rriNEdAjOoe3aSf+recr\n",
       "mJEVdGuXOkoMtArm2P21HGNJfU1JQDQfBOaPXTLX3/FWN6Q6/j86IsaBuAM6j1NftRa/6lFpsjmX\n",
       "NLNSYgrZLmeznISYpVoX//Z2gCnVPyqSFJQvUTOcMh8+WEbM+Oj2IkwUa/+f6AmmXY44PP97u0Ro\n",
       "Cuv04LUzp7JWIyMAPpF+PRxgrnx3pyds/dkChTeRDmIR3oI5JR6BdxB/eJAo827fHXRJR+d8uD3B\n",
       "GsTwkx4leeQ7LLxmV8MsWpgi5iCuzK6R+t4RsC+XmHdWw4q5ItKVXJoB2iD/BndPNRcm0nvAAsJU\n",
       "0w24FErKIdSHBL6QBM3yyFqUMciWP4HjsdxWzbgUQWLXGG2yZfeTm1L/lysSD1JjObVSgSJgIICv\n",
       "qB0A+D65/0zA9Fo2zX1VPPmF1QdmYA1wG5QxN1ma4MYe8JUmphvMuW1v5UEjN6O8fkLNSE/tWudo\n",
       "n8Ktcd+JliU+jbHFO4AthA4H7z8vi2wnwvEQov67bY02ZhxmEWKVM1crwubIKJvAxzqp2R38+86g\n",
       "+BZxKAPTgAAAATpBnspFFSw7/wUaZ06IzM48B8ADaYKZoNsZuCq1WGrGhrGQxWN7cEeLXjMXh2AI\n",
       "G9w/7Oo+omegBcFLghyuDO+Aa9r5sfM6pJRL/iEF3dpHBure74XmOfu/ou+ztLQO+uoNgrOibMus\n",
       "TDbph4pe9SPWCfhdRJNNBh0tn82AWBmV5wg9ZVfi19e3G7M9bhx1Z9XRdcJyUGbqySp/mXkYbC86\n",
       "d7+aXprHUXjBT9akGcxd044Gv+f4O5Di4Es9EU1h9UgOjSzTQrYhyvv2CvqrZo8qynggxsYzEH4B\n",
       "Vj/LABhKSUpvV7N3fEZkQrUVLxh0MalWpj+6t8e94qj1Wyc+BmWLhVCOPUa/IzmoqjlIgmq7K1gZ\n",
       "Ga5MN6UuQaZwRHRul+2ys54yBEgEYXWbpqX8DOHpdIKugVVwZoOahwAAAQEBnul0Q38HSmtEAB1x\n",
       "gKXTGie7dxcJf7UlxD94eUyEFdtDM1oLs0wHDOFUzWSZRWohFInBPjwUp8Q6GIzRprN+hEpjrSo6\n",
       "s4Rf6e3j2ezdPVUSLnYtq7iVEWtoKj2sZ2RNv6HGWRh6n+/2D+MTI31+u+2JsdVgiZr6BSTjFD2m\n",
       "llI8u4TXhJ5P1gGyJT2cL0ME6LknVv44ly4yLjmE8tX3kNr7RGfCeVcR/CKG8klNrPmj5yheWhbZ\n",
       "EoZ3pSIJ44kFfrmK4NNr4P0TvLab7Qh5zZEs08KGoq/e7SMm3IpKLyYnQ+4wEnoaDfFT9s96wj/q\n",
       "Qmd5wIiSYHwrlqIMG7ijrwAAAS4BnutqQ38GrOG2axvocQAa2cfa6PmgfsrfaO5F4MtPO0YEyzRk\n",
       "4k8seNr1osUqMwIWjfWZPi22GTxU/RSsOiVOHHsRQwBZ2lg+qnRXocYB1gUtwnnsd8JV8iW1aEFS\n",
       "G6IYhNCFVz2NWC5VuuG9USuYsBw8avu1Tn4YKgkxzr7OrReeAoZQDW/e3n3MzP53wJXdBJ1LZYjg\n",
       "uh/+5bXIHIgucn2UoOdT3z9kPxkfRzFM9bR88QRE1vuJuQMNNuzD0kqcInoKIfkw8idb3Qrw6ICK\n",
       "rBEOeMWE6a0cN38kFfsyefqGvzJYRNbKYTvOCxZsg0CWSZ3gpH/4oB0J0NaWWMD0XFD2jwahWJ3C\n",
       "xouFp7ImuPgFVN4FzM1rVZPYUZQmbPDTieeg/FS2pgcHD7bHTAAAAVFBmvBJqEFsmUwII//+tSvX\n",
       "51jIAHaUi6jNN9KYWxhTdnk4uHHvQ8ZKTl7r6XmtZIOZEQIHCpdJrOacxShlzUrl+75vmLdQEPWP\n",
       "R7d+IkvYDmfcnpHfbsxzIq3oHBGgNuh/jH4YnVxsRp/f0N9rtLV8ylzxSiK8DsnteIVKcF23Cg78\n",
       "1F3WWxWIpAIKhlWYesF63N+XjdLYdSkR5CsuwpEf19HWGVsu1x9Fud/i1OgHbqncDT5PVXYzCnME\n",
       "QJhXr2rjW7PyYBjJHinYcI9JFc2UCoTs61dTlfwK+Wyk/L1xTylTixOix2QBJBauzkyOcEBox4x0\n",
       "2MnU2L8ZreLFkA03Mt7f66XNVwvRfdafr6OtDzkUTJ3UhADO8MV5x65BAb7V4W1spbo6eRmp+BAL\n",
       "6WYSXkKUFWyrx5w7biK2LUdTWy6zkWjgpjb9f+PgToLmz/mhAAABKEGfDkUVLDv/BYUvg7A0AF5J\n",
       "5H2zbMFiQP/YHU5QeXqmgO+BzEj+G9lq6tW2by4VL9tFDpe6mkQT6Fjx1FTP1VbDLYux5uLGJkPr\n",
       "gygkcZw8Dkkg52KLpAhGUQqOS/DOTbd2OlHKtf7JXoHg9TcGMkLB3CDnuLx5SuTBUJl1hAueOM04\n",
       "SHNh3lPTt6nopW23+4EIC413uY5lxluT/6ES/0sdLAGvAtW40LSVf+zGiouGFWtpZliP298f2wag\n",
       "mYxbTGU0qcfQZnUB2A5AYG7hTVZhRL0M69bIQH1eKNaQ2kwN8QVe9bQQCcCVPoUMAneO7TYsxPL/\n",
       "/3cocuxOeVo3jplkEI7KROF5gcHhYllenucbSQXWvCfrgSGAuAsjrd5Wq32djsYlGOKBAAABHAGf\n",
       "LXRDfwc34cpAAIgbDg/1rTsveWxka3eD0wKfXRvNdo43oyw8qT6ptdbv2J5ifrz6/jUq7I7VUbVB\n",
       "1HGHzZ5qjBuGmG6sJYTOv64UI/PVVt4YR19SI7c6crBVnCz1lJd9vKP9U2TQREBPss9n6fTmFrwx\n",
       "GoYhq/cok5ZfYlozrelwy4RmYy1jLN/nkETA6Fm8/BmozacAzJSxkm+3zUpeCPoMT3NV4NOqEzj/\n",
       "SiKaqppOKujVTRp1IWDscFTlBoDFgoNQMVb3aA6FGeNfzNqxKrfCw/cBfpxTl9GDOj62McXfXFlR\n",
       "7JIUazMmCyPMEmR9iuREAq9sAYUCmk6U3OuYW+/9yb6S24xyTIlO2EmPVG4SWE5J3BveAyA5AAAB\n",
       "EQGfL2pDfwcm5htkAG1XFTSBQTLZfVXYyMJ/D+jvvV37fc7FJf1sBbNHqyDeS5zj4TMfWS/IdvPt\n",
       "N2S2bR6dhPHDUftuLnEw7DbDNLh+Zo9h9FMqwNSMEwZeDUPrp8Y12bZCzJqWItypQiHk+mMdoCPG\n",
       "J8rGaB866QTjaSimNBuX5IfEC9GiBB2tDHeOVsOsC6x846Z80vVw8SXwyDoYYhYMjZlvbK40ug75\n",
       "I6bhpoLgfW7Mwiz4PK6UXA/ywiLfemAzPJB3foOyLcHM3F4X/raM0PY81xWv16CWabpZWQqiUCWj\n",
       "oTL+ZCsQdbxZ72inMXzHHNsrnawxT+zAbfrgobcpDtwuxrn19eAjH7rW/mWG6AAAAZxBmzRJqEFs\n",
       "mUwII//+tSvrv3IWMl4mByADazoxcDOylZgWubHChmcRkZ1jMjpkrGKz7rUDE96KsOywX1FqsIbk\n",
       "IMRDTtazb1lR8bqWPIQ2b0BfnbzVfywZ/bmEbSvGYbD1RMe5FF1MlNbURB6q/4d2ji6CzLLA+jtP\n",
       "Eo8ecN2VejGpQHn4LRKs8pKJyD5BQdjDAJyYabCE/b9XoKYOIVKAup5GlYeCYUP1RWtD6oBG/2T8\n",
       "bmZwAw9V4EzS/RGehq9qnH9jvddCwWKMRsrnd2aD8tLZp/RI5zqG5hF2U1YRGawz4nF5lUMmvys/\n",
       "7t0Ac6mcC9ihwWxWiR4IFQNHN2ufpZ0lD5UfPWfC/uzrXqDJLJmOgLt3mLwYtaoYR7frdVtbsfPl\n",
       "88wwalZTDpfE5rSwiNXV8DK5nupzlcw3Pgb5p58hC0EuUHPm5uLmttNtRtZnCipYeTDj0o0siFtw\n",
       "bekO0WGdSu7TXf0cnQXSlOKGR+2AzfDv8bFzKpP2fgvR6/veIGKvSfe8YltQVeM1e1cl/mY6AYFG\n",
       "mkw1jJqAAAAA+UGfUkUVLDv/BRklXEaSODQAfcCzUgvHzVeluvfOXIPdTaqGsrKuCJeBUBhqtyj4\n",
       "iBuz7kX0t12eruK3ryeOrxk5uz8ArHQiKgFS1QX1fJyjg8pXhOYTFm4quWA9YA9QabYx6i5oTVy/\n",
       "YsrHCwNF9IRHJrVnIwsnHaLdc4k8fUzlp4ltm/DDLfQSeWE+i/zFocJy4ZEXYWGhvi7XLyew6SGQ\n",
       "lfs8teqr8OmZnabSqGat1ljx9a5YhlKYS1ifkWSoLG6wvxAWVmCCi4vmwejL/SvCcbRIbJQLrinp\n",
       "hcm6hQJG1GKPGt4qA0aJ3hesCJRLZpMJ2VCw2dVEEQAAAQcBn3F0Q38GrOGml1Qq0AAiB4O5okpM\n",
       "0VliimwwH+F/vuXB89i3X/P/FUL43zGxs86f+7RNi7+mZGzBzMvBoNXOl7vbTBGDD46QLJHfy25Y\n",
       "SrAjHnFszcKBVDG/eANbSQjE1zY5RoymDuCuC2HSjYugGbgb19AojUNLnUJ5g12XJcez/fgihXEr\n",
       "aD4c8Pj72xq6Xc9qjNxwhxwlePSYwB5gSsT/cWBVslrkRLChAkieLoSXW6a9ETdj6uqMjbAjdJ8T\n",
       "Ek3I91wRyCJcw79BhlySrTae5w5AM04u5lVgCaaiJMlUnBglo8tQv1L0Ohw6N1cVtXviPhsbqCpc\n",
       "mB4ccbVILlNtQt+KcAAAARgBn3NqQ38Gpdyl9KQSiABDF3/fdm3hefSlnf3+Y10UF+wW+GvxLwIM\n",
       "acK85zeDtauNTisObS+mBICqRCIQl9aO10yiVll6gpIs5+7AtjlbAMui3kQv4o0AEWg7dbFLbBC5\n",
       "GtM2sDWKt+8QCdFrAglWid6HwrOdG76lD0jamY8L1kNqWDA8UsP2Yh8f7Hqh2TWWiU/ips1tR6Rz\n",
       "KKV3ByATcd3bPbIJ0foQRhFfYdViLy4ofwPVkN5N+/rabSGgG9l/gAMvQEV8GY3i+hFniAnFY7Nw\n",
       "cdWy+m/gyBmwrdfbRNJMgsi6AGwIzo6NcPcbSUumPzEXs/KsOu78OmzaYz5k/Hd3HZpTEIeAhrjP\n",
       "m2FaiR+o8/3M4lQHAAABPEGbeEmoQWyZTAgh//6qV4UMwAD+kWvZePe3ILMPV9d6gC0xCc6IjcZh\n",
       "OpmoeUyntKGmfWSz90WgS3Vt3gprjczrsdDphxFU28eVOJSRzRPUV0cbSONNUnZ5ZokJdqu36qR1\n",
       "Z7JZTAoCT2e9wUsyc1+4l127m7WAVnaap49wZCa25vajp8ci4a7SnmXF3b1INbY18UKmEai9gBJt\n",
       "oGrP/bP9a2bAH+Dqb+ym1ETxi7Qb2ABdW0sPWBGLOUc1raH+mfbjDMVM37gF5xMUfAWaNx3ORAu4\n",
       "feFxrC6HvKA0/MnriFk7UrTKuahvbp9x0kQ5pJo/dEGSLWRFTjP5NzaQ73MsVYo76VeEKIhgyW1c\n",
       "sMEMm1pmXy8dBDsO+jTIw6LVXF3w1BOLFkRTAQ4BeViB/cYDn/Y2JqudWh0ZrmUAAAEnQZ+WRRUs\n",
       "O/8FkkZYATV8NKxv6TFlOqUqsDMQYgbVKPuPUF89LQ4wcuXjsaQKZ2GX1Sbyk21mmlbvlOztemWF\n",
       "K3st4ElDiCE+vkKeVF36qcwManOjUjfSsLQyCpLX1UvNdsIdC519wvpIrF5px10K8v2yJvux9NAN\n",
       "XgF46IAqwCH3fkoIbIPjxP9YtdWc2CpycmYJiDiHk1UMDeGs6B+lO/g0njov1UBPr98XbSfpAjk4\n",
       "08B2JZ/uFPxT3PRPqwdAOwAMIHwMevrtw+KAJHhMWia2znQmuG0OUANLLFjMVt+/+tvQbHtLw/g1\n",
       "ONPVp4C7gTpbaM/49INv2+REESAz/+Z8RU2EStH3R3ynJ7UQSDeOiOQJQCR55kVzkgBgJxafYnQ+\n",
       "6S4voAAAAPkBn7V0Q38GrOCOOHEAGygswDjGoqfvLLLq06Gpf/LMOxQIGrNKQ8LcIf7rS3PELak0\n",
       "9MD7B54IIdHwV2pbuBYqF6X7TCVqr1poMFFAqCVXFKXVBxSjo2rdZDsV6e67+iRJOFTsRkVOM1iR\n",
       "pKhFwaVBBI2wfi/zWpmeYl82cX3FM61ztDVIiRw2GdrwXR+dEdmU8olN33jvnOzWLMDCa8fMRp2z\n",
       "lJPLruQmfvqqDY+ugWeBHMIDifbmiJDjQJZxKy2ZgI19XUu1yGXYzXmVhXbQvSk/M92qCsa7UXaf\n",
       "yGUdp1Al0N3FPnffDNkH+S8vhaOip0gyovuoKO8AAAFEAZ+3akN/B0fIQACIGcwM5+IGAfOtt0h0\n",
       "vP2bOT7jhbkfCEHlSsBaQTrYtcxygKf9aF0vFbgrGhtVWOb1WIyabxByeVgmyJD0ypXXgPug8fSH\n",
       "9DEL5epC0h5WaLHKBuie9NWgRDZnM1qTdAYTmbVo98qSRPzkqZSQM2eArunv5GXqAO9qRY7RFwGO\n",
       "4oOOyXIH5WlCuQF0JPxK5wHQrBu8LBQVh2WWQl3HfL07e04f8OiPhxXuXXTOEkc4b7faikPRrWFK\n",
       "v0z/VVQ7+Y3C6QMgXTfrZdshrKk62K4Ju1b99gY6McmBOnttDwbQXtgOrDpnStj6yhpW+DjNRXE0\n",
       "tYhlb3rUliqq02K6tXRGNnUUF20JaJQ7aXnLvmKDHrlmIOhLM5nki9lUDo9vcC7GOXMvtgE2Pc0O\n",
       "CQn+R8v5yY7D2lareKexAAABAEGbvEmoQWyZTAh3//6pnGRjLACFMr8YhF7tp7INMCdviPP5QRdz\n",
       "dXRNkWYMsnQ+49P6r4Akt84BmiMBkaBh6cokYedSQ8de62J8/0Mbm45S/lYKJ5lhCtF7sApMtHB4\n",
       "9nWnqntiUwZi8W8glP7YR1HPu4Dm73bU9B35ZkoyZr8cNLEwAsiXjUVm5wE2007xm4eGe7MQf9Gb\n",
       "LP94ATjxkNQ35DctaumJRmbKaEZeDE8XpIr9s5DFithCh6uhDFaRZEiHMLnyc6R+/VpN/zAcnqJd\n",
       "8W6E9Pxja7FMRh5LxQcWBWg9F83ho6zmF/pOay3yGs3J3P+h1F/qnKb/edbkMF4AAADtQZ/aRRUs\n",
       "O/8FlofLACFPQ5WKq2+Bpp47VzDajIQyjyPL/iqbePkygLy2ITiJb/3RnmGE899oBn2kiqZyTAte\n",
       "QcMU0uhDggShfOd3giXr5O6D4efIOh8QnhFWB+CdlOGw4q7GNu/w2zZA1UPsaZJfjTckCD7C5MCs\n",
       "bfr85WFolafVcX+Wnc94ciaxkicNstQEm1CG1lUhr4+VJ0sUteIUdCal5XqcmwgBNm7a811ijL5r\n",
       "b31hNjepwqkghMHtcjGNiSLpQanZLAmnV/Z5ZadxpNsxkiskPbPThTmSRGzM/t5GsXa5sTszxscN\n",
       "+N5hAAABLwGf+XRDfwdEt2QAbVWfqtA+gGvGHGOUuDHVtJP0Ovv4R2Jw2CFp08aFEid5ZneW6zJW\n",
       "t8CjP2+d2Etzt3LF6v/y5RX/FqKa8vmtnrQcXOb+/dc+lN9iBl4Hp2jgK5m4RAg2jbXcKscKK0W6\n",
       "eoEmktQqvbQb6IkXnkr4Db8mHcgloR7fNci102M2Lpg6wgTW/ACc21K+38LNpOVnm7UOE06vCEJN\n",
       "iyC9NuuRqRUIoIHmq5gJoUTvX7QUWTERwxUadT/0be/CDLWqPN/II0tcxuqRawBaL/IXsv03jW5Z\n",
       "i+dRYLxkwrDl3APfDotf/1dA8YLpiqb2Cqt9A3Se9zSAYi+PhGmyt/HC2bW5H35UIafthaY2uAa4\n",
       "s8gwnFnAI6PWxpQAXIjAkxFVyP4n/zHoIAAAATMBn/tqQ38HGmXWQAbVcl3V5uhb/C3zybmQ0pFg\n",
       "5qGkOOCWQRbYt8kRj7e8hr0kOOr/z/719EYRlcq9eQARaj7NRaHgbYhIMvawAmIFuauH2wEGmD2Q\n",
       "UauCjqD2qm0mSufCzGrm6dKHiEcL1d5d3vDVXPGPvonm8RRxd/7V4ZTV3uc5Z1ilD043XsJQGM3y\n",
       "CTcBb0ZpuYNb8zT80sq6/LeNNd/TsiPQzJngNtfbfT2hPWAb9OWbMevR7V6EHreGbgyHyZGU4uEs\n",
       "cwDmIiVvZaJcFmc0F5syedAtLPrVil4GJfMmIzAFws3wagApaynLQwRHheRDf/4ZS/iboPxDvCSI\n",
       "TU4oPCNDVSAVdptArYEi2wvbsf7o5mV5Sk1ig4xJqblQfRdRcQqGmE/oi9GVjHZ3PK3BAAABC0Gb\n",
       "/UmoQWyZTAhv//6njkzkOAEH11/n9d7ieGjjjayjkaNGhtLO2QauFo40jie9o5/1xvunb6bP0lJD\n",
       "qf8xNFglEoIFNRONOV7nSdHdApz5JSHdpLXoeLwxf1MG3cRyKfXDXbNHaeskxMUpLb/n6z+PibXt\n",
       "iUOpNAxiAiEL7fJzCHfj9bFfdq4UqFHiqq5GrgoQHDA1iUJqyeSpMxfvg6c+PTYltCyLdwqRTSW0\n",
       "3a3dJ8TZRHSc0xtQRQd8zJ+PbUY5oq3Y1mfhvCsq0np9JVymEIklWzTClFpiveSTjSLFL247GEai\n",
       "enXSZpYeUyNLBGj3wGaXzY8MaMR1yCUPsnZ6nbAl1CVQ44HqMwAABJZtb292AAAAbG12aGQAAAAA\n",
       "AAAAAAAAAAAAAAPoAAALuAABAAABAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAA\n",
       "AAAAAAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAADwHRyYWsAAABcdGtoZAAA\n",
       "AAMAAAAAAAAAAAAAAAEAAAAAAAALuAAAAAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEA\n",
       "AAAAAAAAAAAAAAAAAEAAAAABsAAAASAAAAAAACRlZHRzAAAAHGVsc3QAAAAAAAAAAQAAC7gAAAgA\n",
       "AAEAAAAAAzhtZGlhAAAAIG1kaGQAAAAAAAAAAAAAAAAAACgAAAB4AFXEAAAAAAAtaGRscgAAAAAA\n",
       "AAAAdmlkZQAAAAAAAAAAAAAAAFZpZGVvSGFuZGxlcgAAAALjbWluZgAAABR2bWhkAAAAAQAAAAAA\n",
       "AAAAAAAAJGRpbmYAAAAcZHJlZgAAAAAAAAABAAAADHVybCAAAAABAAACo3N0YmwAAACzc3RzZAAA\n",
       "AAAAAAABAAAAo2F2YzEAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAABsAEgAEgAAABIAAAAAAAAAAEA\n",
       "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAY//8AAAAxYXZjQwFkABX/4QAYZ2QAFazZ\n",
       "QbCWhAAAAwAEAAADAFA8WLZYAQAGaOvjyyLAAAAAHHV1aWRraEDyXyRPxbo5pRvPAyPzAAAAAAAA\n",
       "ABhzdHRzAAAAAAAAAAEAAAAeAAAEAAAAABRzdHNzAAAAAAAAAAEAAAABAAABAGN0dHMAAAAAAAAA\n",
       "HgAAAAEAAAgAAAAAAQAAFAAAAAABAAAIAAAAAAEAAAAAAAAAAQAABAAAAAABAAAUAAAAAAEAAAgA\n",
       "AAAAAQAAAAAAAAABAAAEAAAAAAEAABQAAAAAAQAACAAAAAABAAAAAAAAAAEAAAQAAAAAAQAAFAAA\n",
       "AAABAAAIAAAAAAEAAAAAAAAAAQAABAAAAAABAAAUAAAAAAEAAAgAAAAAAQAAAAAAAAABAAAEAAAA\n",
       "AAEAABQAAAAAAQAACAAAAAABAAAAAAAAAAEAAAQAAAAAAQAAFAAAAAABAAAIAAAAAAEAAAAAAAAA\n",
       "AQAABAAAAAABAAAIAAAAABxzdHNjAAAAAAAAAAEAAAABAAAAHgAAAAEAAACMc3RzegAAAAAAAAAA\n",
       "AAAAHgAAD/AAAAKwAAABfAAAARwAAAFyAAAByQAAAUAAAAEkAAABcwAAAdkAAAE+AAABBQAAATIA\n",
       "AAFVAAABLAAAASAAAAEVAAABoAAAAP0AAAELAAABHAAAAUAAAAErAAAA/QAAAUgAAAEEAAAA8QAA\n",
       "ATMAAAE3AAABDwAAABRzdGNvAAAAAAAAAAEAAAAsAAAAYnVkdGEAAABabWV0YQAAAAAAAAAhaGRs\n",
       "cgAAAAAAAAAAbWRpcmFwcGwAAAAAAAAAAAAAAAAtaWxzdAAAACWpdG9vAAAAHWRhdGEAAAABAAAA\n",
       "AExhdmY1Ni40MC4xMDE=\n",
       "\">\n",
       "  Your browser does not support the video tag.\n",
       "</video>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create an animation of the traffic density.\n",
    "anim = animation.FuncAnimation(fig, update_plot,\n",
    "                               frames=nt, fargs=(rho_hist,),\n",
    "                               interval=100)\n",
    "# Display the video.\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll see that the result is very similar to the original Lax-Friedrichs method, and with good reason: they're essentially the same! But this is only because we are using a uniform grid. In the finite-volume approach, using the integral form of the equations, we were free to use a spatially varying grid spacing, if we wanted to. \n",
    "\n",
    "The original Godunov method is first-order accurate, due to representing the conserved quantity by a piecewise-constant approximation. That is why you see considerable numerical diffusion in the solution. But Godunov's method laid the foundation for all finite-volume methods to follow and it was a milestone in numerical solutions of hyperbolic conservation laws. A whole industry developed inventing \"high-resolution\" methods that offer second-order accuracy and higher."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Dig deeper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Godunov's method works in problems having waves moving with positive or negative wave speeds. Try it on the green-light problem introduced in [lesson 1](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_01_conservationLaw.ipynb) using the initial condition containing waves traveling in both directions.\n",
    "\n",
    "* Investigate two or three different numerical flux schemes (you can start with van Leer et al., 1987, or Google for other references. Implement the different flux schemes and compare!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MUSCL schemes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Godunov's method is first-order accurate, which we already know is not appropriate for hyperbolic conservation laws, due to the high numerical diffusion. This poses particular difficulty near sharp gradients in the solution.\n",
    "\n",
    "To do better, we can replace the piecewise constant representation of the solution with a piecewise linear version (still discontinuous at the edges). This leads to the MUSCL scheme (for Monotonic Upstream-Centered Scheme for Conservation Laws), invented by van Leer (1979)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruction in space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The piecewise linear reconstruction consists of representing the solution inside each cell with a *straight line* (see Figure 5). Define the cell representation as follows:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "e(x) = e_i + \\sigma_i (x - x_i).\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "where $\\sigma_i$ is the *slope* of the approximation within the cell (to be defined), and $e_i$ is the Godunov cell average. The choice $\\sigma_i=0$ gives Godunov's method.\n",
    "\n",
    "Standard central differencing would give\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\sigma_i = \\frac{e_{i+1} - e_{i-1}}{2 \\Delta x}.\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/cell_boundaries.svg\">\n",
    "\n",
    "#### Figure 5. Piecewise linear approximation of the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we saw with the results [in the second lesson](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/03_wave/03_02_convectionSchemes.ipynb) that this can lead to oscillations near shocks. These [Gibbs oscillations](http://en.wikipedia.org/wiki/Gibbs_phenomenon) will always appear (according to [Godunov's theorem](http://en.wikipedia.org/wiki/Godunov's_theorem)) unless we use constant reconstruction. So we have to modify, or *limit* the slope, near shocks.\n",
    "\n",
    "The easiest way to limit is to compute one-sided slopes\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\Delta e^- = \\frac{e_i - e_{i-1}}{\\Delta x}, \\quad \\Delta e^+ = \\frac{e_{i+1} - e_{i}}{\\Delta x}, \n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/calc_sigma.svg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 6. One-sided slopes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now build the *minmod* slope\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "  \\sigma_i & = \\text{minmod}(\\Delta e^-, \\Delta e^+) \\\\\n",
    "  & = \\begin{cases} \\min(\\Delta e^-, \\Delta e^+) & \\text{ if } \\Delta e^-, \\Delta e^+ > 0 \\\\\n",
    "  \\max(\\Delta e^-, \\Delta e^+) & \\text{ if } \\Delta e^-, \\Delta e^+ < 0 \\\\\n",
    "0 & \\text{ if } \\Delta e^- \\cdot \\Delta e^+ \\leq 0\n",
    "  \\end{cases}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "That is, use the *smallest* one-sided slope in magnitude, unless the slopes have different sign, in which cases it uses the constant reconstruction (i.e., Godunov's method).\n",
    "\n",
    "Once the *minmod* slope is calculated, we can use it to obtain the values at the interfaces between cells.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "e^{R}_{i-1/2} &= e_i - \\sigma_i \\frac{\\Delta x}{2}\\\\\n",
    "e^{L}_{i+1/2} &= e_i + \\sigma_i \\frac{\\Delta x}{2}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $e^R$ and $e^L$ are the local interpolated values of the conserved quantity immediately to the right and left of the cell boundary, respectively.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Index headache"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that for the cell with index $i$, we calculate $e^R_{i-1/2}$ and $e^L_{i+1/2}$.  Look at Figure 5: those are the two local values of the solution are at opposite cell boundaries.\n",
    "\n",
    "However, when we calculate the local flux at the cell boundaries, we use the local solution values on either side of that cell boundary. That is:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F_{i+1/2} = f(e^L_{i+1/2}, e^R_{i+1/2})\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "You can calculate two flux vectors; one for the right-boundary values and one for the left-boundary values. Be careful that you know which boundary value a given index in these two vectors might refer to!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a Python function implementing minmod."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minmod(e, dx):\n",
    "    \"\"\"\n",
    "    Computes the minmod approximation of the slope.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    e : list or numpy.ndarray\n",
    "        The input values as a 1D array of floats.\n",
    "    dx : float\n",
    "        The grid-cell width.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    sigma : numpy.ndarray\n",
    "        The minmod-approximated slope\n",
    "        as a 1D array of floats.\n",
    "    \"\"\"\n",
    "    sigma = numpy.zeros_like(e)\n",
    "    for i in range(1, len(e) - 1):\n",
    "        de_minus = (e[i] - e[i - 1]) / dx\n",
    "        de_plus = (e[i + 1] - e[i]) / dx\n",
    "        if de_minus > 0 and de_plus > 0:\n",
    "            sigma[i] = min(de_minus, de_plus)\n",
    "        elif de_minus < 0 and de_plus < 0:\n",
    "            sigma[i] = max(de_minus, de_plus)\n",
    "        else:\n",
    "            sigma[i] = 0.0\n",
    "    return sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evolution in time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are aiming for second-order accuracy in space, we might as well try for second-order in time, as well. We need a method to evolve the *ordinary* differential equation forwards in time:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "  \\frac{\\partial}{\\partial t} e_i + \\frac{1}{\\Delta x} \\left[ F \\left( x_{i+1/2}, t \\right) - F \\left( x_{i - 1 / 2}, t \\right) \\right] = 0\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "A second-order Runge-Kutta method with special characteristics (due to Shu & Osher, 1988) gives the following scheme:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "e^*_i & = e^n_i + \\frac{\\Delta t}{\\Delta x}\\left( F^n_{i-1/2} - F^n_{i+1/2} \\right) \\\\\n",
    "e^{n+1}_i & = \\frac{1}{2} e^n_i + \\frac{1}{2}\\left( e^*_i + \\frac{\\Delta t}{\\Delta x}\\left( F^*_{i-1/2} - F^*_{i+1/2} \\right) \\right)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Recall that the Rusanov flux is defined as\n",
    "                          \n",
    "$$\n",
    "F_{i+1/2}= \\frac{1}{2} \\left[ F \\left( e_L \\right) + F \\left( e_R \\right)  \\right] -  \\frac{1}{2}  \\max \\left|F'(e)\\right| \\left( e_R - e_L \\right)\n",
    "$$\n",
    "\n",
    "Armed with the interpolated values of $e$ at the cell boundaries we can generate a more accurate Rusanov flux. At cell boundary $i+1/2$, for example, this is:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "F_{i+1/2} = \\frac{1}{2} \\left( F \\left( e^L_{i+1/2} \\right) + F \\left( e^R_{i+1/2} \\right)   -  \\frac{\\Delta x}{\\Delta t} \\left( e^R_{i+1/2} - e^L_{i+1/2} \\right) \\right)\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Now we are ready to try some MUSCL!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def muscl(rho0, nt, dt, dx, bc_values, *args):\n",
    "    \"\"\"\n",
    "    Computes and returns the history of the traffic density\n",
    "    on the road using a MUSCL scheme with a minmod slope limiting\n",
    "    and a Lax-Friedrichs flux.\n",
    "    The function uses a second-order Runge-Kutta method\n",
    "    to integrate in time.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    rho0 : numpy.ndarray\n",
    "        The initial traffic density along the road\n",
    "        as a 1D array of floats.\n",
    "    nt : integer\n",
    "        The number of time steps to compute.\n",
    "    dt : float\n",
    "        The time-step size to integrate.\n",
    "    dx : float\n",
    "        The distance between two consecutive locations.\n",
    "    bc_values : tuple or list\n",
    "        The value of the density at the first and last locations\n",
    "        as a tuple or list of two floats.\n",
    "    args : list\n",
    "        Positional arguments to be passed to the flux function.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    rho_hist : list of numpy.ndarray objects\n",
    "        The history of the car density along the road\n",
    "        as a list of 1D array of floats.\n",
    "    \"\"\"\n",
    "    def compute_flux(rho):\n",
    "        # Compute the minmod slope.\n",
    "        sigma = minmod(rho, dx)\n",
    "        # Reconstruct values at cell boundaries.\n",
    "        rhoL = (rho + sigma * dx / 2.0)[:-1]\n",
    "        rhoR = (rho - sigma * dx / 2.0)[1:]\n",
    "        # Compute the flux.\n",
    "        F = 0.5 * (flux(rhoL, *args) + flux(rhoR, *args) -\n",
    "                   dx / dt * (rhoR - rhoL))\n",
    "        return F\n",
    "    rho_hist = [rho0.copy()]\n",
    "    rho = rho0.copy()\n",
    "    rho_star = rho.copy()\n",
    "    for n in range(nt):\n",
    "        # Compute the flux at cell boundaries.\n",
    "        F = compute_flux(rho)\n",
    "        # Perform 1st step of RK2.\n",
    "        rho_star[1:-1] = rho[1:-1] - dt / dx * (F[1:] - F[:-1])\n",
    "        # Apply boundary conditions.\n",
    "        rho_star[0], rho_star[-1] = bc_values\n",
    "        # Compute the flux at cell boundaries.\n",
    "        F = compute_flux(rho_star)\n",
    "        # Perform 2nd step of RK2.\n",
    "        rho[1:-1] = 0.5 * (rho[1:-1] + rho_star[1:-1] -\n",
    "                           dt / dx * (F[1:] - F[:-1]))\n",
    "        # Apply boundary conditions.\n",
    "        rho[0], rho[-1] = bc_values\n",
    "        # Record the time-step solution.\n",
    "        rho_hist.append(rho.copy())\n",
    "    return rho_hist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set time-step size based on CFL limit.\n",
    "sigma = 1.0\n",
    "dt = sigma * dx / u_max  # time-step size\n",
    "\n",
    "# Compute the traffic density at all time steps.\n",
    "rho_hist = muscl(rho0, nt, dt, dx, (rho0[0], rho0[-1]),\n",
    "                 u_max, rho_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video width=\"432\" height=\"288\" controls autoplay loop>\n",
       "  <source type=\"video/mp4\" src=\"data:video/mp4;base64,AAAAHGZ0eXBNNFYgAAACAGlzb21pc28yYXZjMQAAAAhmcmVlAAA2eW1kYXQAAAKuBgX//6rcRem9\n",
       "5tlIt5Ys2CDZI+7veDI2NCAtIGNvcmUgMTQ4IHIyNjQzIDVjNjU3MDQgLSBILjI2NC9NUEVHLTQg\n",
       "QVZDIGNvZGVjIC0gQ29weWxlZnQgMjAwMy0yMDE1IC0gaHR0cDovL3d3dy52aWRlb2xhbi5vcmcv\n",
       "eDI2NC5odG1sIC0gb3B0aW9uczogY2FiYWM9MSByZWY9MyBkZWJsb2NrPTE6MDowIGFuYWx5c2U9\n",
       "MHgzOjB4MTEzIG1lPWhleCBzdWJtZT03IHBzeT0xIHBzeV9yZD0xLjAwOjAuMDAgbWl4ZWRfcmVm\n",
       "PTEgbWVfcmFuZ2U9MTYgY2hyb21hX21lPTEgdHJlbGxpcz0xIDh4OGRjdD0xIGNxbT0wIGRlYWR6\n",
       "b25lPTIxLDExIGZhc3RfcHNraXA9MSBjaHJvbWFfcXBfb2Zmc2V0PS0yIHRocmVhZHM9OSBsb29r\n",
       "YWhlYWRfdGhyZWFkcz0xIHNsaWNlZF90aHJlYWRzPTAgbnI9MCBkZWNpbWF0ZT0xIGludGVybGFj\n",
       "ZWQ9MCBibHVyYXlfY29tcGF0PTAgY29uc3RyYWluZWRfaW50cmE9MCBiZnJhbWVzPTMgYl9weXJh\n",
       "bWlkPTIgYl9hZGFwdD0xIGJfYmlhcz0wIGRpcmVjdD0xIHdlaWdodGI9MSBvcGVuX2dvcD0wIHdl\n",
       "aWdodHA9MiBrZXlpbnQ9MjUwIGtleWludF9taW49MTAgc2NlbmVjdXQ9NDAgaW50cmFfcmVmcmVz\n",
       "aD0wIHJjX2xvb2thaGVhZD00MCByYz1jcmYgbWJ0cmVlPTEgY3JmPTIzLjAgcWNvbXA9MC42MCBx\n",
       "cG1pbj0wIHFwbWF4PTY5IHFwc3RlcD00IGlwX3JhdGlvPTEuNDAgYXE9MToxLjAwAIAAAA1UZYiE\n",
       "ABH//veIHzLLafk613IR560urR9Q7mffS/DfcIgDmIK8EUV8JDFX7eGduuOz73/AIdcyIgb3io39\n",
       "xyM7hq1HkbIOtWp3FhaPByIRoG4HjU4LmQqJ5Z21ubyVwpfHX252xjtdLlZUbcH8lD/+iMTkgRlf\n",
       "1PScYBJIBRXse15wHlBxdu/CfrdIXZ7F1p4oelRTh8xgEGNbLn9W2vdyQvyQTbdy6v/wMm5hqbcV\n",
       "txyHWni05C2jVkBE1AuVX/jntNdezff0w8CKGf6mHv0GSAtJlMVdV7CMd073XzAuq0qN841Ctrib\n",
       "QHfVATh9hj7qqiTVfR4Rubg+HBwxCWNLvwg+nFCWTn4TzNzUz9bjEdYHcHHTn8nQbj1nRib3AvoQ\n",
       "ND2wLAGW3zC80k9YkgqE3yo3ZMrnLwI5qFqNiGxv9nULUh88CUfn3zMDGPHgvsEFgt+LjdKmc+9X\n",
       "xOErWdjAN/hpvEm9073xzLwuvVDn69nswDTQ4nvIqxrLlmWFLg1fwPpC40bFKhI1N4Exv4hUsEgY\n",
       "wxnmZPJkGJkyQbxDOCZGcZIpmo67JjTi/xz10r/9HiqQVYjFJ1NRVNPwUHP6jtSYIrytx4nMcXig\n",
       "YkukyfWwo539/H22N87LJMru7iJ6LJwN8AXBp3rOPPtwxm5iTb/6FnXV/+u+7QUTvVzktLJQ6ETY\n",
       "AI4vIncE1U+iYY7f8LTelU19cI68Ktg6TaGv8Ah/YxRxSwLalIz1MLb2NNUe6CigS4GcP8Ti2oLI\n",
       "f3j0zvb+5mERb3ju/Hudnrfcf5tzbw/+hTgCLL8YjC9bFmEJdc8zqrrPny+EhZYgMoOHMCggpw2P\n",
       "PXiu1kzR64rgwnWDS+eAum6Wf//hhZmMnWtCqCk4NmTpP04GfNZguaKFkeFepKkRhr74j8APmb4Z\n",
       "2fbVNK7mPaZLvgAV4kvlilymX49aDZMmgsTND6NyGHVykE5nFGuUXyvEssQ3Iexrl4WcTvH8D5DV\n",
       "HRuGuhNOLRvRJr5gNUJvmZX/BTYBWKr7+ZfGmfxsjzunxQvRDzQGQ0NvJwsd9Igc8jklVsOFx+qf\n",
       "nwrbkb0UuQHD9H/PxJGujymFwcWIUQgwC3V9/nbfwLLZFa6KP7inD6s9ZoLtwuMhOaS5MqBGGO1F\n",
       "SIn/F/UYw7zqCTYF/C4KghE2/piBvvwXM+gvYTCSLtuDkiufwrcu44ckrSjjSIe9MVf1icksiSc0\n",
       "WPxbVdc+oFjvl42bPcozuq2Gls11ynQZyJBpGrT/PIuA7c1JdnBFpbfMPXX1ba41KyZQjY67Muvx\n",
       "iH3QGZzPKLGrT1SQOAHzrjT8/KN0LJulWZjUAUuMpPMXz4eNLeB0YY5VM/8Bv+F3SOOGJEgg9/Gr\n",
       "H4Gir2efrChioSI/D4PzHS5KoGScDFjEQgoXNHO4kjlhvlNvcHA90ngBnZukbs1JoK9h3I5ipmRZ\n",
       "IeFZLA3S0k9mVUEE5Ff+tEiR699iSIjTQ/MXrOFMRPdKkxQyF23Wxm1eSYxKnWWWbKk4Cyz6CpgY\n",
       "4ETIQqeoI3oOiWvJSqTvog4kWcdd0ZndifytlXHcO7UshrxlzkmGe7BsMpWTG0l7nLIDzvye/DZt\n",
       "7zqR9KgRueXzRF5f0ZXStHoVw0gyY2rjtxXEn9ra8KUa9ilaiHqEoF1SU3z0WXURhFVq3s14Jjq4\n",
       "jlUzEEe5pDBkH30pAZ49RgyDNMwY4H2DUE1UjHP7GCnb8cvtGLNKbDpiVi71iGgWf9PlPm/KWPOh\n",
       "jG+Z5t102i1F4LJXuYX0BEJXxYg6WPnjyyKjVjEe0Gh0wsUAGFDNsy3SfbewsdLsyOkWppyWyTaP\n",
       "iamN0HsQUgVpTqxJxnyN1ecN9/09dFI1DuQzP5GK2+uUMgLUryfj3pbTSn4k12buJ3Hz/5vErgyq\n",
       "OCBsOxK58k/vycUY19jGJIrm4Ak+N0R8HDvuWJKqdNjd3WPbpsUHDs3w1KVgZWC8r/MKcqVg0NPc\n",
       "ZTDL1x5+gP7+vo8ZNwf8GTw5P5GFljggpBndtvk07f7xagBpN3T3i7gYxG4mKFcwNixQAAADA8wc\n",
       "WsZ/0sn8gmQB0A5D7knCHiuEjsCNlo2DoKWqDQ1+/712XuEZd90/Ed631mB0VXJ/A7KpEBwKcTg6\n",
       "S3+MkDfonXhNvo/HXNnJGgo463xlcWnXpNbbWoOSkTK4apPJ03zBYwT/AuwS8N2nf/QvUqiWbSMN\n",
       "weJ67Wy6a0mB8sqtKgBs6ibzhR01bUNTJCXQ37L6l1MUUJHNtVDN9OpkPr9CLXaZnFlOF5xRPdsP\n",
       "vuMCf/HeJlNMsZsQaBPFM+ZEt8A1PBJMg1/10OJDZB6352f0dl8LdQLzAfNIKEISUSVIsOo4/YWZ\n",
       "hquehe68QQv4cI2wkWfIdfB1xY+UplbyN2mA2ujlabgHRsPQfMcAx5k6ib1o+S+bM4fRKfG3cj5J\n",
       "fWBjgiXZkSwNv/Re03glNwqxwxqxwkiOsb+PJLuvGAt0XEfA5p0bOuq2ZRKweTCQMIgQhaVTWKLA\n",
       "zNl9XKNqdb37PG0bsACL5BlVNiJC1jGfueaODAzea0dH7ULQar7cI8x9tPvJY4VrXOWfCB81SMEj\n",
       "P4WMuBO7GXeeabro1pVNl2pipOO0ZVN7TM0HB+5boUMor0mjVgKs5ugQz+spFnBwXtx3C32p/NAb\n",
       "ge3xpN8pKOmBVxeqo2M5aGNyCOAGwXSmkc5FQn7p7XAOIx3UrUUgcjDCHzMv02fxSXaySYFv6zMf\n",
       "dJLIOH46mZeCDmsNm75HeYEkHdVh3w18MCCmAqqowL/ZFo64LoQooLVH6gYQuan94hMogCof9uGz\n",
       "0uDQ35Ip+tzkZ6yRqhGhWka0nPkfvsx0jn2ZXk9ZXqXso8FId12J86F7ZLQxpfMuvgCoVtSj110s\n",
       "GE9YtAVF7PPUOisIAOFLoKFHjHYSnVHqC6zVVRRzQJZw+ViruAOP0si3TZs1HAOE5aO+BqARCB09\n",
       "Nv9aJ37L0zL5AfkBmtmw2879DqpC+51UKrUKKmfgzMMzn0X//Y5Wjhr58G+ZG//n8qoNwFNfjHfU\n",
       "zrIyvg/KmPXtEJ4cO3NImWP7+ZoNjCTMHTpH/aTbEj+9tuEnZrS5NQp5SPigOrR9Zb+KjZwXYFE/\n",
       "iLRmwPYsEcNEc2jDmn3XFpFpRa7aN75VGXvDa7IUVJfeg7qDtOztD7SZUMH/c5u9AZMFa9Uzlepa\n",
       "t20y6f6vfCmI42ZMfcSWro5k/cHFhBLPaV1+sdZ8hP+jGhdCbPNmgATB1qKdfqXUOJrJUXgp4skM\n",
       "DDaSHpIqn1f26Nd7p9iLF2zJSGZ/KkdpuO/6nwZ9BaHM9Zms4WR0cK+t0WQ3RVpkMw/eZsD74TSv\n",
       "WvHQ1eGXVZc0yXDBx0dMyoDPfoOMqGpmiRgFtC3RGXe092oEfOMSu8fvqp7JrOXUtSapptRJvU4N\n",
       "OCdx6DRoRMwDMtWrs9pD42hqnkEZRjFoPFBQsqebQ6kPTuQ/W69qpMsbhlkVPrHpg2WMIRllDCS2\n",
       "65+48Hw+D2kksWC4X1EECi9Avf0r3akxw4RWWlPSreBllBkA9HH5yXvPKdi1RI96Ctq+Cp60zV9o\n",
       "Hd7t+gtM+W4Jr3rzUx2WaFv7m7fjVZeKU2iR6oiU/v21qJBtnjwbPV27GqBFoZUo0CoXVZ0vHFWU\n",
       "vrJ0GI/5KV565TXxPu8SXZXmmMUMBK/dYuZgs9QCQm1XirXA+s/iSogrPlNuasTlPye+Qi5P3vtD\n",
       "v4aNJoTYSqahpt7jBTBHOaLYUr9lTFOMAiEvwwsPJNBWuZEpDmp9RZCp4NULLkeyRyCSyvOZDNKW\n",
       "J/6DUJHcAy920c2/8G1Z6PLz841smyT3tNK/Y6FFJuSwp3EleGDgr7H87qh8eA4hBXvibDuuyQZO\n",
       "D4phAMNXAOXwkg5TR3C/Bti+i43o3V50V+WIcVdQ/33ux1BfvL1VQev0ajvSvr6TeCOWePhugoyx\n",
       "vO4s/J0wZT8yAS2DSKKs+e1xZukZQ13YMGiLyxZd4Cr+Zw5hLKONCMraUAWGrwKIp1fdPiXC07W8\n",
       "57yAbMdODtf8K7j+5Tv3TFoL6p1U9yFkOI56Em9imsyTzEQgPy3AYq+R3hC2MT4HOL0wGFCrL4pT\n",
       "vPbPaF8A5K1QEa0PNCXBjw2DKbiqd5t8lWyKKkoqRyHcdoebvYDEuG3AnnDvgs6wsRRarTuqDRcc\n",
       "FMnso6gEs+C7I1/+XFa3zGoW5KRXyywjxct8a/17fa0QGGZN8LnDdn2JDmv/X1+wm5YrUJz2Vfe3\n",
       "ARHHhw28jQkaC3MuFbphdJomMF3ceSNVO39sef/WaUD3Ah1l5gTONsjCKdNh33l4ypjnPN6qD5VV\n",
       "7D6BOpBFqtJg1ENwTNxb9rtUyxvUet0Fdpf+7/uFfg6KZgIfeGmfNJqnBAj4cE1icWSrnBvL0rvA\n",
       "wPHBVJEeh98KP1b7A17Rkltq0Ow3cSVI6y43pKJBuSMB5balPKmXpYy+b13teedOQ+vB5fuHxLnT\n",
       "TR1SSjxHlS4AVuX/x4ublH366zsluer05//PYUHa0vl059qbWASddWXlwABBXwAAAtdBmiRsQR/+\n",
       "tSvFysgAdpefR32FldUWIgApB/Mt7rwj4EjuOchcWw33hecWpGq7uzaiWKmqJ9f6le8TRC+jTivr\n",
       "ny4yoVFkG97ABcNs2XzeBIihCRY/O1Ol9CEX3lFMCwaPcUeGgUWH/D3zBJ8ASpfc7siAWM53khw1\n",
       "OsdQEYPRYB9sNxmwmOqPZBotG5KlUzehaYnpXG0vbHk8IqNONKT3970QL8ITBS1E1Vr4By6mxBaT\n",
       "QojZeX0Zc5bwd00aQv+ahWygXIuiyiz3QRZzq5t8jAcanlEI5tqe05f5804Bxrp/k7HpgNZrIy8Q\n",
       "aO+s8DsqwixrNhi8JHHowePGhdT9dS5tThKdV/BJixuhrDfE8myGonm4GALFY5NzeT6ZI6j+UjTh\n",
       "IuWn60xgBL2ml2BugL1MQrtpVi91hrJD57DpBTXJ7pGY453Dk+8gVAr3rcxr52zaaxfe6UtbpzpG\n",
       "EsYNkjwHLIU5TyfN5weyi+/m57pqLuyUseGI/5n3CMQMLYaOpeBHbNUNBq913SkFMPhuq6LMMsyb\n",
       "En/VBbxBGilhZSWUwC9KoICdY3Lx5BKsQEoJaEV3Sb/GdnCMCoDY52UOOkbOYzdDW2XQE7mvV7VM\n",
       "ZgkcsXNTxW+dwY1f+Z9w3WbKRHXzueUlo7r5mTqH5GwZX08D+JUlEjuxSaCUhB2GhacMuQP8bO1a\n",
       "qMATnOjlaz8CBbBIJB1AEX95uM7dZ2diEopb8w9rHs2MWOQEYJvPG/KPsStMux8fmlWnmwYURbs6\n",
       "jQYC/346osXu5O/zdhz2yrHCDDZfmM5UsqYHA1AJ8GFmWzo+JS+8nhLmCbQYVkmg9Zfrej4MN0yJ\n",
       "rBNX1zzPehWIW/N0tJdlia9lK9dx3l3CgO/b9Ucq+J5x480EAo3VuAMYFp28fNWpdIAzLMzZlJcA\n",
       "4QtHQxMQKqaQEMUIh2nxWraGki0o/g5mrX7nGFmUoxi4QBbQAAABi0GeQniHfwWXItAA/pC1Xc5Z\n",
       "B3dYijBs72VldVvidMOjV6qIwksFFOfcuddASkJNZpm8Vd2n+OfGFqgz4KBaNVXMuxkSloA8bof5\n",
       "ya+ZGoj3u66tWARGilojysprSx3mBiyOE7lpwC3j1catZS0EjjCP2EwWvoqsHHCk/YfDKMzze7Bv\n",
       "UKTDxK+bcEUMtD1oi69xShrTgB/1qb//G5DOJ6PbpLlxkIrterk9AgEX3y7YqtED6qfpP/1lKk3O\n",
       "Ij536Uj5apg/CU2+5ZsjAVj+Na+xEQ2TQsioJtVieTTx5YxgeckzmrVIA9FyNhkiXVaNNrlV3Z3o\n",
       "dpSJ3g9i2xgym7EyL8fBSDVVJUUuCZLrXyLvDsEE6YhYPbfDm0lt5jESVJMxH/5zFyLlZo0nuFN3\n",
       "rsDBYYhZftRBetxz/v9Cl7B1avcrmaJLO4GmJPWLhZKD40R+hHYXlsWAxIH3PuBSiUQLcOgOaKHg\n",
       "SQCYn502b8xg60ru1gT6vGTgofKkubgLt9bVgMLeQOOgjPC5AAABAwGeYXRDfwas4I44cQAbKCzA\n",
       "o8OL7zftKrcdB9MYIZccBvYTraNtTQwXqwNNey56XLhQDMX6McBh7xnsiVqMlUaiN9sehhwq4Pz2\n",
       "eAVcSbz9SwdTDELRMr4B6Nb/WqFxPLrYoRwT7qPsZS5ToPv+SfEUhq3ORj4fA/9jNsVf7X31XUNR\n",
       "nW/ZCHWp7t7tnauUx6o+bxcZW13znaChhwXoAaaNyvL/2/IL6m3smk/xj2tTuWyIRfh7brudpZD2\n",
       "3z83cHw6WficodQ18zrZ0ftMkhGw+HWharXyIK4OSUWq8/+1kBsdfoENceFB74XROQ5MQPunv/YE\n",
       "eE3XbQLgN8D+pZgeWoAAAAEyAZ5jakN/B0fIQACIOJZVPjj9IIg0fiW+SHzpIkFuKsOj1+rjCu5O\n",
       "1LKfds+TYlQnsCS7SSVMe5j2uh1IXL0PGYyYjmJi3+tq7KxQ9gPxxqjqGj11Bn6AAUctC6ccOvwF\n",
       "ony8DqYZOxRyuc2D6mdDXakMZ+EqFbcdYMTfj0bQtE8vRM9zpS01F9c6ucwuieNnXL/fCEk2FrXY\n",
       "AO475TnJb+NRCMf6OS+wcFFooHFQ27CkRcV3xakmdPCd+/lyxZ2ALEph20vA8eAuniNBt5/cvUmt\n",
       "8WKEFnvZrJg7JJu6lKXy9lO1p1t+lGX6ByRVMD2tAjrIO3qyxkQ48VURHCGOEiSAWMlgyk+Xwdvv\n",
       "Gv4ebPxW+ohIA45wtaSuskI8ZnwhNCCb97hpIgEQDvbgJSm1MnExAAACXEGaaEmoQWiZTAgj//61\n",
       "K8IQ4AALVH1Zu9dS1LBcrBrvhicxlUPIZfH9lvWIR0z2w2s5AAScRz3do40TPC3/GwDIbr67mQ6Q\n",
       "KI9DEoyXyABoGZ3NvaTwor+hxtAMaBB2CYqj2i8ztLvVjlLF1CZJnpEOZmLZQXsI0bA1n71u29zj\n",
       "J+6WDuQt3YBwX6vJOGo9E4TiZv3SMSCy8xlCQZhYmsr5vmTyChhiv+OWBYasQv7VGDJngkJlMZP2\n",
       "vwi0FoPz6dpkvb2ukZA/f0xJmO5P6fDnzwmz9eNxyPdFko6QgX/tkCmGruFtPzJ5e+BY9tRSZS/w\n",
       "VfQ+5TWvAQqjPh9THUo3cn71s3z7n7lelPgtQASKdketO9MIBfPHX+ub2wH51qAaCNaESDuNpiW4\n",
       "a+7Qc3SPtsY4lZXf7kaQ48lCv+IHAU7e2/VM6ouSZ55Q/hec/sptEURHfqqdvqb9D51foN2t9Dz5\n",
       "RiFeHFM1R/MMwvKso/zqno1b8LGK1vMegAtELKlJR+IUTPgLaTwEM00p0cXZu+nuHkOYF0Yg0iWP\n",
       "+fSRp+YI9nsKx1cGnkd4t7UKRJusbHDyRp4UtVrjqDM1e3SZp5yTTeuyQlMH2lzd5gwJzNbl+iVW\n",
       "ATvMDLqP7QFMIFmQ7uaAexmB3MUOJrMIng9LxEeWbjk3tLfu2v47WnfvpDkV4WdX+ivBdJ/4IzV9\n",
       "+KqM0JRlFW85y01rD+lj84XXECRb6lhb4aja5uR0sgJIx+RqalehrY+xdtVpVqRV7wjON0NfQyGV\n",
       "EX1SdedjAsfZq8cpHt11Y28gmgcAAAFdQZ6GRREsO/8FlofLACFPPGMVVxW44RRiMCmurhlfetBS\n",
       "5m48YvuA13vzXcJ36VW+SZns15Xkl/Hz9GVgUzTYL/t2RXctW3749xQzrK6h6x7/h+X1Kma/WpS9\n",
       "cO6fO6TT+Vg20uFB+dHZqxXKqXWBFq56Fx6IK7gE50bGLKoTJoYWdUHyhx7mq978FstYGinGYzKE\n",
       "RZRk9NNWt5DFn7+X7s1XYkF2xsK37TY+e6pCgLQdAXyEYwPWvLkVlm6hCQ1RgUTGVSlAQ8Hq5PYK\n",
       "Rbn2i/X8CoLkXqPu0AWRvsmdptm9jdsmgoJWucR4+A6xVqVc2NXRccymclFIidO8nYze37Yn7iX6\n",
       "ykrjnEfnYeSuZtS8UUP+5ZcthULZ1zphTdD7fNxWem6yw2nrLEIKCbh1exwycYqj/99EDhsLrl4C\n",
       "WPM5G5+XRsLuql2mvxA/D7etCOMGpyW3k5dMLLb2wQAAAOEBnqV0Q38HRLdkAG1Vn6rEGcJXSa+t\n",
       "AfBjrLJ60tua/COwdyR5H79aGWpUYumd5bywBmUy6s1FB3+ebFaC1IE2O9SIwJk1NhzuFXOmczXv\n",
       "+bad5aFb6VT7lanj6oa+/jLm7GqQHfaCU/iLXXcDnbeE/C/2IwYNiiayOhFpTdsG38vkr4+NA5xe\n",
       "C3yDn1QoHfZfWUm3jFEhyIZqTtzruDhlXFUSsSTEYzs1fnuJKVoaMv8qj4oO3u6T3FZ+IEfdhGE1\n",
       "g4esETfmN+KmNB7btsqTWowfMxZqQqFMomHxXGfDa2UAAAD/AZ6nakN/BxqoAFogBCnvM5XzMKt6\n",
       "hkHHF0zBUiweI90sWlaTr4IM3TkDEksSWYNz2Ldf8/+e8SoOVqj+tnmuDtmexnFhsxAvVeNa/kZD\n",
       "wVwcL1vfDYLKvXDORg48FdorhILiOmVS+LOUdjfoy7AYezqo9CPdeSOIoyKL9n2w9u4nAN6KDoB3\n",
       "vbL04bep5Ox6kikHtoHjmk9nf4GEWK5YItd6h+qt7Y8EWuIUuxACa+zfO3tKZROo7Gos8uo8Q9tp\n",
       "yMfJCWidSRhvqBU1Kd+w2svWo4EYClEro1u/Kw5xfqTgnVi2RXVC5HjchlTvyvsBwho3YTWvquXi\n",
       "eDBdAmPgAAACKkGarEmoQWyZTAgj//61K7972MZa1vfCXTqRwAnW4KYkNh2q6fQoQl1arjS2RuaI\n",
       "fq64jRHQIzqHt2kn/q3nK5iRFXRrlzpKDLQK5tj9tRxjSX1NSUA0HwTmj10y19/xVjekOv4/OiLG\n",
       "gbgDOo9TX7UWv+pRabI5lzSzUmIK2S5ns5yEmKVaF//2doAp1T8qkhSUL1EznDIfPlhGzPjo9iJM\n",
       "FGv/n+gJpl2OODz/e7tEaArr9OC1M6eyViMjAD6Rfj0cYK58d6cnbP3YyC+eh4kHfWrJC94u8H2d\n",
       "90ohQh5r4OJGKyGqZwwOqIGL1cd3YGmuQMD85FOuoMU+pAhTwmf+r1So08tZNuAs4VwgS9Z6cywG\n",
       "K3Ducy7EH5+fuqsXXKqgvu6A3QZw/bv7JEo1wOTdoaHra77RsBYbnslD8CTXMt117QypN+AEZUrL\n",
       "grLfr9Jg/C6RrcPuRI1WsFGQ0OWGk7NsJfDcYfT3f+76PuofHX7/eUK035rCtAE4Hyao9I6NXy/6\n",
       "OVPr7dMxxVFXam0CNFK77K9lvxRdZoUNZWdKH6IxV4TjH3GhZzW0n/slFFRwPOkVjgm2iAvNCpQ7\n",
       "IQJZBfYH8i3/9Pfzo45ejS1OKRkzsczI6ZOXk/FLRLf6YfQJ/c/viNAxIWm9NSpe7vAw14ADVdcM\n",
       "wCsePLaJziMHyukV/WJGIs6MMlbUKcZxQVt7uASOm+X7mLJ2SiHFclX2jqgXpPCQBHYyAAABPkGe\n",
       "ykUVLDv/BRpnTojMzjwHwANpgpmg2xm4KrVYasaGsZDFY3twR4teMxeHYAgb3D/s6j6iZ6AFwUuC\n",
       "HK4M74Br2vmx8zqklEv+IQXd2kcG6t7vheY5+7+i77O0tA766g2Cs6Jsy6xMNukia4UF71c/hAhv\n",
       "xzFAjS4vfs9pJ7lxOUY/jlNg1/mB1gse+/Ehgk/ciEpmthrf85yNn2CTRbBWD4LGEA8wyH6jCz0l\n",
       "i6zzfoLPFgU6yzkXKSJc6axB5nFhiAE6kaijtn4rxE9q2zi3/OZ8WW9Kpv44DcfydoK9hcIhQnOd\n",
       "dVbm14lE3PFLSBw6UzxKCNJuCZF+6vPh6gG00FWUcdb6B7Wz8080BrVR3BJpc3n8PkRj3V9q6Sdu\n",
       "Is8CKOi3gfF6Wpr1XOr2ErpOR9kMWcxLbk4qcu22wQAAARMBnul0Q38HSmtEAB1xgKXTGie7dxcJ\n",
       "f7UlxD94eUyEFdtDM1oLs0wHDOFUzWSZRWohFInBPjwUp8Q6GIzRprN+hEpjrYGOf0lsamgNOemQ\n",
       "d5S0Z/43d3+1+U1J5D+MjKAh4g2299SxQ/vYeSOqOBP3292IrjdZw9SbJ5eUlx/vbweKWTWUEtb7\n",
       "Xb0zHXb+wcZW9bpe+IINURqTZcDn9u48UTocXOjrJOCvJKHjcPmc+Llyl6aFW9MtE8sSlsFqe49y\n",
       "kV6Hb8GavAPDDlKeq9LIScX9n6BFfQA498UF4ocy1JCCLKojeWEOH19de+h6vP9G/W1ACjNK+R/A\n",
       "MWTfhzezUKHxiItMmpVKUIEu3JVU9Cr2yQAAAOcBnutqQ38GrOG2axvocQAa2cfa6PmgfsrfaO5F\n",
       "4MtPO0YEyzRk4k8seNr1osUqMwIWjfWZPi22GTxU/RSsOiVOHHsRQwBZ2lg+qnRXocYB1gUtwnns\n",
       "d7rzvtkIHbNGx1kmF/uacEkARr055kxNzoWZg1+YCNDvCZCcKukOqX+mPnlRi1Y7LotYvEh0wx+6\n",
       "DjKF8fuwCREMo7/LZKt3N7TGBPxkSyQ5lUsb47OfRf3csEgcqM3Ms/djFhJnt0NH+EYhzbj0QSgw\n",
       "21SL2vuslED8oDRd6aHxcq2MVT127xPn30qEWKiVCcAAAAF1QZrwSahBbJlMCCP//rUr1+dYyAB2\n",
       "lIuozTfSmFsYU3Z5OLhx70PGSk5e6+l5rWSDmRECBwqXSazmnMUoZc1K5fu+b5i3UBD1j0e3fiJL\n",
       "2A5n3J6R327McyKt6BwRoDbof4x+GJ1cbEaf39Dfa7S1fMpc7fWT/fJKGfcQU9L1d4KMi/TOjMbQ\n",
       "XJUrBchoT2j9Ve4QdcDLYO7IDgVk7jfmuzO+tZZJ6li16a0NmhjnqT5ikHjNqGbTU54iC3sjsJ/C\n",
       "eHzDUOwYoCOMVRHWkc3agjipVl6DylBJfCx3gvzn7DvbmAeX17PaTjCyyst/wK+n8lB6AlHVs2Nb\n",
       "9wmHHRUZW3V6UvJLHuDu2Ar76pynpH/jlk/nj5SBYGo1Nt36GOeNeG/lQ2WZM1Q8Y9AQfvzzAyND\n",
       "Ia2QD4MfEySXdf1l4DdqA/T+xgfVC1kFWr9JSO9/twJkzySm005d6L59Q0J1qyht677bjB3Cimqk\n",
       "V1wmLichwpZVQQAAAT9Bnw5FFSw7/wWFL4OwNABeSeR9s2zBYkD/2B1OUHl6poDvgcxI/hvZaurV\n",
       "tm8uFS/bRQ6XuppEE+hY8dRUz9VWwy2LsebixiZD64MoJHGcPA5JIOdhBWYHvsL29Rs2GwWlAbO3\n",
       "6JImiFqGIVw9TcGM7E7rXqBXK7Sa/b7tnuotMrVGY0kjHbvSPp/6Y+YKa3PV39tSfU8EblolrNTa\n",
       "mXt/jD+Mc0TneR9LWlUGEeAe5TqNcI79h3WMRm67saqOrvv/dWLj+0Jy7m924Gz1jP3rHW/UuqLN\n",
       "xabFPlYrcmTPQtjQ6FneCTkb9ttdPGNOKYb3UozDoJ2Rc6xdqj1MBVw3UhnlecHZPAuOl1EBcXpi\n",
       "pNivFY2dptcPbOW2ItiahTrrunXjZO2fJD1akUqPWIGMno7Okuta6zteAXEyoOmxAAABEAGfLXRD\n",
       "fwc34cpAAIgbDg/1rTsveWxka3eD0wKfXRvNdo43oyw8qT6ptdbv2J5ifrz6/jUq7I7VUbVB1HGH\n",
       "zZ5p8f9PsRUXZDpekxtB9Vp/P3Wq2FDlzOtPtQU2EX5HDb1LEn1zZkcGZzjb4cbnGGK14PG1iW5w\n",
       "oVfpJgHpeFCaKs7yv1vguo/wViGN8wErRQYn0MVdi4IGTm8Q0FnwNc3GeXZOjIFUwe+D5P09aEdJ\n",
       "SVOu6BqPkrDV6iVYxv2nWcGTJX5WFTGv0ly4kKJGDVlSO5o+0jXbNZBXN0t2nY8Y2auG6dxjQFsc\n",
       "TMfamuMU5admhIjCdg+kLcxA4NQgm+i4JaPwrEp48/8ubAahU9GBAAABIQGfL2pDfwcm5htkAG1X\n",
       "FTSBQTLZfVXYyMJ/D+jvvV37fc7FJf1sBbNHqyDeS5zj4TMfWS/IdvPtN2S2bR6dhPHDUftuLmee\n",
       "7Y0C4qO7+N6v43lWBrDtgcbkhsJUoHOVoVugV6OxDJI6ecx/tEUUqf+NW6l3qKlnBjo+h9sp0FP7\n",
       "88210wCy7NzvlR4/lCp5CVC5Ot7ZEWMGfsxRSSbLB1ELf0HYinE9xTfrEidzLDc8Nmqjj7Pgmnpz\n",
       "TqVqIMQdOSYL2PI/i2rWNkbb2A2yfWeVYK6NuFdW2XCGxHX41ZAo0lnial+GxT/2TWY7iogenYPO\n",
       "hX+K+sCzCNTKoT4xsD5KsnLrk6kcaCtYEtvYZ17u/AsCcWdVIWOkbFSTXzUrDnAAAAHVQZs0SahB\n",
       "bJlMCCP//rUr679yFjJeJgcgA2s6MXAzspWYFrmxwoZnEZGdYzI6ZKxis+61AxPeirDssF9RarCG\n",
       "5CDEQ07Ws29ZUfG6ljyENm9AX5281X8sGf25hG0rxmGw9UTHuRRdTJTW1EQeqv+Hdo4ugsyywPo7\n",
       "TxKPHnDdlXoxqUB5+C0SrPKSicg+QUHYwwCcmGmwhP2/V6CmDiFSgLqeRpWHgmFD9UVrQ+qARv9k\n",
       "/G5mcAMPVeBM0v0Rnoavapx/Y73XQsFijEbK53dmg/LS2af0SXaNjSjkjCfxfmwlcKd+8vNAMN1W\n",
       "5xWR5VuCYLR09iNwAm/esD2uf5llAyv640UeQz4tSxDRFT5oT52EJea2ma3QT3h9ardERxFbXCJ5\n",
       "Zpk1kNN5sV0e8ydLhKgJpv7o7jY83G1j2/ki1pN6KSHJH/jHNVbVBRlhrvd9mipjmFF2aHBTyHlA\n",
       "eOANGHwZ3/QajdUKNCoxY5AbVE5GbqfmW+sE8aBEJ9rHqIK5XPY+/hVmsTam6l/t52l1ktXXCeoO\n",
       "omqkOCP4D53pBdIrsry+0s2U2V5HCmo3l/0yWUfXU2QCYoNEvgJM5MxOyPOZpaw+pySHFMrbeaCF\n",
       "1+7W/gddmAAAASFBn1JFFSw7/wUZJVxGkjg0AH3As1ILx81Xpbr3zlyD3U2qhrKyrgiXgVAYarco\n",
       "+Igbs+5F9Lddnq7it68njq8ZObs/AKx0IioBUtUF9Xyco+SeJwDC5HYpn+9uIe71bHLX2WGuIM9n\n",
       "mZqPw6Q6EKhEiij2Ebp8IggmgLZVxExKGPTNiXCYH9Nfg9QNpOa1uR2MaOTc77hoATNj+w4yOwwt\n",
       "6iMdaj1rOK972ijuGpmUG5SFf+UHVrlLzc9/qhPTqRKOfMLAJVp4ugj8wP1cW4kWTDJZLTSHVdsq\n",
       "U6GFg8fTHNFkx6w+mwdKmlUDIlE5dHcj3z6MGUoJ3KtMNIx6gzXSjyQ4ItqWNjkuJ+4SL6AJIyjy\n",
       "HIj0Bt+SFBZP8QjoJqfxAAABKAGfcXRDfwas4aaXVCrQACIHg7miSkzRWWKKbDAf4X++5cHz2Ldf\n",
       "8/8VQvjfMbGzzp/7tE2Lv6ZkbMHMy8Gg1c6Xu9tMEYMPjpAskd/LblhKsCMecWzNwoFUMb94RfwL\n",
       "dh4k+5fOEfK7gPPxu6/G0z7YRbznLoyA3/65067hJ5EODjlf0gb9e7CCoKRLOSYaYB/BjoqaurwT\n",
       "zqGi8cYU0/p9n0aH85iSJBgMwVZWp1ymhX59dO8yIkxyXBzOWiJ/G0ZAqcAOvKJiKwKhWJ3sacrQ\n",
       "Z3eUVmaMvjqGE/hD2g26Drt6Rtn9CmkzvRJf/yFp7ITTafFODo1QymTaXXxmj1ywemGhqkDvw3p/\n",
       "d0xPFNEQul8aL3cHzzKHuDelWMPv76kAyxfJEfYWAAABJwGfc2pDfwal3KX0pBKIAEMXf992beF5\n",
       "9KWd/f5jXRQX7Bb4a/EvAgxpwrznN4O1q41OKw5tL6YEgKpEIhCX1o7XTKJWWXqCki4RmXCjd60m\n",
       "KFUkz4DkTe3icjUioRWmVxq6GK/eXjDW81c7mtfMT0s1Ncl+h857V7UTZLBVW6GVt3VOltKXzsPz\n",
       "eZtqMf744WnKPByEyKqvYsTYrRwKrGpWZZnbnQtRii7DKOHAVFV/j1Z5rJDs23u0meBD6PopaSJc\n",
       "olMwGS6oYXfJQiUNoP8nz++1b7NiyHxtt4ce4XSRSfmNpTozoMqCM6nHnl/UNxAp7bTJwWS1OLlq\n",
       "3m0ol6elk5l5TPTF7c8MBeRnPDGb2bL0mB86b54038i2QKfyG1x8SQr6f5AAAAFbQZt4SahBbJlM\n",
       "CCH//qpXhQzAAP6Ra9l497cgsw9X13qALTEJzoiNxmE6mah5TKe0oaZ9ZLP3RaBLdW3eCmuNzOux\n",
       "0OmHEVTbx5U4lJHNE9RXRxtI401SdnlmiQl2q7fqpHVnsllMCgJPZ73BSzJzX7iXXbubtYBWdpqn\n",
       "j3BkJrbm9qOnxyLhrtKeZcXdvUg1tjXxQqYUcnYrYfYb0qm686XeXHxyyohqZQ0jm9ffPB87SDYd\n",
       "7Hx0jfcv9XhSnIc+vS0sxWMWcw1HWNX/RhlWxlcy59xzzi0n9XglJEfg6WbVZMA8T1pPZSGJngyZ\n",
       "BMEDxYIeNVGTlDCtMq3utI4i63AZc3EeCXdeGMMg9ss8OJRWH3dNegcH9FbJkoEhUIrpbqVDjIcw\n",
       "jeBKW9CRk5VzmgughTmPMsGNj+qxUhBZKPeJvA1b3JfjTQ00jOrDUh+jGx+mfa3x4nxS388AAAE7\n",
       "QZ+WRRUsO/8FkkZYATV8NKxv6TFlOqUqsDMQYgbVKPuPUF89LQ4wcuXjsaQKZ2GX1Sbyk21mmlbv\n",
       "lOztemWFK3st4ElDiCE+vkKeVF36qcwManOjUjfSsLQyCpLX1UvNdsIdC519wvpIrF5px10K8v2y\n",
       "Jvux7ezmQITHFTW/Rhr+m9aYfrLdhpFYYXwMKK6HM+11rn6FHjggx3ucZSHwqah25QLUoGpkp2Bp\n",
       "2osUq9NHn8lfmyy+fNdrYs81vMitNv38gM8mX14OTAdY/0drkkomFyiI715Oxw+qmhQvXrSdLywV\n",
       "b88lNs+wFA81MKHpWnfVv7mrfeU1qZEUv/0sdkfQViUHxN8PcAqqTNWaDzoKSK6DgmnW25L72bYT\n",
       "vGy4w2cX8NU3tNEJ0IomqMAsO9MBr7EpyIvP7y2AAAAA5AGftXRDfwas4I44cQAbKCzAOMaip+8s\n",
       "surToal/8sw7FAgas0pDwtwh/utLc8QtqTT0wPsHnggh0fAuN9kzgt+xrj8p3LaMzltwxZJ6RrD8\n",
       "RTnTFUFyO8u2W8z8KZY5/Fet9JVdcKg8T30ggL5CzOtnFRBmFkoYhzczQ38tYazhj06alsZieBMk\n",
       "PV6Y3SyucYwkioNLapxBv/y1D/OVbpQHFmfRfjPVeVk3eoQNpk98q0LLEJfYj3AZ9i5ZB15F+u06\n",
       "JGhZRw2OYRaBmYhKVCsDPUkpPn3dY7wXRuvnkHMAcwRvIQAAATcBn7dqQ38HR8hAAIgZzAzn4gYB\n",
       "8623SHS9KQtZAScftE5hwpYOy4YQ2+xC+9fWwtRp3oe5GSOfl1YMIAnlXR2JVRQ//1VIcqdyHGuw\n",
       "obKmvUI6dfmqvqgY+kCrqydyZtFWNuTC3cAA+O44aoa013qouzkFln0+FRblRTdmTeCZsk95zhDz\n",
       "0pFYUjwj56zAKL4f0oqLVdX+WveHglwm6irGJIGWkpWKaASZX10q+0h3+DtaE9yJXZCYy3G1HdKS\n",
       "fuNFoG0/RKCnGfZHQz+mFEygYubTFqqwxANbm6BlsbaaJgfd2k6E7UHvhrie+bUUjtDcanIPHRfq\n",
       "VjrAlpfijhecH8bN/obndd/wUnnhI9geNVJq1MiSbBHdsURO2RKxSOL434OV8kXv+nv8TU/chnbQ\n",
       "ZYoAwp6bpwAAARdBm7xJqEFsmUwId//+qZxkYywAhTK/GIRe7aeyDTAnb4jz+UEXc3V0TZFmDLJ0\n",
       "PuPT+q+AJLfOAZojAZGgYenKJGHnUkPHXutifP9DG5uOUv5WCieZYQrRe7AKTLRwePZ1p6p7YlMG\n",
       "YvFvIJT+2EdRz7uA5u921PQd+WYVklFLA4IT3ghSiIHNBCoCFJdMITaQ866rVS/s0nXckkmxhbu4\n",
       "XLCg95v89PBi94doS1KHTknHDOp4VLEBdcckvpuLCJWiVH11mF0oqdLBRs6UOXVbJB3s0XA9tjLb\n",
       "yjeo2rbIfVhMHBR2KlmIPUocKGinoSYDElKQDDyZxJMohOLfHtszBy9jJBC/tRBfNVJinzFm2Rf6\n",
       "X06wbsAAAAEHQZ/aRRUsO/8FlofLACFPQ5WKq2+Bpp47VzDajIQyjyPL/iqbePkygLy2ITiJb/3R\n",
       "nmGE899oBn2kiqZyTAteQcMU0uhDggShfOd3ixvZzThaZ/8YVsZPEUPmQNm4vXWdtLigr+oY/69S\n",
       "bIYfYj1m9ArY80pcYZMxh+7PCDaxjHb1CfUMaC1Bkvqk2ul/VF5xrEP/pxoHVav78BvjEP/5D/E5\n",
       "FNLDYkKOUCI+bX+6q1ADsSxPxpb7xAbS97TIYIEC0pT8AH1SEMLRY5UiiT1WkVJCjLg3xVn0PJgL\n",
       "NXkjbaMiQHB9/E/Sa6Dl9QpnVekMAg356le23QCE/ayTogKwagvRI9lNeVEAAADqAZ/5dEN/B0S3\n",
       "ZABtVZ+q0D6Aa8YcY5S4MdW0k/Q6+/hHYnDYIWnTxoUSJ3lmd5brMla3wKM/b53YS3O3csXq//Ll\n",
       "Ff8Wopry+a2etBxc5v791z6U32IHC7XDCyAt2QsyR+8lDk695W5BwsTZ6q0sTag4CzEwvXAP7TFr\n",
       "7V7oOcFf/SQsA2jQ+AJSn0mwxjmk9uVcAfHpeYFcyonHIWprxfOlbADyxzddAQS9L4uXDFHUG2xL\n",
       "6JB1EZ28ni2CsAFPlKMqZTi8fo+GE8PIAXxNMu17hDiXNpUCy43OP9GTfoVnAPr3fCCTvC+4AAAB\n",
       "WQGf+2pDfwcaZdZABtVyXdXm6Fv8LfPJuZDSkWDmoaQ44JZBFti3yRGPt7yGvSQ46v/P/vX0RhGV\n",
       "yr15ABFqPs1FoeBtiEgy9rACYgW5q4fbAQaYPZBRq4KOoPaqbSZK58LPfl6oiC6aYexxkrl4mBV1\n",
       "cn0+LmOub1tbmexnNih/FKnOCZ1vbY8nr8uZu8nnD+n03MBfuVGXpDrOJsodk7S8uOTIq4sfox34\n",
       "Jwj95tOkkiF9VQqJduCGBYrQ21SuF17mBFrNeomoF2QvCu3PVrPUZA+kK1ti71Gn4HzUzFWVYmMI\n",
       "rlX2+dIkqwrnjwFc3k5H1U6xpgMh0MBizq2VOPDd9TBhpGSE7HyGT5qwpp3+6dDbwYz4fmH4Rkwz\n",
       "wdxNZPC5u9U7GhEBcZn0xyEvgxO9CwuQ7s1mAByJ3pSbrSrXu1dUtTP/yjv3Cn+wyPzlL5KBoz89\n",
       "X+zUwQAAASVBm/1JqEFsmUwIb//+p45M5DgBB9df5/Xe4nho442so5GjRobSztkGrhaONI4nvaOf\n",
       "9cb7p2+mz9JSQ6n/MTRYJRKCBTUTjTle50nR3QKc+SUh3aS16Hi8MX9TBt3Ecin1w12zR2+oMB/J\n",
       "9QYybv93CTgIMPRR2+928xsqj1YR7xkzVE0rOboDnLFhtNMbspD8WiL2Msqz0d5ZNFKAGWNBMwRL\n",
       "/N8n0S+gnJKBHOu2ajPVou5afGIsbDKEfTPrJ8PsagRw3nDZiVIJlepPPtd3a+JrPynr1guslIBt\n",
       "CyQGEWcHTDqJfYB1rf1kyZ8stBuu24Pv1ODkgW+OA9f5K0fR22BPBeLuSpbxiuSaJra9o4fTq4Tu\n",
       "VlgFXpPjjLHI8pUlzue0gQAABJZtb292AAAAbG12aGQAAAAAAAAAAAAAAAAAAAPoAAALuAABAAAB\n",
       "AAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAA\n",
       "AAAAAAAAAAAAAAAAAAAAAAACAAADwHRyYWsAAABcdGtoZAAAAAMAAAAAAAAAAAAAAAEAAAAAAAAL\n",
       "uAAAAAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAAAAAAAAAAAAEAAAAAAAAAAAAAAAAAAEAAAAABsAAA\n",
       "ASAAAAAAACRlZHRzAAAAHGVsc3QAAAAAAAAAAQAAC7gAAAgAAAEAAAAAAzhtZGlhAAAAIG1kaGQA\n",
       "AAAAAAAAAAAAAAAAACgAAAB4AFXEAAAAAAAtaGRscgAAAAAAAAAAdmlkZQAAAAAAAAAAAAAAAFZp\n",
       "ZGVvSGFuZGxlcgAAAALjbWluZgAAABR2bWhkAAAAAQAAAAAAAAAAAAAAJGRpbmYAAAAcZHJlZgAA\n",
       "AAAAAAABAAAADHVybCAAAAABAAACo3N0YmwAAACzc3RzZAAAAAAAAAABAAAAo2F2YzEAAAAAAAAA\n",
       "AQAAAAAAAAAAAAAAAAAAAAABsAEgAEgAAABIAAAAAAAAAAEAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\n",
       "AAAAAAAAAAAAAAAY//8AAAAxYXZjQwFkABX/4QAYZ2QAFazZQbCWhAAAAwAEAAADAFA8WLZYAQAG\n",
       "aOvjyyLAAAAAHHV1aWRraEDyXyRPxbo5pRvPAyPzAAAAAAAAABhzdHRzAAAAAAAAAAEAAAAeAAAE\n",
       "AAAAABRzdHNzAAAAAAAAAAEAAAABAAABAGN0dHMAAAAAAAAAHgAAAAEAAAgAAAAAAQAAFAAAAAAB\n",
       "AAAIAAAAAAEAAAAAAAAAAQAABAAAAAABAAAUAAAAAAEAAAgAAAAAAQAAAAAAAAABAAAEAAAAAAEA\n",
       "ABQAAAAAAQAACAAAAAABAAAAAAAAAAEAAAQAAAAAAQAAFAAAAAABAAAIAAAAAAEAAAAAAAAAAQAA\n",
       "BAAAAAABAAAUAAAAAAEAAAgAAAAAAQAAAAAAAAABAAAEAAAAAAEAABQAAAAAAQAACAAAAAABAAAA\n",
       "AAAAAAEAAAQAAAAAAQAAFAAAAAABAAAIAAAAAAEAAAAAAAAAAQAABAAAAAABAAAIAAAAABxzdHNj\n",
       "AAAAAAAAAAEAAAABAAAAHgAAAAEAAACMc3RzegAAAAAAAAAAAAAAHgAAEAoAAALbAAABjwAAAQcA\n",
       "AAE2AAACYAAAAWEAAADlAAABAwAAAi4AAAFCAAABFwAAAOsAAAF5AAABQwAAARQAAAElAAAB2QAA\n",
       "ASUAAAEsAAABKwAAAV8AAAE/AAAA6AAAATsAAAEbAAABCwAAAO4AAAFdAAABKQAAABRzdGNvAAAA\n",
       "AAAAAAEAAAAsAAAAYnVkdGEAAABabWV0YQAAAAAAAAAhaGRscgAAAAAAAAAAbWRpcmFwcGwAAAAA\n",
       "AAAAAAAAAAAtaWxzdAAAACWpdG9vAAAAHWRhdGEAAAABAAAAAExhdmY1Ni40MC4xMDE=\n",
       "\">\n",
       "  Your browser does not support the video tag.\n",
       "</video>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create an animation of the traffic density.\n",
    "anim = animation.FuncAnimation(fig, update_plot,\n",
    "                               frames=nt, fargs=(rho_hist,),\n",
    "                               interval=100)\n",
    "# Display the video.\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This MUSCL scheme does not show any of the oscillations you might see with MacCormack or Lax-Wendroff, but the features are not as sharp. Using the _minmod_ slopes led to some smearing of the shock, which motivated many researchers to investigate other options. Bucketloads of so-called _shock-capturing_ schemes exist and whole books are written on this topic. Some people dedicate their lives to developing numerical methods for hyperbolic equations!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Challenge task"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Go back to Sod! Calculate the shock-tube problem using the MUSCL scheme and compare with your previous results. What do think?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Godunov, S.K. (1959), \"A difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics,\" _Math. Sbornik_, Vol. 47, pp. 271–306.\n",
    "\n",
    "* van Leer, Bram (1979), \"Towards the ultimate conservative difference scheme, V. A second-order sequel to Godunov's method,\" _J. Comput. Phys._, Vol. 32, pp. 101–136\n",
    "\n",
    "* van Leer, B., J.L. Thomas, P.L. Roe, R.W. Newsome (1987). \"A comparison of numerical flux formulas for the Euler and Navier-Stokes equations,\" AIAA paper 87-1104 // [PDF from umich.edu](http://deepblue.lib.umich.edu/bitstream/handle/2027.42/76365/AIAA-1987-1104-891.pdf), checked 11/01/14.\n",
    "\n",
    "* Shu, Chi-Wang and Osher, Stanley (1988). \"Efficient implementation of essentially non-oscillatory shock-capturing schemes,\" _J. Comput. Phys._, Vol. 77, pp. 439–471 // [PDF from NASA Tech. Report server](http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19870013797.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "###### The cell below loads the style of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "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",
       "<link href='https://fonts.googleapis.com/css?family=Source+Code+Pro' 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: \"Source Code Pro\";\n",
       "    font-size: 90%;\n",
       "}\n",
       "\n",
       ".alert-box {\n",
       "    padding:10px 10px 10px 36px;\n",
       "    margin:5px;\n",
       "}\n",
       "\n",
       ".success {\n",
       "    color:#666600;\n",
       "    background:rgb(240,242,229);\n",
       "}\n",
       "</style>\n",
       "\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",
       "    MathJax.Hub.Queue(\n",
       "                      [\"resetEquationNumbers\", MathJax.InputJax.TeX],\n",
       "                      [\"PreProcess\", MathJax.Hub],\n",
       "                      [\"Reprocess\", MathJax.Hub]\n",
       "                     );\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 15,
     "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 3 (MOOC)",
   "language": "python",
   "name": "py36-mooc"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}