{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# 3. Slug Test - Falling Head"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = [5, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction and Conceptual Model\n",
    "\n",
    "A well partially penetrates a sandy unconfined aquifer that has a saturated depth of 32.57 ft. The top of the screen is located 0.47 ft below the water table and has 13.8 ft in length. The well and casing radii are 5 and 2 inches, respectively. The slug displacement is 1.48 ft. Head change has been recorded at the slug well. This slug test, taken from the AQTESOLV examples (Duffield, 2007), was reported in Batu (1998). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figs/Falling_head.png\" style=\"width:400pt\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt(\"data/falling_head.txt\", skiprows=2)\n",
    "to = data[:, 0] / 60 / 60 / 24  # convert time from seconds to days\n",
    "ho = (10 - data[:, 1]) * 0.3048  # convert drawdown from ft to meters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters and model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "rw = 5 * 0.0254  # well radius, m\n",
    "rc = 2 * 0.0254  # well casing radius, m\n",
    "L = 13.8 * 0.3048  # screen length, m\n",
    "b = 32.57 * 0.3048  # aquifer thickness, m\n",
    "zt = 0  # top of aquifer, m\n",
    "zst = -0.47 * 0.3048  # top of screen, m\n",
    "zsb = zst - L  # bottom of screen, m\n",
    "zb = zt - b  # bottom of aquifer, m\n",
    "H0 = 1.48 * 0.3048  # initial displacement in well, m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "convert measured displacement into volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "slug: 0.00366 m^3\n"
     ]
    }
   ],
   "source": [
    "Q = np.pi * rc**2 * H0\n",
    "print(f\"slug: {Q:.5f} m^3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create a multi-layer model. For this, we divide the second and third layers into 0.5 m thick layers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# z = np.array([zt, zt - 0.1, zst, zsb, zb])\n",
    "z = np.hstack(\n",
    "    (\n",
    "        zt,\n",
    "        np.linspace(zt - 0.01, zst, 5)[:-1],\n",
    "        np.linspace(zst, zsb, 5)[:-1],\n",
    "        np.linspace(zsb, zb, 10),\n",
    "    )\n",
    ")\n",
    "nlay = len(z) - 1\n",
    "nlay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.      , -0.01    , -0.043314, -0.076628, -0.109942, -0.143256,\n",
       "       -1.194816, -2.246376, -3.297936, -4.349496, -4.969256, -5.589016,\n",
       "       -6.208776, -6.828536, -7.448296, -8.068056, -8.687816, -9.307576,\n",
       "       -9.927336])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self.neq  5\n",
      "solution complete\n"
     ]
    }
   ],
   "source": [
    "ml = tft.Model3D(\n",
    "    kaq=10,\n",
    "    z=z,\n",
    "    Saq=[0.1] + (nlay - 1) * [1e-4],\n",
    "    kzoverkh=1,\n",
    "    tmin=1e-5,\n",
    "    tmax=0.01,\n",
    "    phreatictop=True,\n",
    ")\n",
    "w = tft.Well(\n",
    "    ml,\n",
    "    xw=0,\n",
    "    yw=0,\n",
    "    rw=rw,\n",
    "    tsandQ=[(0, -Q)],\n",
    "    layers=range(6, 6 + 5),\n",
    "    rc=rc,\n",
    "    wbstype=\"slug\",\n",
    ")\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x15b7b8230>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAESCAYAAACiv++BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABEQklEQVR4nO3deVxVdf7H8de5OzuiCaKoWBoqbmkgppkjhkultqlDZWbZYts4M6WV2q+psX0aG2ccbbEmS7NscRfBXURA3PdExQVckAso673f3x/oNVJTEDkX+Dwfj/sgzvle7vucxLf33O85R1NKKYQQQghRYQa9AwghhBA1lZSoEEIIUUlSokIIIUQlSYkKIYQQlSQlKoQQQlSSlKgQQghRSVKiQgghRCWZ9A7gTpxOJ0ePHsXHxwdN0/SOI4QQQgdKKfLy8ggODsZg+P33mlKiv3L06FFCQkL0jiGEEMINZGRk0KRJk98dIyX6Kz4+PkDZjvP19dU5jRBCCD3k5uYSEhLi6oTfIyX6K+cP4fr6+kqJCiFEHXc1H+vJxCIhhBCikqREhRBCiEqSEhVCCCEqSUpUCCGEqCQpUSGEEKKSpESFEEKISpJTXIQQABw4eYbDpwtc31tMBjo19cdslH9rC3E5UqJC1HF7s/L4aNleFmw9dtG6kAAPnv9DSwZ3aoxJylSIi2hKKaV3CHeRm5uLn58fdrtdLrYgar3juYX8feFOftp8lPN/C7Rs6I3RUHaCeWZuITlnSwBoXt+Tsf3C6BveSK+4QlSbinSBvBMVog4qdTh5/MsUthy2A9CvTUNeaZVBCBmuMcWYmJ3Xjn+sz+XAqbM89dVGvn48km43NdArthBuR0pUiDpo2ur9bDlsx99m4OdeJ2i67Q1YvKvcGAvwsMnG0E7D+dvpPny5vZiX525h8Qu342WVvzqEADmcW44czhV1wd6sPAZMXkNb526+DJiBz5n0shU2P2jRCwzGsu+z0+HoRgCU0cpM+vPamfsZHtWc/xsYrlN6Ia4/OZwrhLikUoeTv3y3BU+Hnc+9/onPmWzwqAdRoyFiVFmRnqcU7F8BK99BO5TIQ/zALqM/XyT2oV+7RnRtUV+37RDCXUiJClGHfLomnc0ZOfzb9gX+jmxo0AoejwebL8fyj3Es6xfXWIvRQlhoD0wt7oB1H0PceCZYvmF1YTte/r7ssK6HxajfxgjhBqREhagjDp8+ywdxexhgWE9/1oFmhMFT2X4mg0/Xf8qyg8tQlP90p4l3E0aEj2BgxBNY9y7FcmA1k23TGHzqNf4Zv5ex/cJ02hoh3EOlTvyaMmUKzZs3x2azERkZyYYNG353/Jw5cwgLC8Nms9GuXTsWLlxYbr1SigkTJtCoUSM8PDyIjo5m79695cZkZ2cTGxuLr68v/v7+jBw5kvz8/Eu+3r59+/Dx8cHf378ymydErfTZmgP4lWbztvVzANK7PsGo7f9l6PyhxB2MQ6Fo6tOU5r7Nae7bHB+zD4fzD/O39X8jZm4/vmnfDyw+dFC7eNy4gK/WHyS3sETnrRJCXxUu0dmzZzNmzBgmTpzIxo0b6dChAzExMRw/fvyS49etW8ewYcMYOXIkaWlpDBo0iEGDBrFt2zbXmHfffZfJkyczdepUkpKS8PLyIiYmhsLCQteY2NhYtm/fTlxcHPPnz2fVqlWMGjXqotcrKSlh2LBh9OjRo6KbJkStlVtYwuzkg0wyT8dH5XEmqB3P5G8m8VgiRs3IXS3uYu49c1lw7wLmDZ7HvMHzWPbAMsZGjKWRVyNOFZ7i71unsqDrQwD8xTyH4OJ0vk3OuMIrC1HLqQqKiIhQo0ePdn3vcDhUcHCwmjRp0iXHP/jgg2rAgAHllkVGRqonn3xSKaWU0+lUQUFB6r333nOtz8nJUVarVX3zzTdKKaV27NihAJWcnOwas2jRIqVpmjpy5Ei5n/3SSy+phx56SH3++efKz8+vQttmt9sVoOx2e4WeJ4S7++/KfWr4uDeVmuirnG80UK8sHa3CZ4SrPnP6qIzcjN99brGjWL2f/L4KnxGuus7sqg7/b6BSE33VytduU90mxauSUkf1bIQQ1aQiXVChd6LFxcWkpqYSHR3tWmYwGIiOjiYxMfGSz0lMTCw3HiAmJsY1Pj09nczMzHJj/Pz8iIyMdI1JTEzE39+fLl26uMZER0djMBhISkpyLUtISGDOnDlMmTLlqranqKiI3Nzccg8hapsSh5MZaw/wpHE+AEvaDeDnoysxaAYm9ZhEE58mv/t8s8HMC7e8QMcbOpJfks84PwslmpHbjVvxs+9k0bbM6tgMIdxShUr05MmTOBwOAgMDyy0PDAwkM/PSv0iZmZm/O/781yuNadiwYbn1JpOJgIAA15hTp07x6KOPMmPGjKs+x3PSpEn4+fm5HiEhIVf1PCFqkoVbj9EgdztRxh0cMVt540zZRRUeb/c4nQM7X9XPMBlMTOoxCW+zN2nZO/mkZQQAo0zz+WT1fpScbi7qqFpzReknnniCP/7xj9x+++1X/Zxx48Zht9tdj4wM+XxH1C5KKT5Znc4o03wcwCvNWpJXeob2DdrzVIenKvSzmvg04dWurwLw39JjbLJauMuwnpOH95Fy8PR1SC+E+6tQiTZo0ACj0UhWVla55VlZWQQFBV3yOUFBQb87/vzXK4357cSl0tJSsrOzXWMSEhJ4//33MZlMmEwmRo4cid1ux2Qy8dlnn10ym9VqxdfXt9xDiNpkQ3o2uUd308+wgSVenmx05OJp8uTtHm9jNpgr/PPuanEX/UP741BO3gluhlFzMtK0iE9W778O6YVwfxUqUYvFQufOnYmPj3ctczqdxMfHExUVdcnnREVFlRsPEBcX5xofGhpKUFBQuTG5ubkkJSW5xkRFRZGTk0NqaqprTEJCAk6nk8jISKDsc9NNmza5Hm+88QY+Pj5s2rSJwYMHV2Qzhag1PlubzuPGhWiaYnrDxgCMCB9BiG/lP7p46daXsBltbKOIRJuNocblbNixj0OnzlZVbCFqjAofzh0zZgzTp0/niy++YOfOnTz99NOcOXOGESNGAPDII48wbtw41/gXXniBxYsX88EHH7Br1y5ef/11UlJSePbZZwHQNI0XX3yRN998k59//pmtW7fyyCOPEBwczKBBgwBo3bo1ffv25YknnmDDhg2sXbuWZ599lqFDhxIcHOwaEx4e7no0btwYg8FAeHg49erVu9b9JESNcyq/iLSde3nAuJIVnh7sowgvsxfDwoZd08+t71Gf+1rdB8D0ho3w1IqINSxjbtrhqogtRI1S4RIdMmQI77//PhMmTKBjx45s2rSJxYsXuyYGHTp0iGPHLtzct1u3bnz99ddMmzaNDh068N133/Hjjz8SHn7hAtYvvfQSzz33HKNGjeLWW28lPz+fxYsXY7PZXGNmzpxJWFgYvXv3pn///nTv3p1p06Zdy7YLUavN33KMWMMSrFoJ028ouw/o0JuH4mf1u8Izr+zRto9iMphIMTpIs1p41LSEBRtlgpGoe+QuLr8id3ERtcngj1cy7eTD7PUo4slGDbEZbSy+bzH1ParmwvGvr3ud7/d+T/dixX+OZPBC8TM88uRf6dwsoEp+vhB6qUgX1JrZuUKIC/Ydz8f/2Gpu0OxMDygrzfta3VdlBQrwWPhjGDQDaywaOy1m7jeuYu7GI1X284WoCaREhaiFfkg7zP3GlWyyWkixGjEZTDza9tEqfY2mvk3p27wvANP9/bjNsJ3kzVsoKnVU6esI4c6kRIWoZZxORcLGXUQbNvKpX9mhqIE3DiTI69KnoV2Lx9s9DsAyT08OmY30KVlBws5LX0dbiNpISlSIWiYpPZsu+cs5bVKs8vQEYHjb4dfltVrWa0mPxj1QGsz18eZ+40rmbpRZuqLukBIVopYpO5S7ip98vHBq0CWwC6F+odft9c6f7vKTjzdNDFnY96wh+0zxdXs9IdyJlKgQtUhhiYM9WzfQzrCfud7eANzb8t7r+pq3N7md+rb6ZBsNrPL0YLC2kvlbjl7X1xTCXUiJClGLLNuZRX/HcpJsVo6aTfiYfejTrM91fU2zwczAmwYC8L2PN3cZ17Nw4y/X9TWFcBdSokLUIou3ZDDYuIa5PmXvQge0GIDNZLvCs67d+Xe7az08OGMqJujIMo7mFFz31xVCb1KiQtQSZ4tLKdmdgMmYR7xX2YSi859XXm/NfJvRJbALTg1+9PHiPuMquc+oqBOkRIWoJVbuPsGdrGG+txclmkab+m0ICwirttc//270B29vuhq2s3bzzmp7bSH0IiUqRC2xZPNBog2pzPXxAuDem67vhKLf6tOsDz5mH46aTSR7WGl0NI5Me2G1ZhCiukmJClELFJY4KN0Tx0Gbg30WCzajjX4t+lVrBpvJxoAWAwCY6+PFXYb1LN527ArPEqJmkxIVohZYsfsE0Wod873K3oX2btYbX0v130RhcMuye/eu8PSgrXEXazftqPYMQlQnKVEhaoFlWw5wh2EjS7zLJhQNCB2gS47WAa1p7tucIoOBld42go8u4XiuHNIVtZeUqBA1XGGJg9LdS9nhCdlGI/Ws9ega3FWXLJqm0b9FfwAWeHnR35jEku0yS1fUXlKiQtRwq/eepLdzLQvOndZyZ/M7MRvMuuXpH1pWous9bNxo3Mu6TVt1yyLE9SYlKkQNt2xzOrcZ01znhp6f3KOXZr7NCK8fjkPTWObtQaMjSziRV6RrJiGuFylRIWqw4lInpbsXk+Jp4IzBQCOvRnS4oYPesegXWjYzeKG3JwMM61m6Qw7pitpJSlSIGmz9/lP8wbGWRd5ls3L7hfbDoOn/a903tC8aGmk2G40s+0nZLId0Re2k/2+bEKLSErakc6tpM6s8PIALn0fqraFnQyKCIgBY5OVJg4zF5BaW6JxKiKonJSpEDeV0Kgp2LmGtl4lig8aNfjfSql4rvWO5uA7penlxp5bE8l3HdU4kRNWTEhWihkrLyCGqeB0Lz50b2r9FfzRN0znVBdHNojFpJvZYLfhZ0lm/WS68IGofKVEhaqj4bYfoZNrEBlvZrc76Na/ey/xdiZ/Vj+5NugOw1NsTj18WUVji0DmVEFVLSlSIGkgpxemtcSR5G3BqGm0C2hDiG6J3rIvENI8BYKmXJ71UIut+OalzIiGqlpSoEDXQ3uP5dMxfzdJfXWDBHd3R5A4sBjPpFjMNrPtYvWmP3pGEqFJSokLUQEu3HqaLOZVkmxVw3xL1tnhzW+OyQ7rx3h6o3QtxOJXOqYSoOlKiQtRAx7bEk+qtzh3KbU2Ij/sdyj3vfMEv8fKke+k6Ug+e1jmREFVHSlSIGuZITgE3Z69wHcqNCe2rc6Lfd/6Q7gGLmUbWXSzf8ovekYSoMlKiQtQwS7cdJcKccuFQbjP3PJR73q8P6SZ4WynYthCl5JCuqB2kRIWoYQ6krWCLdwlOTaNtQBua+DTRO9IV/XqW7q2Fq9l5LE/nREJUDSlRIWqQ02eKCTm+7MKs3NAYnRNdnTtCLhzSbWLbQfyWA3pHEqJKSIkKUYPE78wi0pRcYw7lnudl9qJ74x4ArPIycXrrYp0TCVE1pESFqEF2pK1hp3fBuUO5rWvEodzzzs/SXerlSZvcFWRkn9U5kRDXTkpUiBqioNhBg4wlvzqU696zcn/rjpA7sGgmDljMhNq2ELc1Q+9IQlwzKVEhaoiVe04QYdhQ4w7lnudl9nLN0k30MnBs01KdEwlx7aREhaghNqWtJ90np+wCC/VurlGHcs87PxFqqZcnzU/GczK/SOdEQlwbKVEhaoAShxPPXxa7DuX2qWGHcs/r2aQnZs1IusVMS2saCTuO6h1JiGsiJSpEDZCcnk1nEkk+d9uzmGY149SW3/Kx+HBb8G0AJHsr9m9crnMiIa6NlKgQNcD6tE0c8zmBQ9No7X+TW9727GqdnxC11MuTwKNLyS8q1TmREJUnJSqEm1NKYdw1nzjPc7NyWwzQOdG1uSPkDkyakf0WM60sKazcdVzvSEJUmpSoEG5uy2E77RxrSfIoO5Tbp1kfnRNdGx+LD90adQVgi3cROzau1DmREJUnJSqEm1u1cRvZ3kdxaBo3+7agmW8zvSNdsztD+wFlh3TrHVhEUalD50RCVE6lSnTKlCk0b94cm81GZGQkGzZs+N3xc+bMISwsDJvNRrt27Vi4cGG59UopJkyYQKNGjfDw8CA6Opq9e/eWG5OdnU1sbCy+vr74+/szcuRI8vPzXet3795Nr169CAwMxGaz0aJFC1577TVKSkoqs4lCuAWlFKXbf2KZlwcAd95Ysw/lnteraS9MmoF9FgutTEms3XtC70hCVEqFS3T27NmMGTOGiRMnsnHjRjp06EBMTAzHj1/6c41169YxbNgwRo4cSVpaGoMGDWLQoEFs27bNNebdd99l8uTJTJ06laSkJLy8vIiJiaGwsNA1JjY2lu3btxMXF8f8+fNZtWoVo0aNcq03m8088sgjLF26lN27d/PRRx8xffp0Jk6cWNFNFMJt7MnKJ7xoda05lHuer8WXqKBIALb7nGVzylqdEwlRSaqCIiIi1OjRo13fOxwOFRwcrCZNmnTJ8Q8++KAaMGBAuWWRkZHqySefVEop5XQ6VVBQkHrvvfdc63NycpTValXffPONUkqpHTt2KEAlJye7xixatEhpmqaOHDly2ax/+tOfVPfu3a962+x2uwKU3W6/6ucIcT1NXbhezXkvWIXPCFf3zb1b7zhV6se9P6rwGeFq0LRW6r8TR6iSUofekYRQSlWsCyr0TrS4uJjU1FSio6NdywwGA9HR0SQmJl7yOYmJieXGA8TExLjGp6enk5mZWW6Mn58fkZGRrjGJiYn4+/vTpUsX15jo6GgMBgNJSUmXfN19+/axePFievbsedntKSoqIjc3t9xDCHdSsOVn4s4dyo258W6d01StXx/SvdG4ng3p2XpHEqLCKlSiJ0+exOFwEBgYWG55YGAgmZmZl3xOZmbm744///VKYxo2bFhuvclkIiAg4KLX7datGzabjZYtW9KjRw/eeOONy27PpEmT8PPzcz1CQmruuXei9jlw8gwtzy53Hco9f2Pr2sLX4stt52bp7vLJJzXl0v8QF8Kd1brZubNnz2bjxo18/fXXLFiwgPfff/+yY8eNG4fdbnc9MjLkrhLCfSzftIt8n0NlF1jwbUFT36Z6R6pyMS3uAmCJlyfmPT/jdCqdEwlRMaaKDG7QoAFGo5GsrKxyy7OysggKCrrkc4KCgn53/PmvWVlZNGrUqNyYjh07usb8duJSaWkp2dnZF73u+XeTbdq0weFwMGrUKP785z9jNBovyma1WrFarVfabCF0kbf5Z9Z6nzuUe9M9Oqe5PnqF9MKiGUm3QKghkbSMHDo3q6d3LCGuWoXeiVosFjp37kx8fLxrmdPpJD4+nqioqEs+Jyoqqtx4gLi4ONf40NBQgoKCyo3Jzc0lKSnJNSYqKoqcnBxSU1NdYxISEnA6nURGRl42r9PppKSkBKfTWZHNFEJ3x+wFNMtLuHDbs+Y167ZnV8vb4s1tjcp+z/d455Kccuk5DkK4qwq9EwUYM2YMw4cPp0uXLkRERPDRRx9x5swZRowYAcAjjzxC48aNmTRpEgAvvPACPXv25IMPPmDAgAHMmjWLlJQUpk2bBoCmabz44ou8+eabtGzZktDQUMaPH09wcDCDBg0CoHXr1vTt25cnnniCqVOnUlJSwrPPPsvQoUMJDg4GYObMmZjNZtq1a4fVaiUlJYVx48YxZMgQzGZzVewrIarN8rQ9KO8DODV/wv1uIsSn9n5e3/fGu1l+dA1LvD25e9dPKBWDpml6xxLiqlS4RIcMGcKJEyeYMGECmZmZdOzYkcWLF7smBh06dAiD4cIb3G7duvH111/z2muv8corr9CyZUt+/PFHwsPDXWNeeuklzpw5w6hRo8jJyaF79+4sXrwY27k7VkBZST777LP07t0bg8HAfffdx+TJky9siMnEO++8w549e1BK0axZM5599ln+9Kc/VWrHCKGnU2k/k+Z9bkLRTQN1TnN99QzpiVUzctAMTdQath3JpV0TP71jCXFVNKWUfJJ/Tm5uLn5+ftjtdnx9ffWOI+qorNxC1n14JxObZ6M0jSX3LSHYO1jvWNfVmGWjiTuyisdz7FibfspTg2vHRSVEzVSRLqh1s3OFqOniN+4m3ycdpWm0929Z6wsU4M6bys6BXezlibbje+Tf9qKmkBIVws2c3vgDS12zcgfpG6aa3N74dmyaicNmM8GsYesRu96RhLgqUqJCuJFMeyHBeUvZZLOiAX3P3cC6tvM0e9Ir5A4AdvjksT5pnb6BhLhKUqJCuJGE1O2c8DkMQET9djT0bHiFZ9QeA869617s5YVhx1w5pCtqBClRIdxIbtoPLDl3KLdfq/t0TlO9ugV3w9do46TJSH3TWjZn5OgdSYgrkhIVwk0csxfQ4OxS9lgtmDAQ3Sz6yk+qRcxGM3c2K7s+8CafQlI2rNY5kRBXJiUqhJtYnrKVIz5ll8jsHnQrfta6d65k/5Zl58TGeXpi2CmHdIX7kxIVwk2cSfuORd6eAAxodb/OafRxS8NbaGj2Ic9owMucyKZDp/WOJMTvkhIVwg0cySnApyieI2YTHpqJniGXvw9ubWY0GOl77s4uyT4lpCYm6JxIiN8nJSqEG0hYn0q6zykA/tDkdjxMHjon0s/5Q7orPD0w7P0eh9weTbgxKVEh3EDRpm9Y4lV2KLd/HZuV+1ttAtrQzNaAIoMBiy2Z9fuOX/lJQuhESlQIne06ZsdHLeeUyYif0UZU8KVvK1hXaJrGgJZl/5BY4w3b183XOZEQlyclKoTOEtetZKNvIQD9Q/tjNsit++46d+eaRA8bHoe+p7DEoXMiIS5NSlQIHTmdCnbOJsGz7DPQe25+QOdE7iHEN4SOvjfh1DTOeO9k1faDekcS4pKkRIXQUeqBk2i2JAoNBprbbqBt/bZ6R3Ib97T5IwBLfCwcXPedzmmEuDQpUSF0tG3tPFZ6awAMvPlBNE3TOZH7uLP5nZgxsM9iwZr9M/aCEr0jCXERKVEhdFLicGI4OIcUDxsacFfLQXpHcit+Vj96BUUCcNT3MCtSd+icSIiLSYkKoZO1Ow9h99oNwK3+YQR5BemcyP0MbPsQAIu8PTmZ8rXOaYS4mJSoEDrZv/Y7lvpYgAtlIcqLCo6insGDbKMRQ+FSMu2FekcSohwpUSF0YC8oQTv5AwfNZmwYiW7WR+9IbslsMDOgRX8AtvvYWb52rc6JhChPSlQIHcRvSCPD9ygA0Y174Gn21DmR+7onbAgAyz09Kdg8Q+7sItyKlKgQOsjZ8KXrji0D28bqnMa9hQWEcZMtkGKDRqklkbSDJ/WOJISLlKgQ1WxfVh5OFU++wUCw2Y+IoAi9I7k1TdO4v+3DACz1MbJl5fc6JxLiAilRIarZ+pULWO1XCsB9YUMxaPJreCV3tRyEGQO7rRacR2bJZQCF25DfXiGqkcOpKN73PzbabBgUDJTL/F0VP6sf0UFdAcjwOcjyjXLOqHAPUqJCVKO1Ow5w3GsnALfVb0egV6DOiWqO+9o/BsBib0+OJc7QN4wQ50iJClGN9q/+H4t8bADc3/5xndPULLcG3UqwyZd8g4HikiUcPX1W70hCSIkKUV3sZ0sozP2JbKORAM1Gj5AeekeqUQyagfvPne6S6FvEmlVxOicSQkpUiGqzbM1aNvvaARjccqDcN7QSBrYeikHBRpsN+87Pym4lJ4SOpESFqAZKKU6nTWWtR9mh3MHnTtkQFdPQsyG3BYQDcMJjG2t3yn1Ghb6kRIWoBin7jnHamozSNG7xbkEz32Z6R6qxHuw4CoCFPlb2Lf9E5zSirpMSFaIa7Ez4nPm+VgBib3la5zQ1W48mtxNo9MJuNFJY8BNZ9gK9I4k6TEpUiOss+0wx+fbvyDYaqa950KtZb70j1WhGg5GhbcoulbjKr4iVCYt0TiTqMilRIa6z5SviWe+XD8CQsAdlQlEVuLdNLCalsc1q5fTuaThkgpHQiZSoENeRUoqT2/5Dms2GUcH94cP1jlQrBNgC6BMYCUC6917Wbd2rcyJRV0mJCnEdJe0+SIZn2SXqejXoxA2eN+icqPaI7TwagDgvD/at+o/OaURdJSUqxHW0K+G/LPU+N6Goy3M6p6ld2t/QgZbmGyg2aJySKxgJnUiJCnGdHMs5S3bBPAoMBpqb6tE5sIvekWoVTdN4uNMTACzzVaxcPFvnRKIukhIV4jpZuXgOCeduefZw+8fQNE3nRLVP35aD8MbEEbOJ7IzpFBTLLdJE9ZISFeI6KCxxcOLQfzloNuOFibvCHtQ7Uq3kYfLggRvvASDF9wTxa9bonEjUNVKiQlwH8avXkOp7AoAhLe7B0+ypc6La64+dnsaoIMXDxuHUD1FKTncR1UdKVIgqppTiYOo/SPWwYVIXZpGK6yPIK4g+Dco+b97nsZ0NO37ROZGoSypVolOmTKF58+bYbDYiIyPZsGHD746fM2cOYWFh2Gw22rVrx8KFC8utV0oxYcIEGjVqhIeHB9HR0ezdW/68r+zsbGJjY/H19cXf35+RI0eSn5/vWr9ixQoGDhxIo0aN8PLyomPHjsycObMymyfENUne8Qv7PLcDEH3DrTT0bKhzotrvsa5/BWCZt41t8f/QOY2oSypcorNnz2bMmDFMnDiRjRs30qFDB2JiYjh+/Pglx69bt45hw4YxcuRI0tLSGDRoEIMGDWLbtm2uMe+++y6TJ09m6tSpJCUl4eXlRUxMDIWFha4xsbGxbN++nbi4OObPn8+qVasYNWpUuddp374933//PVu2bGHEiBE88sgjzJ8/v6KbKMQ12Rr/DxK8yk5rebzrSzqnqRtaN2hDR1sTHJrGEbWMg8dz9I4k6gpVQREREWr06NGu7x0OhwoODlaTJk265PgHH3xQDRgwoNyyyMhI9eSTTyqllHI6nSooKEi99957rvU5OTnKarWqb775Riml1I4dOxSgkpOTXWMWLVqkNE1TR44cuWzW/v37qxEjRlz1ttntdgUou91+1c8R4tcOZJ1WEz9qpcJnhKtHZvXTO06dsvLAMhU+I1xFfdpGffX5O3rHETVYRbqgQu9Ei4uLSU1NJTo62rXMYDAQHR1NYmLiJZ+TmJhYbjxATEyMa3x6ejqZmZnlxvj5+REZGekak5iYiL+/P126XDjPLjo6GoPBQFJS0mXz2u12AgICLru+qKiI3Nzccg8hrsWqBf9lia8JgJGRf9Y5Td3SvWkvGmve5BkNHDk1k5wzRXpHEnVAhUr05MmTOBwOAgMDyy0PDAwkMzPzks/JzMz83fHnv15pTMOG5T9XMplMBAQEXPZ1v/32W5KTkxkxYsRlt2fSpEn4+fm5HiEhIZcdK8SVnMgt5Njpr8g3GAgx+NC9WS+9I9UpBs3AyHP3Gl3mX8zShXLxBXH91crZucuXL2fEiBFMnz6dtm3bXnbcuHHjsNvtrkdGRkY1phS1TfyCGSz0KzvZf1TnZzFotfLXy63d0/aP1MPCMZOJjIMfc7a4VO9Iopar0G95gwYNMBqNZGVllVuelZVFUFDQJZ8TFBT0u+PPf73SmN9OXCotLSU7O/ui1125ciV33303//jHP3jkkUd+d3usViu+vr7lHkJURl5BMYeOTuOUyUhDbAwIe0DvSHWS1Wjl0TYPAbDcP49lcTKxUFxfFSpRi8VC586diY+Pdy1zOp3Ex8cTFRV1yedERUWVGw8QFxfnGh8aGkpQUFC5Mbm5uSQlJbnGREVFkZOTQ2pqqmtMQkICTqeTyMhI17IVK1YwYMAA3nnnnXIzd4W43uKWfE+cfwEAT3R4Qu4ZqqOhHZ/ERxk5aDazd+cHlDicekcStVlFZy3NmjVLWa1WNWPGDLVjxw41atQo5e/vrzIzM5VSSj388MNq7NixrvFr165VJpNJvf/++2rnzp1q4sSJymw2q61bt7rGvP3228rf31/99NNPasuWLWrgwIEqNDRUFRQUuMb07dtXderUSSUlJak1a9aoli1bqmHDhrnWJyQkKE9PTzVu3Dh17Ngx1+PUqVNXvW0yO1dURkFxqXr3/UgVPiNc9fi8kyosLdQ7Up33rzVvqPAZ4WrQtFZqUUK83nFEDVORLqhwiSql1Mcff6yaNm2qLBaLioiIUOvXr3et69mzpxo+fHi58d9++61q1aqVslgsqm3btmrBggXl1judTjV+/HgVGBiorFar6t27t9q9e3e5MadOnVLDhg1T3t7eytfXV40YMULl5eW51g8fPlwBFz169ux51dslJSoqY/6S+ar/9JtV+Ixw9en6D/WOI5RS9iK7ivi8vQqfEa4+fD9aORxOvSOJGqQiXaApJReaPC83Nxc/Pz/sdrt8PiquSqnDyUcf9eKLBtn4KCPLYtfJdXLdxHsJL/NlxkLaFBbzTMf/0TMyQu9IooaoSBfI9EEhrkHcinjW+JadZhV7471SoG5kZLeXsSjYYbOQtmYiTqe8XxBVT0pUiEoqcTjZvukNfrFY8FIGHo58Ue9I4lcCbAHcG9wbgHU++1mdtF7nRKI2khIVopKWJSwk3q/sdmfDb3oQX4t8BOBunr59Ah5KY6fVQuq68fJuVFQ5KVEhKqG41MnmrW9x2GzGTxkZHvknvSOJSwiwBTCk6QAAVvoeZsXaVTonErWNlKgQlbA07gfi/HIAeKL1cPks1I090X0c3k4D+y1mUpIn4pB3o6IKSYkKUUFFpQ427Xqb4yYTDZSZYV3kptvuzNfiy8M33Q/Act9MElbF6ZxI1CZSokJU0MKF37DU7wwAT7V/EovRonMicSWPdv0z/k4jh81mUtL+T65iJKqMlKgQFZBbUMTm/e9x2mikkbJxb8fH9I4kroKn2ZPHwh4GIM43m0WLv9U5kagtpESFqICfvvuQhX5ldwZ5MeLPco3cGiQ24jkaKisnTCY27p1EXkGx3pFELSAlKsRVOnryNGn2LykwGAgzBNCv9RC9I4kKsBgt/DnirwAs8Cvlx+8/1DmRqA2kRIW4Sj999wrLvI0AvNb7XTRN0zmRqKh+rR8kzBBAocFAas4XHDuVo3ckUcNJiQpxFXb+sp9ElqM0jTs8b6ZDcOSVnyTcjqZpjO/9LgDxPiZ+/G6szolETSclKsQVKKVYuOBPpHmYsSh4pe9HekcS16B9cCR3eLQCYJ1awa5fDugbSNRoUqJCXMGKdSuI99wDwANBvWnk00TnROJavdpvMhYFmzzMzJv/LHIzK1FZUqJC/I6ColKWp/yVDLMJf6eR53q/pXckUQWCfBozpFEfAJZ4/8LS5fN1TiRqKilRIX7HnDlvs8C3EIA/d3oRL7OXzolEVXmu91s0dJrJMplYsW08+YVyyouoOClRIS7j4NFjrMybSbFBo53hBgZ2GK53JFGFPEwejO36KgCLfUuZPWuCzolETSQlKsRlzJn7DBs8TZiV4s3+/5ZTWmqhPq3vI8IcQqmmEV/4I/sOHNQ7kqhhpESFuITla+NY7LELgAdu+AMt6ofpnEhcL2/eNRWrE7Z6mPn25ydlkpGoEClRIX4j72whi1P/SpbJREOHiT/FvKt3JHEdNfJtyiNN7gFgsfch5i3+WudEoiaREhXiN/43cwyLfMuuj/tK19ewmWw6JxLX29O9X6ep08Zpo5HF+//O8ezTekcSNYSUqBC/kpyWxIJzVybqZWtF7zb36R1JVAOzwcykPh9jUIrV3ga+/OYJvSOJGkJKVIhzCotL+G7NMxyymKjnMPC3gZ/pHUlUo/ZNunJ/QA8A5tl2sDThB50TiZpASlSIc774+mUW+xQBMLbzWPxsfjonEtXt5QH/pInDQrbJyE87J5CTl693JOHmpESFANK2pjCvZBFOTeN2cwv6dximdyShA4vRwlu9/4FBKVZ5w6dfPa53JOHmpERFnXemoJD/rXiCgxYT9Rwafx/8ud6RhI5uaXY7g/2jAPjRuoV5S2fqnEi4MylRUedN+/IJ4s7Nxh1/63j8PAJ0TiT09srdU2jusJFjNDJn/1sczjymdyThpqRERZ22KH42c82pANztdQt92j2gcyLhDixGCx8O+BSbU5HmYWTanGE4nHIRBnExKVFRZ2WeOM7Xe98gx2ikucPK64Om6R1JuJGWge15unksAPN8TvLl7Nf1DSTckpSoqJMcDidTZj/IJg8DVqfi/b6fYDFZ9Y4l3MyIO8YSoQVRqmnMyf+WlM2JekcSbkZKVNRJ074awzzvkwA82XQINwd31DeQcEuapvHB/d/QwKGRYTExbe0osk/n6B1LuBEpUVHnxK36jlmlS3FoGj2MzXj8D6/pHUm4MX/PBvy9+/uYlCLRC/45816cDqfesYSbkBIVdcrBw+l8snMC2SYjTUvNfDDkW7nFmbiiqFZ3MqrRYAB+8j7OtK//qnMi4S6kREWdUVRSwkc/DmGHzYiXUzF5wJd4mD31jiVqiKfufIMehhAcmsY3xYuIXy2XBRRSoqKOUErxwadDWOZTAMC4Ni9wY1C4zqlETaJpGh8MmUOzUjPZJiP/2fkq+/bv0juW0JmUqKgTpn/zMnNsewC43zuCgZFylw5RcR4WLz7qNwMfh2K31cg7i4dw2p6jdyyhIylRUevNW/Y5XxQuoFTT6EoQE+79RO9Ioga7Kbg9b3Z5HZNSrPdyMumruyktdegdS+hESlTUamnbE/l3+nvkGg20LLXy8bCfZCKRuGZ/aH8/oxsPAWCRdw4ffP6QzomEXqRERa118PAvvLPmCQ5bjDQshX/fNxebRSYSiarxeJ/xDLR1AGCWeSvTZo3VOZHQg5SoqJVOZp/g9Xn3st2m4eVUvNfzXwT5N9U7lqhl3njwSyKcDSjVND4tmMes+R/qHUlUMylRUeucPVvA+Fn9SfF0YnEq3ujwGrfc1FPvWKIWMmgGpjy0gDalnpw1GPj38U9ZuEJunVaXSImKWqWkpITXvriTNV6FGJXi5Ruf5M5bhuodS9RiNrMn04ctokWJidNGAx/ue4u1KYv1jiWqSaVKdMqUKTRv3hybzUZkZCQbNmz43fFz5swhLCwMm81Gu3btWLhwYbn1SikmTJhAo0aN8PDwIDo6mr1795Ybk52dTWxsLL6+vvj7+zNy5Ejy8/Nd6wsLC3n00Udp164dJpOJQYMGVWbTRA1WWlrK+E/6EeeZA8BTDe/lwduf0zeUqBN8PQP4730/ElyikWU28lbaGJI3r9Q7lqgGFS7R2bNnM2bMGCZOnMjGjRvp0KEDMTExHD9+/JLj161bx7Bhwxg5ciRpaWkMGjSIQYMGsW3bNteYd999l8mTJzN16lSSkpLw8vIiJiaGwsJC15jY2Fi2b99OXFwc8+fPZ9WqVYwaNcq13uFw4OHhwfPPP090dHRFN0vUcA6Hk/Gf9GOBZxYAD3nfzlP939A5lahLguo141/9ZnJDKWRYjLy+4RlSt67SO5a43lQFRUREqNGjR7u+dzgcKjg4WE2aNOmS4x988EE1YMCAcssiIyPVk08+qZRSyul0qqCgIPXee++51ufk5Cir1aq++eYbpZRSO3bsUIBKTk52jVm0aJHSNE0dOXLkotccPny4Gjhw4BW3pbCwUNntdtcjIyNDAcput1/xucJ9OEodatx/7lThM8JV+Ixw9ea3T+gdSdRh2w8kq16flP1Z7D+tjUrdtlrvSKKC7Hb7VXdBhd6JFhcXk5qaWu6dnsFgIDo6msTES99nLzEx8aJ3hjExMa7x6enpZGZmlhvj5+dHZGSka0xiYiL+/v506dLFNSY6OhqDwUBSUlJFNqGcSZMm4efn53qEhIRU+mcJfZSWlvLaJ/2Z53EUgKEekbz6gNxcW+inTbMu/Kv3pzQoVRyyGBi/7kk2bpN3pLVVhUr05MmTOBwOAgMDyy0PDAwkMzPzks/JzMz83fHnv15pTMOGDcutN5lMBAQEXPZ1r8a4ceOw2+2uR0ZGRqV/lqh+RUVFvPRJb+bZjgAwxBbBqw/K1YiE/to0j+Bff7hQpGPXP82q5Hl6xxLXQZ2enWu1WvH19S33EDVDXr6dFz+7nTiPbABiPbvx2pBPdU4lxAVtQyOZ8ofPCCqBY2YDE7aMZeHKL/WOJapYhUq0QYMGGI1GsrKyyi3PysoiKCjoks8JCgr63fHnv15pzG8nLpWWlpKdnX3Z1xW114nTmTz/VS/WeJ7FqBRP1RvA2Af+q3csIS7SJjSC6Xd9R9NijVMmA2/+8g7fLv5I71iiClWoRC0WC507dyY+Pt61zOl0Eh8fT1RU1CWfExUVVW48QFxcnGt8aGgoQUFB5cbk5uaSlJTkGhMVFUVOTg6pqamuMQkJCTidTiIjIyuyCaKG2/HLRp75tg8pHiVYnIq/BD/C6Hve1juWEJfVPOhmPr9/IS2LTeQZDbx77BP+NftFvWOJqlLRWUuzZs1SVqtVzZgxQ+3YsUONGjVK+fv7q8zMTKWUUg8//LAaO3asa/zatWuVyWRS77//vtq5c6eaOHGiMpvNauvWra4xb7/9tvL391c//fST2rJlixo4cKAKDQ1VBQUFrjF9+/ZVnTp1UklJSWrNmjWqZcuWatiwYeWybd++XaWlpam7775b3XHHHSotLU2lpaVd9bZVZEaWqH4J679Xfaa3UeEzwlXXT9uqH1dP0zuSEFfNnn9KxU6LcM0iH//JQFVaUqp3LHEJFemCCpeoUkp9/PHHqmnTpspisaiIiAi1fv1617qePXuq4cOHlxv/7bffqlatWimLxaLatm2rFixYUG690+lU48ePV4GBgcpqtarevXur3bt3lxtz6tQpNWzYMOXt7a18fX3ViBEjVF5eXrkxzZo1U8BFj6slJeq+vpo/SUV9Wlag0dPbqg3bl+kdSYgKKy4pVC981tdVpE9N7aZy807rHUv8RkW6QFNKKd3eBruZ3Nxc/Pz8sNvtMsnITThKHbz1VSxz2YZD02hVbOIf93xL08CWekcTotLemT2KrwvW4dQ0woo0JvaaTnhL+WjKXVSkC+r07Fzh3o6fOsJTn97GHG07Dk0jqrgeXzy0QgpU1HgvD5nGS8GP4uV0ssuqGL3qMb5fNkXvWKISpESFW0rcvISRc2NYbzuDUSliLRH8d+QKvD389I4mRJWIvfMv/Kvrv2haDNkmA28e/g9vffkwTodD72iiAqREhVtRSvHRt8/z4sYxHLBo1Ct18nro84wd9imaQf64itqlS+te/G/IMm4t8qFU05ilNvHYJ105eHSP3tHEVZK/lYTbyDp1hKemdefTguWcNRgIKzLx397/Y1DPUVd+shA1VIBvIJ+MXM0QU2dMSpFqK2TEosF8H/9vvaOJqyATi35FJhbp56cV0/j3vskcNWsYlOIubub1P36N2WLVO5oQ1WbJ+m/4cNtbrt+Dvs5QJgz9Gi9PH72j1SkV6QIp0V+REq1+9vxs3pgdS5wxA6Vp3FDq5LkWzzK419N6RxNCF8ezjzL+uwdZZ7UD0LwYXugwgeiIB3ROVndIiVaSlGj1Wrj2C/61430yLGXfdyvy47WB/yMkMFTfYEK4gf/8MI6Zp3/GbjRgVIoB3MyrQ77A08Nb72i1npRoJUmJVo+sUxm89cMIVpgyUVrZ5KHhN9zPyHv+T+9oQriVPQe38PfFj5NqKwAgpBhGtXpO5glcZ1KilSQlen0pp5Pp8yfw9YkfOGUqm9MWVejDK3d9TvPGN+ucTgj35HQ4mDL3L3ybt5QcY9nvTY/i+rw68DMaN2yhc7raSUq0kqREr5/Vm+YxZcNEtltLAGhc4uThxo8S2++vOicTomY4eGQ378x/nNW2HAB8HU4Get7Gi/dOxmKx6RuulpESrSQp0ap3+PgvvPfzU6wwHcOpaVicij7OUMY+MAN/3/p6xxOixpkT9y++ODCVgxYNgKbFMOKmp7i/12idk9UeUqKVJCVadez5p/nHD8+wpGQL+ecOQd1aaOPJbm8T2a63zumEqNnOFpzho++eZkFJKrnnfr/aF9l4KuI1enQcqHO6mk9KtJKkRK9dUXEBU376Kz/lriDbVPYv5dBixdAmw/mjHLoVokr9krGDfy56htWWk5RqGppSdC0J4Pk/vEv4jV31jldjSYlWkpRo5RUVF/CfeWOZnxNP1rnyDCpx0t/7Dp659wOs8pmNENfNmo0LmLHh/0jyKJvFa1SKbqUNefqOv9PuJinTipISrSQp0Yo7W5jHf+a/yiJ7gqs865U6iTa257nBH1PPr4HOCYWoO35e/imz93zMFlvZReyNShFVegOjbn+dTq166pyu5pASrSQp0at34vQRpix8mYSiTZw2lpVnQKmTPxja8tRdHxF4Q7DOCYWom5TTyXdxU/jpwKdsPlemmlJ0LvHlj52ep0+XoTondH9SopUkJXplW35Zy+er3mIdhzhrKCvPhqVOehrCefKuf0h5CuEmysr038xP/4yNHiWu5WHFZvo3uZeH+vwVs0muTX0pUqKVJCV6aU6ng7mrp/Ljni/ZYj6D0srKs3mxkzs8e/D43ZPw862nc0ohxOUsWfMNP2ybTJItj9Jzv783lCp6eXTm8Tv/j0YNmusb0M1IiVaSlGh5h7L28vmyN1hTsIlM84XlHQtM9AoayEP9X8FisegXUAhRIalbVzIrcRLrTRmuqx+ZlaJzaQCD2j5G/67D0c6VbF0mJVpJUqJQXFrInOUfE3fgBzaZc3Gc+4XydjiJLG7AwI4v0KvrvTqnFEJci6yTx/hy8WusO5vEPuuF0gwu0ejm2Z6H7hjLjcHhOibUl5RoJdXVElVKsWLTXH7e9Bkp6iA5xgu/VDcXQqRnBLF3TiQ4sKmOKYUQVc3pcPLj8mnE7/sfydYcCgxl7041pQgv8aRHYG/++Ic/4+ddt2bZS4lWUl0qUaUUKbvj+WHDVJKLd5c7XBtQ6qRzaSC9bx5Ovx4PYzh32EcIUXsdOLKXb+MnkVqQwg7bhVqwOBUdSv25o1l/7r/9OTxttf8G4VKilVTbS1QpxdptC1iQ9jkbi/dy1Hzhf73N6aRjkSddGkQz5M6X8JeJQkLUSUopViXPY9GmqWzRDpJhufCPaJtT0a7Un9sa9+GBO57D1zNAx6TXj5RoJdXGEi0qPsu89Z+zat88tqojnDRdWGdxKsKLzHT0ieLeO/5Cs8ZyWyUhxAXFJSUsWPkZa36ZzRZjJplmo2udWSnCSjy5pV4E9972NC0atdUxadWSEq2k2lKiew9v4uekT9mUncwuUz6FhgufcdqcTtoV2WjnE8nd3Z/npmZhOiYVQtQUhUWFLFj5KevTf2Cr8ShHflWoAE1LDLQ1N6dnq8H06TIUi7nmXupTSrSSamqJnso9xoL1M0g5vIJdzqMcM5df36DUQduSerSrfzt33/60TBASQlyTkpIS4td/x7rds9np+IXdVuU6fxzAw6loVepFW98O9Ok0lM4te9WoU2ekRCupppToCftRlibPJPXwSvaWZnDQ5Cj3B9ioFDcXadxkCOXW0LuIue1hPDw8dEwshKitlFJs253CstTP2ZWXwk7LGU7/ZjKin0PRyulPG/8O9Gp/P51a9sSgue+ERSnRSnLHElVKsfNQMiu3/MD2E6nsd2aSYb74f1nTYgc3ORrQOiCSXl1GcHOLNjqkFULUdfln84lPnEXagYX8UvILu6ylFBrKF6a3Q3Gjw5ubvW+mS4ve9Ow4CE+re/ydC1KileYOJXrCfpQVm+ay9fA60s/uJ92Yj9148WGQkGIHLRz1aOnTgW7hD9I5/HY5FUUI4XayThwjPmkm24+tIN15iL0Wx0WlalSKZiUmmhuDaFW/A93a9Kd9i9swGk2X+anXl5RoJVV3iR7LPsjabfPZfng9B87uJ0PLJesSf2YsTkWLYo0Q1YCb6nUmqt19dLg5UkpTCFHjnDh1nBXJ37L9yAoOlaSTbi7gpMl40TgPp6JZqZUm5iBaBbQnIqwP7W/sjtl4/S81KiVaSderRAtLzrJpz0o27V/N/uwdHC4+ymHjWdctxH4rpNhBE4c3jS3NCGvUnR633E9wYOMqyyOEEO6iqKiYDVuWkbp3AQdzd3BUO8l+i/Oid6tQdlpNkxITjQ0BNPW+kbDgLnRtG0OjgOZVmklKtJKqokS37U9i3Y75pGfv5FjRUTK1fDJNTtc1aH9NU4rGJU6alHoRbAkhtH4nbm0zgNY3dpR3mUKIOivbnk1i2s/szFjD4TN7Oaplc8ji5MwlihXKJi4FO6wEGuvTxLsFYcG3MKDrCEwm8yXHX4mUaCVVRYmO++we5hvTL1ru43ASUmIgUPkT5NGcmxreSpe2fQkNaVmjpn4LIYQe8vLzSNm2jO0HV5Fh38Vx53EyTQUcMRnKnZ0A4Ol0su7hLRgvcZj4alSkC/T51LYWa1GvPa2Pp9NQ+XKDpTEh9drQunl3Orbujoet5p58LIQQevLx9qFX18H06jrYtUwpxcHD+9i4K54DWWkcO7OfE+oURoyVLtCKkneiv+IOs3OFEELoqyJdIB+8CSGEEJUkJSqEEEJUkpSoEEIIUUlSokIIIUQlSYkKIYQQlSQlKoQQQlSSlKgQQghRSXKxhV85f8psbm6uzkmEEELo5XwHXM1lFKREfyUvLw+AkJAQnZMIIYTQW15eHn5+fr87Rq5Y9CtOp5OjR4/i4+NzTdezzc3NJSQkhIyMDLny0W/Ivrk82TeXJ/vm8mTfXF5l941Siry8PIKDgzFc5qL358k70V8xGAw0adKkyn6er6+v/KG+DNk3lyf75vJk31ye7JvLq8y+udI70PNkYpEQQghRSVKiQgghRCVJiV4HVquViRMnYrVa9Y7idmTfXJ7sm8uTfXN5sm8urzr2jUwsEkIIISpJ3okKIYQQlSQlKoQQQlSSlKgQQghRSVKiQgghRCVJiQohhBCVJCVaxaZMmULz5s2x2WxERkayYcMGvSNVu0mTJnHrrbfi4+NDw4YNGTRoELt37y43prCwkNGjR1O/fn28vb257777yMrK0imxft5++200TePFF190LavL++bIkSM89NBD1K9fHw8PD9q1a0dKSoprvVKKCRMm0KhRIzw8PIiOjmbv3r06Jq4eDoeD8ePHExoaioeHBzfeeCN/+9vfyl0gvS7tm1WrVnH33XcTHByMpmn8+OOP5dZfzb7Izs4mNjYWX19f/P39GTlyJPn5+RUPo0SVmTVrlrJYLOqzzz5T27dvV0888YTy9/dXWVlZekerVjExMerzzz9X27ZtU5s2bVL9+/dXTZs2Vfn5+a4xTz31lAoJCVHx8fEqJSVFde3aVXXr1k3H1NVvw4YNqnnz5qp9+/bqhRdecC2vq/smOztbNWvWTD366KMqKSlJ7d+/Xy1ZskTt27fPNebtt99Wfn5+6scff1SbN29W99xzjwoNDVUFBQU6Jr/+3nrrLVW/fn01f/58lZ6erubMmaO8vb3VP//5T9eYurRvFi5cqF599VU1d+5cBagffvih3Pqr2Rd9+/ZVHTp0UOvXr1erV69WN910kxo2bFiFs0iJVqGIiAg1evRo1/cOh0MFBwerSZMm6ZhKf8ePH1eAWrlypVJKqZycHGU2m9WcOXNcY3bu3KkAlZiYqFfMapWXl6datmyp4uLiVM+ePV0lWpf3zcsvv6y6d+9+2fVOp1MFBQWp9957z7UsJydHWa1W9c0331RHRN0MGDBAPfbYY+WW3XvvvSo2NlYpVbf3zW9L9Gr2xY4dOxSgkpOTXWMWLVqkNE1TR44cqdDry+HcKlJcXExqairR0dGuZQaDgejoaBITE3VMpj+73Q5AQEAAAKmpqZSUlJTbV2FhYTRt2rTO7KvRo0czYMCAcvsA6va++fnnn+nSpQsPPPAADRs2pFOnTkyfPt21Pj09nczMzHL7xs/Pj8jIyFq/b7p160Z8fDx79uwBYPPmzaxZs4Z+/foBdXvf/NbV7IvExET8/f3p0qWLa0x0dDQGg4GkpKQKvZ7cxaWKnDx5EofDQWBgYLnlgYGB7Nq1S6dU+nM6nbz44ovcdttthIeHA5CZmYnFYsHf37/c2MDAQDIzM3VIWb1mzZrFxo0bSU5OvmhdXd43+/fv5z//+Q9jxozhlVdeITk5meeffx6LxcLw4cNd23+p37Havm/Gjh1Lbm4uYWFhGI1GHA4Hb731FrGxsQB1et/81tXsi8zMTBo2bFhuvclkIiAgoML7S0pUXFejR49m27ZtrFmzRu8obiEjI4MXXniBuLg4bDab3nHcitPppEuXLvz9738HoFOnTmzbto2pU6cyfPhwndPp69tvv2XmzJl8/fXXtG3blk2bNvHiiy8SHBxc5/eN3uRwbhVp0KABRqPxolmUWVlZBAUF6ZRKX88++yzz589n+fLl5e7TGhQURHFxMTk5OeXG14V9lZqayvHjx7nlllswmUyYTCZWrlzJ5MmTMZlMBAYG1tl906hRI9q0aVNuWevWrTl06BCAa/vr4u/YX//6V8aOHcvQoUNp164dDz/8MH/605+YNGkSULf3zW9dzb4ICgri+PHj5daXlpaSnZ1d4f0lJVpFLBYLnTt3Jj4+3rXM6XQSHx9PVFSUjsmqn1KKZ599lh9++IGEhARCQ0PLre/cuTNms7ncvtq9ezeHDh2q9fuqd+/ebN26lU2bNrkeXbp0ITY21vXfdXXf3HbbbRedCrVnzx6aNWsGQGhoKEFBQeX2TW5uLklJSbV+35w9exaDofxf10ajEafTCdTtffNbV7MvoqKiyMnJITU11TUmISEBp9NJZGRkxV7wmqZFiXJmzZqlrFarmjFjhtqxY4caNWqU8vf3V5mZmXpHq1ZPP/208vPzUytWrFDHjh1zPc6ePesa89RTT6mmTZuqhIQElZKSoqKiolRUVJSOqfXz69m5StXdfbNhwwZlMpnUW2+9pfbu3atmzpypPD091VdffeUa8/bbbyt/f3/1008/qS1btqiBAwfW2tM4fm348OGqcePGrlNc5s6dqxo0aKBeeukl15i6tG/y8vJUWlqaSktLU4D68MMPVVpamjp48KBS6ur2Rd++fVWnTp1UUlKSWrNmjWrZsqWc4uIOPv74Y9W0aVNlsVhURESEWr9+vd6Rqh1wycfnn3/uGlNQUKCeeeYZVa9ePeXp6akGDx6sjh07pl9oHf22ROvyvpk3b54KDw9XVqtVhYWFqWnTppVb73Q61fjx41VgYKCyWq2qd+/eavfu3TqlrT65ubnqhRdeUE2bNlU2m021aNFCvfrqq6qoqMg1pi7tm+XLl1/y75jhw4crpa5uX5w6dUoNGzZMeXt7K19fXzVixAiVl5dX4SxyP1EhhBCikuQzUSGEEKKSpESFEEKISpISFUIIISpJSlQIIYSoJClRIYQQopKkRIUQQohKkhIVQgghKklKVAghhKgkKVEhhBCikqREhRBCiEqSEhVCCCEq6f8Ba4LErguRU2UAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r = np.linspace(-10, 10, 100)\n",
    "h = ml.headalongline(r, 0, 0.01).squeeze()\n",
    "plt.plot(h[0])\n",
    "plt.plot(h[1])\n",
    "plt.plot(h[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x15b85d4c0>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAAEWCAYAAAD2AJlUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2uElEQVR4nO3de1iUdf4//ufMwMwwCAM4wgCiiFpIClOQhNnmbiSYtbpbLVafVfn2tV3L3frRZlEGdmVf1FyXrUw2+2rqdnDbbft+rtalXIo+HVAL81Ae0oI4OcNBYTjOwMz9+wNmcBKUwRnuOTwf13Vfys09N6+ZK3n2vt8niSAIAoiIiPyIVOwCiIiIxhrDj4iI/A7Dj4iI/A7Dj4iI/A7Dj4iI/A7Dj4iI/A7Dj4iI/A7Dj4iI/A7Dj4iI/A7Dj4iI/M6owm/Lli2Ij4+HUqlEeno6Dh48OOy177zzDtLS0hAWFobg4GDodDrs3r3b4Zrly5dDIpE4HNnZ2aMpjYiI6LICnH3Bnj17kJeXh5KSEqSnp6O4uBhZWVk4deoUIiMjL7o+IiICTz31FBITEyGXy/Hee+8hNzcXkZGRyMrKsl+XnZ2NHTt22L9WKBQjrslqtaKhoQEhISGQSCTOviUiIvIRgiCgvb0dMTExkEov0b4TnDR79mzhoYcesn9tsViEmJgYoaioaMT3uPbaa4U1a9bYv162bJmwaNEiZ0uxq62tFQDw4MGDBw8eAgChtrb2krnhVMvPbDajsrIS+fn59nNSqRSZmZmoqKi47OsFQcCHH36IU6dOYcOGDQ7fKy8vR2RkJMLDw/Gzn/0M69atw/jx44e8j8lkgslkcrgvANTW1iI0NNSZt0RERD7EaDQiLi4OISEhl7zOqfBrbm6GxWJBVFSUw/moqCicPHly2Ne1tbUhNjYWJpMJMpkML7/8Mm699Vb797Ozs/HLX/4SU6ZMwXfffYcnn3wSCxYsQEVFBWQy2UX3KyoqwjPPPHPR+dDQUIYfERFdtgvM6T6/0QgJCcHhw4fR0dGBsrIy5OXlISEhAfPmzQMALFmyxH7trFmzkJycjKlTp6K8vBy33HLLRffLz89HXl6e/Wtb0hMREY2EU+Gn0Wggk8lgMBgczhsMBmi12mFfJ5VKMW3aNACATqfDiRMnUFRUZA+/H0tISIBGo8GZM2eGDD+FQuHUgBgiIqILOTXVQS6XIzU1FWVlZfZzVqsVZWVlyMjIGPF9rFarQ5/dj9XV1aGlpQXR0dHOlEdERDQiTj/2zMvLw7Jly5CWlobZs2ejuLgYnZ2dyM3NBQAsXboUsbGxKCoqAtDfP5eWloapU6fCZDJh79692L17N7Zu3QoA6OjowDPPPIM777wTWq0W3333HVavXo1p06Y5TIUgIiJyFafDLycnB01NTSgoKIBer4dOp0Npaal9EExNTY3D3IrOzk48+OCDqKurQ1BQEBITE/HXv/4VOTk5AACZTIajR49i586daG1tRUxMDObPn49nn32WjzaJiMgtJIJtnoAXMxqNUKvVaGtr42hPIiI/NtI8GJPRnkR0aec6zTjeYBS7DPIAKoUMKRPDIJNytSp3YvgRiUwQBPzi5c/wQ0uX2KWQh1h7RxKW3zhF7DJ8GsOPSGQtnWZ78CVqL70qBfk2Y3cvGtp68MnpZoafmzH8iERW1dwJAIgNC0LpIz8RuRoS06Ga8/jly5/jcG0rBEHgQv1uxP38iERmC78pmmCRKyGxJUWHIlAmQUunGXXnu8Uux6cx/IhEZgu/eI1K5EpIbMpAGWZE949QPFLXKm4xPo7hRySyanvLb5zIlZAnSJkYBgA4Utsqah2+juFHJLLBx55s+RGQEhcGADhS2yZuIT6O4UckIqtVQHULW340SBenBgAcq29Dn8UqcjW+i+FHJCJDew96eq2QSSWYGB4kdjnkARI04xCiCEB3rwXfGjrELsdnMfyIRFTV1N/qmxShQqCM/xwJkEolSB5o/XHQi/vwXxuRiKoGHnnGj2d/Hw3ioBf3Y/gRiajaPs2Bc/xokG3Qy2GGn9sw/IhEZBvpmcDwowvoBsLvW0M7usx94hbjoxh+RCKqYsuPhhAVqoQ2VAmrAByr45QHd2D4EYmkz2JFzbn+Ba25tBn9WAoHvbgVw49IJA2tPei1CJAHSBGj5jQHcqSLCwfAye7uwvAjEoltpOfkCBWk3LiUfsTW8uOgF/dg+BGJpKqpfwIzH3nSUGbFqiGRAPWt3WhqN4ldjs9h+BGJpLqF/X00vBBlIKZN6F/yjvP9XG9U4bdlyxbEx8dDqVQiPT0dBw8eHPbad955B2lpaQgLC0NwcDB0Oh12797tcI0gCCgoKEB0dDSCgoKQmZmJ06dPj6Y0Iq/xPffxo8uwL3LNQS8u53T47dmzB3l5eSgsLMShQ4eQkpKCrKwsNDY2Dnl9REQEnnrqKVRUVODo0aPIzc1Fbm4u3n//ffs1GzduxAsvvICSkhIcOHAAwcHByMrKQk9Pz+jfGZGH4wR3uhwdJ7u7jdPht3nzZqxYsQK5ublISkpCSUkJVCoVtm/fPuT18+bNwy9+8QvMmDEDU6dOxcMPP4zk5GR8+umnAPpbfcXFxVizZg0WLVqE5ORk7Nq1Cw0NDXj33Xev6M0ReSpznxV15/nYky5NZ9/eqBWCIIhbjI9xKvzMZjMqKyuRmZk5eAOpFJmZmaioqLjs6wVBQFlZGU6dOoWf/OQnAICqqiro9XqHe6rVaqSnpw97T5PJBKPR6HAQeZOac12wCoBKLkNkiELscshDXa0NgTxACmNPn72PmFzDqfBrbm6GxWJBVFSUw/moqCjo9fphX9fW1oZx48ZBLpdj4cKFePHFF3HrrbcCgP11ztyzqKgIarXafsTFxTnzNohEZ3/kOT4YEgmnOdDQAmVSzIwJBQAcrj0vcjW+ZUxGe4aEhODw4cP44osv8NxzzyEvLw/l5eWjvl9+fj7a2trsR21treuKJRoD9t3bJ/CRJ10ad3Z3jwBnLtZoNJDJZDAYDA7nDQYDtFrtsK+TSqWYNm0aAECn0+HEiRMoKirCvHnz7K8zGAyIjo52uKdOpxvyfgqFAgoFHxWR97JNcJ8ynuFHl8ZBL+7hVMtPLpcjNTUVZWVl9nNWqxVlZWXIyMgY8X2sVitMpv5Jm1OmTIFWq3W4p9FoxIEDB5y6J5E3sW1iy5GedDm28DveYIS5zypuMT7EqZYfAOTl5WHZsmVIS0vD7NmzUVxcjM7OTuTm5gIAli5ditjYWBQVFQHo759LS0vD1KlTYTKZsHfvXuzevRtbt24FAEgkEjzyyCNYt24dpk+fjilTpuDpp59GTEwMFi9e7Lp3SuRBqls4x49GZlKECmGqQLR29eKk3ojkgY1u6co4HX45OTloampCQUEB9Ho9dDodSktL7QNWampqIJUONig7Ozvx4IMPoq6uDkFBQUhMTMRf//pX5OTk2K9ZvXo1Ojs78cADD6C1tRVz585FaWkplEqlC94ikWfpNltwtq1/DivDjy5HIpEgZWIYPv62CUdqWxl+LiIRfGDyiNFohFqtRltbG0JDQ8Uuh+iSTpw1YsGfP4E6KBCHC27laE+6rM37vsULZafxy+tisflXOrHL8WgjzQOu7Uk0xi5c2YXBRyOhs+3tx0EvLsPwIxpjtjU9E/jIk0YoZeBR53dNnTD29IpbjI9g+BGNsQsnuBONxPhxCsRF9G94fKyO8/1cgeFHNMZsIz3jNSqRKyFvYmv9cb6fazD8iMZYlf2x5ziRKyFvwsnursXwIxpDxp5eNHeYAbDlR85JuSD8fGCQvugYfkRjyNbfpxmnQIgyUORqyJvMjFFDJpWgqd0EvZF7nV4phh/RGLIvaM1WHzkpSC7D1VEhADjlwRUYfkRjqLq5f082jvSk0Rh89MkRn1eK4Uc0hqqaOwBwKyMaHdtkd+7td+UYfkRjqGpgN25uZUSjYWv5Hatrg8XKQS9XguFHNEYEQUBVE1t+NHrTI0OgksvQabbgu4H/lmh0GH5EY+R8Vy+MPX0AgMkRDD9ynkwqwaxY26PPVnGL8XIMP6IxYhvpGa1WIkguE7ka8la2ye4c8XllGH5EY2RwmgNbfTR6KVzpxSUYfkRj5MKtjIhGyxZ+J/Xt6Om1iFuMF2P4EY2RKm5lRC4Qo1ZCM04Bi1XANw2c7zdaDD+iMVLFrYzIBSQSyQWLXDP8RovhRzQGBEG4YCsjhh9dGe7sfuUYfkRjoLHdhC6zBVIJMCmC63rSlbH1+x2paxW1Dm/G8CMaA7ZHnhPDVZAH8J8dXZnk2DAAwA8tXTjXaRa3GC81qn+FW7ZsQXx8PJRKJdLT03Hw4MFhr922bRtuuukmhIeHIzw8HJmZmRddv3z5ckgkEocjOzt7NKUReSROcyBXUqsC7QOn2PobHafDb8+ePcjLy0NhYSEOHTqElJQUZGVlobGxccjry8vLcc899+Cjjz5CRUUF4uLiMH/+fNTX1ztcl52djbNnz9qPN998c3TviMgDVTP8yMU42f3KOB1+mzdvxooVK5Cbm4ukpCSUlJRApVJh+/btQ17/+uuv48EHH4ROp0NiYiJeffVVWK1WlJWVOVynUCig1WrtR3h4+LA1mEwmGI1Gh4PIk33P8CMXS2H4XRGnws9sNqOyshKZmZmDN5BKkZmZiYqKihHdo6urC729vYiIiHA4X15ejsjISFx99dVYuXIlWlpahr1HUVER1Gq1/YiLi3PmbRCNOU5wJ1cbHPTSBkHgDg/Ocir8mpubYbFYEBUV5XA+KioKer1+RPd4/PHHERMT4xCg2dnZ2LVrF8rKyrBhwwZ8/PHHWLBgASyWoVcvyM/PR1tbm/2ora115m0QjSmLVcAP57iVEbnWjOgQBMokONdpRu25brHL8ToBY/nD1q9fj7feegvl5eVQKpX280uWLLH/fdasWUhOTsbUqVNRXl6OW2655aL7KBQKKBSKMamZ6Eo1tHbD3GdFoEyC2PAgscshH6EIkCEpOhRH6tpwuK4Vk8ZzCo0znGr5aTQayGQyGAwGh/MGgwFarfaSr920aRPWr1+PDz74AMnJyZe8NiEhARqNBmfOnHGmPCKPZJvcPilCBZlUInI15EvY7zd6ToWfXC5Hamqqw2AV2+CVjIyMYV+3ceNGPPvssygtLUVaWtplf05dXR1aWloQHR3tTHlEHmlwmsM4kSshX8MRn6Pn9GjPvLw8bNu2DTt37sSJEyewcuVKdHZ2Ijc3FwCwdOlS5Ofn26/fsGEDnn76aWzfvh3x8fHQ6/XQ6/Xo6OjfhbijowOPPfYY9u/fj+rqapSVlWHRokWYNm0asrKyXPQ2icQzGH58LEWuZWv5fd3Qhl6LVdxivIzTfX45OTloampCQUEB9Ho9dDodSktL7YNgampqIJUOZurWrVthNptx1113OdynsLAQa9euhUwmw9GjR7Fz5060trYiJiYG8+fPx7PPPst+PfIJHOlJ7jJlfDBClAFo7+nDKX07Zg7s8k6XN6oBL6tWrcKqVauG/F55ebnD19XV1Ze8V1BQEN5///3RlEHkFbi6C7mLVCpBysQwfHqmGUfqWhl+TuAig0Ru1GuxovZ8/zB0hh+5Qwp3eBgVhh+RG9We64LFKiAoUIaoEOXlX0DkJF1c/2pYR7i3n1MYfkRuZJvmMHm8ClJOcyA3SJnY3/L7trEdHaY+kavxHgw/Ijeqah5Y2YWPPMlNIkOViFErIQjA1/Vs/Y0Uw4/Ijaqa+6f0MPzInWxTHg6z32/EGH5EblQ90PLjNAdyJ6704jyGH5Eb2aY5JDD8yI240ovzGH5EbtLTa0FDW/80B7b8yJ1mxaohlQANbT1oNPaIXY5XYPgRuUnNuS4IAhCiDMD4YLnY5ZAPC1YEYHpkCID+/f3o8hh+RG7yfdPgyi4SCac5kHvZJrsfrj0vciXegeFH5Ca2OX7x3MCWxsDgoBe2/EaC4UfkJlVNXNOTxk7KxDAAwJG6VlitgrjFeAGGH5GbVLUw/GjsXK0NgTJQivaePvt/ezQ8hh+Rm3A3BxpLgTIpZsZwkeuRYvgRuUGHqQ9N7SYAnOZAY4crvYwcw4/IDWwb2EYEy6EOChS5GvIXXOll5Bh+RG7AR54kBt3AoJfjZ40w9VnELcbDMfyI3MDW8uM0BxpLcRFBiAiWo9ci4MTZdrHL8WgMPyI3sK/pOYHhR2NHIpHY9/fjo89LY/gRuUEVJ7iTSDjoZWRGFX5btmxBfHw8lEol0tPTcfDgwWGv3bZtG2666SaEh4cjPDwcmZmZF10vCAIKCgoQHR2NoKAgZGZm4vTp06Mpjcgj2B97alQiV0L+hoNeRsbp8NuzZw/y8vJQWFiIQ4cOISUlBVlZWWhsbBzy+vLyctxzzz346KOPUFFRgbi4OMyfPx/19fX2azZu3IgXXngBJSUlOHDgAIKDg5GVlYWeHq5OTt6ntcuM8129ANjyo7FnW+nl++ZOtA38d0gXczr8Nm/ejBUrViA3NxdJSUkoKSmBSqXC9u3bh7z+9ddfx4MPPgidTofExES8+uqrsFqtKCsrA9Df6isuLsaaNWuwaNEiJCcnY9euXWhoaMC77757RW+OSAy2/r6oUAWCFQEiV0P+JiJYjkkR/U8cjta3iluMB3Mq/MxmMyorK5GZmTl4A6kUmZmZqKioGNE9urq60Nvbi4iICABAVVUV9Hq9wz3VajXS09OHvafJZILRaHQ4iDwFpzmQ2Li57eU5FX7Nzc2wWCyIiopyOB8VFQW9Xj+iezz++OOIiYmxh53tdc7cs6ioCGq12n7ExcU58zaI3Kqa4UciGxz0wh0ehjOmoz3Xr1+Pt956C//85z+hVCpHfZ/8/Hy0tbXZj9raWhdWSXRlqlq6ALC/j8Sjs+/t1wpB4A4PQ3Eq/DQaDWQyGQwGg8N5g8EArVZ7yddu2rQJ69evxwcffIDk5GT7edvrnLmnQqFAaGiow0HkKaqaOwCw5UfiuSZGDZlUguYOExraOHBwKE6Fn1wuR2pqqn2wCgD74JWMjIxhX7dx40Y8++yzKC0tRVpamsP3pkyZAq1W63BPo9GIAwcOXPKeRJ5IEARUN/e3/Bh+JBZloAyJ2hAA7PcbjtOPPfPy8rBt2zbs3LkTJ06cwMqVK9HZ2Ync3FwAwNKlS5Gfn2+/fsOGDXj66aexfft2xMfHQ6/XQ6/Xo6Oj//+OJRIJHnnkEaxbtw7//d//jWPHjmHp0qWIiYnB4sWLXfMuicZIU4cJHaY+SCTApPGc40fi4aCXS3N6HHZOTg6amppQUFAAvV4PnU6H0tJS+4CVmpoaSKWDmbp161aYzWbcddddDvcpLCzE2rVrAQCrV69GZ2cnHnjgAbS2tmLu3LkoLS29on5BIjHYWn2xYUFQBMhErob8WUpcGF4/UMOVXoYhEXygN9RoNEKtVqOtrY39fySqv31Ri9X/OIqbpmuw+/50scshP/atoR3z//Q/UMllOLY2CzKpROySxsRI84BrexK50Pec5kAeYuqEcQiWy9BltuB0I3d4+DGGH5ELcSsj8hQyqQSzuMPDsBh+RC5kX92FWxmRB9DFhQPgZPehMPyIXMRqFVA9sJXRFLb8yAPYJruz5Xcxhh+Ri+iNPTD1WREglWBieJDY5RDZlzk7ZWhHt9kibjEehuFH5CK2R56TIlQIkPGfFolPG6pEZIgCFquArxv46PNC/BdK5CJV9g1s+ciTPINEIuHmtsNg+BG5CLcyIk+ks+/w0CpqHZ6G4UfkItVs+ZEHsi9zVtcqah2ehuFH5CK2ll8Cw488iG2uX+25brR0mESuxnMw/IhcoM9iRc25gX38GH7kQUKVgZg6MO+Urb9BDD8iF6hv7UafVYAiQIroUC7ITp6FO7tfjOFH5ALfX7CsmdRPFhAm78HtjS7G8CNygcHBLtzDjzzPhYNefGAjH5dg+BG5wOA0h3EiV0J0sURtKOQyKVq7eu190/6O4UfkAoPhx5YfeR55gBRJMf1723G+Xz+GH5EL2Ba05lZG5Kk42d0Rw4/oCpn6LKg/3w2AWxmR50rhDg8OGH5EV6j2XBesAhAsl2HCOIXY5RANyba339cNRvRarCJXIz6GH9EV+r5pcANbiYTTHMgzxY9XIVQZAHOfFaf07WKXI7pRhd+WLVsQHx8PpVKJ9PR0HDx4cNhrv/nmG9x5552Ij4+HRCJBcXHxRdesXbsWEonE4UhMTBxNaURjjv195A0u3OHhKz76dD789uzZg7y8PBQWFuLQoUNISUlBVlYWGhsbh7y+q6sLCQkJWL9+PbRa7bD3veaaa3D27Fn78emnnzpbGpEoqpr7h45zTU/ydJzsPsjp8Nu8eTNWrFiB3NxcJCUloaSkBCqVCtu3bx/y+uuvvx7PP/88lixZAoVi+P6QgIAAaLVa+6HRaJwtjUgUVc0dALimJ3m+lIlhABh+gJPhZzabUVlZiczMzMEbSKXIzMxERUXFFRVy+vRpxMTEICEhAffddx9qamqGvdZkMsFoNDocRGKpbuaC1uQdbI89zzR1oL2nV9xiROZU+DU3N8NisSAqKsrhfFRUFPR6/aiLSE9Px2uvvYbS0lJs3boVVVVVuOmmm9DePnSnbFFREdRqtf2Ii4sb9c8muhJd5j7ojT0A+NiTPN+EEAViw4IgCMCxev9e5NojRnsuWLAAd999N5KTk5GVlYW9e/eitbUVf/vb34a8Pj8/H21tbfajtrZ2jCsm6mdr9YWpAhGmkotcDdHlDfb7+Xf4BThzsUajgUwmg8FgcDhvMBguOZjFWWFhYbjqqqtw5syZIb+vUCgu2X9INFZsIz2nsNVHXiIlTo1/HTuLw7XnxS5FVE61/ORyOVJTU1FWVmY/Z7VaUVZWhoyMDJcV1dHRge+++w7R0dEuuyeRO9jX9OQ0B/ISg4Ne2PJzSl5eHpYtW4a0tDTMnj0bxcXF6OzsRG5uLgBg6dKliI2NRVFREYD+QTLHjx+3/72+vh6HDx/GuHHjMG3aNADAH/7wB9xxxx2YPHkyGhoaUFhYCJlMhnvuucdV75PILarsWxkx/Mg7zIxVQyoB9MYe6Nt6oFX75+bLTodfTk4OmpqaUFBQAL1eD51Oh9LSUvsgmJqaGkilgw3KhoYGXHvttfavN23ahE2bNuHmm29GeXk5AKCurg733HMPWlpaMGHCBMydOxf79+/HhAkTrvDtEbnX4G4ODD/yDsGKAFwVFYKT+nYcqWuFVu26Litv4nT4AcCqVauwatWqIb9nCzSb+Pj4y26e+NZbb42mDCLRVTP8yAvp4sL6w6+2FVnX+Gf4ecRoTyJv1Nbdi5ZOMwA+9iTvksLtjRh+RKNla/VNCFFgnGJUD1GIRGEb9HK0rg1W66WfzPkqhh/RKNmnOXCkJ3mZq6LGQRkoRYepD98PLM/nbxh+RKNk38qIjzzJywTIpJgV27+57WE/nfLA8CMaJftWRgw/8kL+vsMDw49olDjNgbyZvw96YfgRjYIgCAw/8mq2QS8nzhrR02sRtxgRMPyIRuFcpxntPX2QSIDJ41Vil0PktInhQRgfLEefVcDxs/63LRzDj2gUbK2+GHUQlIEykashcp5EIrE/+vTHfj+GH9EoDK7pyVYfeS9/HvTC8CMaBfb3kS/w50EvDD+iUbBPc+AEd/JiKRP75/pVt3ShtcsscjVji+FHNApVAzu4J0xg+JH3ClPJET8wYOtInX9Ndmf4ETlJEAT7up5s+ZG389dBLww/IicZjCZ091ogk0oQF8EBL+Td/HXQC8OPyEm2hYDjwoMQKOM/IfJu9pZfXetl9171JfyXS+Sk6oH+Pq7pSb4gKToUAVIJmjvMqDvfLXY5Y4bhR+Qk+1ZGDD/yAcpAGWZEhwLob/35C4YfkZO4lRH5mpS4/ikP/tTvx/AjchLn+JGvsS1yfcSP9vZj+BE5wWIVUNPS3+fHlh/5imsnhQEAjtW3oc9iFbeYMTKq8NuyZQvi4+OhVCqRnp6OgwcPDnvtN998gzvvvBPx8fGQSCQoLi6+4nsSiaWhtRtmixVymRQxYUFil0PkEgmacRinCEB3rwXfGjrELmdMOB1+e/bsQV5eHgoLC3Ho0CGkpKQgKysLjY2NQ17f1dWFhIQErF+/Hlqt1iX3JBKLbU3PyeNVkEklIldD5BpSqQTJA0ud+cugF6fDb/PmzVixYgVyc3ORlJSEkpISqFQqbN++fcjrr7/+ejz//PNYsmQJFAqFS+5pMplgNBodDqKxMLibAx95km/xt5VenAo/s9mMyspKZGZmDt5AKkVmZiYqKipGVcBo7llUVAS1Wm0/4uLiRvWziZzF3RzIV9kGvfjLDg9OhV9zczMsFguioqIczkdFRUGv14+qgNHcMz8/H21tbfajtrZ2VD+byFkMP/JVtkEv3xra0WXuE7eYMeCVoz0VCgVCQ0MdDqKxwGkO5KuiQpXQhiphFYBjfrDDg1Php9FoIJPJYDAYHM4bDIZhB7OIcU8idzD3We3LP3ErI/JF9snufjDoxanwk8vlSE1NRVlZmf2c1WpFWVkZMjIyRlWAO+5J5A6157tgsQpQyWWIDBl68BaRNxsc9OL7Lb8AZ1+Ql5eHZcuWIS0tDbNnz0ZxcTE6OzuRm5sLAFi6dCliY2NRVFQEoH9Ay/Hjx+1/r6+vx+HDhzFu3DhMmzZtRPck8gQX7uEnkXCaA/kenR8NenE6/HJyctDU1ISCggLo9XrodDqUlpbaB6zU1NRAKh1sUDY0NODaa6+1f71p0yZs2rQJN998M8rLy0d0TyJPwMEu5OtmTVRDIgHqW7vR1G7CBB9+wiERfGADJ6PRCLVajba2Ng5+Ibd54h9H8dYXtXjop1PxWFai2OUQucWtmz/G6cYO/N9labhlhvc1QEaaB1452pNorB2ra8M7h+oBAMkDj4aIfJGt38/XH30y/Iguw9jTi4feOASzxYr5SVGYn+R9/zdMNFIMPyKCIAh4/O9HUXOuCxPDg/D8XSkc7EI+TWff3qgVPtArNiyGH9El7Kr4Af/+Wo9AmQQv3Xsd1KpAsUsicqurtSGQB0hh7OlD9cD2Xb6I4Uc0jKN1rVj3r/5pOvkLZkA38DiIyJfJA6SYGdM/UMSXF7lm+BENoa27v5+v1yIg65oo5N4YL3ZJRGPGH/r9GH5EPyIIAlb//Qhqz3UjLiIIG9nPR35Gx/Aj8j+vfV6N978xIFAmwZZ7r4M6iP185F9s2xsdbzDC3GcVtxg3YfgRXeBIbSv+z94TAICnbpvBOX3klyaPVyFMFQizxYqTet/cLJzhRzSgrWuwn2/BTC2WzYkXuyQiUUgkEnvrz1cHvTD8iNDfz/fY34+g7nw3JkWosOGuZPbzkV+zDXr5iuFH5Lt2fFaND44bIJdJseXe6xCqZD8f+TedbW8/hh+Rbzpc24qif/f38625fQZmTVSLXBGR+Gz93d81dcLY0ytuMW7A8CO/1tbVi4de7+/nu22WFr++YbLYJRF5BM04BSaGBwHoX9jd1zD8yG8JgoA//P0I6lu7MXm8CuvvZD8f0YV8ebI7w4/81v/9tAr72M9HNKxrGX5EvuWrmvNY/++TAICnb5+BmbHs5yP6sQtbfr62wwPDj/xOa5cZq974Cn1WAQtnReO/2M9HNKRrYkIhk0rQ1G6C3tgjdjkuxfAjvyIIAv7w9oX9fLPYz0c0DJU8AFdFhQDwvSkPDD/yK69+UoX/nGiEPKC/ny+E/XxEl2Sb73e41rdGfI4q/LZs2YL4+HgolUqkp6fj4MGDl7z+7bffRmJiIpRKJWbNmoW9e/c6fH/58uWQSCQOR3Z29mhKIxpW5Q/nsaG0v5+v4PYk9vMRjYBthwe/b/nt2bMHeXl5KCwsxKFDh5CSkoKsrCw0NjYOef3nn3+Oe+65B/fffz+++uorLF68GIsXL8bXX3/tcF12djbOnj1rP958883RvSOiIbR2mfG7Nw6hzyrg9uRo3Jc+SeySiLyCbdDL0bpWWKy+M+jF6fDbvHkzVqxYgdzcXCQlJaGkpAQqlQrbt28f8vo///nPyM7OxmOPPYYZM2bg2WefxXXXXYeXXnrJ4TqFQgGtVms/wsPDR/eOiH7EahXw6N+OoKGtB1M0wSj6Jfv5iEZqemQIVHIZOs0WfNfUIXY5LuNU+JnNZlRWViIzM3PwBlIpMjMzUVFRMeRrKioqHK4HgKysrIuuLy8vR2RkJK6++mqsXLkSLS0tw9ZhMplgNBodDqLhvPrp9yg72d/P99K917Kfj8gJMqnE3kXgS/P9nAq/5uZmWCwWREVFOZyPioqCXq8f8jV6vf6y12dnZ2PXrl0oKyvDhg0b8PHHH2PBggWwWCxD3rOoqAhqtdp+xMXFOfM2yI9U/nAOG0pPAQAK70jCNTHs5yNyli/2+wWIXQAALFmyxP73WbNmITk5GVOnTkV5eTluueWWi67Pz89HXl6e/Wuj0cgApIuc7zTjd298BYtVwB0pMbh3Nvv5iEbDHn51raLW4UpOtfw0Gg1kMhkMBoPDeYPBAK1WO+RrtFqtU9cDQEJCAjQaDc6cOTPk9xUKBUJDQx0OogtZrQIefZv9fESuYBv0cvJsO04b2sUtxkWcCj+5XI7U1FSUlZXZz1mtVpSVlSEjI2PI12RkZDhcDwD79u0b9noAqKurQ0tLC6Kjo50pj8julU++x4cnG6EYmM83TuERDzmIvFKMWomMhPHoswpYvuMLNLZ7/2ovTo/2zMvLw7Zt27Bz506cOHECK1euRGdnJ3JzcwEAS5cuRX5+vv36hx9+GKWlpfjjH/+IkydPYu3atfjyyy+xatUqAEBHRwcee+wx7N+/H9XV1SgrK8OiRYswbdo0ZGVluehtkj/5svocnn+/v59v7c+vQVIMnwwQXQmJRIIt912HKZpg1Ld243/v/BJd5j6xy7oiTodfTk4ONm3ahIKCAuh0Ohw+fBilpaX2QS01NTU4e/as/fo5c+bgjTfewCuvvIKUlBT8/e9/x7vvvouZM2cCAGQyGY4ePYqf//znuOqqq3D//fcjNTUVn3zyCRQKhYveJvmL75o68Ls3+/v5FulisOR69gUTuUJEsBw7ll+PcFUgjta14fdvHvbqeX8SwQeW6jYajVCr1Whra2P/n59q7jDhz/85jTcO1sBiFZCgCcZ//24uH3cSuVjlD+dwz7YDMPdZsXxOPArvSPKo/vSR5gHX9iSv1m22YMtHZzDv+XLs3v8DLFYBmTMisfN/zWbwEblB6uQI/OlXOgDAa59XY/tn1aLWM1r87UBeyWIV8M6hOvzxg2/tW63MilXjydtmIGPqeJGrI/JtC5OjUXc+EUX/Pol1/zqO2LAgZM8cfgS/J2L4kdf55HQT/s/ekzhxtn9ln9iwIKzOvhp3JMdAKvWcxy9EvuyBnySg5lwXXj9Qg0f2fIU3Q2/AtZO8Z1lKhh95jZN6I4r2nsTH3zYBAEKUAVj102lYNiceykCZyNUR+ReJRIJnfn4N6lu7UX6qCf9755f454M3YtJ4ldiljQgHvJDHMxh78McPTuHvlXWwCkCgTIL/umEyfv+z6QgPlotdHpFf6zD14VclFTh+1oiECcF4Z+UchKnE+3c50jxg+JHH6jD14ZWPv8O2T6rQ3du/zutts7RYnZWIeE2wyNURkY3B2IPFWz7D2bYepE+JwK77Z0MRIM7TGIYfea0+ixV7vqzFn/adRnOHCQCQOjkcT942A6mTvadPgcifnNQbcdfWCnSY+rBYF4M/5ehEmQIx0jxgnx95DEEQ8OHJRhT9+yTONPbvGxY/XoXHsxORPVPrUXOJiMhRojYUW//rOuTu+ALvHm7ApAgV8uZfLXZZw2L4kUc4VteG5/Yex/7vzwEAwlWB+P0t03Ff+mTIAzgdlcgb3DR9Ap77xUw8/o9jeOHDM5gYocKv0jxzlSWGH4mq7nwXNr1/Cu8ebgAAyAOkyL0xHg/OmwZ1EDedJfI2OddPQu25brz00Rk8+c4xxKiDMHe6RuyyLsLwI1G0dffi5fIz2PFZNcx9VgDAL66NxaPzr8LEcO8YKk1EQ3t0/lWoPd+F/3e4ASv/Wom3V2YgUetZ4zEYfjQmrFYBx88a8emZZnx2phkHq87BNBB6NyRE4KnbkjBrIndZJ/IFEokEG+9KxtnWHhysPof/teML/POhGxEVqhS7NDuO9iS3qW/txqenm/DJ6WZ8/l0LznWaHb4/PXIcnliQiJ8lRnIwC5EPau0y45dbP8f3TZ24JiYUf/tNBoLdvOYupzrQmGvr7sX+71vw6elmfHqmGVXNnQ7fD5bLcEPCeNw4TYObpmswLXIcQ4/Ix9W0dOEXL3+Glk4zfpYYiVd+nYoAmfsGsTH8yO3MfVYcrm3tb92dacaR2lZcuL2XTCpBykQ15k6fgLnTNNDFhXHkJpEfOlRzHve8sh+mPiv+64ZJeHbRTLf9jy/n+ZHLCYKA040d9pbd/u9b0GW2OFyToAnG3OkazJ2mwQ1TxyNUyRGbRP7uuknhKM7R4cE3DuGv+2swOSIYK36SIGpNDD+6pEZjDz4902wfqGIwmhy+Pz5YjjnTNLhpmgY3TtcgNixIpEqJyJMtmBWNp26bgXX/OoHn9p5AbHgQbpsVLVo9DD8/JwgCmtpNqDnX5XDUDvz547BTBEgxe0oEbpquwY3TNJihDeU2QkQ0IvfPnYLac13YWfED/r89hxEVqhRtyUKGnx/oNltQd37ocKs514WeXuuwr5VIgJkxasyd3t+6u25yOLcPIqJRkUgkKLijfxuk/5xoxIpdX+KdlXNEWaieA158gNUqoKljoPXWcnG4NbabLvl6qQSICQvCpAgVJkWoEDfw56QIFeI1wVxphYhcqsvch5y/7Mex+jZMHq/COyvnYPw4hUvu7dbRnlu2bMHzzz8PvV6PlJQUvPjii5g9e/aw17/99tt4+umnUV1djenTp2PDhg247bbb7N8XBAGFhYXYtm0bWltbceONN2Lr1q2YPn36iOrxhfDrs1jR3tOHtu5eGHt6YezuG/iz/+u27h+fG7i2uxet3b32VVKGE6IIwKTxg6F24d9jwoIQ6Mahx0REP9bUbsIvt36G2nPdSIkLw5sr0qGSX/nDSLeF3549e7B06VKUlJQgPT0dxcXFePvtt3Hq1ClERkZedP3nn3+On/zkJygqKsLtt9+ON954Axs2bMChQ4cwc+ZMAMCGDRtQVFSEnTt3YsqUKXj66adx7NgxHD9+HErl5VcEcHX4CYIAi1VAr0WAuc8Ks6X/6O2zotdihWngT3OfFb0WYYhzF14n2M91mPrsYWbsvjDoetH5o1GTzpJJJYgJUw7ZepsUoYI6KJBz6ojIo3zX1IG7tn6O8129LpsD6LbwS09Px/XXX4+XXnoJAGC1WhEXF4ff/e53eOKJJy66PicnB52dnXjvvffs52644QbodDqUlJRAEATExMTg0UcfxR/+8AcAQFtbG6KiovDaa69hyZIlLnuzl5P67D50mPpgtlgh1sPgYLkMoUGBCFUGQh0UiNCgAIQqAwfOBfT/OfB92/fUQYHQqpVsvRGR16n84Tzu3dY/B/BXaROx4c7kK/ofdbfM8zObzaisrER+fr79nFQqRWZmJioqKoZ8TUVFBfLy8hzOZWVl4d133wUAVFVVQa/XIzMz0/59tVqN9PR0VFRUDBl+JpMJJtNgP5bRaHTmbQyrp9diX2/yx+QyKeQBUgTKJAgc+PvgOcfvKRzODV4brJDZg0w9RICFKAPcuvIBEZGnSZ0cjpfuvQ6/2f0l/vZlHSaGq/D7W0bW5XUlnAq/5uZmWCwWREVFOZyPiorCyZMnh3yNXq8f8nq9Xm//vu3ccNf8WFFREZ555hlnSh+Rf/3+JsikkovCK1Am4SNDIiI3uTUpCusWz8Lmfacw7+oJY/IzvXKqQ35+vkNr0mg0Ii7uyjdMFGO4LRERAfemT8LCWdFQq8ZmdLlTz9g0Gg1kMhkMBoPDeYPBAK1WO+RrtFrtJa+3/enMPRUKBUJDQx0OIiLybmMVfICT4SeXy5GamoqysjL7OavVirKyMmRkZAz5moyMDIfrAWDfvn3266dMmQKtVutwjdFoxIEDB4a9JxER0ZVw+rFnXl4eli1bhrS0NMyePRvFxcXo7OxEbm4uAGDp0qWIjY1FUVERAODhhx/GzTffjD/+8Y9YuHAh3nrrLXz55Zd45ZVXAPTP+H/kkUewbt06TJ8+3T7VISYmBosXL3bdOyUiIhrgdPjl5OSgqakJBQUF0Ov10Ol0KC0ttQ9YqampgVQ62KCcM2cO3njjDaxZswZPPvkkpk+fjnfffdc+xw8AVq9ejc7OTjzwwANobW3F3LlzUVpaOqI5fkRERM7i8mZEROQzRpoHnFRGRER+h+FHRER+xyvn+f2Y7cmtq1Z6ISIi72TLgcv16PlE+LW3twOASya6ExGR92tvb4darR72+z4x4MVqtaKhoQEhISGQSCT2FV9qa2s5AGYY/IxGhp/T5fEzujx+RiPjis9JEAS0t7cjJibGYebBj/lEy08qlWLixIkXnefqL5fHz2hk+DldHj+jy+NnNDJX+jldqsVnwwEvRETkdxh+RETkd3wy/BQKBQoLC6FQKMQuxWPxMxoZfk6Xx8/o8vgZjcxYfk4+MeCFiIjIGT7Z8iMiIroUhh8REfkdhh8REfkdhh8REfkdhh8REfkdnwu/5557DnPmzIFKpUJYWNiQ19TU1GDhwoVQqVSIjIzEY489hr6+vrEt1MN8++23WLRoETQaDUJDQzF37lx89NFHYpflcf71r38hPT0dQUFBCA8Px+LFi8UuySOZTCbodDpIJBIcPnxY7HI8SnV1Ne6//35MmTIFQUFBmDp1KgoLC2E2m8UuTVRbtmxBfHw8lEol0tPTcfDgQbf+PJ8LP7PZjLvvvhsrV64c8vsWiwULFy6E2WzG559/jp07d+K1115DQUHBGFfqWW6//Xb09fXhww8/RGVlJVJSUnD77bdDr9eLXZrH+Mc//oFf//rXyM3NxZEjR/DZZ5/h3nvvFbssj7R69WrExMSIXYZHOnnyJKxWK/7yl7/gm2++wZ/+9CeUlJTgySefFLs00ezZswd5eXkoLCzEoUOHkJKSgqysLDQ2Nrrvhwo+aseOHYJarb7o/N69ewWpVCro9Xr7ua1btwqhoaGCyWQawwo9R1NTkwBA+J//+R/7OaPRKAAQ9u3bJ2JlnqO3t1eIjY0VXn31VbFL8Xh79+4VEhMThW+++UYAIHz11Vdil+TxNm7cKEyZMkXsMkQze/Zs4aGHHrJ/bbFYhJiYGKGoqMhtP9PnWn6XU1FRgVmzZiEqKsp+LisrC0ajEd98842IlYln/PjxuPrqq7Fr1y50dnair68Pf/nLXxAZGYnU1FSxy/MIhw4dQn19PaRSKa699lpER0djwYIF+Prrr8UuzaMYDAasWLECu3fvhkqlErscr9HW1oaIiAixyxCF2WxGZWUlMjMz7eekUikyMzNRUVHhtp/rd+Gn1+sdgg+A/Wt/fcQnkUjwn//8B1999RVCQkKgVCqxefNmlJaWIjw8XOzyPML3338PAFi7di3WrFmD9957D+Hh4Zg3bx7OnTsncnWeQRAELF++HL/97W+RlpYmdjle48yZM3jxxRfxm9/8RuxSRNHc3AyLxTLk72V3/k72ivB74oknIJFILnmcPHlS7DI9zkg/N0EQ8NBDDyEyMhKffPIJDh48iMWLF+OOO+7A2bNnxX4bbjXSz8hqtQIAnnrqKdx5551ITU3Fjh07IJFI8Pbbb4v8LtxrpJ/Riy++iPb2duTn54tdsihG83uqvr4e2dnZuPvuu7FixQqRKvdPXrGf36OPPorly5df8pqEhIQR3Uur1V40ishgMNi/50tG+rl9+OGHeO+993D+/Hn7Hlovv/wy9u3bh507d+KJJ54Yg2rFMdLPyPY/AUlJSfbzCoUCCQkJqKmpcWeJonPmv6OKioqLFiVOS0vDfffdh507d7qxSvE5+3uqoaEBP/3pTzFnzhy88sorbq7Oc2k0GshkMvvvYRuDweDW38leEX4TJkzAhAkTXHKvjIwMPPfcc2hsbERkZCQAYN++fQgNDXX4xeYLRvq5dXV1AcBFux5LpVJ7i8dXjfQzSk1NhUKhwKlTpzB37lwAQG9vL6qrqzF58mR3lymqkX5GL7zwAtatW2f/uqGhAVlZWdizZw/S09PdWaJHcOb3VH19PX7605/anyBcasdxXyeXy5GamoqysjL71CGr1YqysjKsWrXKbT/XK8LPGTU1NTh37hxqampgsVjsc4ymTZuGcePGYf78+UhKSsKvf/1rbNy4EXq9HmvWrMFDDz3kt9uNZGRkIDw8HMuWLUNBQQGCgoKwbds2VFVVYeHChWKX5xFCQ0Px29/+FoWFhYiLi8PkyZPx/PPPAwDuvvtukavzDJMmTXL4ety4cQCAqVOnYuLEiWKU5JHq6+sxb948TJ48GZs2bUJTU5P9e7729Gmk8vLysGzZMqSlpWH27NkoLi5GZ2cncnNz3fdD3TaOVCTLli0TAFx0fPTRR/ZrqqurhQULFghBQUGCRqMRHn30UaG3t1e8oj3AF198IcyfP1+IiIgQQkJChBtuuEHYu3ev2GV5FLPZLDz66KNCZGSkEBISImRmZgpff/212GV5rKqqKk51GMKOHTuG/B3lg7+OnfLiiy8KkyZNEuRyuTB79mxh//79bv153M+PiIj8jv8+aCYiIr/F8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/D8CMiIr/z/wNUufAJSa9izgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "t = np.logspace(-5, -2, 100)\n",
    "h = ml.head(0, 0, t)\n",
    "zc = 0.5 * (z[:-1] + z[1:])\n",
    "plt.plot(zc, h[:, 10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x15b8b3560>]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAESCAYAAACCU7B8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuy0lEQVR4nO3deVhU9eI/8PeZGWbYB5BNFEQwNTdUEHIjMxNcUsrcMhUzM9Nyu12xW1nfbj8ru5mZqZWCZtesrLyZomgqiAsuuAuKIqDIpjKDss+c3x8IRQKCwpwZ5v16nvMIh8858z6cmefth5kzI4iiKIKIiMhMyKQOQEREZEgsPiIiMissPiIiMissPiIiMissPiIiMissPiIiMissPiIiMisKqQM8LL1ej8zMTNjZ2UEQBKnjEBGRBERRREFBATw8PCCT1T2nM/niy8zMhKenp9QxiIjICGRkZKB169Z1jjH54rOzswNQcbD29vYSpyEiIilotVp4enpWdUJdTL74Kv+8aW9vz+IjIjJz9XnKiy9uISIis8LiIyIis8LiIyIis8LiIyIis2IUxbdixQp4e3vD0tISQUFBSEhIkDoSERE1U5IX36ZNmzBv3jwsWrQIx48fh5+fH0JCQpCTkyN1NCIiaoYEqT+BPSgoCL169cIXX3wBoOKdWDw9PfHaa68hIiLivttrtVqo1WpoNJoHvpwh42Yh0m4UQiYAMpkAuUyATKj8F1Vf//17mSBUjBcEyGSo+PfuOgu5AEuFHDIZ302GiKipNaQLJL2Or7S0FMeOHcPChQur1slkMgwaNAgHDx6scZuSkhKUlJRUfa/Vah86x7bT17F4e9JD76cmSrkMlhYyWFrI7y5//VoOS4Wsar3V3XWqynGKiu+tlH9+rfrLuHu2U8j4tm1ERPchafHl5eVBp9PBzc2t2no3NzckJdVcRIsXL8Z7773XqDkcbZTo6G4HvShCpxehF3H3XxF6vQidKEKnr3gvOF3lmLvr9SKqvq5p7lyq06NUp4e2uLxRM9dGpZDBSim/W5Q1FK5CXvFzCxlUinvLs3KcqnKcQnZPyVYWs1LOoiUi02Ny79yycOFCzJs3r+r7yrepeRhjAjwxJuDh3+9TFP8szXK9HsVlehSV6VBctehRUqa7u05fsa5ch6JSHUrK9dXGFf11m/KKr6ttV/bn+HL9n41bUq5HSbkeQNlDH8/9CAL+LNK7BamqqUgVFeutqs14a575utqr0NGd78BDRE1H0uJzdnaGXC5HdnZ2tfXZ2dlwd3evcRuVSgWVSmWIeA0mCALkAiCXCVBCBmulYW63XKdHcbkeRaUVRVhRlNXLs7I4S+4pXT2K7xZryV/GVW1XXrH+rwVe2bOiCBTdHd+YpvVvi38N69So+yQiqiRp8SmVSvj7+2P37t0ICwsDUPHilt27d2PWrFlSRjMpCrkMtnIZbFVNfzpFUUSZTqwoy9I/y/HemerdIr07o60+rvqM9a/bJWUV4Ou4VLirrTC1X9smPx4iMj+S/6lz3rx5mDx5MgICAhAYGIjPPvsMd+7cwZQpU6SORjUQBAFKhQClQgZ7S4tG3/+qfZfw4fYk/Pv3c/BQW2JI15aNfhtEZN4kL76xY8ciNzcX77zzDrKystC9e3dER0ff84IXMg/Tg31w7VYRvj2UhtmbTsDFToUAbyepYxFRMyL5dXwPqzGu4yPjotOLmP7tMew6nw0HawtsntEHvi62UsciIiPWkC6Q/J1biP5OLhOwfHwP+Hk6IL+wDOGRCcgtKLn/hkRE9cDiI6NkpZRjzeQAeDlZI+NmEaauO4LCUsNcC0lEzRuLj4yWs60K614MhKO1BU5d1WDWfxNRrtNLHYuITByLj4xaW2cbfDO5F1QKGf5IysHbW87CxJ+WJiKJsfjI6Pm3ccSycT0gCMDGhHR8ufeS1JGIyISx+MgkhHZxx6LhFe/msmRHMn5JvCpxIiIyVSw+MhnhfdtiWv+Kd3P550+ncCAlT+JERGSKWHxkUhYOeRTDurVEma7iWr+krIf/WCoiMi8sPjIpMpmA/4z2Q6C3EwpKyhG+9giua4qkjkVEJoTFRybH0kKOryb5w9fFBlnaYkyJPAJtcdN/DBMRNQ8sPjJJDtZKRE0JhIudCklZBZix4RhKy3mNHxHdH4uPTJankzUiw3vBWilHfMoNRGw+xWv8iOi+WHxk0rq0UmPFhJ6QywT8nHgN/9l5QepIRGTkWHxk8p7o4IoPwroAAL7Yk4L/Hk6XOBERGTMWHzUL4wK98PrAdgCAt7ecwR9J2RInIiJjxeKjZmPuU+0xqmdr6PQiZn6XiFNX86WORERGiMVHzYYgCFj8bFf0a+eMojIdXow6goybhVLHIiIjw+KjZkWpkGHlCz3R0d0OebdLMTkyAbfulEodi4iMCIuPmh07SwtETQlES7UlLufewbT1R1FcppM6FhEZCRYfNUvuaktETQmEnaUCR9NuYd4PJ6DX8xo/ImLxUTPWwd0Oqyf6w0IuYNvpLHyw7bzUkYjICLD4qFnr4+uMJc/5AQDW7E/F2v2pEiciIqmx+KjZC+vRCm+EdAAAvP/7OWw/fV3iREQkJRYfmYVXB/hiQpAXRBGYs+kEjqXdlDoSEUmExUdmQRAEvDeiM57s6IqScj2mrjuKS7m3pY5FRBJg8ZHZUMhlWP58D/i1ViO/sAzhkQnILSiROhYRGRiLj8yKtVKBbyb3gqeTFTJuFmHquiMoLC2XOhYRGRCLj8yOi50K66YEwtHaAqeuavDafxNRruOH2BKZCxYfmSUfF1t8MzkAKoUMu5NysOh/Z/khtkRmgsVHZsu/jROWjesOQQC+O5yOlfsuSR2JiAyAxUdmLbRLS7wzvBMA4OPoZPyaeE3iRETU1CQrvitXrmDq1Klo27YtrKys4Ovri0WLFqG0lO+kT4Y1pW9bvNSvLQDgjZ9O4kBKnsSJiKgpSVZ8SUlJ0Ov1WL16Nc6ePYulS5di1apVePPNN6WKRGbszaGPYljXlijTiZj+7TEkZxVIHYmImoggGtEz+kuWLMHKlStx+fLlem+j1WqhVquh0Whgb2/fhOmouSsu02HimsM4cuUWWqot8curfeGutpQ6FhHVQ0O6wKie49NoNHBycqpzTElJCbRabbWFqDFYWsjx9aQA+LrY4LqmGOGRCSgoLpM6FhE1MqMpvpSUFCxfvhzTp0+vc9zixYuhVqurFk9PTwMlJHPgYK1E1JRAONuqkJRVgBkbjqO0nNf4ETUnjV58EREREAShziUpKanaNteuXUNoaChGjx6NadOm1bn/hQsXQqPRVC0ZGRmNfQhk5jydrBEZ3gvWSjn2p+Qh4udTvMaPqBlp9Of4cnNzcePGjTrH+Pj4QKlUAgAyMzMxYMAAPPbYY4iKioJM1rAu5nN81FT2JOXgpfVHodOLeH1gO8wb3EHqSERUi4Z0gaQvbrl27RqeeOIJ+Pv7Y8OGDZDL5Q3eB4uPmtLGhHQs/Pk0AODDZ7tiXKCXxImIqCYm8eKWa9euYcCAAfDy8sInn3yC3NxcZGVlISsrS6pIRPcYH+iF1wa2AwD869cz2JOcI3EiInpYCqluOCYmBikpKUhJSUHr1q2r/YzPp5AxmfdUe1zLL8LPx69h5nfH8cP03ujSSi11LCJ6QJLN+MLDwyGKYo0LkTERBAEfPtsN/do5o7BUh/DII8i4WSh1LCJ6QEZzOQORMVMqZPjyhZ7o6G6HvNslCI9MQH4h316PyBSx+Ijqyd7SAlFTAtFSbYlLuXfw8vpjKC7TSR2LiBqIxUfUAO5qS0RO6QU7lQIJV25i/o8nodfzz/NEpoTFR9RAHd3tsXqiPyzkAn4/dR2Lt5+XOhIRNQCLj+gB9GnnjCXP+QEAvo5LRVR8qsSJiKi+WHxEDyisRyu8EVLxbi7vbT2H6DO8BpXIFLD4iB7CqwN88XyQF0QRmP19Io6l3ZI6EhHdB4uP6CEIgoD/G9EZT3Z0RUm5Hi+tO4LLubeljkVEdWDxET0khVyG5c/3QLfWatwqLEN45BHk3S6ROhYR1YLFR9QIrJUKrJncC55OVki/WYip646isLRc6lhEVAMWH1EjcbFTIWpKIBysLXAyIx+vb0xEuY4fYktkbFh8RI3I18UW30wKgFIhw67zOXj3t7N8/1kiI8PiI2pkAd5OWDa2OwQB2HAoHav2XZY6EhH9BYuPqAkM6doSbw/rBAD4KDoJW05ckzgREVVi8RE1kRf7tcXUfm0BAP/48SQOXMqTOBERASw+oib1r6GPYmhXd5TpREz/9hiSswqkjkRk9lh8RE1IJhPw6Zju6OXtiILickyJTEC2tljqWERmjcVH1MQsLeT4elIAfFxskKkpRnjkERQUl0kdi8hssfiIDMDBWol1UwLhbKvC+etavPrdcZTxGj8iSbD4iAzE08kaa8MDYGUhR9zFPERsPs1r/IgkwOIjMqBurR3w5YSekMsEbD5+FUt3XZQ6EpHZYfERGdgTHV3x77AuAIDPd1/EpiPpEiciMi8sPiIJjA/0wqwn2gEA3vzlDPYm50iciMh8sPiIJDJ/cHs826MVdHoRr353HGeuaaSORGQWWHxEEhEEAR+O6oa+7VqgsFSHKVFHkHGzUOpYRM0ei49IQkqFDCtf8EdHdzvkFpRgStQRaAp5jR9RU2LxEUnM3tICkVN6wd3eEik5tzHt26MoKddJHYuo2WLxERmBlmorRL3YC3YqBRJSb2L+Dyeh1/MaP6KmwOIjMhId3e2xeqI/LOQCtp66jg+jk6SORNQssfiIjEifds74+LluAICvYi9j3YEr0gYiaoZYfERG5pkerfFGSAcAwLu/ncWOs1kSJyJqXlh8REbo1QG+GB/oBVEEXt+YiOPpt6SORNRsGEXxlZSUoHv37hAEASdOnJA6DpHkBEHA+yM7Y2BHV5SU6/HSuqNIzbsjdSyiZsEoiu+f//wnPDw8pI5BZFQUchm+eL4HurVW4+adUoRHJiDvdonUsYhMnuTFt337duzcuROffPJJvcaXlJRAq9VWW4iaK2ulAmsm94KnkxXSbhRi6rqjKCrlNX5ED0PS4svOzsa0adPw7bffwtraul7bLF68GGq1umrx9PRs4pRE0nKxUyFqSiAcrC1wMiMfr3+fCB2v8SN6YJIVnyiKCA8PxyuvvIKAgIB6b7dw4UJoNJqqJSMjowlTEhkHXxdbfDMpAEqFDDHnsvHeb2f5IbZED6jRiy8iIgKCINS5JCUlYfny5SgoKMDChQsbtH+VSgV7e/tqC5E5CPB2wrKx3SEIwPqDafgq9rLUkYhMkiA28n8bc3NzcePGjTrH+Pj4YMyYMfjtt98gCELVep1OB7lcjgkTJmDdunX1uj2tVgu1Wg2NRsMSJLOwZn8q3t96DgCwbFx3jOzeSuJERNJrSBc0evHVV3p6erUXpmRmZiIkJAQ//fQTgoKC0Lp163rth8VH5uj/fjuHtfGpUMplWD81EI/5tJA6EpGkGtIFCgNluoeXl1e1721tbQEAvr6+9S49InP11rBHkaUtwrbTWXh5/VFsntEHj7jZSR2LyCRIfjkDETWcTCbg0zHdEdDGEdricoRHHkG2tljqWEQmwWiKz9vbG6Ioonv37lJHITIJKoUMX08KgI+LDa7lF2HkF/HYm5wjdSwio2c0xUdEDSMIAhxtlFg3JRDeLayRpS1GeOQRLPjpFLTF/BR3otqw+IhMnKeTNbbPDsaUvt4QBGDT0QyELo1F3MVcqaMRGSUWH1EzYKWUY9HTnfH9tMfg5WSNTE0xJq5JwMKfT6OAsz+ialh8RM1IkE8LRM/pj/A+3gCAjQnpCP0sDvsv5kkbjMiIsPiImhlrpQLvjuiMjdMeQ2tHK1zLL8ILaw7jX7+cxu2ScqnjEUmOxUfUTPX2bYEdc4Ix8bE2AIDvDqcj9LNYHLjE2R+ZNxYfUTNmo1Lg/bAu+O9LQWjlYIWrt4rw/NeH8c6WM7jD2R+ZKRYfkRno084ZO+YGY0JQxTsmrT+YhtBlsTh0ue731SVqjlh8RGbCVqXAB890xYapFbO/jJtFGPfVIbz7v7MoLOXsj8wHi4/IzPR7xBnRc/pjfGDF7C/qwBUMWRaHhNSbEicjMgwWH5EZsrO0wOJnu2Ldi4FoqbZE2o1CjP3qIP7vt3MoKtVJHY+oSbH4iMzY4+1dsGNuMMYGeEIUgbXxqRj6eRyOXuHsj5ovFh+RmbO3tMBHz3VD5JRecLe3RGreHYxefRD/3noOxWWc/VHzw+IjIgDAEx1csWNuMJ7zbw1RBL7Zn4qhy+JwLO2W1NGIGhWLj4iqqK0s8MloP6wND4CrnQqX8+5g9KoDWLztPGd/1Gyw+IjoHgM7uiFm7uN4tmcr6EVgdexlDPs8DonpnP2R6WPxEVGN1NYW+HRMd3wzKQAudipcyr2DUSsP4MPtSZz9kUlj8RFRnQZ1ckPM3GCEdfeAXgRW7buEp5fvx8mMfKmjET0QFh8R3ZeDtRKfjeuB1RP94WyrwsWc23h25QEs2ZGEknLO/si0sPiIqN5COrsjZm4wRvh5QKcXsWLPJYxYHo/TVzVSRyOqNxYfETWIo40Sn4/vgVUv9EQLGyWSswsQ9mU8/rMzGaXleqnjEd0Xi4+IHkhol5bYOTcYw7q1hE4vYvkfKRjxxX6cucbZHxk3Fh8RPbAWtiqseL4nVjzfE042SiRlFSBsRTw+jbnA2R8ZLRYfET20Yd0qZn9DurijXC/i890XMXJFPM5laqWORnQPFh8RNQpnWxW+nNATy8f3gKO1Bc5f12LEF/uxbNdFlOk4+yPjweIjokYjCAKe9vPAzrmPI6SzG8r1IpbuuoCwFfFIyuLsj4wDi4+IGp2LnQqrXvDHsnHd4WBtgbOZWjy9fD+++OMiyjn7I4mx+IioSQiCgJHdW2Hn3GA81ckNZToRn+y8gGe+PIDkrAKp45EZY/ERUZNytbPEVxP9sXSsH9RWFjh9TYOnl+/Hij0pnP2RJFh8RNTkBEHAMz1aY+fcYDzZ0RWlOj2W7EjGqJUHcDGbsz8yLBYfERmMm70lvpkcgP+M9oOdpQInr2ow7PP9WLn3Emd/ZDCSF9/vv/+OoKAgWFlZwdHREWFhYVJHIqImJAgCRvm3RszcxzGggwtKdXp8FJ2E51YdRErObanjkRmQtPg2b96MiRMnYsqUKTh58iTi4+Px/PPPSxmJiAzEXW2JyPBeWPJcN9ipFDiRkY+hn8fhq9hL0OlFqeNRMyaIoijJPay8vBze3t547733MHXq1Afej1arhVqthkajgb29fSMmJCJDua4pQsTm09h3IRcA0NPLAZ+M9oOPi63EychUNKQLJJvxHT9+HNeuXYNMJkOPHj3QsmVLDBkyBGfOnKlzu5KSEmi12moLEZm2lmorRE3phY9GdYWtSoHj6fkYsiwO38Rd5uyPGp1kxXf58mUAwLvvvou33noLW7duhaOjIwYMGICbN2/Wut3ixYuhVqurFk9PT0NFJqImJAgCxvbywo65wej/iDNKyvX49+/nMXb1QaTm3ZE6HjUjjV58EREREAShziUpKQl6fcUruP71r39h1KhR8Pf3R2RkJARBwI8//ljr/hcuXAiNRlO1ZGRkNPYhEJGEWjlYYf2Lgfh/z3SFjVKOo2m3MGRZLNbuT4Wesz9qBIrG3uH8+fMRHh5e5xgfHx9cv34dANCpU6eq9SqVCj4+PkhPT691W5VKBZVK1ShZicg4CYKA54O8ENzeGQs2n0J8yg3839ZziD6ThY+f6wZvZxupI5IJa/Tic3FxgYuLy33H+fv7Q6VSITk5Gf369QMAlJWV4cqVK2jTpk1jxyIiE9Ta0Robpgbhu8Pp+H/bziPhyk0MWRaHBaEdMKm3N2QyQeqIZIIke47P3t4er7zyChYtWoSdO3ciOTkZM2bMAACMHj1aqlhEZGQEQcALj7XBjjnB6O3TAkVlOrz72zmM//oQ0m8USh2PTJCk1/EtWbIE48aNw8SJE9GrVy+kpaXhjz/+gKOjo5SxiMgIeTpZ47uXgvD+yM6wspDjcOpNhC6LxbcHr/C5P2oQya7jayy8jo/I/KTfKMQbP53E4dSKV4D38W2Bj0Z1g6eTtcTJSComcR0fEdGD8mphjY3THsO7T3eCpYUMBy7dQOhnsdhwKA0m/n95MgAWHxGZJJlMQHjftoieHYxe3o64U6rDW7+ewcQ1Cbh6i8/9Ue1YfERk0rydbbDp5d54e3jF7G9/Sh5CP4vDxoR0zv6oRiw+IjJ5MpmAqf3aYvvsYAS0ccTtknIs/Pk0Jq1NQGZ+kdTxyMiw+Iio2WjrbINN03vjrWGPQqWQIe5iHkKWxmLTEc7+6E8sPiJqVuQyAS/198G22f3Rw8sBBSXlWLD5NMIjj+C6hrM/YvERUTPl62KLn17pgzeHdoRSIcO+C7kYvDQWPxzN4OzPzLH4iKjZkssEvBzsi22v94OfpwMKisvxz59O4cWoI8jSFEsdjyTC4iOiZq+dqx02v9IbC0I7QimXYU9yLgYv3YfNx65y9meGWHxEZBYUchlmDPDF76/3g19rNbTF5Zj/40lMW38UOVrO/swJi4+IzMojbnbYPKMP3gjpAAu5gF3nc/DU0lj8mniNsz8zweIjIrOjkMsw84l22Ppaf3RtpYamqAxzNp3A9G+PIbegROp41MRYfERktjq42+HnV/tg/lPtYSEXsPNcNp5aug//O5nJ2V8zxuIjIrNmIZfhtScfwf9m9UOnlvbILyzD6xsTMWPDceTd5uyvOWLxEREBeLSlPbbM6ou5g9pDIRMQfTYLg5fGYuupTKmjUSNj8RER3WUhl2H2oEewZVZfPNrSHjfvlGLWfxPx6nfHcIOzv2aDxUdE9DedPdTYMrMvZj/5CBQyAdtOV8z+tp2+LnU0agQsPiKiGigVMsx9qj1+ndkXHd3tcONOKV797jhm/fc4bt4plToePQQWHxFRHbq0UuN/s/rhtYHtIJcJ2HrqOgYv3YfoM1lSR6MHxOIjIroPpUKG+YM74JdX+6C9my3ybpfilQ3HMPv7RNzi7M/ksPiIiOqpW2sH/PZaP7w6wBcyAdhyIhNPLY3FzrOc/ZkSFh8RUQOoFHL8M7Qjfn61L9q52iLvdgle/vYY5m46gfxCzv5MAYuPiOgBdPd0wNbX+uGVxytmf78kXsPgpbHYfT5b6mh0Hyw+IqIHZGkhR8SQjtg8ow98XWyQU1CCqeuOYv4PJ6EpKpM6HtWCxUdE9JB6eDni99f7Y3qwDwQB2Hz8KgYv3Yc9yTlSR6MasPiIiBqBpYUcC4c+ip9e6Q0fZxtka0swJfII3viRsz9jw+IjImpE/m2csG12f7zUry0EAfjx2FWELI3FXs7+jAaLj4iokVlayPHW8E74YXpveLewRpa2GOGRR7Dgp1MoKObsT2osPiKiJtLL2wnbZwfjxb4Vs79NRzMQsjQWsRdypY5m1lh8RERNyEopxztPd8Kml3ujTQtrZGqKMWltAhb+fJqzP4mw+IiIDCCwrRO2z+6P8D7eAICNCekI/SwO8Sl50gYzQ5IW34ULFzBy5Eg4OzvD3t4e/fr1w549e6SMRETUZKyVCrw7ojM2TnsMnk5WuJZfhAnfHMZbv57GnZJyqeOZDUmLb/jw4SgvL8cff/yBY8eOwc/PD8OHD0dWFt/3joiar96+LRA9OxiTercBAGw4lI6Qz2Jx4BJnf4YgiKIoSnHDeXl5cHFxQWxsLPr37w8AKCgogL29PWJiYjBo0KB67Uer1UKtVkOj0cDe3r4pIxMRNboDKXn45+ZTuHqrCAAw8bE2iBjSETYqhcTJTEtDukCyGV+LFi3QoUMHrF+/Hnfu3EF5eTlWr14NV1dX+Pv717pdSUkJtFpttYWIyFT1aeeM6DnBmBDkBQD49lAaQpfF4tDlGxIna74kKz5BELBr1y4kJibCzs4OlpaW+PTTTxEdHQ1HR8dat1u8eDHUanXV4unpacDURESNz1alwAfPdMWGqUFo5WCFjJtFGPfVIbz7v7MoLOVzf42t0YsvIiICgiDUuSQlJUEURcycOROurq6Ii4tDQkICwsLC8PTTT+P69eu17n/hwoXQaDRVS0ZGRmMfAhGRJPo94ozoOf0xPrBi9hd14AqGLItDQupNiZM1L43+HF9ubi5u3Kh7iu7j44O4uDgMHjwYt27dqvb32EceeQRTp05FREREvW6Pz/ERUXMUeyEXEZtPIVNTDEEApvRpizdCOsBKKZc6mlFqSBc0+rOnLi4ucHFxue+4wsJCAIBMVn3SKZPJoNfrGzsWEZFJCW7vgui5wfhg63lsOpqBtfGp2JOcgyXPdUOAt5PU8UyaZM/x9e7dG46Ojpg8eTJOnjyJCxcu4I033kBqaiqGDRsmVSwiIqNhb2mBj57rhsgpveBub4nUvDsYvfog/r31HIrLdFLHM1mSFZ+zszOio6Nx+/ZtDBw4EAEBAdi/fz+2bNkCPz8/qWIRERmdJzq4YsfcYIz2bw1RBL7Zn4qhy+JwLO2W1NFMkmTX8TUWPsdHROZkT1IOIn4+hWxtCWQC8FJ/H8x7qj0sLcz7uT+TuI6PiIga7omOrtg553E827MV9CLwVexlDPs8DonpnP3VF4uPiMjEqK0t8OmY7vhmUgBc7FS4lHsHo1YewIfbk/jcXz2w+IiITNSgTm6ImRuMZ3pUzP5W7buEp5fvx8mMfKmjGTUWHxGRCXOwVmLp2O74aqI/nG1VuJhzG8+uPIAlO5JQUs7ZX01YfEREzcDgzu6ImRuMkd09oNOLWLHnEkYsj8fpqxqpoxkdFh8RUTPhaKPEsnE9sOoFfzjbKpGcXYCwL+Pxn53JKC3nG4NUYvERETUzoV3csXPu4xjerSV0ehHL/0jBiC/248w1zv4AFh8RUbPkZKPEF8/3xJcTesLJRomkrAKErYjHpzEXzH72x+IjImrGhnZtiZi5wRja1R3lehGf776IkSvicTbTfGd/LD4iomauha0KX07wxxfP94CjtQXOX9di5BfxWLbrIsp05jf7Y/EREZmJ4d08sHPu4wjtXDH7W7rrAsJWxOP8da3U0QyKxUdEZEZc7FRY+UJPfD6+BxysLXA2U4sRX+zH8t3mM/tj8RERmRlBEDDCzwM75wbjqU5uKNOJ+E/MBTzzZTySswqkjtfkWHxERGbK1c4SX030x2dju0NtZYEz17QYvjwOK/akoLwZz/5YfEREZkwQBIT1aIWYucEY9KgrynQiluxIxrMrD+BidvOc/bH4iIgIrvaW+HpSAD4d4wd7SwVOXdVg2Of7sXLvpWY3+2PxERERgIrZ37M9WyNm3uMY2NEVpTo9PopOwqhVB5GS03xmfyw+IiKqxs3eEmsmB2DJc91gZ6nAyYx8DP18P1bvuwSdXpQ63kNj8RER0T0EQcDoAE/snBuMAR1cUFqux+LtSXhu1QFcyr0tdbyHwuIjIqJatVRbITK8Fz4e1Q12KgUS0/MxdFkcvo69bLKzPxYfERHVSRAEjOnliR1zg9H/EWeUlOvxwbbzGLP6IC6b4OyPxUdERPXi4WCF9S8G4sNnu8JWpcCxtFsYsiwOa/anQm9Csz8WHxER1ZsgCBgX6IUdc4PRr13F7O/9recw9quDuJJ3R+p49cLiIyKiBmvlYIVvpwbig2e6wEYpx5ErtxC6LBaR8cY/+2PxERHRAxEEAROC2iB6TjD6+LZAcZke7/12DuO/PoS0G8Y7+2PxERHRQ/F0ssaGqUF4f2RnWCvlOJx6E6GfxWH9wStGOftj8RER0UOTyQRM7O2N6NnBCGrrhKIyHd7ZchbPf3MIGTcLpY5XDYuPiIgajVcLa2yc9hjeG9EZVhZyHLp8EyGfxeLbQ2lGM/tj8RERUaOSyQRM7uON6Dn9EejthMJSHd7+9Qwmrj2Mq7ekn/2x+IiIqEm0aWGD719+DIue7gRLCxniU24gZGks/ns4HaIo3eyPxUdERE1GJhMwpW9bbJ8djIA2jrhTqsObv5zGpLUJuJZfJE0mSW6ViIjMSltnG2ya3htvDXsUKoUMcRfzELI0FpuOGH7212TF98EHH6BPnz6wtraGg4NDjWPS09MxbNgwWFtbw9XVFW+88QbKy8ubKhIREUlILhPwUn8fbJvdHz29HHC7pBwLNp9GeOQRXNcYbvbXZMVXWlqK0aNHY8aMGTX+XKfTYdiwYSgtLcWBAwewbt06REVF4Z133mmqSEREZAR8XWzx4yt98ObQjlAqZNh3IReDl8YiW1tskNsXxCaeY0ZFRWHOnDnIz8+vtn779u0YPnw4MjMz4ebmBgBYtWoVFixYgNzcXCiVyhr3V1JSgpKSkqrvtVotPD09odFoYG9v32THQUREjS8l5zb+8eNJ+LjY4NMx3R94P1qtFmq1ul5dINlzfAcPHkTXrl2rSg8AQkJCoNVqcfbs2Vq3W7x4MdRqddXi6elpiLhERNQE2rna4qdXeuPfYV0MdpuSFV9WVla10gNQ9X1WVlat2y1cuBAajaZqycjIaNKcRETUtBRyGayVCoPdXoOKLyIiAoIg1LkkJSU1VVYAgEqlgr29fbWFiIiovhpUsfPnz0d4eHidY3x8fOq1L3d3dyQkJFRbl52dXfUzIiKiptCg4nNxcYGLi0uj3HDv3r3xwQcfICcnB66urgCAmJgY2Nvbo1OnTo1yG0RERH/XZH9UTU9Px82bN5Geng6dTocTJ04AANq1awdbW1sMHjwYnTp1wsSJE/Hxxx8jKysLb731FmbOnAmVStVUsYiIyMw12eUM4eHhWLdu3T3r9+zZgwEDBgAA0tLSMGPGDOzduxc2NjaYPHkyPvzwQygU9e/jhryElYiImqeGdEGTX8fX1Fh8RETUkC4w3OtHm0hlb2u1WomTEBGRVCo7oD5zOZMvvoKCAgDghexERISCggKo1eo6x5j8nzr1ej0yMzNhZ2cHQRAeeD+Vb32WkZFhsn8yNfVjMPX8AI/BWPAYjIMhj0EURRQUFMDDwwMyWd2XqJv8jE8mk6F169aNtr/mcFG8qR+DqecHeAzGgsdgHAx1DPeb6VXi5/EREZFZYfEREZFZYfHdpVKpsGjRIpO+eN7Uj8HU8wM8BmPBYzAOxnoMJv/iFiIioobgjI+IiMwKi4+IiMwKi4+IiMwKi4+IiMwKi4+IiMyKWRXfihUr4O3tDUtLSwQFBd3zCfB/9+OPP6Jjx46wtLRE165dsW3bNgMlvdfixYvRq1cv2NnZwdXVFWFhYUhOTq5zm6ioKAiCUG2xtLQ0UOJ7vfvuu/fk6dixY53bGNM5AABvb+97jkEQBMycObPG8VKfg9jYWDz99NPw8PCAIAj49ddfq/1cFEW88847aNmyJaysrDBo0CBcvHjxvvtt6GPpYdR1DGVlZViwYAG6du0KGxsbeHh4YNKkScjMzKxznw9yX2yqYwAqPsbt73lCQ0Pvu19jOQ8AanxcCIKAJUuW1LpPQ5+HSmZTfJs2bcK8efOwaNEiHD9+HH5+fggJCUFOTk6N4w8cOIDx48dj6tSpSExMRFhYGMLCwnDmzBkDJ6+wb98+zJw5E4cOHUJMTAzKysowePBg3Llzp87t7O3tcf369aolLS3NQIlr1rlz52p59u/fX+tYYzsHAHDkyJFq+WNiYgAAo0ePrnUbKc/BnTt34OfnhxUrVtT4848//hiff/45Vq1ahcOHD8PGxgYhISEoLi6udZ8NfSw15TEUFhbi+PHjePvtt3H8+HH8/PPPSE5OxogRI+6734bcFx/W/c4DAISGhlbLs3Hjxjr3aUznAUC17NevX8fatWshCAJGjRpV534NeR6qiGYiMDBQnDlzZtX3Op1O9PDwEBcvXlzj+DFjxojDhg2rti4oKEicPn16k+asr5ycHBGAuG/fvlrHREZGimq12nCh7mPRokWin59fvccb+zkQRVGcPXu26OvrK+r1+hp/bkznAID4yy+/VH2v1+tFd3d3ccmSJVXr8vPzRZVKJW7cuLHW/TT0sdSY/n4MNUlISBABiGlpabWOaeh9sTHVdAyTJ08WR44c2aD9GPt5GDlypDhw4MA6x0h1HsxixldaWopjx45h0KBBVetkMhkGDRqEgwcP1rjNwYMHq40HgJCQkFrHG5pGowEAODk51Tnu9u3baNOmDTw9PTFy5EicPXvWEPFqdfHiRXh4eMDHxwcTJkxAenp6rWON/RyUlpZiw4YNePHFF+v8ZBBjOweVUlNTkZWVVe13rFarERQUVOvv+EEeS4am0WggCAIcHBzqHNeQ+6Ih7N27F66urujQoQNmzJiBGzdu1DrW2M9DdnY2fv/9d0ydOvW+Y6U4D2ZRfHl5edDpdHBzc6u23s3NDVlZWTVuk5WV1aDxhqTX6zFnzhz07dsXXbp0qXVchw4dsHbtWmzZsgUbNmyAXq9Hnz59cPXqVQOm/VNQUBCioqIQHR2NlStXIjU1Ff3796/6TMW/M+ZzAAC//vor8vPzER4eXusYYzsHf1X5e2zI7/hBHkuGVFxcjAULFmD8+PF1fhpAQ++LTS00NBTr16/H7t278dFHH2Hfvn0YMmQIdDpdjeON/TysW7cOdnZ2ePbZZ+scJ9V5MPmPJTJHM2fOxJkzZ+77t/DevXujd+/eVd/36dMHjz76KFavXo3333+/qWPeY8iQIVVfd+vWDUFBQWjTpg1++OGHev3P0NisWbMGQ4YMgYeHR61jjO0cNGdlZWUYM2YMRFHEypUr6xxrbPfFcePGVX3dtWtXdOvWDb6+vti7dy+efPJJg+d5WGvXrsWECRPu+0Iuqc6DWcz4nJ2dIZfLkZ2dXW19dnY23N3da9zG3d29QeMNZdasWdi6dSv27NnT4M8htLCwQI8ePZCSktJE6RrGwcEB7du3rzWPsZ4DAEhLS8OuXbvw0ksvNWg7YzoHlb/HhvyOH+SxZAiVpZeWloaYmJgGf/bb/e6Lhubj4wNnZ+da8xjreQCAuLg4JCcnN/ixARjuPJhF8SmVSvj7+2P37t1V6/R6PXbv3l3tf+N/1bt372rjASAmJqbW8U1NFEXMmjULv/zyC/744w+0bdu2wfvQ6XQ4ffo0WrZs2QQJG+727du4dOlSrXmM7Rz8VWRkJFxdXTFs2LAGbWdM56Bt27Zwd3ev9jvWarU4fPhwrb/jB3ksNbXK0rt48SJ27dqFFi1aNHgf97svGtrVq1dx48aNWvMY43motGbNGvj7+8PPz6/B2xrsPBj85TQS+f7770WVSiVGRUWJ586dE19++WXRwcFBzMrKEkVRFCdOnChGRERUjY+PjxcVCoX4ySefiOfPnxcXLVokWlhYiKdPn5Yk/4wZM0S1Wi3u3btXvH79etVSWFhYNebvx/Dee++JO3bsEC9duiQeO3ZMHDdunGhpaSmePXtWikMQ58+fL+7du1dMTU0V4+PjxUGDBonOzs5iTk5OjfmN7RxU0ul0opeXl7hgwYJ7fmZs56CgoEBMTEwUExMTRQDip59+KiYmJla94vHDDz8UHRwcxC1btoinTp0SR44cKbZt21YsKiqq2sfAgQPF5cuXV31/v8eSIY+htLRUHDFihNi6dWvxxIkT1R4bJSUltR7D/e6LhjyGgoIC8R//+Id48OBBMTU1Vdy1a5fYs2dP8ZFHHhGLi4trPQZjOg+VNBqNaG1tLa5cubLGfUh9HiqZTfGJoiguX75c9PLyEpVKpRgYGCgeOnSo6mePP/64OHny5Grjf/jhB7F9+/aiUqkUO3fuLP7+++8GTvwnADUukZGRVWP+fgxz5sypOl43Nzdx6NCh4vHjxw0f/q6xY8eKLVu2FJVKpdiqVStx7NixYkpKStXPjf0cVNqxY4cIQExOTr7nZ8Z2Dvbs2VPj/aYyo16vF99++23Rzc1NVKlU4pNPPnnPcbVp00ZctGhRtXV1PZYMeQypqam1Pjb27NlT6zHc775oyGMoLCwUBw8eLLq4uIgWFhZimzZtxGnTpt1TYMZ8HiqtXr1atLKyEvPz82vch9TnoRI/j4+IiMyKWTzHR0REVInFR0REZoXFR0REZoXFR0REZoXFR0REZoXFR0REZoXFR0REZoXFR0REZoXFR0REZoXFR0REZoXFR0REZuX/AxGMYGS+S+uxAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate aquifer parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".....................................\n",
      "Fit succeeded.\n"
     ]
    }
   ],
   "source": [
    "cal = tft.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq\", layers=list(range(1, nlay)), initial=10, pmin=0)\n",
    "cal.set_parameter(name=\"Saq\", layers=list(range(1, nlay)), initial=1e-4, pmin=0)\n",
    "cal.seriesinwell(name=\"obs\", element=w, t=to, h=ho)\n",
    "cal.fit(report=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>optimal</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>kaq_1_17</th>\n",
       "      <td>0.439172</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq_1_17</th>\n",
       "      <td>0.000461</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           optimal\n",
       "kaq_1_17  0.439172\n",
       "Saq_1_17  0.000461"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 0.006 m\n"
     ]
    }
   ],
   "source": [
    "display(cal.parameters.loc[:, [\"optimal\"]])\n",
    "print(f\"RMSE: {cal.rmse():.3f} m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAFBCAYAAAD36+/HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABZTElEQVR4nO3dd1gUVxfA4d+ydBBQqSqCYu9G7D1RMRqjsWGJLWpiSdSgsSRRY4qY2GMs0cQSKxqNJmpsKPZuNGrsglgQ5VNAQWHZne+PDRsRxAWBpZz3eeaBnZ2ZPcsOHO6dO+eqFEVREEIIIUSazEwdgBBCCJGbSaIUQggh0iGJUgghhEiHJEohhBAiHZIohRBCiHRIohRCCCHSIYlSCCGESIckSiGEECIdkiiFEEKIdEiiFNlCpVLxxRdfZHi/sLAwVCoVS5cuzfKYslNISAgqlYqQkBBTh5LjMvLeM/L59u3bF29v71eOryBaunQpKpWKsLCwDO/7xRdfoFKpsj6oPEwSZT6W/MuiUqk4cOBAqucVRcHT0xOVSsVbb71lggjzt1WrVjFr1ixTh2ESBfm9i/xHEmUBYG1tzapVq1Kt37t3L7du3cLKysoEUeV/BTlZvOi9e3l58eTJE3r16pXzQQmRSZIoC4A2bdqwbt06kpKSUqxftWoVtWrVwt3d3USRZa24uDhThyBeQqVSYW1tjVqtNnUor0RRFJ48eWLqMEQOkURZAHTv3p3//e9/7Ny507AuMTGRX3/9lR49eqS5T1xcHCNHjsTT0xMrKyvKly/PtGnTeH6ymYSEBD7++GNcXFwoVKgQb7/9Nrdu3UrzmLdv3+a9997Dzc0NKysrKleuzOLFizP1npK7lffu3cuQIUNwdXWlRIkShuf//PNPGjdujJ2dHYUKFaJt27acP38+xTHu3r1Lv379KFGiBFZWVnh4eNC+ffsU13VedK3V29ubvn37vjC+Zs2asWXLFm7cuGHo/n72etucOXOoXLkytra2FC5cGF9f3zRb/ZmVfC1w2rRpzJ07l9KlS2Nra0urVq24efMmiqLw1VdfUaJECWxsbGjfvj0PHjxIcYzseO+veg162rRpNGjQgKJFi2JjY0OtWrX49ddfU2zTtGlTqlevnub+5cuXx8/Pz/BYp9Mxa9YsKleujLW1NW5ubnzwwQc8fPgw1Xt+66232L59O76+vtjY2PDjjz+m+zOoUqUKf//9N02bNsXW1pYyZcoYYt27dy9169bFxsaG8uXLs2vXrlTH+Ouvv3jzzTdxcHDA3t6eN954gyNHjqTa7vz587z++uvY2NhQokQJvv76a3Q6XZpxGfN7IVIzN3UAIvt5e3tTv359Vq9ezZtvvgnof2FiYmLo1q0b33//fYrtFUXh7bffZs+ePfTv358aNWqwfft2PvnkE27fvs3MmTMN2w4YMIAVK1bQo0cPGjRowO7du2nbtm2qGCIjI6lXrx4qlYoPP/wQFxcX/vzzT/r3709sbCwjRozI1HsbMmQILi4uTJgwwdCiXL58OX369MHPz49vv/2W+Ph45s+fT6NGjfjrr78Mf7Q7derE+fPn+eijj/D29ubevXvs3LmT8PDwVx5E8tlnnxETE8OtW7cMPy97e3sAFi1axLBhw+jcuTPDhw/n6dOn/P333xw9evSF/7hk1sqVK0lMTOSjjz7iwYMHfPfdd3Tt2pXXX3+dkJAQxowZw9WrV5kzZw6jRo3K9D8uz0rvvb+q2bNn8/bbb9OzZ08SExNZs2YNXbp0YfPmzYbzrlevXgwcOJBz585RpUoVw77Hjx/n8uXLfP7554Z1H3zwAUuXLqVfv34MGzaM0NBQfvjhB/766y8OHjyIhYWFYdtLly7RvXt3PvjgAwYOHEj58uXTjfXhw4e89dZbdOvWjS5dujB//ny6devGypUrGTFiBIMGDaJHjx5MnTqVzp07c/PmTQoVKgTok1/jxo1xcHBg9OjRWFhY8OOPP9KsWTNDkgX9P3vNmzcnKSmJsWPHYmdnx8KFC7GxsUkVj7G/FyINisi3lixZogDK8ePHlR9++EEpVKiQEh8fryiKonTp0kVp3ry5oiiK4uXlpbRt29aw38aNGxVA+frrr1Mcr3PnzopKpVKuXr2qKIqinD59WgGUIUOGpNiuR48eCqBMnDjRsK5///6Kh4eHEhUVlWLbbt26KY6Ojoa4QkNDFUBZsmSJUe+tUaNGSlJSkmH9o0ePFCcnJ2XgwIEptr97967i6OhoWP/w4UMFUKZOnZru6zz/PpJ5eXkpffr0MTzes2ePAih79uwxrGvbtq3i5eWVat/27dsrlStXTvd1X1Xyz9HFxUWJjo42rB83bpwCKNWrV1c0Go1hfffu3RVLS0vl6dOnhnXZ8d6N/XwVRVH69OmT6hjJ50myxMREpUqVKsrrr79uWBcdHa1YW1srY8aMSbHtsGHDFDs7O+Xx48eKoijK/v37FUBZuXJliu22bduWar2Xl5cCKNu2bXtp3IqiKE2bNlUAZdWqVYZ1Fy9eVADFzMxMOXLkiGH99u3bU/1MOnTooFhaWirXrl0zrLtz545SqFAhpUmTJoZ1I0aMUADl6NGjhnX37t1THB0dFUAJDQ1VFMX43wtFUZSJEycqkhpSkq7XAqJr1648efKEzZs38+jRIzZv3vzC1svWrVtRq9UMGzYsxfqRI0eiKAp//vmnYTsg1XbPtw4VRWH9+vW0a9cORVGIiooyLH5+fsTExHDq1KlMva+BAwemuN61c+dOoqOj6d69e4rXUavV1K1blz179gBgY2ODpaUlISEhqbrZspuTkxO3bt3i+PHj2f5aXbp0wdHR0fA4uSXy7rvvYm5unmJ9YmIit2/fzvaYXsWzLaWHDx8SExND48aNU5w/jo6OtG/fntWrVxsuFWi1WoKCgujQoQN2dnYArFu3DkdHR1q2bJniXKlVqxb29vaGcyVZqVKlUnTbvoy9vT3dunUzPC5fvjxOTk5UrFjR8DnAf5/J9evXDbHu2LGDDh06ULp0acN2Hh4e9OjRgwMHDhAbGwvofwfr1atHnTp1DNu5uLjQs2fPFLEY+3sh0iZdrwWEi4sLLVq0YNWqVcTHx6PVauncuXOa2964cYNixYoZuoGSVaxY0fB88lczMzN8fHxSbPd8l9T9+/eJjo5m4cKFLFy4MM3XvHfvXqbeV6lSpVI8vnLlCgCvv/56mts7ODgAYGVlxbfffsvIkSNxc3OjXr16vPXWW/Tu3TvbBzeNGTOGXbt2UadOHcqUKUOrVq3o0aMHDRs2THe/u3fvpnjs6OiYZhfbs0qWLJlqHwBPT8801+f0Pw1PnjwhJiYmxbr0fv6bN2/m66+/5vTp0yQkJBjWP3/fX+/evQkKCmL//v00adKEXbt2ERkZmWK07ZUrV4iJicHV1TXN13r+nHz+XHuZEiVKpIrL0dHxpT/7+/fvEx8fn2bXbsWKFdHpdNy8eZPKlStz48aNFEk32fP7Gvt7IdImibIA6dGjBwMHDuTu3bu8+eabODk55cjrJg8sePfdd+nTp0+a21SrVi1Tx34+USS/1vLly9P8g/tsK2rEiBG0a9eOjRs3sn37dsaPH09gYCC7d++mZs2a6b6uVqvNVLyg/2N36dIlNm/ezLZt21i/fj3z5s1jwoQJTJo06YX7eXh4pHi8ZMmSdAfVAC8cXfqi9cpzg7XS8irv/XlBQUH069fPqBj279/P22+/TZMmTZg3bx4eHh5YWFiwZMmSVAOh/Pz8cHNzY8WKFTRp0oQVK1bg7u5OixYtDNvodDpcXV1ZuXJlmq/n4uKS4vHL/il5Xnb87DMrI78XIjX56RQg77zzDh988AFHjhwhKCjohdt5eXmxa9cuHj16lKJVefHiRcPzyV91Oh3Xrl1L8R/spUuXUhwveUSsVqtN8YcqOyS3bl1dXY16LR8fH0aOHMnIkSO5cuUKNWrUYPr06axYsQKAwoULEx0dnWKfxMREIiIiXnrs9Kqb2NnZ4e/vj7+/P4mJiXTs2JFvvvmGcePGYW1tneY+z45aBqhcufJLY3gV2fXen+Xn55fqfb3I+vXrsba2Zvv27Snu/V2yZEmqbdVqNT169GDp0qV8++23bNy4MVU3vY+PD7t27aJhw4YZToLZycXFBVtb21S/R6D/HTQzMzO0Sr28vAytxWc9v29Gfy9ESnKNsgCxt7dn/vz5fPHFF7Rr1+6F27Vp0watVssPP/yQYv3MmTNRqVSGkbPJX58fNfv8jeZqtZpOnTqxfv16zp07l+r17t+/n5m3kyY/Pz8cHByYPHkyGo3mha8VHx/P06dPUzzn4+NDoUKFUnTp+fj4sG/fvhTbLVy40KhWlZ2dXapuRYD//e9/KR5bWlpSqVIlFEVJM+ZkLVq0SLE838LMatnx3p/n4eGR6n29iFqtRqVSpXj9sLAwNm7cmOb2vXr14uHDh3zwwQc8fvyYd999N8XzXbt2RavV8tVXX6XaNykpKdU/CTlFrVbTqlUrNm3alOJWpcjISFatWkWjRo0MXaVt2rThyJEjHDt2zLDd/fv3U7WSjf29EGmTFmUB86Kuz2e1a9eO5s2b89lnnxEWFkb16tXZsWMHmzZtYsSIEYb/TmvUqEH37t2ZN28eMTExNGjQgODgYK5evZrqmFOmTGHPnj3UrVuXgQMHUqlSJR48eMCpU6fYtWtXqnv4MsvBwYH58+fTq1cvXnvtNbp164aLiwvh4eFs2bKFhg0b8sMPP3D58mXeeOMNunbtSqVKlTA3N+e3334jMjIyxQCMAQMGMGjQIDp16kTLli05c+YM27dvx9nZ+aWx1KpVi6CgIAICAqhduzb29va0a9eOVq1a4e7uTsOGDXFzc+PChQv88MMPtG3bNtV1YVPKjvf+Ktq2bcuMGTNo3bo1PXr04N69e8ydO5cyZcrw999/p9q+Zs2aVKlShXXr1lGxYkVee+21FM83bdqUDz74gMDAQE6fPk2rVq2wsLDgypUrrFu3jtmzZ7/wOn52+/rrr9m5cyeNGjViyJAhmJub8+OPP5KQkMB3331n2G706NEsX76c1q1bM3z4cMPtIV5eXil+Jsb+XogXMN2AW5Hdnr09JD3P3x6iKPrh5B9//LFSrFgxxcLCQilbtqwydepURafTpdjuyZMnyrBhw5SiRYsqdnZ2Srt27ZSbN2+meWtBZGSkMnToUMXT01OxsLBQ3N3dlTfeeENZuHChYZuM3h7yove2Z88exc/PT3F0dFSsra0VHx8fpW/fvsqJEycURVGUqKgoZejQoUqFChUUOzs7xdHRUalbt66ydu3aFMfRarXKmDFjFGdnZ8XW1lbx8/NTrl69atQtEo8fP1Z69OihODk5KYDhVocff/xRadKkiVK0aFHFyspK8fHxUT755BMlJiYm3fecEck/x+dvf0mOc926dSnWp/XzzI73/qq3h/z8889K2bJlFSsrK6VChQrKkiVL0r2d4bvvvlMAZfLkyS98nYULFyq1atVSbGxslEKFCilVq1ZVRo8erdy5c8ewTVq/I+lp2rRpmrcAveg4gDJ06NAU606dOqX4+fkp9vb2iq2trdK8eXPl0KFDqfb9+++/laZNmyrW1tZK8eLFla+++kr5+eefU9wekuxlvxeKIreHpEWlKNl4BVkIIUxo9uzZfPzxx4SFhaUaASyEsSRRCiHyJUVRqF69OkWLFpX7BMUrkWuUQoh8JS4ujt9//509e/Zw9uxZNm3aZOqQRB4nLUohRL4SFhZGqVKlcHJyYsiQIXzzzTemDknkcZIohRBCiHTIfZRCCCFEOiRRCiGEEOkocIN5dDodd+7coVChQkaX2RJCCJH/KIrCo0ePKFasGGZm6bQbTXUDp6Ioyt69e5W33npL8fDwUADlt99+e+k+e/bsUWrWrKlYWloqPj4+Rt24/Kzkm+FlkUUWWWSRBVBu3ryZbt4waYsyLi6O6tWr895779GxY8eXbh8aGkrbtm0ZNGgQK1euJDg4mAEDBuDh4WH0PHHJJcJu3rwpU8vkMxqNhh07dhhKkQmRU+Tcy5tiY2Px9PR8aelIkybKN99801BY2xgLFiygVKlSTJ8+HdBPV3TgwAFmzpxpdKJM7m51cHCQRJnPaDQabG1tcXBwkD9WIkfJuZe3vewyXJ66Rnn48OFUswv4+fkxYsSIF+6TkJCQYjaI5JnBNRpNujM1iLwn+fOUz1XkNDn38iZjP688lSjv3r2Lm5tbinVubm7Exsby5MmTNOeUCwwMTHMy3B07dmBra5ttsQrTMXZ+QyGympx7eUt8fLxR2+WpRJkZ48aNIyAgwPA4uU+6VatW0vWaz2g0Gnbu3EnLli2l+0vkKDn38qbkHsaXyVOJ0t3dncjIyBTrIiMjcXBweOEM5VZWVilmQ09mYWEhJ3Q+JZ+tMJXkc0+r1Uo3bC5gYWGBWq1O93lj5KlEWb9+fbZu3Zpi3c6dO6lfv76JIhJCiP8oikJERATR0dGmDkX8y8nJCXd391e6b96kifLx48dcvXrV8Dg0NJTTp09TpEgRSpYsybhx47h9+za//PILAIMGDeKHH35g9OjRvPfee+zevZu1a9eyZcuWHI89IuYJoVFxlHK2w8Mx7dasEKJguXfvHo8ePcLV1RVbW1spamJCiqIQHx/PvXv3APDw8Mj0sUyaKE+cOEHz5s0Nj5OvJfbp04elS5cSERFBeHi44flSpUqxZcsWPv74Y2bPnk2JEiX46aefjL41JKsEHQ9n3Iaz6BQwU0Fgx6r415ZJYYUoyFQqFbGxsbi5uVG0aFFThyPAcEnu3r17uLq6ptsNmx6TJspmzZqhpDN5ydKlS9Pc56+//srGqNIXEfPEkCQBdAp8uuEcTcq5SMtSiAIs+Y+wjKbPXZI/D41Gk+lEKUXRMyg0Ks6QJJNpFYWwKOOGGQsh8jfpbs1dsuLzkESZQaWc7TB77ueuVqnwdpb/IoUQIj+SRJlBHo42BHasivrf/1LUKhWTO1bJc92uETFPOHQtioiYJ6YORQiRi4WEhKBSqQr0SN48dXtIbuFfw5U3H18mQilKUQ9vnIuZgTYJ1HnjxymDkYQQwnh54y97bhN7G4e9E0lR10dlBvZu4FAcHIqBoyc4FgfHEv8unmDnAia+fiGDkYQQImMkUWaGSgWVO0LsHXh0B2IjQKeBRxH65fYL9jO3/i9pOpX8d/HSfy3spU+02ZxI0xuMJIlSiNwjJ+/VTkhI4JNPPmHNmjXExsbi6+vLzJkzqV27tmGbgwcPMm7cOC5fvkyNGjX46aefqFKlCgA3btzgww8/5MCBAyQmJuLt7c3UqVNp06ZNtsadUyRRZkaR0tBlyX+PdTqIuw+xt/VLzG2IvQUxt/Tfx9yER3ch6Sn876p+SYu5tT5xFvb+d/GCwqWgSCn9Y4tX/2VJHoz0bLKUwUhC5C45fXlk9OjRrF+/nmXLluHl5cV3332Hn59fioIwn3zyCbNnz8bd3Z1PP/2Udu3acfnyZSwsLBg6dCiJiYns27cPOzs7/vnnH+zt7bMt3pwmiTIrmJlBITf9Uvy1tLdJSvw3id6E6Jv/fg2Hhzf0X2Nv6RNp1CX9kpZCxaCojz5xFvHRf1+0jD6ZWlgbFWryYKRPN5xDqyh5djCSEPlVTl8eiYuLY/78+SxdutQwP/CiRYvYuXMnP//8s6FVOXHiRFq2bAnAsmXLKFGiBL/99htdu3YlPDycTp06UbVqVQBKly6d5XGakiTKnGJu+W+CK5X281qNvgX6MOyZJVT/9UEoJMTqu3kf3YGw/c/trAInT33SdC7331fnclDIPVV3rn/tkjQp50JYVDzezraSJIXIRXL68si1a9fQaDQ0bNjQsM7CwoI6depw4cIFQ6J8tqZ2kSJFKF++PBcuXABg2LBhDB48mB07dtCiRQs6depEtWrVsjxWU5FEmVuoLV6cSBUF4h/Ag+vPLNfgf9f03bgJsfpWaXQ4XNudcl8rR3AuCy7l/10qgEt5PBxLSoIUIhfKi5dHBgwYgJ+fH1u2bGHHjh0EBgYyffp0PvroI1OHliUkUeYFKhXYFdUvnrVTPqcoEBcF/7sCUVf++xp1Rd8iTYiB2yf0y7Ms7PSJ07UiuFYCt0rgWhnsXU0+MleIgiynL4/4+PhgaWnJwYMH8fLyAvTl3o4fP86IESMM2x05coSSJfXXSR8+fMjly5epWLGi4XlPT08GDRrEoEGDGDduHIsWLZJEKXIJlQrsXfSLV4OUzyUl6Fuf9y/C/cv/fr2kT6aaOLhzSr88y9YZ3CqDe9X/vjqX13cdCyFyRE5eHrGzs2Pw4MF88sknhpmbvvvuO+Lj4+nfvz9nzpwB4Msvv6Ro0aK4ubnx2Wef4ezsTIcOHQAYMWIEb775JuXKlePhw4fs2bMnRRLN6yRR5mfmVv+2GJ87YbVJ+gR67x+4dwHunYfIf/Tr4qMgdK9+SWZmAa4VwL06ePy7uFcBS7ucfT9CFCAejjY5dnlkypQp6HQ6evXqxaNHj/D19WX79u0ULlw4xTbDhw/nypUr1KhRgz/++ANLS/0/0FqtlqFDh3Lr1i0cHBxo3bo1M2fOzJHYc4JKSW/6jnwoNjYWR0dHYmJicHBwePkOBUliPNy/AJHn9cvdc3D3rL77NhWVfrBQsZr/Le5VwdJ011E0Gg1bt26lTZs2Rs9cLkRW0Gg07Nixg1KlSlG6dGmsrY0bhS6y39OnTwkNDaVUqVKpPhdj84G0KMV/LG2heC39kkxR9IOE7p6FiDNw92/910cR/93K8vca/bYqtf56Z/HX9Mco4asfPGSWualthBAiN5BEKdKnUv1b+MALKr713/pHkfqEeeevf5dT8DgSIs/ql1PL9NtZ2kOxmjx2qclN+6oULt8Ad/cSpnkvQgiRCZIoReYUcoNCraBcK/1jRdGX9LtzCm6f/Hc5BYmPIWw/9mH7qQiwB2LtvHEo2whK1tMvRcvISFshRK4liVJkDZXq3yLwxaFiO/06nZb7oWeYsXglNVVXeM3sCmXM7uAQFwanw+D0Cv12ts7gVR9KNtCP3HWvKt21QohcQxKlyD5maq5QktXa11nN6wA48pjXzK7w1WuPKRF7Rt/yjI+CC3/oF9AXSfBqAN6NoFRjcKuqLxMohBAmIIlSZKvnq4zEYM8+5TXULZuDo43+Xs87f8GNg3DjMIQf0Y+yvfynfgGwdvo3aTaF0s30lYakq1YIkUMkUYps9dIqI+ZW/12rbIz+Hs+7ZyDsIIQdgBuH4Gk0XNysX0BfHL50M/BpDqWb64stCCFENpFEKbJdhqqMqM3/u0Wl4TB94ow4rS+AcD0Ewo/qC8OfWaVfANyrQZk3UHk3Q6Uk5cRbEkIUIJIoRY7IdJURtbn+fswSvtB4JGie6Ltnr+/RF4C/e1Z/b+fdvzE/MJM3zWxQP9mgH41bthU4eGT9mxFCFCgyQkLkLRY2+i7Xll/CoAMw6gp0XATV/FFsi2Khe4LZpc3wxzCYUQF+bAJ7JutvVSlYRaiEyLSQkBBUKhXR0dGvdJz4+Hg6deqEg4OD4Xje3t7MmjUrS+LMKRlqUep0Ovbu3cv+/fu5ceMG8fHxuLi4ULNmTVq0aIGnp2d2xSlE2uxdoVpXqNaVpMQEDv06j0Zu8aivB+tH1Eac0S97v4VCHlD+TSjfBko10V8fFULQrFkzatSoYUhgDRo0ICIiAkdHx1c67rJly9i/fz+HDh3C2dn5lY9nKka1KJ88ecLXX3+Np6cnbdq04c8//yQ6Ohq1Ws3Vq1eZOHEipUqVok2bNhw5ciS7YxYibSozou1Ko2syGgbu1rc228/T39dpYacvu3diMazsDN/5wLp+cG49JDwydeRC5CqWlpa4u7ujesXR5deuXaNixYpUqVIlS45nKkYlynLlyvH333+zaNEiYmNjOXz4MOvXr2fFihVs3bqV8PBwrl27RuPGjenWrRuLFi3K7riFeDl7V6jZE/xXwOjr0HM9+L6H1s4dEh/B+Q3w63v6pLm6O5xZA0+iTR21EDmqb9++7N27l9mzZ6NSqVCpVCxdujRF1+vSpUtxcnJi8+bNlC9fHltbWzp37kx8fDzLli3D29ubwoULM2zYMLRaLaBvpU6fPp19+/ahUqlo1qxZmq8fHh5O+/btsbe3x8HBga5duxIZGQlATEwMarWaEyf08+nqdDqKFClCvXr1DPuvWLEi23szjep63bFjx0vnFvPy8mLcuHGMGjWK8PDwLAlOiCxjYQ1lWxAUXY5PD75OVa7TWn2cno5/UyjuBlzaql/MLMDndajSUd9Fay0zzIhXoCigiTfNa1vYGnW/8ezZs7l8+TJVqlThyy+/BOD8+fOptouPj+f7779nzZo1PHr0iI4dO/LOO+/g5OTE1q1buX79Op06daJhw4b4+/uzYcMGxo4dy7lz59iwYYNhSq5n6XQ6Q5Lcu3cvSUlJDB06FH9/f0JCQnB0dKRGjRqEhITg6+vL2bNnUalU/PXXXzx+/NiwX9OmTV/955UOoxJlRibgtLCwwMfHJ9MBCZFdImKeMG7DWXSKGacpw+mkMkx90J0jAzxwCd8O/2zSTzN2Zbt+UVtB2ZZQpROUa23SKcREHqWJh8nFTPPan94xas5YR0dHLC0tsbW1xd3dHYCLFy+m2k6j0TB//nzD3/fOnTuzfPlyIiMjsbe3p1KlSjRv3pw9e/bg7+9PkSJFsLW1NXTjpiU4OJizZ88SGhpqaBX+8ssvVK5cmePHj1O7dm2aNWtGSEgIo0aNIiQkhJYtW3Lx4kUOHDhA69atCQkJYfTo0Zn9KRklw7eHHDt2jMOHD3P37l0A3N3dqV+/PnXq1Mny4ITISqFRcYYKQcm0ClzFC5fm46D5OLh3Ec7/pu+Wjbr8X6EDS3uo0FY/cKhUM/1tK0IUILa2tikaQW5ubnh7e2Nvb59i3b1794w+5oULF/D09EzRdVqpUiWcnJy4cOECtWvXpmnTpvz8889otVr27t1Lq1atcHd3JyQkhGrVqnH16tUXdutmFaN/2+/du0enTp04ePAgJUuWxM3NDYDIyEg+/vhjGjZsyPr163F1dc22YIV4Fc+X0wNQq1R4Oz/TUnStAK7joNlYuPcPnP0Vzv2qn5Pz7yD9Yu8GVTpDdX99sYM8OkBB5AALW33LzlSvnZWHe24ydJVKleY6nU6Xpa/bpEkTHj16xKlTp9i3bx+TJ0/G3d2dKVOmUL16dYoVK0bZsmWz9DWfZ3SiHDJkCFqtlgsXLlC+fPkUz126dIn33nuPoUOHsm7duiwPUois8NJyes9SqcCtsn55YwLcOg5/r9WPkn0cCUfm6hfXylCjh76laS//JIrnqFRGdX+amqWlpWEQTk6qWLEiN2/e5ObNm4ZW5T///EN0dDSVKlUCwMnJiWrVqvHDDz9gYWFBhQoVcHV1xd/fn82bN2f79UnIQKLcvn07+/btS5UkAcqXL8/333+f7c1fIZJFxDwhNCqOUs52Gar4k6FyeslUKvCso19aB8LVXfoRspe2wr3zsOMz2DkByvlBzV76ikDSNSvyEG9vb44ePUpYWBj29vZZ3ip8kRYtWlC1alV69uzJrFmzSEpKYsiQITRt2hRfX1/Dds2aNWPOnDl07twZgCJFilCxYkWCgoKYO3dutsdpdGUeKysrYmNjX/j8o0ePsLKSG7hF9gs6Hk7DKbvpsegoDafsJuh4xkZZezjaUN+naCZL6lnoixZ0XQajLkPbGVDcFxStPnGu6Q4zK8HOifC/axk/vhAmMGrUKNRqNZUqVcLFxSXH7lxQqVRs2rSJwoUL06RJE1q0aEHp0qUJCgpKsV3Tpk3RarUpGmPNmjVLtS7b4lQU4+p6DR06lC1btjBz5kzeeOMNHBz0w+ZjY2MJDg4mICCAt956izlz5mRrwK8qNjYWR0dHYmJiDO9B5B0RMU9oOGV3quuMB8Y2x9nWnK1bt9KmTZtU106y3b2L8NdyfUszPuq/9aWaQq2+UOEtME89PF7kDxqNhh07dlCqVClKly6NtbW1qUMS/3r69CmhoaGUKlUq1edibD4wun9oxowZ6HQ6unXrRlJSkuGemMTERMzNzenfvz/Tpk3L5FsRwjhpj1xVCIuKx7mkCf/xca0Aft/AGxPh8jY4tQyuButnPQndC3au8FovqNUPnKTUoxB5idGJ0srKivnz5/Ptt99y4sQJQ+UEd3d3atWqJa0zkSOMGrlqSuaWUOlt/fLwhr6VeWo5PL4L+6fDgZn6ezLrDNTPpSkjZoXI9TI84sDBwYHXX389O2IR4qXSG7mq0WhMHV5Khb3g9c+h6Rj99cvjP0Hovv+qADmXgzrvQ/XuYGX/8uMJIUzC6ET5/fffG7XdsGHDMh2MEMbI1MhVU1JbQKX2+uX+ZX3CPL1KX9Bg6ygI/krfLVv3A3AqaepohRDPMTpRzpw5M8Xjmzdv4uHhgbn5f4dQqVQZTpRz585l6tSp3L17l+rVqzNnzpx0q/zMmjWL+fPnEx4ejrOzM507dyYwMFAunhcwmZ4I2tRcykGb7/QtzTNr4NiP8L+rcPgHOPLvTCf1PwLP2qaOVGSSkeMjRQ7Jis/D6EQZGhqa4nGhQoXYu3cvpUuXzvSLBwUFERAQwIIFC6hbty6zZs3Cz8+PS5cupVnhZ9WqVYwdO5bFixfToEEDLl++TN++fVGpVMyYMSPTcQiR46wdoO77UHuA/r7MI3Pheoi+3uw/m8CzHjT4SF+Y3UzmV88Lkm/Yj4+Px8YmD/4Tl0/Fx+uL0r/KSHiT3hU9Y8YMBg4cSL9+/QBYsGABW7ZsYfHixYwdOzbV9ocOHaJhw4b06NED0N8k2717d44ePZqjcQuRZczMoFwr/RJ5Hg7P05fJu3kEgo5A0TLQcDhU85eJpnM5RVFwcHAw1Dq1tbXNs/Mv5geKohAfH8+9e/dwcnJCrVZn+lgmS5SJiYmcPHmScePGGdaZmZnRokULDh8+nOY+DRo0YMWKFRw7dow6depw/fp1tm7dSq9evV74OgkJCSQkJBgeJxdN0Gg0uW/wh3glyZ9nnv1ci5SDtrOgyRjMTvyE2cklqP53FX7/CGXPZHR1B6Or2VtfoF3kKsnnXOHChQEMdwUI03NwcKBo0aJp/l0w9m+FyRJlVFQUWq3WUFw9mZubW5pTvAD06NGDqKgoGjVqhKIoJCUlMWjQID799NMXvk5gYCCTJk1KtX7Hjh3Y2uaSWwpEltq5c6epQ8gCtTAvXwmvqBB87m/D5lEE6l0T0O75lvNFWnPAtiUOdrY4SSMzV9m1axegH6/xKi0YkTW0Wm261yiTu2VfxuhE+Xz5OpVKxePHj1Otz877KUNCQpg8eTLz5s2jbt26XL16leHDh/PVV18xfvz4NPcZN24cAQEBhsexsbF4enrSqlUrufczn9FoNOzcuZOWLVvmfGWebNMJkhJIOrcO9aHZWD4Mpeb99fgoW/lF64ddqxG0r1/F1EEWePnz3Mv/0ivL+iyjE6WTk1OK/nZFUahZs2aKxyqVyugK9M7OzqjV6lRdFJGRkS+c5HP8+PH06tWLAQMGAFC1alXi4uJ4//33+eyzzzBLY9CDlZVVmjVoLSws5ITOp/LdZ2thAbX7EVGmC1OmBTJYvYkKZjf50Hwjj4K38/TxBxRqNhxsi5g60gIv3517+Zyxn5XRiXLPnj2ZDiYtlpaW1KpVi+DgYDp06ACATqcjODiYDz/8MM194uPjUyXD5O4NGZIt8rvQB0/ZpG3A79p6tDI7yTDzDVQ2uwHHZsHpn6DeYKj/Idg4mTpUIfIVoxOlVquladOmWdrvHhAQQJ8+ffD19aVOnTrMmjWLuLg4wyjY3r17U7x4cQIDAwFo164dM2bMoGbNmoau1/Hjx9OuXTu5HiDyvf/K95mxXVeb7Ym+tFafZE6xHVjcPwf7psLRhfrbSuoNAqtCpg5ZiHzB6EQ5YMAAoqOjad26Ne3bt+fNN9985Wt8/v7+3L9/nwkTJnD37l1q1KjBtm3bDAN8wsPDU7QgP//8c1QqFZ9//jm3b9/GxcWFdu3a8c0337xSHELkBanL95nRvEM/LHwnwIU/YM9kuH8B9nytL17QZBT49gcLKcYhxKswepotgL///pvff/+d33//nbNnz9KoUSPefvtt2rdvT8mSeaP0lkyzlX9pNBrTTbOVgyJinqRdvk+nhfO/QUigvtoPgEMJaD4OqnWTyaSzUUE59/IbY/NBhkp+VKtWjc8//5xjx45x7do1OnXqxJ9//kn58uWpUaMGEyZM4MSJE68cvBDixV448bSZGqp2hiFH4e05UKgYxN6CTUNhQUO4tA3kWr4QGZbp2ljFihVj0KBBbN26lfv37/P5558TFhZG69atmTx5clbGKITICLU5vNYbhp2CVl+DTWG4fxFW+8PSt+D2SVNHKESekiV9Mfb29nTu3JnOnTuj1Wp58OBBVhxWCPEqLGz0A3tq9tLPg3lkPtw4AIteh6pd9JNMyyTSQrxUphJlcHAwwcHB3Lt3D51OZ1ivUqn4+eefcXFxybIAhRCvyMYJWk7SF2Df841+1pKz6/QDgOp/CI1GyAhZIdKR4a7XSZMm0apVK4KDg4mKiuLhw4eGRVqSQuRiTp7wzgJ4PwS8GkHSU9g/DebU0s+P+cw/vUKI/2S4RblgwQKWLl2abiFyIUQuVqwG9N0MF7fAjs/hYShsHKyfUPrN76CEr6kjFCJXyXCLMjExkQYNGmRHLEKInKJSQcW3YOhRaDFJPyPJ7ZPw0xuwcQg8vm/YNCLmCYeuRRER88SEAQthOhlOlAMGDGDVqlXZEYsQIqeZW+mvUX50Eqrr53nl9Ep9d+zRH1l79DoNp+ymx6KjNJyym6Dj4SYNVwhTMKrr9dnZN3Q6HQsXLmTXrl1Uq1Yt1c21M2bMyNoIhRDZr5A7vDMffPvB1lEQcQb+HE0VnRfVeY+/KItOgU83nKNJOZfU93AKkY8ZlSj/+uuvFI9r1KgBwLlz51Ksl9m8hcjjPOvAwD1wcimanZOolHiD9ZZfsEbbnG+TuhGj2BMWFS+JUhQoRiXKrJ45RAiRi5mpoXZ/HpRoxf55Q+is3kcP8920Up9gctK7eBdtbuoIhchRma7MI4TI39w8PNG+PZduiRO4rCuOsyqWGRbz8PijJzwINXV4QuQYoxLloEGDuHXrllEHDAoKYuXKla8UlBAid/CvXZKZY4bw4N1gYhuOA7UVXNsN8+rDgVmgTTJ1iEJkO6O6Xl1cXKhcuTINGzakXbt2+Pr6UqxYMaytrXn48CH//PMPBw4cYM2aNRQrVoyFCxdmd9xCiBzi4WijvyZZbiy81gX+GA5h+2HXRDi/AdrPBfeqpg5TiGxjVIvyq6++4vLlyzRs2JB58+ZRr149SpYsiaurK+XLl6d3795cv36dhQsXcuTIEapVq5bdcQshTKGoD/T5Q58crR31o2MXNoPdX0NSgqmjEyJbGF2Zx83Njc8++4zPPvuMhw8fEh4ezpMnT3B2dsbHx0dGvApRUKhUUPNdKNMCtoyEi5th31R97dgO86B4LVNHKESWylRR9MKFC1O4cOGsjkUIkZcUcoduK+H8Rv29l/cvwk8t9QUMmo7RFzMQIh+QUa9CiFdTuYN+sugqnUDRwv7p+u7YiDOmjkyILCGJUgjx6uyKQufF0HU52LnAvX/0817unSojY0WeJ4lSCJF1Kr2tb11WfBt0SbDna1jsB1FXTR2ZEJkmiVIIkbXsikLXX+CdhWDlCLdPwI+N4cQSUBSZjUTkOZkazCOEEOlSqaC6P3g31E/bFboXNo/g9rGNdLjpz33FETMVBHasin/tkqaOVoh0GZUoa9asafTtH6dOnXqlgIQQ+YhjCei1EY7MQwmeRPF7IWy1PMUozWD26qrLbCQiTzAqUXbo0MHw/dOnT5k3bx6VKlWifv36ABw5coTz588zZMiQbAlSCJGHmZlBgw85Y1kD698HUcHsJsssv2VRUhumJvnLbCQi1zMqUU6cONHw/YABAxg2bBhfffVVqm1u3ryZtdEJIfINt7K1eF3zFWPUq+hrvoOB5lupb/YPrmYrgKKmDk+IF8rwYJ5169bRu3fvVOvfffdd1q9fnyVBCSHyHw9HG77oWIuvtP3onziSB4o9VczCcF3lB2fWmDo8IV4ow4N5bGxsOHjwIGXLlk2x/uDBg1hbW2dZYEKI/Me/dkmalHMhLKouGutesGsYhO6D3z6A63uhzVSwsjd1mEKkkOFEOWLECAYPHsypU6eoU6cOAEePHmXx4sWMHz8+ywMUQuQvhtlIKKof6LN/OoQEwplV+ltJuiwFt8omjlKI/2Q4UY4dO5bSpUsze/ZsVqxYAUDFihVZsmQJXbt2zfIAhRD5mJkamo4GrwawfgBEXYZFb0Db6VCzp6mjEwLI5H2UXbt2laQohMg63o1g0AHY8D5cC4ZNQyD8ELSZBhYyIlaYllTmEULkDnbO0PNXaP45qMzgrxXwUwv437UUm0llH5HTMtyi1Gq1zJw5k7Vr1xIeHk5iYmKK5x88eJBlwQkhChgzM2j6CXjWgfX9IfIcLGwO7yyACm0IOh7OuA1n0SlIZR+RYzLcopw0aRIzZszA39+fmJgYAgIC6NixI2ZmZnzxxRfZEKIQosAp3RQ+2A+edSEhBtZ05/GW8Xy+4TQ6Rb+JToFPN5yTlqXIdhlOlCtXrmTRokWMHDkSc3Nzunfvzk8//cSECRM4cuRIdsQohCiIHDygz2aoOwgA++Pfs9j8W5x4ZNhEqyiERcWbKkJRQGQ4Ud69e5eqVasCYG9vT0xMDABvvfUWW7ZsydrohBAFm7klvPktdPoZnbkNjdXn+MPycyqpwgBQq1R4O9uaNkaR72U4UZYoUYKIiAgAfHx82LFjBwDHjx/Hysoqa6MTQgiAqp0xG7CLx7Yl8DS7z3rLL3hHfZDJHatInViR7TKcKN955x2Cg4MB+Oijjxg/fjxly5ald+/evPfee1keoBBCAOBeBfsP9/PUqzk2qkRmWszF/+FC0GlNHZnI5zI86nXKlCmG7/39/SlZsiSHDx+mbNmytGvXLkuDE0KIFGyLYN1nPez+Cg7MhENz4N4F6PQz2DiZOjqRT73yxM3169c3TLclhBDZzkwNLb4A96qwcShc3QWLXofua8ClnKmjE/lQpgoOLF++nIYNG1KsWDFu3LgBwKxZs9i0aVOWBieEEC9UpRO8tw0cSsCDa/riBFeDTR2VyIcynCjnz59PQEAAbdq0ITo6Gq1Wf33AycmJWbNmZTiAuXPn4u3tjbW1NXXr1uXYsWPpbh8dHc3QoUPx8PDAysqKcuXKsXXr1gy/rhAiHyhWA94PAc96+vstV3aGIwtAUUwdmchHMpwo58yZw6JFi/jss89Qq9WG9b6+vpw9ezZDxwoKCiIgIICJEydy6tQpqlevjp+fH/fu3Utz+8TERFq2bElYWBi//vorly5dYtGiRRQvXjyjb0MIkV/Yu0Cf36FGT1B0sG0MbB4BWg0gJe/Eq8vwNcrQ0FBq1qyZar2VlRVxcXEZOtaMGTMYOHAg/fr1A2DBggVs2bKFxYsXM3bs2FTbL168mAcPHnDo0CEsLCwA8Pb2zuhbEELkN+ZW0H4uuFSAnRPg5FJ4EMqGst8w6o9wKXknXkmGE2WpUqU4ffo0Xl5eKdZv27aNihUrGn2cxMRETp48ybhx4wzrzMzMaNGiBYcPH05zn99//5369eszdOhQNm3ahIuLCz169GDMmDEpWrfPSkhIICEhwfA4NjYWAI1Gg0ajMTpekfslf57yuRZgdQajKlwa9W/vowrdS7VrXSnBJ4Tjhk6BcRvOUr9UYTwcs3aSeTn38iZjP68MJ8qAgACGDh3K06dPURSFY8eOsXr1agIDA/npp5+MPk5UVBRarRY3N7cU693c3Lh48WKa+1y/fp3du3fTs2dPtm7dytWrVxkyZAgajYaJEyemuU9gYCCTJk1KtX7Hjh3Y2kpFj/xo586dpg5BmJhD6XHUujqDMtxho+V43k8M4IRSAZ0Ca7fuoaxj9lzDlHMvb4mPN678oUpRMn7Ve+XKlXzxxRdcu6af/qZYsWJMmjSJ/v37G32MO3fuULx4cQ4dOpTi9pLRo0ezd+9ejh49mmqfcuXK8fTpU0JDQw0tyBkzZjB16lRDtaDnpdWi9PT0JCoqCgcHB6PjFbmfRqNh586dtGzZ0tA1LwquyDs3ePhzZ6qahZKgmDNKM4gtSgNCRjbJlhalnHt5T2xsLM7OzsTExKSbDzJ1H2XPnj3p2bMn8fHxPH78GFdX1wwfw9nZGbVaTWRkZIr1kZGRuLu7p7mPh4cHFhYWKbpZK1asyN27d0lMTMTS0jLVPlZWVmmW1rOwsJATOp+Sz1YAlPAqw5E2QdzZMhQ/9XHmWP7AgPJqShZtAypVtrymnHt5i7Gf1StN3Gxra5upJAlgaWlJrVq1DOXwAHQ6HcHBwS8sYNCwYUOuXr2KTqczrLt8+TIeHh5pJkkhRMHWuV55qgVs5E5FfXnN6pdmwx/DDSNihTBGhhNlZGQkvXr1olixYpibm6NWq1MsGREQEMCiRYtYtmwZFy5cYPDgwcTFxRlGwfbu3TvFYJ/Bgwfz4MEDhg8fzuXLl9myZQuTJ09m6NChGX0bQogCwqOwPcX8Z8KbU0FlBqeWwepukPDY1KGJPCLDXa99+/YlPDyc8ePH4+HhgeoVujD8/f25f/8+EyZM4O7du9SoUYNt27YZBviEh4djZvZfLvf09GT79u18/PHHVKtWjeLFizN8+HDGjBmT6RiEEAVE3ffByRPW9dOXvVvaBnqsg0JuL99XFGgZHsxTqFAh9u/fT40aNbIppOwVGxuLo6PjSy/eirxHo9GwdetW2rRpI9eJxIvdOgmrukJ8FDiVhJ7rX7lGrJx7eZOx+SDDXa+enp5kYqCsEELkDiVqQf8dUKQ0RIfD4lZwU186U6r4iLRkOFHOmjWLsWPHEhYWlg3hCCFEDijqA/13QvFa8OQhLHubfZuX03DKbnosOkrDKbsJOh5u6ihFLmHUNcrChQunuBYZFxeHj48Ptra2qboZHjx4kLURCiFEdrBzhj5/wLq+cGUHDY4Po5PZANZpm6FT4NMN52hSzgUPRxtTRypMzKhEmZlZQYQQIteztINuq7i36gNcr61nqsVCXIhmnrY9WgXCouIlUQrjEmWfPn2yOw4hhDANtQXadj8wb5qGIea/M9piLUVVjwjUvou3s5S5FK9YcEAIIfIDDydbirb/hq+SegHQ3/xPdpVejYd9poqXiXxGEqUQQgD+tUsy4JOpXGk4A8XMHO/bm2F1d0g0rnC2yL8kUQohxL88HG0o27I/qm6rwdwGru6EFR3hSbSpQxMmJIlSCCGeV64V9N4IVo4QfhiWvgWP75k6KmEikiiFECItJetBvy1g5wqRZ2Gxn75AgShwjLpS3bFjR6MPuGHDhkwHI4QQuYp7VXhvG/zSAR5ch8WtofcmcC5r6shEDjKqReno6GhYHBwcCA4O5sSJE4bnT548SXBwMI6OjtkWqBBCmERRH+i/HZzLQextfbKM+NvUUYkcZFSLcsmSJYbvx4wZQ9euXVmwYIFhWi2tVsuQIUOkyLgQIn9yKAb9/oTl78Ddv/XXLN/9FTzrmDoykQMyfI1y8eLFjBo1KsXck2q1moCAABYvXpylwQkhRK6RXPLOsx4kxOi7Y6/vBSAi5ilXYlRExDw1bYwiW2Q4USYlJXHx4sVU6y9evIhOp8uSoIQQIleycYJeG6B0c9DEwaqu7N2ykmbT9/HDP2qaTd8nxdTzoQyXnejXrx/9+/fn2rVr1Kmj73Y4evQoU6ZMoV+/flkeoBBC5CqWdtB9jb6Y+uU/qX/sI1qqPmK7UkeKqedTGU6U06ZNw93dnenTpxMREQGAh4cHn3zyCSNHjszyAIUQItexsAb/5UT90gfnG1uYa/E9IzWD2KRrhFZRpJh6PpPhrlczMzNGjx7N7du3iY6OJjo6mtu3bzN69OgU1y2FECJfU1ug6bCQX7VNMFfpmGkxny7qENQqlRRTz2cyVXAgKSmJXbt2sXr1asM8lXfu3OHx48dZGpwQQuRmHoXt0babwwptC8xUClMtFhJU6x9pTeYzGe56vXHjBq1btyY8PJyEhARatmxJoUKF+Pbbb0lISGDBggXZEacQQuRK/nW8CS+1hJPLB1Ardju+576C4rZQf4ipQxNZJMMtyuHDh+Pr68vDhw+xsfnvv6Z33nmH4ODgLA1OCCHyAg8nG26V7oG2/jD9iu3j4OBs0wYlskyGW5T79+/n0KFDWFpapljv7e3N7du3sywwIYTIU1QqdM3Ho7a0gb3fws4JoEuCxjLIMa/LcItSp9Oh1WpTrb916xaFChXKkqCEECJPUqmg+afQ/DP94+AvYe93po1JvLIMJ8pWrVoxa9Ysw2OVSsXjx4+ZOHEibdq0ycrYhBAib2o6Gt6YoP9+zzewZzIoimljEpmW4a7X6dOn4+fnR6VKlXj69Ck9evTgypUrODs7s3r16uyIUQgh8p7GI8HMAnaO13fFKgo0/5SI2KeERsVRytlORsfmERlOlCVKlODMmTMEBQVx5swZHj9+TP/+/enZs2eKwT1CCFHgNRwGZmrY/ins+47zdx7S7nwzdIoKMxUEdqyKf+2Spo5SvESGEyWAubk5PXv2pGfPnlkdjxBC5C/1hwIq2D6OylcXMVIdydQkf3SKSsrd5REZvkapVqtp3rw5Dx48SLE+MjJSKvMIIURa6g/huu94AIaa/85o8yBAMZS7E7lbhhOloigkJCTg6+vL+fPnUz0nhBAiNZvGQ/lC0weAIea/84l5EGoVUu4uD8hwolSpVKxfv5527dpRv359Nm3alOI5IYQQqXk42lCxwygm/Zssh5r/zh+V9uDhYG3iyMTLZKpFqVarmT17NtOmTcPf35+vv/5aWpNCCPES/rVL8v7obwn11d86UunaItj9ldw6kstlajBPsvfff5+yZcvSpUsX9u3bl1UxCSFEvuXhaANvjQRnW9g2FvZPB5WZvkiB9MrlShluUXp5eaUYtNO8eXOOHDnCzZs3szQwIYTI1+oNBr9A/ff7pkLIFNPGI14owy3K0NDQVOvKlCnDX3/9RWRkZJYEJYQQBUL9IaDoYMdnsHeK/p7LpqNNHZV4Tqbmo0yLtbU1Xl5eWXU4IYQoGBp8CC2/1H+/5xtubPySiJgnpo1JpGBUi7JIkSJcvnwZZ2dnChcunO7o1ufvrxRCCPESDYdz5uYDql+chdfp6Uw5cYdS7T+Vqj25hFGJcubMmYaZQZ4tiC6EEOLVRcQ84Z0zdRhs1pVPLNYy1nw132xSE1FumlTtyQWMSpR9+vRJ83shhBCvLjQqDp0Cc7UdsFAlMcJ8A5+ZryB0f2l4ayQRMU+kkLoJGZUoY2NjjT6gg4NDpoMRQoiCqJSzHWYq0CkwK6kTanR8ZL6RUie+5MRTHV1PVkKnIIXUTcSowTxOTk4ULlw43SV5m8yYO3cu3t7eWFtbU7duXY4dO2bUfmvWrEGlUtGhQ4dMva4QQuQGHo42BHasilqlAlTM0nblQul+APie+5ouZnsAfSL9dMM5GeyTw4xqUe7ZsyfbAggKCiIgIIAFCxZQt25dZs2ahZ+fH5cuXcLV1fWF+4WFhTFq1CgaN26cbbEJIURO8a9dkiblXAiLisfb2RYPhzbcWaui2IXFBJr/RJKiZr2uiaGQunTB5hyjEmXTpk2zLYAZM2YwcOBA+vXT//e0YMECtmzZwuLFixk7dmya+2i1Wnr27MmkSZPYv38/0dHR2RafEELkFA9HmxQJUOX3DcvO3aKPegdTLX5Eo1GzRWkkhdRzWKZL2MXHxxMeHk5iYmKK9dWqVTP6GImJiZw8eZJx48YZ1pmZmdGiRQsOHz78wv2+/PJLXF1d6d+/P/v370/3NRISEkhISDA8Tr7eqtFo0Gg0Rscqcr/kz1M+V5HTsuvcc7azwKLNt6zaqqWHOpgZFvPpUtMbZ1tzOc+zgLE/wwwnyvv379OvXz/+/PPPNJ/XarVGHysqKgqtVoubm1uK9W5ubly8eDHNfQ4cOMDPP//M6dOnjXqNwMBAJk2alGr9jh07sLWV/8ryo507d5o6BFFAZce5ZwdEV+nF2Zsaqj7eR8PTYzkec4u7TrWy/LUKmvh44+YCzXCiHDFiBNHR0Rw9epRmzZrx22+/ERkZyddff8306dMzHGhGPHr0iF69erFo0SKcnZ2N2mfcuHEEBAQYHsfGxuLp6UmrVq1khG4+o9Fo2LlzJy1btsTCwsLU4YgCJEfOPd2b6P4Yitm5X6lzYx7aOstRyrTIntcqIIy9oyPDiXL37t1s2rQJX19fzMzM8PLyomXLljg4OBAYGEjbtm2NPpazszNqtTpVjdjIyEjc3d1TbX/t2jXCwsJo166dYZ1Op9O/EXNzLl26hI+PT4p9rKyssLKySnUsCwsL+WOaT8lnK0wle889C3jnR9AlofpnI+a/9oEeQeDTPJteL/8z9rPKcK3XuLg4w2jUwoULc//+fQCqVq3KqVOnMnQsS0tLatWqRXBwsGGdTqcjODiY+vXrp9q+QoUKnD17ltOnTxuWt99+m+bNm3P69Gk8PT0z+naEECLvUJtDp5+gfFvQJsDq7hB2wNRR5XsZblGWL1+eS5cu4e3tTfXq1fnxxx/x9vZmwYIFeHh4ZDiAgIAA+vTpg6+vL3Xq1GHWrFnExcUZRsH27t2b4sWLExgYiLW1NVWqVEmxv5OTE0Cq9UIIkS+pLaDLEgh6F67sgJVdodcGKFnP1JHlWxlOlMOHDyciIgKAiRMn0rp1a1auXImlpSVLly7NcAD+/v7cv3+fCRMmcPfuXWrUqMG2bdsMA3zCw8MxM8uySU6EECLvM7eCrsthtT9cD0G3ohPn3vgFlwoN5P7KbKBSFEV5lQPEx8dz8eJFSpYsafQAG1OKjY3F0dGRmJgYGcyTz2g0GrZu3UqbNm3kGqXIUSY79xLjubegHa4PThCj2NJT8zm93mknJe6MZGw+eOWmmq2tLa+99lqeSJJCCJGfRDxR8XrEEE7oyuGoiucXi8ks/W2rlLjLYhnuelUUhV9//ZU9e/Zw7949w6jTZBs2bMiy4IQQQrxYaFQcjxVr+iWOZrnlZGqYXecXi2+4de01PF6rY+rw8o0MtyhHjBhBr169CA0Nxd7eHkdHxxSLEEKInJE868gjbOmdOJbzOi9cVLFUC34X/nfN1OHlGxluUS5fvpwNGzbQpk2b7IhHCCGEkZJnHfl0wzliFXv6aD5ll/NUnB5dhWVvQ7+tUNjL1GHmeRlOlI6OjpQuXTo7YhFCCJFBz8864qR+HZa2hajLsKwdkZ03cC3BSSZ9fgUZ7nr94osvmDRpEk+eyMViIYTIDTwcbajvU1SfCO1doffvUKQ0RN8gfmFbhi/aRsMpuwk6Hm7qUPOkDCfKrl278vDhQ1xdXalatSqvvfZaikUIIYSJOXgQ2XEdNxUXSpndZbXlNxRRYmTS50zKcNdrnz59OHnyJO+++y5ubm6oVKrsiEsIIcQruJbgxOjEz1hr+SVlzO6wwnIy3RM/k0mfMyHDiXLLli1s376dRo0aZUc8QgghskApZzvu4Er3xM9Za/klFcxussIykKL2TU0dWp6T4a5XT09PqWgjhBC5XPKI2Ft40CPxM6IUByqb3cB9Uw94GmPq8PKUDCfK6dOnM3r0aMLCwrIhHCGEEFnFv3ZJDoxtztcDOqH03gQ2ReDOX7CiEzw1bi5GkYmu13fffZf4+Hh8fHywtbVNVdfwwYMHWRacEEKIV+PhaPPvNcmi0HsTLGsHt46T+EsnTjX5CS8PV7lm+RIZTpSzZs3KhjCEEEJkO49q0HsjiYvbYXnnGKzy5w3NaCZ29JVC6unIUKLUaDTs3buX8ePHU6pUqeyKSQghRDaJsKvA4PhP+MUikHpmF1hoPo33N4ymSTkXaVm+QIauUVpYWLB+/frsikUIIUQ2C42K47SuDH0Tx/BYsaaR+jwLzKdxI/KhqUPLtTI8mKdDhw5s3LgxG0IRQgiR3ZILqZ9SytE3cTRxihVN1GepcWgoJCWYOrxcKcPXKMuWLcuXX37JwYMHqVWrFnZ2dimeHzZsWJYFJ4QQIms9W0j9hFKBAZrRLLeeinXYbljbG7r+AuZWpg4zV8lwovz5559xcnLi5MmTnDx5MsVzKpVKEqUQQuRyKQupv475/2rDqq5weRus6wtdloG5panDzDUynChDQ0OzIw4hhBA56L/bRgDHptB9NazuDpe2/pssl0qy/FeGr1E+S1EUFEXJqliEEEKYis/r0G0VqK3g0hb4tR9oNaaOKlfIVKL85ZdfqFq1KjY2NtjY2FCtWjWWL1+e1bEJIYTISWXegO7/JsuLm/nf0h5EPJAKPhlOlDNmzGDw4MG0adOGtWvXsnbtWlq3bs2gQYOYOXNmdsQohBAip5Rpwd7XZpGgWFD05g7+ntmRdUevmToqk8rwNco5c+Ywf/58evfubVj39ttvU7lyZb744gs+/vjjLA1QCCFEzomIeUK/A440VgWw0GIGfurj7Nj8PhFlN+BRxNHU4ZlEhluUERERNGjQINX6Bg0aEBERkSVBCSGEMI3QqDh0CuzVVWegJoAExYJW6hNYbegPSYlExDzh0LWoAjUBdIYTZZkyZVi7dm2q9UFBQZQtWzZLghJCCGEayQUJAPY9kyyL3NrJ7YVdaD5lOz0WHaXhlN0EHQ83bbA5JMNdr5MmTcLf3599+/bRsGFDAA4ePEhwcHCaCVQIIUTe8WxBAq2icFCpweG6P9D05HCK3wthnnk0gzUjSFAs+XTDuQJRIzbDibJTp04cPXqUmTNnGkrZVaxYkWPHjlGzZs2sjk8IIUQOS1mQwBYPRxvOO9hQeucAXlefZhHTGagZSYJiSVhUvCTKtNSqVYsVK1ZkdSxCCCFyiRQFCYAiVVvx3tbR/GQxlSbqs/zMVAYljcLb2daEUeaMVyo4IIQQomDwcLShwzv+vKcZa5h1ZLfHD3hYJ5k6tGxndKI0MzNDrVanu5ibZ6qBKoQQIg/wr12SWWMGE/rmCnSWhXB9cBKWvwNPok0dWrYyOrP99ttvL3zu8OHDfP/99+h0uiwJSgghRO7k4WiDR72W4Pm7PkneOg6/vM3d9qu5HmdFKWe7fHfN0uhE2b59+1TrLl26xNixY/njjz/o2bMnX375ZZYGJ4QQIpcq/hr03Qy/dICIM0TP82NY4qc8UDkS2LEq/rVLmjrCLJOpa5R37txh4MCBVK1alaSkJE6fPs2yZcvw8vLK6viEEELkVu5Vud95PZGKExXMbrLWchJuyv/4dMO5fFWQIEOJMiYmhjFjxlCmTBnOnz9PcHAwf/zxB1WqVMmu+IQQQuRiV5QSdEmcyC3FmdJmd1lnNYni3CUsKt7UoWUZoxPld999R+nSpdm8eTOrV6/m0KFDNG7cODtjE0IIkcuVcrbjFm50SZjIdZ07JVRRrLOcRBnyT9Ueo69Rjh07FhsbG8qUKcOyZctYtmxZmttt2LAhy4ITQgiRuz1bycc/cQLLLQOpYHYTTVAHojqtxrl86trgeY3RibJ3796oVKrsjEUIIUQe9GwlnyNhFXmydyA1E69iveod9tSbS/M3O5s6xFdidKJcunRpNoYhhBAiL0u+JaTnTxFYK5+y0GI6jdTnaXBkEA/c1CT4tCY0Ki5P3j4iFQKEEEJkieQpuuKx5j3NaL7nB1qrj2Pxe3/GavqzVtsMMxV57vYRKWEnhBAiSzw7RVciFgzVDGNtUjPM0PKdxUI+UP+BTlHy3O0juSJRzp07F29vb6ytralbty7Hjh174baLFi2icePGFC5cmMKFC9OiRYt0txdCCJEzkgf2qJPHs6jMuVpvMvOT2gEwzmI1n5mvRKdo89TtIybveg0KCiIgIIAFCxZQt25dZs2ahZ+fH5cuXcLV1TXV9iEhIXTv3p0GDRpgbW3Nt99+S6tWrTh//jzFixc3wTsQQgiR7PkpugAaHuzO/xQHPrdYyUDzrbioYvAu3MTEkRrP5C3KGTNmMHDgQPr160elSpVYsGABtra2LF68OM3tV65cyZAhQ6hRowYVKlTgp59+QqfTERwcnMORCyGESIuHow31fYoapuoK7FiVJbq3CEgchEZR00F9EI/NveFprKlDNYpJW5SJiYmcPHmScePGGdaZmZnRokULDh8+bNQx4uPj0Wg0FClSJM3nExISSEhIMDyOjdV/MBqNBo1G8wrRi9wm+fOUz1XkNDn30texhgf1SxUm/EEtHj5qhMu2D1BdD0FZ0oYk/9VQyN0kcRn7eZk0UUZFRaHVanFzc0ux3s3NjYsXLxp1jDFjxlCsWDFatGiR5vOBgYFMmjQp1fodO3Zga5v/JxwtiHbu3GnqEEQBJefeyx0GHEuPpt616VhHniVhfjPWuY3CzKEYTlY5G0t8vHHXSU1+jfJVTJkyhTVr1hASEoK1tXWa24wbN46AgADD49jYWDw9PWnVqhUODg45FarIARqNhp07d9KyZUssLCxMHY4oQOTcy4SHbYld2hGH+HA63fySwZqPaft2F7rUKpFjIST3ML6MSROls7MzarWayMjIFOsjIyNxd0+/KT5t2jSmTJnCrl27qFat2gu3s7Kywsoq9b8pFhYWckLnU/LZClORc894EVYlaPvwM360mE5ts8sssZjCp388oHnFiTlWkMDYz8qkg3ksLS2pVatWioE4yQNz6tev/8L9vvvuO7766iu2bduGr69vToQqhBAiC4VGxfFAKcS7iZ/yh7Yeliot0yzmo931Deh0pg4vBZOPeg0ICGDRokUsW7aMCxcuMHjwYOLi4ujXrx+grzH77GCfb7/9lvHjx7N48WK8vb25e/cud+/e5fHjx6Z6C0IIITIouThBApYM03zIvKS3AShxdg6sfw80uacggckTpb+/P9OmTWPChAnUqFGD06dPs23bNsMAn/DwcCIiIgzbz58/n8TERDp37oyHh4dhmTZtmqneghBCiAx6tjiBghnTtd05Vu1LMLOA87/B0rbwKJKImCccuhZl0ko+KkVRFJO9ugnExsbi6OhITEyMDObJZzQaDVu3bqVNmzZynUjkKDn3Mi8i5omhOIGHow2EHYCgd+HJQ+Ks3ekW+xFndaWypUassfnA5C1KIYQQBdezxQkA8G4EA4JJKlwGu6d3WWsxibfMDqNTYNz6s5y5+TDHY5REKYQQIncp6sPJVusI0VbHRpXID5Zz+MR8DQo6Osw7RNDx8BwNRxKlEEKIXKdkMQ8GJH3CgqS3ABhq/js/W0yjkPI4x2cfkUQphBAi1/FwtOGbjtX5LqkHwxOH8FSx4HX1af6w/Jyy3MjR2UfydGUeIYQQ+Zd/7ZJUcC9Eh3lwNbEECyxm4mV2j98sJ/A0ygF8euZIHNKiFEIIkWtV9yzMlI5VuUgp2iV+zT5dNWxUiRT+cwjRh3/JkRikRSmEECJXe3aOy7O3XuPszq953ewUnX+3Y4J5eJbeMpIWSZRCCCFyveTbR3r+dASd4s/3vEMClny64RxNyrlka31Y6XoVQgiRJ4RGxaH7t0ROApYAaBUl2wf2SKIUQgiRJyTXh32WWqXC2zl75xaWRCmEECJPeLY+LOiT5OSOVbJ9Wi65RimEECLPeHZgj6E+bDaTRCmEECJP8XC0ybHJnUG6XoUQQoh0SaIUQggh0iGJUgghhEiHJEohhBAiHQVuMI+i6O9WjY2NNXEkIqtpNBri4+OJjY2VWeZFjpJzL29KzgPJeeFFClyifPToEQCenp4mjkQIIURu8OjRIxwdHV/4vEp5WSrNZ3Q6HXfu3KFQoUKoVKo0t6lduzbHjx/PsZiy+vVe5XiZ2dfYfbJquxc9Hxsbi6enJzdv3sTBweGlr5MbybmXPeeesdvKuVewzj1FUXj06BHFihXDzOzFVyILXIvSzMyMEiVKpLuNWq3O0ZM9q1/vVY6XmX2N3SertnvZ8w4ODnn2j5Wce9lz7hm7rZx7Be/cS68lmUwG86Rh6NChefr1XuV4mdnX2H2yaruc/nxykpx72bePMdvKuZd3Xy874y9wXa8i/4qNjcXR0ZGYmJg8+1+9yJvk3MvfpEUp8g0rKysmTpyIlZWVqUMRBYyce/mbtCiFEEKIdEiLUgghhEiHJEohhBAiHZIohRBCiHRIohRCCCHSIYlSCCGESIckSlFgxcfH4+XlxahRo0wdiihAoqOj8fX1pUaNGlSpUoVFixaZOiTxEgWuhJ0Qyb755hvq1atn6jBEAVOoUCH27duHra0tcXFxVKlShY4dO1K0aFFThyZeQFqUokC6cuUKFy9e5M033zR1KKKAUavV2NraApCQkICiKC+d5kmYliRKkevs27ePdu3aUaxYMVQqFRs3bky1zdy5c/H29sba2pq6dety7NixDL3GqFGjCAwMzKKIRX6SE+dfdHQ01atXp0SJEnzyySc4OztnUfQiO0iiFLlOXFwc1atXZ+7cuWk+HxQUREBAABMnTuTUqVNUr14dPz8/7t27Z9gm+frP88udO3fYtGkT5cqVo1y5cjn1lkQekt3nH4CTkxNnzpwhNDSUVatWERkZmSPvTWSSIkQuBii//fZbinV16tRRhg4danis1WqVYsWKKYGBgUYdc+zYsUqJEiUULy8vpWjRooqDg4MyadKkrAxb5BPZcf49b/Dgwcq6deteJUyRzaRFKfKUxMRETp48SYsWLQzrzMzMaNGiBYcPHzbqGIGBgdy8eZOwsDCmTZvGwIEDmTBhQnaFLPKRrDj/IiMjefToEQAxMTHs27eP8uXLZ0u8ImvIqFeRp0RFRaHVanFzc0ux3s3NjYsXL5ooKlFQZMX5d+PGDd5//33DIJ6PPvqIqlWrZke4IotIohQFWt++fU0dgihg6tSpw+nTp00dhsgA6XoVeYqzszNqtTrV4IfIyEjc3d1NFJUoKOT8K5gkUYo8xdLSklq1ahEcHGxYp9PpCA4Opn79+iaMTBQEcv4VTNL1KnKdx48fc/XqVcPj0NBQTp8+TZEiRShZsiQBAQH06dMHX19f6tSpw6xZs4iLi6Nfv34mjFrkF3L+iVRMPexWiOft2bNHAVItffr0MWwzZ84cpWTJkoqlpaVSp04d5ciRI6YLWOQrcv6J56kURWonCSGEEC8i1yiFEEKIdEiiFEIIIdIhiVIIIYRIhyRKIYQQIh2SKIUQQoh0SKIUQggh0iGJUgghhEiHJEohhBAiHZIohcjlQkJCUKlUREdH5/hrq1QqVCoVTk5O6W73xRdfUKNGjRSPk/edNWtWtsYoRHaTRClELtKsWTNGjBiRYl2DBg2IiIjA0dHRJDEtWbKEy5cvZ2ifUaNGERERQYkSJbIpKiFyjhRFFyKXs7S0NOkUTk5OTri6umZoH3t7e+zt7VGr1dkUlRA5R1qUQuQSffv2Ze/evcyePdvQbRkWFpaq63Xp0qU4OTmxefNmypcvj62tLZ07dyY+Pp5ly5bh7e1N4cKFGTZsGFqt1nD8hIQERo0aRfHixbGzs6Nu3bqEhIRkKtYpU6bg5uZGoUKF6N+/P0+fPs2Cn4AQuZMkSiFyidmzZ1O/fn0GDhxIREQEEREReHp6prltfHw833//PWvWrGHbtm2EhITwzjvvsHXrVrZu3cry5cv58ccf+fXXXw37fPjhhxw+fJg1a9bw999/06VLF1q3bs2VK1cyFOfatWv54osvmDx5MidOnMDDw4N58+a90nsXIjeTrlchcglHR0csLS2xtbV9aVerRqNh/vz5+Pj4ANC5c2eWL19OZGQk9vb2VKpUiebNm7Nnzx78/f0JDw9nyZIlhIeHU6xYMUB/HXHbtm0sWbKEyZMnGx3nrFmz6N+/P/379wfg66+/ZteuXdKqFPmWtCiFyINsbW0NSRLAzc0Nb29v7O3tU6y7d+8eAGfPnkWr1VKuXDnD9UN7e3v27t3LtWvXMvTaFy5coG7duinW1a9f/xXejRC5m7QohciDLCwsUjxWqVRprtPpdAA8fvwYtVrNyZMnUw2weTa5CiFSk0QpRC5iaWmZYgBOVqlZsyZarZZ79+7RuHHjVzpWxYoVOXr0KL179zasO3LkyKuGKESuJYlSiFzE29ubo0ePEhYWhr29PUWKFMmS45YrV46ePXvSu3dvpk+fTs2aNbl//z7BwcFUq1aNtm3bGn2s4cOH07dvX3x9fWnYsCErV67k/PnzlC5dOktiFSK3kWuUQuQio0aNQq1WU6lSJVxcXAgPD8+yYy9ZsoTevXszcuRIypcvT4cOHTh+/DglS5bM0HH8/f0ZP348o0ePplatWty4cYPBgwdnWZxC5DYqRVEUUwchhMidVCoVv/32Gx06dMjU/t7e3owYMSJVtSEh8hJpUQoh0tW9e/cMl6KbPHky9vb2WdoiFsJUpEUphHihq1evAqBWqylVqpTR+z148IAHDx4A4OLiYrI6tUJkBUmUQgghRDqk61UIIYRIhyRKIYQQIh2SKIUQQoh0SKIUQggh0iGJUgghhEiHJEohhBAiHZIohRBCiHRIohRCCCHSIYlSCCGESMf/AXbmF0HCkT4XAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tm = np.logspace(np.log10(to[0]), np.log10(to[-1]), 100)\n",
    "hm = w.headinside(tm)\n",
    "plt.semilogx(to, ho / H0, \".\", label=\"obs\")\n",
    "plt.semilogx(tm, hm[0] / H0, label=\"timflow\")\n",
    "plt.xlabel(\"time [d]\")\n",
    "plt.ylabel(\"Normalized head (h/H0)\")\n",
    "plt.title(\"Model results - multi-layer model\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of results\n",
    "Here, the `timflow` performance in analysing slug tests is checked. The solution in `timflow` is compared with the KGS analytical model (Hyder et al. 1994) implemented in AQTESOLV (Duffield, 2007). The parameters of `timflow` and AQTESOLV are similar, even though AQTESOLV only used one layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style type=\"text/css\">\n",
       "</style>\n",
       "<table id=\"T_f6cb8\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th class=\"blank level0\" >&nbsp;</th>\n",
       "      <th id=\"T_f6cb8_level0_col0\" class=\"col_heading level0 col0\" >k [m/d]</th>\n",
       "      <th id=\"T_f6cb8_level0_col1\" class=\"col_heading level0 col1\" >Ss [1/m]</th>\n",
       "      <th id=\"T_f6cb8_level0_col2\" class=\"col_heading level0 col2\" >RMSE [m]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th id=\"T_f6cb8_level0_row0\" class=\"row_heading level0 row0\" >timflow</th>\n",
       "      <td id=\"T_f6cb8_row0_col0\" class=\"data row0 col0\" >0.44</td>\n",
       "      <td id=\"T_f6cb8_row0_col1\" class=\"data row0 col1\" >4.61e-04</td>\n",
       "      <td id=\"T_f6cb8_row0_col2\" class=\"data row0 col2\" >0.006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th id=\"T_f6cb8_level0_row1\" class=\"row_heading level0 row1\" >AQTESOLV</th>\n",
       "      <td id=\"T_f6cb8_row1_col0\" class=\"data row1 col0\" >0.42</td>\n",
       "      <td id=\"T_f6cb8_row1_col1\" class=\"data row1 col1\" >5.70e-04</td>\n",
       "      <td id=\"T_f6cb8_row1_col2\" class=\"data row1 col2\" >-</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n"
      ],
      "text/plain": [
       "<pandas.io.formats.style.Styler at 0x15d60aea0>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = pd.DataFrame(\n",
    "    columns=[\"k [m/d]\", \"Ss [1/m]\", \"RMSE [m]\"],\n",
    "    index=[\"timflow\", \"AQTESOLV\"],\n",
    ")\n",
    "\n",
    "t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
    "t.loc[\"AQTESOLV\"] = [0.4211, 5.70e-4, \"-\"]\n",
    "\n",
    "t_formatted = t.style.format(\n",
    "    {\n",
    "        \"k [m/d]\": \"{:.2f}\",\n",
    "        \"Ss [1/m]\": \"{:.2e}\",\n",
    "        \"RMSE [m]\": lambda x: \"-\" if x == \"-\" else f\"{float(x):.3f}\",\n",
    "    }\n",
    ")\n",
    "t_formatted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "## References\n",
    "\n",
    "* Batu, V. (1998), Aquifer hydraulics: a comprehensive guide to hydrogeologic data analysis, John Wiley & Sons\n",
    "* Duffield, G.M. (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
    "* Hyder, Z., Butler Jr, J.J., McElwee, C.D. and Liu, W. (1994), Slug tests in partially penetrating wells, Water Resources Research 30, 2945–2957."
   ]
  },
  {
   "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.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}