{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Synthetic Pumping Test - Calibration test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import ttim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use observation times from Oude Korendijk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "drawdown = np.loadtxt(\"data/oudekorendijk_h30.dat\")\n",
    "tobs = drawdown[:, 0] / 60 / 24\n",
    "robs = 30\n",
    "Q = 788"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self.neq  1\n",
      "solution complete\n"
     ]
    }
   ],
   "source": [
    "ml = ttim.ModelMaq(kaq=60, z=(-18, -25), Saq=1e-4, tmin=1e-5, tmax=1)\n",
    "w = ttim.Well(ml, xw=0, yw=0, rw=0.1, tsandQ=[(0, 788)], layers=0)\n",
    "ml.solve()\n",
    "rnd = np.random.default_rng(2)\n",
    "hobs = ml.head(robs, 0, tobs)[0] + 0.05 * rnd.random(len(tobs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### See if TTim can find aquifer parameters back"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..................................\n",
      "Fit succeeded.\n"
     ]
    }
   ],
   "source": [
    "cal = ttim.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq0\", layers=0, initial=100)\n",
    "cal.set_parameter(name=\"Saq0\", layers=0, initial=1e-3)\n",
    "cal.series(name=\"obs1\", x=robs, y=0, layer=0, t=tobs, h=hobs)\n",
    "cal.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>layers</th>\n",
       "      <th>optimal</th>\n",
       "      <th>std</th>\n",
       "      <th>perc_std</th>\n",
       "      <th>pmin</th>\n",
       "      <th>pmax</th>\n",
       "      <th>initial</th>\n",
       "      <th>inhoms</th>\n",
       "      <th>parray</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>kaq0</th>\n",
       "      <td>0</td>\n",
       "      <td>59.871288</td>\n",
       "      <td>0.670178</td>\n",
       "      <td>1.119365</td>\n",
       "      <td>-inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>100.000</td>\n",
       "      <td>None</td>\n",
       "      <td>[59.871287575268404]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.000121</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>3.633643</td>\n",
       "      <td>-inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.001</td>\n",
       "      <td>None</td>\n",
       "      <td>[0.00012118965722715004]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      layers    optimal       std  perc_std  pmin  pmax  initial inhoms  \\\n",
       "kaq0       0  59.871288  0.670178  1.119365  -inf   inf  100.000   None   \n",
       "Saq0       0   0.000121  0.000004  3.633643  -inf   inf    0.001   None   \n",
       "\n",
       "                        parray  \n",
       "kaq0      [59.871287575268404]  \n",
       "Saq0  [0.00012118965722715004]  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cal.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse: 0.014201293661013786\n"
     ]
    }
   ],
   "source": [
    "print(\"rmse:\", cal.rmse())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGhCAYAAACphlRxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+qklEQVR4nO3deXhTZfrG8W+SQlmkUagUMRDUsgqIIEsRBAoWRFwQF8SpI1MpiLjgCq44ozJuM8rg0kJcQFR0HBRRKygFVCibVFAR2p9UCFgWgRQECjT5/XFsobSUbslJ0vtzXblsTs5Jn/Zoc/ue9zyvxefz+RAREREJEVazCxARERGpCIUXERERCSkKLyIiIhJSFF5EREQkpCi8iIiISEhReBEREZGQovAiIiIiISXC7AKqm9frZdu2bTRo0ACLxWJ2OSIiIlIOPp+Pffv20bRpU6zWssdWwi68bNu2jWbNmpldhoiIiFTCli1bcDgcZe4TduGlQYMGgPHDR0VFmVyNiIiIlEdeXh7NmjUr+hwvS9iFl8JLRVFRUQovIiIiIaY8Uz40YVdERERCisKLiIiIhBSFFxEREQkpCi8iIiISUhReREREJKQovIiIiEhIUXgRERGRkKLwIiIiIiFF4UVERERCisKLiIiIhBSFFxEREQkpCi8V4Ha7SU9Px+12m11KSNLvT0REqoPCSzm5XC6cTifx8fE4nU5cLpfZJZVbMISGUP79iYhIcLH4fD6f2UVUp7y8POx2Ox6Pp9pWlXa73TidTrxeb9E2m81GTk4ODoejWr6Hv7hcLpKTk/F6vVitVlJTU0lKSgpoDaH8+xMRkcCoyOe3Rl7KISsrq9gHL0BBQQHZ2dkmVVQ+bre7KLgAeL1eRo8eHfARmFD9/YmISHAKSHh55ZVXOOecc6hTpw5dunTh66+/LnP/xYsX06VLF+rUqcO5557La6+9FogyT6ply5ZYrcV/VTabjdjYWJMqKp9gCQ2h+vsTEZHg5PfwMnv2bO6++24efvhh1qxZQ+/evbnsssvYvHlzqftv2rSJwYMH07t3b9asWcNDDz3EnXfeyYcffujvUk/K4XCQmpqKzWYDjA/elJSUoL/kESyhIVR/fyIiEpz8Puele/fudO7cmVdffbVoW9u2bbn66quZPHlyif0ffPBB5s6dy/r164u2jRkzhu+//55ly5aV2D8/P5/8/Pyi53l5eTRr1qxa57wAUFDA4c6dOVhQQGRUFHUaNIDatcv/qFWrYvuXdXytWmAtX+50uVyMHj2agoKCotAQ6DkvhdxuN9nZ2cTGxiq4iIhIMRWZ8xLhz0IOHz7M6tWrmTBhQrHtCQkJLF26tNRjli1bRkJCQrFtAwcOxOVyceTIEWrVqlXstcmTJ/PEE09Ub+GlOXKE2mvXUtv/36l8IiLKFXqSatfmposvZk/9+tTt3JnTzzoLcnKgefNyB6Dq4nA4FFpERKTK/Bpedu3aRUFBATExMcW2x8TEkJubW+oxubm5pe5/9OhRdu3axVlnnVXstYkTJ3LPPfcUPS8ceal2EREwbx4cPlzyceRI6dsr8yjtvfLz4cQBsqNHjceBA6csvQ5wFsDnnx/bWK8etGkD7dode7RtC+eea/ysIiIiQSogn1IWi6XYc5/PV2LbqfYvbTtAZGQkkZGR1VDlKUREwOWX+//7nExBQeVD0KFDsGkT/PQTrF8PGzYYoee774zH8WrXhtatj4WZwmDTsqXxmoiIiMn8Gl6io6Ox2WwlRll27NhRYnSlUJMmTUrdPyIigkaNGvmt1qBns0Hdusajqo4ehV9+ORZmfvrp2NcHD8K6dcbjxO8fG1typKZ1a2MUR0REJED8Gl5q165Nly5dWLBgAUOHDi3avmDBAq666qpSj4mLi+OTTz4ptm3+/PlcdNFFJea7SCVFRECrVsbj6quPbfd6YfPmY2Hm+HCTl2eM2GzYAHPmHDvGYoFzzik5UtO2LTRoEPAfTUREwp/f7zaaPXs2iYmJvPbaa8TFxZGamsq0adP48ccfcTqdTJw4ka1btzJjxgzAuFW6ffv2jB49mlGjRrFs2TLGjBnDu+++y7Bhw075/fzRYbfG8/lg27aSgebHH2H37pMf53CUHKlp1w4aNgxc7SIiEhKC5m4jgBtuuIHff/+dv//97/z222+0b9+ezz77DKfTCcBvv/1WrOfLOeecw2effcb48eN5+eWXadq0KVOmTClXcBE/sVjg7LONx6WXHtvu88HOncUvPRWGm99+A7fbeMyfX/z9evaE66+HYcOMgCMiIlIBWttI/GPPnmOh5vhwc2JzwsIgc+21RjgSEZEaqSKf3wovElhbt8KHH8L778O33xZ/7eKL4brrFGRERGoghReFl9CgICMiIn9SeFF4CT1utxFkPvhAQUZEpAZSeFF4CW2nCjLDh8PNN4POr4hI2FB4UXgJH4VB5v334bj1sLxRUVhvvx3uvBOaNDGxQBERqQ4V+fwO7Mp8IhXlcMBdd8G33/Lus88y3mLhZ8CalweTJ4PTCcnJRvO8CnC73aSnp+N2u/1Tt4iI+I3Ci4QEt9vNXyZM4EWfj3bAVcAyMNZumjbNaIB3zTWQkXHK93K5XDidTuLj43E6nbhcLj9XLyIi1UnhRUJCVlYWXq8XAB8wF+gJrJkyBa64wmiYN2cOxMXBJZcYK4D/uf/x3G43ycnJRe/l9XoZPXq0RmBEREKIwouEhJYtW2K1Fv/X1WazcebQoTB3rrFUwd/+BrVqwddfG4GmQwd4801jdOZPx4egQgUFBWRnZwfixxARkWqg8CIhweFwkJqais1mA4zgkpKSgqNweYF27cDlgk2b4IEHjDuRfvoJRo6Ec8+FF16AvLyThqDY2NhA/0giIlJJuttIQorb7SY7O5vY2NhjwaU0Hg+kpMCLLxrrLAHY7XDbbbwTHc3NDz5IQUFBUQhKSkoKSP0iIlI63Sqt8CKF8vNh1ix47jn4+WdjW+3a7B82jB8HDeLs+PiyQ5CIiASEbpUWKRQZacyF+fFHY27MxRfD4cOc9u67dP/b33A88wzs3m12lSIiUgEKL1IzWK3GJN5vvjG69g4ZAgUFMHUqtGwJr7wCR4+aXaWIiJSDwovUPD17wiefwJdfQvv2xsjL7bdD586Qnl7ut1GjOxERcyi8SM3Vvz+sWWOMvpxxBqxbB/HxMGyYcddSGdToTkTEPJqwKwLw++8waRK8+qpxOSkyEu6/HyZMgPr1i+3qdrtxOp3F+sXYbDZycnI0+VdEpJI0YVekoho1gv/8BzIzjdGX/Hx48klo3Rreecfo4PsnNboTETGXwovI8dq3N+bC/O9/cM45sHUr3HQT9OoFq1YBJ+/2q0Z3IiKBofAiciKLBYYONTr0PvWUcdlo6VLo1g2SknBERJTd7VdERPxKc15ETmXrVpg4EWbONJ43aACPPYb7mmvI3rz51N1+RUTklDTnRaQ6nX02zJhhjL507Qr79sH99+MYNIi+hw4puIiIBJjCi0h5xcVBRga88QbExEBWFlx2GSQlGWspiYhIQCi8iFSE1Qq33AIbN8LddxvzY15/HTp0MCb6ioiI3ym8iFRGVBT8+9+weDGcdx5s2QKXXgpjx8L+/WZXJyIS1hReRKqid2/4/ntjeQEwmtxdcAEsWWJuXSIiYUzhRaSq6tc3lhj48kto3hx++QX69oXx4+HAAbOrExEJOwovItWlf39jfaRbbzU68r74Ilx4ISxbZnZlIiJhReFFpDpFRcG0afDZZ9C0qTGxt1cvePBBOHTI7OpERMKCwouIP1x2GfzwA9x8M3i98Oyz0KULrF5tdmUiIiFP4UXEX844A956Cz76yOgL89NP0L07PPYYHD5sdnUiIiFL4UXE3666yhiFueEGKCiAf/zDCDFr15pdmYhISFJ4EQmE6Gh47z14/31o1AgyM+Gii+Cll4zJvSIiUm4KLyKBdN118OOPxmjMkSNGl96rr4bdu82uTEQkZCi8iARaTAzu//yHjXfeia92bZg7Fzp1gm+/NbsyEZGQoPAiEmAulwtnixa0njKFi44cwRMTYywv0KcPTJ5s3J0kIiInpfAiEkBut5vk5GS8fwaU73w+nDt38sfQocZk3ocegkGDYPt2kysVEQleCi8iAZSVlVUUXAp5vF5W3nGHsTp1vXqwYIGxPpJWqRYRKZXCi0gAtWzZEqu1+H92NpuN2JYtYeRIWLkS2rc3Rl4SEuCRR+DoUZOqFREJTgovIgHkcDhITU3FZrMBRnBJSUnB4XAYO7RrBytWQHKycQv1U09Bv36wdauJVYuIBBeLzxdeTSby8vKw2+14PB6ioqLMLkekVG63m+zsbGJjY48FlxPNng2jRsG+fdC4sfG8b9+A1ikiEigV+fzWyIuICRwOB3379j15cAGjI++aNcb8lx07YMAAeP55NbUTkRpP4UUkmJ13HixdComJxt1I998P119vjMaIiNRQCi8iwa5ePWOBx5dfhlq14L//hW7dYP16sysTETGFwotIEHO73aSnp+PeuhXGjoUlS+Dss+Hnn40A88EHZpcoIhJwCi8iQcrlcuF0OomPj8fpdOJyuaBHD/juO2Pi7v79xiWk++7T7dQiUqPobiORIOR2u3E6ncUa2tlsNnJycoxJvkePGt14n3vOeLFPH+NupJgYkyoWEaka3W0kEuJK68RbUFBAdna28SQiAp591pj/ctppsHgxdO1qjMqIiIQ5hReRIHTSTryxscV3HDYMVq7kyHnnwZYteC++GN57L4CViogEnsKLSBA6ZSfe47i+/ZaYX37hc8B66BDceCNMnGjcWn2cosm/bncgfgQREb/RnBeRIHaqTrzHz42xApOBBwpfvPxyeOcdiIrC5XIVrWZttVpJTU0lKSkpgD+JiEjZKvL5rfAiEsLS09OJj48vtm0EMKN2bWyHD0PbtuSmpHB2374nn/wrIhIENGFXpIYobW7MbJuNXR9+CE2bwvr1NBo8mP5lTf4VEQkxCi8iIexkc2NihgyBVaugRw9q7d/P58D4444rdfKviEiI0GUjkTBw0rkx+flw223wxhsAvAGMtVqZqjkvIhJkNOdF4UXkGJ8PpkzBd889WLxe8rt3J3LePIiONrsyEZEiQTPnZc+ePSQmJmK327Hb7SQmJrJ3796T7n/kyBEefPBBOnToQP369WnatCk333wz27Zt82eZIuHNYoG77sLy6afQoAGRy5cbywz8/LPZlYmIVIpfw8uIESPIzMwkLS2NtLQ0MjMzSUxMPOn+Bw4c4LvvvuPRRx/lu+++43//+x8bN27kyiuv9GeZIjXDoEGwdCm0aAH/938QFwdffWV2VSIiFea3y0br16+nXbt2ZGRk0L17dwAyMjKIi4vj559/pnXr1uV6n5UrV9KtWzd+/fVXmjdvfsr9ddlI5BR27IChQ40gExEBr7wCo0bhdrvJysqiZcuWuoVaRAIuKC4bLVu2DLvdXhRcAHr06IHdbmfp0qXlfh+Px4PFYuH0008v9fX8/Hzy8vKKPUSkDI0bGyMuN91kLPCYnMy6hATOad68+ArWIiJBym/hJTc3l8aNG5fY3rhxY3Jzc8v1HocOHWLChAmMGDHipCls8uTJRXNq7HY7zZo1q1LdIjVCnTowcyb8/e8AdFiwgA99PuoDXq+X0aNHaxkBEQlaFQ4vkyZNwmKxlPlYtWoVABaLpcTxPp+v1O0nOnLkCMOHD8fr9fLKK6+cdL+JEyfi8XiKHlu2bKnojyRSM1ks8Oij/PjooxwErgS+BRyoiZ2IBLeIih4wbtw4hg8fXuY+LVq0YO3atWzfvr3Eazt37iQmJqbM448cOcL111/Ppk2bWLhwYZnXviIjI4mMjCxf8SJSgj05mfgnn2SOz8cFwDLgSqtVTexEJGhVOLxER0cTXY7+EHFxcXg8HlasWEG3bt0AWL58OR6Ph549e570uMLgkpWVRXp6Oo0aNapoiSJSAQ6Hg1unTaNncjKfeL2cD2TUrk3t9etBE3dFJAj5bc5L27ZtGTRoEKNGjSIjI4OMjAxGjRrFkCFDit1p1KZNG+bMmQPA0aNHufbaa1m1ahWzZs2ioKCA3NxccnNzOXz4sL9KFanxkpKSWPLrr+yZO5dDcXHUPnQIBg+GN980uzQRkRL82udl1qxZdOjQgYSEBBISEujYsSMzZ84sts+GDRvweDyA0eJ87ty5uN1uOnXqxFlnnVX0qMgdSiJScQ6Hg15XXEGd9HQYMcK4E2nkSHjiCaNLr4hIkNDyACJSktcLjzwCkycbz0eOhJQUqFXL3LpEJGwFRZ8XEQlhVis8/TS89prx9RtvwOWXg/ooiUgQUHgRkZMbPRrmzoV69WDBArjkEti61eyqRKSGU3gRkbJdfjksXgwxMfD998aijuvWmV2ViNRgCi8icmoXXQTLlkGbNuB2Q69esHCh2VWJSA2l8CIi5XPOOfDtt9C7tzH3ZdAgY4kBEZEAU3gRkfJr2BDmz4cbboAjR+Dmm+Gpp3QrtYgElMKLiFRMnTrwzjvwwAPG80cegeRkI8yIiASAwouIVJzVCs88A1OnGl9Pnw5XXQV//GF2ZSJSAyi8iEjl3X47zJkDdevC559D//6wa5fZVYlImFN4EZGqufJK+OorYz7M8uXGnUi//mp2VSISxhReRKTq4uLgm2+gWTPYsAEuvhh++MHsqkQkTCm8iEiVud1u0nNz+e3DD6FdO6MLb+/eRqAREalmCi8iUiUulwun00l8fDyOHj2YmZwMPXvC3r1w6aXG8gIiItVIq0qLSKW53W6cTider7dom81m49f16zn7nntg3jzjbqTUVEhKMrFSEQl2WlVaRAIiKyurWHABKCgoIGvrVuMupL/9DbxeuPVWY5Xq8Pp/JRExicKLiFRay5YtsVqL/xmx2WzExsZCRITR/2XiROOFhx+Gu+4ywoyISBUovIhIpTkcDlJTU7HZbIARXFJSUnA4HMYOFosx4vLSS8bz//wHRoyA/HyTKhaRcKA5LyJSZW63m+zsbGJjY48FlxO9956xFtKRIxzq1YsVEyZw7gUXnHx/EalRNOdFRALK4XDQt2/fsoPI8OHw2WccjoykzjffUH/IELo0b47L5QpcoSISFhReRCRg3G3acPHhw+wAugBLfD6eTE7G7XabXZqIhBCFFxEJmKysLFb5fPQCfgVaA0u8XrYtXGhyZSISShReRCRgCu9OygJ6AeuBZkCXu++GVatMrU1EQofCi4gEzPF3J7mBvlYrO1u0wLZnD/TrB+npZpcoIiFA4UVEAiopKYmcnBzS09NZ/euvnLl2LcTHw/79MGgQfPRR0b5ut5v09HTNiRGRYnSrtIiY79AhuPFGI7hYrTB9Oi6vl+TkZLxeL1arldTUVJK0xIBI2KrI57fCi4gEh6NHITkZ3ngDgHstFv513J8nm81GTk6O+sKIhCn1eRGR0BMRAS4X3HsvAC/4fDx13MsFBQVkZ2ebU5uIBBWFFxEJHhYLPPccngkTAHgIeBXjD1XRmkkiUuNFmF2AiEgxFgv2yZP5ZutWes6cyRjgdIuFAy+/rEtGIgJo5EVEglSvGTPY8/LLeCMiGO7zcdUbb+DOyjK7LBEJAgovIhK0Go0dy/xx4zgANFq+nOxWrZgxdarZZYmIyRReRCRoud1uLp8yhQTAA/QF2t5xB9vWrTO3MBExlcKLiAStrKwsvF4v3wL9gJ1AV8B+5ZXw22+AGtmJ1EQKLyIStArXQgJYA1wCbAXq5+RAr17M/uc/cTqdxMfH43Q6cblcJlYrIoGi8CIiQev4tZAAsmw2vpk8Gc49F375hYsnTqSV1wuA1+tl9OjRGoERqQEUXkQkqB2/FlJOTg43TJgAX3/NH04nDmAJcOGf+6qRnUjNoD4vIhL0HA5H8R4vTZvimTuX9RdcwEVAOjAYWK5GdiI1gkZeRCQkNe3YkZ/+8x+WAHZgPjD37rvVyE6kBlB4EZGQdfO4cZy3cSO7L7qI+sDgqVNh3jyzyxIRP1N4EZGQdnbLljT85hu4+mrIz4ehQ+H9980uS0T8SOFFREJfZKQRWEaMgKNH4cYb4c03yzxE/WFEQpfCi4iEh1q1YMYMuPVW8Hph5Eh45ZVSd3W5XOoPIxLCLD6fz2d2EdUpLy8Pu92Ox+MhKirK7HJEJNB8Phg/Hl56yXj+zDPwwANFL7vdbpxOJ94/+8MA2Gw2cnJyNNlXxEQV+fzWyIuIhBeLBf79b3j4YeP5gw/CY48ZoYZjSw4cT/1hREKLwouIhB+LBZ58Ep5+2nj+j3/AffeBz1dsyYFCNvWHEQkpCi8iEr4mTjx2+ehf/4LbbsPRtGmxJQdsNhspKSm6ZCQSQtRhV0TC2513wmmnGRN5U1Lgjz9IeuMNBg4cSHZ2NrGxsQouIiFG4UVEwt/f/gb16sFf/gJvvw0HDuB4912FFpEQpctGIlIzDB8OH34ItWvD//5nNLU7eNDsqkSkEhReRKTmuOoq+OQTqFsXPv8cBg+G/fvNrkpEKkjhRURqloQE+OILaNAAFi2CQYPA4zG7KhGpAIUXEal5eveGBQvAbodvv4VLL4U9e8yuSkTKSeFFRGqm7t1h4UJo1AhWroT4eNi1y+yqRKQcFF5EpObq3BnS06FxY8jMhL59ITfX7KpE5BQUXkSkZuvQARYvhqZN4ccfoU8f2LrV7KpEpAwKLyIibdrAkiXQvDls3AiXXAK//mp2VSJyEgovIiIA551nBJhzz4VffjECjBZrFAlKfg0ve/bsITExEbvdjt1uJzExkb1795b7+NGjR2OxWHjxxRf9VqOISBGn0wgwrVrB5s3GJaSffza7KhE5gV/Dy4gRI8jMzCQtLY20tDQyMzNJTEws17EfffQRy5cvp2nTpv4sUUSkuLPPNubAnH8+bNtmBJh168yuSkSO47fwsn79etLS0pg+fTpxcXHExcUxbdo05s2bx4YNG8o8duvWrYwbN45Zs2ZRq1atMvfNz88nLy+v2ENEpEqaNDEa2HXqBDt2GHchffedyUWJSCG/hZdly5Zht9vp3r170bYePXpgt9tZunTpSY/zer0kJiZy//33c/7555/y+0yePLnospTdbqdZs2bVUr+I1HDR0UYfmK5dYfduow/M8uUAuN1u0tPTcbvdJhcpUjP5Lbzk5ubSuHHjEtsbN25Mbhl9FJ555hkiIiK48847y/V9Jk6ciMfjKXps2bKl0jWLiBRzxhnw5Zdw8cXGEgIDBjDvwQdxOp3Ex8fjdDpxuVxmVylS41Q4vEyaNAmLxVLmY9WqVQBYLJYSx/t8vlK3A6xevZqXXnqJN99886T7nCgyMpKoqKhiDxGRahMVBWlp0K8f7N9Pv2efpY/XCxgjxaNHj9YIjEiARVT0gHHjxjF8+PAy92nRogVr165l+/btJV7buXMnMTExpR739ddfs2PHDpo3b160raCggHvvvZcXX3yRnJycipYrIlJ1p50Gn37K73360GjlSj4FhgJfYPyNys7OxuFwmFykSM1R4fASHR1NdHT0KfeLi4vD4/GwYsUKunXrBsDy5cvxeDz07Nmz1GMSExMZMGBAsW0DBw4kMTGRkSNHVrRUEZHqU7cuB999l09iY7kCmAtcC3xmsxEbG2tycSI1i9/mvLRt25ZBgwYxatQoMjIyyMjIYNSoUQwZMoTWrVsX7demTRvmzJkDQKNGjWjfvn2xR61atWjSpEmxY0REzOA47zx2vfYa/wVqAx8CacnJGnURCTC/9nmZNWsWHTp0ICEhgYSEBDp27MjMmTOL7bNhwwY8Ho8/yxARqTYjR4+mx6ZNbI+PpxYwIDUV3n/f7LJEahSLz+fzmV1EdcrLy8Nut+PxeDR5V0T8p6AARo6EmTPBajX+OWKE2VWJhKyKfH5rbSMRkcqw2eCNN4wA4/VCYiLMmGF2VSI1gsKLiEhl2WwwfTokJxsB5pZb4PXXza5KJOwpvIiIVIXVCq++CmPHgs8HSUmQkmJ2VSJhTeFFRKSqrFaYOhXuust4PmYMvPyyuTWJhDGFFxGR6mCxwL//DffdZzwfNw5efNHUkkTClcKLiEh1sVjg2Wdh4kTj+fjx8Nxz5tYkEoYUXkREqpPFAk89BY89Zjx/4AF4+mlzaxIJMwovIiLVzWKBJ56Av//deP7ww8e+FpEqU3gREfGXRx+FyZONrx9/3HgeXn1BRUyh8CIi4k8TJsDzzxtfP/mkMR9GAUakShReRET87d57j9159Mwzxh1JCjAilabwIiISCHfddaz3y7/+BXffrQAjUkkKLyIigTJ27LHuu1OmwO23G8sKiEiFKLyIiARScrKx/pHFYiwrMHq0AoxIBSm8iIgE2siR8NZbxrIC06cb6yEVFJhdlUjIUHgRETFDYiK8/baxMvWbbxorUh89anZVIiFB4UVExCw33gjvvgsREUaQSUxUgBEpB4UXEREzXXcdvP8+vlq14L33ODB0KBw5YnZVIkFN4UVExGSu3bu5+uhRDgP15s1jU1ycAoxIGRReRERM5Ha7SU5OZq7Px1AgHzhn9WoOXnklHD5sdnkiQUnhRUTERFlZWXj/vFX6M+Bq4BBQNy0Nrr9eAUakFAovIiImatmyJVbrsT/FacA1Viu+yEj4+GMYNgzy880rUCQIKbyIiJjI4XCQmpqKzWYDwGazMSw1Fcu8eVCnDsybB9dcA4cOmVypSPCw+HzhtbhGXl4edrsdj8dDVFSU2eWIiJSL2+0mOzub2NhYHA6HsXHhQhgyBA4ehEGDYM4cI9CIhKGKfH4rvIiIBLNFi+Dyy+HAAbj0UuNSUt26ZlclUu0q8vmty0YiIsGsb1/4/HOoXx8WLIArrjCCjEgNpvAiIhLsLrnECDCnnQZffWVcSvrjD7OrEjGNwouISCjo3RvS0qBBA0hPV4CRGk3hRUQkVFx8MXzxhRFgFi2CwYNh/36zqxIJOIUXEZFQEhcH8+dDVBQsWaIAIzWSwouISAhxu92kHzzI9rffBrsdvv4aLrsM9u0zuzSRgFF4EREJES6XC6fTSXx8PE2vvpqPx40zAsw33xh9YPLyzC5RJCAUXkREQkDhAo6F6yB5vV6G/fOfbJ81C04/HZYuVYCRGkPhRUQkBBy/gGOhgoIC1tevD19+CWecAcuWwcCB4PGYVKVIYCi8iIiEgBMXcARjHaTY2Fjo0uVYgMnIgIQE2LvXnEJFAkDhRUQkBJS2gGNKSsqxdZA6dzYa2DVsCCtWKMBIWNPaRiIiIaTUBRyP9/330L8//P47XHSRcVv1GWcEvlCRCtLCjAovIlKTrV1rBJhdu4xLSgsWKMBI0NPCjCIiNVnHjrBwIQUNG8Lq1Rzu0wf27DG7KpFqo/AiIhKGXCtW0HnPHnYCtdetY1enTrB7t9lliVQLhRcRkTBT2BNmrc9HP2AHEL15szECowAjYUDhRUQkzBzfE+ZHoB+wHaj9ww8wYIACjIQ8hRcRkTBzYk+Yn4ABVisF0dGwZs2xu5FEQpTCi4hImCmtJ8zdqanYFi+GmBjIzIQBA9i2bh3p6em43W5zCxapIN0qLSISpkrtCbN+PfTrB9u3kwkMAPZYraSmppKUlGRitVLTqc+LwouIyEnlpqdDfDxNoCjA7LXZyMnJKb3xnUgAqM+LiIic1HqgL/Ab0An4Cji9oIDs7GwTqxIpP4UXEZEapmXLlmRZrfTDCDAXYASYVg0blrq/2+3W3BgJKgovIiI1TOGE3mybrViAaZqYaCwpcByXy4XT6SQ+Ph6n04nL5TKjZJFiNOdFRKSGKpzQ28Ziocnw4ZCbCx06wMKFEB2N2+3G6XQW9YwB484lzY0Rf6jI53dEgGoSEZEg43A4joWQRYugb19Ytw7i42HhwmLN7goV/Dk3RuFFzKTLRiIiAq1bGwHmrLOKAkzrhg2LNbsDY+QlNjbWnBpF/qTwIiIihtatIT29KMA0TUxkxgsvFGt2l5KSolEXMZ0uG4mIyDGFAaZfP1i3jptef51+q1ezcc+e4s3uREykkRcRESnuxBGYv/yFvuefr+AiQUPhRURESjp+DswPPxiTeHfuNLsqEcDP4WXPnj0kJiZit9ux2+0kJiayd+/eUx63fv16rrzySux2Ow0aNKBHjx5s3rzZn6WKiMiJWrUyAkzTpgowElT8Gl5GjBhBZmYmaWlppKWlkZmZSWJiYpnH/N///R+9evWiTZs2LFq0iO+//55HH32UOnXq+LNUEREpTatWxy4hKcBIkPBbk7r169fTrl07MjIy6N69OwAZGRnExcXx888/07p161KPGz58OLVq1WLmzJmV+r5qUici4gcbNxp9YH77Ddq3NxrZnXmm2VVJGAmKhRmXLVuG3W4vCi4APXr0wG63s3Tp0lKP8Xq9fPrpp7Rq1YqBAwfSuHFjunfvzkcffXTS75Ofn09eXl6xh4iIVLPCS0iFIzD9+2sERkzjt/CSm5tL48aNS2xv3Lgxubm5pR6zY8cO9u/fzz//+U8GDRrE/PnzGTp0KNdccw2LFy8u9ZjJkycXzamx2+00a9asWn8OERH50/EBZt06BRgxTYXDy6RJk7BYLGU+Vq1aBYDFYilxvM/nK3U7UNSG+qqrrmL8+PF06tSJCRMmMGTIEF577bVSj5k4cSIej6fosWXLlor+SCIiUl4KMBIEKtykbty4cQwfPrzMfVq0aMHatWvZvn17idd27txJTExMqcdFR0cTERFBu3btim1v27Yt33zzTanHREZGEhkZWc7qRUSkygon8f7ZyI7+/eGrrzQHRgKmwuElOjqa6OjoU+4XFxeHx+NhxYoVdOvWDYDly5fj8Xjo2bNnqcfUrl2brl27smHDhmLbN27ciNPprGipIiLiLyd04qV//6LVqEX8zW9zXtq2bcugQYMYNWoUGRkZZGRkMGrUKIYMGVLsTqM2bdowZ86couf3338/s2fPZtq0aWRnZzN16lQ++eQTxo4d669SRUSkMk7oxEt8POzaVWwXt9tNeno6brfbpCIlHPm1z8usWbPo0KEDCQkJJCQk0LFjxxK3QG/YsAGPx1P0fOjQobz22ms8++yzdOjQgenTp/Phhx/Sq1cvf5YqIiKVUUaAcblcOJ1O4uPjcTqduFwuk4uVcOG3Pi9mUZ8XERETbNhgXEL67Tfo0IFtb79NswsvLLoRA4xVqXNycrRGkpQqKPq8iIhIDVI4AtOkCaxbR9TQoZxxXHABKCgoIDs726QCJZwovIiISPUoXMyxSRNO++UXvgIaHfeyzWYjNjbWpOIknCi8iIhI9TkuwFwARQHGZrORkpKiS0ZSLSp8q7SIiEiZCgNM375ckJtLzrnnkjdnDk07djS7MgkTGnkREZHqd8IlpKY33wy//252VRImFF5ERMQ/jp/E+/33RiM7BRipBgovIiLiP23aFA8wAwYowEiVKbyIiIh/FQaYmBjIzFSAkSpTeBEREf87McBceins3m12VRKiFF5ERCQw2rY1Fm9s3BjWrDFGYBRgpBIUXkREJHDatTNGYAoDjEZgpBIUXkREJLCODzDffWcEmD17zK5KQojCi4iIBF67dsYlpDPPVICRClN4ERERc5x/vjECc+aZsHo1JCTA3r1mVyUhQOFFRETMc/755L7zDodPPx1WrTJGYBRg5BQUXkRExDQul4uzBw6ky9697AQjwGgERk5B4UVEREzhdrtJTk7G6/XyA9Af2AWwciUMHAgej7kFStBSeBEREVNkZWXh9XqLnq8D4oEjUVGwYoUCjJyUwouIiJiiZcuWWK3FP4Z+stnY/cEH0KgRLF8OgwZBXp5JFUqwUngRERFTOBwOUlNTsdlsANhsNlJSUohJSIAvv4SGDSEjwxiBUYCR41h8Pp/P7CKqU15eHna7HY/HQ1RUlNnliIjIKbjdbrKzs4mNjcXhcBx7Yc0a6N/f6P8SFwdpaaC/62GrIp/fCi8iIhK8jg8wPXsaAaZBA7OrEj+oyOe3LhuJiEjwuvBC4xLSGWfA0qVw2WWwb99Jd3e73aSnp+N2uwNYpASawouIiAS3zp1hwQI4/XT49lsYPLjUAONyuXA6ncTHx+N0OnG5XIGvVQJCl41ERCQ0HN+Bt1cv+PxzOO00wBhxcTqdxW69ttls5OTkFJ9HI0FLl41ERCT8XHSRMQJjt8M33xgjMPv3AyV7xgAUFBSQnZ1tRqXiZwovIiISOgoDTFQUfP01XH45/PFHqT1jbDYbsbGxJhUq/qTwIiIioaVr12MBZskSuPxyHGecUWrPGF0yCk+a8yIiIqFp+XJjEce8POjbF+bNw71nT+k9YyToac6LiIiEv+7d4YsvjL4vixbBFVfgaNiQvn37KriEOYUXEREJXT16HAsw6ekwZAgcOGB2VeJnCi8iIhLaCpcOOO00I8BccYUCTJhTeBERkdDXs6cxAnPaabBwIVx1FRw8aHZV4icKLyIiEh4K1z6qX99YUkABJmwpvIiISPi4+GKj8279+sbt1FdfDYcOmV2VVDOFFxERCS+9e8NnnxkBZv58BZgwpPAiIiLh55JL4NNPoV49Yy7M0KEKMGFE4UVERMJTnz7HAkxaGgwbBvn5Zlcl1UDhRUREwtefnXepW9e4lKQAExYUXkREJLz162cEmDp1jJGYa69VgAlxCi8iIhL+4uOPBZh58+C66+DwYbOrkkpSeBERkZqhf3/45BMjwHzyCVx/vQJMiFJ4ERGRmmPAAPj4Y4iMNP55ww0KMCFI4UVERGqWhIRjAeajjzh41VUsWrAAt9ttdmVSTgovIiJS8wwcCB99REFEBHXT0vg9IYHzmjfH5XKZXZmUg8KLiIjUSO727bmyoIB8YBjwts/H7cnJGoEJAQovIiJSI2VlZfGZz8c1wGHgOuAtr5f/27DB5MrkVBReRESkRmrZsiVWq5XPMEZeDgM3AF2nTIGjR80tTsqk8CIiIjWSw+EgNTUVm83GPOB6q5UCm416c+fCzTcrwAQxhRcREamxkpKSyMnJIT09nam//ortww8hIgLefRduuQUKCswuUUoRYXYBIiIiZnI4HDgcjsIn8P77RgO7WbPAaoU33gCbzdwipRiNvIiIiBxv6FCYPdsILDNnwt/+phGYIKPwIiIicqJrroH33jMCzIwZkJSkABNEFF5ERERKc+21xtwXmw3eegtGjQKv1+yqBIUXERGRk7vuuuJzX5KTFWCCgMKLiIhIWW64Ad5+2wgwLheMHl0UYNxuN+np6erKG2C620hERORUbrwRfD5ITITp08FqxXXRRSSPGYPX68VqtZKamkpSUpLZldYIfh152bNnD4mJidjtdux2O4mJiezdu7fMY/bv38+4ceNwOBzUrVuXtm3b8uqrr/qzTBERkVMbMcKY+2KxQGoqR5KT8f45AuP1ehk9erRGYALEr+FlxIgRZGZmkpaWRlpaGpmZmSQmJpZ5zPjx40lLS+Ptt99m/fr1jB8/njvuuIOPP/7Yn6WKiIic2l/+Am++ic9iYQww9biXCgoKyM7ONquyGsVv4WX9+vWkpaUxffp04uLiiIuLY9q0acybN48NZSx6tWzZMv7617/St29fWrRoQXJyMhdccAGrVq0qdf/8/Hzy8vKKPURERPzm5pvZ88ILeIHbgSl/brbZbMTGxppYWM3ht/CybNky7HY73bt3L9rWo0cP7HY7S5cuPelxvXr1Yu7cuWzduhWfz0d6ejobN25k4MCBpe4/efLkostSdrudZs2aVfvPIiIicryG48fzzciReIE7gJcsFlJee+1Yp17xK7+Fl9zcXBo3blxie+PGjcnNzT3pcVOmTKFdu3Y4HA5q167NoEGDeOWVV+jVq1ep+0+cOBGPx1P02LJlS7X9DCIiIidzyeuvs/e55wC40+cj6YcfjEm94ncVDi+TJk3CYrGU+Si8xGOxWEoc7/P5St1eaMqUKWRkZDB37lxWr17NCy+8wNixY/nyyy9L3T8yMpKoqKhiDxERkUBoeN99MG2a8eSll+DeexVgAqDCt0qPGzeO4cOHl7lPixYtWLt2Ldu3by/x2s6dO4mJiSn1uIMHD/LQQw8xZ84cLr/8cgA6duxIZmYmzz//PAMGDKhouSIiIv51661G35fRo+Hf/zY68j77rHFXkvhFhcNLdHQ00dHRp9wvLi4Oj8fDihUr6NatGwDLly/H4/HQs2fPUo85cuQIR44cwWotPiBks9mKbkcTEREJOoWdd2+7DZ5/3mho989/KsD4id/mvLRt25ZBgwYxatQoMjIyyMjIYNSoUQwZMoTWrVsX7demTRvmzJkDQFRUFH369OH+++9n0aJFbNq0iTfffJMZM2YwdOhQf5UqIiJSdWPGwMsvG18/+yw89JAuIfmJXzvszpo1izvvvJOEhAQArrzySqZOnVpsnw0bNuDxeIqev/fee0ycOJGbbrqJ3bt343Q6eeqppxgzZow/SxUREam6sWONEZg77jBGXqxWePJJjcBUM4vPF16xMC8vD7vdjsfj0eRdERExx5QpcNddxtePPAJ//7sCzClU5PNbCzOKiIhUtzvvNCbvgjHy8sQT5tYTZhReRERE/OHuu+GFF4yvn3hCAaYaKbyIiIj4yz33wJ+N7Jg0Cf7xD1PLCRcKLyIiIv50333wzDPG1489Bk89ZW49YUDhRURExN8eeAAmTza+fuSRY19LpSi8iIiIBMKECcdGXR56yOgFI5Wi8CIiIhIoDz10bN7Lgw8a3XilwvzapE5ERETA7XaTlZVFy5YtcTzyiNHI7vHH4f77jUZ299xjdokhRSMvIiIifuRyuXA6ncTHx+N0OnG5XMbE3ccfN3a491548UVTaww16rArIiLiJ263G6fTWWxxYZvNRk5ODo6zzzYCTOFlpJdeMprb1VDqsCsiIhIEsrKyigUXgIKCArKzs43lAp54Ah5+2HjhrrvghPX/pHQKLyIiIn7SsmVLrNbiH7U2m43Y2FjjicVijLxMmGA8v+MOeOWVAFcZehReRERE/MThcJCamorNZgOM4JKSkoLD4Ti2k8UCTz9t9IIBuP12eO01E6oNHZrzIiIi4mdut5vs7GxiY2OLB5fj+XxGgCm8fTolBZKTA1ekySry+a1bpUVERPzM4XCcPLQUsliMxnVeL/zrXzB6NNhskJQUmCJDiC4biYiIBAuLxRh5uesu4/moUfD66+bWFIQ08iIiIhIEijWy+/e/jRGY//wHbr3VGIH561/NLjFoaORFRETEZCUa2b3+utH35fbbjbkwI0fCzJlmlxk0NGFXRETERKdsZHf77fDqq8YlpRkz4C9/MbFa/1GTOhERkRBxykZ2U6cak3d9PuPS0TvvmFRp8FB4ERERMdEpG9lZrUbjulGjjHkwiYkwe7YJlQYPhRcRERETlauRndVqNK5LSjICzE03wfvvm1Sx+TTnRUREJAiUq5Gd12vcffTGG8YdSO+9B9deG9hC/URN6kREREJMuRrZWa0wfboRYt56C2680dh2zTWBKTJI6LKRiIhIKLFaweUy5r4cPQo33AAffWR2VQGl8CIiIhJqbDbj0tFNNxkB5rrrYO5cs6sKGIUXERGRUGSzGZeORowwAsy118K8eWZXFRAKLyIiIqGqMMDccAMcOQLDhsGnn5pdld8pvIiIiISyiAh4+23j0tHhw8bk3c8/N7sqv1J4ERERCXHu3FwWjRrFgcGDjQAzdCh88UXp+7rdpKen43a7A1xl9VF4ERERCWGFizr2S0jgjM8/J+fCCyE/H666CubPL3XfogUgXS6Tqq4aNakTEREJUaUt6ljHamX3gAHUnT8f6tSBTz6BAQPKXgDyVP1lAkALM4qIiNQApS3qeMjrZcV998EVV8ChQ8Y/Fy4sewHIEKPwIiIiEqJOtqjjeW3bwgcfwJAhRoAZMoT2u3aVvQBkCFF4ERERCVFlLuoYGQn//S8MHgwHD3LmLbcw9957y14AMkRozouIiEiIK3NRx0OHjLuP0tKgXj12zJjBT40alb0ApAkq8vmt8CIiIhLuDh06dvdR/fpGkOnVy+yqitGEXRERETmmTh1j8cYBA+CPP+Cyy2DpUrOrqjSFFxERkZqgbl34+GPo3x/274dBgyAjw+yqKkXhRUREpKaoV89YfbpfP9i3DwYOhOXLza6qwhReREREapJ69YzGdX36QF4eJCTAypVmV1UhCi8iIiI1Tf36MG8e9O5tBJhLL4VVq8yuqtwUXkRERGqi006Dzz6Diy8Gj8cIMN99BwT/4o0KLyIiIjXVaafB559Dz56wdy8MGMBHjz8e9Is3qs+LiIhITZeXZ0zezcjgdyAeWPvnS4FavFF9XkRERKT8oqIgLY28Nm1oBHwFdPjzpWBcvFHhRURERMBuZ99//8tKIBojwJxPycUbg2E+jMKLiIiIAHD2+efz80svsQo4E1gIzH7ssaJLRi6XKyjmw2jOi4iIiBSzdd06ooYNo0FWFjRuDIsW4W7QAKfTidfrLdqvOufDaM6LiIiIVNrZHTrQICMDOnWCHTugXz+2fvVVseAC5s2HUXgRERGRkho2hC+/hI4dYft2LnrgAVpbLMV2OXE+TKAovIiIiEjpGjWCr76CDh2w7djBarud1lYjOthsNlJSUvx+C3VpIgL+HUVERCR0REcbAaZfP+r/+CM/NGnC6uef5+w+fUwJLqCRFxERETmVM8+EhQuhXTsicnPpPmECjoIC08pReBEREZFTa9zYCDBt2kCrVkagMYkuG4mIiEj5xMTAokXQoAHUq2daGX4deXnqqafo2bMn9erV4/TTTy/XMT6fj0mTJtG0aVPq1q1L3759+fHHH/1ZpoiIiJRXTIypwQX8HF4OHz7Mddddx2233VbuY5599ln+9a9/MXXqVFauXEmTJk249NJL2bdvnx8rFRERkVDh1/DyxBNPMH78eDp06HDqnTFGXV588UUefvhhrrnmGtq3b89bb73FgQMHeOedd/xZqoiIiISIoJqwu2nTJnJzc0lISCjaFhkZSZ8+fVi6dGmpx+Tn55OXl1fsISIiIuErqMJLbm4uADExMcW2x8TEFL12osmTJ2O324sezZo183udIiIiYp4Kh5dJkyZhsVjKfKxatapKRVlOaD/s8/lKbCs0ceJEPB5P0WPLli1V+t4iIiIS3Cp8q/S4ceMYPnx4mfu0aNGiUsU0adIEMEZgzjrrrKLtO3bsKDEaUygyMpLIyMhKfT8REREJPRUOL9HR0URHR/ujFs455xyaNGnCggULuPDCCwHjjqXFixfzzDPP+OV7ioiISGjx65yXzZs3k5mZyebNmykoKCAzM5PMzEz2799ftE+bNm2YM2cOYFwuuvvuu3n66aeZM2cOP/zwA7fccgv16tVjxIgR/ixVREREQoRfO+w+9thjvPXWW0XPC0dT0tPT6du3LwAbNmzA4/EU7fPAAw9w8OBBxo4dy549e+jevTvz58+nQYMG/ixVREREQoTF5/P5zC6iOuXl5WG32/F4PERFRZldjoiIiJRDRT6/g+pWaREREZFTUXgRERGRkKLwIiIiIiHFrxN2zVA4hUfLBIiIiISOws/t8kzFDbvwUrj6tJYJEBERCT379u3DbreXuU/Y3W3k9XrZtm0bDRo0KHVJga5du7Jy5coy36OsfSrz2onb8/LyaNasGVu2bAmKO6LK8zsJ1PtV5Fidy5Jq6rks63Wdy6q/n85l1ehclu9c+nw+unTpwsaNG7Fay57VEnYjL1arFYfDcdLXbTbbKf9lLmufyrx2su1RUVFB8R9WeX4ngXq/ihyrc1lSTT2XZb2uc1n199O5rBqdy/Kfs9q1a58yuEANnLB7++23V2mfyrxWnu9ppuquryrvV5FjdS5LqqnnsqzXdS6r/n46l1Wjc1n17ScKu8tGoUCN9MKHzmX40LkMHzqX4a/GjbwEg8jISB5//HGthh0GdC7Dh85l+NC5DH8aeREREZGQopEXERERCSkKLyIiIhJSFF5EREQkpCi8iIiISEhReBEREZGQovASAg4cOIDT6eS+++4zuxSpgn379tG1a1c6depEhw4dmDZtmtklSSVt2bKFvn370q5dOzp27MgHH3xgdklSBUOHDuWMM87g2muvNbsUKSfdKh0CHn74YbKysmjevDnPP/+82eVIJRUUFJCfn0+9evU4cOAA7du3Z+XKlTRq1Mjs0qSCfvvtN7Zv306nTp3YsWMHnTt3ZsOGDdSvX9/s0qQS0tPT2b9/P2+99Rb//e9/zS5HykEjL0EuKyuLn3/+mcGDB5tdilSRzWajXr16ABw6dIiCgoJyLf0uweess86iU6dOADRu3JiGDRuye/duc4uSSuvXrx8NGjQwuwypAIWXKliyZAlXXHEFTZs2xWKx8NFHH5XY55VXXuGcc86hTp06dOnSha+//rpC3+O+++5j8uTJ1VSxlCUQ53Pv3r1ccMEFOBwOHnjgAaKjo6upejleIM5loVWrVuH1emnWrFkVq5bSBPJcSuhQeKmCP/74gwsuuICpU6eW+vrs2bO5++67efjhh1mzZg29e/fmsssuY/PmzUX7dOnShfbt25d4bNu2jY8//phWrVrRqlWrQP1INZq/zyfA6aefzvfff8+mTZt455132L59e0B+tpomEOcS4Pfff+fmm28mNTXV7z9TTRWocykhxifVAvDNmTOn2LZu3br5xowZU2xbmzZtfBMmTCjXe06YMMHncDh8TqfT16hRI19UVJTviSeeqK6SpQz+OJ8nGjNmjO/999+vbIlSTv46l4cOHfL17t3bN2PGjOooU8rBn/9dpqen+4YNG1bVEiVANPLiJ4cPH2b16tUkJCQU256QkMDSpUvL9R6TJ09my5Yt5OTk8PzzzzNq1Cgee+wxf5Qrp1Ad53P79u3k5eUBxqq3S5YsoXXr1tVeq5StOs6lz+fjlltuIT4+nsTERH+UKeVQHedSQlOE2QWEq127dlFQUEBMTEyx7TExMeTm5ppUlVRWdZxPt9tNUlISPp8Pn8/HuHHj6Nixoz/KlTJUx7n89ttvmT17Nh07diyagzFz5kw6dOhQ3eVKGarr7+zAgQP57rvv+OOPP3A4HMyZM4euXbtWd7lSjRRe/MxisRR77vP5Smwrj1tuuaWaKpKqqMr57NKlC5mZmX6oSiqjKueyV69eeL1ef5QllVDVv7NffPFFdZckfqbLRn4SHR2NzWYrkf537NhR4v8SJPjpfIYPncvwoXNZcym8+Ent2rXp0qULCxYsKLZ9wYIF9OzZ06SqpLJ0PsOHzmX40LmsuXTZqAr2799PdnZ20fNNmzaRmZlJw4YNad68Offccw+JiYlcdNFFxMXFkZqayubNmxkzZoyJVcvJ6HyGD53L8KFzKaUy8U6nkJeenu4DSjz++te/Fu3z8ssv+5xOp6927dq+zp07+xYvXmxewVImnc/woXMZPnQupTRa20hERERCiua8iIiISEhReBEREZGQovAiIiIiIUXhRUREREKKwouIiIiEFIUXERERCSkKLyIiIhJSFF5EREQkpCi8iIiISEhReBEREZGQovAiIiIiIUXhRURERELK/wPPChFFZHfP3AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hm = ml.head(robs, 0, tobs, 0)\n",
    "plt.semilogx(tobs, hobs, \".k\")\n",
    "plt.semilogx(tobs, hm[0], \"r\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "correlation matrix\n",
      "[[ 4.49139190e-01 -2.49689096e-06]\n",
      " [-2.49689096e-06  1.93916909e-11]]\n"
     ]
    }
   ],
   "source": [
    "print(\"correlation matrix\")\n",
    "print(cal.fitresult.covar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit with `scipy.least_squares` (not recommended)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..........................................\n",
      "      layers    optimal       std  perc_std  pmin  pmax  initial inhoms  \\\n",
      "kaq0       0  59.870969  0.660661  1.103474  -inf   inf  100.000   None   \n",
      "Saq0       0   0.000121  0.000004  3.587586  -inf   inf    0.001   None   \n",
      "\n",
      "                        parray  \n",
      "kaq0       [59.87096932924763]  \n",
      "Saq0  [0.00012119340200370668]  \n",
      "[6.60660677e-01 4.34791694e-06]\n",
      "[[ 4.36472531e-01 -2.42866203e-06]\n",
      " [-2.42866203e-06  1.89043817e-11]]\n",
      "[[ 1.         -0.84548788]\n",
      " [-0.84548788  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "cal = ttim.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq0\", layers=0, initial=100)\n",
    "cal.set_parameter(name=\"Saq0\", layers=0, initial=1e-3)\n",
    "cal.series(name=\"obs1\", x=robs, y=0, layer=0, t=tobs, h=hobs)\n",
    "cal.fit_least_squares(report=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibrate parameters in multiple layers\n",
    "Example showing how parameters can be optimized when multiple layers share the same parameter value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self.neq  1\n",
      "solution complete\n"
     ]
    }
   ],
   "source": [
    "ml = ttim.ModelMaq(\n",
    "    kaq=[10.0, 10.0],\n",
    "    z=(-10, -16, -18, -25),\n",
    "    c=[10.0],\n",
    "    Saq=[0.1, 1e-4],\n",
    "    tmin=1e-5,\n",
    "    tmax=1,\n",
    ")\n",
    "w = ttim.Well(ml, xw=0, yw=0, rw=0.1, tsandQ=[(0, 788)], layers=1)\n",
    "ml.solve()\n",
    "hobs0 = ml.head(robs, 0, tobs, layers=[0])[0]\n",
    "hobs1 = ml.head(robs, 0, tobs, layers=[1])[0]"
   ]
  },
  {
   "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>layers</th>\n",
       "      <th>optimal</th>\n",
       "      <th>std</th>\n",
       "      <th>perc_std</th>\n",
       "      <th>pmin</th>\n",
       "      <th>pmax</th>\n",
       "      <th>initial</th>\n",
       "      <th>inhoms</th>\n",
       "      <th>parray</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>kaq0</th>\n",
       "      <td>0</td>\n",
       "      <td>59.870969</td>\n",
       "      <td>0.660661</td>\n",
       "      <td>1.103474</td>\n",
       "      <td>-inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>100.000</td>\n",
       "      <td>None</td>\n",
       "      <td>[59.87096932924763]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.000121</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>3.587586</td>\n",
       "      <td>-inf</td>\n",
       "      <td>inf</td>\n",
       "      <td>0.001</td>\n",
       "      <td>None</td>\n",
       "      <td>[0.00012119340200370668]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      layers    optimal       std  perc_std  pmin  pmax  initial inhoms  \\\n",
       "kaq0       0  59.870969  0.660661  1.103474  -inf   inf  100.000   None   \n",
       "Saq0       0   0.000121  0.000004  3.587586  -inf   inf    0.001   None   \n",
       "\n",
       "                        parray  \n",
       "kaq0       [59.87096932924763]  \n",
       "Saq0  [0.00012119340200370668]  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cal.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "....................................................................................................................\n",
      "Fit succeeded.\n"
     ]
    },
    {
     "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>layers</th>\n",
       "      <th>optimal</th>\n",
       "      <th>std</th>\n",
       "      <th>perc_std</th>\n",
       "      <th>pmin</th>\n",
       "      <th>pmax</th>\n",
       "      <th>initial</th>\n",
       "      <th>inhoms</th>\n",
       "      <th>parray</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>kaq0_1</th>\n",
       "      <td>[0, 1]</td>\n",
       "      <td>9.998813</td>\n",
       "      <td>3.275126e-04</td>\n",
       "      <td>0.003276</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>30.0</td>\n",
       "      <td>20.000</td>\n",
       "      <td>None</td>\n",
       "      <td>[9.998812542977381, 9.998812542977381]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.100010</td>\n",
       "      <td>1.041096e-07</td>\n",
       "      <td>0.000104</td>\n",
       "      <td>0.00001</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.001</td>\n",
       "      <td>None</td>\n",
       "      <td>[0.10000993019005114, 9.999506259756452e-05]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Saq1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.000100</td>\n",
       "      <td>1.175972e-09</td>\n",
       "      <td>0.001176</td>\n",
       "      <td>0.00001</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0.001</td>\n",
       "      <td>None</td>\n",
       "      <td>[9.999506259756452e-05]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>c1</th>\n",
       "      <td>1</td>\n",
       "      <td>9.999649</td>\n",
       "      <td>9.921790e-05</td>\n",
       "      <td>0.000992</td>\n",
       "      <td>0.10000</td>\n",
       "      <td>200.0</td>\n",
       "      <td>1.000</td>\n",
       "      <td>None</td>\n",
       "      <td>[9.999649350868319]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        layers   optimal           std  perc_std     pmin   pmax  initial  \\\n",
       "kaq0_1  [0, 1]  9.998813  3.275126e-04  0.003276  0.00000   30.0   20.000   \n",
       "Saq0         0  0.100010  1.041096e-07  0.000104  0.00001    0.2    0.001   \n",
       "Saq1         1  0.000100  1.175972e-09  0.001176  0.00001    0.2    0.001   \n",
       "c1           1  9.999649  9.921790e-05  0.000992  0.10000  200.0    1.000   \n",
       "\n",
       "       inhoms                                        parray  \n",
       "kaq0_1   None        [9.998812542977381, 9.998812542977381]  \n",
       "Saq0     None  [0.10000993019005114, 9.999506259756452e-05]  \n",
       "Saq1     None                       [9.999506259756452e-05]  \n",
       "c1       None                           [9.999649350868319]  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cal = ttim.Calibrate(ml)\n",
    "cal.set_parameter(\n",
    "    name=\"kaq0_1\", layers=[0, 1], initial=20.0, pmin=0.0, pmax=30.0\n",
    ")  # layers 0 and 1 have the same k-value\n",
    "cal.set_parameter(name=\"Saq0\", layers=0, initial=1e-3, pmin=1e-5, pmax=0.2)\n",
    "cal.set_parameter(name=\"Saq1\", layers=1, initial=1e-3, pmin=1e-5, pmax=0.2)\n",
    "cal.set_parameter(name=\"c1\", layers=1, initial=1.0, pmin=0.1, pmax=200.0)\n",
    "cal.series(name=\"obs0\", x=robs, y=0, layer=0, t=tobs, h=hobs0)\n",
    "cal.series(name=\"obs1\", x=robs, y=0, layer=1, t=tobs, h=hobs1)\n",
    "cal.fit(report=False)\n",
    "display(cal.parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGhCAYAAACphlRxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABno0lEQVR4nO3deVyU5f7/8dcMqwKOIIIrromouKeiadmCSy5pnfJnBzONk5VHbbHSTuU5nW+2Z9Zps7QyO21mJz2mdTKXUtxJzX0FF0RTBwUEgfv3xy2jI4ugDMPA+/l43I+5576v+57PMMq8ue/rvi+LYRgGIiIiIh7C6u4CREREREpD4UVEREQ8isKLiIiIeBSFFxEREfEoCi8iIiLiURReRERExKMovIiIiIhH8XZ3AWUtLy+Pw4cPExQUhMVicXc5IiIiUgKGYXD69Gnq1auH1Vr8sZVKF14OHz5Mw4YN3V2GiIiIXIHk5GQaNGhQbJtKF16CgoIA883XqFHDzdWIiIhISaSlpdGwYUPH93hxKl14yT9VVKNGDYUXERERD1OSLh/qsCsiIiIeReFFREREPIrCi4iIiHgUhRcRERHxKAovIiIi4lEUXkRERMSjKLyIiIiIR1F4EREREY+i8CIiIiIeReFFREREPEq5hJe3336bJk2a4O/vT6dOnVixYkWx7ZctW0anTp3w9/enadOmvPvuu+VRpoiIiHgAl4eXL774ggkTJvDUU0+xceNGevbsSb9+/UhKSiq0/b59++jfvz89e/Zk48aNTJ48mXHjxjF37lxXl3pZRw/uYcuv8zl6cI+7SymVilJ3WdRxxJ7Jyj3HOWLPvKpaymI/V7OPkm5blu0u16Y83k9p2+e3+y35ZJHtr7Tu8ny/IlK2LIZhGK58ga5du9KxY0feeecdx7KoqChuu+02pk6dWqD9E088wXfffce2bdscy8aMGcNvv/3GqlWrLvt6aWlp2Gw27HZ7mQ7MuPqr16m96R28LAa5hoU9Tf4fLWIGldn+XWXnqu9otu/fF9U9nGu6Dz6/1gIXD4DlmLdgXDRvTvnrLZe0BcNpP5fMn2+/dflXRO6eidUCOYaV31s8RHTsiFK9l++3HGH2j2upZznOYSOUuFuupV+buqXaR/5+Xv9xF8b56h6+5ZpS7+dq9lHSbQtr1z+6YLuFm53bPXLLNfSLrue8r82Hee2iNo/GtqBfdN38T5OFm4/w6g87Hesfi21B/7b1uHh4NKd/KufXWCywYNNhXlq0w7HtE30jGdCunmNwNctF21qw8N1vh5i6cLuj/VO3tmRw+wZYLPltzb1/m3iQf8zfxsW/oCzA3we34o5ODbFaLHy9Ppmnv/3dsa/nh7Zh2LURlx3Y7Yu1SUz6ZjN5BlgtMHVoNHddG1HsNmWxbXGO2DPZdzydJqEB1LVVu+p2JWlbmn2JuFppvr9dGl6ys7OpXr06X331FUOGDHEsHz9+PImJiSxbtqzANr169aJDhw688cYbjmXz5s3jzjvvJCMjAx8fH6f2WVlZZGVlOZ7nD6ldluHl6ME91Hj/WqpZc8tkf2LKMPzIxJdM/Dhr+JKBn2M+k/PrDF/OOubNx1NGEDuMBuwyGpCFr7vfhlRA+UHIarFgtVjM5xbzOUBGdsH/y6GBvvh6WbFYLHhZLVgtYLWa23ud30eeYbDz6JkC217bOJjqvt54W81tvb0seFmtF547PVrPrzefWy0Wth1J48etRx0h7LYO9ejWtBY+Xla8vaz4WC34eFlZuec4s37d72g37qbmDGxXHz9vK77eVny9zEcfLyvfbEhm8rwtRYas0oQwhRwpD6UJL96uLOT48ePk5uYSHh7utDw8PJyUlJRCt0lJSSm0fU5ODsePH6duXee/PKdOncrf//73si38EscObCXYYnDSCHRabmAp0dDd7mIYBhanv1uN/L+Xz/9yN84/u9DGgnHR8ovXGY6/wi9tw0X7ubiNleJzcXVLFtXJAk7DFfwYcw0L+6nHbmsjdtKYXRbzMZUQ58MEFzmXm8fpszkFlgf5e+PjVbKzqFezj6K2DfTzxsfLclE7gzNZV9fO+3y7nCLaBPh64WW1kJNnFPplXs3HC2/rxZ+wKf/vHQPIzTPIyskrsK2Pl/mlbFy0sYFBnmFu42qGYb5snmHAZf4d5jt+JvuKX2/t/pNXvO2lDGDexsPM23j4su3e+Gk3b/y0+7L7zDPgibmbeXvpHgJ8vbFaYcuhNKf1T87dzMakU4QG+lHN1wt/Hy+q+Xjx28GTfLn2oCMwjb/5Gm7v2IBAP28C/Lzx9S7+33xJgo/CkZSWS8NLvku/4A3DKPZLv7D2hS0HmDRpEo888ojjef6Rl7JUu1ErvDAItlz4iyvHsPJH/DrCGzQr09cqS0cP7iF0Rie8LBd+ebujbrOOjlz0nUuOYeHksO+oXas2nMuAc5mQnXFh3vGYCefSyUzZQbXdC53262UxaMYhmhmH6MPKC99R1YIhvA2EtzYf67SB2i3BpxpH7Jn0eGEJYcYfNLGmsC+vDscsofzwcK8S/9LM38fF38FeFkuJ9lHUtj8+4rxtWbYrqs3/Hr2eurZqRa5f8tj1V/x+lj/eu8jTFIW1/+VJs71hGBgGHLZn0vPFnwuNHVbgx0d7YRhwy2vLndpYgQXjr6N2oD+GYYal/NCUd/5FU9IyufPdhALbzbz3WmoF+JFrGOQZBnl5F8KWYRjkGgbHz2TxyBe/FTiV9eygVgT5+ZCbZ5CTZ5Cbl3f+0bjwmFv48oMnM/jfttQC77NjRE0C/LzJyTU4l5vHyYxs9hxLL9AuwNeLXMMgOyePy+XCA39kFLnOAD5fm1zs9gYw7X+7mPa/XY5lvl5Wqvt5EeDrfT7QeBHgZ86nns5i/YGTjp/T0I71uTkqnBrVfKjh74Otmg8/bU9xnB4sy9NwUrm5NLyEhobi5eVV4ChLampqgaMr+erUqVNoe29vb2rVqlWgvZ+fH35+fmVXdCHCGzRjTdspdNz0d7wteeQYVja0fZYuFTi4QMWp26zj7wXriOpV4n1Usx/CeH0RFi78lZ+HBeuQd+H0ETj6O6RsgeM7IfMk7F9hTvksVqh1DXXDW7Mswk69lJ/wsphHb9a3nUJdW/8S11LXVo2pQ6OZ/M0Wcg0DL4uF54e2KVH4Kem2Zdnucm3K4/2UtL3l/OmZBsHVeeH2C+3y5bdvVjsIwKlN/rpWdW3F1twwpHqh290QGXbZ9wuQnZNXYNur+bI9Ys9kyfaCge5fd3csUVDND6EAObl5ZOfmkXQig37TVhQIaG/e3YFAPx+O2DOZNHdzgRD2524ReFmtZGbnknkul0OnMh3h42I+XhbO5ZpbZ+fmkZ2Rx6mMc8W+TwOYu+EQczccKrJN/hGiRVtSqB9cjZAAP0IDfakV4EdIgK85H+hHzWo+WK2F/wGsozhVQ7l02O3UqRNvv/22Y1mrVq0YPHhwkR1258+fz9atWx3LHnjgARITE93aYRfMIwjHD2wntFHLCn3E5VIVpe6rrmPDJxjzJ2AxcjEsXlgGToOOl3T6zcmCY9svhJmj56eMP4rZsRUe3gK2+qUq54g9k/3HM2gcWr3UvyRLum1Ztrtcm/J4P6Vtn9+uuq+VjOy8Qttfad3l+X4v54u1SSUKRCVtV5K2JdlXcUfKQgP9yMjK5Ux2DulZOZzJyjGfZ5nPNx86xUcrDxSoKzI8CAODtMwcTmZkF3ra8XKsFggJuBBqagX6Ehrox6FTGfxva6rjFNezg1oxsnuTy+5PgadiqDAddsG8VDouLo53332XmJgY3n//fWbMmMHvv/9Oo0aNmDRpEocOHeKTTz4BzEul27Rpw/333098fDyrVq1izJgx/Pvf/+b222+/7Ou5MrxIBWA/BCf2QkjTkocNw4AzR80ws30+rP+oYJt2/w/6vQT++jcj7lHWgbYkbUuyr9IEpkv3Xdwpwvw23acuKXAE6OFbWpCTZ/DHmSz+OJPNifRsjqeb8/bM4o/wXCrAz4sGNatTr6Y/9WpWo17NatQ//1ivpj/Ldx7jb98W3bFZyk+FCi9g3qTupZde4siRI7Rp04bXX3+dXr3MUwYjR45k//79LF261NF+2bJlPPzww/z+++/Uq1ePJ554gjFjxpTotRRepFj2QzCtDRiF/LVXLQR6PgLX3gc++utLJN+VHmkqSfApbTg6l5vHyfRsjp/J5o/0LDPYnMnmt4On+C6x+E7OJWEBXr2zHdc2DqFezWp4FXF6SspehQsv5UnhRS5rwycwfwIYuYAVOo+EfSvgj/OdEIPqwQ1PQPs/g1e59GkXqbTK4pRmSV/n0iM9VmDOX7pyLsfg8KlMDp/K5NCpsxw+lckReyYHT2aSU0wvZz9vK01CA2gWFkiz2oE0qx1As9qBNK0dQHVfb8fr6pRT2VB4UXiRy7n09FNuDvz2b1j6AqQdNNuENIMbn4JWQ8CqYcBEKrrSHsU5dCqD614oeFVb09oBHDyRSXZu0f1x6tn8CfDzZleqeRWqBfjnkNbc3bXx1b+RKkrhReFFrtS5s7BuJqx45UIn3zpt4aZnoPnNRd4/RkQqhtIexSkq8ORfxr7n2Bn2Hktnz7Ez7Ek1H/9IL/qeQC3rBNGxUTDR9W1E17fRIjwIX2+rjtCUgMKLwotcrazTsOptWPkmZJ82l0V0h5ufhYhu7q1NRMpUaQPPyfRs/vPbIaZ8t/WybX29rITV8OPgSXMcLAswdWgbhnVpdLVlVzoKLwovUlbS/4BfXoM1MyD3/DAU1/SBm542O/ie2GOeXirlZdYi4tmK6mPzj9vacPBkJlsO2dl08BRphdxRG+CWVuHcEFmbrk1CaFY7EIvFUuWPzii8KLxIWbMfgmUvwsZPz3f0vYjFCgPfKHjPGRGp1C7Xx8YwDL5NPMTDX/xW7H5qBfhSt6Y/vx9Kc9yj5oXbq94l2wovCi/iKsd3w49Pww7noQqweMGEzToCI1LFlOReOoUdoRl1XRN+P5zGhqSThd6ozwL86+4OxLaqg3cJx13zdAovCi/iSvuWw8cDCy6/ZwE06Vn+9YhIhVbcEZqsnFw+W53E3+cX3n8muLoPN0WFE9sqnJ7X1Kaar1d5ll6uKsyo0iKVUkgz81TRpTe68/Z3Tz0iUqHddW0EvVrULvQIjZ+3F33b1OG5BVudjs5YMEeqP5lxjq/XH+Tr9Qfx97HS65raxLauw00twwgO8HW0r2r9ZXTkReRKON3o7ryQpnDvIggqfNBREZGiFHZ05vaODVi7/yQ/bE3hh9+PcuhUpqO91QLXNg4htnUdsnNyeXnxDo8f4kCnjRRepDzk3+jOpzp8NRLsSRDWCkb+F6qHuLs6EfEwxfWfMQyDrUfS+OH3o/yw9SjbjqQVuZ9Lx5DyFAovCi9S3k7shZn94EwK1OsAI77TII8i4jLJJzL4cetRvl6fzNYjpwusf3ZgK0Z2b4zFg26sqfCi8CLukLodZvWDzBPQqAfc/TX4Vnd3VSJSiRU2Mne+FuGBDO8SwZCODbBV8yn32kqrNN/fVeP6K5HyENYS4uaBXw048Ct8GQc5We6uSkQqsbq2arxwezRe54+wWIHOjYPx97Gy8+gZpszfStfn/8fEr34jMfkUhmFwxJ7Jyj3HOWLPLH7nFZiOvIiUtaQEmD0EzmVA1EC44yONTi0iLnVpfxl75jm+3XiIOasPsPPoGUe7ejZ/jtjPYlDxOvfqtJHCi7jbnp/hszshNxvaDoPb3tHI1CJS7gzDYP2Bk8xZncR/Nx0pMFK2Ffh10o0VonOvThuJuFuz3vCnj8w77276HBY+BpXr7wQR8QAWi4XOjUN4/a72vDm8fYH1ecCbP+0mI7vwMZgqKoUXEVdpeSsMeQ+wwLoP4X/PKsCIiNu0bVATayEXH322JomeL/7Me8v2eEyIUXgRcaW2f4IBr5vzv74BK15xbz0iUmXVtVVj6lDnzr1/6tSAiJDq/JGezdTvt9PzxZ95f/ke9h47U6E79arPi0h5WPkW/PCUOd/3Bej2gHvrEZEq69LOvedy85i38RBvLdlN0okMp7blOcK1OuwqvEhFtPQFWDrVnL/lOajX3hwnSSNRi0gFcC43j49W7uP//rvdabkV+OXJ3tSr6dr7VqnDrkhFdP0TEDPWnP/xaXNk6mltzHGSRETczMfLSut6tgLL84AH52xk3/H08i+qCAovIuXFYoGul5wuMvLMAR7th9xSkojIxZqEBhTaqTcx+RR9Xl/Oi4u2s+fYabf3h1F4ESlPJ/cWXGbkmmMjiYi42aWder0sFh6LbcH1LWqTnZvHO0v3cNOryxk+YzXdpy7hi7VJbqlTt/0UKU8hzcBiNY+4OFggpKnbShIRudhd10bQq0Vtp069hmHw1fpkHv96s6OdAUz+Zgu9WtQu95vc6ciLSHmy1YeBb5g3r8vn5Q25GgNJRCqOurZqxDSr5QglFouFBsEFO+zmGgb7j2cUWO5qCi8i5a3jCJiwGUbMh4ZdIfec2e+lcl34JyKVTGH9YbwsFhqHuvYqpMIovIi4g60+NO1ljnnk7Q/7lsFv/3Z3VSIiRSqsP8zzQ9u4ZVwk3edFxN1+eR3+NwWqBcNDayGwtrsrEhEp0qU3uSsrus+LiCeJGQt1oiHzJCx60t3ViIgU69L+MO6g8CLibl4+MOhN8yqkLV/Dzh/cXZGISIWm8CJSEdTrAN0eNOcXPAxZp91bj4hIBabwIlJR9J4MNRtB2kFY8k93VyMiUmEpvIhUFL4BMHCaOb/6PTi4zq3liIhUVAovIhVJsxuh3f8DDPjur5CT7e6KREQqHIUXkYom9v+gei1I3Qor33B3NSIiFY7Ci0hFE1AL+r5ozi97ERI/06jTIiIXUXgRqYii74CwVubQAd8+AK+3hg2fuLsqEZEKQeFFpCJKOwyp2y5aYJjjH+kIjIiIwotIhXRiD+aA8xcxcuHEXreUIyJSkbg0vJw8eZK4uDhsNhs2m424uDhOnTpV7DbffPMNffr0ITQ0FIvFQmJioitLFKmYQpqZd9x1YoGQpm4pR0SkInFpeBk+fDiJiYksWrSIRYsWkZiYSFxcXLHbpKen06NHD1544QVXliZSsdnqw8A3wOJ1YVm1YAjQoI0iIt6u2vG2bdtYtGgRCQkJdO3aFYAZM2YQExPDjh07iIyMLHS7/HCzf/9+V5Um4hk6joBmN5l9X74dA+nHYONsuHa0uysTEXErlx15WbVqFTabzRFcALp164bNZmPlypVl9jpZWVmkpaU5TSKVhq0+XHMz9HrcfL78FTh31r01iYi4mcvCS0pKCmFhYQWWh4WFkZKSUmavM3XqVEefGpvNRsOGDcts3yIVRscRUKM+nD4M6z9ydzUiIm5V6vAyZcoULBZLsdO6deaYLBaLpcD2hmEUuvxKTZo0Cbvd7piSk5PLbN8iFYaPP/R6zJz/5TXIznBvPSIiblTqPi9jx45l2LBhxbZp3LgxmzZt4ujRowXWHTt2jPDw8NK+bJH8/Pzw8/Mrs/2JVFjt/wy/vA6nkmDdh9D9r+6uSETELUodXkJDQwkNDb1su5iYGOx2O2vWrKFLly4ArF69GrvdTvfu3UtfqUhV5+1r9n35biz8Mg063Qt+ge6uSkSk3Lmsz0tUVBR9+/YlPj6ehIQEEhISiI+PZ8CAAU5XGrVs2ZJ58+Y5np84cYLExES2bt0KwI4dO0hMTCzTfjIiHqvdMAhuAhnHYc377q5GRMQtXHqflzlz5hAdHU1sbCyxsbG0bduW2bNnO7XZsWMHdrvd8fy7776jQ4cO3HrrrQAMGzaMDh068O6777qyVBHP4OUDNzxpzq+cDmd1dZ2IVD0WwzCMyzfzHGlpadhsNux2OzVq1HB3OSJlLy8X/tUV/tgFvZ+C6x93d0UiIletNN/fGttIxNNYvS46+vIWZJ5yazkiIuVN4UXEE7UeCrWjIMsOq/7l7mpERMqVwouIJ7Jaofckcz7hHcg44d56RETKkcKLiKdqORDCoyH7NPz8f7BvOdgPubsqERGXU3gR8VRWK/SebM6v/QA+HgjT2sCGT9xbl4iIiym8iHiyOm2dnxt5MH+CjsCISKWm8CLiyU7uLbjMyIUThSwXEakkFF5EPFlIM+CSgU4tXhDS1C3liIiUB4UXEU9mqw+DpuMUYAa8bi4XEamkFF5EPF3HEfBgAnhXM5/bGri3HhERF1N4EakMwlqaIQZgtcYBE5HKTeFFpLLoej9ggV0/wPHd7q5GRMRlFF5EKotazeCaWHN+zXvurUVExIUUXkQqk25jzMfEz+Cs3b21iIi4iMKLSGXStDfUbgnZZ2Djp+6uRkTEJRReRCoTiwW6nj/6svo9yMt1bz0iIi6g8CJS2bS9C6oFw6kDsHORu6sRESlzCi8ilY1vdeh4jzmf8I57axERcQGFF5HKqEu8OUzA/hWQssXd1YiIlCmFF5HKyNYAogaa87ppnYhUMgovIpVVtwfMx81fQfof7q1FRKQMKbyIVFYNu0Ld9pBzFtbPcnc1IiJlRuFFpLKyWC4cfVn7AeSec289IiJlROFFpDJrPQQCwuD0Edj6H3dXIyJSJhReRCozbz+49j5z/tc3YN9ysB9yb00iIldJ4UWksut8L1i9IWUTfDwQprWBDZ+4uyoRkSum8CJS2eWeg7ycC8+NPJg/QUdgRMRjKbyIVHYn9hRcZuTCib3lX4uISBlQeBGp7EKaARbnZRYvCGnqlnJERK6WwotIZWerD4Om4xRgBk4zl4uIeCCFF5GqoOMIGPOr2XEXoE5b99YjInIVFF5Eqoo6rSFqkDm/cbZ7axERuQoKLyJVScc483HTV3Au0721iIhcIYUXkaqkyQ1gi4AsO2yb7+5qRESuiMKLSFVitUKHP5vzulGdiHgohReRqqb9cMAC+1foXi8i4pEUXkSqmpoNodmN5vzGT91bi4jIFVB4EamK8jvuJn4GuTnFtxURqWAUXkSqosj+UL0WnD4Ce35ydzUiIqWi8CJSFXn7Qdth5rw67oqIh3FpeDl58iRxcXHYbDZsNhtxcXGcOnWqyPbnzp3jiSeeIDo6moCAAOrVq8eIESM4fPiwK8sUqZryTx3tXARnUt1bi4hIKbg0vAwfPpzExEQWLVrEokWLSExMJC4ursj2GRkZbNiwgaeffpoNGzbwzTffsHPnTgYNGuTKMkWqprAoqN8Z8nLgt3+7uxoRkRKzGIZhuGLH27Zto1WrViQkJNC1a1cAEhISiImJYfv27URGRpZoP2vXrqVLly4cOHCAiIiIAuuzsrLIyspyPE9LS6Nhw4bY7XZq1KhRNm9GpLJa/zHMHwe1roGxa8Fiufw2IiIukJaWhs1mK9H3t8uOvKxatQqbzeYILgDdunXDZrOxcuXKEu/HbrdjsVioWbNmoeunTp3qOC1ls9lo2LDh1ZYuUnW0GQo+AfDHLkhKcHc1IiIl4rLwkpKSQlhYWIHlYWFhpKSklGgfZ8+e5cknn2T48OFFprBJkyZht9sdU3Jy8lXVLVKl+AVB6yHmvAZrFBEPUerwMmXKFCwWS7HTunXrALAUcgjaMIxCl1/q3LlzDBs2jLy8PN5+++0i2/n5+VGjRg2nSURKIb/j7u/z4Gyae2sRESkB79JuMHbsWIYNG1Zsm8aNG7Np0yaOHj1aYN2xY8cIDw8vdvtz585x5513sm/fPpYsWaJAIuJKDbtCaAs4vhOWvQjdHgRbfXdXJSJSpFKHl9DQUEJDQy/bLiYmBrvdzpo1a+jSpQsAq1evxm6307179yK3yw8uu3bt4ueff6ZWrVqlLVFESsNigbBWZnhZ9RYkvA0D34COI9xdmYhIoVzW5yUqKoq+ffsSHx9PQkICCQkJxMfHM2DAAKcrjVq2bMm8efMAyMnJ4Y477mDdunXMmTOH3NxcUlJSSElJITs721WlilRt9kOw9T8Xnht5MH+CuVxEpAJy6X1e5syZQ3R0NLGxscTGxtK2bVtmz3buFLhjxw7sdjsABw8e5LvvvuPgwYO0b9+eunXrOqbSXKEkIqVwYg9wyR0TjFyNOC0iFVapTxuVRkhICJ9+WvyotRffZqZx48a46LYzIlKUkGZgsZpHXBysENLUbSWJiBRHYxuJVHW2+mYfl4t/HXQfq067IlJhKbyIiNk59+EtcE2s+fzsKbeWIyJSHIUXETHZ6kP3v5rzv/8Hzp11bz0iIkVQeBGRCxpdBzXqQ5Yddv3g7mpERAql8CIiF1it0OZ2c37zl+6tRUSkCAovIuKs7Z3m487FkHnKraWIiBRG4UVEnIW3gdpRkJsN275zdzUiIgUovIiIM4sF2v7JnN+kU0ciUvEovIhIQdHnw8v+XyDtsHtrERG5hMKLiBRUMwIiYgADNn/t7mpERJwovIhI4fKPvuiqIxGpYBReRKRwrYeA1RtSNkPqdndXIyLioPAiIoWrHgLNbzHndfRFRCoQhRcRKVr+PV82fQV5ecW3FREpJwovIlK0yH7gGwT2JEhe7e5qREQAhRcRKY5PNYgaaM7r1JGIVBAKLyJSvPwb1v0+D3Ky3VuLiAgKLyJyOU2uh8BwyDwJe35ydzUiIgovInIZVq8LI01ruAARqQAUXkTk8vJvWLfje8g67d5aRKTKU3gRkcur1wFqNYecTNi2wN3ViEgVp/AiIpdnsUD0+Xu+6KojEXEzhRcRKZnoO8zHvUvh9FG3liIiVZvCi4iUTK1m0OBaMPJgy1x3VyMiVZjCi4iUnE4diUgFoPAiIiXXeghYvODwRji+293ViEgVpfAiIiUXWBua3WjO6+iLiLiJwouIlI5jpOkvwTDcW4uIVEkKLyJSOpH9wac6nNwHh9a7uxoRqYIUXkSkdPwCoeWt5ryGCxARN1B4EZHSy7/q6PdvIDfHvbWISJWj8CIipdesN1SvBenHzJvWiYiUI4UXESk9Lx9oPdSc11VHIlLOFF5E5Mq0vct83LYAstPdW4uIVCkKLyJyZRp0huAmcC4dti90dzUiUoUovIjIlbFYIPpP5rxOHYlIOVJ4EZErl3/Dut0/wbb5YD/k3npEpEpQeBGRKxd6DdgagpELX/wZprWBDZ+4uyoRqeQUXkTkytkPgf3ghedGHsyfoCMwIuJSLg0vJ0+eJC4uDpvNhs1mIy4ujlOnThW7zZQpU2jZsiUBAQEEBwdz8803s3r1aleWKSJX6sQe4JLxjYxcOLHXLeWISNXg0vAyfPhwEhMTWbRoEYsWLSIxMZG4uLhit2nRogVvvfUWmzdv5pdffqFx48bExsZy7NgxV5YqIlcipBlYLvk1YvGCkKbuqUdEqgSLYbhmWNht27bRqlUrEhIS6Nq1KwAJCQnExMSwfft2IiMjS7SftLQ0bDYb//vf/7jpppsKrM/KyiIrK8upfcOGDbHb7dSoUaNs3oyIFG3DJ/DdOBxHYAa9CR1HuLUkEfE8+d/3Jfn+dtmRl1WrVmGz2RzBBaBbt27YbDZWrlxZon1kZ2fz/vvvY7PZaNeuXaFtpk6d6jgtZbPZaNiwYZnULyIl1HEEPPArWL3N5w26uLceEan0XBZeUlJSCAsLK7A8LCyMlJSUYrddsGABgYGB+Pv78/rrr/Pjjz8SGhpaaNtJkyZht9sdU3JycpnULyKlEN4amt9szm/52r21iEilV+rwMmXKFCwWS7HTunXrALBYLAW2Nwyj0OUX6927N4mJiaxcuZK+ffty5513kpqaWmhbPz8/atSo4TSJiBu0ucN83Pw1uOZstIgIAN6l3WDs2LEMGzas2DaNGzdm06ZNHD16tMC6Y8eOER4eXuz2AQEBNG/enObNm9OtWzeuueYaPvzwQyZNmlTackWkvET2A+9qcHIfHN4I9Tu6uyIRqaRKHV5CQ0OLPIVzsZiYGOx2O2vWrKFLF/Mc+OrVq7Hb7XTv3r1Ur2kYhlOnXBGpgPwCzQDz+zewZa7Ci4i4jMv6vERFRdG3b1/i4+NJSEggISGB+Ph4BgwY4HSlUcuWLZk3bx4A6enpTJ48mYSEBA4cOMCGDRu47777OHjwIH/6059cVaqIlJXo86eOtnwDeXnurUVEKi2X3udlzpw5REdHExsbS2xsLG3btmX27NlObXbs2IHdbgfAy8uL7du3c/vtt9OiRQsGDBjAsWPHWLFiBa1bt3ZlqSJSFprfDH42OH0Ykkp2VaGISGm57D4v7lKa68RFxAW+fQgSP4VO98LAae6uRkQ8RIW4z4uIVFHRt5uPW/8DuefcW4uIVEoKLyJSthr3goDakHkC9vzs7mpEpBJSeBGRsuXlDa2HmPNb5rq3FhGplBReRKTs5d+wbvsCOJfp3lpEpNJReBGRstewC9giIPsM7Fzs7mpEpJJReBGRsmexQJuh5rzGOhKRMqbwIiKu0eb8VUc7f4CzdvfWIiKVisKLiLhGnWgIbQG5WbD9v+6uRkQqEYUXEXENi8V5pGkRkTKi8CIirpM/1tHepZB+3K2liEjlofAiIq5TqxnUbQ9GLmz91t3ViEglofAiIq6Vf/Rls25YJyJlQ+FFRFyr9flLppNWgv2ge2sRkUpB4UVEXMtWHyK6m/NbvnFvLSJSKSi8iIjr5Y80rRvWiUgZUHgREddrdRtYvODIb3B8t7urEREPp/AiIq4XEArNepvzGmlaRK6SwouIlI/8G9Zt+RoMw721iIhHU3gRkfLR8lbw8oPjOyFls7urEREPpvAiIuXDvwa0iDXn1XFXRK6CwouIlB/HqaNvIC/PvbWIiMdSeBGR8tOiD/gGgT0Z1r4P9kPurkhEPJDCi4iUH59qENbSnP/+CZjWBjZ84t6aRMTjKLyISPmxH4KD6y48N/Jg/gQdgRGRUlF4EZHyc2IPcMll0kYunNjrlnJExDMpvIhI+QlpBpZLfu1YvCCkqXvqERGPpPAiIuXHVh8GvoHTr55+L5jLRURKSOFFRMpXxxEwYTPUOB9YfAPdW4+IeByFFxEpfzUbQKd7zflNX7i3FhHxOAovIuIebf9kPu5dBmmH3VuLiHgUhRcRcY/gxhARAxiwWcMFiEjJKbyIiPu0vdN81KkjESkFhRcRcZ/WQ8DLF45ugZQt7q5GRDyEwouIuE+1YLjm/EjTm790by0i4jEUXkTEvdreZT5u+gryct1bi4h4BIUXEXGvFn3A3wanD8P+X9xdjYh4AIUXEXEvbz+z7wuo466IlIjCi4i4X9th5uPW7yA7w721iEiFp/AiIu7XsCvUjIDs07BjoburEZEKzqXh5eTJk8TFxWGz2bDZbMTFxXHq1KkSb3///fdjsViYNm2ay2oUkQrAar2o466uOhKR4rk0vAwfPpzExEQWLVrEokWLSExMJC4urkTbfvvtt6xevZp69eq5skQRqSjyw8vu/0H6cffWIiIVmsvCy7Zt21i0aBEffPABMTExxMTEMGPGDBYsWMCOHTuK3fbQoUOMHTuWOXPm4OPjU2zbrKws0tLSnCYR8UCh10C9jmDkwpa57q5GRCowl4WXVatWYbPZ6Nq1q2NZt27dsNlsrFy5ssjt8vLyiIuLY+LEibRu3fqyrzN16lTHaSmbzUbDhg3LpH4RcQPHqSNddSQiRXNZeElJSSEsLKzA8rCwMFJSUorc7sUXX8Tb25tx48aV6HUmTZqE3W53TMnJyVdcs4i4WZvbweIFh9bD8V3urkZEKqhSh5cpU6ZgsViKndatWweAxWIpsL1hGIUuB1i/fj1vvPEGH330UZFtLuXn50eNGjWcJhHxUIG1oflN5rw67opIEbxLu8HYsWMZNmxYsW0aN27Mpk2bOHr0aIF1x44dIzw8vNDtVqxYQWpqKhEREY5lubm5PProo0ybNo39+/eXtlwR8TRt74JdP5injnpPhhL+ISMiVUepw0toaCihoaGXbRcTE4PdbmfNmjV06dIFgNWrV2O32+nevXuh28TFxXHzzTc7LevTpw9xcXHce++9pS1VRDxRZH/wDYRTByB5NUR0c3dFIlLBuKzPS1RUFH379iU+Pp6EhAQSEhKIj49nwIABREZGOtq1bNmSefPmAVCrVi3atGnjNPn4+FCnTh2nbUSkEvOtDlGDzHl13BWRQrj0Pi9z5swhOjqa2NhYYmNjadu2LbNnz3Zqs2PHDux2uyvLEBFP0+78VUebv4bdP4H9kHvrEZEKxWIYhuHuIspSWloaNpsNu92uzrsiniovF15qAmfP/2FjscLAN6DjCPfWJSIuU5rvb41tJCIVz+mUC8EFwMiD+RN0BEZEAIUXEamITuwpuMzIhRN7y78WEalwFF5EpOIJaWaeKrqYxQtCmrqnHhGpUBReRKTisdU3+7hw0T1eBk4zl4tIlafwIiIVU8cRMOZXsJ6/HZVhqM+LiAAKLyJSkdVpDeFtzPn542BaG9jwiXtrEhG3U3gRkYrLfgiO/Hbhua46EhEUXkSkIjuxB7jkVlS66kikylN4EZGKS1cdiUghFF5EpOJyXHV00a+qG/+mq45EqjiFFxGp2DqOgIe3QL0O5vPsM+6tR0TcTuFFRCo+W33oMcGc3zgHcs+5tRwRcS+FFxHxDJH9IaA2nEmBnYvdXY2IuJHCi4h4Bm9faD/cnF//kVtLERH3UngREc/R8R7zcff/4FSye2sREbdReBERz1GrGTTpBRiwcba7qxERN1F4ERHPkn/0ZcNsyM1xby0i4hYKLyLiWaIGQrUQOH3YPH0kIlWOwouIeBZvP3XcFaniFF5ExPPknzratRjSDru3FhEpdwovIuJ5areARj3MUaY3furuakSknCm8iIhncnTc/QTyct1bi4iUK4UXEfFMrQaBf02wJ8PK6WA/5O6KRKScKLyIiGfyqQZ125nz/5sC09qYR2FEpNJTeBERz2Q/BPuWX3hu5MH8CToCI1IFKLyIiGc6sQcwnJcZuXBir1vKEZHyo/AiIp4ppBlYLv0VZoWQpm4pR0TKj8KLiHgmW30Y+AZOv8ZaDTKXi0ilpvAiIp6r4wh4eAt0H2c+T14NOdnurUlEXE7hRUQ8m60+3Pg3CKwDp4/A5q/cXZGIuJjCi4h4Pm8/6PaAOb9yOuTlubceEXEphRcRqRw63wu+QXBsO+z+0d3ViIgLKbyISOXgb4POI835X99wayki4loKLyJSeXR9AKw+cOBXOLjO3dWIiIsovIhI5WGrD23vNOd19EWk0lJ4EZHKpftfzcdt8+GPPe6tRURcQuFFRCqXsCi4pg9gwMo33V2NiLiAwouIVD49xpuPiZ/BmVT31iIiZU7hRUQqn0bdoX4nyM2CNe+7uxoRKWMuDS8nT54kLi4Om82GzWYjLi6OU6dOFbvNyJEjsVgsTlO3bt1cWaaIVDYWy4WjL2tmQNYZ99YjImXKpeFl+PDhJCYmsmjRIhYtWkRiYiJxcXGX3a5v374cOXLEMS1cuNCVZYpIZdRygDnC9NlTsPFTd1cjImXI21U73rZtG4sWLSIhIYGuXbsCMGPGDGJiYtixYweRkZFFbuvn50edOnVcVZqIVAVWL4gZC/99BFa9BdeOBi8fd1clImXAZUdeVq1ahc1mcwQXgG7dumGz2Vi5cmWx2y5dupSwsDBatGhBfHw8qalFd7jLysoiLS3NaRIRAaD9cKgeCvZk+P1bd1cjImXEZeElJSWFsLCwAsvDwsJISUkpcrt+/foxZ84clixZwquvvsratWu58cYbycrKKrT91KlTHX1qbDYbDRs2LLP3ICIezqcadB1jzv/6BhiGe+sRkTJR6vAyZcqUAh1qL53WrTNvy22xWApsbxhGocvz3XXXXdx66620adOGgQMH8v3337Nz507++9//Ftp+0qRJ2O12x5ScnFzatyQildm1o8GnOhzdDHt/dnc1IlIGSt3nZezYsQwbNqzYNo0bN2bTpk0cPXq0wLpjx44RHh5e4terW7cujRo1YteuXYWu9/Pzw8/Pr8T7E5EqpnoIdBwBq981j740u9HdFYnIVSp1eAkNDSU0NPSy7WJiYrDb7axZs4YuXboAsHr1aux2O927dy/x6/3xxx8kJydTt27d0pYqImLq9qB5yfTepXA4Eeq1d3NBInI1XNbnJSoqir59+xIfH09CQgIJCQnEx8czYMAApyuNWrZsybx58wA4c+YMjz32GKtWrWL//v0sXbqUgQMHEhoaypAhQ1xVqohUdsGNoM1Qc15DBoh4PJfe52XOnDlER0cTGxtLbGwsbdu2Zfbs2U5tduzYgd1uB8DLy4vNmzczePBgWrRowT333EOLFi1YtWoVQUFBrixVRCq77uPMx9/nwcn9bi1FRK6OxTAqV/f7tLQ0bDYbdrudGjVqFNkuNzeXc+fOlWNlIp7Jx8cHLy8vd5dRNj65zey02+V+6P+Su6sRkYuU9PsbXHiTuorKMAxSUlIuO0yBiFxQs2ZN6tSpU+yVgh6hx3gzvGycDTc8aXbmFRGPU+XCS35wCQsLo3r16p7/y1jEhQzDICMjw3GjSI/vON/0BqgTDSmbYcVr0CIWQpqBrb67KxORUqhS4SU3N9cRXGrVquXuckQ8QrVq1QBITU0lLCzMs08hWSzQYwLMHQ2r3jQnixUGvmFeTi0iHsGlHXYrmvw+LtWrV3dzJSKeJf//TKXoJ9bgWufnRh7MnwD2Q24pR0RKr0qFl3w6VSRSOpXq/8ypAwWXGblwYm/51yIiV6RKhhcRqcJCmgGXhDGLF4Q0dUs5IlJ6Ci+VxNKlS7FYLGV6FVXjxo2ZNm1ame1PpEKw1YdB03H69RfzkDrtingQhRfxWMuWLaNTp074+/vTtGlT3n33XXeXJJ6i4wh4eAu06Gc+37EQzp11b00iUmIKL1fhiD2TlXuOc8Se6e5SKq3s7OxCl+/bt4/+/fvTs2dPNm7cyOTJkxk3bhxz584t5wrFY9nqw5B3ITAc/tgNv7zu7opEpIQUXq7QF2uT6PHCEobPWE2PF5bwxdokl75eVlYW48aNIywsDH9/f6677jrWrl1boN2vv/5Ku3bt8Pf3p2vXrmzevNmx7sCBAwwcOJDg4GACAgJo3bo1CxcuLHENr732GtHR0QQEBNCwYUMefPBBzpw5A0B6ejo1atTg66+/dtpm/vz5BAQEcPr0aQAOHTrEXXfdRXBwMLVq1WLw4MHs37/f0X7kyJHcdtttTJ06lXr16tGiRYtCa3n33XeJiIhg2rRpREVFcd999zFq1CheeeWVEr8fEarVhH4vmvO/vAbHdrq1HBEpGYWXK3DEnsmkbzaTd35ghTwDJn+zxaVHYB5//HHmzp3Lxx9/zIYNG2jevDl9+vThxIkTTu0mTpzIK6+8wtq1awkLC2PQoEGOy1sfeughsrKyWL58OZs3b+bFF18kMDCwxDVYrVamT5/Oli1b+Pjjj1myZAmPP/44AAEBAQwbNoxZs2Y5bTNr1izuuOMOgoKCyMjIoHfv3gQGBrJ8+XJ++eUXAgMD6du3r9MRlp9++olt27bx448/smDBgkJrWbVqFbGxsU7L+vTpw7p16yrH5bxSflrdBtfEQm42LJgAeXnurkhELkPh5QrsO57uCC75cg2D/cczXPJ66enpvPPOO7z88sv069ePVq1aMWPGDKpVq8aHH37o1PbZZ5/llltuITo6mo8//pijR486Ru1OSkqiR48eREdH07RpUwYMGECvXr1KXMeECRPo3bs3TZo04cYbb+S5557jyy+/dKy/7777WLx4MYcPHwbg+PHjLFiwgFGjRgHw+eefY7Va+eCDD4iOjiYqKopZs2aRlJTE0qVLHfsJCAjggw8+oHXr1rRp06bQWlJSUggPD3daFh4eTk5ODsePHy/xexLBYoH+r4BPdTjwKyTOcXdFInIZCi9XoEloANZLrrT0slhoHOqam9/t2bOHc+fO0aNHD8cyHx8funTpwrZt25zaxsTEOOZDQkKIjIx0tBk3bhz//Oc/6dGjB88++yybNm0qVR0///wzt9xyC/Xr1ycoKIgRI0bwxx9/kJ6eDkCXLl1o3bo1n3zyCQCzZ88mIiLCEZDWr1/P7t27CQoKIjAwkMDAQEJCQjh79ix79uxxvE50dDS+vr6XrefSe4/kjzFaqe5JIuUjuBH0nmzO//A3OHPMvfWISLEUXq5AXVs1pg6Nxuv8l6SXxcLzQ9tQ11bNJa9X1JeyYRgl+qLOb3Pfffexd+9e4uLi2Lx5M507d+bNN98sUQ0HDhygf//+tGnThrlz57J+/Xr+9a9/Ac53Xb3vvvscp45mzZrFvffe63j9vLw8OnXqRGJiotO0c+dOhg8f7thHQEDAZeupU6cOKSkpTstSU1Px9vbW0A9yZbo+YI57dPYULJ7s7mpEpBgKL1formsj+OXJ3vw7vhu/PNmbu66NcNlrNW/eHF9fX3755RfHsnPnzrFu3TqioqKc2iYkJDjmT548yc6dO2nZsqVjWcOGDRkzZgzffPMNjz76KDNmzChRDevWrSMnJ4dXX32Vbt260aJFC8fpoYv9+c9/JikpienTp/P7779zzz33ONZ17NiRXbt2ERYWRvPmzZ0mm81W4p8HmEeYfvzxR6dlP/zwA507d8bHx6dU+xIBwMvbHOPIYoXNX8KeJe6uSESKoPByFeraqhHTrJbLjrjkCwgI4IEHHmDixIksWrSIrVu3Eh8fT0ZGBqNHj3Zq+49//IOffvqJLVu2MHLkSEJDQ7ntttsAs8/K4sWL2bdvHxs2bGDJkiUFwk9RmjVrRk5ODm+++SZ79+5l9uzZhd5XJTg4mKFDhzJx4kRiY2Np0KCBY93dd99NaGgogwcPZsWKFezbt49ly5Yxfvx4Dh48WKqfyZgxYzhw4ACPPPII27ZtY+bMmXz44Yc89thjpdqPiJP6naDLX8z5BQ/DOd0GQaQiUnjxEC+88AK33347cXFxdOzYkd27d7N48WKCg4MLtBs/fjydOnXiyJEjfPfdd47+I7m5uTz00ENERUXRt29fIiMjefvtt0v0+u3bt+e1117jxRdfpE2bNsyZM4epU6cW2nb06NFkZ2c7Ourmq169OsuXLyciIoKhQ4cSFRXFqFGjyMzMpEaNGqX6eTRp0oSFCxeydOlS2rdvz3PPPcf06dO5/fbbS7UfkQJu/BvUqA8n98Oyl9xdjYgUwmLkd6ioJNLS0rDZbNjt9gJfiGfPnmXfvn00adIEf39/N1VY+c2ZM4fx48dz+PDhEnW8lYqvyv3f2f5f+Hw4WL3h/uUQ3trdFYlUesV9f19KR16kzGRkZPD7778zdepU7r//fgUX8Vwtb4WWAyAvB+ZP0L1fRCoYhRcpMy+99BLt27cnPDycSZMmubsckavT7yXwDYKDa2D9rMu3F5Fyo/AiZWbKlCmcO3eOn376qVR37hWpkGz14aanzfn//R1OpxTfXkTKjcKLiEhRrr3PvAIpyw7fP+HuakTkPIUXEZGiWL1gwDSweMHWb2HnYndXJCIovIiIFK9uW4h50Jz/72OQne7eekRE4UVE5LJumAS2CLAnwc/Pu7sakSpP4UVE5HJ8A2DAa+Z8wjtw5Df31iNSxSm8iIiUxDW3QOuhYOTC/PGQl+vuikSqLIWXSmLp0qVYLBZOnTpVZvts3Lgx06ZNK7P9iXi8vi+Anw0Ob4RFT4D9kLsrEqmSFF7EIx05coThw4cTGRmJ1WplwoQJ7i5JqoKgcIjsa86vmQGvt4YNn7i3JpEqSOHlatgPwb7l+uvLhbKzswtdnpWVRe3atXnqqado165dOVclVZb9EGz68qIFBnw3Xr8DRMqZwsuV2vAJTGsDHw80H13811dWVhbjxo0jLCwMf39/rrvuOtauXVug3a+//kq7du3w9/ena9eubN682bHuwIEDDBw4kODgYAICAmjdujULFy4scQ2vvfYa0dHRBAQE0LBhQx588EHOnDkDQHp6OjVq1ODrr7922mb+/PkEBARw+vRpAA4dOsRdd91FcHAwtWrVYvDgwezfv9/RfuTIkdx2221MnTqVevXq0aJFi0Jrady4MW+88QYjRozAZrOV+D2IXJUTe4BLx7LN0/1fRMqZwsuVsB8yO+wZ5wdrM/LMwdtc+NfX448/zty5c/n444/ZsGEDzZs3p0+fPpw4ccKp3cSJE3nllVdYu3YtYWFhDBo0iHPnzgHw0EMPkZWVxfLly9m8eTMvvvhiqW7jb7VamT59Olu2bOHjjz9myZIlPP744wAEBAQwbNgwZs1yHgNm1qxZ3HHHHQQFBZGRkUHv3r0JDAxk+fLl/PLLLwQGBtK3b1+nIyw//fQT27Zt48cff2TBggVX+iMTKXshzcBSyK/NHybDrv+Vfz0iVZTCy5U4sedCcMln5MKJvS55ufT0dN555x1efvll+vXrR6tWrZgxYwbVqlXjww8/dGr77LPPcssttxAdHc3HH3/M0aNHmTdvHgBJSUn06NGD6OhomjZtyoABA+jVq1eJ65gwYQK9e/emSZMm3HjjjTz33HN8+eWFQ+j33Xcfixcv5vDhwwAcP36cBQsWMGrUKAA+//xzrFYrH3zwAdHR0URFRTFr1iySkpJYunSpYz8BAQF88MEHtG7dmjZt2lzpj02k7Nnqw8A3zDvuAmCF2i3hXCb8+y5I/Myt5YlUFQovV6Kwv74sXhDS1CUvt2fPHs6dO0ePHj0cy3x8fOjSpQvbtm1zahsTE3OhzJAQIiMjHW3GjRvHP//5T3r06MGzzz7Lpk2bSlXHzz//zC233EL9+vUJCgpixIgR/PHHH6Snm3cc7dKlC61bt+aTT8xTaLNnzyYiIsIRkNavX8/u3bsJCgoiMDCQwMBAQkJCOHv2LHv27HG8TnR0NL6+vqWqTaTcdBwBEzbDPQvg4S1w/wpoexfk5cC3D8DyV8C49NSSiJQlhZcrcelfXxYvGDjNXO4CxvlfhBaLpcDyS5cVJr/Nfffdx969e4mLi2Pz5s107tyZN998s0Q1HDhwgP79+9OmTRvmzp3L+vXr+de//gXgOC2V/xr5p45mzZrFvffe63j9vLw8OnXqRGJiotO0c+dOhg8f7thHQEBAiWoScRtbfWjS03z09oXb3oUeE8x1S56DhY/pPjAiLqTwcqUu/utrwmbzuYs0b94cX19ffvnlF8eyc+fOsW7dOqKiopzaJiQkOOZPnjzJzp07admypWNZw4YNGTNmDN988w2PPvooM2bMKFEN69atIycnh1dffZVu3brRokULx+mhi/35z38mKSmJ6dOn8/vvv3PPPfc41nXs2JFdu3YRFhZG8+bNnSZ1uhWPZrXCLX+Hfi8BFlj7AXw5wjydJCJlTuHlalz815cLBQQE8MADDzBx4kQWLVrE1q1biY+PJyMjg9GjRzu1/cc//sFPP/3Eli1bGDlyJKGhodx2222A2Wdl8eLF7Nu3jw0bNrBkyZIC4acozZo1IycnhzfffJO9e/cye/Zs3n333QLtgoODGTp0KBMnTiQ2NpYGDRo41t19992EhoYyePBgVqxYwb59+1i2bBnjx4/n4MGDpf655B+5OXPmDMeOHSMxMZGtW7eWej8iZabr/fCnj8DLD7YvgE8GQ8aJy24mIqWj8OIhXnjhBW6//Xbi4uLo2LEju3fvZvHixQQHBxdoN378eDp16sSRI0f47rvvHP1HcnNzeeihh4iKiqJv375ERkby9ttvl+j127dvz2uvvcaLL75ImzZtmDNnDlOnTi207ejRo8nOznZ01M1XvXp1li9fTkREBEOHDiUqKopRo0aRmZlJjRo1Sv0z6dChAx06dGD9+vV89tlndOjQgf79+5d6PyJlqvVtEDfPvBNv8mqY2QdOJbm7KpFKxWIYlatnWVpaGjabDbvdXuAL8ezZs+zbt48mTZrg7+/vpgorvzlz5jB+/HgOHz6sjreVhP7vXIGjW2HOHZB2CALrwJ+/hjrR7q5KpMIq7vv7Ui498nLy5Eni4uKw2WzYbDbi4uJKNPbOtm3bGDRoEDabjaCgILp160ZSkv5yqegyMjL4/fffmTp1Kvfff7+Ci1Rt4a1g9A9QOwrOpMCs/rB3mburEqkUXBpehg8fTmJiIosWLWLRokUkJiYSFxdX7DZ79uzhuuuuo2XLlixdupTffvuNp59+Wn/teYCXXnqJ9u3bEx4ezqRJk9xdjoj72RrAqEXQqAdkpcGnt8Pmry+/nYgUy2WnjbZt20arVq1ISEiga9eugHklTExMDNu3bycyMrLQ7YYNG4aPjw+zZ8++otfVaSORsqf/O1fp3FmYdz9s/dZ8Hvt/0H2sW0sSqWgqxGmjVatWYbPZHMEFoFu3bthsNlauXFnoNnl5efz3v/+lRYsW9OnTh7CwMLp27cq3335b5OtkZWWRlpbmNImIVCg+/nDHLOg6xnz+w1OwaDLk5RW/nYgUymXhJSUlhbCwsALLw8LCSElJKXSb1NRUzpw5wwsvvEDfvn354YcfGDJkCEOHDmXZssLPFU+dOtXRp8Zms9GwYcMyfR8iImXCaoW+L8At/zCfJ/wL5o6GnCz31iXigUodXqZMmYLFYil2WrduHVDwjrBQ/F1h887/FTJ48GAefvhh2rdvz5NPPsmAAQMKvacIwKRJk7Db7Y4pOTm5tG9JRKR8WCzQYzwMeR+s3vD7N2Y/mLN2d1cm4lG8S7vB2LFjGTZsWLFtGjduzKZNmzh69GiBdceOHSM8PLzQ7UJDQ/H29qZVq1ZOy6OiopzuLnsxPz8//Pz8Sli9iEgF0O4uCKwNX8TB/hUwsx8Mfguyz5hjp7n4xpcinq7U4SU0NJTQ0NDLtouJicFut7NmzRq6dOkCwOrVq7Hb7XTv3r3QbXx9fbn22mvZsWOH0/KdO3fSqFGj0pYqIlJxNbsR7l0Ic/4Eqb/DjN7mcovVHDvNhUOOiHg6l/V5yb+La3x8PAkJCSQkJBAfH8+AAQOcrjRq2bIl8+bNczyfOHEiX3zxBTNmzGD37t289dZbzJ8/nwcffNBVpYqIuEfddnDnJVdWGnnw3Xg4pVPgIkVx6X1e5syZQ3R0NLGxscTGxtK2bdsCl0Dv2LEDu/3C+d4hQ4bw7rvv8tJLLxEdHc0HH3zA3Llzue6661xZqlzkhhtuYMKECSVu/9FHH1GzZk3H8ylTptC+ffurqmH//v1YLBYSExOLbNO4cWOmTZt2Va9zpSwWS7FXwZX2ZyhVWG5hHXbz4MNbYPV7kHW63EsSqehKfdqoNEJCQvj000+LbVPYbWZGjRpVYFwcESk7mzdvZuzYsaxZs4aQkBDuv/9+nn766SI704sLhTQzTxUZl1w2ffoIfP84LPkndIiDrn+B4MZgPwQn9qhvjFRpLg0vIuI+586dw8fHp8DytLQ0brnlFnr37s3atWvZuXMnI0eOJCAggEcffdQNlVZxtvpmH5f5E8DIBYsX9HvRXLf6Pfhjl3lZ9ep3ILwNpGwGDPWNkSqtyo8qbRgGGdk5bplKc3PjG264gb/+9a9MmDCB4OBgwsPDef/990lPT+fee+8lKCiIZs2a8f333zttt2zZMrp06YKfnx9169blySefJCcnx7E+PT2dESNGEBgYSN26dXn11VcLvHZ2djaPP/449evXJyAggK5du7J06dJS/ZxnzZpFVFQU/v7+tGzZssBo1mvWrKFDhw74+/vTuXNnNm7cWKL9ZmRkMGrUKIKCgoiIiOD99993Wn/o0CHuuusugoODqVWrFoMHD2b//v2O9WvXruWWW24hNDQUm83G9ddfz4YNG5z2sWvXLnr16oW/vz+tWrXixx9/LNV7B/j000/p3LkzQUFB1KlTh+HDh5OamgqY/wabN2/OK6+84rTNli1bsFqt7NmzBwC73c5f/vIXwsLCqFGjBjfeeCO//fabo33+6bqZM2fStGlT/Pz8Cv03NmfOHM6ePctHH31EmzZtGDp0KJMnT+a1114r1b9JKUMdR8CEzXDPAvOxS7w5PbQG7v7a7Nxr5EHKJuD8Z2TkwXfjYM9S3exOqpwqf+Ql81wurZ5Z7JbX3vqPPlT3LflH8PHHH/P444+zZs0avvjiCx544AG+/fZbhgwZwuTJk3n99deJi4sjKSmJ6tWrc+jQIfr378/IkSP55JNP2L59O/Hx8fj7+zNlyhTA7CD9888/M2/ePOrUqcPkyZNZv369U5+Ve++9l/379/P5559Tr1495s2bR9++fdm8eTPXXHPNZeueMWMGzz77LG+99RYdOnRg48aNxMfHExAQwD333EN6ejoDBgzgxhtv5NNPP2Xfvn2MHz++RD+TV199leeee47Jkyfz9ddf88ADD9CrVy9atmxJRkYGvXv3pmfPnixfvhxvb2/++c9/0rdvXzZt2oSvry+nT5/mnnvuYfr06Y799e/fn127dhEUFEReXh5Dhw4lNDSUhIQE0tLSrqgvS3Z2Ns899xyRkZGkpqby8MMPM3LkSBYuXIjFYmHUqFHMmjWLxx57zLHNzJkz6dmzJ82aNcMwDG699VZCQkJYuHAhNpuN9957j5tuuomdO3cSEhICwO7du/nyyy+ZO3cuXl5ehdayatUqrr/+eqdbDPTp04dJkyaxf/9+mjRpUur3J2XAVr/gaSCrFa65xZw2fgr/eeiSjQyYPRiqBUPDbtAoBiK6mx2BvX11ikkqrSofXjxJu3bt+Nvf/gaYN+d74YUXCA0NJT4+HoBnnnmGd955h02bNtGtWzfefvttGjZsyFtvvYXFYqFly5YcPnyYJ554gmeeeYaMjAw+/PBDPvnkE2655RbADEgNGjRwvOaePXv497//zcGDB6lXrx4Ajz32GIsWLWLWrFk8//zzl637ueee49VXX2Xo0KEANGnShK1bt/Lee+9xzz33MGfOHHJzc5k5cybVq1endevWHDx4kAceeOCy++7fv7/jSrQnnniC119/naVLl9KyZUs+//xzrFYrH3zwgaMvx6xZs6hZsyZLly4lNjaWG2+80Wl/7733HsHBwSxbtowBAwbwv//9j23btrF//37Hz+X555+nX79+l63tYhf34WratCnTp0+nS5cunDlzhsDAQO69916eeeYZx60Fzp07x6effsrLL78MwM8//8zmzZtJTU11hI5XXnmFb7/9lq+//pq//OUvgBmSZs+eTe3atYusJSUlhcaNGzsty7/3UkpKisJLRdW0d+F9Y7yrQeZJ2Pm9OeUvq9kAju8638gCff4Puj1o3igvX0nDjUKQVDBVPrxU8/Fi6z/6uO21S6Nt27aOeS8vL2rVqkV0dLRjWf4XUP7piG3bthETE+PUCbNHjx6cOXOGgwcPcvLkSbKzs4mJiXGsDwkJcbqUfcOGDRiGQYsWLZxqycrKolatWpet+dixYyQnJzN69GhHyALIycnBZrM56mzXrh3Vq1d3rL+4puJc/DOxWCzUqVPH8f7Xr1/P7t27CQoKctrm7NmzjlMxqampPPPMMyxZsoSjR4+Sm5tLRkYGSUlJjtoiIiKcAl1Ja7vYxo0bmTJlComJiZw4ccJxN+mkpCRatWpF3bp1ufXWW5k5cyZdunRhwYIFnD17lj/96U+O93LmzJkCP/PMzEzHewFo1KhRscEl36Udc/NPF6nDbgVWWN+YgdOg3f+DI5sgaSUcWAVJqyDzxEXBBcCAxZNhyXNQszEEN4JzmbBvubkOC/R8BDr8GaqFgL/tQsjZ8AnMH2+GpivpZ6PgIy5Q5cOLxWIp1akbd7q086XFYnFalv/Fk//FWNhQDBd/SZWkf0NeXh5eXl6sX7++wGmIwMDAEm0P5qmjiwfpBBz7u5p+FoX9TPJfMy8vj06dOjFnzpwC2+V/wY8cOZJjx44xbdo0GjVqhJ+fHzExMWRnZxdZW2m/4NPT0x23C/j000+pXbs2SUlJ9OnTx/E6APfddx9xcXG8/vrrzJo1i7vuussR6PLy8qhbt26hfY0uvkw9ICDgsvXUqVOnwPhi+YGvqLtfSwXRcQQ0uwlO7IWQphfCQINO5tT9r2b/l98+K+QUE2ZgObbNnJwYsOJVcwIzGFUPAb8g87Uczc73szm+C2rUB79A8D0/+V3y6BsIm76ABRMuH3wUcKSUPONbW65Iq1atmDt3rlOIWblyJUFBQdSvX5/g4GB8fHxISEggIiICgJMnT7Jz506uv/56ADp06EBubi6pqan07Nmz1DWEh4dTv3599u7dy913311knbNnzyYzM5Nq1aoBkJCQcCVv2UnHjh354osvHB1cC7NixQrefvtt+vfvD0BycjLHjx93qi0pKYnDhw87TputWrWqVHVs376d48eP88ILLzgGDs0f/+ti/fv3JyAggHfeeYfvv/+e5cuXO72XlJQUvL29C5zyKa2YmBgmT55MdnY2vr6+APzwww/Uq1fvqvct5aCwvjEXs1qLOMVkhZELICcTdv8ECW8X3Nbbzxwo0siF9GPmVIABK6eXvm4jD777KyT+GwJqgW+QGXRO7oddP+I4AtTtAWhzB1SrCf41zaNAXqX4qlIQqhIUXiqxBx98kGnTpvHXv/6VsWPHsmPHDp599lkeeeQRrFYrgYGBjB49mokTJ1KrVi3Cw8N56qmnsFovXITWokUL7r77bkaMGMGrr75Khw4dOH78OEuWLCE6OtrxpV+cKVOmMG7cOGrUqEG/fv3Iyspi3bp1nDx5kkceeYThw4fz1FNPMXr0aP72t7+xf//+AlfeXIm7776bl19+mcGDB/OPf/yDBg0akJSUxDfffMPEiRNp0KABzZs3Z/bs2XTu3Jm0tDQmTpzoCFAAN998M5GRkY73n5aWxlNPPVWqOiIiIvD19eXNN99kzJgxbNmyheeee65AOy8vL0aOHMmkSZNo3ry50+mpm2++mZiYGG677TZefPFFIiMjOXz4MAsXLuS2226jc+fOJa5n+PDh/P3vf2fkyJFMnjyZXbt28fzzz/PMM8/otFFlUdQppsY9zPW1o2D1u87hxuIFf91oHnHJOGGeejq20xz5mouPQFqg9RBzWdYZczym/Mf8+ZzMomtLWllM4YYZqi4NVr6BZpCpVtMMM47588/z5w+vh9Xv4whCA9+ATvdc/uelwONxFF4qsfr167Nw4UImTpxIu3btCAkJcQSEfC+//DJnzpxh0KBBBAUF8eijjzrd8RjMTq7//Oc/efTRRzl06BC1atUiJiamRMEFzNMh1atX5+WXX+bxxx8nICCA6Ohox1U7gYGBzJ8/nzFjxtChQwdatWrFiy++yO23335V77969eosX76cJ554gqFDh3L69Gnq16/PTTfd5DgSM3PmTP7yl7/QoUMHIiIieP75552u+LFarcybN4/Ro0fTpUsXGjduzPTp0+nbt2+J66hduzYfffQRkydPZvr06XTs2JFXXnmFQYMGFWg7evRonn/++QI3abRYLCxcuJCnnnqKUaNGcezYMerUqUOvXr1KfarHZrPx448/8tBDD9G5c2eCg4N55JFHeOSRR0q1H6ngijrFBEWHm/w2+Ud36kTDufSC7S7X5yU3B/7YDW93wzn4WKHP8+aRlOwzcHQrbP6y4PYBteHcWcg+f3fh/GCUdrAUPwAD5o+DZS9AUF0ICDMHwwwIg8Aw8zUCw+DASvj5eXTvHM9iMSrZjR3S0tKw2WzY7fYCpwrOnj3Lvn37aNKkCf7+/m6qUKRov/76KzfccAMHDx6sUP1P9H+nkrIfKjzcXGm7S234pPjgYz8E09oUPAI0YbP5Ork5cNYOZ0+ZU+b5x7P2C/P5jycPwJHEktdWnDZ/gtqRULMh2BqYU4364FXwpo86alN2ivv+vpSOvIhUAFlZWSQnJ/P0009z5513VqjgIpXY5frPlLbdpYo7+pO/3+KOAHl5m/1jAi5/ZWOhQQgrDJsDGHAm1ezDcyYV0lPhzDE4dQDSDhXc15avCnkBi3kE5+JAk3YYNn+N4zTVoOlFH7VRyClTCi8iFcC///1vRo8eTfv27QsMXiri0S4XfC4XcErzOoUFoZbFnN4uNPBYzPvhZNnNkb3tB80pNwtOHzan5NWF7MwwOyRv+tI83VarOYReA7WugV0/lOyqKykxnTYSkcvS/x3xGKU9xXW5U1tgXn6ecfx8mDk/Ja+GbfOvsEgr3L8M6rZ1XlzFj87otJGIiFRNpT3FVZIjP1ar2bk3MMy8nw6AfShs/2/B01Q3PQ3px80BNY/vMk9NXXpXZPLgvZ4QGG4Otlmnjdl3Z8MnqONwySi8iIhI1XYlfXqKOk11aeA4sQ+md8D5qisAC5w5ak57fnJelX8zQKsPRPYzLwUXJwovIiIiV6IkR21CmpgdeS8NOW1uNy8VP7rZ7BOz4/tLNjTg2zGABcKioGFXiIiBiK5Qs5HzGFVQ5U45KbyIiIhcqZIctSkq5DS81pyu6QM7FxfsOFyzIZxKgtSt5rR+lrkqsI4ZYhp2Mx+PbIL/PlKlOgQrvIiIiLhacSGnuFNQZ1LNzsFJCebj4UQ4kwJb/2NOlzLy4LvxZliqxEdgFF5ERETcraijM4FhEDXQnMAcXPPQBkhOgKTVcOAXyE6/ZGd5sHAidL4XGvcEn8p3haD18k2kqrnhhhsct+4viY8++shpZOMpU6bQvn37q6ph//79WCwWEhMTi2zTuHFjpk2bdlWvc6UsFgvffvttketL+zMUEcFWH5r0LP6IiU81c4yqno/C3V/Cg6uBQsYk2/FfmHMHvNQEPhsG62aZN9UDs3/MvuXmo4fSkReRKubs2bOMGTOG9evXs23bNgYMGFBsEBORCqxmQ+cOwVihSzzkZpv9aE4fhp3fmxNAjQYXjRF1mbsCV2AKLyKV1Llz5/DxKTgWS25uLtWqVWPcuHHMnTvXDZWJSJkq6pSTYcDRLbBzEez8AQ6uvWRwS8O8JDusFTQo+cj0FYFOGxmGeb7QHVMpbm58ww038Ne//pUJEyYQHBxMeHg477//Punp6dx7770EBQXRrFkzvv/e+XK7ZcuW0aVLF/z8/Khbty5PPvkkOTk5jvXp6emMGDGCwMBA6taty6uvvlrgtbOzs3n88cepX78+AQEBdO3alaVLl5bqxzxr1iyioqLw9/enZcuWvP2285D3a9asoUOHDvj7+9O5c2c2btxYov1mZGQwatQogoKCiIiI4P3333daf+jQIe666y6Cg4OpVasWgwcPZv/+/Y71a9eu5ZZbbiE0NBSbzcb111/Phg0bnPaxa9cuevXqhb+/P61ateLHH38s1XsH+PTTT+ncuTNBQUHUqVOH4cOHk5qaCoBhGDRv3pxXXnnFaZstW7ZgtVrZs2cPAHa7nb/85S+EhYVRo0YNbrzxRn777TdH+/zTdTNnzqRp06b4+flR2A20AwICeOedd4iPj6dOnTqlfi8iUgEVdsrJYjGHKug1Ee77Ee4qbOgRAz64GT4aABtmm4NeQoU/taQjL+cy4Pl67nntyYfBN6DEzT/++GMef/xx1qxZwxdffMEDDzzAt99+y5AhQ5g8eTKvv/46cXFxJCUlUb16dQ4dOkT//v0ZOXIkn3zyCdu3byc+Ph5/f3+mTJkCwMSJE/n555+ZN28ederUYfLkyaxfv96pz8q9997L/v37+fzzz6lXrx7z5s2jb9++bN68mWuuueaydc+YMYNnn32Wt956iw4dOrBx40bi4+MJCAjgnnvuIT09nQEDBnDjjTfy6aefsm/fPsaPH1+in8mrr77Kc889x+TJk/n666954IEH6NWrFy1btiQjI4PevXvTs2dPli9fjre3N//85z/p27cvmzZtwtfXl9OnT3PPPfcwffp0x/769+/Prl27CAoKIi8vj6FDhxIaGkpCQgJpaWlX1JclOzub5557jsjISFJTU3n44YcZOXIkCxcuxGKxMGrUKGbNmsVjjz3m2GbmzJn07NmTZs2aYRgGt956KyEhISxcuBCbzcZ7773HTTfdxM6dOwkJCQFg9+7dfPnll8ydOxcvL69S1ykilVi9jual1AXu+GvA/hXmtPAx874yhxOp0Hf7NSoZu91uAIbdbi+wLjMz09i6dauRmZl5YWHWGcN4toZ7pqwzJX5f119/vXHdddc5nufk5BgBAQFGXFycY9mRI0cMwFi1apVhGIYxefJkIzIy0sjLy3O0+de//mUEBgYaubm5xunTpw1fX1/j888/d6z/448/jGrVqhnjx483DMMwdu/ebVgsFuPQoUNO9dx0003GpEmTDMMwjFmzZhk2m82x7tlnnzXatWvneN6wYUPjs88+c9r+ueeeM2JiYgzDMIz33nvPCAkJMdLT0x3r33nnHQMwNm7cWOTPpFGjRsaf//xnx/O8vDwjLCzMeOeddwzDMIwPP/ywwPvPysoyqlWrZixevLjQfebk5BhBQUHG/PnzDcMwjMWLFxteXl5GcnKyo833339vAMa8efOKrO366693/AwLs2bNGgMwTp8+bRiGYRw+fNjw8vIyVq9ebRiGYWRnZxu1a9c2PvroI8MwDOOnn34yatSoYZw9e9ZpP82aNTPee+89wzDMn7uPj4+Rmppa5Ote6p577jEGDx582XaF/t8REc+z/mPDmBJsfgdNCTafnzxgGMtfMYw3ry3i+6qmYZw66PLSivv+vpSOvPhUN4+AuOu1S6Ft2wuDeHl5eVGrVi2io6Mdy8LDwwEcpyO2bdtGTEwMlovuxNijRw/OnDnDwYMHOXnyJNnZ2cTExDjWh4SEEBkZ6Xi+YcMGDMOgRYsWTrVkZWVRq9blh6k/duwYycnJjB49mvj4eMfynJwcbDabo8527dpRvfqFn8fFNRXn4p+JxWKhTp06jve/fv16du/eTVBQkNM2Z8+edZyKSU1N5ZlnnmHJkiUcPXqU3NxcMjIySEpKctQWERFBgwYNSl3bxTZu3MiUKVNITEzkxIkT5OWZf/kkJSXRqlUr6taty6233srMmTPp0qULCxYs4OzZs/zpT39yvJczZ84U+JlnZmY63gtAo0aNqF27dqnrE5Eqoqj+MT0fhesegfUfmSNgO8mDr0bCDU9C095w+ojb7+ar8GKxlOrUjTtd2vnSYrE4LcsPKflfjIZhOAWX/GX5bY0S9LnJy8vDy8uL9evXFzgNERgYWKLtwTx11LVrV6d1+fsrSR1FKexnkv+aeXl5dOrUiTlz5hTYLv8LfuTIkRw7doxp06bRqFEj/Pz8iImJITs7u8jaLv2ZXk56ejqxsbHExsby6aefUrt2bZKSkujTp4/jdQDuu+8+4uLieP3115k1axZ33XWXI9Dl5eVRt27dQvsaXXyZekCAZ/xbFhE3KuqGeRYLXBNb+Kmlg2vg06EQUNsceNLNp5QUXiqxVq1aMXfuXKcQs3LlSoKCgqhfvz7BwcH4+PiQkJBAREQEACdPnmTnzp1cf/31AHTo0IHc3FxSU1Pp2bNnqWsIDw+nfv367N27l7vvvrvIOmfPnk1mZibVqlUDICEh4UrespOOHTvyxRdfODq4FmbFihW8/fbb9O/fH4Dk5GSOHz/uVFtSUhKHDx+mXj2zb9SqVatKVcf27ds5fvw4L7zwAg0bNgRg3bp1Bdr179/f0Zn2+++/Z/ny5U7vJSUlBW9vbxo3blyq1xcRKbHC7vbbe5IZWDZ+CunHLrQ18sx2bribr642qsQefPBBkpOT+etf/8r27dv5z3/+w7PPPssjjzyC1WolMDCQ0aNHM3HiRH766Se2bNnCyJEjsVov/LNo0aIFd999NyNGjOCbb75h3759rF27lhdffJGFCxeWqI4pU6YwdepU3njjDXbu3MnmzZuZNWsWr732GgDDhw/HarUyevRotm7dysKFCwtceXMl7r77bkJDQxk8eDArVqxg3759LFu2jPHjx3PwoHm5YPPmzZk9ezbbtm1j9erV3H333Y4ABXDzzTcTGRnJiBEj+O2331ixYgVPPfVUqeqIiIjA19eXN998k7179/Ldd9/x3HPPFWjn5eXFyJEjmTRpEs2bN3c6PXXzzTcTExPDbbfdxuLFi9m/fz8rV67kb3/7W6FB6HK2bt3qOIVlt9tJTEws9oaAIlKFdBwBEzbDPQvMx14Tod+LcMfMgm2NXPMUVDlTeKnE6tevz8KFC1mzZg3t2rVjzJgxjB49mr/97W+ONi+//DK9evVi0KBB3HzzzVx33XV06tTJaT+zZs1ixIgRPProo0RGRjJo0CBWr17tOIpwOffddx8ffPABH330EdHR0Vx//fV89NFHNGnSBDBPP82fP5+tW7fSoUMHnnrqKV588cWrfv/Vq1dn+fLlREREMHToUKKiohg1ahSZmZmOIzEzZ87k5MmTdOjQgbi4OMaNG0dYWJhjH1arlXnz5pGVlUWXLl247777+L//+79S1VG7dm0++ugjvvrqK1q1asULL7xQZDgbPXo02dnZjBo1ymm5xWJh4cKF9OrVi1GjRtGiRQuGDRvG/v37HX2dSqN///506NCB+fPns3TpUjp06ECHDh1KvR8RqaQKu/Q6vI15quhiFi+z70w5sxhX0+GgAkpLS8Nms2G32wucKjh79iz79u2jSZMm+PtXvrEexPP9+uuv3HDDDRw8ePCKQomr6P+OiACw4ZPCB5AsA8V9f19KfV5EKoCsrCySk5N5+umnufPOOytUcBERcSjqaqVyptNGIhXAv//9byIjI7Hb7bz00kvuLkdEpGglGUDSxRReRCqAkSNHkpuby/r166lf332/EEREPIHCi4iIiHiUKhleKlkfZRGX0/8ZEalIqlR4yb8ba0ZGhpsrEfEs+f9nLr2jsYiIO1Spq428vLyoWbOmY+yb6tWrl/pW7yJViWEYZGRkkJqaSs2aNTVStYhUCFUqvADUqVMHuDB4oYhcXs2aNR3/d0RE3M2l4eXkyZOMGzeO7777DoBBgwbx5ptvOg0kd6mijoS89NJLTJw48aprslgs1K1bl7CwMM6dO3fV+xOp7Hx8fHTERUQqFJeGl+HDh3Pw4EEWLVoEwF/+8hfi4uKYP39+kdscOXLE6fn333/P6NGjuf3228u0Ni8vL/1CFhER8UAuCy/btm1j0aJFJCQk0LVrVwBmzJhBTEwMO3bsIDIystDtLj00/Z///IfevXvTtGn5j50gIiIiFY/LrjZatWoVNpvNEVwAunXrhs1mY+XKlSXax9GjR/nvf//L6NGji2yTlZVFWlqa0yQiIiKVl8vCS0pKitPovPnCwsJISUkp0T4+/vhjgoKCGDp0aJFtpk6dis1mc0wlHelYREREPFOpTxtNmTKFv//978W2Wbt2LVB451vDMEp8efLMmTO5++67ix3FdtKkSTzyyCOO53a7nYiICB2BERER8SD539sluSlmqcPL2LFjGTZsWLFtGjduzKZNmzh69GiBdceOHSvRiLkrVqxgx44dfPHFF8W28/Pzw8/Pz/E8/83rCIyIiIjnOX36NDabrdg2pQ4voaGhhIaGXrZdTEwMdrudNWvW0KVLFwBWr16N3W6ne/ful93+ww8/pFOnTrRr165U9dWrV4/k5GSCgoIKPcJz7bXXOo4MFaW4Nley7tLlaWlpNGzYkOTkZGrUqFFsLeWhJD+T8tpfabbVZ1lQVf0si1uvz/Lq96fP8urosyzZZ2kYBp06daJevXqXrc1lVxtFRUXRt29f4uPjee+99wDzUukBAwY4XWnUsmVLpk6dypAhQxzL0tLS+Oqrr3j11VdL/bpWq5UGDRoUud7Ly+uy/5iLa3Ml64paXqNGjQrxH6skP5Py2l9pttVnWVBV/SyLW6/P8ur3p8/y6uizLPln5uvri9V6+e64Lh3baM6cOURHRxMbG0tsbCxt27Zl9uzZTm127NiB3W53Wvb5559jGAb/7//9vzKv6aGHHrqqNleyriSv6U5lXd/V7K802+qzLKiqfpbFrddnefX702d5dfRZXv3yS1kMDRdb7tLS0rDZbNjt9grxV4FcOX2WlYc+y8pDn2XlV6VGla4o/Pz8ePbZZ506Gotn0mdZeeizrDz0WVZ+OvIiIiIiHkVHXkRERMSjKLyIiIiIR1F4EREREY+i8CIiIiIeReFFREREPIrCiwfIyMigUaNGPPbYY+4uRa7C6dOnufbaa2nfvj3R0dHMmDHD3SXJFUpOTuaGG26gVatWtG3blq+++srdJclVGDJkCMHBwdxxxx3uLkVKSJdKe4CnnnqKXbt2ERERwSuvvOLucuQK5ebmkpWVRfXq1cnIyKBNmzasXbuWWrVqubs0KaUjR45w9OhR2rdvT2pqKh07dmTHjh0EBAS4uzS5Aj///DNnzpzh448/5uuvv3Z3OVICOvJSwe3atYvt27fTv39/d5ciV8nLy4vq1asDcPbsWXJzc0s09LtUPHXr1qV9+/YAhIWFERISwokTJ9xblFyx3r17ExQU5O4ypBQUXq7C8uXLGThwIPXq1cNisfDtt98WaPP222/TpEkT/P396dSpEytWrCjVazz22GNMnTq1jCqW4pTH53nq1CnatWtHgwYNePzxx0s0QruUXnl8lvnWrVtHXl4eDRs2vMqqpTDl+VmK51B4uQrp6em0a9eOt956q9D1X3zxBRMmTOCpp55i48aN9OzZk379+pGUlORo06lTJ9q0aVNgOnz4MP/5z39o0aIFLVq0KK+3VKW5+vMEqFmzJr/99hv79u3js88+4+jRo+Xy3qqa8vgsAf744w9GjBjB+++/7/L3VFWV12cpHsaQMgEY8+bNc1rWpUsXY8yYMU7LWrZsaTz55JMl2ueTTz5pNGjQwGjUqJFRq1Yto0aNGsbf//73sipZiuGKz/NSY8aMMb788ssrLVFKyFWf5dmzZ42ePXsan3zySVmUKSXgyv+XP//8s3H77bdfbYlSTnTkxUWys7NZv349sbGxTstjY2NZuXJlifYxdepUkpOT2b9/P6+88grx8fE888wzrihXLqMsPs+jR4+SlpYGmKPeLl++nMjIyDKvVYpXFp+lYRiMHDmSG2+8kbi4OFeUKSVQFp+leCZvdxdQWR0/fpzc3FzCw8OdloeHh5OSkuKmquRKlcXnefDgQUaPHo1hGBiGwdixY2nbtq0rypVilMVn+euvv/LFF1/Qtm1bRx+M2bNnEx0dXdblSjHK6vdsnz592LBhA+np6TRo0IB58+Zx7bXXlnW5UoYUXlzMYrE4PTcMo8Cykhg5cmQZVSRX42o+z06dOpGYmOiCquRKXM1ned1115GXl+eKsuQKXO3v2cWLF5d1SeJiOm3kIqGhoXh5eRVI/6mpqQX+SpCKT59n5aHPsvLQZ1l1Kby4iK+vL506deLHH390Wv7jjz/SvXt3N1UlV0qfZ+Whz7Ly0GdZdem00VU4c+YMu3fvdjzft28fiYmJhISEEBERwSOPPEJcXBydO3cmJiaG999/n6SkJMaMGePGqqUo+jwrD32WlYc+SymUG6908ng///yzARSY7rnnHkebf/3rX0ajRo0MX19fo2PHjsayZcvcV7AUS59n5aHPsvLQZymF0dhGIiIi4lHU50VEREQ8isKLiIiIeBSFFxEREfEoCi8iIiLiURReRERExKMovIiIiIhHUXgRERERj6LwIiIiIh5F4UVEREQ8isKLiIiIeBSFFxEREfEoCi8iIiLiUf4/hvzhw/0OT/0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogx(tobs, hobs0, \".C0\", label=\"obs layer 0\")\n",
    "plt.semilogx(tobs, hobs1, \".C1\", label=\"obs layer 1\")\n",
    "\n",
    "hm = ml.head(robs, 0, tobs)\n",
    "plt.semilogx(tobs, hm[0], \"C0\", label=\"modelled head layer 0\")\n",
    "plt.semilogx(tobs, hm[1], \"C1\", label=\"modelled head layer 1\")\n",
    "\n",
    "plt.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate data for head measured in well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "tobs2 = np.hstack((tobs, np.arange(0.61, 1, 0.01)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "self.neq  1\n",
      "solution complete\n"
     ]
    }
   ],
   "source": [
    "ml = ttim.ModelMaq(kaq=60, z=(-18, -25), Saq=1e-4, tmin=1e-5, tmax=1)\n",
    "w = ttim.Well(ml, xw=0, yw=0, rw=0.3, res=0.02, tsandQ=[(0, 788), (0.6, 0)], layers=0)\n",
    "ml.solve()\n",
    "rnd = np.random.default_rng(2)\n",
    "hobs2 = w.headinside(tobs2)[0] + 0.05 * rnd.random(len(tobs2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "....................................................................................................................................................................................................................................................................................................................................................................................................................\n",
      "Fit succeeded.\n"
     ]
    }
   ],
   "source": [
    "cal = ttim.Calibrate(ml)\n",
    "cal.set_parameter(name=\"kaq0\", layers=0, initial=100)\n",
    "cal.set_parameter(name=\"Saq0\", layers=0, initial=1e-3)\n",
    "cal.set_parameter_by_reference(name=\"res\", parameter=w.res[:], initial=0.05)\n",
    "cal.seriesinwell(name=\"obs1\", element=w, t=tobs2, h=hobs2)\n",
    "cal.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGhCAYAAACphlRxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCxklEQVR4nO3de1xUdf7H8dfMqHgDvJCAcfGGmqmVmrfNFSlveSkrV9fNclPEXLe1fqZZu6Vt5WrabrtmiVFZ2aZb2sXUdAvNTUu0qLQiKBFGxGuClqEy5/fHEQRBZJRhODPv5+NxHjpnzsx8jueB8+H7+Zzv12YYhoGIiIiIRdi9HYCIiIiIO5S8iIiIiKUoeRERERFLUfIiIiIilqLkRURERCxFyYuIiIhYipIXERERsZRa3g6gqrlcLnJycggMDMRms3k7HBEREakEwzA4duwYzZs3x26veGzF55KXnJwcIiMjvR2GiIiIXITs7GwiIiIqPMbnkpfAwEDAPPmgoCAvRyMiIiKVkZ+fT2RkZPH3eEV8LnkpKhUFBQUpeREREbGYyrR8qGFXRERELKVakpdFixbRsmVL6tatS9euXdm8eXOFx2/atImuXbtSt25dWrVqxXPPPVcdYYqIiIgFeDx5Wb58OVOnTuWhhx7i888/p0+fPgwePJisrKxyj9+9ezc33ngjffr04fPPP+fBBx/knnvu4c033/R0qCIiImIBNsMwDE9+QI8ePejSpQvPPvts8b4rrriCm2++mTlz5pQ5fsaMGbzzzjt88803xfsmTZrEF198wdatWy/4efn5+QQHB5OXl6eeFxEREYtw5/vboyMvJ0+eZMeOHQwYMKDU/gEDBrBly5ZyX7N169Yyxw8cOJDt27dz6tSpMscXFBSQn59fahMRERHf5dHk5dChQxQWFhIaGlpqf2hoKLm5ueW+Jjc3t9zjT58+zaFDh8ocP2fOHIKDg4s3zfEiIiLi26qlYffc254Mw6jwVqjyji9vP8DMmTPJy8sr3rKzs6sgYhEREampPDrPS0hICA6Ho8woy4EDB8qMrhQJCwsr9/hatWrRtGnTMscHBAQQEBBQdUGLiIhIjebRkZc6derQtWtXNmzYUGr/hg0b6N27d7mv6dWrV5nj169fT7du3ahdu7bHYhURERFr8HjZ6L777uP555/nhRde4JtvvuHee+8lKyuLSZMmAWbZ54477ig+ftKkSezZs4f77ruPb775hhdeeIGkpCSmTZvm6VBFRETEAjy+PMCoUaM4fPgwjz76KPv27aNjx46sWbOG6OhoAPbt21dqzpeWLVuyZs0a7r33Xp555hmaN2/OP//5T2699VZPhyoiIiIW4PF5Xqqb5nkRERHxDKfTSXp6OjExMRdc+dldNWaeFxEREfENSUlJREdHExcXR1RUFE8++aTXYvG5VaVFRESk6jidTrZs2cLEiRNxuVyAOYXJ9OnTsdlsXulJ1ciLiIiIlKtotGXUqFHFiUtJM2bMwOl0VntcSl5ERESkDKfTWWq0pTwul4uMjIxqjMqkspGIiIgAZ0tEAB9++GG5ictAoDnwImC322nTpk21xghKXkRERASzRBQfH8/5bkK+GXgYuAbIA94AHp47t8rvOqoMJS8iIiJ+LiUlhQkTJpT73HXAAqD7mcfHgSRgxn33eW0CWfW8iIiI+LGkpCS6d+9eZn8osAzYjJm4HAceA6KB7xISeGjBguoMsxRNUiciIuKnnE4nUVFRZUpFo4FngCZAIbAEeAQ4ANhsNrKysrw6SZ3KRiIiIn7q6aefLpW41AcWAXeeefwZEH/mzyLz5s3zSp9LSUpeRERE/JDT6WRBidJPFPAOcBXmaMtfgceB0yVeM2/evBqxULKSFxERET+Unp5ePOrSGVgHhAO5wCjgoxLHJiQk8Oc//9nrIy5FlLyIiIj4oYYNGwLQFdgANAa+BIYAQxISmHL99QD06tWrxiQtRZS8iIiI+BGn08nTTz/N/Pnz6QCsx0xcPsZMXI7Z7TVqlKU8Sl5ERET8RMmJ6JoBazDvKPoEGASccDhIXLy4RicuoORFRETELxStVWQYBrUwZ8iNBtIwR1yOA88uXMj48eO9GWalaJI6ERERP/D0008Xr1U0C+iDOc3/MODImWOaNm3qldjcpeRFRETEx5W8Lbo78MCZ/ROA9DN/t9ls9OrVywvRuU/Ji4iIiI8rui3aASQCDuBVzNJRkSVLltT4XpciSl5ERER83KJFiwBzpOUqzDLR1DPPDRs2jOzsbEv0uhRR8iIiIuLDUlJSeOONN6iHuT4RmD0vhwG73c6iRYssM+JSRMmLiIiID9u8eTMA4zFn0M0EnsPscUlMTLRc4gK6VVpERMSn9enTBztw35nHc4FTwLvvvMPQoUO9F9gl0MiLiIiID9u0aRNDgJaYpaKXgDvvvNOyiQto5EVERMRnPfnkk0yfPp23zjx+ASiw2Xjssce8GNWl08iLiIiID3I6ncyYMYOmwI1n9r0EGIZBRkaG9wKrAkpeREREfFDR3C7DgdpAKvA15h1Gbdq08Wpsl0rJi4iIiA+KiYnBbrcz/MzjVWf+nDt3riXvMCpJyYuIiIgPev/996nlctH/zON3gXnz5jFt2jRvhlUl1LArIiLiY4pWkO4FNAAOAl/YbLzz2996ObKqoZEXERERH5Oeno7L5SL2zOONgMsHGnWLKHkRERHxMUX9LtedebwRcDgclm/ULaLkRURExMdEREQw9vbb6XLmcQpw++23W75Rt4iSFxERER/jdDr58JVXaAacBr4CXn31VZxOp5cjqxpKXkRERHxMeno6VxsGYM7t8gtQWFionhcRERGpmWJiYuhkswHwxZl96nkRERGRGisiIoI7fvUrANIxE5fFixer50VERERqrnZ28yt+9IMPkpmZyfjx470cUdXxaPLy448/MnbsWIKDgwkODmbs2LEcPXq0wteMGzcOm81WauvZs6cnwxQREfE5p7/7DoCQnj19ZsSliEeTlzFjxpCamsq6detYt24dqampjB079oKvGzRoEPv27Sve1qxZ48kwRUREfMpLzz5LrdxcADredBNJSUlejqhq2QzjTDtyFfvmm2/o0KEDn3zyCT169ADgk08+oVevXnz77be0a9eu3NeNGzeOo0eP8tZbb1XqcwoKCigoKCh+nJ+fT2RkJHl5eQQFBV3yeYiIiFiJ0+lkUFQUOw2DPKARZs9LZmZmjR6Byc/PJzg4uFLf3x4bedm6dSvBwcHFiQtAz549CQ4OZsuWLRW+duPGjTRr1oy2bdsSHx/PgQMHznvsnDlzistSwcHBREZGVtk5iIiIWE16ejohZ8Yl9p3Z50u3SYMHk5fc3FyaNWtWZn+zZs3IPTOUVZ7BgwezbNkyPvzwQxYsWEBKSgpxcXGlRldKmjlzJnl5ecVbdnZ2lZ2DiIiI1cTExBB65jbpg2f2+dJt0nARycusWbPKNNSeu23fvh0A25l/vJIMwyh3f5FRo0YxZMgQOnbsyLBhw1i7di3fffcd7733XrnHBwQEEBQUVGoTERHxVxEREfxpzBjATF587TZpgFruvmDKlCmMHj26wmNatGjBl19+yf79+8s8d/DgQUJDQyv9eeHh4URHR5Oenu5uqCIiIn6pd9u2APQcMoTM557zqcQFLiJ5CQkJISQk5ILH9erVi7y8PLZt20b37t0B+PTTT8nLy6N3796V/rzDhw+TnZ1NeHi4u6GKiIj4p4Nmwaj5VVeBjyUu4MGelyuuuIJBgwYRHx/PJ598wieffEJ8fDxDhw4tdadR+/btWbVqFQDHjx9n2rRpbN26lczMTDZu3MiwYcMICQlhxIgRngpVRETEp/y8Zw8AR2vX9nIknuHReV6WLVtGp06dGDBgAAMGDKBz58688sorpY5JS0sjLy8PMOtyX331FTfddBNt27blzjvvpG3btmzdupXAwEBPhioiIuITkpKS2PLuuwD88dFHfW6OF/DgPC/e4s594iIiIr7E6XQSHR3NZy4XVwEDgQ8sMMcL1JB5XkRERKR6paen43K5uOzM44P43hwvoORFRETEZ8TExGC32ylqtMjH9+Z4ASUvIiIiPiMiIoLExMSzX+52u8/N8QLqeREREfE5Rt262AoK2PfJJ4SXWKanJlPPi4iIiB+zuVwAhF9+uZcj8QwlLyIiIr7mTPKCw+HdODxEyYuIiIivKSw0/7T75te8b56ViIiIvyrZyqrkRURERGq8opIRqGwkIiIiFlBUMgKNvIiIiIgFaORFRERELEUjLyIiImIpJUdelLyIiIhIjaeykYiIiFiKykYiIiJiKSobiYiIiKUoeRERERFL8fGlAUDJi4iIiG/x8UUZQcmLiIiIb9HIi4iIiFhK0ciLkhcRERGxBJWNRERExFJUNhIRERFLUdlIRERELEVlIxEREbEUlY1ERETEUjTyIiIiIpaikRcRERGxFDXsioiIiKWobCQiIiKWorKRiIiIWIpGXkRERMRSNPIiIiIilqKGXREREbEUlY1ERETEUlQ2EhEREUtR2UhEREQsRWWjS/P444/Tu3dv6tevT6NGjSr1GsMwmDVrFs2bN6devXrExsaya9cuT4YpIiLiO1Q2ujQnT55k5MiR3H333ZV+zbx583jqqadYuHAhKSkphIWF0b9/f44dO+bBSEVERHyEH4y81PLkm8+ePRuAl156qVLHG4bBP/7xDx566CFuueUWAJYuXUpoaCivvfYaCQkJZV5TUFBAQUFB8eP8/PxLD1xERMSqNPJSvXbv3k1ubi4DBgwo3hcQEEDfvn3ZsmVLua+ZM2cOwcHBxVtkZGR1hSsiIlLzqGG3euXm5gIQGhpaan9oaGjxc+eaOXMmeXl5xVt2drbH4xQREamx/KBs5HbyMmvWLGw2W4Xb9u3bLykom81W6rFhGGX2FQkICCAoKKjUJiIi4rf8oGzkds/LlClTGD16dIXHtGjR4qKCCQsLA8wRmPDw8OL9Bw4cKDMaIyIiIuXwg5EXt5OXkJAQQkJCPBELLVu2JCwsjA0bNnDNNdcA5h1LmzZtYu7cuR75TBEREZ/iByMvHj2zrKwsUlNTycrKorCwkNTUVFJTUzl+/HjxMe3bt2fVqlWAWS6aOnUqTzzxBKtWrWLnzp2MGzeO+vXrM2bMGE+GKiIi4hv8oGHXo7dKP/zwwyxdurT4cdFoSnJyMrGxsQCkpaWRl5dXfMz06dM5ceIEkydP5scff6RHjx6sX7+ewMBAT4YqIiLiG/ygbGQzDMPwdhBVKT8/n+DgYPLy8tS8KyIi/icpCSZMgKFD4d13vR1Npbnz/e27Y0oiIiL+yA/KRr57ZiIiIv7ID8pGSl5ERER8ie42EhEREUvRyIuIiIhYikZeRERExFLUsCsiIiKWorKRiIiIWIrKRiIiImIpKhuJiIiIpahsJCIiIpaispGIiIhYikZeRERExFI08iIiIiKWooZdERERsRSVjURERMRSVDYSERERS9HIi4iIiFiKRl5ERETEUtSwKyIiIpaispGIiIhYispGIiIiYikqG4mIiIilqGwkIiIilqKykYiIiFiKRl5ERETEUjTyIiIiIpaihl0RERGxFJWNRERExFJUNhIRERFLUdlIRERELKVo5EVlIxEREbEEjbyIiIiIpahhV0RERCxFDbsiIiJiKSobiYiIiKWobCQiIiKWorLRpXn88cfp3bs39evXp1GjRpV6zbhx47DZbKW2nj17ejJMERER36GRl0tz8uRJRo4cyd133+3W6wYNGsS+ffuKtzVr1ngoQhERER/jByMvtTz55rNnzwbgpZdecut1AQEBhIWFeSAiERERH6eGXe/YuHEjzZo1o23btsTHx3PgwIHzHltQUEB+fn6pTURExG+pbFT9Bg8ezLJly/jwww9ZsGABKSkpxMXFUVBQUO7xc+bMITg4uHiLjIys5ohFRERqED8oG7l9ZrNmzSrTUHvutn379osOaNSoUQwZMoSOHTsybNgw1q5dy3fffcd7771X7vEzZ84kLy+veMvOzr7ozxYREbE8Pygbud3zMmXKFEaPHl3hMS1atLjYeMoIDw8nOjqa9PT0cp8PCAggICCgyj5PRETE0vygbOR28hISEkJISIgnYinX4cOHyc7OJjw8vNo+U0RExLJUNro0WVlZpKamkpWVRWFhIampqaSmpnL8+PHiY9q3b8+qVasAOH78ONOmTWPr1q1kZmayceNGhg0bRkhICCNGjPBkqCIiIr5BIy+X5uGHH2bp0qXFj6+55hoAkpOTiY2NBSAtLY28vDwAHA4HX331FS+//DJHjx4lPDycfv36sXz5cgIDAz0ZqoiIiG/wg5EXm2EYhreDqEr5+fkEBweTl5dHUFCQt8MRERGpXldeCV9/DR9+CP36eTuaSnPn+9t30zIRERF/5AdlIyUvIiIivsQPyka+e2YiIiL+SCMvIiIiYikaeRERERFL8YMZdn33zERERPyRykYiIiJiKSobiYiIiKWobCQiIiKWorKRiIiIWIrKRiIiImIpGnkRERERS9HIi4iIiFiKGnZFRETEUlQ2EhEREUtR2UhEREQsRWUjERERsZSikReVjURERKTGMwxzA428iIiIiAUUJS6gkRcRERGxgKKSEWjkRURERCygqFkXlLyIiIiIBZRMXlQ2EhERkRpPZSMRERGxFI28iIiIiKVo5EVEREQsRQ27IiIiYilKXkRERMRSispGNpu5+SglLyIiIr7CDxZlBCUvIiIivqMoefHhO41AyYuIiIjvKCobaeRFRERELEEjLyIiImIpGnkRERERS1HDroiIiFiKykYiIiJiKSobiYiIiKVo5EVEREQsRSMvlyYzM5Px48fTsmVL6tWrR+vWrXnkkUc4efJkha8zDINZs2bRvHlz6tWrR2xsLLt27fJUmCIiIr5DDbuX5ttvv8XlcrF48WJ27drF3//+d5577jkefPDBCl83b948nnrqKRYuXEhKSgphYWH079+fY8eOeSpUERER3+AnZSObYRhGdX3Yk08+ybPPPssPP/xQ7vOGYdC8eXOmTp3KjBkzACgoKCA0NJS5c+eSkJBwwc/Iz88nODiYvLw8goKCqjR+ERGRGu3TT6FnT2jRAnbv9nY0bnHn+7tax5Xy8vJo0qTJeZ/fvXs3ubm5DBgwoHhfQEAAffv2ZcuWLeW+pqCggPz8/FKbiIiIX1LZqGp9//33/Otf/2LSpEnnPSY3NxeA0NDQUvtDQ0OLnzvXnDlzCA4OLt4iIyOrLmgREREr8ZOykdvJy6xZs7DZbBVu27dvL/WanJwcBg0axMiRI5kwYcIFP8Nms5V6bBhGmX1FZs6cSV5eXvGWnZ3t7imJiIj4Bj+526iWuy+YMmUKo0ePrvCYFi1aFP89JyeHfv360atXLxITEyt8XVhYGGCOwISHhxfvP3DgQJnRmCIBAQEEBARUMnoREREf5icjL24nLyEhIYSEhFTq2L1799KvXz+6du3Kiy++iP0CmWDLli0JCwtjw4YNXHPNNQCcPHmSTZs2MXfuXHdDFRER8S9+MvLisbPLyckhNjaWyMhI5s+fz8GDB8nNzS3Tu9K+fXtWrVoFmOWiqVOn8sQTT7Bq1Sp27tzJuHHjqF+/PmPGjPFUqCIiIr7BTxp23R55qaz169eTkZFBRkYGERERpZ4reXd2WloaeXl5xY+nT5/OiRMnmDx5Mj/++CM9evRg/fr1BAYGeipUERER3+AnZaNqneelOmieFxER8Vtr1sCQIdC1K5xz80xNV2PneREREREP8pOykW+fnYiIiD8patj18bKRkhcRERFfoZEXERERsRQ/adhV8iIiIuIrNM+LiIiIWIrKRiIiImIpKhuJiIiIpahsJCIiIpaikRcRERGxFI28iIiIiKWoYVfO5XQ6SU5Oxul0ejsUERGRslQ2kpKSkpKIjo4mLi6O6OhokpKSvB2SiIhIaSobSRGn08nEiRNxncloXS4XCQkJGoEREZGaRWUjKZKeno7L5aJeiX2FhYVkZGR4LSYREZEy/KRsVMvbAVhBTEwMdpuNXYbBIWA1sNZup02rVt4OTURE5CyVjaRIREQE/37sMVoC1wKzgW0uFxE9e0J8PLz9Nvz0k5ejFBERv+cnIy9KXirpNw8+SM6OHXx7//2cGDQIGjaEffvg+efh5puhaVMYPBieeQYyM70droiI+CM/GXlR2cgNzbt0oXmXLuaDggL46CNYvRrefRd274Z168xtyhS48koYOtTcevaEWvqnFhERD1PDrlQoIAD694enn4bvv4evv4Z58+DXvzaH63btgrlzoU8fCA2F22+H11+HH3/0duQiIuKrVDaSSrPZ4Ior4P77YdMmOHgQ/v1v+N3voHFjOHIEli2D3/4WLrsMYmNh/nz45hswDG9HLyIivsJPyka+fXbe0rgxjB4Nr74KBw7A5s0wY4ZZSiosNBOc+++HDh2gTRv4059gwwazFCUiInKxVDaSKlGrFlx3Hfztb7BzJ/zwA/zrXzBwINSpYz7+5z9hwAAICYFbb4UXX4T9+70duYiIWE3RyIvKRlKlWrY0G3rXrYPDh+Gtt2DCBAgLg+PHYeVKuOsu83GPHvDXv8Lnn6u8JCIiF6aRF/G4hg3hpptgyRLYuxe2b4dZs6BbN/P5bdvg4YehSxeIiICJE+Gddyo1p4wWkRQR8UNq2JVqZbdD167wyCOQkgI5OWfnkGnQwHy8ZImZ7DRtCjfeCIsWwZ49Zd5Ki0iKiPgpP2nYtRmGb9Uj8vPzCQ4OJi8vj6CgIG+HUzUKCswm36I5Zc6dBK9TJ3M+mSFDcEZEEN2qVfEikgAOh4PMzEwiIiKqN24REale06fDk0/C//2feVerhbjz/e3bqZmvCAgwG3r/+U+zwbfkHDJ2O3z1FcyZA9ddR7POnXnJ5eI3QPCZl2sRSRERP+EnZSNN+2o1Npt5i3WHDmaGfeSI2fy7ejWsW0edH39kLDAWOA1sBtbYbLS32cymX5vNu/GLiIjn+EnZyLfPzh80aQJjxsBrr5lzynz0EV8OGsTXmJlpP+BJwyAsNhbatoWpU+G//4WTJ70atoiIeICfjLwoefEltWpBnz50XruWoOxsPlm2jB8ffdQsOdWpAxkZ5nIG/fubc8rcdhu89JKZ9IiIiPX5yciLykY+KiIigogxY87uOH7cHHFZvRreew9yc+HNN83NZoPu3Yubfrn6apWXRESsSPO8iE9p2NC87fr55805ZVJSzNuyu3Y1e2E+/RT+8hdzTpnISEhIMO9s+vnni/5IzTUjIlLNVDYSn2W3mxPhzZplToy3d+/ZOWTq1zcfJybC8OHmnDJDhsCzz0JWVqU/QnPNiIh4gZ+UjXz77KRymjc3lyh46y1zyYJ168wlDKKj4ZdfYM0amDzZfHzVVfDQQ7B169kfknM4nU4mTpxYPNeMy+UiISFBIzAiIp6mspH4pbp1zUUj//Uv2L3bXEzyb38zF5e02+HLL+GJJ6B3b3P9pTvugBUrIC+v+C3S09NLTZIHmmtGRKRaaGFG8Xs2G1x5JcyYAZs3m3clvfoqjB4NjRrBoUPwyiswapR591JcHDz1FFc4HNjPyfodDgdt2rTxznmIiPgLPxl50d1GUnlNm8Lvfmdup0/Dli3m3UurV8M330ByMiQnEwYcadaMlw4e5B3DYIvdzsLFi7U8gYiIp6lh99JkZmYyfvx4WrZsSb169WjdujWPPPIIJy8wOdq4ceOw2Wyltp49e3oqTLlYtWrBr38N8+bB11+XnkOmdm2CDxzgT4bBB8BP9eszft06WLpUc8qIiHiSnzTsemzk5dtvv8XlcrF48WLatGnDzp07iY+P56effmL+BRaLGjRoEC+++GLx4zp16ngqTKkqrVvDPfeY27FjsGGDOSKzZg32/fvhjTfMzWaDHj3MOWWGDoXOnTWnjIhIVVHZ6NIMGjSIQYMGFT9u1aoVaWlpPPvssxdMXgICAggLC/NUaOJpgYFwyy3m5nKZt2O/956ZzHz2GXzyibn9+c8QEWHeij10qNkzU79+8ds4nU7S09OJiYlRyUlEpDJUNqp6eXl5NGnS5ILHbdy4kWbNmtG2bVvi4+M5UEGpoaCggPz8/FKb1CB2uzl77+zZsGMHOJ1n55CpV898vHgxDBtm9tQMHQrPPcfrTz6peWJERNzlJ2Ujm2EYRnV80Pfff0+XLl1YsGABEyZMOO9xy5cvp2HDhkRHR7N7927+8pe/cPr0aXbs2EFAQECZ42fNmsXs2bPL7M/LyyMoKKhKz0Gq2IkTsHGjOSLz7ruQnV3q6VRg9ZntM7udH/bs0QiMiEhFbr4Z3n7b/CUxPt7b0bglPz+f4ODgSn1/u528nC9ZKCklJYVu3boVP87JyaFv37707duX559/3p2PY9++fURHR/P6669zyy23lHm+oKCAgoKC4sf5+flERkYqebEawzDnlFm9mrxlywjctavUsOABwDVwIGETJpgLTVby2qr0JCJ+Zdgw8xfC55+H8eO9HY1b3Ele3O55mTJlCqNHj67wmBYtWhT/PScnh379+tGrVy8SExPd/TjCw8OJjo4mPT293OcDAgLKHZERi7HZoFMn6NSJY2PH0jYqioGGwRBgENAM4P33za1WLejb92zT73nmj0lKSiqe6ddut5OYmMh4i/0wi4i4RQ275QsJCSEkJKRSx+7du5d+/frRtWtXXnzxxTITl1XG4cOHyc7OJjw83O3XijVFRETwxJIlJCQk8EphIQF2Oyvvu48bXS6z8TctDT74wNzuvRfatTvb9HvddVC79nmXKBg4cKBGYETEd6lh99Lk5OQQGxtLZGQk8+fP5+DBg+Tm5pKbm1vquPbt27Nq1SoAjh8/zrRp09i6dSuZmZls3LiRYcOGERISwogRIzwVqtRA48ePJzMzk+TkZDL27OHGJ5+EBQvg22/hu+/g73+H6683R2HS0uCpp8y7lS67DEaNIv+ZZ2isJQpExN/4ScOux26VXr9+PRkZGWRkZJT5Tbdkm01aWhp5Z9bFcTgcfPXVV7z88sscPXqU8PBw+vXrx/LlywkMDPRUqFJDRURElD9KEhMDU6eaW15eqTllOHgQVqygA2afzFbONv1+Y7e7tUSB+mVExHL8pGxUbXcbVRd3Gn7ExxQWQkqKmci89x6kppZ6+niTJjQcPdosL8XGmrdqn4f6ZUTEkuLizKVaXn/dXHfOQjx6t1FNp+RFimVn8+Orr1L47rs0/ewzbCXuSqNePbjhBjORGTIELr+8+Cmn00l0dHSplbEdDgeZmZkagRGRmq1vX/joI1ixAkaO9HY0bnHn+9u3x5XEv0VG0njmTEK2bMF25Ig5IjNpkjmr74kT5twyCQnm4y5d4OGHYds20tPSSiUuoH4ZEbEINeyK+JD69c0Rlmefhawss6T02GPQq5d5m/bnn8Nf/wo9evDr0aN5ARgBNDzzcofD4Va/jIiIV6hhV8RH2Wxw1VXm9tBDZpPv2rXmyMz77+M4dIjfA78HTgIfAQ1HjiSiZNlJRKQm8pOGXd8+O5HKuOwyuOMOs0Z88GDx/DGnWrakDnAD0PP1183J8K64Au6/HzZtglOnvB25iEhpKhuJ+KE6dcxu/aeeovYPP5hzyCxYAP36mXPKfPstzJ9v3q3UrBn89rewbBkcPuztyEVE/KZspLuNRCrr6FFYv/7snDIlExa7HXr3Lr57yRkcTHpGhuaIEZHqdc01Zk/funUwcKC3o3GL7jYS8YRGjeA3v4GXX4b9++Hjj2HmTHNNJpcL/vc/eOAB6NSJ01FR7IqLY2JUFC8995y3IxcRf1E08qKykYiU4XCYIy1PPAFffgmZmfDMM5zo149fgBbAFGCNYTDy7rs5MXCgucprTo5XwxYRH+cnDbu620ikKkRHw+TJfHLFFQxNTuZ6YAgwFLgczHLT+vXmsV27wtCh7L/2Wr6uW5eYdu1UWhKRquEnDbtKXkSqUExMDL/Y7bzrcvHumX1d7XY+uPdegjdvhm3bYMcO2LGDUMAFrAWiJk/mhrlzoWHD87+5iMiF+EnDrm+fnUg1i4iIIDExEceZ33ocDgd3JyYSPH8+fPop7NvHkfnzWQkcA8KBu4AbFi3CaNrUbLD7179g927AXKogOTkZp9PprVMSESvxk7KR7jYS8QCn00lGRgZt2rQpUxJKTk4mLi6OOkBfzpaXWp/zHj82b86SnBzeBT612Xh2yZJyF4fU6tciUqx1a/jhB9i6FXr29HY0btHdRiJeFhERQWxsbLnJRExMDHa7nZPABmAq0M5uJzc5GZ58Evr2xXA4aJyTw3RgM5BrGNSLj+fwwoXw44/F75WUlER0dDRxcXFER0eTlJRUPScoIjWTykYi4gnllZYWJyYSFhsL06bBxo38b+VKRgPLgCNAE2CMYdD0j380ZwTu25ejf/4zT8XHFy8i6XK5SEhIUIlJxJ/5ScOuykYiXlJRacnpdBIdHY3L5cIB9ASG22zcGxND7e++K3Xs98DqM9tHwPvJycTGxpZ6L5WVRPzE5ZebUzJ89pk5YZ2FqGwkYgEVlZZKjs4UAp84HDRdsoTaaWlmM+/ChfwSG0sBZq/MnzBLUIeA7vPmwQsvQG6uykoi/kYNu9akkRfxJRWNzgAsfeYZ3r7nHga7XAzFvHuppG2cHZX5HLNElZmZqREYEV8VGgoHDsBXX0HHjt6Oxi3ufH8reRGxuOIEp1UrIg4eNNdeWr0atm8vddxe4D2g12OP0WnqVGjQwBvhiognhYSY667t2gUdOng7GrcoeVHyIkLOjh08cu213GgY9AdKTX8XEGCunn1mIUmio70UpYhUqSZNzDsSv/kG2rf3djRuUc+LiNC8a1d6LlnCSIeDpsBgu51dcXHQogUUFMDatfCHP5iPO3eGBx+ELVvO3mopItbjJwszauRFxMeV6ZsxDPO3sqLy0scfn23yA2jaFAYP5nCvXnwdGUnLa65Rj4yIVQQGwvHjkJFhTlhnISobKXkRqbwjR2DdOjORWbsWjh4tfuo05iR59UeOpMdf/wpt24LN5q1IReRCGjSAn38270ps0cLb0bhFyYuSF5GLc/o0B956i5dGjmQoUKbdr00bs09m6FDo0wfq1AE0l4xIjVG3rlkW3rMHoqK8HY1b1PMiIhenVi12NW3KDOBKoBXwR+B9wFWrljkU/Y9/wA03mHc13HYbH911F92iojSXjEhNoHlerEkjLyKXpuTsvkUcDgd7du7k8qJemffeg/37i593cXZOmbV2O29nZhIRGVnmfTU6I+JhtWqZTbs5ORB+7sxPNZtGXkTkopW79tLixVzevj2MGAFJSeZ/jNu2sfuOO9iB+R9JT+AxYIfLxWVdusCkSWai8/PPmulXpLr4ycKMGnkRkXJdaHbfomOio6MJc7m4ERgK9AfqlzjGCAhgTUEB72JOkucE7HY7//73v+ndu7dGYUSqimGcTVoOHDAXcbUQjbyIyCWraO2lksckJiay3+HgeeBWh4MVixbBmjUweTJERWErKGAI8ByQjblMwaMuF38fNYqWUVGlRmGcTifJyclaGVvkYpSco0nzvFiLRl5Eqt95R2kMg9wPPuDpAQMYYhj0Akr+l3oQWGuzMWTRIt47dYrfT52Ky+XCbreTmJjI+PHjq/lMRCzs5Elz9mwwZ9lt1Mir4bhLt0oreRGpUZKSkkhISKBRYSGDgGHAQKBRiWNOAR9xdiHJ3ZVcRFKNwCJn/PIL1Ktn/j0/35ywzkJUNhKRGmX8+PFkZmby7IoV/NtuZzRwGRALPGWzkd+8ObWB64G/A+nArsJCjPvug+RkOHWq3PdVI7BICSXLRmrYtRaNvIjUbEWjMIWFhcV3Mg0cOJDro6IYbBgMBfoCtUu+KCgIBg0yJ8cbPBhCQs57S3dlRmtEfNKxY+bPCpiz7BaNwliEykZKXkRqtPJ6ZEomNY3sdt5ISOD6EyfMOWUOHjz7YpsNevbkhw4duDkpia/Oee/k5GRiY2Or7VxEaoy8vLN9LgUFxTNgW4WSFyUvIpZUbuOvywUpKWcXkkxNLfWaLM72yXxkt/Ptnj0aeRH/dOSIubAqmKXWWrW8G4+blLwoeRHxXdnZ5q3Yq1dz+v33qVWiH+Z0nTrUGjjQLC8NGQKXX+7FQEWq2aFDZ+d2cbkst4iqkhclLyL+4eefObRiBQWrVhGakkKtfftKP3/NNWcXkuzWrUwTo+5UEp+yfz+EhZl/t+BXu5IXJS8i/scw4Msvz5aXPv209H/gzZqZozFDh0L//iStWMHEiRM1r4z4jn37oHlzM0kveeeRRdSYW6WHDx9OVFQUdevWJTw8nLFjx5KTk1PhawzDYNasWTRv3px69eoRGxvLrl27PBmmiPgCmw2uugoeegi2bjV/C126FEaONOe7OHAAXnwRbr0Vo2lToidMYIrLRSvA5XKRkJCgmX3F2oruvPPx2XXBw8lLv379WLFiBWlpabz55pt8//333HbbbRW+Zt68eTz11FMsXLiQlJQUwsLC6N+/P8eOHfNkqCLiay67DO64A1asMHsBPvgA7r0X2rTBduoUNwBPA98DXwNzCgs5+MYb551TRqTG85NFGaGay0bvvPMON998MwUFBdSuXbvM84Zh0Lx5c6ZOncqMGTMAKCgoIDQ0lLlz55KQkHDBz1DZSEQuJPejj3gyNpYhhkEfzplTplGjs3PKDBp09u4NkZouMxNatjTnd/n5Z29H47YaUzYq6ciRIyxbtozevXuXm7gA7N69m9zcXAYMGFC8LyAggL59+7Jly5ZyX1NQUEB+fn6pTUSkImG//jUdlixhgMPBZcBom430Xr0gJASOHoXXX4fbbzf7ZPr0gblzYedOSzZBih8pGnlR2ejSzZgxgwYNGtC0aVOysrJ4++23z3tsbm4uAKGhoaX2h4aGFj93rjlz5hAcHFy8RUZGVl3wIuKzipYseCs5mflZWcRs2QK5ubBlCzz4IHTubPYQ/O9/8MAD0KmT+VvtlCmwbp25joxITVLU8+IHZSO3z3DWrFnYbLYKt+3btxcff//99/P555+zfv16HA4Hd9xxBxeqVNnOuTfdMIwy+4rMnDmTvLy84i07O9vdUxIRPxUREUFsbOzZ26QdDujVCx5/HL74AvbsgUWL4MYbzdV69+yBZ54xlyho2hRuvhmefx5ycnA6nSQnJ6vpV7zHjxp23e55OXToEIcOHarwmBYtWlC3bt0y+51OJ5GRkWzZsoVevXqVef6HH36gdevWfPbZZ1xzzTXF+2+66SYaNWrE0qVLLxifel5ExCN++gk+/PDsrdjn3Dm5HXOW3zU2GwmLFzM+Pt4rYYof+/pruPJKM7G+wPd0TeTO97fbcweHhIQQEhJyUYEV5UkFBQXlPt+yZUvCwsLYsGFDcfJy8uRJNm3axNy5cy/qM0VEqkSDBjBsmLkZhrlMwerVnFy5kjqpqXQDugGzDIN9Eyfy04cf0mDUKLjhBmjY0MvBi19Q2ejSbdu2jYULF5KamsqePXtITk5mzJgxtG7dutSoS/v27Vm1ahVgloumTp3KE088wapVq9i5cyfjxo2jfv36jBkzxlOhioi4x2YzZ+/9y1/4+KmnCAPuAt4EjgHhQIPXX4cRI8zfggcNgoULYfdur4YtPs6PykYeS17q1avHypUruf7662nXrh133XUXHTt2ZNOmTQQEBBQfl5aWRl5eXvHj6dOnM3XqVCZPnky3bt3Yu3cv69evJzAw0FOhiohctJiYGA7a7bwI3AaEAIPsdo7ddRe0agUnT8L778Mf/witWnG8ZUvy//AHsxH49GkvRy8+RfO8WJd6XkSkuiUlJZGQkEBhYSEOh4PFixebSw0YBqSlwerV7FuyhMu++650rb5JE7P5d+hQGDgQGjf21imIL9ixw1zDKzISsrK8HY3btLaRkhcRqWZOp5OMjAzatGlTZpFHp9NJdHQ0QS4Xg4AhwGCg1PR3Dgdcd13x+kvOhg1Jz8jQopFSedu2QY8eEB1tTlhnMR5t2BURkbIiIiLOm2Skp6fjcrk4Crx+ZnMAKf/8J9fs3WvevbRrF2zaZG7Tp1MAfAn8zWZj9KJF/H7SpOo6FbEqNeyKiEhViYmJwX7uF4rDwWUjRsDf/mbO3rt7NyxcyC+xsRQArYE/Ae8bBrfdfTc/33ijubDk/v2aU0bKp4ZdERGpKhERESQmJuI486VS1BdTaqSmRQv4wx/Y+vDDNAVuBpYA+4BAoP7atXDXXRAWRk5kJMlxcdwcFUXS88+f93OV5PgZP2rY9f0zFBGpAYqWI0hOTiYzM9Ns6C1HTEwMJ+x23gYmApcD3e128u67j5OdOwPQHXgU2G4YDIqP5/jvfgfvvGNOpHdGUlIS0dHRxMXFER0dTVJSkqdPUbxNZSMREalqZZYjOM8xJUdp7A4HCYmJBC9YwMf/+AfhwHhgJXAcM7lp+NprcNNN5pwyN97Ij48/zmPx8bjOfJm5XC4SEhI0AuPrtDCjiIh4y/lGaWJiYjhgt/MCcCvm3UqD7XaO/f73ZtmpoADWrqXxn//MbsPgC+BxoBdgFBaSkZHhrVOS6qCRFxER8abyRmnOHZUpdDi4LTGRwBdegB9+MO9YmjuXgu7dKQQ6Aw8CW4D9wLULF8KKFXD0aPWfkHieHzXsap4XERGLqWhOmSKvPv007997L4MNg8FAqenvatWCPn3MyfGGDoW2bcu8f3p6uuaYsZq1a80V0Lt0MSessxhNUqfkRUTkbJLTogURWVlnV8T+5pvSB8bEFCcyL6anM2HyZFwuF3a7ncTExPM2F0sN89575nXs1g1SUrwdjduUvCh5ERE5v++/N7/oVq+GjRvh1Knip/KB94HVwFrgiMNBZmamRmCs4N13Yfhwc5bdTz7xdjRuc+f7Wz0vIiL+pnVruOceWL8eDh+GN9+E3/+egsaNCQJGAkuBXGBzYSEnZ82CL74w12qSmkvzvIiIiF8IDIRbboEXXuBgaio9bDYeBXZgfkH0AlolJcHVV0NUFNx9tzlqc+KEV8OWcuhuIxER8TcRUVFMXLKERx0OugFRdjv/u+MOsxRRrx44nfDcczB0KIWNG3Oif3/zcXa2t0MX8Kt5XrQwo4iIFBs/fjwDBw4sezfTiROwcSNfP/kkDZKTiS4ooN5//wv//a/5/FVXnb176dpr/eILtMbxo5EXJS8iIlJKuStk16uHs1MnOm3ahAu4EhgKDAN62+3YvvjC7It5/HG47DLzlt2hQ2HAADhP86Vuya5ifjTPi++nZyIiUiXS09OLlxzYBcwFrgO2rFwJr7wCo0ZBcDAcPAhLl8LIkeaSBTfcAH//O6SnF7+X1l7yADXsioiIlBYTE4P9nC9Gh8NBdNeucPvt8PrrZuKSnAz/93/Qrh2cPg0ffAD33WdOhteuHccmTuS1+HjsWnupavlR2cj3z1BERKrEucsTOBwOFi9eXLrkU7s2xMbC/Pnw7bfw3XfmqMsNN5jPffcdgUuW8IFhcAh4HbgdaKS1ly6dH5WNNEmdiIi4pTLLE5QrPx82bOCn5cv56T//oVmJp1zAqW7dCBgxwuyV6dQJbLaqDt23JSXBhAnmv9+773o7GrdpkjoREfGY8haNrJSgILj1VhqsWMG7iYn8ym7nr8DnmF9GAdu3w0MPmXcuRUfD5MmwZo3mlKksPxp5UfIiIiLVbnx8PMv37KFPcjKXZWdDVlbxHDLUrWvOHfPsszBkiNn0O3w4JCbC3r3eDr3m8qOGXd0qLSIiXlHmluyEBHP7+Wez6bdoIUmn0yyDnCmFHGvTBmPIEILGjDEXIfSDL+tKUcOuiIiIl9Svb464PPusOSKTmgqPPcaBVq1wAYEZGQQ9/bS5AGF4ONx1F6xcCceOeTty71LZSEREpAaw2eCqq3DeeSfhmZmEAXcC/8FcAZsDB+DFF+HWW83yUv/+8PTT5srZ/saPyka+f4YiImJ5RRPkHQReBn4DhACp8+fD1KnQpg2cOmUuV1D0+Ior4P77YdMmnLt3k5yc7NtzyahsJCIiUnOUN0Gey+EgZNSos7P3pqXBggXQrx/UqmXOMzN/PsTG0qBVK3Lj4pgZFcWrTz/tpbPwMD9amFHJi4iI1HiVmiCvbVtzJt8PPzRn+l2+nJ9uvZVDQGPgt8ArhsFvp06loHt3mDsXdu0CX5nuzI9GXnS3kYiIWMJ5V7wuT6NG8JvfsO2yy7jhzTfpAQzBXEzyKsCRkgIpKfDAA9CixdkVsfv2NW/VtiI17IqIiNQ87k6QFxMTA3Y7W4E/A1cDLe12fnz8cXPl64AAyMyEhQth0CAICYERI8zZavft89yJeIIadkVERKyvvHLTnxMTafzgg/Dee3D4MLzzDkycCM2bw08/wVtvmdPsN28O114Ls2fD9u1nRzZqKpWNREREfEOF5aYGDWDYMBg2DGd2NvvWrKFtejrBmzfDtm1m0rJ9O8yaBWFh5vwzQ4eaC002bOi1cyqXH5WNlLyIiIjPKzOb7zmSkpKYOHEiLpcLu91OYmIi499+G9auNWf5Xb8ecnPNclJSEtSpY97VVNQr06JF9Z3M+ahsJCIi4h+cTmdx4gLgcrlISEjAefo0/P738OabcOgQrF/Psbvu4kR4OJw8Ce+/D3/8I7RsCR07ms2///sfnD7tnRPxo5EXJS8iIuLXiibAK6mwsJCMjIyzOwICSMrKotFLL1F/3z462Gx8OnKkeXeSw2Hecj13LvTpA82awe9+B//+N/z4Y/WdiB+NvKhsJCIifq1oArySCYzD4aBNmzbFj88dnfnGMPjVypVkZmYS0aCBOQqzejWsWWMmLK+9Zm4OB/zqV2fLS+3bm0seeIIfNez6/hmKiIhUoDIT4FU4OtO4MYweDa++aq61tHkzzJgBV15pjoZ89BFMnw4dOpjLFvzpT7BhAxQUVO2JqGwkIiLiP8aPH09mZibJyclkZmYyfvz4Us+XtzzBuaMzgLkswXXXwd/+Bjt3wg8/wL/+xS+//jWu2rXNx//8JwwYYM4pc+ut8MILZjMw5gjPRa/B5EdlI4+e4fDhw4mKiqJu3bqEh4czduxYcnJyKnzNuHHjsNlspbaePXt6MkwREZEKJ8Cr1PIE5WnZkqR69Wjwv/8RdOoUt9hspF13nXnb9fHjsHIljB8P4eEcbNmSJVFR/F9cHNFRUSQlJbl3AiobVY1+/fqxYsUK0tLSePPNN/n++++57bbbLvi6QYMGsW/fvuJtzZo1ngxTRETkgi40OlOekr0yPwGrDIMrt27F+emn5vIEjzwCXbsCcFlmJrMNg8+ALMPAFR/PoRdeMCfOqww/WpjRow279957b/Hfo6OjeeCBB7j55ps5deoUtWvXPu/rAgICCAsLq9RnFBQUUFCibpifn3/xAYuIiFTgQvPFnOu8vTI//EBEbCx06wazZvHxG2/wwsiRDAP6A5cD8YZhjspMngxxcWbD75AhEB1d/odp5KXqHTlyhGXLltG7d+8KExeAjRs30qxZM9q2bUt8fDwHDhw477Fz5swhODi4eIuMjKzq0EVERC5KZXtlonv25CW7nRFAU2AgsNBm43REhNnYu3Yt/OEP0KIFp9q3Z8+YMRxYtersaAuoYbcqzZgxgwYNGtC0aVOysrJ4++23Kzx+8ODBLFu2jA8//JAFCxaQkpJCXFxcqdGVkmbOnEleXl7xlp2d7YnTEBERcVtle2VKHlcAfOBwUG/JEmplZZmNv3/7G1x3HS6bjdppaUT/+980u+UWfmnUCMaOheXL+flMT2nesWPVfJbVz2YYhuHOC2bNmsXs2bMrPCYlJYVu3boBcOjQIY4cOcKePXuYPXs2wcHBrF69Glsl73Pft28f0dHRvP7669xyyy0XPD4/P5/g4GDy8vIICgqq1GeIiIh4ktPpLH9tJTeOczqdXBMVxQDDYAgwGGhczns8ZLPRasmSSvXk1CTufH+7nbwcOnSIQ4cOVXhMixYtqFu3bpn9TqeTyMhItmzZQq9evSr9mTExMUyYMIEZM2Zc8FglLyIi4ouSk5OJi4srfuwAegOJw4fjeucdOpzZ/ztgucNhTqDnRn+Ot7nz/e12w25ISAghISEXFVhRnnS+ElB5Dh8+THZ2NuHh4Rf1mSIiIr7g3JmAC4EtDgdf/u53jHrnHVoBUcBmzk6gZ6XkxR0e63nZtm0bCxcuJDU1lT179pCcnMyYMWNo3bp1qVGX9u3bs2rVKgCOHz/OtGnT2Lp1K5mZmWzcuJFhw4YREhLCiBEjPBWqiIhIjXe+/pnevXtjt9v5AdiImdSUO4GeD/FY8lKvXj1WrlzJ9ddfT7t27bjrrrvo2LEjmzZtIiAgoPi4tLQ08vLyAPMf+6uvvuKmm26ibdu23HnnnbRt25atW7cSGBjoqVBFREQsoby5Zi56Aj0Lc7vnpaZTz4uIiPijyjYF11Qe7XkRERGRmsfdCfSszPen4RMRERGfouRFRERELEXJi4iIiFiKkhcRERGxFCUvIiIiYilKXkRERMRSlLyIiIiIpSh5EREREUtR8iIiIiKWouRFRERELEXJi4iIiFiKz61tVLTOZH5+vpcjERERkcoq+t6uzHrRPpe8HDt2DIDIyEgvRyIiIiLuOnbsGMHBwRUeYzMqk+JYiMvlIicnh8DAQGw2W7nHXHvttaSkpJz3PSp63t3nzt2Xn59PZGQk2dnZF1zy25Mu9G9QHe9V2ddV5riLuWaV3V8TrllVXq+LfT93XqOfMf2MubNf18z911zKz1hFz3vzZ8wwDI4dO0bz5s2x2yvuavG5kRe73X7BJcEdDkeF/+AVPe/uc+c7PigoyKs/pBf6N6iO96rs6ypz3MVcM3f3e/OaVeX1utj3c+c1+hnTz9jF7Pf3a1ZdP2MVPe/tn7ELjbgU8cuG3T/84Q8X/by7z13os7ylKuO62Peq7Osqc9zFXDN393tTVcd0Me/nzmv0M6afsYvZ723evmbV9TNW0fNW+RnzubJRTZefn09wcDB5eXle/Q1DKk/XzFp0vaxH18xaasL18suRF28KCAjgkUceISAgwNuhSCXpmlmLrpf16JpZS024Xhp5EREREUvRyIuIiIhYipIXERERsRQlLyIiImIpSl5ERETEUpS8iIiIiKUoeanhfv75Z6Kjo5k2bZq3Q5ELOHbsGNdeey1XX301nTp1YsmSJd4OSS4gOzub2NhYOnToQOfOnfnPf/7j7ZDkAkaMGEHjxo257bbbvB2KnMfq1atp164dMTExPP/88x75DN0qXcM99NBDpKenExUVxfz5870djlSgsLCQgoIC6tevz88//0zHjh1JSUmhadOm3g5NzmPfvn3s37+fq6++mgMHDtClSxfS0tJo0KCBt0OT80hOTub48eMsXbqUN954w9vhyDlOnz5Nhw4dSE5OJigoiC5duvDpp5/SpEmTKv0cjbzUYOnp6Xz77bfceOON3g5FKsHhcFC/fn0AfvnlFwoLCyu1tLt4T3h4OFdffTUAzZo1o0mTJhw5csS7QUmF+vXrR2BgoLfDkPPYtm0bV155JZdffjmBgYHceOONvP/++1X+OUpeLtJHH33EsGHDaN68OTabjbfeeqvMMYsWLaJly5bUrVuXrl27snnzZrc+Y9q0acyZM6eKIpbquGZHjx7lqquuIiIigunTpxMSElJF0fun6rhmRbZv347L5SIyMvISo/Zf1Xm9xDMu9Rrm5ORw+eWXFz+OiIhg7969VR6nkpeL9NNPP3HVVVexcOHCcp9fvnw5U6dO5aGHHuLzzz+nT58+DB48mKysrOJjunbtSseOHctsOTk5vP3227Rt25a2bdtW1yn5PE9fM4BGjRrxxRdfsHv3bl577TX2799fLefmq6rjmgEcPnyYO+64g8TERI+fky+rruslnnOp17C80WabzVb1gRpyyQBj1apVpfZ1797dmDRpUql97du3Nx544IFKvecDDzxgREREGNHR0UbTpk2NoKAgY/bs2VUVst/zxDU716RJk4wVK1ZcbIhyDk9ds19++cXo06eP8fLLL1dFmHKGJ3/GkpOTjVtvvfVSQ5QLuJhr+PHHHxs333xz8XP33HOPsWzZsiqPTSMvHnDy5El27NjBgAEDSu0fMGAAW7ZsqdR7zJkzh+zsbDIzM5k/fz7x8fE8/PDDnghXqJprtn//fvLz8wFz1dWPPvqIdu3aVXmsYqqKa2YYBuPGjSMuLo6xY8d6Ikw5oyqul3hXZa5h9+7d2blzJ3v37uXYsWOsWbOGgQMHVnkstar8HYVDhw5RWFhIaGhoqf2hoaHk5uZ6KSqpSFVcM6fTyfjx4zEMA8MwmDJlCp07d/ZEuELVXLOPP/6Y5cuX07lz5+La/iuvvEKnTp2qOly/V1X/Lw4cOJDPPvuMn376iYiICFatWsW1115b1eFKOSpzDWvVqsWCBQvo168fLpeL6dOne+SOSyUvHnRunc8wjIuq/Y0bN66KIpILuZRr1rVrV1JTUz0QlVTkUq7Zddddh8vl8kRYch6X+v+iJ+5cEfdc6BoOHz6c4cOHezQGlY08ICQkBIfDUea3iQMHDpTJWKVm0DWzHl0za9H1sr6adA2VvHhAnTp16Nq1Kxs2bCi1f8OGDfTu3dtLUUlFdM2sR9fMWnS9rK8mXUOVjS7S8ePHycjIKH68e/duUlNTadKkCVFRUdx3332MHTuWbt260atXLxITE8nKymLSpElejNq/6ZpZj66Zteh6WZ9lrmGV37/kJ5KTkw2gzHbnnXcWH/PMM88Y0dHRRp06dYwuXboYmzZt8l7AomtmQbpm1qLrZX1WuYZa20hEREQsRT0vIiIiYilKXkRERMRSlLyIiIiIpSh5EREREUtR8iIiIiKWouRFRERELEXJi4iIiFiKkhcRERGxFCUvIiIiYilKXkRERMRSlLyIiIiIpSh5EREREUv5f4vLYvk8s0H5AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hm = w.headinside(tobs2)\n",
    "plt.semilogx(tobs2, hobs2, \".k\")\n",
    "plt.semilogx(tobs2, hm[0], \"r\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "artesia",
   "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.10"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}