example_semi_confined_agu2019.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Semi-confined flow between two polders\n",
    "Consider one-dimensional flow between two polders with different water levels, as shown in the figure below; the semi-confining layer at the left has a resistance $c_0$ and the semi-confining layer on the right has a resistance $c_1$. The two polders are separated by a dike. The width of the dike is neglected (the vertical black line in the figure). The water level in the polder on the left is $h_0^*$ while the water level in the polder on the right is $h_1^*$.\n",
    "The semi-confined aquifer extends to infinity on both sides. The transmissivity is the same on both sides and equal to $T$. \n",
    "The governing differential equations are\n",
    "\\begin{equation}\\begin{split}\n",
    "\\frac{\\text{d}^2(h-h_0^*)}{\\text{d}x^2} &= \\frac{h-h_0^*}{\\lambda_0^2}  \\hskip 2em  x\\le 0 \\\\\n",
    "\\frac{\\text{d}^2(h-h_1^*)}{\\text{d}x^2} &= \\frac{h-h_1^*}{\\lambda_1^2}  \\hskip 2em  x\\ge 0 \n",
    "\\end{split}\\end{equation}\n",
    "where $\\lambda_0=\\sqrt{c_0T}$ and $\\lambda_1=\\sqrt{c_1T}$. There are separate solutions for the left and right sides:\n",
    "\\begin{equation}\\begin{split}\n",
    "h=h_0^* +  A_0\\text{e}^{-x/\\lambda_0} + B_0\\text{e}^{x/\\lambda_0}  \\hskip 2em x\\le 0 \\\\\n",
    "h=h_1^* + A_1\\text{e}^{-x/\\lambda_1} + B_1\\text{e}^{x/\\lambda_1}  \\hskip 2em x\\ge 0 \\\\\n",
    "\\end{split}\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAADnCAYAAACTx2bHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAALJklEQVR4nO3df8yvdV3H8ddbDhiHyBUzDgUj0Wphcg66nFokLRmbKCpbfxSZGNG0UbMcmXOLnbKRra2z8g+LUMioWWCskLWdKMwochGHY3oaqyOu4zo0XCoRBsG7P743yvkVHzj3/b2u+z6Px1/n/nF2vf+4v/fz+7mu674+1d0BAP5/z5l6AABYDwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwSTpamqs6vqsqo6o6pOnnoegGdCMFlVVfWaqvrwEb58fJKfTvKmJP+1vKnWl6raUVU7pp4DONCmqQdgw9ma5J4jfO2MJB9KsjfJyUm+vKyh1pltUw8AHEowWW1bk+yvqk8k+fYkP9rdf5Ek3X1rVW3p7v2TTgjwLDgly2rbmuTB7j4vyU8lufSpXxRLYL0STFZNVR2f5JuS/PrKpzYl+eJ0EwGsHsFkNZ2d5N7ufmLl43OS/NOE8wCsGsFkNW1Ncu9TPj4nye6qOqmqbqiqa6vq0iP8X4BZE0xW09Yku5/y8XdnscK8JMlN3X1FkounGAzgaLlLllXT3e886OOzkqSqTk/yqZVPP77suQBWgxUmy7Avyekr//YzB6xLVpgsw0eTvL+qLkryZ1MPA/BsCCZrrrsfTvLWqecAOBpOjwHHFJsA8GwJJrDh2ASAtbC0U7Lf9gsf2xEPlYandcKpL9r2xGNf8Yv86NgEgFXnGiawEdkEgFVX3T31DMBTVNUdSdLd5087yfpVVfcmubG7f62qLkny+u524xlHxTVMYEOxCQBrRTCBjcYmAKwJwQQ2miNtAnBWVV1XVTdNNBfrnGACG81hNwHo7r3dfflEM7EBuEsW2FCOtAkAHC0rTAAYIJjAMaGqTqmqDyQ5t6rePfU8rD9OyQLHhO7+QpK3TT0H65cVJgAMEEwAGODReDAzx59yxr7nHP91X3/aZTt2TT0LrAO77v/Vi96xjANZYQLAACtMmBkPX4d5ssIEgAGCCQADBBMABggmAAwQTAAYIJgAMEAwAWCAYALAAMEEgAGCCQADBBMABggmAAzYtMyDbd++fX+SU5d5TFhvzjzzzCTJ9u3b7YwAT++Bq6++essyDrTsFaZYArCaltYVp2QBYIBgAsAAwQSAAYIJAAMEEwAGCCYADBBMABggmAAwQDABYIBgAsAAwQSAAYIJAAMEEwAGCCYADBBMABggmAAwQDABYIBgAsAAwQSAAYIJAAMEEwAGCCYADBBMABggmAAwQDCPcffff3/27Nkz9Rgwe/fdd1/27t079RhMaNnBfGDJx+NpbN68OTt37szNN9+cRx55ZOpxSLJly5Zs2bJl6jE4yIknnphbbrklt912Wx599NGpx+FrltaV6u5lHYuZqqrNSd6X5I1Jfry7d048EsxSVX1jkt9K8j1J3tzdn5x4JJZIMPmqqrogyXVJ/jTJu7r74YlHglmqqh/KIpy/k+SXu/uxiUdiCVzD5KtWVpbnJHleknuq6hUTjwSz1N1/nOTcJC9LcldVnT3xSCyBYHKA7v5id785ybuT3FJV762qE6aeC+amu/89yeuSfCDJx6vqZ6vK79QNzClZjqiqtiS5NsnpSX6suz818UgwS1X1wiQ3JHksyWXd/bmJR2INeDfEEXX3/iQXZ3Gt5i+r6qqqOm7isWB2uvtfk7w6yZ8n+YequqyqauKxWGVWmAypqhckuT6LN1lv6W5/kAaHUVXnJPlwks8m+cnu/o+JR2KVWGEypLs/m+QHkvxJkr+vqiu8g4ZDdffuJC9PsifJvVX1holHYpVYYfKMVdWLk/xekv1JfmLl5gfgIFX1vVlc2/xEknd095cmHomjYIXJM9bdn07yiiR3J9m18jdpwEG6+84k25L8TxarzfOnnYijYYXJUamql2ex2rw7yZXd/Z8TjwSzVFWvzeKu848keU93exblOmOFyVFZeTTYS5M8mGR3VV048UgwS919WxYPBjk9iztpXzbxSDxDVpismqr6wSQfTPKxJFd5tB4cauVmuR9OsiOLP9m6prv/d9qpGGGFyarp7tuzeAd9UhbXa1418UgwO73wB1mcmTkvyZ1V9Z0Tj8UAwWRVdfeXuvstSa5K8tGquqaqnjv1XDA33b0vyYVZ3EV7Z1Vd6dF68+aULGumqk5N8ttJXpDFVki7Jx4JZqmqviOLm+ceymKLvX+beCQOw7sZ1kx3P5DkTUl+I8ntVfUuj9aDQ3X3fUm+L8kdSe6uqks9GGR+rDBZiqo6M8mHkjw3i0fr/cvEI8EsVdVLs3i03meSvL27H5x4JFZYYbIUK7s3vCbJH2Wxf+DbvIOGQ3X3P2axz+bnsrh57qKJR2KFFSZLV1XflcX1mi8kuby7Pz/xSDBLVfXqLDY92Jnknd390LQTHdusMFm67t6T5FVJ/i7JPVV1+sQjwSx198eTbE1yXBaPoTx54pGOaVaYTGrl2ua+7n586llgzqpqa5Ld7Zf2ZAQTAAY4JcvsVNW3VNXNVXVPVf3zygPegYN4rSyXFSazUlWbstj55D3dfWtVbU5ynJsd4EBeK8tnhcncvDHJnu6+NUm6+7+TPFFVN1TVtVV16bTjwWwc7rXy/Kq6rqpumna0jUkwmZttSe466HOXJLmpu69IcvHyR4JZOuS10t17u/vyiebZ8ASTudmf5MVPflBVz89i/8Ann63pblpYONxrhTUkmMzN9UlOrapPV9WuJK9Msi+LaCZ+ZuFJ1+fQ1wpryE0/zF5VnZTk/Um+kuRvuvvGiUeCWaqqU5L8SpILkvxud18z8UgbimACwACntwBggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIABggkAAwQTAAZsWubBqqqXeTwANr7urmUcp7o1DACejlOyADBAMAFggGACwADBhJmpqh1VtWPqOYADLfUuWWDItqkHAA5lhQkAAwQTAAYIJgAMEEw4hlXVyVV11tRzwHogmHBse1GSu6rqD6vq7CN9U1X9VVVdsPLv91bVby5tQpgJd8nCPD2vql67pGO9PcnrkvxtVe1K8vPd/cmDvufqJL9UVd+c5NwkFy9pNpgNwYT5qSTfmuTKJR7zhCRfTnJekh9JckAwu/uvq6qS/FyS87v78SXOBrMgmDA/neQz3b3mK8yq+oYsVo9vTfL7Sd7X3Z8/zPe9JMlpSR7s7ofWei6YI9cw4dj2wiRPJHlJd//MEWJ5WpIbk7whycNVdeGSZ4RZsL0XzExV3ZEk3X3+tJMkVbU5ye1JfrG7d1bV92exCn3lxKPB0gkmzMycggl8jVOyADBAMAFggLtkYX52TT0AcCjXMAFggFOyADBAMAFggGACwADBBIABggkAAwQTAAYIJgAMEEwAGCCYADBAMAFggGACwADBBIAB/wfkmpvnQSuJDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "plt.plot([-0.1,  1.1], [0, 0], 'k', lw=1)\n",
    "plt.fill([-0.1, 1.1, 1.1, -0.1], [0.8, 0.8, 1, 1], color='grey')\n",
    "plt.plot([-0.1,  0.5], [1.2, 1.2], 'C0')\n",
    "plt.plot([0.5,  1.1], [1.1, 1.1], 'C0')\n",
    "plt.plot([0.5, 0.5], [0.8, 1.3], 'k')\n",
    "plt.text(0.25, 1.25, '$h^*_0$', ha='left')\n",
    "plt.text(0.75, 1.15, '$h^*_1$', ha='center')\n",
    "plt.ylim(-0.2, 1.4)\n",
    "plt.plot([0.5, 0.5], [-0.05, -0.15], 'k')\n",
    "plt.annotate(\"$x$\", va='center',\n",
    "            xy=(0.5, -0.1), \n",
    "            xytext=(0.6, -0.1), \n",
    "            arrowprops=dict(arrowstyle=\"<-\", shrinkA=0, shrinkB=0))\n",
    "plt.annotate(\"$c_0$\", va='center',\n",
    "            xy=(0.2, 0.8), \n",
    "            xytext=(0.3, 0.6), \n",
    "            arrowprops=dict(arrowstyle=\"-\", shrinkA=0, shrinkB=0))\n",
    "plt.annotate(\"$c_1$\", va='center',\n",
    "            xy=(0.8, 0.8), \n",
    "            xytext=(0.9, 0.6), \n",
    "            arrowprops=dict(arrowstyle=\"-\", shrinkA=0, shrinkB=0))\n",
    "plt.axis('off');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary conditions are that the head approaches $h_0^*$ far to the left and $h_1^*$ far to the right\n",
    "\\begin{equation}\\begin{split}\n",
    "&x\\to-\\infty \\hskip 2em h=h_0^* \\\\\n",
    "&x\\to+\\infty \\hskip 2em h=h_1^* \\\\\n",
    "\\end{split}\\end{equation}\n",
    "and continuity requires that the head and flow are continuous at $x=0$\n",
    "\\begin{equation}\\begin{split}\n",
    "&h(x=0^-)=h(x=0^+) \\\\\n",
    "&Q_x(x=0^-)=Q_x(x=0^+) \n",
    "\\end{split}\\end{equation}\n",
    "where $0^-$ and $0^+$ stand for just to the left and just to the right of zero, respectively. Application of the first two boundary conditions gives $A_0=0$ and $B_1=0$, after which application of the two continuity conditions gives the following two equations \n",
    "\\begin{equation}\\begin{split}\n",
    "h_0^*+B_0 &= h_1^* + A_1 \\\\\n",
    "-T\\frac{B_0}{\\lambda_0} &= T\\frac{A_1}{\\lambda_1}\n",
    "\\end{split}\\end{equation}\n",
    "Solving for the two remaining unknowns $A_1$ and $B_0$ gives\n",
    "\\begin{equation}\n",
    "A_1 = \\frac{\\lambda_1}{\\lambda_0+\\lambda_1}(h_0^*-h_1^*)  \\hskip 2em \n",
    "B_0 = -\\frac{\\lambda_0}{\\lambda_0+\\lambda_1}(h_0^*-h_1^*)\n",
    "\\end{equation}\n",
    "and hence\n",
    "\\begin{equation}\\begin{split}\n",
    "h&=h_0^*  -\\frac{\\lambda_0}{\\lambda_0+\\lambda_1}(h_0^*-h_1^*) \\text{e}^{x/\\lambda_0}  \\hskip 2em x\\le 0 \\\\\n",
    "h&=h_1^* + \\frac{\\lambda_1}{\\lambda_0+\\lambda_1}(h_0^*-h_1^*) \\text{e}^{-x/\\lambda_1} \\hskip 2em x\\ge 0 \\\\\n",
    "\\end{split}\\end{equation}\n",
    "Note that for the case of equal resistance values (so that $\\lambda_0=\\lambda_1$), the head in the aquifer at $x=0$ is $(h_0^*+h_1^*)/2$ as is to be expected from symmetry when the aquifer properties are the same on the left and right sides of the dike. \n",
    "\n",
    "The discharge vector is obtained as\n",
    "\\begin{equation}\n",
    "Q_x=-T\\frac{\\text{d}h}{\\text{d}x}\n",
    "\\end{equation}\n",
    "which gives\n",
    "\\begin{equation}\\begin{split}\n",
    "Q_x&=\\frac{T(h_0^*-h_1^*)}{\\lambda_0+\\lambda_1} \\text{e}^{x/\\lambda_0}  \\hskip 2em x\\le 0 \\\\\n",
    "Q_x&=\\frac{T(h_0^*-h_1^*)}{\\lambda_0+\\lambda_1} \\text{e}^{-x/\\lambda_1} \\hskip 2em x\\ge 0 \\\\\n",
    "\\end{split}\\end{equation}\n",
    "The total flow $Q_0$ from the polder on the left to the polder on the right is \n",
    "\\begin{equation}\n",
    "Q_0=Q_x(x=0)=\\frac{T(h_0^*-h_1^*)}{\\lambda_0+\\lambda_1}\n",
    "\\end{equation}\n",
    "Note again that 95% of the $Q_0$ has leaked from the polder into the aquifer over a distance of $3\\lambda_0$ to the left of the dike, and 95% of the $Q_0$ has discharged and into the polder on the right over a distance of $3\\lambda_1$. In other words, three times the leakage factor away from a disturbance, the head in the aquifer is approximately equal to the head above the semi-confining layer. \n",
    "\n",
    "An example is presented below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xd9X3/8ddHW7JsyZJl2bIsbwM2eIDMCGEEMDbTEAgjCXUIKRmkaZo0DfnRlIQmLZSmUFoyHGYSwggjOGyzCWBsGbzxkBcWtmV5L219fn/cI3otJFnr6ki67+fjcR733O/5nnveOlf2R2ebuyMiItJeCWEHEBGR3kkFREREOkQFREREOkQFREREOkQFREREOiQp7ADdadCgQT5y5MiwY4iI9CqLFi3a4e55TdvjqoCMHDmSkpKSsGOIiPQqZrapuXbtwhIRkQ5RARERkQ5RARERkQ5RARERkQ4JtYCY2X1mtt3Mlrcw3czsLjMrNbOlZnZ81LTZZrY2GGZ3X2oREYHwt0AeAGa2Mv08YFwwXA/8CsDMcoCbgZOAE4GbzWxgTJOKiMhhQi0g7v4msKuVLrOA33nEfCDbzIYCM4B57r7L3XcD82i9EImISBfr6deBDAM2R70vC9paao+JJ98vY9POQ6QkJZCcaCQnJpCcmEBKUgKpSQmkJCaQmpxAalIiqUnBa3JkWlpyYjBE+plZrGKKiHSrnl5Amvvf1ltp//QHmF1PZPcXRUVFHQrxzNKtvLpqe4fmjZZgkJ6cSHpKEhkpiVFD5H1mahL9UpPITEsiMzWJAWlJ9E9LZkB6EgPSkhmQnkx2euQ1LTmx03lERDqjpxeQMmB41PtCYEvQfmaT9teb+wB3nwPMASguLu7Q07Pu+8o0Ghqc2oYG6uqdmroGausbqK5roKa+gZq6yFBd10B1XT3VtQ1U1dVTVRt5X1lTT3VdA1W19RyqqaeyNtJ2sLqOyqBtx4FqDlTXRYaqOuoaWo+alpxATkYK2Rkp5PRLITczhdx+qeRmppCXmUpe/8gweEAquf1SSUzQlo+IdK2eXkDmAt82s0eIHDDf6+5bzexF4N+iDpyfC/wolkESEozUhERSk6BfaiyXBO5OdV0D+6vq2FdVG3mtrGVv1LDnUA27D9Wy+2ANuw7VsHnzIXYeqOFAdd2nPi8xwRjcP5UhWWkUZKdTkJVG4cAMCgemUzgwg6KcDNJTtEUjIu0TagExs4eJbEkMMrMyImdWJQO4+6+B54DzgVLgEHBtMG2Xmf0rsDD4qFvcvbWD8b2KmX1y7CSvf/uqVVVtZGumYn812/dXs31fFdv2VbFtbzXb9lWycss+Xl5ZTnVdw2Hz5Q9IZUROP0YN6sfovH6MyctkXH4mwwdmkKCtFxFphsXTM9GLi4tdN1OMbOHsOFBD2e5DbN5dyeZdh9i44yAbdx5kw46D7DhQ80nf9ORExg7O5Ogh/ZlQMICJBVlMKBhAZmpP33gVka5iZovcvbhpu/4XiENm9skxkqlFn758Zu+hWkorDrBu+wFWl+9n9bb9vLpqO39aVBbMD2PzMplUmM2UomyKRwxkfH5/HWcRiTMqIPIpWRnJnDBiICeM+L/i4u5s31/Nii17WVa2jyVle3hjzXaeeD9SVPqnJVE8YiCnjMnllNGDmFAwQAVFpI9TAZE2MTPyB6SRPyCNs47OByJFZfOuSko27WLhxt0s2LCT11ZXAJCVnsxp4wZx5lGDOX38IAb3TwszvojEgAqIdJiZUZSbQVFuBp8/vhCA8n1VvLtuJ38t3cEbayp4ZulWzGDq8GymTxjCjIn5jM7LDDm5iHQFHUSXmGlocFZu3cdrq7bz0spyln28F4CJBQO4aHIBF00uYFh2esgpReRIWjqIrgIi3ebjPZW8sHwbf1myhcWb92AGnxmTyxXFw5kxcYiurhfpoVRAUAHpST7aeYinPviYx9/fzOZdlQxIS+ILxcO55uQRjBzUL+x4IhJFBQQVkJ6oocF5b8Mu/rjgI55ftpW6BufMo/K4/rTRnDImVzefFOkBVEBQAenptu+r4uEFm/n9/I3sOFDD5MIsvnHGGGZMHKKr4UVCpAKCCkhvUVVbzxPvlzHnzfVs2nmIo4f05/vnHsU5xwzWFolICFRAUAHpbeobnL8s2cKdL69h485DTC7M4v+dfwwnjc4NO5pIXGmpgIT9SFuRFiUmGJdMHcbL3zuD/7h8EhX7q7lyzny+9dAiNu86FHY8kbinAiI9XlJiAlcUD+eV75/J96eP57VVFZz9ize48+U1VNfVhx1PJG6pgEivkZ6SyN+dPY7Xf3AmM48dwp0vr+W8/36L+et3hh1NJC6pgEivkz8gjbuunsoD106jtr6Bq+bM5+anl1NZo60Rke6kAiK91plHDeal757BV08dxYPvbuKC/3mLpWV7wo4lEjdCLSBmNtPMVptZqZnd2Mz0O8xscTCsMbM9UdPqo6bN7d7k0lOkpyTyLxdN4KGvnURlTT2X/vIdfvPGOuLp7EKRsIR2Gq+ZJQJrgOlAGZHH017t7itb6P93wFR3/2rw/oC7t+u2rjqNt2/bW1nLj55cynPLtnHuhHxu/8JkstKTw44l0uv1xNN4TwRK3X29u9cAjwCzWul/NfBwtySTXikrPZm7v3g8/3LhBF5dtZ2L/uevrNq2L+xYIn1WmAVkGLA56n1Z0PYpZjYCGAW8GtWcZmYlZjbfzC5paSFmdn3Qr6SioqIrcksPZmZ89bOjePTrJ1NdV89lv3yH11ZtDzuWSJ8UZgFp7p4ULe1Puwp43N2jT7MpCjapvgjcaWZjmpvR3ee4e7G7F+fl5XUusfQaJ4zI4ekbPsvIQf247sGFPPD2hrAjifQ5YRaQMmB41PtCYEsLfa+iye4rd98SvK4HXgemdn1E6c2GZKXxp2+cwjnH5POTv6zk58+u1MF1kS4UZgFZCIwzs1FmlkKkSHzqbCozOwoYCLwb1TbQzFKD8UHAqUCzB98lvmWkJPHrL5/AVz4zkt++tYEbn1hGfYOKiEhXCO2Z6O5eZ2bfBl4EEoH73H2Fmd0ClLh7YzG5GnjED//T8RjgN2bWQKQI3trS2VsiCQnGzRdNYEBaEne9WsqB6jruuHIKKUm6DEqkM3Q3Xokrv31zPT9/7kOmT8jnl186nuREFRGRI+mJp/GKdLu/PX00P714IvNWlvMPjy7W7iyRTghtF5ZIWGZ/ZiRVtfX8+/OrSElK4D8vn6wnHop0gAqIxKWvnzGGqtoG7nh5DQPSkrn5ogl62qFIO6mASNz6ztlj2VdVy71/3UDhwHS+dtrosCOJ9CoqIBK3zIybzj+GLXsq+flzH1KQnc75xw0NO5ZIr6GD6BLXEhKMO66cwtTh2Xz30cUs2rQr7EgivYYKiMS9tORE7pk9jaFZaXzjD+9Tvq8q7EgivYIKiAiQ0y+FOdcUc6Cqjm899D41dQ1hRxLp8VRARAJHDenP7V+YxKJNu/nXZ3RjA5EjUQERiXLhpAKuP300v5+/iScWlYUdR6RHUwERaeKfZhzFSaNy+PHTy9mw42DYcUR6LBUQkSaSEhO448opJCcm8N1HPqC2XsdDRJqjAiLSjILsdP7988expGwvd8xbE3YckR5JBUSkBecfN5Qrigv51RvreHfdzrDjiPQ4KiAirbj5oomMyMngn55YwqGaurDjiPQoKiAireiXmsRtl01i865KfvGSdmWJRAu1gJjZTDNbbWalZnZjM9O/YmYVZrY4GL4WNW22ma0Nhtndm1ziyUmjc/nyyUXc9/YGPvhod9hxRHqM0AqImSUCdwPnAROAq81sQjNdH3X3KcFwTzBvDnAzcBJwInCzmQ3spugSh34482iGDEjjh08s1VXqIoEwt0BOBErdfb271wCPALPaOO8MYJ6773L33cA8YGaMcorQPy2Zn196LGvKD/Cr19eFHUekRwizgAwDNke9LwvamrrMzJaa2eNmNryd82Jm15tZiZmVVFRUdEVuiVNnHZ3PhZOG8svXS9m861DYcURCF2YBae7xb00fUP0XYKS7TwJeBh5sx7yRRvc57l7s7sV5eXkdDisCcNMFx5Bgxr8992HYUURCF2YBKQOGR70vBLZEd3D3ne5eHbz9LXBCW+cViYWhWenc8LkxPL98G2+X7gg7jkiowiwgC4FxZjbKzFKAq4C50R3MLPrxcBcDjX/2vQica2YDg4Pn5wZtIjH3tdNGMzwnnZ/+ZYVucyJxLbQC4u51wLeJ/Mf/IfCYu68ws1vM7OKg23fMbIWZLQG+A3wlmHcX8K9EitBC4JagTSTm0pIT+fEFE1hTfoA/zN8UdhyR0Jh7s4cO+qTi4mIvKSkJO4b0Ae7ONfcuYMWWvbzxT59jQFpy2JFEYsbMFrl7cdN2XYku0gFmxg9nHs3uQ7Xc8+b6sOOIhEIFRKSDjivM4oJJQ7nnrxvYcaD6yDOI9DEqICKd8P3p46mua+B/Xy0NO4pIt1MBEemE0XmZXFFcyEPvbdLFhRJ3VEBEOuk7Z4/DzLjz5bVhRxHpViogIp00NCuda04ewZ8Xf8xHO7UVIvFDBUSkC/ztaaNJNOPXb+pGixI/VEBEusCQrDQuO6GQx0vKKN9XFXYckW6hAiLSRb55xhjq3fmtrguROKECItJFinIzuHhyAQ+99xG7DtaEHUck5lRARLrQt84cQ2VtPfe/vSHsKCIxpwIi0oXG5fdnxsR8fvfuJipr6sOOIxJTKiAiXey6z45mb2UtT35QFnYUkZhSARHpYtNGDmRiwQAeeHsj8XS3a4k/KiAiXczMuPbUUazdfoC/6qmF0oepgIjEwEWThzIoM4X7394YdhSRmAm1gJjZTDNbbWalZnZjM9O/Z2YrzWypmb1iZiOiptWb2eJgmNt0XpEwpSYl8sWTRvDqqu1s2HEw7DgiMRFaATGzROBu4DxgAnC1mU1o0u0DoNjdJwGPA/8RNa3S3acEw8WI9DBfPrmI5ETjwXc2hh1FJCbC3AI5ESh19/XuXgM8AsyK7uDur7l7493p5gOF3ZxRpMMG90/jwkkFPL6ojEM1dWHHEelyYRaQYcDmqPdlQVtLrgOej3qfZmYlZjbfzC5paSYzuz7oV1JRUdG5xCLtdPWJRRyoruPZpVvDjiLS5cIsINZMW7PnPJrZl4Fi4Pao5qLgIe9fBO40szHNzevuc9y92N2L8/LyOptZpF2mjRzI6Lx+PLpw85E7i/QyYRaQMmB41PtCYEvTTmZ2DnATcLG7f/LgaXffEryuB14HpsYyrEhHmBlXTRtOyabdrC3fH3YckS4VZgFZCIwzs1FmlgJcBRx2NpWZTQV+Q6R4bI9qH2hmqcH4IOBUYGW3JRdph88fX0hyomkrRPqc0AqIu9cB3wZeBD4EHnP3FWZ2i5k1nlV1O5AJ/KnJ6brHACVmtgR4DbjV3VVApEcalJnK9An5PPF+GdV1uj+W9B1JR+pgZmnAhcBpQAFQCSwHnnX3FZ1ZuLs/BzzXpO1fosbPaWG+d4DjOrNske505bQinlu2jXkry7lwUkHYcUS6RKtbIGb2E+Bt4BTgPSK7kx4D6oBbzWyemU2KdUiR3u60sYMYlp3OIwu0G0v6jiNtgSx095+0MO2/zGwwUNS1kUT6noQE4/ITCrnr1bVs21vFkKy0sCOJdFqrWyDu/uwRpm9395KujSTSN10ydRju8JclnzrZUKRXatNBdDMrNrOnzOz94L5Uy8xsaazDifQlowb1Y3JhFn9e/HHYUUS6xBEPogceAn4ALAMaYhdHpG+bNWUYtzyzktLt+xk7uH/YcUQ6pa2n8Va4+1x33+DumxqHmCYT6YMunDyUBIM/f6DdWNL7tbWA3Gxm95jZ1Wb2+cYhpslE+qDB/dM4dewgnl7ysZ5WKL1eWwvItcAUYCZwUTBcGKtQIn3ZJVOGsXlXJe9/tDvsKCKd0tZjIJPdXRfuiXSBGccO4aY/L+PPH2zhhBE5YccR6bC2boHMb+ZhTyLSAZmpSZxzTD7PLttKXb3OSZHeq60F5LPA4uDxszqNV6STLpw0lF0Ha1iwcVfYUUQ6rK27sGbGNIVInDl9fB5pyQm8uHwbnxkzKOw4Ih1ypHthZQJEn7rb9DTexj4i0nYZKUmcMT6PF1eU09Cgs7GkdzrSLqynzewXZna6mfVrbDSz0WZ2nZm9iLZORDrkvGOHsm1fFUvK9oQdRaRDjnQvrLOBV4CvAyvMbK+Z7QT+AAwBZrv747GPKdL3fO7owSQnGi8s3xZ2FJEOOeIxkOae2SEinZeVnsxnxgzihRXbuPG8ozGzsCOJtEuYj7TFzGYGZ3aVmtmNzUxPNbNHg+nvmdnIqGk/CtpXm9mM7swt0lVmHjuETTsPsWqbnpcuvU9oBcTMEoG7gfOACcDVzVxrch2w293HAncAtwXzTiDyDPWJRI7B/DL4PJFeZfqEfBIM7caSXinMLZATgVJ3X+/uNcAjwKwmfWYBDwbjjwNnW2Q7fxbwiLtXu/sGoDT4PJFeZVBmKtNG5qiASK90pNN4c1obOrnsYUD08z3LgrZm+7h7HbAXyG3jvI0/w/VmVmJmJRUVFZ2MLNL1Zh47hNXl+1lfcSDsKCLtcqQtkEVASfBaAawB1gbjizq57OaOGDY9Ib6lPm2ZN9LoPsfdi929OC8vr50RRWJvxsQhALy4ojzkJCLtc6TTeEe5+2jgReAidx/k7rlE7sT7ZCeXXQYMj3pfCDR9SMInfcwsCcgCdrVxXpFeoSA7ncnDs3lh+dawo4i0S1uPgUwLTucFwN2fB87o5LIXAuPMbJSZpRA5KD63SZ+5wOxg/HLgVY88RGEucFVwltYoYBywoJN5REIzc+IQlpTt5eM9lWFHEWmzthaQHWb2z2Y20sxGmNlNwM7OLDg4pvFtIls3HwKPufsKM7vFzC4Out0L5JpZKfA94MZg3hXAY8BK4AXgBnev70wekTDNmJgPwIs6mC69iLXlqWjBAfObgdODpjeBn7p7r7qVaHFxsZeUlIQdQ6RZM+54k6yMZB77+ilhRxE5jJktcvfipu1tuhtvUCj+vstTicgnZh47hLteXUvF/mry+qeGHUfkiNq0C8vM8szsdjN7zsxebRxiHU4knsw8dgjuMG+lzsaS3qGtx0AeAlYBo4CfAhuJHAQXkS5y9JD+jMjN4IUVOg4ivUNbC0iuu98L1Lr7G+7+VeDkGOYSiTtmxsxjh/BO6Q52H6wJO47IEbW1gNQGr1vN7AIzm0rk2gsR6UIXTSqgrsF5TteESC/Q1gLyMzPLAr4P/CNwD/APMUslEqcmFgxg7OBMnv5A18VKz9fWs7CeCUb3Ap+LXRyR+GZmXDKlgP98aQ1luw9RODAj7EgiLWrrWVjjzewVM1sevJ9kZv8c22gi8eniyZH7gs5doq0Q6dnaugvrt8CPCI6FuPtSIrceEZEuVpSbwfFF2dqNJT1eWwtIhrs3vddUXVeHEZGIS6YOY3X5fj7cui/sKCItas+9sMYQ3DLdzC4HdJqISIxccNxQEhOMPy/+OOwoIi1qawG5AfgNcLSZfQx8F/hmzFKJxLnczFROHzeIuYu30NBw5PvViYShTQUkeOzsOUAecLS7f9bdN8Y0mUicu/T4QrbureLtdTvCjiLSrDadxmtmqcBlwEggKfJYcnD3W2KWTCTOnTshn+yMZB5ZuJnTxulpmtLztHUX1tPALCIHzg9GDSISI2nJiVw6dRgvrdjGzgPVYccR+ZQ2bYEAhe4+M6ZJRORTrj6xiPvf3shTH3zM104bHXYckcO0dQvkHTM7rqsWamY5ZjbPzNYGrwOb6TPFzN41sxVmttTMroya9oCZbTCzxcEwpauyifQk4/P7c3xRNg8v+Ii2PPxNpDu1WkDMbJmZLQU+C7xvZquD/8wb2zvqRuAVdx8HvBK8b+oQ8DfuPhGYCdxpZtlR03/g7lOCYXEnsoj0aFdNK2JdxUEWbdoddhSRwxxpF9aFMVruLODMYPxB4HXgh9Ed3H1N1PgWM9tO5CywPTHKJNIjXTBpKLc8s5KHF2ymeGRO2HFEPtHqFoi7b2pt6MRy8919a7CMrcDg1jqb2YlACrAuqvnnwdbQHcFZYi3Ne72ZlZhZSUVFRScii4SjX2oSF00u4NllW9hbWXvkGUS6SVuPgbSbmb1sZsubGWa183OGAr8HrnX3hqD5R8DRwDQghyZbL9HcfY67F7t7cV6eToWU3ulLJxVRVdvAYws3hx1F5BNtPQur3YILD5tlZuVmNtTdtwYFYnsL/QYAzwL/7O7zoz678TYq1WZ2P5FnlIj0WccOy+LEkTk8+O5GvvrZUSQmWNiRRGK3BXIEc4HZwfhsIteZHMbMUoCngN+5+5+aTBsavBpwCbA8pmlFeoBrTx1J2e5K5q0sDzuKCBBeAbkVmG5ma4HpwXvMrNjM7gn6XAGcDnylmdN1HzKzZcAyYBDws+6NL9L9pk/IZ1h2Ove/vSHsKCJADHdhtcbddwJnN9NeAnwtGP8D8IcW5j8rpgFFeqCkxARmf2YE//bcKlZs2cvEgqywI0mcC2sLREQ64MriItKTE3ng7Y1hRxFRARHpTbIykrnshGE8vWQLFft1fywJlwqISC/z1VNHUVvfwH06FiIhUwER6WVG52Vy/nFD+f27m3RhoYRKBUSkF7rhzLEcqK7jd+9sDDuKxDEVEJFeaELBAM4+ejD3vb2Bg9V1YceROKUCItJLfetzY9l9qJaHF3wUdhSJUyogIr3UCSMGcsroXOa8uZ6q2vqw40gcUgER6cX+7qyxbN9fzR/f01aIdD8VEJFe7DNjB3Hq2Fzufq2UAzoWIt1MBUSkl/vBjKPZebCGe9/SdSHSvVRARHq5KcOzmTlxCL99az27DtaEHUfiiAqISB/wjzPGc6imjl++Vhp2FIkjKiAifcDYwf257PhCfjd/E2W7D4UdR+KECohIH/EP08eTYPDvz68KO4rECRUQkT6iIDudb54xlmeXbuXddTvDjiNxIJQCYmY5ZjbPzNYGrwNb6Fcf9TTCuVHto8zsvWD+R4PH34rEva+fMZph2en89C8rqKtvCDuO9HFhbYHcCLzi7uOAV4L3zal09ynBcHFU+23AHcH8u4HrYhtXpHdIS07kxxcew6pt+/mjbnEiMRZWAZkFPBiMPwhc0tYZzcyAs4DHOzK/SF83Y+IQTh2byy9eWqPTeiWmwiog+e6+FSB4HdxCvzQzKzGz+WbWWCRygT3u3njZbRkwrKUFmdn1wWeUVFRUdFV+kR7LzLj5ookcrK7jZ8+sDDuO9GExKyBm9rKZLW9mmNWOjyly92Lgi8CdZjYGsGb6eUsf4O5z3L3Y3Yvz8vLa+VOI9E7j8/vzrTPH8OQHH/Pa6u1hx5E+KmYFxN3PcfdjmxmeBsrNbChA8Nrsb7i7bwle1wOvA1OBHUC2mSUF3QqBLbH6OUR6qxvOGsvYwZnc9OQy3SdLYiKsXVhzgdnB+Gzg6aYdzGygmaUG44OAU4GV7u7Aa8Dlrc0vEu9SkxK57bJJbN1Xxe0v6NoQ6XphFZBbgelmthaYHrzHzIrN7J6gzzFAiZktIVIwbnX3xh26PwS+Z2alRI6J3Nut6UV6iRNGDGT2KSN58N1NzF+va0Oka1nkD/r4UFxc7CUlJWHHEOlWB6vruOCut6ipa+D5vz+drIzksCNJL2Nmi4Lj0YfRlegifVy/1CTuvGoq2/dX8//+vIx4+qNRYksFRCQOTBmezT9MH8+zS7fy+KKysONIH6ECIhInvnHGGE4alcPNc1ewYcfBsONIH6ACIhInEhOMO66cQkpSAt/4/SIO6tRe6SQVEJE4UpCdzv9cPZW12/fzwyeW6niIdIoKiEicOW1cHv844yieWbqVe/+q56hLx6mAiMShb54xhvOOHcK/P7+Kv67dEXYc6aVUQETikJlx+xcmMzYvk2/+YRGrt+0PO5L0QiogInEqMzWJ+6+dRkZqIl+5fwHb9laFHUl6GRUQkThWkJ3OfV+Zxr7KWq59YKFuuijtogIiEucmFmRx95eOZ035fr724EIqa+rDjiS9hAqIiHDmUYP5rysm896GXVz/+xKq61RE5MhUQEQEgFlThnHb5yfx1tod3PDQB9TWN4QdSXo4FRAR+cQV04bzr7Mm8vKH5XzzD+9TVastEWmZCoiIHOaaU0Z+UkSuvV8H1qVlKiAi8inXnDKSO66czIKNu/jSPe+x+2BN2JGkBwqlgJhZjpnNM7O1wevAZvp8zswWRw1VZnZJMO0BM9sQNW1K9/8UIn3bpVML+fWXT+DDrfv4/K/e0R185VPC2gK5EXjF3ccBrwTvD+Pur7n7FHefApwFHAJeiuryg8bp7r64W1KLxJnpE/J56Gsnsbeylkt/+bYeiyuHCauAzAIeDMYfBC45Qv/Lgefd/VBMU4nIp0wbmcNT3/oMuf1SuObe93h4wUe6i68A4RWQfHffChC8Dj5C/6uAh5u0/dzMlprZHWaW2tKMZna9mZWYWUlFRUXnUovEqRG5/XjyW6dy8uhcfvTkMn7w+FJdcChYrP6SMLOXgSHNTLoJeNDds6P67nb3Tx0HCaYNBZYCBe5eG9W2DUgB5gDr3P2WI2UqLi72kpKSdv8sIhJR3+D89ytr+Z9X13JUfn9++aXjGZ2XGXYsiTEzW+TuxU3bY7YF4u7nuPuxzQxPA+VBEWgsBttb+agrgKcai0fw2Vs9ohq4HzgxVj+HiPyfxATje9PHc/9XplG+r4oL7vorD723Sbu04lRYu7DmArOD8dnA0630vZomu6+iio8ROX6yPAYZRaQFZx41mOf//nSKRw7kpqeW89UHFrJ9n+7mG2/CKiC3AtPNbC0wPXiPmRWb2T2NncxsJDAceKPJ/A+Z2TJgGTAI+Fk3ZBaRKEOy0njw2hP56cUTeWfdTs75rzf443sf0dCgrZF4EbNjID2RjoGIxMb6igPc9NRy3l2/k+IRA/n5pcdx1JD+YceSLtLtx0BEJH6Mzsvkj397Ev/5hcmsqzjA+Xe9xY//vJxduoK9T1MBEZEuYWZcfkIhr37/TK45eQR/XPARZ9z+Gr95Y51O+e2jtAtLRM8eYWsAAAzNSURBVGKidPt+fvbsh7y+uoK8/qnccOYYrj6piNSkxLCjSTu1tAtLBUREYmrBhl3850urWbBhF/kDUrnus6O4+sQi+qclhx1N2kgFBBUQkbC4O2+X7uTu10p5d/1O+qcl8aWTRvDlk4soHJgRdjw5AhUQVEBEeoIlm/fwmzfX8cLybQCcc0w+15wyglPHDCIhwUJOJ81RAUEFRKQn+XhPJQ/N38QjCzez62ANw7LTueyEQi4/vpCiXG2V9CQqIKiAiPREVbX1zFtZzp8WlfHW2grcYcrwbC6aXMCFk4aSPyAt7IhxTwUEFRCRnm7LnkrmLtnC3MVbWLl1HxApJudOzOfcCfmMycskcgcj6U4qIKiAiPQmpdsP8MLyrby0spylZXsBGJadzhlH5XHG+DxOHpVLVobO5OoOKiCogIj0Vlv3VvLqqu28sbqCt0t3cLCmHjM4tiCLk0fnUDwyh+IRA8nNbPHRQNIJKiCogIj0BTV1DXzw0W7eXb+Td9ft5IOP9lBT3wDAqEH9mFyYxeTh2UwqzOaYof3JSEkKOXHvpwKCCohIX1RVW8/yj/dSsmk3izbtZmnZHsr3VQOQYJGiMqEgi6PyMxmf35/x+f0ZnpNBok4ZbrOWCohKs4j0amnJiZFdWCNzPmkr31fFks17WLl1Hyu27OP9Tbv5y5Itn0xPSUpgVG4/xgzux8jcyFCUm0FRTgb5A9JUXNpIBURE+pz8AWmcO3EI5078v6dqH6iuY235ftaWH2BdRWT4cOt+XlpRTl3UM0ySEoyC7HQKstMoyEqnIDud/Kw0hgyIDIMHpJLbL4WkRN2LNpQCYmZfAH4CHAOc6O7N7lcys5nAfwOJwD3u3vjgqVHAI0AO8D5wjbvrvtEi0qLM1CSmFg1katHAw9rr6hvYureKjTsPUra7ks27DrF5dyVb91Qyf/1OyvdXU9/kIVlmkJORQm5mCrn9UsnNTCGnXwoDMyKv2RnJZKUnk52RQlZ6MgPSkhiQnkxyHys6YW2BLAc+D/ympQ5mlgjcTeSJhWXAQjOb6+4rgduAO9z9ETP7NXAd8KvYxxaRviYpMYHhORkMz2n+6vf6BmfHgWrK91WxbW8VFQeqqdhfzfb91ew6UMPOg9Ws2LKPXQdr2FtZ2+qy0pIT6J+WTP/UJDLTkuiXkkS/1CQyUxPplxoZT09OJCMlMqQlJ5KekkhaUuQ1NSmBtORE0pITSE1KJLXxNSmBlMSEbr8VTCgFxN0/BI50QdCJQKm7rw/6PgLMMrMPgbOALwb9HiSyNaMCIiJdLjHByB+QRv6ANCYVtt63rr6BPZW17DlUy97KWvYcqmFfVS37KuvYW1nLgeo69lfVsq+qjgNVdRysrqNs9yEO1tRRWVPPwep6Kms7/uyUpAQjJSmB5MQEUoKikpxoJCUmcO/sYkbk9uvwZze7vC79tK41DNgc9b4MOAnIBfa4e11U+7BuziYi8ilJiQkMykxlUCeuR2locKrq6jlUU09lTT1VtZGiUllTT1VdA1W1kbaaugaqgyEyHmmrrW+gtt6p/mQ8MqQld/1zWGJWQMzsZWBIM5Nucven2/IRzbR5K+0t5bgeuB6gqKioDYsVEQlPQoKRkZLUK65fiVlCdz+nkx9RBgyPel8IbAF2ANlmlhRshTS2t5RjDjAHIteBdDKTiIgEevIpAQuBcWY2ysxSgKuAuR658vE14PKg32ygLVs0IiLShUIpIGZ2qZmVAacAz5rZi0F7gZk9BxBsXXwbeBH4EHjM3VcEH/FD4HtmVkrkmMi93f0ziIjEO93KREREWtXSrUx68i4sERHpwVRARESkQ1RARESkQ1RARESkQ+LqILqZVQCbOjj7ICLXoPQ0ytU+ytU+ytU+fTXXCHfPa9oYVwWkM8yspLmzEMKmXO2jXO2jXO0Tb7m0C0tERDpEBURERDpEBaTt5oQdoAXK1T7K1T7K1T5xlUvHQEREpEO0BSIiIh2iAiIiIh2iAgKY2RfMbIWZNZhZcZNpPzKzUjNbbWYzotpnBm2lZnZjVPsoM3vPzNaa2aPBrei7KuejZrY4GDaa2eKgfaSZVUZN+3XUPCeY2bIg5112hOcIdzDXT8zs46jlnx81rV3rr4tz3W5mq8xsqZk9ZWbZQXuo66uZnDFfF60se7iZvWZmHwb/Bv4+aG/3dxqDbBuD72KxmZUEbTlmNi/49zXPzAYG7RZ8X6XB9318jDIdFbVOFpvZPjP7bhjry8zuM7PtZrY8qq3d68fMZgf915rZ7HaFcPe4H4BjgKOA14HiqPYJwBIgFRgFrAMSg2EdMBpICfpMCOZ5DLgqGP818M0YZf4F8C/B+EhgeQv9FhC5bb4BzwPnxSDLT4B/bKa93euvi3OdCyQF47cBt/WE9dVked2yLlpZ/lDg+GC8P7Am+N7a9Z3GKNtGYFCTtv8AbgzGb4z6Ts8Pvi8DTgbe64Z1lwhsA0aEsb6A04Hjo3+X27t+gBxgffA6MBgf2NYM2gIB3P1Dd1/dzKRZwCPuXu3uG4BS4MRgKHX39e5eAzwCzAr+Wj0LeDyY/0Hgkq7OGyznCuDhI/QbCgxw93c98tvyu1jkaUW71l9XL9zdX/LIc2UA5hN5emWLQlpf3bIuWuLuW939/WB8P5Fn7wxrZZaWvtPuMovIvys4/N/XLOB3HjGfyFNLh8Y4y9nAOndv7e4WMVtf7v4msKuZ5bVn/cwA5rn7LnffDcwDZrY1gwpI64YBm6PelwVtLbXnAnui/tNqbO9qpwHl7r42qm2UmX1gZm+Y2WlR+cuayRkL3w42je9r3Gym/esvlr5K5C+wRmGvr0ZhrItmmdlIYCrwXtDUnu80Fhx4ycwWmdn1QVu+u2+FSPEDBoeQq9FVHP5HXNjrC9q/fjqVL24KiJm9bGbLmxla+2uvuf3f3oH2rs55NYf/4m4Fitx9KvA94I9mNqAr8rQx16+AMcCUIMsvGmdrYfndlauxz01AHfBQ0BTz9dWeHyGEZX46hFkm8ATwXXffR/u/01g41d2PB84DbjCz01vp263r0SLHNi8G/hQ09YT11ZqY/FtM6nCcXsbdz+nAbGXA8Kj3hcCWYLy59h1ENg2Tgq2Q6P5dktPMkoDPAydEzVMNVAfji8xsHTA+yB+926bdedqaKyrfb4FngrftXX9dnis4KHghcHawW6pb1lc7tLaOuoWZJRMpHg+5+5MA7l4eNb2t32mXcvctwet2M3uKyK6fcjMb6u5bg10w27s7V+A84P3G9dQT1legveunDDizSfvrbV1Y3GyBdNBc4CozSzWzUcA4IgdZFwLjLHLGVQqRTdm5wX9QrwGXB/PPBp7u4kznAKvc/ZNdLWaWZ2aJwfjoIOf6YBN2v5mdHBw3+ZsY5Gk8dtDoUqDxrJB2rb8Y5JoJ/BC42N0PRbWHur6a6JZ10ZLg57wX+NDd/yuqvb3faVfn6mdm/RvHiZwQsTxYfuOZQtH/vuYCfxOcbXQysLdxV06MHLYXIOz1FaW96+dF4FwzGxjsdjs3aGubrjgboLcPRL7wMiJ/lZYDL0ZNu4nImROriTojh8hZDWuCaTdFtY8m8gtSSmTzNrWLsz4AfKNJ22XACiJne7wPXBQ1rZjIL/M64H8J7j7QxZl+DywDlga/qEM7uv66OFcpkf27i4Ph1z1hfTWTM+bropVlf5bILoulUevp/I58p12ca3Tw/SwJvqubgvZc4BVgbfCaE7QbcHeQaxlRZ1PGIFsGsBPI6sy/gS7I8TCR3WW1RP7/uq4j64fI8cHSYLi2PRl0KxMREekQ7cISEZEOUQEREZEOUQEREZEOUQEREZEOUQEREZEOUQERCYmZpQe3UklsxzzfNrNrY5lLpK10Gq9ISMzsBiJ3C/7vdsyTAbztkduwiIRKWyAiXczMpgU31UsLrqheYWbHNtP1SwRXCpvZmcHWyGNmtsbMbjWzL5nZAos8E2MMgEeuqN9oZt15B1yRZsXNvbBEuou7LzSzucDPgHTgD+6+PLpPcNuS0e6+Map5MpFn0+wi8lyGe9z9RIs85OnvgO8G/UqI3JE5lrfEEDkiFRCR2LiFyH2uqoDvNDN9ELCnSdtCD+7fFNzg8aWgfRnwuah+24GjuzStSAdoF5ZIbOQAmUSe8pfWzPTKZtqro8Ybot43cPgfe2nB/CKhUgERiY05wI+JPIPktqYTPfL0t0Qza664HMl4/u9uryKhUQER6WJm9jdAnbv/EbgVmGZmZzXT9SUid8Ntr1OBlzsRUaRL6DRekZCY2VTge+5+TSznEYkVbYGIhMTdPwBea8+FhEQOvv84RpFE2kVbICIi0iHaAhERkQ5RARERkQ5RARERkQ5RARERkQ5RARERkQ75//j/MKJL5FBBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# parameters\n",
    "h0star = 1\n",
    "h1star = -1\n",
    "k = 10\n",
    "H = 10\n",
    "c0 = 400\n",
    "c1 = 400\n",
    "# solution h\n",
    "T = k * H\n",
    "lab0 = np.sqrt(c0 * T)\n",
    "lab1 = np.sqrt(c1 * T)\n",
    "x = np.linspace(-1000, 1000, 101)\n",
    "h = np.zeros_like(x)\n",
    "h[x < 0]  = h0star - lab0 / (lab0 + lab1) * \\\n",
    "            (h0star - h1star) * np.exp(x[x < 0] / lab0)\n",
    "h[x >= 0] = h1star + lab1 / (lab0 + lab1) * \\\n",
    "            (h0star - h1star) * np.exp(-x[x >= 0] / lab1)\n",
    "# basic plot head\n",
    "plt.plot(x, h)\n",
    "plt.xlabel('x (m)')\n",
    "plt.ylabel('head (m)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": [
     "hide_input"
    ]
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8d1fd9a540ad4f63a082c1177e1694d8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=100, description='c0', max=600, min=50, step=50), IntSlider(value=400, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def makeplot(c0=100, c1=400):\n",
    "    # parameters\n",
    "    h0star = 1\n",
    "    h1star = -1\n",
    "    k = 10\n",
    "    H = 10\n",
    "    # solution h, Qx\n",
    "    T = k * H\n",
    "    lab0 = np.sqrt(c0 * T)\n",
    "    lab1 = np.sqrt(c1 * T)\n",
    "    x = np.linspace(-1000, 1000, 101)\n",
    "    h = np.zeros_like(x)\n",
    "    h[x < 0]  = h0star - lab0 / (lab0 + lab1) * (h0star - h1star) * np.exp(x[x < 0] / lab0)\n",
    "    h[x >= 0] = h1star + lab1 / (lab0 + lab1) * (h0star - h1star) * np.exp(-x[x >= 0] / lab1)\n",
    "    Qx = np.zeros_like(x)\n",
    "    Qx[x<0] = T * (h0star - h1star) / (lab0 + lab1) * np.exp(x[x<0] / lab0)\n",
    "    Qx[x>=0] = T * (h0star - h1star) / (lab0 + lab1) * np.exp(-x[x>=0] / lab1)\n",
    "    # solution psi\n",
    "    psi = np.zeros((2, len(x)))\n",
    "    psi[1] = -Qx\n",
    "    xg = x\n",
    "    zg = [0, H]\n",
    "    # plot\n",
    "    plt.figure(figsize=(10, 8))\n",
    "    plt.subplot(211)\n",
    "    plt.xlim(-1000, 1000)\n",
    "    plt.ylabel('head')\n",
    "    plt.plot(x, h)\n",
    "    plt.axvline(-3 * lab0, linestyle=':', color='C2')\n",
    "    plt.axvline(3 * lab1, linestyle=':', color='C2')\n",
    "    plt.axvline(0, linestyle=':', color='C1')\n",
    "    plt.text(-3 * lab0, -1, r'$-3\\lambda_0$', ha='left', va='bottom')\n",
    "    plt.text(3 * lab1, -1, r'$3\\lambda_1$', ha='right', va='bottom')\n",
    "    plt.subplot(212)\n",
    "    plt.title('streamlines')\n",
    "    plt.contour(xg, zg, psi, 20, colors='C0', linestyles='-')\n",
    "    plt.axvline(-3 * lab0, linestyle=':', color='C2')\n",
    "    plt.axvline(3 * lab1, linestyle=':', color='C2')\n",
    "    plt.text(-3 * lab0, 0.2, r'$-3\\lambda_0$', ha='left', va='bottom')\n",
    "    plt.text(3 * lab1, 0.2, r'$3\\lambda_1$', ha='right', va='bottom')\n",
    "    plt.xlim(-1000, 1000)\n",
    "    plt.xlabel('x (m)')\n",
    "    plt.ylabel('z (m)')\n",
    "    \n",
    "interact(makeplot, c0=(50, 600, 50), c1=(50, 600, 50));"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}