{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 7-2: Mixing and dye-dilution discharge\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lab gives an example of how to use a bulk injection of flourescent dye, measure its concentration over time, and determine discharge in a river.\n",
    "\n",
    "This uses the third of four attempts at Glen Aulin in June 2005.  Use the excel file and the worksheet from the main module to assess the entire situation more qualitatively.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing python packages you'll need for this lab:\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sec_after_inj</th>\n",
       "      <th>min_after_inj</th>\n",
       "      <th>DateTime</th>\n",
       "      <th>Rwt_in_stream_ug_per_L</th>\n",
       "      <th>dt_sec</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15:39:02</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15:39:12</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>15:39:22</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sec_after_inj  min_after_inj  DateTime  Rwt_in_stream_ug_per_L  dt_sec\n",
       "0              0            0.0  15:39:02                     0.0     NaN\n",
       "1              0            0.0  15:39:12                     0.0     0.0\n",
       "2              0            0.0  15:39:22                     0.0     0.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('Dye_data_GA050630c.csv', comment='#')\n",
    "df.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the weight of dye in the stream as a function of minutes after the injection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f4a7342eb80>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAFzCAYAAACeg2ttAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA820lEQVR4nO3deXhU5fk+8PuZLGQCBKKACSiCIAJllSTQstpqFUUQRA0EWmVNQMSFFtx+WqgKFVxQCKKgtiyuiIi02hZQYgXCIgglQb8JIEwSQMMWsuf5/TGTNEISssyZdzJzf65rLiaTOXNuptTcOfMuoqogIiIi8hSb6QBERETkX1g+iIiIyKNYPoiIiMijWD6IiIjIo1g+iIiIyKNYPoiIiMijAk0HKK9Zs2bapk0b0zGIiIjIDXbu3HlSVZtf+LhXlY82bdpgx44dpmMQERGRG4jI4Yoe58cuRERE5FEsH0RERORRLB9ERETkUSwfRERE5FEsH0RERORRLB9ERETkUSwfRERE5FEsH0RERORRLB9ERETkUSwf5HUyMjIwcOBAZGZmmo5CREQWYPkgrzNnzhwkJSVh9uzZpqMQEZEFRFVNZygTFRWl3NvFf9ntduTl5V30eEhICHJzcw0kIiKiuhCRnaoadeHjvPJBXiMtLQ2jR49GYKBzv8Pg4GDExcUhPT3dcDIiInInlg/yGpGRkcjLy0NRUREAoKCgAAUFBYiIiDCcjIiI3Inlg7xGcXExNm7ciNDQUHz00UcIDQ3F+vXrcfLkSdPRiIjIjVg+yGssWbIEp06dwrJly3DHHXdg06ZNKCkpwV133YXCwkLT8YiIyE1YPsgrZGZm4rHHHsONN96Ie+65BwAQExOD119/HZs3b8ZDDz1kOCEREblLoOkARAAwY8YM5OXlYdGiRRCRssfHjh2LPXv2YMGCBejevTsmTpxoMCUREbkDr3yQcRs3bsTKlSsxc+ZMdOjQ4aLvz507F7/97W8xdepUJCUlGUhIRETuxHU+yKiCggJ069YNhYWF2LdvH+x2e4XPy87ORu/evXH69GkkJyejdevWHk5KREQ1xXU+yCvNnz8fqampePXVVystHgAQHh6Ojz/+GLm5ubjjjjtw/vx5D6YkIiJ3YvkgY9LT0zFnzhyMGDECgwcPvuTzO3XqhFWrVuGbb77B+PHj4U1X7YiIqPpYPsgIVcUDDzyAgIAAvPTSS9U+bsiQIXjmmWfwzjvvYN68edYFJCIiy3C2Cxmxbt06rF+/Hs8//zyuuuqqGh07a9Ys7N27F4899hi6dOmCIUOGWJSSiIiswAGn5HE5OTno3LkzwsLCsGvXLgQFBdX4Nc6fP4/+/fvju+++w7Zt29CpUycLkhIRUV1wwCl5jTlz5uDIkSNITEysVfEAULYEu91ux7Bhw5Cdne3mlEREZBWWD/Ko/fv3Y8GCBbjvvvvQr1+/Or1W69at8eGHH+LQoUMYNWoUiouL3ZSSiIisxPJBHqOqmDJlCho3buy2waL9+vXDokWL8Nlnn2HWrFlueU0iIrIWywd5zN/+9jd8+eWXmDdvHpo3b+621504cSKmTp2K+fPn429/+xsyMjIwcOBAZGZmuu0cRETkPhxwSh6RnZ2N6667Du3atcNXX30Fm829vbewsBC//e1v8fXXX+O2227D2rVrMXnyZCxevNit5yEiouqrbMApywd5REJCApYuXYqdO3eiR48elpwjJCQE+fn5FT6em5tryTmJiKhynO1Cxmzfvh2vvfYapk2bZlnxAJwrppZfKTU0NBRxcXFIT0+37JxERFRzLB9kqeLiYsTHxyMiIgKzZ8+29FyRkZG4+uqry77Oy8tDWFgYIiIiLD0vERHVDMsHWWrx4sXYvXs3XnzxRYSFhVl+vqysLAwbNgwAcMMNN3DQKRGRF+KYD7JMRkYGOnbsiN69e+Ozzz6DiHjkvKqKa6+9Fq1bt8bGjRs9ck4iIroYx3yQR2VkZKB79+7Izc3Fq6++6rHiAQAigvvuuw+bNm3C//3f/3nsvEREVD0sH2SJyZMn48SJE+jWrRs6dOjg8fP//ve/h81mw1tvveXxcxMRUdX4sQu5ld1uR15e3kWPm5jueuutt+Lbb7/FoUOHEBAQ4NFzExERP3YhD0lLS0NsbGzZ1yanu44bNw5Hjx7Fv/71L4+fm4iIKsfyQW4VGRlZdj8oKMjodNfbb78dl19+OZYvX+7xcxMRUeUCTQcg33Po0CEAwLJly7B161ZkZGQYydGgQQOMGTMGiYmJ+PHHH3H55ZcbyUFERD/HKx/kdiNGjAAADBkyBIsWLcKaNWuMZRk/fjwKCgqwcuVKYxmIiOjnWD7I7VJTU9G8eXOEh4ebjoKuXbsiKioKy5YtgzcNriYi8mcsH+R2qamp6Nixo+kYZcaNG4e9e/di9+7dpqMQERFYPsgCqampuO6660zHKDNq1CiEhIRg2bJlpqMQERFYPsjNfvrpJ5w4ccKrykfTpk0xYsQIrFq1yuNrjRAR0cVYPsitUlNTAcCrygfgHHh66tQprF271nQUIiK/Z3n5EJEAEdktIuutPheZV1o+vGnMBwAMGjQIbdq04ZofRERewBNXPqYDOOCB85AXSElJQVBQENq2bWs6ys/YbDbcd999+Pe//122DgkREZlhafkQkSsB3AbgDSvPQ94jNTUV7dq1Q2Cg961fd++99wIAN5sjIjLM6isfLwH4I4CSyp4gIpNEZIeI7Dhx4oTFcchq3jbNtrzWrVvjpptuwptvvomSkkr/SRIRkcUsKx8iMgTAcVXdWdXzVHWpqkapalTz5s2tikMeUFRUhO+//97rBpuWN27cOBw5cgQbN240HYWIyG9ZeeWjL4ChInIIwDsAfi0iKyw8HxmWnp6OwsJCry4fw4YNQ3h4OAeeEhEZZFn5UNVHVfVKVW0DIBbARlUdY9X5yDxvnWZbXkhICOLi4rBmzRpkZ2ebjkNE5Je4zge5TX0oH4Dzo5f8/HysWrXKdBQiIr/kkfKhqptVdYgnzkXmpKSkoFmzZl6/dX3Pnj3Rs2dPfvRCRGQIr3yQ23jbni5VGTduHHbt2oVvvvnGdBQiIr/D8kFu483TbC80evRoBAcH48033zQdhYjI77B8kFtkZ2fj+PHj9ebKx2WXXYbhw4djxYoVyM/PNx2HiMivsHyQW9SXwabljRs3Dj/99BM+/vhj01GIiPwKywe5RX0sH7/5zW/QunVrDjwlIvIwlg9yi9TUVAQGBuKaa64xHaXaAgICcO+99+Lzzz/HDz/8YDoOEZHfYPkgt0hJSUG7du0QFBRkOkqN3HvvvVBVvP3226ajEBH5DZYPcov6NM22vLZt2+LXv/41li9fzs3miIg8hOWD6qy4uBjff/99vZlme6Fx48YhPT0dX3zxhekoRER+geWD6uzQoUMoKCiol1c+AGDEiBFo0qQJB54SEXkIywfVWUpKCoD6NdOlPLvdjtGjR+P9999H3759kZmZaToSEZFPY/mgOquP02wvVLrZ3Ndff43Zs2ebjkNE5NNYPqjOUlNTcfnll6NZs2amo9SK3W5HdHQ0AEBVkZiYCBGB3W43nIyIyDexfFCdpaSk1OurHmlpaRg9enTZNOGQkBDExcUhPT3dcDIiIt/E8kF1Vl+n2ZaKjIxEWFgYiouLAQB5eXkICwtDRESE4WRERL6J5YPq5PTp08jKyqq302xLZWVlIT4+HkOHDkVgYCCOHj1qOhIRkc8KNB2A6jdfGGwKAGvWrAEAbNu2DevWrcOtt95qOBERke/ilQ+qk/o+zfZCMTEx6NmzJxITE6GqpuMQEfkklg+qk9TUVAQEBNSrDeWqIiKIj4/H3r17sXXrVtNxiIh8EssH1UlqairatWuH4OBg01HcZvTo0WjcuDGWLFliOgoRkU9i+aA6qe8zXSrSqFEjjB07Fu+++y5+/PFH03GIiHwOywfVWnFxMb777jufKx8AMHnyZOTn5+Ptt982HYWIyOewfFCtHT58GPn5+fV+mm1FunXrhl/96ldYsmQJB54SEbkZywfVmq9Ms61MQkICvvvuO2zcuNF0FCIin8LyQbXma9NsLzRy5EhcdtllHHhKRORmLB9Ua6mpqQgPD6+3G8pdSkhICO677z6sXbsWGRkZpuMQEfkMlg+qtdTUVHTs2BEiYjqKZSZPnoyioiIsW7bMdBQiIp/B8kG15ovTbC907bXX4sYbb8TSpUvLNp4jIqK6YfmgWjlz5gwyMjJ8vnwAzoGnP/zwAzZs2GA6ChGRT2D5oFopnenii9NsL3T77bcjMjKSA0+JiNyE5YNqxden2ZYXFBSECRMm4O9//zsOHTpkOg4RUb3H8kG1kpKSgoCAALRr1850FI+YOHEiRARLly41HYWIqN5j+aBaSU1NRdu2bX1qQ7mqXHXVVRgyZAiWLVuGgoIC03GIiOo1lg+qldJptv4kPj4ex48fx9q1a01HISKq11g+qMZ8eUO5qtx8881o06YNEhMTTUchIqrXWD6oxo4cOYK8vDy/Kx82mw2TJ0/G5s2by5aWJyKimmP5oBrzp2m2Fxo3bhyCgoLw2muvmY5CRFRvsXxQjfnTNNsLtWjRAnfeeSfeeustnD9/3nQcIqJ6ieWDaiwlJQVNmzZF8+bNTUcxIj4+HqdOncJ7771nOgoRUb3E8kE1Vrqniy9vKFeVAQMGoFOnThx4SkRUSywfVGP+OM22PBFBfHw8tm/fjl27dpmOQ0RU77B8UI2cPXsWDofDL8d7lDd27FjY7XYOPCUiqgWWD6oRfx5sWl54eDhiY2OxcuVKnDlzxnQcIqJ6heWDasSfp9leKCEhATk5OVixYoXpKERE9QrLB9VIamoqbDab32woV5WoqChcf/31WLJkCVTVdBwionqD5YNqJCUlBW3btkWDBg1MRzGudODpt99+i//85z+m4xAR1RssH1QjpdNsyWnUqFEICwvDkiVLTEchIqo3WD6o2kpKSvDdd99xvEc5jRo1wtixY/H+++/j5MmTpuMQEdULLB9UbT/88ANyc3N55eMC8fHxyM/Px8KFCzFw4EBkZmaajkRE5NVYPqjaSndyZfn4uS5duqBfv354+eWXkZSUhNmzZ5uORETk1Vg+qNo4zbZidrsdSUlJOHPmDEpKSpCYmAgRgd1uNx2NiMgrsXxQtaWmpqJJkyZo0aKF6SheJS0tDffcc0/Z16GhoYiLi0N6errBVERE3ovlg6otJSXFrzeUq0xkZCTCw8PL3pfc3FyEhYUhIiLCcDIiIu/E8kHVxmm2lcvKysKECRPQtGlTXHXVVRx0SkRUhUCrXlhEQgB8CaCB6zwfqOpTVp2PrHXu3DkcO3aM4z0qsWbNGgBA27Zt8dhjj+GDDz4wnIiIyHtZeeUjH8CvVbU7gB4AbhGRPhaejyx08OBBAJzpcilTp05FeHg4/vznP5uOQkTktSwrH+p0zvVlkOvGDTDqKU6zrZ6wsDBMnz4d69atw549e0zHISLySpaO+RCRABH5BsBxAP9U1W0VPGeSiOwQkR0nTpywMg7VQemGcu3btzcdxes98MADCAsL49UPIqJKWFo+VLVYVXsAuBJAjIh0qeA5S1U1SlWjmjdvbmUcqoPU1FS0adMGISEhpqN4vfDwcEybNg0ffvgh9u/fbzoOEZHX8chsF1U9BWAzgFs8cT5yP850qZkHH3wQoaGheOaZZ0xHISLyOpaVDxFpLiJNXfftAG4EkGLV+cg6JSUlLB811KxZM0yZMgXvvvtu2WBdIiJysvLKRySATSKyF0AynGM+1lt4PrLI0aNHkZuby2m2NfTII4+gQYMGePbZZ01HISLyKlbOdtmrqj1VtZuqdlFV7rZVT5Xu6cIrHzVzxRVXYPLkyVixYgXS0tJMxyEi8hpc4ZQuidNsa+8Pf/gDAgMDMXfuXNNRiIi8BssHXVJqair3Kqmlli1bYvz48Xjrrbdw5MgR03GIiLwCywddUulgU24oVzszZ84EAMybN89wEiIi78DyQZfEmS5107p1a9x7771YtmwZHA6H6ThERMaxfFCVcnJy8MMPP7B81NGsWbNQVFSE559/3nQUIiLjWD6oSqVrVHCabd1cc801GDNmDF577TVkZWWZjkNEZBTLB1WJ02zd57HHHkN+fj4WLFhgOgoRkVEsH1SllJQUiAg3lHODDh06IDY2FosXL8bJkydNxyEiMoblg6pUuqGc3W43HcUnPP744zh//jxeeukl01GIiIxh+aAqcaaLe3Xu3Bl33nknXnnlFZw6dcp0HCIiI1g+qFLHjh3Dnj17cOWVV5qO4lOeeOIJnDlzBgsXLjQdhYjICJYPqtSsWbPKdrQl9+nevTuGDRuGF198EWfOnDEdh4jI41g+6CJ2ux0ighUrVgAAtmzZAhHhuA83evLJJ3Hq1CksWrTIdBQiIo9j+aCLpKWlYfTo0QgKCgIAhISEIC4uDunp6YaT+Y5evXph8ODBeOGFF5CTk2M6DhGRR7F80EUiIyMRFhaGoqIiAEBBQQE3lrPAk08+iZMnT2LJkiWmoxAReRTLB1UoKysL1113HZo1a4b4+HhkZmaajuRzfvnLX+LGG2/E888/j9zcXNNxiIg8huWDKrRmzRq0bdsWV199NRYtWoQ1a9aYjuSTnnzySWRlZeH11183HYWIyGNYPqhSDocDLVu2NB3Dpw0YMAADBgzAc889h/79+/MKExH5hUuWDxGJEpGHROR5EZktIneLyGWeCEdmsXx4xpNPPonMzEx89dVXmD17tuk4RESWq7R8iMi9IrILwKMA7ABSARwH0A/AP0XkbRFp7ZmY5GkFBQU4ceIEy4fF7HY7brrpJgCAqiIxMZHTmonI5wVW8b2GAPqqaoUj4USkB4BrARyxIBcZVnr5n+XDWmlpaZgxYwY+/PBD5OfnIygoCHfffTfmz59vOhoRkWUqLR+qWunqRyLyoKq+ZEki8goOhwMAy4fVSqc1FxYWwmazlf3Jac1E5MtqO+D0YbemIK/D8uE5WVlZiI+PxyeffAKbzYZNmzaZjkREZKmqPnapirg1BXkdlg/PKT+N+dFHH8UzzzyDr7/+Gr/85S8NpiIisk5tr3yoW1OQ1zl27BgCAwPRrFkz01H8yqxZsxAZGYnp06ejpKTEdBwiIktUNdvlrIicKfdn6e0sAP467OMcDgciIyNhs3EpGE9q1KgR5s2bh+Tk5LKN/YiIfE2lP1lUtbGqhpX7M6zc17X9uIbqCYfDgVatWpmO4Zfi4uIQExODWbNm4dy5c6bjEBG5XVVXPnaIyMsicouIhHgyFJnHBcbMsdlsePnll5GRkYHnnnvOdBwiIrer6pp6HwAfARgE4AsR2SAi00Wkg0eSkVEsH2b16dMHY8aMwYIFC5Cenm46DhGRW1X1sUuRqm5W1Vmq2hvAeABnAfxZRHaLyGKPpSSPOn/+PE6dOsXyYdjcuXMREBCAP/7xj6ajEBG5VbVHE6pqhqouV9W7AfQCsNK6WGRSRkYGAE6zNa1Vq1aYNWsWPvjgA3zxxRem4xARuU11Npb7RETWlb8BeBtAFMeC+Cau8eE9ZsyYgdatW2P69OkoLi42HYeIyC2qc+UjDcA5AK+7bmcAZAHo4PqafAzLh/ew2+14/vnnsWfPHixfvtx0HCIit6hO+eipqqNV9RPXbQyAGFWdCuB6i/ORASwf3uWuu+5C//798fjjj+P06dOm4xAR1Vl1ykdzEWld+oXrfumylwWWpCKjHA4HQkJC0LRpU9NRCICI4KWXXsLJkycxZ84c03GIiOqsOuXjEQBJIrJJRDYD2ALgDyLSEM6xH+Rjjh07hpYtW0KEW/h4i+uvvx7jxo3DwoULcfDgQdNxiIjq5JLlQ1U3ALgWwIOu23Wq+qmq5qjqS5amIyO4xod3euaZZxASEoJHHnnEdBQiojqpzmyX3wG4B0B3AN0A3O16jHwUl1b3TldccQWeeOIJrF+/Hp9//rnpOEREtVadj12iy936A3gawFALM5FBqsorH15s+vTpaNeuHR566CEUFhaajkNEVCvV+dhlWrnbRAA9AQRbH41MOHv2LHJyclg+vFSDBg2wYMEC/Pe//8WSJUtMxyEiqpXa7Jd+Hs4xIOSDOM3W+w0dOhS/+c1v8NRTT+HHH380HYeIqMZqusLpegCpAD62PhqZwPLh/Uqn3p4+fRpPP/206ThERDVWnSsf8wEscN2eAzBAVWdZmoqMYfmoH7p06YL4+HgkJiZi8+bNGDhwIDIzM03HIiKqlkrLh7gWeVDVL8rdvlLVoxc+h3xHafmIjIw0nIQu5U9/+hMaN26MsWPHIikpCbNnzzYdiYioWqq68rFJRKaVX90UAEQkWER+LSJvA/i9tfHI0xwOBxo3bozGjRubjkKXcNVVV+HUqVM4evQoSkpKkJiYCBGB3W43HY2IqEpVlY9bABQDWC0iDhH5r4ikAfgOwCgAL6rqWx7ISB5Uuropeb+0tDTExsbCZnP+3zgkJARxcXFIT083nIyIqGqBlX1DVfMALAawWESC4NzPJVdVT3koGxnANT7qj8jISDRt2hSqCgDIy8tDUFAQIiIiDCcjIqpatabaqmqhqmawePg+lo/6JSsrCwkJCVixYgVsNhvWr1/PxceIyOtVeuWD/E/p6qZcWr3+WLNmTdl9VcXYsWPx8MMP45VXXjGYioioaiwfVOann35CQUEBr3zUU2PGjMGePXswf/58dO/eHRMmTDAdiYioQtVZZGxedR6j+o9rfNR/c+fOxc0334wpU6YgKSnJdBwiogpVZ8zHTRU8NtjdQcg8lo/6LyAgAKtXr0abNm1w55134siRI6YjERFdpKpFxhJE5FsAHUVkb7lbOoC9notInsLy4RvCw8Oxbt065OXlYfjw4Th//rzpSEREP1PVlY9VAG6Hcx+X28vdeqnqGA9kIw/j6qa+o2PHjli1ahV2796N8ePHl03HJSLyBpWWD1U9raqHAHwPoAOAk6p6WFV/qs4Li8hVIrJJRA6IyH4Rme6eyGQVh8OByy67DCEhIaajkBvcdtttePbZZ/HOO+9g3jwO0yIi71GdMR/pcK5oukNEtovIAhEZVo3jigA8oqqdAPQBMFVEOtchK1mMa3z4npkzZ2LUqFF47LHHsH79etNxiIgAVKN8qOpyVR0H4AYAKwDc5frzUsdlqOou1/2zAA4A4AISXoxLq/seEcEbb7yBnj17YvTo0Thw4IDpSERE1Zpq+4aI/AdAIpzrgowEEF6Tk4hIGwA9AWyr4HuTRGSHiOw4ceJETV6W3IxXPnxTaGgo1q5dC7vdjmHDhiE7O9t0JCLyc9X52OVyAAEATgH4Cc6xH0XVPYGINALwIYAHVfXMhd9X1aWqGqWqUc2bN6/uy5KbFRcXIzMzk6ub+qirrroKa9aswaFDhzBq1CgUFxebjkREfqw6H7sMV9XeAP4CoCmATSJytDov7tqQ7kMAK1V1zaWeT+acOHECxcXFvPLhw/r27YvFixfjs88+w8yZM03HISI/dsnl1UVkCID+AAbA+XHLRgBbqnGcAFgG4ICqvlDHnGQxrvHhHyZMmIA9e/ZgwYIF6N69O8aOHWs6EhH5oep87DIYwC4Ad6pqR1W9T1WXV+O4vgDGAvi1iHzjut1al7BkHZYP//HCCy/ghhtuwMSJE7FhwwYMHDgQmZmZtXqtjIyMOh1PRP6pOh+7TFXVd1XVAQAiEiAicdU4LklVRVW7qWoP122DO0KT+7F8+I+goCC89957iIyMxF133YUtW7Zg9uzZtXqtOXPmICkpqdbHE5F/kspWPhSRMABT4Zweuw7APwHcD2AGgG9UtTprfdRIVFSU7tixw90vS9Xw9NNPY/bs2cjPz0dQUJDpOGQxu92OvLy8Cr/XqFGjSx5/7ty5Ch8PCQlBbm5unbIRke8QkZ2qGnXh41WN+fgbgGwAXwOYAOAPAIIBDFPVb6wISeY4HA60aNGCxcNPpKWlYcaMGfjwww+Rn5+PwMBAtG/fHgMHDkTDhg0vefy5c+ewefNmHDx4EIBzOu/w4cMxf/58q6MTkQ+oqnxco6pdAedaHwBOAmjtWjCMfAzX+PAvkZGRCAsLQ2FhIUJCQlBQUIAbbrgBixcvrvZrJCQk4ODBgwgICEBeXh7CwsIQERFhYWoi8hVVlY/C0juqWiwi6Swevovlw/9kZWUhPj4ekyZNwtKlS5GRkVHj49u0aQNVxW233Vbj44nIf1U15qMYQE7plwDsAM677quqhrk7DMd8mHPFFVdg2LBhWLp0qekoVI8888wzeOKJJ5CdnY2mTZuajkNEXqayMR9V7WoboKphrltjVQ0sd9/txYPMKSwsxPHjx3nlg2osJiYGAMBfGoioJqqzzgf5uNI1Gri0OtVUVJTzF5rk5GTDSYioPmH5IK7xQbUWHh6Oa6+9Ftu3bzcdhYjqEZYPYvmgOomOjuaVDyKqEZYPYvmgOomJicGxY8fK/h0REV0KywfB4XAgICAAzZs3Nx2F6qHSQae8+kFE1cXyQXA4HIiMjITNxn8OVHM9evRAYGAgx30QUbXxpw1xgTGqE7vdjq5du7J8EFG1sXwQywfVWXR0NHbs2IGSkhLTUYioHmD5IBw7dozlg+okJiYGp06dwvfff286ChHVAywffi43NxfZ2dksH1Qn0dHRADjolIiqh+XDz5VuBsbyQXXRuXNnhIaGctwHEVULy4efK12bgUurU10EBgaiV69evPJBRNXC8uHnuMAYuUt0dDR2796NwsJC01GIyMuxfPg5lg9yl5iYGOTl5WHfvn2moxCRl2P58HMOhwMNGjRAeHi46ShUz5UOOuW4DyK6FJYPP1e6xoeImI5C9Vzbtm1x+eWXs3wQ0SWxfPg5LjBG7iIi3OGWiKqF5cPPsXyQO8XExGD//v3IyckxHYWIvBjLh59j+SB3io6ORklJCXbt2mU6ChF5MZYPP3b27FmcPXuW5YPchoNOiag6WD78GKfZkrtdccUVuPrqqznug4iqxPLhx1g+yArR0dG88kFEVWL58GNcWp2sEBMTg/T0dJw4ccJ0FCLyUiwffoxXPsgKpeM+duzYYTgJEXkrlg8/5nA40KhRIzRu3Nh0FPIhvXr1gojwoxciqhTLhx/jNFuyQuPGjdGpUycOOiWiSrF8+DGWD7JKTEwMtm/fDlU1HYWIvBDLhx9j+SCrREdH48SJEzhy5IjpKETkhVg+/JSqsnyQZWJiYgBwsTEiqhjLh5/Kzs5GXl4eywdZolu3bggODua4DyKqEMuHn+I0W7JScHAwevTowSsfRFQhlg8/xfJBVouJicGOHTtQXFxsOgoReRmWDz/F8kFWi46ORk5ODlJSUkxHISIvw/Lhp1g+yGocdEpElWH58FMOhwPh4eGw2+2mo5CP6tChA8LCwjjolIguwvLhpzjNlqxms9kQFRXFKx9EdBGWDz/F8kGeEBMTg7179yIvL890FCLyIiwfforlgzwhOjoahYWF2LNnj+koRORFWD78UElJCTIyMlg+yHIcdEpEFWH58EMnT55EUVERywdZrlWrVoiIiOCgUyL6GZYPP3Ts2DEAnGZL1hORsh1uiYhKsXz4Ia7xQZ4UHR2N1NRUnD592nQUIvISLB9+iOWDPKl03MeOHTsMJyEib8Hy4YdKy0dERIThJOQPoqKiAIDjPoioDMuHH3I4HGjRogWCg4NNRyE/cNlll6F9+/Yc90FEZVg+/BDX+CBPi4mJ4ZUPIirD8uGHWD7I06Kjo3H06NGyj/yIyL+xfPghlg/ytNJBp7z6QUSAheVDRJaLyHER2WfVOajmioqKkJWVxfJBHtWjRw8EBASwfBARAGuvfLwF4BYLX59qISsrC6rK8kEeFRoaiq5du3LQKREBsLB8qOqXAH6y6vWpdrjGB5kSHR2N5ORkqKrpKERkGMd8+BkurU6mxMTE4NSpU/j+++9NRyEiw4yXDxGZJCI7RGTHiRMnTMfxebzyQaZw0CkRlTJePlR1qapGqWpU8+bNTcfxeQ6HAzabDS1atDAdhfxM586dYbfbOe6DiMyXD/Ish8OBiIgIBAQEmI5CfiYwMBC9evXilQ8isnSq7WoAXwO4TkSOish4q85F1edwONCqVSvTMchPRUdHY9euXSgsLDQdhYgMsnK2yyhVjVTVIFW9UlWXWXUuqj4uMEYmxcTEIC8vD/v2cfkfIn/Gj138DMsHmRQdHQ2Ag06J/B3Lhx/Jz8/Hjz/+yPJBxlxzzTW47LLLOOiUyM+xfPiRjIwMAJxmS+aISNliY0Tkv1g+/AjX+CBvEBMTg3379iEnJ8d0FCIyhOXDj3B1U/IGMTExKCkpQd++fZGZmWk6DhEZwPLhR3jlg7xB6aDTvXv3Yvbs2YbTEJEJLB9+xOFwICgoCJdffrnpKOSn7HY7IiIiAACqisTERIgI7Ha74WRE5EksH36kdJqtiJiOQn4qLS0No0ePRmBgYNlj7dq1w9q1a82FIiKPY/nwI1zjg0yLjIxEWFgYSkpKEBwcDAA4dOgQbrnlFsTExGD58uU4f/684ZREZDWWDz/CpdXJG2RlZSE+Ph7bt2/HlClTMHjwYLzyyivIycnB+PHj0apVKzz44IM4cOCA6ahEZBFRVdMZykRFRemOHTtMx/BZTZo0wb333ouXX37ZdBSii6gqtmzZgiVLluCDDz5AYWEhBg0ahISEBNxxxx1lV0qIqP4QkZ2qGnXh47zy4SfOnTuHM2fO8GMX8loiggEDBmDVqlU4evQo5s6di8OHD+Oee+5B69at8fjjj+Pw4cPIyMjAwIEDOU2XqB5j+fATXN2U6pMWLVpg5syZ+P7777Fhwwb07t0bc+fORdu2bdGnTx9s2bIFTz/9tOmYRFRLLB9+gmt8UH1ks9kwePBgfPzxxwgKCoKq4siRI1BVvPbaaxARhISEmI5JRDXE8uEnWD6ovktPT8fo0aPL1gSx2Zz/+SouLkZsbCy++OILeNMYNiKqHMuHn+DS6lTflU7Tzc/PL7vaERsbi2nTpuHzzz/HoEGD8Itf/AILFy7EqVOnzIYloiqxfPgJh8OB0NBQhIWFmY5CVGul03S3bt2K+Ph45Ofn44UXXsCxY8fw5ptvIiwsDNOnT0fLli0xfvx4cPYckXfiVFs/ERsbi507d+K7774zHYXIUrt370ZiYiJWrlyJ8+fPIyoqCvHx8YiNjUXDhg2RkZGB2NhYvPvuu2VLvRORNTjV1s9xdVPyFz179sTSpUvhcDjw6quvIjc3FxMmTECrVq3wwAMP4KGHHkJSUhI3tSMyiFc+/ET79u0RHR2N1atXm45C5FGqiq+++goDBw5ESUnJRd8PCQlBbm6ugWREvo9XPvyYqnJpdfJbIoJ+/frh6NGjGDFiRNmmdjabDSNHjkR6errhhET+h+XDD5w+fRq5ubn82IX8WmRkJFq0aIGSkhIEBQWhpKQE//znP5GTk2M6GpHfYfnwA1zjg8ipdLZMcnIy7rjjDpw/fx59+vTB1q1bTUcj8iuBpgOQ9Vg+iJzWrFlTdv+jjz7CwYMHMXjwYNxwww1YuXIlRowYYTAdkf/glQ8/wPJBVLEOHTpg69at6NGjB0aOHIkXX3yRq6QSeQDLhx8oXd00MjLScBIi79O8eXNs3LgRI0aMwMMPP4wHHngAxcXFpmMR+TSWDz/gcDjQpEkTNGzY0HQUIq9kt9vx3nvvYcaMGXj11VcxfPhwDkQlshDLhx/gAmNEl2az2fD8889j0aJF+PTTTzFw4EBkZmaajkXkk1g+/ADLB1H1TZkyBR9//DEOHDiAPn36YP/+/aYjEfkclg8/wPJBVDNDhgzBl19+ifz8fPTt2xcbN25ERkZGna6G1PV4Il/C8uHjSkpKkJGRwfJBVEO9evXC1q1b0apVK9xyyy0YM2ZMnfaEmTNnDveUIXLh3i4+7sSJE2jRogUWLlyIadOmmY5DVO+EhIQgPz//oseDgoLwzjvvXPL42NhYFBYWVvi63FOGfF1le7twkTEfxzU+iOomPT0dDz/8MN5///2fTcEtLCzEnXfeWePXa9CgAUaOHIn58+e7MyZRvcLy4eNYPojqJjIyEk2bNoWqokGDBigoKMDIkSPxxBNPVPs1/vznP+ODDz4AAOTn58PhcCAiIsKqyERej+XDx7F8ENVd6Z4wkyZNwtKlS5GRkYFu3bpV+/iioiIkJCTg7rvvRlxcHDZt2oRHHnkE8+bNK9tll8if8F+9jystH/wti6j2yu8Js2jRojodf+jQITzyyCN44YUXsG/fPrzzzjsIDw93S06i+oKzXXzcwYMHERgYiOzsbNNRiAhAYGAgXn75ZbzxxhvYtGkTYmJicODAAdOxiDyK5cPHbdmyBUVFRZzeR+Rlxo8fj02bNuHMmTPo3bs3Pv30U9ORiDyG5cNH2e12iAgOHz4MAEhMTISIwG63G05GRKX69u2L5ORktG/fHrfffjvmzZvHXXXJL7B8+Ki0tDRERf1vanVoaCji4uKQnp5uMBURXah169ZISkrC3XffjVmzZmHMmDFc/4N8HsuHj7LZbNi7dy8A52JGeXl5CAsL48BTIi8UGhqK1atX49lnn8Xq1avRv39/HD161HQsIsuwfPioGTNmoLCwEKNGjcLWrVsRHx/PPSWIvJiI4NFHH8XHH3+MgwcPIioqCl9//bXpWESWYPnwQZs3b8aKFSvw6KOPYtWqVejevTsWLVr0s+l+ROSdbr/9dmzduhWNGjXCoEGD8Oabb9ZpUzqTG9pxMz2qDMuHjykoKMCUKVPQpk0bPP7446bjEFEtdO7cGdu3b8eAAQMwbtw43HzzzbXelM7khnbcTI8qw43lfMzcuXPx6KOP4pNPPsGQIUNMxyGiOrDb7cjLy7vocZvNhsGDB1d57N///neUlJRc9LgnNrSrLDc30/M/lW0sxysfPuTQoUOYPXs27rjjDhYPIh+QlpaG0aNHIzg4GIBzXEh4eDi6dOmCzMzMKm9dunRBeHg4RKTs9Vq3bo2VK1daNp03KysLzz33HJo1a3bR9zp16oQvv/zSkvNSPaSqXnPr1auXUu0NHTpUQ0ND9fDhw6ajEJGbxMfHq81m05CQELXZbJqQkFDjYxs0aKAANDg4WAFox44d9aWXXtLs7Ow65yspKdHNmzfrPffco0FBQQpAb7jhBr3pppvUZrOVnROA2mw2vfPOO3XLli1aUlJS53OT9wOwQyv4ec8rHz5i3bp1WLduHZ566im0bt3adBwicpPSTe1qM2ut9Nht27ZhypQpuOWWW/DWW2+hadOmePDBB9GyZUuMGzcOycnJNc516tQpLFy4EL/4xS8waNAgfPbZZ5g6dSoOHDiAjRs3olGjRoiPj8f27dsxZcoU3HzzzZgxYwY2btyI/v37IyoqCn/961+Rn59f43OTD6iokZi6WXHlY/Xq1QpA33//fbe/trc4d+6cXn311dq5c2ctKCgwHYeI6oFdu3bppEmTtGHDhgpAe/Xqpa+//rqeO3eu7DkOh0MHDBigGRkZZY8lJyfr+PHjNTQ0VAFoTEyMvvnmm5qTk1Ot8547d06XLFminTp1UgB6xRVX6FNPPfWzc1R0XqqfUMmVD+OFo/zNivJReskvODjY7a/tLWbNmqUA9IsvvjAdhYjqmdOnT+uiRYu0S5cuCkCbNGmi06ZN0/3792tCQoLabDadMGGCvvHGGxoVFaUANDQ0VCdOnKg7d+6s9XlLSkr0s88+01tvvbXsv9G/+93vdOfOnWXnrclHTOSdKisfPjvbpfwgqwt509+5rv773/+ie/fuiIuLw1tvvWU6DhHVU6qKr776CkuWLMHKlSsrfI6IYOHChRg7diyaNGnitnMfPHgQr7zyCl599dUKv89ZMvWX3812Wb16NQICAn72WEBAAN5//31DidxPVTF16lQ0btwYf/nLX0zHIaJ6TETQr18/rFixAt9++y169OhR9ktcQEAAbrrpJhw7dgz333+/W4sHAHTo0AGvvPIKDhw4gJ49e/7sl8cuXbogKSnJrecj83y2fMTGxlZYPkaOHGkokfutXLkSmzdvxnPPPYcWLVqYjkNEPqJLly7o06cPRAQhISFQVbRv3x6RkZGWnrdjx47o3bs3RARBQUEAgH379iEmJgbDhg3DP/7xjwrXLqH6x2fLBwAUFRWVNejAwEAUFRUZTuQ+2dnZeOSRRxATE4OJEyeajkNEPqYus2zccd7k5OSyWTIzZ87E1q1bMXjwYLRv3x7z5s3D8ePHPZKHrGHpmA8RuQXAywACALyhqnOrer4VK5wOHToUn3zyCT7++GMMHTrUra9t0tSpU7FkyRIkJyfj+uuvNx2HiMhSBQUF+Oijj5CYmIgvvvgCQUFBGDlyJBISEtCvX78qx/mROR4f8yEiAQAWARgMoDOAUSLS2arzVSY9PR0A8I9//MNnNjlKTk5GYmIipk6dyuJBRH4hODgY99xzDzZv3oz9+/cjISEBGzZswIABA9C1a1e8+uqrOH36NIC6bWjHYz3DsisfIvJLAE+r6s2urx8FAFV9rrJjrLjyUdqGRQSTJ0/G0qVLMXnyZCxevNit5/GU4uJi9O7dG8eOHUNKSorbB34REdUXOTk5eOedd5CYmIidO3ciNDQUo0ePRnZ2Nj766CP87ne/w7x582r0mjNnzsRf//pXvzvWqp+LlV35sLJ8jARwi6pOcH09FkBvVb2/smM8NdW2VH2cvrVo0SLcf//9WLVqFUaNGmU6DhGRV9ixYwf69OmD4uJi01HqNXf/XKysfAS67QwVnLOCxy5qOiIyCcAkAG5dFnzx4sWYMmVKlc/p2rUrVqxYgZEjRyIkJMRt57ZKZmYmHnvsMfzmN79BbGys6ThERF4jKioKP/zwAx544AGsW7cOBQUFCAoKQo8ePTB8+PBLXiU+ffo01qxZgz179qCwsNDvjg0NDcXw4cMxf/78Ko9zFyvLx1EAV5X7+koAjgufpKpLASwFnFc+3HXyhIQE3H///RdNy2rQoAEKCwvRp08fnDhxAmPHjsWDDz6IcePGYdKkSWjfvr27IrjdjBkzkJeXh0WLFnFwFRHRBSIjI9GsWTMUFRUhJCQEBQUFiIqKwqOPPlqt448cOYJdu3b55bF5eXkICwtDREREtY6tKyun2iYDuFZE2opIMIBYAOssPN9FKpoPvm3bNsTHx+OKK65ASkoK/vWvf2HQoEF44YUXcO211+Lmm2/GRx995HXTcjdu3IiVK1fij3/8I6677jrTcYiIvJI7NuLjsdazeqrtrQBegnOq7XJVfaaq51s54BQANmzYgMGDB1f4PIfDgWXLlmHp0qU4evQoWrVqhYkTJ2LChAlo1aoVAOeo4NjYWLz77rs1bod1Ofbw4cPo1KkTmjVrhtTUVNjt9hodT0REZIKR5dVVdYOqdlDVdpcqHp5w8803V/q9li1b4sknn0R6ejrWrl2Lrl274k9/+hOuvvpqjBgxAp9//jlmz56NpKQkzJ49u8bnnjNnTq2Pvfvuu5Gbm4uuXbuyeBARUb3HjeWqkJaWhtdee63SfVNEBL17967yNbZt21bh+epybH2cpUNERP7H7zaWGzduXIWPT5o0qdqvcc0112DevHlIT0/Hr371q7K9Ymw2GyIiIjBgwACEhYVVeRswYAAiIiJgs9lqfWzpeUNDQxEXF1e2cBoREVF9ZOVsF6OWLVuG5cuXX/T4a6+9VuPXatOmDbp164atW7eWjSgePnx4tRdkSUhIwNKlS+t8rKdHIxMREVnBZ698lOeO2SH+OBqZiIjICj475oOIiIjM8rsxH0REROSdWD6IiIjIo1g+iIiIyKNYPoiIiMijWD6IiIjIo1g+iIiIyKNYPoiIiMijWD6IiIjIo1g+iIiIyKNYPoiIiMijWD6IiIjIo7xqbxcROQHgsAUv3QzASQtely7G99pz+F57Dt9rz+D77Dmeeq+vVtXmFz7oVeXDKiKyo6KNbcj9+F57Dt9rz+F77Rl8nz3H9HvNj12IiIjIo1g+iIiIyKP8pXwsNR3Aj/C99hy+157D99oz+D57jtH32i/GfBAREZH38JcrH0REROQlfL58iMgtIpIqIt+LyCzTeXyViCwXkeMiss90Fl8mIleJyCYROSAi+0VkuulMvkpEQkRku4jscb3XfzKdydeJSICI7BaR9aaz+DIROSQi34rINyKyw0gGX/7YRUQCABwEcBOAowCSAYxS1f8aDeaDRGQAgHMA/qqqXUzn8VUiEgkgUlV3iUhjADsB3MF/0+4nIgKgoaqeE5EgAEkApqvqVsPRfJaIPAwgCkCYqg4xncdXicghAFGqamxNFV+/8hED4HtVTVPVAgDvABhmOJNPUtUvAfxkOoevU9UMVd3lun8WwAEArcym8k3qdM71ZZDr5ru/rRkmIlcCuA3AG6azkPV8vXy0AvBDua+Pgv+hJh8hIm0A9ASwzXAUn+X6GOAbAMcB/FNV+V5b5yUAfwRQYjiHP1AAn4vIThGZZCKAr5cPqeAx/uZC9Z6INALwIYAHVfWM6Ty+SlWLVbUHgCsBxIgIP1K0gIgMAXBcVXeazuIn+qrq9QAGA5jq+tjco3y9fBwFcFW5r68E4DCUhcgtXOMPPgSwUlXXmM7jD1T1FIDNAG4xm8Rn9QUw1DUW4R0AvxaRFWYj+S5Vdbj+PA7gIziHKHiUr5ePZADXikhbEQkGEAtgneFMRLXmGgS5DMABVX3BdB5fJiLNRaSp674dwI0AUoyG8lGq+qiqXqmqbeD87/RGVR1jOJZPEpGGrsHqEJGGAH4LwOOzFH26fKhqEYD7AXwG58C891R1v9lUvklEVgP4GsB1InJURMabzuSj+gIYC+dvht+4breaDuWjIgFsEpG9cP4i809V5RRQqu+uAJAkInsAbAfwqar+w9MhfHqqLREREXkfn77yQURERN6H5YOIiIg8iuWDiIiIPIrlg4iIiDyK5YOIiIg8iuWDyI1EpKmITCn3dUsR+cDDGVaLyF4ReUhEHhSR0Dq+XryI/O4Sz4kSkYW1fP2fZRSRDaXra9SV67WrzF7BMf+pw/kGicivanlscxHx+JRHIhM41ZbIjVz7raw3tbOviEQA2KaqV7u+PoQa7l4pIgGqWmxRxIrOdwgW7LApIoEAdgG43rXmj+VE5GkA51R1fi2PfxPAG6r6lVuDEXkZXvkgcq+5ANq5Fv96XkTaiMg+ABCRe0VkrYh8IiLpInK/iDwsIrtFZKuIXOZ6XjsR+Ydr06ctItLxwpOISIyI/Md17H9E5DrXtz4H0MJ1/qcAtIRzoaxNruN+KyJfi8guEXnftUcMROSQiPw/EUkCcNcF53paRGa47m8WkXkisl1EDopIf9fjg0Rkvet+QxFZLiLJrnzDXI8HiMh8EfnWdWVmmog8UEHGQyLSzHX/YRHZ57o96HqsjYgcEJHXRWS/iHzuWoH0Qr8GsKu0eLiyvygiX7qOjxaRNSLynYj8udzf91y5v9NmEflARFJEZKWISAUZo1zPawMgHsBDrve/v+tqxoeu9yJZRPq6jhko/1skbre4VpwEsBZAXBX/voh8g6ryxhtvbroBaANgX0VfA7gXwPcAGgNoDuA0gHjX916Ec5M4APg3gGtd93vDudT0hecJAxDoun8jgA8rOf8hAM1c95sB+BJAQ9fXMwH8v3LP+2Mlf6enAcxw3d8MYIHr/q0A/uW6PwjOKz4A8CyAMa77TQEcBNAQQAKce9KU5r7swozlvwbQC8C3rmMbAdgP5y6+bQAUAejhev57pee7IPefAEwr9/VmAPNc96fDuc9TJIAGcO4Ddbnre+fK/Z1Ow7knlA3OFXz7VfC+RgHYfOF75fp6VbljWsO5LD4AfALn5l5w/d1K35NWAL41/e+YN96svgWCiDxpk6qeBXBWRE7D+UMIcP6Q7ea6EvErAO+7fskGnD8cL9QEwNsici2cOzUHVePcfQB0BvCV67WD4fyBWurdav4dSjez2wlnEbjQb+HcJGyG6+sQOH/w3ghgibquRKjqT5c4Tz8AH6lqDgCIyBoA/eHcnyldVb+5RI5IOLdVKK90b6dvAexX1QzXa6fBuQnljxc8f7uqHnU95xvXeZIukbu8GwF0Lve/ZZjrKsdXAF4QkZUA1pSeA8BxOK8EEfk0lg8iz8ovd7+k3NclcP7/0QbglDq3ca/KHDiLzHDX5f7N1Ti3wLk/yahKvp9TjdcA/pe5GBX/N0QA3KmqqT970PkTuCaDzKSK75V/H4sBVPSxSy6cxaei48q/96VfV/R3ufA8pc8pwv8+tr7wHOXZAPxSVXMveHyuiHwK59WjrSJyo6qmuF7rwucS+RyO+SByr7NwfqxSK6p6BkC6iNwFOH9gi0j3Cp7aBMAx1/17q5lnK4C+ItLe9dqhItKhtlmr8BmAaeXGR/R0Pf45gHhxDgRF6RgXVP6efQngDlfOhgCGA9hSgxwHALSvRf7qOATnx0IAcGe5xy/8u3wO5+aWAAAR6eH6s52qfquq8wDsAFA6rqcDDOwwSuRpLB9EbqSqP8L5scY+EXm+li8TB2C8OHed3A9gWAXP+QuA50TkKwABVbzWUgB/F5FNqnoCzqKyWpw7tW7F/37ouUPpVY05cH4MtFecg23nuB5/A8AR1+N7AIy+MOPPXkx1F4C34Nx5cxucs0B21yDP3wEMqMXfozr+BOBlEdkC5xWRUp8AGF464BTAAwCiXANs/wvngFQAeND1b2QPnFc6/u56/AYAn1qUmchrcKotEdWZiNwJYKiq/t50lvJE5CM4B9J+ZzpLdYjIlwCGqWq26SxEVuKVDyKqExEZCuAZAK+ZzlKBWXAOPPV6ItIcwAssHuQPeOWDiIiIPIpXPoiIiMijWD6IiIjIo1g+iIiIyKNYPoiIiMijWD6IiIjIo1g+iIiIyKP+P8X5fq1ilpIGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create the figure and label the axes\n",
    "plt.figure(figsize=(9,6))\n",
    "\n",
    "plt.xlabel('time after injection (minutes)')\n",
    "plt.ylabel('Rwt (ug/L)')\n",
    "plt.plot(df.min_after_inj,df.Rwt_in_stream_ug_per_L,'k-*')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above graph shows us the timeseries of measurements made downstream.  If we presume that the dischage is constant over the time period of measurement (which is generally reasonable for about 5 minutes), and that all of the tracer passes our measurement location, then the total area under the tracer curve above times the discharge that has diluted it should give us the orignal amount we put in:\n",
    "$$ V = Q\\int{RC(t)dt} $$\n",
    "where V is the total mass of tracer injected, Q is the discharge, and RC is the concentration of tracer measured at each time step t."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "348.9\n",
      "62883.34766408714\n"
     ]
    }
   ],
   "source": [
    "# the measured amount of tracer we injected, in micrograms of active ingredient\n",
    "# here we weighted 109.7 grams of a 20% solution.  The sensor measures micrograms per liter\n",
    "V = 109.7*0.2*1000000 \n",
    "intRC = np.sum(df.Rwt_in_stream_ug_per_L*df.dt_sec)\n",
    "print(intRC)\n",
    "Q = V/intRC\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above gives us discharge in liters per second"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "62.88334766408714\n"
     ]
    }
   ],
   "source": [
    "# We can convert this to cubic meters per second\n",
    "Qcms = Q/1000\n",
    "print(Qcms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "The above gives our estimate of discharge in cms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2220.706557752938\n"
     ]
    }
   ],
   "source": [
    "# We can also convert this to cubic feet per second\n",
    "Qcfs = Q*0.0353147\n",
    "print(Qcfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above is discharge is cfs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}