{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spreading out"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You've reached the fourth module of the open course [**\"Practical Numerical Methods with Python\"**](https://openedx.seas.gwu.edu/courses/course-v1:MAE+MAE6286+2017/about), titled *Spreading out: Parabolic PDEs*. We hope that you are enjoying the ride of [#numericalmooc](https://twitter.com/hashtag/numericalmooc) so far!\n",
"\n",
"We introduced finite-difference methods for partial differential equations (PDEs) in the [second module](https://github.com/numerical-mooc/numerical-mooc/tree/master/lessons/02_spacetime), and looked at convection problems in more depth in [module 3](https://github.com/numerical-mooc/numerical-mooc/tree/master/lessons/03_wave). Now we'll look at solving problems dominated by diffusion.\n",
"\n",
"Why do we separate the discussion of how to solve convection-dominated and diffusion-dominated problems, you might ask? It's all about the harmony between mathematical model and numerical method. Convection and diffusion are inherently different physical processes.\n",
"\n",
"The titles of the course modules are meant to spark your imagination: \n",
"\n",
"* _Riding the wave_—imagine a surfer on a tall wave, moving fast towards the beach ... convection implies transport, speed, direction. The physics has a directional bias, and we discovered that numerical methods should be compatible with that. That's why we use _upwind_ methods for convection, and we pay attention to problems where waves move in opposite directions, needing special schemes.\n",
"\n",
"* _Spreading out_—now imagine a drop of food dye in a cup of water, slowly spreading in all directions until all the liquid takes a uniform color. [Diffusion](http://en.wikipedia.org/wiki/Diffusion) spreads the concentration of something around (atoms, people, ideas, dirt, anything!). Since it is not a directional process, we need numerical methods that are isotropic (like central differences)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import Image\n",
"Image(url='http://upload.wikimedia.org/wikipedia/commons/f/f9/Blausen_0315_Diffusion.png')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Parabolic PDEs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You already met the simplest parabolic PDE—the [1-D diffusion equation](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_03_1DDiffusion.ipynb)—in module 2. Its main feature is that it has a second-order derivative in space. Here it is again:\n",
"\n",
"$$\n",
"\\frac{\\partial u}{\\partial t} = \\alpha \\frac{\\partial^2 u}{\\partial x^2}\n",
"$$\n",
"\n",
"Check out the article on [parabolic PDEs](http://en.wikipedia.org/wiki/Parabolic_partial_differential_equation) in Wikipedia. Now compare with the diffusion equation above, with the two independent variables here being $x, t$. You'll see that with no mixed derivatives, and only one second-order derivative (in the spatial variable $x$), it satisfies the condition of a parabolic PDE. Work it out on a piece of paper if you need to.\n",
"\n",
"In the previous module, discussing hyperbolic conservation laws, we learned that solutions have characteristics: information travels along certain paths on space-time phase space. In contrast, parabolic equations don't have characteristics, because any local change in the initial condition will eventually affect the entire domain, although its effect will be felt at smaller intensity with longer distances. This is typical of diffusion processes.\n",
"\n",
"In this first lesson of the module, we first review the 1D diffusion equation and then take a deeper look into the issue of boundary conditions. In the next notebook, we'll introduce _implicit discretizations_ for the first time, which will take us to the land of linear solvers. In the third lesson we'll graduate to two dimensions—more boundary condition and stability issues will come up. We'll then study 2D implicit methods, and go into Crank-Nicolson method: perhaps the most popular of them all. _Enjoy!_"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Heat conduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Heat conduction is a diffusive process. Let's remind ourselves of the heat equation in one spatial dimension:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial T}{\\partial t} = \\alpha \\frac{\\partial^2 T}{\\partial x^2}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Here, $\\alpha$ is the thermal diffusivity, a property of the material, and $T$ is the temperature.\n",
"\n",
"In the [third lesson of module 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_03_1DDiffusion.ipynb), we discretized the diffusion equation with a forward-time, centered-space scheme, subject to the following stability constraint:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\alpha \\frac{\\Delta t}{(\\Delta x)^2} \\leq \\frac{1}{2}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Let's look into it more deeply now, using a 1D temperature-evolution problem."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Problem set up"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Say we have a graphite rod, with [thermal diffusivity](http://en.wikipedia.org/wiki/Thermal_diffusivity) of $\\alpha=1.22\\times10^{-3} {\\rm m}^2/{\\rm s}$, length $L=1{\\rm m}$, and temperature $T=0{\\rm C}$ everywhere. At time $t=0$, we raise the temperature on the left-side end, $x=0$, to $T=100{\\rm C}$, and hold it there. *How will the temperature evolve in the rod?*\n",
"\n",
"As usual, start by importing some libraries and setting up the discretization. We'll begin by using a spatial grid with 51 points and advance for 100 time steps, using a forward-time/centered scheme."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![graphite-rod](./figures/graphite-rod.png)\n",
".\n",
"#### Figure 1. Graphite rod, temperature fixed on ends."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy\n",
"from matplotlib import pyplot\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Set the font family and size to use for Matplotlib figures.\n",
"pyplot.rcParams['font.family'] = 'serif'\n",
"pyplot.rcParams['font.size'] = 16"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Set parameters.\n",
"L = 1.0 # length of the rod\n",
"nx = 51 # number of locations on the rod\n",
"dx = L / (nx - 1) # distance between two consecutive locations\n",
"alpha = 1.22e-3 # thermal diffusivity of the rod\n",
"\n",
"# Define the locations along the rod.\n",
"x = numpy.linspace(0.0, L, num=nx)\n",
"\n",
"# Set the initial temperature along the rod.\n",
"T0 = numpy.zeros(nx)\n",
"T0[0] = 100.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remember the forward-time, centered-space discretization? You should work it out on a piece of paper yourself (if you can't do it without looking it up, it means you need to do this more!).\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{T_{i}^{n+1}-T_{i}^{n}}{\\Delta t}=\\alpha\\frac{T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n}}{\\Delta x^2}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"To obtain the temperature at the next time step, $T^{n+1}_i$, from the known information at the current time step, we compute\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"T_{i}^{n+1}=T_{i}^{n}+\\frac{\\alpha\\Delta t}{\\Delta x^2}(T_{i+1}^{n}-2T_{i}^{n}+T_{i-1}^{n})\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Check the [third notebook of module 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_03_1DDiffusion.ipynb), if you need to refresh your memory!\n",
"\n",
"The following function implements this numerical scheme:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def ftcs(T0, nt, dt, dx, alpha):\n",
" \"\"\"\n",
" Computes and returns the temperature along the rod\n",
" after a provided number of time steps,\n",
" given the initial temperature and thermal diffusivity.\n",
" The diffusion equation is integrated using forward \n",
" differencing in time and central differencing in space.\n",
" \n",
" Parameters\n",
" ----------\n",
" T0 : numpy.ndarray\n",
" The initial temperature along the rod as a 1D array of floats.\n",
" nt : integer\n",
" The number of time steps to compute.\n",
" dt : float\n",
" The time-step size to integrate.\n",
" dx : float\n",
" The distance between two consecutive locations.\n",
" alpha : float\n",
" The thermal diffusivity of the rod.\n",
" \n",
" Returns\n",
" -------\n",
" T : numpy.ndarray\n",
" The temperature along the rod as a 1D array of floats.\n",
" \"\"\"\n",
" T = T0.copy()\n",
" sigma = alpha * dt / dx**2\n",
" for n in range(nt):\n",
" T[1:-1] = (T[1:-1] +\n",
" sigma * (T[2:] - 2.0 * T[1:-1] + T[:-2]))\n",
" return T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are all set to run! First, let's use a time step `dt` that satisfies the stability constraint."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Set the time-step size based on CFL limit.\n",
"nt = 100 # number of time steps to compute\n",
"sigma = 0.5\n",
"dt = sigma * dx**2 / alpha # time-step size\n",
"\n",
"# Compute the temperature along the rod.\n",
"T = ftcs(T0, nt, dt, dx, alpha)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Note\n",
"-----\n",
"In the function `ftcs`, we *copy* the array `T0`; it ensures that the value of `T0` remains unchanged for us to reuse.\n",
"\n",
"Confused? Take a look at [Lesson 0](./04_00_Python_Function_Quirks.ipynb) for a more in-depth explanation.\n",
" \n",
"-----\n",
"\n",
"Now plot the solution!"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEbCAYAAACV0PCVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XecVNX9//HXZzu9F+lFBZEqFiwYLFijiTH2FhOD+SZffzFf49dYYzQmRmMSY/JVUaNRY+wlsRsVewEUERBQBKTD0pcFtn1+f9y77jhum52ZvbOz7+fD+7g755575zPHZT577z33HHN3REREopYTdQAiIiKghCQiIhlCCUlERDKCEpKIiGQEJSQREckISkgiIpIRIk1IZraLmT1vZup7LiLSykWWkMzsBOAdYGgD9fLN7Fozm29mc8zsbTM7qI66F5rZPDObbWYfmNm30xG7iIikXpRnSL8AJgNvNVDvFuAUYKK7jwT+BrxkZmNjK5nZL4ArgOPcfTRwCfCImR2d8shFRCTlLKqRGswsz90rzOwe4Bx3t1rqDAM+Ac5z97/FlM8Flrj7seHrzsAK4CZ3vyqm3jPAIHffM72fRkREkhXZGZK7VzSi2gmAAa/Glb8CHGFm7cPXRwFt66g3wsyGJxOriIikX6b3shsNVAFfxJUvBvKAETH1qsvj68VuFxGRDJUXdQAN6A6UuntlXPmWcN0tph7A1gbqfYWZTQGmAOQWtR8/dECf5KLNElVVVeTkZPrfKs1DbVFDbVFDbVFj4cKFxe7eIxXHyvSEVJev3W9qSj13nwpMBSjcZTd/8KV3GTegS7KxtXjTpk1j0qRJUYeREdQWNdQWNdQWNcxsaaqOlekpvhhoa2a5ceUdwvX6mHqx5XXVq9fDM5YlHKCIiKRGpiek2QQx9o8rHwxUEPTAq64HMKiWerHb6/Xvj1ZRWtaYvhYiIpJqmZ6QngAcmBRXfgjwortX3zN6Hiito948d5/f0BsV5kLJzgqemb0qqYBFRKRpMjohufsCgns8l5pZdwAz+z7B6A6Xx9TbBFwL/MTMhoT1DgeOBH7emPdqXxDcbtJlOxGRaETWqcHMbiQYqWFA+HpWuGlfdy+LqXoB8EvgLTMrJ+hJd4S7z4o9nrtfb2Y7gKfNrAKoBE5y9+caE0+7PCO/IJfpSzayaF0JQ3u0b3gnERFJmSgfjL3Y3ce6e1d3t/DnsXHJCHcvd/cr3H2Yu4909/3d/Y06jvkndx/h7qPdfZy7P9nYeHIMvjl6F0BnSSIiUcjoS3bN7ZR9gr4Tj81cQXllVcTRiIi0LkpIMfYa0IWhPdpRXLKTV+evjTocEZFWRQkphpl9eZaky3YiIs1LCSnOd/bqR16O8eqCdazdsiPqcEREWg0lpDjd2xdy2B49qaxyHv1gedThiIi0GkpItai+bPfIjOVENV+UiEhro4RUi4N360GvjoUsLt7G+4s3RB2OiEiroIRUi7zcHL47vh8AD01X5wYRkeaghFSHk/cOLts9O2cVW3aURxyNiEj2U0Kqw8Bu7ZgwpCs7yqv416yVUYcjIpL1lJDqceo+AwB4cHr8DOoiIpJqSkj1OGpkbzq1yWfOii3MWbE56nBERLKaElI9ivJzOWFcXwD++b7OkkRE0kkJqQGn7Rtctntq1krNJisikkZKSA0Y1rsD4wZ0pmRnBU9rNlkRkbRRQmqE06o7N+iynYhI2ighNcI3x+xC+8I8PvhiEwvXbI06HBGRrKSE1AhtC/I4fmwfQJ0bRETSRQmpkaov2z3x4Qp2lFdGHI2ISPZRQmqkUf06sWefjmwqLeeFuaujDkdEJOsoISXg1LALuC7biYiknhJSAr41tg9t8nN59/MNfL6uJOpwRESyihJSAjoW5XPs6F0ATUshIpJqSkgJOm3fYFqKR2cup6yiKuJoRESyhxJSgvYa0IXderZn/bYy/vPJmqjDERHJGkpICTIzdW4QEUkDJaQm+M64vhTk5fDmZ8Us21AadTgiIllBCakJurQr4OiRvXGHh2eoc4OISCooITVR9WyyD89YRkWlOjeIiCRLCamJJgzpyuDu7VizZSfTFqyLOhwRkRZPCamJzIxT9gm6gD84XZ0bRESSlfEJycz2NrPnzOwTM/vYzN43s5Pi6uSb2bVmNt/M5pjZ22Z2ULpjO3GvfuTlGK/MX8vqzTvS/XYiIlktoxOSmQ0CXgaKgVHuPgr4G/CwmR0XU/UW4BRgoruPDOu8ZGZj0xlfjw6FTB7RiyqHR9S5QUQkKRmdkIBjgI7AH9y9AsDdbwO2AKcDmNkwYApwvbuvC+vcCXwOXJfuAKufSXpoxjKqqjzdbycikrUyPSFVhOu86gIzM4K4c8OiEwADXo3b9xXgCDNrn84AJ+7anb6d27B843be/Kw4nW8lIpLVMj0hPQjMB64ws/ZmlgNcBhQCt4V1RgNVQHzPgsUEiWxEOgPMyVHnBhGRVMira4OZlQNNuQa1yt0HNj2kGu6+xcwOA+4muI9UAmwGJrv7a2G17kCpu8dP47olXHer6/hmNoXgch89evRg2rRpTYqzX3kVBrwwZzX/euFVOhZak46TKUpKSprcFtlGbVFDbVFDbZEedSYkggRwZ4LHM+DkpocTd7Dg/tDLwLNAV2BHePzHzexMd3+ugVjq5e5TgakAw4YN80mTJjU51mdWT+fl+WtZVTSA478xtMnHyQTTpk0jmbbIJmqLGmqLGmqL9KgvIa129ysTPaCZHZtEPPGuBToDP3X37WHZg2Z2CvB3M+tDkDjbmllu3FlSh3C9PoXx1OnUfQfw8vy1PDR9GVMOHkJwq0tERBqrvntIzzfxmE3drzajgOUxyajaQqAHMBiYTfA5+sfVGUzQKeKTFMZTp0OG9aBXx0I+L97Ge4s3NMdbiohklToTkrtf2pQDNnW/OqwFdjGz+DO5gQT3tzYCT4Q/T4qrcwjwortvTWE8dcrLzeHkvcPODZqWQkQkYUn1sgt7vaXTLQTPIV0TdvfGzA4BvgM85O7F7r6A4D7QpWbWPazzfWAocHma4/uKk/fujxk8O2c1m0rLmvOtRURavHoTipn1DYfqed/Mrq2lyuFm9h8z2z0dwbn7o8BRwP7APDObA9xMkGi+F1P1AuAR4K2wzg+BI9x9Vjriqkv/rm05aNfulFVU8cSHK5rzrUVEWryGznBOAMYBLxEMxxNvLrANmGZmfVMcGwDu/oK7H+Lue7j7SHcf7e43uvvOmDrl7n6Fuw8L6+zv7m+kI56GnBaO3PDAe1/grpEbREQaq6GEdCxwgbtf7u6L4ze6+wp3/xZwL/C/6QiwpTl8j1706FDIp2tLePdzdW4QEWmshhJSX+CORhznGuAbyYfT8hXk5Xx5lnTfu0uiDUZEpAVpKCFV1TICwte4eynB8D0CnLHfAPJyjBfmrtG0FCIijZTKXnL5KTxWi9arYxFH7tmbyirngfeWRh2OiEiL0FBCWmxmDV6KM7ODAX3zxjhr/2A4vwfeX0ZZhU4eRUQa0lBCuhl4wMwm1lUhnJn1H8CNqQyspdtvcFd279We4pKdPD93ddThiIhkvHoTkrtPI3jo9DUz+9DMbjWzq83sl+HPHwKvATfHjL4tgJlx9v6DALjvnSVRhiIi0iI0eA/J3X8FnEgwId75wFXAL8OfDTje3X+fziBbqhPG9aVDYR7Tl2xk3sotDe8gItKKNapTg7s/4e6jCbqBHxgu/dx9rLs/k84AW7J2hXmcOL4foC7gIiINSaiXnbuvcvd3wmVluoLKJmdOCDo3PPnhSjZvL484GhGRzFVnQjKzq5tywKbul6127dmeA3ftxvbySh6duTzqcEREMlZ9Z0jHN/GYTd0va1V3brj/3aVUVWl8OxGR2tQ3Y2xvM7uURkwFHqddEvFkpcOG96RPpyIWF2/jzc+KOXj3HlGHJCKScepNSMB1TTimHrqJk5ebwxkTBnLjCwu4950lSkgiIrWo75JdfhOXfmmMt8U6ZZ/+FOTm8PL8tSzbUBp1OCIiGae+Kcwrm7honJxadG9fyDGjeuMO/3hPU5yLiMRL9xTkEuPsAwYB8ND0L9hR3uAg6iIirYoSUjMa178zI/t2ZGNpOU/PXhV1OCIiGUUJqRmZGWdPGARofDsRkXhKSM3suDF96NQmn4+Wb+ajZZuiDkdEJGMoITWzNgW5nLJPfwDufUdTSImIVEs4IZnZAWZ2mZldF74+yMz0MGwCztxvIGbw79kr2bCtLOpwREQyQqMTkpm1M7PngDeBXwPnhpuOBz42s0Epjy5LDejWlkm796CsooqHZyyLOhwRkYyQyBnS9UAX4DhgCLAWwN3/l2COpN+mPLosVjN531IqNb6diEhCCekYYLK7P+PuS4AvH4B19/uB3VMcW1b7xu49GNC1LSs2befV+WujDkdEJHKJJKRyd99az/aOyQbTmuTkGGeFcyXd+646N4iIJJKQSs3sW7VtMLMjAfVhTtBJe/ejMC+H1xeuY3HxtqjDERGJVCIJ6TrgcTN72cyuBbqb2S/M7J/AU8A1aYkwi3VuW8C3xvYBgrmSRERas0YnJHd/DDgLGAZcTjCq92+Ag4Bz3P3faYkwy1V3bnh4xjJKyyqiDUZEJEIJPYfk7g8A/YFRwCHheoC7P5SG2FqFkX07sdeAzmzdUcFTs1ZGHY6ISGQSeQ5pnZmtAAa7+1x3fy1cp73PspmdaGavm9lMM/vczGaY2Vkx2/PN7Fozm29mc8zsbTM7KN1xpUr1WdK97yylGZpTRCQjJXKGlAPs7e6fpyuY2pjZzwguEZ7u7uMJLhkuBA6LqXYLcAow0d1HAn8DXjKzsc0Za1MdPao33dsX8MmqLby3eEPU4YiIRCKRhDTX3eucM8HMvpmCeOKPOYjggdzz3X05gLuXAz8H/hLWGQZMAa5393VhnTuBz2naFOzNrjAvlzPDLuC3v7Yo4mhERKKRSEK608x+bmZ17ZOOXnZnAZvcfXpsobuvdPcZ4csTAANejdv3FeAIM2ufhrhS7pz9B9EmP5dXF6xj/uotUYcjItLsEklIZwEXAavDezQvxi7A0DTEdwCwJLyH9EZ4j+htM/t+TJ3RBKNGxM8LvhjIA0akIa6U69Ku4MtRwKe+1qxXRUVEMkJeAnUnALNiXreJ227Jh/M1/YFBBJfoTiAYP+9E4J9mtou7Xwd0B0rdPX5O8OrTjG51HdzMphBc7qNHjx5MmzYtpcEnas+8KnIMnpq1ggM7bqBbm2hmBykpKYm8LTKF2qKG2qKG2iI9EklIn7n7xLo2mtmHKYgnXhHQDrjY3VeHZY+Y2anAZWb2x3r2bTBBuvtUYCrAsGHDfNKkSUmGm7w3t3zIU7NWMreiN1dNiubkbtq0aWRCW2QCtUUNtUUNtUV6JPIn+HENbK8zWSWheuy8WXHlHwJtCS7HFQNtzSw3rk6HcL0+DXGlzZSDhwDw4PQv2FSquZJEpPVIZKSG5Q1UuSDJWGozP1zHx1kZUz47XPePqzMYqAA+SUNcabNnn04cvHsPSssquU8zyopIK5LIg7EH1LcA30tDfNXDEY2OKx8JbAfmAk8ADkyKq3MI8GIDI5RnpB+FZ0n3vL2EHeXxt8ZERLJTIveQ3iT44m9ODwEXAr82s2+6e4mZTQS+C1zj7tuABWY2FbjUzJ529+KwF95Q4Mxmjjcl9h/ajVF9O/Hxis08MnP5l9NUiIhks0QS0iLgR3Fl7YDhBNOY35yqoKq5e6WZHQX8DphrZjuAncB/u/sdMVUvAH4JvGVm5QT3no5w9/h7Ty2CmfGjbwzlJw98wB2vf87p+w4gNycdnRhFRDJHIgnpOnd/uZbyf5nZ3cCfgEdSE1YNd98A/LCBOuXAFeGSFY4a2ZsBXdvyxYZSnpuzim+O7hN1SCIiaZVIp4Z76tm2jq/f55Ek5OYYPwzvJd3+2ucadFVEsl4inRr61LL0NbORZnYlwTNDkkInje9Ht3YFfLxiM28valG910VEEpbIc0jLgWVxyxfAR8AvCEbklhQqys/lewcMAuA2DboqIlkukXtIy/j6AKqVBMP5vOfu+hM+Dc7afyC3vraINz4tZu7KzezZp1PUIYmIpEUiCekxd78rbZFIrTq3LeDUfQbwt7cWc/trn/Pn08ZFHZKISFokcsnuqdoKzWw/M7vNzAanKCaJ84OJg8nLMZ75eBXLNpRGHY6ISFokkpDqGsh0NbAZuD/5cKQ2fTu34fgxfaiscu56c3HU4YiIpEUiCanWJzPdfam7XwK0iInwWqop36gZdHXDNg26KiLZp957SGZ2HDWjfPcLh+j5WjWgH5Cf4tgkxvDeHZk0rAfTFqzj3neWcOHhu0cdkohISjXUqWFX4Ojw5w4xP8cqJ5idtd7RFCR5P/rGUKYtWMff317C+QcPpU1B/IwbIiItV72X7Nz9j+7e3937A59U/xy3DHH3w9z9rWaKudXab3BXxvTvzMbSch6esSzqcEREUiqRe0jfTlsU0ihmxn+F95LueONzKiqrIo5IRCR1EhnLrt7Z4szsueTDkYZMHtGbwd3bsXzjdp6ds7rhHUREWohEHozFzHYBTgWGAIVxm/dOVVBSt9wc44cTh3DZEx9z27RFHDd6F8w0NYWItHyNTkhmtg/wH6CMoIPDunBTd4LkpD/Xm8l39urLH15ayLxVW3jzs2Im7tYj6pBERJKWyD2k64Hz3L0HMR0cCCbpuxz4SzoClK8rys/l3AMHARp0VUSyRyIJqbu7f20CPnevcvffAoelLixpyJkTBtKuIJe3PlvPx8s3Rx2OiEjSEklIO2NfmFmHmJ8Lgd1SFZQ0rFObfE7fbwAAf37l04ijERFJXiIJaauZnWPBHfQZwONmdqyZHQs8AaxKS4RSpx8ePIQ2+bm8NG8NH3yxMepwRESSkkhCug04ExgKXBuu/wX8G9gfuCjl0Um9enYo4gcHBYOs/+65+ZrmXERatESeQ3rE3Se7+2fu/gUwGvgm8B1gN3d/M11BSt2mfGMIndvm897iDbz+aXHU4YiINFmjE5KZvR8ugwHcvcTdn3P3J91d34QR6ViUz48nDQXghufnU1WlsyQRaZkSuWS3B/ATd9eEPBnm7P0H0btjEXNXbuGZj3UrT0RapkQS0mx3n17XRjPbKwXxSBMU5edy4eFBJ8ebXlxAuca4E5EWKJGE9KSZnVzP9juTDUaa7rvj+zGkezuWrC/loekaCVxEWp5EEtLuwJ/M7EMze8DMpsYuQP80xSiNkJebw8+PHAbAzS9/yvayyogjEhFJTCKDq54NrCUYu25iLds71FImzejokb0Z1bcTH6/YzN1vL+bHk3aNOiQRkUZL5AxpXh0T9H05gV+6gpTGMTMuOWo4ALdNW8Tm0vKIIxIRabxEEtL5DWw/JZlAJDUO2q07B+7ajS07KrhVA6+KSAuSyIOx71f/bGZ9zGx0+HNOuH1h6sOTpvjfI4OzpLvfWszqzTsijkZEpHESOUPCzE41s4XAMqB6htj7zez31kyzxJnZG2bmZjaoOd6vJRrTvzNHj+zNzooqDbwqIi1GIiM1nAb8HfiIYCy7reGmS4DBwC9SHt3XYzgROKiObe3N7C9mtsDM5pnZi2a2Z7pjylQXHTGMHIOHpi9jcfG2qMMREWlQImdIlwCHu/tJ7n41UArg7ssIeuB9N/Xh1TCzAuC3wLN1VHkEGAeMc/cRwHvANDPrm864MtWuPdtz0vj+VFY5N724IOpwREQalEhCKnD3N2rb4O7bSKwLeVP8hGDai6+NFmFmk4GjgCvdvTQsvhbIBS5Lc1wZ66eH70ZBXg5Pz17FnBWaxE9EMlsiCanIzHrUtsHMepLG55DMrCtwMXUnlxOBcuDLEcfdvQx4K9zWKvXp3IZz9h8IwA0v6CxJRDJbQkMHAW+Y2VlmNhTIMbNeZnYk8AzwcFoiDFwF3O/uS+rYPhpYGSahWIuBXmHCbJV+PGlXOhTm8frCdby9SIOyi0jmSuQy22XACIKODQ4YsDLc9hxB0kg5M9sVOJlgtPG6dKemk0WsLeG6G8EoE/HHngJMAejRowfTpk1LKtZMNXmA8fincMXD07lyQhENdYgsKSnJ2rZIlNqihtqihtoiPRqdkNx9B3CUmR0FHE7wJV8MvOTuL6YpPoAbgOvdvSk3Qer95nX3qcBUgGHDhvmkSZOa8BaZb5/9K3j9xlf5fHMZZT334Mg9e9dbf9q0aWRrWyRKbVFDbVFDbZEeCT2HBODuz7v7z939XHe/OJ3JyMwmAiOBWxuoWkzt97Cqy9anMq6Wpl1hHhccGkxPceMLC6jUJH4ikoESfTC2s5ldZWbPmdlH4fpKM+ucpvgmE/SUm25ms8xsFvCjcNuzYdkxwGygT9g1PNZgYI27f+1yXWtz2r4D6NelDZ+tLeHxD5ZHHY6IyNck8mDsXsAigntFowkSxRjgauAzMxub6uDc/Sp3H+ruY6sX4LZw8zFh2bPA40A+cEBMvAXh68dSHVdLVJCXw0VH7A7An/7zKTvKNT2FiGSWRM6Q/go8DQx0977uPtLd+wADCR5W/b90BNgY4WXDF4BrzaxtWHw5UAX8Jqq4Ms3xY/oyvHcHVmzazj/e+yLqcEREviKRhDTI3c9x9xWxhe6+HPgeMCiFcX2NmR1T2yW7mConEVy6m2VmnxCcHU2Kj7c1y80xLg4n8bvllU/ZsC2+l7yISHQSSUh1zovt7lXA0tgyM+vY1KDqeI9nw0t0vd3d3H1EeAmvevtWd/+Ju+/u7nu4+2R3n5vKGLLBocN7csDQbmwqLefXT8+LOhwRkS8lkpD+GXZo+EpXcTPLM7MrCbtPx5iWbHCSembGb04YRWFeDo9/uII3Pl0XdUgiIkBiCekYguF7VpvZu+Fo2u8Cq4H/Bk4Ly140sxeBoWmIV1JgUPd2/L/Dgm7glz8xh+1l6uAgItFLJCFNAGYRTFVeDrQJ158AC8PXsUuzzI8kTTPl4CEM792BLzaU8qf/aG5FEYleIkMHfebuExtb2cw+bEI80kzyc3O4/sTRnPB/b3Hnm4s5bkwfRvbtFHVYItKKJXKGdFyCx060vjSzsf07c87+g6isci59/GMqKquiDklEWrFGJ6Swe3edwvtGja4vmeHnRw6jT6ciPl6xmXveXhJ1OCLSiiU0qZ6ZHQwcAvQmGKkh1rhUBSXNp31hHr8+YSTfv2cGN724sMGBV0VE0qXRCcnMriYYNmg7sJFgCopYaZugT9Lr0OG9OHb0LjwzexVXPDmH7w3W4Ksi0vwSuYf0A+B4d2/n7v3cvX/sQtDbTlqoXx43go5Feby2cB3vrlI3cBFpfokkpC/c/el6th+RbDASnZ4dirj82GAOxAfm72SjhhUSkWaWSEK618wOqWf7X5MNRqJ18t792W9wV7aWwXXP6oRXRJpXIr3sbge+HY7ScL+ZTY1dCDo7SAtmZvz2O6PIy4FHZy7n7c+Kow5JRFqRROZDuhq4ABgFTAKOjlvUqSELDOnRnuOH5gNw6RMfa94kEWk2iVyy+xFwrDo1ZL9jBuczrFcHlq4v5eaXP406HBFpJRJJSEvd/bl6tqtTQ5bIyzF+e+IozGDq65/zyaotUYckIq1AIgnpCTM7tJ7t6tSQRfYa0IWzJwykssr5xWOzqazSs0kikl6JjNQwFLjAzJYBnwGlcdvVqSHL/PzIYbwwdw0fLd/Mve8s4dwDB0cdkohksUTOkM4GqoC+wDdQp4as16Eon2u/PRKAG19YwIpN2yOOSESyWSIJaV58RwZ1ash+k0f04uiRvSktq+TKJ+fgrkt3IpIeiSSk8xvYfkoygUjm+tXxe9KhKI9X5q/l4RnLog5HRLJUIg/Gvl/9s5n1MbPR4c854XZNO5qlenYs4pfH7QnAlU/NZfbyTRFHJCLZKJEzJMzsVDNbCCwDqruA329mvzczTVmexb47vh+n7zeAsooqfnTfTNaX7Iw6JBHJMomM1HAa8HfgI+BaYGu46RJgMPCLlEcnGeWXx41g3IDOrNy8gwv++aFmmBWRlErkDOkS4HB3P8ndrybs9u3uywh64H039eFJJinMy+XWM8bTvX0Bby9az40vLIg6JBHJIokkpAJ3f6O2De6+jQRnn5WWqXenIv56+l7k5hi3v/45T89eGXVIIpIlEklIRWbWo7YNZtYTPYfUauw3pBuXHxPMnfS/j85mweqtDewhItKwRBLSk8AbZnaWmQ0Fcsysl5kdCTwDPJyWCCUjnXvgIL49tg+lZZX86P6ZbN5eHnVIItLCJZKQLgOWEHRsWAiMBlYCzwJrgKtSHZxkrmDupNHssUtHFhdv438emkWVxrsTkSQk8hzSDnc/CjgG+ANBYvoDcLS7f9PdNed1K9OmIJfbzxxPpzb5vDx/Lbe88lnUIYlIC1ZvRwQzqz7r+czdHwBw9+eB59MdmLQMA7q15eZTx3LuPdP508sLGdWvI4cO7xV1WCLSAjV0hvRfgIVLJMxsrJndYWYzzewjM5tnZn+O72BhZu3N7C9mtiCs86KZ7RlV3K3JpGE9uWjy7rjDTx+cxZLibVGHJCItUEMJabW7/8rd/9Es0dTuQaArcLC7jwEmE0wG+JaZtYmp9wgwDhjn7iOA94BpZta3uQNujX48aVeOGNGLrTsqOP++mWzbWRF1SCLSwjSUkBp9l9rMLk0ylvpcEj7rhLuvAG4EdiO4n4WZTQaOAq509+p5mq4Fcgk6Y0ia5eQYN508hiE92rFgzVYueWy2RgYXkYQ0lJAKzKy/mQ1oaAHOSFOMo909/m559dOYXcL1iUA58GZ1hbCTxVvhNmkGHYrymXrWeNoV5PL07FXc9ebiqEMSkRakoYQ0gqCr9+JGLHukI8A6eu/tTnD29nr4ejSwspa6i4Fe4YO70gx27dmBm04eA8Bvn5vP24uKI45IRFoKq++yipmtAm5rzHGAKe7eJ1WB1flGZrnAB8D77v7DsGwhsNPdR8XV/TVwOTDC3b82gaCZTQGmAPTo0WP8ww/r2V6AkpIS2rdvn9QxHl1YxtOfl9MhH64+oA3d2iQ0sHzGSEVbZAu1RQ21RY1DDjlkprvvnYpjNTT+3GrjSfjhAAAVQklEQVR3/1VjDmRmR6Ugnsa4EqgAftaIuvX2DnT3qcBUgGHDhvmkSZOSDi4bTJs2jWTbYuLBzua73+eNT4u5fX4u/zhvAj06FKYmwGaUirbIFmqLGmqL9EhZpwZ3n5BkLA0ys3OBkwkexi2J2VRM7WPpVZetT3ds8lW5OcafTx3Hbj3bs3BNCadOfYc1W3ZEHZaIZLDGdmqodVDV5mRmZwEXAYe6+9q4zbOBPmZWEFc+GFhTS31pBl3aFfDglAkM792BReu2ccrt77By0/aowxKRDNVQQmoPvAbc0Ayx1MnMzqRmPqbVYdk3w3tAAI8D+cABMfsUhK8fa+ZwJUa39oX884cT2LNPR5asL+WUqe+wbENpwzuKSKtTb0Jy90HuPsTdz22ugOKZ2RnAHcA9wOFmdmaYoI4D+oRxvgi8AFxrZm3DXS8HqoDfNHvQ8hVd2hXwwHkTGNOvE8s2bOfUqe+ydL1GcxCRr2oJXZ9uAYoIHoa9L2aZElfvJIJLd7PM7BOCs6NJ4YO0ErFObfO577z92GtAZ1Zs2s4pt7/L5+tKGt5RRFqNjE9I7t7V3a2O5eqYelvd/Sfuvru77+Huk919boShS5yORfnc+4P92HdQV1Zv2cEpU9/ls7Wa3E9EAhmfkCS7tC/M457v78P+Q7qxbutOTrn9Xeav3hJ1WCKSAZSQpNm1Lcjj7nP3YeJu3Vm/rYzTpr7L3JWbow5LRCKmhCSRKMrP5Y6z9+bQ4T3ZWFrO6Xe8x+zlm6IOS0QipIQkkSnKz+W2M8dzxIhebN5ezhl3vMfMpRujDktEIqKEJJEqyMvhr2fsxbGjdmHrzgrOvus9pi/ZEHVYIhIBJSSJXH5uDjefOpZvje3DtrJKzr7rfd5ZpNGeRFobJSTJCHm5Ofzh5LF8Z6++bC+v5Jy73+fed5Zokj+RVkQJSTJGbo7x+++O4ez9B1JWUcVVT83lh/fOZMO22qbEEpFso4QkGSUnx7jmWyP56+l70aEoj/98soajb35dE/2JtAJKSJKRjh29C8/9dCJ7D+zCmi07OePO97jxhfmUV1ZFHZqIpIkSkmSsfl3a8uCUCfy/w3bDgL++uoiTbnuHL9ZrtHCRbKSEJBktLzeH/5m8Ow9O2Z8+nYqYtWwTx/z5DZ6apTFzRbKNEpK0CPsO7sqzP53IUXv2pmRnBT99cBYXPfwRJTsrog5NRFJECUlajM5tC7j1zL34zQmjKMrP4bEPlvPNP7+hIYdEsoQSkrQoZsbp+w3g3/99EMN7d2DJ+lJOvPVtbn9tEVVVemZJpCVTQpIWabdeHXjyJwfyvQMGUV7p/Pa5+Zxz9/uaHl2kBVNCkharKD+Xq4/fkzvP3psubfN549NiDr1pGlc+OYfVm3dEHZ6IJEgJSVq8w0f04vkLD+bbY/tQUeXc9+5SDr7xVa59eh7FJTujDk9EGkkJSbJCr45F/OnUcbx44cEcM6o3ZRVV3PXmYg6+4VVueH4+m0o1/JBIplNCkqyyW68O/N8Z43nm/x3E4Xv0pLSskv+btoiJv3uVm//zKVt3lEcdoojUQQlJstKefTpx5zn78MSPD+CgXbuzdWcFf/zPQibe8Cq3vbaI0jI9vySSaZSQJKuNG9CF+8/bjwenTGCfQV3YVFrO9c/N5+AbpnH3W4vZUV4ZdYgiElJCklZhwpBuPHz+/tz7/X0Z068TxSU7+dW/53HI76dx15uLWbdVnR9EopYXdQAizcXMOHj3HkzcrTv/+WQtN724gPmrt3Lt0/O47pl5HLhrd741ti9H7tmLDkX5UYcr0uooIUmrY2ZMHtGLw4b35MV5q3l05gpeW7iWNz4t5o1Pi7nsiRwO36Mnx4/pS06lRn8QaS5KSNJq5eQYR43chaNG7sKm0jKem7OaJz9cwftLNvDsx6t59uPVtMmD4zZ+xLfG9mXCkG7k5ljUYYtkLSUkEYKBW0/bdwCn7TuAlZu28/TslTw1ayVzV27h4RnLeXjGcnp2KOS4MX341tg+jOrbCTMlJ5FUUkISidOncxumHDyUKQcP5YGnX2F1QV+e+mglS9eXctebi7nrzcUM7NaWCYO7MX5QF/Ye2IXB3dspQYkkSQlJpB592udw+qRh/Gzy7sxatomnZq3k6dlBclq6vpSHZiwDoFu7AvYaGCSnvQd1ZWTfjhTm5UYcvUjLooQk0ghmxrgBXRg3oAtXHLsHH6/YzMylG5m+ZAMzl26kuKSMl+at4aV5awAoyMthTL9OjB/Ylb0HdmH8wC50aVcQ8acQyWxKSCIJysvN+TI5nTdxCO7O0vWlzFi6kZlLNzBjyUY+XVvC9CUbmb5k45f7De3RjuG7dGRwt3YM6t6OweHSpW2+LveJkGUJycx6An8E9g6LPgYudPfl0UUl2c7MGNQ9SDLfHd8PgE2lZXzwRZCQZi7ZyEfLN7Fo3TYWrdv2tf07FuV9mZyqE9WgMGl1aqPnoaT1yJqEZGYFwEvAQmBPwIG/Aa+a2Th3L4kyPmldOrct4NDhvTh0eC8Ayiqq+GTVFhatK2FJ8TYWry8N1sXb2LKjgo+Wb+aj5Zu/dpxu7Qro16UN3dsX0q19QbgupPuXPwfrLm0L1CVdWrysSUjAOcBo4AR3rwAws0uAFcB/ATdGGJu0cgV5OYzp35kx/Tt/pdzdKS4pY8n6IDktKd7GkvXb+HzdNpauL2X9tjLWb2t46owcg67tCujWLkhS3doX0r4wj/aFubQrzKN9YR7tqpeC+LJcOhTmU5Sfo0uHEqlsSkgnAl+4++fVBe6+2szmhduUkCTjmBk9OhTSo0Mh+wzq+pVt7s6aLTtZuXk7xVt3sn5b2ZfrdSU7v/x5fclONpaWU1xSRnFJGaxpWiw5FszCW5CXQ2FeDgV5ORTk5lCQF5RtL9nOXYveC8tit+eQl2Pk5uSQl2vk5lj42sg1Izc3eJ1jYXluTlCeE3z+HDNyDHLMsHBdXWZfKQMzMIzwv2A7NeUWlhP3ujrRVufb2LT71RxsdZR/dZ/PN1XSedmmWtuxrpSejlxvdb5by5RNCWk0weW6eIuBw5o5FpGkmRm9OxXRu1NRg3XLK6vYGCaq9SVlbNhWxtadFWwLl5Ivf6788ueSnRWUllV++fPOiipKyyopLat7BPSFG4tT+RFbtnffijqCrJNNCak7MLOW8i1AWzNr4+7bYzeY2RRgSvhyp5nNSXOMLUV3QN88AbVFDbVFDbVFjWGpOlA2JaS61HlO6+5TgakAZjbD3feuq25roraoobaoobaoobaoYWYzUnWsbJoPqRjoUEt5B6A0/uxIREQySzYlpNnAoFrKBxM8jyQiIhksmxLS48BAMxtUXWBmvYA9gMcasf/U9ITVIqktaqgtaqgtaqgtaqSsLcw9OyYgCx+MnQF8ApwBVAF3AQcBejBWRCTDZc0ZkruXAZOBSmAeQWLqCByqZCQikvmy5gxJpLmZ2a+By4Fz3f2eiMMRSSkz2wW4GzjS3ZvlCdysOUOqjZn1NLN/mNmCcHnUzPo1ct98M7vWzOab2Rwze9vMDkp3zOnS1LYws13M7Fdm9r6ZfRi2x+NmNqo54k6HZH4vYo7RD/ifNIXYbJJtCzMbY2ZPmdkH4e/GAjO7IZ0xp0uS3xe7mNmdYRvMNrO5ZnaZmbXI0XHN7ATgHWBoE/e/0MzmhW3xgZl9u1E7untWLkAB8BHwCMHzVrnA34FPgfaN2P82gpEfeoSvzwO2A2Oj/mzN2RYx7dA/fF0UHqcUGBX1Z2vu34uY49wLPE0wiO/3ov5cUbQFcACwEjgwpuwnwJKoP1tztgXBH/YfAnOAbmHZuPD74vdRf7Ymtsd7wG7APUGaSGjfXxA8hjM0fD0ZKAeObnDfqD94Ghv0h+GXxZCYst4E95gubmDfYQSdIr4fVz4XeCbqz9bMbXEbcF5c2dDweLdE/dmasy1i6u8FLAKObOEJKZnfCyO4T3txXHl+Y754Mm1Jsi1GhPv+LK78KWBV1J+tie2RF64TSkhAZ2AbcE1c+TPA3Ib2z+ZLdrUOtkrQ4eHEBvY9geAf3Ktx5a8AR5hZ+1QG2gySaYv/JpjGI9bKcN0lZRE2n2TaotofCO4d7Ux9eM0qmbY4CBhOcJb4JXcvd/fnUh1oM0imLSrCdfzIN9VnWi2OhzMmNMFRQFtq/+4cYWbD69s5mxPSaIKBVeMtBhq6/zGa4Azpi1r2zSP4i6glaXJbuHuFu1fFFe8erqclH1qzS+b3gvBaeBvgoRTHFYVk2uKAcN0pvIc0N7xf8Gsza5PSKJtHMv9GFgIPAOdXPwdpZocSXKq6JaVRZr7R4Tq+LRfHba9VNiek7sDWWsq/HGy1gX1L3T1+2OMt4bpbCuJrTsm0RW2mEFy+vC/ZwCLQ5LYIb1D/DrjIw+sQLVwyvxf9w/U/gevcfU/gTOB7BJeqWppk/42cAzwLfGpmK4EnCWarvja1YWa87uE6vi0b9d2ZzQmpLsl0X8yuyUea8HnCv/xOAU5295Z+ySpWY9rivwiug7+Z7mAi1pi2qJ4T4y53fx/A3WcTJOzJZvaNdAXXzBpsCzMrJLhEtS8wyN37AJOAS83s8vSG12I06rsmmxNSMoOtFhP8VRR//bf6eOtTEF9zSsnAs2Y2hqB32fHuPi+F8TWnJrWFmXUGLgUuSWNszS2Z34vqv4BnxZV/GK73STK25pZMW/yA4J7axe6+AsDdPwB+D1xrZmNTHWwGq56SI74tG/Xdmc0JKZnBVmcTtE3/uPLBBDcwP0k2uGaW9MCzZjaa4DLEqe7+dupCa3ZNbYsJBP/vHzGzWWY2C7gz3HZNWHZVSiNNv2R+L+aH6/jvkMo6yjNdMm1RfY/p07jyhQRnBi0tOSdjdrgeFFc+OG57rVraL00iGj3Yqpn1MrPYtniCoBvnpLhjHgK86O61XWvOZMm0RXUyego4q/pyVfgg4O1pjjsdmtQW7v68u/d397HVC8GzaQBXhWXXNMsnSJ1kfi+eJUg+8TepR4br6akONs2SaYu14XpA3DEHhuuWdkWl0cysWziOaLXnCZ5RnBRX9RBgnrvPpz5R93dPYz/6AoJs/BBBz7gcgmEwvvKgG3AgwT+sW+P2vw1YAHQPX3+flv1gbJPaguCvv3XArQQ3rauXC4FpUX+25v69iDvWJFr2c0jJ/hv5A7AK2C183Tfc98WoP1tztgXBX/9bgBeBDmHZAOAzgufV2kT9+ZJol3uo4zmk8HPvAJ6LK/9F+J0xJHx9OI18MDZrZ4x19zIzmwz8keBZAid4kjp+sNUSYDPBP6xYFwC/BN4ys3KCa+ZHuHv8NfOMl2Rb/Iqg58yPwiXWa2kLOk1S8HuBmfUk+PKpfh7tGjO7kOAB4pTNnpluKWiLiwnuGTxrZpUED8U+RvDvpkVJpi3cfbGZ7QtcDUw3szKCtngBuNZb4OSgZnYjQbf1AeHr6u+9fT0YyBqCP9A3UPNcIgDufr2Z7QCeNrMKggR+kjfi+TQNrioiIhkhm+8hiYhIC6KEJCIiGUEJSUREMoISkoiIZAQlJBERyQhKSCIikhGUkEREJCMoIYmISEZQQhIRzOxqM/PqgWPNLH5g4aYcs3/M8dzMrk5BqJLFsnboIGm9Yob2GUAwzfpHBKMuFxKMsfUkcHvckDCY2e+Aw9x97wTeqzPBuH5PtsRhpeJ5MGhsqo61DBgLYGYaEkYapDMkyTruvjb8Yv1X+Hqsu49x9+HAT4FvAx+Z2fC4Xdfy9WnrG9KZYOy21jTnjUhaKCFJq+LBxGmHAZsIBgVtG7PtJnf/TmTBibRySkjS6oSjFV9FMHz+DwDM7K9m9kV4r2NQdV0z29/MXjOzD83sIzN7zsxOCLedQDAvENRM0jcrvIyHmV1mZu+b2Uwz+9jM/mlmfWOOvXdYv8zM7jGzi83sHTNbaWa3x80zg5l1NrNbzWypmc0Olz+a2eCYOl3N7I6wzsLw/Y9ualuF8W0wsyVmdoyZvWpmq83sCTPraGYHmtnzZrbCzB4xs05NfS+RyOfb0KIlXQv1z+XShmAG2Gdjyr5HMO3AoPB1B2AjcEb42oAbiJkHimBmzFrnRCI4Cxsd/pxLMLXBB0BuXL0lBNMZHBe+HhHGNiWmTgEwA3idcH4eYDeCy4wXhq8Lw+O/C3QMy74bHuuQBtrq6nra6h6CKReuDF/3CtvlfoJpuwF6h3V+XccxHLg66t8JLZm96AxJWiUP5qgppmZWz9oMI7hHtDjcxwmSyqONfJsJ7j473LcSuB0YB9TWaWKNu/87rDuPYIrwSTHbzwLGA5d52BnD3T8FphIknOo644Ar3H1LWOdRgkSW7BxF7YFbwmOuAd4ETgPuCMtWA28QzAwq0iTqZSetmTWwfQGwBnjSzG4BHnL3hcBfGnn8Lmb2JLArQdKovgQ3BHgvru7CuNcbCM5Eqk0O11+ZGtzdr4h5eTjBmcjbcceaA5xtZvnuXt7I2OOtd/dNcfHFl60nmPJbpEl0hiStUtiZoRuwtK467r4V2A94Avg5sCC8J3NgI44/BniVYArrsR70+jsm3FxYyy6lca+rCC7zVesOlLr7znretjthQoq5nzWL4KxlA0EX+KaKj8/rKMtFpIl0hiSt1WSCL89n66vk7kuB88Mpyk8ErgOeM7NB7r6hnl1PIUg817l7RT31GqsYaGtmhfUkpWKCRDY+vEQo0qLoDElaHTMrAq4huDd0Vz31RpnZZRDcc3L3+4GfEXR2GBRWq74EZuE+481sd2rOgmIfCO2dRNgvhevxcTFeYmYXxdTJI+gUEVtnnJndnsR7izQLJSRpVcxsb+AVgpv0x4SdG+rSDbjIzHYL9zXgAGA18ElYZw2wHegXvr4ZmAA8E77+n3DfAuCSJEK/D5gJXGdm7cJjjiQYJeLFuDo3mVmHsE5Xgs4IC5J4b5HmEXU3Py1aUr0APYFZBPdNPPx5FkHPtTcJ7gd1iNvnrwSjNDgwDziP4J7MTcDH4f5zgaeBUXH7nk9wL2oOwf2morB8CvAZQTJ4OXxfD9/nemBoeNyyMNYnwv3eBkrCZRZQEJZ3Bm4L32tW+FkOjYulus4XBEMmzQB+0og2u5paun0T3AfbEMY4iyBJP9GIsqFxx1G3by0NLuauIaZEWrtw4NNfuntDPQ+benwHfuXuV6fj+JIddMlORCA4G1uTrtG+CS5tljS0j7RuOkMSEZGMoDMkERHJCEpIIiKSEZSQREQkIyghiYhIRlBCEhGRjKCEJCIiGUEJSUREMsL/B8liaD7cd1STAAAAAElFTkSuQmCC\n",
"text/plain": [
"