{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An odd potential\n",
    "\n",
    "So far, the perturbing potentials we have applied have had even symmetry: here we will look at a simple potential with odd symmetry.  This is useful to illustrate more features of perturbation theory, and also to help you think about more aspects of quantum mechanics.  \n",
    "\n",
    "We start, as usual, with the basis set, integration routines, and then define a potential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import libraries and set up in-line plotting.\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as pl\n",
    "import numpy as np\n",
    "# This is a new library - linear algebra includes solving for eigenvalues & eigenvectors of matrices\n",
    "import numpy.linalg as la\n",
    "\n",
    "# Define the eigenbasis - normalisation needed elsewhere\n",
    "def eigenbasis_sw(n,width,norm,x):\n",
    "    \"\"\"The eigenbasis for a square well, running from 0 to a, sin(n pi x/a)\"\"\"\n",
    "    fac = np.pi*n/width\n",
    "    return norm*np.sin(fac*x)\n",
    "\n",
    "# We will also define the second derivative for kinetic energy (KE)\n",
    "def d2eigenbasis_sw(n,width,norm,x):\n",
    "    \"\"\"The eigenbasis for a square well, running from 0 to a, sin(n pi x/a)\"\"\"\n",
    "    fac = np.pi*n/width\n",
    "    return -fac*fac*norm*np.sin(fac*x)\n",
    "\n",
    "# Define the x-axis\n",
    "width = 1.0\n",
    "num_x_points = 101\n",
    "x = np.linspace(0.0,width,num_x_points)\n",
    "dx = width/(num_x_points - 1)\n",
    "\n",
    "# Integrate two functions over the width of the well\n",
    "def integrate_functions(f1,f2,size_x,dx):\n",
    "    \"\"\"Integrate two functions over defined x range\"\"\"\n",
    "    sum = 0.0\n",
    "    for i in range(size_x):\n",
    "        sum = sum + f1[i]*f2[i]\n",
    "    sum = sum*dx\n",
    "    return sum\n",
    "\n",
    "# Now set up the array of basis functions - specify the size of the basis\n",
    "num_basis = 10\n",
    "# These arrays will each hold an array of functions\n",
    "basis_array = np.zeros((num_basis,num_x_points))\n",
    "d2basis_array = np.zeros((num_basis,num_x_points))\n",
    "\n",
    "# Loop over first num_basis basis states, normalise and create an array\n",
    "# NB the basis_array will start from 0\n",
    "for i in range(num_basis):\n",
    "    n = i+1\n",
    "    # Calculate A = <phi_n|phi_n>\n",
    "    integral = integrate_functions(eigenbasis_sw(n,width,1.0,x),eigenbasis_sw(n,width,1.0,x),num_x_points,dx)\n",
    "    # Use 1/sqrt{A} as normalisation constant\n",
    "    normalisation = 1.0/np.sqrt(integral)\n",
    "    basis_array[i,:] = eigenbasis_sw(n,width,normalisation,x)\n",
    "    d2basis_array[i,:] = d2eigenbasis_sw(n,width,normalisation,x)\n",
    "    \n",
    "# Define a function to act on a basis function with the potential\n",
    "def add_pot_on_basis(Hphi,V,phi):\n",
    "    for i in range(V.size):\n",
    "        Hphi[i] = Hphi[i] + V[i]*phi[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A weak perturbation\n",
    "\n",
    "The potential will be very simple: a triangular shape at the bottom of the square well.  We'll start by making the magnitude of the potential small with respect to the lowest energy (remember that when we use words like \"small\" they have to be defined relative to an appropriate scale - in this case the energies of the unperturbed system)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFppJREFUeJzt3XGMnPWd3/H3tyZRFTWcDyGBwQbSFif0FDh8OqCXXDNp\ncsJ2KJyuh1HucqF3RDJQrlFLMFwTHRsR2jMBNT1FdSjJJe4/9UW+6M7nsAGLMroqIjQoa0MOG+NL\nkXaJ7UvSpqZV/vCKb/+YWfZhmLVn55l55pmZ90uymNl9dp4fj3Y/8+zvu9/fLzITSdL0+DujHoAk\nqVoGvyRNGYNfkqaMwS9JU8bgl6QpY/BL0pQpHfwRsTkijkTEyxFx7wrHNCJiLiK+HxHNsueUJPUv\nyvwdf0SsAV4CPgy8CnwX+GhmHi4csxb4NnB9Zi5ExPmZ+eNyw5Yk9avsHf81wLHMfCUzTwN7gJs6\njvkt4M8ycwHA0Jek0Sob/BcD84XnC+2PFV0OnBcRT0fEcxHxOyXPKUkq4ZySX9/LPNHbgE3Ah4B3\nAM9ExHcy8+WS55Yk9aFs8L8KbCg830Drrr9oHvhxZv4M+FlE/BVwFfCm4I8IFw2SpD5kZqzm+LJT\nPc8Bl0fEZRHxduAWYF/HMX8BvD8i1kTEO4BrgRe7vVhm+i+T+++/f+RjqMs/r4XXwmtx5n/9KHXH\nn5mLEXEX8ASwBvhKZh6OiO3tzz+amUci4lvA88DrwGOZ2TX4JUnDV3aqh8ycBWY7PvZox/OHgYfL\nnkuSVJ6duzXUaDRGPYTa8Fos81os81qUU6qBa5AiIusyFkkaFxFBVlzclSSNGYNfkqaMwS9JU8bg\nl6QpY/BL0pQx+CVpyhj8kjRlDH5JmjIGvyRNGYNfkqaMwS9JU6b06pySpOplwlNP9fe1Br8kjZnv\nfQ8+9Sk4fry/r3eqR5LGxPw8fPzj8JGPwC23wAsv9Pc6Br8k1dypU/DpT8Mv/iJccgkcPQrbt8M5\nfc7ZGPySVFOLi/ClL8G73w2vvgqHDsHnPgfvfGe513WOX5JqJhO++U3YsQPWrYPHH4errx7c6xv8\nklQjc3PLhdvPfx62boVY1f5aZ+dUjyTVwPw83HprK+h/8zfh+edbRdxBhz4Y/JI0UsXC7YYN8NJL\ncMcd/Rdue2HwS9IILC7Crl2wcSMsLMDBg63C7bnnDv/czvFLUoWWCrf33LNcuN20qdoxGPySVJGl\njtsTJ+Dhh4dTuO2FUz2SNGRLHbdbt8K2bcMt3PbC4JekIenWcXv77cMt3PbC4JekAess3C513FZR\nuO1F6eCPiM0RcSQiXo6Ie89w3C9HxGJE/EbZc0pSHWXC/v3w3vfC3r0wOwu7d8P69aMe2ZuV+oUj\nItYAXwQ+DLwKfDci9mXm4S7H7QS+BYxoVkuShqfYcTvKwm0vyt7xXwMcy8xXMvM0sAe4qctxvw/s\nBX5U8nySVCvFjts6FG57UTb4LwbmC88X2h97Q0RcTOvNYFf7Q1nynJI0csXC7fr1rY7bMkslV6ls\n8PcS4l8A7svMpDXNU+P3QUk6s25LJT/4YH0Kt70o+970KrCh8HwDrbv+ol8C9kTr957zgS0RcToz\n93W+2MzMzBuPG40GjUaj5PAkaTCKHbcXXTT4pZJ71Ww2aTabpV4jWjfifX5xxDnAS8CHgB8C/wP4\naGdxt3D8V4G/zMxvdPlclhmLJA1LseN2WEsl9ysiyMxVjabUVE9mLgJ3AU8ALwJ/mpmHI2J7RGwv\n89qSNGqde9yOQ+G2F6Xu+AfJO35JdXHqFOzc2ZrLv+OO1k5YdZ3Dr/yOX5ImSbHjtrjHbV1Dv19j\n8IdHkjRcnYXb2dnRFG6rYvBLmmrj1HE7KE71SJpKxY7bm2+enMJtLwx+SVPl1Cn4zGdaUzl1Wiq5\nSga/pKlQ7Lidn2/tcfvAA/DOd456ZNWbovc4SdOoLh23dWLwS5pYnR230zKHfzZO9UiaON06bm+4\nwdBfYvBLmhjFpZI3bBivpZKrZPBLGntLHbfjvFRylXwflDS2lva43bHDwu1qGPySxlKxcDstHbeD\n4lSPpLEyqUslV8nglzQWioXbSy9tddxauO2PwS+p1rrtcTutHbeD4nulpFqy43Z4DH5JtVNcKvmR\nR2DLFufwB8mpHkm10W2pZP9aZ/AMfkkj163jdtqWSq6SwS9pZIqF24WFyd3jtm58P5VUuaXC7Y4d\nsG6dhduqGfySKtW5VLJz+NVzqkdSJZYKt3bcjp7BL2moXCq5fgx+SUPRrePWwm09+J4raaDsuK0/\ng1/SwLhU8nhwqkdSaS6VPF5KB39EbI6IIxHxckTc2+Xzvx0RhyLi+Yj4dkRcWfackurBwu14KhX8\nEbEG+CKwGfhHwEcj4oqOw34A/JPMvBJ4APjPZc4pafQ6O24PHnSP23FS9n35GuBYZr4CEBF7gJuA\nw0sHZOYzheOfBdaXPKekESkWbtetaz3etGnUo9JqlQ3+i4H5wvMF4NozHH8b8HjJc0oaAQu3k6Ns\n8GevB0bEB4HfA9630jEzMzNvPG40GjQajRJDkzQI8/OtefwDB+D+++ETn3AOf5SazSbNZrPUa0Rm\nz9n91i+OuA6YyczN7ed/ALyemTs7jrsS+AawOTOPrfBaWWYskgbr1CnYubM1l3/nna0F1dzusH4i\ngsxc1e9eZf+q5zng8oi4LCLeDtwC7OsY1CW0Qv9jK4W+pPpYXIRdu2DjRve4nVSlfmHLzMWIuAt4\nAlgDfCUzD0fE9vbnHwX+EPh5YFe0JgRPZ+Y15YYtadA6O25nZ+24nVSlpnoGyakeaXTm5uDuu10q\neRyNYqpH0hgr7nFrx+30MPilKdTZcXv0qB2308Tgl6bISnvcWridLr6/S1PAjlsVGfzShCsWbu24\nFTjVI00sC7daicEvTZhi4Xb9epdK1lsZ/NKE6LbHrUslqxvvAaQx5x63Wi2DXxpjLpWsfjjVI42h\n4h6327ZZuNXqGPzSGCkWbi+5pNVxe/vtFm61Oga/NAY6C7cHD9pxq/55nyDV2FLhdseOVsethVsN\ngsEv1VSxcOtSyRokp3qkmil23Fq41TAY/FJNdFsq2cKthsHgl0asuMdtcalkO241LN5LSCPiHrca\nFYNfGgE7bjVKTvVIFSp23LpUskbF4Jcq0K3j1qWSNSoGvzRE3ZZKtuNWo+b9hjQEdtyqzgx+acDm\n5lqF2+PHW4XbLVucw1e9ONUjDUix4/bmm1uFW/9aR3Vk8EslnToFn/mMHbcaHwa/1Kdi4XZ+3qWS\nNT5KB39EbI6IIxHxckTcu8Ixf9z+/KGIsMSlsZYJ+/fDlVfC17/eKuLu3t2625fGQalfRiNiDfBF\n4MPAq8B3I2JfZh4uHLMV+IeZeXlEXAvsAq4rc15pVObm4O67XSpZ463sHf81wLHMfCUzTwN7gJs6\njrkR2A2Qmc8CayPigpLnlSq11HG7ZYtLJWv8lQ3+i4H5wvOF9sfOdsz6kueVKuEet5pEZb99s8fj\nOu+Lun7dzMzMG48bjQaNRqOvQUllLS7CY4/BZz8L11/f6rhd7+2KaqDZbNJsNku9RmT2mt1dvjji\nOmAmMze3n/8B8Hpm7iwc8yWgmZl72s+PAB/IzJMdr5VlxiINwlLhdseO1lLJDz9sx63qLSLIzFVN\nOpa9438OuDwiLgN+CNwCfLTjmH3AXcCe9hvFTztDX6oDl0rWtCg1x5+Zi7RC/QngReBPM/NwRGyP\niO3tYx4HfhARx4BHgTtLjlkaqKWOW5dK1rQoNdUzSE71qGqnTsHOna0mrDvugHvvtflK46efqR47\ndzV1XCpZ084/StPU6Nzj1qWSNa0Mfk0FC7fSMqd6NNHc41Z6K4NfE6mz4/all9zjVlpi8GuiLC7C\nrl2wceObC7fnnjvqkUn14f2PJkJn4XZ21sKttBKDX2PPwq20Ok71aGwVC7culSz1zuDX2HntteXC\n7aWXulSytFoGv8bGUsftxo2wsNAq3D7wgB230mp5j6Tay2x12d5zD6xbZ8etVJbBr1qbm2sVbo8f\nh4cecg5fGgSnelRLS0slb926XLi94QZDXxoEg1+1Uuy43bDBjltpGAx+1UKx43apcGvHrTQc3kdp\npIp73C4VbjdtGvWopMlm8Gtk7LiVRsOpHlXOpZKl0TL4VRmXSpbqweDX0K20x62FW2k0vNfS0LjH\nrVRPBr+GYqnj9sQJeOQR2LLFOXypLpzq0UB1dtweOuRf60h1Y/BrIDo7bo8etXAr1ZXBr1KKhdti\nx61LJUv15f2Y+rJUuC123Fq4lcaDwa9VK3bcfv7zzuFL46bUVE9EnBcRByLiaEQ8GRFruxyzISKe\njoi/jojvR8S/KnNOjc5S4daOW2m8lZ3jvw84kJkbgafazzudBv51Zv4CcB3wLyPiipLnVYVcKlma\nLGWD/0Zgd/vxbuDXOw/IzBOZebD9+P8Ch4GLSp5XFbDjVppMZe/ZLsjMk+3HJ4ELznRwRFwGXA08\nW/K8GiI7bqXJdtbgj4gDwIVdPvXp4pPMzIjIM7zO3wP2Ap9s3/m/xczMzBuPG40GjUbjbMPTgBX3\nuHWpZKl+ms0mzWaz1GtE5opZffYvjjgCNDLzRESsA57OzPd0Oe5twH5gNjO/sMJrZZmxqJyFhdY8\n/pNPwswM3Habc/jSOIgIMnNVt2dl5/j3Abe2H98K/HmXQQXwFeDFlUJfo/Paa63Av+oqO26laVE2\n+P8I+LWIOAr80/ZzIuKiiPhm+5j3AR8DPhgRc+1/m0ueVyUtFW43bnxz4daOW2nylZrqGSSneqpR\nLNyuW9daOdPCrTS++pnq8Rf6KeIet5LARdqmQnGP223b7LiVpp3BP8G67XF7++0WbqVpZ/BPoMVF\n2LVruXB78KAdt5KWee83QTo7bmdnLdxKeiuDf0J0LpXsHL6klTjVM+aKhdulpZJvuMHQl7Qyg39M\nFQu3l15qx62k3hn8Y6bbUskPPGDHraTeeX84JlwqWdKgGPxjYG4O7r7bjltJg+FUT40t7XG7dat7\n3EoaHIO/hoqF2/Xr3eNW0mAZ/DXSrXD74IN23EoaLO8ha8DCraQqGfwj5lLJkqrmVM+IuFSypFEx\n+CvWuVTy0aMulSypWgZ/RToLt0tLJdtxK6lq3mcOWWarWLu0x62FW0mjZvAPUedSyRZuJdWBUz1D\nUOy4tXArqW4M/gEqFm43bLBwK6meDP4BKO5xOz/vHreS6s170RKKHbdLhdtNm0Y9Kkk6M4O/T0uF\n2+PH7biVNF6c6lmlzj1uX3jBwq2k8WLw96izcOtSyZLGVd/BHxHnRcSBiDgaEU9GxNozHLsmIuYi\n4i/7Pd+oFAu3LpUsaRKUueO/DziQmRuBp9rPV/JJ4EUgS5yvUpmwfz+8972wdy/MzsLXvtbaGEWS\nxllk9pfFEXEE+EBmnoyIC4FmZr6ny3Hrga8BDwL/JjP/2Qqvl/2OZdDm5t5cuN2yxTl8SfUUEWTm\nqhKqzB3/BZl5sv34JHDBCsf9B+Ae4PUS56pEseP25ptbHbf+tY6kSXPG0mREHAAu7PKpTxefZGZG\nxFtu1yPiBuBvM3MuIhplBjpMp07BQw+15vLvvLPVceuqmZIm1RmDPzN/baXPRcTJiLgwM09ExDrg\nb7sc9ivAjRGxFfi7wLkR8V8y8+PdXnNmZuaNx41Gg0ajcfb/gxIWF+HLX4bPfhauv75VuHUOX1Kd\nNZtNms1mqdcoM8f/EPCTzNwZEfcBazNzxQJvRHwA+FQd5viXOm537IALL4RHHnGpZEnjqZ85/jJ/\nhf5HwNcj4jbgFWBbexAXAY9l5ke6fM3Iq7dzc3D33S6VLGl69X3HP2jDvuOfn281YD35JMzMwCc+\nYfOVpPFX9V/1jAX3uJWkN5vY4C/ucbuw0CrculSyJE3g6pzFpZIvusg9biWp00QFf7Fw61LJktTd\nREz1FDtub7nFPW4l6UzGOvi77XHrUsmSdGZjGfzFwu3SUsmf+5zLLEhSL8bq3tjCrSSVNzbBv7TH\nrYVbSSqn9lM9S3vcbt0K27ZZuJWksmob/HbcStJw1C74O/e4PXjQjltJGqRa3T/v379cuJ2dtXAr\nScNQq9U5r7giXSpZklahn9U5axX8p0+nc/iStApjvyyzoS9Jw1er4JckDZ/BL0lTxuCXpClj8EvS\nlDH4JWnKGPySNGUMfkmaMga/JE0Zg1+SpozBL0lTxuCXpClj8EvSlOk7+CPivIg4EBFHI+LJiFi7\nwnFrI2JvRByOiBcj4rr+hytJKqvMHf99wIHM3Ag81X7ezX8EHs/MK4ArgcMlzjkVms3mqIdQG16L\nZV6LZV6LcsoE/43A7vbj3cCvdx4QET8H/Gpm/glAZi5m5v8pcc6p4Df1Mq/FMq/FMq9FOWWC/4LM\nPNl+fBK4oMsx7wJ+FBFfjYjvRcRjEfGOEueUJJV0xuBvz+G/0OXfjcXjsrWNV7etvM4BNgH/KTM3\nAf+PlaeEJEkV6HvrxYg4AjQy80RErAOezsz3dBxzIfBMZr6r/fz9wH2ZeUOX16vHHpCSNGZWu/Vi\nmc0O9wG3Ajvb//3zLoM5ERHzEbExM48CHwb+utuLrXbgkqT+lLnjPw/4OnAJ8AqwLTN/GhEXAY9l\n5kfax10FfBl4O/A3wO9a4JWk0ek7+CVJ46nyzt2I2BwRRyLi5Yi4d4Vj/rj9+UMRcXXVY6zK2a5F\nRPx2+xo8HxHfjogrRzHOYevle6J93C9HxGJE/EaV46tSjz8fjYiYi4jvR0Sz4iFWpoefj/Mj4lsR\ncbB9Lf7FCIZZiYj4k4g4GREvnOGY3nMzMyv7B6wBjgGXAW8DDgJXdByzlVbDF8C1wHeqHGPNrsU/\nBn6u/XjzJF6LXq5D4bj/BuwH/vmoxz3C74m1tOpk69vPzx/1uEd4LWaAf790HYCfAOeMeuxDuh6/\nClwNvLDC51eVm1Xf8V8DHMvMVzLzNLAHuKnjmDcawzLzWWBtRHTrERh3Z70WmflMLtdDngXWVzzG\nKvTyPQHw+8Be4EdVDq5ivVyL3wL+LDMXADLzxxWPsSq9XIvjwLntx+cCP8nMxQrHWJnM/O/A/z7D\nIavKzaqD/2JgvvB8of2xsx0ziYHXy7Uoug14fKgjGo2zXoeIuJjWD/2u9ocmtTDVy/fE5cB5EfF0\nRDwXEb9T2eiq1cu1eAz4hYj4IXAI+GRFY6ujVeVmmT/n7EevP7Cdf9o5iT/oPf8/RcQHgd8D3je8\n4YxML9fhC7T6PzIigrd+f0yKXq7F22g1RX4IeAfwTER8JzNfHurIqtfLtfi3wMHMbETEPwAORMRV\nmfnakMdWVz3nZtXB/yqwofB8A613pjMds779sUnTy7WgXdB9DNicmWf6VW9c9XIdfgnY08p8zge2\nRMTpzNxXzRAr08u1mAd+nJk/A34WEX8FXAVMWvD3ci1+BXgQIDP/JiL+J/Bu4LlKRlgvq8rNqqd6\nngMuj4jLIuLtwC20GsGK9gEfB2gv4fzTXF4TaJKc9VpExCXAN4CPZeaxEYyxCme9Dpn59zPzXdnq\nAN8L3DGBoQ+9/Xz8BfD+iFjTXvfqWuDFisdZhV6uxRFaTaG057PfDfyg0lHWx6pys9I7/sxcjIi7\ngCdoVe2/kpmHI2J7+/OPZubjEbE1Io7RWtvnd6scY1V6uRbAHwI/D+xq3+2ezsxrRjXmYejxOkyF\nHn8+jkTEt4DngddpNUtOXPD3+H3x74CvRsQhWjexOzLzf41s0EMUEf8V+ABwfkTMA/fTmvbrKzdt\n4JKkKePWi5I0ZQx+SZoyBr8kTRmDX5KmjMEvSVPG4JekKWPwS9KUMfglacr8f9pB/XNxLB43AAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10876ef50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the potential in the square well\n",
    "def square_well_potential(x,V,a):\n",
    "    \"\"\"Potential for a particle in a square well, expecting two arrays: x, V(x), and potential height, a\"\"\"\n",
    "    for i in range(x.size):\n",
    "        V[i] = a*(x[i]-width/2.0)\n",
    "    # Plot to ensure that we know what we're getting\n",
    "    pl.plot(x,V)\n",
    "    \n",
    "# Declare space for this potential (Vdiag) and call the routine\n",
    "Vdiag = np.linspace(0.0,width,num_x_points)\n",
    "square_well_potential(x,Vdiag,1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will build the matrix elements of the potential and the overall Hamiltonian.  Look carefully at the potential, and ask yourself if the form of the matrix makes sense: how does the symmetry affect the form ? What effect will this have on the energies of the perturbed system ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Full Hamiltonian\n",
      "   4.935   -0.180   -0.000   -0.014    0.000   -0.004    0.000   -0.002    0.000   -0.001\n",
      "  -0.180   19.739   -0.195   -0.000   -0.018    0.000   -0.006   -0.000   -0.002    0.000\n",
      "   0.000   -0.195   44.413   -0.199    0.000   -0.020   -0.000   -0.006    0.000   -0.003\n",
      "  -0.014   -0.000   -0.199   78.957   -0.200   -0.000   -0.021   -0.000   -0.007   -0.000\n",
      "  -0.000   -0.018   -0.000   -0.200  123.370   -0.201    0.000   -0.021   -0.000   -0.007\n",
      "  -0.004    0.000   -0.020   -0.000   -0.201  177.653   -0.201   -0.000   -0.022    0.000\n",
      "  -0.000   -0.006   -0.000   -0.021    0.000   -0.201  241.805   -0.202    0.000   -0.022\n",
      "  -0.002    0.000   -0.006   -0.000   -0.021   -0.000   -0.202  315.827   -0.202    0.000\n",
      "   0.000   -0.002    0.000   -0.007   -0.000   -0.022    0.000   -0.202  399.719   -0.202\n",
      "  -0.001    0.000   -0.003   -0.000   -0.007    0.000   -0.022    0.000   -0.202  493.480\n",
      "Perturbation matrix elements:\n",
      "   0.000   -0.180    0.000   -0.014   -0.000   -0.004    0.000   -0.002    0.000   -0.001\n",
      "  -0.180    0.000   -0.195   -0.000   -0.018   -0.000   -0.006    0.000   -0.002   -0.000\n",
      "  -0.000   -0.195    0.000   -0.199   -0.000   -0.020   -0.000   -0.006    0.000   -0.003\n",
      "  -0.014   -0.000   -0.199   -0.000   -0.200   -0.000   -0.021   -0.000   -0.007   -0.000\n",
      "  -0.000   -0.018   -0.000   -0.200    0.000   -0.201    0.000   -0.021    0.000   -0.007\n",
      "  -0.004   -0.000   -0.020    0.000   -0.201   -0.000   -0.201    0.000   -0.022    0.000\n",
      "  -0.000   -0.006   -0.000   -0.021    0.000   -0.201    0.000   -0.202    0.000   -0.022\n",
      "  -0.002    0.000   -0.006   -0.000   -0.021   -0.000   -0.202    0.000   -0.202    0.000\n",
      "   0.000   -0.002    0.000   -0.007    0.000   -0.022    0.000   -0.202    0.000   -0.202\n",
      "  -0.001   -0.000   -0.003   -0.000   -0.007    0.000   -0.022    0.000   -0.202    0.000\n"
     ]
    }
   ],
   "source": [
    "# Declare space for the matrix elements\n",
    "Hmat2 = np.eye(num_basis)\n",
    "\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "print \"Full Hamiltonian\"\n",
    "for n in range(num_basis):\n",
    "    # Loop over basis functions phi_m (the ket in the matrix element)\n",
    "    for m in range(num_basis):\n",
    "        # Act with H on phi_m and store in H_phi_m\n",
    "        H_phi_m = -0.5*d2basis_array[m] \n",
    "        add_pot_on_basis(H_phi_m,Vdiag,basis_array[m])\n",
    "        # Create matrix element by integrating\n",
    "        Hmat2[m,n] = integrate_functions(basis_array[n],H_phi_m,num_x_points,dx)\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % Hmat2[m,n],\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print\n",
    "    \n",
    "print \"Perturbation matrix elements:\"\n",
    "# Output the matrix elements of the potential to see how large the perturbation is\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "for n in range(num_basis):\n",
    "    # Loop over basis functions phi_m (the ket in the matrix element)\n",
    "    for m in range(num_basis):\n",
    "        # Act with H on phi_m and store in H_phi_m\n",
    "        H_phi_m = np.zeros(num_x_points)\n",
    "        add_pot_on_basis(H_phi_m,Vdiag,basis_array[m])\n",
    "        # Create matrix element by integrating\n",
    "        H_mn = integrate_functions(basis_array[n],H_phi_m,num_x_points,dx)\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % H_mn,\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we diagonalise, as usual.  Look at the eigenvalues: do you understand why they take these values ? Is this simple, first order perturbation expansion good enough for this problem ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Eigenvalues and eigenvector coefficients printed roughly\n",
      "[   4.9326   19.7399   44.4136   78.9571  123.3702  177.653   241.8054\n",
      "  315.8274  399.719   493.4807]\n",
      "[  9.9993e-01   1.2165e-02  -3.3858e-05  -1.9471e-04  -6.9229e-07\n",
      "   2.2987e-05  -6.0505e-08   5.2552e-06   1.0167e-08  -1.6928e-06]\n",
      "[  1.2165e-02  -9.9989e-01   7.8840e-03   1.8058e-05   1.7738e-04\n",
      "  -5.2925e-07   2.5236e-05  -5.7217e-08  -6.4761e-06   1.6305e-08]\n",
      "[  6.0935e-05  -7.8838e-03  -9.9995e-01  -5.7465e-03  -1.0818e-05\n",
      "   1.5022e-04  -3.8187e-07   2.3694e-05   4.6815e-08  -6.5389e-06]\n",
      "\n",
      "    Diag  Perf Square  Difference\n",
      "   4.933        4.935      -0.002\n",
      "  19.740       19.739       0.001\n",
      "  44.414       44.413       0.000\n",
      "  78.957       78.957       0.000\n",
      " 123.370      123.370       0.000\n",
      " 177.653      177.653       0.000\n",
      " 241.805      241.805       0.000\n",
      " 315.827      315.827       0.000\n",
      " 399.719      399.719       0.000\n",
      " 493.481      493.480       0.000\n"
     ]
    }
   ],
   "source": [
    "# Solve using linalg module of numpy (which we've imported as la above)\n",
    "eigval, eigvec = la.eigh(Hmat2)\n",
    "# This call above does the entire solution for the eigenvalues and eigenvectors !\n",
    "# Print results roughly, though apply precision of 4 to the printing\n",
    "print\n",
    "print \"Eigenvalues and eigenvector coefficients printed roughly\"\n",
    "np.set_printoptions(precision=4)\n",
    "print eigval\n",
    "print eigvec[0]\n",
    "print eigvec[1]\n",
    "print eigvec[2]\n",
    "print\n",
    "\n",
    "print \"    Diag  Perf Square  Difference\"\n",
    "for i in range(num_basis):\n",
    "    n = i+1\n",
    "    print \"%8.3f     %8.3f    %8.3f\" % (eigval[i],n*n*np.pi*np.pi/2.0,eigval[i] - n*n*np.pi*np.pi/2.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we plot the eigenstates, with a small potential it's hard to see the change relative to the unperturbed eigenstates, so we'll also plot the difference between perturbed and unperturbed eigenstates.  What do these differences look like ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAADSCAYAAACxUP3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4VMX3h9+B0DsSQgsE6R3pYKHYwIK9N+wd/VrADj/E\nLhZEsSEICHZREAsoIAkQeoBQkkCAEAikV1L3/P6YG1iWlM1mN5sN8z7PfZK9d+7cue3cz8ycOaNE\nBIPBYDAYDAaDwVA81bxdAIPBYDAYDAaDobJjRLPBYDAYDAaDwVAKRjQbDAaDwWAwGAylYESzwWAw\nGAwGg8FQCkY0GwwGg8FgMBgMpWBEs8FgMBgMBoPBUApGNHsJpdRtSqm/vF0OX0ApdY1SKkYpla6U\n6uPB44xQSsV4Kn934SvlNBiqCkqpmUqpl+x+P6yUOqqUSlNKNVFKnauUirRs1FhvlrWqUFmuqVJq\nh1LqAm8d31C5MKLZgyil9iulsqyXvnCZDiAi34jIpd4uoysopcYppVZX4CHfBR4RkQYiElaBxz3j\nUUoFKaVsSiljKwxVEjs7naaUSlZKhSilHlRKqcI0IvKwiEy10tcApgEXikhDEUkGpgDTLRv1m3fO\npHxY12GUt8thh1uuaXnPS0R6ish/ru7vTSrhPfV5/LxdgCqOAFeIyL/eLkhlQilVXUQKnEyrgLbA\nThePVU1EbK7sazgFVXoSg8EnOWGnlVINgBHAh8Bg4J4i0rcAagO77NaVx0Y5bQ89jFAJ3nOllJ+I\n5FOOa+pApTgvL+HyudvdB4M9ImIWDy1ANDCqmG3jgNV2vy8B9gApwMfAKuBeu+33oA1IEvAn0NZu\nmw14EIgAkoEZ1vpaVn497NL6A1lAM+v3FcBWa78QoJdd2kDgZ+AYkAB8BHQFsoF8IB1IstI2AuZa\nafcDLwLK7lxDgPesfKYAHa1zTAHigW+LuEa1gAzr/DKASGt9N2ClVeYdwJV2+8wBZgJLrX1Ou/5A\nU2A2EGtdz1+s9SOAGOAp4ChwGBhnt9/lwBYgFTgITLLbFmSV807ggHVOL9htrwN8bR1vJzABiLHb\n3gr4ybp++4DHHfadY+0bDjxrv28R5/e+Vf5UYBvQAxgIxBXeEyvdtcBW6/9BwEZrnzjgXWv9Qeu8\n0q1lsJPP48NAJJBm3e8OwFrrfn8L1PD2+2kWs4gUbaet96UA6G79ngO8CnQCMu3eiX+AKCttlvW8\n10Dbw1mWDTlk7VvNymscp9vDmugetQPW+zcTqG2lH2HlUZxdqoNu+d5vvV+r7fYdAqxB28qtwPBi\nrsE8u3NIB56x1o+1bE4ysALoWsJ1tAGPA3vR9u9tB3tTms14BP0N21fWa2rlcb+Vf5pV5nOKO68i\nyl7Sd3A/uleh8Fq7ascnA99b+6ehv139rW0TgR8cyvQh8KH1v9vOvaR7ap3rBPR34zhQ3SrbISvv\n3RSjac6UxesFqMoL2hhfWMy2cViiGWiGFitXo11mxgO5wD3W9qvQAqSLtf1FIMQuLxvwG9AQLXSP\nAZda22YBU+3SPgostf4/B22EB6Jro3daZa5hvSxhaGNcBy1gh1n73YWd4LfWzQV+AeoB7dAVgHvs\nzjXPOnY1dCvNQuB5a3vNwryLuVY24Gzr/xpog/ocuqdkpPUyd7a2z0F/OIZav2sVkd/v1vEbWXmc\nb60fYZVzsnX+Y9AfyEbW9uFYFRCgF/rjdpX1O8gq52fWteqNrlx0sba/iTZQjYDWaKN00NpWDdgE\nvGSVpz36w3OJ3b6rgMZAG7SxPVjMtboULX4bWr+7AC2s/8OB0XZpfwH+Z/2/FrjN+r8uJ8VxO+u8\n7A20M8/jL0B9oDuQA/xrXaOGVjnu9Pb7aRaziBTfuIEWsA9a/88Gplj/F/VOnJKH9fzPRNtOfyAU\neMDaNo7T7eH7wCLrHa+PtuevW+lLs0sfW+9XSyu/IWib2hotykdb6S6yfjdz5joAndENDxdax33W\neu+LrPBa1+Qf6xwC0d+Ae61tztiMv6x9a7lwTW9AC7tCEdoBS5QXd3/t8i32O+i4P+Wz45PRQnS0\ndZzXgbV2z1QmUN/6XR0tkAe589xLuKd+1vb9wGbr3GpZ9+sgJ78hbbG+xWfq4vUCVOXFegDT0TW6\nwqXQiIzjpGi+096AWOsOclJ0/lH4v/W7mvWCBVq/bdiJTuA7YKL1/4VAlN22EOB26/+ZWB8Cu+27\ngQuAoWjxXa2I8zpRdut3dbQwsq+xPgCssEt/wCGPr9ECs7UT19FeNJ8PHHHYvgCr1RctmueUkFdL\ndO27URHbRqBr5fYfw6OFhquI9B8A71n/B1nlbGW3PRS40fp/L3Cx3bZ7sVoo0N3AjtfneeAru30v\nsdt2P8W0NKMrEXusPKs5bJsIzLf+b2o9QwHW71Voo97MYZ/C87K/Js48j0Pttm8EnrX7/S7wvjff\nTbOYpXCheNG8lpMV+9nAq9b/Rb0T9sIqAF1hrm23/RbgX+v/U+whWkBlYCdG0PZ3n/V/sXbJevey\nsGsZtUszEZjrsO5PiqmwOl4H4GXsegCtch6i+NZqm4OdehhYbv3vjM0YUVx5nLimf2HXquvM/bXb\nXtx38PwiylEeOz4Z+NtuW3cgy+73auAO6/+Lsb7b7jz3Eu7pBXbpx9lt72g9axdiegcRETMQ0MMI\nuiWyid0yq4h0rdAPrj32v9sBH1qDVJKBRGt9a7s0cXb/Z6FbK0C7MdRVSg1SSgUBfdC11sJ8ny7M\n18q7DVpYBqINgDP+wM3QLcAH7NYddCifY7SHCegXdr01OvluJ44D+lo55nXAWg/6mpcUWSIQ7VKS\nWsz2RIdzPnEtlVKDlVIrlFLHlFIpaJeYsxz2L+4+OJbb8f62crgPzwPNi9n3YHEnJyIrgBno1qej\nSqnPLD9NgG+AK5VSdYEbgf9E5Ki17V50K8QupdR6pdTlxR0D557Ho3b/Hy/id30MhspNG3Q3fFlp\nh7aHR+zekU/RLYSF2L/P/ujenU126f9A29VCirNLzdAt1XuLKccNDnblXLRPtjO0xM7WiFZRMZy0\ntUXhaKcK0zpjM0qy26Vd0zYUfQ2cobjvYFHnWR47DqfawSygtt0g6wVoMQxwK9peF+brrnMv7p4W\neR9EJAp4Ei34jyqlFiqlWjp5rCqJGQhYOTgMXFn4wxr81sZu+0F0K8fCsmYsIgVKqe/RL+MxYLGI\nZNrl+5qIvO64n1JqKNC2mEEq4vA7Ad19GMTJwTFtOdWgnLKPJdYesI51LrBcKbVKRPaVckqHgUCl\nlLJeeNBGZXcp+xUSAzRVSjUqQTgXxwJgOtr1JVcp9T6nfthK4ghasBeWM9ChTNEi0rmEfdty6rUt\nFhH5CPhIKeWP9qF7FnhFRA4ppdahfZlvBz6x2ycKbahRSl0H/KiUasrp9xrK8TwaDL6AUmogWiAF\nu7B7DLrn7awSGh3s36sEdEWyu4gcKeOxEtCtkB3RrgL2HATmicgDTubl+K4fRruhASe+S4HosSDF\n4WinCtM6YzOKsjWFlHZNY9DXoKz5FpatyO9gEZTHjpdWjh+BaUqp1mhXzSF2+brr3J25p47f6oXA\nQqvx5TPgLXTv+BmJaWn2PM6MXF0K9FJKXaWU8kP7utm3BnwKvKCU6g6glGqklLqhDMdcANyMFkUL\n7NZ/ATxktUIrpVQ9pdTlSqn6aNeCI8CbSqm6SqnaSqlh1n5HgTZW6CUsUf098JpSqr5Sqh3wP2B+\nsQVU6galVGHFIAX9ojrTqr0OXUOfoJSqoZQagR7E8W0x534K1kfpD+ATpVRjKw9nY3DWB5ItwTwI\nfT1LM4SFfA88bx2zNfCY3b7rgXSl1ASlVB2lVHWlVE+l1IAi9m2DHmxTJEqpAVaLeA30dcpGu6MU\nMhfdbdsTPcizcL/bLZEN2r++8H7EW3872OVR1ucRTr0vZ+pIdkPlRQEopRoqpa5Aj3mYJyLh9tud\nwbIxfwPvKaUaKKWqKaU6FGdnLCH0BfBB4TuolGqtlLrEiWPZgK+sY7W0bMdQpVRNtP29Uil1ibW+\nttIx3lsXk91RTn3PvwcuV0qNsuzJ02h7sqaEIj1j2alA9Nic76z1rtgM+/Ms7Zp+aR27n/Ut66iU\nKmxccDwvR0r6DjpSHjte2rcpHt0zPAftmrPHA+depnuqlOpspa2FFu6O35MzDiOaPc9idWqc5p+s\n9WItiEgC2pn/bXTLQTe0H2iOtX0Runb3rVIqFdiOHvCFXV44/D6xTkTWo33mWqIFY+H6TWj/2Bno\nbshIrBqkZYyvRNdgD6Jrszdau/6DHswVp5Q6Zq17HO2jtg/tm/UN2g/wtPJYDADWKaXSgV+B8SKy\nv6gL6HAueVa5xqAF3Qy0H1hECcdy5A50y/hutFEZX9SxiuARYIpSKg3tG/adw/aS9p2CbnmPRhvA\nH9CDPQsrHVcAfdHXLx74HD1oDuD/0C4o0WifxLklHKuhtW8S2qc+AXjHbvvP6BagX0Qk2279pcAO\n6368D9wsIjkikgW8BoQo3TU4yIXn0XGdM/fIYKhIFlvv9UF0l/o0wN5lzPGZLe35vRM9GK8wWsQP\nnGwIKer5n4ge4LzOeqeWod2lnDneM+h3cAPa7eENtP/zIfQAvBfQvYwH0SKpuO/+G8BL1nv+lGVT\nb0dHTYpHRw+6UkoOQ/YrejDcFmAJWtC78g0rimKvqYj8iLZTC9ADw38GmhR1Xo6ZlvAdLKpM5bHj\nRd13x98L0P7DCxzWu+XcXbintaw84tGNaM3Q78cZS2FIMNd21rXJuWifHQE+F5HpRaSbjhY5WWgn\n8y0uH/QMQGkfpxjgVhFZ5e3yGNyPUuph9CDBkV44diQ6KoCJH24wGNyCUsoGdHTCxa7K4E07bvAO\n5W1pzkOHrOqB9r95VCnVzT6BUuoy9IvUCe3DOrOcx6ySWF1oja1ukBes1eu8WSaD+1BKtVB6Wthq\nSqku6Jirv5S2nwfKcS16/IcRzAaDwVAGKosdN3iPcg0EFJE4rGgBIpKhlNqFHjxhP1PSWHR4MUQk\n1BKGAXaj9g2aoejulZpo14erRSTHu0UyuJGaaL++9mgf7oXYDcSrCJRSK9GT09xRkcc1GAxnBGeC\ny5XX7bjBu5TLPeOUjHQ4s1XoyR8y7NYvBt4QkTXW7+XoGMKb3HJgg8FgMBgMBoPBw7gl5Jw1yvRH\n4Al7wWyfxOH3aUpdKXUm1FINBkMVRUTOqKggxmYbDAZfxhWbXe7oGVbYkp/QM40tKiJJLKfGMmxD\nMXEeK3JWl8qwTJo0yetlMOdrztmcc/mXMxVvX3fzbJtzNudszteVxVXKJZqVUgqYBewUkQ+KSfYb\nVhgzpdQQIEWMP7PBYDAYDAaDwYcor3vGueiYf9uUUoVh5F7AmrFMRD4TkaVKqcuUUlHoOL7OTpds\nMBgMBoPBYDBUCsobPSMYJ1qrReSx8hynqjJixAhvF6FCOdPOF8w5GwxViTPx2TbnXPU50863PLgt\nekZ5UUpJZSmLwWAwlAWlFHIGDgQ0NttgMPgirtpsM422wWAwGAwGg8FQCkY0GwwGg8FgMBgMpWBE\ns8FgMBgMBoPBUApGNBsMBoPBYDAYDKVgRLPBYDAYDAaDwVAKRjQbDAaDwWAwGAylYESzwWAwGAwG\ng8FQCkY0GwwGQxVEKTVaKbVbKRWplJpYTJrp1vYwpdQ5DtuqK6W2KKUWV0yJDQaDoXJjRLPBYDBU\nMZRS1YEZwGigO3CLUqqbQ5rLgI4i0gl4AJjpkM0TwE7AzGBiMBgMGNFsMBgMVZFBQJSI7BeRPOBb\n4CqHNGOBrwFEJBRorJQKAFBKtQEuA74EzqiZDg0Gg6E4jGg2GAyGqkdrIMbu9yFrnbNp3geeBWye\nKqDBYDD4Gn7lzUAp9RVwOXBMRHoVsX0E8Cuwz1r1k4hMLe9xDQaDwVAszrpUOLYiK6XUFWh7vsWy\n38UyefLkE/+PGDGCESNKTG4wGAxeYeXKlaxcubLc+SiR8rmrKaXOBzKAuSWI5qdEZGwp+Uh5y2Iw\nGAzeQCmFiFQaNwal1BBgsoiMtn4/D9hE5C27NJ8CK0XkW+v3bmAEMB64A8gHagMN0Y0ddzocw9hs\ng8Hgk7hqs8vtniEiq4HkUpJVmo+JwWAwnAFsBDoppYKUUjWBm4DfHNL8BtwJJ0R2iojEicgLIhIo\nIu2Bm4F/HQWzwWAwnImU2z3DCQQYppQKA2KBZ0RkZwUc12AwGM5IRCRfKfUY8BdQHZglIruUUg9a\n2z8TkaVKqcuUUlFAJnB3cdlVTKkNBoOhclMRonkzECgiWUqpMcAioHNRCY1/nMFg8AXc5R/nSUTk\nD+APh3WfOfx+rJQ8VgGr3F86g8Fg8D3K7dMMoJQKAhYX5dNcRNpooL+IJDmsN/5xBoPBJ6lsPs0V\ngbHZBoPBV/GaT3NpKKUClFLK+n8QWqgnlbKbwWAwGAwGg8FQaXBHyLmFwHCgmVIqBpgE1IATXYHX\nAw8rpfKBLPTAkkqLTYQN6elsTk9HAdWV4qwaNbioSRMa+hVxubKyYNs22LMHcnIgLw+Ugg4doFs3\naNMGqrlWNymwFXAs8xiH0w9zJOMImbmZ5BbkkluQS22/2jSs1ZCGtRrSqkEr2jVuR83qNct38pUU\nEcjIgIQESEyEpCRITYX0dL1kZcHx45CdrS9/fj4UFOj9lNKLnx/UrAm1akHt2lCvnl4aNoTGjfXS\npAn4++v/XbxlpXK8oICVKSmsS0tDATWrVaNB9eqMbtqUTnXreuaglYGEBAgOhn379HLsmL4Rdevq\nC9+/PwwaBIGB+oYZDAaDwVDJcIt7hjvwdlfflvR0psfGsjQxEf8aNRjaqBEKLaIP5eSwJi2N8xs1\n4taAAG5JTKTa/PmwaBHs36/FcbduWgD4+WnFFhUFu3drtXfZZXDDDTB6tE5TBAW2AsKOhrFq/yq2\nxG1h+7Ht7E7YTaNajWjZoCUt67ekQa0G1KxekxrVapBTkENaThqp2anEpsdyKO0QAfUC6Nm8JwNb\nDWRg64GcG3guTeo0qdDr6AppaVpHRUfrJSYGDh3SS1wcHD2qRWyzZtC0qV4aNdKCt0EDfUnr1NEa\nrGZNqF5dL0pp4Syib0lOjl6ysyEzUy/p6ZCSAsnJWozHx+v1/v7QqhW0bAmtW0O7dtC2LQQF6fpQ\nQEDZtF1IaiqvHTjA6tRUzqlfnwsaNcJPKXJFiM/LY0liIs1r1OCm5s0Z37o19YuqoPkaR47A3Lmw\neDFs3w7nngudOkH79tCihb4ZWVn6om/cCKGh+gbeeivcfTd07ertM3Aa455hMBgMvoOrNvuMF83J\neXm8HB3ND/HxTGjblmubNaN9nTqnpUvLzWXpn38yPSkJcnP5NC6O3pddBn36QI0axR/g6FEtrn/8\nUQuDBx6AZ54Bf39Ss1NZHLGYn3f9zMr9K2lRvwUjgkYwoNUAejXvRXf/7tSrWc+p88i35ROTGkPY\n0TA2xG4gNDaU9bHr6RXQizEdx3B116vp2bynq5ep3IhAbCyEh+tl506IiNBLRobWUe3ba1Hatq1u\ncGzdWovWgACoX7/iypqbqxtCjxyBw4e1eD94EA4c0HWkvXt1y3aHDlrXFS49e0KXLlr3FRKbk8PE\nvXtZlZrKlKAgrmnWjMZFPC8FIqxJTeXTw4cJTk3l406duKJZs4o7aXcSFQXvvAM//ADXXw/XXgsj\nRuhaTUmIwK5dMGcOzJunL/Crr8LIkRVR6nJhRLPBYDD4DkY0u8AfiYncs2cPVzdrxmvt29O0KPEr\nAsuWwfPPA2CbMoVZffrw4oED3NuyJVPbt6e6s02OMTHY3nid/AXfsOT85jze7yj9u47kum7XcUmH\nS2jZoKUbzw6y87P578B/LI1cyk+7fuKsOmdxe+/bua3XbW4/lj02m9ZNGzbApk2wdSuEhelG+J49\noUcP6N5dC8wuXbQw9rUe+bQ0iIzUXjm7d+tKwI4dWlh36gR9+0K1EfH8cvYeHmjRikkd2zrdevxP\ncjIPRUTQt359ZnXpUrRbUGUkNVW/J99/D488Ao8/rpvsXSEvT1c0X3xR9+K89ZZ+eCopRjQbDAaD\n72BEcxmZdeQIL+7bx489enBe48ZFJ0pKgocfhs2b4fXX4brrTji7HsvN5ZadO2ni58c33btTqxQn\n2ISsBD5e/zFfbP6CfgXNeTekHh23HKTaF1/AJZe4+/ROwyY2/jvwH/PC5vHz7p+5rNNlPDH4CQa1\nHlTuvFNTYe3ak8v69dpNdcAAvfTtqxvkW7Rww4lUco4f1wL6nX0xLKkXw9mzerH3zwa0b69ddocN\n00vXriX7TWcXFPBkVBQb0tP5o3dvmtesxP7qIvDLLzB+vHZFeust/QC4g5wc+PRTeO01ePRRLaIr\nYSXCiGaDwWDwHYxodhIRYeqBA8yOi+OP3r3pUtzgq7/+gnvv1b7Ir7+unWYdyLHZuH3XLhLz8ljU\ns2eRLYKH0g4xbc00vg77muu6Xcfjgx+nd0BvvfHvv+H++7Vo/uADPTKtAkjJTmHW5lnM2DCDNg3b\nMHn4ZEa1H4Vysrk3JQVWrYKVK+G//3Rr64ABWgwOHQpDhrjewOjr2ER4eu9e/k5KYmnv3rSrXZvc\nXN0KHRqqKxUhIdqH+rzzYPhwuOAC6NdP+2HbIyJM3r+fhceO8Xfv3gQV8Qx6nePH4aGH9Ml9/rk+\nGU9w+DDcdZd2OJ8/H84+2zPHcREjmg0Gg8F3MKLZSV7ct4+lSUn80asXLWrVOj2BzQaTJsHXX8Ps\n2XDhhSXmVyDCY5GRbEhLY0XfvjSwhHNiViKvr36dOWFzGNdnHE8NfYrWDVufnkFamu7GDguDX3/V\nI84qiAJbAQt3LGTKqim0qN+C10a9xvntzj8tXV6eFnt//62XXbu0OB45Uou+/v11VAoDPBUVxfq0\nNBb36kWTEnzdjxyB1at1pWPlSu3vPXy4ftwuvVS7eBTWYT46dIi3Y2L4r2/fIv3tvcb+/dpfuWtX\n+OILz1f6bDb48EN44w39fo4Z49njlQEjmg0Gg8F3MKLZCWYfOcLUAwdY168f/kV1d2dnw7hxetTX\nokXQvLlT+YoID0REEJeby3ddOzFj/XTeWfMON3S/gVeGv0KL+qX4JYjolua339aDp847r+wnVw7y\nbfks2L6Al1e8zKDWg3j7orepm9uepUthyRL45x/o2FE3iF9yiRbMRiSfzoxDh/j48GHWnHNOiYK5\nKI4ehX//heXLdSdHjRpaPF92mRbSs5IP8fnhw6zp169y+DivWgU33QTPPQdPPFGxTulr1mhXqUmT\ndCt3JcCIZoPBYPAdjGguhRXJydy8cyer+vala1EtYgkJMHasbumdPbv0kf4O5NpsDA5dxYHY5QzL\nDWPaJdPo0qxL2Qr5119w550wY4Z2C6lARGDrjuNM/HUaK46/T/Utj3B5gxe56vLajB7tdP3hjGVJ\nQgIPREQQcs455W4NFtF+0X/+Cb//roOuDDtXyHogkhptsvl7QC/nB596gt9/1yHhFiyAiy7yThn2\n7tU1iiuv1JVNTwXWdhIjmg0Gg8F3MKK5BCKysjh/yxYWdu/OqKIGKMXH6+a8yy7T/stl/AAnZiXy\nv7/+xz8xG5BzZjClQzfua9XKtcKGhel4zu++C7fd5loeTiKiw+f+8IMOVJCeDldfDeeNieWHjCfZ\nEb+Nz6/4nOFBwz1aDl9ne0YGo8LCWNKrF4MbNnR7/qmpuj61aLGNHwZu46yU+vyvVkeuu073AFQo\n332nB/z99hsMHlzBB3cgMVFXdHv1gk8+8apwNqLZYDAYfAcjmosh12ZjyObN3N+yJQ+3LsKnuFAw\njx2rY8KWsQVv8Z7FPLjkQW7scSNTR00lNr8a523ZwrLevenboIFrhd65Ey6+GKZO1S16bmbPHli4\nEL79Vnuk3HCDDqc7cOCpumPR7kU8tvQxrux8Je9e8q7TMaPPJLILChi4eTNPt2nDuJaeC+NXyNGs\nPPqFbqbn2vaETW9OixZw443aU6JDBw8ffO5cHVLujz+gd28PH8xJ0tO1H0vfvvDxx16LXWhEs8Fg\nMPgORjQXwwv79rE9M5PfevY8PTpEQgKMGuWSYE7PSeeJP59g5f6VzL5q9imtsXPj4ngnJoYN/fpR\n2zEkgrNEROiu79dfh9tvdy0PO+LitFCeP18PQrvxRrjlFh0GraTTTs1O5Yk/n2BNzBq+ufYbBrYe\nWO6yVCX+FxXFoZwcvu/e3enoI+UlNC2Nsdu3s6nfAPaur8X33+uegsBAPZnezTfr2Qzdym+/wYMP\nasfrbt3cnHk5SUvTwnnAAJg+3SvC2Yhmg8Fg8B2MaC6C4JQUbti5k60DBhDgOPAvK0sL5uHD4c03\ny/Sh3Rq3lZt+vInzAs/jg9Ef0KDWqS3KIsL14eGcXacO75Sn+W/nTl3GuXNdiuWck6O1zuzZOvrF\nVVdp/T1y5OnhzUrjh/AfeOyPxxg/aDzPn/881ZR3fUgrA38nJXHvnj2EDRhQ9MQ4HuTl6Gi2pKez\nuFcvlFLk5+soHAsW6JDJAwbAHXfo4Bblnk1x1SrdHfH777o7ojKSmqormVdcoQcIVjBGNBsMBoPv\nYESzA2n5+fTduJEPOnZkrON0xAUFWgTUq6cFqZOCWUT4dOOnvLLyFT4c/SG39rq12LTxubn02biR\nhd27M7y4yVOcISQErrkGli7VSsgJtm/XEcAWLtS96HffrbMob0SwQ2mHuPWnW6lXsx7zrplHs7o+\nOs2zG0jMy6PPhg183a0bF7prIo8yUOh29HCrVtzv0Kx8/DgsXqwf7ZAQXVkaN07XD8vcCLt1q66w\nLVxYavhFr3P0qA7tMmmSjulcgRjRbDAYDL6D10SzUuor4HLgmIj0KibNdGAMkAWME5EtRaRxqwF+\nNCKCHJuNL7t2PX3jk0/Ctm06PIGTM60dzzvOg0seJOxoGD/c8AOdz+pc6j5LEhJ4PCqKHQMHUs9V\nNw3Q8ZsffhiCg4ud1CEzU/sof/65ngfinnu0WA4Kcv2wRZFvy+eFf17gu/Dv+O767xjSZoh7D+Aj\n3Lt7N/WBxmGgAAAgAElEQVSrV+fDTp28VobwzEyGb9nCxv79i534JC5Otz5/9ZX2X7/nHi2gnXLf\nOHpU+++884725/EFdu2CESP0SVegyDei2WAwGHwHl222iJRrAc4HzgG2F7P9MmCp9f9gYF0x6cRd\nbExLk4DgYEnKzT1944wZIt27iyQnO53fwZSD0v+z/nLLj7dIZm5mmcpyS3i4vLh3b5n2KZLp00V6\n9RLJyDhl9c6dIuPHizRtKnLFFSKLF4vk55f/cKWxaNci8X/bX2Zvme35g1Uy1qSkSKuQEEnNy/N2\nUeTV6Gi5Zvv2UtPZbCLr1oncf79I48YiV18tsnRpCc9KdrbIsGEikya5tbwVwsqVIv7++uWoICz7\nVW576kuLO222wWAwVCSu2my3uGcopYKAxVJES7NS6lNghYh8Z/3eDQwXkaMO6cQdZbGJMGzzZh5o\n1Yp7HKMZhIRoJ881a5wONbAmZg3Xf389Tw19iqeHPl3mwV6xOTn02bCBdf360bG4KbudQUQ3E2Zl\nUfDNtyz5XTFjhnbFuO8+PRt3BU4mCMCu+F2M/XYsYzuP5e2L36Z6tXK0pvsI+TYbAzdvZkJgILcE\nBHi7OGQXFNBjwwY+6dyZS5s2dWqf9HTdK/HZZzpq24MP6hnjT0x9LqJXpKbqeIRejoHsEl99pVvI\n168HV6PYlAHT0mwwGAy+g6s2uyK+hq2BGLvfh4A2njrYnLg4AMa1cJiFrzBkxOzZTgvm73Z8x9Xf\nXs2ssbN4ZtgzLkVHaF2rFhPbtuXJqKgy73sKSpH65kzi1kbzTvO3eeMN3c1+4ICOTFfRghmgm383\nQu8LJexoGFcsvIL0nPSKL0QFM/PwYZr4+XFzJZntpXb16nzQsSPjIyPJtdmc2qdBA13J2rgRvv9e\nB2rp3FkPHFy/Hj25zqZNeqpqXxTMoCuYF1ygfZSMsDMYDAaDG6iIlubFwJsiEmL9Xg5MEJHNDulk\nkt2o9xEjRjBixIgylSM5L49u69fze+/e9LdvXcrL01EoLrrIqZH1IsKbwW8yc+NMFt+ymD4t+pSp\nHI7k2mz03rCBdzt04ArHQYlOEB2tZ9meNw9uveAQ74UMoubCud6bjc2BfFs+j/7+KBsOb+D3W3+n\nZQPPxyv2Bkdzc+m5YQOr+vale3lHVboREeGK7dsZ3rgxE9q2dSmPpCTdOLv6vQ3Mjr+c4LfXMvrR\nDs66/FdOcnLg/PN1ZfmZZ9ya9cqVK1m5cuWJ3//3f/9nWpoNBoPBR/Bq9Awn3DNWisi31m+PuWc8\nGRlJjggzOzsM0nv2WR2+bfHiUlvOCmwFjP9jPGsOrWHJLUto3bCICVFc4O+kJB6OiGDnoEHUcrL1\nbsMGPTHgP//o3vLHH4c2bdCxcu+4A7ZsqTTzW4sIbwS/weebPmfpbUvp7t/d20VyOw/t2UP96tV5\nt8Kn4SudyKwshmzezPaBA2lVq5ZrmaSkIP36seGGt3l+4/Xs3g2PPqrdN846y73lrTBiYvRgxu++\n0y3PHsK4ZxgMBoPvUJndM34D7gRQSg0BUhwFszvYf/w4844e5f8cw0UsW6YdOJ3oas7Jz+HWn29l\nV8IuVo1b5TbBDHBJ06Z0r1ePTw8fLjGdiJ4yedQouO46HUErOhreessSzKA3jhunw2o52SXvaZRS\nvHD+C7w68lVGfj2StTFrvV0ktxKRlcVPCQm84A0/GCfoVLcu97RsyZT9+13LQATuvx81ZgyD3rqe\nf/7RUQ4jI/VU3Y88ov/3OQIDYdYsXclMTvZ2aQwGg8Hgw7gj5NxCYDjQDDgKTAJqAIjIZ1aaGcBo\nIBO429E1w0pTrlaLu3btIqh2bf6vffuTK+Pj9fS6c+eWGn4qIzeDa767hoa1GvLNtd9Q26+2y2Up\njm0ZGVwSFkbk4ME08PM7ZVtBAfz8s54AMC8PJk7UM7sVO2dGXp4OvHvddfD0024va3n4I/IP7lx0\nJwuuXcDFHS72dnHcwo3h4ZxTvz7PV1LRDDp2dJfQUNcGnc6cqeMVrl0LtU999uPitJvzZ59pb4cJ\nE2CIr0UafPJJiI3VTtwemDHQtDQbDAaD7+C1kHPuWihH+KLt6enSPDj41BBgNpuOwTZhQqn7pxxP\nkaFfDpV7f71X8gs8G6/t1vBwmRIdfeJ3bq7InDkiXbqIDBok8uuvIgUFTmYWHa1Da23Y4Imilov/\n9v8n/m/7y4/hP3q7KOVmfWqqtAoJkcyKiOVXTl6NjpZbwsPLttOePSLNmons3l1isowMkY8+EgkK\nErngApHff9evmU9w/LhI794is2Z5JHsqYcg5dEPFbiASmFhMmunW9jDgHGtdILACCAd2AOOL2df9\nF9JgMBgqAFdtttcN+4mClMMAj922Td47ePDUlZ98ItK/v0hOTon7JmQmSP/P+stjvz8mBTZn1arr\nRGVlyVmrV0tsRo58+aXI2WeLDB8usny5iwJk4UKRbt20KKhkbD68WVq820Lmhc3zdlFcxmazyagt\nW+TT2FhvF8Up0vPyJCA4WLakpTm3Q16erq3NmOH0MfLyRBYs0Bq0Tx/9CPpAfUJkxw5dOYiIcHvW\nlU00A9WBKCAI3fO3FejmkKbIGPpAC6Cv9X99YI/jvmJEs8Fg8GFctdk+Gk/qJGtSU9makcHD9lOc\nRUfDyy/D/Pklzvh3LPMYo+aOYlT7UUwfM51qyvOXo61fHXomNKfblIN89x3MmQMrV2rvEZd6jW+6\nCXr0gFdecXNJy885Lc9h+R3Lmbh8InO2zvF2cVxieXIyMTk53OMYwrCSUt/PjxfatePF6Gjndnjz\nTWjUSM846SR+fnDLLXqG7dde064bXbvq6Bu5uS4WvCLo0UPbhXHjtD9U1WYQECUi+0UkD/gWuMoh\nzVjgawARCQUaK6UCRCRORLZa6zOAXYAzc0gaDAZDlcbnRfMr0dG8EhRE7cJpqm02HaN14kT9JS+G\n+Mx4Lpx7IWM7j+Wti95yKQZzWcjL06Kic2ewzW0HY+KYsziH888vZ8ZKwSef6Hh0ISFuKas76dG8\nB//c+Q8v/fsSX2z6wtvFKRMiwv/t38/koCBq+FC84gdbtSI8M5OQ1NSSE27eDNOn6wfThfNTCi6/\nXM/uPmuWHm/bqZN+HHNyXCy8p3nsMa36P/zQ2yXxNEXFx3cc2VxqDH0rMtI5QKjbS2gwGAw+hu8o\ngSJYm5pK1PHj3Gk/M9unn0J2Njz1VLH7JWQlcNG8ixjbeSxTRk7xqGAuKIBvvoHu3XXD97x58N/P\ntbinTQumxcSUnoEz+PtrpTJuHGRmuidPN9K1WVdW3LWCV/97lVmbZ3m7OE6zMiWFY3l53FRJwvo5\nS61q1XihXTumHjhQfKLcXP28vPeeXVgW17ngAvj7bz3ObulSHXFjxgz9KlYqqlXTlYTXX9ezulRd\nnB2h52j8TuynlKoP/Ag8YbU4GwwGwxmNX+lJKi+vHTjAc23bnmwFjI7WbgrBwVC96Cmdk44ncfG8\nixnTcQxTR031mGAWgUWL4KWXoHFjHXlg1KiT258JDKTXhg0817Yt/u6YQeKaa+Cnn/QB33+//Pm5\nmU5ndWL5ncsZ9fUo/Kr5cVffu7xdpFJ59cABXmzXjuoe7oXwBHe1aMGrBw6wKT391Il+CnnrLR2O\n7bbb3HrcwYNhyRI92+CUKdr747nn9FTvtd0fkMY1OnTQkxzdfTf891+xtsLHiUUP6CskEN2SXFKa\nNtY6lFI1gJ+A+SKyqLiDTJ48+cT/rkxIZTAYDBWB44RULuOKI7QnFso4qGRzWpq0CgmR44UjkGw2\nkYsuEnnrrWL3Sc1OlQGfD5Cn/3pabB4c9v/PP3psVZ8+JUcYeHD3bnlh7173HTg+XiQgQGT9evfl\n6WZ2xe+Slu+2lG+2fePtopRIcEqKtF+7VnKdDmVS+fggJkau2b799A07d4qcdZbIgQMeL8OGDSKX\nXy7Spo0em1vKuNyKo6BAhwD54AO3ZEflGwjoB+xFDwSsSekDAYdwciCgAuYC75dyDLdcO4PBYKho\nXLXZXjfuJwpSRgN83fbt8r59xIz587VKzc0tMn1mbqZcMPsCeXjJwx4TzJs2iVx8sUiHDjqiQGl6\na19WljRdvVqSiymzS8ybV+J1qAxsP7pdWrzbQn7e+bO3i1Isl27dKp/7SMSM4sjMz5eA4GDZnp5+\ncmVBgciwYWWKluEOQkNFLr1Uh6ubNUtH4PA6u3fryoNj5B0XqGyiWReJMejIF1HA89a6B4EH7dLM\nsLaHAf2sdecBNktob7GW0UXkX+7rZjAYDN7AVZvtlmm03UFZAuWHZ2Zy4dat7B0yhHrVq0NSkh4Z\n/+uvespcB3Lyc7j6u6vxr+vPnKvnuD1Kxr592itixQo9OP/++0uYlMSBu3btolOdOrzkOJOhq4jA\npZfqcBwTJ7onTw+w6fAmxnwzhgXXLeCisy/ydnFOYUNaGteFhxM1eDA1fWgAYFG8eeAA2zMz+aa7\nNa35xx/DwoXaLcEL5xYSot+V2FjtvnHjjV4pxkmmTIFNm7QvVTnccMzkJgbQY1gyMrQv//HjJ4O0\nKKUDOdWrB3Xr6v990OvLYKgyuGqzfVI037FrF93r1j05O9u992prNH36aWkLbAXc8tMt5Nvy+f6G\n7/Gr5j437sREmDpVTzj4xBN67GH9+mXLY3dmJhds3cq+wYOp7+emsu3bpysP69bpEVmVlNUHVnPt\n99fy282/MTRwqLeLc4Lrd+xgeOPGPO6GAXLeJi0/nw6hoaw55xw6paZCr15aMBeKaC/xzz/w4ota\nWEydCldc4SURkZOjZw197TW49lqXszGiueqTm6tNa2SkXg4c0JW/Q4f05LNJSZCerkVx7dpQp44O\n1KJ7dPX+WVl6rHa1atCsmV4CAvRY3LZtIShIR1jq3BnOOsvbZ2xwiuxs/RDExuolI0Pf7Nxc/RA0\naqQHNrVurcdTlFUkGDzCGSOaY7Kz6bNxI/sGD6ZxjRqwahXcfjuEh0PDhqekFREe+f0RIpIiWHrr\nUmr51XJLWbOzdWSAt96CG27QY4rsA3iUFY+ItHfe0crkjz8qdZPGH5F/MO7XcSy/Yzm9Anp5uzhE\nZWUxdMsW9hf2YlQBXo6OJiEvj5mTJ+uv8uuve7tIgBYSixdr8dyokR40eN55XijI6tU68HR4uC6I\nCxjRXLXIytIRGUNDYcsW2LZNC+U2bXRYxU6d9KvUurVeAgKgaVOtjZzpOTl+HBIS9BIXBzExeomO\n1kFd9uzRrdF9+ujlnHP01PUdOlRqc171EdEPwt9/w4YN+iGJioIWLfTD0bq11iE1a+ru5uxsSEmB\n5GQtrPfu1Tamd289anrIEDj3XP3gGCqUM0Y0P7t3LwUivNexow5+3LcvvPpqka1Er6x4hd8jf2fF\nXStoWKthEbmVDREdUuu55/Qz/9ZbJYaCdpq1qanctmsXkYMHuy9SQynXpjLx7Y5veXbZs6y+ezVB\njYO8WpZHIyJo4ufH1LPP9mo53MnR3Fy6hoQQMX48/qGhuimsElFQAAsW6MA3PXrAG2/oBvEK5b77\ndG+Vi/GbjWj2bdLTdd1pxQo92VR4OPTsqXVNv35auHbvXnERYETgyBEIC9OTCBUK+KwsGDoURozQ\nS9++VTX4SyVCRPfafvONjqeZk6NdIIcO1Q9Hz55Qy8kGOZtNt0Zv3apv6Lp1sH69vpGXXQZXXQXd\nunn2fAzAGSKa0/Lzab9uHZsHDKBd7dowbRosW1Zka+rH6z/mw9APCb4nmOb1yh9nd906+N//9Psy\nbRqMHFnuLE9h2ObNPNWmDde7MybwypVw112wc6cWBJWYj0I/YsaGGQTfHYx/PX+vlCEhN5fO69ez\na9AgAtwRBrCykJPDfe+8Q9uBA3nl0ku9XZpiycnRYdZffx3GjNHuxm3bVtDBExK0Klq2TCukMmJE\ns28hArt2we+/ax20cSMMGKDt+ogRMHCg7lmvbMTGwpo1J8V9XBxcdJF+X0aPhpYtvV3CKsTRo/D5\n59r/snp1uOMOLWp79HBvc//x47rHfOlSHTa2eXMdCvTWW6GVmYjTU7hss10ZPWi/AKOB3UAkMLGI\n7SOAVE6Own6pmHxKHe047eBBuTk8XP+IjdUj3/fsOS3dTzt/klbTWsm+pH2l5lkaBw+K3HqrSKtW\nIrNnlx4Rw1V+OnZMhmza5P6Mb71V5Lnn3J+vB3jxnxdlwOcDJD0nvfTEHuD/oqPlvt27vXJsjzJ1\nqoTffbcEBAefDNFYiUlNFXnpJZGmTUWefVYkKamCDvzppyLnnVd8jMgSoBJGz/D04ozNrkzYbCIb\nN2pz2LmzSGCgyMMPiyxeLJKR4e3SuUZsrMhXX4nceKNIkyYigweLvPGGDgxjcJEdO0TuvVekcWOR\nBx7QIVw9GKL2FPLzRf79Vx+/SROR667TMWwr6vhnEK7a7HK1NCulqqNDGl2EDoq/AbhFRHbZpRkB\nPCUiY0vJS0oqS57NRofQUH7u0YMBDRvqmli7dqf5Z4YcDOGa767hz9v/pF/Lfi6fW1YWvP02fPQR\nPPKIDkThSf/9AhG6hIYyt1s3hrnoV1kkR47ovu7gYPf4kngQEeH+xfcTmx7Lbzf/Ro3qToYgcQPH\nCwoIWreOlX370q2St8qXiYMHdRfixo1cnpbG1c2acb+PtF4cPgyTJ+vAFhMn6hmwne0FdYmCAt0f\n/8QTulWpDJiW5spLZKTuWV+wQLcw33CD9ljr379q+Qfn5ekGy19+0Uvz5nDzzXDTTdC+vbdLVzSH\nsrNZk5bGtowMtmVmcjA7m9SCAtLy8ykQobGfH01q1KBVzZr0rV+fvvXrM6hhQ93T7G4iIrSP2MqV\n2tg89JAeqekt0tP1NMIzZ2rb9OyzuvW5KvWCehGvtDQDQ4E/7X4/BzznkGYEsNiJvEqsFSyIi5Ph\nmzfrHytW6GYCh+aBXfG7JOCdAPkz8k8n6hlFY7OJLFigs7/pJpH9+13Oqsx8VNxkFOXl/fdFLrnE\nJ2qrufm5Mmb+GLln0T0enYDGkU9jY+WKbdsq7HgVxo03ikyeLCIi/yYlSZd166TAB54De8LDRa68\nUqR9ex3/3KPFDw0VadFCJDm5TLthWporFRkZumfwvPNEmjcXGT9e31ofe/RdJj9fZOVKkYceEmnW\nTOT883V89NRUb5dMZGt6ury4d6/0Wb9ezlq9Wq7evl1e2bdPfjx2TDanpcnerCyJz8mRxNxc2ZuV\nJRvT0uSXY8dkcnS0XL19uzQPDpae69fLC3v3ysa0tPIXKC5O5P779YV67TWRdO/0dBaLzSaybJme\nvK1NG/09z8rydql8HldtdnmN5vXAF3a/bwc+ckgzHEhEB89fCnQvJq9iT85ms8mAjRvlt/h4PStC\nr14i339/Spq49Dhp/0F7mb1ltqvXUDZv1kb2nHNEVq92ORuXycjPl2bBwRKZmenejHNzRbp1E/nt\nN/fm6yHSc9Kl/2f9ZfKKyRVyPJvNJt1CQ+XfCvMDqCD+/VekXbsTBtZms0m/DRtkcXy8d8vlIitW\niPTrp2fb9Oj7ed99Ik88UaZdjGiuHISHizz2mO7ZvvxykZ9/rtTzPFUIOTkiv/wictVVIo0aidx9\nt8jatRVbgcjKz5evjxyRIZs2SZs1a2RiVJQEp6RIvguFyLfZZE1Kijy/d68ErV0rAzZulK8OH5as\nsrqe5eWJfPihFstPP12BfmDlYONGkauv1v6iH30kcvy4t0vks3hLNF/nhGhuANS1/h8DRBSTl0ya\nNOnEsmLFihMntyYlRTqsXatbyGbOFBk+/JQ3PjM3UwZ9MUgmrZjk0sWLjxd58EE9A/Xnn+taureY\nGBUlT0ZGuj/jP//UUxVmZ7s/bw/gjkqQs/ydmCg916+v0JZtj5OXJ9Kzp8iPP56yeu6RI3Lx1q1e\nKlT5KSjQk14GBopcf71IVJQHDnLsmP6Q7txZbJIVK1acYq+MaPYeBQUiv/4qMnKk7iR4+WWRmBhv\nl6pycvSoyFtviXTsqNueZs70bMNqdkGBfBQTI61CQuTSrVvl1/h4yXPjwKB8m02WJCTIZWFhEhAc\nLNMOHnROPG/YoGfOHTVK17R8jY0bda2wbVuRr7/2rmjxUbwlmoc4uGc8TxGDAR32iQaaFrG+2JO7\nOTxcT5mdlKT72rZsObEtvyBfrvn2Grnj5zvKLHry80U+/ljE31/k8ccrR0XzwPHj0nT1aknzxDzD\nV14p8uab7s/XQ+w8tlOav9Nclu9d7tHjXB4WJl/4+JTZp/HRR/qD4PBOZBcUSEBwsIT76sgni6ws\n3ZN61lm6kaiM3hSl8957et5vJ22KEc0VT1aWyCefiHTqJDJggMg33+hWVUPpFBSILF+uGy2bNhV5\n8kmRvXvdmL/NJrMPH5bANWvk8rAw2eQON4pS2JaeLldt2yatQ0Lks9jYoluxs7NFXnhB64h583zf\nXyc4WGToUJHevXXDmMFpvCWa/YC9QBBQE9gKdHNIE8DJ0HaDgP3F5FXkiR3KzpYmq1dLSl6efrPv\nv/+U7U/9+ZSMmDNCcvLLZi2Dg0X69tWN1pXNlfW67dtlxqFD7s84MlKrjMOH3Z+3h1gRvUL83/aX\n8GOeaQ2IyMyUZsHBZe/aq8wkJOiaYDH+8a/s2ycPFxF1xhc5ckR7UwQE6Aqw2+qaOTkiXbqILFni\nVHIjmiuO5GRdYQoIEBk7VuS//3xf+3iT/ftFJkzQn4brrhMJCSlfflvT02XYpk0yaONGWZOS4p5C\nloHQ1FQ5b/NmGbhx46k+z2Fhuvftqqu04agq2GzaD6lTJ5HLLhPZtcvbJfIJvCKa9XEZg46gEQU8\nb617EHjQ+v9RYIclqNcAQ4rJp8gTe2nfPnksIkI/CM2a6f4li882fiZdPuoiSVnONxEfOSJy550i\nrVtXwKAiF1mVnCxdQ0M9M2BrwgSRcePcn68HmRc2T4I+CJK49Di35z0+IkKec2cTS2Xg8cd1LK1i\nOJydLY1Xr5akKuTsuXWr7p7v3l3kr7/clOnSpfpD5ETzpRHNnicxUeTFF7W4u/NOHRnM4D7S00Wm\nTxc5+2yRc8/VofjK4klxPD9fnomKEv/gYPksNtarA45tVkt3QHCwjI+IkIyPP9b6Yc6cyvnRdwc5\nOSLTpunzfOIJD3S/VS28JprdtRRlgI/n50vz4GDZnZmpa1Dvvnti27K9yyTgnQCJTHTO/9fe5//Z\nZ0UqoLfIZWw2m/RZv17+Skx0f+apqdrxzxMxoT3IpBWTZPAXgyUr132jhlPz8qTJ6tVysCoNpti5\nU6uKY8dKTHZreLi8e/BgBRWqYrDZ9ICnDh1ErriiyBDuZcfB7hSHEc2eo1AsN22qexX2lT/8vqEE\n8vJEvv1W98T27KndXkrriNuUlibdQ0Pl+h075Fgl8pFJOHpUbv/kE+n83XcSeqbUso4d0/GlAwJE\nvvzSc5NL+DhVUjTPOXJERoeF6XArdoPYCn1dV+1f5dTFCQ7WLj8jR5Y4tqdSMevwYbk8LMwzmX/6\nqciIET5V47bZbHLLj7fITT/cJAU29xiB6TExckNVM6ROirzQ1FQJWrvWpdHrlZ3sbJG339Z1h6ee\nKmeDS2EPVykRR4xodj9paSKvvqrv4333iURHe/RwBgdsNpE//tCtzp066RB+jp1TBTabvL5/v/gH\nB8v8uLjKNZh682Ydp/KJJ+T72FhpHhwsk6Ojq6TNK5KNG7W/88CBeoIWwylUOdFcGB7r92PHtOK1\nogAkZCZIhw87OBVV4dgxHV6nMrtiFEdWfr74BwdLhLvDz4mcjKzwyy/uz9uDHM87LkO/HCov//ty\nufOy2WzSZd06+a8qdWGVMULK4I0bZZGPhp9zhri4k/7On31WjgHmjz2mlxIwotl95OToXsGAAD2h\naUSERw5jcBKbTUevHDFCm5evv9afkPicHBkdFibnbtpU+Xrr5szRld1vvz2x6lB2tozaskVGbdki\nRytRa7hHKSjQtZ0WLXSIsIQEb5eo0lDlRPO61FQ5e+1aKZg168TUtjn5OTJ89nCZuGxiiRejoEB/\nJP39Rf73v8oR0N0VJkRFyf88EX5ORDt+duzoc8PNj2YclaAPgmR+2Pxy5bMsMVF6VaUwc3l5Ij16\n6AEhTjLPx8PPOcumTSfjr//3nwsZxMfrD3AJA2yMaC4/BQW6caN9e5ExY/S4LUPlYsUKkQsuEGlz\nSYo0W7ZGnomMktzK1P2fl6f9eTt1KtLpPd9mk5f27ZPWISGyuio1mJRGcrIe69K8uY6rW5numZeo\ncqL5jp075Z3ISB3E2xI39/56r4xdOLbE7vnNm0UGD9a9Er6uB/ZlZclZq1dLpqciO4wZo0Nr+Rjb\nj24X/7f9Zc3BNS7ncfX27fJpVQoz9/nn+mtWhkpAdkGBNA8Olj2e6M2oZNhsWpAFBorcfLNImd25\n331XO0oXgxHN5WP1at2L3L+/btU0VF6+jD0sjf8Nls7j4qVXL5FFiypJL25ysp759pJLSo0f+3tC\ngjQPDpaPPRGlqjKzZYvIsGF6hqiNG71dGq9SpURzfE6ONF69WhKmTBG57TYREXlvzXvSe2ZvScsu\negRfaqquYDZvXrV83y8PC5MvPRUiLjxcN8dXhgDVZWTJniXS8t2Wsj95f5n3LYyFne6JWNjeID1d\npGVLHbC/jDy3d69nJtOppGRk6MkvmjbV/rJO9ypnZ+uwAsuWFbnZiGbX2LtXhzkLDBSZP7/q2O2q\nSF5BgYyPiJDO69bJrowMsdn0pDK9eokMGSKyyrkhRp4hMlKkc2c9X7qTdj0qK0u6h4bKQ3v2SM6Z\n9OAZlw0RqWKi+c0DB2Tcpk36y7Z/v/wZ+ae0eLdFkQLJZhP54Qftt3zPPaWO1/E5liYkyDkbNnjO\njeCBB/TsED7ItDXTpPfM3pKeU7YprV7Yu1fGVyVHyVdeOVG5LCv7rQpERlWKU+0E+/aJXHONdgVw\nute8z2sAACAASURBVKXshx/0LGJFXCsjmstGWprI88/rQX5Tp56Y6d1QSUnOzZWLtm6V0WFhkuww\nGjA/X2TuXJGgID0OuZjw8J4jJEQ7wM+cWeZdU/Py5Ipt22T45s2SUIVCcDqFvcvGp5+ecbMKVhnR\nnG+zSbs1a2TDhAkiEybI7vjd4v+2v6w+sPq0k967V3sY9Oihu/eqIgU2m5y9dq2s9VSQ+CNHdOXE\nB+M4FbrsXLXwKqcjahTOiLe7qrgkHDp0onLpKmO3bat6MyI6ybJlIt266R7dUucEsNm039ecOadt\nMqLZOWw2PRFbq1Yid9yhH19D5WZfVpZ0Cw2VJyIiSow8kZ0t8sEHWoPdc08F3dvvvtPjDZYudTmL\nfJtNno2Kkk7r1klkVfkulIWwMO3a16ePyMqV3i5NhVFlRPPi+HgZuHq1iL+/JB+Ols4fdZYvN315\nysnm5JycQvfNN08Pg1PVeOfAAbnDk7HyJk/Wjp4+SE5+jpz/1fny4j8vOpX+m7g4ucjXnd3tuece\nkYklD4wtjT8TE6WvJ3szKjm5uXpOgMIpuUscOBwSItKmjYjDx9WI5tLZskWHL+vfX2Tt2jLtavAS\na1NSpGVIiEyPiXF6n+RkbZKaNtXxtT02J8K0afpddJM9/zQ2VgKCg8+sAYKF2Gy6AtK2rcj117t3\nTvVKiqs2u3B6a6+jlBIRYcy2bdw0fz53tAzgcv9ldG3WlQ9Gf3Ai3erV8NBDEBQEM2ZA+/beK3NF\nkZiXR4d164gcPBj/mjXdf4CMDOjcGRYtgkGD3J+/h4nPjGfQl4N448I3uLnnzSWmPW/zZp4ODOQa\nf/8KKp0H2bYNLr4YIiKgUSOXs7GJ0GX9er7u2pVh5cjH14mLg+efh7/+grfegttvB6WKSHjDDXDO\nOfDCCydWKaUQkaJSV1kKbXZpJCfDyy/DDz/A1Klwzz1QvbrrxxURCqSAnPwc8m35+FXzo0b1GtSo\nVgNV5A0zuMKi+Hjuj4hgdpcuXNGsWZn3j4mBF1+EZcvg//5P33c/PzcUzGaDZ5+FP//US2CgGzLV\n/JWUxB27dvFRp07c1Ly52/L1GbKyYNo0+OADuPtufQObNHEpq3xbPsnHk0nJTiElO4WcAv2+FtgK\n8KvmR50adajjV4dGtRvRrG4z6tao6+aTKRlXbXalEs17s7IYvG4dBx9+mEnvjmZr8k6W3rYUv2p+\nJCbCxIn6HfnwQ7j22mI+aFWUu3fvpmvdukxs29YzB/jyS5g7F1at8skLu+3oNi6aexFLb1vKgFYD\nik6TkcFl27axf8gQ/KpVq+ASeoAxY/Qyfny5s3ovJoYtGRnM69bNDQXzbdatg8ceg9q14aOPtD4+\nhagoGDIEdu4E68NqRPPp2GzapDz3HFxzDbz2GjRtWnKe+bZ8opOj2Z2wm6ikKA6mHiQmLYbY9FgS\nsxJJPJ5ISnYKALWq18Kvmh/5tnzybHkU2ApoXLsxTes05ay6ZxHYMJCgxkEENQ6iu393ejbvSfN6\nZ6AQcoEZhw7xxsGD/NarF/0bNChXXps2wdNPQ2Ki1mOXXFKOzHJzYdw4rch//bX0B8oFwjIyuGL7\ndv7Xpg1PuVGQ+xRxcfDKK7oh7dlntUGsU+e0ZDaxEZEYwY5jOwg/Fs6uhF3sT9lPTFoM8ZnxNKzV\nkMa1G9OodiPq+NWherXqVFfVybflk52fzfH846RkpxCfGU81VY1WDVrRrnE72jVqR8emHenWrBvd\n/LvRsWlH/Kq5o8Z1kiohmidERpI/fz5j6xzjvrrLCL0vlCa1mzJ/vr5vN96oWyoaNvR2aSue0LQ0\nbtm5k6jBg6nmCVGbnw99+sCbb8KVV7o//wpg0e5FPP7H46y/bz0tG7Q8bfvDERG0qFmTSUFBFV84\nd7N8ue5y2bkT3ND7kJSXR4fQUCIHDaKZJ3ozfIyCAvjqK3jpJbjuOm13Tvk+P/kk5OXBxx8DRjQ7\nsn07PPIIZGfDJ5/AwIGnp8nOz2bLkS2sj13P1qNb2XJkC3sS99Cyfku6NOtCp6adaNeoHYGNAmnd\noDXN6jajaZ2mNK7dmBrVa5yWX4GtgJTsFJKOJxGfFU9MagwHUg8QnRxNeHw4249tp0a1GgxqPYhh\ngcMYFjiMIW2GUNuvtrsvjc9iE+G5ffv4LSGBP3r3pn0RQskVRE7qry5dtHju2rWMmWRk6NayevVg\nwYIiRZy7OJidzZht27i4SROmdexIdR9sSHILO3fqbqLQUHj5ZXLvvI11xzbzz75/WHtoLetj19Ok\nThP6BPShu393uvt3p33j9gQ2CqRVg1ZOC10RITMvk8PphzmQcoD9KfuJTIpkV8IudsbvJC4jjl7N\ne9GvZT8GthrIsMBhdD6rc7l6lqqEaPZfvpy/35jCxZeEs+ru/6iR0p2HH4akJPjss6IN75mCiNB/\n0yZeb9+e0Wed5ZmDLFkCEybobn+39KNVPFP/m8qSiCWsHLfylI9hen4+7datY8fAgbSqVcuLJXQD\nNhv076/dA264wW3Zjtu1ix716vGsp3ozfJCkJN3g8sMP8OqrcO+9lmtBYqL+6oeEQOfORjRbZGTo\nrvivv4YpU+D++0+6YqRmpxJ8MJiV+1ey+uBqth/bTtdmXRnUahD9Wvajb4u+9Gzekzo1PCOGRIRD\naYcIjQ1lbcxaQmJC2Bm/k3PbnsulHS5lbJexnN3kbI8c2xfItdm4e/duDmRn81uvXjStcXrFpLzk\n5Gi3yjfe0O5PkyY52fufkACXXw69e8PMmRXyfUrOy+PqHTtoXrMm87p2pXZ5fIp8mGOZxwj56UPa\nvfs5zWISWXBZIOm33cDQjiMY1HpQhfTepOWkERYXxqYjm1gfu541MWvIyM3g3LbnMjJoJCODRtIr\noBfVlPM9yF4TzUqp0cAHQHXgSxF5q4g004ExQBYwTkS2FJFGLvzoIxrvnMytj8wm/Jcr+fBDrQvG\nj/dZDedWvjh8mCWJifzaq5dnDiACI0fCbbfpr50PIiLc/NPN1Kpei6+v/vpETXTm/7N33nFVlu0D\n/z6AuAVUcO+ZRpoDt+E2zbJhqVmpaePNsnobppWWv7JcmZWZM0tTy17NssyRFiiKOEBQBEEQWbL3\nOJxz/f540MxAGWeA3t/P53w4437u+3rOOVzneq77GtHR7E9JYdudd9pYQjPw7be6h9PHx6yhNL7p\n6Uw8c4YQS+1mVGJOndJ3KHNz9R/93r3RA599feHHHyuk0Vwe3VzCY/9hNO/YoetqT09YtAhc6ufj\nE+XD3vC97A3fy5mEM/Rq0ot7WtzDwKb96FGtNTXTsiE1VY+lzMnRvff29vqtalVwdtatqvr19bh9\nM38vU3NT2R++n93nd7MzZCdNajfh4TseZoL7hNvKgE4vKOChwEDqODiw6Y47qG5hA/HyZf1idPt2\n3XB+5pkb/MZHRekxHVdifKyom3KNRp4MDuZyfj477rwTZwtcSFREUnNT+T7oezad3oR/nD8j2o7g\nwY4Pcm9SXZwWLdcV4ksv6R+cs7NNZIxOj8brohcHLhzgQMQBUnJTGNZ6GCPajGBE2xE0rNXwhsfb\nxGjWNM0eOAcMBaKBY8AEETl7zZhRwAwRGaVpWi/gUxHpXcRc8n+PjyBi+EC8F8ymXTv9x0k5vf4m\ny2ikmY8Pp3r0oHk1C20p+vrqyikkRN8Gq4RkG7IZsH4AE+6cwGt9X0NE6OLnxydt2zKkjEkNFYbc\nXH1/c+NGGDDArFOLCD2OH+fD1q0ZYYFYwcqOCGzapOdWjBwJC97NwW1AB9iyBa1fvwplNJdHN5fk\n2MLjRUSIjNR/P8+dg/nLYkhz+5VfQ3/l2Nn9DM5rzH2mtnhkOtE0MR/7i1FIZCSSkIDm5IRWv75u\nFNeooW+3OzrqsTFGo/5dT03VMwkTEvRFW7bUs7/d3fVwsi5doF07sxhSRpMR74ve/HDmB7YGbaWz\na2cmd53MuE7jqOlYOXVhSYjJy2NUQAD9nJxY3q6dVUMRAgL0SKeEBPjkExg69LoBwcEwYoQ+6JVX\nrCbXtZhEeOX8ef5ITeVXd3eaWeq318aICAciDrDSbyW/h/3O8DbDmeQ+iRFtR/w7hOnUKT3GZtcu\neOIJeP75MsTblE1GE/pnIoC9pl39vkamRrInbA97wvewL3wfbVzaMLrdaO5rfx/dG3f/lxfaVkZz\nH2CuiIwsfDyr8MQ+umbMSuCAiGwtfBwM3CMi8dfNJcOm3MuZvbtY/qnGgw9Wynw0i/NiaCjODg7M\nt2TZkPHjoXNnPZapkhKVFkWvNb1Ye/9aarv24elz5wj28Kj82fWLFukhATt2WGT6NYW7GTsstZtx\nC5CerocgfPMNfDfyG4aeX4ndEZ+KZjSXVTd7Aq1udmzh8/LxxyY+XHece0d+S13jHzik2NHYsTVG\nO2fiazkR27IlCW5uXHZyIqVaNbLs7ckCjNfM46hpODk44OTgQD0HB5pXq0bzqlVpXb067jVr4l6z\nJs4ODroBHRmpJ2KePq1bXMeP68Z1//4wcKCeGNu+fbl/PPIK8tgVuov1p9ZzOOowT3V5iv/0/A9t\n67Yt17wVjTNZWYwKCODZxo2Z1by5TfSjiO5xfu01PfpiyRJo0wbw89Pzaz76CJ56yupy/VNGYUlU\nFMujo9nl7o57rVo2lcecZORlsP7UelYcW4GDnQPP93ieie4TcaleAgdTdLS+67lunW40P/MMjB2r\nXwSXgcyCAoKyswnIzOR8Tg4X8/KIzM3lcn4+aUYjqQUFGEWwQzd6jSJUtbOjpp0dTg4OuFapgpuj\nI25VHJDceGIT/QmK2k9ORhgPtOrHAx3uZ0irIVSvUt1mRvMjwAgRmV74eBLQS0RevGbMz8ACETlc\n+Hgf8KaIHL9uLnnuxUw+/r+at2WiX0kJyspiqL8/F3v3poqlKkCEh+ul566pDlAZOXTxEA9ufZBe\nQ3cyuH4TXqnsmdDJybqX2cvLYlf1WUYjzQt3M25Vj4q5CAqCmTOMzIwdyv3nDlY0o7nMuhloCYy8\n0bGFz8uQD94mtV5jwho3xVilCm3z82lZrRrNXVxo7uZG42rVcHN0xLVKFVwcHKhlb09Ne3uq2Nnp\ndU/RY2nTjEbSCgpIMBiIys0lMi+PsJwcTmdlcTozE1dHR/o7OTHQyYl7nJ1pV7363wbexYvg7Q0H\nDsCvv+o/2Pfdp2eO9+5dbgM6IjWClX4rWXtyLX2b9WVWv1n0adanXHNWBP5KTWVcUBCL27ThiYY3\n3sq2Brm5eqWzxYvhw2EHmLbvMezWroH777e1aFfZHB/PzPPn2dqpE4Mq+a5lXGYcy48uZ9XxVQxq\nNYgXPV5kQPMBZbtwys+Hn3+G1av1pMHRo2HCBH3r4AY5RJdyczmQmopXWhpeaWlE5uZyR40a3FWr\nFu2qV6dFtWq0qFqVBo6OVy+sq15j94gIOSYTWYUGdYLBQILBQExeHlF5eVzMzSUiN5eQ7ExSCgpw\nzE/ElBZG18wCjkx9t0w6u7yRwiW1uK8XrMjjGtRdxNKl+n1PT088PT3LLNitSueaNWlfvTo7EhMZ\nZymDtnVrPUvj/ff1GJlKSr/m/Zg96GP+m5zM8na3gIfoww/1Ug4W3AaraW/P4w0asCo21rK7GZWc\ngwcPcvDgQbr3zeWJ0+jBDBWLsurmUmHvf4YuDdMZdCmGUYMHM2jkyBIfq2kaGlDN3p5q9vY0cHSk\nPfyr5rhJhJDs7Ks/rO9FRFDd3p776tXj/nr1GNisGfYTJ8LEibrb0t9fL0c2ZYr+Yz5hgn6/bdl0\nQEvnlnw09CPm3jOX9afWM/F/E2lWpxlvD3ybYa2HVcrdqy3x8bx0/jyb7riDYRUkFKtaNb084TNu\nO3B44RnG1fieMcmePGmCilIhdEKDBjR0dOSxM2dY2qYNkyrAxUZpiU6P5iPvj9h0ehMT7pzA0WlH\naVO3TfkmdXTUf5seflgPWP/hB32HYOJEPU9q9GgYPBhp1YrjmZn8lJTEL0lJROXmMsjFhQFOTjzf\nuDHuNWuWqhyspmnUsLenhr09ro6OtLt+QEqKrg9OnOC3PXvYERJCIhqRzqWvO351zXJ6mnsD867Z\nxnsLMF2bNFK4BXhQRLYUPi42PKOiVPKo6GyJj2dVbCx/dO1quUUSE3XjzMdHjxmspHx88SIbQg7Q\nPG4zuybuwt6ukmZAR0ToFTOCgsDCivpMVhZDLL2bcQtgNBkZvvURvBs+T/7gERXN01xm3YwennHD\nYwuft4nOFhH8MzP5JSmJ7YmJxObn85ibG4+7udG9du2/jVgRPfZy40Y9ebZLl7+3j8uR0FVgKuD7\noO+Z/9d86teoz/ue7zOo1SAznZ1lEREWXLzIVzEx/Ozuzl0VLczg66/1DkO//IKvsTszZ+rVUD/9\nFPr2tbVwfxOUlcXogACmNmrEOy1aVIoLp/jMeD70+pBvA75l6t1Teb3v6zSo1cCyiyYmwu7dhHp7\ns6lKFb7r2xepXp2Hs7MZU68evTt3xr5Jk/LtBplMEB8PYWF6YkVwsL5LHhCgh3S5u0P37kQ36MY7\n27sRRCe++KoKPXvaJjzDAd3HMgSIAXy5cbJJb2BZcYmAymguGfkmE819fDjYtSsdLZmst2ABnDih\nXzVWQkwitD16lI0dO/DuznF0bdiVxcMX21qssjFpku4pmzfPKst5njzJC02aWG434xbgv7//l59z\nquPRfgKbOt9Z0YzmMuvmkhxbeHyF0NnnsrPZHB/Pt/Hx1HFw4JlGjZjYoAFO15ZjyM3VA2e//JKr\nmYvTp5er6L/RZGRz4Gbe+/M9Wji1YNGwRdzd6PpOOBWHfJOJ50NCOJmZyS/u7hWv9ObSpXp8xp49\nV3fTTCa9JPOsWXDPPbrzsqJE2cXl5TEmMJBONWqwukMHHCuogyEjL4MlPkv4zPcznrjrCWb1n3XT\nyhLmIM9kYntCAqtiYwnKymKCmxuPG4308PVFO3ny77yErCy94kOLFtCggZ4Y7Oysh1ldqaRjNOoV\ndrILq+0kJOi3mBi9uoqTk54g3LGjHsJ4xx16gHzLlmTl2DF/Pqxdq1dqef55fUpblpy7l79LE60V\nkQWapj0LICJfFY75HBgJZAFTROREEfNUCAVcWZgdHk620cgyS3qBs7P1L+APPxTW2Kpc/JaUxNsX\nLuDXvTvJOcn0WtOLd+95lye7PGlr0UrH8eN6jGZICJSzO1dJ2Xr5Ml/FxFh2N6MSs/7kej7wXoD0\n+JqNnTrT19m5QhnNUD7dXNSxRcxfoXS2SYQ/UlL4KjaWfSkpPO7mxsymTWl3fVKSn5+ebbZnDzz7\nLLz6ql7SrowYjAbWnFjD+3+9z7DWw/hg8Ac0c6ogll0hifn5PBwUhHNhSblaFamGq4heW3bHDr2H\nfREls7Ky9AqPX3yhX++8/nqZc83MSpbRyKSzZ0k0GPixc2fcKlBjKKPJyLqT63j34LsMaTWE+YPm\n08rF8iF38fn5fBkdzcqYGDrXrMmzjRvzQP36/4hF/gfp6XpeQmSkbginpOi3nBzdWDYVxufUrKl/\n6HXqgKurfmvUSP++FPFlENGjtGbO1AtNLV78z03aMpcJFZEKcdNFUZSUC9nZUtfLS7IKCiy70Lp1\nIv37i5hMll3HAowJCJDV0dFXHwfGB4rrQlfxifKxoVSlxGQSGTRI5MsvrbpsntEoDQ8dkrOZmVZd\ntzLgHektrgtdZU34Seni6ysmk0kK9ZfN9ag1bxVZZ0fn5srssDCp7+0tYwICxCsl5d+DLlwQee45\nkbp1RV57TeTy5XKtmZ6bLu/88Y7U+7ievH/wfcnOzy7XfObidEaGtPLxkbfCwsRY0fS4wSAybZpI\nz54iCQk3HR4RIfLYYyLNmols3ChiNFpBxptgNJlkTliYtDh8WE5lZNhaHBER+SviL7l75d3Sf11/\nOR5z3Cprns3MlKlnz4qzl5c8ExwsZ2z42xEWJjJqlEjHjiL79xc9pqw62+aK96ogFVgBV1RG+/vL\n2pgYyy5SUCDi7i6yY4dl1zEzkTk54uLlJZnXXVTsDN4pjZc0lqi0KBtJVkp27RLp0EH/cbEys8PC\nZGZIiNXXrchEpERIw8UN5deQX2Xs6dPyVeFFmTKaKyZZBQXy5aVL0trHR/qfOCG/JCaK6XrD8eJF\nkf/8Rzee584VSUsr15oRKRHyyPePSMtlLWX72e3/Xs+K/BAfL67e3rIxLs5mMhRLdrbI2LEiQ4eK\npKeX6lBvb5EePUR69RI5dMhC8pWSLfHxUt/bW7bGx9tMhriMOHn8x8el2dJmsvn0Zqt8946lpclD\np0+Lq7e3vHfhgiTk5Vl8zeLIzhaZN0+kXj2Rjz4SuZEoymi+DfklMVF6+PlZfqHfftMNt/x8y69l\nJuaEhcmMYgy+BV4LpPtX3SUrP8vKUpWSggKRzp1tdsESkZMjdYu48LhdycjLkLu+vEuWHl4qUYUX\nZRmFFzPKaK7YGIxG2RwXJ3f5+srdx47J9suX/+11DQ8XeeIJETc3kU8/vfEvbgnYH75fOn7eUcZ8\nN0YiUiLKNVdpyTca5dXQUGnp4yN+pTRIrUJSkki/fiITJpT5fTYaRTZsEGnSRPc+X7hgXhHLwon0\ndGnl4yMvh4ZKvhXd4AXGAvnC9wupv7C+vLHnDcnMs7yX93Bqqoz095emhw/Lsqgom/9O7Nwp0rq1\nyMMPi0RG3ny8MppvQwpMJmlx+LAcK6dn5KaYTCJDhoisWGHZdczEldCCoGK2h0wmkzz+4+Myftt4\nm3qBbsqaNSIDBtg0NGZMQICssfRuRiXAaDLKg1selKk7porJZJJ3wsPlhXPnrr6ujObKgdFkkh0J\nCdLt2DG5y9dXthVlPAcEiIwcKdKuncj//leu/79cQ67835//J/U+ricLvReKwWj5HaOonBwZcOKE\n3OvvL0kV0dFx8aJIp04ir7xilviKzEyR997TNwreeEOkqEgca5Kcny+j/f2l3/HjEpWTY/H1AuMD\npfea3tJvbT85HX/a4ut5paTI0FOnpMXhw7IyOlpybRwjExKih2J06CCyZ0/Jj1NG823KhxERMuXs\nWcsvdOKESMOGpd5GswWb4+Jk0MmTNxyTnZ8tPVf1lPl/zreSVKUkM1N3oRw5YlMxfktMlG7HjlXs\niwsrMHvfbBmwboDkGnKLvChTRnPlwmQyyc8JCdLDz0/u9PWV7+Pj/208//67Hprm6Sni71+u9c4n\nnZchG4ZI96+6y6nYU+Wa60Zsu3xZ3Ly95YOIiIoXvywicvKkSNOmIosWmX3q6GiRp58220ZBuTCa\nTPJhRIS4eXvLDxYK18g15Mq7f7wr9RfWl5XHVorRZFnj9Yqx3NLHR1ZHR1vVk14U6ekib76ph2Is\nWlT6z1sZzbcpl/PyxOmvvyTRGh6FSZNE3nnH8uuUk/4nTpRIUUWnR0vTpU3lxzM/WkGqUvLee/qe\no40xmkzSxsdHfFJTbS2Kzdjov1FaLWsllzP1RLEt8fHied1FmTKaKycmk0l2JSZKTz8/6Xz0qGyJ\nj5eCa41Ng0FPwnVz05MGS5CsdqO11p5YK64LXWXO/jmSa8g1wxnopBkMMvXsWWnj4yNHLb3zWFZ2\n7xZxdRXZutWiy1zZKGjdWmTLFtsmCx5JS5O2R47I5LNnJc2MeSnHoo9J5y86ywObH5BLaZfMNm9R\n/JmSIoNPnpRWPj6yJibG5say0Siyfr1I48YiTz4pUtaNUGU038ZMOnNGFpUkiKe8REbqe2CXLPtP\nWh78MzKkyaFDJf7H9ov2k/oL68uJmBMWlqwUxMTo73N4uK0lERGRxRcvyqQzZ2wthk3wifIR14Wu\n/9j2HFDERZkymis3JpNJfk1MlF5+fnLH0aPybWysGK7VIUlJIi++qBt9K1bo+QZlJCY9RsZ8N0bc\nV7ibRe/8lJAgTQ8flmnBwZJug4ThErFqlUiDBnoGn5XYv1+ke3c9YXDvXqst+y8yDAaZFhwszQ4f\nlp3luOgS0b3Ls/fNFrdFbrIpYJPFdgBNJpPsS04Wz5MnpY2Pj6yrAMayiIiXl/6Z9ukjcvRo+eZS\nRvNtjE9qqrT28bHOdtysWSKTJ1t+nTLybHCwvFfKjJDvA7+X5p80l9iMWMsIVVqmT9dLYFUQkvLz\nxdnLS+Jtud9pAyJSIqTR4kbyy7lfrj4XUMxFmTKabw1MJpPsSUoSz0LP2opLlyT7WgPZ31/PM7j7\nbpHDh8u1zjenvhHXha4y78A8yS8o/U5hVE6OPBoYKG2PHJE/kpPLLItFMRhEZs4Uad9e5JocAGth\nNIps3izStq2eluPra3URrrI/OVnaHjki4wIDJTq39LsMp2JPifsKd7l/8/0W+626EvPv4ecnHY8e\nla+vv3i0EaGhIg89JNK8ucimTeZJ81FG822MyWSSbseOya+JiZZfLC1N9xicqECe2UJSDQZx9vKS\nmDIopHkH5kmv1b1sX1v19Gl9K9jW2SzX8fTZs/JhRIStxbAaablpcueKO+UTn0/+8XxxF2XKaL71\nOJSaKvcFBIibt7fMDQ+XuCsXjSaT/svduLHIlCki5YhZvZR2SUZ8O0J6rOohZxNKlpuSbjDInLAw\nqevlJXPCwv5p1FckUlNFRozQS8rZ2KjPzxdZuVL/yMaO1UM4bEF2QYG8VfjZvRsefrX6zo0oMBbI\nAq8FUn9hffn65NcW8S5fKc3Y/sgR6eHnJz8WlSBrAy5fFnnpJT1uecECvaScuVBG823OmpgYGV3O\nZJUSs2KF3nCjAvxTXcvyqCh5NDCwTMeaTCaZsG2C7StqjBypZ7FUME6kp0vzw4f/Ge95i1JgLJBR\nm0bJsz8/+4/vwo0uypTRfOtyNjNTng0OFmcvL5l05ox4paTo34u0NJFXX9VDNj7/vMwhGyaTSVb4\nrpB6H9eT5UeWF6t/0gwGWRgZKQ0PHZInzpyRi1aozFBmgoL0cgYzZtikxnxxZGeLLF2q+30esls1\n5gAAIABJREFUe0zEVlFnETk58nhQkDQ6dEg+vUG5tvDkcOm/rr94fu1pkbKFIVlZ8tr58+Lq7S33\nBwTIn1e+2zYmM1Nk/nzdWH7xxXJdlxaLMppvc7IKCqSel5eEm/NSrDgMBpE77tALI1YQjCaTdDhy\nRA6Ww0ObnZ8tvVb3knkH5plRslLw2296masKGgbR5/hx2VHOmLzKwMzfZsrQb4b+a8v80xtclCmj\n+dYnIS9PFl+8KO2PHJHOR4/KoshI3XANDNQrbHTtWq5OGyGJIeKx2kNGbhz5j+33iJwceSssTOp5\necn4oKAK03WuWLZt0y8k1q+3tSTFkpGhey5dXUXGj9dtfFvgl54uDxc2BnknPFxiCy/Ir4Tv1F9Y\nXxYdWmTWyhhpBoNsiI2VISdPiqu3t7x+/ryEZlWMngU5OSLLlumFusaPFzl/3nJrKaNZIf8NDZX/\nhoZaZ7ErneoqSB3Q35OSrrY0Lg+xGbHS4pMW8l3Ad2aSrIRUwAuR69kYFydDblLKr7Lz+dHP5Y7P\n75CUnH9efBlNJml35EjR7ZhFGc23EyaTSQ4kJ8u04GCp6+Ul/U+ckCWRkRL0/fdiKmdKf35Bvrzz\nxzvi9kkbmXniVxl08qTU8/KSGSEhEmYNh0h5yMvTczFatBCxRtMtM5Ce/rfx/MgjekU8WxCSlXV1\nN2PoST/p/dMc6biiq9nKE8bn5ck3sbHySGCg1PnrLxkTECCb4+JsXmP5Crm5evhMs2Yi999f7gqP\nJaKsOlvTj7U9mqZJRZGlsnIhJ4eex48T2acPNe3tLbuYCIwcCaNGwcyZll2rBNwXEMBDrq5MbdSo\n3HOdjj/NkG+GsGP8Dvo262sG6UrAF1/Ajh2wZw9omnXWLCX5JhMtjxxhb5cudK5Z09bimJ1fQ39l\n2s5pHJp6iFYurf75WlIS71y4gF/37mhFfD6apiEiFfODsxBKZ+v/E78nJ/NLUhK7k5MxmUzcExFB\nz99/p0ffvnR+8kmca9S44RwmESJycwnKyuJoejp7U1IIzMxAUk/S3zGbrYNn4lK1lpXOqIyEh8OE\nCeDqCl9/DfXr21qiUpGZCatWwZIl0LUrzJoF/ftbXxXvufAXj/+5ghrNxpJcpTHda9dmqIsL/Zyc\n6FyzJq5VqhSpf66lwGQiPDeXYxkZHEtP51B6OqHZ2Qx2cWFU3bo85OpK3SpVrHRGNyYvD9auhY8+\ngs6dYd486NXLOmuXVWeX2WjWNK0usBVoAUQAj4pIahHjIoB0wAgYRMSjmPluewVsDsaePs3IunV5\nrkkTyy8WFASenhAcDPXqWX69YgjNzqbfyZNE9u5NdTNdLOw+v5spP03Be4o3beq2McucxZKSAh07\nwr594O5u2bXKyfsREUTn5fFVhw62FsWs+Mf5M+zbYeycsJPeTXv/6/WR/v5MaNCApxo2LPJ4ZTQr\nRITg7GwOp6fjFxWFX1gYwS4u2FWtSvNatahfpQpVNQ1HOzuMIiQXFJBsMBCdl0fdKlXoXLMm3WrV\nYpiLC32cnMjNz+D5Xc9zKu4Umx/eTJeGXWx9iv9GBLZs0R0ns2frfyvoRX9JyM2FDRtg8WL9J+2N\nN+CBB8DSPqgCUwHz/5zPqhOrWD1mNfe1v4/MggK80tLYm5KCb3o6Z7KzsQNaV6+Oi4MDzg4OVLOz\nI89kItdkIrWggMi8PGLy8mjk6EiP2rXxqFOH3oU3Rzs7y55EKcjI0C9SPvkEunSBd9+1nrF8BVsY\nzQuBRBFZqGnam4CLiMwqYtwFoLuIJN9kPqWAzcCBlBReCA0lqGfPm16RmoUZM/S/n39u+bWK4eXQ\nUKrb27OgdWuzzvvlsS/59Oin+Dztg0t1F7PO/Q9eeQVycmDlSsutYSbi8/Pp6OtLWK9eFcZbUV4u\npV+i79q+LB6+mEc7P/qv14OzsvA8dYrIPn2oWswPjzKaFUUhv/1G6ttvE9mxIymvvEJeq1bkiWAP\n1K1ShboODjSuWpU6Dg5FHy/CxoCNvLrnVd4Z+A4verxoHb1eEhIS4D//gcBA2LgRune3tURmw2iE\nn36ChQvh8mV46SWYOhXq1DH/WhGpETz+v8epWaUmG8ZuoFHtondLRYTLBgMXcnJILSggtaCAXJOJ\nanZ2VLOzo46DAy2qVaNZ1aoVykC+lrg4fVN15UoYPFj36N99t21kKbPOLktMR6GiDAYaFN5vCAQX\nM+4CUK8E85UjOkVxBZPJJHf6+sqepCTrLJiYqAeElbFqRXlJNxjExctLIi2USf7q7lflnvX3mLV7\n1z84e1ZPEbZQq1VL8OSZM/KxNZrpWIHUnFRxX+EuC70XFjvmP+fOyTs3aTSDimlWFEd+vl5do0ED\nkSeeECllHXkRkdCkUOm5qqeM3jT6amdKm2Eyifzwg56t9dprevbWLYyPj15pw8VF5IUXzPtTtzVw\nq7gudDV7sl9F49Qpvb2Ds7PI88+LhITYWqKy6+zyXI40EJH4wvvxQIPi7HJgn6ZpfpqmTS/HeooS\noGkaM5s0YXl0tHUWrFcP3n5b95bawOu0IS6OIS4uNK9WzSLzLxy2kHo16jF151RMYjLv5CL6duac\nOeDmZt65LchLTZvyRXQ0BSYzvx9WxmA08MgPjzCg+QBe6/takWNSDQY2X77Mc40bW1k6xS1DlSrw\nwgsQEgKtWuke2ZdfhtjYEk/Rtm5bvKd64+7mTtevurIvfJ8FBb4BoaFw770wdy78+CMsWgQW0r0V\nhd699QiUgACoWxeGDtWjEjdv1sM5ykJWfhbTd05nzh9z+PXxX3mt72vYaRXTO1xW8vP1923AABg9\nGtq3h/PnYcUKaNfO1tKVnRuGZ2iathfdi3w9c4ANIuJyzdhkEalbxByNRCRW0zRXYC/wooh4FTFO\n5s6de/Wxp6cnnp6epTkXRSHZRiMtjhzB5+67aXuTJBSzYDDo2RMffABjx1p+vUJMInT09WVthw4M\ncHa22Do5hhyGfDMEz5aefDjkQ/NNvGOHHgfo76//sFYi+p84wSvNmvGwq6utRSkTIsKUn6aQnJPM\n/x77Hw52RW+PL42K4nhGBps6dfrH8wcPHuTgwYNXH7/33nsqPENRMuLi4OOP9eDZp56CN9+EYmLl\ni2Jf+D6e2vEUk9wnMX/wfBztHS0obCHp6bqB/OWXurwvv1zpdJa5yM/XVfeaNXDihJ7/OHkydOtW\nsnBu/zh/xv84Ho8mHnx+7+fUrlrb4jJbk+BgPbnv22/15L7//EePCy8mAslm2Co8o2Hh/UYUE55x\n3TFzgf8W85pZXe+3O7PCwmSGNfdA9u8XadnSvC17bsKOhATp6ednlWLsCVkJ0m55O/ny2JfmmTA7\nW3+/9u41z3xW5vv4eOl7/LitxSgzs/fNFo/VHpKZl1nsGIPRKM0PHxbftLSbzocKz1CUlpgYvd2Z\ns7PI9OkiwcElPvRy5mUZvWm09FzVU0KTLFhmNDdX5JNP9C6lkyaJXLxoubUqIRcuiLz7rkirVnrF\n0A8+ECkukstkMskyn2VSf2F92ei/0apyWpqEBJEvvhDp00eP2nnzTZt0TS8VZdXZ5VGYC4E3C+/P\nAj4qYkwNoHbh/ZrAIWB4MfNZ8O25/YjOzRUXLy9JsmYd5XHjRObOtdpyA06ckC1WjAU+n3ReGi1u\nJNvPbi//ZO+/L/LQQ+Wfx0YYjEZp5eMjh1NTbS1Kqfns6GfSbnm7m8aGbo6Lk4ElbBevjGZFmbl8\nWdebrq4iY8aI/PpriboLmkwmWX5kuWXaK6eliSxZohfOHTPGdn2nKwkmk4i3t8hzz+kfY48eIh9/\n/HfsbnxmvNy78V7xWO0h55Ms2LHDiiQl6f1rRo8WcXLSm5H8/HOFad1wU2xhNNcF9gEhwB7AufD5\nxsCuwvutgVOFt0DgrRvMZ9l36DZk8tmz8kFEhPUWjIwUqVtXJCzM4ksdTUuTFocPi8HKxdn9ov3E\ndaGreEV6lX2SiAj9fSpDQlBFYnlUlDx0+rStxSgV24K2SeMljSU8+caJfSaTSbodOyY7S9gBURnN\ninKTlSWyZo1ucbVsqV9Yl6Almn+cv3T6opNM2DbhX015Sk1QkMh//6snJ48fX2malFQkDAaRfftE\nnnlGpFEjkaaev0nNdxvJxLWzJDO7kliUxRAaqnfsGzpUpHZtkQcfFNm4UW8SU9koq85WzU1uYQIz\nMxkWEMCFXr2oZulCk1f48EM4elSv12NBHgsKok+dOrzcrJlF1ymKPWF7eGL7E/zx5B90dutc+gnu\nvx88PPQEykpMZkEBrY4etV7sfDn548IfjN82nt2TdtOtUbcbjj2YksJzISGc8fDArgSBiqrknMKs\nHD8O69fDtm3QvDmMG6cn4HXuXGTgbLYhm9f2vMavob+y8aGN9G/ev+RrRUbCzp16EGp0NEyaBM8/\nDy1bmu98bkNyC3KZte8tNvtvY0TWNwTvHsS5c3pi3JAh+t+uXSterO+1JCTAwYPwxx+wf7/eBGbU\nKP02YgRU5h5XVq/TbG6UArYM9wYE8IirK0+boVNeicjLg7vu0pNG7r/fIktE5OTQ4/hxLvTuTW0b\naZxNAZt4a/9beE3xooVzi5If+NNPeiKNvz9UrWo5Aa3E2+HhJBcUsKJ9e1uLckP8YvwYtWkU34/7\nHs+Wnjcdf19AAA/Ur8/0ElbNUEazwiIUFOhWy48/wu7detL10KF6J4iePXVd6/h3IuAvIb8w/efp\nTLt7Gu/e8y5V7K9L1hOBiAjw84MjR/Q5L1/WDfKJE/W5K7IVV0k4HX+aSdsn0a5uO1aNWUXd6nqN\nhKQkvYfVn3/CX39BVJT+MV65desGLVrYpj9MdrZecvvECfDx0W/x8TBwoF5TedAgvRFJRSkTXl6U\n0awokv0pKbwYGkpgz54l8piZhQMH9HTioCCoZf72ry+HhuJoZ8fCNhbu1HcTlh9dzue+n+M1xYsG\ntYqruHgNmZm6p2j9el0L3QLE5eVxx7FjhHh44OpohSz+MnAu8RyeGzxZOXolD3R84Kbjz2RlMfjU\nKSJ69y7xDo0ymhUWR0Qv+fbHH3DsmH4LCYHGjaF1a90j7eREZlWNHy7sIj83k0dajqaeqSpcugQX\nL8KFC1C9OvTooVtpw4frJfCstRN5i2MSE58e+ZQPvT9k4dCFTO46+YbNaBITwddXvx07BqdO6d3y\nOneGTp2gTRto21Z3+jdpAg0alO+aJjdXN9SjovRrp5AQOHdOr3gRGak3pu3aVb8m69tXl+FW/Woo\no1lRJCJC9+PHeb9lS+6rX996Cz/5pF57ePFis06bmJ9Pe19fAnr0oGkFqA867+A8fjr3EwefOohT\nNacbD37jDYiJ0btn3UJMCw6mWbVqzK2A27mRqZEM/Hog8+6Zx5S7p5TomKnBwbSsVo13S3E+ymhW\n2IT8fN0YDg/X/6anQ0YGkpXF8eRAdkf/iWenUfTr8xhaixa6G7MS1YSvTFxMu8jUn6aSU5DDtw9+\nS2uXsnWoTU7WPb7nzul1jc+f1w3cmBjdU123Lri46H/r1NGvgapV0ysAmkz6tZXRCFlZ+i0jQz8u\nMVE3mps00a+vmjfXayd36KAbyx06/GPT4pZHGc2KYtkSH8/y6GgO3X239VqwXr4Md94Je/fqezpm\n4p0LF7icn89XHTqYbc7yICLM3D2Tk3En2f34bmo6FhPkFRCgB7IFBurugluIc9nZDDh5krBevWwW\nLlMUMRkxDFw/kJd6vcRLvV4q0TEROTl0P36c0FK2Ca9IRrOmaXWBrUALIAJ4VERSixg3ElgG2ANr\nROTjwucXAfcB+UAYMEVE0oo4XunsCs65xHM8sf0JnKo5sfb+tTR3am5rkW45RIQN/ht4fe/rvNr7\nVV7v93qxdd/Li8GgG78pKfotLU03hHNz9dc0Tb/Z2+vxxjVqQO3aUL++fnNyunXCK8qLMpoVxWIU\noZOvL1+2b89gF5ebH2Au1qyB1avh8GGz7PGkFRTQ5sgRfLt3p3X16mYQ0DyYxMS0ndO4mHaRXyb+\nQjWH6zzgBQXQpw88+yxMm2YbIS3M+KAgutWuzRvNK8aPckJWAp4bPJnkPom3BrxV4uP+ExKCk4MD\nC1qXzktUwYzmhUCiiCzUNO1NwEVEZl03xh44BwwFooFjwAQROatp2jBgv4iYNE37COD64wvnUDq7\nElBgKmDhoYV8cuSTEoUMKEpObEYsz+16jojUCL4Z+w1dGprPQaSwLGXV2bdW30ZFkdhrGrNbtGB+\nZKR1F546Vd87+vRTs0z3RXQ0o+rVq1AGM4CdZsfqMatxrenKI98/Qr4x/58Dli3TL/effto2AlqB\nOS1asDQqimyj0daikJyTzIiNI3igwwOlMpij8/LYcvkyrzRtakHprML9wIbC+xuAotp0egDnRSRC\nRAzAFuABABHZK3K1Z/xRoNK/IbczDnYOzB4wm31P7GO573JGfTeKi2kXbS1WpUZE+Mb/G7qs7MJd\nbnfhO81XGcy3Ccpovk2Y6OZGZG4u3qn/2qW1HHZ2urf5ww/1BJZykGU08umlS7xVQTyZ12NvZ883\nY7+hin0Vxm8bj8Fo0F84fx4++kj3uN/C3h33WrXo6+TE6thYm8qRkpPCsG+HMajlID4Y/EGpjl0c\nFcVTDRviVvkD+xqISHzh/XigqHigJkDUNY8vFT53PVOBX80rnsIWdGnYBd9pvvRv1p/uq7qz0m8l\npqvXRoqSEpkayX2b72Opz1J+n/Q78wfPp6pD5a+EpCgZymi+TahiZ8dbzZvzf9b2NrdtC3Pm6GEJ\nprIr6FUxMQxwcuKOClwYsop9FbY8vAWDycD4H8djMOTp5z17tp4GfYszp0ULFl28SK6NvM2puakM\n+3YY97S4h8XDF5dqC/pyfj4b4uJ43QZ1v8uCpml7NU07XcTtH3UerxTxL2KKm8ZVaJo2B8gXke/M\nJLbCxlSxr8KcgXP4c/KfbPDfwMD1Awm6HGRrsSoFRpORZUeW0X1Vd/o07YPvdF/ubnS3rcVSWJmK\nk7WjsDhPNmzI/MhIfNPT8ahTx3oLv/QS/PADfPklvPBCqQ/PNhpZHBXFL+7uFhDOvFR1qMq2cdsY\n98M4vn7Gg6ezq2I3c6atxbIK3WvXpkutWqyPi+P5JkU5LS1Hck4yIzeOpH/z/iwZvqTUMZtLo6IY\n7+ZG40pSO1tEhhX3mqZp8ZqmNRSROE3TGgGXixgWDVx7hdAM3dt8ZY7JwChgyI3kmDdv3tX7np6e\neHp6lkB6ha3p5NqJQ1MP8ZXfV3hu8OSZbs8wZ+AcalSp+E2KbMGx6GM8v+t5aletzeGnD9O+XsWu\nS6/4NwcPHuTgwYPlnkclAt5mrIiOZmdiIrvNWNGiRAQHQ//+esX0du1KdejHFy/il5HBD53L0H3P\nRuQHBpDbpwdvvX8PS1/85bbZvjuans4jQUGEeHhQ3UoFPhOyEhj27TCGtBpSag8zQGxeHp2PHeNU\njx40L2MZwwqYCJgkIh9rmjYLcC4iEdABPRFwCBAD+PJ3IuBIYAlwj4gk3mAdpbNvAWIyYnj191c5\ncukIS0cs5cGOD6pEwUKSspOYvX82O0N28tGQj3iyy5PqvblFUImAihIxrVEjQnNyOJCSYt2FO3aE\nd97RW7QaDCU+LNVgYHFUFPMrYA3gYjEYcJzyNDUWLCauSR3Gbh1LtiHb1lJZhV516tCjdm2+iI62\nynqxGbF4bvBkTPsxZTKYAeZHRjKlYcMyG8wVkI+AYZqmhQCDCx+jaVpjTdN2AYhIATAD+B04A2wV\nkbOFx38G1AL2app2UtO0FdY+AYX1aFy7MVse2cL6B9bz7oF3Gb5xOIGXA20tlk0xGA18dvQzOq3o\nRFWHqpx94SxPdX1KGcwK5Wm+HfkuPp5PL13iSLdu1lUCJpPerrV3b3jvvRIdMic8nLj8fNZ27Ghh\n4czIvHl6i9rffqNAjEz9aSqRaZH8POFn6lS1YliMjTiblcU9p04R4uGBcylqHZeWiNQIhn07jMld\nJjNn4JwyzXE+O5veJ05wrlcv6pVD1orkabYWSmffehiMBr70+5L/++v/GNtxLO8Pep+GtRraWiyr\nISLsCt3F63tfp1mdZiwZvgT3BhU/LFBReqzuadY0bZymaUGaphk1Tet2g3EjNU0L1jQttLBmqMLG\njHdzI1+E/yUWu/NqGezs4Ouv4auv9DCNmxCXl8fKmJgK2WmuWHx89NjtdetA03Cwc+DrsV/TqX4n\nhn4zlMRsK7/nNuCOmjUZU68eC6Oibj64jAReDmTA+gHM7DWzzAYzwLsREbzctGm5DGaF4lahin0V\nXur1EudmnMOpqhN3rriTd/54h9RcK1ZdshFekV4M/Hogb+57kyXDl/D7pN+Vwaz4F+UJzzgNPAj8\nVdyAwgL6nwMjgU7ABE3T7ijHmgozYKdpLGjVijnh4RSUo6JFmWjUCFau1MM00v7VZOwffHDxIk9W\npm3z5GQYPx5WrYLGja8+bafZsWL0Coa3GU6/df2ISI2wnYxWYl7LlnwVE0NMXp7Z5z508RBDvhnC\nwqELmeExo8zznMrI4EBqKi9X/rrMCoVZcanuwqLhi/B7xo+YjBjafdaOD/76gPS8dFuLZnYORx3m\n3k338tSOp3im2zMEPBfAqHajVCiGokjKbDSLSLCIhNxkWLEF9BW2ZUTdujSqWpV1cXHWX3zsWD1M\nY/JkKGZ792xWFlsuX2Z2Ba3L/C9E4Kmn4JFH4IF/f8U1TeP/Bv8fL3q8SP91/fGP87eBkNajWbVq\nTG3UiPciIsw67/az23lw64N8M/YbJrhPKPM8IsIb4eHMbt6cWhWo9bdCUZFo6dyStQ+s5fDUwwQn\nBdP609bM3j+b+Mz4mx9cgRER9oXvY9CGQTz+v8d5oMMDBM8I5okuT2BvZ50EZkXlxNKJgCUtoK+w\nMpqmsbhNG969cIGUUiTmmY0lSyAmRv97HSLCS+fP83aLFrhWlkYTS5ZAQgIsWHDDYTM8ZvDJiE8Y\n9u0wdp/fbSXhbMPs5s35KTGRExkZ5Z5LRPjE5xNe/O1Ffnv8N0a0HVGu+bYnJhKTl8dz1+wIKBSK\nomlXrx3fPvgtvtN9SctN444v7mD6zumcijtla9FKRY4hh7Un1tL1q6689NtLTOk6hZAZITzX4zkc\n7SvJb43CptzQxaJp2l6gqCyA2SLycwnmL1WWiKr5aV26167NQ66uzLlwgRXtrVx3smpV+P576NUL\nPDxg4MCrL21PTCQuP58XKotB4+0NixaBry+UwMgf13kcjWo3YtwP43h7wNu84FH62tWVAZcqVfiw\ndWv+ExLC4W7dsCvjdmeBqYBXdr/CgYgDHH76MM2dyrf7kGU08sr582zo2JEqdmXzG5ir5qdCUZlo\n7dKaL0Z/wVzPuaw+vpoxm8fQwqkF07tN5+FOD1PLsZatRSySwMuBrDu5jo0BG+nZpCeLhi1iWOth\nKgRDUWrKXT1D07QDwH9F5EQRr/UG5onIyMLHbwEmEfm4iLEqE9sGJBsMdPL1Zdddd9G9dm3rC7B7\nNzz9tG5wNmlCttFIJ19fvu7YEU8XF+vLU1oiIqBPHz3BcUTpvJ/hKeGM/m40w1oPY+mIpTjY3Xph\nAiYR+p88ydSGDZlWhougpOwkHtv2GA52Dmx9ZCtO1ZzKLdOc8HDCc3PZ3KlTuee6gqqeobgdKTAV\nsPPcTtafWo/3RW/u73A/E+6cwOBWg23uub2YdpEfz/zI5sDNxGTEMLnrZCZ3nUzbum1tKpeiYlBW\nnW0uo/k1ETlexGvFFtAvYqxSwDZiXWwsX8XE4FMOb2C5+Phj2LoV/vqLuQkJnMvOZktlaGSSkQH9\n+ulGfxm7/qXmpjJ+23jyjHlsfWQrbjXdzCyk7TmZkcHIgADOeHiUqkrF6fjTjN06lofveJgFQxaY\nJdYwNDubPidO4N+zJ03M2P1PGc2K2534zHi+O/0d285u40zCGe5tey+j241mcKvBNKrdyOLrF5gK\n8IvxY0/YHnaF7iIsOYyxHccyrtM4hrYeqmKVFf/A6kazpmkPAsuB+kAacFJE7tU0rTGwWkRGF467\nF1gG2ANrRaTIoE+lgG3HFW/g5IYNecYWIREiMG0agSIMmjqVEz160KyiV8wwmeDBB8HNTa+WUY6L\nDaPJyNyDc/nG/xu2PboNjyYeZhS0YvBiaCh5JhOrOnQo0fhNAZt4+feXWTZiGY/f9bhZZBARRgQE\nMMzFhdfNnGCqjGaF4m9iM2L56dxP7Anbw8GIgzSq3YgBzQfQo3EPejTuQSfXTuXyRIsIcZlxnL58\nmiOXjnA0+ig+UT40rdOU4W2GM7LtSO5pcQ9V7FUpSUXR2MzTbC6UArYtpzMzGezvz9Fu3WhdvbrV\n18/PzaX39u28EB/P0y+/bPX1S4WI7lkOCIA9e0oUx1wStp/dzjO/PMM7A9/hRY8Xb6l4u1SDgS5+\nfqxs355769UrdlyOIYeZu2dyMOIgP4z7gS4Nzdfu/YvoaL6Ji8P77rvLHMtcHMpoViiKxmgycjLu\nJD5RPvjF+nEs+hhhKWE0qd2EdvXa0cKpBa41XHGt6YpzNWcc7R2vGtTZhmxyDDmk5qYSmxlLbGYs\nkamRBCcGU8W+Cp1dO9OrSS96N+1N76a9reLRVtwaKKNZUW6WRkXxY0ICf3btioOZjYqb8XZ4OP6p\nqeycNAntySdh1iyrrl8q5s6FnTvhwAFwdjbr1GHJYYz/cTyNajVi3QPrqF+jvlnntyUHUlKYdPYs\np3r0KLIqStDlICb+byKdXDux6r5V1K5qvhj74Kws+p88yeFu3Whfo4bZ5r2CMpoVipJjMBqISI0g\nNDmUqLQoErITuJx1mbS8NAxGA/nGfAShRpUa1HCoQZ2qdWhUuxGNajWimVMzOtbveEvpRoX1UUaz\notyYRBjm788gZ2fetmIXviNpaYwNDORUjx40TEqCe+6BF16AiuhxXrZM7/jn5aWHZlhOgDE3AAAQ\nK0lEQVSAfGM+s/fPZmvQVtaMWVPu8moViTfCwgjNyeF/nTtf9aSbxMSyI8tY4L2ABUMW8PTdT5vV\ny24wmeh78iRTGjbkP00sU/FSGc0KhUJReVBGs8IsXMrNpfvx4/zs7o5HnToWXy/VYKDniRN81Lo1\nD7u66k9evKgbzrNmwbPPWlyGEvPll/DRR7rBbIWmK3vD9jLt52kMbz2cJSOWUKeq5T8PS5NnMtHr\n+HFmNGnCtMaNCUsOY/rP08k35rNh7Aba1G1j9jXfuXCBY+np/HbXXRYLeVFGs0KhUFQeyqqzrbsH\nr6jwNK1WjRXt2/NoUBCxFmiBfC0FJhOPnTnDyLp1/zaYQTdI9+2DDz7QPbu2RgQ+/BAWL4Y//rCK\nwQwwrM0wTj9/Gk3TcP/SnZ+Cf6KyGylV7ezY1KkTb4WH84LXZ/Ra04tR7Ubx5+Q/LWIw/5iQwNdx\ncazv2PGWihFXKBQKhfVRnmZFkbwfEcHPSUkc7NqVmvaWKdXzUmgo57Kz2eXuXnQM9cWLMHKk3nJ7\n0SKwcpw1oBvMr72mJ/zt2QONbJNosj98PzN+m0EblzYsv3c5rV1a20QOc3DgwgGmeK8hrvF49t7Z\nkQEN2llkHb/0dO49fZrf77qLbhauQa48zQqFQlF5UJ5mhVl5p0ULOteoweNnzmC0wA/jyuho9qak\nsLVTp+KTDps317vtHTsGEydCTo7Z5bghGRnw6KPg4wN//mkzgxlgSOsh+D/nz4DmA/BY7cGbe98k\nJSfFZvKUhZCkEMZuGcvUnVNZ1G0sC9p35bnIFNIKCsy+1qXcXMYGBrK6fXuLG8wKhUKhuD1QRrOi\nSDRNY1WHDqQbjbwYGorJjIbzhrg43ouM5Oc778T5Zs0u6tbVPbx2dnrL7bP/6otjGYKCoGdPff0/\n/tD/2hhHe0fe7P8mAc8HkJqbSvvP2/Ox98dk5WfZWrQbEpEawTM/P0PftX3p26wvZ184y7jO43i5\naVMGOTvzUGAgWUaj2daLzcvj3tOnmdm0KWOvDftRKBQKhaIcKKNZUSyOdnb8r3NngrKyeOzMGXLN\nYNgsiYri3QsXONClC21LWvqrWjXYtEmvjTxwIKxZo4dNWAKTCVavBk9PPRHxq6/09SsQjWs35qsx\nX+E1xQu/WD9aftqSuQfmkpidaGvR/sHZhLM8/dPTdF/VHbeabgTPCOaNfm9QzUF/PzVNY1nbtrSo\nVg3PU6e4nJ9f7jVDsrPpd/Ik493ceK1Zs3LPp1AoFArFFVRMs+Km5JlMTA4OJio3l5/c3UvVCvkK\nRhHmhIfzU1ISe+66q+wd/86cgQkTwNUVPvkE3N3LNk9RBATAc8/9bTibc24Lci7xHIsPL2bb2W08\nfMfDTO82HY8mHjZJfCswFfBr6K985vsZp+NP81yP53ip10vUrV68p15EmBcRwab4eH676y7albGO\n8rH0dO4PDGR+y5ZMs3JnSxXTrFAoFJUHVXJOYVFMIrwVHs7Wy5dZ2rYtD9avX2KjzD8zk+nnzlHD\nzo5tnTtTv7wd9AoKdA/w++/DQw/BW2+Vr6LFuXOwdCls3w7z58P06bZJOiwncZlxfH3qa1afWE0t\nx1o87v44D93xEG3rtrXouiKCb7Qv353+jq1BW2np3JIZHjMY12kcVR2qlnie1TExzLlwgfmtWjGt\nUSPsS/j9yjUa+Tgqis8uXWJdx47cX9/6TQ+U0axQKBSVB2U0K6zC/pQUZoaG0qhqVRa1bk2XWrWK\nNZ4jc3P5PDqaDXFxLGjdmikNG2JnTu9nSopeCm7dOujbV6/pPHx4ydpaZ2XpHf1Wr9YT/Z5/HmbM\n0D3YlRyTmPQ21EE/sD14Ow1qNWBkm5F4tvSkX/N+Zqn3HJMRg1ekF7vDdrMnbA+1HWsz0X0iE+6c\nQLt6Za+GcSojgxmhoeSaTCxv144+deoU+/0yivB7cjIvnz9P55o1+bRtW5rbKJRGGc0KhUJReVBG\ns8JqGEwmvoyJYVFUFFU0jXvr1qWfkxMA+SYTMfn57EhM5EJuLuNcXZnbsiUNyutdvhHZ2bB1K6xa\nBadP6wl8/ftDmzZ6m2tnZ91IvnQJoqLgyBE4elQf9+ij8OSTYIHWyhUBo8mIzyUf9ofv52DkQY5F\nH6OFcwvc3dy50+1O2ri0udqe1qmaE472jjjaO2IwGsjIzyA9L524zDgupFzgQuoFghKC8IvxI68g\njz7N+jCyzUhGtB1hVm+2iLApPp63L1zAXtO4r149hrq4UNPeHqMI6UYjvycnszMxkQaOjnzQqhX3\n2cC7fC3KaFYoFIrKg9WNZk3TxgHzgI5ATxE5Ucy4CCAdMAIGEfEoZpxSwJUMESEwK4vdycn4ZWRg\nr2lU0TRcHBwYU78+9zg5FV9OzlKkpemeY29v3UhOSYHUVN0obtpUv3XtCoMHw21YiiyvII+ziWcJ\nvBzI6fjTRKRFEJsRS0xGDBn5GRiMBvKMeVSxq0KdqnWoXbU2bjXdaOXcilbOrehYvyM9m/SkhVML\ni8dMiwins7L4JSmJP1NTMYhgB1Szs2OQiwtj69enTfXqFpWhpCijWaFQKCoPtjCaOwIm4Cvgvzcw\nmi8A3UUk+SbzKQWsUCgqJcpoVigUispDWXW2Q1kXFJHgKwuXgNvqx0ShUCgUCoVCcWthjb1zAfZp\nmuanadp0K6ynUCgUCoVCoVCYlRt6mjVN2ws0LOKl2SLycwnX6CcisZqmuQJ7NU0LFhGvogbOmzfv\n6n1PT088PT1LuIRCoVBYj4MHD3Lw4EFbi1EkmqbVBbYCLYAI4FERSS1i3EhgGWAPrBGRj697/b/A\nIqD+zcLrFAqF4nag3NUzNE07wA1imq8bOxfIFJElRbym4uMUCkWlpCLFNGuathBIFJGFmqa9CbiI\nyKzrxtgD54ChQDRwDJggImcLX28GrAY6UExOitLZCoWislJWnW2u8IwiF9Y0rYamabUL79cEhgOn\nzbSmQqFQKP7N/cCGwvsbgLFFjPEAzotIhIgYgC3AA9e8vhR4w6JSKhQKRSWjzEazpmkPapoWBfQG\ndmma9lvh8401TdtVOKwh4KVp2ingKPCLiOwpr9AKhUKhKJYGIhJfeD8eaFDEmCZA1DWPLxU+h6Zp\nDwCXRCTAolIqFApFJaM81TO2A9uLeD4GGF14PxzoWmbpFAqFQvEvbpBvMufaByIimqYVFUNRZFyF\npmnVgdnAsGufLk4OlYeiUCgqA+bKQ1EdARUKhaKcVLCY5mDAU0TiNE1rBBwQkY7XjekNzBORkYWP\n30Kvu78L2A9kFw5tih7z7CEil6+bQ+lshUJRKbF1TLNCoVAoKgY7gacK7z8F7ChijB/QTtO0lpqm\nOQKPATtFJFBEGohIKxFphR620e16g1mhUChuR5TRrFAoFLcWHwHDNE0LAQYXPv5HvomIFAAzgN+B\nM8DWK5UzrkO5khUKhaIQFZ6hUCgU5aQihWdYC6WzFQpFZUWFZygUCoVCoVAoFBZCGc0KhUKhUCgU\nCsVNUEazQqFQKBQKhUJxE5TRrFAoFAqFQqFQ3ARlNCsUCoVCoVAoFDdBGc0KhUKhUCgUCsVNUEaz\nQqFQKBQKhUJxE5TRrFAoFAqFQqFQ3IQyG82api3SNO2spmn+mqb9T9M0p2LGjdQ0LVjTtFBN094s\nu6i3HgcPHrS1CFbldjtfUOesUNxK3I7fbXXOtz632/mWh/J4mvcAnUWkCxACvHX9AE3T7IHPgZFA\nJ2CCpml3lGPNW4rb7Yt6u50vqHNWKG4lbsfvtjrnW5/b7XzLQ5mNZhHZKyKmwodHgaZFDPMAzotI\nhIgYgC3AA2VdU6FQKBQKhUKhsAXmimmeCvxaxPNNgKhrHl8qfE6hUCgUCoVCoag0aCJS/Iuathdo\nWMRLs0Xk58Ixc4BuIvJwEcc/DIwUkemFjycBvUTkxSLGFi+IQqFQVHBERLO1DNZE6WyFQlGZKYvO\ndrjJhMNu9LqmaZOBUcCQYoZEA82uedwM3dtc1Fq31Q+OQqFQVGaUzlYoFLcb5ameMRJ4HXhARHKL\nGeYHtNM0raWmaY7AY8DOsq6pUCgUCoVCoVDYgvLENH8G1AL2app2UtO0FQCapjXWNG0XgIgUADOA\n34EzwFYROVtOmRUKhUKhUCgUCqtyw5hmhUKhUCgUCoVCYeWOgCVpdKJp2vLC1/01TbvbmvJZgpud\ns6Zpjxeea4CmaYc0TbvLFnKak5I2tNE0raemaQWapj1kTfksQQm/256FuzKBmqYdtLKIZqcE3+36\nmqbt1jTtVOE5T7aBmGZD07R1mqbFa5p2+gZjbin9BUpv3w56W+lspbMLX1c6+2aIiFVugD1wHmgJ\nVAFOAXdcN2YU8Gvh/V7AEWvJZ8Nz7gM4Fd4feTuc8zXj/gB+AR62tdxW+JydgSCgaeHj+raW2wrn\nPA9YcOV8gSTAwdayl+OcBwB3A6eLef2W0l+l+JxvqfO+3fS20tlKZ18zRunsm8xpTU9zSRqd3A9s\nABCRo4CzpmkNrCijubnpOYuIj4ikFT4srklMZaKkDW1eBLYBCdYUzkKU5JwnAj+KyCUAEUm0sozm\npiTnHAvUKbxfB0gSPc+hUiIiXkDKDYbcavoLlN6+HfS20tlKZ19B6eyb6C5rGs0laXRS1JjKrIxK\n29zlaYpuElOZuOk5a5rWBP2f9cvCpyp7YH1JPud2QF1N0w5omuanadoTVpPOMpTknFcDnTVNiwH8\ngZlWks1W3Gr6C5TehltfbyudrXT2FZTOvonuumGdZjNT0n+y62t/VuZ/zhLLrmnaIPTOiv0sJ45V\nKMk5LwNmiYhomqb9f3t3rBpFGIVh+P0INumElCJEi5QBBRERtLDRaxDBSrwBCwut9BZEglhqJWgh\n1jZaqoUGjFqkEkQEkRQRj8XsNkYyA8nOMpP3KXe2OIcdPv7ZmfkPO3/zoenS8yHgBM2e5ovAqySv\nq+rjTCubnS493wTeVNX5JMdpdtpZraqfM65tnsaUX2Bu72okuW1m/5+ZbWbv0Oeiucugk3+/c2Ty\n2VB1Gu4yeYlkjWZ64m63EoagS88ngcdN9rIEXEyyXVVD3cO7S8+bwLeq2gK2krwEVoGhBnCXns8A\ndwCq6lOSL8AKzf7tYzS2/AJzG8af22a2mT1lZrdkV5+PZ3QZdPIMuAKQ5DTwo6q+9ljjfmvtOclR\n4Alwuao25lDjfmvtuaqOVdVyVS3TPCN3fcDhC93O7afA2SQLSRZpXjp433Od+6lLz+vABYDJc2Ir\nwOdeq+zX2PILzO2DkNtmtpk9ZWa3ZFdv/zRX1e8k00EnC8CDqvqQ5Nrk+P2qep7kUpIN4Bdwta/6\nZqFLz8At4DBwb3IVv11Vp+ZV81517HlUOp7b60leAO+AP8BaVQ02gDv+zneBh0ne0lyg36iq73Mr\neo+SPALOAUtJNoHbNLdwR5lfYG5zAHLbzDazJ8fN7A7Z5XATSZIkqUWvw00kSZKkIXLRLEmSJLVw\n0SxJkiS1cNEsSZIktXDRLEmSJLVw0SxJkiS1cNEsSZIktfgLayFKnPGe0gAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106d53550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define a figure to take two plots\n",
    "fig = pl.figure(figsize=[12,3])\n",
    "# Add subplots: number in y, x, index number\n",
    "ax = fig.add_subplot(121,autoscale_on=False,xlim=(0,1),ylim=(-2,2))\n",
    "ax.set_title(\"Eigenvectors for changed system\")\n",
    "ax2 = fig.add_subplot(122,autoscale_on=False,xlim=(0,1),ylim=(-0.05,0.05))\n",
    "ax2.set_title(\"Difference to perfect eigenvectors\")\n",
    "for m in range(4): # Plot the first four states\n",
    "    psi = np.zeros(num_x_points)\n",
    "    for i in range(num_basis):\n",
    "        psi = psi+eigvec[i,m]*basis_array[i]\n",
    "    if eigvec[m,m] < 0:  # This is just to ensure that psi and the basis function have the same phase\n",
    "        psi = -psi\n",
    "    ax.plot(x,psi)\n",
    "    psi = psi - basis_array[m]\n",
    "    ax2.plot(x,psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Making the potential larger\n",
    "\n",
    "We will now make the magnitude of the potential larger, so that we're no longer in the weak perturbation regime.  Look at the change to the diagonal of the Hamiltonian, and to the overall eigenvalues.  Are they consistent ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Full Hamiltonian\n",
      "   4.935   -3.603    0.000   -0.288    0.000   -0.079    0.000   -0.033    0.000   -0.017\n",
      "  -3.603   19.739   -3.891   -0.000   -0.368   -0.000   -0.112   -0.000   -0.049    0.000\n",
      "   0.000   -3.891   44.413   -3.970   -0.000   -0.400   -0.000   -0.129    0.000   -0.059\n",
      "  -0.288   -0.000   -3.970   78.957   -4.003   -0.000   -0.417   -0.000   -0.138   -0.000\n",
      "   0.000   -0.368   -0.000   -4.003  123.370   -4.019    0.000   -0.426   -0.000   -0.144\n",
      "  -0.079   -0.000   -0.400   -0.000   -4.019  177.653   -4.029   -0.000   -0.432    0.000\n",
      "   0.000   -0.112   -0.000   -0.417    0.000   -4.029  241.805   -4.035    0.000   -0.436\n",
      "  -0.033   -0.000   -0.129   -0.000   -0.426   -0.000   -4.035  315.827   -4.039    0.000\n",
      "   0.000   -0.049    0.000   -0.138   -0.000   -0.432    0.000   -4.039  399.719   -4.042\n",
      "  -0.017   -0.000   -0.059   -0.000   -0.144    0.000   -0.436    0.000   -4.042  493.480\n",
      "Perturbation matrix elements:\n",
      "   0.000   -3.603   -0.000   -0.288   -0.000   -0.079   -0.000   -0.033    0.000   -0.017\n",
      "  -3.603   -0.000   -3.891   -0.000   -0.368   -0.000   -0.112    0.000   -0.049   -0.000\n",
      "  -0.000   -3.891   -0.000   -3.970   -0.000   -0.400   -0.000   -0.129    0.000   -0.059\n",
      "  -0.288   -0.000   -3.970   -0.000   -4.003   -0.000   -0.417   -0.000   -0.138   -0.000\n",
      "  -0.000   -0.368   -0.000   -4.003    0.000   -4.019    0.000   -0.426    0.000   -0.144\n",
      "  -0.079   -0.000   -0.400   -0.000   -4.019    0.000   -4.029   -0.000   -0.432    0.000\n",
      "   0.000   -0.112    0.000   -0.417    0.000   -4.029    0.000   -4.035    0.000   -0.436\n",
      "  -0.033    0.000   -0.129   -0.000   -0.426   -0.000   -4.035    0.000   -4.039    0.000\n",
      "   0.000   -0.049    0.000   -0.138    0.000   -0.432    0.000   -4.039    0.000   -4.042\n",
      "  -0.017   -0.000   -0.059   -0.000   -0.144    0.000   -0.436    0.000   -4.042    0.000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEMdJREFUeJzt3W2MXNddx/Hfj6R9EUQJlaM43hS5oFAoKqVFuIGCMgg7\nMn1BFRBUPHVt8qICNY2lSjQtUr01wq4tNULIIaqjNuIFJKp4dGz3wUBHrVAbsOLdmqZWY6gluxRD\nW8qTKnmS/HkxM96b8az3zsw99/H7kVba2ZnMOTna/fu3Z//3XEeEAADt9h1VTwAAkB7FHgA6gGIP\nAB1AsQeADqDYA0AHUOwBoAMWLva2P2r7iu1zma+90vZp21+2/Snbty46DgBgfkUk+8cl7Z742kOS\nTkfED0j629FjAEBFXMRFVba3S3oqIl43enxe0j0RccX2Vkn9iPjBhQcCAMwl1Z797RFxZfT5FUm3\nJxoHAJBD8j/QxvBXB85kAIAK3Zzofa/Y3hoR/2b7Dkn/PvkC2/wDAABziAjP+t+kSvbHJS2PPl+W\n9FfTXhQRfERo//79lc+hLh+sBWvBWqx/XL0aWlkJ3XZb6PHHQy++OH9GXjjZ235C0j2Stti+JOn9\nkj4o6WO275d0UdIvLzoOAHTJ2pq0Z490xx3S2bPS0tJi77dwsY+IX9ngqZ2LvjcAdM1gIB06JB09\nKh05Ii0vS5550+Z6qfbsMYNer1f1FGqDtVjHWqzrylpk0/wzz0h33lncexfSZz/XwHZUNTYA1Mlg\nIB08uJ7m9+zZOM3bVszxB1qSPQBUaHJvvsg0n8VBaABQgcFAOnBA2rVLete7pJMn0xV6iWQPAKVL\nuTe/EZI9AJRkMJA+8AFp585y0nwWyR4ASjBO81u3pt2b3wjJHgASmkzzp06VX+glkj0AJFN1ms8i\n2QNAweqS5rNI9gBQoHGa37ZNWl1d/EybopDsAaAAV68O0/yuXdKDD0onTtSn0EskewBY2Nra8MCy\npaXy+uZnRbIHgDll9+b37Rum+ToWeolkDwBzKetMm6KQ7AFgBlVeBbsIkj0A5NS0NJ9FsgeATTQ1\nzWeR7AHgBrJ9801L81kkewCYYpzms33zTS30EskeAK4zuTdfp4uj5kWyB4CRyTR/8mQ7Cr1EsgcA\nSdXcPapMJHsAndaGTps8SPYAOqvJffOzItkD6JyupPkskj2ATllbk/burcfdo8pEsgfQCYOBdODA\nsNPmgQe6keazSPYAWq/tnTZ5kOwBtFYX9+Y3QrIH0Epd6rTJg2QPoFVI89OR7AG0Bml+YyR7AI03\n7rQhzW+MZA+g0Ujz+ZDsATTStBMqKfQbI9kDaJzs3aO62jc/K5I9gMbIXgXbhrtHlYlkD6ARuAp2\nMSR7ALVG33wxSPYAaotOm+KQ7AHUDp02xSPZA6iVyTTflht+V41kD6AWpqV5Cn1xSPYAKkeaT49k\nD6Ay0zptKPRpkOwBVIJOm3IlLfa2L0r6b0kvSBpExI6U4wGov8FAOnhQOnpUOnJkWPDtqmfVfqmT\nfUjqRcQ3E48DoAFI89UpY8+ef7OBjsueacNVsNUoI9n/je0XJH04Ih5LPB6AmuFMm3pIXezfHBFf\ns32bpNO2z0fEZ8dPrqysXHthr9dTr9dLPB0AZWFvvhj9fl/9fn/h93FELD6bPAPZ+yX9b0R8aPQ4\nyhobQLmyaf7YMdJ8kWwrImb+ZzPZnr3tW2x/1+jz75R0r6RzqcYDUD1OqKyvlNs4t0v6Sw9/b7tZ\n0p9ExKcSjgegQnTa1Ftp2zjXDcw2DtAKg4F06NBwb/7wYfbmU5t3G4craAHMjU6b5uBsHAAzo2++\neUj2AGZCmm8mkj2AXLh7VLOR7AFsanV1mOaXljhvvqlI9gA2NE7z994r7dsnnThBoW8qkj2Aqbh7\nVLuQ7AG8BPeCbSeSPYBrxml+2zbSfNuQ7AFcl+bZm28fkj3QcezNdwPJHuioaSdUUujbi2QPdBAn\nVHYPyR7okOzePGfadAvJHugIOm26jWQPtBydNpBI9kCrkeYxRrIHWih73jxpHhLJHmgdzpvHNCR7\noCWm9c1T6DFGsgdagL55bIZkDzQY94JFXiR7oKHYm8csSPZAw5DmMQ+SPdAgpHnMi2QPNACdNlgU\nyR6oOTptUASSPVBTpHkUiWQP1BBpHkUj2QM1wnnzSIVkD9QEJ1QiJZI9UDHOm0cZSPZAhSb35iny\nSIVkD1RgWqcNhR4pkeyBktFpgyqQ7IGS0DePKpHsgRJkO21WV9myQflI9kBCdNqgLkj2QCLZNM8J\nlagayR4oWPa8+XGap9CjaiR7oECcN4+6ItkDBaDTBnVHsgcWRN88moBkD8yJe8GiSUj2wBzYm0fT\nkOyBGZDm0VQkeyAn0jyaLFmyt73b9nnbz9l+T6pxgNTotEEbJEn2tm+SdFTSTklflfSPto9HxJdS\njAekQqcN2iJVst8h6UJEXIyIgaQnJb010VhA4UjzaJtUe/ZLki5lHl+W9KZEYwGFIs2jjVIV+8jz\nopWVlWuf93o99Xq9RNMBNjcYSIcOSUePSocPDwu+XfWs0HX9fl/9fn/h93FErro825vad0taiYjd\no8fvlfRiRBzOvCZSjA3MI5vmjx0jzaO+bCsiZo4hqfbsz0i6y/Z22y+X9DZJxxONBcxt3DfP3jza\nLsk2TkQ8b/udkj4p6SZJH6ETB3Wzuirt3Ts8b569ebRdkm2cXAOzjYOKXL063Jt/5BHpyBFpeZm9\neTTHvNs4XEGLTllbGxb3paVhmucWgegKzsZBJ2TvBbtvH/eCRfeQ7NF63AsWINmjxbJpnnvBoutI\n9milyatg2bJB15Hs0SqTaf7kSQo9IJHs0SKcNw9sjGSPxuOESmBzJHs0GidUAvmQ7NFI0/bmKfTA\nxkj2aJxs3zydNkA+JHs0xrS+eQo9kA/JHo1A3zywGJI9am1apw2FHpgdyR61RacNUBySPWqHvnmg\neCR71AppHkiDZI9a4F6wQFoke1SONA+kR7JHZcZpftcu0jyQGskeleCESqBcJHuUik4boBoke5SG\nvXmgOiR7JEeaB6pHskdSpHmgHkj2SILz5oF6IdmjcJxQCdQPyR6FmZbmKfRAPZDsUQj65oF6I9lj\nIXTaAM1Assfc6LQBmoNkj5ll9+ZJ80AzkOwxk3Ga37aNThugSUj2yGWy0+bECQo90CQke2yKvnmg\n+Uj22NDVq9d32lDogWYi2WMqOm2AdiHZ4yXomwfaiWSPa0jzQHuR7EGaBzqAZN9xpHmgG0j2HUWa\nB7qFZN9BpHmge0j2HcLdo4DuItl3RPZMG86bB7qHZN9y0860odAD3UOybzFOqAQwliTZ216xfdn2\n2dHH7hTjYLrBQDpwYL3ThhMqAaRK9iHp4Yh4ONH7YwN02gCYJuWevRO+NyZMO6GSQg9gLOWe/QO2\n3y7pjKR3R8S3Eo7VaWtr0vLy+t48RR7ApLmLve3TkrZOeep3JT0q6cDo8e9J+pCk+ydfuLKycu3z\nXq+nXq8373Q6aTCQDh6Ujh6VjhwZbt+Y36eAVun3++r3+wu/jyNi8dncaAB7u6SnIuJ1E1+P1GO3\nWXZv/tgx0jzQFbYVETPHulTdOHdkHt4n6VyKcboo2zfP3jyAvFLt2R+2/aMaduV8RdI7Eo3TKfTN\nA5hX8m2cDQdmGye38d78I48M9+aXl9mbB7pq3m0crqCtucm+edI8gHlwNk5NTV4Fe/IkhR7A/Ej2\nNcRVsACKRrKvEe4eBSAVkn1NTHbaUOQBFIlkXzHOmwdQBpJ9hei0AVAWkn0Fpt0LlkIPICWSfcmy\naZ57wQIoC8m+JHTaAKgSyb4E9M0DqBrJPqFpe/MUegBVINknQqcNgDoh2ReMThsAdUSyLxDnzQOo\nK5J9AcYnVGavgqXQA6gTkv2C2JsH0AQk+zlN65un0AOoK5L9HOibB9A0JPsZZDttuAoWQJOQ7HOi\n0wZAk5HsNzHtvHkKPYCmIdnfQDbNc0IlgCYj2U8xrW+eQg+gyUj2EzhvHkAbkexHOG8eQJuR7EXf\nPID263Sy57x5AF3R2WTPmTYAuqRzyZ7z5gF0UaeSPWkeQFd1ItlzQiWArmt9sqfTBgBanOzpmweA\nda1M9qR5AHipViX77Jk2pHkAWNeaZM+ZNgCwscYne/bmAWBzjU727M0DQD6NTPakeQCYTeOSPWke\nAGbXmGRPpw0AzK8RyZ57wQLAYmqd7CdPqOResAAwn9ome06oBIDi1C7Zc948ABRv7mJv+5dsf9H2\nC7bfOPHce20/Z/u87XvzvufamrRjh/T008O9+T17JHveGQIAxhZJ9uck3SfpM9kv2n6tpLdJeq2k\n3ZL+yPYNx+l633y/3696CrXBWqxjLdaxFoubu9hHxPmI+PKUp94q6YmIGETERUkXJO3Y6H2yaf7s\nWWnv3u6leb6R17EW61iLdazF4lLs2W+TdDnz+LKkqbvuXU7zAFCmG3bj2D4taeuUp94XEU/NME5M\n++I4zVPkASAtR0ytw/nfwP60pHdHxDOjxw9JUkR8cPT4E5L2R8TTE//dYgMDQEdFxMyb3UX12WcH\nPi7pT20/rOH2zV2S/mHyP5hnsgCA+SzSenmf7UuS7pZ00vbHJSkinpX0MUnPSvq4pN+ORX99AAAs\nZOFtHABA/SW/gtb27tHFVc/Zfs8Gr/nD0fNrtt+Qek5V2WwtbP/aaA2+YPvvbf9IFfMsQ57vi9Hr\nftz287Z/ocz5lSnnz0jP9lnb/2S7X/IUS5PjZ2SL7U/YXh2txZ4Kppmc7Y/avmL73A1eM1vdjIhk\nH5Ju0rDPfrukl0lalfRDE695i6RTo8/fJOnzKedU1UfOtfgJSd89+nx3l9ci87q/k3RC0i9WPe8K\nvy9ulfRFSXeOHm+pet4VrsWKpEPjdZD0DUk3Vz33BGvx05LeIOncBs/PXDdTJ/sdki5ExMWIGEh6\nUsOLrrJ+XtIfS1IMO3ZutX174nlVYdO1iIjPRcR/jR4+LamtTal5vi8k6QFJfybpP8qcXMnyrMWv\nSvrziLgsSRHx9ZLnWJY8a/E1Sa8Yff4KSd+IiOdLnGMpIuKzkv7zBi+ZuW6mLvZLki5lHk+7wGra\na9pY5PKsRdb9kk4lnVF1Nl0L20sa/qA/OvpSW/+4lOf74i5Jr7T9adtnbP9GabMrV561eEzSD9v+\nV0lrkh4saW51M3PdTH3Ecd4f0Mk2zDb+YOf+f7L9M5J+U9Kb002nUnnW4g8kPRQRYdu6/nukLfKs\nxcskvVHSz0q6RdLnbH8+Ip5LOrPy5VmL90lajYie7e+XdNr26yPifxLPrY5mqpupi/1XJb0q8/hV\neulRCtNec+foa22TZy00+qPsY5J2R8SNfo1rsjxr8WOSnhzWeW2R9HO2BxFxvJwplibPWlyS9PWI\n+Lakb9v+jKTXS2pbsc+zFj8p6fclKSL+2fZXJL1G0plSZlgfM9fN1Ns4ZyTdZXu77ZdreBrm5A/r\ncUlvlyTbd0v6VkRcSTyvKmy6Fra/V9JfSPr1iLhQwRzLsulaRMT3RcSrI+LVGu7b/1YLC72U72fk\nryX9lO2bbN+i4R/kni15nmXIsxbnJe2UpNEe9Wsk/Uups6yHmetm0mQfEc/bfqekT2r4l/aPRMSX\nbL9j9PyHI+KU7bfYviDp/yTtTTmnquRZC0nvl/Q9kh4dJdpBRGx4YmhT5VyLTsj5M3J+dOzIFyS9\nKOmxGF682Co5vy8OSnrc9pqGYfV3IuKblU06EdtPSLpH0pbRxav7NdzOm7tuclEVAHRA7W5LCAAo\nHsUeADqAYg8AHUCxB4AOoNgDQAdQ7AGgAyj2ANABFHsA6ID/B8tmCmxudUz4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10884fc10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Vdiag2 = np.linspace(0.0,width,num_x_points)\n",
    "square_well_potential(x,Vdiag2,20.0)\n",
    "# Declare space for the matrix elements\n",
    "Hmat3 = np.eye(num_basis)\n",
    "\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "print \"Full Hamiltonian\"\n",
    "for n in range(num_basis):\n",
    "    # Loop over basis functions phi_m (the ket in the matrix element)\n",
    "    for m in range(num_basis):\n",
    "        # Act with H on phi_m and store in H_phi_m\n",
    "        H_phi_m = -0.5*d2basis_array[m] \n",
    "        add_pot_on_basis(H_phi_m,Vdiag2,basis_array[m])\n",
    "        # Create matrix element by integrating\n",
    "        Hmat3[m,n] = integrate_functions(basis_array[n],H_phi_m,num_x_points,dx)\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % Hmat3[m,n],\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print\n",
    "    \n",
    "print \"Perturbation matrix elements:\"\n",
    "# Output the matrix elements of the potential to see how large the perturbation is\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "for n in range(num_basis):\n",
    "    # Loop over basis functions phi_m (the ket in the matrix element)\n",
    "    for m in range(num_basis):\n",
    "        # Act with H on phi_m and store in H_phi_m\n",
    "        H_phi_m = np.zeros(num_x_points)\n",
    "        add_pot_on_basis(H_phi_m,Vdiag2,basis_array[m])\n",
    "        # Create matrix element by integrating\n",
    "        H_mn = integrate_functions(basis_array[n],H_phi_m,num_x_points,dx)\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % H_mn,\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Eigenvalues and eigenvector coefficients printed roughly\n",
      "[   4.0834   19.9755   44.5683   79.0523  123.4335  177.6979  241.8388\n",
      "  315.8539  399.7402  493.6551]\n",
      "[  9.7302e-01   2.3030e-01  -1.3321e-02  -4.2000e-03   2.8058e-04\n",
      "   4.6696e-04  -2.4300e-05   1.0562e-04  -4.0734e-06   3.3954e-05]\n",
      "[  2.2955e-01  -9.6072e-01   1.5575e-01   7.1874e-03  -3.6951e-03\n",
      "  -2.1336e-04   5.0959e-04  -2.2948e-05   1.2990e-04  -6.5208e-06]\n",
      "[  2.2642e-02  -1.5451e-01  -9.8109e-01  -1.1424e-01   4.3154e-03\n",
      "   3.0819e-03  -1.5354e-04   4.7700e-04  -1.8763e-05   1.3122e-04]\n",
      "\n",
      "    Diag  Perf Square  Difference\n",
      "   4.083        4.935      -0.851\n",
      "  19.976       19.739       0.236\n",
      "  44.568       44.413       0.155\n",
      "  79.052       78.957       0.095\n",
      " 123.434      123.370       0.063\n",
      " 177.698      177.653       0.045\n",
      " 241.839      241.805       0.033\n",
      " 315.854      315.827       0.027\n",
      " 399.740      399.719       0.021\n",
      " 493.655      493.480       0.175\n"
     ]
    }
   ],
   "source": [
    "# Solve using linalg module of numpy (which we've imported as la above)\n",
    "eigval, eigvec = la.eigh(Hmat3)\n",
    "# This call above does the entire solution for the eigenvalues and eigenvectors !\n",
    "# Print results roughly, though apply precision of 4 to the printing\n",
    "print\n",
    "print \"Eigenvalues and eigenvector coefficients printed roughly\"\n",
    "np.set_printoptions(precision=4)\n",
    "print eigval\n",
    "print eigvec[0]\n",
    "print eigvec[1]\n",
    "print eigvec[2]\n",
    "print\n",
    "\n",
    "print \"    Diag  Perf Square  Difference\"\n",
    "for i in range(num_basis):\n",
    "    n = i+1\n",
    "    print \"%8.3f     %8.3f    %8.3f\" % (eigval[i],n*n*np.pi*np.pi/2.0,eigval[i] - n*n*np.pi*np.pi/2.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "So for this problem first order perturbation theory is not good enough.  We could use higher order perturbation theory here, or a different approach.  We will return to this problem in a notebook on the variational theorem."
   ]
  },
  {
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}