{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Leaky Aquifer Recovery Test - Hardixveld Example\n",
    "**This test is taken from MLU examples.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import ttim as ttm\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = [5, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Introduction and Conceptual Model\n",
    "\n",
    "This test was taken in Hardinxveld-Giessendam, Netherlands, in 1981 to quantify the head-loss at each pumping well by clogging and to assess the transmissivity variation in the area.\n",
    "\n",
    "The hydrogeological conceptualization can be described as the following:\n",
    "* The first ten meters depth is an aquitard\n",
    "* Followed by the first aquifer from 10 to 37 m depth, this is also the test aquifer.\n",
    "* A new aquitard is present from 37 m depth to 68 m depth\n",
    "* A final aquifer is from 68 to 88 m depth.\n",
    "* Below 88 m depth the formations are considered an aquiclude\n",
    "\n",
    "Five pumping wells are screened in the first aquifer. The drawdown of one of them is available in the MLU documentation (Carlson & Randall, 2012).\n",
    "The provided pumping well was operated for 20 minutes at 1848 m3/d. Drawdown was recorded for a total of 50 minutes during and after pumping. The radius of the pumped well is 0.155 m.\n",
    "\n",
    "In this notebook, we reproduce the work of Xinzhu (2020) and demonstrate the use of TTim to fit a recovery test and quantify the skin resistance in the well and the hydraulic parameters of the aquifer.\n",
    "\n",
    "The following figure summarises the hydrogeological conceptual model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAE6CAYAAABeYvgTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHYElEQVR4nO3dd1gUV/s38O+isEsXpItSxN4BC2hUbGDhQRN7CZaQqLG3SIqgUbGgMepjN5hfTGyxY1QwlmgiBhVLJOqjAVEBRVRAUOp5//Bl47qgu7jrCnw/17XXxZw5c+aes8veOzNnZiRCCAEiIiLSGj1dB0BERFTRMdkSERFpGZMtERGRljHZEhERaRmTLRERkZYx2RIREWkZky0REZGWMdkSERFpGZMtERGRljHZviMuXbqEESNGwMXFBTKZDCYmJnB3d8eiRYvw8OFDXYenccnJyQgNDcWFCxfeyvo6duyIjh07qlRPIpHA1dUVJd1c7bfffoNEIoFEIsGmTZs0Ft+mTZsgkUiQmJio9rKhoaGQSCQaiwUAEhMTIZFIEB4eXuL88PDwMsdbFs7Ozhg+fLh8+vjx45BIJDh+/PhbWf+rqNP/L29HaYo/YyW9VFm+PNDG5/ZdVlXXARCwfv16jB07FvXq1cP06dPRsGFD5Ofn4+zZs1izZg1Onz6N3bt36zpMjUpOTsbs2bPh7OyM5s2b6zocBaampkhISMDRo0fRuXNnhXnfffcdzMzMkJmZqaPoCADc3d1x+vRpNGzYUNehaE3fvn0xdepUpXJra2sdRENvislWx06fPo0xY8aga9eu2LNnD6RSqXxe165dMXXqVBw6dEiHEVY+tWrVgqmpKb777juFZJuVlYUdO3ZgyJAhWL9+vQ4jrHhycnJgZGSkcn0zMzO0adNGixHpnq2tbYXfxsqEh5F1bP78+ZBIJFi3bp1Coi1mYGCA//znP/LpoqIiLFq0CPXr14dUKoWNjQ0+/PBD3LlzR2G5jh07onHjxoiNjcV7770HIyMjuLq6YsGCBSgqKlKo+/jxY0ydOhWurq7yNnv06IGrV6/K6+Tl5WHu3Lny9VpbW2PEiBFIS0tTaMvZ2Rm9evXC7t270bRpU8hkMri6umL58uXyOsePH0fLli0BACNGjJAfHgsNDZXHXtIh3+HDh8PZ2VmhbPbs2WjdujUsLS1hZmYGd3d3bNy4scRDwOoYOXIkdu3ahcePH8vLtm7dCgAYOHBgicucOnUKnTt3hqmpKYyMjODt7Y0DBw4o1YuJiUHbtm0hk8ng4OCA4OBg5Ofnl9jmtm3b4OXlBWNjY5iYmMDX1xdxcXFvtG3aEh0djYCAADg6OkImk8HNzQ2ffPIJHjx4oFCv+PDh+fPn0bdvX1hYWKB27doAgPz8fMyYMQN2dnYwMjJCu3bt8Oeffyqt6+XDyA8ePEDNmjXh7e2t0Jfx8fEwNjbGsGHDAACTJk2CsbFxiUcmBgwYAFtbW4Xly9r/qm7Hm1B1mwH135tLly6hX79+MDc3h6WlJaZMmYKCggJcu3YNfn5+MDU1hbOzMxYtWqSwfPH7snnzZkyZMgV2dnYwNDREhw4dVP7clqfPvDqYbHWosLAQR48ehYeHB2rWrKnSMmPGjMFnn32Grl27Yt++ffj6669x6NAheHt7K/3jpKamYsiQIRg6dCj27duH7t27Izg4GJs3b5bXycrKQrt27bB27VqMGDEC+/fvx5o1a1C3bl2kpKQAeJ7gAwICsGDBAgwePBgHDhzAggULEB0djY4dO+Lp06cK671w4QImTZqEyZMnY/fu3fD29sbEiRPl5//c3d0REREBAPjyyy9x+vRpnD59Gh999JHafZiYmIhPPvkE27dvx65du/D+++9j/Pjx+Prrr9Vu60UDBw5ElSpVsGXLFnnZxo0b0bdvX5iZmSnVP3HiBDp16oSMjAxs3LgRW7ZsgampKfz9/bFt2zZ5vfj4eHTu3BmPHz/Gpk2bsGbNGsTFxWHu3LlKbc6fPx+DBg1Cw4YNsX37dvzwww/IysrCe++9h/j4+DfaPlUVFRWhoKBA6fXyDzYAuHnzJry8vLB69WpERUVh1qxZOHPmDNq1a1fij4n3338fbm5u2LFjB9asWQMACAoKQnh4OD788EPs3bsXH3zwAd5//308evTolXFaWVlh69atiI2NxWeffQbg+d5yv379UKtWLXn7I0eORE5ODrZv366w/OPHj7F3714MHToU+vr6AN6s/8u6HS8SQpTY98U/JFXdZkD996Z///5o1qwZdu7ciaCgIHzzzTeYPHkyevfujZ49e2L37t3o1KkTPvvsM+zatUtp+c8//xz//PMPNmzYgA0bNiA5ORkdO3bEP//888ptfhc+81ojSGdSU1MFADFw4ECV6v/9998CgBg7dqxC+ZkzZwQA8fnnn8vLOnToIACIM2fOKNRt2LCh8PX1lU/PmTNHABDR0dGlrnfLli0CgNi5c6dCeWxsrAAgVq1aJS9zcnISEolEXLhwQaFu165dhZmZmcjOzlZYNiIiQml9HTp0EB06dFAqDwwMFE5OTqXGWVhYKPLz88WcOXNE9erVRVFR0WvbLGndjRo1kq/P09NTCCHElStXBABx/PjxEmNv06aNsLGxEVlZWfKygoIC0bhxY+Ho6CiPZcCAAcLQ0FCkpqYq1Ktfv74AIBISEoQQQiQlJYmqVauK8ePHK8SXlZUl7OzsRP/+/eVlISEhQtP/ygkJCQLAa1/F8b6sqKhI5Ofni1u3bgkAYu/evUrxzpo1S2GZ4s/35MmTFcp//PFHAUAEBgbKy44dOyYAiGPHjinUXbhwoQAgdu/eLQIDA4WhoaG4dOmSQh13d3fh7e2tULZq1SoBQFy+fFkI8Wb9r852lOZVff7DDz+ovc0vUuW9WbJkicIyzZs3FwDErl275GX5+fnC2tpavP/++/Ky4vfF3d1d4f8vMTFR6Ovri48++khpXcXU6fPyiHu25cixY8cAQGk0YqtWrdCgQQP8+uuvCuV2dnZo1aqVQlnTpk1x69Yt+fTBgwdRt25ddOnSpdT1RkZGolq1avD391f4hd28eXPY2dkpjQht1KgRmjVrplA2ePBgZGZm4vz586purkqOHj2KLl26wNzcHFWqVIG+vj5mzZqF9PR03L9//43aHjlyJM6ePYvLly9j48aNqF27Ntq3b69ULzs7G2fOnEHfvn1hYmIiL69SpQqGDRuGO3fu4Nq1awCev4edO3eGra2tQr0BAwYotHn48GEUFBTgww8/VOhzmUyGDh06qD0KV5Swl6SKiRMnIjY2Vuk1ceJEpbr379/H6NGjUbNmTVStWhX6+vpwcnICAPz9999K9T/44AOF6eLP95AhQxTK+/fvj6pVVRteMn36dPTs2RODBg3C999/jxUrVqBJkyYKdUaMGIE//vhD/p4AQEREBFq2bInGjRsDeLP+18R2FNcvqe979Oih9jar+9706tVLYbpBgwaQSCTo3r27vKxq1apwc3NT+D4pNnjwYIWRxk5OTvD29pb3TUk0/Zl/13CAlA5ZWVnByMgICQkJKtVPT08HANjb2yvNc3BwUPrQV69eXameVCpVOOyblpaGWrVqvXK99+7dw+PHj2FgYFDi/JcPX9vZ2SnVKS4r3gZN+PPPP9GtWzd07NgR69evh6OjIwwMDLBnzx7MmzdP6fC2utq3b486depg7dq12L59OyZNmlTipQqPHj2CEKLU9wX4d7vT09Nf2T/F7t27BwDyc9sv09NT73fyiRMn4OPjo1CWkJCgdA78ZY6OjvD09FQqf/mLr6ioCN26dUNycjK++uorNGnSBMbGxigqKkKbNm1KfC9e7q/iPnq5L6pWrVriZ7kkxZfGHDhwAHZ2dgrnLYsNGTIE06ZNw6ZNmxAWFob4+HjExsZi1apV8jpv0v+a2A7g+ajjkvr+Za/b5rK8N5aWlgrTBgYGMDIygkwmUyov6fx3aZ/xixcvlrodmv7Mv2uYbHWoSpUq6Ny5Mw4ePIg7d+7A0dHxlfWL/1FTUlKU6iYnJ8PKykrtGKytrZUGV73MysoK1atXL3VUtKmpqcJ0amqqUp3iMlW+bGQyGTIyMpTKX07qW7duhb6+PiIjIxW+BPbs2fPadahqxIgR+PLLLyGRSBAYGFhiHQsLC+jp6cnPcb8oOTkZAOTvTfXq1V/ZP8WK6//888/yPZA34eHhgdjYWIWy4h8CmvDXX3/h4sWL2LRpk0I/3bhxo9RlXv7hUvzZSE1NRY0aNeTlBQUFKv9IS0lJwaefformzZvjypUrmDZtmsLgPOD5+xUQEID/+7//w9y5cxEREQGZTIZBgwbJ67xJ/2tiO9Txum0uy3vzpkr7jL/q/1/Tn/l3Tfn+qVABBAcHQwiBoKAg5OXlKc3Pz8/H/v37AQCdOnUCAIUBTgAQGxuLv//+W+maUFV0794d169fx9GjR0ut06tXL6Snp6OwsBCenp5Kr3r16inUv3LlitIv2J9++gmmpqZwd3cHAPnI65J+VTs7O+P69evIzc2Vl6Wnp+OPP/5QqCeRSFC1alVUqVJFXvb06VP88MMPKm796wUGBsLf3x/Tp09X+OJ8kbGxMVq3bo1du3YpbE9RURE2b94MR0dH1K1bFwDg4+ODX3/9Vf4rHng+UO7FQVQA4Ovri6pVq+LmzZsl9rkqezwvMjU1VVq+tCMVZVGcOF8eUb927VqV2ygegf7jjz8qlG/fvl2lw96FhYUYNGgQJBIJDh48iLCwMKxYsaLEATwjRoxAcnIyfvnlF2zevBl9+vRBtWrV5PPfpP/fdDvUoco2a+K9UdeWLVsUrgi4desW/vjjj1feWEbTn/l3Dfdsdax4hODYsWPh4eGBMWPGoFGjRsjPz0dcXBzWrVuHxo0bw9/fH/Xq1cPHH3+MFStWQE9PD927d0diYiK++uor1KxZE5MnT1Z7/ZMmTcK2bdsQEBCAmTNnolWrVnj69ClOnDiBXr16wcfHBwMHDsSPP/6IHj16YOLEiWjVqhX09fVx584dHDt2DAEBAejTp4+8TQcHB/znP/9BaGgo7O3tsXnzZkRHR2PhwoXyaylr164NQ0ND/Pjjj2jQoAFMTEzg4OAABwcHDBs2DGvXrsXQoUMRFBSE9PR0LFq0SGkUcM+ePbF06VIMHjwYH3/8MdLT0xEeHl7iJVRl5eDgoNKeclhYGLp27QofHx9MmzYNBgYGWLVqFf766y9s2bJF/oX35ZdfYt++fejUqRNmzZoFIyMj/Pe//0V2drZCe87OzpgzZw6++OIL/PPPP/Dz84OFhQXu3buHP//8E8bGxpg9e7bGtvNN1a9fH7Vr18bMmTMhhIClpSX279+P6Oholdto0KABhg4dimXLlkFfXx9dunTBX3/9hfDw8BJHgL8sJCQEJ0+eRFRUFOzs7DB16lScOHECo0aNQosWLeDi4iKv261bNzg6OmLs2LFITU3FiBEjFNp6k/5/0+0odu/ePcTExCiVm5mZyW/moco2a+K9Udf9+/fRp08fBAUFISMjAyEhIZDJZAgODi51mfL2mVebLkdn0b8uXLggAgMDRa1atYSBgYEwNjYWLVq0ELNmzRL379+X1yssLBQLFy4UdevWFfr6+sLKykoMHTpU3L59W6G9F0fVvqikEb2PHj0SEydOFLVq1RL6+vrCxsZG9OzZU1y9elVeJz8/X4SHh4tmzZoJmUwmTExMRP369cUnn3wi/ve//8nrOTk5iZ49e4qff/5ZNGrUSBgYGAhnZ2exdOlSpVi2bNki6tevL/T19QUAERISIp/3/fffiwYNGgiZTCYaNmwotm3bVmLs3333nahXr56QSqXC1dVVhIWFiY0bNyqNlC3LaOTSlDaS+uTJk6JTp07C2NhYGBoaijZt2oj9+/crLf/777+LNm3aCKlUKuzs7MT06dPFunXrShzdu2fPHuHj4yPMzMyEVCoVTk5Oom/fvuLIkSPyOtocjbx48eIS5y9evFgp3vj4eNG1a1dhamoqLCwsRL9+/URSUpLSe1scb1pamlK7ubm5YurUqcLGxkbIZDLRpk0bcfr0aeHk5PTK0chRUVFCT09PYT1CCJGeni5q1aolWrZsKXJzcxXmff755wKAqFmzpigsLCxxO8va/6puR2nwitHIbdu2VXub3/S9CQwMFMbGxkpxvvz/Uvy+/PDDD2LChAnC2tpaSKVS8d5774mzZ88qLFva51aVPi+PJEK84dX/RC9wdnZG48aNERkZqetQiOgtO378OHx8fLBjxw707dtX1+G8U3jOloiISMuYbImIiLSMh5GJiIi0rNzs2YaFhaFly5YwNTWFjY0NevfurXAHGOD5XXJCQ0Ph4OAAQ0NDdOzYEVeuXNFRxERERM+Vm2R74sQJfPrpp4iJiUF0dDQKCgrQrVs3hUsmFi1ahKVLl2LlypWIjY2FnZ0dunbtiqysLB1GTkRElV25PYyclpYGGxsbnDhxAu3bt4cQAg4ODpg0aZL8CRi5ubmwtbXFwoUL8cknn+g4YiIiqqzK7U0tim/nV3wPz4SEBKSmpqJbt27yOlKpFB06dMAff/xRarLNzc1VuFNRUVERHj58iOrVq5d4H1wiIqr4hBDIysqCg4ODRu7LXC6TrRACU6ZMQbt27eRP6Si+F+eLT1Mpni7pqRTFwsLCyvddSYiISGtu37792vvWq6JcJttx48bh0qVLOHXqlNK8l/dGhRCv3EMNDg7GlClT5NMZGRmoVasWzp69DRMT1W+tRkREFceTJ5nw9Kyp9KCVsip3yXb8+PHYt28ffvvtN4VfG8WPdEpNTVV4dNf9+/eV9nZfJJVKS7yXromJGUxNmWyJiCozTZ1OLDejkYUQGDduHHbt2oWjR48q3FQcAFxcXGBnZ6dwc+28vDycOHEC3t7ebztcIiIiuXKzZ/vpp5/ip59+wt69e2Fqaio/R2tubg5DQ0NIJBJMmjQJ8+fPR506dVCnTh3Mnz8fRkZGGDx4sI6jJyKiyqzcJNvVq1cDgNLzECMiIjB8+HAAwIwZM/D06VOMHTsWjx49QuvWrREVFaWxY+5ERERlUW6vs9WWzMxMmJub4+rVDJ6zJSKqpLKyMlG/vjkyMjLUeg5xacrNOVsiIqLyismWiIhIy5hsiYiItIzJloiISMuYbImIiLSMyZaIiEjLmGyJiCqwMWMGolkzW9SrZ4YuXZoiOjpS1yFVSky2REQV2KRJXyE29jauXcvE4sUbMH78EDx8mK7rsCodJlsiogqsXr1GMDAwAABUrVoV+fl5SE29q+OoKh8mWyJ6Ky5ciMWkScPRqpUTXFykaNSoOoYN64kzZ05qdb1PnmRh7twZGDSoG5o0sUaNGhIsWRJaav2//orDyJG94e7ugNq1jdC+fX18880cPH2aU+oyP/20ATVqSFCnjonG2tSkceOGwNVVBj8/D3h7d0KDBk2U6hQWFqJpUxusW/fNK9tasiQUNWpo5kk4lQmTLRFp3aJFX8Hfvw1SUu5g2rQ52Lz5EEJCliIl5Q769u2IXbt+1Nq6Hz1Kx48/rkNeXi78/Hq/su716/EICPDG7duJCA1dhu+/j0RAwEB8880cjB07qMRlUlLu4uuvp8HOzkFjbWraypU/4vr1J/jpp8Po0KFbiY+Ni4n5DenpaejR4/23ElNlU24eREBE5VN4eAi+/XYuvvpqMUaPnqYwLyBgIDp1aoQvvvgUXbr0gpmZucbX7+johPj4R5BIJHj48AF++mlDqXV37/4Jz549w/r1O+HsXBsA0K5dJ9y7l4Iff1yHx48foVo1C4VlZs4cjdat26NaNUscOPCzRtrUhqpVq6JDh27YuHE5XFzqoHPnHgrzDxz4Gc2aecLR0UnrsVRG3LMlIq05dy4G3347F/36BSolWgCQSqUYMuRjZGZm4NSpX7USg0QiUfkB4Pr6+gCglPTNzatBT09Pfu6z2M6dmxETcwLz56/SWJsvKz5sGx9/CR9/3A/165ujUSNLhIZOQUFBAW7cuIYhQ/xQt64pWrd2xqpVi17ZXlFRIRITbyiUCSFw6NBu9OjxgUL5kSMH0LVrc7i4SNGmjQvWrAl/ZdtUOu7ZEr2Bu3eT8PDhA52s29LSCjVq1NLJulW1bNnXkEgkmD59Tql1atVyBQCkpNxRmieEQGFhoUrrqlr1zb/O+vULxIYNyzBz5hh88cVCVK9ujdOnT2Dz5rUYPvxTGBkZy+s+eHAfISGTEBy8AA4Ojhpp81VGj+6P998fiqFDP8HJk9FYtWoRCgrycfLkEQQGjsUnn0zDnj0/Yd68z+Ds7IYePd7H/fupiI39HT4+fjAwkOKXX3bhjz+OITh4gULbZ8/+gXv3UhSS7cmTv2LkyAB4eHhh1aqtKCwsxOrVi5CWdk/NXiWAyZaozO7eTUKHDg3e2iCXlxkaGuHEib/f2YSbmZmB336LQocO3V4ZY3b2EwAoMemcPn0C/fr5qLS+mJgE1KzpXKZYi9Ws6Yx9+05j1Kg+8PauLS8fNWoCZs9eplA3OHgsateuh8DAMRpr81WGDPkYn3wyBQDQvn0XnDgRhYiIldiwYRe6d+8DAPD27ogjRyKxe/eP8nOvGzYsw9SpIyGRSODiUgdr1mxHo0bNFNqOjPwZDRo0gatrHXnZwoVfwNraFlu2REMmkwEAOnb0RevWzirHTP9isiUqo4cPH+Dp0xyEhobC2dn5ra47MTERoaGhePjwwTubbK9evYyCggLUr6888vVF586dBgA0aNBUaV7Tph745ZdYldZna1vyACV13L6diMBAf1hb22Ldup9Rvbo14uLO4Ntv5yI7+wmWLNkIADhwYCeOHNmPw4fjXnuIWtU2X6dLl14K03XqNEB8/EX4+HSXl1WtWhXOzm64c+cWAMDGxg67d79+tPfBg7swcOBI+XROTjYuXoxFYOBYeaIFABMTU3Tt6o8dO75XKWb6F5Mt0RtydnZG/fr1dR3GOycrKxMAUL26dal1njzJwt69W1CrlguaNfNUmm9sbIJGjZqrtD5NHEaeP38mnjzJRHT0Bfmedps27WFpaYUpU0aib98P0bSpB7744lOMGDEetrYOyMh4DADIz88DAGRkPIa+vr58eVXa9PLq8NrYLCwsFab19Q1gaGikkAwBwMDAAE+eZKq8zXFxf+Lu3SSFQ8iPHz9CUVERrK3tlOrb2CiX0esx2RKRVtjbPz+Peft2Yql1Vq1ahCdPsvD11ytK3EN824eRr1y5gDp1Giod0m7WrCUA4Nq1v+Do6IS0tHtYu3YJ1q5dotRGw4YW8PUNwHff7VG5TVWSrbb88stOuLrWRf36jeVl1apZQCKRIC0tVan+/fvKZfR6TLZEpBUNGjSBs3Nt7NmzBdOnf610ecvOnZuxYsV8+Pv3R//+gSW28bYPI9vaOuDatb+Qnf0Exsb/3qCi+FC3vb0jrK3tsGPHMaVl//vfBYiJOYEffjgIS0srtdrUpV9+2YlevforlBkZGaN581Y4eHAXvvxysXzv+cmTLERH79dFmOUeky1RORYVdQ7xf2e91XVevXYHMTFXYWVtBmsrc1S3NEP16qawtDSFVKoPmcwAvt08IJFIsGjRenz4YQ/07NkSY8c+HyWblpaKvXu3IipqH/r3H45Fi9aVui4TE9MSDy+r6+jRg8jJyUZ29vO+un49HpGRz6+J7dy5BwwNjQAAQUGTMHJkbwwc2BVBQZNhaWmF8+djsHJlGOrWbQgfn+4wMDCAt3dHpXVs374JenpVlOap0qau/PXXBSQm3kTPnh8ozZsx42sMGeKHQYO64pNPpqKwsBCrVi2EkZExHj9+qINoyzcmW6JybOmy3YDE8vUV37K7SZsBAG3b+mD//jNYvnweFi/+CunpaSgqKoKjoxO2bIlG+/Zd3ko8wcFj5IOGACAycgciI3cAUDz83K3bf7Bt26/4738XICRkIjIzM+DgUBNDh36CceOCX3tNbEm00aam/PLLTjg6OqFpUw+lee3bd8XGjXuwaNGXGDNmAKyt7RAYOBbPnj3F0qWzdRBt+SYRQghdB/EuyczMhLm5Oa5ezYCpqZmuw6F32OXL5+Hn54FNmza99QFSV69exfDhwxH08bewd3B7q+u+ciUJsWevw8bGHDbW1VC9+vM9W1NTI0gkgKFMig+HdS51+fHjh2H//m3Yu/cPjey1Utl17Ph8zzokRPncc2WXlZWJ+vXNkZGRATOzN88F3LMlKsc+eL8dmjRx13UYapk//7+IjT2FceOG4PDh8yrf1IE07/jxeF2HUGkw2RLRW2VqaoaYmARdh0H0VvHeyERERFrGZEtERKRlTLZERERaxmRLRESkZUy2REREWsZkS0REpGVMtkRERFrGZEtERKRlFTLZrlq1Ci4uLpDJZPDw8MDJk69/eDIREZG2VLhku23bNkyaNAlffPEF4uLi8N5776F79+5ISkrSdWhERFRJVbjbNS5duhSjRo3CRx99BABYtmwZDh8+jNWrVyMsLEzldnJyslGlShVthUkVwLNnTwEAubm5ePr06Vtdd25urjyGnJzst7puospA0/9XFeqpP3l5eTAyMsKOHTvQp08fefnEiRNx4cIFnDhxQmmZ3Nxc+RcX8PypPzVr1nwr8RIR0btNU0/9qVCHkR88eIDCwkLY2toqlNva2iI1NbXEZcLCwmBubi5/MdESEZGmVbjDyAAgkUgUpoUQSmXFgoODMWXKFPl08Z7t119/DZlMptU4qfx79OgRsrP/PdzUsGFDra0rPl7xcWjGxsawsLDQ2vqIKrNnz57hq6++0lh7FSrZWllZoUqVKkp7sffv31fa2y0mlUohlUqVyg0MDEosJ3qRnZ2dwnSzZs20tq4XT3cQkXYVFRVptL0KdRjZwMAAHh4eiI6OViiPjo6Gt7e3jqIiIqLKrkLt2QLAlClTMGzYMHh6esLLywvr1q1DUlISRo8erevQiIiokqpwyXbAgAFIT0/HnDlzkJKSgsaNG+OXX36Bk5OTrkMjIqJKqsIlWwAYO3Ysxo4dq+swiIiIAFSwc7ZERETvIpX2bF+8NEZVX375JSwtLdVejoiIqKJRKdkuW7YMXl5eMDAwUKnRU6dOYdy4cUy2REREUOOc7e7du2FjY6NSXVNT0zIHREREVNGodM42IiIC5ubmKje6du3aUm8iQUREVNmotGcbGBioVqODBw8uUzBEREQV0Rtd+vPkyROlW1pp4ukIREREFYnal/4kJCSgZ8+eMDY2hrm5OSwsLGBhYYFq1arxpuhEREQlUHvPdsiQIQCA7777Dra2tqU+TYeIiIieUzvZXrp0CefOnUO9evW0EQ8REVGFo/Zh5JYtW+L27dvaiIWIiKhCUnvPdsOGDRg9ejTu3r2Lxo0bQ19fX2F+06ZNNRYcERFRRaB2sk1LS8PNmzcxYsQIeZlEIoEQAhKJBIWFhRoNkIiIqLxTO9mOHDkSLVq0wJYtWzhAiugFixcvxrNnzzBq1Cj07t1bXv7gwQOFH6evsnLlSoXHQR4+fBgrV64EAOTl5cHMzAzTp0/XaNxEpH1qJ9tbt25h3759cHNz00Y874yOdfNgYsyHIpHq5j3NRNqDDNy8eVOhvKioCGlpaSq18fKRoWfPniksK9MHOtfPffNgieiVnmTnabQ9tZNtp06dcPHixQqfbInKaseOHZg6dap8Wk9PD9bW1iotW6VKFYVpmUwGa2trpKenK91AhojKD7WTrb+/PyZPnozLly+jSZMmSgOk/vOf/2gsOKLy6OXEamVlhf3795epLV9fX/j6+sLf31/lvWMieveonWxHjx4NAJgzZ47SPA6QIiIiUqZ2suWhLCIiIvVwBBAREZGWqZRsly9fjmfPnqnc6Jo1a5CVlVXmoIiIiCoSlZLt5MmT1UqeM2bM4GAOIg2qXr06bKwtYWVhqutQiKgMVDpnK4RA586dUbWqaqd4nz59+kZBEZGiTZs2wSD3FqS5SboOhYjKQKXsGRISolajAQEBsLS0LFNAROWVUw0rGJla8rNPREq0kmyJKqOIxaORZfaersMgoncQRyMTERFpmdrX2RLR27dgwQI8eZQCC2MgdFJfXYdDRGpisiUqB37//XekpaXB1spc16EQURkw2RJpyIywn/Ag52eYm5uXeDtTIqq81D5nO2fOHOTk5CiVP336lF8wVKmdvfwPzpw5g7i4OF2HQkTvGLWT7ezZs/HkyROl8pycHMyePVsjQREREVUkaidbIQQkEolS+cWLF7V2fWFiYiJGjRoFFxcXGBoaonbt2ggJCUFenuLDfZOSkuDv7w9jY2NYWVlhwoQJSnWIiIjeNpXP2VpYWEAikUAikaBu3boKCbewsBBPnjyRP35P065evYqioiKsXbsWbm5u+OuvvxAUFITs7GyEh4fLY+jZsyesra1x6tQppKenIzAwEEIIrFixQitxERERqULlZLts2TIIITBy5EjMnj0b5ub/joo0MDCAs7MzvLy8tBKkn58f/Pz85NOurq64du0aVq9eLU+2UVFRiI+Px+3bt+Hg4AAAWLJkCYYPH4558+bBzMxMK7ERERG9jsrJNjAwEADg4uICb29v6Ovray0oVWRkZCgctj59+jQaN24sT7QA4Ovri9zcXJw7dw4+Pj4ltpObm4vc3Fz5dGZmpvaCJiKiSkntS386dOiAoqIiXL9+Hffv31d6mHz79u01Flxpbt68iRUrVmDJkiXystTUVNja2irUs7CwgIGBAVJTU0ttKywsjAO7iIhIq9ROtjExMRg8eDBu3boFIYTCPIlEgsLCQpXbCg0NfW2ii42Nhaenp3w6OTkZfn5+6NevHz766COl9b+stAFdxYKDgzFlyhT5dGZmJmrWrKnqJhC9Fd26dUP242RYGBa9vjIRvXPUTrajR4+Gp6cnDhw4AHt7+1cmstcZN24cBg4c+Mo6zs7O8r+Tk5Ph4+MDLy8vrFu3TqGenZ0dzpw5o1D26NEj5OfnK+3xvkgqlUIqlaofPNFbNH78eD5ij6gcUzvZ/u9//8PPP/8MNze3N165lZUVrKysVKp79+5d+Pj4wMPDAxEREdDTU7xqycvLC/PmzUNKSgrs7e0BPB80JZVK4eHh8caxEr1O3+6tkJ5fHSYmJroOhYjeMWon29atW+PGjRsaSbaqSk5ORseOHVGrVi2Eh4cjLS1NPs/Ozg7A88NsDRs2xLBhw7B48WI8fPgQ06ZNQ1BQEEci01sxdlg3PmKPiEqkUrK9dOmS/O/x48dj6tSpSE1NRZMmTZRGJTdt2lSzEeL5HuqNGzdw48YNODo6KswrPm9cpUoVHDhwAGPHjkXbtm1haGiIwYMHyy8NIiIi0hWJeHmUUwn09PQgkUiUBkTJG/n/89QdIPUuyszMhLm5Oc7s/homxjJdh0PljLb2bAcMGIAHafdhY2mCyO9maGUdRPSvJ9nP0LrPV8jIyNDI0VGV9mwTEhLeeEVEVHY5OTnIznmKHCMDXYdCRGWgUrJ1cnLSdhxE5V6nwXNx70EGrK2tsX//fl2HQ0TvELUHSO3bt6/EcolEAplMBjc3N7i4uLxxYERERBWF2sm2d+/eJZ6/ffG8bbt27bBnzx5YWFhoLFAiIqLySu1H7EVHR6Nly5aIjo5GRkYGMjIyEB0djVatWiEyMhK//fYb0tPTMW3aNG3ES0REVO6ovWc7ceJErFu3Dt7e3vKyzp07QyaT4eOPP8aVK1ewbNkyjBw5UqOBEhERlVdq79nevHmzxGHQZmZm+OeffwAAderUwYMHD948OiIiogpA7WTr4eGB6dOnK9zFKS0tDTNmzEDLli0BPL+l48s3nyAiIqqs1D6MvHHjRgQEBMDR0RE1a9aERCJBUlISXF1dsXfvXgDAkydP8NVXX2k8WCIiovJI7WRbr149/P333zh8+DCuX78OIQTq16+Prl27yh8O0Lt3b03HSVSpffbZZyh8cgcmehm6DoWIykDtZAs8v8zHz88Pfn5+mo6HiErQrl07PmKPqBxTKdkuX74cH3/8MWQyGZYvX/7KuhMmTNBIYETlzYLPBiFDv77SwzmIiFR6EIGLiwvOnj2L6tWrv/LuUBKJRD4iubzigwjoTWjzEXvcsyV6e3T+IAI+lIDo7bt69SpE9m0YiwdoVJcj/YnKmzKdswWAvLw8JCQkoHbt2qhatczNEJEKii+3s7Uyx9GfvtR1OESkJrWvs83JycGoUaNgZGSERo0aISnp+WGtCRMmYMGCBRoPkKi8+PPiTcTExODcuXO6DoWI3jFqJ9vg4GBcvHgRx48fh0z27znNLl26YNu2bRoNjqg8mblwCyZNmoTQ0FBdh0JE7xi1j//u2bMH27ZtQ5s2bSCRSOTlDRs2xM2bNzUaHBERUUWg9p5tWloabGxslMqzs7MVki8RERE9p3aybdmyJQ4cOCCfLk6w69evh5eXl+YiIyIiqiDUPowcFhYGPz8/xMfHo6CgAN9++y2uXLmC06dP48SJE9qIkYiIqFxTe8/W29sbv//+O3JyclC7dm1ERUXB1tYWp0+fhoeHhzZiJCIiKtfKdIFskyZN8P3332s6FiIiogqpTMm2qKgIN27cwP3791FUVKQwr3379hoJjIiIqKJQO9nGxMRg8ODBuHXrFl6+rbJEIkFhYaHGgiOi57Zu3Qr93CTI8u7oOhQiKgO1k+3o0aPh6emJAwcOwN7enpf7EL0FxsbGMKhqBGlVPhyDqDxSO9n+73//w88//ww3NzdtxENUbh396UutPvWHiMovtUcjt27dGjdu3NBGLERERBWSSnu2ly5dkv89fvx4TJ06FampqWjSpInSg7KbNm2q2QiJCD/99BOeZSbD3OAZhvftoOtwiEhNKiXb5s2bQyKRKAyIGjlypPzv4nkcIEWV3U8//YQtW7a8tl69evUQHh6uUDZt2jRcu3atxPrp6ekoKiqCrZU5ky1ROaT2w+Mri+PXDWBoKNV1GFSOPHv2DFevXkVaWtpr61pYWODu3bsKZffv33/tss/ygV+v8nNJpG1Pnxa9vpIaVEq2Tk5OGl0pUUUkk8lgZGQEc3Pz19Y1MDBQSrYGBgavXdbMzOyNYiQi3SjTTS10KTc3F61bt8bFixcRFxeH5s2by+clJSXh008/xdGjR2FoaIjBgwcjPDwcBgYGuguYKpVOnTqhU6dOZVr2448/1nA0RPSuKHfJdsaMGXBwcMDFixcVygsLC9GzZ09YW1vj1KlTSE9PR2BgIIQQWLFihY6iJSIiKsOlP7p08OBBREVFKQ0sAYCoqCjEx8dj8+bNaNGiBbp06YIlS5Zg/fr1yMzM1EG0REREz5WbZHvv3j0EBQXhhx9+gJGRkdL806dPo3HjxnBwcJCX+fr6Ijc3F+fOnSu13dzcXGRmZiq8iIiINKlMyfbx48fYsGEDgoOD8fDhQwDA+fPnlQZ8aIoQAsOHD5ffKrIkqampsLW1VSizsLCAgYEBUlNTS207LCwM5ubm8lfNmjU1GjsREZHayfbSpUuoW7cuFi5ciPDwcDx+/BgAsHv3bgQHB6vVVmhoKCQSyStfZ8+exYoVK5CZmfna9ku6T3Px9b+lCQ4ORkZGhvx1+/ZttbaBiIjoddQeIDVlyhQMHz4cixYtgqmpqby8e/fuGDx4sFptjRs3DgMHDnxlHWdnZ8ydOxcxMTGQShWvL/T09MSQIUPw/fffw87ODmfOnFGY/+jRI+Tn5yvt8b5IKpUqtUtERKRJaifb2NhYrF27Vqm8Ro0arzxcWxIrKytYWVm9tt7y5csxd+5c+XRycjJ8fX2xbds2tG7dGgDg5eWFefPmISUlBfb29gCeD5qSSqXw8PBQKy4iIiJNUjvZymSyEgcRXbt2DdbW1hoJ6mW1atVSmDYxMQEA1K5dG46OjgCAbt26oWHDhhg2bBgWL16Mhw8fYtq0aQgKCuKNAIiISKfUPmcbEBCAOXPmID8/H8Dz86RJSUmYOXMmPvjgA40HqKoqVargwIEDkMlkaNu2Lfr374/evXuXeJkQERHR26T2nm14eDh69OgBGxsbPH36FB06dEBqaqr8MO7b4OzsrPBQhGK1atVCZGTkW4mBiIhIVWonWzMzM5w6dQpHjx7F+fPnUVRUBHd3d3Tp0kUb8REREZV7aifbxMREODs7v9E9YImIiCoTtc/Zurq6ol27dli7dq38hhZERERUOrWT7dmzZ+Hl5YW5c+fCwcEBAQEB2LFjB3Jzc7URHxERUbmndrJ1d3fH4sWLkZSUhIMHD8LGxgaffPIJbGxsMHLkSG3ESEREVK6V+UEEEokEPj4+WL9+PY4cOQJXV1d8//33moyNiIioQihzsr19+zYWLVqE5s2bo2XLljA2NsbKlSs1GRsREVGFoPZo5HXr1uHHH3/E77//jnr16mHIkCHYs2cPnJ2dtRAeERFR+ad2sv36668xcOBAfPvtt2jevLkWQiIiIqpY1E62SUlJr3xkHRERESlSKdleunQJjRs3hp6eHi5fvvzKuk2bNtVIYERERBWFSsm2efPmSE1NhY2NDZo3bw6JRKJwb+LiaYlEgsLCQq0FS0REVB6plGwTEhLkj89LSEjQakBEREQVjUrJ1snJSf73rVu34O3tjapVFRctKCjAH3/8oVCXiIiIynCdrY+PT4n3RM7IyICPj49GgiIiIqpI1E62xedmX5aeng5jY2ONBEVERFSRqHzpz/vvvw/g+WCo4cOHQyqVyucVFhbi0qVL8Pb21nyERERE5ZzKydbc3BzA8z1bU1NTGBoayucZGBigTZs2CAoK0nyERERE5ZzKyTYiIgIA4OzsjGnTpvGQMRERkYrUvoNUSEiINuIgIiKqsNROtgDw888/Y/v27UhKSkJeXp7CvPPnz2skMCIioopC7dHIy5cvx4gRI2BjY4O4uDi0atUK1atXxz///IPu3btrI0YiIqJyTe1ku2rVKqxbtw4rV66EgYEBZsyYgejoaEyYMAEZGRnaiJGIiKhcUzvZJiUlyS/xMTQ0RFZWFgBg2LBh2LJli2ajIyIiqgDUTrZ2dnZIT08H8Pw2jjExMQCe3zP5xYcTEBER0XNqJ9tOnTph//79AIBRo0Zh8uTJ6Nq1KwYMGIA+ffpoPEAiIqLyTu3RyOvWrUNRUREAYPTo0bC0tMSpU6fg7++P0aNHazxAIiKi8k7tZKunpwc9vX93iPv374/+/ftrNCgiIqKKRKVke+nSJZUbbNq0aZmDISIiqohUSrbNmzeHRCJ57QAoiUSCwsJCjQRGRERUUaiUbBMSErQdBxERUYWlUrJ1cnLSdhxEREQVltqX/gDADz/8gLZt28LBwQG3bt0CACxbtgx79+7VaHAvO3DgAFq3bg1DQ0NYWVnJn7FbLCkpCf7+/jA2NoaVlRUmTJigdO9mIiKit03tZLt69WpMmTIFPXr0wOPHj+XnaKtVq4Zly5ZpOj65nTt3YtiwYRgxYgQuXryI33//HYMHD5bPLywsRM+ePZGdnY1Tp05h69at2LlzJ6ZOnaq1mIiIiFSh9qU/K1aswPr169G7d28sWLBAXu7p6Ylp06ZpNLhiBQUFmDhxIhYvXoxRo0bJy+vVqyf/OyoqCvHx8bh9+zYcHBwAAEuWLMHw4cMxb948mJmZaSU2IiKi11F7zzYhIQEtWrRQKpdKpcjOztZIUC87f/487t69Cz09PbRo0QL29vbo3r07rly5Iq9z+vRpNG7cWJ5oAcDX1xe5ubk4d+5cqW3n5uYiMzNT4UVERKRJaidbFxcXXLhwQan84MGDaNiwoSZiUvLPP/8AAEJDQ/Hll18iMjISFhYW6NChAx4+fAgASE1Nha2trcJyFhYWMDAwQGpqaqlth4WFwdzcXP6qWbOmVraBiIgqL7WT7fTp0/Hpp59i27ZtEELgzz//xLx58/D5559j+vTparUVGhoKiUTyytfZs2flt4f84osv8MEHH8DDwwMRERGQSCTYsWOHvD2JRKK0DiFEieXFgoODkZGRIX/dvn1brW0gIiJ6HbXP2Y4YMQIFBQWYMWMGcnJyMHjwYNSoUQPffvstBg4cqFZb48aNe+0yzs7O8sf4vbjnLJVK4erqiqSkJADPn0Z05swZhWUfPXqE/Px8pT3eF0mlUkilUrXiJiIiUofayRYAgoKCEBQUhAcPHqCoqAg2NjYAgLt376JGjRoqt2NlZQUrK6vX1vPw8IBUKsW1a9fQrl07AEB+fj4SExPl1wB7eXlh3rx5SElJgb29PYDng6akUik8PDzU3UQiIiKNKdN1tsWsrKxgY2OD1NRUjB8/Hm5ubpqKS4GZmRlGjx6NkJAQREVF4dq1axgzZgwAoF+/fgCAbt26oWHDhhg2bBji4uLw66+/Ytq0aQgKCuJIZCIi0imVk+3jx48xZMgQWFtbw8HBAcuXL0dRURFmzZoFV1dXxMTE4LvvvtNaoIsXL8bAgQMxbNgwtGzZErdu3cLRo0dhYWEBAKhSpQoOHDgAmUyGtm3bon///ujduzfCw8O1FhMREZEqJOJ1Txf4/8aOHYv9+/djwIABOHToEP7++2/4+vri2bNnCAkJQYcOHbQd61uRmZkJc3NzLFy4EIaGhroOh4iIdODp06f47LPPkJGRoZGjoyqfsz1w4AAiIiLQpUsXjB07Fm5ubqhbt65W7xpFRERUEah8GDk5OVk+GtjV1RUymQwfffSR1gIjIiKqKFROtkVFRdDX15dPV6lSBcbGxloJioiIqCJR+TCyEALDhw+XX5P67NkzjB49Winh7tq1S7MREhERlXMqJ9vAwECF6aFDh2o8GCIioopI5WQbERGhzTiIiIgqrDe6qQURERG9HpMtERGRljHZEhERaRmTLRERkZaV6ak/lUHHunkwMeZvESKiyuhJdp5G22M2ISIi0jImWyIiIi1jsiUiItIyJlsiIiItY7IlIiLSMiZbIiIiLWOyJSIi0jImWyIiIi1jsiUiItIyJlsiIiItY7IlIiLSMiZbIiIiLWOyJSIi0jImWyIiIi1jsiUiItIyJlsiIiItY7IlIiLSMiZbIiIiLWOyJSIi0jImWyIiIi1jsiUiItKycpNsr1+/joCAAFhZWcHMzAxt27bFsWPHFOokJSXB398fxsbGsLKywoQJE5CXl6ejiImIiJ4rN8m2Z8+eKCgowNGjR3Hu3Dk0b94cvXr1QmpqKgCgsLAQPXv2RHZ2Nk6dOoWtW7di586dmDp1qo4jJyKiyq5cJNsHDx7gxo0bmDlzJpo2bYo6depgwYIFyMnJwZUrVwAAUVFRiI+Px+bNm9GiRQt06dIFS5Yswfr165GZmanjLSAiosqsXCTb6tWro0GDBvi///s/ZGdno6CgAGvXroWtrS08PDwAAKdPn0bjxo3h4OAgX87X1xe5ubk4d+5cqW3n5uYiMzNT4UVERKRJVXUdgCokEgmio6MREBAAU1NT6OnpwdbWFocOHUK1atUAAKmpqbC1tVVYzsLCAgYGBvJDzSUJCwvD7NmztRk+ERFVcjrdsw0NDYVEInnl6+zZsxBCYOzYsbCxscHJkyfx559/IiAgAL169UJKSoq8PYlEorQOIUSJ5cWCg4ORkZEhf92+fVsr20pERJWXTvdsx40bh4EDB76yjrOzM44ePYrIyEg8evQIZmZmAIBVq1YhOjoa33//PWbOnAk7OzucOXNGYdlHjx4hPz9faY/3RVKpFFKp9M03hoiIqBQ6TbZWVlawsrJ6bb2cnBwAgJ6e4o64np4eioqKAABeXl6YN28eUlJSYG9vD+D5oCmpVCo/r0tERKQL5WKAlJeXFywsLBAYGIiLFy/i+vXrmD59OhISEtCzZ08AQLdu3dCwYUMMGzYMcXFx+PXXXzFt2jQEBQXJ94aJiIh0oVwkWysrKxw6dAhPnjxBp06d4OnpiVOnTmHv3r1o1qwZAKBKlSo4cOAAZDIZ2rZti/79+6N3794IDw/XcfRERFTZlYvRyADg6emJw4cPv7JOrVq1EBkZ+ZYiIiIiUk252LMlIiIqz5hsiYiItIzJloiISMuYbImIiLSMyZaIiEjLmGyJiIi0jMmWiIhIy5hsiYiItIzJloiISMvKzR2k3hYhBADgSc4zHUdCRES6UpwDinPCm5IITbVUQfzzzz+oXbu2rsMgIqJ3wM2bN+Hq6vrG7XDP9iWWlpYAgKSkJJibm+s4mvIhMzMTNWvWxO3bt/mEJTWw39THPisb9pv6MjIyUKtWLXlOeFNMti8pfmauubk5P5RqMjMzY5+VAftNfeyzsmG/qe/l56iXuR2NtEJERESlYrIlIiLSMibbl0ilUoSEhEAqleo6lHKDfVY27Df1sc/Khv2mPk33GUcjExERaRn3bImIiLSMyZaIiEjLmGyJiIi0jMmWiIhIy5hsASQmJmLUqFFwcXGBoaEhateujZCQEOTl5SnUS0pKgr+/P4yNjWFlZYUJEyYo1als5s2bB29vbxgZGaFatWol1mG/KVu1ahVcXFwgk8ng4eGBkydP6jqkd8pvv/0Gf39/ODg4QCKRYM+ePQrzhRAIDQ2Fg4MDDA0N0bFjR1y5ckU3wb4jwsLC0LJlS5iamsLGxga9e/fGtWvXFOqw3xStXr0aTZs2ld/sw8vLCwcPHpTP12R/MdkCuHr1KoqKirB27VpcuXIF33zzDdasWYPPP/9cXqewsBA9e/ZEdnY2Tp06ha1bt2Lnzp2YOnWqDiPXvby8PPTr1w9jxowpcT77Tdm2bdswadIkfPHFF4iLi8N7772H7t27IykpSdehvTOys7PRrFkzrFy5ssT5ixYtwtKlS7Fy5UrExsbCzs4OXbt2RVZW1luO9N1x4sQJfPrpp4iJiUF0dDQKCgrQrVs3ZGdny+uw3xQ5OjpiwYIFOHv2LM6ePYtOnTohICBAnlA12l+CSrRo0SLh4uIin/7ll1+Enp6euHv3rrxsy5YtQiqVioyMDF2E+E6JiIgQ5ubmSuXsN2WtWrUSo0ePViirX7++mDlzpo4iercBELt375ZPFxUVCTs7O7FgwQJ52bNnz4S5ublYs2aNDiJ8N92/f18AECdOnBBCsN9UZWFhITZs2KDx/uKebSkyMjIUbkB9+vRpNG7cGA4ODvIyX19f5Obm4ty5c7oIsVxgvynKy8vDuXPn0K1bN4Xybt264Y8//tBRVOVLQkICUlNTFfpQKpWiQ4cO7MMXZGRkAPj34Srst1crLCzE1q1bkZ2dDS8vL433F5NtCW7evIkVK1Zg9OjR8rLU1FTY2toq1LOwsICBgQFSU1PfdojlBvtN0YMHD1BYWKjUJ7a2tpWyP8qiuJ/Yh6UTQmDKlClo164dGjduDID9VprLly/DxMQEUqkUo0ePxu7du9GwYUON91eFTrahoaGQSCSvfJ09e1ZhmeTkZPj5+aFfv3746KOPFOZJJBKldQghSiwvz8rSb69SWfpNHS9ve2Xvj7JgH5Zu3LhxuHTpErZs2aI0j/2mqF69erhw4QJiYmIwZswYBAYGIj4+Xj5fU/1VoR+xN27cOAwcOPCVdZydneV/Jycnw8fHB15eXli3bp1CPTs7O5w5c0ah7NGjR8jPz1f65VPeqdtvr1KZ+k0VVlZWqFKlitIv4/v371fK/igLOzs7AM/31Ozt7eXl7MPnxo8fj3379uG3336Do6OjvJz9VjIDAwO4ubkBADw9PREbG4tvv/0Wn332GQDN9VeF3rO1srJC/fr1X/mSyWQAgLt376Jjx45wd3dHRESE0jMMvby88NdffyElJUVeFhUVBalUCg8Pj7e6XdqmTr+9TmXqN1UYGBjAw8MD0dHRCuXR0dHw9vbWUVTli4uLC+zs7BT6MC8vDydOnKjUfSiEwLhx47Br1y4cPXoULi4uCvPZb6oRQiA3N1fz/aWBwVvl3t27d4Wbm5vo1KmTuHPnjkhJSZG/ihUUFIjGjRuLzp07i/Pnz4sjR44IR0dHMW7cOB1Grnu3bt0ScXFxYvbs2cLExETExcWJuLg4kZWVJYRgv5Vk69atQl9fX2zcuFHEx8eLSZMmCWNjY5GYmKjr0N4ZWVlZ8s8SALF06VIRFxcnbt26JYQQYsGCBcLc3Fzs2rVLXL58WQwaNEjY29uLzMxMHUeuO2PGjBHm5ubi+PHjCt9hOTk58jrsN0XBwcHit99+EwkJCeLSpUvi888/F3p6eiIqKkoIodn+YrIVzy9bAVDi60W3bt0SPXv2FIaGhsLS0lKMGzdOPHv2TEdRvxsCAwNL7Ldjx47J67DflP33v/8VTk5OwsDAQLi7u8svz6Dnjh07VuLnKjAwUAjx/DKWkJAQYWdnJ6RSqWjfvr24fPmyboPWsdK+wyIiIuR12G+KRo4cKf8/tLa2Fp07d5YnWiE02198xB4REZGWVehztkRERO8CJlsiIiItY7IlIiLSMiZbIiIiLWOyJSIi0jImWyIiIi1jsiUiItIyJlsiIiItY7IlegdJJBLs2bNH12Fo3PDhw+VPjnrT7Xvx6VTLli3TSHxE2sJkS/SWvJho9PX1YWtri65du+K7775DUVGRQt2UlBR0795dpXbLW2L28/NTa/tKM23aNKSkpCg82YboXcVkS/QWFSeaxMREHDx4ED4+Ppg4cSJ69eqFgoICeT07OztIpVIdRqo9UqlUI9tnYmICOzs7VKlSRUOREWkPky3RW1ScaGrUqAF3d3d8/vnn2Lt3Lw4ePIhNmzbJ6724t5qXl4dx48bB3t4eMpkMzs7OCAsLA/Dvc4X79OkDiUQin7558yYCAgJga2sLExMTtGzZEkeOHFGIxdnZGfPnz8fIkSNhamqKWrVqKT3H+c6dOxg4cCAsLS1hbGwMT09PhecT79+/Hx4eHpDJZHB1dcXs2bMVfjSoIjExERKJBNu3b8d7770HQ0NDtGzZEtevX0dsbCw8PT1hYmICPz8/pKWlqdU20buCyZZIxzp16oRmzZph165dJc5fvnw59u3bh+3bt+PatWvYvHmzPKnGxsYCACIiIpCSkiKffvLkCXr06IEjR44gLi4Ovr6+8Pf3R1JSkkLbS5YsgaenJ+Li4jB27FiMGTMGV69elbfRoUMHJCcnY9++fbh48SJmzJghP+R9+PBhDB06FBMmTEB8fDzWrl2LTZs2Yd68eWXqh5CQEHz55Zc4f/48qlatikGDBmHGjBn49ttvcfLkSdy8eROzZs0qU9tEOqeZBxUR0esEBgaKgICAEucNGDBANGjQQD4NQOzevVsIIcT48eNFp06dRFFRUYnLvlj3VRo2bChWrFghn3ZychJDhw6VTxcVFQkbGxuxevVqIYQQa9euFaampiI9Pb3E9t577z0xf/58hbIffvhB2NvblxpDSX2QkJAgAIgNGzbIy7Zs2SIAiF9//VVeFhYWJurVq6fUppOTk/jmm29KXSfRu6CqblM9EQGAEAISiaTEecOHD0fXrl1Rr149+Pn5oVevXujWrdsr28vOzsbs2bMRGRmJ5ORkFBQU4OnTp0p7tk2bNpX/LZFIYGdnh/v37wMALly4gBYtWsDS0rLEdZw7dw6xsbEKe7KFhYV49uwZcnJyYGRkpNK2lxSLra0tAKBJkyYKZcWxEZU3TLZE74C///4bLi4uJc5zd3dHQkICDh48iCNHjqB///7o0qULfv7551Lbmz59Og4fPozw8HC4ubnB0NAQffv2RV5enkI9fX19hWmJRCI/TGxoaPjKmIuKijB79my8//77SvNkMtkrly3Ji7EU//B4uezlUdtE5QWTLZGOHT16FJcvX8bkyZNLrWNmZoYBAwZgwIAB6Nu3L/z8/PDw4UNYWlpCX18fhYWFCvVPnjyJ4cOHo0+fPgCen39NTExUK66mTZtiw4YN8vW8zN3dHdeuXYObm5ta7RJVRky2RG9Rbm4uUlNTUVhYiHv37uHQoUMICwtDr1698OGHH5a4zDfffAN7e3s0b94cenp62LFjB+zs7FCtWjUAz0cV//rrr2jbti2kUiksLCzg5uaGXbt2wd/fHxKJBF999ZXae4WDBg3C/Pnz0bt3b4SFhcHe3h5xcXFwcHCAl5cXZs2ahV69eqFmzZro168f9PT0cOnSJVy+fBlz5859064iqlA4GpnoLTp06BDs7e3h7OwMPz8/HDt2DMuXL8fevXtLvV7UxMQECxcuhKenJ1q2bInExET88ssv0NN7/u+7ZMkSREdHo2bNmmjRogWA5wnawsIC3t7e8Pf3h6+vL9zd3dWK1cDAAFFRUbCxsUGPHj3QpEkTLFiwQB6nr68vIiMjER0djZYtW6JNmzZYunQpnJyc3qCHiComiRBC6DoIIqochg8fjsePH2v0jlfOzs6YNGkSJk2apLE2iTSNe7ZE9FZFRkbCxMQEkZGRb9TO/PnzYWJiojTCmuhdxD1bInpr7t+/j8zMTACAvb09jI2Ny9zWw4cP8fDhQwCAtbU1zM3NNRIjkTYw2RIREWkZDyMTERFpGZMtERGRljHZEhERaRmTLRERkZYx2RIREWkZky0REZGWMdkSERFpGZMtERGRlv0/tZI6opee0zkAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "##Now printing the conceptual model figure:\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "# sky\n",
    "sky = plt.Rectangle((-20, 0), width=50, height=20, fc=\"b\", zorder=0, alpha=0.1)\n",
    "ax.add_patch(sky)\n",
    "\n",
    "# Aquifer:\n",
    "ground = plt.Rectangle(\n",
    "    (-20, -37),\n",
    "    width=50,\n",
    "    height=27,\n",
    "    fc=np.array([209, 179, 127]) / 255,\n",
    "    zorder=0,\n",
    "    alpha=0.9,\n",
    ")\n",
    "ax.add_patch(ground)\n",
    "\n",
    "# Aquifer 2:\n",
    "ground2 = plt.Rectangle(\n",
    "    (-20, -88),\n",
    "    width=50,\n",
    "    height=20,\n",
    "    fc=np.array([209, 179, 127]) / 255,\n",
    "    zorder=0,\n",
    "    alpha=0.9,\n",
    ")\n",
    "ax.add_patch(ground2)\n",
    "\n",
    "# Confining bed:\n",
    "confining_unit = plt.Rectangle(\n",
    "    (-20, -10),\n",
    "    width=50,\n",
    "    height=10,\n",
    "    fc=np.array([100, 100, 100]) / 255,\n",
    "    zorder=0,\n",
    "    alpha=0.7,\n",
    ")\n",
    "ax.add_patch(confining_unit)\n",
    "\n",
    "confining_unit2 = plt.Rectangle(\n",
    "    (-20, -68),\n",
    "    width=50,\n",
    "    height=31,\n",
    "    fc=np.array([100, 100, 100]) / 255,\n",
    "    zorder=0,\n",
    "    alpha=0.7,\n",
    ")\n",
    "ax.add_patch(confining_unit2)\n",
    "\n",
    "well = plt.Rectangle(\n",
    "    (-1.5, -37), width=3, height=37, fc=np.array([200, 200, 200]) / 255, zorder=1\n",
    ")\n",
    "ax.add_patch(well)\n",
    "\n",
    "# Wellhead\n",
    "wellhead = plt.Rectangle(\n",
    "    (-2, 0), width=4, height=10, fc=np.array([200, 200, 200]) / 255, zorder=2, ec=\"k\"\n",
    ")\n",
    "ax.add_patch(wellhead)\n",
    "\n",
    "# Screen for the well:\n",
    "screen = plt.Rectangle(\n",
    "    (-1.5, -37),\n",
    "    width=3,\n",
    "    height=27,\n",
    "    fc=np.array([200, 200, 200]) / 255,\n",
    "    alpha=1,\n",
    "    zorder=2,\n",
    "    ec=\"k\",\n",
    "    ls=\"--\",\n",
    ")\n",
    "screen.set_linewidth(2)\n",
    "ax.add_patch(screen)\n",
    "pumping_arrow = plt.Arrow(x=2, y=7, dx=5, dy=0, color=\"#00035b\")\n",
    "ax.add_patch(pumping_arrow)\n",
    "ax.text(x=7, y=7, s=r\"$ Q = 1848$ m$^3$/d\", fontsize=\"large\")\n",
    "\n",
    "# last line\n",
    "line = plt.Line2D(xdata=[-200, 1200], ydata=[0, 0], color=\"k\")\n",
    "ax.add_line(line)\n",
    "\n",
    "\n",
    "ax.set_xlim([-20, 30])\n",
    "ax.set_ylim([-88, 20])\n",
    "ax.set_xlabel(\"Distance [m]\")\n",
    "ax.set_ylabel(\"Relative height [m]\")\n",
    "ax.set_title(\"Conceptual Model - Hardixveld Example\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = 27  # aquifer thickness, m\n",
    "zt = -10  # upper boundary of aquifer, m\n",
    "zb = zt - H  # lower boundary of the aquifer, m\n",
    "rw = 0.155  # well screen radius, m\n",
    "Q = 1848  # constant discharge rate, m^3/d\n",
    "t0 = 0.013889  # time stop pumping, d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data for the recovery test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is saved in a text file where the first column is the time data in days and in the second is the drawdown in m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt(\"data/recovery.txt\", skiprows=1)\n",
    "to = data[:, 0]\n",
    "ho = data[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conceptual model\n",
    "\n",
    "Here we create a two aquifer leaky model in ModelMaq, so we have to define the top of the first aquitard layer, followed by the tops and bottoms of the aquifer layers. Here we ignore storage (```Sll```) of the aquitard layers.\n",
    "\n",
    "The well is defined with skin resistance (```res```) and the pumping schedule with the start and shutdown of the pump."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self.neq  1\n",
      "solution complete\n"
     ]
    }
   ],
   "source": [
    "ml1 = ttm.ModelMaq(\n",
    "    kaq=[50, 40],\n",
    "    z=[0, zt, zb, -68, -88],\n",
    "    c=[1000, 1000],\n",
    "    Saq=[1e-4, 5e-5],\n",
    "    topboundary=\"semi\",\n",
    "    tmin=1e-4,\n",
    "    tmax=0.04,\n",
    ")\n",
    "w1 = ttm.Well(ml1, xw=0, yw=0, rw=rw, res=1, tsandQ=[(0, Q), (t0, 0)], layers=0)\n",
    "ml1.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters to be calibrated are the hydraulic conductivity and specific storage of the first layer, and the skin resistance of the well. The parameters of the aquitards and the second aquifer are kept fixed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".....................................................................................................................................................................................................................................................................................................................................................................................\n",
      "Fit succeeded.\n"
     ]
    }
   ],
   "source": [
    "ca1 = ttm.Calibrate(ml1)\n",
    "ca1.set_parameter(name=\"kaq0\", initial=50, pmin=0)\n",
    "ca1.set_parameter(name=\"Saq0\", initial=1e-4, pmin=0)\n",
    "ca1.set_parameter_by_reference(name=\"res\", parameter=w1.res[:], initial=1, pmin=0)\n",
    "ca1.seriesinwell(name=\"obs\", element=w1, t=to, h=ho)\n",
    "ca1.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "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>kaq0</th>\n",
       "      <td>44.529551</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq0</th>\n",
       "      <td>0.000006</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>res</th>\n",
       "      <td>0.016205</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        optimal\n",
       "kaq0  44.529551\n",
       "Saq0   0.000006\n",
       "res    0.016205"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 0.005511886856381756\n"
     ]
    }
   ],
   "source": [
    "display(ca1.parameters.loc[:, [\"optimal\"]])\n",
    "print(\"RMSE:\", ca1.rmse())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdUAAAE/CAYAAAAQZlkTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABND0lEQVR4nO3deVxU9f7H8ddhmGEVEEEBBcR9CXBDQ1xbLMlyuRppueu1otCflVmWe2qZpt00SyvL7i0rl9LMXdw1cklNcwtCE3FLUREYZs7vD2RyBJHBGWZGP8/HYx4yZ87yGWB88/2e7/keRVVVFSGEEELcMRd7FyCEEELcLSRUhRBCCCuRUBVCCCGsREJVCCGEsBIJVSGEEMJKJFSFEEIIK5FQFUIIIaxEQlUIIYSwEglVIYQQwkokVIVDmz9/PoqioCgKycnJRV5XVZVatWqhKArt2rWz6rEVRWHs2LEWb5eWloaiKMyfP79U6xU+XFxcqFixIg8++CCrV68uW9FW1q5dO7Pva3Z2NmPHji32Z1FeCr9f/fr1K/b18ePHm9ZJS0uz2nH79etH9erVy7Ttzd/HW9myZQuDBg2iadOmuLm5Wf09CNuTUBVOoUKFCnzyySdFlm/cuJHjx49ToUIFO1RlHS+++CLbt29n8+bNvPvuuxw9epT4+Hg2bdpk79KKyM7OZty4cXYNVSj4ffj222+5fPmy2XJVVZk/fz4+Pj52quzOrFu3jrVr1xIWFkbLli3tXY4oAwlV4RQSEhJYtGgRWVlZZss/+eQTYmNjCQsLs1Nldy4sLIz777+fuLg4Bg4cyJdffonBYCj2jwhRoHPnzqiqytdff222fP369aSmppKQkGCnyu7Mm2++SVpaGkuWLOGxxx6zdzmiDCRUhVPo2bMnAF999ZVp2aVLl1i0aBEDBgwodpsLFy7w/PPPU7VqVXQ6HTVq1GDUqFHk5uaarZeVlcXgwYOpVKkS3t7ePProoxw5cqTYfR49epRevXpRuXJl3NzcqF+/PrNmzbLSuyzQrFkzADIzM82Wnz59miFDhlCtWjV0Oh0RERGMGzeO/Px8s/U+/PBDoqOj8fb2pkKFCtSrV4/XX3/d9PrYsWNRFKXIcQu72m/V3ZiWlkZgYCAA48aNK9INe/bsWf79738TGhqKm5sbgYGBxMXFsXbt2rJ+K27J19eXrl278umnn5ot//TTT4mLi6NOnTrFbvfpp58SHR2Nu7s7/v7+dO3alUOHDhVZb/78+dStW9f0M/7iiy+K3V9eXh4TJ06kXr16pvfcv39/zp49W6b35eIi/yU7O1d7FyBEafj4+NC9e3c+/fRThgwZAhQErIuLCwkJCcyYMcNs/ZycHNq3b8/x48cZN24cUVFRbN68mcmTJ7N3715+/PFHoKC7sEuXLmzbto3Ro0cTExPD1q1b6dixY5EaDh48SMuWLQkLC2PatGkEBQWxatUqkpKSOHfuHGPGjLHKe01NTQUwC4bTp0/TvHlzXFxcGD16NDVr1mT79u1MnDiRtLQ0PvvsMwC+/vprnn/+eV588UXeffddXFxcOHbsGAcPHrzjuoKDg1m5ciWPPvooAwcOZNCgQQCmoO3duze7d+/mrbfeok6dOly8eJHdu3dz/vz5Oz52cQYOHMiDDz7IoUOHqF+/PhcvXmTx4sXMnj272GNOnjyZ119/nZ49ezJ58mTOnz/P2LFjiY2NJSUlhdq1awMFgdq/f386d+7MtGnTuHTpEmPHjiU3N9cs9IxGI507d2bz5s2MGDGCli1b8ueffzJmzBjatWvHL7/8goeHh03eu3BgqhAO7LPPPlMBNSUlRd2wYYMKqAcOHFBVVVVjYmLUfv36qaqqqg0bNlTbtm1r2m7OnDkqoH7zzTdm+3v77bdVQF29erWqqqr6008/qYA6c+ZMs/XeeustFVDHjBljWvbII4+o1apVUy9dumS27gsvvKC6u7urFy5cUFVVVVNTU1VA/eyzz0p8b4Xrvf3226per1dzcnLUvXv3qrGxsWpwcLCamppqWnfIkCGqt7e3+ueff5rt491331UB9bfffjPV4ufnV+Jxx4wZoxb30S/8Xt943LZt25p9X8+ePVvk+1LI29tbHTZsWInHtgZATUxMVI1GoxoREaG+/PLLqqqq6qxZs1Rvb2/18uXL6tSpU83ey99//616eHio8fHxZvtKT09X3dzc1F69eqmqqqoGg0ENCQlRmzRpohqNRtN6aWlpqlarVcPDw03LvvrqKxVQFy1aZLbPlJQUFVBnz55tWnbz97E0bn4PwjlIX4NwGm3btqVmzZp8+umn7N+/n5SUlFt2/a5fvx4vLy+6d+9utrywq3LdunUAbNiwAYCnn37abL1evXqZPc/JyWHdunV07doVT09P8vPzTY/4+HhycnLYsWNHmd7Xq6++ilarxd3dnUaNGnHgwAGWLVtmNtJ0+fLltG/fnpCQELNjF7aoN27cCEDz5s25ePEiPXv25Pvvv+fcuXNlqqksmjdvzvz585k4cSI7duxAr9eXarsb309+fj5qKW/xXNj1vGDBAvLz8/nkk0948skn8fb2LrLu9u3buXbtWpERw6GhoTzwwAOm34fDhw9z6tQpevXqZdZFHh4eXmTg0PLly/Hz8+Pxxx83q79Ro0YEBQXZfTCXsA8JVeE0FEWhf//+fPnll8yZM4c6derQunXrYtc9f/48QUFBRc4dVq5cGVdXV1P34Pnz53F1daVSpUpm6wUFBRXZX35+Pv/5z3/QarVmj/j4eIAyB9jQoUNJSUlhy5YtvPvuu+j1ejp37mzWhZmZmcmyZcuKHLthw4Zmx+7duzeffvopf/75J//617+oXLkyLVq0YM2aNWWqzRILFy6kb9++zJs3j9jYWPz9/enTpw+nT5++5TZpaWlF3lPhHwilUXj+ctKkSezevZuBAwcWu17h9zI4OLjIayEhIWa/D1D051/csszMTC5evIhOpyvyHk6fPl2uf9AIxyHnVIVT6devH6NHj2bOnDm89dZbt1yvUqVK7Ny5E1VVzYL1zJkz5OfnExAQYFovPz+f8+fPmwXrzUFQsWJFNBoNvXv3JjExsdhjRkRElOk9VatWzTQ4KS4ujqCgIJ555hnGjBnDBx98AEBAQABRUVG3fM8hISGmr/v370///v25evUqmzZtYsyYMXTq1IkjR44QHh6Ou7s7ALm5ubi5uZm2u9MQCAgIYMaMGcyYMYP09HR++OEHRo4cyZkzZ1i5cuUt605JSTFbVrdu3VIfMzQ0lIceeohx48ZRt27dW16GUvizzcjIKPLaqVOnzH4foOjPv7hlAQEBVKpU6ZbvzZkv8xJlJ6EqnErVqlV55ZVX+P333+nbt+8t13vwwQf55ptvWLp0KV27djUtLxzF+eCDDwLQvn173nnnHf773/+SlJRkWu9///uf2f48PT1p3749e/bsISoqCp1OZ823Zebpp59m3rx5zJ07l1deeYXw8HA6derEihUrqFmzJhUrVizVfry8vOjYsSN5eXl06dKF3377jfDwcFO38r59+4iJiTGtv2zZstvuszCEr127VuJ6YWFhvPDCC6xbt46tW7fecj2dTmf6g6KsXnrpJTw8POjRo8ct14mNjcXDw4Mvv/zSbL2TJ0+yfv1602mCunXrEhwczFdffcXw4cNNf5D9+eefbNu2zeyPl06dOvH1119jMBho0aLFHb0HcfeQUBVOZ8qUKbddp0+fPsyaNYu+ffuSlpZGZGQkW7ZsYdKkScTHx/PQQw8B0KFDB9q0acOIESO4evUqzZo1Y+vWrSxYsKDIPmfOnEmrVq1o3bo1zz33HNWrV+fy5cscO3aMZcuWsX79equ9x7fffpsWLVowYcIE5s2bx/jx41mzZg0tW7YkKSmJunXrkpOTQ1paGitWrGDOnDlUq1aNwYMH4+HhQVxcHMHBwZw+fZrJkyfj6+trCtD4+Hj8/f0ZOHAg48ePx9XVlfnz53PixInb1lWhQgXCw8P5/vvvefDBB/H39ycgIICKFSvSvn17evXqRb169ahQoQIpKSmsXLmSbt26We37UpwOHTrQoUOHEtfx8/PjzTff5PXXX6dPnz707NmT8+fPM27cONzd3U0jt11cXJgwYQKDBg2ia9euDB48mIsXLzJ27Ngi3b9PPfUU//3vf4mPj2fo0KE0b94crVbLyZMn2bBhA507dzb7g640zp49a+r+3r9/PwA//fQTgYGBBAYG0rZtW4v2J+zA3iOlhCjJjaN/S3Lz6F9VVdXz58+rzz77rBocHKy6urqq4eHh6muvvabm5OSYrXfx4kV1wIABqp+fn+rp6ak+/PDD6u+//17sKNfU1FR1wIABatWqVVWtVqsGBgaqLVu2VCdOnGi2DhaM/p06dWqxr/fo0UN1dXVVjx07pqpqwcjbpKQkNSIiQtVqtaq/v7/atGlTddSoUeqVK1dUVVXVzz//XG3fvr1apUoVVafTqSEhIeqTTz6p7tu3z2zfP//8s9qyZUvVy8tLrVq1qjpmzBh13rx5tx39q6qqunbtWrVx48aqm5ubCqh9+/ZVc3Jy1GeffVaNiopSfXx8VA8PD7Vu3brqmDFj1KtXr5b4fbAU10f/luRWI2fnzZunRkVFqTqdTvX19VU7d+5sGjl983q1a9dWdTqdWqdOHfXTTz9V+/btazb6V1VVVa/Xq++++64aHR2turu7q97e3mq9evXUIUOGqEePHjWtV9rRv4Uj3It7WDp6WNiHoqqlHGonhBBCiBLJ6F8hhBDCSiRUhRBCCCuRUBVCCCGsREJVCCGEsBIJVSGEEMJKJFSFEEIIK5HJH0pgNBo5deoUFSpUKPb+k0IIIe4Nqqpy+fJlQkJCSrzvrYRqCU6dOkVoaKi9yxBCCOEgTpw4QbVq1W75uoRqCQonxD5x4gQ+Pj52rubO6PV6Vq9eTYcOHdBqtfYuRwinIZ8dAZCVlUVoaOhtb5QgoVqCwi5fHx+fuyJUPT098fHxkf8YhLCAfHbEjW53KlAGKgkhhBBWIqEqhBBCWImEqhBCCGElEqpCCCGElUioCiGEEFYio3+FEKIYGZeukXruKtV83SxaPyLAi2Bfj3JbTzgWCVUhhNVYEgSlWdcawVKWfSxMSee1xfsxquCiwJMRCvGlXl9lSpf6PNk4GIz5BQ+DHox6lu/9k/dWHUSjGtApBoa2r87DdQPAqL++jgGCIln4e67Z8Sd3iyQhJqxM71+UL0VVVdXeRTiqrKwsfH19uXTp0l1xneqKFSuIj4+Xa+2c0J2GS1m3t2S7m4OopCAodt1moQWhcj1gluxKZcry/WhUA25KPq88XJP4BgFgyPsnqG78ujCYbvh6d+oZftp34nqI5fNIvUo0DPIsst6N+8jJzWX7kdNoMKDFgKuSjxYDDap4oFOMN2yXD0Y9xnw9V69dwxUDrhjQKgaLfjY3M3iHUPv8VIzqP9dDahSFLSPbS4vVjkqbB9JSFfeU8upSs+ZxLAkroCCYDHnXH/l8vzuVqSsO4Eo+boqB4Q9E8Ei9SteDIe+GUPlnGwx5/PJHJsv3pONKPjrFQMcGlYgsEkgF61/LuYbXb38x29WAlvyCxzIDeXu80SkGs23y9Xm0uXSFnbqCsNKSj+tyA+qPBhT++Ru/K9D1xp7XTdcfFmgCNLnxf7lj1x8lcAfaa4p54Wzx67sAFW47NbiC0UXLNYOCAQ16NOSjQY8rAT5euLu5gaKBs4fQXDmFVtWTi860tUFVSTuXLaHqBCRUizFr1ixmzZqFwXBnf3EKxzovZHE4ARiNN4SNvujXptZRnik4Nh/+i6+3H0ejFoRRQpMqxIRWKGEb8zC78fXc3FzC087wrTYfVwzoyMd1mYH8LVpc1XyzYDPtB/POp85A5xvDacv1x200A5rd2Klx5PqjGB5Ap+KC6FTRRa5AcCnvT2FQFfS4oseVfDR4ebij07mDxhVctKDR3fC1Flxcry/Tcv6ake1pWQUBphYGmSsP31eVYH+f6+trC7bX6Ez7uJirMmnlUfJU1+vBp8GAhkn/iqay303buWg5m20g4ZNd5Kmu6NWCsDQqrvw0/AGC/Lyvr68h89I14qasx3jDj0ejKGwZdL0FajTAeH8AvJUcclWd2XrVAzxL900TdiWhWozExEQSExNNzf075UjBUp5u3cWXX3JQmX2tL906xuJC6p+vc3JyqHLsNF+45qNVCsJJuywf/Q43tGp+0e0K96caLX7frYHWN4bR/uuPMnAD7i9ujP6l0u8jX3UhHw15N4STr/f11tENYfJPWGi5kAsp6ZfRX29N6a8HzAMNqxJU0eemUNOSpVeYvi4V/fXj5KsFwTK2ayMqVvAyO8a5a0b6f7GXPPX6vtGgomXRi22pcj2EMq4YaPXuZgzqP29eoyhsGV76LtC8S9dIKibEHu7YHkrYhx/Q1D2d1xcfwKCq18+pGqgY1RGKOXUSCAzpWqlgfVQ0isKkbvcRFBhgtl6wrweTu0Wa9lu4nun9uGhA4waGXMY+GsGwlReKX084NAlVGytT66gsVLWYQMn9Z1luNhWvHkdJ3wYYiwZRfm4xoXTjvvQ37bfkQMzX59Lq78vs0OX/0x24PB9+zLf+ey8Fd6BdceF03sIdKZrrQaL7p4Wj0ZrC6YrBhaPncgvCSHU1tbKaRARSybeCecuocNsbwqng4WZa52KuyhvLjtwQQAVn7v7TuzkBPl7X19XeYj86Mq7kE/d2ctHW0b9LDqfcS9d4rphAah9ffCD5APUrpJsHRtf7qNis6O96APBM17Ai4VIl5J91g91gUrfoWwdQKdw2xEqQEBNGmzqBpJ3Lpqqvjj1b15d6/eoBnrc8xm3X03nCtVwer+9Hs0bRt92fcDwyUKkEdzpQKeN6d093lw24oUdHPm5KPs+1DqWCq1pCSN0q5G4KyhvXNept8B0oJy5acHX7J5xuCipcry9zcb2+XnFhpCt+++tfX8xVmbjyeEEX3fXWlAEt03o1o1IF7yJhVPzxtVDCfRThn595kRC7g0EmC1PSiwSDJX+YlXX7smyXcelaqYOgNOtasj9r1FScch3kN70hZJ2EweuhalPbHktYRAYqOYDUc1cxqjDO9XM8lLx/XtheHkdXzAJIddGSnWfA09sHRaO7HhSFr7ve0OLRFWznclMrzNXNPGCKC7EbXj93zcigL/eRq7qauh2NN3bxFdbm4grlcAN4PyDGvWhIVLrPur0Gd9I6upXStoKsvX1Ztgv29Sj1/kuzriX7s+U+yo3u+nnTvGz71iHKTELVhiICvHBRYI2xKRoM5BWMcyS+UThenp7Ft4aKPG5oxRXboituXV3B+Zkb5Ov1rC3HS2oCgJ5da5TYxVfe7jSc7HmcOw2Gsm7vVIF0N9BeD1W9hKqzklC1ocJWy/8tTjILFq975CLu8goxS5RXSEgYiTLReRX8m3fVvnWIMpNQtTFHDJbyJOEihAWkper0JFTLgQSLEKJUTOdUpaXqrOQuNUII4Sh03gX/Sqg6LQlVIYRwFNL96/QkVIUQwlHIJTVOT0JVCCEchfb66F+9dP86KwlVIYRwFNJSdXoSqkII4SjknKrTk1AVQghHIZM/OD0JVSGEcBTSUnV6EqpCCOEo5Jyq05NQFUIIRyGjf52ehKoQQjgKOafq9CRUizFr1iwaNGhATEyMvUsRQtxLTKEq3b/OSkK1GImJiRw8eJCUlBR7lyKEuJeYBipdBVW1by2iTCRUhRDCURQOVFKNkJ9r31pEmUioCiGEoygcqARyWY2TklAVQghHoXEFja7gaxms5JQkVIUQwpHIBBBOTUJVCCEciVxW49QkVIUQwpFIS9WpSagKIYQjkakKnZqEqhBCOBKZqtCpSagKIYQjkVmVnJqEqhBCOBJT96+0VJ2RhKoQQjgS6f51ahKqQgjhSGSgklOTUBVCCEcil9Q4NQlVIYRwJDL5g1OTUBVCCEciLVWnJqEqhBCOREb/OjUJVSGEcCSm0b/SUnVGEqpCCOFIZPSvU5NQFUIIR6KT61SdmYSqEEI4Eq1MU+jMJFSFEMKRyEAlpyahKoQQjkQGKjk1CVUhhHAkN7ZUVdW+tQiLSagKIYQjKZz8QTWAIc++tQiLSagKIYQjKRz9C3Je1QlJqAohhCPRaMFFW/C1nFd1Ond9qC5fvpy6detSu3Zt5s2bZ+9yhBDi9mQCCKflau8CbCk/P5/hw4ezYcMGfHx8aNKkCd26dcPf39/epQkhxK1pvSDnkkwA4YTu6pbqzz//TMOGDalatSoVKlQgPj6eVatW2bssIYQomU4mgHBWDh2qmzZt4vHHHyckJARFUVi6dGmRdWbPnk1ERATu7u40bdqUzZs3m147deoUVatWNT2vVq0af/31V3mULoQQZaeT2785K4cO1atXrxIdHc0HH3xQ7OsLFy5k2LBhjBo1ij179tC6dWs6duxIeno6AGox13gpimLTmoUQ4o5p5Ublzsqhz6l27NiRjh073vL16dOnM3DgQAYNGgTAjBkzWLVqFR9++CGTJ0+matWqZi3TkydP0qJFi1vuLzc3l9zcXNPzrKwsAPR6PXq9/k7fjl0V1u/s70OI8maPz47G1QMXIP9aFqp8Zh1CaX/+Dh2qJcnLy2PXrl2MHDnSbHmHDh3Ytm0bAM2bN+fAgQP89ddf+Pj4sGLFCkaPHn3LfU6ePJlx48YVWb569Wo8PT2t+wbsZM2aNfYuQQinVJ6fnZgLWYQAB/emkPqXb7kdV9xadnbpuuJLFapNmjSx6OCKovDDDz+Ync+0tnPnzmEwGKhSpYrZ8ipVqnD69GkAXF1dmTZtGu3bt8doNDJixAgqVap0y32+9tprDB8+3PQ8KyuL0NBQOnTogI+Pj23eSDnR6/WsWbOGhx9+GK1Wa+9yhHAa9vjsaH74ES6m0LBOderHxpfLMUXJCnsub6dUobp3715eeuklvL29b7uuqqpMmTLFrBvVlm4+R6qqqtmyJ554gieeeKJU+3Jzc8PNza3Icq1We9cE0d30XoQoT+X62XEr+L9WY8hFI59Xh1Dan32pu39feeUVKleuXKp1p02bVtrdlllAQAAajcbUKi105syZIq1XIYRwKnL7N6dVqtG/qampBAYGlnqnBw8eJDw8vMxFlYZOp6Np06ZFznOsWbOGli1b2vTYQghhU3L7N6dVqpaqpQEZGhpapmJuduXKFY4dO2Z6npqayt69e/H39ycsLIzhw4fTu3dvmjVrRmxsLB9//DHp6ek8++yzd3TcWbNmMWvWLAwGw52+BSGEsJxMU+i0yjT6Nycnh3379nHmzBmMRqPZa6U9f1kav/zyC+3btzc9LxxE1LdvX+bPn09CQgLnz59n/PjxZGRkcN9997FixYo7biUnJiaSmJhIVlYWvr4y8k4IUc4Kb/8m0xQ6HYtDdeXKlfTp04dz584VeU1RFKu27tq1a1fsBA43ev7553n++eetdkwhhLA73fVBodJSdToWz6j0wgsv0KNHDzIyMjAajWYP6S4VQggrkGkKnZbFoXrmzBmGDx8uI2yFEMJWZJpCp2VxqHbv3p3k5GQblCKEEAKQS2qcmMXnVD/44AN69OjB5s2biYyMLHJBbFJSktWKsxcZ/SuEsCutdP86K4tD9X//+x+rVq3Cw8OD5ORks9mLFEW5K0JVRv8KIexK7qfqtCwO1TfeeIPx48czcuRIXFwc+s5xQgjhnOSSGqdlcSrm5eWRkJAggSqEELZSeE7VmA/5efatRVjE4mTs27cvCxcutEUtQggh4J/RvyCtVSdjcfevwWDgnXfeYdWqVURFRRUZqDR9+nSrFSeEEPckVx24uBa0VPOywaOivSsSpWRxqO7fv5/GjRsDcODAAbPXbr4NmxBCiDLSeUHOJRkB7GQsDtUNGzbYog6HIpfUCCHsTns9VOVaVacio42KkZiYyMGDB0lJSbF3KUKIe5VMVeiUShWq3bp1Iysrq9Q7ffrppzlz5kyZixJCiHueVm7/5oxK1f37/fffc/bs2VLtUFVVli1bxoQJE6hcufIdFSeEEPeswgkgZPSvUylVqKqqSp06dWxdixBCiEJamf/XGZUqVMsyOKlq1aoWbyOEEOI6mVTfKZUqVNu2bWvrOoQQQtyocAIIGajkVGT0rxBCOCKdDFRyRhKqxZg1axYNGjQgJibG3qUIIe5VpnOqV+xbh7CIhGox5DpVIYTdVQgq+DfrL/vWISwioSqEEI6oUq2Cf88fs28dwiISqkII4YhMoXocVNW+tYhSszhUMzMz6d27NyEhIbi6uqLRaMweQgghrMAvHBRNwejfyxn2rkaUksUT6vfr14/09HTefPNNgoOD5c40QghhC646qBgOF/4oaK36hNi7IlEKFofqli1b2Lx5M40aNbJBOUIIIUwq1boeqscgorW9qxGlYHH3b2hoKKr07wshhO351yz4VwYrOQ2LQ3XGjBmMHDmStLQ0G5QjhBDCpFJhqB63bx2i1Czu/k1ISCA7O5uaNWvi6emJVqs1e/3ChQtWK85e5CblQgiHIJfVOB2LQ/W999676wcnJSYmkpiYSFZWFr6+vvYuRwhxryoM1b9TwZAPGov/yxblrEyjf4UQQpQDn6rg6g75OXApHfxr2LsicRsWn1N9+umnmTt3LkeOHLFFPUIIIQq5uNwwWEnOqzoDi0PV29ubadOmUa9ePUJCQujZsydz5szh999/t0V9Qghxb6t0vXUq51WdgsWh+tFHH/H7779z6tQppk+fjq+vLzNnzqRhw4YEBwfbokYhhLh3yWAlp1LmuX8rVKhAxYoVqVixIn5+fri6uhIUFGTN2oQQQkioOhWLQ/XVV1/l/vvvJyAggDfeeIO8vDxee+01MjMz2bNnjy1qFEKIe9eNE+sLh2fx6N+pU6cSGBjImDFj6Ny5M/Xr17dFXUIIIeCfUL10EvTXQOth33pEiSxuqe7Zs4dRo0bx888/06ZNG4KCgkhISODDDz/k0KFDtqhRCCHuXZ6VwN0XUOFCqr2rEbdhcahGR0eTlJTE4sWLOXv2LKtWrcLT05OkpCTuu+8+W9QohBD3LkWROYCdSJmm59izZw/JyckkJyezefNmsrKyaNSoEe3bt7d2fXYh0xQKIRxKpVpwareEqhOwOFQrVqzIlStXiI6Opl27dgwePJg2bdrg4+Nji/rsQqYpFEI4lOvnVTPTfsMYdY1gXzmv6qgsDtUFCxbcdSEqhBCObNslP1oC6Uf3kTBlPZO7RZIQE2bvskQxLD6n2qlTJ1Ognjx5kr/++svqRQkhhCiQcekak3fqAaiunMaowuuLD5Bx6ZqdKxPFsThUjUYj48ePx9fXl/DwcMLCwvDz82PChAkYjUZb1CiEEPes1HNX+cNYMLFOoJKFL1cwqCpp57LtXJkojsXdv6NGjeKTTz5hypQpxMXFoaoqW7duZezYseTk5PDWW2/Zok4hhLgnRQR4cU3x4LgxmJouGbR22c9PakuqB3jauzRRDItD9fPPP2fevHk88cQTpmXR0dFUrVqV559/XkJVCCGsKNjXg8ndIln9fQzPufxAR00KrZ/4twxWclAWh+qFCxeoV69ekeX16tXjwoULVilKCCHEPxJiwjjn8xx89QMd3ffh0qiyvUsSt1CmyR8++OCDIss/+OADoqOjrVKUEEIIcwF1YsGnKi76bPhjg73LEbdgcUv1nXfe4bHHHmPt2rXExsaiKArbtm3jxIkTrFixwhY1CiGEUBSo1wl+/ggOLYe6He1dkSiGxS3Vtm3bcuTIEbp27crFixe5cOEC3bp14/Dhw7Ru3doWNQohhACo/3jBv4d/BEO+fWsRxSrTNIUhISEyIEkIIcpbWCx4+MO1C/DnVqjR1t4ViZuUKlT37dtX6h1GRUWVuRghhBAl0LhCvXjY8yX8vlxC1QGVKlQbNWqEoiioqoqiKKblqqoCmC2TSeiFEMKG6j9REKqHlsOjb4OLxWfxhA2V6qeRmprKH3/8QWpqKosWLSIiIoLZs2ezd+9e9u7dy+zZs6lZsyaLFi2ydb3lYtasWTRo0ICYmBh7lyKEEOYi2oKuAlw+xb6f18t0hQ6mVC3V8PBw09c9evTg/fffJz4+3rQsKiqK0NBQ3nzzTbp06WL1Isub3KVGCOGwtO6kV4ojLGMl25Z/xjvf58oE+w7E4n6D/fv3ExERUWR5REQEBw8etEpRQgghipdx6RqT0wsm4Omu2YRWzZMJ9h2IxaFav359Jk6cSE5OjmlZbm4uEydOpH79+lYtTgghhLnUc1dZbWjKX2olApQsHtdslwn2HYjFl9TMmTOHxx9/nNDQUNMMSr/++iuKorB8+XKrFyiEEOIfEQFeqIqGBfkPM1L7NQM0K1libCsT7DsIi0O1efPmpKam8uWXX/L777+jqioJCQn06tULLy8vW9QohBDiusIJ9t9efJWh6mIauPzJ3LY5MsG+gyjT5A+enp78+9//tnYtQgghSiEhJow2dQK5tHw9Hke/4oGLi4F/2bssQRnOqYaEhNCrVy8+/vhjjhw5YouahBBC3EawrwdBDw8reHJ4BfydZs9yxHUWh+q0adPw8fFh+vTp1KtXj+DgYJ566inmzJnDoUOHbFGjEEKI4lSuBzXag2rkyubZbDt+TkYB25nF3b89e/akZ8+eAGRmZrJhwwaWL1/Oiy++iNFolBmVhBCiPN3/PPyxAeOuBQze1pxriodct2pHZTqneuXKFbZs2cLGjRtJTk5mz549REZG0ratzEMphBDlKaNyHNnGYGq6ZNBPs4pZhi68vvgAbeoEyuAlO7C4+7dFixYEBgby5ptvkp+fz+uvv87p06fZvXs37733ni1qFEIIcQup568xM78bAENcl+PLFblu1Y4sDtWjR4/i6elJjRo1qFGjBrVq1cLPz88GpQkhhLidiAAvflRjOWQMw0fJ5jnXZWgURa5btROLQ/XChQts2LCBuLg41q5dS9u2bQkKCiIhIYE5c+bYokYhhBC3EOzrwaRu0UwzJADQT7OS6fGVpevXTsp0z6CoqCiSkpJYtGgRP/30Ex07dmTx4sUkJiZauz4hhBC3kRATxoRXhpMV2BR3RU/nS/+1d0n3LItDdc+ePbz33nt07twZf39/7r//fvbv38/QoUP54YcfbFGjEEKI2wj288TnsQkFT3Z/ARf+sG9B9yiLR//GxMTQuHFj2rZty+DBg2nTpg0+Pj62qE0IIYQlqsdBrYfg2FpYPxG6f2rviu45FofqhQsXJESFEMJRPTgGjq2DA4ug2QCo3sreFd1TLO7+lUAVQggHFhxVEKYAP74MBr1967nHWByqBoOBd999l+bNmxMUFIS/v7/ZQwghhJ098AYGD384e4isjf+xdzX3FItDddy4cUyfPp0nn3ySS5cuMXz4cLp164aLiwtjx461QYnlb9asWTRo0ICYmBh7lyKEEBZb+NsVXsvqDoDLxrf5YXOKnSu6d1gcqv/973+ZO3cuL7/8Mq6urvTs2ZN58+YxevRoduzYYYsay11iYiIHDx4kJUV+EYUQziXj0jVeW7yfbw1t2GWsjbeSg2b1GzLRfjmxOFRPnz5NZGQkAN7e3ly6dAmATp068eOPP1q3OiGEEBZJPXcVowoqLryp749BVXhMs4O/9620d2n3BItDtVq1amRkZABQq1YtVq9eDUBKSgpubm7WrU4IIYRFIgK8cFEKvj6oVudzwyMA1P15FORcsmNl9waLQ7Vr166sW7cOgKFDh/Lmm29Su3Zt+vTpw4ABA6xeoBBCiNIL9i249ZtGKUjW6YYELnuGorn8F5nfvSzdwDZm8XWqU6ZMMX3dvXt3QkND2bp1K7Vq1eKJJ56wanFCCCEslxATRps6gaSdy6Z6gCcp2+GBHf2pcuwb+r0dQceuz8j9Vm3EopaqXq+nf//+/PHHP9NftWjRguHDh0ugCiGEAwn29SC2ZiUABm1047P8gm7gSdq5TFr8s7RYbcSiUNVqtSxZssRWtQghhLCywoFL7+QnkGasQohygVGaL/hxX4YEqw2U6Zzq0qVLbVCKEEIIayscuHQNd17RD8GoKjzpupFdP80nbsp6Fqak27vEu4rF51Rr1arFhAkT2LZtG02bNsXLy8vs9aSkJKsVJ4QQ4s4UDlx6ffEBUtR6fGh4nETXH5iincu+3Bq8vvgAbeoEyv1XrcTiUJ03bx5+fn7s2rWLXbt2mb2mKIqEqhBCOJjCgUs/7stgyo/daelykMYux5ihm8VTeW+Sdi5bQtVKLA7V1NRUW9QhhBDChoJ9PXgsKphJKw6RpE9khe51YlyOkOS6lOoBD9u7vLuGxedUhRBCOKfCruBTBDFKXzCvwIuuS9Cd3Ma24+dk4JIVlKqlOnz48FLvcPr06WUuRgghhG39cw1rC7J3ncXz4DcYvhnA0Ny3OK9UZHK3SLmG9Q6UKlT37Nlj9nzXrl0YDAbq1q0LwJEjR9BoNDRt2tT6FQohhLCqYF8Pgn09OO07mRMHtlPX5QQf6P7D03mvy8ClO1Sq7t8NGzaYHo8//jjt2rXj5MmT7N69m927d3PixAnat2/PY489Zut6hRBCWMkfl+BZ/TAuqx60cPmdV1wXYlBV0s5l27s0p2XxOdVp06YxefJkKlasaFpWsWJFJk6cyLRp06xanBBCCNuJCPDiT4J5WT8EgCGuP9JRk0L1AE87V+a8LA7VrKwsMjMziyw/c+YMly9ftkpRQgghbK9w4NJatQUf5Rf0NL7v8THBuX/auTLnVaYZlfr37893333HyZMnOXnyJN999x0DBw6kW7dutqhRCCGEjSTEhLFlZHui+r5HbrVYtPlXubbgSU5nnrJ3aU7J4lCdM2cOjz32GM888wzh4eGEh4fz9NNP07FjR2bPnm2LGoUQQthQsK8HsbWr8FO9KZxUA/C4/CfHZvXgm51/3H5jYcbiUPX09GT27NmcP3+ePXv2sHv3bi5cuMDs2bOLTFkohBDCOWRcusbwH/9iUN7LXFXdaOVygCvLXpdrVy1U5skfvLy8iIqKIjo6WsJUCCGcXOHdbH5Xwxiufw6AAa4/cXXHfPsW5mRkRiUhhBCmu9kArDI2Z0Z+wRiZGjve4PyBNXaszLlIqAohhDCNBNYoBck6M/9f/GCIxUXNx/XbPqzYkGzfAp2EhKoQQgjgn5HAH/RsjKIovKIfwi5jbXyVbBomDyYz44S9S3R4EqpCCCFMgn098PfWYVQhFx2D814i3RhIuHIGvn6ajPMX7F2iQ5NQFUIIYebG86sX8KG/fgRZqidVLv3K3hkJfPOz3AL0ViRUhRBCmLn5/OpxtSpD9P9HrupKR83PZP8wkoyLMj9wce6JUO3atSsVK1ake/fu9i5FCCGcQuH51Tceqw/AdmNDXtY/C0A/15XkbnrfnuU5rHsiVJOSkvjiiy/sXYYQQjiVYF8PHosKNnUFLzO2ZKL+aQCq754M+7+zY3WO6Z4I1fbt21OhQgV7lyGEEE7n5q7gz4yPcTjimYIXlwzh/K8r2Hb8nMy8dF2pblJuS5s2bWLq1Kns2rWLjIwMlixZQpcuXczWmT17NlOnTiUjI4OGDRsyY8YMWrdubZ+Ci2EwGNDr9fYuo0R6vR5XV1dycnIwGAz2Lkc4IK1Wi0ajsXcZwgElxITRpk4gaeeyqR7gSXCFjrD4GhxYhMfifkzNe51fqc3kbpEkxITZu1y7snuoXr16lejoaPr378+//vWvIq8vXLiQYcOGMXv2bOLi4vjoo4/o2LEjBw8eJCys4IfXtGlTcnNzi2y7evVqQkJCbFa7qqqcPn2aixcv2uwY1qKqKkFBQZw4cQLl+l+cQtzMz8+PoKAg+R0RRQT7ehDs62F6nvHAexzZd5S2Lvv4VDeVHnmjeX2xQps6gWbr3WvsHqodO3akY8eOt3x9+vTpDBw4kEGDBgEwY8YMVq1axYcffsjkyZMB2LVrl1Vqyc3NNQvnrKwsoKCVV1xLNDMzk6ysLAIDA/H09HTo/4hUVeXq1at4eXk5dJ3CPlRVJTs7m7Nnz2IwGKhSpYq9S3IYhZ99R++NKm/HzubwbN4w/qebRGOXYyzQTaF77hiOZ2YR4Gn3aLG60v78Hfqd5+XlsWvXLkaOHGm2vEOHDmzbts3qx5s8eTLjxo0rsnz16tV4enqaLVMUheDgYIKCgtBqtU7xgdPpdE5Rp7APrVZLhQoVyMjIYPfu3aiqau+SHMqaNTL/7Y0u5kIObvTPe4VvdeOp7fIXX+om8fMvCisO+dm7PKvLzi7dJUQOHarnzp0r9q/mKlWqcPr06VLv55FHHmH37t1cvXqVatWqsWTJEmJiYoqs99prrzF8+HDT86ysLEJDQ+nQoQM+Pj5m6+bm5pKeno6/vz8eHo7f1aGqKpcvX6ZChQrSUhW3pNVquXz5Mg888ABubm72Lsch6PV61qxZw8MPP4xWq7V3OQ5FG3aSN74/yDN5r/Gd2zgiXDKp/veH5Hf6ATwq2rs8qyrsubwdhw7VQjeHgKqqFgXDqlWrSrWem5tbsf+RaLXaIh8mg8GAoihoNBpcXBx/ELXRaAQKvpfOUK+wD41Gg6IouLq6SoDcpLj/B+51ve6PoH39INLOZePm2gy+7YJy9hDahU9Bn+/B7e656qK0P3uH/t81ICAAjUZTpFV65swZOecjhBAOINjXg9ialagcXh/6LC1oof61C/73FOTde7MuOXSo6nQ6mjZtWuRcxpo1a2jZsqWdqrr7JScnoyiKU4xqLk/9+vUrcrmXM6hevTozZswwPVcUhaVLl9qtHnEXq1wfnlkMugrw5xZY+DToc+xdVbmye/fvlStXOHbsmOl5amoqe/fuxd/fn7CwMIYPH07v3r1p1qwZsbGxfPzxx6Snp/Pss8/arKZZs2Yxa9YsuZ5TmJk5c6YM3hHidqo2gae/xbigKy7H15Pzv2dwf/p/4Kqzd2Xlwu4t1V9++YXGjRvTuHFjAIYPH07jxo0ZPXo0AAkJCcyYMYPx48fTqFEjNm3axIoVKwgPD7dZTYmJiRw8eJCUlBSbHUM4H19fX/z8/OxdhhAOb+GZqvTOHk6OqsU9dQ0n5vUCQ769yyoXdg/Vdu3aoapqkcf8+fNN6zz//POkpaWRm5vLrl27aNOmjf0KtpGMS9fKbaqv3NxckpKSqFy5Mu7u7rRq1arYPyC2bt1KdHQ07u7utGjRgv3795te+/PPP3n88cepWLEiXl5eNGzYkBUrVtxRXYXdq+PGjaNy5cr4+PgwZMgQ8vLyTOvc3JUJ0KhRI8aOHWt6rigKH330EZ06dcLT05P69euzfft2jh07Rrt27fDy8iI2Npbjx4+bthk7diyNGjXio48+IjQ0FE9PT3r06GHWBX5z92+7du1ISkpixIgR+Pv7ExQUZFYHwO+//06rVq1wd3enQYMGrF27tsTu12XLluHn52caWLZ3796Cm0W/8oppnSFDhtCzZ0/T823bttGmTRs8PDwIDQ0lKSmJq1ev3ua7LYRtZFy6xmuL97PV2JAh+uHkqq6Enl7DtW8G3hPBavdQFbAwJZ24KevpNXcncVPWszAl3abHGzFiBIsWLeLzzz9n9+7d1KpVi0ceeYQLF8xvPvzKK6/w7rvvkpKSQuXKlXniiSdM17kmJiaSm5vLpk2b2L9/P2+//Tbe3t53XNu6des4dOgQGzZs4KuvvmLJkiXFXjt8OxMmTKBPnz7s3buXevXq0atXL4YMGcJrr73GL7/8AsALL7xgts2xY8f45ptvWLZsGStXrmTv3r0kJiaWeJzPP/8cLy8vdu7cyTvvvMP48eNNYwCMRiNdunTB09OTnTt38vHHHzNq1KgS99emTRsuX77Mnj17ANi4cSMBAQFs3LjRtE5ycjJt27YFYP/+/TzyyCN069aNffv2sXDhQrZs2VLkvQlRXlLPXcV4/SzJRmM0ifqh6FUNHoeXkv3tELYdzbyr5wmWULWzwr/qCn8JjSq8vviAzX7prl69yocffsjUqVPp2LEjDRo0YO7cuXh4ePDJJ5+YrTtmzBgefvhhIiMj+fzzz8nMzGTJkiUApKenExcXR2RkJDVq1KBTp05W6UHQ6XR8+umnNGzYkMcee4zx48fz/vvvm1pupdW/f3+efPJJ6tSpw6uvvkpaWhpPP/00jzzyCPXr12fo0KEkJyebbZOTk8Pnn39Oo0aNaNOmDf/5z3/4+uuvS7wmOioqijFjxlC7dm369OlDs2bNWLduHVAwacjx48f54osviI6OplWrVrz11lsl1u3r60ujRo1MtSUnJ/N///d//Prrr1y+fJnTp09z5MgR2rVrB8DUqVPp1asXw4YNo3bt2rRs2ZL333+fL774gpyce2uAiHAMN97gHGCtsSlJ+UMxKho8f/+OU18MotWUtTZvPNiLhKqd3fhXXSGDqpJ2zjZD0Y8fP45erycuLs60TKvV0rx5cw4dOmS2bmxsrOlrf39/6tata1onKSmJiRMnEhcXx5gxY9i3b98tjzlp0iS8vb1Nj/T0W3+YoqOjzWavio2N5cqVK5w4ccKi9xkVFWX6uvDyq8jISLNlOTk5Zhd0h4WFUa1aNbNjG41GDh8+XKrjAAQHB3PmzBkADh8+TGhoKEFBQabXmzdvftva27VrR3JyMqqqsnnzZjp37sx9993Hli1b2LBhA1WqVKFevXpAwRSd8+fPN/v+PvLIIxiNRlJTU297LCGs7ea72mgUhUYdnuHFvETyVRe6azbxlmYeoxbvuytbrHYf/euIynP0b+FfdTcGq0ZRqB7geeuN7kDh6NWyTqhRuM6gQYN45JFH+PHHH1m9ejWTJ09m2rRpvPjii0W2efbZZ3nyySdNz8tyk4PC47q4uBQZgVvc1Is3XqhduG1xy0pqAReuU9L35eYLwhVFMe3T0klKCrVr145PPvmEX3/9FRcXFxo0aEDbtm3ZuHEjf//9t6nrt7D+IUOGkJSUVGQ/hTecEKK83XxXm9RzV5lsuB9FVZmp/YCnXJMBSDvbAihoXEQEeN0VE/FLS7UY5Tn6t7i/6iZ1u89mv1y1atVCp9OxZcsW0zK9Xs8vv/xC/fr1zdbdsWOH6eu///6bI0eOmFpIAKGhoTz77LMsXryYl156iblz5xZ7TH9/f2rVqmV6uLre+m+5X3/9lWvX/vnrdceOHXh7e5takIGBgWRkZJhez8rKslqLLD09nVOnTpmeb9++HRcXF+rUqVOm/dWrV4/09HQyMzNNy0rzO1V4XnXGjBm0bdsWRVFo27YtycnJZudTAZo0acJvv/1m9v0tfOh098YlDMIxFU4KEezrYWo8LDfG8n/6RAyqwlOuyVRc9xKtpqwtt/Ek5UFC1QEkxISxZWR7vhp8P1tGtrfp/Qi9vLx47rnneOWVV1i5ciUHDx5k8ODBZGdnM3DgQLN1x48fz7p16zhw4AD9+vUjICDANPp12LBhrFq1itTUVHbv3s369euLhHJZ5OXlMXDgQA4ePMhPP/3EmDFjeOGFF0xTKz7wwAMsWLCAzZs3c+DAAfr27Wu1e4C6u7vTt29ffv31VzZv3kxSUhJPPvmkWfetJR5++GFq1qxJ37592bdvH1u3bjUNVCqpBVt4XvXLL780nTtt06YNu3fvNjufCvDqq6+yfft2EhMT2bt3L0ePHuWHH34otsdACHu5sfHwg7ElL+UnYsSFehnfM1kzFwWjzceTlBfp/nUQN9+r0JamTJmC0Wikd+/eXL58mWbNmrFq1SoqVqxYZL2hQ4dy9OhRoqOj+eGHH0ytH4PBQGJiIidPnsTHx4dHH32U9957745re/DBB6lduzZt2rQhNzeXp556yuwylddee40//viDTp064evry4QJE6zWUq1VqxbdunUjPj6eCxcuEB8fz+zZs8u8P41Gw9KlSxk0aBAxMTHUqFGDqVOn8vjjj+Pu7l7itu3bt2f37t2mAK1YsSINGjTg1KlTZn+8REVFsXHjRkaNGkXr1q1RVZWaNWuSkJBQ5rqFsAXzLuEHOJpSi1qb/48nXTfioqiM0P8bg+pC2rlsp+4GVlSZIuaWsrKy8PX15dKlS0XuUpOTk0NqaioRERG3/Q/SERiNRrKysvDx8XHYCfX79evHxYsX7TKF3tixY1m6dCl79+616XG2bt1Kq1atOHbsGDVr1rTpscrC2X6vy4Ner2fFihXEx8fLhPpWlHHpGpPeeYv3XGfhqhhZbGjFq/nPsWnkgw4ZqiXlwY2kpSqEDS1ZsgRvb29q167NsWPHGDp0KHFxcQ4ZqEKUp2BfD1p1+TdDl2qY4fofumm20DS0Ahjbsu34OacduCShKoQNXb58mREjRnDixAkCAgJ46KGHmDZtmr3LEsIhFHQJj+SP3XWpszmJ8FM/8eP0fzFUn4hRcWVyt0gSYsLIuHTNaUYIS/dvMW68pObIkSPS/SvuGc72e10epPu3fFzY/T1e3w/ATclntaEpL+iTMCg6RnSsy9s//Y5RBRcFU9CWt9J2/8r/rsWQCfWFEKJ8/e4bd32uYC0dNLv4SDsdVzWXKdcDFWw/45w1SKgKIYSwu4gALzapjeivf4Vrqo72ml/5RPsubmqu2Xq2nHHOGiRUhRBC2F3htaw71Uj65r3KFdWdVpoDfKF7Gy/+aZlqFAVPnUu53dXLUjJQSQghhEP451rWFuTktcB7aS+a5/7Of3WT6JP3KleVCnRpHELX2dvsfo71VqSlKoQQwmEUTm8YUL819P0BPCrSyOU424Nn8H3/OizZ85dDn2OVUBVCCOGYQhpDvx/BKxCvvw9S48enCFD/NlvF0c6xSqgWY9asWTRo0ICYmBh7l2J37dq1Y9iwYfYuwyrS0tJQFMXmsyYJIayoSkPo/xNUCMHz0lG+cRtPVc6aXrblXb3KQkK1GHf7JTXFBWVycjKKonDx4kWz5YsXL2bChAnlV1w5utV7FkI4mIDa0H8F+IVTXcnkW7fxRCgZprt6AQ4zcEkGKokS+fv7l8tx9Hq9XFgvhLg1/4iCFusXnQk5f5RVvpO51ONb1l+AuCnrHWbgkrRU7zH9+/dn48aNzJw5E0VRUBSFtLQ02rdvDxTcDUVRFPr16wcUbdVWr16diRMn0qdPH7y9vQkPD+f777/n7NmzdO7cGW9vbyIjI/nll19KrENRFObMmUPnzp3x8vJi4sSJACxbtoymTZvi7u5OjRo1GDduHPn5+abtxo4dS1hYGG5uboSEhJjdnFtRlCKT8fv5+TF//vwixy/pPX/33XdERkbi4eFBpUqVeOihh7h69Wppvr1CCFvyrVoQrFUi0eWco9K33fhqyVKHGrgkoWpNqgp5V+3zKOVskzNmzCA2NpbBgweTkZFBRkYGoaGhLFq0CIDDhw+TkZHBzJkzb7mP9957j7i4OPbs2cNjjz1G79696dOnD8888wy7d++mVq1a9OnTh9vNgDlmzBg6d+7M/v37GTBgAKtWreKZZ54hKSmJgwcP8tFHHzF//nzeeustoCDs3nvvPT766COOHj3K0qVLiYyMLOUPx9yt3nNGRgY9e/ZkwIABHDp0iOTkZLp163bb9yKEKCfegdBvGVSLwSXnb77UvkUL5ZDpZXsPXJLuX2vSZ8OkEPsc+/VToPO67Wq+vr7odDo8PT3Nbr5d2M1buXJl/Pz8StxHfHw8Q4YMAWD06NF8+OGHxMTE0KNHD6DgxtmxsbFkZmaWeIPvXr16MWDAANPz3r17M3LkSPr27QtAjRo1mDBhAiNGjGDMmDGkp6cTFBTEQw89hFarJSwsjObNm9/2PRdHo9EU+56PHz9Ofn4+3bp1Izw8HKDMwS2EsBGPitB7KblfJuB9Yguf66bwrH4YycbGdh+4JC1VYbGoqCjT11WqVAHMg6dw2ZkzZ0rcT7Nmzcye79q1i/Hjx+Pt7W16FLaos7Oz6dGjB9euXaNGjRoMHjyYJUuWmHUNW0N0dDQPPvggkZGR9OjRg7lz5/L333/ffkMhRPly88atzyL+qtwOd0XPXO10HtfssPvAJWmpWpPWs6DFaK9jl9ehbhhQpCjKLZcZjcYS9+PlZd6yNhqNjBs3jm7duhVZ193dndDQUA4fPsyaNWtYu3Ytzz//PFOnTmXjxo1otVoURSnSTavX6y16bxqNhjVr1rBt2zZWr17Nf/7zH0aNGsXOnTuJiIiwaF9CCBvTulN1yHdc+2YwHoeX8L7uA1LSKxG3uLbdBi5JS9WaFKWgC9Yej+tBVho6nQ6DwVBkGVBkeXlq0qQJhw8fplatWkUehber8/Dw4IknnuD9998nOTmZ7du3s3//fgACAwPJyMgw7e/o0aNkZ9/63Mqt3rOiKMTFxTFu3Dj27NmDTqdjyZIl1n67Qghr0GjxSPgEmvRFUY003z+Gfi4/AfYZuCQt1WLceD/Vu1H16tXZuXMnaWlpeHt74+/vT3h4OIqisHz5cuLj4/Hw8MDb27tc6xo9ejSdOnUiNDSUHj164OLiwr59+9i/fz8TJ05k/vz5GAwGWrRogaenJwsWLMDDw8N07vOBBx7ggw8+4P7778doNPLqq6+WeJlOce/5t99+Y926dXTo0IHKlSuzc+dOzp49S/369cvr2yCEsJSLBh6fyV85WqoenMdo7QIqkM1MQzcMKqSdyy63m5tLS7UYd/vkDy+//DIajYYGDRoQGBhIeno6VatWZdy4cYwcOZIqVarwwgsvlHtdjzzyCMuXL2fNmjXExMRw//33M336dFNo+vn5MXfuXOLi4oiKimLdunUsW7aMSpUqATBt2jRCQ0Np06YNvXr14uWXX8bT89bd4sW9Zx8fHzZt2kR8fDx16tThjTfeYNq0aXTs2LFcvgdCiDJSFFw6TGBafsGAyf/TLuJN1y9xVdRyHbikqHKtwC2VdKf3nJwcUlNTiYiIwN3d3U4Vlp7RaCQrKwsfHx9TV6oQN3O23+vyoNfrWbFiBfHx8TJBiRNYmJLO70vfZYz2cwD+qNaFGv0/Ac2ddcyWlAc3kv9dhRBC3DUSYsL494i3OdryXVTFhRonl8J3/SE/97bbWoOEqhBCiLtKsK8HtTsMRnnyC1SNDg79wNXlr5XLsSVUhRBC3JUWXommz7WXOGgMp/3OpixMSbf5MSVUhRBC3HUyLl3jtcX72WyM5LG8tzij+pXL5TUSqndIxnmJu4n8Pou7Req5q6aJ9tXrUVce8wJLqJZR4SjAkiYXEMLZFP4+yyhX4ewiArxwuWlOnPKYF1gmfygjjUaDn5+faX5bT09P0/R8jshoNJKXl0dOTo5cUiOKUFWV7Oxszpw5g5+fHxqNxt4lCXFHgn09mNwtktcXH8CgqqYbmtt6EggJ1TtQeAeW200c7whUVeXatWt4eHg4dPgL+/Lz8yvxzkJCOJOEmDDa1Akk7Vw21QM8y2VWJQnVYpR2mkJFUQgODqZy5coWT9xe3vR6PZs2baJNmzbStSeKpdVqpYUq7jrBvh7lNkUhSKgWKzExkcTERNMMGrej0Wgc/j8jjUZDfn4+7u7uEqpCCGEjcnJNCCGEsBIJVSGEEMJKJFSFEEIIK5FzqiUovBA+KyvLzpXcOb1eT3Z2NllZWXJOVQgLyGdHwD85cLsJUiRUS3D58mUAQkND7VyJEEIIR3D58uUSB7DK/VRLYDQaOXXqFBUqVDC7tjMmJqZMNzC3dLvSrl+a9bKysggNDeXEiRMl3gvwblTWn5ctlVdN1j7One7P0T47pVn3Xv7sgON9fuz12VFVlcuXLxMSElLiBDrSUi2Bi4sL1apVK7Jco9GU6cNl6XalXd+S/fr4+Nxz/zGU9edlS+VVk7WPc6f7c7TPjiXr3oufHXC8z489PzulucRSBiqVQWJiYrlsV9r1y1rPvcIRvz/lVZO1j3On+3O0z05Z9n2vcbTvj6N/dqT79x5ROJHFpUuXHOqvTiEcnXx2hCWkpXqPcHNzY8yYMbi5udm7FCGcinx2hCWkpSqEEEJYibRUhRBCCCuRUBVCCCGsREJVCCGEsBIJVSGEEMJKJFSFEEIIK5FQFWYuX75MTEwMjRo1IjIykrlz59q7JCGcxokTJ2jXrh0NGjQgKiqKb7/91t4liXIml9QIMwaDgdzcXDw9PcnOzua+++4jJSWFSpUq2bs0IRxeRkYGmZmZNGrUiDNnztCkSRMOHz6Ml5eXvUsT5UTm/hVmNBoNnp6eAOTk5GAwGG57qyMhRIHg4GCCg4MBqFy5Mv7+/ly4cEFC9R4i3b9OZtOmTTz++OOEhISgKApLly4tss7s2bOJiIjA3d2dpk2bsnnzZouOcfHiRaKjo6lWrRojRowgICDAStULYV/l8fkp9Msvv2A0GuXWkfcYCVUnc/XqVaKjo/nggw+KfX3hwoUMGzaMUaNGsWfPHlq3bk3Hjh1JT083rdO0aVPuu+++Io9Tp04B4Ofnx6+//kpqair/+9//yMzMLJf3JoStlcfnB+D8+fP06dOHjz/+2ObvSTgWOafqxBRFYcmSJXTp0sW0rEWLFjRp0oQPP/zQtKx+/fp06dKFyZMnW3yM5557jgceeIAePXpYo2QhHIatPj+5ubk8/PDDDB48mN69e1u7bOHgpKV6F8nLy2PXrl106NDBbHmHDh3Ytm1bqfaRmZlJVlYWUHB3jk2bNlG3bl2r1yqEo7HG50dVVfr168cDDzwggXqPkoFKd5Fz585hMBioUqWK2fIqVapw+vTpUu3j5MmTDBw4EFVVUVWVF154gaioKFuUK4RDscbnZ+vWrSxcuJCoqCjT+doFCxYQGRlp7XKFg5JQvQspimL2XFXVIstupWnTpuzdu9cGVQnhHO7k89OqVSuMRqMtyhJOQrp/7yIBAQFoNJoif1WfOXOmyF/fQghz8vkR1iChehfR6XQ0bdqUNWvWmC1fs2YNLVu2tFNVQjgH+fwIa5DuXydz5coVjh07ZnqemprK3r178ff3JywsjOHDh9O7d2+aNWtGbGwsH3/8Menp6Tz77LN2rFoIxyCfH2FzqnAqGzZsUIEij759+5rWmTVrlhoeHq7qdDq1SZMm6saNG+1XsBAORD4/wtbkOlUhhBDCSuScqhBCCGElEqpCCCGElUioCiGEEFYioSqEEEJYiYSqEEIIYSUSqkIIIYSVSKgKIYQQViKhKoQQQliJhKoQd5Hk5GQUReHixYvlfmxFUVAUBT8/vxLXGzt2LI0aNTJ7XrjtjBkzbFqjELYmoSqEk2rXrh3Dhg0zW9ayZUsyMjLw9fW1S02fffYZR44csWibl19+mYyMDKpVq2ajqoQoPzKhvhB3EZ1OR1BQkN2O7+fnR+XKlS3axtvbG29vbzQajY2qEqL8SEtVCCfUr18/Nm7cyMyZM01dp2lpaUW6f+fPn4+fnx/Lly+nbt26eHp60r17d65evcrnn39O9erVqVixIi+++CIGg8G0/7y8PEaMGEHVqlXx8vKiRYsWJCcnl6nWKVOmUKVKFSpUqMDAgQPJycmxwndACMckoSqEE5o5cyaxsbEMHjyYjIwMMjIyCA0NLXbd7Oxs3n//fb7++mtWrlxJcnIy3bp1Y8WKFaxYsYIFCxbw8ccf891335m26d+/P1u3buXrr79m37599OjRg0cffZSjR49aVOc333zDmDFjeOutt/jll18IDg5m9uzZd/TehXBk0v0rhBPy9fVFp9Ph6el52+5evV7Phx9+SM2aNQHo3r07CxYsIDMzE29vbxo0aED79u3ZsGEDCQkJHD9+nK+++oqTJ08SEhICFJz3XLlyJZ999hmTJk0qdZ0zZsxgwIABDBo0CICJEyeydu1aaa2Ku5a0VIW4y3l6epoCFaBKlSpUr14db29vs2VnzpwBYPfu3aiqSp06dUznO729vdm4cSPHjx+36NiHDh0iNjbWbNnNz4W4m0hLVYi7nFarNXuuKEqxy4xGIwBGoxGNRsOuXbuKDB66MYiFEEVJqArhpHQ6ndngImtp3LgxBoOBM2fO0Lp16zvaV/369dmxYwd9+vQxLduxY8edliiEw5JQFcJJVa9enZ07d5KWloa3tzf+/v5W2W+dOnV4+umn6dOnD9OmTaNx48acO3eO9evXExkZSXx8fKn3NXToUPr27UuzZs1o1aoV//3vf/ntt9+oUaOGVWoVwtHIOVUhnNTLL7+MRqOhQYMGBAYGkp6ebrV9f/bZZ/Tp04eXXnqJunXr8sQTT7Bz585bjjC+lYSEBEaPHs2rr75K06ZN+fPPP3nuueesVqcQjkZRVVW1dxFCCOenKApLliyhS5cuZdq+evXqDBs2rMgsUUI4E2mpCiGspmfPnhZPNzhp0iS8vb2t2tIWwl6kpSqEsIpjx44BoNFoiIiIKPV2Fy5c4MKFCwAEBgbabd5iIaxBQlUIIYSwEun+FUIIIaxEQlUIIYSwEglVIYQQwkokVIUQQggrkVAVQgghrERCVQghhLASCVUhhBDCSiRUhRBCCCuRUBVCCCGs5P8B+eVSRea42hIAAAAASUVORK5CYII=",
      "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",
    "hm1 = w1.headinside(tm)\n",
    "plt.loglog(to, -ho, \".\", label=\"obs - pumping well\")\n",
    "plt.loglog(tm, -hm1[0], label=\"ttim results\")\n",
    "plt.xlabel(\"time [d]\")\n",
    "plt.ylabel(\"drawdown [m]\")\n",
    "plt.legend()\n",
    "plt.title(\"Model Results - Model 1\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall a good fit has been observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing single layer model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative of simpler model, we will test the assumption that we can neglect the effect of the semi-confining top layer and the underlying aquifer. Thus, we simulate the aquifer as confined.\n",
    "\n",
    "Therefore this model is a single layer semi-confined aquifer that only includes the top aquitard and the first aquifer. In the well, we are not considering the ```rc``` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self.neq  1\n",
      "solution complete\n"
     ]
    }
   ],
   "source": [
    "ml0 = ttm.ModelMaq(kaq=50, z=[zt, zb], Saq=1e-4, topboundary=\"conf\", tmin=1e-4, tmax=0.04)\n",
    "w0 = ttm.Well(ml0, xw=0, yw=0, rw=rw, res=1, tsandQ=[(0, Q), (t0, 0)], layers=0)\n",
    "ml0.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibration\n",
    "\n",
    "In the calibration we repeat the procedure for model 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..............................................................................................................................................................................................................................................................................................................................................................................\n",
      "Fit succeeded.\n"
     ]
    }
   ],
   "source": [
    "ca0 = ttm.Calibrate(ml0)\n",
    "ca0.set_parameter(name=\"kaq0\", initial=50, pmin=0)\n",
    "ca0.set_parameter(name=\"Saq0\", initial=1e-4, pmin=0)\n",
    "# ca0.set_parameter(name=\"c0\", initial=1000, pmin=10, pmax=10000)\n",
    "ca0.set_parameter_by_reference(name=\"res\", parameter=w0.res[:], initial=1)\n",
    "ca0.seriesinwell(name=\"obs\", element=w0, t=to, h=ho)\n",
    "ca0.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "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>kaq0</th>\n",
       "      <td>48.653125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq0</th>\n",
       "      <td>0.000076</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>res</th>\n",
       "      <td>0.022092</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        optimal\n",
       "kaq0  48.653125\n",
       "Saq0   0.000076\n",
       "res    0.022092"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RMSE: 0.00873670055844304\n"
     ]
    }
   ],
   "source": [
    "display(ca0.parameters.loc[:, [\"optimal\"]])\n",
    "print(\"RMSE:\", ca0.rmse())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdUAAAE/CAYAAAAQZlkTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABNtElEQVR4nO3deVhUZf/H8fdhmGEVEEEBRcR9CXBDQ1wrNclyeTTScs/HikJ/VmZZ4pZapmmlmVpZ1lNWLqWZSyruGrmkprkFoYkrKSrbMHN+fyCTI4gMDsyMfl/XNZdw5izfGRg/3Pe5z30UVVVVhBBCCHHHnGxdgBBCCHG3kFAVQgghrERCVQghhLASCVUhhBDCSiRUhRBCCCuRUBVCCCGsREJVCCGEsBIJVSGEEMJKJFSFEEIIK5FQFXZt4cKFKIqCoigkJiYWel5VVWrXro2iKLRv396qx1YUhXHjxlm8XUpKCoqisHDhwhKtV/BwcnKiYsWKPPjgg6xdu7Z0RVtZ+/btzd7XzMxMxo0bV+TPorwUvF8DBw4s8vkJEyaY1klJSbHacQcOHEiNGjVKte3N72NRDAYDM2bM4OGHH6ZatWq4u7vToEEDRo8ezaVLl0p1XFH+JFSFQ6hQoQIff/xxoeWbNm3ixIkTVKhQwQZVWccLL7zAjh072LJlC++88w7Hjh0jJiaGzZs327q0QjIzMxk/frxNQxXyfx++/fZbrly5YrZcVVUWLlyIl5eXjSorvaysLMaNG0dISAgzZ85k1apVDB06lHnz5hEdHU1WVpatSxQlIKEqHEJsbCxLliwhIyPDbPnHH39MVFQU1atXt1Fld6569ercf//9REdHM2TIEL744gsMBkORf0SIfN26dUNVVb7++muz5Rs2bCA5OZnY2FgbVVZ6bm5uJCcn89FHH9GrVy/at2/PyJEjmTdvHocOHWLJkiW2LlGUgISqcAh9+vQB4KuvvjItu3z5MkuWLGHw4MFFbpOens5zzz1H1apV0el01KxZkzFjxpCTk2O2XkZGBkOHDqVSpUp4enry8MMPc/To0SL3eezYMfr27UvlypVxcXGhQYMGzJ4920qvMl/z5s0BOHv2rNnyM2fOMGzYMKpVq4ZOpyM0NJTx48eTl5dntt6HH35IREQEnp6eVKhQgfr16/Paa6+Znh83bhyKohQ6bkFX+626TFNSUvD39wdg/Pjxhbphz58/z3//+1+Cg4NxcXHB39+f6Ohofv7559K+Fbfk7e1Njx49+OSTT8yWf/LJJ0RHR1O3bt0it/vkk0+IiIjA1dUVX19fevToweHDhwutt3DhQurVq2f6GX/++edF7i83N5dJkyZRv35902seNGgQ58+ft/g1aTQaKlWqVGh5ixYtADh58qTF+xTlz9nWBQhREl5eXvTq1YtPPvmEYcOGAfkB6+TkRGxsLDNnzjRbPzs7mw4dOnDixAnGjx9PeHg4W7ZsYcqUKezbt48ff/wRyO8u7N69O9u3b2fs2LFERkaybds2unTpUqiGQ4cO0apVK6pXr8706dMJCAhgzZo1xMfHc+HCBRISEqzyWpOTkwHMguHMmTO0aNECJycnxo4dS61atdixYweTJk0iJSWFTz/9FICvv/6a5557jhdeeIF33nkHJycnjh8/zqFDh+64rsDAQFavXs3DDz/MkCFDePrppwFMQduvXz/27NnDm2++Sd26dbl06RJ79uzh4sWLd3zsogwZMoQHH3yQw4cP06BBAy5dusTSpUuZM2dOkcecMmUKr732Gn369GHKlClcvHiRcePGERUVRVJSEnXq1AHyA3XQoEF069aN6dOnc/nyZcaNG0dOTg5OTv+2Q4xGI926dWPLli2MGjWKVq1a8ddff5GQkED79u359ddfcXNzu+PXuWHDBgAaNWp0x/sS5UAVwo59+umnKqAmJSWpGzduVAH14MGDqqqqamRkpDpw4EBVVVW1UaNGart27UzbzZ07VwXUb775xmx/b731lgqoa9euVVVVVX/66ScVUGfNmmW23ptvvqkCakJCgmlZ586d1WrVqqmXL182W/f5559XXV1d1fT0dFVVVTU5OVkF1E8//bTY11aw3ltvvaXq9Xo1Oztb3bdvnxoVFaUGBgaqycnJpnWHDRumenp6qn/99ZfZPt555x0VUH///XdTLT4+PsUeNyEhQS3qo1/wXt943Hbt2pm9r+fPny/0vhTw9PRUR4wYUeyxrQFQ4+LiVKPRqIaGhqovvfSSqqqqOnv2bNXT01O9cuWKOm3aNLPX8s8//6hubm5qTEyM2b5SU1NVFxcXtW/fvqqqqqrBYFCDgoLUpk2bqkaj0bReSkqKqtVq1ZCQENOyr776SgXUJUuWmO0zKSlJBdQ5c+aYlt38PpbUqVOn1CpVqqjNmzdXDQaDxduL8ifdv8JhtGvXjlq1avHJJ59w4MABkpKSbtn1u2HDBjw8POjVq5fZ8oKuyvXr1wOwceNGAJ588kmz9fr27Wv2fXZ2NuvXr6dHjx64u7uTl5dnesTExJCdnc3OnTtL9bpeeeUVtFotrq6uNG7cmIMHD7JixQqzkaYrV66kQ4cOBAUFmR27oEW9adMmIL+r8NKlS/Tp04fvv/+eCxculKqm0mjRogULFy5k0qRJ7Ny5E71eX6Ltbnw9eXl5qCW8xXNB1/OiRYvIy8vj448/5vHHH8fT07PQujt27CArK6vQiOHg4GAeeOAB0+/DkSNHOH36NH379jXrIg8JCaFVq1Zm265cuRIfHx8effRRs/obN25MQEDAHQ/mSk9PJyYmBlVVWbx4sVkrWdgv+SkJh6EoCoMGDeKLL75g7ty51K1blzZt2hS57sWLFwkICCh07rBy5co4OzubugcvXryIs7NzoXNZAQEBhfaXl5fH+++/j1arNXvExMQAlDrAhg8fTlJSElu3buWdd95Br9fTrVs3sy7Ms2fPsmLFikLHLugSLDh2v379+OSTT/jrr7/4z3/+Q+XKlWnZsiXr1q0rVW2WWLx4MQMGDGDBggVERUXh6+tL//79OXPmzC23SUlJKfSaCv5AKImC85eTJ09mz549DBkypMj1Ct7LwMDAQs8FBQWZ/T5A4Z9/UcvOnj3LpUuX0Ol0hV7DmTNn7ugPmn/++YeOHTvy999/s27dOmrWrFnqfYnyJedUhUMZOHAgY8eOZe7cubz55pu3XK9SpUrs2rULVVXNgvXcuXPk5eXh5+dnWi8vL4+LFy+aBevNQVCxYkU0Gg39+vUjLi6uyGOGhoaW6jVVq1bNNDgpOjqagIAAnnrqKRISEvjggw8A8PPzIzw8/JavOSgoyPT1oEGDGDRoENeuXWPz5s0kJCTQtWtXjh49SkhICK6urgDk5OTg4uJi2u5OW7V+fn7MnDmTmTNnkpqayg8//MDo0aM5d+4cq1evvmXdSUlJZsvq1atX4mMGBwfz0EMPMX78eOrVq1eoNVmg4GeblpZW6LnTp0+b/T5A4Z9/Ucv8/PyoVKnSLV9baS/z+ueff3jooYdITk5m/fr1hIeHl2o/wjYkVIVDqVq1Ki+//DJ//PEHAwYMuOV6Dz74IN988w3Lly+nR48epuUFozgffPBBADp06MDbb7/Nl19+SXx8vGm9//3vf2b7c3d3p0OHDuzdu5fw8HB0Op01X5aZJ598kgULFjB//nxefvllQkJC6Nq1K6tWraJWrVpUrFixRPvx8PCgS5cu5Obm0r17d37//XdCQkJM3cr79+8nMjLStP6KFStuu8+CEL7dNZPVq1fn+eefZ/369Wzbtu2W6+l0OtMfFKX14osv4ubmRu/evW+5TlRUFG5ubnzxxRdm6506dYoNGzaYThPUq1ePwMBAvvrqK0aOHGn6g+yvv/5i+/btZn+8dO3ala+//hqDwUDLli3v6DUUKAjUP//8k3Xr1tGkSROr7FeUHwlV4XCmTp1623X69+/P7NmzGTBgACkpKYSFhbF161YmT55MTEwMDz30EACdOnWibdu2jBo1imvXrtG8eXO2bdvGokWLCu1z1qxZtG7dmjZt2vDss89So0YNrly5wvHjx1mxYoVplKY1vPXWW7Rs2ZKJEyeyYMECJkyYwLp162jVqhXx8fHUq1eP7OxsUlJSWLVqFXPnzqVatWoMHToUNzc3oqOjCQwM5MyZM0yZMgVvb29TgMbExODr68uQIUOYMGECzs7OLFy4sESXbFSoUIGQkBC+//57HnzwQXx9ffHz86NixYp06NCBvn37Ur9+fSpUqEBSUhKrV6+mZ8+eVntfitKpUyc6depU7Do+Pj688cYbvPbaa/Tv358+ffpw8eJFxo8fj6urq2nktpOTExMnTuTpp5+mR48eDB06lEuXLjFu3LhC3b9PPPEEX375JTExMQwfPpwWLVqg1Wo5deoUGzdupFu3bmZ/0N1OVlYWnTt3Zu/evcycOZO8vDyz8/T+/v7UqlXLgndG2ISNB0oJUawbR/8W5+bRv6qqqhcvXlSfeeYZNTAwUHV2dlZDQkLUV199Vc3OzjZb79KlS+rgwYNVHx8f1d3dXe3YsaP6xx9/FDnKNTk5WR08eLBatWpVVavVqv7+/mqrVq3USZMmma2DBaN/p02bVuTzvXv3Vp2dndXjx4+rqpo/8jY+Pl4NDQ1VtVqt6uvrqzZr1kwdM2aMevXqVVVVVfWzzz5TO3TooFapUkXV6XRqUFCQ+vjjj6v79+832/cvv/yitmrVSvXw8FCrVq2qJiQkqAsWLLjt6F9VVdWff/5ZbdKkieri4qIC6oABA9Ts7Gz1mWeeUcPDw1UvLy/Vzc1NrVevnpqQkKBeu3at2PfBUlwf/Vucm0f/FliwYIEaHh6u6nQ61dvbW+3WrZtp5PTN69WpU0fV6XRq3bp11U8++UQdMGCA2ehfVVVVvV6vvvPOO2pERITq6uqqenp6qvXr11eHDRumHjt2zLReSUb/Fvw+3OoxYMCAYrcX9kFR1RIOtRNCCCFEsWT0rxBCCGElEqpCCCGElUioCiGEEFYioSqEEEJYiYSqEEIIYSUSqkIIIYSVyOQPxTAajZw+fZoKFSoUef9JIYQQ9wZVVbly5QpBQUHF3txAQrUYp0+fJjg42NZlCCGEsBMnT56kWrVqt3xeQrUYBRNinzx5Ei8vLxtXc2f0ej1r166lU6dOaLVaW5cjhMOQz44AyMjIIDg4+LY3SpBQLUZBl6+Xl9ddEaru7u54eXnJfwxCWEA+O+JGtzsVKAOVhBBCCCuRUBVCCCGsREJVCCGEsBIJVSGEEMJKJFSFEEIIK5HRv0IIUYS0y1kkX7hGNW8Xi9YP9fMg0Nut3NYT9kVCVQhhNZYEQUnWtUawlGYfi5NSeXXpAYwqOCnweKhCTInXV5navQGPNwkEY17+w6AHo56V+/7i3TWH0KgGdIqB4R1q0LGeHxj119cxQEAYi//IMTv+lJ5hxEZWL9XrF+VLUVVVtXUR9iojIwNvb28uX758V1ynumrVKmJiYuRaOwd0p+FS2u0t2e7mICouCIpct3lwfqhcD5hlu5OZuvIAGtWAi5LHyx1rEdPQDwy5/wbVjV8XBNMNX+9JPsdP+09eD7E8OtevRKMA90Lr3biP7Jwcdhw9gwYDWgw4K3loMdCwihs6xXjDdnlg1GPM03MtKwtnDDhjQKsYLPrZ3MzgGUSdi9Mwqv9eD6lRFLaO7iAtVhsqaR5IS1XcU8qrS82ax7EkrID8YDLkXn/k8f2eZKatOogzebgoBkY+EErn+pWuB0PuDaHy7zYYcvn1z7Os3JuKM3noFANdGlYirFAg5a+flZ2Fx+9/M8fZgJa8/McKA7l7PdEpBrNt8vS5tL18lV26/LDSkofzSgPqjwYU/v0bvwfQ48ae183XHxZoCjS98X+549cfxXAFOmiKeOJ80es7ARVuOzW4gtFJS5ZBwYAGPRry0KDHGT8vD1xdXEDRwPnDaK6eRqvqyUFn2tqgqqRcyJRQdQASqkWYPXs2s2fPxmC4s784hX2dF7I4nACMxhvCRl/4a1PrKNcUHFuO/M3XO06gUfPDKLZpFSKDKxSzjXmY3fh8Tk4OISnn+FabhzMGdOThvMJA3lYtzmqeWbCZ9oN551M3oNuN4bT1+uM2mgPNb+zUOHr9UQQ3oGtRQXS68CJnILCE96cwqAp6nNHjTB4aPNxc0elcQeMMTlrQ6G74WgtOzteXabmYZWRHSkZ+gKkFQeZMx/uqEujrdX19bf72Gp1pH5dyVCavPkau6nw9+DQY0DD5PxFU9rlpOyct5zMNxH68m1zVGb2aH5ZGxZmfRj5AgI/n9fU1nL2cRfTUDRhv+PFoFIWtT19vgRoNMMEXAE8lmxxVZ7ZeDT/3kr1pwqYkVIsQFxdHXFycqbl/p+wpWMrTrbv48ooPKrOv9SVbx1hUSP37dXZ2NlWOn+Fz5zy0Sn44aVfkod/pglbNK7xdwf5Uo8Wvuw3Q5sYwOnD9UQouwP1FjdG/XPJ95KlO5KEh94Zw8va83jq6IUz+DQst6TmQlHoF/fXWlP56wDzQqCoBFb1uCjUtGXqFGeuT0V8/Tp6aHyzjejSmYgUPs2NcyDIy6PN95KrX940GFS1LXmhHleshlHbVQOt3tmBQ/33xGkVh68iSd4HmXs4ivogQ69ilAxSzDx+gmWsqry09iEFVr59TNVAxvAsUcerEHxjWo1L++qhoFIXJPe8jwN/PbL1Abzem9Awz7bdgPdPrcdKAxgUMOYx7OJQRq9OLXk/YNQnVMlaq1lFpqGoRgZLz77KcTCpeO4GSuh0wFg6ivJwiQunGfelv2m/xgZinz6H1P1fYqcv7tztwZR78mGf9114CrkD7osLpooU7UjTXg0T3bwtHozWF01WDE8cu5OSHkepsamU1DfWnkncF85ZRwbY3hFP+w8W0zqUclddXHL0hgPLP3L3frwV+Xh7X19XeYj860q7mEf1WYuHW0X+LD6ecy1k8W0QgdYgpOpC8gAYVUs0Do8d9VGxe+HfdD3iqR/VC4VIl6N91A11gcs+IWwdQCdw2xIoRG1mdtnX9SbmQSVVvHXu3bSjx+jX83G95jNuup3OHrBwebeBD88YRt92fsD8yUKkYdzpQKe16d08vp424oEdHHi5KHs+2CaaCs1pMSN0q5G4KyhvXNerL4B0oJ05acHb5N5xuCiqcry9zcr6+XlFhpCt6++tfX8pRmbT6RH4X3fXWlAEt0/s2p1IFz0JhVPTxtVDMfRTh3595oRC7g0Emi5NSCwWDJX+YlXb70myXdjmrxEFQknUt2Z81aipKuQ7ym9EIMk7B0A1QtVnZHktYRAYq2YHkC9cwqjDe+TPclNx/n9hRHkdXzAJIddKSmWvA3dMLRaO7HhQFzzvf0OLR5W/ndFMrzNnFPGCKCrEbnr+QZeTpL/aTozqbuh2NN3bxFdTm5AzlcAN4HyDStXBIVLrPur0Gd9I6upWStoKsvX1ptgv0divx/kuyriX7K8t9lBvd9fOmuZm2rUOUmoRqGQr188BJgXXGZmgwkJs/zpGYxiF4uLsX3Roq9LihFVdki66odXX552dukKfX83M5XlLjB/TpUbPYLr7ydqfhZMvj3GkwlHZ7hwqku4H2eqjqJVQdlYRqGSpotfzf0nizYPG4Ry7iLq8Qs0R5hYSEkSgVnUf+v7nXbFuHKDUJ1TJmj8FSniRchLCAtFQdnoRqOZBgEUKUiOmcqrRUHZXcpUYIIeyFzjP/XwlVhyWhKoQQ9kK6fx2ehKoQQtgLuaTG4UmoCiGEvdBeH/2rl+5fRyWhKoQQ9kJaqg5PQlUIIeyFnFN1eBKqQghhL2TyB4cnoSqEEPZCWqoOT0JVCCHshZxTdXgSqkIIYS9k9K/Dk1AVQgh7IedUHZ6EahFmz55Nw4YNiYyMtHUpQoh7iSlUpfvXUUmoFiEuLo5Dhw6RlJRk61KEEPcS00Cla6Cqtq1FlIqEqhBC2IuCgUqqEfJybFuLKBUJVSGEsBcFA5VALqtxUBKqQghhLzTOoNHlfy2DlRyShKoQQtgTmQDCoUmoCiGEPZHLahyahKoQQtgTaak6NAlVIYSwJzJVoUOTUBVCCHsiUxU6NAlVIYSwJzKrkkOTUBVCCHti6v6VlqojklAVQgh7It2/Dk1CVQgh7IkMVHJoEqpCCGFP5JIahyahKoQQ9kQmf3BoEqpCCGFPpKXq0CRUhRDCnsjoX4cmoSqEEPbENPpXWqqOSEJVCCHsiYz+dWgSqkIIYU90cp2qI5NQFUIIe6KVaQodmYSqEELYExmo5NAkVIUQwp7IQCWHJqEqhBD25MaWqqrathZhMQlVIYSwJwWTP6gGMOTathZhMQlVIYSwJwWjf0HOqzogCVUhhLAnGi04afO/lvOqDueuD9WVK1dSr1496tSpw4IFC2xdjhBC3J5MAOGwnG1dQFnKy8tj5MiRbNy4ES8vL5o2bUrPnj3x9fW1dWlCCHFrWg/IviwTQDigu7ql+ssvv9CoUSOqVq1KhQoViImJYc2aNbYuSwghiqeTCSAclV2H6ubNm3n00UcJCgpCURSWL19eaJ05c+YQGhqKq6srzZo1Y8uWLabnTp8+TdWqVU3fV6tWjb///rs8ShdCiNLTye3fHJVdh+q1a9eIiIjggw8+KPL5xYsXM2LECMaMGcPevXtp06YNXbp0ITU1FQC1iGu8FEUp05qFEOKOaeVG5Y7Krs+pdunShS5dutzy+RkzZjBkyBCefvppAGbOnMmaNWv48MMPmTJlClWrVjVrmZ46dYqWLVvecn85OTnk5OSYvs/IyABAr9ej1+vv9OXYVEH9jv46hChvtvjsaJzdcALysjJQ5TNrF0r687frUC1Obm4uu3fvZvTo0WbLO3XqxPbt2wFo0aIFBw8e5O+//8bLy4tVq1YxduzYW+5zypQpjB8/vtDytWvX4u7ubt0XYCPr1q2zdQlCOKTy/OxEpmcQBBzal0Ty397ldlxxa5mZJeuKL1GoNm3a1KKDK4rCDz/8YHY+09ouXLiAwWCgSpUqZsurVKnCmTNnAHB2dmb69Ol06NABo9HIqFGjqFSp0i33+eqrrzJy5EjT9xkZGQQHB9OpUye8vLzK5oWUE71ez7p16+jYsSNardbW5QjhMGzx2dH88CNcSqJR3Ro0iIopl2OK4hX0XN5OiUJ13759vPjii3h6et52XVVVmTp1qlk3alm6+Rypqqpmyx577DEee+yxEu3LxcUFFxeXQsu1Wu1dE0R302sRojyV62fHJf//Wo0hB418Xu1CSX/2Je7+ffnll6lcuXKJ1p0+fXpJd1tqfn5+aDQaU6u0wLlz5wq1XoUQwqHI7d8cVolG/yYnJ+Pv71/inR46dIiQkJBSF1USOp2OZs2aFTrPsW7dOlq1alWmxxZCiDIlt39zWCVqqVoakMHBwaUq5mZXr17l+PHjpu+Tk5PZt28fvr6+VK9enZEjR9KvXz+aN29OVFQU8+bNIzU1lWeeeeaOjjt79mxmz56NwWC405cghBCWk2kKHVapRv9mZ2ezf/9+zp07h9FoNHuupOcvS+LXX3+lQ4cOpu8LBhENGDCAhQsXEhsby8WLF5kwYQJpaWncd999rFq16o5byXFxccTFxZGRkYG3t4y8E0KUs4Lbv8k0hQ7H4lBdvXo1/fv358KFC4WeUxTFqq279u3bFzmBw42ee+45nnvuOasdUwghbE53fVCotFQdjsUzKj3//PP07t2btLQ0jEaj2UO6S4UQwgpkmkKHZXGonjt3jpEjR8oIWyGEKCsyTaHDsjhUe/XqRWJiYhmUIoQQApBLahyYxedUP/jgA3r37s2WLVsICwsrdEFsfHy81YqzFRn9K4SwKa10/zoqi0P1f//7H2vWrMHNzY3ExESz2YsURbkrQlVG/wohbErup+qwLA7V119/nQkTJjB69GicnOz6znFCCOGY5JIah2VxKubm5hIbGyuBKoQQZaXgnKoxD/JybVuLsIjFyThgwAAWL15cFrUIIYSAf0f/grRWHYzF3b8Gg4G3336bNWvWEB4eXmig0owZM6xWnBBC3JOcdeDknN9Szc0Et4q2rkiUkMWheuDAAZo0aQLAwYMHzZ67+TZsQgghSknnAdmXZQSwg7E4VDdu3FgWddgVuaRGCGFz2uuhKteqOhQZbVSEuLg4Dh06RFJSkq1LEULcq2SqQodUolDt2bMnGRkZJd7pk08+yblz50pdlBBC3PO0cvs3R1Si7t/vv/+e8+fPl2iHqqqyYsUKJk6cSOXKle+oOCGEuGcVTAAho38dSolCVVVV6tatW9a1CCGEKKCV+X8dUYlCtTSDk6pWrWrxNkIIIa6TSfUdUolCtV27dmVdhxBCiBsVTAAhA5Ucioz+FUIIe6STgUqOSEK1CLNnz6Zhw4ZERkbauhQhxL3KdE71qm3rEBaRUC2CXKcqhLC5CgH5/2b8bds6hEUkVIUQwh5Vqp3/78Xjtq1DWERCVQgh7JEpVE+Aqtq2FlFiFofq2bNn6devH0FBQTg7O6PRaMweQgghrMAnBBRN/ujfK2m2rkaUkMUT6g8cOJDU1FTeeOMNAgMD5c40QghRFpx1UDEE0v/Mb616Bdm6IlECFofq1q1b2bJlC40bNy6DcoQQQphUqn09VI9DaBtbVyNKwOLu3+DgYFTp3xdCiLLnWyv/Xxms5DAsDtWZM2cyevRoUlJSyqAcIYQQJpUKQvWEbesQJWZx929sbCyZmZnUqlULd3d3tFqt2fPp6elWK85W5CblQgi7IJfVOByLQ/Xdd9+96wcnxcXFERcXR0ZGBt7e3rYuRwhxryoI1X+SwZAHGov/yxblrFSjf4UQQpQDr6rg7Ap52XA5FXxr2roicRsWn1N98sknmT9/PkePHi2LeoQQQhRwcrphsJKcV3UEFoeqp6cn06dPp379+gQFBdGnTx/mzp3LH3/8URb1CSHEva3S9dapnFd1CBaH6kcffcQff/zB6dOnmTFjBt7e3syaNYtGjRoRGBhYFjUKIcS9SwYrOZRSz/1boUIFKlasSMWKFfHx8cHZ2ZmAgABr1iaEEEJC1aFYHKqvvPIK999/P35+frz++uvk5uby6quvcvbsWfbu3VsWNQohxL3rxon1hd2zePTvtGnT8Pf3JyEhgW7dutGgQYOyqEsIIQT8G6qXT4E+C7Rutq1HFMvilurevXsZM2YMv/zyC23btiUgIIDY2Fg+/PBDDh8+XBY1CiHEvcu9Erh6AyqkJ9u6GnEbFodqREQE8fHxLF26lPPnz7NmzRrc3d2Jj4/nvvvuK4sahRDi3qUoMgewAynV9Bx79+4lMTGRxMREtmzZQkZGBo0bN6ZDhw7Wrs8mZJpCIYRdqVQbTu+RUHUAFodqxYoVuXr1KhEREbRv356hQ4fStm1bvLy8yqI+m5BpCoUQduX6edWzKb9jDM8i0FvOq9ori0N10aJFd12ICiGEPdt+2YdWQOqx/cRO3cCUnmHERla3dVmiCBafU+3ataspUE+dOsXff/9t9aKEEELkS7ucxZRdegBqKGcwqvDa0oOkXc6ycWWiKBaHqtFoZMKECXh7exMSEkL16tXx8fFh4sSJGI3GsqhRCCHuWckXrvGnMX9iHX8lA2+uYlBVUi5k2rgyURSLu3/HjBnDxx9/zNSpU4mOjkZVVbZt28a4cePIzs7mzTffLIs6hRDinhTq50GW4sYJYyC1nNJo43SAn9RW1PBzt3VpoggWh+pnn33GggULeOyxx0zLIiIiqFq1Ks8995yEqhBCWFGgtxtTeoax9vtInnX6gS6aJNo89l8ZrGSnLA7V9PR06tevX2h5/fr1SU9Pt0pRQggh/hUbWZ0LXs/CVz/QxXU/To0r27okcQulmvzhgw8+KLT8gw8+ICIiwipFCSGEMOdXNwq8quKkz4Q/N9q6HHELFrdU3377bR555BF+/vlnoqKiUBSF7du3c/LkSVatWlUWNQohhFAUqN8VfvkIDq+Eel1sXZEogsUt1Xbt2nH06FF69OjBpUuXSE9Pp2fPnhw5coQ2bdqURY1CCCEAGjya/++RH8GQZ9taRJFKNU1hUFCQDEgSQojyVj0K3HwhKx3+2gY129m6InGTEoXq/v37S7zD8PDwUhcjhBCiGBpnqB8De7+AP1ZKqNqhEoVq48aNURQFVVVRFMW0XFVVALNlMgm9EEKUoQaP5Yfq4ZXw8FvgZPFZPFGGSvTTSE5O5s8//yQ5OZklS5YQGhrKnDlz2LdvH/v27WPOnDnUqlWLJUuWlHW95WL27Nk0bNiQyMhIW5cihBDmQtuBrgJcOc3+XzbIdIV2pkQt1ZCQENPXvXv35r333iMmJsa0LDw8nODgYN544w26d+9u9SLLm9ylRghht7SupFaKpnraarav/JS3v8+RCfbtiMX9BgcOHCA0NLTQ8tDQUA4dOmSVooQQQhQt7XIWU1LzJ+DppdmMVs2VCfbtiMWh2qBBAyZNmkR2drZpWU5ODpMmTaJBgwZWLU4IIYS55AvXWGtoxt9qJfyUDB7V7JAJ9u2IxZfUzJ07l0cffZTg4GDTDEq//fYbiqKwcuVKqxcohBDiX6F+HqiKhkV5HRmt/ZrBmtUsM7aTCfbthMWh2qJFC5KTk/niiy/4448/UFWV2NhY+vbti4eHR1nUKIQQ4rqCCfbfWnqN4epSGjr9xfx22TLBvp0o1eQP7u7u/Pe//7V2LUIIIUogNrI6bev6c3nlBtyOfcUDl5YC/7F1WYJSnFMNCgqib9++zJs3j6NHj5ZFTUIIIW4j0NuNgI4j8r85sgr+SbFlOeI6i0N1+vTpeHl5MWPGDOrXr09gYCBPPPEEc+fO5fDhw2VRoxBCiKJUrg81O4Bq5OqWOWw/cUFGAduYxd2/ffr0oU+fPgCcPXuWjRs3snLlSl544QWMRqPMqCSEEOXp/ufgz40Ydy9i6PYWZCluct2qDZXqnOrVq1fZunUrmzZtIjExkb179xIWFka7djIPpRBClKe0ytFkGgOp5ZTGQM0aZhu689rSg7St6y+Dl2zA4u7fli1b4u/vzxtvvEFeXh6vvfYaZ86cYc+ePbz77rtlUaMQQohbSL6Yxay8ngAMc16JN1flulUbsjhUjx07hru7OzVr1qRmzZrUrl0bHx+fMihNCCHE7YT6efCjGsVhY3W8lEyedV6BRlHkulUbsThU09PT2bhxI9HR0fz888+0a9eOgIAAYmNjmTt3blnUKIQQ4hYCvd2Y3DOC6YZYAAZqVjMjprJ0/dpIqe4ZFB4eTnx8PEuWLOGnn36iS5cuLF26lLi4OGvXJ4QQ4jZiI6sz8eWRZPg3w1XR0+3yl7Yu6Z5lcaju3buXd999l27duuHr68v999/PgQMHGD58OD/88ENZ1CiEEOI2An3c8XpkYv43ez6H9D9tW9A9yuLRv5GRkTRp0oR27doxdOhQ2rZti5eXV1nUJoQQwhI1oqH2Q3D8Z9gwCXp9YuuK7jkWh2p6erqEqBBC2KsHE+D4eji4BJoPhhqtbV3RPcXi7l8JVCGEsGOB4flhCvDjS2DQ27aee4zFoWowGHjnnXdo0aIFAQEB+Pr6mj2EEELY2AOvY3DzhfOHydj0vq2ruadYHKrjx49nxowZPP7441y+fJmRI0fSs2dPnJycGDduXBmUWP5mz55Nw4YNiYyMtHUpQghhscW/X+XVjF4AOG16ix+2JNm4onuHxaH65ZdfMn/+fF566SWcnZ3p06cPCxYsYOzYsezcubMsaix3cXFxHDp0iKQk+UUUQjiWtMtZvLr0AN8a2rLbWAdPJRvN2tdlov1yYnGonjlzhrCwMAA8PT25fPkyAF27duXHH3+0bnVCCCEsknzhGkYVVJx4Qz8Ig6rwiGYn/+xfbevS7gkWh2q1atVIS0sDoHbt2qxduxaApKQkXFxcrFudEEIIi4T6eeCk5H99SK3BZ4bOANT7ZQxkX7ZhZfcGi0O1R48erF+/HoDhw4fzxhtvUKdOHfr378/gwYOtXqAQQoiSC/TOv/WbRslP1hmGWK64B6O58jdnv3tJuoHLmMXXqU6dOtX0da9evQgODmbbtm3Url2bxx57zKrFCSGEsFxsZHXa1vUn5UImNfzcSdoBD+wcRJXj3zDwrVC69HhK7rdaRixqqer1egYNGsSff/47/VXLli0ZOXKkBKoQQtiRQG83ompVAuDpTS58mpffDTxZO5/JS3+RFmsZsShUtVoty5YtK6tahBBCWFnBwKW382JJMVYhSElnjOZzftyfJsFaBkp1TnX58uVlUIoQQghrKxi4lIUrL+uHYVQVHnfexO6fFhI9dQOLk1JtXeJdxeJzqrVr12bixIls376dZs2a4eHhYfZ8fHy81YoTQghxZwoGLr229CBJan0+NDxKnPMPTNXOZ39OTV5bepC2df3l/qtWYnGoLliwAB8fH3bv3s3u3bvNnlMURUJVCCHsTMHApR/3pzH1x160cjpEE6fjzNTN5oncN0i5kCmhaiUWh2pycnJZ1CGEEKIMBXq78Uh4IJNXHSZeH8cq3WtEOh0l3nk5Nfw62rq8u4bF51SFEEI4poKu4NMEMEafP6/AC87L0J3azvYTF2TgkhWUqKU6cuTIEu9wxowZpS5GCCFE2fr3GtaWZO4+j/uhbzB8M5jhOW9yUanIlJ5hcg3rHShRqO7du9fs+927d2MwGKhXrx4AR48eRaPR0KxZM+tXKIQQwqoCvd0I9HbjjPcUTh7cQT2nk3yge58nc1+TgUt3qETdvxs3bjQ9Hn30Udq3b8+pU6fYs2cPe/bs4eTJk3To0IFHHnmkrOsVQghhJX9ehmf0I7iiutHS6Q9edl6MQVVJuZBp69IclsXnVKdPn86UKVOoWLGiaVnFihWZNGkS06dPt2pxQgghyk6onwd/EchL+mEADHP+kS6aJGr4udu4MsdlcahmZGRw9uzZQsvPnTvHlStXrFKUEEKIslcwcOlntSUf5eX3NL7nNo/AnL9sXJnjKtWMSoMGDeK7777j1KlTnDp1iu+++44hQ4bQs2fPsqhRCCFEGYmNrM7W0R0IH/AuOdWi0OZdI2vR45w5e9rWpTkki0N17ty5PPLIIzz11FOEhIQQEhLCk08+SZcuXZgzZ05Z1CiEEKIMBXq7EVWnCj/Vn8op1Q+3K39xfHZvvtn15+03FmYsDlV3d3fmzJnDxYsX2bt3L3v27CE9PZ05c+YUmrJQCCGEY0i7nMXIH//m6dyXuKa60NrpIFdXvCbXrlqo1JM/eHh4EB4eTkREhISpEEI4uIK72fyhVmek/lkABjv/xLWdC21bmIORGZWEEEKY7mYDsMbYgpl5+WNkau58nYsH19mwMscioSqEEMI0Elij5CfrrLz/8IMhCic1D+dv+7NqY6JtC3QQEqpCCCGAf0cCf9CnCYqi8LJ+GLuNdfBWMmmUOJSzaSdtXaLdk1AVQghhEujthq+nDqMKOegYmvsiqUZ/QpRz8PWTpF1Mt3WJdk1CVQghhJkbz6+m48Ug/SgyVHeqXP6NfTNj+eYXuQXorUioCiGEMHPz+dUTalWG6f+PHNWZLppfyPxhNGmXZH7gotwTodqjRw8qVqxIr169bF2KEEI4hILzq68/0gCAHcZGvKR/BoCBzqvJ2fyeLcuzW/dEqMbHx/P555/bugwhhHAogd5uPBIeaOoKXmFsxST9kwDU2DMFDnxnw+rs0z0Rqh06dKBChQq2LkMIIRzOzV3Bnxof4UjoU/lPLhvGxd9Wsf3EBZl56boS3aS8LG3evJlp06axe/du0tLSWLZsGd27dzdbZ86cOUybNo20tDQaNWrEzJkzadOmjW0KLoLBYECv19u6jGLp9XqcnZ3Jzs7GYDDYuhxhh7RaLRqNxtZlCDsUG1mdtnX9SbmQSQ0/dwIrdIGlWXBwCW5LBzIt9zV+ow5TeoYRG1nd1uXalM1D9dq1a0RERDBo0CD+85//FHp+8eLFjBgxgjlz5hAdHc1HH31Ely5dOHToENWr5//wmjVrRk5OTqFt165dS1BQUJnVrqoqZ86c4dKlS2V2DGtRVZWAgABOnjyJcv0vTiFu5uPjQ0BAgPyOiEICvd0I9HYzfZ/2wLsc3X+Mdk77+UQ3jd65Y3ltqULbuv5m691rbB6qXbp0oUuXLrd8fsaMGQwZMoSnn34agJkzZ7JmzRo+/PBDpkyZAsDu3butUktOTo5ZOGdkZAD5rbyiWqJnz54lIyMDf39/3N3d7fo/IlVVuXbtGh4eHnZdp7ANVVXJzMzk/PnzGAwGqlSpYuuS7EbBZ9/ee6PK2/Hz2TyTO4L/6SbTxOk4i3RT6ZWTwImzGfi52zxarK6kP3+7fuW5ubns3r2b0aNHmy3v1KkT27dvt/rxpkyZwvjx4wstX7t2Le7u7mbLFEUhMDCQgIAAtFqtQ3zgdDqdQ9QpbEOr1VKhQgXS0tLYs2cPqqrauiS7sm6dzH97o0s5kI0Lg3Jf5lvdBOo4/c0Xusn88qvCqsM+ti7P6jIzS3YJkV2H6oULF4r8q7lKlSqcOXOmxPvp3Lkze/bs4dq1a1SrVo1ly5YRGRlZaL1XX32VkSNHmr7PyMggODiYTp064eXlZbZuTk4Oqamp+Pr64uZm/10dqqpy5coVKlSoIC1VcUtarZYrV67wwAMP4OLiYuty7IJer2fdunV07NgRrVZr63Lsirb6KV7//hBP5b7Kdy7jCXU6S41/PiSv6w/gVtHW5VlVQc/l7dh1qBa4OQRUVbUoGNasWVOi9VxcXIr8j0Sr1Rb6MBkMBhRFQaPR4ORk/4OojUYjkP9eOkK9wjY0Gg2KouDs7CwBcpOi/h+41/W9P5QODQJIuZCJi3Nz+LY7yvnDaBc/Af2/B5e756qLkv7s7fp/Vz8/PzQaTaFW6blz5+ScjxBC2IFAbzeialWickgD6L88v4X692743xOQe+/NumTXoarT6WjWrFmhcxnr1q2jVatWNqrq7peYmIiiKA4xqrk8DRw4sNDlXo6gRo0azJw50/S9oigsX77cZvWIu1jlBvDUUtBVgL+2wuInQZ9t66rKlc27f69evcrx48dN3ycnJ7Nv3z58fX2pXr06I0eOpF+/fjRv3pyoqCjmzZtHamoqzzzzTJnVNHv2bGbPni3Xcwozs2bNksE7QtxO1abw5LcYF/XA6cQGsv/3FK5P/g+cdbaurFzYvKX666+/0qRJE5o0aQLAyJEjadKkCWPHjgUgNjaWmTNnMmHCBBo3bszmzZtZtWoVISEhZVZTXFwchw4dIikpqcyOIRyPt7c3Pj4+ti5DCLu3+FxV+mWOJFvV4pq8jpML+oIhz9ZllQubh2r79u1RVbXQY+HChaZ1nnvuOVJSUsjJyWH37t20bdvWdgWXkbTLWeU21VdOTg7x8fFUrlwZV1dXWrduXeQfENu2bSMiIgJXV1datmzJgQMHTM/99ddfPProo1SsWBEPDw8aNWrEqlWr7qiugu7V8ePHU7lyZby8vBg2bBi5ubmmdW7uygRo3Lgx48aNM32vKAofffQRXbt2xd3dnQYNGrBjxw6OHz9O+/bt8fDwICoqihMnTpi2GTduHI0bN+ajjz4iODgYd3d3evfubdYFfnP3b/v27YmPj2fUqFH4+voSEBBgVgfAH3/8QevWrXF1daVhw4b8/PPPxXa/rlixAh8fH9PAsn379uXfLPrll03rDBs2jD59+pi+3759O23btsXNzY3g4GDi4+O5du3abd5tIcpG2uUsXl16gG3GRgzTjyRHdSb4zDqyvhlyTwSrzUNVwOKkVKKnbqDv/F1ET93A4qTUMj3eqFGjWLJkCZ999hl79uyhdu3adO7cmfR085sPv/zyy7zzzjskJSVRuXJlHnvsMdN1rnFxceTk5LB582YOHDjAW2+9haen5x3Xtn79eg4fPszGjRv56quvWLZsWZHXDt/OxIkT6d+/P/v27aN+/fr07duXYcOG8eqrr/Lrr78C8Pzzz5ttc/z4cb755htWrFjB6tWr2bdvH3FxccUe57PPPsPDw4Ndu3bx9ttvM2HCBNMYAKPRSPfu3XF3d2fXrl3MmzePMWPGFLu/tm3bcuXKFfbu3QvApk2b8PPzY9OmTaZ1EhMTadeuHQAHDhygc+fO9OzZk/3797N48WK2bt1a6LUJUV6SL1zDeP0sySZjBHH64ehVDW5HlpP57TC2Hzt7V88TLKFqYwV/1RX8EhpVeG3pwTL7pbt27Roffvgh06ZNo0uXLjRs2JD58+fj5ubGxx9/bLZuQkICHTt2JCwsjM8++4yzZ8+ybNkyAFJTU4mOjiYsLIyaNWvStWtXq/Qg6HQ6PvnkExo1asQjjzzChAkTeO+990wtt5IaNGgQjz/+OHXr1uWVV14hJSWFJ598ks6dO9OgQQOGDx9OYmKi2TbZ2dl89tlnNG7cmLZt2/L+++/z9ddfF3tNdHh4OAkJCdSpU4f+/fvTvHlz1q9fD+RPGnLixAk+//xzIiIiaN26NW+++WaxdXt7e9O4cWNTbYmJifzf//0fv/32G1euXOHMmTMcPXqU9u3bAzBt2jT69u3LiBEjqFOnDq1ateK9997j888/Jzv73hogIuzDjTc4B/jZ2Iz4vOEYFQ3uf3zH6c+fpvXUn8u88WArEqo2duNfdQUMqkrKhbIZin7ixAn0ej3R0dGmZVqtlhYtWnD48GGzdaOiokxf+/r6Uq9ePdM68fHxTJo0iejoaBISEti/f/8tjzl58mQ8PT1Nj9TUW3+YIiIizGavioqK4urVq5w8edKi1xkeHm76uuDyq7CwMLNl2dnZZhd0V69enWrVqpkd22g0cuTIkRIdByAwMJBz584BcOTIEYKDgwkICDA936JFi9vW3r59exITE1FVlS1bttCtWzfuu+8+tm7dysaNG6lSpQr169cH8qfoXLhwodn727lzZ4xGI8nJybc9lhDWdvNdbTSKQuNOT/FCbhx5qhO9NJt5U7OAMUv335UtVpuP/rVH5Tn6t+CvuhuDVaMo1PBzv/VGd6Bg9GppJ9QoWOfpp5+mc+fO/Pjjj6xdu5YpU6Ywffp0XnjhhULbPPPMMzz++OOm70tzk4OC4zo5ORUagVvU1Is3XqhdsG1Ry4prAResU9z7cvMF4YqimPZp6SQlBdq3b8/HH3/Mb7/9hpOTEw0bNqRdu3Zs2rSJf/75x9T1W1D/sGHDiI+PL7SfghtOCFHebr6rTfKFa0wx3I+iqszSfsATzokApJxvCeQ3LkL9PO6KifilpVqE8hz9W9RfdZN73ldmv1y1a9dGp9OxdetW0zK9Xs+vv/5KgwYNzNbduXOn6et//vmHo0ePmlpIAMHBwTzzzDMsXbqUF198kfnz5xd5TF9fX2rXrm16ODvf+m+53377jaysf/963blzJ56enqYWpL+/P2lpaabnMzIyrNYiS01N5fTp06bvd+zYgZOTE3Xr1i3V/urXr09qaipnz541LSvJ71TBedWZM2fSrl07FEWhXbt2JCYmmp1PBWjatCm///672ftb8NDp7o1LGIR9KpgUItDbzdR4WGmM4v/0cRhUhSecE6m4/kVaT/253MaTlAcJVTsQG1mdraM78NXQ+9k6ukOZ3o/Qw8ODZ599lpdffpnVq1dz6NAhhg4dSmZmJkOGDDFbd8KECaxfv56DBw8ycOBA/Pz8TKNfR4wYwZo1a0hOTmbPnj1s2LChUCiXRm5uLkOGDOHQoUP89NNPJCQk8Pzzz5umVnzggQdYtGgRW7Zs4eDBgwwYMMBq9wB1dXVlwIAB/Pbbb2zZsoX4+Hgef/xxs+5bS3Ts2JFatWoxYMAA9u/fz7Zt20wDlYprwRacV/3iiy9M507btm3Lnj17zM6nArzyyivs2LGDuLg49u3bx7Fjx/jhhx+K7DEQwlZubDz8YGzFi3lxGHGiftr3TNHMR8FY5uNJyot0/9qJm+9VWJamTp2K0WikX79+XLlyhebNm7NmzRoqVqxYaL3hw4dz7NgxIiIi+OGHH0ytH4PBQFxcHKdOncLLy4uHH36Yd999945re/DBB6lTpw5t27YlJyeHJ554wuwylVdffZU///yTrl274u3tzcSJE63WUq1duzY9e/YkJiaG9PR0YmJimDNnTqn3p9FoWL58OU8//TSRkZHUrFmTadOm8eijj+Lq6lrsth06dGDPnj2mAK1YsSINGzbk9OnTZn+8hIeHs2nTJsaMGUObNm1QVZVatWoRGxtb6rqFKAvmXcIPcCypNrW3/B+PO2/CSVEZpf8vBtWJlAuZDt0NrKgyRcwtZWRk4O3tzeXLlwvdpSY7O5vk5GRCQ0Nv+x+kPTAajWRkZODl5WW3E+oPHDiQS5cu2WQKvXHjxrF8+XL27dtXpsfZtm0brVu35vjx49SqVatMj1UajvZ7XR70ej2rVq0iJiZGJtS3orTLWUx++03edZ6Ns2JkqaE1r+Q9y+bRD9plqBaXBzeSlqoQZWjZsmV4enpSp04djh8/zvDhw4mOjrbLQBWiPAV6u9G6+38ZvlzDTOf36anZSrPgCmBsx/YTFxx24JKEqhBl6MqVK4waNYqTJ0/i5+fHQw89xPTp021dlhB2Ib9LeDR/7qlH3S3xhJz+iR9n/Ifh+jiMijNTeoYRG1mdtMtZDjNCWLp/i3DjJTVHjx6V7l9xz3C03+vyIN2/5SN9z/d4fD8YFyWPtYZmPK+Px6DoGNWlHm/99AdGFZwUTEFb3kra/Sv/uxZBJtQXQojy9Yd39PW5grV00uzmI+0MnNUcpl4PVCj7GeesQUJVCCGEzYX6ebBZbcwg/ctkqTo6aH7jY+07uKg5ZuuV5Yxz1iChKoQQwuYKrmXdpYYxIPcVrqqutNYc5HPdW3jwb8tUoyi465zK7a5elpKBSkIIIezCv9eytiQ7tyWey/vSIucPvtRNpn/uK1xTKtC9SRA95my3+TnWW5GWqhBCCLtRML2hX4M2MOAHcKtIY6cT7AicyfeD6rJs7992fY5VQlUIIYR9CmoCA38ED388/jlEzR+fwE/9x2wVezvHKqFahNmzZ9OwYUMiIyNtXYrNtW/fnhEjRti6DKtISUlBUZQynzVJCGFFVRrBoJ+gQhDul4/xjcsEqnLe9HRZ3tWrNCRUi3C3X1JTVFAmJiaiKAqXLl0yW7506VImTpxYfsWVo1u9ZiGEnfGrA4NWgU8INZSzfOsygVAlzXRXL8BuBi7JQCVRLF9f33I5jl6vlwvrhRC35hua32L9vBtBF4+xxnsKl3t/y4Z0iJ66wW4GLklL9R4zaNAgNm3axKxZs1AUBUVRSElJoUOHDkD+3VAURWHgwIFA4VZtjRo1mDRpEv3798fT05OQkBC+//57zp8/T7du3fD09CQsLIxff/212DoURWHu3Ll069YNDw8PJk2aBMCKFSto1qwZrq6u1KxZk/Hjx5OXl2fabty4cVSvXh0XFxeCgoLMbs6tKEqhyfh9fHxYuHBhoeMX95q/++47wsLCcHNzo1KlSjz00ENcu3atJG+vEKIseVfND9YqYeiyL1Dp2558tWy5XQ1cklC1JlWF3Gu2eZRwtsmZM2cSFRXF0KFDSUtLIy0tjeDgYJYsWQLAkSNHSEtLY9asWbfcx7vvvkt0dDR79+7lkUceoV+/fvTv35+nnnqKPXv2ULt2bfr378/tZsBMSEigW7duHDhwgMGDB7NmzRqeeuop4uPjOXToEB999BELFy7kzTffBPLD7t133+Wjjz7i2LFjLF++nLCwsBL+cMzd6jWnpaXRp08fBg8ezOHDh0lMTKRnz563fS1CiHLi6Q8DV0C1SJyy/+EL7Zu0VA6bnrb1wCXp/rUmfSZMDrLNsV87DTqP267m7e2NTqfD3d3d7ObbBd28lStXxsfHp9h9xMTEMGzYMADGjh3Lhx9+SGRkJL179wbyb5wdFRXF2bNni73Bd9++fRk8eLDp+379+jF69GgGDBgAQM2aNZk4cSKjRo0iISGB1NRUAgICeOihh9BqtVSvXp0WLVrc9jUXRaPRFPmaT5w4QV5eHj179iQkJASg1MEthCgjbhWh33JyvojF8+RWPtNN5Rn9CBKNTWw+cElaqsJi4eHhpq+rVKkCmAdPwbJz584Vu5/mzZubfb97924mTJiAp6en6VHQos7MzKR3795kZWVRs2ZNhg4dyrJly8y6hq0hIiKCBx98kLCwMHr37s38+fP5559/br+hEKJ8uXji0n8Jf1duj6uiZ752Bo9qdtp84JK0VK1J657fYrTVscvrUDcMKFIU5ZbLjEZjsfvx8DBvWRuNRsaPH0/Pnj0Lrevq6kpwcDBHjhxh3bp1/Pzzzzz33HNMmzaNTZs2odVqURSlUDetXq+36LVpNBrWrVvH9u3bWbt2Le+//z5jxoxh165dhIaGWrQvIUQZ07pSddh3ZH0zFLcjy3hP9wFJqZWIXlrHZgOXpKVqTYqS3wVri8f1ICsJnU6HwWAotAwotLw8NW3alCNHjlC7du1Cj4Lb1bm5ufHYY4/x3nvvkZiYyI4dOzhw4AAA/v7+pKWlmfZ37NgxMjNvfW7lVq9ZURSio6MZP348e/fuRafTsWzZMmu/XCGENWi0uMV+DE0HoKhGWhxIYKDTT4BtBi5JS7UIN95P9W5Uo0YNdu3aRUpKCp6envj6+hISEoKiKKxcuZKYmBjc3Nzw9PQs17rGjh1L165dCQ4Opnfv3jg5ObF//34OHDjApEmTWLhwIQaDgZYtW+Lu7s6iRYtwc3Mznft84IEH+OCDD7j//vsxGo288sorxV6mU9Rr/v3331m/fj2dOnWicuXK7Nq1i/Pnz9OgQYPyehuEEJZy0sCjs/g7W0vVQwsYq11EBTKZZeiJQYWUC5nldnNzaakW4W6f/OGll15Co9HQsGFD/P39SU1NpWrVqowfP57Ro0dTpUoVnn/++XKvq3PnzqxcuZJ169YRGRnJ/fffz4wZM0yh6ePjw/z584mOjiY8PJz169ezYsUKKlWqBMD06dMJDg6mbdu29O3bl5deegl391t3ixf1mr28vNi8eTMxMTHUrVuX119/nenTp9OlS5dyeQ+EEKWkKDh1msj0vPwBk/+nXcIbzl/grKjlOnBJUeVagVsq7k7v2dnZJCcnExoaiqurq40qLDmj0UhGRgZeXl6mrlQhbuZov9flQa/Xs2rVKmJiYmSCEgewOCmVP5a/Q4L2MwD+rNadmoM+Bs2ddcwWlwc3kv9dhRBC3DViI6vz31FvcazVO6iKEzVPLYfvBkFezm23tQYJVSGEEHeVQG836nQaivL456gaHRz+gWsrXy2XY0uoCiGEuCstvhpB/6wXOWQMocOuZixOSi3zY0qoCiGEuOukXc7i1aUH2GIM45HcNzmn+pTL5TUSqndIxnmJu4n8Pou7RfKFa6aJ9tXrUVce8wJLqJZSwSjA4iYXEMLRFPw+yyhX4ehC/TxwumlOnPKYF1gmfygljUaDj4+PaX5bd3d30/R89shoNJKbm0t2drZcUiMKUVWVzMxMzp07h4+PDxqNxtYlCXFHAr3dmNIzjNeWHsSgqqYbmpf1JBASqneg4A4st5s43h6oqkpWVhZubm52Hf7Ctnx8fIq9s5AQjiQ2sjpt6/qTciGTGn7u5TKrkoRqEUo6TaGiKAQGBlK5cmWLJ24vb3q9ns2bN9O2bVvp2hNF0mq10kIVd51Ab7dym6IQJFSLFBcXR1xcnGkGjdvRaDR2/5+RRqMhLy8PV1dXCVUhhCgjcnJNCCGEsBIJVSGEEMJKJFSFEEIIK5FzqsUouBA+IyPDxpXcOb1eT2ZmJhkZGXJOVQgLyGdHwL85cLsJUiRUi3HlyhUAgoODbVyJEEIIe3DlypViB7DK/VSLYTQaOX36NBUqVDC7tjMyMrJUNzC3dLuSrl+S9TIyMggODubkyZPF3gvwblTan1dZKq+arH2cO92fvX12SrLuvfzZAfv7/Njqs6OqKleuXCEoKKjYCXSkpVoMJycnqlWrVmi5RqMp1YfL0u1Kur4l+/Xy8rrn/mMo7c+rLJVXTdY+zp3uz94+O5asey9+dsD+Pj+2/OyU5BJLGahUCnFxceWyXUnXL2099wp7fH/KqyZrH+dO92dvn53S7PteY2/vj71/dqT79x5RMJHF5cuX7eqvTiHsnXx2hCWkpXqPcHFxISEhARcXF1uXIoRDkc+OsIS0VIUQQggrkZaqEEIIYSUSqkIIIYSVSKgKIYQQViKhKoQQQliJhKoQQghhJRKqwsyVK1eIjIykcePGhIWFMX/+fFuXJITDOHnyJO3bt6dhw4aEh4fz7bff2rokUc7kkhphxmAwkJOTg7u7O5mZmdx3330kJSVRqVIlW5cmhN1LS0vj7NmzNG7cmHPnztG0aVOOHDmCh4eHrUsT5UTm/hVmNBoN7u7uAGRnZ2MwGG57qyMhRL7AwEACAwMBqFy5Mr6+vqSnp0uo3kOk+9fBbN68mUcffZSgoCAURWH58uWF1pkzZw6hoaG4urrSrFkztmzZYtExLl26REREBNWqVWPUqFH4+flZqXohbKs8Pj8Ffv31V4xGo9w68h4joepgrl27RkREBB988EGRzy9evJgRI0YwZswY9u7dS5s2bejSpQupqammdZo1a8Z9991X6HH69GkAfHx8+O2330hOTuZ///sfZ8+eLZfXJkRZK4/PD8DFixfp378/8+bNK/PXJOyLnFN1YIqisGzZMrp3725a1rJlS5o2bcqHH35oWtagQQO6d+/OlClTLD7Gs88+ywMPPEDv3r2tUbIQdqOsPj85OTl07NiRoUOH0q9fP2uXLeyctFTvIrm5uezevZtOnTqZLe/UqRPbt28v0T7Onj1LRkYGkH93js2bN1OvXj2r1yqEvbHG50dVVQYOHMgDDzwggXqPkoFKd5ELFy5gMBioUqWK2fIqVapw5syZEu3j1KlTDBkyBFVVUVWV559/nvDw8LIoVwi7Yo3Pz7Zt21i8eDHh4eGm87WLFi0iLCzM2uUKOyWhehdSFMXse1VVCy27lWbNmrFv374yqEoIx3Ann5/WrVtjNBrLoizhIKT79y7i5+eHRqMp9Ff1uXPnCv31LYQwJ58fYQ0SqncRnU5Hs2bNWLdundnydevW0apVKxtVJYRjkM+PsAbp/nUwV69e5fjx46bvk5OT2bdvH76+vlSvXp2RI0fSr18/mjdvTlRUFPPmzSM1NZVnnnnGhlULYR/k8yPKnCocysaNG1Wg0GPAgAGmdWbPnq2GhISoOp1Obdq0qbpp0ybbFSyEHZHPjyhrcp2qEEIIYSVyTlUIIYSwEglVIYQQwkokVIUQQggrkVAVQgghrERCVQghhLASCVUhhBDCSiRUhRBCCCuRUBVCCCGsREJViLtIYmIiiqJw6dKlcj+2oigoioKPj0+x640bN47GjRubfV+w7cyZM8u0RiHKmoSqEA6qffv2jBgxwmxZq1atSEtLw9vb2yY1ffrppxw9etSibV566SXS0tKoVq1aGVUlRPmRCfWFuIvodDoCAgJsdnwfHx8qV65s0Taenp54enqi0WjKqCohyo+0VIVwQAMHDmTTpk3MmjXL1HWakpJSqPt34cKF+Pj4sHLlSurVq4e7uzu9evXi2rVrfPbZZ9SoUYOKFSvywgsvYDAYTPvPzc1l1KhRVK1aFQ8PD1q2bEliYmKpap06dSpVqlShQoUKDBkyhOzsbCu8A0LYJwlVIRzQrFmziIqKYujQoaSlpZGWlkZwcHCR62ZmZvLee+/x9ddfs3r1ahITE+nZsyerVq1i1apVLFq0iHnz5vHdd9+Zthk0aBDbtm3j66+/Zv/+/fTu3ZuHH36YY8eOWVTnN998Q0JCAm+++Sa//vorgYGBzJkz545euxD2TLp/hXBA3t7e6HQ63N3db9vdq9fr+fDDD6lVqxYAvXr1YtGiRZw9exZPT08aNmxIhw4d2LhxI7GxsZw4cYKvvvqKU6dOERQUBOSf91y9ejWffvopkydPLnGdM2fOZPDgwTz99NMATJo0iZ9//llaq+KuJS1VIe5y7u7upkAFqFKlCjVq1MDT09Ns2blz5wDYs2cPqqpSt25d0/lOT09PNm3axIkTJyw69uHDh4mKijJbdvP3QtxNpKUqxF1Oq9Wafa8oSpHLjEYjAEajEY1Gw+7duwsNHroxiIUQhUmoCuGgdDqd2eAia2nSpAkGg4Fz587Rpk2bO9pXgwYN2LlzJ/379zct27lz552WKITdklAVwkHVqFGDXbt2kZKSgqenJ76+vlbZb926dXnyySfp378/06dPp0mTJly4cIENGzYQFhZGTExMifc1fPhwBgwYQPPmzWndujVffvklv//+OzVr1rRKrULYGzmnKoSDeumll9BoNDRs2BB/f39SU1Ottu9PP/2U/v378+KLL1KvXj0ee+wxdu3adcsRxrcSGxvL2LFjeeWVV2jWrBl//fUXzz77rNXqFMLeKKqqqrYuQgjh+BRFYdmyZXTv3r1U29eoUYMRI0YUmiVKCEciLVUhhNX06dPH4ukGJ0+ejKenp1Vb2kLYirRUhRBWcfz4cQA0Gg2hoaEl3i49PZ309HQA/P39bTZvsRDWIKEqhBBCWIl0/wohhBBWIqEqhBBCWImEqhBCCGElEqpCCCGElUioCiGEEFYioSqEEEJYiYSqEEIIYSUSqkIIIYSVSKgKIYQQVvL/ERqQZ448xsAAAAAASUVORK5CYII=",
      "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",
    "hm1 = w1.headinside(tm)\n",
    "plt.loglog(to, -ho, \".\", label=\"obs - pumping well\")\n",
    "plt.loglog(tm, -hm1[0], label=\"ttim results\")\n",
    "plt.xlabel(\"time [d]\")\n",
    "plt.ylabel(\"drawdown [m]\")\n",
    "plt.legend()\n",
    "plt.title(\"Model Results - Model 2\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit is still very good. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis and summary of values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "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>k [m/d]</th>\n",
       "      <th>Ss [1/m]</th>\n",
       "      <th>res</th>\n",
       "      <th>RMSE [m]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>MLU</th>\n",
       "      <td>51.53</td>\n",
       "      <td>0.000816</td>\n",
       "      <td>0.022</td>\n",
       "      <td>0.007560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTim-single layer</th>\n",
       "      <td>48.653125</td>\n",
       "      <td>0.000076</td>\n",
       "      <td>0.022092</td>\n",
       "      <td>0.008737</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TTim-two layers</th>\n",
       "      <td>44.529551</td>\n",
       "      <td>0.000006</td>\n",
       "      <td>0.016205</td>\n",
       "      <td>0.005512</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     k [m/d]  Ss [1/m]       res  RMSE [m]\n",
       "MLU                    51.53  0.000816     0.022  0.007560\n",
       "TTim-single layer  48.653125  0.000076  0.022092  0.008737\n",
       "TTim-two layers    44.529551  0.000006  0.016205  0.005512"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ta = pd.DataFrame(\n",
    "    columns=[\"k [m/d]\", \"Ss [1/m]\", \"res\"],\n",
    "    index=[\"MLU\", \"TTim-single layer\", \"TTim-two layers\"],\n",
    ")\n",
    "ta.loc[\"TTim-single layer\"] = ca0.parameters[\"optimal\"].values\n",
    "ta.loc[\"TTim-two layers\"] = ca1.parameters[\"optimal\"].values\n",
    "ta.loc[\"MLU\"] = [51.530, 8.16e-04, 0.022]\n",
    "ta[\"RMSE [m]\"] = [0.00756, ca0.rmse(), ca1.rmse()]\n",
    "ta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both TTim models agree with each other with similar parameters. MLU parameters are higher than the ones in TTim."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "* Carlson F, Randall J (2012) MLU: a Windows application for the analysis of aquifer tests and the design of well fields in layered systems. Ground Water 50(4):504–510\n",
    "* Duffield, G.M., 2007. AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
    "* Newville, M.,Stensitzki, T., Allen, D.B., Ingargiola, A. (2014) LMFIT: Non Linear Least-Squares Minimization and Curve Fitting for Python.https://dx.doi.org/10.5281/zenodo.11813. https://lmfit.github.io/lmfit-py/intro.html (last access: August,2021).\n",
    "* Yang, Xinzhu (2020) Application and comparison of different methodsfor aquifer test analysis using TTim. Master Thesis, Delft University of Technology (TUDelft), Delft, The Netherlands."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}