MWW11/MWW11_Cooper-Jacob_1.2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Developed by: Fritz.Kalwa@tu-dresden.de / Thomas.Reimann@tu-dresden.de\n",
    "\n",
    "<br>Last change: 2020 11 09\n",
    "\n",
    "<br>Current state: funktional und nett\n",
    "\n",
    "<br>ToDo:\n",
    "\n",
    "- Inhalt kurz und Lernziele ergänzen (FK)\n",
    "- Inhalte gliedern (Überschriften) FK\n",
    "- Optimieren(?)\n",
    "- Output hübscher (?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "    \n",
    "### MWW11 - Hydrogeologische und Hydrochemische Methoden\n",
    "### Thema: Pumpversuchsauswertung\n",
    "# Pumpversuchsauswertung nach Cooper-Jacob  </center>\n",
    "_Mit diesem Notebook soll Ihnen die Pumpversuchsauswertung nach dem Ansatz von Cooper-Jacob nähergebracht werden. Zum Schluss sollten Sie vor allem..._\n",
    "- ...die grundlegenden Zusammenhänge zwischen den gesuchten Parametern S & T und den bekannten Parametern Q & r in Abhängigkeit von der Zeit t veranschaulichen können.\n",
    "- ...die Limitierungen des Ansatzes von Cooper-Jacob benennen können.\n",
    "- ...die Reaktionen an unterschiedlichen Messstellen miteinander vergleichen können.\n",
    "- ...Transmissivität und Speicherkoeffizient bei gegebenen Messwerten bestimmen können (iterativ __und__ rechnerisch), unter Berücksichtigung des des kritischen Zeitpunktes $t_{min}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notwendige Funktionen / Bibliotheken einladen\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.special\n",
    "import os\n",
    "import math\n",
    "from ipywidgets import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Die vereinfachte Brunnenfunktion\n",
    "Für eine Pumpversuchsauswertung entwickelt Charles V. Theis 1935 einen Lösungsansatz. Hierbei wurde die radiale Grundwasserströmung auf ein Zentrum mit der Koordinate *0* hin analytisch berechnet. Dabei wird die Absenkung bei der Pumprate *Q* in einer Messstelle mit dem Abstand zum Zentrum *r* in Abhängigkeit von der Pumpzeit *t* bei bekannter Transmissivität *T* bestimmt:\n",
    "\n",
    "<br>\n",
    "<center>$\\large s(r,t)=\\frac{Q}{4\\pi T}W(u) $</center>\n",
    "<br>\n",
    "\n",
    "*W(u)* ist die Brunnenfunktion:\n",
    "\n",
    "<br>\n",
    "<center>$\\large W(u) = \\int_{u }^{+\\infty} \\frac{e^{-\\tilde u}}{\\tilde u}d\\tilde u = -0,5772 - ln(u) + u - \\frac{u^2}{2*2!} + \\frac{u^3}{3*3!} - \\frac{u^4}{4*4!}+...$</center>\n",
    "<br>\n",
    "\n",
    "mit der dimensionslosen Variable *u*:\n",
    "\n",
    "<br>\n",
    "<center>$\\large u = \\frac{Sr^2}{4Tt} $</center>\n",
    "<br>\n",
    "\n",
    "wobei *S* den Speicherkoeffizienten darstellt. Die Vereinfachung von Cooper-Jacob geht darauf zurück, dass *u* mit steigender Zeit *t* immer mehr gegen 0 läuft. Deshalb kann die Brunnenfunktion *W(u)* ab einem gewissen Zeitpunkt vereinfacht werden:\n",
    "\n",
    "<br>\n",
    "<center>$\\large W(u) \\approx -0,5772 - ln(u)$</center>\n",
    "<br>\n",
    "\n",
    "Diese vereinfachte Funktion nimmt bei halblogarithmischer Darstellung die Form einer Geraden an, wir sprechen im Folgenden deshalb von der \"logarithmischen Geraden\". Die Theis'sche Funktion *W(u)*, die logarithmische Gerade und der Parameter *u* in Abängigkeit von der Zeit sind in der folgenden Grafik geplottet. Sehr deutlich wird dabei, dass die Cooper-Jacob-Vereinfachung sich erst bei ausreichend großen Zeiten der Theis'schen Brunnenfunktion annähert."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "12773d0bb3264640a6e3acaeb992b0de",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=10.0, continuous_update=False, description='r:', min=1.0, step=1.0), F…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.uplot(r, T, S)>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#PLOT u as function of t\n",
    "\n",
    "def uplot(r, T, S):\n",
    "    tmax = 100000\n",
    "    u_max = 1\n",
    "    r_max = 100000\n",
    "    u = [u_max for x in range(r_max)]\n",
    "    t = [u_max for x in range(r_max)]\n",
    "\n",
    "    #W(u) as full computaton W and simplified computation W2\n",
    "    W = [u_max for x in range(r_max)]\n",
    "    W2 = [u_max for x in range(r_max)]\n",
    "    for x in range(1,tmax,1):\n",
    "            u[x] = r**2*S/4/T/x\n",
    "            t[x] = x\n",
    "            W[x] = scipy.special.exp1(u[x])\n",
    "            W2[x] = -0.577216-np.log(u[x])\n",
    "            uu = r**2*S/4/T/0.02\n",
    "\n",
    "    #PLOT\n",
    "    fig = plt.figure(figsize=(9,6))\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(t, u, 'green')\n",
    "    ax.plot(t, W,'brown')\n",
    "    ax.plot(t, W2,'blue', ls=':')\n",
    "    3\n",
    "    #plt.yscale(\"log\")\n",
    "    plt.xscale(\"log\")\n",
    "    plt.axis([1,1E5,1E-1,1E+1])\n",
    "    ax.set(xlabel='t', ylabel='u bzw W(u)')\n",
    "    plt.title('Die Brunnenfunktion nach Theis und ihre Vereinfachung',fontweight='bold', fontsize=16)\n",
    "    ax.grid(which=\"both\")\n",
    "    plt.legend(('u','Brunnenfunktion nach Theis','vereinfachte Brunnenfunktion'),loc=0)\n",
    "    plt.axvline(x=uu, ls='--', color='red')\n",
    "    plt.show()\n",
    "    \n",
    "# Computation\n",
    "interact(uplot,\n",
    "         r = widgets.FloatSlider(value=10, min=1, max=100, step=1, description='r:', continuous_update=False, disabled=False),\n",
    "         T= widgets.FloatLogSlider(value=0.0001,base=10,min=-6, max=2, step=0.1, continuous_update=False, readout=True,readout_format='.2e'),\n",
    "         S= widgets.FloatLogSlider(value=0.0001,base=10,min=-6, max=0, step=0.1, continuous_update=False, readout=True,readout_format='.2e'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Der kritische Zeitpunkt und die logarithmische Gerade\n",
    "Die rote, gestrichelte Linie repräsentiert den Zeitpunkt, ab dem die Abweichung der logarithmischen Geraden zu der Theis'schen Brunnenfunktion als vernachlässigbar angesehen werden kann. Dies ist dann der Fall, wenn *u* besonders klein ist ($\\leq 0.02$), was mithilfe folgender Formel überprüft werden kann:\n",
    "\n",
    "<br>\n",
    "<center>$\\large u \\leq 0.02 = \\frac{r^2 S}{4 t_{min} T} \\Rightarrow  t_{min} = \\frac{r^2 S}{0.08 T} $</center>\n",
    "<br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Für jeden Zeitpunkt, der größer ist als $t_{min}$ kann also bedenkenlos auf die logarithmische Gerade zurückgegriffen werden. Dies bedeutet einerseits, dass ein Pumpversuch für die Auswertung nach Cooper-Jacob möglichst lange laufen sollte und andererseits, dass alle Abweichungen vom \"Normalverhalten\" zu früheren Zeitpunkten ignoriert werden. Das sorgt dafür, dass mit Cooper-Jacob eine sehr gute Aussage über die Aquifereigenschaften auf großskaliger (regionaler) Ebeen getroffen werden kann. Auf der anderen Seite kann mit den hier vorgestellten Methoden keine Aussage über Aquifereigenschaften des direkten Brunnenumfeldes getroffen werden (z.B. die Anbindung Brunnen-Aquifer). \\\n",
    "Wenn die eben präsentierte vereinfachte Brunnenfunktion genutzt wird, vereinfacht sich auch die Berechnung der erwarteten Absenkung zu einer Funktion, die eine logarithmische Gerade darstellt:\n",
    "\n",
    "<br>\n",
    "<center>$\\large s(r,t)=\\frac{2.3 \\, Q}{4\\pi T}log\\frac{2.25Tt}{r^2S} $</center>\n",
    "<br>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Der Einfluss der einzelnen Parameter\n",
    "Im folgenden Graphen können Sie an vier Parametern drehen:\n",
    "\n",
    "__r__: Der Abstand der betrachteten Messstelle zum Zentrum des Pumpbrunnens\\\n",
    "__T__: Die Transmissivität\\\n",
    "__S__: Der Speicherkoeffizient\\\n",
    "__Q__: Die Pumprate\\\n",
    "\\\n",
    "<em> <font color=blue>Aufgabe 1: Spielen Sie ein bisschen mit den Parametern und versuchen Sie ein Gefühl dafür zu bekommen, welchen Einfluss diese auf den Verlauf der logarithmischen Geraden, den Verlauf von *u* und die Lage von $t_{min}$ haben. Notieren Sie, auf welche Art und Weise die einzelnen Parameter die Funktionen beeinflussen und überlegen Sie, worauf dieses Verhalten zurückzuführen sein könnte. Was könnte das für die Praxis bedeuten?</font> </em>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "096fa647a8cf4530ba3c56a76cc8c235",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=10.0, continuous_update=False, description='r [m]:', max=1000.0, min=1…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.splot(r, T, S, Q)>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#PLOT s and u as function of t\n",
    "\n",
    "def splot(r, T, S, Q):\n",
    "    tmax = 1000000\n",
    "    s_max = 1\n",
    "    r_max = 1000000\n",
    "    s = [s_max for x in range(r_max)]\n",
    "    u = [s_max for x in range(r_max)]\n",
    "    t = [s_max for x in range(r_max)]\n",
    "    \n",
    "\n",
    "    for x in range(1,tmax,1):\n",
    "        if x>0:\n",
    "            s[x] = (2.3*Q)/(4*3.1416*T)*math.log((2.25*T*x)/(r**2*S))\n",
    "            u[x] = r**2*S/4/T/x\n",
    "            t[x] = x\n",
    "            uu = r**2*S/4/T/0.02\n",
    "\n",
    "    #PLOT\n",
    "    fig = plt.figure(figsize=(9,6))\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(t, u, 'g-')\n",
    "    ax.plot(t, s, 'b-')\n",
    "    3\n",
    "    #plt.yscale(\"log\")\n",
    "    plt.xscale(\"log\")\n",
    "    plt.axis([1,1E5,1E-1,1E+1])\n",
    "    ax.set(xlabel='t [s]', ylabel='s [m] bzw. u [-]')\n",
    "    plt.title('Absenkung s und Parameter u',fontweight='bold', fontsize=16)\n",
    "    ax.grid(which=\"both\")\n",
    "    plt.legend(('u','Cooper-Jacob s'),loc=1)\n",
    "    plt.axvline(x=uu, ls='--', color='red')\n",
    "    plt.show()\n",
    "    \n",
    "# Computation\n",
    "\n",
    "interact(splot,\n",
    "         r = widgets.FloatSlider(value=10, min=1, max=1000, step=1, description='r [m]:', continuous_update=False, disabled=False, readout_format='.1f'),\n",
    "         Q = widgets.FloatLogSlider(value=0.0001, base=10, min=-4, max=0-2, step=0.001, description='Q [m³/s]:', continuous_update=False, readout_format='.2e'),\n",
    "         T = widgets.FloatLogSlider(value=0.0001,base=10,min=-6, max=2, step=0.1,readout=True,description='T [m²/s]:', continuous_update=False,readout_format='.1e'),\n",
    "         S = widgets.FloatLogSlider(value=0.0001,base=10,min=-6, max=0, step=0.1,readout=True,description='S [ - ]:', continuous_update=False,readout_format='.1e'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Die Bestimmung von Transmissivität und Speicherkoeffizient\n",
    "Im oben dargestellten Graphen sind die Parameter *T* und *S* bekannt (weil wir sie selber festgelegt haben). Üblicherweise führen wir einen Pumpversuch aber ja durch, um genau diese Parameter zu bestimmen. Bekannt ist nur der Abstand der Messstelle zum Pumpbrunnen und die Fließrate (so sie denn ordentlich gemessen wird).\n",
    "Nehmen wir z.B. die folgenden Daten, die dem Lehrbuch von Schwartz und Zhang entnommen wurden:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DATA VARNUM (2016) (Q = 1138 m3/day)\n",
    "\n",
    "#R4 => r = 90 m\n",
    "#m_time = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325] # time in minutes\n",
    "#m_ddown = [2E-05,0.0303,0.05141,0.0661,0.07509,0.08408,0.09307,0.10115,0.10556,0.10997,0.11162,0.11878,0.12319,0.12577,0.13002,0.13245,0.13396,0.13822,0.14064,0.14307,0.14457,0.147,0.14851,0.14818,0.15108,0.15307,0.15322,0.15612,0.15902,0.161,0.16115,0.16313,0.16328,0.16343,0.16618,0.16618,0.16618,0.16893,0.16893,0.17077,0.17352,0.17077,0.17352,0.17352,0.17614,0.17602,0.17589,0.1776,0.17564,0.17552,0.17722,0.17527,0.17514,0.17685,0.1768,0.1795,0.17945,0.1794,0.18118,0.18113,0.18108,0.18103,0.18098,0.1791,0.17897,0.18068,0.18056,0.18227,0.18214,0.18202,0.18189,0.18177,0.17981,0.18152,0.18179,0.18207,0.18509,0.18262,0.18564,0.18592,0.18344,0.18647,0.18399,0.18427,0.18773,0.18662,0.18733,0.18622,0.18602,0.18673,0.18837,0.18633,0.18613,0.18777,0.18779,0.18782,0.18784,0.18787,0.18789,0.18792,0.18794,0.18797,0.18799,0.18802,0.19092,0.19107,0.19122,0.19137,0.19152,0.19167,0.19182,0.19197,0.19395,0.19227,0.19244,0.19262,0.19279,0.19297,0.19314,0.19332,0.19349,0.19367,0.19384,0.19402,0.19552,0.1952,0.19304,0.19455,0.19422,0.19207,0.19357,0.19325,0.19109,0.1926,0.19257,0.19255,0.19252,0.1925,0.19247,0.19245,0.19242,0.1924,0.19237,0.19235,0.19252,0.1927,0.19562,0.19305,0.19322,0.1934,0.19632,0.1965,0.19392,0.19685,0.19405,0.194,0.19395,0.19665,0.1966,0.19655,0.1965,0.19828,0.1964,0.19635,0.19635,0.19635,0.19818,0.19818,0.19635,0.19818,0.19635,0.19818,0.19818,0.19635,0.19801,0.19783,0.19766,0.19748,0.19731,0.19713,0.19971,0.19953,0.19936,0.19827,0.19661,0.19862,0.19879,0.19897,0.19731,0.19932,0.20041,0.19783,0.19801,0.20002,0.20063,0.19942,0.20003,0.19973,0.19668,0.20188,0.19792,0.19853,0.19823,0.19702,0.19793,0.19793,0.19977,0.19702,0.19793,0.19793,0.19702,0.19518,0.19793,0.19977,0.19808,0.20007,0.19747,0.20037,0.20052,0.19792,0.19898,0.19913,0.19928,0.19943,0.19938,0.19842,0.19837,0.19923,0.19827,0.19913,0.20183,0.19812,0.19898,0.19893,0.19804,0.19807,0.19901,0.19628,0.20181,0.19817,0.20094,0.20097,0.20191,0.20193,0.19953,0.19988,0.19932,0.20242,0.20002,0.20312,0.20072,0.20198,0.20417,0.20543,0.20142,0.20198,0.20347,0.20312,0.20002,0.20242,0.19932,0.19897,0.20137,0.20193,0.20122,0.19867,0.19887,0.19907,0.20293,0.20222,0.19967,0.20262,0.20007,0.20027,0.20376,0.20358,0.20249,0.20232,0.19939,0.20197,0.20179,0.20162,0.20144,0.20402,0.20233,0.20248,0.20263,0.20187,0.20293,0.20033,0.20232,0.20247,0.20262,0.20277,0.20294,0.20312,0.20329,0.20347,0.20364,0.20473,0.20399,0.20417,0.20434,0.20452,0.20159,0.20417,0.20399,0.20382,0.20364,0.2053,0.20512,0.20312,0.20477,0.20277,0.20477,0.20312,0.20329,0.20347,0.20547,0.20382,0.20399,0.20417,0.20617,0.20452,0.20417] # drawdown in meters\n",
    "\n",
    "#R12 => r = 40 m\n",
    "#m_time =  [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325] # time in minutes\n",
    "#m_ddown = [2E-05,0.02022,0.04591,0.0716,0.09342,0.11433,0.12882,0.14332,0.15139,0.16313,0.17396,0.18203,0.18827,0.1936,0.19878,0.2012,0.20729,0.21247,0.21489,0.22007,0.22249,0.22583,0.22826,0.23068,0.23358,0.23648,0.23938,0.24228,0.24243,0.24533,0.24915,0.2493,0.2522,0.25235,0.2551,0.2551,0.25785,0.25785,0.2606,0.2606,0.26335,0.26335,0.2661,0.2661,0.26597,0.26585,0.26847,0.2656,0.26822,0.27177,0.26797,0.27152,0.27139,0.27402,0.27397,0.27392,0.27387,0.27382,0.27652,0.27647,0.27642,0.27637,0.27907,0.27627,0.27614,0.27877,0.27589,0.27577,0.27839,0.27735,0.27814,0.2771,0.28064,0.2796,0.27712,0.2774,0.28042,0.27795,0.27822,0.2785,0.28152,0.27905,0.28207,0.28235,0.28307,0.28012,0.28175,0.2843,0.2841,0.28482,0.28462,0.28442,0.28422,0.28402,0.28404,0.28407,0.28684,0.28412,0.28689,0.28692,0.28694,0.28422,0.28699,0.28702,0.28717,0.28732,0.28747,0.28762,0.28777,0.28792,0.28807,0.28822,0.28837,0.28852,0.29144,0.28887,0.29179,0.28922,0.28939,0.28957,0.28974,0.28992,0.29009,0.29027,0.28994,0.29237,0.28929,0.29172,0.29139,0.29107,0.29074,0.29042,0.29009,0.28977,0.28974,0.28972,0.28969,0.29333,0.28964,0.28687,0.28959,0.29323,0.28954,0.28952,0.29336,0.29353,0.29371,0.29022,0.29406,0.29423,0.29441,0.29458,0.29109,0.29493,0.29488,0.29483,0.29478,0.29473,0.29468,0.29463,0.29458,0.29453,0.29723,0.29718,0.29443,0.29443,0.29443,0.29443,0.29443,0.29718,0.29443,0.29443,0.29443,0.29718,0.29701,0.29683,0.29666,0.29648,0.29631,0.29613,0.29596,0.29578,0.29561,0.29543,0.29561,0.29853,0.29596,0.29613,0.29631,0.29648,0.29666,0.29683,0.29701,0.29993,0.29688,0.29658,0.29628,0.29598,0.29568,0.29538,0.29508,0.29478,0.29723,0.29693,0.29418,0.29418,0.29418,0.29418,0.29693,0.29693,0.29418,0.29418,0.29693,0.29418,0.29433,0.29723,0.29738,0.29478,0.29768,0.29783,0.29798,0.29813,0.29828,0.29843,0.29838,0.29833,0.29828,0.29823,0.29818,0.29813,0.29808,0.29803,0.29798,0.29793,0.29796,0.29798,0.29801,0.29803,0.30081,0.30083,0.29811,0.29813,0.29816,0.30093,0.29853,0.29888,0.29923,0.29958,0.29993,0.30303,0.30338,0.30098,0.30133,0.30168,0.30408,0.30098,0.30063,0.30303,0.29993,0.30233,0.29923,0.29888,0.29853,0.29818,0.30113,0.30133,0.29878,0.29898,0.30193,0.30213,0.30233,0.30253,0.30273,0.30018,0.30001,0.30258,0.30241,0.29948,0.29931,0.29913,0.30171,0.29878,0.29861,0.29843,0.30133,0.29873,0.29888,0.29903,0.29918,0.29933,0.29948,0.29963,0.30253,0.29993,0.30011,0.30028,0.30046,0.30063,0.30081,0.30373,0.30391,0.30133,0.30151,0.30443,0.30151,0.30133,0.30116,0.30098,0.30081,0.30338,0.30046,0.30028,0.30286,0.29993,0.30286,0.30303,0.30321,0.30338,0.30356,0.30373,0.30391,0.30408,0.30426,0.30443,0.30408] # drawdown in meters\n",
    "\n",
    "#R14 => r = 160 m\n",
    "#m_time = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325] # time in minutes\n",
    "#m_ddown = [0.000,0.00739,0.01474,0.02577,0.03568,0.04284,0.04999,0.05349,0.05973,0.06322,0.06671,0.07387,0.07736,0.08085,0.0842,0.08387,0.08629,0.08964,0.09298,0.09265,0.096,0.09934,0.09901,0.09869,0.1025,0.10632,0.10647,0.11029,0.11044,0.11059,0.1144,0.11455,0.1147,0.11485,0.1176,0.1176,0.1176,0.12127,0.12127,0.12127,0.12127,0.12127,0.12494,0.12494,0.12481,0.12469,0.12456,0.1281,0.12431,0.12419,0.12773,0.12394,0.12748,0.12369,0.1273,0.12725,0.1272,0.12715,0.1271,0.12705,0.13067,0.12695,0.13057,0.13052,0.1304,0.13027,0.13015,0.13002,0.1299,0.12977,0.13331,0.13319,0.1294,0.12927,0.12955,0.13349,0.13376,0.13404,0.13064,0.13092,0.13486,0.13514,0.13541,0.13569,0.13549,0.13529,0.13509,0.13489,0.13469,0.13449,0.13429,0.13409,0.13389,0.13369,0.13371,0.1374,0.13743,0.13379,0.13748,0.1375,0.13753,0.13755,0.13758,0.1376,0.13775,0.1379,0.13805,0.1382,0.13835,0.1385,0.14232,0.13514,0.13895,0.1391,0.13928,0.13945,0.13963,0.1398,0.13998,0.14015,0.14033,0.1405,0.14068,0.14085,0.14053,0.1402,0.14354,0.14322,0.14289,0.1389,0.13858,0.14192,0.14159,0.14127,0.14124,0.14122,0.14119,0.14117,0.14114,0.13745,0.13743,0.14107,0.14104,0.14102,0.1412,0.14137,0.14154,0.14172,0.13823,0.14207,0.14225,0.14242,0.13893,0.14277,0.14272,0.14267,0.14262,0.14257,0.13885,0.1388,0.14242,0.14237,0.13865,0.1386,0.1386,0.14227,0.14227,0.14227,0.14227,0.14227,0.14227,0.14227,0.14227,0.14502,0.14209,0.14467,0.14174,0.14157,0.1414,0.14122,0.14379,0.14087,0.14069,0.14052,0.14069,0.14087,0.14104,0.14122,0.14139,0.14157,0.14449,0.14192,0.14484,0.14227,0.14472,0.14167,0.14137,0.14382,0.14352,0.14322,0.14292,0.14262,0.14232,0.13927,0.14202,0.14202,0.14202,0.14202,0.14202,0.14202,0.14202,0.14202,0.14202,0.14202,0.14217,0.14232,0.14247,0.14262,0.14277,0.14292,0.14307,0.14322,0.14337,0.14352,0.14347,0.14342,0.14337,0.14332,0.14327,0.14322,0.14317,0.14679,0.14307,0.14302,0.14304,0.14307,0.14309,0.14312,0.14314,0.14317,0.14319,0.14322,0.14324,0.14327,0.14362,0.14397,0.14432,0.14467,0.14502,0.14537,0.14939,0.14607,0.14642,0.14677,0.14642,0.14607,0.14572,0.14537,0.14502,0.14467,0.14432,0.14397,0.14362,0.14694,0.14347,0.14367,0.14387,0.14407,0.14427,0.14447,0.14467,0.14487,0.14507,0.14894,0.1451,0.14492,0.14474,0.14457,0.14439,0.14789,0.14405,0.14387,0.14369,0.14352,0.14367,0.14382,0.14397,0.14779,0.14427,0.14442,0.14457,0.14472,0.14487,0.14502,0.14886,0.14537,0.14555,0.14572,0.14956,0.14974,0.14625,0.15009,0.15026,0.15044,0.15026,0.15009,0.14625,0.14607,0.14589,0.14572,0.14554,0.14904,0.14886,0.14502,0.1452,0.14537,0.14921,0.14572,0.1459,0.14607,0.14625,0.14642,0.14659,0.14677,0.14642]   # drawdown in meters\n",
    "\n",
    "#R15 => r = 2,5 m\n",
    "#m_time = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325] # time in minutes\n",
    "#m_ddown = [0.04219,0.20264,0.24483,0.26044,0.28501,0.29492,0.30849,0.31657,0.32556,0.33272,0.33988,0.34612,0.35236,0.35769,0.36103,0.36345,0.36588,0.3683,0.37439,0.37407,0.37558,0.38075,0.38409,0.3856,0.38667,0.38957,0.3943,0.39354,0.39644,0.403,0.40315,0.40147,0.40345,0.4036,0.40635,0.4091,0.40635,0.4091,0.4091,0.4091,0.41277,0.41277,0.41277,0.41552,0.41264,0.41527,0.41514,0.41502,0.41764,0.41752,0.41739,0.42094,0.41623,0.41702,0.42064,0.42059,0.42054,0.42324,0.42319,0.4213,0.42309,0.42487,0.42574,0.42294,0.42556,0.42544,0.42531,0.42519,0.42415,0.42677,0.42848,0.42652,0.42823,0.43085,0.43113,0.42865,0.42893,0.4292,0.43223,0.42975,0.43094,0.42847,0.43058,0.43085,0.43157,0.43137,0.433,0.4328,0.43077,0.4324,0.42945,0.432,0.43364,0.43344,0.43438,0.4344,0.43443,0.43445,0.43448,0.43359,0.43361,0.43455,0.43824,0.4346,0.43475,0.43215,0.4323,0.43245,0.43535,0.4355,0.43565,0.43947,0.43595,0.43519,0.43811,0.43829,0.43846,0.43864,0.43881,0.43899,0.43916,0.43934,0.43951,0.43785,0.44119,0.43904,0.43871,0.44114,0.43806,0.44049,0.43741,0.43984,0.43951,0.43919,0.43916,0.43914,0.43911,0.43909,0.43906,0.43904,0.43901,0.43624,0.44263,0.4426,0.44278,0.43929,0.44313,0.4433,0.44348,0.43999,0.44016,0.444,0.44418,0.44069,0.44064,0.44425,0.44237,0.44232,0.43952,0.44222,0.444,0.44395,0.44665,0.44202,0.44202,0.44202,0.44202,0.44477,0.44202,0.44202,0.44477,0.44385,0.44202,0.44385,0.44368,0.44625,0.44608,0.4459,0.44573,0.4428,0.44538,0.4452,0.44228,0.44485,0.44503,0.4452,0.44538,0.4428,0.44573,0.4459,0.44608,0.44625,0.44918,0.4466,0.4463,0.446,0.4457,0.4454,0.44785,0.4448,0.4445,0.4442,0.4439,0.44177,0.44177,0.44177,0.44177,0.44177,0.44177,0.44177,0.44177,0.44177,0.44177,0.44177,0.44375,0.4439,0.44405,0.44695,0.44435,0.44175,0.44465,0.4448,0.44495,0.4451,0.4478,0.44775,0.44495,0.4449,0.44485,0.4448,0.44659,0.4447,0.44465,0.4446,0.44463,0.44465,0.44468,0.44287,0.44473,0.44292,0.44294,0.44297,0.44299,0.44302,0.44337,0.44372,0.44407,0.44442,0.44477,0.44512,0.44547,0.44582,0.44617,0.45019,0.44984,0.44582,0.44547,0.44512,0.44477,0.44442,0.44407,0.44372,0.44337,0.44302,0.44322,0.44342,0.44362,0.44107,0.44402,0.44605,0.44625,0.44645,0.44665,0.44685,0.44668,0.4465,0.44908,0.4489,0.44873,0.44855,0.44563,0.44545,0.44528,0.44785,0.44525,0.4454,0.44555,0.4457,0.44585,0.446,0.44432,0.4463,0.44462,0.4466,0.44494,0.44512,0.4453,0.44547,0.44565,0.44582,0.44599,0.44617,0.44635,0.44377,0.44359,0.44617,0.44599,0.44582,0.44564,0.44547,0.44529,0.44237,0.44494,0.44844,0.44494,0.44512,0.44529,0.44547,0.44565,0.44582,0.44599,0.44342,0.44634,0.44377,0.44617] # drawdown in meters\n",
    "\n",
    "#B1 => Pumping well\n",
    "#m_time = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325] # time in minutes\n",
    "#m_ddown = [1.06975,0.9955,1.034,1.0505,1.067,1.08075,1.08716,1.0945,1.10916,1.10366,1.11741,1.12016,1.12291,1.13116,1.13391,1.13391,1.14491,1.14216,1.15316,1.13941,1.15041,1.15041,1.15041,1.16141,1.15316,1.16141,1.16141,1.15866,1.16416,1.16966,1.16966,1.16966,1.17791,1.17241,1.17516,1.16966,1.17791,1.18066,1.18341,1.17791,1.17791,1.17516,1.17791,1.18341,1.18066,1.18341,1.17791,1.17791,1.18341,1.18341,1.18891,1.18341,1.18616,1.18066,1.18066,1.17791,1.18066,1.18891,1.18616,1.18066,1.18616,1.19166,1.19716,1.19166,1.19166,1.18616,1.19166,1.19716,1.19716,1.19166,1.19991,1.19441,1.19166,1.18891,1.18341,1.19166,1.18891,1.18891,1.19441,1.19716,1.18616,1.19441,1.19166,1.18891,1.18616,1.19166,1.18891,1.19441,1.19166,1.19441,1.19166,1.19441,1.18891,1.19991,1.19991,1.20266,1.20266,1.20633,1.19716,1.19716,1.19991,1.19991,1.19716,1.20266,1.19716,1.19716,1.19991,1.19716,1.19991,1.20633,1.20633,1.20633,1.20266,1.20908,1.21183,1.20908,1.20633,1.20908,1.20908,1.20633,1.20266,1.20633,1.20266,1.20266,1.20633,1.20908,1.21183,1.20908,1.21183,1.21183,1.20908,1.21183,1.20908,1.21183,1.20908,1.20908,1.20908,1.21458,1.20633,1.20908,1.21183,1.21458,1.21183,1.21458,1.21183,1.21733,1.20908,1.20908,1.21458,1.21458,1.21183,1.20908,1.21183,1.21183,1.21183,1.21458,1.21458,1.21183,1.21733,1.21733,1.21458,1.20908,1.21458,1.21733,1.21733,1.21733,1.21733,1.21733,1.21733,1.21458,1.21183,1.21458,1.21458,1.21183,1.21183,1.21733,1.22008,1.20908,1.21458,1.22008,1.21733,1.22008,1.21733,1.22283,1.22283,1.22283,1.22833,1.22008,1.21183,1.21733,1.22008,1.20908,1.21183,1.22008,1.22008,1.22558,1.22558,1.21458,1.21458,1.21733,1.22008,1.21733,1.20908,1.22008,1.21733,1.21733,1.21183,1.21458,1.21733,1.20908,1.21458,1.21183,1.22008,1.20908,1.22008,1.21458,1.22283,1.21733,1.22283,1.22008,1.21733,1.21733,1.22283,1.21458,1.21733,1.22283,1.21183,1.22008,1.21183,1.21733,1.22283,1.21458,1.21183,1.22008,1.21733,1.22283,1.21733,1.21733,1.22283,1.21458,1.21458,1.22283,1.21733,1.22008,1.21733,1.21733,1.22283,1.21733,1.21458,1.22008,1.21183,1.22008,1.21183,1.22833,1.22558,1.22283,1.21733,1.21733,1.22283,1.22008,1.21458,1.22283,1.22008,1.21733,1.21458,1.21458,1.22283,1.21458,1.22008,1.21458,1.21733,1.21733,1.22558,1.22008,1.22008,1.22558,1.22008,1.22008,1.22008,1.21733,1.21458,1.21458,1.22833,1.21183,1.22558,1.21458,1.21733,1.21458,1.21458,1.21733,1.22008,1.21733,1.21733,1.22008,1.22008,1.21733,1.22283,1.22833,1.22008,1.22283,1.22008,1.22833,1.21733,1.22008,1.22283,1.21733,1.22558,1.22558,1.22283,1.21733,1.22833,1.22008,1.22008,1.22558,1.22008,1.22008,1.22283,1.22283,1.21183,1.22558,1.22833,1.21458,1.22558,1.22283,1.22558] # drawdown in meters\n",
    "\n",
    "\n",
    "\n",
    "# DATA FROM FETTER (2000) (TABLE 5.1), TRANSFERED FROM FEET TO METER\n",
    "# Q = 42400 ft3/d, r = 824 ft\n",
    "#m_time = [3,5,8,12,20,24,30,38,47,50,60,70,80,90,100,130,160,200,260,320,380,500] # time in minutes\n",
    "#m_ddown = [0.093,0.216,0.401,0.648,0.987,1.110,1.264,1.449,1.573,1.635,1.758,1.881,1.943,2.066,2.159,2.313,2.560,2.621,2.837,2.991,3.146,3.362]   # drawdown in meters\n",
    "\n",
    "\n",
    "\n",
    "#DATA FROM SCHWARTZ & ZHANG (2003)\n",
    "#Q = 500 m3/d, r = 300 m\n",
    "m_time = [1,1.27,1.61,2.04,2.59,3.29,4.18,5.3,6.72,8.53,10.83,13.74,17.43,22.12,28.07,35.62,45.2,57.36,72.79,92.37,117.21,148.74,188.74,239.5,303.92,385.66,489.39,621.02,788.05,1000] # time in minutes\n",
    "m_ddown = [0.03,0.05,0.09,0.15,0.22,0.31,0.41,0.53,0.66,0.8,0.95,1.11,1.27,1.44,1.61,1.79,1.97,2.15,2.33,2.52,2.7,2.89,3.07,3.26,3.45,3.64,3.83,4.02,4.21,4.39]   # drawdown in meters\n",
    "m_t = np.array(m_time)*60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGcCAYAAAA70rGtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7xcVX338e8vl8P9EG2qIpfngDUWgxIDjzrV6olBxUJtvQFW+zRiTa3F1hZRVKzx4anRItZSvGsaramtoVQpXpAn5hQp09aAQQhBvJ0iRUHEGMIlJ+T8+sfak7My7JnZc2b2XNb5vF+v/Tpz1l5777XX7Nnzm7X23svcXQAAAKmZ1+8CAAAAlIEgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgpwAzO87MPmBm3zaznWb2kJndbmbXmNkFZjbW7zIiMLMxM/No2mtmu8xs0sy+amavM7MDO9zGm8xsjZm9qVvlHlZmtj6q67EO1nOsmX3KzG7LPl87zOy7ZvZFM/v9WaxvTVSu8dmWq5eyY3dNNo2XtI3jzOxzZnZnVs8/MrMrzOyENtYR122zaTwnfyn71Wtm9tFsf/41Snuqmf21mV1vZndl9ft9M9tgZr/aYD2rzOw/zez+7LvlX83s9AZ5T8/m78zy/6eZrcrJ96dm9m9m9hMzmzKzB8zsFjO7yMweXZf3mmw/PtxhlQwud2dqMkn6P5IelORNpvP7XU6mfe/XWIv3yiXdJGmsg21MZuuZ7Pf+9nuStD6q11nVqaQnSvpFk/fr/89inWui5cf7XU8FyzwelXlNCesfjY7d+um3Z1m3zaZlw/petNj/4yU9nO3PC6L085vUxYOSnl63nvc0yb+6Lu8fNMn7nrq8E03yfkvSvCjvC7L0PZJ+td91W8ZES04TZrZC0jpJByocCB9QOCGPSFqscIB8XNID/SojmvovdzdJh0p6nqTrsvQTJH3ZzA7oW8kQ+xOFL2BJekP2elTS0yS9VSEoTVanLYtt+HVJ/yt7famkwyUdK+mPJP130ZW4+xp3t3hSOEfeEGX7kLtv7VK5B82fSZqvEDBeHaW7pC9LOlXhnHOspK9n8w6UdEEto5mdqBAUSdK2LO9TJf04S/srM3tslvexCt89yuY/VeHH3LYs7a1m9tSoHJ+XtFLSYyQdIukVkqayecuy5WuulvRfkhZISrNlut9R1iBPkv5dMxHwpW0sd7ykv1M4cUxJulvSZZKeWpdvfbT+l0j6rKT7FA7kt0sySa+R9N0sfULSk+vWMU/hJPVNSbsUfjHcJOnNkhbU5X2dpC2S7pW0Oyvf1ZJ+L8pzrKTPSLpd0kOSdki6OSvrY6J8tQ/tTQpB3v1ZGc6u2+Z4tI/vlnSupO9ly2yV9KKc+jtV0lVZOacUTiZ/I2lxgbofi7Y3WTfvEEl3RPNfH817t6SqpLuybd4v6dvZ+zCSsy/102S0rqMkfUTSD7N1/VzSVyQ9p648a6Llf0vSh7Lt/0LS1yQ9scD+rorWsVrSX0j6UXYsXKdH/npsuZ9R3pHs/dqicPw9qHAsXtzgGP7f2f/3SvqZwjH/uAL7cFW0jscXyG+SXivp2qyudmfHyIYGdXuapA9K+kle3Ur6gmZ+zR6cpa2Mln9xlrYwqy+X9LUs7RhJ/yDp1ux93iPpnmyfnt/k/X6JpE9leb2uHuunNWrdQrmmQL2tiPJ/tsvnykujdX9L0gGzeS+y/B9RCJh+mtXnfZL+U9IfSrIOzy3nKZzbHpT0bwrH7GRtPQX281CFz5ZLem/dvMNy8p8clfHWKP3iKP3VUfo7o/Q3ZmlvjNIuiPL+bpR+cYtyXx/lXVo3731Z+n3Kjv+Upr4XYFAnhSg4PokcWXC5Z2cfsrwT0YOSfj3Kuz6a99Oc/F/MSfu+suBFIcD5QoNtuaR/qZ0UFKL5Rvkui8q0rUm+E7I8B2v/ALB+ujRa33iU/vOcvFOSjovyn9tkvT9QFGg1qP+xKP9kzvy3RfO/FKXf2mS7f5uzL/XTZJbnSQ3eS5e0V9KZ0TbXtKibW1UXqObsz6oW6/iZpMPb2c8s34EKQUTDfc05hn+Sk7dlV5OkT0b5f5qtc7Wk43PymqSNjfahQd3mlWtf3Ur64yj9eTnLX5SlVaK0t2Vpz2xSn3slrWhQpnvq8q5vsp41ah3kvKtAPY8o/LCpLfOGLp0rXx6t8z5JS+rmF34vsvwPFdlPtX9ueWdOnp3ZtN/x02RfT42WPa3g90Et/9VR+jVR+tOi9JdG6Z/J0j4Tpb0kyrs8Sr+mwfZrLTm71eDzKOn0aD3Pb7VPwzb1vQCDOkl6evxBqJu3te6Dck80b3uWNpkdhCMKTYR3Z+nfjvKuj9Zxk6SjJf123br/n0LT/eVR2rOy5c+K0t4j6VGSDpP0V1H6S7O8f6PoJJSV65jsA/D7WZ5fipb7a4Vg5lEKv0YukHR0li/ue/6j7IO0WNI/RunLs7zjUdqUwi/YwxVareq/MI7O8rhCy8cxkg6QdGaU95IW79tYlHcyZ/5Lovnb6tKPz+p6oaQnKPwirX1ZPTrKO9lk/V/N5u3I9v0ASb8SHRc/1UzL0Jq4rJJOlHSEpFui9EqL/V0V5f25pOdKerSkzVH6K9vdT0lvietJ4cv8YElPlvSWBsfwDZKOU+jSvStKP6LFPiyL3vf66WZlgUeW94xo3h2Snq9w/D1BUWtGO3Wr0H1ZS3tXlrYp+39a0r9naW+N8j0jSztC0oslHZm91wdr/y+NLzQo088kvVDSQZKekvNZWdOiziqaaVW6Xdlns0n+wxRaLuK63SvprCjPj7L0D7ZxnjxO4VivrfNVOXkKvxdZ/t9R+MwcqnCMPiUq272a+eEW11erc8vhmvnxOSXpN7K0D8Z1UmB/40CpVZ0vUOiuquV/WTTvO1F6HIjFLYhfzdLils74s/CEKP3Wum3HdVOb/kU5LTUK591angta1cGwTX0vwKBO2j/I2VE3LzfIUTi5552o66fHZfnXR2m1QOOgug/uQVn66ij9lVnahgLb+miW9880c9L+jMJ1EC/Q/r/y52nmF9FtCgHWq5S14ET56k+WedP5Wd74wxa3GMVfBLUyvq7Aem9t8b6NRXknc+a/LJofBzmnKAQod2vmosJ4ekaUdzJv/dl7l7ds/fTMLP+aKO2caD3vj9LParG/q6K874/Sz6l/L9rZT+3fivPcJttfH+U7PUq/rH5/W+zHMklXKP9X/H2aCbDjL7Dfa7K+tupWM0HZ1ZrplvqOQlfeHoXg5UtZnl9Imp8tNyLpHQrnhF16ZNm3NyjT23PKPB7NX9Nk305Q+LJ3hRahR7R45SxTO1c8pPDZu00z55gXSnqswrnBJa0qeI4cUejKrJV5XZfeizMVuubvVQjE6uv0sbM4t8QtMHHgebCiALvAPn8oWs9BTfIdoP1/mF5SN79RkHNKlN4qyPmVKL1IkOMKn/15dXnj75zCl2UMy8SFx41NRq8PN7PH1f5x92UeLrar95iC6/6lRttz9wejtLuj/6ei9NoFs0W2V9vWhxWa+acV+nI/qPDhucvMzs+2XZt3h0LA9g6FL5WbzOwmMzt6FtuNfSd6fX/0unbh5WzX2474Vs4fSpKZPUuhLl4o6ZcVLiqsd1CBdT+6wbL1ZlM3RTRdR5v7+dgo7ZZubL8Zd9/q7i9WqJuVCl+Au7LZhyp8SZVZrs3Z34qkZyh8+V0r6RsKv8iflU1S6BrYm72+ROHHwIkKLUr1Gh033ypY9v1kt+lfpdDCer9Cl8n2FsscqBA4SNJmd79SoT7vUgjoLlf44WMKXepXFizOxZJOyl5vVwisW2l1jJ6lcI3TcxX2Me87Kq9OW73Hi6O022sv3P0BhVa1ovLO+/tnMDtM4QLkl2RJH1b4URm7K3q9KHo9Gr2+exZ5JUnuPpF9R40q/Ji9M5v1QoXr/+YMgpwG3P1uhYvdat5cYLH4QLva6+5AyA66ee6+LWfZhwumNdresxps74xsfx5y9zMUvoifLelsSf+hEDC9x8yOzPJdqdBN9CSFZvj/q/Br6gTN3B0Qb/eoBtt9S05590SvvcX+vKPBeosGko9gZocqdK/VXJH9fYVmPgvvU7iA0BRO/nnyyi7N/PKUpO82ef+/lLNsq7opotU62tnP+MR6fJe2n8vM9p2s3f1+d/+6u58naW2UrfZ8j7LKVQtyDlG4RkcKAc412etzFLo34rxS6DKWwjUPv6YQNMRfPo08mJPWtM6yu2yulvR4hX16ubv/R4FtHa6ZgPZwSXL3HyhcBLxLIaB7QTb/ve5+T6sVmtnLNBPUPCjpjCxgaKXVe3FW9PqNCq0lpv3v3JrNeuN9OrL2wswOVns/nH4SvV5cP9PMFit0UT0vS3q3u/+RZ00mkW9Gr58cvV6ak6edvPtx9/vc/WqFH7g1S+qyxftxlxJDkNPcOxRaPiTpT83sXWZ2lJktNLP6A0Xu/l2FZmBJen720LhF2XSymf25wq+Ubol/cV1iZiea2YHZA79eamZXSnqOFE5KZnaOwgf8RoWD/sZsWVO4I0hm9jcKv6R3KTRt/pPCCVwKwU/9dj9lZk/MtrvEzF5tZtdq5lbVdlylmZPVm83sVDM72MyOMLMVZvZRhesi2pKtY1zhC+KILPkWSZ/OXsfB5C5JD5vZaQr99nlqv/wW14JDaV8r3Kbs3yea2V+a2WPM7LDsvfmzaH4/tLOfV0SvP2RmTzezg7L3+Lwul+tSCw/9e4WZPT77fB0n6UVRnlprRVyuvzCzldn7O2Zm7+ygDF+PXr8s+/uNbJKk32yQt1an0wrXphwi6aJZliFuUTjezEZq/5jZ4Qqfx1oXxSp3/2rB9d6tcJuwJFWy89gihbs4J6J8eyV9rtXKsvfmU1HSm9z95oJlaSU+RneGzdlrFB4n0ImqZgLL08zslKxO36MQmBYVBxPxrdjKWrq/oXAN47TC3ZtrGqyndjGxJJ2fHb9PUbiDTArXD30+e/15zTym5A1m9pSsRa92LvRsfTKzZWa2Nvu8Pio699WOaSncvBI7scH+paHf/WWDPkn6fc1cmd5oii88fq6aPzxwIsq7Pkofj9JraZNR2qoofVWWNk8h4GhWtvEs7wVN8typmWt/ml1T8qYsz8EKH4Zm2x3L8o5HaWui/YnT10fpb2mx3jUt3q+xFsu7wnUWY9Eyz9bM9Qi1aa/C7ah578+lOetcn807XuHLqtG24/d0TYP1x+mrWuzvI46LnPQ17e6nZnd31ViD9PEW+/DZBtupTVs0cyfUbO6uKlS3Ct0Y+z4TUXpcP/do/9uYP5FTjtuUX1e5ZYrmH6j8O/PGJf1eizpq9bk4TeEHRKvPxnZJi1qsq9kdnfG0qt33QuEawPr1PKCZC4/3HWdq/9ySd3fVfdm03/HTZN8P08wF3/W3kK/JWX/u8Znlb/YwwD+oy1voYYBqfC1ObbpW0sK6ddduIb9f0qGt6mDYJlpyWnD3TypEuh9WOHk9qPCh+77CL7q3KFykXMv/rwr91J9RuLZlj0I3xrcV7nB6exfLNq3Qv3qOQtfTLoWA7L8UfvWdo5lm3k2S/l7hhL1L4cvtx8r6v33m2p/3KnwQahemPpCt448V7riSh2bp5ygETjdmeR5UuMX7nxW6wmp9wO3u018qtCx8RSFYeFihibiqcBL5TLurzMp3u8IzOVYrPDtmMtrmtQon11sV6u8WhW6+axusc41Cvf00p/zbFS6i/YhCfUwpXKh6i8Kv39e3Wf6uaWc/3f0hhSb38xSesXF/tsz3FFr3uumDCifaqsItzrsVLpDdnqWvdPeHs3J5VubXKVwAv1Ohjm9X562km6PXcZ18I3o9kZWh5k8lfVTh83K/wo+OU2az8azOz1So7/qun5bXgrRY95cUutP+WSFQqz1/6gqFVqqLs6y/KmmjmS1osrpFTeZ1xN03KNTpDxWOgS0KLXr1rQ+zWfeFCq0fdygcY9dp//eq5bU57n6fwl2kknSmmc36fXH3tys8B22Lwvlzl0L36G+6+8fq8n5M4X26RjPPQ9si6TXZemp+IOlvFT47v1A4z/9cYV/PlXSKu+/r2svKf0b27z+4+y4lxvb/vAIAkB4zO0bhgaI3ZP8vUHgcxoVZln9097MaLR+t58kKP1rnS3qhu3+tpCKXzsxeoHCZwF6Fh9UWvZh/aNCSAwCYC54q6fpscMvbFVoBawHOXSrYyp4FAp/M/n1H10vZW7WbST6RYoAjldySY2aTCv2deyU97O4nl7YxAAAaMLMnKXR/nqzwCIVphW6xLys81Tq5O4vQmyDnZC9wSyIAAEA30V0FAACSVHZLzg81M0zAx9z94zl5Vivc8aIDDzzwpGOOOaY+S9ft2bNHCxe282iE2S3bKm+z+Y3m5aXXp7X6vyzUazk63U7R5Yvka7duZ5s2DHXLMdvYoNRrszyd1Gte2lyr214ds7fddts97v7LhQqdp8z70yU9Pvv7GIVbjZ/TLP+SJUu8FzZv3tyTZVvlbTa/0by89Pq0Vv+XhXotR6fbKbp8kXzt1u1s04ahbjlmGxuUem2Wp5N6zUuba3Xbq2NW0hYf1OfkuPud2d+7FZ7P8PTmSwAAAHRHaUGOmR2SDVQmMztEYWyUbj36GwAAoKlmT7Xs1GMl/XP2QMgFkv7ei4+1AgAA0JHSghwPo9ye2DIjAABACbiFHAAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAJIkgBwAAdEW1Km3YcIyq1X6XJCDIAQAAHatWpZUrpXXrjtXKlRqIQIcgBwAAdGxiQpqakqanTVNT4f9+I8gBAAAdGx+XRkakefOmNTIS/u83ghwAANCxSkXatEk6++xJbdoU/u+3Bf0uAAAASEOlIu3efbsqleP6XRRJtOQAAIBEEeQAAIAkEeQAAIAkEeQAAIAkEeQAAIAkEeQAAIAkcQs5AABoqVoNTzEehIf8FUWQAwAAmqqNSzU1FZ5qfNFFo0MR7NBdBQAAmqqNS7V3b/i7deuifhepEIIcAADQVG1cqvnzw99ly3b0u0iF0F0FAACaqo1LVbsmZ/funf0uUiEEOQAAoKVKZWbQzYmJvhalMLqrAABAkghyAABAkghyAABAkghyAABAkghyAABAkghyAABAkghyAABAkkoPcsxsvpl9y8yuLHtbAACgfdWqtHZt+JuSXjwM8E8kbZc02oNtAQCANtQPvrlpU79L1D2ltuSY2VGSTpP0yTK3AwAAZqd+8M1heZpxEWV3V31Q0lskTZe8HQAAMAv1g2+Oj/e7RN1j7l7Ois1Ol/Qb7v4GMxuX9GZ3Pz0n32pJqyVp8eLFJ23cuLGU8sQmJyc1NjZW+rKt8jab32heXnp9Wqv/y0K9lqPT7RRdvki+dut2tmnDULccs40NSr02y9NJvealDXvdbts2qq1bF2nZsh1aunTnwByzK1asuN7dT25YkFbcvZRJ0lpJd0ialPQTSQ9I+myzZZYsWeK9sHnz5p4s2ypvs/mN5uWl16e1+r8s1Gs5Ot1O0eWL5Gu3bmebNgx1yzHb2KDUa7M8ndRrXtpcq9teHbOStngHsUhp3VXu/jZ3P8rdxySdJenr7v7qsrYHAAAQ4zk5AAAgSb24hVzuPiFpohfbAgAAkGjJAQAAiSLIAQAASSLIAQAASSLIAQAASerJhccAAKA/qtUwVMPo6GhSTzMugiAHAIBExYNvLlhwopYvlyqVfpeqd+iuAgAgUfHgm3v2WFKDbxZBkAMAQKLiwTcXLnS6qwAAQBoqFWnTpto1OTeqUlne7yL1FEEOAAAJq1TCNDGxs99F6Tm6qwAAQJIIcgAAQJIIcgAAQJIIcgAAQJIIcgAAQJIIcgAAQJIIcgAAQJIIcgAAQJIIcgAAGFLVqrR2bfiLR+KJxwAADKF4hPGRkTB8w1waYbwIWnIAABhC8QjjU1OacyOMF0GQAwDAEIpHGB8Z0ZwbYbwIuqsAABhC8Qjj4+N0VeUhyAEAYEjVRhhHPrqrAABAkghyAABAkghyAABAkghyAABAkghyAABAkghyAABAkghyAABAkghyAAAYQAy+2TkeBggAwIBh8M3uoCUHAIABw+Cb3UGQAwDAgGHwze6guwoAgAHD4JvdQZADAMAAYvDNztFdBQAAkkSQAwAAkkSQAwAAkkSQAwAAkkSQAwAAkkSQAwAAkkSQAwAAkkSQAwBAjzH4Zm/wMEAAAHqIwTd7h5YcAAB6iME3e4cgBwCAHmLwzd6huwoAgB5i8M3eIcgBAKDHGHyzN+iuAgAASSLIAQAASSLIAQAASSLIAQAASSLIAQAASSLIAQAASSLIAQAASSLIAQCgixh8c3DwMEAAALqEwTcHCy05AAB0CYNvDhaCHAAAuoTBNwdLad1VZnagpGskHZBt5zJ3f1dZ2wMAoN8YfHOwlHlNzm5Jz3P3XWa2UNK1ZvYVd//3ErcJAEBfMfjm4CgtyHF3l7Qr+3dhNnlZ2wMAAIhZiEVKWrnZfEnXS/oVSR9y97fm5FktabUkLV68+KSNGzeWVp6ayclJjY2Nlb5sq7zN5jeal5den9bq/7JQr+XodDtFly+Sr926nW3aMNQtx2xjg1KvzfJ0Uq95aXOtbnt1zK5YseJ6dz+5UKHzuHvpk6RFkjZLOqFZviVLlngvbN68uSfLtsrbbH6jeXnp9Wmt/i8L9VqOTrdTdPki+dqt29mmDUPdcsw2Nij12ixPJ/WalzbX6rZXx6ykLd5B/NGTu6vcfYekCUmn9mJ7AAAApQU5ZvbLZrYoe32QpFMk3VrW9gAAAGJl3l11hKRPZ9flzJP0eXe/ssTtAQAA7FPm3VXflvS0stYPAADQDGNXAQBQULUaHvQ3OjrK04yHAEEOAAAFxINvLlhwopYv56F/g46xqwAAKCAefHPPHmPwzSFAkAMAQAHx4JsLFzrdVUOA7ioAAAqIB98cHb1RlcryfhcJLRDkAABQUG3wzYmJnf0uCgqguwoAACSJIAcAACSJIAcAACSJIAcAACSJIAcAACSJIAcAACSJIAcAACSJIAcAAIWxqTZsOEbVar9Lgm4hyAEAzHm1wTfXrTtWK1eKQCcRBDkAgDmvNvjm9LRpakoMvpkIghwAwJxXG3xz3rxpjYyIwTcT0XTsKjN7dIF1TLv7ji6VBwCAnqsNvrlu3aTOPvs4VSr9LhG6odUAnXdmkzXJM1/SMV0rEQAAfVCpSLt3365K5bh+FwVd0irI2e7uT2uWwcy+1cXyAAAAdEWra3KKNNjRqAcAAAZO05Ycd3+o9trMHiXp6HgZd78hzgMAADAoWnVXSZLM7EJJqyR9X5JnyS7peeUUCwAAoDOFghxJZ0h6grtPlVkYAACAbin6nJybJS0qsyAAAADdVLQlZ62kb5nZzZJ21xLd/cWllAoAgC6qVsNTjMfHxTNw5pCiQc6nJb1P0k2SpssrDgAA3VUbl2pqKjzVeNOmfpcIvVI0yLnH3S8ptSQAAJSgNi7V3r3aNy4VrTlzQ9Eg53ozWyvpCu3fXXVDKaUCAKBLauNS1Vpyxsel3btbLYUUFA1yak89fmaUxi3kAICBVxuXKr4mh1HG54ZCQY67ryi7IAAAlKVSoYtqLmp6C7mZnd5qBUXyAAAA9FqrlpyLzOy/1XwU8vdIurJ7RQIAAOhcqyDnLkkfaJHnu10qCwAAQNe0GqBzvEflAAAA6KqiwzoAAAAMFYIcAACQJIIcAACQpEJBjpm9wswOy15fYGaXm9nycosGAEBr1aq0dm34C8SKPvH4ne6+0cyeLemFkt4v6SOSnlFayQAAaCFv8E0e+oeaot1Ve7O/p0n6iLt/UdJIOUUCAKCYvME3gZqiQc5/m9nHJJ0h6ctmdkAbywIAUIra4Jvz588MvgnUFO2uOkPSqZLe7+47zOwISeeVVywAAFrLG3wTqCk6QOcDki6P/v+xpB+XVSgAAIpi8E00QpcTAABIEkEOAABIEkEOAABIEkEOAABIEkEOAABIEkEOAABIEkEOAABIEkEOAGCgMQAnZqvoE48BAOg5BuBEJ2jJAQAMLAbgRCcIcgAAA4sBONEJuqsAAAOLATjRCYIcAMBAYwBOzBbdVQAAIEkEOQAAIEkEOQAAIEmlBTlmdrSZbTaz7Wa2zcz+pKxtAQAA1CvzwuOHJZ3r7jeY2WGSrjezq939lhK3CQAAIKnElhx3/7G735C9vk/SdklHlrU9AACAmLl7+RsxG5N0jaQT3H1n3bzVklZL0uLFi0/auHFj6eWZnJzU2NhY6cu2yttsfqN5een1aa3+Lwv1Wo5Ot1N0+SL52q3b2aYNQ91yzDbWzna2bRvV1q2LtGzZDi1durOr9dosTyf1mpc2iHXbybKDcsyuWLHienc/uVCh87h7qZOkQyVdL+mlrfIuWbLEe2Hz5s09WbZV3mbzG83LS69Pa/V/WajXcnS6naLLF8nXbt3ONm0Y6pZjtrGi27nuOveDDnKfPz/8ve667tZrszyd1Gte2qDVbafLDsoxK2mLdxCDlHp3lZktlPRPkja4++VlbgsAMFwYlwplK/PuKpP0KUnb3f0DZW0HADCcGJcKZSuzJedZkn5X0vPMbGs2/UaJ2wMADJHauFQXXhj+MnQDuq20W8jd/VpJVtb6AQDDj3GpUCaeeAwAAJJEkAMAAJJEkAMAAJJEkAMAAJJEkAMAAJJEkAMAAJJEkAMAAJJEkAMAKEW1Kq1dG/4C/VDawwABAHNXtSqtXBnGpBoZ4YnG6A9acgAAXcfgmxgEBDkAgK5j8E0MArqrAABdVxt8c2IiBDh0VaEfCHIAAKVg8E30G91VAAAgSQQ5AAAgSQQ5AAAgSQQ5AKo8xD0AAAv6SURBVAAgSQQ5AAAgSQQ5AAAgSQQ5AIC2MS4VhgHPyQEAtIVxqTAsaMkBALSFcakwLAhyAABtYVwqDAu6qwAAbWFcKgwLghwAQNsYlwrDgO4qAACQJIIcAACQJIIcAACQJIIcAACQJIIcAACQJIIcAACQJIIcAMB+GJcKqeA5OQCAfRiXCimhJQcAsA/jUiElBDkAgH0YlwopobsKALAP41IhJQQ5AID9MC4VUkF3FQAASBJBDgAASBJBDgAASBJBDgAASBJBDgAASBJBDgAASBJBDgDMIYxLhbmE5+QAwByxbduozjuPcakwd9CSAwBzxNatixiXCnMKQQ4AzBHLlu1gXCrMKXRXAcAcsXTpTsalwpxCkAMAcwjjUmEuobsKAAAkiSAHAAAkiSAHAAAkiSAHAAAkiSAHAAAkiSAHABKxbdsoQzYAEW4hB4AEVKvSueeeqIcfnhmyAZjraMkBgARMTEh79sxjyAYgQpADAAkYH5cWLpxmyAYgQncVACSgUpEuvvhG7dy5fN+QDbTmYK4jyAGARCxdupMWHCBSWneVma0zs7vN7OaytgEAANBImdfkrJd0aonrBwAAaKi0IMfdr5F0b1nrBwAAaMbcvbyVm41JutLdT2iSZ7Wk1ZK0ePHikzZu3FhaeWomJyc1NjZW+rKt8jab32heXnp9Wqv/y0K9lqPT7RRdvki+dut2tmnDULccs40NSr02y9NJvealzbW67dUxu2LFiuvd/eRChc7j7qVNksYk3Vw0/5IlS7wXNm/e3JNlW+VtNr/RvLz0+rRW/5eFei1Hp9spunyRfO3W7WzThqFuOWYbG5R6bZank3rNS5trddurY1bSFu8gDuE5OQAwBKpVMWQD0CZuIQeAAVetSitXhicZM2QDUFyZt5B/TlJV0pPM7A4ze21Z2wKAlE1MhACHIRuA9pTWkuPuryxr3QAwl4yPhxacWkvO+Li0e3e/SwUMPq7JAYABV6mELqoLLwx/K5V+lwgYDlyTAwBDoFIhuAHaRUsOAABIEkEOAABIEkEOAABIEkEOAABIEkEOAPTZtm2jPM0YKAF3VwFAH1Wr0rnnnqiHH555mjF3UQHdQUsOAPTRxIS0Z888nmYMlIAgBwD6aHxcWrhwWvPnzzzNGEB30F0FAH1UqUgXX3yjdu5crvFxuqqAbiLIAYA+W7p0Jy04QAnorgIAAEkiyAEAAEkiyAEAAEkiyAEAAEkiyAGAElWr0oYNx/A0Y6APCHIAoCTVqrRypbRu3bFauZJhG4BeI8gBgJJMTISnGE9PG08zBvqAIAcASjI+Hp5iPG/eNE8zBvqAIAcASlKphAE3zz57koE3gT7giccAUKJKRdq9+3ZVKsf1uyjAnENLDgAASBJBDgAASBJBDgAASBJBDgDMUrUqrV3L82+AQcWFxwAwC9u2jeq888JzcEZGxN1TwACiJQcAZmHr1kWampL27hUP+gMGFEEOAMzCsmU7NDIizZ8vHvQHDCi6qwBgFpYu3alNm0ILzvg4XVXAICLIAYBZqlQIboBBRncVAABIEkEOAABIEkEOAOTgGTjA8OOaHACoU61KK1fyDBxg2NGSAwB1JibEM3CABBDkAECd8XHxDBwgAXRXAUCdSkU8AwdIAEEOAOTgGTjA8KO7CgAAJImWHABzTrUauqJGR0e53gZIGEEOgDklvj18wYITtXw53VJAquiuAjCnxLeH79lj3B4OJIwgB8CcEt8evnCh010FJIzuKgBzSnx7+OjojapUlve7SABKQpADYM6p3R4+MbGz30UBUCK6qwAkhYE1AdTQkgMgGQysCSBGSw6AZDCwJoAYQQ6AZDCwJoAY3VUAksHAmgBiBDkAhkZtOIZmAQwDawKoIcgBMBS4qBhAu7gmB8BQ4KJiAO0iyAEwFLioGEC76K4CMBDi623ycFExgHYR5ADou23bRnXeeTPX21x00WhusMNFxQDaQXcVgL7bunXRftfbbN26qN9FApAAWnIA9EStO2p09JGtNMuW7dDIyExLzrJlO/pRRACJIcgBULr49u8FC07U8uX7dzstXbpzv+ttdu9mdHAAnaO7CkDHqlVpw4ZjGo78Hd/+vWeP5d7+XalIb3sb19wA6B6CHAAdqbXSrFt3rFauVG6gE9/+vXChc/s3gJ6guwpAR2qtNNPTtu8hffWtMfHt36OjN6pSWd6HkgKYa8zd+12GfczsPknf6cGmDpf0ix4s2ypvs/mN5uWl16fV/79Y0j1NS9od1Gs5OqnXdpYvki8nz2GHSE9cotAyPC199zbpvvubrLNI2jDULcdsY4NSr83ydFKveWlzrW57dcw+yd0Pa17UJtx9YCZJW3q0nY/3YtlWeZvNbzQvL70+Led/6nWO1ms7yxfJ127dzjZtGOqWY3bw67Xd+itar9Tt8Byzc/WanH/p0bKt8jab32heXnp9Wif71wnqtRydbrfo8kXytVu3naT1AsdsOQalXpvl6aRei267DINSt0NxzA5ad9UWdz+53+VIDfVaDuq1PNRtOajX8lC35ei0XgetJefj/S5AoqjXclCv5aFuy0G9loe6LUdH9TpQLTkAAADdMmgtOQAAAF1BkAMAAJJEkAMAAJJEkAMAAJI0FEGOmR1nZp8ys8v6XZYUmNkhZvZpM/uEmb2q3+VJBcdpeczst7Pj9Ytm9oJ+lycVZna8mX3UzC4zsz/sd3lSkp1nrzez0/tdlpSY2biZfSM7bsdb5S89yDGzdWZ2t5ndXJd+qpl9x8y+Z2bnN1uHu//A3V9bbkmHW5v1/FJJl7n76yS9uOeFHSLt1CvHaXvarNsvZMfrKkln9qG4Q6PNet3u7q+XdIYknvHSxCy+y94q6fO9LeVwarNuXdIuSQdKuqPVunvRkrNe0qlxgpnNl/QhSS+S9GRJrzSzJ5vZU8zsyrrpMT0oYwrWq2A9SzpK0o+ybHt7WMZhtF7F6xXtWa/26/aCbD4aW6826tXMXizpWkmbelvMobNexb/LTpF0i6S7el3IIbVexY/Zb7j7ixSCyHe3WnHpQY67XyPp3rrkp0v6XvbLd0rSP0j6LXe/yd1Pr5vuLruMKWinnhWi36OyPEPRZdkvbdYr2tBO3VrwPklfcfcbel3WYdLuMevuV7j7r0mi67qJNut1haRnSvodSa8zM86zTbQZJ0xn838u6YBW6+5XxR+pmZYEKXzpHtkos5n9kpl9VNLTzOxtZRcuIY3q+XJJLzOzj6h/468Ms9x65TjtikbH7BslnSLp5Wb2+n4UbMg1OmbHzewSM/uYpC/3p2hDLbde3f0d7v4mSX8v6RPRFzOKa3TMvjQ7Xv9O0qWtVrKgpMK1YjlpDR+97O4/k8SJrX259ezu90t6Ta8Lk5BG9cpx2rlGdXuJpEt6XZiENKrXCUkTvS1KUpp+l7n7+t4VJTmNjtnLFX6oF9Kvlpw7JB0d/X+UpDv7VJaUUc/loF7LQ92Wg3otB/Vanq7Ubb+CnG9KeqKZHWtmI5LOknRFn8qSMuq5HNRreajbclCv5aBey9OVuu3FLeSfk1SV9CQzu8PMXuvuD0s6R9JVkrZL+ry7byu7LCmjnstBvZaHui0H9VoO6rU8ZdYto5ADAIAkcVsbAABIEkEOAABIEkEOAABIEkEOAABIEkEOAABIEkEOAABIEkEOAABIEkEOgJ4xs0Vm9oYG88bM7EEz29piHRvM7F4ze3k5pQSQCoIcAL20SFJukJP5vrsva7YCd3+VeHQ+gAIIcgD00nslPcHMtprZRc0ymtkhZvYlM7vRzG42szN7VEYAiVjQ7wIAmFPOl3RCq9aazKmS7nT30yTJzA4vtWQAkkNLDoBBdZOkU8zsfWb26+7+i34XCMBwIcgBMJDc/TZJJykEO2vN7M/7XCQAQ4buKgC9dJ+kw4pkNLPHS7rX3T9rZrskrSqzYADSQ5ADoGfc/Wdm9m9mdrOkr7j7eU2yP0XSRWY2LWmPpD/sSSEBJIMgB0BPufvvFMx3laSrSi4OgIRxTQ6AQbFX0uFFHgYo6bmSHupJqQAMLXP3fpcBAACg62jJAQAASSLIAQAASSLIAQAASSLIAQAASfofiSsjYgdAQ7EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot Data\n",
    "fig = plt.figure(figsize=(9,6))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9)  # adjust plot area\n",
    "ax.plot(m_t, m_ddown,'bo', markersize=3)\n",
    "plt.yscale(\"linear\")\n",
    "plt.xscale(\"log\")\n",
    "plt.axis([0.1,1E5,1E-2,5E+0])\n",
    "    \n",
    "ax.set(xlabel='t [s]', ylabel='s [m]')\n",
    "plt.title('Gemessene Daten',fontweight='bold', fontsize=16)\n",
    "ax.grid(which=\"both\", color='grey',linewidth=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Fließrate beträgt *500 m³/Tag* und die Messstelle befindet sich in *300 m* Entfernung zum Pumpbrunnen. Um die Vereinfachung der Brunnenfunktion nach Cooper-Jacob auf diese Messwerte anzuwenden, plotten wir beides zusammen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "172ea01d46254d2da881c51a41caf0c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(BoundedIntText(value=10, description='r [m]:', max=1000, min=1), BoundedFloatText(value=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.splot(r, T, S, Q)>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Plot Function\n",
    "def splot(r, T, S, Q):\n",
    "    tmax = 1000000\n",
    "    s_max = 1\n",
    "    r_max = 1000000\n",
    "    s = [s_max for x in range(r_max)]\n",
    "    t = [s_max for x in range(r_max)]\n",
    "    \n",
    "    for x in range(1,tmax,1):\n",
    "        if x>0:\n",
    "            s[x] = (2.3*Q)/(4*math.pi*T)*math.log((2.25*T*x)/(r**2*S))\n",
    "            t[x] = x\n",
    "            uu = r**2*S/4/T/0.02\n",
    "\n",
    "    #PLOT\n",
    "    fig = plt.figure(figsize=(9,6))\n",
    "    ax = fig.add_subplot(1, 1, 1)\n",
    "    ax.plot(t, s, 'green')\n",
    "    ax.plot(m_t, m_ddown,'bo', markersize=3)\n",
    "    3\n",
    "    #plt.yscale(\"log\")\n",
    "    plt.xscale(\"log\")\n",
    "    plt.axis([10,1E5,1E-1,5E+0])\n",
    "    ax.set(xlabel='t [s]', ylabel='s [m]')\n",
    "    plt.title('Auswertung nach Cooper-Jacob',fontweight='bold', fontsize=16)\n",
    "    ax.grid(which=\"both\", color='grey',linewidth=0.5)\n",
    "    plt.legend(('Cooper-Jacob','Messwerte'),loc=0)\n",
    "    plt.axvline(x=uu, ls='--', color='red')\n",
    "    plt.show()\n",
    "# Computation\n",
    "\n",
    "interact(splot,\n",
    "         r = widgets.BoundedIntText(value=10, min=1, max=1000, description='r [m]:', disabled=False),\n",
    "         Q = widgets.BoundedFloatText(value=0.0001, min=1E-4, max=1E-1, description='Q [m³/s]:', disabled=False),\n",
    "         T = widgets.BoundedFloatText(value=0.0001, min=1E-5, max=1E-1, description='T [m²/s]:', disabled=False),\n",
    "         S = widgets.BoundedFloatText(value=0.0001, min=1E-7, max=5E-1, description='S [ - ]:', disabled=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Iteratives Vorgehen__ \\\n",
    "<font color=blue> <em>Aufgabe 2: Versuchen Sie, durch richtige Abschätzung der Parameter die logarithmische Gerade möglichst gut an die Messwerte anzupassen.</font> </em> \\\n",
    "Nach Einstellung der richtigen Pumprate und des Abstandes zur Messstelle können wir durch fröhliches Rumprobieren (wissenschaftlich gesprochen: iterativ) versuchen, die Parameter *T* und *S* zu bestimmen.\n",
    "Denken Sie daran: Für die Anpassung der Geraden an die Messwerte dürfen Sie nur die Werte verwenden, die sich __rechts__ der roten Linie befindet, die den kritischen Zeitpunkt darstellt, ab dem die vereinfachte Brunnenfunktion gilt. \\\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__\"Klassische\" Berechnung nach Cooper-Jacob__ \\\n",
    "Ein solches iteratives Verfahren ist grundsätzlich zulässig. Es gibt aber auch einen anderen Weg, der oft schneller und auch anwendbar ist, wenn man gerade kein eigens dafür programmiertes Worksheet zur Hand hat und/oder die Werte nur auf Papier vorliegen (beispielsweise in einer Klausur).\n",
    "So lässt sich die logarithmische Gerade der angepassten Brunnenfunktion auch von Hand durch die vorab ausgewählten Messwerte zeichnen.  \\\n",
    "Der Grund, warum wir nun aus unserer vereinfachten Brunnenfunktion nicht sofort Transmissivität und Speicherkoeffizient ablesen können, ist, dass wir __zwei__ unbekannte Parameter (*T* und *S*) in __einer__ Gleichung haben:\n",
    "\n",
    "<br>\n",
    "<center>$\\large s(r,t)=\\frac{Q}{4\\pi T}log\\frac{2.25Tt}{r^2S} $</center>\n",
    "<br>\n",
    "\n",
    "__Hintergrund: Subtraktion zweier Absenkungen__ \\\n",
    "Das Problem lässt sich durch einen kleinen Trick umgehen, denn wenn die Absenkung an zwei unterschiedlichen Zeitpunkten bekannt ist, kann die Gleichung durch Subtraktion auf eine Unbekannte (*T*) reduziert werden:\n",
    "\n",
    "<br>\n",
    "<center>$\\large s_1=\\frac{Q}{4\\pi T}log\\frac{2.25Tt_1}{r^2S}     \\hspace{35pt}         s_2=\\frac{Q}{4\\pi T}log\\frac{2.25Tt_2}{r^2S} $</center>\n",
    "<br>\n",
    "\n",
    "<br>\n",
    "<center>$\\large s_2-s_1= \\Delta s=\\frac{2.3 \\, Q}{4\\pi T} log\\frac{t_2}{t_1} $</center>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Bestimmung der Transmissivität__\\\n",
    "Jetzt suchen wir uns zwei beliebige Zeitpunkte auf der logarithmischen Geraden und bestimmen die Differenz in der Abweichung $\\Delta s $. Im Idealfall liegen diese genau eine Größenordnung auseinander (z.B. die Zeitpunkte 10² und 10³ Sekunden), denn dadurch wird die Gleichung weiter vereinfacht, da: $ log\\frac{1000}{100}=1 $ . Die Transmissivität lässt sich nun ganz einfach mit folgender Formel bestimmen:\n",
    "\n",
    "<br>\n",
    "<center>$\\large T = \\frac{2.3 Q}{4 \\pi \\Delta s} $</center>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Bestimmung des Speicherkoeffizienten__\\\n",
    "Um den Speicherkoeffizienten zu bestimmen, muss die vereinfachte Brunnengleichung erneut angepasst werden. Da *T* inzwischen bekannt ist, lässt sich *S* sehr einfach bestimmen, indem der Schnittpunkt der Geraden mit der x-Achse abgelesen wird ($t_0$). Dieser befindet sich sich zwar meistens auf der linken Seite des kritischen *t*'s, was aber in diesem Fall kein Problem darstellt. \\\n",
    "Wenn also der Schnittpunkt mit der x-Achse bekannt ist, ergibt sich:\n",
    "\n",
    "<br>\n",
    "<center>$\\large s = 0 = \\frac{Q}{4\\pi T}log\\frac{2.25Tt_0}{r^2S} $</center>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Da der Faktor $\\frac{Q}{4\\pi T}$ nur aus bekannten Konstanten aufgebaut ist, kann er nicht zu *0* werden. Deshalb kann er vernachlässigt werden und es ergibt sich:\n",
    "\n",
    "<br>\n",
    "<center>$\\large log\\frac{2.25Tt_0}{r^2S} = 0 \\qquad \\Rightarrow \\qquad \\frac{2.25Tt_0}{r^2S} = 1 $</center>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nach Umstellung ergibt sich eine Formel, durch die sich nach Ablesen des Zeitpunktes *t_0* (dem Schnittpunkt der Geraden mit der x-Achse) einfach die Transmissivität errechnen lässt:\n",
    "\n",
    "<br>\n",
    "<center>$\\large S = \\frac{2.25Tt_0}{r^2} $</center>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Überprüfung des kritischen Zeipunktes $t_{min}$__\\\n",
    "Schlussendlich muss bei jeder händischen Bestimmung von *S* und *T* überprüft werden, ob die für die Bestimmung der logarithmischen Geraden verwendeten Messwerte allesamt oberhalb des kritischen Zeitpunktes $t_min$ liegen, bei dem gilt: $u \\leq 0.02 $. \\\n",
    "Dies erfolgt durch einfaches Einsetzen der bestimmten Werte in die schon bekannte Gleichung:\\\n",
    "\\\n",
    "<br>\n",
    "<center>$\\large t_{min} = \\frac{r^2 S}{0.08 T} $</center>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Falls bei der händischen Festlegung der logarithmischen Geraden Messwerte berücksichtigt wurden, die kleiner als das errechnete $t_{min}$ sind, geht das ganze Spiel nochmal von vorne los... Das ist insbesondere dann von Belang, wenn die Bearbeitung mit Tabellenkalkulationsprogrammen (Excel, LibreOffice Calc, etc.) gearbeitet wird."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=blue> <em>Aufgabe 3: Rechnen Sie mithilfe der oben dargestellten Formeln und manuell abgelesenen Werten die Transmissivität, den Speicherkoeffizienten und den kritischen Zeitpunkt aus. Sind Ihre Ergebnisse identisch mit denen, die Sie oben eingestellt haben?</font> </em>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}