{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The wave equation in 1D\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{\\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{1}{3}$ on the wave equation.\n", "\n", "The problem reads\n", "$$ \\drondtt\\rho = c^2 \\drondxx\\rho, \\qquad t>0, \\quad x\\in(0,2\\pi),$$\n", "\n", "where $c$ is a constant scalar. In this session, two different kinds of boundary conditions will be considered:\n", "\n", "* periodic conditions $\\rho(0)=\\rho(2\\pi)$,\n", "\n", "* Homogeneous Dirichlet conditions $\\rho(0)=\\rho(2\\pi)=0$.\n", "\n", "The problem is transformed into a one order system:\n", "\n", "$$\n", "\\begin{aligned} &\\drondt \\rho + \\drondx q = 0, && t>0, \\quad x\\in(0, 2\\pi),\\\\ &\\drondt q + c^2 \\drondx \\rho = 0, && t>0, \\quad x\\in(0, 2\\pi).\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The scheme $\\DdQq{1}{3}$\n", "\n", "The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximation of the solution on discret points of $(0,2\\pi)$ at discret instants.\n", "\n", "The spatial mesh is defined by using a numpy array. To simplify, the mesh is supposed to be uniform.\n", "\n", "First, we import the package numpy and we create the spatial mesh. One phantom cell has to be added at each bound for the treatment of the boundary conditions." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADnhJREFUeJzt23+sX3ddx/HnixYmP3Q/S1vb1bu4RlM0YeSbEQKSBfaji0AX3R+bURuD6T/MgIvRIomDQSIYZcSIJM2GqYiMZbjQiFLLNuKP4Oi3YwbKGK1zpK0bK3SilegyePvHPdP7uX57f/T77T33Xp6P5OZ+zzmf9ry/y9Ln/Z5zbqoKSZKe94K+B5AkLS+GQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGmv7HuBsXHLJJTU1NdX3GJK0ohw6dOhbVbVuvnUrMgxTU1MMh8O+x5CkFSXJNxayzktJkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1JhIGJJsT/JYkqNJdo84fl6ST3bHH0oyNev4liSnk/zGJOaRJJ29scOQZA3wYeB6YBtwc5Jts5a9FXimqi4H7gA+MOv4B4G/HncWSdL4JvGJ4UrgaFU9XlXPAncDO2at2QHs7V7fC7wxSQCS3AD8C3B4ArNIksY0iTBsAo7N2D7e7Ru5pqqeA74DXJzkZcBvAe+ZwBySpAno++bzu4E7qur0fAuT7EoyTDI8efLkuZ9Mkn5ArZ3A33ECuHTG9uZu36g1x5OsBc4Hvg28Grgxye8BFwDfT/JfVfVHs09SVXuAPQCDwaAmMLckaYRJhOEgsDXJZUwH4CbgF2at2QfsBL4A3Ag8UFUF/MzzC5K8Gzg9KgqSpKUzdhiq6rkktwD7gTXAR6vqcJLbgWFV7QPuAj6W5Chwiul4SJKWoUz/4L6yDAaDGg6HfY8hSStKkkNVNZhvXd83nyVJy4xhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGhMJQ5LtSR5LcjTJ7hHHz0vyye74Q0mmuv3XJDmU5Mvd9zdMYh5J0tkbOwxJ1gAfBq4HtgE3J9k2a9lbgWeq6nLgDuAD3f5vAW+uqp8GdgIfG3ceSdJ4JvGJ4UrgaFU9XlXPAncDO2at2QHs7V7fC7wxSarqS1X1r93+w8CLk5w3gZkkSWdpEmHYBBybsX282zdyTVU9B3wHuHjWmp8HHq6q/57ATJKks7S27wEAkryC6ctL186xZhewC2DLli1LNJkk/eCZxCeGE8ClM7Y3d/tGrkmyFjgf+Ha3vRm4D/jlqvrnM52kqvZU1aCqBuvWrZvA2JKkUSYRhoPA1iSXJXkRcBOwb9aafUzfXAa4EXigqirJBcBngN1V9Q8TmEWSNKaxw9DdM7gF2A88CtxTVYeT3J7kLd2yu4CLkxwFbgWef6T1FuBy4HeSPNJ9vXzcmSRJZy9V1fcMizYYDGo4HPY9hiStKEkOVdVgvnX+5rMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktSYSBiSbE/yWJKjSXaPOH5ekk92xx9KMjXj2Du7/Y8luW4S80iSzt7YYUiyBvgwcD2wDbg5ybZZy94KPFNVlwN3AB/o/uw24CbgFcB24I+7v2+iNmzYQJL/97Vhw4ZJn2rJz7ea39tSn281v7elPt9qfm9Lfb6lfm8Aqarx/oLkNcC7q+q6bvudAFX1uzPW7O/WfCHJWuApYB2we+bamevmOudgMKjhcLiYGc94bNz33/f5VvN7W+rzreb3ttTnW83vbanPN8lzJTlUVYP51k3iUtIm4NiM7ePdvpFrquo54DvAxQv8s5KkJbRibj4n2ZVkmGR48uTJvseRpFVrEmE4AVw6Y3tzt2/kmu5S0vnAtxf4ZwGoqj1VNaiqwbp16yYwtiRplEmE4SCwNcllSV7E9M3kfbPW7AN2dq9vBB6o6Ytj+4CbuqeWLgO2Al+cwEySpLM0dhi6ewa3APuBR4F7qupwktuTvKVbdhdwcZKjwK38303nw8A9wFeBzwJvq6rvjTvTbOvXr1/U/pV0vtX83pb6fKv5vS31+Vbze1vq8y31e4MJPJXUh8U+lSRJWtqnkiRJq4hhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGmOFIclFSQ4kOdJ9v/AM63Z2a44k2dnte0mSzyT5WpLDSd4/ziySpMkY9xPDbuD+qtoK3N9tN5JcBNwGvBq4ErhtRkB+v6p+ErgCeG2S68ecR5I0pnHDsAPY273eC9wwYs11wIGqOlVVzwAHgO1V9d2qehCgqp4FHgY2jzmPJGlM44ZhfVU92b1+Clg/Ys0m4NiM7ePdvv+V5ALgzUx/6pAk9WjtfAuSfA7YMOLQu2ZuVFUlqcUOkGQt8AngD6vq8TnW7QJ2AWzZsmWxp5EkLdC8Yaiqq890LMk3k2ysqieTbASeHrHsBHDVjO3NwOdnbO8BjlTVh+aZY0+3lsFgsOgASZIWZtxLSfuAnd3rncCnR6zZD1yb5MLupvO13T6SvA84H3jHmHNIkiZk3DC8H7gmyRHg6m6bJIMkdwJU1SngvcDB7uv2qjqVZDPTl6O2AQ8neSTJr445jyRpTKlaeVdlBoNBDYfDvseQpBUlyaGqGsy3zt98liQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJahgGSVLDMEiSGmOFIclFSQ4kOdJ9v/AM63Z2a44k2Tni+L4kXxlnFknSZIz7iWE3cH9VbQXu77YbSS4CbgNeDVwJ3DYzIEl+Djg95hySpAkZNww7gL3d673ADSPWXAccqKpTVfUMcADYDpDkZcCtwPvGnEOSNCHjhmF9VT3ZvX4KWD9izSbg2Izt490+gPcCfwB8d8w5JEkTsna+BUk+B2wYcehdMzeqqpLUQk+c5JXAj1fVryeZWsD6XcAugC1btiz0NJKkRZo3DFV19ZmOJflmko1V9WSSjcDTI5adAK6asb0Z+DzwGmCQ5Ilujpcn+XxVXcUIVbUH2AMwGAwWHCBJ0uKMeylpH/D8U0Y7gU+PWLMfuDbJhd1N52uB/VX1kar60aqaAl4HfP1MUZAkLZ1xw/B+4JokR4Cru22SDJLcCVBVp5i+l3Cw+7q92ydJWoZStfKuygwGgxoOh32PIUkrSpJDVTWYb52/+SxJahgGSVLDMEiSGoZBktQwDJKkhmGQJDUMgySpYRgkSQ3DIElqGAZJUsMwSJIahkGS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhqGQZLUMAySpIZhkCQ1DIMkqWEYJEkNwyBJaqSq+p5h0ZKcBL6xxKe9BPjWEp9zHM57bjnvueW858aPVdW6+RatyDD0IcmwqgZ9z7FQzntuOe+55bz98lKSJKlhGCRJDcOwcHv6HmCRnPfcct5zy3l75D0GSVLDTwySpIZhmEeS7UkeS3I0ye6+55lPko8meTrJV/qeZT5JLk3yYJKvJjmc5O19zzSXJD+U5ItJ/qmb9z19z7QQSdYk+VKSv+x7loVI8kSSLyd5JMmw73nmk+SCJPcm+VqSR5O8pu+ZxuWlpDkkWQN8HbgGOA4cBG6uqq/2OtgckrweOA38aVX9VN/zzCXJRmBjVT2c5IeBQ8ANy/W/b5IAL62q00leCPw98Paq+seeR5tTkluBAfAjVfWmvueZT5IngEFVrYTfCyDJXuDvqurOJC8CXlJV/9b3XOPwE8PcrgSOVtXjVfUscDewo+eZ5lRVfwuc6nuOhaiqJ6vq4e71fwCPApv6nerMatrpbvOF3dey/skqyWbgZ4E7+55lNUpyPvB64C6Aqnp2pUcBDMN8NgHHZmwfZxn/w7WSJZkCrgAe6neSuXWXZR4BngYOVNWynhf4EPCbwPf7HmQRCvibJIeS7Op7mHlcBpwE/qS7XHdnkpf2PdS4DIN6l+RlwKeAd1TVv/c9z1yq6ntV9UpgM3BlkmV7uS7Jm4Cnq+pQ37Ms0uuq6lXA9cDbusujy9Va4FXAR6rqCuA/gWV/L3I+hmFuJ4BLZ2xv7vZpQrpr9Z8CPl5Vf9H3PAvVXS54ENje9yxzeC3wlu6a/d3AG5L8Wb8jza+qTnTfnwbuY/qS7nJ1HDg+45PjvUyHYkUzDHM7CGxNcll3U+kmYF/PM60a3c3cu4BHq+qDfc8znyTrklzQvX4x0w8lfK3fqc6sqt5ZVZuraorp/3cfqKpf7HmsOSV5afcgAt0lmWuBZfuEXVU9BRxL8hPdrjcCy/LhicVY2/cAy1lVPZfkFmA/sAb4aFUd7nmsOSX5BHAVcEmS48BtVXVXv1Od0WuBXwK+3F23B/jtqvqrHmeay0Zgb/e02guAe6pqRTwCuoKsB+6b/pmBtcCfV9Vn+x1pXr8GfLz74fFx4Fd6nmdsPq4qSWp4KUmS1DAMkqSGYZAkNQyDJKlhGCRJDcMgSWoYBklSwzBIkhr/A57/wrLrWAp+AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import pylab as plt\n", "\n", "def mesh(N):\n", " xmin, xmax = 0., 2.*np.pi\n", " dx = (xmax-xmin)/N\n", " x = np.linspace(xmin-.5*dx, xmax+.5*dx, N+2)\n", " return x\n", "\n", "x = mesh(10)\n", "plt.plot(x, 0.*x, 'sk')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To simulate this system of equations, we use the $\\DdQq{1}{3}$ scheme given by\n", "\n", "* three velocities $v_0=0$, $v_1=1$, and $v_2=-1$, with associated distribution functions $\\fk{0}$, $\\fk{1}$, and $\\fk{2}$,\n", "\n", "* a space step $\\dx$ and a time step $\\dt$, the ration $\\lambda=\\dx/\\dt$ is called the scheme velocity,\n", "\n", "* three moments\n", " $$ \\mk{0}=\\sum_{i=0}^{2} \\fk{i}, \\quad \\mk{1}= \\lambda \\sum_{i=0}^{2} v_i \\fk{i}, \\quad \\mk{2}= \\frac{\\lambda^2}{2} \\sum_{i=0}^{2} v_i^2 \\fk{i},$$\n", " \n", "and their equilibrium values $\\mke{0} = \\mk{0}$, $\\mke{1} = \\mk{1}$, and $\\mke{2} = c^2/2 \\mk{0}$.\n", "\n", "* a relaxation parameter $s$ lying in $[0,2]$.\n", "\n", "In order to prepare the formalism of the package pylbm, we introduce the three polynomials that define the moments: $P_0 = 1$, $P_1=\\lambda X$, and $P_2=\\lambda^2/2 X^2$, such that\n", "$$ \\mk{k} = \\sum_{i=0}^2 P_k(v_i) \\fk{i}.$$\n", "\n", "The transformation $(\\fk{0}, \\fk{1}, \\fk{2})\\mapsto(\\mk{0},\\mk{1}, \\mk{2})$ is invertible if, and only if, the polynomials $(P_0,P_1,P_2)$ is a free set over the stencil of velocities.\n", "\n", "The lattice Boltzmann method consists to compute the distribution functions $\\fk{0}$, $\\fk{1}$, and $\\fk{2}$ 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", " $$\\mks{2}(t,x)=(1-s)\\mk{2}(t,x)+s\\mke{2}(t,x).$$\n", "* m2f: \n", "$$\n", " \\begin{aligned}\\fks{0}(t,x)&\\;=\\mk{0}(t,x)-2\\mks{2}(t,x)/\\lambda^2, \\\\ \\fks{1}(t,x)&\\;=\\mk{1}(t,x)/(2\\lambda)+\\mks{2}(t,x)/\\lambda^2, \\\\ \\fks{2}(t,x)&\\;=-\\mk{1}(t,x)/(2\\lambda)+\\mks{2}(t,x)/\\lambda^2.\\end{aligned}\n", " $$\n", "* transport: \n", "$$\n", " \\begin{aligned} \\fk{0}(t+\\dt, x)&\\;=\\fks{0}(t,x), \\\\ \\fk{1}(t+\\dt, x)&\\;=\\fks{1}(t,x-\\dx), \\\\ \\fk{2}(t+\\dt, x)&\\;=\\fks{2}(t,x+\\dx). \\end{aligned}\n", " $$\n", "* f2m:\n", "$$\n", " \\begin{aligned}\\mk{0}(t+\\dt,x)&\\;=\\fk{0}(t+\\dt,x)+\\fk{1}(t+\\dt,x)+\\fk{2}(t+\\dt,x), \\\\ \\mk{1}(t+\\dt,x)&\\;=\\lambda\\fk{1}(t+\\dt,x)-\\lambda\\fk{2}(t+\\dt,x), \\\\ \\mk{2}(t+\\dt,x)&\\;=\\tfrac{1}{2}\\lambda^2\\fk{1}(t+\\dt,x)+\\tfrac{1}{2}\\lambda^2\\fk{2}(t+\\dt,x).\\end{aligned}\n", "$$\n", "\n", "The moments of order $0$, $\\mk{0}$, and of order $1$, $\\mk{1}$, being conserved during the relaxation phase, the equivalent equations of this scheme read at first order\n", "\n", "$$\n", "\\begin{aligned}&\\drondt\\mk{0} + \\drondx\\mk{1} = \\grandO(\\dt),\\\\ &\\drondt\\mk{1} + 2\\drondx\\mke{2} = \\grandO(\\dt). \\end{aligned}\n", "$$\n", "\n", "We implement a function equilibrium that computes the equilibrium value $\\mke{2}$, the moment of order $0$, $\\mk{0}$, and the velocity $c$ being given in argument." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def equilibrium(m0, c):\n", " return .5*c**2*m0\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create three vectors $\\mk{0}$, $\\mk{1}$, and $\\mk{2}$ with shape the shape of the mesh and initialize them. The moments of order $0$ and $1$ should contain the initial value of the unknowns $\\rho$ and $q$, and the moment of order $2$ the corresponding equilibrium value.\n", "\n", "We create also three vectors $\\fk{0}$, $\\fk{1}$ and $\\fk{2}$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def initialize(mesh, c, la):\n", " m0 = np.sin(mesh)\n", " m1 = np.zeros(mesh.shape)\n", " m2 = equilibrium(m0, c)\n", " f0 = np.empty(m0.shape)\n", " f1 = np.empty(m0.shape)\n", " f2 = np.empty(m0.shape)\n", " return f0, f1, f2, m0, m1, m2\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Periodic boundary conditions\n", "\n", "We implement the four elementary functions f2m, relaxation, m2f, and transport. In the transport function, the boundary conditions should be implemented: we will use periodic conditions by copying the informations in the phantom cells." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def f2m(f0, f1, f2, m0, m1, m2, la):\n", " m0[:] = f0 + f1 + f2\n", " m1[:] = la * (f2 - f1)\n", " m2[:] = .5* la**2 * (f1 + f2)\n", "\n", "def m2f(f0, f1, f2, m0, m1, m2, la):\n", " f0[:] = m0 - 2./la**2 * m2\n", " f1[:] = -.5/la * m1 + 1/la**2 * m2\n", " f2[:] = .5/la * m1 + 1/la**2 * m2\n", "\n", "def relaxation(m0, m1, m2, c, s):\n", " m2[:] *= (1-s)\n", " m2[:] += s*equilibrium(m0, c)\n", "\n", "def transport(f0, f1, f2):\n", " # periodic boundary conditions\n", " f1[-1] = f1[1]\n", " f2[0] = f2[-2]\n", " # transport\n", " f1[1:-1] = f1[2:]\n", " f2[1:-1] = f2[:-2]\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We compute and we plot the numerical solution at time $T_f=2\\pi$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVdXV//HPYgAHFQRhokiPYKMICvoYE4OKgk9QNNFgi2g0RBNj7EAUQVSENLArxYgVsT0SfxpFBTH2ITQBEYIaBlAQu9JZvz/2Id7BGabccm75vl+v+5p7zz1l3VFm3bX3PnubuyMiIrJNnbgDEBGR7KLEICIi5SgxiIhIOUoMIiJSjhKDiIiUo8QgIiLlKDFIzjOzZ8xswA7ev9PMhlbzXDPM7LzURSeSe5QYJCuZ2ftm1qs6+7r7ce4+KTrubDP753bvn+/u16UjzlygZCc1pcQgIiLlKDFI1ttWBZjZn83sUzN7z8yOS3h/hpmdZ2b7A3cCh5nZV2b2WfT+PWZ2ffS8iZk9ZWZronM9ZWYtqxnHcDN7xMzuN7MvzWy+me1jZkPMbLWZLTezYxP238vMpprZJ2a21Mx+lcS5djOziWa2ysxWmNn1ZlZU1e/HzG4AfgTcGv1ObrVgTHSdL6Jrd0rmv5HkFyUGyRWHAouBZsAfgYlmZok7uPsi4HzgNXff1d0bV3CeOsDfgDZAa2AdcGsN4jgeuA9oAswGno3O2QIYAdyVsO9koAzYCzgZGGlmR9XyXPcAm4H2QDfgWCCxeajC34+7XwW8DFwY/U4ujI49AtgH2A34ObC2Br8DyXNKDJIrPnD38e6+BZgENAf2qOlJ3H2tuz/m7t+4+5fADcCPa3CKl939WXffDDwClACj3H0TIRG0NbPGZtYKOBwY5O7r3X0OMAE4qxbn2gP4X+Bid//a3VcDY4BTa/n72QQ0BPYDzN0XufuqGvwOJM/VjTsAkWr6cNsTd/8mKhZ2relJzGxnwh/VPoRv6gANzawo+qNalY8Snq8DPk44bl30c1dClfBJlHy2+QDoXstz1QNWJRRJdYDlCcdX+/fj7i+a2a3AbUAbM3scuNzdv6jwE0vBUcUg+aaq6YIvA/YFDnX3RoQmFQCr/JBaWQnsbmYNE7a1BlbU4lzLgQ1AM3dvHD0auXvHah7/nd+Ju9/s7gcDBxCalK6oRVySp5QYJN98BLQ0s/qVvN+Q8G38MzPbHRiWjiDcfTnwKnCjmRWbWRfgXOD+WpxrFfAc8Bcza2RmdcxsbzOrbhPYR8D3t70wsx5mdqiZ1QO+BtYDW2sal+QvJQbJNy8CC4APzezjCt4fCzQAPgZeB/6RxlhOA9oSqocngGHu/nwtz3UWUB9YCHwKPEroR6iOm4CToxFLNwONgPHReT4gdDz/qZZxSR4yLdQjIiKJVDGIiEg5SgwiIlKOEoOIiJSjxCAiIuXk5A1uzZo187Zt28YdhohITpk1a9bH7l5S1X45mRjatm1LaWlp3GGIiOQUM/ugOvupKUlERMpRYhARkXKUGEREpJyc7GMQEcmkTZs2UVZWxvr16+MOpVqKi4tp2bIl9erVq9XxSgwiIlUoKyujYcOGtG3blu3Wh8o67s7atWspKyujXbt2tTpHSpqSzOzuaJnAtyt538zs5mh5w3lmdlDCewPMbEn0GJCKeEREUmn9+vU0bdo065MCgJnRtGnTpKqbVPUx3ENY+KQyxwEdosdA4A6AhGmPDwUOAYaZWZPKTiIiEpdcSArbJBtrSpqS3H2mmbXdwS79gHs9TOX6erRcYXOgJzDN3T8BMLNphATzUCrikhrYsAEWLIB//xtWrYK1a2HbzLu77w7Nm0O7dtClCxQXxxuriKRVpvoYWlB+GcKyaFtl27/DzAYSqg1at26dnigLydat8Oab8NRT8OyzMHcubNpUfh+zb5PDNkVF0KkTHHssHH88HHYY1FVXlUg+yZl/0e4+DhgH0L17dy0iUVtlZXD33eHxwQfhD/0PfgCXXgoHHwz77gt77RWqhDp1QmL49FNYuRLefRdmzYLXX4exY+FPfwr7nnMOnHtuqChEJOdl6j6GFUCrhNcto22VbZdU+/e/4bzzwh/vYcOgQwe4/35YswZmzoRRo+CUU0JTUbNmISlAqBp23z1UCT/9KdxwA7zwAnz8MUyZAt26wY03Qvv2cOaZsGhRvJ9TJM/dfvvtdOrUiTZt2nDLLbek5RqZqhimAhea2WRCR/Pn7r7KzJ4FRiZ0OB8LDMlQTIXhk0/gmmvgzjtDk8/558Mll8D3v1/1sTvSqFFIJKecAsuXw803w+23w4MPwtlnh2Sxxx4p+QgiWeXii2HOnNSes2vXUIVX4bHHHmPatGnMnj2bjz/+mM6dO3PBBRdQN8XNuakarvoQ8Bqwr5mVmdm5Zna+mZ0f7fI0sAxYSlhr9jcAUafzdcBb0WPEto5oSZI7TJwI++wDd9wBAwfCsmVwyy3JJ4XttWoVmpXefz80Sd1/f7juzTeHvgwRSYmbb76Z0aNHU69ePZo3b069evXYmo5/Y+6ec4+DDz7YZQdWrHA/7jh3cD/iCPe5czN7/Xfece/d+9vrL1uW2euLpNjChQvjDsE3btzou+22239fr1y50jt27Fjp/hXFDJR6Nf7Gaq6kfPP3v4f+gBkzQnUwfXroN8ikffeFZ56Bv/0tlNxdusBDGoEskoyFCxfyxRdfsGzZMrZu3cqQIUO46KKL0nItJYZ8sWULXHUVnHACtG0b/iBfeOG3nciZZhb6GubPhwMPhNNPh4sugo0b44lHJMfNnj2bM844g9NOO40uXbrQunVrBg4cmJZr5cxwVdmBL7+E/v3Dt/Rzz4Vbb82em9Batw5Vy5VXhs61OXPgiSegadO4IxPJKXPmzKFv3770798/7ddSxZDrVqyAH/0InnsudDJPmJA9SWGbevVgzBh44AF4441w38S//x13VCI5Zc6cOXTt2jUj11LFkMsWL4ZeveDzz8MdzH12NF1VFjj99FBB9OsH//M/IZl16xZ3VCI5YcaMGRm7liqGXDVvHhxxRGiznzkz+5PCNj/8Ibz2Guy8Mxx5ZLiLWkSyihJDLpo1C3r2DE00M2eGm2NyyT77hLibNQsVz8yZcUckIgmUGHLN/PlhArtGjeDll8PQ0FzUpk2Iv1Ur+MlPQt+DiGQFJYZcsq1PobgYXnwx9yeta948zLu0xx6hKWz27LgjEhGUGHLHihVwzDHh+QsvpH5ai7jstVf4PA0bQu/esHRp3BGJFDwlhlzw+edw3HHw2Wfwj3/AfvvFHVFqtWkDzz8f5lXq0wdWr447IpGCpsSQ7TZuDNNdL1oEjz2Wv8M799knTOexYgX07Qtffx13RCIFS4khm7nDb34T+hMmTvy2KSlfHXYYTJ4MpaUwYIBmZhWJiRJDNrvpppAQrroKzjor7mgyo1+/MIX3Y4/BtdfGHY1IVlmyZAk9e/akU6dOXHbZZey9995puY7ufM5Wzz4Ll10GJ50EI0bEHU1mXXopLFgQPnfHjvDzn8cdkch/xbVOz5YtWzjrrLO47bbbOOigg/jd735Hx44dUxtIRIkhG733Hpx2Wpg++95745shNS5mYd6nxYvDpIAHHpi792uIpMj//d//ccABB3DQQQcBsP/++9O4ceO0XCslicHM+gA3AUXABHcftd37Y4Ajo5c7A99z98bRe1uA+dF7/3H3E1IRU85avx5OPjm0rz/+OOy6a9wRxWOnneDhh0Nn+ymnhKkzdt457qhEqrMCZ1rMnj273CR6c+fOpVevXmm5VtJfRc2sCLgNOA44ADjNzA5I3MfdL3H3ru7eFbgFeDzh7XXb3iv4pABhDYV//Qvuuw/S1H6YM1q2DMuEvv02/O53cUcjEqumTZvyzjvvAPDGG29w7733cuCBB6blWqloozgEWOruy9x9IzAZ6LeD/U8DtJxXRSZPDp3NQ4bA8cfHHU126N07dL7ffTdMmhR3NCKx+cUvfkFpaSmdO3fm8ccfp2nTprRv3z4t10pFU1ILYHnC6zLg0Ip2NLM2QDvgxYTNxWZWCmwGRrn7/1Vy7EBgIEDr1q1TEHaWef99OP/8MGSz0DqbqzJ8OLzyClxwAXTvHjqkRQpMs2bNeCOaU2z58uXMmDGDOmnqf8x0r+apwKPuviVhWxt37w6cDow1swrbT9x9nLt3d/fuJSUlmYg1czZvhjPPDP0KDzwAdTUmoJyiInjwwTBx4CmnwFdfxR2RSKzmzp1LlzSu5Z6KxLACaJXwumW0rSKnsl0zkruviH4uA2YAeXpr7w6MHBm+Ed9xR+5PjJcue+4JDz0URipdcknc0YjEqm/fvowfPz5t509FYngL6GBm7cysPuGP/9TtdzKz/YAmwGsJ25qY2U7R82bA4cDCFMSUO159NdzIdeaZcMYZcUeT3Y48MqwdPWECPP103NGI5K2kE4O7bwYuBJ4FFgFT3H2BmY0ws8RRRqcCk93dE7btD5Sa2VxgOqGPoXASw+efh2TQpg3cdlvc0eSG4cOhc2c47zz45JO4oxHJSylpzHb3p4Gnt9t2zXavh1dw3KtA51TEkJMuvBCWLw8L1jRqFHc0uWGnncJNfz16hN/fgw/GHZEUCHfHzOIOo1rKf/+uuQK7pTaLPPVUGKN/1VVhJJJUX9euMGxY6HN45JG4o5ECUFxczNq1a5P+g5sJ7s7atWspLi6u9TksFz7o9rp37+6lpaVxh1F7n38ehlw2aRLWb65fP+6Ics/mzfCDH8CyZeEGuD33jDsiyWObNm2irKyM9evXxx1KtRQXF9OyZUvq1atXbruZzYpGge6QxkXG4YorYNUqeOIJJYXaqls3NCl16wYDB8KTT4Y5lkTSoF69erQroBGDakrKtBdfhPHjwwyiPXrEHU1u228/uPHGsMCP+hpEUkZNSZn09ddhRE1REcydq0nhUmHLFjj88NCk9M47sPvucUckkrWq25SkiiGTrr46TKk9caKSQqoUFcG4cWHo6pVXxh2NSF5QYsiU118PK7L95jdwxBFxR5NfunQJixpNnAgzZ8YdjUjOU1NSJmzeDAcfHL7VLlwIDRvGHVH++frrsLBRcXFYXmunneKOSCTrqCkpm9x2G8ybFyoGJYX02GUXuP320M8wenTc0YjkNCWGdFu1CoYOhT59wvrNkj7HHQf9+8MNN8C778YdjUjOUmJItyuugA0b4JZbNM4+E8aOhQYNwtoWOdhMKpINlBjSacaMsL7CoEGQppWWZDt77gmjRsH06TBlStzRiOQkdT6ny6ZNYU6fdetgwYLwLVYyY8uWcPPgmjWhz2GXXeKOSCQrqPM5bjfdFEYg3XyzkkKmFRWF33tZWageRKRGlBjSoawsrBtw/PHQt2/c0RSmH/4QTj8d/vSncFOhiFRbShKDmfUxs8VmttTMBlfw/tlmtsbM5kSP8xLeG2BmS6LHgFTEE7shQ8K9CzfdFHckhe2PfwyT7V12WdyRiOSUpBODmRUBtwHHAQcAp5nZARXs+rC7d40eE6JjdweGAYcChwDDzKxJsjHF6s03wzoLl12m9Zvj1qIF/OEPYRbbadPijkYkZ6SiYjgEWOruy9x9IzAZ6FfNY3sD09z9E3f/FJgG9ElBTPFwDwvV77knDP5O4SRxuPRS+P734fe/DwMCRKRKqUgMLYDlCa/Lom3b+5mZzTOzR82sVQ2PzQ1TpsCrr8L11+sO52xRXAxjxsCiRVpXW6SaMtX5/Hegrbt3IVQFk2p6AjMbaGalZla6Zs2alAeYtHXrwuyeXbvC2WfHHY0kOv54OPZYuPbaMF+ViOxQKhLDCqBVwuuW0bb/cve17r4hejkBOLi6xyacY5y7d3f37iUlJSkIO8XGjIH//Cf8LCqKOxpJZAZ/+Qt88UWo5kRkh1KRGN4COphZOzOrD5wKTE3cwcyaJ7w8AVgUPX8WONbMmkSdzsdG23LLqlVhJbETT4SePeOORirSqROccw7ceiv8+99xRyOS1ZJODO6+GbiQ8Ad9ETDF3ReY2QgzOyHa7SIzW2Bmc4GLgLOjYz8BriMkl7eAEdG23HL11WE+pD/9Ke5IZEdGjIB69cJIJRGplKbESNa8eaFf4ZJLQnOFZLfhw0Nfw2uvwf/8T9zRiGSUpsTIlMGDoXHjUDVI9rv88jCc+PLLNfuqSCWUGJIxfTo880xommiS2/flFYxddw1NSq+8Em58E5HvUFNSbbnDIYfARx+FRWGKi+ONR6pv82Y48EDYuDHMfFu/ftwRiWSEmpLS7ZFHoLQUrrtOSSHX1K0bBgosXQp33RV3NCJZRxVDbWzaBPvvDzvvDLNn676FXOQORx8Nb78dhq/qTnUpAKoY0mn8+PDHZNQoJYVcZRbuPVmzJiwHKiL/pcRQU19+GYY7/vjHYfF5yV2HHgonnRSalT7+OO5oRLKGEkNN/fWvsHp1mOvfLO5oJFnXXw9ff62V3kQSKDHUxJo18Oc/w8knhxFJkvsOOADOOitMlbF8edX7ixQAJYaa+OMf4ZtvwkgkyR/Dh4fO6BEj4o5EJCsoMVTXypXhW+UvfgH77Rd3NJJKbdrA+efD3XfD4sVxRyMSOyWG6ho5MtwYdc01cUci6XDVVdCgAQwdGnckIrFTYqiODz6AcePg3HPDMpGSf773vbAM6COPwKxZcUcjEislhuq47jqoU0cT5eW7yy+H3XdXVSgFT4mhKkuWwD33wAUXQMuWcUcj6dSoEVxxBTz9dJiWW6RAKTFUZfhw2GmnML225L8LL4SSEhg2LO5IRGKTksRgZn3MbLGZLTWz7/wFNbNLzWyhmc0zsxfMrE3Ce1vMbE70mLr9sbF6+2146CG46CLYY4+4o5FM2HVXGDQIpk2Dl1+OOxqRWCQ9iZ6ZFQHvAscAZYQlOk9z94UJ+xwJvOHu35jZBUBPd+8fvfeVu+9ak2tmbBK9n/0Mnn8e3nsvtD1LYfjmG9h77zAsefr0uKMRSZlMTqJ3CLDU3Ze5+0ZgMtAvcQd3n+7u30QvXweyv7F+1ix4/PEwUkVJobDsvDMMGQIzZsCLL8YdjUjGpSIxtAAS5xIoi7ZV5lzgmYTXxWZWamavm9mJlR1kZgOj/UrXrFmTXMTVMXRoSAiXXJL+a0n2GTgQWrQII5RycGp6kWRktPPZzM4EugN/StjcJiptTgfGmtneFR3r7uPcvbu7dy8pKUlvoK+8EpbsHDQojFSRwlNcHG56e+UVeO65uKMRyahUJIYVQKuE1y2jbeWYWS/gKuAEd9+wbbu7r4h+LgNmAN1SEFNyhg4Nnc2//W3ckUicfvlLaN1aVYMUnFQkhreADmbWzszqA6cC5UYXmVk34C5CUlidsL2Jme0UPW8GHA4sJE4zZoQOxyFDYJddYg1FYrbTTuFLwptvhnsbRApESpb2NLP/BcYCRcDd7n6DmY0ASt19qpk9D3QGVkWH/MfdTzCzHxASxlZCkhrr7hOrul5aRyX17AnvvgvLlmktZwnLuO63HzRuHNb41hocksOqOyqpbiou5u5PA09vt+2ahOe9KjnuVULCyA4zZsBLL8HNNyspSFCvXmhKOvtsePJJOLHS8REieSMlFUOmpa1iOPLIMO2yqgVJtHkzdOwYmpbmzAnzZonkoEzex5AfXnopVAyDByspSHl164YpMubPh0cfjTsakbRTxbDNUUfBokWhWmjQILXnlty3ZQt07hz6GObPV9UgOUkVQ03MnBlGIg0erKQgFSsqClXDwoWqGiTvqWIAOPro8A9e1YLsyJYt0KVLeK6qQXKQKobqevnlMB/OoEFKCrJjRUVhhJKqBslzqhh69YIFC1QtSPWoapAcpoqhOl5+GV54Aa68UklBqkdVgxSAwq4YevUKi/EsWxamWhapDlUNkqNUMVTln//8tlpQUpCaUNUgea5wK4ZjjoF588LqbEoMUlOqGiQHqWLYkVdeCUt2qlqQ2lLVIHmsMCuGY4+FuXND34Km1pbaUtUgOUYVQ2VefRWmTQvVgpKCJENVg+SpwqsYeveG2bND34ISgyRLVYPkEFUMFXnttbB+r6oFSRVVDZKHUpIYzKyPmS02s6VmNriC93cys4ej998ws7YJ7w2Jti82s96piKdS114LJSVwwQVpvYwUmJNPhgMOCP9/bd0adzQiSUs6MZhZEXAbcBxwAHCamR2w3W7nAp+6e3tgDDA6OvYAwhrRHYE+wO3R+dJj2DC44w5VC5JaiVXDI4/EHY1I0lKxtOchwFJ3XwZgZpOBfsDChH36AcOj548Ct5qZRdsnu/sG4D0zWxqd77UUxPUdFz98GHPmALek4+xS0PznsHM7+CVwu2ttaEmLrl1h7Nj0XycVTUktgOUJr8uibRXu4+6bgc+BptU8FgAzG2hmpWZWumbNmhSELZJCZtCmLXzzDej/T8lxqagYMsLdxwHjIIxKqs05MpFppYBtaQpdfhuevzAvNDGJpMq6dWEJYu+d9oo0FRXDCqBVwuuW0bYK9zGzusBuwNpqHiuSGzRCSdJp/Hg47jh48820XyoVieEtoIOZtTOz+oTO5Knb7TMVGBA9Pxl40cMNFFOBU6NRS+2ADkD6P7VIumwboTRiRLjHQSQV1q+H0aPhiCPg0EPTfrmkE0PUZ3Ah8CywCJji7gvMbISZnRDtNhFoGnUuXwoMjo5dAEwhdFT/A/itu+tfk+QuVQ2SDhMmwMqVMHx4Ri5XeHc+i6Rb4t3Q89TXIElavx723js8Xnopqf4F3fksEhdVDZJKEyeGamHYsIwNg1bFIJIOqhokFTZsgPbtoU2bsBRxkolBFYNInIqKwjc8VQ2SjLvvhrKyjFYLoIpBJH22boXOncNzVQ1SUxs2QIcO0KpVWIo4BYlBFYNI3OrUUdUgtXfPPbB8ecarBVDFIJJeqhqkNjZuDNXCXnuFxcVSlBhUMYhkA1UNUhv33AP/+U8s1QKoYhBJP1UNUhMbN8I++8Aee8Drr6c0MahiEMkWqhqkJu69Fz74ILZqAVQxiGSGqgapjk2bQrVQUgJvvJHyxKCKQSSbqGqQ6rj3Xnj//VirBVDFIJI5qhpkRzZuhH33hWbNwtTaaUgMqhhEso2qBtmRv/0tVAsjRsS+NKwqBpFMUtUgFdk2J1KrVvDKK2lLDKoYRLKRqgapyPjxYU6kLKgWQBWDSOapapBE69aFtRY6dIAZM9KaGDJSMZjZ7mY2zcyWRD+bVLBPVzN7zcwWmNk8M+uf8N49Zvaemc2JHl2TiUckJ6hqkER33QWrVmVNtQBJVgxm9kfgE3cfZWaDgSbuPmi7ffYB3N2XmNlewCxgf3f/zMzuAZ5y9xr961DFIDlPVYMAfP01fP/74f+F559P++Uy1cfQD5gUPZ8EnLj9Du7+rrsviZ6vBFYDJUleVyS3qWoQgNtvh9WrQ7WQRZKtGD5z98bRcwM+3fa6kv0PISSQju6+NaoYDgM2AC8Ag919QyXHDgQGArRu3frgDz74oNZxi2QFVQ2F7csvoV076NEDnnkmI5dMWcVgZs+b2dsVPPol7uchw1SaZcysOXAfcI67b402DwH2A3oAuwODKjkcdx/n7t3dvXtJiQoOyQOqGgrbLbfA2rVw7bVxR/IdyVYMi4Ge7r4q+sM/w933rWC/RsAMYGRl/Qlm1hO43N37VnVd9TFI3lDVUJg+/zxUC4cfDn//e8Yum6k+hqnAgOj5AODJCgKpDzwB3Lt9UoiSybZmqBOBt5OMRyS3qGooTGPHwqefZl3fwjbJVgxNgSlAa+AD4Ofu/omZdQfOd/fzzOxM4G/AgoRDz3b3OWb2IqEj2oA50TFfVXVdVQySV1Q1FJZPP4W2beHoo+HxxzN66epWDHWTuYi7rwWOrmB7KXBe9Px+4P5Kjj8qmeuL5IVtVUP//vDII3DqqXFHJOn0l7/AF19kZd/CNpoSQyQbnHwydOoE11wT5uSX/PTRR6EZqX//b6vELKTEIJIN6tSBG26AJUvCer+Sn0aOhPXr4brr4o5kh5QYRLLF8cfDYYeFJoZ16+KORlLt/ffhjjvg3HPDvEhZTIlBJFuYwahRsGIF3HZb3NFIqg0bFgYWXHNN3JFUSYlBJJsccQT06QM33hjGukt+ePttuO8++N3voEWLuKOpkhKDSLYZORI++QT+/Oe4I5FUufpqaNgQBlU6uUNWUWIQyTbduoVRK2PGhFEskttefx2efBKuvBKaNo07mmpRYhDJRtddF0av3HBD3JFIMtxh8GD43vfg97+PO5pqU2IQyUYdOsB558Gdd8J778UdjdTWc8/BSy/B0KGw665xR1NtSgwi2eqaa6Bu3dA+Lblny5bQfNS2LQwcGHc0NaLEIJKt9toLLrsMHnwQ3nor7mikpu67L8x9NWoU1K8fdzQ1ktQkenHRJHpSML78Etq3h/32S/tC8ZJC33wTmgNbtgydz1ny3y1T026LSDo1bBjuhJ45M6Pz9kuSxoyBlSvDkOMsSQo1oYpBJNtt3hwmXHOH+fOhXr24I5Id+eijUOX16gVPPBF3NOWoYhDJF3Xrwh//CIsXw4QJcUcjVbn22jDUePTouCOptaQSg5ntbmbTzGxJ9LNJJfttMbM50WNqwvZ2ZvaGmS01s4ej1d5EZHt9+8KPfxzm2/nii7ijkcq88w6MGwe//jXss0/c0dRashXDYOAFd+8AvBC9rsg6d+8aPU5I2D4aGOPu7YFPgXOTjEckP5mF9uo1a0L1INlp0CDYeeeQwHNYsomhHzApej6JsG5ztUTrPB8FbFvotkbHixSc7t3hjDPCCmD/+U/c0cj2pk+HqVPhD3+AkpK4o0lKsolhD3dfFT3/ENijkv2KzazUzF43s21//JsCn7n75uh1GZD90w6KxGnkyFA9XHll3JFIos2bw5QXbdrk1NQXlalyzWczex7Ys4K3rkp84e5uZpUNcWrj7ivM7PvAi2Y2H6jRnMJmNhAYCNC6deuaHCqSP1q3Ds0Vw4fDBReEfgeJ3113hRFjjz0GDRrEHU3SkhquamaLgZ7uvsrMmgMz3H3fKo65B3gKeAxYA+zp7pvN7DBguLv3ruq6Gq4qBW3dOth/f9htN5g1K4wRM5aKAAANVklEQVRakvisXRtuZuvWDZ5/PqvvW8jUcNWpwIDo+QDgyQoCaWJmO0XPmwGHAws9ZKTpwMk7Ol5EttOgQehnmDcPxo+POxoZOjSMFLvppqxOCjWRbGIYBRxjZkuAXtFrzKy7mW0bcL0/UGpmcwmJYJS7L4zeGwRcamZLCX0OE5OMR6Qw/PSncOSRYYK9tWvjjqZwzZ0bmpF++1vo1CnuaFJGdz6L5Kr580Pzxa9/rTWi4+AOPXvCwoXw7rvQpMLbuLKK7nwWyXedO8NvfhPWbJg7N+5oCs+UKWEOq5EjcyIp1IQqBpFc9umnoeOzY0fNvppJX38dZrwtKQlTohcVxR1RtahiECkETZrAjTeGb66TJlW9v6TGsGFQVga33pozSaEmlBhEct2558Lhh8Pll8PHH8cdTf6bPRvGjg19Oz/4QdzRpIUSg0iuq1MnjIz54ouQHCR9tmwJy3Q2axYqtTylxCCSDzp2DNNkTJoEL74YdzT56/bbobQ0VAx51uGcSJ3PIvli3bowUqlOnXDzW3Fx3BHll7KycMf5D38ITz+dkx396nwWKTQNGoShq0uWhCGUkloXXRSakm6/PSeTQk0oMYjkk1694MwzYdQoWLQo7mjyx5NPhmU6hw2Ddu3ijibtlBhE8s1f/gING8J554VvuJKczz4LU1507gyXXhp3NBmhxCCSb773vTCh26uvwl//Gnc0ue+ii+DDD2HiRKhXL+5oMkKJQSQfnXFGmGjv6qvh7bfjjiZ3PfEE3HcfXHUV9OgRdzQZo8Qgko/MQkf0brvBWWfBpk1xR5R7Vq8ON7F16xYSQwFRYhDJVyUlMG5cuFP3+uvjjia3uIcV8j7/HO69F+rXjzuijFJiEMlnJ54Iv/gF3HBDmOxNqueBB+Dxx+G66/JqnYXq0g1uIvnus8/CiJpdd4V//Ssv1iROq7KykAw6dYKXXsqrSfIycoObme1uZtPMbEn08zv3iJvZkWY2J+Gx3sxOjN67x8zeS3ivazLxiEgFGjcOI2reeSdMmyGV27IFzj479Mncc09eJYWaSLYpaTDwgrt3AF6IXpfj7tPdvau7dwWOAr4BnkvY5Ypt77v7nCTjEZGKHHssXHJJmCb60UfjjiZ7XXcdvPAC3HwztG8fdzSxSTYx9AO2TQI/CTixiv1PBp5x92+SvK6I1NSoUXDoofDLX8LSpXFHk32efx5GjAijuH75y7ijiVWyiWEPd18VPf8Q2KOK/U8FHtpu2w1mNs/MxpjZTpUdaGYDzazUzErXrFmTRMgiBap+/bAcZb16cMopsH593BFlj5Ur4fTTwyR5BTAXUlWqTAxm9ryZvV3Bo1/ifh56sSvtyTaz5kBn4NmEzUOA/YAewO7AoMqOd/dx7t7d3buXlJRUFbaIVKR16zD8cs4cuPjiuKPJDps3w6mnhuU6H30Udtkl7ohiV7eqHdy9V2XvmdlHZtbc3VdFf/hX7+BUPweecPf/3mmTUG1sMLO/AVplRCTdfvITGDQIRo+GI44I35QL2TXXwMsvhzuc998/7miyQrJNSVOBAdHzAcCTO9j3NLZrRoqSCWZmhP4J3bsvkgnXXx/WFRg4EBYujDua+Dz1VFiJ7Ve/CrPSCpB8YhgFHGNmS4Be0WvMrLuZTdi2k5m1BVoBL213/ANmNh+YDzQDdHumSCbUrQuTJ4d7G044AdaujTuizFu4MFRL3bqFSQflv3SDm0ghe+016NkzLGr/3HMFM3soa9fCIYeEfoW33oJWreKOKCO0gpuIVO2ww8LNbzNmhDUHcvCLYo1t2AA/+xmsWBEW4CmQpFATVXY+i0ieO/PMsNrbyJFh1NLVV8cdUfps3RrubH7ppTAf0qGHxh1RVlJiEJHQGV1WBkOHQosWcM45cUeUHldeGfpWRo/WaKwdUGIQkXBD1/jxsGpVGKHTpEmYmTWfjB4dlj298EK44oq4o8lq6mMQkaB+fXjsMejeHfr3h3/8I+6IUueWW2DwYDjtNBg7tuDvbK6KEoOIfKthw5AQOnaEk04KE8rlunHjwrrNJ50EkyYV7IypNaHEICLlNW4chq62bx/ukv5//y/uiGrvppvC8pz/+7/w0EOFMxw3SUoMIvJdzZqFIaydOoVv2rk2Vbd7GGV18cXw05+G1dh2qnSOTtmOEoOIVKxp09CUdMgh8POf587dwZs3hw7mq66CM86Ahx9WUqghJQYRqdxuu4VmpRNPDN++L7oorHKWrb78Evr1C1NnX3llmEm2rgZf1pQSg4js2M47wyOPwKWXhtE9vXtDNq6J8s474Ya1Z5+FO+8Mw1Pr6E9cbei3JiJVKyoK9wDcfTe88gocdBC8+mrcUX1ryhTo0QM+/jgkhl//Ou6IcpoSg4hU3znnhIRQrx786Efwhz+EuYfi8umnoR+hf//QUf6vf8HRR8cXT55QYhCRmunWLawAd845YS2DHj3gn//MbAzuoXmrU6dQLVx7LcycCS1bZjaOPKXEICI116gRTJgAf/87fPZZqB7OOAM++CD91549G446KoyUKimB118Pq7DpHoWUUWIQkdrr2zfMzDp0aJhOo337MNfSsmWpv1ZpaRhxdNBBMG9eGHk0axYcfHDqr1XgkkoMZnaKmS0ws61mVuniD2bWx8wWm9lSMxucsL2dmb0RbX/YzOonE4+IxGCXXWDECFiyJHT63ntvSBC9e4fmnm++qf25P/kkTGnRo0d4zJwZmo2WLoULLtD0FmmS1ApuZrY/sBW4C7jc3b+zrJqZFQHvAscAZcBbwGnuvtDMpgCPu/tkM7sTmOvud1R1Xa3gJpLFVq4MM7VOnAjLl0NxcWj6OfLI8O3+wAPD7K3bT2TnDh9+GJqKZs2CadPCCKitW0Nfwq9+BQMGhHsrpFaqu4JbSpb2NLMZVJ4YDgOGu3vv6PWQ6K1RwBpgT3ffvP1+O6LEIJIDtmwJ02o89VR4LF367XsNGkDz5iFpAHz1VUgKGzd+u0/XrnD88WFN6oMP1oyoKVDdxJCJWwJbAMsTXpcBhwJNgc/cfXPC9haVncTMBgIDAVq3bp2eSEUkdYqKwtDRo4+GMWPCTXGzZsGCBaGqSEwEDRrAXnuFRYK6dAkjnxo1ijf+AlZlYjCz54E9K3jrKnd/MvUhVczdxwHjIFQMmbquiKRISQn06RMektWqTAzu3ivJa6wAElfbbhltWws0NrO6UdWwbbuIiMQoE8NV3wI6RCOQ6gOnAlM9dG5MB06O9hsAZKwCERGRiiU7XPUkMysDDgP+n5k9G23fy8yeBoiqgQuBZ4FFwBR3XxCdYhBwqZktJfQ5TEwmHhERSV5KRiVlmkYliYjUXHVHJenOZxERKUeJQUREylFiEBGRcpQYRESknJzsfDazNUBt5vdtBnyc4nAyLdc/Q67HD/oM2UKfoebauHtJVTvlZGKoLTMrrU6PfDbL9c+Q6/GDPkO20GdIHzUliYhIOUoMIiJSTqElhnFxB5ACuf4Zcj1+0GfIFvoMaVJQfQwiIlK1QqsYRESkCkoMIiJSTkEkBjPrY2aLzWypmQ2OO57aMLO7zWy1mb0ddyy1YWatzGy6mS00swVm9vu4Y6opMys2szfNbG70Ga6NO6baMLMiM5ttZk/FHUttmdn7ZjbfzOaYWc7NqGlmjc3sUTN7x8wWRUsbZ42872MwsyLgXeAYwvKhbwGnufvCWAOrITM7AvgKuNfdO8UdT02ZWXOgubv/y8waArOAE3Ppv4OZGbCLu39lZvWAfwK/d/fXYw6tRszsUqA70Mjd+8YdT22Y2ftAd3fPyRvczGwS8LK7T4jWqdnZ3T+LO65tCqFiOARY6u7L3H0jMBnoF3NMNebuM4FP4o6jttx9lbv/K3r+JWFtjkrX+M5GHnwVvawXPXLqm5WZtQR+AkyIO5ZCZWa7AUcQrT/j7huzKSlAYSSGFsDyhNdl5NgfpHxjZm2BbsAb8UZSc1EzzBxgNTDN3XPtM4wFrgS2xh1Ikhx4zsxmmdnAuIOpoXbAGuBvUZPeBDPbJe6gEhVCYpAsYma7Ao8BF7v7F3HHU1PuvsXduxLWKD/EzHKmWc/M+gKr3X1W3LGkwA/d/SDgOOC3UVNrrqgLHATc4e7dgK+BrOr7LITEsAJolfC6ZbRNMixql38MeMDdH487nmREpf90oE/csdTA4cAJUfv8ZOAoM7s/3pBqx91XRD9XA08QmoxzRRlQllBtPkpIFFmjEBLDW0AHM2sXdfKcCkyNOaaCE3XcTgQWuftf446nNsysxMwaR88bEAY0vBNvVNXn7kPcvaW7tyX8O3jR3c+MOawaM7NdogEMRE0wxwI5M1rP3T8ElpvZvtGmo4GsGoRRN+4A0s3dN5vZhcCzQBFwt7sviDmsGjOzh4CeQDMzKwOGufvEeKOqkcOBXwDzozZ6gD+4+9MxxlRTzYFJ0Ui3OsAUd8/ZIZ85bA/gifBdg7rAg+7+j3hDqrHfAQ9EX1aXAefEHE85eT9cVUREaqYQmpJERKQGlBhERKQcJQYRESlHiUFERMpRYhARkXKUGEREpBwlBhERKef/AzeQmKkgNlduAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xnc1WP+x/HXp7u7BVHqRtopW0VRjDE/E6L8JmKGyR7DNMxgyFKNpYSmZisMZlKWDJJtGD+GomQQ7rSppGTprpDstPf5/XF9G+fOfXcvZ/me5f18PM6jc77L+X5OdD7nc13X97rM3REREdmiTtwBiIhIdlFiEBGRcpQYRESkHCUGEREpR4lBRETKUWIQEZFylBgkZ5lZazP72syKUvBed5vZDamISyTX1Y07AJGqmNl7wK7ApoTNe7n7B8AOsQSVA8zsbqDM3a+OOxbJLUoMkiuOc/cpcQchUgjUlCQ5y8zampmbWd3o9TQzu97MXjKzr8zsWTNrlnD8Q2b2oZl9YWbTzaxjNa9zdvSeo83sczNbamY/jLYvM7OPzax/wvE7mdkEM1tlZu+b2dVmVqeW71XfzP5kZh+Y2Udm9jczaxjt62FmZWZ2WXTeSjM7J9o3ADgduDJqbvtXtH2QmS2P/n4WmdlRyf+XkHyjxCD55jTgHGAXoB5wecK+p4EO0b43gPtq8L6HAHOBpsD9wESgO9AeOAP4q5ltada6BdgJ2AP4MXBWFFNt3msksBfQJdrfArg24b12i67VAjgXuNXMmrj72Ojz/cHdd3D348xsb+BCoLu7NwJ6Ae/V4O9ACoQSg+SKf0a/sD83s39u47i73P1td18DTCJ8oQLg7ne6+1fuvg4YBhxgZjtV8/rvuvtd7r4JeBBoBQx393Xu/iywHmgfdYSfAgyJrvUe8GfgzFq8lwEDgEvd/VN3/woYEb3/Fhuicze4+1PA18DelXyGTUB9YD8zK3b399z9nWp+fikg6mOQXHFCNfsYPkx4/i1R53T0hX0jcDJQAmyOjmkGfFGN9/0o4fkaAHffetsO0fsVA+8n7Huf8Iu+pu9VAmwHzAw5AgADEkdhrXb3jQmv//uZt+buS8zsEkJS7GhmzwAD3X1FRcdL4VLFIIXiNKAv0JPQ9NI22m6VnVBLnxB+xbdJ2NYaWF7L91oDdHT3xtFjJ3ev7kis702d7O73u/uPovgcGFWLuCTPKTFIoWgErANWE36Fj0jHRaLmoUnAjWbWyMzaAAOBf9TivTYDdwCjzWwXADNrYWa9qvkWHxH6OYjO3dvMjjSz+sBaQtLZXNnJUriUGKRQTCA06SwHFgAz0niti4BvgKXAfwgdzHfW8r0GAUuAGWb2JTCFyvsQtjae0J+wpV+mPqEz+xNCk9suwJBaxiV5zLRQj4iIJFLFICIi5SgxiIhIOUoMIiJSjhKDiIiUk5M3uDVr1szbtm0bdxgiIjll5syZn7h7SVXH5WRiaNu2LaWlpXGHISKSU8zs/aqPUlOSiIhsRYlBRETKUWIQEZFycrKPQUQkkzZs2EBZWRlr166NO5RqadCgAS1btqS4uLhW5ysxiIhUoaysjEaNGtG2bVsSpkDPSu7O6tWrKSsro127drV6j5Q0JZnZndHSgm9Wst/M7GYzW2Jmc83swIR9/c1scfToX9H5IiJxWrt2LU2bNs36pABgZjRt2jSp6iZVfQx3A723sf9YwpKKHQgrUt0OYGY7A0MJSx0eDAw1syYpiklEJGVyISlskWysKWlKcvfpZtZ2G4f0BSZ4mMp1hpk1NrPmQA9gsrt/CmBmkwkJ5oFUxCU1sG4dzJ8P77wDK1fC6tWwZebdnXeG5s2hXTvYf39o0CDeWEUkrTLVx9ACWJbwuizaVtn27zGzAYRqg9atW6cnykKyeTO89ho8+SQ88wzMmQMbNpQ/xuy75LBFURF06gTHHAPHHQeHHgp11VUlkk9y5l+0u48FxgJ069ZNi0jUVlkZ3HlneLz/fvii/+EPYeBAOOgg2Htv2H33UCXUqRMSw2efwYoV8PbbMHMmzJgBY8bAH/8Yjj3nHDj33FBRiEjOy9R9DMuBVgmvW0bbKtsuqfbOO3DeeeHLe+hQ6NAB/vEPWLUKpk+HkSPh5JNDU1GzZiEpQKgadt45VAk//SnceCM89xx88glMmgRdu8Lvfw/t28MZZ8DChfF+TpE8d9ttt9GpUyfatGnDLbfckpZrZKpieAK40MwmEjqav3D3lWb2DDAiocP5GLTUYGp9+ilcey387W+hyef88+HSS2GPPao+d1t23DEkkpNPhmXL4Oab4bbb4P774eyzQ7LYddeUfASRrHLJJTB7dmrfs0uXUIVX4ZFHHmHy5MnMmjWLTz75hM6dO3PBBRdQN8XNuakarvoA8Aqwt5mVmdm5Zna+mZ0fHfIUYf3bJYTFzX8NEHU6Xw+8Hj2Gb+mIliS5w/jxsNdecPvtMGAALF0Kt9ySfFLYWqtWoVnpvfdCk9Q//hGue/PNoS9DRFLi5ptvZtSoURQXF9O8eXOKi4vZnI5/Y+6ec4+DDjrIZRuWL3c/9lh3cD/8cPc5czJ7/bfecu/V67vrL12a2euLpNiCBQviDsHXr1/vO+20039fr1ixwjt27Fjp8RXFDJR6Nb5jNVdSvvnXv0J/wLRpoTqYOjX0G2TS3nvD00/DXXeFknv//eEBjUAWScaCBQv48ssvWbp0KZs3b2bIkCFcfPHFabmWEkO+2LQJrroKjj8e2rYNX8gXXvhdJ3KmmYW+hnnz4IAD4LTT4OKLYf36eOIRyXGzZs3i9NNP59RTT2X//fendevWDBgwIC3XypnhqrINX30F/fqFX+nnngt//Wv23ITWunWoWq68MnSuzZ4Njz0GTZvGHZlITpk9ezZ9+vShX79+ab+WKoZct3w5/M//wLPPhk7mceOyJylsUVwMo0fDfffBq6+G+ybeeSfuqERyyuzZs+nSpUtGrqWKIZctWgQ9e8IXX4Q7mHtva7qqLHDaaaGC6NsXfvCDkMy6do07KpGcMG3atIxdSxVDrpo7Fw4/PLTZT5+e/Ulhix/9CF55BbbbDo44ItxFLSJZRYkhF82cCT16hCaa6dPDzTG5ZK+9QtzNmoWKZ/r0uCMSkQRKDLlm3rwwgd2OO8KLL4ahobmoTZsQf6tW8JOfhL4HEckKSgy5ZEufQoMG8PzzuT9pXfPmYd6lXXcNTWGzZsUdkYigxJA7li+Ho48Oz597LvXTWsRl993D52nUCHr1giVL4o5IpOApMeSCL76AY4+Fzz+Hf/8b9tkn7ohSq00bmDIlzKvUuzd8/HHcEYkUNCWGbLd+fZjueuFCeOSR/B3euddeYTqP5cuhTx/45pu4IxIpWEoM2cwdfv3r0J8wfvx3TUn56tBDYeJEKC2F/v01M6tITJQYstlNN4WEcNVVcNZZcUeTGX37him8H3kErrsu7mhEssrixYvp0aMHnTp14rLLLmPPPfdMy3V053O2euYZuOwyOPFEGD487mgya+BAmD8/fO6OHeHnP487IpH/imudnk2bNnHWWWdx6623cuCBB3LRRRfRsWPH1AYSUWLIRu++C6eeGqbPnjAhvhlS42IW5n1atChMCnjAAbl7v4ZIivzzn/9kv/3248ADDwRg3333pXHjxmm5VkoSg5n1Bm4CioBx7j5yq/2jgSOil9sBu7h742jfJmBetO8Ddz8+FTHlrLVr4aSTQvv6o4/CDjvEHVE86teHBx8Mne0nnxymzthuu7ijEqnOCpxpMWvWrHKT6M2ZM4eePXum5VpJ/xQ1syLgVuBYYD/gVDPbL/EYd7/U3bu4exfgFuDRhN1rtuwr+KQAYQ2FN96Ae++FNLUf5oyWLcMyoW++CRddFHc0IrFq2rQpb731FgCvvvoqEyZM4IADDkjLtVLRRnEwsMTdl7r7emAi0Hcbx58KaDmvikycGDqbhwyB446LO5rs0KtX6Hy/80645564oxGJzZlnnklpaSmdO3fm0UcfpWnTprRv3z4t10pFU1ILYFnC6zLgkIoONLM2QDvg+YTNDcysFNgIjHT3f1Zy7gBgAEDr1q1TEHaWee89OP/8MGSz0DqbqzJsGLz0ElxwAXTrFjqkRQpMs2bNeDWaU2zZsmVMmzaNOmnqf8x0r+YpwMPuvilhWxt37wacBowxswrbT9x9rLt3c/duJSUlmYg1czZuhDPOCP0K990HdTUmoJyiIrj//jBx4Mknw9dfxx2RSKzmzJnD/mlcyz0ViWE50CrhdctoW0VOYatmJHdfHv25FJgG5OmtvdswYkT4RXz77bk/MV667LYbPPBAGKl06aVxRyMSqz59+nDHHXek7f1TkRheBzqYWTszq0f48n9i64PMbB+gCfBKwrYmZlY/et4MOAxYkIKYcsfLL4cbuc44A04/Pe5ostsRR4S1o8eNg6eeijsakbyVdGJw943AhcAzwEJgkrvPN7PhZpY4yugUYKK7e8K2fYFSM5sDTCX0MRROYvjii5AM2rSBW2+NO5rcMGwYdO4M550Hn34adzQieSkljdnu/hTw1Fbbrt3q9bAKznsZ6JyKGHLShRfCsmVhwZodd4w7mtxQv3646a979/D3d//9cUckBcLdMbO4w6iW8r+/a67AbqnNIk8+GcboX3VVGIkk1delCwwdGvocHnoo7mikADRo0IDVq1cn/YWbCe7O6tWradCgQa3fw3Lhg26tW7duXlpaGncYtffFF2HIZZMmYf3mevXijij3bNwIP/whLF0aboDbbbe4I5I8tmHDBsrKyli7dm3coVRLgwYNaNmyJcXFxeW2m9nMaBToNmlcZByuuAJWroTHHlNSqK26dUOTUteuMGAAPP54mGNJJA2Ki4tpV0AjBtWUlGnPPw933BFmEO3ePe5octs++8Dvfx8W+FFfg0jKqCkpk775JoyoKSqCOXM0KVwqbNoEhx0WmpTeegt23jnuiESyVnWbklQxZNLVV4cptcePV1JIlaIiGDs2DF298sq4oxHJC0oMmTJjRliR7de/hsMPjzua/LL//mFRo/HjYfr0uKMRyXlqSsqEjRvhoIPCr9oFC6BRo7gjyj/ffBMWNmrQICyvVb9+3BGJZB01JWWTW2+FuXNDxaCkkB7bbw+33Rb6GUaNijsakZymxJBuK1fCNddA795h/WZJn2OPhX794MYb4e23445GJGcpMaTbFVfAunVwyy0aZ58JY8ZAw4ZhbYscbCYVyQZKDOk0bVpYX2HQIEjTSkuyld12g5EjYepUmDQp7mhEcpI6n9Nlw4Ywp8+aNTB/fvgVK5mxaVO4eXDVqtDnsP32cUckkhXU+Ry3m24KI5BuvllJIdOKisLfe1lZqB5EpEaUGNKhrCysG3DccdCnT9zRFKYf/QhOOw3++MdwU6GIVFtKEoOZ9TazRWa2xMwGV7D/bDNbZWazo8d5Cfv6m9ni6NE/FfHEbsiQcO/CTTfFHUlh+8MfwmR7l10WdyQiOSXpxGBmRcCtwLHAfsCpZrZfBYc+6O5dose46NydgaHAIcDBwFAza5JsTLF67bWwzsJll2n95ri1aAG/+12YxXby5LijEckZqagYDgaWuPtSd18PTAT6VvPcXsBkd//U3T8DJgO9UxBTPNzDQvW77QaDv1c4SRwGDoQ99oDf/jYMCBCRKqUiMbQAliW8Lou2be1nZjbXzB42s1Y1PDc3TJoEL78MN9ygO5yzRYMGMHo0LFyodbVFqilTnc//Atq6+/6EquCemr6BmQ0ws1IzK121alXKA0zamjVhds8uXeDss+OORhIddxwccwxcd12Yr0pEtikViWE50Crhdcto23+5+2p3Xxe9HAccVN1zE95jrLt3c/duJSUlKQg7xUaPhg8+CH8WFcUdjSQygz//Gb78MlRzIrJNqUgMrwMdzKydmdUDTgGeSDzAzJonvDweWBg9fwY4xsyaRJ3Ox0TbcsvKlWElsRNOgB494o5GKtKpE5xzDvz1r/DOO3FHI5LVkk4M7r4RuJDwhb4QmOTu881suJkdHx12sZnNN7M5wMXA2dG5nwLXE5LL68DwaFtuufrqMB/SH/8YdySyLcOHQ3FxGKkkIpXSlBjJmjs39CtcemlorpDsNmxY6Gt45RX4wQ/ijkYkozQlRqYMHgyNG4eqQbLf5ZeH4cSXX67ZV0UqocSQjKlT4emnQ9NEk9y+L69g7LBDaFJ66aVw45uIfI+akmrLHQ4+GD76KCwK06BBvPFI9W3cCAccAOvXh5lv69WLOyKRjFBTUro99BCUlsL11ysp5Jq6dcNAgSVL4O9/jzsakayjiqE2NmyAffeF7baDWbN030IucoejjoI33wzDV3WnuhQAVQzpdMcd4ctk5EglhVxlFu49WbUqLAcqIv+lxFBTX30Vhjv++Mdh8XnJXYccAieeGJqVPvkk7mhEsoYSQ0395S/w8cdhrn+zuKORZN1wA3zzjVZ6E0mgxFATq1bBn/4EJ50URiRJ7ttvPzjrrDBVxrJlVR8vUgCUGGriD3+Ab78NI5EkfwwbFjqjhw+POxKRrKDEUF0rVoRflWeeCfvsE3c0kkpt2sD558Odd8KiRXFHIxI7JYbqGjEi3Bh17bVxRyLpcNVV0LAhXHNN3JGIxE6JoTrefx/GjoVzzw3LREr+2WWXsAzoQw/BzJlxRyMSKyWG6rj+eqhTRxPl5bvLL4edd1ZVKAVPiaEqixfD3XfDBRdAy5ZxRyPptOOOcMUV8NRTYVpukQKlxFCVYcOgfv0wvbbkvwsvhJISGDo07khEYpOSxGBmvc1skZktMbPvfYOa2UAzW2Bmc83sOTNrk7Bvk5nNjh5PbH1urN58Ex54AC6+GHbdNe5oJBN22AEGDYLJk+HFF+OORiQWSU+iZ2ZFwNvA0UAZYYnOU919QcIxRwCvuvu3ZnYB0MPd+0X7vnb3HWpyzYxNovezn8GUKfDuu6HtWQrDt9/CnnuGYclTp8YdjUjKZHISvYOBJe6+1N3XAxOBvokHuPtUd/82ejkDyP7G+pkz4dFHw0gVJYXCst12MGQITJsGzz8fdzQiGZeKxNACSJxLoCzaVplzgacTXjcws1Izm2FmJ1R2kpkNiI4rXbVqVXIRV8c114SEcOml6b+WZJ8BA6BFizBCKQenphdJRkY7n83sDKAb8MeEzW2i0uY0YIyZ7VnRue4+1t27uXu3kpKS9Ab60kthyc5Bg8JIFSk8DRqEm95eegmefTbuaEQyKhWJYTnQKuF1y2hbOWbWE7gKON7d123Z7u7Loz+XAtOArimIKTnXXBM6m3/zm7gjkTj94hfQurWqBik4qUgMrwMdzKydmdUDTgHKjS4ys67A3wlJ4eOE7U3MrH70vBlwGLCAOE2bFjochwyB7bePNRSJWf364UfCa6+FextECkRKlvY0s/8FxgBFwJ3ufqOZDQdK3f0JM5sCdAZWRqd84O7Hm9kPCQljMyFJjXH38VVdL62jknr0gLffhqVLtZazhGVc99kHGjcOa3xrDQ7JYdUdlVQ3FRdz96eAp7badm3C856VnPcyIWFkh2nT4IUX4OablRQkKC4OTUlnnw2PPw4nVDo+QiRvpKRiyLS0VQxHHBGmXVa1IIk2boSOHUPT0uzZYd4skRyUyfsY8sMLL4SKYfBgJQUpr27dMEXGvHnw8MNxRyOSdqoYtjjySFi4MFQLDRum9r0l923aBJ07hz6GefNUNUhOUsVQE9Onh5FIgwcrKUjFiopC1bBggaoGyXuqGACOOir8g1e1INuyaRPsv394rqpBcpAqhup68cUwH86gQUoKsm1FRWGEkqoGyXOqGHr2hPnzVS1I9ahqkBymiqE6XnwRnnsOrrxSSUGqR1WDFIDCrhh69gyL8SxdGqZaFqkOVQ2So1QxVOU///muWlBSkJpQ1SB5rnArhqOPhrlzw+psSgxSU6oaJAepYtiWl14KS3aqWpDaUtUgeawwK4ZjjoE5c0LfgqbWltpS1SA5RhVDZV5+GSZPDtWCkoIkQ1WD5KnCqxh69YJZs0LfghKDJEtVg+QQVQwVeeWVsH6vqgVJFVUNkodSkhjMrLeZLTKzJWY2uIL99c3swWj/q2bWNmHfkGj7IjPrlYp4KnXddVBSAhdckNbLSIE56STYb7/w/9fmzXFHI5K0pBODmRUBtwLHAvsBp5rZflsddi7wmbu3B0YDo6Jz9yOsEd0R6A3cFr1fegwdCrffrmpBUiuxanjoobijEUlaKpb2PBhY4u5LAcxsItAXWJBwTF9gWPT8YeCvZmbR9onuvg5418yWRO/3Sgri+p5LHjyU2bOBW9Lx7lLQ/OewXTv4BXCba21oSYsuXWDMmPRfJxVNSS2AZQmvy6JtFR7j7huBL4Cm1TwXADMbYGalZla6atWqFIQtkkJm0KYtfPst6P9PyXGpqBgywt3HAmMhjEqqzXtkItNKAdvUFPb/TXj+3NzQxCSSKmvWhCWIvVfaK9JUVAzLgVYJr1tG2yo8xszqAjsBq6t5rkhu0AglSac77oBjj4XXXkv7pVKRGF4HOphZOzOrR+hMfmKrY54A+kfPTwKe93ADxRPAKdGopXZAByD9n1okXbaMUBo+PNzjIJIKa9fCqFFw+OFwyCFpv1zSiSHqM7gQeAZYCExy9/lmNtzMjo8OGw80jTqXBwKDo3PnA5MIHdX/Bn7j7vrXJLlLVYOkw7hxsGIFDBuWkcsV3p3PIumWeDf0XPU1SJLWroU99wyPF15Iqn9Bdz6LxEVVg6TS+PGhWhg6NGPDoFUxiKSDqgZJhXXroH17aNMmLEWcZGJQxSASp6Ki8AtPVYMk4847oawso9UCqGIQSZ/Nm6Fz5/BcVYPU1Lp10KEDtGoVliJOQWJQxSAStzp1VDVI7d19NyxblvFqAVQxiKSXqgapjfXrQ7Ww++5hcbEUJQZVDCLZQFWD1Mbdd8MHH8RSLYAqBpH0U9UgNbF+Pey1F+y6K8yYkdLEoIpBJFuoapCamDAB3n8/tmoBVDGIZIaqBqmODRtCtVBSAq++mvLEoIpBJJuoapDqmDAB3nsv1moBVDGIZI6qBtmW9eth772hWbMwtXYaEoMqBpFso6pBtuWuu0K1MHx47EvDqmIQySRVDVKRLXMitWoFL72UtsSgikEkG6lqkIrccUeYEykLqgVQxSCSeaoaJNGaNWGthQ4dYNq0tCaGjFQMZrazmU02s8XRn00qOKaLmb1iZvPNbK6Z9UvYd7eZvWtms6NHl2TiEckJqhok0d//DitXZk21AElWDGb2B+BTdx9pZoOBJu4+aKtj9gLc3Reb2e7ATGBfd//czO4GnnT3Gv3rUMUgOU9VgwB88w3ssUf4f2HKlLRfLlN9DH2Be6Ln9wAnbH2Au7/t7ouj5yuAj4GSJK8rkttUNQjAbbfBxx+HaiGLJFsxfO7ujaPnBny25XUlxx9MSCAd3X1zVDEcCqwDngMGu/u6Ss4dAAwAaN269UHvv/9+reMWyQqqGgrbV19Bu3bQvTs8/XRGLpmyisHMppjZmxU8+iYe5yHDVJplzKw5cC9wjrtvjjYPAfYBugM7A4MqOR13H+vu3dy9W0mJCg7JA6oaCtstt8Dq1XDddXFH8j3JVgyLgB7uvjL64p/m7ntXcNyOwDRgRGX9CWbWA7jc3ftUdV31MUjeUNVQmL74IlQLhx0G//pXxi6bqT6GJ4D+0fP+wOMVBFIPeAyYsHVSiJLJlmaoE4A3k4xHJLeoaihMY8bAZ59lXd/CFslWDE2BSUBr4H3g5+7+qZl1A8539/PM7AzgLmB+wqlnu/tsM3ue0BFtwOzonK+ruq4qBskrqhoKy2efQdu2cNRR8OijGb10dSuGuslcxN1XA0dVsL0UOC96/g/gH5Wcf2Qy1xfJC1uqhn794KGH4JRT4o5I0unPf4Yvv8zKvoUtNCWGSDY46STo1AmuvTbMyS/56aOPQjNSv37fVYlZSIlBJBvUqQM33giLF4f1fiU/jRgBa9fC9dfHHck2KTGIZIvjjoNDDw1NDGvWxB2NpNp778Htt8O554Z5kbKYEoNItjCDkSNh+XK49da4o5FUGzo0DCy49tq4I6mSEoNINjn8cOjdG37/+zDWXfLDm2/CvffCRRdBixZxR1MlJQaRbDNiBHz6KfzpT3FHIqly9dXQqBEMqnRyh6yixCCSbbp2DaNWRo8Oo1gkt82YAY8/DldeCU2bxh1NtSgxiGSj668Po1duvDHuSCQZ7jB4MOyyC/z2t3FHU21KDCLZqEMHOO88+Nvf4N13445GauvZZ+GFF+Caa2CHHeKOptqUGESy1bXXQt26oX1acs+mTaH5qG1bGDAg7mhqRIlBJFvtvjtcdhncfz+8/nrc0UhN3XtvmPtq5EioVy/uaGokqUn04qJJ9KRgfPUVtG8P++yT9oXiJYW+/TY0B7ZsGTqfs+S/W6am3RaRdGrUKNwJPX16RuftlySNHg0rVoQhx1mSFGpCFYNIttu4MUy45g7z5kFxcdwRybZ89FGo8nr2hMceizuaclQxiOSLunXhD3+ARYtg3Li4o5GqXHddGGo8alTckdRaUonBzHY2s8lmtjj6s0klx20ys9nR44mE7e3M7FUzW2JmD0arvYnI1vr0gR//OMy38+WXcUcjlXnrLRg7Fn71K9hrr7ijqbVkK4bBwHPu3gF4LnpdkTXu3iV6HJ+wfRQw2t3bA58B5yYZj0h+Mgvt1atWhepBstOgQbDddiGB57BkE0Nf4J7o+T2EdZurJVrn+Uhgy0K3NTpfpOB06wannx5WAPvgg7ijka1NnQpPPAG/+x2UlMQdTVKSTQy7uvvK6PmHwK6VHNfAzErNbIaZbfnybwp87u4bo9dlQPZPOygSpxEjQvVw5ZVxRyKJNm4MU160aZNTU19Upso1n81sCrBbBbuuSnzh7m5mlQ1xauPuy81sD+B5M5sH1GhOYTMbAAwAaN26dU1OFckfrVuH5ophw+CCC0K/g8Tv738PI8YeeQQaNow7mqQlNVzVzBYBPdx9pZk1B6a5+95VnHM38CTwCLAK2M3dN5rZocAwd+9V1XU1XFUK2po1sO++sNNOMHNmGLUk8Vm9OtypJq8xAAANT0lEQVTM1rUrTJmS1fctZGq46hNA/+h5f+DxCgJpYmb1o+fNgMOABR4y0lTgpG2dLyJbadgw9DPMnQt33BF3NHLNNWGk2E03ZXVSqIlkE8NI4GgzWwz0jF5jZt3MbMuA632BUjObQ0gEI919QbRvEDDQzJYQ+hzGJxmPSGH46U/hiCPCBHurV8cdTeGaMyc0I/3mN9CpU9zRpIzufBbJVfPmheaLX/1Ka0THwR169IAFC+Dtt6FJhbdxZRXd+SyS7zp3hl//OqzZMGdO3NEUnkmTwhxWI0bkRFKoCVUMIrnss89Cx2fHjpp9NZO++SbMeFtSEqZELyqKO6JqUcUgUgiaNIHf/z78cr3nnqqPl9QYOhTKyuCvf82ZpFATSgwiue7cc+Gww+Dyy+GTT+KOJv/NmgVjxoS+nR/+MO5o0kKJQSTX1akTRsZ8+WVIDpI+mzaFZTqbNQuVWp5SYhDJBx07hmky7rkHnn8+7mjy1223QWlpqBjyrMM5kTqfRfLFmjVhpFKdOuHmtwYN4o4ov5SVhTvOf/QjeOqpnOzoV+ezSKFp2DAMXV28OAyhlNS6+OLQlHTbbTmZFGpCiUEkn/TsCWecASNHwsKFcUeTPx5/PCzTOXQotGsXdzRpp8Qgkm/+/Gdo1AjOOy/8wpXkfP55mPKic2cYODDuaDJCiUEk3+yyS5jQ7eWX4S9/iTua3HfxxfDhhzB+PBQXxx1NRigxiOSj008PE+1dfTW8+Wbc0eSuxx6De++Fq66C7t3jjiZjlBhE8pFZ6IjeaSc46yzYsCHuiHLPxx+Hm9i6dg2JoYAoMYjkq5ISGDs23Kl7ww1xR5Nb3MMKeV98ARMmQL16cUeUUUoMIvnshBPgzDPhxhvDZG9SPffdB48+Ctdfn1frLFSXbnATyXeffx5G1OywA7zxRl6sSZxWZWUhGXTqBC+8kFeT5GXkBjcz29nMJpvZ4ujP790jbmZHmNnshMdaMzsh2ne3mb2bsK9LMvGISAUaNw4jat56K0ybIZXbtAnOPjv0ydx9d14lhZpItilpMPCcu3cAnotel+PuU929i7t3AY4EvgWeTTjkii373X12kvGISEWOOQYuvTRME/3ww3FHk72uvx6eew5uvhnat487mtgkmxj6Alsmgb8HOKGK408Cnnb3b5O8rojU1MiRcMgh8ItfwJIlcUeTfaZMgeHDwyiuX/wi7mhilWxi2NXdV0bPPwR2reL4U4AHttp2o5nNNbPRZla/shPNbICZlZpZ6apVq5IIWaRA1asXlqMsLoaTT4a1a+OOKHusWAGnnRYmySuAuZCqUmViMLMpZvZmBY++icd56MWutCfbzJoDnYFnEjYPAfYBugM7A4MqO9/dx7p7N3fvVlJSUlXYIlKR1q3D8MvZs+GSS+KOJjts3AinnBKW63z4Ydh++7gjil3dqg5w956V7TOzj8ysubuvjL74P97GW/0ceMzd/3unTUK1sc7M7gK0yohIuv3kJzBoEIwaBYcfHn4pF7Jrr4UXXwx3OO+7b9zRZIVkm5KeAPpHz/sDj2/j2FPZqhkpSiaYmRH6J3Tvvkgm3HBDWFdgwABYsCDuaOLz5JNhJbZf/jLMSitA8olhJHC0mS0GekavMbNuZjZuy0Fm1hZoBbyw1fn3mdk8YB7QDNDtmSKZULcuTJwY7m04/nhYvTruiDJvwYJQLXXtGiYdlP/SDW4iheyVV6BHj7Co/bPPFszsoaxeDQcfHPoVXn8dWrWKO6KM0ApuIlK1Qw8NN79NmxbWHMjBH4o1tm4d/OxnsHx5WICnQJJCTVTZ+Swiee6MM8JqbyNGhFFLV18dd0Tps3lzuLP5hRfCfEiHHBJ3RFlJiUFEQmd0WRlccw20aAHnnBN3ROlx5ZWhb2XUKI3G2gYlBhEJN3TdcQesXBlG6DRpEmZmzSejRoVlTy+8EK64Iu5ospr6GEQkqFcPHnkEunWDfv3g3/+OO6LUueUWGDwYTj0Vxowp+Dubq6LEICLfadQoJISOHeHEE8OEcrlu7NiwbvOJJ8I99xTsjKk1ocQgIuU1bhyGrrZvH+6S/r//izui2rvpprA85//+LzzwQOEMx02SEoOIfF+zZmEIa6dO4Zd2rk3V7R5GWV1yCfz0p2E1tvqVztEpW1FiEJGKNW0ampIOPhh+/vPcuTt448bQwXzVVXD66fDgg0oKNaTEICKV22mn0Kx0wgnh1/fFF4dVzrLVV19B375h6uwrrwwzydbV4MuaUmIQkW3bbjt46CEYODCM7unVC7JxTZS33go3rD3zDPztb2F4ah19xdWG/tZEpGpFReEegDvvhJdeggMPhJdfjjuq70yaBN27wyefhMTwq1/FHVFOU2IQkeo755yQEIqL4X/+B373uzD3UFw++yz0I/TrFzrK33gDjjoqvnjyhBKDiNRM165hBbhzzglrGXTvDv/5T2ZjcA/NW506hWrhuutg+nRo2TKzceQpJQYRqbkdd4Rx4+Bf/4LPPw/Vw+mnw/vvp//as2bBkUeGkVIlJTBjRliFTfcopIwSg4jUXp8+YWbWa64J02m0bx/mWlq6NPXXKi0NI44OPBDmzg0jj2bOhIMOSv21ClxSicHMTjaz+Wa22cwqXfzBzHqb2SIzW2JmgxO2tzOzV6PtD5pZvWTiEZEYbL89DB8OixeHTt8JE0KC6NUrNPd8+23t3/vTT8OUFt27h8f06aHZaMkSuOACTW+RJkmt4GZm+wKbgb8Dl7v795ZVM7Mi4G3gaKAMeB041d0XmNkk4FF3n2hmfwPmuPvtVV1XK7iJZLEVK8JMrePHw7Jl0KBBaPo54ojw6/6AA8LsrVtPZOcOH34YmopmzoTJk8MIqM2bQ1/CL38J/fuHeyukVqq7gltKlvY0s2lUnhgOBYa5e6/o9ZBo10hgFbCbu2/c+rhtUWIQyQGbNoVpNZ58MjyWLPluX8OG0Lx5SBoAX38dksL69d8d06ULHHdcWJP6oIM0I2oKVDcxZOKWwBbAsoTXZcAhQFPgc3ffmLC9RWVvYmYDgAEArVu3Tk+kIpI6RUVh6OhRR8Ho0eGmuJkzYf78UFUkJoKGDWH33cMiQfvvH0Y+7bhjvPEXsCoTg5lNAXarYNdV7v546kOqmLuPBcZCqBgydV0RSZGSEujdOzwkq1WZGNy9Z5LXWA4krrbdMtq2GmhsZnWjqmHLdhERiVEmhqu+DnSIRiDVA04BnvDQuTEVOCk6rj+QsQpEREQqluxw1RPNrAw4FPg/M3sm2r67mT0FEFUDFwLPAAuBSe4+P3qLQcBAM1tC6HMYn0w8IiKSvJSMSso0jUoSEam56o5K0p3PIiJSjhKDiIiUo8QgIiLlKDGIiEg5Odn5bGargNrM79sM+CTF4WRarn+GXI8f9BmyhT5DzbVx95KqDsrJxFBbZlZanR75bJbrnyHX4wd9hmyhz5A+akoSEZFylBhERKScQksMY+MOIAVy/TPkevygz5At9BnSpKD6GEREpGqFVjGIiEgVlBhERKScgkgMZtbbzBaZ2RIzGxx3PLVhZnea2cdm9mbcsdSGmbUys6lmtsDM5pvZb+OOqabMrIGZvWZmc6LPcF3cMdWGmRWZ2SwzezLuWGrLzN4zs3lmNtvMcm5GTTNrbGYPm9lbZrYwWto4a+R9H4OZFQFvA0cTlg99HTjV3RfEGlgNmdnhwNfABHfvFHc8NWVmzYHm7v6GmTUCZgIn5NJ/BzMzYHt3/9rMioH/AL919xkxh1YjZjYQ6Abs6O594o6nNszsPaCbu+fkDW5mdg/woruPi9ap2c7dP487ri0KoWI4GFji7kvdfT0wEegbc0w15u7TgU/jjqO23H2lu78RPf+KsDZHpWt8ZyMPvo5eFkePnPplZWYtgZ8A4+KOpVCZ2U7A4UTrz7j7+mxKClAYiaEFsCzhdRk59oWUb8ysLdAVeDXeSGouaoaZDXwMTHb3XPsMY4Argc1xB5IkB541s5lmNiDuYGqoHbAKuCtq0htnZtvHHVSiQkgMkkXMbAfgEeASd/8y7nhqyt03uXsXwhrlB5tZzjTrmVkf4GN3nxl3LCnwI3c/EDgW+E3U1Jor6gIHAre7e1fgGyCr+j4LITEsB1olvG4ZbZMMi9rlHwHuc/dH444nGVHpPxXoHXcsNXAYcHzUPj8RONLM/hFvSLXj7sujPz8GHiM0GeeKMqAsodp8mJAoskYhJIbXgQ5m1i7q5DkFeCLmmApO1HE7Hljo7n+JO57aMLMSM2scPW9IGNDwVrxRVZ+7D3H3lu7elvDv4Hl3PyPmsGrMzLaPBjAQNcEcA+TMaD13/xBYZmZ7R5uOArJqEEbduANIN3ffaGYXAs8ARcCd7j4/5rBqzMweAHoAzcysDBjq7uPjjapGDgPOBOZFbfQAv3P3p2KMqaaaA/dEI93qAJPcPWeHfOawXYHHwm8N6gL3u/u/4w2pxi4C7ot+rC4Fzok5nnLyfriqiIjUTCE0JYmISA0oMYiISDlKDCIiUo4Sg4iIlKPEICIi5SgxiIhIOUoMIiJSzv8Dr7d+3DDbzusAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# parameters\n", "c = 1 # velocity for the transport equation\n", "Tf = 2.*np.pi # final time\n", "N = 128 # number of points in space\n", "la = 1. # scheme velocity\n", "s = 2. # relaxation parameter\n", "# initialization\n", "x = mesh(N) # mesh\n", "dx = x[1]-x[0] # space step\n", "dt = dx/la # time step\n", "f0, f1, f2, m0, m1, m2 = initialize(x, c, la)\n", "plt.figure(1)\n", "plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\\rho$')\n", "plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')\n", "plt.title('Initial moments')\n", "plt.legend(loc='best')\n", "# time loops\n", "nt = int(Tf/dt)\n", "m2f(f0, f1, f2, m0, m1, m2, la)\n", "for k in range(nt):\n", " transport(f0, f1, f2)\n", " f2m(f0, f1, f2, m0, m1, m2, la)\n", " relaxation(m0, m1, m2, c, s)\n", " m2f(f0, f1, f2, m0, m1, m2, la)\n", "plt.figure(2)\n", "plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\\rho$')\n", "plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')\n", "plt.title('Final moments')\n", "plt.legend(loc='best')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Anti bounce back conditions\n", "\n", "In order to take into account homogenous Dirichlet conditions over $\\rho$, we introduce the bounce back conditions. At edge $x=0$, two points are involved: $x_0=-\\dx/2$ and $x_1=\\dx/2$. We impose $\\fk{1}(x_0)=-\\fk{2}(x_1)$. And at edge $x=2\\pi$, the two involved points are $x_{N}$ and $x_{N+1}$. We impose $\\fk{2}(x_{N+1})=-\\fk{1}(x_{N})$.\n", "\n", "We modify the transport function to impose anti bounce back conditions. We can compare the solutions obtained with the two different boundary conditions." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVdXV//HPYgAHFQRhokiPYKMICvoYE4OKgk9QNNFgi2g0RBNj7EAUQVSENLArxYgVsT0SfxpFBTH2ITQBEYIaBlAQu9JZvz/2Id7BGabccm75vl+v+5p7zz1l3VFm3bX3PnubuyMiIrJNnbgDEBGR7KLEICIi5SgxiIhIOUoMIiJSjhKDiIiUo8QgIiLlKDFIzjOzZ8xswA7ev9PMhlbzXDPM7LzURSeSe5QYJCuZ2ftm1qs6+7r7ce4+KTrubDP753bvn+/u16UjzlygZCc1pcQgIiLlKDFI1ttWBZjZn83sUzN7z8yOS3h/hpmdZ2b7A3cCh5nZV2b2WfT+PWZ2ffS8iZk9ZWZronM9ZWYtqxnHcDN7xMzuN7MvzWy+me1jZkPMbLWZLTezYxP238vMpprZJ2a21Mx+lcS5djOziWa2ysxWmNn1ZlZU1e/HzG4AfgTcGv1ObrVgTHSdL6Jrd0rmv5HkFyUGyRWHAouBZsAfgYlmZok7uPsi4HzgNXff1d0bV3CeOsDfgDZAa2AdcGsN4jgeuA9oAswGno3O2QIYAdyVsO9koAzYCzgZGGlmR9XyXPcAm4H2QDfgWCCxeajC34+7XwW8DFwY/U4ujI49AtgH2A34ObC2Br8DyXNKDJIrPnD38e6+BZgENAf2qOlJ3H2tuz/m7t+4+5fADcCPa3CKl939WXffDDwClACj3H0TIRG0NbPGZtYKOBwY5O7r3X0OMAE4qxbn2gP4X+Bid//a3VcDY4BTa/n72QQ0BPYDzN0XufuqGvwOJM/VjTsAkWr6cNsTd/8mKhZ2relJzGxnwh/VPoRv6gANzawo+qNalY8Snq8DPk44bl30c1dClfBJlHy2+QDoXstz1QNWJRRJdYDlCcdX+/fj7i+a2a3AbUAbM3scuNzdv6jwE0vBUcUg+aaq6YIvA/YFDnX3RoQmFQCr/JBaWQnsbmYNE7a1BlbU4lzLgQ1AM3dvHD0auXvHah7/nd+Ju9/s7gcDBxCalK6oRVySp5QYJN98BLQ0s/qVvN+Q8G38MzPbHRiWjiDcfTnwKnCjmRWbWRfgXOD+WpxrFfAc8Bcza2RmdcxsbzOrbhPYR8D3t70wsx5mdqiZ1QO+BtYDW2sal+QvJQbJNy8CC4APzezjCt4fCzQAPgZeB/6RxlhOA9oSqocngGHu/nwtz3UWUB9YCHwKPEroR6iOm4CToxFLNwONgPHReT4gdDz/qZZxSR4yLdQjIiKJVDGIiEg5SgwiIlKOEoOIiJSjxCAiIuXk5A1uzZo187Zt28YdhohITpk1a9bH7l5S1X45mRjatm1LaWlp3GGIiOQUM/ugOvupKUlERMpRYhARkXKUGEREpJyc7GMQEcmkTZs2UVZWxvr16+MOpVqKi4tp2bIl9erVq9XxSgwiIlUoKyujYcOGtG3blu3Wh8o67s7atWspKyujXbt2tTpHSpqSzOzuaJnAtyt538zs5mh5w3lmdlDCewPMbEn0GJCKeEREUmn9+vU0bdo065MCgJnRtGnTpKqbVPUx3ENY+KQyxwEdosdA4A6AhGmPDwUOAYaZWZPKTiIiEpdcSArbJBtrSpqS3H2mmbXdwS79gHs9TOX6erRcYXOgJzDN3T8BMLNphATzUCrikhrYsAEWLIB//xtWrYK1a2HbzLu77w7Nm0O7dtClCxQXxxuriKRVpvoYWlB+GcKyaFtl27/DzAYSqg1at26dnigLydat8Oab8NRT8OyzMHcubNpUfh+zb5PDNkVF0KkTHHssHH88HHYY1FVXlUg+yZl/0e4+DhgH0L17dy0iUVtlZXD33eHxwQfhD/0PfgCXXgoHHwz77gt77RWqhDp1QmL49FNYuRLefRdmzYLXX4exY+FPfwr7nnMOnHtuqChEJOdl6j6GFUCrhNcto22VbZdU+/e/4bzzwh/vYcOgQwe4/35YswZmzoRRo+CUU0JTUbNmISlAqBp23z1UCT/9KdxwA7zwAnz8MUyZAt26wY03Qvv2cOaZsGhRvJ9TJM/dfvvtdOrUiTZt2nDLLbek5RqZqhimAhea2WRCR/Pn7r7KzJ4FRiZ0OB8LDMlQTIXhk0/gmmvgzjtDk8/558Mll8D3v1/1sTvSqFFIJKecAsuXw803w+23w4MPwtlnh2Sxxx4p+QgiWeXii2HOnNSes2vXUIVX4bHHHmPatGnMnj2bjz/+mM6dO3PBBRdQN8XNuakarvoQ8Bqwr5mVmdm5Zna+mZ0f7fI0sAxYSlhr9jcAUafzdcBb0WPEto5oSZI7TJwI++wDd9wBAwfCsmVwyy3JJ4XttWoVmpXefz80Sd1/f7juzTeHvgwRSYmbb76Z0aNHU69ePZo3b069evXYmo5/Y+6ec4+DDz7YZQdWrHA/7jh3cD/iCPe5czN7/Xfece/d+9vrL1uW2euLpNjChQvjDsE3btzou+22239fr1y50jt27Fjp/hXFDJR6Nf7Gaq6kfPP3v4f+gBkzQnUwfXroN8ikffeFZ56Bv/0tlNxdusBDGoEskoyFCxfyxRdfsGzZMrZu3cqQIUO46KKL0nItJYZ8sWULXHUVnHACtG0b/iBfeOG3nciZZhb6GubPhwMPhNNPh4sugo0b44lHJMfNnj2bM844g9NOO40uXbrQunVrBg4cmJZr5cxwVdmBL7+E/v3Dt/Rzz4Vbb82em9Batw5Vy5VXhs61OXPgiSegadO4IxPJKXPmzKFv3770798/7ddSxZDrVqyAH/0InnsudDJPmJA9SWGbevVgzBh44AF4441w38S//x13VCI5Zc6cOXTt2jUj11LFkMsWL4ZeveDzz8MdzH12NF1VFjj99FBB9OsH//M/IZl16xZ3VCI5YcaMGRm7liqGXDVvHhxxRGiznzkz+5PCNj/8Ibz2Guy8Mxx5ZLiLWkSyihJDLpo1C3r2DE00M2eGm2NyyT77hLibNQsVz8yZcUckIgmUGHLN/PlhArtGjeDll8PQ0FzUpk2Iv1Ur+MlPQt+DiGQFJYZcsq1PobgYXnwx9yeta948zLu0xx6hKWz27LgjEhGUGHLHihVwzDHh+QsvpH5ai7jstVf4PA0bQu/esHRp3BGJFDwlhlzw+edw3HHw2Wfwj3/AfvvFHVFqtWkDzz8f5lXq0wdWr447IpGCpsSQ7TZuDNNdL1oEjz2Wv8M799knTOexYgX07Qtffx13RCIFS4khm7nDb34T+hMmTvy2KSlfHXYYTJ4MpaUwYIBmZhWJiRJDNrvpppAQrroKzjor7mgyo1+/MIX3Y4/BtdfGHY1IVlmyZAk9e/akU6dOXHbZZey9995puY7ufM5Wzz4Ll10GJ50EI0bEHU1mXXopLFgQPnfHjvDzn8cdkch/xbVOz5YtWzjrrLO47bbbOOigg/jd735Hx44dUxtIRIkhG733Hpx2Wpg++95745shNS5mYd6nxYvDpIAHHpi792uIpMj//d//ccABB3DQQQcBsP/++9O4ceO0XCslicHM+gA3AUXABHcftd37Y4Ajo5c7A99z98bRe1uA+dF7/3H3E1IRU85avx5OPjm0rz/+OOy6a9wRxWOnneDhh0Nn+ymnhKkzdt457qhEqrMCZ1rMnj273CR6c+fOpVevXmm5VtJfRc2sCLgNOA44ADjNzA5I3MfdL3H3ru7eFbgFeDzh7XXb3iv4pABhDYV//Qvuuw/S1H6YM1q2DMuEvv02/O53cUcjEqumTZvyzjvvAPDGG29w7733cuCBB6blWqloozgEWOruy9x9IzAZ6LeD/U8DtJxXRSZPDp3NQ4bA8cfHHU126N07dL7ffTdMmhR3NCKx+cUvfkFpaSmdO3fm8ccfp2nTprRv3z4t10pFU1ILYHnC6zLg0Ip2NLM2QDvgxYTNxWZWCmwGRrn7/1Vy7EBgIEDr1q1TEHaWef99OP/8MGSz0DqbqzJ8OLzyClxwAXTvHjqkRQpMs2bNeCOaU2z58uXMmDGDOmnqf8x0r+apwKPuviVhWxt37w6cDow1swrbT9x9nLt3d/fuJSUlmYg1czZvhjPPDP0KDzwAdTUmoJyiInjwwTBx4CmnwFdfxR2RSKzmzp1LlzSu5Z6KxLACaJXwumW0rSKnsl0zkruviH4uA2YAeXpr7w6MHBm+Ed9xR+5PjJcue+4JDz0URipdcknc0YjEqm/fvowfPz5t509FYngL6GBm7cysPuGP/9TtdzKz/YAmwGsJ25qY2U7R82bA4cDCFMSUO159NdzIdeaZcMYZcUeT3Y48MqwdPWECPP103NGI5K2kE4O7bwYuBJ4FFgFT3H2BmY0ws8RRRqcCk93dE7btD5Sa2VxgOqGPoXASw+efh2TQpg3cdlvc0eSG4cOhc2c47zz45JO4oxHJSylpzHb3p4Gnt9t2zXavh1dw3KtA51TEkJMuvBCWLw8L1jRqFHc0uWGnncJNfz16hN/fgw/GHZEUCHfHzOIOo1rKf/+uuQK7pTaLPPVUGKN/1VVhJJJUX9euMGxY6HN45JG4o5ECUFxczNq1a5P+g5sJ7s7atWspLi6u9TksFz7o9rp37+6lpaVxh1F7n38ehlw2aRLWb65fP+6Ics/mzfCDH8CyZeEGuD33jDsiyWObNm2irKyM9evXxx1KtRQXF9OyZUvq1atXbruZzYpGge6QxkXG4YorYNUqeOIJJYXaqls3NCl16wYDB8KTT4Y5lkTSoF69erQroBGDakrKtBdfhPHjwwyiPXrEHU1u228/uPHGsMCP+hpEUkZNSZn09ddhRE1REcydq0nhUmHLFjj88NCk9M47sPvucUckkrWq25SkiiGTrr46TKk9caKSQqoUFcG4cWHo6pVXxh2NSF5QYsiU118PK7L95jdwxBFxR5NfunQJixpNnAgzZ8YdjUjOU1NSJmzeDAcfHL7VLlwIDRvGHVH++frrsLBRcXFYXmunneKOSCTrqCkpm9x2G8ybFyoGJYX02GUXuP320M8wenTc0YjkNCWGdFu1CoYOhT59wvrNkj7HHQf9+8MNN8C778YdjUjOUmJItyuugA0b4JZbNM4+E8aOhQYNwtoWOdhMKpINlBjSacaMsL7CoEGQppWWZDt77gmjRsH06TBlStzRiOQkdT6ny6ZNYU6fdetgwYLwLVYyY8uWcPPgmjWhz2GXXeKOSCQrqPM5bjfdFEYg3XyzkkKmFRWF33tZWageRKRGlBjSoawsrBtw/PHQt2/c0RSmH/4QTj8d/vSncFOhiFRbShKDmfUxs8VmttTMBlfw/tlmtsbM5kSP8xLeG2BmS6LHgFTEE7shQ8K9CzfdFHckhe2PfwyT7V12WdyRiOSUpBODmRUBtwHHAQcAp5nZARXs+rC7d40eE6JjdweGAYcChwDDzKxJsjHF6s03wzoLl12m9Zvj1qIF/OEPYRbbadPijkYkZ6SiYjgEWOruy9x9IzAZ6FfNY3sD09z9E3f/FJgG9ElBTPFwDwvV77knDP5O4SRxuPRS+P734fe/DwMCRKRKqUgMLYDlCa/Lom3b+5mZzTOzR82sVQ2PzQ1TpsCrr8L11+sO52xRXAxjxsCiRVpXW6SaMtX5/Hegrbt3IVQFk2p6AjMbaGalZla6Zs2alAeYtHXrwuyeXbvC2WfHHY0kOv54OPZYuPbaMF+ViOxQKhLDCqBVwuuW0bb/cve17r4hejkBOLi6xyacY5y7d3f37iUlJSkIO8XGjIH//Cf8LCqKOxpJZAZ/+Qt88UWo5kRkh1KRGN4COphZOzOrD5wKTE3cwcyaJ7w8AVgUPX8WONbMmkSdzsdG23LLqlVhJbETT4SePeOORirSqROccw7ceiv8+99xRyOS1ZJODO6+GbiQ8Ad9ETDF3ReY2QgzOyHa7SIzW2Bmc4GLgLOjYz8BriMkl7eAEdG23HL11WE+pD/9Ke5IZEdGjIB69cJIJRGplKbESNa8eaFf4ZJLQnOFZLfhw0Nfw2uvwf/8T9zRiGSUpsTIlMGDoXHjUDVI9rv88jCc+PLLNfuqSCWUGJIxfTo880xommiS2/flFYxddw1NSq+8Em58E5HvUFNSbbnDIYfARx+FRWGKi+ONR6pv82Y48EDYuDHMfFu/ftwRiWSEmpLS7ZFHoLQUrrtOSSHX1K0bBgosXQp33RV3NCJZRxVDbWzaBPvvDzvvDLNn676FXOQORx8Nb78dhq/qTnUpAKoY0mn8+PDHZNQoJYVcZRbuPVmzJiwHKiL/pcRQU19+GYY7/vjHYfF5yV2HHgonnRSalT7+OO5oRLKGEkNN/fWvsHp1mOvfLO5oJFnXXw9ff62V3kQSKDHUxJo18Oc/w8knhxFJkvsOOADOOitMlbF8edX7ixQAJYaa+OMf4ZtvwkgkyR/Dh4fO6BEj4o5EJCsoMVTXypXhW+UvfgH77Rd3NJJKbdrA+efD3XfD4sVxRyMSOyWG6ho5MtwYdc01cUci6XDVVdCgAQwdGnckIrFTYqiODz6AcePg3HPDMpGSf773vbAM6COPwKxZcUcjEislhuq47jqoU0cT5eW7yy+H3XdXVSgFT4mhKkuWwD33wAUXQMuWcUcj6dSoEVxxBTz9dJiWW6RAKTFUZfhw2GmnML225L8LL4SSEhg2LO5IRGKTksRgZn3MbLGZLTWz7/wFNbNLzWyhmc0zsxfMrE3Ce1vMbE70mLr9sbF6+2146CG46CLYY4+4o5FM2HVXGDQIpk2Dl1+OOxqRWCQ9iZ6ZFQHvAscAZYQlOk9z94UJ+xwJvOHu35jZBUBPd+8fvfeVu+9ak2tmbBK9n/0Mnn8e3nsvtD1LYfjmG9h77zAsefr0uKMRSZlMTqJ3CLDU3Ze5+0ZgMtAvcQd3n+7u30QvXweyv7F+1ix4/PEwUkVJobDsvDMMGQIzZsCLL8YdjUjGpSIxtAAS5xIoi7ZV5lzgmYTXxWZWamavm9mJlR1kZgOj/UrXrFmTXMTVMXRoSAiXXJL+a0n2GTgQWrQII5RycGp6kWRktPPZzM4EugN/StjcJiptTgfGmtneFR3r7uPcvbu7dy8pKUlvoK+8EpbsHDQojFSRwlNcHG56e+UVeO65uKMRyahUJIYVQKuE1y2jbeWYWS/gKuAEd9+wbbu7r4h+LgNmAN1SEFNyhg4Nnc2//W3ckUicfvlLaN1aVYMUnFQkhreADmbWzszqA6cC5UYXmVk34C5CUlidsL2Jme0UPW8GHA4sJE4zZoQOxyFDYJddYg1FYrbTTuFLwptvhnsbRApESpb2NLP/BcYCRcDd7n6DmY0ASt19qpk9D3QGVkWH/MfdTzCzHxASxlZCkhrr7hOrul5aRyX17AnvvgvLlmktZwnLuO63HzRuHNb41hocksOqOyqpbiou5u5PA09vt+2ahOe9KjnuVULCyA4zZsBLL8HNNyspSFCvXmhKOvtsePJJOLHS8REieSMlFUOmpa1iOPLIMO2yqgVJtHkzdOwYmpbmzAnzZonkoEzex5AfXnopVAyDByspSHl164YpMubPh0cfjTsakbRTxbDNUUfBokWhWmjQILXnlty3ZQt07hz6GObPV9UgOUkVQ03MnBlGIg0erKQgFSsqClXDwoWqGiTvqWIAOPro8A9e1YLsyJYt0KVLeK6qQXKQKobqevnlMB/OoEFKCrJjRUVhhJKqBslzqhh69YIFC1QtSPWoapAcpoqhOl5+GV54Aa68UklBqkdVgxSAwq4YevUKi/EsWxamWhapDlUNkqNUMVTln//8tlpQUpCaUNUgea5wK4ZjjoF588LqbEoMUlOqGiQHqWLYkVdeCUt2qlqQ2lLVIHmsMCuGY4+FuXND34Km1pbaUtUgOUYVQ2VefRWmTQvVgpKCJENVg+SpwqsYeveG2bND34ISgyRLVYPkEFUMFXnttbB+r6oFSRVVDZKHUpIYzKyPmS02s6VmNriC93cys4ej998ws7YJ7w2Jti82s96piKdS114LJSVwwQVpvYwUmJNPhgMOCP9/bd0adzQiSUs6MZhZEXAbcBxwAHCamR2w3W7nAp+6e3tgDDA6OvYAwhrRHYE+wO3R+dJj2DC44w5VC5JaiVXDI4/EHY1I0lKxtOchwFJ3XwZgZpOBfsDChH36AcOj548Ct5qZRdsnu/sG4D0zWxqd77UUxPUdFz98GHPmALek4+xS0PznsHM7+CVwu2ttaEmLrl1h7Nj0XycVTUktgOUJr8uibRXu4+6bgc+BptU8FgAzG2hmpWZWumbNmhSELZJCZtCmLXzzDej/T8lxqagYMsLdxwHjIIxKqs05MpFppYBtaQpdfhuevzAvNDGJpMq6dWEJYu+d9oo0FRXDCqBVwuuW0bYK9zGzusBuwNpqHiuSGzRCSdJp/Hg47jh48820XyoVieEtoIOZtTOz+oTO5Knb7TMVGBA9Pxl40cMNFFOBU6NRS+2ADkD6P7VIumwboTRiRLjHQSQV1q+H0aPhiCPg0EPTfrmkE0PUZ3Ah8CywCJji7gvMbISZnRDtNhFoGnUuXwoMjo5dAEwhdFT/A/itu+tfk+QuVQ2SDhMmwMqVMHx4Ri5XeHc+i6Rb4t3Q89TXIElavx723js8Xnopqf4F3fksEhdVDZJKEyeGamHYsIwNg1bFIJIOqhokFTZsgPbtoU2bsBRxkolBFYNInIqKwjc8VQ2SjLvvhrKyjFYLoIpBJH22boXOncNzVQ1SUxs2QIcO0KpVWIo4BYlBFYNI3OrUUdUgtXfPPbB8ecarBVDFIJJeqhqkNjZuDNXCXnuFxcVSlBhUMYhkA1UNUhv33AP/+U8s1QKoYhBJP1UNUhMbN8I++8Aee8Drr6c0MahiEMkWqhqkJu69Fz74ILZqAVQxiGSGqgapjk2bQrVQUgJvvJHyxKCKQSSbqGqQ6rj3Xnj//VirBVDFIJI5qhpkRzZuhH33hWbNwtTaaUgMqhhEso2qBtmRv/0tVAsjRsS+NKwqBpFMUtUgFdk2J1KrVvDKK2lLDKoYRLKRqgapyPjxYU6kLKgWQBWDSOapapBE69aFtRY6dIAZM9KaGDJSMZjZ7mY2zcyWRD+bVLBPVzN7zcwWmNk8M+uf8N49Zvaemc2JHl2TiUckJ6hqkER33QWrVmVNtQBJVgxm9kfgE3cfZWaDgSbuPmi7ffYB3N2XmNlewCxgf3f/zMzuAZ5y9xr961DFIDlPVYMAfP01fP/74f+F559P++Uy1cfQD5gUPZ8EnLj9Du7+rrsviZ6vBFYDJUleVyS3qWoQgNtvh9WrQ7WQRZKtGD5z98bRcwM+3fa6kv0PISSQju6+NaoYDgM2AC8Ag919QyXHDgQGArRu3frgDz74oNZxi2QFVQ2F7csvoV076NEDnnkmI5dMWcVgZs+b2dsVPPol7uchw1SaZcysOXAfcI67b402DwH2A3oAuwODKjkcdx/n7t3dvXtJiQoOyQOqGgrbLbfA2rVw7bVxR/IdyVYMi4Ge7r4q+sM/w933rWC/RsAMYGRl/Qlm1hO43N37VnVd9TFI3lDVUJg+/zxUC4cfDn//e8Yum6k+hqnAgOj5AODJCgKpDzwB3Lt9UoiSybZmqBOBt5OMRyS3qGooTGPHwqefZl3fwjbJVgxNgSlAa+AD4Ofu/omZdQfOd/fzzOxM4G/AgoRDz3b3OWb2IqEj2oA50TFfVXVdVQySV1Q1FJZPP4W2beHoo+HxxzN66epWDHWTuYi7rwWOrmB7KXBe9Px+4P5Kjj8qmeuL5IVtVUP//vDII3DqqXFHJOn0l7/AF19kZd/CNpoSQyQbnHwydOoE11wT5uSX/PTRR6EZqX//b6vELKTEIJIN6tSBG26AJUvCer+Sn0aOhPXr4brr4o5kh5QYRLLF8cfDYYeFJoZ16+KORlLt/ffhjjvg3HPDvEhZTIlBJFuYwahRsGIF3HZb3NFIqg0bFgYWXHNN3JFUSYlBJJsccQT06QM33hjGukt+ePttuO8++N3voEWLuKOpkhKDSLYZORI++QT+/Oe4I5FUufpqaNgQBlU6uUNWUWIQyTbduoVRK2PGhFEskttefx2efBKuvBKaNo07mmpRYhDJRtddF0av3HBD3JFIMtxh8GD43vfg97+PO5pqU2IQyUYdOsB558Gdd8J778UdjdTWc8/BSy/B0KGw665xR1NtSgwi2eqaa6Bu3dA+Lblny5bQfNS2LQwcGHc0NaLEIJKt9toLLrsMHnwQ3nor7mikpu67L8x9NWoU1K8fdzQ1ktQkenHRJHpSML78Etq3h/32S/tC8ZJC33wTmgNbtgydz1ny3y1T026LSDo1bBjuhJ45M6Pz9kuSxoyBlSvDkOMsSQo1oYpBJNtt3hwmXHOH+fOhXr24I5Id+eijUOX16gVPPBF3NOWoYhDJF3Xrwh//CIsXw4QJcUcjVbn22jDUePTouCOptaQSg5ntbmbTzGxJ9LNJJfttMbM50WNqwvZ2ZvaGmS01s4ej1d5EZHt9+8KPfxzm2/nii7ijkcq88w6MGwe//jXss0/c0dRashXDYOAFd+8AvBC9rsg6d+8aPU5I2D4aGOPu7YFPgXOTjEckP5mF9uo1a0L1INlp0CDYeeeQwHNYsomhHzApej6JsG5ztUTrPB8FbFvotkbHixSc7t3hjDPCCmD/+U/c0cj2pk+HqVPhD3+AkpK4o0lKsolhD3dfFT3/ENijkv2KzazUzF43s21//JsCn7n75uh1GZD90w6KxGnkyFA9XHll3JFIos2bw5QXbdrk1NQXlalyzWczex7Ys4K3rkp84e5uZpUNcWrj7ivM7PvAi2Y2H6jRnMJmNhAYCNC6deuaHCqSP1q3Ds0Vw4fDBReEfgeJ3113hRFjjz0GDRrEHU3SkhquamaLgZ7uvsrMmgMz3H3fKo65B3gKeAxYA+zp7pvN7DBguLv3ruq6Gq4qBW3dOth/f9htN5g1K4wRM5aKAAANVklEQVRakvisXRtuZuvWDZ5/PqvvW8jUcNWpwIDo+QDgyQoCaWJmO0XPmwGHAws9ZKTpwMk7Ol5EttOgQehnmDcPxo+POxoZOjSMFLvppqxOCjWRbGIYBRxjZkuAXtFrzKy7mW0bcL0/UGpmcwmJYJS7L4zeGwRcamZLCX0OE5OMR6Qw/PSncOSRYYK9tWvjjqZwzZ0bmpF++1vo1CnuaFJGdz6L5Kr580Pzxa9/rTWi4+AOPXvCwoXw7rvQpMLbuLKK7nwWyXedO8NvfhPWbJg7N+5oCs+UKWEOq5EjcyIp1IQqBpFc9umnoeOzY0fNvppJX38dZrwtKQlTohcVxR1RtahiECkETZrAjTeGb66TJlW9v6TGsGFQVga33pozSaEmlBhEct2558Lhh8Pll8PHH8cdTf6bPRvGjg19Oz/4QdzRpIUSg0iuq1MnjIz54ouQHCR9tmwJy3Q2axYqtTylxCCSDzp2DNNkTJoEL74YdzT56/bbobQ0VAx51uGcSJ3PIvli3bowUqlOnXDzW3Fx3BHll7KycMf5D38ITz+dkx396nwWKTQNGoShq0uWhCGUkloXXRSakm6/PSeTQk0oMYjkk1694MwzYdQoWLQo7mjyx5NPhmU6hw2Ddu3ijibtlBhE8s1f/gING8J554VvuJKczz4LU1507gyXXhp3NBmhxCCSb773vTCh26uvwl//Gnc0ue+ii+DDD2HiRKhXL+5oMkKJQSQfnXFGmGjv6qvh7bfjjiZ3PfEE3HcfXHUV9OgRdzQZo8Qgko/MQkf0brvBWWfBpk1xR5R7Vq8ON7F16xYSQwFRYhDJVyUlMG5cuFP3+uvjjia3uIcV8j7/HO69F+rXjzuijFJiEMlnJ54Iv/gF3HBDmOxNqueBB+Dxx+G66/JqnYXq0g1uIvnus8/CiJpdd4V//Ssv1iROq7KykAw6dYKXXsqrSfIycoObme1uZtPMbEn08zv3iJvZkWY2J+Gx3sxOjN67x8zeS3ivazLxiEgFGjcOI2reeSdMmyGV27IFzj479Mncc09eJYWaSLYpaTDwgrt3AF6IXpfj7tPdvau7dwWOAr4BnkvY5Ypt77v7nCTjEZGKHHssXHJJmCb60UfjjiZ7XXcdvPAC3HwztG8fdzSxSTYx9AO2TQI/CTixiv1PBp5x92+SvK6I1NSoUXDoofDLX8LSpXFHk32efx5GjAijuH75y7ijiVWyiWEPd18VPf8Q2KOK/U8FHtpu2w1mNs/MxpjZTpUdaGYDzazUzErXrFmTRMgiBap+/bAcZb16cMopsH593BFlj5Ur4fTTwyR5BTAXUlWqTAxm9ryZvV3Bo1/ifh56sSvtyTaz5kBn4NmEzUOA/YAewO7AoMqOd/dx7t7d3buXlJRUFbaIVKR16zD8cs4cuPjiuKPJDps3w6mnhuU6H30Udtkl7ohiV7eqHdy9V2XvmdlHZtbc3VdFf/hX7+BUPweecPf/3mmTUG1sMLO/AVplRCTdfvITGDQIRo+GI44I35QL2TXXwMsvhzuc998/7miyQrJNSVOBAdHzAcCTO9j3NLZrRoqSCWZmhP4J3bsvkgnXXx/WFRg4EBYujDua+Dz1VFiJ7Ve/CrPSCpB8YhgFHGNmS4Be0WvMrLuZTdi2k5m1BVoBL213/ANmNh+YDzQDdHumSCbUrQuTJ4d7G044AdaujTuizFu4MFRL3bqFSQflv3SDm0ghe+016NkzLGr/3HMFM3soa9fCIYeEfoW33oJWreKOKCO0gpuIVO2ww8LNbzNmhDUHcvCLYo1t2AA/+xmsWBEW4CmQpFATVXY+i0ieO/PMsNrbyJFh1NLVV8cdUfps3RrubH7ppTAf0qGHxh1RVlJiEJHQGV1WBkOHQosWcM45cUeUHldeGfpWRo/WaKwdUGIQkXBD1/jxsGpVGKHTpEmYmTWfjB4dlj298EK44oq4o8lq6mMQkaB+fXjsMejeHfr3h3/8I+6IUueWW2DwYDjtNBg7tuDvbK6KEoOIfKthw5AQOnaEk04KE8rlunHjwrrNJ50EkyYV7IypNaHEICLlNW4chq62bx/ukv5//y/uiGrvppvC8pz/+7/w0EOFMxw3SUoMIvJdzZqFIaydOoVv2rk2Vbd7GGV18cXw05+G1dh2qnSOTtmOEoOIVKxp09CUdMgh8POf587dwZs3hw7mq66CM86Ahx9WUqghJQYRqdxuu4VmpRNPDN++L7oorHKWrb78Evr1C1NnX3llmEm2rgZf1pQSg4js2M47wyOPwKWXhtE9vXtDNq6J8s474Ya1Z5+FO+8Mw1Pr6E9cbei3JiJVKyoK9wDcfTe88gocdBC8+mrcUX1ryhTo0QM+/jgkhl//Ou6IcpoSg4hU3znnhIRQrx786Efwhz+EuYfi8umnoR+hf//QUf6vf8HRR8cXT55QYhCRmunWLawAd845YS2DHj3gn//MbAzuoXmrU6dQLVx7LcycCS1bZjaOPKXEICI116gRTJgAf/87fPZZqB7OOAM++CD91549G446KoyUKimB118Pq7DpHoWUUWIQkdrr2zfMzDp0aJhOo337MNfSsmWpv1ZpaRhxdNBBMG9eGHk0axYcfHDqr1XgkkoMZnaKmS0ws61mVuniD2bWx8wWm9lSMxucsL2dmb0RbX/YzOonE4+IxGCXXWDECFiyJHT63ntvSBC9e4fmnm++qf25P/kkTGnRo0d4zJwZmo2WLoULLtD0FmmS1ApuZrY/sBW4C7jc3b+zrJqZFQHvAscAZcBbwGnuvtDMpgCPu/tkM7sTmOvud1R1Xa3gJpLFVq4MM7VOnAjLl0NxcWj6OfLI8O3+wAPD7K3bT2TnDh9+GJqKZs2CadPCCKitW0Nfwq9+BQMGhHsrpFaqu4JbSpb2NLMZVJ4YDgOGu3vv6PWQ6K1RwBpgT3ffvP1+O6LEIJIDtmwJ02o89VR4LF367XsNGkDz5iFpAHz1VUgKGzd+u0/XrnD88WFN6oMP1oyoKVDdxJCJWwJbAMsTXpcBhwJNgc/cfXPC9haVncTMBgIDAVq3bp2eSEUkdYqKwtDRo4+GMWPCTXGzZsGCBaGqSEwEDRrAXnuFRYK6dAkjnxo1ijf+AlZlYjCz54E9K3jrKnd/MvUhVczdxwHjIFQMmbquiKRISQn06RMektWqTAzu3ivJa6wAElfbbhltWws0NrO6UdWwbbuIiMQoE8NV3wI6RCOQ6gOnAlM9dG5MB06O9hsAZKwCERGRiiU7XPUkMysDDgP+n5k9G23fy8yeBoiqgQuBZ4FFwBR3XxCdYhBwqZktJfQ5TEwmHhERSV5KRiVlmkYliYjUXHVHJenOZxERKUeJQUREylFiEBGRcpQYRESknJzsfDazNUBt5vdtBnyc4nAyLdc/Q67HD/oM2UKfoebauHtJVTvlZGKoLTMrrU6PfDbL9c+Q6/GDPkO20GdIHzUliYhIOUoMIiJSTqElhnFxB5ACuf4Zcj1+0GfIFvoMaVJQfQwiIlK1QqsYRESkCkoMIiJSTkEkBjPrY2aLzWypmQ2OO57aMLO7zWy1mb0ddyy1YWatzGy6mS00swVm9vu4Y6opMys2szfNbG70Ga6NO6baMLMiM5ttZk/FHUttmdn7ZjbfzOaYWc7NqGlmjc3sUTN7x8wWRUsbZ42872MwsyLgXeAYwvKhbwGnufvCWAOrITM7AvgKuNfdO8UdT02ZWXOgubv/y8waArOAE3Ppv4OZGbCLu39lZvWAfwK/d/fXYw6tRszsUqA70Mjd+8YdT22Y2ftAd3fPyRvczGwS8LK7T4jWqdnZ3T+LO65tCqFiOARY6u7L3H0jMBnoF3NMNebuM4FP4o6jttx9lbv/K3r+JWFtjkrX+M5GHnwVvawXPXLqm5WZtQR+AkyIO5ZCZWa7AUcQrT/j7huzKSlAYSSGFsDyhNdl5NgfpHxjZm2BbsAb8UZSc1EzzBxgNTDN3XPtM4wFrgS2xh1Ikhx4zsxmmdnAuIOpoXbAGuBvUZPeBDPbJe6gEhVCYpAsYma7Ao8BF7v7F3HHU1PuvsXduxLWKD/EzHKmWc/M+gKr3X1W3LGkwA/d/SDgOOC3UVNrrqgLHATc4e7dgK+BrOr7LITEsAJolfC6ZbRNMixql38MeMDdH487nmREpf90oE/csdTA4cAJUfv8ZOAoM7s/3pBqx91XRD9XA08QmoxzRRlQllBtPkpIFFmjEBLDW0AHM2sXdfKcCkyNOaaCE3XcTgQWuftf446nNsysxMwaR88bEAY0vBNvVNXn7kPcvaW7tyX8O3jR3c+MOawaM7NdogEMRE0wxwI5M1rP3T8ElpvZvtGmo4GsGoRRN+4A0s3dN5vZhcCzQBFwt7sviDmsGjOzh4CeQDMzKwOGufvEeKOqkcOBXwDzozZ6gD+4+9MxxlRTzYFJ0Ui3OsAUd8/ZIZ85bA/gifBdg7rAg+7+j3hDqrHfAQ9EX1aXAefEHE85eT9cVUREaqYQmpJERKQGlBhERKQcJQYRESlHiUFERMpRYhARkXKUGEREpBwlBhERKef/AzeQmKkgNlduAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xnc1WP+x/HXp7u7BVHqRtopW0VRjDE/E6L8JmKGyR7DNMxgyFKNpYSmZisMZlKWDJJtGD+GomQQ7rSppGTprpDstPf5/XF9G+fOfXcvZ/me5f18PM6jc77L+X5OdD7nc13X97rM3REREdmiTtwBiIhIdlFiEBGRcpQYRESkHCUGEREpR4lBRETKUWIQEZFylBgkZ5lZazP72syKUvBed5vZDamISyTX1Y07AJGqmNl7wK7ApoTNe7n7B8AOsQSVA8zsbqDM3a+OOxbJLUoMkiuOc/cpcQchUgjUlCQ5y8zampmbWd3o9TQzu97MXjKzr8zsWTNrlnD8Q2b2oZl9YWbTzaxjNa9zdvSeo83sczNbamY/jLYvM7OPzax/wvE7mdkEM1tlZu+b2dVmVqeW71XfzP5kZh+Y2Udm9jczaxjt62FmZWZ2WXTeSjM7J9o3ADgduDJqbvtXtH2QmS2P/n4WmdlRyf+XkHyjxCD55jTgHGAXoB5wecK+p4EO0b43gPtq8L6HAHOBpsD9wESgO9AeOAP4q5ltada6BdgJ2AP4MXBWFFNt3msksBfQJdrfArg24b12i67VAjgXuNXMmrj72Ojz/cHdd3D348xsb+BCoLu7NwJ6Ae/V4O9ACoQSg+SKf0a/sD83s39u47i73P1td18DTCJ8oQLg7ne6+1fuvg4YBhxgZjtV8/rvuvtd7r4JeBBoBQx393Xu/iywHmgfdYSfAgyJrvUe8GfgzFq8lwEDgEvd/VN3/woYEb3/Fhuicze4+1PA18DelXyGTUB9YD8zK3b399z9nWp+fikg6mOQXHFCNfsYPkx4/i1R53T0hX0jcDJQAmyOjmkGfFGN9/0o4fkaAHffetsO0fsVA+8n7Huf8Iu+pu9VAmwHzAw5AgADEkdhrXb3jQmv//uZt+buS8zsEkJS7GhmzwAD3X1FRcdL4VLFIIXiNKAv0JPQ9NI22m6VnVBLnxB+xbdJ2NYaWF7L91oDdHT3xtFjJ3ev7kis702d7O73u/uPovgcGFWLuCTPKTFIoWgErANWE36Fj0jHRaLmoUnAjWbWyMzaAAOBf9TivTYDdwCjzWwXADNrYWa9qvkWHxH6OYjO3dvMjjSz+sBaQtLZXNnJUriUGKRQTCA06SwHFgAz0niti4BvgKXAfwgdzHfW8r0GAUuAGWb2JTCFyvsQtjae0J+wpV+mPqEz+xNCk9suwJBaxiV5zLRQj4iIJFLFICIi5SgxiIhIOUoMIiJSjhKDiIiUk5M3uDVr1szbtm0bdxgiIjll5syZn7h7SVXH5WRiaNu2LaWlpXGHISKSU8zs/aqPUlOSiIhsRYlBRETKUWIQEZFycrKPQUQkkzZs2EBZWRlr166NO5RqadCgAS1btqS4uLhW5ysxiIhUoaysjEaNGtG2bVsSpkDPSu7O6tWrKSsro127drV6j5Q0JZnZndHSgm9Wst/M7GYzW2Jmc83swIR9/c1scfToX9H5IiJxWrt2LU2bNs36pABgZjRt2jSp6iZVfQx3A723sf9YwpKKHQgrUt0OYGY7A0MJSx0eDAw1syYpiklEJGVyISlskWysKWlKcvfpZtZ2G4f0BSZ4mMp1hpk1NrPmQA9gsrt/CmBmkwkJ5oFUxCU1sG4dzJ8P77wDK1fC6tWwZebdnXeG5s2hXTvYf39o0CDeWEUkrTLVx9ACWJbwuizaVtn27zGzAYRqg9atW6cnykKyeTO89ho8+SQ88wzMmQMbNpQ/xuy75LBFURF06gTHHAPHHQeHHgp11VUlkk9y5l+0u48FxgJ069ZNi0jUVlkZ3HlneLz/fvii/+EPYeBAOOgg2Htv2H33UCXUqRMSw2efwYoV8PbbMHMmzJgBY8bAH/8Yjj3nHDj33FBRiEjOy9R9DMuBVgmvW0bbKtsuqfbOO3DeeeHLe+hQ6NAB/vEPWLUKpk+HkSPh5JNDU1GzZiEpQKgadt45VAk//SnceCM89xx88glMmgRdu8Lvfw/t28MZZ8DChfF+TpE8d9ttt9GpUyfatGnDLbfckpZrZKpieAK40MwmEjqav3D3lWb2DDAiocP5GLTUYGp9+ilcey387W+hyef88+HSS2GPPao+d1t23DEkkpNPhmXL4Oab4bbb4P774eyzQ7LYddeUfASRrHLJJTB7dmrfs0uXUIVX4ZFHHmHy5MnMmjWLTz75hM6dO3PBBRdQN8XNuakarvoA8Aqwt5mVmdm5Zna+mZ0fHfIUYf3bJYTFzX8NEHU6Xw+8Hj2Gb+mIliS5w/jxsNdecPvtMGAALF0Kt9ySfFLYWqtWoVnpvfdCk9Q//hGue/PNoS9DRFLi5ptvZtSoURQXF9O8eXOKi4vZnI5/Y+6ec4+DDjrIZRuWL3c/9lh3cD/8cPc5czJ7/bfecu/V67vrL12a2euLpNiCBQviDsHXr1/vO+20039fr1ixwjt27Fjp8RXFDJR6Nb5jNVdSvvnXv0J/wLRpoTqYOjX0G2TS3nvD00/DXXeFknv//eEBjUAWScaCBQv48ssvWbp0KZs3b2bIkCFcfPHFabmWEkO+2LQJrroKjj8e2rYNX8gXXvhdJ3KmmYW+hnnz4IAD4LTT4OKLYf36eOIRyXGzZs3i9NNP59RTT2X//fendevWDBgwIC3XypnhqrINX30F/fqFX+nnngt//Wv23ITWunWoWq68MnSuzZ4Njz0GTZvGHZlITpk9ezZ9+vShX79+ab+WKoZct3w5/M//wLPPhk7mceOyJylsUVwMo0fDfffBq6+G+ybeeSfuqERyyuzZs+nSpUtGrqWKIZctWgQ9e8IXX4Q7mHtva7qqLHDaaaGC6NsXfvCDkMy6do07KpGcMG3atIxdSxVDrpo7Fw4/PLTZT5+e/Ulhix/9CF55BbbbDo44ItxFLSJZRYkhF82cCT16hCaa6dPDzTG5ZK+9QtzNmoWKZ/r0uCMSkQRKDLlm3rwwgd2OO8KLL4ahobmoTZsQf6tW8JOfhL4HEckKSgy5ZEufQoMG8PzzuT9pXfPmYd6lXXcNTWGzZsUdkYigxJA7li+Ho48Oz597LvXTWsRl993D52nUCHr1giVL4o5IpOApMeSCL76AY4+Fzz+Hf/8b9tkn7ohSq00bmDIlzKvUuzd8/HHcEYkUNCWGbLd+fZjueuFCeOSR/B3euddeYTqP5cuhTx/45pu4IxIpWEoM2cwdfv3r0J8wfvx3TUn56tBDYeJEKC2F/v01M6tITJQYstlNN4WEcNVVcNZZcUeTGX37him8H3kErrsu7mhEssrixYvp0aMHnTp14rLLLmPPPfdMy3V053O2euYZuOwyOPFEGD487mgya+BAmD8/fO6OHeHnP487IpH/imudnk2bNnHWWWdx6623cuCBB3LRRRfRsWPH1AYSUWLIRu++C6eeGqbPnjAhvhlS42IW5n1atChMCnjAAbl7v4ZIivzzn/9kv/3248ADDwRg3333pXHjxmm5VkoSg5n1Bm4CioBx7j5yq/2jgSOil9sBu7h742jfJmBetO8Ddz8+FTHlrLVr4aSTQvv6o4/CDjvEHVE86teHBx8Mne0nnxymzthuu7ijEqnOCpxpMWvWrHKT6M2ZM4eePXum5VpJ/xQ1syLgVuBYYD/gVDPbL/EYd7/U3bu4exfgFuDRhN1rtuwr+KQAYQ2FN96Ae++FNLUf5oyWLcMyoW++CRddFHc0IrFq2rQpb731FgCvvvoqEyZM4IADDkjLtVLRRnEwsMTdl7r7emAi0Hcbx58KaDmvikycGDqbhwyB446LO5rs0KtX6Hy/80645564oxGJzZlnnklpaSmdO3fm0UcfpWnTprRv3z4t10pFU1ILYFnC6zLgkIoONLM2QDvg+YTNDcysFNgIjHT3f1Zy7gBgAEDr1q1TEHaWee89OP/8MGSz0DqbqzJsGLz0ElxwAXTrFjqkRQpMs2bNeDWaU2zZsmVMmzaNOmnqf8x0r+YpwMPuvilhWxt37wacBowxswrbT9x9rLt3c/duJSUlmYg1czZuhDPOCP0K990HdTUmoJyiIrj//jBx4Mknw9dfxx2RSKzmzJnD/mlcyz0ViWE50CrhdctoW0VOYatmJHdfHv25FJgG5OmtvdswYkT4RXz77bk/MV667LYbPPBAGKl06aVxRyMSqz59+nDHHXek7f1TkRheBzqYWTszq0f48n9i64PMbB+gCfBKwrYmZlY/et4MOAxYkIKYcsfLL4cbuc44A04/Pe5ostsRR4S1o8eNg6eeijsakbyVdGJw943AhcAzwEJgkrvPN7PhZpY4yugUYKK7e8K2fYFSM5sDTCX0MRROYvjii5AM2rSBW2+NO5rcMGwYdO4M550Hn34adzQieSkljdnu/hTw1Fbbrt3q9bAKznsZ6JyKGHLShRfCsmVhwZodd4w7mtxQv3646a979/D3d//9cUckBcLdMbO4w6iW8r+/a67AbqnNIk8+GcboX3VVGIkk1delCwwdGvocHnoo7mikADRo0IDVq1cn/YWbCe7O6tWradCgQa3fw3Lhg26tW7duXlpaGncYtffFF2HIZZMmYf3mevXijij3bNwIP/whLF0aboDbbbe4I5I8tmHDBsrKyli7dm3coVRLgwYNaNmyJcXFxeW2m9nMaBToNmlcZByuuAJWroTHHlNSqK26dUOTUteuMGAAPP54mGNJJA2Ki4tpV0AjBtWUlGnPPw933BFmEO3ePe5octs++8Dvfx8W+FFfg0jKqCkpk775JoyoKSqCOXM0KVwqbNoEhx0WmpTeegt23jnuiESyVnWbklQxZNLVV4cptcePV1JIlaIiGDs2DF298sq4oxHJC0oMmTJjRliR7de/hsMPjzua/LL//mFRo/HjYfr0uKMRyXlqSsqEjRvhoIPCr9oFC6BRo7gjyj/ffBMWNmrQICyvVb9+3BGJZB01JWWTW2+FuXNDxaCkkB7bbw+33Rb6GUaNijsakZymxJBuK1fCNddA795h/WZJn2OPhX794MYb4e23445GJGcpMaTbFVfAunVwyy0aZ58JY8ZAw4ZhbYscbCYVyQZKDOk0bVpYX2HQIEjTSkuyld12g5EjYepUmDQp7mhEcpI6n9Nlw4Ywp8+aNTB/fvgVK5mxaVO4eXDVqtDnsP32cUckkhXU+Ry3m24KI5BuvllJIdOKisLfe1lZqB5EpEaUGNKhrCysG3DccdCnT9zRFKYf/QhOOw3++MdwU6GIVFtKEoOZ9TazRWa2xMwGV7D/bDNbZWazo8d5Cfv6m9ni6NE/FfHEbsiQcO/CTTfFHUlh+8MfwmR7l10WdyQiOSXpxGBmRcCtwLHAfsCpZrZfBYc+6O5dose46NydgaHAIcDBwFAza5JsTLF67bWwzsJll2n95ri1aAG/+12YxXby5LijEckZqagYDgaWuPtSd18PTAT6VvPcXsBkd//U3T8DJgO9UxBTPNzDQvW77QaDv1c4SRwGDoQ99oDf/jYMCBCRKqUiMbQAliW8Lou2be1nZjbXzB42s1Y1PDc3TJoEL78MN9ygO5yzRYMGMHo0LFyodbVFqilTnc//Atq6+/6EquCemr6BmQ0ws1IzK121alXKA0zamjVhds8uXeDss+OORhIddxwccwxcd12Yr0pEtikViWE50Crhdcto23+5+2p3Xxe9HAccVN1zE95jrLt3c/duJSUlKQg7xUaPhg8+CH8WFcUdjSQygz//Gb78MlRzIrJNqUgMrwMdzKydmdUDTgGeSDzAzJonvDweWBg9fwY4xsyaRJ3Ox0TbcsvKlWElsRNOgB494o5GKtKpE5xzDvz1r/DOO3FHI5LVkk4M7r4RuJDwhb4QmOTu881suJkdHx12sZnNN7M5wMXA2dG5nwLXE5LL68DwaFtuufrqMB/SH/8YdySyLcOHQ3FxGKkkIpXSlBjJmjs39CtcemlorpDsNmxY6Gt45RX4wQ/ijkYkozQlRqYMHgyNG4eqQbLf5ZeH4cSXX67ZV0UqocSQjKlT4emnQ9NEk9y+L69g7LBDaFJ66aVw45uIfI+akmrLHQ4+GD76KCwK06BBvPFI9W3cCAccAOvXh5lv69WLOyKRjFBTUro99BCUlsL11ysp5Jq6dcNAgSVL4O9/jzsakayjiqE2NmyAffeF7baDWbN030IucoejjoI33wzDV3WnuhQAVQzpdMcd4ctk5EglhVxlFu49WbUqLAcqIv+lxFBTX30Vhjv++Mdh8XnJXYccAieeGJqVPvkk7mhEsoYSQ0395S/w8cdhrn+zuKORZN1wA3zzjVZ6E0mgxFATq1bBn/4EJ50URiRJ7ttvPzjrrDBVxrJlVR8vUgCUGGriD3+Ab78NI5EkfwwbFjqjhw+POxKRrKDEUF0rVoRflWeeCfvsE3c0kkpt2sD558Odd8KiRXFHIxI7JYbqGjEi3Bh17bVxRyLpcNVV0LAhXHNN3JGIxE6JoTrefx/GjoVzzw3LREr+2WWXsAzoQw/BzJlxRyMSKyWG6rj+eqhTRxPl5bvLL4edd1ZVKAVPiaEqixfD3XfDBRdAy5ZxRyPptOOOcMUV8NRTYVpukQKlxFCVYcOgfv0wvbbkvwsvhJISGDo07khEYpOSxGBmvc1skZktMbPvfYOa2UAzW2Bmc83sOTNrk7Bvk5nNjh5PbH1urN58Ex54AC6+GHbdNe5oJBN22AEGDYLJk+HFF+OORiQWSU+iZ2ZFwNvA0UAZYYnOU919QcIxRwCvuvu3ZnYB0MPd+0X7vnb3HWpyzYxNovezn8GUKfDuu6HtWQrDt9/CnnuGYclTp8YdjUjKZHISvYOBJe6+1N3XAxOBvokHuPtUd/82ejkDyP7G+pkz4dFHw0gVJYXCst12MGQITJsGzz8fdzQiGZeKxNACSJxLoCzaVplzgacTXjcws1Izm2FmJ1R2kpkNiI4rXbVqVXIRV8c114SEcOml6b+WZJ8BA6BFizBCKQenphdJRkY7n83sDKAb8MeEzW2i0uY0YIyZ7VnRue4+1t27uXu3kpKS9Ab60kthyc5Bg8JIFSk8DRqEm95eegmefTbuaEQyKhWJYTnQKuF1y2hbOWbWE7gKON7d123Z7u7Loz+XAtOArimIKTnXXBM6m3/zm7gjkTj94hfQurWqBik4qUgMrwMdzKydmdUDTgHKjS4ys67A3wlJ4eOE7U3MrH70vBlwGLCAOE2bFjochwyB7bePNRSJWf364UfCa6+FextECkRKlvY0s/8FxgBFwJ3ufqOZDQdK3f0JM5sCdAZWRqd84O7Hm9kPCQljMyFJjXH38VVdL62jknr0gLffhqVLtZazhGVc99kHGjcOa3xrDQ7JYdUdlVQ3FRdz96eAp7badm3C856VnPcyIWFkh2nT4IUX4OablRQkKC4OTUlnnw2PPw4nVDo+QiRvpKRiyLS0VQxHHBGmXVa1IIk2boSOHUPT0uzZYd4skRyUyfsY8sMLL4SKYfBgJQUpr27dMEXGvHnw8MNxRyOSdqoYtjjySFi4MFQLDRum9r0l923aBJ07hz6GefNUNUhOUsVQE9Onh5FIgwcrKUjFiopC1bBggaoGyXuqGACOOir8g1e1INuyaRPsv394rqpBcpAqhup68cUwH86gQUoKsm1FRWGEkqoGyXOqGHr2hPnzVS1I9ahqkBymiqE6XnwRnnsOrrxSSUGqR1WDFIDCrhh69gyL8SxdGqZaFqkOVQ2So1QxVOU///muWlBSkJpQ1SB5rnArhqOPhrlzw+psSgxSU6oaJAepYtiWl14KS3aqWpDaUtUgeawwK4ZjjoE5c0LfgqbWltpS1SA5RhVDZV5+GSZPDtWCkoIkQ1WD5KnCqxh69YJZs0LfghKDJEtVg+QQVQwVeeWVsH6vqgVJFVUNkodSkhjMrLeZLTKzJWY2uIL99c3swWj/q2bWNmHfkGj7IjPrlYp4KnXddVBSAhdckNbLSIE56STYb7/w/9fmzXFHI5K0pBODmRUBtwLHAvsBp5rZflsddi7wmbu3B0YDo6Jz9yOsEd0R6A3cFr1fegwdCrffrmpBUiuxanjoobijEUlaKpb2PBhY4u5LAcxsItAXWJBwTF9gWPT8YeCvZmbR9onuvg5418yWRO/3Sgri+p5LHjyU2bOBW9Lx7lLQ/OewXTv4BXCba21oSYsuXWDMmPRfJxVNSS2AZQmvy6JtFR7j7huBL4Cm1TwXADMbYGalZla6atWqFIQtkkJm0KYtfPst6P9PyXGpqBgywt3HAmMhjEqqzXtkItNKAdvUFPb/TXj+3NzQxCSSKmvWhCWIvVfaK9JUVAzLgVYJr1tG2yo8xszqAjsBq6t5rkhu0AglSac77oBjj4XXXkv7pVKRGF4HOphZOzOrR+hMfmKrY54A+kfPTwKe93ADxRPAKdGopXZAByD9n1okXbaMUBo+PNzjIJIKa9fCqFFw+OFwyCFpv1zSiSHqM7gQeAZYCExy9/lmNtzMjo8OGw80jTqXBwKDo3PnA5MIHdX/Bn7j7vrXJLlLVYOkw7hxsGIFDBuWkcsV3p3PIumWeDf0XPU1SJLWroU99wyPF15Iqn9Bdz6LxEVVg6TS+PGhWhg6NGPDoFUxiKSDqgZJhXXroH17aNMmLEWcZGJQxSASp6Ki8AtPVYMk4847oawso9UCqGIQSZ/Nm6Fz5/BcVYPU1Lp10KEDtGoVliJOQWJQxSAStzp1VDVI7d19NyxblvFqAVQxiKSXqgapjfXrQ7Ww++5hcbEUJQZVDCLZQFWD1Mbdd8MHH8RSLYAqBpH0U9UgNbF+Pey1F+y6K8yYkdLEoIpBJFuoapCamDAB3n8/tmoBVDGIZIaqBqmODRtCtVBSAq++mvLEoIpBJJuoapDqmDAB3nsv1moBVDGIZI6qBtmW9eth772hWbMwtXYaEoMqBpFso6pBtuWuu0K1MHx47EvDqmIQySRVDVKRLXMitWoFL72UtsSgikEkG6lqkIrccUeYEykLqgVQxSCSeaoaJNGaNWGthQ4dYNq0tCaGjFQMZrazmU02s8XRn00qOKaLmb1iZvPNbK6Z9UvYd7eZvWtms6NHl2TiEckJqhok0d//DitXZk21AElWDGb2B+BTdx9pZoOBJu4+aKtj9gLc3Reb2e7ATGBfd//czO4GnnT3Gv3rUMUgOU9VgwB88w3ssUf4f2HKlLRfLlN9DH2Be6Ln9wAnbH2Au7/t7ouj5yuAj4GSJK8rkttUNQjAbbfBxx+HaiGLJFsxfO7ujaPnBny25XUlxx9MSCAd3X1zVDEcCqwDngMGu/u6Ss4dAAwAaN269UHvv/9+reMWyQqqGgrbV19Bu3bQvTs8/XRGLpmyisHMppjZmxU8+iYe5yHDVJplzKw5cC9wjrtvjjYPAfYBugM7A4MqOR13H+vu3dy9W0mJCg7JA6oaCtstt8Dq1XDddXFH8j3JVgyLgB7uvjL64p/m7ntXcNyOwDRgRGX9CWbWA7jc3ftUdV31MUjeUNVQmL74IlQLhx0G//pXxi6bqT6GJ4D+0fP+wOMVBFIPeAyYsHVSiJLJlmaoE4A3k4xHJLeoaihMY8bAZ59lXd/CFslWDE2BSUBr4H3g5+7+qZl1A8539/PM7AzgLmB+wqlnu/tsM3ue0BFtwOzonK+ruq4qBskrqhoKy2efQdu2cNRR8OijGb10dSuGuslcxN1XA0dVsL0UOC96/g/gH5Wcf2Qy1xfJC1uqhn794KGH4JRT4o5I0unPf4Yvv8zKvoUtNCWGSDY46STo1AmuvTbMyS/56aOPQjNSv37fVYlZSIlBJBvUqQM33giLF4f1fiU/jRgBa9fC9dfHHck2KTGIZIvjjoNDDw1NDGvWxB2NpNp778Htt8O554Z5kbKYEoNItjCDkSNh+XK49da4o5FUGzo0DCy49tq4I6mSEoNINjn8cOjdG37/+zDWXfLDm2/CvffCRRdBixZxR1MlJQaRbDNiBHz6KfzpT3FHIqly9dXQqBEMqnRyh6yixCCSbbp2DaNWRo8Oo1gkt82YAY8/DldeCU2bxh1NtSgxiGSj668Po1duvDHuSCQZ7jB4MOyyC/z2t3FHU21KDCLZqEMHOO88+Nvf4N13445GauvZZ+GFF+Caa2CHHeKOptqUGESy1bXXQt26oX1acs+mTaH5qG1bGDAg7mhqRIlBJFvtvjtcdhncfz+8/nrc0UhN3XtvmPtq5EioVy/uaGokqUn04qJJ9KRgfPUVtG8P++yT9oXiJYW+/TY0B7ZsGTqfs+S/W6am3RaRdGrUKNwJPX16RuftlySNHg0rVoQhx1mSFGpCFYNIttu4MUy45g7z5kFxcdwRybZ89FGo8nr2hMceizuaclQxiOSLunXhD3+ARYtg3Li4o5GqXHddGGo8alTckdRaUonBzHY2s8lmtjj6s0klx20ys9nR44mE7e3M7FUzW2JmD0arvYnI1vr0gR//OMy38+WXcUcjlXnrLRg7Fn71K9hrr7ijqbVkK4bBwHPu3gF4LnpdkTXu3iV6HJ+wfRQw2t3bA58B5yYZj0h+Mgvt1atWhepBstOgQbDddiGB57BkE0Nf4J7o+T2EdZurJVrn+Uhgy0K3NTpfpOB06wannx5WAPvgg7ijka1NnQpPPAG/+x2UlMQdTVKSTQy7uvvK6PmHwK6VHNfAzErNbIaZbfnybwp87u4bo9dlQPZPOygSpxEjQvVw5ZVxRyKJNm4MU160aZNTU19Upso1n81sCrBbBbuuSnzh7m5mlQ1xauPuy81sD+B5M5sH1GhOYTMbAAwAaN26dU1OFckfrVuH5ophw+CCC0K/g8Tv738PI8YeeQQaNow7mqQlNVzVzBYBPdx9pZk1B6a5+95VnHM38CTwCLAK2M3dN5rZocAwd+9V1XU1XFUK2po1sO++sNNOMHNmGLUk8Vm9OtypJq8xAAANT0lEQVTM1rUrTJmS1fctZGq46hNA/+h5f+DxCgJpYmb1o+fNgMOABR4y0lTgpG2dLyJbadgw9DPMnQt33BF3NHLNNWGk2E03ZXVSqIlkE8NI4GgzWwz0jF5jZt3MbMuA632BUjObQ0gEI919QbRvEDDQzJYQ+hzGJxmPSGH46U/hiCPCBHurV8cdTeGaMyc0I/3mN9CpU9zRpIzufBbJVfPmheaLX/1Ka0THwR169IAFC+Dtt6FJhbdxZRXd+SyS7zp3hl//OqzZMGdO3NEUnkmTwhxWI0bkRFKoCVUMIrnss89Cx2fHjpp9NZO++SbMeFtSEqZELyqKO6JqUcUgUgiaNIHf/z78cr3nnqqPl9QYOhTKyuCvf82ZpFATSgwiue7cc+Gww+Dyy+GTT+KOJv/NmgVjxoS+nR/+MO5o0kKJQSTX1akTRsZ8+WVIDpI+mzaFZTqbNQuVWp5SYhDJBx07hmky7rkHnn8+7mjy1223QWlpqBjyrMM5kTqfRfLFmjVhpFKdOuHmtwYN4o4ov5SVhTvOf/QjeOqpnOzoV+ezSKFp2DAMXV28OAyhlNS6+OLQlHTbbTmZFGpCiUEkn/TsCWecASNHwsKFcUeTPx5/PCzTOXQotGsXdzRpp8Qgkm/+/Gdo1AjOOy/8wpXkfP55mPKic2cYODDuaDJCiUEk3+yyS5jQ7eWX4S9/iTua3HfxxfDhhzB+PBQXxx1NRigxiOSj008PE+1dfTW8+Wbc0eSuxx6De++Fq66C7t3jjiZjlBhE8pFZ6IjeaSc46yzYsCHuiHLPxx+Hm9i6dg2JoYAoMYjkq5ISGDs23Kl7ww1xR5Nb3MMKeV98ARMmQL16cUeUUUoMIvnshBPgzDPhxhvDZG9SPffdB48+Ctdfn1frLFSXbnATyXeffx5G1OywA7zxRl6sSZxWZWUhGXTqBC+8kFeT5GXkBjcz29nMJpvZ4ujP790jbmZHmNnshMdaMzsh2ne3mb2bsK9LMvGISAUaNw4jat56K0ybIZXbtAnOPjv0ydx9d14lhZpItilpMPCcu3cAnotel+PuU929i7t3AY4EvgWeTTjkii373X12kvGISEWOOQYuvTRME/3ww3FHk72uvx6eew5uvhnat487mtgkmxj6Alsmgb8HOKGK408Cnnb3b5O8rojU1MiRcMgh8ItfwJIlcUeTfaZMgeHDwyiuX/wi7mhilWxi2NXdV0bPPwR2reL4U4AHttp2o5nNNbPRZla/shPNbICZlZpZ6apVq5IIWaRA1asXlqMsLoaTT4a1a+OOKHusWAGnnRYmySuAuZCqUmViMLMpZvZmBY++icd56MWutCfbzJoDnYFnEjYPAfYBugM7A4MqO9/dx7p7N3fvVlJSUlXYIlKR1q3D8MvZs+GSS+KOJjts3AinnBKW63z4Ydh++7gjil3dqg5w956V7TOzj8ysubuvjL74P97GW/0ceMzd/3unTUK1sc7M7gK0yohIuv3kJzBoEIwaBYcfHn4pF7Jrr4UXXwx3OO+7b9zRZIVkm5KeAPpHz/sDj2/j2FPZqhkpSiaYmRH6J3Tvvkgm3HBDWFdgwABYsCDuaOLz5JNhJbZf/jLMSitA8olhJHC0mS0GekavMbNuZjZuy0Fm1hZoBbyw1fn3mdk8YB7QDNDtmSKZULcuTJwY7m04/nhYvTruiDJvwYJQLXXtGiYdlP/SDW4iheyVV6BHj7Co/bPPFszsoaxeDQcfHPoVXn8dWrWKO6KM0ApuIlK1Qw8NN79NmxbWHMjBH4o1tm4d/OxnsHx5WICnQJJCTVTZ+Swiee6MM8JqbyNGhFFLV18dd0Tps3lzuLP5hRfCfEiHHBJ3RFlJiUFEQmd0WRlccw20aAHnnBN3ROlx5ZWhb2XUKI3G2gYlBhEJN3TdcQesXBlG6DRpEmZmzSejRoVlTy+8EK64Iu5ospr6GEQkqFcPHnkEunWDfv3g3/+OO6LUueUWGDwYTj0Vxowp+Dubq6LEICLfadQoJISOHeHEE8OEcrlu7NiwbvOJJ8I99xTsjKk1ocQgIuU1bhyGrrZvH+6S/r//izui2rvpprA85//+LzzwQOEMx02SEoOIfF+zZmEIa6dO4Zd2rk3V7R5GWV1yCfz0p2E1tvqVztEpW1FiEJGKNW0ampIOPhh+/vPcuTt448bQwXzVVXD66fDgg0oKNaTEICKV22mn0Kx0wgnh1/fFF4dVzrLVV19B375h6uwrrwwzydbV4MuaUmIQkW3bbjt46CEYODCM7unVC7JxTZS33go3rD3zDPztb2F4ah19xdWG/tZEpGpFReEegDvvhJdeggMPhJdfjjuq70yaBN27wyefhMTwq1/FHVFOU2IQkeo755yQEIqL4X/+B373uzD3UFw++yz0I/TrFzrK33gDjjoqvnjyhBKDiNRM165hBbhzzglrGXTvDv/5T2ZjcA/NW506hWrhuutg+nRo2TKzceQpJQYRqbkdd4Rx4+Bf/4LPPw/Vw+mnw/vvp//as2bBkUeGkVIlJTBjRliFTfcopIwSg4jUXp8+YWbWa64J02m0bx/mWlq6NPXXKi0NI44OPBDmzg0jj2bOhIMOSv21ClxSicHMTjaz+Wa22cwqXfzBzHqb2SIzW2JmgxO2tzOzV6PtD5pZvWTiEZEYbL89DB8OixeHTt8JE0KC6NUrNPd8+23t3/vTT8OUFt27h8f06aHZaMkSuOACTW+RJkmt4GZm+wKbgb8Dl7v795ZVM7Mi4G3gaKAMeB041d0XmNkk4FF3n2hmfwPmuPvtVV1XK7iJZLEVK8JMrePHw7Jl0KBBaPo54ojw6/6AA8LsrVtPZOcOH34YmopmzoTJk8MIqM2bQ1/CL38J/fuHeyukVqq7gltKlvY0s2lUnhgOBYa5e6/o9ZBo10hgFbCbu2/c+rhtUWIQyQGbNoVpNZ58MjyWLPluX8OG0Lx5SBoAX38dksL69d8d06ULHHdcWJP6oIM0I2oKVDcxZOKWwBbAsoTXZcAhQFPgc3ffmLC9RWVvYmYDgAEArVu3Tk+kIpI6RUVh6OhRR8Ho0eGmuJkzYf78UFUkJoKGDWH33cMiQfvvH0Y+7bhjvPEXsCoTg5lNAXarYNdV7v546kOqmLuPBcZCqBgydV0RSZGSEujdOzwkq1WZGNy9Z5LXWA4krrbdMtq2GmhsZnWjqmHLdhERiVEmhqu+DnSIRiDVA04BnvDQuTEVOCk6rj+QsQpEREQqluxw1RPNrAw4FPg/M3sm2r67mT0FEFUDFwLPAAuBSe4+P3qLQcBAM1tC6HMYn0w8IiKSvJSMSso0jUoSEam56o5K0p3PIiJSjhKDiIiUo8QgIiLlKDGIiEg5Odn5bGargNrM79sM+CTF4WRarn+GXI8f9BmyhT5DzbVx95KqDsrJxFBbZlZanR75bJbrnyHX4wd9hmyhz5A+akoSEZFylBhERKScQksMY+MOIAVy/TPkevygz5At9BnSpKD6GEREpGqFVjGIiEgVlBhERKScgkgMZtbbzBaZ2RIzGxx3PLVhZnea2cdm9mbcsdSGmbUys6lmtsDM5pvZb+OOqabMrIGZvWZmc6LPcF3cMdWGmRWZ2SwzezLuWGrLzN4zs3lmNtvMcm5GTTNrbGYPm9lbZrYwWto4a+R9H4OZFQFvA0cTlg99HTjV3RfEGlgNmdnhwNfABHfvFHc8NWVmzYHm7v6GmTUCZgIn5NJ/BzMzYHt3/9rMioH/AL919xkxh1YjZjYQ6Abs6O594o6nNszsPaCbu+fkDW5mdg/woruPi9ap2c7dP487ri0KoWI4GFji7kvdfT0wEegbc0w15u7TgU/jjqO23H2lu78RPf+KsDZHpWt8ZyMPvo5eFkePnPplZWYtgZ8A4+KOpVCZ2U7A4UTrz7j7+mxKClAYiaEFsCzhdRk59oWUb8ysLdAVeDXeSGouaoaZDXwMTHb3XPsMY4Argc1xB5IkB541s5lmNiDuYGqoHbAKuCtq0htnZtvHHVSiQkgMkkXMbAfgEeASd/8y7nhqyt03uXsXwhrlB5tZzjTrmVkf4GN3nxl3LCnwI3c/EDgW+E3U1Jor6gIHAre7e1fgGyCr+j4LITEsB1olvG4ZbZMMi9rlHwHuc/dH444nGVHpPxXoHXcsNXAYcHzUPj8RONLM/hFvSLXj7sujPz8GHiM0GeeKMqAsodp8mJAoskYhJIbXgQ5m1i7q5DkFeCLmmApO1HE7Hljo7n+JO57aMLMSM2scPW9IGNDwVrxRVZ+7D3H3lu7elvDv4Hl3PyPmsGrMzLaPBjAQNcEcA+TMaD13/xBYZmZ7R5uOArJqEEbduANIN3ffaGYXAs8ARcCd7j4/5rBqzMweAHoAzcysDBjq7uPjjapGDgPOBOZFbfQAv3P3p2KMqaaaA/dEI93qAJPcPWeHfOawXYHHwm8N6gL3u/u/4w2pxi4C7ot+rC4Fzok5nnLyfriqiIjUTCE0JYmISA0oMYiISDlKDCIiUo4Sg4iIlKPEICIi5SgxiIhIOUoMIiJSzv8Dr7d+3DDbzusAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def transport(f0, f1, f2):\n", " # anti bounce back boundary conditions\n", " f1[-1] = -f2[-2]\n", " f2[0] = -f1[1]\n", " # transport\n", " f1[1:-1] = f1[2:]\n", " f2[1:-1] = f2[:-2]\n", "\n", "\n", "# parameters\n", "c = 1 # velocity for the transport equation\n", "Tf = 2*np.pi # final time\n", "N = 128 # number of points in space\n", "la = 1. # scheme velocity\n", "s = 2. # relaxation parameter\n", "# initialization\n", "x = mesh(N) # mesh\n", "dx = x[1]-x[0] # space step\n", "dt = dx/la # time step\n", "f0, f1, f2, m0, m1, m2 = initialize(x, c, la)\n", "plt.figure(1)\n", "plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\\rho$')\n", "plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')\n", "plt.title('Initial moments')\n", "plt.legend(loc='best')\n", "# time loops\n", "nt = int(Tf/dt)\n", "m2f(f0, f1, f2, m0, m1, m2, la)\n", "for k in range(nt):\n", " transport(f0, f1, f2)\n", " f2m(f0, f1, f2, m0, m1, m2, la)\n", " relaxation(m0, m1, m2, c, s)\n", " m2f(f0, f1, f2, m0, m1, m2, la)\n", "plt.figure(2)\n", "plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\\rho$')\n", "plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')\n", "plt.title('Final moments')\n", "plt.legend(loc='best')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }