{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The Fourier Grid Hamiltonian method\n", "\n", "## Source\n", "\n", "This is a quick and dirty implementation of the \"FGH\" method as described in either\n", "\n", "- Tannor's book \"Introdoction to Quantum Mechanics: a time-dependent perspective\" (my warmest recommendation!)\n", "\n", "- BALINT-KURTI et al., INTERNATIONAL REVIEWS IN PHYSICAL CHEMISTRY, VOL. 11, No. 2, 317-344\n", "\n", "\n", "The method computes the matrix elements of the kinetic energy operator in position representation analytically. For position-dependent potentials, the potential energy operator is diagonal in position representation. Both terms are added to obtain the symmetric Hamiltonian matrix, which is diagonalised to obtain both eigenvalues and enigenvectors (wavefunctions in position representation). This notebook aims at an example computation and a sample visualization." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#%matplotlib inline\n", "\n", "from itertools import cycle\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.linalg import eigh\n", "from scipy.integrate import simps\n", "\n", "%load_ext Cython" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define cythonized function that computes the Hamiltonian matrix:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%cython\n", "\n", "# cython: boundscheck=False\n", "# cython: cdivision=True\n", "# cython: wraparound=False\n", "\n", "import numpy as np\n", "cimport numpy as np\n", "\n", "cpdef H_cont(double L, double[:] V):\n", " \"\"\"Compute the Hamiltonian NxN matrix in atomic units (hbar = 1; m = 1)\n", " according to Tannor's book, sec. 11.6.4\n", " \n", " Parameters\n", " ----------\n", " L : number\n", " Length of x-interval.\n", " V : array\n", " Sampled potential.\n", "\n", " Returns\n", " -------\n", " Hij : array\n", " Hamiltonian.\n", " \"\"\"\n", " \n", " \n", " cdef int N = len(V)\n", " cdef np.ndarray[dtype=double, ndim=2] Hij = np.zeros([N, N])\n", " cdef double K, pi = np.pi\n", " cdef int i, j\n", " \n", " K = pi/(L/N) # pi/dx\n", " \n", " for i in range(N):\n", " for j in range(i+1):\n", " if i == j:\n", " Hij[i, j] = 0.5*K**2/3. + V[i]\n", " else:\n", " Hij[i, j] = K**2/pi**2 * (-1.)**(j-i)/(j-i)**2\n", " Hij[j, i] = Hij[i, j] # use Hermitian symmetry\n", " \n", " return Hij" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define potential and sample:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pot = lambda x: x**4 - 20*x**2 # potential function callable\n", "#pot = lambda x: 0.5*x**2 # instructice example: harmonic oscillator potential\n", "\n", "N = 2**10 # number of samples in discretization\n", "L = 10. # length of x support\n", "x_vals = np.linspace(-L/2, +L/2, N, endpoint=False)\n", "\n", "V_sampled = pot(x_vals)\n", "H_sampled = H_cont(L, V_sampled)\n", "\n", "# diagonalize the Hamiltonian matrix:\n", "E, psi = eigh(H_sampled)\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAD/CAYAAABvuWSAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFoFJREFUeJzt3V+MXVd1x/Hfug4YUaidjOgfUJUZpzRGQsiecfwnUoDi\nIPWhElImKW/0wX8A9QmpqIr6WAkMcZUiFanpOFLlJxL/kYx4aojVQiG2y8SDoK1N8diWJUql+E8C\nilInzOrD3Xd85s69c889f+7Z55zvB1353p07M8cz48XeZ+29lrm7AKCtOlVfAABUiSAIoNUIggBa\njSAIoNUIggBa7b4yP7mZzUu6I2nW3Z8p82sBQBalzQTNbKckd/eXJd0xsx1lfS0AyKrM5fBn1Z0F\nStKypMdL/FoAkEmZQXCrpFuJ11Mlfi0AyITECIBWKzMxclvSA+H5Vkk3+99gZpzZAyri7pbn4+3d\nv+16+1dp337d3afzfL2ylBkEX5Q0J+mspG2SXhr0pjfvrujE0g3dXXF1cv1Iyndm4Vl95tCXqr6M\nVOp0rRLXW6ZB13pgz3T+T/z2r/SeHX+R6q1vLX3zwfxfsBylLYfd/aIkmdl+SbfdfWnQ+77z01/o\nqR1/oN9QyAGoH+uke0Ss1H2C7n5s1Htev/uOnr9wTQd2T+v44nVtMot+Rggg6Gyq+gpyiyJEv+e+\njr7z01/o6s3/i3pG+PDs3qovIbU6XavE9Zap1Gs1S/eIWOVBsGPdx+t339HcB9+nz809qLdXXCsR\nxsLtc/uqvoTU6nStEtdbplKvtQHL4WiurhcIn79wTZ+be1Dv7liUgRBAAjPB4vWWxiRLgBpgJlis\n5NK4NyOMdWkMQMwEy1SXZAnQaswEy1GnZAnQap1N6R4RizII9pAsASLHcngy+pMlBEIgEiyHyzco\nWcKMEIhEziAYii8nXx8Jfx5KjM2b2X4z+/K4Y2lEHwST2D4DRKY3Sxn1GCDUFTjRN3zYzP5b0pXw\nnmSF+ttmtjPF2FiV7GsTBNk+A0Qox0wwBKwrfcMH3f3D7n42vE5WqL+qboX6UWNjVbKvTRBMYvsM\nEIniEyPb+pa0gyrUb0k5lkrtgiDbZ4CIFLxFxt2PhhniA2G5XLraBcEekiVABArMDpvZITN7Iry8\nJWlG6yvUv6busnfU2LpK9sOUWk9wEpLJkuOL19WJfE8S0ChD/r395tYVrdzqv903+DMknv+7uvfz\nJOkhSf8g6UeSdml9hfq0YyPVdiYokSwBKjdk5rdp6sN614f/ZPUx8EPN5iXN9WZ/ofr8Z8P4z919\nqVeRPlmhPu1Y2r9C7WeCPclkyczUZmaEwCTk+Hfm7qckneobWxjwvnUV6tOOpVHrmWAPyRKgIpwY\niQvJEmDCODscJ06WABPSuS/dI2KNC4IkS4AJYiYYN06WACXjnmC8SJYAE8BMMH4kS4ASMROsDwqz\nAiVgJlgPFGYFymFmqR4xa0UQTGL7DFAc61iqR8xaFQTZPgMUi5lgjbF9BsiPIFhTbJ8BikEQrDmS\nJUA+BMGGIFkCZGQpHxFrfRAkWQJk1+l0Uj2GGdB3+FB4HEmM0Xd4UkiWAOPJsxzu7zscXr8UCqtu\nM7NPRdt3OE+0jhXJEmB8eYLggL7D23SvX/ByeB1f3+GM0XqsyFwlkiXAGAq8J+juC4kS+bPqNlmK\nsu9wlmg9VmSOAWeNgdHKyA6HCdTiOM2S8hg7COaM1rXAWWMgnZK2yOx396fD83j7Diejdez7gPKg\nrzEw3LB/+3f/5z/09i//M9Wn6Pt8h9z9aHi+X9ILKrnvcJ7i/6OitStFZD6z8Ozq84dn92r73L4c\nl1Ss3rnv3ozwwO5pHV+8rk1mivxMOLDGpcVXdPnVc4V/3mHFETZ/6KPa/KGPrr5+c+nk+o9N9B12\n99Mh6B0xs7+SdL+kp8Ika5f19RNOO5bq7+AZtoKEaL0Qnu9XN8DtcvdjIRO8GpmTY/0XZmb+/Plr\nY3/9qrz59ooO7J7WiaUburviBELU1oE903L3XL/BZua/e/DE6DdK+t9jT+X+emXJmh0+YmY/N7Ob\n6maAC+0IHyuSJcBaTTg2N/ZyOGx5WZfkKLIjfKySyRJmhMDwe4J1womRDDhrDAScHW4fzhoD9zRh\nOUwQzIGzxmg7gmCLcdYYoMcIxMkStBszQawiWYI2IghCEskStBdBEOuQLEGbEASxBskStA77BDEI\nyRK0BTNBbIizxmi6TsdSPWJGECwJhVnRBswEkQrbZ9BUZukeMSMIloztM2iyJswE81SWxpiS22dm\npjZTqh+1l+dXOBRbviJpKlGkeV7dfiGz7v7MOGNZMROcELbPoImyJkZCsWV399OSHjKzaZtAo/WB\nf4c8H4zxkSxBk+S4J/hpdVvxSt3Z4Kc1gUbrgxAEK0KyBE2QY4vMTa1txPaQ0jdVL7SdL0GwAiRL\n0BQ5EiMn1Q18Cn++NrGL7kNipGIkS1BnwzK/v762pF9f+/HQj3P3q2b2Qri/d0fdZe2URrfuTd3O\nNy1mghUiWYK6G3YP8P0zO/T7f/znq4/1H2c71W3Je1HSlpAgeVHdxukKf35X3ebracYyIwhGgGQJ\n6irrcjgEv1tm9oSk5xJjvczx0Na9RbfzZTkckWSy5PjidZbGiF6eX9Ew++sfS9W6t8h2vswEI0Gy\nBHVEAQWUgsKsqIsmHJsjCEaGZAnqhAIKKA3JEtQBM0GUjsKsiBkzQZSKwqyIHTNBTAxnjREjZoKY\nCLbPIFZskcHEsX0GMWE5jIli+wxi0/ogGMpj957Pm9n+UWPIj2QJYtHqe4Lh8PLj4XnpJbCxHskS\nVK31M8GE0ktgYy2SJYhBa2eCZrYzzPB6tqrkEtgYjmQJqtLmmeD9hV4FMiNZgipt6liqR8zGDoJh\nFni2b/iO1pfA7h/LVQIbGyNZgirkWQ6H3MF86CHcG0uVYC0y6ZplJrjNzJ4ws0OSpkLC41squQQ2\n0iFZgknKuRx+2t1PSZoxsx0pEqylJF3HriwdLlohCG4JY0tmtqu/3PWgsX5nFp5dff7w7F5tn9s3\n/t8CkrqzQenejPDA7mkdX7yuTWaKfEWCkl1afEWXXz1X+OfN+nsVZn8XJMndj4axI5L+Obyll2Cd\nGjHWS7pmLrGfuby+uy9IWki8zlQC+zOHvpT1ErAButghafvcvjUTjG8f+0YhnzdH0uMRSR5mdY+7\n+zNKn2Cl7zA2RrIEk5Jzi8zNRHOleXVbaU4cQbDBSJagbJbyfwPcVHcpK3WTqI8oXYK18KQr3eZa\noL+LncQ9QhRj2PaX1y4v6rXLixt96ElJvazwVnXvDy5L2iXprLrJ1JfCf087lglBsOEGJUtOLN3Q\n3RUnECK3YUvdD2yf0we2z62+/tl3Ftb8d3e/amZ3wjL4gXBPcGAyNWvSNS2CYIvQ1xhFy/M7lEic\nnhowNuh9G45lxT3BluCsMcrQ2rPDqDfOGqMobT47jJpi+wyKxEwQtcX2GRShY5bqETOCYMtx1hh5\nEARRayRLkFfvd2jUI2YEQUgiWYJsSIygEUiWICsSI2gUkiUYFzNBNBLJEqTFPUE0DskSjIOZIBqN\nZAlG2WSW6hEzgiAGIlmCNEiMoPFIlmAjLIfRGv3JEgIhJGaCaIlByRJmhJDyHZsLfYP3hy5zvbFa\n9B1Gi7F9BklZZ4KhKvSToXfwbJV9hwmCSI3tM+iXNTvs7i+7+xfDy5lQIv+z6jZRku71GB411us7\nnBlBEJmwfQZS/sRIWM5+Pryk7zDqge0z6Ml7YiQ0WPqCmW2Z2EX3IQgiM5IlyBoEw7293r28ZUmH\nJd0WfYdRR3Sxa69hS90bPzmvGz+5sNGHPi6p15i413f4JdF3GHUzqK/x8cXr2mQ0eG+DYT/jBz+2\nRw9+bM/q63Pf+mb/W56T9Gdm9pC6md7TEn2HUXPJZMnM1GZmhC2Q9Ufs7m9IytxjmL7DiA7Jkna6\nzyzVI2YEQRSKZEm7cGwOGIKTJe1AtzlgAE6WtAczQWAETpY0G+X1gQ2QLGk+lsNACiRLmqu1y+Fw\n5GXezOYTY6XW/EL9UZi1edrcY+Rpdz8laSZlHbDcNb9QbxRmbaZW3hMMs78LkuTuR1PWActd8wvN\nwfaZ5mhlEJT0iKSpMNvrLXPT1gFDy7F9plna3GjpprtflFZnhvwKY2xsn6m/JswEsxRQuKnu8lbq\nLncf0eCaX64UNb/OLDy7+vzh2b3aPrcvwyWhbpLVZ+Y++D796Uc/SPWZEl1afEWXXz1X+OeNfJKX\nSpYgeFJSLyvcqwO2rIw1vz5z6EsZLgFNkVwaH9g9rRNLN3R3xQmEBds+t2/NBOPbx75RyOeNfQ9g\nGmMvh939qrrZ3nlJD7j76USNr9X6XoPGirxwNAvJknra1En3iFmmeoKJWl6nBowNeh8wFIVZ66uj\n7D+gMJG6I2k29BqpROQxGm1DsqRecvQdjmYfMUEQ0eCscf3kyA5Hs4+YIIjocLKkPnIUUBi0t7gS\nBEFEi2RJ/FpbQAEoGydL6iHHTLC/x3Cu3sF50G0O0aOLXbw2DflRXFp8RZcWN9yc/aKkORXUOzgP\ngiCixsmSuA07F/yRXY/qI7seXX195tjfrfnv7n7RzOZi2EdMEEQtcLIkTnm+/bHsI+aeIGqFwqxx\nobw+MEEUZo2PpXzEjCCIWmL7TBzYIgNUgO0z8WhzjxEgCpw1rlabK0sDleOscfW4JwhEgGRJdZgJ\nAhEhWTJ5nZSPmMV+fUAqJEuqwUwQiBDJksnhniAQGZIlk8UWGSBSJEsmg83SQORIlpTLUv4vZgRB\nNBbJkvIxEwRqgmRJOTqyVI+YEQTReCRLysNMEKgRkiXFKzMImtn+8DiSGJsPY18ed2wYgiBah8Ks\nxSlri0wou/9kaM4+a2Y7+hq23zaznSnGRjZ2JwiiVSjMWqyyssPu/rK7fzG8nAk9SJIN26+q27B9\n1NjIxu4EQbQW22fyK/ueYFjOfj68HNSwfUvKsaEIgmglts8Uo+x9gu7+jKQvmNmW4q56LbrNofXo\na5zdsG5/F8//m5Yu/GDDjzWzQ5J6/7dj4fmyu59N3NdbUndJe1jrG7a/Fj5m1NiGjd0Jgmg1+hrn\nM2yWN7vnMc3ueWz19T/9/dfXvcfdFzb41I9LWgzPt0q6oG6D9l1a37A97dhALIcBkSzJqsR7gs9J\n2tabLbr76V6D9mTD9rRjG30hZoJAQjJZcnzxOkvjEcqqEOPub0ha15x9UMP2tGPDMBMEApIl46Oe\nINBQnDVOqQFRMFMQTBxJOThgLNPRFSAWnDVOr5WltELqejkcSbma8jjLyKMrQGxIlozW5gIKXwt/\npj3OMvLoChArTpYM14DV8PhB0N0vSlo2s1u6dzQl7XEWoFZIlozQgCiYZTm8Rd2d21+RtGBmM4Vf\nFRAhkiXrdcxSPWKWZZ/gYUlfdfc3zGxZ0pNKd5xl4NGVMwvPrj5/eHavts/ty3BJQLnqfrLk0uIr\nuvzqucI/bw3+6iOZj/n/aGb2l5IW3P318PqgpB9J2uXux0ImePXoSnKsf+e2mfnz56/l/TsAE7Xi\n0lvvrOjA7mmdWLqhuytei0CYdGDPtNw911WbmS9eez3Ve+emt+T+emXJck/wqKTDZvaEmR1092NF\nHF0B6oTCrF1N2CKT6dhcKG/TP5br6ApQF8ml8fMXrtV6RphX5Lf7UuHECJBD27fPNCA5TBAEsmL7\njGRmqR4xIwgCBWjr9pk2nxgBELT5rDHLYQCrWnnWuMQoGGoQzJvZfGKMvsNA7NqULCl5i8zT7n5K\n0gx9h4GaaFuypKx7gmH2d0Hq7k2m7zBQQ21IlpS4Gn5E0lSY2fWWtPQdBuqiLcmSkrfI3AxVq3oz\nw1K+ezRaAkqUXBo38WTJsPh2/off04Uffn/Exw7vO6xuwZXl8N/uqDszvCP6DgP11N/FTqpH9ZlR\nhv0V9j76ce199OOrr7/5t19Z954RfYdPSuplhXt9h5dF32GgfgYlSxqzfaakm4LuflXdzO68pAfK\n7Ds8dimtIlFKC23SC3pb3n3faj3Cd1U0HSyqlNbPfvlmqvf+0e+9tzmltABk08TtMxybAzC2jknv\nfVenEUtjjs0ByKwJhVmpIgMgsyYkS1gOA8itzmeNWQ4DyKX2yZIGREGCIBCJOp41bkKjJYIgEIG6\nnjXmniCAQtUtWdKA1TBBEIhRXZIlbJEBULg6JUtYDgMoVezJEpbDAEpTh2QJM0EApYs5WcIWGQAT\nE+VZ4washwmCQA3EWpi1zBhoZl82sydCGf7eGH2HgbaLaftMxyzVY1yhKrS7+2lJD5nZNH2HAcS3\nfaa8qeCnda/R0pXwmr7DAO6JYftMicvhm1rbMe4hpe8xTN9hoOli2T5T4haZk+oGPoU/Xyvqmvul\narlpZjt7TZDD63l1p5uz7v7MOGMAilN1X+Nh21/+7Xv/oh98/183/tgN+g67+1UzeyHc37uj7rJ2\nSlX0HQ43KJ+T9Ifh9epNRzObCa81Ymybme0Y1foOQDb9fY2zJCOyGPZlHvvEJ/XYJz65+vrrX/2b\nde/ZqO9wiCFz7n7MzA67+2kzuyppTpPuOxwyLFcSQ6XcnASQTZXJkrKWw2HlecvMnlB3EtYbK7zv\ncKrlsNbe29yqEm5OAsgvmSyZmdpc+oywzNMgYXtM/9ixrGPDkBgBGqKKZEmbzg4nv423tf5G5J0U\nYxvenARQjEmeLGnAqblMy+EXVeDNyTMLz64+f3h2r7bP7Ut5SQA2kkyW/PXCC7py8VzxXyT2CJeC\n+YhNlmGbyz9KOtRbo5vZQXUTIDO9tXfasb7P7c+fv1bc3wbAGisuvfXOig7sntbxxevaZKaOSQf2\nTMvdc4UwM/NfvbWS6r3vf08n99crS5rs8Cl3n0repHT3Y+7+cjKwpR2rs0uLr1R9CanV6VolrrdM\n1398vrSTJW26JwhJl18tYTlRkjpdq8T1lqVj0s8unluXLCkKQRBALfQnS4rShKKqaRMjpfnI77y/\n6ktI7Qe/tbk211una5W43jL1rrU3//uvX75e2OeOfZaXxsjESKlf3CyG2rhAKxWQGLkmKe208rq7\nT+f5emWpNAgCQNW4J9gQRZYbR70kCpb0Xhdegr7JKguCsf4AzOxQeBxJjEX9CxQOij8enucuN17y\nte4M37v5xFiU39/ENRzc6LqqvNbwsz+ReF1KCfomqyQIxvoDCL9QL4USP9vM7FM1/AWKvaLP0+5+\nStKMme2I9fsbrmE5XMPVWK+VKk/5VTUTjPUHsE33rmU5vI76FygUvH05MZS2ys/EhdnfBUly96Oh\nxFHM39+vhT9nIr9WqjzlUFUQHPSDqpy7LyROt8xK+pEiDirB/RV+7XE9ImkqzKB6S8cov7+hdt2y\nmd1KXEuU14p8SIwMEJY4i7FXwg6zwLN9w7FX9LmZKI45r7UViqJhZlvUrZj0FUkLZjZT8SVthCpP\nOVS1Wbr/BxXbD2C/uz8dng/6pRqrh0GJtoV/nFPqzrB2SPqWcpYbL9FN3WujeEfdmeGgf6AxfH8P\nS/qqu79hZsuSnlS8vwulVXlqg6pmgi+q+41X+PO7FV3HOmZ2yN2Phuf7Jb2g9dc6aGziQnGLXmGL\nLWEsd7nxEp3Uve/bVnXvD8b6/XWF4BK+x7eHXFel1xpm03OhDH1pJeibrLLN0qPKbFUh/EK8qO4v\n/P2SnnL3s1nKhGGw8H27LWlXb7Yd6/c33Le8IumBPCXjEDdOjABoNRIjAFqNIAig1QiCAFqNIAig\n1QiCAFqNIAig1QiCAFqNIAig1f4f0wupeoHVDfMAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.imshow(H_sampled, interpolation=\"none\"); plt.colorbar(); plt.show()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "E[0] = -95.5531875434\n", "E[1] = -95.5531875434\n", "E[2] = -86.7633541611\n", "E[3] = -86.7633541611\n", "E[4] = -78.1354749183\n", "E[5] = -78.1354749183\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAENCAYAAAABh67pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4U1X6wPHvTdJ9p9BS9pYdCrRQkEWhUtxwGQUFf+4L\nCIyMOoqKjArOjKMOKs6ooyLLqDMjiyLuCrQWZNNCW7CsBUpLoWzdt7RNcn5/3CS0TbrRtGnT83me\nPm1uzr33BNo3J+ee+76KEAJJkiTJdWmc3QFJkiSpZclAL0mS5OJkoJckSXJxMtBLkiS5OIcEekVR\noms9ftX8fXa1bdMVRYlTFOVpR5xTkiRJapxmB3pFUeKA9bU2P6IoSjpw3NwmGhBCiHigQFGUqOae\nV5IkSWqcZgd6c/A+XmvzLCFEfyFEgvnxTKDA/PMJYEpzzytJkiQ1TkvN0UfUmqYJBPKqPR/cQueV\nJEmSatG1xEGFEK8DKIoyxTy1I0mSJDmJw0f0iqLMVhRlmvlhHhAO5AOdzNsCgVxHn1eSJEmyz1Ej\neqXaz0mo8/AAfYH3gT1ADJAARACbbQ6gKDIXgyRJ0mUQQij1Pe+IVTfTgVGWUbwQIhWYad5+TAiR\nat5mWaGTb3lsp7Mu+7V48WKn90G+Nvn65Otzva/GaPaIXgjxOfB5rW0f2mm3ornnkiRJkppO3hkr\nSZLk4mSgbyWxsbHO7kKLceXXBvL1tXeu/voaQ2nsHE9LUxRFtJW+SJIktReKoiAauBjbIuvoJUmS\nmqpPnz5kZmY6uxttVu/evTl58uRl7StH9JIktQnmkamzu9Fm1fXv05gRvZyjlyRJcnEy0EuSJLk4\nGeglSZJcnAz0kiRJLk4GekmSJBcnA70kSVIDCgsLmTFjBqNHj66xxLGwsJCYmBhWrKiZ4SUhIYGG\nFBYWEh8f7+iu2iUDvSRJUgMCAgKYM2cOERER9OnTx7o9Ly+PFStWMGvWLOu2pUuXMnny5EYdMzg4\nmJSUlJbocg0y0EuS1C4oiuKwr8sRERHB8eM1q6ampKQQFRVV43Hfvn0bfcyoqCjWrl17Wf1pChno\nJUmSGiE8PJwTJ05YH8fHxzNt2rQabdauXWuzLT4+nvj4eBYuXEh8fDwzZ86s8byiKBQVFbVcx5GB\nXpKkdqK1c7jbYwnKGRkZBAfblr4uKCio8biwsJDRo0cTFxdHfHw8cXFxvPbaazXaxMTEsGfPnsvu\nU2PIQC9JktRIMTExHD9+3GbKxiI/P7/G44CAAPz9/cnIyCAmJgagxhw/QGBgoM0bhKPJQC9JktRI\n4eHhLF++3GZ6xqJTp041HhcWFlJYWMhnn33GqFGjAGwuvhYUFBAYGNgyHTaT2SslSZIaKSYmhoiI\niDqfrx2wly9fTlBQEH379iUvL4+EhATryN4iKSmJRYsWtUh/LTpc9sqTJ0/y0UcfERwczPz581v8\nfJIkNY4rZK9MTU3l+PHjTJ8+vdH7PPfcc7zyyisNtpPZK5vg/PnzLFmyhNWrVzu7K5IkuZioqCgy\nMjIa3T4lJcVmFU5L6HCBvnv37gBkZ2c7uSeSJLmiBQsWNPrO2Pz8fLsXdR2tw03dGI1GPDw8MBqN\nVFRU4O7u3uLnlCSpYa4wddOS5NRNE2i1WsLCwgA4c+aMk3sjSZLU8jpcoAfo0aMHIKdvJEnqGGSg\nlyRJcnEdMtDLC7KSJHUkHTLQyxG9JEnN1dhc8m0hN32HDvSnT592ck8kSWqvGpPuuK3kpu/QgV6O\n6CVJuhwpKSkNBvC2lJteBnpJkqQGWKZWFi5cCFzKUrl06dI6p2baUm76DpnUrFu3bgDk5ORgMBjQ\n6TrkP4MktStPPAGpqc0/TlQUvPVW0/axpBlevnw5gPXGpTvuuIOCggIKCwtZvnw5ffv2JTw8nOjo\n6Dpz0/v7+7Nw4UJeffVVmxG/JTd9Y6Z7mqJDjujd3d0JCQnBaDRy7tw5Z3dHkqQ2LioqinXr1llH\n4JaiI5YUBsuXL2fOnDlMmzbNOv2Sl5dX4xjOzE3fYYeyPXr04Pz582RnZ1uXW0qS1HY1dRTuaFu2\nbOHDDz8kJSWF6OhoCgsLrRdkk5KSePrppwGsSc3s5aYHbHLTR0dHW9sUFBQQFBTk8L536ECfnJws\nV95IktQozz33HGvXrq0R4C0JySxz6/7+/tYRee2A7czc9B060IO8ICtJUuNERUURFRVFQkKCTcbJ\n0aNHk5eXh7+/v7X4yIwZM/j888+tuektI/76KIqCv7+/w/suA70M9JIkNVJhYaHdClOzZ8+2jtif\ne+45AKKjo5t0E1RL5qZ3SKBXFCVaCJFS7fF0oAAYKYRYWtc2Z5JpECRJaqqAgAACAgLsbrc3Yrfk\npm9oFY0lN72jV9tYNHvVjaIoccD6ao+jASGEiAfyFUWJrrWtQFGUls+03wA5opckqTU09s7Ylgry\n4IBAbw7ex6ttmok6cgfIAKbU2nbCvM2pZBoESZI6Cketo6+e9CEQqL6ANBgIsLPNqapP3ciqNpIk\nubIOecMUgI+PD0FBQVRWVnLx4kVnd0eSJKnFOCrQVx8S5wOWOwUCgYuo0zbVt+U66LzNIufpJUnq\nCBy1vLL61M06YBSQAEQAm83bY+xsq2HJkiXWn2NjY4mNjXVQ9+zr3r07v/32G9nZ2TXuTpMkSWqr\nEhMTSUxMbNI+zQ705mWToxRFmSaE2CCESFEUZZR5NU6+ECLV3C6m9rbaqgf61iBH9JIktTe1B8Ev\nvfRSg/s4YtXN50KIYCHEhmrbVggh4oUQK+rb5mwy0EuSdLlkhal2olevXgBkZWU5uSeSJLU3ssJU\nOyEDvSRJl0NWmGpHevfuDUBmZqaTeyJJUlsmK0y1Y9Xn6I1GI1qt1sk9kiSpXnWtxKtrFYq99k1c\nsQINV5gqKioiKSmJ5ORka84bWWGqjfD09CQ0NBSj0ciZM2ec3R1JktqohipM+fv7ExERUaOqlGXU\nbyErTDlR7969OXfuHFlZWfTs2dPZ3ZEkqT5NHY1fxui9LvVVmLKnLVWY6tAjepAXZCVJahxLhan4\n+HhSU1PJyMiwKUBSnaUAicXy5ctZv369dbomISHBZuomKSnJpuqUI8gRvbwgK0lSI9RXYcqieoLE\nmTNntpkKU3JEL0f0kiQ1Ul0VpkCd2klJSeHkyZOA+sZgKRTeGC1ZYUoGenOglyN6SZIaEhAQYHMB\n1WL27Nn8+OOPNZ63VJhqiKXCVH1TQc2htJVc7IqiCGf0JSUlhZEjRxIZGclvv/3W6ueXJEmlKIqs\nDVGPuv59zNvrvU1Xjujl1I0kSS6uwwf6Tp064ePjQ1FRUYusX5UkSXK2Dh/oFUWRo3pJklxahw/0\nIC/ISpLk2mSg59JaejmilyTJFclAj7wgK0mSa5OBHnl3rCRJTScrTLUzckQvSVJTyQpT7YwM9JIk\nNYWsMNUOde/eHY1Gw5kzZ6isrHR2dyRJamNkhSkX4ObmRrdu3cjOzub06dOEh4c7u0uSJNXyww8/\ncPbs2WYfp2vXrlx//fVN2qehClOFhYXs2bPHml8+Li5OVphqiywXZC2Z5yRJkiwaqjC1bt06+vbt\ny4IFC3jttdcAalSbAllhqk3o06cPO3bskIFektqopo7CHa2+ClOzZ88Gas7Lt6UKUzLQm1lyTJ84\nccLJPZEkqS2yVJiqHuBrpxVet26ddURfO2AvX76coKAg+vbtS15eHgkJCTbVpJKSkli0aJHD+y4D\nvZllXr4phQIkSeo4Gqow9fnnn7Nw4UJyc3Px9/eXFabaIhnoJUlqSF0VpixBfsaMGdYLtm2pwlSH\nLzxikZWVRe/evenatSs5OTlO64ckdVSuWngkISGhwVU0hYWF7N27t952zSk8IgO9mdFoxMvLi6qq\nKsrKyvDy8nJaXySpI3LVQO8ossKUA2i1WusdsnLljSRJrkQG+mrkyhtJklyRDPTVyAuykiS5Ihno\nq5GBXpIkVyQDfTVy6kaSJFckA301ckQvSZIrkoG+muqBXi7zkiSpPrLCVDsVHByMn58fRUVFNpnn\nJEmSquvwFaYURXnV/H12tW3TFUWJUxSl4YQPTqIoipy+kSSpQbLClOoRRVHSgeMAiqJEA0IIEQ8U\nKIpimxGojZCBXpKk2mSFKftmCSE2VHs8E9hk/vkEMAVIbaFzN4tl5Y0M9JLUBm2Jtb99SmLj29fV\nth6NrTB14sQJ+vbty+TJkztEhamIWtM0gUD1Se/gFjpvs1lG9HKJpSRJFg1VmNqzZw+KohATE0Ny\ncjLQASpMCSFeB1AUZYqiKHGN3W/JkiXWn2NjY4mNjXV43xoip24kqQ1r6mj8MkbvdamvwlRcXBwZ\nGRn8/e9/57333gMuvRlYOKrCVGJiIomJiU3qu8MDvfkCbK556iYPCAfyAUtdrUAg196+1QO9s8ip\nG0mS7GmowlR4eDjPPPMMM2bMYN26dQQGBtbY31EVpmoPgl966aUG+94SI/ok1Hl4gL7A+8AeIAZI\nACKAzS1wXoewfJQ6efIkRqMRrVbr3A5JktQm1FdhauHChcydO5fw8HDrIHHGjBmuW2FKCJEKzFQU\nZTpwTAiRat6GeRon3/K4LfL29qZbt25UVVWRlZXl7O5IktSG1FVh6s477+T48eMsXbrUWjM2Ojpa\nVpiqzdmFR6qLjY1l69atbNq0iWuuucbZ3ZGkDsFVC4+0hQpT8s5YO/r16wdAenq6k3siSVJ719g7\nYx29pLI6Gejt6N+/PwDHjh1zck8kSZKaTwZ6O+SIXpIkVyIDvR1yRC9JkiuRF2PtKC0txdfXF3d3\nd8rKyuQSS0lqBa56MdZR5MVYB/Px8aFbt25UVlZy6tQpZ3dHkiSpWWSgr4Ocp5ckyVXIQF8HOU8v\nSVJ9ZIUpFyBH9JIk1afDV5hyBXJEL0lSXWSFKRchR/SSJFnIClMuyhLoT5w4IbNYSlIb8ET2GVLL\ny5t9nCgvL97q0a1J+zRUYcpi6dKlzJkzB39//w5RYard8/HxISwsTC6xlCSpwQpToAbypKQka2Up\ny6jfwuUqTLmK/v37k5OTw7Fjx2z+QyRJal1NHYU7Wn0VpgD27NnDmDFjrI87depUY39HVZi6HHJE\nXw/L9M3Ro0ed3BNJkpzNUmEqPj6e1NRUMjIyrKP5lJQUYmJiaty5aq/C1Pr1663TNQkJCTZTN0lJ\nSTZVpxxBjujrMWjQIACOHDni5J5IkuRs9VWYysjIIC8vj6SkJIKCgpg1axYzZ8503QpTrsQS6A8f\nPuzknkiS1BbUVWFq2rRpxMTEWKdnQH1jkBWmamlLSc0s0tPTGTBgAL169SIzM9PZ3ZEkl+aqSc3a\nQoUpGejrYTAY8Pb2pqqqipKSEnx8fJzdJUlyWa4a6B1FZq9sITqdznqHrLwgK0lSeyUDfQMGDhwI\nyHl6SZLaLxnoGyAvyEqS1N7JQN8AGeglSWrvZKBvgAz0kiS1d3LVTQMKCwsJDAzE09OTkpISmdxM\nklqIXHVTP7nqpgUFBAQQFhaGXq8nKyvL2d2RJKmNkBWmXIycvpEkqTZZYcrFyEAvSVJ1ssKUC5KB\nXpI6tvZeYUoG+kaQgV6S2ojYOr6a0v4yWIqFWJKUVa8w1alTJ+vz8+bN4+TJkwB1VpiKi4sjPj6e\nuLg4XnvttRptLBWmHE2mKW4ES6A/ePCgk3siSZIzREVF8eGHH9qtMBUdHU1GRgYJCQk1UgzbqzAF\nyApTbVXPnj3x8/Pj4sWLnD9/npCQEGd3SZI6psQWbl+PhipMrV27FkVRiImJISoqSlaYam8URSEy\nMhKAtLQ0J/dGkiRnqK/CVHh4OLNnz2bWrFm8//77gKww1S5FRkaya9cu0tLSHF6hXZKktq++ClOW\naR1/f3/rlM2MGTPaTIUpGegbSY7oJUmqq8LUNddcQ1JSEhkZGdYLrNHR0U26CcolKkwpijIdKABG\nCiGW2nm+TaZAsEhISCAuLo5x48axc+dOZ3dHklyOq6ZA6DAVphRFiQbChRAbFEWZDSQJIVJrtWnT\ngf78+fOEhobi5+dncxFGkqTmc9VA7yjtIdfNTNTRPMAJYEornddhQkJC6NKlC8XFxZw6dcrZ3ZEk\nSWq01gr0gUBetcfBrXReh5Lz9JIktUfyYmwTREZG8tNPP5GWlsbUqVOd3R3JgYQQpKWlsWvXLlJT\nUzly5Ajnzp0jL08dn7i5uRESEkLPnj0ZPHgwV1xxBePHj6dz585O7rkkNay1An0+YLl7IBDItddo\nyZIl1p9jY2OJjY1t6X41iRzRuxYhBMnJyaxatYqvvvqK7OzsettnZWWxZ88evvjiCwA0Gg3jx49n\n+vTp3Hvvvda7JSWpJSUmJpKYmNikfVrzYuwoIcQKRVGeBja3t4uxADt37mTChAlER0eTnJzs7O5I\nl8loNLJmzRpef/11UlMv/RqGhoYSFxfHyJEjGTp0KGFhYQQHB6MoCpWVlZw9e5asrCxSU1PZvXs3\nO3bsoKqqCgAPDw9mzJjBwoULGTJkiLNeWrvWp08fMjMznd2NNqt3797WPDrVNeZiLEKIVvkCZgFx\nwKw6nhdtXUFBgQCEh4eHMBgMzu6O1EQmk0l8+umnYsCAAQIQgOjUqZN44oknxJ49e4TRaGzS8QoL\nC8WaNWvEDTfcIBRFEYBQFEXcddddIj09vYVehdSS3nvvPQGIBx54wNldaTRz7Kw//jbUoLW+2kOg\nF0KInj17CkAcOXLE2V2RmiAtLU1MmjTJGuAjIiLEypUrhV6vd8jxT5w4IX7/+98LNzc3AQh3d3ex\nePFiUV5e7pDjS61j3rx5AhBvvPGGs7vSaI0J9DLXTRNZ5ul/++03J/dEagyj0chf//pXoqKi2Lp1\nK507d2b58uUcPnyYhx56CA8PD4ecJzw8nHfffZf09HTuu+8+Kisreemllxg+fDi7d+92yDmklrd/\n/34Ahg8f7uSeOJYM9E00YsQIgBpzu1LblJmZSWxsLC+88AIGg4G5c+dy5MgRZs+ejZubW4ucs3fv\n3nz00Uds3bqVIUOGkJ6ezpVXXsnLL7+M0WhskXNKjiGEkIFeUllSispA37b9+OOPjBgxgu3btxMW\nFsbmzZt57733bFLHtpSJEyeSnJzMU089hdFo5Pnnn+f666+3LteU2p7MzEyKi4sJDQ11uVTkMtA3\nkSVrXUsV8ZWaRwjBsmXLmDp1KoWFhdx8883s37+fKVNa/2ZsDw8PXn/9dX788UdCQkLYsmULV1xx\nBYcOHWr1vkgNc9XRPMhA32T9+vXDx8eH06dPc+HCBWd3R6rGYDAwZ84cnnzySUwmEy+++CIbN250\n+k1N1157LUlJSURFRXHs2DHGjh3bpKyGUuuQgV6y0mg0cp6+DaqoqGDmzJl8+OGHeHl5sW7dOl56\n6SU0mrbxK96rVy+2b9/O7bffTlFREVOnTmXDhg3O7pZUjQz0Ug1ynr5tKSkp4aabbmLDhg0EBASw\nZcsW7rjjDmd3y4aPjw9r167lD3/4A5WVldxxxx2sWrXK2d2SzPbt2wfAsGHDnNwTx5OB/jLIefq2\no6SkhOuuu44tW7YQEhLC1q1bGT9+vLO7VSeNRsM//vEPlixZgslk4uGHH2bFihXO7laHV1RUxNGj\nR3F3d2fo0KHO7o7DyUB/GSwjehnonau8vJxbbrmFnTt30rNnT7Zv326dVmvLFEVh8eLFvPnmmwA8\n8sgjfPLJJ07uVcdm+XQ+bNgw3N3dndwbx5OB/jIMHToUrVbLkSNHKC0tdXZ3OqSKigqmT5/OTz/9\nRFhYGAkJCfTv39/Z3WqSP/7xj7z66qsIIXjggQdYt26ds7vUYVlyV40cOdLJPWkZMtBfBk9PT4YM\nGYIQQt4h6wRGo5F77rmH77//ns6dO7Nlyxb69evn7G5dlmeffdY6jXP33XezadMmZ3epQ9q7dy8A\no0aNcnJPWoYM9JdJXpB1nqeffprPPvuMgIAANm3a1O6zRb744os8/fTTGAwGpk+fLn+nnMAS6OWI\nXqpBXpB1jn/+858sW7YMNzc3vvjiC+sbbnumKAqvvvoq//d//0dJSQlTp04lKyvL2d3qMEpLSzl8\n+DA6nc4lV9yADPSXzRJgZF761rNx40aeeOIJAFauXMnVV1/t5B45jkajYfXq1UyaNImcnBymTp1K\nQUFBwztKzZaamooQgsjISDw9PZ3dnRYhA/1lGjlyJIqisG/fPioqKpzdHZeXmprKXXfdhRCCP//5\nz9x7773O7pLDeXh48MUXXzB48GAOHDjAXXfdJROhtQJXvxALMtBfNn9/fwYNGkRVVZX1jjqpZVy8\neJFbb72V8vJy7r//fp5//nlnd6nFBAUF8d133xEcHMz333/PCy+84OwuuTxXvxALMtA3y5gxYwD4\n9ddfndwT12UwGLjzzjvJzMxk9OjRvP/++yhK/VXT2rs+ffqwdu1atFotr7zyCuvXr3d2l1yaq1+I\nBRnom2X06NGADPQtaeHChcTHxxMSEsKGDRtcdg61tri4OF5//XUAHnjgAfmpsYWUlZVx6NAhtFqt\nS+a4sZCBvhksI/qkpCQn98Q1rVmzhjfeeAOdTsdnn31Gjx49nN2lVvX4449z7733UlZWxm233UZh\nYaGzu+Ry9u7di9FoZPjw4Xh7ezu7Oy1GBvpmGD58OO7u7hw+fFj+ETpYeno6s2fPBmDZsmVcddVV\nTu5R61MUhQ8++IDo6GhOnDjBrFmzLPWVJQexlHm84oornNyTliUDfTN4eHgwYsQIhBDWeT6p+Swp\nh0tKSpg5cyaPPvqos7vkNJaUy35+fnz22We89957zu6SS7EE+rFjxzq5Jy1LBvpmktM3jvf000+T\nkpJCREQEy5cvd/mLrw3p168fH374IaDmx5F3zjqOHNFLjSJX3jjWxo0befvtt3Fzc2Pt2rX4+/s7\nu0ttwsyZM5kzZw6VlZXMmDGD4uJiZ3ep3cvOzubMmTMEBgYyYMAAZ3enRclA30yWlTe//PKLk3vS\n/mVmZvLggw8C8Pe//52YmBgn96htWbZsGcOHDyc9PZ158+Y5uzvtXvXRfFupRNZSXPvVtYKBAwcS\nGBjI6dOnOXXqVKufXwi4cEH93p4ZjUbuvvtuCgoKuPnmm3n88ced3aU2xzJf7+3tzX//+1/WrFnj\n7C41X2mp+uUEHWXaBmSgbzaNRsO4ceMA2LFjR6ueu6oKbr0VQkJg3Dhozwt/li5dyo4dOwgLC2P1\n6tUdfl6+LgMHDmTZsmUAzJs3j+zsbCf3qBm+/hpCQ6FLF3DCTWGWT+GufiEWZKB3iCuvvBKA7du3\nt+p5ly6Fr76C+++HPXvgj39s1dM7TGpqKi+++CIAq1atIjg42Mk9attmz57NTTfdREFBAQ888AAm\nk8nZXWq6CxfgnnugXz8YMgQeeABa8U2rqqqKPXv2AJeus7k0IUSb+FK70j4lJiYKQERFRbXaOUtL\nhejcWYgbb1QfP/mkEBqNEMePN/1Yf845KyYcSRc7S0oc28lGKC8vF5GRkQIQ8+bNa/Xzt1dnz54V\nXbp0EYBYtmyZs7vTdIsXC6EoQhw8KERGhhBLlghRUNBqp09KShKA6N+/f6uds6WYY2f98bWhBq31\n1Z4DfWlpqdDpdEKj0YjCwsJWOefq1er/XmKi+jg7W/27efHFph1nQ36BIHmf0CTvE533p4kCg8Hh\nfa3PggULrH9wJU54o2nPNm7cKADh4eEh0tLSnN2dpiksFGLDBqed/s033xSAePDBB53WB0dpTKCX\nUzcO4O3tzahRozCZTNYLPC1t/Xro0wcmTlQfd+8OU6bAxx837cLsG+cv0NfdnV0D+nHRYOS9C7kt\n0l97tm7dyhtvvIFGo+Hjjz/Gx8en1c7tCn73u9/x8MMPU1FRwT333ENlZWXrnNgEJADNWQDg7w+3\n3eagDjXdtm3bAJho+QNycTLQO8iECROA1rkgW1gImzfD9OlQ/ZrlnXfCyZPQ2PxXZ6uq2FFaxkPB\nnRjj481EXx8+ystvldvsi4qKuP/++xFCsGjRog5xQawlLFu2jIiICFJTU3n55Zdb56RbgDjga/Pj\nM8Dm1jm1Iwgh+PnnnwEZ6KUmckSg/7GomOfPnCXXYKi33U8/qStubrml5vbrrjMf58fGnW9LcQkA\n1/v7AvB/QYEcrqhgX7m+Sf2+HM888wyZmZmMGjXKeiFWajo/Pz9Wr14NwN/+9rfWuWt2A+ALXGt+\n/BfgNqCVPlA016FDh8jNzaVbt26Eh4c7uzutQgZ6B7EE+t27d2NoIFDbs7esjJuOZ/DyufM8lFX/\n6oPERPDygtrLf7t3h8jIxgf6TUXFdNZpifLyAmB6QAAK8FVhUZP73xSJiYl88MEHuLm58e9//xs3\nN7cWPZ+rmzhxIvPnz8dgMPDggw9SVVXVciczAV8CNwCWjNFxQCmwzwHHFwJaeAqq+rRNR1nGKwO9\ng4SGhtKvXz9KS0svq47sn8+ex1+rZUFIZ74qLOJgPaPqn36CCRPAw8P2ueuug+3bG3cPyq6yMib6\n+KAx/7J3cdMx0svLOtJvCWVlZcyaNQuA559/nsjIyBY7V0fyyiuvEB4eTmpqKq+++mrLnWg/cBa4\nudo2y4CjMTeH5+RAXaU38/OhZ09YsaLu/RcDrzWmo3XraNM2IAO9Q1mKVf/0009N2i/XYODbwiIe\n6dyJp0K6oAHW1FEY+uJFdQ4+Ntb+seLi1AFRQ6l3Co1GjlVUMqpWDu4pfr7sKi2luIVqlb744osc\nP36cYcOGsXDhwhY5R0fk6+vLypUrAfjLX/7Cb7/91jIn2mr+HlttWw8gjMYF+kcfhago+88FBqoj\n+rqmP03A+0AzZqeEEGzdqr6IjpT6WgZ6B5o8eTIACQkJTdrv68IijMDtgQF0dXNjrI83m4vsj6ot\nfwN1BXrzTbrs3Fn/OVPLygEY6e1VY/s1/r4YgK0ljr8t/ddff2XZsmVoNBpWrlyJu7u7w8/RkV19\n9dXMmzePqqoqHnjggZaZwukDPAj0rLZNQR3VN5TXzxLE67pBSVFg/Pi6A30ycB64yfzYiPrG04SZ\nnpMnT3IB+X0sAAAgAElEQVT69Gk6derEkCFDGr9jOycDvQNZRvTbt29v0lK3LwuL6OnmxkjzXHms\nry9JZWWU2BlVJyWBVgt1lbcMDFRvNGwo0CeXq4E+2qtmab4JPj54KgoJDp6+qays5OGHH8ZkMvHk\nk09ak8FJjvXaa6/Ru3dvkpOTWbp0qeNP8DtglZ3ttwLXUf+Sy6wsOH/e9uJSdRMmQGYmnD5t+5zl\nd/pq8/fvUT9ZbGuo05dY5uevvPJKl09kVl2LvFJFUV41f59dbdt0RVHiFEV5uiXO2RaEhoYyZMgQ\nysrKGp222CQEW0tKudbf13phaIKPN0Ygxc48fVKSesHVy8vmKavx42HXLqjvzvj95Xq66nSE1roQ\n6qnRMNbH2+Ej+ldeeYW0tDT69evHSy+95NBjS5f4+fmxwjzHvWTJEtLS0lrnxPcD/0Qd3dfFUrOh\nvjf58ePV7/ZGKr8C3YFu5sexgBZ1TX8jbd6srgO1fPruKFrqLe0RRVHSgeMAiqJEo969FQ8UKIpS\nxyRd+9fU6ZsDej35RiMTfXyt26LN0ykp5ukVCyHUnDa1/04O6fWEHzjEjIxMTEIwfrx6XevIkbrP\ne6SigkGedq7mApN8fUgpL6fA4Jh5+rS0NOsa7xUrVrh0bc62YMqUKTzyyCNUVVXx0EMPXdYqsBbx\n66/g7g7Dh6sXdXsD91HzU0BUFHh7q6N6m/2B6r/7vubHW22b2iOEYMuWLQBcc801l/EC2q+WCvSz\nhBD9hRCWaDcTsFxdPAFMaaHzOp0l0Df2gqxl5DzR99JdoV11OkJ0OlLKawb6jAzIy7MN9C/knONk\nZRXrCwr5tqjYOiiqb0n/UX0FA+0t2wEm+foigO0OSB9rNBp56KGHqKqqYu7cuUyaNKnZx5QatnTp\nUnr06EFSUhJvvfWWs7uj0mrV1QIeHrAQyAI+QQ36Fu7u6i/5ggW2+38N/LXWtklAElDW8OnT0tI4\nd+4c3bp1Y/DgwZf3Gtqplgr0EbWmaQKBvGrPu2x6wkmTJqEoCjt37qS8VqC2Z1tJKT3d3OjtfmkK\nRVEUor08bQK9Odke1etxFBqNbCwo5PEunemi0/JJXj4DBkBQENRVCyXXYCDXaGRAHSP6sT7euCsK\nW0uaP0//1ltvkZSURI8ePXjttWaui5Mazd/fn+XLlwPwwgsvkJ6e7uQeAa+8At99BxeAH4CnUYP8\niFrt6hiAMBAYWmvbJKAK2NXw6S2j+SlTpnSY9fMWupY4qBDidQBFUaYoihLX2P2WLFli/Tk2NpbY\nupaWtGGdOnUiKiqKlJQUduzYwZQpdX94EUKwraSUKX6+Nr94UV5evHnhIgYh0JmfS0pS/waqLz2P\nLy7BCEwPDKDYaOSLwiIEgjFjlDqXWB7Rq+uY6xrRe2k0XOHtTWIz5+mPHTvGCy+8AMD7778vywK2\nshtuuIH77ruPjz/+mIcffpjExMTLvwB5HHgFeBbo38yO/Yg6XTMTGNbMY01AXdNfx3tDdZb5+fr+\nJtuDxMREEhMTm7ZTQ1nP7H0Bs4FZ5i/Lz5OrPTfN/PPT5udeqfb8dGCBnWO2XHq3VvbMM88IQDz1\n1FP1tjtarhck7xPvX7ho89zqi7mC5H3iaLneui02VogxY2q2W5B9Wnik7BeVJpP46GKeIHmfSC0t\nEy++qKYtLi62Pa+9Y9f2/OkcoUneJwovM5ul0WgUsbGxAhB33333ZR1Dar7c3FwRGhoqAPHuu+9e\n/oGWC/Uv9XA9bfYJIRYLIer+tVJ9LYS4TQhhvPzuNFVFRYXw9vYWgDhz5kzrnbgV0FLZK4UQHwoh\nVpi/LD9b5uOTUNMeAfQF9gBrgQjztohqz7ukqVOnAvDdd9/V226neQ78SjtZGweap1WOmO8iNJlg\n717b+fkD+goGe3rgpihM8FUvciaVlTFmjLqPvZt0j1RU4KYohHvUvY491s8HE7DjMkf1K1asIDEx\nkS5durSdOeIOqFOnTvzrX/8C4NlnnyXT3kXOxvgZCAHqq6G9H3gJ8xKMetyEmi+nFVc37tq1i7Ky\nMoYOHUpYWFjrnbiNcPg/tRAiFZipKMp04JgQItW8DfM0Tr7lsasaP348/v7+HDp0iJMnT9bZbmdp\nGQFaDYPtzJVbplUs0yxHj0Jxcc35eYC0cj1DPdW18OHu7vhrNKSW6633pNibpz+ir6Cvu7t1Ssie\ncT4+uCnKZS2zzM7OZoH5Ytrbb79N586dm3wMyXGmTZvG7bffTklJCY888sjlZSf9GbiK+pdPDjJ/\nr2e1V6Pl50NKigMOpHKVaZvL1SLvqeZR/ufCPFdv3rZCCBEvhKgnkYVrcHNzsy7f+v777+tst7O0\njHHel3LNVNdJp6OLTsth84jeMjKvfSH2VFUVkeZAr1EURnh5kVpeTpcuEB5eR6CvqKjzQqyFt0bD\naG+vJs/TCyGYO3cuxcXF/O53v2PGjBlN2l9qGe+88w6dOnVi06ZN/Pvf/27azqeAk6iBvj7m0X7u\njlxWrFjB+vXrKSsrU3PbfPQRnDljfz8DkFFr29NPwzXXqGuKS4EuQBO7Xd23334LwPXXX3/5B2nH\nOs6tYa3MMn1TV6AvMBg5oNdbp1vsGejhYR3RJyeDpycMGnTp+QPmG6oiq93dOsLLk33lekxCcMUV\ntjlvjEJwrKKyzgux1U3y9WFPHXfo1uXTTz/l22+/JSAggH/9618dbnVDWxUaGso//vEPAJ588knO\n1BV07fnZ/L2hHGD+YAw1kp2QTUlJCUeOHOGzzz5DpKaqNWF31bE0Zj7qevjqHzRGjYLcXPVu2sPA\nRcCv8V2uLjs7m9TUVHx8fNrlAg9HkIG+hVhGDvHx8ej1tne47i4rQwDj66mqNNDTwzpHn5wsGD4c\ndNXWSR0wH9cydQMQ5e1FicnEicpKxoyBU6fg4qGfwaimZMisrKRSCAbqGv74HuvrixH1k0djXLhw\ngcceewyAN954g27dujWwh9Sa7r77bm688UYKCgr4/e9/3/gpnOuBz4Hh9TcTQnAu6Bxdcrswe/Zs\nrrvuOjIyMjhmucu1rjtihwC5wLlq28w5PlK3bGHT25sAqOpbT+6en4HHsJuC4ZtvvgHUm6Q8PT1t\nG3QAMtC3kG7duhEVFUVZWZk1LWp1O0tL0QBjvOvOZTDIw5PzBgN5GZ+z5Ko4Ro2smdMgTa/HR6Op\nsQZ/mPkX+UC5niuugJ7BWQSnTIJDat4TyyeEAft+D6KeHAnAeB9vtEBiI9fTP/bYY+Tm5hIXF8dD\nDz3UqH2k1qMoinWZ65dffsm6desat2MnYBpquoF6ZGVlsSNyB/r79Pj4+DBy5Eh8fHzYm5MDoaFq\nCuI1qEszqrMssayecHP4cE717s2Xp04ReD4QoQi+O1LP4obDwNvAMdunLIH+pptusn2yg5CBvgVZ\npm+++uorm+d2lJQxwssTX23dfz2WlTe7LyhMHPgTt43+osbzaeUVDPX0qDHHP8S8T5peT3Q0PDjp\nI7WGbJ+7gUureAbm/wzn6k/T4KvVEuPduLw3X331FWvWrMHb25vly5fLKZs2qkePHrz+unrpbP78\n+Vy4cMFhx965cycnx5yk54tqakutVktUVBRHvbwoGTdOzU75JvBhrR3tBXovL7bdcAM+VVWM8huF\nvque1MOp5OTk2D/5lebv22tuLisrIz4+HoAbb7yxGa+ufZOBvgXdeuutAGzYsAFTtQxjVULwS1lZ\nvdM2AIPM8+jf5k0k82IvxgTXvI6dptdbL8Ra+Gq19HF3I61cj5cXzJzwJQfPjwPfPgAcragkUKuh\ni06B9PcafA2TfH1IKiun1Fj36L+goIC5c+cCajm7iIiIOttKzjdr1iwmT57MxYsXefzxxx1yzOLi\nYtLT0xk5cmSNimGRffogNBrSo6PVfPIHgNq1ZjoDXakR6HNzcznWtStjKivRntDiHumOp6en3U/H\ngLriJxibQJ+QkIBer2f06NF07dq1uS+z3ZKBvgXFxMTQq1cvzpw5UyObZVJpGaUmE7G+dQT689sg\nN4lwD3fcFIWUvCrW7L4L/7LNoD8PwIUqA+cNBoZ62c45Rnp6kqavgPIchnTdy4bdN1ozWR7RVzDA\nwwMl/AHI/gr0F+t9DbG+PlQJwa568t4sWLCAnJwcxo0bx/z58+v/R5GcTlEUPvzwQ7y9vfn000/t\nfuJsqrS0NIQQjBhRM59BqL8/fiYT6SEh6sqaMuzfDTsRqHZbx8GDBwGIeuEF2Ana/2mJjo7myJEj\nlNr7XVRQ75Kt9T6wYcMGAG6++WabXToSGehbkKIoTJs2DYDPP//cuj2+pAQFuNrP1/6OKc/Ar3PQ\nKQr9PNw5bqog+eKdKMIIp78GLl2IrT2iB3UVzmG9nsoz6trhDb9MtWayPFyhZ5CnpzqVIwyQvaHe\n1zDB1wcN8FMd0zebN2+2FhFZuXIl2nqmoqS2IyIigr/97W8AzJ07lwJ7Fc2MQMPpmgDYv38/3bp1\ns7lnQgkJoX9MDCcKCjDtM4827AX6tcAHlx4eOnSIHj16qGkz3IDOEB0djclkYt++OorTXgmkY72o\nW1VVxcaNGwG4/fbbG/dCXJQM9C2seqC3rHLYUlxCtJcXwTo7qYZKMiD3F+h9J6BO3+T6VeDTYzh4\nhUHhIUCdtoGaSystIj09MQDp7r0p6PQg+7OG88svUGQ0cqbKwGAPDwiKgpBJYKq/CpG/VstVvj5s\nLCy07WpJCbNnqyUHFi9e3OEyArZ38+fPZ/z48eTk5PDUU0/ZNtiFeiG2gTTA58+f5+zZswwfbn9Z\nTp8+faioqKD0F/NgoYHCTqWlpeTk5DBgQM3bcLt06UKPHj3qDvTTUVcHmcdP8fHx5OfnM3To0A7/\nuykDfQsbP348oaGhZGRkkJqaSonRyK7SMuLqGs1nmpck9FJvNOpe5YExrIIRo4CbjsBI9UJaWrme\nIK2WMDtvFpZRfprvCPyvXYWfn4Zff4XD5hU3gzw91AtjUxJhwKMNvobpgQEc1FdwqNYy0UWLFpGZ\nmUlUVBRPP+2y9WRcllarZeXKlXh4eLBq1So2bdpUs4FlUU4d1cws9u/fj6Iolwq9rwX+cOn53r17\nA5AVngUvYw3EdcnKygLUN4jahg0bxvnz5zl//rztjhGoq4PMM6Lr168H4I477qj/hB2ADPQtTKvV\ncttttwHqqP6HomKqhOAG/zru/shcA50vXTzVZXuCG4SOrAS3S/sc0OsZ6ulhd3XLQE8PtKhvBhqN\nunz5l18uBXp7KRfqMy0gQO1/waVR/Y4dO3jnnXfQarWsWrWqxgU4qf0YNGiQNWvsI488QnFxsfqE\nEVgP3Ei9NyoJIThw4AB9+/bFx7K44DfgPay1XP39/QkICOBApwOwqOE+ZWVlodPp7N6HMWTIEBRF\nabD4uZy2qUkG+lZg+UX79NNP+bygkM46dTrERuEhKNhnnbYBKElTg7Im4tJoWghBmr7C7rQNqOUA\n+3t4WKd3xoyB/fvht1I9OiCiEXfFVtfd3Y2rfHz4OC8fkxDo9XoefvhhhBA8++yzREdHN+l4Utuy\nYMECRo4cSWZmJs8995y6MRE4C9xZz47A6dOnKSgoYOjQaoniB6K+UVRLbta7d2+ysrIadZNWVlYW\n3bt3V6/3HMmBt98B8zJQX19fIiIirBd/65KQkEBeXh6DBw+u2bcOSgb6VhAbG0u3bt04ceoUX+UX\ncGtAgP2EYh5dIPp16HXpo+apn9WgfJIK67YzVQYKjEa7F2ItIr08rYH+iivAYIBfz1fQz0PNdNlU\nj3TuRHpFJT+VlLB48WKOHDnCoEGDrPnmpfZLp9OxatUqdDod7777rrqE8V+o8/MNLD1PS0tDq9Uy\nqHpuDsuPb/8I5vrAvXr1orS0lLy8PNuDWGRB1cYqcnJy6NWrl7ptchA85lsjfUJkZCQFBQWctldA\n3OyTTz4BYObMmfW/gA5CBvpWoNVqufvuu+HqyZQpCncEBthv6NkZBj+lXnRFzee0d6sWrxIdh/SX\nAv0+c+WpYXWM6EGdpz9eUUmZycTIkQZuvPFbDpSfwePMaQ4fPtzk13B7YACdtFr+fPgoS5cuRaPR\nsGrVqg57S7mrGTFihHU0P+uhWRh8DGoOmnqK0JtMJg4cOED//v0v/R5UFYNivj/j+xPwww8A1sBt\nmX+3axXopuvQVmrV9oXAGU9Q0tUc3WaDBg1Cq9XWOX1TVFTET5//hBtu3HfffY16/a5OBvpWct99\n98H0O9Dk5HCVe8Pz2UIIvvlmJ3fd9U+6l1/gl3zz8jdjJckX1EAd7VXHX6GpisjMDxDAwXI927at\nZdi4VHJ9felVVsratWtJSkq61D77S9h+p/rOUgdPjYZZAX5s07khBgzk2WefZdy4cY19+VI78Kc/\n/YmhQ4dy9NhRFnVdpOaWr0dWVhYlJSWXLsIaSiHxBjj6KIRUwilf9eMk0LlzZzw9PTl16lTdBxwG\nikkh5GIIPXv2hBTzJH/cyUt1NAFPT08GDBjAgQMHatyIaBG/LJ7j+uMsGrSI8PDwpvwTuCwZ6FtJ\nSXg4jB6Nad0avjOnTK2LEIJvv/2W5OTNFBQE0qeshCNGI78dPQplp0g+vYMBSjl+da1Zz0shMnsV\nAF8eOMCxY8f4rUpNsjZn4lUMGDCA77//nowMc25Y/QXIWgt5e+wfz+zCm29AYSG+f3qeF198sWn/\nAFKbZ1l9o9FoeOONN2oOBuzYv38/bm5u9O9vri346xy4uAuuXAuPZYBxKUxUU14qTypcfeTq+gO9\n+f2iv74/Hu7u8MVqdcMgPezeDdWCemRkJKWlpXbrPfwz/p+c5jTzy+fbTXLWEclA3wpMQvDsmbP4\nVlbCurWsXLmy3vZ79uxh7969lJWNZ8OGe5k1fiwmjYYPEhMp14aR7B1JdFXtBN7VXNxB38psPBTY\nlnOW8PBwxGA1/0iE3pfp06cTHBzMhg0b1ALmvW4HjTuc/G+dh/zxxx9Z/c9/oln+ASWRw/iomfVk\npbZpzJgxPPnkk5hMJh588EEqKirstqusrOTAgQMMHToUd3d3OLVR/f2JXKxeY9JsAH5TA70RWAE9\nCnpw8eJFNUe9HcZwI1W6KnoV94LjK2CfGwSWQ9SNkJeH9a4/oH///ri7u9tM35w4cYLEnxN5ze01\nOmd2tk2g1kHJQN8K3rpwkW0lpfylawgeJhM//PDDpdE0QGWhNZNkfn4+mzdvpl+/fiQkTGHMGIUr\nA9Wi2se9vPlixw4ydaGMLtpR91TLhZ3ofHoxRAiO+fkxfvx4SsPL4YwbJ/bqcHd3Z/r06ZSWlqpr\np90DocfvIOMTMNjeCpmfn2/NRvmX4ZFc4+fL49ln2FRU7Nh/KKlN+POf/0z//v05cOAAixcvttsm\nLS2NyspKRo4cCRW5kDQXgqJhqGXVTqJaxb5zZ9h6HErA4yp1cJCdnW33mDnnc7jQ5QIhZ7vAwVfB\n2Buu9oQpcbB4MfhdWufp5ubG4MGDOXToEAaDwbr9gw/U22sr7qiAUajr+e2frkORgb4FCSF4/2Iu\nC07ncGuAP4/36snMmTMRQlh/IQHY+xh8H4Uwmfj2229RFIXJk28iNVVh7Fh1eWNfd3cuRPRjwxk1\ne9+Ugk1QamdULwRc2I7oPJ6w7FPkBAbh26sXqR7FaPb4Yike37VrVyZMmEBqairHjx+H/o9CZR5k\nfmrzGh599FHOnDnD+PHjefbpp/m0Ty8Genhw84mTvHn+AlWXU5pOajsygVuxJgTz8vJi9erVaDQa\n/v73v1uzP1aXnJxsvVMVk0G992Psv0Fjvv701lvw7rvqz1s7ARA4+CM0Gk2dF2SzsrI4OuAoHmOO\nQckJ+F8hfKZA796wZAn06KE2LD4O5xKJjIykoqKC9PR0APR6vfXT8tz5c+ET1LX8v2/2v1C7JwO9\ngwghMAlBgcHIvrJyPriYy4T048w7dZrr/fz4X69eKEJh3rx5AKxcuVL9WFxZCMe/gbKbOPplOiXb\nS7il6y1krgmgrwEmTVKPf0uAP2ne3vzcJ5wuxkqGVaTD2W1qkqg84Iz5q/gY6M9yQTOYnseOIRSF\nuzJPUWgyMfSiP9u/Q/3D/g0mhUwiODiYb775hsrAsRA4DE5/A0Wot7+nwsZXN7Lj0x308erDx//6\nGK1WS7BOR0L/CK7z8+Wp0zn0PXCYBafP8EVBIQfL9VyoMmCQwb/1CNTMkCbUaZLqX/YUg9glEO8J\nxDSB6CsQmwXGdIFBCKqEYMz48SxavBjh4cE98+Zx+Nw5TldWcbKikq0nT5JcVIRvTAy/lpWz3ehL\nQvQn/OzWl5Syco7qKzjTrx+FEyZgMAlYFwSR6WiL3qJXWAAnsrMpMhq5UGUgu7KS4xUVZFZWknYq\nm5Q7DqB7XkC/udD9dxgVQb7BQFZlJQfL9fxaWkbC3lf4MXkpWUHeXAgL46v0dA6W63lnwxfkAsPH\nj2fkmDGIQQK+Af7SKv8LbZpyWYWCW4CiKKK1+nLjZyf46BF1uZciLn3tH6TnluUnrX83AoEQEHXQ\nk80P91XbV9tn7+ByrlxtW/I+5oAXSff2tz1xDIhfBSNHjiQ1NZVPPvmEe2IuwLo1sNi2uOseBQYX\ng48PpOsrGHToCCbgie/O8+aSWBRDrfQHY4Gfy+H8Nj5PSOfkuQo+u+kWDEmQMC8Cn1INGlFtDf2V\nkPXfLFavXs2oUaO4afJI8OoK2zX2y8ZNoEYaWCEE3xUVE/9jCS/MClX/xcyHFwrsHl7G7f/MRKco\naBXQoqBRYPQ+Lz5+rJfN4fcOK+fed2xHe6NTvfjPfMv/16X+7xlRxl3/sm0/JsWL/z3a29z+0vak\nEWXc+b5t+yuSvVkzt7fN9qToMmZ8kGm3/fpHbI//S3QZt6/ItF7/E0KNwWOTvdk4u49N+13Rpdy0\n4qTa1rIPMGGvNz88bJvqeXt0KdesOlGjLeb2CbP62rTfFl1C7MoTNdoCTP7Fl/h56vFzOlfxnxsK\nePvOi5wKqz/v0eVwr1QonDiU+YtP8cl1F6jUuDe4jwbw0WgQQImdVTWNpQG8NBp0CijmquaW3x4F\nSL2+P53ydRjcBAbtpd/dSRuOcTHY9rxbb4ugc55typFJG45zMfjSu+q6Pr2YVFeKkxagKApCiHpv\njrGTVcv1je7hzcGb9AjF/B+vqIGpoLuRB4M7ocEc0M3PB/bXsvd+tfSfoliehKIwEy90DeHSJgUf\njYaB7u5cfM5AsE6rpiiwNOim/qfMnz+fWbNm8ebrf+ful4pRhg5l7x/3kpGdwZSpUwgMCeSxP4LR\nA94130Db39ODLyP6sD3nLMF5aZy9eyhh4WHgibrW2VM9PjovzmmjSDu+m8mTJ3NHeG/W5xZSeb+g\ntEJh6QqY9iBMuEFt36tXL8aNG8euXbsYNGgQ/fppYChUflnJs48/y9mTZ7l67NXMvnc2SkjN3yVF\nUbgxwJ8bR/pT9YAgr8pAvtFIhUlQYRKI3oI/dOmMQQiMCIxCfQPt0kvLgZttyysW9zFye2CgzfYu\nfbSkTqvW3tyNvN7223cO15I0s9ym/cVe9tsH99Wy837bi8u5Petuv/Xhau0tx+9hsLavHlA699eS\nMLdalS5ztM/vYeLB4E412gJ0Gqxl0xPFNbap7Y081qVzzeMrEDhMy5YFttdLCrsb+VNoCJb74yz7\neE1S+O8neeT2NZHf24iiwCw6WX/fq7ctKSjkrdeXoi8u5q4ZMxg1YgQ7EhMZPngwUUOH4qYouCsK\nbopClRCUmEyUmr9KTCb05Sa+X1KM3z1u/PHUNnTlZ9l9LopxMaPpGhSEh6LgrmjILSxg0/btDBg+\nAr+QEEpMJhQgQKshQKvFX6vFT6PBV6PBV6tBl/4+FZnruNBvMRt2nKZzWBjvrliBd1AQL7z8Miad\nG+UmE+XChEGogzagxptw6j3leBVp0FYqaKp9+rk2xA+9n+2gM/PaSs6VGGy2X9vFF73/pfYhbm0w\nrAoh2sSX2pWOQa/Xi7CwMPHIZIT4L+JM0gdiyZIlYtOmTUIIIfLzhdBohHjhBdt9TSaTWL58uXj7\n7beFyWSye/wvvvhCvPzyy6KsrKzGdqNRiNBQIe68s2b7qqoq8e6774o33nhDFBcXCyGEmD9/vgBE\nv379RFFRUfNftNRurV27VgDC09NT/OMf/xCvv/660JeXN/1Ap74UpvXB4t2//V7s2LGjxlO7d+8W\nS5YsEQUFBY07VlWZEF/2FeLLvmLNf1eL559/Xri7u4tFixY1vV/tnDl21htf5Ry9E3h4ePDHP/6R\nkxfgh6Nh/DexkC5dulgr1P/wg7pk+IYbbPdVFIWxY8eSm5trvQhVXXFxMb/99hvR0dF41bqhSqOB\nqVPhu++geiJKnU7HtGnT0Ov1rF27ls8++4x33nkHNzc31qxZg59fPVmtJJc3Y8YMHn74YaKjo8nP\nz+e6aC88tl8H5edsG+fm1ljvXkP3m1BuPQUBQ9UFANUcO3aMoKAgAgLquGsc4NNP4aqr1OPrvOCK\nD6H8NP2Cy9BqtVxzzTU88cQTzXilrksGeieZM2cOyTmd+fz0jZSX67n99tutGSC//BJCQtRkZPYM\nGTIEf39/du/ebfPczp07EUIwduxYu/vOmAFFRVA7I23Xrl259dZbyc7O5uuvv0an0/H+688wakQD\nycOlDmHBggVcc801GHK2EZ71OKL8LGjtJMe75x6YPNn+QRQN6Lzo378/J0+etK7Rr6ys5ELWPob3\ntZ0mq8Fkgu3bwXIjV+jVcNMRlq7eRlJSEjExMXWu0e/oZKB3Ep1Ox6OPPkq3bt1ISUmhS5cuAFRW\nwvffw803Q103vmq1WsaMGUNGRkaNYskluRns2bOHYcOGERQUZHffuDgICoK1dm4k6d69O7t37yY8\nPJznn3mUB8Leg8Tr1XXSUod19OhRNmz4nCvD0nnuip84nWvgnQO3qPdfVHf+PGzeDOPH13u8AQMG\nYLptgtsAABH/SURBVDKZOGK+AerYsWNMDNzCxJJHoaqk7h2nTgU3N1i3zropLaOI9evXs23bNgIC\nAli7di3nztn5pNHByUDfmoTAYDCQkpLCe++9h06n48cff+SLL77gu+++A9SRdmEh/O539R9q1KhR\nBHobKd1yOyL7Gyg8hPeP/RjotY+JE+0tmVG5ucHtt8PGjVC9epzJZOLee+/lhx9+4JdffsHDrytf\nnr4G0/ldmL4eAunvq7lMpA7j4sWLfPnll6z59D/c1+M/XOP/KeU+w5j0V3h80Rv8YE5YZrVmDRiN\ncNdd9R63V69eBAUFscecv+Zk8gai/Peh9H0Q3OpZrRIUBNdfrwZ68/TQwoULMZlM3Hfffdx///3m\n+ggr2LNnj908OB1Vh1xemZOVjrbkMCAQwmT9LjQ+VPqoJccsfRFCoBiKcS/Zp65EEJf2MWl80ftF\n12gLoFQV4FWyF4QBbdVFlLJTeBSnkl8ZwPpT16PX6wkLC+Pmm29mzZo1PPnkk0RGRpKamsodd2jZ\nsQOys9WgXJ/kvb8SnnYdOndfjG7BeOsPsjv0f0y8pv5CC8nJMGoU/OMf8Nhj6rZFixbxyiuvEBQU\nRFJSEp07dyY+Pp4LR7cwtcvX9PTKxoA7GYHzKQz5P7y8vHB3d0ej0aDRaHDXZ+BmzKu2LEkDioLB\nszcm9842fdCWZ6Ktsk1Za/DsifDoYqf9SbSVlk8Wl35PDF69MbnbaV92Am3VRZuc5QavcEzuXWy2\n68pOoK26YLO9yisco5vt8XVlx9BWXrDZXuUVgcFu+3S0lbZVkSo9IzC6h9ppfxRdrfZCCCq9+mJ0\n72qz3b08HV1lzZGsQFDp2R9DreMLIfAoP1qjvUBQVVFGrrEHOUVunD59mgsXLlg/PU7p9D0a/37Q\nfx5//stfWbx4MX5+fmzfvl0tIWgywZAh6t2rDeTIAdi9ezc/b9nA1YMhongFPl46PG49BB7B9e/4\n3/+q00PbtvGTwcDkyZPx9fXl+PHjhISEUFRUxLGND9DDuJ0zVX3RBA7EIzAcd+9OmAJGoPHtiUaj\njm8VRUFRFHT6k2gNBVxaHqcyevXC5NbJpgva8kw0Vfk22y3tg4ODWzWra2OWV3bIQL966Twe7P6+\nzfbT+m6sOPWIzfZuHtnM7rXistubhMK5ilCyTFGc7TKXYcOGER4ejqIoVFRUMHDgQDIzM3nttdX8\n6U8P8Pjj8PrrDb8OIQTJ373KiPwX0GmMJCmzGDnj/UYV6B47FvLz4eBBWLlyOXPmzEGr1fLdd99x\n7bXXWtsVFBRw6OBBijN+oHPZVvbl9SVLb7vm/JbQjUT7p9ps//LcLaQW2daiuyVkI9EBsn1ba//V\nuVs4YriSbt26ERERQWRkpM3FeCEEd999N59++ik9evTgl19+odu+ferUyn/+A3ffbXPc2gwGAzmf\nRtFTewC98EE35Qd0oVc2uB+lpdCjB6Y//IGRX33Fvn37+Mtf/sLzzz9/qX/HV1N+6EN0Rb/hzqWp\noA1nb+O34hE2h7wt9HOG+9umPP7i7G3st9P+1tANjPDfX2f7e+65h759be9raCky0NfhxOFkPIr3\nAgqKYl41ryiYdH5U+kVb+mP9rhhL8Cg9aG0vLIvvdb5U+Qyu0RZAYypHV3YMFA3/3965BkdRZXH8\nf/OYPCEhTAATFhJEgglCwvPDrkUBWRZYpMSwuIhlWbtY7hdLWRG1YAlqFQRkRao2ggValuKHRaDc\nFaEiK7hVsmtEgWgQw0J4E5IJJCaYZDKZ+e+H7klm0h0mgRky6ZxfVVe6b5++99z0ndOnb997rsdm\nR2xyBhIHJHVpgPfs2YNFixbBZtuA1taV+OEHoCdrGV+/8DXczgak3jfLdGlBMz76SPsw+/zzJ/Dm\nm5Phdruxfft2LFu27JbXtbW1obm5Gc3NzWhtbYXH44HH40Hkz6cR4awGPW7dJ9KmnbXGjUFbJw8U\nAGxNPxo8UABojc/qQr4CUS5fD1crxRk32tTDjWk5Y+pBu+KN8gBgaz6DKFetify9Bo9bKYXo5rMG\nj14pBVfcKLhtQwzpUU1nEemTv/c+ueIy4fGR96ZHNVUiwkTe+0bimwZ432D0NySfJuCOyzCXbz6P\niNbrfmm2mDjEpY5DVIK//ma0tLQgPz8fR44cQV5eHr5ctw7xGzYAJSWALfCkKABwt9TBcboEg0fn\nIzre+NbXJdXV+OvOnVixYgUyMjJQXl7esYxhZ5w34Gy4hOZGB5ojhsKJAYahh9E//4AoZzU6h7p0\nxo81bYsxTadM265XPj09vWt9QkB3DH2vj5/3buhH4+g74/F4OGfOEgL1TEv7z10p0+0mx4xppFLn\nCMRx9erVd6VcwTo4HA6OHj2aADhz5kzDvI1Qce7cOcbHxxMA9+/ff1fKDGcg4+j7BkopDBv2NoAB\nuHr1T9i5c2fIy/zmm69x6dJikBm4//6P8eqrr4a8TMFa2O12HDhwAMOGDcOhQ4ewcOFCtLQYZzwH\nk7a2Njz22GNoamrC4sWLMddssolgJNCT4G5t6McefUkJqRQ5c2Y5ATAhIYE//vhjyMo7evQok5KS\nCID33XeAAPnBByErTrA4J0+eZGpqKgFw/vz5bL6dWbPdZNWqVQTA9PR01tbWhqycvgS64dH3uoFv\nV6SfGvrDh8mBA8nx48mGBg+XLFlCAMzOzuaNGzeCXl5JSQkTExMJgAUFBWxoaOWMGWRkJPnWW2QX\nURUEwR+Hg3Q62w/LysqYkpJCAJw+fTrr6uqCXuSuXbsIgEopfvHFFx0nPvuM7MdhOsTQhzFXr5Ir\nV5JRUWR2Nnnpkpb+008/MTs7u/0HE0zv6L333mNUVBQBcOnSpXTqP9SGBnLePK01zJhB7ttHulxB\nK1awEleukEVFpN1Ovvii36mysjKmpaURAHNycnjhwoWgFXvkyBHGxsYSADdu3NhxoqqKjI0l772X\n/Ogjv4dPf6E7hj4oo26UUnkkj/scFwCoBzCR5OtdpXXKg8HQpTvsebsWQ48d0B4wHn0MvIdoSbTj\n0oT50ExexxbzUw0yyv7RkQDtAdmckIozEwp8k0ECcQ3VGHNil/dJqn3MJ+GIGIoDyUtQXg6cOKEN\nPX7ySWDzyiok79vZnkldXR22FhfjdGMjqmbPxt69e/2/4ldWAq+9Brhc/tvo0cDmzYb6Ok+dwrcL\nF6K8ogJtAMZPnIhfTp8OlZkJPPMMAE2XbduAV14BWFODJba9SBsZjZQhUUgaFIG4OKAtaTCq8uYh\nJkaLm+MlttGB4d8f0A58IiU6B9hxZcI8dB4IFNPgQPp3B9D5djsH2HF5/Dy/NFKTH/7dfkO9WhLt\nuDzht4Z8Yhoc+MV3/uvykmi/v95jX/07ywNAc4IdF8fP98vDKz/y+31amq8+CXacf+AhQ7mxjQ5k\nln9iyKc5MRXnH3jIoH9sQw0yT+4z5NOcmIpz4x4y6B/XWINRJz8x5NOUkIrKcQsM5cY11mD0D/80\n5NOUkIozOf4z9UjAfq0cv/rsL0hxVMBecwoAcO6+2Tg8dyMcaRP88qivr8e7774Lh+M64uMHYsGC\nRRgx4l60tXU001vtk9qgnejoju3GjWs4eHA/XK4mjBs3FnPnzoLNptrPj7r4BWbvfRopjtNwxiej\nNnMKLk5ZhDMzjUOfE2sqkfq/I4b0m/ZMOLKMwzt7Kg8A+fnAPfeYngoJd2XUDYBZAM74HOcBeETf\nX6Yf+6Y9BSDXJJ8QPe+MPBhT2tmWkwBLMcUsmVNw5/KHAX6NKRwxgpw5kywsJCsqdIVKzfM/rnvf\n06ZN4yWvy0+Sx4+Tw4eTmZnkmDFkTg6Zm0s+8YShrseOHePvs7J4AeBVgE2JiWRSEpmQQD74oEHe\n6ST/veG/pvp8hamm9R2L4h7JT8VXfUz+cJjpE1z5wwHk78dJfo8cfowFXIkiZqPcVM58a2NkZCsT\nEz1MTiZTU8m0NHLECM0JHzuWHDeOzMsjp04lp0zRmnJOjta0U1MbqdQlAlW02Ro4YICHsbFadFff\nciLh4lx8yu34I48hl1vwjM/5jvv3ON43VfR9PG6qf0/lAa0n6W6i207caguWR19C8jf6fhGAz0ge\nUkrNAjARwOBOaXkkN3XKg8HQpTucO9WCiKorWvztCKVtCkBMDDzD0nR9fDZXKyJvONrloXR5mw0Y\nPNhPFgBUmwsRNxu0fT3/dW8U4bXC1YhIMokE6XZ3hJNsz0ThdGUlfj1vHi5evIiUlBQUFxfj0Ucf\n7dZY+WvXrmH9+vUoLi6G2+1GRkYGdu/ejUmTJgX+B7W2alEIdXfL4yZcLqCFMWixD4fTqTVpQPv7\n5qbVWLH4D35pAEBbDNqGprdnS+r/J2cLIquv+lXXK+8elm54A4hobUFkTZVfmlK6/NCO+9WO04ko\nh788AMBma7+/ftc4nYisqTKUy2gbeE8aNm1aixUr1nZc43QiouaaMR+bDZ6h9xjyUa1ORDj8x10r\npec/dJip/pG1JvFabB3ynfVXjhpjPjYbMLTzHABdvrZjDsC6v72BVc/8Wcs/dYhRf5PmFljGg9df\nX4+1a9fA4/Fg5MiR2LBhAxYtWtStCX21tbV46aWX2pcGXLZsGbZt858M6PEYX2pdLq35+kY/2LJl\nLZ59dq2m581GRF43zq9gQiLcduMM5Z7KA5o3Hx8fsIpB466No4dmxL3726B77NC8/SIAWzulrTfJ\nI1QPvLCgsLDwtq6rrq7mnDlzCGgdQNOmTeOHH35oOmbZ5XLxyy+/5NNPP824uDgCYEREBJ977rn2\nOPOh4Hbr1leQ+t0+paWlzM3NbW+/WVlZ3Lx5My9fvmwqf/bsWa5atYqDBg0iANpsNm7durXLtRe6\ng9XvH7rh0YfhUiiCL0OGDMGnn36KHTt2YM2aNSgtLcXSpUsRHR2NCRMmID09HZGRkaiurkZZWRlu\n3uyY8v3www+jsLAQubm5vVgDoT8zdepUHD16FO+88w6KiopQUVGB5cuXY/ny5cjIyEBWVhYGDhyI\n+vp6VFRU+C0cnp+fjy1btiA7W0Jl3ykBu26UUk8B7d+clL5fSfKQj4xv1816AAepddMUAMiE1nXj\nl0aTrptgVUoQBKE/wTtdM5bk9m6U41vILgCTABwCMArAQT19sklatxUVBEEQbo87DoGge+iTlFKP\nAAD1YZb6R9c6kidInuicdqflCoIgCN0jbKJXCoIghAKl1As0mbvTn5CgZr2EUuqF3tZBCIxSqkAp\nNcuq90sp9ZS+FfW2LqFA70XI7209QoFSKk9vnwWBZMPS0PekAn0Rizc+yxgOpVQetKFrnwOoV0pZ\naviS3g4P6t/hRimluljVWwhTXia5B0BmoLYZloYePaiAED5Y0HA8Ci1sBwBUwnoP51HoqFOlfmwZ\n9NAsn8N/sIgl0J3grwGA5KZA3z3DztD3tAJ9DSs3PljPcCQD8F3YNsCCpn0LkttJete8nAjgm97U\nJwQM6m0FQsgUAIP13o+A3YphZ+jRwwr0QSzb+PqB4bAkehfVt1ZyqnSHyjvXx6ojTq77jHK8ZTd3\nuM6MvU7yuFIqXylVoHfj9AluNcGsnzQ+KxmOOgAp+n4ygOu9qEsomUXy5d5WIsiMUkp5J2sOVkrl\nWqA9+nId2lszoHUvTgbQpZ3sFUMfYLZtjyoQbgSYYNbnG193ZkrDOoajq8l/lkEp9ZR3lrpSapbe\nrdjn8TqHentN6mV1QsFuAF4vPhnA0VsJh904et0QFpDcpHfdnCW5t7f1CiZ641sJ4Hd9zdAHQjcc\n2/X9Pm84lFLLAJyDFrZjRyD5voT+8XwXtDeXQdDa46FbXyWEC3rbrAMwOZBjFXaGHuhZBYTwQQyH\nIIQnYWnoBUEQhOARjqNuBEEQhCAihl4QBMHiiKEXBEGwOGLoBUEQLI4YekEQBIsjhl4QBMHiiKEX\nBEGwOOEa60YQeh09UNQoaCE5JgNYT7Khd7UShJ4jHr0gmKCUytDjpdRDi+fzdzHyQl9FDL0gmEDy\nvL47CcC/rBaTSOhfiKEXBBP0UMuAFsyswedYEPoc0kcvCObk65FUD+pLIt4IdIEghCsS1EwQBMHi\nSNeNIAiCxRFDLwiCYHHE0AuCIFgcMfSCIAgWRwy9IAiCxRFDLwiCYHHE0AuCIFgcMfSCIAgW5/8w\nGUbXKIKrNwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# quick and dirty visualization:\n", "\n", "N_plot_min = 0 # quantum number of first eigenfunction to plot\n", "N_plot = 6 # number of eigenfunctions to plot\n", "\n", "WF_scale_factor = (np.max(V_sampled) - np.min(V_sampled))/N_plot\n", "plt.plot(x_vals, V_sampled, ls=\"-\", c=\"k\", lw=2, label=\"$V(x)$\")\n", "\n", "style_cycler = cycle([\"-\", \"--\"]) # line styles for plotting\n", "color_cyler = cycle([\"blue\", \"red\", \"gray\", \"orange\", \"darkturquoise\", \"magenta\"])\n", "\n", "for i in range(N_plot_min, N_plot_min+N_plot):\n", " # physically normalize WF (norm = 1)\n", " WF_norm = simps(np.abs(psi[:,i])**2, x=x_vals)\n", " psi[:,i] /= np.sqrt(WF_norm)\n", " # higher energy --> higher offset in plotting\n", " WF_plot = WF_scale_factor*np.abs(psi[:,i])**2 + E[i] # also try plotting real part of WF!\n", " plt.plot(x_vals, WF_plot, ls=style_cycler.next(), lw=1.5, color=color_cyler.next(),\n", " label=\"$\\psi_{}(x)$\".format(i)) \n", " print(\"E[%s] = %s\"%(i, E[i]))\n", "\n", "plt.xlabel(\"$x$\")\n", "plt.legend(loc=\"best\")\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "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.12" } }, "nbformat": 4, "nbformat_minor": 0 }