Pellet.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'fipy.solvers.petsc.linearGMRESSolver.LinearGMRESSolver'>\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import fipy as fp\n",
    "import scipy.constants as sc\n",
    "\n",
    "import time\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "#plt.rcParams[\"animation.html\"] = \"jshtml\"\n",
    "#from matplotlib.animation import FuncAnimation\n",
    "#from IPython.display import HTML\n",
    "\n",
    "def update_progress(process,progress):\n",
    "    import sys\n",
    "    barLength = 100 # Modify this to change the length of the progress bar\n",
    "    status = \"\"\n",
    "    if isinstance(progress, int):\n",
    "        progress = float(progress)\n",
    "    if not isinstance(progress, float):\n",
    "        progress = 0\n",
    "        status = \"error: progress var must be float\\r\\n\"\n",
    "    if progress < 0:\n",
    "        progress = 0\n",
    "        status = \"Halt...\\r\\n\"\n",
    "    if progress >= 1:\n",
    "        progress = 1\n",
    "        status = \"Done...\\r\\n\"\n",
    "    block = int(round(barLength*progress))\n",
    "    text = \"\\r\"+process+\": [{0}] {1:.0f}% {2}\".format( \"#\"*block + \"-\"*(barLength-block), progress*100, status)\n",
    "    sys.stdout.write(text)\n",
    "    sys.stdout.flush()\n",
    "    \n",
    "print(fp.DefaultSolver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Te = 4000 #eV\n",
    "n0 = 1e8 #line integrated density\n",
    "n_hot = 5e19 #ambient density [m^-3]\n",
    "M = 2. * sc.m_p #atomic mass\n",
    "\n",
    "LL = 12\n",
    "#Braginskii Eq. 2.5e\n",
    "tee = 12. * np.pi**1.5 /2**0.5 * sc.m_e**0.5 * sc.epsilon_0**2 * (sc.e * Te)**1.5 / (sc.e**4 * n_hot * LL)\n",
    "SQ  = 3. * n0 * Te * 1.6e-19 / tee \n",
    "tau = 1./ (3. * n0) * SQ # Dnestrovskiy Eq. 65\n",
    "\n",
    "#Analytical solution\n",
    "n_ana   = lambda x,t: np.sqrt(3.*M/(8.*t**3*np.pi*tau))*np.exp(-3.*M*x**2/(8.*t**3*tau))\n",
    "#T = tau * t\n",
    "#v  = 1.5 * x / t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our mesh will span $[-L_x/2;L_x/2]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lx = 400.; dt = 0.5e-6; steps = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#base = 30; N =200\n",
    "#nonUniform_x = Lx/(base-1.)*(np.logspace(0.,1.,N+1,base=base,endpoint=True)-1.)\n",
    "#mesh  = fp.Grid1D(dx=np.diff(nonUniform_x))\n",
    "\n",
    "mesh = fp.Grid1D(nx=1000, Lx=Lx) - 0 * Lx/2.; \n",
    "x = mesh.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "''Shallow water'' equations for expanding plasma take the form\n",
    "\\begin{equation}\n",
    "\\frac{\\partial n}{\\partial t} + \\frac{1}{\\partial x} (nv) = 0,\n",
    "\\end{equation}\n",
    "\\begin{equation}\n",
    "\\frac{\\partial v}{\\partial t} + v\\frac{\\partial v}{\\partial x} = - \\frac{T(t)}{nM}\\frac{\\partial n}{\\partial x}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = fp.CellVariable(mesh=mesh, hasOld = True)\n",
    "v = fp.CellVariable(mesh=mesh, hasOld = True)\n",
    "v_face = fp.FaceVariable(mesh=mesh, rank=0, value=0.)\n",
    "\n",
    "T = fp.Variable(value = 0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cont_eq = (fp.TransientTerm(var=n) + fp.ConvectionTerm(coeff=v.faceValue*[[1.0]], var=n) == 0)\n",
    "#moti_eq = (fp.TransientTerm(var=v) + fp.ConvectionTerm(coeff=v.faceValue*[[0.5]], var=v) + T / M * n.grad[0]/n == fp.DiffusionTerm(coeff=1.e6,var=v)) # T * (fp.numerix.log(n.faceValue)).divergence #\n",
    "\n",
    "#cont_eq = (fp.TransientTerm(var=n) + fp.ConvectionTerm(coeff=v_face*[[1.0]], var=n) == 0)\n",
    "#moti_eq = (fp.TransientTerm(var=v) + fp.ConvectionTerm(coeff=v_face*[[0.5]], var=v) + T / M * n.grad[0]/n == fp.DiffusionTerm(coeff=1.e6,var=v)) # T * (fp.numerix.log(n.faceValue)).divergence #\n",
    "\n",
    "cont_eq = (fp.TransientTerm(var=n) + fp.ConvectionTerm(coeff=v_face*[[1.]], var=n) == 0)\n",
    "moti_eq = (fp.TransientTerm(var=v) + fp.ConvectionTerm(coeff=v_face*[[1.]], var=v) - \\\n",
    "                    fp.ImplicitSourceTerm(coeff=v_face.divergence, var=v) + T / M * n.grad[0]/n == fp.DiffusionTerm(coeff=1.e6,var=v))\n",
    "\n",
    "eqn = cont_eq & moti_eq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rough boundary conditions for $v$ and $n$ on the right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.faceGrad.constrain(0.  , where=mesh.facesLeft)  #correct\n",
    "n.faceGrad.constrain(0.  , where=mesh.facesRight) #arbitrary \n",
    "v.constrain(0.           , where=mesh.facesLeft)  #correct\n",
    "v.faceGrad.constrain(0.  , where=mesh.facesRight) #arbitrary \n",
    "\n",
    "#Initial conditions == analytical around finite density\n",
    "t = 10.e-6\n",
    "T.setValue(tau*t)\n",
    "v.setValue(1.5*x/t, where= x < 10.)\n",
    "v.setValue(0., where= x >= 10.)\n",
    "n.setValue(n0*n_ana(x,t)+1e-15*n0/300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check convergence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6188891535743.795\n",
      "2.945613471800911e+18\n",
      "4.479535579484638e+23\n",
      "4.493617561019341e+22\n",
      "1.0926114997987734e+22\n",
      "2.6101651991604704e+21\n",
      "6.3433172575606e+20\n",
      "1.557947349073216e+20\n",
      "3.822709259039078e+19\n",
      "9.3989248208487e+18\n",
      "2.305900214123529e+18\n",
      "5.612070827616478e+17\n",
      "1.3471000795458498e+17\n",
      "3.1673003724301332e+16\n",
      "7404961891972778.0\n",
      "1804493661598282.2\n",
      "454029500619669.8\n",
      "115181525629723.89\n",
      "29244019658586.12\n",
      "6891873137491.468\n",
      "1216826843027.6611\n",
      "487031370988.6832\n",
      "450870278727.06464\n",
      "342068069441.8376\n",
      "114776606072.20105\n",
      "173491166168.29446\n",
      "68185862514.93431\n",
      "96390732769.00186\n",
      "32970059233.913746\n",
      "59075440802.921196\n",
      "18544272537.058613\n",
      "33974079007.941\n",
      "10975715093.802519\n",
      "20339508984.982826\n",
      "6326305756.408315\n",
      "11714800818.758322\n",
      "3756028003.838876\n",
      "6864207513.677132\n",
      "2197163294.2535133\n",
      "3971136310.2591634\n",
      "1299263551.323517\n",
      "2311657162.4683156\n",
      "764252159.0824708\n",
      "1340123474.1037154\n",
      "451159815.2799031\n",
      "778529258.4939299\n",
      "265826209.45168963\n",
      "451643786.9495594\n",
      "156795516.97767588\n",
      "262188442.43688115\n",
      "92413852.91909589\n",
      "152128439.18404034\n",
      "54480674.581204176\n",
      "88286633.83028062\n",
      "32105678.22230987\n",
      "51226018.56567497\n",
      "18919139.500352442\n",
      "29723702.029231265\n",
      "11145826.19650641\n",
      "17245327.97415274\n",
      "6565448.723071438\n",
      "10005353.837815005\n",
      "3866599.4143496775\n",
      "5804504.937045605\n",
      "2276797.0880771717\n",
      "3367301.5594638246\n",
      "1340421.4151246992\n",
      "1953330.764352148\n",
      "789018.6486195371\n",
      "1133054.2051285587\n",
      "464365.8338244098\n",
      "657210.0289319128\n",
      "273252.2067188667\n",
      "381186.596479136\n",
      "160767.36991860933\n",
      "221080.30585875548\n",
      "94572.53143895179\n",
      "128215.89256904215\n",
      "55624.49049494311\n",
      "74355.41711918484\n",
      "32711.66019552794\n",
      "43118.37155696443\n",
      "19234.275890100424\n",
      "25002.925457755016\n",
      "11308.035657283195\n",
      "14497.661562166046\n",
      "6647.1852939364735\n",
      "8405.887164936057\n",
      "3906.8697477871697\n",
      "4873.573095664956\n",
      "2295.9464348367537\n",
      "2825.462636804642\n",
      "1349.079028990978\n",
      "1637.9839880760887\n",
      "792.6052464065117\n",
      "949.5275548883478\n",
      "465.6090169483002\n",
      "550.4058619656377\n",
      "273.48410447684034\n",
      "319.03357818093633\n"
     ]
    }
   ],
   "source": [
    "n.updateOld()\n",
    "v.updateOld()\n",
    "\n",
    "res = 1e100\n",
    "for i in range(100):\n",
    "    v_face[:] = v.arithmeticFaceValue\n",
    "    v_face[..., mesh.facesLeft.value] = 0.\n",
    "    res = eqn.sweep(dt=dt)\n",
    "    print(res)\n",
    "t += dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb2534d1f40>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAEDCAYAAAA/TmJVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAuj0lEQVR4nO3de5xVdb3/8ddn7xlAHZEIArkophxDOCZCXkJ0NFIgT6YnO3iOdlOxftXJ2yntYsewMjPLjiVaolmmZqCpUWnKDqwwAbkMoog4JoJyGS4O95n9+f2x10ybYS57Zvbea62938/HYx7O3nvtvd6sNXz98J3vZy1zd0REREREyk0i7AAiIiIiImFQISwiIiIiZUmFsIiIiIiUJRXCIiIiIlKWVAiLiIiISFlSISwiIiIiZSnUQtjMZpjZejOryWHbH5jZ4uBrpZltKUJEERHJonFbREqJhXkdYTM7FagH7nX3UZ143xeA0e7+6YKFExGR/WjcFpFSEuqMsLvPBeqynzOzI83sD2a20Mzmmdl7WnnrBcD9RQkpIiLNNG6LSCmpCDtAK+4EPuPuL5vZicBPgDOaXjSzw4EjgKdDyiciIvvSuC0isRSpQtjMqoD3Aw+ZWdPTPVtsNgX4jbs3FjObiIjsT+O2iMRZpAphMks1trj7ce1sMwX4XHHiiIhIBzRui0hsReryae6+DXjVzM4HsIz3Nr1uZkcD7wD+FlJEERHJonFbROIs7Mun3U9mcDzazNaY2cXAfwEXm9kSYDlwTtZbLgAe8DAvdSEiUsY0botIKQn18mkiIiIiImGJ1NIIEREREZFiCa1Zrl+/fj5s2LBOv2/79u0cdNBB+Q/UBVHJEpUcEJ0sUckByhLlHND1LAsXLtzo7v0LECmSNGaXXg5QlijngOhkiUoOKMCY7e6hfI0ZM8a7Ys6cOV16XyFEJUtUcrhHJ0tUcrgrS2uiksO961mABR7S+BnGl8bs/IlKDndlaU1UcrhHJ0tUcrjnf8zW0ggRERERKUsqhEVERESkLKkQFhEREZGyFLU7y4lIjDU2NlJXV8fevXvb3a53796sXbu2SKna11GWyspK+vbtSzKZLGKqeMjlfOtci0iUqRAWkbypq6ujV69e9OvXDzNrc7s9e/YwaNCgIiZrW3tZ3J36+nrq6uro379sLhCRs1zOt861iESZlkaISN7s3buXqqqqdovgODEzqqqqOpzhLleldL51rkXKkwphEcmrUiiKspXanyffSun4lNKfRURyE6tCeNaiNaRe2x12DBEREREptmW/of/6Z/L6kbFaIzz8iU9w9O6twFlhRxERERGRYti1DWZfDUsfZGDf48G/Cnn6DU6sCmG3BEkaw44hIiIiIsXw+t9h5iWw9XWovpaa9Ps4LY/LmGK1NCJtFVR4Q9gxRCTizj33XL72ta8xfvx4Bg4cyJ/+9KewI0mBLFy4kNNPP735cU1NDSeffHKIiUQkLxobIPVdmDERcPjUH6D6GjyR38sbxmpGOJ2o1IywSExc/9hyXli7rdXXdu3aRa9e6zr9mccM6s03/m1kh9vV1NQwbtw45s2bx6xZs7jvvvuYMGFCp/cnuWvrfHf1XENu53vEiBGsXLmy+fF1113HtGnTurQ/EYmIza/BrKnw+nw49j9g8veg1yEF2VWsZoQ9UaFCWETatWPHDrZu3coVV1wBQENDA3369AHg61//eojJpBAOPPBAevXqxZYtW1i0aBGbN29mwoQJOtcicbX0IZh+Cqx/Ac77GZx3Z8GKYIjZjLCbCmGRuGhvJq+2tpZhw4YVZL/Lly9nzJgxzXcHW7p0KaNGjeLNN9+koUFLqwqlrfNdyHPd5JhjjuHFF19k2rRp3HDDDTrXInGU1RDH0BMzBfA7hhV8t/GaEbYKKlyFsIi0raamhuOOO6758dKlSzn22GN5/vnn93leSsfIkSO5++67cXfGjRuncy0SN6//PTMLvOwhqL4WPjm7KEUwxKwQTicrqED/yheRti1btmyfIqimpoZRo0axePFiFUclauTIkfz0pz9tXhuscy0SE/s0xNHcEEeyeAsWYrU0AqugQksjRKQdt9xyyz6PV69eDcCqVasYPnx4GJGkwC666CIuuuii5sc61yIxsF9D3M3Qq3fRY8SqEPZEJUnSYccQkRi66667wo4gRaJzLRJxSx+C312Z+f68n8Gx54cWJVaFMImklkaIiIiIxNE+DXEnBQ1xh4caKVaFsCcqtTRCREREJG7+8SzMuhS2roHqr8D4q4q6Frgt4SfojEQFFZbG02ksEas+PxEREZHy09gA826GP98EhwyBT/8Bhp4QdqpmsSqEPVkJQLpxL8lEz5DTiIiIiEib9mmImxLcIa74DXHtiVUhTCJTCO/du4dkpQphERERkUjKboj797vgXz8abp42xKwQzsRt3Ls35CAiIiIisp9dW2H2/0SqIa49sSqELZG5ZWpjo64cISIiIhIp/3gWZl0CW9+IVENcezrsODOzoWY2x8xWmNlyM/tiK9tUm9lWM1scfF1XiLBNhbCndeUIERERkUhobIDUjXD3JMAyDXHVX458EQy5zQg3AFe5+yIzOxhYaGZPuvsLLbab5+5n5z9iFsvU7ZoRFpGuuP7666mrq6NPnz5cf/31rFmzhgsuuICPfexjzJ8/nwcffDDsiJJH2ef70ksv1bkWKYTNr2Uui/b6s5FtiGtPhzPC7r7O3RcF378NrAAGFzpYq4JLpnnaQ9m9iMTXG2+8wd69e+nTpw/z588HYMmSJZx55plcccUVVFREf+ZCctfyfOtcixTA0l/D9FNg/YpMQ9x5d8SqCIZOrhE2s2HAaODZVl4+2cyWAGuBq919eSvvnwpMBRgwYACpVKpTYbet3wDAc889S6/eqzv13kKor6/v9J+hlHNAdLJEJQeUV5bevXuzZ8+eDrfbs2cPtbW1Bctx2WWXMXz4cJ599llqa2v5wQ9+wKOPPsp1111HXV0dK1asoLa2llQqxRlnnEFtbS3bt29vM1NdXR0rV64sWF7puoULF3L11VczZ84cAGpqarj00ksZMWIEt956Kxs2bOD1119nyZIlnHrqqQCYWZiRReJv11b43dWw7Ndw2Mlw7h2RbohrT86FsJlVATOBy919W4uXFwGHu3u9mU0GHgGGt/wMd78TuBNg7NixXl1d3amwz25aDBvh+NGjedfQozr13kJIpVJ09s9QyjkgOlmikgPKK8vatWsZNGhQh9vV1tYybNiwguVYvXo1H/zgB7ntttuYNWsWjz32GCeeeCIPPfQQmzZtYvz48QwbNoxNmzZx9NFHU1VVxVFHHdVmph49enD88ccXLK903YgRI/b5R8p1113HtGnTWLJkCTfffDObNm1i9OjRPP/880yZMoWNGzcycODAEBOLxFx2Q9zpX4VTrozFWuC25JTczCrJFMH3ufuslq9nF8buPtvMfmJm/dx9Y/6i0rxGOK1mOZHo+/018OayVl8auGsX9OrV+c8c+K8w6cZ2N9mxYwdbt27liiuuAKChoYE+ffpw1VVX7bftjBkzqK2tpV+/ftx8882dzxNzZlYLvA00Ag3uPrbLH9bG+e7yuYaczveBBx5Ir1692LJlC6tXr2bz5s1MmDCBCRMm7LdtOZ9rkW5rbIC534O5N8EhQyN3h7iu6rAQtszvkO4CVrj7LW1sMxB4y93dzE4gs/Z4U16TQvNtldONKoRFpHXLly9nzJgxJJOZq8wsXbqUUaNGhZwq0k7P+6RFkR1zzDG8+OKLTJs2jRtuuCHsOCIlp9fOt+CeybFtiGtPLjPC44CLgGVmtjh47ivAYQDuPh34KPBZM2sAdgJT3D3/HW2W+R8bns77R4tInrUzk/dmAZdG1NTUcNxxxzU/Xrp0Keeccw4AX//615k2bVpB9lv22jjfhTzXTUaOHMndd9+NuzNu3DhA51okb5b+mrELLoeKykjfIa6rOiyE3f0ZoN3OAne/DbgtX6Ha0jwjrKURItKGZcuWceKJJzY/rqmpYdSoUbz55ps0NOjSiy048ISZOXBH0MfRrKMG51yaIwvdGAnwrne9i5tuuonf/va31NbWsmHDBjZt2rTffnPJUozGyHJqoO2MqGSJSg4IN0uyYTv/snI6A9bPZWvV0awcdTW7N/WDkI9Nvo9JvFY3BzPCKoRFpC233LLvCq7VqzNXmEmlUvvMFAsA49x9rZm9C3jSzF5097lNL3bU4JxLc2ShGyMBrrzySq688srmxytWrOC0007bb7+5ZClGY2Q5NdB2RlSyRCUHhJjlH/Nh1heaG+JqGsdw2hkfKH6OVuT7mHR4HeEoabqzHCqERaSTFi9erEK4BXdfG/x3PfAwEP/OF3SuRbqssQHmfCfrDnF/hNO+hDfVXyUoVjPC1nzVCK0RFpHOWbVqFcOH73dVx7JlZgcBCXd/O/j+TOCbIcfKC51rkS7YXAszL4U1f4f3XgCTbiqZhrj2xKoQ1uXTRKSr7rrrrrAjRM0A4OHg5hIVwK/c/Q/hRsoPnWuRTlryIPzuqkydVYINce2JVSFsweWQXIWwiEi3uPtq4L1h5xCREO3amimAlz2UuUPceXdCn8PCTlVU8SqEgxlhtDRCREREpOv+MR9mXRo0xH0Nxl8JJbwWuC2xKoSbTpC7ZoRFREREOi37DnF9Dss0xA19X9ipQhOrQvifd5bTjLCIiIhIp5RpQ1x7YlYIa42wSNS5O0EDVkkoxE0yS0kpnW+daylpZdwQ1554FcLBGmEtjRCJpsrKSurr66mqqiqJ4sjdqa+vp7KyMuwokVRK51vnWkqWGuLaFa9CuHlGWEsjRKKob9++1NXV8fbbb7e7XV1dHT169ChSqvZ1lKWyspK+ffsWMVF85HK+da5FQvSP+ZmlENvKuyGuPTEthDUjLBJFyWSS/v37d7jdypUrC34b21xFKUvc5HK+o3R8o5RFpKAaGzLNcHO/p4a4DsSrEG5aGqEZYREREZH91b0Ks6YGDXH/CZO+W/YNce2JVyHcNJ2vGWERERGRfWU3xH10Boz697ATRV7MCuHg8mmuGWERERERoEVD3PvhvDvUEJejmBXCWiMsIiIi0kwNcd0Sy0IYXT5NREREylnLhriLn4AhY8NOFTuxKoQTpsuniYiISJmrexVmXQprnss0xE2+CXoeHHaqWIpVIWzJpqtGaEZYREREyow7LH0Qfne1GuLyJF6FcPPSCM0Ii4iISBnZtRUevxJqfqOGuDyKZyGsGWEREREpF6/9LXNt4G1vwBlfg1PUEJcvsSqEE8Hl01wzwiIiIlLiLN0IT38L5t2shrgCiVUh/M/Lp6kQFhERkRJW9yrHLb4Wtr2khrgCilUhnNDSCBERESllWQ1xBzWm1RBXYLEqhE1LI0RERKRU7dySuUNc0BD33KBPc7KK4IJKhB2gMxK6s5yIiIiUotf+BtPHw/KH4YyvwycfZ3ev/mGnKnkdFsJmNtTM5pjZCjNbbmZfbGUbM7MfmdkqM1tqZscXIqwlMxPYphlhERERKQWNDZmGuHsmZ64EcfGTcOrVuipEkeSyNKIBuMrdF5nZwcBCM3vS3V/I2mYSMDz4OhG4PfhvXiUSBmhphIiIiJSAutWZy6KteQ6O+y+Y9F01xBVZh4Wwu68D1gXfv21mK4DBQHYhfA5wr7s7MN/M+pjZocF780bXERYREZHYc4clD8Dsq8GS8NG7YdR5YacqS51qljOzYcBo4NkWLw0GXs96vCZ4Lq+FcCKRiasZYREREYmlnVvgd1dCzUw4fBycewf0GRp2qrKVcyFsZlXATOByd9/W8uVW3uKtfMZUYCrAgAEDSKVSuScFduyoZzKwYf1bnX5vIdTX1ytHC1HJEpUcoCxRzgHRylJsZpYEFgBvuPvZYecRKXmv/TW4Q9zaTEPcKVdoLXDIciqEzaySTBF8n7vPamWTNUD2P2eGAGtbbuTudwJ3AowdO9arq6s7FXbz5jr4O/R/Z1/e18n3FkIqlaKzf4ZSzgHRyRKVHKAsUc4B0coSgi8CK4DeYQcRKWmNe+HP34V534c+h2ca4oaMCTuVkNtVIwy4C1jh7re0sdmjwMeDq0ecBGzN9/pggEQy+FeTlkaIiHSLmQ0BPgT8LOwsIiWtbjXMmAhzvwfvvQA+M09FcITkMiM8DrgIWGZmi4PnvgIcBuDu04HZwGRgFbAD+FTekwKJ4IYauJrlRES66YfAl4A2W9S7u5wNorX0JCpZopIDlKWgOdwZ8NYchr98B25JVh7zP2zocwr8bWHxs3RTVHJA/rPkctWIZ2h9DXD2Ng58Ll+h2tLULKcZYRGRrjOzs4H17r7QzKrb2q67y9kgWktPopIlKjlAWQqWY+cWePwKeHFWc0PcyC40xJXUMcmTfGeJ1S2WE82XT1MhLCLSDeOAD5vZZKAX0NvMfunuF4acSyT+mhri3l6nhrgYiNctlpOZuLp8mohI17n7te4+xN2HAVOAp1UEi3RT4154+ga450OQqIBPP6E7xMVAzGaEE6TdtDRCREREoqNuNcy8FN5YAMddCJNu1B3iYiJehbAZaQxTs5yISF64ewpIhRxDJJ7cYcn9MPt/MjO/598DI88NO5V0QswKYdhDQjPCIiIiEq6mhrjls3SHuBiLVSFsZjgGac0Ii4iISEjUEFcyYlUIAzSSyPwqQkRERKSYsu8Q945hmYY43Rwj1mJXCKe1NEJERESKbb+GuO9Cz6qwU0k3xa4QdjXLiYiISLGoIa6kxa4QbtSMsIiIiBTDPg1xp8B5d8AhQ8JOJXkUu0LY0XWERUREpMBq/wIPX5ZpiPvAdTDucjXElaDYFcKaERYREZFCsXQDPDUNnrkl0xB38RMwWA1xpSp2hXBmjbAKYREREcmzutWMfv5aeHsljL4QJqohrtTFrhDWVSNEREQkr7Ia4g5IuxriykgsC2FDhbCIiIjkwc7N8PiVzQ1xCw79JCerCC4bibADdJaa5URERCQvav8Ct58CKx7NNMR94lF29+ofdiopohjOCOs6wiIiItINjXshdaMa4iSGhbDpFssiIiLSRXWrYeYl8MZCNcRJ/AphV7OciIiIdJY7LP4V/P5LukOcNItdIZzW5dNERESkM3ZuDu4Q9zAMGw/nTtcd4gSIYSGs6wiLiIhIzmr/ArOmQv2b8IFvwLgv6g5x0ix2hXCaBOjyaSIiItKexr2Q+g7MuwX6vlsNcdKq2BXCmhEWERGRdm16BWZdqoY46VDsCmHdUENERERatU9DXAWc/3MY+ZGwU0mExa4Q1oywiIiI7GfnZnjscnjhETXESc5iVwjrOsIiIiKyj9pnYNZlaoiTTotdIawZYREREQHUECfdluhoAzObYWbrzaymjderzWyrmS0Ovq7Lf8x/cq0RFhERkU2vwF1nwrzvZxriLpurIlg6LZcZ4XuA24B729lmnrufnZdEHUhjVGhGWEREpDy5w+L7YPaXIFkJH7sXjjkn7FQSUx0Wwu4+18yGFSFLTjQjLCIiUqbUECd5lq81wieb2RJgLXC1uy9vbSMzmwpMBRgwYACpVKrTO+qN0bB3T5fem2/19fXK0UJUskQlByhLlHNAtLIUk5n1AuYCPcn8v+A37v6NcFOJtCO7IW7C/8L7/1sNcdJt+SiEFwGHu3u9mU0GHgGGt7ahu98J3AkwduxYr66u7vzO5iaorEjSlffmWyqVUo4WopIlKjlAWaKcA6KVpch2A2cEY3cl8IyZ/d7d54cdTGQf+zXEPQmDjw87lZSIbhfC7r4t6/vZZvYTM+vn7hu7+9mtSWMktEZYRKRb3N2B+uBhZfCla1NKtGx6BWZeAmsXweiLYOKNukOc5FW3C2EzGwi85e5uZieQuRLFpm4na4PWCIuI5IeZJYGFwFHAj9392Ravd3s5W5SWnkQlS1RyQISzuDPwzacY/vJPSScqeGnkl9l4yPvhbwuKmyNkUckSlRyQ/ywdFsJmdj9QDfQzszXAN8jMHODu04GPAp81swZgJzAlmGkoCDfDdEMNEZFuc/dG4Dgz6wM8bGaj3L0m6/VuL2eL0tKTqGSJSg6IaJamhriXHoFh40meewejDhlc/BwREJUsUckB+c+Sy1UjLujg9dvIXF6tKDQjLCKSX+6+xcxSwESg1WvGixRF7TMwayrUv6WGOCkK3VlORKQMmVl/YG9QBB8ATAC+G3IsKVeNezli9S8gNVMNcVJUsSuE05ZQISwi0n2HAj8P1gkngF+7++MhZ5JyFDTEHb52ERz/cTjrO2qIk6KJXSHsGKbGZhGRbnH3pcDosHNIGXOH538Jv/8yJCupGfllRn34K2GnkjITw0I4ocuniYiIxNl+d4i7g43Pvxx2KilDsSuEMVOznIiISFy12RCnQliKL3aFcFpXjRAREYmfxr0w59vwzA/gnUfCJX+CQVqdI+GKXSHsGAldR1hERCQ+Nr0CMy+Gtc9nGuIm3gg9Dgo7lUgcC+EECc0Ii4iIRF+Lhjg+di8cc07YqUSaxa8QNi2NEBERibydm+GxL8ILv21uiKOId4gTyUUMC2FdPk1ERCTSXp0HD18WNMRdD+//gu4QJ5EUv0KYBKY1wiIiItHTuBfmfAue+aEa4iQWYlcIg5GkMewQIiIikk0NcRJDsSuE05bQ0ggREZGoyG6Iq+gBH/sFHPPhsFOJ5CR2hTC6xbKIiEg07KiDxy/PNMQdcWqmIa73oLBTieQsdoWwmy6fJiIiErrmhrj18MFvwslfgEQi7FQinRK7QhhMhbCIiEhYGvZA6ttZDXFPqiFOYit2hXBmRlhLI0RERIpu4yqYdUnQEPcJmPgdNcRJrMWuEAbDXDPCIiIiReMOz/8iaIjrCf/xSxjxb2GnEum22BXCmhEWEREpoh11mTvErXhUDXFScmJXCGuNsIiISJG8OhdmXQbbN6ghTkpS7AphzQiLiIgUWMOezB3i/nJrpiHuAjXESWmKZyFsnlmvZBZ2HBERkdKycVXmDnHrFqshTkpe7AphyBS/7mnMkiFnERERKRFqiJMyFL9C2DJrk9KNjSQTKoRFRES6TQ1xUqZiVwh7UAg3NjaSrAw5jIiISNypIU7KWOwK4aZ1wel0Y8hBREREYqxhD+9+5eeQejhoiPsTDDou7FQiRRW/QpimQrgh5BwiIiIxFTTEHbZuMYz5JJz1bTXESVnq8HcfZjbDzNabWU0br5uZ/cjMVpnZUjM7Pv8xs3fYtDRC1xIWERHpFHdY+HO4YzxseY2akdfAv92qIljKVi6LgO4BJrbz+iRgePA1Fbi9+7Ha5jQ1y6kQFhERydmOOvj1x+Gx/4Yh74PP/pWN/U8OO5VIqDoshN19LlDXzibnAPd6xnygj5kdmq+A+9EaYRGRbjOzoWY2x8xWmNlyM/ti2JmkgF6dC7ePg5d+n2mIu+gRXRVChPysER4MvJ71eE3w3LqWG5rZVDKzxgwYMIBUKtXpnTUEM8HPzv8bvar6dPr9+VRfX9+lP0Op5oDoZIlKDlCWKOeAaGUpsgbgKndfZGYHAwvN7El3fyHsYJJH+9wh7ii44H41xIlkyUch3Nrt3Vq9B7K73wncCTB27Fivrq7u9M4er3kcgDFjjqf/oYd1+v35lEql6MqfoVRzQHSyRCUHKEuUc0C0shSTu68jmLBw97fNbAWZSQwVwqUi+w5xaogTaVU+CuE1wNCsx0OAtXn43DY0LY3QGmERkXwws2HAaODZFs93+7d4UZpxj0qWgudw59B1T3LUqp+RTvTgpZHXsPHgk+GvzxU/SydEJUtUckB0skQlB+Q/Sz4K4UeBz5vZA8CJwNZgpqEwgjXCjVojLCLSbWZWBcwELnf3bdmv5eO3eFGacY9KloLm2FGXaYZb+RgccRrJc+9gVO+223aickwgOlmikgOikyUqOSD/WToshM3sfqAa6Gdma4BvAJUA7j4dmA1MBlYBO4BP5S1dq4Ey/X3eqEJYRKQ7zKySTBF8n7vPCjuPdNPqP8PDnwnuEDcNTv687hAn0oEOC2F3v6CD1x34XN4SdaDpFstaGiEi0nVmZsBdwAp3vyXsPNINDXtgzg3wlx+pIU6kk2J3ZzlrLoQ1Iywi0g3jgIuAZWa2OHjuK+4+O7xI0mkbX4aZlwQNcZ+Cs76lhjiRTohdIdzULOcqhEVEuszdn6H1q/5IHLjDonvhD9dARU/4j/tgxNlhpxKJnfgVwpoRFhGRctbUELci0xDHuXdAOw1xItK2+BXCwcL/RjXLiYhIuVFDnEhexa8QDi6f5o1qlhMRkTKhhjiRgohdIewkAUinG0JOIiIiUgT7NcR9G3ocGHYqkZIQu0K4+TrCKoRFRKSUucOin8MfroWKXmqIEymA+BXCicyMMI0qhEVEpETtqINHvwAvPg7vroaPTFdDnEgBxK8QtmBphAphEREpRdkNcWfeACd9Tg1xIgUSu0LYtTRCRERKUcMeeHoa/PX/oN9w+M8H4ND3hp1KpKTFrhBumhF2XT5NRERKxcaXYebFsG6JGuJEiih+hXBCM8IiIlIiWjbETfkVvOdDYacSKRvxK4S1NEJEREqBGuJEQhe/Qji4aoTrFssiIhJXq1NBQ9xGNcSJhCh+hTBNa4Q1IywiIjHTsId3v3IPpB4JGuIeVEOcSIjiVwg3/Ys5vTfcHCIiIp2xYSXMuoTD1i2BsZ+GM7+lhjiRkMXv9zBNV41Ip0MOIiIikgN3WHA33HEqbHmdZaO+Amf/QEWwSATEbkbYmu4sp2Y5ERGJun0a4k6Hj9zOpkUvhZ1KRAKxK4SbrhqhQlhERCJtn4a4b8FJ/y9Y3qdCWCQqYlgI66oRIiISYfvdIe7XcOixYacSkVbErxBO6KoRIiISURtWZu4Q9+ZSNcSJxEDsCmELrhphWhohIiJR4Q4L78ncIa7yAN0hTiQm4lcIN99ZTksjREQkAlo2xJ07HQ4eGHYqEclB7Aphb7pqhKsQFhGRkL0yBx75LOzY1KIhTkTiIHaFsAXNcloaISIioWnYA09/M2iI+xc1xInEVPwK4YSuGiEiIiHapyHuYjjzBjXEicRUTr+/MbOJZvaSma0ys2taeb3azLaa2eLg67r8Rw3ohhoiIhIGd1gwI3OHuK1rYMr9cPYtKoJFYqzDGWHLrEX4MfBBYA3wnJk96u4vtNh0nrufXYCM+0jqhhoiIt1mZjOAs4H17j4q7DyRt30TPPbfaogTKTG5zAifAKxy99Xuvgd4ADinsLHalkwajW5aGiEi0j33ABPDDhELr8yB298PLz+RaYi7cJaKYJESkcsa4cHA61mP1wAntrLdyWa2BFgLXO3uy1tuYGZTgakAAwYMIJVKdTrwzu3baSDJlrpNXXp/PtXX14eeIUo5IDpZopIDlCXKOSBaWYrJ3eea2bCwc0Raw+6sO8QdDf/1kBriREqMuXv7G5idD5zl7pcEjy8CTnD3L2Rt0xtIu3u9mU0GbnX34e197tixY33BggWdDpxKpThxzvksG3Q+J1z2k06/P59SqRTV1dWhZohSDohOlqjkAGWJcg7oehYzW+juY/OfqHiCQvjxtpZGtJi8GPPAAw90eh/19fVUVVV1J2bedCbLgdvXMGLF9zm4fjVvDJrEK0d+inSyZ9FzFJqyRDcHRCdLVHJA17OcfvrprY7ZucwIrwGGZj0eQmbWt5m7b8v6fraZ/cTM+rn7xk4nzUEjSS2NEBEpMHe/E7gTMpMXXfkHQ+z+0eMOC++GZ76SaYKbcj+D3zOZwcXOUSTKEt0cEJ0sUckB+c+SSyH8HDDczI4A3gCmAP+ZvYGZDQTecnc3sxPIrD3elLeULTRaQs1yIiKSX9s3Ze4Q99Lv1BAnUiY6LITdvcHMPg/8EUgCM9x9uZl9Jnh9OvBR4LNm1gDsBKZ4R2suuiFNUoWwiIjkzytz4OHPwM46OOvbcOJndYc4kTKQ0w013H02MLvFc9Ozvr8NuC2/0drWSEJLI0REusHM7geqgX5mtgb4hrvfFW6qEDTshqe+CX+7TQ1xImUodneWA2i0JKZCWESky9z9grAzhE53iBMpe7EshJ0kuJZGiIhIFzQ1xP0haIi74AE4elLYqUQkBLEshNOWhHQ67BgiIhI32Q1xR54BH7ldDXEiZSyehTAJTM1yIiLSCe+oWwy3XxY0xH0HTvyMGuJEylw8C2HT0ggREclR0BD33qW3Qf/3wIW/gYH/GnYqEYmA2BbC5loaISIiHdjwUtAQt4w3Bk1i8CdnqCFORJrFshB2S2pphIiItM0dFsyAP361uSHu5XUHMFhFsIhkiW8h7Lp8moiItKKthrh1qbCTiUjExLIQTlsFCc0Ii4hIS688HdwhbrMa4kSkQ7EshBsTlVQ07Ak7hoiIREX2HeL6vwcunKmGOBHpUEwL4R5U+PawY4iISBRkNcTxvksyd4irPCDsVCISA7EthHv53rBjiIhImPZriHsQjp4YdioRiZFYFsLpRA8qXUsjRETK1vZN8Ojn4aXZcOQHgoa4AWGnEpGYiWUh7IkeVKJmORGRsqSGOBHJk1gWwulkJZVaGiEiUl7UECcieRbLQtgTPeiBCmERkbKxT0PcpXDmNDXEiUi3xbMQTvagUoWwiEjpc4cFdwUNcQepIU5E8iqWhTDJnvSggb2NaSqTWhcmIlKStm8M7hCnhjgRKYxYFsJW2YtKa2Trrt0ccpB+NSYiUnLUECciRRDPQrhnFQA7t29TISwiUkr2a4ibBQNHhZ1KREpULAvhRK+gEK7fBu/Sr8lERErC+hdh5iXwlhriRKQ4YlkIJ3seDMCeHdtCTiIiIt22T0NclRriRKRoYlkIVxyQmRHevePtkJOIiEi3ZDfEHTUBzvmJGuJEpGhiWQj3OrA3kFkjLCIiMbXqKXjks5mGuIk3wgmXqSFORIoqloVwVZ++AOx+e1PISUREpNMadsOfrof5P1ZDnIiEKpaFcO93HgpAw7YNIScREZHOeHPV8+x58NMctnc1a4ZfxJCPfU8NcSISmlgWwj0PGQhAuv6tkJOIiEguPJ1m8cO3MGLpjWznAKamv8S27WfwgIpgEQlRTouxzGyimb1kZqvM7JpWXjcz+1Hw+lIzOz7/UbNU9GSr9SaxbU1BdyMiUso6GtvzZevGtSz7/ocYvWwaK3q+l12XzOPwk85l0T+2sGtvY6F2KyLSoQ4LYTNLAj8GJgHHABeY2TEtNpsEDA++pgK35znnfjYedBQDt79I/e6GQu9KRKTk5Di2d9u2Vxey97aTeU/935l35FUc++UnGDx0GOOH92dPQ5rbU69Qu3E72zWWi0gIclkacQKwyt1XA5jZA8A5wAtZ25wD3OvuDsw3sz5mdqi7r8t74kCPEZM58rkbeP3bI6hLVOJmOABWqF3u57B0mtq54Xc4RyUHRCdLVHKAskQ5B8DbfSdBdXXYMcKQy9jeLfN/eT0ffu0WXk0cxuaPPMj4Y09qfm388H6cfnR/bn3qZW596mUA+lX14KCeFSQTRsIs76P59h07OGjRn/P8qfHNAcoS5RwQnSxRyQEwuMfuvA7ZuRTCg4HXsx6vAU7MYZvBwD6FsJlNJTNjzIABA0ilUp2MC/X19aRSKeyA0awZcCGJLatpTDuZMthbf1MbT3eXm2NWvMI76jkgOlmikgOUJco5AHbZAV0ai0pAh2N7d8fsbclBPNHzLBJjPk1F3S7eaPH+Cw933te7F3W70mzZ5azfmWZP427SDukCjNtVPdMkbWf+PzimOUBZopwDopMlKjkAqpIN+R2z3b3dL+B84GdZjy8C/q/FNr8DTsl6/BQwpr3PHTNmjHfFnDlzuvS+QohKlqjkcI9OlqjkcFeW1kQlh3vXswALvIPxM8pfuYzt2V8as/MnKjnclaU1UcnhHp0sUcnhnv8xO5ffTa4BhmY9HgKs7cI2IiISHRq3RaTs5VIIPwcMN7MjzKwHMAV4tMU2jwIfD64ecRKw1Qu4PlhERLotl7FdRKSkdbhG2N0bzOzzwB+BJDDD3Zeb2WeC16cDs4HJwCpgB/CpwkUWEZHuamtsDzmWiEhR5XRDDXefTabYzX5uetb3Dnwuv9FERKSQWhvbRUTKSTSuXyQiIiIiUmQqhEVERESkLKkQFhEREZGypEJYRERERMqSZfrcQtix2QbgtS68tR+wMc9xuioqWaKSA6KTJSo5QFlaE5Uc0PUsh7t7/3yHiSqN2XkVlRygLK2JSg6ITpao5IA8j9mhFcJdZWYL3H1s2DkgOlmikgOikyUqOUBZopwDopWlFEXp+EYlS1RygLJEOQdEJ0tUckD+s2hphIiIiIiUJRXCIiIiIlKW4lgI3xl2gCxRyRKVHBCdLFHJAcrSmqjkgGhlKUVROr5RyRKVHKAsrYlKDohOlqjkgDxnid0aYRERERGRfIjjjLCIiIiISLepEBYRERGRshSbQtjMJprZS2a2ysyuCWH/tWa2zMwWm9mC4Lm+Zvakmb0c/PcdBdr3DDNbb2Y1Wc+1uW8zuzY4Ti+Z2VkFzvG/ZvZGcFwWm9nkQucIPnuomc0xsxVmttzMvhg8X9Tj0k6Ooh8XM+tlZn83syVBluuD54t9TNrKEcrPSvD5STN73sweDx4X/e9POQpz3NaY3W6WMMYnjdn7Z4nEmN1BlrD+H1+8MdvdI/8FJIFXgHcDPYAlwDFFzlAL9Gvx3E3ANcH31wDfLdC+TwWOB2o62jdwTHB8egJHBMctWcAc/wtc3cq2BcsRfP6hwPHB9wcDK4N9FvW4tJOj6McFMKAq+L4SeBY4KYRj0laOUH5Wgn1cCfwKeDysvz/l9kXI4zYas9vLEsb4pDF7/8+OxJjdQZaiH5fg84s2ZsdlRvgEYJW7r3b3PcADwDkhZ4JMhp8H3/8c+EghduLuc4G6HPd9DvCAu+9291eBVWSOX6FytKVgOYIs69x9UfD928AKYDBFPi7t5GhLIc+Pu3t98LAy+HKKf0zaytGWgv6smNkQ4EPAz1rss6h/f8pQFMftshqz28nSlkKOTxqz988SiTG7gyxtKViWYo/ZcSmEBwOvZz1eQ/s/uIXgwBNmttDMpgbPDXD3dZD5ywW8q4h52tp3GMfq82a2NPgVXNOvK4qWw8yGAaPJ/As2tOPSIgeEcFyCXyctBtYDT7p7KMekjRwQzs/KD4EvAems56L096dUhX0sNWa3L7RxW2P2PhkiMWa3kwWKf1x+SBHH7LgUwtbKc8W+7ts4dz8emAR8zsxOLfL+c1XsY3U7cCRwHLAO+H4xc5hZFTATuNzdt7W3aSHztJIjlOPi7o3ufhwwBDjBzEa1F7tQWdrIUfRjYmZnA+vdfWGubylUljIU9rHUmN220MZtjdktPigiY3Y7WYp6XMIYs+NSCK8BhmY9HgKsLWYAd18b/Hc98DCZqfe3zOxQgOC/64sYqa19F/VYuftbwV+eNPBT/vkriYLnMLNKMgPZfe4+K3i66MeltRxhHpdg/1uAFDCREH9WsnOEdEzGAR82s1oyv5o/w8x+SUT+/pS4UI+lxuy2hTU+acxuW1TG7JZZQjguRR+z41IIPwcMN7MjzKwHMAV4tFg7N7ODzOzgpu+BM4GaIMMngs0+Afy2WJna2fejwBQz62lmRwDDgb8XKkTTD2bgXDLHpeA5zMyAu4AV7n5L1ktFPS5t5QjjuJhZfzPrE3x/ADABeJHiH5NWc4RxTNz9Wncf4u7DyIwbT7v7hUTk70+JC23c1pjdvpDGJ43Z+2eJxJjdXpZiH5dQxmzPU4dfob+AyWS6O18Bvlrkfb+bTFfiEmB50/6BdwJPAS8H/+1boP3fT+ZXEnvJ/Ovn4vb2DXw1OE4vAZMKnOMXwDJgafADeWihcwSffQqZX38sBRYHX5OLfVzayVH04wIcCzwf7LMGuK6jn9MCHZO2coTys5K1j2r+2YFc9L8/5fhFSOM2GrM7yhLG+KQxe/8skRizO8gS2rhNkcZs3WJZRERERMpSXJZGiIiIiIjklQphERERESlLKoRFREREpCypEBYRERGRsqRCWERERETKkgphERERESlLKoRFREREpCz9fyclYPuM+L0HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,4))\n",
    "ax1 =  fig.add_subplot(121)\n",
    "ax2 =  fig.add_subplot(122)\n",
    "\n",
    "x_plot = x.value\n",
    "p_n,      = ax1.plot(x_plot, n.value, '-', label = '$n_i$')\n",
    "p_v,      = ax2.plot(x_plot, v.value, '-',label = '$v_i$')\n",
    "a_n,      = ax1.plot(x_plot, n0*n_ana(x.value,t), '-', label = '$n_i^{an}$')\n",
    "a_v,      = ax2.plot(x_plot, 1.5*x.value/t, '-',label = '$v_i^{an}$')\n",
    "\n",
    "ax1.grid()\n",
    "ax2.grid()\n",
    "ax1.legend(framealpha=0.5)\n",
    "ax2.legend(framealpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_list = []; n_list = []; u_list = []; n_sweeps = []\n",
    "\n",
    "def update(i):\n",
    "    update_progress('Progress',1.*i/steps)\n",
    "    global t,tau,T\n",
    "    n.updateOld()\n",
    "    v.updateOld()\n",
    "    T.setValue(t*tau)\n",
    "    res = 1e100\n",
    "    isweep = 0\n",
    "    while res > 1.:\n",
    "        isweep+=1\n",
    "        v_face[:] = v.arithmeticFaceValue\n",
    "        v_face[..., mesh.facesLeft.value] = 0.\n",
    "        res = eqn.sweep(dt=dt)    \n",
    "        if isweep > 100:\n",
    "            print(i,'No converge',res)\n",
    "            break\n",
    "    n_sweeps.append(isweep)\n",
    "    #print(isweep, res)\n",
    "    t_list.append(t)\n",
    "    n_list.append(n.value.copy())\n",
    "    u_list.append(v.value.copy())\n",
    "    t = t + dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Progress: [####################################################################################################] 100% Computation walltime  320.7130470275879\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "for i in range(steps):\n",
    "    update(i)\n",
    "print('Computation walltime ', time.time()-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fb2540241c0>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAapUlEQVR4nO3df2zc9Z3n8ed7ftjjxI6dECfYCWmApF04KKHrFtrssiw/bmngStGpd+1u2dzCXbjTItFupRa60omquhPXFuju3h5SKIhs22uPPUAgutU1CmS5arvQBEJILkCgpDTYtQ3BjhNsxz/e98d8J54kdvwdZ8bf+cy8HpLlme/MxO+PBa988v5+vp+vuTsiIhKeVNIFiIjI3CjARUQCpQAXEQmUAlxEJFAKcBGRQGXm84ctXbrUV69ePZ8/UkQkeDt37nzX3dtPPj6vAb569Wp27Ngxnz9SRCR4Zvbr6Y6rhSIiEigFuIhIoBTgIiKBUoCLiARKAS4iEigFuIhIoBTgIiKBCiLAt+3r5YHtbyZdhohIVQkiwLe/1s+D//dXSZchIlJVggjwdMoYn5hMugwRkaoSO8DNLG1mL5nZ09Hzu83sHTPbFX1tqFSR6ZQxMak7B4mIFCtlL5Q7gH3AoqJj97v7d8pb0qkyKWNCt34TETlBrBm4ma0Erge+V9lypqcZuIjIqeK2UL4LfBU4uRF9u5ntNrOHzWxxWSsrkk4Z4wpwEZETzBrgZnYD0OfuO0966QHgfGAd0APcO8PnN5nZDjPb0d/fP6ci0ynDHSYV4iIix8WZga8HPmNmB4AfA1eZ2Q/cvdfdJ9x9EngQ+MR0H3b3ze7e5e5d7e2n7EceSyZlAOqDi4gUmTXA3f0ud1/p7quBzwPPuPsXzayj6G03AXsqVCOpQoBrBi4ictyZ3JHnW2a2DnDgAHBbOQqaTkYBLiJyipIC3N23A9ujxzdXoJ5ppVP5fyjoRKaIyJQwrsTMT8A1AxcRKRJGgKfzZSrARUSmhBHgph64iMjJggjwwknM8UltaCUiUhBEgKejAFd+i4hMCSrANQMXEZkSVICrBy4iMiWIANel9CIipwoiwAuX0o9PKMBFRAqCCHBdSi8icqogAjytFoqIyCnCCnDNwEVEjlOAi4gEKogAz6S0F4qIyMmCCPBoLyttJysiUiSQAM+XqXtiiohMCSLApzazUoCLiBQEEeCp49vJai8UEZGC2AFuZmkze8nMno6eLzGzrWa2P/q+uFJFZtKFAK/UTxARCU8pM/A7gH1Fz+8Etrn7WmBb9LwitBuhiMipYgW4ma0Erge+V3T4RmBL9HgL8NmyVlZEd+QRETlV3Bn4d4GvAsVT4OXu3gMQfV823QfNbJOZ7TCzHf39/XMqUhfyiIicatYAN7MbgD533zmXH+Dum929y9272tvb5/JHFPXAFeAiIgWZGO9ZD3zGzDYAOWCRmf0A6DWzDnfvMbMOoK9SRRZaKFpGKCIyZdYZuLvf5e4r3X018HngGXf/IvAUsDF620bgyUoVefyemNqNUETkuDNZB34PcK2Z7QeujZ5XRGEvFN3QQURkSpwWynHuvh3YHj1+D7i6/CWdKspv9cBFRIoEcSXm8d0I1UIRETkuiADXMkIRkVMFFeDqgYuITAkiwKP8VgtFRKRIEAFuZmRSxrh2sxIROS6IAIf81Zi6kEdEZEowAd6QTnFsXDNwEZGCYAK8MZtmVAEuInJcMAHekE4xOj6RdBkiIlUjmABvzKiFIiJSLJgAb1CAi4icIJgAb8yk1AMXESkSTIBrBi4icqJgArwxk+aYLuQRETmupO1kk9SQSTE4PJZ0GSIVMTg8xvjEJGbG4gVZLLoLlcjphBPgWkYoNerZV/v4s0d+efz51zf8DpuuOD/BiiQU4bRQsuqBS23a2z0IwN3/6kJy2RTdAyMJVyShCCbAdSm91KruwRGWLGzg360/l6ZsWvveS2yzBriZ5czsBTN72cz2mtk3ouN3m9k7ZrYr+tpQyUIbtIxQalTPwDAdrTkgv/e9tk2WuOL0wEeBq9z9iJllgZ+b2U+j1+539+9UrrwpjZm0ZuBSk7oHRjhnyQIAUmZMagYuMc0a4O7uwJHoaTb6mvf/whoyKUa1jFACMjI2wVf+/mXeP3rstO97692jXH7eEiCagSvAJaZYPXAzS5vZLqAP2Oruz0cv3W5mu83sYTNbPMNnN5nZDjPb0d/fP+dCCxfyuP55KYF4vXeIn+zu4d0jo4xNTM74tW5VG5++uAPIz8DVQpG4Yi0jdPcJYJ2ZtQFPmNlFwAPAN8nPxr8J3AvcMs1nNwObAbq6uub8X2ZjJv93zbGJSRoz6bn+MSLzprCa5L5/s46LVrTG+kw6pRaKxFfSKhR3HwC2A9e5e6+7T7j7JPAg8InylzfleICrDy6B6BkcBjh+gjKO/EnMSlUktSbOKpT2aOaNmTUB1wCvmllH0dtuAvZUpMJIQxTgWokioegZHKExk2LJwobYn0kZmoFLbHFaKB3AFjNLkw/8R939aTP7vpmtI99COQDcVrEq0QxcSvfMq730Hh5N7Oe/8NYhOlpzJV0Wr5OYUoo4q1B2A5dOc/zmilQ0g0Lfe2RMl9PL7PqHRrnlkR1Jl8GnLzq7pPfrJKaUIpi9UJob86UeGR1PuBIJwcH3PwDg3s9dwvo1SxOrY2lz/PYJaB24lCaYAG/J5UsdGlGAy+x6BvMrQC7oWMTZJZxETFo6ZUxqBi4xBbMXSksuCyjAJZ7ugfwKkM62cMIbIKVVKFKCAGfg2hO8Hh0dHS9pZvqbQx/QlE3T2pStYFXll9YqFClBgAGuGXi9+V+/fJuvPfZKyZ9bs6w5uBsjaBWKlCKYAF+ok5h16+WDg7Q0ZrjjmrUlfe7SVdPu7lDVtApFShFMgGfTKZqyabVQ6lDPwDCrzlrAv//985IupeLSKdO1DhJbMCcxAZpzGbVQ6lDP4AgdrU1JlzEvtB+4lCKoAG/JZRhSC6XuvDMwzIrAVpPMldaBSymCaaEAtDRqBl7tBofH+PoTr3C0TH/RuudPXHe0aQYucrKgAnzJwgb6hpLb20Jmt/PXh/jJ7h4+vLyZpmx5tv3t+tBirljbXpY/q9qlzNB9SySuoAK8o62Jlw8OJl2GnEZhD+y/u+WyoK6ArBbplNaBS3xB9cA7W3McOnqM4WPa0Kpa9QwOk0kZ7S2NSZcSJLVQpBRBBXhhJUJho3ypPt0DIyxflCOdCusCmmphpr1QJL6gWiidbYUAH+G89uaEq6ldg8NjPPtq35yuCNzbPRjc/iPVJK1VKFKCoAK80FPtPTyScCW17e/+6QD3bn19zp//k8tWlbGa+qIWipQiqABfpP1Q5sXbhz6gvaWRx/7jp+b0ec3A5y6/DjzpKiQUQQV4s3YknBc9gyOsaGti1VkLki6l7qRTaDMriS3OTY1zZvaCmb1sZnvN7BvR8SVmttXM9kffK75zUGMmTUMmpasxK6x7cFiz6ISohSKliLMKZRS4yt0vAdYB15nZ5cCdwDZ3Xwtsi55XnK7GLJ27Mzo+EetrZGyCnoH62Xuk2uhSeilFnJsaO3AkepqNvhy4EbgyOr4F2A58rewVnqQll+GIArwk/+Un+/jez98q6TMr6uTS9WqjGbiUIlYP3MzSwE5gDfC37v68mS139x4Ad+8xs2UzfHYTsAlg1aozX53QksuqB16inW+/z3lLF/Kvf3dlrPdn08ZNl66ocFUynfyl9ApwiSdWgLv7BLDOzNqAJ8zsorg/wN03A5sBurq6zvi/zGa1UErWMzDC+jVL+fM/XJN0KTKLdEotFImvpCsx3X2AfKvkOqDXzDoAou995S5uOi25jO7KU4LxiUn6hkZ0UjIQaqFIKeKsQmmPZt6YWRNwDfAq8BSwMXrbRuDJCtV4gnwLRQEeV+/QKJM+dRWrVDetA5dSxGmhdABboj54CnjU3Z82s18Aj5rZrcDbwOcqWOdxLblMXfbA//sz+9k9h50YD0e/qw7tDBiEdArNwCW2OKtQdgOXTnP8PeDqShR1OoUWirsHd8fxuZqcdP5q235amxpY2txQ8uc/vnoxl6xsK39hUnYpbWYlJQjqSkzIn8ScdDh6bILmxuDKn5N3j4wyNuHccfUabv7k6qTLkQpKmeFOXU1QZO6C2k4W8j1woK7Wgr8zkN8+VxfX1L7CNrxaSihxBBjg9bcfSs9gfvdFnYisfccDXG0UiSG4HsTxDa1qbCnh671DvPbboWlfe+71fkC7/NWDVNQ20UoUiSO4AK/VLWVv+/5O3nr36Iyvt7c00tqUnceKJAnp6N/EmoFLHMEFeHNjPsRqqYUyMen85tAH/PFlq7hl/epp37O0uVEntepAYQauHrjEEVyAF3rgtXQSs39olPFJ54KORaxZ1pJ0OZKgQg9cl9NLHAGfxKydAO+ObtK8Qj3uuqeTmFKK4AJ8YUMGs9pqofQM5FeZaJmgTJ3EVIDL7IJroaRSRktjhr9+5g0e+Mc3eWjjx7niw+1Jl8XdT+3lkX86cEZ/RqcCvO5pBi6lCC7AAb79uUvYfXCAv332TV777VBVBPgLbx3i/PaF3PDRzjl9/pwlC2hdoFUm9S6tk5hSgiAD/I/+xdlce8Fy/sf2N6umldI9OMz1F3fw5Ws/nHQpErBUSuvAJb7geuAFqZTR3JDhcBWczBw+NsHAB2O6UlLOWJTf2tBKYgk2wKF6bu5QWEWiKyXlTKkHLqUIsoVSkOT9MQ++/wH3/ex1jk1M8v4HxwCtIpEzp1UoUoqgA7w5l9z9Mbf+v14ef+kdzl26EDO4dFUbF3QsSqQWqR3ZdD7AxyYU4DK7oAO8JZfh0NFjifzsnsERGjMpnvnKH+gSdymbbLQZyrjOYkoMce6JeY6ZPWtm+8xsr5ndER2/28zeMbNd0deGypd7oiTvj9k9MExHa07hLWWViQJ8bEIBLrOLMwMfB77i7i+aWQuw08y2Rq/d7+7fqVx5p9fcmFwLpWdwRD1vKbtCC+XYuFooMrs498TsAXqix0Nmtg9YUenC4liUy3B4eIz/s/e38/6zf/3e0aq4gEhqS4NaKFKCknrgZraa/A2OnwfWA7eb2Z8CO8jP0t+f5jObgE0Aq1atOtN6T9DRmuPYxCS3fX9nWf/cuNYsa07k50rtyqqFIiWIHeBm1gw8BnzJ3Q+b2QPANwGPvt8L3HLy59x9M7AZoKurq6z/LvzTT67m8vPPSuSy45QZaxXgUmaFAFcLReKIFeBmliUf3j9098cB3L236PUHgacrUuFppFLG75ytpXtSOxoyhWWEmoHL7OKsQjHgIWCfu99XdLyj6G03AXvKX55Ifcmk1EKR+OLMwNcDNwOvmNmu6NjXgS+Y2TryLZQDwG0VqE+krmQzCnCJL84qlJ8D0y12/ofylyNS33QlppQi6M2sRGpNg1ahSAkU4CJVRMsIpRQKcJEqMhXgaqHI7BTgIlVk6lJ6zcBldgpwkSpiZmRSphaKxKIAF6ky2XRKAS6xKMBFqkw2beqBSywKcJEq05DRDFziUYCLVBm1UCQuBbhIlcmohSIxKcBFqkw2neKYZuASgwJcpMo0pFOMaR24xKAAF6ky2XSK8QRuUiLhUYCLVJn8MkLNwGV2CnCRKpNNp3QpvcSiABepMg0ZncSUeBTgIlWmKZtm+NhE0mVIABTgIlWmJZdlaGQ86TIkAHFuanyOmT1rZvvMbK+Z3REdX2JmW81sf/R9ceXLFal9LbkMh0fGki5DAhBnBj4OfMXdLwAuB/7czC4E7gS2uftaYFv0XETOUEsuw5HRcdy1lFBOb9YAd/ced38xejwE7ANWADcCW6K3bQE+W6EaRepKc2MGdziqPrjMoqQeuJmtBi4FngeWu3sP5EMeWDbDZzaZ2Q4z29Hf33+G5YrUvpZcFoAj6oPLLGIHuJk1A48BX3L3w3E/5+6b3b3L3bva29vnUqNIXWnJZQAYUh9cZhErwM0sSz68f+juj0eHe82sI3q9A+irTIki9aU5CvDDmoHLLOKsQjHgIWCfu99X9NJTwMbo8UbgyfKXJ1J/FkUBfmRUAS6nl4nxnvXAzcArZrYrOvZ14B7gUTO7FXgb+FxFKhSpM82N+R64Wigym1kD3N1/DtgML19d3nJEZKoHrhm4nJ6uxBSpMoUA1yoUmY0CXKTKLGzIYKYWisxOAS5SZVIpo7khw5BOYsosFOAiVagll1EPXGalABepQs25jFooMisFuEgVaslltQ5cZqUAF6lCzY1qocjsFOAiVUg9cIlDAS5ShXRXHolDAS5ShVp0ElNiUICLVKG2BVlGxyc5qhOZchoKcJEq1NnaBEDP4EjClUg1U4CLVKGO1hwA3QPDCVci1UwBLlKFOtsKM3AFuMxMAS5ShZYvymEG3QNqocjMFOAiVaghk6K9uVEzcDktBbhIlepoa9IMXE5LAS5SpTpbc3RrBi6nEeemxg+bWZ+Z7Sk6dreZvWNmu6KvDZUtU6T+dLQ20TMwgrsnXYpUqTgz8EeA66Y5fr+7r4u+/qG8ZYlIZ1uO4bEJBod1RaZMb9YAd/fngEPzUIuIFCksJXxHa8FlBmfSA7/dzHZHLZbFM73JzDaZ2Q4z29Hf338GP06kvhQu5unRiUyZwVwD/AHgfGAd0APcO9Mb3X2zu3e5e1d7e/scf5xI/dHFPDKbOQW4u/e6+4S7TwIPAp8ob1ki0t7cSDZtdGs/FJnBnALczDqKnt4E7JnpvSIyN6mUsXxRTvuhyIziLCP8EfAL4CNmdtDMbgW+ZWavmNlu4A+BL1e4TpG61NnaxJO7utVGkWllZnuDu39hmsMPVaAWETnJJ88/ixcOHOKxnQe5/aq1SZcjVUZXYopUsS9f+2EWL8iqDy7TUoCLVLnOtib1wWVaCnCRKle4pF7kZApwkSrX2ZZfifJG35D2RZETKMBFqtyHzlrI0Og419z3HC++/X7S5UgVUYCLVLk/uWwV9//bSwB4s+9owtVINVGAi1S5XDbN9Rd3YqaNreRECnCRAOgWazIdBbhIIDramujRenApMuuVmCJSHTpbc2zb18cffPtZAC5Z2cZff+HShKuSJGkGLhKIjZ9azYaLz+bSc9pY0JDh6d3djE1MJl2WJEgzcJFAXH7eWVx+3lkA/OiFt7nr8VfoPTzCysULEq5MkqIZuEiApm72oJ54PVOAiwSoM7rdmvZIqW9qoYgEqCOagf/9joO83js07Xs+tmoxV1+wfD7LknmmABcJUHNjhktWtvLPv3qPf/7Ve6e8PuHOspZGnleA1zQFuEignrz992Z87b6tr/M3z+zn2PgkDRl1SmtVnFuqPWxmfWa2p+jYEjPbamb7o++LK1umiJRiRVsOd+g9rJOctSzOX82PANeddOxOYJu7rwW2Rc9FpEp0tGqVSj2Ic0/M58xs9UmHbwSujB5vAbYDXytnYSIyd51t+VUqL779PouaZu+UZlLGeUubSaWs0qVVjb7DIxz64Ni8/byVixfQ3FjervVc/7Tl7t4D4O49ZrasjDWJyBnqbGsikzLu+emr3PPTV2N95r/edDF/fNmqCldWHY6OjnPFt59lZGz+rmR95M8+zpUfKW9UVvwkppltAjYBrFpVH/9xiCRtQUOG//2fPkVPzHXif/Hoy7zRd6TCVVWPdwaGGRmb5D/8/rl8bNX8nMK7sHNR2f/MuQZ4r5l1RLPvDqBvpje6+2ZgM0BXV5fuByUyT9ad08a6c9pivfc7P3utri4KKuyrft1FZ/O7H1qScDVzN9f1RU8BG6PHG4Eny1OOiCShs62prvYaL9wkunCyN1RxlhH+CPgF8BEzO2hmtwL3ANea2X7g2ui5iASqozVHdx2tWOkZHCZlsKylMelSzkicVShfmOGlq8tci4gkpLOtif6hUa697x+TLmVe9A2NsnxRjkw67IucdCWmiHDDRzt4s/8oE5P1sb/42uXNrF+zNOkyzpgCXERYs6yFv9HdfYIT9r8fRETqmAJcRCRQCnARkUApwEVEAqUAFxEJlAJcRCRQCnARkUApwEVEAmXu87dBoJn1A7+e48eXAu+WsZwQaMz1QWOuD2cy5g+5e/vJB+c1wM+Eme1w966k65hPGnN90JjrQyXGrBaKiEigFOAiIoEKKcA3J11AAjTm+qAx14eyjzmYHriIiJwopBm4iIgUUYCLiAQqiAA3s+vM7DUze8PM7ky6nnIxs4fNrM/M9hQdW2JmW81sf/R9cdFrd0W/g9fM7I+SqXruzOwcM3vWzPaZ2V4zuyM6XstjzpnZC2b2cjTmb0THa3bMBWaWNrOXzOzp6HlNj9nMDpjZK2a2y8x2RMcqO2Z3r+ovIA28CZwHNAAvAxcmXVeZxnYF8DFgT9GxbwF3Ro/vBP5b9PjCaOyNwLnR7ySd9BhKHG8H8LHocQvwejSuWh6zAc3R4yzwPHB5LY+5aOx/AfxP4OnoeU2PGTgALD3pWEXHHMIM/BPAG+7+K3c/BvwYuDHhmsrC3Z8DDp10+EZgS/R4C/DZouM/dvdRd38LeIP87yYY7t7j7i9Gj4eAfcAKanvM7u5HoqfZ6Mup4TEDmNlK4Hrge0WHa3rMM6jomEMI8BXAb4qeH4yO1arl7t4D+cADlkXHa+r3YGargUvJz0hresxRK2EX0AdsdfeaHzPwXeCrQPFdkmt9zA78zMx2mtmm6FhFxxzCTY1tmmP1uPaxZn4PZtYMPAZ8yd0Pm003tPxbpzkW3JjdfQJYZ2ZtwBNmdtFp3h78mM3sBqDP3Xea2ZVxPjLNsaDGHFnv7t1mtgzYamavnua9ZRlzCDPwg8A5Rc9XAt0J1TIfes2sAyD63hcdr4nfg5llyYf3D9398ehwTY+5wN0HgO3AddT2mNcDnzGzA+RbnleZ2Q+o7THj7t3R9z7gCfItkYqOOYQA/yWw1szONbMG4PPAUwnXVElPARujxxuBJ4uOf97MGs3sXGAt8EIC9c2Z5afaDwH73P2+opdqeczt0cwbM2sCrgFepYbH7O53uftKd19N/v/XZ9z9i9TwmM1soZm1FB4D/xLYQ6XHnPSZ25hndzeQX7HwJvCXSddTxnH9COgBxsj/jXwrcBawDdgffV9S9P6/jH4HrwGfTrr+OYz398j/M3E3sCv62lDjY/4o8FI05j3Af46O1+yYTxr/lUytQqnZMZNfJfdy9LW3kFOVHrMupRcRCVQILRQREZmGAlxEJFAKcBGRQCnARUQCpQAXEQmUAlxEJFAKcBGRQP1/qdwB6uNMIloAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(n_sweeps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pickle\n",
    "#pickle.dump((t_list,n_list,u_list),open( \"save.p\", \"wb\" ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import pickle\n",
    "#t_list,n_list,u_list = pickle.load(open( \"save.p\", \"rb\" ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_func(i):   \n",
    "    p_n.set_ydata(n_list[i])\n",
    "    p_v.set_ydata(u_list[i])\n",
    "    a_n.set_ydata(n0*n_ana(x.value,t_list[i]))\n",
    "    a_v.set_ydata(1.5*x.value/t_list[i])\n",
    "    ax1.set_ylim((0,max(n_list[i])*1.2))\n",
    "    ax2.set_ylim((0,max(1.5*x.value/t_list[i])*1.2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#animation2 = FuncAnimation(fig, plot_func, frames=steps, repeat=False)\n",
    "#animation2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAEDCAYAAADgABhKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABVf0lEQVR4nO3dd3gV1drG4d+bHnoPvQkoEJUmRVSCKCIWBEVQD6JyxAIeUFRs2FCPWLChKApSRAErqKAiGjuCIFWko/QWWoD09f2xh/MFDBBCktlJnvu69pW9156Z/WQmDG8ma9Yy5xwiIiIiIpJ7QvwOICIiIiJS2KjIFhERERHJZSqyRURERERymYpsEREREZFcpiJbRERERCSXqcgWEREREcllKrJFRMR3ZjbGzLaZ2ZJsLn+1mf1hZkvN7N28zicicqJM42SLiIjfzOw8IBEY75yLPc6y9YEpwPnOuV1mVsk5ty0/coqIZJeuZIuIiO+cc98DCZnbzOwUM/vCzOaZ2Q9mdpr31s3Aq865Xd66KrBFJOioyBYRkWA1CrjDOdccuBt4zWtvADQws5/MbLaZdfItoYjIUYT5HUBERORIZlYCOBt438wONUd6X8OA+kAcUB34wcxinXO78zmmiMhRqcgWEZFgFALsds41yeK9DcBs51wqsNbMlhMouufmYz4RkWNSdxEREQk6zrm9BAro7gAWcKb39idAe6+9AoHuI2v8yCkicjQqskVExHdm9h7wC3CqmW0wsz7AdUAfM1sILAW6eIt/Cew0sz+Ab4F7nHM7/cgtInI0GsJPRERERCSX6Uq2iIiIiEguK3Q3PlaoUMHVrl37hNfbv38/xYsXz/1AORAsWYIlByhLMOeA4MkSLDkg51nmzZu3wzlXMQ8iBa2cnLcLw7HOC8GSJVhyQPBkCZYcoCy5meOY52znXKF6NG/e3OXEt99+m6P18kKwZAmWHM4pS1aCJYdzwZMlWHI4l/MswG8uCM6l+fnIyXm7MBzrvBAsWYIlh3PBkyVYcjinLFnJi3O2uouIiIiIiOQyFdkiIiIiIrlMRbaIiIiISC4rdDc+ikjhlJ6eTkJCAqmpqUddplSpUmzatCkfUx3d8bKEh4dTrlw5QkND8zFVwXG8461jLSLBTkW2iBQICQkJREVFUaFCBcwsy2VSUlKoWrVqPifL2rGyOOdITEwkISGBihWL1EAi2Xa8461jLSLBTt1FRKRASE1NpUSJEkctsAsSM6NEiRLHvCpf1BWW461jLVJ0qcgWkQKjoBdcmRWm7yWvFJZ9VFi+DxE5Mcctss0syszmmNlCM1tqZo957Y+a2UYzW+A9Omda534zW2Vmy83sokztzc1ssffey+adecws0swme+2/mlntTOv0NrOV3qN3rn73IiIiIiJ5IDt9spOB851ziWYWDvxoZjO8915wzj2XeWEzawT0BBoDVYGvzayBcy4dGAn0BWYD04FOwAygD7DLOVfPzHoCw4AeZlYOeARoAThgnplNc87tOrlvW0REREQk7xz3SrY3oU2i9zLce7hjrNIFmOScS3bOrQVWAS3NrApQyjn3izdDznjgikzrjPOefwB08K5yXwTMdM4leIX1TAKFuYiIiIhI0MpWn2wzCzWzBcA2AkXvr95b/c1skZmNMbOyXls1YH2m1Td4bdW850e2H7aOcy4N2AOUP8a2RER80bVrVx566CHOPfdcKleuzNdff+13JMkj8+bNo3379v97vWTJEtq0aeNjIhEpSLI1hJ/X1aOJmZUBPjazWAJdP4YSuKo9FHgeuAnI6g4Pd4x2crjO/5hZXwLdUIiJiSE+Pv4Y303WEhMTc7ReXgiWLMGSA5QlmHNA/mQpVaoUKSkpALzy02ZW7Uj6xzLOOczW5mj79SpEcUfbKsddbsGCBTRs2JAJEybwxRdf8MYbb1CvXr1/LJeSksK6deuOua2EhARWrFiRo7xFyWOfLuWPTXsPa0tKSiIqanOOt9moaikeuazxMZdp2LDhYcfn4YcfZujQoTn+TBEpWk5onGzn3G4ziwc6Ze6LbWZvAp95LzcANTKtVh3Y5LVXz6I98zobzCwMKA0keO1xR6wTn0WuUcAogBYtWri4uLgjFzmu+Ph4crJeXgiWLMGSA5QlmHNA/mTZtGnT/8YiLrV4P1GJ/1wmUHhF5Wj7pUqVonbt2sdc5sCBA+zfv5+hQ4cSGhpKuXLlqF69OrVr12bIkCGHFWDr1q077vYiIiJo1qxZjvJK3itWrBhRUVHs3r2bNWvWsGvXLi644IJ/HGsRkawct8g2s4pAqldgRwMXAMPMrIpz7tBlhK7AEu/5NOBdMxtO4MbH+sAc51y6me0zs9bAr8D1wCuZ1ukN/AJcBXzjnHNm9iXwVKauKB2B+0/yexaRAu5oVyCzU9iejKVLl9K8efP/zdy3aNEiYmNj2bJlC2lpaXn2uUVdVsc7r4/1IY0aNeLPP/9k6NChPPHEEzrWIpJt2bmSXQUYZ2ahBPpwT3HOfWZmE8ysCYHuG+uAWwCcc0vNbArwB5AG9PO6mwDcBowFogmMKnJolJLRwAQzW0XgCnZPb1sJZjYUmOst97hzLiHn366ISM4tWbKEJk2a/O/1okWL6NKlC7///vth7VJ4NG7cmLfffhvnHG3btmXGjBk61iKSLcctsp1zi4CmWbT3OsY6TwJPZtH+GxCbRXsS0P0o2xoDjDleThGRvLZ48WJatWr1v9dLliwhNjaWF198kW7duvmYTPJK48aN6d27N3PnBq71LFiwQMdaRLLlhPpki4gUZcOHDz/s9Zo1awBYtWoV9evX9yOS5LFevXrRq9f/X1PSsRaR7FKRLSJykkaPHu13BMknOtYikl3ZGidbRERERESyT0W2iIiIiEguU5EtIiIiIpLLVGSLiIiIiOQyFdkiInJcZhZlZnPMbKGZLTWzx7JYxszsZTNbZWaLzEzTWYpIkaXRRUREJDuSgfOdc4lmFg78aGYznHOzMy1zMYFZfusDrYCR3lcRkSJHV7JFROS4XECi9zLce7gjFusCjPeWnQ2UMbMq+ZlTRCRY6Eq2iIhki5mFAvOAesCrzrlfj1ikGrA+0+sNXtvmLLbVF+gLEBMTQ3x8/GHvlypVipSUlKNmSUlJYd26dSf8PeSF7GRJSEhgxYoVeZ4lMTHxH/vSD8GSA4InS7DkAGXJrxwqskVEcsFjjz1GQkICZcqU4bHHHmPz5s3ceOONXH755cyePZvJkyf7HfGkOefSgSZmVgb42MxinXNLMi1iWa12lG2NAkYBtGjRwsXFxR32/qZNm6hatepRs6xbt47atWufSPxclfl4d+rUiQceeOCYxzoiIoJmzfK+i3p8fDxH7ks/BEsOCJ4swZIDlCW/cqi7iIjISdq4cSOpqamUKVOG2bMDXZSXLVvGFVdcwZ133klYWOG6nuGc2w3EA52OeGsDUCPT6+rApvxJlX+OPN6F+ViLSM7pbCAicgK6du1K48aN+e6771i5ciXvvPMO7777Li+99BLbt29n/fpAb4lly5bRp08fAMyyusBbsJhZRSDVObfbzKKBC4BhRyw2DehvZpMI3PC4xzn3j64iBcW8efO4++67+fbbbwFYsmQJN998Mw0bNjzseBe2Yy0iuUNFtojICViyZAlt27blhx9+4KOPPmLixInExsby3HPPsXPnTpo2bQrAX3/9RYMGDdixYweVK1f2OXWuqAKM8/plhwBTnHOfmdmtAM6514HpQGdgFXAAuNGvsLmhYcOGh/Wjfvjhhxk6dCgLFy487Hj/8MMPhe1Yi0guUJEtIgXPjPtgy+J/NFdOSoKoqJxts/LpcPHTx1zkwIED7NmzhzvvvBOAtLQ0ypQpw6BBg/6x7DPPPENISAgVKlTgueeey1mmIOKcWwQ0zaL99UzPHdAv1z88i+N9UscasnW8ixUrRlRUFLt372bNmjXs2rWLCy64gAsuuOCw5S677LJCdaxFJHeoT7aISDYtXbqU5s2bExoaCsCiRYuIjY31OZXkpUaNGvHnn38yZMgQnnjiCb/jiEgBoivZAKu+psHyUZD0FYRHQ2QJKFkVSlWF0tWhTC0I0e8jIkHjKFcgt+TxiBNLliyhSZMm/3u9aNEiunTpAsCQIUMYOnRonn12kZbF8c7rY31I48aNefvtt3HO0bZtW0DHWkSyR0U2sH/zCspu/5WMnT9hqQcwl3H4AuHFoVJDqBwLNc+G2m0DxbeIFCmLFy+mVav/n8BwyZIlxMbGsmXLFtLS0nxMJnmlcePG9O7dm7lz5wLoWIsURumplNi3BojL1c2qyAbGp1/EsMRa3itHNMnUj9pHw+L7ODUygUahG6iTtI4Kiz8idN7YwGJl68CpnaHhZVCjla50ixQBw4cPP+z1mjVrgMD4qpmvcEvh0atXL3r16vW/17///ruOtUhhsmkBTO1Pkx2rocMVUKxcrm1aRTZwQcNK7Nq0ljr1GnAwJZ39yWlsT0xm854kPtpzkOe27udASjohZBAbup6uZdfSzi2l9pw3CZn9KpSIgdgroVlvqHSa39+OiOSzBQsW0K1bN79jSD7QsRYpJFKT4Lth8NNLULwCf542gNhcLLBBRTYA9WNKcnbVMOJa1szy/YwMx98JB/hzy14WrK/P1DU7eWJje6Iz9nNx5GL+xe+c/uubhMx+LXBVu8VNgaI7NDyfvxMR8cOqVauoX7++3zEkH+hYixQC6+fA1H6wYwU0+RcH2z/OT7MWkNu3sR+3yDazKOB7INJb/gPn3CNmVg6YDNQG1gFXO+d2eevcD/QB0oH/OOe+9NqbA2OBaALjqQ5wzjkziwTGA82BnUAP59w6b53ewENenCecc+NO+rs+QSEhRu0KxaldoTidYqsAkJicxuzVO/nqjwbc8Mc52MFruSbyJ27Y9h0VP74FN2sodnZ/aHY9RBTP78giko9Gjx7tdwTJJzrWIgVYyn745gmYPTJwb92/PuJHdyYPvLGIrXuS6NEphTLFInLt47JzJTsZON85l2hm4cCPZjYD6AbMcs49bWb3AfcBg82sEdATaAxUBb42swbOuXRgJNAXmE2gyO4EzCBQkO9yztUzs54EZhHr4RXyjwAtAAfMM7Nph4p5P5WIDOOCRjFc0CiGtPQMfl2bwIfzGnLO4s60yZjPoAMzOP2L+3DfPYOdcye07AvhJzGmq4iIiIjkzNrvYdodsGsdnHUzu85+gKEz/+aj+b9St0Jx7mwelasFNmRjnGwXkOi9DPceDugCHLqqPA64wnveBZjknEt2zq0lMPNXSzOrApRyzv3iTVgw/oh1Dm3rA6CDBeamvQiY6ZxL8ArrmQQK86ASFhpC23oVGN6jCXMevJAOl/2LO4s9RbfkR/k1uRbMHELGK83g93cgI93vuCIiIiJFQ9Je+HQgjLsMLAR3w+d8UvVOOoyYx7QFm+jfvh7TB5zLaeVCc/2js9Un25tGdx5QD3jVOfermcU45zYDOOc2m1klb/FqBK5UH7LBa0v1nh/Zfmid9d620sxsD1A+c3sW62TO15fAFXJiYmKIj4/Pzrd1mMTExBytl5UawANNHQu3n87Da0+j/J7FPLBnEqdP7cfeWc+zqv6t7C19ar5kORnBkgOUJZhzQP5kKVWqFCkpKcdcJiUlhXXr1uVpjuzKTpaEhITDpu0WEZFctHImfDoA9m2Gs+9gw5kDefDzNXy3YgFNapTh6StP57TKpfLs47NVZHtdPZqYWRngYzM7Vt9wy2oTx2jP6TqZ840CRgG0aNHCxcXFHSNe1uLj48nJesdyPnAnMO+vVjwx4xwq/j2dR/ZPpNnv9+Ka9cYueDTLoWLyIktOBEsOUJZgzgH5k2XTpk1UrVr1mMusy6cJSrIjO1kiIiJo1qxZ/gQSESkqDiTAF/fDoklQsSFpV41j7F/leX7Eb4QYPHpZI3q1qU1oSFZlZu45odFFnHO7zSyeQJeNrWZWxbuKXQXY5i22gcDF3EOqA5u89upZtGdeZ4OZhQGlgQSvPe6IdeJPJHMwaF6rHJNuacP3K+tz24xz6bh9HH3mT8AtnUbYZcMhVsNBiWSHc45AT7KCL9BrTo6lsBxvHWuRfPTHVPj8bjiYAOfdyx/1b2bwJytYvHEZ559WiaFXxFKtTHS+RDlun2wzq+hdwcbMooELgD+BaUBvb7HewFTv+TSgp5lFmlkdoD4wx+tass/MWnv9ra8/Yp1D27oK+Mbrt/0l0NHMyppZWaCj11bgmBntGlRkyh0dKdf1Ga4NfZYlB8vBBzeS8t71sH+n3xFFglp4eDiJiYmFomBxzpGYmEh4uIb5PJrCcrx1rEXySeI2mHJ94FGyMkk3fsN/k7tx2cjf2LznIK9c05TRvVvkW4EN2buSXQUY5/XLDgGmOOc+M7NfgClm1gf4G+gO4JxbamZTgD+ANKCf190E4Db+fwi/Gd4DYDQwwcxWEbiC3dPbVoKZDQXmess97pxLOJlv2G8hIcZVzatzYaPevPTVWXw9dwQDln9E8ss/EXnla9DgIr8jigSlcuXKkZCQwL59+466TEJCAhERuXt3eE4dL0t4eDjlyuXuxAeFyfGOt461iADgHCyaAl8MDgzR1+Fhfo65lvsn/clfOw9wdYvqPNC5Ya6PHJIdxy2ynXOLgKZZtO8EOhxlnSeBJ7No/w3+Oda3cy4Jr0jP4r0xwJjj5SxoSkeH83CXM1l61nMMfO8c+u95jobvXk1Si9uwYlnuVpEiLTQ0lIoVKx5zmRUrVgRNH+dgylIQHe94B9P+DaYsIkXKno3w2Z2w8kuo3pI9HV/g8dnpfPj5fGqXL8a7N7fi7FMq+BZPMz76rHHV0rw4sBdvzGrOvB8e5V+/jaRB9NfQtAGUreV3PBEREZHg4hzMHwdfDYGMNNxF/2Va1KU8Pm45ew6mcnvcKfynQ32iwnN/WL4Tcdw+2ZL3wkND6N8xlma3jeHJ4vdR/MAGDo5oS+ofn/sdTURERCR47FoH47sEhuarciabr/uGG/9swYDJi6leNppP7ziHezud5nuBDSqyg0qjqqW4+857ebrCMFallid0ynUkfPF04Dc2ERERkaIqIwNmvw6vtYGN88m45EVGn/IyHcb8xZy1CTx8aSM+ur0tDavk3bjXJ0rdRYJMZFgoHU+vxdZyU1n/0e10nv1fNm5YSLXrR0NEMb/jiYiIiOSvHSthan9YPxvqd2TFWUO556sdLNywjPanVmToFbFULxt8NZKuZAepC86ozZkDPmBC8Ruosn4GW16MI33X+uOvKCIiIlIYpKfBD8NhZFvY/icpl49kWLnHuXjsGjbuPsjL1zRlzA1nBWWBDSqyg1q1ssXocecLTKj7NMX3/83eEeexb93vfscSERERyVtblsBbHWDWY9DgIuZe8gUdZ1Vh5Hdr6Na0Gl/f1Y7Lz6wa1BNWqcgOchFhIfTufSvfnzuRpDQIGXsxG+ZN9zuWiIiISO5LS4Fvn4JR7WDvRvZ3Gc09IXfTfeIaHPDuv1vxbPczfRn3+kSpyC4gLrmgA9t7fM4mKhIz7V+s+OotvyOJiIiI5J6N8wLF9XfDcLHdmNFuKu0+L81Hv2/ktrhT+HLgeZxdz79xr0+UiuwC5IxGjShx29csCWtMg58HsXjSIxp5RETyhZnVMLNvzWyZmS01swFZLBNnZnvMbIH3eNiPrCJSwKQeDIx5/dYFcHA3O7pMoM/evtz20TqqlYnm0/7nMDhIhuU7ERpdpICpEhND8bu+4KcR19H2zxeZOzqBFje9hIXo9yURyVNpwCDn3HwzKwnMM7OZzrk/jljuB+fcpT7kE5ECqPTupTDyLkhYTUazG5hYqg///XgTsJMhlzbihrNrExoSvP2uj0WVWQFUqnhxzrrrA34u24WzNozn5xE3kZqW5ncsESnEnHObnXPzvef7gGVANX9TiUiBlZwI0++h6YIHICONvy+bRNf1VzPki/W0rFOOr+48jz7n1CmwBTboSnaBFREeRps7xvLb6Dtou/EdfnrxAM3vmEBUZKTf0USkkDOz2kBT4Ncs3m5jZguBTcDdzrmlR9lGX6AvQExMDPHx8SeUITEx8YTXySvKErw5IHiyBEsO8D9L2YQFnLr8VSKTt7O6YideC72OqR9kUCx8N7ecEUnrKvtZtXAOq/IxU17sExXZBZiFhNDi3yNY8E5J2q4eyewXriT2P1MoUSw4x4sUkYLPzEoAHwIDnXN7j3h7PlDLOZdoZp2BT4D6WW3HOTcKGAXQokULFxcXd0I54uPjOdF18oqyBG8OCJ4swZIDfMxycDd89SAsegfK12dx3Kv8e2YaWw84rmxWnYcuaUjZ4v6MGpIX+0TdRQo6M5r0eprFje+hddIPLH7xSnbv2+93KhEphMwsnECBPdE599GR7zvn9jrnEr3n04FwMys4QwGISN75czq81hoWvEdSqwE8GDOSy6amkuHgnT6teP7qM30rsPOKiuxC4vTuD7GsyQO0SfmZxS93Z9ueRL8jiUghYoEZH0YDy5xzw4+yTGVvOcysJYH/Y3bmX0oRCTr7d8AHfWDSNbhi5fghbjLnzDuPSQu2c0u7ujxxTjTn1C+cv4uru0gh0vCKwawNgXPnP8WsV3oS238SMWVK+B1LRAqHtkAvYLGZLfDaHgBqAjjnXgeuAm4zszTgINDTOY0zKlIkOQdLP4Lp90DSXva2uZe7N7Xnqxm7OL1aacbeeBax1UoTH7/V76R5RkV2IVPn8sFsCEmnw2/DmDniWs7o/x4xZYr7HUtECjjn3I/AMW/zd86NAEbkTyIRCVr7tsBnd8Hyz3FVmzG11gM8+FM6GW4vD13SkBvOrk1YaOHvTKEiuxCqfukDbMhI58L5z/HliOto2n8ilVRoi4iISF5yDha8C1/eD2nJbGv9ELeuasX8b/dxXoOKPHlFLDXKFZ3BGQr/rxFFVPXLh7CxyZ1clPYtv7zah217D/odSURERAqr3X/DO1fC1NvJqNiIt2Lf4ezvG7NuVzIv9mjCuBvPKlIFNuhKdqFW7YpH2Zx2gC5L3uCdEXfQ8Y4RVCoZ5XcsERERKSwyMmDeGJj5CDjH2paP8e+lZ7B65UG6NavGQ5c0olwhGzUku457JdvMapjZt2a2zMyWmtkAr/1RM9toZgu8R+dM69xvZqvMbLmZXZSpvbmZLfbeeznTXeiRZjbZa//Vm+jg0Dq9zWyl9+idq999EVDlymFsq9+Df6W8z0cj7mfX/hS/I4mIiEhhsHM1jLsUPh9EatUWDDvlbdp/X58UBxP6tGT41U2KbIEN2buSnQYMcs7NN7OSwDwzm+m994Jz7rnMC5tZI6An0BioCnxtZg2cc+nASAIzfM0GpgOdgBlAH2CXc66emfUEhgE9zKwc8AjQAnDeZ09zzu06uW+7CDGj0jUj2Tl2L7f+PYaXRpbipv4PUTIq3O9kIiIiUhBlpMPs1+CbJ3Gh4Sxu/iR9Fp5KwoFUbjmvLgMvaEB0RKjfKX133CvZzrnNzrn53vN9wDKg2jFW6QJMcs4lO+fWAquAlmZWBSjlnPvFG9JpPHBFpnXGec8/ADp4V7kvAmY65xK8wnomgcJcTkRIKOWvH0dCTFv673uJN954iaTUdL9TiYiISEGzbRmM7ghfPURSzfO4u9KbXP5THWJKRzG1X1vu79xQBbbnhPpke904mgK/Ehgztb+ZXQ/8RuBq9y4CBfjsTKtt8NpSvedHtuN9XQ/gnEszsz1A+cztWayTOVdfAlfIiYmJydHc83kxZ31O5VWW0Pp3UHfPdu5I+C9PDHe0O6sFYSFHH5GrKOyTnAiWLMGSA4InS7DkgODKIiJy0tJT4ccX4bthuKhS/HDG09y+sA5pGY4HOzfkxrZFY1i+E5HtItvMShCYTnegc26vmY0EhhLoxjEUeB64iazHUXXHaCeH6/x/g3OjgFEALVq0cDmZez4v5qzPqTzNcnYrdr92Affue55X/3qFe3pfRehRCu0is09OULBkCZYcEDxZgiUHBFcWEZGTsmkBTO0PWxezt14X/rO7J/FzHOfWL8tTXU8vcqOGZFe2fuUws3ACBfZE59xHAM65rc65dOdcBvAm0NJbfANQI9Pq1YFNXnv1LNoPW8fMwoDSQMIxtiU5VawcZW6ehkWV5IZ19/Dc+7PQhGwiIiLyD6lJ8PVj8Ob5uMRtTD31GZov68nChDBe6HEm429qqQL7GLIzuogBo4FlzrnhmdqrZFqsK7DEez4N6OmNGFIHqA/Mcc5tBvaZWWtvm9cDUzOtc2jkkKuAb7x+218CHc2srJmVBTp6bXIySlenxE2fUDYshSuWDmDkF/P9TiQiIiLBZP0ceONc+HE420/pRhcbzoCF1bn0jKp8fVc7ujatjjdInBxFdrqLtAV6AYvNbIHX9gBwjZk1IdB9Yx1wC4BzbqmZTQH+IDAyST9vZBGA24CxQDSBUUVmeO2jgQlmtorAFeye3rYSzGwoMNdb7nHnXEJOvlE5QkxjIq57l1MmXEnCz/2ZXGYCPdrU8zuViIiI+CllP3zzBMweSUapaoytM5zHl1Smetloxt10Ou0aVPQ7YYFx3CLbOfcjWfeNnn6MdZ4Ensyi/TcgNov2JKD7UbY1BhhzvJxy4qxuHHR5lTaf3MKO6f2ZVfptOjSqctz1REREpBBa+z1MuwN2reOvutfS++/O/L09hJvPrcOdFzagWITmMDwRug20iAtr0pOUuCFcFvoLaybdw+9/awhyERGRIiVpL3w6EMZdRpoznqnyAu3+uJTipcowtd85PHhJIxXYOaA9JkS0G8TBXeu5eeFYnnq7MmX6PUadCsX9jiUiIiJ5beVM+HQAbt9m/qjdm95rLyBxZzj3X9yAPufU0bB8J0F7TsCM6C7DOVDrAu51oxnx5hts35fsdyoRERHJKwcS4KNbYOJVJIeW4L6yw7nkz4s4rUYMXw48j1vanaIC+yTpSrYEhIRS7NqxHHz9Ah5NeIb7Rlfm8lj1zxYRESlsKmz/GV69GXcwgV+r96HP2naER0bzfPdGdGtWTaOG5BL9iiL/L7Ik0b3fJyK6OIN3PcykhdvJyNAY2iIiIoVC4jaYcj2xS4exP7Iit0Q/T89VHbjw9JrMuqsdVzbXsHy5SUW2HK5MTSL/NZmqoXu5fc+zPDdjkd+JRERE5GQ4Bwsnw6stcctn8EHxazhz02CWptdk7I1n8WLPppQvEel3ykJHRbb8U/UWhHZ7nbNCVlBv9gNMnvOX34lEREQkJ/ZshHd7wMd92RVdm54hz3LPzsu44Zz6zLzrPOJOreR3wkJLRbZkyWK7sbrWNXQL/ZEN057g51U7/I4kIiIi2eUczBsLr7XGrf2eyeVvp/mmQewrcQoPt4nioUs1LF9eU5EtR7W+dg9SGl/FoLApfPjOCFZtS/Q7koiIiBzPrnUwvgt8OoAtxU+lc+ozPLz1PO69uBFT+7elTulQvxMWCSqy5ejMiLjiVZKrnMWTvMozo98lYX+K36lEREQkKxkZMPt1eK0NGRvm8UapO2i9aSDlajTgqzvP49Z2pxCuYfnyjfa0HFt4FJHXvYeVqMjQpKcYPPYrktPS/U4lIvnMzGqY2bdmtszMlprZgCyWMTN72cxWmdkiM2vmR1aRImnHSnj7YvhiMGtLNCXuwNOMTDyPZ686k3f6tKJWeU0yl99UZMvxlahI5L8mUz4sidu2PsKQD37DOQ3tJ1LEpAGDnHMNgdZAPzNrdMQyFwP1vUdfYGT+RhQpgtLT4IfhMLItaVuX8XTUnbTffDtNYmP5+q52dG9RQ8Py+URFtmRP5dMJ6/Y6zUJWcdaSJ3jr+zV+JxKRfOSc2+ycm+893wcsA6odsVgXYLwLmA2UMTPNaiWSV7Ysgbc6wKzHWFKiDW32/pdPOY+3b2zJy9c0pYKG5fOVimzJvsZX4M69h+5h37N55gvEL9/mdyIR8YGZ1QaaAr8e8VY1YH2m1xv4ZyEuIicrLQW+fQpGtSM5YT33h93N5dv6clnbpnx153m017B8QUFjt8gJsfYPkLZ1KQ+umMit79WiRr/bOKViCb9jiUg+MbMSwIfAQOfc3iPfzmKVLPuWmVlfAl1KiImJIT4+/oRyJCYmnvA6eUVZgjcHBE+W3MpRcu9KTl3+CiX2/8V34ecyYM+/KFmyNA+1jqBuyW3M/eX4F8CCZZ9A8GTJixwqsuXEhIQQduUoUkd14LmdL9Lv7eq8esdVlI4O9zuZiOQxMwsnUGBPdM59lMUiG4AamV5XBzZltS3n3ChgFECLFi1cXFzcCWWJj4/nRNfJK8oSvDkgeLKcdI7Ug/DtU7jfR3AwogK3cx+zDjZhQKf63Hxu3RMaNSRY9gkET5a8yKHuInLiIksSft1kikeG8UjiUO6d+CPpGboRUqQws8CdU6OBZc654UdZbBpwvTfKSGtgj3Nuc76FFCms/voZRraFn1/m66iLaLXnKXZVbc8XA8/j9rh6GpYvSOlKtuRMuTqE9RjPKRO6ctVfjzNsRiUeuKSx36lEJO+0BXoBi81sgdf2AFATwDn3OjAd6AysAg4AN+Z/TJFCJDkRZj0Gc0axJ6oad6Q9xMKkMxlyVUO6N6+uUUOCnIpsybm67Qjp9F8unHEvf/7yDB9WfpIrm1f3O5WI5AHn3I9k3ec68zIO6Jc/iUQKudXfwLQBuD3r+TjiMh7c3ZULzqzL15c2omJJjRpSEKjIlpPTsi8Zmxdzx4IJDPikFnUr3knTmmX9TiUiIlIwHdwNXz0Iv7/D9sia3JryCJsjz+DVG2I5/7QYv9PJCThuJ56jzfJlZuXMbKaZrfS+ls20zv3ejF/LzeyiTO3NzWyx997LXh8/zCzSzCZ77b96w0MdWqe39xkrzax3rn73cvLMCLn0edKqnsWw0JE8N/4DtuxJ8juViIhIwfPndHitNRkL3mN8aDfO2fs4p7fuyFd3tVOBXQBlp6f80Wb5ug+Y5ZyrD8zyXuO91xNoDHQCXjOzUG9bIwkM2XRoRrBOXnsfYJdzrh7wAjDM21Y54BGgFdASeCRzMS9BIiySsGvfJbR4eZ5NfZp7xn9DUqqmXhcREcmW/Tvggz4w6Ro2JEdzedJjTCxxI5Nui+PRyxtTIlIdDwqi4xbZx5jlqwswzltsHHCF97wLMMk5l+ycW0vgBpiW3qxfpZxzv3j99sYfsc6hbX0AdPCucl8EzHTOJTjndgEz+f/CXIJJiUqEX/suMaH76Lf9cR79ZIGmXhcRETkW52DJh7hXW5K+9BNepQcX7X+cThd24tM7zlH3ywLuhMZ8OWKWr5hDQzN5Xw9NL3S0Gb+qec+PbD9sHedcGrAHKH+MbUkwqtaM0C6v0DpkGactfJp3fv3b70QiIiLBad8WmHQdfHATa1LLc3HSk3xX+UamDmhP//PrExGmYfkKumz//eHIWb6OMWzM0Wb8OtZMYDlZJ3O2k5o5DIJnxiEIniw5yxFD3epduGHDVO779BVGbb6YBmVDj79anmTJG8GSJVhyQPBkCZYcEFxZRCSIOAcL3sV9eT/pKUkMz7iOd1Mu5Z4rGnPNWTUJCdGwfIVFtorso8zytdXMqjjnNntdQQ7N43m0Gb82eM+PbM+8zgYzCwNKAwlee9wR68Qfme9kZw6D4JlxCIInS45znHsOqROuYui6t7llSW0u+08fqpSO9idLHgiWLMGSA4InS7DkgODKIiJBYvff8OlAWD2LpWGN6X/wJuo3bMIXXWKpXDrK73SSy7IzusjRZvmaBhwa7aM3MDVTe09vxJA6BG5wnON1KdlnZq29bV5/xDqHtnUV8I3Xb/tLoKOZlfVueOzotUkwCw0jvMfbuFLVGJb2LA+M+4rkNN0IKSIiRVRGBlU3zsC91oaUtT/zcNqN3MSjDL72Ekb1aq4Cu5DKToefQ7N8nW9mC7xHZ+Bp4EIzWwlc6L3GObcUmAL8AXwB9HPOHaqwbgPeInAz5Gpghtc+GihvZquAu/BGKnHOJQBDgbne43GvTYJddFki/jWZsuEpDNjxKI9/NF83QoqISNGzczWMu5QGK1/nt7RTOP/g0yQ3uYmZd7Xn4tOraNbGQuy43UWOM8tXh6Os8yTwZBbtvwGxWbQnAd2Psq0xwJjj5ZQgVKkhYVe+SZPJ17Fq8eO8U+MVerWp7XcqERGRvJeRDrNfw33zJEkZoTyc2pe5pS/mmevO4OxTKvidTvKBbl2VvNXwUjLa3cdVod/z1+fPMXed/hAhIiKF3LZluNEd4auH+D49lvOThrG/5oV8cWc7FdhFiEY3lzwX0m4wqZsWcf/Kidwxvg41BvRT/zMRESl80lPhxxdx3w0jkWI8kNKftZUv4s0rz2THyt+JCj/50bak4NCVbMl7ISGEXzWKtLL1eSp9OI+M+0w3QoqISOGyaQFuVBx8+wQz0s+iY+qznN7pJj7pdw6x1Ur7nU58oCJb8kdkSSJ7TaJYRCh37niUJz6cqxshRUSk4EtNgq8fw715Pru3b6Rvyp1MrPEIkwZeSt/zTiEsVKVWUaUjL/mnXF0ieo6lQchG2iwZwsTZ6/xOJCIiknPr5+BePxd+HM6H6edyWcZwLuzWh3f6tKJW+eJ+pxOfqciW/HXK+XDh43QOncOO6U/ym26EFBGRgiZlP3xxP250R7Yl7KJXyn3En/YIHw+6hO4tamhYPgF046P4IOTs/qRsWsjApe8zaHxdqg+4SzdCiohIwbD2ezKm3kHI7nVMSL+QcdE3cH/3s7igUYzfySTI6Eq25D8zIq54haSKZ/J4+ks8MfZj3QgpIiLBLWlvYEr0cZexcU8yVycPYVWLR/nkrotUYEuWVGSLP8KjifrXe4RFlWDQzkd4+qPZuhFSRESC04qvSH+1FRnzxvFG2iXcWuJl7r3lJh7vEkvJqHC/00mQUpEt/ildjcjrJlIzdCdxi+/j3dlr/U4kIiLy/w4k4D7qC+92Z+2+ULqnPcb+do/y0YAOtKhdzu90EuRUZIu/arbGLnmedqGL2D99iGaEFBGR4PDHVNJHtCR90Qe8lNaV+yuO4Kn+N3LXhQ2IDNOkMnJ8KrLFdyEtbiC56Y30Df2MT8a/yKbdB/2OJCIiRVXiNtzk62HK9SzfX5zuGU9R+uJHmHRbO06tXNLvdFKAaHQRCQqRlzzDwa1/MGTja9z3dh2e7n+9pp8VEZH84xwsmkL69HvJSE7khdQerDjlBkZ0a0q1MtF+p5MCSFeyJTiERRB97UQoXoF7dg/lySnf60ZIkSBjZmPMbJuZLTnK+3FmtsfMFniPh/M7o0iO7NlI+sSr4eO+LDxYkZ4hz3Fq90d488Y2KrAlx1RkS/AoUZGof02iUmgily6/j7e/X+F3IhE53Fig03GW+cE518R7PJ4PmURyzjmYN5a0Ea1IWfUdj6X2YmLjN3hz0HV0aVJNk8rISVGRLcGlahPCrniVViF/Evn1/Xy/YrvfiUTE45z7HtDdyVI47FpH2tjL4dMBzEmqQe/IF2nf+xGe79GccsUj/E4nhYD6ZEvQsTO6k7JpIdfNfoWh7z5DrTse9TuSiGRfGzNbCGwC7nbOLc1qITPrC/QFiImJIT4+/oQ+JDEx8YTXySvKErw5IIssLoNqGz+n9uoJJLsQnkrtw85qHenTIJKMTUuJ35RPOXykLPmTQ0W2BKWIjo9xcNMS7vt7NPeOqUvHpqf5HUlEjm8+UMs5l2hmnYFPgPpZLeicGwWMAmjRooWLi4s7oQ+Kj4/nRNfJK8oSvDngiCzbV5D68e2Eb5rLN+lNeKv0HdxzdQea1iybvzl8piz5k0PdRSQ4hYQSfc1Y0krW4MHEp/h44UYyMnQjpEgwc87tdc4les+nA+FmVsHnWCKQnob7YTjpI9tyYNMy7k67nYXnjuLtgV3zpcCWoklFtgSv6DJEXz+FUmHp9N/zHK/NzHJAAxEJEmZW2bw7xcysJYH/Y3b6m0qKuuKJ60gddT426zG+TG3CHeVG8u/+93Nnx1M1qYzkKXUXkeBW8VTCu79F7KRrWf3jYL6qNpqOsVX8TiVSJJnZe0AcUMHMNgCPAOEAzrnXgauA28wsDTgI9HQai1P8kpaC+/5Zmv32PLspwWPpAzn9wt6MaVubsFBdY5S8d9yfsqzGRTWzR81sY6axUDtneu9+M1tlZsvN7KJM7c3NbLH33suZrnZEmtlkr/1XM6udaZ3eZrbSe/TOte9aChQ7rTOra11D19CfWDDlCVZs3ed3JJEiyTl3jXOuinMu3DlX3Tk32jn3uldg45wb4Zxr7Jw70znX2jn3s9+ZpYjaOI+Ukedi3z/D1PQ23BvzJncNuJebz6urAlvyTXZ+0saS9bioL2QaC3U6gJk1AnoCjb11XjOzQ3+LGUngTvL63uPQNvsAu5xz9YAXgGHetsoRuErSCmgJPGJm6jhVRG2sfTUH61/KoJB3eXPMm+w5kOp3JBERCTapB8n4cggZb15Awo5t3O7uY2H9gbx5S0dqVyjudzopYo5bZJ/guKhdgEnOuWTn3FpgFdDSzKoApZxzv3h/OhwPXJFpnXHe8w+ADt5V7ouAmc65BOfcLmAmx58EQQorM6KveoPksvV5KOk5Hh//Kem6EVJERA7562dSRrQh5JeXmZQWx9BaYxhy10Da1wwnJESTykj+O5k+2f3N7HrgN2CQVwhXA2ZnWmaD15bqPT+yHe/regDnXJqZ7QHKZ27PYp3DnOx4qxA84zRC8GQJlhzgZfnlN6IaDOKMuXdxy6aHuGuk0bVx/v9xI1j2S7DkgODJEiw5ILiyiBRqyYmkz3yU0N/eZKurxNCQh7nkqp6MOLMqZsZyv/NJkZXTInskMBRw3tfngZuArH5VdMdoJ4frHN54kuOtQvCM0wjBkyVYcsARWU6tTL3xXbliywtsbT6WHq3r+pfFR8GSA4InS7DkgODKIlJorf6GlI/vICxxI2PSOrG04QD+26U55UtE+p1MJGdD+Dnntjrn0p1zGcCbBPpMQ+Bqc41Mi1YnMOvXBu/5ke2HrWNmYUBpAt1TjrYtKerqtsN1fpb2oQtJ/nwwP6/e4XciERHJTwd3k/bx7TChK+v3ptM37AlqXPsyz193tgpsCRo5KrK9PtaHdAUOjTwyDejpjRhSh8ANjnOcc5uBfWbW2utvfT0wNdM6h0YOuQr4xuu3/SXQ0czKejc8dvTaRAht2YfklrdzfehX/DDhCdbu2O93JBERyQ9/Tifl5bOwhe/xWtrljD1jAsPvvoULG8X4nUzkMMftLnKUcVHjzKwJge4b64BbAJxzS81sCvAHkAb0c86le5u6jcBIJdHADO8BMBqYYGarCFzB7ultK8HMhgJzveUed85l9wZMKQIiOz3BgR1ruHvNWIaMrsbgOwZSuli437FERCQv7N9B6uf3EP7HR6zOqMnzxZ7hpqu6cns9TSoqwem4RbZz7posmkcfY/kngSezaP8NiM2iPQnofpRtjQHGHC+jFFEhoRTrOYbENzry4I5neWpsVR695RrCNQaqiEjh4Rws/YiUTwdhyfsYnnYVSS3/wyudYomO0IyNErxUjUjBFlGcEr0/wKLL0H/rQ7zw4Xd+JxIRkdyybwspE6+BD27ij4Nl6VfiBdr3fY4HLj9TBbYEPRXZUvCVqkKxGz6iXGgSnZcMZOL3f/idSEREToZzuN/fIeXls8hYNYun068j/pyJjBj4L5rW1Lx0UjCoyJbCoXIsYT3G0ijkb2Jm9uO7P7f4nUhERHJi998kj+2KTe3H78lVGVh2BFf0e5qBHRsREaayRQoO/bRKoRF66kWkdvwvF4TO5+/37mT5ln1+RxIRkezKyMDNeYvUV1qR/tfPPJZxE4s6TOTVO67mtMql/E4ncsJOZsZHkaATefatJG5bRa8Fb/Lymw9RZsAwYkpF+R1LRESOZedqkj68nahNs/kl/XTerzKIQVd3pHaF4n4nE8kxXcmWQqfE5cPYU6cz/0kfy9tvPEdicprfkUREJCsZ6WT89DJpr55NysZFDHG3seGSibx06xUqsKXAU5EthU9IKKWvfZvdFc/irsThvDJ6NKnpGX6nEhGRzLYt4+Dr5xMycwjfpjbmsRpjuP2uR7i2dS1CQszvdCInTUW2FE7hUZS56X0OlKxDv62P8Op7nxCYSFRERHyVnkrat8+QPvIcDm5dzX02kAPdxvNcn05UKR3tdzqRXKMiWwqv6LKUuXkqFlmSnivvYtyMH/1OJCJStG1awMFXzyPsuyf5PO0snqs3gXsGPUCXptUx09VrKVxUZEvhVro6Jfp8QunQVNrOvoXPfl3idyIRkaInNYnUrx4lY9T57Nu5iXvDBlPsmnE81as95UtE+p1OJE9odBEp9CymMWHXTaL2hK7s+fwmZpf+iNanVfc7lohI0bB+Dgfev5Vie1czJa0dy84YzEOXt6RUVLjfyUTylK5kS5EQfsq5pF4+kmYhKzj4Xm/+2LDT70giIoVbyn6SP7uXjNEdSdizh0GRj1D9xjE8cnVbFdhSJOhKthQZxZp1Z/e+HbT/9j6mj76B4v3epVaFkn7HEhEpdMrsWsTBl/oTvX8949MvZOtZg3miUzOiI0L9jiaSb3QlW4qUMu1uY2erwXR23zPv9b5s23vQ70giIoVH0l6SPv4PTRYOYcu+FAYV/y9n9H2Ley4/SwW2FDm6ki1FTvlO97P1wG66LX6D914bQOcBr1I6Wn+6FBE5GW7FlyR9/B8iDm7jzbRLSD73Pv7bIZaIMF3Pk6JJP/lS9JgR020Ym0/pwTVJk5n22n0kpab7nUok6JnZGDPbZmZZDtNjAS+b2SozW2RmzfI7o/jgQAIHJ/8be/dq/t4fxuAyzxPd6t/0v+gMFdhSpOmnX4omM6pcN5KN1TrRa99bTH79CdI0K6TI8YwFOh3j/YuB+t6jLzAyHzKJj9zSTzj4YgvC/viIVzO68VOHDxn2nxupXlLlhYj+FUjRFRJKtRsnsLHCOfTa8QLvvf0iGRmaFVLkaJxz3wMJx1ikCzDeBcwGyphZlfxJJ/kqcRv737kOe783q5NK8lClV7h0wAhuancaoZoSXQRQn2wp6sIiqNb3fTaM6EzP9U8weUIkPa+/XTOPieRMNWB9ptcbvLbNRy5oZn0JXO0mJiaG+Pj4E/qgxMTEE14nrxSpLM5RcWs8tVe8RXh6EsMzepBYvyudakSxdvFc1uZXjhMQLFmCJQcoS37lUJEtElGMardPZeMrnblqzRA+eDeMq67tq0Jb5MRl9Y8myz8POedGAaMAWrRo4eLi4k7og+Lj4znRdfJKkcmyZyOJH95Bib9nMS+jPp/UuI9+V19K5dJR+ZvjBAVLlmDJAcqSXzmO210kqxtdzKycmc00s5Xe17KZ3rvfu+lluZldlKm9uZkt9t572bwKxswizWyy1/6rmdXOtE5v7zNWmlnvXPuuRY5gUaWp1v9zthY/lS4r7ueTKWP8jiRSEG0AamR6XR3Y5FMWyS3OkTbnbZJfPouQv37kWbuRjd0+5vF/X5llgS0iAdnpkz2Wf97och8wyzlXH5jlvcbMGgE9gcbeOq+Z2aGBMUcS+NPgoZtiDm2zD7DLOVcPeAEY5m2rHPAI0ApoCTySuZgXyW0WXYaq/aeztVh9Ov9xL59+MM7vSCIFzTTgem+UkdbAHufcP7qKSAGyax373ryEsOkDmZdSi+dOeZubBg3j8iY19Nc+keM4bpF9lBtdugCHKpBxwBWZ2ic555Kdc2uBVUBL78aXUs65X5xzDhh/xDqHtvUB0MG7yn0RMNM5l+Cc2wXM5Nh3tYuctJBiZana/wu2RdWh4+JBfDF1ot+RRIKGmb0H/AKcamYbzKyPmd1qZrd6i0wH1hA4978J3O5TVDlZGRmk/PQaKa+0go3z+W/orRzo+TEPX38J5UtE+p1OpEDIaZ/smENXJ5xzm82sktdeDZidablDN72kes+PbD+0znpvW2lmtgcoz9FvoPmHk72BBoKn4z0ET5ZgyQH5n8WaDSHt1yHEzR/A2M1bqH1ac9+yHE2w5IDgyRIsOSC4suQW59w1x3nfAf3yKY7kle0r2DvlVkptn8c36U2YGzuEfl3aUSpKk3aJnIjcvvHxaDe9HOtmmJysc3jjSd5AA8HT8R6CJ0uw5AB/sqS0OZutIzpyzeb/8mPV4XS4/HrfsmQlWHJA8GQJlhwQXFlEsiU9jaTvXyT0+2FkZITzZOQA2l99B4PrVfQ7mUiBlNNxsrceGvvU+7rNaz/aTS8bvOdHth+2jpmFAaUJdE/RDTTiq4hSFanU/ys2Rdbl3HkDmfXhm35HEhHJG1uWsHdEO6K+G8rXaU0Y13QKd939CGerwBbJsZwW2dOAQ6N99AamZmrv6Y0YUofADY5zvK4l+8ystdff+voj1jm0rauAb7w/OX4JdDSzst4Njx29NpF8E1mqItX+8yXro08jbtE9fDP5Zb8jiYjknrQUDnz5OOmvn0dywt8MLXYfVfq+z4ArziE6IvT464vIUR23u4h3o0scUMHMNhAY8eNpYIqZ9QH+BroDOOeWmtkU4A8gDejnnEv3NnUbgZFKooEZ3gNgNDDBzFYRuILd09tWgpkNBeZ6yz3unDvWTGMieSKiRFlqDviCVS9fTtwfD/PBlptB3QBEpIBzG+axb8otlNq7kk8yzmFrm0cZfGFzIsI0GbRIbjhukX2MG106HGX5J4Ens2j/DYjNoj0Jr0jP4r0xgAYsFt+FR5filAGfs+yVbly9axTx40rQ7vqHNYSViBQ8qQdJ/GIoxeaNZL8rw4ulH6XHtf/misol/U4mUqjo11WRbAqLKs5pA6cxJ6IlcWuH89Obg3AZGX7HEhHJNrfuJ/a+0IoS817l/Yz2fBX3CQ8OHMipKrBFcp2KbJETEBoeyf7Wg5lXrjPnbBrNnFd6k5aa4ncsEZFjS05k70cDsbGd2Z14gCfKD6P1gAn0bn8moSH6i5xIXlCRLXKCLCSMZv0n8mv1G2i1axqLX7yCpAOJfscSEclS+spv2PdCC0osHMsE15lfO33Gg/1voVb54n5HEynUVGSL5ICFhNDq3y8x+9R7OTPxZ9a+eBH7du/wO5aIyP87uJvdk/oSOrEr2w7AsKovcuFdb9P97NN0P4lIPlCRLXISWl/zIL+d9SynJC9j5yvns2vzOr8jiYiQ+sdnJL7QnBLL3udt68qfV8zgvr69qVw6yu9oIkWGimyRk9Ty0ptZ2v4tKqZtJXXU+Wz8c47fkUSkiApP2cOu8b0In3Id65OK8WKdkXS5exSXNKujq9ci+UxFtkguaBrXjb+6fEiGg7KTLmXVj+/7HUlEihLnSFnwPmf80o/iqz9nVGhPtvSYzt039KRc8Qi/04kUSccdJ1tEsqdRs3NYX24m68d3p97Mm1m6bRWNu94HunokInlp3xZ2Tu5H+Q1f80dGXb477VFu6taZklHhficTKdJ0JVskF9WofQqV7pjFb1FtaLzoaZa8eTMuPdXvWCJSGDnHgTnjOfhic4qv/46RETcw+4xhDLi2iwpskSCgIlskl5UtW5YzB01jZrlriN30PqteuJiUfTv9jiUihcnuv9nxxmUUm34Hi1OrM/bMidxw93BOq6CuISLBQkW2SB6IiginQ/+RzDjlQWrtm0/Ci2eTsGa+37FEpKDLyCDxx9dJeqkl0ZvnMCL6VqJunsGt3S4iOiLU73QikomKbJE8EhJiXNzrXma3m4ClJRM9vhN/fzfB71giUkC5navZ/uqFlPh6ML+l1+ODVu9zy93/5Ywa5fyOJiJZUJEtksfOO/8Sdv1rJiusLjW/7c/KCQMgPc3vWCJSUGSks2fWcFJGtCFyxx+MKDmQyv2m07tzO8JD9d+4SLDSv06RfHBa/fpUHziTL4pdRv3VY1n7YkfS9mzxO5aIBLmMLX+w/aXzKP3DY/yUcTozzvuE2+58lHoxpfyOJiLHoSJbJJ+UL12SDoPG81HNB6i6dxH7XmpNwpKZfscSkWCUnsrOGU+Q8fq5hOz+i1fL30/9/0yjR4dWhIZoWFCRgkBFtkg+Cg8NodtNg/mx/WQS0qMp80F3/vrgQXUfEZH/Sd3wOzuGn035X59lJq346aLp3N5/MDXKF/c7moicABXZIj7oENcB+n7L1+HtqbVkBBte6kD67g1+xxIRP6Umsf2TBwh563wyErfxRpXHaT7oIy4/+wxNiS5SAKnIFvHJKdUqc+49HzCp2oOU3bOM/S+3IWH+J37HEhEfJK/9hZ3Pt6Liglf53OJY1OVLbrllAJVKRvkdTURySEW2iI+iI0LpefO9/Hj+h2xML0e5ab35e8yNkLTX72gi/2BmncxsuZmtMrP7sng/zsz2mNkC7/GwHzkLlJT9bJl8J+HjLubgwUTerPUc7e6ZwgXNTvM7mYicJBXZIkHgonbnEn17PO8X60G1vz5m53MtSFz+rd+xRP7HzEKBV4GLgUbANWbWKItFf3DONfEej+dryALmwPJvSHiuBZWXjeHj0E6s7zGLm2+8mdLFNCW6SGFwUkW2ma0zs8XeFYvfvLZyZjbTzFZ6X8tmWv5+7wrIcjO7KFN7c287q8zsZfM6n5lZpJlN9tp/NbPaJ5NXJJjVjilL10Gv82GTMexNNUq8dwWbJt0JqQf9jiYC0BJY5Zxb45xLASYBXXzOVDAl7WXjhFso9l5X9iSlM7bBa1w8+B3aNKrtdzIRyUVhubCN9s65HZle3wfMcs497f058T5gsHfFoyfQGKgKfG1mDZxz6cBIoC8wG5gOdAJmAH2AXc65embWExgG9MiFzCJBKSw0hKu7dmPxma35+L276frnGHY8N5NiV46gWIM4v+NJ0VYNWJ/p9QagVRbLtTGzhcAm4G7n3NKsNmZmfQmc94mJiSE+Pv6EwiQmJp7wOnnlRLJEb51LveUjqZyewMSQS0lpch11yhVjzs8/5nuWvBQsOSB4sgRLDlCW/MqRG0X2kboAcd7zcUA8MNhrn+ScSwbWmtkqoKWZrQNKOed+ATCz8cAVBIrsLsCj3rY+AEaYmTnnXB7kFgkap9etSr17JjD+/Qm0W/EkFd7twqa63ana/VmILnv8DYjkvqyGtzjyXDwfqOWcSzSzzsAnQP2sNuacGwWMAmjRooWLi4s7oTDx8fGc6Dp5JTtZ3IEENrw3gBrrp7EiozpfNnmbqy6/nMiw0HzPkh+CJQcET5ZgyQHKkl85TrbIdsBXZuaAN7yTZoxzbjOAc26zmVXylq1G4Er1IRu8tlTv+ZHth9ZZ720rzcz2AOWBzFfOT/qKCATPb1IQPFmCJQcU3Sw1q9Xmm+iXCF/8Lj1Xf8juZ75iVf2bSaxyDon79xfJfVIQckBwZcklG4AamV5XJ3C1+n+cc3szPZ9uZq+ZWYUj/tpZ5Oz67QNCZtxN5bS9TC7Wkyb/eoJrq1X0O5aI5LGTLbLbOuc2eYX0TDP78xjLHu0qyLGujmTnyslJXxGB4PlNCoInS7DkgKKdJQ5IubwTkz/7nNPnP0yLFc+xadev/FXnWtrEXZpvOY4lWI5PsOSA4MqSS+YC9c2sDrCRQPe/azMvYGaVga3OOWdmLQnc97Mz35MGiYy9W9nwbj9qbpnJUleb5a1GclWnTpqxUaSIOKkbH51zm7yv24CPCdwYs9XMqgB4X7d5ix/tKsgG7/mR7YetY2ZhQGkg4WQyixREEWEhXHfFZUTd9i1vl+xLqW1zaTGnP9s+HAzJ+/yOJ0WAcy4N6A98CSwDpjjnlprZrWZ2q7fYVcASr0/2y0DPItm9zzm2/zSe/S80J2bzt0wqeSMl+n1Ht84Xq8AWKUJyXGSbWXEzK3noOdARWAJMA3p7i/UGpnrPpwE9vRFD6hDopzfH61qyz8xae6OKXH/EOoe2dRXwTZE8YYt4GlQpww13PcMPnb5kumtLpcWvs/fZM0j8dTxkZPgdTwo559x051wD59wpzrknvbbXnXOve89HOOcaO+fOdM61ds797G/i/Je2az1/jbiUijPvYHVGFWa1+5Aed71ArUpl/I4mIvnsZK5kxwA/elcs5gCfO+e+AJ4GLjSzlcCF3mu8O8ynAH8AXwD9vJFFAG4D3gJWAasJ3PQIMBoo790keReBkUpEijQz4+I2TYg49y7ebvgWa1LKUmLGHWx98TyS1/zkdzyRosk5Nn/zOskvt6Tijrm8V+52Kg+Mp/P5cZoSXaSIynGfbOfcGuDMLNp3Ah2Oss6TwJNZtP8GxGbRngR0z2lGkcIsOsy4sUd3Vm3tyKjJr9Bl55tEju/M5ph2xFzxJCFVTvc7okiRELp/C3+9eCG19sxlLrHs6zScnq1bqrgWKeI046NIAVcvpjR9//MQa675kbHFelN8y1x441y2jr0el7DW73gihVdGOuunP0+Luf+h3O4lTKk8iHp3f8P5bVqpwBaRPBknW0R80Oa0GrRq8BJf/HYLCV89y1VrPyPt5c/ZUf9qKncejJWt7XdEkULjwKY/2DGxLzX3L+ZHmhDa5WWubvaPP+6KSBGmIlukEAkJMTq3bExys9FM++E3wn54jktWTCZ95WS21+1K5c73YxXq+R1TpOBKT2PttP9SdeFLlHIRfFh7CMWqN+NiFdgicgR1FxEphCLDQunevhWdH5jEZ+1n8FFoJ8qunkrGiLPYPOY6MrZkOdO1iBzD3nW/s/7Zs6mz8Dlmh7VgbY9vufLGu4kO13+lIvJPupItUohFhoXSLa4VKee8w2c/L+Dg9y/T5a8ZhLz+GZvKt6F8h4FEntYRQlQkiBxVWgorP3yEOsveIMUVZ+qpT3FR91uJCs/dKdFFpHBRkS1SBESEhdDtvGaktR3DrN//ZOu3r9NxxzQip/RgZ1QtItreTslWvSCiuN9RRYJKwopfOPjBbdRPWcs3Ee2p2vNFutSt7XcsESkAdPlKpAgJCw3hohaN6HX3S2zo/StvVXqADQdCKTlrMPufbsDG9/6jriQigEs5wJ/jB1L63YsJTd7D56e/xLmDP+I0Fdgikk26ki1SBJkZLerG0OL2wazb3o/x33xKpWUTaP/nREKWj2NLqTMo3qYPJZt319VtKXK2Lv6GjKn9OS1tIzOjL6bedcO5pHpVv2OJSAGjIlukiKtdsQS1e1xDctrVfDNvGdt+HEvb3Z9T+csBHPjqAXbU6kzMOdcTWfcc9d2WQi09aR/L37mbRhsmsd5VYmbLN+lwcXdCQjTmtYicOBXZIgIEbpK8uFUstHqOVVsfZvy3n1F+xWTi1k4jct377AqvxIEGV1Clba/AbJKabEMKkY2/fU749IGclr6dr0p1Jfb657iwYgW/Y4lIAaYiW0T+oV5MKer1vJb0jGuYs2I9q394n+obPqPtktGELB3F1qjaJNe7lCqtryK8WhO/44rkWEriLlZO+A+Nt05jLVX587x3uPD8SzRjo4icNBXZInJUoSFGm9Nq0ua0QRxMGcisBX+ya+5kTtk2k+aLXyV0yQgSwmIoVbI5SVUdUaecC6E6rUjBsPbH9yk1615OzdjNV+Wvpdn1T3NemdJ+xxKRQkL/G4pItkRHhNKpZWNo+ThJqY/w/eLlbPntEypvmkWbhJlEvTud/VaCLRVaU6xhRyo364yVqeF3bJF/OLhrK2sm9KdxwlestFqsvPAtOp7Twe9YIlLIqMgWkRMWFR5K+2aNoFkjUtPv4+0Pp1Mu+W9K/j2LM7fNo/L2r+H7e9kaUZPE6u0offpFVGh4HkTpKqH4yDmWfzOOmB+GUN/tZ2aVf9Oq1+PUL64RdEQk96nIFpGTEh4aQoNKJYiLux24nc27D/DF/F/Zv+xLKm//mearJxG1ZgLpU0PYHFmXA5VbUvq086gU2x4rWdnv+FJE7Nuxnr/H30bjvT+wLKQ+yZ1f5sIWZ/sdS0QKMRXZIpKrqpQpRpXz28P57UnPcCzfsI2/F8aTtvYnKu2aR+y6Dyn217vwJWwNq0pC2TMJrdGcSqe2pkyd5hBRzO9vQQoT51gy/XVqzh3KKS6FWTXv4OzrhhAdFel3MhEp5FRki0ieCQ0xGtWMoVHNHkAPnHOs3bqbVYt+Inn1j5TfOZ9Ttv1CzPYZMB/SCGFLRG32lI0lpGoTytZtSqW6TQkpXtbvb0UKoISNq9gy8VZiD8xlcWhjwru9SofGTf2OJSJFhIpsEck3ZkbdymWpW/lS4FIA9ien8fuqlWxf/jMZG3+n3O4l1NvyDeW2ToPfA+vtCCnPjmL1SCl3KhFVYyld6wwq1mpEWDH18ZZ/chnpLPzkBeovepZazhFfbzBn97yXiHD9lyci+UdnHBHxVfHIMJo2bgiNG/6vbX9SKkvWLCdh7QKSNy0latdyKiau4dR984n8OxVmB5ZLsNLsjKjO/hK1cGXrkJQcycZKkZSvcSpRJcpqwpwiaMvapeyedAtNkhezILwppa9+jbj6jfyOJSJFkIpsEQk6xaPCiW0UC41iD2vfuXc/f67+g8T1i0ndvoqwPWspdeBvKu+YTeWd0wMLrR8OQCLRJIRWZG9EDMnFqpBRqhrh5WpQrEJNilesRZlKNShesowK8UIiIy2N+e8/Seyfr1CMcH5o/ChnXzmA0NAQv6OJSBGlIltECozypYpTvulZ0PSsw9qdc2xN2MV3X3xCjdIhZOxcA3s3EnlgM6VStlLlwArK79wDaw/f3kEXwe6QMuwLLcuBiPKkRJUno1hFQkpUwkrFEFmiAlEly1GsTEVKlqlIyVJlCFHRFnT+/nMeyR/eRovU5cyLbkPla1/j3Jp1/Y4lIkVcgSiyzawT8BIQCrzlnHva50giEkTMjJjy5ahUrS5t4uKyXGZv4j62b1zDvq1/kZKwnvS9W2H/NsIO7iAqeQelkjZS+sBSyuzcS6i5LLeR5kLYbcXZbyU5EFqCpLBSpISVIjW8FC6iOC6yJBZZEossxa7EDCDrLAXV8c7FFpiL/CWgM3AAuME5Nz+v8qSmJDP/3UdounYUiVaM2c2epdWl/8ZC9IuQiPgv6ItsMwsFXgUuBDYAc81smnPuD3+TiUhBUqpESUqdeiaceuYxl0tNTSVhx2YO7NpM0p6dJO3bSer+XaQfSICDu7Gk3YQm7yEidQ/RaXupkLKR6Iz9FHcHibTU/21nfugZwMC8/abyUTbPxRcD9b1HK2Ck9zXX7duygr+HDaRV+lrmlmxP7V4jaB1TPS8+SkQkR4K+yAZaAqucc2sAzGwS0AVQkS0iuS48PJyKVWpClZonvK5LS+bAvj0cSNzFtvm/50E6X2XnXNwFGO+cc8BsMytjZlWcc5tzM8ji76dy8bLB7LZSzG8zgrMu6pWbmxcRyRUFociuBqzP9HoDR1wZMbO+QF/vZaKZLc/B51QAduQoYe4LlizBkgOUJSvBkgOCJ0uw5ACoANfkJEutXE+SO457Lj7KMtWAfxTZuXDergC7d8D1BB6+CrKfu6DIEiw5IHiyBEsOUJas5DTHUc/ZBaHIzurW/8M6TDrnRgGjTupDzH5zzrU4mW3klmDJEiw5QFmCOQcET5ZgyQHBlSWXHPdcnM1lAo0ned4Opv2rLMGbA4InS7DkAGXJrxwF4e6QDUCNTK+rA5t8yiIiUlRl51ys87WIiKcgFNlzgfpmVsfMIoCewDSfM4mIFDXZORdPA663gNbAntzujy0iUlAEfXcR51yamfUHviQwbNQY59zSPPiok+puksuCJUuw5ABlyUqw5IDgyRIsOSC4spy0o52LzexW7/3XgekEhu9bRWAIvxvzMFIw7V9l+adgyQHBkyVYcoCyZCXXc1jgJnAREREREcktBaG7iIiIiIhIgaIiW0REREQkl6nIJjBVsJktN7NVZnZfPn/2OjNbbGYLzOw3r62cmc00s5Xe17J59NljzGybmS3J1HbUzzaz+719tNzMLsqHLI+a2UZv3ywws855ncXMapjZt2a2zMyWmtkArz1f98sxcvixT6LMbI6ZLfSyPOa15/vPyjGy5Pt+8bYdama/m9ln3mtf/v0UNTpn65ydabtBcc4+TpZ83S86Zx8zT/6es51zRfpB4Aae1UBdIAJYCDTKx89fB1Q4ou0Z4D7v+X3AsDz67POAZsCS43020MjbN5FAHW+fheZxlkeBu7NYNs+yAFWAZt7zksAK7/Pydb8cI4cf+8SAEt7zcOBXoLUfPyvHyJLv+8Xb/l3Au8Bnfv77KUoPdM7WOfvwbQfFOfs4WfJ1vxzjPKlzdj6fs3UlO9NUwc65FODQVMF+6gKM856PA67Iiw9xzn0PJGTzs7sAk5xzyc65tQRGD2iZx1mOJs+yOOc2O+fme8/3AcsIzFiXr/vlGDmOJi/3iXPOJXovw72Hw4eflWNkOZo8y2Jm1YFLgLeO+Lx8//dTxOicnb3P1jk7/89PQXHe1jk7a36cs1VkH30a4PzigK/MbJ4FphkGiHHe2LLe10r5mOdon+3XfupvZou8P00e+jNOvmQxs9pAUwK/efu2X47IAT7sE+9PbAuAbcBM55xv++QoWSD/98uLwL1ARqa2YPv3Uxj5vS91zj62In/OziIL5PN+0Tk7Sy+Sz+dsFdknMA1wHmnrnGsGXAz0M7Pz8vGzT4Qf+2kkcArQBNgMPJ9fWcysBPAhMNA5t/dYi+Zllixy+LJPnHPpzrkmBGbwa2lmsceK7UOWfN0vZnYpsM05Ny+7q+RFjiLK732pc/bRFflz9lGy5Pt+0Tn7cH6ds1Vk+zwNsHNuk/d1G/AxgT9HbDWzKgDe1235lecYn53v+8k5t9X7x5kBvMn//6kmT7OYWTiBE+RE59xHXnO+75escvi1Tw5xzu0G4oFO+PyzkjmLD/ulLXC5ma0j0F3hfDN7hyD691OI6Zx9uKD5mSvq5+yjZfHzvK1z9v/4cs5Wke3jtO1mVtzMSh56DnQElnif39tbrDcwNT/yeI722dOAnmYWaWZ1gPrAnLwMcugH39OVwL7J0yxmZsBoYJlzbnimt/J1vxwth0/7pKKZlfGeRwMXAH/iw8/K0bLk935xzt3vnKvunKtN4JzxjXPuXwTRv59CTOfswwXNz1xRPmcfK0t+7xeds//Jt3O2y6U7Ngvyg8A0wCsI3D36YD5+bl0Cd68uBJYe+mygPDALWOl9LZdHn/8egT/TpBL4ra3PsT4beNDbR8uBi/MhywRgMbDI+4GvktdZgHMI/EloEbDAe3TO7/1yjBx+7JMzgN+9z1wCPHy8n1MfsuT7fsm0/Tj+/051X/79FLUHOmfrnP3/2w2Kc/ZxsuTrfjnGeVLnbJe/52xNqy4iIiIiksvUXUREREREJJepyBYRERERyWUqskVEREREcpmKbBERERGRXKYiW0REREQkl6nIFhERERHJZSqyRURERERy2f8BnI2ShZrToEYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_func(steps-100)\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": true,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}