{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f5ac2f81-49ec-48e8-9968-4d5f3c882a61",
   "metadata": {},
   "source": [
    "# Using Intel oneMKL PARDISO as an external sparse direct solver for PyFVTool\n",
    "\n",
    "MW 250721\n",
    "\n",
    "This is an example of how to use an external sparse solver via a PyFVTool interface.\n",
    "\n",
    "An external solver can be supplied as an optional argument to `pyfvtool.solvePDE()`, *e.g.*\n",
    "\n",
    "```python\n",
    "pf.solvePDE(phi, eqnterms,\n",
    "            externalsolver = solveur)\n",
    "```\n",
    "\n",
    "where `solveur` is a function that will be called instead of `scipy.sparse.linalg.spsolve(A, b)`. It should provide the same interface as the original `spsolve`.\n",
    "\n",
    "In this example, we use `pyfvtool.solvers.oneMKL_pardiso.spsolve_oneMKL_pardiso(A, b)`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec93bdd4-096c-46e1-8258-60118a4e29f2",
   "metadata": {},
   "source": [
    "## PyFVTool MKL PARDISO interface\n",
    "\n",
    "PyFVTool includes a minimal interface to the [Intel oneMKL PARDISO solver](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-2/onemkl-pardiso-parallel-direct-sparse-solver-iface.html)\n",
    "\n",
    "Beyond `numpy`, `scipy` and `matplotlib`, it requires that the Intel oneMKL runtime (package `mkl`) be installed in the PyFVTool environment. On Windows, this may already be the case, since NumPy already depends on it, but we found that it needs to be installed explicitly on our Linux system:\n",
    "\n",
    "```\n",
    "conda install mkl\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c827dee8-ae1e-4d35-894a-49bca4c3caf1",
   "metadata": {},
   "source": [
    "## Test drive the solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e77a132-0f37-491c-a58e-6d5a39b0c60e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import sparse\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4ae46233-1a77-4522-b271-1bee7dbbd448",
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "import pyfvtool as pf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "70b1f63a-54e4-4cd3-a23e-105726fcc489",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyfvtool.solvers.oneMKL_pardiso import spsolve as spsolve_oneMKL_pardiso\n",
    "\n",
    "# additional import to get MKL version printed\n",
    "from pyfvtool.solvers.oneMKL_pardiso import oneMKL_pardiso_solver_instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3b468728-8a4e-4614-8d7d-bbdc0da11c74",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python 3.12.9 | packaged by conda-forge | (main, Mar  4 2025, 22:37:18) [MSC v.1943 64 bit (AMD64)]\n",
      "NumPy version          :  2.3.0\n",
      "SciPy version          :  1.15.2\n",
      "oneMKL runtime version :  2025.0.1\n",
      "PyFVTool version       :  0.4.2\n"
     ]
    }
   ],
   "source": [
    "print('Python', sys.version)\n",
    "print('NumPy version          : ', np.__version__)\n",
    "print('SciPy version          : ', sp.__version__)\n",
    "print('oneMKL runtime version : ', oneMKL_pardiso_solver_instance.get_MKL_version_string())\n",
    "print('PyFVTool version       : ', pf.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fd00e51e-2bc8-4e55-abdf-1445977a2c72",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = sparse.rand(10, 10, density=0.5, format='csr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "29c6c0a7-6f4f-4728-ac02-6e6f3b5c3ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.random.rand(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "66613514-0885-456a-932f-c40d10a84297",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = spsolve_oneMKL_pardiso(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6d20b1a3-09af-43d4-9d05-3c93a6980e39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.90190165,  0.7385448 , -2.45810562,  1.8050738 ,  1.65217373,\n",
       "        0.98397367,  1.01222876, -0.27855745, -4.20957481, -3.05200821])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba7f0702-7886-4df5-92fd-585a3af98374",
   "metadata": {},
   "source": [
    "## PyFVTool finite-volume (from cylindrical2D_convection notebook)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df8e1325-198a-47a6-8e71-06b3604a6ed5",
   "metadata": {},
   "source": [
    "#### Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b56afb82-9f77-406f-829c-321b741b3c16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualization routine (imshow-based)\n",
    "def phi_visualize():\n",
    "    print(f't = {t:.1f} s')\n",
    "    # avoid ghost cells\n",
    "    plt.imshow(phi.value, origin = 'lower',\n",
    "               extent = [zz[0], zz[-1], rr[0]*rzoom, rr[-1]*rzoom])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1cd0611d-f786-4fa2-a278-9b8cb158b03b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate simple finite-volume integral over r\n",
    "def integral_dr(phi0):\n",
    "    v = phi0.cellvolume\n",
    "    c = phi0.value\n",
    "    return (v*c).sum(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44b77e8a-d5fb-4229-8d1b-741447fb32fd",
   "metadata": {},
   "source": [
    "#### FVM settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b70cedd2-baae-40e6-a123-65007b605c93",
   "metadata": {},
   "outputs": [],
   "source": [
    "Lr = 7.5e-05 # [m] radius of cylinder\n",
    "Lz = 0.3 # [m] length of cylinder\n",
    "umax = 2*9.4314e-3 # [m s^-1] max flow velocity = 2 time average flow velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b331b3d2-e52d-4852-8233-5b6794632859",
   "metadata": {},
   "outputs": [],
   "source": [
    "# regular grid parameters\n",
    "Nr = 40\n",
    "Nz = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7a1af730-d1d4-49c1-b567-a82aa5f4014a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial condition parameters (cell indices)\n",
    "loadix0 = 20\n",
    "loadix1 = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2addc814-63d0-4e00-b3b1-f85fc3fd8b25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# timestep parameters\n",
    "deltat = 0.01 # [s] per time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fbd7f090-6a9f-4510-9c64-364ab2e72c43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualization parameters\n",
    "rzoom = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0c228f5-ad35-40e0-817e-55ea640b5cfe",
   "metadata": {},
   "source": [
    "#### 2D cylindrical mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "7831e5d0-d129-482d-9b26-590cc5f97a5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "msh = pf.CylindricalGrid2D(Nr, Nz, Lr, Lz)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32216205-b7a3-49eb-9fca-feb9cd98a637",
   "metadata": {},
   "source": [
    "#### Set up Poiseuille flow velocity field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "0264a95e-db00-4441-b988-9be0c8bf26f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "rr = msh.cellcenters.r\n",
    "zz = msh.facecenters.z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "399ea439-8a21-45ff-afe9-034f8e410164",
   "metadata": {},
   "outputs": [],
   "source": [
    "uu =  umax*(1 - (rr**2)/(Lr**2)) # does not depend on zz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0f0dfe6b-95f0-4532-a6c9-8d63e863661b",
   "metadata": {},
   "outputs": [],
   "source": [
    "u =  pf.FaceVariable(msh, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e141e300-0b88-4fb8-a3d4-251e9fa98f4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "u.rvalue[:] = 0\n",
    "u.zvalue[:] = uu[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "149e0f7b-f14e-4507-98b1-729df06739bc",
   "metadata": {},
   "source": [
    "#### Solution variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "efb9f15b-93ef-496f-8941-e507faaed92e",
   "metadata": {},
   "outputs": [],
   "source": [
    "bc = pf.BoundaryConditions(msh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "2563b8d6-a57b-4ed5-a958-5ffbc478fc66",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi = pf.CellVariable(msh, 0.0 , bc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae0736a3-107c-4f6a-b093-b5d06cb8d54b",
   "metadata": {},
   "source": [
    "#### Initial condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "89b7fda9-871d-4ad3-afff-2723f108ff79",
   "metadata": {},
   "outputs": [],
   "source": [
    "t=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "bc57e0e8-7fa2-45fb-b422-9b6a38e98c57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial condition\n",
    "for i in range(loadix0, loadix1):\n",
    "    phi.value[:, i] = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5c596a4c-c983-48cb-888f-5338c66b106b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 0.0 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAACkCAYAAACTpfSdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUUUlEQVR4nO3df0yV593H8c+Rc+DMPoVZsAd/IIUtT8GxZvWwVeQhbk8sTPeri0ttu/7Yps1YzCM/YlaVJrUuk66znWGKrBTabYnUVLvUpKTCnkxCC12DHpetkrqtKDwtpwaWgraZ/LqePyzHnZ6j4+B9Kr14v5L7j3Od731x3d9cCR/vcx90GWOMAAAALDbnWi8AAAAg3gg8AADAegQeAABgPQIPAACwHoEHAABYj8ADAACsR+ABAADWc1/rBThlYmJC77zzjq6//nq5XK5rvRwAADAFxhidO3dOCxcu1Jw58bsPY03geeedd5SRkXGtlwEAAKahr69Pixcvjtv81gSe66+/XpL0X1ojtzxXPZ/bd6Oeb//fq55n0v/83236v//+wLH5AACwwZhG9YqaQ7/H48WawDP5MZZbHrldDgSeOYlKvt65W2uJ/5Eot2vUsfkAALDCh//BVbwfR+GhZQAAYD0CDwAAsB6BBwAAWI/AAwAArEfgAQAA1iPwAAAA6xF4AACA9Qg8AADAegQeAABgPQIPAACwHoEHAABYj8ADAACsR+ABAADWI/AAAADrEXgAAID1CDwAAMB6BB4AAGC9aQWe2tpaZWVlyev1yu/3q729/Yr1bW1t8vv98nq9ys7OVl1dXUTNe++9p40bN2rBggXyer3Kzc1Vc3PzdJYHAAAQJubAc+DAAZWXl6uqqkqBQEBFRUVavXq1ent7o9b39PRozZo1KioqUiAQ0LZt27Rp0yYdOnQoVDMyMqLbb79dp0+f1sGDB/Xmm2+qvr5eixYtmv6VAQAAfMgd6wlPPvmk1q9frw0bNkiSdu/erSNHjmjfvn2qrq6OqK+rq9OSJUu0e/duSVJubq66urq0a9curV27VpLU2Niof/zjH+ro6JDH45EkZWZmTveaAAAAwsR0h2dkZETHjh1TcXFx2HhxcbE6OjqintPZ2RlRX1JSoq6uLo2OjkqSDh8+rIKCAm3cuFE+n095eXnauXOnxsfHL7uWCxcuaHh4OOwAAACIJqbAMzAwoPHxcfl8vrBxn8+nYDAY9ZxgMBi1fmxsTAMDA5Kkt956SwcPHtT4+Liam5v18MMP64knntBPf/rTy66lurpaKSkpoSMjIyOWSwEAALPItB5adrlcYa+NMRFj/67+X8cnJiZ044036qmnnpLf79ddd92lqqoq7du377Jzbt26VUNDQ6Gjr69vOpcCAABmgZie4UlLS1NCQkLE3ZyzZ89G3MWZlJ6eHrXe7XYrNTVVkrRgwQJ5PB4lJCSEanJzcxUMBjUyMqLExMSIeZOSkpSUlBTL8gEAwCwV0x2exMRE+f1+tba2ho23trZqxYoVUc8pKCiIqG9paVF+fn7oAeXCwkL97W9/08TERKjm1KlTWrBgQdSwAwAAEIuYP9KqrKzU008/rcbGRnV3d6uiokK9vb0qLS2VdPGjpvvvvz9UX1paqjNnzqiyslLd3d1qbGxUQ0ODNm/eHKr50Y9+pMHBQZWVlenUqVN66aWXtHPnTm3cuNGBSwQAALNdzF9LX7dunQYHB7Vjxw719/crLy9Pzc3Noa+R9/f3h/1NnqysLDU3N6uiokJ79+7VwoULVVNTE/pKuiRlZGSopaVFFRUVuuWWW7Ro0SKVlZXpoYcecuASAQDAbOcyk08Qf8INDw8rJSVFX9a35HZ5rno+d7pPLx0/4sDKLnqwr1C9t73v2HwAANhgzIzqqF7U0NCQkpOT4/Zz+L+0AACA9Qg8AADAegQeAABgPQIPAACwHoEHAABYj8ADAACsR+ABAADWI/AAAADrEXgAAID1CDwAAMB6BB4AAGA9Ag8AALAegQcAAFiPwAMAAKxH4AEAANYj8AAAAOsReAAAgPUIPAAAwHoEHgAAYD0CDwAAsB6BBwAAWI/AAwAArDetwFNbW6usrCx5vV75/X61t7dfsb6trU1+v19er1fZ2dmqq6u7bO1zzz0nl8ulO+64YzpLAwAAiBBz4Dlw4IDKy8tVVVWlQCCgoqIirV69Wr29vVHre3p6tGbNGhUVFSkQCGjbtm3atGmTDh06FFF75swZbd68WUVFRbFfCQAAwGXEHHiefPJJrV+/Xhs2bFBubq52796tjIwM7du3L2p9XV2dlixZot27dys3N1cbNmzQD37wA+3atSusbnx8XN/97nf16KOPKjs7e3pXAwAAEEVMgWdkZETHjh1TcXFx2HhxcbE6OjqintPZ2RlRX1JSoq6uLo2OjobGduzYofnz52v9+vVTWsuFCxc0PDwcdgAAAEQTU+AZGBjQ+Pi4fD5f2LjP51MwGIx6TjAYjFo/NjamgYEBSdKrr76qhoYG1dfXT3kt1dXVSklJCR0ZGRmxXAoAAJhFpvXQssvlCnttjIkY+3f1k+Pnzp3Tvffeq/r6eqWlpU15DVu3btXQ0FDo6Ovri+EKAADAbOKOpTgtLU0JCQkRd3POnj0bcRdnUnp6etR6t9ut1NRUvfHGGzp9+rS+8Y1vhN6fmJi4uDi3W2+++aY+85nPRMyblJSkpKSkWJYPAABmqZju8CQmJsrv96u1tTVsvLW1VStWrIh6TkFBQUR9S0uL8vPz5fF4lJOToz//+c86ceJE6PjmN7+pr3zlKzpx4gQfVQEAgKsW0x0eSaqsrNR9992n/Px8FRQU6KmnnlJvb69KS0slXfyo6e2339ZvfvMbSVJpaan27NmjyspKPfjgg+rs7FRDQ4OampokSV6vV3l5eWE/49Of/rQkRYwDAABMR8yBZ926dRocHNSOHTvU39+vvLw8NTc3KzMzU5LU398f9jd5srKy1NzcrIqKCu3du1cLFy5UTU2N1q5d69xVAAAAXIHLTD5B/Ak3PDyslJQUfVnfktvluer53Ok+vXT8iAMru+jBvkL13va+Y/MBAGCDMTOqo3pRQ0NDSk5OjtvP4f/SAgAA1iPwAAAA6xF4AACA9Qg8AADAegQeAABgPQIPAACwHoEHAABYj8ADAACsR+ABAADWI/AAAADrEXgAAID1CDwAAMB6BB4AAGA9Ag8AALAegQcAAFiPwAMAAKxH4AEAANYj8AAAAOsReAAAgPUIPAAAwHoEHgAAYD0CDwAAsN60Ak9tba2ysrLk9Xrl9/vV3t5+xfq2tjb5/X55vV5lZ2errq4u7P36+noVFRVp3rx5mjdvnlatWqXXX399OksDAACIEHPgOXDggMrLy1VVVaVAIKCioiKtXr1avb29Uet7enq0Zs0aFRUVKRAIaNu2bdq0aZMOHToUqjl69Kjuvvtu/eEPf1BnZ6eWLFmi4uJivf3229O/MgAAgA+5jDEmlhNuu+02LVu2TPv27QuN5ebm6o477lB1dXVE/UMPPaTDhw+ru7s7NFZaWqo//elP6uzsjPozxsfHNW/ePO3Zs0f333//lNY1PDyslJQUfVnfktvlieWSonKn+/TS8SNXPc+kB/sK1Xvb+47NBwCADcbMqI7qRQ0NDSk5OTluPyemOzwjIyM6duyYiouLw8aLi4vV0dER9ZzOzs6I+pKSEnV1dWl0dDTqOR988IFGR0d1ww03XHYtFy5c0PDwcNgBAAAQTUyBZ2BgQOPj4/L5fGHjPp9PwWAw6jnBYDBq/djYmAYGBqKes2XLFi1atEirVq267Fqqq6uVkpISOjIyMmK5FAAAMItM66Fll8sV9toYEzH27+qjjUvS448/rqamJr3wwgvyer2XnXPr1q0aGhoKHX19fbFcAgAAmEXcsRSnpaUpISEh4m7O2bNnI+7iTEpPT49a73a7lZqaGja+a9cu7dy5U7///e91yy23XHEtSUlJSkpKimX5AABglorpDk9iYqL8fr9aW1vDxltbW7VixYqo5xQUFETUt7S0KD8/Xx7PpYeLf/7zn+snP/mJXn75ZeXn58eyLAAAgCuK+SOtyspKPf3002psbFR3d7cqKirU29ur0tJSSRc/avrXb1aVlpbqzJkzqqysVHd3txobG9XQ0KDNmzeHah5//HE9/PDDamxs1E033aRgMKhgMKjz5887cIkAAGC2i+kjLUlat26dBgcHtWPHDvX39ysvL0/Nzc3KzMyUJPX394f9TZ6srCw1NzeroqJCe/fu1cKFC1VTU6O1a9eGamprazUyMqLvfOc7YT/rkUce0fbt26d5aQAAABfF/Hd4Zir+Dg8AAJ88M/Lv8AAAAHwSEXgAAID1CDwAAMB6BB4AAGA9Ag8AALAegQcAAFiPwAMAAKxH4AEAANYj8AAAAOsReAAAgPUIPAAAwHoEHgAAYD0CDwAAsB6BBwAAWI/AAwAArEfgAQAA1iPwAAAA6xF4AACA9Qg8AADAegQeAABgPQIPAACwHoEHAABYj8ADAACs577WC3CKMUaSNKZRyTgw4cSIhs9NODDRRSPnRzRmRh2bDwAAG4zp4u/Gyd/j8WJN4BkcHJQkvaJmZyZ8V5r3n85MddFbTk4GAIBVBgcHlZKSErf5rQk8N9xwgySpt7c3rg2z3fDwsDIyMtTX16fk5ORrvZxPNHrpHHrpDProHHrpnKGhIS1ZsiT0ezxerAk8c+ZcfBwpJSWFzeeA5ORk+ugQeukceukM+ugceumcyd/jcZs/rrMDAADMAAQeAABgPWsCT1JSkh555BElJSVd66V8otFH59BL59BLZ9BH59BL53xcvXSZeH8PDAAA4Bqz5g4PAADA5RB4AACA9Qg8AADAegQeAABgPQIPAACw3owNPLW1tcrKypLX65Xf71d7e/sV69va2uT3++X1epWdna26urqImkOHDmnp0qVKSkrS0qVL9bvf/S5ey59RnO7ls88+K5fLFXH885//jOdlXHOx9LG/v1/33HOPbr75Zs2ZM0fl5eVR69iTzvRytu5JKbZevvDCC7r99ts1f/58JScnq6CgQEeOHImom4370uk+sien1stXXnlFhYWFSk1N1ac+9Snl5OToF7/4RUSdI3vSzEDPPfec8Xg8pr6+3pw8edKUlZWZ6667zpw5cyZq/VtvvWXmzp1rysrKzMmTJ019fb3xeDzm4MGDoZqOjg6TkJBgdu7cabq7u83OnTuN2+02r7322sd1WddEPHr5zDPPmOTkZNPf3x922CzWPvb09JhNmzaZX//61+YLX/iCKSsri6hhTzrXy9m4J42JvZdlZWXmZz/7mXn99dfNqVOnzNatW43H4zHHjx8P1czGfRmPPrInp9bL48ePm/3795u//OUvpqenx/z2t781c+fONb/61a9CNU7tyRkZeL70pS+Z0tLSsLGcnByzZcuWqPU//vGPTU5OTtjYD3/4Q7N8+fLQ6zvvvNN89atfDaspKSkxd911l0Ornpni0ctnnnnGpKSkOL7WmSzWPv6rlStXRv0lzZ685Gp7ORv3pDFX18tJS5cuNY8++mjo9Wzcl/HoI3vyklh7+e1vf9vce++9oddO7ckZ95HWyMiIjh07puLi4rDx4uJidXR0RD2ns7Mzor6kpERdXV0aHR29Ys3l5rRBvHopSefPn1dmZqYWL16sr3/96woEAs5fwAwxnT5OBXvykqvtpTS79qTkTC8nJiZ07ty5sP+lerbty3j1UWJPToqll4FAQB0dHVq5cmVozKk9OeMCz8DAgMbHx+Xz+cLGfT6fgsFg1HOCwWDU+rGxMQ0MDFyx5nJz2iBevczJydGzzz6rw4cPq6mpSV6vV4WFhfrrX/8anwu5xqbTx6lgT15ytdc92/ak5Ewvn3jiCb3//vu68847Q2OzbV/Gq4/syUum0svFixcrKSlJ+fn52rhxozZs2BB6z6k96Y6p+mPkcrnCXhtjIsb+Xf1Hx2Od0xZO93L58uVavnx56P3CwkItW7ZMv/zlL1VTU+PUsmeceOwf9uRFV3vds3VPStPvZVNTk7Zv364XX3xRN954oyNzfpI53Uf25CVT6WV7e7vOnz+v1157TVu2bNFnP/tZ3X333Vc150fNuMCTlpamhISEiOR29uzZiIQ3KT09PWq92+1WamrqFWsuN6cN4tXLj5ozZ46++MUvWvsvl+n0cSrYk5c4fd2270np6np54MABrV+/Xs8//7xWrVoV9t5s25fx6uNHsSev3MusrCxJ0uc//3m9++672r59eyjwOLUnZ9xHWomJifL7/WptbQ0bb21t1YoVK6KeU1BQEFHf0tKi/Px8eTyeK9Zcbk4bxKuXH2WM0YkTJ7RgwQJnFj7DTKePU8GevORqe/lRtu9Jafq9bGpq0ve+9z3t379fX/va1yLen237Ml59/Cj25NT3jzFGFy5cCL12bE/G9Ijzx2Tya20NDQ3m5MmTpry83Fx33XXm9OnTxhhjtmzZYu67775Q/eRXqSsqKszJkydNQ0NDxFepX331VZOQkGAee+wx093dbR577DHrv2ppTHx6uX37dvPyyy+bv//97yYQCJjvf//7xu12mz/+8Y8f+/V9XGLtozHGBAIBEwgEjN/vN/fcc48JBALmjTfeCL3PnnSul7NxTxoTey/3799v3G632bt3b9hXpd97771QzWzcl/HoI3tyar3cs2ePOXz4sDl16pQ5deqUaWxsNMnJyaaqqipU49SenJGBxxhj9u7dazIzM01iYqJZtmyZaWtrC733wAMPmJUrV4bVHz161Nx6660mMTHR3HTTTWbfvn0Rcz7//PPm5ptvNh6Px+Tk5JhDhw7F+zJmBKd7WV5ebpYsWWISExPN/PnzTXFxseno6Pg4LuWairWPkiKOzMzMsBr2pDO9nK170pjYerly5cqovXzggQfC5pyN+9LpPrInp9bLmpoa87nPfc7MnTvXJCcnm1tvvdXU1taa8fHxsDmd2JMuYz58IhUAAMBSM+4ZHgAAAKcReAAAgPUIPAAAwHoEHgAAYD0CDwAAsB6BBwAAWI/AAwAArEfgAQAA1iPwAAAA6xF4AACA9Qg8AADAev8P9qIg47vhAD4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phi_visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "73f8ace3-dd40-4b75-9410-6ef40c8cda29",
   "metadata": {},
   "outputs": [],
   "source": [
    "initInt = phi.domainIntegral()\n",
    "# print(initInt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "f27415f5-5a94-41c8-9f7e-a35e57f61ccf",
   "metadata": {},
   "outputs": [],
   "source": [
    "phiprofs = []\n",
    "phiprofs.append((t, integral_dr(phi)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b32d38-17c1-4e69-8657-f1ec723d8fea",
   "metadata": {},
   "source": [
    "### Solve the convection PDE with time-stepping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "98594df8-a52a-416d-820a-5ee83c3a514c",
   "metadata": {},
   "outputs": [],
   "source": [
    "exect0 = time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "6ca90110-c48c-46c2-83c1-c6ace789e13f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The convection term only needs to be constructed once, since it\n",
    "# will be constant during all time steps.\n",
    "# convectionterm = pf.convectionTerm(u) # really ugly results?\n",
    "convectionterm = pf.convectionUpwindTerm(u) # numerical diffusion\n",
    "\n",
    "def step_solver(Nstp):\n",
    "    global t\n",
    "\n",
    "    for i in range(Nstp):\n",
    "        # Transient term needs to be re-evaluated at each time step\n",
    "        transientterm = pf.transientTerm(phi, deltat, 1.0)\n",
    "        eqnterms = [transientterm,\n",
    "                    convectionterm]\n",
    "        # solve the PDE using oneMKL PARDISO as the external solver,\n",
    "        #    via `spsolve_oneMKL_pardiso`\n",
    "        pf.solvePDE(phi, eqnterms,\n",
    "                    externalsolver = spsolve_oneMKL_pardiso)\n",
    "        t += deltat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "50f94ce2-cd02-4849-a4b0-c287ee87205e",
   "metadata": {},
   "outputs": [],
   "source": [
    "step_solver(200)\n",
    "phiprofs.append((t, integral_dr(phi)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "ad8521f0-72d0-4a1b-b529-34ac0563aa26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# print(t, initInt, pf.domainInt(phi))\n",
    "# test conservation of mass\n",
    "assert np.isclose(initInt, phi.domainIntegral())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "ec265066-a81f-4e47-a6d6-d63c50c0e622",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 2.0 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAACkCAYAAACTpfSdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8cElEQVR4nO29fbBdZXn//blf1tr7nJNzTkgIASSJwU4FS/1VkrZGyvD7Q0O17aiPHa1vtR1kmspUIONUqE4FOgNWraUIAaWgtQ7I+NIRZzJK+rQylMT6gNFphSlUkKAkYEKS87r3Wuu+7+ePe6219z4vISeck5ycXB9mzd577bXXWfvOmtkX1/f6XpcKIQQEQRAEQRCWMPpEX4AgCIIgCMJCIwGPIAiCIAhLHgl4BEEQBEFY8kjAIwiCIAjCkkcCHkEQBEEQljwS8AiCIAiCsOSRgEcQBEEQhCWPPdEXMF9473nuuecYHBxEKXWiL0cQBEEQhKMghMDo6Chnn302Wi9cHmbJBDzPPfcca9asOdGXIQiCIAjCMfDss89yzjnnLNj5l0zAMzg4CMDv8BYsycs+n7KWyTf+H37++553X/gD3j60mzVW0adSjJp7BPpsMcaVT/0/PPevazlz5zj2qX34wyOEIgdpdi0IgiCcohTk/Afb69/xhWLJBDyVjGVJsOrlBzwEQ/+4pv/AMna8+BuE/n5+q/EUr0p+yemmxXJtaagEjTqqAGhFoTnvzFGe+pWE0dEVDA8MkD4/ijo4gh+fIGQZIS8geAmABEEQhFOH8idvoctRlkzAM+8ET/rTfbwiOZsXn1/J13/19ex41av5nbOf4neGnuDX0n2sNi0GdQqBlwx6VpgG71z5n+hNnv/37F/l8JNDDP1vg+Gnh2j8/DDqxcOE8fEY+DgnQY8gCIIgzCMS8MxGCLhf7if1ntOzsyH0cShdzk/6z2JN80VeYQ8yrCdoBodWGvMSp2uohAvTFs2VO+kzOd8OFzDeGiIdS0hG+jATLWi3Y5YHdzy+oSAIgiCcMkjAcwSCD1AUmPGM5qEmjV8anjltBbvScwEYbf6cdfYgK02LfmVeUuLSaAZUwTnpQVYPj/KzM/oZP2wxWT/9Cmxi0aPjhIkJQpaXmR6RuARBEATh5SIBz5EIHj82jtl3gGHAtpZx+FA/P37hVTy5bhX/Z/UaNg7/jF9vPssaO8IK7ehXCQlmxsCnoSxnm4I39D+JfoXn35uv5sfLz2HyjD6WnbGMwZ/30bdvAnNglDAySphsicQlCIIgCPOABDxHIgRCluEPHUaFwEAIwCDBWMaSQZ5IV7E8mWSFGWNQt+hXbRLlMChmErmM0izTDc61bej7KawEHxQ/9ucwkfVh2gbTatCczFFZDnkRgx0fAMn0CIIgCMKxIgHPURCcRxUFqpVjxxzpYUNyUHNgaBn/0zyDAdvGqIBPXmCVmWRY+1kzPRqFVop+VbDKjnBm3yhPD05yeEVKa8xgWwm6GCAlVqzrcUMoa3tE4hIEQRCEY0MCnpciBIJzhPEJlDY0AV0MYSdTxsb6+Onhs9l79hA/XXU6vzH8c85vPscrk/2sMhmDStNQdlrg01SWVSbj19J9MAwDps3/17eOPYMryYZTWqc1WbY3oe+FPuz+MfTIGGFyktBqE4oy8CmvTRAEQRCEl0YCnqPBO3wGenQUBaTAMjVE0A2CNozbAZ6yKxlM2iwzLQZ0m6Y6RKIdifL4LonLKA2BaGcnw6XPw7IobeVe81yxAp2nmExjsgTdbmLyAgoHrszu+BAzPYIgCIIgHBUS8MyBEEKUtrIcM1mQTCQkYwZ72HB4WT8/61vBoG3RVDkaD3YER8GgDiSEOtNTYVA0lWNIt1iZjLOqb5z9g8vIhi12XGMnDbrdIHUhhktaoZQiKBUzPVXgI5keQRAEQTgiEvAcLd4RcgiTk6AVVimWuYDJ+jAtw8REHz8bW82BM/v5+YrlvGZoBb/a3Msr0/2cacYZVI5+bWgqiy7/ayjLCp3j7EHoh0QX9NuMnzTOZKR/kGKZJRtMGRiyNH+ZYg820aPjMD4B7QyyDJwTF5cgCIIgvAQS8MyF4AlFAZMtUAoDNI0imCbBGLw1jNoBnjGOfpvTrzMGdEZT5SS6TRI8CQFdJnkSZegHVpkcOIhrKnzQFEHzeFCM+WWowqILjc4TVOFReYGq5C3nYpZHiYtLEARBEI6EBDxzIYQYYDiHch7yAt0qMJOeZNyQjCmKUcvoQB97+4YYSloMmhZNnWEIQAY6p4knUR3bukGRKs+Ayhg2E5zeGGeor8X4QJN8mSEfVLQnNTpL0HkT7T3K99bwBAcS9AiCIAjCzEjAM1dCmVmZnARAa0UzBHTRj84TTFsz2ernqYmEg6v6eOG0ZTw/OMSvNJ9nTXKAM80Yw7pgQHkSFaWtRGmGdcDZ0XhOFUh1wU+SnGfSFYz29VH0aYr+hKLf0BxIsAdT9OgkTEyi6gaFXqzrgiAIgjADEvDMldKmDkC7DWOxr06qFMEogrEEowkm4ZAd4OfW0TQFTZ3TVHFL1CQGj1YKTRw50VSwXDuwozgUPigKb8id4ede03INlNeoQqOLBJX7KHEVXTU8IUMFRfAayfYIgiAIQgcJeI6R4EOUlZyDwqFyh5l02JbBTgTshKKYsBwe7+OXjWUsT4ZZZlo0dR67Mes2eEe/jpIWgAYMgabKGTQthuwkyxuT7O8bYLw/oehX5MsUWUtj2gk69+jCoUIo7eqeACicBD2CIAiC0IUEPMdCiEXCoSigBahoNk/KwEO5Riw0ziyt1gD/27aMZA1eHOrn4LIBDjX6Odse5AwzxnDIGdCKBIVB0a8Cq8wk8MvSueVITcH/2tM5mA7iGgmuqSmaCX19mkbTYg8n6FGLShJCq0XI8hiIVQ6u+poFQRAE4dREAp5jZaq0VfbIsVpH55ZOCVoTjKZlUg6YZSTa0zAFDVVEeUvn0bkVHFpBojRNysyMmcRxABcUuTcUXlN4zahbRttZlFPowqCLGGApV2WcPMqX2R6lo4NLmhQKgiAIpzj6pQ+ZzrZt21i/fj3NZpMNGzbw0EMPHfH4Bx98kA0bNtBsNjn33HO54447ph1z6NAhrrjiCs466yyazSbnn38+27dvP5bLO/7U0laBygt022FbHtMibpOafCLh8GSTA+0BXsgH+WUxxItuGaM+ZSIoWiGQB48jZmIMgRTHgM4YtpMsS9oMNtvovgLXFyj6oOiHvF/j+iy+LyE0UlQjhTSBJEEZjTIGlAalXuJLCIIgCMLSZc4Znvvuu4+rrrqKbdu2cdFFF/H5z3+eN7/5zTz22GOsXbt22vFPP/00b3nLW7j88sv5yle+wsMPP8yHPvQhVq1axTve8Q4AsizjTW96E2eccQZf//rXOeecc3j22WcZHBx8+d9wIenO8mQ6Sls6DpFoAARQ3qILRStPGMkHeSKzjGYNRgebHO7vYyRt8orkICv1BIM6p6GIn1cwrNs4ewhX1vgkypNoz3N2iFbah08NrmFwDUWzqUkbFnM4QY9blG0RquaEher07BEHlyAIgnAKMueA57Of/SyXXXYZH/zgBwG4+eab+e53v8vtt9/OTTfdNO34O+64g7Vr13LzzTcDcP755/PII4/wmc98pg547r77bl588UV27txJkiQArFu37li/0/Gl6s2TF/UupRRGK1KtCFoRtImPxtDWTX5pPal2Ud7SOalyJIkjCR6Nw5RBz4D2QAvsQXzQ5MGQB03uNc97TeYaKGdQTqMc5RZ79CgfAxtVBjjBB5T2sV+PUhL0CIIgCKcUc5K0sizj0UcfZfPmzT37N2/ezM6dO2f8zK5du6Ydf+mll/LII4+Q5zkA999/P5s2beKKK65g9erVXHDBBdx44424KnsyA+12m5GRkZ7thBF83Lyv+/RQOHTmMJnHtgOmDaalUG1Nq5XEIuasn4P5AIdcP4dcP6M+pRU0OVB9c6NCdG7pnGWmxZBtsyzJSNOC0PC4ZilvNRWuqXBNi28khDRBJQlYC8agjC6lrfKfXCQuQRAE4RRiThme/fv345xj9erVPftXr17Nvn37ZvzMvn37Zjy+KAr279/PWWedxVNPPcW//du/8d73vpft27fz5JNPcsUVV1AUBX/9138943lvuukmrr/++rlc/oISOzB7yPLataWVIgWgCSEWGitnaBdNns817dwyUaSMDTSYaDaYSBqcaQ+zPEzSrwsSYrAzqDMcI7hE4/tiwKJU4BfGM2b6CTYhWIVPDC7VNFJN0jCYEYuyBoyBzE6Xt8S2LgiCIJwiHJNLS03JDoQQpu17qeO793vvOeOMM/jCF76AMYYNGzbw3HPP8elPf3rWgOfaa69l69at9euRkRHWrFlzLF/n5VPZ1KckpJTWaKVItCIoCNrG5oRa09Yph03A6IBVjqRr08qjQ8CoAg00lYtBjx3BoXAo8qBxQfOcV0x4hfIW5RXKU8pcCcoFdAioStKq5C2kV48gCIJwajGngOf000/HGDMtm/PCCy9My+JUnHnmmTMeb61l5cqVAJx11lkkSYIxnflS559/Pvv27SPLMtI0nXbeRqNBo9GYy+UfH6ZKW86hco/OAyYLmHbAZArTVhQtw2SWMJo3OZT3cziZYNBMMuDbpDiMDiQqWsoNAYMnVY6mKlhm2gzYjEZSMNlwsXg5VbhGtWlMalCphdyi0oTgHZhY3xOUlpoeQRAE4ZRhTgFPmqZs2LCBHTt28Pa3v73ev2PHDt761rfO+JlNmzbx7W9/u2ffAw88wMaNG+sC5Ysuuoh77rkH7z1aR8nmiSee4Kyzzpox2FmUVFker3ulLaXKLA9AE7AQNMor2j5h3CmeLTTtwjLpEib6U8YbUdpaacYY1C2aKqaOBlSOM2O4pLf0SqnAIT1ASzcIVuOtxluFTxRporFWo41BaQ02g7aGvCAUBQo68lb9PQRBEARhaTFnSWvr1q28//3vZ+PGjWzatIkvfOEL7Nmzhy1btgBRavrFL37Bl7/8ZQC2bNnCrbfeytatW7n88svZtWsXd911F/fee299zj//8z/nc5/7HFdeeSV/8Rd/wZNPPsmNN97Ihz/84Xn6mseJbmlLq1gzoxVqopS2lCLoUtrSGm8UQVvausl+HTDaY7Uj0ZW8VZCoAkPM9CTKM0CGN2OlrGXIg6EIGu81h70i8ynKabRTcfaWs7EpoYvNCFVppVc+TJe3pEGhIAiCsESZc8Dzrne9iwMHDnDDDTewd+9eLrjgArZv317byPfu3cuePXvq49evX8/27du5+uqrue222zj77LO55ZZbaks6wJo1a3jggQe4+uqree1rX8srXvEKrrzySj760Y/Ow1c8QZQBRSVr4eOwT10ETB7QOZg26AxUpskzy2SeMJY3OFz0MWwmmdANmirH6AAUaBWzL1pFaauhc/pMTr/NaKY5Y4nHpQHXCFHWSjvylm4bTJ7EuV+FJbjStl51ZK58YSJvCYIgCEsQFcLS+HUbGRlheHiY/8tbsSo5sRejVJSzEouyFpWmqP4+Qn8Tv6yPYrhBttzSWm5on6ZonxbITvOo0zKWD49z1uAoawcOck7jIKuTw6yyIyw3EzH4IZAHw3hIOeCW8Xy+nD3tFTwzsYJnR5dz4NAy8kMNkoOGxiFF42CgcTiQjjiSkQwz1kZNtFGTbUK7DVlOKIrYQLG7OSFI4CMIgiAsOEXI+R7f4vDhwwwNDS3Y35FZWgtBKW3hqkER1LO2tFJYExsSeqPwppS2jCYzCSO6v3RueaxyGOWjrKU8aGiqOG19gAynJ/BW19JW5iyFMxz0mqKWtBR44nOfoHxAe1DltPdQDTyFTnNCkbcEQRCEJcYxzdISjo7gQ5S2ID52S1vOR9dWDjoHnSl0pnG5ppVbJouE8aLBhE+Z8A1aPiUPBh90PWrCKI9WnkQ5Grqg32Y0bIFNCkLi8Qm4FHwKLo1FzD4xhMQQrIlNCa0tmxKqWHekum4JaU4oCIIgLBEkw7OQBA9BEYpy7IRW0SmlFKariDkOklCAJgspo17hvSbzhra3tPoSWiFuy804AypDK48hMKRb5DYWMfug8GUwtB9oqyYoC6oacQFBKxITs0zaaNRkGeyYPM4DK7tfE5S4twRBEIQlgwQ8C0UIMdDwAYWLz7McVCltaV1LW6GUttCxKWFmLBOmwYs6ztzquLaic8voQJMcrTxNcpabcVxQ5GklbRkKr3nRaXKvUN7EpoSVzOUNKiTgAtrHxoSxgDkQfAyZgkPkLUEQBGHJIJLWcSKEUDYk7Li3VOHRVVPC0rlVyVsh02SltDXp4jbhU1o+IQumzOh0/vnSMiBq6IKmKWjaAps4sB5vAz4BnypconBp7NUTEg3l6AllLRiNMgZ0DL56EHlLEARBOImRDM9C0tWMUKkY5HRnedAaY6LEFLQlKEVQpbSFpRUULwRwQVH4sjg5tbSSNEpbuk2Cq4eLLjcT5ImhaGqKoCm85kBQtIBMWYLS8fxKd8ZcmHgNqiyqjh2YFeQqBjnO98pbIm0JgiAIJyES8Cw0U+ZsKe2iFVxplNFoo7FalzU2MQCpApFgoG0aHDKeRHusjgXKiS5Iy4aERrfReJoqZ1C3cEaTJzY2JPQaFxQHnCJ3paTlVN2UULtQNib0GBcLrFUIhOCjvFVlo7z06REEQRBObkTSOt7UVnAfsyfOQ+na0kWIWyVt5THTkueGVmFpFQmTPmXCNWj5JNbrBIMv/xmN8tHGrgsS5UiNo2EcSeIgCfgkEBKivGXBJQpvVena0tG5pXWcrl5u1YBX1e3gEnlLEARBOMmQDM/xIsyS5Sl782DKqerGdrI9SgGWHHgRahdW7k18DJZhM86AzkhUdII1Vc6wmezJ8hQ+jp5oBciCJca5Ub6qHFyJAlNOd4/BTbSph/J6g3O9Iygk0yMIgiCcREjAczwoHVvxuY9djZVC5VksFG5HacvoMugp3VpR2lIEY8hNyogOJMZhVXRuaUK0p5cjJxJV0NQ5Dk1uqzlbhsKX8lYpbVXyFl6hHeXsLVs3I6wcW7VzC3onrEvQIwiCIJxkiKR1PKns3dWcLR/iTKtyUz7Ugz6VK2WtAlShoNC4wpAVlpZLaDtL21vapbTlUHg0LmgMHkOs97HKkZqC1Dis9WADIQl4C8FW8lYpbRlNsBpM3FTp1qoKrKc1JhQEQRCEkwTJ8Bwvuvvy6K4sT1FAPiXLY7qyPGXGB2UoSDlUnq6Stxyx8/KwsdG1VUlbOmfQtGgnlqIcO+GCwntFO0DuLSpoCIo4/8JUvQ+hlNmUio0Sg2p3nFuu1OWkMaEgCIJwEiEBz/FkJpt6UdRZFKU1WmuM0THoMXS5tmIGprCWEdPEaB+zNuVoCUPAEOdtJaogwdGv2wwbQ24N7YalCJrcGQ45TV50Obe8QlXztpwFF2JTQudRoTNza0bnlpLGhIIgCMLiR/SJE0XwnWaEtcTlUd6jnEcVpaxVgHKguqUtp8kKQ7uwtMvxE9UAURc6DQnrWVva0TAFqS5IjMNU0paN0lZnK6Utq8GWTQm1LjdTSlqq17lVIc4tQRAEYREjGZ4TQZgly2MMqpWhbZnlsXH+lTdR2oqjJwxOJYxC/HyJDwqXaJzV5UT1jmtrmWnRsgnt1JJ5S+4M3isyryg8qKBjlifE59Vt0fFzEZsSls+D1ijy6eMnpIhZEARBWKRIwHO8merYCgrlHOSKoPMobU1qTOnSsrWcRY9rqzCBMdMkMZ7EOBIVGxMa5cGU2Z2ycLlfZwybSYpEx2yQMxQuZopcV0PCuFXSlkG5mHXSITYkrJsSUhC8mXnmlgQ9giAIwiJEJK0TRfCx6NeXsla5r5615T3KhbIbcpS3qk2X0pYvpa3cmdKxFaUtH6Jby3f985oyIEq0i9Z24zHGgwl4G2JAVUlbVXBVObe0JphqqrqJdTu6PPfUmVuCIAiCsAiRDM+JYIYsTxwmWs2yigGFNgZTZnmibbwMRMoRFIVOGFMBo0OPvAXgUDRVjlaxmLmh85jpSSajpb1pO1keD3nZl6cuYK6krQDdp1bQkbYA5Vz5umxKCIhzSxAEQVhsSMBzIplaywMopaNNvZyzhYlFxN3ztUIdABlymzBqPEZ3ZXBUbEqIgSYx6Elw9Jt2HEeRdhxbhdeMO1VKW6ZX2qr3BUzp1qrnbfmACr7TlNCYMuiRIEcQBEFYfIikdaLozn5Ujq3u52VjQipJy0+XtZQDCoV3mtzFACb3pnZs+aBxaHyI2SRDiNKWKu3sxmG1R5kAppS1DHVTwrpmyPb2BZrWhLCSt9TUR5G7BEEQhMWBZHhONFOzPFrHZoRZlLa0iRbxWE8TszzeTsnymIQxPV3a8iiWoWiWji2tPA1V0GdyBpMWLWdpN2IR80R3lsdXoyeAoFHBvLS0BQQcClN/JyliFgRBEBYLEvCcSI7k2FLRpo7RKKtLq7rC2k49T6eupwx6SsdWauKsrThHPWBMqJsTNnUeZ2wZTVbZ1H2s5ZksSteW1yin0E5RONBOl5mlgPEe5ZNeaYveoIeyGXPdlFCCHkEQBOEEc0yS1rZt21i/fj3NZpMNGzbw0EMPHfH4Bx98kA0bNtBsNjn33HO54447Zj32q1/9Kkop3va2tx3LpZ3cTK1/ca5rvpYvuyF35m3pWuJSBKdwpayVdUlbcap6nLVVERsSBqzy2HL4qNYBZQLBBIKm7P9TbZ2CabTuSFvlRPV63laF0r1NCQVBEAThBDPnDM99993HVVddxbZt27jooov4/Oc/z5vf/GYee+wx1q5dO+34p59+mre85S1cfvnlfOUrX+Hhhx/mQx/6EKtWreId73hHz7HPPPMMH/nIR7j44ouP/RudjFSyVs+MrdKxlRdgygxPKW11OiKXwYiFYDtZnqqAWZeBTfWvnChXPzZ0QZ/JWJa0aTlL5gxFoWkXGleUWR7fyfIor8GDCmUBcyinqJePVbfoOsvjZbK6IAiCsHiYc8Dz2c9+lssuu4wPfvCDANx8881897vf5fbbb+emm26advwdd9zB2rVrufnmmwE4//zzeeSRR/jMZz7TE/A453jve9/L9ddfz0MPPcShQ4eO7RudbFSyVvCAKQuVy1oeVdSOLTJTd2C2ZT2Pt8TnNhYUF8bQrqQt7cs5W3FyurYBNGUQ5KO0ZQxDSYvMx4xQ4XRZAK1whUIXUNSOLR0zSt50Zm35csq7SFuCIAjCImdOklaWZTz66KNs3ry5Z//mzZvZuXPnjJ/ZtWvXtOMvvfRSHnnkEfI8r/fdcMMNrFq1issuu+yorqXdbjMyMtKznez0WLrLZoS1Y8v5LrdWfNRdri0VSmnLl9KW19F6Xjq14pytzj+3LoeNarqkLR1QOoDucmzVHZ4ppS4FptyqhoSVawtE2hIEQRAWJXPK8Ozfvx/nHKtXr+7Zv3r1avbt2zfjZ/bt2zfj8UVRsH//fs466ywefvhh7rrrLn70ox8d9bXcdNNNXH/99XO5/MXLTMXLIWZ6cBqyPDYjzAw6MZi2xicaWw3/TEqZy2ictbRN6MnymC57VRPVkba0o6Fz+kzOQNKm7SxZYWIzwkLhCoMqFEU1wNTpso6ozPKEgC7t87W0lSHSliAIgrDoOCaXlprSXyWEMG3fSx1f7R8dHeV973sfd955J6effvpRX8O1117L1q1b69cjIyOsWbPmqD+/KOmyqAdmqOWxBtWO9Ty2rN3xtrSp26qWR+ONpW08o7qq54l9dwy+nLEV17+q5RmwbfKga8dW0e3YKhS63kIpcWl0YepsUyVrTWtIiEhbgiAIwuJgTgHP6aefjjFmWjbnhRdemJbFqTjzzDNnPN5ay8qVK/nJT37Cz372M/7gD/6gft+Xco61lv/5n//hVa961bTzNhoNGo3GXC7/5KGq56maD5YztpTzZSNC32lE6Mut7owcwCm817igcF5TeENRSlsuaLxSddCjVYjjJ7ofS2kraKJry6jasVXJWp0RFyrWGKku15bS1OMloJS2vHRhFgRBEE4Ycwp40jRlw4YN7Nixg7e//e31/h07dvDWt751xs9s2rSJb3/72z37HnjgATZu3EiSJJx33nn813/9V8/7H//4xxkdHeUf/uEfTv6szdHSLWtBDHRU7LocAx0Xt7xAZ4aQGHTmsa1SzkootziKwllDyySM6aqjcifLAx3HllFxonqqC/ptRttaWqmlKAztQuMLhStCnemp5a2uieozZnmYUsAs0pYgCIJwApmzpLV161be//73s3HjRjZt2sQXvvAF9uzZw5YtW4AoNf3iF7/gy1/+MgBbtmzh1ltvZevWrVx++eXs2rWLu+66i3vvvReAZrPJBRdc0PM3li9fDjBt/ylBt0W9akSoFKEoUEp1anlsx60VmxFS1vF0bOqFTZg0AWtiLU+PVV332tT7TE7hDVliomurMLhCkxfRsaWKmEGqRlpop9FFKIMej/Y2ZqG8j7O1gkcFTfAepRRBI9KWIAiCcMKYc8Dzrne9iwMHDnDDDTewd+9eLrjgArZv3866desA2Lt3L3v27KmPX79+Pdu3b+fqq6/mtttu4+yzz+aWW26Z1oNHYHqWp94dMyf1MSHU1vC6GWE15dwT1SQPwSu8VzivKIKm8Bof4mwtHxS+628ZfC1vWe3R2peOLQg6ylqoblmrlLbKye21pKXKWVvdc7aqAmyRtgRBEIQThAphafwv9sjICMPDw/xf3opVyYm+nGOnCkKURhkT62OMQSUWkhTVSAnNlNBs4PsTimUpxYAhW6bJlymyQUW+DIplgWLQoQYK+gbaDPW3GG60WN6YZDiZZJlp09BFzPYAbW8Zdw1G8iaHsj4Otfo4PNlkcqKBG7PocUMyprDjimQMkrFAOu6xEx474TATObqVoyYzVJZDlhPyHIqCkBext1AZ+NQBTyjrfJbGLSgIgiAcA0XI+R7f4vDhwwwNDS3Y35FZWouNqRZ1B0prgvMo46AoUHkMgnSmMW1HsAqbKHyiMEk56bxycBlD2yaMlzb11Li6IaFWgQRXZ3ZSHQeLZtaQJ6bsvmzwTYV3CpcbVAG63CrHlnIBVURpS7kkZqR8iL2EyucBUM5Fx1mV5RFpSxAEQThOHNMsLeE4U7rW8CFmScoxDrWLy3caEmoX4giIUuLCq1Laiq6tKGsp8lLaqjCUHZnLQEipgJkmbVE/diQtYsBSNSLscmupbudWhTQiFARBEE4AkuFZrFTFy6as4fE+ZkOcA+9RhYPcoXOHyTU+05gk4DKFaUfHlkkUPtF4a8iMZcKkMcujo2urakhodcdCnihPqguapqBtC9pJQZ4YfGLwDYUrNCpX6Ebsy1P16FGFRqUGXXiUt+A8WFsHZpLlEQRBEE4kEvAsRqbIWnjdNV/LR8eW1mA0KusMFvWJwrZDlz09ureKRFMkhpZNSKyrbeqJ8ljl0CF0FSw7GtrRtDmZN2SpIXcGV5jSoh4dW7p0bOkcdBEdW95F11bwHpyNjq3Spo7pykqF0GVTlyBHEARBWHhE0jqZqLI8PZIWXZJWVxPC0qnVPWOrdmz1yFq90haUA0aJAZACVNmMkFLSQkNQUdbqaUZYbp0GhL2SVtVxe9au3Efo1i0IgiAILwfJ8Cxmpvbk6e687FyUtZxD5QaVe3TuMZnCpArfDvhU4dsKnShc21Aknnae0LKOCZNidSxejuS1Y8uoQMMUpK6gYQsa1pInDpdqfK7wmcGn4HLQObi8S9ZyBlX4cpK6jc0InSNYG6/bmPK7henNCAVBEARhgZCAZ7HyEp2Xg6nqeApU5dhKNMGqWLtTdV+uhosmCt+OtTwtk8TBorpqSOjQwZDgMESZyypD0xRkNqdtLVliKBJDkXpcQ0dZK48SV3ctj3e6y7EVYuPBMvAJ3sTOy9KMUBAEQTjOiKR1stEta1E2HSyLgpUPEKpGhKWsFSqZS8UZW0HFGVtV80FUKWv13gqVW6tuRKgCSvtS0ipdWz1NCEtZK2pgsTu0mkHegtiMUJd/T2mUOLcEQRCEBUYyPIud2WQtX87W8gYKF7MquUdnHpMqXB7QOZgsZnd0AiHV+MyQWUvLOFKTRMeWKu3owfTIWrp0bCXG0UgKssLgEhedX0lZHJ0qfKZwCegEdKLxaSxe1mURM86D8b19eVS0rQe6sjmS5REEQRAWCAl4FjOzyVrOobQiFC42I7QOlTt0bjC5ry3qvt09VJT42mqKxNK2npZ1pDqJjq3Spm5LjUkTaGhHoR0NU9A2ltQ68sTh09KxlUeXlmtQNiSsZm1pdKEJhSaUAc+cHVsS9AiCIAjziEhaJztVV+NQdjYupazuuVqquxFhUARfjrfyupS0OjLXTOiuJoRaB5QOhK5mhN1NCGMjwvi6R8qa4sCa1bGl5JYUBEEQ5h/5dTkZCGX2I8SiZabJWj5mUYq46dyXhcQBU0pbsV+OQpUuqzyz5IWhXVgyH7fCG4pg6noeXRY0W+Wio8t4rPXoxIMNUdZKQ1kYXW5WRZkr0XirCYkBawg2PsbZYDoGQ8aU9Txq5loesakLgiAI84RIWicb3bKWMXHGVuFQugBr0LkjFLp0bQVcotBtMKWkpa3CJZoiN7RzWzYiTKJbq5qvpR2agFEhdl42joavZK2C3BpcomNNUK7wKfiCWuLSRbTC60TjS5u68h4KN737clXLM5NjSxAEQRDmCcnwLHaOVMcyddp4lf1x3Y4tykLnjrRF6d6qGhF23FpqWhNCoG5CaLSvmxBWslbcps7Zig6t2ISQjlur3GZvPCiOLUEQBGFhkIDnZCF0inqnyVre98paLtSSli5KSauIhcWqAJVrQq5xRRwbkVWbN+RBk3uDJwYetTW9S9YyxqNtLD4OlihtWTrSlqWWtYLRYDXB6ihrGV1KWaZL0tKz1/KIrCUIgiDMAyJpnWxUs7WqJoS+bO5XlG6tYgZZK1foDEwKIYvOLZdqXK7JsmhRT3SUrqz30ZLeJSkl2mG1IdGO1Lha1vKpxhcKlRtCQpS2coVLY4ClrcKnunRtGZT1ZSNCjzKaEDQq6Dh7C2Itjzi2BEEQhAVAMjwnM1Nma6nuGVsu9DQgVD7UUlYlbQWvCEHhy9la1YwtYIZGhJ0Bo/VsLVXO11JTpC0Vt2DoNCLUdDUenOLa6m5EKAiCIAgLgPzKnAx0ZTd6ZC3ozNbqlrXKYEcXoUveonZrKQcUCgqNc5rCa3IfpayifIz1PJ3bQxOwymOmylq6lLVsOUi0em5Lx1Y57iIYTTCqR8bCmFi3M1ttj1jUBUEQhHlCflFOJuoi5U6xcogNdaIs5OJQUWprukdnU63psfuyzhRUFvXSsdUuLC2X1LU8VbanKlquGhTGWh6HNg4Sj089PiHKWmX3ZVfO8wpW1fb02qJuuup2Kot6FQSJRV0QBEFYACTgOVnxZaPBbqqZWnVRc1SdOjLWFFkrqLJfYSlt9Wyzz9bqlrWUDqWE1SVrqSkztkxnxtb0+Vr6yJkcyfIIgiAI84D8mpzsVHU8znWNavCdGp7Co10oty63Vpes5b0md1HOckFTVJIWvTb1apCo0R5bylrKdLm1umStYOLmrcIbFV1aWoFRvZJW2ZG5R9aaLcsjCIIgCMeIBDwnC3WvnVns6SF0ui4XDuXK+VouoPPeGp74PDYJpFD4XFMUmqzuvFxa1Lvs6UaFso7HkWpHYhymDHqwgZDEoKeWtZKpspaKspbVHVnLaJTWnYzPbMXLYlEXBEEQXibHFPBs27aN9evX02w22bBhAw899NARj3/wwQfZsGEDzWaTc889lzvuuKPn/TvvvJOLL76Y0047jdNOO403vvGN/OAHPziWSzt16LKNVwXMoQp+oJazKhmrakTYkbOq+VrRrTWzrDVzgKFUQBGDIFTZhFCFKGWpjqRVu7U0oGeQteoT9spaszYmFARBEIRjZM4Bz3333cdVV13Fxz72MXbv3s3FF1/Mm9/8Zvbs2TPj8U8//TRvectbuPjii9m9ezd/9Vd/xYc//GG+8Y1v1Md873vf493vfjf//u//zq5du1i7di2bN2/mF7/4xbF/s1OB7joe38n0qHI6eY+05YjZniI+V0U5XNQBLlrTi3JzQceZWsSBolObEFbSltYeY6JTCw3BBIINXXJW3IKJXZePVtYq/5gULwuCIAjzhgphbh3dfvu3f5sLL7yQ22+/vd53/vnn87a3vY2bbrpp2vEf/ehHuf/++3n88cfrfVu2bOHHP/4xu3btmvFvOOc47bTTuPXWW/njP/7jo7qukZERhoeH+b+8FauSuXylk4fqx74KBJQuh3EaVGLBWlSSQCMlNFJCM8H1p7g+SzFgyAc0eb8iH1AUA5APBFx/wA14VH9Boy+nv9lmIM1ZlrbptxmpdljtMCqQl1JXyyWM5ynjWcpYq0GrleAmLGrSYMY1dkJhJ8BOQDIRSCYCdsJjJx1mskC3C3SrgHaGygvIckKeQ1EQSqdZcK623NdNCKeO0hAEQRBOeoqQ8z2+xeHDhxkaGlqwvzOnDE+WZTz66KNs3ry5Z//mzZvZuXPnjJ/ZtWvXtOMvvfRSHnnkEfI8n/EzExMT5HnOihUrZr2WdrvNyMhIz7bkmeGHvqcfT/dx5X4VomurakJI6HZuqVrmopS1qiaEM0latVOLbqdWmf/pkrDQndeVxEX92JXFEVlLEARBOE7MKeDZv38/zjlWr17ds3/16tXs27dvxs/s27dvxuOLomD//v0zfuaaa67hFa94BW984xtnvZabbrqJ4eHheluzZs1cvsrJTZgyemHq8NC6iDnUTQirOh7tiPsc5aZQruy67KOEFTeND7p0amnc1OCnGiiqAtr4WMdj4uZNiG6trs1bFaUtW8pb3bJW1X9HZC1BEARhgTimouWp//cdQjji/5HPdPxM+wE+9alPce+99/LNb36TZrM56zmvvfZaDh8+XG/PPvvsXL7C0qCSeLyvC5aD89PqeFQR6vqduuty16aK2ITQOU1RGApXDhUtOy/3WNOJDQiN9iR1HY9HWV/b0ytbenfH5WCIbi2jy67LcZBoTxPCo3VrCYIgCMIcmdPw0NNPPx1jzLRszgsvvDAti1Nx5plnzni8tZaVK1f27P/MZz7DjTfeyL/+67/y2te+9ojX0mg0aDQac7n8pYkPsQi46scDnUxPqBoNVlkfSvmKWt6qXhO6ZmsFpslaPmi06m10qFQo+weGmHRREHQoZauOUyu6tip3Fh3Za6bZWkqXFxQD4oDU6wiCIAgvnzn9L3OapmzYsIEdO3b07N+xYwdveMMbZvzMpk2bph3/wAMPsHHjRpKkU1z86U9/mr/5m7/hO9/5Dhs3bpzLZQndVHO1yiBHed+xoju6holW+1QMfjx1HU/hDC5UwY+uXVo9s7VKt5bRPnZcVl1urXJwKFWXZVN1XO44tmp7et2Hp1PDM03WApG1BEEQhJfFnDWCrVu38o//+I/cfffdPP7441x99dXs2bOHLVu2AFFq6nZWbdmyhWeeeYatW7fy+OOPc/fdd3PXXXfxkY98pD7mU5/6FB//+Me5++67eeUrX8m+ffvYt28fY2Nj8/AVlxizDRLtrt0p52tVtnRdeHQ5VFQXZR1PKW1Fi3pZx1NEWcuVnZczbyiCni5rqThItA54qjqesoYn2E4dT1XL48tgp+66bHXMTJWSFkpFx1lZyxP/UNVwsOs2FVlLEARBOAbmJGkBvOtd7+LAgQPccMMN7N27lwsuuIDt27ezbt06APbu3dvTk2f9+vVs376dq6++mttuu42zzz6bW265hXe84x31Mdu2bSPLMv7wD/+w52994hOf4LrrrjvGr3bqEkJsDEjp0Apdc7WqZoQqqJ4mhPiOpSrGTFVWp3xEoafIS5r4dyrHVjx56cbqbjxYS1qha9/U5oNTnFtax87RgiAIgjAPzLkPz2LllOjDUzG1H48xKGPio7WQJrEfT5oQmim+L8H1Jbh+S9Gvyfs1+UDZj6cfioFA0R/w/Q7dX5A2i7ofT3+SxX48xtUBj0fRKhJazjKRp4xnCROtBlnL4ictesJgJhV2XGEnwY4H7CQkEx47GbATDtMq0JMFOpvSj6dwsR9PUdSjMmIhtvTjEQRBWIosyj48wiKl6rhcFS5PkbZUZU/vsaN31fSUz/GK4DTOdep46q7Ls/XkUQGjQ2eQaGVPr2t2ojurrt+pXmsFVk+v4xF7uiAIgrAASMCzFOkuXIauR6Y0Hqy2qDspH+Usgormramzteit46mlLDrztWoXVunS6pmvNbUBIXTcWlMlLkEQBEGYRyTgWWrM1HG56stTdlwmhF5behn8VE6tEGLHZXeErssV3YFP5dQK5YYqMz2VPV3TcWiZsvlg1XunDHakH48gCIKwEMgvx8nIbCMmZuu4XBUsFzHQ0dUg0W6Jqwx24iDRuDmvuoIePWuWx2jf03GZbjmrHDPRsaRX+8uOy6WcFYeK6in9eMqnkvkRBEEQXiYS8Jzk9IyY6N7fFfyoacEQ0+StOstTPg91Hx5VO7ZmonJqVdJWfNHl1poiaYXaoVWfYLqcJf14BEEQhHlGAp6Tmdqt1OmAXI2YKF/UElYtbQXKgaIdGauq5aGq3/Eda3psxKzKmVpTsjtTbOpKETM8U6zpHSmr67Gq46kaEFYnmEnCqv+A3K6CIAjCsSG/IEsJ3zX6YUqgUw0PrZ1aM3Vcdgo8BK/wLg4MnVa4PCXbo1TAao/R5eR0XXVcDj0dl9HU9TxTOy7H93Un06O7NpA6HkEQBOFlI78aS4WeYmU/7T01RcrqlrR6tyr90itrzVa0XFHP04Jep9aMLq1OoNPj2KrP1fX6SBkfQRAEQThK5NdkqeG7ZK4p1nRVTamvpS3KwaHTNz9DDU/3LC3o9OKZzZoeuup5evd11fLUWZ4p1nSpzREEQRDmEQl4ljo9slZV01PKW777sSsAKvvyVE6t3qnp07M9Pdb0UtYKlaxV9eKpJS1V29SrbE+V8Qm6K/CBWZ1aUrgsCIIgzBUJeJYiXUXLPY8wXdqCKdkdVbq0ek8ZpljSK/TU7E79hzrZnJ6t5/0ut9bUoOUITi1BEARBmCsS8CxVpkYsVQaH7nqdLmlrBpt6dw3P0aIqGQvqYGamTst1sFMVLdcffqk/ILesIAiCMHfk12MJ0dN8sKLqrNzl1uruxVPZ0btt6pU13Zddl2tr+izBTyVpdWSt6oRdjQe75Kxua3pVtFwXMk+ZqSUIgiAI84EEPEuYMDX4iTt7Mj1T5S01LV46SpcW051ataTV/eerrI/m5Qc0ku0RBEEQjhL5xThZ6QpmpnZbntp8sH6cWstTP68+2PUY1AzJItXzCJ0anh7qoCe+11PL0/0+dGV1up732NKP0ItHEARBEI4S+RVZKkztvdOFqoKf2ok1vSfP1Fqe7qnpocel1bllqmCn6rhcDRBVM9jQZ7OmSy8eQRAE4XggvyQnO0cIdHqP6834VDZ0NSWrQ1DTmg8eDd1urfikK/NzpOaDU+mZqSU1PIIgCML8IAHPUmaWwaJTVaiZVCmgq+Hg3N1a8cRzO7yWryTQEQRBEOYZCXiWOl21O911O0caMdFdxDxbhmdmt1b5pOquDDPKWrN2W34J1JGOkSBJEARBOAIS8Cw1/EtIXN0NB2d7f8rT7qBnpuaDqitF1BN3dJQxQRAEQTihSMBzCtNTwwNHDoRegm631rQYZ0o9z5EI0ntHEARBWAAk4FmqzDAxffZj6Q12eiSt6uPqqAuYgRn774SpMlf9fIbanRnmafWcXgIjQRAEYQ5IwHOKomZqSghTAp+5BRVqtupnQRAEQTjB2BN9AfNF1VW4IH9Z0szJRRWQ6NJKXlbYhIAKOvbd8R68JzhPcAFXOFxhcbmmyA0uUzircBacBq8CXnsCDhUczuc4l1PYuGnrULqoryD3nsJpisLjCo9rB3zL4dsWNWmgpXEtoK1QbdBZQGUB8gC5h8KhC0dwBcoV4HKUz8E7gs8heEIoILjyeaizV93Pyx3Ha+EFQRCEeaIgB2aZDjCPLJmA58CBAwD8B9tP8JUcR7p76HigikNaJ+ZyBEEQBOFYOXDgAMPDwwt2/iUT8KxYsQKAPXv2LOiCLXVGRkZYs2YNzz77LENDQyf6ck5qZC3nD1nL+UHWcf6QtZw/Dh8+zNq1a+vf8YViyQQ8uhw/MDw8LDffPDA0NCTrOE/IWs4fspbzg6zj/CFrOX/oBR4jJEXLgiAIgiAseSTgEQRBEARhybNkAp5Go8EnPvEJGo3Gib6UkxpZx/lD1nL+kLWcH2Qd5w9Zy/njeK2lCgvtAxMEQRAEQTjBLJkMjyAIgiAIwmxIwCMIgiAIwpJHAh5BEARBEJY8EvAIgiAIgrDkkYBHEARBEIQlz6INeLZt28b69etpNpts2LCBhx566IjHP/jgg2zYsIFms8m5557LHXfcMe2Yb3zjG7zmNa+h0Wjwmte8hn/5l39ZqMtfVMz3Wn7pS19CKTVta7WW9hCvuazj3r17ec973sOrX/1qtNZcddVVMx4n9+T8rOWpek/C3Nbym9/8Jm9605tYtWoVQ0NDbNq0ie9+97vTjjsV78v5Xke5J49uLf/jP/6Diy66iJUrV9LX18d5553H3//93087bl7uybAI+epXvxqSJAl33nlneOyxx8KVV14ZBgYGwjPPPDPj8U899VTo7+8PV155ZXjsscfCnXfeGZIkCV//+tfrY3bu3BmMMeHGG28Mjz/+eLjxxhuDtTZ8//vfP15f64SwEGv5xS9+MQwNDYW9e/f2bEuZua7j008/HT784Q+Hf/qnfwq/8Ru/Ea688sppx8g9OX9reSrekyHMfS2vvPLK8Ld/+7fhBz/4QXjiiSfCtddeG5IkCT/84Q/rY07F+3Ih1lHuyaNbyx/+8IfhnnvuCf/93/8dnn766fDP//zPob+/P3z+85+vj5mve3JRBjy/9Vu/FbZs2dKz77zzzgvXXHPNjMf/5V/+ZTjvvPN69v3Zn/1ZeP3rX1+/fuc73xl+93d/t+eYSy+9NPzRH/3RPF314mQh1vKLX/xiGB4envdrXczMdR27ueSSS2b8kZZ7ssPLXctT8Z4M4eWtZcVrXvOacP3119evT8X7ciHWUe7JDnNdy7e//e3hfe97X/16vu7JRSdpZVnGo48+yubNm3v2b968mZ07d874mV27dk07/tJLL+WRRx4hz/MjHjPbOZcCC7WWAGNjY6xbt45zzjmH3//932f37t3z/wUWCceyjkeD3JMdXu5awql1T8L8rKX3ntHR0Z4p1afafblQ6whyT1bMZS13797Nzp07ueSSS+p983VPLrqAZ//+/TjnWL16dc/+1atXs2/fvhk/s2/fvhmPL4qC/fv3H/GY2c65FFiotTzvvPP40pe+xP3338+9995Ls9nkoosu4sknn1yYL3KCOZZ1PBrknuzwcr/3qXZPwvys5d/93d8xPj7OO9/5znrfqXZfLtQ6yj3Z4WjW8pxzzqHRaLBx40auuOIKPvjBD9bvzdc9aed09HFEKdXzOoQwbd9LHT91/1zPuVSY77V8/etfz+tf//r6/YsuuogLL7yQz33uc9xyyy3zddmLjoW4f+SejLzc732q3pNw7Gt57733ct111/Gtb32LM844Y17OeTIz3+so92SHo1nLhx56iLGxMb7//e9zzTXX8Cu/8iu8+93vflnnnMqiC3hOP/10jDHTIrcXXnhhWoRXceaZZ854vLWWlStXHvGY2c65FFiotZyK1prf/M3fXLL/53Is63g0yD3ZYb6/91K/J+HlreV9993HZZddxte+9jXe+MY39rx3qt2XC7WOU5F78shruX79egB+/dd/neeff57rrruuDnjm655cdJJWmqZs2LCBHTt29OzfsWMHb3jDG2b8zKZNm6Yd/8ADD7Bx40aSJDniMbOdcymwUGs5lRACP/rRjzjrrLPm58IXGceyjkeD3JMdXu5aTmWp35Nw7Gt577338id/8ifcc889/N7v/d6090+1+3Kh1nEqck8e/f0TQqDdbtev5+2enFOJ83GisrXddddd4bHHHgtXXXVVGBgYCD/72c9CCCFcc8014f3vf399fGWlvvrqq8Njjz0W7rrrrmlW6ocffjgYY8InP/nJ8Pjjj4dPfvKTS95qGcLCrOV1110XvvOd74Sf/vSnYffu3eFP//RPg7U2/Od//udx/37Hi7muYwgh7N69O+zevTts2LAhvOc97wm7d+8OP/nJT+r35Z6cv7U8Fe/JEOa+lvfcc0+w1obbbrutxyp96NCh+phT8b5ciHWUe/Lo1vLWW28N999/f3jiiSfCE088Ee6+++4wNDQUPvaxj9XHzNc9uSgDnhBCuO2228K6detCmqbhwgsvDA8++GD93gc+8IFwySWX9Bz/ve99L7zuda8LaZqGV77yleH222+fds6vfe1r4dWvfnVIkiScd9554Rvf+MZCf41FwXyv5VVXXRXWrl0b0jQNq1atCps3bw47d+48Hl/lhDLXdQSmbevWres5Ru7J+VnLU/WeDGFua3nJJZfMuJYf+MAHes55Kt6X872Ock8e3Vrecsst4dd+7ddCf39/GBoaCq973evCtm3bgnOu55zzcU+qEMqKVEEQBEEQhCXKoqvhEQRBEARBmG8k4BEEQRAEYckjAY8gCIIgCEseCXgEQRAEQVjySMAjCIIgCMKSRwIeQRAEQRCWPBLwCIIgCIKw5JGARxAEQRCEJY8EPIIgCIIgLHkk4BEEQRAEYckjAY8gCIIgCEue/x9XPxlssE0uNQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phi_visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2048546e-02c5-440e-84b7-91abcf43732e",
   "metadata": {},
   "outputs": [],
   "source": [
    "step_solver(300)\n",
    "phiprofs.append((t, integral_dr(phi)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "79e2c6d4-bc8c-4e8e-9faf-60719118371f",
   "metadata": {},
   "outputs": [],
   "source": [
    "step_solver(500)\n",
    "phiprofs.append((t, integral_dr(phi)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "f69518d9-b522-4a0c-ae99-15dc36398106",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.999999999999831 2.1205750411731096e-10 2.1205750411729896e-10\n"
     ]
    }
   ],
   "source": [
    "print(t, initInt, phi.domainIntegral())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "14ea960f-2d57-4c5f-912d-758cd3b03617",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t = 10.0 s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAACkCAYAAACTpfSdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABM8UlEQVR4nO29f7BcdX3//3y93+fs7r35BSQhBAkY2s8IljqVpNXIl6F/aKi2HXHsaLVS61imVKYCmU6B2lagM2C1nU4RAsKkWuuAjNpWnMlY0mnhgxDriNFphRY/FcVCIk0quUlu7u457/fr+8f7x3mfs2dv7k3uTW7ufT1m1rN79uzZs++7sq+8nq/n60XMzBAEQRAEQVjEqFN9AYIgCIIgCPONBDyCIAiCICx6JOARBEEQBGHRIwGPIAiCIAiLHgl4BEEQBEFY9EjAIwiCIAjCokcCHkEQBEEQFj3Zqb6AucJai5deegkrVqwAEZ3qyxEEQRAEYQYwMw4dOoRzzz0XSs1fHmbRBDwvvfQSNmzYcKovQxAEQRCE4+BHP/oRzjvvvHk7/6IJeFasWAEA+P/wNmTI5+7ERKAsh1q5HFh7FiY3rMTExhyHNlqMXXAIr1/339i88gf4me6LOEdP4kyt0aMMGTQ0DUeqhi0Ocx8vlcC3+6/C/33lNfjGf58P+4PlWPFDYPl/lxjbNwn1vxPgw4fBU33YQQmwBaQptiAIgrDIKFHga9gZf8fni0UT8AQZK0OOjOYw4AFA0FAlgfqM8SMKfLAD7Fc4mo1hd3EG/v2MC7Fm+RGs6hzF8ryPrjLoqhJdVSBTFuNqgHHdx7gaYJnqQ5HFVNbBPnsGaHwMK1cT/tdoTHZzmNWEqVcvRza5FlmfoQoGGYAMgyxDlQxVWlBhoQoLGpTu1i+AQQEu/LYsgaIAGwuwBVt2QRMggZMgCIKwcPA/SfNdjrJoAp75hI2BnZwEGYPsyFGcsW8cK/9jDGZZF2Y8g+kuR9FZgZczwj5NsBlgc4LpAKZHKMeBYhmjXMGwK0p0V/RxxvKjWNWdwrK8j/9z1n7grP2wryGUVsGyggXBsntsWKFfZjha5Dgy1UH/aA57pItsQiOfIHQOAt2DjO6EQedgiWyiD3W4D5qcAk/1gX4fXJbgwmWK2Bj/wSTwEQRBEJYGEvDMBGYXMJQlcOQI8D9uN8EtYAYARAApUJ6Bsgw01gMtG4ddsQzlmWPor84xuVpjak0HU2ty/HhNB/2zMmw4o8T5vZ/gp3r/gw2dA1irJ3CWmsIyZaMwN8XAIc7wP2YZflSsxv+bWofnDp+N519ZjQP/uxyD/+mgOKBQHtAwOaGrCR0iKADEDPZyGFkLNgAoyfYIgiAIwhJAAp65xhgX/JQl0B+A8gzZIQXWBCYCWQU9IPQnu5g4mOPfVizHfy4/GyuXTeGMsaNY1TmKVfkUlmV9jOkCPVUgJ+OyPDbDYdPFK8UYDhddaGXR6RXon6EwpTOUY4T+GRr5EY1ssoNsajmyKYbuW6i+hS68FDYwoMLdUJSgogRKA4SgriydFGaMywZ5OUxkMUEQBOF0RQKeuYIZYANmAlsGGQua6gNHJkE/ydH9cQfdbgcruh1wL4cdy2G7GmVPw3Z6KHtjONQ9Cz/pEkyPYHpA2QPMGLvbuAV6BnmvRKdbYqxToJuVeNVZB5GttiBifxmEkhUKozEwGoeLDP0iw2CQwRzNQEc70JMK2RFCNgnkhxn5YaBzxCI7YpAfLqEmC6jJPqg/APoD8GAAFC4YYmNcICTBjyAIgnAaIQHPXBMCH2vABYCpqeFjvPylFUFr7SSwTgfU7QC9LnjZGMyKLorlOQarMvRXKQxWZe52Ro4jqwrwSmBZZ4C1Y4dxTm8C6zsHsSY7hBX6KJapPjQYA9aYtF28Ysaxv1yBl/pn4KXJVXjp8Eq8MjGOqYNdZAc1OgcJnQmNzkGF7phCPqGR5Qr6iAYpBVIKTH2EcjIGQDBgq1zQQyRBjyAIgrCgkYDnVMAMwIKtcpkZY5yUpAhEBFIK2h9KDJDNoAoF1SfoowrFoQ6mluf44fJx/PeyMzA2NsDyXh8rO32s6ExhPBugqwwUuQxMaTWOmhyTZQeHii6YCVluUI6XKAFwplCOEQYrCFOTGbKjGtnRDrIpGyUx3TeggYUalK1yGJelk/OCHGa9JMbcLonFdRAEQRCE+UcCnlNFlMAs2ChQWYIHAxBNAodzUJZBZxmybge9PAN3cieFdTPYbgYzpmG6CmWvC9PtoegRXu4RXuoBpgeYHsN0ATtmwR0L6hpkHYO8UyLTFst6A6xadhSKqqDDWIXCODlssshQFBqmr4G+hprKoKcI+ighOwpkk0A2yciOMvKjFtmkhZ4soadKqKOFq1/qF84aXxROEgtyWJtVPqyJIAiCIMwDEvCcaqIEBpcdAeoyWOhLQAqkCKQ1Mq2RZxnQyUHdLtDJwWNd2PEOyuUdlMs0Bis0BssJxQqNYrlGsSJDsdLALFcYX9bHWF7grN4kzuxO4qz8CJZnfXSphPZZoSmb+wLpcfxkMIYDU8vwyuQYDh/pYepwDnVYIz9MyA8ptz2s0TmikR/OkB3JoCYzqHwAmtIgrcGqABUDsFIgFGCmuiwWPqsEPYIgCMI8IAHPQofZBwKVBMZwTjAqCUwKBNewScH9QcnkgAVUqaAHCnoK0FOE8kiGcpnGkfEODo+PY29vFXpjA4x1CoznBXpZgY42yMj16SnZFT73jesBZJmgtMsWWUsolILNgXKMUCwHBkc19FGFbCpD1u9CTzH0lHEuscKA+gZqUAJFCVU6SQxlCbbWudpC9meUHAZIRkgQBEE4LiTgOR2IP+wWXFqAlAsIBgWgpkBaA774WWUZdCdHJ8/AucsC2V7mpTAN09WuGWIvQ9nNYXpjGHQJR3uA6XoprAPYngXnDHQNVG6hM4Mss+h0SvS6BdQZ7posuwaJxij0SwVTathCAX0NGhBUvwPddwGXngKyo4CeYuRH2dUHHbXIpgzUlIGaKqF8x2gqSt8xOrHLGwNYW5fEAHGLCYIgCMdEAp7TifCDHiSwsDvcSeQvAC4QCjJYlgF5BnhHGHdzcK8LHsthxjMU4xnKZQrFuEKxDCiWa5TjQLlMwyyzKJeVoLEC3bzEil4fyzt9LM/76OkiFkhbVuhbjSlfIH140MWhfgdHjnYxdTQHT2bQRxSySUI2qZAdAfIjCvmkRnbUIjuSITuaQ00Vrg4oz3wdkAK0coFPUYLg7PdEHDNf4hYTBEEQpkMCnsVEIn+BfBBglZPAgIat3EHM7nWGoUzmZnUVCmoQsjIEM6lRHlYoxjK80uviYG8cebdEp1Oik5XoZAa5srEA2jKhsAql0SiM85uRZnDHwjDASsFmgOkQzBihmCJkRwm6r70rrAPdZ+i+geob1ygxusNKUGnqjRJNIoGFLJDIYYIgCEKCBDyLjTQLxATAVBJYUQJ9Z31n7YqJkXsZzGeAuJMDeQbbybwU5txgpqdgOk4KM90MptuD6QFTHWCyC5guw3bcjXMGcgvKLVRmoZRF3ilB3dLVIDG5GMsoFFbBlgQuFFAqL4Np6AGg+rmTw/pOBsumAN1nZH2GnrLu1jdQ/TBA1dvkCz84NQ2GgkNMbPKCIAhLEgl4FjOzkMBIuWaIQQJTWQaltXOCeUdYtMaP5TC9DOWYRjnuevgU4+TvA+UyhhlTMGMWZsyAukDWKWJ36LG8QK5MLSNkrELfZJgsckwVGfr9HP1+Bj6qQVPaZYAmXSbI2eIVsqPsCqSP5siOli4bdLTwwU/mgx8X+FBJziFmDEBuphgpL4cBIokJgiAsciTgWco0HWCqLoGB2dXLWLcN8pdiBhl3U0ZDlco7wnw2pu/GY5RjBNPTMD2LI90ck10DnVtkuUGWGWTKZX8C1ioYJpSl9pfG4CxIYQpMAIdJ9F1COU7OCdZXyPraN0js1qWwpEmiSmeGGQtYM60zDP6z+zv1dRMEQRBOKyTgWeqkDjADpBIYSIHVoC6BRTeYrmSwPAN3M3Ang+3oyg3WJSeHdQhlN4Pp5rBd5wYrO8Cgw7AdwOYM7jA4s0DGIM0gbUEK0GMlaMxlpZgJhQUGRoENASWBggxWENRAQw0IeuCkMNVPJLC+c4WpgXU2+YGXwvq+LihIYWXpnGFpMJS6wwBxiAmCIJyGSMAjONIf7DYJrOkAU+SCH/JbrUG+HoiSeiAng/nu0D3tZLAxQtnz8teYK1x2w1IVzBiDuxbcBXTXIM8N8rxErg1ybaF9RsgyuQDIKAxK1xW6LDIMBsp3hlZQU4TMF147SzxDT7nC6GzKQk/lzhLfd5Z4GpSVJb50tn/SBigIbIybKWYM4LNhbLkqDg/rI4GPIAjCgkQdz4u2b9+OjRs3otfrYdOmTXjiiSemPf7xxx/Hpk2b0Ov1cOGFF+K+++4bOuaVV17Bddddh/Xr16PX6+Hiiy/Gzp07j+fyhPmA6xkNtlxzRrHxjQQHhesPVJSgQQHqD6COFlCTBfRkifxIifywReewRecQozPB6BxkdA4CnYOE7k8U8lc01Cs57E+6mHqlh8OvjOPgxDL878Q4Dkwsw08Oj2NisodDR7s42u+gKDSsUVGhg2LYjoXtMsoeoxxnFMuAwUrCYBWhfwZh6kyNqbMyTK3uoL+6i/6aMRRrl8GsWQGzZiXsmSuBs1YBZ64CnbkK6oxVUCtXQK1YDrV8mbstG4fqdaG6XTf8NXP1TpRlgNLuRlTdBEEQhFPGrDM8Dz/8MG644QZs374dl112GT71qU/hrW99K5555hmcf/75Q8c///zzeNvb3oZrrrkGn/vc5/Dkk0/iQx/6ENauXYt3vvOdAIDBYIC3vOUtOPvss/HFL34R5513Hn70ox9hxYoVJ/4JhbkjKYIGkZPALDubuyrBSdEzT1KUwJBl0FpB5zmQaXTyLDrBuKthc+cEs10VJTDTcXU6pqPdTLAOXEPEDqPIAM7ZSWEZA7q6EQHUNUAX4GVwNUGWAONuZAhUuJsqCWqgoAbeFTZwEpiOW38LFvmBHbbIN+uCQsdoscgLgiAsKIh5dv+1fcMb3oBLL70U9957b9x38cUX46qrrsKdd945dPxNN92ERx55BM8++2zcd+211+I73/kOdu/eDQC477778IlPfAL/8R//gTzPj+uDTExMYNWqVfhFvB0ZHd85hDmiTf5KgiEngakYDFGeV7VAndwFQb4zdNnTMD2Fskcoe+QHoxLKHlw9UI9hOgzuuuJm6vjO0Lq6UTIg1VoFYxRM6baxK3RBbhp9P2zhZLApdltfB+Ss8CPs8MZU3aHbaoAk8BEEQRii5AKP4cs4ePAgVq5cOW/vM6sMz2AwwNNPP42bb765tn/r1q146qmnWl+ze/dubN26tbbvyiuvxI4dO1AUBfI8xyOPPIItW7bguuuuw5e//GWsXbsW733ve3HTTTdBa9163n6/j36/Hx9PTEzM5qMI80mjASJbrjnACPXmhwznAAsuMGI3C4wsg2zqBCOYglAWgBq44Ef1CbpLvg8QwXYUTM4wmQVlvvhZe5eZD3zY1/+wJYAJUAzWzv0FBpgYUARWAGcEm7tCa90n6B4h6yvovoIqMui+AQ38rDCf+UFpYuaHjIkBEBnjMmLGuM/aCIIAcYUJgiDMF7MKePbv3w9jDNatW1fbv27dOuzbt6/1Nfv27Ws9vixL7N+/H+vXr8f3v/99/PM//zN+4zd+Azt37sT3vvc9XHfddSjLEn/yJ3/Set4777wTt91222wuXziZtMlfxIBB4gBzTRChlHeBVVkfrRSQusDyDJzrKIHZjpO/bIdgOgqmA5Rd8o/Ju7/8rcOwGWAzJ4NBA6zZXY8CoBjwGSIzTjAMFJZc0FW6Ia2qJFAJ6IECFYAuNFSQwQZeBosSmIUqOEpgQ/b4WqfoejaIxB4vCIIwLxyXS4saBZjMPLTvWMen+621OPvss3H//fdDa41NmzbhpZdewic+8YmRAc8tt9yCbdu2xccTExPYsGHD8XwcYb5Jgx8kv9lESRPEEfIXuS15Ozxluh4INezwZU/DdgllV7lePT14e7zb2i7HWiDuMDi3QA5QZqE0u75AxJXpyrrhqNYQ2CgUhXJ2+IKgBu6m++Trf7wU1lfQ/SCD2XoNUN/X//ip8RSmxhvj63+8BGYMwOxdYQzSqOaGhUUUZ5ggCMKMmVXAs2bNGmith7I5L7/88lAWJ3DOOee0Hp9lGVavXg0AWL9+PfI8r8lXF198Mfbt24fBYIBOpzN03m63i263O5vLFxYa4Yd6OvnLB0VN+QthywxYQBkGGMgtwxgNKgHli5RVAX8j2IKgOq4I2nYINidX/JwrGM0w2vUBAnH1pgBgvQQGfzHKN0Fkf1VEAAGsAdbks0kE3VHQXUAPCKrP0APtAh+f8aG+l7rSGWHGbxsyGJhBM5XB0vUVBEEQZhfwdDodbNq0Cbt27cI73vGOuH/Xrl14+9vf3vqaLVu24Ctf+Upt36OPPorNmzfHAuXLLrsMDz74IKy1UMr9C/a5557D+vXrW4MdYZGRND8EUJe/FAHGglVZyV+UOMC0gsrcFHidaXCmkQf5K9fOCdZRMLmKMpjp+LqcKH85SczmoZMzu/vayWCs4VxgimMjB868Q6wHGIZ3gsG5wEqCMgAVcINYC7gapIELvIIEporKCaYGFnpgQYV3ghU2OsFgTF0Ga3SJHpLBgNGF0bX1FgRBWDrMWtLatm0brr76amzevBlbtmzB/fffjxdeeAHXXnstACc1vfjii/jsZz8LwDmy7r77bmzbtg3XXHMNdu/ejR07duChhx6K5/zd3/1dfPKTn8T111+P3/u938P3vvc93HHHHfjwhz88Rx9TOC1oND8ERstfpBInWEP+it2gtap1g4YPhJz8lVUBUBiMmkpfwQbfJZiOy9Zw7gajsmYgs84Gr9yt+gguE2Ss6wQNHwBFCax00pcqkEhhLvMTrfB9W3WEDlmgxAZPPgPEpb/fkMFcMTS5Iu2mDObXcmi9BUEQFjmzDnje/e5348CBA7j99tuxd+9eXHLJJdi5cycuuOACAMDevXvxwgsvxOM3btyInTt34sYbb8Q999yDc889F3fddVfswQMAGzZswKOPPoobb7wRr3vd6/CqV70K119/PW666aY5+IjCoiB2FUQy+8tLYUH+0to5oOJLVLxPzK7WLMwCY4Cs9k4wn5kxBGVcUGJKV7BsS4AKAucuAHKuLQZn5LI8isEKrRJY1OKUl7oyhvUeNSZyBrFUAvOOMJtrqELBdizUQEF3tXeCZVXxc3CCWTsc/LB1W2unl8HcoLThdRYEQViEzLoPz0JF+vAsMdJC+LaCZy9/peMvoJXrgqw1kLlMUM0BliuX/ek4F5jNqZK9coLpIMph4RYlsCzU73DchmAHcFZ3Yh8EWYAsgYIEZuFqjgrnBFNlkL6c7BVkMFUAumCoAUMVTgJTgyB9+W1LPRDbkAFKpK+W3kCA1AMJgnDyWZB9eARhwTBK/vKBkJv/dYxAyNvgSWmo4P7K9LD85W+VDT50gfZBUDexwXfI1wC5rE2s9VEAaxszQfHqGS4bZF2GCSYEPhQLrV3tD/lu0MEC76fTD6wLgLwERqWFip2gW4Kf0oCsAYyq6n+SnkAigwmCsFiRgEdYXDSkryHnl4L7kYdXnKwbSRHdX172Gpa/GGQ1jFUgo7wMBi+DMUzp635M6N0DUO4lsKTwOWZ+KA3YkgiI4CQyzWAm2KCLUXjON0TUgNWAzRR0xlA5weYKemDBmaoCnzIDEgmMvAQWMz7N4mcf/CBmfnwAGTI/kMyPIAinJxLwCIuP2o/wCOcXUG9+qDVYqar5oXIZIPJZH5Vpl/XRelj+8sGG6XrpKw8uMPjnvAyWeQksY3AGH7j4GiDlYwviKIvZHC6wsuwzQO6x8vVFqkRsiOiyQUH28g6wgn0DRHaZoNK7wKIDrMr8qIYERg0JjJoSGFB3gvnH7X8DQRCEU48EPMLiZlrpyz8mBaAAKaoGoBKBlao5v6IF3tcAcZ4ldUCNLtCdRP5Kuj+7bb0OqBYIaXLZHeVs8ENhA7v6n3YLfL3nUKz7GYQaIDcPLNT/jLLAj5LAYHmkE8yta0tTxObfQBAE4RQhAY+wNJmp9KVU3fmF0sle8KYsW5fB4KUvsgzD2mdo3M1YAhmGMaFg2ctfuSt8JkO+4Nk5tyr5yxU9A6gKn+EvQLFriJi5DuY2XKkbHhYdZKwoOsL0wEltqlDgTPmMjwIV2mV9QuGzsT74KWPhMxntMjvT1P8Aup75gWR+BEE49UjAIyxdZiV9lZXbSylAK1cInTi/KNNQ2jU/dBkg3/ww1+COgs2qzE8lfRFMzPaQd36Rn/1V2dlZ+1IfP/uLY00PwLn7HEECQwiyjO82XboAS8VCaAVVcswGRemrTBxgRSV/oUwcYKXLAFHpJK6q/scHSqOGowLTS2BDfw9BEIS5RQIeQQDapS8OeR0/7BQAqdI915S+GvZ3aO07QGvoEASlHaBjDZC3v3dVUvtDSf2Pd34F2St3spdNLPBhCKq73BBcINrfYX2tT7C/h0aIZSV/Vfb3KgAKMpgbhZGMwygMqLSthdBp8DOcBZpGAgPECSYIwrwiAY8gjKI56wvVj3Sr64uo1nvQucAskHnJC5UTTLGbG+cyISrW5hjfsNDdZ9evxzrHlzVwM70ydhJWIoGlRc/VBVRXw8o7rjRgmaGYXE4raHNElWVeEbRy2SOtAVbKNVrUBKVdYTcVFpT5TFdpQEoBVrvMj1FgpUE6yfqEwmeyMZBslcCk/kcQhHlCAh5BOBZtgU9T+vIZHy7LOOYiSl9h7EUifbkMkJO+slyDM+UyP5mqOb9sx0tfWZXxMTn5bE+QvWi48WHo8JwEQc4RRiDtpsUTo6ol8p2mnfSVur+4ygSV7BxgQfoqnfylBhZkKvcXSluTv1Caat5XyAD5Amj22R8yIasmjRAFQZgfJOARhJnS/JFl4wqDDRAcXyAFLsrY7JCBOOuLfRDUan33+zjPGvU/PgDK3JDTEACZDtUCIJvUAbm6n8T+TpTY312mKSR/iOEzSt76bqsgyAU9zgUWpC8qg+QVrO/ay18WVHLl/irtsPxVVnVAo4IfyQIJgjBfSMAjCCdCmv0Bhh1f5OzcDONqfIA4YstJPX7eV5C4gNj0MEpfDBArH5x4+YtdmxxXo+PngfmAxeY+RrAEtj7jw75hoaqcX4n1zH8G71xzb+GyLaDosWIKddwEJjcPDH7LSkEP2MthBGgCKYJSCqytk7y0Ako1LH8Z64JC696f/fwv56TzwU9YTyS1P3G9JfgRBOHYSMAjCHPBdLIXnJ2cARelKEokLw0MCp/xcdNEg+RFmXOEqZD1CbKXl8Bc9kf57A/BZsHxFSQwwGZe/soq+StKXwpV00MgBkBMALRrveMCL5/1sZWdXhlKskCJ9FW77wqfKUhffuscX24MhrPB26r/jy905rJ0mR/rs0KWZyZ/xSBOJDBBEOpIwCMIc8m0shdQd3w1ZC8fBHGL44tUo+tzkL206/rMeXB8haGnfhvnetWt7zZzgU/N+u4DoHT6Ozc+W3B9VfU/zu0V63+KYH/nxAbPUKWuBp/GAKje+Zlban9are9aH0P+EgeYIAjDSMAjCPPNsWQvYNjxpXXMZABeBiudJT7O+gJqPf1UODA4vZgApih3eZ0q7rOZO56tk7ysZkC7YMbJVQ3nVyz88Q+VOx/rsMMFVYiyl9vJ0QHm9nFwfhGgiHw3a/IFz97qT+R6HJXG1T1R6RosErnAy5hK/gpjLxTgbWgifwmCMIQEPIJwshgpeykAppK9SPkppF728uMtWFGt2DmVvXhU1if3kpd3fLn7obFhNePLSV+V28umdvda5geVbORdYEY7x1fo+xOdX7aSvqrZX+x7AfniZ79PFS3Sl7Gx6DkUPrtuzzZmf8hnf8Du+NAFOgxApVrzQz0sfzWlr/TvJAjCokICHkE42RxL9vJ295rs5Rse1gadjpK9tA+AouurCn44rzo+c5a4vTKu3F5Z6PjcrPuBywKFAMhb32vSV6j5AVwAlNT+qNLfD7U+XvoKwZD20pdKOz83pK+hrs+mMeg0HXzaOgLDF0BrtLu/3CK3/50EQTitkYBHEBYCzdleAJrzvQAXaJAf4DlS9mLtkhdAzQGmmGEZiCoa+yaKwafODGICscvYWC+DxQQJA4r9jC+GC3aUjwvSjouB0AsI7rz+MoAsVeJ89+Uh9xe8pDYsfUHBSV9KgUvj659c0EjWeunMuoJnny1jv74UnV8c13RI/pLARxAWJRLwCMJCIf1hHerubFy2x7ggAMx12csXOnPq9tIK8Bkf8uMtVHB7Bdkrq7I+oeDZZXx8wXOWjLWI96mSvVQqfVUzvkL9D/uejJaq+qCQ8TG2KXuF2V++4DnIXaX221T6YijjMj/kM0Bkqp4/sDa6vmrSl+/7Qw3pK24BcX4JwiJFAh5BWIg0Z3u1OL0ANAacVk4vxEGnLXU/Snmbu6/70bpV9hoZAIUOz+mQ0wxglUx7V1yr/4mFzuF/Y8NDfzNp08Nmx+dmAIRE8qqkLxf0sJv5Vdoh15fr9mzd9PdQ7NwmfR3L+eUeSAG0IJxmSMAjCKcDqeQFNJxePhuUOr2UAhsDhMaGAJh9IW+WAWXS5NAyGFlQoGBRyV5VF0IA8JIQgtRVaVnMlesLDEAnrrA44qLlcwUnWIglkrghXIICuQyRl72iWFfLJikoIqjCwiYSGBSi64tKE9eQjKkKxK3xp6a69BWcX0PSF7vXpRkfkb8EYcEjAY8gnC5MK3kB0elFKjY4JCKwUV4vcmMtuCxjBshlezK3P9MgnchePuvDWtXcXpxVWR+bVbJX5fZC0uMntaKjNt+L0/jNd4BmRYnjyylNMQOUjLxQJcceQDHr42Uvis4vX/BccrvjK8heIfNjqr4/QTKkpOFhTfpqc30B7c4vCYIEYUEgAY8gnI60SF7hbpS8wmBTAGG4aWp1j8NNSdWcXiPt7tMEQKbTCIBS2cvP9qrX/ZBLkiiOkpernQ6dqVHJXkndTwiEanU/bQFQAahS1WSvYHun0gVCYchpbdZXKntZ68d2NOZ9Tev6mqb4ufl3EwThpKKOfcgw27dvx8aNG9Hr9bBp0yY88cQT0x7/+OOPY9OmTej1erjwwgtx3333jTz285//PIgIV1111fFcmiAsTYas7i3TxluyFRwyG4mtG2HAZ8iAlKYaBOpvtX45foSELpJp6gWSLssuMHEjKYB0QKkLZFz35hjchI/SKH6OGSJV9QkKmaRmlqnqKu3HboT7uQL7jtSc++LtzNc05Zm38/usl898QSXWf61qAaOrmfI3CtJXEmRS4z+xafAjCMJJZdYZnocffhg33HADtm/fjssuuwyf+tSn8Na3vhXPPPMMzj///KHjn3/+ebztbW/DNddcg8997nN48skn8aEPfQhr167FO9/5ztqxP/zhD/H7v//7uPzyy4//EwnCUuV4JS/lGx0GyauR8SHf6yc4vaBdpkf5Ymebqyrrk1ND8vJZH113e1V9fYLcRbWAJtb2ADECisXP3qTGGqAskb5qslcqf8EVMyeyV5Xxabi9woyvuZC9AO/48tk3kb0E4ZRCzLP7f9ob3vAGXHrppbj33nvjvosvvhhXXXUV7rzzzqHjb7rpJjzyyCN49tln475rr70W3/nOd7B79+64zxiDK664Ah/4wAfwxBNP4JVXXsE//MM/zPi6JiYmsGrVKvwi3o6M8tl8JEFY/KSZhSTrQElmAk2Xlw+A2iSv8DhKXtNZ3ZvdnWtWdydvxZleiexVC35UveYnXr+Xu6LjywBgn0mqBT9VpkmVgEoDoLTLs+Eht1eQvmAawU9ZAsw1u3vs+VOTvUbY3QMS/AhLnJILPIYv4+DBg1i5cuW8vc+sJK3BYICnn34aW7dure3funUrnnrqqdbX7N69e+j4K6+8Et/85jdRFEXcd/vtt2Pt2rX44Ac/OKNr6ff7mJiYqN0EQRgBt//AcvIjPFLyYjsseaUTzktfDGxMMgzU1cmogXusB9Z3UuZYZByCDwqP4yT2RPYKUpcld2uJBTiVvnyA1JS8WNeDrEoCoyQrlUpd5KbSJ7do5Q+1T1mWyFxVC4BW2Qtol73cg/oHCg4zQRDmlFlJWvv374cxBuvWravtX7duHfbt29f6mn379rUeX5Yl9u/fj/Xr1+PJJ5/Ejh078O1vf3vG13LnnXfitttum83lC8LSpmWWFxAkr5bGhkDi8nLzvKBK17cncXmRdsEAJUXOpJzspTOX6XGuL2pkfCgGIrGfTxqgJEEL6xDQkM/4oLK5J9vgBiMLGAWQb3ZokxEXLuvDSb8fV/SsTJC8VHR6uQJn9gFdkvWxLbJXaGSYyF2hz8/0bi8My14AapNh07+fIAjHxXG5tKjxrw9mHtp3rOPD/kOHDuF973sfHnjgAaxZs2bG13DLLbdg27Zt8fHExAQ2bNgw49cLwpJl6IczGWQKIMzyIpUMM1UELstK8gqdnZPGhlDkevyEWiCfDZm2u3PWVvPjAyA9Xc1P0+aOWt1PGGMR3F5kER1exI16n1T6KikJfJzTywU6oe4nqwabhg7P1sle6XBTssNZMQpBkA9+Zuz2gt8nbi9BOCFmFfCsWbMGWuuhbM7LL788lMUJnHPOOa3HZ1mG1atX47vf/S5+8IMf4Fd/9Vfj89a6/8NnWYb//M//xE/91E8Nnbfb7aLb7c7m8gVBaGOGTQ2hqPqhBtwsLws4/zsAaBcUKQX4RoacITY5bKL8/1o0tfWquaAlQIFh/RlCY8R4qfA9ANNXNt8sGXHB8CfQ1dPu/QmWOTl/dQ3pe7r7qn69pnpvAL7J4fCngTEuOIQBWbhicd8kkkKTQ8CtNVsXcMZCZ2l0KAgnyqwCnk6ng02bNmHXrl14xzveEffv2rULb3/721tfs2XLFnzlK1+p7Xv00UexefNm5HmOiy66CP/2b/9We/6P/uiPcOjQIfzVX/2VZG0E4WSQyl3AUJaBlHU/7KTAfqKo+5H2tSrGgqj0k9zJ1fV4Zxf8KIvg9HLjLMg7vZKMj/YuL13P+FTW8yTL4/v7wFvU4xaJ3EX14CcGPVTFeKQrl5cy5KQv44ucE8krurxMyPTodsmr2d8n1DpFuavK+oQhsKHYGS3FzqQxWvZqZn3Sv6MgCEPMWtLatm0brr76amzevBlbtmzB/fffjxdeeAHXXnstACc1vfjii/jsZz8LwDmy7r77bmzbtg3XXHMNdu/ejR07duChhx4CAPR6PVxyySW19zjjjDMAYGi/IAjzzEzkLsDZ24Eod7mdlcOLQuAzwuFF2jUwrM3y8jcbtxS3qctrOABy3Z2VroIfDlJXQ/aqBUAhxvOF0TbIXQ2Le132olgD5Op+VOXyMlVfopGSl6/3GZK8bDa6wWGY7g7tZK8490PXg5/wNxj5txSEpc2sA553v/vdOHDgAG6//Xbs3bsXl1xyCXbu3IkLLrgAALB371688MIL8fiNGzdi586duPHGG3HPPffg3HPPxV133TXUg0cQhAXI8chdfo4XQbuhXD5YqiZvOZKpWLXn4jwvKFhYEAUJyR1pQVBB8oKXvKgueRFcMBRmf0Y5q+UCQpEzV288RHhP60aDuSRMPEEqe6maSEcMcAlAs1uP0sl8w5JX0NkMAO1mfTUlL2YXXCadnKPs1ZS8AF+cLkGPIARm3YdnoSJ9eARhnmnKJ22dhYH2fj7p9PZkcnu0eKdjLDL/ONcu26PJS16+U7KXvFymx2d+Gr18ai6vpKdPLdOTurzCfwVDp+c40gJVsXPM9HAt8+McX6icXgX7Y5LGhmnWp9Hbp83lxaHAuVHoPLK/D5DIkNLjRzi9OFl9eGSWliAIM2NUJ2df2xNSOU25i8g/ju6usnJ3Ne3tWkV7O6UOr8x1c3Yzvaju8Mqa9T6uu3Pl6qq7vIaCn4bFPcY+SeATnF5umCnFHkH15oZUa2pIhhObu6/9CbJXI/iBTSSvtMGhtcd2eY2UvBKnl0hegiABjyAIx8EIqcvdrctdrHxrijZ3l9Yuo4GG2uSlH9i648nVC1WyURSP/IstCETsXFf+OUvuWtn384vvFWSu1OKVEuq3/V1ukbr8xFOwZv/+4YXVNQBAqpZxLC6q3oZLgLyjjQEX3ITnmV2d07FcXqnk5QMckbwEoUICHkEQjo8Rzq5Q4+Pu2qqZoSI/e6rh7tJuthdrFd1d0NplOXw/HyoUKNPgwsQiZ8oUlKahAufQu6eSuUZkfFQledXmdzWKmwHv8EJ1LPn4Lri82JDP+LB3eoV+Pun9kOFRMfvTzPjAsJO6Gg6v4SLnZIq7Hwzb5vKq9/bRw3KXuLyEJYQEPIIgnBgjnF1AXe6iYKtuDi819VofbtT6EJEf41DJXZXFXdUkr1jjkzY01A25S3Oto3Nr8BMSMG1BUCh0bq31oVjbY6yv+SlD7Q+1yF0MKlUMfpy9PavZ22cS/MC6IaZgPw6EyGWJNGbe2BAQyUtY1EjAIwjC3DITZxcpV4uiEJ1d8JkfMs6pBHipS6sohwFe7sKwuwuoNwa0FNxc7qCa3EUuCLHpGZrJDuXeINY0N4qcY9bHv3G8Hk7lK/dRorwGJyMpkHd6ASo5HlAgtmBNLkBk5Vxd4Zoalznk8LIKpP01JVPbXVdsuOCIZtDY0L1Igh5hUSEBjyAIc09b0IMqqxALa23SyND/wLLfOqeXAVvtJpMnvXzImKqhYZC7fKaHCrdVg2R+l04LnDm6uFzmZ4TclY6vCEXOgRbJC14WsyEg0i6rE4qb2bu7bEZVXx9DoE41ud1lfnyRc5u7K/T0GZHxCdPbKbNxjtfIURYAKPj2ARxT8pLgRzjNkYBHEIT5oenqAlqdXbWhpanUVeviTKMHl4bp5N7ZFaztrFwHZxUcXi3BD2v2MtcMgh+N6OqKnyyxt8c5Xl4pIusLkW1weCXNDG1LJ2eTdHQuGWS81DUq+Glxd8WtMWA9LHfBsuuRlHRdDHIXUAWisdA5+ZtJ8COc7kjAIwjC/NPM+ACzl7pm6uwKs7syVGO+AADKS1p1iSnIXKBqZhcRfJFOIneF5oRB5mrT1MKxvs4HadNDoBYpMTPALrhyji6G8o4vC3JSl05kOmvdmyoCNLlWiCVA3iE2E7kLMHFNwVyf4yVyl7DIkYBHEISTQ9PVBdSCHmAmUpdy2QqV7LMWZLzUlXZ6LrWf2q4qZ1fpipqVd3S57A+S4mZU4yo0YoFzaGZYK272MtYoZxfQEiAlgRNritke9kXNNkvGWOQEVTKsL3C2HeUdXfXZXWx1zPqwl7dcHx9dzfCy2j+nXHZnlgXOIncJiwEJeARBOLmMaGAIpM6hhtRVK75tSF3aukaGyrm7ppO6SFMMfmx0eJGXuYZt7ZWzC4ncRbWhpdGqThiu9YF7zH6wadrEMNxgyTU0bEhdVSfnIHl5p1dRyV22EfxUHZyzdqmr0cmZdKOTcwh+yEtbIncJiwgJeARBOHXMooFhc15Xm6sL0KiEKQcZ3+m5hJtnRda9jQFCB0IFwCpXVAyQc8n7GV1kmk0E2UthdRlpKJvT9nG9syt2eEbd3cWJ1OUSXM5GX5e83OuCw4s0+SJp8l2ufXNC/waxkaHyaShf9M2UaG7KDjUzjM44RSJ3CYsCCXgEQTi1zKKBIch1boYxw64u7ep/oAistB/F4DJAoZkhWwsqE1eXcVkgm3lnlM/sqDJIXhylLpVIXZRKXWrY2UWpzNUoao73k/3k3V1kAasI0G5ERZC62Gd8Qv+gUODMGcOaIHWpmO0ho2rNDJElc7uMreSu6YqbpxtfIf18hNMQCXgEQVgYjAp8UP2o1gKfZo2PMU6S0drZ2UvlZK20e3MYVloEicvV92gvabGuXF02T+SuOLCU6x2cE6kr1OU0Le2hzidmgpJgJ6Z3gmqXyF1sG92bY1PDytUVB5fmFAOdehNDX+dTuqwOGQbXpK6kzseo2tBSslzJX0HqAqrgp9liABL8CAsbCXgEQVhYnIijC/CZCQAhONJwGZLm+1jvyALA0K4nDfmCXnJW8KqJoRsTEV1dxktMvolhLYvjg5hQr8NtzQtRHT9k9krkruqkPlDgEBg56ctqihKbCu+n/fs3pC6wPzKVusJ6+7aHNalLufepSV3TNDIUhIWOBDyCICw8ZuvoapW5VNVYT2swux9sUs5qRd6WzVbFbAdKN6MryF2k2c3JKl22x01ETySuLEhNLttD02R7QlEzpY6uhsxFqGqBwvNMPrtDACl2gYiGuy6DmtxFhsAlV9etfRG0VkDusz3aNSuMrq5SVVkepUDa3Wdr/YcclrpgbSxuFqlLOF2QgEcQhIXLichcxrQPKU3mdEUJTCug1FHuoqR5YRhU2mpnz5I5Xb6mJ9T5OPmrXuNDiZ09NC1Mg580mRNkLgplNNoFPuxlrlB4bZP5XUHyqtnZfW2PMgwuq9oeyp3sBa29pFV1b46SF7OTukLXZvbDTYnqUpe1gCaRuoQFjQQ8giAsfNpkrvDUKJmrebxlJ3NZ5TITBi7oCX1pyFbSkyVQCSBD4uhSgGFnWoIFk3IBBrn2fkQuRZMWLLtLqK6d/WnA7RIXwj6uZ3woND2MH5qSk/i77LNbtWaGzlXmJDA4MxuTC+zgghTn7PL9eVjV5bXQp0cr1x37GK4ukbqEhYwEPIIgnB7M1M0FHLfMBbIg1i5ToRSINWAYlDHYOlcUWwXKfG1LRrCGQN5VZTN33zm5ksyPD8JCtgcqiQeSbA+AVkt7M+ipyVyW3bkUoBT5LZItQSkGaz+rS/sC55Kcm0v7bE5J7rNaWxU2K+uCI2NBZFxARKOzPWD2tnspbBYWHhLwCIJwejFbG/uxOjYbC07dXNZWbq7S+B98HYeScuYiiabMxZmC8rO52Du5rK7b2VlN4+hKs0KNImhKipFZ+YJoGwIh7+bSicxl3IDS0NCwmtPFrsbHAFQSlGFn07fVnC6U1sl4LY4uJw3amnsrOrrSGp+kgWFa4xP/NkB7A8Pa31UCH2FukYBHEITTk5k2LQz72poWBjcXEwA3VsE1MqxwjQT9+xhyLQCz6WUukLeLexdXGsTEWh0bzl1lb8LjoY/qg54YCPnz1OZ0wQUf8JKZSx5VMldwbEWZCwrOhwag9G+Szumyvu4prF14vq15oUJjbpknzOoKjQ+T8RTSwFA42UjAIwjC6ctsOjUjZER80GMMmMh1X4bxP7rsmhYC3oatAc0gZjBrl8mwoe7HyUmtMpd10pbN4LMvSUGzRWyTbDWgLNVGVIQimqoGCMO1PsG6bpOgh+ALtL0P3vpmhjRC5iqdlMeGoRSDDIEUVU0LNTmZSylXw5PKXDbU9vhePb7WZ0jmCs4vv55OkpNsj3BqkIBHEITTm9nW9qQSF1BJM0r52VwMtuS1qJbBpFr5bs0+SzSdjb0kL3v5njmhU7JywQ+FxoXhcYuNvZbxSTNFcEEUh1ghZIkU+Rld8LZ0/9wIG7sqq/oeyhKZy8tXZNnJXCFb05zNFazsM5W5YuCDoRofCXyE+UQd+5Bhtm/fjo0bN6LX62HTpk144oknpj3+8ccfx6ZNm9Dr9XDhhRfivvvuqz3/wAMP4PLLL8eZZ56JM888E29+85vxjW9843guTRCEpUrzx7BmYef6vvAjnL42zI7yP9LV1mUxMGTXtpUl2/iAwVg3osL4EQ/W34zP9NTuIwYmrsMyx/vpLTQwHKIhkaVT22OX5+R+Wi9Umwjva47iTaVbF3mxdlIXK+WkKXLPu0DRb4m8U01V+9sIgSUwdAylj6nl52mEU08QZsKsA56HH34YN9xwAz7ykY9gz549uPzyy/HWt74VL7zwQuvxzz//PN72trfh8ssvx549e/CHf/iH+PCHP4wvfelL8ZjHHnsM73nPe/Av//Iv2L17N84//3xs3boVL7744vF/MkEQlh7HE/RYG7sIc5BnYoDjHrO19aCnNO5xUbqeNoVp3CzUwEIP3FYV7O6X7G+o3y8wtC9MTE8DJCQB0lAQFIKeWmDjRmPEHkFhEnxGyX3A5n50RkawuS/CzhW4o2A7GpzXb8gzIM/AmXZ9fDINdHIgy0Bhn9JAlrnsmA4T7HUtm0YxYFIxwCFFVeCT7K8+J0ngIxwXxDy7HOEb3vAGXHrppbj33nvjvosvvhhXXXUV7rzzzqHjb7rpJjzyyCN49tln475rr70W3/nOd7B79+7W9zDG4Mwzz8Tdd9+N3/zN35zRdU1MTGDVqlX4RbwdGeWz+UiCICxGmj+KyQ9n7QcVABQlWQeFusSl3Y+yclvn4PJNCr3EBaWATMdGhS5bouJ8rjiTK2lYWM3iam6rLAyryskVszUhsAGGa3uSjFDcJsES2UriUiGQMuybFibb6ORy/XaUz1yRz1yhrLo1x6aFqcwVgkbmqmlh2KaBZvj5CYXNQFXYHDs2p5m4lr4+InOd9pRc4DF8GQcPHsTKlSvn7X1mleEZDAZ4+umnsXXr1tr+rVu34qmnnmp9ze7du4eOv/LKK/HNb34TRVG0vmZychJFUeCss84aeS39fh8TExO1myAIQuREJK7m6/w08ZESl/+RJxv2JXJWvFWylcvacE3aqj3m4X1DwUwbDWmrKXe1S1zUsqWGdZ6qGwHQVEldNOIWA0efwUm2lD4fLn1E1kZkLmGumFXR8v79+2GMwbp162r7161bh3379rW+Zt++fa3Hl2WJ/fv3Y/369UOvufnmm/GqV70Kb37zm0dey5133onbbrttNpcvCMJSYzbW9XQQqX9t7Nmj4QtsjVOSFNWSKwSvMPnGe8FiHrs4s3LudybAEphVZRWP22AhJ5fkUPCdkV2nZArHpo0K02xP+rGpuqaadZ1cwKL8HncfsER+eGqy9TGKG8XBbmiqAsKoDpB3cQVpypJbL1s9hvEZMePXjXwLAHJurtipWWtf6Dxc0Ay2M7Owh/UWhBEcl0urGYlz4niY6fFt+wHg4x//OB566CE89thj6PV6I895yy23YNu2bfHxxMQENmzYMKPrFwRhCXE8Li5F9Z49gO/QTC6A0X5yOzNgyLu4vMRlratXsapq4Jfa1xWBjXNtwTqJizLv9LYE67M6rCg641mxz8xUTq5jBT5DQQ9V8Z8litkjq3zGKe3SrBGHojoJyz9fwlnYLTmJS5ObtaXqEhd763rNwh46NRufslI0bGHX+vidXLW/sQQ+wjCzCnjWrFkDrfVQNufll18eyuIEzjnnnNbjsyzD6tWra/v//M//HHfccQf+6Z/+Ca973eumvZZut4tutzubyxcEQahIsj1t+2bcqNAd7OdHOfmKyf0QExEY1mU6iMHkeuSQdc+R9T1zDPnmhBxncoWGhDGoYcQMUK1oOZQj8TGCHlTHu2DCPxOyRkAt4wTlAxH287GY/YBQVK/xASBpd8HsZ3ONlLis9Q0IFZw/Xg3P5QqXGQKhBsfM9gjCCGZVw9PpdLBp0ybs2rWrtn/Xrl1405ve1PqaLVu2DB3/6KOPYvPmzcjzqrj4E5/4BP70T/8UX/3qV7F58+bZXJYgCMKxaftX/6jCWP+4VlwbCmstx4JcpK6utuLd0ni7ugFKC+VvwbquCl8gXHK8X3NxmVA8jMS1xdV09NTe7mt7Wut70pqeJDOUTnJP7erVBPiGdT11cWW+MDvXceQGa3cfmXNusb+Fx9Aa5BsZQvnCb18UXqvr0boqGlftTq7qs42o65HaHqHBrCWtbdu24eqrr8bmzZuxZcsW3H///XjhhRdw7bXXAnBS04svvojPfvazAJwj6+6778a2bdtwzTXXYPfu3dixYwceeuiheM6Pf/zj+OM//mM8+OCDePWrXx0zQsuXL8fy5cvn4nMKgiAMy1tAzOpU8pav62l2Z/avr83i0hpQ3D6E1Cr3Y8/s5C4LFwR5mctaN3uLrQKFIKPRodlalwUKHZorB1cYF0FV+iZxcIVPFwqWw1iKcGhoaOgySOTv+wLpZD8p9sNH3WuUCXFH2qGZ3WsUOUlLk2tUGGp5jKvtIZNke5hBZFyDR8uxLojY1iWuoS7NmL3EJfKW4Jl1wPPud78bBw4cwO233469e/fikksuwc6dO3HBBRcAAPbu3VvrybNx40bs3LkTN954I+655x6ce+65uOuuu/DOd74zHrN9+3YMBgP82q/9Wu29PvrRj+LWW289zo8mCIIwQ2Yyh6stYxACI01R0oIxVRGusb5Yl0DkZTHjpC3y4ynCKAg3g8sX9IaAo1FIDFsFKuwVKeJQMO2vKdGwgszVFvSEY9m/bZTOlHs9B8lKJbO4fCFQqDEKA02dpAXXhZptlLZYKfc52Tm84jwuCz+N3tcmGVRSV1PiSmo+j0vikqBH8My6D89CRfrwCIIwK9oCmFG9evz9muySSi6knDxDXq4JDfbSXj2+d0+Yth5lIEWz6tXjxlBghLUcqBUxt6k6ITZI+vRU/XoSa7xJt8O9eig+5ijRhY7TqgwNGv2+pDePs+ZX8h+sb+xYs/639OvxWbYoKwK14Eb69Zy+nKw+PDJLSxCEpck08pa7O83wUf/6mm3dwGVD/HNO3mJn6fY/0szehaR8dbA/znp5SrECs4VlcodYt98VUDsHl8uq+EJin2GpWdeDXDUq8Am1yqGg2Wd0wlq413qZyvjYzhBsUlAdymkUIU6DZ6WglEs5WQIoPKHYzWb1waGbxdWwrgMA8bB1nW3NUj8r63r4e6avT//uwpJjVkXLgiAIi57pGhQC7Q0Ka7O3wjbJTnDy2CbBjglb6+dscbINGRcvcTHXuyenW4aXjhBlrvpnajwOv/1pMJQGR6m8RSGIqrJIsQA6doP2wQwlc7jCvC0Fl9Ea1aCwrUlheqnNTFzbLK62YmZBaCAZHkEQljbN5oTTYTnatefUsh7OTa6uRwGwrtAFyh0FBYIFx2245rRIOVTouEwUYgPEUO/TKnEB9doegg9OuAp6wqkR6miqa4Fmf9nkDi1dTY+LsZSvMbKupgeqKqhuuRx3mRYudWVG1/SEye3MPs3EsX5HmhQKo5CARxAEoa0jM4Dg3orSFuB/hKcJemAqWYwqlxFY++CAwdpndlj7YaCqCozY1QO51jft8lbYWu2fi9IWorwVApWavBU+X8jgJNIW/MsRg6Qk6PHZH2V800LjA6/k3FQieax8obWXt0ogFm77HkQoG0GPOoa8FTozh/WcTdCT/k3jGkgx81JDJC1BEARg2h+/WkEsUBXNAu19esL52O+LkpYPdFJZK0pdQaaq5m+l87OQyFa1WVtp4fFM5C0gqdnxDxMZK529xb6mp20uV03aCo4tQqwDinIWuaLsdPYWh4LuIHvNRN5KB7wmiLwlzBTJ8AiCIIziRO3qCvUZXSGrEG4mZJJ0zH64iuBQ6Bu2QK0js8+NhO7MTOSSIUDM3LRZ1mMmh6p9rQSJq1HTE4qFOZyA0/v1jsyA7wodMkapVT3467VfR8A1a1QhlZRY1uPArxZ5y2/Fsi7MBAl4BEEQAm31PDN1brXV89jQZEa7OpbktGnMQeFE4TED1lfyACqauix7g1eGpK4HUMzOup68ntn3z4mZjyojE+t60o/ZkLYoXlJ1peFeeO9QcUNEvtbI1/H4TI/1n8JdmM+8EINLuF8fQ8PyVhhOCreuDPis2YiaHi9zyeBR4ViIpCUIgpByvD98bdJW+HGN2xHSVptzixPHlj8mlavI1+sEp1bq3qrN2hola40ika/C46a01XRu1UZX6Op4qMSxFaQvgssAUSV51W6Al7Pc/ZhBa5W51HCAeszPJz97SxX5ywuCIByL6azqoUHeKJKgptZIr+UWgp3qPmJgFGp2kDQDrAc7LfU8tqWeZ5ranmbWp9Wq3qjliY9D80Py9TwqCXCSIGdkPU9o0BiCmLZ6niQgmnO7uszeWvSIpCUIgtDkGNJW275preozlLYAeKeSBWeAKp205eppGCgRhS6AGpZ1IAhD4Zjqfdz+aFNP3htoD3SqVwGhMHk6aQsI1vTqveDHqkdpK4Nzc2UKqnTSXGJwr7/n0BqFMR5+KU/Urh7+frU3kpqexYxkeARBEGbJcbm23Av9drS0RfE5VNIWY1bSVjObM7Rvus/WVoMdMzkt0pbf3xxpMZ20VR3XkLaatElbcb8anrLuaS0kxwwzPcKiRQIeQRCENo7jX/ojpa3QfRlItqOlLUQpygdFscanCnCG5mDFoIaHAp+aVR3DdT6tNT6NACbdN/RcTeai1iCoOjYEOYmkFbZh30yt6m2o5PmGtDUjRNpatIikJQiCMFOms6n7LswAapbp2XVhZsA4yYutF6fIi0ekhierg518lZE7RwgyUNnlowU9KWxuSltwL2/N7sSP7k+dNiT0o8SqA+CdYf7sYe5XPAfDSVKo9scOzGyQOtUYOD6rug86o1V9tp2YAZG2FimS4REEQRjFXGZ5gFlneeJ9AFUX5XqWp9mQsJnZGcryYERGZ9T+mWR50P44zfIgkbNiQ8KGpDXtvC0g2R7DnaVmk9GRn8GlgvylBUEQZsMxHFtD2EYGwQcdtcAoreXxhOAkdmDmltodNB7HY9AuX6U1PM06HyTHxweNS6f6/fCYhwKdpJYnGXERB4yG4Cc9Zwh82pBaHmEOkIBHEARhOo43y5PIK7W+PM2J6kOZHesa7zEDsVA5Ob42Ub2q5ams66hqefz9tlqeY05UH0Uzc+P78MQREyHISXvvKKqeD8drX5OjqDZNnZWqj52YiU29jdnU8rTtl1qeRYcEPIIgCHNBWy3IyGPTQGjE66KUlchQrVkd1IOVEUXL9XM3tseiTdZKT9eUtlA/PgY5QMzyVK8LTi7Ug4xUykr3YUSWJz6eWTNCyfIsPSTgEQRBmC1tstaJ0rSooy5rhcfN4GVU1oaaQU1bcIR6IDVbaoEOjhH0JDJWeJweWzUVRJXFQfJcM0BpeTztcFFhySMBjyAIwlxyrJ48o8ZNtJ0ndF1u9O0JtvWqMzOGOizPxKI+I1mrrY6HkoAlHJPsTy3ozU7N9cdUBTNNizpQyVpIj00yOHNpURdZa9EjAY8gCMKxmGu3VqDWlLClngeI2Z3aNSQ9dVKouT8tXB55oce+zJGkwUu6Dy1urfh8kumJBc1UexyOG3JptV7DiOckWBEaSMAjCIJwshnq1HyM+p9asFN/qrWOB8NBzjGzOSeozHFL0DNEm6Q1W0bV8RzPqaSOZ0khAY8gCMJccbyFy9Mc07SnA3VZqkns0RPPgfr9GQRHrdCxg5rm87UsT/prkwQas7anj7y+hj09fUqyPQIk4BEEQTg+ZhPczOq800cfNOL5tgBoqHB5LjhG7NDs1jzd45GFy+F+s3A5EAKak5GhkWBp0bBoRkuEfyWVKOb2/9yCIAgA2n/p3Q8v1X7VVfxVJ6DyZLMCxbSHAizBzVkgP0Ii3AhstH+dAliBSYHZP8cEZgUGwYatdVkRywS25Gqlwz7jet+wP2V1c8+zv9BaMXGTJHCqdXBOe//4LVmAjS+WNu4YZQAyDFW6gmtlGVT6fkKlrbahENvauIUxbgp97EfULAQ3Sd2Te46T+yEw5TTzlbrsatmwY7QIEOaFEgWAGda9nQCLJuA5cOAAAOBr2HmKr0QQhEVJ23+Lw755SvYIwlLiwIEDWLVq1bydf9EEPGeddRYA4IUXXpjXBVvsTExMYMOGDfjRj36ElStXnurLOa2RtZw7ZC3nBlnHuUPWcu44ePAgzj///Pg7Pl8smoBHeU131apV8uWbA1auXCnrOEfIWs4dspZzg6zj3CFrOXeo2Qx9PZ7zz+vZBUEQBEEQFgAS8AiCIAiCsOhZNAFPt9vFRz/6UXS73VN9Kac1so5zh6zl3CFrOTfIOs4dspZzx8laS+L59oEJgiAIgiCcYhZNhkcQBEEQBGEUEvAIgiAIgrDokYBHEARBEIRFjwQ8giAIgiAseiTgEQRBEARh0bNgA57t27dj48aN6PV62LRpE5544olpj3/88cexadMm9Ho9XHjhhbjvvvuGjvnSl76E1772teh2u3jta1+Lv//7v5+vy19QzPVafuYznwERDd2mpqbm82Occmazjnv37sV73/tevOY1r4FSCjfccEPrcfKdnJu1XKrfSWB2a/l3f/d3eMtb3oK1a9di5cqV2LJlC/7xH/9x6Lil+L2c63WU7+TM1vJrX/saLrvsMqxevRpjY2O46KKL8Jd/+ZdDx83Jd5IXIJ///Oc5z3N+4IEH+JlnnuHrr7+ely1bxj/84Q9bj//+97/P4+PjfP311/MzzzzDDzzwAOd5zl/84hfjMU899RRrrfmOO+7gZ599lu+44w7Osoy//vWvn6yPdUqYj7X89Kc/zStXruS9e/fWbouZ2a7j888/zx/+8If5b/7mb/jnfu7n+Prrrx86Rr6Tc7eWS/E7yTz7tbz++uv5z/7sz/gb3/gGP/fcc3zLLbdwnuf8rW99Kx6zFL+X87GO8p2c2Vp+61vf4gcffJD//d//nZ9//nn+27/9Wx4fH+dPfepT8Zi5+k4uyIDnF37hF/jaa6+t7bvooov45ptvbj3+D/7gD/iiiy6q7fud3/kdfuMb3xgfv+td7+Jf+qVfqh1z5ZVX8q//+q/P0VUvTOZjLT/96U/zqlWr5vxaFzKzXceUK664ovVHWr6TFSe6lkvxO8l8YmsZeO1rX8u33XZbfLwUv5fzsY7ynayY7Vq+4x3v4Pe9733x8Vx9JxecpDUYDPD0009j69attf1bt27FU0891fqa3bt3Dx1/5ZVX4pvf/CaKopj2mFHnXAzM11oCwOHDh3HBBRfgvPPOw6/8yq9gz549c/8BFgjHs44zQb6TFSe6lsDS+k4Cc7OW1locOnSoNqV6qX0v52sdAflOBmazlnv27MFTTz2FK664Iu6bq+/kggt49u/fD2MM1q1bV9u/bt067Nu3r/U1+/btaz2+LEvs379/2mNGnXMxMF9redFFF+Ezn/kMHnnkETz00EPo9Xq47LLL8L3vfW9+Psgp5njWcSbId7LiRD/3UvtOAnOzln/xF3+BI0eO4F3velfct9S+l/O1jvKdrJjJWp533nnodrvYvHkzrrvuOvz2b/92fG6uvpPZrI4+iRBR7TEzD+071vHN/bM952JhrtfyjW98I974xjfG5y+77DJceuml+OQnP4m77rprri57wTEf3x/5TjpO9HMv1e8kcPxr+dBDD+HWW2/Fl7/8ZZx99tlzcs7TmbleR/lOVsxkLZ944gkcPnwYX//613HzzTfjp3/6p/Ge97znhM7ZZMEFPGvWrIHWeihye/nll4civMA555zTenyWZVi9evW0x4w652JgvtayiVIKP//zP79o/+VyPOs4E+Q7WTHXn3uxfyeBE1vLhx9+GB/84AfxhS98AW9+85trzy217+V8rWMT+U5Ov5YbN24EAPzsz/4sfvzjH+PWW2+NAc9cfScXnKTV6XSwadMm7Nq1q7Z/165deNOb3tT6mi1btgwd/+ijj2Lz5s3I83zaY0adczEwX2vZhJnx7W9/G+vXr5+bC19gHM86zgT5Tlac6Fo2WezfSeD41/Khhx7Cb/3Wb+HBBx/EL//yLw89v9S+l/O1jk3kOznz7w8zo9/vx8dz9p2cVYnzSSLY2nbs2MHPPPMM33DDDbxs2TL+wQ9+wMzMN998M1999dXx+GClvvHGG/mZZ57hHTt2DFmpn3zySdZa88c+9jF+9tln+WMf+9iit1oyz89a3nrrrfzVr36V/+u//ov37NnDH/jABzjLMv7Xf/3Xk/75ThazXUdm5j179vCePXt406ZN/N73vpf37NnD3/3ud+Pz8p2cu7Vcit9J5tmv5YMPPshZlvE999xTs0q/8sor8Zil+L2cj3WU7+TM1vLuu+/mRx55hJ977jl+7rnn+K//+q955cqV/JGPfCQeM1ffyQUZ8DAz33PPPXzBBRdwp9PhSy+9lB9//PH43Pvf/36+4oorasc/9thj/PrXv547nQ6/+tWv5nvvvXfonF/4whf4Na95Ded5zhdddBF/6Utfmu+PsSCY67W84YYb+Pzzz+dOp8Nr167lrVu38lNPPXUyPsopZbbrCGDodsEFF9SOke/k3KzlUv1OMs9uLa+44orWtXz/+99fO+dS/F7O9TrKd3Jma3nXXXfxz/zMz/D4+DivXLmSX//61/P27dvZGFM751x8J4nZV6QKgiAIgiAsUhZcDY8gCIIgCMJcIwGPIAiCIAiLHgl4BEEQBEFY9EjAIwiCIAjCokcCHkEQBEEQFj0S8AiCIAiCsOiRgEcQBEEQhEWPBDyCIAiCICx6JOARBEEQBGHRIwGPIAiCIAiLHgl4BEEQBEFY9Pz/mTAkPDppG2MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "phi_visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "81ba399b-f716-4e59-b8b7-9b20fc74c86b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAHACAYAAABwEmgAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABL4UlEQVR4nO3deXgT5cI28HuyNN1oS/cWSik7iCIWlBY5gGKx4u4RFJVF8IUPkc2Vg0dxeUWOG4ejoEcpuAAigsflRaAeZZFFBIoiVNkKZWmpLdB9Teb7I03atGlJMpNOkrl/15Ur6WTmyZNhoDfPNoIoiiKIiIiIFKJRugJERESkbgwjREREpCiGESIiIlIUwwgREREpimGEiIiIFMUwQkRERIpiGCEiIiJFMYwQERGRohhGiIiISFEMI0RERKQorwoj27Ztw2233Yb4+HgIgoD//Oc/in/e+vXrMXLkSERGRkIQBBw4cMCtdSIiIvI1XhVGysvL0a9fP7z99tse83nl5eUYPHgwXn311TapExERka/RKV0BZ6SnpyM9Pb3F92tqavDss89i5cqVuHTpEvr27YuFCxdi2LBhbvk8AHjooYcAACdPnnTpM4iIiNTOq8LI5UycOBEnT57Ep59+ivj4eHzxxRe4+eabcfDgQXTv3l3p6hEREZEdXtVN05rjx49j9erVWLt2LYYMGYKuXbviiSeewPXXX4/ly5crXT0iIiJqgc+Ekf3790MURfTo0QPBwcHWx9atW3H8+HEA5q4UQRBafUyfPl3hb0JERKQuPtNNYzKZoNVqsW/fPmi1Wpv3goODAQAdOnRAdnZ2q+W0b9/ebXUkIiKi5nwmjPTv3x9GoxEFBQUYMmSI3X30ej169erVxjUjIiKi1nhVGCkrK8OxY8esP+fk5ODAgQMIDw9Hjx498MADD2DcuHF444030L9/fxQWFuL777/HlVdeiVtuuUXWz+vUqRMA4MKFC8jNzcW5c+cAAH/88QcAIDY2FrGxsVK+LhERkTqIXuSHH34QATR7jB8/XhRFUaypqRGfe+45sXPnzqJerxdjY2PFu+66S/z111/d8nmiKIrLly+3u8/zzz8v/QsTERGpgCCKoqhABiIiIiIC4EOzaYiIiMg7MYwQERGRorxiAKvJZMK5c+fQrl07CIKgdHWIiIjIAaIoorS0FPHx8dBoWm7/8Iowcu7cOSQkJChdDSIiInLB6dOn0bFjxxbf94ow0q5dOwDmLxMSEqJwbYiIiMgRJSUlSEhIsP4eb4lXhBFL10xISAjDCBERkZe53BALDmAlIiIiRTGMEBERkaIYRoiIiEhRXjFmhIiIqC0YjUbU1tYqXQ2vodfrodVqJZfDMEJERKoniiLy8/Nx6dIlpavidcLCwhAbGytpHTCGESIiUj1LEImOjkZgYCAX2HSAKIqoqKhAQUEBACAuLs7lshhGiIhI1YxGozWIREREKF0drxIQEAAAKCgoQHR0tMtdNhzASkREqmYZIxIYGKhwTbyT5bxJGWvDMEJERITLL8xF9slx3hhGiIiISFEMI0RERKQohhEiIiIvNWHCBAiC0OyRnp6OESNG2D1m165dEAQB+/fvx8mTJyEIAnQ6Hc6ePWuzX15eHnQ6HQRBwMmTJ936PTibRiaFZdWoqjVKKsNPq0F0iL9MNSIiIjW4+eabsXz5cpttO3bswL333otTp04hMTHR5r2MjAxcffXVuOaaa6whIz4+Hh999BHmzp1r3e/DDz9Ehw4dkJub6/bvwDAig7V7T+PJz3+VpawnR/bEo8O7yVIWERH5PoPBgNjYWJttd9xxB6Kjo7FixQo8//zz1u0VFRVYs2YNXnnlFZv9x48fj+XLl9uEkRUrVmD8+PF46aWX3PsFwG4aWfx6phgAoNUIMOg0Lj10GvNo5IP1ZRERkXJEUURFbYUiD1EUJddfp9Nh3LhxWLFihU15a9euRU1NDR544AGb/W+//XZcvHgRP/74IwDgxx9/xIULF3DbbbdJrotD9W2TT/FxIsx/0NOHd8Psm3q4VMYnu0/h2f/8Zi2LiIiUU1lXietWXafIZ/809icE6h1f8+Sbb75BcHCw9ef09HSsXbsWDz/8MF577TVs2bIFw4cPB2Duorn77rvRvn17mzL0ej0efPBBZGRk4Prrr0dGRgYefPBB6PV6eb7UZTCMyECGEOuWsoiIyPcNHz4cS5cutf4cFBQEAOjVqxdSU1ORkZGB4cOH4/jx49i+fTs2b95st5xJkyYhJSUFr7zyCtauXYtdu3ahrq6uTb4Dw4iMpKz7wrV2iIg8R4AuAD+N/Umxz3ZGUFAQunWzP9Zw0qRJmD59Ot555x0sX74ciYmJuPHGG+3u27dvX/Tq1Qv3338/evfujb59++LAgQPOVt8lDCMysDRmCHA9UViOZcMIEZHyBEFwqqvEU40ePRozZ87EqlWr8OGHH+KRRx5pdcXUhx9+GNOmTbNpaWkLDCMyYDcNERF5ouDgYIwZMwZ/+9vfUFxcjAkTJrS6/yOPPIJ7770XYWFhbVI/C86mkRG7aYiIyNNMmjQJFy9exIgRI9CpU6dW99XpdIiMjIRO17ZtFWwZkYW5OUNKnmg4lk0jRETkmBUrVlx2n5SUlBanC3fu3LnVqcRXX321LFONL4ctIzJgNw0REZHrGEZkxG4aIiIi5zGMyMDSmtHaCOXL4WwaIiJSK4YRGci5ampb9M0RERF5EoYRT1HfqMIoQkSkDP5n0DVynDeGERk0dNO4XgaHjBARKcNy/5WKigqFa+KdLOdNyn1sOLVXBnJmaQZzIqK2pdVqERYWhoKCAgBAYGCgpDGAaiGKIioqKlBQUICwsDBotVqXy3I6jGzbtg2vvfYa9u3bh7y8PHzxxRe48847Wz1m69atmDNnDg4dOoT4+Hg89dRTmDp1qqt19jjWlhEpy8ELHMBKRKSU2NhYALAGEnJcWFiY9fy5yukwUl5ejn79+mHixIm45557Lrt/Tk4ObrnlFjzyyCP45JNPsGPHDkybNg1RUVEOHe9N2E1DROSdBEFAXFwcoqOjUVtbq3R1vIZer5fUImLhdBhJT09Henq6w/u/++676NSpExYtWgQA6N27N/bu3YvXX3/dZ8KIKMcKrJYBrOynISJSjFarleWXKznH7QNYd+3ahbS0NJttI0eOxN69e1tMn9XV1SgpKbF5eDTmByIiIpe5PYzk5+cjJibGZltMTAzq6upQWFho95gFCxYgNDTU+khISHB3NWXBFViJiIic1yZTe5uOSrZ0RbQ0Wnnu3LkoLi62Pk6fPu32OkphaRiRNIDVsgIrW1mIiEhl3D61NzY2Fvn5+TbbCgoKoNPpEBERYfcYg8EAg8Hg7qrJRs5xHnKu5kpEROQN3N4ykpKSgszMTJttmzdvxoABAyQtkOKJ2E1DRETkPKfDSFlZGQ4cOIADBw4AME/dPXDgAHJzcwGYu1jGjRtn3X/q1Kk4deoU5syZg+zsbGRkZGDZsmV44okn5PkGHoCLnhEREbnO6W6avXv3Yvjw4daf58yZAwAYP348VqxYgby8PGswAYCkpCRs2LABs2fPxjvvvIP4+HgsXrzYZ6b1AvIGCIYRIiJSG6fDyLBhw1odI7FixYpm24YOHYr9+/c7+1FeR8rywVx6mIiI1Io3ypNBw2wa11mO5QBWIiJSG4YRGcg6m4ZZhIiIVIZhREZyzKZhFiEiIrVhGJGBPN00HDNCRETqxDAih/o0Im0Aq21ZREREasEwIgM5B51yACsREakNw4iMJI0Zka8aREREXoVhRAaWGTCSxoxYBrCyYYSIiFSGYUQGsq7AKl9RREREXoFhRE6SVlFlRw0REakTw4gMLINO5emmYdsIERGpC8OIDNhNQ0RE5DqGERlxNg0REZHzGEZk0LACq/S79rKXhoiI1IZhRAbspiEiInIdw4iMZOmmYdMIERGpDMOILOSbTUNERKQ2DCMyYDcNERGR6xhGZCSpm4bLwRMRkUoxjMhAltk0nNxLREQqxTAiA1GWO+XVl8WOGiIiUhmGERnIGR/YTUNERGrDMCIj3iaPiIjIeQwjMrD20kgYwcoVWImISK0YRmQgazeNjGURERF5A4YRGbGbhoiIyHkMIzKwzKaRZ50Rto0QEZG6MIwQERGRohhGZCTtRnnsqCEiInViGJFBw5pnUmbT2JZFRESkFgwjMpBz1VSuwEpERGrDMCIjad00ZmwZISIitWEYkYEsAYJDRoiISKUYRmQgywqs9WmEDSNERKQ2DCMehuuMEBGR2jCMyMAy6FTSCqzspiEiIpViGJFBQzeN62VYB7BKrg0REZF3YRiRgawBgmmEiIhUhmFERtIWPWM/DRERqRPDiBzk6KaxrMAqvTZERERehWFEBrKuwMrZNEREpDIMIzKSNJtGtloQERF5F4YRGcgym4bdNEREpFIMIzKQM0Cwl4aIiNSGYURW0jtqeNdeIiJSG4YRGVgGncrRTUNERKQ2DCMyYDcNERGR6xhGZCTHbBqGESIiUhuGERk0zKbhCqxERETOYhiRgaUxg+uMEBEROY9hxMNwBVYiIlIbhhE5cDYNERGRy1wKI0uWLEFSUhL8/f2RnJyM7du3t7r/ypUr0a9fPwQGBiIuLg4TJ05EUVGRSxX2RNZuGilhxLrOCBERkbo4HUbWrFmDWbNmYd68ecjKysKQIUOQnp6O3Nxcu/v/+OOPGDduHCZNmoRDhw5h7dq1+PnnnzF58mTJlfdF7KUhIiK1cTqMvPnmm5g0aRImT56M3r17Y9GiRUhISMDSpUvt7r9792507twZM2bMQFJSEq6//npMmTIFe/fulVx5T2GdTSNhGCq7aYiISK2cCiM1NTXYt28f0tLSbLanpaVh586ddo9JTU3FmTNnsGHDBoiiiPPnz+Pzzz/HqFGjWvyc6upqlJSU2Dw8mXUJdxkCBZeDJyIitXEqjBQWFsJoNCImJsZme0xMDPLz8+0ek5qaipUrV2LMmDHw8/NDbGwswsLC8K9//avFz1mwYAFCQ0Otj4SEBGeq2ebk7FphNw0REamNSwNYmy7QJYpii4t2HT58GDNmzMBzzz2Hffv2YePGjcjJycHUqVNbLH/u3LkoLi62Pk6fPu1KNducpHVG2E1DREQqpXNm58jISGi12matIAUFBc1aSywWLFiAwYMH48knnwQAXHXVVQgKCsKQIUPw8ssvIy4urtkxBoMBBoPBmaopSpYVWDmbhoiIVMqplhE/Pz8kJycjMzPTZntmZiZSU1PtHlNRUQGNxvZjtFotAN9Z4Is3yiMiInKd0900c+bMwQcffICMjAxkZ2dj9uzZyM3NtXa7zJ07F+PGjbPuf9ttt2H9+vVYunQpTpw4gR07dmDGjBm49tprER8fL9838QDydNMwjRARkbo41U0DAGPGjEFRURFefPFF5OXloW/fvtiwYQMSExMBAHl5eTZrjkyYMAGlpaV4++238fjjjyMsLAw33HADFi5cKN+3UJjIFViJiIhc5nQYAYBp06Zh2rRpdt9bsWJFs22PPfYYHnvsMVc+yqtIWmfEMmaEDSNERKQyvDeNh2EWISIitWEYkUHDbBrXy2A3DRERqRXDiAwsq6ZKGsBqKYv9NEREpDIMIx6GUYSIiNSGYUQG1sYMdtMQERE5jWFEBg1ZRHpHDXtpiIhIbRhGPAzHjBARkdowjMiAi54RERG5jmFEBjIMGWmYTSOxLkRERN6GYUQOst4pT8ayiIiIvADDiIwECX0tlmOZRYiISG0YRmRg7aaRMmZElpoQERF5H4YRGcg5A4azaYiISG0YRmQkaQBr/cGMIkREpDYMIzKQp5uGHTVERKRODCMyaOhZkTKAtWlZRERE6sAw4mFEdtQQEZHKMIzIwBIguIoqERGR8xhGZGDpWpFlACsbRoiISGUYRjwMswgREakNw4gMrC0jMqzASkREpDYMIzKS40Z5bBohIiK1YRjxMJxNQ0REasMwIgPLEu6SFj1jLw0REakUw4gMrCuwSln0rP5YzqYhIiK1YRiRgZwBglmEiIjUhmFERnJ00/CuvUREpDYMIzKQY9Aph4wQEZFaMYzIoGGdEQmFWFpGJNeGiIjIuzCMeBj20hARkdowjMhAztk0REREasMwIgM5umm4zggREakVw4gH4owaIiJSE4YRWciwAqtMNSEiIvI2DCMysHbTSBkz0ijJsGGEiIjUhGHEAzGLEBGRmjCMyMA6m4bdNERERE5jGJGB9a69EspoHGQ4gJWIiNSEYcQDMYoQEZGaMIzIQJ5uGg5gJSIidWIYkUFDeOCgESIiImcxjMhA7jEectwFmIiIyFswjMhIruXg2U1DRERqwjAig4Yb5bmOvTRERKRWDCNysN4oT54VWImIiNSEYcQDsZuGiIjUhGFEBuymISIich3DiAysK7DKNYCVs2mIiEhFGEY8ELtpiIhITRhGZNDQTSNhACs7aoiISKUYRmQgWmfTuF6GbTcNERGRejCMeCDetZeIiNTEpTCyZMkSJCUlwd/fH8nJydi+fXur+1dXV2PevHlITEyEwWBA165dkZGR4VKFPREHnBIREblO5+wBa9aswaxZs7BkyRIMHjwY7733HtLT03H48GF06tTJ7jGjR4/G+fPnsWzZMnTr1g0FBQWoq6uTXHlPwW4aIiIi1zkdRt58801MmjQJkydPBgAsWrQImzZtwtKlS7FgwYJm+2/cuBFbt27FiRMnEB4eDgDo3LmztFr7OPbSEBGRmjjVTVNTU4N9+/YhLS3NZntaWhp27txp95ivvvoKAwYMwD/+8Q906NABPXr0wBNPPIHKykrXa+1hrLNppCwHDzaNEBGROjnVMlJYWAij0YiYmBib7TExMcjPz7d7zIkTJ/Djjz/C398fX3zxBQoLCzFt2jRcuHChxXEj1dXVqK6utv5cUlLiTDXbnqWbRkIRvDUNERGplUsDWJu2AIii2GKrgMlkgiAIWLlyJa699lrccsstePPNN7FixYoWW0cWLFiA0NBQ6yMhIcGVarY5SWNGGr3mgFgiIlITp8JIZGQktFpts1aQgoKCZq0lFnFxcejQoQNCQ0Ot23r37g1RFHHmzBm7x8ydOxfFxcXWx+nTp52pZpuTOzxwzAgREamJU2HEz88PycnJyMzMtNmemZmJ1NRUu8cMHjwY586dQ1lZmXXbkSNHoNFo0LFjR7vHGAwGhISE2Dw8mXU2jZQVWNlPQ0REKuV0N82cOXPwwQcfICMjA9nZ2Zg9ezZyc3MxdepUAOZWjXHjxln3Hzt2LCIiIjBx4kQcPnwY27Ztw5NPPomHH34YAQEB8n0TBTUMYHW9DNtuGiIiIvVwemrvmDFjUFRUhBdffBF5eXno27cvNmzYgMTERABAXl4ecnNzrfsHBwcjMzMTjz32GAYMGICIiAiMHj0aL7/8snzfwsdwBVYiIlITp8MIAEybNg3Tpk2z+96KFSuabevVq1ezrh1fYgkPnE1DRETkPN6bRgbWdgxJK7A2HMx2ESIiUhOGEQ/EXhoiIlIThhEZyDGbhoiISK0YRmQkddyH5XguekZERGrCMOKJmEWIiEhFGEYkajwNV2onjeV4ZhEiIlIThhGJGg82lbqKKldhJSIiNWIY8UCcTUNERGrCMCJR49wgXzcN0wgREakHw4hENmNGZJpNQ0REpCYMIzKSus6I5Xh20xARkZowjEjkjtzALEJERGrCMCKR6I5BI0RERCrCMCJR48GmkseMWMpkPw0REakIw4gHYhYhIiI1YRiRyGbRM4llcTYNERGpEcOIjCSvwMpBI0REpEIMIx6I3TRERKQmDCMSsZuGiIhIGoYRidwym4YrjRARkYowjHggdtMQEZGaMIxIZNtNI3EAa33TCrMIERGpCcOIRDYLsMrUTUNERKQmDCOepD6NcAVWIiJSE4YRidwRHBhFiIhITRhGJGI3DRERkTQMIzKSbQArm0aIiEhFGEYkck9wYBohIiL1YBiRqvHUXqndNOynISIiFWIYkchmBVaJZVlXYGXDCBERqQjDiAdiFiEiIjVhGJHIZgVWif0sUo8nIiLyRgwjEtlM7ZVYFrtpiIhIjRhGPBDv2ktERGrCMCJR4xVY5ZpNw5YRIiJSE4YRiWxXYJWro4aIiEg9GEY8EFtGiIhITRhGJJIzOFi7aThmhIiIVIRhRCJLcJBjVi47aYiISI0YRmQiR5DgAFYiIlIjhhGpGByIiIgkYRiRyJJF5Fg9VWBHDRERqRDDiEzYTUNEROQahhGJ3BEcOJuGiIjUhGFEIs6mISIikoZhRCJLy4gc4z0s407YTUNERGrCMOKBmEWIiEhNGEYksgYH9rEQERG5hGFEIstde+WdTcO2ESIiUg+GEQ/EKEJERGrCMCKRdQCrHLNpuM4IERGpEMOITGSZTcOBJ0REpEIMIzKRs2WEHTVERKQmLoWRJUuWICkpCf7+/khOTsb27dsdOm7Hjh3Q6XS4+uqrXflYj+SWFViZRYiISEWcDiNr1qzBrFmzMG/ePGRlZWHIkCFIT09Hbm5uq8cVFxdj3LhxuPHGG12urCeyrsAqQ1nspCEiIjVyOoy8+eabmDRpEiZPnozevXtj0aJFSEhIwNKlS1s9bsqUKRg7dixSUlJcrqwnk+WuvZYVWCWXRERE5D2cCiM1NTXYt28f0tLSbLanpaVh586dLR63fPlyHD9+HM8//7xDn1NdXY2SkhKbh6diNw0REZE0ToWRwsJCGI1GxMTE2GyPiYlBfn6+3WOOHj2KZ555BitXroROp3PocxYsWIDQ0FDrIyEhwZlqtilLbmA3DRERkWtcGsDatEtCFEW73RRGoxFjx47FCy+8gB49ejhc/ty5c1FcXGx9nD592pVqti0Z0whXYCUiIjVxrKmiXmRkJLRabbNWkIKCgmatJQBQWlqKvXv3IisrC9OnTwcAmEwmiKIInU6HzZs344Ybbmh2nMFggMFgcKZqinFHcGAUISIiNXGqZcTPzw/JycnIzMy02Z6ZmYnU1NRm+4eEhODgwYM4cOCA9TF16lT07NkTBw4cwHXXXSet9h6A3TRERETSONUyAgBz5szBQw89hAEDBiAlJQX//ve/kZubi6lTpwIwd7GcPXsWH330ETQaDfr27WtzfHR0NPz9/Ztt93ayzqZh0wgREamI02FkzJgxKCoqwosvvoi8vDz07dsXGzZsQGJiIgAgLy/vsmuO+BK3zKZhRw0REamIIHrBaMmSkhKEhoaiuLgYISEhSlfHxrGCUox4cxvCAvU48Fza5Q9oxU1vbsXRgjKsmnwdUrtFylRDIiIiZTj6+5v3ppHIetdeGcqS4/42RERE3oZhxAN5fFMVERGRjBhGJLLOppFjACs4gJWIiNSHYUQidtMQERFJwzAiEzmDBGfTEBGRmjCMSOSO4MBuGiIiUhOGEYkagoN8i54RERGpCcOITOTIEZYi2DBCRERqwjAikVtWYGU/DRERqQjDiESWMSNyzaaJRDEHjRARkaowjMhEjm6a26q+xl7//4c+Ox4D6qqlF0hEROQFnL5RHtmSrRGj4gIerPwEABB9ehPw0Z1Av/uAsvNAdUmjHetTjzX9WH7WAFo9oNE1PHQGIKA9EBgOtIsDIrqZtxEREXkQhhGZCFI7ag6tR7BYjjNiJOJ05dDm7gRyd8pTOQtBC0T3AbqPAPreA8ReKW/5RERELmAYkYnkbprzhwAAXxlTMXDkBAw8vxYoLwRC4gD/UJhbQOqbYew1x5iMgKkOMNUCxjrz67oqoPIiUFEEFJ8GqoqB8wfNjx/fArreAIyYD8T1k1h5IiIi1zGMSCRbN01BNgDgiKkjerbvC6TeKFPB9UQRKM0DTu4Afv8ayP4aOP49kLMduOFZIHUGoOEQIiIianv87SORLLNpRNEaRo6KHd0zmUYQgJB44Kp7gdEfATOygN63mVtSvnseWDcJMNa64YOJiIhaxzAikfVGeVL6acrOA1WXYIQGx8R4eSp2Oe07A6M/Bm5bDGj0wKH1wJfTOa2YiIjaHMOIJ6hvFcnXxqIafm23AqsgAMnjgTGfmAe3/vopsOOfbfXpREREABhGJJMlOPz5OwDgtDbRXGZbt070vBm45TXz6/++CJzZ17afT0REqsYwIpElOEiaTVPfMnJa10mGGrlowMPAFXcDohH4ajrHjxARUZthGJGJpDBiaRnR1beMyFAfpwkCcMvrQGAEUHAY2P+RErUgIiIVYhiRSHJwEEWgoEkYUWoMaVAEMPQZ8+utC4GacoUqQkREasIwIpF1No2rk3tL84DqYkDQ4pyuo3wVc1XyBPNMm7LzwO6lSteGiIhUgGFEJi5309SPF0F4FxgFv/qNCk6v1fkBw581v965mK0jRETkdgwjkkkMDvXjRRDdq6FEpZf66HuPuXWkqhj4dY3ClSEiIl/HMCJRQzeNiywtI1G9pd/fRi4aDTDwEfPrn/7tAemIiIh8GcOITFxegbVRy4hl3IlH/Orv/yCgDwT+zAZOble6NkRE5MMYRiSSFBwazaRBVG+bzYoLCAOuGmN+/fMHilaFiIh8G8OIRJK6aYrPADWlgEYHRHSzFiJ6RtsIMGCi+fmPjebxI0RERG7AMCIXV9KIpYsmohug85N25193iL0KiOwJGKuB7G+Urg0REfkohhGJJN1Hxjp4tZfNZo/opgHM85WvvNf8+rfPla0LERH5LIYRiSy5waVWDevgVfN4EcHaTeNBrrzH/HxiC1BWoGhViIjINzGMSGQdM+LKbJomLSMur+LqTuFdgA4DANEEHPpC6doQEZEPYhiRidMxwmQC/vzD/Lppy4jH9NPUu+Iu8/Pv/6dsPYiIyCcxjEjk8syX4tNAbTmg0ZtbHzxZz3Tz86kdnFVDRESyYxiRytpN4+RxlvEikd0Brd61MtpKRFcgsgdgqgOOfad0bYiIyMcwjMjE6fEedmbSWFdg9bBeGgBAj5vNz39sVLYeRETkcxhGJHI5NzSZSWNbpgemEUtXzdHNgLFO2boQEZFPYRiRSHS1m8Zey4indtMAQMdrgYD2QNUl4PRPSteGiIh8CMOIEkwmoPCI+bW9lhEPbBiBVgd0G2F+ffx7ZetCREQ+hWFEIpe6VC6dAmorAK0f0D6peZmeGEYAoMtw8/OJH5StBxER+RSGEYlcWvTMOpOmh7nFAS6UoYQuw8zP57KAyouKVoWIiHwHw4hMnIoRBYfNz026aCxleGrDCEI7mAOUaAJO/qh0bYiIyEcwjEjkUnAoqG8ZaXKDPGuZHttPg4bWkePsqiEiInkwjEhkCQ5O9bD8WT+TpmnLiCfeKK8p67iRLYpWg4iIfAfDiEwcDiMmI1B41Py6ScuIh48YMes8GBC0wIXjwKVcpWtDREQ+gGFEIksrhsMrsF48CdRVATp/oH1nm7cEb2ga8Q8F4q82vz61U9GqEBGRb2AYkcrZ4GBZ7CyyB6DRtlCkJ6cRAImDzc+ndihbDyIi8gkMIxJZgoPD3TQtjBcBvKSbBmgURtgyQkRE0jGMyMThINHKTBprL42HN4yg03UABKDoGFB6XunaEBGRl2MYkcjp4NDKDfKsZbpenbYR0B6I6Wt+za4aIiKSiGFEImsYcaSfxljXcE8au2uMeE1HDZCYan5mVw0REUnkUhhZsmQJkpKS4O/vj+TkZGzfvr3FfdevX4+bbroJUVFRCAkJQUpKCjZt2uRyhT2VQzHiYg5grAH0gUBYYvMyvKWbBjBP8QUYRoiISDKnw8iaNWswa9YszJs3D1lZWRgyZAjS09ORm2t/zYlt27bhpptuwoYNG7Bv3z4MHz4ct912G7KysiRX3hM4lRtsZtK0fOo9fjYNAHSqbxkpOARUXFC2LkRE5NWcDiNvvvkmJk2ahMmTJ6N3795YtGgREhISsHTpUrv7L1q0CE899RQGDhyI7t2745VXXkH37t3x9ddfS668J3BqBVbreJE+dt/2ok4aIDjKHKoAIHe3snUhIiKv5lQYqampwb59+5CWlmazPS0tDTt3OtZcbzKZUFpaivDwcGc+2uM5FCSsN8izf08ar+qmARqNG+EgViIicp1TYaSwsBBGoxExMTE222NiYpCfn+9QGW+88QbKy8sxevToFveprq5GSUmJzcNTOddNY5nW2/JMGqfLVBIXPyMiIhm4NIBVaNInIYpis232rF69GvPnz8eaNWsQHR3d4n4LFixAaGio9ZGQkOBKNduEpRXjst/fWGtelwNouWUEXtY0YmkZyfsVqC5Tti5EROS1nAojkZGR0Gq1zVpBCgoKmrWWNLVmzRpMmjQJn332GUaMGNHqvnPnzkVxcbH1cfr0aWeqqYjLRrGi44CpFvALBkLthyun7vzrCUI7mr+LaATO/Kx0bYiIyEs5FUb8/PyQnJyMzMxMm+2ZmZlITU1t8bjVq1djwoQJWLVqFUaNGnXZzzEYDAgJCbF5eC4HWzEsy8BH9bxs6vCSdhGzToPMzxzESkRELnK6m2bOnDn44IMPkJGRgezsbMyePRu5ubmYOnUqAHOrxrhx46z7r169GuPGjcMbb7yBQYMGIT8/H/n5+SguLpbvWyiooZvmMjs6MF7E6wawAkCnFPNz7i5l60FERF5L5+wBY8aMQVFREV588UXk5eWhb9++2LBhAxITzYt45eXl2aw58t5776Gurg6PPvooHn30Uev28ePHY8WKFdK/gYcQLtdRY71Bnv3xIg6V4YksYeTMz+ZxMVq9svUhIiKv43QYAYBp06Zh2rRpdt9rGjC2bNniykd4DWsjhqMtI63ck6Zh/KoXNY1E9QL8w4CqS0D+r0CHZKVrREREXob3ppHIodxQV90wk+Yy03oBLxszotFw3AgREUnCMCKRZen2VhtGio6ZZ5wYQoCQ+BZ388JOGjNLGOF9aoiIyAUMIzJpdQCr5Z40Ub1a3dGyVok39dIAaDSIdbcXVp6IiJTGMCKRQ797/3RgvEjjMl2vjjLi+wNaA1BRaF5PhYiIyAkMIxJZgkOrM2EsLSOXCSNe202jMzQMXM1lVw0RETmHYUQmDnfTOFCGV82mseAgViIichHDiESXDQ61VcDFHPPr6D7ur5BSuPgZERG5iGFEJi22jBQeAUQTENAeCG755oCAF3fTAEDCtQAE4MIJoPS80rUhIiIvwjAikxbHjFi7aHpfds14r51NAwABYUDMFebXbB0hIiInMIxIdNng8Kdjg1dtyvS++TRmHDdCREQuYBiRyLroWUuNHo4sA19PACDoi7Dn4hr8ceEPeSrYlqzjRjijhoiIHMcw4m4Fh83Pl5lJAwBl4lkEdl6Cny6txr1f34uNORvdXDmZWcJI/kGgulTZuhARkddgGJHI0k0j2GsaqSkHLp0yv3agZeRI7UpodOXmciHi+Z3PI6c4x2afwspCHPzzIE4Un0BxdbFnTQMO7QCEdjIP2D3zs9K1ISIiL+HSXXupQatZ4M/6rpagKCAostVyfvnzFxSZDkIUNXigw9s4alqOn/N/xpwtc/DaX17DkYtHsO7oOuzJ32NznJ/GD+382sEkmmAUjTCJJogQYRJN9fWTHlbsBS3LgF3Le9afw7VAaEdgxxPAHv/W92006LfpZzQ9xu5xguP7Xu47CBAQoAtAoD4QgbpA6+tQQyhiAmMQHRiNmMAYJIUmob1/e7vniYiIXMMwIhO7Q0Ysy8A70EWz8vBKAEBdcX+EdY7HwmsX4q9f/xXHLh3DXV/d1ehzBEQHRqOitgKltaWoMdWgqKpIhm8gI60GEGuBmlqla+IWkQGR6BneEwNjBuLa2GvRO6I3dBr+VSIichX/BZXIuhy8vTRiGS9ymS6awspCZOZmAgBqLgyGCBFRgVFYfvNyPPvjszhYeBCxQbG4q9tduKvbXYgLjgMAVBurUVhZiLKaMmgFLTQaDbSCFgIE6//27XYfOfP96ltWbGb4iLDZ1vg9segYsHI0RJ0/8D9bIGq1LR/XqNWmaVn2WnSavmfzuQ2FN3+vyf723jOJJlTVVaGirgKVdZWoqK1AeV05LlZdxPmK8yioKEBeWR7OlZ9DYWUhCs8WYsfZHQCAcP9w3Nz5Ztza5Vb0jewr+ZwTEakNw4hErXaDODiT5vMjn6POVIcwTXeUVsdbu366hHbBqlGrUGuqhU7QNfslZ9Aa0CG4g5Tqyy84AdC3AyovAhXFQMcBStdIVuW15Th26Rh+K/wNe/L24OfzP+NC1QWs+n0VVv2+Cn0i+uB/rvofDE8YDo3AIVlERI5gGJFJ6900LYeROlMd1h5ZCwDopB+B03b20Wv0kuvXZjQaIGEQcORb8+JnPhZGgvRB6BfVD/2i+uGB3g+gzlSHXed24ZsT3+D73O9xuOgwZv0wC93CuuHJAU8itUOq0lUmIvJ4/K+bRA3dNE3iSFUJUFwfLaJbHjOy5fQWFFQUINw/HHHa69xSxzanosXPdBodhnQcgoV/WYjNf92MR658BMH6YBy7dAxTvpuCOVvmoKjSw8b0EBF5GIYRqVrqpbHMpAmONd+XpgWf/vEpAODu7ndDK5hbQDxquq4rEutbA3J3eena9q5p798eM66ZgU1/3YSH+jwEraBF5qlM/PXrv2LnWS4ER0TUEoYRiawrsDZ9w4Fl4E8Un8BPeT9BI2gwusfoy926xnvEXQ3o/IGKooZQpiIhfiF4auBTWHPrGnQL64bCykJM+W4K3tr3lnXKNRERNWAYkUmzIOHA4NU1v68BAAztOBRxwXHWMry+MUHnV38XXwA525Sti4J6hvfE6lGrMabnGABAxm8ZeHzL46isq1S4ZkREnoVhRKIWg8NlloGvqK3AV8e/AgDc1/M+2zLlqpySugwzP+dsVbQaSvPX+ePZQc9i4ZCF0Gv0+C73O0zaNInjSIiIGmEYkaghODRpGrHMpInuY/e4b058g7LaMiSGJGJQ/CCbMry+ZQRoFEa2A8Y6RaviCW7pcgveT3sfoYZQHCw8iImbJqKwslDpahEReQSGEZnYdNNUXgJK88yvo3o221cURaz5w9xFM7rHaOt6FD4zZgQwjxvxDwWqi4G8A0rXxiMkxyTjk/RPEBsUi5ziHDy86WEGEiIiMIxIZrcVw9IqEtIR8A9p9nZWQRaOXDwCf60/7uh2R/MyfaGjRqMFkv5ifn3iB2Xr4kE6h3ZGRloGAwkRUSMMIxLZnU1zmWXgP/3dPJ13VJdRCDWEWrdbyvCJbhqgoavmhLrHjTSVEJJgE0imfTcN5bXlSleLiEgxDCMyselisc6kaT54tfF9aCyzLOyW4Qu6DDc/n/4JqC5Tti4eJiEkAR+kfYBw/3BkX8jGrB9modbomzcWJCK6HIYRiSytGELjthHLGiN2loFfd2Qd6kx16BfVD70jbN+3lOErDSMI7wK0TwKMNcCJLUrXxuMkhiRiyY1LEKALwO683Zi3Yx7XISEiVWIYkchucGihZaTxfWiatorYFuojcUQQgB4jza+PblK2Lh7qisgrsGjYIugEHb7N+Rbv/vKu0lUiImpzDCMysXaxlBcB5QXm103WGNl6eivOV5xHuH84RnYe2XIZ9USjEWJNjRtq24YsYeTIZsDE//Xbk9ohFc+lPAcAWPrLUnyb863CNSIialsMI1LVt2JYg4SliyYsEfALstl19R+rAZjvQ+On9WtWVMMAVhFFy5bhaOpg/JE8AOeefRam6mp31N79EgcDfsFAWT6Q/4vStfFYd3W/CxOumAAA+PuOv+PgnweVrRARURtiGJGoWYdKgf170jS+D829Pe5ttczu/7cKBa+9DmNxMcTaWhR/vg5nZ86CaDTKV/G2ojMAXesHsh5hV01rZl0zC0M7DkW1sRozfpiB/PJ8patERNQmGEZkYh3A2kIYWZW9CgDwl45/QXxwvP0yBAGD8n5D983mcSXRTz2FjkuXQPD3R9mWLShc2vp4AlNVFUq//wGF7/0bhe++h9L//hd1hR6whkX3+q6a3/9P2Xp4OK1Gi4V/WWi9ud6M72egorZC6WoREbmdTukKeDvrWFNrN0394NVGM2nOl5/H+qPrAQDj+oxrsaygovN4fJ95DZLw8eMQ8fBEAEDcC/Nx7ulnUPjOOwi4+moEXz/Ytg5GIy599hn+XPRPGIuLm5Vr6N0b7YYPQ1BKCvSJidAEBMBUUYG6ggLUnstD7blzqM07h7q8PBhLyyAY/KAx+EMTGAhNcDA0wUHQBgdD8A8AjHUQa2sbHnVGQDRBNJkAo6nR9CKh4VFXBeGXUOCXHKDw70BgeP17zp1rSezMmxbszqUWmu7kUFl2v0uT/ex+np1tb9SkYs0fp1FZdwhrt9+HW7qOsp2t5UQdmn2mI9/Z3n6ufp7dYx09pw58pqN/rq7UQSNAYzBA8DOY/074+0Oo/1njb4CmXTtoQ0Ig+Pu3cC0RkaMYRiQSm858+fMP83OjZeCXH1qOWlMtrom+BgNjB9otx1haiiErFiC4rgoXOvdEr8cft74XescdqNifhUtr1uDcE08g6Yv10MfFAQAqfv4Z+f/7Cqp/N4cgXXwcggYOBCCg6vBhVB87hursbFRnZ6NwyVL5vrjT6sfPZH+uYB28x53WV0dRiEWK1YMuT9DroQkNhTYkBLrwcOhiY6GPjYEuJhb6uFj4de4Mv06dIPg1HydGRGYMIxLZNIxUXgQq6rtFIroBMC9y9vkR8y/gqf2m2i+jrg5n5zyO0PNnUOgfgqwJT2Jwk3+4Yv42F1UHD6Lq8GGcHPsAwu66E1WHDqNsq3l1U01oKKJmPIb2Y8ZA0DX8sdZduICyrdtQ9sMPqDp0CLV5eeZZLRoNdFFR0MfFQR8fB318PHRxcdCGhEKsqYapqgqmigqYysphKiuDqbwcpspKCDodBL0egp9f/WsdIGgAjQaCRjC/hgiIojmoiTC3lpw/DBz/LxAQAfGq0W27mIq9qdIObWu+T7Pw2WJZLtah0YGnSnKxN/9nAMCAmAHoFNLJ8c+zt58LdTDv4uLnOVoHh8+zC59nd5uDn2c0Qaypgam6CmJ1DcTq6obXlZUwlpUBRiPE2loYCwthLCxEzYkTzcsBAK0WfgkJ8OvaFf49e8D/yisRcNVV0EVE2N+fSGUYRmQiCAJQdNz8Q7t4wBAMAPjw0IeoNlbjqqirMChuULPjjGVlOPfkUyjfvh11ej+8MOhh3BTW/B8ojcGADov/idxJk1B7KrehlUMQEDZmNKJmzoSufftmx+nCwxF2150Iu+tOADB3rZhM5jDRlk3LlReB17oDphxg8l0tLpVPDeIA7Nm/GO8ffB8Zmt+wbORs9I/ur3S1qJ4oiuaQXlwMY0kJjMUlqCsqRF3+edSezzc/nzuHmpwcmMrLUXPyJGpOnkTZf/9rLUPfoQMCBw5EUGoKAgcNgj46WsFvRKQchhGJbP5DVXjU/BxpbhW5UHXBenfeKVdNsfnlL4oiSjdtRsHrr6P2zBkIBgN2PfA4jl0Ix4gWWg38OnZE0rr1KF73OaoOZ0PfIR4ht90GQ1KSw/UV9Po2HaphFdAe6DYCOPItcPBz4Ma/K1ELrzO9/3TkFOfgu9zvMOuHWVg1ahU6BHdQuloE839AtMHB0AYHQ9+h5T8TURRRV1CAmuPHUX38BKoOH0blr7+i5sQJ1J49i+KzZ1H8n/8AAAw9eqDdiBvRLi0Nhp49ORaFVINhRCKbbpqi+jBS30Xz8eGPUVlXiT4RfTCkwxDrMZUHD+L8gldRuX8/AEAXF4cOb76B82cNwI85rd61VxschPDx493wTdrAVfeaw8iBlcCwuYCWl9/laAQN/vf6/8XZjWeRfSEb0/87HR+nf4xgv2Clq0YOEgQB+pgY6GNiEJSaat1uLC1F5S+/ouKn3Sjfucs8xuvIEVQfOYLCJUuh79QJ7W4agdBRo+Dfp4+C34DI/Ti1VyaCAKDomPmHiO4oqSnB6t/Ni5xZWkVEoxHnX12Ik/eORuX+/RACAhA5bRq6fvM1Avv3V6bFoi31ug0IjARK88yhhBwSqA/E4hsWIyogCscuHcNT255CnalO6WqRRNp27RB8/WBEP/44ktZ9ju47dyD+HwsRPOJGCAYDanNzcWFZBnLuvgcn7rwLFz76GHUXLypdbSK3YBiRyDLwTQCAQksY6YYvj32J8tpydAvrhmEJwwAA+S+9hAsrVgAwz5DpuvFbRM14DJog80wTa4usj9yaphmdH3DNQ+bXezOUrYuXiQ2KxeIbFsOgNWD72e2Yv3O+/UGX5LV07dsj9PbbkfD22+ixcwc6LFqEdjffDEGvR/Xvv+P8K6/g6F+G4szMWSj/aQ///MmnMIzIRBBNwAXzAFYxois+++MzAMD9ve6HRtCg+OuvcenTNYAgIP711xG/8FXoY2LsluXT/8QkTwAgAMe/bxjwSw7pG9kXrw99HVpBiy+Pf4m39r2ldJXITTRBQQi5eSQ6LnoL3bdvQ8zfn4X/FVcAtbUo3bQJuePHI+f2O3BxzWcwVXBhPPJ+DCMyaW/8E6irAjR6HKorw8mSkwjQBWBUl1EwlpXj/D/+AQCIfPRRhN46ym4Zqhis1r6zeSArAOx6R9GqeKNhCcMwP3U+APP6NR8c/EDZCpHbacPCEP7AA0ha9zmS/vMFwu6/D0JAAKqPHkX+88/j6LDhOL/wH6g5c0bpqhK5jGFEIktLaXT1afOL8C7IPPM9APPS70H6IFxcuRLGPwuhT+yEyP95pMWyGt8oz6ddP8v8vP8j4NJpRavije7sdifmJM8BAPxz/z/x/q/vK1wjaiv+vXoh7vnn0X3rFkQ/8zT0nTrBVFKCC8uX43jaSJx57DGU72EXDnkfhhGJLDNfomvqf6lGdsd/c83rCIxIHAGxrg4XPzUv8R459f85tAqjz/870vl6oPMQwFQL/Pim0rXxShP7TsT0q6cDABZnLcbSX5RcXZfamjYkBBETJqDrxm/R8d2lCBo8GDCZUJr5HXLHjUfOXXfj0vovvPdu36Q6DCMyiarJBQCcDo3BqZJT0Ak6DOkwBGVbtqAuLw/a9u0Rckt664WooJfGatgz5uf9HwMXWli1klo1pd8UzLxmJgBgyYElWLhnIYwmL7yzM7lM0GjQbtgwdFr2Abr83zcIu28MBH9/VP/+O/L+9jccG34D/lz8L9T9+afSVSVqFcOIRNZumvqWkV31S2f0i+5n7qJZZb5bb9hf74HGYGi1LMvN0Hy9YQSAuXWky3Bz68jXs1TQHOQek6+cjCcGPAEA+CT7E8z8YSbv9KtShq5dETd/Prpv+QHRTzwOXWwsjBcuoHDJEhy94Uace/ppVP52SOlqEtnFMCKR5VdoVI158NiuGvO9aVLiUlB9IgflO3cBGg3CxtzneJlq+b086g1A5w/kbAV+Wa10bbzW+CvG4/Whr8OgNWDrma146NuHkFOco3S1SCHasDBETJ6Mbpmb0eGtNxHQvz9QW4viL7/Cyb/+FSfHPoCSjZsg1tYqXVUiK4YRGRhQg7Da8zAC+KnYvAprSnwKLq0xLwUfPHQo/DpefglvNUymsRHRFRj6tPn1xrmc6ivByM4jsWzkMoT7h+PIxSMY880YrDuyjgMZVUzQ6xGSno7Oq1eh89rPEHLbbYBOh8r9+3F21izzLJxXF6LqyBGlq0rEMCKVKAJJQj40EHGoXThKa8vQzq8degd1td5vov19Yxwqq2HNMxX9Akl9DOiQDFRdAlbfB1QVK10jr9Uvqh8+u/UzXBd7HSrrKjF/13xM/346cktyla4aKSzgyivR4bV/oNt//4uI/zcV2ogIGIuKcGHFCuTcfgdy/novLnz8CWrPn1e6qqRSDCMSiRCRJOQBAHaFmRcxuy72OlRkfgdjcTF08XEIuv5658pUURaBVg/ct8p8p+PCI8Dq+4HKS0rXymvFBMXg32n/xpzkOdBpdNh2Zhvu/PJOvLXvLZTUlChdPVKYPiYa0TNnovuWH9BxyTsIHnEjoNOh6rffcP5//xfHhg7DyfvuR9GyDNSc5rR7ajsuhZElS5YgKSkJ/v7+SE5Oxvbt21vdf+vWrUhOToa/vz+6dOmCd99916XKeiJRBLrUh5Ef/cynMyU+BRfXmFdgbX/vvRC0WofKUl03jUW7WGDsp4BfO+DUDiDjZuDiSaVr5bU0ggYT+07EutvWYXD8YNSaapHxWwZuWnsTFu5ZiLNlZ5WuIilM0OvR7oYbkPD22+i+bSti5j5jHlsCoPLAARS89hqO35SG4yNvRt78+SjZuIn3xSG3cjqMrFmzBrNmzcK8efOQlZWFIUOGID09Hbm59puCc3JycMstt2DIkCHIysrC3/72N8yYMQPr1q2TXHlP0VeTg3M6LQ6YSiFAQGp+iPmOvDodQu++x+FyBFXN7W0irh8wcQMQHAv8mQ0sSQF2/guo4zoJruoS1gVLRyzF2ze8je7tu6OirgKfZH+C9HXpmLhxIj774zMUVhYqXU1SmC48HOHjx6Pz6lXotnUrYv7+LAKvuw7QaFBz6hQufbrGPMYkdTBO3H4Hzv1tHi6sXInKX36BqapK6eqTjxBEJ0e4XXfddbjmmmuwdGnDIku9e/fGnXfeiQULFjTb/+mnn8ZXX32F7Oxs67apU6fil19+wa5duxz6zJKSEoSGhqK4uBghISHOVNft3vn+KO7bOhxfhAn4Z3gYUsOT8fSyYlT//jvaj70fsc8953BZr2/6A2//cAzjUxLxwh193VhrD3bpNLD+f4Dcneafg6KB5PFAr1vNgUW1zUfSiKKIXed2Yfmh5didt9vmvS6hXTAwdiCuiLgCPcJ7oGtoV/jr/BWqKXkKY0kJKn7+GeW7dqN89y7UHLMzwFyrhb5jB/glJsKvc2fzc0ICdFFR0EVFQRseDkHD0QBq5ujvb50zhdbU1GDfvn145plnbLanpaVh586ddo/ZtWsX0tLSbLaNHDkSy5YtQ21tLfR6vTNVkNXx33ah9MJ5WCfoirAO2BAsGU0U698WG+bxiqL156DDv+NElREba9rjilMmPPJ/f6L69xPQhoYi8tFHnaqP5fdsaXUdzlxU61oREcCdaxF4aDVCfnoDurI8YNtrwLbXYDKEojaiJ2ojesEYHA9TQHuY/MNh8g+DqNFD1PoBGh1EjR7Q6iEKTf8RbBJkmgUboZX3vF8SIvBizydQkPgnthbsxNY/d+JI6XGcKD6BE8W2C8+114ch2j8SUYYItPcLQ5AuEMG6IATpghCoDYCfxg96jQ46QQedRlv/rINe0EGoP+8CLK195nMpCOafG7cACkKTn637+t7590odDcC9Q82Pi8XQHMmBcOwkhOOnoDl2CkJxKWpP5aL2VC7KtzXvrhc1GiAsBGL7ECAoEGKAPxAUADEgAAgKAAICIBr8AL0W0OnMD70O0Osh6rTm11qt+eLR1F9LGqH+YhJa3C5a32/zM9YKD6pMC3+/Ona7CpFxSW1cGTOnwkhhYSGMRiNimtxtNiYmBvn5+XaPyc/Pt7t/XV0dCgsLERcX1+yY6upqVDdaxrikxD0D736dNwW9/pA21/46AEAYXkJ9QMEJCH5+6Pj2v6CLiHCqLMvlsX7/Wazfr/Z+/XjosBDpmj24Vbsb12sOIqi6GIZze2A4t0fpynm1OAD9AMwAcEmjwT5/A/b6G3DEzw9H/PS4pNXiYu0lXKy9hD9KjylcW/I4XesfN4loX6ZF3AURcReAuIvm58gSEe3LgJByQGMyARcuQbhwSeFKkyN+GNMf976wSpHPdiqMWDT9X4soiq3+T8be/va2WyxYsAAvvPCCK1VzikmrQY2uvsFDaNTwUV8tsfHrRtus+zY6RitoERIaheBrBiBy2v+DoUsXp+szpEcUVu05jdIqLkZk5ofNuB6bxeuhN9ahC86ih3Aa3YTTiEQx2qMUYUIZQlAOHeqgRx30MEIHI/Sog1D/p9O8/aNxz6TYaLv9fQQfnmrtbwQGl9dhcHkdgHKIMAeU8zoNzuu0yNdpcVGrQZlGgzKNgFKNBuUaATWCgDoBqIOAWkFArQDrszmW264mbInqYqOTLDZ734P+50gOqQgCjgcJOJ4ANP2bpjGJCCkHwsqB0HIgsFpEQDXMjxrzc2A14FcnQmsE9EZA1+ShNwIaEyCIjR4wP2vqLx5N/XaIDa8t+7iL269UheouaF2KBLJw6pMjIyOh1WqbtYIUFBQ0a/2wiI2Ntbu/TqdDRAstB3PnzsWcOXOsP5eUlCAhIcGZqjrknvUHZC9TioGdw7H32RFKV4NULgDm1hMiUpd+Cn62UyOL/Pz8kJycjMzMTJvtmZmZSE1NtXtMSkpKs/03b96MAQMGtDhexGAwICQkxOZBREREvsnpYc5z5szBBx98gIyMDGRnZ2P27NnIzc3F1KlTAZhbNcaNG2fdf+rUqTh16hTmzJmD7OxsZGRkYNmyZXjiiSfk+xZERETktZzuIBozZgyKiorw4osvIi8vD3379sWGDRuQmJgIAMjLy7NZcyQpKQkbNmzA7Nmz8c477yA+Ph6LFy/GPfc4vv4GERER+S6n1xlRgievM0JERET2Ofr7m6vREBERkaIYRoiIiEhRDCNERESkKIYRIiIiUhTDCBERESmKYYSIiIgUxTBCREREimIYISIiIkUxjBAREZGiGEaIiIhIUU7fm0YJlhXrS0pKFK4JEREROcrye/tyd57xijBSWloKAEhISFC4JkREROSs0tJShIaGtvi+V9woz2Qy4dy5c2jXrh0EQZCt3JKSEiQkJOD06dO8Ad9l8Fw5h+fLcTxXjuO5chzPlXPcdb5EUURpaSni4+Oh0bQ8MsQrWkY0Gg06duzotvJDQkJ4sTqI58o5PF+O47lyHM+V43iunOOO89Vai4gFB7ASERGRohhGiIiISFGqDiMGgwHPP/88DAaD0lXxeDxXzuH5chzPleN4rhzHc+Ucpc+XVwxgJSIiIt+l6pYRIiIiUh7DCBERESmKYYSIiIgU5VNhZMmSJUhKSoK/vz+Sk5Oxffv2VvffunUrkpOT4e/vjy5duuDdd99tts+6devQp08fGAwG9OnTB1988YW7qt/m5D5fK1asgCAIzR5VVVXu/BptwplzlZeXh7Fjx6Jnz57QaDSYNWuW3f189dqS+1z58nUFOHe+1q9fj5tuuglRUVEICQlBSkoKNm3a1Gw/XluOnStfvracOVc//vgjBg8ejIiICAQEBKBXr1546623mu3n1utK9BGffvqpqNfrxffff188fPiwOHPmTDEoKEg8deqU3f1PnDghBgYGijNnzhQPHz4svv/++6Jerxc///xz6z47d+4UtVqt+Morr4jZ2dniK6+8Iup0OnH37t1t9bXcxh3na/ny5WJISIiYl5dn8/B2zp6rnJwcccaMGeKHH34oXn311eLMmTOb7eOr15Y7zpWvXlei6Pz5mjlzprhw4UJxz5494pEjR8S5c+eKer1e3L9/v3UfXltmjpwrX722nD1X+/fvF1etWiX+9ttvYk5Ojvjxxx+LgYGB4nvvvWfdx93Xlc+EkWuvvVacOnWqzbZevXqJzzzzjN39n3rqKbFXr14226ZMmSIOGjTI+vPo0aPFm2++2WafkSNHivfdd59MtVaOO87X8uXLxdDQUNnrqjRnz1VjQ4cOtfsL1levLXecK1+9rkRR2vmy6NOnj/jCCy9Yf+a11bKm58pXry05ztVdd90lPvjgg9af3X1d+UQ3TU1NDfbt24e0tDSb7Wlpadi5c6fdY3bt2tVs/5EjR2Lv3r2ora1tdZ+WyvQW7jpfAFBWVobExER07NgRt956K7KysuT/Am3IlXPlCF+8ttx1rgDfu64Aec6XyWRCaWkpwsPDrdt4bdln71wBvndtyXGusrKysHPnTgwdOtS6zd3XlU+EkcLCQhiNRsTExNhsj4mJQX5+vt1j8vPz7e5fV1eHwsLCVvdpqUxv4a7z1atXL6xYsQJfffUVVq9eDX9/fwwePBhHjx51zxdpA66cK0f44rXlrnPli9cVIM/5euONN1BeXo7Ro0dbt/Hass/eufLFa0vKuerYsSMMBgMGDBiARx99FJMnT7a+5+7ryitulOeopnf0FUWx1bv82tu/6XZny/Qmcp+vQYMGYdCgQdb3Bw8ejGuuuQb/+te/sHjxYrmqrQh3XAe+em3J/b18+boCXD9fq1evxvz58/Hll18iOjpaljI9ndznypevLVfO1fbt21FWVobdu3fjmWeeQbdu3XD//fdLKtNRPhFGIiMjodVqmyW0goKCZknOIjY21u7+Op0OERERre7TUpnewl3nqymNRoOBAwd69f8yXDlXjvDFa8td56opX7iuAGnna82aNZg0aRLWrl2LESNG2LzHa8tWa+eqKV+4tqScq6SkJADAlVdeifPnz2P+/PnWMOLu68onumn8/PyQnJyMzMxMm+2ZmZlITU21e0xKSkqz/Tdv3owBAwZAr9e3uk9LZXoLd52vpkRRxIEDBxAXFydPxRXgyrlyhC9eW+46V035wnUFuH6+Vq9ejQkTJmDVqlUYNWpUs/d5bTW43LlqyheuLbn+HoqiiOrqauvPbr+uZBkG6wEsU5mWLVsmHj58WJw1a5YYFBQknjx5UhRFUXzmmWfEhx56yLq/Zarq7NmzxcOHD4vLli1rNlV1x44dolarFV999VUxOztbfPXVV31iipwouud8zZ8/X9y4caN4/PhxMSsrS5w4caKo0+nEn376qc2/n5ycPVeiKIpZWVliVlaWmJycLI4dO1bMysoSDx06ZH3fV68td5wrX72uRNH587Vq1SpRp9OJ77zzjs1U1EuXLln34bVl5si58tVry9lz9fbbb4tfffWVeOTIEfHIkSNiRkaGGBISIs6bN8+6j7uvK58JI6Ioiu+8846YmJgo+vn5iddcc424detW63vjx48Xhw4darP/li1bxP79+4t+fn5i586dxaVLlzYrc+3atWLPnj1FvV4v9urVS1y3bp27v0abkft8zZo1S+zUqZPo5+cnRkVFiWlpaeLOnTvb4qu4nbPnCkCzR2Jios0+vnptyX2ufPm6EkXnztfQoUPtnq/x48fblMlry7Fz5cvXljPnavHixeIVV1whBgYGiiEhIWL//v3FJUuWiEaj0aZMd15XvGsvERERKconxowQERGR92IYISIiIkUxjBAREZGiGEaIiIhIUQwjREREpCiGESIiIlIUwwgREREpimGEiIiIFMUwQkRERIpiGCEiIiJFMYwQkWKSkpKwceNGpatBRApjGCEiRfz6668oKirC8OHDla4KESmMYYSIJDl58iQEQWj2GDZsWKvHffnllxg5ciQMBoPd9wVBwHvvvYdbb70VgYGB6N27N3bt2oVjx45h2LBhCAoKQkpKCo4fP+6Gb0VEbYlhhIgkSUhIQF5envWRlZWFiIgI/OUvf2n1uK+++gp33HFHq/u89NJLGDduHA4cOIBevXph7NixmDJlCubOnYu9e/cCAKZPny7bdyEiZQiiKIpKV4KIfENVVRWGDRuGqKgofPnll9Bo7P9/5+zZs0hKSkJ+fj7Cw8Pt7iMIAp599lm89NJLAIDdu3cjJSUFy5Ytw8MPPwwA+PTTTzFx4kRUVla65wsRUZvQKV0BIvIdkyZNQmlpKTIzM1sMIoC5VWTw4MEtBhGLq666yvo6JiYGAHDllVfabKuqqkJJSQlCQkIk1p6IlMIwQkSyePnll7Fx40bs2bMH7dq1a3VfR7poAECv11tfC4LQ4jaTyeRKlYnIQ3DMCBFJtm7dOrz44ov47LPP0LVr11b3LSsrww8//IDbb7+9jWpHRJ6OLSNEJMlvv/2GcePG4emnn8YVV1yB/Px8AICfn5/dbpiNGzeie/fu6NKlS1tXlYg8FFtGiEiSvXv3oqKiAi+//DLi4uKsj7vvvtvu/l9++aVDXTREpB6cTUNEbcZoNCI6Ohrffvstrr32WqWrQ0Qegi0jRNRmioqKMHv2bAwcOFDpqhCRB2HLCBERESmKLSNERESkKIYRIiIiUhTDCBERESmKYYSIiIgUxTBCREREimIYISIiIkUxjBAREZGiGEaIiIhIUQwjREREpKj/D46VLx9XdE8+AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for ix, (tprof, phiprof) in enumerate(phiprofs):\n",
    "    if ix == 2:\n",
    "        lbl1 = 'FVM'\n",
    "    else:\n",
    "        lbl1 = None\n",
    "    plt.plot(phi.domain.cellcenters.z, phiprof, \n",
    "             label=lbl1)\n",
    "plt.xlabel('z / m')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "976a4cd4-5b88-4f0b-a184-dbe4dcde8aeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "exect1 = time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "1488bee9-d24f-411f-8886-ba0a105f309c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time  6.2673749923706055 s\n"
     ]
    }
   ],
   "source": [
    "print('Elapsed time ', exect1 - exect0, 's')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c15e591f-994b-44ba-b069-06581db6d907",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}