{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9667595f-10e1-4605-b8e0-c42f65094697",
   "metadata": {},
   "source": [
    "# Cartesian 1D: Solve 1D diffusion with the explicit solver\n",
    "\n",
    "A tutorial example adapted from the FiPy 1D diffusion example\n",
    "\n",
    "see: https://pages.nist.gov/fipy/en/stable/generated/examples.diffusion.mesh1D.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "eb280c09-a91f-4433-8f6d-1254aefb57bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyfvtool as pf\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.special import erf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "97430bc0-7240-472f-a8db-7494b2c7fb04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the domain\n",
    "L = 5.0  # domain length\n",
    "Nx = 100 # number of cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0ecbec02-e1dd-4012-8bca-0bfda3ac0bb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "meshstruct = pf.Grid1D(Nx, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "81face1d-06c0-4170-80f2-765ee4659f7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "BC = pf.BoundaryConditions(meshstruct) # all Neumann boundary condition structure\n",
    "BC.left.a[:] = 0 \n",
    "BC.left.b[:]=1 \n",
    "BC.left.c[:]=1 # left boundary\n",
    "BC.right.a[:] = 0 \n",
    "BC.right.b[:]=1 \n",
    "BC.right.c[:]=0 # right boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "664b1c8f-c769-4312-a034-80eb2c6b8d0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = meshstruct.cellcenters.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d201db49-956f-4778-9294-d321976e29e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "## define the transfer coeffs\n",
    "D_val = 1.0\n",
    "alfa = pf.CellVariable(meshstruct, 1)\n",
    "Dave = pf.FaceVariable(meshstruct, D_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "38b12682-bf11-4843-87c9-8883bda4b1de",
   "metadata": {},
   "outputs": [],
   "source": [
    "## define initial values\n",
    "c_old = pf.CellVariable(meshstruct, 0, BC) # initial values\n",
    "c = pf.CellVariable(meshstruct, 0, BC) # working values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1b4c90ca-68ea-46c3-abbf-9488701b0f0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "## loop\n",
    "dt = 0.001 # time step\n",
    "final_t = 0.5\n",
    "for t in np.arange(dt, final_t, dt):\n",
    "    # step 1: calculate divergence term\n",
    "    RHS = pf.divergenceTerm(Dave*pf.gradientTerm(c_old))\n",
    "    # step 2: calculate the new value for internal cells\n",
    "    c = pf.solveExplicitPDE(c_old, dt, RHS)\n",
    "    c_old.update_value(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3aaaa61d-398e-4c71-87b5-3fab3c243f1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytic solution\n",
    "c_analytic = 1-erf(x/(2*np.sqrt(D_val*t)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "06e4c664-37b2-41a7-88cb-ddc12fa106cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJ+0lEQVR4nO3deVxU9eLG8c8w7Ai4gygi7guuuIHi0mKaWd5rqZW7Vtrq1mKWqVm2mnVTW1yyMrXFtPqZxa0Ml9xI1ITMHVEUsQRBZZk5vz9ErgQYKHBYnvfrdV6XOfM9M8+M3Obhe86cYzEMw0BERETEJA5mBxAREZGKTWVERERETKUyIiIiIqZSGRERERFTqYyIiIiIqVRGRERExFQqIyIiImIqlRERERExlaPZAQrCbrdz4sQJPD09sVgsZscRERGRAjAMg3PnzuHn54eDQ/7zH2WijJw4cQJ/f3+zY4iIiMg1OHbsGHXq1Mn3/jJRRjw9PYFLL8bLy8vkNCIiIlIQycnJ+Pv7Z3+O56dMlJHLu2a8vLxURkRERMqYfzrEQgewioiIiKlURkRERMRUhS4jERER9OvXDz8/PywWC6tXr/7HbX7++WeCg4NxdXWlfv36vPPOO9eSVURERMqhQh8zkpqaSuvWrRk5ciQDBgz4x/GHDx/m1ltv5b777uPjjz9m06ZNPPjgg9SoUaNA24uIVASGYZCZmYnNZjM7ikiBWa1WHB0dr/u0G4UuI3369KFPnz4FHv/OO+9Qt25d5s6dC0CzZs3YsWMHr732msqIiAiQnp5OfHw858+fNzuKSKG5u7tTq1YtnJ2dr/kxiv3bNL/88gu9evXKse6WW25h0aJFZGRk4OTkVNwRRERKLbvdzuHDh7Farfj5+eHs7KyTO0qZYBgG6enpnD59msOHD9OoUaOrntjsaoq9jJw8eRIfH58c63x8fMjMzCQxMZFatWrl2iYtLY20tLTs28nJycUdU0TEFOnp6djtdvz9/XF3dzc7jkihuLm54eTkxNGjR0lPT8fV1fWaHqdEvk3z95ZvGEae6y+bPXs23t7e2YvOvioi5d21/kUpYrai+N0t9t9+X19fTp48mWNdQkICjo6OVKtWLc9tpkyZQlJSUvZy7Nix4o4pIiIiJin2MhISEkJ4eHiOdd9//z3t27fP93gRFxeX7LOt6qyrIiJilvXr12OxWDh79uxVx/344480bdoUu91eMsGAI0eOYLFYiIqKAgqe9bIePXowfvz4q4755ptvaNu2bbG/rkKXkZSUFKKiorJf/OHDh4mKiiI2Nha4NKsxbNiw7PFjx47l6NGjTJw4kZiYGBYvXsyiRYuYPHly0bwCERExxYgRI7BYLFgsFpycnKhfvz6TJ08mNTW1wI9Rr1697Me4vNSpU4f09HSqV6/OrFmz8txu9uzZVK9endDQ0FzbX7nUq1eviF7t1T3xxBNMnTrV1N1toaGhxMfH4+3tXaDxq1at4vnnn8++Xa9evexvvl522223YbFY+OSTT4oyai6Fftd27NhB27Ztadu2LQATJ06kbdu2TJs2DYD4+PjsYgIQGBjI2rVrWb9+PW3atOH555/nrbfe0td6RUTKgd69exMfH8+hQ4eYNWsW8+fPL/QfmzNnziQ+Pj572blzJ87OzgwZMoQPPvgg+zjDKy1ZsoShQ4fyzTffZG+3bds2AP773/9mr9u+fXuRvM6r2bx5M/v37+euu+4q9ue6GmdnZ3x9fQv8bayqVav+4wXsAEaOHMl//vOf6413dUYZkJSUZABGUlJSkT7u6tWrjYEDBxrR0dFF+rgiIgV14cIFIzo62rhw4YLZUQpt+PDhxh133JFj3ZgxYwxfX1/DbrcbDRo0MF599dUc9+/Zs8ewWCzGgQMHDMMwjICAAOONN97I8/F3795tAMb69etzrI+IiDAAY8+ePTnWHz582ACMnTt3Zq9bv3690aFDB8PZ2dnw9fU1nnzySSMjIyP7/osXLxqPPPKIUaNGDcPFxcXo0qWLsW3btuz7f/rpJwMw/vrrr3zfh0ceecS48847c63/6quvjHbt2hkuLi5GYGCgMX369OznnjFjhlGrVi0jMTExe3y/fv2MsLAww2azGYZhGIAxf/58o3fv3oarq6tRr14949NPP8339eaVdePGjUa3bt0MNzc3o3LlykavXr2MP//80zAMw+jevbvx2GOPZf8M5FguO3LkiAEYBw8ezPP1X+13uKCf3xX68O3v5syh86efcnjSJLOjiIhkMwyD1NTUEl+MPGYgCsvNzY2MjAwsFgujRo1iyZIlOe5fvHgxYWFhNGjQ4B8fq2XLlnTo0CHPx+jYsSNBQUFX3f748ePceuutdOjQgV27drFgwQIWLVqUY9fPE088wRdffMHSpUv59ddfadiwIbfccgt//vlngV9zREQE7du3z7Huu+++Y8iQITz66KNER0fz7rvv8sEHH/DCCy8AMHXqVOrVq8eYMWOASycIjYiI4KOPPsqxq+fZZ59lwIAB7Nq1iyFDhnD33XcTExNToFxRUVHceOONtGjRgl9++YWNGzfSr1+/PM/yu2rVKurUqZNjluqygIAAatasyYYNGwr8nhTaVatKKVFcMyMb7rvPMMDY7+JSpI8rIlJQef1VmZKSkuuv1JJYUlJSCpX97zMjW7duNapVq2YMHDjQMAzDOHHihGG1Wo2tW7cahmEY6enpRo0aNYwPPvgge5uAgADD2dnZ8PDwyF7efPPN7PsXLFhgeHh4GOfOnTMMwzDOnTtneHh4GO+++26uPH+fKXj66aeNJk2aGHa7PXvMvHnzjEqVKhk2m81ISUkxnJycjGXLlmXfn56ebvj5+RmvvPKKYRgFmxnx9vY2PvzwwxzrwsLCjBdffDHHuo8++sioVatW9u2DBw8anp6expNPPmm4u7sbH3/8cY7xgDF27Ngc6zp16mSMGzcuz9f796x333230aVLl3xzXzkzYhhXn6Vq27atMX369Dzv08zIdWo5dSrpQMO0NA589ZXZcUREypxvvvmGSpUq4erqSkhICN26dcs+vqBWrVr07duXxYsXZ4+9ePFirmMrHn/88ewvRkRFReX4EsTdd9+N3W5n5cqVAKxcuRLDMBg8ePA/ZouJiSEkJCTHMRRdunQhJSWFuLg4Dh48SEZGBl26dMm+38nJiY4dOxZ49gHgwoULuU72FRkZycyZM6lUqVL2ct999+U47X/9+vV57bXXePnll+nXrx/33ntvrscOCQnJdbuwMyNFwc3NrVgvV1DsZ2AtzbwDAthWsyYdExKIe+MNGt5+u9mRRERwd3cnJSXFlOctrJ49e7JgwQKcnJzw8/PLdcqGMWPGMHToUN544w2WLFnCoEGDcj1P9erVadiwYZ6P7+3tzZ133smSJUsYPXo0S5Ys4c477yzQKR8Mw7jqSTeNfE7Amdd2V1O9enX++uuvHOvsdjszZszg3//+d67xVxaXiIgIrFYrR44cITMzE0fHf/5YLmg2Nze3Ao0riD///JMaNWoU2eP9XYWeGQFIv+MOAGr/8ovJSURELrFYLHh4eJT4ci3XxPHw8KBhw4YEBATkee6oW2+9FQ8PDxYsWMC3337LqFGjCv0co0ePZtOmTXzzzTds2rSJ0aNHF2i75s2bs3nz5hzHwmzevBlPT09q165Nw4YNcXZ2ZuPGjdn3Z2RksGPHDpo1a1bgfG3btiU6OjrHunbt2rFv3z4aNmyYa7l8TMjKlStZtWoV69ev59ixYzm+ZnvZli1bct1u2rRpgXK1atWKH374ocCvw9nZOc/jSS5evMjBgwezv0VbHCp8Gbm8q6aRdtWIiBQ5q9XKiBEjmDJlCg0bNsy126EgunfvTsOGDRk2bBgNGzakW7duBdruwQcf5NixYzzyyCP8/vvvrFmzhueee46JEyfi4OCAh4cH48aN4/HHH2fdunVER0dz3333cf78+QIXHrh08dcrCw3AtGnT+PDDD5k+fTp79+4lJiaGlStX8swzzwAQFxfHuHHjePnll+natSsffPABs2fPzlU+PvvsMxYvXswff/zBc889x7Zt23j44YcLlGvKlCls376dBx98kN27d/P777+zYMECEhMT8xxfr149IiIiOH78eI4xW7ZswcXF5Zr+7QqqwpcR74AAomrWBCDujTdMTiMiUv6MHj2a9PT0a5oVuWzUqFH89ddfhXqM2rVrs3btWrZt20br1q0ZO3Yso0ePzi4EAC+99BIDBgxg6NChtGvXjgMHDvDdd99RpUqVAj/PkCFDiI6OZt++fdnrbrnlFr755hvCw8Pp0KEDnTt3Zs6cOQQEBGAYBiNGjKBjx47ZxeLmm2/m4YcfZsiQITl20c2YMYMVK1bQqlUrli5dyrJly2jevHmBcjVu3Jjvv/+eXbt20bFjR0JCQlizZk2+u4JmzpzJkSNHaNCgQY5dMsuXL+fee+8t1gs5WgyjCL7LVcySk5Px9vYmKSmpWE4Nv/H++wl8/31W+Pgw6W/X0RERKU4XL17k8OHDBAYGXvMVT0u7TZs20aNHD+Li4nJdxb28eOKJJ0hKSuLdd98tsse0WCx8+eWX9O/fv8ges7BOnz5N06ZN2bFjB4GBgXmOudrvcEE/vyv8zAhA0AsvUN/RkcmnTuXa7yciItcmLS2NAwcO8OyzzzJw4MByW0Tg0nlDAgIC8jzmoiw7fPgw8+fPz7eIFBWVEaByjRrcfMstwKX9cyIicv2WL19OkyZNSEpK4pVXXjE7TrHy9vbm6aefxmq1mh2lSHXs2JFBgwYV+/OojGQZOHAgDsD+pUvNjiIiUi6MGDECm81GZGQktWvXNjtOmWMYhqm7aEqSykiW2/v04TDw8eHDHPj6a7PjiIiIVBgqI1kq16jByaxv1RybM8fkNCIiIhWHysgVLp8ArY5OgCYiIlJiVEaukOMEaNpVIyIiUiJURq7gHRDALu2qERERKVEqI3+TnnXksnbViIiIlAyVkb8Jevrp7F01+3WtGhERU/Xo0YPx48df12McOXIEi8VCVFRUkWSSoqcy8jfeAQG827o1HYBlkZFmxxERkUIYMWJErnNz+Pv7Ex8fT1BQkDmh5B+pjOSh8qRJ7AA++/xzs6OIiMh1slqt+Pr65nuBODGfykgebr/9dpydnYmOjmbv3r1mxxERKZXWrVtH165dqVy5MtWqVeO2227j4MGDwP92jaxatYqePXvi7u5O69at+eWK4/HOnDnD3XffTZ06dXB3d6dly5YsX7483+ebOXMmLVu2zLU+ODiYadOmMX36dJYuXcqaNWuwWCxYLBbWr1+f526avXv30rdvX7y8vPD09CQsLCw7u5Q8lZE8eHt783DnziwCjk2YYHYcEamoUlPzXy5eLPjYCxf+eew1xUtl4sSJbN++nR9++AEHBwf+9a9/Ybfbs8dMnTqVyZMnExUVRePGjbn77rvJzMwELl3tNTg4mG+++YbffvuN+++/n6FDh7J169Y8n2/UqFFER0ezffv27HW7d+9m586djBgxgsmTJzNw4EB69+5NfHw88fHxhIaG5nqc48eP061bN1xdXfnxxx+JjIxk1KhR2bmk5GnOKh//btaMLhERHPz5ZwzDwGKxmB1JRCqaSpXyv+/WW+H//u9/t2vWhPPn8x7bvTusX/+/2/XqQWJizjGGUeh4AwYMyHF70aJF1KxZk+joaCplZZ88eTJ9+/YFYMaMGbRo0YIDBw7QtGlTateuzeTJk7O3f+SRR1i3bh2fffYZnTp1yvV8derU4ZZbbmHJkiV06NABgCVLltC9e3fq168PgJubG2lpafj6+uabe968eXh7e7NixQqcnJwAaNy4caFfvxQdzYzko+XUqaQBDdLTObBmjdlxRERKnYMHD3LPPfdQv359vLy8si8zHxsbmz2mVatW2T/XqlULgISEBABsNhsvvPACrVq1olq1alSqVInvv/8+x/Z/d99997F8+XIuXrxIRkYGy5YtY9SoUYXKHRUVRVhYWHYREfNpZiQfXv7+bPP1pePJk8TNmUOjCnLlRBEpRVJS8r/v75eqz/qAz5PD3/7uPHLkmiNdqV+/fvj7+/P+++/j5+eH3W4nKCiI9PT07DFXfuBfnmG+vBvn9ddf54033mDu3Lm0bNkSDw8Pxo8fn2P7vJ7TxcWFL7/8EhcXF9LS0nLN0PwTNze3Qo2X4qcychUZ/fvDO+9Qd+tWDLsdy9//Dy0iUpw8PMwfm48zZ84QExPDu+++S1hYGAAbN24s1GNs2LCBO+64gyFDhgCXSsr+/ftp1qxZvts4OjoyfPhwlixZgouLC4MHD8bd3T37fmdnZ2w221Wft1WrVixdupSMjAzNjpQS+nS9ipZTp3KRS7tq9q9ebXYcEZFSo0qVKlSrVo333nuPAwcO8OOPPzJx4sRCPUbDhg0JDw9n8+bNxMTE8MADD3Dy5Ml/3G7MmDH8+OOPfPvtt7l20dSrV4/du3ezb98+EhMTycjIyLX9ww8/THJyMoMHD2bHjh3s37+fjz76iH379hUqvxQdlZGr8KpTh11ZB0Ed17VqRESyOTg4sGLFCiIjIwkKCmLChAm8+uqrhXqMZ599lnbt2nHLLbfQo0cPfH19c52wLC+NGjUiNDSUJk2a5DrQ9b777qNJkya0b9+eGjVqsGnTplzbV6tWjR9//JGUlBS6d+9OcHAw77//vmZJTGQxjGs4hLqEJScn4+3tTVJSEl5eXiX63Jsfeohq8+ezomZNpp08qW/ViEiRunjxIocPHyYwMBBXV1ez45QJhmHQtGlTHnjggULPxkjRu9rvcEE/vzUz8g9avvACbVxcmJ6QwO7du82OIyJSoSUkJDBnzhyOHz/OyJEjzY4jRURl5B94Vq7MrVnfkf/0009NTiMiUrH5+Pjw0ksv8d5771GlShWz40gRURkpgIEDB+ICnFyyBOOKMwuKiEjJMgyD06dPc88995gdRYqQykgB9O3Vi6PAovh49ml2REREpEipjBRApSpVOFy7NgDxb75pchoREZHyRWWkgIy77gKg/o4d2lUjIkWuDHyxUSRPRfG7qzJSQK2mTOE8EJCZScyyZWbHEZFy4vK5Lc7nd5E7kVLu8u/u9ZynRaeDLyCPmjX5xd+fkGPHOPX22zQfOtTsSCJSDlitVipXrpx98Th3d3edz0jKBMMwOH/+PAkJCVSuXBnr36+XVAgqI4VgGTQIXnuNhpGRulaNiBSZy5e7T7jaxe5ESqnKlStn/w5fK52BtRAunDmDvXp1PIA9H31Ey6yLO4mIFAWbzZbntVRESisnJ6erzogU9PNbMyOF4FatGm907co7Gzdy286dvK4yIiJFyGq1XtdUt0hZpf0MhRQ4aRJ/cOlsrHZ9q0ZEROS6qYwUUu/evfH09CQuLo4tW7aYHUdERKTMUxkpJFdXVyaGhvIZkPrYY2bHERERKfNURq5B3w4duBNo8euv2HSwmYiIyHVRGbkGrZ94giTAz25nzzvvmB1HRESkTFMZuQbOnp7sadAAgLPvvWdyGhERkbJNZeQauQ4fDkDzvXvJTEszOY2IiEjZpTJyjVpPmsSfFgs1DYNdb71ldhwREZEyS2XkGjm5u7O3cWMAUhcvNjmNiIhI2aUych08Ro9mO7A6NlancBYREblGKiPXodWECdxWsyZvnD/PDz/8YHYcERGRMkll5Do4Ojpy5513ArBy5UqT04iIiJRNKiPXadCgQXgBzitXknbunNlxREREyhyVkevUtUsX9litvHvhAlGvvmp2HBERkTJHZeQ6OVitHAkKAiDjo49MTiMiIlL2qIwUgSoPPghAqyNHOH/mjMlpREREyhaVkSIQNHo0x61WvICol14yO46IiEiZojJSBCxWKwfatbt0Y8UKc8OIiIiUMddURubPn09gYCCurq4EBwezYcOGq45ftmwZrVu3xt3dnVq1ajFy5EjOlLPdGT6PPQZAm7g4zp04YXIaERGRsqPQZWTlypWMHz+eqVOnsnPnTsLCwujTpw+xsbF5jt+4cSPDhg1j9OjR7N27l88++4zt27czZsyY6w5fmjS5+26OODnhCmx/802z44iIiJQZhS4jc+bMYfTo0YwZM4ZmzZoxd+5c/P39WbBgQZ7jt2zZQr169Xj00UcJDAyka9euPPDAA+zYseO6w5cmFgcH/jtkCP7AG9HRZscREREpMwpVRtLT04mMjKRXr1451vfq1YvNmzfnuU1oaChxcXGsXbsWwzA4deoUn3/+OX379s33edLS0khOTs6xlAWhkydzAvjuu+/4888/zY4jIiJSJhSqjCQmJmKz2fDx8cmx3sfHh5MnT+a5TWhoKMuWLWPQoEE4Ozvj6+tL5cqV+c9//pPv88yePRtvb+/sxd/fvzAxTdO8eXNatWpFRkYGq7/4wuw4IiIiZcI1HcBqsVhy3DYMI9e6y6Kjo3n00UeZNm0akZGRrFu3jsOHDzN27Nh8H3/KlCkkJSVlL8eOHbuWmKZ4JCyM74AWU6aYHUVERKRMcCzM4OrVq2O1WnPNgiQkJOSaLbls9uzZdOnShccffxyAVq1a4eHhQVhYGLNmzaJWrVq5tnFxccHFxaUw0UqNm//1LwLmzSPzzBlOR0dTo3lzsyOJiIiUaoWaGXF2diY4OJjw8PAc68PDwwkNDc1zm/Pnz+PgkPNprFYrcGlGpbwJuPFGYtzdcQSin3/e7DgiIiKlXqF300ycOJGFCxeyePFiYmJimDBhArGxsdm7XaZMmcKwYcOyx/fr149Vq1axYMECDh06xKZNm3j00Ufp2LEjfn5+RfdKSpGEG24AwOvbb01OIiIiUvoVajcNwKBBgzhz5gwzZ84kPj6eoKAg1q5dS0BAAADx8fE5zjkyYsQIzp07x9tvv82kSZOoXLkyN9xwAy+//HLRvYpSpvEzz8A339A6KYkT27fj16GD2ZFERERKLYtRBvaVJCcn4+3tTVJSEl5eXmbHKZDd3t60Sk7mp9tvp+eaNWbHERERKXEF/fzWtWmKydk+fQCo+d//mpxERESkdFMZKSbNp03jR2Du+fP8sW+f2XFERERKLZWRYlK9eXNe7d2bhcByXclXREQkXyojxejuu+8GYPny5eXya8wiIiJFQWWkGPXv358AFxd67dtHzJdfmh1HRESkVFIZKUZeXl6sqF6dt4BT5firzCIiItdDZaS4Ze2qaRgZid1mMzmMiIhI6aMyUszaPP00yYC/zcae994zO46IiEipozJSzFyrVGFP/foA/DV/vslpRERESh+VkRLgOnIkAM337iXjwgWT04iIiJQuKiMloPWkSZyxWKhpGETNmWN2HBERkVJFZaQEOLq5Ed28OenAvtWrzY4jIiJSqqiMlBDXF17AFxj3++9c0K4aERGRbCojJaT97bfjXa8eKSkpfP3112bHERERKTVURkqIxWLJPj3814sXm5xGRESk9FAZKUFDb7+dzcD7333HX4cOmR1HRESkVFAZKUHNOnemuqsrrsDu6dPNjiMiIlIqqIyUsBM9egDgqeNGREREAJWREtfwmWcAaHP2LCe2bzc5jYiIiPlURkpY7S5d2O3lhQOwb8YMs+OIiIiYTmXEBEm33QaAz48/mpxERETEfCojJgiaPp0MoPmFC+zXsSMiIlLBqYyYoEqjRqxs2pRBwMcbNpgdR0RExFQqIyZxmj6dT4EPP/sMwzDMjiMiImIalRGT9OvXj0qVKnHkyBE2b95sdhwRERHTqIyYxN3dndG9e/M0kDBpktlxRERETKMyYqJhbdvyAtBt2zbSU1PNjiMiImIKlRETtZ40iQQHB6oZBr/Onm12HBEREVOojJjI6uLC723aAGBbutTcMCIiIiZRGTGZz+TJALSLi+NsbKzJaUREREqeyojJGg8axEEXF9yAXc89Z3YcERGREqcyYjKLgwNxWVfyrbR6talZREREzKAyUgo0eu45zgF7zp7l6KFDZscREREpUSojpYBfSAgDwsIYCXyycqXZcUREREqUykgpMXjECAA++ugjnR5eREQqFJWRUmLAgAG4ODvjFBPD3m+/NTuOiIhIiVEZKSW8vb35IiCAXcDpadPMjiMiIlJiVEZKkZp33QVAs507yUxLMzmNiIhIyVAZKUVaP/UUf1os+Nrt7HztNbPjiIiIlAiVkVLE2dOTvUFBAKS9/77JaUREREqGykgpU33SJADaHT1KUlycyWlERESKn8pIKdN06FAOOTvjDkQ9+6zZcURERIqdykgpY3FwIK5nz0s/f/21yWlERESKn8pIKdTohRe4DbjhzBkO6fTwIiJSzqmMlEK1goNJv/lmbMDHH39sdhwREZFipTJSSg0fPhyAj5cuxbDbTU4jIiJSfFRGSqn+/fszw9mZHw4dYre+5isiIuWYykgp5eHhwU116uAP/PXWW2bHERERKTYqI6WYx4MPAtA6JoaLSUkmpxERESkeKiOlWMtHHyXeaqWKYRA5c6bZcURERIqFykgp5uDkxP5OnQBw/Ogjk9OIiIgUD5WRUi7guecACD59mvidO01OIyIiUvRURkq5gF69+M3TE0fg92eeMTuOiIhIkVMZKQMSBw9mATB3714MwzA7joiISJFSGSkDgl9/ncnu7nx19Chbt241O46IiEiRUhkpAzw9PbnzzjsBWLJkiclpREREipbKSBkxYvhwQoDgJUu48NdfZscREREpMiojZUT3sDA+s1q5PyODyGnTzI4jIiJSZK6pjMyfP5/AwEBcXV0JDg5mw4YNVx2flpbG1KlTCQgIwMXFhQYNGrB48eJrClxROTg5cSA0FACX5ctNTiMiIlJ0Cl1GVq5cyfjx45k6dSo7d+4kLCyMPn36EBsbm+82AwcO5IcffmDRokXs27eP5cuX07Rp0+sKXhHVnzEDgHZnznB82zaT04iIiBQNi1HI74p26tSJdu3asWDBgux1zZo1o3///syePTvX+HXr1jF48GAOHTpE1apVrylkcnIy3t7eJCUl4eXldU2PUV7sqlyZ1klJ/HTzzfT8/nuz44iIiOSroJ/fhZoZSU9PJzIykl69euVY36tXLzZv3pznNl999RXt27fnlVdeoXbt2jRu3JjJkydz4cKFfJ8nLS2N5OTkHItckpr1rZp669dj2O0mpxEREbl+hSojiYmJ2Gw2fHx8cqz38fHh5MmTeW5z6NAhNm7cyG+//caXX37J3Llz+fzzz3nooYfyfZ7Zs2fj7e2dvfj7+xcmZrnWetYsUoDAjAyi5s0zO46IiMh1u6YDWC0WS47bhmHkWneZ3W7HYrGwbNkyOnbsyK233sqcOXP44IMP8p0dmTJlCklJSdnLsWPHriVmueTh68uuxo05Avy8apXZcURERK5bocpI9erVsVqtuWZBEhIScs2WXFarVi1q166Nt7d39rpmzZphGAZxcXF5buPi4oKXl1eORf7HZf586gNPb93K2bNnzY4jIiJyXQpVRpydnQkODiY8PDzH+vDwcEKzvnb6d126dOHEiROkpKRkr/vjjz9wcHCgTp061xBZgm+4geYtWnDhwgWW62u+IiJSxhV6N83EiRNZuHAhixcvJiYmhgkTJhAbG8vYsWOBS7tYhg0blj3+nnvuoVq1aowcOZLo6GgiIiJ4/PHHGTVqFG5ubkX3SioQi8XC6NGjcQL2vPGG2XFERESuS6HLyKBBg5g7dy4zZ86kTZs2REREsHbtWgICAgCIj4/Pcc6RSpUqER4eztmzZ2nfvj333nsv/fr146233iq6V1EBDb3rLg4B8/fvZ98XX5gdR0RE5JoV+jwjZtB5RvK2rXZtOp44wU9t2tBz506z44iIiORQLOcZkdLFev/9ALTatYuLSUkmpxEREbk2KiNlWJunnuKk1Uo1w2CHLp4nIiJllMpIGWZ1ceGPkBAAXD7+2OQ0IiIi10ZlpIyrP2sWAMF//knsxo0mpxERESk8lZEyrk737kRVqYIDsDOPCxWKiIiUdioj5cDpyZNpCTwYFUVmZqbZcURERApFZaQc6DZpEierV+fEiROsXbvW7DgiIiKFojJSDri4uDBixAgAli5YYG4YERGRQlIZKSfuGzGCD4Gl69ZxfNs2s+OIiIgUmMpIOdG4RQtae3tTCdj31FNmxxERESkwlZFyJPWeewBoFBFBZlqayWlEREQKRmWkHGn7wgv8ZbHgb7MR+dJLZscREREpEJWRcsS1ShX2tGkDgP2dd8wNIyIiUkAqI+VMnRkzAOhw8iQnIiNNTiMiIvLPVEbKmfr9+rHbywtH4PcnnzQ7joiIyD9SGSmH/hw9mieByTEx2Gw2s+OIiIhclcpIOdT5xRdZWLUqO3VGVhERKQNURsohV1dXRo4cCcACnZFVRERKOZWRcuqBBx7gLmDyt98Su3Gj2XFERETypTJSTjVq1IgpVatyA3DgiSfMjiMiIpIvlZFyLGP0aACCtmwh7dw5k9OIiIjkTWWkHGs3YwanHByoaRhsf/pps+OIiIjkSWWkHHN0c2NfWBgAHh9+aHIaERGRvKmMlHONX3uNTKBtcjJ/rF5tdhwREZFcVEbKOd/27fnVzw+AE88+a3IaERGR3FRGKgCnxx5jA/D+gQOc04GsIiJSyqiMVABtHn+c+5o04ZOLF/n444/NjiMiIpKDykgFYLFYGDduHADz5s3DMAyTE4mIiPyPykgFMWLECPzd3bl1715+1SniRUSkFFEZqSC8vb1ZVq8erwAXX3zR7DgiIiLZVEYqEL9ZswDodPw4cVu2mJxGRETkEpWRCqTBv/7FrsqVcQT2TZxodhwRERFAZaTCSbv/fgBabtnChbNnzQ0jIiKCykiFEzxjBvFWKzUNg22TJ5sdR0RERGWkorG6urL/ppsAqLZ8ub7mKyIiplMZqYCC5s4lGYg6f54tP/1kdhwREangVEYqoKpNmzJl2DCGAm++957ZcUREpIJTGamg7s/6Ns0XX3zB8ePHTU4jIiIVmcpIBdW6dWu6detG48xMfhg/3uw4IiJSgamMVGCzbrqJvcCtX3zB+TNnzI4jIiIVlMpIBRb6xBPEOTpS3TDYptkRERExicpIBWZ1ceFQ374A+H36KXabzeREIiJSEamMVHBt336bc0Dj9HR2vPSS2XFERKQCUhmp4Dzr1CGqTZtLN954w9QsIiJSMamMCIFz52IHOp45wx9r1pgdR0REKhiVEaFO9+7sqFWLBGDdvHlmxxERkQpGZUQuWbCAusATEREkJCSYnUZERCoQlREBoMPtt9OqQwfS0tJ45513zI4jIiIViMqIAGCxWJgwYQIW4I/XX+fC2bNmRxIRkQpCZUSy3Xnnnfzs6srHyclsffRRs+OIiEgFoTIi2ZycnKBXLwD8V6zAlpFhciIREakIVEYkh7YLFpAMNMjIYOu0aWbHERGRCkBlRHKo5OdHVOfOALjNm4dhGCYnEhGR8k5lRHJpumAB6UDbc+fY/d57ZscREZFyTmVEcqnZpg3bGzcGIHX6dHPDiIhIuacyInnye/11AKwnT/L7rl0mpxERkfJMZUTyFHjbbUzo2ZPOwGv/+Y/ZcUREpBy7pjIyf/58AgMDcXV1JTg4mA0bNhRou02bNuHo6Eiby1eJlVLtrlmzAPjoo4+Ij483OY2IiJRXhS4jK1euZPz48UydOpWdO3cSFhZGnz59iI2Nvep2SUlJDBs2jBtvvPGaw0rJCg0NpUuXLrimp7NmwgSz44iISDllMQr53c1OnTrRrl07FixYkL2uWbNm9O/fn9mzZ+e73eDBg2nUqBFWq5XVq1cTFRVV4OdMTk7G29ubpKQkvLy8ChNXrtP6d96hzbhxWAD7oUNUCQw0O5KIiJQRBf38LtTMSHp6OpGRkfTKOkvnZb169WLz5s35brdkyRIOHjzIc889V6DnSUtLIzk5Occi5uh+330kurriDfw6ZozZcUREpBwqVBlJTEzEZrPh4+OTY72Pjw8nT57Mc5v9+/fz1FNPsWzZMhwdHQv0PLNnz8bb2zt78ff3L0xMKUIWq5XErBLS+qefOHfqlMmJRESkvLmmA1gtFkuO24Zh5FoHYLPZuOeee5gxYwaNs85bURBTpkwhKSkpezl27Ni1xJQi0uG11zjm6Eh1w2D7Aw+YHUdERMqZgk1VZKlevTpWqzXXLEhCQkKu2RKAc+fOsWPHDnbu3MnDDz8MgN1uxzAMHB0d+f7777nhhhtybefi4oKLi0thokkxsrq4EHv33fh/9BFNv/6ai8nJuOrYHRERKSKFmhlxdnYmODiY8PDwHOvDw8MJDQ3NNd7Ly4s9e/YQFRWVvYwdO5YmTZoQFRVFp06dri+9lJgO8+ZxysEBP7udLQ89ZHYcEREpRwo1MwIwceJEhg4dSvv27QkJCeG9994jNjaWsWPHApd2sRw/fpwPP/wQBwcHgoKCcmxfs2ZNXF1dc62X0s3Z05M/br+daqtXs/frr+mamVngY4BERESuptCfJoMGDeLMmTPMnDmT+Ph4goKCWLt2LQEBAQDEx8f/4zlHpGwKfu89QiIi2PHnn3ivWMGQIUPMjiQiIuVAoc8zYgadZ6T0eOmll5gyZQrNmjVjz549WK1WsyOJiEgpVSznGRF58MEHqVy5MkZMDD9d5SR3IiIiBaUyIoXi5eXFwl692Av4z5qFLSPD7EgiIlLGqYxIod388sukAE3S0vjlqafMjiMiImWcyogUmle9ekR17w5AjXnzsGdmmpxIRETKMpURuSatlywhmazZkSefNDuOiIiUYSojck28AwOJ6tYNgOqaHRERkeugMiLXrPUHH/xvdkTHjoiIyDVSGZFrdnl2JBH4Zvly7Ha72ZFERKQMUhmR69Lq449p4+XFSydO8Nlnn5kdR0REyiCVEbkulf39uX/SJACmTZtGpo4dERGRQlIZkes2fvx4qlWtStM//uCnrGIiIiJSUCojct28vLxYduONrAGazJtH2rlzZkcSEZEyRGVEikS3+fNJcHCgrs3G5jFjzI4jIiJliMqIFAm36tX54667AGj2+eekJiaanEhERMoKlREpMp0WLiTO0RFfu52tw4aZHUdERMoIlREpMk6VKhE7ciQAbdat42xsrMmJRESkLFAZkSLV6T//4ZCzM1UNg8ghQ8yOIyIiZYDKiBQpq4sLiY89xi7gP9u2ER8fb3YkEREp5VRGpMh1eOklxnXqxJq0NKZPn252HBERKeVURqTIWRwceOW11wBYuHAh0dHRJicSEZHSTGVEikXXrl0Z2K8fk+12TvTubXYcEREpxVRGpNi89MADvAjcdOwYO99+2+w4IiJSSqmMSLEJ7NuXX5o1A8D61FPYbTaTE4mISGmkMiLFqsmKFaQCrVJT2Tx5stlxRESkFFIZkWJVo1UrInv0AKD222/rInoiIpKLyogUu+Dly0lwcCAwM5NNQ4eaHUdEREoZlREpdh6+vuy/914Agr76ioRjx0xOJCIipYnKiJSIzgsXsrxGDcIMg2eef97sOCIiUoqojEiJsDo7479qFX9w6URoO3fuNDuSiIiUEiojUmK6du3K4MGDMQyD1++/H8NuNzuSiIiUAiojUqJefvll5lmtfLhjB5uffNLsOCIiUgqojEiJqlu3LkFduuAA+L/xBhf++svsSCIiYjKVESlx7T//nHirlbo2G78MHGh2HBERMZnKiJQ49xo1ODpuHAAd//tf4rZsMTmRiIiYSWVETNHpzTfZ4+VFJSB2wACz44iIiIlURsQUFgcHXBcvJhMIPXGCrTNmmB1JRERMojIipmk0YACbg4M5AyydN4/z58+bHUlEREygMiKmavf11/SsVYsFp08ze/Zss+OIiIgJVEbEVJVq1WL6228D8Morr7Bv3z6TE4mISElTGRHT/etf/6JP797cmp7OoRtu0JlZRUQqGJURMZ3FYmH+1KmsAPqcOMGmSZPMjiQiIiVIZURKhXpdu7K1e3cAGr31Fn8dOmRyIhERKSkqI1JqdFqzhsPOzvjY7ezp3dvsOCIiUkJURqTUcPH2JnXuXAC67d/PjldeMTeQiIiUCJURKVWCxo1jQ8uWANSYOpWUhASTE4mISHFTGZFSp+26dZywWgnIzOTL4cPNjiMiIsVMZURKnUp+fhx//nn6AMO/+46tW7eaHUlERIqRyoiUSh2mTKHG0KEYhsHo0aNJT083O5KIiBQTlREptd544w1q1KjB2b17WX3XXWbHERGRYqIyIqVWtWrVWDprFr8B//7qK/YsWmR2JBERKQYqI1Kq9bn/fmICAnAEPMaNIzUx0exIIiJSxFRGpNRr/tNPnHRwoH5GBttuusnsOCIiUsRURqTU8w4M5MSMGQD03LWLX19/3eREIiJSlFRGpExo98wzbGjeHIAaTz5J0rFjJicSEZGiojIiZUa7H34g1tERf5uNn2691ew4IiJSRFRGpMzw8PUl6c03eRsY/NtvfPzxx2ZHEhGRIqAyImVKywcf5Mz06aQB48aN4+DBg2ZHEhGR66QyImXOM888Q7du3TifksKKG28kPTXV7EgiInIdrqmMzJ8/n8DAQFxdXQkODmbDhg35jl21ahU333wzNWrUwMvLi5CQEL777rtrDixitVr5+OOPWe3kxNSjR9l4441mRxIRketQ6DKycuVKxo8fz9SpU9m5cydhYWH06dOH2NjYPMdHRERw8803s3btWiIjI+nZsyf9+vVj586d1x1eKi5/f398JkwAoMfWrex46SWTE4mIyLWyGIZhFGaDTp060a5dOxYsWJC9rlmzZvTv35/Zs2cX6DFatGjBoEGDmDZtWoHGJycn4+3tTVJSEl5eXoWJK+XchhYtCIuO5pSDA5Zdu6gZFGR2JBERyVLQz+9CzYykp6cTGRlJr169cqzv1asXmzdvLtBj2O12zp07R9WqVfMdk5aWRnJyco5FJC/tIyI44OKCj91OXFgYmRcvmh1JREQKqVBlJDExEZvNho+PT471Pj4+nDx5skCP8frrr5OamsrAgQPzHTN79my8vb2zF39//8LElArErVo1rKtWkQK0O3uWDT16mB1JREQK6ZoOYLVYLDluG4aRa11eli9fzvTp01m5ciU1a9bMd9yUKVNISkrKXo7pbJtyFYG33spvWcePdNy6lW91dV8RkTKlUGWkevXqWK3WXLMgCQkJuWZL/m7lypWMHj2aTz/9lJv+4WJnLi4ueHl55VhErqbznDl8ExZGCDBowgT++OMPsyOJiEgBFaqMODs7ExwcTHh4eI714eHhhIaG5rvd8uXLGTFiBJ988gl9+/a9tqQi/+CWH36gclgY586dY8CAAaTq/CMiImVCoXfTTJw4kYULF7J48WJiYmKYMGECsbGxjB07Fri0i2XYsGHZ45cvX86wYcN4/fXX6dy5MydPnuTkyZMkJSUV3asQAZycnFi5ciW+vr5U+u03fmrfHsNuNzuWiIj8g0KXkUGDBjF37lxmzpxJmzZtiIiIYO3atQQEBAAQHx+f45wj7777LpmZmTz00EPUqlUre3nssceK7lWIZKlVqxarFizgB+C2339nvWbiRERKvUKfZ8QMOs+IFFbE4MF0W7kSO7Dt6afp/MILZkcSEalwiuU8IyJlRbfly9nQogUOQIsXX2TfF1+YHUlERPKhMiLlk8VC523b2FmlCp6A+6BBJMbEmJ1KRETyoDIi5ZaTuzv1tm3jqJMT/jYbcZ07k3bunNmxRETkb1RGpFyr0rAhti+/JAn4IzmZsfffj13fsBERKVVURqTcq9+3L7+99x5DrFY+WLGCp59+2uxIIiJyBZURqRC63Hcf72edJv6Vl1/mi0mTTE4kIiKXOZodQKSkDB8+nONHj1Lvuef415w5bHFyovNLL5kdS0SkwtPMiFQoU555Bv+mTXEAWr/8MnvefdfsSCIiFZ7KiFQoFgcHQnbuZJuPD25AnXHj+OPzz82OJSJSoamMSIXj6OpK0J497PH0pIphUGXgQA59+63ZsUREKiyVEamQ3GvUoO6ePfzu5kYNw8Dttts4FhFhdiwRkQpJZUQqLO+AAGr8+isHXFyoarczbeBAjh8/bnYsEZEKR2VEKrRqTZvi+csvjPbz44NTp7jppptISEgwO5aISIWiMiIVnk/btryweTP+/v78/vvvDA8N5VRUlNmxREQqDJURESAgIIAff/yRLr6+vH/wIKkdO3Ji+3azY4mIVAgqIyJZGjZsyCfLl4PVSv2MDDJCQ4nbtMnsWCIi5Z7KiMgV6vboAT//zFFHRwIyM6F7d2LXrzc7lohIuaYyIvI3dbp0wXnzZg47OVHHZsPpxhvZ/+WXZscSESm3VEZE8lCrQwc8tm/ngIsLtex2agwYwM6sC+2JiEjRUhkRyUfN1q2ptncvuz09OWwY3PLQQ6xZs8bsWCIi5Y7KiMhVVGnQgIaHD/P6zTdzOi2Nf//737z//vtmxxIRKVdURkT+gXu1anywdi2jRo3Cbrez7/77WR8Whj0z0+xoIiLlgsqISAE4OjqycOFC/jNmDK8BPTZuZFu9epw/fdrsaCIiZZ7KiEgBWSwWHn7/fTaMGkU60Pn4cWIDAojXydFERK6LyohIIYUtWkTM22+TaLHQ9MIFHDp3JvqDD8yOJSJSZqmMiFyD1g89xIX16/nDxQUfu536I0eyYfRos2OJiJRJKiMi18i/Wzd8Dxxgq48PrsDixYsZM2YMFy5cMDuaiEiZojIich286tShQ1wcK4YMYanFwqJFiwgNDeXg/v1mRxMRKTNURkSuk4OjI4M/+ojvv/+eGjVqEBsVxcWmTdny+ONmRxMRKRNURkSKyE033cSvv/7K3Nq1aWG30/m119jYuDEp8fFmRxMRKdVURkSKUJ06dRj8xx/81LkzdqDr/v2cCQhgr65rIyKSL5URkSLm5O5Oz19+YdfcucRZrQRkZNBkzBh+7tmTTB3cKiKSi8qISDFp+9hjVDpwgI116+IIdF+/noUNGhAVFWV2NBGRUkVlRKQYVa5Xjy5HjrDxwQf5zWrlyfh42rdvz9SpU7l48aLZ8URESgWVEZFiZrFY6DpvHtVjY+l1553YbDZefPFFVvv6svutt8yOJyJiOpURkRLi6+fHZ599xhdffMGDXl4MTkqi1WOP8Uu9epyKjDQ7noiIaVRGRErYv//9b17YtYv1zZphB0KOHsWjfXsi+vYlIyXF7HgiIiVOZUTEBJXr1aNHdDQxH37IrkqVqAR0W7uWuGrV2D51KobdbnZEEZESozIiYqIWQ4fS8uxZfh49mgSLhcD0dBxffJGePXqwZcsWs+OJiJQIlRERkzlYrXRfuBDnI0f4sXNnpjo58fOGDYSEhHBv//4c/r//MzuiiEixUhkRKSUq163LDb/8wrsHDzJy5EgcHByot2YNAbfdxi8BARz66iuzI4qIFAuVEZFSxt/fn8WLF7Nr1y56+vvjAITExlL/jjvY6u/PgS++MDuiiEiRUhkRKaWCgoK4KTaWmBUr2OTnhx3oFBdHwzvvJLJmTaLmzMEwDLNjiohcN5URkVKu2aBBdDl+nD8+/5yN/v7YgeDTp/l10iTat2/PsmXLyMjIMDumiMg1UxkRKSOaDhhA19hYjoaH81NQEPOcnfn1118ZMmQIt/n5EREWxonNm82OKSJSaCojImVM4E030XPPHr47fpznn38eHx8fBiYm0m3jRny7dGGHjw87nn0WW1qa2VFFRArEYpSBnc7Jycl4e3uTlJSEl5eX2XFESpW0tDS2PvMMbosW0eGvv7LXJzg48Hvbtvg99RQNBwwAi8XElCJSERX081szIyJlnIuLC91efZUOf/7JgW+/5cfgYE5bLNS02+kWGYn1rrsIbteOOXPmcOzYMbPjiojkopkRkXIoLSWFX198EfvSpXx38iTPZ51e3hnY4ulJSlgY9Z98ktrdupkbVETKtYJ+fquMiJRziYmJrFy5kpUrV+K9YQNfX3HffldX4oODqT5iBE2HD8fBycm0nCJS/qiMiEgu8b/9xu8vvIDnunW0OXsWxyvuO2Ox8NENN1Bz1ChuvPFGfHx8TMspIuWDyoiIXNXpffuIeeMNHNaupcWxY1QBmgMxWfc/FhDAnVWq4NK7N41GjaJyo0YmphWRskhlREQKLP38eXYvXMinx44R/t//EhUVxSfA3VeMOezszPEGDbB2706dAQOoc8MNWBx0DLyI5E9lRESu2enTp4l+/XXsX31F7YMHaZyenuN+G9CgWjVadOpEx44d6Vm1Kk06daJm+/YqKCKSTWVERIrM6X372P/BB1wMD6favn1cTE2l8xX/6fgF6AyctVg4Urky5+rXx9q6NVXDwqjbuzfuvr6mZRcR86iMiEixSbt4kV27d7Nt2za2btnCxFWraHHhAs55jD0I9Kpfn2bNmtGkSRN6JydTo1EjfMLC8GnfXt/gESnHVEZEpERdOHuWg//3f5wJD8e2cydesbHUSU5mm93OHVeMSwSqZf2cDsQ7OXHGy4vzPj6ktmhB0p13UrduXfz9/fGtUQOrc14VR0TKgmItI/Pnz+fVV18lPj6eFi1aMHfuXMLCwvId//PPPzNx4kT27t2Ln58fTzzxBGPHji3w86mMiJRdCXFx7N2/n3379nEwOpp/f/op1f/6i7rp6bj8bey3wK1X3D4N2Bwc+NPNjRRPT9KqVcPu44NDrVo4NG+O4003UbNmTWrWrIm7m5tOeS9SyhRbGVm5ciVDhw5l/vz5dOnShXfffZeFCxcSHR1N3bp1c40/fPgwQUFB3HfffTzwwANs2rSJBx98kOXLlzNgwIAifTEiUnZkXLzIiR07SNy2jXN79mD/4w/222x87OREbGwsSXFxnM06c2xe/g+47YrbZ4AMBweSnJ1JdXMj3d2djEqVsHt7kxwYyIkbbsDb25vKlSvjFxeHR40auNWsiYevLx6+vlhdXYv7JYtUOMVWRjp16kS7du1YsGBB9rpmzZrRv39/Zs+enWv8k08+yVdffUVMTEz2urFjx7Jr1y5++eWXAj2nyohIxWPLzOTU7t38+dtvnNu3j7TDh7HHxeGQkIDL2bP86ujIK1Yrp06dwpKWxoWrPNbXwO1X3D4PuP1tzAUg1WJhi7s7T9evT6VKlfDw8OCZmBhcHBywubhgd3XFyFpwdeW8jw9Hu3bFxcUFFxcXAn77DWdHR6zu7ljd3LC6ueHo7o6jhwcOXl5Y6tfHyckJZ2dnnC5cwNHRESd3d5zc3HBwdtY3kaTcKejnt2O+9+QhPT2dyMhInnrqqRzre/XqxebNm/Pc5pdffqFXr1451t1yyy0sWrSIjIwMnPI4eC0tLY20Ky5/npycXJiYIlIOWB0d8WvXDr927fK8PwR4CDAMg9Rz5zj2228kHz7M+aNHSTt+nMzTp7GfOYPlr7847erK7TVqkJSURMrZs5zauxcPm41KhpFdStwAN8OA1FT27NmT/TyrAY98MkYA45Yty759CqiZz9gdQIcrbh8CAv82JjNriXFw4AYvL6xWK46Ojqz66y/q2WzYAbvFgs1iwZ61JDg780jDhlitVhwcHHj28GECLl7EsFgw4NL/WiwYDg6kODnxYps2ODg44ODgwPDffycgJQWyxpE1FiDd0ZG3O3XCYrFgsVi4/fffqXf2bPb9QI7xC7N21VssFm74/XfqnTmDBTCyxl25zSddumCzWrFYLITu20fgqVPZ9/G38V+GhJCWddxQ8P79NDx5MvebmzX+/zp25LyrKxaLhVYHD9I4Li6ffw34vkMHznlc+pdtfuQIzY4ezXfsj+3a8ZenJwBN4uIIOnQo37ERrVuTWLkyAA2OH6f1gQP5jt0UFMSpapeOoAo4eZLgffvyHbu1eXOO16gBQO2EBDpd8Qf+33dPRjZpQmzWN9h8z5wh5Lff8n3cqIYNOVy7NgDDhg2jXT7/fytuhSojiYmJ2Gy2XKeJ9vHx4WRevyDAyZMn8xyfmZlJYmIitWrVyrXN7NmzmTFjRmGiiUgFZbFYqOTlRaXQUAgNzXNMT2BUHusNw+BiSgop8fGcP3WKC6dPU9tm4/vKlUlNTSUlJYUd//0v9pQUjNRUSE3FcuEClrQ0LGlpnHB351+Bgdl/QB389VcS0tNxtNlwtNtxtNtxNgyc7HaSrVbcnZ3JyMi49IdYHnkcLy92O2fPns1eXxXwy+f1O6alsXPnzuzbtYCW+Yw9BXz33XfZt58AOuYzNgX4LDY2+/ZwIO9395KwKz6g+wL5H0UIvfbt42LWz6HA1S7X+O/oaM5k/Twf6H6VscP27OFy4lf/Yey4PXuyzzY8HbjaUYyTdu8mMuvnx7lUgvMzbdcuIrJ+fgh49CpjX961i2+zfh4BTLzK2P/s2sUXWT/fxaV/u/wsjoriw6yf+wJPXWXsyqgoLu/n6Ny5c9koI5dZ/tbCDMPIte6fxue1/rIpU6YwceL//lmSk5Px9/e/lqgiIvmyWCy4enri6ukJjRvnPWjIkKs+xt1Xvfd/bgBSs342DIPMixe5cPEiGefPY0tLy/5fW1oanobB79Wrk5mZic1mw75/P9Hnz2NPT8dus2FkZGDPyMDIzCTTwYFvmzfHbrdjs9m4sGsXW86dw7DbMex2sNkwMjPBMMi0Wlnavj12ux273c7FXbvYcPbspf8m22xgGGAYGIaB3WLhrdBQjKzbmbt3s/7MmcsvIMdiGAYvd++ePdaydy/rT5/GACyXx/G///ZPCwvDlrVLyjU6mvWnTmWPycEwmBAaSmbWzIj3vn2sP3Ei15jLxnbuzAWXS4dF1zxwgPVZMyN/PxrBAgzr0IHUrJmRuocPsz6vmZGs7e5q147eWbsYAmJjWX+VmZHb27QhLGtmpM7x46zfvz/P5we4uWVL2mbNjPjFx7P+KjMjYc2b06TmpXk3n4QE1kdH5zu2Y9Om1MmaGamemMj6q8yMtGrUiKezZkaaN2+e77jiVqgyUr16daxWa65ZkISEhHwvquXr65vneEdHR6pVq5bnNpf3v4qIlEcWiwUnNzec3Nxwq1Llnzdo1argD96v37UH+5uepWBsD40tNWOLU6GOlnJ2diY4OJjw8PAc68PDwwnNZ3o0JCQk1/jvv/+e9u3b53m8iIiIiFQshT50e+LEiSxcuJDFixcTExPDhAkTiI2NzT5vyJQpUxg2bFj2+LFjx3L06FEmTpxITEwMixcvZtGiRUyePLnoXoWIiIiUWYU+ZmTQoEGcOXOGmTNnEh8fT1BQEGvXriUgIACA+Ph4Yq846CkwMJC1a9cyYcIE5s2bh5+fH2+99VaBzzEiIiIi5ZtOBy8iIiLFoqCf3zrDjoiIiJhKZURERERMpTIiIiIiplIZEREREVOpjIiIiIipVEZERETEVCojIiIiYiqVERERETGVyoiIiIiYqtCngzfD5ZPEJicnm5xERERECury5/Y/ney9TJSRc+fOAeDv729yEhERESmsc+fO4e3tne/9ZeLaNHa7nRMnTuDp6YnFYrnux0tOTsbf359jx47pWjfFTO91ydF7XXL0Xpccvdclq6jfb8MwOHfuHH5+fjg45H9kSJmYGXFwcKBOnTpF/rheXl765S4heq9Ljt7rkqP3uuTovS5ZRfl+X21G5DIdwCoiIiKmUhkRERERU1XIMuLi4sJzzz2Hi4uL2VHKPb3XJUfvdcnRe11y9F6XLLPe7zJxAKuIiIiUXxVyZkRERERKD5URERERMZXKiIiIiJhKZURERERMVSHLyPz58wkMDMTV1ZXg4GA2bNhgdqRyJyIign79+uHn54fFYmH16tVmRyq3Zs+eTYcOHfD09KRmzZr079+fffv2mR2rXFqwYAGtWrXKPiFUSEgI3377rdmxKoTZs2djsVgYP3682VHKnenTp2OxWHIsvr6+JZqhwpWRlStXMn78eKZOncrOnTsJCwujT58+xMbGmh2tXElNTaV169a8/fbbZkcp937++WceeughtmzZQnh4OJmZmfTq1YvU1FSzo5U7derU4aWXXmLHjh3s2LGDG264gTvuuIO9e/eaHa1c2759O++99x6tWrUyO0q51aJFC+Lj47OXPXv2lOjzV7iv9nbq1Il27dqxYMGC7HXNmjWjf//+zJ4928Rk5ZfFYuHLL7+kf//+ZkepEE6fPk3NmjX5+eef6datm9lxyr2qVavy6quvMnr0aLOjlEspKSm0a9eO+fPnM2vWLNq0acPcuXPNjlWuTJ8+ndWrVxMVFWVahgo1M5Kenk5kZCS9evXKsb5Xr15s3rzZpFQiRSspKQm49CEpxcdms7FixQpSU1MJCQkxO0659dBDD9G3b19uuukms6OUa/v378fPz4/AwEAGDx7MoUOHSvT5y8SF8opKYmIiNpsNHx+fHOt9fHw4efKkSalEio5hGEycOJGuXbsSFBRkdpxyac+ePYSEhHDx4kUqVarEl19+SfPmzc2OVS6tWLGCX3/9le3bt5sdpVzr1KkTH374IY0bN+bUqVPMmjWL0NBQ9u7dS7Vq1UokQ4UqI5dZLJYctw3DyLVOpCx6+OGH2b17Nxs3bjQ7SrnVpEkToqKiOHv2LF988QXDhw/n559/ViEpYseOHeOxxx7j+++/x9XV1ew45VqfPn2yf27ZsiUhISE0aNCApUuXMnHixBLJUKHKSPXq1bFarblmQRISEnLNloiUNY888ghfffUVERER1KlTx+w45ZazszMNGzYEoH379mzfvp0333yTd9991+Rk5UtkZCQJCQkEBwdnr7PZbERERPD222+TlpaG1Wo1MWH55eHhQcuWLdm/f3+JPWeFOmbE2dmZ4OBgwsPDc6wPDw8nNDTUpFQi18cwDB5++GFWrVrFjz/+SGBgoNmRKhTDMEhLSzM7Rrlz4403smfPHqKiorKX9u3bc++99xIVFaUiUozS0tKIiYmhVq1aJfacFWpmBGDixIkMHTqU9u3bExISwnvvvUdsbCxjx441O1q5kpKSwoEDB7JvHz58mKioKKpWrUrdunVNTFb+PPTQQ3zyySesWbMGT0/P7Jk/b29v3NzcTE5Xvjz99NP06dMHf39/zp07x4oVK1i/fj3r1q0zO1q54+npmeu4Jw8PD6pVq6bjoYrY5MmT6devH3Xr1iUhIYFZs2aRnJzM8OHDSyxDhSsjgwYN4syZM8ycOZP4+HiCgoJYu3YtAQEBZkcrV3bs2EHPnj2zb1/e7zh8+HA++OADk1KVT5e/pt6jR48c65csWcKIESNKPlA5durUKYYOHUp8fDze3t60atWKdevWcfPNN5sdTeSaxcXFcffdd5OYmEiNGjXo3LkzW7ZsKdHPxQp3nhEREREpXSrUMSMiIiJS+qiMiIiIiKlURkRERMRUKiMiIiJiKpURERERMZXKiIiIiJhKZURERERMpTIiIiIiplIZEREREVOpjIiIiIipVEZERETEVCojIiIiYqr/ByS/Xx8jb/NAAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(1)\n",
    "plt.clf()\n",
    "plt.plot(x, c.value, 'k', label = 'PyFVTool (explicit)')\n",
    "plt.plot(x, c_analytic, 'r--', label = 'analytic')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2586549b-a449-405b-a29f-761968b8f1e2",
   "metadata": {},
   "source": [
    "## Comparison with 2D cylindrical\n",
    "\n",
    "Now, solve the same problem in a 2D cylindrical system (a rod) and compare to the 1D analytic solution (along line near the center of the rod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1419383e-99bd-4090-9f9f-0f999c099c05",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the domain\n",
    "L = 5.0  # domain length\n",
    "N = 100 # number of cells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3960abd8-428f-41b5-9739-394b8a7e9902",
   "metadata": {},
   "outputs": [],
   "source": [
    "meshstruct = pf.CylindricalGrid2D(N, N, L, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2548114d-1559-45f4-a2d1-93ec1f4e8f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "BC = pf.BoundaryConditions(meshstruct) # all Neumann boundary condition structure\n",
    "BC.bottom.a[:] = 0.0 \n",
    "BC.bottom.b[:] = 1.0 \n",
    "BC.bottom.c[:] = 1.0 # bottom boundary\n",
    "BC.top.a[:] = 0.0 \n",
    "BC.top.b[:] = 1.0 \n",
    "BC.top.c[:] = 0.0 # top boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5296df4c-f1aa-4091-81fe-fbd42587a0c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = meshstruct.cellcenters.r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "df0db336-95af-4caa-932e-a8d4cad9f0d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "## define the transfer coeffs\n",
    "D_val = 1.0\n",
    "alfa = pf.CellVariable(meshstruct, 1.0)\n",
    "Dave = pf.FaceVariable(meshstruct, D_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a7b41cf5-ddcc-454a-991e-7c80d7ff7ae4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## define initial values\n",
    "c_old = pf.CellVariable(meshstruct, 0.0, BC) # initial values\n",
    "c = pf.CellVariable(meshstruct, 0.0, BC) # working values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5da03b92-7083-4f5e-b3f7-ac85e2838795",
   "metadata": {},
   "outputs": [],
   "source": [
    "## loop\n",
    "dt = 0.001 # time step\n",
    "final_t = 100*dt\n",
    "for t in np.arange(dt, final_t, dt):\n",
    "    # step 1: calculate divergence term\n",
    "    RHS = pf.divergenceTerm(Dave*pf.gradientTerm(c_old))\n",
    "    # step 2: calculate the new value for internal cells\n",
    "    c = pf.solveExplicitPDE(c_old, dt, RHS)\n",
    "    c_old.update_value(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "935d74b7-370b-485a-a3fc-f98aa29a4564",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQ70lEQVR4nO3dd3wUdf7H8demh4QE6YkECL23oAgKgiKKyok/laYgUhSxoBwW9BTl9FA5FVFBVMDGARaKh4iiHEVQKRJpoUkglGAIQkJLIZnfH0OWhCSQTXZ3Nrvv5+Mxj92d/c7MZ/c4953vfOc7NsMwDEREREQs4md1ASIiIuLbFEZERETEUgojIiIiYimFEREREbGUwoiIiIhYSmFERERELKUwIiIiIpZSGBERERFLBVhdQEnk5uZy6NAhKlasiM1ms7ocERERKQHDMDhx4gTR0dH4+RXf/1EuwsihQ4eIiYmxugwREREphf3791OrVq1i3y8XYaRixYqA+WEiIiIsrkZERERKIj09nZiYGPvveHHKRRjJOzUTERGhMCIiIlLOXGqIhQawioiIiKUURkRERMRSDoeRlStX0qtXL6Kjo7HZbCxYsOCS26xYsYK4uDhCQkKoV68e7733XmlqFRERES/k8JiRU6dO0bp1a+677z7uuOOOS7ZPTEzk5ptvZvjw4Xz22WesXr2akSNHUq1atRJtLyIijsvJySE7O9vqMsTL+fv7ExAQUOZpNxwOIz179qRnz54lbv/ee+9Ru3ZtJk2aBEDTpk1Zv349//73vxVGRERc4OTJkxw4cADDMKwuRXxAhQoViIqKIigoqNT7cPnVND///DM9evQosO7GG29k+vTpZGdnExgY6OoSRER8Rk5ODgcOHKBChQpUq1ZNE0WKyxiGQVZWFkeOHCExMZGGDRtedGKzi3F5GDl8+DA1atQosK5GjRqcPXuW1NRUoqKiCm2TmZlJZmam/XV6erqryxQR8QrZ2dkYhkG1atUIDQ21uhzxcqGhoQQGBrJv3z6ysrIICQkp1X7ccjXNhck8r+uwuMQ+YcIEIiMj7YtmXxURcYx6RMRdStsbUmAfTqjjomrWrMnhw4cLrEtJSSEgIIAqVaoUuc3YsWNJS0uzL/v373d1mSIiImIRl4eRjh07snTp0gLrvv/+e9q3b1/seJHg4GD7bKuadVVERLzB8uXLsdlsHD9+3OXHKunUG57C4TBy8uRJ4uPjiY+PB8xLd+Pj40lKSgLMXo1BgwbZ248YMYJ9+/YxevRoEhISmDFjBtOnT2fMmDHO+QQiIlLuDR48GJvNhs1mIzAwkHr16jFmzBhOnTpV4n3UrVvXvo+8pVatWmRlZVG1alVeeumlIrebMGECVatWpVOnToW2z7/UrVvXSZ/24jZu3Mhdd91FjRo1CAkJoVGjRgwfPpydO3eWan979+5l6NChxMbGEhoaSv369Rk3bhxZWVkF2uT/rBUrVqR58+Y89NBD7Nq1y1kfrVgOh5H169fTtm1b2rZtC8Do0aNp27Ytzz//PADJycn2YAIQGxvL4sWLWb58OW3atOGf//wnkydP1mW9IiJSwE033URycjJ79uzhpZdeYsqUKQ7/4Tp+/HiSk5Pty8aNGwkKCuKee+7ho48+KvJy55kzZzJw4EAWLVpk327t2rUA/PDDD/Z169atc8rnvJhFixZx1VVXkZmZyaxZs0hISODTTz8lMjKS5557rlT73L59O7m5uUybNo2tW7fy5ptv8t577/HMM88Uapv3eX///Xf+9a9/kZCQQOvWrfnxxx/L+tEuzigH0tLSDMBIS0uzuhQREY925swZY9u2bcaZM2esLsUh9957r3HbbbcVWDds2DCjZs2aRm5urlG/fn1j4sSJBd7fvHmzYbPZjN27dxuGYRh16tQx3nzzzSL3v2nTJgMwli9fXmD9ypUrDcDYvHlzgfWJiYkGYGzcuNG+bvny5cYVV1xhBAUFGTVr1jSeeuopIzs72/5+RkaG8cgjjxjVqlUzgoODjauvvtpYu3at/f3//e9/BmAcO3asyBpPnTplVK1a1ejdu3eR7x87dqzE3wVgzJ8/v8j9GIZhvPbaa0ZsbOxFP69hGEZOTo7RtWtXo06dOsbZs2eL3NfF/s2V9Pfbt+9Ns3UrPP44vPOO1ZWIiLiEYRicOnXKksUo46RroaGhZGdnY7PZGDJkCDNnzizw/owZM+jcuTP169e/5L5atmzJFVdcUeQ+rrzySlq0aHHR7Q8ePMjNN9/MFVdcwe+//87UqVOZPn16gVM/Tz75JF999RUff/wxv/32Gw0aNODGG2/kr7/+KtHn/e6770hNTeXJJ58s8v1KlSo55bsASEtLo3Llypds5+fnx6hRo9i3bx8bNmwo0b5L5aJRxUO4rGfk448NAwyjWTPn7ldExCIX/pV68uRJA7BkOXnyZInrvrBn5NdffzWqVKli9OnTxzAMwzh06JDh7+9v/Prrr4ZhGEZWVpZRrVo146OPPrJvU6dOHSMoKMgICwuzL2+99Zb9/alTpxphYWHGiRMnDMMwjBMnThhhYWHGtGnTCtVzYU/BM888YzRu3NjIzc21t3n33XeN8PBwIycnxzh58qQRGBhozJo1y/5+VlaWER0dbbz22muGYVy6Z+TVV181AOOvv/666HdVku+Ci/SM7N6924iIiDA++OCDYj9vfgkJCQZgzJ07t8j9qWekrP72NwgMhG3bzEVERCyzaNEiwsPDCQkJoWPHjnTp0oW3334bgKioKG655RZmzJhhb5uRkcFdd91VYB9PPPGE/SKL+Pj4AhdU9O/fn9zcXObOnQvA3LlzMQyDfv36XbK2hIQEOnbsWGD+lquvvto+9f4ff/xBdnY2V199tf39wMBArrzyShISEkr0+Y0S9iSV9LsoyqFDh7jpppu46667GDZsmEN1uXLuGt8OI5UqQd5U9V98YWkpIiKuUKFCBU6ePGnJUqFCBYdq7datG/Hx8ezYsYOMjAzmzZtH9erV7e8PGzaMOXPmcObMGWbOnEnfvn0LHaNq1ao0aNDAvlSqVMn+XmRkJHfeeaf9FMfMmTO58847SzR9hGEYF53As7gf7KK2K06jRo0Ac8DppZTku7jQoUOH6NatGx07duT9998vUU2APUzFxsaWeBtH+XYYAejTx3z8/HNr6xARcQGbzUZYWJgli6N/SYeFhdGgQQPq1KlT5DxUN998M2FhYUydOpVvv/2WIUOGOPx9DB06lNWrV7No0SJWr17N0KFDS7Rds2bNWLNmTYHeizVr1lCxYkUuv/xyGjRoQFBQED/99JP9/ezsbNavX0/Tpk1LdIwePXpQtWpVXnvttSLfzz8/iaPfxcGDB+natSvt2rVj5syZJZ41NTc3l8mTJxMbG2u/itYVFEZ0qkZEpFzw9/dn8ODBjB07lgYNGtCxY0eH93HttdfSoEEDBg0aRIMGDejSpUuJths5ciT79+/nkUceYfv27SxcuJBx48YxevRo/Pz8CAsL48EHH+SJJ55gyZIlbNu2jeHDh3P69OkSB56wsDA+/PBDvvnmG/72t7/xww8/sHfvXtavX8+TTz7JiBEjSvVdHDp0iK5duxITE8O///1vjhw5wuHDhwvNjg5w9OhRDh8+zJ49e/j666/p3r07a9euZfr06fj7+5foc5SGwohO1YiIlBtDhw4lKyurVL0ieYYMGcKxY8cc2sfll1/O4sWLWbt2La1bt2bEiBEMHTqUf/zjH/Y2r7zyCnfccQcDBw6kXbt27N69m++++47LLrusxMe57bbbWLNmDYGBgQwYMIAmTZrQv39/0tLSCk3aVtLv4vvvv2f37t0sW7aMWrVqERUVZV8u1L17d6KiomjZsiVPP/00TZs2ZdOmTXTr1q3En6E0bEZJR8xYKD09ncjISNLS0pw6NbxhGBw6dIjLvv6aCi+9ZF7mq5lhRaQcy8jIIDExkdjY2FLfQdWTrV69mq5du3LgwIFCd4T3NZ7yXVzs31xJf799umfk1ltvpVatWnwRGAj79yuIiIh4qMzMTHbv3s1zzz1Hnz59fDqIeON34dNhpE6dOgDs2LMHnHALZBERcY3Zs2fTuHFj0tLSih3g6Su88bvw6V/gxo0bA5y/+dDZs7BmjYUViYhIUQYPHkxOTg4bNmzg8ssvt7ocS3njd+HTYSTvmu4dO3ZAdjbUqwdXX62rakRERNxIYQTYtWsXuf7+0KqV+YauqhEREXEbnw4jeRPrZGZmsn///vMToCmMiIiIuI1Ph5GAgAAaNGgAnBs3kjcB2tatOlUjIiLiJj4dRuCCcSOVKsGNN5pvqHdERETELRRGzoUR+xU1eXc91L1qRERE3MLnw0ihy3v/9jcICjJP05Twts8iIlI+dO3alccee6xM+9i7dy82m434+Hin1CQKI4V7RipVgqlTYd06aNLEusJERMRygwcPpnfv3gXWxcTEkJycTIsWLUq931GjRhEXF0dwcDBt2rQp9P7y5cux2WzYbDb8/PyIjIykbdu2PPnkkyQnJ5f6uJ5KYeRcGNm7dy8ZGRnmyiFDoH17cPD21yIi4v38/f2pWbMmAQEBpd6HYRgMGTKEvn37XrTdjh07OHToEOvWreOpp57ihx9+oEWLFmzevLnUx/ZEPh9GqlevTkREBIZh8Mcff1hdjoiIT1qyZAnXXHMNlSpVokqVKtx6660F/pucd2pk3rx5dOvWjQoVKtC6dWt+/vlne5ujR4/Sv39/atWqRYUKFWjZsiWzZ88u9pjjx4+nZcuWhdbHxcXx/PPP88ILL/Dxxx+zcOFCey/F8uXLizxNs3XrVm655RYiIiKoWLEinTt3vuhvyuTJk3nooYeoV6/eRb+X6tWrU7NmTRo1akS/fv1YvXo11apV48EHH7zoduWNz4cRm81WeNwIwPr1Zg/JG29YVJmIiBOdOlX8ktcrXJK2Z86UrK3D5Z1i9OjRrFu3jh9//BE/Pz9uv/12cnNzC7R79tlnGTNmDPHx8TRq1Ij+/ftz9uxZwLx7bFxcHIsWLWLLli3cf//9DBw4kF9//bXIYw4ZMoRt27axbt06+7pNmzaxceNGBg8ezJgxY+jTpw833XQTycnJJCcn06lTp0L7OXjwIF26dCEkJIRly5axYcMGhgwZYq/LmUJDQxkxYgSrV68mJSXF6fu3Sun7mLxIo0aNWLduXcEwkpAAM2dC06YwerR1xYmIOEN4ePHv3XwzfPPN+dfVq8Pp00W3vfZaWL78/Ou6dSE1tXA7w3CovDvuuKPA6+nTp1O9enW2bdtWYGzGmDFjuOWWWwB48cUXad68Obt376ZJkyZcfvnljMl39/VHHnmEJUuW8MUXX9ChQ4dCx6xVqxY33ngjM2fO5IorrgBg5syZXHvttfYei9DQUDIzM6lZs2axtb/77rtERkYyZ84cAgMDgfNDAFyhybnxjHv37qV69eouO447+XzPCFww10ievKtqEhLMSdBERMRl/vjjDwYMGEC9evWIiIggNjYWgKSkpALtWuXdtgOIiooCsPcQ5OTk8PLLL9OqVSuqVKlCeHg433//faF95Dd8+HBmz55NRkYG2dnZzJo1iyFDhjhUe3x8PJ07d7YHEVczzgU9mxeNa1TPCEVcUQMQGQk33QRff23OOfLiixZVJyLiBCdPFv+ev3/B1xfr/ve74G/YvXtLXVJ+vXr1IiYmhg8++IDo6Ghyc3Np0aIFWVlZBdrl/8HP+zHOO5Xz+uuv8+abbzJp0iRatmxJWFgYjz32WKF9XHjc4OBg5s+fT3BwMJmZmYV6aS4lNDTUofZllXBu2om6deu69biupDBCEXON5LnrrvNh5IUXdHWNiJRfYWHWty3G0aNHSUhIYNq0aXTu3BmAn376yeH9rFq1ittuu4177rkHMEPKrl27aNq0abHbBAQEcO+99zJz5kyCg4Pp168fFSpUsL8fFBRETk7ORY/bqlUrPv74Y7Kzs13eO3LmzBnef/99unTpQrVq1Vx6LHfSaRqgYcOGABw5coRjx46df+Nvf4PgYNi+XadqRERc5LLLLqNKlSq8//777N69m2XLljG6FGP1GjRowNKlS1mzZg0JCQk88MADHD58+JLbDRs2jGXLlvHtt98WOkVTt25dNm3axI4dO0hNTSU7O7vQ9g8//DDp6en069eP9evXs2vXLj799NOCp/4vsHv3buLj4zl8+DBnzpwhPj6e+Pj4Qr04KSkpHD58mF27djFnzhyuvvpqUlNTmTp1agm/lfJBYQQIDw8nOjoauKB3JCLCPFUDMHeuBZWJiHg/Pz8/5syZw4YNG2jRogWPP/44EydOdHg/zz33HO3atePGG2+ka9eu1KxZs9CEZUVp2LAhnTp1onHjxoUGug4fPpzGjRvTvn17qlWrxurVqwttX6VKFZYtW8bJkye59tpriYuL44MPPrhoL8mwYcNo27Yt06ZNY+fOnbRt25a2bdty6NChAu0aN25MdHQ0cXFxvPLKK3Tv3p0tW7bQrFmzkn0p5YRO05zTqFEjDh06xM6dOwv+Y7zrLnMQ67mBUiIi4nzdu3dn2wV3SzfyXZFTt27dAq8BKlWqVGBd5cqVWbBgwUWPszz/lUD5jvPnn3/ywAMPFHqvWrVqfP/990Vuk1+rVq347rvvLnrsS9WRX9euXQsdw5upZ+ScYseN9OtnnqYZOdKCqkRExJVSUlJ44403OHjwIPfdd5/V5fgs9YycU+QVNVB4lLmIiHiNGjVqULVqVd5//30uu+wyq8vxWQoj5xQ510h+GRmwbBn07KmrakREvIQvnQrxZDpNc07eaZpdu3YVmn6Y7GyoUwduuQV+/92C6kRERLyXwsg5devWJSAggNOnTxcazUxgIJy79l1X1YiIiDiXwsg5gYGB9nsRFBo3AtCnj/n4+ecO33NBRMTddPpB3MUZ/9YURvK56LiRW26B0FDYswd++83NlYmIlIz/uUH3F5sCXcSZTp+7qWJZZp/VANZ8GjduzKJFi4ruGQkLg1tvhS++MHtH4uLcX6CIyCUEBARQoUIFjhw5QmBgIH4X3ktGxEkMw+D06dOkpKRQqVIlexAuDYWRfIq9vDdPnz7nw8grr+iqGhHxODabjaioKBITE9m3b5/V5YgPqFSpEjVr1izTPhRG8rlkGLn5ZqhQwbxL5aZN0Lq1+4oTESmhoKAgGjZsqFM14nKBgYFl6hHJozCST14Y2bNnD1lZWQQFBRVsUKECfPIJtGgB5y4FFhHxRH5+foSEhFhdhkiJ6GRiPlFRUYSHh5Obm8sff/xRdKM77lAQERERcSKFkXxsNpt98rOL3fpZREREnEdh5AIlCiO//GL2kDzzjJuqEhER8V4KIxcoURhJSYF588zxIxdOHS8iIiIOURi5QInCyI03QmQkHDwIq1e7qTIRERHvpDBygRKFkeBg6N3bfK571YiIiJSJwsgFGjZsCMDRo0c5evRo8Q379jUfv/wScnLcUJmIiIh3Uhi5QFhYGDExMcAleke6d4fKleHPP2HFCjdVJyIi4n0URopQolM1gYHwf/9nPtepGhERkVJTGClCieca6dcPrrhCN80TEREpA00HX4QSh5Hrr4e1a91QkYiIiPdSz0gRNAuriIiI+yiMFCEvjOzevZuzZ89eeoO0NHMCNN0hU0RExGEKI0WIiYkhNDSU7Oxs9u7de/HGhgFt2sC998L337ujPBEREa+iMFIEPz8/+3wjlzxVY7NBr17mc11VIyIi4jCFkWI4NG6kXz/zceFCOHPGhVWJiIh4H4WRYjgURq66CmrXhhMn4NtvXVyZiIiId1EYKYZDYcTP7/z08HPmuLAqERER71OqMDJlyhRiY2MJCQkhLi6OVatWXbT9rFmzaN26NRUqVCAqKor77rvv4vd98QAOX96bF0YWLTJ7SERERKREHA4jc+fO5bHHHuPZZ59l48aNdO7cmZ49e5KUlFRk+59++olBgwYxdOhQtm7dyhdffMG6desYNmxYmYt3pUaNGgFw+PBh0tPTL71Bu3bQoAFkZMCaNS6uTkRExHs4HEbeeOMNhg4dyrBhw2jatCmTJk0iJiaGqVOnFtn+l19+oW7dujz66KPExsZyzTXX8MADD7B+/foyF+9KkZGR1KhRAyhh74jNZs41cuAA3Hiji6sTERHxHg6FkaysLDZs2ECPHj0KrO/RowdriukN6NSpEwcOHGDx4sUYhsGff/7Jl19+yS233FLscTIzM0lPTy+wWMHhUzUdO0J0tAsrEhER8T4OhZHU1FRycnLsPQZ5atSoweHDh4vcplOnTsyaNYu+ffsSFBREzZo1qVSpEm+//Xaxx5kwYQKRkZH2JSYmxpEynaZM08Ln5Di5GhEREe9UqgGsNputwGvDMAqty7Nt2zYeffRRnn/+eTZs2MCSJUtITExkxIgRxe5/7NixpKWl2Zf9+/eXpswyK1UYiY+HG244PxGaiIiIXJRDd+2tWrUq/v7+hXpBUlJSCvWW5JkwYQJXX301TzzxBACtWrUiLCyMzp0789JLLxEVFVVom+DgYIKDgx0pzSVKFUbCw+GHH8DfH1JSoHp1F1UnIiLiHRzqGQkKCiIuLo6lS5cWWL906VI6depU5DanT5/Gz6/gYfz9/QGzR8WT5YWRXbt2kZubW7KNGjSAuDjzNM1XX7mwOhEREe/g8Gma0aNH8+GHHzJjxgwSEhJ4/PHHSUpKsp92GTt2LIMGDbK379WrF/PmzWPq1Kns2bOH1atX8+ijj3LllVcS7eGDPWNjYwkMDOTMmTOOnSrKmx5eE6CJiIhckkOnaQD69u3L0aNHGT9+PMnJybRo0YLFixdTp04dAJKTkwvMOTJ48GBOnDjBO++8w9///ncqVarEddddx6uvvuq8T+EiAQEB1K9fn+3bt7Njxw77Z7ykPn3giSdg1So4eBAuv9y1hYqIiJRjNsPTz5UA6enpREZGkpaWRkREhFuP3bt3bxYuXMjkyZN55JFHSr7hNdfA6tXw+uswerTrChQREfFQJf391r1pLqFJkyYAbN++3bEN+/c3H2fPdnJFIiIi3kVh5BLywojDc43cdRdcfz3cfz94fueTiIiIZRweM+JrSt0zUr26eYmviIiIXJR6Ri4h7/LegwcPckJ34xUREXE6hZFLuOyyy+wTujncOwJw+DBMngyl2VZERMQHKIyUQNOmTYFShpGHH4ZRo+Cjj5xblIiIiJdQGCmBUo8bgYIToGkgq4iISCEKIyVQpjByyy1QsSLs2wc//+zkykRERMo/hZESyAsjCQkJjm8cGgq9e5vP//Mf5xUlIiLiJRRGSiBvzMju3bvJzs52fAd5E6B98QWcPevEykRERMo/hZESqFWrFhUqVCA7O5vExETHd9C9O1SpAikpsGyZ8wsUEREpxxRGSsDPz88+30ipxo0EBpozsgYG6hJfERGRCyiMlFCZBrECPPcc/PknPPqoE6sSEREp/zQdfAmVaRArQHS0E6sRERHxHuoZKaEyTXx2oSNHyr4PERERL6EwUkL5T9MYpZ28LC0NrroKYmLg2DEnViciIlJ+KYyUUMOGDbHZbBw/fpyUlJTS7SQyEk6dgsxMmDfPuQWKiIiUUwojJRQSEkJsbCxQhnEjAAMGmI+zZjmhKhERkfJPYcQBZb6iBs5PgLZ8ORw8WPaiREREyjmFEQc4ZRBr3bpw9dXmTfPmzHFOYSIiIuWYwogDnNIzAnD33eaj7lUjIiKiMOKIMs81kueuuyAgAH77TTOyioiIz9OkZw7ICyNJSUmcOnWKsLCw0u2oalUYNw4aNoQ6dZxYoYiISPmjnhEHVK1alSpVqgCwc+fOsu3sH/+Avn0hNNQJlYmIiJRfCiMOcupMrCIiIqIw4iinDWIF89Lel1+GV18t+75ERETKKYURBzltECvAli3m6ZqJEyE7u+z7ExERKYcURhzk1J6R66+HGjXg6FH47ruy709ERKQcUhhxUN6YkZ07d5KTk1O2nQUEnJ+R9bPPyliZiIhI+aQw4qA6deoQHBxMZmYmiYmJZd/hPfeYjwsXmnf1FRER8TEKIw7y9/d37riRdu2gaVPIyNCdfEVExCcpjJRCs2bNANi2bVvZd2azne8d0akaERHxQQojpeDUMAIwYACEh0NMDJR1HIqIiEg5o+ngS8HpYaRuXThyBEJCnLM/ERGRckQ9I6WQF0YSEhLIzc11zk4VRERExEcpjJRC/fr1CQwM5NSpU+zfv995OzYM+P13SEpy3j5FREQ8nMJIKQQGBtKoUSPAiadqAEaPhjZt4O23nbdPERERD6cwUkpOHzcCcO215uOsWRrIKiIiPkNhpJRcEkZ69oTKlSE5GX780Xn7FRER8WAKI6XkkjASHAz9+pnPP/nEefsVERHxYAojpZQ/jBiG4bwd33uv+ThvHpw44bz9ioiIeCiFkVJq2LAh/v7+pKenc+jQIeft+IoroHFjOHMGvvzSefsVERHxUAojpRQcHEyDBg0AJ5+qsdlg0CDz+ddfO2+/IiIiHkphpAxcMm4EYPBg+OYb+OIL5+5XRETEAymMlEHz5s0BF4SR6Gi4+WYI0Gz9IiLi/RRGysBlPSP55eSYM7OKiIh4KYWRMsgLI1u3bnXuFTV5/vlP8yZ6v/7q/H2LiIh4CIWRMmjUqBF+fn4cO3aMlJQU5x9g5044cEBzjoiIiFdTGCmD0NBQ6tWrB7joVE3enCNz5kBmpvP3LyIi4gEURsrIpeNGunWDyy+HY8dg0SLn719ERMQDKIyUkUvDiL8/DBxoPv/oI+fvX0RExAMojJSRy6+oGTzYfPz2Wzh82DXHEBERsZDCSBm5PIw0bgwdO5qX+H72mWuOISIiYiHNqlVGTZo0ASAlJYXU1FSqVq3q/IM89BC0aQM33OD8fYuIiFhMPSNlFBYWRt26dQFISEhwzUHuvhumTIHWrV2zfxEREQspjDiBW2ZiFRER8VIKI06QfyZWlzEMWL3aPGWTkeG644iIiLiZxow4gctumJefYUD//rB/P3TpAn37uu5YIiIiblSqnpEpU6YQGxtLSEgIcXFxrFq16qLtMzMzefbZZ6lTpw7BwcHUr1+fGTNmlKpgT9SiRQsANm/e7LqD+Pmdn5F15kzXHUdERMTNHA4jc+fO5bHHHuPZZ59l48aNdO7cmZ49e5KUlFTsNn369OHHH39k+vTp7Nixg9mzZ9uvQvEGzZo1w2azkZKS4pp71OTJm3Nk6VI4eNB1xxEREXEjm+Hg7WY7dOhAu3btmDp1qn1d06ZN6d27NxMmTCjUfsmSJfTr1489e/ZQuXLlUhWZnp5OZGQkaWlpRERElGofrtawYUN2797Njz/+yHXXXee6A3XpAqtWwb/+BWPHuu44IiIiZVTS32+HekaysrLYsGEDPXr0KLC+R48erFmzpshtvv76a9q3b89rr73G5ZdfTqNGjRgzZgxnzpwp9jiZmZmkp6cXWDxd3qmaLVu2uPZAeb0jM2ea40hERETKOYfCSGpqKjk5OdSoUaPA+ho1anC4mKnK9+zZw08//cSWLVuYP38+kyZN4ssvv+Shhx4q9jgTJkwgMjLSvsTExDhSpiVatmwJuHjcCMBdd0FYGOzaZV5dIyIiUs6VagCrzWYr8NowjELr8uTm5mKz2Zg1axZXXnklN998M2+88QYfffRRsb0jY8eOJS0tzb7s37+/NGW6ldt6RipWhD59oE4d+Osv1x5LRETEDRy6tLdq1ar4+/sX6gVJSUkp1FuSJyoqissvv5zIyEj7uqZNm2IYBgcOHKBhw4aFtgkODiY4ONiR0iyX1zOyZcsWcnNz8fNz4RQub75phhJXHkNERMRNHPo1CwoKIi4ujqVLlxZYv3TpUjp16lTkNldffTWHDh3i5MmT9nU7d+7Ez8+PWrVqlaJkz9SgQQOCgoI4efLkRa8scorISAURERHxGg7/oo0ePZoPP/yQGTNmkJCQwOOPP05SUhIjRowAzFMsgwYNsrcfMGAAVapU4b777mPbtm2sXLmSJ554giFDhhAaGuq8T2KxwMBA++XKLh83kicry7zMV0REpBxzOIz07duXSZMmMX78eNq0acPKlStZvHgxderUASA5OblAz0B4eDhLly7l+PHjtG/fnrvvvptevXoxefJk530KD5H/VI3LZWRAbCz06AHuOJ6IiIiLlGo6+JEjRzJy5Mgi3/voo48KrWvSpEmhUzveyC0zseYJCYEOHWD+fJg+3RxHIiIiUg5p4IETubVnBGDoUPPx008hM9M9xxQREXEyhREnyusZ2b59O9nZ2a4/4I03QnQ0HD0KX3/t+uOJiIi4gMKIE9WuXZuKFSuSnZ3Nzp07XX/AgIDzM7JOn+7644mIiLiAwogT2Ww2944bARgyxHz8/nvYt889xxQREXEihREnc/u4kfr1oWtX8z41333nnmOKiIg4UamuppHiub1nBOC118yra84FIRERkfJEYcTJ3HaPmvyuuMJ9xxIREXEynaZxsrwwsmfPngJT4LvNqVPuP6aIiEgZKIw4WbVq1ew3Ddy2bZv7DpydDffcAzVqwMGD7juuiIhIGSmMuEDeIFa3jhsJDISkJLNnZOZM9x1XRESkjBRGXMCScSMA999vPn74IeTkuPfYIiIipaQw4gJuv7w3zx13wGWXmfON+MC9gERExDsojLiAJZf3AoSGwqBB5vP333fvsUVEREpJYcQFmjVrBsCff/7JkSNH3Hvw4cPNx6+/huRk9x5bRESkFBRGXCA8PJx69eoBsGnTJvcevHlzuPpqc8zIjBnuPbaIiEgpKIy4SJs2bQD4/fff3X/wJ5+EV16BYcPcf2wREREHKYy4SOvWrQGLwsjf/gZPPWXOOSIiIuLhFEZcJK9nJD4+3tI6REREPJ3CiIvkhZGEhASysrKsKWLuXOjWzZwMTURExEMpjLhITEwMlSpVIjs7273Twuf3/vuwfLku8xUREY+mMOIiNpvN+lM1Dz5oPn74IVjVOyMiInIJCiMuZOkgVoDbboOoKPjzT1iwwJoaRERELkFhxIUs7xkJDDx/ee+UKdbUICIicgkKIy6Uf64RwzCsKeL++8HfH1asAKvGroiIiFyEwogLNW3alICAAI4dO8b+/futKaJWLejVy3z+3nvW1CAiInIRCiMuFBwcbL9PjaXzjYwcCddcA127WleDiIhIMRRGXMzSaeHz3HADrFoF//d/1tUgIiJSDIURF8u7okYzsYqIiBRNYcTFLL+iJr+jR+HVV+HXX62uRERExE5hxMXyekb27NlDenq6tcU8+yw8/TS8/rq1dYiIiOSjMOJiVapUoVatWgBs2rTJ2mJGjjQf582DAwesrUVEROQchRE38IhBrACtWkGXLpCTo8t8RUTEYyiMuIFHjRt55BHz8f33ISPD2lpERERQGHELj7qipndvcyK0I0fg88+trkZERERhxB3yeka2bNnC2bNnrS0mIOD83XzfecfaWkRERFAYcYt69eoRHh5ORkYGO3futLocGD4cKlaExo3hzBmrqxERER+nMOIGfn5+tGrVCvCAQawA1apBcjJ8+imEhlpdjYiI+DiFETfxqEGsAGFhVlcgIiICKIy4TV4Y2bhxo7WFXGjLFvj6a6urEBERHxZgdQG+om3btgD89ttvGIaBzWazuCLgf/+D664zT9skJUFIiNUViYiID1LPiJu0bNmSwMBAjh49yr59+6wux9S5M9SubV7mO2uW1dWIiIiPUhhxk+DgYPsg1g0bNlhczTkBAecnQXvzTTAMa+sRERGfpDDiRnFxcQCsX7/e4kryGTbMHMy6dSv88IPV1YiIiA9SGHGj9u3bAx4WRipVgiFDzOdvvmlpKSIi4psURtwoL4xs2LABw5NOiYwaBTYbfPstJCRYXY2IiPgYhRE3at68OUFBQRw7dozExESryzmvfn3429/Mq2r++MPqakRExMcojLhRUFCQ/aZ5HnWqBmDqVPPy3ltvtboSERHxMQojbpb/VI1HiYrSPCMiImIJhRE388gravLLzTVnZM3MtLoSERHxEQojbuaxg1jz9OgBt90Gn31mdSUiIuIjFEbcrFmzZoSEhJCWlsYfnjhY9KabzMd//9vsJREREXExhRE3CwwM9NxBrAD33w8REbB9O3zzjdXViIiID1AYsYBHTn6WJyICRowwn7/2mrW1iIiIT1AYsYDHXlGTZ9QoCAyEn36CX36xuhoREfFyCiMWyLuiZsOGDeR64riM6Gi45x7z+cSJ1tYiIiJeT2HEAk2bNiU0NJQTJ06wa9cuq8sp2pgx5uOBA7rMV0REXEphxAIBAQG0bdsW8OBTNc2aQXy8eZomONjqakRExIuVKoxMmTKF2NhYQkJCiIuLY9WqVSXabvXq1QQEBNCmTZvSHNarePzkZwCtW5s30BMREXEhh8PI3Llzeeyxx3j22WfZuHEjnTt3pmfPniQlJV10u7S0NAYNGsT1119f6mK9iUdfUXOh48dh+XKrqxARES9lMxycBrRDhw60a9eOqVOn2tc1bdqU3r17M2HChGK369evHw0bNsTf358FCxYQHx9f4mOmp6cTGRlJWloaERERjpTrsbZt20bz5s0JDw/n+PHj+Pv7W11S0bZvhw4dzOf79kGlSpaWIyIi5UdJf78d6hnJyspiw4YN9OjRo8D6Hj16sGbNmmK3mzlzJn/88Qfjxo0r0XEyMzNJT08vsHibxo0bExYWxsmTJ9m5c6fV5RSvUSOIiYH0dJgyxepqRETECzkURlJTU8nJyaFGjRoF1teoUYPDhw8Xuc2uXbt4+umnmTVrFgEBASU6zoQJE4iMjLQvMTExjpRZLvj7+9sHsa5du9biai7Czw/GjjWfv/kmnD5tbT0iIuJ1SjWA1XbBoEbDMAqtA8jJyWHAgAG8+OKLNGrUqMT7Hzt2LGlpafZl//79pSnT41111VUA/OLpE4v17Qv16kFqKnzwgdXViIiIl3EojFStWhV/f/9CvSApKSmFeksATpw4wfr163n44YcJCAggICCA8ePH8/vvvxMQEMCyZcuKPE5wcDAREREFFm/UsWNHAH7++WeLK7mEgAB46inz+cSJkJVlbT0iIuJVHAojQUFBxMXFsXTp0gLrly5dSqdOnQq1j4iIYPPmzcTHx9uXESNG0LhxY+Lj4+mQNzDSR+X1jGzevJkTJ05YXM0l3HsvREXBwYPw6adWVyMiIl7E4dM0o0eP5sMPP2TGjBkkJCTw+OOPk5SUxIhzN1cbO3YsgwYNMnfu50eLFi0KLNWrVyckJIQWLVoQFhbm3E9TzkRHR1O7dm1yc3NZt26d1eVcXHCwOSurvz/s2WN1NSIi4kVKNqI0n759+3L06FHGjx9PcnIyLVq0YPHixdSpUweA5OTkS845Iud17NiRpKQkfv75Z6677jqry7m4Bx6A3r3N8SMiIiJO4vA8I1bwxnlG8rz11ls89thj3Hrrrfz3v/+1uhwRERGncck8I+J8eYNYf/nlF8pBLjxv2zbzvjUiIiJlpDBisTZt2hAcHExqaiq7d++2upyS+fxzaNEC7r8fcnOtrkZERMo5hRGL5V2hBOXgEt883btDxYqweTPMn291NSIiUs4pjHiA/KdqyoXKlWHUKPP5+PHqHRERkTJRGPEA5Wbys/weewwiImDTJliwwOpqRESkHFMY8QB5YWTTpk2cPHnS4mpKqHJlePRR8/mLL6p3RERESk1hxANER0cTExNTPiY/y+/xx82xI5s2wcKFVlcjIiLllMKIhyiXp2ryxo5UqQKePp29iIh4LIURD1HuBrHmefJJ2LsXzt0CQERExFEOTwcvrpG/Z8QwDGw2m8UVlVDFilZXICIi5Zx6RjxEmzZtCAoKIjU1lT/++MPqchxnGOZVNYsWWV2JiIiUMwojHiI4OLj8TX6W34wZcPvt5hU2WVlWVyMiIuWIwogHKZeDWPP06wc1akBiIkyfbnU1IiJSjiiMeJByHUbCwuAf/zCf//OfcPq0tfWIiEi5oTDiQTp16gSYk58dP37c2mJKY/hwqFMHkpNhyhSrqxERkXJCYcSDREdH07BhQ3Jzc/npp5+sLsdxwcHwwgvm8wkTID3d0nJERKR8UBjxMNdeey0AK1assLiSUrrnHmjSBP76C954w+pqRESkHFAY8TDlPowEBJhjRlq3hnNjYERERC7GZhiGYXURl5Kenk5kZCRpaWlERERYXY5L7d+/n9q1a+Pv78+xY8eoWB4nFTMMc/FT1hUR8WUl/f3Wr4WHiYmJITY2lpycHFavXm11OaVjsymIiIhIiekXwwOV+1M1ec6cgVdfhfvus7oSERHxYAojHshrwsjevfDMM/DRR1Bee3lERMTlFEY8UF4YWbduHadOnbK4mjJo2hSGDTOfP/GEOY5ERETkAgojHqhu3brUrl2bs2fPsmbNGqvLKZsXXjBnZ/35Z/jqK6urERERD6Qw4oFsNpv3nKqJioIxY8znY8fqJnoiIlKIwoiH8powAmYYqVEDdu+G996zuhoREfEwCiMeKi+MrF27ljNnzlhcTRmFh8P48ebzf/0LMjOtrUdERDyKwoiHql+/PtHR0WRlZfHLL79YXU7ZDR0Kjz4KK1ea97ARERE5R2HEQ3nVuBEAf3946y1o1MjqSkRExMMojHgwrwojF9qxQ5f6iogIoDDi0fLCyM8//0xGRobF1TjRo4+ac5AsWGB1JSIi4gEURjxY48aNqVGjBpmZmaxdu9bqcpwnMtLsFfn738GbQpaIiJSKwogHs9lsdO3aFYBly5ZZW4wzPf00XH45JCbCm29aXY2IiFhMYcTD3XDDDQB89913FlfiRGFh5g30AF5+GQ4csLYeERGxlMKIh7vxxhsBc76Rv/76y+JqnGjAAOjUCU6dgscft7oaERGxkMKIh6tVqxbNmjUjNzeXH374wepynMdmg6lTzUt+v/wSliyxuiIREbGIwkg5kNc74lWnagBatTKvrKlcGU6etLoaERGxiMJIOXDTTTcBZhgxvG1ujhdfNOccufNOqysRERGLKIyUA507dyYkJISDBw+ybds2q8txrooVoWpVq6sQERELKYyUA6GhofYJ0LzuVE0ew4B582DYMM3MKiLiYxRGygmvHTeS5+BB6N8fpk+Hr76yuhoREXEjhZFyIi+MrFixgtOnT1tcjQvUqmVOhgbwyCNw7Ji19YiIiNsojJQTTZs2JSYmhszMTFauXGl1Oa4xdiw0bgyHD8MTT1hdjYiIuInCSDlhs9m8/1RNSAh88IH5fPp08KYp8EVEpFgKI+WI14cRgM6d4cEHzefDh4M3npISEZECFEbKkeuvvx4/Pz8SEhJISkqyuhzXeeUVcwzJnj2wcKHV1YiIiIspjJQjl112GR06dAC8vHckIgJmzoRvvzWvsBEREa+mMFLO5J+N1at17w7nPquIiHg3hZFyJm/cyA8//EB2drbF1bjJ/v3w2WdWVyEiIi6iMFLOtG/fnurVq5OWlsby5cutLsf19u+Hli1h8GBYt87qakRExAUURsoZf39/brvtNgDmz59vcTVuEBMDPXtCTg4MGgRnzlhdkYiIOJnCSDl0++23A7BgwQJyc3MtrsYN3n0XoqJg+3Z45hmrqxERESdTGCmHrrvuOipWrEhycjJr1661uhzXq1zZnAQNYNIk+N//LC1HREScS2GkHAoODuaWW24BfORUDZinau6/33x+332Qnm5tPSIi4jQKI+XU//3f/wEwb948DMOwuBo3+fe/ITYW9u0zn4uIiFdQGCmnevbsSXBwMLt372br1q1Wl+MeFSvCxx/Dww9r7IiIiBdRGCmnwsPDueGGGwAfOlUD5r1r3n7bvKmeiIh4BYWRcizvqhqfCiP55eSYd/n1lcnfRES8VKnCyJQpU4iNjSUkJIS4uDhWrVpVbNt58+Zxww03UK1aNSIiIujYsaP3T2XuJr169cLPz4+NGzeSmJhodTnud+ed5qDW55+3uhIRESkDh8PI3Llzeeyxx3j22WfZuHEjnTt3pmfPnsXeRXblypXccMMNLF68mA0bNtCtWzd69erFxo0by1y8r6tWrRqdO3cGzDlHfM4995iPr74KP/xgbS0iIlJqNsPBSzE6dOhAu3btmDp1qn1d06ZN6d27NxMmTCjRPpo3b07fvn15voR/0aanpxMZGUlaWhoRERGOlOv1Jk+ezKhRo+jcuTMrV660uhz3GzECpk2DmjUhPh5q1LC6IhEROaekv98O9YxkZWWxYcMGevToUWB9jx49WLNmTYn2kZuby4kTJ6hcuXKxbTIzM0lPTy+wSNF69+4NwE8//URKSoq1xVjhjTegeXM4fBj694ezZ62uSEREHORQGElNTSUnJ4caF/z1WaNGDQ4fPlyifbz++uucOnWKPn36FNtmwoQJREZG2peYmBhHyvQptWvXJi4uDsMwWLhwodXluF+FCvDFFxAWZs7M+o9/WF2RiIg4qFQDWG02W4HXhmEUWleU2bNn88ILLzB37lyqV69ebLuxY8eSlpZmX/bv31+aMn3GHXfcAZjfr09q2hRmzDCfv/02HDpkbT0iIuIQh8JI1apV8ff3L9QLkpKSUqi35EJz585l6NChfP7553Tv3v2ibYODg4mIiCiwSPHuvvtuAP73v/+xb98+i6uxSJ8+5kDWn3+G6GirqxEREQc4FEaCgoKIi4tj6dKlBdYvXbqUTp06Fbvd7NmzGTx4MP/5z3/s91QR56lduzbdunUDYNasWRZXY6Enn4RWrayuQkREHOTwaZrRo0fz4YcfMmPGDBISEnj88cdJSkpixIgRgHmKZdCgQfb2s2fPZtCgQbz++utcddVVHD58mMOHD5OWlua8TyH27/yTTz7xnXvVXMzPP8OYMaDvQkTE4zkcRvr27cukSZMYP348bdq0YeXKlSxevJg6deoAkJycXGDOkWnTpnH27FkeeughoqKi7MuoUaOc9ymEO+64g9DQUHbs2MG6deusLsdaf/4J118Pr78OEydaXY2IiFyCw/OMWEHzjJTM3XffzX/+8x8efvhh3n77bavLsda775o31LPZYOFC6NXL6opERHyOS+YZEc+Wd6pm9uzZZGVlWVyNxUaONCdEMwwYMAC2bLG6IhERKYbCiBe5/vrriYqK4ujRo3z77bdWl2Mtmw0mT4Zu3eDkSbNn5MgRq6sSEZEiKIx4kYCAAPtlvp988onF1XiAwEBzQrT69WHvXrjjDvD1HiMREQ+kMOJlBg4cCMB///tf/vrrL4ur8QBVqsB//wsRERAVpatrREQ8kMKIl2nVqhWtW7cmOzubuXPnWl2OZ2ja1LzUd/ZsCA62uhoREbmAwogXyhvI+umnn1pciQdp1gz8zv1zz82F1autrUdEROwURrzQgAED8Pf35+eff2br1q1Wl+NZzp6Fe+6Bzp1h/nyrqxERERRGvFLNmjW57bbbADTfyIX8/SE8/Pwlv2vWWF2RiIjPUxjxUo8++ihgnqo5duyYxdV4EJsNpkyBW2+FjAzzcdMmq6sSEfFpCiNeqkuXLrRq1YrTp08zY8YMq8vxLAEBMGcOdOwIx47BDTfAjh1WVyUi4rMURryUzWaz946888475OTkWFyRhwkLg8WLoU0bSEkx72WTmGh1VSIiPklhxIsNGDCAypUrs3fvXhYtWmR1OZ6nUiX4/nvzSpsjR2D7dqsrEhHxSQojXiw0NJThw4cDMHnyZIur8VDVqsHSpbBkCfTsaXU1IiI+SWHEy40cORI/Pz+WLVvGFt0srmjR0eY9bPLs3AnJydbVIyLiYxRGvFzt2rW5/fbbAXPsiFzC7t1mMOnaFQ4etLoaERGfoDDiA/IGsn7yySe6zPdS/PzMq2127oQuXWDfPqsrEhHxegojPqBz5860bt2aM2fO8MEHH1hdjmerVw9WrjQf9+wxA8mePVZXJSLi1RRGfIDNZmPUqFEAvP7665w6dcriijxcnTqwYgU0agRJSebU8RpvIyLiMgojPuKee+4hNjaWlJQUpk6danU5nq9WLVi+3Lzs99AhM5D89pvVVYmIeCWFER8RGBjI888/D8Crr77KyZMnLa6oHIiKglWroFMniI2FBg2srkhExCspjPiQe+65hwYNGpCamqora0qqcuXz85BERFhdjYiIV1IY8SEBAQGMGzcOgIkTJ5Kenm5xReVEhQpQvfr51xMnwnPPQW6udTWJiHgRhREf079/fxo3bsxff/2lWVlL4/ff4ckn4aWXoE8f0GBgEZEyUxjxMf7+/rzwwguAeWXN8ePHLa2n3GndGmbOhMBA+Oorc2DrgQNWVyUiUq4pjPigPn360Lx5c44fP86bb75pdTnlz+DBsGwZVK0KGzfCFVfA2rVWVyUiUm4pjPggPz8/XnzxRQDefPNNUlNTLa6oHLrmGjOAtGgBhw+bk6N99JHVVYmIlEsKIz7q9ttvp23btpw4cYKnnnrK6nLKp9hYWL0aevWCzEyrqxERKbcURnyUn5+f/fLeGTNmsHr1aosrKqciImDBAli0yDx9k0dX2oiIlJjCiA/r1KkTQ4cOBeDBBx8kOzvb4orKKT8/uOWW869TU82Brl9+aV1NIiLliMKIj3v11VepUqUKmzdv5u2337a6HO8wcaJ5L5u77oIhQ+DECasrEhHxaAojPq5KlSq8+uqrAIwbN44Duky17F56CcaOBZvNvAy4bVv45RerqxIR8VgKI8J9991Hx44dOXnyJI8//rjV5ZR/gYHwr3+ZN9qrXRv++MO8+uaFF0CnwkREClEYEfz8/Hjvvffw9/fnyy+/ZMmSJVaX5B26dDFnbB0wAHJy4MUXYcIEq6sSEfE4CiMCQKtWrXj00UcBGD58OEePHrW4Ii9RqRLMmgX/+Q+0agWjRlldkYiIx1EYEbvx48fTsGFDDhw4wODBgzEMw+qSvEf//uZsrZGR5mvDgBEjzFM5IiI+TmFE7MLDw/n8888JDg5m0aJFTJo0yeqSvItfvv+7zZoF06ZBt25mUNHAYRHxYQojUkCbNm144403AHjqqadYt26dxRV5qV694MEHzStu5syBJk3g1Vc1k6uI+CSFESnkwQcf5M477yQ7O5u+ffuSlpZmdUneJzISpkyB9euhUyc4dQqeftocV7JggXkaR0TERyiMSCE2m40PPviA2NhYEhMTGTZsmMaPuEq7dvDTT/DJJ1CjBuzcCePHK4yIiE9RGJEiVapUiblz5xIYGMiXX37JSy+9ZHVJ3stmg4EDzSDy7LPmDK5540tOnYJt26ytT0TExRRGpFhXXHEFkydPBuD555/nvffes7giLxcRYc7eev3159dNmgQtWkC/fuYU8yIiXkhhRC5qxIgRjBs3DoCRI0fypW7+5l67d5unbObOhZYt4c47zYnURES8iMKIXNK4ceN48MEHMQyDu+++mx9//NHqknzHzJlm+LjrLvN0zldfQZs20LMnLFtmdXUiIk6hMCKXZLPZePvtt7nrrrvIysqid+/erF+/3uqyfEerVvD557B5szknic0GS5bAZ59ZXZmIiFMojEiJ+Pv78+mnn3L99ddz8uRJunfvzjL9Ze5ezZub08rv2gWPPAKjR59/77ffzMGve/daVp6ISGkpjEiJBQcHM3/+fDp37kxaWho33XQTn+mvc/erXx8mTzYHtuZ5913zTsH16sHNN8PChXD2rHU1iog4QGFEHFKxYkW+//57+vTpQ3Z2NgMHDuTll1/WPCRW690bbrjBHOz67bfm61q1zN6T+HjNWyIiHs1mlINfkfT0dCIjI0lLSyMiIsLqcgTIzc3l6aefZuLEiQAMGzaMKVOmEBgYaHFlPm73bnj/ffjoIzhyxFxXty788UfBe+OIiLhBSX+/9V8nKRU/Pz9ee+013nnnHfz8/Pjwww/p2LEjCQkJVpfm2xo0gNdeg4MH4b//Na/CGT78fBDJzIRrr4VXXjEDioiIB1DPiJTZokWLuPfee/nrr78ICQnh1Vdf5eGHH8ZPf4l7nkWLzJv05WnRAm691Vyuugr8/a2rTUS8jnpGxG1uvfVWNm/ezI033khGRgajRo2iR48e7N+/3+rS5EIdO8K0adC9uxk8tmwxe0muuca8N86iRVZXKCI+SGFEnCI6Oppvv/2WKVOmEBoayo8//kiTJk14/vnnSU9Pt7o8yVOlCtx/PyxdCikpMGuWOXdJpUpw9Kh5pU6eBQvg0Udh3rzz409ERFxAp2nE6Xbu3Ml9993HmjVrAKhWrRrPPfccDzzwAEFBQRZXJ0U6exbWrjV7Tmw2c13//jBnzvk2jRtD587m0qEDNGyoQbEiclEl/f1WGBGXMAyD+fPnM3bsWHbu3AlAvXr1+Pvf/87AgQOpWLGixRXKJS1ZYg6CXbmy8E36/PwgPR3CwszXv/9u9q7Urn0+zIiIz1MYEY+QnZ3NjBkzeOGFFzh8+DAA4eHhDBo0iJEjR9K8eXOLK5QS+esvWL0aVq0yH3Ny4Jdfzr9/1VXw669mIGnTxlxatDBnjW3WzLwjsYj4HIUR8SinTp1ixowZTJkyhe3bt9vXd+rUiTvvvJPbb7+dunXrWlegOMYwzveAGIY5AHbdOsjOLty2fn1z/pM8n39uhpZGjSAmRlfwiHgxhRHxSIZhsGzZMqZMmcLChQvJycmxv9emTRtuv/12unfvTvv27TW+pLzJzISEBNi40Txts3WrubRvD19/bbYxDKhcGY4fN18HBpqndurVg9hYcyzKkCHn95mTo7AiUo65NIxMmTKFiRMnkpycTPPmzZk0aRKdO3cutv2KFSsYPXo0W7duJTo6mieffJIRI0aU+HgKI97p4MGDzJs3j3nz5rFy5Upyc3Pt74WEhNChQwc6d+5Mp06daNWqFdHR0dg0HqH8ycqCvGCZkWEOjN250+wtycoq2LZnT1i8+PzrqlXNbaOjzy9RUeZlyE2aQNeu59vm760REY/gsjAyd+5cBg4cyJQpU7j66quZNm0aH374Idu2baN27dqF2icmJtKiRQuGDx/OAw88wOrVqxk5ciSzZ8/mjjvucOqHkfIrNTWV//73vyxatIiVK1eSmppaqM1ll11Gy5YtadmyJQ0bNiQ2Nta+hIeHW1C1lElODhw6BHv2QGKiOSNsvXpw333m+ydOXHysyc03wzffnH992WUQGmoGmKpVzdd5S+vWcM8959uuXQsVKkDFiucX3cpAxOlcFkY6dOhAu3btmDp1qn1d06ZN6d27NxMmTCjU/qmnnuLrr78uME34iBEj+P333/n5559LdEyFEd9iGAY7duxg1apVrFq1ivXr17Nz584Cp3QuVLlyZaKioqhZs6Z9qVKlCpdddpl9qVSpEuHh4QUWnQryYIYBf/5pBpb8S3Kyuf7KK+GZZ8y2p0+fv7KnKL16nT9VBGZoycgo2CYkBMLD4frrC17SPGCAeelzhQrmEhpqLiEhZngaMOB822++gdxcCA42e3TyluBgc9+xsefbHj9uXpUUGAgBAeainh3xMi4JI1lZWVSoUIEvvviC22+/3b5+1KhRxMfHs2LFikLbdOnShbZt2/LWW2/Z182fP58+ffpw+vTpIm+slpmZSWZmZoEPExMTozDiwzIyMkhISGDz5s1s2bKFPXv2kJiYSGJiIseOHSv1fv39/QkJCSmwBAUFERQURGBgoP0xICCAgIAA/P39CQgIwM/PD39//wKPeYvNZrM/FrcABZ7nvc7/eOH6S60rDW857WUzDCJPniTszBnCzpwhPCOD0IwMKmRmUiEjgz8rV+aXFi0A8MvJYeynnxKSlUVwVhZBF4TcrXXr8sFtt9lfv/buuwSdPVvkcXdffjnv3Hmn/fU/33+fimfOFNk2qXp13ujf3/76uZkzqXLBhIA5Nhu5fn4kV6lSoO1DX31F1ePHMc69n2uz2Z8fDw9nWu/e9rb9v/+e6sePYwBGXrtzj2eCg/nollvsbW/96Sdq/vUX5Gub94OQHRDAZzfdZG97/fr1RB85gpHv30zec8Nm4z89etjXX/3779ROSTHfu6AtwJfdupFzbizQlVu3UvfclXb2H6N8bRdecw1Z5/5oaLtjB/UOHSrq6wXg244dOR0SAkDL3btpeOBAsW2XXnEFJ84F2KaJiTTdt6/YtsvateP4ud+eRklJNE9MLLbtqtatSa1UCYB6Bw/SOv/g7QusadGCP6tUAaDO4cO027Gj2LZrmzXjYLVqANRKSeGKi9wHbEPjxiTVrAlAzaNH6Xjhpfn5xDdoQOLllwMwaNAg2rVrV2zb0ihxZ4LhgIMHDxqAsXr16gLrX375ZaNRo0ZFbtOwYUPj5ZdfLrBu9erVBmAcOnSoyG3GjRtnYP67LLCkpaU5Uq74iOPHjxubNm0yli5danz66afGv//9b2PMmDHGfffdZ/Tu3du49tprjVatWhm1a9c2KleubAQHBxf570uLby4BYFwGRm0wmoFR/4L37wPjITCeAGMcGK+C8RYY086ty9/2GzB+AWMjGFvB2AVGEhh/gvHdBW33g2EUs/x+QdvtF2n7xwVt11+k7eEL2q68SNsTF7T99iJtcy5o++VF2hpghORr+/El2lbJ13bKJdrG5Gs78RJtm+Rr++Il2sbla/vkJdp2ydf2oUu0vSlf28GXaHtHvrZ9LtF2UL62t1yi7Yh8bWfPnu30/z6npaUZcOnf7wBK4cK/pgzDuOhfWEW1L2p9nrFjxzJ69Gj767yeEZGiREZG2seSlFR2djYnT57kzJkzZGRkFFiysrLIzs4mKyvL/jwnJ4ezZ8/al9zcXHJycsjJybE/NwwDwzDIzc0lNzfX/vrCBSjwPO91/scL119qXWnaiOOygbRzC0Ag8Ey+91dfYvv8bafl5uKXm4tfTg5+ubn4573OzcXw8+OZfBMDfpuayg/Z2fjl5mIzDPwMA1tuLn6GwVl/f57J99/Hjfv2sT0jAwwD27nF79y/hxx/f55p0sTeNnHXLlJPnTrfFsznQK6fH8+0bWtvm7pjB9+kpdlf2/L+O37u8ZmrrrK/l7FjB98dPVqobZ4xV11F7rmeEb8dO/gh/+0G8rW1AQ9fdRXZ53rQw3bt4sfk5AL7yv8rMvzKK8k41zNS9Y8/WFZMz4gNuKd9e06f6xm5PDGR/yUlFW54rpa74uK48dz/HnWTkvjfnj1F7hegV5s2XHOuZyTm4EH+t2tXkccHuKFlS9qd6xmJTk5m+UV6Rjo3a0bj6tUBqJGSwvJt24pte2WTJtQ61zNSNTWV5RfpGWnVsCHPnOsZadasWbHtXM2hMFK1alX8/f3tk1flSUlJoUaNGkVuU7NmzSLbBwQEUOXc/wgXCg4OJjg42JHSRBwSGBhoH0siIkXr7UDbng607XHpJnbXO9D2OgfadvOAtl3LWVtXcujGEkFBQcTFxbF06dIC65cuXUqnTp2K3KZjx46F2n///fe0b9++yPEiIiIi4lscvsvV6NGj+fDDD5kxYwYJCQk8/vjjJCUl2ecNGTt2LIMGDbK3HzFiBPv27WP06NEkJCQwY8YMpk+fzpgxY5z3KURERKTccnjMSN++fTl69Cjjx48nOTmZFi1asHjxYurUqQNAcnIySfnOvcXGxrJ48WIef/xx3n33XaKjo5k8eXKJ5xgRERER76bp4EVERMQlSvr77fBpGhERERFnUhgRERERSymMiIiIiKUURkRERMRSCiMiIiJiKYURERERsZTCiIiIiFhKYUREREQspTAiIiIilnJ4Ongr5E0Sm56ebnElIiIiUlJ5v9uXmuy9XISREydOABATE2NxJSIiIuKoEydOEBkZWez75eLeNLm5uRw6dIiKFStis9nKvL/09HRiYmLYv3+/7nXjYvqu3Ufftfvou3Yffdfu5ezv2zAMTpw4QXR0NH5+xY8MKRc9I35+ftSqVcvp+42IiNA/bjfRd+0++q7dR9+1++i7di9nft8X6xHJowGsIiIiYimFEREREbGUT4aR4OBgxo0bR3BwsNWleD191+6j79p99F27j75r97Lq+y4XA1hFRETEe/lkz4iIiIh4DoURERERsZTCiIiIiFhKYUREREQs5ZNhZMqUKcTGxhISEkJcXByrVq2yuiSvs3LlSnr16kV0dDQ2m40FCxZYXZLXmjBhAldccQUVK1akevXq9O7dmx07dlhdlleaOnUqrVq1sk8I1bFjR7799lury/IJEyZMwGaz8dhjj1lditd54YUXsNlsBZaaNWu6tQafCyNz587lscce49lnn2Xjxo107tyZnj17kpSUZHVpXuXUqVO0bt2ad955x+pSvN6KFSt46KGH+OWXX1i6dClnz56lR48enDp1yurSvE6tWrV45ZVXWL9+PevXr+e6667jtttuY+vWrVaX5tXWrVvH+++/T6tWrawuxWs1b96c5ORk+7J582a3Ht/nLu3t0KED7dq1Y+rUqfZ1TZs2pXfv3kyYMMHCyryXzWZj/vz59O7d2+pSfMKRI0eoXr06K1asoEuXLlaX4/UqV67MxIkTGTp0qNWleKWTJ0/Srl07pkyZwksvvUSbNm2YNGmS1WV5lRdeeIEFCxYQHx9vWQ0+1TOSlZXFhg0b6NGjR4H1PXr0YM2aNRZVJeJcaWlpgPkjKa6Tk5PDnDlzOHXqFB07drS6HK/10EMPccstt9C9e3erS/Fqu3btIjo6mtjYWPr168eePXvcevxycaM8Z0lNTSUnJ4caNWoUWF+jRg0OHz5sUVUizmMYBqNHj+aaa66hRYsWVpfjlTZv3kzHjh3JyMggPDyc+fPn06xZM6vL8kpz5szht99+Y926dVaX4tU6dOjAJ598QqNGjfjzzz956aWX6NSpE1u3bqVKlSpuqcGnwkgem81W4LVhGIXWiZRHDz/8MJs2beKnn36yuhSv1bhxY+Lj4zl+/DhfffUV9957LytWrFAgcbL9+/czatQovv/+e0JCQqwux6v17NnT/rxly5Z07NiR+vXr8/HHHzN69Gi31OBTYaRq1ar4+/sX6gVJSUkp1FsiUt488sgjfP3116xcuZJatWpZXY7XCgoKokGDBgC0b9+edevW8dZbbzFt2jSLK/MuGzZsICUlhbi4OPu6nJwcVq5cyTvvvENmZib+/v4WVui9wsLCaNmyJbt27XLbMX1qzEhQUBBxcXEsXbq0wPqlS5fSqVMni6oSKRvDMHj44YeZN28ey5YtIzY21uqSfIphGGRmZlpdhte5/vrr2bx5M/Hx8falffv23H333cTHxyuIuFBmZiYJCQlERUW57Zg+1TMCMHr0aAYOHEj79u3p2LEj77//PklJSYwYMcLq0rzKyZMn2b17t/11YmIi8fHxVK5cmdq1a1tYmfd56KGH+M9//sPChQupWLGivecvMjKS0NBQi6vzLs888ww9e/YkJiaGEydOMGfOHJYvX86SJUusLs3rVKxYsdC4p7CwMKpUqaLxUE42ZswYevXqRe3atUlJSeGll14iPT2de++91201+FwY6du3L0ePHmX8+PEkJyfTokULFi9eTJ06dawuzausX7+ebt262V/nnXe89957+eijjyyqyjvlXabetWvXAutnzpzJ4MGD3V+QF/vzzz8ZOHAgycnJREZG0qpVK5YsWcINN9xgdWkipXbgwAH69+9Pamoq1apV46qrruKXX35x6++iz80zIiIiIp7Fp8aMiIiIiOdRGBERERFLKYyIiIiIpRRGRERExFIKIyIiImIphRERERGxlMKIiIiIWEphRERERCylMCIiIiKWUhgRERERSymMiIiIiKUURkRERMRS/w+3zxTtp4gB0AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(2)\n",
    "plt.clf()\n",
    "plt.plot(r, c.value[1,:], 'k', label='PyFVTool Cyl2D')\n",
    "plt.plot(r, c_analytic, 'r--', label='analytic 1D')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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
}