Example_Direct_Step.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Example_Direct_Step.ipynb",
      "provenance": [],
      "authorship_tag": "ABX9TyPd8yT4jhkaTsbbyMdAde4v",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/mqquiroz/ICV411/blob/main/Example_Direct_Step.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y0VGG8mIKjc_"
      },
      "source": [
        "GVD example Direct Step (Moglen 6.1)\n",
        "\n",
        "[Video Youtube](https://www.youtube.com/watch?v=L5oBguZC_2c)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "o1jk8Tx8BSBO"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.optimize import fsolve,newton"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "15iLHlL-BarQ"
      },
      "source": [
        "#Manning dor a trapezoidal or rectangular channel\n",
        "def manning(y,Q,b,S,n,Co,m):\n",
        "  A = (b+m*y)*y\n",
        "  R = A/(b+2*y*(1+m**2)**(1/2))\n",
        "  return Q-(Co/n)*A*R**(2/3)*S**(1/2)\n",
        "\n",
        "#Critical flow depth\n",
        "def critT(y):\n",
        "  A = (b+m*y)*y\n",
        "  B = b+2*m*y\n",
        "  return Q/np.sqrt(g*A**3/B)-1"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UXkX6AOhCodM",
        "outputId": "766573a9-26fa-456e-cbfc-96ebc2c6da48"
      },
      "source": [
        "#Conditions for the M2 profile using channel properties from example 6.1\n",
        "\n",
        "g = 9.81\n",
        "Q = 30\n",
        "b = 10\n",
        "n = 0.03\n",
        "m = 0\n",
        "So = 0.0050\n",
        "x0 = 0 #Starting point\n",
        "\n",
        "if g>10:\n",
        "  Co = 1.486\n",
        "else:\n",
        "  Co = 1\n",
        "\n",
        "#Specify depth (crit and normal)\n",
        "yc = newton(critT,0.2)\n",
        "yo = fsolve(manning,1,args=(Q,b,So,n,Co,m))\n",
        "print('Critical depth',np.around(yc,decimals=3))\n",
        "print('Normal depth',np.around(yo,decimals=3))"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Critical depth 0.972\n",
            "Normal depth [1.265]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vQMa8A7lEKG6"
      },
      "source": [
        "#Specify depths\n",
        "yd = yc #downstream\n",
        "yu = np.around(yo,decimals=6) #upstream\n",
        "yu = 1.2619 #needed to match Moglen"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rmJS5UIIEu_-"
      },
      "source": [
        "#Calculate profile [small number if points]\n",
        "\n",
        "N = 8 #number of calculation points\n",
        "y = np.array([yd,1,1.05,1.1,1.15,1.2,1.25,yu])\n",
        "\n",
        "#N = 100\n",
        "#y =np.linspace(yd,yu,N)\n",
        "\n",
        "dx = np.zeros(N-1)\n",
        "x = np.zeros(N)\n",
        "\n",
        "A = (b+m*y)*y\n",
        "R = A/(b+2*y*(1+m**2)**(1/2))\n",
        "Sf = (Q*n/(Co*A*R**(2/3)))**2\n",
        "E = y+Q**2/(2*g*A**2)\n",
        "\n",
        "dx = (E[1:]-E[:-1])/(So-0.5*(Sf[1:]+Sf[:-1]))\n",
        "x[0] = x0\n",
        "x[1:] = dx\n",
        "x = np.cumsum(x)"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SyCNrLOkHeK1",
        "outputId": "6ac1f575-2610-4e2b-b2c8-0e2194a6cd67"
      },
      "source": [
        "print(np.around(x,decimals=2))"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "[   0.     -0.2    -1.8    -5.78  -13.87  -30.72  -77.92 -136.11]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "id": "xDAigmSOHkEr",
        "outputId": "e6dd52ba-2519-42a8-aca2-9febfe29b3dc"
      },
      "source": [
        "plt.plot(x,y,linestyle='-',marker='o',color='b')\n",
        "plt.ylabel('depth, y [m]')\n",
        "plt.xlabel('distance, x [m]')\n",
        "plt.show()\n",
        "print('distance to the upstream bc =', np.around(x[N-1],decimals=2))"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5xcdX3/8dc79wQIARJuCcnKvdwvK1AQCkTut0CkrUYRtU2xtd7AH9C0grZB/IlCabUklcsvdo2oBEUEBFEKiFwSIOHShIZLEsIlARIkxIRgPr8/vmfdyTKzO7vJmTOz834+HucxM98zM+fDIbvvPed7vt+jiMDMzKyzfkUXYGZm9ckBYWZmZTkgzMysLAeEmZmV5YAwM7OyBhRdwKY0cuTIaGlpKboMM7OGMWfOnNciYlS5dX0qIFpaWpg9e3bRZZiZNQxJiyqt8ykmMzMrywFhZmZlOSDMzKwsB4SZmZXlgDAzs7KaPiDa2qClBfr1S49tbUVXZGZWH/rUZa491dYGkyfD6tXp9aJF6TXApEnF1WVmVg+aOiCmTOkIh3arV8MXvwi77gojRnQsgwcXU6OZWVGaOiAWLy7fvmwZHHbYhm1DhqSg2HLLDYOj2mXIkPz/e8zMNqWmDoixY9Nppc622w6uvx5Wrqy8rFgBzz/f8Xzduq63NXjwhoHR06AZMgSkfPaDmVk5TR0QU6du2AcBMGwYfPObcNJJ1X9PBKxZ0xEeb77Zdbi0L4sWdQTMO+90vY1BgyqHRzVhM3SoA8bMeqapA6K9I3rKlHS6aezYFBo97aCW0i/goUNhhx16V0tpwFQbMosXdzxfs6br7x84sLojlUphM2yYA8as2agv3ZO6tbU1mnWyvjVr3hsq1R7JrFwJv/99198/YEDPTol1DprNNnPAmNUjSXMiorXcutyOICRdB5wKLIuIfcqsnwRcCAh4C/h0RMzN1r2Qtf0BeLdS8dZhyJC0bLdd7z6/dm3lQKnU/tJLHc87Xw3WWf/+vevcr1XAtLVt/JGkWV+T5ymmG4B/B2ZUWP888GcRsULSScB04NCS9cdExGs51mclBg+GbbdNS2+8807XRyzl1s2f3/H87be7/v7+/Xt/BdmIEbD55pUDxuNhzMrLLSAi4l5JLV2sf6Dk5YPAmLxqsfwNGgSjRqWlN9atq+6UWOl7nnmm4/mqVV1/f79+lftYfvzj8uNhpkxxQFhzq5dO6k8Bt5e8DuBOSQFMi4jplT4oaTIwGWDs2LG5Fmn5GTgQRo5MS2+8+271fS7t71u4MD2+9Vb571y8OB0ZDRrU+/8us0ZWeEBIOoYUEB8oaf5ARCyVtC1wl6T5EXFvuc9n4TEdUid17gVbXRowALbZJi091dJSfjxMRDrlduqpMHEinHBCuprLrFkUOlmfpP2A7wJnRMTr7e0RsTR7XAbcDBxSTIXWDKZOfe8v/mHD4IILUjDccQecdVY6ffahD8HMmfC73xVTq1ktFRYQksYCs4CPRcQzJe2bSdqi/TlwPPBkMVVaM5g0CaZPh3HjUkf2uHHp9Te+AddeC6+8AnffDeeeCw88AB/5SAqLU09NI+5ff73bTZg1pNzGQUiaCRwNjAReBS4BBgJExDWSvgtMBNoP7t+NiFZJO5OOGiCdAvt+REytZpvNPA7CamP9enjwQbjpprQsWpSusDr66HSUceaZvR8saVaErsZBeKCcWS9FwGOPdYTFggXpCOTww1NYnHVW6t8wq2cOCLMaePppmDUrhcXjj6e2gw5K/RgTJ8IeexRbn1k5DgizGnv2Wbj55hQWDz6Y2vbaKwXFWWfB/vt76hGrDw4IswK9+CL85CcpLO69N/Vj7LxzR1gcckgayGdWBAeEWZ1Yvhx++tMUFnffnUaQjx6dOrcnToQjj0yd3ma14oAwq0MrV8Ktt6Z+i9tvTzPyjhoFZ5yRwuLYYz2K2/LngDCrc2+/nUJi1qwUGm+9leaNOu20FBbHH+9R3JYPB4RZA1mzJp1+uummdDrqjTdSOJx8cuqzOOUUGD686Cqtr3BAmDWodetSx/ZNN6Wrol55JZ12Ov74FBann967+afM2jkgzPqA9evht79NYTFrVsco7mOO6RjFvf32RVdpjcYBYdbHRMCjj3YMzCsdxd1++ey4cUVXaY3AAWHWh0VsOIp77tzUfvDBHWHhUdxWiQPCrIksXNgxivuhh1Lb3nunoJg4Efbbz6O4rYMDwqxJLVnSMYr7vvtSP8Yuu3SExfvf71Hczc4BYWYsW7bhKO5334UxYzpGcX/gAx7F3YwcEGa2gZUr4Wc/S/0Wd9zRMYp7woR0dOFR3M3DAWFmFa1ateEo7lWrYMSINIr7rLPSvbiHDi26SsuLA8LMqrJmDfzylx2juFes6BjFPXFiGsW9xRZFV2mbkgPCzHps3Tr47//uGMX96qsweDAcd1wKi9NPh623LrpK21gOCDPbKH/4w4ajuBcv7hjFPXFi6rvwKO7G5IAws00mAubM6RiY98wzaVzFEUd03Ivbo7gbhwPCzHIRAU891REW8+al9tbWjrEWu+9ebI3WNQeEmdXEwoUdYfHww6lt771TUEycCPvu61Hc9cYBYWY1t2RJx5Qf992XjjZ23XXDUdwOi+I5IMysUK++mi6bnTVrw1Hc7X0WHsVdHAeEmdWNFSs2HMW9dm3HKO6JE9OVUR7FXTtdBURu03RJuk7SMklPVlg/SdI8SU9IekDS/iXrTpS0QNJCSRflVaOZ1d5WW8E556RJBF97DW68MU3tMXMmnHgibLddWv/Tn8Lvf9/xubY2aGlJkwu2tKTXlq/cjiAkHQWsAmZExD5l1h8O/E9ErJB0EnBpRBwqqT/wDHAc8CLwCPDhiHi6u236CMKsca1ZA3fdlfosbrklHWlstlkaxb3ddnDttRsGxrBhMH06TJpUXM19QWGnmCS1ALeWC4hO79sKeDIiRkv6U1JYnJCtuxggIr7W3fYcEGZ9w7p1cM89HaO4ly0r/75x4+CFF2pZWd9TyCmmHvoUcHv2fDSwpGTdi1lbWZImS5otafby5ctzLNHMamXgwDSlxzXXwEsvVb7aafHi2tbVbAoPCEnHkALiwt58PiKmR0RrRLSOGjVq0xZnZoXr3x/Gji2/bsyY2tbSbAoNCEn7Ad8FzoiI17PmpcBOJW8bk7WZWZOaOjX1OXQmwSuv1L6eZlFYQEgaC8wCPhYRz5SsegTYTdL7JA0C/hK4pYgazaw+TJqUOqTHjUuhMG4cfOlL6Sqoww9PI7ht08vzKqaZwNHASOBV4BJgIEBEXCPpu8BEYFH2kXfbO0oknQxcBfQHrouIqdVs053UZs3loYfSPSr69YPbbktzQFnPeKCcmfVZCxaku9699lq64um444quqLE0wlVMZma9ssce8MADsMsu6Whi5syiK+o7HBBm1vB23DHd/e5P/xQ+8hG46qqiK+obHBBm1ieMGAG/+EWa/O8LX4CLLkozyFrvOSDMrM8YMgR++EM47zz4+tfh3HPTqGzrnQFFF2Bmtin17w/f+Q7ssANccgksXw4/+lGa18l6xkcQZtbnSPDlL8O0aem00/jx6Son6xkHhJn1WZMnpwn/Hn883ZRo0aLuP2MdHBBm1qdNmAB33pmm5Dj8cHjiiaIrahwOCDPr8446Kt0XG+DIIzueW9ccEGbWFPbdNw2o2377NNr6Jz8puqL654Aws6Yxbhzcfz8ccEC6//W0aUVXVN8cEGbWVEaOhLvvTve/Pu88+OpXPaCuEgeEmTWdzTZLp5g+/vE0VuJv/xb+8Ieiq6o/HihnZk1p4EC4/vrUJ/H1r6f7Xre1pdHYlvgIwsyalgSXXw5XXgmzZqXTTitXFl1V/XBAmFnT+/zn4fvfT1c5HXUUvPRS0RXVBweEmRnw4Q/Dz38Ozz+fBtQtWFB0RcVzQJiZZY47Du65B1avhiOOgIcfLrqiYjkgzMxKHHww/OY3MHw4HHMM3HFH0RUVxwFhZtbJbrul/ojdd4fTToPvfa/oiorhgDAzK2P77dNtTI86Cs45B664ouiKas8BYWZWwfDhcNtt8Od/Dl/6Epx/PqxfX3RVteOBcmZmXRg8GGbOhG23hW99C159Fa67DgYNKrqy/FUMCEm/6+azAl6OiN03bUlmZvWlXz+4+up0G9MpU9JtTG+6CTbfvOjK8tXVKaZnI2J4F8sWwNu1KtTMrEgS/MM/wLXXpsn+jj02BUVf1lVATKzi8xXfI+k6ScskPVlh/Z6SfitpraQLOq17QdITkh6XNLuKOszMauKTn4Sbb053pjviCLjqKmhpSUcZLS1pPqe+QlHlPLeShlNySioi3ujm/UcBq4AZEbFPmfXbAuOACcCKiLiiZN0LQGtE9Og2462trTF7tvPEzPL3wANpYN3q1Ru2DxsG06fDpEnF1NVTkuZERGu5dd1exSTpbyS9AswD5mRLt7+FI+JeoGKIRMSyiHgEWNfdd5mZ1ZvDD4ctt3xv++rVqZ+iL6jmKqYLgH16+tf8RgrgTkkBTIuI6ZXeKGkyMBlg7NixNSrPzAxeeaV8++LFta0jL9WMg3gWWN3tuzatD0TEQcBJwN9lp6vKiojpEdEaEa2jRo2qXYVm1vQq/U3aV/5WrSYgLgYekDRN0tXtS55FRcTS7HEZcDNwSJ7bMzPrjalTU59DqaFDU3tfUM0ppmnAr4AngNzHEEraDOgXEW9lz48Hvpr3ds3Meqq9I3rKlHRaKQIOO6xxOqi70+1VTJIei4gDe/zF0kzgaGAk8CpwCTAQICKukbQ9qbN7OCl4VgF7Ze+/OfuaAcD3I6KqPPZVTGZWpM9+Fr7zHXj66TTRXyPo6iqmagLiMuAF4GfA2vb27i5zLYIDwsyK9OqrsMsucMopcOONRVdTnY26zBX4MFk/BD24zNXMrNlst12a0O+HP4Q5c4quZuN1GxAR8b4yy861KM7MrNGcfz5ss02alqPRVQwISQd19+Fq3mNm1kyGD0+d1nfeCb/6VdHVbJyujiCul7SVpK0rLcC1tSrUzKxRfPrTsNNOcNFF6cqmRtXVZa5bkvob1MV7+vhchmZmPTdkCHzlKx0T+511VtEV9U7Vk/U1Al/FZGb14t13Yb/90hHEE0/AgDq9PdvGXsVkZmY9NGBAGlE9fz7MmFF0Nb3jgDAzy8mECXDIIXDJJbBmTdHV9JwDwswsJxJcfjm8+GIaYd1oqrkfxCxJp0hymJiZ9dAxx8AJJ8Bll8GbbxZdTc9U80v/O8BHgP+VdLmkPXKuycysT7nsMnj9dfjmN4uupGeqGUn9y4iYBBxEmpPpl5IekPQJSQPzLtDMrNEddBD8xV/At76V5mtqFFWdNpK0DXAu8FfAY8C/kgLjrtwqMzPrQ/75n1NH9b/8S9GVVK+aPoibgfuAYcBpEXF6RNwYEX8PbJ53gWZmfcFuu8Ff/RVMmwbPPVd0NdWp5gji6ojYKyK+FhEvl66oNLjCzMze68tfhv7902WvjaCaPohf16IQM7O+bscd4XOfg7Y2mDev6Gq650tXzcxq6MILYcst04yv9c4BYWZWQ1ttlULi1lvh/vuLrqZrPQ4ISTtIGpxHMWZmzeCzn4Uddqj/6cB7cwTxPWC+pCs2dTFmZs1g2LDUUf2b38BttxVdTWU9DoiI+CCwM3D9pi/HzKw5fPKTsOuucPHFsH590dWUV+1Auf6SdpQ0VtJYYKeIeCrn2szM+qyBA9OguSeegO9/v+hqyqtmoNzfA6+SRk3/PFtuzbkuM7M+7+yz4cAD4Z/+Cd55p+hq3quaI4jPAXtExN4RsW+27Jd3YWZmfV2/fvC1r8ELL8D06UVX817VBMQSoMEmqTUzawzHHw9HH53malq1quhqNlQxICR9UdIXgeeAeyRd3N6WtXdJ0nWSlkl6ssL6PSX9VtJaSRd0WneipAWSFkq6qKf/UWZmjUJKRxHLlsFVVxVdzYa6OoLYIlsWk/ofBpW0VTNJ3w3AiV2sfwP4LLDB5bKS+gPfBk4C9gI+LGmvKrZnZtaQDjsMDj44zdXUrx+0tKTpOIo2oNKKiPgKgKSzI+JHpesknd3dF0fEvZJauli/DFgm6ZROqw4BFkbEc9m2fgCcATzd3TbNzBpRWxs89VTHoLlFi2Dy5PR80qTi6qqmD+LiKts2ldGkfo92L2ZtZmZ90pQp6V4RpVavLn6+popHEJJOAk4GRku6umTVcODdvAurlqTJwGSAsWPHFlyNmVnPLV7cs/Za6eoI4iVgNrAGmFOy3AKckGNNS4GdSl6PydrKiojpEdEaEa2jRo3KsSwzs3xU+tu26L95KwZERMyNiP8H7ArMJN1q9FHg1ohYkWNNjwC7SXqfpEHAX5JCycysT5o6Nc3PVGrYsNRepIqnmEocB0wDngUEvE/S30TE7V19SNJM4GhgpKQXgUuAgQARcY2k7UlHKMOB9ZI+D+wVEb+T9BngF0B/4DpP62FmfVl7R/SUKamDGuDKK4vtoAZQdDPXrKT5wKkRsTB7vQvw84jYswb19Uhra2vMnj276DLMzHrtvvvgqKNg1iw488z8tydpTqXbR1dzFdNb7eGQeQ54a5NUZmZmGzj00HR66e67i66kulNMsyXdBvwQCOBs4BFJZwFExKwc6zMzayqDBqUjiHoIiGqOIIaQZnP9M1KfwnJgKHAacGpulZmZNanx42H+fFha8frN2uj2CCIiPlGLQszMLBk/Pj3++tfw0Y8WV0c194PYXdLd7ZPuSdpP0j/mX5qZWXPaf3/YeuviTzNVc4rpP0lTa6wDiIh5pLEJZmaWg3794JhjUkB0c6FpvnVU8Z5hEfFwp7a6mWrDzKwvGj8eliyBhQu7f29eqgmI17KxDwEg6UPAy7lWZWbW5Nr7IX71q+JqqCYg/o40knpPSUuBzwOfzrUqM7Mmt9tuMGZMsf0Q1VzF9BzwQUmbAf0iwoPkzMxyJsGxx8LPfw7r16d+iVrrarrvsrcVlQRARHwrp5rMzIx0mmnGDJg3Dw44oPbbr+aWo62kU0qjs+U84KD8SzMza27t/RBFnWbqarrvr2S3HR0DHBQR50fE+cDBgO/MY2aWs9GjYY89iuuoruas1nbAOyWv38nazMwsZ8ceC/feC+vW1X7b1QTEDOBhSZdKuhR4CLghz6LMzCwZPx5WrYKHO49Gq4FuAyIipgKfAFZkyyci4mt5F2ZmZmlEtVRMP0RVF05FxKMR8a/Z8ljeRZmZWbL11nDggcX0QxRwZa2ZmfXE+PHw29/C6tW13a4Dwsyszh17LLzzDtx/f22364AwM6tzRx4JAwfWvh/CAWFmVuc22wwOO8wBYWZmZYwfD48+CitW1G6bDggzswYwfny6edA999Rumw4IM7MGcMghMGxYbU8zOSDMzBrAoEFw1FF9JCAkXSdpmaQnK6yXpKslLZQ0T9JBJev+IOnxbLklrxrNzBrJ+PEwfz689FJttpfnEcQNwIldrD8J2C1bJgP/UbLu9xFxQLacnl+JZmaNY+3a9Dh6NLS0QFtbvtvLLSAi4l7gjS7ecgYwI5IHgRGSdsirHjOzRtbWBpdd1vF60SKYPDnfkCiyD2I0sKTk9YtZG8AQSbMlPShpQu1LMzOrL1OmvHeqjdWrU3teur0ndUHGRcRSSTsDv5L0REQ8W+6NkiaTTlExdqzvY2RmfdPixT1r3xSKPIJYCuxU8npM1kZEtD8+B9wDHFjpSyJiekS0RkTrqFGj8qvWzKxAlf7+zfPv4iID4hbgnOxqpsOANyPiZUlbSRoMIGkkcATwdIF1mpkVburUNA6i1LBhqT0vuZ1ikjQTOBoYKelF4BJgIEBEXAPcBpwMLARWk25KBPAnwDRJ60kBdnlEOCDMrKlNmpQeL7gAXnkFRo2CK6/saM+DIiK/b6+x1tbWmD17dtFlmJnlZsmSdFpp+nT467/e+O+TNCciWsut80hqM7MGsv326RakS5fmvy0HhJlZAxk4ELbd1gFhZmZljB7tgDAzszIcEGZmVpYDwszMyho9Gl5/HdasyXc7DggzswYzOpu1Lu9pvx0QZmYNxgFhZmZl7bhjesy7H8IBYWbWYNqPIBwQZma2gREjYOhQB4SZmXUi1eZSVweEmVkDckCYmVlZDggzMytrxx3TZa553rHBAWFm1oBGj4a1a+GNN/LbhgPCzKwB1eJSVweEmVkDckCYmVlZDggzMytrhx3SowPCzMw2MGhQ/rcedUCYmTWoHXd0QJiZWRmjR+c75bcDwsysQeU9mtoBYWbWoEaPhuXL04C5POQaEJKuk7RM0pMV1kvS1ZIWSpon6aCSdR+X9L/Z8vE86zQza0Ttl7q+/HI+35/3EcQNwIldrD8J2C1bJgP/ASBpa+AS4FDgEOASSVvlWqmZWYNZsCA97rwztLRAW9um/f5cAyIi7gW6minkDGBGJA8CIyTtAJwA3BURb0TECuAuug4aM7Om0tYG//Zv6XkELFoEkydv2pAoug9iNLCk5PWLWVuldjMzA6ZMgTVrNmxbvTq1bypFB8RGkzRZ0mxJs5cvX150OWZmNbF4cc/ae6PogFgK7FTyekzWVqn9PSJiekS0RkTrqFGjcivUzKyejB3bs/beKDogbgHOya5mOgx4MyJeBn4BHC9pq6xz+viszczMgKlTYdiwDduGDUvtm8qATfdV7yVpJnA0MFLSi6QrkwYCRMQ1wG3AycBCYDXwiWzdG5L+GXgk+6qvRkSOt8UwM2sskyalx499LHVSjxuXwqG9fVNQ5Hm/uhprbW2N2bNnF12GmVnNjBgB554LV13Vu89LmhMRreXWFX2KyczMNsKgQQ06ktrMzPI1aBC8804+3+2AMDNrYA4IMzMrywFhZmZlOSDMzKwsB4SZmZU1eLADwszMyvARhJmZleWAMDOzshwQZmZWlgPCzMzeo60N7rgDnnyyAW85amZm+WhrS7cYXb06ve6Ltxw1M7NemDKlIxza+ZajZmbWFLccNTOzXmiGW46amVkv1OKWow4IM7MGNGkSTJ+ebjUqpcfp0zftLUdzvSe1mZnlZ9KkTRsInfkIwszMynJAmJlZWQ4IMzMrywFhZmZlOSDMzKwsRUTRNWwykpYDi8qsGgm8VuNyNpZrzl+j1QuuuRYarV7YuJrHRcSociv6VEBUIml2RLQWXUdPuOb8NVq94JprodHqhfxq9ikmMzMrywFhZmZlNUtATC+6gF5wzflrtHrBNddCo9ULOdXcFH0QZmbWc81yBGFmZj3kgDAzs7L6XEBIOlvSU5LWS3rPZV+SxkpaJemCkrYTJS2QtFDSRfVQr6TjJM2R9ET2eGzJuoOz9oWSrpakeqg5W3dxVtcCSSeUtBe2jzuTdICkByU9Lmm2pEOydmX7c6GkeZIOKrLOziT9vaT52b7/vyXtZfd5PZB0vqSQNDJ7Xbf7WNI3sv07T9LNkkaUrKvnfZzfz1ZE9KkF+BNgD+AeoLXM+h8DPwIuyF73B54FdgYGAXOBvYquFzgQ2DF7vg+wtGTdw8BhgIDbgZPqYR8De2X7bzDwvmy/9i96H5ep/872fQacDNxT8vz2bL8eBjxU5L/lTjUfA/wSGJy93rarfV50vVltOwG/IA1eHdkA+/h4YED2/OvA1xtgH+f6s9XnjiAi4n8iYkG5dZImAM8DT5U0HwIsjIjnIuId4AfAGflXmlSqNyIei4iXspdPAUMlDZa0AzA8Ih6M9C9kBjChVvVmtVXax2cAP4iItRHxPLCQtH8L3cdlBDA8e74l0L6fzwBmRPIgMCLb3/Xg08DlEbEWICKWZe2V9nk9uBL4P6T93a5u93FE3BkR72YvHwTGZM/reR/n+rPV5wKiEkmbAxcCX+m0ajSwpOT1i1lbPZkIPJr9chhNqrFdPdVbaV/W2z7+PPANSUuAK4CLs/Z6q7PU7sCRkh6S9N+S3p+112XNks4gHfXO7bSqLust45OkIx2o75pzra0h7ygn6ZfA9mVWTYmIn1b42KXAlRGxqsan7Htbb/tn9yYd7h6fR21dbLfXNdeDruoHxgNfiIibJP05cC3wwVrWV043NQ8Atiadlnk/8ENJO9ewvPfopt5/oMb/ZqtRzb9rSVOAd4G2WtZWjxoyICKiNz/MhwIfyjr3RgDrJa0B5pDOlbYbAyzd+Co79LJeJI0BbgbOiYhns+aldBz6Qg71Qq9rXkrlfZnrPu6sq/olzQA+l738EfDd7HlX9eeum5o/DczKTis+LGk9aYK2wmquVK+kfUnn6udmf4yNAR7NLgao230MIOlc4FRgfLavoeCau5FvbUV3suS1UKGTOlt3KR2d1AOA50j/oNs7efYuul5SiM0Fzirz3s6d1CfXwz4G9mbDzrznSJ1odbGPS+r8H+Do7Pl4YE72/BQ27EB9uKgay9R8HvDV7PnupNMKqrTPi663U+0v0NFJXc/7+ETgaWBUp/a63cd5/2wV/h+Yww47k3Qebi3wKvCLMu/5Y0Bkr08GniFdDTClHuoF/hF4G3i8ZGm/cqUVeDKr99/JRsQXXXO2bkpW1wJKrq4qch+Xqf8DpCPHucBDwMFZu4BvZzU+QYU/MAqqeRDwX9n/90eBY7vb5/WydAqIet7HC7Pgbf95u6YR9nGeP1ueasPMzMpqmquYzMysZxwQZmZWlgPCzMzKckCYmVlZDggzMyvLAWFmZmU5IKxPk3Rp+9Tukr4qqavRyhMk7VW76npP0rmSlkv6bvfv3uBz35D0Sul092aVNORUG2a9ERFf7uYtE4BbSaNpG8GNEfGZnnwgIr4k6e28CrK+xUcQ1udImiLpGUn3k+5b0d5+g6QPZc8vl/R0dnOYKyQdDpxOmuX1cUm7SPprSY9ImivpJknDSr7nakkPSHqu/TuzdRcq3cxprqTLs7ZdJN2hdOOn+yTt2U39X5B0XfZ8X0lPtm+7i8+cK+knku6S9IKkz0j6oqTHlG6OtHUvd6c1MR9BWJ8i6WDgL4EDSP++HyVNq1H6nm1I04XsGREhaURErJR0C3BrRPw4e9/KiPjP7Pm/AJ8C/i37mh1IU3bsCdwC/FjSSaS5+A+NiNUlv5SnA+dFxP9KOhT4DvDHOwSW8bb7IcgAAAHqSURBVK/APZLOJE3x8DcRsbqK//x9SDeaGkKaNuLCiDhQ0pXAOcBVVXyH2R85IKyvORK4uf0XavZLv7M3gTXAtZJuJZ1WKmefLBhGAJuT7o7W7icRsR54WtJ2WdsHgevbtx0Rb2T3ITkc+FHJNPODu/oPiIj12ayi84BpEfGbrt5f4tcR8RbwlqQ3gZ9l7U8A+1X5HWZ/5ICwphMR72ZTT48HPgR8hvJ/0d8ATIiIudkv7KNL1q0ted7VDUb6ASsj4oAelrkbsArYsQefKa1pfcnr9fhn3XrBfRDW19wLTJA0VNIWwGmd35D9Vb9lRNwGfAHYP1v1FrBFyVu3AF6WNBCYVMW27wI+UdJXsXVE/A54XtLZWZsk7Z89P1PS18rUtyVwNXAUsE1pH4dZLTkgrE+JiEeBG0lTed8OPFLmbVsAt0qaB9wPfDFr/wHwpaxjdxfgn0jTgf8GmF/Ftu8g9UfMlvQ40H4p6STgU5Lmku4v3n7P4F2A35X5qiuBb0fEM6R+j8slbdvd9s02NU/3bVYQSf9FuvXp8l589lzSvRR6dJlr9tlLgVURcUVPP2vNxUcQZgWJiI/2JhwyvwdO6s1AOeCjpJtRmXXJRxBmZlaWjyDMzKwsB4SZmZXlgDAzs7IcEGZmVtb/Bwf2p0nvPQiJAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        },
        {
          "output_type": "stream",
          "text": [
            "distance to the upstream bc = -136.11\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "id": "RTbNXNhGIHQh",
        "outputId": "f61b8269-d83f-4864-835e-27a9c77a2f3f"
      },
      "source": [
        "eta = -x*So\n",
        "WSE = eta+y\n",
        "\n",
        "plt.plot(x,eta,'k-',label='bed')\n",
        "plt.plot(x,WSE, linestyle='-',color='b',label='WSE')\n",
        "plt.legend()\n",
        "plt.ylabel('Elevation, [m]')\n",
        "plt.xlabel('distance, x [m]')\n",
        "plt.xlim(x[N-1],x[0])\n",
        "plt.show()\n",
        "print('distance to theupstream bc =',np.around(x[N-1],decimals=2))"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZhU1Z3/8feXXZAdRPbFoMgmS7MIDNHghmNUIpmo0agYiUk0A2acxDi/0ZhMYkaTqNGMQUXDRHFHGeOa4I4QuhGQRZAoQreo7KBE1u/vj3OLruqtqpuqrurqz+t56umqc++tOpei+XDuOfccc3dERESq0iDbFRARkdynsBARkaQUFiIikpTCQkREklJYiIhIUo2yXYF06tChg/fq1Svb1RARqTOKioo2u3vHZPtlLCzMrDswC+gEODDD3W8vs48BtwNnAruBS919cbTtEuA/ol1/7u5/TPaZvXr1orCwMH0nISKS58zsw1T2y2TLYj/wQ3dfbGYtgSIze8ndV8btMxHoGz1GAf8DjDKzdsANQAEhaIrMbK67b8tgfUVEpBIZ67Nw942xVoK77wJWAV3L7HYOMMuDBUAbM+sMnA685O5bo4B4CTgjU3UVEZGq1UoHt5n1AoYCC8ts6gpsiHtdHJVVVl7Re081s0IzK9y0aVO6qiwiInEy3sFtZkcCTwDT3H1nut/f3WcAMwAKCgo0d4mIVNu+ffsoLi7miy++yHZVMqZZs2Z069aNxo0b1+j4jIaFmTUmBMWD7v5kBbuUAN3jXneLykqAk8qUv5KZWopIfVdcXEzLli3p1asXYdxNfnF3tmzZQnFxMb17967Re2TsMlQ00uk+YJW7/6aS3eYC37JgNLDD3TcCLwCnmVlbM2sLnBaViYik3RdffEH79u3zMigAzIz27dsfVsspky2LscDFwDtmtiQq+wnQA8Dd7waeJQybXUsYOntZtG2rmf0MWBQdd5O7b81gXUWknsvXoIg53PPLWFi4+xtAlbXzMD/69yvZNhOYWZ3P3LYNNmyAbt0gz793EZFalVfTfbz/PvToAV27wqRJcPPN8PLLsGtXtmsmIlK1devWMXDgwFo/NlV5Nd1Hv37wve/BwoXh8dRTobxBAxgwAEaNKn307w8NG2a3viIidUVehUWLFnD11eEBsGUL/O1vpeHxxBNw771h25FHQkFBYoB06ZK9uouI7N+/n29+85ssXryYAQMGMGvWLFatWsU111zDZ599RocOHXjggQfo3LkzRUVFTJkyBYDTTjst43XLq7Aoq317mDgxPADc4b33SsNj4UL49a9h//6wvXv3xPAYPhyaN89e/UWk9k2bNo0lS5Yk37EahgwZwm233ZZ0v9WrV3PfffcxduxYpkyZwl133cWcOXN4+umn6dixI4888gjXX389M2fO5LLLLuPOO+9k/PjxXHvttWmtb0XyOizKMoNjjw2Piy8OZV98AW+/HYJjwYLw8/HHw7aGDWHQoBAco0eHn8cdFy5riYikW/fu3Rk7diwAF110Eb/4xS9Yvnw5p556KgAHDhygc+fObN++ne3btzN+/HgALr74Yp577rmM1q1ehUVFmjWDE08Mj5hPPkm8fDV7NvzhD2Fb69YwYkRiC+Soo7JTdxFJv1RaAJlSdnhry5YtGTBgAG+99VZC+fbt22uzWkCejYZKl06d4KtfhZ//HF56KQzJXbkSZs6E88+HzZvDSKuzzw779ukDF1wAt90Gb70VWisiItW1fv36Q8Hw0EMPMXr0aDZt2nSobN++faxYsYI2bdrQpk0b3njjDQAefPDBjNet3rcsUtGgARx/fHhcdlko+/xzWLy4tPXx5pvw8MNhW+PGcMIJia2Pvn1174eIVO24447jrrvuYsqUKfTv35+rr76a008/nR/84Afs2LGD/fv3M23aNAYMGMD999/PlClTMLNa6eC2cF9cfigoKPBsLn700UeJneeLFoVQAWjXDkaOLA2PkSNDB7yIZN+qVas4/vjjs12NjKvoPM2syN0Lkh2rlkUadekSbgacNCm8PnAgXL6KhceCBfDCC2FUFsCXvlQaHqNHh9ZIkybZq7+ISGUUFhkUG001aBB8+9uhbNcuKCwsDZB58yB2ubFpUxg6NPHyVe/eunwlItmnsKhlLVvCySeHB4RWRnFxYutjxgy4PVqtvGPH0stXo0eHkVht2mSv/iJSPyksssws3AzYvTtMnhzK9u2D5csT+z/+/OfSY/r1S2x9DBoUOtVFRDJFYZGDGjcOl6OGDoUrrwxlO3aEDvPYjYPPPgt//GPYdsQRMGxY6Y2Do0aF8NHlKxFJF4VFHdG6NZxySnhAuHy1bl1i6+POO8P0JQBHH53Y+hgxIlwCExGpCYVFHWUWOr979w43CgLs3QvLlpW2PhYuhKefLt2/7My7AwZo5l2RXDB9+nR69uzJtGnTADj99NPp3r0790Yzn/7whz+ka9eurFu3jnnz5mFmNGvWjEcffZTevXvTq1cvWrZsScPoF3r8+PHccccdaa2jwiKPNGkSZtItKICrrgplW7cmTl0yZw7cd1/Y1qJF6cy7sUtYmnlXpPaNHTuWRx99lGnTpnHw4EE2b97Mzp07D22fP38+Z599Nh999BHLli2jQYMGFBcX06JFi0P7vPzyy3To0CFjdcxYWJjZTOAs4FN3L7cqh5ldC3wzrh7HAx2jJVXXAbuAA8D+VG4YkYq1awdnnBEeEC5frV2bePnqt78NneoQVhksO/Nu3N9HEcmAMWPGMH36dABWrFjBwIED2bhxI9u2baN58+asWrWKSZMm0blzZxpEM5l269atVuuYyZbFA8CdwKyKNrr7LcAtAGb2VWB6mXW2T3b3zRmsX71kFqYe6dsXLroolH3xBSxZkhggTzwRtjVsCAMHJrY++vXTzLuSv6ZNC78P6TRkSJg7rjJdunShUaNGrF+/nvnz53PiiSdSUlLCW2+9RevWrRk0aBAXXngh48aN4/XXX2fChAlcdNFFDB069NB7nHzyyYcuQ11yySWHwiddMrkG92tm1ivF3S8AZmeqLlK1Zs1CEIweXVq2aVNieDzySLj/A6BVq/Iz73bqlJ26i+SLMWPGMH/+fObPn88111xDSUkJ8+fPp3Xr1owdO5Zu3bqxevVq5s2bx7x585gwYQKPPfYYEyZMAOrwZahUmVlz4AzgqrhiB140Mwf+4O4zqjh+KjAVoEePHpmsar3SsSOcdVZ4ABw8CGvWJHae/+pXYUoTgJ49E4fuDh0ahvSK1DXZmqF87NixzJ8/n3feeYeBAwfSvXt3fv3rX9OqVSsui2Ywbdq0KRMnTmTixIl06tSJp5566lBYZFrWwwL4KvBmmUtQ49y9xMyOAl4ys3fd/bWKDo6CZAaEiQQzX936qUGDcPmpXz+49NJQtnt34sy7b70VWiAAjRpVPPOuLl+JVGzMmDHceuut9OnTh4YNG9KuXTu2b9/OihUruOeee1i8eDFHH300Xbp04eDBgyxbtozBgwfXWv1yISzOp8wlKHcviX5+amZzgJFAhWEh2dO8OYwbFx4xGzcmXr6aNQt+//uwrU2bxPAYORIy2GoWqVMGDRrE5s2bufDCCxPKYmtvFxYWcsUVV7Bnzx4ARo4cyVVXlV6Qie+zGDx4MLNmVdhdXGMZnaI86rN4pqLRUNH21sAHQHd3/zwqawE0cPdd0fOXgJvc/flkn5ftKcqlvAMHYNWqxABZvjxc1gI45pjyM+82bZrdOkv9oynKszhFuZnNBk4COphZMXAD0BjA3e+OdpsEvBgLikgnYE60vGAj4KFUgkJyU2w01cCBcPnloeyzzxJn3n3lFXjoobCtSZPyM+/26aOpS0SyTYsfSU4oO/NuUVHoE4FwqWrkyNIO9JEjNfOupJdaFlr8SOqIbt3C47zzwuv9+8vPvPvcc6ULRx13XGLrY/Bgzbwrh8fdsTxuwh5uw0BhITmpUaNwI9OQIfCd74SyHTtKL18tWADPPx860CHcKzJ8eGKA9Oihy1eSmmbNmrFlyxbat2+fl4Hh7mzZsoVmzZrV+D10GUrqLHf48MPE1sfixeGOdAg3CpadebdVq+zWWXLTvn37KC4u5ovYX5481KxZM7p160bjMk3wVC9DKSwkr+zbV37m3TVrwjYz6N+//My7jdS+lnpMYSES2bYtcebdhQthy5awrXnzMPPuP/0TnHQSjBkTykTqC4WFSCXc4f33S1sfCxaEy1cHDoRO8lGjQnCcdBKceKLCQ/KbwkKkGnbtgjfeCPd8vPJKGLp74EC476NseGjOK8knCguRw7BzJ7z5ZgiOl18O4XHwYGl4nHxyCI/RoxUeUrcpLETSaOfO8i2PWHiMHp3Y8jiM0YkitU5hIZJBO3YkhsfixSE8mjZNDI/RoxUektsUFiK1aMcOeP310vB4++3S8DjxxNLwGDVK4SG5RWEhkkXbt5dvebiXhkesz2PUKM2yK9mlsBDJIdu3l295uIdWRtmWh8JDapPCQiSHbduWGB5LlpSGx5gxpeExcqTCQzJLYSFSh8TC4+WXQ3gsXZoYHrHLViNGKDwkvRQWInXY1q2JLY9YeBxxRPmWR5Mm2a2r1G1ZDwszmwmcBXxa0bKqZnYS8DRhWVWAJ939pmjbGcDtQEPgXne/OZXPVFhIvtq6FV57LTE8IITH2LGl4TFihMJDqicXwmI88Bkwq4qw+Dd3P6tMeUNgDXAqUAwsAi5w95XJPlNhIfXFli2J4bFsWShv3jwxPAoKFB5StayvlOfur5lZrxocOhJY6+7vA5jZw8A5QNKwEKkv2reHSZPCA0rDI9bncf31oTwWHrE+j4ICrSgoNZPtmfxPNLOlwEeEVsYKoCuwIW6fYmBUZW9gZlOBqQA9evTIYFVFclfZ8Ni8ObHl8ZOfhPIWLcq3PBQekopshsVioKe7f2ZmZwJPAX2r+ybuPgOYAeEyVHqrKFI3degAX/taeABs2lR5eIwbVxoew4crPKRiWQsLd98Z9/xZM/u9mXUASoDucbt2i8pEpIY6doTzzgsPKA2P2GWr664L5UcemRgew4YpPCTIWliY2dHAJ+7uZjYSaABsAbYDfc2sNyEkzgcuzFY9RfJR2fD49NPElsePfxzKy4bH8OFahra+ytjXbmazgZOADmZWDNwANAZw97uBycB3zWw/8A/gfA9Ds/ab2VXAC4ShszOjvgwRyZCjjoLJk8MDQni8+mrF4RFbgjbW8lB41A+6KU9Ekvrkk8TwWLUqlLdsmRgeQ4cqPOqarN9nkQ0KC5HaEQuPWJ/Hu++G8lh4xIbqDhmi8Mh1CgsRqTUff5zY8oiFR6tWpS2Pk08O4dGwYRYrKuUoLEQkazZuTAyP1atDeatWMH586WUrhUf2KSxEJGfEwiN22WrNmlDeunVieJxwgsKjtiksRCRnffRRYsujbHjE+jwGD1Z4ZJrCQkTqjJKSxPB4771Q3qZNYstD4ZF+CgsRqbNKSkqD45VXYO3aUN62bfnwaNAga9XMCwoLEckbxcWJfR5//3soj4VH7LLVoEEKj+pSWIhI3tqwIfGyVXx4fPnLcOaZcO65YVoTqZrCQkTqjfXrS8Nj3jxYty70bZx0UpjCZNIk6NQpy5XMUQoLEamX3MOys48/Do89FkZaNWgQLldNnhymbe/cOdu1zB0KCxGp99xh+fLS4Fi1CszCTLqTJ4dZd7t2zXYts0thISJSxsqVpcGxfHkoGzMGvv71EBzdu1d9fD5SWIiIVOHdd+GJJ0JwLF0aykaNKg2OXr2yWr1ao7AQEUnRe++FFsfjj8PixaGsoCAEx+TJ0KdPduuXSWkJCzP7Wgqf9YW7P1udymWKwkJEDtf775cGx6JFoWzo0NLg6Ns3u/VLt3SFxRbgacCqeI/x7n5M9auYfgoLEUmndevCparHH4cFC0LZCSeUrirYr19Wq5cW6QqLP7n7RUk+qMJ9zGwmcBbwqbsPrGD7N4EfEYJoF/Bdd18abVsXlR0A9qdyIqCwEJHM2bChNDjefDOUDRxYGhwDBmS3fjWV9T4LMxsPfAbMqiQsxgCr3H2bmU0EbnT3UdG2dUCBu2+uzmcqLESkNpSUwJNPhuB4/fUwRPf440NofP3rIUSsqusxOSStYWFmDYF/BnoBhxZJdPffJDmuF/BMRWFRZr+2wHJ37xq9XofCQkTqgI0bYc6cEByvvgoHD8Kxx5b2cZxwQm4HR6phkeqUW/8HXAq0B1rGPdLlcuC5uNcOvGhmRWY2taoDzWyqmRWaWeGmTZvSWCURkeQ6d4bvfS9MM/LRR3D33eF+jV/+MnSMH3ssXHcdFBWFFkhdlWrLYpm7D672m6fQsjCzk4HfA+PcfUtU1tXdS8zsKOAl4Gp3fy3Z56llISK5YtMmePrpcB/HX/8KBw5A796lfRwjRuRGiyPdLYvnzOy0w6xTOWY2GLgXOCcWFADuXhL9/BSYA4xM92eLiGRSx47w7W/DCy/AJ5/AffeF0VO33RZu/vvSl+DOO2H37mzXNDWphsUCYI6Z/cPMdprZLjPbeTgfbGY9gCeBi919TVx5CzNrGXsOnAYsP5zPEhHJpvbtYcoUePbZEBwPPABHHw1XXw09esBPfwqbq9VDW/tSDYvfACcCzd29lbu3dPdWVR1gZrOBt4DjzKzYzC43syvN7Mpol/8k9IH83syWmFns+lEn4A0zWwr8Dfizuz9f3RMTEclFbdvCJZeE4bdvvBHmprrxxhAaV18d7u3IRan2WbwGnOTuBzNfpZpTn4WI1EUrV8Itt8CDD4bRVP/yL3DttaGDPNPS3WfxPvCKmV1nZtfEHodXRRERAejfH+6/P0w1Mn06PPMMDBsGp58eOsdzYRRVqmHxAfBXoAmZGTorIlLvdesWWhjr14eht0uXwimnhJFTjz4K+/dnr26adVZEJEd98QX87/+GAHnvvTD77Q9/CJddBkcckZ7PSMtlKDO7MYUPSrqPiIhUX7NmcMUVYYW/J58Mw3G//33o2RN+9jPYsiX5e6RLsokEiwkjoSrdBbjC3XNi7kW1LEQkn7mHuah+9aswDLd58xAm06eHAKmJdHVw30NiH0XZx5HRPiIikmFmMH48/PnPsGxZuBP8rrvgmGPg8sthx44Mfrb6LERE6q716+E3vwl3g/foAbNnhzvEU5XuobMiIpKDevQIU4i89lq4R2PcOLj55vA8nRQWIiJ5YMwYWLIEJk0Ks9yedlqYPj1dFBYiInmiTRt45BG45x6YPx8GDw4d4elQo7Aws++Z2TfMrFHyvUVEpLaYhdlui4qgSxf4538Oo6X27Dm8961py8KAcYRZY0VEJMccfzwsXAhXXRX6NEaPhtWra/5+NQoLd7/L3a9297Nr/tEiIpJJzZrB734XFmFavz7MN/XEEzV7r5QuI5lZU+A8yq/BfVPNPlZERGrL2WeH+zK+9rUwPfqwYWHVvupItWXxNHAOsB/4PO4hIiJ1QNeuYYnXBg3CDXzVHVqbagd1N3c/o9q1ExGRnNGjB9x6K3znOzBjBlx5ZfJjYlJtWcw3s0HVrZiZzTSzT82swmVRLbjDzNaa2TIzGxa37RIzey96XFLdzxYRkfKuuCJMe37ttdVblS/VsBgHFJnZ6ugf9XfMbFkKxz0AVNUimQj0jR5Tgf8BMLN2wA3AKGAkcIOZtU2xriIiUgkzuPfe8Pzyy1M/LtXLUBOrXSPA3V8zs15V7HIOMMvDBFULzKyNmXUGTgJecvetAGb2EiF0ZtekHiIiUqpnz3A5Ku2Xodz9Q6AN8NXo0SYqO1xdgQ1xr4ujssrKyzGzqWZWaGaFmzZtSkOVRETy39Sp0K8ai0ukFBZm9q/Ag8BR0eNPZnZ1TSqYbu4+w90L3L2gY8eO2a6OiEidYAZDh6a+f6qXoS4HRrn75+FD7FfAW8DvqlvBMkqA7nGvu0VlJYRLUfHlrxzmZ4mISJwBA1LfN9UObgMOxL0+EJUdrrnAt6JRUaOBHe6+EXgBOM3M2kYd26dFZSIikibDh6e+b6oti/uBhWY2J3p9LnBfsoPMbDahhdAhWqL1BqAxgLvfDTwLnAmsBXYDl0XbtprZz4BF0VvdFOvsFhGR9DijGnfPpbxSXnQPxLjo5evu/na1a5ZhWilPRKR6Ul0pr8qWhZm1cved0X0P66JHbFs7/W9fRKR+SHYZ6iHgLKAIiG+CWPS6T4bqJSIiOaTKsHD3s6Kf1ZyfUERE8kmq91n8NZUyERHJT8n6LJoBzQmjmdpSOly2FZXcUS0iIvknWZ/Fd4BpQBdCv0UsLHYCd2awXiIikkOS9VncDtxuZle7++HerS0iInVUSjflufvvzGwg0B9oFlc+K1MVExGR3JHqGtw3EO7E7k+463oi8AagsBARqQdSnRtqMjAB+NjdLwNOAFpnrFYiIpJTUg2Lf7j7QWC/mbUCPiVxtlgREcljqU4kWGhmbYB7CKOiPiNMUS4iIvVAqh3c34ue3m1mzwOt3D2VNbhFRCQPpHoH91wzu9DMWrj7OgWFiEj9kmqfxa8J05OvNLPHzWxydHe3iIjUA6lehnoVeNXMGgJfAa4AZhKm/RARkTyXassCMzsCOA+4EhgB/DGFY84ws9VmttbMflzB9t+a2ZLoscbMtsdtOxC3bW6q9RQRkfRL9aa8R4GRwPOEOaFejYbSVnVMQ+Au4FSgGFhkZnPdfWVsH3efHrf/1cDQuLf4h7sPSfVEREQkc1IdOnsfcIG7H6jGe48E1rr7+wBm9jBwDrCykv0vIKzRLSIiOSbVy1CvA9eZ2QwAM+trZmclOaYrsCHudTGVTGtuZj2B3sC8uOJmZlZoZgvM7NwU6ykiIhmQaljcD+wFxkSvS4Cfp7Ee5wOPl2m59IwWEb8QuM3MjqnoQDObGoVK4aZNm9JYJRERiUk1LI5x9/8G9gG4+25K17aoTAmJU4J0i8oqcj4wO77A3Uuin+8Dr5DYnxG/3wx3L3D3go4dOyapkoiI1ESqYbE3Gg3lANH/8vckOWYR0NfMeptZE0IglBvVZGb9gLbETR9iZm3NrGn0vAMwlsr7OkREJMNS7eC+gTASqruZPUj4x/vSqg5w9/1mdhXwAtAQmOnuK8zsJqDQ3WPBcT7wsLt73OHHA38ws4OEQLs5fhSViIjULkv8N7qKHc3aA6MJl58WuPvmTFasJgoKCrywsDDb1RARqTPMrCjqH65SlS0LMxtWpmhj9LOHmfVw98U1raCIiNQdVbYszOzlKo51d/9K+qtUc40bN/ZTTz2VgoICRowYwYgRIzj66KOzXS0RkZyVassi5ctQdUGHDh28S5curFixgoMHww3m3bp1OxQcI0aMoKCggDZt2mS5piIiuSFdl6H+PRoyi5l93d0fi9v2C3f/yeFXNX169epFYWEhn3/+OW+//TaLFi069JgzZ86h/fr27ZsQIEOHDqV58+ZZrLmISG5LdhlqsbsPK/u8ote5oKoO7m3btlFYWJgQICUl4baPhg0bMmDAgIQAGTRoEI0bN67N6ouI1Lq0XIYys7fdfWjZ5xW9zgXVHQ21cePGhPBYtGgRW7duBaBp06YMGTIkIUCOO+44GjRIeaJeEZGcl66wyJuWRSrcnQ8++CAhPIqKivj8888BaNmyJcOHD08IkJ49e2KW7GZ2EZHclK6wOAB8Tri34ghgd2wT0Mzdc+o6TSbuszhw4ADvvvtuQoAsXbqUvXv3AtChQ4eE0VcagSUidUm9HA1VWzfl7d27l2XLlrFo0aJD/SAagSUidZHCopZVNAJr7dq1h7ZrBJaI5CKFRQ7QCCwRyXUKixylEVgikksUFnVEshFYRx55ZLkRWL169dIILBFJC4VFHRY/Ait2GWvJkiUagSUiaaewyDN79+7lnXfeSWiBaASWiBwuhUU9oBFYInK4FBb1lEZgiUh15ERYmNkZwO2EZVXvdfeby2y/FLgFKImK7nT3e6NtlwD/EZX/3N3/mOzzFBYVSzYC64QTTig3Aqthw4ZZrrWI1Iash4WZNQTWAKcCxcAi4IL4tbSjsChw96vKHNsOKAQKAAeKgOHuvq2qz1RYpKbsCKzCwkKKior47LPPAI3AEqlP0rKexWEaCax19/ejCj0MnAOsrPKo4HTgJXffGh37EnAGMDtDda1XzIw+ffrQp08fvvGNbwBhBNbq1asTWh933HGHRmCJCJDZsOgKbIh7XQyMqmC/88xsPKEVMt3dN1RybNeKPsTMpgJTAXr06JGGatdPDRs2pH///vTv359LLrkEqHgE1osvvqgRWCL1UCbDIhX/B8x29z1m9h3gj0C11vV29xnADAiXodJfxfqrSZMmDB8+nOHDh3PllVcCFY/A0iqEIvkvk2FRAnSPe92N0o5sANx9S9zLe4H/jjv2pDLHvpL2Gkq1tWjRgnHjxjFu3LhDZWVHYL366qs89NBDADRo0KDCEVhNmjTJ1imISA1ksoO7EeHS0gTCP/6LgAvdfUXcPp3dfWP0fBLwI3cfHXVwFwGxxZUWEzq4t1b1mergzh0agSVSN2R9NFRUiTOB2whDZ2e6+3+Z2U1AobvPNbNfAmcD+4GtwHfd/d3o2CnAT6K3+i93vz/Z5ykscpe7s27dunJzYGkElkh25URY1DaFRd1S0QgszYElUrsUFlInaQ4skdqlsJC8kcocWPEtkGHDhmkElkiKFBaS16qaA0sjsERSp7CQeic2Ais+RLZsCaOzNQJLpGIKC6n3NAJLJDmFhUgFko3Aat++fUJ4aASW5DuFhUiKUh2BFetELygooG3btlmutUh6KCxEDkOyEVhf+tKXys2B1aJFiyzWWKRmFBYiabZt2zaKiooSAqS4uBjQCCypuxQWIrVg48aN5YbwagSW1CUKC5Es0AgsqWsUFiI5IpURWGXnwOrcuXOWay31hcJCJIeVHYFVWFjI8uXLD43A6tq1a7k5sDQCSzJBYSFSx+zevbvcCKz33nvv0HaNwJJMUFiI5AGNwJJMU1iI5KmPP/643CqEsRFYTZo0KTcCq1+/fhqBJZXKibAwszOA2wkr5d3r7jeX2X4N8G3CSnmbgCnu/mG07QDwTrTrenc/O9nnKSykPq1TdxwAAAyzSURBVEplBNawYcMSAqR3794agSVADoSFmTUkrMF9KlBMWIP7AndfGbfPycBCd99tZt8FTnL3b0TbPnP3I6vzmQoLkSA2Aiv+HpAlS5awZ88eQCOwpFQuhMWJwI3ufnr0+joAd/9lJfsPBe5097HRa4WFSBrt3buX5cuXl5sD68CBA4BGYNVXuRAWk4Ez3P3b0euLgVHuflUl+98JfOzuP49e7weWEC5R3ezuT1Vy3FRgKkCPHj2Gf/jhh2k/F5F8VZ0RWAUFBQwbNkwjsPJMqmHRqDYqk4yZXQQUAF+OK+7p7iVm1geYZ2bvuPvfyx7r7jOAGRBaFrVSYZE80bx5c8aOHcvYsWMPlZUdgfX6668ze/ZsIIzA6t+/f0ILZPDgwRqBVQ9kMixKgO5xr7tFZQnM7BTgeuDL7r4nVu7uJdHP983sFWAoUC4sRCS92rZtyymnnMIpp5xyqKzsCKy5c+dy//33AxqBVV9k8jJUI0IH9wRCSCwCLnT3FXH7DAUeJ1yuei+uvC2w2933mFkH4C3gnPjO8Yqoz0Kkdrg7H374YbkRWLt27QI0AqsuyXqfRVSJM4HbCENnZ7r7f5nZTUChu881s78Ag4CN0SHr3f1sMxsD/AE4CDQAbnP3+5J9nsJCJHsOHjxY4RxYGoGV23IiLGqbwkIkt6Q6Ait+FcJ27dpludb1i8JCRHJSshFYxxxzTELrQyOwMkthISJ1xvbt28vNgbVhwwZAI7AyTWEhInXaxx9/XG4Vws2bNwMagZVOCgsRySvJRmC1aNGC4cOHJ3Si9+nTRyOwklBYiEjeSzYCq127duVGYHXp0iXLtc4tCgsRqZf27t3LihUrEgJk+fLlh0ZgdenSpdwcWPV5BJbCQkQksnv3bpYsWZIQIGvWrDm0vT6PwFJYiIhUQSOwAoWFiEg1ffLJJ+VWIcz3EVgKCxGRw1Qf5sBSWIiIZEB15sCK/czlEVgKCxGRWhI/B1bsRsK6MgJLYSEikkV1ZQSWwkJEJMdUdwTWoEGDaNq0aUbrpLAQEakDko3AGjx4cEKAHH/88WkdgaWwEBGpg+JHYMX6PwoLCxPmwCo7Autw5sBSWIiI5ImDBw+yZs2ahNbH22+/nZY5sHIiLMzsDOB2wrKq97r7zWW2NwVmAcOBLcA33H1dtO064HLgAPADd38h2ecpLESkvti3b1+5VQhrMgIr62FhZg2BNcCpQDGwCLjA3VfG7fM9YLC7X2lm5wOT3P0bZtYfmA2MBLoAfwGOdfcDVX2mwkJE6rOajMA68sgjUwqLRhms90hgrbu/D2BmDwPnACvj9jkHuDF6/jhwp4ULb+cAD7v7HuADM1sbvd9bGayviEid1rx5c8aMGcOYMWMOlZUdgfXmm2/y8MMPA2EEVqoyGRZdgQ1xr4uBUZXt4+77zWwH0D4qX1Dm2K4VfYiZTQWmAvTo0SMtFRcRyRdt2rRhwoQJTJgw4VBZ/Aism266KaX3ST1WcpS7z3D3Ancv6NixY7arIyKS8zp16sRZZ53FT3/605SPyWRYlADd4153i8oq3MfMGgGtCR3dqRwrIiK1JJNhsQjoa2a9zawJcD4wt8w+c4FLoueTgXkeetznAuebWVMz6w30Bf6WwbqKiEgVMtZnEfVBXAW8QBg6O9PdV5jZTUChu88F7gP+N+rA3koIFKL9HiV0hu8Hvp9sJJSIiGSObsoTEanHUr3Pos53cIuISOYpLEREJCmFhYiIJKWwEBGRpPKqg9vMNgEfZrselegAbM52JQ5TPpwD5Md56BxyQz6cw3Hu3jLZTpmc7qPWuXvO3sJtZoWpjDjIZflwDpAf56FzyA35cg6p7KfLUCIikpTCQkREklJY1J4Z2a5AGuTDOUB+nIfOITfUm3PIqw5uERHJDLUsREQkKYWFiIgkpbDIADP7upmtMLODZlYQV36qmRWZ2TvRz6/EbRsela81szui5WWzprJziLZdF9VztZmdHld+RlS21sx+XPu1rpyZDTGzBWa2xMwKzWxkVG7Rn/daM1tmZsOyXdeqmNnVZvZu9N38d1x5hd9JLjOzH5qZm1mH6HWd+S7M7Jboe1hmZnPMrE3ctjrzXVTrd9bd9UjzAzgeOA54BSiIKx8KdImeDwRK4rb9DRgNGPAcMDFHz6E/sBRoCvQG/k6Ygr5h9LwP0CTap3+2v4u4er8Y+zMFzgReiXv+XPTnPhpYmO26VnEOJwN/AZpGr4+q6jvJdn2TnEt3wvIFHwId6uB3cRrQKHr+K+BXde27qO7vrFoWGeDuq9x9dQXlb7v7R9HLFcAR0QJPnYFW7r7Aw7c4Czi3FqtcTmXnAJwDPOzue9z9A2AtMDJ6rHX39919L/BwtG+ucKBV9Lw1EPsezgFmebAAaBN9H7nou8DN7r4HwN0/jcor+05y2W+Bfyd8LzF15rtw9xfdfX/0cgFhNU+oW99FtX5nFRbZcx6wOPrF7woUx20rjspyUVdgQ9zrWF0rK88V04BbzGwDcCtwXVSe6/WOdyzwT2a20MxeNbMRUXldOgfM7BxCq3ppmU116jziTCG0iKBunUO16ppX033UJjP7C3B0BZuud/enkxw7gNB0PS0TdUvV4ZxDLqrqfIAJwHR3f8LM/oWwSuMptVm/VCQ5h0ZAO8IlmhHAo2bWpxarl7Ik5/ETsvx3PxWp/H6Y2fWE1TwfrM26ZYPCoobcvUb/0JhZN2AO8C13/3tUXEJpM5boecnh1TC5Gp5DCeF6c0x8XSsrrxVVnY+ZzQL+NXr5GHBv9Lyq86l1Sc7hu8CT0aXKv5nZQcJEdjl1DlD5eZjZIMK1/KXRGI5uwOJowEFOnUey3w8zuxQ4C5gQfSeQY+eQRLXqqstQtSgaMfFn4Mfu/mas3N03AjvNbHQ0CupbQK7+z34ucH7U19Ib6EvonF8E9DWz3mbWhLCe+tws1rOsj4AvR8+/ArwXPZ8LfCsaiTMa2BF9H7noKUInN2Z2LKFTcjOVfyc5x93fcfej3L2Xu/ciXPoY5u4fU4e+CzM7g9Dncra7747bVGe+C6r7O5vtHvl8fACTCL8Ee4BPgBei8v8APgeWxD1iI1oKgOWE0Ql3Et1dn2vnEG27PqrnauJGbRFGs6yJtl2f7e+hzPmMA4oIIz4WAsOjcgPuiur8DnEjv3LtQQiHP0V/TxYDX0n2neT6A1hH6WiouvRdrCVc74/9Ht9dF7+L6vzOaroPERFJSpehREQkKYWFiIgkpbAQEZGkFBYiIpKUwkJERJJSWIiISFIKC8lrZnajmf1b9PwmM6vqDulzzax/7dWu5szsUjPbZGb3Jt874bhbzOzj2J+JSKo03YfUG+7+n0l2ORd4BlhZC9VJh0fc/arqHODu15rZ55mqkOQvtSwk75jZ9Wa2xszeIKzJESt/wMwmR89vNrOV0eI1t5rZGOBswsy0S8zsGDO7wswWmdlSM3vCzJrHvc8dZjbfzN6PvWe07UcWFrFaamY3R2XHmNnzFha8et3M+iWp/3Qzmxk9H2Rmy2OfXcUxl5rZU2b2kpmtM7OrzOwaM3vbwqJP7Wr4xykCqGUhecbMhhPmuBlC+Pu9mDDNR/w+7QnTmfRzdzezNu6+3czmAs+4++PRftvd/Z7o+c+By4HfRW/TmTCFSD/CfDqPm9lEwnoAo9x9d9w/0DOAK939PTMbBfyeMD9VZW4HXjGzSYSpI77jifMPVWYgYYGtZoTpKH7k7kPN7LeE+cZuS+E9RCqksJB880/AnNg/rlEAlLUD+AK4z8yeIVx6qsjAKCTaAEcSVnaLecrdDwIrzaxTVHYKcH/ss919q5kdCYwBHrPSlXKbVnUC7n4wmtF0GfAHj5t0MomX3X0XsMvMdgD/F5W/AwxO8T1EKqSwkHrH3fdHU2JPACYDV1Hx//QfAM5196XRP94nxW3bE/e8qvXSGwDb3X1INavZF/gM6FKNY+LrdDDu9UH0uy6HSX0Wkm9eA841syPMrCXw1bI7RP/bb+3uzwLTgROiTbuAlnG7tgQ2mllj4JspfPZLwGVxfRvt3H0n8IGZfT0qMzM7IXo+ycx+WUH9WgN3AOOB9vF9IiLZorCQvOLui4FHCFORP0eYs7+slsAzZrYMeAO4Jip/GLg26hQ+Bvh/hOnM3wTeTeGznyf0XxSa2RIgNjz1m8DlZraUsPZ6bJ3jY4CdFbzVb4G73H0NoZ/kZjM7Ktnni2SSpigXyRIz+xNhqddNNTj2UsJ6D9UaOhsdeyPwmbvfWt1jpf5Sy0IkS9z9opoEReQfwMSa3JQHXERYhEskZWpZiIhIUmpZiIhIUgoLERFJSmEhIiJJKSxERCSp/w/hpPOq+k34hwAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        },
        {
          "output_type": "stream",
          "text": [
            "distance to theupstream bc = -136.11\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}