{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Steady-state heat diffusion in 1D cylindrical mesh\n",
    "\n",
    "MW 230817, 240503\n",
    "\n",
    "This example will evaluate the analytic solution for a simple steady-state heat transfer problem, and use PyFVTool to solve the same problem by the finite-volume method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import exp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# explicity import all required routines from pyfvtool\n",
    "from pyfvtool import CylindricalGrid1D\n",
    "from pyfvtool import CellVariable\n",
    "from pyfvtool import transientTerm, diffusionTerm\n",
    "from pyfvtool import harmonicMean\n",
    "from pyfvtool import constantSourceTerm\n",
    "from pyfvtool import solvePDE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem statement\n",
    "\n",
    "The steady-state heat equation in 1D cylindrical geometry (total cylinder radius $R$), with an internal heat source (strength $S$) and Dirichlet (constant-value) outer boundary. This can model, for example the temperature (temperature $T$ relative to outer temperature) profile in a electrically resistive wire through which a current passes.\n",
    "\n",
    "$$ 0 = k \\nabla^2 T + S $$ \n",
    "\n",
    "which, under the circumstances mentioned above, should have solution\n",
    "\n",
    "$$ T(r) = \\frac{S}{4k} (R^2 - r^2) $$\n",
    "\n",
    "*Do not forget that we are in cylindrical coordinates and that $\\nabla^2$ has the corresponding form.*\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System parameters\n",
    "\n",
    "It is possibly wise to not set all parameters to 1 or even integer values for testing purposes but rather use something else."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "R = 2.1\n",
    "k_val = 3.7 # heat transfer coefficient\n",
    "S_val = 4.2 # source strength\n",
    "T_outer = 0.0  # (outer) boundary temperature"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite-volume solution with PyFVTool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define 1D cylindrical grid with radius $R$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nr = 50\n",
    "Lr = R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = CylindricalGrid1D(Nr, Lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create coefficients of diffusion and source terms in the form of `CellVariable`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = CellVariable(mesh, k_val) # heat transfer coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = CellVariable(mesh, S_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create solution variable `T`, initialize to arbitrary 0.0 value.\n",
    "\n",
    "Subsequently, specify a boundary condition for this variable: outer wall will be kept at 0.0. (Dirichlet boundary condition)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = CellVariable(mesh, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# switch the right (=outer) boundary to Dirichlet: fixed temperature\n",
    "T.BCs.right.a = 0.0\n",
    "T.BCs.right.b = 1.0\n",
    "T.BCs.right.c = T_outer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diffusion term requires the face values of the diffusion coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "k_face = harmonicMean(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation is defined as a list of the different matrix equation terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eqnterms = [-diffusionTerm(k_face),\n",
    "             constantSourceTerm(S)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Go solve. Since no transient term is part of the equation, this will yield directly the steady-state solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "solvePDE(T, eqnterms);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrieve solution temperature profile from solution `CellVariable`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "rnum, Tnum = T.plotprofile()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Comparison with the analytic solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_an = np.linspace(0, R)\n",
    "T_an = (S_val/(4*k_val))*(R**2 - r_an**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG0CAYAAADU2ObLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWwklEQVR4nO3dd3hUZdrH8e+k00uAkJAIwQK4iEJQAYkEgSC4CCJFWAQECyolRCyIFQurLoouggVQ9AVEmhWErCYQiiwlqCg2DCaEZCGgoYdkct4/TmZMmUD6TGZ+n+uaK8w5zznzxHGYm6fct8UwDAMRERERN+Xl7A6IiIiIVCUFOyIiIuLWFOyIiIiIW1OwIyIiIm5NwY6IiIi4NQU7IiIi4tYU7IiIiIhb83F2B5wtLy+PQ4cOUa9ePSwWi7O7IyIiIqVgGAYnTpwgJCQEL6/zj914fLBz6NAhwsLCnN0NERERKYfU1FRCQ0PP28bjg5169eoB5n+s+vXrO7k3IiIiUhrHjx8nLCzM/j1+Ph4f7NimrurXr69gR0REpIYpzRIULVAWERERt6ZgR0RERNyagh0RERFxax6/ZkdERMRqtZKTk+PsbkgRfn5+F9xWXhoKdkRExGMZhkFGRgZ//vmns7siDnh5eREeHo6fn1+F7qNgR0REPJYt0GnWrBm1a9dWclkXYkv6m56ezkUXXVSh90bBjoiIeCSr1WoPdAIDA53dHXGgadOmHDp0iNzcXHx9fct9Hy1QFhERj2Rbo1O7dm0n90RKYpu+slqtFbqPgh0REfFomrpyXZX13mgaq6pYrZCYCOnpEBwM3brB1q1/PY+MBG/v81/jqI2IiIiUiYKdqrB6NUyZAgcP/nXM29sMZmxCQ+HVV2Hw4JKvKdoGFBCJiIiUkYKdyrZ6NQwZAoZR+HjR+ca0NLPdypXmc0fXFGwzeHDpAyLb6ykoEhGpHvo716VpzU5lslrNYMQwsALJQCqQDhwGjgHHgdPAOcMgzzDM9vnXFGM7FhNjBjxDhhQOdOCvgGj16r+OrV4NrVpBz54wcqT5s1Wrwm1s/U1IgGXLzJ8VXAAmIuKRSvt3biUaO3YsFoul2KNfv3707t3b4TXbtm3DYrGwe/duDhw4gMViwcfHh7S0tELt0tPT8fHxwWKxcODAgSr7HaqTRnYqU2KiPRjJBFqX4hLLwYP4AXWBOvk/6xZ8bhjUTU2l3ujRNDYMmgCB8NdPwyAQ8IuJgYED4eOPK3eUSP9aEREpWUmj+UX/zq0CN954I++8806hY1u2bGHo0KH8/vvvtGzZstC5RYsWcdVVV9GpUyd7EBMSEsJ7773H9OnT7e0WL15MixYtSElJqZJ+O4OCncqUnm7/oxWoDeTmP/JKuMQAsvMfR8937zNnzvvS9VJTaRIaSvCxY4QaBi2AFkBo/s8WhkEI4B8TA3l5MGxY5U6biYh4mgKj+cUYBlgs5sj8wIFV8o9Ef39/mjdvXujYwIEDadasGe+++y5PPvmk/fjp06dZvnw5zz//fKH2Y8aM4Z133ikU7Lz77ruMGTOGZ555ptL77CwKdipTcLD9jyHAqQKnDMwAKLfAw/b8bH7bkwV+FnycArIwp8EyMYMi289jmIHUCeBERgbJF+hi09RUWo0YwcWGwcVQ6BFsGHjZPpylDYhAoz8i4pkKjOY7ZBiQmmq2i4qqli75+PgwevRo3n33XZ544gn71u0VK1Zw7tw5/vGPfxRqf/PNN/PGG2+wefNmunfvzubNmzl27BgDBgxQsCMliIw0Rz7S0ooFCRbM/9iF/oNbLNCihflnB9fY2zRpAkeOOHzJPOBPzMDnyJgxHFq8mDTgIJCW/7D9ORs4AhzJzWWHg3sFAOGGwcWpqbQZPZrLDYO/AZcD9WyNiv5r5eOPNR0mIp6pwGh+pbQro88++4y6devan/fr148VK1Ywbtw4XnrpJRISEujZsydgTmENHjyYRo0aFbqHr68vo0aNYtGiRXTv3p1FixYxatSoCmUrdkUKdiqTt7f5JT9kiBkQOApebGyJkl591fzp6Bpbm9dfh9hYhwGRF9DYYqFxaCiX3n47LF7s8OUMzFGgg5gLp/cXeRzAHGHal//4rMi02UWYQc/fgL8ZBn9LTeXyJ56g7qxZmg4TEc9UYDS/UtqVUc+ePZk/f779eZ06dQBo27Yt3bp1Y9GiRfTs2ZP9+/eTmJjIhg0bHN5n/PjxdO3aleeff54VK1awbds2cnNzq6TPTmN4uKysLAMwsrKyKu+mq1YZRmioYZhhgPnw9i78PCzMbHe+awq2WbXKMCwW81Gwje3YqlWGkZtr3qNom4JtmzZ1eC4HjP1gbABjHhiTwegFRrAZJzl8WMBoC8ZIMGaDkQBGVsHXCgszjBUrHPenYL9FRJzgzJkzxg8//GCcOXOmfDcozd+5YWFmu0o2ZswYY+DAgSWeX7hwoVGrVi0jKyvLmDFjhtGqVSsjLy/Pfj45OdkAjKSkJMMwDKNz585GVFSUcfXVVxuGYRhJSUkGYCQnJ1d638vifO9RWb6/tfW8KgweDAcOQHw8LF1q/jx9uvDz5OTCoxqOrinYZvBgc6TENu1lExr61wiKbWQJ/hoVsik4ShQaWuy8D+busT4WC/c2bcqrwH+AQ5hTZInAG8Ak4AYgCDPi+RFYCjwARAENgMuAEYbBS6mpxN95JycvtK3etuVdW+FFpCYpzd+5c+Y4Zcp+2LBheHt7s3TpUhYvXswdd9xx3tIL48aNIyEhgXHjxlVjL6tRVURiNUmVjOxUpdxcw4iPN4ylS82fjv7FUNFRog8/vPC/Vho3NjLAWAvGM2DcAsZFJYwAeYHREYyJYCwF43cw8greLz7ecZ9DQzXyIyJVpsIjOzYX+ju3ClxoZMcwDGP8+PFGo0aNDC8vL+P3338vdK7oyE5OTo5x5MgRIycnxzAM9xvZ0Zqdmsbb+8Kr+gcPNhcPl7Qg2DZK5GgdzZw5f40SnW8d0ZQpBD35JP2AfgVeOhPYnf/YBfwXSAGS8h9z89u1ALoB1wHd3nyTjh98UHwBmXZ+iUhNcKG/c51k/PjxLFy4kOjoaC666KLztvXx8aFJkybV1LPqZzGM862idX/Hjx+nQYMGZGVlUb9+fWd3p3pdKHBwtLA4LMwMiAYONDOElmIX2UFgG7AF2IoZ9BRd+lYP6IE5RdYLaE9+em+LxQzCkpNLv/NLRKQUzp49S3JyMuHh4QQEBDi7O+LA+d6jsnx/K9jx5GCnNM4XENkyh4Lj0Z/lyx3uIjsN7MAMfLb4+bHl3Dn+LPKyTYCemMHPDcClTz2F5emniwdWtteqwiylIuKeFOy4vsoKdlxqgfKmTZsYMGAAISEhWCwWPvroo/O2X716NX369KFp06bUr1+frl27sn79+urprKewTZuNGGH+LDjyc6FF00OHOly8VxvoYbEw3WLhs/vuIxNzyuslzCmxOpjTYSuAe4E2wEUzZzLOMFiJmWDRTgudRUTkAlwq2Dl16hRXXnklc+fOvXBjzOCoT58+rF27ll27dtGzZ08GDBhAUlJSFfdU7Cq6i2zgQLyBTsA0YC3wB7AZmIm5w8sPOJiXxzvAUMxRnyjgRWAvYBTMUuqEgnwiIuLaXHYay2KxsGbNGgYNGlSm6/72t78xfPhwnnjiiVK11zRWNSlpOsxqveDanzMNG7L5jz/4AjMY+rFIkzCgP9BvwAB6ffopdR3cA9BUl4gUomks1+eW01gVlZeXx4kTJ2jcuHGJbbKzszl+/Hihh1SDkqbDSpGnolZMDH2A2ZjZnfdj7urqj1niIhV4Exj06ac0AQYCizEzRgOa6hIR8XBuFezMnj2bU6dOMWzYsBLbzJo1iwYNGtgfYWFh1dhDcehCU10zZhRKhNgauB/4HDOgWQtMDAggHLP+1yfAWKAZ0BuYBxzSVJeIiMdym2Bn2bJlPPXUUyxfvpxmzZqV2G769OlkZWXZH6mpqdXYSynR+db+nGf0p5bFQj+LhX9PmMB+4BvgSaADZlX5LzEDI1ten5cef5z9t95avFKxLaePAh4REbfjFsHO8uXLGT9+PB9++CG9e/c+b1t/f3/q169f6CEuoiI7vwYOxIIZ5DyFGfT8grmIuUt+023AQ5s3cwlwLTAHsNci1lSXiIjbqvHBzrJlyxg7dixLly7lpptucnZ3pCqdb/QnMrJYza9LgAcxg5yDwOt169IL83/6/wJTMUd8egELgT801SUiUkhUVBQxMTEVuseBAwewWCzs2bOnUvpUHi4V7Jw8eZI9e/bY/4MkJyezZ88eUlJSAHMKavTo0fb2y5YtY/To0cyePZsuXbqQkZFBRkYGWVlZjm4v7qCcC51bWCzcd+ed9uKm/wa6Yhbv+gq4E7O46SBg+YsvclpTXSIiZTZ27Nhiu6jDwsJIT0+nffv2zukULhbs7Ny5k44dO9KxY0cAYmNj6dixo30beXp6uj3wAXjzzTfJzc3l/vvvJzg42P6YMmWKU/ovTlaKqS4wg5qJmBmcfwOeB64AcoCPgdvWrSMIGIeZ78e+IV5TXSIiZebt7U3z5s3x8XFiOc7KrE5aE9W4qudyYSVVhs/NPW819+/AeLR2baNVkartl4LxPBgHValdxK1UWtVzJ1m3bp1x3XXXGQ0aNDAaN25s3HTTTcavv/5qGMZfVc1XrVplREVFGbVq1TI6dOhgbN261X59ZmamcdtttxktWrQwatWqZbRv395YunRpodfo0aOHMWXKFMMwDOPpp5822rdvX6wfnTp1Mh5//HHjySefLPR3J2DEx8cXq7BuGIaxd+9eo3///ka9evWMunXrGt27d7f3vaDKqnruUiM7IpWinFNd7S0Wnrv7bn4DEjFHdupgLnR+FLgIs5zFCiB71SpzSktTXSJuxTAMTp065ZSHUcYcv6dOnSI2NpYdO3bw5Zdf4uXlxS233EJeXp69zYwZM5g2bRp79uzhsssuY8SIEeTmmqWYz549S0REBJ999hl79+7l7rvv5vbbb2f79u0OX2/cuHH88MMP7Nixw37s22+/JSkpibFjxzJt2jSGDRvGjTfeSHp6Ounp6XTr1q3YfdLS0rj++usJCAjgq6++YteuXYwbN87erypxwXDIzWlkxwM5GpEJCzOPx8cXOn4CjHfAuL7Iv1YaWyzGRDC+dTBCZFgs5v1sI0oi4pIcjRqcPHmy2OhEdT1OnjxZod/n8OHDBmB899139tGUBQsW2M9///33BmDs27evxHv079/feOCBB+zPC47sGIZh9OvXz7j33nvtz2NiYoyoqCj78zFjxhgDBw4sdM+iIzvTp083wsPDjXPnzl3wd9LIjkh5lWFXV13MBIUbMUd4ZgAtvLw4ZhjMxdzq3h1YgpnQEDBDHtuuLtC6HhGpEvv372fkyJG0bt2a+vXrEx4eDlBobWuHDh3sfw4ODgbg8OHDAFitVp577jk6dOhAYGAgdevWZcOGDYWuL+quu+5i2bJlnD17lpycHJYsWcK4cePK1O89e/YQGRmJr69vma6rCCeuFhJxIttUl6Pjr75qTkVZLIXqdV1isfAs8PTEifzntddYAHwEbMl/xGBOfd2DmeWZ9HRzOmvKlMLTXaGh5muoTpeIy6lduzYnT5502muXxYABAwgLC+Ptt98mJCSEvLw82rdvz7lz5+xtCgYUlvx/xNmmuWbPns0rr7zCnDlzuOKKK6hTpw4xMTGFrnf0mv7+/qxZswZ/f3+ys7O59dZby9TvWrVqlal9ZVCwI1KUbVeXoyBlzhy8Gzem72uv0RczKeEC4C3MXD4vAi8BfYF7P/2Um5Ytw7vo/W3relSYVMTlWCwW6tSp4+xuXNDRo0fZt28fb775JpGRkQBs3ry5TPdITExk4MCBjBo1CjCDoF9++YV27dqVeI2Pjw9jxozhnXfewd/fn9tuu61QkObn54f1AqPXHTp0YPHixeTk5FTb6I6msUQcKeVUVzDwOJAMrAGiMSffvwAGLltGODALOFrw3o62sIuIlEGjRo0IDAzkrbfe4tdff+Wrr74iNja2TPe45JJLiIuLY+vWrezbt4977rmHjIyMC15355138tVXX7Fu3bpiU1itWrXi22+/5aeffiIzM5OcnJxi10+cOJHjx49z2223sXPnTn755Rfef/99fvrppzL1vywU7IiUpAy7unwwExKut1j4BZjWrRuBmBXZHwXCgAnAj7Z7a12PiFSAl5cXH3zwAbt27aJ9+/ZMnTqVl156qUz3ePzxx+nUqRN9+/YlKiqK5s2bF0sI6Mill15Kt27daNOmDddee22hc3fddRdt2rShc+fONG3alC1bthS7PjAwkK+++oqTJ0/So0cPIiIiePvtt6t0lMdiGGXc6+Zmjh8/ToMGDcjKylKdLCkbR+txwsJgzhzIzubsyJEsx6zBtafAZf0wS1X0BixLl4K/v9b1iDjB2bNnSU5OJjw8nICAAGd3p8YwDIO2bdtyzz33lHk0qazO9x6V5ftbIzsi5XW+qa7gYAKAMcBuIAEYCFiAdZjTXR2AhR98wFmVphCRGuLw4cO8/PLLpKWlcccddzi7O6WmBcoiFVHSri7bup60NCyGQQ+gB/Ar8BqwCNgL3PnJJ0wH7gMmAYG26w3DnCKLiTHLXHgXW+YsIlLtgoKCaNKkCW+99RaNGjVydndKTSM7IlWhhGzNlwCvWSwcBF6KjuYi4AjwNGaG5qmYu7oAresREZdjGAZHjhxh5MiRzu5KmSjYEakq5ylM2nDVKqaNHct+4AOgI3Aac31Pa8x8PfZ9CbZ8Pa1aQc+eMHKk+bNVK01ziYiUgoIdkap0gXU9PsBwYBfmdvUozOrr7wDtgCHAzoQE1eESqUIevk/HpVXWe6PdWNqNJc5itZqjM2lphTI1f42Zm+eTAk17A9OBnpiLnO0sFnNtUHKy1vWIlJHVauXnn3+mWbNmBAYGXvgCqXZZWVkcOnSISy65pNjW9LJ8f2uBsoizlFCaogvwscXCXsPgxSuvZOk33/Af4D+YdbieAm4gP+gpuK7H0UJpESmRt7c3DRs2tNeKql27tr2kgjhfXl4eR44coXbt2vj4VCxcUbAj4kznKU3Rfs4c3svOZubIkfwLsyzFZsxRnkjMRc09be3T082fVqsZ+KSnQ3CwuStMIz4iJWrevDnwV3FMcS1eXl5cdNFFFQ5CNY2laSxxBSUFKQkJ5mJkIA34J2YdLluZvh6YIz1R8fFw7JiSE4qUk9VqdVjaQJzLz88PLy/Hy4vL8v2tYEfBjrgyB+t6DmIGPW/zV9AT9be/8dT339Oj6PW2fw2p6KiIuBllUBZxFw7y9YQCc4H9wP2An48PCd9/TxTmWp5tBa9X0VEREQU7Ii6vhHw9oWFhzF21il+XLOE+wA+IB7phFiX93tawaHJCEREPo2BHpCY4T76eMKuV1zFLUdyJ+aH+GLgCGAv8brtHwUXMysQsIh5EwY5ITWGrwzVihPnTtssqOBiAMMx1PN9jJiM0gMXAZUAMcDggQJmYRcQjaYGyFihLTVdCcsIdmIkIv8x/XjcggAfOniUWKPR/uhYxi0gNpAXKIp6khKKjVwP/sViIAyJat+bk2bM8DVyMucDZvslWi5hFxM0p2BFxB+cpOtp71Sp2LFjACswprUxgEuaans8wp7u0iFlE3JmCHRF3cZ5FzJaMDIZgrueZBzTFrKo+AOgDfGO7hxYxi4gbUrkIEXdiW8RcVP4iZh/gXmAkZrHRVzDX9HQExgHP+PkRvHq1MjGLiFvRAmUtUBZPUMIi5mTMRczL85/X8ffnkexsHgBqFbxei5hFxMVogbKIFFbCIuZw4AOLhS3ANZdcwqnsbB4H2gBLyV/PA1rELCI1moIdEU9xnkXM3VatYtubb7IUuAhIBf4BRAHf2tppEbOI1FAKdkQ8yXkWMXv973+MAH4EnsGcxtqEuZ5nMvCH7R5axCwiNYyCHRFPc4FMzLWAxzCDnqFAHvBvzG3rC4G8oCBlYhaRGkULlLVAWcRUwiLmLzHz8uzLf371xRfz+v79XF30ei1iFpFqpAXKIlJ2JSxi7oWZh2c2UK9WLXbs38+1wF3AkYLXaxGziLgoBTsi8pcSFjH7hoURu2oVP73/PqMxd2ktwNy1tQBzqgvQImYRcUkKdkSksPMsYg4+d47FwGbgKsxFy3cB12NmZ7azLWIWEXEByqAsIsVdIBPzdZhV1f8NPA5swQx+HsJc3Fwrvx1WqznKk55uXhsZ+deCaBGRaqKRHREpvchIs3SExYIPMBX4AbgZyAWeB9p7e7P+9Gnt2BIRl6FgR0RKz8Ei5ouAj4E1QAvgN6uVG2+6iZG33kpGwfpaYO70GjJEAY+IVCsFOyJSNiUsYh4UFsa+//s/pkyejBewDGgLvEmRBcygHVsiUq2UZ0d5dkTKp6T1OAkJ7OrZk7uB3flNe2Du2rqk4PXx8Y7XBYmIlEJZvr+1QFlEyqekRczp6UQA24G5wAxgI9ABswxFDOCd305EpDpoGktEKlf+TiwfzMDmO8zEhGeAaZg7ub4v0E41tkSkqrlUsLNp0yYGDBhASEgIFouFjz766ILXbNy4kYiICAICAmjdujVvvPFG1XdUREpWYMcWQGsgDngbqI854tMReHbTJnI+/FA7tkSkyrlUsHPq1CmuvPJK5s6dW6r2ycnJ9O/fn8jISJKSknj00UeZPHkyq1atquKeikiJHOzYsgB3Yo7o/B3IAR5/8kmuHj6c3dqxJSJVzGUXKFssFtasWcOgQYNKbPPwww/zySefsG/fPvuxCRMm8M0337Bt2zaH12RnZ5OdnW1/fvz4ccLCwrRAWaSyrV4NU6ZAwWAmLAzjlVdYduYMk8eM4WheHt7Ag8CTQICtncVijg4lJysJoYg45DGFQLdt20Z0dHShY3379mXnzp3k5OQ4vGbWrFk0aNDA/ggLC6uOrop4nhLKTlhuvZWRoaH8kJfHMMAK/BPoDCTZrlWNLRGpRDU62MnIyCAoKKjQsaCgIHJzc8nMzHR4zfTp08nKyrI/UlNTq6OrIp7JtmNrxAjzp22UJj2dZsByYDXQDHOK6xrMHVu5tuu1Y0tEKkGN33puyV8TYGOblSt63Mbf3x9/f/8q75eInIdtJxZwC9AduBdYBTwBfAq8B7RVjS0RqQQ1emSnefPmZGRkFDp2+PBhfHx8CAwMdFKvROSCiuzYagqsAP4PaIhZZLSjxcKrSUnkrVypHVsiUiE1Otjp2rUrcXFxhY5t2LCBzp074+vr66ReicgFlbBj6x+YeXmigbOGQUxsLL2HDuV37dgSkQpwqWDn5MmT7Nmzhz179gDm1vI9e/aQkpICmOttRo8ebW8/YcIEfv/9d2JjY9m3bx+LFi1i4cKFTJs2zRndF5GyKKHGVmhYGF+sXMm8uXOpbbEQD1wBvAPYt46qxpaIlIFLbT1PSEigZ8+exY6PGTOGd999l7Fjx3LgwAESEhLs5zZu3MjUqVP5/vvvCQkJ4eGHH2bChAmlfk3VxhJxsvPU2Pq1Z0/GAFvzmw7CTE7YpOD1qrEl4pHK8v3tUsGOMyjYEXFRy5bByJFYgdnAY5jJCIMxFy/3trVbutTc7SUiHsVj8uyIiBvL34nlDTyEWWaiLZAO9MFMRHiuQDsRkZIo2BER11Rkx1ZHYBdwT/7pfwFdfH35sWlT84AKiopICRTsiIhrcrBjqzbwBvAREAgk5eTQ6eqreWvCBIyWLbU9XUQcUrAjIq6rhB1bA8PC+Pbtt+nduzdnzpzhnjff5Na0NI4WbKTt6SKSTwuUtUBZxPWVsGMrLyeHV5o2ZXpWFjlACPA+cIPtOhUUFXFbWqAsIu6lhBpbXlu28EBWln3x8iHMXVpPYhYYVUFREQEFOyJSk+UXCrUtXr4TM/HgTKAXZvBTsJ2IeCYFOyJScxXYdl4bM+HgEqAusBG4ClhfsJ12bIl4JAU7IlJzFdmeDjASc5TnSuAIcCMwfd06clesUEFREQ+lBcpaoCxSs61ebe66gr9qZgFngQeAefnPrwOWAWEFr7UFSStXmju/RKTG0AJlEfEcJWxPDwgL4/VVq/jwgw+ob7GwBXNa6/OCjVRQVMQjKNgRkZpv8GA4cMAsCrp0qfkzORkGD2ZoUBC7DYMI4Bjwd8xSEzm2a7VjS8Tt+Ti7AyIilcK2Pb2o9HQuBrYADwOvYpaa2A4sxywsamsnIu5JIzsi4t7yd2L5A3OAVUA9IBHolP+zYDsRcT8KdkTEvRXZsTUY2Am0BzKAnsDLDRtidO9uttf2dBG3o2BHRNybg4KilwFfA//AzLT8wJ9/MmzECE4sWaLt6SJuSMGOiLg/Bzu26gDvh4Yy98478fX1ZeXKlVw9ahQ/HDxY+FoVFBWp8ZRnR3l2RDxHCQVFv96yhaE9enDQaqUOsAC4reB1Kigq4nLK8v2t3Vgi4jlK2LHVJSeH3VYrI4AvgRHANsxdW75QeHu6ox1fIuLSNI0lIpKeTlPMOlqP5h96DbOC+uEi7USk5lGwIyKSv+3cG3gO+Bhze/omoDNmra2C7USkZlGwIyJSZHv6zcB/MXdtpQLdgSWNG5vtRKTGUbAjIuJge3pbzCzL/TGLio46dowHHnqI3Nxc5eIRqWEU7IiIgMPt6Q2BT0JDmXHrrQC8/PLL9IuI4GhYmHLxiNQg2nqureciUlAJ29NXrlzJmFGjOJ2dTTjwEdDBdk3+aBArV5pBk4hUubJ8fyvYUbAjIqVhtfJtSAiDDh8mGagNLAaG2M4rF49ItSrL97emsURESiMxkQ6HD7MDc0v6aWAo8DiQB4Vz8YiIS1GwIyJSGvk5dgKBdUBs/uFngeGYwU/BdiLiOhTsiIiURoEcOz7AbGAhZobllcD1QFqRdiLiGhTsiIiURpFcPADjgP9gjvbsAq728mJnrVrmSW1PF3EZCnZERErDQS4eMEd0/gtcDqTn5XF9z56smDbN3I6u7ekiLkHBjohIaTnIxQPQOiyMbe+/T79+/Thz5gzDZs/mmYMHKbTVNS0NhgxRwCPiBNp6rq3nIlJWJeTisZ47x7TAQOacPAmY1dMXArVs12l7ukilKcv3t0819UlExH14e0NUVPHDW7fyysmTtAPuB5YBv2EmIGwOhbenO7heRKqGprFERCpL/rbzu4ENQCPM+lrXAt87aCci1UPBjohIZSmw7bwnZqBzKZACXAd85aCdiFQ9BTsiIpWlyPb0S4FtQHcgC+gLLG7c2GwnItVGwY6ISGVxsD09EIgDbgNygbHHjvHkzJkYhqFcPCLVRMGOiEhlcrA9PQBYEhrKo/kV0WfOnMmYnj0517KlcvGIVANtPdfWcxGpCiVsT1+wYAET7rkHa14eUcBqzIXMwF/JCleuNIMmESlRWb6/Fewo2BGR6mS1sqF5c4ZkZnICaAusBcJt55WLR6RUyvL9rWksEZHqlJhIdGYmm4FQ4EegC2bJCaBwLh4RqRQKdkREqlN+jp0OmFvTrwIOY25VX+ugnYhUnMsFO/PmzSM8PJyAgAAiIiJIvMC/bpYsWcKVV15J7dq1CQ4O5o477uDo0aPV1FsRkTIqkGMnBNgE3AicBm4G3nHQTkQqxqWCneXLlxMTE8OMGTNISkoiMjKSfv36kZKS4rD95s2bGT16NOPHj+f7779nxYoV7NixgzvvvLOaey4iUkpFcvHUAz4BRgNWYBzwXIMGGN27O6+PIm7GpYKdl19+mfHjx3PnnXfSrl075syZQ1hYGPPnz3fY/uuvv6ZVq1ZMnjyZ8PBwunfvzj333MPOnTtLfI3s7GyOHz9e6CEiUm0c5OLxBd4FHslv8lhWFhOnTMFqtSoXj0glcJlg59y5c+zatYvo6OhCx6Ojo9m6davDa7p168bBgwdZu3YthmHwv//9j5UrV3LTTTeV+DqzZs2iQYMG9kdYWFil/h4iIhfkIBePBZgVFsZr48djsViYN28ew667jrPKxSNSYS4T7GRmZmK1WgkKCip0PCgoiIyMDIfXdOvWjSVLljB8+HD8/Pxo3rw5DRs25N///neJrzN9+nSysrLsj9TU1Er9PURESmXwYDhwAOLjYelS82dyMpMWLGD58uX4+fiwevt2+qal8UfB69LSYMgQBTwiZeAywY6NxZZUK59hGMWO2fzwww9MnjyZJ554gl27dvHFF1+QnJzMhAkTSry/v78/9evXL/QQEXEKb2+IioIRI8yf+Xl1hg4ezPqGDamPuYA5Ejhou8aWGi0mRlNaIqXkMsFOkyZN8Pb2LjaKc/jw4WKjPTazZs3iuuuu48EHH6RDhw707duXefPmsWjRItK1bVNEaqrERKIyM0kEgoHvga75PwHl4hEpI5cJdvz8/IiIiCAuLq7Q8bi4OLp16+bwmtOnT+PlVfhX8M7/l5GHJ4YWkZqsQC6ebZhZlg9ijvB87aCdiJyfywQ7ALGxsSxYsIBFixaxb98+pk6dSkpKin1aavr06YwePdrefsCAAaxevZr58+fz22+/sWXLFiZPnsw111xDSEiIs34NEZGKKZBjpyWwGTPL8h9AL8wq6kXbiUjJfJzdgYKGDx/O0aNHmTlzJunp6bRv3561a9fSsmVLANLT0wvl3Bk7diwnTpxg7ty5PPDAAzRs2JAbbriBF154wVm/gohIxdly8aSlgWEQCPwHGAxsAG4ClgYGMiQy0qndFKkpVAhUhUBFxBWtXm3uugL7ouRzwO3Ah4CXlxdvvPEGd911V4kV1kXcWZUWAl2/fj15eXnl7pyIiJSCg1w8fsDS0FDu7tOHvLw87r77bl64/XYz945y8YiUqMwjO97e3qSnp9OsWbOq6lO10siOiLg0B6M2hpcXM2bMYNasWQA8CLyAmZgQsGdmZuVKM2gScUNl+f4uc7Dj5eVFRkaGgh0REWeyWvlXYCAPZmUBMB54E7BPXlks5rqf5GRNaYlbqtJpLBERcQGJiUzLymIB5l/kC4HhQLbtvHLxiNiVK9iZO3cu69evJzMzs7L7IyIipZGfY2c8sAJzPc8q4O/AKQftRDxZubaez58/n2effRaLxUKLFi3o1KkTERERdOrUiU6dOhGs3A8iIlWrwN+zg4HPgUGYW9T75j9vUKSdiKcq95qd3NxckpKS2L17t/2RmpqKxWIhKCiIQ4cOVVWfK5XW7IhIjWS1mruu8nPxAGwHbgT+BCKA9SEhBKakaM2OuKWyfH+XeWTHVpQzJCSEkJAQbrrpJvu5Y8eOsXPnTvbs2VPW24qISFl4e8Orr5q5eCwWMAyuBeKBPsAuIMrHh/9kZpZYX1DEU2g3lkZ2RKQmW70apkyBg/a66PzQvDm9s7NJ/+MPLrvsMr788ktCg4OVeFDcSpWO7Kxbt44GDRqUu3MiIlKJBg+GgQMLBTKXR0ayKTmZXr168fPPP3N9RARfenkRnpHx13WhoebIkPLwiAdQuQiN7IiIm/r999/p1aUL+zMyaAF8CbSxnVTiQanhlGdHRERoGRrKJouFdkAacD3wne2k7d+5MTHmYmcRN6ZgR0TEXSUmEpKezkbgKuAwEAXstJ1X4kHxEAp2RETcVX5CwabAV8C1wDGgF7DVQTsRd6VgR0TEXRVIKNgIiAN6AMcxEw9udtBOxB0p2BERcVeRkeauq/zFyPWAtZgjOycxExBuatrUbCfixio92OnduzetW7eu7NuKiEhZ2RIPgj3gqQ18ipl48BTQ7/hxErRmR9xcpQc7t9xyC2PGjKns24qISHkMHmxuL2/Rwn6oFvBxixb0veoqTmdn079/f7766itzV1ZCAixbZv7ULi1xE8qzozw7IuIJrNZiGZTP5uQwePBg1q1bR4CfH5/Wr0/vzMy/rlHiQXFhZfn+rpRgx3YLW92smkTBjoh4suzsbG697jo+37WLAOAjzMXLgBIPikurtqSCCxcupH379gQEBBAQEED79u1ZsGBBRW4pIiLVyN/Hh1UZGQwAzgIDgXW2k0o8KG6i3MHO448/zpQpUxgwYAArVqxgxYoVDBgwgKlTp/LYY49VZh9FRKSqJCbin5bGSmAQkJ3/8zPbeSUeFDdQ5kKgNvPnz+ftt99mxIgR9mM333wzHTp0YNKkSTz77LOV0kEREalC+QkF/YAPgduA1cBgYBUwoEg7kZqo3CM7VquVzp07FzseERFBbm5uhTolIiLVpEBCQV/gA2AokAPcipmXp2g7kZqm3MHOqFGjmD9/frHjb731Fv/4xz8q1CkREakmRRIP+gJLgSGYAc9gYIMSD0oNV+5pLDAXKG/YsIEuXboA8PXXX5Oamsro0aOJjY21t3v55Zcr1ksREakatsSDQ4aYAY9h4IMZ8ORi7s4amJXFZwkJ9OrVy6ldFSmvcm8979mzZ+lewGIxk1W5KG09FxEBVq+GKVPg4EH7oXOhoQxp3pxPd+6kVq1arFu3jh7duxfL14O3txM7Lp6q2vPs1GQKdkRE8jlIPJidm8stt9zCunXrqBMQwBf16tH9yJG/rlHiQXESBTtloGBHROT8zp49y81duhD3zTfUBTYAXW0nlXhQnKQs398VWrNz9uxZvv32Ww4fPkxeXl6hczfffHNFbi0iIi4iwNeXj44cYQDwFWa19DjgGjDz8FgsZuLBgQM1pSUuqdzBzhdffMHo0aPJLFhHJZ/FYsGqbJsiIu4hMZHahw7xCdAf2AREA18CEVA48WBUlPP6KVKCcm89nzhxIkOHDiU9PZ28vLxCDwU6IiJuJD+hYB3gc+A6IAvoA+xx0E7E1ZQ72Dl8+DCxsbEEBQVVZn9ERMTVFEgoWBcz0WAX4A/MgOcHB+1EXEm5g50hQ4aQkJBQiV0RERGXVCTxYH3gC8wprEygN7Dftg1dxAWVezfW6dOnGTp0KE2bNuWKK67A19e30PnJkydXSgermnZjiYiUwurVZuJBsFdDPwpEAXuBlk2bkrhrF2FhYU7qoHiaatl6vmDBAiZMmECtWrUIDAzEYtt+iLlA+bfffivPbaudgh0RkVJykHgwIySEyLw8fs3I4NJLL2XTpk00b9pUiQelylVLsNO8eXMmT57MI488gpdXuWfDnE7BjohIGThIPJiSlkZkZCQpKSm0v+giEnJyCCy4WFmJB6UKVEuw07hxY3bs2MHFF19crk66CgU7IiIV9+uvvxJ59dVk/PknnTG3pdv/RlXiQakCZfn+LveQzJgxY1i+fHl5LxcRETdySXg4//H3JxDYCfwdOG07afs3dUyMOTIkUs3KnVTQarXy4osvsn79ejp06FBsgbIqnYuIeJDERP72v/+xAegJJAK3AJ8A/qDEg+JU5Q52vvvuOzp27AjA3r17C50ruFhZREQ8QP4anU7AOsz8OxuA4cAKwLdIO5HqVO5gJz4+vjL7ISIiNVmBhILdMEd0bgI+BsYC7wHeRdqJVJcKbaNKTExk1KhRdOvWjbS0NADef/99Nm/eXO57zps3j/DwcAICAoiIiCAxMfG87bOzs5kxYwYtW7bE39+fiy++mEWLFpX79UVEpByKJB7sBazE/Bf1UmASYISGKvGgOEW5g51Vq1bRt29fatWqxe7du8nOzgbgxIkTPP/88+W65/Lly4mJiWHGjBkkJSURGRlJv379SElJKfGaYcOG8eWXX7Jw4UJ++uknli1bRtu2bcv1+iIiUk7e3ub2crAHPH8H/g+wAPOBJ7p0Ub4dcQ6jnK666ipj8eLFhmEYRt26dY39+/cbhmEYSUlJRlBQULnuec011xgTJkwodKxt27bGI4884rD9unXrjAYNGhhHjx4t1+sZhmFkZWUZgJGVlVXue4iISL5VqwwjNNQwzCXJhgHG/EaNDMAAjFdeecXZPRQ3UZbv73KP7Pz0009cf/31xY7Xr1+fP//8s8z3O3fuHLt27SI6OrrQ8ejoaLZu3erwmk8++YTOnTvz4osv0qJFCy677DKmTZvGmTNnSnyd7Oxsjh8/XughIiKVZPBgOHAA4uNh6VKIj2fCkSM8++yzAEydOpX33nvP3IKekADLlpk/tSVdqlC5FygHBwfz66+/0qpVq0LHN2/eTOvWrct8v8zMTKxWa7Eq6kFBQWRkZDi85rfffmPz5s0EBASwZs0aMjMzue+++zh27FiJ63ZmzZrF008/Xeb+iYhIKXl7F9te/uijj3L06FFeeeUVxt1xBw1jY7n56NG/GijLslShco/s3HPPPUyZMoXt27djsVg4dOgQS5YsYdq0adx3333l7lDRbeuGYZS4lT0vLw+LxcKSJUu45ppr6N+/Py+//DLvvvtuiaM706dPJysry/5ITU0td19FRKR0LBYL//rXvxgTFYU1L49hR4+SULBBWppZaHT1aif1UNxZuUd2HnroIbKysujZsydnz57l+uuvx9/fn2nTpjFx4sQy369JkyZ4e3sXG8U5fPhwsdEem+DgYFq0aEGDBg3sx9q1a4dhGBw8eJBLL7202DX+/v74+/uXuX8iIlIxXobBgl9+4U/MLek3A/FABJireywWM8vywIFayCyVqkJbz5977jkyMzP573//y9dff82RI0d45plnynUvPz8/IiIiiIuLK3Q8Li6Obt26Obzmuuuu49ChQ5w8edJ+7Oeff8bLy4vQ0NBy9UNERKpIYiI+aWl8gJll+QRwI/CT7XzBLMsilajcwU5KSgqGYVC7dm06d+7MNddcQ926de3nyiM2NpYFCxawaNEi9u3bx9SpU0lJSWHChAmAOQU1evRoe/uRI0cSGBjIHXfcwQ8//MCmTZt48MEHGTduHLVq1SrvryYiIlUhP3tyAPAR5ohOJma25VQH7UQqS7mnscLDw0lPT6dZs2aFjh89epTw8HCs5VhZP3z4cI4ePcrMmTNJT0+nffv2rF27lpYtWwKQnp5eKJCqW7cucXFxTJo0ic6dOxMYGMiwYcPsq/5FRMSFFMieXB+zrEQk5shONGY9rSZF2olUBoth2MrRlo2Xlxf/+9//aNq0aaHjv//+O5dffjmnTp2qlA5WtbKUiBcRkQqwWqFVK3Mxcv5XTypwXf7Pa4EvW7Sgzu+/a82OXFBZvr/LPLITGxsLmCvrH3/8cWrXrm0/Z7Va2b59O1dddVVZbysiIu7OlmV5yBBzMbJhEIZZMPQ6YDswLCiIj/Ly8FWwI5WozMFOUlISYG4J/+677/Dz87Of8/Pz48orr2TatGmV10MREXEfgwfDypUwZQocPAhAW+CzZs3o9eefrN29m7vvvptFixaVmHZEpKzKPY11xx138Oqrr9b4qR9NY4mIOIHVau66Sk831+hERvLZunUMGjQIq9XK9OnTef6ZZ4q10fSW2JTl+7vcwY67ULAjIuI6Fi1axPjx4wF4rWFDJhUsP6Qsy1JAWb6/K5RnR0REpDKNGzeOZ0eMAGDKn3/yYcGTyrIs5aRgR0REXIfVyqObNnE/Zpn024GvbOdsExExMSocKmWiYEdERFxHYiKWtDReBYYA54BBwB7beWVZlnJQsCMiIq4jP3uyN/A+0AOzrEQ/INlBO5HSqPRgp3fv3rRu3bqybysiIp6gQPZkW1mJK4AMoC9wxEE7kQspd7mIktxyyy1kZmZW9m1FRMQTREaau67ysyw3BL4AugK/AAOAr1q0oHZkpDN7KTVMpQc7999/f2XfUkREPIWDLMshwHqgG2aW5X+0aMFKzKkukdLQmh0REXEttizLLVrYD7UFPm7WDH9fXz7673+ZOnUqHp4mTsqgTEkFbXWxSuPll18uV4eqm5IKioi4KAdZlj9ctYrhw4cDMHv2bGKnTFGWZQ9VZRmUe/bsWap2FouFr7766sINXYCCHRGRmuVf//oXDz74IAAfNm7M0GPH/jqpLMseQ+UiykDBjohIzWIYBpP69+f1L77AH/gP0N120lY8dOVKBTxuTuUiRETEbVny8nj1u++4GcgGBgI/2U4qy7I4UKFgJzExkVGjRtG1a1fS0tIAeP/999m8eXOldE5ERKSYxES809JYBlwDHMNMOvg/23llWZYiyh3srFq1ir59+1KrVi2SkpLIzs4G4MSJEzz//POV1kEREZFC8rMn1wY+BVpjZlceAJxy0E6k3MHOs88+yxtvvMHbb7+Nr6+v/Xi3bt3YvXt3pXRORESkmALZk5sB64BAYAcwArA6aCeerdzBzk8//cT1119f7Hj9+vX5888/K9InERGRktmyLOcvRr4M+ATwxxzpmQIYoaFmOxEqEOwEBwfz66+/Fju+efNm1cYSEZGqY8uyDPaApxuwBLAArwOv9emjfDtiV+5g55577mHKlCls374di8XCoUOHWLJkCdOmTeO+++6rzD6KiIgU5iDL8q3ACw0aADD13Xf59NNPndQ5cTUVyrMzY8YMXnnlFc6ePQuAv78/06ZN45lnnqm0DlY15dkREanBimRZNrp35+5772XBggXUqVOHxMREOnbs6OxeShWo8qSCOTk5REdH8+abbxIaGsoPP/xAXl4el19+OXXr1i13x51BwY6IiHvJycmhf//+/Oc//6FFixZs37qVFr/9ppISbqYs39/lqnru6+vL3r17sVgs1K5dm86dO5eroyIiIpXN19eXFStW0K1bN/bt28eASy5hU04O9n+Kq6SExyn3mp3Ro0ezcOHCyuyLiIhIpWjYsCGfT5pEUyApJ4eRFNiSnpYGQ4bA6tXO66BUq3KN7ACcO3eOBQsWEBcXR+fOnalTp06h8zWl6rmIiLghq5Xw55/nY6An5pb0B4GXwcywbLGYJSUGDtSUlgcod7Czd+9eOnXqBMDPP/9c6JzFVohNRETEGRIT4eBBugKLgduAV4BLgXuhcEmJqCjn9VOqRbmDnfj4+Mrsh4iISOUpUCpiOPAr8BgwCQgHbnTQTtyXqp6LiIj7KVIq4lFgDOa6nWHAdyW0E/ekYEdERNxPkZISFuAtoAdwAvg78L+QEJWU8BAKdkRExP04KCnhB6zGXLeTAgyuV4/s3FwndVCqk4IdERFxTw5KSjQGPm3enIZ16rD1p5+4++67qUAhAakhFOyIiIj7GjwYDhyA+HhYuhTi42lz8CAfrlmDt7c37733Hi+99JJZdiIhAZYtM39arRe4sdQkFaqN5Q5ULkJExDPNnTuXSZMmYbFY+KhxY24+evSvk8qy7PLK8v2tkR0REfFI999/PxOiozEMg38cPfrXDi1QlmU3o2BHREQ8kiUvj9e+/54bgJPAAOCw7aRt0iMmRlNabkDBjoiIeKbERHzT0lgBXAL8DgwGsm3nC2ZZlhpNwY6IiHim/OzJjTFrZzUAtgATAMNBO6m5FOyIiIhnKpA9uS3wIeaX4rvAv0poJzWTgh0REfFMRbIsRwNz8k89jDnaQ1iYsiy7AQU7IiLimRxkWZ4I3IM5jTUS+D421mwnNZqCHRER8VxFsixbgH8DUf7+nAQGzp3LsWPHnNlDqQQuF+zMmzeP8PBwAgICiIiIILGUq+C3bNmCj48PV111VdV2UERE3EuRLMu+8fGs+P13WrZsyf79+7ntttvIVQ2tGs2lgp3ly5cTExPDjBkzSEpKIjIykn79+pGSknLe67Kyshg9ejS9evWqpp6KiIhb8faGqCgYMQKiomgSFMTHH39M7dq1iYuL45FHHlFJiRrMpcpFXHvttXTq1In58+fbj7Vr145BgwYxa9asEq+77bbbuPTSS/H29uajjz5iz549pX5NlYsQEZGSrFixgmHDhgHwXqNG3P7HH3+dVEkJp6qR5SLOnTvHrl27iI6OLnQ8OjqarVu3lnjdO++8w/79+3nyySdL9TrZ2dkcP3680ENERMSRoUOHMuPWWwG4648/2FnwpEpK1BguE+xkZmZitVoJCgoqdDwoKIiMjAyH1/zyyy888sgjLFmyBB8fn1K9zqxZs2jQoIH9ERYWVuG+i4iIm7Jamfn11wzAzKw8CLB/I6mkRI3hMsGOjSV/+5+NYRjFjgFYrVZGjhzJ008/zWWXXVbq+0+fPp2srCz7IzU1tcJ9FhERN5WYiFdaGv8HtAPSgFtRSYmapnTDIdWgSZMmeHt7FxvFOXz4cLHRHoATJ06wc+dOkpKSmDhxIgB5eXkYhoGPjw8bNmzghhtuKHadv78//v7+VfNLiIiIe8kvFVEf+Bi4BtgK3A+8jblVvWA7cU0uM7Lj5+dHREQEcXFxhY7HxcXRrVu3Yu3r16/Pd999x549e+yPCRMm0KZNG/bs2cO1115bXV0XERF3VaBUxKXAMswvzoXAvBLaietxmZEdgNjYWG6//XY6d+5M165deeutt0hJSWHChAmAOQWVlpbGe++9h5eXF+3bty90fbNmzQgICCh2XEREpFxsJSXS0sAwuBH4J/AQEAP8DYhSSQmX51LBzvDhwzl69CgzZ84kPT2d9u3bs3btWlq2bAlAenr6BXPuiIiIVBpbSYkhQ8ySEobBNGAPsBQYCux89FFaqqSES3OpPDvOoDw7IiJyQatXw5QpcPAgAGeASF9fduXk0KlTJzZv3kytWrWc20cPUyPz7IiIiLisIiUlasXHs/qnn2jSpAm7d+/m3nvvxcPHDlyaS01jiYiIuCxbSYl8F2GWOerTpw+LFy/m6quv5v7773da96RkGtkREREppxtuuIEXXngBgJiYGDZv3Kj6WS5Ia3a0ZkdERCrAMAxGjBjB8uXLae7lxa68PEJsJ1U/q8pozY6IiEg1sVgsLLz5ZtoDGXl5DAHO2U6qfpZLULAjIiJSEVYrdR5+mDVAQ2AbZg4eQPWzXISCHRERkYpITISDB7kEWIJZQmI+8I7tvOpnOZ2CHRERkYooUBerP/B0/p/vBXaW0E6ql4IdERGRiihSF2sGcDNmZfTBwOES2kn1UbAjIiJSEbb6WRazBroX8B5wGZAKDAdyQ0NVP8uJFOyIiIhUhK1+FtgDngbAGqAukABM79TJbCdOoWBHRESkogYPhpUroUUL+6HLgXcDAwH41yefsFrbz51GwY6IiEhlKFI/i/h4bv3f/3jggQcAGDt2LD///LNz++ihlEFZGZRFRKQK5eTk0KtXLxITE2nfvj1fb9lCnd27zd1ZwcHmWh5NcZVZWb6/VQhURESkCvn6+rJ8+XI6duzI3r17mRAczHunT2OxNVBJiSqnaSwREZEqFhwczIf334838H+nT/NGwZMqKVHlFOyIiIhUNauV6996i3/mP50C/Nd2TiUlqpyCHRERkaqWX1LiAcxEgznAECDTdl4lJaqUgh0REZGqll8qwgIsAi7FTDj4D8DqoJ1ULgU7IiIiVa1AqYgGwCqgFrABmFlCO6k8CnZERESqWpGSElcAb+WfmgmsBQgLU0mJKqJgR0REpKo5KCkxCrMyuu3PBx59VPl2qoiCHRERkergoKTEK8DVfn78AQxZsIDs7Gyndc+dKdgRERGpLkVKSvjHx7Pyxx9p3Lgxu3btYtq0ac7uoVtSuQiVixARESf7/PPP+fvf/w7Ahx9+yNChQ53cI9dXlu9vjeyIiIg42U033cTDDz8MwPjx4/nlxx8hIQGWLTN/KtlghWhkRyM7IiLiAnJzc7nhhhtITEzkKl9ftuXkEGA7qfpZxWhkR0REpIbx8fFh2e230xTYk5NDTMGTqp9VIQp2REREXIHVSouZM/k/zEzLbwJLbedUP6tCFOyIiIi4gvz6WdHAY/mH7gZ+tJ1X/axyU7AjIiLiCgrUxXoS6AmcAoYCp0toJ6WjYEdERMQVFKiL5Y05hRUE7AUmltBOSkfBjoiIiCsoUj+rObAM84v6HeBdUP2sclKwIyIi4goc1M/qyV9V0e8D9sbGqn5WOSjYERERcRUO6mdNB/oGBHAGGPrmm5w6dcpp3aupFOyIiIi4kiL1s7zi43k/OZmQkBB+/PFHJk6ceMFbSGHKoKwMyiIiUgNs3LiRG264gby8PN5//31GjRrl7C45lTIoi4iIuJkePXrw5JNPAjBhwgR+3rdP9bNKSSM7GtkREZEawmq10qdPH+Lj4z2+fpZGdkRERNyQt7c3/zdypL1+1rSCJ1U/q0QKdkRERGoKq5WQp5/mvfynrwP20Eb1s0qkYEdERKSmyK+fdSPwYP6h8cAB23nVz3JIwY6IiEhNUaAu1nPAtcCfwAggp4R24oLBzrx58wgPDycgIICIiAgSzxOdrl69mj59+tC0aVPq169P165dWb9+fTX2VkREpBoVqIvlC3wANAC+Bh4voZ24WLCzfPlyYmJimDFjBklJSURGRtKvXz9SUlIctt+0aRN9+vRh7dq17Nq1i549ezJgwACSkpKqueciIiLVoEj9rFbAwvxTLwDrQfWzHHCprefXXnstnTp1Yv78+fZj7dq1Y9CgQcyaNatU9/jb3/7G8OHDeeKJJxyez87OJjs72/78+PHjhIWFaeu5iIjUDKtXm7uuwL4o+X5gHpi7tN5+m5A773RW76pNjdx6fu7cOXbt2kV0dHSh49HR0WzdurVU98jLy+PEiRM0bty4xDazZs2iQYMG9kdYWFiF+i0iIlKtHNTPmg1c6evLEWDU0qVYtRurEJcJdjIzM7FarQQFBRU6HhQUREZGRqnuMXv2bE6dOsWwYcNKbDN9+nSysrLsj9TU1Ar1W0REpNoVqZ8VEB/P8j17qFOnDvHx8bzwwgvO7qFL8XF2B4qy5M9D2hiGUeyYI8uWLeOpp57i448/plmzZiW28/f3x9/fv8L9FBERcSpvb4iKsj9tA8ydO5c77riDJ554ghtuuIEuXbo4rXuuxGVGdpo0aYK3t3exUZzDhw8XG+0pavny5YwfP54PP/yQ3r17V2U3RUREXNaYMWO47bbbsFqtjBw5kqxjx1Q/CxcKdvz8/IiIiCAuLq7Q8bi4OLp161bidcuWLWPs2LEsXbqUm266qaq7KSIi4rIsFgtvvPEGrVq1Ijk5mfvDwqBnTxg50vzZqpVHlpNwmWAHIDY2lgULFrBo0SL27dvH1KlTSUlJYcKECYC53mb06NH29suWLWP06NHMnj2bLl26kJGRQUZGBllZWc76FURERJyqQYMGLL37bryBJadP837Bkx5aP8ulgp3hw4czZ84cZs6cyVVXXcWmTZtYu3YtLVu2BCA9Pb1Qzp0333yT3Nxc7r//foKDg+2PKVOmOOtXEBERcS6rla7z5vFU/tP7gF9t5zy0fpZL5dlxhrLs0xcREXF5CQnQsydW4AZgE3A1sBnwK9guPr7QAueapkbm2REREZFKkF8Xyxv4P6ARsAMolmrXg+pnKdgRERFxJwXqYoUBC/L//CLwZQnt3J2CHREREXdSpH7WYOBuwABuB46Ax9XPUrAjIiLiTry94dVXzT/nBzyvAO2AdGA8YLzyitnOQyjYERERcTdF6mfVBpZhLlD+FHjdg9brgHZjaTeWiIi4L6sVEhPNxcjBwby2Zw9Tpk7F39+fHTt2cMUVVzi7h+VWlu9vBTsKdkRExEMYhsGAAQP4/PPPad++PTt27CAgIMDZ3SoXbT0XERGRYiwWC4sWLaJZs2bs3buX6dOnO7tL1UIjOxrZERERD7N27Vp7Pcn1a9cSXauWfaqLyMgasXhZIzsiIiJSov79+3PfffcBMPbvfyfTzYuFKtgRERHxQC91705bID0vz56HB3DLYqEKdkRERDyN1Urthx5iKeALrAEW2c65YbFQBTsiIiKeJjERDh6kI/Bs/qEpwC+284YBqalmOzegYEdERMTTFEgq+AAQBZwCRgE5JbSryRTsiIiIeJoCRUC9gfeAhsB/gWdKaFeTKdgRERHxNEWKhYYBb+Sfeg7YAm5VLFTBjoiIiKdxUCx0OGZV9DzM6azjzz9fI/LtlIaCHREREU9UpFgowFyglbc3B4BJcXHO6lmlUwZlZVAWERFPVqRY6GaLhR433EBeXh7Lly9n2LBhzu6hQyoEWgYKdkRERAp77LHHeO6552jUqBF79+4lJCTE2V0qRuUiREREpNyefPJJOnXqxB9//MG4ceOo6eMiCnZERESkEF9fX95//30CAgJYv349b8ybBwkJsGyZ+bOGZVZWsCMiIiLFXH755fzzn/8E4IFJk/i5BhcLVbAjIiIiDk1q0YJewBnD4HYg13aihhULVbAjIiIixVmteE2dyjtAA8zsyrNs52pYsVAFOyIiIlJcfrHQMOD1/EMzgZ228zWoWKiCHRERESmuQBHQkcBQzGms24EzJbRzVQp2REREpLgCRUAtwHwgGPgReKSEdq5KwY6IiIgUV6RYaCCwKP/Ua8B/oMYUC1WwIyIiIsU5KBZ6I3Bv/umxwB/PPlsjioUq2BERERHHHBQLfQm41MeHNGDi+vVO61pZqDaWamOJiIicX5FioV/7+nLd9deTl5fHBx98wPDhw6u9SyoEWgYKdkRERMruiSee4JlnnqFx48Z8//33NG/evFpfX4VARUREpEo99thjdOzYkWPHjnHXXXe5dLFQBTsiIiJSZn5+fixevBhfX18+++wzFr/zjssWC1WwIyIiIuVyxRVXMHPmTACm3HknqS5aLFTBjoiIiJTbtIsv5lrguGEwHrBPZrlQsVAFOyIiIlI+Vis+sbEsBgKAOOBN2zkXKhaqYEdERETKJ79YaBvgn/mHpgG/2c67SLFQBTsiIiJSPgWKgE4CegCnMLMr55XQzhkU7IiIiEj5FCgC6gW8A9QBEoFXS2jnDAp2REREpHyKFAsNB2bnn3oUs0I6TZuai5WduB3d5YKdefPmER4eTkBAABERESReYJ5v48aNREREEBAQQOvWrXnjjTeqqaciIiIezkGx0LuBaOAsMAbIPXIERo1y6nZ0lwp2li9fTkxMDDNmzCApKYnIyEj69etHSkqKw/bJycn079+fyMhIkpKSePTRR5k8eTKrVq2q5p6LiIh4qCLFQi3AQqAB8F/MwqF2TtqO7lK1sa699lo6derE/Pnz7cfatWvHoEGDmDVrVrH2Dz/8MJ988gn79u2zH5swYQLffPMN27ZtK9VrqjaWiIhIJbAVC01Lg6lTee/IEcYAvsBOoIOtncViTn0lJ5sjQ+VUI2tjnTt3jl27dhEdHV3oeHR0NFu3bnV4zbZt24q179u3Lzt37iQnJ8fhNdnZ2Rw/frzQQ0RERCrI2xuioswRniNHuB24GcjBnM46Z2vnhO3oLhPsZGZmYrVaCQoKKnQ8KCiIjIwMh9dkZGQ4bJ+bm0tmZqbDa2bNmkWDBg3sj7CwsMr5BURERMS+zdyCmWAwEGgMFBtaqMbt6C4T7NhY8hc42RiGUezYhdo7Om4zffp0srKy7I/U1NQK9lhERETsCmwzbw58jZlZucl52lU1n2p7pQto0qQJ3t7exUZxDh8+XGz0xqZ58+YO2/v4+BAYGOjwGn9/f/z9/Sun0yIiIlKYbTt6WhoYBpcUPW9bsxMZWW1dcpmRHT8/PyIiIoiLiyt0PC4ujm7dujm8pmvXrsXab9iwgc6dO+Pr61tlfRUREZESONiObmd7PmdOhRYnl5XLBDsAsbGxLFiwgEWLFrFv3z6mTp1KSkoKEyZMAMwpqNGjR9vbT5gwgd9//53Y2Fj27dvHokWLWLhwIdOmTXPWryAiIiJFtqPbhYaaxwcPrtbuuMw0FsDw4cM5evQoM2fOJD09nfbt27N27VpatmwJQHp6eqGcO+Hh4axdu5apU6fy+uuvExISwmuvvcatt97qrF9BREREwAxoBg40d12lp5trdCIjq3VEx8al8uw4g/LsiIiI1Dw1Ms+OiIiISFVQsCMiIiJuTcGOiIiIuDUFOyIiIuLWFOyIiIiIW1OwIyIiIm5NwY6IiIi4NQU7IiIi4tYU7IiIiIhbc6lyEc5gSyB9/PhxJ/dERERESsv2vV2aQhAeH+ycOHECgLCwMCf3RERERMrqxIkTNGjQ4LxtPL42Vl5eHocOHaJevXpYipair6Djx48TFhZGamqq6m65OL1XNYver5pD71XNUdPeK8MwOHHiBCEhIXh5nX9VjseP7Hh5eREaGlqlr1G/fv0a8T+O6L2qafR+1Rx6r2qOmvReXWhEx0YLlEVERMStKdgRERERt6Zgpwr5+/vz5JNP4u/v7+yuyAXovapZ9H7VHHqvag53fq88foGyiIiIuDeN7IiIiIhbU7AjIiIibk3BjoiIiLg1BTsiIiLi1hTsVNC8efMIDw8nICCAiIgIEhMTz9t+48aNREREEBAQQOvWrXnjjTeqqadSlvcqISEBi8VS7PHjjz9WY48906ZNmxgwYAAhISFYLBY++uijC16jz5XzlPX90mfLOWbNmsXVV19NvXr1aNasGYMGDeKnn3664HXu8tlSsFMBy5cvJyYmhhkzZpCUlERkZCT9+vUjJSXFYfvk5GT69+9PZGQkSUlJPProo0yePJlVq1ZVc889T1nfK5uffvqJ9PR0++PSSy+tph57rlOnTnHllVcyd+7cUrXX58q5yvp+2eizVb02btzI/fffz9dff01cXBy5ublER0dz6tSpEq9xq8+WIeV2zTXXGBMmTCh0rG3btsYjjzzisP1DDz1ktG3bttCxe+65x+jSpUuV9VFMZX2v4uPjDcD4448/qqF3UhLAWLNmzXnb6HPlOkrzfumz5RoOHz5sAMbGjRtLbONOny2N7JTTuXPn2LVrF9HR0YWOR0dHs3XrVofXbNu2rVj7vn37snPnTnJycqqsr56uPO+VTceOHQkODqZXr17Ex8dXZTelnPS5qpn02XKurKwsABo3blxiG3f6bCnYKafMzEysVitBQUGFjgcFBZGRkeHwmoyMDIftc3NzyczMrLK+erryvFfBwcG89dZbrFq1itWrV9OmTRt69erFpk2bqqPLUgb6XNUs+mw5n2EYxMbG0r17d9q3b19iO3f6bHl81fOKslgshZ4bhlHs2IXaOzoula8s71WbNm1o06aN/XnXrl1JTU3lX//6F9dff32V9lPKTp+rmkOfLeebOHEi3377LZs3b75gW3f5bGlkp5yaNGmCt7d3sZGBw4cPF4uEbZo3b+6wvY+PD4GBgVXWV09XnvfKkS5duvDLL79UdvekgvS5qvn02ao+kyZN4pNPPiE+Pp7Q0NDztnWnz5aCnXLy8/MjIiKCuLi4Qsfj4uLo1q2bw2u6du1arP2GDRvo3Lkzvr6+VdZXT1ee98qRpKQkgoODK7t7UkH6XNV8+mxVPcMwmDhxIqtXr+arr74iPDz8gte41WfLiYuja7wPPvjA8PX1NRYuXGj88MMPRkxMjFGnTh3jwIEDhmEYxiOPPGLcfvvt9va//fabUbt2bWPq1KnGDz/8YCxcuNDw9fU1Vq5c6axfwWOU9b165ZVXjDVr1hg///yzsXfvXuORRx4xAGPVqlXO+hU8xokTJ4ykpCQjKSnJAIyXX37ZSEpKMn7//XfDMPS5cjVlfb/02XKOe++912jQoIGRkJBgpKen2x+nT5+2t3Hnz5aCnQp6/fXXjZYtWxp+fn5Gp06dCm3jGzNmjNGjR49C7RMSEoyOHTsafn5+RqtWrYz58+dXc489V1neqxdeeMG4+OKLjYCAAKNRo0ZG9+7djc8//9wJvfY8tq3JRR9jxowxDEOfK1dT1vdLny3ncPQeAcY777xjb+POny2LYeSvNhIRERFxQ1qzIyIiIm5NwY6IiIi4NQU7IiIi4tYU7IiIiIhbU7AjIiIibk3BjoiIiLg1BTsiIiLi1hTsiIiIiFtTsCMiIiJuTcGOiIiIuDUFOyJSJaKiooiJiSnxeVmvd4bS9MEV+iki5+fj7A6IiGdYvXo1vr6+zu5GmRTtc1RUFFdddRVz5swpsY2IuB4FOyJSonPnzuHn51cp92rcuHGl3Kc6labPzvy9cnNz8fHRX+MiF6JpLBGxi4qKYuLEicTGxtKkSRP69OnDF198Qffu3WnYsCGBgYH8/e9/Z//+/YWuO3XqFKNHj6Zu3boEBwcze/Zsh/e2TfeU5p6l7evEiRPt93nssccwDMPeJjs7m8mTJ9OsWTMCAgLo3r07O3bssJ9fuXIlV1xxBbVq1SIwMJDevXtz6tQph30eO3YsGzdu5NVXX8VisWCxWDhw4ECxaawLvWZUVBSTJ0/moYceonHjxjRv3pynnnrqgr/vgQMHsFgsrFy5kuuvvx5/f3/WrFlTpv9mIp5KwY6IFLJ48WJ8fHzYsmULb775JqdOnSI2NpYdO3bw5Zdf4uXlxS233EJeXp79mgcffJD4+HjWrFnDhg0bSEhIYNeuXSW+RmnuWZa+bt++nddee41XXnmFBQsW2M8/9NBDrFq1isWLF7N7924uueQS+vbty7Fjx0hPT2fEiBGMGzeOffv2kZCQwODBgwsFSwW9+uqrdO3albvuuov09HTS09MJCwsr1u58r1mw33Xq1GH79u28+OKLzJw5k7i4uPP+rnv27AHghRde4PHHH+f7778nOjq6TP+9RDyWISKSr0ePHsZVV1113jaHDx82AOO7774zDMMwTpw4Yfj5+RkffPCBvc3Ro0eNWrVqGVOmTCl074LPz3fPC7W3nW/Xrp2Rl5dnP/bwww8b7dq1MwzDME6ePGn4+voaS5YssZ8/d+6cERISYrz44ovGrl27DMA4cODAeV/jQr9DwWMXek1b++7duxe6x9VXX208/PDDJfbDMAzjqaeeMurUqWMkJyeft52IFKeRHREppHPnzoWe79+/n5EjR9K6dWvq169PeHg4ACkpKfbz586do2vXrvZrGjduTJs2bUp8jQvds7S6dOmCxWKxP+/atSu//PILVquV/fv3k5OTw3XXXWc/7+vryzXXXMO+ffu48sor6dWrF1dccQVDhw7l7bff5o8//ijT6zv6vc73mjYdOnQodF1wcDCHDx8+77337NnDzTffTKtWrSrURxFPpGBHRAqpU6dOoecDBgzg6NGjvP3222zfvp3t27cD5uJloMRpn/O50D0rg61fBYMh23GLxYK3tzdxcXGsW7eOyy+/nH//+9+0adOG5OTkKntNm6K7tywWywWn8L755huioqLK3TcRT6ZgR0RKdPToUfbt28djjz1Gr169aNeuXbHRj0suuQRfX1++/vpr+7E//viDn3/+udz3LK2Cr2l7fumll+Lt7c0ll1yCn58fmzdvtp/Pyclh586dtGvXDjCDjOuuu46nn36apKQk/Pz8zrvo18/PD6vVWuL50rxmeRw/fpwDBw7QsWPHct9DxJNpz6KIlKhRo0YEBgby1ltvERwcTEpKCo888kihNnXr1mX8+PE8+OCDBAYGEhQUxIwZM/DycvxvqdLcs7RSU1OJjY3lnnvuYffu3fz73/+27wSrU6cO9957Lw8++CCNGzfmoosu4sUXX+T06dOMHz+e7du38+WXXxIdHU2zZs3Yvn07R44cOW9Q0qpVK7Zv386BAweoW7dusW3nF3rN8vrmm2/w8vLiiiuuKPc9RDyZgh0RKZGXlxcffPABkydPpn379rRp04bXXnut2HTKSy+9xMmTJ7n55pupV68eDzzwAFlZWRW6Z2mMHj2aM2fOcM011+Dt7c2kSZO4++677ef/+c9/kpeXx+23386JEyfo3Lkz69evp1GjRtSvX59NmzYxZ84cjh8/TsuWLZk9ezb9+vUr8fWmTZvGmDFjuPzyyzlz5ozDKa/zvWZ5ffPNN7Rt25aAgIBy30PEk1mM8ky4i4g4maNsxiIijmjNjoiIiLg1BTsiIiLi1jSNJSIiIm5NIzsiIiLi1hTsiIiIiFtTsCMiIiJuTcGOiIiIuDUFOyIiIuLWFOyIiIiIW1OwIyIiIm5NwY6IiIi4NQU7IiIi4tYU7IiIiIhb+3/smBvOEmPjaQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(rnum, Tnum, 'ro', label='FVM')\n",
    "plt.plot(r_an, T_an, 'k', label='analytic')\n",
    "plt.ylabel('rel. temp. $T$')\n",
    "plt.xlabel('radial position $r$')\n",
    "plt.legend(); # semicolon avoids output of '<matplotlib.legend.Legend...>' descriptor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantitative testing of agreement between FV result and analytic solution\n",
    "\n",
    "The tolerance is manually tuned. This should help detect code changes that break good agreement with experiment and theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytic solution evaluated at cell positions\n",
    "T_an_rnum = (S_val/(4*k_val))*(R**2 - rnum**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normalized error\n",
    "norm_err = (Tnum-T_an_rnum)/T_an_rnum.max() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAG0CAYAAAAFP+K7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/lklEQVR4nO3de3RU9b3//9ckmckNMuWimeQLhYBdQAQUkqJQbtaaYKu22irlLAP12B6wh3ITISgtVFsI1IK2ClhMb7YVjgS66K8ViUIiSIgCAQNJsUIiWJLGgM6gSK6f3x80U8dcmBlmGyY8H2vNwtn7vT/7PTMM++WePZ+xGWOMAAAAYImIzm4AAACgKyNsAQAAWIiwBQAAYCHCFgAAgIUIWwAAABYibAEAAFiIsAUAAGChqM5uAFJzc7NOnTql7t27y2azdXY7AADAD8YYnT17VsnJyYqIaP/8FWHrMnDq1Cn17du3s9sAAABBOHnypPr06dPuesLWZaB79+6SLrxYCQkJndwNAADwh8fjUd++fb3H8fYQti4DLR8dJiQkELYAAAgzF7sEiAvkAQAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsxAzyV5CmZqPXK86o5ux5Xd09RqNSeioywkYNNdRQ027N5dgTNdRc6t/pz9plEbbWrFmjn/3sZ6qqqtK1116rJ554QuPGjWu3vrCwUPPmzdORI0eUnJysBQsWaMaMGT41eXl5+uEPf6hjx45p4MCB+ulPf6o777wzoP1u3rxZzzzzjPbv36/Tp0+rpKRE119/vc8YdXV1mj9/vp5//nl9/PHHuvnmm7VmzZoOf5CyM2w7XKUf/6VMVe7z3mVJzhgtuT1Vk4YmUUMNNdS0qrkce6KGmkup6Sw2Y4zpzAY2btyorKwsrVmzRl/60pf0zDPP6Nlnn1VZWZk+//nPt6qvqKjQ0KFD9b3vfU/Tp0/Xa6+9pu9///t6/vnn9c1vflOSVFRUpHHjxumxxx7TnXfeqS1btuhHP/qRdu/erRtuuMHv/T733HOqqKhQcnKyvve977UZth544AH95S9/0W9/+1v16tVLDz74oM6cOaP9+/crMjLSr+fA4/HI6XTK7XZb8tuI2w5X6YE/HNCnX+iWrL/23pGSRA011FDjrZk0NIl/O6jpUjVWBC5/j9+dHrZuuOEGjRw5UmvXrvUuGzJkiL7xjW9o+fLlreoXLlyorVu3qry83LtsxowZOnTokIqKiiRJkydPlsfj0YsvvuitmTRpknr06KHnn38+4P1WVlYqJSWlVdhyu9266qqr9Nxzz2ny5MmSpFOnTqlv377629/+pszMTL+eAyvDVlOz0dgVO3yS/ifZJCUmREuyqdpDDTXUUCO5nDEqfOgmTfjZTv7toKZL1LicMdq98Msh/0jR3+N3p36MWF9fr/379ys7O9tneUZGhvbs2dPmNkVFRcrIyPBZlpmZqdzcXDU0NMhut6uoqEhz585tVfPEE08Evd+27N+/Xw0NDT79JCcna+jQodqzZ0+7Yauurk51dXXe+x6Px+99Bur1ijPt/mMpSUZStaeu3fXUUEPNlVdT5T6v54oq+beDmi5TU+U+r9crzmj0wF4d1lqlU7+NWFtbq6amJiUmJvosT0xMVHV1dZvbVFdXt1nf2Nio2traDmtaxgxmv+314nA41KNHj4DGWb58uZxOp/fWt29fv/cZqJqz7f9jCQDteefMuc5uAQipzjweXhZTP9hsvqf1jDGtll2s/tPL/Rkz0P3662LjLFq0SG6323s7efLkJe+zPVd3j7FsbABdV7+ecZ3dAhBSnXk87NSw1bt3b0VGRrY6C1RTU9PqrFMLl8vVZn1UVJR69erVYU3LmMHst71e6uvr9f777wc0TnR0tBISEnxuVhmV0lNJzhi1F/1sklwJ0XIlUEMNNdRcqElyxihrdH/+7aCmy9QkOS9MA9FZOjVsORwOpaWlKT8/32d5fn6+xowZ0+Y2o0ePblW/fft2paeny263d1jTMmYw+21LWlqa7Ha7zzhVVVU6fPhwQONYKTLCpiW3p0pSq7+ILfeX3nGtlt5BDTXUUHPBkttT5YiK4N8OarpMzZLbUzt1vq1O/xhx3rx5evbZZ/XrX/9a5eXlmjt3rk6cOOGdN2vRokWaOnWqt37GjBl65513NG/ePJWXl+vXv/61cnNzNX/+fG/N7NmztX37dq1YsUJ///vftWLFCr388suaM2eO3/uVpDNnzujgwYMqKyuTJB09elQHDx70nhFzOp26//779eCDD+qVV15RSUmJ7r33Xg0bNkxf+cpXrHzaAjJpaJLW3jtSLqfvKVSXM8b7dVhqqKGGmk/WSPzbQU3XqulMnT71g3RhctGVK1eqqqpKQ4cO1erVqzV+/HhJ0ne+8x1VVlaqoKDAW19YWKi5c+d6JzVduHBhq0lNN23apMWLF+v48ePeSU3vuusuv/crSb/97W913333tep3yZIlWrp0qSTp/Pnzeuihh/SnP/3JZ1LTQC56t3qerRaX2yy+1FBDzeVfczn2RA01l8sM8mEzzxY+u7AFAABCx9/jd6d/jAgAANCVEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKXRdhas2aNUlJSFBMTo7S0NO3atavD+sLCQqWlpSkmJkYDBgzQunXrWtXk5eUpNTVV0dHRSk1N1ZYtWwLerzFGS5cuVXJysmJjYzVx4kQdOXLEp6a6ulpZWVlyuVyKj4/XyJEjtWnTpiCeBQAA0BV1etjauHGj5syZo0ceeUQlJSUaN26cbr31Vp04caLN+oqKCn31q1/VuHHjVFJSoocfflizZs1SXl6et6aoqEiTJ09WVlaWDh06pKysLN1zzz0qLi4OaL8rV67UqlWr9NRTT+mNN96Qy+XSLbfcorNnz3prsrKydPToUW3dulWlpaW66667NHnyZJWUlFjwbAEAgLBjOtmoUaPMjBkzfJYNHjzYZGdnt1m/YMECM3jwYJ9l06dPNzfeeKP3/j333GMmTZrkU5OZmWm+/e1v+73f5uZm43K5TE5Ojnf9+fPnjdPpNOvWrfMui4+PN7///e99xunZs6d59tln233Mn+Z2u40k43a7/d4GAAB0Ln+P3516Zqu+vl779+9XRkaGz/KMjAzt2bOnzW2Kiopa1WdmZmrfvn1qaGjosKZlTH/2W1FRoerqap+a6OhoTZgwwae3sWPHauPGjTpz5oyam5u1YcMG1dXVaeLEie0+7rq6Onk8Hp8bAADomjo1bNXW1qqpqUmJiYk+yxMTE1VdXd3mNtXV1W3WNzY2qra2tsOaljH92W/LnxfrbePGjWpsbFSvXr0UHR2t6dOna8uWLRo4cGC7j3v58uVyOp3eW9++fdutBQAA4a3Tr9mSJJvN5nPfGNNq2cXqP73cnzFDUbN48WK9//77evnll7Vv3z7NmzdPd999t0pLS9vtf9GiRXK73d7byZMn260FAADhLaozd967d29FRka2OotVU1PT6oxSC5fL1WZ9VFSUevXq1WFNy5j+7Nflckm6cIYrKSmpzZpjx47pqaee0uHDh3XttddKkq677jrt2rVLTz/9dJvfkpQufBwZHR3dzrMCAAC6kk49s+VwOJSWlqb8/Hyf5fn5+RozZkyb24wePbpV/fbt25Weni673d5hTcuY/uw3JSVFLpfLp6a+vl6FhYXemnPnzkmSIiJ8n8bIyEg1Nzdf/AkAAABdXyBX3W/bts00NTUFe9F+mzZs2GDsdrvJzc01ZWVlZs6cOSY+Pt5UVlYaY4zJzs42WVlZ3vrjx4+buLg4M3fuXFNWVmZyc3ON3W43mzZt8ta89tprJjIy0uTk5Jjy8nKTk5NjoqKizN69e/3erzHG5OTkGKfTaTZv3mxKS0vNlClTTFJSkvF4PMYYY+rr680111xjxo0bZ4qLi83bb79tHn/8cWOz2cxf//pXv58Dvo0IAED48ff4HVDYioiIMP/6178uqbG2PP3006Zfv37G4XCYkSNHmsLCQu+6adOmmQkTJvjUFxQUmBEjRhiHw2H69+9v1q5d22rMF154wQwaNMjY7XYzePBgk5eXF9B+jbkw/cOSJUuMy+Uy0dHRZvz48aa0tNSn5q233jJ33XWXufrqq01cXJwZPnx4q6kgLoawBQBA+PH3+G0z5t9Xl/shIiJC1dXVuvrqq6060XZF8ng8cjqdcrvdSkhI6Ox2AACAH/w9fl8W30YEAADoqgIOW0899ZReeukl75xWAAAAaF/AHyP26tVLp0+fls1m0//7f/9PI0eOVFpamkaOHKmRI0f6TJMA//AxIgAA4cff43fA82wdOXJEjY2NKikp0YEDB3TgwAE9++yzOnnypGw2mxITE3Xq1KlLah4AAKCrCChstcycnpycrOTkZH3ta1/zrjtz5oz27dungwcPhrRBAACAcBZQ2OroE8eePXsqIyOj1Y87AwAAXMkCukD+xRdflNPptKoXAACALiegM1uZmZlW9QEAANAlMc8WAACAhQhbAAAAFiJsAQAAWIiwBQAAYKGAJzVtcf78eb355puqqalRc3Ozz7o77rjjkhsDAADoCoIKW9u2bdPUqVPb/H1Em82mpqamS24MAACgKwjqY8SZM2fq7rvvVlVVlZqbm31uBC0AAID/CCps1dTUaN68eUpMTAx1PwAAAF1KUGHrW9/6lgoKCkLcCgAAQNdjMx394GE7zp07p7vvvltXXXWVhg0bJrvd7rN+1qxZIWvwSuDxeOR0OuV2u5WQkNDZ7QAAAD/4e/wO6gL5P/3pT3rppZcUGxurgoIC2Ww27zqbzUbYAgAA+LegwtbixYv16KOPKjs7WxERTNUFAADQnqCSUn19vSZPnkzQAgAAuIig0tK0adO0cePGUPcCAADQ5QT1MWJTU5NWrlypl156ScOHD291gfyqVatC0hwAAEC4CypslZaWasSIEZKkw4cP+6z75MXyAAAAV7qgwtbOnTtD3QcAAECXxBXuAAAAFgrqzNajjz7a4fof/ehHQTUDAADQ1QQVtrZs2eJzv6GhQRUVFYqKitLAgQMJWwAAAP8WVNgqKSlptczj8eg73/mO7rzzzktuCgAAoKsI2TVbCQkJevTRR/XDH/4wVEMCAACEvZBeIP/BBx/I7XaHckgAAICwFtTHiL/4xS987htjVFVVpeeee06TJk0KSWMAAABdQVBha/Xq1T73IyIidNVVV2natGlatGhRSBoDAADoCoIKWxUVFaHuAwAAoEsKKmxJ0vnz5/Xmm2+qpqZGzc3NPuvuuOOOS24MAACgKwgqbG3btk1ZWVk6ffp0q3U2m01NTU2X3BgAAEBXENS3EWfOnKl77rlHVVVVam5u9rkRtAAAAP4jqLBVU1OjefPmKTExMdT9AAAAdClBha1vfetbKigoCHErAAAAXY/NGGMC3ejcuXO6++67ddVVV2nYsGGy2+0+62fNmhWyBq8EHo9HTqdTbrdbCQkJnd0OAADwg7/H76AukP/Tn/6kl156SbGxsSooKJDNZvOus9lshC0AAIB/CypsLV68WI8++qiys7MVERHSX/wBAADoUoJKSvX19Zo8eTJBCwAA4CKCSkvTpk3Txo0bQ90LAABAlxPUx4hNTU1auXKlXnrpJQ0fPrzVBfKrVq0KSXMAAADhLqiwVVpaqhEjRkiSDh8+7LPukxfLAwAAXOmCCls7d+4MdR8AAABdEle4AwAAWMjvM1vz5s3TY489pvj4eM2bN6/DWq7ZAgAAuMDvsFVSUqKGhgbvfwMAAODigvq5HoQWP9cDAED48ff4HdQ1Wy+//HK765555plghgQAAOiSggpbX/va1/Tggw+qvr7eu+y9997T7bffrkWLFoWsOQAAgHAXVNh69dVX9Ze//EVf/OIXdeTIEf31r3/V0KFD9eGHH+rQoUOh7hEAACBsBRW2brjhBpWUlGj48OFKS0vTnXfeqQcffFA7duxQ3759Q90jAABA2Ap6nq2jR4/qjTfeUJ8+fRQVFaW///3vOnfuXCh7AwAACHtBha2cnByNHj1at9xyiw4fPqw33njDe6arqKgo1D0CAACEraDC1pNPPqk///nP+uUvf6mYmBhde+21ev3113XXXXdp4sSJAY+3Zs0apaSkKCYmRmlpadq1a1eH9YWFhUpLS1NMTIwGDBigdevWtarJy8tTamqqoqOjlZqaqi1btgS8X2OMli5dquTkZMXGxmrixIk6cuRIq3GKior05S9/WfHx8frc5z6niRMn6uOPPw7wWQAAAF2SCcJ7773X7rqCgoKAxtqwYYOx2+1m/fr1pqyszMyePdvEx8ebd955p83648ePm7i4ODN79mxTVlZm1q9fb+x2u9m0aZO3Zs+ePSYyMtIsW7bMlJeXm2XLlpmoqCizd+/egPabk5NjunfvbvLy8kxpaamZPHmySUpKMh6Px2dfCQkJZvny5ebw4cPmrbfeMi+88II5f/6838+B2+02kozb7Q7kqQMAAJ3I3+N3p09qesMNN2jkyJFau3atd9mQIUP0jW98Q8uXL29Vv3DhQm3dulXl5eXeZTNmzNChQ4e8H2FOnjxZHo9HL774ordm0qRJ6tGjh55//nm/9muMUXJysubMmaOFCxdKkurq6pSYmKgVK1Zo+vTpkqQbb7xRt9xyix577DG/H3NdXZ3q6uq89z0ej/r27cukpgAAhBFLJzVtUVZWpm3btmnr1q0+N3/V19dr//79ysjI8FmekZGhPXv2tLlNUVFRq/rMzEzt27fP+3NC7dW0jOnPfisqKlRdXe1TEx0drQkTJnhrampqVFxcrKuvvlpjxoxRYmKiJkyYoN27d3f4uJcvXy6n0+m98Q1OAAC6Lr9/G/GTjh8/rjvvvFOlpaWy2WxqOTlms9kkSU1NTX6NU1tbq6amJiUmJvosT0xMVHV1dZvbVFdXt1nf2Nio2tpaJSUltVvTMqY/+235s62ad955x/s8SNLSpUv1+OOP6/rrr9fvf/973XzzzTp8+LC+8IUvtPkYFi1a5PNj3i1ntgAAQNcT1Jmt2bNnKyUlRf/6178UFxenI0eO6NVXX1V6eroKCgoCHq8lpLUwxrRadrH6Ty/3Z8xLrWlubpYkTZ8+Xffdd59GjBih1atXa9CgQfr1r3/dbv/R0dFKSEjwuQEAgK4pqLBVVFSkRx99VFdddZUiIiIUERGhsWPHavny5Zo1a5bf4/Tu3VuRkZGtzmLV1NS0OqPUwuVytVkfFRWlXr16dVjTMqY/+3W5XJLUYU1SUpIkKTU11admyJAhOnHiRAePHAAAXCmCCltNTU3q1q2bpAvB5dSpU5Kkfv366ejRo36P43A4lJaWpvz8fJ/l+fn5GjNmTJvbjB49ulX99u3blZ6eLrvd3mFNy5j+7DclJUUul8unpr6+XoWFhd6a/v37Kzk5udVjfuutt9SvXz+/ngMAANDFBfNVx7Fjx5otW7YYY4yZMmWKmTRpktm9e7eZOnWqufbaawMaq2UKhtzcXFNWVmbmzJlj4uPjTWVlpTHGmOzsbJOVleWtb5n6Ye7cuaasrMzk5ua2mvrhtddeM5GRkSYnJ8eUl5ebnJycdqd+aG+/xlyY+sHpdJrNmzeb0tJSM2XKlFZTP6xevdokJCSYF154wfzjH/8wixcvNjExMebtt9/2+zlg6gcAAMKPv8fvoMLWtm3bTF5enjHGmGPHjpkhQ4YYm81mevfubV555ZWAx3v66adNv379jMPhMCNHjjSFhYXeddOmTTMTJkzwqS8oKDAjRowwDofD9O/f36xdu7bVmC+88IIZNGiQsdvtZvDgwd5+/d2vMcY0NzebJUuWGJfLZaKjo8348eNNaWlpq3GWL19u+vTpY+Li4szo0aPNrl27Anr8hC0AAMLPZz7P1pkzZ9SjR48OL2xH2/ydpwMAAFw+/D1+BzX1gySdP39eb775pmpqarzfymtxxx13BDssAABAlxJU2Nq2bZuysrJ0+vTpVutsNpvf82wBAAB0dUF9G3HmzJm65557VFVVpebmZp8bQQsAAOA/ggpbNTU1mjdvXrtzYQEAAOCCoMLWt771raBmigcAALjSBPVtxHPnzunuu+/WVVddpWHDhnknE20RyCzy4NuIAACEI0u/jfinP/1JL730kmJjY1VQUNDqNwkJWwAAABcEFbYWL16sRx99VNnZ2YqICOqTSAAAgCtCUEmpvr5ekydPJmgBAABcRFBpadq0adq4cWOoewEAAOhygvoYsampSStXrtRLL72k4cOHt7pAftWqVSFpDgAAINwFFbZKS0s1YsQISdLhw4d91vHbiAAAAP8RVNjauXNnqPsAAADokrjCHQAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBCQX0bUZI++OAD5ebmqry8XDabTUOGDNH9998vp9MZyv4AAADCWlBntvbt26eBAwdq9erVOnPmjGpra7V69WoNHDhQBw4cCHWPAAAAYctmjDGBbjRu3Dhdc801Wr9+vaKiLpwca2xs1He/+10dP35cr776asgb7co8Ho+cTqfcbrcSEhI6ux0AAOAHf4/fQYWt2NhYlZSUaPDgwT7Ly8rKlJ6ernPnzgXe8RWMsAUAQPjx9/gd1MeICQkJOnHiRKvlJ0+eVPfu3YMZEgAAoEsKKmxNnjxZ999/vzZu3KiTJ0/q3Xff1YYNG/Td735XU6ZMCXWPAAAAYSuobyM+/vjjstlsmjp1qhobGyVJdrtdDzzwgHJyckLaIAAAQDgL6pqtFufOndOxY8dkjNE111yjuLi4UPZ2xeCaLQAAwo+l12ydOHFCxhjFxcVp2LBhGj58uDdotXUtFwAAwJUqqLCVkpKi9957r9Xy06dPKyUl5ZKbAgAA6CqCClvGGNlstlbLP/zwQ8XExFxyUwAAAF1FQBfIz5s3T5Jks9n0wx/+0OcaraamJhUXF+v6668PaYMAAADhLKCwVVJSIunCma3S0lI5HA7vOofDoeuuu07z588PbYcAAABhLKCwtXPnTknSfffdpyeffJJvzgEAAFxEUPNs/eY3vwl1HwAAAF1SUBfIAwAAwD+ELQAAAAsRtgAAACxE2AIAALBQ0GFr165duvfeezV69Gj985//lCQ999xz2r17d8iaAwAACHdBha28vDxlZmYqNjZWJSUlqqurkySdPXtWy5YtC2mDAAAA4SyosPWTn/xE69at0/r162W3273Lx4wZowMHDoSsOQAAgHAXVNg6evSoxo8f32p5QkKCPvjgg0vtCQAAoMsIKmwlJSXp7bffbrV89+7dGjBgwCU3BQAA0FUEFbamT5+u2bNnq7i4WDabTadOndIf//hHzZ8/X9///vdD3SMAAEDYCurnehYsWCC3262bbrpJ58+f1/jx4xUdHa358+dr5syZoe4RAAAgbNmMMSbYjc+dO6eysjI1NzcrNTVV3bp1C2VvVwyPxyOn0ym3282PewMAECb8PX4H9THifffdp1deeUWxsbFKT0/XqFGjCFoAAABtCCpsnT59Wl/72tfUp08fPfjggzp48GCI2wIAAOgaggpbW7duVXV1tZYsWaL9+/crLS1NqampWrZsmSorK0PcIgAAQPi6pGu2Wrz77rt6/vnn9etf/1r/+Mc/1NjYGIrerhhcswUAQPix9JqtT2poaNC+fftUXFysyspKJSYmXuqQAAAAXUbQYWvnzp363ve+p8TERE2bNk3du3fXX/7yF508eTKU/QEAAIS1oObZ6tOnj06fPq3MzEw988wzuv322xUTExPq3gAAAMJeUGHrRz/6ke6++2716NEj1P0AAAB0KUGFrf/5n/8JdR8AAABdkt9ha968eXrssccUHx+vefPmdVi7atWqS24MAACgK/A7bJWUlKihocH73+2x2WyX3hUAAEAX4fe3EXfu3KnPfe5zkqTf/e53euWVV7Rz506f244dO/Tb3/424CbWrFmjlJQUxcTEKC0tTbt27eqwvrCwUGlpaYqJidGAAQO0bt26VjV5eXlKTU1VdHS0UlNTtWXLloD3a4zR0qVLlZycrNjYWE2cOFFHjhxpsydjjG699VbZbDb9+c9/9v/BAwCALi2oqR9SUlJUW1vbavmZM2eUkpIS0FgbN27UnDlz9Mgjj6ikpETjxo3TrbfeqhMnTrRZX1FRoa9+9asaN26cSkpK9PDDD2vWrFnKy8vz1hQVFWny5MnKysrSoUOHlJWVpXvuuUfFxcUB7XflypVatWqVnnrqKb3xxhtyuVy65ZZbdPbs2VZ9PfHEE5zVAwAArZkg2Gw2869//avV8srKShMXFxfQWKNGjTIzZszwWTZ48GCTnZ3dZv2CBQvM4MGDfZZNnz7d3Hjjjd7799xzj5k0aZJPTWZmpvn2t7/t936bm5uNy+UyOTk53vXnz583TqfTrFu3zme7gwcPmj59+piqqiojyWzZsuUij9qX2+02kozb7Q5oOwAA0Hn8PX4H9G3ElgvjbTabfvSjHykuLs67rqmpScXFxbr++uv9Hq++vl779+9Xdna2z/KMjAzt2bOnzW2KioqUkZHhsywzM1O5ublqaGiQ3W5XUVGR5s6d26rmiSee8Hu/FRUVqq6u9tlXdHS0JkyYoD179mj69OmSpHPnzmnKlCl66qmn5HK5/HrcdXV1qqur8973eDx+bQcAAMJPQGGr5cJ4Y4xKS0vlcDi86xwOh6677jrNnz/f7/Fqa2vV1NTU6id+EhMTVV1d3eY21dXVbdY3NjaqtrZWSUlJ7da0jOnPflv+bKvmnXfe8d6fO3euxowZo69//ev+PmwtX75cP/7xj/2uBwAA4SugsLVz505J0n333acnn3wyZD+a/OlrnYwxHV7/1Fb9p5f7M+al1mzdulU7duzo8NuZbVm0aJHP9Bkej0d9+/YNaAwAABAegrpA/je/+U1Iglbv3r0VGRnZ6ixWTU1Nuz9o7XK52qyPiopSr169OqxpGdOf/bZ8JNhRzY4dO3Ts2DF97nOfU1RUlKKiLmTXb37zm5o4cWK7jzs6OloJCQk+NwAA0DUF/UPUklRWVqZt27Zp69atPjd/ORwOpaWlKT8/32d5fn6+xowZ0+Y2o0ePblW/fft2paeny263d1jTMqY/+01JSZHL5fKpqa+vV2FhobcmOztbb775pg4ePOi9SdLq1av1m9/8xu/nAQAAdGHBXH1/7NgxM3z4cGOz2UxERISx2Wze/46IiAhorA0bNhi73W5yc3NNWVmZmTNnjomPjzeVlZXGGGOys7NNVlaWt/748eMmLi7OzJ0715SVlZnc3Fxjt9vNpk2bvDWvvfaaiYyMNDk5Oaa8vNzk5OSYqKgos3fvXr/3a4wxOTk5xul0ms2bN5vS0lIzZcoUk5SUZDweT7uPR3wbEQCAK4K/x++gwtZtt91mvv71r5uamhrTrVs3U1ZWZnbt2mVGjRplXn311YDHe/rpp02/fv2Mw+EwI0eONIWFhd5106ZNMxMmTPCpLygoMCNGjDAOh8P079/frF27ttWYL7zwghk0aJCx2+1m8ODBJi8vL6D9GnNh+oclS5YYl8tloqOjzfjx401paWmHj4WwBQDAlcHf47fNmH9fXR6A3r17a8eOHRo+fLicTqdef/11DRo0SDt27NCDDz4Y8AXjVzqPxyOn0ym32831WwAAhAl/j99BXbPV1NSkbt26SboQvE6dOiVJ6tevn44ePRrMkAAAAF1SQFM/tBg6dKjefPNNDRgwQDfccINWrlwph8OhX/3qVxowYECoewQAAAhbQYWtxYsX66OPPpIk/eQnP9Ftt92mcePGqVevXtq4cWNIGwQAAAhnQV2z1ZYzZ86oR48e/BhzELhmCwCA8OPv8TuoM1tt6dmzZ6iGAgAA6DL8Dluf/HmZi1m1alVQzQAAAHQ1foctf6dz4GNEAACA//A7bLX8CDUAAAD8F/RvI+7atUv33nuvxowZo3/+85+SpOeee067d+8OWXMAAADhLqiwlZeXp8zMTMXGxurAgQOqq6uTJJ09e1bLli0LaYMAAADhLKiw9ZOf/ETr1q3T+vXrZbfbvcvHjBmjAwcOhKw5AACAcBdU2Dp69KjGjx/fanlCQoI++OCDS+0JAACgywgqbCUlJentt99utXz37t38XA8AAMAnBBW2pk+frtmzZ6u4uFg2m02nTp3SH//4R82fP1/f//73Q90jAABA2ApqBvkFCxbI7Xbrpptu0vnz5zV+/HhFR0dr/vz5mjlzZqh7BAAACFuX9NuI586dU1lZmZqbm5Wamqpu3bqFsrcrBr+NCABA+PH3+B3wx4gNDQ266aab9NZbbykuLk7p6ekaNWoUQQsAAKANAYctu92uw4cP87M8AAAAfgjqAvmpU6cqNzc31L0AAAB0OUFdIF9fX69nn31W+fn5Sk9PV3x8vM/6VatWhaQ5AACAcBdU2Dp8+LBGjhwpSXrrrbd81vHxIgAAwH8EFbZ27twZ6j4AAAC6pKCu2QIAAIB/CFsAAAAWImwBAABYiLAFAABgoUuaQR4AAAAdYwZ5AAAACzGDPAAAgIWYQR4AAMBCzCAPAABgIWaQBwAAsBBTPwAAAFgoqDNbkvTBBx8oNzdX5eXlstlsGjJkiO6//345nc5Q9gcAABDWgjqztW/fPg0cOFCrV6/WmTNnVFtbq9WrV2vgwIE6cOBAqHsEAAAIWzZjjAl0o3Hjxumaa67R+vXrFRV14eRYY2Ojvvvd7+r48eN69dVXQ95oV+bxeOR0OuV2u5WQkNDZ7QAAAD/4e/wOKmzFxsaqpKREgwcP9lleVlam9PR0nTt3LvCOr2CELQAAwo+/x++gPkZMSEjQiRMnWi0/efKkunfvHsyQAAAAXVJQYWvy5Mm6//77tXHjRp08eVLvvvuuNmzYoO9+97uaMmVKqHsEAAAIW0F9G/Hxxx+XzWbT1KlT1djYKOnCbyY+8MADysnJCWmDAAAA4SzgsNXQ0KDMzEw988wzWr58uY4dOyZjjK655hrFxcVZ0SMAAEDYCjhs2e12HT58WDabTXFxcRo2bJgVfQEAAHQJQV2zNXXqVOXm5oa6FwAAgC4nqGu26uvr9eyzzyo/P1/p6emKj4/3Wb9q1aqQNAcAABDuggpbhw8f1siRIyVJb731ls86m8126V0BAAB0EUGFrZ07d4a6DwAAgC4p4Gu2GhoadNNNN7U6owUAAIDWAg5bn/w2IgAAADrGtxEBAAAsxLcRAQAALMS3EQEAACzEtxEBAAAsFNQ1W5K0a9cu3XvvvRozZoz++c9/SpKee+457d69O2TNAQAAhLugwlZeXp4yMzMVGxurAwcOqK6uTpJ09uxZLVu2LKQNAgAAhLOgwtZPfvITrVu3TuvXr5fdbvcuHzNmjA4cOBCy5gAAAMJdUGHr6NGjGj9+fKvlCQkJ+uCDDy61JwAAgC4jqLCVlJSkt99+u9Xy3bt3a8CAAQGPt2bNGqWkpCgmJkZpaWnatWtXh/WFhYVKS0tTTEyMBgwYoHXr1rWqycvLU2pqqqKjo5WamqotW7YEvF9jjJYuXark5GTFxsZq4sSJOnLkiHf9mTNn9IMf/ECDBg1SXFycPv/5z2vWrFlyu90BPwcAAKCLMkFYsWKFSU1NNXv37jXdu3c3u3btMn/4wx/MVVddZX75y18GNNaGDRuM3W4369evN2VlZWb27NkmPj7evPPOO23WHz9+3MTFxZnZs2ebsrIys379emO3282mTZu8NXv27DGRkZFm2bJlpry83CxbtsxERUWZvXv3BrTfnJwc0717d5OXl2dKS0vN5MmTTVJSkvF4PMYYY0pLS81dd91ltm7dat5++23zyiuvmC984Qvmm9/8ZkDPgdvtNpKM2+0OaDsAANB5/D1+BxW2jDHm4YcfNrGxscZmsxmbzWZiYmLM4sWLAx5n1KhRZsaMGT7LBg8ebLKzs9usX7BggRk8eLDPsunTp5sbb7zRe/+ee+4xkyZN8qnJzMw03/72t/3eb3Nzs3G5XCYnJ8e7/vz588bpdJp169a1+3j+7//+zzgcDtPQ0NBuzacRtgAACD/+Hr+Dnvrhpz/9qWpra/X6669r7969eu+99/TYY48FNEZ9fb3279+vjIwMn+UZGRnas2dPm9sUFRW1qs/MzNS+ffvU0NDQYU3LmP7st6KiQtXV1T410dHRmjBhQru9SZLb7VZCQoKiotqfwqyurk4ej8fnBgAAuqagw5YkxcXFKT09XaNGjVK3bt0C3r62tlZNTU1KTEz0WZ6YmKjq6uo2t6murm6zvrGxUbW1tR3WtIzpz35b/gykt9OnT+uxxx7T9OnT233MkrR8+XI5nU7vrW/fvh3WAwCA8HVJYStUPv0TP8aYDn/2p636Ty/3Z8xQ1UiSx+PR1772NaWmpmrJkiXt9i5JixYtktvt9t5OnjzZYT0AAAhfQf1cT6j07t1bkZGRrc4U1dTUtDqj1MLlcrVZHxUVpV69enVY0zKmP/t1uVySLpzhSkpK6rC3s2fPatKkSerWrZu2bNniM/dYW6KjoxUdHd1hDQAA6Bo69cyWw+FQWlqa8vPzfZbn5+drzJgxbW4zevToVvXbt29Xenq6N+S0V9Mypj/7TUlJkcvl8qmpr69XYWGhT28ej0cZGRlyOBzaunWrYmJiAnkKAABAV2f9tfoda5mCITc315SVlZk5c+aY+Ph4U1lZaYwxJjs722RlZXnrW6Z+mDt3rikrKzO5ubmtpn547bXXTGRkpMnJyTHl5eUmJyen3akf2tuvMRemfnA6nWbz5s2mtLTUTJkyxWfqB4/HY2644QYzbNgw8/bbb5uqqirvrbGx0e/ngG8jAgAQfiyf+iGUnn76adOvXz/jcDjMyJEjTWFhoXfdtGnTzIQJE3zqCwoKzIgRI4zD4TD9+/c3a9eubTXmCy+8YAYNGmTsdrsZPHiwycvLC2i/xlyY/mHJkiXG5XKZ6OhoM378eFNaWupdv3PnTiOpzVtFRYXfj5+wBQBA+PH3+G0z5t9Xl6PTeDweOZ1O77QRAADg8ufv8fuy+DYiAABAV0XYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBCUZ3dAKzR1Gz0esUZ1Zw9r6u7x2hUSk9FRtg6uy0AAD4zl8ux8LI4s7VmzRqlpKQoJiZGaWlp2rVrV4f1hYWFSktLU0xMjAYMGKB169a1qsnLy1Nqaqqio6OVmpqqLVu2BLxfY4yWLl2q5ORkxcbGauLEiTpy5IhPTV1dnX7wgx+od+/eio+P1x133KF33303iGchdLYdrtLYFTs0Zf1ezd5wUFPW79XYFTu07XBVp/YFAMBn5XI6FnZ62Nq4caPmzJmjRx55RCUlJRo3bpxuvfVWnThxos36iooKffWrX9W4ceNUUlKihx9+WLNmzVJeXp63pqioSJMnT1ZWVpYOHTqkrKws3XPPPSouLg5ovytXrtSqVav01FNP6Y033pDL5dItt9yis2fPemvmzJmjLVu2aMOGDdq9e7c+/PBD3XbbbWpqarLg2bq4bYer9MAfDqjKfd5nebX7vB74wwECFwCgy7vcjoU2Y4z5TPf4KTfccINGjhyptWvXepcNGTJE3/jGN7R8+fJW9QsXLtTWrVtVXl7uXTZjxgwdOnRIRUVFkqTJkyfL4/HoxRdf9NZMmjRJPXr00PPPP+/Xfo0xSk5O1pw5c7Rw4UJJF85iJSYmasWKFZo+fbrcbreuuuoqPffcc5o8ebIk6dSpU+rbt6/+9re/KTMz06/nwOPxyOl0yu12KyEhwd+nrpWmZqOxK3a0+svVwibJ5YzR7oVf5iNFAECX9FkeC/09fnfqma36+nrt379fGRkZPsszMjK0Z8+eNrcpKipqVZ+Zmal9+/apoaGhw5qWMf3Zb0VFhaqrq31qoqOjNWHCBG/N/v371dDQ4FOTnJysoUOHttu/dCG0eTwen1sovF5xpt2/XJJkJFW5z+v1ijMh2R8AAJeby/FY2Klhq7a2Vk1NTUpMTPRZnpiYqOrq6ja3qa6ubrO+sbFRtbW1Hda0jOnPflv+vFiNw+FQjx49/O5fkpYvXy6n0+m99e3bt93aQNScbf8vVzB1AACEm8vxWNjp12xJks3mexrPGNNq2cXqP73cnzFDVfNpF6tZtGiR3G6393by5MkOx/PX1d1jQloHAEC4uRyPhZ0atnr37q3IyMhWZ4FqampanVFq4XK52qyPiopSr169OqxpGdOf/bpcLkm6aE19fb3ef/99v/uXLnwcmZCQ4HMLhVEpPZXkjFF7Mc8mKcl54auvAAB0RZfjsbBTw5bD4VBaWpry8/N9lufn52vMmDFtbjN69OhW9du3b1d6errsdnuHNS1j+rPflJQUuVwun5r6+noVFhZ6a9LS0mS3231qqqqqdPjw4Xb7t1JkhE1Lbk+VpFZ/yVruL7k9lYvjAQBd1mV5LDSdbMOGDcZut5vc3FxTVlZm5syZY+Lj401lZaUxxpjs7GyTlZXlrT9+/LiJi4szc+fONWVlZSY3N9fY7XazadMmb81rr71mIiMjTU5OjikvLzc5OTkmKirK7N271+/9GmNMTk6OcTqdZvPmzaa0tNRMmTLFJCUlGY/H462ZMWOG6dOnj3n55ZfNgQMHzJe//GVz3XXXmcbGRr+fA7fbbSQZt9sd1HP4aS+WnjI3LnvZ9Fv4/3lvNy572bxYeiok4wMAcLn7LI6F/h6/Oz1sGWPM008/bfr162ccDocZOXKkKSws9K6bNm2amTBhgk99QUGBGTFihHE4HKZ///5m7dq1rcZ84YUXzKBBg4zdbjeDBw82eXl5Ae3XGGOam5vNkiVLjMvlMtHR0Wb8+PGmtLTUp+bjjz82M2fOND179jSxsbHmtttuMydOnAjo8Yc6bBljTGNTs9nzdq35c8m7Zs/btaaxqTlkYwMAEA6sPhb6e/zu9Hm2ELp5tgAAwGcnLObZAgAA6OoIWwAAABYibAEAAFiIsAUAAGAhwhYAAICFCFsAAAAWImwBAABYiLAFAABgIcIWAACAhaI6uwFILZP4ezyeTu4EAAD4q+W4fbEf4yFsXQbOnj0rSerbt28ndwIAAAJ19uxZOZ3Odtfz24iXgebmZp06dUrdu3eXzWYL2bgej0d9+/bVyZMn+c3FMMDrFT54rcILr1f4CLfXyhijs2fPKjk5WRER7V+ZxZmty0BERIT69Olj2fgJCQlh8ZcWF/B6hQ9eq/DC6xU+wum16uiMVgsukAcAALAQYQsAAMBChK0uLDo6WkuWLFF0dHRntwI/8HqFD16r8MLrFT666mvFBfIAAAAW4swWAACAhQhbAAAAFiJsAQAAWIiwBQAAYCHCVphbs2aNUlJSFBMTo7S0NO3atavD+sLCQqWlpSkmJkYDBgzQunXrPqNOIQX2ehUUFMhms7W6/f3vf/8MO74yvfrqq7r99tuVnJwsm82mP//5zxfdhvdW5wj0teJ91XmWL1+uL37xi+revbuuvvpqfeMb39DRo0cvul1XeG8RtsLYxo0bNWfOHD3yyCMqKSnRuHHjdOutt+rEiRNt1ldUVOirX/2qxo0bp5KSEj388MOaNWuW8vLyPuPOr0yBvl4tjh49qqqqKu/tC1/4wmfU8ZXro48+0nXXXaennnrKr3reW50n0NeqBe+rz15hYaH+93//V3v37lV+fr4aGxuVkZGhjz76qN1tusx7yyBsjRo1ysyYMcNn2eDBg012dnab9QsWLDCDBw/2WTZ9+nRz4403WtYj/iPQ12vnzp1Gknn//fc/g+7QHklmy5YtHdbw3ro8+PNa8b66fNTU1BhJprCwsN2arvLe4sxWmKqvr9f+/fuVkZHhszwjI0N79uxpc5uioqJW9ZmZmdq3b58aGhos6xXBvV4tRowYoaSkJN18883auXOnlW0iSLy3wg/vq87ndrslST179my3pqu8twhbYaq2tlZNTU1KTEz0WZ6YmKjq6uo2t6murm6zvrGxUbW1tZb1iuBer6SkJP3qV79SXl6eNm/erEGDBunmm2/Wq6+++lm0jADw3gofvK8uD8YYzZs3T2PHjtXQoUPbresq762ozm4Al8Zms/ncN8a0Wnax+raWwxqBvF6DBg3SoEGDvPdHjx6tkydP6vHHH9f48eMt7ROB470VHnhfXR5mzpypN998U7t3775obVd4b3FmK0z17t1bkZGRrc6K1NTUtPq/gBYul6vN+qioKPXq1cuyXhHc69WWG2+8Uf/4xz9C3R4uEe+t8Mb76rP1gx/8QFu3btXOnTvVp0+fDmu7ynuLsBWmHA6H0tLSlJ+f77M8Pz9fY8aMaXOb0aNHt6rfvn270tPTZbfbLesVwb1ebSkpKVFSUlKo28Ml4r0V3nhffTaMMZo5c6Y2b96sHTt2KCUl5aLbdJn3VidenI9LtGHDBmO3201ubq4pKyszc+bMMfHx8aaystIYY0x2drbJysry1h8/ftzExcWZuXPnmrKyMpObm2vsdrvZtGlTZz2EK0qgr9fq1avNli1bzFtvvWUOHz5ssrOzjSSTl5fXWQ/hinH27FlTUlJiSkpKjCSzatUqU1JSYt555x1jDO+ty0mgrxXvq87zwAMPGKfTaQoKCkxVVZX3du7cOW9NV31vEbbC3NNPP2369etnHA6HGTlypM9XaKdNm2YmTJjgU19QUGBGjBhhHA6H6d+/v1m7du1n3PGVLZDXa8WKFWbgwIEmJibG9OjRw4wdO9b89a9/7YSurzwt0wN8+jZt2jRjDO+ty0mgrxXvq87T1uskyfzmN7/x1nTV95bNmH9faQYAAICQ45otAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAsRtgAAACxE2AIAALAQYQsAAMBChC0AAAALEbYAAAAsRNgC0CVNnDhRc+bMafd+oNt3Bn96uBz6BNCxqM5uAAA+C5s3b5bdbu/sNgLy6Z4nTpyo66+/Xk888US7NQAuP4QtAJet+vp6ORyOkIzVs2fPkIzzWfKn5858XI2NjYqK4jACXAwfIwK4bEycOFEzZ87UvHnz1Lt3b91yyy3atm2bxo4dq8997nPq1auXbrvtNh07dsxnu48++khTp05Vt27dlJSUpJ///Odtjt3ycZs/Y/rb68yZM73jLF68WMYYb01dXZ1mzZqlq6++WjExMRo7dqzeeOMN7/pNmzZp2LBhio2NVa9evfSVr3xFH330UZs9f+c731FhYaGefPJJ2Ww22Ww2VVZWtvoY8WL7nDhxombNmqUFCxaoZ8+ecrlcWrp06UUfb2VlpWw2mzZt2qTx48crOjpaW7ZsCeg5A65UhC0Al5Xf/e53ioqK0muvvaZnnnlGH330kebNm6c33nhDr7zyiiIiInTnnXequbnZu81DDz2knTt3asuWLdq+fbsKCgq0f//+dvfhz5iB9FpcXKxf/OIXWr16tZ599lnv+gULFigvL0+/+93vdODAAV1zzTXKzMzUmTNnVFVVpSlTpui///u/VV5eroKCAt11110+Ye2TnnzySY0ePVrf+973VFVVpaqqKvXt27dVXUf7/GTf8fHxKi4u1sqVK/Xoo48qPz+/w8d68OBBSdKKFSv0wx/+UEeOHFFGRkZAzxdwxTIAcJmYMGGCuf766zusqampMZJMaWmpMcaYs2fPGofDYTZs2OCtOX36tImNjTWzZ8/2GfuT9zsa82L1LeuHDBlimpubvcsWLlxohgwZYowx5sMPPzR2u9388Y9/9K6vr683ycnJZuXKlWb//v1GkqmsrOxwHxd7DJ9cdrF9ttSPHTvWZ4wvfvGLZuHChe32YYwxS5cuNfHx8aaioqLDOgCtcWYLwGUlPT3d5/6xY8f0X//1XxowYIASEhKUkpIiSTpx4oR3fX19vUaPHu3dpmfPnho0aFC7+7jYmP668cYbZbPZvPdHjx6tf/zjH2pqatKxY8fU0NCgL33pS971drtdo0aNUnl5ua677jrdfPPNGjZsmO6++26tX79e77//fkD7b+txdbTPFsOHD/fZLikpSTU1NR2OffDgQd1xxx3q37//JfUIXIkIWwAuK/Hx8T73b7/9dp0+fVrr169XcXGxiouLJV24eF5Sux+7deRiY4ZCS1+fDGMty202myIjI5Wfn68XX3xRqamp+uUvf6lBgwapoqLCsn22+PS3F20220U/Qj106JAmTpwYdG/AlYywBeCydfr0aZWXl2vx4sW6+eabNWTIkFZnf6655hrZ7Xbt3bvXu+z999/XW2+9FfSY/vrkPlvuf+ELX1BkZKSuueYaORwO7d6927u+oaFB+/bt05AhQyRdCDlf+tKX9OMf/1glJSVyOBwdXnTucDjU1NTU7np/9hkMj8ejyspKjRgxIugxgCsZ39kFcNnq0aOHevXqpV/96ldKSkrSiRMnlJ2d7VPTrVs33X///XrooYfUq1cvJSYm6pFHHlFERNv/L+nPmP46efKk5s2bp+nTp+vAgQP65S9/6f0mZHx8vB544AE99NBD6tmzpz7/+c9r5cqVOnfunO6//34VFxfrlVdeUUZGhq6++moVFxfrvffe6zAU9e/fX8XFxaqsrFS3bt1aTftwsX0G69ChQ4qIiNCwYcOCHgO4khG2AFy2IiIitGHDBs2aNUtDhw7VoEGD9Itf/KLVx1k/+9nP9OGHH+qOO+5Q9+7d9eCDD8rtdl/SmP6YOnWqPv74Y40aNUqRkZH6wQ9+oP/5n//xrs/JyVFzc7OysrJ09uxZpaen66WXXlKPHj2UkJCgV199VU888YQ8Ho/69eunn//857r11lvb3d/8+fM1bdo0paam6uOPP27zI8eO9hmsQ4cOafDgwYqJiQl6DOBKZjPBXPAAAFe4tmZzB4C2cM0WAACAhQhbAAAAFuJjRAAAAAtxZgsAAMBChC0AAAALEbYAAAAsRNgCAACwEGELAADAQoQtAAAACxG2AAAALETYAgAAsBBhCwAAwEKELQAAAAv9/5pkjz9feXmqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(rnum, norm_err, 'o')\n",
    "plt.ylabel('error relative to maximum $T$')\n",
    "plt.xlabel('radial position $r$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ True  True  True  True  True  True  True  True  True  True  True  True\n",
      "  True  True  True  True  True  True  True  True  True  True  True  True\n",
      "  True  True  True  True  True  True  True  True  True  True  True  True\n",
      "  True  True  True False False False False False False False False False\n",
      " False False False  True]\n"
     ]
    }
   ],
   "source": [
    "print(np.abs(norm_err) < 0.0001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.all(np.abs(norm_err) < 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 4
}