main_case12.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import flopy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "# ! Custom modules\n",
    "from utils import pht3d_fsp\n",
    "from utils.tools import calculate_MAS_sums, clear_directory\n",
    "from utils.pht3d_utils import copy_pht3d_datab, run_pht3d_program\n",
    "\n",
    "clear_directory()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_folder = 'case12'\n",
    "model_name = 'case12_modflow'\n",
    "model_ws = './models_folder/' + data_folder\n",
    "\n",
    "Lx = 0.212\n",
    "Ly = 1.0\n",
    "\n",
    "nrow = 1\n",
    "ncol = 212\n",
    "nlay = 1\n",
    "\n",
    "delr = Lx / ncol\n",
    "delc = Ly / nrow\n",
    "\n",
    "top = 0.0\n",
    "botm = -0.00038\n",
    "\n",
    "laytyp = 1\n",
    "\n",
    "hk = 1.0\n",
    "\n",
    "nper = 2\n",
    "perlen = [0.08333, 1.03939]\n",
    "nstp = [260, 3132]\n",
    "steady = [True, True]\n",
    "\n",
    "model = flopy.modflow.Modflow(modelname=model_name, model_ws=model_ws, exe_name='./bin/mf2005')\n",
    "dis = flopy.modflow.ModflowDis(model, nlay=nlay, nrow=nrow, ncol=ncol, \n",
    "                               delr=delr, delc=delc, top=top, botm=botm,\n",
    "                               nper=nper, perlen=perlen, nstp=nstp, steady=steady)\n",
    "\n",
    "ibound = np.ones((nlay, nrow, ncol), dtype=np.int32)\n",
    "ibound[0, 0, -1] = -1\n",
    "strt = np.zeros((nlay, nrow, ncol), dtype=np.float32)\n",
    "bas = flopy.modflow.ModflowBas(model, ibound=ibound, strt=strt)\n",
    "\n",
    "lpf = flopy.modflow.ModflowLpf(model, hk=hk, laytyp=laytyp, ipakcb=53, ss=0)\n",
    "\n",
    "chd_spd = {0: [[0, 0, ncol-1, 0.0, 0.0],]}\n",
    "chd = flopy.modflow.ModflowChd(model, stress_period_data=chd_spd)\n",
    "\n",
    "wel_spd = {0: [[0, 0, 0, 0.00048],],}\n",
    "wel = flopy.modflow.ModflowWel(model, ipakcb=66, stress_period_data=wel_spd, unitnumber=78)\n",
    "\n",
    "lmt = flopy.modflow.ModflowLmt(model, output_file_name='mt3d_link.ftl')\n",
    "\n",
    "stress_period_data = {}\n",
    "for kper in range(nper):\n",
    "    for kstp in range(nstp[kper]):\n",
    "        stress_period_data[(kper, kstp)] = [\"save head\", \"save drawdown\", \n",
    "                                            \"save budget\", \"print head\", \"print budget\",]\n",
    "oc = flopy.modflow.ModflowOc(model, stress_period_data=stress_period_data)#, compact=True)\n",
    "\n",
    "pcg = flopy.modflow.ModflowPcg(model=model)\n",
    "\n",
    "model.write_input()\n",
    "success, mfoutput = model.run_model(silent=True, pause=False, report=True)\n",
    "if not success:\n",
    "    raise Exception('MODFLOW did not terminate normally.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R2 score:  1.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAADaCAYAAADZqVrgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5HUlEQVR4nO3deVxUdfv/8dfMMKwCLmwiiKi5KygmgiBatyTuW+KSS2qpuWDaXVl3afX7ZtltBe6ZVpZb5ZJ7YbmA4IaouOSuoIK4AiL7nN8fJHekFiDDMMP1fDzm8YgzZ2beM/HB6zoz5xqVoigKQgghhKjS1IYOIIQQQgjDk4JACCGEEFIQCCGEEEIKAiGEEEIgBYEQQgghkIJACCGEEEhBIIQQQgikIBBCCCEEUhAIIYQQAikIqryZM2eiUqm4efOmoaMA/8sjhPifX3/9lbZt22JjY4NKpWLDhg0ArFmzhubNm2NlZYVKpeLIkSNlXkMjR46kXr16ZcpX0tt26tSJFi1alOkx9KVTp0506tTJ0DEqBSkIhBCiElMUhYEDB6LVatm4cSOxsbEEBQVx48YNhg0bRoMGDdi+fTuxsbE0atSIMWPGEBsbW+rHeeedd1i/fr0enoEwFmaGDiCEEOLxrl27xu3bt+nbty/PPvts0fa9e/eSl5fHCy+8QFBQUNF2a2tr3NzcSv04DRo0KJe8wnjJEQIBwPXr1xk8eDD29vY4OzszatQo0tLSiq5XFIUFCxbg7e2NlZUVNWrUYMCAAVy4cKHY/URGRtK7d2/c3NywtLSkYcOGjB079pFvSWzZsgVvb28sLCzw9PTkv//9r96fpxCVSXR0NM8++yy2trZYW1vj7+/Pli1biq6fOXNm0T/ub7zxBiqVinr16jFy5EgCAgIACA0NRaVSFR32ftxbBitXrsTPz49q1apRrVo1vL29Wbp0adH1jzrsP3/+fDp27IiTkxM2Nja0bNmS2bNnk5eX90TP++DBgwQGBmJtbU39+vX56KOP0Ol0xfZJT0/ntddew9PTE3Nzc+rUqcOUKVPIzMwsU0ZFUZg9ezYeHh5YWlrSpk0btm3b9kTPw9TIEQIBQP/+/QkNDWX06NEkJCQwffp0AJYtWwbA2LFj+frrr5k8eTIff/wxt2/f5v3338ff35+jR4/i7OwMwPnz5/Hz82PMmDHY29tz6dIlPv30UwICAkhISECr1QKF74n27t0bPz8/Vq9eTUFBAbNnz+b69euGeQGEqGC7d++mS5cutGrViqVLl2JhYcGCBQvo2bMnq1atIjQ0lDFjxuDl5UW/fv2YNGkSQ4YMwcLCAjs7O9q1a8eECRP48MMP6dy5M3Z2do99rHfffZcPPviAfv36MW3aNOzt7Tl+/DiXL1/+24znz59nyJAhRf8oHz16lP/7v//j999/L/rbUFopKSkMHTqUadOmMWPGDNavX8/06dNxdXVl+PDhANy/f5+goCCuXLnCW2+9RatWrThx4gTvvvsuCQkJ7Nixo6joKWnG9957j/fee4/Ro0czYMAAkpKSeOmllygoKKBx48Zlei4mRxFV2owZMxRAmT17drHtr7zyimJpaanodDolNjZWAZQ5c+YU2ycpKUmxsrJSXn/99Ufet06nU/Ly8pTLly8rgPLTTz8VXefr66u4uroqWVlZRdvS09OVmjVrKvJrKaqC9u3bK05OTkpGRkbRtvz8fKVFixaKm5ubotPpFEVRlIsXLyqA8sknnxS7/c6dOxVA+eGHH4ptf7CmH7hw4YKi0WiUoUOH/m2eESNGKB4eHo+9vqCgQMnLy1OWL1+uaDQa5fbt2yW+7QNBQUEKoOzfv7/Y9mbNminPPfdc0c+zZs1S1Gq1cvDgwWL7/fjjjwqgbN26tVQZ79y5o1haWip9+/Yttv/evXsVQAkKCvrH7FWBvGUgAOjVq1exn1u1akV2djapqals3rwZlUrFCy+8QH5+ftHFxcUFLy8vdu3aVXS71NRUxo0bh7u7O2ZmZmi1Wjw8PAA4deoUAJmZmRw8eJB+/fphaWlZdFtbW1t69uyp/ycrhIFlZmayf/9+BgwYQLVq1Yq2azQahg0bxpUrVzh9+nS5PFZkZCQFBQVMmDCh1LeNj4+nV69e1KpVC41Gg1arZfjw4RQUFHDmzJky5XFxcaFdu3bFtrVq1arY0YrNmzfTokULvL29i/3Nee6551CpVMX+5pQkY2xsLNnZ2QwdOrTY4/r7+xf9fRLyloH4Q61atYr9bGFhAUBWVhbXr19HUZSitwX+qn79+gDodDqCg4O5du0a77zzDi1btsTGxgadTkf79u3JysoC4M6dO+h0OlxcXB66r0dtE8LU3LlzB0VRqF279kPXubq6AnDr1q1yeawbN24AlPqDhomJiQQGBtK4cWPCw8OpV68elpaWHDhwgAkTJhSt59L6698aKPx78+f7u379OufOnSt6i/GvHnwmqaQZH7yW8jfn70lBIP6Rg4MDKpWKqKiookLhzx5sO378OEePHuXrr79mxIgRRdefO3eu2P41atRApVKRkpLy0H09apsQpqZGjRqo1WqSk5Mfuu7atWtA4borD46OjgBcuXIFd3f3Et9uw4YNZGZmsm7dumJd9JEjR8ol199xcHDAysrqsZ9TePDalDTjgyLkcX9zyjp/wdTIWwbiH/Xo0QNFUbh69Spt27Z96NKyZUuAog/5/LVoWLx4cbGfbWxsaNeuHevWrSM7O7toe0ZGBps2bdLzsxHC8GxsbPD19WXdunXFOmOdTsd3332Hm5sbjRo1KpfHCg4ORqPRsHDhwlLd7lHrWVEUlixZUi65/k6PHj04f/48tWrVeuTfnAf/gJc0Y/v27bG0tGTFihXFtsfExPzjByurEjlCIP5Rhw4dePnll3nxxRc5dOgQHTt2xMbGhuTkZKKjo2nZsiXjx4+nSZMmNGjQgDfffBNFUahZsyabNm0iMjLyofv84IMP6Nq1K126dGHatGkUFBTw8ccfY2Njw+3btw3wLIWoWLNmzaJLly507tyZ1157DXNzcxYsWMDx48dZtWpVuU3srFevHm+99RYffPABWVlZRacXnzx5kps3b/Lee+898nZdunTB3NycwYMH8/rrr5Odnc3ChQu5c+dOueT6O1OmTGHt2rV07NiRV199lVatWqHT6UhMTOSXX35h2rRp+Pr6ljhjjRo1eO211/h//+//MWbMGJ5//nmSkpKYOXOmvGXwJ3KEQJTI4sWLmTdvHnv27GHQoEF0796dd999l8zMzKIPCGm1WjZt2kSjRo0YO3YsgwcPJjU1lR07djx0f126dGHDhg2kp6cTGhrK1KlT6d+/P6NGjaropyaEQQQFBfHbb79hY2PDyJEjGTRoEGlpaWzcuJHQ0NByfaz333+f5cuXc/nyZYYOHUqfPn346quv8PT0fOxtmjRpwtq1a7lz507RaY/e3t5ERESUa7ZHsbGxISoqipEjR/LFF1/QvXt3Bg4cSEREBG5ubkVHCEqT8f3332fWrFn88ssv9OrVi7lz57Jo0SI55fBPVIqiKIYOIYQQQgjDkiMEQgghhJCCQAghhBBSEAghhBACKQiEEEIIgRQEQgghhEAKAiGEEEJgJIOJdDod165dw9bWttyGdQhRFSmKQkZGBq6urqjVlbcfkDUvRPkp6bo3ioLg2rVrpZrBLYT4e0lJSaX+spuKJGteiPL3T+veKAoCW1tboPDJ2NnZGTiNEMYrPT0dd3f3ojVVWcmaF6L8lHTdG0VB8OCQoZ2dnfxxEKIcVPbD8LLmhSh//7TuK++biEIIIYSoMFIQCCGEEMK0CoLsvAJDRxBCVKDsvAJ0Ovl+NiHKg8kUBLvP3KDTJ7uIOXfT0FGEEBXkrXUJjP7mIDfv5Rg6ihBGz2QKgkW7zpOSns3QpfuZvf138gp0ho4khNCjSzcz2ZyQzM7TNwgJjyL6rDQDQjwJkykIlo5sy+B27igKLNh1noGLY0m6fd/QsYQQelLPwYaNEzvwlFM1bmTkMGzZfj6WZkCIMjOZgsDa3IxZ/Voxf0gbbC3NiE+8S7eIKDYfu2boaEIIPWniYsfGiQEM8a2LosDCXed5fpE0A0KUhckUBA90b1WbrZMDaVO3OhnZ+UxcGc+ba49xPzff0NGEEHpgZa7hw74tWTC0DXaWZhxJuku38Cg2HpVmQIjSMLmCAMC9pjXfj/VjYueGqFSw+mASPedGcyo53dDRhBB60q1lbbaGBdLWowYZOflMXhXP6z8elWZAiBIyyYIAwEyj5rXnGrNitC9Othacv5FJ7/l7WR57CUWR05SEMEVuNaxZ/XJ7Jj9T2Ax8f+gKPeZGc+JamqGjCVHpmWxB8IB/Qwe2T+nIs02cyM3X8e5PJ3hpeRx3MnMNHU0IoQdmGjVTgxuzckx7nO0suHAjk77zY/hq70VpBoT4GyZfEADUtDHnyxFtmdGzGeYaNTtOXSckPIp9F24ZOpoQQk/8GtRiW1hH/tXUidwCHe9tOslLyw9xW5oBIR6pShQEUPilDi928GTdK/7Ud7AhJT2bIUv28WnkGfLlNCUhTFJNG3OWDG/LzKJmIJWQ8D3EnJeZBUL8VZUpCB5oUceeTZMCGNjWDZ0CEb+eZdAX+7h6N8vQ0YQQeqBSqRjZwZMNEzrQwNGG6+k5DP1yP3N+OS3NgBB/UuUKAgAbCzNmD/AifJA31SzMOHT5DiGf72FbQrKhowkh9KSZqx2bJgUQ2rZwgNnc384R+sU+rtyRmQVCQBUtCB7o7V2HrZMD8XKvTnp2PuNXHOat9QnyJUlCmChrczM+HtCKuYNbY2thRtzlO4SER7FVmgEhqnZBAFC3ljU/jvNjXFADAFbuT6TXvGhOp2QYOJkQQl96ermyNSyQ1n8MMHtlxWGmr0sgK1eaAVF1VfmCAECrUfNmSBO+Hd0OR1sLzly/R6950Xy377KcpiSEiXowwOyVTg1QqWDVgcJm4PcUGWAmqiYpCP4k8ClHtoUFEtTIkZx8Hf/ZcJzx3x3m7n05TUkIU6TVqHm9axO+G+2Lo60FZ1Pv0WveXr6VAWaiCpKC4C8cqlnw1cin+U/3pmg1KrafSKFbeBQHL902dDQhhJ50aOjA9rBAOjd2JDdfxzs/nWDcd3HSDIgqRQqCR1CrVYwJrM+68R2oV8uaa2nZhC6OJXzHWQp00jUIYYpqVbNg2cineadHM7QaFT+fKBxgtl8GmIkqQgqCv9HSzZ7NkwPp16YOOgU+23GGwUv2kZwmMwuEMEUqlYrRAZ6sf6UDng42JKdlM3jJPj6TAWaiCpCC4B9UszDj04HefBbqhY25hgMXbxMSHsUvJ1IMHU0IoSct6tizeVIA/dsUDjAL//UsQ5bs55oMMBMmTAqCEurb2o0tkwNpWceeu/fzePnbOGb8dFxmFghhomwszJgz0IvPQwsHmB24VNgMbD8uzYAwTVIQlEI9BxvWjvfn5Y71Afgm9jJ95u/lXKrMLBDCVPVpXYctkwPwcrMnLSuPcd/F8Z8NMsBMmB4pCErJ3EzNW92a8vWLT+NQzZzfUzLoMTea1QcS5TQlIUyURy0bfhjnz9g/moHv9iXSZ/5ezl6XZkCYDikIyqhTYye2hgUS+JQD2Xk63lyXwMRV8aRl5Rk6mhBCD8zN1Ezv1pTlo9oVNQM950Wzcr80A8I0SEHwBJxsLfnmxXZMD2mCmVrFlmPJdAuPIu7yHUNHE0LoScdGjmwL60jHRo5k5+l4a30CE1YeJu2+NAPCuElB8ITUahVjgxrw43h/6ta05urdLAYujmX+znMys0AIE+Voa8HXI5/mrW6FzcDWhBS6RUQRd1kGmAnjJQVBOfF2r86WyQH08nKlQKfwyc+neeHL/VxPzzZ0NCGEHqjVKl7u2IC14/3xqPWgGdjH3F9lgJkwTlIQlCNbSy3hg7z5ZEArrM01xF64RdfP9/DrqeuGjiaE0BMv9+psnhRAH+/CZmBO5BmGfrmPlDRpBoRxkYKgnKlUKp5v686mSQE0d7Xjzv08Rn9ziPc2nSAnX05TEsIU2Vpq+XxQa+Y874W1uYZ9F24TEr6HHSelGRDGQwoCPWngWI11r/jzYod6AHy19xJ958dw/sY9wwYTQuhNfx83Nv+pGRiz/BAzN56QmQXCKEhBoEcWZhpm9GzO0hFtqWljzsnkdHrOjeaHQ0lympIQJqr+H83AmABPAL6OuUTfBTGcS5VmQFRuUhBUgGebOrMtLBD/BrW4n1vAv388RtjqI2Rky2lKQpgiCzMN/+nRjK9GPk0tG3NO/dEMfH9QmgFReZWpIFiwYAGenp5YWlri4+NDVFTUY/ddt24dXbp0wdHRETs7O/z8/Pj555/LHNhYOdtZ8u1oX/79XGM0ahUbj16jW0QUR5LuGjqaECUi6770OjdxYltYIB0a1iIrr4DX1x5j0qp40qUZEJVQqQuCNWvWMGXKFN5++23i4+MJDAwkJCSExMTER+6/Z88eunTpwtatW4mLi6Nz58707NmT+Pj4Jw5vbDRqFRM6N+T7sX7UqW5F0u0sBiyMYdHu8+jkNCVRicm6LzsnO0u+HeXLG10LZxZsPpZM94go4hNlgJmoXFRKKY9f+fr60qZNGxYuXFi0rWnTpvTp04dZs2aV6D6aN29OaGgo7777bon2T09Px97enrS0NOzs7EoTt9JKy8rjrfUJbDmWDEDgUw7MGeiFk62lgZMJU1bWtVTR694U1zzA4cQ7TF4Vz5U7WZipVUwNbsS4jg1Qq1WGjiZMWEnXU6mOEOTm5hIXF0dwcHCx7cHBwcTExJToPnQ6HRkZGdSsWfOx++Tk5JCenl7sYmrsrbTMG9yaj/u3xFKrJursTUI+j2LX6VRDRxOimIpY91VhzQO0qVuDrWGB9GhVm3ydwuztpxm2bD+pMsBMVAKlKghu3rxJQUEBzs7OxbY7OzuTklKy7wifM2cOmZmZDBw48LH7zJo1C3t7+6KLu7t7aWIaDZVKRejTddk8KYAmLrbcysxl5FcH+b8tJ8nN1xk6nhBAxaz7qrLmAewstcwd3JrZ/VthpdWw99wtQsKj2CnNgDCwMn2oUKUqfnhLUZSHtj3KqlWrmDlzJmvWrMHJyemx+02fPp20tLSiS1JSUlliGo2GTrZsmNCBEX4eACyJukj/hTFcvJlp4GRC/I8+131VW/MqlYqBTxcOMGta245bmbm8+NVBPth8UgaYCYMpVUHg4OCARqN5qCtITU19qHv4qzVr1jB69Gi+//57/vWvf/3tvhYWFtjZ2RW7mDpLrYb3erfgi2E+VLfWknA1jR4RUaw7fMXQ0UQVVxHrviqueYCGTtVY/4o/I/3rAbA0urAZuCADzIQBlKogMDc3x8fHh8jIyGLbIyMj8ff3f+ztVq1axciRI1m5ciXdu3cvW9IqIri5C9vCAvH1rElmbgFTvz/K1DVHuJeTb+hoooqSda9flloNM3s158vhbalhreX41XR6zI1mbdwVmVkgKlSp3zKYOnUqX375JcuWLePUqVO8+uqrJCYmMm7cOKDw0N/w4cOL9l+1ahXDhw9nzpw5tG/fnpSUFFJSUkhLSyu/Z2FiattbsfKl9kzt0gi1CtbFX6VHRBQJV+Q1E4Yh617//tXMmW1hHWlfvyb3cwuY9sNRXl0jA8xEBVLKYP78+YqHh4dibm6utGnTRtm9e3fRdSNGjFCCgoKKfg4KClKAhy4jRowo8eOlpaUpgJKWllaWuEbtwMVbit+HOxSPNzYrDd/aonyx+7xSUKAzdCxhpJ5kLVXkuq/Kaz6/QKfM/fWMUn/6FsXjjc1Kx9m/KUcS7xg6ljBiJV1PpZ5DYAimek5ySaXdz+ONtcfYfqLwPdygRo7893kvHG0tDJxMGBtjWUvGklOf4i7fZvKqI1y9Wziz4PWujRkTUF9mFohS08scAmEY9tZaFr7Qhv/r2wILMzW7z9wgJDyKqLM3DB1NCKEnPh412To5kJAWLuTrFD7c+jsjvz7IjYwcQ0cTJkoKAiOhUqkY6uvBxokBNHKuxs17OQxbeoBZ206RVyAzC4QwRfbWWhYMbcOsfoUDzPacuUFI+B72nJFmQJQ/KQiMTGMXWzZODGCob10AFu++wIBFsSTeum/gZEIIfVCpVAxuV5dNEwNo7GzLzXu5DF92gFlbT8kAM1GupCAwQpZaDf/XtyWLXmiDnaUZR5Pu0i0iip+OXDV0NCGEnjzlbMtPEzswrH3hALPFey7w/KIYLt+SAWaifEhBYMS6tqjNtikdebpeDe7l5BO2+gj//uEomTKzQAiTZKnV8EGfFiwe5oO9lZajV9LoHhHNhnhpBsSTk4LAyNWpbsWql9oz+dmnUKvgh7gr9JwXzfGrcr63EKbquT8GmLWrV5N7OflMWXOEad9LMyCejBQEJsBMo2Zql0asfKk9LnaWXLiRSb8FMSyLviiTzoQwUa7VrVj5ki9T/lXYDKw9fIUec6UZEGUnBYEJaV+/FtvCAunSzJncAh3vbz7J6G8OceuenKYkhCky06iZ8q9GrH7Zj9r2lly8mUnfBXtZKs2AKAMpCExMDRtzvhjmw/u9m2Nupua331MJCY8i5txNQ0cTQuhJO8+abAsLJLiZM3kFCh9sPsmorw9KMyBKRQoCE6RSqRjuV4+fJnSgoVM1UjNyGLp0P5/8/LvMLBDCRFW3NmfxMB8+6NMCczM1O0/foGt4FHulGRAlJAWBCWta246NEzswuJ07igLzd54ndHEsSbdlZoEQpkilUjGsvQcbJ3bgKadq3MjI4YWl+5m9XZoB8c+kIDBx1uZmzOrXinlDWmNracbhxMKZBVuOJRs6mhBCT5q42LFxYgCD29VFUWDBrvM8v0iaAfH3pCCoInq0cmXr5EBa161ORnY+E1Ye5s21x8jKLTB0NCGEHliZa5jVryULhhYOMDuSdJdu4VFsOnrN0NFEJSUFQRXiXtOa78f6MaFzA1QqWH0wiZ7zojmVnG7oaEIIPenWsjZbwwLx8ahBRk4+k1bF88aPx7ifKzMLRHFSEFQxWo2afz/XhBWjfXGyteBc6j16z9/L8thLcpqSECbKrYY1a15uz6RnGqJSwZpDSfSYG82JazKzQPyPFARVlH9DB7aFBfJMEydy83W8+9MJXv42jjuZuYaOJoTQAzONmmnBjVkxxhdnOwsu3Mik7/wYvt4rMwtEISkIqrBa1SxYOqIt7/ZohrlGTeTJ63SLiGLfhVuGjiaE0BP/Bg5sC+vIv5o6kVugY+amk7y0/BC3pRmo8qQgqOJUKhWjAjxZ94o/9R1sSE7LZsiSfXwaeYZ8OU1JCJNU08acJcPbMrNnYTOw41QqIeF7iD0vzUBVJgWBAKBFHXs2TQrgeR83dApE/HqWwUv2cfVulqGjCSH0QKVSMbKDJxsmdKC+ow3X03MY8uU+5vxyWpqBKkoKAlHExsKMT573InyQN9UszDh46Q7dwqPYflxmFghhqpq52rF5UgAD27qhKDD3t3OEfrGPK3dkZkFVIwWBeEhv7zpsnRyIl3t10rLyGPfdYd5en0B2nswsEMIUWZubMXuAFxGDW2NrYUbc5cJmYFuCNANViRQE4pHq1rLmx3F+jAtqAMCK/Yn0mhfNmesZBk4mhNCXXl6ubA0LxNu9OunZ+YxfcZjp6xJkgFkVIQWBeCytRs2bIU34dnQ7HKpZcOb6PXrOjWbF/stympIQJsq9pjU/jPNjfKfCAWarDiTSe340p1OkGTB1UhCIfxT4lCPbpwQS1MiRnHwdb68/zisrDpN2P8/Q0YQQeqDVqHmjaxO+HeWLo21hM9BrXjTf7pNmwJRJQSBKxKGaBV+NfJr/dG+KVqNi2/EUQsL3cPDSbUNHE0LoScBTDmwPC6Rz48Jm4J0Nxxn3XRx378vMAlMkBYEoMbVaxZjA+qwd70+9WtZcS8smdHEs4TvOUqCTrkEIU1Q4wOx/zcDPJ64TEh7FgYvSDJgaKQhEqbVyq87myYH0a10HnQKf7TjDkCX7SE6TmQVCmKIHzcD6Vzrg+ccAs0FfxPL5jjPSDJgQKQhEmVSzMOPTUG8+HeiFjbmG/RdvExIexS8nUgwdTQihJw8GmPVvUzjA7PMdhQPMrskAM5MgBYF4Iv3auLF5ciAt69hz934eL38bx4yfjsvMAiFMVDULM+YM9OLzUG9szDUc+KMZ+FmaAaMnBYF4Yp4ONqwd789LgZ4AfBN7mT7z93IuVU5TEsJU9Wldhy2TA2nlZk9aVh5jv43jnQ3SDBgzKQhEuTA3U/N292Z89eLT1LIx5/eUDHrO3cuag4lympIQJqqegw0/jvNnbMf6AHy7r7AZOCsDzIySFASiXHVu7MS2KYEEPuVAVl4Bb6xNYOKqeNKyZGaBEKbI3EzN9G5N+WZUOxyq/dEMzItm5X5pBoyNFASi3DnZWvLNi+14M6QJZmoVW44l0z0iirjLdwwdTQihJ0GNHNkW1pHApxzIztPx1voEJq6UZsCYSEEg9EKtVjEuqAE/jPPDvaYVV+5kMXBxLPN3npPTlIQwUY62FnzzYjumP2gGEpLpFh5F3GWZWWAMpCAQetW6bg22TA6kl5crBTqFT34+zbCl+7menm3oaEIIPVCrVYwNasDa8f541LLm6t0sBi7ex7zfZIBZZScFgdA7O0st4YO8+WRAK6y0GmLO3yIkPIrffr9u6GhCCD3xcq/O5kkB9PEubAb++8sZXvhyPylp0gxUVlIQiAqhUql4vq07mycH0Ky2Hbczcxn19SHe33SSnHw5TUkIU2RrqeWzUG/++7wX1uYaYi/cIiR8D7+ekmagMpKCQFSoBo7VWD/Bnxc71ANg2d6L9FsQw4Ub9wwbTAihFyqVigE+bmyeFEBzVzvu3M9j9DeHmLnxhMwsqGSkIBAVzsJMw4yezVk6oi01bcw5cS2dHnOj+eFQkpymJISJqu9YjXWv+DM6oHCA2dcxl+i3IIZzqdIMVBZSEAiDebapM9vCAvGrX4v7uQX8+8djhK0+Qka2nKYkhCmyMNPwTo9mfDXyaWramHMyOZ2ec6P5/qA0A5WBFATCoJztLPlujC//fq4xGrWKjUev0T0imiNJdw0dTQihJ52bOLE9LJAODWuRlVfA62uPMXn1EdKlGTCoMhUECxYswNPTE0tLS3x8fIiKinrsvsnJyQwZMoTGjRujVquZMmVKWbMKE6VRq5jQuSHfj/WjTnUrEm/fZ8DCGBbtPo9OTlOqNGTdi/LkZGfJt6N8eb1rYTOw6eg1ukdEEZ8oA8wMpdQFwZo1a5gyZQpvv/028fHxBAYGEhISQmJi4iP3z8nJwdHRkbfffhsvL68nDixMl49HDbaGBdK9ZW3ydQofbfudEV8dIDVDTlMyNFn3Qh/UahWvdGrID+P8cKthRdLtLJ5fFMuCXeekGTAAlVLKN258fX1p06YNCxcuLNrWtGlT+vTpw6xZs/72tp06dcLb25vPP/+8VCHT09Oxt7cnLS0NOzu7Ut1WGB9FUVhzMImZm06QnafDoZo5cwZ6E9TI0dDRjF5Z11JFr3tZ81VPenYeb61LYPOxZAACGjrw6UAvnOwsDZzM+JV0PZXqCEFubi5xcXEEBwcX2x4cHExMTEzZkgrxFyqVikHt6rJ5UgBNXGy5eS+XEcsO8OHWU+Tm6wwdr8qRdS8qgp2llrmDWzO7f+EAs+hzNwkJj2Ln6VRDR6sySlUQ3Lx5k4KCApydnYttd3Z2JiUlpdxC5eTkkJ6eXuwiqp6GTrZsmNCBEX4eAHyx5wL9F8Zw6WamgZNVLRWx7mXNCyhsBgY+7c6mSR1oWtuOW5m5vPjVQT7YLAPMKkKZPlSoUqmK/awoykPbnsSsWbOwt7cvuri7u5fbfQvjYqnV8F7vFnwxzIfq1loSrqbRPSKK9fFXDB2tytHnupc1L/6soZMt61/xZ6R/PQCWRl+k/8IYLkozoFelKggcHBzQaDQPdQWpqakPdQ9PYvr06aSlpRVdkpKSyu2+hXEKbu7CtrBA2nnWJDO3gFfXHGXq90e4l5Nv6GgmryLWvax58VeWWg0zezXny+FtqWGt5fjVdLpHRLE2TpoBfSlVQWBubo6Pjw+RkZHFtkdGRuLv719uoSwsLLCzsyt2EaK2vRWrXmrP1C6NUKtg3eGr9IiIIuFKmqGjmbSKWPey5sXj/KuZM9vCOtK+fk3u5xYw7YejvLpGmgF9KPVbBlOnTuXLL79k2bJlnDp1ildffZXExETGjRsHFFb6w4cPL3abI0eOcOTIEe7du8eNGzc4cuQIJ0+eLJ9nIKoUjVrF5GefYs1YP1ztLbl06z79Fu7ly6gLcpqSHsm6F4bkYm/JijHtmdalERq1ivXxV+keEcWxK3cNHc2klPq0QygcUDJ79mySk5Np0aIFn332GR07dgRg5MiRXLp0iV27dv3vQR7xPqOHhweXLl0q0ePJKUjiUe7ez+WNtcf4+UThN6cFNXJkzkAvHKpZGDhZ5fUka6ki172sefE4hy7dJmz1Ea7ezcJMreL1ro0ZE1Aftbr8Psdmakq6nspUEFQ0+eMgHkdRFFbsT/zjU8g6HG0t+GygNwFPORg6WqVkLGvJWHIKw0i7n8eb646x7Xjh51o6NnJkzvNeONpKM/AoeplDIERlo1KpeKG9BxsnBtDIuRo3MnIYtmw/H237nbwCmVkghCmyt9ayYGgbPuzbEgszNXvO3CAkPIo9Z24YOppRk4JAmITGLrb8NCGAIb51URRYtPs8AxbFknjrvqGjCSH0QKVSMcS3LpsmBdDY2Zab93IYvuwAs2SAWZlJQSBMhpW5hg/7tmTh0DbYWZpxNOku3SOi2Hj0mqGjCSH0pJGzLT9N7MCw9oUDzBbvucDzi2K4fEtmFpSWFATC5IS0rM22KR1p61GDjJx8Jq+K598/HOV+rpymJIQpstRq+KBPCxa94IO9lZajV9LoHhHNT0euGjqaUZGCQJikOtWtWP1yeyY/0xCVCn6Iu0KPudGcuCYzC4QwVV1buLA1LJCn69XgXk4+YauPMO37o2TKzIISkYJAmCwzjZqpwY1ZOaY9LnaWXLiRSd/5MXy19yJGcHKNEKIM6lQvHGAW9uxTqFWw9vAVes6N5vhVaQb+iRQEwuT5NajFtrBA/tXUmdwCHe9tOsmYbw5xOzPX0NGEEHpgplHzapdGrHqpPbXtLblwM5N+C2JYGi3NwN+RgkBUCTVszFky3If3ezfH3EzNr7+n0vXzPcScv2noaEIIPfGtX4utkwMJblbYDHyw+SSjvj7IrXs5ho5WKUlBIKoMlUrFcL96bHilAw0cbUjNyGHol/v578+nyZeZBUKYpBo25iwe5sMHfzQDO08XziyIOSfNwF9JQSCqnGaudmyaFMCgp91RFJi38xwDF8eSdFtmFghhilQqFcP86rFxYgcaOlUrbAaW7mf2dhlg9mdSEIgqydrcjI/6t2LekNbYWppxOPEu3SKi2HIs2dDRhBB60sTFjk0TAxjcrnCA2YJd56UZ+BMpCESV1qOVK1snB9K6bnUysvOZsPIw09cdIyu3wNDRhBB6YGWuYVa/lswf0gZbSzPiE+/SLTyKTTLATAoCIdxrWvP9WD8mdG6ASgWrDiTRc140p5LTDR1NCKEn3VvVZltYID5/DDCbtCqeN348VqUHmElBIASg1aj593NN+G60L062FpxLvUfv+XtZHntJTlMSwkS51bBmzcvtmfTHALM1h5LoOTeak9eqZjMgBYEQf9KhoQPbwgJ5pokTufk63v3pBC9/G8cdmVkghEky06iZFtyYFWN8cbaz4PyNTPos2Ms3MVWvGZCCQIi/qFXNgqUj2vJuj2aYa9REnrxOt4go9l+4ZehoQgg98W/gwLawjjz7RzMwY+MJXlpetZoBKQiEeASVSsWoAE/WveJPfQcbktOyGbxkH59FnpGZBUKYqJo25nw5oi0zexY2AztOXSckPIrY81WjGZCCQIi/0aKOPZsmBTDAxw2dAuG/nmXIkv1cvZtl6GhCCD1QqVSM7ODJ+gn+1He0ISU9myFf7uPTX0x/gJkUBEL8AxsLM/77vBfhg7ypZmHGgUu36RYexfbjKYaOJoTQk+au9myeFMDAtm4oCkT8do5BX+zjyh3TnVkgBYEQJdTbuw5bJgfg5WZPWlYe476L4z8bEsjOk5kFQpgia3MzZg8obAZsLcw4dPkO3cKj2JZgmgPMpCAQohQ8atnwwzh/xgbVB+C7fYn0nreXM9czDJxMCKEvhc1AIN7u1UnPzmf8isO8tT7B5AaYSUEgRCmZm6mZHtKU5aPa4VDNgtPXM+g5N5oV+y9XudOUhKgq6tay5odxfozvVDjAbOX+RHrPj+Z0iuk0A1IQCFFGHRs5si0skI6NHMnJ1/H2+uO8suIwaffzDB1NCKEHWo2aN7o24dtRvjjaWnDm+j16zYvmu32m0QxIQSDEE3C0teDrkU/zdremaDUqth1PoVtEFIcu3TZ0NCGEngQ8VTjArFPjwmbgPxuOM/67w9y9b9wzC6QgEOIJqdUqXupYn7Xj/alXy5qrd7MYuDiWiF/PUqAz/q5BCPEwh2oWLBvxNP/pXtgMbD+RQrfwKA5cNN5mQAoCIcpJK7fqbJ4cSN/WddAp8GnkGYYs2UdymswsEMIUqdUqxgTWZ934DtSrZc21tGwGfRFL+A7jbAakIBCiHFWzMOOzUG8+HeiFtbmG/RdvExIeReTJ64aOJoTQk5Zu9myeHEi/NoXNwGc7zjB4yT6uGdkAMykIhNCDfm3c2DI5kBZ17Lh7P4+Xlh9ixk/HZWaBECaqmoUZnw705rNQL2zMNRy4eJtuEVH8csJ4BphJQSCEnng62LBufAfGBHgC8E3sZfrM38u5VNM5TUkIUVzf1oXNQCs3e+7ez+Plb+N410iaASkIhNAjczM1/+nRjK9efJpaNub8npJBz7l7WXMw0SROUxJCPKyegw0/jvPn5Y6FA8yW/9EMnK3kA8ykIBCiAnRu7MS2sEACGjqQlVfAG2sTmLQqnvRsmVkghCkyN1PzVremfDOqHQ7V/mgG5kWz6kDlbQakIBCigjjZWbJ8VDve6NoEM7WKzceS6RYexeHEO4aOJoTQk6BGjmwNCyTwKQey83RMX5fAxJXxpGVVvmZACgIhKpBarWJ8pwb8MM4P95pWXLmTxfOLYpm/8xw6IzxNSQjxz5xsLfnmxXZMDylsBrYkFDYDcZcr18wCKQiEMIDWdWuwZXIgPb1cKdApfPLzaYYt28/19GxDRxNC6IFarWJsUAN+HO9P3ZoPBpjtY95vlWdmgRQEQhiInaWWiEHezB7QCiuthr3nbhESHsXO31MNHU0IoSfe7tXZMjmA3t6FzcB/fznDC1/uJyXN8M2AFARCGJBKpWJgW3c2Tw6gWW07bmfm8uLXB3l/00ly8iv/aUpCiNKztdTyeag3/32+cIBZ7IVbhITv4ddThh1gJgWBEJVAA8dqrJ/gz4sd6gGwbO9F+i2I4cKNe4YNJoTQC5VKxQAfNzZPCqC5qx137ucx+ptDzNx4wmDNgBQEQlQSFmYaZvRsztIRbalhreXEtXR6zI3mx7grlfY0JSHEk6nvWI11r/gzqkPhALOvYy7Rd34M5w3QDEhBIEQl82xTZ7aFdcSvfi3u5xbw2g9HmbLmCBkys0AIk2RhpuHdns1YNrItNW3MOZmcTo+IaL4/lFShzYAUBEJUQi72lnw3xpd/P9cYjVrFT0eu0T0imiNJdw0dTQihJ880cWZbWCD+DWqRlVfA6z8eI2z1kQobYCYFgRCVlEatYkLnhnw/tj11qluRePs+AxbGsGj3eZlZIISJcraz5NvR/2sGNh69RveIKOIrYICZFARCVHI+HjXZGhZI95a1ydcpfLTtd0Z8dYDUDMOfpiSEKH8PmoEfxvnhVsOKpNuFA8wW7tJvM1CmgmDBggV4enpiaWmJj48PUVFRf7v/7t278fHxwdLSkvr167No0aIyhRWiqrK30jJvSGtm9WuJpVZN1NmbdAuPYveZGxWWQda9EBWrzR8DzLq3KmwGPt7+O8OXHSBVTwPMSl0QrFmzhilTpvD2228THx9PYGAgISEhJCYmPnL/ixcv0q1bNwIDA4mPj+ett95i8uTJrF279onDC1GVqFQqBrery6aJATRxseXmvVxGLDvAh1tPkZuv0+tjy7oXwjDsrbTMG9yaj/sXNgPR524WDjA7Xf4DzFRKKT/C6OvrS5s2bVi4cGHRtqZNm9KnTx9mzZr10P5vvPEGGzdu5NSpU0Xbxo0bx9GjR4mNjS3RY6anp2Nvb09aWhp2dnaliSuEScrOK+DDradYHnsZgFZu9kQMak09B5u/vV1Z11JFr3tZ80I87FxqBhNXxvN7SuHXKI8J8OT1rk0wN/v73r6k66lURwhyc3OJi4sjODi42Pbg4GBiYmIeeZvY2NiH9n/uuec4dOgQeXmP/uRkTk4O6enpxS5CiP+x1Gp4v3cLFg/zwd5Ky7EraXSPiGJD/NVyf6yKWPey5oX4Zw2dbNkwoQMj/esB8GX0RfovjOHizcxyuf9SFQQ3b96koKAAZ2fnYtudnZ1JSUl55G1SUlIeuX9+fj43b9585G1mzZqFvb190cXd3b00MYWoMp5r7sK2sEDaedYkM7eAtYfLf4hRRax7WfNClIylVsPMXs1ZMrwt1a21JFxN4/Dl8jkDoUwfKlSpVMV+VhTloW3/tP+jtj8wffp00tLSii5JSUlliSlEleBa3YpVL7Xnja5NmDPQ62/X4pPQ57qXNS9E6XRpVjiz4M2QJvT3cSuX+zQrzc4ODg5oNJqHuoLU1NSHuoEHXFxcHrm/mZkZtWrVeuRtLCwssLCwKE00Iao0jVrF+E4N9HLfFbHuZc0LUXq17a0YF1R+675URwjMzc3x8fEhMjKy2PbIyEj8/f0feRs/P7+H9v/ll19o27YtWq22lHGFEBVN1r0QVYRSSqtXr1a0Wq2ydOlS5eTJk8qUKVMUGxsb5dKlS4qiKMqbb76pDBs2rGj/CxcuKNbW1sqrr76qnDx5Ulm6dKmi1WqVH3/8scSPmZaWpgBKWlpaaeMKIf6krGupote9rHkhyk9J11Op3jIACA0N5datW7z//vskJyfTokULtm7dioeHBwDJycnFzk329PRk69atvPrqq8yfPx9XV1ciIiLo379/+VQ0Qgi9k3UvhOkr9RwCQ0hLS6N69eokJSXJOclCPIH09HTc3d25e/cu9vb2ho7zWLLmhSg/JV33pT5CYAgZGYVDGORUJCHKR0ZGRqUuCGTNC1H+/mndG8URAp1Ox7Vr17C1tX3saU4PKiBj6ygkd8UxxsxQvrkVRSEjIwNXV1fU6sr73WYlWfMg/08rkjFmBuPMXd6ZS7rujeIIgVqtxs2tZOdZ2tnZGc3/9D+T3BXHGDND+eWuzEcGHijNmgf5f1qRjDEzGGfu8sxcknVfeVsEIYQQQlQYKQiEEEIIYToFgYWFBTNmzDC6aWeSu+IYY2Yw3twVwVhfG2PMbYyZwThzGyqzUXyoUAghhBD6ZTJHCIQQQghRdlIQCCGEEEIKAiGEEEJIQSCEEEIITKggWLBgAZ6enlhaWuLj40NUVJShIxWZOXMmKpWq2MXFxaXoekVRmDlzJq6urlhZWdGpUydOnDhR4Tn37NlDz549cXV1RaVSsWHDhmLXlyRnTk4OkyZNwsHBARsbG3r16sWVK1cMlnnkyJEPvfbt27c3aOZZs2bx9NNPY2tri5OTE3369OH06dPF9qmMr3VlJOv+yRjjmi9Jbln3ZWMSBcGaNWuYMmUKb7/9NvHx8QQGBhISElLs29cMrXnz5iQnJxddEhISiq6bPXs2n376KfPmzePgwYO4uLjQpUuXonnuFSUzMxMvLy/mzZv3yOtLknPKlCmsX7+e1atXEx0dzb179+jRowcFBQUGyQzQtWvXYq/91q1bi11f0Zl3797NhAkT2LdvH5GRkeTn5xMcHExmZmbRPpXxta5sZN0/OWNc8yXJDbLuy0RvX8Bcgdq1a6eMGzeu2LYmTZoob775poESFTdjxgzFy8vrkdfpdDrFxcVF+eijj4q2ZWdnK/b29sqiRYsqKOHDAGX9+vVFP5ck5927dxWtVqusXr26aJ+rV68qarVa2b59e4VnVhRFGTFihNK7d+/H3sbQmRVFUVJTUxVA2b17t6IoxvFaVway7suXMa75R+VWFFn3ZWX0Rwhyc3OJi4sjODi42Pbg4GBiYmIMlOphZ8+exdXVFU9PTwYNGsSFCxcAuHjxIikpKcXyW1hYEBQUVKnylyRnXFwceXl5xfZxdXWlRYsWBn0uu3btwsnJiUaNGvHSSy+RmppadF1lyJyWlgZAzZo1AeN+rSuKrHv9M/bfQ1n3pWf0BcHNmzcpKCjA2dm52HZnZ2dSUlIMlKo4X19fli9fzs8//8ySJUtISUnB39+fW7duFWWszPmBEuVMSUnB3NycGjVqPHafihYSEsKKFSv47bffmDNnDgcPHuSZZ54hJycHMHxmRVGYOnUqAQEBtGjRoijTgwyPy2To3IYm617/jPn3UNZ92RjFtx2WxF+/IlVRlL/92tSKFBISUvTfLVu2xM/PjwYNGvDNN98UfdClMuf/s7LkNORzCQ0NLfrvFi1a0LZtWzw8PNiyZQv9+vV77O0qKvPEiRM5duwY0dHRD11nbK+1IVTmdWMq694Yfw9l3ZeN0R8hcHBwQKPRPFQdpaamPlRpVRY2Nja0bNmSs2fPFn3quLLnL0lOFxcXcnNzuXPnzmP3MbTatWvj4eHB2bNnAcNmnjRpEhs3bmTnzp3FvurXVF5rfZJ1r3+m9Hso675kjL4gMDc3x8fHh8jIyGLbIyMj8ff3N1Cqv5eTk8OpU6eoXbs2np6euLi4FMufm5vL7t27K1X+kuT08fFBq9UW2yc5OZnjx49Xmudy69YtkpKSqF27NmCYzIqiMHHiRNatW8dvv/2Gp6dnsetN5bXWJ1n3+mdKv4ey7kse0uitXr1a0Wq1ytKlS5WTJ08qU6ZMUWxsbJRLly4ZOpqiKIoybdo0ZdeuXcqFCxeUffv2KT169FBsbW2L8n300UeKvb29sm7dOiUhIUEZPHiwUrt2bSU9Pb1Cc2ZkZCjx8fFKfHy8AiiffvqpEh8fr1y+fLnEOceNG6e4ubkpO3bsUA4fPqw888wzipeXl5Kfn1/hmTMyMpRp06YpMTExysWLF5WdO3cqfn5+Sp06dQyaefz48Yq9vb2ya9cuJTk5uehy//79on0q42td2ci6f3LGuOb/Kbes+7IziYJAURRl/vz5ioeHh2Jubq60adOm6FSOyiA0NFSpXbu2otVqFVdXV6Vfv37KiRMniq7X6XTKjBkzFBcXF8XCwkLp2LGjkpCQUOE5d+7cqQAPXUaMGFHinFlZWcrEiROVmjVrKlZWVkqPHj2UxMREg2S+f/++EhwcrDg6OiparVapW7euMmLEiIfyVHTmR+UFlK+++qpon8r4WldGsu6fjDGu+X/KLeu+7OTrj4UQQghh/J8hEEIIIcSTk4JACCGEEFIQCCGEEEIKAiGEEEIgBYEQQgghkIJACCGEEEhBIIQQQgikIBBCCCEEUhAIIYQQAikIhBBCCIEUBEIIIYRACgIhhBBCAP8fPXwBQ98vSD0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 600x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head = flopy.utils.binaryfile.HeadFile(model_ws+'/'+model_name+\".hds\").get_alldata()\n",
    "head_official = flopy.utils.binaryfile.HeadFile(\"./official_examples/ex12_pht3d/mf2k/heads.dat\").get_alldata()\n",
    "\n",
    "plt.figure(figsize=(6,2))\n",
    "plt.subplot(121)\n",
    "plt.plot(head[-1, 0, 0])\n",
    "plt.title(\"head\")\n",
    "plt.subplot(122)\n",
    "plt.plot(head_official[-1, 0, 0])\n",
    "plt.title(\"official head\")\n",
    "\n",
    "print(\"R2 score: \", r2_score(np.ravel(head_official), np.ravel(head)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'./models_folder/case12\\\\pht3d_ph.dat'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Create PHT3D Reactive Transport Model, supported by PHT3D-FSP\n",
    "## Define \"spec\", calling PHT3D-FSP (adopt function variables as needed)\n",
    "spec = pht3d_fsp.create(\n",
    "    xlsx_path=\"./\", \n",
    "    xlsx_name=f\"./data/{data_folder}/pht3d_species.xlsx\",\n",
    "    nlay=nlay, \n",
    "    nrow=nrow, \n",
    "    ncol=ncol,\n",
    "    ph_os=2,\n",
    "    ph_temp=25,\n",
    "    ph_asbin=1,\n",
    "    ph_eps_aqu=0,\n",
    "    ph_ph=0,\n",
    "    # ph_print=0,\n",
    "    ph_cb_offset=0,\n",
    "    ph_surf_calc_type=\"-no_edl\",\n",
    "    write_ph=\"yes\"\n",
    ")\n",
    "# Move the pht3d_ph.dat file into the models_folder\n",
    "source = 'pht3d_ph.dat'\n",
    "destination = os.path.join(model_ws, 'pht3d_ph.dat')\n",
    "shutil.move(source, destination)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: [0,\n",
       "  0,\n",
       "  0,\n",
       "  0.0,\n",
       "  2,\n",
       "  1e-06,\n",
       "  5e-05,\n",
       "  0.0098456,\n",
       "  0.01,\n",
       "  0.0,\n",
       "  4.26,\n",
       "  4,\n",
       "  0.0,\n",
       "  0.0,\n",
       "  0.0],\n",
       " 1: [0, 0, 0, 0.0, 2, 0.0, 0.0, 0.0099398, 0.01, 0.0, 4.26, 4, 0.0, 0.0, 0.0]}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define SSM and RCH data (for the SSM package) for each stress period\n",
    "ssm_data = {}\n",
    "itype = flopy.mt3d.Mt3dSsm.itype_dict() # ... to check the key words\n",
    "# print(itype)\n",
    "# {'CHD': 1, 'BAS6': 1, 'PBC': 1, 'WEL': 2, 'DRN': 3, 'RIV': 4, 'GHB': 5, 'MAS': 15, 'CC': -1}\n",
    "\n",
    "wel_spec={}\n",
    "\n",
    "for i in range(1):\n",
    "    for key in spec.keys():\n",
    "        wel_spec[key]=spec[key][0,0,0]\n",
    "\n",
    "    wel_spec['tracer'] = 0.000001\n",
    "    wel_spec['u6'] = 0.00005\n",
    "    wel_spec['na'] = 0.0098456\n",
    "    wel_spec['n5'] = 0.01\n",
    "    wel_spec['f'] = 0.0\n",
    "\n",
    "    wel_spec['ph'] = 4.26\n",
    "    wel_spec['pe'] = 4\n",
    "\n",
    "    wel_spec['saoh'] = 0.0\n",
    "    wel_spec['sboh'] = 0.0\n",
    "    wel_spec['scoh'] = 0.0\n",
    "\n",
    "    ssm_per=[0,0,0,0.0,itype['WEL']]\n",
    "    for key in spec.keys():\n",
    "        ssm_per.append(wel_spec[key])\n",
    "    \n",
    "    ssm_data[i] = ssm_per\n",
    "\n",
    "for i in range(1,2):\n",
    "    for key in spec.keys():\n",
    "        wel_spec[key]=spec[key][0,0,0]\n",
    "\n",
    "    wel_spec['tracer'] = 0.0\n",
    "    wel_spec['u6'] = 0.0\n",
    "    wel_spec['na'] = 0.0099398\n",
    "    wel_spec['n5'] = 0.01\n",
    "    wel_spec['f'] = 0.0\n",
    "\n",
    "    wel_spec['ph'] = 4.26\n",
    "    wel_spec['pe'] = 4\n",
    "\n",
    "    wel_spec['saoh'] = 0.0\n",
    "    wel_spec['sboh'] = 0.0\n",
    "    wel_spec['scoh'] = 0.0\n",
    "\n",
    "    ssm_per=[0,0,0,0.0,itype['WEL']]\n",
    "    for key in spec.keys():\n",
    "        ssm_per.append(wel_spec[key])\n",
    "    \n",
    "    ssm_data[i] = ssm_per\n",
    "\n",
    "ssm_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DSP: setting dmcoef for component 2 to zero, kwarg name dmcoef2\n",
      "DSP: setting dmcoef for component 3 to zero, kwarg name dmcoef3\n",
      "DSP: setting dmcoef for component 4 to zero, kwarg name dmcoef4\n",
      "DSP: setting dmcoef for component 5 to zero, kwarg name dmcoef5\n"
     ]
    }
   ],
   "source": [
    "model_name = \"case12_mt3dms\"\n",
    "icbund = np.abs(ibound)\n",
    "\n",
    "## Initiate model object\n",
    "mt = flopy.mt3d.Mt3dms(model_name, model_ws=model_ws, exe_name='./bin/mt3dms',\n",
    "                      ftlfilename='mt3d_link.ftl',modflowmodel=model, namefile_ext='nam_pht3d')\n",
    "\n",
    "prsity = 0.42\n",
    "save_times = np.loadtxt(\"./data/case12/save_times.txt\")\n",
    "\n",
    "exec(f'btn = flopy.mt3d.Mt3dBtn(mt, \\\n",
    "                                nper=nper, perlen=perlen, nstp=nstp, dz=0.00038, \\\n",
    "                                nlay=nlay, ncol=ncol, nrow=nrow, \\\n",
    "                                laycon=[1], prsity=prsity, icbund=icbund, nprs=1, mxstrn=50000, tsmult=1, timprs=save_times, \\\n",
    "    ncomp={pht3d_fsp.create.ncomp}, mcomp={pht3d_fsp.create.mcomp}, {pht3d_fsp.create.sconc_btn})') # call BTN package in this way to invoke PHT3D-FSP variables\n",
    "\n",
    "adv = flopy.mt3d.Mt3dAdv(mt, mixelm=-1, percel=0.75, mxpart=5000, nadvfd=0)\n",
    "dsp = flopy.mt3d.Mt3dDsp(mt, al=0.000075, trpt=0.1, trpv=0.1, multiDiff=True,dmcoef=0.0,)\n",
    "gcg = flopy.mt3d.Mt3dGcg(mt, isolve=1, cclose=0.00001)\n",
    "exec(f'ssm = flopy.mt3d.Mt3dSsm(mt, stress_period_data=ssm_data, mxss=1000)')\n",
    "mt.write_input()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Manually add the PHC reaction package to the nam file and push out as pht3d.nam\n",
    "s='PHC               39  pht3d_ph.dat\\n'\n",
    "namfiletxt=open(model_ws + '/'+ mt.namefile, 'r').read()\n",
    "pht3d_nam = namfiletxt+s\n",
    "file = open(model_ws + '/pht3d.nam','w')\n",
    "file.write(pht3d_nam)\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "copy_pht3d_datab(source_folder=data_folder)\n",
    "\n",
    "run_pht3d_program(work_dir=data_folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracer = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D001.UCN').get_alldata()\n",
    "u6 = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D002.UCN').get_alldata()\n",
    "na = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D003.UCN').get_alldata()\n",
    "n5 = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D004.UCN').get_alldata()\n",
    "f = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D005.UCN').get_alldata()\n",
    "ph = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D006.UCN').get_alldata()\n",
    "saoh = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D008.UCN').get_alldata()\n",
    "sboh = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D009.UCN').get_alldata()\n",
    "scoh = flopy.utils.binaryfile.UcnFile(model_ws+'/PHT3D010.UCN').get_alldata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAADUCAYAAADQgpstAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAv20lEQVR4nO3deXxU5bnA8d+ZyWSfBLInJIGwhC0xQthFxQUUioq4g4ptrXXBW4pd9N5exV7r1qu1rXWtWlu0aC1arIiALC5cEAJI2EPYs5KQZLIxmcy894+TxUACScjMmcw834/nM5OZszznZXzmzHveRVNKKYQQQridyegAhBDCX0jCFUIID5GEK4QQHiIJVwghPEQSrhBCeIgkXCGE8BBJuEII4SGScIUQwkMk4QohhIdIwhVCCA8xLOF+8cUXXHPNNSQlJaFpGh999JFbj7do0SI0TWuzJCQkuPWYQgjxXYYl3NraWrKysnjxxRc9dsyRI0dSVFTUsuTm5nrs2EIIEWDUgadPn8706dM7fL+hoYFf/epXvPPOO1RWVpKRkcEzzzzDlClTun3MgIAAuaoVQhjGa+twv//97/P111+zZMkSduzYwU033cTVV19NXl5et/eZl5dHUlISaWlp3HrrrRw8eLAHIxZCiLPTvGF4Rk3T+PDDD5k1axYA+fn5DBkyhOPHj5OUlNSy3pVXXsm4ceN48sknu3yMTz/9lLq6OtLT0ykpKeGJJ55g79697Nq1i+jo6J46FSGE6JBXXuFu3boVpRTp6emEh4e3LOvXryc/Px+Aw4cPn3ET7PRl/vz5LfucPn06N9xwA5mZmVx55ZV88sknALz99tuGnKMQwv8YVod7Ni6XC7PZTE5ODmazuc174eHhAPTr1489e/acdT99+/bt8L2wsDAyMzPPq4pCCCG6wisT7qhRo3A6nZSWlnLxxRe3u47FYmHYsGHdPobdbmfPnj0d7l8IIXqaYQm3pqaGAwcOtPx96NAhtm/fTlRUFOnp6cydO5c777yT5557jlGjRlFWVsaaNWvIzMxkxowZXT7ez372M6655hpSU1MpLS3liSeewGazMW/evJ48LSGE6JgyyNq1axVwxjJv3jyllFINDQ3q0UcfVQMGDFAWi0UlJCSo66+/Xu3YsaNbx7vllltUYmKislgsKikpSc2ePVvt2rWrB89ICCHOzitaKQghhD/wylYKQgjhizxeh+tyuSgsLMRqtaJpmqcPL4QQPU4pRXV1NUlJSZhMHV/HejzhFhYWkpKS4unDCiGE2x07dozk5OQO3/d4wrVarYAeWERERKe3czgcrFy5kmnTpmGxWNwVnteTcmglZdFKyqKVEWVhs9lISUlpyW8d8XjCba5GiIiI6HLCDQ0NJSIiwq8/UFIOraQsWklZtDKyLM5VTSo3zYQQwkMk4QohhIdIwhVCCA/xyrEUTvfF/hP8fvV+KirMLC7ajNmkoaFhMoGpqc7EpGmYNL0OpflRA0IDzcRag4i1BjEoNpzM5EjirMHGnpAQwi/1ioRbXmsn52gloHGwuuK89zciMYLrLkzixuxkosODznt/QgjRGb0i4Y5Pi+ZPt2WxJWcro0aNwmQ241J6Y2OlwHXao0LhUuB0KeoaGjlRbafYZmd/cTX7SqrZXWRjd5GNF1bnceek/sy/bDDWYP++syuEcL9ekXCT+oQQGxZP42HF9IyE82rqUV5j57NdJbz7zRF2Fth4df1BPtxawHM3Z3HxkNgejFoIIdryu5tm0eFBzBmfysfzJ/PnO8fQPzqU0mo7d775DS+s3o/TJWP5CCHcw+8SbjNN07hyRDyfLbiE28aloBS8sDqP7/9lMzX2RqPDE0L4IL9NuM2CLWaemn0Bz92URbDFxBf7T3D7nzdRWddgdGhCCB/j9wm32Q3Zybx3z0T6hFrYfqySW1/bSGn1KaPDEkL4EEm435GV0of37plIrDWIvcXVzHl9Eydr5UpXCNEzJOGeZmiClX/8eCIJEcEcKK1h3pvfUH3KYXRYQggfIAm3HQNiwlh89ziiwgLJLaji7re3cMrhNDosIUQvJwm3A4PjrPz1B+OwBgWw6dBJfvredlzSZEwIcR4k4Z5FRr9IXp83hkCziU93FvPMir1GhySE6MUk4Z7DhIHRPHvjBQC8+sVB3tl0xOCIhBC9lSTcTpg1qh8Lp6YD8Oi/drFuX6nBEQkheiNJuJ304OWDuTE7GadL8cA7W9lTZDM6JCFELyMJt5M0TePJ6zOZNCia2gYnd7+9hRPVdqPDEkL0IpJwuyAwwMTLc7MZGBNGQWU99/xNmosJITpPEm4XRYZaeOOusUSGWNh2tJJffLADpaS5mBDi3CThdkNaTBgv3z6aAJPGsm8LeXHNAaNDEkL0ApJwu2nSoBj+Z1YGAM+t2s8nO4oMjkgI4e0k4Z6H28al8sPJaQA89I/t7DheaWxAQgivJgn3PP3njOFcNjSWUw4Xd7+9haKqeqNDEkJ4KUm458ls0vjDbaMYGm+ltNrO3W9voa5BZowQQpxJEm4PsAZb+PO8MUSHBbKr0MbC976VgW6EEGeQhNtDUqJCefWObALNJlbsKuZ/V+4zOiQhhJc5r4T71FNPoWkaCxYs6KFwercxA6J45sZMAF5al88/c44bHJEQwpt0O+Fu3ryZ1157jQsuuKAn4+n1rh+VzAOXDQLgkaW5bD580uCIhBDeIqA7G9XU1DB37lxef/11nnjiibOua7fbsdtbxxyw2fRBXxwOBw5H56euaV63K9sY5T+mDORASTWf7S7lnr9u4Z/3jielb2iP7Ls3lYO7SVm0krJoZURZdPZYmupGv9R58+YRFRXF7373O6ZMmcKFF17ICy+80O66ixYt4vHHHz/j9XfffZfQ0J5JQt7I7oQ/7DJzvFYjIUTx0wwnwd36ehNCeLu6ujrmzJlDVVUVERERHa7X5RSwZMkStm7dyubNmzu1/iOPPMLChQtb/rbZbKSkpDBt2rSzBnY6h8PBqlWrmDp1KhaLpathG2Lipae44ZVNFFfbWV6ZwCtzLyTAfH73KXtjObiLlEUrKYtWRpRF8y/3c+lSwj127Bg/+clPWLlyJcHBwZ3aJigoiKCgoDNet1gs3SqM7m5nhJRovbnYza/+H+vzyvjtqnwevWZEj+y7N5WDu0lZtJKyaOXJsujscbp0uZWTk0NpaSnZ2dkEBAQQEBDA+vXr+cMf/kBAQABOpwxVeLoLkvvw/M0XAvDm14d4d9NRYwMSQhimSwn3iiuuIDc3l+3bt7csY8aMYe7cuWzfvh2z2eyuOHu1GZmJPNQyRc9ONhwoMzgiIYQRulSlYLVaycjIaPNaWFgY0dHRZ7wu2pp/+WDyT9Tw0fZC7l2cw0cPXMTA2HCjwxJCeJD0NPMQTdN4+oYLGJ3aB9upRn749hYq6xqMDksI4UHnnXDXrVvXYZMw0Vawxcyrd4yhX58QDpXVcv87W3E4XUaHJYTwELnC9bBYaxB/njeGsEAzG/LLefRfu2SKHiH8hCRcAwxPjOD3t45C0+Dv3xzlza8PGx2SEMIDJOEa5MoR8fzn9OEA/OaT3azdW2pwREIId5OEa6C7L07jljEpuBQ8+Pdt7CuuNjokIYQbScI1kKZp/M+sDManRVFjb+QHf9lMWY393BsKIXolSbgGCwww8crt2QyIDqWgsp4f/y0He6P02BPCF0nC9QJ9wwL587yxWIMDyDlSwSP/zJWWC0L4IEm4XmJwXDgvz83GbNJYuq2Al9blGx2SEKKHScL1IpOHxLDo2pEA/PazfazYWWRwREKIniQJ18vcMaE/8yb2B+Cn733LzoIqgyMSQvQUSbhe6L9njuCS9FjqHU5++PZmSmynjA5JCNEDJOF6oQCziRfnjGJwXDglNjt3v72F+gZpuSBEbycJ10tFBFt4Y94Y+oZayC2o4qF/bMflkpYLQvRmknC9WP/oMF69YwwWs8by3GJe+DzP6JCEEOdBEq6XG5cWxW9mZQLwh8/z+HiHtFwQoreShNsL3Dw2hXsuGQjAwx/u4rAMuSBEryQJt5f45dXDuHJ4HA2NLt7YZ6aoSlouCNHbSMLtJcwmjRduHcXQ+HBsDo0fL95GXUOj0WEJIbpAEm4vEh4UwKu3jyLcothTXM1P35OWC0L0JpJwe5l+fUK4e6gTi1njs10lPLdqn9EhCSE6SRJuL5RmhSdn6WMu/GltPh9uO25wREKIzpCE20vNujCJ+6cMAuCXH+SSc6TC4IiEEOciCbcX+9m0oVw1Mp4Gp4v7FufImAtCeDlJuL2YyaTx/M0XMjTeSmm1nfsWy2wRQngzSbi9XFhQAK/ekU1EcABbj1ayaNluo0MSQnRAEq4PGBATxh9uG4Wmwd+/Oco7m44YHZIQoh2ScH3ElKFx/GzaUAAWLdvFlsMnDY5ICHE6Sbg+5P4pg5iRmYDDqbjvna1yE00ILyMJ14domsZvb8xiaLyVE9V27l2cQ0Ojy+iwhBBNJOH6mLCgAF67MxtrcADbjlbyzIq9RockhGgiCdcH9Y8O47mbsgB446tDMvuvEF5CEq6PmjYyoWUM3Z//YweHy2oNjkgIIQnXh/38qqGMHdCXansj97+zlVMO6RQhhJEk4fowi9nEH28bTXRYILuLbDz+sXSKEMJIknB9XEJkMC/cemFLp4hPc6U+VwijSML1AxcPieXeS/WRxR5emkthZb3BEQnhnyTh+omFU9PJSo6kqt7Bwve345SZIoTwOEm4fsJiNvH7W0cRGmhm48GTvLI+3+iQhPA7knD9yICYMB6/Vp8p4ner9rP9WKWxAQnhZyTh+pkbs5P53gWJNLoUC5Zso75BmooJ4SmScP2Mpmk8eX0mCRHBHC6v47mVMgmlEJ4iCdcPRYZYeGp2JgBvfH2InCMylKMQniAJ109dNiyOG0YnoxT8/IMd0gtNCA+QhOvHHp05gjhrEAdP1PK7VfuNDkcInycJ149Fhlp48nq9auH1Lw+y7ahMtS6EO3Up4T711FOMHTsWq9VKXFwcs2bNYt8+uenSm105Ip7rR/XDpeAXH+yQAcuFcKMuJdz169fzwAMPsHHjRlatWkVjYyPTpk2jtlaG/uvNHrtmBNFhgeSV1vDm14eMDkcInxXQlZVXrFjR5u+33nqLuLg4cnJyuOSSS9rdxm63Y7fbW/622WwAOBwOHA5Hp4/dvG5XtvFF7iiHMIvGL69K5xdLd/L71fu5ekQs/fqE9Nj+3cXrPxPleZg2v47p8BeAhopJx5VxM2rodNB6tjbP68vCg4woi84eS1NKdbtT/YEDBxgyZAi5ublkZGS0u86iRYt4/PHHz3j93XffJTQ0tLuHFj1MKfjjLjP51RqZfV3cPUyqFs5HUsVGRh39MwGuhjPeqwgdyPaU72ML7W9AZMId6urqmDNnDlVVVURERHS4XrcTrlKK6667joqKCr788ssO12vvCjclJYWysrKzBnY6h8PBqlWrmDp1KhaLpTsh+wR3lsP+kmque2kjjS7Fq7eP4vKhsT26/57mrZ8J7eBazH+/GQ2Fq/9kXOPuhcBQtINrMeW8hdZQgzIH4Zz+W1TWnB45preWhRGMKAubzUZMTMw5E26XqhS+a/78+ezYsYOvvvrqrOsFBQURFBR0xusWi6VbhdHd7XyNO8phZHIUP5ycxqtfHOTpFfu5bFgCgQHe35DFqz4TVQXw0T2AggvnYrr2j5hMZv29IVfApAfh4/9A27+CgH//BxRtg+nPQkBgjxzeq8rCYJ4si84ep1v/Nz344IMsW7aMtWvXkpyc3J1dCC81//LBxIQHcqislsUbjxgdTu+z6lGor4DEC+F7z0Nzsm1mjYdb/w6X/QrQIOct+Pst0CA3nv1BlxKuUor58+ezdOlS1qxZQ1pamrviEgaxBlt4aNpQAH7/eR6VdWfWQYoOFOTAzg8ADa79I1iC21/PZIJLfw5z/wGWMMhfA3+9Duqki7Wv61LCfeCBB1i8eDHvvvsuVquV4uJiiouLqa+XGQR8yc1jUhiWYKWq3sELq/OMDqf3WPOE/ph1KyRecO71h0yFO/8FwX3g+Gb4y/egtsytIQpjdSnhvvzyy1RVVTFlyhQSExNblvfee89d8QkDmE0av/reCAAWbzzCIZli/dxK9+hXqpoJpjzc+e1SxsIPVkB4ApTulitdH9flKoX2lrvuustN4QmjTB4Sw2VDY2l0KRlnoTM2vqw/DpsJfQd0bdu44XDXvyE8Hkp26km3XrpZ+yLvvwUtDNNcl/vxjkL2FtsMjsaL1VfCjqZfeRPu794+YobAncsgNAaKd8DfZsOpqh4LUXgHSbiiQxn9IpmRmYBS8PxKucrt0J6PofEUxA6H1And30/cMJi3DEKioHArLL5RWi/4GEm44qwWTk3HpMHK3SV8K3OgtS/3ff3xgptA085vX/Ej4c6Pmm6kfQP/+D44G883QuElJOGKsxocZ2XWqH4APCd1uWeyFcGhpp6WGTf0zD4Ts2DOexAQDHmfwb9/ove9Fr1et3ua9WpKwdGNcORrqDoOjjp9UUpvqK6ZQGt6NJn1xRIGgaEQGAaBVgiP1W9yhCfojdmDrEafldssuCKdf20v5Iv9J8g9XkVmcqTRIXmPXUsBBSnju36z7GxSJ8CNb8J7t8O2xRDRDy77z57bvzCE/yXcmlJY9iDsX3HudbsiPB5i0iF6MMQO1XsaJWbpSbqXS40O5dqsJD7cVsBL6w7w8u3ZRofkPXY0VSdk3tTz+x72Pb232r8XwPpn9M9X5o09fxzhMf6VcJWCD34Ah78EkwWGz4SYofpVqyVEr39TCpQLXE79UTnB1QgNTVfBDTVwyga1J6CmRE/gdlvT8xJ93800M8SNgORs6H8RDLgYIhKNO//zcN+UQXy4rYAVu4o5UFrD4Lhwo0MyXlkeFG3X/51HXu+eY4z5PpzMhw1/hH89AH3T9M+T6JX8K+Hu+1RPiOYg+NHnkJDZM/s9ZYPyPCg7AGX79UbwBTlQUwwlufqS8xd93eghkHaxnnwHToHQqJ6Jwc3S461cOTye1XtKeHV9Pr+9KcvokIyX+4H+OOhyCItx33GufFxP7vtXwJI5cM9aiEhy3/GE2/hPwlUKVi/Sn098oOeSLUBwBPTL1pfvshXqiffYJv3GStG3emIuz4MtbwKavs2QqTD4SkgadeZgJ17k/ssGsXpPCR9uK+CnU9NJ6gWDlLuNUpD7D/25O6oTvstkhtmvwxvT4MQeveXCXf8Gs4wK1tv4T8KtOARl+8AcCJN/6pljRiTpy/Br9L/rK+DIBj35HvoCSndBwRZ9WfeU3v5y8BV68h10hX5jzouMTu3LhIFRbDx4kje+OsR/zxxhdEjGKdyq/9QPCIFhM9x/vOAIuPUdeG0KHNuoXzxc9Rv3H1f0KP9JuEf+T39MGq1/eI0Q0le/ETLse/rftkI4sFpf8tdB/Un9qqn5yilplJ58B0/Vr4TNxv9z/fjSQWw8eJL3Nx/jp1PTCQ8yPiZDNFcnDJvhuRYq0YNg1kt6y4X/e1FvGTHiWs8cW/QI//m/5WhTwj2fnkA9LSIJRt+pL04HHN8CB1ZB3iq9e2fhNn354rd6Q/iBU9AGXEpIg3FtMi8dEsvAmDAOltXyz5zjzJs0wLBYDONyws5/6s/dXZ1wuuHX6IOYN99ES8yCvjJVT2/hRwl3o/6YOtHYODpitkD/ifpyxaNQXQL5nzdd/a7RqyN2f0TA7o+YBqjil/Xqh0GX6y0ggjzTasBk0pg3aQCPLdvFXzYc5o4J/TGZzrN3VW9z+Eu9RUpwH73qp4tK60rZULiBgpoCGl2NWAOt9AnqQ0xIDElhSSSFJxFqOUtzwiseg2Pf6PcGPvwx3PVJ989FeJR/JNzaMv1GFUDKOGNj6SxrPFw4R19cTv3mW/5aXAc+h+ObMTXffNv0it7ELXUCDLpMT8AJWfog125yQ3Yy//vZPg6V1bI+7wSXDY1z27G80o6mKp+Rs7o0NY5LuXh9x+u8/O3LOJXzrOv2DepLUriefPuF92t9bE7Is1+Dlyfrv9y+eh4mLuj++QiP8Y+Ee3yz/hg7vNc0w2rDZNa/KFLG4bxoISs//oCrhoQQcGQ9HPgcKo/oV12Hv4TPfw2h0TBgst70bMBkiB12/n38vyM8KICbx6bwxleHeOvrw/6VcB2nYM8y/XnmzV3a9MlNT/LePn1UsZHRIxkZPRKL2YLNbqPSXsmJ+hMU1BRQ3VBNhb2CCnsFu8p3tbuvyKBIkgYPJ6nsIInfvkSiuYGCBhMDTw4ktU8qEYERaD34by56hn8k3JMH9ce4YcbG0UMazaGoYTMgc5b+wsmDerVD/lo4uB7qymH3v/QF9CH/BkxuTcKxQ887Ac+bOIA3vz7EF/tP+FdHiLzP9I4uEf26VD21s2wn7+/Te6X9etKvuX5Ixx0lbA02imqKKKgpoKhWfyysKaSwppDjNcepbqimyl5FFVXsCWuqejiwBIB3V7wLQJgljMSwRPqF9yMxLLHlajkpLInE8ESig6MlIRvAPxJu5VH9sY+P3lyIGqgvY+/Wb74VbG264v1Kr7uuK4PdH+kLQFisnnxTJ0HqeIgb2eUWEKnRoVwxLI7Ve0p5d9NRHr3GT5qIbVusP2be1OlqG6UUT3/zNArFzIEzz5psASICI4iIimBo1NB2369pqKGwtlBPyhV5FH39HIUuOwdDYjgZFsLJUyepddRyoPIAByoPtLuPIHNQSyJuTswJYQktS3xoPIHmnplJWLTyj4Rb0TT7bJ9UY+PwBLNFT6Kp4+GSn0Fjg95mtCUBb9K7Je/6UF8AAsP1ZmcpTdslj4Xgcw9QM3dCf1bvKeWDnGP8/KqhhAR6b6eNHtHcjA9g1B2d3mx/xX6+PfEtgaZAfpp9/m3AwwPDSQ9MJ71vOqRcCgFx8P4duKjA+cPVNCaOoKi2iKKaotbE3HS1XFhTSGldKXanncO2wxy2He7wONHB0W2ScEJoQpu/Y0JiCDD5RwrpKf5RWpVNCdcfm88EBOo31FInwCU/h0Z70xXwV3oD+mObwV4Fh9brCwCaPgZEyjh9u5Rxeh/+036CXjIkluS+IRyvqOffOwq5aUyK58/Pk7a/q4+vkToJYgZ3erO1x9YCMClpEnGhbqjvHnEtrmHXYtq7DO3fP8Fyz1oGRg5kYOTAdld3OB0U1xWfkYiLa4spriumuLYYu9NO+alyyk+Vd1iPbNJMxITEtJuMm/+ODonGpMkosM18P+Eq5ftVCl0RENTa/AzA5YITe/Xke3ST3tSo4pDeC650F+S8pa8XEqV3xOg3Wu880m80ZmsCc8an8uyKfSzedNS3E67TAVuaymJ0569uAdYcXQPAZamX9XRULZxXPU1j3ucEluTC//0JJi/ocF2L2UKKNYUUa/v/XkopKu2VegL+ThJuXkrqSiipK6HR1UhpXSmldaXsYEe7+wowBRAfGk98aHxLMo4LjWtZ4kPj/epK2ffPsu6kPsIXQKQPJ4TuMpkgfoS+jPmB/lp1iT7bwNGNenvPou16L7j8z/WlmTWJH8ZnUW+JYFvBQHYfTGbEQB/9Utv9L7Ad1+u/R87u9GbFtcXsObkHDY1Lki9xX3zhcezsN4fRR1+H9c/qdcyR/bq1K03T6Bvcl77BfRkePbzddVzKRXl9+RkJuaSupOX5ifoTNLoaKagpoKCmoOPjoREdEt0mCbck5ZCmx7A4rBZrr7/R5/sJt/Kw/mhNBEuwoaH0GtZ4vUdT8xgQjXZ9NtmCrXrPt4Kt+lVxdSFB1YU8ZAbMwF+f0qse+o3We0AlZELCBe4dScsTlIINf9Cfj7unS5+jrwu+BiArNouYEPeWw7Goi7jQtQPT8U3w2X/CzW+77VgmzURsaCyxobFk0v5AUI2uRsrqy9pcHRfXFbdcFZfWlXKi7gSNSl+vrL6M3eW7OzxmSEBIm6vjNgm5+Wo51Ls/a36QcJurE/zghpm7BASdORqavUYf/axwK2X7N1J78Bv6m0r16oiKQ61dX0H/skvI/M5ygZ6Y3dg5o0ft/UQ/14AQGPPDLm26s3wnAKPjR7sjsrY0E86rn8H0xuV6i5QDn+u9EQ0SYApoqUboiEu5OHnqZJskXFJX0pKMm5/bGmzUN9ZzxHaEI7YjHe5PQyMqOIogRxAr168kISyB2NBYvVojLJ6E0ATiw+IJs4S545TPyfcTbksLBR/9qWuUoHAYcBEMuIjoifO59XdfUFZaxP9e5OTKyEIo3gnFufqIWtVF+pK3snX7wHCIz2ibiOOG6wPBexOnA1Y9qj+f+ACERXdp8z3lewA6/Gne4+IzYNyPYdPL8Okv4L4N+heml2q+8RYTEsOI6I6bFtY31rdJwO0u9aU0uhopP1UOQGFBYYf7C7eEtyThNo+h8aRFppFsTe7xcwV/SLjNV7j+2ELBQzRNY+74VB7/uIbfHrByxYI5rXVt9moo2a0PxlOcqy+lu/V69WMb9aVlRyb9yjduuN5KovkxepBxY79u+KP+pREaAxf9pEubOpwO9lfoE2+OjBrpjujad9kj+lxr5Qf0UcUufshzx3aTkIAQUiNSSY3o+JeqS7moOFVBoa2Q5V8uJ3VEKicbTlJSW9Jy5VxSW0K1o5oaRw01VTXkV+WfsZ9bht7Cryb8yi3n4fsJt+qY/ig3zNxq9uhknlmxl30l1eQcqWDMgKYu1EHW1nbBzZyNejIozv1OIt6h95A7ma8ve//dur7JAjFDmhLwd5JxePduCnVaQQ6sbRpzduqvuzysZ35VPg6XA6vF6rYrpnYFR8K0J2Dpj2D9b/UbaH5QpWbSTESHRBMREMFBy0FmDJmBxXLmF3Wto7Yl+bZ5bLp6TotMc1uMfpBwj+uPkR78wPuhyBAL12Yl8f6W4yzeeKQ14bbHHKB3s44bBhc0DW+olN4ho3S3PkVRy+Me/Wq4dLe+fEeAJZRLLPGYnZ9Cwkh9zIiYIfqX6/nOnHFiP7x7qz6f3cjr9UGEuqj5BtDw6OGev7ueeRPkvA1HvoLP/gtu+Ztnj+/FwixhZ22n7E5+kHCbmqNIwnW72yf05/0tx1meW8yj1zQQFdaFrqGaBuFx+jJwSuvrSum/Ukr3tk3GJ/ahOero6zgEOw7RphloQLA+e3LMEH2m25h0/Xn0YH3C0LNRSr+6/niB3iU6PgNmvtCtsSeaE+7Z6ibdRtNgxrPwysX6YDsH17UtV2EI3064p2x6LyrQBxsRbnVBch8y+0WSW1DFP7Yc48eXDjr/nWqa/nO4TyqkT2t93eXEUbqfbZ+9Q3ZKKOayffoEnuUHoPGU3oytZOeZ+4tMaU3Ekcn6LBxBEfo2ZXmwb3nrdgkXwB0fQUifboXefMPMkIQLED9SH1/jm1dh+S/gvq9lHjSD+XbCtTVd3Qb38dgA3f7u9gmp/PKfubz7zVF+dPFA9w1ObjJD9GCK+ozFdfEMzM11dS6n3pW7LE9PwGX7W5/XletXy1XH9NHVOhIQDBPuhykPd/sOf6OrkX0V+wADEy7oN9B2fqDP5/fN6zDxfuNiET6ecKU6weOuyUriiU/2cKS8jq8OlHFJuocnwjSZW0dPS7+q7Xu15U3T2e+HE/ugplSfScNu05NsZAokj9EHFg/pe15h5FfmY3faCbOEddiF1iNC+uozRHz8H/pEpZk36tU2whA+nnCbWyhIwvWU0MAAbhidzF82HOadTUc8n3DPJixaXzwwr92ek03tb6OGGz94y6jbYcubehft1Y/DrD8ZG48f6yVdfbqpuUpB6m89as54vQnS6j2lFFedMjgaYxh6w+x0JjPM+F/9+fbF+mSlwhC+nXBbqhQk4XpSeryVcWlROF2KJZuPGh2OITzew+xcUsZCVlPTtuU/10eJEx7n4wlXOj0YZW7TVe6Sb47R6PSv/7mdLqd33DA73ZWLINCqD0i/fbHR0fgl3064UqVgmKszEogOC6TYdopPdxYbHY5H7T25l/rGeqwWK/2tXtSl3BoPU36pP1/53/ownMKjfDfhulxSpWCgoAAzd0zUk81L6/JRShkckedsLtZnic6Oz8Z8vj3eetr4e/X2xacqYfnPjI7G7/huwi0/AE67PqSeVCkY4q5JAwgLNLOnyMaavaVGh+Mxm0v0hDsmYYzBkbTDbIHr/gSmAL0H2q6PjI7Ir/huwi36Vn9MyDz/fvWiW/qEBnL7BP0q9/lV+3G5fP8qt9HVyNaSrQCMTRhrcDQdSLwAJi/Un3/ykN4+WXiEDyfc7fpj0oVGRuH3fnzpIKxBAewqtPHxjo7HJ/UVe0/upcZRg9ViZWjf9qc59wqX/Axih+vjRSybr48hIdzOhxNu0xVuYpaxcfi5qLBA7p2ij6nw7Ip91DU0GhyRe60+ok+jnp3ghfW33xUQBLNfBXOgPn7ExpeMjsgv+GbCdbm+k3AvNDQUAT+4KI1+fUIoqKzn+ZX7jQ7HbRqcDSzNWwrArEGzjA2mMxKz4Kon9ecr//vs40uIHuGbCbfikN4/3hwEsV78s85PhASaeWJWBgBvfn2IzYdPGhyRe6w6sooKewVxoXFcmnKp0eF0zti79Q4Rygnv39V6oSLcwjcTbvM3dUKGDEfnJS4bFsfsUf1wKXjgna2UVvtWl986Rx2vfPsKADem30iAqZcMU6JpcM0LkDpJH8r0b9fr89EJt/C9hOts1OehAn3Ue+E1/mdWBkPiwimttjPvzc2crG0wOqQe4VIufrPpNxy2HSYuNI45w7o+O4ShAoJgzhJIGqUPYfnWdDiw2uiofFK3Eu5LL71EWloawcHBZGdn8+WXX/Z0XN2X+74+HmpIFIy+0+hoxHeEBQXw2p1jiLUGsafIxk2vbGBvsc3osM5LYU0hP1nzE5blL0ND4+mLnyYyKNLosLouOFIfbL3/ZL06bvEN8OnDUF9pdGQ+pcu/e9577z0WLFjASy+9xEUXXcSrr77K9OnT2b17N6mpBk1UpxTUlsGO92D1Iv21CfefezoV4XFpMWH8/UcTmPvnjeSfqOXaP37NDdnJ3DC6Hxn9Igm2eNedfZdy4XA5cDgdVDdUU36qnBN1J8ivymdz8WY2FW3CqZwEmgJ5bNJj3tv2tjNC+sAdS2HFI7DlDX2q9W//rs/nNnSGfpMtyNqt6YaETlNd7HM5fvx4Ro8ezcsvv9zy2vDhw5k1axZPPfXUGevb7XbsdnvL3zabjZSUFMrKyoiI6NwsqBu3vsJru9/C2diI2WwGFBqgUHqLBFdDSztCBRAShYpK06fdbtLRaSo6f/pd3Ud7r3e1i+vp+1BKYbPZsFqtZ0xM2JU4OjrtnjiXzpSp06U4UW0/o5mYyaRh1jTQ/wNO//+77b6dThdmc0c/1Lrwb4sCXCicKK0Rhb6gnXvgnVDnMGIaZhPsMm5mXKUUFZWV9O3Tp0cmrLzQvoW7av5MirPtaG+nCKLK1Be7FoRTM+PC1LIoTWv6P9NgChqdjQSYA+hOOJVpMxlz0y+6tI3NZiMmJoaqqqqz5rUuXeE2NDSQk5PDww8/3Ob1adOmsWHDhna3eeqpp3j88cfPeH3lypWEhoZ26rhHT25iu8kBgQCnteM0Q9MbrVQtlPt2xX9xlQ8MCKOBqZ0ZbJxd2YcJPDUWmXKZUc5wVGM4roZYnKf64awZSnVDHPowMJUeiqQjGoeqq3pkTzkM5i2eZIppO9eZNzDetIcErYJg7AS7eslnr5tNvlccGcDy5cu7tE1dXV2n1utSwi0rK8PpdBIfH9/m9fj4eIqL2/9HeOSRR1i4cGHL381XuNOmTev0FW5xSRJD8ldw6PARBqQN0q9yNRNoJrSAYFRwJARH6M3AmrT3La918euuK/vo8PUuXG10Zh+NjY1s3bqV7Oxs/Rvczcc71/odlmmHL7e/D4Wizu6kot5BVZ2DRqVQLj2ZKpfCpVqvm5v34XK5yMvLY8iQIZg66GTQcWmc+Y5ZM2PWAjBrFgJMlqbnAQRolpbXPT7deSc5nU6+/fZbsrKymn4F9pRRwPfZApgb6wiylxF0qgyTqwFNOdGU6zuP3jEMp8vl4tDhw6QNGIDJ1PXbVIP6pTNgRNeqhmy2zt2L6FbblTN+yirV4QcxKCiIoKAzL2MsFgsWS+eabKUkjyUh/kKWVy7nyokzOr2dL3I4HFTnVjM5ebJflwPoZbH8pJ0ZE66QsnA44Ph2ZlyQ5OaySHfjvnuGw+GgaPlyRl3luVzR2eN0Kf3HxMRgNpvPuJotLS0946pXCCFEW11KuIGBgWRnZ7Nq1ao2r69atYpJkyb1aGBCCOFrulylsHDhQu644w7GjBnDxIkTee211zh69Cj33nuvO+ITQgif0eWEe8stt1BeXs6vf/1rioqKyMjIYPny5fTv37mpRJqbEnW2krmZw+Ggrq4Om83m1/V1Ug6tpCxaSVm0MqIsmvPZuZp9drkd7vk6fvw4KSkyA4MQwvccO3aM5OTkDt/3eMJ1uVwUFha223D/bJqbkx07dqzTzcl8kZRDKymLVlIWrYwoC6UU1dXVJCUlnbUpmseHNDKZTGf9BjiXiIgIv/9AgZTDd0lZtJKyaOXpsoiMPPcYGr43WpgQQngpSbhCCOEhvSbhBgUF8dhjj7Xba82fSDm0krJoJWXRypvLwuM3zYQQwl/1mitcIYTo7SThCiGEh0jCFUIID5GEK4QQHiIJVwghPKRXJFyvniXYTRYtWoSmaW2WhISElveVUixatIikpCRCQkKYMmUKu3btMjDinvPFF19wzTXXkJSUhKZpfPTRR23e78y52+12HnzwQWJiYggLC+Paa6/l+PHjHjyL83eucrjrrrvO+IxMmDChzTq+UA6gT9U1duxYrFYrcXFxzJo1i3379rVZpzd8Lrw+4TbPEvxf//VfbNu2jYsvvpjp06dz9OjRc2/cy40cOZKioqKWJTc3t+W9Z599lueff54XX3yRzZs3k5CQwNSpU6murjYw4p5RW1tLVlYWL774Yrvvd+bcFyxYwIcffsiSJUv46quvqKmpYebMmTidXZoxzVDnKgeAq6++us1n5PS5uHyhHADWr1/PAw88wMaNG1m1ahWNjY1MmzaN2tralnV6xedCeblx48ape++9t81rw4YNUw8//LBBEXnGY489prKystp9z+VyqYSEBPX000+3vHbq1CkVGRmpXnnlFQ9F6BmA+vDDD1v+7sy5V1ZWKovFopYsWdKyTkFBgTKZTGrFihUei70nnV4OSik1b948dd1113W4jS+WQ7PS0lIFqPXr1yules/nwquvcJtnCZ42bVqb1882S7AvycvLIykpibS0NG699VYOHjwIwKFDhyguLm5TLkFBQVx66aU+Xy6dOfecnBwcDkebdZKSksjIyPC58lm3bh1xcXGkp6fzox/9iNLS0pb3fLkcqqr02YmjoqKA3vO58OqE251Zgn3F+PHj+etf/8pnn33G66+/TnFxMZMmTaK8vLzl3P2xXDpz7sXFxQQGBtK3b98O1/EF06dP55133mHNmjU899xzbN68mcsvvxy73Q74bjkopVi4cCGTJ08mIyMD6D2fC48Pz9gdXZkl2FdMnz695XlmZiYTJ05k0KBBvP322y03RvyxXJp159x9rXxuueWWlucZGRmMGTOG/v3788knnzB79uwOt+vt5TB//nx27NjBV199dcZ73v658OorXJkluFVYWBiZmZnk5eW1tFbwx3LpzLknJCTQ0NBARUVFh+v4osTERPr3709eXh7gm+Xw4IMPsmzZMtauXdtmXO3e8rnw6oQrswS3stvt7Nmzh8TERNLS0khISGhTLg0NDaxfv97ny6Uz556dnY3FYmmzTlFRETt37vTp8ikvL+fYsWMkJiYCvlUOSinmz5/P0qVLWbNmDWlpaW3e7zWfC4/cmjsPS5YsURaLRb3xxhtq9+7dasGCBSosLEwdPnzY6NDc6qGHHlLr1q1TBw8eVBs3blQzZ85UVqu15byffvppFRkZqZYuXapyc3PVbbfdphITE5XNZjM48vNXXV2ttm3bprZt26YA9fzzz6tt27apI0eOKKU6d+733nuvSk5OVqtXr1Zbt25Vl19+ucrKylKNjY1GnVaXna0cqqur1UMPPaQ2bNigDh06pNauXasmTpyo+vXr53PloJRS9913n4qMjFTr1q1TRUVFLUtdXV3LOr3hc+H1CVcppf70pz+p/v37q8DAQDV69OiWpiC+7JZbblGJiYnKYrGopKQkNXv2bLVr166W910ul3rsscdUQkKCCgoKUpdcconKzc01MOKes3btWgWcscybN08p1blzr6+vV/Pnz1dRUVEqJCREzZw5Ux09etSAs+m+s5VDXV2dmjZtmoqNjVUWi0WlpqaqefPmnXGOvlAOSql2ywFQb731Vss6veFzIePhCiGEh3h1Ha4QQvgSSbhCCOEhknCFEMJDJOEKIYSHSMIVQggPkYQrhBAeIglXCCE8RBKuEEJ4iCRcIYTwEEm4QgjhIZJwhRDCQ/4fY2GLzUQLQiQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4, 2))\n",
    "plt.plot(u6[56, 0, 0, :], label=\"1.5h\")\n",
    "plt.plot(u6[112, 0, 0, :], label=\"3h\")\n",
    "plt.plot(u6[465, 0, 0, :], label=\"12.5h\")\n",
    "\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAADFCAYAAACW0gNvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2EElEQVR4nO3dd3wVVd748c9tuek3vZFKCcVQQ1c6RCn2VVZYQVdsCwqiu6vrPj/xkV18dtV1XWRRRBRXQRcVUFGJ0ntIggRCCQkkAVJIvem5ZX5/DAmEEEife2/O+/Wa1713MnPnO4fJl5Mz55xRSZIkIQiCINgktdIBCIIgCE0TSVoQBMGGiSQtCIJgw0SSFgRBsGEiSQuCINgwkaQFQRBsmEjSgiAINkyrdADNYbVauXjxIh4eHqhUKqXDEQRBaDNJkigrKyMkJAS1uun6sl0k6YsXLxIWFqZ0GIIgCO0uOzub0NDQJn9uF0naw8MDkE/G09Oz2fuZTCa2bt1KXFwcOp2uo8KzC6IsrhBlIRPlcIUSZWE0GgkLC6vPb02xiyRd18Th6enZ4iTt6uqKp6enuAhFWdQTZSET5XCFkmVxsyZcceNQEATBhokkLQiCYMNEkhYEQbBhdtEmbTNqyuDkd5B9ECqLoKoYLCZQqUGtBpUG1BpQ60DvDnoPeXE2gEcIeIaAdyQYwuTtBUEQbkIk6eYqSIOPZkB5btu/y8kd/PtAt1gIHwkRo8EjqO3fKwiCwxFJujmMOfDxXXKCNoRBzP1yrdjZC7R6kCxgtV5+tYClFmor5Jp3TZlc4y67CKUXoPgc1JbDhcPycug9+RhB/SH6DugzHYIHgRi0IwgCIkk3z8GVcpL1i4ZHvwc3v9Z/l8UEhemQdwyyD0HWfshNubLs+rt8nAEz5cVLDOIRhK5MJOmbsVrh2Jfy+wkvty1BA2h0ENBHXvr/Sl5XUQBp8XBqC6RthYLTsO01eYkaC8PmQe/poBH/XILQ1Yjf+ps5fwhKs8HJA6Jv75hjuPnBoIfkpboUUjfD0c/h3G44u0tePLtB7KMQOxfcAzomDkEQbI7oYnAzKRvk1753gs6l44/nbIAhD8Mj38KiFBjzPLj6gfECbF8Kb/WDDY9BdkLHxyIIguJEkr6ZrAPya5/pnX9sr3CY9P9gcSrc+z6EDgOrCY5tgNWTYdVEOPpfMNd2fmyCIHQKkaRvxGqFwjPy+4C+ysWh1cPAmTDvJ3hiBwycBRonuJAIX82Dt/vDzr9D+SXlYhQEoUOIJH0jxgtgrpIHp3hFKB2NLGQw3PtveO64fCPTPVDuGrh9KfzjFtj4O8g5qnSUgiC0E5Gkb6TgtPzqE2V7PSvcA2DcH2DRMbhvFYQMAUsNHPkU3hsDa6bJNyCtFqUjFQShDWws89iYuqYO317KxnEjWicY8CD0fwDOJ8CBf0PqJsjcKy+GcBj+OPR/SOlIBUFoBZGkb6QgTX71s+EkXUelgrDh8lJ6AQ6vhsNroDQL4v8H7Y5lDDCMhIJeENxP6WgFQWgmkaRvpNCOkvTVDN3kXiFjfw8p/4UDK1HlHyeqYBu8Nwp6TIQRT0PPyWKiJzuRWpjKJ6mfcDjvMGW1ZXjpvQh0DSTYPZho72j6ePch2icaP5c2DrYSbI5I0jdSV5O25eaOG9G5wJA5MPhhzOk7uPTNawSVJqFK3wbp2+SboUMehkGz5blIBJu0+/xuntn2DBbpyv2FClMFF8ovQD58x3f16/1c/Ojt05ve3r3p49OH3t69CfcMR6sWv+r2SvzLNaW2Qu7dAfZXk76WSoUUcRuHui9k2uhb0CWtgaRPoCQTti2F7X+FXnFyQu8VJw9dF2zCmeIzPL/zeSyShTHdxvBozKP4ufhRWlNKbmUu2cZsThWf4lTRKTKNmRRUFVBwoYC9F/bWf4deo6enV0+ivaPp7dObaO9oor2jcVW7KnhmQnOJJN2UsstTkurcwNVH2Vjak1cE3P4Xufte6iZIWgtZ++D0D/LiHggDH5KXgD5KR9vlfXDsA6rMVYwIGsE/J/wT3Q3+A600VZJWksapIjlpnyw+SVpxGlXmKo4XHud44fEG2we5BmGoNZDxSwZ9/foS7R1NhEcEGrWmo09LaAGRpJtSni+/Ouo8GU6uV+YLKUiD5E/gyGdQngd735aXoAHyTHz9fyXmu1ZASXUJ8efiAXgu9rkbJmgAV50rA/0HMtB/YP06q2QluyybU0WnOF18mlPFpzhddJqLFRfJrcwll1xOHT9Vv31drfvqGndPr554O3t3zEkKNyWSdFPK8+RXR03SV/PrBVP+Fyb+j1ybTv4UzsRD7lF5if8fiBonJ+y+M+SnzQgdbnP6ZmqttfT16Us/39b1yFGr1ER4RhDhGUFcZFz9+rLaMlIvpbJx70Z0ITrOlJ4hrTiNakv1dWvdPs4+9PTqSU+vnvTw6lH/atAb2nSOws2JJN2UistDrLtCkq6j0ckTSfW9EyoKIfVrOPqF/LiwjO3y8q2z3Cuk393yQwqcPZWO2mF9d1a+IXh/r/tRtfNDIDycPBgSMIRcfS7TRkxDp9NhsVrkWnexXOs+XXSa08VyrbuouohDuYc4lHuowff4u/g3SNp1rx5O4j/y9iKSdFPqatJuXShJX83NV57Hetg8KMqQZwM8+rk8wOfkt/KicYIek+SE3XsquHgpHbXDqDRVcqpIboYYFzauU46pUWuINEQSaYjk9sgr0/JWmipJL0nnTMkZ+bVUfs2tyOVS1SUuVV3iQM6BBt8V6BrYIHFHGaKIMkSJmncriCTdlPo26UBl47AFPt3lIehjfy8/USZ1ExzfKPcjP/29vKh10H089LsLet0OHqLc2iK1MBWLZCHANYAgN2XvB7jqXOnv35/+/v0brC+vLSe9NJ30knTSitNIL5Hf51flk1eZR15lHnsv7m2wj4+zD5GekfVJu24JcQsRNyybIJJ0U+qTtL+ycdgSlUp+FmNQf7l3yKWTVxL2pRNyO/YZ+UYX3WLl5pDoO+TtxTMbW+RogTxJ1tU3AW2Nu5N7oxuVAKU1pWSUZlypeZec4WzpWfIr8ymqLqKouoik/KQG+zipnQj3DK9P2pGekXQ3dCfSEImbzq0zT8vmiCTdlApRk74hlUqevjWgL4x/ES6dkid0OvUdXEyWp1G9kAjb/yI/VSb6djlhR43tnIcn2Lmjl+QkPcBvgMKRtJxBb2BwwGAGBwxusL7CVME54znOlp6tX84Zz5FZmkmttZYzJWc4U3Km0fcFuAYQ5RlFpEGugUd4RhDuEU6Ie0iXGKTj+GfYWnVzM3fVNumW8u8N434vL2W5cPpHuadI+nZ5UNDhD+VF6wwRt8pD03tMgIB+opZ9DUmS+OXSLwAM8Le/JN0UN50bt/jewi2+tzRYb7FauFhxUU7apec4a7ySxIuqi8ivzCe/Mp+DuQcb7KdVaenm0Y1wj3A5cXuGE+ERQZhnmEM1n7QpSS9btow//elPLFy4kLfffvum2+/du5dx48YRExPDkSNH2nLojiVJV3XBE80dLeYRJD+LMXYumKrg7O4rg2WMFyD9Z3kBcA+Sk3WPiXKbdlfqTdOEvMo8CqoK0Kq0re56Z080ag1hHmGEeYQxNnRsg5+V1pQ2qH2fKz1HVlkW2WXZ1FhqyDRmkmnMZPeF3Q3206q1hLqHNkreEZ4RBLkG2VUCb3WSTkhI4P3332fAgOb9T19aWsqcOXOYNGkSeXl5rT1s56gxynMzg6hJt5XOBaLj5EV6U27Hrps75Nxe+YEFv6yTF5Dbr7tPkPtlh4/okn2yM0oyAIjwjMBZ66xwNMoy6A3Xbfe2SlbyK/Prk3SWMYvMskyyjdlkl2VTa63lnPEc54znGn2nTq0jzCOsPnmHe4YT4hpCsbUYi9WCDtuaFqFVSbq8vJzZs2ezatUqli5d2qx9nnzySWbNmoVGo2Hjxo033Lampoaampr6z0ajEQCTyYTJZGp2nHXbtmQfAEouoAMkJ3fMKh20dH8b1OqyaG/ePWFoTxj6BJirUZ0/hCpjB+qM7ajyUiD38rLvHSSVBil4EFLEaKTwW5HC2idp20xZNOFMsdwuG+4R3qEx2no53Iyvky++fr4M8RvSYL3FaiGvKo/ssmyyyrLIMmaRXS6/P19+HpPVREZpBhmlGY2+819f/Itu7t3kJO4RLtfw3eX3ga6B7VoDb265qyRJklr65XPnzsXHx4d//OMfjB8/nkGDBt2wuWPNmjWsWLGC/fv3s3TpUjZu3HjD5o4lS5bw6quvNlr/2Wef4era8ZPC+Jaf5La0v1KuD+Tnfn/v8OMJMieTEf+yYwSUHce3/CRutQ2f2WhFTalrJAXufShw70uRey/MGsebJGhz5WYO1R5irH4scS5xN99BaDarZKXEWkKhtZAiaxGF1kIKLYUUWgvlmjRNP8lIgwZvtTe+al98Nb74qH3k92pfDGoDGlXLEnhlZSWzZs2itLQUT8+mB4W1uCa9fv16kpKSSEhIaNb2aWlpvPjii+zevRuttnmHe+mll1i8eHH9Z6PRSFhYGHFxcTc8mWuZTCbi4+OZMmUKOl3z/4RRpdZCGrgGRDFt2rRm72fLWlsWSjKVnkeVtQ915l5UmXtRl5zDuzID78oMeuVvQULuYWINHY50ecEr4qY3Im29LDb9vAnyYPLgyUzr3nHXn62XQ2cymUz8uPVH+t/an5yqHLnmfbkGnl2WXV8DL7AWUGAtAHPD/cPcw9h016YWHbOuheBmWpSks7OzWbhwIVu3bsXZ+eZtZRaLhVmzZvHqq68SHR3d7OPo9Xr0en2j9TqdrlUXU4v3qy4CQO0RiNrBLt7WlqEi/KLkZchs+XPpebkd+9xuyNyLqigD8lPR5KdC0kfyNm4Bclt22AgIGwnBA+SnrV+HrZZFpjETgJ4+PTslPlsth86mVqmJ8Iqgp3/PRj+zWC3kVebVN59c2wbezaNbi8uwudu3KEknJiaSn59PbGzsleAtFnbt2sXy5cupqalBo7lS5S8rK+Pw4cMkJyezYMECAKxWK5IkodVq2bp1KxMnTmxJCJ2jbiCLuGloWwyhMHCmvACU5cH5Q5B1ALIPyf2zK/LhxDfyAqDRy09YDx0qD7DpNgTcbPcBB+W15eRXyddfpCFS2WCEehq1hhD3EELcQxgZPLLBz6ySlXJTeYcdu0VJetKkSaSkpDRY9+ijj9KnTx/++Mc/NkjQAJ6eno22X7FiBdu2bWPDhg1ERUW1MuwOVt/9TgxksWkegVcmhAIwVUPOkctJ+6C8VBZC9gF5uUzr6ssIbSjqXccgfLj8pHU3X2XO4Rp1vRH8XPzEJEV2Qq1S4+nUcRONtShJe3h4EBMT02Cdm5sbvr6+9etfeuklLly4wNq1a1Gr1Y22DwgIwNnZudF6m1I/A57oI21XdM4QPlJeQO7vXpguJ+uLSfIIyNxjqCoLCaIQdv9yZV/vSDlZd4uVl+CB8pzbnexs6VkAogw2WoEROl27jzjMyckhKyurvb+2c3X1GfAchUoFfj3lZfDldm1TNeYLR0j96T/EeNegvpgsTxRVfE5ejn91eV+NPBqy2+XEHToU/PtABw+CqEvSkZ6RHXocwX60OUnv2LGjweePPvrohtsvWbKEJUuWtPWwHatuSLho7nA8OmekbrGc9c+j77Rp8o3hqpKr5hu5XOMuz4W8FHlJ+lje19lLHh0Z8yt5HhJN+8+qUNfcIWrSQh0xd8e1JOmqyZVEc0eX4OJ1eWj6hCvrSi9cmSTqQqKcxKtL4PjX8uLTQ36aTd8Z7RqKaO4QriWS9LWqS8BSK78XzR1dl6GbvPS7S/5sMcvt2ic2y8+CLEqHz2fDiKflB/u2QzOIxWqp734nkrRQR610ADanrqlDb5BvRAkCyE0bYcMhbiks/AVGPyuvP/hv2PrndjnExfKLmKwm9Bo9wW7B7fKdgv0TSfpaXekBtELr6D0g7jW4933584EVkLS2zV971ig3dUR4RqBWiV9NQSauhGvVt0eLJC3cxMCZMOFyLXrrn6GioE1fJ9qjhesRSfpa5SJJCy1w23Py9KrVpfDTkjZ9leh+J1yPSNLXEkPChZbQaGHam/L7I5/J84u0kqhJC9cjkvS1RPc7oaXCR0DkGJAscGhVq77CbDVzougEANHezZ+MTHB8Iklfq1w8gFZohZFPy6+JH0FtZYt3Ty9Jp8pchbvOnR5ePdo3NsGuiSR9rcJ0+dUQpmwcgn2JvkOey7q6BE5tafHudQ+ejfGLET07hAbE1XC12goouvxInUAbngBKsD1qDfT/lfw+dWOLd3fEp4ML7UMk6avlnwQk+aahaJMWWqrfPfJrWjzUtGx+4aOXjgI0euCqIIgkfbW8Y/Jr4C3KxiHYp6D+4NMdzNWQ9mOzdyuqLqqfWGmAn6hJCw2JJH21vOPyq0jSQmuoVNDvbvn9yea3S/9w9gcA+vj0wcvZqwMCE+yZSNJXq0/Soj1aaKWek+XXjB1gtTZrl41nNgJwT897OiQkwb6JJF1HkiC/Lkn3UzYWwX6FDgedG1QWXGk+u4EThSc4UXQCnVrH9KjpnRCgYG9Ekq5zMQmqikGtA7/eSkcj2CutE0TeKr/P2H7DTa2SlTcOvwHAxPCJoqlDuC6RpOvs/af82v9XYopSoW26X354QPqNk/TqlNUcyj2Ei9aFZwY/0wmBCfZITPoPkH8CUjfL70eLXxahjeqe8JK1H8y1cu36KgVVBSxPXs6XaV8C8MLQF4jwjOjsKAU70XWTtNUiPxU8fRv8+DIgQa840bNDaDWL1UKttRaTZzAmNx9M1aVUpv9AoSGYwqpCzhnPcST/CIdyD2GymlChYv6g+TwQ/YDSoQs2TCVJkqR0EDdjNBoxGAyUlpbi6enZrH3Ssn5h0Y+/ASRU9WvlU1UhocaKiiunXqPSU6D2x6rSoALUahUaNeg0avRaeVGr5W9qSZFJXH/bpr7jetu3ZNsm45AkKisrcXF1QdWgRFr23deLpck4mlzd/O9or3O/dt/a2lqcnJwab9vCOKySFZPVhMlqwio1rzcHyCMLFwxawKiQUc3ep72ZTCa2bNnCtGnT0Ol0isVhC64ui0ozHMwo5HReGRdLq8ktraay1ozJImGyWKk1W7n2cugZ6M67s4a06JjNzWsOW5OuNdWS5QRclZAavr/2mXQSkN9wlfXyYmr/+JRSUlGidAg2o7Km5RMhtYROknBRafHxDMPH2YdQj1B6e/dmdMhoenj1QKVS3fxLhE6TVwWL/3uUH47nYbK0rO7qpO2423sOm6TDg3rxctAjnMs6T3hEBCq1FqtKrkNLKh0mJ0/MWnesKg2SJPfAk/9ZJGotVsqrzJRWm7hYUk36pQoKymvqv9tVp+GewaFM7R+Ek+bKP871fulUXP8XsSXbNnf/pr5DhQqz2cy+ffsYPXo0Wq225TFcZ3VTx2puvB15vjf6DrPZzK5duxg7diw6bRM1yBacr5PGCSeNEzq1Tl40OrTZh1GtuQPc/GHOZnmgi2CTJElizb5M/vaLBrOUC0B3PzcGhXsR6u1KsMEZd70WnUaNk1aFVq1Go77y76kC3PQdl0odNkl7uHlx/8Rn5T9hJrT9z7mc0iq2pOTy6cFMMi5V8MkO2HG0nD9P78fttwS1S8wdyWQykaXNIsYvRvxpazJxUnOS7obuHVcWIYNB4yTf9yjKAF8x/agtkiSJN7ae4t3t6YCKsb18+f3tfekfalA6tHqiC14zBRtceOy2KH56bhxvPjCQQE892UVVPPlJIi99dZSqWovSIQq2ROcMwYPk9+cTFA1FaNpnh7IuJ2i4K9zCBw8PsakEDSJJt5hareL+2FC2vzCeJ8d1R6WCdYeyufvdPZzJL1M6PMGWhA6VXy8kKhuHcF1peWW89m0qAM9P7smkbpJN3icQSbqVXJ20vDS1L5/8dgR+7npO55Vz77v72J12SenQBFvRLVZ+PX9Y2TiERiRJ4uWNx6g2WRnTy48nxtjucyVFkm6j23r58f3CMQyP9KGsxswjaxL47GCW0mEJtqAuSeemgLnmxtsKnWrHqUscOluEXqvm/+4fUN+91haJJN0O/D30fDJvOPcN7obFKvGnr1N4K/50i/pTCw7IOxJcfcFqkhO1YBMkSeJvP54C4JHRkYR4uSgc0Y2JJN1O9FoNbz44kOcmy096fufnNF79JhWrVSTqLkululKbFu3SNmN/RiEncoy4OWl4erzt97ppU5JetmwZKpWKRYsWNbnNV199xZQpU/D398fT05NRo0bx44/Nf2qFPVGpVCyc3IvX7paHln+07xy/33AUs6X5I9EEB9Pt8s1D0S5tM/5zIBOAe4d0w8u18ahTW9PqJJ2QkMD777/PgAE3ftzPrl27mDJlClu2bCExMZEJEyZw5513kpyc3NpD27yHR0Xyj5kD0ahVfJl0nvmfJVFjFl30uiRRk7YpecZqfjyeB8BvRtrHpFatGsxSXl7O7NmzWbVqFUuXLr3htm+//XaDz3/961/ZtGkT33zzDYMHD77uPjU1NdTUXLnRYjQaAXkQgsnU/DHaddu2ZJ/2MiMmEGfNQJ79/Bd+PJ7Hb9cksGLWQFydlBk/pGRZ2JpOLYvAAegAitIxGfPBxbvjj9lMXfGa+CoxG4tVYki4Fz18XRqVQWeWRXOP1aqMMX/+fKZPn87kyZNvmqSvZbVaKSsrw8fHp8ltli1bxquvvtpo/datW3F1dW1xvPHx8S3ep7083lvFByfV7E0v5O5//MyTfS24KjjOU8mysDWdVRaT9IG41+SRsOk9Lnna3oNmu9I1se6oBlDRXVPIli2Nn0PZmWVRWdm8uWNanC7Wr19PUlISCQmtG0X15ptvUlFRwYMPPtjkNi+99BKLFy+u/2w0GgkLCyMuLq7Zs+CB/D9VfHw8U6ZMUWwo9DRgQnYJj61N4ly5mbXZ3qyZOwRfd32nxmELZWErOrssNKbNcGwDI7ppsI6Z1uHHa66udk1kF1eStX8PahU8/+BE/K76HVSiLOpaCG6mRUk6OzubhQsXsnXrVpydW/70knXr1rFkyRI2bdpEQEBAk9vp9Xr0+sZJTKfTtaoAW7tfexnW3Z/PnxzFw6sPciK3jNkfHuY/80YQbOj8rj9Kl4Ut6bSyCB0GxzagyTmCxgbLvqtcE1tPFAAwsrsvwd7u192mM8uiucdp0Y3DxMRE8vPziY2NRavVotVq2blzJ++88w5arRaLpembY59//jmPPfYYX3zxBZMnT27JYR1C32BPvnhyFCEGZ9IvVfDAyv1kFlYoHZbQGeqHhx+m0UTEQqeJT5VvGE7tH6xwJC3ToiQ9adIkUlJSOHLkSP0ydOhQZs+ezZEjR9Borp2jWbZu3ToeeeQRPvvsM6ZP77pPRO7u784XT40i0teV88VVPLByP2l5Yr4PhxcYIz/guLIQSjKVjqZLKqmsJTmrGICJfZr+K94WtShJe3h4EBMT02Bxc3PD19eXmJgYQG5PnjNnTv0+69atY86cObz55puMHDmS3NxccnNzKS0tbd8zsROh3q588eQoegd6kF9Ww4Pv7SflfNcsiy5D5wxB/eX3or+0InanFWCVoFeAO91sfIThtdp9xGFOTg5ZWVfmrnjvvfcwm83Mnz+f4ODg+mXhwoXtfWi7EeDpzPonRjIw1EBxpYlZqw6QcK5I6bCEjlTfXzpJ2Ti6qB2n5InPxvf2VziSlmtzZ7AdO3Y0+PzRRx/d8OeCzNvNif/MG8FjHx/m0NkiHl59kPceHsq4aPu7iIRmCB0KCavkdmmhU0mSxK7Ls1OOi7avpg4Qc3coysNZx8ePDmdctD/VJiuPfZTA18nnlQ5L6Ah1NemcX8DSdQaP2IIz+eVcKqtBr1UzLMp2BhM1l0jSCnNx0rBqzlDuGhiC2Srx3Oe/sHJnuphBz9H49ABnA5irIe+40tF0KQcyCgEYGumNXnv9zg22TCRpG+CkVfP2zEE8fnni8de/P8mr36RiETPoOQ61GkKGyO/FPB6d6kCGfL9nZJSvwpG0jkjSNkKtVvHy9H78eXpfQJ5B75l1SVSbxMRMDkM8TqvTSZJUX5Me2UMkaaEdzBvTnXceGoxOo2JLSi5zPjxEcUWt0mEJ7UHMiNfp0vLLKayoxVmnZmCol9LhtIpI0jboroEhfPzocDz0Wg6dLeK+f+8j41K50mEJbVWXpC+dgurmzdsgtM2hs3JTR2yEN05a+0x39hl1FzC6px8bnh5NNy8XzhZUcO+KfexPL1Q6LKEt3APAEA5IcNFx51O3JUmXRxnGhttfr446IknbsN5BHmycfyuDwrworTLx8OqDfHE4W+mwhLYIrWvyEP2lO0NyVgkAgyNEkhY6iL+HnvVPjGTGgGDMVok/bDjK69+fFM9OtFdi5GGnKaqo5WyBPInZkDCRpIUO5KzT8M6vB/PMxJ4ArNyZzuNrD2OsFoMi7M7VzzwUfeE7VN2ESj383TC42u9UrCJJ2wm1WsXzcb35x8yBOGnV/Hwyn3ve3Uu6uKFoX4IHgkoD5blgvKh0NA6trj16iB23R4NI0nbn3sGhbHhqFMEGZzIuVXDP8r38fCJP6bCE5nJyhcB+8nvRLt2hkjJLABhix+3RIJK0XRoQ6sXmBbcxLNKbshoz89Ye5l8/p4l2ansh+kt3OLPFypHsEkDUpAWF+Hvo+XTeSB4eGYEkwZvxp/ndp0mU15iVDk24mfp2aZGkO8rJ3DKqTBY89Fp6BVz/UVn2QiRpO+akVfPaPTG8fl9/nDRqfjiey72indr21Q0Pv5gEZjGatCPU3TQcFO6FWq1SOJq2EUnaAfx6eDjrnxxJgIeetPxy7vrXHr49Km5K2Sy/3uDqB6ZK0S7dQZIu94+296YOEEnaYQwJ9+bbZ29jZHcfKmotLPgsmSWbj1NrtiodmnAttRqixsrvM3YqG4uDqu/ZYec3DUEkaYcS4OHMfx4bwdPjewDyTHoz39/PxZIqhSMTGuk+Xn7N2KFkFA6poLyGzMJKAAaFeSkbTDsQSdrBaDVq/nhHHz6YMxRPZy3JWSVMf2c3u88UKB2acLXu4+TXC4ehRjwxvj3VDQXvFeCOwcV+B7HUEUnaQU3uF8i3z4whppsnxZUmHlubxPfZKvEgAVvhHSkvVjOc3a10NA7FUQax1BFJ2oGF+7qy4anRPDQ8HEmCH85rmPdJEkVifmrb0CtOfj39vbJxOJikzLr2aC9lA2knIkk7OGedhmX39edv98WgU0vsOVPI9Hd219c2BAX1niq/nvoBrOIGb3swWawcPV8KiJq0YGfuHRzC4hgLkb6u5JRWM/O9/azZe1Y88FZJEbeB3hMq8sXow3ZyMkcexOLprKWHv30PYqkjknQXEuIGXz01kun9gzFZJF79JpUFnyVTJmbTU4bWCXpOlt+f/FbZWBxE3V+Ig8O97X4QSx2RpLsYD2cty2cN5pU7+6FVq/guJYe7l+/lZK54nJMi+t4pvx77UjR5tANHu2kIIkl3SSqVikdvjeLzJy/PpldQwT3v7mVD4nmlQ+t6ek8FvQFKsyFzr9LR2L0rg1i8lA2kHYkk3YXFRnjz3bNjGBvtT7XJygv//YUXvzxKtcmidGhdh84FbrlHfv/LekVDsXeXymrILqpCpXKMQSx1RJLu4nzcnPjokWEsnhKNSgXrE7K5b8U+MgsrlA6t6xj4kPx6/GuoKlE0FHtWV4uODvDAw9n+B7HUEUlaQK1W8eykXqz97XB83JxIzTEy4197+PF4rtKhdQ3hIyGgH5gqIPEjpaOxW47Y1AEiSQtXGdPLn++evY3YCG/Kqs08+Ukif/kuFZNF3NDqUCoVjFogvz/4npi+tJUOn3O8m4bQxiS9bNkyVCoVixYtuuF2O3fuJDY2FmdnZ7p3787KlSvbclihAwUbXFj/xEgeHxMFwKrdZ3no/QPkGasVjszB9f8VuAdB2UVIXqt0NHanqtbC0fMlAIyI8lU2mHbW6iSdkJDA+++/z4ABA2643dmzZ5k2bRpjxowhOTmZP/3pTzz77LN8+eWXrT200MF0GjUvT+/Hyt8MwUOv5XBmMXf+a0/9ROpCB9DqYewL8vvty6C6VNl47ExydjEmi0SwwZkwHxelw2lXrUrS5eXlzJ49m1WrVuHtfeM/LVauXEl4eDhvv/02ffv2Zd68efz2t7/ljTfeaFXAQue5IyaYb565jehAd/LLapj53gG+OJytdFiOK/YR8IuGygLY9helo7ErBzOKABge5YNK5RiDWOpoW7PT/PnzmT59OpMnT2bp0qU33Hb//v3ExcU1WHf77bezevVqTCYTOl3ju7A1NTXU1NTUfzYa5YEWJpMJk6n5o+Pqtm3JPo6qtWXRzeDE548P5w9fHiP+RD5/2HCUY+dLePGOaHQa+7ylYcvXhWrKX9CuewAOvYe5+0SkHpM67Fi2XA4tdTBDnoo3NtzQqvNRoiyae6wWJ+n169eTlJREQkJCs7bPzc0lMDCwwbrAwEDMZjMFBQUEBwc32mfZsmW8+uqrjdZv3boVV1fXloZMfHx8i/dxVK0ti2kG0ISq+OG8hrUHstifeo5Hoq2423FPJ1u9Lvr7T6H7pXgs/32MPdF/ptw5pEOPZ6vl0FxmKySe0wAqqrNS2HIppdXf1ZllUVlZ2aztWpSks7OzWbhwIVu3bsXZ2bnZ+13750fdpD5N/Vny0ksvsXjx4vrPRqORsLAw4uLi8PT0bPZxTSYT8fHxTJky5bo19q6kPcpiBhCfms/vv0whzQgr09344OEhdPd3a99gO5jNXxfmiVjX3ok+J5mJF/6FedaX4Nuz3Q9j8+XQTImZxZgOJuDr5sQj901pVXOHEmVR10JwMy1K0omJieTn5xMbG1u/zmKxsGvXLpYvX05NTQ0ajabBPkFBQeTmNuxvm5+fj1arxdf3+ndh9Xo9er2+0XqdTteqAmztfo6orWUxbWA3egZ5Mu/jw2QVVfLgqkOsmjOU4VE+7Rhl57DZ60Kng998CR/ejqowDd3H0+BXH0KPCR10OBsth2ZKzJaT3YjuPjg5ObXpuzqzLJp7nBY1Kk6aNImUlBSOHDlSvwwdOpTZs2dz5MiRRgkaYNSoUY3+hNi6dStDhw616wujK4sO9OCr341mUJgXpVUmfvPBQTYduaB0WI7FzRce3QIhQ6CqCD65B75ZCNViIqxrHTx7+aZhpP1VFJqjRUnaw8ODmJiYBoubmxu+vr7ExMQAclPFnDlz6vd56qmnyMzMZPHixZw4cYIPP/yQ1atX88ILL7TvmQidys9dz/onRjI1Johai5WF64+wYscZMT91e3IPgEe+g2Hz5M+JH8GKkZD8KVjF/CoAZouVxHNykh7R3bH6R9dp99vzOTk5ZGVl1X+Oiopiy5Yt7Nixg0GDBvHaa6/xzjvvcP/997f3oYVO5qzT8O6sITwxtjsAf/vhFK9/f1Ik6vbk5ArT35STtXcUGC/Apt/BilGQuhm6eFmn5hipqJUn+e8d6KF0OB2iVV3wrrZjx44Gnz/66KNG24wbN46kpKS2HkqwQWq1ij9N60uAh56l353gvV0ZGKtNLL2nPxoHmXTdJkTeBr/bD4dWwZ63oOAUfPGw3Bwy+RXoPl7pCBWxO03uejc8ytdhJvm/ln12dBVszrwx3fnb/QNQq2DdoWwWrk+m1izm/GhXOhe49VlY+AuM/QPo3OBiEqy9W14udL2K0M7TlwAY19tf4Ug6jkjSQrt5cFgYy2cNQadR8e3RHJ7+T6JI1B3B2QATX5aT9YinQa2DjB2wagJ8MQcK0pSOsFMYq031TwYf10skaUFolmn9g/lg7jCcdWp+PpnPM+uSxCx6HcXdH6a+Ds8kXp6TWgWpm+DdEfDd81BZpHSEHWrfmQLMVonufm6E+7Z8kJu9EElaaHfjov1ZNWcoTlo1Px7P47nPj2AWibrjeEfAvSvh6X3QexpIFkj4AP41RH510J4gO045flMHiCQtdJAxvfx57zex9U0ff9hwFIu1a/dE6HCB/eChdTD3Wwi4BaqK5Rr1e2Mh64DS0bUri1XipxN5AEzoHaBwNB1LJGmhw0zoE8DyWUPQqFV8lXyBP29MEd3zOkPUGHhyF0x7A5y9IO8YfHgHfPeCwwyGOXyuiILyWgwuOkb1cMz+0XVEkhY61O23BPHPXw+q7/Xx9k9d46aW4jRaGP44PJsMg38DSJCwSh4Mc+oHpaNrs++PyVNNTO4baLezMTaXY5+dYBNmDAjhtXvkEan//DmNzw5m3WQPod24+sDd78KcTeAdKQ+GWTcTvnrSbmvVVqvE1svP35waE6RwNB1PJGmhU8weEcGzE+WZ3P68MYX41DyFI+piuo+Hp/fD6GdApYaj69F+MA6f8lNKR9ZiB84WcrG0Gg+9ltt6+SkdTocTSVroNM9NiebBoaFYJXhmXRKJmeJxXJ3KyRXilsKj34NXOKrSbG5L+yvq7X8Bi/1M/L/h8HkAZgwMwVnXeFI3RyOStNBpVCoVf7m3PxN6+1NtsjLv4wTSL5UrHVbXEz4SntqLdcCvUSGh2fcP+cZicabSkd1UWbWJLcdyAHhgaKjC0XQOkaSFTqXTqHl39hAGhhoorjQx98ND5IsnkXc+Z08sdy4nIXIBkrMBLhyGlWPkwTA2bGPyBapNVnr4uzE4zEvpcDqFSNJCp3N10rL6kWFE+rpyvriKR9YkUFZtP39uO5KL3sMxz9sBocOgplQeVv7dC2Cyvf84LVaJD/acBWDOqEiHe+BsU0SSFhTh567n498Ox9fNidQcI0//J0nM86EUQ5jcTn3rQvlzwipYPRkKzigb1zW2Hs8ls7ASL1ddl2nqAJGkBQVF+Lqx5tFhuDpp2HOmgD9+eVQMdlGKRgdT/hdmfwmufpCbIo9U/OVzpSMD5Fr0P3+W+9j/ZkQErk5tnmXZbogkLShqQKgXK2bLoxK/Tr7A//1gf13CHEqvyfDUHogcA6YK+PoJ2DgfaisUDeuLw9mczC3D4KJj3pgoRWPpbCJJC4ob3zuA1+/rD8DKnen83w/i6S6K8gyWB7+Mf0nuU33kP/D+BMhLVSScwvIa3twq/+f97KReeLm27WGz9kYkacEmPDA0jD9P7wvAv3ek8/82HccqJmRSjloD41+EOZvBPUh+EsyqCXB4Tac+skuSJH6/4SgF5bVEB7rz8MiITju2rRBJWrAZ88Z05y/3xqBSwScHMnl2fTJVtY45zabdiBoDT++FnpPBXA3fLoK1d0H+yU45/BtbT7HtZD5OWjXvPDQYJ23XS1ld74wFmzZ7RARvzxyERi1PcXrvir0cv1iqdFhdm5sfzPovTHkNtM5wdhesvBW+/yMYL3bIIa1Wib//eJJ3t6cD8Nrdt9AnyLNDjmXrRJIWbM7dg7rx2bwR+LnrOZlbxt3L9/LKpmNi0IuS1Gr5+YrzD0Lv6WA1w8GV8PYA2LRAfr5iOzWDnM4r4+EPD9Yn6Ben9mHmsPB2+W571HX6sQh2ZUR3X75fOIZXNh9jS0ouH+/P5LNDWUzpF8iUfoHcEmIg0tetS/75qyjvSHjoM0jfBrvegMy9kPyJvHhFQN87oVssBA+Ut1XffG6NylozZwsqSMoqYevxXPacKUCSwFmn5q/39ue+IV2nT/T1qCQ7uI1uNBoxGAyUlpbi6dn8P3lMJhNbtmxh2rRp6HS6DozQ9tlzWew7U8Cb8acbTcikUavwd9fjrFOj12pw0qrRalTcbByaJEkUF5fg7e3VaNRae41iu/bX6tpfsmt/6xr9ErZ4f+nGP7/O8SRJosxoxKMZv1NNHa+f+QT31m5mhCkBZ2obbGNFRZnKHaPKExNarKjrFwtqLBKYrGC9TgrydnUi1NsFZ23nTKBklSSKi4vw9vZB3ZprwK+nPCVsCzQ3r4matGDzRvf0Y3RPP45dKGVLSg570wvJyC+nrMZMbqubQFScKxdt3aCCyrJW732aMDYyHxeqmaA+wm3qY/RTn6OPKhtnlQmDVIZBusH3qy4v16oGclodVoupAV+A1nYHN3dcU5xI0oLdiOlmIKabAZBrgfllNeQba6gxW6gxW6k1W6ltxgNvLWYLiUmJxA6JRXNVTa1tf1NKXJttrq2QXZuLGtXiG/38Jp+v3aOFx7NYzCQcSmDY8GFoNdqbfv/1KpgNV42Tcz6QZDWjqylGV1OCtqYEtbUWJCtqyQKSFZVkwUmjwk2vxU2vxVWnVnQuDrPZQlJSIkOGxKJtTe3d2dD+QV0mkrRgl1QqFYGezgR6Ord4X5PJhOmcxJR+AXbX9NOeTCYTZaclxvT066BysJ+npkgmEzkZKqQ+08DGrglx10UQBMGGiSQtCIJgw0SSFgRBsGF20SZd153JaGzZ041NJhOVlZUYjcYu3fYIoiyuJspCJsrhCiXKoi6f3awXtF0k6bIyuQtPWFiYwpEIgiC0r7KyMgyGpnuH2MVgFqvVysWLF/Hw8GhRNx2j0UhYWBjZ2dktGgTjiERZXCHKQibK4QolykKSJMrKyggJCUGtbrrl2S5q0mq1mtDQ1g8N9fT07PIXYR1RFleIspCJcriis8viRjXoOuLGoSAIgg0TSVoQBMGGOXSS1uv1vPLKK+j1eqVDUZwoiytEWchEOVxhy2VhFzcOBUEQuiqHrkkLgiDYO5GkBUEQbJhI0oIgCDZMJGlBEAQbJpK0IAiCDXPYJL1ixQqioqJwdnYmNjaW3bt3Kx1Sh1uyZAkqlarBEhR0ZeJ1SZJYsmQJISEhuLi4MH78eI4fP65gxO1n165d3HnnnYSEhKBSqdi4cWODnzfn3GtqanjmmWfw8/PDzc2Nu+66i/Pnz3fiWbSPm5XFI4880ug6GTlyZINtHKEsli1bxrBhw/Dw8CAgIIB77rmHU6dONdjGHq4Lh0zSn3/+OYsWLeLll18mOTmZMWPGMHXqVLKyspQOrcPdcsst5OTk1C8pKSn1P/vb3/7GW2+9xfLly0lISCAoKIgpU6bUT2BlzyoqKhg4cCDLly+/7s+bc+6LFi3i66+/Zv369ezZs4fy8nJmzJiBxWLprNNoFzcrC4A77rijwXWyZcuWBj93hLLYuXMn8+fP58CBA8THx2M2m4mLi6Oi4sqDDO3iupAc0PDhw6Wnnnqqwbo+ffpIL774okIRdY5XXnlFGjhw4HV/ZrVapaCgIOn111+vX1ddXS0ZDAZp5cqVnRRh5wCkr7/+uv5zc869pKRE0ul00vr16+u3uXDhgqRWq6Uffvih02Jvb9eWhSRJ0ty5c6W77767yX0ctSzy8/MlQNq5c6ckSfZzXThcTbq2tpbExETi4uIarI+Li2Pfvn0KRdV50tLSCAkJISoqil//+tdkZGQAcPbsWXJzcxuUi16vZ9y4cQ5fLs0598TEREwmU4NtQkJCiImJccjy2bFjBwEBAURHR/P444+Tn59f/zNHLYvSUvnp8D4+PoD9XBcOl6QLCgqwWCwEBgY2WB8YGEhubq5CUXWOESNGsHbtWn788UdWrVpFbm4uo0ePprCwsP7cu2K5NOfcc3NzcXJywtvbu8ltHMXUqVP59NNP2bZtG2+++SYJCQlMnDiRmpoawDHLQpIkFi9ezG233UZMTAxgP9eFXUxV2hrXzjstSZKij4zvDFOnTq1/379/f0aNGkWPHj34+OOP628MdcVyqdOac3fE8pk5c2b9+5iYGIYOHUpERATfffcd9913X5P72XNZLFiwgKNHj7Jnz55GP7P168LhatJ+fn5oNJpG/8vl5+c3+h/T0bm5udG/f3/S0tLqe3l0xXJpzrkHBQVRW1tLcXFxk9s4quDgYCIiIkhLSwMcryyeeeYZNm/ezPbt2xvMS28v14XDJWknJydiY2OJj49vsD4+Pp7Ro0crFJUyampqOHHiBMHBwURFRREUFNSgXGpra9m5c6fDl0tzzj02NhadTtdgm5ycHI4dO+bw5VNYWEh2djbBwcGA45SFJEksWLCAr776im3bthEVFdXg53ZzXXTK7clOtn79ekmn00mrV6+WUlNTpUWLFklubm7SuXPnlA6tQz3//PPSjh07pIyMDOnAgQPSjBkzJA8Pj/rzfv311yWDwSB99dVXUkpKivTQQw9JwcHBktFoVDjytisrK5OSk5Ol5ORkCZDeeustKTk5WcrMzJQkqXnn/tRTT0mhoaHSTz/9JCUlJUkTJ06UBg4cKJnNZqVOq1VuVBZlZWXS888/L+3bt086e/astH37dmnUqFFSt27dHK4snn76aclgMEg7duyQcnJy6pfKysr6bezhunDIJC1JkvTuu+9KERERkpOTkzRkyJD6bjeObObMmVJwcLCk0+mkkJAQ6b777pOOHz9e/3Or1Sq98sorUlBQkKTX66WxY8dKKSkpCkbcfrZv3y4BjZa5c+dKktS8c6+qqpIWLFgg+fj4SC4uLtKMGTOkrKwsBc6mbW5UFpWVlVJcXJzk7+8v6XQ6KTw8XJo7d26j83SEsrheGQDSmjVr6rexh+tCzCctCIJgwxyuTVoQBMGRiCQtCIJgw0SSFgRBsGEiSQuCINgwkaQFQRBsmEjSgiAINkwkaUEQBBsmkrQgCIINE0laEATBhokkLQiCYMNEkhYEQbBh/x+lBVt4IrHG5wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4, 2))\n",
    "plt.plot(ph[56, 0, 0, :], label=\"1.5h\")\n",
    "plt.plot(ph[112, 0, 0, :], label=\"3h\")\n",
    "plt.plot(ph[465, 0, 0, :], label=\"12.5h\")\n",
    "\n",
    "plt.grid()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "work",
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}