{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The heat equation in 2D\n",
    "\n",
    "$$\n",
    "\\renewcommand{\\DdQq}[2]{{\\mathrm D}_{#1}{\\mathrm Q}_{#2}}\n",
    "\\renewcommand{\\drondt}{\\partial_t}\n",
    "\\renewcommand{\\drondx}{\\partial_x}\n",
    "\\renewcommand{\\drondtt}{\\partial_{tt}}\n",
    "\\renewcommand{\\drondxx}{\\partial_{xx}}\n",
    "\\renewcommand{\\drondyy}{\\partial_{yy}}\n",
    "\\renewcommand{\\dx}{\\Delta x}\n",
    "\\renewcommand{\\dt}{\\Delta t}\n",
    "\\renewcommand{\\grandO}{{\\mathcal O}}\n",
    "\\renewcommand{\\density}[2]{\\,f_{#1}^{#2}}\n",
    "\\renewcommand{\\fk}[1]{\\density{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\fks}[1]{\\density{#1}{\\star}}\n",
    "\\renewcommand{\\moment}[2]{\\,m_{#1}^{#2}}\n",
    "\\renewcommand{\\mk}[1]{\\moment{#1}{\\vphantom{\\star}}}\n",
    "\\renewcommand{\\mke}[1]{\\moment{#1}{e}}\n",
    "\\renewcommand{\\mks}[1]{\\moment{#1}{\\star}}\n",
    "$$\n",
    "\n",
    "In this tutorial, we test a very classical lattice Boltzmann scheme $\\DdQq{2}{5}$ on the heat equation.\n",
    "\n",
    "The problem reads\n",
    "$$\n",
    "\\begin{gathered} \\drondt u = \\mu (\\drondxx+\\drondyy) u, \\quad t>0, \\quad (x, y)\\in(0,1)^2,\\\\ u(0) = u(1) = 0, \\end{gathered}\n",
    "$$\n",
    "\n",
    "where $\\mu$ is a constant scalar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The scheme $\\DdQq{2}{5}$\n",
    "\n",
    "The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximatation of the solution on discret points of $(0,1)^2$ at discret instants.\n",
    "\n",
    "To simulate this system of equations, we use the $\\DdQq{2}{5}$ scheme given by\n",
    "\n",
    "* five velocities $v_0=(0,0)$, $v_1=(1,0)$, $v_2=(0,1)$, $v_3=(-1,0)$, and $v_4=(0,-1)$ with associated distribution functions $\\fk{i}$, $0\\leq i\\leq 4$,\n",
    "* a space step $\\dx$ and a time step $\\dt$, the ration $\\lambda=\\dx/\\dt$ is called the scheme velocity,\n",
    "* five moments\n",
    "  $$ \\mk{0}=\\sum_{i=0}^{4} \\fk{i}, \\quad \\mk{1}= \\sum_{i=0}^{4} v_{ix} \\fk{i}, \\quad \\mk{2}= \\sum_{i=0}^{4} v_{iy} \\fk{i}, \\quad \\mk{3}= \\frac{1}{2} \\sum_{i=0}^{5} (v_{ix}^2+v_{iy}^2) \\fk{i}, \\quad \\mk{4}= \\frac{1}{2} \\sum_{i=0}^{5} (v_{ix}^2-v_{iy}^2) \\fk{i},$$\n",
    "  \n",
    "  and their equilibrium values $\\mke{k}$, $0\\leq k\\leq 4$.\n",
    "* two relaxation parameters $s_1$ and $s_2$ lying in $[0,2]$ ($s_1$ for the odd moments and $s_2$ for the odd ones).\n",
    "\n",
    "In order to use the formalism of the package pylbm, we introduce the five polynomials that define the moments: $P_0 = 1$, $P_1=X$, $P_2=Y$, $P_3=(X^2+Y^2)/2$, and $P_4=(X^2-Y^2)/2$, such that\n",
    "$$ \n",
    "\\mk{k} = \\sum_{i=0}^4 P_k(v_{ix}, v_{iy}) \\fk{i}.\n",
    "$$\n",
    "\n",
    "The transformation $(\\fk{0}, \\fk{1}, \\fk{2}, \\fk{3}, \\fk{4})\\mapsto(\\mk{0},\\mk{1}, \\mk{2}, \\mk{3}, \\mk{4})$ is invertible if, and only if, the polynomials $(P_0,P_1,P_2,P_3,P_4)$ is a free set over the stencil of velocities.\n",
    "\n",
    "The lattice Boltzmann method consists to compute the distribution functions $\\fk{i}$, $0\\leq i\\leq 4$ in each point of the lattice $x$ and at each time $t^n=n\\dt$.\n",
    "A step of the scheme can be read as a splitting between the relaxation phase and the transport phase:\n",
    "\n",
    "* relaxation: \n",
    "$$\n",
    "    \\begin{aligned}\\mks{1}(t,x,y)&=(1-s_1)\\mk{1}(t,x,y)+s_1\\mke{1}(t,x,y),\\\\ \\mks{2}(t,x,y)&=(1-s_1)\\mk{2}(t,x,y)+s_1\\mke{2}(t,x,y),\\\\ \\mks{3}(t,x,y)&=(1-s_2)\\mk{3}(t,x,y)+s_2\\mke{3}(t,x,y),\\\\ \\mks{4}(t,x,y)&=(1-s_2)\\mk{4}(t,x,y)+s_2\\mke{4}(t,x,y).\\end{aligned}\n",
    "$$\n",
    "\n",
    "* m2f:\n",
    "$$\n",
    "    \\begin{aligned}\\fks{0}(t,x,y)&\\;=\\mk{0}(t,x,y)-2\\mks{3}(t,x,y), \\\\ \\fks{1}(t,x,y)&\\;=\\tfrac{1}{2}(\\phantom{-}\\mks{1}(t,x,y)+\\mks{3}(t,x,y)+\\mks{4}(t,x,y)), \\\\ \\fks{2}(t,x,y)&\\;=\\tfrac{1}{2}(\\phantom{-}\\mks{2}(t,x,y)+\\mks{3}(t,x,y)-\\mks{4}(t,x,y)), \\\\ \\fks{3}(t,x,y)&\\;=\\tfrac{1}{2}(-\\mks{1}(t,x,y)+\\mks{3}(t,x,y)+\\mks{4}(t,x,y)), \\\\ \\fks{4}(t,x,y)&\\;=\\tfrac{1}{2}(-\\mks{2}(t,x,y)+\\mks{3}(t,x,y)-\\mks{4}(t,x,y)).\\end{aligned}\n",
    "$$\n",
    "\n",
    "* transport: \n",
    "$$\n",
    "    \\begin{aligned} \\fk{0}(t+\\dt, x,y)&\\;=\\fks{0}(t,x,y), \\\\ \\fk{1}(t+\\dt, x,y)&\\;=\\fks{1}(t,x-\\dx,y), \\\\ \\fk{2}(t+\\dt, x,y)&\\;=\\fks{2}(t,x,y-\\dx), \\\\ \\fk{3}(t+\\dt, x,y)&\\;=\\fks{3}(t,x+\\dx,y), \\\\ \\fk{4}(t+\\dt, x,y)&\\;=\\fks{4}(t,x,y+\\dx). \\end{aligned}\n",
    "$$\n",
    "\n",
    "* f2m:\n",
    "$$\n",
    "    \\begin{aligned}\\mk{0}(t+\\dt,x,y)&\\;=\\fk{0}(t+\\dt,x,y)+\\fk{1}(t+\\dt,x,y)+\\fk{2}(t+\\dt,x,y)\\\\&\\;\\phantom{=}+\\fk{3}(t+\\dt,x,y)+\\fk{4}(t+\\dt,x,y), \\\\ \\mk{1}(t+\\dt,x,y)&\\;=\\fk{1}(t+\\dt,x,y)-\\fk{3}(t+\\dt,x,y), \\\\ \\mk{2}(t+\\dt,x,y)&\\;=\\fk{2}(t+\\dt,x,y)-\\fk{4}(t+\\dt,x,y), \\\\ \\mk{3}(t+\\dt,x,y)&\\;=\\tfrac{1}{2}(\\fk{1}(t+\\dt,x,y)+\\fk{2}(t+\\dt,x,y)+\\fk{3}(t+\\dt,x,y)+\\fk{4}(t+\\dt,x,y)), \\\\ \\mk{4}(t+\\dt,x,y)&\\;=\\tfrac{1}{2}(\\fk{1}(t+\\dt,x,y)-\\fk{2}(t+\\dt,x,y)+\\fk{3}(t+\\dt,x,y)-\\fk{4}(t+\\dt,x,y)).\\end{aligned}\n",
    "$$\n",
    "\n",
    "The moment of order $0$, $\\mk{0}$, being conserved during the relaxation phase, \n",
    "a diffusive scaling $\\dt=\\dx^2$, yields to the following equivalent equation\n",
    "$$\n",
    "\\drondt\\mk{0} = \\bigl(\\tfrac{1}{s_1}-\\tfrac{1}{2}\\bigr) \\bigl(\\drondxx(\\mke{3}+\\mke{4})+\\drondyy(\\mke{3}-\\mke{4})\\bigr) + \\grandO(\\dx^2),\n",
    "$$\n",
    "\n",
    "if $\\mke{1}=0$.\n",
    "In order to be consistent with the heat equation, the following choice is done:\n",
    "$$\n",
    "\\mke{3}=\\tfrac{1}{2}u, \\qquad \\mke{4}=0, \\qquad s_1 = \\frac{2}{1+4\\mu}, \\qquad s_2=1.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using pylbm\n",
    "\n",
    "pylbm uses Python dictionary to describe the simulation. In the following, we will build this dictionary step by step.\n",
    "\n",
    "### The geometry\n",
    "\n",
    "In pylbm, the geometry is defined by a box and a label for the boundaries. We define here a square $(0, 1)^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+----------------------+\n",
      "| Geometry information |\n",
      "+----------------------+\n",
      "    - spatial dimension: 2\n",
      "    - bounds of the box: [0. 1.] x [0. 1.]\n",
      "    - labels: [0, 0, 0, 0]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD8CAYAAACcoKqNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAETxJREFUeJzt3X9s3PV9x/HnKz8tQ1oQdqE4AadaUHBJBezqda0goLLOeCKRotL8oOqYUKJ2AQRJJzEVQRvUP0ZXEsVka50N0RY1afpDq1WCmNalYqoaYqfQgBNl9dJ0mJTGbmnVxQpJ8Ht/fA96uZxzX+Lv3X2TvB6S5ft+vx/u++Kce93nvnf3PUUEZnZ+m9LoAGbWeC4CM3MRmJmLwMxwEZgZLgIzw0VgZrgIzAwXgZkB0xq145aWlmhvb2/U7s3OC7t37x6NiNZq4xpWBO3t7QwMDDRq92bnBUm/TDPOTw3MzEVgExsfH2f9+vXMnz+fpqYm5syZw9q1azly5Eijo1nGXAQ2ofvvv581a9bQ0dFBT08Pt99+Oxs3buS2225jfHy80fEsQw07RmD5Njg4SE9PD0uWLOG73/3u2+vnzp3Lvffey9atW1mxYkUDE1qWPCOwirZs2UJEcN999520fuXKlTQ3N/PUU081KJnVQtUikPSEpMOSXp5guyRtlDQkaY+k67OPafXW39/PlClT6OzsPGl9U1MT1157Lf39/Q1KZrWQZkbwJNB1mu23AvOKP6uAf558LGu0Q4cO0dLSwsyZM0/Z1tbWxujoKMeOHWtAMquFqkUQEc8Bvz3NkMXA1yOxE7hI0nuzCmiNMTY2VrEEIJkVvDXGzg1ZHCxsA14pWR4urvtV+UBJq0hmDVxxxRWprvyyy/6RX//aL1fV3/8BR5C+UGHbbgAuvvjL+HhzY1x66QW89tpnM7u+LA4WqsK6imdEjYjeiChERKG1teq7HgFcAg0zCxgDTlTY9gegGZdA42R9v8iiCIaBOSXLs4FDGVyvNdTlJH3+atn648Brxe12rsiiCPqATxVfPfgQ8PuIOOVpgZ1trin+3lm2/qckZbCgvnGspqrO7SRtAW4CWiQNAw8D0wEi4ivAdqAbGCKZS/5NrcJaPV0KdAK7gK0kLwqNAs8DV+IiOLdULYKIWF5lewCrM0tkOdIFXERycPDnJMcFOoGb8XvRzi0+2mOnMQX4cPHHzmWudTNzEZiZi8DMcBGYGS4CM8NFYGa4CMwMF4GZ4SIwM1wEZoaLwMxwEZgZLgIzw0VgZrgIzAwXgZnhIjAzXARmhovAzHARmBkuAjPDRWBmuAjMDBeBmeEiMDNcBGaGi8DMcBGYGS4CMyNlEUjqkrRf0pCkBypsv0LSDkkvSNojqTv7qGZWK1WLQNJUYBNwK9ABLJfUUTbsQWBbRFwHLAP+KeugZlY7aWYEncBQRByIiGPAVmBx2ZgA3lW8/G7gUHYRzazWpqUY0wa8UrI8DPxZ2ZjPA/8u6R7gAuCWTNKZWV2kmRGowrooW14OPBkRs4Fu4BuSTrluSaskDUgaGBkZeedpzawm0hTBMDCnZHk2p0797wK2AUTET4AmoKX8iiKiNyIKEVFobW09s8Rmlrk0RdAPzJM0V9IMkoOBfWVj/hf4KICkq0mKwA/5ZmeJqkUQESeAu4FngX0krw4MSlonaVFx2FpgpaSfAVuAOyOi/OmDmeVUmoOFRMR2YHvZuodKLu8FPpJtNDOrF7+z0MxcBGbmIjAzXARmhovAzHARmBkuAjPDRWBmuAjMDBeBmeEiMDNcBGaGi8DMcBGYGS4CM8NFYGa4CMwMF4GZ4SIwM1wEZoaLwMxwEZgZLgIzw0VgZrgIzAwXgZnhIjAzXARmhovAzHARmBkuAjMjZRFI6pK0X9KQpAcmGPMJSXslDUr6ZrYxzayWplUbIGkqsAn4C2AY6JfUFxF7S8bMA/4e+EhEvC7pPbUKbGbZSzMj6ASGIuJARBwDtgKLy8asBDZFxOsAEXE425hmVktpiqANeKVkebi4rtRVwFWSfixpp6SuSlckaZWkAUkDIyMjZ5bYzDKXpghUYV2ULU8D5gE3AcuBf5F00Sn/UURvRBQiotDa2vpOs5pZjaQpgmFgTsnybOBQhTHfj4jjEfELYD9JMZjZWSBNEfQD8yTNlTQDWAb0lY35N+BmAEktJE8VDmQZ1Mxqp2oRRMQJ4G7gWWAfsC0iBiWtk7SoOOxZ4DeS9gI7gL+LiN/UKrSZZavqy4cAEbEd2F627qGSywGsKf6Y2VnG7yw0MxeBmbkIzAwXgZnhIjAzXARmRg6LYHx8nPXr1zN//nyampqAx0jepnCswcnMGmkc+AnQAzwCPMbatWs5cuRIJteeuyK4//77WbNmDR0dHfT09AAdwPPAN0luDLPz0bPFn1agG+hg48aN3HbbbYyPT/5+kasiGBwcpKenhyVLlvC9732PlStXAl3AXwIHgZcbms+sMQ6TPBheTfIO/z8FunjsscfYsWMHW7dunfQeclUEW7ZsISK47777yrZcD0wH9jQglVmjvVT8/aGT1q5cuZLm5maeeuqpSe8hV0XQ39/PlClT6OzsLNsyHbiMUz/0aHY+OERyNoCTTwPS1NTEtddeS39//6T3kKsiOHToEC0tLcycObPC1lnAGHCizqnMGu0PQDOVPhrU1tbG6Ogox45N7mB6ropgbGxsghKAP94Ix+sVxywnjgNTK25JXllL7juTkasiaG5u5o033phg61szgen1imOWE9OBNytuOXr0KJDcdyYjV0Vw+eWXMzo6OkEZTDw9Mju3Tfy0+NVXX6WlpYUZM2ZMag+5KoIPfvCDjI+Ps2vXrrItx4HXgMsbkMqs0S4nOU3oqyetPXr0KC+++CKFQmHSe8hVESxduhRJbNiwoWzLT0nKYEEDUpk12jXF3ztPWrt582bGxsa44447Jr2HXM2zFyxYwOrVq3n88cdZsmQJ3d3dJO+meh64EheBnZ8uJfl6kV0kXysyDxhlzZp+Fi5cyIoVKya9h1wVAcCGDRtob2+nt7eXp59+GphJciPcTM4mMGZ11AVcBOwGfg40c88997Bu3TqmTJn8/ULJ6Qbrr1AoxMDAQNVx0hfqkMbs7BPxcNUxknZHRNWDCH6INTMXgZm5CMwMF4GZ4SIwM1wEZoaLwMxwEZgZLgIzw0VgZqQsAkldkvZLGpL0wGnGfVxSSJr85yLNrG6qFoGkqcAm4FaSLxlYLqmjwrhZwL0kHxU0s7NImhlBJzAUEQci4hjJ5yAXVxj3CPAocDTDfGZWB2mKoA14pWR5mLLzKku6DpgTET/IMJuZ1UmaIlCFdW9/dlnSFGA9sLbqFUmrJA1IGhgZGUmf0sxqKk0RDANzSpZnc/I3jcwiOZfSjyQdJPk6lr5KBwwjojciChFRaG1tPfPUZpapNEXQD8yTNFfSDJIvX+t7a2NE/D4iWiKiPSLaSU6stigiqp91xMxyoWoRRMQJ4G6SkwfuA7ZFxKCkdZIW1TqgmdVeqnMWRsR2YHvZuocmGHvT5GOZWT35nYVm5iIwMxeBmeEiMDNcBGaGi8DMcBGYGS4CM8NFYGa4CMwMF4GZ4SIwM1wEZoaLwMxwEZgZLgIzw0VgZrgIzAwXgZnhIjAzXARmhovAzHARmBkuAjPDRWBmuAjMDBeBmeEiMDNcBGaGi8DMcBGYGSmLQFKXpP2ShiQ9UGH7Gkl7Je2R9ENJV2Yf1cxqpWoRSJoKbAJuBTqA5ZI6yoa9ABQi4gPAd4BHsw5qZrWTZkbQCQxFxIGIOAZsBRaXDoiIHRExVlzcCczONqaZ1VKaImgDXilZHi6um8hdwDOVNkhaJWlA0sDIyEj6lGZWU2mKQBXWRcWB0ieBAvClStsjojciChFRaG1tTZ/SzGpqWooxw8CckuXZwKHyQZJuAT4HLIyIN7KJZ2b1kGZG0A/MkzRX0gxgGdBXOkDSdcBXgUURcTj7mGZWS1WLICJOAHcDzwL7gG0RMShpnaRFxWFfAi4Evi3pRUl9E1ydmeVQmqcGRMR2YHvZuodKLt+ScS4zqyO/s9DMXARm5iIwM1wEZoaLwMxwEZgZLgIzw0VgZrgIzAwXgZnhIjAzXARmhovAzHARmBkuAjPDRWBmuAjMDBeBmeEiMDNcBGaGi8DMcBGYGS4CM8NFYGa4CMwMF4GZkfIrz+x8NQ48DwwAvwMuAN4P3AzMaGAuy5qLwE7jWZIimA98GBgpLv8K+BSeUJ47XAQ2gcMkd/qrgaUl6y8GngFeBj7QgFxWC650m8BLxd8fKlt/PTAd2FPfOFZTqYpAUpek/ZKGJD1QYftMSd8qbn9eUnvWQa3eDgEC2srWTwcuK263c0XVIpA0FdgE3Ap0AMsldZQNuwt4PSL+BFgP/EPWQa3e/gA0U/nZ4yxgDDhR10RWO2lmBJ3AUEQciIhjwFZgcdmYxcDXipe/A3xUkrKLafV3HJg6wbZpJWPsXJCmCNqAV0qWhzl1vvj2mIg4AfweuCSLgNYo04E3J9h2omSMnQvSFEGlR/Y4gzFIWiVpQNLAyMhImnxceukFqcZZ1k43/T/d0warh6zvF2n+ksPAnJLl2Zx6pOitMcOSpgHvBn5bfkUR0Qv0AhQKhVOKopLXXvtsmmGWsQcfPM4Xv/hFnnvuY9xwww1vrz969CiXXPIoN954I88883ADE1qW0swI+oF5kuZKmgEsA/rKxvQBf128/HHgPyMi1R3d8mnp0qVIYsOGDSet37x5M2NjY9xxxx0NSma1UHVGEBEnJN1N8jazqcATETEoaR0wEBF9wL8C35A0RDITWFbL0FZ7CxYsYPXq1Tz++OMsWbKE7u5u9u3bx8aNG1m4cCErVqxodETLkBr1wF0oFGJgYKAh+7Z03nzzTTZs2EBvby8HDx6kpaWFpUuXsm7dOi688MJGx7MUJO2OiELVcS4Cs3NX2iLwW4zNzEVgZi4CM8NFYGY08GChpBHglymGtgCjNY5zpvKcDZxvMvKcDdLnuzIiWqsNalgRpCVpIM1Rz0bIczZwvsnIczbIPp+fGpiZi8DMzo4i6G10gNPIczZwvsnIczbIOF/ujxGYWe2dDTMCM6ux3BRBnk+QmiLbGkl7Je2R9ENJV9YrW5p8JeM+Likk1e1oeJpskj5RvP0GJX2zXtnS5JN0haQdkl4o/n2765jtCUmHJb08wXZJ2ljMvkfS9We8s4ho+A/Jx5v/B3gfyVfo/AzoKBvzt8BXipeXAd/KUbabgebi5c/UK1vafMVxs4DngJ1AIS/ZgHnAC8DFxeX35Om2I3ku/pni5Q7gYB3z3Uhy/viXJ9jeTfIlEyI57/zzZ7qvvMwI8nyC1KrZImJHRIwVF3eSnMWpXtLcdgCPAI8CR3OWbSWwKSJeB4iIwznLF8C7ipffTR3P4x4Rz1HhTF8lFgNfj8RO4CJJ7z2TfeWlCPJ8gtQ02UrdRdLS9VI1n6TrgDkR8YM65oJ0t91VwFWSfixpp6SuuqVLl+/zwCclDQPbgXvqEy2Vd/pvc0J5OftkZidIrYHU+5X0SaAALKxporLdVlj3dj5JU0i+a+LOegUqkea2m0by9OAmkpnUf0m6JiJ+V+NskC7fcuDJiPiypD8nORPXNRExXvt4VWV2n8jLjOCdnCCV050gtUHZkHQL8DlgUUS8UYdcb6mWbxZwDfAjSQdJnkv21emAYdq/6/cj4nhE/ALYT1IM9ZAm313ANoCI+AnQRPI+/zxI9W8zlXod+KhyUGQacACYyx8P2ry/bMxqTj5YuC1H2a4jOeg0L4+3Xdn4H1G/g4Vpbrsu4GvFyy0kU91LcpTvGeDO4uWri3c01fHv287EBwv/ipMPFu464/3U638oxf9wN/DfxTvU54rr1pE8wkLSxN8GhoBdwPtylO0/gF8DLxZ/+vJ025WNrVsRpLztBDwG7CX55tVlebrtSF4p+HGxJF4EPlbHbFtIvoP+OMmj/13Ap4FPl9x2m4rZX5rM39XvLDSz3BwjMLMGchGYmYvAzFwEZoaLwMxwEZgZLgIzw0VgZsD/A4wzxVoWY8tjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylbm\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "xmin, xmax, ymin, ymax = 0., 1., 0., 1.\n",
    "dico_geom = {\n",
    "    'box': {'x': [xmin, xmax], \n",
    "            'y': [ymin, ymax], \n",
    "            'label': 0\n",
    "           },\n",
    "}\n",
    "geom = pylbm.Geometry(dico_geom)\n",
    "print(geom)\n",
    "geom.visualize(viewlabel=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The stencil\n",
    "\n",
    "pylbm provides a class stencil that is used to define the discret velocities of the scheme. In this example, the stencil is composed by the velocities $v_0=(0,0)$, $v_1=(1,0)$, $v_2=(-1,0)$, $v_3=(0,1)$, and $v_4=(0,-1)$ numbered by $[0,1,2,3,4]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+---------------------+\n",
      "| Stencil information |\n",
      "+---------------------+\n",
      "    - spatial dimension: 2\n",
      "    - minimal velocity in each direction: [-1 -1]\n",
      "    - maximal velocity in each direction: [1 1]\n",
      "    - information for each elementary stencil:\n",
      "        stencil 0\n",
      "            - number of velocities: 5\n",
      "            - velocities\n",
      "                (0: 0, 0)\n",
      "                (1: 1, 0)\n",
      "                (2: 0, 1)\n",
      "                (3: -1, 0)\n",
      "                (4: 0, -1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAE/CAYAAAAjcYRfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEUNJREFUeJzt3H2QXXV5wPHvkyUJhCUsLyGETTAIBgREVlLkxRGqUiLD1EBlKjhOqC+pU3GITUdh7AjSmc7UCuI0WBoggwUKtUWQ8YVIlJRGXs0bEiI0IoElkBhgyS7kbZNf/9iLs6tMdpN72bP77Pczs+Oeveee88yZ65dzz70nUUpBkrIZVfUAkvR2MG6SUjJuklIybpJSMm6SUjJuklIybhp2IuLwiOiKiKba8uKI+GzVc2loMW7aLRHxgYh4MCJei4hXIuIXEfEntccujoglb/cMpZTnSinNpZQdA1k/Ir4UES/VZl4QEWPf7hlVPeOmAYuI8cAPgX8BDgRaga8DW6uca1ci4mzgMuDDwFTgnfTMrOSMm3bHNIBSyu2llB2llM2llJ+WUh6PiHcD1wOn1t4ydgBExNiI+GZEPBcR6yPi+ojYp/bYmRHRHhFzI2JDRLwYEX/15s4iYp+IuDoi1tbOupbU/jY1IkpE7DWAmWcBN5VSVpVSXgX+Abi4wcdFQ5Bx0+54GtgREd+NiI9GxAFvPlBKWQ18Hnio9paxpfbQP9ETxROBo+g52/tar20eCuxf+/tngOt6bfebwEnAafScKX4Z2LmbMx8HrOy1vBKYGBEH7eZ2NMwYNw1YKWUT8AGgADcAv4uIeyJi4lutHxEBfA74UinllVJKJ/CPwCd6rbYduKqUsr2U8mOgCzg6IkYBnwYuLaW8UDtTfLCUsrtvgZuB13otv/n7fru5HQ0zAzmtl36vdoZ2MUBEHAPcClwLXPgWq08AxgFLezoHQABNvdZ5uZTS3Wv5DXqCdDCwN/CbOkfuAsb3Wn7z9846t6shzjM37bFSyq+Bm4Hj3/zTH6yyEdgMHFdKaan97F9KaR7A5jcCW4Aj6xxzFfDeXsvvBdaXUl6uc7sa4oybBiwijqld/J9cW55Czxnbw7VV1gOTI2IMQCllJz1vX78VEYfUntNa+wRzl2rPXQBcExGHRURTRJy6B1/j+HfgMxFxbO1a3t/TE2QlZ9y0OzqB9wOPRMTr9ETtCWBu7fGf03Om9FJEbKz97SvAGuDhiNgELAKOHuD+/g74FfAY8Ao9H07s1mu2lHIv8A3gfmBt7eeK3dmGhqfwH6uUlJFnbpJSqjtuETElIu6PiNURsSoiLm3EYJJUj7rflkbEJGBSKWVZROwHLAVmllKebMSAkrQn6j5zK6W8WEpZVvu9E1hNz7fNJakyDb3mFhFTgTbgkUZuV5J2V8PuUIiIZuBOYE7tNp0/fHw2MBtg3LhxJ02bNq1Rux72duzYQVNTU/8rjgDR3cXOnTuJMeP7X3mE8PXR14oVKzaWUib0t15DvgoSEaPp+adwFpZSrulv/ba2trJ8+fK695tFR0cHLS0t/a84EqxfTFdXF81Hnlv1JEOGr4++ImJpKWV6f+s14tPSAG4CVg8kbJI0GBpxze104FPAhyJiRe3nnAZsV5L2WN3X3EopS+j5lx4kacjwDgVJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3DQmf+9w9HHvsdTS/8385fPoKzjnnNlat2lD1WBrGGhK3iFgQERsi4olGbE8jz403Lmf8+LFceN4hjG9u4ic/WcPZZ9/Kli3dVY+mYapRZ243AzMatC2NQL/4xad5+OHPcsPVR/OjW6cB8MILnTz55O8qnkzDVUPiVkp5AHilEdvSyHTaaVN+//u27QWAUaOCSZOaqxpJw5zX3DSkdL2+g89/+VkA5s49lUmT9qt2IA1bew3WjiJiNjAboLW1lY6OjsHa9ZDX2dlZ9QhDwssvb+Yvz1/G0sffYNas47n88um+TvD1sacGLW6llPnAfIC2trbS0tIyWLseFkb68Vi7toNzzrmTp59+g7/960O5+vq/qHqkIWWkvz72xKDFTdqV005bwLp1nRzeOpYtW3cyZ869AFx00Xs4+eTWiqfTcNSQuEXE7cCZwMER0Q5cUUq5qRHb1siwbl3PW6/nXtjKd27eAPR8x+3EEw81btojDYlbKeXCRmxHI1cpV/T8sn4xXV1dNB95brUDadjz01JJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRk3SSkZN0kpGTdJKRm3Cs2adTetrdcwceI8Dj74G8yYcSvLl79Y9ViV2LKlmy9+8cccctyDTDhuGaefvoBHHmmveqzKXHvtw5xwwr/S1HQVBxzwba68cnHVIw07DYlbRMyIiKciYk1EXNaIbY4Ea9d2cMYZ7+CTnzyWgw4ax8KFv2HmzP+seqxKzJlzL/PmPcbECWM496wWHnroec466xY2bnyj6tEqsXTpixx44D5MmTK+6lGGrb3q3UBENAHXAWcB7cBjEXFPKeXJered3eLFFwPQ0dHBM89s5qST5tPevont23cwenRTtcMNog0bXmfBguWMGhX87L9OYNw+29i7ZTu33vo48+Y9ypVXnln1iIPullvOA2DmzDtYu/a1iqcZnhpx5nYysKaU8kwpZRtwB/CxBmx3RJg371Hmzv05F154JwBz5546osIGsGrVBrZv38nhh+/PIRPGADB9+iQAVqx4qcrRNIzVfeYGtALP91puB96/qydEdxesX9yAXQ9///0fK/ifh3r+yzz5sLGcftxrI+7YrH96AwDNe2+HV1fQtHkz++6YCsBLz68bccejj60be/6369mRfRz2QCPiFm/xt/JHK0XMBmYDHDH5YLq6uhqw6+Hvh7ccxSuvdvHgL7v55Bd+w8c/u4oVi47nHZPHVj3aoBm/7w4AOru2s3nzdrZs2cLLr/Rcazv4wFEj+rXS3d0NwLZt20b0cdgTjYhbOzCl1/JkYN0frlRKmQ/MB2hrayvNR57bgF0PX5s3b2fMmCaamkbR3dHBee9ppvkr32TTpq1s2Pk+jjvyiKpHHDQnNXcxevS3eH7ddjaNPoV9x2/j8d9uA9o56bT30nzkn1Y9YmX22vcO4DXGHDiN5iPPrHqcYaURcXsMeFdEHAG8AHwCuKgB203tkUde4KKL7uSDH3wH48aN4tFHX2LTpq1MmDCO971vUtXjDaqJE5u5+OITueGGZXz44ys55qixfP/Hr9LcPIZLLjm56vEqceONy1iy5DmWLev5atDdd/+aZ5/tYObMY5g585iKpxse6o5bKaU7Ii4BFgJNwIJSyqq6J0vusMP2Y9q0g7jvvmfo7NzKhAn7csEFx/K1r53B/vvvXfV4g+7b357B6NGj+N4dK1nz282ccspkrr76z5gwYd+qR6vEkiXP8d3vrvz98sqV61m5cj1Tp7YYtwGKUv7o8tjbrq2trSxfvnzQ9ztUdXR00NLSUvUYQ8P6xXR1dTHSL1v05uujr4hYWkqZ3t963qEgKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjpiHl9rs2sN9RS4n4OnPm3Fv1OBrG6opbRFwQEasiYmdETG/UUBqZ2ts38TeX/R977VX1JMqg3jO3J4DzgQcaMItGsFIKs2bdzWETx/Cxsw+oehwlUFfcSimrSylPNWoYjVzXXvswS5Y8x23feTdjx3q1RPXzDYAq98QTG7j88p9x1VVncuLx3VWPoyT6jVtELAIOfYuHvlpK+cFAdxQRs4HZAK2trXR0dAx4yOw6OzurHqFSt922jG3bdrBo0RruX9jBr558HYC7715NxA6uuOL0iies1kh/feypfuNWSvlII3ZUSpkPzAdoa2srLS0tjdhsGiP5eIwduzelwKJFa/v8fe3aTSxf/rsRfWze5DHYfV7cUOWuvPJMSrmi5+elM7jo/IMAuPTS97N48cXVDqdhq96vgpwXEe3AqcCPImJhY8aSpPrU9YFCKeUu4K4GzSIB8G/fmMptd15S9Rga5nxbKikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYpJeMmKSXjJikl4yYppbriFhH/HBG/jojHI+KuiGhp1GCSVI96z9zuA44vpZwAPA1cXv9IklS/uuJWSvlpKaW7tvgwMLn+kSSpfo285vZp4CcN3J4k7bG9+lshIhYBh77FQ18tpfygts5XgW7gtl1sZzYwu7a4NSKe2P1x0zoY2Fj1EEOIx6Mvj0dfRw9kpSil1LWXiJgFfB74cCnljQE+55ellOl17TgRj0dfHo++PB59DfR49Hvm1s9OZgBfAc4YaNgkaTDUe81tHrAfcF9ErIiI6xswkyTVra4zt1LKUXv41Pn17Dchj0dfHo++PB59Deh41H3NTZKGIm+/kpRSZXHz1q2+IuKCiFgVETsjYkR+MhYRMyLiqYhYExGXVT1P1SJiQURs8GtTEBFTIuL+iFhd+//Jpf09p8ozN2/d6usJ4HzggaoHqUJENAHXAR8FjgUujIhjq52qcjcDM6oeYojoBuaWUt4NnAJ8ob/XR2Vx89atvkopq0spT1U9R4VOBtaUUp4ppWwD7gA+VvFMlSqlPAC8UvUcQ0Ep5cVSyrLa753AaqB1V88ZKtfcvHVLrcDzvZbb6efFq5EpIqYCbcAju1qvrq+CDGCIhty6lcVAjscIFm/xNz/KVx8R0QzcCcwppWza1bpva9xKKR/Z1eO1W7fOpefWrfQv5P6OxwjXDkzptTwZWFfRLBqCImI0PWG7rZTy/f7Wr/LT0jdv3fpzb90S8Bjwrog4IiLGAJ8A7ql4Jg0RERHATcDqUso1A3lOldfcvHWrl4g4LyLagVOBH0XEwqpnGky1D5cuARbSc7H4e6WUVdVOVa2IuB14CDg6Itoj4jNVz1Sh04FPAR+q9WJFRJyzqyd4h4KklIbKp6WS1FDGTVJKxk1SSsZNUkrGTVJKxk1SSsZNUkrGTVJK/w+3RUYyFNITOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dico_sten = {\n",
    "    'dim':2,\n",
    "    'schemes': [\n",
    "        {'velocities': list(range(5))}\n",
    "    ],\n",
    "}\n",
    "sten = pylbm.Stencil(dico_sten)\n",
    "print(sten)\n",
    "sten.visualize();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The domain\n",
    "\n",
    "In order to build the domain of the simulation, the dictionary should contain the space step $\\dx$ and the stencils of the velocities (one for each scheme). \n",
    "\n",
    "We construct a domain with $N=10$ points in space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+--------------------+\n",
      "| Domain information |\n",
      "+--------------------+\n",
      "    - spatial dimension: 2\n",
      "    - space step: 0.1\n",
      "    - with halo:\n",
      "        bounds of the box: [-0.05 -0.05] x [1.05 1.05]\n",
      "        number of points: [12, 12]\n",
      "    - without halo:\n",
      "        bounds of the box: [0.05 0.05] x [0.95 0.95]\n",
      "        number of points: [10, 10]\n",
      "        \n",
      "    +----------------------+\n",
      "    | Geometry information |\n",
      "    +----------------------+\n",
      "        - spatial dimension: 2\n",
      "        - bounds of the box: [0. 1.] x [0. 1.]\n",
      "        - labels: [0, 0, 0, 0]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARIAAAEICAYAAACTenveAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFGlJREFUeJzt3X+QXWV9x/H3JwmBCYlAsiAUWBYw2kYGQSNa/IEUGAFbsYNVqEigaApaOx21IwwMzeA4oo7FsWotIuKPkR8yrUaLI5Ef2hEDLiIJ0AZCAE1Z8hNpcPmRkG//OM/qzbK79+59zrn37N3Pa2Znzz3nPPf7zdnNd8859zzPo4jAzCzHjG4nYGZTnwuJmWVzITGzbC4kZpbNhcTMsrmQmFk2FxLrGkn9kp6WNLPbuVgeF5JpSNKjkp6RtE3SbyXdIel8SR39fYiIX0fE3Ih4oZNxrXwuJNPXX0TEPOAQ4HLgY8BXu5uSTVUuJNNcRDwVEcuBdwNLJB0haS9J35C0SdJjki4ZOVuRdI6kn0m6Ip3NrJN0bFr/G0kbJS0ZeX9Jb5N0j6T/S9uXNWwbkBSSZqXXt0v6eHr/bZJultTX4UNibXAhMQAi4i5gPfAm4F+AvYDDgOOAs4FzG3Z/HbAKWAB8G7gOeC3wMuAs4AuS5qZ9f5fa7w28DbhA0jsmSOWvU6z9gNnAR0v451nFXEis0ePAfIqzk4siYltEPAp8Fnhvw36PRMTX0r2N64GDgcsi4rmIuBl4nqKoEBG3R8TqiNgZEauAaymK03i+FhEPRsQzwA3AUSX/G60CLiTW6EBgFsWZwGMN6x9L20ZsaFh+BiAiRq+bCyDpdZJuS5dJTwHnAxNdrjzRsDw88j5Wby4kBoCk11IUi+8C2yluwo7oB/63zbf+NrAcODgi9gK+DCgjVashF5JpTtJLJP05xX2Ob0XEvRSXFJ+QNE/SIcCHgW+1GWIesDUinpV0DMU9EOsxs7qdgHXN9yXtAHYCDwD/THG2APAhihuu64Bnga8AV7cZ5wPAZyV9AfgJRZHaOyNvqyF5YCMzy+VLGzPL5kJiZtlcSMwsmwuJmWWr7ac2fX19MTAw0O00zHra3XffvTki9s19n9oWkoGBAQYHB7udhllPk/RY872a86WNmWVzITGzbKUUEklXp3Eo7htn+3skrUpfd0h6VRlxzaweyjojuQY4eYLtjwDHRcSRwMeBK0uKa2Y1UMrN1oj4qaSBCbbf0fByJXBQGXHNrB66cY/kPOCHXYhrZhXp6Me/ko6nKCRvHGf7UmApQH9/fwczM7McHTsjkXQkcBVwWkRsGWufiLgyIhZHxOJ9981+RsbMOqQjhURSP/DvwHsj4sFOxDSzzinl0kbStcBbgD5J64F/AnYDiIgvA5dSjDj+JUkAOyJicRmxzaz7yvrU5swm298HvK+MWGZWP36y1cyyuZCYWTYXEjPL5kJiZtlcSMwsmwuJmWVzITGzbC4kZpbNhcTMsrmQmFk2FxIzy+ZCYmbZXEjMLJsLiZllq+1Mex11z4Xw7BN/eL3H/nD05dW2dUzHLKNtTbiQQPFD3HPgD69/92j1bR3TMctoWxO+tDGzbC4kZpbNlzZQXJM2nk7usX/1bR3TMctoWxOKiG7nMKbFixfH4OBgt9Mw62mS7i5jIPZOTSIuSZ+XtDZNJP7qMuKaWT10ahLxU4CF6Wsp8K8lxTWzGiilkETET4GtE+xyGvCNKKwE9pZ0QBmxzaz7OvWpzYHAbxper0/resIQyxhiWcfaOWbvxZzqOlVINMa6F93llbRU0qCkwU2bNnUgLTMrQ6cKyXrg4IbXBwGPj97Jk4ibTU2dKiTLgbPTpzevB56KiKEOxTazinVqEvGbgFOBtcAwcG4Zcc2sHjo1iXgAHywjlpnVjx+Rh+xu3MPD23lmeDu7MUxf35zK2zlmj8X0MAI9opVu3MueH7Pp8PAO2LaYORIr3vQwJ510eEu/RJs3D7NuaCuSWL269XY5bR2z2pg/WPccQvzkrm0cfvhuzJmz25j7Ljtu1AoPI2DDw9uRxOzZM5kxQ2zZMtxSuy1bhpHE3LmzJ9Uup61jVhyT4vdAEsPD21uO2Qt8RtKqZbPHXr95B3euOJwZM8TOncGCBa395VuwYA6/XHEBWybZLqetY1Ybc99fns2MGeL0+Rs46bi59PW1HHbKc+9fyL5G3bx5mC1bhlmwYM6kro3bbeeYPRazi/dIyur960JiNo3VahgBM5veXEjMLJsLiZllcyExs2wuJGaWzYXEzLK5kJhZNhcSM8vmR+TBT7Y6Zndjuvdvj8js/Tv0cNFbdMXpQ5PqLbpixcO/76Mz2R6q7bR1zGpj/u1125BEhHv/2iSN9P5tp7fojBnipS+d21YP1XbaOma1MUd+D9z718Y3Qe/f1Su2ttX7d+fOYMOGp9vqodpOW8esNubp89OZzOxw79+6cO9fx5w2Md37tzru/WtWvVr1/pV0sqQ1aZLwC8fY3i/pNkn3pEnETy0jrpnVQ3YhkTQT+CLFROGLgDMlLRq12yXADRFxNHAG8KXcuGZWH2WckRwDrI2IdRHxPHAdxaThjQJ4SVreizFm2TOzqauMT23GmiD8daP2WQbcLOlDwJ7AiSXENbOaKOOMpJUJws8EromIgyhm3PumpBfF9iTiZlNTGYWklQnCzwNuAIiInwN7AC/6lN2TiJtNTWUUkl8ACyUdKmk2xc3U5aP2+TVwAoCkP6EoJD7lMOsR2YUkInYAfwf8CPhvik9n7pd0maS3p90+Arxf0r3AtcA5UdcHWMxs0sqaRPwm4KZR6y5tWH4AeEMZscysftzXBvyIvGN2N6aHEegRHkbAMUuI6WEErG0eRsAxR9p5GAFrzsMIOGaTdh5GoIY8jIBjTpuYHkagOh5GwKx6tRpGwMymNxcSM8vmQmJm2VxIzCybC4mZZXMhMbNsLiRmls2FxMyy+RF58JOtjtndmO792yPc+9cxS4jp3r/WNvf+dcyRdu79a825969jNmnn3r815N6/jjltYrr3b3Xc+9eserXq/dtsEvG0z7skPSDpfknfLiOumdVD9j2ShknET6KYLOsXkpankeNH9lkIXAS8ISKelLRfblwzq49OTSL+fuCLEfEkQERsLCGumdVEGYVkrEnEDxy1z8uBl0v6maSVkk4uIa6Z1UQZH/+2Mon4LGAh8BaKuYH/S9IREfHbXd5IWgosBejv7y8hNTPrhE5NIr4e+F5EbI+IR4A1FIVlF55E3Gxq6tQk4t8FjgeQ1EdxqbOuhNhmVgOdmkT8R8AWSQ8AtwH/GBFbcmObWT34gTSzaaysB9Lc1wb8iLxjdjemhxHoER5GwDFLiOlhBKxtHkbAMUfaeRgBa87DCDhmk3YeRqCGPIyAY06bmB5GoDr+1MaserUaRsDMpjcXEjPL5kJiZtlcSMwsmwuJmWVzITGzbC4kZpbNhcTMsvkRefCTrY7Z3Zju/dsj3PvXMUuI6d6/1jb3/nXMkXbu/WvNufevYzZp596/NeTev445bWK692913PvXrHq16v3byiTiab93SgpJ2YmbWX1kF5KGScRPARYBZ0paNMZ+84C/B+7MjWlm9dKpScQBPg58Gni2hJhmViMdmURc0tHAwRHxgxLimVnNlFFIJpxEXNIM4ArgI03fSFoqaVDS4KZNm0pIzcw6oROTiM8DjgBul/Qo8Hpg+Vg3XD2JuNnUVPkk4hHxVET0RcRARAwAK4G3R4Q/2zXrEZ2aRNzMelgpj8hHxE3ATaPWXTrOvm8pI6aZ1Yf72oAfkXfM7sb0MAI9wsMIOGYJMT2MgLXNwwg45kg7DyNgzXkYAcds0s7DCNSQhxFwzGkT08MIVMfDCJhVr1bDCJjZ9OZCYmbZXEjMLJsLiZllcyExs2wuJGaWzYXEzLK5kJhZNj8iD36y1TG7G9O9f3uEe/86Zgkx3fvX2ubev4450s69f6059/51zCbt3Pu3htz71zGnTUz3/q2Oe/+aVa9WvX+bTSIu6cOSHpC0StItkg4pI66Z1UOnJhG/B1gcEUcCN1LMAWxmPaIjk4hHxG0RMXL7eyXFbHxm1iM6Mon4KOcBPywhrpnVRBkf/044ifguO0pnAYuB0Y/kjGxfCiwF6O/vLyE1M+uETkwiDoCkE4GLKeb9fW6sN/Ik4mZTU+WTiANIOhr4N4oisrGEmGZWI52aRPwzwFzgO5J+JWn5OG9nZlNQRyYRj4gTy4hjZvXkvjbgR+Qds7sxPYxAj/AwAo5ZQkwPI2Bt8zACjjnSzsMIWHMeRsAxm7TzMAI15GEEHHPaxPQwAtXxMAJm1avVMAJmNr25kJhZNhcSM8vmQmJm2VxIzCybC4mZZXMhMbNsLiRmls2PyIOfbHXM7sZ0798ekdn7997XXIIEv3zh3En1Fv3x0EeQxM4VF0y6h2o7bR2z2pifX3cpIO5f8VH3/rXJKXr/wuzZM9vqLbrnJHsN57R1zGpjgpg9e6Z7/9oEJuj9u3HomeIv2OrJ9RaNoeB3Tz/fVg/Vdto6ZrUxFz3/NJI4fP6Gadf714UEimvSxtPJPfZvuWlf3xwOmzOf4eHtHLF/66fQ7bZzzN6LmfP7VxcuJJB9Y2vOnOJ6uI/J3dRrt51j9ljMKXZjdSydmkR8d0nXp+13ShooI66Z1UOnJhE/D3gyIl4GXAF8KjeumdVHRyYRT6+/npZvBE6QNNZUn2Y2BXVqEvHf75Mm1HoKWFBCbDOrgTIKSSuTiLc00bikpZIGJQ1u2rSphNTMrBOyx2yV9KfAsoh4a3p9EUBEfLJhnx+lfX4uaRbwBLBvTBDcY7aaVa9OY7Y2nUQ8vV6Slt8J3DpRETGzqSX7OZKI2CFpZBLxmcDVI5OIA4MRsRz4KvBNSWuBrRTFxsx6RKcmEX8W+KsyYplZ/fjJVsjrxt1uW8d0zDLa1oQLCeR14263rWM6Zhlta8LDCJhZNhcSM8vmSxvI68bdblvHdMwy2taEJxE3m8bq9ECamU1zLiRmls2FxMyyuZCYWTYXEjPL5kJiZtlcSMwsmwuJmWVzITGzbC4kZpbNhcTMsrmQmFk2FxIzy+ZCYmbZsgqJpPmSVkh6KH3fZ4x9jpL0c0n3S1ol6d05Mc2sfnLPSC4EbomIhcAt6fVow8DZEfFK4GTgc5L2zoxrZjWSW0gaJwf/OvCO0TtExIMR8VBafhzYCOybGdfMaiS3kLw0IoYA0vf9JtpZ0jHAbODhzLhmViNNx2yV9GNgrEEkL55MIEkHAN8ElkTEznH2WQosBejv75/M25tZFzUtJBFx4njbJG2QdEBEDKVCsXGc/V4C/CdwSUSsnCDWlcCVUIzZ2iw3M6uH3EubxsnBlwDfG71Dmlj8P4BvRMR3MuOZWQ3lFpLLgZMkPQSclF4jabGkq9I+7wLeDJwj6Vfp66jMuGZWI56Owmwa83QUZlYbLiRmlq22lzaStgFrup3HKH3A5m4n0cD5TKxu+UD9cnpFRMzLfZM6z/27poxrtzJJGqxTTs5nYnXLB+qXk6RSbkT60sbMsrmQmFm2OheSK7udwBjqlpPzmVjd8oH65VRKPrW92WpmU0edz0jMbIpwITGzbF0tJLlDNUq6RtIjuX14JJ0saY2ktZJeNMqbpN0lXZ+23ylpoGHbRWn9GklvbSd+G/l8WNID6XjcIumQhm0vNByP5WXk02JO50ja1BD7fQ3blqSf8UOSloxuW1E+VzTk8qCk3zZsK/0YSbpa0kZJ942zXZI+n/JdJenVDduqOD7N8nlPymOVpDskvaph26OSVqfj09rHwxHRtS/g08CFaflC4FNj7PNyYGFa/iNgCNg7vb4GeGdmDjMpBlo6jGLQpXuBRaP2+QDw5bR8BnB9Wl6U9t8dODS9z8wO5HM8MCctXzCST3r9dAU/p1ZyOgf4whht5wPr0vd90vI+Veczav8PAVdXfIzeDLwauG+c7acCPwQEvB64s6rj02I+x47EAU4ZySe9fhTom0y8bl/a1GGoxmOAtRGxLiKeB65LeY2X543ACZKU1l8XEc9FxCPA2vR+leYTEbdFxHB6uRI4KDNmdk4TeCuwIiK2RsSTwAqKsXs7mc+ZwLWZMScUET8Ftk6wy2kUQ2lEFGPy7J3G8Kni+DTNJyLuSPGghN+hbheSMoZq/EQ6PbtC0u5t5HAg8JuG1+vTujH3iYgdwFPAghbbVpFPo/Mo/tKN2EPSoKSVkl5UmCvO6fT0s7hR0sGTbFtFPqTLvkOBWxtWV3GMmhkv5yqOz2SN/h0K4GZJd6sYtbCpyh+RV7VDNV4EPEFRXK4EPgZcNtkUx1g3+jPx8fZppe1ktfyeks4CFgPHNazuj4jHJR0G3CppdUTkjpHbSk7fB66NiOcknU9xBvdnLbatIp8RZwA3RsQLDeuqOEbNdPJ3qGWSjqcoJG9sWP2GdHz2A1ZI+p90hjOuys9IIuLEiDhijK/vARtSgRgpFJMaqjEihtKp4nPA12jvsmI9cHDD64OAx8fbR9IsYC+K08ZW2laRD5JOpCjGb0//fuD3l39ExDrgduDozHxayikitjTk8RXgNa22rSKfBmcw6rKmomPUzHg5V3F8WiLpSOAq4LSI2DKyvuH4bKQY3bD5/6uybzpN8obQZ9j1Zuunx9hnNsWcOf8wxrYD0ncBnwMubyOHWRQ3uA7lDzfuXjlqnw+y683WG9LyK9n1Zus68m+2tpLP0RSXdwtHrd8H2D0t9wEPMcFNyJJzOqBh+S+BlWl5PvBIym2ftDy/6nzSfq+guHGoqo9Rer8Bxr+5+TZ2vdl6V1XHp8V8+inu6R07av2ewLyG5TuAk5vGKiPhjH/oglQkHkrf56f1i4Gr0vJZwHbgVw1fR6VttwKrgfuAbwFz28zjVODB9J/z4rTuMoq/9gB7AN9JB/4u4LCGthendmuAU0o6Ls3y+TGwoeF4LE/rj03H4970/bwSf1bNcvokcH+KfRvwxw1t/yYdu7XAuZ3IJ71exqg/LlUdI4qznqH0u7qe4nLhfOD8tF3AF1O+q4HFFR+fZvlcBTzZ8Ds0mNYflo7NvenneXEr8fyIvJll6/anNmbWA1xIzCybC4mZZXMhMbNsLiRmls2FxMyyuZCYWbb/B+vFa7F9tVEIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 10\n",
    "dx = (xmax-xmin)/N\n",
    "dico_dom = {\n",
    "    'box': {'x': [xmin, xmax], \n",
    "            'y': [ymin, ymax], \n",
    "            'label': 0\n",
    "           },\n",
    "    'space_step': dx,\n",
    "    'schemes': [\n",
    "        {'velocities': list(range(5)),}\n",
    "    ],\n",
    "}\n",
    "dom = pylbm.Domain(dico_dom)\n",
    "print(dom)\n",
    "dom.visualize(view_distance=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The scheme\n",
    "\n",
    "In pylbm, a simulation can be performed by using several coupled schemes. In this example, a single scheme is used and defined through a list of one single dictionary. This dictionary should contain:\n",
    "\n",
    "* 'velocities': a list of the velocities\n",
    "* 'conserved_moments': a list of the conserved moments as sympy variables\n",
    "* 'polynomials': a list of the polynomials that define the moments\n",
    "* 'equilibrium': a list of the equilibrium value of all the moments\n",
    "* 'relaxation_parameters': a list of the relaxation parameters ($0$ for the conserved moments)\n",
    "* 'init': a dictionary to initialize the conserved moments\n",
    "\n",
    "(see the documentation for more details)\n",
    "\n",
    "The scheme velocity could be taken to $1/\\dx$ and the inital value of $u$ to \n",
    "\n",
    "$$ u(t=0,x) = \\sin(\\pi x)\\sin(\\pi y).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+--------------------+\n",
      "| Scheme information |\n",
      "+--------------------+\n",
      "    - spatial dimension: 2\n",
      "    - number of schemes: 1\n",
      "    - number of velocities: 5\n",
      "    - conserved moments: [u]\n",
      "        \n",
      "    +----------+\n",
      "    | Scheme 0 |\n",
      "    +----------+\n",
      "        - velocities\n",
      "            (0: 0, 0)\n",
      "            (1: 1, 0)\n",
      "            (2: 0, 1)\n",
      "            (3: -1, 0)\n",
      "            (4: 0, -1)\n",
      "\n",
      "        - polynomials\n",
      "                    \n",
      "            ⎡   1   ⎤\n",
      "            ⎢       ⎥\n",
      "            ⎢  X    ⎥\n",
      "            ⎢  ──   ⎥\n",
      "            ⎢  LA   ⎥\n",
      "            ⎢       ⎥\n",
      "            ⎢  Y    ⎥\n",
      "            ⎢  ──   ⎥\n",
      "            ⎢  LA   ⎥\n",
      "            ⎢       ⎥\n",
      "            ⎢ 2    2⎥\n",
      "            ⎢X  + Y ⎥\n",
      "            ⎢───────⎥\n",
      "            ⎢     2 ⎥\n",
      "            ⎢ 2⋅LA  ⎥\n",
      "            ⎢       ⎥\n",
      "            ⎢ 2    2⎥\n",
      "            ⎢X  - Y ⎥\n",
      "            ⎢───────⎥\n",
      "            ⎢     2 ⎥\n",
      "            ⎣ 2⋅LA  ⎦\n",
      "\n",
      "        - equilibria\n",
      "                    \n",
      "            ⎡  u  ⎤\n",
      "            ⎢     ⎥\n",
      "            ⎢ 0.0 ⎥\n",
      "            ⎢     ⎥\n",
      "            ⎢ 0.0 ⎥\n",
      "            ⎢     ⎥\n",
      "            ⎢0.5⋅u⎥\n",
      "            ⎢     ⎥\n",
      "            ⎣ 0.0 ⎦\n",
      "\n",
      "        - relaxation parameters\n",
      "                    \n",
      "            ⎡0.0⎤\n",
      "            ⎢   ⎥\n",
      "            ⎢0.4⎥\n",
      "            ⎢   ⎥\n",
      "            ⎢0.4⎥\n",
      "            ⎢   ⎥\n",
      "            ⎢1.0⎥\n",
      "            ⎢   ⎥\n",
      "            ⎣1.0⎦\n",
      "\n",
      "    - moments matrices\n",
      "                \n",
      "        ⎡1   1    1     1     1  ⎤\n",
      "        ⎢                        ⎥\n",
      "        ⎢   10         -10       ⎥\n",
      "        ⎢0  ──    0    ────   0  ⎥\n",
      "        ⎢   LA          LA       ⎥\n",
      "        ⎢                        ⎥\n",
      "        ⎢         10         -10 ⎥\n",
      "        ⎢0   0    ──    0    ────⎥\n",
      "        ⎢         LA          LA ⎥\n",
      "        ⎢                        ⎥\n",
      "        ⎢    50   50    50    50 ⎥\n",
      "        ⎢0  ───  ───   ───   ─── ⎥\n",
      "        ⎢     2    2     2     2 ⎥\n",
      "        ⎢   LA   LA    LA    LA  ⎥\n",
      "        ⎢                        ⎥\n",
      "        ⎢    50  -50    50   -50 ⎥\n",
      "        ⎢0  ───  ────  ───   ────⎥\n",
      "        ⎢     2    2     2     2 ⎥\n",
      "        ⎣   LA   LA    LA    LA  ⎦\n",
      "\n",
      "    - inverse of moments matrices\n",
      "                \n",
      "        ⎡                  2        ⎤\n",
      "        ⎢               -LA         ⎥\n",
      "        ⎢1   0     0    ─────    0  ⎥\n",
      "        ⎢                 50        ⎥\n",
      "        ⎢                           ⎥\n",
      "        ⎢                  2      2 ⎥\n",
      "        ⎢    LA          LA     LA  ⎥\n",
      "        ⎢0   ──    0     ───    ─── ⎥\n",
      "        ⎢    20          200    200 ⎥\n",
      "        ⎢                           ⎥\n",
      "        ⎢                  2      2 ⎥\n",
      "        ⎢          LA    LA    -LA  ⎥\n",
      "        ⎢0   0     ──    ───   ─────⎥\n",
      "        ⎢          20    200    200 ⎥\n",
      "        ⎢                           ⎥\n",
      "        ⎢                  2      2 ⎥\n",
      "        ⎢   -LA          LA     LA  ⎥\n",
      "        ⎢0  ────   0     ───    ─── ⎥\n",
      "        ⎢    20          200    200 ⎥\n",
      "        ⎢                           ⎥\n",
      "        ⎢                  2      2 ⎥\n",
      "        ⎢         -LA    LA    -LA  ⎥\n",
      "        ⎢0   0    ────   ───   ─────⎥\n",
      "        ⎣          20    200    200 ⎦\n",
      "\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "import sympy as sp\n",
    "\n",
    "def solution(x, y, t):\n",
    "    return np.sin(np.pi*x)*np.sin(np.pi*y)*np.exp(-2*np.pi**2*mu*t)\n",
    "\n",
    "# parameters\n",
    "mu = 1.\n",
    "la = 1./dx\n",
    "s1 = 2./(1+4*mu)\n",
    "s2 = 1.\n",
    "u, X, Y, LA = sp.symbols('u, X, Y, LA')\n",
    "\n",
    "dico_sch = {\n",
    "    'dim': 2,\n",
    "    'scheme_velocity': la,\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': list(range(5)),\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, X/LA, Y/LA, (X**2+Y**2)/(2*LA**2), (X**2-Y**2)/(2*LA**2)],\n",
    "            'equilibrium': [u, 0., 0., .5*u, 0.],\n",
    "            'relaxation_parameters': [0., s1, s1, s2, s2],\n",
    "        }\n",
    "    ],\n",
    "    'parameters': {LA: la},\n",
    "}\n",
    "\n",
    "sch = pylbm.Scheme(dico_sch)\n",
    "print(sch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The simulation\n",
    "\n",
    "A simulation is built by defining a correct dictionary.\n",
    "\n",
    "We combine the previous dictionaries to build a simulation. In order to impose the homogeneous Dirichlet conditions in $x=0$, $x=1$, $y=0$, and $y=1$, the dictionary should contain the key 'boundary_conditions' (we use pylbm.bc.Anti_bounce_back function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+------------------------+\n",
      "| Simulation information |\n",
      "+------------------------+\n",
      "        \n",
      "    +--------------------+\n",
      "    | Domain information |\n",
      "    +--------------------+\n",
      "        - spatial dimension: 2\n",
      "        - space step: 0.1\n",
      "        - with halo:\n",
      "            bounds of the box: [-0.05 -0.05] x [1.05 1.05]\n",
      "            number of points: [12, 12]\n",
      "        - without halo:\n",
      "            bounds of the box: [0.05 0.05] x [0.95 0.95]\n",
      "            number of points: [10, 10]\n",
      "            \n",
      "        +----------------------+\n",
      "        | Geometry information |\n",
      "        +----------------------+\n",
      "            - spatial dimension: 2\n",
      "            - bounds of the box: [0. 1.] x [0. 1.]\n",
      "            - labels: [0, 0, 0, 0]\n",
      "        \n",
      "    +--------------------+\n",
      "    | Scheme information |\n",
      "    +--------------------+\n",
      "        - spatial dimension: 2\n",
      "        - number of schemes: 1\n",
      "        - number of velocities: 5\n",
      "        - conserved moments: [u]\n",
      "            \n",
      "        +----------+\n",
      "        | Scheme 0 |\n",
      "        +----------+\n",
      "            - velocities\n",
      "                (0: 0, 0)\n",
      "                (1: 1, 0)\n",
      "                (2: 0, 1)\n",
      "                (3: -1, 0)\n",
      "                (4: 0, -1)\n",
      "\n",
      "            - polynomials\n",
      "                        \n",
      "                ⎡   1   ⎤\n",
      "                ⎢       ⎥\n",
      "                ⎢  X    ⎥\n",
      "                ⎢  ──   ⎥\n",
      "                ⎢  LA   ⎥\n",
      "                ⎢       ⎥\n",
      "                ⎢  Y    ⎥\n",
      "                ⎢  ──   ⎥\n",
      "                ⎢  LA   ⎥\n",
      "                ⎢       ⎥\n",
      "                ⎢ 2    2⎥\n",
      "                ⎢X  + Y ⎥\n",
      "                ⎢───────⎥\n",
      "                ⎢     2 ⎥\n",
      "                ⎢ 2⋅LA  ⎥\n",
      "                ⎢       ⎥\n",
      "                ⎢ 2    2⎥\n",
      "                ⎢X  - Y ⎥\n",
      "                ⎢───────⎥\n",
      "                ⎢     2 ⎥\n",
      "                ⎣ 2⋅LA  ⎦\n",
      "\n",
      "            - equilibria\n",
      "                        \n",
      "                ⎡  u  ⎤\n",
      "                ⎢     ⎥\n",
      "                ⎢ 0.0 ⎥\n",
      "                ⎢     ⎥\n",
      "                ⎢ 0.0 ⎥\n",
      "                ⎢     ⎥\n",
      "                ⎢0.5⋅u⎥\n",
      "                ⎢     ⎥\n",
      "                ⎣ 0.0 ⎦\n",
      "\n",
      "            - relaxation parameters\n",
      "                        \n",
      "                ⎡0.0⎤\n",
      "                ⎢   ⎥\n",
      "                ⎢0.4⎥\n",
      "                ⎢   ⎥\n",
      "                ⎢0.4⎥\n",
      "                ⎢   ⎥\n",
      "                ⎢1.0⎥\n",
      "                ⎢   ⎥\n",
      "                ⎣1.0⎦\n",
      "\n",
      "        - moments matrices\n",
      "                    \n",
      "            ⎡1   1    1     1     1  ⎤\n",
      "            ⎢                        ⎥\n",
      "            ⎢   10         -10       ⎥\n",
      "            ⎢0  ──    0    ────   0  ⎥\n",
      "            ⎢   LA          LA       ⎥\n",
      "            ⎢                        ⎥\n",
      "            ⎢         10         -10 ⎥\n",
      "            ⎢0   0    ──    0    ────⎥\n",
      "            ⎢         LA          LA ⎥\n",
      "            ⎢                        ⎥\n",
      "            ⎢    50   50    50    50 ⎥\n",
      "            ⎢0  ───  ───   ───   ─── ⎥\n",
      "            ⎢     2    2     2     2 ⎥\n",
      "            ⎢   LA   LA    LA    LA  ⎥\n",
      "            ⎢                        ⎥\n",
      "            ⎢    50  -50    50   -50 ⎥\n",
      "            ⎢0  ───  ────  ───   ────⎥\n",
      "            ⎢     2    2     2     2 ⎥\n",
      "            ⎣   LA   LA    LA    LA  ⎦\n",
      "\n",
      "        - inverse of moments matrices\n",
      "                    \n",
      "            ⎡                  2        ⎤\n",
      "            ⎢               -LA         ⎥\n",
      "            ⎢1   0     0    ─────    0  ⎥\n",
      "            ⎢                 50        ⎥\n",
      "            ⎢                           ⎥\n",
      "            ⎢                  2      2 ⎥\n",
      "            ⎢    LA          LA     LA  ⎥\n",
      "            ⎢0   ──    0     ───    ─── ⎥\n",
      "            ⎢    20          200    200 ⎥\n",
      "            ⎢                           ⎥\n",
      "            ⎢                  2      2 ⎥\n",
      "            ⎢          LA    LA    -LA  ⎥\n",
      "            ⎢0   0     ──    ───   ─────⎥\n",
      "            ⎢          20    200    200 ⎥\n",
      "            ⎢                           ⎥\n",
      "            ⎢                  2      2 ⎥\n",
      "            ⎢   -LA          LA     LA  ⎥\n",
      "            ⎢0  ────   0     ───    ─── ⎥\n",
      "            ⎢    20          200    200 ⎥\n",
      "            ⎢                           ⎥\n",
      "            ⎢                  2      2 ⎥\n",
      "            ⎢         -LA    LA    -LA  ⎥\n",
      "            ⎢0   0    ────   ───   ─────⎥\n",
      "            ⎣          20    200    200 ⎦\n",
      "\n",
      "        \n"
     ]
    }
   ],
   "source": [
    "dico = {\n",
    "    'box': {'x': [xmin, xmax],\n",
    "            'y': [ymin, ymax],\n",
    "            'label': 0},\n",
    "    'space_step': dx,\n",
    "    'scheme_velocity': la,\n",
    "    'schemes': [\n",
    "        {\n",
    "            'velocities': list(range(5)),\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, X/LA, Y/LA, (X**2+Y**2)/(2*LA**2), (X**2-Y**2)/(2*LA**2)],\n",
    "            'equilibrium': [u, 0., 0., .5*u, 0.],\n",
    "            'relaxation_parameters': [0., s1, s1, s2, s2],\n",
    "        }\n",
    "    ],\n",
    "    'init': {u: (solution, (0.,))},\n",
    "    'boundary_conditions': {\n",
    "        0: {'method': {0: pylbm.bc.AntiBounceBack,}},\n",
    "    },\n",
    "    'parameters': {LA: la},\n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico)\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run a simulation\n",
    "\n",
    "Once the simulation is initialized, one time step can be performed by using the function one_time_step.\n",
    "\n",
    "We compute the solution of the heat equation at $t=0.1$. On the same graphic, we plot the initial condition, the exact solution and the numerical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAEICAYAAABs2F48AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXu0bFtd3/n5zblqn8tLQcHI+xrpGB9DsQMYTVpRiRANahzaoqBixAxtHeqwEzVGW3Rox9Z0BxUCvlGMooBG7YAYldu+RURAjTpUuAj44vIQvJezd9Wcv/7j95tzzbWqau/aZ9c5u85lfs/YZ60116pVa8361q++v8ecS1SVjo6Ojo7DQbjsC+jo6OjomKIb5o6Ojo4DQzfMHR0dHQeGbpg7Ojo6DgzdMHd0dHQcGLph7ujo6DgwdMPc0TGDiNwuIrdf9nV0vPuiG+brABFRETm1QNy//Coit96Yq6rv+1h/36ffyPc9JIjIbWd9Ph0dl4nhsi+go+MA8fGXfQEd797ohrmjYwZV/bPLvoaOd2/0UMYBQkT+oYg8V0TeICLHIvLXIvKjIvIBG479ByLyrSLyChF5sx//ehH5HhF5yOzY5wIv881vKCEX/3vsvq/Nj3+EiLxARN4mIneKyK+LyCeJyFP9fZ86O15F5LYt53rupvCPn+tFIvJaEXmXiLxDRH5NRJ4yO+5WD2F8TPNeOn/PbTFmEbkiIl8jIq8Rkbv8fX5FRP7XDcfe6ud9rq8/X0TuEJGr/ln9i40d3NFBV8wHBxF5AvCTwAL4WeBPgYcAnwZ8koh8rKq+snnJpwFfhBncXwdOgA8GngY8UUQepapv8mP/qy8/D/j/gNua89y+72sTkf8J+A3gvYGXAK8CHuHX8ZKze2NnPBv4H8AvA3/p7/eJwPNE5ANU9ev9uLcD3wg8FXi4rxfcftobiMgR8FLMqP8R8CzgnsCnAz8uIo9U1a/d8NKHAy8HXgs8D3gv4DOBnxaRx6nqyza8puPdHara//b8B6j/Pf2Uv7f7Mbc2r7sf8DbgDuCDZuf8YODvgFfO2h8MXNlwDZ8AJODZs/bHlms75z1dy7X9vL/Xl8/aP6Xpo6du6LvbtlzDc+d95u3vv+HYI+AXgSXw4Nm+24z6W+/1duD2Wdu/8/d+MTA07e/jxyvwUU37rc09fsPsXI8v57psrva/w/y79Au4O/41X8hd/m5tXvfl3vYlW877n3z/B+14Ha8BXjtru1bDfK5rw5S0Ykoxbjj+tn0Z5lOu+dP8+M/d9N6nvG6TYf4TIAP/cMPxX+Dv8wNNWzHMt2+5/9cDd1w2V/vfYf71UMZ1hKrKtn0ew3z4rPkjfflhW8rZ/oEvPxBz3RERAZ6Muecfhinb2Lzm5JyXvQ3nvbYP9+1fVdW04fjb8FjvRSEiDwO+GqumeBhwj9khD77g+e+DhWDepKp/tOGQX/Llh2/Y96ot9/8Gxj7t6JigG+bDwnv78gvPOO7ezfr/A3wFFlt9KfAm4F2+76msG/8bdW3v6cu/3nLcX134igAR+ftYDPd+wK9g4ZO/xcI4t2Lx9CsXfJtyL3+5ZX9pv++GfW/f8poVPfnesQXdMB8W/taXH6aqrznrYBF5H+DLgN/H4pvvnO3/rMu6tub4v7dl//tuaVe283KT4ftK7Efj81X1ue0Ov//PO/0yd0K5l23X/MDZcR0dF0L/xT4s/KYv/5cdj//72Gf48xuM8kN8/xzFrY4b9u3z2n7Xl/9URDa912O3vO5twEPnjX6OR244/hG+fNGGfdtCJak555nwvv0z4MFeaTLHx/rylRv2dXScG90wHxZ+EHN9v0FEHjPfKSJhVm98uy8nxk9E7g18L5uV51t8+bDreW2q+kbgvwPvB3zp7NhPYbvRfDnwMBH5hFn717E5LHO7Lx/bNorI47GSwU24lj74AUCAb5/19f2Br2+O6ei4MG6KUIaI/AFWDXDbRY4TkZcAz1fVH9rhPW8Hnqaqv3DuC75GqOpbROTTgZ8CflNEfhH4A6wa4GFYsui9gVv8+L8SkecDTwJeJSI/j8VD/xlwFasbnqvMP8bi0E8SkRPgz7HwwfNU9fX7ujbHl2B1zM9wQ/tqTOH+S6wO+okb3uo/YuVkPy0iPw68FfgozMDfxrrS/s/A5wMvEJEX+b19CPAE4CewmuE5fhH4DOAnReTFWEz+9ar6vG3379f1z7FSv1f76+7p53kf4NtU9VdPeX1Hx+647LKQ6/WH1Qr/yAVefzvwuGt8rXJKOVZz/o2lX1jS6plYidZV4B3YoIbnAZ86O/aewLdggz2uYtn+Z2FG8rZN1wE8GjNOf4sZVgUeu+O97XxtfvwjgBdiavtOzFB/EpaYXCuX89d8MvAKP/9bgOdjavm5m/oMM9y/hIVB3gn8KvCpbCkNxMI4/ydWyrdkVqLHhnI5b78F+Fospv+u5r0+a0s/KfDcLf248bO5zt+Jp2JVMjfsPfvftf2Jf2B3O3hJ1yNU9SlnHbvl9bdzgxXzuxN8KPYPsiFp13F94H3+NFX9p5d9LR2n46aIMfvcBY8TkaeLyE+IyA+LyDtF5A9E5FEbjnsCpmw+U0T+TkRe7ftvE5Gn+fr7i8gvichbfA6D/yIim7L+HR3XDSLyIJ/n480i8joR+TJvf7GI/N/NcT8uIj/g66dyV0QeKiI/6ed8i4g8U0Q+EHgO8JH+ndhWxtdxALgpDPMMn4y5tvcFfgZzqydQ1Z/DXNUfV9V7q+qHbTiPAP8BeBA2KOKhWPijo+OGQEQCFmt/NTYI5uOBr/DE5b8CPkdEPk5EnoyFn768vJQt3PXE5P+LjSy81c/7fFX9Q2xOld/w70QXIQeMm9Ew/6qqvlhtNNXzsNFu54aq/qmq/ndVPVbVN2MDNfYyEq2jY0c8GniAqn6Tqp6o6muxaponqepfYYb0h4DvwIaVvxPO5O5jMIP9b1X1TlW9qj0pedPhpqjKmKEdMXYXcIuIDKq6Os9JfHDGd2J1uffBfqTetrer7DgVHld+7iVfxmXj4cCDZmGFiI1gBFO+zwT+uDWuZ3D3oViFybm+Dx2HhZtRMe+Ks7Ka/8GP+VBVfQ/gKZiL2NFxo/AG4HWqet/m7z6q+om+/1uAPwQeOBvFeRp334DVgW8SXXfPTP/dEHdnw/zXwK0ex9uE+2BTVb5dRB4M/NsbdmUdHYaXA+8Qka8WkXuISBSRDxGRR4vIR2P12Z/rf9/lPIXTuftybO6ObxWRe4nILSLyT3zfXwMP8bmlOw4Yd2fD/AJfvkVENg2V/Ubgf8Zqef8bNgF8R8cNg+dJnogNAnodNtf192Fzb/ww8KWq+iYPY3w/8IM+m+BW7jbnfAQ2eOiNjINsfgkbFPRXInLHdb/BjmvG3baOuaOjo+Nmxd1ZMXd0dHTclOiGuaOjo+PA0A1zR0dHx4GhG+aOjo6OA0M3zB0dHR0HhlNH/h3JFb2Fe92oa7k27GNIyE1QmPJO3naHqj6gbD/+Y++lb3nrpmd8nh+/85rjl6rqE/ZyspsEnduHg5bb++Q13LzcPtUw38K9+Aj5eNuQGUuacRsStu8DYLZf5ucKW4T7/Lg5tr3uWpDz6fu3lRXOXrdWfphn2zo7vt0/29e+5y/oCyeT2N/x1hW//nMXevhzxS0Pet3993Kimwid2w0OiNv75DXcvNzeba4MkUpICb7uhBSRkUTNurU35GteU0lZSN6eq90PU4KufYH2OIJ6Trp2uyWot1eSFvI58aRtb8mYtXlNrueRnKfn0twQOm/90mSU4z4dwsXRub3Wfpnc7rw2nG2YnbhVOcTozU7UlrCxkNiJ2pI0yJScLeEnZG1ft05mnXwhtpD3NFJvUwdNu0x+6XV9mXX0MlXHfQ0p1Y+rBFa182qGGNCU67WKKuSMBoCI+PNSNQe2EViBY/bn8r1bonN7fXnJ3O68NpxvdjkJaypCYhzVQoyjAilklSnpKWpDpiTVVnFMSL6BrL7UTSTdxQPc4NlJQ0715Xi8TslbFcGGdm8TgJSAaOdsyawN+UXQlCAEa8sZlbDu+s2gwLKP2twfOrcPgtud14adDLMEgRiNjDFaRxdiDcNI1hBGpRHDSEQRtN2OI/lVBKJMtwPj68qnHGRK1rJaiXz+m5fy+VcCtvu0unKiNOT0fYXkSes2ACnXY6VdVx335YzmXMksKxndQREEUBWEhG4RD1mVq53AF0bnNgfF7c5rw44x5hIvCyNxi5qQYG5eaSvrIZhSKIQd4pSs0pAxNkQtJBXGtnJcaJSEL+bb54J//jIjr6iOP+qqDXnLvqYtGIkrgUMYyS1pJG82YqIKyQiKqJEzRhBFUkJD8GOS9/tmy6wIy2v5xnZM0bl9UNzuvDbsEGOextAmxA3RiFvcvCFWsmpZdzJrbBTGEKZkDTIlqpN0dAGpxwITQtfPcP5Zbvps5z/ElaxlW8f1SlA/JutI6obQkosb59viCQ9V6ytXExoUWdkx9cvobp+9hbmGkpJ9GUQsLieBTdJCgasaN9xkx87o3LZjDojbndeG3WPMNRkilbiE0f2ztjC6dUVFxClhESEPYaoa4qgkNExJPSGn0BzXtLMhJrcDedfURJaqHFCZkLuQtKxrURjJ2tRVRyiEVQXJExIrjG5gHN1JcxUjkOx9S6wzbU+CZISrejM+gOYA0bl9MNzuvDbsFmN20kmMNe5W1cRiGF27si6CLqLF00JAo8XhClFzDFUpaHT1UJREpJLUflVHFVEIq9K0w/lcv4mLNyWoEde3J+rBtqW4fMnXs6uNNB6nKXtsDiRlJKnH3EDcFdScYSn+hcg1maLW2fblSMmUxdbbEE66srgwOrcPi9ud14YzDbO4WqjKoSWut2sbe4t+/BCcnPY6HVwlBCEPMpJR8G1XFGFK1pao5haObbBhvV74hptpmCDl4IawdV2Leye1zcjs63HcDiv7AhayW/lVUQwBFTWyB9uvOSN4v4lnsEtJFu7ZBZdNqkiQjUnsrMJVXZz18XWcgs7tw+N257Xh/D5Dm/RoYmgMER2CBfijVLWgQUxJVDI7kRsFkWMhspCrqvBjWwXRLtfcPSaE3VRuVN07WCMtTAmKykRpmHIwModEVRolkGZKw6qVjKxCIFt4TYLH6gIkAcmIRlglNIgZCB37dP3Chbm8UGBJVxZ7Ref2pXO789qw88i/oiQkuKIYXFUUF28I6CIaYYeALkIlr0YhL8r6SNaqIiI15mbrOKmprl8lbCHwNtcPTnH3piSYu3hWLsRIzBqTg5BkdPuSqRJz9UxRSFYkCyG5isiKLsVIn/0vCLLKkMW7Vaxwv2S1/dI1RlueMpRWEa7mriwujM7tg+J257VhN8M8H5Zas9GhunjqsThTFhZrM2I6eau7h68zEjgybsfRtdNGYWwi77lcvYI5eb1NlC3kNdJqKC6giQQjrJjAcuLKOOBpfL042cWIzRBglb0fFcHcZVG1fm0z26fMl2BJkv5MzQujc/uguN15bditXK4sZcw+j0XyhdAjcXUITdxtJLAGUxR5aBRFAB0KYanHFSJPyBtHkk7cvVZVzJct5u6dE3Xi7lHcOxrXrxDWFYSoJ0wKUUdlAU5SaeJnIlZkn524Q/CRVcHicrM+rX1tWRW213t2l+9C6Nw+OG53Xht2r8oopURFSQzm8ukiWjIkCroIpiYGc++qkoiQFiORjbyja1dUBgHyQFUQ5W+MywGi66qi9eLOQ96yvqYqpMbVqouXR6UQVuJZawgrHV3A7O5ssn0x4O6gElYANi9DSD6ZS6kIyNnqO1dYqZH39WlVGVm7y7cPdG4fFrc7rw07hjK8k4uSmCRJGvesKoUpcXOcxt9G8sqoKgJ1XyVvURVBG/IyVRqC/7xvUBab0CgKbRIhqm2b2g+7OiFLAiQ3r89ANterVRrm4pmyyFkIRkvy4EolKKreh9neh6ZfzYUsam37jXRlsSd0bh8UtzuvDeecxKghbRzLhzROkyFtfG1C3DiSsh4XGzI3SsMUhlYiV0I7WetnV0gbWqnAZgI3qkIBss1+VRLa5uJ5xrqExHzkU3D3TbJnp4uqAM9OCyQfPSWF1IoVEBmBNai5xWr9Qwz2hYmCaLBkSenfM6D0sqK9onP7ILjdeW042zDXoaNNkqSNw8XiBsoYb4tbiBstG52HaVKkVRJ5YSqCJllSlUUbl4taCaqiNZNt17gtAEAjPfAyn3ISI1Qhc3XfXF1krxeSLPZW5UtVY3HUUVDBw2ZFVRQCq09oY4OhzEfV2MxU1hJ3Np/v2m1AVxYXRef2wXG789qw+yRGPkRVY7D4W0mIlLjbEMbYWxSLuzXETUdUUtt6oyQWIzlzxIgb1N1AJ6aARq2uHbVd/RJHMssp5NWWrNnjXComFdy1szicEbEQVla2JNt1Tchd4nSuLjSBJLsGDUII6q6cc1JAUrBaUKwcy0ZiRVTVRqHlPBJ4A7IKxz0Wd3F0bh8UtzuvDTvHmCcTgU8y1+7arcXipmpCJ+tT1y4P2mSxtbp7tt6QNI7ELWQVJ7OEcR7YCXmlUQ6M5FWMvKi3OZnJRloVaUgNIDbzVom9ibmCIcvEtRx/7MW+bNjxQf2+VKrbavFMXetTuweP1W1Br/fcEzq3OSRud14brinGrJN4XBN/iy1pWxKPpM2TWJwphYnCKErC9xEawg65ElaikdXIqwSPw4nMyDuDqtS4W85GUFVTF5pkJLIE/Ee/xtZIOLl9W+wLRhpJluPoeebo4UF362xbxy+zgrQzkzVJqLOgCKtLcPlE5LOAr8V64S+Ap6jqHed4/T8CngvcA3gx8OWqqiLyGcDTgQ8EHqOqr9jzpZ91YZ3bB8Dty+L1oeHsuTJaopaSInf18mAlRHkQc/m82D4P1D/1uFvr4uUrkBvS5qMSd1N0oVVByKBVQUhQQsiEoPadidlILEr0pYg2j15bJ3BRFFkLiYVUyKtCSsHn+hbyIpg76EpDVwJObonuDmaMcNlGRAUnsUZqsgXRxsWDIkHCICZWNCBDIKhPOJ7DZGiwZbHX3b6swnG6sbNwicgAfAfwQap6h4h8G/ClmEHdFc8G/jXwm5hhfgLwEuD3gU8Dvnuf13waOrcPj9uXwetDxPmn/YTGNcH/7BdxMuKpZqIbN69kqlsl4fuIairClYREUxESlCBG3mFI9h0KmSHkSlYRJYZMKDG5cskNgXPj5pXtlAODSiXzKgRyNgKvVkbMHJzA+AxiGcrkXSSh/XHPqKkEh8biSmrTL74u5jLWuOIkOSKnjvor97HU04+5Diif+L1E5C3AewB/CiAi7w88C3gAcBfwhar6R5MXizwQeA9V/Q3f/mHgU4GXqOofetsNupUGndsHw+1L4vXB4XxPMGncEq2kbUg4IS0Tl2/q4rXbOhI3uJqIRtroy+CKYoiZ6KRdxEJec/tGMuuEtHNkFbKriJXYPLJlW5KiIZOy3W/OwZIgIiT1+03mBhb+GnmnhK3rft9WRqTTL3e08qO2H4uKWMteb4CqcJxvrLJQ1aWIfDHwe8CdwJ8AX+K7vwf4IlX9ExH5COA/Ax83O8WDgTc222/0tstD5/ZBcfsyeF0gIt8CfC5wP1W99zW8/v2A5wPvBbwS+BxVPRGRjwaeAXwo8CRVfeFZ59q9B7xD68xWdWKWUV2sKYomCTJu64S4OthfJe6QCdHduCFXFy+IcjSkkbyuIoqaiD7mtCXvJlVRyAuQNFR1kVUQsQL55JN852CF9KY07L6yZFRDnfd7fPhlcfU8IaPj/df1RlnUPvNZxsqTmlXKo3nOjsWd7I/A9xeRNqb7Par6PfODRGQBfDHw4cBrge8C/p2IPAP4KOAFjeK9suF9Nt3Udktzo9C5fTDc3jOvz4ufBZ6JCY5rwf8F/CdVfb6IPAf4Aix09+fAU4F/s+uJdqhjnsaFiFKTImWZPQZXFcPg9ZxtPG7hGepocbfq4g0Ki2zuXczERSJGUxFHQ6pKIobMUUwMIRNQFjExSCaIkbe0G3m3z8qWnbAZYeWkzRpYaWCZYm0/GSIpB1IOrFLgJGRyDqQkJIloCpZQIVgMMTHOzBWUsBKyJ1YCnr0etD4AInutZ8hYP2aL75Fncc8tbp8Cq7y3JMkdqvqoHY57JICq/hmAiPwE8DXAdwJvV9VHtgeLSAR+xzd/BiPpQ5pDHoIlEC8HndsHx+098/pcUNXfhPVwmog8AHgO8DBv+gpV/bXZMYJ5iJ/tTT+E5V6eraq3+zHbP7wZzvXTNH8M+/jU35miECZKohbRl+3AxMUrxJWoxGhuXgyZhauIKEbWK3FViXoUVhPCLkI6N3mXEidEDqKscqgKZQk18ZJVPEEdyFkpaW1N9kUUoJQLmYqwUqH2nnXSTyUjz1rGWoNslJaTz0KFkxtP4DcBHyQiD1DVNwP/DPhDVX2HiLxORD5DVV/gJP1QVX01bswLROSdIvKPgd/C3MbvutE3sQmd23AI3L4OvN7JGzwD34Ep4V8VkYcBL8Wqh1q8NyZOVr59oTDd7pMYjRv2IRR3pdQktjG58mFtc/s2xd0KaZ24Q0wsYmLhimIQVxWSnLyJIJmFq4qRvJmz43CBjDBIJquw1MAgTl4JrDSiKgRRUg4sRUl57IOcheRaQcsoJ6CMsLISIhoiW0it3n+JvYW2H+2v7evTHy0Fq3xjkySq+hci8o3AL4vIEng95qIBPBl4toh8HbDAYm2v3nCaL2Ysl3uJ/yEi/xIz0g8A/puIvEpVH3/97sbQuX1Y3L4OvN7VGzwNj8MESdl+DxG5j6q+szlmr2G6c1RltErC15v5Ztt6zjLLVi0pWvj+wUc8LTxD7XG34uLFmLnH0ZIhJqIo91qcmJsXEkNI3BJXLIK5eVfCiiEkIk5eSTUeFxjjci2Sk7bE35ZqqiIRWOXIcR7M7cuRq2FglSMn2dw+AdKQWCX7NQ9BScHLkWSMy2kytZARZFFcPe/CkijJ1h+SC8ldVYSm3vOUwSVwaYoZVX0O5tbN21+Hlb6d9fpXAB+yof2ngJ/axzWeG53bB8Pty+L1GQjAR6rqu9pGEXkp8PeAVwBfCNxXRAZXzRcK0+38BJMxNkR189pfxOryyejijW6gNirDazmDqYkwc/GGRkmUWNtRXDFI5h5xWdXDPeJJJW2UbOR10rbuXkRJzY9Z1lBJvNRI0sBSIwsrxGTRqI2iTlaiLHIi+C95DIU4gRyVrJ40yeVexVXEOFx1U/9M+q+4fWHW31uQEU7SwRH45kPn9kFx+0B5/fNYvf63A4jII1V1zaMTkZcBn455i58H/PS1vuH5059e46nNr+BYElNiTFO3j4bAeGJEohfXiyVDCnGjZ6OPQuIorNzFyyxCqkpiIbZeiBtQFmFFdHevLOfIGkhi7l5CCFnJIiw01YlTlhoJKKsYCDkSRDlJkWWInGDu4hCN6PZABi+yjzo+qb0h8Zrb60NVJ6VEVU1wZjVGi6Q3NpRxt0fn9kFw+7J47QOmPhu4p4i8Efg+VX068GXAs0TkNZjN/GXgizac4quB54vINwO/C3y/n/fRmDd4P+CJIvKNqvrBp13LDk8wkXHZ/M3nDijDMtcK7j1brYO7ez7qSQaPuQ2WoV4MFne71+KERTTi3ntxzEIyV+KKQRL3jscsKnmXRl6UhbREtoQKQBzrfTxuBkmFjCmLZYgsdSC5wjjOC5YaWeZIkMxKI8dp4Mizx4s8sAzRiuCDEiRaNltCTZ7oqgTlBF1BViWoWF1nxFy9bP2lqS3Ol7U+nvT/DKrC8vCUxc2Fzu2D4/Zl8lpVvwr4qg3tdwCfucPrXws8ZkP7bzOtRjoT53vmX6nv9P6cB/dp/iYunxTXx5Mis6GopWxongwpxD1yBVGUxCKsuEVWLMQy2UeFvOSqKoBJLK78CldVIYGomaCWlY6Y6gnZiv6XGgll4oE0cBSTJVaCsAjmDmoUQsjkMqWi39uYqS5KYdof8z5r+7OojEm/b4BCHSzQcQF0bh8UtzuvDeceYFLXZwSuWexJXE7HD624fkIdihp85FMpG6punsfdNhH3FlcTV8KSI1kRyBxJmqkKI21skqLJrz9pIIupiogdn8XG85MhBlteCas60QvAcR44ilYJcxwyMQe0DBAIimYjr87utcQepekbc5F1Qtw6dHVTf2+A6o2vyrjbonP7YLjdeW04f4y5juKhdvhGVTEjsbWvK4oywqkkQ4Ywxt0GSRuJu5AVt8iJ1YCKxd9uCSeVrAtZ1cstCZIykeCSAUhG5nxEkMxSBydx+TX3ybo97ZxFWIRk5Uh+fTkmFB+N1SgLnPzmEqsPez2tf2RK4jOqMUYIq9QJvFd0bh8Atzuv4TxVGb6shfdBml/Ndr2JxW2YN0B8KOrgI5+OhsSRF9gfRSsbukdcciWsuHc8nhD33vFqde3uFY4rcQOZW8KSyDQ5MlEVLoNKciQRWJAsJoew1IEo5uadqHXL1bzwL4u5eoNkQlpwMqzqZOFHw/iUX1XIg2WxVe2pDpqxaR5LHC6Nrt+0GF/GqoAz4svlvbrLtwd0bh8UtzuvDTsr5nZi6zKnQJ1bYParqUVBTNrHMfR1esMSd3N1MUiyBIjXdpZkSFEThbgLj8EtGMuJjkg1BhdZn+xlgc+6JUpAiKpVPQQNILCQMemwkDSqiiAMIdkoKg0MklhJYAhlSK3NSRCaeQLa+xff1uImt/1S+9H2l14+c4QUNpVjx8XRuX043O68Npx/onzHpuTIpuTJWD6DuXpixA1hnHM2YENRB0+OBLT+mi9CS1hbHkliQfL428qz16mSdtFkrKfXD0sC0UmMWhJlyQAKRzIqhOIypiDkLCwksSK6S5oZ1GbqqvcQMiK2XYesCthk57KeBNnQd5v6eSt0nLymYw/o3D4MbndeA9dUx1z+ZK3YfqzrZBw/77WdNf7mBfdDyCyiZatLCdGR13PeI55wJay4EpbcImV5Ul28BcXdSxyV0VFeSmREHi/XhUHFUjNJbBLxpRfbn0hioRGyKYulJpIKVz1pQoArYXzczVGJyUXLYtsTG4Tsw2/FkyVEtQnJE3XY6qSP2r5rFMeuyKkTeK/o3D4Ibnde7zq7XC20X+9n5az/AAAgAElEQVSw+Y/bmqrAlyVrDfWXuMwxW36pgydIyqinUUWsajJkVBNT4h5Jthp/GWv9C8p6AhYCQe3JECgsJXCk9lifI0mUyWijKEesyCKWYJHEUqJfnyuLbO/f3o/M7rfc/6RfTuk/a2smftlWVqSgPRZ3MXRuHxy3O68NpxvmTT9c0k4mbk0Tl6X5mxJYZ6SlmXO2uHiZSK7lQcV9C+Tq1lncbbWmJhYNaaMImz7aCCS1985AQkFzJfBSViSsFnQhK5La0NWoWuN9KwILyazKSKyQCTlQnjZR/rTEAaWNxU37RWfkniRH2i4PYt+82YfTlcUF0Ll9oNzuvIZrCWXMMVcPbZvH4NoPEmmeY4aX5Mz+6lBUVxRHJcZWiN3E26IoC5s5thJ30XzrgghZy2gpvwbfXuAJDSdwVFMzVgeqHEkiu7IY5ypYv96ilKL4XAn1S20xuBKL29pP1wIF7QS+vujcvvHc7rwGrsEwz4vu2/bJH36czyMw1nf6Exw8IRLddVo0meppYqRkr0+4JSw5wtvcxStq4hYRIkIoy1ZX+GVm7EuRUVMUnuVe+tyGi+LuBVgSIR+RsJFUi7CqiZJFSKzUJiAv159UWIllr7Pfq8XdPEGyoY/GvhPWivJ3/jA69oXO7QPhduf1xRXz3PWb7ms22h/V8mvsaiOg459nsotrF320k5UKjcNSo2gdhBRhjbhxU/ZXPV0iWNaa8v4Wd4uqJMlElfH9dXz/9vrK36b7kfrf2A9rV9O6eE1N6s7oyuK6o3P7ErjdeQ3A+aLsux69yQUsm9IMl19znZwoUmJz4xKY1HEauZq424y4YcO/0l6JLuKvZ3reSurpdUyubXbtdj8bwmjX4trt2s/++PkL/50DInKbiPyxiLzK/97nnK9/PxH5LRH5ExH5cRE58vaPFpFXishKRD79XBe1D3RuHw6398Xrc3L7kHBtc2VM2mfr821YS46UXS1pKyFK4sGVQ5sUKe0LRjevxN1a4g5N3jpKIKk/zJIIkqq6WLjUSWoJloVmm2zc3y8hnhwpX5rps9daEreKaZokkd36aNf+LlBBLk9ZPNknvL8W7O2BlXtD5/bhcPtyeX0wuPa6lJlymJQbnaIqgErgTRifCLy5kL4d9RQZb6BVDXaeQJT19flx8/KjbY/uqUNht1zX1vvaoirqSKkN+3ZG2tPfHiAiDxCRF4nIb/vfP9lwjGAPrCyPb/8h4FMBVPV2VX0N09Lcy0Hn9m73db24vS9e74nbl4FLe0746CKdQuQmBmfb1xCLvQa071dicduwy31cFyjI5blqPygiCXgR8M2qqlzCAysPFZ3bF8Dl8vpgcGMN844f8CaSzttaqR9mrlFREGvnkEDW8Wc0iJB0TJLM32+5w3VtxA0i8ikC57w4z5OEn6yqbxKR+2CG+XOAH+YSHlh5UOjc3hv2yOubFjfWMO9YBpM2fH/nbT6xla2rTr7ySfNGAqeZOig1oOV8572GrbgR5T7KPpMbOz9JWFXf5Mt3isiPYk9s+GHs+39DH1h5UOjc3g/2y+ubFpc29rFMVHLahCUJm/S7kGZn8lwQ7fvZ+2/vpl3u43pB0n7+dn4/kUFE7u/rC+BfAL/vu8sDK8uxjwRQ1cer6iNV9Wke8ngZ9sBKuOADKw8VndsXw754fR5uHxqu3TAroxOq2ETaW/atvfQ0wvpjcvKWBzK2JEmMaiCTSapkb0maq4po1+fHZaY5gm0kLNdz2oMiN95X2w9NX4jq1n27QBQkyV7+zoErwEv9oZSvAt4EfK/v+zLgUSLyGhH5H2x+WCXYAyu/UkT+FIs51wdW+gMwPwP4bhH5g/P1yB7Rub3bfV0Hbu+T1zdzdcfuoQzd0rs6W59vA6jNUFX+yq6s4n/B/4RMeW6ZzZJly8CSwduVJYGlWllRclfPCuKzlQtJqpnpNu7WEjf5qCjbhqXatIllPoGlDlXRJJXabtc3u2aV8Vbb+6yzu+zQR7v2d4MbHYtT1TuBf7Rl3w1/YOXe0Ll9UNzuMebzxph37bBTVIUqZG8byTsSIvmfEcWXEoBkj2VHfM5ZIfikLagNRTUvrRB4/WJb4mY14iYgKU5O8SdAuAvH9Dom1za7drufDZy7FuWwSz/rze2qHRw6tw+D253XwB6Sf6KK1AkE5vua5ta7Kb/E5UNHxj9XFjUG5g+XTGKPzLGhpUpSIYrWzzmVN2gJPMOEuGh1FzNMlEMi0Mbgirpor6/8bbofrf+N/bAGdbdtB2W8DV1ZXF90bl8Otzuvr8EwS/l1rB0vtX3yh//CZoHsk2pne2KCqj2uZiVK0sAqB5sPVjJLjfZInKwsQ7Tx/GTIRyzwR+IoLCWDGmkXfkE25LTMwFWK+ccZuIDq4iWKi+cTixPsmWhEruYFV/NRvZalRpZ5aNbtb5VDvf6Ug5HY71Oz+DfDvtiyoY/GPtWmT3f8ILQTeN/o3D4AbndeA/sol9vk2lWC+9N1S5zK25IKg7etu3zij2G3h0gG1TqePzvBoj+dYSkBfCLwWodZ43KG1BA3+3ZREy1xi6u3VHvwe0I48YdXLnUwpbPBzSsxOPX7GqcfA1HxL7uc3k/XAKG7fNcdnds3nNud14bTDfNGN0XdxaP5kBh/ESeqo4lLzZIkWY24KZdYl7B092qpkYWmGhvL4gkLhODJkhNJPgG4EXiJzciVPWkSN1x6IW3S0cUrxD1hmhQxwo6kTmrXZddXEiXBr188Brc5ObKpXyYKonX9Nrh/mjd8EF1ZXAyd24fJ7c5rYBfFnLN1qJN2jrl7Ig2htflwyNRY3JTE9nTeVQ4MEljlaI+6cdcqolVZLHVwPxJ7jpnQENimN8w+aUv5bNefi1ZqOOfEjZxoZEl0EktVFOVasgqrHMl+vStXGe396Ox+y/23/XJa/1mb1j4nnzJ/QVcWF0Pn9kFy+7J4LSK3AQ8EyiCpT1DVvznH698PeD7wXsArgc9R1RMR+UrgacAKeDPwr1T19aed6/yhjPrrqDamvfxClkxDHtetJhF7PHyJTyV7PPkqBCTZTFXLFOtMVsd5vKTjvBgfsy7ij9mJ9vyybM8xW8qqPp0hqj/94ZTx/3PXriiIE43cpVeMwDpwV77CUgeu6sBxXnCcB96VjjjOAyc5cpIHlimyzIFl8lhcCmgqMTiBJFbortN+qett3zWxuJ0/h64s9ovO7cvn9uXz+nrMnPi7wKNU9S4R+WLg2zijtPR8hrlRFa3LsrZOoyo8FleSJfbrCzkHNGRTFoj9Sov9Ui+cWCVREoOddKmjE7eQCK4QytMZbCJwJSMbJ15pS4ayBk6w97Ba0liJu2yXTWIkMyqgVTaXr6qjHFwIyJgU8fufq4ptfbepn09D6Ip5f+jcPhhuHxqvReQBwHOAh3nTV6jqr82OEWzmxM/2ph8Cng48W1Vf1hz6m8BTznrPnQ2zZK0/eMXtq+5f8yFMPzSdkLj8clr2GlIOpKyssj3GZuVZ4cHnhl3mSAhqv6ABTnS83KUmEGoNJgGiChl7UgM6nZil1m+WwnoCV/OiJkOWjWt34sureVEz1avcZKyr6xf8HoJ/IZtESEva+iWe99GW/vT+PhWXryzuNujcPiBuXz6vr/fMiV8AvOSsi9jNMFeSWixOM15xLkj2jq7rjGPV/fk4tl2eQBzIi8BqNZ7+ZIh+euFqGIwYMRC8xOhKWFVFUR77Xmo9yxOGl8RK1oWMJ4/oZB6CpX8BEsLVfFSJm1Sqi7fUyDvTLSw1cpwHjvPAnemI42Su3vHKXb4UOVlFVimwWkVTFqvgbp67em1/lG139STr6PZldRdZJ/29DQKEbpgvjs7tg+L2deD1wcycKCJPAR4FfMxZF33+GHMuv4CM8bjZL2bryqji8TojOBmLx4m5SDkoKQeW2Lyvq+wxOf/1DjqqiqslLgdcVeUII7HF5Y78CcS5krWd3rBO3uKjnUom+sRLiNpkSFETx+7qFbVjrp4tlymOWescmvrOsaxIGrdvW/8UpVXdvbOUcoH25N/e0bl9+dzeP68PYuZEEXkc8O+Bj1HV47Ou5drmyijBfOZu3pTIksXcxJI+dhJqFnIQxIvyUw6IGIlPsqmHIMpx8svz1y88MZKC1KRJECWz8pFT/vwyv9b2iQxlgpZ2BFRx7bJKTYaUuFsh7kkeOE4DKw2cuKpIOdhoqhzIfg9Znbza3KsWFTEl8cai+0LkTf29Bb2saE/o3D4obl8Gr0VkAO6rqnc0Myf+gu8uMyd+ux/7SFV9lao+fnaOl2EzJz6fZuZEEflw4LuBJ+xa5bGbYS5lLVnXkiBbfzVzib/5/lwUBWMW25XFKhl5l07glSgnKXIUAjiBs1giJAcjC8GIGFXJIiSxWF7QTPQLbJ/OUKY3TCoTVVFqSI2wY9ztOA+sNHKcBpYaOEmmJlIOLD32tkphXVFk8S/tLO6WZ/2j0/4r/SnF9Wv7fRO6Yt4POrcPi9uXx+syc+ICmw77F5jOnPgsn1VxAH6ZzbMnfjXwfBH5ZqwS4/u9/duBewMv8HDIn6vqJ592MWcb5jYm1PxZDEmb2JsSktgcK2GMPwURsigWGhMQ0JUAgaQWezsJ2QvyBQEWObEMpi6OYuI4DyxCImtgCImFJK6EhZUReTJlEVb1oZJlOUdJjpTlMo+jnibuXY7cmY6qi3eSIncur3CSI8sUubocWCWPvy0jOQmaArry+NtKkJXf/2o9HheS9VftO+/LeR9P+n8GoSvmC6Nz++C4fVm8vs4zJz7uvNdzDXXMJW6krjLKdok9qf+qurKotY1io3wSkMQGD4mRPedAmTo1DYmQAyfAItt0iEdxRVZhEHvS78rHPi0lVgLnLAS0PoK9uKPzBEnJdOdS6+nEzSq8Kx1Ze44TF2+lYeLmrUoMzucQKDWsJKlDsIqKohJ0JOqYuda1fty1VM6UxY7HduyGzu3L53bnNXCeqow6WoeavZ64K+7SyMy1qXGokqX1ek8yaApkyaRkKgNglYyYWYVliOTgpUAhEdKClQYG/0kdQmJFsMJ7SQTRqjI2PfG3kLbMC1BIm7DEh7l4gWWOnHgZUSHtMsUJcVMKpCRVUZTESPniFgKv9UXTP5P+K1nrWpCvZxL50Oo9b0p0bh8ctzuvz6OY3R0x4jYuXwKJxYVxV0+UsHLS4ZNmDZBx5RHFi/JBNZAkkrNafA2IITLEhAKLkDkOmUEyJ8OKQRJDyByFRJDMQrKRNySCu3mnPdW3TgTu6sFIbHManPisWiu1sqE27lZcvJQDy2UcR0ItA6w8/rYSZCnVpQvL4vLZnySdunzFBSyuXlUXenYGu3wJbiBE5J7AC4D3x/TTz6rq15zzHNuGrX408AzgQ4EnqeoL93rxp6Fz+3C4fQm8PkTsZJhVdXSYdEySqLt946/k6Oa1bg7iH5KUdY/HASRTF5BJBCu6dyyDklWIOZBjQpKyksCgFrcbQmaFDVVdaXDybo7BFbRPkyjkLaOzTvJgo55UasytZKinasKJm6yuk9w8zia3pGRjf1RXb+Iuj+1tv2+DcGku339U1ZeJyBHwiyLyz1X1zIL5BtuGrf458FTg3+z9ik9B5/ZhcfsSeX1QOFeMuY6QKkmS6p6MMaaSsW0/MA3F9VHqdLLJkiEqWAzL64ZSsHVVCBLRKDa8tVxwyBYDi8KQc1UR10LeQtSsodZvViIXsnqGumSpU5JKXG2IayRuCStNv0xdvWnxPaOScLKeOeoPLkVZqOpd2MNUcZX7SvxxUHsYtnq7H3Mpeqlz+0C43RUzsOvscm1WNSkE+0XUZMuwoqoFEVMPYeXj+gHUCByyoIPaKKkIGjGXj4AmNaKqkKNNAJNzIIRMCEbKoyERQ0ZEWQQjbSxL/zTLhDFlvd6GSl2W9RKPK8X0S094pBw4WbniyF72tIx2bcXFK8Q9cTfPR0OFJb5t6yGBrKw/Sjbb3D21tmx9Wpaezh9dvlNK5sJq667rDhG5L/BEbLgq7HfY6o1B5/ZBcvsyeX0oOPcAE1FXFqXu00m3piyKm+NuXnH7SNL8aroTGZzIig1r1YwESBLIQYy8Pk4zBkHE3lNECTkgmNoQ0Ql5N6FOAO7zAai3qQrLNJK3qIjy1IacSoa6ibulKXHNzZOpkpi5fRNF4f3X1neOczTsoCz25/KdZ9hqKcb/MeA7vUQI9jRs9VLQuX043N4vr29a7FiV4b9srZunHp/T9kNqiJtAREEaIjumBeSCz9kyxuU0QFSSUBMyJYOdXGmoOnnF3ifpOOtW+fZvUhXabNdH5iiVzGUmrdUqGml1JG1x71oX76x5A6xdZ23NdonJ6brbxylTPIrH8faEnYetOr4H+BNVfUbTduFhq5eCzu2D4vaeeX3T4hxVGbnp1NLJQNv5VU2MJCYpIjIhbEg2U1aBJChj8NXjdPZWAQ2KBs+YK4SgiARyNBUhokSxxEkhM/gXZ4b6gEwna3lkTl1P40xa81FPuvJYm45KglY1JLuvdSLPi+0Ladf7cNK/p436q/2486e3N/iopvfEJv5ucaFhq5eKzu2D4nYvl9vBMGtVER4XShlCQCQTRNCUQYKVDlUXREb3LWOv9dibTxdgZUURNNocBBpAg1g8LihE0JXF/FSAoOQhg6sIiZZNF0/5loy3yGbijvcjlSPZ6zPL0xk0lUfnYDNp1fkBiorwL+hyVEpWNiRjfK3E3hKEpXopkZcYeTwuJEVW6uuZsMrIKnscbhr3VC2dOL8R7PgbCBF5CDYRyx8Br/SwxTNV9fu44LBVEXk08FPA/YAnisg3quoHX8/76dzm8Lh9Cbw+RJx/ovxS79kkTCSMNYwiRuKcQMRKh0IC9V/BgMfc6jkFCc2vqwA+kUtJrJSJXlQDCKYyPAtu5IUcxrKnCXmlSJVyC6PbVyZmqYStM2hRFUQtpl81menWdVvJpK3WvLb1nO7yhaIwqtrQSUJkzd07A+EGE1hV3wgb48TsYdjqb+MVHpeCzu2D4faN5vUhYsdJjIysAo2rVzrc9kuQMfjvsbeQLHsNeEC/rDfnVhtzr0FNdSBIMIIa87E/EfvCeGxPg1pipcTeglaTcZaqKO+ruaz7l6UQV11FKA15bUlL0CxTIrdxuEnB/frcAeYe61jfOSOuFhW3BeX1HRdE5/ZBcbvz2rB78i8rpISkjEpyxZCRlIFAIKPL8usMMZg7Jc2oJ41q7p66qxfw+QSgzC0gCyduEHQF6srEyO3rAsSRrCo+/SLj9vZ7aQRfxuZBgKok6lBSJ2ctkE9UIoclNftc1mssbqmNu+fELSVESzV1sVJ38RRJ5upJUlhZ/5KSu3ynxOIUZNUJfGF0bh8WtzuvgZ3qmL2TNGPBsZmqSGpkKb+k5Zc92cQrVu+p9tBKwMSCWuyrlBKJu3Xq7mAwNZHVzleK+DUb90R0VAQWjPMT+7bP9LUGnS3b6QuhJmrqIIJmGba4e6PCcJdu1RC5Ie6oNOzPnjPv6xtURSXuVvWgbjg6rhmd2wfI7c5ruMYYs8XfMiKC5mwfuoRaKmPunrlrlcBBx/kEgp/La0SNsK4wsA9fAwRtFIUAER/6amVIRWGU0qUqGDZGQcs9+CGFtNooidrmhNRGWaSWyI2SWI3EncbjZsQtWWuPv1lfZciuKFJeJ/CpnwX0es89onP7MLjdeQ2cI8ZcsqilzlCDd3J2MpYPRAR8xJQ9rX1UFUEbly2AqjphTUVIcdnKfiexODE1Y4pE7FgjLkZcRkFBu2wxUxVS1gtxcyGwVPKOgwum5C0Kw0Y4taQdidwSt4yEmsfi2hm3pOnnWilwCrqy2AM6tw+O253X55nEKOtYUiSWMFBAQvDOD0gQQlI0W1ZDEjV+VpRCjni7VCURBiOkBghxbK+uoPgyYnE5GONu0hxTttvl5Eamy1ZFgJPT2831YxaXmykIbVRFLiQdEyIl7lb3nXjMLastVxlW2eJvq2QByVJS5EmpbRB3tTsuhs7tw+J257VhhyeY5HGpM5eklBflbMmFVYYhwCojwWNj0bPVAVBTFuCZ6lLfqaO7p1qILK4iGvI2Lh5hbIeGvLCZuPV+msMaMk9VxaggpuTVSuxJtroprq8qI48xt0rghrgWh8vTMqL6l6f9vu0+Vl1ZXAid24fH7c5r4DzP/IuuHpyoiFgWO/t0LpKtJMcJHJINLa1hpeLCeZ1nJW5QJ+/o0mkARD2jTZO9ZkJeaIh9reTVhsQbyaugPhpJG8JqIeioKDa6eG3craiJrEiJv2XP/ntsczIK7dQRUoqkdMr+jp3QuX1g3O68hnM8wUSTVW0aWR1RYSk2WqqwMliiRLM6OQWiICkYOaMQBvFZuBQEcvTaTRErOXIyaiyP6dEpSQuRYaoqGtIWt7BFnUQFpqRt3T2dkVqL+6d124g8knVMqKiXCdmxwdVDS1zbzsgyQXXxkhE3JVdoyUadnZYoUfoIqX2gc/uwuN15DVzLM/+yf5oewJfsJUUrr/+M7rt5osTUCFYL6q5dxnYXd65kN9TFiZakR5qSVoNMCAynuXwbPtzTXL12vSiGlryeyCjr5flmYaWNCmkTIUZkPENtMbs8ZqpXiZoUaYapbkyKbCKxYrG7jv2hc/vyud15DewyV4a7JnUyWs11xizAPgycWO76qR8vwXke3T8TRWJhqY4qI1Pjam1SZOLelRFRbHDxKnllsr35hvyQWlNZtpkQuSVrdfN0qiBo421VcbhyUIy07tKRsfhkcZlzduLmqixMSZQEibXrKXXMu0x01LEdndscILc7r+E8VRmqaEo2rj5G+8xVjSfBstZkRWOwGbecwDa3rKBDqF8AGWxegFqfGaWJtXmbu4oTgs5jcqWdQtxWNmy6kenmGoFnsbeW2OXXv66X4ypxi8oYk0ji67XQviiJdgRUzuhq5evJSWzHnVaVYefrM4pfFJ3bB8btzmvgvNN+hjC6eyQg2ocAlMJLUaXUXtq2k9bbEc9eO3El2PSDpZ5TRFxhSB1RNSZFjKBFZZS3tWPmzNxwD3MutIoCRhUBlagjsbUqiVFt6ITULWEnpPWYmjh527aRxGVZ3L4zVINiX4aOi6Nz+3C43XkN7Foul638R3CFkZy4pKb43dkSjayyGl1AI3F0wo4EL6+RKKPCCH5MqzrKucO6S7eTi7f13vyla66fkxQmRK1KoiFrSVTU0qCk43pD1uraFeIW5ZDySNyUqpKoE71snS9jB+PdcTo6tw+Q253XcK4nmMRKRA3BCOxKQhMQxEcbOTGjIlkqSdXbbfSU+WriZMbnBijbEhhfV0hZSF2uqaqJmbo4B1o1Yct2n9ZkxUheZfKASZiSFSbDT6VZnyiJbKVYlZyuLLSJ0Y39vgVKd/n2gc7tw+J25zWwa4w5K0JCJVSCVncuK+XRChp9wL8EI/NEOZRZxAV/Zo693tuKa9e+ph5b0LbTKIoWYb1pDRs4sfY8sjYGlnXa7uuyqb1tq65wIaaTtbh5fu4Sd7PtXBMl2xN/fs5e73lhdG4fGLc7r4FzT2JkhfUkLx8KwTveVUI5zsk7umLBMrUTUoZxvVULQdZISkNSbcgsm8g7O379HrYTor68Jc0mUs/3N8Qr62vu2mxby3wAuVEa9T12cOVUYdmVxd7QuX0Y3O68BnaKMSuQfY4AEBJIML4WEpcPrCRRcKJOFEFoEh4NqaEeJxvIWknu7RNankbS82JO6omqyGvtlWh5Rri2vSVhbjLRsxFQOnuNtufc9mVDLePdce3o3D5AbndewzlG/hUfyR6jkyrxFKycqIWEaZJ4tn9NDYTxXLMDT7+usItvtyPOrILYYiBnr1srA5q7bDPFMHHp1p5/dloog9Fl7Lh2dG4fFrc7r4HzhDLWfnXHzrNRTi3RTu9Ylf2Rbu2Lcw04NZZ7/pOd49hrf19VJXeXbz/o3N71ZOc49tret/PaIKdOLSnyZuD1N+5yOk7Bw1X1AWVDRH4OuP+ezn2Hqj5hT+e6KdC5fVCo3N4zr+Em5faphrmjo6Oj48Zjj4Gsjo6Ojo59oBvmjo6OjgPDTWmYReQDROR3ReSdIpJF5Ov3cM5bRURF5PxToXZ0dHTsETerEfoq4DZV/fDLvpCOjo6OfeOmVMzAw4E/uOyL6Oi40djk0Z3Xy+te4eHjpjPMIvJLwMcCzxSRvxORHxWRb/Z9jxWRN4rI/y4ifyMifykin9+89pM8BPIOEXmDiDz9km6jo2MCEXmQiLxIRN4sIq8TkS/z9qeLyAtF5EdE5B3AU7e0XRGRZ4jIX/jfM0Tkip+jfC++WkT+CvjBS7zVjh1w0xlmVf044FeAL1XVewMns0PeF3hP4MHAFwDPEpH7+b47gc8F7gt8EvDFIvKpN+TCOzq2QEQC8LPAqzHefjzwFSLyeD/kU4AXYrz9L1va/j3wj4FHAh8GPAb4uuZt3hd4L8zb/NfX8XY69oCbzjDvgCXwTaq6VNUXA38HfACAqt6mqr+nqllVXwP8GPAxl3itHR0AjwYeoKrfpKonqvpa4HuBJ/n+31DV/+q8fdeWtidjvP8bVX0z8I3A5zTvkYFvUNXj5hwdB4q7Y6zpLarajum8C7g3gIh8BPCtwIcAR8AV4AU3/Ao7OqZ4OPAgEXl70xYxz/D1wBs2vGbe9iCmIxlf720Fb1bVq3u41o4bgLujYj4NPwr8DPBQVX1P4DlwTc+H6OjYJ94AvE5V79v83UdVP9H3b56GbYq/wAx8wcO8bdvxHQeMdzfDfB/grap6VUQeA3z2ZV9QRwfwcuAdnpy7h4hEEfkQEXn0Oc7xY8DXicgDROT+wP8B/Mh1udqO6453N8P8vwHfJCLvxIj7E5d8PR0dqGoCnogl7l4H3AF8H5bE3hXfDLwCeA3we8Arva3jJkSfxKijo6PjwPDuppg7Ojo6Dh7dMHd0dHQcGLph7ujo6DgwdMPc0dHRcWA4daiFoR0AACAASURBVIDJkVzRW7jXjbqWa8M+qpBvgvznO3nbHe2jpR7/sffSt7x1Pw+t/J3XHL/0Znz8zkVw//eKeutDFzsdq5dEkH2862UV6cs53vl3XnNcub1PXvu5b0pun2qYb+FefIR8vG3Mn+rbPHRy05OEJ9jxScJr6E8Srqu/oC+cPJ/ujrcmfuulDznlgnfH4oF/ts9nrN0UuPWhC17+0ocCkGb9nhuTmJnuS7PPd23/zJzmLbyZHzfHOR57eibO+pbELUY0zL5/8+PC7MxRtu8Ps9fGxkbEB/5p5fY+eQ03L7d3G5ItUo2tBF93Yysio4Fs1q29+TCa11SDWz6c9lztfpga37Ufhz3qgbUnJTfbrfH19mqAi2H1L7e07e0XPmvzmlzPIzlPz6W5MdZ56w+Coiy1P+b9okiaqyHOZJJqNbYJrYY1odVYJtXJs7IzkHRct+PFjxVvl0k7QNZxPc0M13z7IoizH4F2O4iutQdfRsmzdj9OpsY+4kZZrb0Y8CAyrhOIIqOx1jwxzgWd14azDbMb5aqKY/RmN8KtMY7FQLsRbg1wkKnhbY35xBC3r1s31Dox9lvIe5rB3qZ8m3aZqFhdX2Ydvzaq477G4KofV42zqp1XM8SAplyvVVQhZzSA0dqIqTmwzTgryrH2x7xfBMUoF0NcDEIxyK0xPtFivM0ItwY4qUwMbzG4CZkZ4jAabDdQWUfj1B6bdLPOzafo37BFZxcDC3OjnEEhMhrg0B6rOhpladZVie12ea2/7shvI6iSEVfemaywqLcYNhrnzmvD+SYxkrCmkCXGUQnHOKrrYohlatApSlqmBlhbNT0x4BsMsS91kwHeJbqxgb/SGF715Xi8Tg1zVbsb2r1NAFICop2zNdTaGHYRNCUI7uzljEpYD2tsuIWrXVnsDakY4kYhLxtlvHRj3BriYkiXbkTNKAc3yKPhHdVzqIbVzjOStbY355ojbzHWLVrDWhBnCrg14FFyY2DzRCWHup6Jatvl2EVrzDEjPRpnfz+wILdaeCiIuFHfLpw6rw07GWYJAjGaoY3RjEgxmsMwGuIQRhUdw2hkRdB2O46GXUXMN2q3A+PrymcYZGqIy2o10ue/+erFVePa7tMaphClMby+rxjwpHUbgJTrsdKuq477ckZzroZaVjKGOsTSJqqCkNjGUUVZ3gxZywNHJrPU5IY4uyI2g3ysoyFeEshqxnipsRrfrIETjWSC77OvVNL5ttgxGvx1m4132YbROO9ikOcI8zDETDW3RtfUsBllM85ufGVlBtePLduBzJEkN9RmrBeS/LzKQnNV1FeK0XajjFgIZyEQiGvX3Xlt2DHGXGLBYTTKRSlLsBBGaSvrIZgKLsZ4iFNDLI2hjY0RLgZYGNvKcaFRyb6Yb58L/vnLzDCL6ihYVRvDXPY1bcEMdDXOIYyGW9JomLMZXVQhec5a1AxvjCCKpISG4Mck7/fNlllVuHotv0YdE5RkXlatRnmJhSvMsIa6LAZ5qUNVvJnAVV1MDHExwgAnOlQjXAzwZOlGNzcKOjcK267t/J9ziR2PKnncnqy7Cp4uzVgfSfTXmqFdSKyGOrMkaK6KOyFmoElufBUks1Qlm9IAEbIr5qTKsOG2Oq8NO8SYp/HhiVEO0YxyCWEMsRpiLetuqDU26nkIU0McZGqE3QCP4Q3qscDEWNfPcP5Zbvps5z/E1RCXbR3Xq/H1Y7KOBrsx1pJLiMK3xZN5qtZXrpQ1KLKiErS+X30LC3tIsny9iFjMWQKbZLMCy07gC6HEl0uFRGuUl9UYmxG9qkM1rld10RjhYrADSYWrejQxxMWYFyNcDPBSzeiNx8pku6zbMsyue/1zDzNyF0VcDbQb29LWbi8kVYPdGutFo4qDKLeImBLWTJJAlMyCVVXlSbInO1djZYu4l+jGeSHi/W3x/blm7rw27B5jrok+qUaZMIY2rC2MIYuikOPUGCNCHsJUEcdRJWuYGuyJ4RWa45p2NsSbdzDMa0o5S1XFqEwMdzHAZV2Lek7Wpq6oQzHGrhhaA11i15IyxDFUYmGQCCR73xLHT9tjbRaLW3cFO86PkuhLjEY5qRnmk2pghxqyKAq5GOSreVGN61VdTBTxUmM1yMsc63FLjRPDm0vsWhsjPavsGK93PbQxjy2XcEQx2EFcDYvFhFvDbQbYjPUipGqgR8Nsx2eRetwtYTkx0ABBLcQRUZInFgOjUQ7iZYK6XlpX743Oa9g1xuwGVWKsMeWqlBfDGLYo6yLoIlqsOAQ0Woy5GOEcQ1XBGl0ZF5UcqQbYFOOokAsfVZp2OF9YYxK+mBpfM8q+PVHGti0lnJFKDNqVdBqP05Q97gySMpLU48kgHubQnGEpbuxzTRSqdbYZ/pRMNW+5jewqruNiMINsib4SUy5K+U5d1LDFnflKXb8rX6mGeakDV3WoRvhqXlTVu8rWttJIUmGVYzXAZT2r/a2yK+imDUbVrDNjvQnFCEujjMsyuEEGGEKatA0hVYNd1gdJLCQxhEYxh2U11rfkFQtZVbV9z3DsKntFCkJUZSlWXbGwL4qFMdxQzxV+Qee14cweEFfCVRW3RtnbtY0rRz9+CG547XU6uAIOQh5kNLSCb7taDlND3BphC3mMbbBhvV74hpuZJPekto2hjMZAZ7VjtCjmcb/EcTus7MelGHIrLSxqOKCiZsiD7decEbzfxKszSrkhHrUI7hKoIkE2Fmhok1jquDZk8hhbZmqUT5iq5KUOLDVy4sY3UQxz5DgvbD1HjvPgRtkM73EaqlpeuRrOGwxzymHNKOsGAw3rETmY0r01yLLBOMcQ1wxz2R5yUc0DV+KKkN1Yky1pF8xgE0o5YCBqrglBwJOpRtoTSeMFi4WNIh7TlwyzYEbnteH8PdAm9Jr4MENEh2DJqyhVCWsQU8nVULuRbtRxjsVIC7kqZj+2Vcftci2UwYSdm0rpaugC1gwyTI0vKhMVbarYDHVIVBVdgsSmos0VM0MsVpaUMKuc1S46CUhGNMIqoUHsx0/HPl2/cFn7NlqSpBP4WrFpmHXyqotJgs5DF1d14YZ5qEa6/N2VjqpCflc+miji4zRUg3ySh2p4lznW9WKAs5qHlPI0xmz7m2vfEIMVGQ+wr+RUNcdgWrsYYCmGWdQMsa8fhRUnOTJIthpklEFDDY0kAks3plkCWQILSVxVJbOqhvoWWZLI1o+iBISg0nzpSp/PR852XsM5Rv4VlSzB1fLgirmEL4aALqIZ4yGgi1ANs0YhL8r6aIirQo5jbNnWcYPNely5GOdtYQ04JZQxNXDz8IWVwjEa3RpvhpBkDGkkU9wWxjC1LFmRLITkCjkruhQz6Nn/giCrjKWpPURUyuZKchXQaOOl9JRh4opw0gl8YSSUE1WWCktMLV9141vCF1d14eumju/KR1UhLzVyZ7pSFXExxKscWGnkJHmMWYVlih5HFlYpesm7GebkS1Uh59Egq04N8bbxUVAp5Ovq+skNdFBvs8EhUo0zDNFjzKIs6nrmKCYGD2cMknlXXDCIKeqlRhaSqoK+p0aueKgDIEmwKo6gXva3AoGlZc836g/ovC7YrQfmQ65rpUWo4Qv1OLOpZosjm9F1w1xDGfg6o3GOjNtxDFtoo543GeZzhTEK5obZ28QSxRsMsxlkDSW8YQLYjLEY+d0ol/yLtfnrxQ25mNFmCLDK3o9qmiR6iV30UrlaHri9ftVicbtNwtOxHblZlqRcwmqT58q4hCyu+nKVA8d5wXEaWKqFL66mgZXvSzlwkiMpW8XGMkUzwllY5VANsZW2h2qENXt7uUg/xtZPuZmSYhEquS0qqEg1zJBCrkZbRElZiL4/qVdehIyqsAqBQTNDMbghkdWUdA5Lq0Zx1UwGAiw11tGEJxoJkr1NiybZOhdI57Vht3K5spTpoJFxQIiHL4bGODeGWGtow9RyHhq1HECHYoypxxUjPTHMcTTAk1BGq5jnyxbz0IUb4UkogxK6YBpzDlTjK6KeDCxGeFTN4AZYmtiwiA0gyW6Uh+CjBoPFnGd9WvvaMoZsqmXOCFdzJ/BFUGuYGYdPlzBGGTBy4uGL1igf58ENs8WUj7Op5JMUuZoWpBxYqRnm42qMA6tkBjrnQEqNIXbDXAywFuvlIbXRKDftc1TelwSgrxdjHbS2hSCNYYYcMyFYxYaqWNgjm+CKmkk5s3KRcISQZSy1SySyCAtZmBebrfwuYGV3udRvSx77tylPnNcyd14bdq/KKGVyRSUPFs7QRbREXxR0EUwpDxa6qCo5QlqMRtoM8xi2KAqaAHmgquPyN8acAdF1xdxGKM5jmMv6mmKWGjOu4Ys8quCwEq/IgLDSMbyRPVSTbF8MeKhDCSsAm3MkJJ+oqFS7ZKOqrLAyOu/r06oyVOlJkj0gqdpgElfLJXRxVRfcla+4Oh64Kx9Vo3xnulJV8kkeuGt1ZBUYOXK8Gqo6TtmWOQdyFlYrI21WQZMp42p4UyE5NdQ1cnNmkE81zH5IMcpVyJRtJcXRSEtQJFp8+f9v7+1jbmvTurDfdd1r73NemJkgDrU6A41pMQKJsTLFEFukQlAGCiU1hTigtECiDEFCKHSiKTRADNIKJB3BUaMQSwb6MYEm1iEQKMGCMoI1EZKWSFHAdhxbGz7ec5697uvqH9fHfa21136e/ZzzvOe8JOtKnmevtfbXWve+92//7t/1BVJMk4BZwSy4aR2NTd5opDhNDYdu0sasjKfccCDBo2bRFydtGbkhZPp8OAQByzjsOmckzNZ3dJ/XZldKGQ4gwZIXDsAiPSQLXoKytKW2PICZBmNm5H0JzMGYWQswY8miCzM4Y81bVtiyFiefaj2mRlrVwTace1KeLwDEfuErizb5wlizCHlYkP0YUYdlCqqPodj7oIyrySOxErl8Ibov+R7EOlzGKIw5svoifbrKGSdtC1C+6aYjn1yyuOltwY7nuUHEdGPpzUqlCEE7DxCutz7nciVXnNBAud2yOv+JBmMuMmCu3gg2FwlGDNjm4gw4MNuLiQwWzaTupOyYPCsQPAN9wpFnq4OhJlswyRhDKqsRUEbBbNk+r83uWcSoAHIboXHalo6+qh0vQLkNwM3HtQLUhUUbe9bFpBrhc4qMQU/NudJgbINzmdwKAGKV3WKpaPKFfxFC7vWsPnZpgsQjL4IxAx55QUD3zEAKwNb0bANkrKX5ctVXH6Q2FqRsjsAY3zvMwor2QPyHsgj96spn2nI4+mYxHXkByjKZE9Ade6fO7tgj9M6YZ4a6niyzs2QhYC4MWW1llUDsqEUVqIH7MWZ2OS2OtwHUg+DYMVV2YFYADSoKYdPhtFmY3NSAxiP55UYKdPCMWUxXPqHhxM3SteHask44oKPTsibIlu3z2uxuYM606OIArBpzC4mDhpbcLoBys0gLmZYOv8qS5eCTqjgCkzVXzbmNpZCSLqvK0S20onoJRRe6nfIA6pQmnDmLx8KRkL1V/GCkzozM8GP/YgVjDnBWL9ZkiX72RdRWqvBVUF7Vql6bgPB01+Key0bdZWeHoS0nUz53/j2VaQnK3UD5JAbIp97QO6d8sWDJJzZQFsq5FGx5IZ0FSMOnePF93MWYNacPOR3w41KAWsiITM5h21e2FHVbLTZTV5TAbFJbRHAAwLTKNHwqdpxJjTFr1NDwVHTiVXGmUW51+Zns8xq4TxEjT7/WxqYth7MvNOWJh67cyDTlAsr9iARs2y4s+TCAVxoMlFld4nDQJUCbjuVZHndnBw+gpluAWSsQi2u4Sr6epSFv2KzMLxDN/kUSO68FcIcG7cxZO0BeOV2ZwKwuUzjeEkCdvfyihRpalmGDqlqGpcSDL1/H7iR5frPCRZ5+rS3jlZ/IAU/c6feb/Zi68lPXlEO+OPWGV08HzN1C5G5uJkh3htwLGHefUz7PeLbblL+KD4PKcWAloV3DOQJwy/Zi1enyma1M1YhRkKGZ7XvGCjkQuHlYXbMojakz5maAOitjZsas3eKe1aSKCJkDA0/kgMaCJ3oAQ3BSxgFiRGfrGvZ5DeAeGvOiyP0iKgMjBnmhMy+Zsi62l7KFTFoiNDSlDNsuANwGKAcQkwM18ahxvABmKqwYA5gVcOeLH3OgNq1PTeNNwAYAsqpyoSuTyRwstJBNxiqMbILDHs+qyU5CkjGtXs/G1K7BdegLti/5Hsas6L0z5sKWl0WI2GOWOWOUA5RPrifPHmkh3ZhyOvdcsqBu2/kjPtNwOIux5Lq/iAgqTHmBZ0mJfbcCM41birkW865prgTVV3XxvaV4K7HniAImJ9v1hZ3a2GZSzMpgZbDXAzlpw0F7jmWwZntpyo4va9vntdkzacy60JqLthxAW8LjbB8LKWPozD5JKnsOluz3xZILrKBJEoypacZoWgiQllO8nTGHpizpFTfmrJ0GSBNnXGboxuYpovGFIPvxQB/fDmlDVZHm0rdLFrav44dKAapV94qD9S6zJd/uvX4oWzr+KGOaU1v2OOXq6OueKBKg3Csod7Z5IYMpJ0vunj0aIZe+Isss0zUwZwnaKxhzasxYAjQXkFZKh7olpRLQjNwIAaSeQxDs2UWHCsxz97RuYZzIxsWKH6n9cKHhRK61Ky0cgLfZPq/N7q6VUUE4wuVcxpDJwuNkIpMzPJFEJuSfuqZc5Qt5BEgBZDmGpqzQgyY7pkmTHZOH8LDHY7YmFzOZ7LzPZ3AWgtEA6GXGVcSWihDkwCZ1OIvW2ZajqgRq48umFBmAnh0Iu85wJFpIEly+AOKbw5N5p6EMmhgcxfSFF2nvFqFxLmlYWNHOLJ7HonXUqAQ3ebjc0ULmXFN+VY542ic86ROe9AOezqYxnzrj1FvKFyIEedoWgEw3NBjxPNixgTRGPHyJ+lmHatYMVGAbnJdsmYacQcGSMQB6xigSxgBNwaLtnGNbjzq2m2IGIJ3R/bsnCvQWNTSG/PBqOwIwJv1EDx6tYb38xKM0JMd/eTH7vDa7f9lPoCy7gfDuxge9SLWOkLiVfFFBOe5Ds0kAZ8nUjCETe24/K6ape/KhYGJJICbPVIqg95ijtdFk1h0o+10Yk1IC9cycGVizT15hB2d4dTyxyWMMmlDnkECNAbtpC5lEy7j4NtEyCWDh+KNbs/7sOnZm8VAWnUkAJFuO8ptRZ3lWS7EOptzd2beWLxYsuTvIeto+u5+COgykA3iDTYfW3AcQVwYNjNsKzhpqWoKzJkMOoFax75tHy1mZb7XVIGOwaCtqT/Z97ABhxNJrdzLhsgbBOERnsfHghlnFa47wchz9+5c9ES8Q531em92vg0lZcmsCcgHYBSBjmcW3kC/qvg5QZmfKzQC5+W0EvE9NLCuJFIcWwGySxgBqXQDy2mrhmJmslE3sU1eoTzLAsrFI7Hq7+vV2kzgCmzUbm5lVoE5HSyOo6vKHq1loXR3HGiu+GPcNU4xSkS/SiOibAfxJAL9NVd/wDM//3QDeC+AjAfwMgC9S1Rsi+hQA3w7g9wH4AlX97x/wtC9abY4aHUk6Rt++WZr/cUoap94wdweeSK0O+WIByjRqqeS+AW1KGZmo5PuqC8acYF315gu2iO0PHkWAkkUaZds2YLBl901CATTykDkAsHNVaIKzdgYg5noRhrBi7gyihnliNGHM5GNGgs5jHKt+X8f97Bpe0rx+vdn1P00OFlm1LYsODeZ8xpaLg2/s6wKUdbK/BOVJ0hPcJkn5gklxnPoAZmfIwZSzWHcB5i3GXMsqRnGZLnZLnpLavZSXsHpygH/5GBASqHLWtB8RPyFjuLNRx/XndmHNOWZeQS+6kCtF26nbtThVvCxm8T8B+G8A/B/P+PxvAfBtqvpeIvouAF8C4DsB/FMAXwzgax7iJO9jA5xjjpjjb45SnKDBloXRxXXTkC/C0ddXoDw7ELvMlVKGYLktAEUGqYYMVmSObKhQgLkCdOjKC415kCaP4jSAZoCx/I5SALNavL6o+XhojlWmg3M3kiLd5L7e2ciQ2Heok7HmHDPhHMutcd6ylzivX1d2RRzzUvNEo3T4xa24vpxsePJY5ao1Hzz6opmmnPLFpMBBTLpognboaM0Y8nHqyZIbe7UrthY4h2YhOux5+3E8OjVcslhajYljy61gQnH8ZhpfxLkzbli8xgGhU4N2tokKNn08lq1qrINngrjT0HgGgSfN5ibiccwssHEU067N+1I0/QuShoBw8xKYhar+FIARpeNGRB8F4LsAfIwf+ipV/burxxCAPwLgT/ih7wbwDQC+U1X/T3/M5Q/vgU1Us9Tn6FKyil3uE572iFtueOqxyvPcMM9sccontugLcU25sON2U8D5NLRk29YFY86U/64LB2AkOV0dxxwkqWxLG7rzIgrKa4mnzOgOaVv5qflZcpVLgDB0UvTWFpXunvZmfp6ueNonKxdKPcfxRiccisZskpGeSRova14DABH9XgB/A8AfAPDnVPW/eobXeBeMbHQAX6mq7y/3NQAfAPArqvrZt73OvX6aRg++IWlkjHFly7mPEj85JgMYC/kiQJmaojWTMBoLDs6QoxzhozYnCB95XoBxtMS5DzCH1zjAmN2jHOz7BKRTUZQ8+IIhYp5qgKHdfmQIQITCGUNWRKxohv8tximiTXAWjaFMd/iu7b1eZ0u+74Ax4Z8goo8B8H4AH7d6zG8H8K9Udfb9Xwbwlhd4jhdtsOWRkj2r/1CrL8Mjg8+TR1TYInY8umfIFlW+GGA82POosQIZ+8GSuS/BOLer1rxlVL9vFZiHn4M89FMaRks08jroAKhMKYvL9xcmK3GbfhYBlOy7ECvKLozORnJylaE8dGXczpYBvOx5/f8A+EoA/+GzPJmIPh7AFwD4BAC/C8APE9HvUc2mnX8WwM8DeNNdr3V9EaOxYwBTHQulmeoCiC5JGluacgCyg/LUOg6t4+BseVrVhz2ytbg5OGMewCxXaMw2cSIo/qSMiRyYyb6QqtbfrAvjRFYaMV9DomGmmJQTDnONL8G43pAxqI5F6MoLBw08nGm8z+2tpfCQzOLNRPSBsv8eVX3PPV/j0wF8fDn/NxHRG1X118pjtr6Vt0HNa261VdNwWo0efLWgfTj9ovZFlOisySMVlLloyjQPUKao2R33Sd3W0ohB81jqzXdozKkrF+ksnp+xymRfPfN/wFZ4MMcezUhU4E4jEaQT4PXHFTAmzcixMFCWRSOAI8/elYXTiboc9/OLeeB5fS9T1Q8C+CARfdb6PiL6QhhoHwH8PQBfXgA37HMBvFdVnwL4RSL6BQCfBOAnieitAD4LwDcD+Oq7zuUeURmVJft2qaVcY5WjglyGyx38/smz+Q4efeGacsgXrQleOZ4weTWrDz/cmITBHRN3PG4zDmwSxiOes+XNoqEkjW6/a6tdimtDzHDyRPnGkzQ84cniVj1mlQD0qWPuNmmYFZ39y0lDc9ZuTFhAKJUQ7TnhBBQbD5KxRIwKfsmeb0kuAeBOktsjN+5hH1LVtz3nazCAT1bVV+tBIno/gN8BW8J9GYCPIKLJWfNbAfzqc77vc9voYM25fQqnn/Ii9Xr2+hfSm9W+OHEmj/BswMZzaMxFvuhAe6rp9Gs3A4RHlcIA7JAyondksOdyu7IMkYvb7P5jAMpthMFJL99ZJxK2cgUAj+DwcgisllglCvDJE04A6Ik9urlhBtIJ2EhxI5NVoBOP+ab4sfNxVt7+icaDz+sHMSL6OACfD+APqeqJiP4ygHcA+J7VQ98C4KfKfl0RfjuArwXwxmve8+oOJkP3REoYCy9wWT7Fcmr8cmth0P6zz8aUeSVfTIUlh458bDMmErzSTsmMX2k3CcjZ0Rej/XqY5euXX+r8ElJ2MI5uDABwKEw6mPdMioN0sE+YxvGLzpCmEHWHoMS1kjPkkYq9NT5rL7oSLarN3eYAVH15zOKC/RCArwDwrQBARL9fVf+hqv7R+iAi+lEAfxwWmfGnAPzAiz7RaqPqGY9ld0YS0IItR9w7lLL+xah7gUyxzpC4tXxRWHLoyAHCfNJRMGtWY6cRpSFqv+7RHDisNFUg6FhxsTNaMsmMmLLwVq7S1J2Aqp5Y5WFzfRAJsmBjgKJxsDcK9mu1OH+b710Jk1IJN1yOY8hDUec66mWs7TWY1w+xGvw0AJ8I4Kd9RfgKgA9uPG5zRUhEnw3gg6r6D4joU695w/u7P+svM48swKpnLfTllDX8k3enH7XRVYF5gHJ0Tzhyx5Fnly+snU2w5APZdoAyQ3HgOdusx+3aMrjdv4AsakW+tWdQ+0ktg2lull7KpJZMwA03gFfaMhC374UnkDSvUtewAOgzScfTsBdhcsmUcSsYV1NYrd8XbUT0F2HOuw8jol8G8NdU9Rtgy7x3E9E/gs2rHwfwpzde4usAvJeIvgnAzwL46/66/w6A9wH4bQD+AyL6L1X1E17r66nWC2OWNciodR4RdxgPUKYByKuQuHTsRe0UWYP1kCqobgd7Lp3YUdmynAPacNAbgFoonKbkwBhaM3d4UpVC2qhmp2TnySDzQZeMVyILJw3Qj+tWoQwb7CLL8cpbTjnjLnsN5vWtq0EieidsJQcAb1fVrRUcAfhuVX3X6rmfB+DrffdLYQz5o8tDYkX4OQA+h4jeDuAxTOb7W6r6hZfO64oOJkuHXzr+VnUxIuX4LJnEIzF0cinDM/pocj15suiLw2Sa8ocfbnBoBspvODzNQtwTdbyhPc0eY9FfLFqmD5CWbBzZMMA54yd1sKMTt2xJv+52zCSY3St/9MiIg0w4cYMCOLGCyYqgd+J0DOocgjNBZ/P8s5bloQAivtTsNfGEzsZ4Mf4rU7UiMi/aVPVrYUuy9fEPwZZ7dz3/n8B0t/Xxn4ZN5Bdm3RlbxC2H4y8iCSKG+SRtkXqdMcveMMEkDHInnssZnkTCJ5cpuskX0ROy3ag3UFCXNgTZGzIYs6oDumSMc4ZCVDmjzJV0yEfLYTuzIAAAIABJREFUtyBQfXxnDakVMhVQbvDyADAW7I7s+B1g2HdYyGpr9AOA2WQ8ZbUSp8xovQ0ZSEZ0Sx3fEcu8zZhf5LxW1XcDePcdD/sRAD9ARN+mqh8koo8E8EZVfR+MTAAAiOhVAN9LRH8J5vz7WAB/X1V/EsC7/DGfCuBrbgNl4L49/yI20ufB2nGF8reQMyiW9e7wW6VZR0jc2tEXoHx0dhws+cAzHtOMA1mUxjGAGZKMGcBCZ+4RSxmMmbztulrERbO1HVgsocVKF4bjY8KxWaiPMOHgTg5tVhJRolyoX9uIwggWvByP9ZjV8QwGvRj3TaPUu3d7dqsx7cGY43jE4xoDDJaJLH4FJa84OAoQRf2LwY5rhxstMoYmYCdjjsa9XVxXVqArSDzdNMAZ2ATmlMKyXK2n9DfXynzFyt3i9RlWWY67AS555AUoztsIBTFcvkAWPorKilBn5T4+oudjZ+M5ihnVcT+3lzeviehfh/lC3gRAiOirAHy8qv4cEf15AD9ERAwL2HongF+qz1fVf0xE3w/g52A9B9654SC8yu6dYJLbK3DOCI2F5qwDkELWIGSaNXtWX4TEpYThmvIWKEfrmkd8wpFmMARH6ivGbIDcyi9yp/EFFLIJ0szdAyGrVQEBGtvtI54jIg6ABb0fvYXOUxY0YWgkv7BCxYBZV9caujqVsdFYchZQxlrGuCLBpEaK7PbsJoWhhYxhx8eyXMufpb4FOIe+jEyfXv8FSHNx7C1AeQ5gdtY8iwGyFW4ZIF1ljA1gJqYRB98s7EKJSmsqm5w0q3/zR53wlDA8FI5odf4UoIxFQkqECtbxqeMW49jLGN/q/HuJ81pV/y9cWLWp6vcB+L4rXuObYZEXl+7/MQA/dtfr3F9jzgw1JJhsMuYVQNvxc7Yc2Xvh6Jt4aMoT9U1QPtCMx3Rj8c1k2vJjvkkgPtCcpxvOvyixcsIEoBtQy9E7+E4O0MFUvZCKh1QIEQ7eHVj8/KTZQoxdIw/WDAd2hONF7hofWgL0HdEYYQoLSdrt4WxkqEXG2hqUBz4ugSlAbFXKswLy4tjQjlNPDlDuDszBllWBuQ8glrEKJNGRWwD4sgwAEQjNGv02n1z5Q+/HPCwui3D5PIzzDY2ZAqjZr0+cbHl5XJQxGc1lKzgPbVmukCj2eW12fVSG3w4NiwojrNtFZ96oiUGeZj15Vt9x6jh68sixWUjcK+2ERzzjDe3pApTf0J6kbPHh/DRBmSF4zCc0LB1/C8acCQQjf/+Anl7ik07WcNKzlAAr8n3wHwdRtkzDfsDNNGch/OM0ViqqgEwWoaFqTEQtomhozH3IGstEExoRL3foy/FetQzjbs9mOS8wIjJG0lHLdP3obi3Cy+pxGX1BI1Y5QuZm9W2vPjirRV90NU25D6bMNz1lCzp1oEuiHQUwbzBmAsY8KX6KRff15mIxcxbZ0sYj25QY5GGt3OKlCTIbiwYcqGeYE5DV+1e6nNFjbBRdJcsczNIW4xnRL3Xc17bPa7OrGXP9ZY56GVk3Y8UINdjx4rjm/Vm6MzRlZ84TdXPuedxyOPqCKQcoH1xfPmCEyh3RU19uOC9kdIAvqcgmW/P4INOSbWl1KGlPB+qDMTNh4p5ZYBN1zMSYONLFbSJyqYFRr598P2JMF+OS46ijToaP913ceZcyHsYqSKQvIpbhqKwZY6kTMkaw5SprrFjzYNS6SBqp8kVqyV0NlF1XptjOyAxdsOY05jHZ2H4ktHEBbUshUbDX3nCdzpkxKM6LMkMw2PKy/Oi43rFyCDkDgy0jpMO4LXLRHTN7n9fPUijfbcvxt+UYHKFhMBmDDJSZRz1lhqVZT+74Y2gy1QNXMLbbI3Uc0F1bnj0yoycgH7AxcWHncAKjOUBDzUF4wgQoFq3WQw7pbJlNB+qY0VxuEUxqVejyGlhAZPu1EaYV8qdzB9/G2G2N8yVTUFbB2+35bamBjuxQUc90W4AzCiBhAHEBZ4qY4wrKi4y+8ldAmURcVx7gnOxZBLoFyoA9h3kANODSmrcCFgdlCBCMlIakYrWY49zCeejXEZ3idQD1uP4qZ/gYxRhilP3cGuct2+e12TPEMccfnSWSjJhljNoQHrec2rInk0wsODSLxIjwuKPHKr/SbvCIZzziEx5T3N6kfHFASBkdx8j88zA5A+lxujXzDgBOKllE5eSJJDfUcdAGiLHmk3Z0JTxxhyAYeFQKgR9Db24WoWHdSAjiqeXkjkA091h3ZEr2Yozq2BU2fZVFdMBuD2aZig1GdWD1rI0x4pcziaT79ipGucYqj+QSj8SIvy6g2f9CvhABnebBmgOMxVGxl9kc4m/sNkb056TuIM0MEgdnZkCcRSdAMqhpfk+4xzxUSKcRvVlYdF57dH0X5NhEevZi/DJZ5wrA3ec1gGury2USyfmAraNezhgz/DYiMoBkmVE/OVgou/MvMvoGQ57T0TeY8hKUj2TqVaORxxIW2x3AgWAdFZwxn4hxVAHIGbPaeTZSHDFDyLSxA3WcqPn5OWsWe/96PbS63rj+xbjcMn52rBQ1uhAyp9i1uOc1ATIjbWt5vQ7p0rqv45YKi0yWuWLLI3ojMvk8RjnD4gpTXoNy70NfVlklmETjU7JNVpMVsPyNJy69QsKxSAolzWw+Skemju28niLfxPXWcViPz8b4YTXeW9x/n9dmtwPz1g+XA/RSI8WZzrzWUbECMCaUesohXwgaJEPfQppgSEoWpinPZ0z5UAC5EW3mGDUAXe297UtpEz3A+URzloA80IyulpbdVFPLnsE4kGCOLEMWsHB2Uok/DY17oTMvx6WGyo0wuvNBp/jSVduZxXOZbhzLtOGikYZmGkv0sYzHojYyxX69vwK2a8RD5vB9LeAc+vGaKXdJQNZLGrOajGaTzZiz1hJNIgbeFZT9/YdcoSlrxPnr6npTW8ZaZ0ZKPjFewNDoL/74reF5n9cAnkXKWNuaGddjri9XkAKVHn1AFravf5lm7Wz5GPpxgHbRkhspDtAFKB/KBMjC334/CDmzDlDfN3BuakzdYpwVR+oQZ82jDsf5+cYqoMXMzR8s05dDZ744Ts848Nr3CfxaWo3FzbDhZI4OXvHghb6MBZBVh1mmWQdbFhkgvdaSRWEVhwoo9/ILLTJWVL1DW7MfAxZ7HtoAZ399i2tWZ+I+N8VkCVIq57m6jqoto3yNUmdejlMdu/vZPq+BZwDmZfQFLY4v/uCP8xoZI3bZQDmiMZrLAocShbF0+kVkxg0e8wlH+DGXL4IpPyZCA4HjtnJmP02BAb5APR3UNOmT/9wfQspg4IQGyNGylYhx4DmdgAfu3gOO8/y7EmayyAzxazVN2UBZN8ZojB0loF89lxXQ3UnyoDZil5clKmtiiTnBKLtbZ/RFAa8ar3zu5AspQyw0riswdwuJE5MxUr5wpqyneQD2ulZGaM5M5uBz6ZHsgoDm+rI/jmZAmUGw8oZa46rz/JZ6Mup1pSTjWrODcoC0FrYMRA/FMa532j6vATwAY17LGsv7yk4ljME0nUkzdPx5lEbIFs0z+SwMbqRcN9JMsGvAGSi3rcgGdRcHAdAR52zatEkWnQRNaby/jvev5xd/W9eTOnMZh7OzqfLF5sL6CrvgoN/tYWwkSsR+vXO1nct7PZM5Fn/iLxS39U9Wt6Enr0FZNz54YQ9zcxYtasxZafm6tH7POB86P1eMh1sqOC3uW4/DYMtYyBn3tn1e3xOYr/0h25I3YpdKKYgzWcBBkEJ3HrcAFjHKBpxFU16BMm+dbNQ9dHBusGWYEKzYkIOzJGD7+6ucn9vq3O16NiTiZ5Etrhlnxb7ke0C7is0BZTmIJQD7/sIptpIBsghRkTQSzarGHIBcNOUKyrpRXY5YFuCsrmfnCYjY5Kzv5xKGi8PjR6XIMevryffTKu/g6qXeneO8z2sAz1orY3F8tb3eB7B2/MVdFZAT7MKp5qy4Ovzi+AFDwghNuYLyVGIyGjG6MwxGg7X+NXA+OI3vas7Dg4qFS/n7dZA7/uIHYdlXsAJ0XQ0sHYB03RhdO971IbuT5EFMLvwS1u4mugKmNTinJrs6tnAU5nEdQOiyRrBadTkDKqkpL0C5asyqOU+0wwrdBziH5qxiTRyaxznH+3kGnyYgl1hl3bgGnGvPi3Hw7bqqkAuT+9J4h+3z+nmkjBUrXoTS3cKYASQ4b9nodr29nqkZfQ2DXK6ZcivxnbHdVez+wpxDCol3Yw+jO39fMfnjwnldvK4LjDl15Y37rjIlYGcWr4nVmFsjk/WDwzYoud0ytUdluK16ysAy2iIANLYrU67olyESVrQ+mLNNamfJwVOqs3DxvmotpvTyyW9e1wU2nREsGLHhV9s+rwE8RFTGM9pY/t8C0kVftv1bZv0DWn2/0Jkv2TXX8ZqYYp/Ar1OjW4A7req8sf8ibC2d3Pa+hSm/sNm9z2sALxqYrwSvLQBeH6u/+7xa9le2vD4upTwqE6H7xFw/w6I17j6vTXtBIH0Led/tRdu1H/kGENKaQVcisI5ZvgSkqstVlwjQ2vnr+fudSb3h3LvLXsDU3uf1iwbmKx0EW4Ho62OCskJbTcqusgnOfTVBpUzy9Vy45hwu2rN6o+9puxb3OrJrP4qt7Nl1wSpijIw+XqZhVwfeba9bJYvVd0G3Sste4dOwx133sOexfV6/RCmjFtG+ZNHxoJPVU651lV9LCwAe6aOXnRXXXMdrYoo9rOh1ajWj86Jd2ULswa2+34VM0/FYu3mhU3uf1wCeB5hXThCqS6E7HCS3xTeO1jMXPOUlKNhym8J55w49EjBasuMalQF/XFfNVFDBMtv5EsBmEfVbwn02r+tiuJFevO9ao3Wa9m4PYhEZBCCzOtPucGxvxqznfQSCXm6GwDyki+jN17ttC1u9C+9abS+oC2AlJkQRo3y9+l6XWpVl84vL38vNr8UFx3bWjAEykuk+ts/r+wDzRW1rtb0ZUrRqzeN31S4HUmoURDfrruS3jBMmP644gXFSC5nrLmNw/NRarFBGZ1RNuYJy94w/2wdOaiVBo1bGSafM7+9KeXz0gFt1uchhWrUgunaMrh1vN1JYf7bdntv4AkWrgEIr8FlGJI2/9bHN+iilKBhlc1S7k5itSlw3gLM0a/EgIiton9EZFQwLKBMR0JoXwmKgsYXKRRKBA7aW89iqd7O4Bqyuc+vHaUXALwHypfEG9nkddj/GfO0S4xbGXDNKa2nAkQob7ePZOvJ5jz6go4O8yL0BNntBIkRSCPlJRqXvlVVQFjVQ7vDa5KhdfF2ewPI8Fue2One7ng08fRZWfOU4706Sh7NL4ZlnVmqhnKXQrwErgCwqbWbVQG+KwLCMO2DIChylO9WAWmAp1BWceeNcKygz+63/ES+qROb7cYRtIgEafr6LgmTra8MGOF/p8L5mnPd5/QAa86Ik4Nl95XBduW9VoIo/HW1ounezDp25gz1tWtGV0EgTw3q8QQXnlS1AGQbKnva/YMXRzj705VGmcJxf/G1dj+a/MQ5nphjJBs9iui/5XmuLVHsuWDbuXG0nqwzgXd0XfwHGcVu1XmfNeUtsadUBtgWcz6yAcr4W1e1yu3hPlNvVHyrTj0HYuPYcrzyVRYmCe9k+rwE8AzAvs5mGrryopOWkQhXwfOdFIW1V7wlGiq6MWdhqHZPgpM3aPYnixM1qVUAAOeIAb/ekwMkTPjqsuhw8XbshqstFosqoLgcg5YuOkC+8aD7Y+v2h4Ykc8ESOeS4nbTjJVLbtbxbO8+/CBtCLgurIGgS0MUZjTLWM6T0/i90ezEYmp1UtDMuMVVJnvqbtKivI+94FMw7Gac0QaPTBZHgPSAdBtgw9NEU2TvUewtRdBCAGxNdvnp5NvV9OGAn5gk3KADcDbDY5A8zQqQFTs3PxQvp2nuTH1n/L69LF9Y3qkeQTnFaA3DJLVq9elbyseU1E7wDwdb776wD+jKr+b/d8jXcB+BKY6+orVfX95b4G4AMAfkVVP/u213n+qIwt2SLBu1SfwjjWlTBpaUGz+uvKELIGqayaNSvEwbN555ETMaBR98LfPDVns74KiQtQFixBOWSMkzZEY84bb8x60slY/IaEUVsP9RTh7P2iCtmiBdHWOD3juO/M4rW1Wgcl2XKwjgDneHAwSwcvKnqsVlCLVOhGvtLkbJyqjbMdFAHQCAr1gkQkgLY2HISttIOo8gXxEpSdRUcPQA1m3QyMQ9JYnudSF68ddwCftrQajzJOdezuZS93Xv8igD+sqv8vEX0mgPcA+IPXPpmIPh7AFwD4BAC/C8APE9HvUU1H158F8PMA3nTXa90OzJtL8CiujQJAGGxvwaiL5rpyAEYFKmtDY061k0sHJ204aE/dV8idcSCwOwJvqHtxewPnEzRXh1110cEkLAC565AvApRvsHT4GRgPwO5q52XnN1rljDY6uOj42xqXBTuussaGtLFZtAY7MD+P0caxyPSMmi3AkDMW3WkKYMWxtbMPKzBTlxGUffXoerNSdKoOaYGTBSc4t0B7tdoXW9ETNSLDATpBuYCzpWivHX/RJo6y4/36/C85A+v9tBovICJcdJHBuzjtVSjqy5zXqvq/lt2fAvDW2CGiLwTwlQCOAP4egC8vgBv2uQDeq6pPAfwiEf0CgE8C8JNE9FYAnwXgmwF89V3ncjdjjoaQUexkZesfRSpgvag+5RW2gl0OgLbO07MwJrKW5wfqKRs0aLLmk04xW61HH6GAsy2bxAsSZe0LLH1pJ4345DUoN9xowwnNAZqSLce5iFpLe/HznZ1B1+vR1fXG9ddxuW387JjmmG92q4jX3J0kz2VWJ+UW0Fh9OGdhc34bzr0BXpTLe3VJI7LtQvqASwcUaNdGbQ4qGrJPdy/ZKSNGFMC65x9qz78KyM0jPULSCHBuBuDKcRvXU7bzenAO2HUc1uOzMX5Yjfdm8N7rZ15/CYD/GQCI6OMAfD6AP6SqJyL6ywDeAeB7Vs95CwzQw37ZjwHAtwP4WgBvvObN7y9lJPMbfcFyMLOAtt/6skSZgNBeO6F3xsxshblJceotlz5PZZzSUzkksgqRt5Bq1ptPrEffiebsPNK8dOfhltoWa9ki2PGNNvymPjJw1gm/KY9w0glPdMJTOeCpTHi1H/FUJtxIw41MOPWGkzBO3XXmztA+Gnaie9NKXY5LbtexKzrztfYymAUR/RiA3wngVT/0Gar6wXs8/3cDeC+AjwTwMwC+SFVviOirAXwpgBnAvwDwn6rqLz3kud9lEXPbsCzr2liswUM2QPCU5maFgwJ8TUc2vMzf1Ub52XKzPncMeGlLjkGxx4SDD/BiRPZCi84m7ZZQs8qMA4Aj/O4wmabcCHpwjbkxdCLoRJDJgFoa/I/sWupf6uUKberHFChjwxzt4sr4IaoyXqkxP+y8fjMRfaDsv0dV33Pr+xP9+zBg/nf90KcB+EQAP022WnkFwNac31qIKRF9NoAPquo/IKJPveak7wfMhTEvyxiutlEYs+vM4QiMFvAiDGUx1gwyBkrGQg8OmuEEbGwvetIhUByoAc5+o/OIFbm3Zjpbv9Y1HE6UcQN7D4uTbgnKp3pbnH6Cwe5nMTkjmb+wfxlpOPyKvlzH5dLYbY3z5c/ipTKLd6jqB+5+2KZ9C4BvU9X3EtF3wb4A3wngZwG8TVV/k4j+DIC/CGMpL8Rq1cAs64pl6zBKJ1f8qcsDcJ0WySiDNec+E7Sph78FS8WQvJhBECgs5ljtRLLziDFmZ8Jbq6gqVXAB5ADnAGXmTScfKMB3gy2X/WX0hua2vfWqZRyG429rnDft4ef1h1T1bZfuJKJ3Avgy3307gDcD+GsAPlNV/2U8DMB3q+q7Vs/9PABf77tfCmPIH10e8lYAvwrgcwB8DhG9HcBjAG8ior+lql946byuBmYSTTIXkkZKGysddQCSLgA6mINFZgBdGF0Us1iLptkjHib/UpykgdkFYQZudJzuSTtAyPhiMNDUpnaPEp2FfmZsciSNgPFEDunoOxXZ4sZvn8ghozBmKdEYKWuwXwP7j81YQSwAOX+g1mN0YTx9vO/8TF5HGjMRfRSA7wLwMX7oq1T1764eQwD+CIA/4Ye+G8A3APhOVf3R8tCfAnBx0j601XkSwBFRGQnOJfKgAnJ1mIV0QQHU/hdVZhP01BueNnsHQECNDJQhIxxOPd55xmA6Ysx68eNd49RCO57aAOo2QBkByq0CNFCjMMaPCRKIF2C9jnMuESsZYuhjGuNYAfmuYmAvcl6r6rsBvBsAiOhjAPyPsFXc/14e9iMAfoCIvk1VP0hEHwngjar6PgDvy/MmehXA9xLRX4I5/z4WwN9X1Z8E8C5/zKcC+JrbQBm4FphLqUBr5AjPpqDSnj22bWCpu/zFsR+OBoYcGPM8Xv5mav7yhCc8Geg1Bnv43COeky1bD8A545ije/YJLT/wA40Xb9BF8aGTg3sH4YkcE5S7UsoXJ234tf4YJ214KhOeyoTf6Ec87SZjPJ1dzugNN3PD3Bnz3Iw1z+wShssYdTxi32WM6K/m4Sku/+hivC/ZOuTuBdvfIKIO4H8A8E2qqgC+A8aEf8In+PsBfNzqeb8dwL9S1fiAqgZXLfW919pGe7Hz9mETdzRu7rxSdBYwE3ozgNQGWwUybGmvBJ1HAhWDIJMCIBDB+qoSwB2AMrgraFZz/sG/R52HtNH9+1aBuLDl7GgdFqFzRBYS51oyqMgXTJBjMxCeCNIIciSXLkLSAGSCSxy+3QCd/DpDsmkmY6ApmAXM4uFx0dOzn7WLq92Ituwlz+v/AjZH/7JLFrOqvk1Vf46I/jyAHyIiBnAC8E4AC6lNVf8xEX0/gJ+DSXLv3HAQXmX315gl2B2G1rxig3WZrgrXog28IUivtAhDWNGFcYItf2ZxvdmZKetgzE9CcwbwRBVHGECb5nz07tqSQFxLd2ZhIs/kiyiLGw+Pq46+YMpPXcYIJm8yht2eehsRGcIldjloBVKDr5LFenxGq6Hl+F5jD8gs7qPDvUNVf4WI3ggD5i+COUE+HcDH04gYeBMRvVFVf62e8sbrLS7Yvd9vA/CHn+E6nstCA43ektmdpjDm+EvGGM12M67XMvpq/LI0T+Rju1pVAtiAHBOAGSY7hJta2UCZrKMJoSFF61InI3OothhzcwmjOPqSKTsoh5yyjreWtjz/kGlSW8/rRpExlnLGol1c5hS8FI35alPVL4VJElv3fR+A77viNb4ZFnlx6f4fA/Bjd73Os9XKKAV4lsvzJUiTDG90ptixAbMwgTzhpAuDyAD6RowZMymedj89f/7BnX6dKR2CTArB7FmB3ptPz5dPUXyoZveFbCFK6egLTTlA+UYmPO0TZmXcOGPuYl2UuzDEr0HUgVnLtWow5CVAbyaUBEhvjffm5/GgE/hWHW7xtqq/4re/RkTfCwsH+h7Yp/TJqvpqfTwRvR/A74AF1n8ZgI8goslZc2hw8dhPB/DnYLGkT5//sq6zChitaKLVgbVojUaagJxRGVqW+YqllKEFoEFQVchEvtw3cFZx/cPfm8j7/XVyLcRXVLyxmnJg1pLJF7JFgHQCchugbMw4GD8lINdzX/wViQaLH6bl+Jz3wtSVhn8LQD/svP4ta9cBc2lxs3bwXWSEPpfyfgm2jBGh4ax57gbMJwfnmRQ3veHIDDg4C5mTT9iAEIzsySdE6GQ6NevI2qqdR6J0Z1daMOaIjzYwHpryU5kwa8PTPuGkjJtuTLkL4+S68tz5nC0L+Q/SSlOW1fjocvxiPClkjTruF+xFO/+IaALwEar6ISI6APhsAD/sd/8QgK8A8K3+2N+vqv9QVf/o6jV+FMAfh0Vm/CkAP+DH/20AfwXAH7tPlMfzWrQSi47sARqLbuikTkadNXvGG8hJh8cmh3QXUWyZ7acKhGNNFdJM5pLJNOwRywwnzWzhdhqsGQnKuiVzlfoXCc7h5HNnYIIxm2wBDlB2GSMdlOW8C1vGYtsiU8L5R1z15fOxs/EcHe9z3C/Y6yRc7qXa3cC81YpG1fVRLbqygjt5rObQVpkIQurpps4wZvukuy/rblg82cQKhh+k48TGnI+t46lMOHCHKGPijgN1POKDhci5o/DAczZMjdu1heMvbk8yMvoW0oU0/EY/pnxx0xt+4/QIN9Jw6g1PThPm7tryqUE6QTtDZ9eWZ3PYUAd4Pteaudt45dj5WK7HeDH+Z5/LS2EWjwC830G5wUD5r/p9Xwng3UT0j2Dz6scB/OmN1/g6AO8lom+CRWL8dT/+rQDeAOC/cznkn6rq57xWFwJ4l3R4C7NIhIC6H6NjYvs7cMfUOroQpAmoqXeiNqQKgBUl0IRlGQwB4Nl+UEDd9wDYPNBmq0oQj7nQ/LsVLDnnhq+2tuSurBYHZKJIRIcsgBmQY9luhH509twAOTh7nmIbqTlLA3RSyAS7yElBPh6tCSYWTM3HK8bOx3LR3T515g1la2fMAJ4pjnk5QbLleeqq6ozRWXPG7VrMp5Vzi4B1nxzC2earTx0sjBsAB7FSn8c2Q5QwkXWxnj3C/kQtwVnElobNu1eH1LJ2/kUUh0Qcs4OyKOHVfrTj0hbyxay8kDDm0Je9PkbEZ6NTphfGCgEJvgOER1SGno3jVaFy8O/fC2YWqvobsHjOrfs+hCvC21T1n8Dkj/XxT3/uE3xOC8YcckY4skY8szu5KGKZXdLgEWURTsBYGVUJw3Rley9ttgq0WH6yDtdkr8Nw1i0O6L7KTD9NJpgUzRkYJUWD/UaoG5tTL5JHkiHnbZEwiqSRf3Fd+bhx3cQ+Piw2PnW88lYWWvNt9jLm9evRro/KyIh5ZGTGYikeeupq2Z4aa0QgeCwzBNDOEBL0Emw/d5t1ooQTN4h7nYU7uB8wK2PyT27ijhlsSSXUxySAbsZLBiBHzYsA5A5z6pl8wThJw42HyAUgn3qgSNN7AAAbjUlEQVRbgHLvjN4p2XI4/eJHKcD5bCzK+CzGLyIyMtmksObNzwTgfh2I73bZWokWGBlpNSligAz5LUhtVZjg7FqyA7Symr6rBmqAfazURugyNbKlv68i5TBWVQInMuGz6YMMaZAiAEAFZQfeuG1L8qORMJLFisLJZ0AMRiaXLBJKPIkkJZsIocukEgCr8UlwLuOY6e6+OrkYMrfPawD3Ycy+nDJQLnJG9wlXQuSIFDw7oMIn4wQrk6k2KdXBSJXRqUFETTsG0Lhhah0K4MCCpyyYSHAzzZioY2LBkTuYBAdnN4cMzZFb9asscu/M2ADa6nXceMW4WS0krmrKIV90YZxObWT5nRiYXVueCXSilCv4FHKG/VHXpZwR8kYsVZM5X1iurmxf8j2MJUsujDmX4yI48oxD6+ieTDRNAov3axAIdLaiQ0IAdbJbVq9saN+PkPLUZSzEY30lxS1+yE0SrOGUscpaOIs3bJk2PXRj0ABhTW15AG/KF7E9+X0HQCYdcsZBc1VABwFPCm4d0ySYmuDQOg6t48gzJpKUgYIxx4/fXdEZ+7y+EphVR6uYSBrJUDjVwgCHhFGX8DYJkVlT1P1XHQC6MWdA0MGWUOJ2YtPtmjCkdVBXzMSY1DTpiQWzs5pZuTCcyx987ZQSwByZhzcyWUafUurJEX2xZMoOyt1iliERt7yK5a5jUOWMkDEWUtA4Xsf98oeyL/kewmqXjUaacoatuuwH/0ZGLHOkZzMbQaGIYYbNf20uNwAQHTHKFnkRRMVlPbJj3H2BxL66JGSM8mKFJbhd6fKokATmAtIZcUFIPTmjL0K+aAWUGyBFvrCYZb9tzpZJcyxaiWEOotTCCZiOv3HyF9tN7fMawD015sz+CwdgLr2HfhrRCBWMYsKRT0Yr7WaOPiWUugFioXCwTDqmCIqX/BgnFtN3G2GSEW/6LMAcICzKGZucIB1A7NEXEYHROyUoawFlA+gKxlTGZSljLBNLsCxa5ON852eBfcn3kJaJJjSW3RN1MCZMzqIbC0jYkykIwgJSD1IWpFSApgBCcx7vkbKGT3zTl30flN8V0KhDsyV9XTTCAoxriri0EdonAb4BygWMh6aMBSgnOFNoywDlj5SYjMHi42Y/bJP7f8LxV8f5lkvY5zWurS5XIwa6fzCi0G63PCOZcKbrz16zAkA4QVgIOll4UHzwEJiq19VDiwjS7NdY8ktgH+xx6vblIMXBnTBRMKWG4dT4ybyMZDGjFVRozZEocnJnXhfGzexsWjyk79Ts3EK+CFC+cQnDve18gu/bNndbxvKsGalhUobaMbExjVsPVRlyxm3V5fYl33MZExlwqLosNuPoDuX4e9RmCAjH1s3R2wg3bQx8J0AOHgUxM/SoLlfYPADUP28DSjojK7Yvs+aPt/Tyo14ZZJUyKnY5KV+U5HRgriw546yrlMEuX8T+YTBlkzXUWTKASYHJozCmbpEYk8kXj1rHkTuOzcbsUZsX43ikGQeaM2rK2DWffyj7vAbwDAkmpM6aI6bZAfWMNccS3iWMkDTQqTBCn0X+Kw2FpWyrgBjoxBCm/FUGgMbmNNFwnojpe5MD9l0FurO4vde6UD+mSjj1AczBkKMjifSIviiacl+CskkYtGTJK0ljwZZ9/Grs8qg/cgdr8Ofu9vwWLG6wZluJTSxg0WR/M1sSU2OFiEKbFeHipsZdmoGY9cK211KXJoQUrB61EwEWZV6YvOHSIA2pMBO17qExV0DO/YaMca6MGSumnJoyL0G5ShjcXL5o4vKOs2WOVUaktEuO5dY4b1/EPq+Bq6My/Ce7ShgK9xpXACqg7F1+QQWk3Za/iO4QAYbmrAw0taQndzZGdEZ3Fq1KGdBOZD0As9uEv/IWY9ayn+2gFAnUUSVunpsBsg5ADumiyhd31cSw47o6VvZDb9ZzSQO3lC+1Jd9Vn95ut9gACylRAxE5YA6sSdlidFXQxZxcUUs8lvHUAECMqQIJziZf+NzzCnHo9jtszVjLxxwgTFj4HlC1Zb+tAB0F6lBuM1PPIzW0aMwRhxyPGbLFUlNegDKPmGWiETYYTr/Jf8gmr5ExcV+M44h6WYL02vZ5bXaPqIwyOxJAAFRgSaY8ABoRp1kGm7tVgQujDkR9CXUN2t6KM140QoWYFURsgf4lRGcmTqAGcFa020475AzbjnZQud1Hlbh1Rp/OriPrYMmojNiXq+cgvU4kCUA+H8PF+N6R9RfOw92e35i0pPEPP0XUEm6kBjzBml1Cm6zyLLiZS0/gIAYMcPb4ZHSP0BD7PmSsMhlAi68q1QF75AUAmTdQQHn9yVdwrt2uK2NOwC6yRjLoVXEiZSxAGc0SX5gV3MRkjNbT6Rds2fT4WtN6xIQDWOjNm7bPawBXALMmQ3bNswusv5iAiaBdALJJOUR7GtJExOXSyC4CYCFzEXqjJSW0e2xoA3S2SRFpoDIJ4AyZ2qiTC/cOA8jU0MvXQ4l/4rHH0XlEO2VygM4Rm4zCkP3H5zRWARYSR0M7Dl25A3xSD5Pz8DnXmqOqmG0LeBbQLK4xLzV9jTXv2YUANO8T+HmsuQ8kMkUPNONADY+JIER4zKdcab3aDvm80+RRGsxgQv6o986YgYVzWDl+xBV8itWWQicgm0mIl/ZUDE052HLsR6jpFVIGsATlMwY9OTt2n3tEZIDVQuLC0TcFIBtTng6mK7cmeDR1HKcZBzbG/Lid8Ni15Vf4Bo94xmM+4TGd8JhufGyHxsw5/rS6iH1eA89SKD9imYszkHjE5xIZQEsHyH+5uQNR/I6BhacaOtJVnWYAXqQonIZRxEiVfZJpRngYMAPCI6Rv2QIoaHhcwpA0ouhQgnFWh0Oy40wUmWN7JVXMtDjGRc7g+rjugJxsujj8FhJRYc532B5W9HAWSQ/DMeWp2dzRYUlNnTuOIBx6T1/F1DqmPpxY0qOwpRhzFTanX2HMYM/689/cYNFjH0tgju8GVsAc35fYXbDm1W2pdZGMOPYnHZpzsGRSc/TxYMoBypF6HUz54E6/iYekcfB07AwzLON7l+3z+uoiRjpimdc1M9TuJ+9rZk4+cpCyyAwALujHdnlttXoSoW8JomyiPz6AmbwYjOvWmXEVunKUIcTdjDneV6W4swUDlNVZjqIAs92igq/QEqSrxrxIJjmvi0GxrK0yRgFlvSvJRGNMd3seawQrLQsMHVRHQaMDdZzQvFYLQRx4pJTZnUu7p+6yhhKs/vKkHnlk01PF5h3DCx85aBKXfACXAeP3OZ3FuIIxAwOIy7b69yMBOME6al8Uppyasjn6iIp8wZ5M4kzZEnEsRO5ABsgN4qAcev3Q722cFY3Oz98uYp/XwH2cf6JA76AuUOrOhsVb3zAYAj0F87SmCeKxy5HRp01NytBSxYpdTvXUUTr4JGKCzih6mUd/BAtoA4iVSm1aDLDevpYyI8S1PCBZcrIUB95cXpYaGHwaX5bYTp35pEXKcFCO8LiTGnOe1eULBXWTMagrMNv4ort36C7n37xTi+c1hunKBwiO1CE4oRPbMlzmqBRgbc5KMtNEVhI2ygOcmiUhESl6iejprSUY64nTedwPY34hfuDL3EsgBhAlJq6OY/YL03rcvy+pNydQa6ZYgy2jL1Kt29QzIirkCytSJHhlOmWW37F1vNJOeMQzXmk3+LB2g0d8wmMyOeNAHY/phCN1HFyD3giUy+vd5/VVccz+8aoAtWB3MOauBoTBEoO1dmMFFss8Oi0YEfZuDxEmF7/ePkeiMpZ46FDEfKrAW/joYLsmNPsL+77XHzgzXd3W0pyIL8R5tlWw+i0pY7BnlyvmAtIFlAeLtj90Y8TU9WxMx3jjMmu+Mm17t8vGDg+LMDkdrPlAs1UjJPaSs9Z24Sm3xevMzgqYLGyOYNJaOJOV2AgIZMhlM+Wc92wqD4+jrLkcIXMUxGY9f6vR8jYZcuyHE7KAtP2paclRkGgyhkwOxq1F9IVpyo10kXp99NtHPOORg/D4mxdseR02x1vwvM9rAM+oMZu2LCAi67DQAStbqL4dscoFnFlHrQz21/L4ZwNjZ89ABuCzFrZMALyNGchC7GLiRVjeWdjQ5jX4QwKQlc61vGDM8aVIBg0H6cKSE4h1AdpnoBwRGa4t21gJIM6WoxPytRqzmiNxt4exCOtqJDhgRncw7mADZteaO4z1WVaV2ewtnVgYvQmYkDXGAZQ6MM26agtBiR2EfJVZQToKIOX81HNSsWVnGvMA4ABqo6S6uKU2WDK3nmnW02Q68tTEtWQZmnIB5QNFiNzQls8AmkYY4q22z2sA99CYI0KAPJxF2WeOONAG2BCZB3r2eMnCmFmLHMGmoxoYG0PORKC43wGaHHRVYGyb7LEGyjBQxiDLqLfVVpObYjtAOeOKV0vLEpMcwBzsmec1IA+QrqAcWX5rnblWk6Myzlcxh51ZPLd5yYhRAS0dgJU5G8B0YggRHrXZGjg4OM9qTq4TebVCFhA1sDPdkDVmAGCPjffi+ogooXA2VwUrgLrKbfW22poxF1AGMMCY1CItfDudex7dNE0jzfrQOhobQ06m7ABcQXmd5Wf1l2WMYRlXzoiMW2yf1/coYiQ6wuW8c68C1k7dK7AQE7grVGyqUx/acLBgafDjIy2UJ824Sm7jeMoc8asf2UvAwrmx0NbuA8yFIQNIb3BqyoVFL518I32WKnv2H6Rw9oWmnPfduJ4sarezALOYtjx3E9sjXE5Kt4rNDwX2/N2eyxpRgoYx5DFxPoyfgiF4oooP00jVtrC5I89eKtYiEaKhwoEPljk6WZ2Vp70tskm7ktdccZnDgVmEE6SHUxoDuIFzgK6W895lkMqYgewyAoraFsiOLK2NRqrBjokUjzxOeSI79ridLPqCJDXlKIb/hvY09WTTl29KuNzJw+U8koNs3Ddtn9cArupgIuNWV8vtCJ0TMcfZLMDEwCwgdt03OgAzACX3gFPWeFXX1hKMNUCanCEXYC7yRba5SYZQzvkKKQPIObzBmAc7XgKzJmgvIjFK4kgyaNGFpmyyxwBl05hlGSKXf7Ic9w0jKOiuJJTdbrVGQ0qrrK6H1gx1h6Dpp1HD9uRNgiOLTZTAXkQLAGZuaMLoNHpZBouelNBFMDNnYpMBsw6glhFbD2AFzrdcUIBwgLIfomDGDsRLYLbU6eb3H0oY3JEdmD155HELIJbUlCe2FUUF5SPNOFLPOtd8xprL+J9dwsub10T0uQC+EfZJzwC+SlV/4h7PJ1i3+LcD+E0AX6yqP1PufxOAnwfwPlX9itte6/qefy00MfEGkWQRGuKlikjcqWHgzF1y0jkWDxYMDFD2rg5ZdMWZNUg9WgObmUzx6S7iNXMEbrmWNTArBkBvArPpe9zH8aE/11TrC/JF1ZSDKYtNPgpwDn3ZNeeSAXPrdezM4vmNyy2T4oCendcPVkQZHaY3Rzz9Yz6hqeCE5s2ACSwts9tmFczEmLmBuqKzFclqZKUDRAmtt1UGquR2RjHl92dM6NsWURXnRldvewInMI+i9pEpO7URb3xoo+GEseMRDveoWZ3lA3e80m4yzjsLPm04/45Uw+d0Md6b9nLn9Y8A+EFVVSL6fQC+H8DvvcfzPxPAx/rfHwTwnX4b9o0A/pdrXujqDibaxeMwZeBeU+BElgkYiMvmBFSJbg4ENAJ19lAdAk/e9sa9w9bqRgHybgoht2UXBl0CcIA0sGTMZWLq1q9xndUVkKuUoSvA1pA2NPcNpAcQD2ehegicPZadGVdQtn0BnTqQ8kU3UO7dVx/dMipvlTLUnrfbc1kD4ejOhoOKzyHrH9mZcNKeTX1POirPnbThxC2PWRccxqvtaJ1uxErIRosyq/E9ZWJK1gIPMI5tWD0YoNR3UVrIrqobc7uEiHIF5EiFZvHcrdH1u9ZOju0jzxaF4UBsxYiMRb/CN8auIfiwAOYCykeXKz6cn6Z8YbfGsg8QHAg4Em33+7OLe2nzWlV/vex+OAqNI6L/DMB/DOt7+T5V/fqNl/hcAN+jpkH+FBF9BBH9TlX950T0ibBu8X8HwJ0d6e/f808cqdw5ReLhcrPHNjfXJdwJaEzbgvdDtvCw5ZQqwnOnTrw1HHp9CcjRjeFMwtiUMzZA7TYZo24HG67A7E662I6uEjxrYdjVyWcgDY++oGACMmKW0+FXy6puOT4uAPQeiP/sRhvA0EghqhCXM5oqOgSPyWKboxiPEGVoHRfQPrlLq4Mxk2SoaDRlmLzpsBTgjv0uvNgHkIANDJAGbpeYgQHGAcCxne2eWBb7UwHmYMQMXQIzTL4IhlzZcYPgMd34KsNilqOF1ELG2PhOto3Sny9zXhPR5wH4CwD+NQCf5cc+A8aCPwk21D9IRJ+iqj++evpbAPyzsv/LAN5CRP83gP8awBcB+LRrzuPuWhm+7M5CyypZDQ6AAQ0cNF3WUH+8FWMx5mthF2odHZSTBWsjsGt3GfyeQIwCxM6qsSFfJDDTYn/7gvwhGS8c+1iAdAXilDB0yY5RteRk086KFQbILleYauWasjhAy9geLDmcf3Zcb4tjnnfG/DzGYDCR68UeXxseYAVOKWUIDlY+DkwCgTX0zfAvBg7acSKXNpTQ2Rj0RB2ztnQWWs/Jc2CevYjMGpwrSMf9l6/HdeUVW64gDGABxJUR1+2JolN4ZEFq6shMgsde+yKaIT/mkztQZ49ftu0jPCMwATpkI7ocx/yw8/rNRPSBsv8eVX3PpQer6vsAvI+IPgUmPXw6gM/wv5/1h70BBtRrYL4UcvDlAP62qv4zuuT0XNn1URmq0N6tZkTzaMRI02aLyIAotLHl/Ts4W91kgk6c4E4Th1fC45Op6Mh+zGWQBfiu9eYyFAbKlRJfGKJiZ+C80pUraAezze14XIJyMOjhICXfziSSYMk1u08EOs++3R2g7XG3RmUAd1eg2+1Oa7BuNiDgERQny+hIELuhDlFGY8WNGmc+Us/u6h2MJ3KweGclPNGDdWKHdcaJTuzBmuNx0QgYQPagNP2Zx3GfxH0lXYieA9q6a0+w07iOqPRW2z3F8axpgcGWG9W0anUwPuXjDIhlALOz5GMwZhIcIXjkkRgHKA4EHFzGuChlAA89rz+kqhelAyJ6J4Av8923q+qvAoCq/jgR/ZtE9GYYmvwFVf0rtz0XxpA/ujzkrQB+FcAnA/j3iOjLYaB+JKJfV9X//NJ53a/sJ/OQMtABtKEHeVAxqWnFamcOClnDj8MZSoAyMXnNAN+P9utEmS04HH4GvsGg423tMWvU3biGNc5VtgwMhgwkCA/QVoxWWgOMK2BXMF4AsuvF5MBcjw2AjtuQNO6YnKrAPN/+mN2uMvbojAb1tk8KEHDwIhedFB0mXwQodiI0nYw5s/Wr7B51FJ3XE5gdpDuNLu0nHd3gA6zrfmzb7RKIt1jzuodeBV4AyXrjWN0PcA1ADjCu9S6YFI/pJkt6xnOCHQcY533QTL8+wEPk4BXlbmONL3heq+q7AbwbAIjo3yIicuffHwBwBPAvAbwfwDcS0X+rqr9ORG8BcKrP9ef/IICvIKL3wpx+/5+q/nMA7yiP+WIAb7sNlIFrw+XEQttsziqoOyijhzw8XMLNgJjmIW8YQDcH4wHe8RxqNNgz+2Mqo47X5nO54ir54uK1+VPPZA3N+gQVhJMl1+JDrodl2FvXsV2AOGWLAOVgxV0GKPeeLDmLGF2ql6HYpYznNEYspwUC9fA5jYmesoaVppwNZP1YB+GADiHGEz2gKztYe+SFg+uNTgugTpCOWwfdAOzYBjD2N5x9d15bAHJew9hfbNNSC65AfHQph70mdYBwZcgtQXxOMD/kMUlQPmCExzUf97OSn8DLntf/EYA/SUQnAK8C+Hx35P0QEX0cgJ90KeLXAXwhgA+unv+3Ycz5F2Dhcv/Js57IPTqYtARZZTZwdpasHQC7OyVAt6mlMzsAqx8nL4UIwFkywbtR5j55NblFmnUAdpxTMuUVc76HVaZst/U+TUfcAGbXlTMGEEsgBhap1bROs459sTDDBF5nzVr05zHuF2yPyngQa2QRD6E1x+rOlDmXNdSkjgZjznHbycDWHIIGtievOtcdfA+5T6lNB0MezJgXiS3JyrHNmK+xs7oUReoIALbt2hBAEoQBLIC47odkMarHabJkCzmMYvkmXzQf59CWY9w37eVGZXwLgG+5cN93wGKUb3u+AnjnHY/5mwD+5l3ncp3GLApChxIn+KZUYbPaHte8mAWxAfWCFUeFfIp4Hnu+HwvZoj4nHxtWj6Ow5WrXzOENvDvrtVf1XdHlcd+mreP1WMo8AboOxHXyiaambPuSTsCLjj970Z0xP4AxGAeCg7KB8YEIPZJN3CN8UoGQgWVXa9AawHkqQBoguwW6pjWPx/YCuHl8BcrVrgHorQ7xrTBku+YC0jTqVwQDjudU8K6sGoDVC4nHYdRcZgCHWDxjaMp2boQDmoP0pWvZ5zVw7yJGFgSP7qFxzA4qzoDjcQ7MQ2aw2p60AFwe25UJM50BcI2c1wLUFz2cd2lYdxxftLbZAuz1/QVUY/tMiljta5flc0SGs+82ppzvCxv73R7EjDkTmoNyaHSjpoN9NuJgHZ+QSRoy5AcqXdhB6GUuLkC6ADfKa+X2BRCWW5gHbzEODEAGlr32KvDGfWtmPRh32S5AbPte19qfd/RrtoxK05QDnC8y5bB9XgO4SmNWAOL1LwBCB7yMYQJ0gFE4COEgvGC7XJx5BbCBfBxtAHECuB9ffKxXhp5cZWvAXjBmOTueICorMK3HK8DW2her7D5dPUfra972Q7I7/57LGnmlLIeUAwHd57uxOotrBpBNdQDzrAzoUAhpuBrsMYTCkl2O8MiPpWRByWbWDHmLMT/zda4cg0tw1rPjQ4PW1XF/HC0XpiFVxGMqQ85t8IIpM2gzhnmf12ZXZ/7FtDTJrCeoKiyyYmG08hGv7j9jujxea/XA28+L717aXW3XREFc8byzELe1HLFiwwu5Ys2UbwuXU4We9gn8vDbAGWA0TGSsOExW0kAP52Dcv2KpfTWLJT/D5e36cWt7yICxu74ll0LXeAWc68et5Yg1G673rx19m6AM7PPa7Xop44xRFs4QqXxpty9F9NKH8gx29qPwDHa7lnvvF7vHY5/9fRX7ku+hbA0SbbXXy2c63THdZA24zzE916D/LHZZy32W17r+Yi4C7x22z2szui2JgYj+BYBfenGns9st9m+o6kfFDhH9HQBvfqDX/pCq/rEHeq3fErbP7deV5dx+4HkN/Bad27cC82677bbbbi/eHlCk3W233Xbb7SFsB+bddtttt9eZ7cC822677fY6sx2Yd9ttt91eZ7YD82677bbb68x2YN5tt912e53ZDsy77bbbbq8z24F5t9122+11Zjsw77bbbru9zuz/B+NHUBsoSySWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sympy as sp\n",
    "import pylab as plt\n",
    "%matplotlib inline\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import pylbm\n",
    "\n",
    "u, X, Y = sp.symbols('u, X, Y')\n",
    "\n",
    "def solution(x, y, t, k, l):\n",
    "    return np.sin(k*np.pi*x)*np.sin(l*np.pi*y)*np.exp(-(k**2+l**2)*np.pi**2*mu*t)\n",
    "\n",
    "def plot(i, j, z, title):\n",
    "    im = axarr[i,j].imshow(z)\n",
    "    divider = make_axes_locatable(axarr[i, j])\n",
    "    cax = divider.append_axes(\"right\", size=\"20%\", pad=0.05)\n",
    "    cbar = plt.colorbar(im, cax=cax, format='%6.0e')\n",
    "    axarr[i, j].xaxis.set_visible(False)\n",
    "    axarr[i, j].yaxis.set_visible(False)\n",
    "    axarr[i, j].set_title(title)\n",
    "\n",
    "# parameters\n",
    "xmin, xmax, ymin, ymax = 0., 1., 0., 1.\n",
    "N = 128\n",
    "mu = 1.\n",
    "Tf = .1\n",
    "dx = (xmax-xmin)/N # spatial step\n",
    "la = 1./dx\n",
    "s1 = 2./(1+4*mu)\n",
    "s2 = 1.\n",
    "k, l = 1, 1 # number of the wave\n",
    "\n",
    "dico = {\n",
    "    'box': {'x':[xmin, xmax], \n",
    "            'y':[ymin, ymax], \n",
    "            'label': 0},\n",
    "    'space_step': dx,\n",
    "    'scheme_velocity': la,\n",
    "    'schemes':[\n",
    "        {\n",
    "            'velocities': list(range(5)),\n",
    "            'conserved_moments': u,\n",
    "            'polynomials': [1, X/LA, Y/LA, (X**2+Y**2)/(2*LA**2), (X**2-Y**2)/(2*LA**2)],\n",
    "            'equilibrium': [u, 0., 0., .5*u, 0.],\n",
    "            'relaxation_parameters': [0., s1, s1, s2, s2],\n",
    "        }\n",
    "    ],\n",
    "    'init': {u: (solution, (0., k, l))},\n",
    "    'boundary_conditions': {\n",
    "        0: {'method': {0: pylbm.bc.AntiBounceBack,}},\n",
    "    },\n",
    "    'generator': 'cython',\n",
    "    'parameters': {LA: la},    \n",
    "}\n",
    "\n",
    "sol = pylbm.Simulation(dico)\n",
    "x = sol.domain.x\n",
    "y = sol.domain.y\n",
    "\n",
    "f, axarr = plt.subplots(2, 2)\n",
    "f.suptitle('Heat equation', fontsize=20)\n",
    "\n",
    "plot(0, 0, sol.m[u].copy(), 'initial')\n",
    "\n",
    "while sol.t < Tf:\n",
    "    sol.one_time_step()\n",
    "\n",
    "sol.f2m()\n",
    "z = sol.m[u]\n",
    "ze = solution(x[:, np.newaxis], y[np.newaxis, :], sol.t, k, l)\n",
    "plot(1, 0, z, 'final')\n",
    "plot(0, 1, ze, 'exact')\n",
    "plot(1, 1, z-ze, 'error')\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}