{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6ec88f25",
   "metadata": {},
   "source": [
    "# Poisson equation example\n",
    "\n",
    "- Written by Ali A. Eftekhari (last update 2021)\n",
    "- Ported to Python and full rewrite by Gavin M. Weir June, 2023.\n",
    " \n",
    "The generalized form of the equations solved in this package looks like, \n",
    "\n",
    "$$ \\alpha \\frac{\\partial \\varphi}{\\partial t} + \\nabla\\cdot\\left(\\vec{u}\\varphi\\right) + \\nabla\\cdot\\left(-D\\nabla\\varphi\\right) + \\beta \\varphi = \\gamma $$\n",
    "\n",
    "with boundary condition,\n",
    "\n",
    "$$ a\\nabla\\varphi\\cdot \\vec{e} + b\\varphi = c $$.\n",
    "\n",
    "\n",
    "The inhomogeneous Poisson equation\n",
    "\n",
    "$$ \\nabla^2 \\varphi + s\\left(\\vec{x}\\right) = 0 $$\n",
    "\n",
    "can be generalized to a simple 1D case,\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\frac{\\partial^2 \\varphi}{\\partial ^2 x}  + s\\left(x\\right) &= 0 \\\\\n",
    "    \\varphi\\left(x_L\\right) &= 0 \\\\\n",
    "    \\frac{\\partial \\varphi}{\\partial x}|_{x_R} &= 0\n",
    "\\end{align}    \n",
    "$$\n",
    "\n",
    "The corresponding equation in our form has \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    D = 1.0, \\vec{u} &= \\vec{0}  \\\\\n",
    "    \\alpha = 0, \\beta = 0, \\gamma &= s  \\\\\n",
    "    a_L = 0, b_L = 1, c_L &= 0 \\\\\n",
    "    a_R = 1, b_R = 0, c_R &= 0 \\\\    \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "see this link\n",
    "http://scicomp.stackexchange.com/questions/8577/peculiar-error-when-solving-the-poisson-equation-on-a-non-uniform-mesh-1d-only\n",
    "\n",
    "Strange behavior when change the number of grids from even to odd\n",
    "Wrong results does not always mean that the code has bugs.\n",
    "\n",
    "Wrong use of the code can also give you wrong results.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a57b83ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f43c27c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# explicit imports as an alternative to `import pyfvtool as pf`\n",
    "from pyfvtool import Grid1D\n",
    "from pyfvtool import cellLocations, CellVariable\n",
    "from pyfvtool import faceLocations, FaceVariable\n",
    "from pyfvtool import BoundaryConditions\n",
    "from pyfvtool import diffusionTerm\n",
    "from pyfvtool import constantSourceTerm\n",
    "from pyfvtool import boundaryConditionsTerm\n",
    "from pyfvtool import solvePDE\n",
    "from pyfvtool import visualizeCells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b6982425-7882-4006-8fb1-8854f8caea51",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the domain and create a mesh structure\n",
    "L = 20      # domain length\n",
    "# Nx = 10000  # number of cells (original)\n",
    "Nx = 100000  # number of cells (test)\n",
    "m = Grid1D(Nx, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7cbde30b-9b4d-407a-829d-72d3da49ed8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution variable\n",
    "c = CellVariable(m, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fa10196f-b006-4423-af87-21b53b4e56a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Boundary Conditions\n",
    "\n",
    "# Left boundary\n",
    "c.BCs.left.a, c.BCs.left.b, c.BCs.left.c = 0.0, 1.0, 0.0\n",
    "\n",
    "# Right boundary (Neumann zero-flux is the default in this package)\n",
    "c.BCs.right.a, c.BCs.right.b, c.BCs.right.c = 1.0, 0.0, 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1c0bef95-54a6-4d8e-87d4-21f555fd8ab7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the transfer coeffs\n",
    "D_val = 1.0\n",
    "D = FaceVariable(m, D_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "13571a64-0013-4201-88d0-2fd93bdd8696",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define source term\n",
    "def rho(x):\n",
    "    return -1.0*((x>=-1.0)*(x<=0))+(x>0)*(x<=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7cb9682c-0140-42c5-8afd-f95ae26cb29c",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = m.cellcenters.x-0.5*L # shift the domain to [-10,10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ecb0458b-a78b-49d0-94e3-9eb8cdb66e83",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = solvePDE(c, [-constantSourceTerm(CellVariable(m, rho(x))),\n",
    "                  diffusionTerm(D)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ce6da8d0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGwCAYAAABb3Do8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXyklEQVR4nO3deVxUVeM/8M/IMizKKBIMJAGa4oILgrKYWypqYpmVmEXaVzFzy7AyygWxMq3UNrcy0fRRewK0J00Fxe0RDBXX1ExRUEEUdQZU9vP7w4f5OTIwMzjDMPB5v1739XLuPffMOXNl+HDuPfdKhBACRERERFSlRqZuABEREVFdx8BEREREpAUDExEREZEWDExEREREWjAwEREREWnBwERERESkBQMTERERkRaWpm5AfVFeXo5r166hSZMmkEgkpm4OERER6UAIgfz8fLi5uaFRo6rHkRiYDOTatWtwd3c3dTOIiIioBrKystCiRYsqtzMwGUiTJk0APPjAHRwcTNwaIiIi0oVSqYS7u7vq93hVGJgMpOI0nIODAwMTERGRmdF2OQ0v+iYiIiLSgoGJiIiISAsGJiIiIiItGJiIiIiItGBgIiIiItKCgYmIiIhICwYmIiIiIi0YmIiIiIi0YGAiIiIi0oKBiYiIiEgLswtM+/btw9ChQ+Hm5gaJRILNmzdr3Wfv3r3w8/ODjY0NWrZsieXLl1cqExcXh/bt20MqlaJ9+/ZISEgwQuuJiIjIHJldYLp79y46d+6M7777TqfyGRkZeO6559CzZ0+kp6fjo48+wtSpUxEXF6cqk5KSgrCwMISHh+P48eMIDw/HiBEjcOjQIWN1g4iIiMyIRAghTN2ImpJIJEhISMCwYcOqLDNjxgz89ttvOHPmjGrdhAkTcPz4caSkpAAAwsLCoFQq8ccff6jKDBo0CM2aNcOGDRt0aotSqYRMJoNCoeDDd4nqmZIS4No1Y9Rbgnv3ciGTlRq+cqJ6yNHREU2aNDFonbr+/rY06LvWQSkpKQgJCVFbN3DgQKxatQolJSWwsrJCSkoK3n333UpllixZUmW9RUVFKCoqUr1WKpUGbTcR1Q1CAP7+wIkTxqjdCoAcwAsAthrjDYjqlRUrVmD8+PEmee96H5hycnLg4uKits7FxQWlpaW4efMmXF1dqyyTk5NTZb3z58/H3LlzjdJmIqo7ior+f1iSSgGJ5PHqE6L8oT+2rAFYwNKyGywtdz1exUQNgIWFhcneu94HJuDBqbuHVZyFfHi9pjKPrntYVFQUIiMjVa+VSiXc3d0N0VwiqqNu3AAe52zAlStX4O/vj+vXr6N9+/Zo3XoXtmyRY/bsOZg1a47hGkpEBmd2F33rSy6XVxopys3NhaWlJZo3b15tmUdHnR4mlUrh4OCgthARVaWsrAyvv/46rl+/jo4dO2L//v2Qy+WmbhYR6ajeB6agoCAkJiaqrdu5cyf8/f1hZWVVbZng4OBaaycR1W8rVqzA3r17YW9vj7i4ODg6Opq6SUSkB7M7JVdQUIB//vlH9TojIwPHjh2Do6MjnnrqKURFReHq1atYu3YtgAcz4r777jtERkYiIiICKSkpWLVqldrst3feeQe9evXCggUL8MILL2DLli1ISkrCgQMHar1/RFT/3LlzB7NnzwYAfP7552jdurWJW0RE+jK7EabDhw/D19cXvr6+AIDIyEj4+vqqvoyys7ORmZmpKu/l5YVt27Zhz5496NKlC+bNm4dvvvkGL730kqpMcHAwNm7ciNWrV6NTp06IjY3Fpk2bEBAQULudI6J6aeHChcjLy0O7du0wYcKEStvN9+YuRA2HWd+HqS7hfZiI6qfCQsDW9sG/lUr9L/rOz8+Hu7s7FAoF4uPj8eKLL6q2TZgArFgBzJ0L/O9vPiKqZbr+/ja7ESYiInOyatUqKBQKeHt744UXXlDb9ri3KCCi2sPARERkJOXl5fj6668BANOnT0ejRvzKJTJX/OklIjKSPXv24NKlS5DJZHj99ddN3RwiegwMTERERrJ69WoAwMiRI2FbcSEUEZklBiYiomrUdFqMQqFAXFwcAODNN980ynsQUe1hYCIiMoLNmzfj/v37aNu2Lbp3727q5hDRY2JgIiLSkT6z2uLj4wEAYWFhVT6XkrPkiMwHAxMRkYEVFBRg586dAIDhw4ebuDVEZAgMTEREBrZ9+3YUFhaiVatW6Nixo6mbQ0QGwMBERGRgFafjXnzxxSpPxxGReWFgIiIyoNLSUvzxxx8AgGHDhum0D2fJEdV9DExERAZ0+PBh3LlzB02bNuUDvInqEQYmIiIDSkxMBAD069cPlpaW1Zbl2Toi88HARERkQBWz40JCQkzcEiIyJAYmIiIDUSqVSElJAQAMGDDAxK0hIkNiYCIiMpDk5GSUlZXh6aefhpeXl8778aJvorqPgYmIqBr6hJndu3cD4OgSUX3EwEREZCD79+8HAPTu3dvELSEiQ2NgIiLSUXWz2pRKJY4fPw4AeOaZZx67PiKqWxiYiIgMICUlBeXl5fDy8sKTTz5p6uYQkYExMBERGcCBAwcAAD179jRxS4jIGBiYiIgMoOL6JV1Pxz2Ms+SI6j4GJiKix1RcXIxDhw4B4AgTUX3FwERE9JiOHj2KwsJCODk5wdvb29TNISIjYGAiInpMFaNLQUFBkOgx9Y2z5IjMBwMTEdFjSktLAwB069bNxC0hImNhYCIiekwMTET1HwMTEVE1tM1gu3PnDv7++28AgL+/v1Heg4hMj4GJiOgxHDlyBADg6ekJJycnE7eGiIyFgYmISEeaLtLm6TiihoGBiYjoMTxOYOIsOSLzYZaBaenSpfDy8oKNjQ38/PxUd9jVZMyYMZBIJJWWDh06qMrExsZqLFNYWFgb3SEiM8YRJqKGwewC06ZNmzBt2jR8/PHHSE9PR8+ePTF48GBkZmZqLP/1118jOztbtWRlZcHR0RGvvPKKWjkHBwe1ctnZ2bCxsamNLhGRmbp+/TqysrIgkUjg5+dn6uYQkRGZXWBatGgRxo4di3HjxqFdu3ZYsmQJ3N3dsWzZMo3lZTIZ5HK5ajl8+DBu376NN998U62cRCJRKyeXy2ujO0Rkxo4dOwYAaNOmDZo0aVLjejhLjqjuM6vAVFxcjCNHjiAkJERtfUhICA4ePKhTHatWrUL//v3h4eGhtr6goAAeHh5o0aIFQkNDkZ6eXm09RUVFUCqVagsRNSzHjx8HAHTu3NnELSEiYzOrwHTz5k2UlZXBxcVFbb2LiwtycnK07p+dnY0//vgD48aNU1vftm1bxMbG4rfffsOGDRtgY2ODHj164Pz581XWNX/+fMhkMtXi7u5es04RkdmqCEydOnUycUuIyNjMKjBVePRZTUIInZ7fFBsbi6ZNm2LYsGFq6wMDA/H666+jc+fO6NmzJ3755Re0adMG3377bZV1RUVFQaFQqJasrKwa9YWIzNeJEycA1HyEibPkiMyHpakboA8nJydYWFhUGk3Kzc2tNOr0KCEEfvrpJ4SHh8Pa2rraso0aNUK3bt2qHWGSSqWQSqW6N56I6pWioiKcPXsWAE/JETUEZjXCZG1tDT8/PyQmJqqtT0xMRHBwcLX77t27F//88w/Gjh2r9X2EEDh27BhcXV0fq71EVH/99ddfKC0tRbNmzdCiRYvHqosXfRPVfWY1wgQAkZGRCA8Ph7+/P4KCgrBy5UpkZmZiwoQJAB6cKrt69SrWrl2rtt+qVasQEBAAHx+fSnXOnTsXgYGBaN26NZRKJb755hscO3YM33//fa30iYjqrqrCzMMXfOtySQARmTezC0xhYWHIy8tDTEwMsrOz4ePjg23btqlmvWVnZ1e6J5NCoUBcXBy+/vprjXXeuXMH48ePR05ODmQyGXx9fbFv3z50797d6P0hIvPEGXJEDYvZBSYAmDhxIiZOnKhxW2xsbKV1MpkM9+7dq7K+xYsXY/HixYZqHhHVUw8PJD3uBd9EZF7M6homIqK6QAhhkFsK8EwekflgYCIi0lN2djby8vJgYWGh9lxKIqq/GJiIiPR0+vRpAEDr1q0N8sxJzpIjqvsYmIiI9PTXX38BANq1a2filhBRbWFgIiLS05kzZwAA7du3N3FLiKi2MDAREempYoSJgYmo4WBgIiLSk6FOyXGWHJH5YGAiItLDjRs3kJeXB4lEAm9vb1M3h4hqCQMTEVE1Hp3BVjG65OXlBTs7O6O8BxHVPQxMRER64Aw5ooaJgYmISEcSCWfIETVUDExERHrgDDmihomBiYhID4Y8JcdZckTmg4GJiEhHd+7cQXZ2NgBew0TU0DAwERHp6OzZswCAJ598Eg4ODgarl7PkiOo+BiYiIh39888/AMD7LxE1QAxMREQ6unjxIgDg6aefNnFLiKi2MTAREemoYoSpdevWJm4JEdU2BiYiIh1duHABgOFGmDhLjsh8MDAREemoYoSJp+SIGh4GJiKiajw8g02pVAAAWrVqZbT3IKK6iYGJiEgPLVq0gK2trambQUS1jIGJiEhnghd8EzVQDExERHrg9UtEDRMDExGRHgwZmDhLjsh8MDAREenBGKfkeNE3Ud3HwEREpAeekiNqmBiYiIj0YOhbChCReWBgIiLSkZvbk7CzszN1M4jIBBiYiIh0xNElooaLgYmISEeGDkycJUdkPswyMC1duhReXl6wsbGBn58f9u/fX2XZPXv2QCKRVFrOnj2rVi4uLg7t27eHVCpF+/btkZCQYOxuEJEZeHgGm7Eu+OYsOaK6z+wC06ZNmzBt2jR8/PHHSE9PR8+ePTF48GBkZmZWu9+5c+eQnZ2tWh6eGpySkoKwsDCEh4fj+PHjCA8Px4gRI3Do0CFjd4eIzAhPyRE1XGYXmBYtWoSxY8di3LhxaNeuHZYsWQJ3d3csW7as2v2cnZ0hl8tVi4WFhWrbkiVLMGDAAERFRaFt27aIiopCv379sGTJEiP3hojMSevWvKUAUUNlVoGpuLgYR44cQUhIiNr6kJAQHDx4sNp9fX194erqin79+iE5OVltW0pKSqU6Bw4cWG2dRUVFUCqVagsR1T+3b99W/btly5YmbAkRmZJZBaabN2+irKwMLi4uautdXFyQk5OjcR9XV1esXLkScXFxiI+Ph7e3N/r164d9+/apyuTk5OhVJwDMnz8fMplMtbi7uz9Gz4iorsrIyFD9297e3oQtISJTsjR1A2pC8sjUEiFEpXUVvL294e3trXodFBSErKwsfPnll+jVq1eN6gSAqKgoREZGql4rlUqGJqJ66MKFCwC6GqVuzpIjMh9mNcLk5OQECwuLSiM/ubm5lUaIqhMYGIjz58+rXsvlcr3rlEqlcHBwUFuIqP65ePGi0d+Ds+SI6j6zCkzW1tbw8/NDYmKi2vrExEQEBwfrXE96ejpcXV1Vr4OCgirVuXPnTr3qJKL66cEIExE1dGZ3Si4yMhLh4eHw9/dHUFAQVq5ciczMTEyYMAHAg1NlV69exdq1awE8mAHn6emJDh06oLi4GOvWrUNcXBzi4uJUdb7zzjvo1asXFixYgBdeeAFbtmxBUlISDhw4YJI+ElHdwcBERIAZBqawsDDk5eUhJiYG2dnZ8PHxwbZt2+Dh4QEAyM7OVrsnU3FxMd577z1cvXoVtra26NChA7Zu3YrnnntOVSY4OBgbN27EzJkzMWvWLLRq1QqbNm1CQEBArfePiOoWBiYiAgCJEDx7bghKpRIymQwKhYLXMxHVE7du3ULz5l4AFACAoiLA2tpw9X/wAfDFF8D06cCXXxquXiLSna6/v83qGiYiotr08OQQY+AsOSLzwcBERFSFB4HJ+KmG4/xEdR8DExFRFYw9wkRE5oOBiYioCo8GJp5CI2q4GJiIiKrAESYiqsDARESkgRCCgYmIVBiYiIg0yMvLg0KhMOp7VJzi40XfRHUfAxMRkQYVo0tPPtnCxC0horqAgYmISIOKwNSyZUsTt4SI6gIGJiIiDSoCU6tWrUzcEiKqCxiYiIg0YGAioocxMBERafD3338D4Ck5InqAgYmI6BHl5eU4e/YsAKB16zZGex/OkiMyHwxMRESPyMzMxP3792FtbQ1PT09TN4eI6gAGJiKiR5w5cwYA0Lp1a1haWqrW89EoRA0XAxMR0SMqAlO7du1M3BIiqisYmIiIHsHARESPYmAiInoEAxMRPYqBiYjoIUKIWgtMnCVHZD4YmIiIHnLjxg3cunULEokE3t7epm4OEdURDExERA+pGF3y9PSEra2tiVtDRHUFAxMR0UN4/RIRacLARET0EAYmItKEgYmI6CEMTESkCQMTEdFDTp8+DQBo3749AOPOYOMsOSLzwcBERPQ/N2/exLVr1wAAPj4+Jm4NEdUlDExERP9z8uRJAEDLli3RpEmTStv5LDmihouBiYjof06cOAEA6Nixo4lbQkR1DQMTEdH/VIwwderUycQtIaK6hoGJiOh/KkaYajsw8aJvorqPgYmICEBZWRlOnToFoPYCE6+JIjIfZhmYli5dCi8vL9jY2MDPzw/79++vsmx8fDwGDBiAJ554Ag4ODggKCsKOHTvUysTGxkIikVRaCgsLjd0VIqojLly4gPv378PW1hatWrUydXOIqI4xu8C0adMmTJs2DR9//DHS09PRs2dPDB48GJmZmRrL79u3DwMGDMC2bdtw5MgR9O3bF0OHDkV6erpaOQcHB2RnZ6stNjY2tdElIqoDKq5f6tChAywsLEzcGiKqayxN3QB9LVq0CGPHjsW4ceMAAEuWLMGOHTuwbNkyzJ8/v1L5JUuWqL3+7LPPsGXLFvznP/+Br6+var1EIoFcLjdq24mo7jLV9UtEZB7MaoSpuLgYR44cQUhIiNr6kJAQHDx4UKc6ysvLkZ+fD0dHR7X1BQUF8PDwQIsWLRAaGlppBOpRRUVFUCqVagsRmS/eUoCIqmNWgenmzZsoKyuDi4uL2noXFxfk5OToVMdXX32Fu3fvYsSIEap1bdu2RWxsLH777Tds2LABNjY26NGjB86fP19lPfPnz4dMJlMt7u7uNesUEdUJR44cAQB07dpVbX1tzGDjLDmius+sAlMFySNTS4QQldZpsmHDBkRHR2PTpk1wdnZWrQ8MDMTrr7+Ozp07o2fPnvjll1/Qpk0bfPvtt1XWFRUVBYVCoVqysrJq3iEiMqnc3FxkZWVBIpGonao3Ns6SIzIfZnUNk5OTEywsLCqNJuXm5lYadXrUpk2bMHbsWPz73/9G//79qy3bqFEjdOvWrdoRJqlUCqlUqnvjiajOqhhdatOmjcZHohARmdUIk7W1Nfz8/JCYmKi2PjExEcHBwVXut2HDBowZMwb/+te/MGTIEK3vI4TAsWPH4Orq+thtJqK6ryIw+fv7V1uOI0JEDZdZjTABQGRkJMLDw+Hv74+goCCsXLkSmZmZmDBhAoAHp8quXr2KtWvXAngQlt544w18/fXXCAwMVI1O2draQiaTAQDmzp2LwMBAtG7dGkqlEt988w2OHTuG77//3jSdJKJadfjwYQCAn5+fiVtCRHWV2QWmsLAw5OXlISYmBtnZ2fDx8cG2bdvg4eEBAMjOzla7J9OKFStQWlqKSZMmYdKkSar1o0ePRmxsLADgzp07GD9+PHJyciCTyeDr64t9+/ahe/futdo3IjINXUeYiKjhkgjB+RmGoFQqIZPJoFAo4ODgYOrmEJGOrl+/DrlcDolEAqVSicaNG6ttz8sDnJwe/LusDGhkwAsZZs0CPvkEmDQJ+O47w9VLRLrT9fe3WV3DRERkaBWjS97e3pXCkrHxmigi88HAREQNWlpaGgBev0RE1WNgIqIGLSUlBcCD+7EREVWFgYmIGqzy8nJVYKru1iRERAxMRNRg/fXXX1AqlbCzs+NDd4moWgxMRNRgVTy0OyAgAJaWmu+ywmfJERHAwEREDVhFYDLV6TjOkiMyHwxMRNRgVQSmHj166FSeAYeo4WJgIqIG6caNG6oHbHOGHBFpw8BERA3Svn37AAAdOnRAs2bNTNwaIqrrGJiIqEHavXs3AODZZ581cUt40TeROWBgIqIGqS4FJiKq+xiYiKjBuXr1Ks6ePQuJRILevXubrB28iJzIfDAwEVGDk5ycDADo2rUrr18iIp0wMBFRg1NxOq5fv34mbgkRmQsGJiJqUIQQ2LlzJwBev0REumNgIqIGJT09HVevXoW9vb1Jr196GGfJEdV9DExE1KD8/vvvAIABAwbAxsZGa3mGGSICGJiIqIGpCEyhoaEmbglnyRGZEwYmImowcnJykJaWBgB47rnn9N6fAYeo4WJgIqIGIz4+HgDQrVs3uLq6mrg1RGROGJiIqMHYsGEDACAsLMzELSEic8PAREQNQmZmJg4cOACg7gUmXlhOVPcxMBFRg/DLL78AAHr27IkWLVqYuDUP8JooIvPBwERE9Z4QAj/99BMAYNSoUSZuDRGZIwYmIqr39u3bhzNnzsDe3p6BiYhqhIGJiOq9ZcuWAQBee+01ODg4mLg1RGSOGJiIqF7LyMjAr7/+CgCYMGGCiVtDROaKgYmI6rX58+ejrKwMISEh8PX11Xv/2pjBxllyRHUfAxMR1Vvnz59HbGwsAGD27NmmbYwGnCVHZD4YmIioXhJCYMqUKSgpKcGgQYPQo0cPUzeJiMyYWQampUuXwsvLCzY2NvDz88P+/furLb937174+fnBxsYGLVu2xPLlyyuViYuLQ/v27SGVStG+fXskJCQYq/lEVAt+/PFH7NixA9bW1vjmm29M3RwiMnN6BSalUokZM2YgJCQEw4cPx1dffYVLly4ZqWmabdq0CdOmTcPHH3+M9PR09OzZE4MHD0ZmZqbG8hkZGXjuuefQs2dPpKen46OPPsLUqVMRFxenKpOSkoKwsDCEh4fj+PHjCA8Px4gRI3Do0KHa6hYRGdDu3bsxZcoUAMC8efPQunVrE7eIiMydRAjdLzccMmQIDh06hIEDB8LCwgJHjx7FP//8g88++wyRkZHGbKdKQEAAunbtqpomDADt2rXDsGHDMH/+/ErlZ8yYgd9++w1nzpxRrZswYQKOHz+OlJQUAA8ek6BUKvHHH3+oygwaNAjNmjVTPXtKG6VSCZlMBoVCYdBpy+fP30Ju7j2N2/Q4dDrvZy7bdNlem+0xRf/rSltqu3/V7WtrW4yjR3fjk08+QVFREYYNG4a4uDg0alTzwfTcXMDFpeJ9a1yNRvPmAbNnA6++Cmj4+iKiRzg6Ak2aGLZOnX9/Cz04ODiIffv2qa07ePCgaNGihVi2bJk+VdVIUVGRsLCwEPHx8Wrrp06dKnr16qVxn549e4qpU6eqrYuPjxeWlpaiuLhYCCGEu7u7WLRokVqZRYsWiaeeeqrKthQWFgqFQqFasrKyBAChUChq0rUqtWu3Vzz4mubChYv2pUgA3QUAMXz4cFFYWPjYP4PXr///+g0tJsbUnxcXLua1rFhh+J9DhUIhdPn9balPCgsPD8fgwYMRHByMLl26oGvXrujcuTO+//57vPvuu0a/x8nNmzdRVlYGl4o/9/7HxcUFOTk5GvfJycnRWL60tBQ3b96Eq6trlWWqqhN4MFV57ty5NeyJ7iwtBYD7Rn8fqqtqNo2qIc6+EsIagDVathyOuXOn4LXXXoOkjn8QISHA0qXAnTumbgmRebCwMN176xWY/Pz80L9/fxw9ehTp6enYsGEDrl69Cmtra5SWlmLixIno2LEjOnbsiGeeecZYba70JSiEqPaLUVP5R9frW2dUVJTaaUilUgl3d3ftjdfTiRO9DV4nUX00bBiwZcuD0/Cvv27q1ugmIADIzjZ1K4hIF3oFpt27d+P+/ftYt24dbGxsAAB5eXnYsmULxo8fj4KCAqxcuRLnzp3DvXuar7t5HE5OTrCwsKg08pObm1tphKiCXC7XWN7S0hLNmzevtkxVdQKAVCqFVCqtSTeIiIjIzOh1JeSXX36Js2fPwsfHB59++im2bt2K3bt3Y/Xq1ejSpQvWrl2L9PR05OfnG6Wx1tbW8PPzQ2Jiotr6xMREBAcHa9wnKCioUvmdO3fC398fVlZW1Zapqk4iIiJqWPQKTC4uLkhLS8Po0aPxyy+/YOjQoQgLC8Pdu3exZs0aVTkLI55kjIyMxI8//oiffvoJZ86cwbvvvovMzEzV9VNRUVF44403VOUnTJiAy5cvIzIyEmfOnMFPP/2EVatW4b333lOVeeedd7Bz504sWLAAZ8+exYIFC5CUlIRp06YZrR9ERERkPvQ6JQcAtra2mDVrFmbNmoXCwkKUlpaicePGxmibRmFhYcjLy0NMTAyys7Ph4+ODbdu2wcPDAwCQnZ2tdk8mLy8vbNu2De+++y6+//57uLm54ZtvvsFLL72kKhMcHIyNGzdi5syZmDVrFlq1aoVNmzYhICCg1vpFRIbxv0sU62x9RGSe9LoPE1XNWPdhIiLdVFz0vXw58NZbhqv3+nVALn/wb35bEtU/uv7+NstHoxARPcrYdxCo43coICIjY2AiIiIi0oKBiYiIiEgLBiYiIiIiLRiYiKhe4YXZRGQMDExEREREWjAwEVG9wFlsRGRMDExEREREWjAwEREREWnBwEREVA1eRE5EAAMTEdUzDDhEZAwMTEREOuBF5UQNGwMTEdULDDREZEwMTERERERaMDARERERacHARERERKQFAxMR1SucJUdExsDARET1Ai/6JiJjYmAiIiIi0oKBiYiIiEgLBiYiIiIiLRiYiIiqwYvIiQhgYCKieoYBh4iMgYGJiOoFY8+S4yw8ooaNgYmIiIhICwYmIiIiIi0YmIiIiIi0YGAiIiIi0oKBiYjqFc6SIyJjYGAionqBs9iIyJgYmIiIiIi0MKvAdPv2bYSHh0Mmk0EmkyE8PBx37typsnxJSQlmzJiBjh07wt7eHm5ubnjjjTdw7do1tXJ9+vSBRCJRW0aOHGnk3hAREZG5MKvANGrUKBw7dgzbt2/H9u3bcezYMYSHh1dZ/t69ezh69ChmzZqFo0ePIj4+Hn///Teef/75SmUjIiKQnZ2tWlasWGHMrhCRmeA1UUQEAJamboCuzpw5g+3btyM1NRUBAQEAgB9++AFBQUE4d+4cvL29K+0jk8mQmJiotu7bb79F9+7dkZmZiaeeekq13s7ODnK53LidICIiIrNkNiNMKSkpkMlkqrAEAIGBgZDJZDh48KDO9SgUCkgkEjRt2lRt/fr16+Hk5IQOHTrgvffeQ35+frX1FBUVQalUqi1EZHrGGhHiReVEDZvZjDDl5OTA2dm50npnZ2fk5OToVEdhYSE+/PBDjBo1Cg4ODqr1r732Gry8vCCXy3Hq1ClERUXh+PHjlUanHjZ//nzMnTtX/44QkVEw0BCRMZl8hCk6OrrSBdePLocPHwYASDR8IwohNK5/VElJCUaOHIny8nIsXbpUbVtERAT69+8PHx8fjBw5Er/++iuSkpJw9OjRKuuLioqCQqFQLVlZWXr2nIiIiMyFyUeYJk+erHVGmqenJ06cOIHr169X2nbjxg24uLhUu39JSQlGjBiBjIwM7N69W210SZOuXbvCysoK58+fR9euXTWWkUqlkEql1dZDRERE9YPJA5OTkxOcnJy0lgsKCoJCocCff/6J7t27AwAOHToEhUKB4ODgKverCEvnz59HcnIymjdvrvW9Tp8+jZKSEri6uureESIiIqq3TH5KTlft2rXDoEGDEBERgdTUVKSmpiIiIgKhoaFqM+Tatm2LhIQEAEBpaSlefvllHD58GOvXr0dZWRlycnKQk5OD4uJiAMCFCxcQExODw4cP49KlS9i2bRteeeUV+Pr6okePHibpKxEREdUtZhOYgAcz2Tp27IiQkBCEhISgU6dO+Pnnn9XKnDt3DgqFAgBw5coV/Pbbb7hy5Qq6dOkCV1dX1VIxs87a2hq7du3CwIED4e3tjalTpyIkJARJSUmwsLCo9T4S0ePhfZOIyBhMfkpOH46Ojli3bl21ZcRD35aenp5qrzVxd3fH3r17DdI+IjIdzpIjImMyqxEmIiIiIlNgYCIiIiLSgoGJiKgavCaKiAAGJiIiIiKtGJiIqF6ouOibz5IjImNgYCIiIiLSgoGJiIiISAsGJiIiIiItGJiIiIiItGBgIiIiItKCgYmI6gVjz5IjooaNgYmIiIhICwYmIiIiIi0YmIiIqsFTfEQEMDARERERacXARESkAz4ahahhY2AionqBs+SIyJgYmIiIiIi0YGAiIiIi0oKBiYiIiEgLBiYiIiIiLRiYiIiIiLRgYCKieoGz5IjImBiYiIiIiLRgYCIiIiLSgoGJiKgaPMVHRAADExEREZFWDExEVC8Y+1lvfJYcUcPGwERE9QpPoRGRMTAwEREREWnBwERERESkhVkFptu3byM8PBwymQwymQzh4eG4c+dOtfuMGTMGEolEbQkMDFQrU1RUhClTpsDJyQn29vZ4/vnnceXKFSP2hIiIiMyJWQWmUaNG4dixY9i+fTu2b9+OY8eOITw8XOt+gwYNQnZ2tmrZtm2b2vZp06YhISEBGzduxIEDB1BQUIDQ0FCUlZUZqytERERkRixN3QBdnTlzBtu3b0dqaioCAgIAAD/88AOCgoJw7tw5eHt7V7mvVCqFXC7XuE2hUGDVqlX4+eef0b9/fwDAunXr4O7ujqSkJAwcONDwnSEig+MsNiIyJrMZYUpJSYFMJlOFJQAIDAyETCbDwYMHq913z549cHZ2Rps2bRAREYHc3FzVtiNHjqCkpAQhISGqdW5ubvDx8am23qKiIiiVSrWFiEyPs+SIyBjMJjDl5OTA2dm50npnZ2fk5ORUud/gwYOxfv167N69G1999RXS0tLw7LPPoqioSFWvtbU1mjVrprafi4tLtfXOnz9fdS2VTCaDu7t7DXtGREREdZ3JA1N0dHSli7IfXQ4fPgwAkGgYcxdCaFxfISwsDEOGDIGPjw+GDh2KP/74A3///Te2bt1abbu01RsVFQWFQqFasrKydOwxERERmRuTX8M0efJkjBw5stoynp6eOHHiBK5fv15p240bN+Di4qLz+7m6usLDwwPnz58HAMjlchQXF+P27dtqo0y5ubkIDg6ush6pVAqpVKrz+xKReeIpPiIC6kBgcnJygpOTk9ZyQUFBUCgU+PPPP9G9e3cAwKFDh6BQKKoNNo/Ky8tDVlYWXF1dAQB+fn6wsrJCYmIiRowYAQDIzs7GqVOnsHDhwhr0iIiIiOobk5+S01W7du0waNAgREREIDU1FampqYiIiEBoaKjaDLm2bdsiISEBAFBQUID33nsPKSkpuHTpEvbs2YOhQ4fCyckJL774IgBAJpNh7NixmD59Onbt2oX09HS8/vrr6Nixo2rWHBHVfXyWHBEZk8lHmPSxfv16TJ06VTWj7fnnn8d3332nVubcuXNQKBQAAAsLC5w8eRJr167FnTt34Orqir59+2LTpk1o0qSJap/FixfD0tISI0aMwP3799GvXz/ExsbCwsKi9jpHRAbBU2hEZAxmFZgcHR2xbt26asuIh74tbW1tsWPHDq312tjY4Ntvv8W333772G0kIiKi+sdsTskRERERmQoDExEREZEWDExEREREWjAwEVG9wFlsRGRMDExEVK9wlhwRGQMDExEREZEWDExERNXgiBURAQxMRERERFoxMBFRvcBHoxCRMTEwEREREWnBwERE9QqvOSIiY2BgIiIiItKCgYmIiIhICwYmIiIiIi0YmIioXuAsNiIyJgYmIiIiIi0YmIioXuEsOSIyBgYmIiIiIi0YmIiIqsERKyICGJiIiIiItGJgIqJ6gc+SIyJjYmAiIiIi0oKBiYjqFV5zRETGwMBEREREpAUDExEREZEWDExEREREWjAwEVG9wFlsRGRMDExEREREWjAwEVG9wllyRGQMDExERNVgACMigIGJiIiISCuzCky3b99GeHg4ZDIZZDIZwsPDcefOnWr3kUgkGpcvvvhCVaZPnz6Vto8cOdLIvSEiQ+KjUYjImCxN3QB9jBo1CleuXMH27dsBAOPHj0d4eDj+85//VLlPdna22us//vgDY8eOxUsvvaS2PiIiAjExMarXtra2Bmw5ERERmTOzCUxnzpzB9u3bkZqaioCAAADADz/8gKCgIJw7dw7e3t4a95PL5Wqvt2zZgr59+6Jly5Zq6+3s7CqVNRYhBMrKylBaWlor70fUEDRtCnh4APb2QGGh4eotK3tQr52dYeutCSsrK1hYWJi2EUQNlNkEppSUFMhkMlVYAoDAwEDIZDIcPHiwysD0sOvXr2Pr1q1Ys2ZNpW3r16/HunXr4OLigsGDB2POnDlo0qRJlXUVFRWhqKhI9VqpVGp9fyEE7ty5gxs3bqCsrExreSLS3QsvAP36PQhOGRmGq7e0FFi+/MEpOUPWW1NNmzaFXC6HhOcIiWqV2QSmnJwcODs7V1rv7OyMnJwcnepYs2YNmjRpguHDh6utf+211+Dl5QW5XI5Tp04hKioKx48fR2JiYpV1zZ8/H3PnztW7D3fu3IGDgwMcHBxgaWnJLz0iA7G2Bm7fBpydHyyGUlwMFBU9CExeXoarV19CCNy7dw+5ubkAAFdXV9M1hqgBMnlgio6O1ho80tLSAEBjuBBC6Bw6fvrpJ7z22muwsbFRWx8REaH6t4+PD1q3bg1/f38cPXoUXbt21VhXVFQUIiMjVa+VSiXc3d2rfO+ysjIoFAo88cQTcHJy0qm9RKS7ijNVVlbAIz/ij6Xi60UiMWy9NVFxbWVubi6cnZ15eo6oFpk8ME2ePFnrjDRPT0+cOHEC169fr7Ttxo0bcHFx0fo++/fvx7lz57Bp0yatZbt27QorKyucP3++ysAklUohlUq11lWhpKQEQgjY29vrvA8R0aPs7OwAPPhOYWAiqj0mD0xOTk46jbgEBQVBoVDgzz//RPfu3QEAhw4dgkKhQHBwsNb9V61aBT8/P3Tu3Flr2dOnT6OkpMQoQ948BUdEj4PfIUSmYTb3YWrXrh0GDRqEiIgIpKamIjU1FREREQgNDVW74Ltt27ZISEhQ21epVOLf//43xo0bV6neCxcuICYmBocPH8alS5ewbds2vPLKK/D19UWPHj2M3i8iIiKq+8wmMAEPZrJ17NgRISEhCAkJQadOnfDzzz+rlTl37hwUCoXauo0bN0IIgVdffbVSndbW1ti1axcGDhwIb29vTJ06FSEhIUhKSuJwN5EZ4qNMiMgYTH5KTh+Ojo5Yt25dtWWEhm/L8ePHY/z48RrLu7u7Y+/evQZpHxEREdVPZjXCRHVXbGys2qNlLC0t0aJFC7z55pu4evWq3vUtXboUsbGxhm+oBgcPHkR0dLTGx+z06dMHffr0qVG9j7Mv8GCyw5gxY1Sv9+zZA4lEgj179uhVj7E/y7Vr1+KJJ55Afn6+0d7DXOl6zHbt2oXGjRvX6GeFiGoHAxMZ1OrVq5GSkoLExERERERgw4YN6NmzJ+7evatXPbUdmObOnasxMC1duhRLly6tlXZo07VrV6SkpFQ5c7Mqxvws7927h48++ggzZsyo9kavtcGcr4Xu168funfvjo8++sjUTSGiKjAwkUH5+PggMDAQffv2xZw5c/DBBx8gIyMDmzdvNnXTaqR9+/Zo3769qZsBAHBwcEBgYCAcHByM9h4lJSV6PbJnzZo1yMvL0zihgvQzadIkrF+/HllZWaZuChFpwMBkYkII3L17t84smq4BexyBgYEAgMuXLwMACgsLERUVBS8vL1hbW+PJJ5/EpEmT1EZ3PD09cfr0aezdu1d1is/T01O1XalU4r333lOrY9q0aZVGsSQSCSZPnoyff/4Z7dq1g52dHTp37ozff/9dVSY6Ohrvv/8+AMDLy0v1fhWnUDSdVps7dy4CAgLg6OgIBwcHdO3aFatWrarxZ1dSUoIPPvgAcrkcdnZ2eOaZZ/Dnn39WKqfp9M7FixcxcuRIuLm5QSqVwsXFBf369cOxY8e0fpYV9f3888+YPn06nnzySUilUvzzzz86t33ZsmUYOnQomjZtqrZel88eAMaMGaN2bCtER0dXmj5fUefq1avh7e0NW1tb+Pv7IzU1FUIIrFjxBV54wQutWzfGs88+q1c/HlZeXo5PPvlE9R4uLk3x6qud8K9/fa1W7sCBA+jXrx+aNGkCOzs7BAcHY+vWrTV6TwAYOnQoGjdujB9++KHGdRCR8ZjVRd/10b1799C4cWNTN0OloKDAoDfXrPil9cQTT0AIgWHDhmHXrl2IiopCz549ceLECcyZMwcpKSlISUmBVCpFQkICXn75ZchkMtXpsIqbhN67dw+9e/fGlStX8NFHH6FTp044ffo0Zs+ejZMnTyIpKUntF+3WrVuRlpaGmJgYNG7cGAsXLsSLL76Ic+fOoWXLlhg3bhxu3bqFb7/9FvHx8ap7b1U3qnTp0iW89dZbeOqppwAAqampmDJlCq5evYrZs2fr/RlFRERg7dq1eO+99zBgwACcOnUKw4cP1+maoOeeew5lZWVYuHAhnnrqKdy8eRMHDx5UBdDqPssKUVFRCAoKwvLly9GoUSONjyDS5MqVKzh58iTefvttjdu1ffY18fvvvyM9PR2ff/45JBIJZsyYgSFDhmD06NE4deoi3n//O1hZKTBvXiReeuklHDt2TO/7Fi1cuBDR0dGYOXMmevXqhXv3SpCcfBYFBXdUZfbu3YsBAwagU6dOWLVqFaRSKZYuXYqhQ4diw4YNCAsL07tv1tbWqtAVExOj9/5EZGSCDEKhUAgAQqFQaNx+//598ddff4n79++rrS8oKBAA6sxSUFBQo/6vXr1aABCpqamipKRE5Ofni99//1088cQTokmTJiInJ0ds375dABALFy5U23fTpk0CgFi5cqVqXYcOHUTv3r0rvc/8+fNFo0aNRFpamtr6X3/9VQAQ27ZtU60DIFxcXIRSqVSty8nJEY0aNRLz589Xrfviiy8EAJGRkVHp/Xr37q2xHRXKyspESUmJiImJEc2bNxfl5eU67yuEEGfOnBEAxLvvvqu2fv369QKAGD16tGpdcnKyACCSk5OFEELcvHlTABBLliyp9j2q+iwr6uvVq1e1+1el4rilpqZW2qbrZz969Gjh4eFRaf85c+aIR7+eAAi5XK72f3Tz5s0CgOjSpYvIyCgXaWlCXL0qxJIlSwQAceLECb37FRoaKrp06aJ6XVgoRFqaEEeO/P8ygYGBwtnZWeTn56vWlZaWCh8fH9GiRQvV/4NHj5k2H3/8sWjUqFG1P4dVfZcQUc1o+/1dgafkTMzOzg4FBQV1Zql47EJNBQYGwsrKCk2aNEFoaCjkcjn++OMPuLi4YPfu3QCgNvMLAF555RXY29tj165dWuv//fff4ePjgy5duqC0tFS1DBw4UONspL59+6pdjOzi4gJnZ2fVKcKa2L17N/r37w+ZTAYLCwtYWVlh9uzZyMvLUz0YVVfJyckAHjwA+mEjRoyApWX1A8COjo5o1aoVvvjiCyxatAjp6ekoLy/XrzMAXnrpJb33AYBr164BQJUjUsb47Pv27as2AtquXTsAwODBg9VGkirW1+S9unfvjuPHj2PixInYsWMHlEql2va7d+/i0KFDePnll9VGhy0sLBAeHo4rV67g3Llzer8v8OCzLC8v1/mB4kRUe3hKzsQkEkm9er7c2rVr0a5dO1haWsLFxUXt8TJ5eXmwtLTEE088obaPRCKBXC5HXl6e1vqvX7+Of/75B1ZWVhq337x5U+118+bNK5WRSqW4f/++Lt2p5M8//0RISAj69OmDH374AS1atIC1tTU2b96MTz/9VO96K/osl8vV1ltaWmps+8MkEgl27dqFmJgYLFy4ENOnT4ejoyNee+01fPrppzrPWqvpI4Aq+vrow6wrGPqzBx6ExIdZW1tXu76wsFDv94iKioK9vT3WrVuH5cuXw8LCAl269MI77yxA167+uH37NoQQGj83Nzc3ANDp/7ImFZ/l43xGRGQcDExkUO3atYO/v7/Gbc2bN0dpaSlu3LihFpqEEMjJyUG3bt201u/k5ARbW1v89NNPVW43po0bN8LKygq///67WlCo6SzAilCRk5ODJ598UrW+tLRUp1+6Hh4eWLVqFQDg77//xi+//ILo6GgUFxdj+fLlOrWhps8mq/isb926VePQZWNjg6KiokrrHw2+tcnS0hKRkZGIjIzEnTt38McfSfjww48wefJADB+ehWbNmqFRo0bIzs6utG/FqFtN/x/eunXrsfYnIuPhKTmqNf369QOASndrj4uLw927d1XbgapHIkJDQ3HhwgU0b94c/v7+lRZNM660qbgIWpe/6ituyvnwY3Pu379f6RE9uqqYgbd+/Xq19b/88ote0/sBoE2bNpg5cyY6duyIo0ePqtY/7qhOVdq2bQvgwfMYa8rT0xO5ubm4fv26al1xcTF27Njx2O0zhKZNm2L48Jfx8suToFDcwqVLl2Bvb4+AgADEx8erfa7l5eVYt24dWrRogTZt2tTo/S5evIjmzZvDxcXFUF0gIgNhYKJaM2DAAAwcOBAzZszA3LlzkZSUhEWLFuHNN9+Er68vwsPDVWU7duyI48ePY9OmTUhLS8PJkycBANOmTYO3tzd69eqFRYsWISkpCTt37sSPP/6IESNG4NChQ3q3q2PHjgCAr7/+GikpKTh8+HCVM9SGDBmCgoICjBo1ComJidi4cSN69uxZaeaZrtq1a4fXX38dS5YswYwZM5CYmIjFixfj/fff13q/pRMnTqBXr1749ttvsX37duzevRszZ87EiRMnMGDAALX+afosqzJ27FhYWlpqvf4nICAAtra2SE1N1b3DjwgLC4OFhQVGjhyJbdu2IT4+HiEhISgrK6txnVWJiYmBpaWl1kchDR06FFFRUYiLi8O+ffuwfv3P2LhxCVxdPdC6dWsAwPz585GXl4e+ffvi119/xW+//YbnnnsOp06dwpdfflntqN3atWthaWmJtWvXVtqWmpqK3r1713jUj4iMh6fkqNZIJBJs3rwZ0dHRWL16NT799FM4OTkhPDwcn332mVromDt3LrKzsxEREYH8/Hx4eHio/rrfv38/Pv/8c6xcuRIZGRmwtbXFU089hf79+9dohKlPnz6IiorCmjVr8MMPP6C8vBzJyckaH2vy7LPP4qeffsKCBQswdOhQPPnkk4iIiICzszPGjh1bo89l1apVcHFxQWxsLL755ht06dIFcXFxGDlyZLX7yeVytGrVCkuXLkVWVhYkEglatmyJr776ClOmTFGVq+qzrEpZWRnKysq03lfK2toaL7/8MrZs2YLPPvtMrz5X8PLywpYtW/DRRx/h5ZdfhqurKyIjI3Hjxg3MnTu3RnVWpby8XKd+9e3bF3Fxcfjxxx+hVCrh4iJHt24DEBExS3XtXO/evbF7927MmTMHY8aMQXl5OTp37ozffvsNoaGhOrXj0Qv0L1y4gJMnTyI6Ovqx+klExiER2r49SCdKpRIymQwKhULjyEBhYSEyMjLg5eVV5UWyRObm8OHD6NatG1JTUxEQEGDStly+DNy4Abi5PVgMpbAQOHUKaNQI0POpNHqZNWsW1q5diwsXLlQ7Q5LfJUSGpe33dwWekiOiGvP398eIESMwb948UzfFrN25cwfff/89PvvsM623kyAi02BgIqLH8tVXX6Fbt2463ZmcNMvIyEBUVBRGjRpl6qYQURX4pwwRPZYWLVpgzpw5pm6GWfP19YWvr6+pm0FE1eAIExEREZEWDExEREREWjAwEREREWnBwERERESkBQMTERERkRYMTERERERaMDARERERacHARAZ14sQJvPnmm6rHNjRu3Bhdu3bFwoULcevWLVU5T09Prc/cqi+io6MrPUzV09MTY8aM0auegwcPIjo6Gnfu3DFc44iISCe8cSUZzA8//ICJEyfC29sb77//Ptq3b4+SkhIcPnwYy5cvR0pKChISEkzdzDohISGh2mcWaXLw4EHMnTsXY8aMQdOmTY3TMKrSI5mXiBoYBiYyiJSUFLz99tsYMGAANm/eDKlUqto2YMAATJ8+Hdu3b6/VNpWVlaG0tFStLXVFbdzV+f79+7C1tTX6+xARNQQ8JWdiQgB379adRYia9eOzzz6DRCLBypUrNQYUa2trPP/885XWb9++HV27doWtrS3atm2Ln376SW37jRs3MHHiRLRv3x6NGzeGs7Mznn32Wezfv1+t3KVLlyCRSLBw4UJ88skn8PLyglQqRXJyMgBgy5Yt6NSpE6RSKVq2bImvv/5a46kyIQSWLl2KLl26wNbWFs2aNcPLL7+Mixcv6vQ5bN26FV26dIFUKoWXlxe+/PJLjeUePSVXXl6OTz75BN7e3rC1tUXTpk3RqVMnfP311wAenNZ7//33AQBeXl6QSCSQSCTYs2ePqr7Q0FDEx8fD19cXNjY2mDt3rk5tJiIi7TjCZGL37gGNG5u6Ff9fQQFgb6/fPmVlZdi9ezf8/Pzg7u6u837Hjx/H9OnT8eGHH8LFxQU//vgjxo4di6effhq9evUCANV1T3PmzIFcLkdBQQESEhLQp08f7Nq1C3369FGr85tvvkGbNm3w5ZdfwsHBAa1bt8b27dsxfPhw9OrVC5s2bUJpaSm+/PJLXL9+vVKb3nrrLcTGxmLq1KlYsGABbt26hZiYGAQHB+P48eNwcXGpsj+7du3CCy+8gKCgIGzcuBFlZWVYuHChxvd51MKFCxEdHY2ZM2eiV69eKCkpwdmzZ1XXK40bNw63bt3Ct99+i/j4eLi6ugIA2rdvr6rj6NGjOHPmDGbOnAkvLy/Y63sg64mahn4iomoJMgiFQiEACIVCoXH7/fv3xV9//SXu37+vtr6gQIgHX/F1Yyko0L/vOTk5AoAYOXKkzvt4eHgIGxsbcfnyZbXPyNHRUbz11ltV7ldaWipKSkpEv379xIsvvqhan5GRIQCIVq1aieLiYrV9unXrJtzd3UVRUZFqXX5+vmjevLl4+EcgJSVFABBfffWV2v5ZWVnC1tZWfPDBB9X2KSAgQLi5uakdY6VSKRwdHcWjP2oeHh5i9OjRqtehoaGiS5cu1db/xRdfCAAiIyOj0jYPDw9hYWEhzp07V20d9dmlS0KkpQlx5Yph671//0G9R48att6aquq7hIhqRtvv7wo8JWdidnYPRnXqymJnV3t979KlC5566inVaxsbG7Rp0waXL19WK7d8+XJ07doVNjY2sLS0hJWVFXbt2oUzZ85UqvP555+HlZWV6vXdu3dx+PBhDBs2DNbW1qr1jRs3xtChQ9X2/f333yGRSPD666+jtLRUtcjlcnTu3Fl1+kuTu3fvIi0tDcOHD4eNjY1qfZMmTSq9jybdu3fH8ePHMXHiROzYsQNKpVLrPo/q1KkT2rRpo/d+RESkHU/JmZhEov8psLrGyckJdnZ2yMjI0Gu/5s2bV1onlUpx//591etFixZh+vTpmDBhAubNmwcnJydYWFhg1qxZGgNTxamqCrdv34YQQuOptEfXXb9+vcqyANCyZcsq+3L79m2Ul5dDLpdX2qZp3aOioqJgb2+PdevWYfny5bCwsECvXr2wYMEC+Pv7a90fqNz3hoaz2IjImMxqhOnTTz9FcHAw7OzsdJ5WLYRAdHQ03NzcYGtriz59+uD06dNqZYqKijBlyhQ4OTnB3t4ezz//PK5cuWKEHtRPFhYW6NevH44cOWLwz23dunXo06cPli1bhiFDhiAgIAD+/v7Iz8/XWP7Ri7ibNWsGiUSi8TqinJwctddOTk6QSCQ4cOAA0tLSKi2bN2+usp0V7/NonZreRxNLS0tERkbi6NGjuHXrFjZs2ICsrCwMHDgQ9+7d07o/ULnvRERkOGYVmIqLi/HKK6/g7bff1nmfhQsXYtGiRfjuu++QlpYGuVyOAQMGqP3CnTZtGhISErBx40YcOHAABQUFCA0NRVlZmTG6US9FRUVBCIGIiAgUFxdX2l5SUoL//Oc/etcrkUgqzbo7ceIEUlJSdNrf3t4e/v7+2Lx5s1q7CgoK8Pvvv6uVDQ0NhRACV69ehb+/f6WlY8eO1b5P9+7dER8fj8LCQtX6/Px8vfvdtGlTvPzyy5g0aRJu3bqFS5cuAYDqc3h4BI6IiGqHWZ2Sq5gmHRsbq1N5IQSWLFmCjz/+GMOHDwcArFmzBi4uLvjXv/6Ft956CwqFAqtWrcLPP/+M/v37A3gwquHu7o6kpCQMHDjQKH2pb4KCgrBs2TJMnDgRfn5+ePvtt9GhQweUlJQgPT0dK1euhI+Pj07X8zwsNDQU8+bNw5w5c9C7d2+cO3cOMTEx8PLyQmlpqU51xMTEYMiQIRg4cCDeeecdlJWV4YsvvkDjxo3V7j7eo0cPjB8/Hm+++SYOHz6MXr16wd7eHtnZ2Thw4AA6duxYbVifN28eBg0apLrvVFlZGRYsWAB7e3u199Fk6NCh8PHxgb+/P5544glcvnwZS5YsgYeHB1q3bg0AqsD29ddfY/To0bCysoK3tzeaNGmisc69e/eiX79+mD17NmbPnq3TZ1UflJUBRUWGq6+kxHB1EZH5MqvApK+MjAzk5OQgJCREtU4qlaJ37944ePAg3nrrLRw5cgQlJSVqZdzc3ODj44ODBw9WGZiKiopQ9NC3ck0u0q1vIiIi0L17dyxevBgLFixATk4OrKys0KZNG4waNQqTJ0/Wu86PP/4Y9+7dw6pVq7Bw4UK0b98ey5cvR0JCQrUXYT9s0KBBiIuLw+zZsxEWFga5XI6JEyfi2rVr+Pnnn9XKrlixAoGBgVixYgWWLl2K8vJyuLm5oUePHujevXu171Nx086ZM2eqvc/9+/e13hOpb9++iIuLw48//gilUqkaCZ01a5bqIvY+ffogKioKa9aswQ8//IDy8nIkJydXurVCBSEEysrKUF5ertPnVF/k5j5YiIgMSSKE+d21JDY2FtOmTdP6TK2DBw+iR48euHr1Ktzc3FTrx48fj8uXL2PHjh3417/+hTfffFMt/ABASEgIvLy8sGLFCo11R0dHa/wlqFAoND7yorCwEBkZGapnrJFplZSUoEuXLnjyySexc+dOUzeHDEChADIyHowwGYOTE+DhYZy69cHvEiLDUiqVkMlkVf7+rmDyEaaqgsfD0tLSdJ4ppImmuzlru0BWW5moqChERkaqXiuVSr1u2ki1a+zYsRgwYABcXV2Rk5OD5cuX48yZM6o7aZP5k8mALl1M3Qoiqq9MHpgmT56MkSNHVlvG09OzRnVXTOfOyclRm3Kdm5urmjoul8tRXFyM27dvo1mzZmplgoODq6xbKpXWyWeUkWb5+fl47733cOPGDVhZWaFr167Ytm2b6ro1IiKi6pg8MDk5OcHJyckodXt5eUEulyMxMVH1sNPi4mLs3bsXCxYsAAD4+fnBysoKiYmJGDFiBAAgOzsbp06dwsKFC43SLqp9v/zyi6mbQEREZszkgUkfmZmZuHXrFjIzM1FWVoZjx44BAJ5++mk0/t8D2dq2bYv58+fjxRdfhEQiwbRp0/DZZ5+hdevWaN26NT777DPY2dlh1KhRAACZTIaxY8di+vTpaN68ORwdHfHee++hY8eOHH0gIiIiAGYWmGbPno01a9aoXleMGj08U+jcuXNQKBSqMh988AHu37+PiRMn4vbt2wgICMDOnTvVpmIvXrwYlpaWGDFiBO7fv49+/fohNjYWFhYWBu+DGV5jT0R1CL9DiEzDLGfJ1UXarrIvKyvD+fPnjXoKkojqv7y8POTm5qJNmzZG+aOOqKExm1lyDYWFhQVkMhlu3LiBoqIiODg4wNLSko+zICKdCCFw79495ObmomnTpgxLRLWMgakWyeVy2NraIjc3lze6JKIaadq0qU4PdCYiw2JgqkUSiQRNmzaFTCZDWVmZzo/2ICICACsrK44sEZkIA5MJSCQSWFpawtKSHz8REZE5aGTqBhARERHVdQxMRERERFowMBERERFpwcBEREREpAUDExEREZEWnKZlIBU3TOf9lYiIiMxHxe9tbQ8+YWAykPz8fACAu7u7iVtCRERE+srPz4dMJqtyO58lZyDl5eW4du0amjRpYtDHnSiVSri7uyMrK6vaZ9yYs/rex/reP6D+95H9M3/1vY/sX80JIZCfnw83Nzc0alT1lUocYTKQRo0aoUWLFkar38HBoV7+EDysvvexvvcPqP99ZP/MX33vI/tXM9WNLFXgRd9EREREWjAwEREREWnBwFTHSaVSzJkzB1Kp1NRNMZr63sf63j+g/veR/TN/9b2P7J/x8aJvIiIiIi04wkRERESkBQMTERERkRYMTERERERaMDARERERacHAVAd8+umnCA4Ohp2dHZo2baqxTGZmJoYOHQp7e3s4OTlh6tSpKC4urrbeoqIiTJkyBU5OTrC3t8fzzz+PK1euGKEHutuzZw8kEonGJS0trcr9xowZU6l8YGBgLbZcP56enpXa++GHH1a7jxAC0dHRcHNzg62tLfr06YPTp0/XUot1d+nSJYwdOxZeXl6wtbVFq1atMGfOHK3/H+v6MVy6dCm8vLxgY2MDPz8/7N+/v9rye/fuhZ+fH2xsbNCyZUssX768llqqn/nz56Nbt25o0qQJnJ2dMWzYMJw7d67afar6OT179mwttVo/0dHRldoql8ur3cdcjh+g+ftEIpFg0qRJGsvX9eO3b98+DB06FG5ubpBIJNi8ebPa9pp+F8bFxaF9+/aQSqVo3749EhISDNpuBqY6oLi4GK+88grefvttjdvLysowZMgQ3L17FwcOHMDGjRsRFxeH6dOnV1vvtGnTkJCQgI0bN+LAgQMoKChAaGgoysrKjNENnQQHByM7O1ttGTduHDw9PeHv71/tvoMGDVLbb9u2bbXU6pqJiYlRa+/MmTOrLb9w4UIsWrQI3333HdLS0iCXyzFgwADVcwrrirNnz6K8vBwrVqzA6dOnsXjxYixfvhwfffSR1n3r6jHctGkTpk2bho8//hjp6eno2bMnBg8ejMzMTI3lMzIy8Nxzz6Fnz55IT0/HRx99hKlTpyIuLq6WW67d3r17MWnSJKSmpiIxMRGlpaUICQnB3bt3te577tw5tePVunXrWmhxzXTo0EGtrSdPnqyyrDkdPwBIS0tT61tiYiIA4JVXXql2v7p6/O7evYvOnTvju+++07i9Jt+FKSkpCAsLQ3h4OI4fP47w8HCMGDEChw4dMlzDBdUZq1evFjKZrNL6bdu2iUaNGomrV6+q1m3YsEFIpVKhUCg01nXnzh1hZWUlNm7cqFp39epV0ahRI7F9+3aDt72miouLhbOzs4iJiam23OjRo8ULL7xQO40yAA8PD7F48WKdy5eXlwu5XC4+//xz1brCwkIhk8nE8uXLjdBCw1q4cKHw8vKqtkxdPobdu3cXEyZMUFvXtm1b8eGHH2os/8EHH4i2bduqrXvrrbdEYGCg0dpoKLm5uQKA2Lt3b5VlkpOTBQBx+/bt2mvYY5gzZ47o3LmzzuXN+fgJIcQ777wjWrVqJcrLyzVuN6fjB0AkJCSoXtf0u3DEiBFi0KBBausGDhwoRo4cabC2coTJDKSkpMDHxwdubm6qdQMHDkRRURGOHDmicZ8jR46gpKQEISEhqnVubm7w8fHBwYMHjd5mXf3222+4efMmxowZo7Xsnj174OzsjDZt2iAiIgK5ubnGb+BjWLBgAZo3b44uXbrg008/rfaUVUZGBnJyctSOl1QqRe/evevU8aqKQqGAo6Oj1nJ18RgWFxfjyJEjap89AISEhFT52aekpFQqP3DgQBw+fBglJSVGa6shKBQKANDpePn6+sLV1RX9+vVDcnKysZv2WM6fPw83Nzd4eXlh5MiRuHjxYpVlzfn4FRcXY926dfi///s/rQ96N6fjV6Gm34VVHVNDfn8yMJmBnJwcuLi4qK1r1qwZrK2tkZOTU+U+1tbWaNasmdp6FxeXKvcxhVWrVmHgwIFwd3evttzgwYOxfv167N69G1999RXS0tLw7LPPoqioqJZaqp933nkHGzduRHJyMiZPnowlS5Zg4sSJVZavOCaPHue6drw0uXDhAr799ltMmDCh2nJ19RjevHkTZWVlen32mn4mXVxcUFpaips3bxqtrY9LCIHIyEg888wz8PHxqbKcq6srVq5cibi4OMTHx8Pb2xv9+vXDvn37arG1ugsICMDatWuxY8cO/PDDD8jJyUFwcDDy8vI0ljfX4wcAmzdvxp07d6r9I9Pcjt/DavpdWNUxNeT3p6XBaiI10dHRmDt3brVl0tLStF63U0HTXxJCCK1/YRhiH13UpL9XrlzBjh078Msvv2itPywsTPVvHx8f+Pv7w8PDA1u3bsXw4cNr3nA96NPHd999V7WuU6dOaNasGV5++WXVqFNVHj02xjpemtTkGF67dg2DBg3CK6+8gnHjxlW7b104htXR97PXVF7T+rpk8uTJOHHiBA4cOFBtOW9vb3h7e6teBwUFISsrC19++SV69epl7GbqbfDgwap/d+zYEUFBQWjVqhXWrFmDyMhIjfuY4/EDHvyROXjwYLUzDo8yt+OnSU2+C439/cnAZCSTJ0/GyJEjqy3j6empU11yubzShWu3b99GSUlJpUT98D7FxcW4ffu22ihTbm4ugoODdXpffdSkv6tXr0bz5s3x/PPP6/1+rq6u8PDwwPnz5/Xet6Ye55hWzAb7559/NAamihk9OTk5cHV1Va3Pzc2t8hgbmr79u3btGvr27YugoCCsXLlS7/czxTHUxMnJCRYWFpX+Eq3us5fL5RrLW1paVhuITWnKlCn47bffsG/fPrRo0ULv/QMDA7Fu3TojtMzw7O3t0bFjxyr/b5nj8QOAy5cvIykpCfHx8Xrvay7Hr6bfhVUdU0N+fzIwGYmTkxOcnJwMUldQUBA+/fRTZGdnq/4D7dy5E1KpFH5+fhr38fPzg5WVFRITEzFixAgAQHZ2Nk6dOoWFCxcapF0P07e/QgisXr0ab7zxBqysrPR+v7y8PGRlZan9QBnb4xzT9PR0AKiyvV5eXpDL5UhMTISvry+AB9cq7N27FwsWLKhZg/WkT/+uXr2Kvn37ws/PD6tXr0ajRvqf3TfFMdTE2toafn5+SExMxIsvvqhan5iYiBdeeEHjPkFBQfjPf/6jtm7nzp3w9/ev0f9nYxJCYMqUKUhISMCePXvg5eVVo3rS09NNfqx0VVRUhDNnzqBnz54at5vT8XvY6tWr4ezsjCFDhui9r7kcv5p+FwYFBSExMVFtdH/nzp2GHSAw2OXjVGOXL18W6enpYu7cuaJx48YiPT1dpKeni/z8fCGEEKWlpcLHx0f069dPHD16VCQlJYkWLVqIyZMnq+q4cuWK8Pb2FocOHVKtmzBhgmjRooVISkoSR48eFc8++6zo3LmzKC0trfU+PiopKUkAEH/99ZfG7d7e3iI+Pl4IIUR+fr6YPn26OHjwoMjIyBDJyckiKChIPPnkk0KpVNZms3Vy8OBBsWjRIpGeni4uXrwoNm3aJNzc3MTzzz+vVu7hPgohxOeffy5kMpmIj48XJ0+eFK+++qpwdXWtc328evWqePrpp8Wzzz4rrly5IrKzs1XLw8zpGG7cuFFYWVmJVatWib/++ktMmzZN2Nvbi0uXLgkhhPjwww9FeHi4qvzFixeFnZ2dePfdd8Vff/0lVq1aJaysrMSvv/5qqi5U6e233xYymUzs2bNH7Vjdu3dPVebR/i1evFgkJCSIv//+W5w6dUp8+OGHAoCIi4szRRe0mj59utizZ4+4ePGiSE1NFaGhoaJJkyb14vhVKCsrE0899ZSYMWNGpW3mdvzy8/NVv+cAqL4vL1++LITQ7bswPDxcbRbrf//7X2FhYSE+//xzcebMGfH5558LS0tLkZqaarB2MzDVAaNHjxYAKi3JycmqMpcvXxZDhgwRtra2wtHRUUyePFkUFhaqtmdkZFTa5/79+2Ly5MnC0dFR2NraitDQUJGZmVmLPavaq6++KoKDg6vcDkCsXr1aCCHEvXv3REhIiHjiiSeElZWVeOqpp8To0aPrTF8edeTIEREQECBkMpmwsbER3t7eYs6cOeLu3btq5R7uoxAPptPOmTNHyOVyIZVKRa9evcTJkydrufXarV69WuP/10f//jK3Y/j9998LDw8PYW1tLbp27ao27X706NGid+/eauX37NkjfH19hbW1tfD09BTLli2r5Rbrpqpj9fD/vUf7t2DBAtGqVSthY2MjmjVrJp555hmxdevW2m+8jsLCwoSrq6uwsrISbm5uYvjw4eL06dOq7eZ8/Crs2LFDABDnzp2rtM3cjl/FbQ8eXUaPHi2E0O27sHfv3qryFf79738Lb29vYWVlJdq2bWvwgCgR4n9XuhERERGRRrytABEREZEWDExEREREWjAwEREREWnBwERERESkBQMTERERkRYMTERERERaMDARERERacHARERERKQFAxMRkR7GjBmDYcOG6bVPdHQ0JBIJJBIJlixZ8ljvHxsbq6pr2rRpj1UXEemOgYmI6pyahBJDu3TpEiQSCY4dO2aQ+jp06IDs7GyMHz/+seoJCwtDdnY2goKCDNIuItKNpakbQETUEFhaWkIulz92Pba2trC1tYW1tbUBWkVEuuIIExGZnb/++gvPPfccGjduDBcXF4SHh+PmzZuq7X369MHUqVPxwQcfwNHREXK5HNHR0Wp1nD17Fs888wxsbGzQvn17JCUlQSKRYPPmzQAALy8vAICvry8kEgn69Omjtv+XX34JV1dXNG/eHJMmTUJJSYne/ZBIJFixYgVCQ0NhZ2eHdu3aISUlBf/88w/69OkDe3t7BAUF4cKFC3rXTUSGxcBERGYlOzsbvXv3RpcuXXD48GFs374d169fx4gRI9TKrVmzBvb29jh06BAWLlyImJgYJCYmAgDKy8sxbNgw2NnZ4dChQ1i5ciU+/vhjtf3//PNPAEBSUhKys7MRHx+v2pacnIwLFy4gOTkZa9asQWxsLGJjY2vUn3nz5uGNN97AsWPH0LZtW4waNQpvvfUWoqKicPjwYQDA5MmTa1Q3ERkOT8kRkVlZtmwZunbtis8++0y17qeffoK7uzv+/vtvtGnTBgDQqVMnzJkzBwDQunVrfPfdd9i1axcGDBiAnTt34sKFC9izZ4/qNNmnn36KAQMGqOp84oknAADNmzevdCqtWbNm+O6772BhYYG2bdtiyJAh2LVrFyIiIvTuz5tvvqkKezNmzEBQUBBmzZqFgQMHAgDeeecdvPnmm3rXS0SGxREmIjIrR44cQXJyMho3bqxa2rZtCwBqp646deqktp+rqytyc3MBAOfOnYO7u7taEOrevbvObejQoQMsLCw01q2vh9vp4uICAOjYsaPausLCQiiVyhrVT0SGwREmIjIr5eXlGDp0KBYsWFBpm6urq+rfVlZWatskEgnKy8sBAEIISCSSGrehurofp66KNmlaV9P6icgwGJiIyKx07doVcXFx8PT0hKVlzb7C2rZti8zMTFy/fl01qpOWlqZWpmIWWllZ2eM1mIjqBZ6SI6I6SaFQ4NixY2pLZmYmJk2ahFu3buHVV1/Fn3/+iYsXL2Lnzp34v//7P53DzYABA9CqVSuMHj0aJ06cwH//+1/VRd8VIzrOzs6wtbVVXVSuUCiM1lciqvsYmIioTtqzZw98fX3VltmzZ8PNzQ3//e9/UVZWhoEDB8LHxwfvvPMOZDIZGjXS7SvNwsICmzdvRkFBAbp164Zx48Zh5syZAAAbGxsAD+6b9M0332DFihVwc3PDCy+8YLS+ElHdJxFCCFM3gojI1P773//imWeewT///INWrVoZtO7o6Ghs3rzZYHcNBx7ca6pLly6P/agVItINR5iIqEFKSEhAYmIiLl26hKSkJIwfPx49evQweFiqcPLkSTRu3BhLly59rHrWr1+Pxo0bY//+/QZqGRHpgiNMRNQgrV27FvPmzUNWVhacnJzQv39/fPXVV2jevLnB3+vWrVu4desWgAf3d5LJZDWuKz8/H9evXwcANG3aFE5OTgZpIxFVj4GJiIiISAuekiMiIiLSgoGJiIiISAsGJiIiIiItGJiIiIiItGBgIiIiItKCgYmIiIhICwYmIiIiIi0YmIiIiIi0+H8e9M4CXLUmWgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualization\n",
    "plt.figure()\n",
    "plt.plot(x, c.value, 'k-', label='Potential distr. (num. sol.)')\n",
    "plt.plot(x, rho(x), 'b-', label='Charge distr.')\n",
    "plt.xlabel('Length [m]')\n",
    "plt.ylabel(r'$\\varphi$')\n",
    "plt.legend(fontsize=12, loc='best');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9dc1e12",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}