Theis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "eb507583-4c7a-42f2-a957-556b0e35a2b9",
   "metadata": {},
   "source": [
    "# Theis Pumping Test Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc105bea-1927-4724-9ac8-9541cb436f48",
   "metadata": {},
   "source": [
    "This notebook will solve the equations for pumping test analysis in a confined aquifer. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73aacb35-c4cd-4269-9aac-ff65622b4247",
   "metadata": {},
   "source": [
    "After the definition of variables, the value u is computed:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f82d478-6a5f-41af-97c3-a054256d5aa6",
   "metadata": {},
   "source": [
    "$$ u = \\frac{r^2 \\cdot S} {4 T t} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc13fb62-299a-4955-842b-94f1cf34cdd3",
   "metadata": {},
   "source": [
    "Python libraries for operating system functions, for numerical computation, for scientific constants and functions, for mathematical functions and for plotting are loaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 323,
   "id": "27bb613e-87d4-4a84-ae8a-2ad36c50e309",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy\n",
    "import scipy\n",
    "import math\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a57114-9bd3-469b-8b89-08819cd60697",
   "metadata": {},
   "source": [
    "The parameters for the pumping test are loaded: $h_=$ is the water level at the start, here 10 m. $Q_p$ is the pumping rate in $\\frac{m^3}{s}$, note per second to stay within SI units (!), $S_s$ is the specific storativity (default is $10^{-4}$), $k_f$ is the hydraulic conductivity in $m/s$, $H_a$ is the aquifer thickness in m, $T_a$ the transmissivity calculated as $H_a \\cdot k_f$. Finally, we define a range of values for the radius from the borehole in m (avoiding the value 0 as the drawdown is not defined at r=0) and we define a time range in seconds (to stay within SI units) from 1 minute (60 seconds) to a full day (86400 seconds). \n",
    "\n",
    "Comments are added with #."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 428,
   "id": "8377f25d-b180-45fc-8070-a6453766721d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "h0 = 10 # water level at start given in meters\n",
    "Qp = 0.02 # pumpint rate in m³/s\n",
    "Ss = 1.4E-4 # specific storage \n",
    "kf = 8E-5 # hydraulic conductivity in m/s\n",
    "Ha = 60 # thickness of the aquifer in m\n",
    "Ta = kf*H # transmissivity\n",
    "rb = numpy.arange(1, 100, 1) # distance from the borehole\n",
    "ts = numpy.arange(60, 86401, 60) # time after start of pumping test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d810db-1f47-4cff-8715-cd3d7972b012",
   "metadata": {},
   "source": [
    "According to theory, we first calculate $u$. Note, that here we take the element that has the index 20 (the index starts with 0, hence, it is element number 21) for $r_b$ and the element with index 60 of time $t_s$, again the sixty-first (as we start from 0).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 429,
   "id": "220a21ae-873b-48b8-af54-254c0e0f6214",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "u = rb[20]**2*Ss/(4*Ta*ts[60])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ff4153e-ed5e-47fc-a8b4-c47c749a9d12",
   "metadata": {},
   "source": [
    "After the value of $u$ has been calculated for a given combination of parameters, the value of $u$ is inserted into the well function. Note that the function for the natural logarithm is log2() from the math library and that factorial is calculated with factorial() from the same library:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc021c5e-6edb-4c03-915d-87cfc26c099e",
   "metadata": {},
   "source": [
    "$$ W(u) = -0.5772-ln(u)+u-\\frac{u^2}{2 \\cdot 2!}+\\frac{u^3}{3\\cdot 3!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b3107a4-1f40-4670-b194-4e7260c5b643",
   "metadata": {},
   "source": [
    "In the code, an explicit line break \\ is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 430,
   "id": "da319d2b-f6e4-4d38-b219-cc2641b3821f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wu=-0.5772-math.log2(u)+u \\\n",
    "   -u**2/(2*math.factorial(2)) \\\n",
    "   +u**3/(3*math.factorial(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d28c885e-807a-4a7f-93fe-aa306586aa3b",
   "metadata": {
    "tags": []
   },
   "source": [
    "Now, the drawdown can be calculated simply by using the function:\n",
    "$$s = \\frac{Q}{4 \\cdot \\pi \\cdot T} \\cdot W(u)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 431,
   "id": "096230c4-4890-4042-9a69-c8a53125981b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "s = Qp/(4*math.pi*Ta)*wu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2487c85b-0ae3-4fa6-ab63-1d22a96629e1",
   "metadata": {},
   "source": [
    "Finally, the computed drawdown is abstracted from the initial water level $h_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bf728c6-a252-46d9-b9ec-f7599bbc1f2c",
   "metadata": {},
   "source": [
    "$h = h_0-s$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12114556-6270-464d-8a80-21ba99e63e32",
   "metadata": {},
   "source": [
    "We will now make the program more modular by introducing functions for the calculation of u (called fu), w (called (fw), s (called fs) and h (falled fh). Finally, all funcions are integrated into a function Theis (called ftheis). Please note that a function is defined with **def**, then the function name that can be any string that is not a python key word or that has not been used before as a variable name. In the brackets, there are the variables needed in the function, please note, some rules. You should keep the order, in which they are listed. When you call the function, you transfer parameters into the function. You can define default values (such as Ha=100) in the definition, if you want, the the variable does not have to be called. If you do so, the variables without default values must come first (like the pumping rate $Q_p$ in the function ftheis(Qp). Since in this function, all other variables have default values, you could call the function with ftheis(0.0138), but you could also call it with ftheis(0.0138,10,100,1E-4,10,1E-4,3600), then the drawdown is given for a pumping rate of 0.0138 m³/s, an initial level of 10 m, an aquifer thickness of 100 m, a hydraulic conductivity of 1E-4, at a distance of 10 m, with a Specific Storativiy of 1E-4 and after 3600 seconds.  \n",
    "\n",
    "Values are explicitely returned with the command **return**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 432,
   "id": "884df16f-96a2-4e58-a9ee-57cbc8070dbf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fu(Ha,kf,rb,Ss,ts):\n",
    "    Ta = Ha*kf\n",
    "    u = rb**2*Ss/(4*Ta*ts)\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 433,
   "id": "b6dfcbb8-e61c-4af7-9d48-888873a215ee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fw(u):\n",
    "    w=-0.5772-math.log2(u)+u-u**2/(2*math.factorial(2))+u**3/(3*math.factorial(3))\n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 434,
   "id": "e4887602-8579-48b7-856c-bb2eaf54b57d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fs(w):\n",
    "    s = Qp/(4*math.pi*Ta)*w\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 435,
   "id": "2e52a49c-d2fd-4814-b58b-f72dbb3af633",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def fh(h0,s):\n",
    "    h = h0-s\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 436,
   "id": "5ead8c16-adb2-4c42-b138-bdd13a09b208",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ftheis(Qp,h0=0,Ha=100,kf=1E-4,rb=10,Ss=1E-4,ts=3600):\n",
    "    u=fu(Ha,kf,rb,Ss,ts)\n",
    "    w=fw(u)\n",
    "    s=fs(w)\n",
    "    h=fh(h0,s)\n",
    "    return h"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "344e70a3-4a6e-4fcc-9e0e-3075918fbc3d",
   "metadata": {},
   "source": [
    "Finally, we write the code for producing the table that is similar to the Excel file. We measure the length of the array with radii and call it lr, the same for the time. We now initialize empty tables to receive the data. What is new: We also define a counter cr, and ct and initialize it. Now the for-loops are defined, using ri as the i-th element of the array rb and ti as the i-th element of the array ts. Inside the loop, it is checked whether are are within the indices. We use the functions fu, fw, fs, and fh as well as ftheis, that have been defined before. The values u,w,s, and h are assigned to the arrays mu, mw, ms, mh, and mht at the locations given by the indices cr and ct. I used the functions one after the other and I used the integral function ftheis - just one way is necessary. We stop the loop with the command **break** if we should exceed it by any chance. Finally, one we got out of the first loop, we have to reset the counter for the inner for-loop with **ct=0** and we must increment the counter for the outer loop with **cr=cr+1**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 437,
   "id": "ee7b9225-72c2-4c1a-94e6-4ece69a7be99",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "lr = len(rb)\n",
    "lt = len(ts)\n",
    "mu = numpy.zeros((lr,lt), dtype=numpy.float64)\n",
    "mw = numpy.zeros((lr,lt), dtype=numpy.float64)\n",
    "ms = numpy.zeros((lr,lt), dtype=numpy.float64)\n",
    "mh = numpy.zeros((lr,lt), dtype=numpy.float64)\n",
    "mht = numpy.zeros((lr,lt), dtype=numpy.float64)\n",
    "cr = 0\n",
    "ct = 0\n",
    "for ri in rb:\n",
    "    for tj in ts:\n",
    "        if cr <= lr-1 and ct <= lt-1:\n",
    "            u = fu(Ha,kf,ri,Ss,tj)\n",
    "            w = fw(u)\n",
    "            s = fs(w)\n",
    "            h = fh(h0,s)\n",
    "            mu[cr,ct] = u\n",
    "            mw[cr,ct] = w\n",
    "            ms[cr,ct] = s\n",
    "            mh[cr,ct] = h\n",
    "            mht[cr,ct] = ftheis(Qp,h0,Ha,kf,ri,Ss,tj)\n",
    "            ct = ct+1\n",
    "        else: \n",
    "            break # the break command exits the for-loop \n",
    "    ct = 0  # we leave the inner loop and need to reset the counter ct\n",
    "    cr = cr+1 # we increase the counter for the distance cr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e49e17cb-d932-49e0-8e2c-b5ef023188df",
   "metadata": {},
   "source": [
    "Finally, we plot the data with the library **matplotlib** that contains a function **pyplot**, that has been imported and labelled with a shorter expression **plt**. We can therefore call the function **plt.plot**, it requires a time series. We plot the water levels in the matrix **mh** at the distance of 10 m (first index i in m[**i**,j], and we use all the times from the first to 1430 with the index range m[10,**1:1439**]. Please note, the times are listed with a time step of 60 seconds, hence, the correspond to minutes. At the location of the index, we find the water level **h**, that has been computed.The graph gets an xlabel and a ylabel and - of course - we shot it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 438,
   "id": "3494f5a8-6eba-4a5c-9834-3804089cdcd4",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAG2CAYAAACKxwc0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2xUlEQVR4nO3deXjTZb7//1eaJt1oAwVaKJTNhUWUVRF1WNRBXGE8KqO4oKMXHkFBZESOembcYByXUfQ4rj9wGQePR3H8uWEdCoIgIgURUTaBIhbK1r1N2uT+/pE2NLSFpiQ0H/p8XFeu5HN/tnfvK5oX92ezGWOMAAAALCqmuQsAAAA4FoQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgaYQZAABgac0aZr788ktdfvnlysjIkM1m0wcffBA03xijP//5z8rIyFBCQoJGjBihH374oXmKBQAAUalZw0xpaan69eun559/vt75f/3rX/X000/r+eef16pVq9ShQwf99re/VXFx8XGuFAAARCtbtDxo0mazacGCBRo7dqwk/6hMRkaGpk6dqhkzZkiS3G630tPT9fjjj2vixInNWC0AAIgWsc1dQEO2bdum3bt3a9SoUYG2uLg4DR8+XMuXL28wzLjdbrnd7sC0z+fTgQMH1LZtW9lstojXDQAAjp0xRsXFxcrIyFBMzJEPJEVtmNm9e7ckKT09Pag9PT1dO3bsaHC92bNn66GHHopobQAA4PjYuXOnOnfufMRlojbM1Dh8NMUYc8QRlpkzZ2ratGmB6cLCQnXp0kU7d+5USkpKxOoEAADhU1RUpMzMTCUnJx912agNMx06dJDkH6Hp2LFjoD0/P7/OaE1tcXFxiouLq9OekpJCmAEAwGIac4pI1N5npnv37urQoYOysrICbR6PR0uWLNE555zTjJUBAIBo0qwjMyUlJdqyZUtgetu2bVq7dq1SU1PVpUsXTZ06VbNmzdIpp5yiU045RbNmzVJiYqKuu+66ZqwaAABEk2YNM99++61GjhwZmK451+Wmm27SvHnzdO+996q8vFx33HGHDh48qCFDhujzzz9v1PEzAADQMkTNfWYipaioSC6XS4WFhZwzAwCARYTy+x2158wAAAA0BmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYGmEGAABYWosJM+4qb3OXAAAAIqDFhJn9JZ7mLgEAAERAiwkz+0rczV0CAACIgBYUZhiZAQDgRNSCwgwjMwAAnIhaTJjZX8zIDAAAJ6IWE2b2lVY0dwkAACACWk6Y4TATAAAnpBYTZvIKCDMAAJyIWkyY2b6/RMaY5i4DAACEWYsJM2Uen/IKOW8GAIATTYsJM5K0Jb+kuUsAAABh1qLCzKY9xc1dAgAACLMWFWbW7Cxo7hIAAECYtagwk7PjYHOXAAAAwqzFhBl7jE15hRXaeaCsuUsBAABh1GLCTP/OrSVJi37Kb95CAABAWLWYMDOyV3tJUtaGPc1cCQAACKcWE2ZG9EqTJH39834Vllc2czUAACBcWkyY6dY2ST3Tk1XlM/rX2l3NXQ4AAAiTFhNmJOn3Z2VKkt5emcujDQAAOEG0qDBz5YDOinfE6KfdxVqxdX9zlwMAAMKgRYUZV6JD1wz2j8488+/NjM4AAHACiOowU1VVpQceeEDdu3dXQkKCevTooYcfflg+n6/J27xjxMlyxsbom20HuEwbAIATQFSHmccff1wvvviinn/+ef3444/661//qieeeELPPfdck7fZwRWvm8/tJkn604c/qNzjDVO1AACgOUR1mFmxYoXGjBmjSy+9VN26ddNVV12lUaNG6dtvvz2m7d51/inKcMXrl4PleurzjWGqFgAANIeoDjPnnXee/v3vf2vTpk2SpO+++07Lli3TJZdc0uA6brdbRUVFQa/DJcXF6pGxfSVJry7bpkU/cSM9AACsKqrDzIwZM3TttdeqV69ecjgcGjBggKZOnaprr722wXVmz54tl8sVeGVmZta73AW90zXhnG6SpLvf+U5b8osj8ScAAIAIi+ow88477+itt97S22+/rZycHL3++ut68skn9frrrze4zsyZM1VYWBh47dy5s+FlL+mlAV1aq7C8Ujf9f6u0p6giEn8GAACIIJuJ4uuTMzMzdd9992nSpEmBtkcffVRvvfWWfvrpp0Zto6ioSC6XS4WFhUpJSakz/0CpR//x9+Xatq9UPdol6a1bhyijdULY/gYAABC6o/1+1xbVIzNlZWWKiQku0W63H9Ol2YdLTXLqjVvOUqfWCfp5X6mufnGFft5bErbtAwCAyIrqMHP55Zfrscce08cff6zt27drwYIFevrpp/W73/0urPvJTE3Uu7cPVY/2SdpVUK6x//OVsrkHDQAAlhDVh5mKi4v14IMPasGCBcrPz1dGRoauvfZa/fd//7ecTmejthHKMNX+Erdue+Nb5eQWyGaTJo88WXddcIoc9qjOfAAAnHBC+f2O6jATDqF0hiR5qnx66P//Qf9YmStJOi0jRU9f0189OyRHulQAAFDthDlnpjk4Y2P02O9O15xrB8iV4NAPvxbpsueWavanP6rEXdXc5QEAgMMQZhpwRb8MZd09TBf2TlOl1+ilJT9r5JOL9b+rdqrKG74TkAEAwLHhMNNRGGO06Kd8PfLRBm3fXyZJ6to2UXeef4rG9s9QLOfTAAAQdpwzU8uxhpka7iqv3li+Q39fslUHSj2SpC6piZpwTjddPbizkuMd4SoZAIAWjzBTS7jCTI1Sd5Xe/HqHXv7y50CoSXLadfXgTF1/dhednMaJwgAAHCvCTC3hDjM1yjxVej9nl+Yt364t+Ydustcvs7WuGtRZl5/RUa0TG3f5OAAACEaYqSVSYaaGMUZfbdmvecu3K3tjvrw+f3c67TEa2au9Lu7bUef3TlMKh6EAAGg0wkwtkQ4zte0tdutfa3fpvZxd+jGvKNDusNt07sntNPq0Djq/d5rSkuMjWgcAAFZHmKnleIaZ2jb8WqRP1+fp0/W7gw5DSVLvjikafmp7DTu1nQZ3TZUzliuiAACojTBTS3OFmdq25Bdr4Q97tPCH3Vr3S2HQvESnXWf3aKuze6TqrO5tdVpGCo9PAAC0eISZWqIhzNS2r8StZZv36ctNe/Xl5n3aV+IOmp/gsGtQ1zY6q3uqzuyWqn6ZLiU6Y5upWgAAmgdhppZoCzO1+XxGG/KKtHzrPn2z7aBWbT+gwvLKoGVibNKp6ck6o7NL/TJbq1/n1urZIZnRGwDACY0wU0s0h5nD+XxGm/NL9M22/fpm+0Gt2nZAu4sq6iwXFxuj0zJSdFqGS707pqhXx2T1TE9WUhwjOACAEwNhphYrhZn67Cmq0Hc7C7Tul0J990uBvttZoKKKug+8tNmkrqmJ6tXBH256d0xRz/RkZaYmyh5ja4bKAQBoOsJMLVYPM4czxmj7/jKt+6VAG/KK9GNesX7KK1J+sbve5Z32GHVvl6ST0pJ0cvtWOimtlU5q30o92idxLg4AIGoRZmo50cJMQ/aXuPXT7mL9WB1wfswr0ta9JXJXNfyE706tE9SjfZK6t0tSl9REdUlNVNe2/s8JTvtxrB4AgGCEmVpaSpipj9dn9GtBubbsLdHW/BJt3VuiLfkl2rq3NPBcqYakJcepa9tEdUlNUte2idWfE9W5TaLatXLKZuPQFQAgcggztbTkMHMkB0o92lodcnYcKFPu/jLtOFCqHfvLVFzPOTm1OWNjlOGKV6c2CcpwJSijdYI6tUlQp9b+zx1d8Yp3MLIDAGi6UH6/OWmihUpNcio1yX8vm9qMMSooq9SOA2Xasb+0OuQcCjv5xW55qnzavr9M2/eXNbj9dq3i1Kl1vDJaJyg9JV4dXPFKT4lTenK80qqnW3H1FQAgDPg1QRCbzaY2SU61SXKqf2brOvM9VT7tKarQLwfL9WtB9auwvNZ0hcorvdpX4ta+Ere+O+yOx7UlOe1KT4mvfsUpPaU66FRPpyXHq12ykxOVAQBHxK8EQuKMjVFmaqIyUxPrnV8zsrOroFy7qsPOniK38osqtKe4QrsLK5Rf5Faxu0qlHq9+3leqn/eVHnGfiU672rWKU7tWTv97cpzaJTn9763iDs1LjlNyXCzn8wBAC0OYQVjVHtnp28nV4HKl7irlF7v94aa4QnuKKrS70K09xRX+4FPk1p6iCrmrfCrzeJV7oEy5Bxo+rFXDGRsTFHTaJjmVWl1Pm0SH2iTWnnbKleDgPjwAYHGEGTSLpLhYdY+LVfd2SQ0uY4xRibtK+0s8gcNWe0s82lfs1v5St/YVH2rfV+JRibtKniqffi2s0K+Fde+cXB+bTWqd4AiEG3/YOTSdmugPPqlJDrWuDj+uBAePkwCAKEKYQdSy2WxKjncoOd6hbkcIPTUqAufq+APPvhK39pd6VFDm0YHSSh0s8+hAqUcHyzw6WOpRUUWVjJEOllXqYFmlpCMf7qotyWmXK8GhlOpw0zrREQg6rgSHXLWCT82rdfXyjAQBQHgRZnDCiHfY1bmN/144jVHp9amg7FDIqRN6qoPPgbJK/+dSj4rd/svWSz1elXq8jR4Bqi05LrZOCEqOj60ObofeU+ppS46PVVwsl70DQG2EGbRYDnuM2ifHqX1yXKPXqfL6VFxRpcLyShWWV6qg+r2wvFJFNW1lnkBbYXlVoL2kOggVu6tU7K7SroLyJtXttMcEgk3tkFM7+KRUt7WKC57fKi5WSXF2JTljFcMIEYATBGEGCEGsPSZwgnOoKr2+QLA5/FVcUaWiCv+7/1UZ9F5SURUYFfJ4fdpf6tH+o9zF+WgSnXYlxcUGBRz/59jq9trzY5XotAc+H1rO35bgsHMVGYBmQ5gBjhOHPUZtW8WpbavGjwTV5vMZlXjqDztF1Z9L6oSh2iGpUqUer7w+/02/yzxelXm82tvAQ0pDEWOTkpzBAacmFCU57UqMi1Wiw65Ep10JTv8yCQ67Ep3+kOR/xSqh+nNS9WdnLCdaAzg6wgxgETExNqXEO5QS75CU0KRtGGPkrvKpxF2lUndV9bu31udDbWWe2m3+ZUprtZW6vSr1+E+i9plDh8/CKTbGFgg6/iB0KOjUnq75XHvZQHAKzPOPICU47Ip3xshpj2E0CThBEGaAFsRmsyneYVe8w38jwmPl8xmVV3qDQlAg7FQHnzK3t3oUqCowGlRe6V+23ONVWaW/vdzj3055pVeVXv/oUZXPqKh65CncbDYdCjcOu+IdMUpw1p62B+YnOO2Kc8QETcfH2hXvtNfaRox/nVrbSHDYFRcbw/lJQIQRZgA0WUyMLXA4KS2M2630+oICULmnnkBU/bm01udyj3+0qPY6NWGr3ONVRdWhoGTMoUNtkRYIOrXCU4KzOkBVT8fF+qfrvvvXjXPEKD42+L32srWnGXVCS0OYARB1HPYYuRJi5EpwhH3blV6fKiq9qqj0v5dX+oNPeaW3ur2m7dD8iqBlgterqKp+D2rzyVPlC+zTv45PBaoM+99TH5tNioutLxjVCj81ISk2RnGO4OkjvgeFqRg5q/fjjPWHKIfdRpDCcUeYAdCiOOwxcthjlBwf2f14fSYo4Lir/AGpvPKw4FQdktxV/pBU37v7qPN9qqjyyvgHnWTMoQBV2LQ7ADSZzea/fUBNyPGHKv+0s/Zne60QVGcZ/3pOe4ziHLXf7fUsWz1ttweWrWmP5U7dLQZhBgAiwF7rENzxYIxRpdeoosrrDzeNDETuKv9oU933o2/DXemVx+sLHLrz16Hq7fpUrPCf6xSKGJuCRo6CgpTDrjh73SDltMfIEXsoFDnsNjntdjlibbXaai9nq6et1nKBbRxqi41h9CrcCDMAcAKw2Wxyxvp/WBXhUafD+XxGHq8vEI481WHGU/0KfPb6g9ahZWsv462zvLvKH5Zq2txB8/3zarZX015z6wHJf5VdzQhVNLHZFAg/NWGnJvjUDj0Ou03OWHtQYKq9nKN2kAq0xRzWZqunLTho1Wy39merPXaFMAMAOCYxMTbFx/hPZJbCf55TKLw+ExSO3IeHoypfPSGoVgDz+lRZZVTpPTTf3+ar1VY9P6jt0OfKKhNYx1PdZg5lLBmjQHDTsd/mKSJibP6bhDrt9Qce/7yaz7ZAOKvvc01wio2x1fnsiI2Ro4HPnrKSRtdLmAEAnDDsMTb/5fHO6HmGmTFGXp//MKCnVsCpCUjumiBU5T9kV7et9nLBQcpda37N+nXbfPJ4jTzVV/N5am/b51+n9oiW5B/VCgSuZuJzlzV6WcIMAAARZLPZFGu3KdauqApZtfl8JhBsKmuFnMpAMDKBgBT8+UjzGvhcZepu31f3c3mpTTsbWT9hBgCAFi4mxqa4GLviYiUd+/00w6KoqEiu+xu3LNetAQAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASwv52UwVFRV67rnnlJ2drfz8fPl8wU/UzMnJCVtxAAAARxNymLnllluUlZWlq666SmeddZZsNlsk6gIAAGiUkMPMxx9/rE8++UTnnntuJOoBAAAIScjnzHTq1EnJycmRqAUAACBkIYeZp556SjNmzNCOHTsiUQ8AAEBIQj7MNHjwYFVUVKhHjx5KTEyUw+EImn/gwIGwFQcAAHA0IYeZa6+9Vrt27dKsWbOUnp7OCcAAAKBZhRxmli9frhUrVqhfv36RqAcAACAkIZ8z06tXL5WXl0eiFgAAgJCFHGb+8pe/6J577tHixYu1f/9+FRUVBb0AAACOJ5sxxoSyQkyMP/8cfq6MMUY2m01erzd81YVBUVGRXC6XCgsLlZKS0tzlAACARgjl9zvkc2ays7ObXBgAAEC4hRxmhg8fHok6AAAAmiTqn5q9a9cuXX/99Wrbtq0SExPVv39/rV69urnLAgAAUSLkkZnj6eDBgzr33HM1cuRIffrpp0pLS9PWrVvVunXr5i4NAABEiagOM48//rgyMzM1d+7cQFu3bt2aryAAABB1ovow04cffqjBgwfr6quvVlpamgYMGKBXXnnliOu43W4uFwcAoAWJ6jDz888/6+9//7tOOeUULVy4ULfffrvuuusuvfHGGw2uM3v2bLlcrsArMzPzOFYMAACOt5DvM7Nnzx5Nnz5d//73v5Wfn6/DVw/nfWacTqcGDx6s5cuXB9ruuusurVq1SitWrKh3HbfbLbfbHZguKipSZmYm95kBAMBCInqfmQkTJig3N1cPPvigOnbsGNEHTXbs2FF9+vQJauvdu7fee++9BteJi4tTXFxcxGoCAADRJeQws2zZMi1dulT9+/ePQDnBzj33XG3cuDGobdOmTeratWvE9w0AAKwh5HNmMjMz6xxaipS7775bX3/9tWbNmqUtW7bo7bff1ssvv6xJkyYdl/0DAIDoF3KYeeaZZ3Tfffdp+/btESgn2JlnnqkFCxbon//8p/r27atHHnlEzzzzjMaPHx/xfQMAAGsI+QTgNm3aqKysTFVVVUpMTJTD4Qiaf+DAgbAWeKx40CQAANYT0ROAn3nmmabWBQAAEHYhh5mbbropEnUAAAA0SaPCTFFRUWCI52h31OVQDgAAOJ4aFWbatGmjvLw8paWlqXXr1vXeW8YYI5vNFtab5gEAABxNo8LMokWLlJqaKknKzs6OaEEAAAChCPlqJqvhaiYAAKwnlN/vqH7QJAAAwNEQZgAAgKURZgAAgKURZgAAgKURZgAAgKU16tLsAQMG1Htvmfrk5OQcU0EAAAChaFSYGTt2bITLAAAAaBruMwMAAKJOxO8zU1BQoFdffVUzZ87UgQMHJPkPL+3ataspmwMAAGiykJ+avW7dOl144YVyuVzavn27brvtNqWmpmrBggXasWOH3njjjUjUCQAAUK+QR2amTZumCRMmaPPmzYqPjw+0X3zxxfryyy/DWhwAAMDRhBxmVq1apYkTJ9Zp79Spk3bv3h2WogAAABor5DATHx+voqKiOu0bN25U+/btw1IUAABAY4UcZsaMGaOHH35YlZWVkiSbzabc3Fzdd999+o//+I+wFwgAAHAkIYeZJ598Unv37lVaWprKy8s1fPhwnXzyyUpOTtZjjz0WiRoBAAAaFPLVTCkpKVq2bJkWLVqknJwc+Xw+DRw4UBdeeGEk6gMAADiikG+at337dnXr1i1C5YQfN80DAMB6InrTvB49eui8887TSy+9FLhhHgAAQHMJOcx8++23Gjp0qB599FFlZGRozJgxevfdd+V2uyNRHwAAwBGFHGYGDhyoJ554Qrm5ufr000+VlpamiRMnKi0tTbfcckskagQAAGhQWB40mZOToz/84Q9at26dvF5vOOoKG86ZAQDAeiL+oElJ2rlzp/7617+qf//+OvPMM5WUlKTnn3++qZsDAABokpAvzX755Zf1j3/8Q1999ZV69uyp8ePH64MPPrDUFU4AAODEEXKYeeSRR/T73/9ezz77rPr37x+BkgAAABov5DCTm5srm80WiVoAAABCFvI5MzabTUuXLtX111+voUOHateuXZKkN998U8uWLQt7gQAAAEcScph57733dNFFFykhIUFr1qwJ3F+muLhYs2bNCnuBAAAARxJymHn00Uf14osv6pVXXpHD4Qi0n3POOcrJyQlrcQAAAEcTcpjZuHGjhg0bVqc9JSVFBQUF4agJAACg0UIOMx07dtSWLVvqtC9btkw9evQIS1EAAACNFXKYmThxoqZMmaKVK1fKZrPp119/1T/+8Q9Nnz5dd9xxRyRqBAAAaFDIl2bfe++9Kiws1MiRI1VRUaFhw4YpLi5O06dP1+TJkyNRIwAAQIOa/GymsrIybdiwQT6fT3369FGrVq3CXVtY8GwmAACsJ5Tf75BHZmokJiZq8ODBTV0dAAAgLBoVZq688spGb/D9999vcjEAAAChalSYcblcka4DAACgSRoVZubOnRvpOgAAAJok5EuzAQAAoglhBgAAWBphBgAAWBphBgAAWFpIYaayslIjR47Upk2bIlUPAABASEIKMw6HQ+vXr5fNZotUPQAAACEJ+TDTjTfeqNdeey0StQAAAIQs5McZeDwevfrqq8rKytLgwYOVlJQUNP/pp58OW3EAAABHE3KYWb9+vQYOHChJdc6d4fATAAA43kIOM9nZ2ZGoAwAAoEmafGn2li1btHDhQpWXl0uSjDFhKwoAAKCxQg4z+/fv1wUXXKBTTz1Vl1xyifLy8iRJt956q+65556wFwgAAHAkIYeZu+++Ww6HQ7m5uUpMTAy0jxs3Tp999llYiwMAADiakM+Z+fzzz7Vw4UJ17tw5qP2UU07Rjh07wlYYAABAY4Q8MlNaWho0IlNj3759iouLC0tRAAAAjRVymBk2bJjeeOONwLTNZpPP59MTTzyhkSNHhrU4AACAown5MNMTTzyhESNG6Ntvv5XH49G9996rH374QQcOHNBXX30ViRoBAAAaFPLITJ8+fbRu3TqdddZZ+u1vf6vS0lJdeeWVWrNmjU466aRI1AgAANAgmwnxBjG5ubnKzMys926/ubm56tKlS9iKC4eioiK5XC4VFhYqJSWlucsBAACNEMrvd8gjM927d9fevXvrtO/fv1/du3cPdXMAAADHJOQwY4ypd1SmpKRE8fHxYSmqIbNnz5bNZtPUqVMjuh8AAGAdjT4BeNq0aZL8Vy89+OCDQZdne71erVy5Uv379w97gTVWrVqll19+WWeccUbE9gEAAKyn0WFmzZo1kvwjM99//72cTmdgntPpVL9+/TR9+vTwVyj/qM/48eP1yiuv6NFHH43IPgAAgDU1OszUPC375ptv1rPPPntcT6adNGmSLr30Ul144YVHDTNut1tutzswXVRUFOnyAABAMwr5PjNz586NRB0Nmj9/vnJycrRq1apGLT979mw99NBDEa4KAABEi5DDjOQ/f+Xdd99Vbm6uPB5P0Lz3338/LIVJ0s6dOzVlyhR9/vnnjT65eObMmYHzeyT/yExmZmbYagIAANEl5KuZ5s+fr3PPPVcbNmzQggULVFlZqQ0bNmjRokVyuVxhLW716tXKz8/XoEGDFBsbq9jYWC1ZskRz5sxRbGysvF5vnXXi4uKUkpIS9AIAACeukEdmZs2apb/97W+aNGmSkpOT9eyzz6p79+6aOHGiOnbsGNbiLrjgAn3//fdBbTfffLN69eqlGTNmyG63h3V/AADAekIOM1u3btWll14qyT8KUlpaKpvNprvvvlvnn39+WM9XSU5OVt++fYPakpKS1LZt2zrtAACgZQr5MFNqaqqKi4slSZ06ddL69eslSQUFBSorKwtvdQAAAEcR8sjMb37zG2VlZen000/XNddcoylTpmjRokXKysrSBRdcEIkagyxevDji+wAAANYRcph5/vnnVVFRIcl/5ZDD4dCyZct05ZVX6sEHHwx7gQAAAEcS8lOzrYanZgMAYD0RfWp2zWMFNm3a1OQCAQAAwiXkMNOqVSs99dRT6tWrlzIyMnTttdfqxRdf1E8//RSJ+gAAAI6oyYeZdu/ercWLF2vx4sVasmSJNm3apLS0NOXl5YW7xmPCYSYAAKwnooeZaiQnJ6tNmzZq06aNWrdurdjYWHXo0KGpmwMAAGiSkMPMjBkzdPbZZ6tdu3Z64IEH5PF4NHPmTO3Zs0dr1qyJRI0AAAANCvkwU0xMjNq3b6+7775bY8aMUe/evSNVW1hwmAkAAOsJ5fc75PvMrFmzRkuWLNHixYv11FNPyW63a/jw4RoxYoRGjBgR9eEGAACcWI75PjPfffednnnmGb311lvy+Xz1Psm6OTEyAwCA9UR0ZEbyj87UXMm0dOlSFRUVqX///ho5cmSTCgYAAGiqkMNMmzZtVFJSon79+mnEiBG67bbbNGzYMEY9AABAswg5zLz55puEFwAAEDVCDjOXXXZZJOoAAABokibfNA8AACAaEGYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClEWYAAIClRXWYmT17ts4880wlJycrLS1NY8eO1caNG5u7LAAAEEWiOswsWbJEkyZN0tdff62srCxVVVVp1KhRKi0tbe7SAABAlLAZY0xzF9FYe/fuVVpampYsWaJhw4Y1ap2ioiK5XC4VFhYqJSUlwhUCAIBwCOX3O/Y41RQWhYWFkqTU1NQGl3G73XK73YHpoqKiiNcFAACaT1QfZqrNGKNp06bpvPPOU9++fRtcbvbs2XK5XIFXZmbmcawSAAAcb5Y5zDRp0iR9/PHHWrZsmTp37tzgcvWNzGRmZnKYCQAACznhDjPdeeed+vDDD/Xll18eMchIUlxcnOLi4o5TZQAAoLlFdZgxxujOO+/UggULtHjxYnXv3r25SwIAAFEmqsPMpEmT9Pbbb+tf//qXkpOTtXv3bkmSy+VSQkJCM1cHAACiQVSfM2Oz2eptnzt3riZMmNCobXBpNgAA1nPCnDMTxTkLAABECctcmg0AAFAfwgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0wgwAALA0S4SZF154Qd27d1d8fLwGDRqkpUuXNndJAAAgSkR9mHnnnXc0depU3X///VqzZo1+85vf6OKLL1Zubm5zlwYAAKKAzRhjmruIIxkyZIgGDhyov//974G23r17a+zYsZo9e/ZR1y8qKpLL5VJhYaFSUlIiWSoAAAiTUH6/Y49TTU3i8Xi0evVq3XfffUHto0aN0vLly+tdx+12y+12B6YLCwsl+TsFAABYQ83vdmPGXKI6zOzbt09er1fp6elB7enp6dq9e3e968yePVsPPfRQnfbMzMyI1AgAACKnuLhYLpfriMtEdZipYbPZgqaNMXXaasycOVPTpk0LTBcUFKhr167Kzc09amfgyIqKipSZmamdO3dyyO4Y0I/hQ1+GD30ZHvRj+BhjVFxcrIyMjKMuG9Vhpl27drLb7XVGYfLz8+uM1tSIi4tTXFxcnXaXy8UXK0xSUlLoyzCgH8OHvgwf+jI86MfwaOwgRFRfzeR0OjVo0CBlZWUFtWdlZemcc85ppqoAAEA0ieqRGUmaNm2abrjhBg0ePFhDhw7Vyy+/rNzcXN1+++3NXRoAAIgCUR9mxo0bp/379+vhhx9WXl6e+vbtq08++URdu3Zt1PpxcXH605/+VO+hJ4SGvgwP+jF86MvwoS/Dg35sHlF/nxkAAIAjiepzZgAAAI6GMAMAACyNMAMAACyNMAMAACzthA4zL7zwgrp37674+HgNGjRIS5cube6Sosrs2bN15plnKjk5WWlpaRo7dqw2btzY4PITJ06UzWbTM888E9Tudrt15513ql27dkpKStIVV1yhX375JcLVR7fZs2fLZrNp6tSpgbaSkhJNnjxZnTt3VkJCgnr37h30AFWJvpSkbt26yWaz1XlNmjRJlZWVmjFjhk4//XQlJSUpIyNDN954o3799degbdCPfsXFxZo6daq6du2qhIQEnXPOOVq1alXQMj/++KOuuOIKuVwuJScn6+yzz1Zubm5gfkvsyy+//FKXX365MjIyZLPZ9MEHHwTmhfM7ePDgQd1www1yuVxyuVy64YYbVFBQcBz+whOQOUHNnz/fOBwO88orr5gNGzaYKVOmmKSkJLNjx47mLi1qXHTRRWbu3Llm/fr1Zu3atebSSy81Xbp0MSUlJXWWXbBggenXr5/JyMgwf/vb34Lm3X777aZTp04mKyvL5OTkmJEjR5p+/fqZqqqq4/SXRJdvvvnGdOvWzZxxxhlmypQpgfZbb73VnHTSSSY7O9ts27bNvPTSS8Zut5sPPvggsAx9aUx+fr7Jy8sLvLKysowkk52dbQoKCsyFF15o3nnnHfPTTz+ZFStWmCFDhphBgwYFbYN+9LvmmmtMnz59zJIlS8zmzZvNn/70J5OSkmJ++eUXY4wxW7ZsMampqeaPf/yjycnJMVu3bjUfffSR2bNnT2AbLbEvP/nkE3P//feb9957z0gyCxYsCMwL53dw9OjRpm/fvmb58uVm+fLlpm/fvuayyy47Xn/mCeWEDTNnnXWWuf3224PaevXqZe67775mqij65efnG0lmyZIlQe2//PKL6dSpk1m/fr3p2rVrUJgpKCgwDofDzJ8/P9C2a9cuExMTYz777LPjVXrUKC4uNqeccorJysoyw4cPDwozp512mnn44YeDlh84cKB54IEHjDH0ZUOmTJliTjrpJOPz+eqd/8033xhJgX+o0I9+ZWVlxm63m48++iiovV+/fub+++83xhgzbtw4c/311ze4DfrS1Akz9WnKd3DDhg1Gkvn6668Dy6xYscJIMj/99FP4/5AT3Al5mMnj8Wj16tUaNWpUUPuoUaO0fPnyZqoq+hUWFkqSUlNTA20+n0833HCD/vjHP+q0006rs87q1atVWVkZ1NcZGRnq27dvi+zrSZMm6dJLL9WFF15YZ955552nDz/8ULt27ZIxRtnZ2dq0aZMuuugiSfRlfTwej9566y3dcsstDT5ctrCwUDabTa1bt5ZEP9aoqqqS1+tVfHx8UHtCQoKWLVsmn8+njz/+WKeeeqouuugipaWlaciQIUGHVOjLxmnKd3DFihVyuVwaMmRIYJmzzz5bLpeLvm2CEzLM7Nu3T16vt87DKNPT0+s8tBJ+xhhNmzZN5513nvr27Rtof/zxxxUbG6u77rqr3vV2794tp9OpNm3aBLW3xL6eP3++cnJyNHv27Hrnz5kzR3369FHnzp3ldDo1evRovfDCCzrvvPMk0Zf1+eCDD1RQUKAJEybUO7+iokL33XefrrvuusBD/ehHv+TkZA0dOlSPPPKIfv31V3m9Xr311ltauXKl8vLylJ+fr5KSEv3lL3/R6NGj9fnnn+t3v/udrrzySi1ZskQSfdkYTf0O7t69W2lpaXW2l5aWRt82QdQ/zuBYHP4vOWNMg/+6a+kmT56sdevWadmyZYG21atX69lnn1VOTk7I/dbS+nrnzp2aMmWKPv/88zr/Eq4xZ84cff311/rwww/VtWtXffnll7rjjjvUsWPHekdyarS0vqzttdde08UXX6yMjIw68yorK/X73/9ePp9PL7zwwlG31RL78c0339Qtt9yiTp06yW63a+DAgbruuuuUk5Mjn88nSRozZozuvvtuSVL//v21fPlyvfjiixo+fHiD222JfVmfY/0O1teH9G3TnJAjM+3atZPdbq+TbvPz8+uM1kC688479eGHHyo7O1udO3cOtC9dulT5+fnq0qWLYmNjFRsbqx07duiee+5Rt27dJEkdOnSQx+PRwYMHg7bZ0vp69erVys/P16BBgwJ9tWTJEs2ZM0exsbEqLS3Vf/3Xf+npp5/W5ZdfrjPOOEOTJ0/WuHHj9OSTT0qiLw+3Y8cOffHFF7r11lvrzKusrNQ111yjbdu2KSsrK/AvYol+rO2kk07SkiVLVFJSop07d+qbb75RZWWlunfvrnbt2ik2NlZ9+vQJWqd3796Bq5noy4Yd63ewQ4cO2rNnT53t7t27t8X3bVOckGHG6XRq0KBBysrKCmrPysrSOeec00xVRR9jjCZPnqz3339fixYtUvfu3YPm33DDDVq3bp3Wrl0beGVkZOiPf/yjFi5cKEkaNGiQHA5HUF/n5eVp/fr1LaqvL7jgAn3//fdBfTV48GCNHz9ea9euldfrVWVlpWJigv+Ts9vtgX8h05fB5s6dq7S0NF166aVB7TU/Ips3b9YXX3yhtm3bBs2nH+tKSkpSx44ddfDgQS1cuFBjxoyR0+nUmWeeWed2DJs2bQo8yJe+rF84voNDhw5VYWGhvvnmm8AyK1euVGFhYYvu2yZrrjOPI63m0uzXXnvNbNiwwUydOtUkJSWZ7du3N3dpUeM///M/jcvlMosXLw66FLasrKzBdQ6/mskY/yWInTt3Nl988YXJyckx559//gl/6WZjHH410/Dhw81pp51msrOzzc8//2zmzp1r4uPjzQsvvBBYhr7083q9pkuXLmbGjBlB7ZWVleaKK64wnTt3NmvXrg363rrd7sBy9KPfZ599Zj799FPz888/m88//9z069fPnHXWWcbj8RhjjHn//feNw+EwL7/8stm8ebN57rnnjN1uN0uXLg1soyX2ZXFxsVmzZo1Zs2aNkWSefvpps2bNGrNjx46wfgdHjx5tzjjjDLNixQqzYsUKc/rpp3NpdhOdsGHGGGP+53/+x3Tt2tU4nU4zcODAOpcct3SS6n3NnTu3wXXqCzPl5eVm8uTJJjU11SQkJJjLLrvM5ObmRrZ4Czg8zOTl5ZkJEyaYjIwMEx8fb3r27GmeeuqpoEuO6Uu/hQsXGklm48aNQe3btm1r8HubnZ0dWI5+9HvnnXdMjx49jNPpNB06dDCTJk0yBQUFQcu89tpr5uSTTzbx8fGmX79+Qfc9MqZl9mV2dna937GbbroprN/B/fv3m/Hjx5vk5GSTnJxsxo8fbw4ePHh8/9gThM0YY47bMBAAAECYnZDnzAAAgJaDMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAMAACyNMAO0cIsXL5bNZlNBQcFx37fNZtMHH3xw3Pf75z//Wf379z/u+wUQGdw0D2hBRowYof79++uZZ54JtHk8Hh04cEDp6enH/Wm9u3fvVps2bRQXF3dc91tSUiK3213nmTrHYvv27erevbvWrFlDUAKOs9jmLgBA83I6nerQoUOz7Lu59tuqVSu1atWqWfYNIPw4zAS0EBMmTNCSJUv07LPPymazyWazafv27XUOM82bN0+tW7fWRx99pJ49eyoxMVFXXXWVSktL9frrr6tbt25q06aN7rzzTnm93sD2PR6P7r33XnXq1ElJSUkaMmSIFi9efMSaah9m2r59u2w2m95//32NHDlSiYmJ6tevn1asWHHUbbz00ku67LLLlJiYqN69e2vFihXasmWLRowYoaSkJA0dOlRbt24NrHP4YaYJEyZo7NixevLJJ9WxY0e1bdtWkyZNUmVlZb211mjdurXmzZsnSYGnzg8YMEA2m00jRowILDd37lz17t1b8fHx6tWrl1544YWgfps8ebI6duyo+Ph4devWTbNnzz7i3wwgGGEGaCGeffZZDR06VLfddpvy8vKUl5enzMzMepctKyvTnDlzNH/+fH322WdavHixrrzySn3yySf65JNP9Oabb+rll1/W//3f/wXWufnmm/XVV19p/vz5Wrduna6++mqNHj1amzdvDqnO+++/X9OnT9fatWt16qmn6tprr1VVVdUR13nkkUd04403au3aterVq5euu+46TZw4UTNnztS3334rSZo8efIRt5Gdna2tW7cqOztbr7/+uubNmxcIKo3xzTffSJK++OIL5eXl6f3335ckvfLKK7r//vv12GOP6ccff9SsWbP04IMP6vXXX5ckzZkzRx9++KH+93//Vxs3btRbb72lbt26NXq/AKQT+qnZAIId/iRvYw49Ibjmab1z5841ksyWLVsCy0ycONEkJiaa4uLiQNtFF11kJk6caIwxZsuWLcZms5ldu3YFbfuCCy4wM2fObLAeSWbBggXGmENPxH711VcD83/44Qcjyfz4449H3MYDDzwQmF6xYoWRZF577bVA2z//+U8THx8fmP7Tn/5k+vXrF5i+6aabTNeuXU1VVVWg7eqrrzbjxo2rt9YaLpcr8JT5mvrXrFkTtExmZqZ5++23g9oeeeQRM3ToUGOMMXfeeac5//zzg56eDiA0nDMDoI7ExESddNJJgen09HR169Yt6DyT9PR05efnS5JycnJkjNGpp54atJ2mnGR7xhlnBD537NhRkpSfn69evXo1ap309HRJ0umnnx7UVlFRoaKiIqWkpNS7jdNOO012uz1o399//31ItR9u79692rlzp/7whz/otttuC7RXVVXJ5XJJ8h/i+u1vf6uePXtq9OjRuuyyyzRq1Khj2i/Q0hBmANThcDiCpm02W71tPp9PkuTz+WS327V69eqgQCAp5BNta++n5uqqmv2Esk6o2znS31czbQ67+LP2OTX1qVn/lVde0ZAhQ4Lm1fTTwIEDtW3bNn366af64osvdM011+jCCy8MOoQH4MgIM0AL4nQ6g07aDZcBAwbI6/UqPz9fv/nNb8K+/WjQvn175eXlBaY3b96ssrKywLTT6ZSkoP5NT09Xp06d9PPPP2v8+PENbjslJUXjxo3TuHHjdNVVV2n06NE6cOCAUlNTI/CXACcewgzQgnTr1k0rV67U9u3b1apVq7D9WJ566qkaP368brzxRj311FMaMGCA9u3bp0WLFun000/XJZdcEpb9NKfzzz9fzz//vM4++2z5fD7NmDEjaDQnLS1NCQkJ+uyzz9S5c2fFx8fL5XLpz3/+s+666y6lpKTo4osvltvt1rfffquDBw9q2rRp+tvf/qaOHTuqf//+iomJ0bvvvqsOHTqodevWzffHAhbD1UxACzJ9+nTZ7Xb16dNH7du3V25ubti2PXfuXN14442655571LNnT11xxRVauXJlg1dMWc1TTz2lzMxMDRs2TNddd52mT5+uxMTEwPzY2FjNmTNHL730kjIyMjRmzBhJ0q233qpXX31V8+bN0+mnn67hw4dr3rx5gUu5W7Vqpccff1yDBw/WmWeeqe3bt+uTTz5RTAz/ewYaizsAAwAASyP6AwAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAAS/t/32Et/r7hyQYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(mh[10,1:1439])\n",
    "plt.xlim([0,1440])\n",
    "plt.xticks(numpy.arange(0,1440,240), numpy.arange(0, 1440, 240))\n",
    "plt.ylim([0,10])\n",
    "plt.xlabel(\"time in minutes\")\n",
    "plt.ylabel('water level in m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0538122b-8d81-49d9-9a54-883a6c49ac7f",
   "metadata": {},
   "source": [
    "In the same manner, the matrix mh (or mth, when the full Theis function ftheis is used, just or testing), can be used to plot the water level at different distances. For this purpose we use a range of the first index in mh[0:100,1439]. The first index i with 0:100 specifies the distances r, the second index j specifies the fixed time: As we increment with steps of 60 seconds when we define the variable ts = arange(60,86400,1439), each increment of the index corresponds to 60 seconds or one minute. At j=1439 we are at 1 day after the start of pumping (please not that the index starts with 0, hence index 1439 corresponds to 1440 minutes = 1 day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 439,
   "id": "d6745275-3250-42a1-bf43-c85c5c0f2c05",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAG2CAYAAACKxwc0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4o0lEQVR4nO3deXxU9b3/8fdkm+whCWSDsIPIvqogslysXKVWfvZapS4grWKLCqIVqXpBq8SlKvbaakWLC263FbzUDaJCAFFklU3ZISyJgYRksk6Syff3xyQDYc3AhJkTXs/HYx5zzvecM/l8iTBvv9+z2IwxRgAAABYV5O8CAAAAzgVhBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWJpfw8zSpUt17bXXKi0tTTabTR999FG97cYYzZgxQ2lpaYqIiNCwYcO0efNm/xQLAAACkl/DTGlpqXr16qWXXnrppNufeeYZPf/883rppZe0atUqpaSk6Gc/+5mKi4vPc6UAACBQ2QLlQZM2m03z58/X6NGjJblHZdLS0jR58mRNnTpVkuR0OpWcnKynn35aEyZM8GO1AAAgUIT4u4BT2b17t3Jzc3XVVVd52ux2u4YOHaoVK1acMsw4nU45nU7Pek1NjQoKCpSYmCibzdbodQMAgHNnjFFxcbHS0tIUFHT6iaSADTO5ubmSpOTk5HrtycnJ2rt37ymPy8jI0GOPPdaotQEAgPNj3759atWq1Wn3CdgwU+f40RRjzGlHWKZNm6YpU6Z41ouKitS6dWvt27dPsbGxjVYnAADwHYfDofT0dMXExJxx34ANMykpKZLcIzSpqame9ry8vBNGa45lt9tlt9tPaI+NjSXMAABgMQ05RSRg7zPTrl07paSkKDMz09NWWVmprKwsDRo0yI+VAQCAQOLXkZmSkhLt2LHDs757926tX79eCQkJat26tSZPnqyZM2eqU6dO6tSpk2bOnKnIyEj9+te/9mPVAAAgkPg1zKxevVrDhw/3rNed6zJ27Fi98cYbevDBB1VeXq7f//73OnLkiC699FItWrSoQfNnAADgwhAw95lpLA6HQ3FxcSoqKuKcGQAALMKb7++APWcGAACgIQgzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0ggzAADA0kL8XcDpVFdXa8aMGXrnnXeUm5ur1NRUjRs3To888oiCgshhAACcD64ao/Iql8qc1SqrdKms0qXyqmqVOl2169X13ssrXSfsV17pUlnV0e0f/m6Q0ppF+KS+gA4zTz/9tF555RW9+eab6tatm1avXq3bb79dcXFxmjRpkr/LAwAgoBhjVFFVo9LKapU5Xe73yvqho7TSHUqOfS+va68LJZ5j3W0VVTU+r7XUWe2zzwroMPPNN9/ouuuu06hRoyRJbdu21XvvvafVq1f7uTIAAM6ds9qlUqdLpc5qldaGjlKnO4CUOF2eIFK3vS5klDqPho+jx7i3GdN49dpsUlRYiCLCghUZFqzIsJDa96PrEWHBigoLVsQx2yJCj+5bd2x6QqTP6groMDN48GC98sor2rZtmzp37qzvv/9ey5cv16xZs055jNPplNPp9Kw7HI7zUCkA4EJQ7apRqdOlkspqlVRUq8RZGyyc1So+ZrmkXkCpDRzHLJfUBpYqV+Mlj7pwEWU/GiSi7CGKDA1WpD1YUZ6wUX+fiLBgRdvrQsnRQBJlD5E9JEg2m63Raj5bAR1mpk6dqqKiInXp0kXBwcFyuVx68sknNWbMmFMek5GRoccee+w8VgkACGSuGqMS59HgUVwbQtxhpErFFXUBo0oltUGj1Hk0rHjWndVyVvt+ukWS7CFBiraHuMNGbZiItIcoujZkRNuPhpGosGBF2kPcQcMerMjQ2na7O5REhYUoIjRYQUGBFzoaS0CHmQ8++EBz587Vu+++q27dumn9+vWaPHmy0tLSNHbs2JMeM23aNE2ZMsWz7nA4lJ6efr5KBgD4SF0IKa6oqn13BwxH7Xpd2CiuqAsoR9uLj9leVunyeW1hnvDhDg8x4SFHA0Vt6Ii2H22Lrt3v2LYo+9GRkpBgLmo5FzZjGnN27dykp6froYce0sSJEz1tTzzxhObOnasff/yxQZ/hcDgUFxenoqIixcbGNlapAIBjVLtqVFwbPI59d7+qjns/GlDq2ksq3OeE+FJdAKn3Cj8aMGLC3aMd7rZgRdtDFWUP9ux3bGgJJXw0Om++vwN6ZKasrOyES7CDg4NVU9M4w3wAgKNXxDgqquQor3K/V1TXLrvf6wWV8ioV1+5TXFElR3m1yqt8F0TsIUGKCQ9RTHioomtDR13AiA13B47jt9WFkxh7qDuI2INlDwn2WU0ILAEdZq699lo9+eSTat26tbp166Z169bp+eef1/jx4/1dGgAEtCpXjRzlVSqqDSBFdcueNnfo8ISV8mMDS5XPTkyNDAv2BJFj32PrlutCR217XUCpCyUx4aEKC2EUBKcX0NNMxcXFevTRRzV//nzl5eUpLS1NY8aM0X//938rLCysQZ/BNBMAq6p21XiCSGFZpSeQFJVXqaisSoXHrh8TVIrKq3xynkhwkM0TOmIj3CEj9phQUtcWEx6i2Ii6kHLsPiGcC4Kz5s33d0CHGV8gzADwt8rqGhWWV6qwrKr2VekOImVVKiyv1JGyYwOKe7+i2qmccxVjdweN2IhQxdUFkohQxdWGj7ja5br22IgQz3pkWHBAXoaLC0OTOWcGAAJJTY1RcUW1CsoqdaSsUoVllTpSWlW7fDSY1LUXlbu3nesoSbT9aOiIiwhVs8ijy7HHr4cfXWZkBBcKwgyAC5IxRo6Kah0prVR+aaWOlFa6Q0qpO5C43+teR9drznIsO8imY8JImJpFhqpZ7XJcRKjiI0MVFxmqZhFhiqsNJ81qwwpXzgCnR5gB0CRUu2pUUFapgtJKFZS4A0rBsUHlmFd+qXtUpfosk0lUWLCaRYYpPipU8ZFh7uXI0GPe65Zr1yPCFBMeckHdxAw4nwgzAAJSTY1RYXmV8kucOlxSqfxSp/JrQ0p+iXu5oLRSh0udKih1T/OcjaiwYMVHhSkhyh0+jr6HKj6qLpAcG1xCucQXCDCEGQDnTZWrRgWllTpU7NTh2pByuMSpw8VO5Ze6lw/VLheUVsrl5ciJzSZPIEmIClPiMe/xnmW74qNCPaElPJRgAlgdYQbAOTHGqLCsSnnF7iByqKTC/e5Zd9aGF3dA8VZcRKgSo8PUPMquhKgwNY8JU0KUXc2j3cEkISpMidHuwNIsMkzBTOUAFxzCDICTqnbVKL+0UnkOp/KKK/RT7XtesVN5DqcOFVd4woo3N1gLDrK5Q0m0O5C0iLareczRcNI8xq7E2u0JUWHcMA3AGRFmgAuMMUYFpZX6yeHUT44K/eSoUK6jNqw4KvRTcYXyHO5pIG9meeIjQ5UUE64WMXbPq3l0mHs5OlzNY9zBJT4yjBNhAfgUYQZoQqpdNcordiqnqEK5Re6QkltUrpyi+qGlsrphzzcLDrKpeXSYkmPDlRRjV4sY93tSrF1JMXVtdjWPtjOCAsBvCDOARbhqjA4VO3WwqFw5hRU6WFiug0Xlyi2q0MEid2g5VNzw0ZTEKHdISY61KyUuXEkx4Z715NhwJcXalRhl5xwUAAGPMAMEiLLKah04Uq79heXuoFJYroOFFTpQu5xbVNGg+6KEBtuUHBuulNhwpcSFKzXOHVJS4txtdUGFy4sBNBWEGeA8Ka6o0v4j5bWvMu0/Ul4bXsp04Ei5jjTgPinBQTalxLoDSmqzCKXVhpWUuAilNXMHluZRds5JAXBBIcwAPuKsdunAkXJlF5Rp35Fy7Sso074Cd2jZd6SsQTd1iwkPUctmEe5XfITSmrlfLZuFKzUuQsmx4Uz7AMBxCDOAF4rKq7Q3v1R788uUXVDmWd5XUKYcR4XO9Az6+MhQtYqPVKv4CLWKrwst7vWW8RGKDQ89Px0BgCaEMAMcp7CsUrsPl2pPfql2Hz4aWPbml55xKigiNFitEyKVnhCp9IQIpccfXW4VH6loO3/lAMDX+JcVF6SKKpf25Jdq16FS7T7sft91uES7D5eecTqoRYxdbRPdIaVNQpTa1C0nRioxKkw2G9NAAHA+EWbQpOWXOLUjr0Q7DpVo16FS7TxUop2HSrT/SPlpp4SSY+1qmxilds2j1CYxSu2aR6pNYpRaJ0QqitEVAAgo/KsMyzPGKK/YqW0/FWv7TyXanleiHXnF2pFXctppodjwELVvEa32zaPUvkWU2jWPVtvmkWqbGEVgAQAL4V9sWEphWaV+zC3W1txibf2pWNtyi7Xtp2I5KqpPur/NJrVsFqEOLaLVMSlaHVpEq0OLKHVIimZKCACaCMIMAlKVq0a7DpXqhxyHfsh16IecYv2Y41BesfOk+wfZpLaJUeqUHK1OSTHqlFwXXKIVEcbN4QCgKSPMwO9KnNX6IcehLQcd2nywSJsPOrT9pxJVuk7+/KBW8RHqkhKjzskxuiglRp2SYtS+RZTCQwktAHAhIszgvHJUVGnzAYc2HSjSxgNF2nSgSLvzS096Mm60PURdUmJ0cWqsuqTGqEtKrDonRyuGe7EAAI5BmEGjqahyafNBh77fV6gN+wu1YX+Rdh0uPem+KbHh6pYWq25pseqaFquuqXFqFR/BbfkBAGdEmIFPGGO0N79M6/Yd0brsQq3NPqIfc4pP+mDEls0i1KNlnHq0ilP3lnHqnharxGi7H6oGADQFhBmclYoqlzYeKNLqPUe0Zm+B1mYXqqC08oT9mkeHqVerZuqV3kw9W8WpZ6tmSogK80PFAICmijCDBimuqNLqvUe0cleBVu0p0Mb9RSecoBsWHKTuLWPVp3W8+rRupt7pzdSyWQSXPwMAGhVhBidVXFGllbsK9M2ufH23u0CbDxbp+Bmj5tFh6t8mQf3bxqtvm3h1S4uVPYQrigAA5xdhBpLc00Zr9x7R1zsP6+sd+dp4oEiu49JLm8RIXdI2QZe0S9CAtglqkxjJqAsAwO8IMxcoY4x25JVo6fbDWrrtkFbuzldFVf1po7aJkRrYobkua5+gS9slKiUu3E/VAgBwaoSZC0hZZbW+3pGvr37M05Ktecopqqi3PSnGrsEdm2tgh0QN6thcLZtF+KlSAAAajjDTxB0sLNcXP/ykL3/I0ze78lVZfXT0JSwkSJe2S9CQTi00pHMLdU6OZtoIAGA5hJkmaEdeiRZuztXCzbnasL+o3rZW8REa0SVJw7sk6bL2iTwCAABgeYSZJmJHXrEWfJ+jTzfmaEdeiafdZpP6t4nXiIuTNaJLkjomMfoCAGhaCDMWtje/VB9vyNG/vz+oH3OLPe2hwTYN6tBc/9k9RVdenKwWMdxdFwDQdBFmLKaovEqfbMjRv9bs09rsQk97aLBNQzq10M97pWrExcmK5WGMAIALBGHGAlw1Rst3HNa/1uzXws25npN4g2zS5R2b69qeaRrZLUVxkQQYAMCFhzATwA4VO/W/q/fpve+ytf9Iuae9c3K0ftm3lUb3aankWO79AgC4sBFmAowxRt/uKtDclXu1aHOuqlzuu/DGRYRqdO80/Ve/dHVvGctJvAAA1CLMBIjK6hp9vOGgXlu2W1tyHJ72Pq2b6eZL2+jnPVO5jBoAgJMgzPhZUXmV3vsuW298vUe5DvcdeSNCg/X/+rbUzZe2Vre0OD9XCABAYCPM+ElhWaVmL9ulN77eo9JKlySpRYxd4wa11c2XtlazyDA/VwgAgDUQZs6zorIqvbZ8l+Z8vUclzmpJ7hN6f3tFe13XO032EKaSAADwBmHmPClxVmv20l36x/LdKq4NMV1SYjT5ys4a2S2ZE3oBADhLhJlGVlNj9K+1+/Xswq06VOyUJF2UHKPJV3bSyG4pCgoixAAAcC4IM43ou90Fevzjzdp0wH11UpvESP1h5EW6pnsqIQYAAB8hzDSCPEeFHvt4iz7ZkCNJirGH6J4RHTV2UFvOiQEAwMcIMz5kjNG8tQf0+MdbVFRepSCbdNMlrTXlZ53VPJqHPQIA0BgIMz6SU1SuP87bqMVbD0mSerSM09O/7KmuabF+rgwAgKaNMHOOjDH639X79MTHP6jYWa2w4CBN/lkn3XlFe4UEB/m7PAAAmjzCzDkoq6zWH+dt1EfrD0pyP3rg2f/qqY5JMX6uDACAC4fXYaaiokL/8z//o8WLFysvL081NTX1tq9du9ZnxQWy3YdLddfba7T1p2IFB9n0wFUX6c4h7RXMVUoAAJxXXoeZ8ePHKzMzU//1X/+lSy655IK82duizbm6/3+/V7GzWs2j7frrr/vo0vaJ/i4LAIALktdh5pNPPtGnn36qyy+/vDHqCWg1NUbPZW7VXxfvlCT1bxOvv97cV8mx4X6uDACAC5fXYaZly5aKibnwzglx1RhN/XCD/rVmvyTp9svb6o/XXKxQTvIFAMCvvP4mfu655zR16lTt3bu3MeoJSFWuGt37/jr9a81+BdmkP9/QS9Ov7UaQAQAgAHg9MtO/f39VVFSoffv2ioyMVGhoaL3tBQUFPisuEFRUuXT3u2v1xQ95Cg226S839dHVPVL9XRYAAKjldZgZM2aMDhw4oJkzZyo5uWk/7bmsslp3vrVGy3cclj0kSK/c2k/DL0ryd1kAAOAYXoeZFStW6JtvvlGvXr0ao56AUVHl0th/fKdVe44oMixYr48doIEduGIJAIBA43WY6dKli8rLyxujloBhjPtk31V7jigmPERvjr9EfVvH+7ssAABwEl6fwfrUU0/p/vvv15IlS5Sfny+Hw1Hv1RS8unSX/m/9QQUH2fT3W/sRZAAACGBej8z853/+pyRpxIgR9dqNMbLZbHK5XL6pzE+WbM3T05//KEmafm1XDerQ3M8VAQCA0/E6zCxevLgx6ggIuw6V6J731qnGSDcNSNetl7Xxd0kAAOAMvA4zQ4cObYw6/K64okp3vLVaxRXV6tcmXo9d161JX6kFAEBTEfB3fTtw4IBuueUWJSYmKjIyUr1799aaNWt8+jOMMbrvg/XaeahUKbHhevmWvrKHBPv0ZwAAgMbh9cjM+XTkyBFdfvnlGj58uD777DMlJSVp586datasmU9/zscbcvTFD3kKCwnSq7f1U1IMz1oCAMAqAjrMPP3000pPT9ecOXM8bW3btvXpzyivdCnj0x8kSXcP76ierZr59PMBAEDjCuhppgULFqh///664YYblJSUpD59+mj27NmnPcbpdHp1ufjfl+7UwaIKtWwWoTuHtPdl+QAA4DwI6DCza9cuvfzyy+rUqZMWLlyou+66S/fee6/eeuutUx6TkZGhuLg4zys9Pf2U+x4sLNcrWTslSdOu6aLwUM6TAQDAamzGGOPNAT/99JMeeOABffnll8rLy9Pxh/vyPjNhYWHq37+/VqxY4Wm79957tWrVKn3zzTcnPcbpdMrpdHrWHQ6H0tPTVVRUpNjY2Hr73vveOi34/qAuaZugDyZcxtVLAAAECIfDobi4uJN+fx/P63Nmxo0bp+zsbD366KNKTU1t1ACQmpqqrl271mu7+OKL9eGHH57yGLvdLrvdfsbPXr2nQAu+PyibTfrva7sSZAAAsCivw8zy5cu1bNky9e7duxHKqe/yyy/X1q1b67Vt27ZNbdqc283samqMHvv3FknSjf3T1b1l3Dl9HgAA8B+vz5lJT08/YWqpsdx333369ttvNXPmTO3YsUPvvvuuXn31VU2cOPGcPvdfa/dr44EixdhDdP9VF/moWgAA4A9eh5lZs2bpoYce0p49exqhnPoGDBig+fPn67333lP37t31pz/9SbNmzdLNN9981p/prHbp2YXu0Z57RnRUi5gzT0kBAIDA5fUJwPHx8SorK1N1dbUiIyMVGhpab3tBQYFPCzxXx59A9MWWn/Tbt1YrKcau5VP/Q2EhAX1BFwAAF6RGPQF41qxZZ1tXQPhkY44k6ZoeqQQZAACaAK/DzNixYxujjvOiosqlzC0/SZKu7ZXq52oAAIAvNCjMOBwOzxDPme6oe6ahIH9auu2QSpzVSo0LV5/0eH+XAwAAfKBBYSY+Pl45OTlKSkpSs2bNTnpPFmOMbDabT2+a52sfbzg6xRQUxH1lAABoChoUZr766islJCRIkhYvXtyoBTWWiiqXvvjBPcX0855MMQEA0FQ0KMwMHTr0pMtWsmz7IZVVutSyWYR6pzfzdzkAAMBHLpjLeT7flCtJGtWzcR/BAAAAzq8LJsxkbTskiSkmAACamgsmzFRU1Sg9IUI9eA4TAABNygUTZiRpVI80ppgAAGhiLqgwwxQTAABNT4OuZurTp0+DRzTWrl17TgU1ltYJEeqWFrg39AMAAGenQWFm9OjRjVxG4xvZLYUpJgAAmiCvn5ptNXVP3fxu6z4N6NzK3+UAAIAG8Oap2Wd1zkxhYaFee+01TZs2TQUFBZLc00sHDhw4m487Lzonx/i7BAAA0Ai8fmr2hg0bdOWVVyouLk579uzRHXfcoYSEBM2fP1979+7VW2+91Rh1njOmmAAAaJq8HpmZMmWKxo0bp+3btys8PNzTfvXVV2vp0qU+LQ4AAOBMvA4zq1at0oQJE05ob9mypXJzc31SFAAAQEN5HWbCw8PlcDhOaN+6datatGjhk6IAAAAayuswc9111+nxxx9XVVWVJPe5KNnZ2XrooYf0y1/+0ucFAgAAnI7XYebPf/6zDh06pKSkJJWXl2vo0KHq2LGjYmJi9OSTTzZGjQAAAKfk9dVMsbGxWr58ub766iutXbtWNTU16tu3r6688srGqA8AAOC0vL5p3p49e9S2bdtGKsf3vLnpDgAACAyNetO89u3ba/Dgwfr73//uuWEeAACAv3gdZlavXq2BAwfqiSeeUFpamq677jr985//lNPpbIz6AAAATsvrMNO3b189++yzys7O1meffaakpCRNmDBBSUlJGj9+fGPUCAAAcEo+edDk2rVr9Zvf/EYbNmyQy+XyRV0+wzkzAABYT6M/aFKS9u3bp2eeeUa9e/fWgAEDFBUVpZdeeulsPw4AAOCseH1p9quvvqp33nlHX3/9tS666CLdfPPN+uijjyx1hRMAAGg6vA4zf/rTn3TTTTfpxRdfVO/evRuhJAAAgIbzOsxkZ2fLZrM1Ri0AAABe8/qcGZvNpmXLlumWW27RwIEDdeDAAUnS22+/reXLl/u8QAAAgNPxOsx8+OGHGjlypCIiIrRu3TrP/WWKi4s1c+ZMnxcIAABwOl6HmSeeeEKvvPKKZs+erdDQUE/7oEGDtHbtWp8WBwAAcCZeh5mtW7dqyJAhJ7THxsaqsLDQFzUBAAA0mNdhJjU1VTt27Dihffny5Wrfvr1PigIAAGgor8PMhAkTNGnSJK1cuVI2m00HDx7UO++8owceeEC///3vG6NGAACAU/L60uwHH3xQRUVFGj58uCoqKjRkyBDZ7XY98MADuvvuuxujRgAAgFM662czlZWVacuWLaqpqVHXrl0VHR3t69p8gmczAQBgPd58f3s9MlMnMjJS/fv3P9vDAQAAfKJBYeb6669v8AfOmzfvrIsBAADwVoPCTFxcXGPXAQAAcFYaFGbmzJnT2HUAAACcFa8vzQYAAAgkhBkAAGBphBkAAGBphBkAAGBpXoWZqqoqDR8+XNu2bWusegAAALziVZgJDQ3Vpk2bZLPZGqseAAAAr3g9zXTbbbfp9ddfb4xaAAAAvOb14wwqKyv12muvKTMzU/3791dUVFS97c8//7zPigMAADgTr8PMpk2b1LdvX0k64dwZpp8AAMD55nWYWbx4cWPUAQAAcFbO+tLsHTt2aOHChSovL5ckGWN8VhQAAEBDeR1m8vPzNWLECHXu3FnXXHONcnJyJEm//e1vdf/99/u8QAAAgNPxOszcd999Cg0NVXZ2tiIjIz3tN954oz7//HOfFgcAAHAmXp8zs2jRIi1cuFCtWrWq196pUyft3bvXZ4UBAAA0hNcjM6WlpfVGZOocPnxYdrvdJ0UBAAA0lNdhZsiQIXrrrbc86zabTTU1NXr22Wc1fPhwnxYHAABwJl5PMz377LMaNmyYVq9ercrKSj344IPavHmzCgoK9PXXXzdGjQAAAKfk9chM165dtWHDBl1yySX62c9+ptLSUl1//fVat26dOnTo0Bg1AgAAnJLNeHmDmOzsbKWnp5/0br/Z2dlq3bq1z4rzBYfDobi4OBUVFSk2Ntbf5QAAgAbw5vvb65GZdu3a6dChQye05+fnq127dt5+HAAAwDnxOswYY046KlNSUqLw8HCfFHUqGRkZstlsmjx5cqP+HAAAYB0NPgF4ypQpktxXLz366KP1Ls92uVxauXKlevfu7fMC66xatUqvvvqqevbs2Wg/AwAAWE+Dw8y6deskuUdmNm7cqLCwMM+2sLAw9erVSw888IDvK5R71Ofmm2/W7Nmz9cQTTzTKzwAAANbU4DBT97Ts22+/XS+++OJ5PZl24sSJGjVqlK688sozhhmn0ymn0+lZdzgcjV0eAADwI6/vMzNnzpzGqOOU3n//fa1du1arVq1q0P4ZGRl67LHHGrkqAAAQKLwOM5L7/JV//vOfys7OVmVlZb1t8+bN80lhkrRv3z5NmjRJixYtavDJxdOmTfOc3yO5R2bS09N9VhMAAAgsXl/N9P777+vyyy/Xli1bNH/+fFVVVWnLli366quvFBcX59Pi1qxZo7y8PPXr108hISEKCQlRVlaW/vKXvygkJEQul+uEY+x2u2JjY+u9AABA0+X1yMzMmTP1wgsvaOLEiYqJidGLL76odu3aacKECUpNTfVpcSNGjNDGjRvrtd1+++3q0qWLpk6dquDgYJ/+PAAAYD1eh5mdO3dq1KhRktyjIKWlpbLZbLrvvvv0H//xHz49XyUmJkbdu3ev1xYVFaXExMQT2gEAwIXJ62mmhIQEFRcXS5JatmypTZs2SZIKCwtVVlbm2+oAAADOwOuRmSuuuEKZmZnq0aOHfvWrX2nSpEn66quvlJmZqREjRjRGjfUsWbKk0X8GAACwDq/DzEsvvaSKigpJ7iuHQkNDtXz5cl1//fV69NFHfV4gAADA6Xj91Gyr4anZAABYT6M+NbvusQLbtm076wIBAAB8xeswEx0dreeee05dunRRWlqaxowZo1deeUU//vhjY9QHAABwWmc9zZSbm6slS5ZoyZIlysrK0rZt25SUlKScnBxf13hOmGYCAMB6GnWaqU5MTIzi4+MVHx+vZs2aKSQkRCkpKWf7cQAAAGfF6zAzdepUXXbZZWrevLkeeeQRVVZWatq0afrpp5+0bt26xqgRAADglLyeZgoKClKLFi1033336brrrtPFF1/cWLX5BNNMAABYjzff317fZ2bdunXKysrSkiVL9Nxzzyk4OFhDhw7VsGHDNGzYsIAPNwAAoGk55/vMfP/995o1a5bmzp2rmpqakz7J2p8YmQEAwHoadWRGco/O1F3JtGzZMjkcDvXu3VvDhw8/q4IBAADOltdhJj4+XiUlJerVq5eGDRumO+64Q0OGDGHUAwAA+IXXYebtt98mvAAAgIDhdZj5+c9/3hh1AAAAnJWzvmkeAABAICDMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASyPMAAAASwvoMJORkaEBAwYoJiZGSUlJGj16tLZu3ervsgAAQAAJ6DCTlZWliRMn6ttvv1VmZqaqq6t11VVXqbS01N+lAQCAAGEzxhh/F9FQhw4dUlJSkrKysjRkyJAGHeNwOBQXF6eioiLFxsY2coUAAMAXvPn+DjlPNflEUVGRJCkhIeGU+zidTjmdTs+6w+Fo9LoAAID/BPQ007GMMZoyZYoGDx6s7t27n3K/jIwMxcXFeV7p6ennsUoAAHC+WWaaaeLEifrkk0+0fPlytWrV6pT7nWxkJj09nWkmAAAspMlNM91zzz1asGCBli5detogI0l2u112u/08VQYAAPwtoMOMMUb33HOP5s+fryVLlqhdu3b+LgkAAASYgA4zEydO1Lvvvqv/+7//U0xMjHJzcyVJcXFxioiI8HN1AAAgEAT0OTM2m+2k7XPmzNG4ceMa9Blcmg0AgPU0mXNmAjhnAQCAAGGZS7MBAABOhjADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAsjTADAAAszRJh5m9/+5vatWun8PBw9evXT8uWLfN3SQAAIEAEfJj54IMPNHnyZD388MNat26drrjiCl199dXKzs72d2kAACAA2Iwxxt9FnM6ll16qvn376uWXX/a0XXzxxRo9erQyMjLOeLzD4VBcXJyKiooUGxvbmKUCAAAf8eb7O+Q81XRWKisrtWbNGj300EP12q+66iqtWLHipMc4nU45nU7PelFRkST3HwoAALCGuu/thoy5BHSYOXz4sFwul5KTk+u1JycnKzc396THZGRk6LHHHjuhPT09vVFqBAAAjae4uFhxcXGn3Segw0wdm81Wb90Yc0JbnWnTpmnKlCme9cLCQrVp00bZ2dln/MMIZA6HQ+np6dq3b59lp8uaQh+kptGPptAHiX4EkqbQB6lp9KMp9EFyf9cXFxcrLS3tjPsGdJhp3ry5goODTxiFycvLO2G0po7dbpfdbj+hPS4uztK/1DqxsbGW70dT6IPUNPrRFPog0Y9A0hT6IDWNfjSFPjR0ECKgr2YKCwtTv379lJmZWa89MzNTgwYN8lNVAAAgkAT0yIwkTZkyRbfeeqv69++vgQMH6tVXX1V2drbuuusuf5cGAAACQMCHmRtvvFH5+fl6/PHHlZOTo+7du+vTTz9VmzZtGnS83W7X9OnTTzr1ZCVNoR9NoQ9S0+hHU+iDRD8CSVPog9Q0+tEU+uCtgL/PDAAAwOkE9DkzAAAAZ0KYAQAAlkaYAQAAlkaYAQAAltakw8zf/vY3tWvXTuHh4erXr5+WLVvm75JOa+nSpbr22muVlpYmm82mjz76qN52Y4xmzJihtLQ0RUREaNiwYdq8ebN/ij2FjIwMDRgwQDExMUpKStLo0aO1devWevtYoR8vv/yyevbs6bnp1MCBA/XZZ595tluhD8fLyMiQzWbT5MmTPW1W6MeMGTNks9nqvVJSUjzbrdCHOgcOHNAtt9yixMRERUZGqnfv3lqzZo1ne6D3pW3btif8Lmw2myZOnCgp8OuvU11drUceeUTt2rVTRESE2rdvr8cff1w1NTWefazQl+LiYk2ePFlt2rRRRESEBg0apFWrVnm2W6EPPmOaqPfff9+Ehoaa2bNnmy1btphJkyaZqKgos3fvXn+Xdkqffvqpefjhh82HH35oJJn58+fX2/7UU0+ZmJgY8+GHH5qNGzeaG2+80aSmphqHw+Gfgk9i5MiRZs6cOWbTpk1m/fr1ZtSoUaZ169ampKTEs48V+rFgwQLzySefmK1bt5qtW7eaP/7xjyY0NNRs2rTJGGONPhzru+++M23btjU9e/Y0kyZN8rRboR/Tp0833bp1Mzk5OZ5XXl6eZ7sV+mCMMQUFBaZNmzZm3LhxZuXKlWb37t3miy++MDt27PDsE+h9ycvLq/d7yMzMNJLM4sWLjTGBX3+dJ554wiQmJpqPP/7Y7N692/zzn/800dHRZtasWZ59rNCXX/3qV6Zr164mKyvLbN++3UyfPt3Exsaa/fv3G2Os0QdfabJh5pJLLjF33XVXvbYuXbqYhx56yE8Veef4MFNTU2NSUlLMU0895WmrqKgwcXFx5pVXXvFDhQ2Tl5dnJJmsrCxjjHX7YYwx8fHx5rXXXrNcH4qLi02nTp1MZmamGTp0qCfMWKUf06dPN7169TrpNqv0wRhjpk6dagYPHnzK7VbqS51JkyaZDh06mJqaGkvVP2rUKDN+/Ph6bddff7255ZZbjDHW+F2UlZWZ4OBg8/HHH9dr79Wrl3n44Yct0QdfapLTTJWVlVqzZo2uuuqqeu1XXXWVVqxY4aeqzs3u3buVm5tbr092u11Dhw4N6D4VFRVJkhISEiRZsx8ul0vvv/++SktLNXDgQMv1YeLEiRo1apSuvPLKeu1W6sf27duVlpamdu3a6aabbtKuXbskWasPCxYsUP/+/XXDDTcoKSlJffr00ezZsz3brdQXyf3v7Ny5czV+/HjZbDZL1T948GB9+eWX2rZtmyTp+++/1/Lly3XNNddIssbvorq6Wi6XS+Hh4fXaIyIitHz5ckv0wZeaZJg5fPiwXC7XCQ+jTE5OPuGhlVZRV7eV+mSM0ZQpUzR48GB1795dkrX6sXHjRkVHR8tut+uuu+7S/Pnz1bVrV0v14f3339fatWuVkZFxwjar9OPSSy/VW2+9pYULF2r27NnKzc3VoEGDlJ+fb5k+SNKuXbv08ssvq1OnTlq4cKHuuusu3XvvvXrrrbckWef3Ueejjz5SYWGhxo0bJ8la9U+dOlVjxoxRly5dFBoaqj59+mjy5MkaM2aMJGv0JSYmRgMHDtSf/vQnHTx4UC6XS3PnztXKlSuVk5NjiT74UsA/zuBc2Gy2euvGmBParMZKfbr77ru1YcMGLV++/IRtVujHRRddpPXr16uwsFAffvihxo4dq6ysLM/2QO/Dvn37NGnSJC1atOiE/3s7VqD34+qrr/Ys9+jRQwMHDlSHDh305ptv6rLLLpMU+H2QpJqaGvXv318zZ86UJPXp00ebN2/Wyy+/rNtuu82znxX6Ikmvv/66rr76aqWlpdVrt0L9H3zwgebOnat3331X3bp10/r16zV58mSlpaVp7Nixnv0CvS9vv/22xo8fr5YtWyo4OFh9+/bVr3/9a61du9azT6D3wVea5MhM8+bNFRwcfEL6zMvLOyGlWkXd1RtW6dM999yjBQsWaPHixWrVqpWn3Ur9CAsLU8eOHdW/f39lZGSoV69eevHFFy3ThzVr1igvL0/9+vVTSEiIQkJClJWVpb/85S8KCQnx1Bro/TheVFSUevTooe3bt1vmdyFJqamp6tq1a722iy++WNnZ2ZKs9Xdj7969+uKLL/Tb3/7W02al+v/whz/ooYce0k033aQePXro1ltv1X333ecZwbRKXzp06KCsrCyVlJRo3759+u6771RVVaV27dpZpg++0iTDTFhYmPr166fMzMx67ZmZmRo0aJCfqjo3df9xHtunyspKZWVlBVSfjDG6++67NW/ePH311Vdq165dve1W6cfJGGPkdDot04cRI0Zo48aNWr9+vefVv39/3XzzzVq/fr3at29viX4cz+l06ocfflBqaqplfheSdPnll59wm4Jt27Z5Hpprpb7MmTNHSUlJGjVqlKfNSvWXlZUpKKj+119wcLDn0mwr9UVyB/zU1FQdOXJECxcu1HXXXWe5Ppwzv5x2fB7UXZr9+uuvmy1btpjJkyebqKgos2fPHn+XdkrFxcVm3bp1Zt26dUaSef755826des8l5M/9dRTJi4uzsybN89s3LjRjBkzJuAus/vd735n4uLizJIlS+pdwllWVubZxwr9mDZtmlm6dKnZvXu32bBhg/njH/9ogoKCzKJFi4wx1ujDyRx7NZMx1ujH/fffb5YsWWJ27dplvv32W/Pzn//cxMTEeP4uW6EPxrgvjw8JCTFPPvmk2b59u3nnnXdMZGSkmTt3rmcfK/TF5XKZ1q1bm6lTp56wzQr1G2PM2LFjTcuWLT2XZs+bN880b97cPPjgg559rNCXzz//3Hz22Wdm165dZtGiRaZXr17mkksuMZWVlcYYa/TBV5psmDHGmL/+9a+mTZs2JiwszPTt29dzeXCgWrx4sZF0wmvs2LHGGPflgtOnTzcpKSnGbrebIUOGmI0bN/q36OOcrH5JZs6cOZ59rNCP8ePHe/7badGihRkxYoQnyBhjjT6czPFhxgr9qLs3RmhoqElLSzPXX3+92bx5s2e7FfpQ59///rfp3r27sdvtpkuXLubVV1+tt90KfVm4cKGRZLZu3XrCNivUb4wxDofDTJo0ybRu3dqEh4eb9u3bm4cfftg4nU7PPlboywcffGDat29vwsLCTEpKipk4caIpLCz0bLdCH3zFZowxfhkSAgAA8IEmec4MAAC4cBBmAACApRFmAACApRFmAACApRFmAACApRFmAACApRFmAACApRFmAJzSsGHDNHnyZM9627ZtNWvWLL/V01hmzJih3r17+7sMAGeJMAOgwVatWqU777yzQftaKfg88MAD+vLLL/1dBoCzFOLvAgBYR4sWLfxdQqOIjo5WdHS0v8sAcJYYmQEgSSotLdVtt92m6Ohopaam6rnnnjthn+NHW2bMmKHWrVvLbrcrLS1N9957ryT39NTevXt13333yWazyWazSZLy8/M1ZswYtWrVSpGRkerRo4fee++9ej9j2LBhuvfee/Xggw8qISFBKSkpmjFjRr19CgsLdeeddyo5OVnh4eHq3r27Pv74Y8/2FStWaMiQIYqIiFB6erruvfdelZaWnrLvx08zjRs3TqNHj9af//xnpaamKjExURMnTlRVVdUZP+Mf//iHWrdurejoaP3ud7+Ty+XSM888o5SUFCUlJenJJ5885WcAODuEGQCSpD/84Q9avHix5s+fr0WLFmnJkiVas2bNKff/17/+pRdeeEF///vftX37dn300Ufq0aOHJGnevHlq1aqVHn/8ceXk5CgnJ0eSVFFRoX79+unjjz/Wpk2bdOedd+rWW2/VypUr6332m2++qaioKK1cuVLPPPOMHn/8cWVmZkqSampqdPXVV2vFihWaO3eutmzZoqeeekrBwcGSpI0bN2rkyJG6/vrrtWHDBn3wwQdavny57r77bq/+PBYvXqydO3dq8eLFevPNN/XGG2/ojTfeOO0xO3fu1GeffabPP/9c7733nv7xj39o1KhR2r9/v7KysvT000/rkUce0bfffutVLQDOwN9PugTgf8XFxSYsLMy8//77nrb8/HwTERFR7ynbbdq0MS+88IIxxpjnnnvOdO7c2VRWVp70M4/d93SuueYac//993vWhw4dagYPHlxvnwEDBpipU6caY9xPbQ4KCjrpU5uNMebWW281d955Z722ZcuWmaCgIFNeXn7SY6ZPn2569erlWR87dqxp06aNqa6u9rTdcMMN5sYbbzxlP6ZPn24iIyONw+HwtI0cOdK0bdvWuFwuT9tFF11kMjIyTvk5ALzHyAwA7dy5U5WVlRo4cKCnLSEhQRdddNEpj7nhhhtUXl6u9u3b64477tD8+fNVXV192p/jcrn05JNPqmfPnkpMTFR0dLQWLVqk7Ozsevv17Nmz3npqaqry8vIkSevXr1erVq3UuXPnk/6MNWvW6I033vCcBxMdHa2RI0eqpqZGu3fvPm19x+rWrZtntOf4Gk6lbdu2iomJ8awnJyera9euCgoKqtd2ps8B4B1OAAYgY4zXx6Snp2vr1q3KzMzUF198od///vd69tlnlZWVpdDQ0JMe89xzz+mFF17QrFmz1KNHD0VFRWny5MmqrKyst9/xx9tsNtXU1EiSIiIiTltXTU2NJkyY4Dl/51itW7ducP9OV4M3x5zN5wDwDmEGgDp27KjQ0FB9++23ni/8I0eOaNu2bRo6dOgpj4uIiNAvfvEL/eIXv9DEiRPVpUsXbdy4UX379lVYWJhcLle9/ZctW6brrrtOt9xyiyR38Ni+fbsuvvjiBtfas2dP7d+/X9u2bTvp6Ezfvn21efNmdezYscGfCcDamGYCoOjoaP3mN7/RH/7wB3355ZfatGmTxo0bV2965HhvvPGGXn/9dW3atEm7du3S22+/rYiICLVp00aSe8pl6dKlOnDggA4fPizJHZoyMzO1YsUK/fDDD5owYYJyc3O9qnXo0KEaMmSIfvnLXyozM1O7d+/2nHQrSVOnTtU333yjiRMnav369dq+fbsWLFige+655yz/dAAEOsIMAEnSs88+qyFDhugXv/iFrrzySg0ePFj9+vU75f7NmjXT7Nmzdfnll6tnz5768ssv9e9//1uJiYmSpMcff1x79uxRhw4dPPenefTRR9W3b1+NHDlSw4YNU0pKikaPHu11rR9++KEGDBigMWPGqGvXrnrwwQc9o0A9e/ZUVlaWtm/friuuuEJ9+vTRo48+qtTUVO//UABYgs2czWQ5AABAgGBkBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWBphBgAAWNr/Bypizxanr9gwAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(mh[0:100,1439])\n",
    "plt.xlim([0,100])\n",
    "plt.xticks(numpy.arange(0,100,10), numpy.arange(0, 100, 10))\n",
    "plt.ylim([0,10])\n",
    "plt.xlabel(\"distance in m\")\n",
    "plt.ylabel('water level in m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90ce85a8-598d-41f8-b460-2e2ffd7bb01b",
   "metadata": {},
   "source": [
    "Finally, we try the new function **ftheis()**. First, we call it with only the non-default value for the pumping rate. All other values are set as default. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "id": "03b7ef7f-6d9e-4f6f-a920-f3322d6fe615",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-15.800087838405927"
      ]
     },
     "execution_count": 384,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ftheis(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3e11324-520c-407a-a840-36d9b2ce7f21",
   "metadata": {},
   "source": [
    "If you want to specify just some (and keep the default values for the others), you must use the variable names and assign values to them. If you do that, you could even mix the variables of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "id": "65cca84d-63ea-4221-9f4e-173ea9eb4566",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-18.571619028767373"
      ]
     },
     "execution_count": 385,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ftheis(0.01, kf=0.001, Ha=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73790390-164b-4064-beec-217be93da6be",
   "metadata": {},
   "source": [
    "Finally, you can try to specify all other variables. Or you use all variables in the right order without their names, or you need to assign all variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "id": "02965b5e-e6f1-4c28-99ad-edb748b68ba7",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.3786358870470363"
      ]
     },
     "execution_count": 386,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ftheis(0.01,h0=6,Ha=50,kf=1E-5,rb=50,Ss=1E-4,ts=1440)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}