{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from matplotlib import rc\n", "\n", "import numpy as np\n", "from matplotlib import pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "plt.rcParams[\"text.usetex\"] = True\n", "plt.rcParams[\"font.family\"] = 'serif'" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def setBwStyles(ax):\n", " from cycler import cycler\n", "\n", " ax.set_prop_cycle( cycler('marker', ['.', 'x', 4, '+', '*', '1']) +\n", " #cycler('linestyle', ['-.', '--', ':', '-.', '--', ':']) +\n", " #cycler('markersize', [15, 15, 15, 15, 15, 15]) +\n", " cycler('color', ['k', 'k', 'k', 'k', 'k', 'k']) )" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/andreb/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:35: RuntimeWarning: invalid value encountered in true_divide\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yb13Xw8d8FuCc4RUkUKZGSrGXLJimPeNuUVxM7cWQpjptmNJbTN2naLMtJ2zdp2iaR7HSkTVrTGW/ipl60MxwndkzbqXcsidrU5BBJkRJJkOAEiXXfPx4AIikOcIAPQJ7v56OPADwYFyTxHNx77j1Xaa0RQgghACxmN0AIIUTkkKAghBAiSIKCEEKIIAkKQgghgsIaFJRSJTM5LoQQYm6FLSgopcqBx6Z7XAghxNwLW1DQWlcBndM9LoQQYu5JTkEIIUSQBAUhhBBBEhSEEEIExczliymlbFprR4j33Q5sB0hOTi5ds2ZNWNsmhBDzzd69ezu01jlTeYwKV+0jpdQWjNlF92utK/237dVal453fDxlZWV6z549YWmnEELMV/5zbtlUHhO2noL/RF856rbSiY4LIYQwl+QUhBBCBElQEEIIESRBQQghRJAEBSGEEEESFIQQQgRJUBBCCBEkQUEIIUSQBAUhhBBBEhSEEEIESVAQQggRJEFBCCFEkAQFIYQQQRIUhBBCBElQEEIIESRBQQghRJAEBSGEEEESFIQQQgRJUBBCCBEkQUEIIUSQBAUhhBBBEhSEEEIESVAQYTfg8qC1NrsZQogQSFAQYdXWO8hl33yZt07ZzW6KECIEEhREWJ0818eQx0eLw2l2U4QQIZCgIMKquWsAgCGP1+SWCCFCIUFBhFVTp9FDGPL4ANjX2EVNS88F93N5fHz6p7vZ3+SY0/YJIUaSoCDCKtBTcHmNoPDV5w6x66VjwWPff+0UWmsa7P1UHW1jd32naW0VQkhQEGHW3GX0FFweH1prTtsHGHAZQ0nP7j3Dwy8dx97v4rTdCB4DLi+Hmrv56GPvMuj20trtpNE+gMfr4/Z/e4MXD7ea9l6EWAgkKIiwagr0FDw+7P0unG5vcCipsdM4Nuj2ctreD8CA28Nrx9t4u9bOuZ5BvvrcIb749H6aupwcbe3h0JnuC15Daz3mkJQQYuokKIiwGfJ4OdczBBhBIRAEhtxGT6Gx0wgERlDwBwiXl4YO4/b+IS8nz/Vh73dR39EHQN+gB59v5JqHV462ccf33uCPdXa+/MwBGu0D/OUT++gddIf/TQoxz0hQEGHT4hgMXh7y+GgKBIULego+GgI9BZc3eLlrwEVLt5O+IQ/1HcZ9OwfcXPWdV3h6T1PwuWtajV7C4++epnJvM9979STPH2jhUPPIXsUv9jVz1H/fs92DvCf5CyEuIEFBhE0gCIDRU2ga1lMYdJ/vRTjd3mCAcA7rNRxt7UFr6B/yBHsKR1t7ONczxKm2PgbdXnw+TW27cWzv6S4AdjcYJ/uOfhcAvYNuBt1eHqw8yI/frAfg3145wSd/8t4FK61fO97GV545wN7TnXh9Oni8f8gzyz8dISJTjNkNEPNXIMkcZ7Xg8vpGTE9tHBYw+oY8nPHft61nCLv/ZH7EnycYcHmpbTN6D4EA0Nnv4tpdr/FXN68K3tbabfRMAkGluctIUN/2r29QUpiB26s5draXk+d6OdraS7/Ly1XffpVnPnMVyzKTAPjkT3YD8Mze5gvez7995FLuunRp8LrWmmerz3DL+kWkJcTS2u3EohSL0hJm/LMTwiwSFETYNHcNEGNRLM1IxOXxcdZ/0h50e2m0nw8KtW19ePx5gqNnzyeMj7R0X3A58MX+SEsP7b1GjyEQMEbb9eJx3qm1c8bhxN5v9EoOnelm87+8TnKcFYCzPYPc8b03+OGflXFFUdaE76dyb3MwKNR39POvVSf41f4WvtJzEduvK+L//LwareHzN68EYMOSdDRGL8mnNXsauuh2uvnUNSs41zNIdko8ACfbelmTl8aP3qynurGL796zkYRYo32N9gFsybGkJcRO8tOemNPlJcaqiLXK4ICYmAQFETZNXU6W2BJJiLUaOYWu8zmF08N6CifO9QKQGGuld/D8ME1t+/mTfc/gyOGb4/7gcfxsL073+Kul3zjZARh5i+H6Xecf0zvoYVvFu3z2xuIJ30//kIe69j6e3N1Exet1wdt/uc+YWhvw5WcOYlGwIjsZx4CbjOQ4znQ5iY+xcLpzgMtXZHLPf73Dtk3LuKzAxl89uZ/v3rORZ/Y0cexsL26Pj+/fV4JPa+743hskxFr42h1ruevSpVgtasI2jlbf0c//e6uex989zZbSfL5z9yU8sbuRU219/N/3r+P5g618+7dH+fs713OqvY8tpfkcbOrmyuIsqmrOUVKQwX+8dpL3FWezfkka2SnxpCXGEmiGx6ex97no6BvC6fZySX46bT1D/PpAC7euX8SvD7SybdMyhtxGruj61bk8/k4DsTEWbl2fFwyMR1q6ebnmHBlJcSTHx+Dx+nD7NHFWxV2XLg0GyeauAf5Y18nF+eksTk9gz+ku3lecRatjkMKsJGrb+8nPSOR/T7RTkJlEakIMHq/G7fVhtSgsyviS4vVp4qwW3D4fLo8Px4Cb/IxEDp/pwac1K3KSJwzETpeXMw4naYkx5KTEY+93kZEUh09rYiyKfpeXhBgLMcOCsNYaj09zrmeQvLQEzvUOkZUch9WiaHUMkpMaT6L/y8pYvD5NZ7+LtMQY4qwW7P0uspLjaO8bwpYYF/zb8GmN12f8mw4VzuqVSqkSrXX1OMe2AA6gRGu9a6Ln2bh0qf71TTeHo4liFriJpcuaS5clh35LGgMqlUGViMtnQVusaKZ2IhNCzI7PPXb7Xq112VQeE7aeglKqHNgJlI5xrARAa12llCqaKHiIyORUSbRal9MSswKHJQuUBbSPJN1Pku4l29uNva8LpXqw6Fi0Vmg0CuN/i1JoDZrzX0piLBY8Pl/wPsCIy1al8I7zJWb4/WZKKTXtUt+TtcOiFD498r0Nf0xCrBWtR9aKSoqLweszvu0Of6zFolCAAjTGN1HfsJ+pQhFrVcRYLfi0ZtDtRaGIi7EEnz/WakGhcHm9wbYF2+X/OcRaLXi8esL3pVAoRbB9sVYLXp/xrdmrwevzEWe1BIcJ42OseHw+vD4dfExCrBWLGvkFwuXxBf8mLBbl/7Zv3G5RxmWPz/h78vo0FovRZqtFBYcatf9nFODT59+bRQXub8GrNRZF8DbjZ8n5MUsVeKdgUcbfSeDbuFLg8w17Xkvg73v0z8l4rE8bjzHag/Etf/hrBl7Xfyflfw1F4Hfkf01tPJ8O3D3YzukLW1Dwn/DHm/O3DXjZf7kOKAfGDQqxixdT+PjPZrmFYjrO1ndzoKqJ2n3taJ8me1kKZRdns2S1jbwV6cTGWzl5rpelGYmUPvYRsm39FDq/wZunjGGc29bn8eKRs2Qmx7F+SVpweGf1ohRKCzN54r1GLl+RSU1LD31DHq5dlR28z1/cUMx//qH2gjalxsdwaYGNN091cOWKLP5Yb2esnnNZYQZ7/DOU8tIS2LnlEj7+4/cuuN/Jf7odi1K4PD4G3V4ykuM4ca6XW/7l9RH325ifzt+9fx0f+9F7DHq8fLF8NRaL4p9fPsHmtYt48chZrludw+sn2vni5tX88I06nv/La7j+4T8A8MM/K+PTP9vDvZcXcO/ly/jBa7U8snUjKfExvHGynb9+cj+pCTH84Ss3AsZJv7a9j5dr2qhu7KKho5/2viF8Pk18rJUltkSWZyVRWpjBlUVZrMpNQflPKlprfl9zjsuXZ5KRHMfvj5wlxqq4ac0iANxeHx6vpqa1mzV5aexrdHD1yizaeoeCifMhjxeny4vL46Otd4jEOCuL0xNQKBLjrLi9Pk619bF2cdoFP9NBtzc4BDSaz6fpc3nGHK7RWuP2auJiJBcyHR978utTfoxZOQUbMDxgTJzhE6azn+nj3V/W0nDITnxSDBtvXsb6a5ZgW5Q04n517X1s/pfX+dodq7Em1rM6/Ra0+/wHetWiFF48YsweWr0oNXjCL8hMJsk/nroiK5mGjn76hjxsWJrOGyc7yE2NZ3G6cXIqykmmbli+oSg3hfyMRFbmpHBJfjoDbi//uu1SvvbcId6ps5OeGEu30829lxfw409u4rJvvszqvFSuX53Ds39xFTuePcSptr7g8wWSsYlx1uAYb3FOyoj3+YP7SrhpTS4JsVY+efVyUhNi+YsbinF7fdy8NpfMpDiuKMrk7pJ8/vdEO3duXMJf3rQSpRRPP3AVLo+Pa1Zl8/ifX86GJelkJMfxXx8736m+dlUOr37phhH5EqUUK3NTWZmbOuXfn1KKW9fnBa/fMuxy4D3HWqG0MBOAa1ZlA4yYSRUfYyU+xvh55I4xwyrWahkzIADjBgQAi0WNO36vlCIuRoYf55IkmsWE3ENe3nu+jgOvNhOXYOXKDxZx8Q35xCWM/acTOLm+Vr8PZXVRmlvKwS7jJBtjURRmJQfvuyI7mVirwu3VLM9KCp44CrOTSKo3Lm9Ykg7A8uxkcvxJybLCDOra+4mPsTDk8VGck8xX71iL0+UlOyXeSCLGWPjAxiW8U2fnQ5ct5bXjbVy+IpO0hFhuXb+Iq1caJ73Swkx++/lrqXi9lkd+f4KLFo19wrVaFJ+/aSU5qfG8/5IlZCTHBY89eNua4OVYq4U1ecaJ8ZNXrwDgzo1LAILf2i9fkRm8/7Wrcsb92acnxZLOzGYdCTFVZgUFBxD4ZNiAC7blUkptB7YDFBQUzF3LRNDZ+m5e/tERejoGWXftEq76YDEJyROfpAJrE072HIA0uGn5VRyvawNgaUZicCooQKE/ELi9HgqzkoIzjJZnJZMUZ/xpblhqnGCLspPZvG4Rv/rs1Rxp6eHpPc1ckp/O7oYuinNSSEs4P20zMAsjL90IIuVrF/GNO9cHX/cH941Mc8XFWPjcTav4yOUFJE7wjfaLt1w0yU9MiOg3pwN1Simb/+JTQJH/chFQNfq+WusKrXWZ1rosJ2f8b1Ni9mmtqX7pNL94uBrtgw99qYQb71szaUCA80Gh33Ic31Auq7OXEOcfjlmWkUR87Pk/uYLMYb2DrPPDR4VZSSTHW8lOiWdxeiLpibFcusxGjNXCxmU2inKSiY+x8L5i49t+cU4yY7l2VQ7/sm0j7ysObXQyOyWe5HjpPIuFLZyzj7YAZUqpLVrrSv/NrwClWutqpVSZf4aSQ2YeRQ6P28trjx/jxHvnKC7J4cY/XUN8UuhDGMZaBC/WxAYSXZuwWFQwSbgsMyk4Jm21KP8aBuNYYVYSmclxlK9dxMrcFNbkpbEoLYG4GAtvP3TTiG/wVxZlceDrt3CwuZv//EMtG5amj9mWWKuFD12WP82fhBALUzhnH1UClaNuKx12uSJcry2mZ2jAzQvfP0hrbTdX3LmC0tuXB8fBQ9Xc5cSS0IKyDrEobh3AsKCQSLz/8hJbArFWC4mxVmPVsy2RGKuFH37cmFL9Dx/cEHzOsb69J8RauXxFJvu/vjk41CSEmDn5NAkAnH0ufv1v++ls6eeWT69nVdmiKT+H1prmzgGsycZq31VplwLng0LBsJ5CYaYx5JMQa2VpRuKIlZ9TIQFBiNklnyjBYL+bX/7zPrrbndz+mYtZfnH2tJ6nx+mhd8hDYk493qEcVhUuBiB+WE4hMFwUKEB30aJUptgZEUKEkQSFBc7t8vLC9w/iaBvg/Z/byLI1mZM/aBzBfEJSA+6ejeRnGCf+BH8CuSAzKVjbqDDLOPbwPRtn9gaEELNKgsIC5vX6eOmxw5yt7+a2+zfMKCCAUazMktCKsg7iHSgiPyMRgC0l+RRkJpGRHEdqQgz3lOaPWEglhIgcEhQWsNefOMHpQ3au/+hFFJfkzvj5mjqdWJOMfIJ3YAXL/D2F3LQE3n+JsYArxmqR3oEQEUyCwgJ15I0z1LzZQslthWy4bunkDwhBc9cACan1LEpYxvsuWkV2StzkDxJCRBSpMrUAna3r5vUnT1CwLpMr7iya/AEhauzqRyXWc23BFTz6sbIpT2cVQphPgsICM9Dj4sWKw6RkxLP5z9djmeKmLRNp6DmBVoNsWrRp1p5TCDG3JCgsIFprXvvvYwz2ubntgYtDKlsxlec+564BoCxvSnt6CCEiiASFBeTo2600HOzgyg8WkbNs6uWXJ9LZ78IXX4stdjG5STNPWgshzCFBYYHobnfy5tMnWXqRjY03LZv15z/d2U9M4mlWp8vMIiGimQSFBUD7NK/8tAal4OaPr0PNYh4hYF/rMVTMAGWLSmb9uYUQc0eCwgJQ81YLrae6uWbralIzL9wxazZUtxmFbm9YfkVYnl8IMTckKMxzzl4X7/yiliWrbKy5KnyriGt7DoEnjTVZK8L2GkKI8JPFa/Pc28+dwj3o5fp7LwrLuoHH3z1NYWYibe6jJOlVsjZBiCgnQWEeaznp4Ng7Zym5tYDMJWPvTjZTj7x0nNJiH27VxbK4tWF5DSHE3Jny8JFSavnsN0PMNu3TvPH0CVIy4ym7IzxDOj2Dbrqdbs4OHQVgtU1mHgkR7UIKCkqp+5VS/6mUutt//e7wNkvM1In3ztLR1MdVHywmNn78zehnornT2I/5rOso2pvAhuzVYXkdIcTcCbWnUAc8BCjgM8DsFcwRs87j9vLur+vIKUid1g5qoWruGgBg0HoK78BylmWGZ4hKCDF3xs0pKKV+D9QCLwOZwG6t9bPAs3PUNjFNh147Q1/nEDf/2dqwrEkIaOpyoqy9WOPbGXKUBUtlCyGi10SJ5h0YPYRyYCVQqZTqAnYDVVrr/XPQPjFFg/1u9r7YQMH6LPJnuGnOZJo6B7AmNQDgcS5nqX9THSFE9Bo3KGit9/kvjugdKKUuA4oBCQoRaP/LjQw5Pbzv7uKwv1ZzlxNrUgPaF0uGdQVJcTKZTYhoN+anWCm1AkgfqzfgDxb7LnyUMNtgv5uDf2hmZUkuWUtTwv56zV0DWBMb8DqXkZ+RFvbXE0KE35iJZq11PaD8s44+rZS6dI7bJabhwKtNuAe9lN2xPOyvpbWm2dGFNaHFv/WmDB0JMR9MNny0D4yeg1LqfkADdVrrV+eofSJEQ04PB19tpujSnDnpJTgG3DittSQpjXdgOfmSZBZiXghpENjfc3gMLggQnVrr58LYPhGiQ6814XJ65qSXANDUZSSZFVa8zkKWZUpPQYj5YMqZwVEBIl0pdZP0HMzlGvSw/5Umll+STU7B7G6eMx4jyVxPYcoqXLlZXLEivDOdhBBzY1plLpRSaQBa624JCOY79s5Zhvo9lN5WOGevWd/hwJrQxFVLN/HSF65jZe7cBCMhRHiFWubiv5RSTymlPg3YgK3hbZYIlc+nOfBqE4tWpJFXlD5nr3vYfhhl8XLlEtmPWYj5JNScwmcAlFI3A5sx8gkiAjQc7KCn3cmVd81t5ZG63sNggZJc2WlNiPkkpKDgn5KaqbV+BXhFKXVTeJslQnXglSZSMuMpvixnTl+3w3OUpLil2BJsc/q6QojwCjWnsAkoDgwjAfL1MAK0ne6h5aSDjTctw2Kdu030PF4PQ9Z6FifI/glCzDehzj6qAmxa68fC2RgxNQdebSI23sraq5fM6evubqlBWQe5yHbJnL6uECL8prJOATBmH2mtG8LWIhESZ5+LU3vbWH/NUuIT57bm0JtNewEoXSQL3YWYb0LNKdyPMWT0MlCtlLpbFq2Z69jbZ/F5NOuvnbtewouHW/nfEx00Wffh86RyaV74i+4JIabJNTCth4X6FbMOeBqjjPZngI7JHqCU2gI4gBKt9a4xjj/of95MrXVFyC0WaK058uYZFhenz0lJi4DP/Hc1AEs21OAdKCQ/U0pbCBExhnqh6Y/Q8CY0vAUt1dN6mrBssqOUKgHQWlcppYqUUiVa6+phx8v9xyuVUjuVUkVa67ppvYMF6MzxLrrbnGyao5IWw6mYHnq950jwXimlsoUw02A3NL5rBIHTb0HLftBesMTAksvgqs8B35zy04Zrk51tGMGEYc8xPGxt9j8PGIGnHJDeQogOv95CfHIMxSW5c/7a1sTTAOTFr5nz1xZiQRvsMU7+9W/A6Tfh7CHQPrDGwdJSuOYLsPxqWHYFxAW2xp3FoDDDTXZsQOew61mjjtsxeh+B+44+LsYx0OOifn87F9+YT0ycdc5ed8jjBYygoH0xFKWtnrPXFmJB8rqheQ/U/QHqXjMuay9Y42HZ5XDdg0YQyN8EsbNXkHI6BfFmY5OdSuAB/+UsjCAhQnDsnVZ8vrlNMAO0OAYBsCadxutcRuFi2VRHiFmlNbQfPx8EGt4EVx+gjOGga/4aim6A/MshNiFszQjXoLCDkT2BESd9rXWdv5ZSYBHcBfkEpdR2YDtAQUFBmJoZXbTWHHv3LHlF6WTkJU/+gFnU3DUAyoUl4Qwu+3Usk/0ThJg5p8MIACdfhtpXobfVuD2zCC7ZagSB5ddC0txVIQ5XUHgKCFRKK8JY/IZSyqa1dviDQZnWukIp9YDWunL0E/hnJFUAlJWVSa0loL2xl67Wfq7/6EVz/tqNnQNYE5pRyofXWUi+7LQmxNRpDe3H4MRLRiBofMcYEkqwQfGNUHSjEQgy5q7i8WjTCgpKqTStdc94x7XW1UqpMv8sI8ewmUevAKX+40X+aauPTqcNC9Hxd89ijbGwsnTuE8xNnU6sSUaS2essYJlMRxUiNK5+Izl80h8IupuM2xddbAwJrboFlpaBNTJm8000JfXL4x3CmC1060RPPNbaA6116bDLF/QOxPi8Xh8ndp9j+SXZJCTHzvnrN3UNEJ/SiHcoB+VLZoktfGOaQkS93nNw/Ldw7AWofx28QxCbbPQGrvsKrNoMaXObFwzVRKEpG3gS6Aa24B8CEuZoPGxnsM/NmivzTHn9ps4+rCmNuLrXsig1gfiYuZv5JERUsNfCsd8YgaDpPUBDxgrY9OdGb6DwfRATb3YrJzXRlNSHApeVUnuHTVFFKZUR7oaJkY6/e5bE1FiWrTdn28um3ka8qf1Yh4pkP2YhwMgPtB44HwjaaozbF2+EG78Ga94PuWtBKXPbOUWhDmKVKuON1WEkjksA2YZzjgz2u6k/1MHF1+VjncMS2QH9Qx761EkSgHTLSlZkz+3MJyEihtZwphqOPAc1vzLyA8oChVfDbd+BNX8CtuieLRlqldSH/UXx7gH2aq0fCW+zxHB1+9vxeTSrLl9kyus3dQ1gTTxNkjWdH993B1kpkd8FFmLWaA1nD8Lh5+DIL8BxGiyxsPJmuOEhWH07JM+f9bdTqZJaDOwFnlFK3aS1lp7CHDm1t4207ARyC1NNef2mTifWxNOsybiYtYvnbh9oIUx1rsboERx+DjprjZpCRTfA9TuMHkHi/Nx1MNTho1qt9WNKqcu01t0qysbIopmzz0XzsS4u21yAWT/34+2tWOI72JR3mSmvL8Sc6T4Dh56GA09B+1FjaGj5tXD152HtnXO6iMwsU80p2JRSGihFcgpzom5fO9qnWVk292sTAg51HADgqqWlk9xTiCg01AdHn4cDTxjTR9FGUbk7HoF1d0GKeZ89M0wlp/AdjCTze1rrh8PbLBFwck8b6bmJZOfP3b4Jo9X3HQEdw4acDaa1QYhZ5fMaAeDAk3D01+AegIzlxtDQxm1GmYkFKtScwkvAA1rrBqXUh5VST2mtt4W5bQveQI+LlhNdlN6+3LShIwC7+zgpsYXEWyXBLKJcZx1U/8wYHuptgfh0o8bQxnuN3oEMjYc8fFQB7FJK1Wqtvyo5hblRW92G1phS1iJgyDPEkLWRwvgJF7ALEbk8Q8bwUPVPjd6BssDKzXDbt4yZQ2GsOBqNQg0KWmu9VSl1s39Htr2EsAObmJlTe9vIyEsic4l56wL+eOYgyuLhItslprVBiGlpPw57f2rkCpydxvqBm/4WLr0vYktMRIJQg0KxfxrqK8Ar/vyCCCNnn4vWUw7Th47eajZqGZYtlplHIgq4B6Hml7DnJ9D0rrGeYM0dUPoJWHEDWOZ+8We0CTnRPOqm/wpDW8QwDQftaA0rNmab2o5DHQfxuTJYv2ipqe0QYkLdZ2DPj4yewUAHZBbD5m/Cxo9CSo7ZrYsqE1VJfQoILFrbCXQFDgGXAavC3roFrP5AOykZ8eQUmLNgLdiO3hqjVLZsqiMijdbGnsV/fNSoPaR9cNEdcPn9xiIzyX1Oy0Q9hYe01j1KKQfGzKP6wAH/Ps0iTNwuL001nay9eompQ0fn+s/R5+0g3ncNyfGRUetdCFwDxgKzP1ZA2xFjg5qrPgubPm3q5jTzxURVUusD//v3Vnhk2LGZ7tEsJtBU04nH7WPFpWYPHR0CIC9utantEAKAvjajV7D7hzDoMDapufPfYcMWiJOe7GwJ9evfiD2UlVKXaq33h6E9Aqjf3058UgxLVplbW+Vg+0HQMRSlz/32n0IEtZ+Ad/7dWFvgdRl1h676LBRcJUNEYRBqUPiIUmonUI3kFMLK5/VRf6iDwouzTCmTPdz+9v34BpdQuCjN1HaIBUhrOP02vP3vcOJ3EJMAl90HV34Wslea3bp5LdSg8JTWemvgilLq5jC1Z8FrPdXNUL+Hoo3mzphw+9wc6ajBM7BJNtURc8fnM4LAG9+FM3shMROuf8jIF8gsojkR6pTU4EI1pVQaIPWTw6T+UAeWGMWydeZVY3R5fBxsP4rLNyQzj8Tc8HmNTWve+C6cOwy2QviT7xpTSiVfMKemsp/CA4AdY/hoL/BcGNu1YDUetrN0dQZxCebN9vnnl0/wm/pfQTJGUMiUD6UIE68HDj8LbzwCHScgaxV86FEjeWyVGW9mCPWn3qm1LlNK3ay1fkWmpIZHd7uTrrMDrL/W3IViR1q66fKdJMWSSZ83nSU2qQ0jZpnHBQefhDf+GbrqIXc9bPmJUaraYjW7dQtayKHYPy11n1Lq04ANkGmps6zxiB2Awg3mbu3X1DkAtkaSdRF5aYnEx8iHVMwSnxcOPQOvfcvY1nLxRtj2c2PRmZSgiAgh5xQC01D9Q0l9YW7XgnT6sO2/ZvIAACAASURBVJ30nERsi8wbrvH6NGd620nItTPYf60MHYnZobVRqfS1f4L2Y5B3CXz0aVh1i0wrjTATlbn4PedLW/hvUprzU1Kl/tEs8ri8NB/vYv015lZvPNcziC+uEYD2jjyuukiCgpgBraH2VXj1H6BlH2Svhnt+amxtKT2DiDRRT2GnvyrqBSSnMPvOnHDgdfsiYujImtiI1hbcA0tkOqqYvsY/wivfhNNvQnoB3PUDuGSbJJAj3ERlLsYMCH5dExwT03D6sJ2YWAtLVpu7irmpy4k1sRHf4GLQcTIdVUydvRaqvm4MFyXnwu0PQ+nHIUZ27osGoU5J/fbwq8DNwKawtGgB0lpz+nAH+WsyiIk1N6l72t6LNaEJd3cpgOQUROgGOuH1h+G9x8AaBzf+jVGOIs68TaLE1IXaj1PAo/7LRcDu8DRnYepuc9LTMcil5QVmN4VjnadQVhdep9EWGT4Sk/K4YPdj8L+7YKgHLvtTIyCk5pndMjENoc4+emjY1Xql1E1has+C1HS0E4CC9eatYg5o6K2BOGPRWpzVwqJUWaMgxqE1HPsN/P7vjLUGRTfCLf8IeRvMbpmYgVCHj34PaMDhv2k38Gq4GrXQNB3tJC07gfQc84dqOjwniIlLRbszWZqdiMUi0wXFGDpOwe++YswsylkD91XCynKZXjoPhDp8NO5MJDEzPq+PM8e7WLlpkdlNYcjjxWVtoDB+Nd1KkZ8hQ0diFFe/kTd4+z8gNhFu+45RrM4aa3bLxCwJNSisGH7FXxTvqxjVU2VfhRk419CLa9DLsjXmDx0db2vDEt/GyvTbsKfGszI3xewmiUihNRz5Bfz+b6HnjFGorvwbkGr+lxkxu0INChn+PZt3aK0bMAJCBcYiNgkKM9B0tBMU5F+UYXZTeKvZqFxSsuhS/vr+K8lKkSmEAug4CS98Cer/F/Iuhi0/hoIrzW6VCJNQg8JerfXD/gRzA1Dk36azKHxNWxiaj3WSW5BKQor53e8DbQfQWnHNsssoypFewoLnccFb/2oMF8Umwh2PQNmnpGDdPBdqUCgNBAClVD1QPNm+CkqpLRiJ6RKt9a4JjhdprSum3PJ5wDXo4VxdD5feYv5UVIC63hq0axErssxdVS0iQNN78OvPQ/tRWH833L4TUnLNbpWYA6EWH6kA6rXWPwRsWusyjP0VHGPdWSlVAqC1rgIcgeujjtf5j9eNPr5QnDnhwOfTLFtj/tCRT/tod50k0VeEVWYcLVyDPfDCl+FHt8BQr1G07p6fSEBYQEJdp9ANvOK/vE8plaa1fniCh2wDXvZfrgPKMfZ3Hm4nsBmjp1A1pVbPE01HO4mJtbC42NzSFgANPQ146Gdx3GqzmyLMcvx38JsvQm8rXPEA3PS3EJ9qdqvEHAt1ncJlGCf64VVSb53gITagc9j1EeMRWutqpVSdUqoW2DGlFs8jzUc7WbLKhjXW/GqRh9oPAVCUus7klog55+yC3+2Ag08Zm91sexzyy8xulTBJqDmFcs6XuQDYMpMXVUrZMIaeHgUeU0pVa63rRt1nO7AdoKAgMsbcZ9NAj4uuswOsuWqx2U0BoPrcQbQ3nrXZxWY3Rcylk1Xw689Bfztc/xBc+yWIiTO7VcJEU5l9VB+4opR6eaI7Y5zwAxPvbRh7Ow+3Hfi21tqhlKrGCDIjktH+5HMFQFlZmQ6xnVGj5aSRjjG7KmrAgbaDeAfzKcySWUcLwlAvvPQ3UP1TyFkL9z4BS6Qivgg9KDyklNqJMSTUjTF8tGqC+z8FBPqfRUAVGD0ErfWI5LTWumohTm09c6KL2HgrOQXmj9kOeYdo6D2F13mNlMpeCOrfgF/9H3A0wdV/BTd8DWKlxpUwTKvMhVLq5onu7M8ZlCmlygGH1jqQZH4FKNVa71JKPaiUqgMyF+KU1JaTDhYXp2O1mp9PON55HK/24BvMp0BKZc9fHpexA9rb34PMIvjUS1BwhdmtEhEm1NlHo+se1YbwmAtO9Frr0mGXL1i7sFAM9LjobOln9eWRUSLgUIeRZE70rcCWZP4iOhEG9lqo/BS07ofST8Kt/yT7HIgxTbRH81PA/UAxxvTRwG5rgdlHEw0fiQkE8glLV5u/PgGMoBCrbSxOX4ySKpfzi9aw/3/gt18xitZtfRzW3Wl2q0QEm6in8JDWukcp5QAeGJVolozUDLScdBATZyGn0Px8AsDhjsMoVwGFWTJ0NK84HfDCF+Hws1B4DdxdAelLzW7VpNxuN83NzQwODprdlKiRkJBAfn4+sbEz7+lPtEdzfeB/pdRNSimNMZOoHKic8SsvYGdOdLF4pS0i8gndQ92c7jmNp+d2li2ToDBvNO+FZz5hVDS96e/gmi9ETc2i5uZmUlNTWb58ufRcQ6C1xm6309zczIoVKyZ/wCRCPiv5q6M+prV+BFiQZSlmg7PPyCcsWRUZU1GPdBwBwDWwlMJMGWOOelobeyT/2L+29FMvwXVfjpqAADA4OEhWVpbpAaGiooLi4mIqKyupqKhgx44dOBxjVvYJqqqqCt6nurqazZs3T+u1x3qsw+Ggunp0YQhQSpGVlTVrPauQ92j2V0iVjXZmKBLzCQqFV2YeRb+hPnj+88Zw0apb4UP/BUnm79MxHWYHBICysjLKy8vZssVYq+twOLj//vt55plnxrx/IBjYbMYXvqKiIoqKpjfbfqzH2mw26urqKCm58Dv5bP68Qu0pdGLUKfqWUurDGGsPxDScOeEgJtZCbgTlE7Li8sGXIEEhmrUdg8duNDbCufn/wr1PRm1AiBRVVVUjvq3bbLYxv6kHVFRUUF5ePu7jZ/LaASUlJVRWhnf0PtSgUIYx62grxkK0fWFr0TzXcsJBXnE61hjz8wlaaw51HCLdYlRGXWKTBUxR6eAzRkBwdsGf/cooVWEx/+8r2u3evfuCb+WdnZ3j3Btqa0fO1N+9ezednZ1UVVVRUWHM0K+rq6OioiIYXB544IFxX3v0Y8HoQbz88mQFJWYm1L+cWq31QxjlLroxCuOJKRrsc2M/0xcxQ0et/a3YB+1Y3IUstSUSEwGJbzEFXg+8+FV47tOw+FJ44A1YcZ3ZrZo36urqRgzhVFdXU1ZWFrxcVTVxcefq6mq2bt1KWVkZe/fuBc4PMQWCy3hBZqzHBkwUmGZDqGeBUn9OYYVS6lKgdLIHiAu11vrrHUVIkjmwaG2gZ4lMR402A53w8w/Duz+AK/4CPv48pEVGccX5wOFwkJk5cvjt0UcfZccOo6jzU089NWnS2WazYbPZRgwFlZSUsHfvXsrLyyccXhrrsQGj2zXbQl3R/LBS6jsYuYT3JtlLQYyjtbYbi1VFVD4h1hLL2Y4MSi6WoBA12o7CE/ca003v+j5c9qdmtyhs/v75I9S09Mzqc65bksbXP7B+wvvs2bNnxNBRXZ1RxDmQM9i2bVvwtrFUV1ezadMmwBgKeuCBB4I9j8A3/erqasrLy4O3OxyOYN5ivMfOhalMSX1Ia73VPyVVTMPZum5yClKJiYuM6YGHOg6xynYRjgFNoSSZo8OxF+CH5eAegE+8MK8Dglnq6urYuXMnDoeDyspKKisrqaqq4tFHH538wcOeIxBAsrKyqK6uDp7UN23aFBx62rNnT/Cbf2lp6aSPnQsT9hT8he92YvQQ0jF2T/u21vq5OWjbvOJ1+2hr6GXDDZGxotTj81Bjr+H6xX/CH0FmHkU6reH1R+C1fzRKXH/kfyBtidmtCrvJvtGHQyjJ3D179lBbW0t5eXlwCmrgpA4Ep7ECPPjggyMeG7g+fKYSEHzNiR5bV1c37RlNoZqo9tGHgU3APYHVzUqpdGC7UurT/v2aRYjam3rxenwsLk43uykA1HXX4fQ4sVmNTXUKJKcQuTxD8KvPwaGn4eKtcOf3IDbR7FYtaNu3b7/gtq1bt1JZWTnipD4Vk+UowBhymu7zh2qi4aNM/5BRsOaR1rrbn0+IjOkzUaS1thuAvKLICAqHOw4DYHUXAtJTiFjOLnj8biMg3PS3Rv0iCQgRKZAcDuXkPpaxFqUNN1d5hYmGj0bvljbcpKWzxUhna7tJy04gOT3e7KYARj4hNS6V7p50MpMHSU2QktkRp6sBfn6P8f/dP4RL7jG7RWISo4eEZtNc5RUmCgqb/JvgjHkMkLxCiLTWtNZ1U7A2claYHmo/xMXZF9N02sky6SVEnuY98MRHwOuGj/0Sll9tdovEAjFRULgHoyrqWEU1bga+GpYWzUM9HU6cPS7yIiSfMOAe4JTjFNcvu55n9g2wcVlkrJsQfkefh2c/Dal5cF8lZMvWJWLuTBgUtNZjlrOQ/RSmJpBPiJQk87HOY3i1l7WZ6znjcHLnxvk/iyVqvPeYsSHO0lKjflFKjtktEgvMuInm8QLCZMfEhc7WdhOXYCVjcWSUpg6sZM6JXYXXpyXJHAm0hj98B377ZbjodvjEbyQgmCxaSmfPtlBLZ4sZaK3tJq8oHYvF/HLAAAfaDpGXtJiefqMAnkxHNZnPBy/ugPcq4NL74APfA6t8NM0WTaWzZ5NUQAuzoQE3na39EZNPAHirqZq+nsU0dg4AMh3VVB4XPHe/ERCu+hzc+R8SECKElM4WYXG2vgc0ERMUOgc7GdDtOPvyabQPEGe1kJcmJbNN4eqHJ++Fw5VQ/g245R+l5HUEkdLZIizO1najFCxanmZ2U4Dzi9YG+4yeQn5mYsQMay0ogz3w3x+G2lfhA/9m7KEcAbuNifPGK51dV1dHdXU1u3btmrQo3nwunS2m6VxDD5lLU4hLiIwhgX3nDqG1YrBvCbXtfVIIzwxOBzz+IWjeDVt+DKWfMLtFYpSJSmdXV1dTUlJCeXn5hEM587p0tpge7dO0NfRQXJprdlOCqs8ewufKAR1PbXs/VxVlmd2khWWgEx7/IJyrga0/gzV/YnaLItvvHoKzh2b3OfMuhtu/M+FdJiudDca4/3h1iBZE6Wwxdd3tToYGPBEzdKS15qTjKD5nPgBen5bVzHOprx1++gFjP+WP/I8EhAgVSunsqqoqysvLxz1Rz9vS2WJmztUbi9YiJSicGzhHr6cT7+D5kgmFWZGxdmLe6z0LP7sLuk7DR5+E4pvMblF0mOQbfThMlsytqqpi586dFBUVsXnz5mBvYd6XzhYzd66hl5j4yFm0dsR+BACvc1nwNpmOOgd6z8L/+xPoaYX7noEV15rdIjED5eXlYxa+Wwils8UMnWvoIbcgNWJm9xzpOALagtVzvqyFBIUw62uHn95pBIQ/fVYCwjy2EEpnixnwun10NPey8cZlk995jtTYa7B6F7MqJ5Oa1h5yUuNJjJCtQeelQFLZ0Wj0EAqvMrtFIszmQ+ls6SmESUdzHz6PZtGKyMgnaK05Yj+Cq38J65cYbZLpqGEUmHbacRLu/R/pIYioIUEhTM41GEnm3AhJMp/pO4NjyIF7ID8YFGToKEyGeuHnW+DcEdj2uCSVRVSR4aMwOdfQQ1JaHCkZkbHTWjDJPLiUVYtS2bA0jStljcLsc/XDz7fCmWrY+lNYfavZLRJiSqSnECZtDb0sWpGGipDSBUfsR7CqGHxDeRRkJvGbv7yWrZsiJ98xL3hc8NTHoOld+PBjsPYDZrdIzICUzp5lSqktgAMo0VrvGnWsBNgLBAqHVGmtx64MFYUG+904zg1w0ZV5ZjclqKajhnRrIf2WWJbYZOP3WefzwS//AmpfMUpfb/iw2S0SMySls2eR/6SP1roKcASuD5OptVZa62KMbT93hqMdZmk73QNEzqI1n/ZRY68hzlvAsowkrBEyRXbe0BpefMiodnrz16H042a3SMwCKZ09u7Zh9BLA6A2MmKflDxYBRVrr8UsNRqFz9UZQiJQkc2NPI73uXob6l8qGOuHw+iPw3qPGfgjXfMHs1ohZIqWzZ5cNGP7TGzOjqZQqB6rGOhbN2hp6yMhLIj4xMvL4gSRzZ2euTEOdbbt/BK/9I2y8Fzb/g5S/nkfGK53tcDioqqpi165dE+YYpHT29GzWWk9v+V+E0lobK5kjpJcARlCIs8TT15tFgdQ6mj1HfgkvfAlW3wZ3/rtskDOPTFQ6e8+ePZSVlQXH+McjpbNHcgCBltsA+zj3GzdjopTaDmwHKCgomNXGhVO/Ywhnr5vcwlSzmxJ0pOMIBSkrsWNluQwfzY7T78Bz22HZFbDlJ2CNNbtF89LO93ZyrPPYrD7nmsw17Lh8x4T3max0tsPhwOFwjJv0ldLZF3oKCLyDIvxDREopW+AOSqkizucdLqC1rtBal2mty3JycsLUzNnX3tgLQM6yyAgKXp+Xo51HyYotBqBQgsLM2WuNbTRtBXDvExAnP9P5ZLLS2RUVFdhsNkpKSti1a9e4zyGls4fRWlcrpcr8OQOH1jqQsn8FKB121/AOjpmgrbEXpSA7QoJCfXc9To+TeF8hSkF+hpzAZqS/w9hGU1ngvqchKbxd+YVusm/04TBZMjcw9FNXVzeiYqmUzp6E1rpijNtKh12uA+bN2oSA9sZebHnJxMZHRqG5QJLZ3b+UvLQEEmIjo11Rye2EJ+6F3lb4+POQOXff3kTkGG/9gZTOFmNqb+wltyAyeglgBIWkmCQ6HOkydDQTPh/84jPGvsp3V8Cyy81ukYgwUjpbXKC/e4iBbhc5ERYU1mat5XD9IDeviZy9oqPOK38PNb80pp2uu8vs1ogIJaWzxQjtp/1J5ggJCm6fm+Odx1ltW0tH35AsXJuuff8Nb/0rlP05vO8vzW6NEGElQWEWtTf1goLsZSlmNwWAWkctQ94hcuJWATLzaFqa3oPffAGKboDbd8niNDHvSVCYRW2ne7HlJhGXEBmjckc6jCRzos9Y51GYKQvXpqT7DDx5H6Qt9a9FiIzfqxDhJEFhFrU39kbM0BEY+YTU2FT6+ozlITJ8NAVuJzz5UeP/e5+UqacLyK5du9ixYweVlZUopaioqGDXrl3j1imazEwK45lBgsIsGehx0e8YiqyVzPYjrMteR2OXk4ykWNITZdVtSLSGX30OWg/Ah38IuWvMbpGYQzabjZ07d7JlyxaKiorYvn07Dz744Ih1CJMZXsm0vLw8WE47GkhQmCXBlcwR0lNweV2c6DrBhqwNnLb3S82jqXjzX/xlsP8vXHSb2a0Rc6ysrGxKt4/mcDjCXsk0nCQozJJAUIiUlcwnuk7g8XlYn72e0/YBqY4aqhMvwSvfhA1bpAz2AjXeeoG6ujruueceKisr2bVr14hhoR07drBjh7H6es+ePezZs+eCfQ8ClVXnYve0mZCgMEvaG3tJz02MnHLZ/iTz6vS1tDicUggvFJ318Nz9kHexUfVUZhqJYbZs2RJcUfzggw+OWJMwPN9QXl5OUVHRiJXHgXpGW7Zs4amnnprTdk9VZJzB5oG2xh4WF6Wb3YygI/YjZMRn4HXb8Glk+Ggybic8/THj8rbHpchdhDj7rW8xdHR2q6TGr11D3te+Nq3Hju5FhJorCHe569kkPYVZ4Oxz0dc5RHaE5BNgWJK5cwCQNQqT+u2X4ewhuPsxyFhudmtElBlvX4XAUFE0JZqlpzALOhr7gMhJMjs9Tmodtdy47EZO2/1BQXIK46v+mbFq+bqvwOpbzW6NGGa63+hnKrC7Wl1dHbt27WLLli3U1dVRXV1NdXV1sMcQKIM9/P6BgnmVlZXBiqqBxwUuB/ZOiERKa212GyZVVlam9+zZY3YzxlX9+9O881wtf/7da0lINn/a5/62/Xzsdx/jezd+jzcOLOKJ9xqp+eatKBkjv1DLfvjRLVD4PvjTZ8EiVWTNdvToUdauXWt2M6LOWD83pdRerXVo06b8ZPhoFtib+0jJiI+IgADny2Wvz15PY2c/hVlJEhDGMtBp5BGSc+DDP5KAIAQSFGZFR3Mf2fmRUe8IjJlHOYk55Cblcto+QIEMHV0oUAq7pxW2/hSSs8xukRARQYLCDHncXrrODpAVQUGhxl7Duqx1+Hyaxs4BSTKP5d0fwMmX4NZvQf6UetdCzGsSFGaos6Uf7dNk50dGknnAPUB9Tz3rstZxtmeQIY+PQpmOOlLLPqj6Bqx5P1x+v9mtESKiSFCYoY5mY+ZRpAwfneg6gU/7WJe1joaOfgBWZEtQCBrqhcpPQUquLFATYgwyJXWGOpr7iIm3kp6TaHZTgPNJ5rWZa6k6bASF5RIUznvhS9DVAJ94QSqfignt2rWLkpKS4Paa4d4bOVJIT2GG7M19ZC9NRlki4xvnUftRshKyyE3KpaGjn/gYC4vTEsxuVmTY/wQcfAquf8iYgirEOO655x62bNkSLE3R2dl5QS2j+UqCwgxorelo7iMrQvIJADWdNazNWotSivoOI8lsiZCAZaqOU0YvofAauO7LZrdGRLC6urrgIrSA7du3BwvezXcSFGag1z6Iy+mJmHzCoGeQOkcd67LWAdBg72e5JJnBMwSVn4SYOLi7QtYjiAlVV1ePCAgBnZ2dwZXLY1VHraysHFFFNVpJUJiBSEwye7WXdZnr8Po0jfYBSTIDvPZPcPYg3PUDSF9qdmtElBuvOuroKqrRShLNM9DR3AcKspZGRlCosdcAsC5rHS0OJy6vT5LMp9+Gt74HpZ+ANXeY3RoxRW88fYKOpr5Zfc7sZSlcu3X1uMdLSkr49re/PeK2QLI5UK9ovLpF4+3FEE2kpzAD9uY+bLlJxMZHxnDE0c6j2OJt5CXn0WD3zzxayMNHQ73GquWM5XDLP5ndGhElAgXthm+GU1FRwc6dOy+473jVUaOZ9BRmoKO5l9zCNLObERRYyayUkjUKAC9+Fbqb4JMvQnxk9ObE1Ez0jT6cnnnmGXbt2hXMI9hsNrZv3x48PlZ11LGqqEYjCQrTNOT00NMxyNqrl5jdFMDYk/lU1yk+vv7jANR3DJAYa2VRWrzJLTPJsd/Cvsfhmi9CwRVmt0ZEoYnyAsOPBdYvFBUVUVtbG/Z2hZsMH02TPcKSzCe7TuLRnhEzjxZsddT+Dnj+87DoYrjhq2a3RoioIkFhms7PPIqMNQo1nUaSeW2WUU+9oaN/YQ4daQ3P/xUMdhvTT2PizG6REFFFgsI02Zt7SUiOJdkWGSedGnsNqXGp5Kfk4/H6aOwcWJgzjw48Acd+Azf9HSxaZ3ZrhIg6EhSmyd7ST1Z+csQMz9TYa1iXaSSZzziceHyaFQtt5lFPK/zuISh4H1z1WbNbI0RUkqAwDdqn6WzpJ3NJZOQT3F43J7tOBvMJdR0LsBCe1vDCF8E7BHf9h6xajnLRsE1wJJnNn5cEhWno7RzEPeQla0lknHRPOU7h9rkZ6Mvjqd2Nwemoy7MX0OY6h5+F47+FG/8GsorNbo2YgYSEBOx2uwSGEGmtsdvtJCTMTuFLmZI6DfYW46QbKSuZj3YeBeDtIwm86W7gihWZJMdZyUlZINNR+zvgdw/C0lIZNpoH8vPzaW5upr293eymRI2EhATy8/Nn5bnCFhSUUlsAB1Citb6gOpRSqgQoAtBaR1VN2s4WY+ZR5uLI6CnU2GtIjk2mqT2RtAQP9XYjyRwp+Y6w+92DMNgDd31fho3mgdjYWFasWGF2MxassAwf+U/4aK2rAEfg+igP+INB0TjHI5b9TD+pmQnEJUZGR+uo/Sir0tfQO+ijx+mmoaN/4eQTjv3WGDq6/kHIXWt2a4SIeuHKKWzD6CUA1AHlww/6exG1AFrrXVrraqJIZ0sfWUsj46Tr8Xk43nWcRQnGOHrvkIczDufCmHnkdMBvvgCLNsA1XzC7NULMC+EKCjagc9j1rFHHNwFZSqkSpVRU1Zj1en10nR2ImJlHdd11DHmHSKEweJvXpxdGT+H3fwP97cZsI2us2a0RYl4wc/zDrrWuVkqVK6W2jM4rKKW2A4EKVENKqcNz38QJ/GBWny0b6JjJExzm4yOu33NhQUezzPi9Teobpo4+hv/9mUveX3S7aKoPCFdQcACBXdFtgH3U8VrO9yTqMHoOI4KC1roCqABQSu3RWpeFqa2mm8/vbz6/N5D3F+0Wwvub6mPCNXz0FP6ZRf7/qwCUUoGdKapGHd8dpnYIIYSYgrAEhUDiWClVDjiGJZJf8R+vw5iVtMV/PaqmpAohxHwVtpyCf/hn9G2lEx2fwFTuG43m8/ubz+8N5P1FO3l/oyhZSi5mk1KqZLwpxpMtaIwGk7y/nVrrHUqp7VP80iNExIi42kdKqS3+GUljTlWd7HgkC+G97fT/v32s45HOP1z42DjHQlnQGNEmen9+25VStRiTJ6KOUmq7/9+Yc9ei+bMHIb2/qP38+X8v5bPxu4uooDDZiSOaTywhtj2qTyr+99Y5zuEJFzRGg0neH8A9Wuti//2iij/gVfl7OEX+68OPR+1nDyZ/f35R+fnz/y42+383JTM9b0ZUUGDyE0c0n1hCaXvUnlRCMNmCxvmgJIq/SRdx/m+yjvOzAwOi+bMHk78/iNLPn9a6Wmu9w3+1aIzhzSn97iKjeM95k504ovnEEkrbS/xF7KJ2zH0hC/zOlFKblVLl0XRyGZUDKcGYVj5cNH/2Qnl/EOWfP/+XkQfGODSl312k9RQWNH8dqCqMEiDR9k1sMpMtaIxq/jHbLf6rdsb+Jhrx/EMLL0dbPbJQTfT+ov3z5w9kDwxbDzYtkRYUJjtxRPOJZcK2z5eTymjD/kDHXNAY7Ya9vzrOv6diYMorSSNE+TjfkqP5szfcmO8vmj9//hpygTxBHefLAwVM6XcXaUFhspXQ0Xximey9Rf1Jxf+hKhv24YLzCxbHW9AYNUJ4f1sDFYCj9P1tHzYEVu7/fz589oBJ3180f/7KGXnSr4Pp/+4ibp2CfzpYHUbCJFD7aG9g4dtYx6NFiO+t03886sY0RfTynySfwfj7y8RIulbNo89eqO8v6j5//pP/Voy2b9ZaP+C/fVq/u4gLCkIIIcwTacNHQgghTCRBQQghuVwktQAAAVVJREFURJAEBSGEiFL+mUfl412fDgkKQghhojHKUkylxtQ2jBlH412fMgkKQghhktFFFqdRY2r0yuyxVmpPSaSVuRBCiAXDPy12eAmKbcDL/suBOkXVo9bGgFHcz0EYSE9BiGnwd/G7lFI2pdQzUVoET4SZf4y/aNj18knKUIxZp0hrXTnqXyAglAGbhj3n6OtTJj0FIaZBa13pL562FXg0morfibmjta727+EQ3Jd+Nr/hj16INhuLCqWnIMQ0+fcWf4DoKokg5pj/RL0DYzXxZF8eTK8xJUFBiGnyJwnvB8bc7UoIGFFig+FDSeMwvcaUBAUhpsFfS2aHv/Bd0XjbIIqFbfhsIn+PYUROYXSRxUgoHCm1j4QQQgRJT0EIIUSQBAUhhBBBEhSEEEIESVAQQggRJEFBCCFEkAQFIYQQQRIUhBBCBElQEEIIESRBQQghRND/ByciEmgaoeVRAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "/home/andreb/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: divide by zero encountered in true_divide\n", "/home/andreb/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:23: RuntimeWarning: divide by zero encountered in true_divide\n", "/home/andreb/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:26: RuntimeWarning: divide by zero encountered in true_divide\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEOCAYAAACJlmBtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3iUVRbA4d+dyUwqSUiHhJBCC0VIAVQEQWLHgiLYdq0EF1ddK7hi2XVXBbtYIFbUFVHsDSUUC0USQu8khJoeEhJSZ+buH98EI0xCIJmSzH2fJw+TfHfmO7Q5c9u5QkqJoiiKouicHYCiKIriGlRCUBRFUQCVEBRFURQrlRAURVEUQCUEhBBJQojU5r5XFEVxF50iIQghko77foIQIlUI8VArnj4JCGzhe0VRFLfQ4ROC9dP8m02+TwKQUmYA5ccnCxsWnOR7RVEUt+Dh7ADaSkqZIYQoa/KjScBi6+NcIBXIFkJMOO6pGVLKckfEqCiK0hG4fEKwfsIvl1LmWr9PBbJaeDMPBJomiGAAKeXCZtqnAPFCiMYEcfz3iqIobsHlE4KUMlsIkSaEyADirD9rtzdqKWV6S98riqK4iw4xh2B9k54GxFnnBlpSDgRZHwcCpfaMTVEUpbPoEAnBOkz0qfVx3EmaL8Dak7D+erIEoiiKomCnhGBd8pkqhJjZzPVWLwttumrI2lNIFUIENrk+AUhpnDSWUmY3xoA295DdDr8lRVGUTk+0d7VT6xv4JCnlNCHEYmBa0zdl6/U4KeVCIUQa2gSxetNWFEVxsnbvIUgps6WU06zfxtl4s5+ENs4PfywLVRRFUZzMbnMI1uGgKTYu2VwWqiiKojiX3ZadSilnCSE+FUK0tGfAJutQUhqAr69vcr9+/docz94je7FIC7EBsSdc211UhV4niA3xPb0XN9dD4RYIiALf0DZGqiiKy6gthyOHwFQHnl3AvzsYfJwdlU1lRaWYq3VI4YGQ1ewrPVAipTylN6R2TwhNJoGz0YaE0oBZTZqcdFmodfI4HSAlJUVmZWW1Oa5HfnuENQVrWDxh8QnX7vwom62HjrDsgdGnf4PZyRAUBzd8evqvoSiKa9j3O/w0Aw7sgtAz4Px/Q+/zQQhnR3aCvNw9LPrP55iNiRjqivFLqub6u25BCLH3VF/LHj2EVKBx3iAQyAQQQgRaewoL0HYDgwOXhYZ6h1JSU4KUEnHcX2pUV28WbynEYpHodKf5Fx5/Hqz7UPsk4eHZDhEriuJwpTmQ8Ths+wb8IuCyV2DIDaB3zT2882bOoW5HKGbDYHTmLMY/fTOhEWGn/Xr2mENIB+KaLANtLBmxxPq9U5aFhniHYLKYqKirOOFaVKA39WYLxVV1p3+D+POgoRr2/96GKBVFcYqjJfD9g/DaMNi9FEb/E+7OhuSbXDIZbFy/mbk3v0zVnj4ISw3hZ5bwtzcfalMyADv0EKy9gMbyDwub/Dy5yWOHl4cI8QkBoLimmECvP1e3juqqjQkeOFxDuL/X6d0g5hzQeUDOUogd1aZYFUVxkPpqWP06/PaS9oEu+SY4dzp0CXd2ZDaZTCbmPfkapgMxmI0JCJHJVS9MITiofSr2u17qs5NQb21upbimmN5de//pWlRXbwAOHK4muWfX07uBZxfoMRx2Z0DqE22IVFEUu7OYYcPHsOy/cOQg9L1E+38b2tfZkTUrc+XvrJubRYP3IIymvURfEsrFV087+RNPgdskhBBvrYdQUlNywrXIYwmhpm036X0+ZDwBFQchILJtr6Uoin3sXgKLH4fCTdA9Ca5K13r4Lqqhvp53H5mN5XA/pDEevVcWN7xwFz4+3u1+L7dJCI09BFsJwcfoQZCvse0Joe8lWkLYuQiG3ta211IUpX0VbILFj2nDuoHRcPXbMOAq0LluSbelPywl5+PdNHgnYmzYRdx1fRh7UWsOgjw9bpMQfAw++Hj4UFxdbPN6VFdvDpa3MSGE9NGWnu74QSUERXEVh/Ng6X9h06fgFQAX/BeGTXbp1YBHKqv48OHXEHWDEYZIPPzXcfPsezAY7PuW7TYJAbRhI1s9BNASwvaCyrbdQAjoczFkvgl1VeDp17bXUxTl9B0tgV+ehcy3QaeHEffAOf8A79OcJ3SQL+Z/ScmicqT3UDzqt5EydThJwy91yL3dKiGE+4ZTWF1o81pkoDdLthXZ3KdwSvpeDKtf07ql/S8//ddRFOX01FXBqtdg5SvayqHEv8Do6douYxdWVnaY+f9MR1iS0HsIfLpv4a+PTEWv1zssBvdKCD7hZBfa3vbQI8iHOpOFosq60196ChB9JngFavMIKiEoiuOY6iF7Hvw8E44WQ8JlcN5jENrH2ZGd1BfvLqDkZxN4DcWjYSPnT7+E2L5XODwOt0oIEb4RFFUXYZEWdOLPE0kxwVodo7ySo21LCHoD9L4AdnwP5gbte0VR7MdigS2fw9IntfmCniPg2vnQY6izIzup0sJCFs54H5MuEYOugi7xu/jrg/9wWjxulRDCfcIxSROlNaWE+vy55tOxhFB6lOFxbSzAOuBK2PQJ7PkZeqnq3opiNzlLtSWkBRshbABc/6nL1hw63vsvvkXtRn9Mnsno6rO5YMY1xMT1dGpMbpcQAAqrC09ICN0DvTDoBXtKqtt+o/ix4OkPW75QCUFR7OFgtrbEe8/PEBAN4+fCoGu0yWMXtyd3Dz/953NMxkSMFOE/YB9/uesBZ4cFuFlCiPCNAKDwaCEDQwb+6ZqHXkePIB/2lh5t+40MXtqehG3fwqUvgoex7a+pKIpWfG7pk9qHLe8guOgZSLnVpZeQNvX+029QsysCk2EwwpzFlTNvITTMdUrmu1VCCPfVeggF1QU2r8cE+7KnpB0SAsCA8bDxY8hdDn0uaJ/XVBR3VVkAP8/SJo31Rhj1EJx9F3j5OzuyVtm6YTO/vrgEk9cgjJaDdB9hYvxf7bfB7HS5VULo6tkVo85I4VHbS09jgn1ZlVPa9qWnAPFjwDNA+ySjEoKinJ7qMvjtRVjzJlgaIOkmOHeayxafO15Dg4l3npiNLOiFxdgPocvkmpfuIDAwwNmh2eRWCUEIQbhveLM9hNgQH2oazG1fegpaFzZhHGz9Gi59HoyuecqSorik2iOw+g1Y9SrUVcIZE7W9BEFxzo6s1X77eQVb39mAyXswRlMesZd348Ir2rcYXXtzq4QA2sRysz0E6xGae9q69LTR4Otg/f9g+3dwxjVtfz1F6ewaaiDzLfj1Bagpg37jYMwjEN7f2ZG1Wk1tHe898gocGQDGWAzea7nxxbvw9m6H9xQ7c7uEEOEbwbqidTavNd2LcGZbl56Cth46IFpLCiohKErzTPWw7gOt1ERlPsSNgbGPQmTyyZ/rQpZ8u5jcT/di8U7GWL+LM24dwPBzH3R2WK3mdgkh3EcrX2Frc1r3QG+Meh15pe2w9BS0KopDrtMmw1RJbEU5kcWsFZ1b/rS2qazHcLjqTYgd6ezITsnRqqPMm/6qtRhdBB6B67jl1Xvw8OhYb7EdK9p20N2vOyaLieLq4mOrjhrpdYIeQd7ktddKI4DB12pb6Td+DCPvb7/XVZSOTErt3OJl/4Xi7RAxqENtKmvqq4++pOjHCmsxui2cdfdIBiWPc3ZYp8UtEwLAoaOHTkgIoA0b5bXHXoRGQXHa0NG6/8E593W4f+yK0q6k1HYXL30SDq2D4N4w4V3of6VLn0tgS3FJGQtnvIm0JKH3AO/Irdwy4+9tX6HoRG6bEA5WHSQxLPGE6zEhvqzIKWmfpaeNkm6CL9K0XZVxo9vnNRWlo9m7Cpb+B/b+BgE94IrX4IxrXfIQ+5P58p35FP8isXgNxaN+PSMfvIT+A650dlht1vH+Jtqou681IVQetHk9NsSX2gYLhypqiQxspyPq+l8Bi6ZrqyfiRrfPaypKR7F/DSx7CnKXgW8YXPysdph9B9ld3NShA4f45on/YdInYhDlBPXL5bp/3OfssNqN2yUELw8vQrxDOHT0kM3rvcO0Q212F1W1X0IweEHSX2Dlq2pyWXEfB7K0RJCzBHxC4PwntZMEjb7Ojuy0zHv+Teo2B1qL0a1l9D8n0qe3c4vRtbeONWjXTrr7dedgle0eQu/wLgDsKmzj6WnHS74FpEXbeq8ondnBbPjfNfDWWG2eIPUJuGcDjLi7QyaD7dt2MeeWF6jaFY+QJoKHHORv7zzY6ZIBuGEPASDSL5LNJZttXgvyNRLsa2R3UVX73jQoVqt8uvY9bbVRB+wuK0qL8jfAsqdh5w/aIVFjH4NhaeDZxdmRnRYpJW/++1Xk3ijMxkHoZRYTnr+drsFBzg7Nbtw2ISzeuxizxYzeRrnc3uF+7GrvhABw1lT4YDxs/EQbQlKUzqBgEyx/BrZ/qx1iP2YGDJ/SYQrP2ZK9Zh2Zr/2KyXsgRvMBYi8QXHKt6xWja29umRCO7UWoKT5WErup3mFd+Gr9wfZdaQTa7svwQbByNgy5ocMts1OUPyncoiWCbV9rhRxHPwxn/k1LCh2U2WzmnUdfwVzUC4tnXzwMWdzw3N/w8euYvZxT5ZYJIdJXm9Q9WHXQdkII9+NIrYniyjrC2qOmUSMhtJK9X6TB7sXQ58L2e21FcZSi7fDzM1olX2MXrRT1WVPBu6uzI2uTdStXkzV3HfXegzE27KHPNT04b1zn7xU0ZZeEIIRIsz6Ml1KeUN5PCDFTSjlNCJEmpUy3RwwtieyiJYQDlQdIDj+xVkov60qjnYVV7ZsQAAZeBUv+DSteVglB6VjyN2q1hrZ9o00Oj3wAzroTfDr2mHptbR1vT38Z3dFBYIwBryz+8tI9eHm53zxfuycEIUQqkCGlzBVCfCqESJVSZhzXLE0IMQGY0t73b43uft3RCz37KvfZvN47zLrSqKiSc3qHtO/N9QY4++/avoS83yDmnPZ9fUVpbwfXws/PapPFnv7aoogzp4JvOxSAdLJfFi1hx/w88E5BX7+Dntf158KL3KtX0JQ9eghx1q90INf6+HjX2EgSDmPQGeju1539R/bbvB7iZyTQx2CfiWWA5Ju1Qz+WPwM3f2ufeyhKW+1brRVmzFmirRoa/U9tstg70NmRtVll5RHen/Y6wpSIzhCGMWgDtz55N3q965/JbE/tnhCOGwJKAhbYaJZknaxNklLOau8YWiO6SzR7K/favCaEoHeYH7sL7ZQQDN5wzr1aL2HPrx2usqPSiUkJeb9qiSDvV21DWeoTkHJbh1411NT89z6lalkNeA/Do34zZ/xtJGcOv8zZYbkEu00qCyGSgMVSyuzjrzUmASHE+baGlKxzEGkA0dHRdokv2j+ajTkbm11J1CusCz9szm//lUaNkm+G317SdnLGfK+K3inOJSXsXqLNEexfDX4RcOFT2r/TDriZzJaD+UV8/fh7WEQiHnoPPLtv4fbH7nZ2WC7FnquMUm19+rfOHSClXAiUYmNIydrLSAdISUmR9gguuks0lQ2VHK47TJDXiZNivcP8mF/dQHFVHWFd7HDSkcEbRj0A3z8AO36Afpe0/z0U5WSk1P79/fIsHMoG/yi45DlI/ItWcqWT+PLtjyj+VWDxSkFft54LZ1xBTPx4Z4flcuy2yqhJLyBVSpkhhAiUUpajzSvkWpvGA3PtEcPJRPtrPY99R/bZTAgJ3bTu8bb8SvskBNA+fa1Jh8WPanXg9Qb73EdRjmc2actGV7wEhZuhawxc9op27KuH0dnRtZsD+/bz7b8WYDYkYRClBCXs4bp7Ok8xuvbW7jujrKuMZgohcoQQh5tcWgJgHUKaaO0p5NgaUnKE6C7WhNDMSqP+1oSw9dAR+wWhN2gFv0p3Q9Y79ruPojSqr4Y1b8LsRPj8drCY4Mo58Pe11gqknScZfPjCm3z/+O+YPYaga8hi9IzhXHfPbc4Oy6XZY1I5Azhhh4qUMrnJY4fvPThepF+ktvT0iO2EEOBjIDLQm635dkwIoO1FiD1XO0LwjIkdfnOP4qJqDmvl11fPgeoSiBoGF82EPhd1uh3zO7ftYNms7zF5DsYoCwhJLmRSmvsuJT0VbrlTGcCgN9DNtxt7j9heaQTQv7s/Ww9V2DcQIeDC/8Kckdoy1Itn2vd+ins5kg+rX4Osd6G+CnpfoK1wiz6r0y1kaGgw8eEzc6jPi8JsGAgyk6tfSCMoSH3Iai23TQgAsQGx7KnY0+z1/t38WbKtkJp6M95GO65Pjhik1Yn/fS4MmghRJ+6eVpRTUrJL2w2/cYE2LDTwahhxj/ZvrRPamLWe1bN/pcF7AEbTfsLGSsZff0KRBOUk3DohxAfG83v+781WPU3o5o9Fwo7CSob0sPNmnLGPw/bv4eu7YMrPaoJZOT0H1moTxdu+0UqsJ92k7YzvGuPsyOzCbDYz9+GX0JX1RRp7ozdkcd1zd+Ln1zmWyjqaWyeEuIA46i31HKw6eGzVUVMDuv8xsWz3hODlD5c+Bx9fDytf0coDKEprWMyw43tY9RrsW6VVHh15Pwy/A/xCnR2d3WQsWs6e+TuQ3onoG3KJmxjD2EvVXEFbuHVCiA+MByCnPMdmQojq6k0XTw+25tt5HqFRv0sh4XJYPhP6XQahfRxzX6VjqquC9R/B6tfh8B4IjIaLnoHEGzvsoTStUVNTy4ePzMZUORBhiELns5a/vHS3Wxaja28qIQA5FTmMYcwJ14UQJHT3t+/S0+Nd8qxW9O6z2+D2DHWymnKiI4e0/StZ70JtOUQN1cpL9BsH+s79X/qnrxax9/ND1HsnY6jfzhm3J3HmyAedHVan0bn/9ZyEr8GXCN8Icspzmm3Tv5s/n2Ttx2yR6HUOWJXRJQKueFUbOlryb20FkqKAVn561WuweaF2Pne/cdr5Gj2GOTsyu6uuquTD6XMw1Z+B3iMEr9CN/PWxv2MwuPVbWLtz+z/N+ID4lhNCd3+q683sKamiV5iDuuH9LoWht8OqVyF+jHYWs+KeLGbY+aM2LJT3Kxh8YehkrepoUKyzo3OIhfM+5fCyWhq8ktHXbeKMO0Zx1lmXOzusTsntE0JcYBxrd6zFIi3oxIkbdBonkzfsr3BcQgC44D+wdyV8cQekLYeAKMfdW3G+6jJY94G2max8H/hHQuq/tHInnaD8dGvs3HuIn578AKFPxkNXiTFyG7c9cjc6R/TU3ZTbJ4T4gHhqzbUcqjpEVJcT33TjQ/3wNerZcKCcq5Md+KZs8IZr3oM3x2rDR7csAqOP4+6vOEf+Rm1+YNOnYKqFnudoHw76Xtrp5wea+u6djzj0qw7hORRD/XouenQ80bGqGJ29uc+/sGY0TiznVuTaTAh6nWBQVAAb9pc7OjQI7QsT3oaPJsFXd8KEdzrd7lIFMDdo+wbWpGvLRj28YfC12tBQxEBnR+dQxQcO8OUTC6j3SMRAMaED9zLx76oYnaO4fUKIC9Sqb+8u382oqFE22wzuEcg7v+2hzmTG08PBJyr1uRDGPgZL/qUliNHTHXt/xX4qCyD7fa2wYWW+tnnsgv9C4g1uWdPqs5ffpnRDIA2GwciGLEb+81oSetvnPBTFNrdPCP5Gf8K8w1qcWB4SFUiDWbI9v5LB9t6gZss590LJTq0Anm+INuGsdEwWC+QshbXvaucQSDPEj4XLXtYWD9jYMd/Z7d+1i0VPfU+95yCMlny8BlVz650P2udgKqVFbp8QAHp17cXOwzubvd6YBDYcKHdOQhACLp8NNeXw3QPa+baDJjg+DuX0HcmH9R/C2vehYh/4BMNZd2qlJUJ6OTs6p6itb2DuY6/gVdwLs6E/ZvPvpP73NmIjw5wdmttSCQFICEpg3tZ51JvrMepPrAffLcCLED9P1u8v569nOSFA0GobXfMufHg1fDFF27CWoM6BdWkWs7U38N4fvYHYUXD+v7SlxW686TBr9Vqy31iFh3ciOtM+fEeauOGm6apX4GQqIQD9gvthspjYXb6b/sH9T7guhGBIDydNLDdl8IbrPoYPr4JPboLxc7QzFBTXcnivVmU0+wNrbyBEKzCXdBMExzs7Oqcym83Mm/Ey9cV9kMZeGL3WctOLd2H07jzHdXZkKiEA/YO0JLC9bLvNhAAwOCqQjG1FVNQ0EODtxEqkXv7wly9h/rXweZpW4z7lVufFo2jqqmDrV7BhvraBDJr0BsZ1qpPITtcP3yzhwMIc6r2HYGzYTa/rejHmIlV2wpWohABEdYnC1+DL1tKtXNX7Kpttkntqqz7W7TvM6L5OHuP09IMbPtV6Cd/eq30iHft4pzv5yuVZLNqb/4b5WjJoqIausTDmEThjEnTt6ewIXUJdbS3vTpuNrB6EMETS4JXJTS/cg5+P6hW4GpUQAJ3Q0S+oH9vKtjXbZkh0IHqdIDOvzPkJAbTho2v/Bz88pNW/L9kFV6VryUKxr5JdsOFjbVioYj94+sOga2DI9dBjuNor0sQnH31JxY+HMXsnY6jfxrCpwxky/FJnh6U0QyUEq4SgBBbuXNjsYTk+Rg8GdvcnM++wE6Jrht4Al74AoQmwaBq8eZ428Rw+wNmRdT6H82Dz59pX4SYQOogbY60yeqmWoJVjaiqr+GD665gahqD3ENB1PbfMvlsVo3Nx6m/HKiE4gVpzLXlH8o7tXj7e0Jgg3l+91zkb1JojBAxP085O+GyylhQufEqbV1CfVNum4iBs/RI2fwYH12o/ixqqnTnQ/0rw7+bc+FzUF/M+oWRZPQ1eKXjUb+S8hy6kd/8rnB2W0goqIVglBCUAsK1sW7MJISUmiLd+28PmgxUk9wxyZHgnFzca/rZCK4b33X3aCVrjXtQOTVFar3y/9me35UvYt1L7WbfBWmG5AePVvEAL8g8V8uXj87DokjDoKtB128ptM+7GQ6/mtjoKlRCsYgNi8dR7sq10G+PixtlskxKjTSxn5h12vYQA4BcGNyyEzDch41/w2plw3gwYNlmd0dwcKaFwM2z/Tvsq2Kj9PDRBmxwecJXbbhw7FXOeTUe3NQCLZwq6umxGTr+KhL5XOzss5RSphGDlofOgT9c+bCnd0mybED9P4kJ8ydxTxh3nuuh6cp1Oq5Xf92L49j748WGtVs4FT0Kfi9QwEoCpHvavhu3fw47vtPLSCG1C+Px/a5VFVRJolW3bc/ll1heYjYnoKKbrwH1c//cHnB2WcppUQmhicOhgFu5cSIOlAYPO9ifqoTFBLNpSgMUiXbsue2C0tjR1xw+w+DFt30L0WTDyAeg11v0SQ2mOtmt49xJtqWh9Feg9taG2kQ9oCdTPBVaPdRBSSl7/92sY9nbHZBgMpizGP3MLoeGhzg5NaQOVEJoYHDaYD7d9yM6ynQwIsb1SZ2hsEAuy9rOjsJKEbv4OjvAUCQH9LoHe52vlE357Cf53NUScAef8A/pd1nk3TB0t1eYAcpfD7gxtlRBAYE9td3f8WC0ZqGW6pywzayNrZy8F7zMQlkN0P9vE+JsecnZYSjtQCaGJIaFDAFhfvL7ZhHB2fDAAK3aXuH5CaKQ3aPMISTfBpk+0xLDwVq2kwpDrIPGv2iqljqyyAPaugLwV2klzxdY9JQZfiB0JZ96p9YyC4tyvd9ROauoaePOx2XiW9MLi2R+9LovrX/kbvv4OPElQsSu7JAQhRJr1YbyUcpqN6xOAciBJSjnLHjGcjgjfCCJ8I1hftJ4bEm6w2aZ7oDexIb6s2F3C7SPjHBxhG3kYIfFGGHyd9qk5+31Y/QasnA0Rg7QeQ8I4COvv2m+adVWQvx4OrYOD2XAo+48egNFPmws44xroOQK6J7p1Ebn28uvPq9j6TjZ67yHoTHuJuyyc869UvYLOpt0TghAiFciQUuYKIT4VQqRKKTOaXE8CkFJmCCHihBBJUsrs9o7jdA0JHcKG4g0tthnRK5gvsg/SYLZg6IhL6nR67eCdPhdCVZG243bbN9p5C8ufgoAeEHOO9tXzbK0cgzMShMUC5XlQtB2KrV/5G6B4ByC1NgHREJmonRHRc4Q2HOZGR03aW21tHe/OeAVR3h9pjMPgvZYbX7wLb1WMrlOyx/+cOOtXOpBrfdzUJGCx9XEukAq4TEIYHDqYRXmLKDhaQIRvhM02I+JD+HD1PjbsLyclxgWXn54KvzA4+y7tq7JQW3WTuxx2/aTV6AHwDNB6EBEDIbiXdrJXYLT21ZYdumYT1ByGqgJt/X/5Pu2rYp/2ib9kN5hq/mjvH6Xtwh4wHronaZ/+/dQkpr389tMytv1vNxbvZIwNu+j3lwRGpqpidJ1ZuycEKWV6k2+TgAXHNQkEypp8H9zeMbTFkDBtHmFD8YZmE8JZ8cEIASt2l3b8hNBUl3Bth3PKrdr6/OLtsG81FGzS1udnv68VcGvKMwC8A7UjH727aglC56HNW+gM2hkAplow1UFDjfZVcxiqS6HWRjlxg4/WQwmMhphRENZP2xMQ2ler9KrYXW11NW8/NBtd3WCEoRv4reXGl/+Bt6fay9LZnXJCEELESCnzWtEuCVh8OsNB1jmINIDoaMfutO0b1BcvvRfri9ZzYcyFNtsE+hgZ2D2AFTkl3JPa26HxOYwQEJagfTWyWKCqEMr3ahVWy/fC0RLtDb7x62gxWEzawfGWBi05eHhp4/geXtqqnsBo7ShQn2DtyzcUAntoK4B8gl17/qKT+/DdhdQsrwLvoejrtzJg8lBGnGN7o6bS+bQqIQghJqN92l8MZAshrpJSfn6Sp6U2M2FcDjR+rA4ESo9vYO1lpAOkpKTI1sTYXgw6AwNDBrK2cG2L7c7uFcw7v+3haJ0JX083GbPW6bT6Pf7dIPpMZ0ejtKNd+wr46b/zEDIJvYcBQ9gmbpoxFU+j6hW4k9bOiOYC0wEB3MGJ8wJ/IoRIa0wG1klmhBCNhxEvaPL8OCDjxFdwrmERw9hetp2Kuopm24zqHUqDWbIq54R8pigdhpSSz9/9mOWPLQExFH39Zs65tw9p/75HJQM31GxCEEL8JIR4QwhxFRALSCnlZ1LK6VLK51p4XiowUwiRI4RoWit6CdqLZDdpV+5KK4waDe82HIkksyCz2TYpMV3xNepZtqPIgZEpSvs5XFzEnMnPkb86BKnzJKDXbqa8dy8DBiac/MlKp9TSWMc0/lgF1AtYaIDvvxMAACAASURBVH2Dz0RbVrre1pOsS0y72vh5cpPH6cdfdyWDQgbh7eHN6vzVpPZMtdnG00PPiF4hLNtehJRSHQ6udBj1JgsfvvwOps3+WDyT0ddnc/6MCcTGxTg7NMXJmk0IUsp11oefWb8AEEIkAvGAzYTQGRj0BpLDk1lTsKbFduf1C+OnrYXsLKyib4Taram4vrzcPXz378/AKwkjRZhidzP5/vswenTA/TRKuzvl2VBrolh30oYd3PCI4Ty/9nkKjxYS7htus03jUZpLtxephKC4NLNFMu/p12nI6Q6eg/FkLVc/dytdg11q1bfiZOpjQTOGdxsO0GIvISLAi/7d/NU8guLSfv99HW/e8go1+xPQmavwTizm9jkPqmSgnMBN1kueur5BfQn0DGR1/moui7+s2XZj+oUy5+dcKqobCPBRqzIU12E2m/n0qTcoz+uJxbMfnoa1XPvsFPy6qA1+im2n1UMQQnT6f1E6oePMbmey4uAKLNLSbLvUhHDMFsnSHYUOjE5RWrb6l1W8c2s6pQf7ozcV4TeyhttnP6iSgdKiZnsIQojmjj0SaCuPbG/j7URGRY1iUd4itpZuZWDIQJttBkcFEuHvxaLNBYxPjHJwhIryZ4crq/lgxmwM1YPAGIPZkMkl/7mDyJAAZ4emdAAtDRmFAB8DFcAEXHADmb2NjByJTuhYvn95swlBpxNcOCCcBVn7qa434WNUo3CK45ktkvc++BrL0mL03kPxqN/JgFsHcva5J1SfV5RmNTtkZN2Atl5KuQdYK6Vc1/gFtFzXoZMI9ApkSOgQfjnwS4vtLhwYQW2DhV92FjsoMkX5Q9bufF6Y8hR1K7wwGyLwD9/Ere9M5uxzz3F2aEoH09qPs8nWjVeN5ayTgKX2CsqVjIoaxUvZL7VYDntYTBBdfQws2lzARQO7OThCxV1V15t4efZHhGyU+HifhbFuK6kPjCZmgCpGp5yeVk0qSymfRduMNg2Ia6l0RWczusdogBZ7CR56Hef3D2fJtiLqTGYHRaa4s2Xrc5lzxywCdnTD4hFA9967uf3dvxMzwPbQpqK0RqsSgrXaaTzaUNGnQojz7BqVC4kLiCPSL5Ll+5e32O7igd2orDPx684SxwSmuKXaBjPPPDWX3S+txtN4Jp7mLVzx6GDG35928icrykm0dtlpjpRyOtpcQvMlQDshIQSp0amsyl/VYvXTc3qHEORr5Iv1Bx0YneJO1m3N4c205+myNx4pjMQlHuC2t/9BREyss0NTOonWJoRka68gVggxBEg+2RM6k4tiL8JkMbF0X/PTJga9jksHdSNjayGVtQ0OjE5xB+nPzCHz2Wx0xiS85AYmPDOCi6f81dlhKZ3MqcwhXABci3bwzbN2jcrFDAgeQJRfFIvyFrXY7srE7tSZLPy4RW1SU9rHgdw9zL3lBRry+iCkifgRpdyWfj9B4WrxgtL+Wnti2o/AFCllnhDiaiHEAinlJDvH5jKEEFwYcyHvbXmPstoygrxsn6OcFN2VHkHefLX+IBOS1SY15fTtKTnKb+nvUb27OybjIHQyi4kvTiag6wmV5RWl3bR2yCgdmCWEeFpK+RnwiR1jckkXx16MWZrJ2Nv8/jwhBFcOiWTF7hKKKmsdGJ3SmXyxaAWL7nubI/sS0JmP4JdSxt/mPqSSgWJ3rU0IUko5EcgQQvwEpNgxJpfUp2sfYvxj+GHPDy22u2JIJBYJ32zId1BkSmdhNptJf/gFCheWIb36YhFrGPb4Rdw02W0644qTtTYhxAshzpNSLpFSXoBWz8itCCG4PP5ysgqz2H9kf7PteoX5MTDSny/WHXBgdEpHJqUke8Uq3r31TRoOD0HfUMjQSUbuemM6g+Nsn8WhKPbQ6kllKWXTJTZz7BSPS7s8/nJ0QscXu79osd01yT3YfPAImw641Qpd5TTU1dXz0l0z+f29chqM0dR5/M71c25k6Ngxzg5NcUPNJgQhxAIhhL8QIlEI8ZP1+wVCiE+AxQ6M0WWE+4YzovsIvsr5CrOl+R3JVyZG4mXQ8dGafQ6MTulIahvMvDz3U95P+wCjaRj6ujwqhjVwzyvT8fX2cnZ4iptqqYcwXUp5BChHW2E0yfo1EZjomPBcz1W9r6KouoiVh1Y22ybA28C4M7rz9fqDVNWZHBid0hFUVh5hzt9nYVjrj9kQBgHruO3dNKZPvgK9zu1GYxUX0lK10z1Nfr36uGud/kzl5pwbdS5BXkF8vuvzFttdPzyao/Vmvl5/yEGRKa7OZLbwxNNv8/Hfv8QghmOo30HMDaHcOfN+9Hq9s8NTlFZPKuc2/ca6W9ktGfQGruh1Bcv2L6PgaEGz7RJ7BNIvogsfrdnrwOgUV1VRVkr6354lNC8ai74LltCNTH7vbi4ac6azQ1OUY1qbEK4VQuxqMofwqT2DcnXX9r0WiWT+9vnNthFCcP3waDW5rPDKs2+x4L6fkLqhGBo2Me7RIdz15D+cHZainKC1CWGBlLJ3kzmEO+wZlKvr7tedsdFjWbhzIdUN1c22u2JIJN4GPfNW5TksNsV1fLtiE7NvnoU+Jw6EB+Whm7n9rXuJjFXF6BTX1Nplp581PhZC+ANuf0DrjQk3cqT+CN/mfttsmwBvAxOSo/h6/SG1c9mN7C+r5oUnZnPorR3oPJOQdVlc/t9zeOTJu9GpSWPFhbX6PAQhRJa1ptFCYGgrnpPUwrWZ1l87bBH3xLBEEoIS+GDrBy0uQb31nFgaLBY+XKXmEjo7KSW7tu3gm3vn4FkwACHr6Xl2KZPfvJ+I7qoYneL6WjtkVCalTAFmWXcqt1jLSAiRCrzZQpM0IUQOx01WdyRCCG4bdBt5R/JYvLf5bRmxIb6kJoTzweq91Dao09Q6K7PZzPPTX2Dpc9vAcxCCTM569FzG3TQJTw+1gkjpGFqbEBBCPGD99XZgbEttpZQZQFkLTa6RUsZb23VYqdGpxAbEMnfjXCzS0my728+J5XB1A59nq8NzOhuLRbJk+WreuvUNvCsS0ZkqkIOLmTpnGgPjI50dnqKcklOZQ8iQUi5Bq2NU1cb7JgkhUoUQD7XxdZxKr9MzedBkdpfvZtm+Zc22GxYbxBlRAbz1Wy4Wi3RghIo9VdbU8eSdT7HrwxIsxl4cNa9i0qtXc+ed1zs7NEU5LS2VrmharmIB8E/rr+cD97flplLKWdbeQbB1eKnDujj2Ynp06cGcjXOa7SUIIbjtnFhyi4+yeJs6PKcz+P67pXx8x3uEyLPQ1efjO7qB2196EP8ufs4OTVFOW0sH5My09ghOIIRIPN0bCiEmAEgpFwKlQJyNNmlAGkB0dPTp3sohPHQeTB0ylYd/fZjv93zPuLhxNttdOqgbLy7eyStLdnFB/3CEUKtNOhopJfuKK/j40dfwtyQhDFE0eGYyZfZ9eBgMzg5PUdqspdIVNpOB1eFTvZEQItD6MBdonDuIB7Js3DtdSpkipUwJDQ091Vs53CWxl5AQlMDs7NnUm+tttvHQ6/j7eb3ZcugIGduKHByh0lbb8o/w8Kz3+fHBz/ETZ6Gr20Pl2ZLJs+5XyUDpNFq77PTpJl/PcJKdytZeQEpjb8BqCYCUMhuYaL2WY/2+Q9MJHfcm38uho4da3L185ZDu9Az24eUlO5FSzSV0FDkHivnq4VeIzolAeoSgC1rPNa/dyoO3jsPXs1Wn0CpKh9Daf80CmGt9HAdkttTYOhy08LifJTd5nH4KMXYIZ3U/ixHdRzB341zGxY0j2Dv4hDYeeh13junFQws3snR7EWMTwjlaZ2LOzzlMHhWHv5f6pOlqvvpgIUVLaujidTYetZs49/5U+g663NlhKYpdtHaV0XQp5R7r1xJOY8jIHTw07CFqTDW8sPaFZtuMT4ykR5A3L2XsQkrJr7tKmL10Ny8u3unASJWTKSws4o3JszjwWyAWnS/1IRtJe+8e+g4a4OzQFMVuWjtk9JMQ4scmK46a3YXszuIC4rh5wM18nfM1awvX2mxj0Ou4Z2wfNh2s4LtN+RwqrwHgg1V72V3U1tW8SnuY/9p7fDV9KRZ9CpaadRwcHUjao393dliKYneiNWPZQoixJ5lktpuUlBSZlXXCvLPLqm6o5sqvrsTHw4dPLvsEo954QhuzRXLpK79SXW9mdN9QFmTux+ihI7lnV967ZZgTolYqqhuoKS3gqycWYDIkYagrZnvQQW685xaSors6OzxFOWVCiLXWChOt1tqdyn8qz2g9WvNpdz4XoTk+Bh8ePfNRcipyeH396zbb6HWChy9JYF9ZNZ9k7Seqqzf3jO3N8h3FLNuuViA52v6yah598Hm+fnwNJo8h1Nf8jpjUh/88fbdKBopbaW1C6GodLoqxfv8wkI6NPQQKjIwayVW9r+LdLe+yoXiDzTajeodwTq8QahssRHb14a9nxRAX4suT326lzqRqHjmCxSL57Kff+eYfc+lrHoaw1LLebzPjn7+byRcl4mNUK4gU99LahLBWSjmJPxJAnPVoTXXySzMeTHmQMJ8wHvntEarqT5wbEEIw/eJ+AEQGemP00PH45QPILTnKG8tzHB2u2zGbzbww7XlKFhSA1wAEmUx4cRzP/edOYkJ8nR2eojhFaxNCsrWoXZwQIhaIV+citMzP6MfT5zzN/sr9/GvVv2zuOxgYGcCcG5OZMkrLs+f2CeWywd15fVmOmmC2k6o6E9uzs7VidJVJ6Exl+I+oYuqcaXTtGkCAt1r6q7iv1iaEdGCPlPItINA6UTEFKLdbZJ1ASkQKdyXexaK8RSzYscBmm4sGRvzpE+lj4/rjZdDxyBeb1Oa1dlZaWcOsu55m+ev5WIy9KKtbwXnPXMmNN13l7NAUxSW0dh9CReMqIynlOiGEv5TyWSnlUvuG1/HdOvBWRkaOZFbmLDaXbD5p+9Aunjx8SQK/7ynjk6z9DojQPfz47VI+u/MDwnUj0NUfZGX4QR55ewbx3dSksaI0au2y00RgEiDRdi0nSikvtHNsQMdbdmpLeW05k76dRL2lno8u+Yhufi2fnmWxSK57czVbDx1h0b2jiAz0dlCknc+K7Yf45fl3CCQFIc0c9tzInTPvxaLTEehz4pJgReks7LnsNBWtdEW69dcOfbCNowV6BfLa2NeoM9UxdclUKusrW2yv0wmeu2YwFil58NMN6gyF07Bydwk3T3udrU//iL/+bERdDgUpZqa/Oh1/Py+VDBTFhlNZZdRYumIP0PyZkYpNvbr24sUxL5JXkce9y++lwdzQYvseQT7MGNeflTmlzFuV55AYO4uNufms+O9chh3uhcUjiDr/bHpPv4r7br5YlR1XlBa0dsjoJ6Ar2rGYFWhDRr3tHBvQOYaMmvpq91fMWDGD83qcx3Ojn8Oga35Vi5SSW9/LZGVOKd/dfQ69wro4MNKOp7SqjuxvviPvx2rqvbpjqN1It7+OYMyIRPxUVVLFzdhzyGimlHKolPJCKeVE4I5TD08BuKLXFTw87GGW7l/KQz8/RIOl+Z6CEIKZV5+Bn6cHd/5vHTX1asNacypKS/jo7pfYudwfqfPmkPdarnltKpedN1QlA0VppdauMjq+jpHaOdUG1ydcz7Sh08jYl8G0X6Y1e6gOQJi/Fy9OGsLOokqe+HqLA6PsOJ564jUWPLAUnXEosnYd88LN3POvu+jqq+YJFOVUNPvRyVrVdDLaqWYz+aPktQASAYcMGXVWN/a/EYlkVuYsyuvKeWnMS/gb/W22HdUnlDtH9+LVZbsZHhfEVUlRDo7W9SzanM/MBb9y/aEDBBiGAMVk6TL556yp3BHii4e+tZ1fRVEaNTuHIISIlVLuse5MxjqZ3HgtUUq5zhEBdrY5hON9m/stj654lBj/GN5IfYMI3wib7UxmC9e/9TubDlTwzV0j3HY+4fDRepZsL2LLxx8TURqHyeBPXW0mPdMuJzqyG0NjgpwdoqK4hHadQ2hMANZfY4UQMUKIIUKIB1AH5LSbcXHjeCP1DfKP5jPp20lkFtg+jM5Dr2P2dYn4GPWkfbCWipqWVyl1VnM+/5milxcSUpmEzlKNcXABBy+9hMuGJ6hkoCht1Op+tZQyD3hTSvkc6oCcdnVmtzP53yX/w9/oz+SfJvP+lvdtlq0I9/fi9RuS2FdazV3z12EyW5wQrXNk5hbz1L0zCVpeA179qTOtZvDDY5g89UaenzgYo4caIlKUtmrt/yIhhDgPcMohOe4gPjCe+ZfOZ3SP0Tyb9Sx3LrmTkpqSE9oNjwvmySsH8svOYp7+Yfufrv28s7jTrUT6ZWcx61evYf2/PiWgZig6Uyl1Z5Rx31v/ZFjfSGeHpyidSmsTQhlwPvCUEOJq1DkIduFn9OPF0S8yfdh01hSsYfxX4/kx78cTegvXDYvm5rNjePu3PSzI3AfAzsJKbnpnDWkfZHWa8xQ27C1h8cw5rH6rGJMxlhqxGstNZ3LLraoYnaLYQ2sTQgra6qKJaGUrHDKh7I6EENyQcAOfjPuE7n7deeDnB5i6ZCr7juz7U7sZlyYwsncIj3yxmV93FbO9QCuH8euuEv7x8foOPZz04uKd/PfVT8h87AsiDCPQ1R/gc589hP31eqacl6DKTiiKnbQ2IeRIKaejlbCoQCtyp9hRXGAcH17yIQ8NfYh1Reu48qsreXXdq1Q3VAPaJPNrNyTRK8yPOz5YyzcbDqETMP3ifvywuYCHP9/UIWsgHa06SvX8Twja4I/J0I1iVnHTW7fyz6kTuHZotLPDU5RO7VQOyDkPbbXRECDZjjEpVgadgb/0/wvfXPkNF8RcwNyNc7n484t5f8v71Jpq8fcyMO/WYQT6GFm8tZCewb7ccW4894ztzadrD/Dvb7d2mDMVPlt7gFumvcb8v39KoPFsPBp2syq2Ep8rr8bTaGBM3zA1cawodtaqWkYAQohn0OYO1lhXGjlEZ9+HcCo2FG/g1XWvsjp/NWHeYdw88GbG9xpPQTlcM2clZ/cK4bXrk5BS8p/vtvH2b3u46ayePHH5AJcs6ra94AheHnrCPM28fd/r6EUyevNRtuk3cfd/7iEuzPZGPaXzamho4MCBA9TW1jo7lA7Dy8uLqKgoDIY/10U7nX0IrU4IzqISwokyCzJ5dd2rZBdl42fwY3zv8VzQ4ypiA3oeOwJSSslT32/jzV/3cMPwaJ68YiA6nZYUvt+UT2J0IN0CnHvOwplPLaFX/jrOrexGvVc3PGo3YL4smf2iK0+NH+iSSUyxrz179tClSxeCg4PV338rSCkpLS2lsrKS2NjYP107nYTQYtUvIcRYtLIVcWjnJ2cDT0spPz+1sJX2NDRiKPMunsem4k18uO1D5m+bzwdbP2BYxDCu6HUFqdGp+Bh8+OclCXjodbyxPAeTWfLUVYPYW3qUqf/LJqyLJ+/cPJSBkY47FttikXyz8RCjeodSVlzMpB2ZGL1SkLpyNvE7SbdOZOqoeIfFo7ie2tpaYmJiVDJoJSEEwcHBFBcXt8vrNTsoa11eej5wjZQySEqpRzsoJ14IcXu73F1pk0Ghg5g5aiaLrl7EnUPuJP9oPo/89gijPxnN/cvv5/s933PHmG7cfV4vFmTt56752WTv047Brm0wM2nuKpbvKHJYvKtyS7nn4/XMmP4SS/61EqP3MLzkRiJu7kN50mgu6G+7bIfiXlwhGaSnpxMfH8/ChQtJT09n2rRplJe3fIR8RkbGsTbZ2dmcf/75p3VvW88tLy8nOzvbZvv2/PNqqYcQZF1ZdIx1hdGzQogHT/bCQogkKaXN34EQYgJQDiRJKWedSsDKicJ9w7lj8B1MOWMK64rW8W3utyzbv4yf9v6Eh/AgOSKZi87ux6Ls/SzdHoVRr+f7e0aS9v5abpuXxVPjBzKpyQqeg+U1BPsa8TLo2xxbvclCTYOZAG8Dm7bu4oFD2xA+KSCLWKPL4uUX/kGgj5HLR7X5VorSblJSUkhNTWXChAmA9oY8efJkPv30U5vtGxNBYGAgAHFxccTFnd52LVvPDQwMJDc3l6Qk+xaJaCkhlLZwrcXy10KIVLShphNWIwkhkgCklBlCiLiWEodyaoQQJIUnkRSexIwzZ7CpZBPL9i1j+f7l5FTMwzcWpNkTH9mLb/bmcddlfZi3XM+0zzaSW3KUhy7sR1WdibHPLycm2Jc3bkwmNsS3TTG9smQX6b/kcpduF8aDPRDeg/E3rufoFWPobT5D7SlQXFJGRsafPqUHBgY2+wkdtB7FQw891Ozz23LvRklJSSxcuPBYkrKHlhLCUCFEbnPXgGbnEaxv9mXNXJ7EH0dw5qINQ6mE0M50Qsfg0MEMDh3MP5L/QUlNCVmFWfy4ewXbKzbwxoY3kEjQQ3B/f97PDWfxvB4M69GXBmMtueWhXD67imevSeKigRFYLJJ7Fqynd5gfU0fH46HXYTJbWLq9iFF9QpvtTezaspl79x/A4jMInfkguyILeOnJ+xz8p6EopyYzM/OEN96ysube0iAn58+fkTMzM4mPjycjI4Pc3FzS0tLIzc0lIyODlJQUkpKSmDJlCnPnzrV57+OfC1rPYebMmU5LCNcAgWg7lI83Fnj4NO8ZiFYKo1Hwab6OcgpCvEO4KOYiLoq5CIDqhmp2HN7B1tKtbCvdxuoDmymo+Y3vDi7Fp3H0SAru/92H/6wPpmdABGuLTPx4qAtf5HblmqRe5B+WzF9dRHf/QO46rx/9wrvy05ZiftlZyh0je1H9xUoSc3ti8eqHp3Et86IGMvGsM533h6AorZSbm/unYZvs7GxSUlKOPS4rKyM1NbXZ52dnZ/Pww9pbZOMwU+OwUmNiaS7B2Hpuo5aSUntoMSE0d+aBECLRTvEoDuJj8CExLJHEsD/+KtfsKeVv838hLOgIaWP9OVB5iJ927GJX6SEqqgvQ+1bgaThKCSbe2Kw9x7uHVgv932v/eO3oI8HsmzkJk1cCxvocas7yZ+rfHkStRFBOxb++2cLWQ0fa9TX7d/fn8csGtNimvLycoKA/l1KfO3cu06ZNA2DBggUMHTq0xdcIDAwkMDCQhQsXHhv+SUpKYu7cuaSlpbU4pGTruY2Oj6u9NZsQWjoAp42H45QDjb+rQGzMVQgh0oA0gOhoVa7AUYbFBvPz/ZdhMluOje3flQRfrT/II19sJsTPyM8PjiG3uJz7Fv7OhkOFjO3vz41nR/Duyl2s3FXIDXmScI9BSA8zhbVr+DAigZfPbvk/j6K4kqysrD9N3ubmaiPnjT2CSZMmHfuZLdnZ2ccSRmZmJlOmTDnW42j8hJ+dnU1qauqxn5eXlx+bp2juuY7gsNPHhRCBUspyYAFasTzQ9jdkHN9WSpkOpIO2Mc1RMSrYPJD+iiGRnBUXTJ1JK5gXFxrIwrTz+Wr9IYbHBRHV1YfeBwNZnLGBo8ZQQsz7SL1/DFsNwynP3M85vUIc/dtQOoGTfZK3h9zcXGbOnElcXBwLFy4EtGEaW2P9Lb1GY/IIDg4mOzv72Lj/0KFDycjQ3vKysrKYOHEiAMnJyeTk5LT4XEewy05l67LSN4HJUsqF1p+tlVImWx+noU0ox1nf/Juldiq7ttrKWn6Z9T27i7pgaDjK8KF6Bt0xziXWkisdz7Zt20hISHB2GC1KT08nJyeHhx9++Ngy0/T09GOTv6ejNb2A3NzcZhOErT+3dt+pfLqsSWDhcT9LbvK4xSSguD4pJVu+zGbV9/nU6/yJ0e9l9COX4hutNpcpnZutN/6JEye2aUnoyTa9AQ7pLThsyEjpPIr3lLH0ld8oqfEjoK6U1Eu6EXvtbc4OS1GcpnEiuHEu4FSdbMOZo+YRVEJQWq2+xsTKd35n68YaPEyQ6L+VYU/dgEdQV2eHpihO19Iy1LbqdJPKSsdlsUi2Lc9j9cId1Jo9iKpYx4jbzyQk9XJnh6YoSjtSCUFp0f5tZfw6bwOHyyX+FQc4q3cVfZ+fit6vbSUtFEVxPSohKDYdLjjKio+3sXf7EbxqShhStZIh//wLvsnqsDxF6azUmYTKn9RU1fPLxzuZ/6/VHNxSRHzul1xyRj5nfTRTJQPFbXSk8tftSfUQFECbMF6fsY/1Gfsw1ZnofvBX+hp3Efvyo3i5+LpwRWlv7lr+WvUQ3Jyp3sy6n/bxwYxVZH6XR9eCDQxfN5NzLoui38fzVDJQ3NLplL9uusrInuWv7UklBDdlNlvY/MtBPnx0FSs/302Xyr2krJ3JUI9MBs5PJyRtMuK4Q7sVxV1kZmae8Gn8VMtfl5WVkZGRQXq6tg83NzeX9PT0Y4llypQpzd77+OeC1nNYvHixzee0FzVk5GbMZgs7fy8g6/s8jpTUEtKljr7b3iGoMpfQB+6n67XXInTqc4Li3porf52bm0t5eTkZGRlMmDCh2WGhzlj+WulETA1mtq/MJ/vHfVSW1RIcZmTo0UX4Lf8Gv5Ej6fbE1xgiI50dpqL84YfpULCpfV8zYhBc/EyLTVoqf920fMTChQv/dEpaU52u/LXSOTTUmdny60HWLd5HdUU94TFdGBKYg+GDZ9H7+BA+8xn8L79cFaNTFKuTlb8GjvUQbFHlrxWXU3u0gc2/HGTDkv3UVjUQ2TeQUWN8Yc6/qd++nS4XXUTEjEfwCFGlqRUXdZJP8vbQmvLXGRkZpKamNvsmrcpf25Eqf31qyouq2bj0ANtWHsJUbyF6QDBJY7th+OEDSt95F31QVyIeewz/01wBoSj25OrlrzMyMo4ljPPPP//Ym7Uqf624DCkl+bvLWZ+xnz0bS9DpBH2GhTN4bDQ+BdvJv+8mjuTlETDhasIffBB9QICzQ1aUDik1NdVmETtV/lpxOlODmZzsYjYu3U/R3kq8fA2kXBzDwHMj8dKbKH7hefZ+NB9DZCTR77yN79lnOztkRemUVPlrxWnKC6vZ8tshtq/Mp/ZoA4HhPoy+oS99hkdgMOqp+uUXch9/AlNBAUE3/ZXQe+5B5+Pj7LAV9uRQRgAAEMNJREFUpVNT5a8VhzGbLeRtKGHzLwc5sP0wOp0gdkgIA0ZFEtWnK0InMB0+zKFHn6Hiq68xxsfT86P/4ZOY6OzQFUXpIFRCcHFl+UfZsbqA7avyqT5Sj1+QJ8MvjyNhRDd8AzwBbQ7hyKJFFDz5H8wVFYRM/RvBd9yBzmh0cvSKonQkKiG4oJqqenZlFrFjdT5FeysROkHPAUEMGBVJ9IBgdLo/9gw0FBZR8OS/qcpYgteAAUS//RZe/fo5MXpFUToqVaPARTTUmdm9tojv39jIe9NW8OuCnVgsknOu6c3Nz4zg0jsHEzMo5FgykFJSvnAhuePGcfTX3wh78AFiFnyskoGitANV/lpxuIY6M3s3l7J7bRF7N5dgqrfg42/kjDFR9D2zGyFRfjafV79/P/mPPUb1qtX4pKTQ7T9PYoyJcWzwitKJuWv5a5UQHKyuuoF9W8rIWVd8LAl4dzHQ78xu9EoOo1vvwD8NCTUlzWYOf/ghRS+9jNDpiHjicQInTlTF6BSlnZ1O+eumdY3sWf7annsRVEJwgPLCavZsLGHvphIO7a5AWuSxJBCfHEb3FpJAo7rdu8l/ZAY1Gzbge+4ouj3xBIZu3Rz0O1AU95KZmXnCG++plr+Oj48nIyOD3Nxc0tLSyM3NJSMjg5SUFJKSkpgyZcqfSmK09FzQeg4zZ85UCaGjqa8xcXBXOQe2l7F3cykVRTUABEf6knhBNLFnhBAW43/SJAAg6+speestSt+Yg87Xl+7PzsJ/3DhVjE5R7Ki58tfl5eVkZWWRnZ1NWlpas5vQVPlrN2ZusFCwp4ID2w9zYHsZhXmVSItEb9AR2TuQM8b0IGZQMP4h3qf0ujWbNpH/yAzqdu7E/5JLCH/kn3gEB9vpd6EormXmmplsL9verq/ZL6gf04ZNa7FNS+Wvs7Kyjp2L0NKYvip/7UZqKuvJz6mgIKeCgtwKivZWYjZZEALCYvxJujCaqH5BRMT542HQn/LrW2pqKH71VcrefQ+PkBCiXn+NLuedZ4ffiaIoxztZ+evy8nLKy8ubTQaq/HUnVlfdQMn+Kor3V1K8v5LCPUeODQHpPARh0V0YNCaKbvEBRPYJxNOnbcdOHl2zhvxHH6Vh7z4Cr7mGsIceRN+lS3v8VhSlQznZJ3l7OFn568aqpklJScyaNcvmATmq/PXxLyrEBKAcSJJSzrJxfaaUcpoQIk1KmX7iK/zBUeWvzQ0WyouqOVxQzeGCo5Qe0JLAkZLaY23+3969x7ZV3XEA/548nbRpnOZNk1Z1ShljqMK0wBhIIFyGxBgDknZqtQnxSKVVYxJCeKMP2tAEYsqgUB7pKg06VgYxEmxdqZTwGkIaUjASjIeQ4gJtSdLErpt3nMRnf/jaOK4dP2L33mt/P1KUuDePc+Te+7v3nnO+d5GxEFUrSlDTUIpaUykqV5QkdQUQ8e+PjuL03r3w/ONV5NfXo/aRFiy66qqU/G4ivdB6/HXgVlHgwB04c2f8dRRCCDMASCm7hRAmIYRZShk+X6tZKRqRnzKdJlPj0xhxT2LENYlh1yRG3JM4O+AvAsNDEwjWRgEsqShC5fIl+PE1F6CyvgQV9SUoXpKeKIiR995D/67dmDl9GkvvvBOVf7gPOUWJjTcQUfpFW1/A+OvoNgLoUr52ArAACC8ITVLK7mT/gPRJTHtnMT05C+/kDLyTs5hWPk+OTWN82IuJES8mhr0YH5nGxIgXo2em4J2YmfN78vJzsKSyCBX1JbjwimqU1RSjrGYRjNXFyC9IzZn/fGbcbgy0PYrhI0dQeOEq1O17CkVr1qT97xJRajH+OjojgNC5UZGmxZiVaZMRbymFGjo5ioP3/wezsxJyVsI360M8d7nyDbkoLilAUUkBSiuLsGx1GUqWGlBS7v9YUm6AYXG+KtM3pZQYPnoUA3taMTs6ioqtW1GxpRmCYXREusX46yQFioAQYr0QwhJ+tSCEaAbQDADLaxpw0ZU1ELkCubkCIkcgJ0cgrzAXBYY8FBj8n/OVz4WL8lBcUoC883CGn4zpgQH079qN0XffheHSS7F8zx4YLlqtdrOIiNJSEDwAApNljQBcoRuVsQNIKe3KtnNKnzLQfADwDypfu1H/B0wpJTydnThtexxyZgZVViuW/vY3ELnaLFxElH3SURBeBRAY2TYB6AYAIYRRSumBf1zBqWxvAHDu2u0M4/3uO/Tt2Inxjz5C8ZVXovaRFhQsX652s4iI5kh5QZBSOoQQa4UQFgCekBlGbwO4XNneLIRwA+iNMAMpY8jZWbhfOoTBp5+GyMtDTctuGJuaGDtBpFE2mw0ulwvr1q1DU1MTOjo64PF40NvbGzF3KJbu7m60t7ejq6sr9jdrQFrGECKtLZBSXj7f9kwz+fXX6Nu2HZOffYbF11+Pml0PI7+6Wu1mEdE8jEZjcLGZyWQKri04cCD+Q1bo9FOLxZJUIVELc5NTTHq9GHxmP47f0YjpkydxwRN7UffcsywGRDqwdm3kdVzR/j2cx+PRzdVAJCwIKTTx6ac4fscdGHr2WSy56SaYjv4bpTffzFtERDoRbT2A0+lEU1MT7HY7bDbbnHA6q9UKq9Ufs9HT04Oenp5g7EVAd3c3bDbbeXnq2UIwyygFfBMTGNz3NNyHDiGvqgp1LzyPkuuuU7tZRLrW39aGqS9Tm3ZaePGPUPPQQwn/XGNjI6xW65w46vb2dgCY81yDQJxF6Iri0JiLjo6OtD/1bCFYEBZo7L8f+cPoTpyA8dcbUfXAA8hdHPnRl0SkX+EH8nhXJKc7sjqVWBCSNDs8jNOP74WnsxP5K5Zj+aGXsOiKK9RuFlHGSOZMXg2BeOxwDocDZrM5qSgLtXAMIQkj77wD5y9ugef117H07rtgeuMNFgOiDOHxeGC32+F0OmGz2YKPvnQ4HHPGAAJR1h6PJ/i4SwDB6GyTyTTn5+x2OxwOR1xBdmpJS/x1Kp2v+Ot4zLhcGGhtxfDRt1C4ejVqW1tRdOlP1G4WUcbQevy1Vmk2/joTSSkxfOQIBlrbMDs2hor7fo+Ke+5hGB0RZRQWhBim+/r8YXTvv4+iNWtQ27oHhatWqd0sIqKUY0GIQvp88Lz2Gk4/vhfS50P1Q39C2ebNDKMjoozFghCB95tv0Ld9B8Z7elD806tQ29KCgvp6tZtFRJRWLAgh5MwM3C++iMFn9kMUFKC2dQ9Kb7+dK42JKCuwICgmv/rKH0b3+edYfMMNqNm5E/nVVWo3i4hUYLPZYDabg1NE0/0sY63I+nUIPq8Xp/ftw/HGJkz392PZU0+ibv8zLAZEWaqpqQmNjY2wWCxobGyE2+0+J5soU2X1FcL4J5+gb/sOeHt7UXrrraj6oxV5ZWVqN4uIVOJ0Os95oH1zczMaGhqy4iohK68QfGNj6G9rw7ebNsM3Po76Ax24oP0xFgOiLOdwOCI+0N7tdgdXJEdKObXb7XPSUPUq664QRj/8EP07H8b0qVMo27QJlfffj9zFi9RuFhGF+eC1rzF0YjSlv7OifjGu3ZD8M9otFkvElNNIaah6lDVXCLNnz+L7bdtw4u57IPLysOLlv6Fm5w4WAyIKMpvN54TVBQaWAyF10cLqtBxrHa+suEIY7upCf0sLZt1nUH7vvajY+jvkGAxqN4uI5rGQM/lkmUwmmEymYFIp4H98ZuCqIFS0lFM9y+iCMDM0hP49rRg5dgyFF1+M+hdeQNEll6jdLCLSsM7OTthstuC4gdFoDD5bGYiccup0OoOppnq+UsjItFMpJc6++SYGHn0McnwcFVu3ovzuuyDy89PUSiJKBaadJodpp1FMnzqFvl27MfbBByi67DJ/GF2EWQNERDRXxhQE6fPhzCuvYPCJP0MCqN62DWWbN0HkZM24ORHRgmREQZhyHkffjh2Y+PhjLLr6atS0tKCgbpnazSIi0hVdFwQ5PQ3XX1/E0P79EAYDatvaUHrbrxhGR6RjUkruwwlI5TiwbgvC5Bdf4Pvt2zH1xZcoufFG1OzYjrzKSrWbRUQLYDAY4HK5UF5ezqIQByklXC4XDCmaRq+7guCbmsLQc8/DdfAgcsvKsGzfPiz5+Y1qN4uIUqCurg4nT57E4OCg2k3RDYPBgLq6upT8rrQUBCFEIwAPALOU8pxgj1jboxl3ONC3bTu8x4+j9LbbUG19ELlRVg0Skf7k5+dj5cqVajcja6W8IAghzAAgpewWQpiEEGYppSPe7ZHMjo5h8MkncebwYeTX1qL+4EEsvuZnqW46EVFWS8eczI3wn/0DgBOAJcHtc/hGR+H85S04c/gwyjZvhulf/2QxICJKg3TcMjICcIe8Lk9w+xzeb75FzkoTVvz9ZRTreEk4EZHWaXJQWQjRDCAQHjK16thb/8Oxt9RsUjpVABhSuxFpksl9A9g/vcv0/l2U6A+koyB4ACxVvjYCcCW4HVLKAwAOAIAQoifRPA49yeT+ZXLfAPZP77Khf4n+TDrGEF4FEAgPMgHoBgAhhHG+7UREpK6UF4TAjCEhhAWAJ2QG0dsxthMRkYrSMoag3PIJ/7fL59s+j0S+V48yuX+Z3DeA/dM79i+M5p+HQERE5wezoSllAosOo2xrFEJYhBAPns82pVKM/rUrn5ujfQ+R1mmmIMQ6YOj9gBJH/3R9QFHGhP4SZVtwdToAz3wHVq2ar3+KZiFEL/yLLXVHCNGsfJz78GDoe/+Lo2+63/eUjwW/d5ooCLEOGHo/oMTZfl0fUJS+uaNsTmh1uhbF6B8ANEkpG5Tv0xWl2HUrY3sm5XXodt3uf7H6ptDtvqe8F+uV98a80GOnJgoCUhx3oUHxtF+3B5Q4JLQ6XafMej2Dhn/6d+D/pBM/TAsP0PP+F6tvgI73PSmlQ0ppVV6aIszaTOi908pK5ZTGXWhQPO03K/nvCSXAkjYE3jMhxHohhEVPB5ewWX9m+NcKhdLt/hdH34AM2PeUE5EtETYl9N5p5Qoh60kpbcpBpDzKZa2exVydrmfKPdpG5aULkc9CNU+5ndCViWuD5utbJux7SiHbErIAOClaKQgLjrvQuHnbnykHlHCZvjo9pH9O/NCnBgAJRwZohCXKGbLe9z8gSt/0vu8JIULHDZz4IQMuIKH3TisFIdPjLmL1T/cHFGWnWhuycwEZtDo9jv5tULb16rR/zSG3vSzK54zY/2L0Te/7ngVzD/hOIPn3TjML05QpX074B0YCwXYfB1Y4R9quJ3H2z61s1+V9TNIn5SDZCf//v6XwD7J2Z8L+l0DfdLnvKQf+DfC3f72Ucovy70m9d5opCEREpC6t3DIiIiKVsSAQEemQMqBsifY6GSwIREQqibCyOJGIkI3wDyRHe50wFgQiIhWE52MlERESvsgu0qK7hGhlpTIRUVZRZjuFriLeCKBL+ToQM+EIm+oM+LOZPEgDXiEQJUi5rD8jhDAKITp1ml9Eaabc0zeFvLbEWEkcMWZCSmkP+wgUg7UA1oX8zvDXCeMVAlGCpJR2JftmA4AOPeUW0fkjpXQosdvdUBaHpfLMPnxNQSrWh/AKgSgJUko7/GFielvZSueRcpC2wr8oLNaJg+oRISwIRElQBgTvBRDxoSREwJyV0gi9fRSF6hEhLAhECVKiAKxKZpEp2pOqKLuFzhpSrhTmjCGE52NpIfOL0RVERASAVwhERKRgQSAiIgAsCEREpGBBICIiACwIRESkYEEgIiIALAhERKRgQSAiIgAsCEREpPg/jWO93xZdQEYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def desingularise_our(hu, h, eps):\n", " return hu / np.maximum(np.minimum(h*h/(2.0*eps)+0.5*eps, eps), np.abs(h))\n", "\n", "def desingularise_kp1(hu, h, eps):\n", " return np.sqrt(2.0)*h*hu/np.sqrt(h**4+np.maximum(h**4, eps**4))\n", "\n", "def desingularise_kp2(hu, h, eps):\n", " return h*hu/(h**2+eps**2)\n", "\n", "def desingularise_kp3(hu, h, eps):\n", " return 2*h*hu/(h**2+np.maximum(h**2, eps**2))\n", "\n", "\n", "\n", "\n", "def desingularised_h_our(hu, h, eps):\n", " return np.maximum(np.minimum(h*h/(2.0*eps)+0.5*eps, eps), np.abs(h))\n", "\n", "def desingularised_h_kp1(hu, h, eps):\n", " return np.sqrt(h**4+np.maximum(h**4, eps**4)) / (np.sqrt(2.0)*h)\n", "\n", "def desingularised_h_kp2(hu, h, eps):\n", " return (h**2+eps**2) / h\n", "\n", "def desingularised_h_kp3(hu, h, eps):\n", " return (h**2+np.maximum(h**2, eps**2)) / (2*h)\n", "\n", "\n", "\n", "\n", "x = np.linspace(0, 3.0e-11, 500).astype(np.float32)\n", "h = x.copy()\n", "u = np.ones_like(h)*1.01\n", "hu = h*u\n", "u_true = hu/h\n", "u_our = desingularise_our(hu, h, 1.0e-11)\n", "u_kp1 = desingularise_kp1(hu, h, 1.0e-11)\n", "u_kp2 = desingularise_kp2(hu, h, 1.0e-11)\n", "u_kp3 = desingularise_kp3(hu, h, 1.0e-11)\n", "\n", "\n", "fig = plt.figure()\n", "#setBwStyles(fig.gca())\n", "plt.plot(x, u_kp1, label='$D_1(hu, h)$', markevery=50)\n", "plt.plot(x, u_kp2, label='$D_2(hu, h)$', markevery=50)\n", "plt.plot(x, u_kp3, label='$D_3(hu, h)$', markevery=50)\n", "plt.plot(x, u_true, label='Truth', markevery=50)\n", "plt.plot(x, u_our, label='Our', markevery=50)\n", "plt.ylim([0.6, 1.1])\n", "plt.xlim([0, 3.0e-11])\n", "plt.xlabel('x')\n", "plt.ylabel('Desingularized $hu/h$')\n", "plt.legend()\n", "plt.show()\n", "\n", "\n", "h = x.copy()\n", "u = np.ones_like(h)*1.01\n", "hu = h*u\n", "h_true = h\n", "h_our = desingularised_h_our(hu, h, 1.0e-11)\n", "h_kp1 = desingularised_h_kp1(hu, h, 1.0e-11)\n", "h_kp2 = desingularised_h_kp2(hu, h, 1.0e-11)\n", "h_kp3 = desingularised_h_kp3(hu, h, 1.0e-11)\n", "\n", "fig = plt.figure()\n", "#setBwStyles(fig.gca())\n", "plt.plot(x, h_kp1, label='$D_1(hu, h)$', markevery=50)\n", "plt.plot(x, h_kp2, label='$D_2(hu, h)$', markevery=50)\n", "plt.plot(x, h_kp3, label='$D_3(hu, h)$', markevery=50)\n", "plt.plot(x, h_true, label='Truth', markevery=50)\n", "plt.plot(x, h_our, label='Our', markevery=50)\n", "plt.ylim([0, 3.0e-11])\n", "plt.xlim([0, 3.0e-11])\n", "plt.xlabel('x')\n", "plt.ylabel('Desingularized $h$')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }