Muskingum_routing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we solve a storage flow relationship, assuming the so-called Muskingum equation. An important law we use (as always) is the continuity equation (i.e. how much flow goes into and out of a given river section, how much does this change the storage?)\n",
    "\n",
    "$$\n",
    "\\frac{\\partial{A}}{\\partial{t}} + \\frac{\\partial{Q}}{\\partial{x}} = 0\n",
    "$$\n",
    "\n",
    "Here $A$ [L$^2$] is the wetted cross-section, $t$ is the time [T], $Q$ [L$^3$ T$^{-1}$] is the flow and $x$ is the longitudinal position [L], usually we use meters and seconds for the magnitudes L and T.\n",
    "\n",
    "We can integrate this equation over a certain length $x$ and simply look at the storage over a given channel section $S$, assuming that at any moment in time $S=Ax$, yielding:\n",
    "\n",
    "$$\n",
    "\\frac{\\Delta A\\ x}{\\Delta{t}} + \\Delta{Q} = 0\n",
    "$$\n",
    "or in simpler terms:\n",
    "$$\n",
    "\\frac{\\Delta S}{\\Delta t} = -\\Delta Q\n",
    "$$\n",
    "When integrating over a certain time span $t_1$ until $t_2$ we yield:\n",
    "$$\n",
    "\\frac{S_2-S_1}{\\Delta t}=0.5\\left(I_1+I_2\\right)-0.5\\left(O_1+O_2\\right)\n",
    "$$\n",
    "where $I$ and $O$ represent $Q$ at the inflow and outflow boundaries of the channel stretch $x$ respectively. \n",
    "\n",
    "Ok, so far we have only simplified the continuity equation. Now we need a storage-flow relationship. Muskingum assumes that the storage can be described with a linear relationship of the  weighted average of inflows and outflows in a given channel section. ADD EXPLANATION, WHAT DOES THIS MEAN IN PRACTICE?\n",
    "\n",
    "$$\n",
    "S=k\\left[fI+\\left(1-f\\right)O\\right]\n",
    "$$\n",
    "It is pretty easy to see that this relationship approaches a linear reservoir when $f$ approaches 0. We can also subsitute this in the earlier found discretized continuity equation:\n",
    "$$\n",
    "\\frac{k}{\\Delta t}\\left[f\\left(I_2-I_1\\right)+\\left(1-f\\right)\\left(O_2-O_1\\right)\\right]=0.5\\left(I_1+I_2\\right)-0.5\\left(O_1+O_2\\right)\n",
    "$$\n",
    "This allows us to relate the outflow from a storage section at a future time step to the inflows over the time step and outflow over the previous time step. Reordering then provides:\n",
    "$$\n",
    "S_2 = c_0I_2+c_1I_1 + c_2O_1\n",
    "$$\n",
    "with:\n",
    "$$\n",
    "c_0 = \\frac{-2f+\\Delta t/k}{2\\left(1-f\\right)+\\Delta t/k}, \n",
    "c_1 = \\frac{2f+\\Delta t/k}{2\\left(1-f\\right)+\\Delta t/k}, \n",
    "c_2 = \\frac{2\\left(1-f\\right)-\\Delta t/k}{2\\left(1-f\\right)+\\Delta t/k}\n",
    "$$\n",
    "Ok, this is extremely simple! Let's put this in a little model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class simple_runoff(object):\n",
    "    \"\"\"\n",
    "    A very very simple rainfall runoff response model, computing specific runoff as follows:\n",
    "    q = max(P-I, 0)*rc\n",
    "    where P is precipitation in mm per day, I is interception in mm/day, rc is runoff coefficient [-]\n",
    "    Consequently, it will compute flows at the outlet as follows:\n",
    "    Q(t) = q(t)A * f(t)\n",
    "    It will convert flows to m3/s (instead of day) before returning\n",
    "    where f(t) is a unit hydrograph convolution function, computed assuming a triangular shape,\n",
    "        with a time of concentration given in days\n",
    "    Input:\n",
    "        I: interception threshold [mm/day]\n",
    "        rc: runoff coefficient [-]\n",
    "        A: catchment surface area [km2]\n",
    "        t_c: time of concentration [day]\n",
    "    \"\"\"\n",
    "    def __init__(self, I, rc, A, t_c):\n",
    "        self.I = I\n",
    "        self.rc = rc\n",
    "        self.A = A\n",
    "        self.t_c = t_c\n",
    "        self.triangle_UH()\n",
    "        \n",
    "    def triangle_UH(self):\n",
    "        \"\"\"\n",
    "        Make a unit hydrograph assuming a time of concentration in time steps\n",
    "        \"\"\"\n",
    "        uh_right = np.append(np.linspace(0, 1, self.t_c), np.linspace(1, 0, self.t_c)[1:])\n",
    "        uh_left = np.zeros(len(uh_right))\n",
    "        self.conv = np.append(uh_left, uh_right)/uh_right.sum()\n",
    "\n",
    "    def compute_flows(self, P):\n",
    "        q = np.maximum(P-self.I, 0)*self.rc\n",
    "        _Q = np.convolve(q*A*1e3/86400, self.conv, 'same')\n",
    "        return pd.Series(_Q, index=pd.DatetimeIndex(start=P.index[0], periods=len(_Q), freq=P.index.freqstr))\n",
    "        \n",
    "class Muskingum(object):\n",
    "    def __init__(self, k, f, dt):\n",
    "        self.k = k  # residence time\n",
    "        self.f = f  # weighting factor\n",
    "        self.dt = dt\n",
    "        self.c0 = (-2*f+dt/k)/(2*(1-f)+dt/k)\n",
    "        self.c1 = (2*f+dt/k)/(2*(1-f)+dt/k)\n",
    "        self.c2 = (2*(1-f)-dt/k)/(2*(1-f)+dt/k)\n",
    "    \n",
    "    def compute_O_2(self, I_1, I_2, O_1):\n",
    "        # check for NaNs in O_1\n",
    "        if np.isnan(O_1):\n",
    "            O_1 = 0\n",
    "        return self.c0*I_2+self.c1*I_1+self.c2*O_1\n",
    "\n",
    "    def compute_outflows(self, s_inflows):\n",
    "        s_outflows = pd.Series(index=s_inflows.index, )\n",
    "        for _i in range(1, len(s_inflows)):\n",
    "            s_outflows.iloc[_i] = compute_O_2(s_inflows.iloc[_i-1], s_inflows.iloc[_i], s_outflows.iloc[_i-1])\n",
    "        return s_outflows\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Timestamp('2001-01-01 00:00:00', freq='D')"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "s = pd.Series(index=pd.DatetimeIndex(start='2001-01-01', end='2001-01-30', freq='D'))\n",
    "s2 = pd.Series()\n",
    "P.index[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare some flows from a unit hydrograph principe\n",
    "p = np.array([0., 5., 20., 18., 15., 12., 8., 10., 2., 0., 0.])\n",
    "P = pd.Series(p, index=pd.DatetimeIndex(start='2001-01-01', periods=len(p), freq='D'))  # mm/day\n",
    "\n",
    "# make a model\n",
    "I = 3. # interception threshold [mm/day]\n",
    "rc = 0.5# event runoff coefficient [-]\n",
    "A = 150000 # surface of catchment [km2]\n",
    "t_c = 5\n",
    "runoff_model = simple_runoff(I, rc, A, t_c)\n",
    "Q = runoff_model.compute_flows(P)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0xfcacd30>"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAETCAYAAAAPucBzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xl8VOX1+PHPSSYJIZBAQLYEZHUB\nZJPNXbDuCy6ooLVoVWpdWmsXbe1PW1vrVrfW1n5t1VKLqFVQVFwQ91Z22RdZxBB2CHsgkOT8/rh3\ncIwzyex3JnPer9e8ZubOM889gcycPM99FlFVjDHGmETJ8joAY4wxjZslGmOMMQllicYYY0xCWaIx\nxhiTUJZojDHGJJQlGmOMMQllicYYY0xCWaIxxhiTUJZojDHGJJTP6wDqk5WVpfn5+V6HYYwxaaWy\nslJVNWUaEimdaPLz89m7d6/XYRhjTFoRkX1exxAoZTKeMcaYxskSjTHGmISyRGOMMSahLNEYY4xJ\nKEs0xhhjEsoSjTHGZBgR6SgiH4jIUhFZLCI/do8Xi8hUEVnh3rcM8f4xbpkVIjKmwfOl8g6bBQUF\nasObjTEmMiJSqaoF9bzeHmivqnNFpDkwB7gQuBqoUNX7ReQOoKWq3l7nvcXAbGAgoO57j1XV7aHO\nl9LzaFKVqnLPG0s4o2c7juvWyutwjInO9L/Bmk8grxCaFAbcNw94XOQ897+eWwAiXkduYqSqG4AN\n7uPdIrIUKAFGAKe6xcYBHwK313n7mcBUVa0AEJGpwFnAhFDns0QThU27qnj2v2t4fkYZT48ZxIk9\nWnsdkjGRmf8ivH07FHUCrYWqXc6tIZJdJxEFJKemreH4W6CoJPHxm4b4RGR2wPOnVPWpYAVFpDPQ\nH5gBtHWTEKq6QUTaBHlLCbA24Hm5eyx0MOHHbfzKKioByPNlce24WTxz9SBO6G7JxqSJDQvg9R/D\n4SfA916D7BzneG0tHNgNVbthv5t4/PeHHu+uc3w37FrvPN61HtbNgWvegmz7avFYtaoObKiQiDQD\nXgFuVdVdEl5rNViheq/B2G9DFPyJ5tlrBnHnpEVOshkziOMt2ZhUV1kBL34X8lvApf/8OskAZGVB\nkyLnVhRF3Qtfhleuhf8+Cif/PF4RmwQRkRycJDNeVSe6hzeJSHu3NdMe2BzkreV83b0GUIrTxRaS\njTqLQllFJSJwTEkLxl83hMOLC/j+uFn8b+VWr0MzJrTaGph4vdPyuOw5aBasVyQGx4yE3pfAh/fD\n+nnxrdvElThNl6eBpar6SMBLkwH/KLIxwGtB3v4OcIaItHRHpZ3hHgvJEk0U1lZU0qEon1xfFq2a\n5TH++iF0Km7qJJtVlmxMivrwflj5HpzzIHQclJhznPNHKGgDE8fCwZRa19F80wnAVcBwEZnn3s4B\n7gdOF5EVwOnuc0RkoIj8A8AdBPA7YJZ7u8c/MCAUG94chUue/B852cILY487dGzrniqu+Pt0yioq\nefbqwTYazaSWZVPghdHQ77sw4onEjhxb9T48dxEMvRHOui9x5zEhNTS8OdmsRROFsopKOhU3/cax\n1s3yeP76oXRs2ZTv/3MW01dv8yg6Y+rYuhIm/QDa94Nz/5j44cndhsPgH8D0v8LqDxN7LpMWLNFE\naN+BGrbsrvpWooGvk01py3yuedaSjUkBVXvgxSshyweXPwc5SdpI8Du/gVY94NUbYd+O5JzTpCxL\nNBFau90ZcdYxSKIBOKy5k2xK3GQzw5KN8YoqTL4Ztn4BI5+BFp2Sd+7cpnDxU7BnE0yxEWiZzhJN\nhMq2OYkmWIvGz0k2Q+jQognX/HMWM7+s9zqZMYnx2ROweBKcdjd0G5b885cMgJN/AQtfgkUTGy5v\nGi1LNBHyz6GpL9EAtGnehAljh9K+qAlXPzvTko1Jri8/hql3wdEXwAk/9i6Ok34KJcfCGz9xhlWb\njGSJJkJlFZUU5GZTXJDbYNk2zZsw4fqvk82sNZZsTBLsLIf/XAOtusOFf/V2bbJsH1z0FFRXwWs3\nO915JuNYoonQ2opKOhY3JcylGmhT6CSbdkVNuPoZSzYmwQ7uhxevcr7YLx/vrEHmtdbd4czfw6pp\nMOsfXkdjPGCJJkLBhjY3pE1hE164fiht3WQz25KNSZS3fgHr58JFT8JhR3gdzdcGXgvdvwPv/j/Y\nusLraEySWaKJgKpGlWggINkUNmGMJRuTCHPGwdxxcOJtcPT5XkfzTSJwwROQ08RZNaDmoNcRmSSy\nRBOBLburqKqupVOryBMNuN1oY79ONnO+smRj4qR8Dkz5GXQdBsN/7XU0wRW2h/MedVpcnzzsdTQm\niSzRRMA/4izUHJpwtHWTTZvCJox5ZhZzvgq5KZ0x4dmzBV66Cpq1c+bLZGV7HVFovS6CPpfDRw86\nWwqYjGCJJgLhDm1uSFt3gMBhzfPclo0lGxOlmmp4+Rqo3ObM/G9a7HVEDTv7QWjezulCO1DpdTQm\nCSzRRMC/PUBpy9iX8WhX5CSb1s1yGfPMTOaWWbIxUZj2W2c75vMehQ79vI4mPPkt4MInYdtKeO9u\nr6MxSWCJJgJlFZW0L2xCni8+XRPtipxutFbNchnz9EzWbE29lapNCls8Cf73Jxh0HfS7wutoItP1\nFGd155lPOVsXmEbNEk0E/HNo4ql9UT7PXz8UBH7x8gJqa21CmwnD5qXw6k1QOhjOTNOl+E+7Cw47\nyvk5Km1gTGNmiSYC0Q5tbkhJi3zuOq8nM9dUMO6zNXGv3zQy+3fCC1dCbgFc9i/wNbxKRUrKyXcW\n3qzcCm/eZqsGJJGIPCMim0VkUcCxFwM2QVsjIkG3SXVfW+iWmx3O+SzRhGn/wRo27Qq+PUA8jDy2\nlFOPPIwH3l5mXWgmtNpamPRD2PEVXDbOGTKcztr3hVN/6XQDLnzZ62gyyT+BswIPqOrlqtpPVfsB\nrwD1rYQ6zC07MJyTWaIJU7m7PUC0c2gaIiLcd/Ex5GRn8YtXrAvNhPDpw7D8TTjjXjj8eK+jiY8T\nbnW6AKf81FmnzSScqn4MBO2vFGd9rcuACfE6nyWaMMVjDk1D2hfl8//O68nMLyt4bvpXCTuPSVNf\nfgzv3wvHXAZDfuB1NPGT7YOL/uYM1X71RqfVZrx0ErBJVUOtFaTAuyIyR0TGhlOhJZowhbMPTTxc\nemwppxxxGPe/tezQOY0BYPrfoFlbOP9xb1dkToRW3eDMe+HLj5yRaCZWPhGZHXALKyG4RlN/a+YE\nVR0AnA3cJCInN1ShJZowlVXso2luNq3C2B4gFv4uNF+W8ItX5lsXmnHs3wkrp0Lvi53dKxujY6+G\nHmc6c2u2LPc6mnRXraoDA25hZW8R8QEXAy+GKqOq6937zcAkYHBD9YaVaETkJyKyWEQWicgEEWki\nIl1EZIaIrHBHK+S6ZfPc5yvd1zsH1PNL9/hyETkznHOnCv+Is3C3B4hFhxb5/Pq8o5m+uoLxM6wL\nzQDLpkDNAeh1sdeRJI4IXPBnyGkKE6+H6gNeR5SJvgMsU9WgF8tEpEBEmvsfA2cAi4KVDdRgohGR\nEuBHwEBV7Q1kA6OAB4BHVbUHsB241n3LtcB2Ve0OPOqWQ0R6uu/rhTPa4a8iksKLMn1TIubQ1Oey\ngR05qUdr7ntrGWsrrAst4y2eCEWdoDSsQT7pq7nbNbhhPnz8oNfRNFoiMgH4DDhSRMpFxP/9PYo6\n3WYi0kFEprhP2wKfish8YCbwpqq+3dD5wu068wH5brOqKbABGA74xyOOAy50H49wn+O+fpo7imEE\n8IKqVqnql8BKwmhypYJYtgeIlohw/yV9yBKxiZyZrrICVr0PvS5sfNdmgul5AfS9wlnhuWK119E0\nSqo6WlXbq2qOqpaq6tPu8atV9W91yq5X1XPcx6tVta9766Wq94ZzvgYTjaquA/4IlOEkmJ3AHGCH\nqla7xcqBEvdxCbDWfW+1W75V4PEg7zlERMb6L2BVV1fXfdkTW/ccYN/BmqQmGnAmct557tF8tnob\nz88sS+q5TQpZ9gbUVjvXZzLFaXeBZMGsp72OxMRBOF1nLXFaI12ADkABzmiDuvx/cgf7k0vrOf7N\nA6pP+S9g+Xy+hsJLinit2hyNUYM6cmL31tw3Zal1oWWqRROhZRdonyaLZsZDYXs4+gL4/Dk4YBOY\n0104XWffAb5U1S2qehBntujxQAu3Kw2gFFjvPi4HOsKhEQxFOBODDh0P8p6UtjYJc2hCcbrQjgHg\njokLUFumI7Ps3erMn+l9cWZ0mwUaPNYZbbfwP15HYmIUTqIpA4aKSFP3WstpwBLgA2CkW2YM8Jr7\neLL7HPf199X5dpwMjHJHpXUBeuBcTEp5/hZNPLYHiEZpy6b86tyj+e/KbUyYubbhN5jGY8lroDWN\ne7RZKJ2GQttjYObfbR20NBfONZoZOBf15wIL3fc8BdwO3CYiK3Guwfg7U58GWrnHbwPucOtZDLyE\nk6TeBm5S1Zq4/jQJUlZRSbvCJjTJ8W6Q3BWDO3F8t1bc++aSQ8vhmAyweBK0PgLa9vI6kuQTgcHX\nw6ZF8NX/vI7GxEBSuSumoKBA9+71vn/2sv/7DBReuuE4T+NYW1HJmY99zLGHt+Rf3x+clDk9xkO7\nN8LDR8Ept8OwX3odjTcOVMIjR0PXU51FRE1YRKRSVQu8jsPPVgYIQ7Ln0ITSsbgpvzznaD5ZsZUX\nZ1kXWqO35DVAM2u0WV25TWHAVbD0ddiVFpd0TRCWaBqw/2ANG3ft92TEWTBXDu7EcV1b8fs3l7Ju\nxz6vwzGJtGgitOkFhx3pdSTeGngtaC3MftbrSEyULNE0oHz7PlShUytvBgLUlZUlPDiyD7Wq/HLi\nQhuF1ljtLIe106H3RV5H4r3iLnDEmTDnWaiu8joaEwVLNA1Y6+EcmlA6FjfljrOP4uMvtvCf2bZ/\nR6O0+FXnPhNHmwUzeCzs3eJ2J5p0Y4mmAcnYhyYa3x1yOEO6FPO7N5awYad1oTU6iyc6u0+26uZ1\nJKmh6zBo1d22EEhTlmgaUFZRSZOcLA5rlud1KN/g70KrrrUutEZn+xpYN8daM4GysmDQ9VA+C9bN\n9ToaEyFLNA1I5vYAkTq8VQG3n3UkHy7fwstzrAut0Vg8ybnvZddnvqHfaMgpcCZwmrRiiaYBa5O8\nanOkvndcZwZ3LuaeN5awced+r8Mx8bBoIpQMhJaHex1JamlS5CSbRa84S/OYtGGJph7+7QFS7fpM\nIH8X2sGaWn41ybrQ0t62VbBxQWbPnanPoOuhpgrm/svrSEwELNHUY9veA1QeSP72AJHq3LqA2886\niveXbWbi3HVeh2NisWiic9/zwvrLZao2R0GXk2H2M1CTGtuImIZZoqmHl9sDRGqM24X229cXs2mX\ndaGlrcUTodNxUPStrZqM3+CxsHMtfPGW15GkLRF5RkQ2i8iigGO/EZF1IjLPvZ0T4r1nichyEVkp\nIneEcz5LNPVIxTk0oWRlCQ+M7MOBmlp+ZaPQ0tPmZbB5iY02a8gRZ0NRRxvqHJt/AmcFOf6oqvZz\nb1Pqvigi2cBfcPYk6wmMFpGeDZ3MEk09yrb5twdI/UQD0KV1AT8/8yimLdvMpM+tCy3tLJ7o7CrZ\nc4TXkaS2bB8M/L6zT8/mZV5Hk5ZU9WOcfcIiNRhY6W7pfAB4AWdjzHpZoqlHWUUlbZrnkZ/r3fYA\nkbr6+M4MPLwlv5m8mM3WhZY+VJ3rM4efAM3beh1N6hvwPcjOg1k21DkEn4jMDriNDfN9N4vIArdr\nrWWQ10uAwBV9y91j9bJEU4+yFB/aHEy2OwqtqtpGoaWVTYtg2wobbRaugtbQ+xKYN8HZhdPUVa2q\nAwNu4fQzPgl0A/oBG4CHg5QJNqGwwS8ZSzT1SPU5NKF0PawZPz/zSN5bupnX5tnS6mlh0USQbDja\nus3CNmQsHNzrJBsTM1XdpKo1qloL/B2nm6yucqBjwPNSoMEvGUs0IVRV17Bh1/6UnkNTn2tO6MKA\nTi24e/JiNu+2LrSUpupcn+l6ChS08jqa9NGhP5QOcrrPamu9jibtiUj7gKcXAYuCFJsF9BCRLiKS\nC4wCJjdUtyWaENb5twdI00STnSU8dGlf9h2s4c5Ji6wLLZWt/9xZ38xGm0Vu8FjYthJWf+B1JGlF\nRCYAnwFHiki5iFwLPCgiC0VkATAM+IlbtoOITAFQ1WrgZuAdYCnwkqoubuh8lmhCODSHplV6JhqA\nboc142dnHMHUJZuYPN+60FLW4omQlQNHn+d1JOmn5wgoOMzWP4uQqo5W1faqmqOqpar6tKpeparH\nqGofVb1AVTe4Zder6jkB752iqkeoajdVvTec81miCSGd5tDU59oTu9Lf7ULbsts2jUo5qs7eM92G\nQ36wQT6mXr48OPYa+OJtqPjS62hMCJZoQiirqCTPl3rbA0QqO0t4aGQfKg/U8OtXbRRayimf5cxy\nt9Fm0Rt4jTP/aPbTXkdiQrBEE4J/Mc2srNTbHiBS3ds057bTj+CdxZt4Y8EGr8MxgRZNdOaDHBl0\ntQ8TjsIOcPT5MPc5OFDpdTQmCEs0IZRV7Ev7brNA153Yhb4dW3DXa4vYuse60FJCbS0seRV6nA5N\nCr2OJr0NHgv7d8Cil72OxARhiSYIVU3bOTSh+LKz+OPIPuytquGu14KNWjRJV/YZ7N5gG5zFw+HH\nQ9veMOMp57qXSSmWaILYXnmQPVXVaTuHJpQebZtz6+k9mLJwI29aF5r3Fk8EXz4cEWxtQxMRERh8\nPWxaCGXTvY7G1GGJJoh02h4gUmNP6kqf0iL+32uL2GZdaN6pqYYlr8ERZ0JeM6+jaRyOudTZhdNW\ndU45lmiCaMyJxpedxUMj+7JnfzV3TW5wnpVJlK8+hb1bbLRZPOUWQP+rYOlk2GUt9lRiiSYI/xya\njsX5HkeSGEe2a86Pv9ODNxdsYMpC+0B6YtFEyG0GPc7wOpLGZdC1UFsDc571OhITwBJNEGXbKmnd\nLI+muT6vQ0mYH5zclWNKivh/ry6iYu8Br8PJLDUHnb+6jzwbchrnHzOeKe7qJO/Zz0K1/V6nCks0\nQTjbAzTuLwBfdhYPXdqHXfsPcrd1oSXX6o9g33ZnmXsTf4PHwt7NTjI3KcESTRDpuA9NNI5qV8iP\nhvfg9fnreXuRdaElzeKJkFfkLDtj4q/bcCjuZoMCUoglmjoOVNeyYWfjmqxZnxtO7UavDoX8+tVF\nbLcutMSrroKlbzgLaPrSe3mjlJWV5Qx1XjsD1s/zOhqDJZpvWb9jH7VKo5tDE0pOdhZ/vLQvOyoP\n8pvXrQst4VZOg6qdtiVAovUdDTkFtqpziggr0YhICxF5WUSWichSETlORIpFZKqIrHDvW7plRUT+\nJCIr3b2nBwTUM8Ytv0JExiTqh4pFYx7aHMrR7Qu5ZXgPXpu3nncWb/Q6nMZt8UTIL3Y2OTOJk98C\n+l4OC/8DlRVeR5Pxwm3RPA68rapHAX1xNry5A5imqj2Aae5zgLOBHu5tLM4+1IhIMXA3MARni9C7\n/ckplTSGfWiiceOwbvRsX8idkxaxo9K60BLi4D5Y/pazAGR2jtfRNH6DroeaKpj7L68jyXgNJhoR\nKQROBp4GUNUDqroDGAGMc4uNAy50H48A/qWO6UALd4vQM4GpqlqhqtuBqUDKrb2xtqKS3Ows2jZv\n4nUoSZXjjkLbUXmA376+xOtwGqcV78KBPTZJM1na9oTOJ8Gsp525NeYQEXlGRDaLyKKAYw+5vVYL\nRGSSiLQI8d417k6c80RkdjjnC6dF0xXYAjwrIp+LyD9EpABoG7AD2wagjVu+BFgb8P5y91io43V/\niLEiMltEZldXV4fzM8RVWUUlpcX5jWJ7gEj16lDETcO6M+nzdUxdssnrcBqfRROd3SAPP9HrSDLH\n4LGws8zZGM0E+iff/kN/KtBbVfsAXwC/rOf9w1S1n6oODOdk4SQaHzAAeFJV+wN7+bqbLJhg39Ba\nz/FvHlB9SlUHqupAny/5EyYzZWhzKDcN685R7Zrzq0kLrQstnqr2wBfvOFsPZzfeicAp58hzoLDE\nhjrXoaofAxV1jr2rqv6/7qcDpfE6XziJphwoV9UZ7vOXcRLPJrdLDPd+c0D5jgHvLwXW13M8Zagq\nZdsyO9Hk+pxRaBV7D3DPG9aFFjdfvA3V+2y0WbJl+2Dg92H1h7BludfRJJPP3zPk3sZG+P7vA2+F\neE2Bd0VkTrj1NphoVHUjsFZEjnQPnQYsASYD/pFjY4DX3MeTge+5o8+GAjvdrrV3gDNEpKU7COAM\n91jK2LnvILurqjM60QD0LiniplO7MXHuOt6zLrT4WDwJmreHTsd5HUnmOfZqyM7NtKHO1f6eIfcW\ndpNORO4EqoHxIYqcoKoDcAZ+3SQiJzdUZ7ijzm4BxovIAqAf8AfgfuB0EVkBnO4+B5gCrAZWAn8H\nbgRQ1Qrgd8As93aPeyxllB1aTDOzEw3AzcN7HOpC21l50Otw0tv+XbBiKvS80JlMaJKroDUccxl8\n/hzstuH79XGnnZwHXKkafAc5VV3v3m8GJuGMIq5XWL/1qjrPzYp9VPVCVd2uqttU9TRV7eHeV7hl\nVVVvUtVuqnqMqs4OqOcZVe3u3lJuedVMnEMTSq7P2U5gm3WhxW75FGeYrY02887JP3UWM/30Ua8j\nSVkichZwO3CBqlaGKFMgIs39j3F6phrcstf+vApgLZpvOqa0iB+e0o1X5pbzwbLNDb/BBLdoIhR1\nhNJBXkeSuYq7Qr8rnFWdd67zOhrPicgE4DPgSBEpF5FrgSeA5sBUd+jy39yyHURkivvWtsCnIjIf\nmAm8qaoNDumz4S8B1lZU0qogl2Z59s/id8tp3ZmycAP3v7WMU444LCOHfcdk33ZY9T4MvcHZbth4\n5+Sfw/wX4JOH4bxHvI7GU6o6Osjhp0OUXQ+c4z5ejTNpPyLWoglQVlFprZk68nzZ/Oi0HizftNuW\np4nG0jeg9qCNNksFLQ+HAVc5KwXsKPM6moxiiSZAps+hCeX8vh3o2rqAx6etoLY26PVBE8riidCy\nM3To73UkBuCknzkty48f8jqSjGKJxnWwppb1O/ZbogkiO0u45bTuLNu4m3eXWKsmbHu3OZuc9brY\nus1SRVEJHHsNfD4eKlZ7HU3GsETj2rBjPzW1aokmhPP7OK2ax96zVk3YFv4HtAaOGel1JCbQSbc5\ni5p+ZK2aZLFE47IRZ/XzZWdx83B/q8YmcYZl3nho3w/a9vI6EhOoeTsYdB0seAG2rvA6moxgicaV\nqdsDROKCvh3oYtdqwrNxIWxcAP2/63UkJpgTbgVfPnz0gNeRZARLNK6yikpysoV2hZm1PUAkfNlZ\n3DK8O0s37GLqUmvV1Gve886yJ70v8ToSE0yzw2DIWFj4Mmxe6nU0jZ4lGtfaikpKWzYl2+aJ1OuC\nvh3o3Kopj7+3ghArVJjqA7DgRWfl4KbFXkdjQjn+R5DbDD68v+GyJiaWaFw2hyY8TqumB0s27LI9\na0JZ8S5UboN+V3odialP02IY+kNY8qrT1WkSxhKNy5lDk+91GGlhRD+3VTPNWjVBzRsPzdpBt+Fe\nR2IactyNkFdkrZoEs0QD7Kw8yM59B21oc5icEWg9WLx+F+8ttTXQvmHPZmeDs76X2wZn6SC/JRx/\nMyx7A9Z/7nU0jZYlGmDtdlu1OVIX9uvA4a2a8th7X1irJtCCl5y5M9Ztlj6G3OAknA/u8zqSRssS\nDTaHJhq+7CxuHtadxet3Mc1aNQ5Vp9usZCAcdmTD5U1qaFLoDAxY8Q6sneV1NI2SJRos0UTrov4l\ndCpuymPTrFUDwIZ5sHmJsxy9SS+Dx0LT1vDhH7yOpFGyRIOTaFo2zaGwSY7XoaQV/2oBi9ZZqwZw\n587k2dyZdJTXDE681dnS4avPvI6m0bFEgzOHxq7PROei/iV0LM63EWjVVc7aZkefB/ktvI7GRGPg\ntdCsLXxwr9eRJJyIPCMim0VkUcCxYhGZKiIr3PuWId47xi2zwt36uUGWaLA5NLHIyc7ilmE9WLhu\nJ+9n8i6cy99yNjmzQQDpK7cpnHgbrPkEvvzY62gS7Z/AWXWO3QFMU9UewDT3+TeISDFwNzAEGAzc\nHSohBcr4RFNdU8u67fusRRODiwZYq4Z5z0NhCXQ91etITCyOvRqad4AP/uAM7mikVPVjoKLO4RHA\nOPfxOODCIG89E5iqqhWquh2YyrcT1rdkfKLZsHM/1bY9QExy3BFoC8p38sHy5LZqnvp4FT//z/yk\nnvNbdm+ElVOh7yjIyvY2FhObnCZw8k+h7DPnek1maauqGwDc+zZBypQAawOel7vH6pXxiWZthc2h\niYeLB5RS2jI/qWugPfbeF/xhyjL+M6ecDTv3JeWcQS14EbQW+tpos0ah/1VQ1NG5VpO+rRqfiMwO\nuI2NU73BFoNs8B8p4xONDW2OD3+rZn75Tj5cviXh53v8vRU89t4Kju/WCoAZq+v2AiSJqrNbY8eh\n0Lq7NzGY+PLlwck/h3VznHXr0lO1qg4MuD0Vxns2iUh7APc+WPdEOdAx4HkpsL6hii3RVFTiyxLa\nF9n2ALHyt2oeS/C1mj9NW8Gj733BJQNKGff9wRTl5zB99baEna9e6+bC1uU2d6ax6XcFtOyc7q2a\nSE0G/KPIxgCvBSnzDnCGiLR0BwGc4R6rlyWaikpKWubjy874f4qY5fqyuGlYd+av3cGHXySmVfPn\naSt4ZOoXXDyghAdH9iEnO4tBnYu9SzTz/u1soNXrIm/ObxIjOwdOuR02zIdlb3odTdyJyATgM+BI\nESkXkWuB+4HTRWQFcLr7HBEZKCL/AFDVCuB3wCz3do97rF4Z/+1qc2ji65IBpZS0SMy1mifeX8HD\nU7/g4v4lPDSy76G9g4Z2LWbNtsrkX6c5uB8WvgI9L3CWMTGNyzGXQavuzgi02lqvo4krVR2tqu1V\nNUdVS1X1aVXdpqqnqWoP977CLTtbVa8LeO8zqtrdvT0bzvkyPtHYHJr48rdq5q3dwUdxbNX85YOV\n/PFdN8lc2vcbG9QN7erRdZplb0DVTps701hl++CUO2DzYlgarBfJhCujE82u/QfZXmnbA8TbyGPd\nVk2crtX85YOVPPTOci4KkmTGQyGDAAAgAElEQVQAjm5fSGETHzO+THL32bznoagTdD4puec1ydP7\nYjjsKGdl59oar6NJWxmdaGxoc2Lk+rK4cVg3Pi/bwccrtsZU118/dJLMhf068McgSQYgO0sY3KWY\n6cls0exc58yz6DcasjL6Y9S4ZWXDqXc4Az4WTfQ6mrSV0Z8QSzSJc+mxHelQ1ITHY9iv5skPV/Hg\n28sZ0a8DD1/WL2iS8RvatRVfbt3Lxp37ow05MgteANSZpGkat6NHQNve8OF9UFPtdTRpKaMTjX8O\nTadWlmjiLdeXxU3DuzO3bAefRNGq+dtHq3jg7WWM6NeBRxpIMhBwnSYZ3Wf+uTOHnwDFXRN/PuOt\nrCw49ZdQsQoWvuR1NGkp4xNNC9seIGEOtWoivFbzfx+t4v63lnFB3w48HKK7rK6j2xfSvIkvOd1n\na2c6Xzo2CCBzHHUutO8LHz0ANQe9jibthJ1oRCRbRD4XkTfc511EZIa7VPSLIpLrHs9zn690X+8c\nUMcv3ePLReTMeP8wkSqrsMU0E8m5VtOdOV9t59OV4bVq/u+jVdz31jLO79uBRy7rG/b8puwsYXDn\nYmYkYz7NvH9DTgH0HJH4c5nUIALD7oTta5xBICYikbRofgwsDXj+APCou6T0duBa9/i1wHZV7Q48\n6pZDRHoCo4BeOKt9/lVEPF2BcK0NbU64SweW0r6oSVjzap76+Osk82gEScZvaNdWrN66l027Enid\n5kAlLJoEvS50NssymaPHGc423R8/5Ow/ZMIW1idZREqBc4F/uM8FGA687BYJXFI6cKnpl4HT3PIj\ngBdUtUpVvwRW4uxn4ImaWqV8u03WTLQ8XzY3DuvO7K+289+VoVsbf/94NX+Ysozz+rSPKsnA19dp\nErpKwNLX4cBuW3ImE4nAsF/BzrXw+XNeR5NWwv00Pwb8AvBPj20F7FBV/xCMwKWiDy0j7b6+0y0f\n1fLSibJx134O1tj2AMlwmb9VMy34CLR/fLKae6cs5dw+7Xns8n5RLwfUs0MhzfN8zPgygddp5o13\n1sDqdHzizmFSV7fh0Ok4+PhhZ2UIE5YGP9Eich6wWVXnBB4OUlQbeC2s5aVFZKx/aevq6sQNJSzb\nZkObkyXPl82Np3Zj1prt/G/VN1sb//hkNb9/cynnHtOex2NIMuBcpxnUJYHrnu0oc3Ze7HuFzZ3J\nVP5Wze71MPsZr6NJG+F8Wk4ALhCRNcALOF1mjwEtRMTnlglcKvrQMtLu60U4O7mFtby0qj7lX9ra\n5/PVfTlubA5Ncl02qCPtCr95rcafZM45ph2PjYotyfgN7VrM6i172ZyI6zTz3bkz/UbHv26TPrqc\n7Nw+eRiq9ngdTVpo8JOtqr90F13rjHMx/31VvRL4ABjpFgtcUjpwqemRbnl1j49yR6V1AXoAM+P2\nk0SorKKSbNseIGmcazXdmLmmgs9WbePpT7/k928u5eze7Xh8VH9y4rR69qHrNPHuPqutdbrNupwM\nLTrFt26Tfk67Gyq3wvQnvY4kLcTy6b4duE1EVuJcg3naPf400Mo9fhtwB4CqLgZeApYAbwM3qapn\niweVVVRS0sK2B0imywZ2pG1hHre+OI/fvbGEs3u340+j45dkAHq2d6/TxLv7rOwzZ2hrv+/Gt16T\nnkoHwpHnwv/+BJUebbqXRiL6hKvqh6p6nvt4taoOdpeKvlRVq9zj+93n3d3XVwe8/15V7aaqR6rq\nW/H9USJTZtsDJF2TnGxuPLU7m3dXJSTJAPiysxJznWbeeMhtDkefH996Tfoa/muo2g3/fczrSFJe\nxv45b3NovPHdoYfzr+8PTkiS8RvSpZhVW/ayeXecrtNU7YHFr0LviyDXfmeMq21P6HMZzHgKdm3w\nOpqUlpGJZk9VNdv2HrAWjQeys4STjzgsYUkGErA/zZLX4OBeW3LGfNupv4Tag84kzjQiIkeKyLyA\n2y4RubVOmVNFZGdAmbuiPV9GJhobcda49epQSLM8X/y6z+Y9D8XdoOOQ+NRnGo/iLjBgDMwdBxVf\neh1N2FR1uar2U9V+wLFAJTApSNFP/OVU9Z5oz5eRiabMEk2j5svOYlDnlvGZuFnxJXz1qbMSgDS8\nuKfJQCf/HLJynG0E0tNpwCpV/SpRJ8jIRGMtmsZvSNdWrNy8hy27Y1yTav4EQKCvzZ0xIRS2hyFj\nYcFLsGmJ19FEYxQwIcRrx4nIfBF5S0R6RXuCjEw0ZRWVFDbxUdTUtgdorOKyP01tLcybAN2GQZFn\nqyWZdHDCrZDXHN7/vdeR+Pn8K6y4t7HBCrmr7l8A/CfIy3OBw1W1L/Bn4NVog8nYRGObnTVuvTsU\nUpCbHdt1mjWfwM4yGwRgGta0GI7/ESx/E8pnex0NQLV/hRX39lSIcmcDc1V1U90XVHWXqu5xH08B\nckSkdTTBZG6isW6zRs0/nyamkWfznoe8ImfTK2MaMvQGaNoapv3W60giMZoQ3WYi0s5deR8RGYyT\nL6L6yy3jEk1trVJesc/m0GSAIV1asWLzHrbuieI6zf5dzrDm3hdDTn78gzONT15zOPlnzsKrqz/0\nOpoGiUhT4HRgYsCxG0TkBvfpSGCRiMwH/gSM0ki2yg2QcYlm0+79HKiptRZNBhjatRiIcj7Nkleh\neh/0tyVnTASOvQYKS2HaPRDdd3LSqGqlqrZS1Z0Bx/6mqn9zHz+hqr1Uta+qDlXV/0V7roxLNLY9\nQOboXVIU/XWaz8dD6yOg5Nj4B2Yar5wmcOodsG4OLHvT62hSRuYlGhvanDFysrMY2Lk48pFn21bB\n2unOIACbO2Mi1Xc0tOoB7/8Oaj1bNzilZFyiWVtRSZZAhxbW754JhnQt5otNEV6nmTceJAv6XJ64\nwEzjle2D4XfClmWwMNio4cyTcYmmrKKSDi3yE7rWlkkd/vk0M8NdJaC2xtngrPt3nIl4xkTj6BHQ\nrg988AeoPuB1NJ7LuG9bG9qcWY4pKaJpJNdpVn8Iu9Y5S84YE62sLGdztB1fOeugZbgMTDT7LNFk\nkEPXacIdeTbveWjSAo48J7GBmcav+2nQ6XhnZecDe72OxlMZlWgqD1SzdU+VzaHJMEO6FLN80262\nNXSdZvtXzrDmvqPAl5ec4EzjJQKn3QV7NsHMUBPzM0NGJZq1FfsAG3GWacK+TvPxgyDZcMKPkxCV\nyQiHHwc9zoBPH4N9O7yOxjMZlWhsaHNm6lNaRH5OA9dptq1yFtAc+H0o7JC84EzjN/zXsH8H/O/P\nXkfiGUs0ptFzrtO0ZHp912k+ehCyc+HEnyQvMJMZ2veFXhfD9Cdhz2avo/FERiWatRWVNM/z0cK2\nB8g4Q7u2Yvmm3VTsDTLUdMsXsPAlGHwdNG+b/OBM4zfsTqjeD5887HUknsioRFNWUUnH4qaIzfbO\nOP51z2YGWyXgo/vBl+/sKWJMIrTuDv2vhNnPwI4yr6NJuoxLNNZtlpmOKWnhXqep0322aQksmghD\nfgAFUW21YUx4TrkdEPjwAa8jSbqMSTS1tcpa2/AsY+X6/Ndp6rRoPrwPcpvB8bd4E5jJHEWlMOg6\nmP+8012bQTIm0WzZU0VVda3NoclgQ7oUs2zjbrb7r9NsWABLJ8NxNzo7JBqTaCfdBjlN4YOU2fI5\nKTIm0diIM+OfTzPDP5/mw/ugSREMvdHDqExGKWgNx93kbKq3/nOvo0mazEk0tg9NxutT2oImOVlO\n99m6ubB8Chx3C+S38Do0k0mOuwnyW8K033kahoisEZGFIjJPRGYHeV1E5E8islJEFojIgGjPlTmJ\npqISESix7QEyVq4vi4GHFzuJ5oM/OB/2IT/wOiyTaZoUwYm3wappsOZTr6MZpqr9VHVgkNfOBnq4\nt7HAk9GeJGMSzdqKSjoU5ZPry5gf2QQxtGsxyzfuZseK/zlLzTQp9Dokk4kGXw/N26f6ls8jgH+p\nYzrQQkSi2jsjI751a2uVWV9V0L1NM69DMR4b0rUVCszIGQyDx3odjslUOflwyi9g7QxY8W4izuAT\nkdkBt2C/7Aq8KyJzQrxeAqwNeF7uHotYRiSa/67aytqKfVxybKnXoRiP9alZQhOqmN7mMsgt8Doc\nk8n6XwUtuzjXampr4117taoODLgFWz76BFUdgNNFdpOInFzn9WAz26NqfmVEopkws4yWTXM4s5ct\nL5LRVMn7+D6OzfmK6ZVR/WFmTPxk5zhL02xaCIsnJv30qrrevd8MTAIG1ylSDnQMeF4KrI/mXI0+\n0WzZXcW7izcx8thS8nzZXodjvPTlR/DVpwzt3pZlm/awo9K22DUe630JtOkFH9yb1C2fRaRARJr7\nHwNnAIvqFJsMfM8dfTYU2KmqG6I5X4OJRkQ6isgHIrJURBaLyI/d48UiMlVEVrj3Ld3jIYfEicgY\nt/wKERkTTcCRenlOOdW1yqjBnZJxOpOqVJ2RZoUlDDl+GKph7E9jTKJlZcF3fgMVq+Ht25N55rbA\npyIyH5gJvKmqb4vIDSJyg1tmCrAaWAn8HYh6wpkvjDLVwE9Vda6bAeeIyFTgamCaqt4vIncAdwC3\n880hcUNwhsQNEZFi4G5gIE4/3xwRmayq26MNviG1tcoLs8oY0qWYbofZQICMtnKac+H13Efo26UN\neb4spq+u4Ixe7byOzGS6I85wFnT972PQpqczIi3BVHU10DfI8b8FPFbgpnicr8EWjapuUNW57uPd\nwFKckQcjgHFusXHAhe7jUEPizgSmqmqFm1ymAmfF44cI5bPV2/hqWyVXDLHWTEZTdbomijpB/6vI\n82Vz7OFB1j0zxiun3QVHnAVv3Q6rP/Q6mriL6BqNiHQG+gMzgLb+/jr3vo1bLNSQuLgNlQvX8zPL\naNE0hzPtr9bM9sXbsH4unPJz8OUCznI0SzfuYmflQY+DMwbIyoaL/w6tj4CXxjg7vjYiYScaEWkG\nvALcqqq76isa5JjWc7zuecb6x35XV1eHG963bN1TxbuLN3LJgFKa5NgggIzlb8207AJ9Rx86PKRL\nsXOdZo1dpzEpokkhjJ4AkgUTRsP+nV5HFDdhJRoRycFJMuNV1T8Ob5N/lqh779+jNNSQuLCGyqnq\nU/6x3z5fOJeQgntlTjkHa5TRgzs2XNg0Xktfh40Lnb1Asr/eWbVvxxbudRrrPjMppLgLXPYvqFgF\nL18LtTVeRxQX4Yw6E+BpYKmqPhLw0mTAP3JsDPBawPFgQ+LeAc4QkZbuCLUz3GNxp6pMmFnG4M7F\ndG/TPBGnMOmgttZZoblVDzjm0m+81CQnmwGd7DqNSUFdToKzH4SVU+G9u72OJi7CadGcAFwFDHdX\n+ZwnIucA9wOni8gK4HT3OYQYEqeqFcDvgFnu7R73WNx9tnoba7ZVMnqItWYy2pJJsHkJnHoHZH+7\ndTy0ayuWbLDrNCYFDboWBl0P//szzHve62hi1mDflKp+SvDrKwCnBSkfckicqj4DPBNJgNGYMHMt\nRfk5nN07qvXfTGNQWwMf3g+HHQ29Lg5aZEjXYvQ9mLWmgu/0tFUjTIo56z7Yuhxe/zEUd4NOQ7yO\nKGqNbmWAbXuqeGfRRi4eUGKDADLZwpdh6xcw7JfOpLgg+nVsQa5dpzGpKjsHLh0HhSXw4ndhZ7nX\nEUWt0SWaiXPXcaCmltG2EkDmqqmGj+6HtsfAUeeHLOZcp2nB9C8t0ZgU1bQYrngRqvc7I9EO7PU6\noqg0qkTjHwQw8PCWHNHWBgFkrAUvOEt6DPtVyNaM39CurVi8fhc799l1GpOiDjsSLnnaGT356g8T\nsdJzwjWqRDPjywpWb91rrZlMVnMQPnoAOvSHI89usPiQLq1Qhdk2n8aksiPOgNPvgSWvwccPeh1N\nxBpVopkws4zCJj7O7WODADLW5/+GHWXO8usSagzL1/p3sus0Jk0cfwv0vcIZsr/4Va+jiUijSTTb\n9x7grYUbudhWAshc1VXw8R+hdBB0/05Yb2mSk03/ji2YvtpaNCbFicD5j0HpYJh0A2yY73VEYWs0\nieaVueUcqKlllK0EkLnm/gt2lYfdmvFzrtPsZNd+u05jUpwvDy7/tzNIYMIVsGdzw+9JAY0i0fgH\nAQzo1IKj2hV6HY7xwsF98MnD0Ol46HpqRG8d0rWYWrtOY9JF87bOmmj7KuCFK52WfIprFIlm1prt\nrNpigwAy2uxnYfcGGB5ZawZgQKeW5GZnWfeZSR/t+8KFT0L5THj9Vmfx2BTWKBLNhJllNG/i47w+\nHbwOxXjhwF749BHocjJ0PjHitzfJyaZfpxY2IMCkl14Xwil3wPzn4bMnInprqJ2T65Q5VUR2Biw9\ndle0oUa/PHKK2FF5gDcXbmDUoI7k59oggIw06x+wdwsM+3fUVQzt2oon3l/Brv0HKWyS0/AbjEkF\np9wOW5bC1Lug9ZHOMOjwBN05WVWX1Cn3iaqeF2uYad+imTh3HQeqaxk1yLrNMlLVbvj0Meh2GnQa\nGnU1Q7vYdRqThrKynC60tr3glWthy/Kw3lbPzsmJCTNRFSeDfxBAv44t6NnBBgFkpI8ecC6KDrsz\npmr6u9dpZth1GpNucgtg1ARnRNrzl0NlBYDPv4Gkexsb6u11dk6u6zgRmS8ib4lIr2hDTOuuszlf\nbWfF5j08eEkfr0MxXpjxlLOM+rFXQ+mxMVWVn5tNv452ncakqRYd4fLxMO48+M/VANWqOrChtzWw\nc/Jc4HBV3eNuDfMq0COa8NK6RfP8zDKa5/k4r6+tBJBxFr0Cb/0CjjwXznk4LlUO7VrMwnU72W3z\naUw66jQEzn8cvvworOIhdk4+RFV3qeoe9/EUIEdEWkcTWtommp2VB3lzwQYu7F9C09y0bpiZSK16\nHyb+ADodByOfDrqpWTSGdm3lXqfZHpf6jEm6flfAcTc3WKyenZMDy7RzyyEig3HyRVRN/rRNNJM+\nL6eq2rYDyDjr5sAL33VWtB09AXLy41Z1/04tyckW3l60MW51GpN0p98TTqmgOyeLyA0icoNbZiSw\nSETmA38CRrkbW0YsLZsCziCAtfS1QQCZZesKGH8pFLSG774C+S3iWn1+bjajBnXiuelfUdoyn1tO\ni6o72hhvZTU8zaOBnZP9ZZ4AIpugE0JaJpq5ZTtYvmk3D1xyjNehmGTZtR6euwgkC66aBM3bJeQ0\nv7mgF3sPVPPw1C8QgZuHW7IxJlZpmWgmzCyjWZ6tBJAxKivguYth3w64+g1o1S1hp8rOEh4a2RcU\n/vjuF4gINw3rnrDzGZMJ0i7R7Nx3kDcWrOeSAaUU5KVd+CZSByphwiioWOV0l3Xol/BTZmcJD13a\nF4CH3nEmwFmyMSZ6afdN/dq8dew/aIMAMkLNQWdOwNqZcNk4Zy2zJPEnG8VJNiJw46mWbIyJRlol\nGlXl+Rll9CktondJkdfhmESqrYXJt8CKd+C8R6HniKSHkJ0l/PHSvqgqD77ttGws2RgTubRKNPPW\n7mDZxt3cd7ENAmj03rsL5k9wlpYZ+H3PwsjOEh6+rB8KPPj2cgThh6cm7hqRMY1RWiWaCTPLKMjN\n5vy+NgigUfvv487SMoPHwsk/9zoaJ9lc2hdVeODtZYjADadYsjEmXGmTaHbtP8jr852VAJrZIIDG\n6/PxzpLnvS6Gsx6IeBOzRPFlZ/HIZc41m/vfWoYAP7BkY0xY0uYb+7XP17HvYA1X2CCAxmv5W851\nma7D4KL/c5ZATyG+7CwevcwZjXbfW07LZuzJlmyMaUhaJBpVZfyMMnqXFHJMqQ0CaJS++swZYda+\nD1z+HPhyvY4oKH+yUVX+MGUZgnD9yV29DsuYlJYWiWZ++U6WbdzNvRf19joUkwibFsOEy6GoFK58\nGfKaex1RvXzZWTx2uTNA4N4pSxGB606yZGNMKGmRaCbMKKNpbjYX2CCAxmf7V/DvSyCnqbO0TEFU\nq5AnnS87i8cv7wcKv39zKWDJxphQUj7R7N5/kMnz1zOiXwea217ujcueLc76ZQcr4Zq3oUV6XX/z\nZWfx2Kh+KGrJxph6pHyieW3eevYdrLGVABqbqt0wfqSzWOb3XoW2Pb2OKCo52Vk8Pqo/qp/z+zeX\nIiJce2IXr8MyJqWk1rCeIJ6fUUbP9oX0sUEAjUd1FbxwJWxc6Cwt02mo1xHFJCc7iz+N7s/Zvdvx\nuzeW8MynX3odkjEpJemJRkTOEpHlIrJSRO6or6xm5bBkwy5GD+mEpMh8ChMlVaerbP3nMPF6Z7vZ\nEX+BI870OrK48Cebs3q14543lvDsfy3ZmNTW0HexiOSJyIvu6zNEpHO050pq15mIZAN/AU4HyoFZ\nIjJZVZcEK1+b25T8nGxG9LNBACmvag/sLIdd5c79znXffL5rPVTv/7r86b+DfqO9izcBcrKz+PMV\n/bn5+bn89vUlCHD1CdaNZlJPmN/F1wLbVbW7iIwCHgAuj+Z8yb5GMxhYqaqrAUTkBWAEEDTRaE5T\nzj+mDYXZ1XCwOolhmm9Qhb1bYJebPPy3Q8/Xwv6d33yPZEHz9lBYAu37wVHnQlFHZwhzq+7OVsyN\nUE52Fn8ePYCbn5/Lb153fq1Hxfn6Yp4vy1r4JlbhfBePAH7jPn4ZeEJEJJrtnJOdaEqAtQHPy4Eh\nIUuLMHrRdbBkVaLjMpHKbwmFpU7y6DTUSSBFHZ3EUlTq7ICZnZmjBHN9WTxxxQBucpONP+HEy8p7\nz8aXbYnGxCSc7+JDZVS1WkR2Aq2ArZGeLNmJJtin4xvZUUTGAmMB8gpb0e/0GxvY2dokRdNWbhLp\nCEUlkFvgdUQpLdeXxV+uGMDLc8rZue9gXOvOstaMaZhPRGYHPH9KVZ8KeN7gd3GYZcILJpo3xaAc\n6BjwvBRYH1jA/cd4CqCgoEDlpJ8kLzpj4ijXl8UVQ2xYvvFEtaoOrOf1Br+LA8qUi4gPKAIqogkm\n2aPOZgE9RKSLiOQCo4DJSY7BGGMyXTjfxZOBMe7jkcD70VyfgSS3aNx+vpuBd4Bs4BlVXZzMGIwx\nJtOF+i4WkXuA2ao6GXgaeE5EVuK0ZEZFez6JMkElRUFBge7du9frMIwxJq2ISKWqpsyF1JRfGcAY\nY0x6s0RjjDEmoSzRGGOMSShLNMYYYxLKEo0xxpiESulRZyKiQGUcq8wGalK0vlSOLd71pXJsqV5f\nKscW7/pSObZ41xfv2JqqauosIaGqKXsDtsS5vqdStb5Ujs1+1tSpL5Vjs581Nepy64vrd2est1Tv\nOtsR5/peT+H6Ujm2eNeXyrGlen2pHFu860vl2OJdX7xji/d3Z0xSvetstta/Xo8xxpg6Uu27M9Vb\nNE81XMQYY0wdKfXdmdItGmOMMekv1Vs0UQu2H7aIjHePLRKRZ0Qk7J25QtT3tIjMF5EFIvKyiDSL\npb6A1/4sIntijO2fIvKliMxzb/1irE9E5F4R+UJElorIj2Ks75OA2NaLyKsx1neaiMx16/tURLrH\nUNdwt65FIjLOXSI9nLqeEZHNIrIo4FixiEwVkRXufcsIfs5g9V0qIotFpFZEIuoaCVHfQyKyzP0d\nniQiLWKo63duPfNE5F0RCXsP9mD1Bbz2MxFREWkdS30i8hsRWRfwe3dOLLGJyC3u785iEXkwxthe\nDIhrjYjMC7e+tOD1aIQ6IyXOApYDK4E73GM3u88VaB1mPdnAKqArkAvMB3oC5+Bs5iPABOCHMdZX\nGFDmEX/M0dbnvjYQeA7YE2Ns/wRGRvF/EKq+a4B/AVluuTax/qwBZV4BvhdjfF8AR7tlbgT+GUNd\na4Ej3DL3ANeGGdvJwABgUcCxBwN+l+8AHojg/yJYfUcDRwIfAgMj/L8NVt8ZgM99/EC48YWoK/Dz\n8CPgb7HE5h7viLPC8FeE+fmvJ77fAD+L4jMRrK5hwHtAXiSfh/p+1oDXHwbuiqC+Z4DNdeLrB0wH\n5gGzgcGR/tzxvKVMi0ZEsoG/AGfjfNhHi0hP4L/Ad3B+0cJ1aD9sVT0AvACMUNUp6gJm4mz2E0t9\nu9zYBcgn/N3ngtbn/hs8BPwizHpC1hXB+8Ot74fAPapaC6Cqm+MRn4g0B4YD4bZoQtWnQKFbpohv\nb+IUbl2XAFWq+oVbZqp7rEGq+jHf3hhqBDDOfTwOuDCcukLVp6pLVXV5uHWEUd+7qlrtPp1OmJ+J\nEHXtCnhaQAS7MYb4twN4FOfzEFEffz31RSxEXT8E7lfVKrdMuJ+HemNzv0suw/lDOFz/xPkjPdCD\nwG9VtR9wl/vcMymTaAj9Zf65qq6JsK5g+2GX+J+I02V2FfB2rPWJyLPARuAo4M8x1nczMFlVN4RZ\nT72xAfe6XRmPikhejPV1Ay4Xkdki8paI9IhDfAAXAdPqfElFU991wBQRKcf5v70/yrraATkB3VIj\n+eZOhJFq6///dO/bxFBXon0feCuWCtzu1bXAlThfcLHUdQGwTlXnx1JPHTe7n4lnIunGDOII4CQR\nmSEiH4nIoDjFdxKwSVVXhPuGEIkrmj+8EiaVEk1DX0iRaGiv678CH6vqJ7HWp6rXAB2ApcDlMdSX\nB1xK+Mmqodh+iZP8BgHFwO0x1pcH7FdnyOTfcZrrsdTnN5rI/noLVd9PgHNUtRR4FqcrM5q6anE2\neHpURGYCu4HqIOUaFRG5E+fnHB9LPap6p6p2dOu5OYZ4mgJ3EmOyquNJnD+Y+gEbcLqoouUDWgJD\ngZ8DL7mtkVhF+nkI5VbgITfp/xHn+8AzqZRoGvpCikTI/bBF5G7gMOC2eNQHoKo1wIuE2cUSor41\nQHdgpYisAZqKs7NdVLGp6ga3l7AK54t3cAyxrXePv+IemwT0ibE+RKSVG9ebYdYVqr7NQF9VneEe\nexE4PtrYVPUzVT1JVQcDHwNh/3UZxCYRaQ/g3ofdxZIsIjIGOA+40u1WjofnCf/zEEw3oAsw3/08\nlAJzRaRdtBWq6iZVrXG7f/9O+J+JYMqBie5nbCbOHyhhD1YIxh10cjHO72+sfgj8xE36P8HZLdMz\nqZRo6v0yj1DQ/bBF5BE55YkAAASnSURBVDrgTGC0/1pDjPV1h0P9qucDy2Ko71VVbaeqnVW1M1Cp\nquGMnAoVm//LTXCuC3xrNE8k9eFcQxnuljkF5+J7LPWB04J7Q1X3h1lXffUVicgRbpnTcVqYUdUl\nIm0A3O7G24G/RRBfXYH7ro8BXouhrrgTkbNwfsYLVDWmdQXrdKdeQPifh29R1YWq2ibg81AODFDV\njTHE1z7g6UWE/5kI5tDnwf29ywW2xlAfONeil6lqeYz1gPO7NtF9/B9iS6qx83IkQuANpym6Guev\nGP8IoF4Br68hslEn5+B8Ga4C7nSPVbvP57m3SEZ2fKM+nCT9X2Ahzi/seAJG3UQTX53Xwxp1Vs/P\n+n5AbP8GmsVYXwuclsdC4DOcFkRMPyvOyKmzovhdCRbfRW5s8916u8ZQ10M4iWo5cGsEcU3A6ZI5\niPPFeC3QCpiG0yqaBhTHWN9F7uMqYBPwToz1rcTpsvZ/JsIaKRairlfc37cFOEuqlMQSW53X1xDZ\n5z9YfM+5vyMLcP4AaB9DXbnu52oRMBcYHuvPinNR/4ZIPw/uezvzzVFnS4FT3cenAXOiqTdet5Sa\nsOmOa38MZ9jpM6p6rzjzNX6Bc5F2MzBFVa/zMExjjEkZIjIBOBWn624TcDfOH0mP4/wBvx+4UVXn\neBZjKiUaY4wxjU8qXaMxxhjTCFmiMcYYk1CWaIwxxiRUSiQaiWABSWOMMeklJRKNMcaYxitlEo2I\nNBORaeIsz75QREa4xzuLsyz9393luN8VkXyv4zXGGBOelBje7HadtQCaquoud9+J6UAP4HCcSWUD\nVXWeiLyEs/Dkv72L2BhjTLjC2tApSQT4g4icjLNuUAnQ1n3tS1X1bwQ0B2cWrDHGmDSQSonmSpzF\nLo9V1YPuQnpN3NeqAsrV4Oz9YowxJg2kzDUanD0TNrtJZhhOl5kxxpg053mLxl0auwpnUcrXRWQ2\nzuJ+Ua/8aowxJnV4PhhARPoCf1dn7w9jjDGNjKddZyJyA86S2b/2Mg5jjDGJ43mLxhhjTOOW1BaN\niHQUkQ/cCZiLReTH7vFiEZkqIivc+5bucRGRP4nIShFZICIDAup6W0R2iMgbyfwZjDHGRCbZXWfV\nwE9V9WhgKHCTiPQE7gCmqWoPnF0I73DLn40zabMHMBZ4MqCuh4CrkhW4McaY6CQ10ajqBlWd6z7e\njbPdaAkwAhjnFhuHs8c97vF/qWM60MK/77eqTgN2JzN+Y4wxkfNsMICIdAb6AzOAtqq6AZxkBLRx\ni5Xg7GfuV+4eM8YYkyY8STQi0gx4BbhVVXfVVzTIMRu9YIwxaSTpiUZEcnCSzHhVnege3uTvEnPv\nN7vHy4GOAW8vBdYnK1ZjjDGxS/aoMwGeBpaq6iMBL00GxriPxwCvBRz/njv6bCiw09/FZowxJj0k\ndR6NiJwIfAIsxFmhGeBXONdpXgI6AWXApapa4SamJ4CzgErgGlWd7db1CXAU0AzYBlyrqu8k7Ycx\nxhgTFpuwaYwxJqFSafVmY4wxjZAlGmOMMQllicYYY0xCWaIxxhiTUJZo/n97dSwAAAAAMMjfehT7\nSiIAVqIBYCUaAFaiAWAVhBCyIwlTUZMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xfd67128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.axes()\n",
    "\n",
    "df = pd.DataFrame({'P': P, 'Q': Q})\n",
    "df['P'].plot(ax=ax, secondary_y=True, label='rainfall')\n",
    "df['Q'].plot(ax=ax, label='Discharge')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}