{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Boundary conditions tutorial"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This tutorial aims to demonstrate how users can implement various boundary conditions in Devito, building on concepts introduced in previous tutorials. Over the course of this notebook we will go over the implementation of both free surface boundary conditions and perfectly-matched layers (PMLs) in the context of the first-order acoustic wave equation. This tutorial is based on a simplified version of the method outlined in Liu and Tao's 1997 paper (https://doi.org/10.1121/1.419657).\n",
"\n",
"We will set up our domain with PMLs along the left, right, and bottom edges, and free surface boundaries at the top as shown below.\n",
"\n",
"\n",
"\n",
"Note that whilst in practice we would want the PML tapers to overlap in the corners, this requires additional subdomains. As such, they are omitted for simplicity.\n",
"\n",
"As always, we will begin by specifying some parameters for our `Grid`:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"shape = (101, 101)\n",
"extent = (2000., 2000.)\n",
"nbpml = 10 # Number of PMLs on each side"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will need to use subdomains to accomodate the modified equations in the PML regions. As `Grid` objects cannot have subdomains added retroactively, we must define our subdomains beforehand."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"from devito import SubDomain\n",
"\n",
"class MainDomain(SubDomain): # Main section with no damping\n",
" name = 'main'\n",
" def __init__(self, PMLS):\n",
" super().__init__()\n",
" self.PMLS = PMLS\n",
" \n",
"\n",
" def define(self, dimensions):\n",
" x, y = dimensions\n",
" return {x: ('middle', self.PMLS, self.PMLS), y: ('middle', 0, self.PMLS)}\n",
"\n",
"\n",
"class Left(SubDomain): # Left PML region\n",
" name = 'left'\n",
" def __init__(self, PMLS):\n",
" super().__init__()\n",
" self.PMLS = PMLS\n",
"\n",
" def define(self, dimensions):\n",
" x, y = dimensions\n",
" return {x: ('left', self.PMLS), y: y}\n",
"\n",
"\n",
"class Right(SubDomain): # Right PML region\n",
" name = 'right'\n",
" def __init__(self, PMLS):\n",
" super().__init__()\n",
" self.PMLS = PMLS\n",
"\n",
" def define(self, dimensions):\n",
" x, y = dimensions\n",
" return {x: ('right', self.PMLS), y: y}\n",
" \n",
" \n",
"class Base(SubDomain): # Base PML region\n",
" name = 'base'\n",
" def __init__(self, PMLS):\n",
" super().__init__()\n",
" self.PMLS = PMLS\n",
"\n",
" def define(self, dimensions):\n",
" x, y = dimensions\n",
" return {x: ('middle', self.PMLS, self.PMLS), y: ('right', self.PMLS)}\n",
"\n",
"\n",
"main_domain = MainDomain(nbpml)\n",
"left = Left(nbpml)\n",
"right = Right(nbpml)\n",
"base = Base(nbpml)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And create the grid, attaching our subdomains:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from devito import Grid\n",
"\n",
"grid = Grid(shape=shape, extent=extent,\n",
" subdomains=(main_domain, left, right, base))\n",
"x, y = grid.dimensions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can then begin to specify our problem starting with some parameters."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"density = 1. # 1000kg/m^3\n",
"velocity = 4. # km/s\n",
"gamma = 0.0002 # Absorption coefficient"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We also need a `TimeFunction` object for each of our wavefields. As particle velocity is a vector, we will choose a `VectorTimeFunction` object to encapsulate it."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from devito import TimeFunction, VectorTimeFunction, NODE\n",
"\n",
"p = TimeFunction(name='p', grid=grid, time_order=1,\n",
" space_order=6, staggered=NODE)\n",
"v = VectorTimeFunction(name='v', grid=grid, time_order=1,\n",
" space_order=6)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A `VectorTimeFunction` is near identical in function to a standard `TimeFunction`, albeit with a field for each grid dimension. The fields associated with each component can be accessed as follows:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[[0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" ...\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]]\n",
"\n",
" [[0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" ...\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]\n",
" [0. 0. 0. ... 0. 0. 0.]]]\n"
]
}
],
"source": [
"print(v[0].data) # Print the data attribute associated with the x component of v"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may have also noticed the keyword `staggered` in the arguments when we created these functions. As one might expect, these are used for specifying where derivatives should be evaluated relative to the grid, as required for implementing formulations such as the first-order acoustic wave equation or P-SV elastic. Passing a function `staggered=NODE` specifies that its derivatives should be evaluated at the node. One can also pass `staggered=x` or `staggered=y` to stagger the grid by half a spacing in those respective directions. Additionally, a tuple of dimensions can be passed to stagger in multiple directions (e.g. `staggered=(x, y)`). `VectorTimeFunction` objects have their associated grids staggered by default.\n",
"\n",
"We will also need to define a field for integrating pressure over time:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"p_i = TimeFunction(name='p_i', grid=grid, time_order=1,\n",
" space_order=1, staggered=NODE)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we prepare the source term:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from examples.seismic import TimeAxis, RickerSource\n",
"\n",
"t0 = 0. # Simulation starts at t=0\n",
"tn = 400. # Simulation length in ms\n",
"dt = 1e2*(1. / np.sqrt(2.)) / 60. # Time step\n",
"\n",
"time_range = TimeAxis(start=t0, stop=tn, step=dt)\n",
"\n",
"f0 = 0.02\n",
"src = RickerSource(name='src', grid=grid, f0=f0,\n",
" npoint=1, time_range=time_range)\n",
"\n",
"# Position source centrally in all dimensions\n",
"src.coordinates.data[0, :] = 1000."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"tags": [
"nbval-skip"
]
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAGBCAYAAADYEOPMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcZGV59//P1fs6+wIDs7LJAIJmROIGogkaDZhgNCpuiWJiTDSGJ5GfCTFIFOMSf2pUMHncMIlKMBBFRdlUBA2yybDNDLMwbNOz9VrdXV19PX+cc3qKmurqU9W1nvq+X696dfepc07dZ6qm6qrrvu77NndHREREpFgttW6AiIiINCYFESIiIlISBREiIiJSEgURIiIiUhIFESIiIlISBREiIiJSEgURIiIiUhIFESIiIlISBREiIiJSEgURIiIiUpK2Wjeg3i1btszXrVtX62aIiIhUxa9+9au97r48zr4KIuawbt067rzzzlo3Q0REpCrMbGfcfdWdISIiIiVRECEiIiIlURAhIiIiJVEQISIiIiVRECEiIiIlURAhIiIiJWmIIMLMjjazz5rZ7WY2ZmZuZutiHttiZheb2Q4zGzeze83s/Mq2WEREJPkaIogAjgVeBxwAflrksR8GPgR8DnglcAfwbTP7nXI2UEREpNk0ymRTP3H3lQBm9g7gt+McZGYrgIuAy939E+Hmm83sWOBy4PpKNFZERKQZNEQmwt2nSzz0HKADuCpn+1XAKWa2fl4NExERaWINEUTMw0nABLA1Z/vm8OfG6jZHREQkOZIeRCwBDrq752zfn3W/FCkz7dz/+CCZ6dx/VhERaSZJDyJKYmYXmtmdZnbnwMBArZtTd75z9+O8+rM/4+xP3sITB1O1bo6IiNRI0oOIA8AiM7Oc7VEGYj95uPuV7r7J3TctXx5rNdSmcteuA7S3Gjv3jXHjg0/XujkiIlIjSQ8iNgOdwDE526NaiAeq25xk2PzEEJvWLqG/s41Hnh6pdXNERKRGkh5E/ABIA2/K2X4BcL+7b69+kxrbVGaah54c4qRVCzhuZR9b9gzXukkiIlIjjTJPBGb22vDX3wh/vtLMBoABd7813GcK+Kq7/zGAu+8xs08BF5vZMHAX8HrgbODcql5AQmwbGGViapqTjlrA8PgUP1Z3hohI02qYIAL4ds7fnw9/3gqcFf7eGt6yfRAYAd4LHAE8DLzO3b9bmWYm2+YnBgE4edVC9o1M8s07H2PfyARL+zpr3DIREam2hgki3D23ODLWPu6eAS4LbzJPDz45RGdbCxuW9/Hk4DgAW/aMKIgQEWlCSa+JkDJ7YnCcVYu6aW0xjlvZB8CWp1UXISLSjBRESFEGhidYHmYdjljQRWdbC7v2j9W4VSIiUgsKIqQoe4cnWL4gCCLMjGV9newbmaxxq0REpBYUREhRsjMRAMv6OxkYmahhi0REpFYUREhsqckMwxNTrFhwKIhY3tfBXmUiRESakoIIiW1gOMg4ZGcilvZ2sk+ZCBGRpqQgQmIbGAmGdC7vz+7O6GDf6CTTWtFTRKTpKIiQ2PYMBRmHFf1dM9uW9XWSmXYOptK1apaIiNSIggiJLSqgzM5ERJNMqUtDRKT5KIiQ2AaGJ2htMZb0dsxsW9YX/K4RGiIizUdBhMQ2MDzB0t4OWlsOzS4eFVlqhIaISPNRECGx7RmeeEZXBqg7Q0SkmSmIkNjyrda5qLud1hZjr4IIEZGmoyBCYhsan2JRd/sztrW0GEt7O9g7rO4MEZFmoyBCYhtMpVnQffjq8Uv7Otk3qkyEiEizURAhsbg7g6k0C3MyERB0aQylpmrQKhERqSUFERLL6GSGzLTnDSIWdLcxqMmmRESajoIIiWUoDBIWdOUJIrraGRpXECEi0mwUREgsUaYhfyaifSbIEBGR5qEgQmIpGER0tTM6mWEqM13tZomISA0piJBYoiBiwSw1EQDD4yquFBFpJgoiJJahOTIRgOoiRESajIIIiaVwJiIMIjTMU0SkqSiIkFiGUmnMoL/z8MmmFnQF25SJEBFpLgoiJJbBVJoFXe20ZK3gGVnYE2UiFESIiDQTBRESy9D4VN4pr0E1ESIizUpBhMQy25TXoJoIEZFmpSBCYikURPR2tNJiaOprEZEm0xBBhJmtNrOrzWzQzIbM7BozWxPz2DVm9lUz22VmKTN7xMwuM7PeSrc7SQoFEWYWzFqp7gwRkaaSv5O7jphZD3ATMAG8FXDgMuBmM3u2u48WOLYX+DHQDvwdsAt4HvAPwHHA6yvb+uQYCgsrZ7OgS1Nfi4g0m7oPIoB3AhuAE9x9K4CZ3QdsAd4FfKrAsS8kCBbOcfcbwm03m9kS4CIz63H3sco1PTmGxtN554iILOhuY0gzVoqINJVG6M44F7gjCiAA3H07cBtw3hzHdoQ/h3K2HyS49sPHK8ph0plpxtPT9OWZIyKiTISISPNphCDiJOD+PNs3AxvnOPbHBBmLj5nZRjPrM7OzgfcCXyzUFSKHjE4EGYZCQcRC1USIiDSdRggilgAH8mzfDywudKC7jwMvIrjOzcAwcCPwXeA95W1mco3ECCL6u9o0OkNEpMk0Qk1EycysC/gmsAJ4M0Fh5enAJcAU8KezHHchcCHAmjWxBoEk2uhEBoDeAkFEX2f7zH4iItIcGiGIOED+jMNsGYpsfwycBRzr7tvCbT8xs0HgSjP7orvfm3uQu18JXAmwadMmL7XhSTEyEWQY+rpmf7n0drYyOjmFu2OmUhMRkWbQCN0ZmwnqInJtBB6Y49hTgANZAUTkl+HPE+fZtqYwEmYY+jpbZ92nt7MNdxibVDZCRKRZNEIQcR1whpltiDaY2TqC4ZvXzXHsU8BiMzs2Z/vzw5+Pl6mNiRYVVhbqzojui/YVEZHka4Qg4kvADuBaMzvPzM4FrgUeA66IdjKztWY2ZWaXZB37FYJiyuvN7K1m9lIz+z/AJ4BfEQwTlTmMhPM/9HYUqokIshQjCiJERJpG3QcR4TDMs4FHgK8D3wC2A2e7+0jWrga0knVN7r4DOAO4h2CWy+sJJq+6Evgtd5+uwiU0vCgw6C9UE9ERZSLUnSEi0iwaobASd98FnD/HPjvIM3mUuz8AvK4yLWsOcbozouGfykSIiDSPus9ESO2NTE7R0dZCe+vsL5cowBibVBAhItIsFETInEbGp+gvkIWAQ0GEMhEiIs1DQYTMaXRiqmBXBhzqzlBNhIhI81AQIXMamcjMGUT0hqMzNMRTRKR5KIiQOY1MpAtONAWHRmeoO0NEpHkoiJA5jU5kCi6+BdDSYvR0tCoTISLSRBREyJzi1ERAUFw5qtEZIiJNQ0GEzGl4YmrOTAQExZUjKqwUEWkaCiJkTqMxg4jeTnVniIg0EwURUlBm2hmbnHt0BgTFlSqsFBFpHgoipKCoxiFud4YyESIizUNBhBQUZ92MSK+CCBGRpqIgQgo6FEQUnici2EeFlSIizURBhBQUTWMdTSZVSK/miRARaSoKIqSgsckgiOiJmYlIpTNkpr3SzRIRkTqgIEIKipb27omRiZhZhEsTTomINAUFEVJQlIno7YiXiQAYU12EiEhTUBAhBUWZiO4YQURPuM+YMhEiIk1BQYQUdCgTMXd3RvdMEKFMhIhIM1AQIQUVU1jZ3R7sk0oriBARaQYKIqSgsckpWluMjta5Xyo9ykSIiDQVBRFS0OhEhp6OVsxszn2j7oyUaiJERJqCgggpKDWZmckwzCUaBqpMhIhIc1AQIQWNTk7FKqoEdWeIiDQbBRFS0NhkJtbwTsjuzlAQISLSDBRESEFjxWQi2pWJEBFpJgoipKBiMhFtrS10tLYwllZhpYhIM1AQIQWNTWZiLQMe6e5oVXeGiEiTaJggwsxWm9nVZjZoZkNmdo2ZrSni+BPN7NtmttfMUmb2sJm9t5JtToKxiSm62+N1Z0BQXKnuDBGR5hD/06GGzKwHuAmYAN4KOHAZcLOZPdvdR+c4flN4/C3AO4BB4Digr4LNToSxtDIRIiKSX0MEEcA7gQ3ACe6+FcDM7gO2AO8CPjXbgWbWAnwNuNHdfy/rrpsr19zkGJvIxFoGPNLT0appr0VEmkSjdGecC9wRBRAA7r4duA04b45jzwJOpECgIfmlM9NMZqZjTzYF0NPeplU8RUSaRKMEEScB9+fZvhnYOMexLwp/dpnZHWaWNrM9ZvYZM+suaysTZmbxrSKCCHVniIg0j0YJIpYAB/Js3w8snuPYVeHPbwI3AL8F/BNBbcS/l6uBSZSaCSJUWCkiIodrlJqI+YgCpavc/ZLw91vMrBW43MxOdPcHsw8wswuBCwHWrIk9ACRxRsNuiWILKxVEiIg0h0bJRBwgf8ZhtgxFtn3hzx/lbL8h/Pmc3APc/Up33+Tum5YvX15UQ5NkbCIIBrrbiwgi2lVYKSLSLBoliNhMUBeRayPwQIxjC5kuqUVNYGwmE1Fsd4YKK0VEmkGjBBHXAWeY2YZog5mtA14Y3lfI9wnmlzgnZ/srwp93lqeJyRN1S8Sd9jrYt43x9DTT016pZomISJ1olCDiS8AO4FozO8/MzgWuBR4Droh2MrO1ZjZlZlHtA+6+D/go8Cdm9hEze7mZfQC4BPhq9rBReaZSRmdE+6pLQ0Qk+RqisNLdR83sbOCfga8DBtwIvM/dR7J2NaCVw4OjS4Fh4N3ARcCTwMeBD1e46Q1tpjujyNEZwbGZorpBRESk8TTMu7y77wLOn2OfHQSBRO52J5hsShNOFSHKJnQVWVgJaK4IEZEm0CjdGVIDqZK6M4K4VMuBi4gkn4IImVUpmYjs7gwREUk2BREyq9Rkhs62FlpbDushmlU0kkPdGSIiyacgQmaVSmeKGt4JykSIiDQTBREyq7HJDD1FdGWAhniKiDQTBREyq1Q6Q1eRmYjusLAypVkrRUQST0GEzCo1mSlqZAYwk7lQd4aISPIpiJBZpSYzRS2+BYcKKxVEiIgkn4IImdVYOjPTPRFXZ1sLLabRGSIizUBBhMxqfDJDd3txLxEzo7u9VZkIEZEmoCBCZjWWnpqZgbIY3R1tpDRjpYhI4imIkFmlJqeLmq0y0tOhTISISDNQECGzSk1OFV1YCQoiRESahYIIycvdSaWLH+IJwQgNFVaKiCSfggjJa2Jqmmmn6GmvIcpEqCZCRCTpFERIXuPhtNWldGd0t7epO0NEpAkoiJC8orUvSs1EaO0MEZHkUxAheUWZhFJqIlRYKSLSHBRESF5RYWQpQzxVWCki0hwUREheUXdEqZmIVDqDu5e7WSIiUkcUREheUSahtHki2shMO5OZ6XI3S0RE6kjsIMIC55rZJ8zsy2a2Ntx+ppmtqlwTpRaimoZSCiujwENdGiIiyRZrYQQzWwxcDzwfGAb6gM8CO4F3AvuBv6hQG6UG5jPEsydrOfBFPWVtloiI1JG4mYiPA6uBFwJLAcu678fAy8rcLqmxQ6MzSlmA61AQISIiyRX3E+I84CJ3v93Mcr+a7iIIMCRBUvPKRAQvK3VniIgkW9xMRB/w+Cz3dfHMzIQkQCqctno+NRGa+lpEJNniBhEPA789y31nAr8uT3OkXqTSGVpbjPbW4uPDme4MzVopIpJocbszPg98zswGgX8Pty0ys7cD7wEurETjpHbGJjN0t7diVnwQERVWqjtDRCTZYmUi3P1K4FPAPwBbw80/Aq4EPu3u36hM8w4xs9VmdrWZDZrZkJldY2ZrSjjPB8zMzexnlWhnUoynMyV1ZcAzR2eIiEhyxS69d/cPmNkXgN8CVgD7gB+5+6OValzEzHqAm4AJ4K2AA5cBN5vZs919NOZ5NgB/C+ypVFuTIspElKJ7JhOhmggRkSQravyeu+8E/rVCbSnkncAG4AR33wpgZvcBW4B3EWRJ4vgC8A3gBIq89maTmsyUNOU1HBqdoUyEiEiyzfpBWmxXgbvvmn9zZnUucEcUQISPt93MbiMYfjpnEGFmbwSeC7wBuKZSDU2KVDpT0uJbAD3RjJUqrBQRSbRC38Z3EHQbxFXaJ048JwHX5tm+GfiDuQ4OZ9z8Z+Cv3X1/KcWCzWY+mYiWFqOzrUWFlSIiCVcoiPgjDgURnQS1BEPAt4CngSOA1wH9wIcr2EaAJcCBPNv3A4tjHP9x4BHgK2VsU6Kl0hkWdreXfHx3uJKniIgk16xBhLt/JfrdzD4N3AX8nmet72xmlwL/DWysYBvnxcxeDLwFeK7HXJvazC4kHLa6Zk3RA0ASITVZ+ugMCCacUiZCRCTZ4k429QbgitwP4fDvLwJvLHfDchwgf8ZhtgxFtiuAfwN2m9kiM1tEEDy1hn935h7g7le6+yZ337R8+fL5tr0hpdKlj86AIBOhyaZERJIt7giFPmC2T9MVQG95mjOrzQR1Ebk2Ag/MceyJ4e1P8tx3APhL4NPzal0Cjc2jJgKCTMS4MhEiIokWN4i4BfiImT3o7v8bbTSz04F/DO+vpOuAT5jZhmheCjNbR7Cq6AfmOPalebZ9mqAQ9M85NHmWZEmlM3TNI4jo6WjVEE8RkYSLG0S8h2DJ7zvM7DGCwsqVBKt3bg/vr6QvhY9xrZn9LUHB54eBxwi6KwAws7XANuBSd78UwN1vyT2ZmR0E2vLdJ5CZdianpulpL30qja72VobHNdmUiEiSxZ32ejvwLIIugRsJZqu8kWCipxPdfUelGhg+/ihwNsEIi68TTBi1HTjb3UeydjWCDEPcWg/JY2YZ8I7S/xm721sZV02EiEiiFTPtdZogI/ClyjWn4OPvAs6fY58dxFiW3N3PKk+rkmlsZhnw0jMR6s4QEUk+fWOXw4xPTgPMe3SG5okQEUm2WF81zWw7hWevdHc/pjxNklobSweZiPmNzmjTPBEiIgkXN199K4cHEUuBFwAjBCtsSkJEH/7zy0S0kEpncHc0zbiISDLFCiLc/W35tocTN/2AYOSGJEQURJS6ABcEAUhm2pnMTNPZVsllVUREpFbmVRPh7gcJ1qW4pDzNkXoQ1TLMqzsjLMqM6itERCR5ylFYOQ4cXYbzSJ2IRlXMd+0M0HLgIiJJVvIYPjNrA04GPkQwLbUkxMw8EfPozoiyGNFwURERSZ64ozOmmX10xhDwqrK1SGpuPD3/TESXMhEiIokXNxNxKYcHEePATuD77j5Y1lZJTUXdGfOpiYiO1TBPEZHkijs640MVbofUkZnRGfMYVRFlMZSJEBFJrliFlWZ2k5k9a5b7jjczzRORIKl0hq72FlpaSp/fYaawUpkIEZHEijs64yxgwSz39QNnlqU1UhdSk5l5FVWCMhEiIs2gmCGesxVWHkMwa6UkxNhkhp55LL4FykSIiDSDWT8pzOztwNvDPx240syGc3brJhjmeWNlmie1MB52Z8zHoSGeCiJERJKq0CfFNJAJb5bzd3TbB3wB+OPKNlOqaWxyat6ZCA3xFBFJvlk/Kdz9q8BXAczsZuBP3f2hajVMaieVnn9NRGdbC2bqzhARSbK4QzxfWumGSP1ITWZY1NMxr3OYGT3trcpEiIgkWKGaiLcA33P3feHvBbn718raMqmZsckMqxbNf+XN7g4FESIiSVYoE/EV4AyCuoevzHEeBxREJMRYGYZ4QhhEqDtDRCSxCgUR64Ens36XJjGezsxr3YxId7uCCBGRJCtUWLkz3++SfME8EeXIRLQxpu4MEZHEmt9kAJI409MejM6Y5xBPgO72FsaViRARSaxChZXbmX2Wylzu7seUp0lSS+NT4TLg5aiJaG9l78jkvM8jIiL1qdDXzVuJH0RIQpRjGfBIT0cbY5Nj8z6PiIjUp0I1EW+rYjukTkSFkOUorOxqb2U8PT3v84iISH1STYQ8QzSvQ3kyEa2MTU7N+zwiIlKfYgcRZnacmX3VzB4xs9Hw51fM7NhKNlCqK+rOKNs8ERqdISKSWLFK8M3sLOB6IAV8D3gaWAn8LvB6M3uFu99aqUZK9USZg3J2Z0xPOy0tNu/ziYhIfYmbifgkcDew1t3f4u7/x93fAqwD7gnvrygzW21mV5vZoJkNmdk1ZrYmxnGbzOxKM3vIzMbMbJeZfcPMNIFWHuMz3RnzH+IZdYlEIz5ERCRZ4gYRG4GPuftI9kZ3HwY+BpxU7oZlM7Me4CbgWcBbgTcDxwE3m1nvHIf/Ydi+zwCvBD4APBe408xWV6zRDaqcozOiLhHNWikikkxxv27uBmZb1rEDeLw8zZnVO4ENwAnuvhXAzO4DtgDvAj5V4NiPuftA9gYzuw3YHp73koq0uEGVuyYiOufSeZ9NRETqTdxMxMeAfzCzVdkbzewo4O+Bj5S7YTnOBe6IAggAd98O3AacV+jA3AAi3LYTGACOKnM7G145h3hGgci4iitFRBIpbibiTGAB8KiZ3cGhwsozwt/PCosvIZi98q1lbudJwLV5tm8G/qDYk5nZicAK4MF5titxyj3EEw5lN0REJFniBhEvAqYIVvVcG97g0CqfL87atxKzXC4BDuTZvh9YXMyJzKwN+CJBJuLf5t+0ZIk+8LvaylgToUyEiEgixQoi3D1JIxk+B7wAeJW75wtMMLMLgQsB1qyZcwBIoqQmp+hqbynLkMyuDgURIiJJ1igzVh4gf8ZhtgxFXmZ2OUFw8EfufsNs+7n7le6+yd03LV++vOjGNrJgGfD5D++EQ90ZGp0hIpJMRX1ahEMiVwNdufe5+03lalQem8k/jHQj8ECcE5jZB4G/Af7c3b9exrYlSiqdKcvIDNAQTxGRpIs7Y+UG4BvA6dGm8KeHvztQnk+e/K4DPmFmG9z90bBN64AXEsz7UJCZ/QVwGfBBd/9cBdvZ8FKTmbIUVULWEE91Z4iIJFLcTMS/AmuA9wEPAZMVa1F+XwLeA1xrZn9LELR8GHgMuCLayczWAtuAS9390nDbHwKfBn4A3GRmZ2Sdd8jdY2UymsVYOYOIaIinMhEiIokUN4h4HvA2d/+vSjZmNu4+amZnA/8MfJ0g+3Ej8L6cWTSNICOSXevxinD7K8JbtluBsyrU7IaUmszQVebuDA3xFBFJpmJmrKx29uEZ3H0XcP4c++zgUFdLtO1twNsq1a6kSaUzLOubbXLS4rS1ttDR2qLRGSIiCRV3dMZHgL+JsU6FNLixyamyjc4A6Gpv0YyVIiIJFXeeiK+b2bOAHeGMlbnDKisxS6XUQGoyU5YpryM9HW0zy4uLiEiyxB2d8TbgYiBDsAJmbtdGJWaplBoYK+MQTwhGaKTS02U7n4iI1I+4eet/AL4D/LG7H6xge6TGyjk6A4LiypQyESIiiRS3JmIp8HkFEMmWmXYmp6bL2p0RZCJUEyEikkRxg4ifASdWsiFSe+VcwTMSZCIURIiIJFHc7oz3At8yswMEkzYdtl6Fu6vju8FFBZDlronYOzJRtvOJiEj9iBtEPBj+/FqBfSo57bVUQZQx6C7jEM/udnVniIgkVdxPi0vRCIzEq0R3Rk9Hq2asFBFJqLjzRHxotvvM7CzgLWVqj9TQ2EwmosyFlQoiREQSKW5h5TOY2bFmdqmZbSdYw+J15W2W1MJMd0YZayJ6w8mm3JXIEhFJmthBhJktNLMLzew24GHggwQFlu8GVlWofVJFUSairKMzOlqZdpiYUt2tiEjSFAwizKzFzH7HzL4JPAl8EVgL/Eu4y/vc/Qp3H6pwO6UKKlUTAahLQ0QkgWYNIszsk8DjwP8AryaYsfIVwBrgEnJWy5TGF80sWc7RGVEQMaYRGiIiiVPo0+IvCUZkXA+8zd33RXeYmTq4E2isAjURUUAyNqGpr0VEkqZQd8a/AcPAq4CHzexzZnZ6dZoltVCJmojeKBOh7gwRkcSZNYhw93cCRwBvAu4E3gXcbmYPAn+D5o1InPF0BjPobCtp0E5e3QoiREQSq+CnhbuPu/t/uHtUCxEtB/4BgpqIy83sAjPrqnxTpdLGJjP0tLdiVr5yl56wOyOVVneGiEjSxP7K6e5Puvs/ufvJwOkEIzSOI5gK+8kKtU+qaGwyU9aiSsgqrFQmQkQkcUrKW7v7ne7+5wTzQ5wP3FLORkltpCan6O4oX1cGZAUREwoiRESSZl5fO909TTD08zvlaY7UUtCdUe5MRDg6Y1LdGSIiSVPer53S0FLpTFnXzQDNEyEikmQKImRGajJT1uGdEIz0MNOMlSIiSaQgQmaMTWbKOtEUgJnR067lwEVEkkhBhMyoRHcGQE9nm2oiREQSSEGEzKhEdwYEdRHKRIiIJI+CCJkxNjk1M5qinLrVnSEikkgKImRGKp2hq8w1ERBkIlRYKSKSPA0TRJjZajO72swGzWzIzK4xszUxj+0ys4+b2ZNmljKz283sJZVucyNJZ6ZJZ7wi3Rm9qokQEUmkhggizKwHuAl4FvBW4M0EU27fbGa9MU7xb8A7gUuAVxNM0/1DMzutMi1uPKl0+VfwjKg7Q0QkmcrfAV4Z7wQ2ACe4+1YAM7sP2EKwuuinZjvQzE4F3gj8kbt/Odx2K7AZuBQ4t7JNbwxRd0NFRmeosFJEJJEaIhNB8EF/RxRAALj7duA24LwYx6aBb2YdOwX8J3COmXWWv7mNJ/qQr0gmoqNNQYSISAI1ShBxEnB/nu2bgY0xjt3u7mN5ju0Ajp1/8xrf6ERQs9BbgdEZQWGlaiJERJKmUbozlgAH8mzfDyyex7HR/VWxY+8oAyMTrOzvYsWCzoqMhCjVTBDRWf6XRG9HK2PpDO6OmZX9/KUam5xi94EUQ6k0w+NTDI2nmUhPk3Fn2p3paScz7Uw7wd/utW6yiMhhlvZ2cv5vHF2Tx26UIKKqzOxC4EKANWtiDQCJ5epf7eZzN8/0yLBheS+nHb2IMzYs5ZyTj2Bhd3vZHqtYle7OcIfx9HRFai7i2n1gjJsf2sOtjwxw3+5B9gxP1KwtIiLlsvHIBQoi5nCA/BmH2bIMuceuneVYOJSRmOHuVwJXAmzatKlsXz/f+Pw1nL5+CU8NjfPEwRT3Pz7ET7bs5Zq7H+eS6+7nDaev4X0vO56FPdUPJkYnK5iJ6AxX8pycqkkQ8ejACB/7wUPc8MDTuMPqJd285PjlrF/Wy+olPSzuaaevs43+rja62ltpbTFazTCzQ7+3QIvXM/ReAAAgAElEQVQZ9ZNHEREJtNQww9soQcRmgtqGXBuBB2Ic+3tm1pNTF7ERmAS25j+s/FYt6mbVou5nbHN37ts9yNdu38nXbt/Jd+97ks/84XP4zWOWVqtZAIxNBJmISgQR0SyYoxMZlvaV/fSzcne++vMdfOT6h+hsa+HPzjqW83/jaNYt7amrbhURkUbVKIWV1wFnmNmGaIOZrQNeGN5XyP8A7cAfZB3bBrweuMHda5rTNjNOXb2IT77uVK79sxeysLudt/zfX/Dd+56oajtGZgory58p6AszEaNVLK50d/7xew/yof95gBcft4wbLzqTi845gfXLehVAiIiUSaMEEV8CdgDXmtl5ZnYucC3wGHBFtJOZrTWzKTO7JNrm7ncTDO/8tJm9w8xeRjC8cz3w91W8hjmdfNRC/utPX8BzVi/mL795Dz/furdqjx3NKFmJtTMOZSKqF0T8y81b+defbedtL1jHl96yiRX9XVV7bBGRZtEQQYS7jwJnA48AXwe+AWwHznb3kaxdDWjl8Ot6O/Bl4DLge8Bq4BXufleFm160hd3tfOmtm1i/rJc//4+72TtSnUTJ6GSG9lajo638L4moi2S0SnNF/HzbXj75o0d4zWmruOTVG2lpUeZBRKQSGiKIAHD3Xe5+vrsvcPd+d3+Nu+/I2WeHu5u7fyhne8rd3+/uR7h7l7s/391vqWLzi7Kwu53PvfG5DE9M8f9d8+uqPObYxFRF6iHgUGFlNTIRY5NT/NW37mX90l7+8fdOUQAhIlJBDRNENJvjV/bzvpcfxw0PPM1tVejWGJ3MVGSiKTg0gVU1gogrbn2UJwfH+afXPrtiQZGIiAQURNSxP3rheo5a1M1Hv/8g09OVnehodGKqInNEQFZ3RoWDiD3D41zxk2286pQj2bSuanOIiYg0LQURdayrvZX3vfw47n98iFu3DFT0sUYnM/RU6Jt7FJxUuibiaz/fycTUNBedc0JFH0dERAIKIurceacdxYr+Tr5y246KPs7YxNTMUMxy62xroa3FKpqJGJuc4qpf7OS3N65k/bI4q8OLiMh8KYiocx1tLVxwxlpufWSAbQMjcx9QotHJTEWGd0IwF0allwO/7p4nODiW5h0v3jD3ziIiUhYKIhrAG05fQ2uLcc1duyv2GKMTUxWZaCrS19k2M6FVJVxz1+Mcs7yXTWvnWo9NRETKRUFEA1je38kLj13Gtfc8gVdoJcmxyamK1UQA9HS2zUxoVW6P7R/jlzv28/vPPVqzUYqIVJGCiAZx3qmr2H0gxV275lpvrDSjE5mKZiJ6O9sYmahMd8a19zwOwHmnrarI+UVEJD8FEQ3it09aSUdbC9+776mynzsz7aTSmYrOq9Db0cpYhbozbnjgaZ6zZhFHL+6pyPlFRCQ/BRENor+rnd/csJSbH95T9nOn0uEKnhUqrIRg/YxK1ETsGRrnvt2DvPzElWU/t4iIFKYgooGc/awVbN87yva9o2U9bzT0sqdCQzwhWMmzEqMzbnooCKpeduKKsp9bREQKUxDRQF56QvBBefND5c1GjM4sA954hZU3PrSHoxZ1c8LK/rKfW0REClMQ0UDWLO1hw/JebnmkvLNXRhmCStZEVGKIZ2bauePRfbzk+GUalSEiUgMKIhrMC45Zyq927GcqM122cx7KRFSuO6Ono5Xx9DSZMq4B8sATQwyPT3HGhqVlO6eIiMSnIKLBPH/9UkYnM2x+Yqhs54wyEZWcJ6IvWoSrjF0adzy6D4DfVBAhIlITCiIazPM3BKtT/mL7vrKdc7gqmYggiBgr41wRtz+6jw3Le1mxoKts5xQRkfgURDSYFf1dbFjWyy8e3V+2c46MB0FEf1d72c6Zqzcc+VGuuojMtPO/2/fz/PXKQoiI1IqCiAZ0+vol/O+O/WWbAntkIg1Af1clJ5sKMxFl6s7YNjDC8MSU1soQEakhBREN6DlrFjE0PlW2+SJGxqcwC4ofKyUa+RFlPebr7nD679PWLCrL+UREpHgKIhrQqauDD857dx8sy/mGJ6bo62yr6DDJKMsxXKbujHseO8jC7nbWL+0ty/lERKR4CiIa0HEr+unpaOXexwbLcr6R8Sn6KzgyAw6NzihfJuIgp65eREuL5ocQEakVBRENqLXFOPmohdzzWJkyEeNT9FWwHgKyMhHj6Xmfa3RiikeeHua01erKEBGpJQURDeo5qxfxwBNDTEzNf8jkSNidUUlRkFKO0Rmbnxhi2uG01QvnfS4RESmdgogGdfJRC5nMTLPl6ZF5n2t4Yoq+Cg7vBOhsa6WjraUsNREPPBF045y0SkGEiEgtKYhoUBtXLQDgwSfnP3PlyHi64jURAAu62hguQ03E5ieGWNrbwYr+zjK0SkRESqUgokGtW9pLd3srD5QjiKhCdwaEi3CVKYjYuGqBFt0SEakxBRENqrXFeNaR/TxQhjU0RsanKjrRVKS/q33ehZWTU9Ns2TM8k4kREZHaURDRwDYeuYAHnxya18yVmWlndDJT8dEZUJ7lwLfsGSadcdVDiIjUgYYIIsysxcwuNrMdZjZuZvea2fkxjltgZpeY2c/NbJ+ZHQx/f0012l1pJx65gKHxKR4/mCr5HNGHejW6M/rLUBPx0JPDAGw8sr8cTRIRkXloiCAC+DDwIeBzwCuBO4Bvm9nvzHHcGuDdwK3ABcDrgUeA75jZn1WstVVy4pFRceVwyeeIgohqdGf0lSGI2DYwQluLsVYzVYqI1FzlPznmycxWABcBl7v7J8LNN5vZscDlwPUFDt8ObHD3saxtPzSz1cDfAP9SiTZXy3Er+wDYumeE39q4sqRzRIWOfZ2VHeIJsKAMNRHbBkZYu7SH9tZGiX9FRJKrEd6JzwE6gKtytl8FnGJm62c70N1HcwKIyJ3AqvI1sTYWdLWzckEnW/eUPldEtIJnNWsi5lPDsW1glGOW95WxVSIiUqpGCCJOAiaArTnbN4c/N5ZwzpcAD82nUfXiuBX9bN1TenfG8Hh1ayKmHcYmS5tlM52ZZue+UY5ZoSBCRKQeNEIQsQQ46Id/fd2fdX9sZnYhcAbw0TK0reaOXdHH1j0jJX+7j4KIatVEQOlTXz+2f4x0xtmwTPUQIiL1oOpBhJm93Mw8xu2WCjz2WcBngK+5+zcK7Hehmd1pZncODAyUuxlldeyKPkYnMzw5OF7S8dUdnRHUXZRaF7FtYBRAmQgRkTpRi8LKnwMnxtgvqmU4ACwyM8vJRkQZiP3EYGbPA64DbgLeUWhfd78SuBJg06ZNpXfgV8GxKw4VV65a1F308SNVzEREU2uXOkJj20BQ+3HMMgURIiL1oOpBRFjoWEw9wmagEziGZ9ZFRLUQD8x1AjM7BfghcA9wvrvPfz3qOnFcGERs2TPCS45fXvTx0YJYvR3VqYmA0oOIRwdGWNbXycKeyo8kERGRuTVCTcQPgDTwppztFwD3u/v2Qgeb2XHAj4BHgVe7e+kzM9WhpX2dLO5pL3mExlAqTX9XGy0tlV+HYr41EcHIDNVDiIjUi7qfJ8Ld95jZp4CLzWwYuItg0qizgXOz9zWzG4G17n5s+PcKggCiA/h7YGPOok13u/tE5a+isuYzQuPg2CQLu6vzzX4+NRHuztY9I7zq2UeWu1kiIlKiug8iQh8ERoD3AkcADwOvc/fv5uzXyjOvaSOwNvw9d1+A9cCOsra0Bo5Z0cf3738Sdy96ZcvBVLqKQUTp3Rn7RycZTKU1R4SISB1piCDC3TPAZeGt0H5n5fx9C5D49aKPW9HHf4yl2Tc6ybK+zqKOHUylWVSlGoO+jjZaLHjMYs2MzFB3hohI3WiEmgiZQ/YIjWJVMxPR0mIs7G7n4FgpQUQ4MkOZCBGRuqEgIgGiNTS2lBRETFUtiABY1NPBwVIyEXtG6GxrKWkYq4iIVIaCiAQ4YkEXfZ1tbCsyiHB3BlOTLKhiEBFkIiaLPu7RvaOsX9ZLaxVGkYiISDwKIhLAzDhmeW/R3RmpdIZ0xlnU3VGhlh1uUU97iTURI5qpUkSkziiISIgNy/vYvne0qGOiD/Nqdmcs7C4+iBhPZ3hs/5jqIURE6oyCiIRYv6yXxw+mSBWxQmYtgohFJRRW7tw3xrRrZIaISL1REJEQG8IP2GKyEYNjNchE9HQwNJ4mMx1/SRKNzBARqU8KIhJiQ7goVTFBRDRKolrzRECQiXAvbtbKqGB0gzIRIiJ1RUFEQqxb1gMEi1TFVZPujDBgKaZLY9vACKsWdtFThUXCREQkPgURCdHT0caqhV08WkQmYigMIqo5xHMmiCiiuPLRvaMamSEiUocURCTIhuV9RQURg6k0ZtDfWb1v+AvD4aRx54pwd7btGVE9hIhIHVIQkSDrl/Xy6MAI7vGKFqMpr6uxDHgk6jqJO8zz6aEJRiczGpkhIlKHFEQkyIblvQyPT7F3JN63/INj1Vs3IxJ1Z8QNIjQyQ0SkfimISJANy4sboXEwlWZRlYOIKGiJW1g5E0SoJkJEpO4oiEiQDcuClH/cERr7RiZYWuTS4fPV3tpCX2db/CBizwh9nW2s6K9uO0VEZG4KIhJk1aJuOtpaYhdX7h2ZYFlf9dbNiBSzCNe2gVGOWd6LmRbeEhGpNwoiEqS1xVi3tIdHB+YOItydfSOTVc9EACzr62DvaLwg4tGBkZluGhERqS8KIhJmw7I+Ht07d3fGYCrN1LSzrCZBRCcDwxNz7jc6McUTg+MamSEiUqcURCTMhuW97No3RjozXXC/aARHLbozlvd3sndk7iAiKhDVyAwRkfqkICJh1i/rZWra2X0gVXC/feGH+NLe2mQi9o1MzLkIl0ZmiIjUNwURCRPVD8w1QmMmE9Ffm0zEtMOBOYort+0ZocVg7dKeKrVMRESKoSAiYaL6gbmKK/eN1jYTAcxZF7FtYJQ1S3robGutRrNERKRICiISZlFPB4t72ucc5rl3eAIzWNJbm0wEMGddxLYBrZkhIlLPFEQk0IblfXN3Z4xOsqSng9YqrpsRiYKIQpmIzLSzXat3iojUNQURCbRhWW+sTEQthnfCoREhhTIRTxxMMTE1PTMLp4iI1B8FEQm0fnkvA8MTDI/PPrX0vtFJltZgeCdAX2cbXe0tBTMRW8NMiiaaEhGpXwoiEmjDsmiExuzZiGDK69pkIsyMZX2dBVcb3fL0MADHr1QQISJSrxREJNBx4Qfvlj356yLcnacGx1m5oHaLWi3vLzxr5cNPjbCiv5NFPbXJloiIyNwURCTQ2iU9dLS18PBTQ3nvHxiZYGJqmtVLajf/wvK+TvYMj896/5Y9wxy/sr+KLRIRkWI1RBBhZi1mdrGZ7TCzcTO718zOL+E8G8xszMzczI6tRFvrQVtrC8et6OPhp/NnIqLZLI9e3F3NZj3DUYu72X0ghfvhs1ZOTztbnh5RECEiUucaIogAPgx8CPgc8ErgDuDbZvY7RZ7n88BgeZtWn05Y2T9rJuKx/WMAHL24dpmItUt6GJvM5K2L2H0gRSqdUT2EiEidq/sgwsxWABcBl7v7J9z9Znd/F3AzcHkR53kj8BzgY5VpaX054Yh+nh6a4GCeqaXrIROxJpzKelcY0GR7JCqqPEKZCBGRelb3QQRwDtABXJWz/SrgFDNbP9cJzGwx8CmCYORg2VtYh04IP4Affmr4sPt2H0ixtLeDno62ajdrxpqwHuOxPEHEw2EQcZwmmhIRqWuNEEScBEwAW3O2bw5/boxxjn8CHnL3r5ezYfUsCiIeyhtEjNU0CwGHulJ27js8iPj17kHWLu2hv6u92s0SEZEi1O6raHxLgIN+eAXe/qz7Z2VmLwbeQtCVEYuZXQhcCLBmzZr4La0jRyzoYllfB/fuPjzxsvtAio1HLqhBqw7pam/liAVdebsz7tt9kN9YV/BpFRGROlD1TISZvTwcHTHX7ZYyPFYHcAXwz+7+QNzj3P1Kd9/k7puWL18+32bUhJlx2upF3PPYM4OI6Wnn8QMpjl5S20wEBHURu/Y/c0KsgeEJnhgc59SjF9aoVSIiElctMhE/B06MsV/0FfUAsMjMLCcbEX1V3c/s3gcsBj5jZovCbdGQhH4z63f3w/P9CfGcNYv58YN7GBxLs7An6Bp4/GCKycw0q2s4MiOyZkkPP90y8Ixt94WZk2cfvSjfISIiUkeqHkS4+xjwUBGHbAY6gWN4Zl1EVAtRKMOwETgCeDzPfXcB9wKnFdGWhnLa6uCD+J7dBznz+CCjcu/Mh3Ttv+mvX9bL1b/azdB4mgVh/cO9uwdpMTj5qNp2t4iIyNwaobDyB0AaeFPO9guA+919e4FjLwdemnOLhnheALyjvE2tL88+eiFmcM+uQ10a9+w6SGdbC886ovYf0lGQc9fOAzPbfrVzP8ev7K/pyBEREYmn7oMId99DMDzzYjN7v5mdZWZfAM4GLs7e18xuNLOtWcc+5O63ZN84lAX5hbvfWaXLqIn+rnaedcQCfrb1UJfBPY8d5OSjFtLRVvun/jlrFtHWYvxye9AjNTye5pfb989kTUREpL7V/pMkng8ClwHvBX4IvBB4nbt/N2e/VhpjxEnVvOKkI7hz5wGeHhonnZnm148PzmQAaq2no42Tj1rI/+4IgoifPLKXdMZ5+caVNW6ZiIjE0RBBhLtn3P0yd1/r7p3u/mx3vzrPfme5+7o5zvUVdzd3z513IpFe9ewjcIcf3P8U9z8+yMTUNM9ZUx9BBMDp65dw72ODjKcz/PjBp1nc085z1yyudbNERCSGhggipHTHrujnuBV9/Mcvd/GJGx6mv7ONFxyzrNbNmvGbG5YymZnmI9c/yPfvf5KXnbiS1hardbNERCQGBRFN4KJzTmDrnhFu27qPv/rt41nS21HrJs048/jlnHPSSr52+04WdLXz1+ecUOsmiYhITKofaALnnHQEX3778/jJIwNccMbaWjfnGVpajE++7jSWXf8gr3/ealYs6Kp1k0REJCY7fDZpybZp0ya/885ED+IQERGZYWa/cvdNcfZVd4aIiIiUREGEiIiIlERBhIiIiJREQYSIiIiUREGEiIiIlERBhIiIiJREQYSIiIiUREGEiIiIlERBhIiIiJREQYSIiIiUREGEiIiIlERBhIiIiJREQYSIiIiURKt4zsHMBoCdZTzlMmBvGc9XS7qW+pOU64DkXEtSrgN0LfWoEtex1t2Xx9lRQUSVmdmdcZdYrXe6lvqTlOuA5FxLUq4DdC31qNbXoe4MERERKYmCCBERESmJgojqu7LWDSgjXUv9Scp1QHKuJSnXAbqWelTT61BNhIiIiJREmQgREREpiYKIKjCz1WZ2tZkNmtmQmV1jZmtq3a5CzOwsM/M8t4M5+y02s381s71mNmpmPzazU2rY7qPN7LNmdruZjYVtXpdnvy4z+7iZPWlmqXD/l+TZr8XMLjazHWY2bmb3mtn5dXYt+Z4nN7PT6uFazOy1ZvZfZrYz/Ld+2Mw+amb9OfvFei3Ffe5qcR1mtq7A87GoHq4jfOxzzOwmM3vKzCbMbLeZfcvMNubsF+u9q5bvA3GuJe77Wa2vJU9bfhC287JS2liV15i761bBG9ADbAHuB14DnAf8GtgG9Na6fQXafRbgwJ8DZ2TdNmXtY8DPgN3AG4BXALcSjFk+uobtfhq4HvhheA3r8uz3DeAg8E7gZcA1QAo4LWe/fwQmgIuAlwJXANPA79TRtTjw5Zzn6Qygpx6uBbgD+BbwJuBM4H3hv/0dQEuxr6W4z12NrmNd+Hx8JM/z0VoP1xE+9huAjwOvDa/lzcBmYIhgjgCI+d5V6/eBmNdyFnO8n9XDteS5rifDdl9WShur8Rqr2j9Is96A9wIZ4NisbeuBKeD9tW5fgXZH/+leXmCf88J9Xpq1bSGwH/hMjdrdkvX7O8jzwQucGm5/e9a2NuBh4LqsbSsIPnT/Ief4G4H76uFawvue8SYzy7lqdi3A8jzb3hK2++xiXktxn7saXse68O93zHGuml1HgTadELbpr8K/Y7131en7QO61zPl+Vk/XAiwGniIIEnKDiLr6v6LujMo7F7jD3bdGG9x9O3AbwYuhkZ0LPOHuN0cb3H0Q+B9qdG3uPh1jt3OBNPDNrOOmgP8EzjGzznDzOUAHcFXO8VcBp5jZ+vm3eHYxryWuml2Luw/k2fy/4c+jwp9xX0txn7uyi3kdcdXsOgrYF/6cCn/Gfe+qu/cBDr+WuOrlWj4G3O/u/5Hnvrr6v6IgovJOIkgH5toMbMyzvd58w8wyZrbPzP49pz+00LWtMbO+6jSxaCcB2919LGf7ZoIP2mOz9psAtubZD+rr+fvTsD94LOwffnHO/fV2LWeGPx8Mf8Z9LcV97qol9zoiHzWzqbCW4Lo8/dV1cR1m1mpmHWZ2HEH31lNA9MEV972rLt4H5riWSKH3M6iDazGzFxFkuP5sll3q6v+KgojKWwIcyLN9P0HKql4NAp8kSKOfDXwYeDlwu5mtCPcpdG1Qv9c3V7uXZP086GEesMB+tXYV8G6C5+dCYClwk5mdlbVP3VyLmR0FXAr82N3vzHr8OK+luM9dxc1yHRMEH2DvIqg7uQg4Bfi5mZ2YdXi9XMcvCNr8CPBsgm6ZPVltiPPeVS/vA4WuJc77GdT4Wsysg+D18wl3f3iW3erq/0pbOU4iyePudwN3Z2261cx+AvwS+Avgb2vSMDmMu78568+fmtm1BN9ULgNeVJtW5Rd+S7qWIM389ho3p2SzXYe7Pwn8SdauPzWzHxB8+/sgcEE12xnDm4EFwAaCgOdHZvYid99R01aVZtZraaD3s78GugmKoBuCMhGVd4D80etsUWLdcve7CKL854WbCl1bdH89mqvd+7P2W2RmNsd+dcXdh4Hvceh5gjq4FjPrJui33QCc4+67c9oX57UU97mrmDmu4zDu/hhBNX3u81HT6wBw9wfd/Rdh3/vLgD7gA+Hdcd+76uJ9YI5rybd/7vsZ1PBawq6VDwJ/B3Sa2aKsYcHR361FtLEqrzEFEZW3maBvKtdG4IEqt6VcopR4oWvb5e4j1WtSUTYD682sJ2f7RmCSQ3UDm4FO4Jg8+0H9P3/ZXRc1vRYzaweuBjYRDCn9dc4ucV9LcZ+7iohxHYXkPh81u4583P1g+LhRX3nc9666ex/Icy0Fd8/6vZbXsgHoIuiePJB1gyCzcoCga6yu/q8oiKi864AzzGxDtMGCCYNeGN7XMMxsE8HQqV+Gm64DjjKzM7P2WQD8LvV9bf8DtAN/EG0wszbg9cAN7j4Rbv4BQXXzm3KOv4Cgcnp7FdpatPA5eDWHnieo4bWYWQvBePWzgde4+x15dov7Wor73JVdzOvId9wagm6l7OejZtcxGzNbCTyLYB4IiP/eVXfvA3muJd8+ue9nUNtruYegjib3BkFg8VKCD/76+r9SrXGvzXoDesMn/tcEw2/OBe4FHgX6at2+Au3+BkGf+u8TvGn+FcFkJruAZeE+LcDPgceAPyQYRngLQZpsdQ3b/trw9gWCbxl/Gv59ZtY+/0kQ2b+DIPV5NTAOPDfnXJeH299PMNb8CwQTNL26Hq6F4BvKl4A3hu17a/hamwReXA/XktX2yzh8Aqaji30txX3uanQdnwT+GXgdwZv+nwA7CSb8OaEeriN87O8QpM3PC9v5LuChsJ3Hh/vEeu8q5rmr4bXM+X5WD9cyy/XlzhNRV/9Xqv4P0ow3YA3wXwQzqA0D/02eSYPq6QZcDNxHUNWcDl+wVwJH5uy3BPi/4Qt4jGDyolNr3Haf5XZL1j7dwKcIhoGNE1R2n5XnXK0ERVc7CSq/7wNeWy/XQvDt47bwDTFNMD7+OuD0erkWYEeB6/hQsa+luM9dLa4D+COCuSMOhM/HU8C/kxNA1PI6wsf+G+BXBB+0YwQTEF2R+74U972rlu8Dca6FmO9ntb6WWa7vGUFEMW2sxmtMq3iKiIhISVQTISIiIiVRECEiIiIlURAhIiIiJVEQISIiIiVRECEiIiIlURAhIiIiJVEQIdKEzMxj3HaE+34l+r1emNlnzOy7VXy815jZ03W8vL1ITWieCJEmZGZn5Gz6DsFshB/K2jbh7neb2THAAg9WQqy5sD0PAi/wQ0twV/oxjWAVyGvd/e+r8ZgijUBBhIgQZhp+5u71tlT1Yczss8AZ7v68OXcu7+O+G/gwcJS7j1fzsUXqlbozRKSg3O4MM1sXdnf8iZl91MyeMrNhM7vKzHrM7Fgz+6GZjZjZVjN7a55znmpm15nZATNLmdltZvbiGG3pJFg07N9ztp8Vtuk1ZnaFme03s4Nm9mkzazWz55nZz8xs1Mw2m9k5Occ/z8x+ZGb7wvY8amafz3n4bwGLCNZfEBEURIhI6S4GVhEs+nUJweqAXyToGvke8HsE6xV82cxmli42s+cSLCC0BHgncD7Bmh8/NrPfmOMxzyD4IP/pLPd/GhgN2/JZ4L3htq8RrDXw+wTrDVxjZsvC9vQBPwQywNuAVwKXAm3ZJ3b3vQTdKK+Yo40iTaNt7l1ERPLa5u5RluGHYSbhzcCb3f0qADO7k2D1x9cCm8N9P06weuLZ7j4Z7vdD4H6C1RhfU+AxzyBYkOi+We6/yd3fH/7+IzN7FfAeghVNfxY+1pME9R+vAr5KsGT0YuCv3T37vF/Jc/67wzaICMpEiEjpvp/z90Phzx9GG9z9ALAHWA1gZt3AmcC3gWkzazOzNsCAHwMvmeMxVwFDUfARs02jUQCR087V4c8tBCtAXmFmF5jZamY3ELZBRFAQISKlO5Dz92SB7V3h70sIliT/O4IlmbNv7wEWm1mh96UugmXMi2nTwewNWQFIV/j3IPBS4Ang88AuM7vfzM7Pc/5U1rWIND11Z4hINR0EpoF/IahTOIy7Txc4fh9BTURZufs9wPlhVmQTQb3Ht8zsVHe/P2vXJWEbRAQFESJSRe4+amY/BU4F7pojYLLbGEoAAAEpSURBVMjnIaDDzI52990VaN8UcIeZ/R1BLceJBLUakfXAw+V+XJFGpSBCRKrt/cBPCIox/w14ElgGPBdodfcPFDj2J+HP04GyBBFm9mrgQuC/ge1AL/AXwDBwe9Z+Fj5u7tBPkaalmggRqSp3vwt4HkG3wGeAG4D/HziFQ0HCbMfuAH4J/G4Zm7SFoNbh7wgKM78MTAG/lZPteAHBKI7/LONjizQ0zVgpIg3FzN5GEHQc6e5jVXzcLwAnu/uck2KJNAtlIkSk0VxFMJLi3dV6QDM7gmBSrQ9W6zFFGoGCCBFpKGHx49uBqmUhgHXAX7l7we4WkWaj7gwREREpiTIRIiIiUhIFESIiIlISBREiIiJSEgURIiIiUhIFESIiIlISBREiIiJSkv8HlpgLNwGi9NEAAAAASUVORK5CYII=\n",
"text/plain": [
"