{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "The real and imaginary part of dielectric permittivity of the metals is important to simulate the optical properties of metal films and nanoparticles. Permittivity data is obtained experimentally by ellipsometry and is fitted with analytical models. The most common model for fitting experimental data is with Drude-Lorentz model shown below. \n", "\n", "$$\\epsilon(\\omega)=1-\\frac{f_1\\omega_p^2}{(\\omega^2+i\\Gamma_1\\omega)}+\\sum_{j=2}^{n}\\frac{f_j\\omega_p^2}{(\\omega_{o,j}^2-\\omega^2-i\\Gamma_j\\omega)}$$\n", "\n", "The first term is the Drude part. It represents the response of electron in the Fermi sea/conduction band when it sees external oscillating electric field (these transitions are called as intraband transitions). The Drude term has the plasma frequency ($\\omega_p$, oscillator strengh ($f_1$) and damping term ($\\Gamma_1$). The rest of the terms represent the Lorentz oscillators with specific resonance frequencies ($\\omega_{o,j}$), oscillator strengths ($f_{j}$) and damping terms($\\Gamma_j$) associated with them. They represent electron excitation from one band to another following an external oscillating electric field at certain resonance frequencies (these transitions are called as interband transitions).\n", "\n", "Just out of curiosity, I wanted to do this fitting by myself for a long time now. So I attempted to fit the data with the commonly used fitting method, the least squares method, which involves minimization of the square of the error between the model and experimental data, while changing the parameters. However I was not getting good fitting results even with several attempts, this may be due to the following problems with this method: \n", "- The fits were very sensitive to the initial guess values of the parameters (oscillator strengths, damping terms, etc). \n", "- The method is a local minimization technique. What that means is that the least square algorithm (Levenberg - Marquardt ) is changing the parameters in parameter space and searching for the minimum of sum of squares of difference between model and data. While doing this it can get trapped in a local minimum and report the parameters at that point as if it was as global minimum.\n", "For more information on the problems of least-square method. See the last section of the article here.\n", "\n", "To avoid these problems, researchers use global optimization techniques. See, Rakic et al. and Djurisic et al.. Rakic et al. used simulated annealing method , whic is a global optimization method and does a great job fitting experimental data. I wanted to use Scipy's simulated annealing method to fit the dielectric functions, but found that simulated annealing algorithm is deprecated in scipy 0.14 and they suggest using basin-hopping or differential evolution (DE) algorithms instead. \n", "\n", "I gave DE method a try because I didnot see any research article using DE algorithm for fitting dielectric data. If any of the readers find an article that uses DE to fit complex dielectric function, please point me to the article. DE algorithm is a stochastic method that does not involve finding a gradient but rather involves mutation and recombination of random population of solutions. I dont completely understand the guts of the algorithm, but you can find more information here. \n", "\n", "Below is my code of how I do the fitting of gold dielectric function with python and lmfit module. lmfit is a python module that provides a better and convinient fitting and optimization interface for Scipy optimize/minimze methods. The default minimization technique for lmfit is least square (Levenberg-marquardt) methood, so it needs to be changed to differential evolution. \n", "\n", "You can download my ipython notebook here \n", "\n", "You can use this method for fitting other dielectric data, but have to tweak the bounds a little bit. If you use my implementation for your fitting your experimental data, please let me know or even better cite it as:\n", "\n", "Juluri B.K. \"Fitting Complex Metal Dielectric Functions with Differential Evolution Method\". http://juluribk.com/?p=1597.\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following code, I get the gold's experimental dielectric data from Refractiveindex.info. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Keys in the data: ['REFERENCES', 'DATA']\n", "Reference : S. Babar and J. H. Weaver. Optical constants of Cu, Ag, and Au revisited, Appl. Opt. 54, 477-481 (2015)\n" ] } ], "source": [ "import urllib2 # To make a request and gather data,\n", "\n", "# Query the website\n", "\n", "response = urllib2.urlopen('http://refractiveindex.info/database/main/Au/Babar.yml')\n", "\n", "# Alternative Experimental data sources\n", "# response = urllib2.urlopen('http://refractiveindex.info/database/main/Au/McPeak.yml')\n", "# response = urllib2.urlopen('http://refractiveindex.info/database/main/Au/Lemarchand-3.96nm.yml')\n", "\n", "# Read the response. Response is in YAML format. Will use PyYAML parser,\n", "yaml_data = response.read()\n", "import yaml\n", "data = yaml.load(yaml_data)\n", "\n", "# Lets see what is inside the data,\n", "print\"Keys in the data:\" , data.keys() \n", "\n", "# Lets print the reference for this data,\n", "print \"Reference : \", data['REFERENCES']\n", "\n", "if 'COMMENTS' in data.keys():\n", " print \"Comments : \", data['COMMENTS']\n", "\n", "# Data seems to be stored in 'DATA', if needed uncomment the following line\n", "#print(data['DATA'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Data string contains real and imaginary part of the refractive index (not dielectric permitivitty). I clean up the real and imaginary part of the refractive data and convert them to dielectric permittivities" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Some clean up needed. The required data string holding wavelength, n and k, is in the form of dictionary item \n", "# with key 'data', which is inside a list\\n\",\n", "\n", "cleaned_nk = data['DATA'][0]['data'] \n", "\n", "# Use StringIO to convert String buffer as a file like class\\n\",\n", "import numpy as np\n", "from StringIO import StringIO\n", "\n", "# Read the data into as numpy vectors\\n\",\n", "wave_micron, n, k = np.genfromtxt(StringIO(cleaned_nk), unpack=True)\n", "\n", "wave_exp = wave_micron * 1000 # Convert wavelength from microns to nm\n", "eps_exp = (n + 1j* k)**2 # Convert n,k into dielectric function\n", "\n", "# Lets convert wavelengh in nm to eV and assign it to 'w' \n", "h = 4.135667516E-15 # plancks's constant in eV-sec\n", "c = 299792458E9 # speed of light in vacuum in nm/sec\n", "w_exp = h*c/wave_exp; # Convert wavleength in nm to eV" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I define a function that will plot experimental dielectric data. It will also plot fitted models if the data is provided." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "\n", "def plot_fit(w_exp,eps_exp, w_for_fit = None,eps_fit_result = None):\n", " '''\n", " This functions plots the experimental dielectric function as function of wavelength.\n", " If fitted model data is available it will also plot it\n", " '''\n", "\n", " %matplotlib inline\n", " import matplotlib.pyplot as plt\n", " from matplotlib.ticker import ScalarFormatter\n", " from matplotlib import rcParams\n", " \n", " rcParams.update({'font.size': 18}) # Increase the font size\n", " rcParams['mathtext.default']='regular' # make the mathtext the same size of normal text for better readbility\n", " \n", " # Plot the real part of the dielectric function\n", " fig,ax = plt.subplots(1,2,figsize = (12,5))\n", " ax[0].scatter(w_exp, abs(eps_exp.real), marker = 'o',facecolor = 'none', edgecolor = 'g')\n", " \n", " ax[0].set_xlabel('Energy (eV)')\n", " ax[0].set_ylabel(r'|$\\epsilon^\\prime$|')\n", " \n", "\n", " # Plot the imaginary part of the dielectric function\n", " ax[1].scatter(w_exp, eps_exp.imag, marker = 'o',facecolor = 'none', edgecolor = 'g')\n", " ax[1].set_ylabel(r'$\\epsilon ^ {\\prime \\prime}$')\n", " ax[1].set_xlabel('Energy (eV)')\n", " \n", " # If the fit data is available we will plot it with the actual plots\n", " if (w_for_fit is not None and eps_fit_result is not None):\n", " ax[0].plot(w_for_fit, abs(eps_fit_result.real),'-r')\n", " ax[1].plot(w_for_fit,eps_fit_result.imag,'-r')\n", " \n", " \n", " # Set the grid, log and axis formatter\n", " for axis in ax:\n", " axis.grid('on')\n", " axis.set_xscale('log')\n", " axis.set_yscale('log')\n", " axis.set_xlim([min(w_exp),max(w_exp)])\n", " for x_yaxis in [axis.xaxis, axis.yaxis]:\n", " x_yaxis.set_major_formatter(ScalarFormatter())\n", " \n", " fig.tight_layout()\n", " \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets plot the experimental data to see the data. Data seems to be from 0.1ev to 10eV. There seems to be three transitions (are all of them interband?) above 1eV" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAAFRCAYAAAC/qtYsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0VNXax/HvTiAEEiDUUERCkS4gHRQITQEVRYFrQ/Ha\nRaxXr+gVEi92sFwbAiIqFkRAxYoCYwEE6YIiClKkGyCkkrbfPzLJG2MSAinnzMzvsxZr5pyzZ8+T\nMMmTfXYz1lpEREREREQCXZDTAYiIiIiIiLiBGkciIiIiIiKocSQiIiIiIgKocSQiIiIiIgKocSQi\nIiIiIgKocSQiIiIiIgKocSQiIiIiIgK4vHFkjBlvjJlrjNlujMkyxvx+gvItjTEfGGMOG2MSjTHf\nGGP6FVI2yBhzlzFmizEmxRizyxgz2RhTpbzrFhERd1M+EhEJDMbNm8AaY7KAOGAt0AWIt9Y2LaRs\nM2AVkAY8CxwDbgDaAUOstYvzlX8OGAfMBz4D2niPvwUG2jzfmLKsW0RE3E/5SEQkMLi9cRRlrd3h\nfb4JqFJEMnoPGA50ttZu9J4LAzYDqdbaVnnKtgV+BOZZa0fmOX8b8D/gSmvtO+VRt4iIuJ/ykYhI\nYHD1sLqcRHQi3sQwDPDkJAvv65OAGUALY0zXPC+53Pv4bL6qpgPJwFXlUbeIiPgG5SMRkcDg6sbR\nSWgPhAArCri20vvYJc+5rkAm2UMTcllrjwMbvNfLo24REfEvykciIj7MXxpHDbyPewq4lnOuYb7y\nf1pr0wspX9sYU6Ec6hYREf+ifCQi4sP8pXGUs+rO8QKupeYrk/O8oLIFlS/LukVExL8oH4mI+DB/\nuWOU7H2sVMC10Hxlcp7XLqSuUMDmKV+WdQNgjHHvqhgiIqXIWmucjqGMKR+JiPiAwvKRv/Qc7fU+\nNizgWs65vMMQ9pI9nKBiIeX/tNZmlEPduay1fvdv4sSJfvfepVHvqdZxsq8rbvnilDtRmWuuucax\n/+uy+ufk57cs39/Jz3CAUD5y4T/lo9KtQ/mofP8pH5V+HUXxl8bRj2QPHehVwLUe3sfVec6tAoKB\n7nkLGmNCgY75ypZl3bliYmLweDwFXfJZ0dHRfvfepVHvqdZxsq8rbvnilHPy/9IpTn/N/vQZ9ng8\nxMTElPh9fYTP5yN/pHxUunUoH5Uvp79mf/wMF8np1vBJtPA2AduLuP4ekAG0z3MuHNgJbMlXth3Z\nK/i8n+/8OCALuKK86vZesyK+bOLEiU6HID7A+7vO8XxS0n/+no8mTpxoly5devL/wSIuoHwkRVm6\ndKmdOHFikfnI7ZvAjgYaew/HARWBp73HO6y1s/OUzdk1PB14Bkgge9fwtsD51tov89X9P+A2YAHZ\nu4a39r7Hd9ba/vnKllnd3vL2j/g/aFitoJESIu7n8Xgcv7Ml7meMwfronKNAykdu/rtA5ESUj6Q4\nispHbl+Q4Z9AX+/znN/WD3sfPUBuMrLWbjPGnA08DtxP9l4Qa4DB1tolBdR9J7ADuBE4HzhE9o7h\nE/IXLMu6c5xx6RnEjonl3ivvLayIiIhP8ng8/jBsOGDykYhIIHN1z1GgMMbYr3d8zai5o9h3zz6M\n8ckbqxLAdKdOisOXe44ChTHGTpw4kejoaP1Mi09SPpKi5Nysi42NLTQfqXHkAsYYO2HCBB7Z/QgJ\nryRQuWJlp0MSESk1xUlG4g4aVicigaCom3VqHLmAMcbOWDODKSum8NPYn5wOR0SkTKjnyP3UOBKR\nQFBUPvKXpbx93qPfPcrckXOdDkPklPjBfBIR8fLHrSUkcOizK0UpztYSbl+QIWBcGXclh346BHWd\njkREpHT5yYIMASOA9qQSkQCTM58yNja20DIaVucCGsYgIoFAw+rcT/lIRAKBhtWJiIhIsWhYnYj4\nq+IMq1PPkQvoTp34Oi2dKsWhniP3Uz4SX6d8JMWhniMfoDt1IuKvinOnTkRExA3Uc+QCulMnIoFA\nPUfup3wkIoFAPUciIiIiIiInoMaRiJSYhoSK+A8N8xZfps+uFEX7HPmQmJiY3LXXRUT8ifY58i2a\nHyYi/kr7HPkIjfEWkUCgOUfup3wkIoFAc47kb7bGbaXf6/1o+HRDzpt9HjuP7nQ6JBERERERR6lx\nFICS0pI4b/Z5XNLqElZct4LoxtEMfmswaZlpTocmPkpDpkRExA2Uj6Sk1DgKQBsPbKRW5VqM6z6O\n06ufzvje48nMyuS3w785HZqIiDhMCzKIiL8qzoIMmnPkAuU9xnvzwc0MeWsIW8dtJbRCKAnHE2jy\nXBPW3bSORtUblVscIhJYNOfI/TTnSEQCQVH5SKvVuUR5rlbXpk4b+jTuw4A3BnBes/P46JePGNV2\nlBpGIlImtFqdiIj4CvUcuYATd+qybBZvbniTrXFbaVe3HZe1uwxjdENXTo3H49Ey9HJC6jlyP/Uc\nia9TPpLiUM+R/E2QCeKajtc4HYaIiIiIiGuo58gFdKdORAKBeo7cT/lIRAKB9jkSERERERE5ATWO\nRKTENNlexH9oKW/xZfrsSlG0lLeP0DAG8XWaACvFoWF17qd8JL5O+UiKo6h8pMaRC/hTMjqaepT7\nvryPdfvX0axGMyafO5nTqp3mdFgi4gJqHLmfP+UjEZHCaM6RD/CHYQzWWi5+92KstTw/5Hla1W5F\nv9f7kZSW5HRoIuKg4gxjEBERcQP1HLmAv9yp2x2/m67Tu7Ln7j0EBwUD0PPVnjza/1H6NenncHRS\nljSMQYpDPUfu5y/5SAKX8pEUh3qOpFyEBIdwPPM4xzOPA9kbzSamJRISHOJwZCIiIiIiJ6aeIxfw\npzt1Vy+4mj+O/cGVZ17Jou2L2JewjyXXLKFCkPYbFgl06jlyP3/KRyIihdGCDC7nT8koMyuTF1a9\nkLsgwz297qFKxSpOhyUiLqDGkfv5Uz4SESmMGkcup2Qkvk5jvKU41DhyP+Uj8XXKR1IcATPnyBgz\n3hgz1xiz3RiTZYz5/QTlWxpjPjDGHDbGJBpjvjHGFLhygDEmyBhzlzFmizEmxRizyxgz2RhTYLfI\nydQtIiLiFv6weqqISEECbhNYY0wWEAesBboA8dbapoWUbQasAtKAZ4FjwA1AO2CItXZxvvLPAeOA\n+cBnQBvv8bfAwLy32k6hbt2pyyczK5O9CXupUbkG4SHhTocjIqVAPUfup3wkIoEgYIbVGWOirLU7\nvM83AVWKaBy9BwwHOltrN3rPhQGbgVRrbas8ZdsCPwLzrLUj85y/DfgfcKW19p1Tqdt7Tckoj61x\nW7nwnQtJOJ5AQloCk/pN4o4edzgdloiUkBpH7qd8JCKBIGCG1eU0jE7E21AZBnhyGi/e1ycBM4AW\nxpiueV5yuffx2XxVTQeSgatKULfkc8W8KxjXbRx779nL5ls3M2XFFFbsXuF0WFIEDcERERE3UD6S\nkvKrxtFJaA+EAAX9xb3S+9glz7muQCbZQ+VyWWuPAxu810+1bskjy2axbv86bu5yMwCnVz+doWcM\nZe2+tQ5HJiISGLrP6M6OozucDkNExBGB2jhq4H3cU8C1nHMN85X/01qbXkj52saYCnnKnkzdkkeQ\nCeL06qezeHv2tKzk9GSW715OVESUs4FJkbQykIj/uLT1pVz87sVoeJ34IuUjKalAbRzlrDB3vIBr\nqfnK5DwvqGxB5U+2bsnntYte46oFVzF49mDavdSObg27MfSMoU6HJSISEO7tdS/bjmwj/ni806GI\niJS7Cicu4peSvY+VCrgWmq9MzvPahdQVCtg85U+2bgDGjBlDVFQUABEREXTs2DH37kfO+NlAOWYH\nTG07lUrNK1E3rC5JW5P4+uuvXROfjv9+vH79eu68807XxKNjdxx7PB5mzZoFkPv7Tdxva9xWrLVa\nKVR8kkf7HEkJ+dVqdXkVtVqdMaYnsAyYZK2dkO/aIOALYKy19mXvuS+A/t760vOVXwY0t9ZGnkrd\n3vNaHagUWWsxRgtilSclIykOrVbnfsYYW29yPR7t/yjXnnWt0+GInDTlIymOgFmt7iT8SPawt14F\nXOvhfVyd59wqIBjonregMSYU6Jiv7MnWDWjTvdKwP3E/A98YSMikEBpMacCCnxc4HVLAUCKSoniK\nsemev/HlTck/ueITNYzEZykfSUkFZM+R9/p7wCVApzx7EYWTvRdRSr59jtqRvSrdAmvtiDznxwHP\nAVdZa98+lbq919RzVAr6v96frg26EtsvlvX71zPsnWEsuWYJ7eq2czo0ESGweo60KbmIiHsF0iaw\no4HG3sNxQEXgae/xDmvt7DxlcxJGOvAMkEB2wmgLnG+t/TJf3f8DbgMWkJ2MWnvf4ztrbf98ZU+2\nbiWjEsrMyiT0kVCSH0imYnBFAK7/6Hq6NujKTV1ucjg6/6dhDFIcAdY40qbkIg5QPpLiCKRhdf8E\nHvb+qw1Uz3P8z7wFrbXbgLOB74H7gafIbsQMzt948boT+BfZDZwXgFFkJ6IL8hc8hbo1rK6EgkwQ\nEaERbDq4CchuLG0+tJk6YXUcjkxEAnFYnTYlFxHxTX7Vc+SrdKeudLy76V3u/PxOhrcazsaDG6lW\nqRoLL19IhaBAXZRRxF0Cqecor3JeIOgMa23dU6nbe175SET8XlH5SH81ukRMTAzR0dHqCi6By9pd\nRstaLVm2exn9m/RneOvhahiJuIDH41HPeOFKe1PynsaYCtbajFOoW0Qk4PnbsDqfldM4kpI5q/5Z\n3NbtNka2HXnChtH2I9vpO6svNZ6oQdfpXdmwf0M5Rel/9IevFCU6OjrghtWdBG1KLlKKlI+kpHRb\nXQJSemY6Q98aynVnXce8UfP4ZOsnDH17KJtv3UxEaITT4YlI4NCm5DrWcSlvSu6meHTsjmPPSWxK\nrjlHLmCMsRMnTiRaw+rKzda4rQyePZjtd2zPPXfOzHOY1H8S0VHRzgUm4oc83mF1sbGxmnP092va\nlFxEpJwF0mp1PkvD6spXRGgEh1MOE5ccB0BKegq7j+1Wr5FIGYjWsLqiuG5TchGRQKbGkQSkumF1\nuaXLLfR+rTf3f3U/fWf1JToqmg6RHZwOzSfldF2LyMmx1iYCC4FoY0z7nPPejcOvB7Zaa3/I85I5\nZA+duzNfVTcAlYG3SlA3oK0lxLfpsytF8RRjawkNq3MBDWNwhrWWj7d+zIYDG2heszmj2o4iyBR8\nv+CL377g018/JSI0grHdxlI3rG45R+tuHm26J8UQSEt5a1NyEWcoH0lxFJWP1DhyAc05crdZ62cx\nYekEbu9+O9sOb+OLbV+w6oZV1K5S2JxoEckrEOccGWOWAn29hzmJNudr9xTQiGkFPO59TQiwBoix\n1i4poO4gsnuObgSigENk9yhNsNb+bYGFk6xbjSMR8XtqHLmckpG7NftfM9699F26NszeSP7qBVfT\nqX4n7uyRf1SLiBQlkHqOfJVu1omIPyvOzTrNORI5geT05L8Mo4sMiyQ5/W83ZwOaxniL+A8tECS+\nTPlIilKcBYLUOBI5gZFtRnLjxzeyYf8G5v88n1kbZnFBiwucDktERERESpkaRy6h1YHca/K5k2lf\ntz1XzL+CycsnM2fEHNpHtv9LmeMZxxn7yVjqTa5Hs/81440NbzgUrTN0l1mKUpzVgcQ9lI/Elykf\nSVG0Wp2P0Jwj33fX53fxS9wvvHLBK+xL3Mel713KrItmMaDpAKdDE3ENzTlyP+UjEQkE2gRWpIx9\n8usnPDXoKRpVb0S3ht0Y23Usn/32mdNhlRvdZRYRETdQPpKSUuNIpBREhEaw7ci23OPfDv9GRGiE\ngxGJiJwaDasTEX+lYXU+QsMYfN+X277kivlXMKbDGPYm7mXlHytZef1KalWpBWRvODt5+WTe3Pgm\nFYMrck/Pe7jizCscjlqkfGlYnfspH4lIINA+Ry6nZOQfNuzfwCe/fkJYxTBGdxhNzco1c6899/1z\nzNowi1cueIXEtESu+eAaXrngFYaeMdTBiEXKlxpH7qd8JCKBoKh8VKG8g5GC5ewroVVWfFeHeh3o\nUK9Dgdfm/jSXKedOoVvDbgA82PtB5v00z28aRx6PR59dKVTOpnsiImVN+UhKSnOOXEKb7vm3KhWr\ncDDpYO7xgcQDVKlY5S9lMrIyOHb8WHmHJlLmirPpngSurXFb6fNaH+o8VYc+r/Vha9xWp0MSkQCm\nniORcjD+nPH84/1/8Nvh30g4nsCsDbP49tpvc69PXT2VexbdA0DbOm2Z/4/5nFbtNKfCPWlq2Iv4\nj/IcyZCSnsLg2YO5s8edjGo7irmb53Le7PPYfOvmv91AEikO5SMpSnFGMmjOkQtojHdg+GHPD8zZ\nPIeQ4BCuO+s6mtVsBsDy3csZNXcUX4/5mqY1mhLjieHbXd+y5JolDkcsUro058j9yjsfrdm7hms/\nvJaNt2zMPddhagdmDptJ5wadyy0OEQksmnMk4gJdG3ala8Oufzv//R/fc2nrS3MbS/eefS9PLn+y\nvMMrEY3xFpFTEREawYGkAySmJRIeEk5iWiL7E/dTPbS606GJj1I+kpJS40jEYQ2rNmTuT3NJz0yn\nYnBFVuxeQcOqDXOvr/xjJbM3ziY4KJgbOt1A27ptHYxWRKT0NKvZjOGthtN3Vl+GNh/Kp799yvBW\nw2les7nToYlIgNKwOhfQsLrAlpmVyYi5I9h+ZDsta7Vk6Y6lvHPpOwxsOpClvy/lH+//g3/1+hep\nGak8v+p5Fl+9mPaR7Z0OW+SkaVid+zmRj6y1zP1pLj8d+ok2ddowss1IjNHHRETKjvY5cjk1jiTL\nZvHV9q+IS46jZ6OeREVEAXD+2+dzebvLuar9VQA8uexJfo37lenDpjsYrcipUePI/YwxduLEidpa\nQkT8Us6CDLGxsYXmIy3l7RIxMTHaBySABZkgzm12LpefeXluwwggNSP1L5vJ1qpci9TMVACS0pJ4\n/6f3mb1xNvsT95d3yH+hz64UxePxaClvH6KtJcSXKR9JUYqztYR6jlxAPUdSmBlrZzBlxRSmnj+V\n1IxUrl94PVPPn8rZp59N79d6ExkWSY3KNfhu13d8Nforx+YjaQKsFId6jtxP+Uh8nfKRFIeG1bmc\nkpEUxlrL1NVTeW39awQHBXNH9zu4rN1lTFg6gT+O/cGrw17FGMNLP7zEJ79+widXfOJ0yCKFUuPI\n/ZSPRCQQaClvER9ljOGWrrdwS9db/nJ+X8I+ujXsljtpuWuDrsxYOwOAXfG7+Hbnt0SERnBe8/Oo\nEKQfcxEREZHi0JwjER/Uu3Fvpq6eysGkg6RmpPLk8ic55/Rz+Hbnt3Se1pmPtn5EzNcxDH1rKOmZ\n6WUej8Z4i4iIGygfSUmpcSTig0a3H80FLS7g9GdOp/rj1TEYHh/4OGM/HcuMC2cwZ8Qcvr/ue7Js\nFm9ufNPpcEXEh2iBIBHxV8VZIChg5xwZY8YDnYDOQBSw01rbpIjyLYEngD5ACLAWmGitXVpA2SDg\nDuAmoDFwCHgPmGCtTS6gvMZ4yynJzMok02YSEhwCQO0na7P51s1EhkcCMP6r8VSpWIUbO9/IO5ve\nIT0znYtbXcwZtc5wMmwJUJpz5H7KRyISCIrKR4Hcc/QIEA38ChwBCs0GxphmwHKgO9kNpHuBcOAL\nY8yAAl7yDDAF2ATcBswFbgcWGu1sJ6UoOCg4t2EE0KtRL55c9iRZNoudR3fy7uZ3OaPWGXSZ3oX1\n+9ezK34XvWb2YvXe1Q5GLSIiIuJOgdxzFGWt3eF9vgmoYq1tWkjZ94DhQGdr7UbvuTBgM5BqrW2V\np2xb4EdgnrV2ZJ7ztwH/A6601r6Tr37dqZNScTDpIKPmjmLlnpUYDI8NeIyd8TupEFSBJwc9CcDM\ndTN5/6f3+fTKT0vtfbV0qhSHeo7cT/lIfJ3ykRSHeo4KkNMwOhFvI2gY4MlpGHlfnwTMAFoYY7rm\necnl3sdn81U1HUgGrjrVmEVOpG5YXTxjPPx5758kjE/gjh53cCT1CGfU/P9hdGfUPIOjqUd5ff3r\ndJveja7TuzJ9zXQHoxYRERFxB63xe2LtyZ5jtKKAayu9j12AH7zPuwKZwKq8Ba21x40xG7zXRcpU\nWEhY7vMhzYfw0NKH6NmoJ1VDqvKfpf/h9OqnE/N1DDMunEFwUDA3LryRShUqcXWHq0/p/XSXTkRE\n3ED5SEpKjaMTa+B93FPAtZxzDfOV/9NaW9D6yXuAnsaYCtbajFKMUaRQo9qOYl/CPoa8NYS0zDRG\ntx/N9iPbmdRvEgOaZk+Ze3zg47y67lW6NOjCD3t+oH7V+gxqOghNkRMREZFAErDD6k5CFe/j8QKu\npeYrk/O8oLKFlRcpc3f0uIPdd+3mwL8OMPncyYSFhBGXEpd7PS45jrjkOKJnRfPl9i+5Z9E9XP3B\n1RR37oGW/RURETdQPpKSUs/RieUsvV2pgGuh+crkPK9dSF2hZK+K97flvMeMGUNUVBQAERERdOzY\nMbdrOOcHXcc6Lq3j3lm9efCbBzmaepRdG3Yx96e5ZDXO4ttrv+XolqOkVU/j7l/u5qvtX1FhV/av\niX79+hVa3/r161319enYHccej4dZs2YB5P5+Eykv8anx/GvRv1i9bzWNqzfm6fOepmmNAtddEhHJ\nFbCr1eVV1Gp1xpiewDJgkrV2Qr5rg4AvgLHW2pe9574A+nvrS89XfhnQ3Fobme+8VgeScrfxwEZe\nW/caFsvItiMZ8PoAUh5MyR1Kd+W8K0lOT2bJjiVkZmVybcdreWbwM1QI0j0VOTVarc79/CUfWWsZ\n9OYgGldvzM1dbmbJ70t4efXLbLh5A9VDqzsdnog4TKvVlcyPZA+T61XAtR7ex7ybxqwCgsneEymX\nMSYU6JivbC7tSC7lrX1ke54Z/AzPDn6WsxudTes6rZmyYgrWWtbvX8/CrQvZEreFTbds4vc7fmfz\noc089u1jToctPshTjB3JxT38IR8dSj7E2n1rmXbhNLo27Mq/z/k3TWs0Zfnu5U6HJiIOKk4+Us8R\nxd7n6BKgU559jsLJ3ucoJd8+R+2ADcACa+2IPOfHAc8BV1lr385Xv1/cqRPftv3Idka8N4KfDv1E\n5YqVaVOnDTd3vpnRHUYDsHj7Yh5c8iCd6nfiSOoRhjQfwuj2ozHG4NG+ElIM6jlyP3/JR0dTj3La\n06ex9569VKtUjSybRdfpXXlq0FP0b9Lf6fCkDCkfSXEUlY8CdnyMMWY00Nh7WAeoaIz5j/d4h7V2\ndp7i44EBwCJjzDNAAnADUB84P2+91tpNxpgXgduMMfOAz4DWwDiy90r6S8NIxC2a1mjK2pvWkpye\nTOUKlbnj8zvYdHBT7vVlu5ex4cAGBjUdRM/TevLEsifYn7if+86+z8GoRUT+LiI0gtHtRzN49mCu\n7nA1S3csJbRCKOecfo7ToYmIywVsz5ExZinQ13uY803IaUF6rLX985VvBTzufU0IsAaIsdYuKaDu\nIOBO4EYgCjgEzAEmWGv/thiDMcZOnDiR6Oho3e0Q19hzbA+9Zvaix2k9qFyhMvN/ns/ApgOZ/4/5\nAGz5cwv9Xu/HxL4TSUxLZOgZQ2lTp43DUYsbeTwePB4PsbGx6jlyOX/pOQLIsllMWzONNXvX0Dii\nMXf1uOsve8CJSOAqqucoYBtHbuJPyUj8S1xyHPN/nk9GVgb7EvbxZ8qfvHT+SwBs2L+BztM6M7z1\ncOqH1+edTe/w3oj36Nekn8NRi1tpWJ37KR+JSCDQsDofEBMTo54jcZ1aVWpxQ+cbAPj9yO90n9Gd\ntnXa0qJWC65feD2ta7dm7si5eDwe+l3Qj/sX38/K61c6HLW4TU7PkYhIWdOcIykp9Ry5gO7Uia/Y\nsH8DsV/HciT1CJlZmQxuPpgHej+Ax+Mhsm0kw94dxq/jfnU6THEp9Ry5n/KR+Do1jqQ4NKzO5ZSM\nxBct2raIGxbewIeXfUiDqg24ceGNnFbtNF4Y+kJumeW7l7Plzy20rt2ano16OhituIEaR+6nfCQi\ngUD7HPkAf9hXQgLLuc3OZUKfCVz4zoW0fKEltSrXYvK5k3OvP/z1w1wx7wq+3vk1l827jEe/fdTB\naMVJ2udIRER8hXqOXEB36sTX5R/GsPPoTjpP68zPY3+mTlgdDiQeoPWLrfnxlh9pWK2hc4GKo9Rz\nVDBjzHigE9CZ7BVOd1prmxRRviXwBNCH7NVT1wITrbVLCygbBNwB3ET29hWHgPcoYvVU5SPxZRpW\nJ8WhniMRKVcHkg7QOKIxdcLqABAZHkmj6o04kHTA4chEXOkRIBr4FTjC/28v8TfGmGbAcqA72Q2k\ne4Fw4AtjzIACXvIMMAXYBNwGzAVuBxYaY9RQFRHJR6vViUiJ5b9L17JWS/Yc28OHWz5kWMthLNiy\ngINJB2lRqwUAx44f462Nb5GQlsDg5oNpH9negahFXKOptXYHgDFmE1CliLKPAdWAztbajd7XvAFs\nBl4EWuUUNMa0JXsD8nnW2pF5zv8O/A+4DHinVL8SEYep10hKSj1HLqE5R+JPqodW54PLPuCOz++g\n0qRK/GvRv/jwsg8JDwnnaOpRer7ak8W/L2Z/4n4GvjGQz3/73OmQpQxpzlHRchpGJ2KMCQOGkb1R\n+cY8r08CZgAtjDFd87zkcu/js/mqmg4kA1edaswiIv5Kc45cQGO8xdcVNcY7JT2FyhUr5x5PXj6Z\ndfvX8dYlbwHw2a+fMX7xeNbfvL48QhUHac7RieX0HFlrmxZwrSewDJhkrZ2Q79og4AtgrLX2Ze+5\nL4D+3vrS85VfBpxhra2b77zyUT4ZWRn8cewPIkIjiAiNcDocOQHNOZLi0JwjEXFM3oYRwJGUI7So\n2SL3uEWtFhxJPVLeYYn4ogbexz0FXMs5l3fFkwbAn/kbRnnK1zbGaHh9EbYd3kbbl9rS+7XeNHqm\nEY99+5jTIYlIGVPjSERK7GTu0p3b7Fymr53OD3t+4EDiAe776j4GNxsMwI8HfmToW0Pp9Eon7vz8\nTlLSU8ooYhGflDMX6XgB11Lzlcl5XlDZwspLPlctuIpbutzC7rt2s/W2rUxfO50lvy9xOiwpgnqN\npKR0x8iMEiC2AAAgAElEQVQlYmJiiI6O1g+1+L2+UX15fODjjJw7koS0BIa1HMYzg59hb8JeBr05\niJjoGDrX78zjyx7nnx/9k3cu1XxxX+fxeDSnsnTkLL1dqYBrofnK5DyvXUhdoWSvive35bzHjBlD\nVFQUABEREXTs2DE3N+X8PwbK8Zrla5hwevYIxvpV63NW6lnM/WQu/W/r74r4dKxjHRfv2OPxMGvW\nLIDc32+F0ZwjF9AYb/F1nlIY4z1r/Sw+/+1z3h3xLgDJ6cnUeKIGieMTqRhcsRSiFKdpztGJlfOc\no+bW2sh855WP8mjzYhsm9Z/EJa0vISU9hV4zezGhzwSGtx7udGhSiNLIR+L/NOdIRFwvtELoX+Ye\nxafGE2SCCA4KJjUjlaW/L2XJ70tIzUgtohYRv/Yj2cPkehVwrYf3cXWec6uAYLL3RMpljAkFOuYr\nKwWYedFMbvnkFga8MYA2L7WhY72OXNzqYqfDEpEypJ4jF9CdOhFITEukx4wedG/YnS4NuvDy6pcZ\n3mo4t3e/nf5v9KdScCWMMaSkp7DkmiXUrlLYaCFxK/UcnVhRPUfe6+8BlwCd8uxzFE72Pkcp1tq8\n+xy1AzYAC6y1I/KcHwc8B1xlrX07X/124sSJRGuYd65DSYdYt38dtSrXolP9TmjvXBHf5fEO846N\njS00H6lx5AJqHIlkO5JyhGe+f4b9ifvp07gPV555Jbd/djsWy/NDngfgjs/vIDMrkxfPf9HhaOVk\nqXFUMGPMaKCx93AcUBF42nu8w1o7O0/ZZmT3CKUDzwAJwA1AW+B8a+2X+er+H3AbsAD4DGjtfY/v\nrLX9C4hF+UhE/F5R+UgLMohIiZXWGO8alWvwcL+H/3Lu96O/c0OnG3Lv1g5qOoiXV79c4vcScZF/\nAn29z3NaJjk/CB4gt3Fkrd1mjDkbeBy4HwgB1gCDrbUFLaN2J7ADuBE4HzgE/A+YUEBZQAsEiW/T\nnCMpSk7PUVHUc+QCGsYgvq4sk9GEpRPYeGAjc0bMwRjDP97/B23rtOWilhdx75f3ciDpAP2i+vHU\noKcICwkrkxikZIozjEHcQT1H4uvUOJLiKKrnSI0jF1AyEilcakYqV8y7As8OD+BdCnzA4/R+rTeT\nz53MWfXOYtK3kwgyQVr22+U0rM79lI9EJBCoceRySkYiRbPWcjDpIBZLZFgkM9bO4Lvd3/H6xa8D\nkJSWRM0na5L8QDLBQcEORyuFUePI/ZSPRCQQaClvESlTJxq/W1LGGCLDI6kXXg9jDFUqVuFA4gFy\n/og7lHyISsGVSMtM4/X1r/P0iqdZu29tmcYk4q9iYmLK/GdapKzosytF8Xg8xMTEFFlGPUcuoDt1\n4uvKe4x3UloSPV/tScd6HelYryOvrHmFMR3G8PGvHxNWMYzWtVvz7uZ3eXHoi4xoM+LEFUq5UM+R\n+ykflZ5DSYdYu28ttavU1hLg5UhzjqQ4NKzO5ZSMRE5efGo8L6x6gQNJB4iOiiYxLZE3N77JoqsW\nYYzh+z++Z9TcUey6a5fToYqXGkfup3xUOlbsXsHFcy6mXd12bD+ynf5R/ZkxbIYaSCIuoaW8RcTv\nVA+tzoN9Hsw9fmbFM7Su3Tr3j482ddoQlxIHwPGM4wQHBVMhSL/yRE5ES3mX3LUfXssrF7zCxa0u\nJjk9mV6v9uKjXz7iolYXOR2aSEDTUt4+QnfqxNe5YRjDun3rGPzWYBb8YwFt67Tl31/9mz0Je6hS\nsQofbPkAg+GO7nfw+MDHdffWIeo5cj/lo9IR8t8Q4u+Pp3LFygCM+3QcTWs05a6edzkcmf9zQz4S\n99OCDD5AE2BFSuas+mfx8vkvc/m8y2n4dEMOJB3gtKqnkZmVSfz98fxx9x8s/n0xM9fNdDrUgFOc\nCbAi/uSs+mcxfe10APYl7OPjXz/mrPpnORyViBSHeo5cQHfqRMpGt+ndeG7wc/Rs1BOA6Wums/yP\n5bx20WsORxaY1HPkfspHpePXuF+54J0LSE5P5mjqUR7s/SD3n3O/02GJiFepzTkyxiwFSvJb0wDW\nWtu/BHWIiBRL/ar1WbVnFT0b9cRay6o9q4gMi2TW+lkcSDxA78a96dWol9NhiriK5hyV3Bm1zmDz\nrZv549gf1AitQfXQ6k6HJCKUwZwjY4yH7MZRSe78WWttvxK83u/oTp34OreO8d7y5xb6v96fno16\nEp8az77EfdQIrUFohVA61uvIO5veYVK/SVx71rVOhxoQ1HPkfspH4uvcmo/EXUqt58haG10qEYmI\nlINWtVux/ub1fLX9K0KCQ0jNSGXq6qksGr2IIBPEP8/6J2fPPJsxHcdokQYRERHRnCM30J06kfIx\nbc00vv/je2ZelL0oQ3pmOlUerULKgyla5rscqOfI/ZSPRCQQaLW6UmCMGW+MmWuM2W6MyTLG/H6C\n8i2NMR8YYw4bYxKNMd8YYzScUMRBfRv3ZeHWhXy1/SsOpxzmnkX30PO0nlz07kW0eL4Fl753KfsS\n9jkdpoiIiDhEjaPiewSIBn4FjlDEwhTGmGbAcqA78ARwLxAOfGGMGVDmkYqUM19Zhr5l7Za8OfxN\nbv3kVpo814Rth7exK34XfRv35cPLPqRlrZYMeWsIGVkZTocqIgFm/f71zFo/i292fuN0KD7NV/KR\nuJcaR8XX1Fpbx1p7HnCiW8uPAdWA86y1T1hrXwZ6A3uBF8s4ThEpwuDmg9k6bivx98fzQO8HqBNW\nh/vOvo/WdVrzSP9HiD8ez7bD25wOU4pgjHnBGNPb6Tj8lfbdK3+vrH6FoW8NZfHvi7nuo+u463Nt\nFitSFoqz757mHJ0CY8wmoIq1tmkB18KAOOBba+2gfNf+AzwMdLfW/pDnvMZ4izhg3b51jJg7gp/H\n/kxIcAiJaYlEPRvFmhvXULVSVWqE1tBCDaWotOYcGWNaAcOAKsA8a+2PJQ5OAOUjJySmJVJ/Sn02\n3LyBpjWacuz4Mdq91I6PLv+IjvU6Oh2eiF8qtdXqpFjaAyHAigKurfQ+dgF+KOC6iJSjDvU60D6y\nPee/fT6Dmw1m3s/zOKveWXR8pSOZWZnUCavDgn8soH1ke6dDlTystVuALU7HIVIa4pLjqF6pOk1r\nZN9vrVapGq1qt2Jfwj41jkQcoGF1pa+B93FPAddyzjUsp1hEyoWvDsEJMkHMHTmXi1tezO5ju7m4\n1cWs27+ORVct4tj4Y8T0jeGidy8iMyvT6VBFxE81rNaQisEVeX3961hr+Xbnt6zbv44O9To4HZpP\n8tV8JO5RZONI47pPSRXv4/ECrqXmKyMiDqsQVIGx3cby7OBnaVGrBT1O60HXhl0BGN1hNMnpyRxI\nOuBwlJKXcpP4kwpBFfjoso949LtHqfxIZUbMHcFbl7xFg6oNTvxiESl1JxpW9wIwzLvCmsZ1F0+y\n97FSAddC85XJNWbMGKKiogCIiIigY8eOuTs859wF0bGO3Xycwy3xnMpxo2qN+GH5D3xc+2MuOPcC\nfj70Mwm/JLBp1SYaDGjgeHy+duzxeJg1axZA7u+3UqLcJH7lzMgz+eW2X0hKS6JKxSqa61gCOb+L\nRE5VqS3IYIypb63dl+e4qbV2e6lU7jInWJChJ7AMmGStnZDv2iDgC2CsdwW7nPOaACviEvd8cQ8L\ntiygc4POfLPzG/7R9h98s/MbktOTGdFmBA/3e1gbxp4iJzaBDaTcVBqMMXbixIlER0frj0wR8Tse\njwePx0NsbGyh+ahUGkfGmLuByUAra+1W77kuwJXAfdba9BK/iYucoHEUDhwClllrB+a79hAQSwGr\n1SkZiS/zeDx+9dldvns5u+J3kZaRxv2L7+fN4W8SGR7J2E/H0vv03kzqP8npEH1KcZJRWQi03FQa\ndLNOfJ2/5SMpG0XdrCutxtEC4GPgdeAqa+0s7/mhQH9r7b9K/CYuUlTjyHv9PeASoJO1dqP3XDiw\nGUix1rbKV17JSHyavyaje764hzphdbj/nPuB7E0ar5x/JZtv3exwZL6pvHuOAi03lQblI/F1/pqP\npHQVlY9Ka7W6DdbaV621GcDwnJPW2k+BC0rpPRxljBltjPmPd6+iOkBEzrEx5qp8xccD8cAiY8y/\njTG3At8C9YFxBdWvTffEl/lrIgoLCWNvwt7c4z3H9hAeEu5gRL7JU4xN98qI3+cmCSzWWl7+4WUG\nzx7MqLmjWLdvndMhuY6/5iMpP6XVc/SotfYB7/NF1tpzvc8rAPustXVK/CYOM8YsBfp6D3O+aTkt\nTo+1tn++8q2Ax72vCQHWADHW2iUF1K07dSIutDdhL91ndOfCFhdSL7weL/7wIk8MfIIaoTWoG1aX\nHqf10MTpk+BAz5Hf56bSpnzkbk989wRvb3qbSf0msSt+FzFfx/Ddtd/RsnZLp0MT8SnlMazuXbI3\n5KsB3Ag8BbwL3Ab0sda2K/Gb+DElI/F1/jyMYV/CPmasnUFyejINqzUk9utYujfszi9xv3DO6ecw\nc9hMNZCKyYHGkXLTSVI+crfm/2vOvFHzcvdA+teif1E1pCoToyc6HJl7+HM+ktJTVD4qrSWXJgFL\nye5JuQqYDfzHe+3WUnoPvxYTE6MFGURcqH7V+jzU9yEAGj/bmLcveZtBzQaRkp5Cj1d78Omvn3J+\ni/MdjtLdchZkcIByk/gVYwyZ9v83pc7MytTNGZFSVppLeVcEKlprk40xbYChwA/W2q9L5Q38mO7U\nibhfls2i4n8rcvw/x3OX8r5p4U10qNeBW7vq7+zicGgpb+Wmk6B85G7Pff8cU9dMZUKfCeyK38WU\nFVP4/vrvaVrjr+tDLfxlIXM2z6Fyhcrc0eMO2tVVJ6lIXuWxIAPW2nRrbbL3+U/W2slKPiLiL4JM\nEGfVO4sXVr0AwI6jO/jk10/oXL+zw5FJUZSbxJ/c3v12xp8znvd/fp8tcVv4eszXf2sYvbXxLcZ+\nOpb+TfrTvGZz+r3ej58P/exQxCK+p9R6juTUaZ8j8XWBMsb7t8O/ceE7FxKXHEdyejJPDHyCsd3G\nOh2W6zm1z5GcPPUc+b5u07vx6IBHGdg0e6vFCUsnkJyezORzJzscWfkIlHwkJVNqc468K7aV5Lem\nAWz+ld0Ep5a5FZGT0LxmczbfupkDiQeICI3ggy0f0PaltqSkp3BZu8t4uN/DuUPu5P/l3PiJjY11\nOhQpBs2B9W2ZNpPQCqG5x6EVQkk4nuBgRCLuUZw5sCfVc2SM8ZDdOCrJnT9rre1Xgtf7Hd2pE/E9\nX277kms/vJa3L32bOlXqcNPHN9Evqh+x/dQAKIwTc47k5Cgf+b7nVz7Py6tfZsq5UziUfIh7Ft3D\nZ1d+RpcGXZwOTcQ1ynwpbykZJSMR33P7Z7cTFRHF3T3vBmD13tVc/9H1rL95vcORuZcaR+6nfOT7\nrLVMXzs9d0GGe3vdS9+ovid+oUgAKZcFGaRkYmJinFrqVqTEAvGzGx4Szu743bnHu+N3Ex4S7mBE\n7uXxeDR0WKScGGO4sfONLL56MR9f8XHANYwCMR9J6VLPkQvoTp34ukCcALvn2B56vNqDC864gDph\ndZi6eiqzL5nNuc3OdTo011LPkfspHwWGY8eP8fWOrwkOCqZfVD8qV6zsdEilJhDzkZw8DatzOSUj\nEd+0N2Evr617jZSMFIa3Gk7nBlrWuyhqHLmf8pH/2x2/m+jXo2kS0YTUjFQS0hLwXOOhRuUaTocm\nUm7UOHI5JSMRCQRqHLmf8pH/G71gNE0jmhLbLxZrLTd/fDPVKlXjqXOfcjo0kXKjOUc+QHOOxJfp\nsytF0ZwjEffYFb8rdx6SMYa+UX3ZdWyXw1GVHuUjKSk1jlwiZ18JEfFdW+O28tX2r9ibsNfpUFwl\nOjpajSMRl+jWoBsv/fASxzOOk5iWyIy1M+jWoJvTYYm4hobVuYCGMYj4vke/fZTnVj5Hmzpt+PHA\nj7x20Wtc2PJCp8NyFQ2rKx3GmPFAJ6AzEAXstNY2KaJ8S+AJoA8QAqwFJlprlxZQVvnIz6Wkp3Dl\n/CtZtG0RWTaLy9tdzrQLpxEcFOx0aCLlRnOOXE7JSMS3bTywkcGzB7PupnVEhkeyas8qBs8ezP5/\n7SckOMTp8FxDjaPSYYzJAuLIbuR0AeKttU0LKdsMWAWkAc8Cx4AbgHbAEGvt4nzllY8CxNHUowSb\nYKpWqup0KCLlTnOORKRMBfoY7+1HttO5QWciwyMB6NawGxWDK/Jn8p8ORyZ+qqm1to619jxg3wnK\nPgZUA86z1j5hrX0Z6A3sBV4s4zjFxSJCIwpsGCWmJTJl+RT+/eW/+ezXzxyIrGQCPR9JyalxJCJS\nQm3qtGHlHyvZGrcVgI+3fkywCaZuWF2HIxN/ZK3dUZxyxpgwYBjgsdZuzPP6JGAG0MIY07VMghSf\nlJyeTJ/X+rByz0qqh1Zn7KdjeX7l806HJVKuKjgdgGTLWZBBizKILwr0z22LWi14ctCTdJvejVpV\napGcnsy8UfOoEKRfsZB9J1d3cx3Rnuw5RisKuLbS+9gF+KHcIhJX+2DLB9SuUps5I+ZgjGFkm5F0\nnd6V27rdhjG+MSI20PORlJzmHLmAxniL+If41HgOJB2gUbVGfrXjfGnRnKPSZ4zZBFQpaM6RMeZS\nYC5wi7X2lXzX2gCbgEettf/Jc175KIBNWzON7//4npkXzQSye5JqPFGDlAdTCDIabCT+Q3OORKRM\nqVcgW/XQ6rSo1UINI3GLKt7H4wVcS81XRoQBTQawcOtC3v/pfX47/Bs3LryRYS2H+VTDSPlISkpj\nPkRERPxTsvexUgHXQvOVyTVmzBiioqIAiIiIoGPHjrlDlXL+8NSxfx7v3ribCadP4KnlT3Ew6SCt\nEloxrts4cjgdX3GO169f76p4dOyOY4/Hw6xZswByf78VRsPqXEDDGEQkEGhYXek7wbC6nsAyYJK1\ndkK+a4OAL4Cx3hXscs4rH0mhktOTWbZrGcYYzjn9HEIrhJ74RSIuVFQ+Us+RiEgZ+HrH12w6uIkz\nap3BoKaDfGYys/iVH8keUtergGs9vI+ryy8c8WUHEg/Q7/V+1Khcg8ysTFIyUlh6zVJqVq7pdGgi\npcp3BpGKiGvldF1LtknfTOLaD6/lx4M/cvtnt3PPonucDkkCkLU2EVgIRBtj2uecN8aEA9cDW621\nf1upLiYmRj/T8jcPLX2IoWcMZdk/l7HiuhWc3ehsHv76YafD+ht9dqUoHo+HmJiYIsuo58gltJS3\niH84kHiAycsn88ttvxAZHkl8ajytXmzFTZ1vomXtlk6H5wiPlvIuVcaY0UBj72EdoKIxJmfFuR3W\n2tl5io8HBgCLjDHPAAnADUB94PyC6j/RHw4SmHYc3cFdPe4Csock9W/Sn9kbZ5/gVSLukvO3dmxs\nbKFlNOfIBTTGW8R/bD64mUvfu5Qtt23JPdfr1V48MfAJejfu7WBkztOco9JhjFkK9PUe5iSPnO+r\nx1rbP1/5VsDj3teEAGuAGGvtkgLqthMnTtTNOvmb8V+NZ+vhrbxz6Ttk2SwumXMJPU/ryUN9H3I6\nNJFiy7lZFxsbW2g+UuPIBdQ4EvEfqRmptHyhJTF9Y7iq/VUs3LqQWz+5lZ/H/kyNyjWcDs9Rahy5\nn/KRFCYlPYUr5l+BZ4eHLJvFkOZDeGP4G4QEhzgdmshJKyofqXHkAkpG4us8Ho/uMuex6eAmrph3\nBZsPbaZ5zea8cfEbdD+tu9NhOU6NI/dTPpKiWGv5M/lPjDHUrlI79/zBpIO8seENktOTGdZyGB3r\ndXQsRuUjKQ6tViciUo7a1W3Hxls2kmWzfGrzRBHQHFgpnDGGOmF1/nJuf+J+eszoQf8m/alTpQ6D\n3hzEu5e+y4CmAxyKUqRwxZkDq54jF9CdOhEJBOo5cj/lIzlZDy15iCOpR3hh6AsAzP95PlNWTGHZ\nP5c5HJk4bdWeVSzbvYzvd3/PwaSDhIWE0S+qHx3qdaB/k/6O3jxUz5GIiIiIlLpjx4/RJKJJ7nGT\niCbEp8Y7GJGUls9/+5wN+zdwNPUoK/5YwfHM4zSq1ohalWtRL7wekeGRVK5QmeqVqpNpM6lZuSZx\nKXEs+X0J3+36jm1HttGoWiP2Je6jQdUG7Nyzk8W/LyasYhghwSFccMYFxPaLpX7V+k5/qX+hnqMy\nYowZD3QCOgNRwE5rbZNCyupOnfg0jfH2D4dTDvPToZ+oF16P5jWbl3r96jlyP+UjOVmLti3iuo+u\n470R71E3rC43fnwj3Rp047GBjzkSj/LRqcvKyuKFH17g1XWvcjDpINZa+jbuy/wt8+lzeh8ysjJY\ns28NvU7rhWenh8jwSJKOJ5GSmULDqg35/ejvVK9UnaT0JKy1XNr6UuZsnsPOO3fS7H/N+G+///LK\nmleoVKESSWlJJKYlcjT1KE1rNOXO7ndyceuLOa3aaeXytRaVjzQYvuw8AkQDvwJH+P/lVkVEXMVa\ny4dbPqTlCy255eNbaP1ia8IeDeOidy/iYNJBp8OTcqZNYOVknNvsXB7p/wjXfHAN0a9H0yGyAw/3\ny94cNjMrkx1HdxCXHOdwlFKUuZvn0mBKAyo9UokHlzzIA+c8wLHjx6gUXIkjqUd4tP+j/BL3Cz//\n+TPzRs1j1d5VPD/keRLTEmlfrz3nNTuPg0kHmdh3IolpiSy9ZinpWelMu2AaAKv3riYjK4OI0Ah2\nHt3J2K5jiU+NZ2DTgTSJaMK+hH3cvehumj3XjC7TuvDyDy9TVjdpirMJrBpHZaeptbaOtfY8YJ/T\nwYiUJd2l812JaYkMeWsIl7x3CUlpSfx25DdeHPIitSvXJiMzgyFvDSHLZjkdppSjnAUZRIrr6g5X\ns3XcVnbftZunz3uaisEV2XNsD52ndeacmefQ5Lkm3Lvo3jL7gzcvfXaL52jqUe78/E7OmXkOYz4c\nw7uXvkuf0/sQ3TiaqWumUqdKHf7b/7/sOLqD1IxUIsMiqVShEsnpyWRkZdC1QVeCTTBn1j2TyLBI\naoTWoF+TfmTaTLo26ErFoIo8uORBalepze2f306QCWLammlkkUVEpQiS0pO4pNUl/JHwBzHRMQBc\n0OICfjz4I/csuofGzzbmnR/fKfWvOzo6Wo0jp1hrdzgdg4i4W1JaEvsS9pXLHwwFiU+N5/y3zmdf\nwj6MNbxy/iuEh4RzOPUw6VnpLP9jOev2rSN6VjSJaYmOxCgivummj29iWMth7L5rNzvu3MEX275g\n3s/znA5LgPTMdAbPHkxSWhKtarWiYdWGPPzNw4SHhDO4+WC+/+N7QoJDmLt5Li1rt+TJ5U+y9fBW\n0jPTGb1gNJ3qd+LKBVcSWiGUz377jC+3f8nhlMO8sPIF6oXXo89rfQgPCef9n9/nUPIhko4n0al+\nJ3Yf202wCeaaD68hyATx5sY3qRBUgWPHjxESHMJZ9c8iMyuTSf0mcSTlCDcsvIFL51zKn8l/luv3\nR40jESkxDcE5eY9/9zh1J9el/dT2dJrWid3xu8v1/fcm7KXztM5sOriJNnXbgIGFvy6kSsUqzFw3\nk/jj8cy6aBYhwSEYDGM/GVuu8YmIb1u3fx3Xd7oeYww1K9dkRJsRrNu3rszfV/noxNbvX8+x48eY\nduE0optE06BqA36J+4XLz7ycB5c8SEhwCAObDuTzbZ+zePtiGlZtyNDmQxnYdCBdGnQhLjmOpLQk\nDiYdZFf8LnYf2016Vjrzfp5HRlYGmw9tpnWd1oxoM4K4e+M4fP9hVt2wioP3HmT/PfuJ6RtDVEQU\n3+z6hqysLF5Z8wqNqzdm3k/zqFapGtuObCM1M5WRbUay8eBGes7oydGUo+X2/dFqdSIi5ezLbV8y\nY+0Mfh33K/XD63P3F3cz4I0B3N79dka1HUXdsLoAJKcnk5mVSdVKVUvtvfcl7OPFH17kqWVPERwU\njDGG8JBwnhv8HHd+cSdYyLAZtK/bnls+uYXgoGAS0xOZs3kODao2cGyStZQf7XMkpaFpjaZ88dsX\n3ND5BtIy01jy+xJGtx/tdFiCdzEC71T4UW1HMX3NdA4mHeSjXz4iOCiYcxqdQ7VK1dh0yyZa1m5Z\naD2ZWZkEmSDSMtNIzUilUoVKHEo6RP2q9akQVHATo3ZYbR7q+xAP9X2IRdsW8dq61/jglw/4M/lP\ngoOCuajFRUxdPZW7e93N4eTD7Dq6i0ybScsXWvLV1V9xZuSZJfratc+RSxhjNgFVrLVNC7mu1YFE\nAoC1lo9++Yhpa6YRHhLOnJFzWLdvHYPeHERCWgIj2ozgk62f0LdxX/Ym7OXHgz8SZIIYcsYQZg+f\nTeWKlU/5vQ+nHOaFVS8weflkUjJSCDbBnNfsPM5tfi53fn4nHSI7sPHARlrUapG72tCvcb/iGeNh\n2pppBAUF8flvn/P+yPfp2rDrKcWg1ercT/lISsumg5sYPHswzWs2Z2/CXtrVbcd7I98r9I9mKT8Z\nWRn0fq03bWq3YVjLYbz949v89OdPXHfWdbn7EJWnhOMJLN+9nMkrJvPNzm9Iz0zn1i638tLql3j9\n4teZvnY6W+O2cjT1KC8OfZHrOl1X4vcsKh+pcVQOitM4uuaaa4iKigIgIiKCjh075t61y2nh6ljH\nOvbN48ysTNp1a8cDix/A87WHWpVrsbHyRsafM54vFn9Bzco12V1zN1ERUaxetprgoGCSGybTpk4b\nLg+/nPu+vI/jjY7TuUFnbql9C6dXP/2k3n/1ntU8vudxjqYepcKuCjSo2oBH/vkIV86/koebPMxj\n3z1G265tiQyP5Pa6t5OWmUbMzhjW7lvLaXGnERQURMvOLVm9dzUdUzoyvvd4BvQfcML393g8zJo1\nC4CoqChiY2PVOHI5NY6kNB1JOcKafWuoGlKVrg27EmSCeHrF0zz23WOkZqRyebvLeX7I81SqUMnp\nUANOfGo8sV/HsjVuK+3qtuOhPg8RFhLmdFis2rOKqxdczZHUIzQIb0CvRr2YumYq0y+YzrjPxmGM\nYYwyVGcAACAASURBVECTAbw+/HUiQiNO+X3UOHKYeo7E33m0r0Shlu9ezoj3RpCcnkxCWgILRi3g\n/Bbnc9G7F/HZb59RNaQqQSaI6Rdm/+KPiY7hse8eI6ZvDA9/8zCJaYlc0+Eavt31LVedeRVTVkzh\np7E/EVohtFjvn5GVQf0p9XlswGM8t/I5HhvwGBe9exHfjPmGaz64hvjUeOJS4ri247W8MPSF3N6p\ntMw0Gj/bmMvbXs4bG9/g9m638+zKZ4mKiOK8Zued0vA69Ry5n/KRlKX3f3qfBxY/wMdXfEzNyjW5\nesHVtK3Tlh6n9eC5lc+RaTO5tuO1XN/p+lN+D+Wjwu08upPrPrqO9fvX07RGU6ZfOL3ce4mK43jG\ncW5YeANzNs+hZuWa3NLlFnYe3cmsDbO4uv3VzN8yn8oVKvPm8DcZ1GzQKb2H9jnyAdpXQsT/JKUl\nccmcS5h+4XTmjZrHmXXP5LqF13E45TAfXf4RkWGR9Dm9Dy1rt6RWlVoYDM9+/yxt67Tlu13fYcie\nyFytUjWiIqIY2XYkaZlpzFw7k7TMtGLFkLPKz4g2I9ibsJfIsEia1mjK0LeGsit+F1lk8cFlH/Dq\nRa/+ZdheSHAIn17xKTPXzyQ1I5Up30/hucHP8f6o93l9w+sn9X3wFGNfCRHxf19u+5Lbu99Oi1ot\nqF2lNrHRsczfMp87Pr+D+86+j9joWJ5c9iSz1s9yOlS/k5GVwflvn0+/qH5svnUzY7uOZchbQziS\ncsTp0P7m/9q77/Aoqv2P4++TkJCEUEPvoIBUkQ6ioIioYKEjgoWfFAUV1KtcG+AtNkTRq14FxIKg\nFANyVUSFoIAKShNCEaQoBAgtkN7O749dYgibkJDd7G7yeT3PPHFmzsw5G4f95sxppUuV5oO+H3B/\nu/uJT41n0fZFvLvpXab2nMq3e78lLCiMABPArR/fyr2f3UtiWqJb81flyEdoXQnxZ3p2Xdt7ai8V\nQyvSu3FvWldvzeH4w0SERrAhZgNvrHuDMsFl+GTgJ3Su3ZkB8wcQmxhLudLlGNN2DIt3LubQmUPs\nObGHGRtmcGerO2n5VkuOJhzl7Q1v0+ODHiSlJeWZf3xqPE+teIqTSSfp/l53nrzqSXp+2JO9J/eS\nYTN4uPPDbBy9kVua3OLy+itqXMETXZ+gT+M+xDwSQ4WQCjz57ZMkpCUUaGrV7vlYV0J8h17WiadU\nDqtMdGx01n50bDSJqYlM7j6ZPo37cF3D65jWaxpztsy56DwUj1w7EHeAuJQ4nrjqCaqFV+Ou1nfR\nKKIRG2I2eLtouXrlhleYfsN09p7YSwABrDmwhiMJRxjbfizHk47TtEpT3tv0HrWm1eLFNS/ma02+\n/LysU7c6DzHGDAfqOXcfAIKAac79fdbaOdnS2kmTJtFdswOJFCvHE49zyWuXsHH0RhpUbMCXv33J\nzfNuJsAEcEWNK/iw74c0jmiclf5k0kkmfjORrbFbaVixIb0u6cWcLXM4ePogcSlxpGemc0fLO3ih\n5wv0n9+frnW68kiXR3LNv//8/oSWCqVHgx48tOwhEtMSCQsKY2z7sTx25WNUDK14wc8QcyaGdjPa\n0aJqCzbGbMQYwyUVLyE2MZb1I9fnq893lHN2II058n3qVieeFJsQS+dZnWldvTWVQisRuSOSzrU7\n061et6zvso+2fMTcrXP5fOjnXi5t8XI88TgNpjdg70N7iQiLICU9hWZvNmP+gPm0rdnW28XL08aY\njfSf35+4lDgur3Y5W49uJcNmEBoYSmxiLGWCy3Aq+RTBgcHcdfldPN71cRpWdDmSJYvGHHmBMWYl\n0M25e/aXfPZ/QpS19tpsaRWMxK+pj3fu3lr/FlNWTaFznc6sO7iOce3H8fer/p7v6621LN6xmJFL\nR/KPa/7BmHZjMMbw0pqX+PPMn0y/YbrL61IzUgn/dzhn/n6G0qVKcyblDIMWDOLGRjfyYMcHC/QZ\n9p/aT7M3m3Ftg2sZ2mIot7e8nYELBnJdg+sY3W50vu+jMUe+T/FIPO1k0kkWRi8kOT2Z3o17E5cc\nx/VzrmdCpwmElArhhTUvMLffXCqGVmT1gdVUD69Ov6b98j3LneJR7p749gkW71hM38v6snLfSupV\nqMfcfnMxxve/lhPTEhn52UjmR88H4L529/Hp9k8pXao0wQHBpNt09pzYQ6AJxGK5su6V/Ovaf9Gl\nThcCzPkd5VQ58nEKRuLvFIzytu3oNqJjo2kU0YjW1Vtf1D2GRw4nrFQYb/Z+k2dXPcs/v/8nBsOA\nZgOYdcus82YZyrSZhP87nF0P7KJ2udpYa7nuw+sY03YMA5sPLHD+5Z8vz28P/EbVMlVJy0hj9P9G\n07JqSyZ0npDve6hy5PsUj8QbthzZwju/vENGZgbDWg1j36l9PLL8Efo17cfmI5spV7ocS29fmq8K\nkuJR7qy1LN21lE2HN3FJxUu4veXtLisOvmzZ7mUMWTiEpPQkqpapSnJ6MqkZqVQvU53dJ3dTu2xt\n0jLTiImPwWCoUqYKT131FGM7jD3ns6py5OPUrU5ELuRk0kkGLhjImj/WkJKewv3t7+eF615gxGcj\niAiN4M3eb553zXPfP8cHWz5gZJuR/HzoZ3Yc28GaEWsuar2kUUtHcejMIS6vdjlTf5hKWkYaTSKa\nsGzYMupVqJfntepW5z9UORJvs9ZS+aXKrLxrJa2qtSIjM4Mr372Sx698nL5N+3q7eOIDTqec5s7I\nO/ls52cYY2hRtQXHE48TEx9Dy6otiUuOY3/cfiqGVCSTTE6nnCasVBhj2o3hsSsfo0qZKqoc+ToF\nIxHJr9FLR9OwYkMe7/o4AL8e+ZVBCwexfez2c9LtPrGb9QfXs/P4Tk4knqBG2RqM6zCOsqXLXlS+\nKekp3LPkHhZEL6BNjTZMu34a3+3/jqW7lrL2/9bm6x5qOfJ9ikfibemZ6YT8M4Tkp5KzWopGLBlB\n59qdGdl2pJdLJ75k78m9jFo6iqj9URgMaZlpVAypSLPKzVgfs5565euRmp7K0cSjpGakkmkzMRgu\njbiUXQ/s0lTeIuI5mtmq6NQtX5dtsds4+wfspsObqFam2jlpluxYQudZnYncEckn2z4hLjWOiV0n\nXnTFCBxTq3au3Zl7r7iXn+79iSvrXsnDnR9m3cF16I9pEXGXUgGluKreVTz57ZMkpyez9o+1/G/X\n/+hSpwu7ju9i69GtpGem53q94tG5TiWfYtTSUXSY0YEhC4fw5+k/vV0kt2lQsQFf3/k1m8dspl/T\nfgSaQE4mn+SXw78QHhTO3lN7qV2uNuHB4ZQKKEWTiCYEmIAL/g5UOfIRmjpVRPJjXIdxbDq8iRs/\nupG7Ft/FI8sf4aWeL2Wdt9Zy79J7+Xzo58wfOJ+NozeyIWYDy/csL3TedcrXYd2hdVlrLK35Yw21\nytW64GBerXPkXxSPxNvm9Z/HLzG/UO65cgxaMIi3+rzFEyue4Nr3r6XvJ33pNLNTgZYTKKkybSa3\nfnwr1lqm3zCdJhFN6PFBDxJSE7xdNLdqVqUZHw/4mMQnEnm2+7OEB4VzKvkUxho2HN5ArbK1yLSZ\nxCbG0ialDaVW5T12Td3qfIC6MYhIQSSkJrBk5xIOnznMhpgNnEw5ScdaHZnYdSLWWsKfCyf1qdSs\nSstdi++iW71ujLhiRKHyzbSZDF00lC1HthBoAtl5fCdta7bljZveoE2NNhe8Xt3qfJ/ikfgSay3G\nGF5c8yKr9q8icnAkQQFBTPhqAqeST/Hebe95u4g+7UDcATrO7MjBhw9mTUbQaWYnnr/uebrX7+7d\nwnlQRmYGn+38jInfTGTXiV0EBQQRYAJIy0xjSHPHZA6RQyLVrU5EpLgoE1yGW5rcwjsb3qF62eqM\nbjuaH//8kXuW3EPpUqVpWbUlr/74KtZaomOj+Wr3V7Sr2a7Q+QaYAOb1n8cV1a8gJSOF2bfO5u7L\n76bXnF7sPrHbDZ9MROQvZ1/wbIvdRv+m/QkODMYYw+Dmg89ZTFZcCwoIIiU9hZT0FMDxgis+NZ7g\nwGAvl8yzAgMC6du0Lzsf2MkPI36gfa32WWOOlu5aypGEI3len79J40VE8qCpU4veqn2rqB5enanX\nTwWgZ8OeVH6pMqdTTrNg4AL6ftKXp1c+TYAJ4M3eb9KqWiu35GuMYeW+lawesTprkb1tsdtYFL0o\na5IIERF3alq5KYt3LObOy+8k0ASyaPsiLqt8GakZqSSkJlAhpEJWRUrx6C81ytagd+Pe9JnXh6Et\nhrL89+VEhEXQoVYHbxetyHSq04k1I9ZwOvk003+azowNM/jpz5/yvEaVIx8xefJkTeUtIvlmyb3r\n0yWVLmHzmM2cTjlNeHA4gQGBbs27VEApktKSsvaT0pLyXH/k7FTeIiIXY0KnCXy3/zsavd6IsKAw\nggKCGNJ8CBVfqEigCaRxRGM+HfQp8WnxRMdG0z61/Xlrv5VUs2+dzRvr3uD7A9/TokoLZt86O98L\n6hYn5ULK8XS3p3m629OkpKcQMikk17Qac+QD1MdbRAoqPjWedu+0o0/jPlxd72r++/N/qRhakY/6\nfeTxvF/54RXe2fAOg5sPZs6WORyIO8D97e7npetfIigwKNfrNObI9ykeia/KtJlEx0aTmpHK6eTT\nDF88nO/u/o76Ferz7KpnefPnNwkODKZqmaqcTDrJV8O+olFEI28XW3xUXvFIY45ERPxQeHA4q+5e\nRUJqAm//8jaXVrqU3o16M+zTYXSc2ZEB8wew58Qej+Q9vtN4RrYZyfOrn6du+brM7TeXrbFbeWT5\nIx7JT0QkwATQomoL2tRow4bDG+h3WT8aVGyAMYbq4dWJTYhl17hd/DLqFx7s+CBjPh/j7SKLn1Ll\nSEQKTV2mvKNaeDXe6vMWfS/ry8dbP+bhrx5m0fZFdKjZgfY123PN+9dwIumE2/M1xhBgArin9T2s\nuGsFA5oP4MO+HzJnyxy35yUiklPtcrX56eBPpGWkAbBy30oqhlYkNCiUqKgo+l7Wlx3Hdni5lOKv\nVDnyEVpXQkQuxvHE4zy6/FGWD1/O6ZTTRN8fzbyt8xjcYjAtqrYgal+UR/INDgzmdOrprP24lLhc\nZ0DSOkci4k79m/anVrlatHmnDf3n9+eL376gclhl4lPjAZj761xaVG3h5VJ6169HfmXJjiUe60FQ\nnJW8EVk+Sn84iD/TRCLec+jMIWqUrUGTiCZYLJVCK9E4ojEHTh3gTOoZj03ZOqj5IF5c8yITlk2g\naZWmTPthGo9d+ZjLtGcnm5kyZYpHyiIiJUtgQCALBi5g5d6VHE86zrTrp/Hc6udoML0BEaERACwb\ntszLpfSef6z6B//95b9cUf0Kfjr4E6/0eoVhrYZ5u1h+QxMy+AANgBWRixWfGk/D6Q35qN9HfLbz\nM6L2RbEvbh99GvVhx/EdrB2xltCgUI/kHXMmhhfXvsg3v39DSnoKLau15Lkez9E4orHL9JqQwfcZ\nY+ykSZM0e6r4pf2n9nMm9QyNKjWidKnS3i6OV0THRtPjgx5sHrOZqmWqEh0bTedZnTn08CHN4Mdf\ns6dOmTJFEzKIiOeoS6j3hAeHs2DgAoZFDmPh9oXsPbWXDrU60LBiQ6LuivJYxQgca2gcTTjKJRUv\nYdYts+hSuwvXvH8NRxOOeixP8byzS0uI+Jt6FepxLPpYia0YARyIO0CLqi2oWqYqAM2qNKN86fL6\nXnbq3r37BXtrqVudiIif61a/G39M+IMj8UeoWqZqkf1hkJqRysLohZx6/BShQaFcVe8q1v65lq/3\nfM0dre4okjKIiFzI7hO7GbFkBNtit9Ekogkzb5lJsyrNvF0sj2hepTkbYzayIWYDbWq04bOdn5Ge\nmU6tcrW8XTS/ocqRj9AisOLP9Nx6X3BgMHXK1ynSPANMAAZDYlpiVgtVfGr8eeOctAisiBSVnPEo\nJT2FG+bcwLgO41g4aCGR2yO58aMb2Xb/NsKDw71TSA+qU74OM26eQY8PehBSKoQAE8Cngz/12PjT\n4khjjnyAxhyJiL/62/K/sWr/Ku5rdx/rD61n5b6VrB+53uUfHRpz5PsUj6S4+fXIrwxaOIjtY7dn\nHWv3Tjv+c9N/6FS7kxdL5l4zfpnB5FWTSUhNoH/T/kzrNY0zqWeoHl6dUgFqC8lJi8CKiEepVaDk\nerHni9zb5l5W7FtB2eCyrL5ndbF8Gysi/iFnPKoQUoHYhFjikuMASEhNICY+hgohFbxQOs/4avdX\n/PP7f/LF0C/Y9cAujicd56kVT1G7XG1VjC6CfmMiInLRjDGMajuKUW1H5ZomKS2Jncd3FmGpREQc\n6pSvwx0t7+Dq966md6PeLN+znJsuvYkmEU0AWH1gNc+sfIZTyafo3ag3k7pP8rsKxde/f8197e7j\n8uqXA/DvHv/mto9v83Kp/Jd//d8XEZ+kMUeSnJ7Mir0rSElPoVv9blQKrQTAzmM7ufGjGwkLCvNy\nCUWkJHAVj1694VUid0Sy7eg2HrvyMQY2G4gxhm1Ht9H3k768dsNrNKzYkInfTuTxrx/n5V4vF33B\nL0JSWhK7ju8i0ASe8wJqx7EdWd/BUnAac+QD1MdbRPzZ6ZTTXPP+NZQOLE35kPJsO7qNFXet4NJK\nl9L13a4MaTGEcR3GacyRH1A8kpLk39//m+OJx7MqQ3tP7qXr7K4cfPigl0t2YdGx0dz00U2EBYUR\nEx9DgAngmvrXUCO8Bh9v+5gFAxfQvX53bxfTZ2nMkYh4lMYclWwvr32ZFlVbsGbEGr6840se7Pgg\nf/v6b4AjgA9uPtjLJSzZjDF/N8YsMMb8bozJNMbs9XaZRDylIPEopFQIJ5NPZu2fSDpB6UD/WCNp\neORwnrzqSaLHRrPnwT1UCq1EzbI1qV+hPqvuXqWKUSGociQiIoVy4PQButbpijGOl3Bd63blz9N/\nAnBZ5ctYtH2RN4sn8C+gO/AbcBJQ05AIMKzVML75/RsmLJvAf9b9hwELBjCx60QANh/ezLXvX8tl\n/7mMEUtGcDrltJdLe65tR7cxtOVQACqFVuKmS2+ifoX6PNLlkWK7hlNRUeXIR0yePFlv38VvacxR\nydapVidmbZzFqeRTpGWk8fq61+lYqyMAs26ZxaT3JlG9T3Uvl7JEa2itrWKt7QXEeLswIp5UkHhU\ntUxVfrz3R4IDg9l6dCuv3fAao9qOIuZMDL3m9OKOlnewaNAiLJahi4Z6rtAXoUnlJkTuiAQcXZu/\n3ftt1iQTUjgac+QD1MdbRPxZps1k/LLxzNgwg0ATyDUNrmFe/3lZU3onpCYQHRtNh9odNObIy4wx\nW4Ewa23DXM4rHkmJN+/XeSzcvpBFgxyt3mkZaZR9riwnHz+ZteC1N8z9dS5Pr3ya+NR4rqp7FWv/\nWEutcrX4I+4PBjcfzKs3vJrVgi95y2vMkWarE5FCi4qKUutRCRZgAnjtxtd4/rrnSc9Mp1zpcuec\nLxNchva12nupdCJSkrgjHoUFhXEk/gjWWowxnEg6AUBQYJAbSnhxVu1bxd++/huLBi2idrna3P/5\n/fRp3Id7Wt9DRFgEjSMae61sxY261YmIiFuEBYWdVzESEfE3N1x6A2mZaQxeOJiX175Mjw96MLHr\nREoFlOLDzR/ScWZH2r3TjrfWv0VRtbQu272M0W1H06l2J2qXq83U66fyze/f0LlOZ1WM3EwtRyJS\naGo1EhERX+COeFS6VGlW3LmCN9a/wYG4AzzT7RkGNhtI5PZInl75NLNumUVQYBCjlo6idKnSjLhi\nROELfgEVQioQfSw6a3/3id1UCKng8XxLIlWOREREJMvdd99N/fr1AahQoQKtW7fO+oPz7MRB2td+\ncd8vE1yGDmkdIBS6N3ecf33+6wypNYQeDXsAcGe5O3lr4VtZlSN3l+fbFd+yfM9ygi8J5tKIS3l9\n/utcu/la2nZuy/ub3+fRmo+e043Ql35/vrYfFRXFe++9B5D1/ZYbTcjgAzQAVvxd9i9nKbk2xGzg\n7sV389uJ32hZtSUzbp7B5799zo5jO2hepTkTr5qoCRm8TBMySHHnyXh09+K7aVWtFQ93fhiA2Rtn\ns2TnEhYPWez2vDJtJgMXDORY4jGuqX8Ni7YvokeDHjSs2JD41HhuvPRGrqhxhdvzLSk0IUMRM8YE\nAA8Bo4F6QCwwH3jGWpvozbKJiHhCXHIcfeb2Yer1U7mlyS28v/l9Os/qTI8GPejXtJ/WOvIjkydP\npnv37nrhIZLDhE4TuO7D6zidcprgwGBe+fEVIgdHsiFmAw9/9TCH4w/TtW5XXun1CmVLly1UXusP\nrufXI7+y9f6tBAcG82DHB6n/an32jd9HpdBKbvpEJU9UVFRWi1Ju1HLkAcaY6cADwKfAl0Az5/73\nwHU5X8vpTZ2I+Lvv93/PY988xg//9wMAW45sod077fh51M+0qtaK1IxUSpcqrZYjL1PLkUjhbDu6\njXc3vkuGzWB4q+FUD69Om3fa8MJ1L9C+ZnteWPMCcSlxLBmypMD3TkxL5Lnvn2PH8R2EBIawP24/\n393zHQDWWmq/Upu1I9ZSr0I9d3+sEkctR0XIGNMcR0VokbV2YLbje4HXgCHAPC8VT0TEIyLCIjgQ\nd4D41HjCg8OJTYglw2ZQJawKAEEB3psCt6QzxgzH0YsBoAoQZIx5yrm/z1o7xzslE/E/zas25+Ve\nL2ftf7j5Q7rX787dre8GYOYtMyn7XFkSUhP4cveXHE04Ste6XWlVrVWe9820mdw872Yqh1Wmf9P+\nzP11LusOrmPGLzO4sdGNzPhlBlXCqlC7XG1PfjxBU3l7wu3On6/mOD4DSASGFW1xRDzvQk3UUvw1\nq9KMW5vcSpdZXRi/bDxjPh9D5bDKTF07lTUH1vDQsoe8XcSSbATwrHOrDJTPtn/eNFuTJ0/Wv2nx\nW0X97GZfEwngWOIxAkwA/T7px7QfprH58GZ6ftiTRdF5dy3ednQb+07tY26/uQxpMYRPB39KpdBK\nvPnzm7Sf0Z71h9bzv6H/IzAgsCg+VrEVFRXF5MmT80yjbnVuZoz5CrgWR7eFtBzn1gCNrLVVcxxX\nNwbxa5qQQcDR7WPprqX8dvw3WlZrSevqrfnb13/LmpBh9m2z1a3Oxykeib8r6niUnJ7MVbOv4pKK\nl9C+ZntmbZzF5dUvZ/+p/Xx/z/cEBgSy/uB6bp53M4cfPXze9RmZGSSkJbDnxB5uX3Q728duxxhD\nps2k0euNWDJkCS2qtiiyz1NS5NWtTpUjNzPG/ApUttbWcHFuPjAACLbWpmc7rmAkIsVeXsFIfIPi\nkUjBxafG88a6N4iJj6Fr3a4ciT/CliNbePvmtwFIzUilzL/LsHDgQj7Y8gGlAkoxrv04/jz9J2M+\nH0N6ZjoNKjQgwATQtW5X+jftzyfbPmHr0a18d893lArQKBh305ijohUGpORyLjlbmtNFUxwREZH8\n02x1IgUTHhzO410fz9rfcmQLz373LCPbjqRVtVY8teIpmlZuytgvxvL8dc+TlJbEbR/fRiaZrL5n\nNS2qtmD6T9OZtXEW1lr+8d0/aF6lOV/c8YUqRm6m2eq8IB8tR/2B0jlbju666y4tuqd9v93ftGkT\n48eP95nyaN839qNyLLo3ZcoUtRz5OLUcib+L8pFu3gu2LeCBLx/geNJxrq53NekZ6YzvNJ6+TfsC\ncMeiO/jx4I/seXAP4OiWHPKvEE4+fpKwoDBvFr1EULe6IpSPMUeXWmur5TiuYCR+zVeCkfg2davz\nfYpH4u98LR5l2kwCTAA9P+zJ2PZjue2y2wAY98U45myZQ8wjMYQGhfLLoV/o+WFPjj92HGP0Nelp\n6lZXtNYBPYGOwOqzB40xIUBrIMrVRerGIP5Mz63kJSof3RhERNzB1+JRgHFMDD267Wge+PIBktKS\nSEpPYv62+XSu3Zk277ShVbVWrNi7gpm3zFTFyAeo5cjNjDEtgM1ApLV2QLbjDwDTgWHW2rk5rtGb\nOhEp9tRy5PuMMXbSpEl6WSfiAZHbI3l/8/uOCRk6jKNbvW6s3LeSmDMxtKvZjiaVm3i7iMXe2Zd1\neXXzVuXIA4wxrwHjgEjgS6ApjoVhV1trr3WRXsFI/JqvdWMQ35KfYCS+QS/rxN8pHkl+aMxRETPG\nBADjgVFAfSAW+AR4xlqb6CK9gpH4NQUjyQ+1HPk+xSPxd4pHkh+qHPk4BSMRKQlUOfJ9ikciUhLk\nFY8Cirow4trkyZM1YFlEiqWoqCgmT57s7WJIPikeiUhxlZ94pJYjH6A3deLv1I1B8kMtR75P8Uj8\nneKR5IdajkRERERERC5A6xz5CK1zJP5Mz63k9NXur3j5h5dJTk+mVVIrIo5EeLtIIlICKB5JYalb\nnQ9QNwYRKS5OJp3kxo9u5KeDPxEUEMTwVsPZcHgDQ1sM5bGuj6lbnY9TPBKRkkDd6kTEozR4W84a\n+8VYktOTeebqZ9g3fh+r/1jNsFbDeH/z+94umuSTJmQQf6ZnV/KSnwkZVDnyEQpGIlIcrP1jLV1q\ndyEpPYmaZWtyZ6s7+WL5Fxz74pi3iyb5dLabt4hIcdO9e3fNVucP1I1BRIqLLrO6MKDZAF5Y8wL3\ntbuPyB2R/BH3B9NvmM6dre9Utzofp3gkIiWBFoH1cQpGIlJc/HzoZ3rP7c0V1a9g0+FNWCxv3PQG\nA5oN0FTefkDxSERKAlWOfJyCkfg7rSsh2R08fZDv9n9HeHA4vS7tRXBgMKB1jvyB4pH4O8UjyY+8\n4pGm8hYREbeqVa4Wt7e83dvFEBERKTBNyOAjNCGD+DO9pZO85Gd2IBERd1A8ksJStzofoG4MIlIS\nqFud7zPG2EmTJmlRchEplqKiooiKimLKlCkac+TLVDkSf6c+3pIfqhz5PsUj8XeKR5IfWgRWRERE\nRETkAtRy5AP0pk5ESgK1HPk+xSMRKQnUcuQHNCGDiBRXmpBBRET8hVqOfIDe1Im/Ux9vyQ+1RveG\nhAAADidJREFUHPk+xSPxd4pHkh9qORIREREREbkAtRz5AL2pE5GSQC1Hvk/xSERKArUciYiIiIiI\nXIAqRyJSaJpMREREfIHikRSWKkciIiKSRbOnikhxlZ/ZUzXmyAeoj7eIlAQac+T7FI9EpCTQmCM/\noDd1IlJcaZ0jERHxF2o58gF6Uyf+TutKSH6o5cj3KR6Jv1M8kvxQy5GIiIiIiMgFqOXIB+hNnYiU\nBGo58n2KRyJSEqjlSERERERE5AJUORKRQtNkIiIi4gsUj6SwVDkSERERERFBY448whgTADwEjAbq\nAbHAfOAZa22ii/Tq4y0ixZ7GHBU9xSMRkfNpzFHRewV4GdgKjAMWAA8CS40x+sNARESKiuKRiEgB\nqHLkZsaY5sADwCJr7QBr7Sxr7SPAw8A1wBCvFlDEA9THW8T3KB5JSaR4JIWlypH73e78+WqO4zOA\nRGBY0RZHxPM2bdrk7SKIyPkUj6TEUTySwlLlyP3aAxnAuuwHrbUpwGbneZFi5dSpU94ugoicT/FI\nShzFIyksVY7cryZwzFqb5uLcQaCyMaZUEZfJK7zZtO2pvN1x34u9R0Gvy2/6/KQrid0UvP2Zi+Mz\nLEVO8chJ8ci991A8Klre/szF8RnOiypH7hcGpORyLjlbmmJPwci99/DlYLRv37585eVPFIw8ew8p\nEopHTopH7r2H4lHR8vZ3bnF8hvOiqbzdzBjzK1DZWlvDxbn5QH+gtLU2Pdtx/U8QkRJBU3kXHcUj\nEZHc5RaPSkRzehE7BFxmjAly0ZWhFo4uDunZD+qPBRER8QDFIxGRAlK3OvdbBwQCHbMfNMaEAK2B\nn71RKBERKXEUj0RECkiVI/f7BLDA+BzHRwKhwEdFXiIRESmJFI9ERApIY448wBjzGo6VyCOBL4Gm\nOBbiW22tvdabZRMRkZJD8UhEpGBUOfIAY0wAjjd1o4D6QCyON3jPWGsTvVg0EREpQRSPREQKRt3q\nCsEYE2CMmWCM2WGMSTLGHDDGTAVCrLXTrLWXWWtDrLV1rLWP5gxExphBxpjZxpjNxpg0Y0ymMaZu\nHvmVN8a8bow56MxvqzFmjMc/qBRbxpi/G2MWGGN+dz5/ey/yPjcZY9YaY+KNMceNMfONMfVzSavn\nWMQDXMUk4EXgvxeKR87rFZPEaxSPxFeoclQ4rwAvA1txdFtYADwILDXG5GfGn/uAQUACsBtH33CX\njDHBwNfAaGCeM7+dwJvGmEmF+AxSsv0L6A78Bpwkj2cwN8aYfsD/gNLAo8BLwNXAGmNMjRxp9RyL\neI5ikvgzxSPxDdZabRexAc2BTGBBjuPjnMdvz8c96gABzv/+j/O6urmkvd95fmyO4wtxLPLn8jpt\n2vLagPrZ/nsr8HsBrw8CDgJ7gbBsxy8H0oG3c6TXc6xNmwc2xSRt/r4pHmnzlU0tRxfvdufPV3Mc\nnwEkAsMudANr7R/W2sx85jcUx9u8GTmOv4rjC2FwPu8jksVau6+Qt+gG1ABm2mzddKy1m4EoYLAx\nJvt6anqORTxDMUn8muKR+ApVji5eeyADxzoSWay1KcBm53m3cA6obQNstNam5ji9HkfTczt35SdS\nAGef8x9cnPsJKAc0Bj3HIh6mmCQlneKRuIUqRxevJo7VxXOuOg6OZt3KOd5QFEZFIMR533M4A99x\nHKudixS1ms6f5z2b2Y6dTaPnWMRzFJOkpFM8ErdQ5ejiheHok+pKcrY07sqLC+TnrrxECiKvZzPn\nvwM9xyKeo5gkJZ3ikbiFKkcXLxHHbCiuhOBoknXXGhJn75NXflqvQrwhr2czJEcaPccinqOYJCWd\n4pG4hSpHF+8Qjm4KQS7O1cLRvSHdTXmdBJJw0cRrjCkNVMZ1M7KIpx1y/nTV/eDssbPPpp5jEc9R\nTJKSTvFI3EKVo4u3DggEOmY/aIwJAVoDP7srI+fsQRuANs55+bPr4PzptvxECuDs4O8uLs51AuKA\nXaDnWMTDFJOkpFM8ErdQ5ejifYKjm8L4HMdHAqHAR2cPGGOqG2MuM8aEFiK/eTj6v47KcXw8kOYs\nj4jH5PIcrwJigHuNMWWypb0cx2J+C6y1GdnS6zkW8QzFJCkxFI/Ek4y1BV6AWJyMMa/hWGAvEvgS\naAo8AKy21l6bLd17wJ3ANdbaVdmOX41j5WaAPjjeVryM4+2Gtdb+K1vaIGAtjsXMXgN2ADcBtwH/\nsNZqNWcpMGPMcKCec/cBHGs7THPu77PWzsmW9j1cP8cDcASRzcBMHNOlTsAxrXBba21MtrR6jkU8\nRDFJ/JnikfgMb69C688bjpa3h3H8g0oG/gCmkm1lZme62Tj+YV6d4/gkHKszZzrPZ2Tfd5FfeeB1\nHP1gk3GsIH2/t38P2vx3A1bm8QyuyJHW5XPsPNcbx9oSCcAJYD7QIJc89Rxr0+aBTTFJmz9vikfa\nfGVTy5GIiIiIiAgacyQiIiIiIgKociQiIiIiIgKociQiIiIiIgKociQiIiIiIgKociQiIiIiIgKo\nciQiIiIiIgKociQiIiIiIgKociQiIiIiIgKociQibmCMecEY87sxplQR5dfaGJNpjLm6KPITERH/\noHgkhaXKkUg+GWO6O78Ac9vSvF1GbzDGNAAeBKZYa9MLeG0pY0yMMeZIXoHMGNPQ+Tv+CsBauwmI\nBF4uTNlFRPyR4pFrikfiDkVSqxYpZuYCX7g4nlnUBfERE4E4YE5BL7TWphtj3gMeB/oAi3NJerfz\n56xsx14FVhljbrLWuvr/ISJS3CkenUvxSArNWGu9XQYRv2CM6Q6sAB611k7zcnHOYYwJBVKttRlF\nnG854BAww1o74SLv0QjYCSy11t7q4nwAsBcIA2paa9Oynfsd2GqtveVi8hYR8UeKRy7zVTwSt1C3\nOhEPMMbUdza7TzLG9DHGrDfGJBljDhljXjTGBLq4ppEx5kNns36KMWavM21YjnTvOe9d2RjzrjHm\nCBAP1HKeb2WMWW6MiTfGHHOmr+y8ZrYzTVVjTKoxxuXbNWPMG8aYDGNM3Qt81JtwBAmXb8ry85ms\ntb8B3wM3GmOqubhND6AOMDd7IHL6CrjBGFPmAuUUESmRFI/y/5kUjwTUrU7kYpQxxlR2cTzFWnsm\nx7GbgPuBt4CZwG3Ao8BJ4LmziYwxbXG8BTzhTHsQaI2j7/SVxphuLvpPfw3EAFOAMkCC863X987z\n05336Q186TxmAay1R40xS4B+xpjy1tq4bGUJAYYCX1trD1zgd9HN+XN9zhMF/EzvAlcBw4GpOW51\nj/PnLM73IzAa6IojMImIlCSKR39RPBL3sNZq06YtHxvQHUc/7ty2z7Klre88dgaom+M+vwKHchzb\nDEQDZXIcv815n7uyHXvPeewDF2WcD2QAnXMc/9h5zbvZjvV0HrsvR9o7nMcH5ON3sgo4lsu5gnym\nMBz9xLflSFsBSAJ+ziWPrs57TfD286FNmzZtRbUpHrn8nSgeaXPLpm51IgX3NnCdi+1JF2kX2/Pf\ndkUB1c825RtjWgItgXlAqLPLQWXn28A1QCJwvYt7n/NGy9k14iZgnbX2hxxpz5tFx1r7NY6+0/+X\n49T/AcfIfTBqdlVwvIk7R0E/k7U2EUfAbGqM6ZDtVkOA0rh+Swdw3Pmzaj7KKiJS3Cge/UXxSNxC\n3epECu43a+2KfKb93cWxs1+gETi+lJs696c4N1dcfdnuyrFfBccbr535SHvWTOBfxpjLrbWbjTEN\ncXRNeNXmbxpUCxgXxy/mM80CRgIjgHXOYyNwvKmbm8s9zuatmWVEpCRSPPqL4pG4hSpHIp6V12w9\nJsfPqcCyXNKezHnAWptciHKd9S6OYPF/OPpej3CWZ2Y+r48FWrk4XuDPZK1dZ4zZBgw2xowHLgHa\n4Rj4GufqBkClbOUQEZHcKR4pHkk+qHIk4n1n36JlFuANoCuxQALQxMU5V8ew1h4xxiwF7jDGTMSx\nfsOP1trt+cxzK3C1MaaStTZ7d4aL/Uzv4uhy0Q9ok+1Ybi7NVg4RESkcxaO/KB6VUBpzJOJl1tqN\nOL5MxxjH6t7nMI5VuyvmvMzFfTJwzALU0RjTJcfpR/IowgygIo6+6zXJ/1s6gJXOn51zlOViPhPA\nh0AaMArHQNy9FwhmnZzp1xSgzCIi4oLi0TkUj0ootRyJFFxbY8ywXM5FWmsTLuKew3FMM7rFGPMu\njll1wnC8ieqLY9XvD7Kld9WvGuApoBewzBjzH/6aOrWK87yrvtBfAftxfPmfwTEQNb+WOa+5Cfi8\nkJ8Ja+0xY8xnQH/noUm5ZWyMMcANwDLnAFoRkZJG8egvikfiFqocieTf2S/yIcDtuZxfi+tBrznT\nnRMUnINPrwD+DtwCjMHxJb8XmA18m9f12e6zyxhzNY6+1Q8ByTgWxBsL7MExmDTnNdYYMwt4Fphf\nkC92a22Cc+G+wcaY8TbbongF/EzZzcIRjDJwTBObm6uBusB9+S2viEgxoXh0/rWKR+IWxlpNqiFS\n3DkXwFsPTLTWvuji/GPA8zjWo/ipgPeuB+wAxllrc5vi1O2MMZFALWtthwsmFhERn6B4JL5OlSOR\nYsYYE2qtTcq2b3B0TRgItHX2v86evhSO6VbPWGtbX2SezwGDgcb5nHK1UJxvAH8Gultrv79QehER\nKXqKR+KPVDkSKWaMMTtxdBHYCpQBbsaxcvfH1tqh2dLVB7oAt+IIVEOstfOLurwiIlI8KR6JP9KY\nI5HiZzGOADQcx7/x33EMjH0hR7ruOKYljQWmKBCJiIibKR6J31HLkYiIiIiICFrnSEREREREBFDl\nSEREREREBFDlSEREREREBFDlSEREREREBFDlSEREREREBID/B2vI0dEQw+VkAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_fit(w_exp,eps_exp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets define a Drude-Lorentz model. This will be the model we will try to fit the experimental data with. The model has one Drude term and five Lorentz oscillators" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def drude_lorentz_model(params, w):\n", " '''\n", " This is functional representation of drude lorentz model with the first term being Drude\n", " and the rest terms being lorentz oscillators\n", " '''\n", " \n", " # Get the plasma frequency\n", " omega_p = params['Omega_p'].value\n", " \n", " \n", " # Get the Drude parms\n", " f0 = params['f0'].value\n", " gamma0 = params['Gamma0'].value\n", "\n", " # Get the first Lorentz term parameters,\n", " f1 = params['f1'].value\n", " gamma1 = params['Gamma1'].value\n", " omega1 = params['Omega1'].value\n", " \n", " # Get the second Lorentz term parameters\n", " f2 = params['f2'].value\n", " gamma2 = params['Gamma2'].value\n", " omega2 = params['Omega2'].value\n", " \n", " # Get the third Lorentz term parameters\n", " f3 = params['f3'].value\n", " gamma3 = params['Gamma3'].value\n", " omega3 = params['Omega3'].value\n", "\n", " # Get the fourth Lorentz term parameters\n", " f4 = params['f4'].value\n", " gamma4 = params['Gamma4'].value\n", " omega4 = params['Omega4'].value\n", " \n", " # Get the fifth Lorentz term parameters\n", " f5 = params['f5'].value\n", " gamma5 = params['Gamma5'].value\n", " omega5 = params['Omega5'].value\n", "\n", " # Drude component\n", " epsilon_D = 1 - (f0 * omega_p ** 2 / (w ** 2 + 1j * (gamma0) * w))\n", "\n", " # Lorentz first oscillator\n", " epsilon_L1 = (f1 * omega_p ** 2) / (omega1 ** 2 - w ** 2 - 1j * gamma1 * w)\n", "\n", " # Lorentz SEcond oscillator \n", " epsilon_L2 = (f2 * omega_p ** 2) / (omega2 ** 2 - w ** 2 - 1j * gamma2 * w)\n", "\n", " # Lorentz Third oscillator \n", " epsilon_L3 = (f3 * omega_p ** 2) / (omega3 ** 2 - w ** 2 - 1j * gamma3 * w)\n", " \n", " # Lorentz Fourth oscillator \n", " epsilon_L4 = (f4 * omega_p ** 2) / (omega4 ** 2 - w ** 2 - 1j * gamma4 * w)\n", " \n", " # Lorentz Fifth oscillator \n", " epsilon_L5 = (f5 * omega_p ** 2) / (omega5 ** 2 - w ** 2 - 1j * gamma5 * w)\n", "\n", " # Sum all the terms\n", " epsilon = epsilon_D + epsilon_L1 +epsilon_L2 + epsilon_L3 + epsilon_L4 + epsilon_L5\n", "\n", " return epsilon" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets define a callable function which returns a residual. This function will be used by our minimization algorithm. It is important to note how the residual is calculated as it is little bit different than normal fitting examples shown in lmfit, because we are minimizing complex data." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def complex_residuals(params, model, w, exp_data):\n", " '''\n", " This is the residual function that we will try to minimize.\n", " It takes the params dict that has the parameters that need to be found. \n", " '''\n", " if model == 'Drude':\n", " # Our Drude model\n", " epsilon= drude_model(params,w)\n", " elif model == 'Drude-Lorentz':\n", " # Our Drude model\n", " epsilon= drude_lorentz_model(params,w)\n", " \n", " # Lets calculate our complex residual as the way it is done in page 5264\n", " # Rakic, a D.,et al (1998). Optical properties of metallic films for vertical-cavity optoelectronic devices. Applied Optics, 37(22), 5271–83.\n", " residual = (abs((epsilon.real - exp_data.real)/exp_data.real) + abs((epsilon.imag - exp_data.imag)/exp_data.imag))**2\n", " \n", " # if the residual is being used for least square optimizaiton we should have used\n", " # residual = (abs((epsilon.real - exp_data.real)/exp_data.real) + abs((epsilon.imag - exp_data.imag)/exp_data.imag))\n", " # least square method does the square of the residual in its algorithm. see http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.leastsq.html\n", "\n", " return residual" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is where differential evolution will be used to make a fit to the data. Differential evolution algorithm requires bounds for each parameter. I gave resonable bounds based on looking at the experimental data and hope to see what the algorithm does. It seem to do a resonable job and much better job than least squares method. **The fitting results will vary a little bit each time we run this cell. So run the cell couple of times to see if things change for better.**" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "differential_evolution step 1: f(x)= 182.623\n", "differential_evolution step 2: f(x)= 182.623\n", "differential_evolution step 3: f(x)= 182.623\n", "differential_evolution step 4: f(x)= 50.6673\n", "differential_evolution step 5: f(x)= 20.0853\n", "differential_evolution step 6: f(x)= 20.0853\n", "differential_evolution step 7: f(x)= 20.0853\n", "differential_evolution step 8: f(x)= 20.0853\n", "differential_evolution step 9: f(x)= 15.1101\n", "differential_evolution step 10: f(x)= 15.1101\n", "differential_evolution step 11: f(x)= 11.3224\n", "differential_evolution step 12: f(x)= 8.07834\n", "differential_evolution step 13: f(x)= 3.05303\n", "differential_evolution step 14: f(x)= 3.05303\n", "differential_evolution step 15: f(x)= 2.4401\n", "differential_evolution step 16: f(x)= 2.4401\n", "differential_evolution step 17: f(x)= 2.4401\n", "differential_evolution step 18: f(x)= 2.4401\n", "differential_evolution step 19: f(x)= 2.31825\n", "differential_evolution step 20: f(x)= 2.26033\n", "differential_evolution step 21: f(x)= 2.09022\n", "differential_evolution step 22: f(x)= 2.00165\n", "differential_evolution step 23: f(x)= 1.58789\n", "differential_evolution step 24: f(x)= 1.57863\n", "differential_evolution step 25: f(x)= 1.57519\n", "differential_evolution step 26: f(x)= 1.56108\n", "differential_evolution step 27: f(x)= 1.51249\n", "differential_evolution step 28: f(x)= 1.47025\n", "differential_evolution step 29: f(x)= 1.43602\n", "differential_evolution step 30: f(x)= 1.43602\n", "differential_evolution step 31: f(x)= 1.43602\n", "differential_evolution step 32: f(x)= 1.43602\n", "differential_evolution step 33: f(x)= 1.37272\n", "differential_evolution step 34: f(x)= 1.37272\n", "differential_evolution step 35: f(x)= 1.36562\n", "differential_evolution step 36: f(x)= 1.32582\n", "differential_evolution step 37: f(x)= 1.32201\n", "differential_evolution step 38: f(x)= 1.31609\n", "differential_evolution step 39: f(x)= 1.29776\n", "differential_evolution step 40: f(x)= 1.28905\n", "differential_evolution step 41: f(x)= 1.28905\n", "differential_evolution step 42: f(x)= 1.28905\n", "differential_evolution step 43: f(x)= 1.27764\n", "differential_evolution step 44: f(x)= 1.27764\n", "differential_evolution step 45: f(x)= 1.27625\n", "differential_evolution step 46: f(x)= 1.26961\n", "differential_evolution step 47: f(x)= 1.25183\n", "differential_evolution step 48: f(x)= 1.24244\n", "differential_evolution step 49: f(x)= 1.23059\n", "differential_evolution step 50: f(x)= 1.22122\n", "differential_evolution step 51: f(x)= 1.214\n", "differential_evolution step 52: f(x)= 1.21086\n", "differential_evolution step 53: f(x)= 1.20955\n", "differential_evolution step 54: f(x)= 1.20823\n", "differential_evolution step 55: f(x)= 1.20721\n", "differential_evolution step 56: f(x)= 1.20556\n", "differential_evolution step 57: f(x)= 1.20399\n", "differential_evolution step 58: f(x)= 1.20319\n", "differential_evolution step 59: f(x)= 1.19969\n", "differential_evolution step 60: f(x)= 1.19968\n", "differential_evolution step 61: f(x)= 1.1984\n", "differential_evolution step 62: f(x)= 1.19704\n", "differential_evolution step 63: f(x)= 1.19592\n", "differential_evolution step 64: f(x)= 1.19556\n", "differential_evolution step 65: f(x)= 1.19137\n", "differential_evolution step 66: f(x)= 1.18851\n", "differential_evolution step 67: f(x)= 1.18354\n", "differential_evolution step 68: f(x)= 1.18286\n", "Print exited successfully? : True\n", "Termination Status: Optimization terminated successfully.\n", "[[Fit Statistics]]\n", " # function evals = 66984\n", " # data points = 69\n", " # variables = 18\n", " chi-square = 1.137\n", " reduced chi-square = 0.022\n", " Akaike info crit = -226.448\n", " Bayesian info crit = -186.234\n", "[[Variables]]\n", " Omega_p: 9.97441670 (init= 9)\n", " f0: 0.68683341 (init= 0.5)\n", " Gamma0: 0.01987236 (init= 0.03)\n", " f1: 0.98709783 (init= 0.5)\n", " Gamma1: 1.4876e-05 (init= 0.5)\n", " Omega1: 8.02150038 (init= 5)\n", " f2: 0.07682351 (init= 0.5)\n", " Gamma2: 0.65934600 (init= 0.5)\n", " Omega2: 2.91927926 (init= 3)\n", " f3: 0.16124225 (init= 0.5)\n", " Gamma3: 1.12448236 (init= 1)\n", " Omega3: 3.70495567 (init= 4)\n", " f4: 0.15690286 (init= 0.5)\n", " Gamma4: 1.18873234 (init= 1)\n", " Omega4: 4.50000000 (init= 5)\n", " f5: 0.23132574 (init= 0.5)\n", " Gamma5: 1.88224399 (init= 1)\n", " Omega5: 5.71650394 (init= 6)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0cAAAFRCAYAAAC/qtYsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcjXX/x/HXd2aMsXPbFdn3W/ZCmLGnBQktKpWokNa7\nVMwMWoS0y1JJ/IiksqUsR9mTXYlIhCQJYwxj5vv745xxzz2NMZiZ6zrnvJ+Ph8eZ67q+5zofNeYz\nn+u7GWstIiIiIiIiwS7E6QBERERERETcQMWRiIiIiIgIKo5EREREREQAFUciIiIiIiKAiiMRERER\nERFAxZGIiIiIiAig4khERERERARweXFkjBlkjJlpjNltjEk2xvxygfbVjDGfGWP+MsbEGWO+McZE\nnadtiDHmMWPMdmPMKWPMXmPMKGNM3py+t4iIuJvykYhIcDBu3gTWGJMMHAHWAw2BY9baiudpWwlY\nC5wBXgOOAw8AtYHrrbWL07R/HRgAfAosAGr6jr8F2thU/2Gy894iIuJ+ykciIsHB7cVReWvtHt/X\nW4G8GSSjGUAXoIG1drPvXD5gG5Bgra2eqm0tYAswy1rbLdX5/sAbwJ3W2mk5cW8REXE/5SMRkeDg\n6mF1KYnoQnyJ4WbAk5IsfO8/CUwEqhpjGqV6y+2+19fS3GoCEA/0zIl7i4iIf1A+EhEJDq4uji5C\nHSAcWJXOtTW+14apzjUCkvAOTTjHWnsa2OS7nhP3FhGRwKJ8JCLixwKlOCrje92fzrWUc1ekaf+n\ntTbxPO2LGWPCcuDeIiISWJSPRET8WKAURymr7pxO51pCmjYpX6fXNr322XlvEREJLMpHIiJ+LFCe\nGMX7XnOncy0iTZuUr4ud514RgE3VPjvvDYAxxr2rYoiIZCFrrXE6hmymfCQi4gfOl48CpefogO/1\ninSupZxLPQzhAN7hBLnO0/5Pa+3ZHLj3OdbagPsTHR0dcJ+dFfe91Htc7Psy2z4z7S7U5p577nHs\n/3V2/XHy+zc7P9/J7+EgoXzkwj/KR1l7D+WjnP2jfJT198hIoBRHW/AOHWiazrVrfa/rUp1bC4QC\n16RuaIyJAOqmaZud9z4nJiYGj8eT3iW/FRkZGXCfnRX3vdR7XOz7Mts+M+2c/H/pFKf/zoH0Pezx\neIiJibnsz/UTfp+PApHyUdbeQ/koZzn9dw7E7+EMOV0NX0SFtxXYncH1GcBZoE6qc/mBX4HtadrW\nxruCzydpzg8AkoE7curevmtWxJ9FR0c7HYL4Ad/POsfzyeX+CfR8FB0dbZcuXXrx/4NFXED5SDKy\ndOlSGx0dnWE+cvsmsHcBV/kOBwC5gFd9x3ustVNStU3ZNTwRGAOcwLtreC3gBmvt12nu/QbQH5iN\nd9fwGr7PWG6tbZWmbbbd29fe/nbsN64omN5ICRH383g8jj/ZEvczxmD9dM5RMOUjN/9eIHIhykeS\nGRnlI7cvyHAf0NL3dcpP66G+Vw9wLhlZa3cZY5oBLwPP4N0L4nugg7V2STr3fhTYA/QBbgAO490x\nfEjahtl57xRVulYhtlcsT9351PmaiIj4JY/HEwjDhoMmH4mIBDNX9xwFC2OMXbZnGd1ndufgEwcx\nxi8frEoQ05M6yQx/7jkKFsYYGx0dTWRkpP5Ni19SPpKMpDysi42NPW8+UnHkAsYYO2TIEF7Y9wIn\nxp0gT648TockIpJlMpOMxB00rE5EgkFGD+tUHLmAMcYuHdabAQWXs+WRH50OR0QkW6jnyP1UHIlI\nMMgoHwXKUt5+L98HU1j3+imYPBnO/mPbCRFXC4D5JCLiE4hbS0jw0PeuZCQzW0uo58gFjDF28ODn\naVW4CJGffw4HDsBzz8Gdd0Ku9PbuE3EXjfGWjGhYnf9Qz5H4O+UjyQwNq3O5fyQjjweGDoVff4Vn\nn4W771aRJCJ+T8Pq3E/FkYgEAw2r8zeRkbBkCUyaBNOmQdWqMGECnDnjdGQiIhLgNKxORAKVhtX5\niQs+qVuxwtuTtH07DBoE994LuXPnXIAiF6BhDJIZ6jlyP/Ucib9TPpLMUM+RH8jwSV2zZrBwIUyf\nDp9/DpUrw9tvQ0JCjsYoInIpMvOkTkRExA3Uc+QCF/2kbu1ab0/Sxo3w9NPQuzfk0d5IIuJu6jly\nP/UciUgwUM9RoGncGObOhc8+g6+/hkqV4LXX4NQppyMTEREREfFbKo78WcOG8MUXMG8eLFsGFSvC\nq6/CyZNORyZBRpO3RQKHFmQQf6bvXclIZoZ5qzhyictKRvXqwezZsGCBd/GGSpVg5EgVSSLiCppz\n5F9iYmI0oV1EAlJkZKRWq/MHWT7Ge8sWGDbM25v0+OPQrx/kz5919xcRuQSac+R+mnMkIsFAc46C\nzb//DTNmePdK2rDBO9zuxRfh+PFzTXYc2UHUh1Fc8eoVtJ/Snl///tXBgEVEREREnKfiKJDVquVd\n/tvjgW3bvMPthg/n5OEDtJ/Snluq38Kq+1cReVUkHaZ24EySNpmVS6Mx3iIi4gbKR3K5VBwFg5o1\nYepUWL4cfvqJ8Go1eWbRaQZUu4tyhcoxqPkgkpKT+Pmvn52OVEREHKYFGUQkUGVmDqzmHLlATo/x\n3rlqHhsHdufW3Xkw/fpx4qH7qTC5Phv6bqBsobI5FoeIBBfNOXI/zTkSkWCgOUd+ICef1FW+tiOf\nP9OFHk+UY8N3c0iqXJEPN1agrC2QI58vIsFFq9WJiIi/UM+RCzjxpC7ZJvPRpo/YcWQH1yYU48ZP\nNmPmzIEBA+DRR6FQoRyNR/ybx+PR0r9yQeo5cj/1HIm/Uz6SzMgoH4XldDDiDiEmhHvq3vPfEzcA\nP//sXQK8UiV45BEYOFBFkoiIiIgEDfUcuYDrntTt3OktkhYs+G+RVLCg01GJiJ9Tz5H7uS4fiYhk\nA805kotTpQpMnuxd3W7HDm9P0gsv/M8+SSIiIiIigUbFkZxftWrw0Ufw7bfw449QubJ3M9kTJ5yO\nTFxGy/6KBA4t5S3+TN+7khEt5e0n/GYYw/btMHQoLFoEjz0G/ftDAa1wJ5oAK5mjYXXu5zf5SOQ8\nlI8kMzLKRyqOXMDvktGPP3qLpMWL4fHHvUVS/vwA/J3wN//5+j9s+H0DlYpUYlS7UVxZ8EqHAxYR\nN1Bx5H5+l49ERC6B5hz5Ab8axlCjBkybBkuXwsaN3jlJr7yCPXGCztM7Y63lzevfpHqx6kR9GMXJ\nMyedjlhEHKR9jkRExF+o58gF/P5J3datMGwYSUuX8ELjBJ6btp/QAt7V7Zq814QXW71IVIUoh4OU\n7KRhDJIZ6jlyP7/PRxL0lI8kM9RzJNmrdm34+GOOzp1Jnb2nCalSBUaNIvlkHHFn4ggPDXc6QhER\nERGRC1LPkQsE0pO6u2ffTa5t23l+aRKF1//AlA5leGjiRsLyaeEGkWCnniP3C6R8JCJyPlqQweUC\nKRklJSfx1tq32PD7Bpr+mZf75vxG2Lrv4emnoU8fiIhwOkQRcYiKI/cLpHwkInI+Ko5cLuCT0fr1\nEBPjfX32Wbj/fsid2+moJAtpjLdkhooj9wv4fCQBT/lIMiNo5hwZYwYZY2YaY3YbY5KNMb9coH01\nY8xnxpi/jDFxxphvjDHprhxgjAkxxjxmjNlujDlljNlrjBlljMl7ufcOePXrwxdfwOzZMG8eVKkC\n48bBmTNORyYiImn41eqpIiIXIeg2gTXGJANHgPVAQ+CYtbbiedpWAtYCZ4DXgOPAA0Bt4Hpr7eI0\n7V8HBgCfAguAmr7jb4E2qR+1XcK9g+tJ3erV3p6k7dvh+efhnnsgV67/aZKUnMSBEwcokqcI+cPz\nOxOniGQp9Ry5X9DlIxEJSkEzrM4YU95au8f39VYgbwbF0QygC9DAWrvZdy4fsA1IsNZWT9W2FrAF\nmGWt7ZbqfH/gDeBOa+20S7m371pwJqOVKyE6GnbtgsGD4a67ICyMHUd2cNO0mzhx+gQnzpxgeNRw\nBl470OloReQyqThyv6DNRyISVIJmWF1KYXQhvkLlZsCTUrz43n8SmAhUNcY0SvWW232vr6W51QQg\nHuh5GfcOXk2bwtdfw4cfwkcfeTeXnTyZu2bczoDGAzjwxAG2PbyN0atGs2rfKqejlQxoCI6IiLiB\n8pFcroAqji5CHSAcSO837jW+14apzjUCkvAOlTvHWnsa2OS7fqn3lubNYckSGD8eO3Eikwev56Gf\nCkJSEuUKlaNjlY6sP7je6ShFRILCNROvYc/fe5wOQ0TEEcFaHJXxve5P51rKuSvStP/TWpt4nvbF\njDFhl3hvSREVhVm2jKE9SnL8tRFQuzanp05m1a8rKF+4vNPRSQa0MpBI4Ohaoyudp3dGw+vEHykf\nyeUK1uIoZYW50+lcS0jTJuXr9Nqm1/5i7y2pGcMDT02nao9DPHdTPn54rg/zRu6n4+ZTkJzsdHQi\nIgHvqaZPsevoLo6dPuZ0KCIiOS7swk0CUrzvNb3NdiLStEn5uth57hUB2FTtL/beAPTq1Yvy5csD\nULhwYerWrXvu6UfK+NlgOWYPvFt7HLkr5ybxkeLs/GgJPz/3HJEvvACxsXgKFABjXBOvjj1s3LiR\nRx991DXx6Ngdxx6Ph0mTJgGc+/km7rfjyA6stVopVPySR/scyWUKqNXqUstotTpjTBNgBTDcWjsk\nzbW2wEKgn7V2rO/cQqCV736JadqvACpba0teyr1957U60IVY690rKToawsJg6FC4/now/1xoxFqL\nSee8ZB8lI8kMrVbnfsYYW2pUKV5s9SL31rvX6XBELprykWRG0KxWdxG24B321jSda9f6XtelOrcW\nCAWuSd3QGBMB1E3T9mLvDWjTvQsyBjp1gvXrYdAg+M9//rvana+w/D3ud9pMbkP48HDKjC7D7B9n\nOxx08FAikox4MrHpXqDx503J590+V4WR+C3lI7lcQdlz5Ls+A7gFqJ9qL6L8ePciOpVmn6PaeFel\nm22tvTXV+QHA60BPa+3/Xcq9fdfUc3SxkpNhxgzvZrLFi8PQobTaO4xGZRoRGxXLxt83cvO0m1ly\nzxJql6jtdLQiQnD1HPn1puS9e8PYsd5eehGRABRMm8DeBVzlOxwA5AJe9R3vsdZOSdU2JWEkAmOA\nE3gTRi3gBmvt12nu/QbQH5iNNxnV8H3GcmttqzRtL/beKo4uVVISTJuGjY3Fk7SL5u8vJizS+0C0\n9xe9aVSmEX0b9nU4yMCnYQySGUFWHPnvpuTt2kHu3DB9OuTV+kHiX5SPJDOCaVjdfcBQ359iQKFU\nx/elbmit3QU0A1YDzwAj8RYxHdIWLz6PAk/iLXDeArrjTUQ3pm14CffWsLpLFRoKPXvCDz/waaN8\nJN1zF7RrR9LKFWw7vI3i+Yo7HaFI0AvGYXV+vSn5nDlQsCC0aQNHjmTmryEiEjACqufIX6nnKGtM\n3zqdJ+cNZMS+6rSetpq9ZQtRf/wXhDW+9sJvFpFsF0w9R6nl8AJBVay1JS7l3r7z3nyUnAzPPOMt\nlL78Eq66ChGRQJFRPtKAYpeIiYkhMjJSXcGX4bbat1GtaDVW7FvBivv70nnFn4R26QqNGkFsLFx9\ntdMhigQlj8ejnvHzy+pNyZsYY8KstWcv4d7/FRICr7wCZcrAddfB/Pnw739n9PcQEQkIgTaszm+l\nFEdyeeqVrkf/xv3pWu8OQvs/Aj//DJGR0KEDdOsGP/xwru3uo7tpOaklRUYUodGERmz6fZNzgfs5\n/eIrGYmMjAy6YXUXwd2bkj/6KIwcCa1bw7JlGTYVcQPlI7lc6jmSwJYnjze5P/AAvP02REVBmzYk\nPv8sHRd35f569zOr+yzm7ZhHx//ryLaHt1E4orDTUYtI8PCPTcmnT4du3fD06wctW7pik2Ed6zi9\n440bN7oqHh2749hzEZuSa86RCxhjbHR0NJEaVpf9TpyAN97g7JjRfFYxkVunb4KK3mkA171/HcNb\nDSeyfKSzMYoEGI9vWF1sbKzmHP3zmv9sSr5xI9xwg3evuf79L+Y/gYiIqwTTanV+S8PqckiBAvDc\ncxzdtIZdBRJJbtQQ+vQhYdcO9h3fp14jkWwQqWF1GXHdpuTnVbcuLF8Ob74Jzz57bgNuEZFAouJI\nglLxK6rw96DHaPZ0cZae2MLpq2vy3ld5uDpJS39fipSuaxG5ONbaOGAOEGmMqZNy3rdxeG9gh7X2\nu1Rv+Rjv0LlH09zqASAPMPUy7g1cYGuJChW8BdLixXDvvZCY3roQIs5RPpKMeDKxtYSG1bmAlvJ2\nhrWWuTvmsunQJmrZ4nT+/CfMpEnehP/001CixLm2C39eyPyd8ykcUZh+jftRIl+J8984CHm06Z5k\nQjAt5R3wm5KfPAndu3t7j2bOhHz5LvwekRygfCSZkVE+UnHkAppz5CIHDsCLL8K0adCnDzz1FJP2\nfsGQpUN45JpH2PXXLhbuWsjaB9ZSLO/55kSLSGrBOOfIGLMUaOk7TEm0KX93TzpFTHXgZd97woHv\ngRhr7ZJ07h2Ct+eoD1AeOIy3R2mItfYfCyxc5L0z/7AuMdH7c3LbNpg3D4qr511E/IOKI5dTz5EL\n7d0Lw4fDp5/yeqMkmr/6KfVrRAFw9+y7qV+6Po9em3ZUi4hkJJh6jvzVRT+ssxYGD4YZM7ybxVb8\nx5oTIiKukZmHdSqOXEDFkYvt2sWMHrXpuicvoY8/AY88wlMrYymSpwjPNn/W6ehcQ8MYJDNUHLnf\nJeejd96BF16AuXOhXr2sD0wkk5SPJDO0Wp3IpapUieXDHqDvf2pwdO03JJQvS+43x3JT2TZORyYi\n4h4PPwxvvAHt28OiRU5HIyJyyVQcuUSGqwOJo0a1G0WRuk24rt0+HuhfjoFn6vLvZl28T0rPnAHg\n9NnT9JvXj1KjSlHpjUpM3jTZ4ahzlp7SSUYyszqQuMcl56OuXeGTT+DOO73zNkUcoHwkGdFqdX5C\nw+r80Lp1MGQI/PADDBnCEyU38ePfOxl34zgOxh2k64yuTOo0idYVWzsdqYhraFid+2VJPtqyBTp2\nhMcfh8cey5rARESykIbViWS1hg1h/nyYOhWmTKH/vWMZfzySsvnL0PiKxvRr1I8FPy9wOsoco15P\nETnn3/+GFStgwgR46ilITnY6IgkiykdyuVQciVyOZs1gyRJG3lWRvBM+gDp1YNYsfj6yk8IRhZ2O\nTkTkomXJMO9y5bybxa5cCXfffW4IsoiIkzSszk9oWJ3/+3rX19wx63ZePtOSVh8s5eSZk5R7/UMK\ndukB3q5bRq0cxUebPyJXaC6eaPIEd/z7DqfDFslRGlbnflmej06dgttu877OmgUFCmTdvUVELpH2\nOXI5FUeBYdPvm5i3cx75wvJy354iFBj+ChQuDMOH83qezUzaNIlxN44j7kwc93x2D+NuHEfHKh2d\nDlskx6g4cr9syUdnz0K/ft65mvPnQ8mSWXt/EZGLpOLI5S560z3xD0lJMH06REezLvxPGD6chrf0\nB+Ddde/y3f7veK/Tew4HmTW0r4RkJDOb7ok7ZNvDOmth2DD48EPvZrFVqmT9Z4igfCSZowUZ/EBM\nTIz+MQea0FDvkrY//sg3TcpQo18M3HgjbNzIobhD5M2V93+an00+y/HTx52JVSQbRUZGainvYGeM\nd4XPZ56BFi3gu+/OXdpxZActPmhB8ZHFafFBC3Yc2eFgoCIS7FQciWS3XLmo9/zbVHvE8GUlON76\nOuo++hKP/+u/Q+reXfcuhV4uROnRpWk8oTG/Hf/NwYAvngp7kcCRrfvuPfAAjBvnXer7yy85lXiK\nDlM6cGvNW9ny0Ba61exG+yntiU+Mz57Pl4CnfCQZ0YIMfkJzjoLDd/u/4+NtH5PvDDyy1lB07CS4\n8Ua+f+BGOq0eyLJey6hYpCIxnhi+3fstS+5Z4nTIIllKc47cL8fy0cqV0KULvzzXj065PmHzQ5vP\nXbr63at5/+b3aVCmQfbHISJBSXOOXE7FUZD6+2949VVOvT6adc0r0Xz8l1CmDHFn4ig+sjinnjvl\ndISZpjHekhkqjtwvR/PRjz+S2K4NL9U5xuOfHiR/7gLEnYmj0huVWHHfCir/q3LOxCEBRflIMiOj\nfBSW08GIiE/hwjB0KAvblCX+hWhs7dqY3r1Z1/0arihwxblma35bw5TNUwgNCeWB+g9Qq0QtB4MW\nEckiNWqQa/Va7m9Why86VmL7oD7M272ALtW7qDASEceo58gF1HMU3JKSk7h15q3E/fITTy05TYMV\nv/D3A3dRaeibLD3yPT0+6cGTTZ8k4WwCb659k8V3L6ZOyTpOhy1y0dRz5H5O5CN79CiH2zXjj/wh\nbH91EF3r3oEx+jYRkeyjYXUup+JIkm0yi3Yv4kj8EZonXcGVr06EhQuZ2OZf5Bv4FLc3vg+AV1a8\nws4jO5lw8wSHIxa5eCqO3M+xrSUSEqBnTzh6FGbPhoIFc+6zRSRoZGZrCRVHLqB9jiRd27bx7T1R\nNPwtmTyxL8B99/Helsl4fvXwUZePOHnmJAt+XkDC2QTaVGxDqfylHAtVY7wlI9rnyH84+rAuKQkG\nDIDVq72bxZZy7mea+C/lI8kM9Ry5nHqO5Hwmrp/IwunDGbeqOBF79/N089N0eH4Szco3p/kHzSmZ\nryRF8hRh+d7lLLprkWPzkZSMJDPUc+R+juej1JvFfvUVVKrkXCzil5SPJDNUHLmc48lIXMtay7vr\n3uWDjR/QaEccsV+dpZiNYGqPmiyulYf3Or2PMYZ3vnuHeTvnMe+OeU6HLHJeKo7czzX5aNw4iI2F\nuXOhfn2noxGRAKPiyOVck4zE/ayFOXP47ZFehBUsQqk33oPISL7b/x195/Zlfd/17D22l29//ZbC\nEYVpX7k9YSFalFLcQcWR+7kqH82eDX37wrRp0Lq109GISABRceRyrkpG4hc+Wj+JH9+KZpjHQOUq\nPNvacqpOTbrV7MYtM26hVYVW7D66myIRRZh3xzxyhebK1ng0jEEyQ8WR+7kuH33zDdx6K7z1FnTv\n7nQ04geUjyQztM+RSIDpWe8eBvf8mcLlRtJr/W/EjshFwdaF6VKrNxNvm0in6p1ISk6i/ZT2fLT5\nI+6rd5/TIYuIn4iJiXHPAkEtWsCiRdCxIxw65F2wQUTkEqUsEJSRoO05MsYMAuoDDYDywK/W2goZ\ntK8GjABaAOHAeiDaWrs0nbYhwECgL3AVcBiYAQyx1san095dT+rEbyQlJ5FkkwhPSIQ33uDPF54j\nb7c7yTvsJbjySgYtGkTeXHnp06AP07ZOIzEpkc7VO1OlaBWnQ5cgpJ4j93NtPtqzB9q39/YiDR8O\n2gdJRC5DRvkoJKeDcZEXgEhgJ3AUOG82MMZUAlYC1+AtkJ4C8gMLjTHpDYQeA4wGtgL9gZnAI8Ac\no53tJAuFhoQSHhoO+fLBoEE8MqYdq+N/wl59NccG9GXB2qlUKVqFhhMasvH3jew9tpem7zdl3YF1\nTocuIpJ55cvD8uXw9dfwwANw9qzTEYlIgArmnqPy1to9vq+3AnmttRXP03YG0AVoYK3d7DuXD9gG\nJFhrq6dqWwvYAsyy1nZLdb4/8AZwp7V2Wpr7u/NJnfidP07+QfeZ3dn742qeXXqWO3ZGsLRzXVZ2\nacgLN78GwPsb3ueTHz5h/p3zs+xzNcZbMkM9R+7n+nwUF+ftPcqd27tQQ968TkckLqN8JJmhnqN0\npBRGF+Irgm4GPCmFke/9J4GJQFVjTKNUb7nd9/pamltNAOKBnpcas8iFlMhXAk8vD1uGHeHe1afJ\nu2Y9xXfu59neH8K770JiIlX+VYW/E/7mw40f0nhCYxpNaMSE7yc4HbqIyIXlzw9ffAEFCkC7dvDX\nX05HJCIBJmiLo4tQB+8co1XpXFvje22Y6lwjIAlYm7qhtfY0sMl3XSRb5QvPR2hIKFStyp5xI7i7\nV0Hipn1IYvWqLB7Rl6sKlCVmWQwvtX6JkW1HMnLlSCZvmnzJn6endCKSY8LDYfJkuOYa74INv/3m\ndETiIspHcrm0Wt2FlfG97k/nWsq5K9K0/9Nam3ie9k2MMWHWWg2YlhzRvVZ3Dt5ykBqlR3Htj3G8\nuugUpxd+xZ5nHqR1Re+UuZfbvMx7G96jYZmGfLf/O0oXKE3bim3RFDkRcaWQEBg9GkqXhmbN4Msv\noUYNp6MSkQCgnqMLSxnQfDqdawlp2qR8nV7b87UXyXYDrx3Ivsf2MXP8UcpuP8CXnWvTcOhE7+pP\nGzZwJP4IR+KPEDkpkq93f80TXz3B3Z/dTWbnHlxoWUwRkWzx5JMwbBhERcGq9AZ4SLBRPpLLpZ6j\nC0tZejt3Otci0rRJ+brYee4VgXdVvH8s592rVy/Kly8PQOHChalbt+65ruGUf+g61nGWHH/zDeFN\n76JG8Wd571AuTMumHCiRRFyXXHz1xAr+3v43Zwqd4fGfHmfR7kWE7fX+mIiKijrv/Tdu3Oiev5+O\nXXPs8XiYNGkSwLmfbyJZ7u67oVgxuPlmmDQJbrgBgGMJx3jyqydZd3AdVxW6ilfbv0rFIumuuyQi\nck7QrlaXWkar1RljmgArgOHW2iFprrUFFgL9rLVjfecWAq1890tM034FUNlaWzLNeXevDiQBafOh\nzXyw4QNynTpN/29Pk+/d9/lXn4GY55+HYsW4c9adxCfGs2TPEpKSk7i37r2M6TCGsBA9U5FLo9Xq\n3M+v89Hq1dC5M4wYgb37btp+1JarCl3Fgw0fZMkvSxi7biybHtxEoYhCTkcqIg7LKB/pt5wL24J3\nmFzTdK5d63tNvWnMWqAt3j2RlqecNMZEAHUBT3of4qodySUo1ClZhzEdxngPukDr0msYtXotdatX\n52Cf2/k64guKFr2SrQ9tJSIsgh6f9OClb19icMvBzgYufseTiR3JxT38Nh9dey0sXQrXX8/JfbtY\nH/Y9C3suJDQklEZXNGLhroWs3LeS66tc73SkIuKQzOQj9RyR6X2ObgHqp9rnKD/efY5OpdnnqDbe\nVelmW2tvTXV+APA60NNa+39p7u+/T+okYOw+uptbZ9zK6e1beWExNPnN8POj99DsubEQGsri3Yt5\nbslz1C+8O2NsAAAgAElEQVRdn6MJR7m+8vXcVecujDF4tK+EZIJ6jtwvIPLRb7+R1L4t7xT5mXu+\nPkzBPIVJtsk0mtCIkW1H0qpCK6cjlGykfCSZoZ6jdBhj7gKu8h0WB3IZY573He+x1k5J1XwQ0Br4\nyhgzBjgBPACUBm5IfV9r7VZjzNtAf2PMLGABUAMYgHevpP8pjETcomKRiqzvu574xHjyhOVhzJju\ndHl/PnxSD155hRURa9h0aBNtK7alyZVNGLFiBL/H/c5/mv3H6dBFRP7ryisJXb6S65vVZFVUFX4d\nE83i/d8SERbBdeWuczo6EXG5oO05MsYsBVr6DlP+I6RUkB5rbas07asDL/veEw58D8RYa5ekc+8Q\n4FGgD1AeOAx8DAyx1v5jMQZjjI2OjvbPYQwSsPYf30/T95rQ/2A5ek7dwo95TzK3dwtefdr7Lb/9\nz+1EfRhFdMto4s7E0bFKR2oWr+lw1OJGKcMYYmNj1XPkcgHRc+STHH+SX69vSvzxI8x96V76t3qG\nfOH5nA5LRFwgo56joC2O3CSQkpEEliPxR/j0x09JPnOaMtPn02Kyh0Kde8CwYWwKO0KD8Q3oUqML\npfOXZtrWacy4dQZRFaKcDltcSsPq3C/g8tHZs/Dgg7B5M8ybB8WLOx2RiLhARvlI+xy5RExMjCYs\ni+sUzVuUBxo8QN8m/akd8zb1H8vL9/YAZ/5dk2V3t6BB/qrM7DaTW/Lcwvgbx/PM4mecDllcyOPx\nEBMT43QYEozCwmDCBGjXDq67DvbscToiyWb6XUoul3qOXCDgntRJwNr0+yZil8WSa/9B7vtkN013\nnKLAS6PxVKxIyTpluHn6zewcsNPpMMWl1HPkfgGdj958E0aMgPnzoU4dp6ORbKIFGSQzNKzO5QI6\nGUnA+mrXV7zx9t1MX3UFEcfjGdq5CH81q89bHd8612blvpVs/3M7NYrVoEnZJg5GK26g4sj9Aj4f\nTZ8OjzwCs2ZB8+ZORyMiDtGwOj+gYXXib9pVakeXni9Qo9shetXfS//3tvL6O7/ATz8BMHTZUO6Y\ndQfLfl3GbbNu48VvX3Q4YnGKhtWJa9x2G0ydCl27wuefOx2NiLiQeo5cIOCf1EnA83g8RDZp4h22\n8vLLHO/eiXolZrP6qZ8onq84h+IOUePtGmx5aAtXFLzC6XDFIeo5Sp8xZhBQH2iAd4XTX621FTJo\nXw0YAbTAu3rqeiDaWrs0nbYhwECgL97tKw4DM8hg9dSgyEfr1sFNN8GwYdC7t9PRSBbSsDrJDPUc\niUj2y50bnnwSfviBhBNH+e7VExT/4GNITKRk/pKULVSWQycPOR2liBu9AEQCO4Gj/Hd7iX8wxlQC\nVgLX4C2QngLyAwuNMa3TecsYYDSwFegPzAQeAeYYY4K3UG3YEJYtgxdf9P4JhoJQRDJFPUcuEDRP\n6iRoHEs4RqfnKjJjVVmK/32GFY/fSrdjE9g5YCf5w/Nz/PRxpm6eyokzJ+hQuQN1SmpydDBQz1H6\njDHlrbV7fF9vBfJaayuep+0MoAvQwFq72XcuH7ANSLDWVk/VthawBZhlre2W6nx/4A3gTmvttDT3\nD658dOAAdOgAkZHw2msQomfGIsFAPUd+QHOOJJAUiijEy4/No3G3o3Stt4MrB73ED19XJf9vf/B3\nwt80ea8Ji39ZzO9xv9Nmchu+/PlLp0OWbKQ5RxlLKYwuxFcE3Yx3o/LNqd5/EpgIVDXGNEr1ltt9\nr6+ludUEIB7oeakxB4wyZeCbb2DjRrjzTjhzxumIRMRh6jlygaB7UicBJ6Mx3qcST5EnycCYMTB6\nNKtvuJrx7Yry/p0zAFiwcwGDFg9i44MbczBicYJ6ji4so54jY0wTYAUw3Fo7JM21tsBCoJ+1dqzv\n3EKgle9+iWnarwCqWGtLpDkfnPno1Cm44w6Ii4NPP4UCBc5dOpt8lt+O/0bhiMIUjijsYJCSGZpz\nJJmhniMRcUyeXHkgIgIGDYJNm4jYf4gxjy2E//s/sJaqRatyNOGo02GK+IMyvtf96VxLOZd6xZMy\nwJ9pC6NU7YsZY8KyMD7/lScPzJwJFSpAq1bwxx8A7PprF7XeqUXzD5pTdkxZXvr2JYcDFZHspuJI\nRC5bpp/SXXEFxya+zT09wjn58jDONGvCW+/3pUOlDgBsObSFjlM7Un9cfR798lFOJZ7KvqBF/E9e\n3+vpdK4lpGmT8nV6bc/XPriFhcG4cd45SNddB7/8Qs/ZPXmo4UPse2wfO/rvYML6CSz5ZYnTkUoG\n1Gskl0tPjFwiJiaGyMhI/aOWgNeyfEv2PTCGf1d4ji7L9xI7PJk891TnYP2faDutLTGRMTQo3YCX\nV7zMfV/cx7Su0y58U3E1j8ejOZVZI2Xp7dzpXItI0ybl62LnuVcE3lXx/rGcd69evShfvjwAhQsX\npm7duudyU8r/x4A9XrYMWrcmskQJaN6cvxr9TrVy3hGMpQuUpl5CPWbOm0mr/q3cEa+OdazjTB17\nPB4mTZoEcO7n2/lozpELBO0YbwkYnssZ4334MDzzDCfnfMr7t1dnwGsrwRjiE+MpMqIIcYPiyBWa\nK0vjFWdoztGF5fCco8rW2pJpzisfpfj4Y470vpOt70TT8q7BnEo8RdP3mzKkxRC61OjidHRyHpeV\njyRoaM6RiLhX8eLw3nusfPUxOn62zTve/4cfOJZwjBATQmhIKAlnE1j6y1KW/LKEhLMJF76nSGDa\ngneYXNN0rl3re12X6txaIBTvnkjnGGMigLpp2kpaPXpwaPyr1HoomiED61DznZrULVWXztU7Ox2Z\niGQj9Ry5gJ7UiUDcmTiajbuG/2wuyC2fbGXKNXn4/dHePNziCVpNbkXu0NwYYziVeIol9yyhWN7z\njRYSt1LP0YVlcp+jW4D6qfY5yo93n6NTafY5qg1sAmZba29NdX4A8DrQ01r7f2nub6Ojo4nUMO9z\nji5fRJ6uPTj05EOUe3IYwbx3roi/8/iGecfGxp43H6k4cgEVRyJeR08dZczqMZz6dRf9p/1MuV1H\nePve2myvdyVvXv8mAAO/HEhSchJv3/C2w9HKxVJxlD5jzF3AVb7DAUAu4FXf8R5r7ZRUbSvh7RFK\nBMYAJ4AHgFrADdbar9Pc+w2gPzAbWADU8H3Gcmttq3RiUT5Kz86d0L493H8/PPssqEAS8WsZ5SMV\nRy6gZCT+LtvGeM+dy6F7u3OmSWPKTpwBJUow56c5jF03lvl3zs/6z5NspeIofcaYpUBL32FKMkj5\n7+RJW8QYY6oDL/veEw58D8RYa/+xjJoxJgR4FOgDlAcOAx8DQ6y1/1iMQT1HGTh40LuSXYsW8Prr\nEKKZCW6kOUeSEfUc+QklI/F32ZmMhs1/htpjZ9F5zTGSXhxOt/zzqVWiNp2qdeKpr5/i0MlDRJWP\nYmTbkeQLz5ctMcjlyUwyEnfQw7oL+Ptv6NQJypSBDz+E8HCnI5I0VBxJZqjnyOWUjETOL+FsAnfM\nuoM/Vy7izU/isUX/RYEPp9Nk0W2MajeKeqXqMfzb4YSYEC377XLqOXI/5aNMSEiA22+HuDj49FMo\nUMDpiETkIqk4cjklI5GMWWv54+Qf2MQzlBz7EadfeZEZ3Wtx97urICSEk2dO8q9X/kX8s/GEhoQ6\nHa6ch4oj91M+yqSzZ+Hhh2HDBpg/37vqpoj4DS3lLSLZKmWjtexijKFk/pKUKlIW8+yzfD1pCA2X\n7sBGRsKOHRyOP0zu0NycSTrDhxs/5NVVr7L+4PpsjUkkUMXExGT7v2m/FxYG48Z55yA1awZ79jgd\nkfjoe1cy4vF4iImJybCNeo5cQE/qxN/l9Bjvk2dO0mzCtTy5PoKus35gVNt8hPUfyNxd88mXKx81\nitVg+rbpvN3xbW6teeuFbyg5Qj1H7qd8dAnefBNGjPD2INWpc+704ZOHWX9wPcXyFqN+6fpaAjyH\naM6RZIaG1bmckpHIxTuWcIy31r5F8k/b6ffWGhLyhPPkbUWYOvAbjDGs/m013Wd2Z+9je50OVXxU\nHLmf8tEl+vhjeOQR+OQTaN6cVftW0fnjztQuUZvdR3fTqnwrJt48UQWSiEtklI/CcjoYEZGsUCii\nEM+1eA5aAPeeZUW/mxkX48EUmwp33knN4jU5cuoIAKfPniY0JJSwEP3IE7mQmJgYrZ56sXr0gKJF\noWtXmDCBe/c9zbgbx9G5emfiE+Np+l5TvvjpCzpV7+R0pCJBLWX11Iyo58gF9KRO/J0bhjFsOLiB\nJ19pzZz5hclVqw5Pdy/CTvsneXPl5bPtn2EwDLxmIC+3eVlPbx2iniP3Uz66TOvWwU030eeaw7w+\n8wR5cuUBYMD8AVQsUpHHmjzmcICBzw35SNxPCzL4AU2AFbk89UrXo1+fidTtfZYJv8/jmQEf02Jf\nKEnJSRx75hi/Pf4bi39ZzPsb3nc61KCTmQmwIgGhYUNYtoyY5aFsGNgdrOXgiYPM3TmXeqXrOR2d\niGSCeo5cQE/qRLLBF19wpGdXTj54H+VeegdCQ5nw/QRW/raSDzp94HR0QUk9R+6nfJQ1dm9dTmL7\n1nxbKZzH21uebfk8z1z3jNNhiYhPli3IYIxZClzOT00DWGttq8u4R8BRMhLJHve93Y6hE3dzZZFy\n2ClT6LMumqJ5i1K9WHUOxR2i+VXNaVq2qdNhBg0VR+5njLHR0dGac5QFzv71J2dv7EhI2asI/2gq\nhIc7HZJI0EuZcxQbG5tlxZEHb3F0OcnNWmujLuP9AUfFkfg7t47x3v7ndtp8EMWotUVos2g3j91b\nil9qX0lEWAR1S9Vl2tZpDI8azr317nU61KCg4sj9lI+y2KlTcPvt3tdZsyB/fqcjCnhuzUfiLlrK\n2+WUjMTfuTkZ/XHyDxbtXkS5ldto8OxbjG9XlAFTdxISEsoPh3+g2fvN+Os/f2mRhhyg4sj9lI+y\nwdmz8NBDsGkTzJsHxYs7HVFAc3M+EvdQceRySkYiOWPa5y/Q9LFXuapJBxg/nsSIcPK+mJdTz53S\nMt85QMWR+ykfZRNrYfBgmDkTFi6E8uWdjkgkqGm1uixgjBlkjJlpjNltjEk2xvxygfbVjDGfGWP+\nMsbEGWO+McZoOKGIg+o3u5Xm9xsOnDrM2SbXMHxKH5pc2YRO0ztR9c2qdJ3RlYMnDjodpogEGmNg\n+HDo3x+aN4ctW5yOSETOQ8VR5r0ARAI7gaNksDCFMaYSsBK4BhgBPAXkBxYaY1pne6QiOcxflqGv\nVqwa43tMITLyF6LL/syAx6ZRaMtOWl7Vks9v+5xqRatx/dTrOZt81ulQRSQQDRgAI0dCmzawfPn/\nXNr4+0YmbZzEN79+41BwgcFf8pG4l4qjzKtorS1urW0PXOjR8ktAQaC9tXaEtXYs0Bw4ALydzXGK\nSAY6VO7Ajkd28sK8BA6PiGbyhD/5z+Gq1ChegxdavcCx08fY9dcup8OUDBhj3jLGNHc6jkClffey\n2W23wZQpcMst8MUXAIxbN46OUzuy+JfF3P/F/Tz2pTaLFckOmdl3T3OOLoExZiuQ11pbMZ1r+YAj\nwLfW2rZprj0PDAWusdZ+l+q8xniLOGDDwQ3Ejr6J2dMs5vHHiRvQl/Kvlef7Pt9TIHcBikQU0UIN\nWSir5hwZY6oDNwN5gVnWWo1RyiLKRzlo3Tq46SYSYp6n+F/PsOnBTVQsUpHjp49T+53afHH7F9Qt\nVdfpKEUCUkb5SDOQs14dIBxYlc61Nb7XhsB36VwXkRx0damrMY0aceeVf/LG6yOZv/hV6vWoS91x\ndUlKTqJ4vuLM7jGbOiXrOB2qpGKt3Q5sdzoOkcvSsCEsW0ZY29Y8WzuEis9UAKBg7oJUL1adgycO\nqjgScYCG1WW9Mr7X/elcSzl3RQ7FIpIj/HUITogJYWa3mTRrehtjXu5Eu92GHuOW89UdX3J80HFi\nWsbQaXonkpKTnA5VRAJR1aqwYgVdNpzih57tsUlJfPvrt2z4fQNXl7ra6ej8kr/mI3GPDIsjjeu+\nJHl9r6fTuZaQpo2IOCwsJIx+jfvxQo9xrPtoBE3+ykujwWPh7Fnuuvou4hPjOXTykNNhSirKTRJI\nwq4sR5JnKSfXLGdGvXBun9aVqbdMpUyBMhd+s4hkuQsNq3sLuNm3wprGdWdOvO81dzrXItK0OadX\nr16U9+17ULhwYerWrXtuE7OUpyA61rGbj1O4JZ5LOS59RXWubxHK+99vps3tt/PjG0M48dMJtq7d\nSpnWZRyPz9+OPR4PkyZNAjj38y2LKDdJQKlVtRlsOUK9Ht3pvuQ0pl8Tp0PyWyk/i0QuVZYtyGCM\nKW2tPZjquKK1dneW3NxlLrAgQxNgBTDcWjskzbW2wEKgn28Fu5TzmgAr4hJPLHyCeVs/ZfpHp/gx\n1zFWx/Rm2b5viU+M59aatzI0aqg2jL1ETmwCG0y5KSsYY2x0dDSRkZH6JdMJZ8/CQw/Bxo0wfz4U\nL+50RCIBxePx4PF4iI2Nzd5NYI0xjwP7jTFVU53+lzFmjDEmV1Z8hh/ZgndIXdN0rl3re12X9oKW\nThV/Fkjfu6Pbj+b97lPZOf5lWideSa0RHzC67Sg+7fEpK/atIMYT43SIfseTiaVTs4Ny06WJiYlR\nYeSUsDAYPx46dIBmzWDPHqcj8juBlI8k60VGRl4wH2XVggzNgQeA3caYXgDW2nXA13j3/Aka1to4\nYA4QaYw5t8SVMSY/0BvYkXoZ7xRKRiLu0bRsU7o17sUbz7ej84ECtJ62mtolavN6h9eZvX220+H5\nncwko2yi3CT+xxgYNsy7Yex118HmzU5HJBJUsqo42mStfc9aexboknLSWjsfuDGLPsNRxpi7jDHP\n+/YqKg4UTjk2xvRM03wQcAz4yhjztDHmYeBboDQwIL37q+dI/FmgFvYhRYvyemxHmDQJ3n2X/cf3\nkz88v9Nh+R2neo4IgtwkAWzAABg9Gtq0gW++AcBay9jvxtJhSge6z+zOhoMbHA7SfQI1H0nOyZI5\nR8aYF621z/q+/spa2873dRhw0Frr94NmjTFLgZa+w5T/aCljFT3W2lZp2lcHXva9Jxz4Hoix1i5J\n596acyTiQgdOHOCaidfQq2BLnn7qc3r2CKNz3zEUiShCiXwluPbKa7VJ7EXI6TlHwZCbsprykQst\nWgR33AHjxzOi2E/839b/Y3jUcPYe20vMshiW37ucasWqOR2liF/JKB9lVXE0He+GfEWAPsBIYDrQ\nH2hhra192R8SwJSMxN95PJ6AfVp38MRBJq6fSJnV27ht1AKa9gmlbK2m/HTkJ64rdx3v3/y+CqRM\ncqA4Um66SMpHLvX993DTTTzX4izdX/v63B5IT371JAXCCxAdGe1wgO4RyPlIsk5G+SirllwaDizF\n25PSE5gCPO+79nAWfUZAS5lzpH/QIu5SukBpBrccDC3hpeVF+GZuCQoN/oRToZZr37uW+Tvnc0PV\nG5wO09VSVgdygHKTBIYGDWDZMvo2qUnYmHHw8ttgDEnJSXo4I5LFsnIp71xALmttvDGmJtAR+M5a\nuyxLPiCA6UmdiPsl22RyDQ0jccsthBQqBBMn0nfug1xd6moebqTfszPDoaW8lZsugvKRu02YO5Tr\nHnyR0OYt+OzBKEatHcPq3qupWOR/dxaZ89McPt72MXnC8jDw2oHULqFOUpHUMspHWbUgA9baRGtt\nvO/rH6y1o5R8RCRQhJgQ6pWuz9iHG8Hq1Rye8Drzds6jQekGTocmGVBukkDS+4bBbJz+Gmc2byDy\n2fF8c/vX/yiMpm6eSr/5/WhVoRWV/1WZqA+j+PHwjw5FLOJ/sqw4ksuj1erEnwXL9+70W6fz1o+T\naNv6APbxxxhaox/XXHmN02G5noOr1YkEFGMMt1/3ILW/30fj0g2ofs/jcPz4/7R5fc3rvN/pfe6r\ndx9PX/c0DzV8iPc2vOdQxDkvWPKRZJ+LmnPkW7HtcvrbDWDTruwm6BcHET9Q+V+V2fbwNg7FHaJI\n7tG0ee1zav0yhVOJp7it9m0MjRpKWEhWTeUMHCnzKWNjY50ORTJBc2D9QEQEfPyxd7nvli1hwQIo\nVQqAJJtERFjEf5uGRXDi9AmnIhVxlczMgb2oOUfGGA/e4uhyxoxba23UZbw/4GiMt4j/WbxtLldF\ndeHM0CGYrrfSd25fospHERulAuB8nJhzJBdH+cjPWAvDh3v3Ylu4ECpX5s01bzJ23VhGtxvN4fjD\nPPHVEyy4cwENyzR0OloR18j2pbzl8igZififRxY8QtO9lttiZ8GWLaw7/Qu9v+jNxgc3Oh2aa6k4\ncj/lIz81fjzExMCcOdj69ZmwfsK5BRmeavoULcu3vOAtRIJJjizIIJdHc47EnwXj927+8PysuSoM\nunWDZ55h37F95A/P73RYrqQ5RyLZrE8feOcduP56zKJF9GnQh8V3L2buHXODrjAKxnwkWUs9Ry6g\nJ3Xi74Jx0739x/dz7XvX0q10W4Y9+DGt++RmaN/ptKvUzunQXEs9R+6nfOTnvvnG+8Dm9dfhttvO\n2+z46eMs27OM0JBQospHkSdXnhwMMnsFYz6Si6dhdS6nZCTinw6cOMAHGz6g8YdfU/9oHop+usDp\nkFxNxZH7KR8FgC1boGNHePJJGDjwH5f3HdtH5IeRVChcgYSzCZw4cwLPPR6K5CniQLAizlBx5HJK\nRiJ+7sQJqFwZFi+G2tps8XxUHLmf8lGA+PVXaN8eOneGl14C899/dnfNvouKhSsSGxWLtZYH5z5I\nwdwFGdlupIMBi+QszTnyA5pzJP4s6L93CxSAp56C6GinI3ElzTkSyWFXXQXLl4PHA716QWLiuUt7\nj+09Nw/JGEPL8i3Ze3yvM3Fmg6DPR3LZVBy5RMq+EiLipx5+mLMrvmXNF2M5cOKA09G4SmRkpIoj\nkZxWrJi3N/vPP6FTJzh5EoDGZRrzznfvcPrsaeLOxDFx/UQal2nscLAi7qHiSEQumwp7ePH71xjc\n5BTJQ56nztg6zPlpjtMhSYAyxgwyxsw0xuw2xiQbY365QPtqxpjPjDF/GWPijDHfGGO032AwyJcP\nPvsMSpaEVq3gzz8ZGjWUZJtM0VeKUmJkCSoUrsCj1z7qdKRZRvlILpfmHLmAxniL+LfNhzbTYUoH\nNvRaTcnqDdj0+Xiivr2f35/8nfDQcKfDcw3NOcoaxphk4AiwHmgIHLPWVjxP20rAWuAM8BpwHHgA\nqA1cb61dnKa98lEgshaeew5mzfJuFlu+PH8n/E2oCaVA7gJORyeS4zTnSESyVbCP8d59dDcNyjSg\nZNFycPvtXL1wI7lCc/Fn/J9OhyaBqaK1tri1tj1w8AJtXwIKAu2ttSOstWOB5sAB4O1sjlPcwhh4\n8UXo3x+uuw42baJwROF0C6O4M3GMXjmap79+mgU7/W8FzmDPR3L5VByJiFymmsVrsua3New4sgPu\nu4/4iWMJsyGUyFfC6dAkAFlr92SmnTEmH3Az4LHWbk71/pPARKCqMaZRtgQp7jRgALz6KrRt612s\nIY34xHhafNCCNfvXUCiiEP3m9+PNNW/mfJwiDgpzOgDxSlmQQWNlxR8F+/dt1aJVeaXtKzSe0Jii\neYvyGUf56opRhIXoRyx4n+Tqaa4j6gDhwKp0rq3xvTYEvsuxiMR53btD8eLe17ff9m4a6/PZ9s8o\nlrcYH9/6McYYutXsRqMJjejfuD/G+MeI2GDPR3L5lLldQis5ifi3XnV70aV6Fw6dPET55LmEz10D\ndzgdlTukPPiJjY11OpRgU8b3uj+daynnrsihWMRNoqLgq6/gxhvh99+9PUp4h9RdWfDKc4XQFQWv\n4NTZU1gsBv8ojkQul4bVichlU6+AV6GIQlQtWpXwu+6B+fPh6FGnQ5Lgltf3ejqdawlp2kiwqVvX\nuxfS22/DM8+AtbSu0Jo5O+bwyQ+f8PNfP9NnTh9urnYzIcZ/fl1UPpLLpZ4jEZGsVrQotGsH06fD\nQw85HY0Er3jfa+50rkWkaXNOr169KF++PACFCxembt2654YqpfziqeMAOd6zB0aMIPLll+Gee9h3\n990MKTeEkStH8sfJP6h+ojoDGnt7lVwRbyaON27c6Kp4dOyOY4/Hw6RJkwDO/Xw7Hy3l7QJaOlUk\nAC1YANHRsHat05G4hpbyznrGmK1A3vSW8jbGNAFWAMOttUPSXGsLLAT6+VawSzmvfBSM4uOhRw84\ncwY++QQKpL+8d3xiPCv2rsAYw3XlriMiLCLddiJup6W8RURy2LKquYn/+Ue+WTYZ/bIpDtmCd0hd\n03SuXet7XZdz4Yhr5c0Ls2dDuXLe+UiHDv2jyaG4QzQc35CYZTE8v+R5rpl4DX+d+suBYEWyl4oj\nEblsKV3X4jX8m+HcO7c3W/9diq/G/YcnvnrC6ZAkCFlr44A5QKQxpk7KeWNMfqA3sMNa+4+V6mJi\nYvRvOhiFhcH48d5FGpo1g59//p/Lg5cOpmOVjqy4bwWr7l9Fs7LNGLpsqEPBnp++dyUjHo/ngoug\nac6RS2gpb5HAcCjuEKNWjuKn/j9RssA86i6Yx1Vbp9G3QV+qFavmdHiO8Ggp7yxljLkLuMp3WBzI\nZYx53ne8x1o7JVXzQUBr4CtjzBjgBPAAUBq4Ib37a/XUIGYMxMRAmTLQvDl8/jk0bgzAnr/38Ni1\nj/maGVpVaMWUzVMyuJmI+6T8rp3R6qmac+QCGuMtEji2/bGNrjO6sr3/dti7Fxo0oNmLlXm53Ss0\nv6q50+E5SnOOsoYxZinQ0neYkjxS/rt6rLWt0rSvDrzse0848D0QY61dks69bXR0tB7WCcyZA/fd\nBx98ADfeyKBFg9jx1w6mdZ1Gsk3mlo9vocmVTRjccrDTkYpkWsrDutjY2PPmIxVHLqDiSCRwJJxN\noNpb1YhpGUPPOj1JqFyezp0T+OTFnymSp4jT4TlKxZH7KR/J/1izBjp1gmHDONWrJ3d8egeePR6S\nbRhKTaIAACAASURBVDLXV76eyV0mEx4a7nSUIhcto3yk4sgFlIzE33k8Hj1lTmXrH1u5Y9YdbDu8\njcmLCtC8xV2UG/6m02E5TsWR+ykfyT/s3AkdOkDPntjoaP48dQRjDMXyFjvX5I+TfzB502TiE+O5\nudrN1C1V17FwlY8kM7RanYhIDqpdojabH9pM4uBE7nx8EuW+2+F0SCKZpgUZ5H9UqQIrV8L8+Zje\nvSkeXvh/CqPf436n8YTG/HD4B06eOUnbj9qyePdiBwMWOb/MLMigniMX0JM6kQD2999QtiwcPgwR\nwb0niHqO3E/5SM4rLs67F1JyMsycCfnzAzB4yWCOJhzlrY5vAfDpj58yetVoVty3wsloxQXW7l/L\nin0rWL1vNX+c/IN84fmIKh/F1aWuplWFVoQY5/poMspHWq1ORCQ7FS4MtWrBqlXe/UNERPxR/vze\n1esefBBatoR586BUKY6fPk6FwhXONatQuALHEo45GKhklS9//pJNv2/i74S/WfXbKk4nnaZswbIU\nzVOUUvlLUTJ/SfKE5aFQ7kIk2ST+ledfHDl1hCW/LGH53uXsOrqLsgXLcjDuIGUKlOHX/b+y+JfF\n5MuVj/DQcG6sciOxUbGU/v/27ju6iqIP4/h3kpCEUAOhd0R6kw4iRAGRItIRFClSxQJibxRfBBRU\nsNMEQRQQQVSkKAQEVJBeRJEuhE4C6SSZ9497wRATSEi5Kc/nnD1hd2d3Zy+b+8vstDzFXH2r11HN\nURoxxrwI1AHqAmWBo9bacomk1Zs6ydTUxvsmXn0VYmLgjTeSd9yhQ+DuDmXK3DxtKrgQfoF9Z/dR\nNHdRKhSokOrnV81Rxqd4JDdlLfzvfzBrFixfzirP4zy67FEWdl1I4VyFGfTdIBoUb8D4luNdkj3F\no1sXGxvL+1veZ8m6TyhwOJAoN6hcrj7vhv1EszLNiY6NZmvgVpqUbELA0QCK5C5CaGQo4THhlMhT\ngsNBh8nnlY/QK6EQE8uAIm3Y9+u3fDF4NZWWt+Gl1v/jk62f4OXhRWhUKCFRIQRFBFHetzzDGw6n\nY5WOlMxbMl3uVX2OXGMc4A8cAC7y73CrIpLdtGwJP/6Y9PQ7djiar9Su7RgpKiYm7fIGWGv5Zv83\nVHq/EkO/G0qVD6qQ641cPPDlA5wJPZOm15aMR32O5IaMcbzwGT0a/P2594Q34+4ZR5+lffCf40+t\nIrUYe7djctiY2BiOBB3hfNh51+ZZbmjR3kUUn1wcr3FezJv/At+/fYoRK4IZtTyUkRPX8/c3ZXHb\nvYc/zv3B4u6L2XxyM++1eY+QqBBqFq1J69tacyb0DKOaj6Ls8csc+bs9pyZE896Yzby+BnI93JdT\nE6Lp/dRM7lt9hKcr9iU4IpiW5VtSLn85Ai8H8vSqp7ltym3Um1aPj7Z8RFq9pFGfIxcyxpS11h5x\n/nsP4GOtLZ9IWr2pE8nKIiPBzw9OnIC8eW+c9uJFKFcOXnsNBg6Edu3gkUdgwIA0yVpIVAhdF3Zl\n9aHVeLl7YbFMaT2FcT+Po1tQCaLPnebtt/fi5pXy/lKqOcr4FI8kWX78EXr1gqlT4cEHr9t14tIJ\n2s1vx7mwc1yKvMTguoN5s9WbGKOvgIwgKCKI0QGj+f3k72w/tZ0fev3AoplPM3bqbj7oeRszKocx\n9u6xjF/7OpOPVabhjBWsqexJvrFv0Xnbc6zvu557591Lz+o9iYqOJGjlN3y8pwxh27dQ9KmXqRg2\nkXYthvLl3i/x8vDiwoWTDD5XloY/H6LTUR++KRNOyYcfo3vQNJ68fxwvr3mZ+yvez3cHvsPduOPn\n48fElhPpWaNnmty/hvJ2MRWORIQmTWD8eEdbfafQqFAuRV6iaO6i//7B8PnnsGABLFvmWN+6Fdq3\nhz//vHnBKpmCI4Lp8EUHgiKC2HtmL58+8ClPr36aF2o/we3PjueOo5EcyWupFpQD79Vr8WlwZ4qu\np8JRxqd4JMm2a5fjO+qJJ+CZZxw1S0D7+e2pU6wOY/zHcDHiIv6z/Xmt+Wt0rdrVxRmWKzFXuOvT\nu6hRuAbWWtYfW889Qb5MfGsHP495lG5hn1IqbykqFawEBgKOBJAv3DL8V+i3MZSd1fzYUxhO5Iym\nyhnL3XtDiTIx/NS1LhPLnaBogTL8ef5PvD28CQwJxNfLlwoFK3Ak6AgXwi+QLzyW7vvc6H3Cjyp7\nT+Pt7sWh3FHkLVyKHVHHKFmjKVPct/B9ZXeaVbmPT+7/5LoRElODmtWJSJpSE5wkqFvXUdBxmrBh\nAoUnFabmxzWpM60Ox4OPO3Z8842jKV3c49q0gXHjUjU7Jy+fpO60uuw5s4eqhauCgW8PfItPDh9K\njJ5MVOwVdq75klYDvXiv522EdLgPLlxI1TyISBZQs6ZjqO+5c+HxxyE6GoDtp7YzoM4AjDEUyFmA\nrlW7sj1we5pnR/Ho5nac2sGlyEtMu38a/uX8KZmrGEOmb2f/s315OHI+nu6etCzfkhUHV/DToZ8o\nkacEd9Zox44hHXnkrTtZf1sOTGQkFQ8HczDHZTr3sNR4MgePF9tOmFsMe8/upUqhKnSt2pXzz57n\nwgsX2DxwM2eePcOpkacYft8YfmpZntZdQin+khe1n8nNq/3K8dq9nixokItD+WLpuTWKI5Oiqffl\nOu6c1oig8KB0+3w0Wp2ISHqoWxdWrQJg9cHVzNg2gwNPHKBY7mI8vfJpWnzWguF3DGXIqpW4veeY\nMDbsShgxsTHkGTcOKleGV16BPHlSlI3Ay4F8sOUD3tr4Fu5u7hhjyO2Zmyn3TWH4yuG0/yOWhvti\neejVahxZMxx3N3eW1fEh/+5w7unYhBrr96f4o5CMbfTo0fj7+6tTuyRdyZLw88/QrRt07Ahffkl5\n3/Ks/HslA+sOJComijWH19C7Zm9X51Rw1po4u8J3r9adU5PHEGqi+bDyJdwPu9O0VFPyeuVlz9A9\nVPKr9N8TPOb4ERMbg5tx45WYKCKiI/Dy8OJs6FmK5SmGh1vCRQy/XH682vxVXm3+KqsOruLT7Z+y\n9M+lHDAXcPdx54H6D9B9/xKe/nAkaw8eoc3ExTTbdp7GQbezcOgaahSpkaJ7DwgIuGkBWs3q0oGa\n1YkIu3dju3Zl2TdvMm3rNHJ75mZBtwVsD9xOq7mtuBx1mbGRTWg+92fGj2/Hycsn2X1mN27GjTa3\nt2HhJxdxf/Kp62uVkuFC+AXe3/w+kzZNIjw6HHfjTuvbWnNvhXsZvmI4tYrU4sDxnfz1HowZXJlN\n5dw5cP4AAX0DmLZ1Gl5XYnmt32zOLPmc6i0evPkFE6BmdRmf4pGkyJUr8Nhj8Pvv7J89iZY/9aFC\ngQqcvHyS6oWrs7DbwkT/aJb0Ex0bzV2f3kVVv6p09WtGk7aD6f94Ke66f9i1eYjS0+XIy2w6volJ\nv0xi/dH1XIm5wmP1HuPD3z9kTodPMa+8QoPfA2n7iDsv9vyQR+s8muJrqs+RiyWlcNSnTx/Kli0L\nQP78+aldu/a1t3ZXS7ha17rWM+d6TGwM1etUJn/RstR5tCS58hdiV85dvNj0RVb+tJICOQtwvMBx\nJn4VzNfhF1hRLz9hJcKoWqgqPXP35LnVz/HY0TAaRvoRMuItSucrnazr/37idyacmEBQRBAexzwo\nnqc44/qP46GvH2JsubGM3zCeavWrMWzNZUoGFibipecYfXQ02wK3UfJ8Sdzc3KhUtxJN560nX2Be\nqr7zOS3uaXHT6wcEBDB79mwAypYty5gxY1Q4yuBUOJIUsxYmToQPP+TSV/PZ7BdBHs881C9RHzfj\nxtu/vM34DeOJiI6gZ/WevNfmPbw8vFyd62wnOCKYMevG0OaNRXgXKkqd+QHk8szl6myx+cRmHlny\nCBcjLlI8d3GalGrCx1s/5pcLXSgw7ytaD/Cmeq1WzOk0h/ze+W/5OiocuZhqjiSrC9C8EonadHwT\nXRd2JexKGD98EEzsuHE0evh5HvjyAX74+wfyeObBzbgxvd0n3Hnng6ybNYoX/vmU0c1HM3b9WEKi\nQuhTqw///LKSjz8+Qe3n8rDv8T/w9kja6HHRsdEUm1yM8S3GM+W3KYxvMZ4HvnyA9X3X02dpH4Ij\ngjkffp5hFXoxdfhKzMaNULEiUTFRlHm3DD2r9eSzXZ/xZIMnmRXwNjsmhfDJJ4N4/qEPk/1ZqOYo\n41M8klSzYIFjkIa5c6F1awC+2vcVL/30Et/1+o4COQvwyJJHqFaoGo1KNmLKb1OIsTH0q92PAXVu\nfXROxaPEHQ06yqPLHmXHqR2U9y3PvFLDqdj/Gdi/P9UH/EmJyOhIBn47kAV7F1AgZwGG1hvK0aCj\nFJv6KUMO5KNpn2gi8uVibqe5tLqt1S1dQwMyZAKaV0Ik6wmNCqXzgs5Mv386i7sv5uhtBVmx6A0u\nhF9gWc9lFMlVhGalm1HJrxJl/z7HZW83xpz6kmqFqrHh2AYMjo7Meb3yYqtWIZebF6UCw5i1bRZR\nMVFJysO5sHMAdK3alZOXT1IkVxHK+5an7edtORZ8jFhiWfrgUt7bVgTTqRNUrAiAp7sny3stZ9aO\nWURERzD518m83uk9zIM9if18XrI+h4AkzCshIllMjx6wZAn07Qsffww4+ls+2fBJKhasiJ+PH2P8\nx/D1/q95asVTPHfnc4zxH8ObG99k9o7ZLs16VhQdG027+e24u+zd7H1sL8PqD+PUk/0IfeGZDFUw\nAvDy8OKzTp/xWL3HCIkKYfEfi5m1YxYFX5/EtxVi+GZ2BPki4IEvH2DAsgGEXQlL1eurcJRBXO0A\nK5IZ6dlN2OGgw/jm9KVdxXbULlqbX4tcocEpN7YFbuODzR+QyzMXC7otoHHJxgS8N5JvKsaS1ysv\nQ+oOYemfSzl5+SQHLxxk+rbpPFKrD/NLXqDernN8su0TWnzWgvAr4Te8fkhUCK+seYWcpy8w8qV6\nTPPsSpvZLTl88TAxNoanGz/N9sHb6XC5OMyb55j1Po47it3BS01fon3F9gSODCS/d37eKXOS+7eH\nXit0JYW/v78KR5mIXtZJqrnzTtiwAd59F0aMoJBXAfad3Xdt976z+wiLCmO0/2jaV2xPy/Itebv1\n28zblbwXMHEpHiXsWPAxgiODeemulyiSuwh9zpWg3CV3fmtdzdVZS9Q7973DlPumcPjCYdxwY+Px\nTQy/O4Ko+nWZM+0MDXxuZ/aO2ZR4uwRvbnyTWBt703Mm5WWdmtWlEWNMb6CMc/UJIAfwtnP9iLV2\nXpy0dtSoUfhrdCCRLOV82Hlum3ob2wdvp5xvOTZ8+yF+/YdR88kc3FHsDuZ2mkvFgo6aGpo14/Kz\nT/EMq9hzdg/lfcvT+rbWzNs1jxOXThAcGUyrnZd5bk9+bt9ykC4Lu9C0VFNGNhmZ6PW7LOxCuTNX\n+N+ra9mWLxzPKzF4x7qx/sWH6DlwCr45fR0dqOvVc8xP0vu/I0kFXg6k3vR6VC9cne2B23GzsGtC\nMH2GFOGL13Ylqc13gHN0IPU5yvjUrE7SxMWL0LUrkV4eNPA/wO1l61AgZwGW7F9C45KNaV6m+bXv\nss93fc78PfP5vtf3Ls501nI+7DzlppTj8FOHKZizALH16vFUzX/oO245dYvXdXX2bmh74Ha6LOxC\ncGQwtYrUYu+p3Yz67jJ3H4yhXS9DUKHcBEUE4enuSZ9afXi+6fOU902wJ8s16nPkAsaYtcDV2R6v\nfshX/xMCrLX3xEmrYCSZmtp4J+6jLR8xZt0YGpdqzNZjv/H3a+fwPHP++iG5o6KgQAE4efI/zRus\ntSzdv5SB3w5kYv2X6N9hFOb0ad7a/gH/XP6HKfdNSfC6UTFR+I3JRdDSKrgNG8blvr3ovrAbTx0p\nwn0frYaGDaFtW3j/fShfHr7++trkjfEdDTpK1Q+rck+5e+hVvRc9Z//OV0d/4PwLTzG43uAkfxbq\nc5TxKR5JmrlyBYYNI/rXTSwc/zDn/XLRrmI7giOCuXfevYxoNAJvD28mbpzI/M7z8c3py4ZjGyia\nuyidq3RO8ih3ikeJe+mnl1i6fymvnq7MHXNWMXZSez7v+sW/k5BnYGFXwhi4bCAL9y0EYGjdIRT+\n6DMGrwvltV7F+LGaNwcvHMTduGOx3Fn6TsbdM44mpZrgZv7bUE6FowxOwUgyOwWjG9t7Zi/7zu7j\n9oK3U7vjEHjzTWjW7N8Ev/4KQ4bAjh2JnqP3kt74ePjw8aQ/mNemJH2jFmAwdK3alZkdZv5nlKFY\nG8vMRp48VLo9PguXYIGWc1sypO4QupVrBzNmwOrVjuu2bZtoweiqfBPyceCJAxTOVZjoTRs5/WB7\nFi54jRGNRyT5c1DhKONTPJI0ZS28/bZjWbIEGjQAYNfpXUzbOo2Y2BgervkwR4KOMHLVSDpX6czO\n0zvJ65WXb3t+m6QCkuJR4qy1fLd3CU3aDGLHKwO4e9AbCRYcMrIVf6/gwa8eJDw6nMK5ClP7wGWm\nLLjEsaI+vHhnGCcrl+SKjSYwJBCDoVCuQrxy1ysMazDsuntV4SiDU7M6kWxk2DCoUAFGxClUTJ4M\nhw7BBx8ketjF8It0W9SNFp86hgY/9fwwJracSP9l/SmYsyAftos3etyFC0SULk6zUaV5sOkQfj/5\nO/vP7Wdj/43kzJEz2dke9O0gTl4+Sa0itXh741v8M+EKnV+6jc8e/4ky+cvc8Fg1q8s8VDiSdPHN\nNzBgAEydCj17XrfLWovfW36s7bOWmkVqEhMbw52z7uT5O5+nU5VOLspwFjJ1KixfDitWuDont+xS\n5CUeWfIIy/5chjGGO3yr0nbNcfqsCyYqd06W3OHNBxUuElG4ALHEcinyEj4ePgypN4Tn7nyOQrkK\nqXCU0SkYiWQjs2bBmjWOARCu6twZunaFXr1uevgHr95H25WHKLf5LwB2n95N96+688ewP65Ld2bK\nG0R89w2zXm7DhbALFMtTjMcbPE4erzwJnfamIqMj6fdNPxbtW0SdYnX4fokPv9YqyBu3nWTTo5uS\ndA7VHGV8ikeSbnbtgg4d4KGH4PXXwc3xVj86Nhrv/3kT8UrEtZqi/t/0p3HJxgysO9CVOc78LlyA\nypUdMah6dVfnJsUOXzzMoG8HEXA0AIMhOvoKHU7mZsB+H5psPcOO8rn44J7c/FDkElExUcTaWAyG\nCgUr8NcTf2kobxFJOxrZKhnq1oXff/933VrYuBGaNEnS4Vfq1KbQH0exsY5ReXac2kGRXEWuS/PN\n/m/Y//5oFtV0Z8HeBQRHBfNC0xduuWAEjqFVG5dszIA7BvDbgN/w69iLNoc92HxiM/pjWkSSrWZN\n2LwZ1q93vCAKCQHAw82Du8rcxcs/vUxEdASbjm/iu7++o0mpJvx1/i/2nNlDdGx0oqdVPLpeUEQQ\ng74dRIPpDfj+kcaEtG+dJQpGAOV8y7H6kdXsHLKTzlU64+buzjclQ+h23yVqvOTLV+XD+eDzIGYt\niaVgpDuVClbCzbjxz6V/bnheFY4yCA2dKpJNVKsGkZGO4W0BDh4EDw8oc+OmaVf1a/Milzwtj77d\nnD5L+zBy1UjeavXWtf3WWl7+vD9Nznoz8n8/sX3wdrYFbmPVwVUpznqpfKXYfHKzY46lVq2IWb2S\nkrmL37Qzr+Y5ylwUjyTdFC4MP/4IBQs6hv0+cgSAL7p8wdbAreQdn5fui7rzUfuPeGnNS9wz5x46\nLehEoxmNkjWdQHYVa2N54MsHsNbySaVn8P/5OPdW2ERoVKirs5aqqhaqypddvyTspTDG+o8ld47c\nnIoJZkY9d2oMM7j75CZgehTe/5yiTmQdPNbduO+amtVlAGrGIJLNfPopzJkDa9fCZ5/B99/DwoVJ\nPjy6Syc21y/Or81uY1vgNi5GXqRhiYa80PQFrLW82MGHyYUexsyZA0CfpX1oXqY5/e/on6Jsx9pY\nei3uxa7Tu3A37nz96h7GPlGDEYNnU6dYnZser2Z1GZ/ikbiEtTBlCkyYAF98AXff7dxsMcbw5sY3\nWXd0HUt6LCGHWw5GrBxBUEQQszvOdm2+M7hjwcdoOKMhJ54+gVvnLtCwIY38ljKh5QT8y/q7Ontp\nJiY2hmV/LuOFH1/grwt/kcMtB0/8GssTm2KYNLEjJ/MZljy4RM3qREQyjN69HcN2r1njaFJ3553J\nOtyjURPq/WOZtm0aRfMUZXDdwfz6z6/0+6YfXh5e9PkzJ0tqe2GtZd/Zfaz8eyX1itdLcbbdjBtf\ndPmCO4reQWRMJMbfn8FhVWk9rzV/X/g7xecXkWzKGBg+3NEXs2dPR0HJWTAC2Ht2L12qdMHT3RNj\nDD2q9bhuMllJWA63HERGR3Jl6dewezexTz5BSFQInu6ers5amnJ3c6dTlU78+cSf/NL/F+qXqM87\nDWL5uC4MGbWMS+dP3vB4FY5EJMXUBCeZPDxgzBh45ZVbKhzRoAGhG9dQNHdRJt07iQ6VOrC4+2KW\n7l/KpVNHqXEG/uexiTzj89BoRiMm3TuJmkVqpkrWjTGsPbKWFQ+voEK73jT9x42e1XuyeN/iVDm/\niGRjLVvCL784atf79oXwcACq+FVh6f6lRMdGY61l8R+LqexXmaiYKC6GX7yu36Pi0b+K5SlGl5Kt\nCBrwMMuf6UjP5f0p6FOQBiUauDpr6aZRqUZs7L+RoOeDyPXKGPaV8qbPtM03PEaFowxCbbxFspke\nPRwdkI8cgVq1knds3brk+fMIHtGx/9nlsWETbo2bsPXx3Zx4+gQXn7/IwzUfTp08X72GmwfhV8Id\ng0j88gvhV8JvOP+I+hyJSJKVK+d4aRQVBXfdBceOMaLRCKJjo7n9vdup/lF1fjz0I1X9quI70Zcy\n75ah/vT6HAs6xr6z+9h3dl+W61OTEp9sKsjZO+9gYdFzVC9UnR8e+iHJE+pmJXm98/Kq/2t0X3uG\nh65UvmFa9TnKANTGWySbWr4c5s+/fljvJIqpXo3ubUMp17Irzco04+PfP8Y3py+fry8ERYrAiy+m\nQYYd3vnlHaZtm8aDVbrz1P3/o8YwN7rcPYy37n2LHO45Ej1OfY4yPsUjyTCuThj71lswZw6x97Zi\n39l9RMVEcSniEr2X9mZ93/WUzV+WsevG8uHvH+Lp7knhXIW5GH6RlQ+v5PaCt7v6Llxr40bo1g32\n7gVfX1fnJuM4dw5TqJD6HImIZDht295SwQjAvWEjPi0+lNCoUD7Z+gkVClSg3e3tOLp0Dv0vz6Pr\nwq4cvHAwlTPsMLzRcAbWGcj4TRPZXyE/C0oMZ8/ZPYxcNTJNrici2ZAxMHIkLFgA/fvj9vr/qO5X\nlTrF6rDt1DY6V+5MOd9yGGMomrsoZ0PP8tfjf7F10FaebPgkQ74f4uo7cK2ICBg40DHpqwpG1/Pz\nu+FuFY5EJMXUJNQFGjYk79Y9fNT+IzpV7sSXe75k3OLh5D8VRJ7G/tQvXp+759zNhfALqX5pYwxu\nxo1+tfvRqNsImvwDczvNZd6uWyvoiYgkqnlzx9xwP/4I7dvD+fOUzFuS3078xpWYKwCsPbIW35y+\n5MyRk4CAADpV7sT+c/tdnHEXGz7cMZdUly6uzkmmo8JRBqE+RyKSLK1bw8qVnL98hmdWPcOq3quo\n9WcwOZr58/kfC+hRvQfVC1cn4EhAmlze092TS1GXoHFj+PVXgiODEx0BSX2ORCRFihWDn36CqlWh\nTh26XCxGibwlqDOtDl0WdmH5geX4+fgREuWYSHb+7vlUL5w1Jjq9JfPnE7l6Bd8/8wAHLx5ydW4y\nHfU5ygDUxltEbknt2hx8/WnaHxvPtkHbmHNXHvo8MJoWhZczocUEXl77Ms/f+TztK7ZP9UufCztH\nvWn1eKhEG0b1m02tN0ozsP5gnm78dKLHqM9Rxqd4JBnesmUwcCCxw4eztls9zjvneRu/YTyL/1hM\nwZwFAVjx8ArK5i/r2ry6wv79hDaqS8cBufCq04DfTvzGO63fSfWBeTK7G8Wj7DdchYhIVtGhA6XW\nbed8kfNsOLaBjify0id4Nrs5zUe/f0RIVAgtyrVIk0v7+fjxy6O/8OamNznrFUPJs5FsPL6R9hXb\nU7FgxTS5pqSP0aNH4+/vj7+/v6uzIvJfHTpA7dq49epFi4AAx0Ta+YvwcfuPebHpi1yOusztBW7H\ny8PL1TlNf2FhRHTuwKhWHnw+Zg+FcxVm39l9NJ7ZmE6VO5HLM5erc+hyAQEBN22ppWZ1IpJiahLq\nIh064Pn9ChZ1W8SMCd0IjbjEityBNCjRgPK+5QnoE0DOHDnT7PLF8hTjTOgZjpcvyNSi/WhSsgl3\nz7mbM6Fn0uyakvauFo5EMqzSpSEgAOrVgzvugB9+AKBM/jKc23cuexaMYmLgoYc4X7k0O++vT+Fc\nhQGoWqgq+bzy6XvZyd/f/6bNvFU4EhHJrOrUgUuXaB5emC9+K0W+dz/m7AsX+OmRnxjXYhz5vPOl\n6eWjYqL4at9X1G03kCrHwxnZZCSNSjZi9cHVaXpdERE8PGDcOMd0CEOGwLBhEBaWYNK/L/xNs0+b\nUfDNgjSZ2YR9Z/elc2bTmLXw2GMQEkLsxx+z/dQOtgVuA2DZn8uIjo2mRN4SLs5k5qHCUQahARkk\nM9NbZhdxc4P774fevXHLlx+/Xo+m6xtTN+OGwRBeozJs3w5ASFTIfwZm0IAMIpJm/P1h504ICoI6\ndfDPnfu63ZHRkdw37z46V+nMH8P+oE+tPrT5vM21wRuyhDFjHCP6ff01pQpXYPr902nxWQuKTS7G\n0O+H8nWPrxMdMEf+SwMyZADqACsit2z5cmjXDn79FRo2TPfLP7vqWfbtWM1XEw4x8rOHWHs0gC0D\nt5DbM/d/0mpAhoxP8UgytS++cAxh/cgjMHo05MrF7tO76f5Vd/4Y9se1ZPWm1eP9tu/TqGQj1+U1\ntXz0EUyezLwPh/L8nrcJjQqlS5UuvN36bS5HXaZo7qJ4uGmIgfhuFI9UcyQiKaZaTxdq1cpRNIIb\nDAAAHsFJREFUQHJBwQjgzVZv8sA9Q7kSE0WJy4YN/TYkWDBi1ar0z5yIZC89exLw8ccQGAg1asCq\nVeT3zs/Z0LMERwQDEBoVSmBIIPm987s4s6ngrbdgwgTWf/IyL++byvJey/nrib84H36eV9a8Qsm8\nJVUwugX6xEREMrMcOaBNG5dd3hjDoHqDof4iXvbtAD4F/5Mm/Eo4Ma8874LciUi24+sL8+bBihUw\neDClmjZlULtONJvdjHa3t2PVwVW0rdCWSgUrAbDh2AZeW/saQRFBtLu9HaP8R2X8AoW18Nxz8P33\nsHEjy/a9y9B6Q6lVtBYAb7R4g45fdnRxJjMv1RyJSIqpz5FEV67IH+sWs+SPJVwIv3Bt+5/n/qTa\nh9W4eHivC3MnItnFtXh0332wZw8ULsy4p5Yx82Izcrp789ydzzHt/mkYY9h7Zi+dFnRiYJ2BfNTu\nIzYc38DzqzP4i5zoaOjXDzZsIHzNKnZ6nMfduPPn+T+vJdl/bj8FchZwYSYzN/U5ygDUxltEMrNL\nkZd479Hq1DgayUdD6rD3zF7W9FlDhQIVaDqrKQ9W68Hjdz+PCQ9Xn6MMTvFIsqStW2HgQChQACZP\nhlrOGpaf3+B82Hkmt54MwOGLh2n6aVNOPH3ClblN3D//QM+ekDcvf3wwhjZLu+KTw4fAkEDcjBt3\nl72bYrmL8eXeL1nUbRH+Zf1dneMMS32ORCRNqc9R9jZ502SoWpX7o8vzw0M/8GTDJ3l29bMA7Du7\njwdLtQV3dxfnMvsyxrxojFlkjDlkjIk1xhx2dZ5E0kqC8ahuXdi82TGBbOvW8NBDcPAg3h7eXIy4\neC3ZhfALeLln0DmSvv/eMa9Tmzbw7bc8vGowL9/1MvuG7ePgkwcpkLMAxfMUp2z+sqzru04FoxRQ\n4UhERFLk2KVjlGrYGrNvH1hL09JN+efSPwBU9qvMj5vmQrFiLs5ltjYO8AcOABcBVQ1J9uPhAU8+\nCQcOQKVK0LAhQz7dw64dKxmxYgTvb36frou68kLTFwDYeWon98y5h8rvV6b/N/25FHnJNfkOD4dn\nn4WhQ2HRInjpJXBzY++ZvfSq0QuAAjkL0LZCW8rmL8vIJiOpWqiqa/KaRahwlEFoniPJzNTnKHtr\nVKIRHx5eQGxOb64cP8p7m9+jYQnH6HkzO8xk2rwpDDp/xLWZzN7KW2sLWWtbA4GuzoxIWrppPMqT\nB157DfbvxyePL5unhHH/nF85eGgrU++byqC6gwi8HEjrea15qMZDLO6+GIul1+Je6ZL/a6yFJUug\nalU4ehS2bYO77rq2u5JfJZbsXwI4mjb/dPina4NMSMqoz1EGoDbeIpKZxdpYhq8YTudhH/B28xzY\ne1vxRZcvrg3pHTH3U0K/+gK/ZavV58jFjDF7AB9rbflE9iseSfZy/LhjTqSvv4YePWDoUL5w28dX\nf3zF4u6LAbgSc4U84/Nw8fmL5MyRM+3ztHcvPPUUnD4NU6bAPfcAMH/3fF5d+yohUSHcVfouNh3f\nRIm8JTgefJwe1Xrw7n3vYoy+YpNCfY5EJE2p1jN7czNuTG0zlTvvG8jCaqP5tue318115H32IgXL\nqZmHiKS9ZMejUqVg5kxHgaRECWjfntaPjOGOVbux4eEA10bgzOGeI5VzG4e1sH49dO4M/v7QsSNs\n336tYLTuyDqeXf0sn3f+nK2DthIVE0X7iu2Zet9U1vdbz5Q2U1QwSiUqHImISKrIUb0m3gcS6Osf\nGKg+RyKSsRUvDq++CocPk+fV/3HPb2e4XCQ/27rfxasvNuSVBs/g4ebB3J1zaTijIfWm1eOjLR+R\n4prWS5dg9myoUwcGDYKWLeHwYXj8cUc/KacVf69gcN3BNCrZiJJ5SzLp3kn8eOhHGpdqTMWCFVOW\nB7lOBp/lSkQyA/U5EgDKl3e0kY/v6mz1IiJpLMXxyMODHJ27Uqt9G+YtHUupb9fxxgYvCs6ewuna\n33Ok4FGmDJ1EVKXbGPTdYLw8vOh/R/+kn99a2L8fli93jEC3ZQs0bw5vvOEYSc8t4XqL/N752Xdu\n37X1vy/8TX7v/Cm7V0mQCkciIpI6ypVzvPGMTzVHmUrfvn0pW7YsAPnz56d27drX/uC82mRJ61rP\n6uu5PHNRqXAbeLQNfv7+EBTEmMFNaRxYmkaDxkBwMJN8c/H3h89Ax7+hQgUCzp8HY/CvUQOsJWDn\nToiIwN/TE/76i4AtW+D4cfzz5YN27Qho2RKefRb/Nm0SzM9Pa35i1cFVeN7mSYWCFXhv4Xvcs/Me\n6jauy5ydc3im+DMEBARkiM8ro68HBAQwe/ZsgGvfb4nRgAwZgDrASmYX98tZsq/th3+hWsU78X3N\nk2pFazL9/ul8f+B7+vR6kxXj+jGgz7sakMHFNCCDZHVpGY/6Lu1LzSI1ebrRCDh/nm+Xv8uBLSt4\nuvADcPAgnHBOHmvMv4uPD1SoABUr/rsUKeLYdwOxNpZui7pxLuwcd5e9m8V/LKZFuRaU9y1PSFQI\nbSq04Y5id6TJfWYHNxqQQTVHacAY4wY8BQwGygBngYXAa9baMFfmTUQkLQRHBNPu6y78XSA/Zx7c\nwOygtTSe2ZgW5VowMvgKP4bvdXUWJYlGjx6Nv7+/XniIxDOi0Qhazm3JpchLeLp78s7pT1jy7BK2\n5fDh6ZVPcyrkFE1LN+Wd1u+QxytPiq615cQWdp/ezZ7H9uDp7smTDZ+k7LtlOTL8CAVyFkilO8p+\nAgICrtUoJUY1R2nAGDMFeAL4GvgBqOpc/xloGf+1nN7UiUhm9/PRn3nux+f4ZU4OGDuWXVUKUG9a\nPbb2/YUa5RsRFR6CVw5v1Ry5mGqORFJm75m9zNo+ixgbQ++avSmauyh1ptVhYsuJ1C9en4kbJxIc\nGcw3D36T7HOHXQlj/M/j2X9+P97u3hwNPsr6fusBsNZS8p2SbOq/iTL5y6T2bWU7qjlKR8aYajgK\nQouttd3ibD8MTAUeBL5wUfZERNJEQZ+CHAs+xpUyzclx6BBnS8cQY2MoFOMNefOSw93T1VnMtowx\nvXG0YgAoBOQwxrziXD9irZ3nmpyJZD7VCldjcuvJ19bn7pyLf1l/+tbuC8CMDjPIMz4PoVGh/PD3\nD5wJPUPT0k2pWaTmDc8ba2O5/4v78fPxo0uVLszfPZ/NJzYzfet02tzehulbp1PIpxAl85ZMy9sT\nNJR3Wujp/PluvO3TgTDg4fTNjkjau1kVtWR9VQtV5YFKDzDjwk+s/PEThnw/BD8fP2aue4cIH0+e\nWvGUq7OYnfUHxjoXPyBfnPX/DLM1evRo/U5LppXez65PDh9Oh5y+NqT3ubBzuBk3Oi/ozNu/vM3O\nUztpNbcVi/ctvuF59p7Zy5GgI8zvPJ8Hqz/I1z2+pkDOAnz4+4fUn16fLSe38F2v73B3c0+P28qy\nAgICGD169A3TqOYo9dUHYoDNcTdaayONMTud+0VEspwP2n7A9j2eeK/5mQ/afkDtorWZOn0Ax20w\nIVEhrs5etmWtvTs56W/2h4OI/KtdxXZM2DiBnot7Ur94fWZun0nHyh05GnSUn/v9jLubOwPqDOD+\nL+6nS9Uu/zk+JjaG0CuhRMdGk8MtB27GUW/hZtzImSMnczvNpXrh6ul9W1nW1f6UY8aMSTSN+hyl\nMmPMbsDPWvufcWuNMQuBroCntTY6zna18RaRrOHnn+H552HTJsf6hg2O9Y0bb9jGWzIGxSOR5AuJ\nCuGDzR8QGBJI09JNOR1yml2nd/HJ/Z8AEBUTRa43cvFVt6/4bNdneLh58Hj9x/nn0j8M+X4I0bHR\nlMtfDjfjRtPSTelSpQsL9i5gz5k9rO+3Hg831WWkNvU5Sl8+QGQi+yLipLmUPtkREUlH5crBoUP/\nrl++DHlSNmqTpC+NVieSPLk9c/N80+evre86vYux68cysO5AahapyStrXqGKXxWGLR/GhJYTCL8S\nTscvOxJLLBv6baB64epM+W0KM7fPxFrL6+tfp1qhaix/aLkKRqlMo9W5QBJqjroAXvFrjvr06aNJ\n97Seadd37NjB8OHDM0x+tO7C9Z9+glat8I+OJmD9emaPGgVHj1K2b1/GjBmjmqMMTjVHktkFZJB5\n9xbtXcQTPzzB+fDzNCvTjOiYaIY3Gk6nKp0AeGjxQ/x64lcOPnkQcIxG5z3Om4vPX8Qnh48rs54t\n3KjmSIWjVGaMWQncg2Oo1Cvx9m0EKlhri8TbrmAkmVpGCUaSQeTLB0eOgK8vzJgBv/wCM2eqWV0m\noHgkmV1Gi0exNhY340arua0YVn8YHSt3BODx5Y8zb9c8AkcGkjNHTrae3Eqrua04/9x5zE0miJWU\nU7O69LUZaAU0BDZc3WiM8QZqAwEJHaRmDJKZ6bmV6xQoABcvOgpHly8TEBxMgDr5i0g6yGjx6OoA\nC4PrDuaJH54g/Eo44dHhLNy7kMYlG1NnWh1qFqnJmsNrmNFhhgpGGYBqjlKZMaY6sBNYYq3tGmf7\nE8AU4GFr7fx4x+hNnYhkHXXqwPTpULcujBkDMTEwdqxqjjIBY4wdNWqUXtaJpIElfyxhzs45jgEZ\nGjxO8zLNWXtkLYGXA6lXvB6V/Cq5OotZXoCzz9GNmnmrcJQGjDFTgceBJcAPQBUcE8NusNbek0B6\nBSPJ1DJaMwZxsRYt4IUXoFUreOYZAkJCCChaVH2OMgG9rJPMTvFIkuJGL+s0CWzaGA48A1QD3ge6\nA1OB9okdcLVZnYhIpne1WR3A5cv433GH5s4REZFMQX2O0oC1NhZ427mIZHkq2Mt1fH3/LRxduqSh\nvEUk3SgeSUqp5iiDGD169LUhcUVEMjVfX7hwwfHvy5cJOHpUNUeZiOKRiGRVAQEBN41H6nOUAaiN\nt2R2auMt15kwwVE4evNNaNYMXn8dmjfXgAyZgOKRZHaKR5IU6nMkIiLpJ16fIzWrExGRzEKFowxC\nzRgkM9NbOrmOry+nju/n3rn38s/J/byw7B1GjRrl6lyJSDageCQppQEZMgi1xxeRrOBi+EVe++U1\nOv25n4AjOcgb6c6WXNu5r0lNGOvq3ImIiNyYao5EJMVU6ylXDVs+jLPeMVTxKMqR4UfwDIukQ/2H\nmbNzjquzJkmklgySmenZlRtJyoAMKhxlEApGIpIVbDq+idvLN8DrchjFvQvhGWtYtm4V55afc3XW\nJIk0756IZFX+/v4qHGUWCkaSmenZlauK5ylO0dJVyBF0mYk/vEKIp2G79w7eGveWq7MmItmA4pGk\nlApHIiKSaqa2mcrrO6bgcwW+3jSTYC/LtPun0btWb1dnTURE5KZUOBKRFFOTULmqXvF6bB2yneg8\nuXin6tOUKFGZrlW7ujpbIpJNKB5JSmm0OhERSVUl8pYAvyI0sSUgbz5XZ0dERCTJVHOUQWhABsnM\n1MZb/sPXF44ehTx5kjQ6kIhIalA8kpQy1lpX5yHbM8ZY/T+ISJbSujWEhUGRIvDVVwAYY7DWGhfn\nTG7AGGNHjRqFv7+//sgUkSwnICCAgIAAxowZk2g8UuEoA1DhSDK7gIAA/SEl15sxA7ZuhR49wPls\nqHCU8SkeSWaneCRJcaN4pD5HIiKS+gYMcCwiIiKZiGqOMgC9qROR7EA1Rxmf4pGIZAc3ikcakCGD\n0IAMIpJVaUAGERHJLFRzlAHoTZ1kdmrjLUmhmqOMT/FIMjvFI0kK1RyJiIiIiIjchGqOMgC9qROR\n7EA1Rxmf4pGIZAeqORIREREREbkJFY5EJMU0mIiIiGQEikeSUiociYiIyDUaPVVEsqqkjJ6qPkcZ\ngNp4i0h2oD5HGZ/ikYhkB+pzlAnoTZ2IZFWa50hERDIL1RxlAHpTJ5md5pWQpFDNUcaneCSZneKR\nJIVqjkRERERERG5CNUcZgN7UiUh2oJqjjE/xSESyA9UciYiIiIiI3IQKRyKSYhpMREREMgLFI0kp\nFY5ERERERERQn6M0YYxxA54CBgNlgLPAQuA1a21YAunVxltEsjz1OUp/ikciIv+lPkfp7x1gMrAH\neBxYBDwJfGuM0R8GIiKSXhSPRESSQYWjVGaMqQY8ASy21na11s601o4EngbuBh50aQZF0oDaeItk\nPIpHkh0pHklKqXCU+no6f74bb/t0IAx4OH2zI5L2duzY4eosiMh/KR5JtqN4JCmlwlHqqw/EAJvj\nbrTWRgI7nftFspSgoCBXZ0FE/kvxSLIdxSNJKRWOUl9x4Jy19koC+04AfsYYj3TOk0u4smo7ra6d\nGue91XMk97ikpk9KuuzYTMHV95wVn2FJd4pHTopHqXsOxaP05ep7zorP8I2ocJT6fIDIRPZFxEmT\n5SkYpe45MnIwOnLkSJKulZkoGKXtOSRdKB45KR6l7jkUj9KXq79zs+IzfCMayjuVGWN2A37W2mIJ\n7FsIdAG8rLXRcbbrP0FEsgUN5Z1+FI9ERBKXWDzKFtXp6ewkUNkYkyOBpgwlcDRxiI67UX8siIhI\nGlA8EhFJJjWrS32bAXegYdyNxhhvoDbwuysyJSIi2Y7ikYhIMqlwlPoWABYYHm/7QCAn8Hm650hE\nRLIjxSMRkWRSn6M0YIyZimMm8iXAD0AVHBPxbbDW3uPKvImISPaheCQikjwqHKUBY4wbjjd1g4Cy\nwFkcb/Bes9aGuTBrIiKSjSgeiYgkj5rVpYAxxs0YM8IYs98YE26MOWaMmQR4W2vfttZWttZ6W2tL\nWWufiR+IjDHdjTGfGmN2GmOuGGNijTGlb3C9fMaY94wxJ5zX22OMGZLmNypZljHmRWPMImPMIefz\nd/gWz9PWGLPJGBNijDlvjFlojCmbSFo9xyJpIKGYBLwJfHyzeOQ8XjFJXEbxSDIKFY5S5h1gMrAH\nR7OFRcCTwLfGmKSM+DMU6A6EAn/jaBueIGOMJ7AaGAx84bzen8CHxphRKbgHyd7GAf7AAeAiN3gG\nE2OM6Qx8B3gBzwBvAc2AjcaYYvHS6jkWSTuKSZKZKR5JxmCt1XILC1ANiAUWxdv+uHN7zyScoxTg\n5vz3+87jSieS9jHn/mHxtn+FY5K/BI/TouVGC1A2zr/3AIeSeXwO4ARwGPCJs70WEA18Ei+9nmMt\nWtJgUUzSktkXxSMtGWVRzdGt6+n8+W687dOBMODhm53AWnvcWhubxOv1wvE2b3q87e/i+ELokcTz\niFxjrT2SwlM0B4oBM2ycZjrW2p1AANDDGBN3PjU9xyJpQzFJMjXFI8koVDi6dfWBGBzzSFxjrY0E\ndjr3pwpnh9o6wHZrbVS83VtwVD3XS63riSTD1ef8lwT2/QbkBSqCnmORNKaYJNmd4pGkChWObl1x\nHLOLx591HBzVun7x3lCkhC/g7TzvdZyB7zyO2c5F0ltx58//PJtxtl1No+dYJO0oJkl2p3gkqUKF\no1vng6NNakIi4qRJrWtxk+ul1rVEkuNGz2b83wM9xyJpRzFJsjvFI0kVKhzdujAco6EkxBtHlWxq\nzSFx9Tw3up7mqxBXuNGz6R0vjZ5jkbSjmCTZneKRpAoVjm7dSRzNFHIksK8EjuYN0al0rYtAOAlU\n8RpjvAA/Eq5GFklrJ50/E2p+cHXb1WdTz7FI2lFMkuxO8UhShQpHt24z4A40jLvRGOMN1AZ+T60L\nOUcP2gbUcY7LH1cD589Uu55IMlzt/N0kgX2NgGDgL9BzLJLGFJMku1M8klShwtGtW4CjmcLweNsH\nAjmBz69uMMYUNcZUNsbkTMH1vsDR/nVQvO3DgSvO/IikmUSe43VAIDDAGJMrTtpaOCbzW2StjYmT\nXs+xSNpQTJJsQ/FI0pKxNtkTEIuTMWYqjgn2lgA/AFWAJ4AN1tp74qSbDTwC3G2tXRdnezMcMzcD\ntMfxtmIyjrcb1lo7Lk7aHMAmHJOZTQX2A22BjsDr1lrN5izJZozpDZRxrj6BY26Ht53rR6y18+Kk\nnU3Cz3FXHEFkJzADx3CpI3AMK1zXWhsYJ62eY5E0opgkmZnikWQYrp6FNjMvOGrensbxCxUBHAcm\nEWdmZme6T3H8YjaLt30UjtmZY537Y+KuJ3C9fMB7ONrBRuCYQfoxV38OWjLvAqy9wTO4Jl7aBJ9j\n5752OOaWCAUuAAuBcolcU8+xFi1psCgmacnMi+KRloyyqOZIREREREQE9TkSEREREREBVDgSERER\nEREBVDgSEREREREBVDgSEREREREBVDgSEREREREBVDgSEREREREBVDgSEREREREBVDgSEREREREB\nVDgSkVRgjJlojDlkjPFIp+vVNsbEGmOapcf1REQkc1A8kpRS4UgkiYwx/s4vwMSWK67OoysYY8oB\nTwJjrLXRyTzWwxgTaIw5faNAZowp7/yMVwJYa3cAS4DJKcm7iEhmpHiUMMUjSQ3pUqoWyWLmA8sT\n2B6b3hnJIF4AgoF5yT3QWhttjJkNPA+0B5YmkrSv8+fMONveBdYZY9paaxP6/xARyeoUj66neCQp\nZqy1rs6DSKZgjPEH1gDPWGvfdnF2rmOMyQlEWWtj0vm6eYGTwHRr7YhbPMftwJ/At9baBxLY7wYc\nBnyA4tbaK3H2HQL2WGs73Mq1RUQyI8WjBK+reCSpQs3qRNKAMaass9p9lDGmvTFmizEm3Bhz0hjz\npjHGPYFjbjfGzHVW60caYw470/rESzfbeW4/Y8wsY8xpIAQo4dxf0xizyhgTYow550zv5zzmU2ea\nwsaYKGNMgm/XjDEfGGNijDGlb3KrbXEEiQTflCXlnqy1B4CfgTbGmCIJnKYFUAqYHzcQOa0E7jPG\n5LpJPkVEsiXFo6Tfk+KRgJrVidyKXMYYvwS2R1prL8fb1hZ4DPgImAF0BJ4BLgLjryYyxtTF8Rbw\ngjPtCaA2jrbTdxpjmifQfno1EAiMAXIBoc63Xj87909xnqcd8INzmwWw1p4xxnwDdDbG5LPWBsfJ\nizfQC1htrT12k8+iufPnlvg7knlPs4C7gN7ApHin6uf8OZP/+hUYDDTFEZhERLITxaN/KR5J6rDW\natGiJQkL4I+jHXdiy7I4acs6t10GSsc7z27gZLxtO4F9QK542zs6z9MnzrbZzm2fJZDHhUAM0Dje\n9i+dx8yKs62Vc9vQeGkfcm7vmoTPZB1wLpF9ybknHxztxPfGS5sfCAd+T+QaTZ3nGuHq50OLFi1a\n0mtRPErwM1E80pIqi5rViSTfJ0DLBJaXE0i71P73bVcAUPRqVb4xpgZQA/gCyOlscuDnfBu4EQgD\n7k3g3Ne90XI2jWgLbLbW/hIv7X9G0bHWrsbRdvrReLseBc6ReGfUuArheBN3neTek7U2DEfArGKM\naRDnVA8CXiT8lg7gvPNn4STkVUQkq1E8+pfikaQKNasTSb4D1to1SUx7KIFtV79AC+L4Uq7iXB/j\nXBKS0JftX/HWC+F44/VnEtJeNQMYZ4ypZa3daYwpj6Npwrs2acOgWsAksP1W7mkmMBDoD2x2buuP\n403d/ETOcfXaGllGRLIjxaN/KR5JqlDhSCRt3Wi0HhPv5yRgRSJpL8bfYK2NSEG+rpqFI1g8iqPt\ndX9nfmYk8fizQM0Etif7nqy1m40xe4EexpjhwG1APRwdX4MTOgFQIE4+REQkcYpHikeSBCocibje\n1bdoscl4A5iQs0AoUCmBfQltw1p72hjzLfCQMeYFHPM3/Gqt/SOJ19wDNDPGFLDWxm3OcKv3NAtH\nk4vOQJ042xJTIU4+REQkZRSP/qV4lE2pz5GIi1lrt+P4Mh1iHLN7X8c4Zu32jX9YAueJwTEKUENj\nTJN4u0feIAvTAV8cbdeLk/S3dABrnT8bx8vLrdwTwFzgCjAIR0fcwzcJZo2c6TcmI88iIpIAxaPr\nKB5lU6o5Ekm+usaYhxPZt8RaG3oL5+yNY5jRXcaYWThG1fHB8SaqE45Zvz+Lkz6hdtUArwCtgRXG\nmPf5d+jUQs79CbWFXgkcxfHlfxlHR9SkWuE8pi3wfQrvCWvtOWPMMqCLc9OoxC5sjDHAfcAKZwda\nEZHsRvHoX4pHkipUOBJJuqtf5A8CPRPZv4mEO73GT3ddUHB2Pr0DeBHoAAzB8SV/GPgU+OlGx8c5\nz1/GmGY42lY/BUTgmBBvGHAQR2fS+MdYY8xMYCywMDlf7NbaUOfEfT2MMcNtnEnxknlPcc3EEYxi\ncAwTm5hmQGlgaFLzKyKSRSge/fdYxSNJFcZaDaohktU5J8DbArxgrX0zgf3PARNwzEfxWzLPXQbY\nDzxurU1siNNUZ4xZApSw1ja4aWIREckQFI8ko1PhSCSLMcbktNaGx1k3OJomdAPqOttfx03vgWO4\n1cvW2tq3eM3xQA+gYhKHXE0R5xvA3wF/a+3PN0svIiLpT/FIMiMVjkSyGGPMnziaCOwBcgH345i5\n+0trba846coCTYAHcASqB621C9M7vyIikjUpHklmpD5HIlnPUhwBqDeO3/FDODrGToyXzh/HsKRn\ngTEKRCIiksoUjyTTUc2RiIiIiIgImudIREREREQEUOFIREREREQEUOFIREREREQEUOFIREREREQE\nUOFIREREREQEgP8DCZEImILVPOIAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Note that lmfit should be Version > 0.9.0. We will check this\n", "import lmfit\n", "from distutils.version import StrictVersion\n", "assert StrictVersion(lmfit.__version__) > StrictVersion(\"0.9.0\")\n", "\n", "from lmfit import minimize, Parameters, printfuncs\n", "\n", "# Choose the model we want to fit it with \n", "model = 'Drude-Lorentz'\n", "\n", "params = Parameters()\n", "params.add('Omega_p', value = 9, min = 8.5 , max = 10) # Omega_p has a value in eV, we will add a starting guess here\n", "\n", "# Drude Term\n", "params.add('f0', value=0.5, min = 1E-10, max = 1) # f has no units, we will add a starting guess here \n", "params.add('Gamma0', value= 0.03 , min = 1E-10, max = 0.1 ) # Gamma is damping term and has units in eV, we will add a starting guess here\n", "\n", "# Lorentz first Oscillator term\n", "params.add('f1', value=0.5, min = 1E-10, max = 1) # f has no units, we will add a starting guess here \n", "params.add('Gamma1', value= 0.5, min = 1E-10 , max = 2) # Gamma is damping term and has units in eV, we will add a starting guess here\n", "params.add('Omega1', value = 5 , min = 5 , max = 10) # Omega_o is oscillator central frequency term and has units in eV. we will add a starting guess here\n", "\n", "# Lorentz Second Oscillator term\n", "params.add('f2', value=0.5, min = 1E-10, max = 1) # f has no units, we will add a starting guess here \n", "params.add('Gamma2', value= 0.5, min = 0.01, max = 2 ) # Gamma is damping term and has units in eV, we will add a starting guess here\n", "params.add('Omega2', value = 3 , min = 2.5 , max = 3.5) # Omega_o is oscillator central frequency term and has units in eV. we will add a starting guess here\n", "\n", "# Lorentz Third Oscillator term\n", "params.add('f3', value=0.5, min = 1E-10, max = 1) # f has no units, we will add a starting guess here \n", "params.add('Gamma3', value= 1, min = 0.01 , max = 2) # Gamma is damping term and has units in eV, we will add a starting guess here\n", "params.add('Omega3', value = 4 , min = 3.5, max = 4.5 ) # Omega_o is oscillator central frequency term and has units in eV. we will add a starting guess here\n", "\n", "# Lorentz Fourth Oscillator term\n", "params.add('f4', value=0.5, min = 1E-10, max = 1) # f has no units, we will add a starting guess here \n", "params.add('Gamma4', value= 1, min = 0.01 , max = 2) # Gamma is damping term and has units in eV, we will add a starting guess here\n", "params.add('Omega4', value = 5 , min = 4.5, max = 5.5 ) # Omega_o is oscillator central frequency term and has units in eV. we will add a starting guess here\n", "\n", "# Lorentz Fifth Oscillator term\n", "params.add('f5', value=0.5, min = 1E-10, max = 1) # f has no units, we will add a starting guess here \n", "params.add('Gamma5', value= 1, min = 0.01 , max = 2) # Gamma is damping term and has units in eV, we will add a starting guess here\n", "params.add('Omega5', value = 6 , min = 5.5, max = 6.5 ) # Omega_o is oscillator central frequency term and has units in eV. we will add a starting guess here\n", "\n", "\n", "\n", "# Lets fit to all the data.\n", "w_for_fit = w_exp\n", "eps_for_fit = eps_exp\n", "\n", "# if we were to fit only a part of the data such as fit only below 2 ev\n", "#w_for_fit = w_exp[w_exp<2]\n", "#eps_for_fit = eps_exp[w_exp<2]\n", "\n", "\n", "# Call the minimize function with required parameters and use differential evolution as a method \n", "minimizer_results = minimize(complex_residuals, params, args=(model, w_for_fit, eps_for_fit), method = 'differential_evolution', strategy='best1bin',\n", " popsize=50, tol=0.01, mutation=(0, 1), recombination=0.9, seed=None, callback=None, disp=True, polish=True, init='latinhypercube')\n", "\n", "# If we were to fit ith with least squares methods\n", "# minimizer_results = minimize(complex_residuals, params, args=(model, w_for_fit, eps_for_fit), method = 'leastsq')\n", " \n", "#lets see whether the fit exited successfully?\n", "print \"Print exited successfully? : \", minimizer_results.success\n", "\n", "#lets see the termination status\n", "print \"Termination Status: \", minimizer_results.message\n", "\n", "# lets print the fit report. We dont need lengthy Correlation table\n", "printfuncs.report_fit(minimizer_results, show_correl=False)\n", "\n", "# Caluclate the epsilon based on the fit results\n", "eps_fit_result = np.array([drude_lorentz_model(minimizer_results.params, i) for i in w_for_fit])\n", "\n", "# Lets plot the fit data\n", "plot_fit(w_exp,eps_exp, w_for_fit,eps_fit_result)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }