{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving for eigenvalues and eigenvectors\n",
    "\n",
    "We saw in the previous notebook how to set up a matrix by integrating between pairs of basis functions.  In this notebook we will solve for the eigenvalues and eigenvectors of the Hamiltonian for a square well, exploring what happens when we use a finite basis set, and introduce changes to the potential in the square well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenbasis, adding potential\n",
    "\n",
    "We will start by using the eigenbasis (since it's simple) and will add small changes to the potential.  This will give us some insight into the effect of perturbations (and will feed directly into our studies of perturbation theory).  So we'll first set up the bits and pieces of code we'll need, some taken from the previous notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now set up all the mechanics that we need to create matrices for different Hamiltonians, except for one thing: the potential.  We have an implicit potential already in the infinite square well (we set $V=\\infty$ when $x=0$ or $x=a$).  If we create a potential function, then we can change this and create different Hamiltonian matrices.  So we'll do this next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEVVJREFUeJzt3WuMXGd5wPH/43UigSouUSQqOYYATWkiFQRV3CDaMij5\nsCCBESClLhRxawOSERJFDVQV2YhKiPAlqkxTF6UXvtSqaEUMiola0SkI5WYgCW3tyG6IajsqBEIR\nJank9T79MBPvsN2dObM7M7tznv9PWnkux+M3R2f/fvPuO+PITCRJ7bFruwcgSZoswy5JLWPYJall\nDLsktYxhl6SWMeyS1DIjwx4RixFxMiJORcTN6zzfiYifRMR3+l9/PJ2hSpKa2D3syYhYAA4BNwDn\ngAcj4mhmnlhz6L9k5lumNEZJ0hhGzdj3Aacz8/HMPA8cAfavc1xMfGSSpE0ZFfY9wJmB+2f7jw1K\n4LUR8VBE3B0R10xygJKk8QxdiqEX7VG+Dbw4M5+OiDcCXwJ+ecsjkyRtyqiwnwP2DtzfS2/WflFm\n/nTg9rGI+LOIuCwznxo8LiL8UBpJ2oTMHGu5e9RSzHHgqoi4MiIuBW4Ejg4eEBEviojo394HxNqo\nDwzOr0xuueWWbR/DTvnajnPxsY8lt922/nPXXps88ECdc7FTvzwXq1+bMXTGnpnLEXEQuAdYAO7M\nzBMRcVP/+cPAO4APRcQy8DTw25saiTQjy8uwe4Mrf/fu3vPSPBu1FENmHgOOrXns8MDtzwGfm/zQ\npOlYXoaFhfWfW1gw7Jp/vvN0G3Q6ne0ewo6xHedip87YvS5WeS62xrBvAy/aVYZ9ldfFKs/F1hh2\nlbNTwy5NimFXORcuDA/7hQuzHY80aYZd5ThjV9sZdpVj2NV2hl3luN1RbWfYVY4zdrWdYVc5hl1t\nZ9hVjrti1HaGXeU4Y1fbGXaVY9jVdoZd5bgrRm1n2FWOM3a1nWFXOYZdbWfYVY67YtR2hl3lOGNX\n2xl2lWPY1XaGXeUMC7u7YtQGhl3lDNvu6IxdbWDYVY5LMWo7w65yDLvazrCrHLc7qu0Mu8pxxq62\nM+wqx10xajvDrnLcFaO2M+wqx6UYtZ1hVzmGXW1n2FXKygpEwK4Nrnx3xagNDLtKGTZbB2fsagfD\nrlJGhd1dMWoDw65Shu2IAWfsagfDrlJcilEFhl2lGHZVYNhVyrDPiQF3xagdDLtKccauCkaGPSIW\nI+JkRJyKiJuHHHdtRCxHxNsmO0Rpcgy7Khga9ohYAA4Bi8A1wIGIuHqD4z4DfBWIKYxTmgi3O6qC\nUTP2fcDpzHw8M88DR4D96xz3YeCLwJMTHp80UW53VAWjwr4HODNw/2z/sYsiYg+92N/RfygnNjpp\nwlyKUQVDLnGgWaRvBz6emRkRwZClmKWlpYu3O50OnU6nwctLk2PYtdN1u1263e6WXiMyN253RFwH\nLGXmYv/+J4CVzPzMwDGPsRrzy4Gngd/LzKNrXiuH/VnSLBw/Dh/8YO/X9Xzve3D99fDYY7Mdl7SR\niCAzx/rZ5agZ+3Hgqoi4EngCuBE4MHhAZr5sYAB/BXx5bdSlncIZuyoYGvbMXI6Ig8A9wAJwZ2ae\niIib+s8fnsEYpYlxV4wqGDVjJzOPAcfWPLZu0DPzvRMalzQV7opRBb7zVKW4FKMKDLtKMeyqwLCr\nFD8ETBUYdpXijF0VGHaV4q4YVWDYVcqosO/aBSsrvS9pXhl2lTJqu2OE6+yaf4ZdpYyasYPr7Jp/\nhl2ljNoVA87YNf8Mu0pxxq4KDLtKaRJ2d8Zo3hl2leKMXRUYdpUyalcMGHbNP8OuUpyxqwLDrlIM\nuyow7CrF7Y6qwLCrFHfFqALDrlJcilEFhl2luCtGFRh2leKMXRUYdpVi2FWBYVcp7opRBYZdpThj\nVwWGXaW43VEVGHaV4oxdFRh2leJ2R1Vg2FWKM3ZVYNhVirtiVIFhVynO2FWBYVcp7opRBYZdpThj\nVwWGXaW4K0YVGHaV4oxdFRh2lWLYVYFhVylud1QFI8MeEYsRcTIiTkXEzes8vz8iHo6I70TEgxHx\nuukMVdo6d8WogqGXeEQsAIeAG4BzwIMRcTQzTwwc9k+ZeVf/+F8F/g64ekrjlbbEpRhVMGrGvg84\nnZmPZ+Z54Aiwf/CAzPzZwN1fAFYmO0Rpcgy7KhgV9j3AmYH7Z/uP/ZyIeGtEnAC+ArxvcsOTJsvt\njqpgVNizyYtk5pcy82rgrcCfbHlU0pQ4Y1cFIy5xzgF7B+7vpTdrX1dmfiMiXhYRl2XmU2ufX1pa\nuni70+nQ6XTGGqy0Ve6K0U7X7Xbpdrtbeo3I3HhSHhG7gUeB64EngAeAA4M/PI2IlwOPZWZGxGuA\nuzJz7zqvlcP+LGkWXvIS+PrXe79u5LOfhR/8oPertN0igsyMcX7P0LlLZi5HxEHgHmABuDMzT0TE\nTf3nDwNvB94dEeeBZ4AbNzV6aQZcilEFo5ZiyMxjwLE1jx0euH0bcNvkhyZNnmFXBb7zVKW4K0YV\nGHaV4oxdFRh2leKuGFVg2FWKM3ZVYNhVih8CpgoMu8rI7C2x+MNTtZ1hVxkXLsCuXRAj3uph2DXv\nDLvKaLIMA4Zd88+wqwzDrioMu8postUR3O6o+WfYVUbTGbu7YjTvDLvKcClGVRh2lWHYVYVhVxlN\nPgAMDLvmn2FXGc7YVYVhVxnuilEVhl1luCtGVRh2leFSjKow7CrDsKsKw64y3BWjKgy7ynDGrioM\nu8pwV4yqMOwqw10xqsKwqwyXYlSFYVcZhl1VGHaV4a4YVWHYVYYzdlVh2FVG010xCwvuitF8M+wq\no+mMfVf/u2JlZbrjkabFsKuMpmEHl2M03wy7yjDsqsKwqwzDrioMu8pout0RDLvmm2FXGc7YVYVh\nVxlNtzuCHwSm+WbYVcY4M3Y/CEzzrFHYI2IxIk5GxKmIuHmd598ZEQ9HxCMR8c2IeOXkhyptjUsx\nqmJk2CNiATgELALXAAci4uo1hz0G/FZmvhL4FPAXkx6otFWGXVU0mbHvA05n5uOZeR44AuwfPCAz\n783Mn/Tv3g9cMdlhSlvnrhhV0STse4AzA/fP9h/byPuBu7cyKGkanLGriiaXeTZ9sYh4A/A+4HWb\nHpE0JRcuwHOe0+xYd8VonjUJ+zlg78D9vfRm7T+n/wPTzwOLmfnj9V5oaWnp4u1Op0On0xljqNLW\nuCtG86Db7dLtdrf0GpE5fEIeEbuBR4HrgSeAB4ADmXli4JgXA18D3pWZ923wOjnqz5Km6aMfhSuu\n6P06ynXXwe23936VtlNEkJkxzu8ZOX/JzOWIOAjcAywAd2bmiYi4qf/8YeCTwAuBOyIC4Hxm7hv3\nP0CaJtfYVUWjyzwzjwHH1jx2eOD2B4APTHZo0mQZdlXhO09VhtsdVYVhVxl+VoyqMOwqw10xqsKw\nqwzX2FWFYVcZhl1VGHaVYdhVhWFXGe6KURWGXWU4Y1cVhl1ljLPdcWHB7Y6aX4ZdZThjVxWGXWUY\ndlVh2FWGYVcVhl1luCtGVRh2leGMXVUYdpXhrhhVYdhVhjN2VWHYVYZhVxWGXWUYdlVh2FWGYVcV\nhl1luN1RVRh2leE/jacqDLvK8J/GUxWGXWW4xq4qDLvKMOyqwrCrDMOuKgy7ynBXjKow7CrDGbuq\nMOwqYWUFMmFXwyveDwHTPDPsKuHZPewRzY53xq55ZthVwjjLMGDYNd8Mu0ow7KrEsKsEw65KDLtK\nGGerIxh2zTfDrhLG+QAwcFeM5pthVwkuxagSw64SDLsqGRn2iFiMiJMRcSoibl7n+V+JiHsj4n8j\n4g+mM0xpawy7Khl6qUfEAnAIuAE4BzwYEUcz88TAYT8CPgy8dWqjlLbIsKuSUTP2fcDpzHw8M88D\nR4D9gwdk5pOZeRw4P6UxSlvmrhhVMirse4AzA/fP9h+T5oq7YlTJqEs9J/mHLS0tXbzd6XTodDqT\nfHlpQy7FaF50u1263e6WXmPUpX4O2Dtwfy+9WfumDIZdmiXDrnmxdtJ76623jv0ao5ZijgNXRcSV\nEXEpcCNwdINjG35unjR7hl2VDL3UM3M5Ig4C9wALwJ2ZeSIibuo/fzgifhF4EHgesBIRHwGuycz/\nmfLYpcYMuyoZealn5jHg2JrHDg/c/i9+frlG2nHcFaNKfOepSnDGrkoMu0pwu6MqMewqwRm7KjHs\nKsGwqxLDrhLGDfuzSzE50bfoSbNh2FXCuGGPcJ1d88uwq4RxtzuCyzGaX4ZdJYy7KwacsWt+GXaV\nMO5SDDhj1/wy7CrBsKsSw64SDLsqMewqwbCrEsOuEtwVo0oMu0pwV4wqMewqwaUYVWLYVYJhVyWG\nXSUYdlVi2FWCYVclhl0lGHZVYthVwma2Oy4sGHbNJ8OuEjaz3XH3brc7aj4ZdpXgUowqMewqwbCr\nEsOuEgy7KjHsKsGwqxLDrhLcFaNKDLtKcFeMKjHsKsGlGFVi2FWCYVclhl0lGHZVYthVgmFXJYZd\nJfhP46kSw64S/KfxVIlhVwkuxagSw64SDLsqGRn2iFiMiJMRcSoibt7gmD/tP/9wRLx68sOUtsaw\nq5KhYY+IBeAQsAhcAxyIiKvXHPMm4Jcy8yrg94E7pjTW1uh2u9s9hB1jVudiHsLudbHKc7E1o2bs\n+4DTmfl4Zp4HjgD71xzzFuBvADLzfuAFEfGiiY+0RbxoVxn2VV4XqzwXWzMq7HuAMwP3z/YfG3XM\nFVsfmjQ5fgiYKhk1h8mGrxNNft+b39zw1Vru0UfhW9/a7lHsDLM6F2fPwiWXjPd7Lr0UvvAFeOih\n6YxpLa+LVZ6LrYnMjdsdEdcBS5m52L//CWAlMz8zcMyfA93MPNK/fxJ4fWZ+f81rNf1LQpI0IDPX\nTp6HGjVjPw5cFRFXAk8ANwIH1hxzFDgIHOn/RfDfa6O+mYFJkjZnaNgzczkiDgL3AAvAnZl5IiJu\n6j9/ODPvjog3RcRp4GfAe6c+aknShoYuxUiS5s/E33nqG5pWjToXEfHO/jl4JCK+GRGv3I5xTluT\na6J/3LURsRwRb5vl+Gap4fdHJyK+ExH/GhHdGQ9xZhp8fzw/Ir4cEQ/1z8V7tmGYMxERfxkR34+I\n7w45pnk3M3NiX/SWa04DVwKXAA8BV6855k3A3f3bvw7cN8kx7JSvhufitcDz+7cX23gumpyHgeO+\nBnwFePt2j3sbr4kXAP8GXNG/f/l2j3sbz8UfAZ9+9jwAPwJ2b/fYp3Q+fhN4NfDdDZ4fq5uTnrH7\nhqZVI89FZt6bmT/p372fdu7/b3JNAHwY+CLw5CwHN2NNzsXvAH+fmWcBMvOHMx7jrDQ5FyvA8/q3\nnwf8KDNb+c6CzPwG8OMhh4zVzUmH3Tc0rWpyLga9H7h7qiPaHiPPQ0TsofdN/ezHUbT1Bz9Nromr\ngMsi4p8j4nhE/O7MRjdbTc7FIeCaiHgCeBj4yIzGthON1c0x32Q90kTf0DTnGv83RcQbgPcBr5ve\ncLZNk/NwO/DxzMyICP7/9dEWTc7FJcBrgOuB5wL3RsR9mXlqqiObvSbnYhH4dma+ISJeDvxjRLwq\nM3865bHtVI27OemwnwP2DtzfS+9vlmHHXNF/rG2anAv6PzD9PLCYmcP+V2xeNTkPv0bvfRDQW0t9\nY0Scz8yjsxnizDQ5F2eAH2bmM8AzEfF14FVA28Le5Fy8B/g0QGb+R0R8D3gFvffXVDNWNye9FHPx\nDU0RcSm9NzSt/eY8CrwbLr6zdd03NLXAyHMRES8G/gF4V2ae3oYxzsLI85CZL8vMl2bmS+mts3+o\nhVGHZt8fdwG/ERELEfFcej8o+/cZj3MWmpyL/wRuAOivJ78CeGymo9w5xurmRGfs6RuaLmpyLoBP\nAi8E7ujPVs9n5r7tGvM0NDwPJTT8/jgZEV8FHqH3w8PPZ2brwt7wuvgU8NcR8Qi9ZYg/zMyntm3Q\nUxQRfwu8Hrg8Is4At9BblttUN32DkiS1jP80niS1jGGXpJYx7JLUMoZdklrGsEtSyxh2SWoZwy5J\nLWPYJall/g+xvQq90MU3ygAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a6ea750>"
      ]
     },
     "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] = 0.0\n",
    "    # Let's define a small bump in the middle of the well\n",
    "    i_mid = x.size/2\n",
    "    V[i_mid-1] = a\n",
    "    V[i_mid] = a\n",
    "    V[i_mid+1] = a\n",
    "    # Plot to ensure that we know what we're getting\n",
    "    pl.plot(x,V)\n",
    "    \n",
    "# Declare space for this potential (Vbump) and call the routine\n",
    "Vbump = np.linspace(0.0,width,num_x_points)\n",
    "pot_bump = square_well_potential(x,Vbump,0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving for the eigenvalues\n",
    "\n",
    "Now that we have the potential, we need to build the matrix (remember that we have set $\\hbar = m = 1$):\n",
    "\n",
    "$H_{mn} = \\langle \\phi_m \\vert -\\frac{1}{2}\\nabla^2 + \\hat{V}\\vert\\phi_n\\rangle$\n",
    "\n",
    "By diagonalisation, we will find the eigenvalues.\n",
    "\n",
    "It's worth noting that we act with $\\hat{V}$ by multiplication in position representation:\n",
    "\n",
    "$$\\langle x \\vert \\hat{V} \\vert \\phi_m\\rangle = V(x) \\phi_m(x)$$\n",
    "\n",
    "Then we make matrix elements by integration in position representation; as we've defined a grid on which to represent the system, we'll just sum over the values on the grid (though there are much more accurate ways to do numerical integration).\n",
    "\n",
    "$$H_{ij} = \\langle \\phi_i \\vert \\hat{H} \\vert \\phi_j \\rangle = \\int dx \\phi_i^\\star(x) \\hat{H} \\phi_j(x)$$\n",
    "\n",
    "We will print out the matrix of _just_ the potential first: $\\langle \\phi_m \\vert \\hat{V} \\vert \\phi_n\\rangle$.  This will be useful when thinking about perturbation theory (the change in eigenvalues should be the same as the diagonal elements of the potential matrix, to first order in the potential)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Potential matrix elements:\n",
      "   0.030    0.000   -0.030   -0.000    0.030    0.000   -0.030   -0.000    0.029    0.000\n",
      "   0.000    0.000   -0.000   -0.000    0.000    0.000   -0.000   -0.000    0.000    0.000\n",
      "  -0.030   -0.000    0.030    0.000   -0.030   -0.000    0.029    0.000   -0.029   -0.000\n",
      "  -0.000   -0.000    0.000    0.000   -0.000   -0.000    0.000    0.001   -0.000   -0.001\n",
      "   0.030    0.000   -0.030   -0.000    0.030    0.000   -0.029   -0.000    0.029    0.000\n",
      "   0.000    0.000   -0.000   -0.000    0.000    0.001   -0.000   -0.001    0.000    0.001\n",
      "  -0.030   -0.000    0.029    0.000   -0.029   -0.000    0.029    0.000   -0.029   -0.000\n",
      "  -0.000   -0.000    0.000    0.001   -0.000   -0.001    0.000    0.001   -0.000   -0.002\n",
      "   0.029    0.000   -0.029   -0.000    0.029    0.000   -0.029   -0.000    0.028    0.000\n",
      "   0.000    0.000   -0.000   -0.001    0.000    0.001   -0.000   -0.002    0.000    0.002\n",
      "\n",
      "Full Hamiltonian\n",
      "   4.965   -0.000   -0.030   -0.000    0.030    0.000   -0.030   -0.000    0.029   -0.000\n",
      "  -0.000   19.739    0.000   -0.000   -0.000    0.000    0.000   -0.000   -0.000    0.000\n",
      "  -0.030   -0.000   44.443    0.000   -0.030    0.000    0.029   -0.000   -0.029    0.000\n",
      "  -0.000   -0.000    0.000   78.957    0.000   -0.000   -0.000    0.001   -0.000   -0.001\n",
      "   0.030    0.000   -0.030    0.000  123.400    0.000   -0.029    0.000    0.029   -0.000\n",
      "   0.000    0.000    0.000   -0.000    0.000  177.654    0.000   -0.001    0.000    0.001\n",
      "  -0.030    0.000    0.029   -0.000   -0.029   -0.000  241.834   -0.000   -0.029   -0.000\n",
      "   0.000   -0.000   -0.000    0.001    0.000   -0.001   -0.000  315.829   -0.000   -0.002\n",
      "   0.029   -0.000   -0.029   -0.000    0.029    0.000   -0.029   -0.000  399.747    0.000\n",
      "  -0.000    0.000    0.000   -0.001   -0.000    0.001   -0.000   -0.002    0.000  493.482\n"
     ]
    }
   ],
   "source": [
    "# Declare space for the matrix elements - simplest with the identity function\n",
    "Hmat = np.eye(num_basis)\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]\n",
    "        \n",
    "print \"Potential matrix elements:\"\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "# Calculate and output matrix elements of the potential\n",
    "\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,Vbump,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",
    "print\n",
    "print \"Full Hamiltonian\"\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "# Calculate and store the matrix elements for the 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",
    "        # First the kinetic energy\n",
    "        H_phi_m = -0.5*d2basis_array[m] \n",
    "        # Now the potential\n",
    "        add_pot_on_basis(H_phi_m,Vbump,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",
    "        Hmat[m,n] = H_mn\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": "markdown",
   "metadata": {},
   "source": [
    "Notice that two things have changed compared to the perfect square well: first, the diagonal elements are _slightly_ larger; second, there are now off-diagonal elements.  Does it surprise you that these alternate (i.e. only in odd row and columns) ? Think about the symmetries of the system, particularly of the wavefunctions.\n",
    "\n",
    "What effect will this have on the eigenvalues and eigenvectors ? We'll diagonalise, print out the eigenvalues and plot the first few eigenvectors (as well as looking at their coefficients to get a rough idea of how much they've changed)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   4.9647   19.7393   44.443    78.9571  123.3996  177.6536  241.8344\n",
      "  315.8286  399.7474  493.4821]\n",
      "[  1.0000e+00  -2.0635e-17   7.5696e-04  -7.1521e-17  -2.5115e-04\n",
      "   7.7171e-17  -1.2461e-04   3.6019e-17  -7.3983e-05   5.9905e-17]\n",
      "[  2.2443e-17   1.0000e+00  -2.2204e-15  -2.6578e-06   2.2204e-16\n",
      "   1.4901e-06   5.5511e-17   1.0548e-06   3.8858e-16  -8.1916e-07]\n",
      "[  7.5708e-04  -2.2204e-15  -1.0000e+00  -3.4056e-16   3.7569e-04\n",
      "   1.3964e-16   1.4913e-04   1.0053e-16   8.1990e-05  -2.4266e-16]\n"
     ]
    }
   ],
   "source": [
    "# Solve using linalg module of numpy (which we've imported as la above)\n",
    "eigval, eigvec = la.eigh(Hmat)\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",
    "np.set_printoptions(precision=4)\n",
    "print eigval\n",
    "print eigvec[0]\n",
    "print eigvec[1]\n",
    "print eigvec[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the eigenvalues look close to the perfect well values (we'll check them properly below).  The eigenvector coefficients show a single dominant value (corresponding to the unperturbed eigenvector), with very small contributions from other eigenvectors.  Now print the eigenvalues and calculate the change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Changed Original  Difference\n",
      "   4.965    4.935    0.030\n",
      "  19.739   19.739    0.000\n",
      "  44.443   44.413    0.030\n",
      "  78.957   78.957    0.000\n",
      " 123.400  123.370    0.030\n",
      " 177.654  177.653    0.001\n",
      " 241.834  241.805    0.029\n",
      " 315.829  315.827    0.001\n",
      " 399.747  399.719    0.028\n",
      " 493.482  493.480    0.002\n"
     ]
    }
   ],
   "source": [
    "# Now print out eigenvalues and the eigenvalues of the perfect square well, and the difference\n",
    "print \" Changed Original  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": [
    "Compare the differences in eigenvalues due to the small potential we've added to the diagonal terms of the potential matrix.  How close is the agreement ? You might like to change the magnitude of the potential (higher up - it's passed as an argument to the potential function) and see how this agreement changes.\n",
    "\n",
    "Now we'll plot the eigenvectors (after building them) and look at the change with respect to the eigenvectors of the original system.  Remember that any function in this space can be written as:\n",
    "\n",
    "$$\\vert\\psi\\rangle = \\sum_i c_i \\vert \\phi_i \\rangle$$\n",
    "\n",
    "We'll use this to build the eigenfunctions of the perturbed system.  In this case, $c_i$ are the coefficients in the eigenvector of the matrix, and $\\vert\\phi_i\\rangle$ are the basis functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAADSCAYAAACxUP3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4FFXXwH+XbqOEhBaaAiJgAwFRRPOJBUFAVAQ76mvB\n/tpQbGBBVFBEEPUVUUSwgIAKNlR67whEOglJCC0kgfTd8/1xJrCETbJJdrMp9/c882R25raZzJw5\n995zzjUigsVisVgsFovFYsmdCsFugMVisVgsFovFUtKxSrPFYrFYLBaLxZIPVmm2WCwWi8VisVjy\nwSrNFovFYrFYLBZLPlil2WKxWCwWi8ViyQerNFssFovFYrFYLPlgleYgYYy53RjzW7DbURowxvQx\nxkQbY5KNMRcEsJ4IY0x0oMr3F6WlnRZLWcEYM84Y85LH74HGmHhjTJIxppYxprMxZqsjo3oFs61l\nhZJyT40x/xhjLg9W/ZaShVWaA4gxZpcxJsV56bO30QAi8rWIXBvsNhYGY8wAY8yCYqxyBPCwiJwh\nIuuKsd5yjzGmqTHGbYyxssJSJvGQ00nGmARjzCJjzIPGGJOdRkQGisgbTvrKwEigq4hUF5EE4DVg\ntCOjfgzOlRQN5z5cGex2eOCXe1rU6xKRc0VkfmHzB5MS+D8t9VQKdgPKOAJcLyJ/BbshJQljTEUR\ncfmY1gCNgU2FrKuCiLgLk9dyAib/JBZLqeSYnDbGnAFEAB8AFwP3eklfD6gGbPY4VhQZ5bM8DDBC\nCXjPjTGVRCSLItzTHJSI6woShb52j/+DxRMRsVuANmAncGUu5wYACzx+XwP8CxwGxgLzgPs8zt+L\nCpBDwK9AY49zbuBBYAuQAIxxjld1ymvjkTYMSAFCnd/XA2udfIuA8zzSNgJ+APYBB4APgXOANCAL\nSAYOOWlrABOdtLuAFwHjca2LgPeccl4DmjvXeBjYD3zj5R5VBY4413cE2OocbwXMddr8D9DTI88X\nwDhgtpPnpPsPhAATgBjnfk53jkcA0cBTQDwQCwzwyNcDWAMkAlHAqx7nmjrtvAvY7VzTYI/zpwBf\nOvVtAp4Doj3ONwCmOfdvB/BYjrxfOHk3As965vVyfe877U8E1gNtgA7A3uz/iZPuRmCts98RWOnk\n2QuMcI5HOdeV7GwX+/g8DkSfxyTn/90MWOyU/y1QOdjvp93sJuJdTjvviwto7fz+AngdaAEc9Xgn\n/gS2OWlTnOe9MioPxzsyZI+Tt4JT1gBOlodV0Bm13c77Nw6o5qSPcMrITS6dgo5870Ll6QKPvJ2c\n9y4BlfNX5HIPvvK4hmTgGed4L0fmJAB/A+fkcR/dwGPAdlT+vZND3uQnMx52ZMaOgt5Tp4z7nfKT\nnDa3ze26vLQ9r+/gLnRWIfteF1aODwG+c/Inod+ui5xzg4Dvc7TpA+ADZ99v157X/9S51ufQ70Yq\nUNFp2x6n7Ehy0WnKyxb0BpTlDRXGXXM5NwBHaQZCUWXiBtRk5nEgA7jXOd8b2Aq0dM6/CCzyKMsN\n/AhURxXdfcC1zrnxwBseaR8BZjv7bVEh3AHtjd7ltLmy87KsQ4XxKagCe6mT7248FH7n2ERgOnAa\n0ATtANzrca2ZTt0V0FGaKcALzvkq2WXncq/cwFnOfmVUoD6PzpT8n/Myn+2c/wL9cFzi/K7qpbxZ\nTv01nDK6OMcjnHYOca7/OvQDWcM5fwVOBwQ4D/249XZ+N3Xa+Ylzr85HOxctnfPDUQFVAwhHhVKU\nc64CsAp4yWnPmeiH5xqPvPOAmkBDVNhG5XKvrkWV3+rO75ZAPWd/I9DNI+104L/O/hLgdmf/VI4r\nx02c6/IU0L48j9OB04HWQDowx7lH1Z123BXs99NudhPJfXADVWAfdPYnAK85+97eiRPKcJ7/cajs\nDAOWAQ845wZwsjx8H5jhvOOno/J8mJM+P7k0FvgLqO+U1wmVqeGoUt7NSXeV8zvUl/sAnI0OPHR1\n6n3Wee+9dnide/Kncw2N0G/Afc45X2TGb07eqoW4p31RxS5bCW2Go5Tn9v/1KDfX72DO/BRNjg9B\nFdFuTj3DgCUez9RR4HTnd0VUQe7oz2vP439ayTm/C1jtXFtV5/8VxfFvSGOcb3F53YLegLK8OQ9g\nMtqjy96yhcgAjivNd3kKEOdYFMeVzl+y953fFZwXrJHz242H0omO5A1y9rsC2zzOLQLucPbH4XwI\nPM5HApcDl6DKdwUv13Ws7c7viqhi5NljfQD42yP97hxlfIkqmOE+3EdPpbkLEJfj/GScUV9Uaf4i\nj7Lqo73vGl7ORaC9cs+PYXy24PKSfhTwnrPf1GlnA4/zy4BbnP3twNUe5+7DGaFAp4Fz3p8XgM89\n8l7jce5+chlpRjsR/zplVshx7jlgkrMf4jxDdZ3f81ChHpojT/Z1ed4TX57HSzzOrwSe9fg9Ang/\nmO+m3eyWvZG70ryE4x37CcDrzr63d8JTsaqLdpireZy/FfjL2T9BHqIK1BE8lBFU/u5w9nOVS867\nl4LHyKhHmkHAxBzHfiWXDmvO+wC8jMcMoNPOPeQ+Wu3OIacGAnOcfV9kRkRu7fHhnv6Gx6iuL/9f\nj/O5fQe7eGlHUeT4EOB3j3OtgRSP3wuAO539q3G+2/689jz+p5d7pB/gcb6586x1xc4OIiLWETDA\nCDoSWctjG+8lXQP0wfXE83cT4APHSSUBOOgcD/dIs9djPwUdrQA1YzjVGNPRGNMUuADttWaX+3R2\nuU7ZDVHFshEqAHyxBw5FR4B3exyLytG+nNEenkNf2OWOd/I9PtQDeq9ylrXbOQ56z/OKLNEINSlJ\nzOX8wRzXfOxeGmMuNsb8bYzZZ4w5jJrE1M6RP7f/Q8525/z/Nsjxf3gBqJNL3qjcLk5E/gbGoKNP\n8caYTxw7TYCvgZ7GmFOBW4D5IhLvnLsPHYXYbIxZbozpkVsd+PY8xnvsp3r5fToWS8mmIToNX1Ca\noPIwzuMd+RgdIczG830OQ2d3Vnmk/wWVq9nkJpdC0ZHq7bm0o28OudIZtcn2hfp4yBpRLSqa47LW\nGznlVHZaX2RGXnI7v3vaEO/3wBdy+w56u86iyHE4UQ6mANU8nKwno8owwG2ovM4u11/Xntv/1Ov/\nQUS2AU+iCn+8MWaKMaa+j3WVSawjYMkgFuiZ/cNxfmvocT4KHeWYUtCCRcRljPkOfRn3AT+JyFGP\nct8UkWE58xljLgEa5+KkIjl+H0CnD5ty3DmmMScKlBPyOMraA05dnYE5xph5IrIjn0uKBRoZY4zz\nwoMKlch88mUTDYQYY2rkoTjnxmRgNGr6kmGMeZ8TP2x5EYcq7NntbJSjTTtF5Ow88jbmxHubKyLy\nIfChMSYMtaF7FnhFRGKMMUtQW+Y7gI888mxDBTXGmJuAqcaYEE7+X0MRnkeLpTRgjOmAKkgLC5E9\nGp15q53HoIPne3UA7Ui2FpG4AtZ1AB2FbI6aCngSBXwlIg/4WFbOdz0WNUMDjn2XGqG+ILmRU05l\np/VFZniTNdnkd0+j0XtQ0HKz2+b1O+iFosjx/NoxFRhpjAlHTTU7eZTrr2v35X+a81s9BZjiDL58\nAryNzo6XS+xIc+DxxXN1NnCeMaa3MaYSauvmORrwMTDYGNMawBhTwxjTtwB1Tgb6o0rRZI/j/wMe\nckahjTHmNGNMD2PM6ahpQRww3BhzqjGmmjHmUidfPNDQCb2Eo1R/B7xpjDndGNME+C8wKdcGGtPX\nGJPdMTiMvqi+jGovRXvozxljKhtjIlAnjm9yufYTcD5KvwAfGWNqOmX4GoPzdCDBUZg7ovczP0GY\nzXfAC06d4cCjHnmXA8nGmOeMMacYYyoaY841xrT3krch6mzjFWNMe2dEvDJ6n9JQc5RsJqLTtuei\nTp7Z+e5wlGxQ+/rs/8d+528zjzIK+jzCif+X8urJbim5GABjTHVjzPWoz8NXIrLR87wvODLmd+A9\nY8wZxpgKxphmuckZRxH6HzAq+x00xoQbY67xoS438LlTV31HdlxijKmCyt+exphrnOPVjMZ4D8+l\nuHhOfM+/A3oYY6505MnTqDxZnEeTnnHkVCPUN+db53hhZIbndeZ3Tz9z6m7nfMuaG2OyBxdyXldO\n8voO5qQocjy/b9N+dGb4C9Q0598AXHuB/qfGmLOdtFVRxT3n96TcYZXmwPOTOTFO8zTnuDgbInIA\nNeZ/Bx05aIXagaY752egvbtvjDGJwAbU4QuPssjx+9gxEVmO2szVRxXG7OOrUPvYMeg05FacHqQj\njHuiPdgotDd7i5P1T9SZa68xZp9z7DHURm0Hapv1NWoHeFJ7HNoDS40xycBM4HER2eXtBua4lkyn\nXdehCt0Y1A5sSx515eROdGQ8EhUqj3urywsPA68ZY5JQ27Bvc5zPK+9r6Mj7TlQAfo86e2Z3Oq4H\nLkTv337gU9RpDmAoaoKyE7VJnJhHXdWdvIdQm/oDwLse56ejI0DTRSTN4/i1wD/O/+N9oL+IpItI\nCvAmsMjo1GDHQjyPOY/58j+yWIqTn5z3OgqdUh8JeJqM5Xxm83t+70Kd8bKjRXzP8YEQb8//INTB\neanzTv2Bmkv5Ut8z6Du4AjV7eAu1f96DOuANRmcZo1AlKbfv/lvAS857/pQjU+9AoybtR6MH9ZS8\nw5DNRJ3h1gA/owp9Yb5h3sj1norIVFROTUYdw38Aanm7rpyF5vEd9Namoshxb//3nL8no/bDk3Mc\n98u1F+J/WtUpYz86iBaKvh/lluyQYIXLrL3JiajNjgCfishoL+lGo0pOCmpkvqbQlZYDjNo4RQO3\nici8YLfH4n+MMQNRJ8H/C0Ld21DPaxs/vJRgjOmGOp5WBD4Tkbe9pPEqZ3PLa4x5F/3IZ6A2kfdk\nmywZY15AQ4S50A7t74G9QktpxxjjBpr7YGJXZgimHLcEh6KONGeiIavaoPY3jxhjWnkmMMZ0R1+k\nFqgN67gi1lkmcabQajrTIIOdw0uD2SaL/zDG1DO6LGwFY0xLNObq9PzyBaAdNwFuqzCXHowxFdFR\nsG6ox/2tvsrZfPL+joZQvACNj/uCk6c10M9J3w01ZbKzkpZyT0mR45bgUSRBKCJ7RWSts38EdQDI\n6XHaCw0vhogsA2oaY+oWpd4yyiXo9Fz2lMkNIpIe3CZZ/EgV1K4vCTVvmYGHI15xYIyZi0bVeKQ4\n67UUmY5o+KldjnnSN+i0uyfe5Gy9vPKKyB8ejkXLOO583BuYIiKZjsnUNqcciyUvyoPJVdDluCW4\n+C16htFwZm1R4etJOCeHaGnIiaFXyj0iMhS1XbWUQUQkCg+v5SC1ISKY9VsKjTcZerEPacLxHiIr\nZ15QU4zsyAYNOHGWK7ssiyVXRKRisNsQaEqCHLcEF78ozY6X6VTgCWfE+aQkOX6f1CM1xpSHXqrF\nYimjiEigooL4KhsLVb8x5kUgQ0RyOh/l2QYrsy0WS2mmMDK7yHZqTtiSaehKYzO8JInhxFiGDckl\nzqOvK7KUle3VV18Nehvs9dprttdc9C3A5JShjTh5MSRvcnZPfnmNMQOA7sDt+ZRlZXY5fbbtNZf9\nrbxdr0jhZXaRlGZjjAHGA5tEZFQuyX7ECWNmjOkEHJbjq5BZLBaLJW9WAi2MMU2d+Lv9ULnqSW5y\nNte8TlSNZ9FVS9NylNXfGFPFGHMm0AKNQWuxWCzlmqKaZ3RGY/6tN8Zkh5EbjLNimYh8IiKzjTHd\nnTBXRzkx9qXFYrFY8kBEsowxjwK/oWHjxovIZmPMg875XOVsbnmdoj9EHZv+0PEPlojIwyKyyegq\nopuALOBhKcrQjMVisZQRiqQ0i8hCfBitFpFHi1JPWSUiIiLYTShWytv1gr1mi38QkV/wWJjIOfZJ\njt9e5ay3vM7xFnnUNwzwZVnhckV5fLbtNZd9ytv1FoUiLW7iT4wxdjDDYrGUSowxSOAcAUskVmZb\nLJbSSmFltg1Yb7FYLBaLxWKx5INVmi0Wi8VisVgslnywSrPFYrFYLBaLxZIPVmm2WCwWi8VisVjy\nwSrNFovFYrFYLBZLPlil2WKxWCwWi8ViyQerNFssFovFYrFYLPlglWaLxWKxWCwWiyUfrNJssVgs\nFovFYrHkg1WaLRaLxWKxWCyWfLBKs8VisVgsFovFkg9FVpqNMZ8bY+KNMRtyOR9hjEk0xqxxtpeK\nWqfFYrGUJ4wx3YwxkcaYrcaYQbmkGe2cX2eMaZtfXmNMX2PMRmOMyxjTzuN4U2NMqofM/iiwV2ex\nWCylg0p+KGMC8CEwMY8080Sklx/qslgslnKFMaYiMAa4CogBVhhjfhSRzR5pugPNRaSFMeZiYBzQ\nKZ+8G4A+wCdeqt0mIm29HLdYLJZyS5FHmkVkAZCQTzJT1HosFoulnNIRVWJ3iUgm8A3QO0eaXsCX\nACKyDKhpjKmXV14RiRSRLcV1ERaLxVLaKQ6bZgEuMcasNcbMNsa0LoY6LRaLpawQDkR7/N7jHPMl\nTQMf8nrjTGPMamPMXGPMZQVvssVisZQ9/GGekR+rgcYikmKMuQ6YAZztLeGQIUOO7UdERBAREVEM\nzbNYLJaCMXfuXObOnVtc1YmP6fw1oxcLNBKRBMfWeYYxpo2IJOdMaGW2xWIpDfhLZhsRX+VxHoUY\n0xT4SUTO8yHtTuAiETmU47j4oy0Wi8VS3BhjEJGAmKEZYzoBQ0Skm/P7BcAtIm97pPkYmCsi3zi/\nI4ErgDN9yPs38LSIrM6lfq/nrcy2WCyllcLK7ICbZxhj6hpjjLPfEVXUD+WTzWKxWCzKSqCFE9Wi\nCtAP+DFHmh+Bu+CYkn1YROJ9zAseo9TGmFDHgRBjzFlAC2CHn6/JYrFYSh1FNs8wxkxBRzRCjTHR\nwKtAZQAR+QS4GRhojMkCUoD+Ra0zkGS43cw/fJhVR46QPYpSs1IluoWE0PSUU07OEB8Pq1bBpk2Q\nkgLp6SACTZpAixbQqhXUr1+otuw/up8dCTuISoxiT9IeEtMTSclMISUzhUoVKnFq5VM5tfKphJ0a\nRuMajWlUoxHNQ5pTrVK1otyCEk9KChw8CAkJuiUnw9GjcOQIpKbqvyA9HbKywO0GlwuMgQoVdKtc\nGapW1e3UU+H003WrXh1q1dItJETP+xu3CCuTk5l/+DBZzvNVrUIFutaqxbmnnYbTvyx7iMCOHbBl\ni26xsfqPqFZNb/z550O7drpvOQERyTLGPAr8BlQExovIZmPMg875T0RktjGmuzFmG3AUuCevvADG\nmD7AaCAUmGWMWSMi16HyfKgxJhNwAw+KyOFivWiLxWIpgfjFPMMfBHuqb3lSEu9FR/NbQgLnnHoq\nnatXp5KjwMRlZDD70CEaVKlCv7Awnjx0iFMnToQfflCN7aKL4LzzVPPK1rR27oStW2HjRmjUCPr0\ngRtv1HReyHBlsHTPUv7c8Scr41ayJm4NqVmpNA9pTqPqjWhYvSG1qtXi1MqnckrlU3CLm5TMFI5m\nHCX+aDzRSdHsPryb3Ym7aR7SnLb12nJpo0u56qyraFarWalQxtxu1aV27YLdu3WLidEtNhb27dNN\nBGrXhpo1VcGtXl1v/WmnwSmnHFeIK1c+rihnl+9yQWamKtVpaapkHzmiW1LScUX80CEtKywM6taF\n8HBo0AAaNtT+UJMmcOaZEBqqCnl+7ExNZXhUFD8ePEhIpUp0rVWLU52GJbpc/HpIJ1/6hIYyqHFj\n6lapEqC7XIy43bBgAUyfDjNm6M1v1Uo7kw0baq8mPV1v9tq1sH49NG0K/fvDnXfqTS4lBNI8o6QS\nbJltsVgshaWwMrvcK80HMzMZvGMHPx48yEtNmnBTaCj1vAwxulwulk6fzodRUSxr1IgP4uLodc01\n0LJl3lqTywWLF6vS8N13mv6ll+CKK0hMT2JG5Aymbp7KvF3zaBnakq5nduXi8ItpV78djWs0LrCy\nm56Vzj/7/mF13GoWRC1gzo45VK1UlR4tetCvTT86N+5MBRPchSCPHNG+xKZNEBmp27ZtOhBZs6Yq\no9mKacOGqqzWrw/16qkSe9ppvimqRUEEEhNVSd+7F+LiVHmPjj6u0O/cqXpf8+aqB55zjuqErVrp\nftWqkO52825UFKP27OHR8HBur1uXFqee6qU+4Z+jR/l8714mxccztGlTHmzQgIqloLNzEllZ+qy/\n+SZUrAh9+2qnsU2bvP9xWVmwciV89RV8+y1ceCG8+ip06VJ8bS8kVmm2WCyW0oNVmgvBrIMHuTcy\nklvq1OH1pk2pWbmy94R//AGDBkGlSvDGG/x50UU8sm0brU87jS/POYczKvlo5ZKZiXvSV6S+9gpR\nVVO5/+o0QjtfzS1tbqFb826EnBLiv4tzEBE27d/EjMgZfLvxWw6lHuL2827nofYPcWatM/1eX05i\nYmD1alizRrf161UBPeccaN36uIJ59tlw1lmqEJcmDh2C7dvV4iAyEjZv1s7Azp0Q3uUI+x7eSGM5\njRdrNKdb22rUqpV/mf8cOcLDW7eS6nYzrU0bGlcrReY2M2fCM89oD+fll+HqqwvXw0lPhylTYMgQ\nnZ156y0491y/N9dfWKXZYrFYSg9WaS4gn8fFMXjHDqadey6da9TwnighAR55BJYvh+HD4aabjikA\n6W43j23dyurkZGadf36+0+lJ6Ul8uupTxq0cR63K1RkZdz5dPp5Nhf/cD6+8oradxcCm/ZsYv3o8\nX677kosbXszjHR/nmmbX+MV8IzUVVqzQgfVly/S2ZWSo9Uq7dtC2rZquNm+uA5BlmT/2Habf5o30\n3d+cagvqsnq1WiDUrw8XX6xb586qD3rrc4kI7+3Zw6g9e5h93nmcd/rpxX8RBSE+Hh5/XHtI48bB\nVVf5p9y0NC3vrbfg4YfhxRfV7qaEYZVmi8ViKT1YpdlHRIRhUVF8FhfHr+efT0svU+UA/P033H03\n9O4Nb7+tHmNeyhq6axeT4uP57YILaObFUTAhNYHRy0YzZsUYrjrrKp68+Ek6hndUJXXvXnjsMdiw\nQaezzz/f35ebK6mZqXzzzzeMWDKC0yqfxkuXv0TPs3sWSHlOToZFi2DuXJg3T0eR27SBSy+FTp2g\nY0c1tSiNFgZFYeq+fTy8dSvftG7NlR5Dyy6XjkIvXw5Ll+q927NH79UVV0BEBHToAJ79r2/i43li\n2za+b9OGy2vWLP6L8YWZM+GBB+Cee9ScwpvDbFGJi9PyExJg0iS1hylBWKXZYrFYSg9WafaRV3bu\nZOaBA/x6/vnU9xYeQQSGDYOPPoLx46Fbt3zL/CQ2ltd27WJR27bHImykZqYyaukoRi4ZSa+WvXj+\nsuc5u7bXNV1g8mR48kkYO1btP4sRt7iZETmDN+a/gTGGd69+lyvPvNJr2sxMVfbmzFGLlfXroX17\nVfYuv1xHT0ubeYW/+X7fPp7cto1Z553HhWeckW/6gwdVeZ43TzsfW7bAZZfpQO3VV+tI9F+HE7h1\n0yamljTF2e2G117T92TqVH0AAomIviNDh8Knn6qddAnBKs0Wi8VSerBKsw9M3LuXV3ftYlm7dtTx\nZk6RkaEjZhs2wE8/qQeaj4yKjuazuDgWtL2Qnzd9y4t/vUjH8I681fUtWtT2YVRszRpVAm67Dd54\n43jIh2JCRPh+0/c8P+d5Woe15t2r36VVWCtiY2H2bPjlF/jzT7U7vvpq3Tp3DsygYmllWVIS12/Y\nwJwLLuCCQppTJCToJMcff8Dvv6vJS7duEN7rEJ+GbmZRu7Y0z212pDhJToa77oL9+1Vhrlev+Ope\nuVJngJ56SrcSMJVhlWaLxWIpPVilOR/mHz7MzRs38veFF9LG23Do4cOqtNaoAV9/XeAhUxHh1vXL\n+T1mPWfFfMYH175H58adC9bI/fs1LF3TpjBhgndj1wCTlpnOSz+OZdw/b1F96wOk/fYS13Y9he7d\n4Zprilc3Kk3sSk3l0jVr+PTss7k+NNRv5W7dqh2W2bNh7hmxVOwXzcvx7bitZ2UaN/ZbNQXj0CHV\n5M8/X2dkghEeLyoKevTQYfkxY4JuJG+VZovFYik9WKU5D7anptJ59WomtmrFNSFeIlQcPqxDp506\nwahRBf4Ap2Wl8cb8N/h41f8Iu3g8V9ZrwdizWxausSkpqryfcYaabRSDQuJ2q+PetGkaUlcEruoT\ny5azniQ6axUfX/8xVze7OuDtKK0kZ2Vx6Zo13FevHk82ahSweo4cgdsWbmNFwhEy/3s+TcIrcOON\ncPPNGsmwWNi3T9+Va66Bd94J7ihvUpK+K/Xrw5dfBlVxtkqzxWKxlB6s0pwLmW43ndes4Y66dXm8\nYcOTE2QrzJ07w/vvF1gJWLt3LXdOv5MWIS0Y030Mp51Shw6rVvHWWWdxU1hY4Rqdng79+mnc2qlT\nAxJZQ0QjXXz7LXz/vQ6s9+2rA90XXHD8Nvyy9RcemvUQPVr04N2r3+W0KuXcaNkLAzZvpqIxfNay\nZcAXkXGJ0GP9etqfXp2uUWfyww/a2QkLg1tu0XVBmjULUOWxsdC1qz6br75aIswiSEmBnj119ZkJ\nE4KmOFul2WKxWEoPVmnOhSE7d7IsOZnZ5513skKTmHh8hPmDDwqkBLjcLoYvHM4Hyz7gvWvf4/bz\nbj9W/rKkJHpt2MDa9u29Oxv6Qmam2je73RpZw0/KQGSkDmBPnqxm0/37q7KV17oTh9MO88SvT7A4\nejETb5jIJY0u8UtbygLT9u/n+R07WHPRRZxeTOY0cenptF25khnnnkunGjVwudSZ8LvvtAPUtCnc\neqv+b/1mTpOQoIuM3HYbDB7sp0L9REoKXH+9rrw5YUKx+wOAVZotFoulNGGVZi8sS0qi94YNrPGm\nvGZkqF1my5Zql1kAhTk2OZbbpt1GxQoV+aL3FzSqcfKUfJ7Kuq+kp6vdZvPmGqu2kOUcOADffAMT\nJ+qKdrfeqrrPRRcVrMjpm6czcNZAnuz0JM91fi7oKwsGm7j0dC5cuZKZjvJanPywfz+DvCjrWVnw\n119qlv/jj3DJJeqv17t3EZw2U1LUHOPii2HEiJIxwpyTo0fhuus0Zt/IkcVevVWaLRaLpfRgleYc\nHHW5uHBpKAM/AAAgAElEQVTlSoZ7M5MQ0RjMSUk6t12AUdzft//O3TPu5uH2DzO4y2AqVvCeN9ss\nZEC9ejwcHl74C0lO1phuvXrplLiPuFwafeHzzzUSQ/fueslduxbNvzA6MZr+0/pTvWp1Jt4wkbDT\nCmmCUsoREbpv2ECHM87gtTMDv7KiNwZs3ky1ChX4OBeD5qNHdfX2iRPVFKdfP7j3Xg0T6LPem5Wl\nNjs1aqjdcBBGcX0mIUHNrAYO1PjnxYhVmi0Wi6X0UFiZXeQvoDHmc2NMvDFmQx5pRhtjthpj1hlj\n2ha1Tl94cccOLqle3btd8auvwr//qo2CjwqziPDm/De5Z+Y9TL5xMi9f8XKuCjNA5QoV+KpVK17Z\nuZNdqamFvQx1CPzlF13QYcKEfJNHRenlNWmif6+6Cnbt0ku99tqiB+RoVKMRc++eywV1L+CiTy9i\nZezKohVYSpmwdy/7MzJ4uUmToLXhgxYt+PXQIeYcOuT1/Gmnwe23w2+/wbp1avbbv7/arI8Zo+b8\n+fL442oq9PnnJVthBqhVS8OMDB+uC66UIYwx3YwxkY4cHZRLGq9yNre8xpi+xpiNxhiXMaZdjrJe\ncNJHGmOuCdyVWSwWSylCRIq0AV2AtsCGXM53B2Y7+xcDS3NJJ/5iXXKyhC1cKPvT008++eWXImed\nJRIf73N5yenJctO3N0mnzzpJTFJMgdry+s6dcsOGDQXK45XNm0XCwkSWLDnpVFaWyE8/iXTvLhIS\nIvLooyLr1hW9yvyYtmmahL4TKl+t+yrwlZUgDmZkSN2FC2VlUlKwmyIz9++XlkuXSrrL5VN6l0tk\nzhyRfv1EatQQuftukaVLRdxuL4k//likVSuRxES/tjngrFghEhoqsmpVsVXpyK8iy1NvG1AR2AY0\nBSoDa4FW4oOczSsvcA5wNvA30M6jrNZOuspOvm1ABS/tCvBdtVgslsBQWJld5KEjEVkAJOSRpBfw\npZN2GVDTGFO3qPXm0R4e2bqV15o2JTRnuLY1a+Dpp3Xhkjp1fCpv1+FdXDL+EmpUrcHcu+fS4Azf\nFzwBeKZRIzYcOcKvBw8WKN9JnHOOjvbdfLNGMUDDOg8bpguOvP66Rr+IjoYPPyyeFblvbHUjf9/9\nN0PmDuGZ35/B5XYFvtISwMs7d9InLIyLfFjxL9D0rF2bZqecwqg9e3xKX6GCmuh8843GgG7T5rh9\n+/jxupgKAAsXwiuv6Iht9eqBu4BA0L69+gDcdJMuuVj66QhsE5FdIpIJfAP0zpHGm5ytl1deEYkU\nkS1e6usNTBGRTBHZhSrNHQNwXRaLxVKqKI751nAg2uP3HsBL7Df/8HV8PCkuF/fnXM3v4EG1zfzo\nI2jd2qeylscs59Lxl/Kftv/hs16fUbVSwSNhVKtYkdEtWvD4tm2ku90Fzn8C118PAwdy9Jo+3H9n\nGmefDdu3ww8/aJzlAQOguBeLO7fOuSy/fzmr4lZx03c3cTTjaPE2oJhZk5zM1P37eTNIdsw5Mcbw\nQfPmvBMVxZ60tALlDQuDZ59V5XnYMLV/btwY3nwomqybblEb5hY+rGZZErn5Zu1F3nabGviXbrzJ\n0JyOErmlaeBD3pw0cNIVJI/FYrGUeYprybmcxtZevUeGDBlybD8iIoKIiIgCVZKUlcWgHTuY1qYN\nFT09nVwuNe686Sb9kPrA9M3TeeDnB/i81+f0bNmzQO3ISffatfkkNpaR0dEMLqQNrMulg36jfh3M\n8zvWMPCMpxi+7SNq1y5S0/xCyCkh/HbHb9z/0/1EfBnBT7f+RL3Ty97SgW5nFuPNM88kpHLlYDfn\nGM1PPZWB4eE8vX0737ZpU+D8FSpoIJlu3WB7ZCYVr7yFN5Of4J/x3XjyDLj00pIZMCNfhg3TqB+v\nvqpL0/uRuXPnMnfuXL+WmQe+etsF8r8UEJltsVgsxYHfZHZhbDpybqjdW242zR8D/T1+RwJ1vaQr\nso3Ks9u2yb2bN598YuhQkYgIkcxMn8r5cNmHEj4yXFbGrCxym7LZkZIiIQsWyJ60tALlS0oSGTVK\npGlTkUsvFfnuO5HMA4dFzjxTZOpUv7XPH7jdbnlt7mvSdFRT+ffAv8Fujt/5Ki5OOqxcKS6vBsDB\n5WhWljRZvFjmJiQUraBBg0Suu06SDrtk9GiRZs1EOnQQ+fZbn1+fksW+fSKNG4vMmhXQagisTXMn\n4FeP3y8Ag3Kk8Spnfcyb06b5eeB5j9+/Ahd7aVeA7qbFYrEElsLK7OJQmj0dVDoRIEfA6NRUCVmw\nQGJyKqULF4rUrSsSk78Dn9vtllf+ekVajG4hOxN2Fqk93nh22zZ5MDLSp7SxsSLPPy9Su7bIzTd7\n8f9btkwdA3f6v51FZfzq8VJvRD2/djqCTZrLJU2XLJH5RVVKA8hXcXHSadUqcRdWqf/1V5HwcFU0\nHVwukRkzRC67TDtuo0aJHDnipwYXF/Pni9SrJxIXF7AqAqw0VwK2O3K2Cvk7Ah6Tsz7m/Ru4yON3\ntiNgFeBMJ7/x0q6A3U+LxWIJJIWV2UWO02yMmQJcAYQC8cCrqNc1IvKJk2YM0A04CtwjIqu9lCNF\nact/IiOpU6UKw8466/jBxES48EJd7a9Xrzzzu9wuHvvlMZbFLOOX23+hzmm+OQoWhEOZmZy9bBmL\n27Xj7FyMj7duhXff1dWzb7sNnnpKHf28MmKExpmePx9KkLkAwMzImdz/0/18c/M3XHnmlcFuTpEZ\nvWcPvx06xKzi8LAsJG4R2q5cyWtnnknv0NCCZY6Lg3btNDbh//2f1yTLlumzOX/+8VDIBa0maLzy\nigarnjUrIKHzAh2n2RhzHTAKjYYxXkTeMsY8CPnLWW95neN9gNGo7E4E1ojIdc65wcC9QBbwhIj8\n5qVNRZLZwcbthk2bYMkSfbZ37dLXIC5OzeEqVtQQnaGhuthko0bQqhW0batbrVrBvgKLpfCkpcGG\nDRofYcMGDVcbHa3Pf0aGvgMi6vvSoIGGLL3oIl0w66KLoFq1YF9B0SjXi5tEHj1Kl7Vr2dqxIzWz\nlUcR1TpDQmDs2DzzZ7mzuHvG3cQmxzKz/0yqVw1ctIBhu3ez7siRk2xP166Ft97S1dwefhgefVQf\n1jxxu9U58KKLNHxGCWPernn0/b4v43uNL7JdeDBJzsqixbJl/HbBBVxw+unBbk6ezDp4kOe2b2d9\nhw4n2vXnhYiuptexI7z2Wr7Jt2zR/trUqbpgztNPQ8OAufb6iawsuPxyXTP+ySf9Xrxd3KR04Hbr\nkvNff63Pb61aqgR06qQLr9avr0vPV6mij0xWlkYpio5WpWLjRlUy1q7VZ/7//k+3rl31U2OxlFTS\n02HBAvj7b93WrlU/73btNNpW06baMaxfH6pW1Q6jiD7/sbH6/K9YoZ3MyEhdA+K226BnzyKsNhtE\nCi2zCzM8HYiNIkz13bRhg7y9e/eJBydOFGndWiQlJc+86VnpcuO3N0q3Sd0kJSPvtP7gSFaW1F+0\n6FiM38WLNb5ygwYiI0eKJCcXsMDYWJE6ddRcowSybM8yqfNuHfl+4/fBbkqhGbpzp9y2cWOwm+ET\nbrdbLlu9Wr4oiCnCuHEi7duLZGQUqK6YGJGnnxapVUvkP/8R2batgI0tbnbs0PjNa9f6vWgCaJ5R\nUreiyOziJiVFZMwYNTE691yRt94qmmVbVpaGAR8xQqRHD5Hq1UW6dBEZPlxD6lssJYF9+0TGjxfp\n00ef0U6dRF58UWP1Hz1a+HIPHxaZMEHkqqtEatYUeeYZVUVKE4WV2UEXvMcaUkgBvDwxUcIXLZKj\nWVnHD+7Zox/H1avzzJuamSo9vu4hvaf0lrTMgjnoFYWxe/ZIh7/XSteuIk2aiHz0kUhqahEK/OYb\nkZYt8+0gBIs1cWuk3oh6MmndpGA3pcDsS0+XkAULZFsJvbfeWJCQIE0WL5Y0XxY82bZNDeeL0Ck4\ncEDk5Ze1mDvvFPHRbD84fP65yIUXFriDkB9WaS6ZpKWpYlu/vkjPnl7XhvILqakiv/wi8vDDOgDS\nurXIK6+I+GNdK4ulIMTFaQcxIkIXsOrbV8cQ9+8PTH27d4s89pgOngwcGFDXEb9SbpXmbuvWybg9\ne44fcLu16//KK3nmS81MlW6Tuknf7/pKRpZ/P6B5MXeuyBVdXVLpuyUy6OsE/327+/UTefJJPxXm\nf/6J/0cajGxQ6lYPfG7bNnno39IXCaTHunUyxvO98EZWlnr4jRzplzoTEkRef137q7ffXkJH3Nxu\nndoZMsSvxVqlueTx668iLVqoslwcK6Rm43LpDOJTT4k0aiTSpo3Ia6+JbN1afG2wlC8OHhT59FOR\nK6/Ukd/bb1cH7uIc64mP15nH0FB1GC/p0ZbKpdK8IjFRGuYcUfviC5ELLhDxtoS2Q1pmmnT/urv0\n/a6vZLqK5z87f77I//2fruA9YYLIuKgYudaf08QHDugQx99/+69MP7Np3yapP6K+fL3+62A3xScO\nZGRIrQULZHeRpgGCw9LERGm8eHHey2uPHCly+eX6lfcjiYkib7xxXHkucX2OPXs08syaNX4r0irN\nJYf4eJGbblJZ+9NPwW2LyyWyaJGOxNWtK3LxxSKjR2sbLZaikJqqIWh791bTi759RX74oYiz1n5g\n0yaRrl1FzjtPZMWK4LYlL8ql0nzDhg0yOjr6+AEfPoZpmWly/eTr5aZvbyqWEeYlS0Suvlpt6caP\nPz4rnOZyScPFi2V5YqL/Kps5U6R58xJrpiGiI871RtSTKRumBLsp+fLKjh1yX4kcLvWNq9eulfG5\nGZpt3672FAE0RE5M1JHn2rVFBgzQKksMEybk27kuCFZpLhnMmqWmGIMGlTwxmJmpo9933KHT5r17\ni0yf7rdH0FIOcLtFli4VeeghkZAQVU4nTFBZW5Jwu0UmT1Z17M03dVKzpFHulOZ1yclSb9EiScn+\nb7jdIr165WmWkenKlD7f9JEbvrkh4Arz2rUi11+v03OffOJdMI6OjpZe69f7t+KbbxZ54QX/luln\nNsRvkHoj6sm0TdOC3ZRcOZyZKbUXLJCtRfGWCDLzEhKk2ZIlkplzJNntVg+Od94plnYkJOhrGRKi\nwj4/q5FiIdtM4/XX/VKcVZqDS2qqyCOP6Do2c+cGuzX5k5Sk5vVduugI9PPPq5+qxeKNw4fVTvnc\nc3Vc7I03RKKigt2q/ImKUtvqLl1KXnvLndLc759/5B3PiBnTpomcc456fnjB5XbJHT/cId0mdQuo\n09+WLSL9++taCqNG5T1VkpKVJfUWLZK1BQ6ZkQdxcdq9C0CEAH+yKnaV1Hm3jvyy9ZdgN8Urw3bt\nkjs2bQp2M4pMl9WrZdLevSce/OILkbZti93obP9+kWefVYeRp59Wi6KgsmuXDoP7wX7EKs3BIypK\npGNHNckowWsP5crmzSL//a8+it27i/z1l/bpLJYtW9S5rmZNkVtuEfnzz9L3bGRliQwbpjpRSbIe\nLVdKc+TRoxK6cKEkZX/0ExN1JbN587ymd7vd8uBPD8rlEy6XoxmBGTmMiRF58EG143zzTd9Dx727\ne7fc8s8//m3MZ59pCLESbom/KGqRhL0TJnN3zg12U07gSFaW1Fm4UDaWuqXvTua3gwel1bJlx5f+\n3rtXQxSuWhW0NsXG6oegdm0d6PVnn7HAjByp3jNF/BJZpTk4zJ2rH+Phw0ufMpGTlBSR//1Px37a\ntdOgSH52N7CUEpYv1zBxYWEamai0RKTIi99/10/P+++XjHe1XCnNAzZvlqGeQTYffVTkvvtyTf/c\n789Jh087SGKa/w1/EhJ0ai0kREfRDh4sWP7kzEwJW7hQIv1pBuB265zIe+/5r8wAMWf7HAl7J0xW\nxJQcj4H3o6LkpjISK8rtdkuHlStlWvbS2Lfdpg9qCWDrVpFbb1WlZ8wYv0eB843MTB11//LLIhVj\nlebi53//04/wb78FtRl+x+US+fFHkQ4d1Ox+1qySoWRYAs/Gjaosh4eLfPihSBkYtzmBHTv0mb7n\nnuDb8pcbpTkmLU1qLlggB7O/sMuW6Vc3F231nYXvSKsxreTAUf/OBaemavzPsDDV1z39EQvKyzt2\nyIP+Dm4bGalDeTEx/i03AEzfPF3qjagnkfuDH+A30+WSxv520Awy38fHyyWrVuncXpMmJU4Sr14t\ncs01aqv37bdBUBBWrFDtqwiBTK3SXHy4XDpQ0axZCYzM4kfcbo2G0KqVyBVXlNAQjha/kJCgM9Vh\nYepqUtKcWP1JcrL6e115ZXDNqQorsysUeAnBIPNhTAx31K1LSOXKujj6Qw/pmr5e1jCdsGYCY1eM\n5fc7f6f2qbX9Ur/bDZMnQ6tWMG+eLkf52WdFW0b4kfBwvt2/n30ZGX5pIwAtW8KDD+oaxyWcG865\ngWFXDuPaSdcSkxQT1LZM3b+fptWq0aF64JZSL276hIWxNz2dxSNHwgcfwGmnBbtJJ9C2Lfz2G4wb\nB2+/rUsaL1hQjA1o3x7694fnny/GSi2FIS1Nl+6dPx+WLoWzzw52iwKHMdCnD2zYADfdBJddBkOH\n6nLIlrKBiC7n3qaN/r+3bIFnny2dy1L7yumnw4wZ0Lo1dO4Mu3cHu0UFpDCadiA2fBi1SHYiGhxb\nnW3sWO2Cexmamhk5U+qNqCf/HvDfUMTff4tcdJFOm/nbQ/v+yEh51d/u00eP6sjinDn+LTdAvLPw\nHWk9trUcTCmgjYufcLvd0m7FCvkxUEsnBZEPP/9c+nzySYmf53W5RCZN0se2d+9iXF3w8GGNVbZ0\naaGyY0eaA87hw2p11rdv8GPRBoOoKA0Q1aaNxsK1lG4SEzXYVatWIgsWBLs1weH990UaNhTxdxAx\nXyiszPaH4OwGRAJbgUFezkcAicAaZ3spl3LyvcgPoqPl5mynuX37dC7Di+3p4qjFfrWT3bxZhVXT\npiJTpgTGOWPzkSNSZ+HC4yH0/MXMmbrEdi5RRUoaz/z2jHQe31lSMop/furvQ4ek5dKlx53mygq7\ndsmRBg0kdN482VJKQuilpuo0Ze3aGkos2yQ7oEycqL3iQryDgVaa85OzTprRzvl1QNv88gIhwB/A\nFuB3oKZzvCmQ6iGzP8qlvgLfp8ISG6u2kI88UjJjvhYXbrf6eYeGqqOgpXSyfr2uVvnQQ+WzA+jJ\n5MlqHTd/fvHWGxSlGagIbHOEbGVgLdAqR5oI4EcfysrzAjNdLmm6ZIkszbY1vfdejdOTg837N0vd\nd+v6JZTZgQPqYxgaKvLuu4F/uHuuXy8fB8IGuWdPjflSCnC5XXLr1Fulzzd9JMtVvF/HHuvWyael\nwAa8wPTpI/Laa/Li9u0ysJQZge7fr6up1a6tSnRA+35uty4rPm5cgbMGUmn2Uc52B2Y7+xcDS/PL\nC7wDPOfsDwKGy3GleYMP7SrwfSoM27bp6n6vv17iJ0qKjdWr9Z48+WT57kSURqZMUZ1i4sRgt6Tk\n8PvvOgY6c2bx1RkspfkS4FeP388Dz+dIEwH85ENZeV7gt/Hxctnq1fpjyRJdMjqHs1ZsUqw0HdVU\nJqyZUKCbl5P0dI1EFRqqSnNxzdbPS0iQFoEY6dy2TbWOErGqRP6kZabJlV9eKQ///LC4i+krufHI\nEam7cKGklrUv0B9/6Nc1NVXiHCfafcF2Wy4EkZE623PmmSLffx9A5WntWpXeBXzpA6w0+yJnPwb6\nefyOBOrllddJU9fZrwdESglTmtetU1H/yScBr6rUceiQmqv07x+kyDOWAjN2rJojrFsX7JaUPFas\n0IV+iqszUViZXVRHwHAg2uP3HueYJwJcYoxZa4yZbYxpXZiK3ouO5qmGDdUT77HH1GPIw1krOT2Z\n7pO7c1/b+xhw4YDCVIEI/PijGuXPmaPOJh9+CKGhhSquwHSpUYMalSox6+BB/xbcrBk88ECpcXSq\nWqkqP9zyAwujF/Lu4neLpc5Re/YwMDycahUrFkt9xUJmJjzxBIwcCdWqUa9qVW4MDeWT2Nhgt6zA\ntGwJM2fC//4Hr78OERGwZk0AKrrgAnUKfOmlABReaHyRs7mlaZBH3roiEu/sxwN1PdKdaYxZbYyZ\na4y5rIjtLxSLF8PVV8P776v4spxIrVowezYcPaoOg6mpwW6RJTdE4K234L33VK84//xgt6jk0b69\nBlZ48UUYPTrYrcmdSkXMLz6kWQ00FpEUY8x1wAzAq8/zkCFDju1HREQQEREBwLKkJOIzM+kVGgpf\nfgmVK8Pttx9Lm+XOov+0/rSv354Xu7xYqAv55x/4738hJkYV5W7dClVMkTDG8ER4OKNjYujpb019\n8GA45xz9El16qX/LDgA1qtVg1m2zuHT8pTSp0YR+5/YLWF0HMzP5fv9+Ijt2DFgdQWHcOAgPh969\njx16vGFDuq9fz6DGjalcodQFz6FrV1i9GsaPh+uugx494M03oV49P1YydKi+KwMHqhLthblz5zJ3\n7lw/VponvshZAONjmpPKExExxmQfjwUaiUiCMaYdMMMY00ZEknPmy01mF5Xff4c77oCJE4Mji0sL\np5wC06bBgAHQvTv88gtUqxbsVlly8sorMH26RgWqXz/YrSm5tGql9+jqq+HQIXj1VY0q4g/8JrML\nMzwtx6fnOnHi1N8L5OKk4pFmJxDi5Xiuw+i3b9woI6Ki1Byjfn1dLsche7W/a7+6VjKyCj5HdfCg\nOpeEhYmMHh38aa40l0vqLVoUmNXoJk1SR6dStMzU2ri1EvZOmCzYHTj34rd375a7ypo7+r59al+0\nceNJp65YvVqm5FxauxRy+LDIM88EyN75o4907ttHOxACa56Rr5xFzTP6e/yOREeOc83rpKnn7NfH\nMc/wUv/fQDsvxwtzZ/Plu+/UMWjhwoAUXyZxuUT69dOlxEuReC8XjBkjcvbZxWfmWRbYu1fkwgtF\nHn88cM9zYWV2UYV5JWA7agNXBe8OKnUB4+x3BHblUpbXC4tNS5NaCxZIQkaGyKBBInfffcL5txe+\nLeePO7/Aq/1lZqp9UViYKs0H/Lv2SZF4dccOeSgQDltut8ill6r7dSni162/St136/o1fGA22YuZ\nrExK8nvZQeXBB0WeeMLrqR/27dPFTsoI//4r0qOHLo7y009+KjQzU+S880SmTvUpeYCVZl/krKcj\nYCeOOwLmmhd1BMxWoJ/nuCNgKFDR2T8LNemo6aVdRbjB3vn0U7VhXrvW70WXedLSNALrk08GuyWW\nbKZP13G+7duD3ZLSR0KC+mXfdZeKY38TFKVZ6+U64F/UQ/sF59iDwIPO/iPAP46wXgx0yqUcrxf2\n6o4d6vG/dasOKcXGHjs3bdM0afheQ4lOLNhyfH//rd/DiIiSaZAf6zhsHQrEsPeKFfoWlzIl8dOV\nn0qL0S38vrLjtH375NIypECKiMYzqlNHPYW8kOlySZMytuqhiMjs2Tqic911forv/OefGmfSh7A5\ngVSaxQc56/we45xf5zky7C2vczwEmMPJIedudGT2GmAV0COXNhXt/nrgdosMH663e8sWvxVb7jh0\nSKR1a5H33gt2SyxLluhk3wr/RL4tlxw9qvK8Vy//r5IYNKXZX5s3AZzmckndhQvVVOGGG04Im7Z8\nz3IJfSdUVsX6rvDs3i1yyy0ijRsH2APfD9y+caO8u3t3YAq/806RwYMDU3YAefb3Z+XyCZdLepb/\noj9cvnq1fBMf77fygo7bLXLVVSIffphnsnd275Y7yppJimjkmxEjtH/97LN+6Bv26SPy5pv5Jgu0\n0lwSN38pzS6XyNNP66IdpSTAT4lm924dF/nrr2C3pPwSHy8SHi7y44/BbknpJz1d5NZbRbp08e+y\n24WV2SXaE+i7ffs4//TTab1ypbrK//e/AEQlRnHDtzfwWc/PaFe/Xb7lpKWps1Dbturfs3kz3Hyz\n/wzMA8ETDRsyJiYGl/jqA1QAhg2Djz8udetXDr9qOCGnhHD/T/dnf7SLxNrkZLanpnJjcYVHKQ5m\nzYI9e3QJ9Tz4T/36/HzwIHvL2Jq8VaroyvH//AP79un7PmmSeq8XinffVZf3+Pj801oKTGYm3HMP\nLFmiUQXCc8YEsRSYxo3hiy/gzjv1HbAUL2433H23OrL27Bns1pR+qlRRGX7hhXDFFRAXF9z2lGil\n+cOYGB5v0ACeegqGD4dq1UhOT6bnlJ481ekpep/TO98yfv5ZQ8itWKHb0KFw6qnF0Pgi0qF6depV\nqcLP/g4/B9CwoYbtKyUh6LKpYCowqc8kNu7byFsL3ypyeWNiYhgYHl4qo0h4JTMTnnlGQ8xVrpxn\n0lqVK9MvLIxPgy2BAkS9eqo4TJ2qIcsuvxzWrStEQc2a6RfwlVf83cRyz5EjGtjlwAH44w8ICQl2\ni8oO11wDd92lj67bHezWlC9GjoTERA2NafEPFSrABx9A377QuTP8+28QG1OY4elAbOSY6luemChN\nlyyRrC+/FOnUScTtlixXllw/+Xr5z8z/5LvoxfbtItdfr0tV/lL0xQGDwldxcXJ1oDxijhzR+aPF\niwNTfgDZk7hHGr7XUKZtmlboMg5lZEjNBQskvhQu9JErH3wgcs01PtsdrU9OlvBFiySjjLvbZ2WJ\nfPyxmnk/+miupt65c+iQegxv2JBrEqx5RoHYu1cD+dx7b2CcfCx6Xy+9VCPLWIqHJUtUzuzaFeyW\nlF3Gj9dFUIqquhRWZpfYIbaxMTEMDAuj4osv6vSoMTw/53mOZBxhbI+xmFxsK1JTNbZfx44ajnjD\nhtIb5/PmsDDWHjnClpQU/xd+2mlqpvHUU0WYuw4O4dXDmdFvBg/+/CCr41YXqowv9u6le0gIdapU\n8XPrgkRCArzxBowY4bPd0Xmnn85Zp5zCj4GYzShBVKyo1iqbNulgfKtWMGFCAUbgatXSxU6eeSag\n7SwvbNmisrlHD/jsM6hU1NUCLF6pVAmmTFELow0bgt2ask9qqprEjBsHTZoEuzVll3vvVfndqxfM\nmEctmtwAACAASURBVFH89ZdIpflARgYzDhzg3u++07H4Sy7h8zWfM+PfGUy7ZRpVKnpXdH78EVq3\nVpvlNWvghRegatVibrwfqVaxIvfVr89HMTGBqeCOO9Tge+rUwJQfQC5qcBEf9/iY3t/0Jja5YCvc\nuUX4KDaWR8qSAeWwYTrXfd55Bcr2cIMGjA3U81XCqF1bTfl//ln/du5cgFUFBw6EnTvh118D2say\nzrx50KWLWoYNHVqy/UrKAo0ba1/6gQesmUageeMN9Zu68cZgt6Tsc911upDPI4+o+V2xjvsVZng6\nEBseU33v7N4td61aJRISIrJjh8zfNV/qvFtHIvd7jyOVbYrRsqXIH38UdrC+ZLIrNVVCFiyQI1lZ\ngangzz9FzjrLzytDFB9vzHtDOnzaQVIyfI9H8+vBg3LhihX5mviUGnbu1HfFIxyjr6QHcjGdEozL\npeHK69TROO0+eWVPny5y7rlq75EDrHlGvnz5pVq5lDUZXdJxuTTe7dixwW5J2WX9eg0vVwgRbCkC\nu3dr+OCHHiq4mVdhZXaJG2l2iTAuNpZHpk6FAQPYVctwy9Rb+KrPV7QMbXlC2rQ0Ha3o2BEuuwzW\nr4errgpSwwNEk2rVuKxGDSYHynv/yit1vnrcuMCUH2AGdxlM85Dm3Pfjfdkf8nz5KCaGRxo0yNXE\np9QxeDA8/nih1metUqEC99evz7jYgo3Wl3YqVID77tNZqawsfQW+/DKfEYvevaFmTfUwtPiM2w0v\nvqhmc3Pnlj0ZXdKpUAE+/VTvfzmZVCpWXC64/34dabZLZBcvjRvDwoWwa5eaeyUkFEOlhdG0A7Hh\njFr8fOCAtF+wQCQ0VJJid8l5H50nHyz94KRewuzZIs2aidx4o/Y2yjK/HjwoFyxfHriR0Y0bdQio\nwB5SJYOUjBTp8GkHeXN+/vF0d6WmSu1AjtwXN8uX6xJqycmFLmJPWpqELFggSeXYI2v5cpH27XVE\nbv36PBIuXar3O8fIPHak2SuJiSI9e4pcfrmu7G4JHq+8omHHLf5lzBiRzp3t8uXBJDNT5L//1cAP\nvi4/UFiZXeJGmj+KieHhWbNwDxrEHfMep1PDTjzW8bFj56Oi1GboscdgzBiYNk17G2WZq2vVIsXt\nZnFSUmAqaN0a+vTRYNalkFMqn8KM/jMYt3Ic0zdPzzPtJ7Gx3FmvHqdVrFhMrQsgIuqcNmQInH56\noYsJr1qV/6tZk0nlOBZxhw6wdCncdht07aqxnpOTvSS8+GI1yn3vvWJvY2lj2za45BKNvfzHHxAW\nFuwWlW9eeEEdAufMCXZLyg6HDqn4/eQTHdG3BIdKlVQkv/CCxnKeNStwdZWof/PO1FSWHTpE/2nT\nGNJmP4fTDjOm+xiMMWRkwDvvQLt2cMEFunhBaY2KUVAqGMNDDRowLpBza0OHqkvqrl2BqyOANDij\nAdP7TeeBnx9gffx6r2nS3W7Gx8UxsEGDYm5dgPj5Zzh4UFeHKCIDw8MZFxvrs4lLWaRiRfX3++cf\n/Ri2agXff+/FZGPYMBg1yi54kgc//6yOlo89ppZfJwSpEdF51O3bNXj+ggVqWxcVpSEILAGhWjV4\n6y149lnrFOgv3nwTbrpJ14KwBIisLIiN1fBHixfrFhmpK/e4XCckvecemDlToyUNHRqY59yUlI+k\nMUYGbdtGxtdfc3PV/dxe7WeW/2c5YaeFMW8ePPywhnH58ENdb6C8cSgzk7OWLmXrxRcTFqgwaUOG\n6PDQpEmBKb8YmLJhCoP/Gnzs2TnhXHw84+PimHPhhUFqnR/JytJIGSNGqDFXERERzlm+nM/POYfO\nNWr4oYGln4ULVe7Uq6ezWmef7XHyqadUwXN8AYwxiEgZMZL3DWOM5Px+uN3w2mswfjx8952ONJOe\nrsbMc+fCqlWwerXG/gsN1RVNqlaFpCQ4fFhXOmnaVEdHOnVSO/JGjYJwdWUTEe3MPPSQLn5iKTw7\nd0L79rBxo8oIi59ISoLZs3WJ0FWrdBTjjDM09GetWsc73QcPQkoKnH++yotLL9WwGiEh7N2rC6HU\nqAFffaXZclJomV0Ym45AbIDUmTNH1lx1hYQNry3r9q6TvXtF7rxTpFEjkR9+8HnNhjLL3Zs2yduB\nNOBOThapX19k1arA1VEMDJ4zWLp83kXSs05cuKTL6tUytawYVn7yiUhEhF9fiveiouT2jRv9Vl5Z\nICNDZORIkdq1RV5+WSQlO0jLwYPqLh+pEX2wNs1y4IDIddep/XJcVIYK7b59RWrUELnkEpFXXxX5\n6ae8QwxkZIisXasrGNx9t974Dh1ERowotT4XJY2FC/WbmuJ7wCGLF/r3Fxk6NNitKCNkZIhMmSLS\no4fIGWeIdO8u8v77IvPniyQl5Z4vMVFk3jyR994T6dVL80ZEiIwZIxn7EuSJJzQ4mDeVprAy2x+C\nsxsQCWwFBuWSZrRzfh3QNpc0cuWYMdJ3YKhM/ecHGTtWfdOefbZIPk5liqWJiXLWkiXiCmTv4eOP\nRa68slT3UFxul/Sa0kse+PGBY86TG5KTpUFZWf0uu3OzYoVfiz3orJK4ryytkugnoqNFbr5ZBfDs\n2c7Bt98+5lkVaKW5KHI2t7xACPAHsAX4Hajpce4FJ30kcE0u9R27P0uWiDRuLDLswV3ien6wPp+X\nXSbyv//p8n+FJSNDY9TddZdIrVoiAweK/Ptv4cuziIg+tm+9FexWlF6y/a/LWaRO/3PokMjw4SIN\nG6qy+/XXIocPF768o0dFZs4U6ddPO+t33SV/vL5EQkN1nMlTrQmK0gxUBLYBTYHKwFqgVY403YHZ\nzv7FwNJcypI37uwmD05+Tdq3F+nSJc9Va8slbrdb2q5YIb8cOBC4SjIzNeB1aV173CEpLUnajG0j\nY5drcNJH/v1XXtmxI8it8hNDh4rcemtAih6weXNgZzNKOb/8cjxqT/SWFB2yW7gwoEpzUeRsXnmB\nd4DnnP1BwHBnv7WTrrKTbxtQwUu7xO0WGTVK5P9qrZHoy2/TeOFPPum7C3tBiI0VeeklHU25917t\nyVgKxb//6iD+wYPBbknpJCJC+4OWQpKSor220FDtEK9Z4/869u3TGaozz5Sj7TrLY41nyJ23u44N\nwgZLab4E+NXj9/PA8znSfAz08/gdCdT1UpZc859rpU5dt3zxRake6Awon8bESK88Y2L5gRkzNGJ4\nKQ/Ltu3gNqnzbh35edufUmvBAolOTQ12k4rO3r2qmGzfHpDilxXHbEYpJzVVLQ1q1xaZ1e9LcXW6\nJNBKc2HlbL288nrKYidtpLP/Qo4R6V+BTl7aJYMuXyyLanSTzHrhIu+8U7RRIl9JSBB5/nl9D154\nwQ73FZL77lOTI0vBmDtXO87lOEJn4XG7Rb74QkeWb765eGaNMjNFvv1WstpeJDE1zpFn6k+S9asz\nCy2zK1E0woFoj9970FGO/NI0BE5yPW9YaRpTNhlCQorYqjLMbXXrMmjHDqLS0mhcrVpgKunVSx3M\nvvoKBgwITB3FQLOQZky+cTJ95o+n0/mP0zBQ96s4GToU7r4bzjorIMV3OOMMalaqxO+HDtGtdu2A\n1FHaqVZNfWZvvx0ef+R2XAkTAl1lYeVsONAgj7x1RSRbDscDdZ39BsBSL2WdxBnhY/m6/+W8f+5d\nuCpUwP3Pb7gBNyAY3CII4MYgiB4XUPdBPYeTJvuIEdG/iBO6xNl3fhsErm7GKf/f3nmHR1U1Dfw3\nm0KvoYUAAi/tA7tiwRYFlKIIKiICNlSqvCo2BBV7BxVBBBs2igVEQQVURKSJUl+kg5AEEhJCT92d\n74+zgYApm822JOf3PPtk995T5u5uZufOmTPT7mlaLt9C1H2D2HBZC/Y1ikJy9XO4nwscP37iuSnh\n7XDP6kCOp5JyCIj7b5gI4j6Xc8yBuP+e/Fowbn1zzIwZLmY+B+ZYWM44KGG55hFO7EfKXXRJkOOv\nc54LgkMciJi/DjHSqzhwSRia85xwXBKGUxw4CXM/HGRJGNmEUe4O4dWPHaT86cBRTshSIRPzN0Mh\nU3E/V7KPv1ayVclScGKKkTnV/bnmCK164j0EwsS8L+EC4ShhKOECEUbCXA8n4eoknCzCNBuHKxuH\nZuHQTMSVCa5McGaAK8M8d2XicqajrizUlYnLlYHTlY1LXTjViUtdxx/q/k6diiIgEaiEgUSAIwKX\nIwyRSFwSAWGRqESijghcEoFKJKvXR1DrqXCu+CUClyMcJ+G4JBynhJn3O+e9FvMtcOLAhfzroe5P\nNedvrg/9+CtHru/rie9gzvfJ/T1zv78nvo+5zguEoQhCmPuc5PX9dX9OOedO/M39/Xb/HxzvK8fl\nyhlPOPG5Sy45BaV8UgoNv5pHeHome95+mmONo5Fj25A1WxH3/znu///j2kJztII5p5pbZyiq4EJx\nqWnlcusWF3r8PXYquBqBjh9O1W0J1Jo/n2+f+e5f3wVPKa7R7GnqjVN3KObZr2HdV3nrLfM8NjaW\n2NhYrwUrrVQKC6NP3bpM3rOHZ5s08c8kIia/3803Q69eUKGCf+YJAFc1uYqquzPZvP41Drf5gCrl\nqgRbJO/ZtMmkJNi0yW9TiAiD6tfnnYQEazQXwMKFC1m4cCFntU2nX3nAfx8JeK9n82vzr/FUVUWk\noHnyPPem4yAs+R1Z8juVz2xB1bOaHzcScwzdEz+mJ4xVcZuJgjmQyyw8bhga6yHHxHQ/F8GFA5Vw\nDlcpxx8d2lI18SCNVu1A9jn55+ymZJWL5ISpLccN+FPMb1x64nnuh8stjQtwua86ZwzUfTxXm3/P\nk9swOnHOddyM53gbFWOqi7r+ZdznvOU5H2ruD+Dk6zHjiLpw4MKhTsT916FOHGTj0DTEbYiKZuFw\nZYFmopdnMmVjJpUqZOQySjPBlYG6/+I2StWZAa5sVLNAs0FdiOa8G+r+7ECOG/JhOBzhOBwRiCOC\nsLBIHI5ySFg5HI4Tz3GUQ8LKg6McGhZp/kokLkckKhE4HZVxEo5TwnGGhZEdFo77qshG3M+FbCTX\nDQrHDbzc75+6P9OcGzswRlCYQARKhNvIjBQlAiVcXETgOm7gpx92kRbhJPo0JxHiJEzdJrJm49B0\n93uejUOzQZ2IZiHqBHWCZiPqQtVp3jd1uQ1DV66byBOmfc53L+f7A4KKuA3DnP+DE+/+ie+b46Tv\notPdP9vd5vh3O+c7pCe+T64Tt6vuc7n/RwSVE8eQnLFyfXPF4f4Ucr7bgig03n2UmL1pLDi/JXHR\nVVDJhIRd7v/rXLcG7udm7OMagpPuJMh9i3miVd43xydeZ6zZQNqajQhKpGThLcU1muOB3PmAGmK8\nEgW1aeA+9i9Gjx5dTHHKBgPr16fDmjU8edppRPgro/rFF5tCDm+9BY8+6p85AsCyQ4coX64GsdWr\n029mP77u9TUOCan05J7z+OMmyaqfjdnegVjNKOHExsZy2eWXcc30nqR1GAHfLvTndN7q2ThMXHJ+\n+jdRROqp6l4RiQaSChgrT52d9Om3RbgMP9ItzdSJvuNRUzP62muDLZHHuFRxqR733Kr7WI6xAiaM\nUuSEsRAmJ7zfEQ6H20AsevasTZvg0kthyzaoWtUHFxNkXG5PeLb7PTWG4MmGVZjISY+i0Lkz/LcH\n3Hu572UvlWzaBP36mdSS770HDRoET5aOJ7/05v8Fil/cZCXQXEQai0gk0AuYfUqb2cBtACJyEXAg\n15KgxQvaVKpEiwoVmJWc7N+JXnjBhGn4ex4/8k5CAgPr12dC1/GkpKUweuHoYIvkHUuWmEIQw4b5\nfarcqxmW/Hn8p8fZWa4VPevmGbngS4qjZwvqOxu43f38dmBWruO3iEikiDQBmgMr/HNpPqJCBbM6\nNm2aqahy771w9GiwpfIIhwjhDgflHA4qhoVRKSyMKuHhVA0Pp5r7UT0igmruY1XDw6kUFkaFsDDK\nh4WZEBIvDYCWLaFjR5gwwccXFSQcIkS638fK4eFUcb9/Vd3PK4eHUyEsjEiHo8gG84oVJifz7bcX\n3rbMo2py2F96qak48v33wTWYfUixjGZVzQaGAj8CG4Dpqvq3iAwQkQHuNnOB7SKyFXgXGFxMmS1w\nfAndr7RoYcIznnvOv/P4iZSsLGYnJ3NndDSRYZF8dfNXTFkzhRn/mxFs0YqGqvEwP/tswEJlBtav\nz3t79pBlS4flyadrP+WLv7/CFX0tQ2L8+2NQHD2bX1/30C8BHUVkM3CV+zWqugGY4W7/PTBYVT0N\nEQkul18Oa9ZAejpccIGxciwFMnIkjB1bYu4xgsZzz5lF13Llgi1JiHPggKks8t578Pvvpsyqlzd1\noUhIVQQMFVlKApkuF42WLmXh2WfTqlIl/02UlAStW8Py5SWuFONru3ax9uhRPv6//zt+bPXe1XT8\npCPz+s7jnOhzgihdEfj6a7PzbNUqU+s5QMSuWsWQmBh61qkTsDlLAiviV9D18648f+OPTEhxsur8\n83E4HKitCBhafPSRudl8+WXj7SpFP9y+pkcP6NABhgwJtiShyd9/Q2ws7NxZorf4+J8//jCOtmuv\nhVdfDek7DG8rApbQ4E5LpMPBXdHRTPS3t7lOHbj/fuOOKEG4VJmYkMCg+vVPOn52vbOZ0GUC3ad3\nJ/FICYgSysqCxx4zCiiABjPAoJgY/69mlDDiD8Vzw/QbeL/b+/xwLJJBMTFeL41b/Mwdd8Cvv8Lr\nr8Pdd5uy55Y8GT7ceJudzmBLEpqMHQuDB1uDOV9UYeJE6NrVhHS+9VZIG8zFwRrNJZh7o6P5JDGR\nY/7WdA88AL/9ZoK6SggLUlOpHBbGRXnsbunZpid3nHUHPab3ICM7IwjSFYFJk6BxY7jmmoBP3aNW\nLf4+doyNdt0WgLSsNLpP787gtoM5r/E1LDxwgD7WCx/a5KySHTli4it37gy2RCHJJZeY/cXfhsi+\nzlAiKQm++MIYzZY8SEszKznjx5twjBtuCLZEfsUazSWYxhUqcHHVqkxNSiq8cXGoVAmeecYsdYby\ncmwu3klIKNAL+FTsU0RXiWbgnIGE7BLzoUMmjvmVV4IyfaTDQf969ay3GZO94O5v76ZZzWaMuHQE\nkxMS6F2nDpXDi5uAyOJ3Klc2GwT79oWLLoKffgq2RCGHiPE2v/56sCUJPSZMMCG6tWsHW5IQ5J9/\nzB1XZiYsWwbNmwdbIr9jjeYSzuCYGCbEx/vf8LvjDkhJKRGuiF3p6SwqxAvoEAdTuk9h1Z5VjFk6\nJoDSFYFXXjEe5rPPDpoI99avz6eJiRwt4+u2Ly5+kc0pm/mg2wdkqzJpzx4Gx/g9a4bFV4iYFbNp\n00xVmjfeKDEOgEBxww2we3eQFxQzMszKwFtvmexN8+ebjWVBIi3NJIF48MGgiRC6LFpkbkL79oXP\nPjPOtTKANZpLOJ1q1uRAdjbLDx3y70RhYcaIe+QRE2cbwkxKSKBv3bqFegErR1Zmdu/ZvL70deZu\nmRsg6TwkLs5o62efDaoYjcqX5/Lq1fkssQTEf/uJmX/PZMIfE/jmlm+oEFGBmcnJtKhQgTZl5Eei\nVBEbazxiH35olpTT04MtUcgQHg7//W+QvM3790P//iaf74ABJuvJ/v1G/zVoAFdcEZRMKJ9+Cm3b\nQqtWAZ86tJk40bjfP/7Y3FGUoX0d1mgu4TjcFdzGB2IJvXNno8AmT/b/XF6S4XLx3p49/9oAmB+N\nqjXiy5u/5I5Zd7Bh3wY/S1cERo0yPx6NGgVbEgbXr8/4QKxmhCBrE9dy73f3MrPXTOpXMd+p8fHx\nDLFe5pJL48Ym7/nRo9C+vQlatQDGbl2wwKy6BwRVmD4d2rQxu+zi42H1anj3XbOhbNEi42nu3dvc\n8Dz5pPFGB0i0sWOtl/kksrNh6FCzEvD77ybJdxnDGs2lgDujo/k2OZl9mZn+nUjEuCGeeQYOHvTv\nXF7y1b59nF6pUpHS8LVr2I7Xrn6N66ZeR/KxECjksmoV/PCDyZoRArSvUYN0l4sl/l7NCDESjyTS\nbWo3xnUeR9uYtgCsP3KELWlpdK9VK8jSWYpFpUrGWGvf3lQ+Xbcu2BKFBFWrmgJuEycGYDKXCwYO\nNL8nX38Nb78N1av/u114uGm3erX5nNq1g9RUv4v3889mgfXKK/0+VcngwAGTHWPrVli6FJo1C7ZE\nQcEazaWAqIgIetSuzfuBqOB21lnQpQu8+KL/5/KCCV56AW876zZ6tu7JDdNvINPp55uPglCFhx4y\nJYFDpK6tQ4TBMTGMj8+zknKpJD07ne7Tu3P7Wbdzy+m3HD8+ISGBe6Oj/Ve+3hI4HA5jsD3/PFx1\nFcwNsRCtIDF4MLz/vp8jV1RN5cb1600M88UXF94nJsYY17GxZq+Hnx0348ebvNVlKPIgf7ZvN59R\nq1bw3XdQrVqwJQoaVvOXEobUr8/EhAScgVhCf/ZZE6IRYumb1hw5wj8ZGVwXFeVV/xfav0BUxSgG\nfhfEjBpz5kBCAtxzT3Dmz4fb69bl+/37SfT3akYIoKr0n92fRtUa8VTsU8ePH8rOZmpSEvd6GPpj\nKSHceivMnm1iEwLiYg1tWrSAc86BGf4qnKpqYh7++MPcqFSu7HlfERO2ceGFxnlz+LBfRNy1CxYu\nNHvcyjzLl5t0jUOHwptvGs9/GcYazaWE86tWpW5kJHNSUvw/WUyM+Qd6/HH/z1UExsfHMyA6mnAv\nvYAOcfBJj09YtXcVry551cfSeUBWlknr9+qrIaeYqkdE0LN2bSaXgfRzz//2PFtStvDR9R/hkBPf\npY/37qVDjRrUL6VJ+8s0F18MixebINZHHjGhA2WYIUOMp9UvvPCCKTrz44/eeSxFjPF2+ummlKEf\nMvu8+64JUymKPV8q+fprU91v0iRbLtKNNZpLEffFxDAuUEvoDz9sNmksXRqY+Qphf1YWX+zbxz3F\n9AJWjqzMt72/ZdyKcXz999c+ks5D3n3X3JB07RrYeT3kvpgYJiYkkFWKDYpp66cx+a/JxzNl5OBS\n5e34eIbaDYCll//8x2wQXLrUpKUL0IazUKRrV0hMNM5gn7JoEYwbZ5b4a9TwfhyHwyRQdjqNEe5D\nMjLgvfdsMRPGjTMhND/+aAxnC2CN5lJFzzp1WH/0KBsCUcGtcmWjrB54ICS8Mu/t2UO3qCjqRkYW\ne6wGVRswq9csBnw3gJUJK30gnQfs32/iK8eODdkgujMqV6Z5xYp8tW9fsEXxC0t3L2XY98P4tve3\nRFeJPunc/NRUyjscXF6GY/nKBFFRJjdwZqZZ/i9jm19zCAuDQYN87G1OSTHxDu+/D74IcQoLMznh\nxo83qwQ+4osvzNadli19NmTJQtWsIudU+Dv33GBLFFJ4bTSLSE0RmS8im0Vknojkse0VRGSniKwV\nkVUiUnLqMJdAyjkcDIiODpy3uW9fc6c/dWpg5suHbJeL8fHxDGvQwGdjnlf/PN677j2un3Y9uw7u\n8tm4+fL003DjjXDGGf6fqxgMi4nhrVK4IXB76nZumHEDH3X/iDPrnvmv82/FxTGsQYN8K0z6iyLo\n2U4islFEtojIo570F5ER7vYbReTqXMcXuo+tcj/KVqqQ8uVNQG+rViY/8N69wZYoKPTvD7NmQbIv\nEgqpmgFvusm3K2kxMcYt3KePcTz4gJwNgGWSrCyTv/znn82NSOPGwZYo5CiOp/kxYL6qtgB+cr/O\nCwViVfUcVb2gGPNZPGBg/fpMS0oiNRAFSBwOU1nrscdMztMgMTslhYblynFelSo+Hff6Vtcz/OLh\ndP28KwfT/bhT+++/TUWlZ57x3xw+4rqoKBIyMvijFHngUtNS6fp5V0ZeNpIuzbv86/zmY8f44/Bh\nbi2gwqQfKVTPikgY8DbQCWgN9BaR/yuov4i0Bnq523cCJsiJOwIFbnXr7HNUNQTyMAaYsDCTAu2G\nG+Cyy0Ju03MgqFULunWDKVN8MNjEiabcoD+yLl17rYltvvfeYg+1Zo1JFV0moxHS081NTVKSKTVv\n02rmSXGM5m5Azr/TFKB7AW1Dc725FFKvXDmujYoKTPo5MHXnL7nEbF4LEm+6vYD+4IGLHuDyRpfT\n84ueZDn9dCPy0EMwYgTUru2f8X1IuMPBkFLkbc7IzqDH9B50btaZoRcMzbPN2/Hx3BMdTfmwsABL\nB3imZy8AtqrqTlXNAqYB1xfS/3pgqqpmqepOYCtwYa4xrc4WgSeeMGXyLr/c3NyWMQYMMFstipVM\naM8eU5Tkk0/AX5toX34Z1q6Fb78t1jDvvgt3323umcoUhw+bFYAKFczygq12mi/FMZrrqmpObd1E\noG4+7RSYJyIrRSS08miVUobFxPB2fHxg0s+BUVhvvx3AMlInWH34MNvS0ujhp7tiEeHNzm9SLrwc\nA74b4PtUdHPmwJYtZsNFCaF/dDTfpaSwt4RvlMpJLRdVMYrXrn4tzzaHsrP5NDHR4wqTfsATPRsD\n7M71Os59rKD+9d3tcvfJfZEfukMzRhVH+FLB0KEncjn/9VewpQko7dpBZCT88ksxBnnoIWOJtm7t\nM7n+RblyJq5i2DA4dsyrIY4cgWnTTBRJmSI11VT2a9bMrHj6YF9QaabAvFYiMh+ol8epkblfqKqK\nSH7WxCWqukdEagPzRWSjqv6WV8PRo0cffx4bG0tsbGxB4lnyoW3VqkSXK8fs5GR6BMJ7edppRlkN\nHw5ffun/+XLxVnw8g2Ni/FpsItwRztQbp3LFR1fw7KJnefKKJ30zcEYG3H+/KUlaghRVzYgIetWu\nzcSEBEY3aRJscbzmyV+eZOv+rfx8+88npZbLzYd799KxRg0alC9/0vGFCxeycOFCn8jhAz176jHJ\n41hhejo3fVQ1QUQqA1+JSD9V/SSvhmVGZ+fkH+vc2dzonn9+sCUKCCInvM1XXeXFAD//bDaTTZrk\nc9n+RceOJn/zCy/Ac88VufvUqWZBoUwlyNm/37xvV1xhqv2G6CZ0X+Azna2qXj2AjUA99/No2ZoL\nJgAAIABJREFUYKMHfZ4ChudzTi2+Y3piol7611+Bm/DYMdUmTVQXLAjYlHvS07X6b7/pvoyMwMx3\neI82eaOJfvDXB74Z8IUXVK+7zjdjBZj/HTmidRYv1mPZ2cEWxSsm/jFRm73VTBOPJObbJsvp1CZL\nl+qSAwcKHc+tv7zWp/k9PNGzwEXAD7lejwAeLag/Jrb5sVx9fgAuzGPs24Fx+chW6PtS6vjmG9Xa\ntVWXLQu2JAEjNVW1WjXVxPz/VfImI0O1ZUvVWbP8IleexMerRkWpbtxY5K7nnac6d64fZApVkpNV\nzz5b9eGHVV2uYEsTcLzV2cVxz812K9QcxTrr1AYiUlFEqrifVwKuBtYVY06Lh9xQqxZxGRksD9SG\nrQoVYMwYE2YQiE2ImFjTW+vUoVaAvLT1Ktfj+z7fM+KnEXy/5fviDRYXZ+7sx471jXABpnWlSrSt\nUoVPEhMLbxxifLPxG57+9Wl+6PMDdSrlv7lvZnIy9SMjuTi4aeYK1bPASqC5iDQWkUjMBr/ZhfSf\nDdwiIpEi0gRoDqwQkbCcbBkiEgFch9XZJ+jWDT74AK67DpYtC7Y0AaF6dbPP7sMPi9hxzBho3ty8\nZ4Gifn0YNcqE1BQhlO7PP02WkKuvLrxtqSA52SwddOpkwitLsYfZ53hjaRsjnZrAAmAzMA+o7j5e\nH5jjft4UWO1+rAdGFDCef28ryiBv7t6tPdevD9yELpfqNdeovv6636c6kp2ttRYv1i1Hj/p9rlNZ\nsmuJ1nqllq6IW+H9IL16qY4a5TuhgsAv+/dri2XL1FmCvBQ5n90f8X8U2M7lcmnblSt1ZlKSR+Pi\nP09zoXrW/bozsAmzoW9EYf3d5x53t98IXOM+VgljhK9x6+yxgOQjm0fvTanku+9U69RRXbUq2JIE\nhKVLVZs2VXU6PeyQlGQ8vlu2+FWuPMnKUm3Vqkhu43vuUX3uOT/KFEocOKB67rmqjz5aJj3MOXir\ns32u5L19lGkF7CcOZ2Vp1G+/6bZjxwI36caNRlnGxfl1mnG7d2uPdev8OkdBfLPxG633Wj3duK/o\ny4C6YIFqo0aqQTD4fYnL5dLz/vhDv9m3L9iieMT6xPVa59U6Ondz4T+mi1JTtdmyZZrt4Y+Kv4zm\nUH6UeZ39xReq9eqpbtgQbEn8jsuletZZqvPmedhh2DDVoUP9KlOBzJypesYZqh6Ejx06pFq9umpC\nQgDkCjaHD6u2a6d6331l2mBW9V5n24qApZjK4eHcU78+b8TFFd7YV7RsCQMHmkqBfsKpyti4OB5q\n2NBvcxRGt5bdeOGqF7jm02uIO1SE9zc93ZTaGjcOKlb0n4ABQER4qGFDXtu9u/DGQWbngZ10+qwT\nY64eQ+fmnQtt//ru3TzYoAFhdtnSkh833QSvvGI2Um3bFmxp/IqISYM8ebIHjbdvN1kYnnjC73Ll\ny/XXQ9WqpmJgIUydCldeCdHRhTYt2aSnm/elVStTX8HqNq+wRnMp576YGD5NTGR/gOKMARg50gSJ\nfV/MuN98mJWcTN3ISNoFuaTxnefcyZC2Q7j6k6tJOZbiWadXXoE2bQIb5+dHbqpdm13p6awI4WIn\nSUeTuPqTq3m43cP0ObNPoe03HTvGkkOHuL1eXgktLJZc9OtnYmivuQZKYHx/UejTB+bNM7UvCmTk\nSJMVKDjFgAwiRtc+8QSkpRXYdPJkuKe0J8N1Os0HGBVlMpn4MdtUace+c6Wc+uXKcX2tWkwIZDGK\nChVO1CL1Mmdmfqgqr+zaxfAgeplz8/AlD3Ndi+vo8nkXDmccLrjxli0mvdybbwZGuAAQ7nBwf4MG\nvBqi3uaD6Qfp/FlnerXpxbALh3nU5/XduxlYvz4Vy1yFA4tXDBwIffua4hCHC9EBJZhq1cyGwAIr\nBP7xByxa5NeVRo9p186kBhw3Lt8mq1aZm4BSvQFQ1aSETU01BWasXisW1mguAzzasCFvxcdzJDs7\ncJN26gRt25qiAD7kp9RUDjmddA+hEp8vdXiJs+ueTbdp3UjLyseroWpuIkaMgEaNAiugn7k7OppF\nBw7wdxBLqefF0cyjdP28K5c0vIRnrvSsRHlcejpf7tvHsDKVrNVSbJ56Cs4914RsZGYGWxq/cc89\nxjOr+SWmePxx816ESkW5F1801WoPHMjz9OTJZaAC4AsvmFzZM2f6ryJjGcIazWWAVpUqcWX16rwb\nqNLaObzxhlkKWrvWZ0M+988/jGzUKKRiTUWECV0nUL9KfW6ccSOZzjx+ND/91CzfDvPM21mSqBwe\nzn8bNODFXbuCLcpx0rPT6T69O82jmvNGpzcQD78vr+7eTf/o6IClMbSUEkRgwgQoX95YlvlalSWb\niy82dZh+/TWPk7/+Cjt2wJ13BlyufGnZ0qQHHDPmX6eOHjUVAENJXJ/z8cfw3nsmVDLI4YylBWs0\nlxFGnnYar+/eTZrTGbhJo6PhpZfgrrvAB17u3w4cYHdGBrcEM1YuH8IcYXx0/UeUCy/HrV/dSpYz\nVwz53r2mlOwHH0BERPCE9CNDYmKYm5LCtkLiBwNBpjOTXl/2omaFmrx33Xv5Vvs7lb0ZGXySmMjw\nBg38LKGlVBIeDp9/DuvXGw9nKUTkhLf5JFRN/PCTT4aejnviCRMumJx80uEZM+DSS6HU/rv/9pv5\n3Zkzpwzscgwc1mguI5xZuTJtq1Thg717AzvxXXdBjRp53ukXlef/+YcRjRoRHqKbGCLCIph24zTS\nstPoO7Mv2S73jcLQoeZ9OO+84AroR6qFhzM4JoaXg+xtznJmccuXtwDwSY9PCHN4vu46Ji6OPnXr\nUs8uYVq8pVIlmD0b3nkHvvgi2NL4hX79jB2Wknvv8/z5sG+f2WwWajRpAjffbMI0cjFpkskIUirZ\nuhV69jQrnK1bB1uaUoVoiCwjiYiGiiyllRWHDnHT//7H1gsvJDKQhueOHSa++fffzXKZFwRNdi9I\nz06nx/QeVC9fnU+d3Ql78ilYvdos3ZZiUrKyaLF8OavOP59GQbjWLGcWvb/qTaYzky9v/pLIMM9D\nLFKysmi+fDmrvZRdRFDV0IkZCgBWZxfAqlVmd9mcOXDBBcGWxuf06wfnnAMPPojxMl90kXnRq1ew\nRcubuDg46yz43/+gXj3WroUuXWDnTrNAUKo4cMB8HvffbzapWvLEW50d2taHxadcULUq/1exIlMC\n7W1u0sRsDunf36S+8YLn/vmHRxo2DHmDGaB8eHlm9ppJVtJeDg28k+z3JpV6gxkgKiKC/tHRvBIE\nb3OWM4tbv76V9Ox0vuj5RZEMZoA34uK4sXbtoBj7llLIOeeYGIYbbzThWaWMe+81nlpV4LvvTFq3\nnj2DLVb+NGhgMpy89BIA775rNgCWOoPZ6YRbbzU3bNZg9guhb4FYfMroxo159p9/SA9kbDOYzBHh\n4fD660XuuuzgQVYdOUL/EhSXVT6sHNPnVePni+rSe8+4vDcHlkIeatiQqUlJ7AxgbHNGdgY9v+hJ\nWlYaX978JeXCixZekZSZyYT4eB4vZVlNLEGme3djmd14Y6nLqHHppSbjxK8L1ThEnn469HP/jhgB\nH3/MsS3xTJ1qPppSxxNPmCImXvzOWjwjxL/lFl9zcbVqnFu5MhMSEgI7scNhEny+9ppZuvQQVeWx\n7dsZ3bgxFUpSXqApUwjbtp1rZ6wm05nJTTNuIj07PdhS+Z06kZEMjYnhqZ07AzJfWlYa3ad3J9wR\nzte9vqZ8eNE9xc//8w996talSYUKfpDQUqZ54gmoXbvUZc0RgQED4M8nZ5kD3bsHVyBPqFcP+vdn\n54AXS+cGwBkzzEbU6dNDbzNmKcIazWWQ55s25aVduzgYyLzNAKedZjYE9u1baJWmHOalprI3M5Pb\n69b1s3A+ZMcOePhh+PRTylWuxpc9v6R8eHm6Te3G0czQymXsD4Y3bMgP+/ez/sgRv85zOOMwXT/v\nSlSFKKbdNK3IIRkAO9PS+DQxkZGnneYHCS1lHofDpP1atMjEM5Qi+vVx0WnJUxwc/kzJKcn8yCPE\nLJrKf3uETnpMn7BunVnNnTnT3KRZ/IbXRrOI9BSR/4mIU0TOLaBdJxHZKCJbRORRb+ez+I42lSrR\npWZNXg9GFbc+feD00+Gxxwpt6lJlxPbtPN+0achmzPgXTifcdhs8+iiceSZgsmp8fuPnNKjagPYf\nt/e85HYJpWp4OI81asTIHTv8NkfS0SSunHIlLaJaMKX7FMId3gUnPrVzJ0NiYqgbonmZRaSmiMwX\nkc0iMk9EqufTLk89m19/9/FfROSwiIw7ZazzRGSde6zSU74yWFStCrNmmfLSK1cGWxqfUePnr6hQ\nozyT4rsGWxSPWRVXm08r3cuVy14Itii+4+BBuOEGGDvWxNJb/EpxLJF1QA9gUX4NRCQMeBvoBLQG\neovI/xVjTouPeLpJE8bHx5MY6Fg7EZOOaeZMk5qpAGYkJREuwg0hVP2vUJ55xiyNPfjgSYfDHeG8\n3+19rjjtCi778DJ2HwzNstO+YlD9+qw+coTfDx70+dg7D+zksg8vo0vzLrzT9Z0ipZXLzfojR/hh\n/34eCpGS7PnwGDBfVVsAP7lfn0Qheja//unAKOChPOZ8B+ivqs2B5iLSyYfXUzZp0QImTjQVA1NK\nwU2z0wmjR5Px+NNMmiy4XMEWyDPefRfSBj+E48svTOqMko4q3HGH2fjXt2+wpSkTeG00q+pGVd1c\nSLMLgK2qulNVs4BpwPXezmnxHaeVL89t9erxdDAUR82aJu7qnnvyVVwZLhejduzgxaZNPa7mFnTm\nzTPVlz7/PM9NMSLCyx1fpv85/bn0w0tZl7guCEIGhvJhYTzduDGPbNuGL9OSrd67mss+vIwhbYfw\nzJXPFOu78dj27TzaqBFVQ3sLfTdgivv5FCCv4NGC9Gye/VX1mKr+DmTkHkhEooEqqrrCfejjfOa0\nFJUbbzRGc9++lBgrMz9mzICqVWl1fycqVoSffw62QIVz+LD52bn1vigYPBieey7YIhWf11+HhASf\n1EGweIa/17xjgNwutTj3MUsI8MRpp/H1vn2sOnw48JNffLEJ0ejZEzIy/nX6td27Ob1SJdrXqBF4\n2bwhLs6EZXz2mdlwUgDD2w3npfYv0f7j9szbNi9AAgaefvXqkeFy8Wliok/Gm7tlLld/cjVjrxnL\nsAuLt7Hqu+RktqSlMTQm5NVRXVXNeQMTgbyC+wvSs4X1P/WOJsbdP4d4rM72HS++aOo3l2SDLTsb\nRo+GZ55BHMLgwaaCeKjz8cfQoQPUr49ZCZw1C7ZsCbZY3rNokdlY/8UXYAsyBYwCXSwiMh/IywJ4\nXFW/9WD8IrmYRo8effx5bGwssbGxReluKSI1IyJ4rkkThmzZwuJzzsERaI/u/febUp/Dh8Pbbx8/\n/E96OmN372ZlSamgl5Vlkvr/97/g4Xe29xm9aVC1ATd9cRPPX/U8d59b+vIfhYkwvkULeqxfT7da\ntahWDI/uhD8m8OyiZ/nmlm+4uOHFxZIrzenkv1u38k6LFl7n/V64cCELFy4slhw5FKBnR+Z+oaoq\nInnp1FOPSR7HCurvNVZnF5GICOPuPO88uOQSaN8+2BIVnY8+MpZnhw6A2aby+OOwaxeEatZGVfMT\nM3Gi+0CNGvDAA6bs99SpQZXNKxITTT7mjz4K3Tc9xPCZzlbVYj2AX4Bz8zl3EfBDrtcjgEfzaauW\nwON0ufSClSv1w4SE4AiQmqrarJnq++8fP9Rj3Tp9eseO4MhTVFwu1QEDVK+9VtXpLHL3TcmbtNlb\nzXTY3GGamZ3pBwGDT/+//9b/bt7sVd+M7Awd9N0gbfV2K92astUn8ozesUNvWLfOJ2Pl4NZfxdan\npz6AjUA99/NoYGMebfLVs4X1B24HxuV6HQ38net1b2BiPrL59D0sUyxYoBodrRosvestaWmqDRuq\nLlly0uH771d97LEgyeQB8+ernn66UdfHOXxYtW5d1dWrgyaXV2Rnq3booDpqVLAlKdF4q7N9FZ6R\nn4tyJWYjSWMRiQR6AQXv/rIEFIcI45s3Z8SOHRzIygq8ANWrw7ffmsTzCxfy4/79rD1yhEdCe3PW\nCd56y5QH/+wzr5L7t4hqwYq7V7B5/2au+fQako8l+0HI4PJi06Z8npTEuiKmoEs8kkiHjzsQdyiO\nZf2X8Z+a/ym2LDvS0hgXF8fYZs2KPVaAmI0xbHH/nZVHm4L0bGH9T9LdqroHOCQiF4oJGO+Xz5yW\n4tC+vUl03Lu3CXcoKUycCGefbcLrcjF4MLz/vqmrEYq8/TYMHXpKZrzKlY2LfNSooMnlFc89dyJE\nxhJ4vLG0jZFOD0wcXRqwF/jefbw+MCdXu87AJmArMKKA8fx5U2EphHs3btSBmzYFT4D58/Vogwba\n/Lff9Lvk5ODJURTmzDHeIh94xbOd2fro/Ee18RuNdUXciuLLFmJMiIvTS/78U7NPcvXkz5JdS7TR\n2EY66qdR6nQV3YOfFy6XS7uuWaPP79zpk/Fyg/88zTWBBcBmYB5QXYugZ/Pr7z63E0gBDrt1eSv3\n8fMw2ZG2Am8VIJvP38cyRY7HcOTIYEviGYcOGc/smjV5nu7USfWjjwIskwfs2KFas6bqkSN5nExP\nV23USHXx4kCL5R0ldYUiBPFWZ4vpG3xERENFlrLIgawszly5kkktWtApKiooMtw3cyb7t2zhs7vu\nglBPM7d2rYnpmzUL2rXz2bBfbfiKQXMGMfKykQy7cFjJyRxSCE5Vrlq9mmujoni4gBg8l7oYs3QM\nry55lUnXTuL6Vr5LtjM5IYF3EhJYdu65Xscy54eIoKql48PyEKuzfUBSEpx7rsm60ynEM/s9+yxs\n3GhW1fJgzhxTUfuPP0Kr1sljj5kq5vkmmPjgA1OtduHC0BL8VBIS4Pzz4dNP4aqrgi1NicdbnW2N\nZstxfk5N5ba//2ZN27ZEBbgM54/793Pvpk2smTuX6t9/Dz/9BNWqBVQGj9m0Ca68Et5802T/8DHb\nU7fT68teNKjagPeue4+oisG5ifE1/6Snc/6ff7LgrLM4q3Llf53fd3Qfd82+i31H9zH9pumcVt13\nVfq2HjvGRX/9xaJzzqF1pUo+GzcHazRbvGbRIrj5ZmNthmpY2t69pijV8uXwn7zDpFwuaN7c2NQX\nXRRg+fIhLc0Uol2yBPKNyMrONoWoXn4ZrrsuoPJ5THa2Cenp0MGUZrcUG291dgkps2YJBFfVqEGv\nOnUYsGkTgfwxTMnKov/GjXzYqhXVn3nGaNxrrzWpmUKNHTugY0d44QW/GMwATWs0ZfGdi2lavSln\nTjyTOZvn+GWeQHNa+fK8/p//0Pfvv0l3Ok86983Gbzhz4pm0qd2GRXcu8qnBnO1y0W/jRp5o3Ngv\nBrPFUiwuv9xkErrlFpOJJxQZNcoU0cjHYAazpWPYMJM6OFSYMgUuvLAAgxkgPNy4oYcPNy7pUOTJ\nJ01auccfD7YkZR7rabacRLrTyfl//snDjRpxeyH5hn2BqnLzhg00KleO13M0m8sF/ftDfDx88w1U\nqOB3OTwiLg6uuMLk+BwyJCBT/rrzV+745g46NOnAa1e/RrXyIep99xBVpef//kej8uUZ06wZqWmp\nPDjvQRb9s4gp3adwaaNLfT7nszt3sujgQX4880y/pVW0nmZLsXC5jJfz//7P5N4NJVavNqEjGzea\njdsFcOQINGkCS5cWYqgGAKcTWrY0Wdku9UStdO4M11xjbmBCiblz4d574a+/oE6dYEtTarCeZotP\nKB8WxtTWrXl42zZWHDrk9/le3LWL7WlpPN+kyYmDDgdMnmwUxDXXwIEDfpejUDZuNJp38OCAGcwA\nVzS+gjUD1+AQB60ntGb6+ukBXQXwNSLCuy1bMjM5mUErZ9J6QmsqhldkzcA1fjGYZycnMzEhgY9a\ntQp8HnKLxVMcDlN946uvzCNUUDX5jEePLtRgBpOQYuDA0ChQN2uW+Qm55BIPO7z+Ojz/fGiVOd+x\nA+68E6ZNswZziGA9zZY8+TY5mYGbN7Pk3HM5rXx5v8zxRVISw7dtY9m551I/r4pGLpdZMluwAH74\nAYJVvW35crj+enjpJbNEGSR+3/U7A+cMpEHVBrxxzRu0rNUyaLIUhw37NnD3guf4o3Zv3j4tigHN\nfbeRMjerDh/m6rVrmXvGGbStWtUvc+RgPc0Wn7ByJXTpYuKcW7UKtjTG8hw1ynibPSxOlJRkRN+4\nMXh2nqoJyxgxAnr0KELH++4zf8eN84tcRSItzVj8t90Wet7vUoD1NFt8ynW1avFQw4Zcu24dh/yQ\nR3T5oUMM3rKF2aefnrfBDMb7MmYM9O1rlMfq1T6Xo1BmzTLLpu+/H1SDGeCSRpfw171/0b5Jey79\n8FKGzh3KvqP7gipTUUg8ksjA7wYS+1EsNze5gNlnX8RTScrWY8d8Pld8Rgbd1q/nnebN/W4wWyw+\n4/zzzX6JG24wsQ7B5OhRE4o2dqzHBjMYQ7lXr5OKvAacRYvg4EHo1q2IHUePNhUbg/FbkxtVs6LZ\nooWpNGsJHbzJU+ePBzbnZ8jhcrl04KZNeuWqVXooK8tn464+fFijf/9dZ+/b53mnqVNVa9VSnTjx\nlLJOfiIjQ/WBB1RPO011+XL/z1dE9h3dp/fNvU+jXo7SUT+N0uSjoZvbOulIko5YMEJrvlxTH/jh\nAU05lnL83IS4OG22bJnuTEvz2Xx70tP1zBUr9AU/5GPOD/yUpzmUH1Zn+5G77lLt2TMwui4/HnxQ\ntU8fr7pu2WLU9eHDPpbJQ7p0UZ00ycvOH36oeu65qj78zSsyEyeqtm4dvDewDOCtzraeZku+iAjj\nmjWjeYUKxK5eTaIPdhb/kppKxzVreKNZM64rSi7mW26BxYth/Hjo08e/cc7btsFll8HWrWbzxQUX\n+G8uL6lVsRZvdX6L5XcvZ++RvTQf15xH5j9C3KG4YIt2nF0HdzH8x+G0fLslqWmprLxnJWOuGUPN\nCjWPtxkUE8PQmBgu+esv1vrAs7bl2DEuWbWKm2rX5rEC8kFbLCHN+PGwa5fJjRwMli83uePeeMOr\n7s2aQWwsTJrkW7E8YfVqo7b79fNygNtvh6io4AVm//KLyZYxc6YJEreEFt5Y2v54YL0WIYvL5dKn\nd+zQpkuX6uajR70eZ0ZiotZevFh/3r/fe2GOHlUdPNhUpnr/fVWnb6rFHR/7iSdM+agxY4Lr5Ski\n/xz4R4fNHaY1XqqhPab10Pnb5vuskl5RcLqc+uPWH7Xb1G5a8+Waev/39+vug7sL7Td1716tvXix\n/lKM78aKgwe13u+/6+T4eK/H8Basp9nia/bsMdXqZswI7LwZGapt2qh+/nmxhlm3TrVOHVNIMJB0\n6aL65pvFHGTHDtWoKNVAV8ndssX8ti1YENh5yyDe6my7EdDiMZMTEhixfTtPNG7MkPr1Cfewqtq+\nzEwe3b6dBampfHvGGXkWtigyf/4JQ4eazYKjR5ssG95WeUtPN7uTR482OaJffTV0iwwUwpHMI3y2\n9jPeWfkOKWkp9Gzdk15tenFBzAV+qy7oUhfL45Yz/X/T+WLDF9SrXI9B5w+i9+m9qRTpeV7kn1NT\nuWXDBm6rW5enGjemiodxlGlOJy/u2sWE+Hg+aNWKbkGoJmk3Alr8wurVJi/8Dz/AeecFZs6nnzYb\nEmfPLnaFvH79jNf5qad8JFshLF5stsBs2mTSGheLN9+Er782nl8fVxDNk4MHze/PsGEwaJD/5yvj\n2IqAloCw8ehRhmzZQkpWFq83a8aV1avnm8rrqNPJJ3v38uTOnfStW5fRjRtTtQgbSgrF5YKpU02q\noCNHzMaJm2+G6OjC+6rCli0m+/1778E555h6q7GxvpMvyGzYt4Hp66czY8MM9qftp0PTDnRo0oF2\nDdvRPKo5DvHuh8DpcrJl/xaW7F7C/O3z+Wn7T9SqWItebXrR6/RetKrl/a7/xMxMHtm2jZ8PHOCV\npk25oXZtyuXzg5XlcjEnJYXh27ZxXpUqjPnPf2jgp0wvhWGNZovfmDnTOAh+/dX/yY9//dXs4lu5\nEho0KPZw27eb6LaNG8Hf97Kqpk5M//4+2rPtdJrKrx07+r8KX3q6yZrSpk1oZO4oA1ij2RIwVJXp\nSUk8v2sX+zIzuTYqitjq1Y8bN8lZWcxJSWHRwYNcWq0aLzVtypn+jM1Shd9/hwkTjEemXj2j7Fq3\nNlu569QxZUiTkiAx0QS8/fyz8aL06GGM7ZYlM32bp+w8sJMF2xewYPsCVsSvYN+xfZxZ90xaRLWg\nUdVGNKzWkFoVa1ExoiIVIyqiqhzLOsaxrGMkH0tm96Hd7D60m03Jm1iXtI46lepwYcyFdGjagfZN\n2vu0gh/AbwcO8PiOHaw7coSONWvSqWZNqoaFAeZmbH5qKt/v30/zChV4pkkTrqlZs5AR/Ys1mi1+\n5d13TZnn337zX+rNuDhj4U6ZYgxFHzFkiKlP5e+aLd9/bzKUrlsHblVRfPbsgbZtTXB2ly4+GvQU\nsrLgppvMm/TZZz4U3lIQATeaRaQnMBpoBbRV1b/yabcTOAQ4gSxVzXNXlVXAJZPtaWl8m5LCkoMH\ncbmPVQ4L45oaNehUsybVIyICK5DTaZY0f/nFbOhLSjKPiAhjPNeube7m27c3XpsyWvDiQPoBVu9d\nzdb9W9l9cDe7Du0iNS2VY1nHOJp1FEGOG9A1K9SkYdWGNKrWiGY1m3F2vbMDVpkwKTOTOSkp/JSa\nSoZbP0SIEFu9OtdGReWfrjDA+MtoFpGawHTgNGAncLOq/msXrIh0At4AwoD3VPXlgvq7j38FnA98\npKr35RprIVAPSHMf6qiqyXnMaXV2IHnlFVPebtEi37ttMzJMtdPu3c2Kmw/ZswdOPx3WrPGJ8zpP\nXC4TvfLEEyZbn09ZssQ4V5YsKbCMuFe4XCYPc2qqWVGIjPTt+JZ8CYbR3ApwAe8CwwvGIzWeAAAM\ndklEQVQwmncA56nq/kLGswrYYrGUSPxoNL8CJKvqKyLyKFBDVR87pU0YsAnoAMQDfwC9VfXv/PqL\nSEXgHOB04PRTjOZfKECn52pndXagGTECfvzRlFauV883Y6rCgAGQnGyqEfrBkTBypAnVmDrV50MD\nJo3+pEmwbJmf/CDjx5sJFi+GKlV8M2ZWlnnft241K6QVK/pmXItHBLy4iapuVNXNHjYvm+48i8Vi\nKR7dgCnu51OA7nm0uQDYqqo7VTULmAZcX1B/VT2mqr8DGfnMa3V2KPLCC6Y6abt2ZrdbcXG5TBW8\nv/4yXmw/rbyNHAkrVsCcOb4fe88ecy8xaZIfFw4HD4ZLLzVhK6mpxR/vyBHzOe7ZY26ArMFcYghE\nnmYF5onIShG5JwDzWSwWS2mhrqomup8nAnXzaBMD7M71Os59zJP++bmKPxSRVSIyyguZLf5CxKSi\nGDXKhFMsWeL9WE6n2TW3ejX89BP4sXJmxYrGqB00CA4f9u3Yw4bBPffAWWf5dtyTEDElDtu1M/tl\nkpK8Hysx0Ww4j442GUpsLuYSRYGpDERkPia27VQeV9VvPZzjElXdIyK1gfkislFVf8ur4ejRo48/\nj42NJbYUZTKwWCylh4ULF7Jw4UKfjFWAnh2Z+4WqqojkZeSeekzyOFZQ/1Ppo6oJIlIZ+EpE+qnq\nJ3k1tDo7SNx1lzG6rr/eeEFHjixaPOyBA8bSPHDAhHtU8jw1pLe0b28ctSNG+K7E9qxZsHYtfJLn\nt9PHiJhMTU89ZW5YZs6EVkXMFPTll8azP2iQCcAuo3tqgoGvdHaxs2d4Gv/mbvsUcERVX8/jnI2P\ns1gsJRI/xjRvBGJVda+IRAO/qGqrU9pcBIxW1U7u1yMAl6q+XFh/EbkdOD93TPMpY+d73ursECA+\n3hhg27eb1JkXXVRwe1X4/HN46CGzuW3MGAhgmsb9+82mwBkzTLRDcUhNhTPOMJdz+eW+kc9j3n3X\nePsHDDA3LBUqFNw+IcEYyxs2mADsdu0CI6clXwIe03zq/HkeFKkoIlXczysBVwPrfDSnxWKxlHZm\nA7e7n98OzMqjzUqguYg0FpFIoJe7nyf9T9LdIhImIrXczyOA67A6O3SJiYFvvjGG2003mZRxEyYY\n6zQ327cbQy821uR+mzXLtAtwXvOaNU2YRq9e8M8/3o+TmQk9e5q0/AE3mMEYy2vWmFz/bdqYWPOV\nK02MeG4hZ80yGUlatzZe6VWrrMFcwilO9owewFtALeAgsEpVO4tIfWCyqnYVkabA1+4u4cBnqvpi\nPuNZr4XFYimR+Dnl3AygESenjDuuZ93tOnMi5dz7OXo2v/7uczuBKkAkcADoCOwCFgER7rHmAw/m\npZytzg4xsrNhwQKzoW+2+56penVTzS4728RGdO5sLE1fFpnygrFjjcP199+hWhGzV6qaUOzkZBMh\nEfS0xosWmcqB8+aZeOXKlY0b/NgxYyDfeae5ofFV1g2LT7DFTSwWiyVI2OImlpBC1RhtBw6YHMxN\nmoRU/KyqiVbYtMkkjyhKOv/nnjMO3F9/DUgodtFISDAe5ho1jJEciPLbFq+wRrPFYrEECWs0WyxF\nIzvbhFWLwIcfQlRU4e2ff960XbrU7IO0WLwl2DHNFovFYrFYLB4RHm6SSbRoYdLFLViQf9sdO0zC\nit9+MyEd1mC2BAvrabZYLJZiYj3NFov3zJ9vQn8vuQQ6dIDLLjOxzosXm5DhadNMde8HHrARDxbf\nYMMzLBaLJUhYo9liKR4pKWZj36JF5nHwoDGeL7sMrruu6CmRLZaCsEazxWKxBAlrNFssFkvJwcY0\nWywWi8VisVgsfsIazRaLxWKxWCwWSyFYo9lisVgsFovFYikEazRbLBaLxWKxWCyFYI1mi8VisVgs\nFoulEKzRbLFYLBaLxWKxFII1mi0Wi8VisVgslkLw2mgWkVdF5G8RWSMiX4tItXzadRKRjSKyRUQe\n9V7U0sfChQuDLUJAKWvXC/aaLcVDRGqKyHwR2Swi80Skej7t8tSz+fUXkY4islJE1rr/Xpmrz3ki\nss491pv+v8qSQ1n8bttrLv2UtestDsXxNM8D2qjqWcBmYMSpDUQkDHgb6AS0BnqLyP8VY85SRVn7\nopa16wV7zZZi8xgwX1VbAD+5X59EIXo2v/77gGtV9UzgduCTXEO+A/RX1eZAcxHp5PvLKpmUxe+2\nvebST1m73uLgtdGsqvNV1eV+uRxokEezC4CtqrpTVbOAacD13s5psVgsZYxuwBT38ylA9zzaFKRn\n8+yvqqtVda/7+AaggohEiEg0UEVVV7jPfZzPnBaLxVLm8FVM813A3DyOxwC7c72Ocx+zWCwWS+HU\nVdVE9/NEoG4ebQrSs570vxH4021wx7j75xCP1dkWi8UCgKhq/idF5gP18jj1uKp+624zEjhXVW/M\no/+NQCdVvcf9ui9woarel0fb/AWxWCyWEEdVxZt+BejZkcAUVa2Rq+1+Va15Sv9T9Ww/oK2qDhOR\n1IL6i0gb4Bugo6ruEJHzgRdVtaP7/GXAI6p6XR5yW51tsVhKLN7o7PBCBuxY0HkRuQPoArTPp0k8\n0DDX64ac7MXIPZdXPzgWi8VSkilIz4pIoojUU9W97tCJpDyanapnG7iPAeTbX0QaAF8D/VR1R66x\nGuQz1qlyW51tsVjKFMXJntEJeBi4XlXT82m2ErORpLGIRAK9gNnezmmxWCxljNmYjXq4/87Ko01B\nejbP/u4sGnOAR1V1ac5AqroHOCQiF4qIAP3ymdNisVjKHAWGZxTYUWQLEAnsdx9aqqqDRaQ+MFlV\nu7rbdQbeAMKA91X1xeKLbbFYLKUfEakJzAAaATuBm1X1gKd6toD+ozCZNLbkmq6jqiaLyHnAR0AF\nYK6qDvP7hVosFksJwGuj2WKxWCwWi8ViKSsEtCKgJ4VOROQt9/k1InJOIOXzB4Vds4j0cV/rWhH5\nXUTODIacvsTTgjYi0lZEskXkhkDK5w88/G7HisgqEVkvIgsDLKLP8eC7XU1EvhWR1e5rviMIYvoM\nEfnAHWO8roA2pUp/gdXbZUFvW51tdbb7vNXZhaGqAXlglg23Ao2BCGA18H+ntOmCWQ4EuBBYFij5\ngnjNFwPV3M87lYVrztXuZ+A74MZgyx2Az7k68D+ggft1rWDLHYBrfhyTiQGgFpAChAdb9mJc82XA\nOcC6fM6XKv1VhM+5VF13WdPbVmdbnZ2rjdXZhYwZSE+zJ4VOjifiV9XlQHURySuvaEmh0GtW1aWq\netD9Mr8iMSUJTwva3Ad8ialMVtLx5JpvBb5S1TgAVU0OsIy+xpNrdgFV3c+rAimqmh1AGX2Kqv4G\npBbQpLTpL7B6uyzobauzrc7OwersQnRXII1mTwqd5NWmJCujohZ36U/eRWJKEoVes4jEYP5Z33Ef\nKumB9Z58zs2BmiLyi4isFJNLtyTjyTW/DbQWkQRgDfDfAMkWLEqb/gKrt6H0622rs63OzsHq7EJ0\nV4F5mn2Mp/9kp+b+LMn/nB7LLiJXYiorXuI/cQKCJ9f8BvCYqqo7rVVJz/fqyTVHAOdicppXBJaK\nyDJV3VJwt5DFk2vuBPylqleKyH+A+SJylqoe9rNswaQ06S+wertASonetjo7b6zOtjr7XwTSaPak\n0ElBSfpLIh4Vd3FvIpmMqepV0FJCScCTaz4PmGZ0L7WAziKSpaolNYe3J9e8G0hW1TQgTUQWAWdx\ncsqvkoQn13wH8CKAqm4TkR1AS0xe4dJIadNfYPU2lH69bXW21dk53IHV2QXqrkCGZ3hS6GQ2cBuA\niFwEHFDVxADK6GsKvWYRaYSpytVXVbcGQUZfU+g1q2pTVW2iqk0wMXKDSrDyBc++298Al4pImIhU\nxGw62BBgOX2JJ9e8C+gA4I4TawlsD6iUgaW06S+werss6G2rs63OzsHq7EJ0V8A8zaqaLSJDgR85\nkYD/bxEZ4D7/rqrOFZEuIrIVOArcGSj5/IEn1ww8CdQA3nHfxWep6gXBkrm4eHjNpQoPv9sbReQH\nYC1ms8VkVS2xCtjDz/lZ4CMRWYtZAntEVffnO2iIIyJTgSuAWiKyG3gKs4RbKvUXWL1NGdDbVmdb\nne0+b3W2B7rLFjexWCwWi8VisVgKIaDFTSwWi8VisVgslpKINZotFovFYrFYLJZCsEazxWKxWCwW\ni8VSCNZotlgsFovFYrFYCsEazRaLxWKxWCwWSyFYo9lisVgsFovFYikEazRbLBaLxWKxWCyF8P+C\nMPkTZZiKKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a930a50>"
      ]
     },
     "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.002,0.002))\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",
    "    # Now subtract the unperturbed eigenvector to see the change from the perturbation\n",
    "    psi = psi - basis_array[m]\n",
    "    ax2.plot(x,psi)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the eigenvectors for this changed system are *very* close to the original system.  This is a perfect example of a good system to study with perturbation theory.  The changes in the eigenvectors can be related to matrix elements of the potential and the eigenvalues of the unperturbed system, but we won't do this just yet.\n",
    "\n",
    "In the next notebook, we will explore a more complex perturbation (a truncated quantum harmonic oscillator, or QHO) and the effects of varying the magnitude of the perturbation.  As well as improving your understanding of matrix mechanics, this should give a useful insight into the limitations of perturbation theory and finite basis sets."
   ]
  },
  {
   "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
}