{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# 4. Slug test for confined aquifer - Dawsonville"
]
},
{
"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 timflow.transient as tft\n",
"\n",
"plt.rcParams[\"figure.figsize\"] = [5, 3]"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Introduction and Conceptual Model\n",
"\n",
"This Slug Test, was reported by Cooper Jr et al. (1967), and it was performed in Dawsonville, Georgia, USA. \n",
"\n",
"A fully penetrated well (Ln-2) is screened in a confined aquifer, located between depths 24 and 122 (98 m thick). The volume of the slug is 10.16 litres. Head change has been recorded at the slug well. Both the well and the casing radii of the slug well is 0.076 m."
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load data"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"data = np.loadtxt(\"data/dawsonville_slug.txt\")\n",
"to = data[:, 0]\n",
"ho = data[:, 1]"
]
},
{
"cell_type": "markdown",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"### Parameters and model"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# known parameters\n",
"b = 98 # aquifer thickness in m\n",
"zt = -24 # top of aquifer in m\n",
"zb = zt - b # bottom of aquifer in m\n",
"rw = 0.076 # well radius of Ln-2 Well in m\n",
"rc = 0.076 # casing radius of Ln-2 Well in m\n",
"Q = 10.16 / 1000 # slug volume in m^3 (10.16 l = 0.01016 m^3)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"self.neq 1\n",
"solution complete\n"
]
}
],
"source": [
"ml = tft.ModelMaq(kaq=10, z=[zt, zb], Saq=1e-4, tmin=1e-6, tmax=1e-3, topboundary=\"conf\")\n",
"w = tft.Well(ml, xw=0, yw=0, rw=rw, rc=rc, tsandQ=[(0, -Q)], layers=0, wbstype=\"slug\")\n",
"ml.solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Estimate aquifer parameters\n",
"\n",
"We calibrate hydraulic conductivity and specific storage, as in the KGS model (Hyder et al. 1994)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"...............................\n",
"Fit succeeded.\n"
]
}
],
"source": [
"# unknown parameters: kay, Saq\n",
"cal = tft.Calibrate(ml)\n",
"cal.set_parameter(name=\"kaq0\", initial=10, pmin=0, layers=0)\n",
"cal.set_parameter(name=\"Saq0\", initial=1e-4, layers=0)\n",
"cal.seriesinwell(name=\"obs\", element=w, t=to, h=ho)\n",
"cal.fit(report=False)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" layers | \n",
" optimal | \n",
" pmin | \n",
" pmax | \n",
" initial | \n",
" inhoms | \n",
" parray | \n",
"
\n",
" \n",
" \n",
" \n",
" | kaq0_0_0 | \n",
" 0 | \n",
" 0.420908 | \n",
" 0.0 | \n",
" inf | \n",
" 10.0000 | \n",
" None | \n",
" [[0.42090788883700814]] | \n",
"
\n",
" \n",
" | Saq0_0_0 | \n",
" 0 | \n",
" 0.000017 | \n",
" -inf | \n",
" inf | \n",
" 0.0001 | \n",
" None | \n",
" [[1.7004312750679313e-05]] | \n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" layers optimal pmin pmax initial inhoms \\\n",
"kaq0_0_0 0 0.420908 0.0 inf 10.0000 None \n",
"Saq0_0_0 0 0.000017 -inf inf 0.0001 None \n",
"\n",
" parray \n",
"kaq0_0_0 [[0.42090788883700814]] \n",
"Saq0_0_0 [[1.7004312750679313e-05]] "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"rmse: 0.004409616694582213\n"
]
}
],
"source": [
"display(cal.parameters)\n",
"print(\"rmse:\", cal.rmse())"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAdsAAAFBCAYAAAA/hwURAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAV7dJREFUeJzt3XdcVfX/wPHX4XLZMpQhKoKSCxcG7hQtR85wJK4cmZUjK9K+auWoFDNXjqIsR6lpmuOXmVtyjxypuXMrojgAQeHCPb8/rtxEQLnK5TLez8fjPC733M85533vh8ubc85nKKqqqgghhBDCbKwsHYAQQghR2EmyFUIIIcxMkq0QQghhZpJshRBCCDOTZCuEEEKYmSRbIYQQwswk2QohhBBmJslWCCGEMDNJtkIIIYSZSbIVT01RFMaMGWPydufPn0dRFObNm5frMeU1Pz8/+vTpY+kwTJZXdTBmzBgURcnVfTZp0oQmTZrk6j6Limepjz59+uDn55e7ARUhkmwLuHnz5qEoCoqisH379kyvq6qKj48PiqLQtm1bC0T49KKioozvTVEUNBoNnp6edO7cmePHj1s6vCwdO3aMMWPGcP78eYvF8NtvvxESEoKnpycODg6UL1+eLl26sHbtWovFJERRZ23pAETusLOzY9GiRbzwwgsZ1v/5559cvnwZW1tbC0X27IYMGULt2rXR6XQcPnyYyMhIoqKiOHr0KCVLlrR0eBkcO3aMsWPH0qRJE4ucBUyaNIlhw4YREhLCiBEjcHBw4MyZM2zcuJHFixfz8ssvA+Dr68u9e/fQarV5HqMQRZEk20KidevWLF26lOnTp2Nt/V+1Llq0iKCgIGJjYy0Y3bNp1KgRnTt3Nj6vVKkSAwYM4Mcff+TDDz+0YGT5S2pqKp999hnNmzdn/fr1mV6/fv268WdFUbCzs8vL8IoEvV5PSkqKfLYiE7mMXEh069aNmzdvsmHDBuO6lJQUli1bRvfu3bPcJjExkQ8++AAfHx9sbW2pVKkSkyZN4tGJoJKTk3n//ffx8PCgWLFitG/fnsuXL2e5zytXrvD666/j5eWFra0tVatWZc6cObn3RjEkX4B///33qY49Y8YMqlatioODA25ubgQHB7No0SLj69ndm3rS/a558+bx6quvAtC0aVPj5e+oqCgA/vrrL1q2bIm7uzv29vaUK1eO119/3dS3n63Y2Fji4+Np2LBhlq97enoaf87qnm2fPn1wcnLiypUrhIaG4uTkhIeHB0OHDiUtLS3Dvm7evMlrr72Gs7Mzrq6u9O7dm7///jvH94EXLFhAUFAQ9vb2FC9enK5du3Lp0qWnet8pKSmMGjWKoKAgXFxccHR0pFGjRmzZssVYRlVV/Pz8eOWVVzJtf//+fVxcXHjrrbeM65KTkxk9ejTPPfcctra2+Pj48OGHH5KcnJxhW0VRGDx4MAsXLqRq1arY2to+9nK9n58fbdu2JSoqiuDgYOzt7alevbrxd2T58uVUr14dOzs7goKCOHjwYKZ9bN68mUaNGuHo6IirqyuvvPJKlrdVtm/fTu3atbGzs8Pf359vv/0227hysz5E1uTMtpDw8/Ojfv36/Pzzz7Rq1QqAP/74g7i4OLp27cr06dMzlFdVlfbt27Nlyxb69etHYGAg69atY9iwYVy5coWpU6cay77xxhssWLCA7t2706BBAzZv3kybNm0yxRATE0O9evWMf4A8PDz4448/6NevH/Hx8bz33nu58l7T74e6ubmZfOzZs2czZMgQOnfuzLvvvsv9+/c5fPgwe/bsyfafkpxq3LgxQ4YMYfr06YwcOZIqVaoAUKVKFa5fv06LFi3w8PBg+PDhuLq6cv78eZYvX/5Mx3yYp6cn9vb2/Pbbb7zzzjsUL17c5H2kpaXRsmVL6taty6RJk9i4cSOTJ0/G39+fAQMGAIazt3bt2rF3714GDBhA5cqVWbVqFb17987RMcaNG8cnn3xCly5deOONN7hx4wYzZsygcePGHDx4EFdXV5Nijo+P5/vvv6dbt27079+fhIQEfvjhB1q2bMnevXsJDAxEURR69uzJxIkTuXXrVobP5rfffiM+Pp6ePXsa31/79u3Zvn07b775JlWqVOHIkSNMnTqVU6dOsXLlygzH37x5M7/88guDBw/G3d39ibcPzpw5Q/fu3Xnrrbfo2bMnkyZNol27dkRGRjJy5EgGDhwIQEREBF26dOHkyZNYWRnOizZu3EirVq0oX748Y8aM4d69e8yYMYOGDRty4MAB47GPHDli/H0bM2YMqampjB49Gi8vL7PXh8iGKgq0uXPnqoC6b98+debMmWqxYsXUpKQkVVVV9dVXX1WbNm2qqqqq+vr6qm3atDFut3LlShVQP//88wz769y5s6ooinrmzBlVVVX10KFDKqAOHDgwQ7nu3burgDp69Gjjun79+qne3t5qbGxshrJdu3ZVXVxcjHGdO3dOBdS5c+c+9r1t2bJFBdQ5c+aoN27cUK9evaquXbtWfe6551RFUdS9e/eafOxXXnlFrVq16mOP27t3b9XX1zfT+tGjR6uPfmV8fX3V3r17G58vXbpUBdQtW7ZkKLdixQpjPZnTqFGjVEB1dHRUW7VqpY4bN07dv39/pnJZ1UHv3r1VQP30008zlK1Vq5YaFBRkfP7rr7+qgDpt2jTjurS0NPXFF1/MtM9HP7Pz58+rGo1GHTduXIZjHDlyRLW2ts60PishISFqSEiI8XlqaqqanJycoczt27dVLy8v9fXXXzeuO3nypAqo33zzTYay7du3V/38/FS9Xq+qqqr+9NNPqpWVlbpt27YM5SIjI1VA3bFjh3EdoFpZWan//PPPE+NWVcPvC6Du3LnTuG7dunUqoNrb26sXLlwwrv/2228z/S4FBgaqnp6e6s2bN43r/v77b9XKykrt1auXcV1oaKhqZ2eXYX/Hjh1TNRrNU9dHdt8LkTNyGbkQ6dKlC/fu3WP16tUkJCSwevXqbM/W1qxZg0ajYciQIRnWf/DBB6iqyh9//GEsB2Qq9+hZqqqq/Prrr7Rr1w5VVYmNjTUuLVu2JC4ujgMHDjzV+3r99dfx8PCgVKlSvPzyy8TFxfHTTz9Ru3Ztk4/t6urK5cuX2bdv31PF8rTSzw5Wr16NTqcz23HGjh3LokWLqFWrFuvWreOjjz4iKCiI559/PsctuN9+++0Mzxs1asTZs2eNz9euXYtWq6V///7GdVZWVgwaNOiJ+16+fDl6vZ4uXbpkqKeSJUtSoUKFDJd+c0qj0WBjYwMYzkpv3bpFamoqwcHBGX7nKlasSN26dVm4cKFx3a1bt/jjjz/o0aOH8RbB0qVLqVKlCpUrV84Q44svvgiQKcaQkBACAgJyHG9AQAD169c3Pq9bty4AL774ImXLls20Pv2zj46O5tChQ/Tp0yfDmXmNGjVo3ry58bualpbGunXrCA0NzbC/KlWq0LJlywyxmKM+RNbkMnIh4uHhQbNmzVi0aBFJSUmkpaVlaFj0sAsXLlCqVCmKFSuWYX36pc8LFy4YH62srPD3989QrlKlShme37hxgzt37vDdd9/x3XffZXnMhxvomGLUqFE0atSIu3fvsmLFChYvXmy8rGbqsf/3v/+xceNG6tSpw3PPPUeLFi3o3r17tvc5c0tISAidOnVi7NixTJ06lSZNmhAaGkr37t0f21I8Li6Oe/fuGZ/b2Ng88fJwt27d6NatG/Hx8ezZs4d58+axaNEi2rVrx9GjRx/beMfOzg4PD48M69zc3Lh9+7bx+YULF/D29sbBwSFDueeee+6xcQGcPn0aVVWpUKFClq+nt46+e/cud+/eNa7XaDSZ4nrY/PnzmTx5MidOnMjwz0y5cuUylOvVqxeDBw/mwoUL+Pr6snTpUnQ6Ha+99lqGGI8fP57t8R79PX70GE/ycAIEcHFxAcDHxyfL9emfffp38tHvHhi+t+vWrSMxMZGEhATu3buX5WdcqVIlY1KGnNeHeHaSbAuZ7t27079/f65du0arVq3y7H6LXq8HoGfPntneu6tRo8ZT7bt69eo0a9YMgNDQUJKSkujfvz8vvPACPj4+Jh27SpUqnDx5ktWrV7N27Vp+/fVXvv76a0aNGsXYsWMBsm0E9WgjIVMoisKyZcvYvXs3v/32G+vWreP1119n8uTJ7N69Gycnpyy3e/fdd5k/f77xeUhIiLExzZM4OzvTvHlzmjdvjlarZf78+ezZs4eQkJBst9FoNCa9L1Pp9XoUReGPP/7I8ljpn8OkSZOM9QGGrkrZ9V1esGABffr0ITQ0lGHDhuHp6YlGoyEiIiJTI7quXbvy/vvvs3DhQkaOHMmCBQsIDg7OkMD0ej3Vq1dnypQpWR7v0aRob2+fo/eeLrvPOLv16iMNFnNTTutDPDtJtoVMhw4deOutt9i9ezdLlizJtpyvry8bN24kISEhw9ntiRMnjK+nP+r1ev79998Mf5BOnjyZYX/pLZXT0tKMidFcJkyYwIoVKxg3bhyRkZEmH9vR0ZGwsDDCwsJISUmhY8eOjBs3jhEjRmBnZ4ebmxt37tzJtF36mcXjPGl0nnr16lGvXj3GjRvHokWL6NGjB4sXL+aNN97IsvyHH35obLgDGRuFmSI4OJj58+cTHR39VNs/zNfXly1btpCUlJTh7PbMmTNP3Nbf3x9VVSlXrhwVK1bMtlyvXr0y9Bl/XEJbtmwZ5cuXZ/ny5Rk+/9GjR2cqW7x4cdq0acPChQvp0aMHO3bsYNq0aZli/Pvvv3nppZdyffSrZ5H+nXz0uweG7627uzuOjo7Y2dlhb2/P6dOnM5V7dNuc1od4dnLPtpBxcnLim2++YcyYMbRr1y7bcq1btyYtLY2ZM2dmWD916lQURTG2aE5/fLQ186N/oDQaDZ06deLXX3/l6NGjmY5348aNp3k7WfL396dTp07MmzePa9eumXTsmzdvZnjNxsaGgIAAVFU1Xn709/cnLi6Ow4cPG8tFR0ezYsWKJ8bm6OgIkClZ3759O9MZSmBgIECm7iQPCwgIoFmzZsYlKCgo27JJSUns2rUry9fS78FndQnSVC1btkSn0zF79mzjOr1ez6xZs564bceOHdFoNIwdOzbT56GqqrF+ypcvn+F9P+4yf/oZ2cP727NnT7afxWuvvcaxY8cYNmwYGo2Grl27Zni9S5cuXLlyJcP7S3fv3j0SExOf+D7Nwdvbm8DAQObPn5/h9+vo0aOsX7+e1q1bA4bPo2XLlqxcuZKLFy8ayx0/fpx169Zl2GdO60M8OzmzLYRy0gWjXbt2NG3alI8++ojz589Ts2ZN1q9fz6pVq3jvvfeM92gDAwPp1q0bX3/9NXFxcTRo0IBNmzZleRYzYcIEtmzZQt26denfvz8BAQHcunWLAwcOsHHjRm7dupVr73HYsGH88ssvTJs2jQkTJuT42C1atKBkyZI0bNgQLy8vjh8/zsyZM2nTpo3xDL9r167873//o0OHDgwZMoSkpCS++eYbKlas+MRGXoGBgWg0Gr744gvi4uKwtbXlxRdfZNGiRXz99dd06NABf39/EhISmD17Ns7OzsY/ks8qKSmJBg0aUK9ePV5++WV8fHy4c+cOK1euZNu2bYSGhlKrVq1nPk5oaCh16tThgw8+4MyZM1SuXJn/+7//M37Gjzsb9Pf35/PPP2fEiBGcP3+e0NBQihUrxrlz51ixYgVvvvkmQ4cONSmetm3bsnz5cjp06ECbNm04d+4ckZGRBAQEZLjvm65NmzaUKFGCpUuX0qpVqwz9j8GQjH/55RfefvtttmzZQsOGDUlLS+PEiRP88ssvrFu3juDgYJNizC1ffvklrVq1on79+vTr18/Y9cfFxSXDOOVjx45l7dq1NGrUiIEDB5KammrsX/7wP5HmqA+RjTxv/yxy1cNdfx7n0a4/qqqqCQkJ6vvvv6+WKlVK1Wq1aoUKFdQvv/zS2AUi3b1799QhQ4aoJUqUUB0dHdV27dqply5dytT1R1VVNSYmRh00aJDq4+OjarVatWTJkupLL72kfvfdd8Yypnb9Wbp0aZavN2nSRHV2dlbv3LmT42N/++23auPGjdUSJUqotra2qr+/vzps2DA1Li4uw77Xr1+vVqtWTbWxsVErVaqkLliwIEddf1RVVWfPnq2WL1/e2M1iy5Yt6oEDB9Ru3bqpZcuWVW1tbVVPT0+1bdu26l9//fXYz8AUOp1OnT17thoaGqr6+vqqtra2qoODg1qrVi31yy+/zNA9JruuP46Ojpn2m9X7vnHjhtq9e3e1WLFiqouLi9qnTx91x44dKqAuXrz4sduqqqH70AsvvKA6Ojqqjo6OauXKldVBgwapJ0+efOL7fLTrj16vV8ePH298z7Vq1VJXr1792K4qAwcOVAF10aJFWb6ekpKifvHFF2rVqlVVW1tb1c3NTQ0KClLHjh2b4XcFUAcNGvTEmNNl9T3Mbj/pdfTll19mWL9x40a1YcOGqr29vers7Ky2a9dOPXbsWKZ9/vnnn2pQUJBqY2Ojli9fXo2MjHym+pCuP89GUVUz3n0XQhQZK1eupEOHDmzfvt3srbuf1fvvv88PP/zAtWvXMrWqFsIc5J6tEMJkD3dHAkNL7RkzZuDs7Mzzzz9voahy5v79+yxYsIBOnTpJohV5Ru7ZCiFM9s4773Dv3j3q169PcnIyy5cvZ+fOnYwfP97krjB55fr162zcuJFly5Zx8+ZN3n33XUuHJIoQSbZCCJO9+OKLTJ48mdWrV3P//n2ee+45ZsyYweDBgy0dWraOHTtGjx498PT0ZPr06cbW4ELkBblnK4QQQpiZ3LMVQgghzEySrRBCCGFmRe6erV6v5+rVqxQrVixfDcUmhBAib6mqSkJCAqVKlcowuYk5FLlke/Xq1UwDiQshhCi6Ll26RJkyZcx6jCKXbNOH5Lt06RLOzs7ZltPpdKxfv54WLVrINFNFiNR70ST1XvTodDpWrlzJG2+8kWmqUXMocsk2/dKxs7PzE5Otg4MDzs7O8uUrQqTeiyap96Invc7hybN15QZpICWEEEKYmSRbIYQQwswk2QohhBBmVuTu2QohhLmlpaWh0+ksHUaRp9Vq0Wg0lg4DkGQrhBC5RlVVrl27xp07dywdinjA1dWVkiVLWnxcBUm2Tyk67h7nYhMp5+6It0v+nOVECJG30hOtp6cnDg4OFv8DX5SpqkpSUhLXr18HwNvb26LxSLJ9Ckv2XWTE8iPoVbBSIKJjdcJql7V0WEIIC0pLSzMm2hIlSlg6HAHG6R6vX7+Op6enRS8pSwMpE0XH3TMmWgC9CiOXHyU67t7jNxRCFGrp92hlQvr8Jb0+LH0PXZKtic7FJhoTbbo0VeV8bJJlAhJC5Cty6Th/yS/1IcnWROXcHbF6pO40ioKfu/w3K4QQImtyz9ZE3i72TGvlQa1N3UlU7UjCjlKeJSj5x0KwcXxocTIstumPxQyL8Wdnw6O1jaXfkhBCZCsqKoqmTZty+/ZtXF1dLR1OgSXJ9im0r1wMNt+A9DPc2NMQ+5Q709gakq6dsyEB2zmDnQvYuhgeH17sXR88uoGdq+FRa5c7b0oIIYTZSLJ9Gq6+0G8jpNwFXRKkJEJyQsafUxIfLHcfPH/wmHwXkuMNZQHSkiEpGZKeMltb2xuSrr0bOBR/6LE4OJTIuDiWAAd3w5l3PrmPIYQQRYEk26dh4wA+tZ9tH/q0B8k33vB4P97w8/14SI6D++lL/IPHO3DvzkOPcaCmQeo9SLgHCVdzfmxrO3D0AEf3B48PFidPcPJ66NHLcCYtiVmIPJeXffmTk5MZNmwYixcvJj4+nuDgYKZOnUrt2v/9nduxYwcjRozg1KlTBAYG8v3331OtWjUALly4wODBg9m+fTspKSn4+fnx5Zdf0rp1a7PGXZBIsrUUK43hsrC969Ntr6qGJH3v9oPlFiTdMvycdOvB85v/LYk3DWfPqfcNS9wlw/Ik1naGpFus5IPF27A4l/pvKVZKLmcLkYvyui//hx9+yK+//sr8+fPx9fVl4sSJtGzZkjNnzhjLDBs2jK+++oqSJUsycuRI2rVrx6lTp9BqtQwaNIiUlBS2bt2Ko6Mjx44dw8nJyWzxFkSSbAsqRXlwf9cZ3Hxzto2qGi5tJ954kIBvwN3rhsf0n+9eh7sxhsfkOENivnPBsDyOQwlwLg0uPuBSGlzKGH52LWtYHD3kDFmIHMiuL3/jih5mOcNNTEzkm2++Yd68ebRq1QqA2bNns2HDBn744Qfj2e3o0aNp3rw5APPnz6dMmTKsWLGCLl26cPHiRTp16kT16tUBKF++fK7HWdBJsi1KFMXQOtrWCYqXe3J53T1IuGZIvgnXHixXDY/xV/9bUu/9dwZ97XDW+7K2e5B4fQ3/HLj5PVjKGR5t5b9gIeDxffnNkWz//fdfdDodDRs2NK7TarXUqVOH48ePG5Nt/fr1ja8XL16cSpUqcfz4cQCGDBnCgAEDWL9+Pc2aNaNTp07UqFEj12MtyCTZiuxp7Q1J+XGJWVUNl67jr0DcFYi/DHGX4c4lw2PcpQcJ+T7EnjIsWXHyguLlobg/lEhfnjM8l0vUoghJ78v/cMLN733533jjDVq2bMnvv//O+vXriYiIYPLkybzzzjuWDi3fkGQrno2iGFo/OxSHktWzLpOaYkjCty/AnYtw+/xDyzlDsr4bY1gu7nr0AODqA+4V/1s8KoNHJcMxhShkvF3siehYnZHLj5KmqmgUhfEdq5mtkZS/vz82Njbs2LEDX1/DLSmdTse+fft47733jOV2795N2bKG+8a3b9/m1KlTVKlSxfi6j48Pb7/9Nm+//TYjRoxg9uzZkmwfIslWmJ+1zYOz1mzu49y7DbfOwa2zcPNfuPUv3DxjWO7HGRL0nYtwZmPG7Rw9wSsAPAPAswp4VjU82uTfMwAhciKsdlkaV/TgfGwSfu4OZm2N7OjoyIABAxg2bBjFixenbNmyTJw4kaSkJPr168fff/8NwKeffkqJEiXw8vLio48+wt3dndDQUADee+89WrVqRcWKFbl9+zZbtmzJkIiFJFuRH9i7QWk3KP18xvWqComxhkvPN0/DjVMQexJunDRcnk68Dmevw9mohzZSDJegvapCyRqGxbuGoSW1EAWIt4t9nk3fOWHCBPR6Pa+99hoJCQkEBwezbt063NzcMpR59913OX36NIGBgfz222/Y2BhGwEtLS2PQoEFcvnwZZ2dnXn75ZaZOnZonsRcUkmxF/qUo4ORhWPwaZnwtOcGQdK8ff7D8AzHHDAk4/az42Kr/yjt5gXdN8A6EUoFQOkgSsBAP2NnZMX36dKZPn57ptSZNmqCqhhvIbdu2zXL7GTNmmDW+wkCSrSiYbItBmWDD8rC71yHmKFw7amgZfe2I4cz4bgycXm9Y0hUrZTibLhMMZWpDqVqgyFjVQojcJ8lWFC5OnuD0Ivi/+N+6lCRDAr56CKIPwdWDcOOEoRvTiatwYrWhnKLB2jOA6mleKP/ch3INDf2FhRDiGUmyFYWfjQP41DEs6ZLvQvTfcGU/XN5nWBKiUWKOUJ4jsPJBYywXH/Bt8GBpaOiOJINzCCFMJMlWFE22Tob7wA/fC467Qur5nVzY/gvlNDFYxRw1NMQ6vMSwADiVBL8XoFwjKNc4+xbWQgjxEEm2QqRzKY0aEMrR8zaUbd0aK32y4Yz3wk7Dcnkf3L0GR5cZFjCMiFW+yX+L9P0VQmRBkq0Q2bF1Av+mhgVAdx8u74Xz2+HcVkPyvXMBDsw3LIqVoZWz/0vwXDND4ysrjWXfgxAiX5BkK0ROae0Ml47LNYamIw33fS/shLNb4N/NhkZX6fd//5xgmJyhQgvD8txLhukKhRBFkpWlAwCYNWsWfn5+2NnZUbduXfbu3Ztt2Xnz5qEoSobFzk7GzhUWYOsEFVvAyxEwaA+8/w+0nwEBr4Cts2Fihr9/hmV9YaI//BgKe2cbxpAWQhQpFj+zXbJkCeHh4URGRlK3bl2mTZtGy5YtOXnyJJ6enllu4+zszMmTJ43PFWkdKvIDlzLwfC/DkqYzjPN8ap1huXnacAZ8dgusGWq43FylPQS0l0ZWQhQBFj+znTJlCv3796dv374EBAQQGRmJg4MDc+bMyXYbRVEoWbKkcfHy8srDiIXIAY3WcLm55Th45y8Y/Bc0Gws+dQHF0OVo42iYXgsiX4BtUwwTMwiRD0RFRaEoCnfu3Hmm/SQlJdGpUyecnZ2N+/Pz82PatGm5EmdBYtEz25SUFPbv38+IESOM66ysrGjWrBm7dj06+8t/7t69i6+vL3q9nueff57x48dTtWrVLMsmJyeTnJxsfB4fHw8YZrXQ6XTZHiP9tceVEYWP2erdxQ/qDjIsCdewOrUG5cRqlAs7UK4dMYx0tWkseu9aqNU6oQ/oYBhiUuSJ3Kh3nU6Hqqro9Xr0en1uhZYnXnzxRWrWrGkcz7hevXpcuXKFYsWKPdN7mTt3Ltu2bWP79u24u7tTrFgxAOPnlBf0ej2qqqLT6dBo/muwmNd/2y2abGNjY0lLS8t0Zurl5cWJEyey3KZSpUrMmTOHGjVqEBcXx6RJk2jQoAH//PMPZcpkHu0nIiKCsWPHZlq/fv16HByePDvMhg0bcvhuRGFi/novCW5vYFMsDO87f1H69h7c7x7HKvogRB/EasMn3ChWlUvFGxLtEkyaxtbM8Qh4tnq3tramZMmS3L17l5SUlFyMyvxSU1NJSUkxnowAODg4kJCQ8Ez7PXHiBBUqVDBOzZeQkIBer+f+/fsZjmVOKSkp3Lt3j61bt5Kamponx8yKxe/Zmqp+/frUr1/f+LxBgwZUqVKFb7/9ls8++yxT+REjRhAeHm58Hh8fj4+PDy1atMDZ2Tnb4+h0OjZs2EDz5s3RarW5+yZEvmWZeg8DIPXudayOr0I5ugyrq/vxTDiKZ8JRVJsFqJVfQV+zK6pPfRnBygxyo97v37/PpUuXcHJyKlCNNvv27cuOHTvYsWMHkZGRAPzwww/069ePmzdv4urqyrx58wgPD+fHH39k2LBhXLp0iVatWjF//nyWLl3K2LFjiYuLo2fPnkyZMgWNRsOLL77In3/+CYCbmxshISFs3rwZKysr7OzsjH9/L168yJAhQ4yvtWzZkunTp+Pl5UVcXBzu7u7s2rWL4OBg9Ho9np6eVKxYkZ07dwKwYMECPvroIy5cuJDl+7t//z729vY0btw4Q73odDpWrVqV5TbmYNFk6+7ujkajISYmJsP6mJgYSpbM2YwsWq2WWrVqcebMmSxft7W1xdY281mBVqvN0Zcqp+VE4WKRencrDQ0GGpab/8KRpfD3zyi3z6McXoTV4UVQ3B+efw1qdodicpk5tz1LvaelpaEoClZWVlhZPWgOo6qgS8rFCE2gdcjRP2bTp0/n9OnTVKtWjU8//RSAf/75B8D4XqysrEhKSmLmzJksXryYhIQEOnbsSKdOnXB1dWXNmjWcPXuWTp068cILLxAWFsby5csZPnw4R48eZfny5djY2Bg/l/TPSa/X06FDB5ycnPjzzz9JTU1l0KBBdOvWjaioKNzc3AgMDGTr1q3UqVOHI0eOoCgKBw8eJCkpCScnJ7Zt20ZISMh/n/kjrKysUBTF4n/LLZpsbWxsCAoKYtOmTcZJiPV6PZs2bWLw4ME52kdaWhpHjhyhdevWZoxUiDxWwh+aDIeQ/xlaNR9aBP+sgFv/wsYxsOkzqNQKgl+H8k0hmz80wsJ0STC+lGWOPfIq2Dg+sZiLiws2NjY4ODgYT3Kyuo2n0+n45ptv8Pf3B6Bz58789NNPxMTE4OTkREBAAE2bNmXLli2EhYVRvHhxHBwcsLGxyfbkadOmTRw5coRz587h4+MDwI8//kjVqlXZt28ftWvXpkmTJkRFRTF06FCioqJo3rw5J06cYPv27bz88stERUXx4YcfPu2nlGcs/g0NDw9n9uzZzJ8/n+PHjzNgwAASExPp27cvAL169crQgOrTTz9l/fr1nD17lgMHDtCzZ08uXLjAG2+8Yam3IIT5KIphEoRXZsIHJ+GVWVCmDqhphtmKFnSEGc/Djq8g6ZaloxWFmIODgzHRgqFtjZ+fH05OThnWXb9+Pcf7PH78OD4+PsZECxAQEICrqyvHjx8HICQkhO3bt5OWlsaff/5JkyZNjAn46tWrnDlzhiZNmjz7GzQzi9+zDQsL48aNG4waNYpr164RGBjI2rVrjY2mLl68mOHywO3bt+nfvz/Xrl3Dzc2NoKAgdu7cSUBAgKXeghB5w9YJavU0LNePw19zDYNm3D4HG0bBlvFQowvUeQtKVrN0tAIMl3JHXrXcsXNzd49cgk2/NPvoutxuZdy4cWMSEhI4cOAAW7duZfz48ZQsWZIJEyZQs2ZNSpUqRYUKFXL1mOZg8WQLMHjw4GwvG0dFRWV4PnXqVGPzdCGKLM8q0HoiNBsNR5bBvtmG7kMHfjQsfo2g/iCo0FIuMVuSouToUq6l2djYkJaWlufHrVKlCpcuXeLSpUvGs9tjx45x584d4wmUq6srNWrUYObMmWi1WipXroynpydhYWGsXr2akJCQPI/7aci3UIiCzMYRgnrDW9ug71qo2gEUDZzfBj93hVl14K85oLtn6UhFPubn58eePXs4f/48sbGxedYHtlmzZlSvXp0ePXpw4MAB9u7dS69evQgJCSE4ONhYrkmTJixcuNCYWIsXL06VKlVYsmSJJFshRB5SFPCtD6/Og/cOQ8N3wdbFMEzk6vdhWnXYNhnu3bF0pCIfGjp0KBqNhoCAADw8PLh48WKeHFdRFFatWoWbmxuNGzemWbNmlC9fniVLlmQoFxISQlpaWoZ7s02aNMm0Lj9TVFVVLR1EXoqPj8fFxYW4uLgn9rNds2YNrVu3lq4/RUihqvfkBDi4AHZ9DXEP/njaOhtaMNcfDE4elo0vH8mNer9//z7nzp2jXLlyBaqfbWGXXb3odDqWLVtG9+7dn5gPcoOc2QpRWNkWg3oDYMgB6PAdeFSB5HjYMc1wprt2JCRcs3SUQhQJkmyFKOw0WqgZBgN2QtefDTMOpd6D3bNgWg3443+QEPPk/QghnpokWyGKCisrqNwa3tgEPX81zECUlgx7IuGrmrD+E0i8aekohSiUJNkKUdQoCjzXDF5fB6+thDK1DWe6O6fDVzUgaoLhfq8QItdIshWiqFIU8G8K/TZA96VQsgak3IWoCMM8u3u+g9SCNXtNflDE2pzme/mlPiTZClHUKQpUbAFvbTV0HSpeHhJvwB/DSJ1RmxNbFhJ9x0KD6Rcg6a2Yk5Lks8pP0uvD0r0L8sUIUkKIfEBRDINiVG4LB37k3sZx2Medp/KfA9m7ZTpHG42meXOZ8CM7Go0GV1dX49jADg4OKDIdosWoqkpSUhLXr1/H1dU1w8TxliDJVgiRkUZLdMXuNF/uSn/Nb7yp+Z06VidgRzeSbnfBodVnUCxnU2AWNemz25gyGL8wL1dX1xxP2WpOkmyFEJmci03krmrH1NRX+Tn1RYZpl9BJsx2HY7/AmT8gZBjUHQDWNpYONV9RFAVvb288PT3R6XSWDqfI02q1Fj+jTSfJVgiRSTl3R6wU0KtwjRJ8oBvIwrSWLC67HJtrBw2zDB34EdpMhvJNLB1uvqPRaPLNH3mRP0gDKSFEJt4u9kR0rI7mwT1HjaIQ1iEUmzc3wytfg6MH3DwDP74Cy/rJSFRCPIGc2QohshRWuyyNK3pwPjYJP3cHvF3sDS/U6gGV28CWcbDvezi6DE6vh5dGQXA/mdJPiCzIt0IIkS1vF3vq+5f4L9Gms3eF1l9C/81Q6nnDmMtrhsKcloaJ7YUQGUiyFUI8vVK14I2N0HoS2BSDy3shshFsHgepyZaOToh8Q5KtEOLZWGmgTn8YtAcqtQa9DrZOhG9D4MoBS0cnRL4gyVYIkTtcSkPXRfDqfEMDqhvH4ftmsHGsnOWKIk+SrRAi9ygKVA2FgXugWidQ02D7FMNZbvTflo5OCIuRZCuEyH2OJaDzHAhb8N9Z7uyXYNtk0KdZOjoh8pwkWyGE+VRpBwN3G8Zb1utg06cwtxXcPm/pyITIU5JshRDm5ehuOMMN/cbQYvnSHvjmBTi81NKRCZFnJNkKIcxPUSCwOwzYAT51ISUBlr8By9+E+/GWjk4Is8vRCFLPP/+8STtVFIX/+7//o3Tp0k8VlBCikHLzhT5rYOuXhu5Bh5cYznRfnWfosytEIZWjZHvo0CE++OADnJycnlhWVVUmTJhAcrI09RdCZEFjDU1HGCYwWN7fcP/2hxbQYpyhv67MASsKoRyPjTxs2DA8PT1zVHby5MlPHZAQoojwrQ9vbYVVg+DkGvhjGJzfBq/MBDsXS0cnRK7K0T3bc+fO4eHhkeOdHjt2DF9f36cOSghRRDgUNwyE0TICrLRw/P/gu6YQ84+lIxMiV+Uo2fr6+qKYcGnHx8dH5nIUQuSMokD9gfD6OnDxgVv/Gvrk/r3E0pEJkWueaoq9+/fvc/jwYa5fv45er8/wWvv27XMlMCFEEVMmCN7809BK+d/NsOJNw8QGLSPA2sbS0QnxTExOtmvXrqVXr17ExsZmek1RFNLSZHQYIcRTciwBPZZB1ARDa+V930PMMegyH5xy1mZEiPzI5H6277zzDq+++irR0dHo9foMy9Mm2lmzZuHn54ednR1169Zl7969Odpu8eLFKIpCaGjoUx1XCJEPWWngxY+g22LDIBgXd8J3TWQGIVGgmZxsY2JiCA8Px8vLK1cCWLJkCeHh4YwePZoDBw5Qs2ZNWrZsyfXr1x+73fnz5xk6dCiNGjXKlTiEEPlMpVaGyelLPAfxV2DOy3D4F0tHJcRTMTnZdu7cmaioqFwLYMqUKfTv35++ffsSEBBAZGQkDg4OzJkzJ9tt0tLS6NGjB2PHjqV8+fK5FosQIp/xqGhIuBVfhrRkQ7/cTZ/BI21FhMjvTL5nO3PmTF599VW2bdtG9erV0Wq1GV4fMmRIjveVkpLC/v37GTFihHGdlZUVzZo1Y9euXdlu9+mnn+Lp6Um/fv3Ytm2bqW9BCFGQ2LlA159h01jYMQ22TYLYk9DhW7BxtHR0QuSIycn2559/Zv369djZ2REVFZWhS5CiKCYl29jYWNLS0jJdkvby8uLEiRNZbrN9+3Z++OEHDh06lKNjJCcnZxjNKj7eMA6rTqdDp9Nlu136a48rIwofqfd8rMnHKMUroFnzPsrx31BvnSc1bBEU837mXUu9Fz15XdcmJ9uPPvqIsWPHMnz4cKys8nYeg4SEBF577TVmz56Nu7t7jraJiIhg7NixmdavX78eBweHJ26/YcMGk+MUBZ/Ue35VDK3P/2hycRoOMUdI/SaEXf4fkGDvkyt7l3oX5mJysk1JSSEsLCxXEq27uzsajYaYmJgM62NiYihZsmSm8v/++y/nz5+nXbt2xnXp/Xytra05efIk/v7+GbYZMWIE4eHhxufx8fH4+PjQokULnJ2ds41Np9OxYcMGmjdvnulSuSi8pN7zt6X7L/PxbmtK8RnztBN5TneVpmcjSOs4B9X/xafer9R70aPT6Vi1alWeHc/kZNu7d2+WLFnCyJEjn/ngNjY2BAUFsWnTJmP3Hb1ez6ZNmxg8eHCm8pUrV+bIkSMZ1n388cckJCTw1Vdf4eOT+b9bW1tbbG1tM63XarU5+lLltJwoXKTe85/ouHt8vOoYehUu40nHlDF8ZzONeinHsF7SDdpOhaDez3QMqXdhLiYn27S0NCZOnMi6deuoUaNGpl/MKVOmmLS/8PBwevfuTXBwMHXq1GHatGkkJibSt29fAHr16kXp0qWJiIjAzs6OatWqZdje1dUVINN6IUThci42Eb363/N4nHgtZTg7A1bgcXYF/DYEEq5ByIcyc5DId0xOtkeOHKFWLcO8k0ePHs3wminjJ6cLCwvjxo0bjBo1imvXrhEYGMjatWuNjaYuXryY5/eGhRD5Tzl3R6wUMiRcvaIltf3XsN/f0Eo5ajwkXIXWkw1T+QmRT5j827hly5ZcD2Lw4MFZXjYGntind968ebkejxAi//F2sSeiY3VGLj9KmqqiURTGd6yGt6sDvPQJFCsJa4bB/nlw9zp0ngNae0uHLQTwlBMRCCGEJYTVLkvjih6cj03Cz90Bb5eHkmmd/uDkBb++YZgfd0En6PazzI0r8oUcXZ/t2LGjsX9qTvTo0eOJwy0KIcTT8Haxp75/iYyJNl1Ae3htBdg6w4UdMK8t3L2R90EK8YgcJdtVq1Zx48YN4uPjn7jExcXx22+/cffuXXPHLoQQmfk1hD6rwdEDrh2GOS3h9gVLRyWKuBxdRlZVlYoVK5o7FiGEyB3eNQ2T0f8YapiMfm4r6LUK3CtYOjJRROUo2T5No6jSpUubvI0QQuSaEv7Q70HCjT1pSLivrYSS0k1Q5L0cJduQkBBzxyGEELnPuRT0XQM/hcK1IzCvDby2HEoHWToyUcRIB1YhROHm6A69V0OZ2nD/Dsx/BS7utnRUooiRZCuEKPzsXQ2tlP0aQUoC/NQRzu+wdFSiCJFkK4QoGmyLQfdfoHwT0CXCws5wTubDFnlDkq0QouiwcYBui8H/JdAlwcJX4WyUpaMSRYAkWyFE0aK1h66LoEILSL0Hi8JQzm+1dFSikDM52cbExPDaa69RqlQprK2t0Wg0GRYhhMj3tHYQtgAqvgyp99Es6UGJhBOWjkoUYiaPjdynTx8uXrzIJ598gre391PN9COEEBZnbQtdfoTFPVDObKDe2clwsT74N7Z0ZKIQMjnZbt++nW3bthEYGGiGcIQQIg9Z20LYAvQ/d8X67BbUJV0NA1/41LF0ZKKQMfkyso+PD6qqPrmgEEIUBFo70jr/yPViVVFSEkn9sQM3Tkk/XJG7TE6206ZNY/jw4Zw/f94M4QghhAVo7Znh+D579JWx1t3FemEn1m7eZOmoRCFi8mXksLAwkpKS8Pf3x8HBAa1Wm+H1W7du5VpwQgiRF6Lj7rPgvD3LGcYCmwhqWZ0h6M++XC/3O57lqls6PFEImJxsp06dKo2ihBCFyoWbSagoJGJP75QP+dlmHFWtLpC8tDP0Xw9uvpYOURRwT9UaWQghChPfEg4oqKgoxOPEaykjWGLzGRWSrhgmMei7Fop5WTpMUYCZfM+2V69ezJ07l3///dcc8QghRJ7zdrEjrLweqwcX7eIUF443nw+uZeHWWfipA9y7bdkgRYFm8pmtjY0NERER9OvXj9KlSxMSEkKTJk0ICQmhQgWZmFkIUTDV91IZ2LExV+JS8HN3wNvFHqqugjkvw/V/YGEXw2QGtk6WDlUUQCaf2X7//fecOnWKS5cuMXHiRJycnJg8eTKVK1emTJky5ohRCCHyhLeLHfX9SxgSLUDx8oYEa+cKl/fCL69BaopFYxQF01OPjezm5kaJEiVwc3PD1dUVa2trPDw8cjM2IYSwPK+q0PNX0DrAv5th5QDQ6y0dlShgTE62I0eOpEGDBpQoUYLhw4dz//59hg8fzrVr1zh48KA5YhRCCMsqEwxhP4GVNRxdButGgAzuI0xg8j3bCRMm4OHhwejRo+nYsSMVK1Y0R1xCCJG/PNcMQiNh+RuwJxIc3aHxMEtHJQoIk89sDx48yEcffcTevXtp2LAhpUuXpnv37nz33XecOnXKHDEKIUT+UONVePkLw8+bP4cDP1k2HlFgmJxsa9asyZAhQ1i+fDk3btxgzZo12NjYMGjQIKpUqWKOGIUQIv+o9za8EG74+bd34dR6y8YjCgSTLyOrqsrBgweJiooiKiqK7du3Ex8fT40aNQgJCTFHjEIIkb+8NAoSouHvn2Fpb+izGkoHWToqkY+ZnGyLFy/O3bt3qVmzJiEhIfTv359GjRrh6upqhvCEECIfUhRoPwPuXod/Nxn64PZbDyX8LR2ZyKdMTrYLFiygUaNGODs7myMeIYQoGDRa6DIf5rWB6L9hYWfot8HQcEqIR5h8z7ZNmzbGRHv58mUuX76c60EJIUSBYFsMui81DuuYsqALu09eJjrunqUjE/mMyclWr9fz6aef4uLigq+vL76+vri6uvLZZ5+hl47eQoiippgX9FhGinUxbKL3c3tBHxpN2MiSfRctHZnIR0xOth999BEzZ85kwoQJHDx4kIMHDzJ+/HhmzJjBJ5988lRBzJo1Cz8/P+zs7Khbty579+7Ntuzy5csJDg7G1dUVR0dHAgMD+eknaX4vhLCcaJuyvJb0HsmqNa00+xiuWcjI5UflDFcYmXzPdv78+Xz//fe0b9/euK5GjRqULl2agQMHMm7cOJP2t2TJEsLDw4mMjKRu3bpMmzaNli1bcvLkSTw9PTOVL168OB999BGVK1fGxsaG1atX07dvXzw9PWnZsqWpb0cIIZ7ZudhE9uirMEz3NtNtZvKG9R9cVD05H1v3v3GWRZFm8pntrVu3qFy5cqb1lStX5tatWyYHMGXKFPr370/fvn0JCAggMjISBwcH5syZk2X5Jk2a0KFDB6pUqYK/vz/vvvsuNWrUYPv27SYfWwghckM5d0esFPg/fQMm6sIAGG39I5Xid1g4MpFfmHxmW7NmTWbOnMn06dMzrJ85cyY1a9Y0aV8pKSns37+fESNGGNdZWVnRrFkzdu3a9cTtVVVl8+bNnDx5ki+++CLLMsnJySQnJxufx8fHA6DT6dDpdNnuO/21x5URhY/Ue9H0rPXu7mDN568E8PGqY3yd1h4/qxi6aKJw++NtdB4+4FUtF6MVuSGvv+MmJ9uJEyfSpk0bNm7cSP369QHYtWsXly5dYs2aNSbtKzY2lrS0NLy8vDKs9/Ly4sSJE9luFxcXR+nSpUlOTkaj0fD111/TvHnzLMtGREQwduzYTOvXr1+Pg4PDE2PcsGHDE8uIwkfqvWh6lnp3BEbXghv3FZJte3HjUgwed4+TOr8jWyuN5r7WLfcCFQWOyck2JCSEU6dOMWvWLGNC7NixIwMHDqRUqVK5HmBWihUrxqFDh7h79y6bNm0iPDyc8uXL06RJk0xlR4wYQXh4uPF5fHw8Pj4+tGjR4rF9hXU6HRs2bKB58+ZotVpzvA2RD0m9F01mqfd7jVHnt8L+5mlaxP5A6mu/gY1j7uxbPDOdTseqVavy7HgmJ1uAUqVKmdwQKivu7u5oNBpiYmIyrI+JiaFkyZLZbmdlZcVzzz0HQGBgIMePHyciIiLLZGtra4utrW2m9VqtNkdfqpyWE4WL1HvRlKv1rvWAHkvh+5dQrh1Gu3owvPojWD31NOKiAMtRsj18+HCOd1ijRo0cl7WxsSEoKIhNmzYRGhoKGPrxbtq0icGDB+d4P3q9PsN9WSGEyBeKl4OwhfBjezj+G2z53DCusihycpRsAwMDURQFVVVRFMW4Xn0wefLD69LS0kwKIDw8nN69exMcHEydOnWYNm0aiYmJ9O3bF4BevXpRunRpIiIiAMM92ODgYPz9/UlOTmbNmjX89NNPfPPNNyYdVwgh8oRvfcM4yivegm2ToUQFCOxm6ahEHstRsj137pzx54MHDzJ06FCGDRuWoYHU5MmTmThxoskBhIWFcePGDUaNGsW1a9cIDAxk7dq1xkZTFy9exOqhyy6JiYkMHDiQy5cvY29vT+XKlVmwYAFhYWEmH1sIIfJEza4Qe8qQbP/vHXDzA9/6RMfd41xsIuXcHaU/biGXo2Tr6+tr/PnVV19l+vTptG7d2riuRo0a+Pj48MknnxgvB5ti8ODB2V42joqKyvD8888/5/PPPzf5GEIIYVFNP4bY03D8/2BJT36ru5B3195Er4KVAhEdqxNWu6yloxRmYvKd+iNHjlCuXLlM68uVK8exY8dyJSghhCh0rKygQySUrAFJsTy3qT926n0A9CoyvGMhZ3KyrVKlChEREaSkpBjXpaSkEBERQZUqVXI1OCGEKFRsHKHbz6TYuVPF6iJTtV+jYJjAJU1VOR+bZOEAhbmY3PUnMjKSdu3aUaZMGWPL48OHD6MoCr/99luuByiEEIWKSxniQ+dT7OdQWmr+4gN1KZNSw9AoCn7uTx5oRxRMJifbOnXqcPbsWRYuXGgc1CIsLIzu3bvj6CgdtoUQ4kncK7/A7sCx1Pt7JIOtV3FaLUuD0DelkVQh9lSDWjg6OvLmm2/mdixCCFFk1OswiLvaKzj9NYtp9t+jlO4ESAOpwsrke7Zly5alV69e/PDDD5w9e9YcMQkhRJHg1PozqNACJfUeLO4BCTFP3kgUSCYn2/Hjx2NnZ8cXX3zBc889h4+PDz179mT27NmcPn3aHDEKIUThZKWBTt+De0WIvwK/vAapMhpeYWRysu3Zsyffffcdp06d4sqVK3z55ZcADBw4MMt5boUQQjyGnQt0/dnweGkP/P4BPBidTxQeT3XPNikpie3btxMVFcWWLVs4ePAg1apVy3IiACGEEE/g/hx0ngsLO8PBn8C7JtTpb+moRC4yOdk2aNCAgwcPUqVKFZo0acLw4cNp3Lgxbm4yV6MQQjy1516CZmNhwyewdjh4VgG/FywdlcglJl9GPnHiBI6OjlSuXJnKlStTpUoVSbRCCJEbGrwD1buAPhV+6QV3LgIQHXePnf/GyghTBZjJyfbmzZts3ryZevXqsW7dOho2bEjp0qXp3r07s2fPNkeMQghRNCgKtJ9uuIycdBMWd2fZ7pM0nLCZ7rP30HDCZpbsu2jpKMVTMDnZKopCjRo1GDJkCMuWLeOPP/6gefPmLF26lLffftscMQohRNGhtTfMgevoAdeOYPP7e+gfNJiSMZQLLpOT7YEDB5gyZQrt27enRIkS1K9fn8OHD/POO++wfPlyc8QohBBFi6sPvDofvWJNe81O+mt+N74kYygXTE81XGOtWrUICQmhf//+NG7cGBcXF3PEJoQQRZdfQxKafIrLlpEMt/6Z46ov2/XVZQzlAsrkZHvr1i2cnZ3NEYsQQoiHuDQeyNnT+yh/eQUztdMJTRnHgI4vyRjKBZDJl5El0QohRB5RFMr3jiSlZC1clUQ2lP6WsJolLB2VeAomJ9u0tDQmTZpEnTp1KFmyJMWLF8+wCCGEyEVaO2y6LwJHT7Sxx2HVYBlhqgAyOdmOHTuWKVOmEBYWRlxcHOHh4XTs2BErKyvGjBljhhCFEKKIcy4FXeaDlTX8sxx2zrB0RMJEJifbhQsXMnv2bD744AOsra3p1q0b33//PaNGjWL37t3miFEIIYRvA3h5guHnjaPh382WjUeYxORke+3aNapXrw6Ak5MTcXFxALRt25bff//9cZsKIYR4FrXfgMCeoOph2etw+7ylIxI5ZHKyLVOmDNHR0QD4+/uzfv16APbt24etrW3uRieEEOI/igJtJkOp5+HebVjSE1Kkz21BYHKy7dChA5s2bQLgnXfe4ZNPPqFChQr06tWL119/PdcDFEII8RCtHYT9BA7ucO0I/PYu0XeSZOzkfM7kfrYTJkww/hwWFoavry87d+6kQoUKtGvXLleDE0IIkQWXMoYGU/Pbw5FfmH3QljmprbBSIKJjdcJql7V0hOIRJp3Z6nQ6Xn/9dc6dO2dcV69ePcLDwyXRCiFEXvJ7gbjGYwAYqVlIPatjMnZyPmZSstVqtfz666/mikUIIYQJ/vHpxvK0F7BW9MzUTsebmzJ2cj5l8j3b0NBQVq5caYZQhBBCmKKchxMfp/bjH70v7ko8kTZTsVd0MnZyPmTyPdsKFSrw6aefsmPHDoKCgnB0dMzw+pAhQ3ItOCGEENnzdrFndMdgBiz/gFU2I6lpdZbf/Vfg7fyKpUMTjzA52f7www+4urqyf/9+9u/fn+E1RVEk2QohRB4Kq12WxhW7EH3EFddNfSh/eSX89YOhT67IN0xOtg83jhJCCGF53i72eL/wCihjYMMo+GM4eFWDsvUsHZp4wOR7tkIIIfKpBkOgagfQ6+CXXhAfTXTcPemDmw/k6Mw2PDw8xzucMmWKyUHMmjWLL7/8kmvXrlGzZk1mzJhBnTp1siw7e/ZsfvzxR44ePQpAUFAQ48ePz7a8EEIUGYoCr8yCGyfh+jFi54TRNCac+6pW+uBaWI6S7cGDBzM8P3DgAKmpqVSqVAmAU6dOodFoCAoKMjmAJUuWEB4eTmRkJHXr1mXatGm0bNmSkydP4unpmal8VFQU3bp1o0GDBtjZ2fHFF1/QokUL/vnnH0qXLm3y8YUQolCxcYSwBei/a4r7nb/5RPMjH6X2M/bBbVzRQyaft4AcXUbesmWLcWnXrh0hISFcvnyZAwcOcODAAS5dukTTpk1p06aNyQFMmTKF/v3707dvXwICAoiMjMTBwYE5c+ZkWX7hwoUMHDiQwMBAKleuzPfff49erzcOISmEEEVeCX9ONJyCXlXoYb2JMM0WAOmDa0EmN5CaPHky69evx83NzbjOzc2Nzz//nBYtWvDBBx/keF8pKSns37+fESNGGNdZWVnRrFkzdu3alaN9JCUlodPpsp24Pjk5meTkZOPz+Ph4wDAalk6ny3a/6a89rowofKTei6bCWO9OAS2Zsv5Vhlr/wqfWczmlL8PfVKC0i02hep9PK68/A5OTbXx8PDdu3Mi0/saNGyQkJJi0r9jYWNLS0vDy8sqw3svLixMnTuRoH//73/8oVaoUzZo1y/L1iIgIxo4dm2n9+vXrcXB4csfvDRs25CgOUbhIvRdNha3eY33a8sflc7TS7CPSZirfl/yMgzs2c/DJm4pcZnKy7dChA3379mXy5MnGRkl79uxh2LBhdOzYMdcDfJwJEyawePFioqKisLOzy7LMiBEjMjTwio+Px8fHhxYtWuDs7JztvnU6HRs2bKB58+Zotdpcj13kT1LvRVNhrffWwLUbdUhc0havuDOMsP6RtBYrwVqmQ9XpdKxatSrPjmdyso2MjGTo0KF0797deBpubW1Nv379+PLLL03al7u7OxqNhpiYmAzrY2JiKFmy5GO3nTRpEhMmTGDjxo3UqFEj23K2trZZzrOr1Wpz9KXKaTlRuEi9F02Fsd59SnlDr1/gu6ZYXdmH1caPoN1Xlg6ryDG5n62DgwNff/01N2/e5ODBgxw8eJBbt27x9ddfZxq68UlsbGwICgrK0LgpvbFT/fr1s91u4sSJfPbZZ6xdu5bg4GBT34IQQhQtJfyh8w+AAvvnwV+GBqjSBzfvmHxmm87R0fGxZ5Q5FR4eTu/evQkODqZOnTpMmzaNxMRE+vbtC0CvXr0oXbo0ERERAHzxxReMGjWKRYsW4efnx7Vr1wBwcnLCycnpmeMRQohCqUJzeOkT2PQprPmQTTdL0D/KGr2K9MHNAxYfQSosLIxJkyYxatQoAgMDOXToEGvXrjU2mrp48SLR0dHG8t988w0pKSl07twZb29v4zJp0iRLvQUhhCgYXgiHgFDQ66ix8x081ZsAMg9uHnjqM9vcNHjwYAYPHpzla1FRURmenz9/3vwBCSFEYaQoEPo1idEn8Lh9gm9tptIlZRTJ2Bj74MqAF+Zh8TNbIYQQecjGkcQO87mlOlHT6izjtd8DKhpFkXlwzUiSrRBCFDGeZSvzd71ppKpWdNJs503rNYzvWE3Oas1Ikq0QQhRBTVu9SmLTTwEYof2ZMNdTFo6ocJNkK4QQRZRLyGCo1RNF1cOy1yH2tKVDKrQk2QohRFGlKNBmCvjUheQ4+Lkb3Lsj/W/NIF+0RhZCCGEh1rYQtgC+awI3TxM9pzuNL7+NTtVI/9tcJGe2QghR1Dl5QtdFqNb2eN/Ywf80iwDpf5ubJNkKIYSAUoGcamAY3/4N6z/oInPg5ipJtkIIIQBwDurM1NTOAHxuPYc6ynHpf5tLJNkKIYQAwNvFnlLtR/F7Wj1slDQibabyVUsX6X+bCyTZCiGEMAqr48vz7y7ibonqFFfu0vboe3DvjqXDKvAk2QohhMjA270ETr2XgnNpiD0FS/tAms7SYRVokmyFEEJk5uwN3RaD1hHObiFxVTg7z9yQlslPSZKtEEKIrHnXgE7fo6LgePhHNs0dQ8MJm1my76KlIytwJNkKIYTIVrR3U8andgfgI+uFNFf2Sd/bpyDJVgghRLbOxSYyO7U1P6U2w0pRmaadRTVOS99bE0myFUIIka1y7o5YKQpjUnuzOS0QeyWFH2wm4a+9IWMom0CSrRBCiGx5u9gT0bE6KNa8o3uHf/S+uCvx2C3pSusJ/0f32XvkPm4OSLIVQgjxWGG1y7J9eFO+798U97dWkuZUCufEc3ynnYQtKTKGcg5IshVCCPFE3i721PcvgVfp8hxp8j3xqgO1rU4xTTsLK/QyhvITSLIVQghhEq8Kz/OWLpxk1ZpWmn2Msv4RjQIONlZyDzcbMp+tEEIIk3i72BPaIYyhK+OZoZ1OH+v1uJcqR4evDdPyyTy4mcmZrRBCCJOF1S7LyA9Hci74YwDaXv+WTlZRgMyDmxVJtkIIIZ6Kt4s95doO43LVtwCYYD2bFlb7AJkH91GSbIUQQjwTTfMx/JLWBI2iMkM7k7oyD24mkmyFEEI8E29XB2g7lfVpwdgqOr63mcQ3L1nJPLgPkWQrhBDimXWpW57q7y0jzqsuxZR7tNg/AK6fsHRY+YYkWyGEELnCu4QbLn2XQann4d4t+PEVuHXW0mHlC5JshRBC5B47Z+j5K3gGwN1r3P+hLTGXz1g6KouTZCuEECJ3ORRnVfVZnNOXxC7xCknftWHVtr+K9MQFMqiFEEKIXBUdd4/310RTUh3JL7afUs7qGvoNPXllzSdcV12L5KAXFj+znTVrFn5+ftjZ2VG3bl327t2bbdl//vmHTp064efnh6IoTJs2Le8CFUIIkSPnYhPRq3AVd7qmfMxl1R1/q2gWasfhTlyRHPTCosl2yZIlhIeHM3r0aA4cOEDNmjVp2bIl169fz7J8UlIS5cuXZ8KECZQsWTKPoxVCCJEThjlwDT9fVj3pnvIRV9XiVLC6wkKbcZQgrsgNemHRZDtlyhT69+9P3759CQgIIDIyEgcHB+bMmZNl+dq1a/Pll1/StWtXbG1t8zhaIYQQOZE+B65GMWTcy6oX3VM+5prqRiWryyy2+ZySyp0iNeiFxe7ZpqSksH//fkaMGGFcZ2VlRbNmzdi1a5elwhJCCJELwmqXpXFFD87HJuHn7sDWUzfosQJ+0n5OBasrbHCbQDFCgDKWDjVPWCzZxsbGkpaWhpeXV4b1Xl5enDiRex2hk5OTSU5ONj6Pj48HQKfTodPpst0u/bXHlRGFj9R70ST1bh7uDta4l3UGoGOgN/XLdeHKxZp4bulDsYSLqHNakdpzBbj65nlseV3Xhb41ckREBGPHjs20fv369Tg4PPkSxoYNG8wRlsjnpN6LJqn3vHG77PvUO/0FznEXSfm2GbsrfMhdu9KWDsusLJZs3d3d0Wg0xMTEZFgfExOTq42fRowYQXh4uPF5fHw8Pj4+tGjRAmdn52y30+l0bNiwgebNm6PVanMtHpG/Sb0XTVLveWvp/su0OGDLT9rxVEi9QsMzX3An9CdOW1fCt4QD3i52Zo9Bp9OxatUqsx8nncWSrY2NDUFBQWzatInQ0FAA9Ho9mzZtYvDgwbl2HFtb2ywbU2m12hx9qXJaThQuUu9Fk9S7+UXH3ePjVcfQq250SfmEuTYTCdSdpdgvnflG9wG71WqFsg+uRVsjh4eHM3v2bObPn8/x48cZMGAAiYmJ9O3bF4BevXplaECVkpLCoUOHOHToECkpKVy5coVDhw5x5owMBSaEEAVBeh9cgNs40yPlI7anVcVRSWaudiIvK7sLZR9ci96zDQsL48aNG4waNYpr164RGBjI2rVrjY2mLl68iJXVf/8PXL16lVq1ahmfT5o0iUmTJhESEkJUVFRehy+EEMJE6X1w0xNuIva8rvuQacyitWYvX9tM53PdTX7/uxJtapYqNNP0WXwEqcGDB3PhwgWSk5PZs2cPdevWNb4WFRXFvHnzjM/9/PxQVTXTIolWCCEKhkf74FoBOrQM1g1hXmoLAD7WLkSzfgSNJmxkyb6LFow29xT61shCCCHyl6z64I5cfpQxqb25rHrwsXYhfa3XUUq5ydDlg2hc0aPAn+FKshVCCJHnvF3sjQk0Pfn+fjiaz39XiFZLMEX7NS01f+GjjCFqrxdN6jxfoBOuxS8jCyGEEN4u9rSp4Y2VAr/r69Et5WNuqM4EWF2g2fauvPNFZIG+pCzJVgghRL7w8P3cA2pFXkn+nGN6XzyUeBZqP+fgyun8fem2pcN8KpJshRBC5BthtcuyfXhTPm5Thau40zllNH+k1cZWSWWCdjbHv+vL0t2nLR2mySTZCiGEyFcevqSchB0Dde/ypa4LelWhq2YLlde8yvWLJy0dpkkk2QohhMh30i8pWwEqVsxKC6WXbji3VCeqW53D+cfm3Dy60dJh5pgkWyGEEPlSWO2yrBjUgAddctmur07b5PEc0pcnVZdM2KKzBabRlCRbIYQQ+VZNHzcmPDQIxlXc6ZIymm4pH3NGX7rADO0oyVYIIUS+9nCjKYAUtBxRywOQpqqcj02yZHg5IslWCCFEvvdwo6mHaRQFP/cnz01uaZJshRBCFAiPjqusURTGd6xWIEaWkuEahRBCFBiPjqtcEBItSLIVQghRwDw8rnJBIZeRhRBCCDOTZCuEEEKYmSRbIYQQwswk2QohhBBmVuQaSKmqCkB8fPxjy+l0OpKSkoiPj0er1eZFaCIfkHovmqTei570Oof/8oI5Fblkm5CQAICPj4+FIxFCCJEfJCQk4OLiYtZjKGpepPR8RK/Xc/XqVYoVK4byoGN07dq12bdvX4Zy8fHx+Pj4cOnSJZydnS0Rarayijc/7NfU7XNa/knlnvZ1qffc2W9+rHdTX8uv9W6uOn/WfT/Ntvmt3tPr/NixY1SqVAkrK/PeVS1yZ7ZWVlaUKVMmwzqNRpPtF8zZ2Tlfffng8fFacr+mbp/T8k8q97SvS73nzn7zY70/7Wv5rd7NVefPuu+n2Ta/1nvp0qXNnmhBGkgBMGjQIEuHYBJzxfus+zV1+5yWf1K5p31d6j139psf6/1pX8tvzBnrs+z7abYt6vVe5C4j51R8fDwuLi7ExcXlq/90hXlJvRdNUu9FT17XuZzZZsPW1pbRo0dja2tr6VBEHpJ6L5qk3ouevK5zObMVQgghzEzObIUQQggzk2QrhBBCmJkkWyGEEMLMJNkKIYQQZibJVgghhDAzSba54Ny5czRt2pSAgACqV69OYmKipUMSecDPz48aNWoQGBhI06ZNLR2OyENJSUn4+voydOhQS4ci8sCdO3cIDg4mMDCQatWqMXv2bJP3UeSGazSHPn368Pnnn9OoUSNu3bolffWKkJ07d+Lk5GTpMEQeGzduHPXq1bN0GCKPFCtWjK1bt+Lg4EBiYiLVqlWjY8eOlChRIsf7kDPbZ/TPP/+g1Wpp1KgRAMWLF8faWv6HEaKwOn36NCdOnKBVq1aWDkXkEY1Gg4ODAwDJycmoqmrytHyFPtlu3bqVdu3aUapUKRRFYeXKlZnKzJo1Cz8/P+zs7Khbty579+7N8f5Pnz6Nk5MT7dq14/nnn2f8+PG5GL14WuaudwBFUQgJCaF27dosXLgwlyIXzyIv6n3o0KFERETkUsQiN+RFvd+5c4eaNWtSpkwZhg0bhru7u0nbF/pTsMTERGrWrMnrr79Ox44dM72+ZMkSwsPDiYyMpG7dukybNo2WLVty8uRJPD09AQgMDCQ1NTXTtuvXryc1NZVt27Zx6NAhPD09efnll6lduzbNmzc3+3sT2TN3vZcqVYrt27dTunRpoqOjadasGdWrV6dGjRpmf28ie+au93379lGxYkUqVqzIzp07zf5+RM7kxffd1dWVv//+m5iYGDp27Ejnzp3x8vLKeZBqEQKoK1asyLCuTp066qBBg4zP09LS1FKlSqkRERE52ufOnTvVFi1aGJ9PnDhRnThxYq7EK3KHOer9UUOHDlXnzp37DFGK3GaOeh8+fLhapkwZ1dfXVy1RooTq7Oysjh07NjfDFs8oL77vAwYMUJcuXWrSNoX+MvLjpKSksH//fpo1a2ZcZ2VlRbNmzdi1a1eO9lG7dm2uX7/O7du30ev1bN26lSpVqpgrZJELcqPeExMTSUhIAODu3bts3ryZqlWrmiVekTtyo94jIiK4dOkS58+fZ9KkSfTv359Ro0aZK2SRC3Kj3mNiYozf97i4OLZu3UqlSpVMiqPQX0Z+nNjYWNLS0jJdCvDy8uLEiRM52oe1tTXjx4+ncePGqKpKixYtaNu2rTnCFbkkN+o9JiaGDh06AJCWlkb//v2pXbt2rscqck9u1LsoeHKj3i9cuMCbb75pbBj1zjvvUL16dZPiKNLJNre0atVKWiYWMeXLl+fvv/+2dBjCgvr06WPpEEQeqVOnDocOHXqmfRTpy8ju7u5oNBpiYmIyrI+JiaFkyZIWikqYm9R70ST1XjTll3ov0snWxsaGoKAgNm3aZFyn1+vZtGkT9evXt2Bkwpyk3osmqfeiKb/Ue6G/jHz37l3OnDljfH7u3DkOHTpE8eLFKVu2LOHh4fTu3Zvg4GDq1KnDtGnTSExMpG/fvhaMWjwrqfeiSeq9aCoQ9f5U7Z4LkC1btqhApqV3797GMjNmzFDLli2r2tjYqHXq1FF3795tuYBFrpB6L5qk3oumglDviqqaOOaUEEIIIUxSpO/ZCiGEEHlBkq0QQghhZpJshRBCCDOTZCuEEEKYmSRbIYQQwswk2QohhBBmJslWCCGEMDNJtkIIIYSZSbIVIp+LiopCURTu3LmT58dWFAVFUXB1dX1suTFjxhAYGJjhefq206ZNM2uMQhQEkmyFyEeaNGnCe++9l2FdgwYNiI6OxsXFxSIxzZ07l1OnTpm0zdChQ4mOjqZMmTJmikqIgqXQT0QgREFnY2Nj0SngXF1d8fT0NGkbJycnnJyc0Gg0ZopKiIJFzmyFyCf69OnDn3/+yVdffWW8BHv+/PlMl5HnzZuHq6srq1evplKlSjg4ONC5c2eSkpKYP38+fn5+uLm5MWTIENLS0oz7T05OZujQoZQuXRpHR0fq1q1LVFTUU8U6YcIEvLy8KFasGP369eP+/fu58AkIUXhJshUin/jqq6+oX78+/fv3Jzo6mujoaHx8fLIsm5SUxPTp01m8eDFr164lKiqKDh06sGbNGtasWcNPP/3Et99+y7Jly4zbDB48mF27drF48WIOHz7Mq6++yssvv8zp06dNivOXX35hzJgxjB8/nr/++gtvb2++/vrrZ3rvQhR2chlZiHzCxcUFGxsbHBwcnnjZWKfT8c033+Dv7w9A586d+emnn4iJicHJyYmAgACaNm3Kli1bCAsL4+LFi8ydO5eLFy9SqlQpwHBfde3atcydO5fx48fnOM5p06bRr18/+vXrB8Dnn3/Oxo0b5exWiMeQM1shCiAHBwdjogXw8vLCz88PJyenDOuuX78OwJEjR0hLS6NixYrG+6lOTk78+eef/PvvvyYd+/jx49StWzfDuvr16z/DuxGi8JMzWyEKIK1Wm+G5oihZrtPr9QDcvXsXjUbD/v37MzVaejhBCyHMQ5KtEPmIjY1NhkZNuaVWrVqkpaVx/fp1GjVq9Ez7qlKlCnv27KFXr17Gdbt3737WEIUo1CTZCpGP+Pn5sWfPHs6fP4+TkxPFixfPlf1WrFiRHj160KtXLyZPnkytWrW4ceMGmzZtokaNGrRp0ybH+3r33Xfp06cPwcHBNGzYkIULF/LPP/9Qvnz5XIlViMJI7tkKkY8MHToUjUZDQEAAHh4eXLx4Mdf2PXfuXHr16sUHH3xApUqVCA0NZd++fZQtW9ak/YSFhfHJJ5/w4YcfEhQUxIULFxgwYECuxSlEYaSoqqpaOgghRP6kKAorVqwgNDT0qbb38/PjvffeyzQqlhBFjZzZCiEeq1u3biYPuzh+/HicnJxy9cxciIJMzmyFENk6c+YMABqNhnLlyuV4u1u3bnHr1i0APDw8LDausxD5hSRbIYQQwszkMrIQQghhZpJshRBCCDOTZCuEEEKYmSRbIYQQwswk2QohhBBmJslWCCGEMDNJtkIIIYSZSbIVQgghzEySrRBCCGFm/w9MZhWYJCVafwAAAABJRU5ErkJggg==",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"tm = np.logspace(np.log10(to[0]), np.log10(to[-1]), 100)\n",
"hm = ml.head(0, 0, tm)\n",
"plt.semilogx(to, ho, \".\", label=\"obs\")\n",
"plt.semilogx(tm, hm[0], label=\"timflow\")\n",
"plt.xlabel(\"time [d]\")\n",
"plt.ylabel(\"drawdown [m]\")\n",
"plt.title(\"Model Results - Single-layer model\")\n",
"plt.legend()\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Comparison of results\n",
"\n",
"We now compare the values in `timflow` and add the results of the modelling done in MLU (Hemker & Post, 2014). Results are similar between both models. The RMSE of MLU is slightly better than the one from `timflow`."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
" \n",
" \n",
" | | \n",
" k [m/d] | \n",
" Ss [1/m] | \n",
" RMSE [m] | \n",
"
\n",
" \n",
" \n",
" \n",
" | timflow | \n",
" 0.42 | \n",
" 1.70e-05 | \n",
" 0.004 | \n",
"
\n",
" \n",
" | MLU | \n",
" 0.41 | \n",
" 1.94e-05 | \n",
" 0.004 | \n",
"
\n",
" \n",
"
\n"
],
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = pd.DataFrame(\n",
" columns=[\"k [m/d]\", \"Ss [1/m]\", \"RMSE [m]\"],\n",
" index=[\"timflow\", \"MLU\"],\n",
")\n",
"\n",
"t.loc[\"timflow\"] = np.append(cal.parameters[\"optimal\"].values, cal.rmse())\n",
"t.loc[\"MLU\"] = [0.4133, 1.9388e-05, 0.004264]\n",
"\n",
"t_formatted = t.style.format(\n",
" {\"k [m/d]\": \"{:.2f}\", \"Ss [1/m]\": \"{:.2e}\", \"RMSE [m]\": \"{:.3f}\"}\n",
")\n",
"t_formatted"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### References\n",
"\n",
"* Cooper Jr, H.H., Bredehoeft, J.D. and Papadopulos, I.S. (1967) Response of a finite diameter well to an instantaneous charge of water, Water Resources Research 3, 263–269\n",
"* Duffield, G.M. (2007), AQTESOLV for Windows Version 4.5 User's Guide, HydroSOLVE, Inc., Reston, VA.\n",
"* Hemker, K. en Post V. (2014) MLU for Windows: well flow modeling in multilayer aquifer systems; MLU User's guide. https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fmicrofem.com%2Fdownload%2Fmlu-user.pdf&data=05%7C02%7CMark.Bakker%40tudelft.nl%7Cad7f16364d2d4fd55dbf08de73832eaa%7C096e524d692940308cd38ab42de0887b%7C0%7C0%7C639075204580287861%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=OBoe8seXZUfoat89Dfr4g6lF%2Bn1FdtXqtp%2F18BMXCn0%3D&reserved=0\n",
"* Hyder, Z., Butler Jr, J.J., McElwee, C.D. and Liu, W. (1994) Slug tests in partially penetrating wells, Water Resources Research 30, 2945–2957."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "timflow (3.13.11)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.11"
}
},
"nbformat": 4,
"nbformat_minor": 4
}