{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Roots of the depth \u2013 specific energy function"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In free surface flows, the specific energy is a fundamental quantity governing the physics of relevant phenomena, e.g, [Chow (2009)](http://blackburnpress.stores.yahoo.net/ophy.html \"Open-Channel Hydraulics\"). Recently, Computational Community (e.g [Le Floch & Duc Thanh (2011)](http://dx.doi.org/10.1016/j.jcp.2011.06.017 \"A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime\"), [Castro et al. (2013)](http://dx.doi.org/10.1016/j.jcp.2013.03.033 \"High order exactly well-balanced numerical methods for shallow water systems\"), [Bernetti et al. (2008)](http://dx.doi.org/10.1016/j.jcp.2007.11.033 \"Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry\")) has re-discovered its crucial role in solving some Shallow Water problems on discontinuous bottom.\n",
      "In the work of [Valiani & Caleffi (2008)](http://dx.doi.org/10.1016/j.advwatres.2007.09.007 \"Depth-energy and depth-force relationships in open channel flows: Analytical findings\") the analytical inversion of the depth - specific energy relationship is given.\n",
      "\n",
      "-------------------------------\n",
      "\n",
      "**The reference:**\n",
      "\n",
      "Valiani, A. & Caleffi, V., [Depth-energy and depth-force relationships in open channel flows: Analytical findings](http://dx.doi.org/10.1016/j.advwatres.2007.09.007 \"Depth-energy and depth-force relationships in open channel flows: Analytical findings\")\n",
      "(2008) Advances in Water Resources, 31 (3), pp. 447-454.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "--------------------------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# preparation of the numerical environment\n",
      "%matplotlib inline\n",
      "from numpy import sqrt, cos, pi, arctan, linspace\n",
      "from matplotlib import pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The problem to solve"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For a fixed value of the specific discharge $q$, the classical expressions of the specific energy $E$ as functions of the depth $Y$ is:\n",
      "\n",
      "$$\n",
      " E = Y + \\frac{{q^2 }}{{2g\\,Y^2 }};\n",
      "$$\n",
      "\n",
      "where $g$ the is gravity acceleration, and can be interpreted as a third-order equation in $Y$.\n",
      "\n",
      "The problem to solve is to *analytically* found the two *pysically meaningful* depths $Y_{sb}$ and $Y_{sp}$ (subcritical and supercritical, respectivelly) for given values of $E=E_0$ and $q=q_0$."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The analytical solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The analytical solution described in [Valiani & Caleffi (2008)](http://dx.doi.org/10.1016/j.advwatres.2007.09.007 \"Depth-energy and depth-force relationships in open channel flows: Analytical findings\") is given by the following steps:\n",
      "\n",
      "1. for the given discharge $q_0$, the critical values of the specific energy $E_c$ and depth $Y_c$ are:\n",
      "$$\n",
      " Y_c  = \\sqrt[3]{{\\frac{{q_0^2 }}{g}}};\\quad E_c  = \\frac{3}{2}Y_c;\n",
      "$$\n",
      "\n",
      "2. the non-domensional depth, $\\eta$, and the non-domensional energy, $\\Gamma$, are obtained using their critical values:\n",
      "$$\n",
      "\\eta  = \\frac{Y}{{Y_c }};\\quad \\Gamma  = \\frac{E}{{E_c }};\n",
      "$$\n",
      "\n",
      "3. the non-domensional form of the depth - specific energy relationship becomes:\n",
      "$$\n",
      "\\Gamma = \\frac{2}{3}\\eta + \\frac{1}{3}\\left( {\\frac{1}{\\eta }} \\right)^2;\n",
      "$$\n",
      "\n",
      "4. this equation has two meaningful solutions (see [Valiani & Caleffi (2008)](http://dx.doi.org/10.1016/j.advwatres.2007.09.007 \"Depth-energy and depth-force relationships in open channel flows: Analytical findings\") for the details):\n",
      "$$\n",
      "\\eta _{sb}  = \\frac{{\\Gamma _0 }}{2}\\left[ {1 + 2\\cos \\left( {\\frac{\\pi -2 \\alpha}{3}} \\right)} \\right]; \\quad \\text{subcritical depth}\n",
      "$$\n",
      "$$\n",
      "\\eta _{sp}  = \\frac{{\\Gamma _0 }}{2}\\;\\left[ {1 + 2\\cos \\left( {\\frac{\\pi+2\\alpha}{3}} \\right)} \\right]; \\quad \\text{supercritical depth}\n",
      "$$\n",
      "where $\\alpha  = \\arctan \\left( {\\sqrt {\\Gamma _0^3  - 1} } \\right)$;\n",
      "\n",
      "5. corresponding to the dimensional form:\n",
      "$$\n",
      "Y_{sb}  = \\eta _{sb}Y_c;\n",
      "\\qquad\n",
      "Y_{sp}  = \\eta _{sp}Y_c.\n",
      "$$\n",
      "\n",
      "\n",
      "----------------------------"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The following *python code* is proposed to highlight the simplicity of the method."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def radix(q_0,E_0,g):\n",
      "    \"\"\"\n",
      "    Compute analitically the roots (Y_sb, Y_sp) of the equation:\n",
      "    E_0 - Y - q_0^2/(2 g Y^2) = 0\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    Input:                             |  Output:\n",
      "        q_0 : float                    |     Y_sb : float\n",
      "          [m^2/2] specific discarge    |       [m] subcritical depth\n",
      "        E_0 : float                    |     Y_sp : float\n",
      "          [m] specific energy          |       [m] supercritical depth\n",
      "        g : float                      |\n",
      "          [m/s^2] gravity              |\n",
      "    \"\"\"\n",
      "    \n",
      "    Y_c = (q_0**2/g)**(1.0/3.0)   # critical depth\n",
      "    E_c = 3.0/2.0*Y_c             # crtitical energy\n",
      "\n",
      "    # check the esistence of physical solutions!\n",
      "    if E_0 > E_c:\n",
      "        Gamma_0 = E_0/E_c\n",
      "    else:\n",
      "        print('no physical solutions!')\n",
      "        return 0.0, 0.0\n",
      "\n",
      "    alpha = arctan(sqrt(Gamma_0**3 - 1.0))\n",
      "\n",
      "    # non-dimensional subcritical depth\n",
      "    eta_sb = Gamma_0 / 2.0 * (1.0 + 2.0*cos(( pi - 2.0*alpha)/3.0))\n",
      "    # non-dimensional supercritical depth\n",
      "    eta_sp = Gamma_0 / 2.0 * (1.0 + 2.0*cos(( pi + 2.0*alpha)/3.0))\n",
      "\n",
      "    Y_sb = eta_sb * Y_c # subcritical depth\n",
      "    Y_sp = eta_sp * Y_c # supercritical depth\n",
      "\n",
      "    return Y_sb, Y_sp "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Test of the solution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To test the solution we find the two physical roots of the specific energy function for $q$ = 2.00 [m$^2$/2] and $E$ = 2.50 [m]:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gg = 9.81    # [m/s^2] gravity\n",
      "qq = 2.00    # [m^2/2] specific discarge\n",
      "EE = 2.50    # [m]     specific energy\n",
      "\n",
      "Y_ss = radix(qq,EE,gg)\n",
      "    \n",
      "Y = linspace(0.2,3,100)\n",
      "E = Y + qq**2 / (2.0*gg*Y**2)\n",
      "    \n",
      "plt.figure(2,figsize=(10, 7))\n",
      "plt.title('Speficic Energy vs. Depth')\n",
      "plt.plot(E,Y,label='E = E(Y)')\n",
      "plt.plot((EE,EE),Y_ss,'or',label='solutions')\n",
      "plt.xlabel('E [m]')\n",
      "plt.ylabel('Y [m]')\n",
      "plt.axis((1.0,4.0,0.0,3.0))\n",
      "plt.grid('on')\n",
      "plt.legend()\n",
      "plt.show()    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAHBCAYAAAAo6sxCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcjXX/x/HXMURjm7EzM8iSpcVSUVGmRJZEmyyJcpeU\nre67u5/ipu5UukMpFS3cJVruCGFKaUJuTRhtyNZkFsluFDNj5vr9cd1zmjFnVud8r7O8n4/HeXCd\nc805n3k38XFdn+t7uSzLshARERERR5VzugARERERUVMmIiIi4hfUlImIiIj4ATVlIiIiIn5ATZmI\niIiIH1BTJiIiIuIH1JSJSIlMmDCB2rVr06BBA5KTk6latSrFraizdu1aWrZsaahC8ZbJkyczZMgQ\np8sQCTlqykSCyLp167jyyiuJiIigZs2adO7cmY0bN571++7du5fp06ezfft20tLSiImJIT09HZfL\nVeTXXXXVVWzfvr3Unzd58mQqVKhA1apV3Y8aNWqUtfyAlptFtWrVqFatGi1atGD06NH8+uuvXnn/\n+Ph4YmJi8j1X3H9XEfENNWUiQeL48ePccMMNjB07liNHjpCamsqkSZOoWLHiWb/33r17qVmzJjVr\n1vRCpcVzuVwMHDiQ9PR09+Pw4cNe/5zTp097/T29LTeL48ePc+TIERYvXsyvv/7KJZdc4rXG7Exa\nU1zEGWrKRILEjh07cLlc3H777bhcLipVqkS3bt246KKLAJg3bx6dOnVi9OjRRERE0KpVK1avXu3+\n+mPHjjF8+HAaNGhAdHQ0EydOJCcnh88++4zu3buTlpZG1apVufvuu0lKSqJcuXLk5OQAcPjwYe66\n6y6ioqKoUaMGN910E1DwKExycjI333wzderUoVatWowePdrj92JZVpGNQbly5Zg9ezbnn38+kZGR\njBo1Kt/rb775Jq1bt6ZGjRr06NGDvXv35vval19+mebNm9OiRQsAnn32Wff3/frrr1OuXDn27NnD\nN998Q7169fLVsmjRItq2bVugpq+//pr69evn23fx4sW0adMGgISEBC699FKqV69OvXr1+Otf/1ro\n91dYFmFhYbRu3Zr33nuP2rVrM23aNPd+H3/8MW3btiUyMpJOnTrx/fffu19r3LgxzzzzDBdccAE1\natTg7rvvJiMjg99//52ePXu6/9tWq1aNffv24XK5yMzMZOjQoVSrVo0LL7yQTZs2laheESk7NWUi\nQaJFixaEhYUxbNgw4uLiOHLkSIF9EhISaNasGYcOHeLxxx/n5ptv5ujRowAMGzaMc845h927d5OY\nmMinn37K66+/znXXXcfKlStp0KAB6enpvPnmmwXed8iQIZw6dYqtW7fy22+/8dBDDxXYJzs7mxtu\nuIHzzjuPX375hdTUVAYMGFDm73f58uVs3LiR7777jvfff59PPvkEgCVLlvD000+zePFiDh48yFVX\nXcXAgQPzfe2SJUv45ptv2Lp1K3FxccyYMYPPP/+cnTt3Eh8f797vsssuo2bNmu73Bnj77bcZOnRo\ngXo6duxI5cqV+fzzz93PLViwgMGDBwMwduxYHnzwQY4dO8aePXvo379/mb/3cuXK0bdvX9auXQtA\nYmIiw4cP57XXXuPw4cOMGDGCG2+8kaysrHy1fPrpp+zevZsdO3bw5JNPUrlyZeLi4tz/bY8fP+5u\nLJcuXcrAgQM5duwYN954Y4HGV0R8wBKRoLFt2zZr2LBhVnR0tFW+fHnrxhtvtPbv329ZlmXNnTvX\natCgQb79O3ToYL399tvWr7/+alWsWNE6efKk+7UFCxZY11xzjWVZlvXFF19Y0dHR7td+/vlny+Vy\nWdnZ2VZaWppVrlw56+jRowXqyft169evt2rXrm1lZ2cX+31MmjTJOuecc6yIiAj349prr3W/7nK5\nrK+++sq93b9/f2vq1KmWZVlWjx49rDfeeMP9WnZ2thUeHm7t3bvX/bVffPGF+/W77rrLevTRR93b\nu3btslwul7V7927LsizrmWeesQYPHmxZlmUdOnTICg8Pt3799VePdU+YMMG6++67LcuyrOPHj1uV\nK1d2f+7VV19tTZo0yTpw4ECx3/+ZWdxxxx0Fnn/llVes5s2bW5ZlWffdd581ceLEfK+3aNHCWrNm\njWVZltW4cWNr9uzZ7tdWrFhhNW3a1LKsgv9tcz+zW7du7u0ff/zROvfcc0tVt4iUno6UiQSRli1b\nMnfuXJKTk/nhhx9IS0tj3Lhx7tejoqLy7d+oUSPS0tLYu3cvWVlZ1K9fn8jISCIjI7nvvvs4cOBA\nsZ+ZnJxMjRo1qF69erH7NWrUiHLlSvbHzu23386RI0fcj7xHoADq1avn/n14eDgnTpwA4JdffmHs\n2LHu7yN3Di41NdW9f95Tqvv27cu3HR0dne9zBg8ezLJly/jjjz94//33ufrqq6lbt67HmgcNGsSi\nRYvIzMxk0aJFXHLJJe73fuONN9ixYwetWrWiQ4cOLF++vEQ5FCY1NdX9vf3yyy9MmzbN/T1HRkaS\nkpJCWlqax++5YcOG+V7zJO/3GB4ezqlTp9ynq0XEN8o7XYCI+EaLFi0YOnQoc+bMcT+XtzEB+y/z\nvn37EhMTQ8WKFTl06FCJm6ZcMTExHD58mGPHjhXZmMXExLB3716ys7MJCwsr8j1dLleZh80bNmzI\nxIkTC5yyPPP9c9WvX5/k5GT3dt7fg92kXX755SxatIj58+dz//33F/q+rVq1olGjRqxcuZIFCxYw\naNAg92vNmjVjwYIFAHz44YfceuutHD58mHPPPbfI78fTlZA5OTksW7aM7t27u7/nxx57jEcffbTQ\n98k7V7d3714aNGhQ6Pvr6ksRZ+hImUiQ+Omnn5g+fbq78UpOTmbhwoVcccUV7n1+++03Zs6cSVZW\nFh988AHbt2+nV69e1KtXj+7du/PQQw+Rnp5OTk4Ou3fvZs2aNcV+bv369enZsyf3338/R48eJSsr\ny+PXdejQgfr16/N///d//PHHH5w6dYr169d7fM/SNmRWnmH4++67j6eeeoqtW7cC9gUMH3zwQaFf\n279/f+bOncv27dv5448/+Oc//1lgnzvvvJOpU6fyww8/cPPNNxdZy6BBg3j++edZu3Ytt912m/v5\n+fPnu488Vq9eHZfLVaIGOG8Wp0+fZtu2bQwcODDf7N4999zDq6++SkJCApZl8fvvv7N8+XL30UPL\nsnj55ZdJTU3l8OHDTJkyxT3PV7duXQ4dOsTx48c9fqaImKOmTCRIVK1ala+//pqOHTtSpUoVrrji\nCi6++OJ8V+h17NiRnTt3Urt2bSZOnMiHH35IZGQkAG+99RaZmZnuqxZvu+22fEsunHn0JO/222+/\nTYUKFWjZsiV169Zl5syZBfYLCwtj2bJl7Nq1i4YNGxITE8P777/v8XtxuVy89957+dYpq1atGgcP\nHiy0ltzn+vXrxyOPPMKAAQOoXr06F110Ub5B/TO/tkePHowZM4ZrrrmG888/393E5l1K5Oabb2bv\n3r3cdNNNVKpUyWPNuQYOHMiaNWvo2rVrvrXVPvnkEy688EKqVq3Kgw8+yLvvvuv+jKpVq/LVV18V\nm0VERAR9+/aldu3abNq0yX0K95JLLuG1115j1KhR1KhRg+bNm/PWW2+5v1eXy8WgQYPo3r07TZs2\npXnz5kyYMAGwT3kPHDiQJk2aUKNGDffVl0X99xYR33BZPvon0alTp+jSpQsZGRlkZmbSt29fnn76\n6QL7jRkzhpUrVxIeHs68efNo166dL8oRCXnz5s3jjTfecF+xJ55t27aNiy66iMzMzHxHspo3b87s\n2bO59tprHayubM477zzeeOONgKxdJJT47EhZpUqV+OKLL9iyZQvfffcdX3zxBevWrcu3z4oVK9i1\naxc7d+5kzpw5jBw50lfliIgUavHixWRkZHDkyBEeeeQRbrzxxnwN2aJFi3C5XGpqRMSnfHr6Mjw8\nHIDMzEyys7ML3CZl6dKl7vV+OnbsyNGjR9m/f78vSxIJWZ5OSYltzpw51K1bl2bNmlGhQgVeeeUV\n92uxsbHcf//9zJo1y8EKRSQU+PTqy5ycHNq3b8/u3bsZOXIkrVu3zvd6ampqgUvRU1JSCr3cXETK\nbujQoR4XPRVYuXJloa/lXUw2UP38889OlyAiJeDTI2XlypVjy5YtpKSksGbNGo9/uJ050qZ/yYuI\niEgoMrJOWfXq1enduzcbN24kNjbW/XxUVFS+NYFSUlIKLG6Zu19xCx2KiIiI+IOmTZuya9euUn+d\nz46UHTx40H1PvZMnT7Jq1aoCV1beeOONvPXWWwBs2LCBiIgIj6cu09LS3OsQ6WHmMWnSJMdrCLWH\nMlfmofAYP34Sw4dbNGtm8fXXztcTCg/9nJt/7N69u0y9k8+OlO3bt4+hQ4eSk5NDTk4OQ4YMoWvX\nrsyePRuAESNG0KtXL1asWEGzZs2oXLkyc+fO9VU5UkpJSUlOlxBylLl5ytysb76BF19Mon9/SEyE\nKlWcrig06Oc8cPisKbvooovYvHlzgedHjBiRb/ull17yVQkiIuIHsrNh6lR44QW45BJ44w2nKxLx\nT7r3pXg0bNgwp0sIOcrcPGXue7/8AkOGQFgYbNwIu3cPc7qkkKOf88DhsxX9velsbk4sIiLOWLgQ\nxo6Fv/0N/vpXuzETCQVl7Vt070vxKBjWZgo0ytw8Ze4bx47ZR8cefxxWroS///3PhkyZe1eNGjXc\nC0PrYf5x5qL4Z0unL0VExGvWr4c77oDu3WHTJqhc2emKgtuRI0d0JslB3l5bVacvRUTkrJ0+Df/8\nJ8yebT/69nW6otCgvx+dVVj+Zf3voiNlIiJyVvbsgcGDoWpVe6mL+vWdrkgkMGmmTDzS3Id5ytw8\nZX52LAveegs6doT+/SEurviGTJmLFE5NmYiIlNrRozBwoL3+2GefwYMPQjn9jSJ+qFOnTnz77bfF\n7rd//35at25NZmamgao80/9C4lHee5SKGcrcPGVeNmvWQJs2ULu2vfZYmzYl/1plHjoaN25MeHg4\nVatWdT/GjBnj1c8YNmwYFStWzPcZeW/puGzZMqpXr06bNm1YtWoVdevW5dChQ+7XMzIyaNWqFXPm\nzKFu3bpcc801zJkzx6s1loaaMhERKZGsLHjsMbj9dnj5ZXjxRTj3XKerEn/lcrn4+OOPSU9Pdz9m\nzpzp9c945JFH8n1GYmKi+/VXX32VIUOGANCtWzf69OnD2LFj3a8/+eSTREVFce+99wIwePBg9+0g\nnaCmTDzS3Id5ytw8ZV5yO3dC586wZYv96N27bO+jzMWUzMxMvvjiC7p06eJ+bvr06cTHx7NixQp+\n+OEHZs2axeuvv+5+vUOHDuzZs4fk5GQnSlZTJiIihbMsePNNuPJKe0HYjz+GunWdrkoCRUmXhViw\nYAGRkZEeHzVq1CAlJaXUn7Fz507KlStHgwYN3M9Vq1aNV199lREjRjB8+HAmT55M48aN3a+XL1+e\nZs2asWXLlpJ9g16mdcpERMSjw4fh3nthxw77lkkXXOB0RXImf/77sXHjxhw6dIjy5f9cfeu5555j\n+PDhXvuMYcOG8d5771GpUiX3c/369WPu3Ll89dVX3Hrrrezbt6/A1/Xv35+kpCQSEhIKvNa5c2fu\nu+8+7rjjjmI/X+uUiYiIz61eDUOHwq23wvz5kOfvPAkg3lpwvix9n8vlYsmSJVx77bXeKaKQz3j4\n4Yd54oknCrwWGRlJenq6x6+74IILqFixosfX0tPTiYiI8GqdJaXTl+KR5j7MU+bmKfOCMjPhkUfs\nWyW99hrMmOHdhkyZm2VZ3nn42jvvvJPvCsq8j2rVqhV5+rIwzZo1w7Isj0fKCjuKdfr0aXbt2kWb\n0lxS7EVqykREBIDt2+GKK2DbNvj2W+jRw+mKJNCV9BTe4MGD811Bmfdx/PhxoqOjC33/wj7jnHPO\n4brrrivVPwQSEhJo3LgxMTExJf4ab1JTJh5pLSHzlLl5ytxmWTBnjn115V/+AkuW2GuQ+YIyDy19\n+vTJd9Trlltu8er7u1wunn322XyfUadOHffrI0aM4O233/b4dZ5uJv7OO+8wcuRIr9ZYGhr0FxEJ\nYQcPwj33wM8/28P8rVo5XZGUhv5+LF7nzp2ZNWtWsackf/vtN2JjY9myZQvnnHNOid7b24P+OlIm\nHmnuwzxlbl6oZ/7ZZ9C2LTRpAl9/baYhC/XMxbx169aVaEasTp06bN26tcQNmS/o6ksRkRCTkWGv\nzL9wIcybB926OV2RiIBOX4qIhJRt22DQIGjUCF5/HWrVcroiORv6+9FZOn0pIiKlZlnw6qtw9dUw\nciQsXqyGTMTfqCkTjzT3YZ4yNy9UMj94EPr1s6+wXLvWXqXfW4uKllaoZC5SFmrKRESC2KpV0KYN\ntGgBGzZAy5ZOVyQihdFMmYhIEMrIgEcfhffeg3//G7p2dboi8QX9/egs3ftSRESKlDvMf9559sr8\nNWs6XZGIlIROX4pHmvswT5mbF2yZWxa88oo9zH///fDhh/7XkAVb5uI98fHxZ3V7o6effpp77rnH\nixWZpyNlIiJB4MABGD4cUlNh3Tp7hkwkWMXHxzNkyBCSk5Pdz40fP97BirxDR8rEI92fzjxlbl6w\nZP7pp/bK/K1awX//698NWbBkLuILaspERALUqVPw4IP2EbK334apU8HBO8SIH1mzfDkTrr+eybGx\nTLj+etYsX278PaZOnUp0dDTVqlWjZcuWrF69mszMTMaNG0dUVBRRUVE8+OCDZGZmevz6cuXKsWfP\nHvf2sGHDmDhxIn/88Qc9e/YkLS2NqlWrUq1aNfbt28fkyZMZMmSIe/+lS5dywQUXEBkZyTXXXMP2\n7dvdrzVu3Jhp06bRpk0bIiIiGDBgABkZGQAcPHiQG264gcjISGrWrMnVV19t7GIKnb4Uj+Lj4/Uv\nWsOUuXmBnPnWrTBwIDRtClu2+N/sWGECOfNAsWb5cj4ZO5Ypu3e7n3vsf7+/undvI+/x008/MWvW\nLDZu3Ei9evXYu3cvp0+f5sknnyQhIYFvv/0WgL59+/Lkk0/yxBNPFPueLpcLl8tFeHg4cXFx3HHH\nHflOX7ryLL63Y8cOBg0axJIlS4iNjWX69On06dOHbdu2Ub58eVwuFx988AGffPIJFStWpFOnTsyb\nN48RI0Ywbdo0YmJiOHjwIAAbNmzI996+pCNlIiIBxLLg5ZehSxcYPdo/h/nFWZ/OnJmvmQKYsns3\nq1580dh7hIWFkZGRwY8//khWVhYNGzakSZMmLFiwgH/84x/UqlWLWrVqMWnSJN5+++0S15V7xMrT\nkau8z7333nvccMMNdO3albCwMP72t79x8uRJ1q9f795nzJgx1KtXj8jISPr06cOWLVsAOOecc9i3\nbx9JSUmEhYXRqVOnEtd3ttSUiUf6l6x5yty8QMv8wAHo2xfefBO++gr+8hfnVuYvq0DLPBCV/99p\nuDOFnTpl7D2aNWvG888/z+TJk6lbty4DBw4kLS2NtLQ0GjVq5N6vYcOGpKWllbiukkpLS6Nhw4bu\nbZfLRUxMDKmpqe7n6tWr5/79ueeey4kTJwB4+OGHadasGd27d6dp06ZMnTrV6/UVRk2ZiEgAyDvM\nv349nH++0xWJvzpdsaLH57MrVTL6HgMHDmTt2rX88ssvuFwuHnnkERo0aEBSUpJ7n71799KgQQOP\nXx8eHs4ff/zh3t63b5/7NGJxpxOjoqL45Zdf3NuWZZGcnExUVJTH/fO+X5UqVXjuuefYvXs3S5cu\nZfr06axevbrY79cb1JSJR1pLyDxlbl4gZJ6RAQ89ZA/zz58f+MP8gZB5oOs+ZgyPNW2a77lHmzal\n2+jRxt5jx44drF69moyMDCpWrEilSpUoX748AwcO5Mknn+TgwYMcPHiQJ554It9wfl5t27blnXfe\nITs7m7i4ONasWeN+rW7duhw6dIjjx497/NrbbruN5cuXs3r1arKyspg2bRqVKlXiyiuv9Lh/3lOf\nH3/8Mbt27cKyLKpVq0ZYWBhhYWEl+r7Plgb9RUT81Nat9sr8gTbML87KHcSf+OKLhJ06RXalSvQY\nPbrEQ/7eeI+MjAzGjx/Ptm3bqFChAp06dWLOnDlERkZy/PhxLr74YgD69+/PhAkT3F+X94jVCy+8\nwNChQ5k1axb9+vXjpptucr/WsmVLBg4cSJMmTcjJyeHHH390XwgA0KJFC+bPn8/o0aNJTU2lXbt2\nLFu2jPLlPbc9eb92165djB49mgMHDhAZGckDDzxAly5dSpzd2dC9L0VE/IxlwezZMHEiPP20fZQs\n0GbHxAz9/egs3ftSRCSIHTxoN2EpKVqZXyTUaKZMPNLch3nK3Dx/y/yzz+xh/pYt/X9l/rLyt8xF\n/ImOlImIOCwjAyZMgIUL4d//hq5dna5IRJygmTIREQdt324P8zdqBK+/rmF+KR39/egsb8+U6fSl\niIgDcof5r7oK7rsPFi1SQyYS6tSUiUea+zBPmZvnVOYHD8LNN9tN2dq1cO+9oXN1pX7ORQqnpkxE\nxKDPP7eH+Zs1s4f5W7Z0uiIR8ReaKRMRMSAz0x7mX7AA5s6Fbt2crkiCQY0aNThy5IjTZYSsyMhI\nDh8+XOB5rVMmIuKnfvrJHuaPjrZX5q9Vy+mKJFh4aggkcOn0pXikuQ/zlLl5vs7csuwrKjt3hnvu\ngY8+UkOmn3PzlHng0JEyEREfOHzYHuDfuRO+/BJat3a6IhHxd5opExHxsvh4uPNOuOUW+96VlSo5\nXZGImKSZMhERh2VlweTJ9iD/m29Cjx5OVyQigUQzZeKRZhDMU+bmeTPzXbvs2bHERPuhhswz/Zyb\np8wDh5oyEZGzYFn2/SqvuALuuAOWL4e6dZ2uSkQCkWbKRETK6OhR+xZJ339v30z84oudrkhE/IHu\nfSkiYtC6dfbK/LVqwcaNashE5OypKROPNINgnjI3ryyZnz4NkybBrbfCiy/CSy/Bued6v7ZgpZ9z\n85R54NDVlyIiJfTzzzB4MFSpYg/z16/vdEUiEkx8NlOWnJzMnXfeyW+//YbL5eLee+9lzJgx+faJ\nj4+nb9++NGnSBIBbbrmFCRMmFCxSM2Ui4rCFC2HsWPi//4Nx46CczjOISCH8bp2yChUqMGPGDNq2\nbcuJEye45JJL6NatG61atcq3X5cuXVi6dKmvyhAROSvp6TBqFHz9NXzyCbRr53RFIhKsfPZvvXr1\n6tG2bVsAqlSpQqtWrUhLSyuwn46A+SfNIJinzM0rLvOEBLsJq1gRNm1SQ+YN+jk3T5kHDiMH4JOS\nkkhMTKRjx475nne5XKxfv542bdrQq1cvtm7daqIcEZEiZWfbt0fq0wemToU5c6ByZaerEpFg5/N1\nyk6cOEFsbCwTJkygX79++V5LT08nLCyM8PBwVq5cydixY9mxY0fBIjVTJiKGpKTAkCH2orBvvw0x\nMU5XJCKBpqx9i0+bsqysLG644QZ69uzJuHHjit3/vPPOY9OmTdSoUSN/kS4XQ4cOpXHjxgBERETQ\ntm1bYmNjgT8PzWpb29rW9tlsL1oEw4fHc8stMHt2LGFh/lWftrWtbf/czv19UlISAP/+97/9qymz\nLIuhQ4dSs2ZNZsyY4XGf/fv3U6dOHVwuFwkJCfTv39/9DeUrUkfKjIuPj3f/0IkZyty83Mx//x0e\negg++wwWLIAzJi3Ei/Rzbp4yN8/vrr786quvmD9/PhdffDHt/jcd+9RTT7F3714ARowYwX/+8x9e\neeUVypcvT3h4OO+++66vyhERcVuzfDmfzpxJyv79LAqvy6rkMVx2TW8SE6FaNaerE5FQpXtfikhI\nWbN8OZ+MHcuU3bvdz42p05Rb33yBq3v3drAyEQkWuveliEgJfDpzZr6GDGDmb7tZ9eKLDlUkImJT\nUyYe5R1eFDOUuRnH0jLcv4/P83zYqVPGawlF+jk3T5kHDjVlIhISTp2yb4+0ZVdFj69nV6pkuCIR\nkfw0UyYiQW/rVhg0CJo1g7tuW876x/LPlD3atCk9XtBMmYh4h99dfSki4jTLslfjnzDBXqF/+HBw\nuXpTtQpMfPFFwk6dIrtSJXqMHq2GTEQcpyNl4pHWtTFPmXvXoUPwl79AUhIsXAgtWxbcR5mbp8zN\nU+bm6epLEZH/+eILaNsWmjSBDRs8N2QiIv5GR8pEJGhkZcHkyTB3Lrz5JvTo4XRFIhKKNFMmIiFt\nzx57mD8yEhIToW5dpysSESkdnb4Uj7SujXnKvOzeece+X+WAAbB8eckbMmVunjI3T5kHDh0pE5GA\ndfw4jBoF33wDq1bZc2QiIoFKM2UiEpASEuzTlV27wvTpULmy0xWJiNg0UyYiISEnB559FmbMgJdf\nhltucboiERHv0EyZeKQZBPOUefHS0qBbN1i5EjZuPPuGTJmbp8zNU+aBQ02ZiASEpUuhfXuIjYXV\nqyEmxumKRES8SzNlIuLXTp6Ehx+2r6p85x248kqnKxIRKZpW9BeRoPPjj9ChAxw8aK89poZMRIKZ\nmjLxSDMI5inzP1kWvPKKfaryoYfse1dGRHj/c5S5ecrcPGUeOHT1pYj4lUOHYPhw2LsXvvoKzj/f\n6YpERMzQTJmI+I34eBgyBPr3h6eegooVna5IRKT0tE6ZiAQs3UhcREQzZVIIzSCYF6qZ//wzXH01\nbN5sD/ObbMhCNXMnKXPzlHngUFMmIo559137RuL9+5fuRuIiIsFIM2UiYtyJEzB6NKxfb19Z2b69\n0xWJiHiP1ikTkYCwebPdhLlcsGmTGjIRkVxqysQjzSCYF+yZ5+TA9On2zNgTT9gD/VWqOFtTsGfu\nj5S5eco8cOjqSxHxuf37YdgwOHoUvv4azjvP6YpERPyPZspExKc+/dRuyO66y172okIFpysSEfEt\nrVMmIn4lMxMmTIAFC+wbiV9zjdMViYj4N82UiUeaQTAvmDLftQs6dYLt22HLFv9tyIIp80ChzM1T\n5oFDTZmIeNX8+XDFFTB0KCxZArVqOV2RiEhg0EyZiHhFejo88AB88429KGybNk5XJCLiDK1TJiKO\n2bjRXm+zOnsxAAAgAElEQVSsUiX792rIRERKT02ZeKQZBPMCMfOcHHjuOejVC556CubMgcqVna6q\n5AIx80CnzM1T5oFDV1+KSJn8+qs9N3bihH3KslEjpysSEQlsmikTkVL75BN73bHhw2HSJCivf96J\niLhpnTIR8bnMTHjsMXuQf8ECiI11uiIRkeChmTLxSDMI5vl75rlrj/30EyQmBkdD5u+ZByNlbp4y\nDxxqykSkWLlrj915p9YeExHxFc2UiUih0tNh1Cj7JuLvvaelLkRESkLrlImIV23aBJdcYg/xb9qk\nhkxExNfUlIlHmkEwz18ytyyYMQN69oQnnoA33gistcdKw18yDyXK3DxlHjh09aWIuB04AMOGwaFD\nsGEDNGnidEUiIqFDM2UiAsDq1fYg/x13wD//CRUqOF2RiEhg0jplIlImp0/D5Mkwdy7Mmwfdujld\nkYhIaNJMmXikGQTznMj8l1+gSxf7JuKbN4deQ6afc/OUuXnKPHCoKRMJUR9+CJddBjfdBCtWQN26\nTlckIhLaNFMmEmJOnoQHH4RVq2DhQujQwemKRESCi9YpE5Fi/fij3YQdO2afrlRDJiLiP9SUiUea\nQTDPl5lbFrz2mn2/ygcftG8mXr26zz4uYOjn3Dxlbp4yDxy6+lIkyB07BvfeC9u3w5o10KqV0xWJ\niIgnmikTCWJffw0DB9qr8z/3HJx7rtMViYgEP61TJiJuOTl2EzZtGrz6qn2FpYiI+DfNlIlHmkEw\nz1uZ799vHxlbsgQSEtSQFUU/5+Ypc/OUeeBQUyYSRD77DNq3t9cf+/JLaNTI6YpERKSkNFMmEgSy\nsmDSJPj3v+Gtt6BrV6crEhEJXX63TllycjLXXHMNF1xwARdeeCEzZ870uN+YMWNo3rw5bdq0ITEx\n0VfliAStX36xl7rYvBkSE9WQiYgEKp81ZRUqVGDGjBn8+OOPbNiwgVmzZrFt27Z8+6xYsYJdu3ax\nc+dO5syZw8iRI31VjpSSZhDMK0vmixfbpyr79bNvlVSnjvfrCmb6OTdPmZunzAOHz66+rFevHvXq\n1QOgSpUqtGrVirS0NFrlWSRp6dKlDB06FICOHTty9OhR9u/fT13dhE+kSKdOwd/+Zjdiy5ZBx45O\nVyQiImfLyKB/UlISiYmJdDzjb47U1FRiYmLc29HR0aSkpJgoSYoRGxvrdAkhp6SZ//QTXH65fZXl\n5s1qyM6Gfs7NU+bmKfPA4fOm7MSJE9x666288MILVKlSpcDrZw7CuVwuX5ckErDeegs6d4b77oP3\n34eICKcrEhERb/Hp4rFZWVnccsst3HHHHfTr16/A61FRUSQnJ7u3U1JSiIqK8vhew4YNo3HjxgBE\nRETQtm1bd/efe75c297b3rJlC+PGjfObekJhO/c5T6+fPAkLF8aycSNMnRpPkybgcvlX/YG4fWb2\nTtcTCtvPP/+8/vw2vK0/z838+R0fH09SUhJnw2dLYliWxdChQ6lZsyYzZszwuM+KFSt46aWXWLFi\nBRs2bGDcuHFs2LChYJFaEsO4+Ph49w+dmFFY5lu2wO2320fIZs6EypXN1xas9HNunjI3T5mbV9a+\nxWdN2bp167j66qu5+OKL3ackn3rqKfbu3QvAiBEjABg1ahRxcXFUrlyZuXPn0r59+4JFqimTEGRZ\n8PLLMHkyvPACDBrkdEUiIlISfteUeZOaMgk1R47A8OGQlATvvQfNmztdkYiIlJTfLR4rgS3veXIx\nIzfzDRugXTuIjob//lcNmS/p59w8ZW6eMg8cPh30F5GSy8mBZ5+FadNg9mx7QVgREQkdOn0p4gcO\nHIA774Rjx2DhQt1IXEQkkOn0pUiA+vJL+3Rl27b279WQiYiEJjVl4pFmEHwvOxueeAIGDIDXX4fr\nr4+nQgWnqwot+jk3T5mbp8wDh2bKRBywbx8MHmwve7FpEzRoAPpzU0QktGmmTMSwTz+FoUNhxAiY\nOBHCwpyuSEREvKmsfYuOlIkYcvo0/OMf9v0rFy4ELbAtIiJ5aaZMPNIMgnclJ9tN2ObN9sNTQ6bM\nzVPm5ilz85R54FBTJuJjy5fDZZfBDTfAihVQp47TFYmIiD/STJmIj2Rlwfjx8MEHsGABdOrkdEUi\nImKCZspE/EhSkr3URe3a9unKmjWdrkhERPydTl+KR5pBKLuPPoKOHeG222Dp0pI3ZMrcPGVunjI3\nT5kHDh0pE/GSzEz4+9/tpmzJErj8cqcrEhGRQKKZMhEv2LMHbr8doqJg7lyIjHS6IhERcYrufSni\nkEWL7KNid9wBixerIRMRkbJRUyYeaQaheBkZMHo0/O1v9rIXY8eCy1X291Pm5ilz85S5eco8cGim\nTKQMdu+2T1c2bGhfXRkR4XRFIiIS6DRTJlJKH34II0fChAn2kbKzOTomIiLBR+uUifhYRgY8/DB8\n/PGfq/SLiIh4i2bKxCPNIOS3Z4+9In9Kin260hcNmTI3T5mbp8zNU+aBQ02ZSDFyr64cMsQ+dan5\nMRER8QXNlIkUIjPTPl25dCm89x506OB0RSIiEgg0UybiRUlJ0L8/NGhgn67U2mMiIuJrOn0pHoXy\nDMKSJfZRsQEDzC4GG8qZO0WZm6fMzVPmgUNHykT+JysL/u//4D//sU9Z6t6VIiJikmbKRIDkZHsx\n2Bo14N//hpo1na5IREQCle59KVJGK1bYS1z062cfIVNDJiIiTlBTJh6FwgzC6dPw6KNw773wwQfw\n979DOQf/jwiFzP2NMjdPmZunzAOHZsokJO3bBwMHQoUK9tWVdeo4XZGIiIQ6zZRJyFm9Gu64A+67\nDx57DMLCnK5IRESCidYpEylGTg5MmQIvvwzz50PXrk5XJCIi8ifNlIlHwTaDcPAg9OoFq1bBpk3+\n2ZAFW+aBQJmbp8zNU+aBQ02ZBL3//hfat4e2be1Tlw0aOF2RiIhIQZopk6BlWfD88/DMM/D669Cn\nj9MViYhIKNBMmUgex47B3XfD3r2wYQOcd57TFYmIiBRNpy/Fo0CeQfj2W7j0UqhXD9atC5yGLJAz\nD1TK3Dxlbp4yDxxqyiSovPkmXHcdPP44zJoFFSs6XZGIiEjJaKZMgsLJk/DAA/apyg8/hFatnK5I\nRERCle59KSFr1y64/HLIyICEBDVkIiISmNSUiUeBMoOweDFceaW9Ov/8+VClitMVlV2gZB5MlLl5\nytw8ZR44dPWlBKSsLBg/Hv7zH1i+HC67zOmKREREzo5myiTgpKXB7bfbR8Xmz4eaNZ2uSERE5E+a\nKZOQEB9vL3fRrZt9hEwNmYiIBAs1ZeKRv80gWBb8618wYADMmwf/+AeUC7KfXn/LPBQoc/OUuXnK\nPHBopkz83rFjMGyYfdoyIQEaNnS6IhEREe/TTJn4te++g1tugeuvh2nTtBisiIj4P82USdCZPx+6\ndoVJk+Cll9SQiYhIcFNTJh45OYOQkWGvzv/EE7B6Ndxxh2OlGKW5D/OUuXnK3DxlHjiKnCm76KKL\nin2D2rVrs3r1aq8VJKEtORluuw3q14dvvoHq1Z2uSERExIwiZ8pat27NypUrizwveuONN/Ldd9/5\npLhcmikLDatXw+DBMG4c/P3v4HI5XZGIiEjplbVvKfJI2Zw5c2jUqFGRbzBr1qxSf6hIXrnLXcyY\nAe+8A9de63RFIiIi5hU5U9a5c+di3+Cqq67yWjHiP0zNIBw/bl9d+eGH9nIXodyQae7DPGVunjI3\nT5kHjhIN+i9btox27doRGRlJ1apVqVq1KtWqVfN1bRLktm6171lZpw6sWQMxMU5XJCIi4pwSrVPW\ntGlTFi9ezIUXXkg5B5ZR10xZ8PngA7j/fnj2WbjrLqerERER8R6fzJTlio6O5oILLnCkIZPgcvo0\njB9vN2VxcXDJJU5XJCIi4h9K1GVNnTqVnj178vTTTzNt2jSmTZvG9OnTi/26u+++m7p16xa6tEZ8\nfDzVq1enXbt2tGvXjieffLJ01YvP+GIG4cAB6N4dvv0WNm5UQ3YmzX2Yp8zNU+bmKfPAUaKmbOLE\niVSpUoVTp05x4sQJTpw4QXp6erFfd9dddxEXF1fkPl26dCExMZHExEQmTJhQsqol4HzzDVx6KVx+\nOaxcCbVqOV2RiIiIfynRTNmFF17IDz/8UKYPSEpKok+fPnz//fcFXouPj2fatGksW7as6CI1UxbQ\n3nwTHnkEZs+Gm292uhoRERHf8um9L3v16sUnn3xS6jcvjsvlYv369bRp04ZevXqxdetWr3+GOCcz\nE0aOtIf516xRQyYiIlKUEjVlL7/8Mj179qRSpUpeXRKjffv2JCcn8+233zJ69Gj69et31u8p3nG2\nMwhpaRAbC/v22euPtWrllbKCmuY+zFPm5ilz85R54ChRU3bixAlycnI4deoU6enppKenc/z48bP+\n8KpVqxIeHg5Az549ycrK4vDhw2f9vuKs9evt9cd69YJFi0BL2omIiBSvyCUx9u3bR/369Yt8g5Ls\nU5j9+/dTp04dXC4XCQkJWJZFjRo1PO47bNgwGjduDEBERARt27YlNjYW+PNfAdr27nauku7fpUss\nc+bAI4/E88gjMH68f30/2tb2mduxsbF+VU8obOc+5y/1hMp2Ln+pJ9i2c3+flJTE2Shy0L99+/Zs\n3ry5yDcoap+BAwfy5ZdfcvDgQerWrcvjjz9OVlYWACNGjGDWrFm88sorlC9fnvDwcKZPn87ll19e\nsEgN+vu9jAwYNco+SvbRR9C8udMViYiIOKOsfUuRTVlYWJj79GJhqlWrRmpqaqk/uDTUlJmX91+y\nxUlLs4f4o6Jg3jyoWtWnpQWt0mQu3qHMzVPm5ilz83yyon92dnaZC5LQ8N//wm232VdZPvoouFxO\nVyQiIhKYSrROmdN0pMw/vf663YjNnQu9eztdjYiIiH/w6b0vRfLKzIRx4+CLL2DtWmjRwumKRERE\nAl+5ol7s2bMnP//8s6laxI+cecVOrt9+g+uug5QU2LBBDZk3FZa5+I4yN0+Zm6fMA0eRTdndd9/N\n9ddfz5QpU9xXTUro2rzZXn+sSxf7Csvq1Z2uSEREJHgUO1N24sQJnnjiCT755BOGDBmC63+T3C6X\ni4ceeshMkZopc9y778Lo0fDyy/Zgv4iIiHjms5myChUqUKVKFfdq/uXKFXlwTYJMdjZMmGA3ZZ99\nBm3aOF2RiIhIcCqyKYuLi+Ohhx6iT58+JCYmFrtmmQSP+Ph42rePZdAgOHHCvn9l7dpOVxXctJaQ\necrcPGVunjIPHEU2ZVOmTOGDDz7gggsuMFWP+InUVLj/fnt+bOZMqFDB6YpERESCW5EzZZZluWfI\nnKSZMrM+/xwGDYLJk+1FYUVERKTkfDJT5g8NmZhjWfDSSzBlij1Dds01TlckIiISOjS1LwBkZcF9\n98Hs2fZNxV2ueKdLCjlaS8g8ZW6eMjdPmQcONWXCoUPQvbt9Y/H166FJE6crEhERCT2692WI27oV\nbrwRbr4Znn4awsKcrkhERCSw6d6XUmpxcXDnnfCvf8HQoU5XIyIiEtp0+jIEWRa8+CLcdRcsXuy5\nIdMMgnnK3Dxlbp4yN0+ZBw4dKQsxWVkwdiysWQP//S80bux0RSIiIgKaKQspR45A//72QrDvvgvV\nqjldkYiISPApa9+i05chYvduuPJKaN0ali5VQyYiIuJv1JSFgK++gs6dYfRoeOEFKF+Ck9aaQTBP\nmZunzM1T5uYp88ChmbIgt2ABjBsHb70FPXo4XY2IiIgURjNlQcqy4IknYO5c+PhjuPBCpysSEREJ\nDVqnTNwyMuCee+Cnn2DDBqhXz+mKREREpDiaKQsyR47A9dfDiRPwxRdlb8g0g2CeMjdPmZunzM1T\n5oFDTVkQ2bMHrrgC2reHDz6A8HCnKxIREZGS0kxZkPj6a7jpJnjsMXjgAaerERERCV2aKQthS5bA\nX/4Cb74Jffo4XY2IiIiUhU5fBrhZs2DkSFixwrsNmWYQzFPm5ilz85S5eco8cOhIWYDKyYHx4+Gj\nj2DdOmjSxOmKRERE5GxopiwAZWbCXXdBUpJ9y6SaNZ2uSERERHLp3pch4vhx6NUL/vgDPvtMDZmI\niEiwUFMWQPbtgy5d4Pzz4T//gXPP9d1naQbBPGVunjI3T5mbp8wDh5qyAPHTT3DllXDrrfZwf1iY\n0xWJiIiIN2mmLAB88w3ceCNMmQJ33+10NSIiIlIUrVMWpFatgkGD4I037MZMREREgpNOX/qx99+H\nwYPhww/NN2SaQTBPmZunzM1T5uYp88ChI2V+6pVX7NOVn30GF1/sdDUiIiLia5op8zOWBU8/bZ+u\nXLVKi8KKiIgEGs2UBQHLgr//HeLi7FX669d3uiIRERExRTNlfiI7G+69F9auhS+/dL4h0wyCecrc\nPGVunjI3T5kHDh0p8wNZWTBkCBw4YM+QVanidEUiIiJimmbKHHbqFNx+u32D8Q8+gEqVnK5IRERE\nzobufRmA/vjDXuqiYkVYtEgNmYiISChTU+aQ9HTo2dOeHVuwACpUcLqi/DSDYJ4yN0+Zm6fMzVPm\ngUNNmQOOHYPrr4dWrWDuXCivyT4REZGQp5kyw44ehR494NJL4cUXweVyuiIRERHxJs2UBYAjR6Bb\nN+jYUQ2ZiIiI5KemzJAjR+C66+Cqq+D55/2/IdMMgnnK3Dxlbp4yN0+ZBw41ZQYcPQrdu0NsLEyb\n5v8NmYiIiJinmTIfO37cbsg6dIAXXlBDJiIiEuzK2reoKfOhEyfsof6LLoKXX1ZDJiIiEgo06O9n\nTp6EPn2gZUuYNSvwGjLNIJinzM1T5uYpc/OUeeBQU+YDmZlw663QoAHMng3llLKIiIgUQ6cvvSw7\nGwYNso+Uffih/63ULyIiIr5V1r5Fa8l7kWXBfffBwYOwfLkaMhERESk5nVjzovHj4fvvYcmSwL+5\nuGYQzFPm5ilz85S5eco8cOhImZe88AJ89BGsWwdVqjhdjYiIiAQan86U3X333Sxfvpw6derw/fff\ne9xnzJgxrFy5kvDwcObNm0e7du0KFunnM2XvvgsPP2w3ZI0aOV2NiIiIOMkvl8S46667iIuLK/T1\nFStWsGvXLnbu3MmcOXMYOXKkL8vxidWrYexYWLFCDZmIiIiUnU+bsquuuorIyMhCX1+6dClDhw4F\noGPHjhw9epT9+/f7siSv2roVBgyA996zF4gNJppBME+Zm6fMzVPm5inzwOHooH9qaioxMTHu7ejo\naFJSUhysqOT274feveG55+x7WoqIiIicDcevvjzznKsrAJa+P3kS+vaFIUPgzjudrsY3YtVpGqfM\nzVPm5ilz85R54HD06suoqCiSk5Pd2ykpKURFRXncd9iwYTRu3BiAiIgI2rZt6/5Byz00a2LbsqB3\n73gqV4bHHzf/+drWtra1rW1ta9u/tnN/n5SUxNnw+Yr+SUlJ9OnTx+PVlytWrOCll15ixYoVbNiw\ngXHjxrFhw4aCRfrR1ZdTp9or9a9ZE/hrkRUlPj7e/UMnZihz85S5ecrcPGVunl+u6D9w4EC+/PJL\nDh48SExMDI8//jhZWVkAjBgxgl69erFixQqaNWtG5cqVmTt3ri/LOWtxcfZ6ZAkJwd2QiYiIiHm6\n92UJ7dwJnTrBokXQubOjpYiIiIgf88t1yoLFyZNwyy3w+ONqyERERMQ31JSVwNixcOGF9s3GQ0Xe\n4UUxQ5mbp8zNU+bmKfPAoXtfFmPhQoiPh02bIABW6xAREZEApZmyIuzcCVdeCatWQdu2xj9eRERE\nApBmyrzs9Gl7cdh//EMNmYiIiPiemrJCPPssVKkCDzzgdCXO0AyCecrcPGVunjI3T5kHDs2UefDt\ntzBjBmzeDOXUtoqIiIgBmik7Q1YWXHopPPQQDB1q5CNFREQkiGimzEtefBHq1QveG42LiIiIf1JT\nlkdKCjz1FLz0kpa/0AyCecrcPGVunjI3T5kHDjVleTz0kD3Y37y505WIiIhIqNFM2f98+SXcdRf8\n+COce65PP0pERESCmGbKzoJlwfjx8MQTashERETEGWrKgI8/hvR0GDjQ6Ur8h2YQzFPm5ilz85S5\neco8cIR8U5aTA489BlOmQFiY09WIiIhIqAr5mbLly2HiRN1wXERERLxDM2Vl9MILMG6cGjIRERFx\nVkg3ZVu3wvffw+23O12J/9EMgnnK3Dxlbp4yN0+ZB46QbspmzYIRI6BiRacrERERkVAXsjNl2dnQ\noAFs2ADnnefVtxYREZEQppmyUlq3DqKi1JCJiIiIfwjZpuw//4FbbnG6Cv+lGQTzlLl5ytw8ZW6e\nMg8cIduUrVgBffs6XYWIiIiILSRnyg4dsk9bHj0K5UK2LRURERFf0ExZKWzeDO3aqSETERER/xGS\nbcmmTXDppU5X4d80g2CeMjdPmZunzM1T5oEjJJuyPXvg/POdrkJERETkTyE5U3b77XDTTTBggNfe\nUkRERATQTFmpHD8O1ao5XYWIiIjIn0KyKUtPhypVnK7Cv2kGwTxlbp4yN0+Zm6fMA0dINmXVq8Ox\nY05XISIiIvKnkJwpu+8+uPhiuP9+r72liIiICKCZslKJjobkZKerEBEREflTSDZlTZrATz85XYV/\n0wyCecrcPGVunjI3T5kHjpBsyrp2hdWrISPD6UpEREREbCE5UwbQqRP84x9w/fVefVsREREJcZop\nK6W+fWHJEqerEBEREbGFbFPWvz+8/z4cPOh0Jf5JMwjmKXPzlLl5ytw8ZR44QrYpa9zYvs3SlClO\nVyIiIiISwjNlAPv3Q+vWsHEjnHee199eREREQpBmysqgbl0YPRoeeQT8vzUVERGRYBbSTRnAww/D\njh0wc6bTlfgXzSCYp8zNU+bmKXPzlHngKO90AU6rXBk++gguvxwuuACuu87pikRERCQUhfRMWV7x\n8XD77bB+PTRt6tOPEhERkSCmmbKzFBtrLybbowfs2eN0NSIiIhJq1JTl8cAD8OCD0LkzfPON09U4\nSzMI5ilz85S5ecrcPGUeONSUneH+++HVV6F3b1i2zOlqREREJFRopqwQCQnQrx889pjdqLlcRj9e\nREREAlRZ+xY1ZUXYs8duzBo2tI+eRUcbL0FEREQCjAb9faBJE3u1/8sug3bt4LXXQmeRWc0gmKfM\nzVPm5ilz85R54FBTVoxzzoFJk+CLL+ymrGtX2L3b6apEREQk2Oj0ZSmcPg3PPw/PPAN//at9i6Yq\nVZyuSkRERPyJTl8aUL48/O1vsGEDfPedfXrzmWcgPd3pykRERCTQqSkrg2bNYOFC+y4A331n3wHg\nqafg+HGnK/MezSCYp8zNU+bmKXPzlHngUFN2Flq3hgUL4Msv4ccf7WZtypTgas5ERETEDM2UedH2\n7fDkkxAXBwMGwPDh9lWbIiIiEjr8cqYsLi6Oli1b0rx5c6ZOnVrg9fj4eKpXr067du1o164dTz75\npC/L8bmWLWH+fNi0CWrXhptugvbt4aWX4MgRp6sTERERf+azpiw7O5tRo0YRFxfH1q1bWbhwIdu2\nbSuwX5cuXUhMTCQxMZEJEyb4qhyjGjWyl9HYswemToWvvoLzzoOBA+GzzyAnx+kKi6cZBPOUuXnK\n3Dxlbp4yDxw+a8oSEhJo1qwZjRs3pkKFCgwYMIAlS5YU2C8QTkuWVbly0K2bfVHAnj3QqRM8/LB9\n1ebjj2u9MxEREfmTz5qy1NRUYmJi3NvR0dGkpqbm28flcrF+/XratGlDr1692Lp1q6/KcVyNGjBq\nFCQmwuLFcOCA3aS1bg2PPALr1kF2ttNV/ik2NtbpEkKOMjdPmZunzM1T5oHDZ02ZqwR38G7fvj3J\nycl8++23jB49mn79+vmqHL/Srp09Z5aWBvPm2XcNGDUK6tWDO++EDz7QFZwiIiKhpryv3jgqKork\n5GT3dnJyMtFn3NG7atWq7t/37NmT+++/n8OHD1OjRo0C7zds2DAaN24MQEREBG3btnV3/7nnywNx\nu0MH+OOPeLp2hSZNYvn4Y/jXv+IZOhQ6dYqlTx+oXTue+vXN1rdlyxbGjRvneD6htJ37nL/UEwrb\nZ2bvdD2hsP38888HzZ/fgbKtP8/N/PkdHx9PUlISZ8NnS2KcPn2aFi1a8Pnnn9OgQQM6dOjAwoUL\nadWqlXuf/fv3U6dOHVwuFwkJCfTv39/jNxQoS2J404kTsGoVLFsGy5fbV3P26gVXX22f9oyM9O3n\nx8fHu3/oxAxlbp4yN0+Zm6fMzStr3+LTdcpWrlzJuHHjyM7OZvjw4YwfP57Zs2cDMGLECGbNmsUr\nr7xC+fLlCQ8PZ/r06Vx++eUFiwzBpiyvnBxISIBPPrFnz77+2r7C86qroHNn+9GwodNVioiICPhp\nU+Ytod6UnSkrC779FtautZu0tWvh3HP/bNCuusq+gKBcOacrFRERCT1+uXis+EaFCnDppfDgg/Dh\nh7B/v32q89pr7SNq/fpBrVrQp8+f66RlZJTuM/KeJxczlLl5ytw8ZW6eMg8cPhv0F3NcLjj/fPsx\nfLj93L59djO2di2MGQM//QQXXght2tiPtm3hoosgz7UWIiIi4iCdvgwR6emwZYt92jP3161boUGD\n/I1amzYQE2M3eiIiIlJ6mimTUjt9GnbssBu0vM1aRsafjVpus9a6NVSs6HTFIiIi/k9NmXjN/v3w\n1lvxQKy7Wdu9G5o3z9+otWplH2nTUTXv0GXr5ilz85S5ecrcvLL2LZopkwLq1oXLLoO8/w+fOmWf\n7sw9mvbxx7B9u31atFkzu2E781G3rho2ERGRktKRMjkrx4/Drl2wc2fBR0ZG4Q1brVpq2EREJDjp\n9KX4naNHCzZqO3bYv1qW52ateXP75u0iIiKBSk2ZeJUvZxAsCw4d8nx0bedOex223AatWTP7atCY\nGIiOtn+tUsUnZTlOcx/mKXPzlLl5ytw8zZRJwHC57NOXtWrBFVfkf82y4Lff/mzQdu+GNWsgOdl+\npF8IonkAAAy4SURBVKTYV4HmbdLObNqioyE83JnvTUREpKx0pEwCimXB4cN/Nmm5jVre7dRU+2ha\nYU1b7u+1xIeIiPiCTl+K/I9lwYEDRTduaWkQEZG/aTuzcWvQAM45x+nvRkREAo2aMvGqYJ9ByMmx\n12MrrGlLSYFff4WaNSEqCurUKf5xtg1csGfuj5S5ecrcPGVunmbKREqhXDmoX99+dOjgeZ/Tp+3G\nbN8+e84t97Fvn71eW97nDhywT5kW1rDVrZt/OyLCrkFERCSXjpSJeEFOjr0ESG6Ttn9//qbtzMeJ\nE1C7dtGNW97Huec6/R2KiEhJ6fSlSADJyLCPrhXWtJ3Z1J1zTvGNW5069hpvkZFQqZIW5xURcYqa\nMvEqzSCYV1jmlmXfOaGopi33cfgwHDlif11kZP5HRETJnqtSJXQaOv2cm6fMzVPm5mmmTCRIuVxQ\nvbr9aN68+P0ty75X6ZEj+R9Hj/75+7177XuY5n0ud59Tp/I3ayVt5iIi7BrDwnyfiYhIMNKRMhHJ\nJzMzf7N2ZuNW1HPp6VC1atkaushI+24OJqxZvpxPZ86kfEYGpytWpPuYMVzdu7eZDxeRoKcjZSLi\nFbnza3XqlP5rs7Ph2LGiG7e9ez03eEeP2rNwhTVvuUfiqlSxG7/Cfg0PL/rK1jXLl/PJ2LFM2b3b\n/dxj//u9GjMRcZKOlIlHmkEwL9Qztyz7SFtRR+aOH7evXE1PL/zXkyehcuXCm7ZT667no9RPAYgH\nYv/3+fe2vZ7+/4orsH/VquaO4IWCUP85d4IyN09HykQkoLlcUK2a/WjYsOzvk50Nv/9eeNP2xeYM\nj193OPkUTz9dcP/0dHtOrqijc6X9NTw8dC6mEJGS05EyEQkpE66/nic//bTA8xOvv55/xsUVeN6y\n7CVMijtCV5pfT50q+mhe3qN0JW30yuuf2CJ+Q0fKRERKoPuYMTy2e3e+mbJHmzalx+jRHvd3uexZ\nt0qVoFYt79SQezSvJE3cb7/B7t3F71e+fP4mrXJle9HhMx/h4Z6fL+k+FSvqKJ+Ir+hImXikGQTz\nlLk5a5YvZ9WLL5L866/E1KtHt9GjA3rIP/doXt4m7fff7fm6kyfhjz/+/H1hj+L2yX09K8tuUEvS\n3Hlq8FJS4rn44tgSN4FaCPns6c8W83SkTESkhK7u3Zure/cOmr+s8h7Nq13bt5+VnW2ffi1LY3fw\nICQn528Yi3uPzEz76Jw3jvKVZL9KlXRfWnGOjpSJiIjfyskp2AR648hfYftkZNjLwnijwStqn4oV\n7c/J/VWLLgcX3WZJRETkLOXk2I1ZaU/rlqX5y8y0HxkZ9tG5vE1a3t8X9qu39ilu3woVdPSwtNSU\niVcFy2mdQKLMzVPm5inzgizLPi2c26id+aun50qzz88/x1OrVmyZ3ycz027MTDeLJW0oy5f3v7lD\nzZSJiIgEIJfLbizKl7evmvW2+Hg4mz7YsuwLPLzVJOZelHLo0Nm/T2YmnD7tf81iWelImYiIiASs\nnJziGzhvNJIl3ScjAw4c0OlLEREREceVtW/R6J54FB8f73QJIUeZm6fMzVPm5inzwKGmTERERMQP\n6PSliIiIiBfp9KWIiIhIAFNTJh5pBsE8ZW6eMjdPmZunzAOHmjIRERERP6CZMhEREREv0kyZiIiI\nSABTUyYeaQbBPGVunjI3T5mbp8wDh5oyERERET+gmTIRERERL9JMmYiIiEgAU1MmHmkGwTxlbp4y\nN0+Zm6fMA4eaMhERERE/oJkyERERES/STJmIiIhIAFNTJh5pBsE8ZW6eMjdPmZunzAOHmjIRERER\nP6CZMhEREREv0kyZiIiISABTUyYeaQbBPGVunjI3T5mbp8wDh5oyERERET+gmTIRERERL9JMmYiI\niEgA82lTFhcXR8uWLWnevDlTp071uM+YMWNo3rw5bdq0ITEx0ZflSCloBsE8ZW6eMjdPmZunzAOH\nz5qy7OxsRo0aRVxcHFu3bmXhwoVs27Yt3z4rVqxg165d7Ny5kzlz5jBy5EhflSOltGXLFqdLCDnK\n3Dxlbp4yN0+ZBw6fNWUJCQk0a9aMxo0bU6FCBQYMGMCSJUvy7bN06VKGDh0KQMeOHTl69Cj79+/3\nVUlSCkePHnW6hJCjzM1T5uYpc/OUeeDwWVOWmppKTEyMezs6OprU1NRi90lJSfFVSSIiIiJ+y2dN\nmcvlKtF+Z16dUNKvE99KSkpyuoSQo8zNU+bmKXPzlHngKO+rN46KiiI5Odm9nZycTHR0dJH7pKSk\nEBX1/+3dT0jTfxzH8deWsFqG6aFFPzqYIW2t2leKnUKtQ8xwFXXI6B8ZRIegunXoEMQORYdCqE79\nIcqoU5iXIMtIJKMZhJeCxEkSSB0qlOb2/Z02Gk637Lvt63w+QNj2+ezLhzcvxtvv98vn+9+0Y9XV\n1dGslcCdO3dKvYQFh5oXHzUvPmpefNS8uOrq6ub0vYI1ZZs3b9bHjx81PDysVatW6eHDh3rw4EHG\nnHA4rI6ODu3fv1/9/f1avny5PB7PtGN9+vSpUMsEAACwhYI1ZRUVFero6NCOHTuUSCTU3t4ur9er\nmzdvSpJOnDihlpYWdXd3a+3atVq6dKlu3bpVqOUAAADY2rzY0R8AAKDc2WZH/2PHjsnj8WjDhg0z\nzmGjWWvlqvmLFy9UVVUlwzBkGIYuXrxY5BWWn1gspubmZq1fv15+v1/Xrl3LOo+sWyefmpN1a01O\nTioYDCoQCMjn8+ncuXNZ55Fz6+RTc3JuvUQiIcMw1NramnX8rzNu2kRvb6/57t070+/3Zx1/+vSp\nGQqFTNM0zf7+fjMYDBZzeWUpV817enrM1tbWIq+qvI2NjZnRaNQ0TdP88eOHWV9fbw4NDWXMIevW\nyqfmZN16v379Mk3TNOPxuBkMBs1Xr15ljJNz6+WqOTm33pUrV8wDBw5kretcMm6bM2Vbt25VdXX1\njONsNGu9XDWXpm9Zgn+zcuVKBQIBSVJlZaW8Xq++fPmSMYesWyufmktk3Wput1uS9Pv3byUSCdXU\n1GSMk3Pr5aq5RM6tNDo6qu7ubh0/fjxrXeeScds0Zbmw0WzxORwO9fX1adOmTWppadHQ0FCpl1RW\nhoeHFY1GFQwGMz4n64UzU83JuvWSyaQCgYA8Ho+am5vl8/kyxsm59XLVnJxb68yZM7p8+bKczuyt\n1FwyPm+aMomNZoutoaFBsVhM79+/16lTp7R79+5SL6ls/Pz5U/v27dPVq1dVWVk5bZysW2+2mpN1\n6zmdTg0ODmp0dFS9vb1ZH4pNzq2Vq+bk3DpdXV1asWKFDMOY9ezj32Z83jRl+W40C+ssW7YsfTo8\nFAopHo/r27dvJV7V/BePx7V3714dPHgw648iWbderpqT9cKpqqrSzp079fbt24zPyXnhzFRzcm6d\nvr4+PXnyRLW1tWpra9Pz5891+PDhjDlzyfi8acrC4bDu3r0rSbNuNAvrfP36Nd3lv3nzRqZpZr1H\nAfkzTVPt7e3y+Xw6ffp01jlk3Vr51JysW2t8fDz9EOyJiQk9e/ZMhmFkzCHn1sqn5uTcOpFIRLFY\nTJ8/f1ZnZ6e2bduWznPKXDJesM1j/1ZbW5tevnyp8fFxrV69WhcuXFA8HpfERrOFkqvmjx8/1vXr\n11VRUSG3263Ozs4Sr3j+e/36te7du6eNGzemfzAjkYhGRkYkkfVCyKfmZN1aY2NjOnLkiJLJpJLJ\npA4dOqTt27ezeXgB5VNzcl44qcuS/5pxNo8FAACwgXlz+RIAAKCc0ZQBAADYAE0ZAACADdCUAQAA\n2ABNGQAAgA3QlAEAANgATRkAAIAN0JQBKBuLFi2SYRjpv0uXLk2b09TUpHXr1qmrqyvv405OTioQ\nCMjlcvFYGgAFY5sd/QHgX7ndbkWj0VnnOBwO3b9/Xw0NDXkfd/HixRocHFRtbe2/LhEAZsSZMgAL\nzp8PMmlqatLZs2e1ZcsWeb1eDQwMaM+ePaqvr9f58+dLuEoACw1NGYCyMTExkXH58tGjR1nnpZ5T\nl3rtcrk0MDCgkydPateuXbpx44Y+fPig27dv6/v378VaPoAFjsuXAMrGkiVLcl6+zCYcDkuS/H6/\n/H6/PB6PJGnNmjUaGRlRdXW1pesEgGw4UwZgwXO5XJIkp9OZfp16n0gkSrUsAAsMTRkAzOLP+88A\noJC4fAmgbKTuKUsJhUKKRCJ5f9/hcGTcb5b6DACKgaYMQNmYmprKa96fZ796enrSrxsbG9XY2Jh1\nDAAKjcuXABaUmpoaHT16dE6bx05NTcnp5GcTQGE4TG6YAAAAKDn+5QMAALABmjIAAAAboCkDAACw\nAZoyAAAAG6ApAwAAsIH/AcQcG1ATue6oAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x6805630>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "-------------------------------\n",
      "\n",
      "**Bibliography:**\n",
      "\n",
      "1. Bernetti, R., Titarev, V.A. & Toro, E.F. [Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry](http://dx.doi.org/10.1016/j.jcp.2007.11.033 \"Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry\") (2008) Journal of Computational Physics (227), pp. 3212\u20133243.\n",
      "1. Castro D\u00edaz, M.J., L\u00f3pez-Garc\u00eda, J.A. & Par\u00e9s C. [High order exactly well-balanced numerical methods for shallow water systems](http://dx.doi.org/10.1016/j.jcp.2013.03.033 \"High order exactly well-balanced numerical methods for shallow water systems\") (2013) Journal of Computational Physics, 246, pp. 242\u2013264.\n",
      "1. Chow, V.T., [Open-Channel Hydraulics](http://blackburnpress.stores.yahoo.net/ophy.html \"Open-Channel Hydraulics\"), 2009, Blackburn Press.\n",
      "1. LeFloch P.G. & Duc Thanh M., [A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime](http://dx.doi.org/10.1016/j.jcp.2011.06.017 \"A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime\") (2011) Journal of Computational Physics, 230, pp. 7631\u20137660.\n",
      "1. Valiani, A. & Caleffi, V., [Depth-energy and depth-force relationships in open channel flows: Analytical findings](http://dx.doi.org/10.1016/j.advwatres.2007.09.007 \"Depth-energy and depth-force relationships in open channel flows: Analytical findings\") (2008) Advances in Water Resources, 31, pp. 447-454."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}