TISE-Numerov.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical solution to the 1-dimensional Time Independent Schroedinger Equation\n",
    "Based on the paper \"Matrix Numerov method for solving Schroedinger's equation\" by Mohandas Pillai, Joshua Goglio, and Thad G. Walker, _American Journal of Physics_ **80** (11), 1017 (2012).  [doi:10.1119/1.4748813](http://dx.doi.org/10.1119/1.4748813)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.379434Z",
     "start_time": "2019-01-18T21:41:04.370438Z"
    }
   },
   "outputs": [],
   "source": [
    "# import some needed libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.optimize as opt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.395425Z",
     "start_time": "2019-01-18T21:41:04.384430Z"
    }
   },
   "outputs": [],
   "source": [
    "autoscale = False # set this equal to true to use Pillai's recommended step sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.429402Z",
     "start_time": "2019-01-18T21:41:04.399421Z"
    }
   },
   "outputs": [],
   "source": [
    "# values of constants\n",
    "hbar = 1.0\n",
    "mass = 1.0 # changing the mass will also change the energy scale\n",
    "omega = 1.0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.453387Z",
     "start_time": "2019-01-18T21:41:04.437398Z"
    }
   },
   "outputs": [],
   "source": [
    "# bounds (These are overwritten if autoscale=True)\n",
    "xmin = -5.0 # lower bound of position\n",
    "xmax = 5.0 # upper bound of position\n",
    "n = 100 # number of steps (may be overwritten if autoscale == True)\n",
    "dx = (xmax-xmin)/(n-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.478372Z",
     "start_time": "2019-01-18T21:41:04.460383Z"
    }
   },
   "outputs": [],
   "source": [
    "# the function V is the potential energy function\n",
    "def V(x):\n",
    "    # make sure there is no division by zero\n",
    "    # this also needs to be a \"vectorizable\" function\n",
    "    # uncomment one of the examples below, or write your own.\n",
    "    return 0.5*mass*omega**2*x*x # harmonic oscillator\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.504357Z",
     "start_time": "2019-01-18T21:41:04.481371Z"
    }
   },
   "outputs": [],
   "source": [
    "if (autoscale): \n",
    "    #Emax is the maximum energy for which to check for eigenvalues\n",
    "    Emax = 20.0\n",
    "    #The next lines make some reasonable choices for the position grid size and spacing\n",
    "    xt = opt.brentq(lambda x: V(x)-Emax ,0,5*Emax) #classical turning point\n",
    "    dx = 1.0/np.sqrt(2*Emax) #step size\n",
    "    # bounds and number of steps\n",
    "    n = np.int(0.5+2*(xt/dx + 4.0*np.pi)) #number of steps\n",
    "    xmin = -dx*(n+1)/2 \n",
    "    xmax = dx*(n+1)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.534338Z",
     "start_time": "2019-01-18T21:41:04.510353Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-5.0, 5.0, 100)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xmin, xmax, n #show the limits and number of steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.553327Z",
     "start_time": "2019-01-18T21:41:04.539335Z"
    }
   },
   "outputs": [],
   "source": [
    "#define the x coordinates\n",
    "x = np.linspace(xmin,xmax,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.594301Z",
     "start_time": "2019-01-18T21:41:04.557326Z"
    }
   },
   "outputs": [],
   "source": [
    "#define the numerov matrices\n",
    "B = np.matrix((np.eye(n,k=-1)+10.0*np.eye(n,k=0)+np.eye(n,k=1))/12.0)\n",
    "A = np.matrix((np.eye(n,k=-1)-2.0*np.eye(n,k=0)+np.eye(n,k=1))/(dx**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.626280Z",
     "start_time": "2019-01-18T21:41:04.597299Z"
    }
   },
   "outputs": [],
   "source": [
    "#calculate kinetic energy operator using Numerov's approximation\n",
    "KE = -0.5*hbar**2/mass*B.I*A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.640278Z",
     "start_time": "2019-01-18T21:41:04.632279Z"
    }
   },
   "outputs": [],
   "source": [
    "#calculate hamiltonian operator approximation\n",
    "H = KE + np.diag(V(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.700235Z",
     "start_time": "2019-01-18T21:41:04.647268Z"
    }
   },
   "outputs": [],
   "source": [
    "#Calculate eigenvalues and eigenvectors of H\n",
    "energies, wavefunctions = np.linalg.eigh(H) # \"wavefunctions\" is a matrix with one eigenvector in each column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.715226Z",
     "start_time": "2019-01-18T21:41:04.704233Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0.49999959,   1.49999715,   2.49998984,   3.49997482,\n",
       "         4.49995293,   5.49995153,   6.50014214,   7.50129661,\n",
       "         8.50607825,   9.52150074,  10.56141767,  11.64604745,\n",
       "        12.79708967,  14.03179647,  15.36046274,  16.78765495,\n",
       "        18.31462584,  19.94114087,  21.66644546,  23.48967816,\n",
       "        25.41001613,  27.42670997,  29.53907833,  31.74649046,\n",
       "        34.04834693,  36.44406207,  38.93304868,  41.51470496,\n",
       "        44.18840319,  46.95347979,  49.80922647,  52.75488201,\n",
       "        55.7896247 ,  58.91256508,  62.12273897,  65.41910068,\n",
       "        68.80051628,  72.26575701,  75.81349263,  79.44228487,\n",
       "        83.15058079,  86.93670626,  90.79885933,  94.73510369,\n",
       "        98.74336215, 102.8214102 , 106.96686956, 111.17720203,\n",
       "       115.4497033 , 119.78149707, 124.16952936, 128.61056316,\n",
       "       133.10117332, 137.63774192, 142.21645406, 146.83329415,\n",
       "       151.48404277, 156.16427421, 160.86935467, 165.59444129,\n",
       "       170.33448199, 175.08421624, 179.83817686, 184.59069278,\n",
       "       189.33589301, 194.06771177, 198.77989479, 203.46600697,\n",
       "       208.1194413 , 212.73342912, 217.30105171, 221.81525326,\n",
       "       226.26885512, 230.65457125, 234.96502493, 239.19276637,\n",
       "       243.3302913 , 247.37006004, 251.30451686, 255.12610911,\n",
       "       258.82730543, 262.40061208, 265.83858612, 269.13384327,\n",
       "       272.27905738, 275.26694657, 278.09023772, 280.74159516,\n",
       "       283.21348744, 285.49794108, 287.5860726 , 289.46708686,\n",
       "       291.12733678, 292.52000407, 293.79198764, 294.23037647,\n",
       "       296.4085974 , 296.41538741, 300.43135346, 300.4313561 ])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energies #display the lowest four energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:42.346083Z",
     "start_time": "2019-01-18T21:41:41.669500Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABcHElEQVR4nO3dd3xUVd7H8c+Zkt4bCYSShN57bwooCoioqyKoKCB2UVz7qqu7dll9sCJlWUCwggUVRbrU0GuoCWmk9zbtPH/cEAEpIcnkTpLz5jWvyfTvhMxvzj333HOElBJFURSl4TDoHUBRFEWpXarwK4qiNDCq8CuKojQwqvAriqI0MKrwK4qiNDCq8CuKojQwTiv8Qoh5Qoh0IcT+s64LEkL8JoQ4Wn4e6KzXVxRFUS7MmS3+/wIjz7vuGeB3KWUr4Pfyy4qiKEotEs48gEsI0QL4UUrZsfxyHDBUSpkqhIgA1kop2zgtgKIoivIXplp+vUZSylSA8uIfdrE7CiHuA+4D8Pb27tG2bdtaivinlJQUUlNT6dy5M2azudZfX1GUuiszM5OEhATatm2Lt7e3Lhl27NiRKaUMPf/62i78lSalnA3MBujZs6eMjY2t9QwHDx6kQ4cOTJ48mUcffbTWX19RlLpr+PDhmEwmDh48iBBClwxCiIQLXV/bo3rSyrt4KD9Pr+XXvyLt27enc+fOLF26VO8oiqLUIadPn2bNmjWMHz9et6J/KbVd+L8H7i7/+W7gu1p+/Ss2fvx4Nm/eTHx8vN5RFEWpI7766iscDgfjx4/XO8oFOXM45xJgM9BGCJEkhJgMvAGMEEIcBUaUX3Zpt912GwBffPGFzkkURakrli5dSqdOnWjfvr3eUS7IaX38UsqLfdUNc9ZrOkNUVBR9+/ZlyZIlPP3003rHUS7DarWSlJREaWmp3lGcwsPDg8jISDXYwIUlJCSwadMmXnvtNb2jXJTL7tx1JePHj+exxx7j0KFDtGvXTu84yiUkJSXh6+tLixYtXLJvtTqklGRlZZGUlERUVJTecZSLOLNP8Pbbb9c5ycWpKRsq4dZbb8VgMLBkyRK9oyiXUVpaSnBwcL0r+gBCCIKDg+vt1kx9sWTJEvr06ePSX86q8FdCeHg4V199NYsXL0atWOb66mPRP6M+v7f64MCBA+zZs4cJEyboHeWSVOGvpAkTJnDixAm2bt2qdxRFUVzU4sWLMRqNFYNCXJUq/JV000034eHhweLFi/WOoiiKC3I4HHz++eeMGDGCsLCLTkrgElThryQ/Pz/GjBnDF198gdVq1TuO4qKGDh3KypUrz7nuvffe48EHHyQ1NZXRo0df8vE//vgjL730kjMjKk6yadMmEhISXL6bB1ThvyITJkwgIyOD33//Xe8oiosaP378X470Xrp0KePHj2fmzJlMnTr1ko8fNWoU33//PcXFxc6MqTjB559/jqenJ2PHjtU7ymWp4ZxX4LrrriMwMJDFixczcuT5M04rrmb69Ons3r27Rp+za9euvPfeexe9/ZZbbuGFF16grKwMd3d34uPjSUlJYeDAgdx9993861//AmDmzJns37+fefPmsW/fPsaPH8+2bdvw8vJi6NCh/Pjjj9x66601ml1xHqvVypdffsnYsWPx9fXVO85lqRb/FXBzc+OWW25h2bJlFBUV6R1HcUHBwcH07t2bX375BdBa+7fddhvx8fEEBgbi7u4OaF9Kx44dY9myZdxzzz18+umneHl5AdCzZ082bNig23tQrtzKlSvJysqqE908oFr8V2zChAl89tlnfP/99y47D4eiuVTL3JnOdPeMHTuWpUuXMm/ePFJTUwkN/XN2XIPBwH//+186d+7MtGnTGDBgQMVtYWFhpKSk6BFdqaLFixcTFBTENddco3eUSlEt/is0aNAgmjVrxsKFC/WOorioG2+8kd9//52dO3dSUlJC9+7d8fT0/MuBV0ePHsXHx+cvRb60tBRPT8/ajKxUQ35+PsuXL+e2227Dzc1N7ziVogr/FTIYDEyYMIFff/2VtLQ0veMoLsjHx4ehQ4dy7733VmwVtm7d+pwZXvPy8njsscdYv349WVlZfP311xW3HTlyhI4dO9Z2bKWKvv32W0pLS7nzzjv1jlJpqvBXwcSJE7Hb7WoKB+Wixo8fz549eyrma/H29iYmJoZjx44B8Pjjj/Pggw/SunVr5s6dyzPPPEN6urY8xZo1axg1apRu2ZUrs3DhQlq2bEnfvn31jlJpqvBXQfv27enevbvq7lEuaty4cUgpOXvJ0Icffpj//ve/AMybN69iVbemTZty7NgxwsLCSEtLo6SkhE6dOukRW7lCSUlJrFmzhokTJ9ap6TRU4a+iO++8k507d3Lw4EG9oyh1xLhx42jRosUl73Pq1Cnefffd2gmkVNuZ+bsmTpyod5Qrogp/FY0fPx6j0ciiRYv0jqLUIVOmTLnk7b169aJr1661E0apFiklCxcupH///sTExOgd54qowl9FjRo14pprrmHx4sU4HA694yiKUst2797NgQMH6tRO3TNU4a+GO++8k1OnTrF+/Xq9oyiKUssWLVqE2Wyuk0dYq8JfDWPHjsXHx4f//e9/ekdRFKUW2Ww2Pv/8c0aNGkVQUJDeca6YKvzV4OXlxa233spXX32lpnBQzvHLL7/Qpk0bWrZsyRtvvKF3HKWG/frrr5w+fZpJkybpHaVKVOGvprvvvpvCwkKWLVumdxTFRdjtdh566CF+/vlnDh48yJIlS9Tor3rmv//9LyEhIVx33XV6R6kSVfiraeDAgURFRVWMz1aUbdu20bJlS6Kjo3Fzc+P222/nu+++0zuWUkNycnL47rvvuOOOO+rMFA3nU5O0VZPBYOCuu+7ilVdeITExkaZNm+odSSn3zx8OcDAlv0afs31jP14a0+GS90lOTj7n7yAyMlIt2VmPfPHFF1gsFu6++269o1SZavHXgLvuuqtiTK+iSCn/cl1dOqpTubQFCxbQsWNHunXrpneUKlMt/hoQHR3N4MGDWbBgAc8++6z6kLuIy7XMnSUyMpLExMSKy0lJSTRu3FiXLErNiouLY8uWLbzzzjt1+nOuWvw15O677+bIkSNqk16hV69eHD16lJMnT2KxWFi6dCk33HCD3rGUGrBgwQKMRmOdWXDlYlThryG33HILnp6eaievgslk4oMPPuDaa6+lXbt23HrrrXTooM/Wh1Jz7HY7Cxcu5NprryU8PFzvONWiCn8N8fPz45ZbbmHJkiVqoWyF66+/niNHjnD8+HGef/55veMoNWDVqlUkJSVxzz336B2l2lThr0H33HMP+fn5aky/otRD8+fPJygoiDFjxugdpdpU4a9BQ4YMISoqivnz5+sdRVGUGpSdnc2yZcuYOHEi7u7uesepNlX4a5DBYGDSpEn8/vvv5yyzpyhK3bZkyRIsFku96OYBnQq/EOJxIcQBIcR+IcQSIYSHHjmc4e6770YIwYIFC/SOoihKDZk3bx7dunWrN2sl1HrhF0I0AR4FekopOwJG4PbazuEszZs3Z9iwYcyfP1/N068o9cCePXvYuXNnvWntg35dPSbAUwhhAryAFJ1yOMW9995LQkICa9as0TuKoijVNH/+fNzc3Ljjjjv0jlJjar3wSymTgXeAU0AqkCel/PX8+wkh7hNCxAohYjMyMmo7ZrXceOONBAQEMG/ePL2jKDq49957CQsLo2PHjnpHUaqprKyMRYsWMXbsWIKDg/WOU2P06OoJBMYCUUBjwFsI8ZeViqWUs6WUPaWUPUNDQ2s7ZrV4enpyxx138M0335CTk6N3HKWWTZo0iV9++UXvGEoN+O6778jKymLy5Ml6R6lRenT1DAdOSikzpJRW4Fugvw45nGrKlCmUlZWxePFivaMotWzw4MF1clUm5a/mzJlDs2bNGD58uN5RapQek7SdAvoKIbyAEmAYEKtDDqfq1q0bPXr04LPPPuOhhx6q0xM61Vk/PwOn99Xsc4Z3guvUiloNwcmTJ/ntt9/45z//idFo1DtOjdKjj38r8DWwE9hXnmF2beeoDVOmTGHv3r3s2LFD7yiKolyh+fPnI4SoV6N5ztBlWmYp5UvAS3q8dm0aP348TzzxBHPmzKFnz556x2l4VMtcqSK73c68efMYOXJkvVxcSR2560T+/v7ceuutfP7552oxdkWpQ1auXElycjJTpkzRO4pTqMLvZFOmTKGgoICvvvpK7yhKLRk/fjz9+vUjLi6OyMhI5s6dq3ck5QrNmTOHsLAwRo8erXcUp1ArcDnZgAEDaNOmDbNnz2bSpEl6x1FqwZIlS/SOoFTD6dOn+eGHH5g+fXqdXUz9clSL38mEEEydOpXNmzezf/9+veMoinIZ8+bNw2azMXXqVL2jOI0q/LXg7rvvxs3Njdmz6+XgJUWpNxwOB5999hlXXXUVrVu31juO06jCXwtCQkK45ZZb+N///qdW51IUF/bbb78RHx/PtGnT9I7iVKrw15Jp06aRl5fHl19+qXcURVEu4tNPPyUkJIQbb7xR7yhOpQp/LRk0aBBt27ZV3T2K4qJSU1P5/vvvueeee+rFKluXogp/LRFCcN9997F582b27avhaQQURam2efPmYbfb6/VO3TNU4a9Fd999N+7u7nz66ad6R1GcKDExkauuuop27drRoUMH3n//fb0jKZdht9v57LPPuPrqq2nVqpXecZxOFf5aFBQUxN/+9jcWLlxIYWGh3nEUJzGZTLz77rscOnSILVu28OGHH3Lw4EG9YymXsHLlShISEur9Tt0zVOGvZQ888AD5+fnqIJ96LCIigu7duwPg6+tLu3btSE5O1jmVcikff/wxjRo1qvc7dc9QR+7Wsn79+tGlSxc++ugjpkyZoqZrdqI3t73J4ezDNfqcbYPa8nTvpyt9//j4eHbt2kWfPn1qNIdSc+Lj41mxYgXPP/98vT1S93yqxV/LhBA88MAD7N69m61bt+odR3GiwsJCbr75Zt577z38/Pz0jqNcxOzZsysGXzQUqsWvgwkTJvD3v/+djz76iL59++odp966kpZ5TbNardx8881MmDCBm266SbccyqWVlZUxZ84cxowZUy+nX74Y1eLXgY+PD3fddRdffPEFmZmZesdRapiUksmTJ9OuXTueeOIJveMol/DNN9+QkZHBAw88oHeUWqUKv04eeOABLBYL8+fP1zuKUsP++OMPFi5cyOrVq+natStdu3blp59+0juWcgEff/wxMTExjBgxQu8otUp19eikQ4cODBo0iI8//pgnnnii3q3p2ZANHDgQKaXeMZTL2Lt3Lxs3buTtt9/GYGhYbeCG9W5dzEMPPcTJkyf55Zdf9I6iKA3Ohx9+iIeHR71cU/dyVOHX0U033URERASzZs3SO4qiNCg5OTksWrSICRMmEBwcrHecWqcKv47MZjP3338/K1eu5MiRI3rHUZQGY/78+RQXF/Pwww/rHUUXqvDr7L777sNsNvPhhx/qHUVRGgS73c6HH37IwIED6dq1q95xdKEKv87Cw8P529/+xn//+18KCgr0jqMo9d4vv/zCiRMnGmxrH1ThdwkPP/ww+fn5LFq0SO8oilLvffDBB0RERDToA+tU4XcBffv2pUePHsyaNUsNA6wHSktL6d27N126dKFDhw689NJLekdSyh05coRffvmF+++/H7PZrHcc3ajC7wKEEDzyyCMcOnSIVatW6R1HqSZ3d3dWr17Nnj172L17N7/88gtbtmzRO5YCzJo1C7PZ3KDm5bkQVfhdxO23305YWJhatKMeEELg4+MDaHP2WK1WNQurC8jNzWX+/PmMHz+e8PBwvePoSh256yLc3d154IEH+Oc//8mRI0do3bq13pHqvNOvvUbZoZqdltm9XVvCn3vusvez2+306NGDY8eO8dBDD6lpmV3AvHnzKCoq4rHHHtM7iu5Ui9+F3H///bi5uakDuuoBo9HI7t27SUpKYtu2bezfv1/vSA2a3W5n1qxZDBo0qGKRnIZMtfhdSHh4OLfffjvz58/n1VdfJSAgQO9IdVplWubOFhAQwNChQ/nll1/o2LGj3nEarO+//574+HjeeecdvaO4BF1a/EKIACHE10KIw0KIQ0KIfnrkcEWPPfYYRUVFzJs3T+8oShVlZGSQm5sLQElJCatWraJt27b6hmrg3n//fZo3b87YsWP1juIS9OrqeR/4RUrZFugCHNIph8vp3r07gwYNYtasWdjtdr3jKFWQmprKVVddRefOnenVqxcjRoxg9OjResdqsHbv3s26det4+OGHMZlUJwfo0NUjhPADBgOTAKSUFsBS2zlc2fTp07n55ptZtmwZt9xyi95xlCvUuXNndu3apXcMpdzMmTPx9vZm8uTJekdxGXq0+KOBDGC+EGKXEGKOEML7/DsJIe4TQsQKIWIzMjJqP6WOxo4dS0xMDDNnztQ7iqLUacnJySxZsoTJkycTGBiodxyXoUfhNwHdgY+llN2AIuCZ8+8kpZwtpewppewZGhpa2xl1ZTQamT59Ops3b2bz5s16x1GUOuuDDz7A4XAwffp0vaO4FD0KfxKQJKXcWn75a7QvAuUs99xzD4GBgbz77rt6R1GUOqmwsJBPPvmEm266iaioKL3juJRaL/xSytNAohCiTflVw4CDtZ3D1Xl7ezNt2jSWLVvGiRMn9I6jKHXO/Pnzyc3NVQveX4Beo3oeARYLIfYCXYHXdMrh0h555BGMRqOaxkFRrpDdbue9996jX79+9OunRoufT5fCL6XcXd5/31lKeaOUMkePHK6ucePGjB8/nrlz55KTo35FilJZ3333HSdOnGDGjBl6R3FJasoGFzdjxgyKior4+OOP9Y6iXCG73U63bt3UGP5aJqXkzTffJDo6mhtvvFHvOC5JFX4X17lzZ0aOHMn7779PaWmp3nGUK/D+++/Trl07vWM0OBs2bGDbtm08+eSTGI1GveO4JFX464Cnn36a9PR0FixYoHcUpZKSkpJYsWIFU6ZM0TtKg/PWW28RGhrKpEmT9I7istTxy3XAkCFD6NWrF++88w5TpkxRrZhK2vDlETITC2v0OUOa+jDo1stPmT19+nTeeusttY5yLdu/fz8rVqzg1VdfxdPTU+84Lku1+OsAIQRPPfUUx44dY9myZXrHUS7jxx9/JCwsjB49eugdpcF5++238fLy4sEHH9Q7iktTLf46Yty4cbRs2ZK33nqLm2++Wa3oVAmVaZk7wx9//MH333/PTz/9RGlpKfn5+UycOJFFixbpkqehSExM5PPPP+ehhx4iKChI7zguTbX46wij0ciTTz7J9u3bWbNmjd5xlEt4/fXXSUpKIj4+nqVLl3L11Verol8L/vOf/yCl5PHHH9c7istThb8OufvuuwkPD+f111/XO4qiuJTMzEw+/fRTJkyYQPPmzfWO4/JU4a9DPDw8eOKJJ1i1ahXbtm3TO45SCUOHDuXHH3/UO0a993//93+UlJTwzDN/me9RuYBKFX4hxDtCiA7ODqNc3v33309AQIBq9StKufz8fGbNmsW4cePUcROVVNkW/2FgthBiqxDifiGEvzNDKRfn6+vLo48+yvLlyzlw4IDecRRFd5988gm5ubk8++yzekepMypV+KWUc6SUA4C7gBbAXiHE50KIq5wZTrmwRx99FG9vb9544w29o7gkKaXeEZymPr+3qigpKWHmzJlcc8019OzZU+84dUal+/iFEEagbfkpE9gDPCGEWOqkbMpFBAcHM23aNJYsWaKmbD6Ph4cHWVlZ9bJASinJysrCw8ND7yguY/78+aSlpfHcc8/pHaVOEZX5gAghZgI3AL8Dc6WU2866LU5K2eaiD64BPXv2lLGxsc58iTonJSWF6Oho7rrrLmbPnq13HJdhtVpJSkqqt/MaeXh4EBkZidls1juK7iwWCy1btqRp06Zs3LhRHdtyAUKIHVLKv2wKVfYArv3AC1LK4gvc1rtayZQqady4MVOmTGH27Nk8//zzDWIIm93u4PSxPLJTiygptFJSYMFgEAQ08iIgzIvgSB+8/Nzq9WpLUkpy04rJSc0lN72Y/MwSTG5GPH3NePq6ER7lT2CEV4MoggsWLCAxMZE5c+Y0iPdbkyrb4r/Q0oh5QIKU0lbjqc6jWvwXlpiYSExMDFOmTOGjjz7SO45TSIckfl8mR2PTOXUgi7LiP//c3L1NOGwSa5kdACGgWYdg2vaLIKpzCEZz/RmtXJxv4ci20xzenEpWclHF9R7eZmw2B7by3wGAX6gnUZ1DaNc/guAmPnrEdTqr1Urr1q1p1KgRmzdvVoX/Ii7W4q9s4d+Cti7uXkAAHct/Dgbul1L+WrNxz6UK/8Xdf//9zJ8/n+PHjxMZGal3nBojHZLjuzKI/ekkWclFePqaad4xmKjOoTSK9sPDx4zRaEBKSXG+hdzTxSQeyiZu62kKc8rw9HOj79ho2vaLwGCou0WhtMjK9hUn2b82GYdD0ijKjzZ9wmkU5Yd/mBfuntpGu9Vipyi3jKTDOZzck0lSXDYOmyS6Wyg9r29BaFNfnd9JzZo/fz733nsvK1as4Prrr9c7jsuqbuFfCrwqpTxQfrk98HfgVeBbKWXXmo17LlX4Ly4+Pp5WrVrxwAMP8H//9396x6kRp0/msXZRHFnJhQSGe9Hjuha06hmGwXj5FrzDIUk6nE3sinhSj+dVzKbZuFWA84PXIIdDcnBDMlu/P0lZsZV2AxrT5eqmBDX2rtTjS4us7FmdyN7VSVhKbLTq1YiBf2uFl5+bk5M7n81mo23btgQEBLB9+3bV2r+E6hb+3ecX9zPXXei2mqYK/6VNnjyZxYsXc/LkSSIiIvSOU2WWUhtbvzvB3rVJ+AS40/+mlsT0CKtSi11KybHYdDZ9e4zCnDK6jmhG37HRGE2u3/1TmFPKb/MOknI0lyZtAhj4t1aERFatxV5WbGX3qkR2/pqA2d3IwFta0aZveJ0ulgsXLuSuu+5i+fLljB07Vu84Lq26hf9LIAs4M3TzNiAEuBPYKKXsVYNZ/0IV/ks7fvw4bdq04eGHH+a9997TO06VnD6Rx69zD1CQXUqnIZH0HRuNm2f1J4+1Wuxs+voY+9cnE9bclxGTOxAQ5lUDiZ3jxO4MVi88hN0mGXJ76xor0tmpRaxddJjU43k06xDM8Ent8PSte61/m81G+/bt8fLyYteuXXX6C6w2VLfwewIPAgPR+vg3Ah8BpYCXlLJmV7s4jyr8l3em1X/8+HGaNGmid5xKkw7Jzl8T2Pr9SXyD3Bk+qT0RLQNq/HWO70xnzaLDOBySEfd2IKpzSI2/RnU4HJIty46z67dThDT14dopHQloVLNfUNIh2bcuiU3fHMfd28SIezsQ2SawRl/D2RYsWMCkSZNUa7+Sqlz4yw/cWimlHO6scJejCv/lnTx5ktatWzNt2jQ++OADveNUSkmhhVXzDnLqYDYx3cO46s62FTsrnaEgu5SfP9lHRmIB/W6Mods1zVyixVhWYuO3uQdI2J9Fh8FNGPS3Vk4dkZSZVMDKzw6Qm15Mr+tb0HNUVJ3YAW61Wiv69mNjY13i/87VXazwX/avS0ppB4rV/DyuLSoqinvvvZfPPvuMU6dO6R3nstLi8/ny39tJPpLLkDvacO3UDk4t+gC+QR6Me7I7LXuEsXnZcVbNP4jNar/8A50oN72Yb96MJfFgNkPGt2boHW2cPgw1JNKXvz3bkzZ9wtm+Ip4VH+6htMjq1NesCQsWLODEiRO88sorquhX05X08fcFfgMqBhFLKR91XrQ/VbXFf/CPFBIPZSMAhEAIQGhLGYqzzw1C+9kgMBgEwqidG4zayWgyVJyb3AwYzQZMZiNmNyMmNwNmDxNuHkbcPEy4eRorNfrEGU6dOkXLli259957+eSTT3TJcDlSSg5uTGH9F0fw9nNn5LSOhDX3q/UMO36OZ+v3JwmP9uO6+zvrMtol5WgOP32yD4CR93Wq9W4XKSUHNqSw4csjePu7c920ToQ2c81hnxaLhVatWhEeHs6WLVt0K/w2qx1LiR1LiQ1rmR1rmQ2rRTuOwmZ1YLc6tHObA4fdgcMuy09//iwdEocsn3ep/GccEilBIpEO7bWkdgU9rmte5Z371T1yd0X5qU4pyi0jM7Gw4hd45ktOSq2/E7S+VXnml+6QFZcdDu0/iSpM+eLmacLD24SHtxkvPzc8fd3wDnDHO8Adn0B3fIM88Av1xOxWs4umN2vWjKlTpzJ79myefvpplzuC1Waxs27pEQ5vSqVZ+yBG3NsBD5/an3pACEHP66MIDPdm1fyDfP1GLKMe6lyrBzsd3pzKmkWH8QvxZNRDnXXZ4SyEoOPgJoQ29eWX2fv45q0dDB7fmvYDGtd6lsuZN28ep06dYvbs2TVe9KWUlBXZyMssoTC7lMKcMgpzyyjOL6Mk30JxgZXSQitlRVZsVscVP78QYDAaKhqUwoDWtSYEhvKG59kNUs5qmAJYSmp+q7RSLX4tvPAEmkkp42o8xWXo2cfvcEgcNu0b3G6T2Kz2im91W5kdq8WOtezPVkBZiY2yIiulxVZKC6wUF1gozrdQkm/h/F+1l78bAWFeBEV4ExjhTXBjb0Ka+VaryyM5OZmYmBjGjx/P/Pnzq/nua05eRjE/f7qfrKRCeo5qQS8X6VdOT8hnxUd7sZbZGXFPe6K6hDr19Rx2B1u+O8GuX08R2TaQa6d2xMNb/3l3Sgos/Dr3AEmHc2jXP4LBt7fGVMMNk6oqKSmhZcuWtGjRolpz8kgpKcwpIzOxgOzUIrJTi8hJLSYvvRhL6bnF1Wg24OXnVtFw8/Qx4+5lwt1LO3fzNGF2N+LmYcTkZsTsbjynR8BoNFT0FAgd/86rO6pnDPAO4CaljBJCdAVekVLeUONJL6A+7Nx12B0U51sozCkjP6uE/IwS8jJKyE0rJju1GEvJn1MR+Id6Etrcl/Aof8Jj/Alp6oPxCrqPZsyYwXvvvcf+/ftdYmGKYzu0ETVCwPB72tOik2uNqCnMKeWnj/eRcaqAnqNa0HtUlFM+rCWFFn6doxXXjoObMPC2Vlf0/+psDodk2w8n2PFzAiFNfbhmcgcCwyt3wJgzvfvuuzz55JOsXbuWIUOGVPpxZcVW0k7mk3oij7ST+WScKqC08M99GT6B7gRGeBMQ6olfqCd+IZ74BnvgG+iBu7epXuxHqG7h3wFcDayVUnYrv26flLJTjSe9gPpQ+C9FSklxnoXMpEIyEgvIPFVAWnw+hTllAJjMBiJa+tOkTSBN2gQS1tzvkq3lzMxMoqKiuPbaa/n6669r6238hbXMzoYvj3Doj1TCWvhx7ZQO+IV46pbnUmwWO+uWxHF482madwpm+KT2NdoST0/I55dP91OUX8aQ8W1csjvljPh9mfz+30PYrHYG3dqadgMidCuC+fn5REdH06NHD1auXHnJ+1pKbaQczSUpLofkuBwykwpBal0tQY19CGvuS2gzX0Ka+hLU2NvpgwlcQXUL/1YpZR8hxK6zCv9eKWVnJ2T9i/pe+C+mMKeU0yfySTmWS3JcDtkp2n51d28TzdoF0axDMM07Bl/wQJyXX36Zf/7zn2zfvl2XBSpSj+ex+n+HyE0vpse1zek1JsqlWrcXIqVk/7pkNn55FA8fM0MntKl214/d5iD253h2/JyAt78bI+/rRKOo2t2ZXRVFuWWs+u9Bkg7nENMtlEG3t8bb373Wc5z5O46NjaVHjx7n3CalJOd0MQn7skg4kEXqsVwcdonBJIiI9qdx60AiYvxp1MKvRg4GrIuqW/jnos3F/wxwM/AoYJZS3l/TQS+koRb+8xXnW0iKyybxQDYJB7MpybeAgPAof6K6hBDVJaRi0/xKWko1qazYyublJziwIRmfAHeGTWpf5w4Sykgs4Pf/HiIruZA2fcLpf3PLKo36SYvPZ82iw2Qlac8z8NZWLtGfX1nSIdn12ym2/nACs5uRfuNiaD+gca31WWdkZBAdHX3OlqvD7iD1WB4n92Zycm8m+RklAAQ19qZ5h2CadggiItrfZfZP6K26hd8LeB64Bu3I3ZVok7bVymoXqvD/lXRIMhILiN+bSfy+LDJOFQAQ0Mir/EsglCXL5zLjyRmsWbOGoUOHOjWP3e4gbstptn5/gpJ8C52vakrvG6Jw86ibLa2zW+pGo6D9oMZ0G9Ecn8DLt3pTj+cR+1M8pw5k4eXnViNbDnrKTStm7eLDJB/JJTzan343xdDYCUdXn++JJ57g/fffZ8+uvXjaQzm5J5P4fZmUFdkwmgw0aRNIVOdgmncKwTdIrUp2IdUq/M5QfkRwLJAspRx9qfuqwn95BdmlxO/N5OSeDJLjcnE4JJ6+ZjYf+I1s6ymWfP8Zbp4139p02B0cjU1n248nyc8oIayFH0PGt671sfnOkptWzI5f4onbmoYwQGSbIJq0CSCyTSDeAdqXgHRAVkohyXE5JB3OIeNUAR4+ZroMa0qnoZH1oi9ZSsnhzalsXq59sTdrH0TvMdFO67Y6tPcY90+YwTX9bibYvRl2qwN3LxMtOmlbtk3bB9XZRkVtqm6LvzXwJNpC6xW/bSnl1dUI9ATQE/BThb9mlRVbSTiQxcndmRzbfRrsBhCSpm2DiGwXRJPWgYQ29bmiA83submUHjpE2bHjlB07Rk5iLgmWJiS5t8Zi9sWnOIWY1FWElZ7AHBGBW2Qk5mZN8ereHc/uPTD66D86pDryM0vYuzqJUwezyDl9oYXowGAUNIryI6ZbGO0HNsbsXv+6G6wWO/vWJrFr5SlKi6yENvOlXf8IWvVqVK1uLEupjdRjeSTH5XDqYDZZydr0X57+Jlp1DyeqSwgRrQJcfj+Rq6lu4d8DfALsACoGvEopd1QxTCSwAPg38IQq/M5jKbNyw7DbiPBpxfC+48gtL1pmDyONWvgR0tSXkEgfgiK88Q5wx9PHjDAIbJmZFG3eTNHWbeTvOUhBSi6F3hHkBrQiN6gNxR5hCBw0cssiOiCXCI8ccNhxlBRjTUnBmpSMNSUF7HYwGvHs2BG/66/Db8wYTEFBOv9Wqqcor4yUo7nnrAbmH+JJeIz/FRd7e2ER1pRkbdSMwYDBywtTuOtPm2wpsXFocyqHNqWSlVSofem18KNxqwAiWgbgH+aJT6A7JvO5vw+73UFRbhmFOWVkJxeSkVhIxqkCMpMKkQ6JwSjwCTOyYNmH9B3RkVffesHlfxeurNrDOaWUPS57x8qH+Rp4HfAFnrxQ4RdC3AfcB9CsWbMeCQkJNfXyDc5PP/3EqFGjmDVrFvfcOZWUo7kkH8klIyGfrOQi7LY/j0Y0CAcmWynYbSAldpM7duOf/admdwMRrQJo2jaIVr0aXXKkh6OkhJLduynaupWiDRspPXAATCZ8hg4h+J578OpRY39SdYajtJSCX3+lcN16Sg8cwJKQwPlH9hkDAvDo1AnPrl3wHzUKtxYt9AlbSRmnCjgam0bK0VzSEwoqjooHbWlIg1EgpTYVQWmx9Zyj4d29TIQ09SUixp/GrQMIj/Zn1Ojr2LFjB8ePHycgIKD231A9Ut3C/zKQDiwDys5cL6XMrkKQ0cD1UsoHhRBDuUjhP5tq8VePlJJhw4axf/9+jh07hp+f1i9rLyyiYP0GUn7fRva+k5TY3Sn1CEKEN8UUHo6xUTjuYUH4BmuttzMLmlf1iNvSI0fIW7acvO++w56djVe/voQ+9BBeOgw3rW1lJ06Ss2gReT/8gKOgAFNYGB6dO+HRvj3uUVHaYHOHA3t+PqUHDlCydx9lR4+Cw4FXr14E3Po3/K67DmFy7X5tS6mNjIQCCrJLKcwppTDXgpQSgTY1gaevGz6B2vQlgeFe+AZ5nNOi//XXX7n22muZOXMmjz/+uH5vpJ6obuE/eYGrpZQyugpBXkdbwMUGeAB+aMs3TrzYY1Thr77Y2Fh69erFvx5/nAd79aJg9RqKt2xBWq0YAwLwGTIEn6uuwntAf4y+zp2oy1FSQs4XX5A1Zy72zEx8rr6aRs89h1tk3VlHoLJs2dlkfvABOV98iTAa8b32WgJuuQWvXj0Rhkv3V1vT0slbvpzcr7/GmpiIW1QUYTOewGfYsHrZ/WG32+nRowd5eXkcPnwYd/faP26gvnG5UT0AqsXvfNJup2TPXgrXr+PAgv8RXqKNezY3b4bv1cPwHXY1nl276tKSdJSUkL1oEZkffwIOByH3TyPo3nsxuNW9laHOJx0OchYtJuP//g9HSQkBt/6N0IcfxhQcXKXnKly9mvSZ/8Fy4gSe3bsT/uI/8Gjb1gnJ9TNv3jwmT57MkiVLuP322/WOUy9UqfALIZ6SUr5V/vPfpJRfnXXba1LK56oZaiiq8NcoKSXWxESKNm3Wds5u2YIjLw+MRgzt2vGfDetxHziI977+ymVajdbUVNLeeJOClStxb92axm+9WaeLmjUlhZRnn6N461a8Bw6k0bPP4B4TU+3nlTYbud9+S8b7/4c9L4+QadMImXYfoh58URYWFtKqVStatGjBpk2bXOZvs667WOHXdrpc5ATsvNDPF7rszFOPHj2kcmEOq1WWHDwos5cskUkznpRHhgyVB9u0lQfbtJVHhgyVyc88K/NWrJC23FwppZTPPfecBOTWrVt1Tv5X+WvWyLiBA+XBjp1kxqezpcNm0zvSFcv9/nt5uEdPebhbd5nz1VfS4XDU+GtYs7Nl0t//Lg+2aSuPj7lBlhw6VOOvUdteeOEFCcjNmzfrHaVeAWLlBWrq5Vr8u+Sfc/NU/Hyhy86kWvwaR1ERZSdOUHr4MGWH4yg9fJjSAweQpdoB1KbQULx69cSzZ0+8+/bFLSrqLy2ngoICWrduTXR0dLWmuHUWW04Op1/+JwUrV+LVsyeN330Hc6NGese6LIfFQtprr5G79As8u3en8Ztv4Na0qVNfs2DNGk6/+BL2vDzCnnqKwAl3uNz/Z2UkJibSunVrxo0bx+eff653nHqlql09O6WU3c//+UKXnamhFH5psWDLyMCalo7tdCqWpGSsSUlYEk9hORmP7fTpivsavLxwb9MGj04d8ezcBc8unTFHRlbqgz937lymTJnCF198wa233urMt1QlUkryvvuO06+8isHdncZvvoHP4MF6x7ooa3IySdMfp3TfPoKnTCZ0+vRa22diy84m5dlnKVq3Hp9hw2j8739hrGNDICdOnMjXX39NXFwczZs31ztOvVLVwm9HW2pRAJ7AmUMWBeAhpayVGadcpfBLKcHhAIcDabcjbXawWZE2m3ayWpEWC7KsDEdZmXZeUoIsLcVRXIyjsBB7URGOgkLs+Xk48vOx5+Riy8nBnpWFPTf3L69pDAzE3LQp7lFRuEVF4RYdhUfbtlqRv8yokIs5M3oiJyeHQ4cO4eVV+6s/VUbZiRMkT3+csiNHCJ46hdDHHnO54YxFW7aQPP1xpNVKxBuv4zdiRK1nkFKS87//kfbOu5hCQ4j8z3/w7NKl1nNUxaZNmxgwYADPPvssr732mt5x6h2XHNVTWVUt/Onvvkvut8v+vELKPw+WOetneeayw1FxvTxzu8Oh/Wy3a7fXAOHhgdHPD6O/HwZ/f0zBIRiDAjEFh2BqFIa5USNMjRphbhLptKkO1q9fz5AhQ3jppZd4+eWXnfIaNcFRWkraa6+T++WXePbsQZN3Z2JuFKZ3rD+L7Vtv49aiBZEfzNLG4+uoZN8+kqc/jjU9nUZ//zuBd0506a4fu91Onz59SE1NJS4uDh+f2lv6ErQRb/bcXOx5edjzyhtihYU4iopwFBXjKClGlpTiKC1FljfkpNVS3sCzauc2G9Juh/Jz6bCD3QEOu3Yg25mG4pl6cqamUL6s6wXq0fkH9DV++228+/ap0ntskIU/78cVFMduP/dKIc76MGjrXgJgMFSsewmi/LJAGA1oi2QKhMEIRoPW0jaZEEYTwmREmM3aZbMZg5sbwt0dUX5u8PTE4OGB8PTC6OONwdtbu78LGD9+PMuXL+fQoUO0cPGjQ/N++IHUl17G4OFBk3fexrt/f92yOEpKOP3yy+R99z0+w4fR+I03XWYuInteHinPPEvhmjX4jhhBxL//hdHPNSfMmzNnDlOnTmXx4sXccccdNfrc0mbDmpqKNTERS1IS1uQUbGlp2NLTsKanY8/K1rawL1f/zGYM7u4IDw/ts33mZDYjyj/zWi04UxuMYDRqNcJoLJ/C+s96otWYMz+LiuuAv9amcoF3TsSjdesq/R4aZOFXLi0xMZG2bdty3XXX6bpSV2WVHT9O0mOPYTl+guCpUwl95OFa/xK1JCSQ9OhjlB05QsjDDxHywANV7nJzFikl2fPmk/6f/2Bu1Igm/5mJZ+daWTOp0nJzc2nVqhVt27Zl/fr11VpH15qUpA14OHSYsqNHKTt5AmvCKaT1z2UWMRoxhYZqW9RhYRhDQjAFBmEMCsLo74/R3w+jnx8GHx/t5O2tNdhcpJFWVarwKxf073//mxdeeIFVq1YxbNgwveNclqO4mLTX3yD3q6/w7NKFxu++W2tH/Bb8/jspzzyLMBho/M7b+AwaVCuvW1Ulu3eT/MQMrOnphE1/jKB77tFapi5g+vTpzJo1ix07dtC1a9dKP85eWETJrl2U7NpJyd59lOzbpx2nAiAEbs2b4xYTg3t0FG4tWmCObIpbZBNMjRq53P6h2qAKv3JBpaWldOjQAXd3d3bv3o1bHTkYKP/nn0n9x4sgJWF//zsBt93qtP5sR1ERaW+9Te4XX+DRoQNN3n+/zkwvYc/LI/WFf1Dw22949uxB4zfewC0yUtdMe/fupXv37kyZMoVPPvnkkvd1lJRQHLuDos2bKd62jdKDB7V+c4MB91at8OzcGY+OHfFo1xb3Vq0weLrmms56UYVfuagff/yRMWPG8MYbb/D000/rHafSLEnJnH7xHxRt2oxX375EvPpKjY+dL965i5RnnsGamEjQPfcQ+tijGOrYHDJnhsem/evf4HAQ9tRTBNz6N126qBwOB4MGDeLIkSMcPnyY4POmsJBSYjkZT+H6dRStX0/x9lik1Yowm/Hs0gWv3r3w6tULz86dMXi7xn4VV6YKv3JJ48aN49dff+XgwYN1aiy1lJLcr74i/c23kFYrgXdOJOS++zD6+1frea1paWS89z55y5djbtyYxm+8jlevXjWUWh/W5GRSnn+B4i1b8OjcmfAXX8SzY4dazTB//nzuvfde5s2bxz333AOAtFop3rGTwjWrKVizFuupUwC4t2qJ94CBeA/oj1ePHhhcdNixK1OFX7mkhIQE2rdvzzXXXMOyZcsu/wAXYz19WivU332Hwc+P4MmTCbjl5ite9MWWmUnO55+TNf+/YLMReOedhDz4AMZaHmroLFJK8n/8kbQ338KelUXALbcQcv80zE2c33WVlZVFmzZtaNu2LWt//pniP/6g4PfVFK5bhyM/H+Hmhle/vvgOHYrP4MG1kqm+U4Vfuaw333yTZ555hh9++IHRoy85b94VsTvsJBYkklKUQlpRGqeLT1NiK8Fit2CxW/AweeDr5oufmx9NfZvSMqAlEd4RVeqzLz18mPR33qVo40Ywm/EdPgz/sWPx6tbtolsB9oICimNjyfv2WwrWrAWbDd+RIwmb8YTTp12oLVa7lfj8eI7lHiO9OJ2i3AyafvEHUauPICQkD25N+s0DCYxuS4R3BJG+kYR6htbofpOn7rqL0ytW8PSw4RgOHQKrFWNgID5Dh+Jz1VB8BgxQ3Tc1TBV+5bIsFgvdunWjqKiIAwcO4F2FD6GUklMFp9iZtpM9GXs4nH2YY7nHKLOXnXM/d6M7bkY3zAYzpbZSim3nrmPrbfamS2gXBjQewMAmA4ny/+u8Q5dSeuQIed98Q97y77CXj/pwa94c99atEe7uYBDI0jJK4w5jTdC6FoxBQfjfeCMBt9yMe/QVLzXhUqwOK7vTd/NH8h9sStnE0Zyj2OSfS0UKBL5uvoQUCEZsKGDwTgtmO+xrLljbSbC1rcDbJ4jWQa1pH9SermFd6R7WnQCPgEpnsGVmUrxjJ0WbN5H1+2qMGRkAuMXE4DNkCL5XX4Vnt24uM9KoPlKFX6mUjRs3MmjQIJ544gnefffdSj0mrSiNzamb2ZyymW2nt5FZkgmAn5sf7YLa0TqoNa0DWxPpE0mETwRhnmGYjeeOj7Y6rOSX5XOq4BRHc45yJOcI209v50TeCQCi/aO5qdVN3BBzA4EegZV+Pw6LhZLYWEr27adk314sJ04i7TZwSITRiHurVnh06IBHhw549+ld56c4PppzlG+OfsMPx38g35KPSZjoGtaVrmFdaRnQUtua8onAx+yDQfy5c9ealkbGl0vJW74ckk/jMBvJiAniUDMDm4NzSPa3k+UHLYJa0q9xP/o17kfPRj3xMnvhKCvDlpGB5fhxyo4dp+xIHMW7dlf01QsvL7YVF7PdYuHVn1YQUIen3K5rVOFXKu3+++/ns88+Y9u2bfS4wLq4doedvZl7WZ+0nvVJ6zmScwSAII8g+kb0pWd4T7qHdSfKP+qc4lIVKYUpbEjawPcnvmdvxl7MBjPXRV3H1E5TaeHfolrPXV9IKdmYvJHP9n3GrvRdmAwmhjcbzsgWI+kT0Qcft8rvn5AOB8WxsRSuXqMNnzx06M+pTQyCIh8zJcKKXUgMCPxLDLiX2c95DmNoiDYCp1s3PLt1Z+YP3/OPl19mxYoVXH/99TX63pVLU4VfqbTc3Fzat29PeHg427Ztw2QyUWgp5I+UP1iXuI4NyRvILcvFKIx0C+vGoMhBDGg8gFaBrapd6C/lSM4Rvj7yNcuOLsPisHB91PXc3+V+mvvVnVFINelMwf94z8fsy9xHY+/G3NHuDsbEjCHI48p2al+MPS+P0kOHsSYnYUlMxJaZid1qIbsog4ziDE6SRaIpnzxvkM0b06rrVfRvdy1dQrtgMpg4fPgwXbp0Ydy4cSxdurRGMimVpwq/ckW++eYb7njwDib8YwLG1kZi02KxOWz4u/szqMkghkQOoX+T/vi51f48MJklmSw4sIAv4r7Aarcysf1EpnWedkUt27ruRN4J3tr2Fn+k/EETnyZM7TSVG2Ju+EsXWm1ILEhkfdJ61iauPefvZEDjAaybt46jq45ycMdBGtWBdRXqG1X4lcsqshaxI20HfyT/wcbkjZwq0PpoIz0jGRE9gsGRg+ka1hWTwTUOfc8syWTWrlksO7qMYM9gZvScwaioUS49I2V1FVuL+XjPxyw6uAhPkycPdn2Q29rehtngGnPKnNkyXJ+0nl+P/UqpKEVIQZewLgxsMpD+jfvTPrg9RoPaoVsbVOFX/qLYWszezL3sSNvBttRt7M3Yi03a8DB60Cu8F518OvGPif+gXUQ71q1bh8HFJiM7Y3/mfl7b+hr7MvcxqMkgXur3Eo2861/rcvvp7bz4x4skFSZxU6ubeLTbowR7Xvni7bUhPj6ejp060nNMT275+y1sTN7IgawDAPi6+dKrUS96R/SmR6MetApopb4InEQV/gbOIR2cyj/F/qz97MvYx77MfRzKOoRN2jAIA22D2tIvoh99G/elW1g33I3atAQLFixg0qRJvP/++zz66KM6v4uLszvsLDm8hPd3vo/JYOKpXk9xY8sb60Xrv9hazMwdM/ki7gua+jbllf6v0DP8r+tnuwopJSNGjGDr1q3s37+/4kjw7NJstqVuY0vqFrakbiG5MBkAX7MvnUI70SmkE51DO9M+uD0hniF6voV6QxX+BsIhHdpOt/yTnMg9wYm8E8RlxxGXE0eJrQQAT5Mn7YPb0zW0Kz0a9aBrWFd83Xwv+HxSSkaPHs2aNWvYu3cvLVu2rM23c8US8xN5cdOLxKbFcnXTq3m5/8tXNPzT1ezL2MczG54hsSCRie0n8ki3R/A0ufZEZLNnz2batGl88sknTJs27aL3SylMYUfaDnama8d8HM89jkNqix2FeobSJqgNrQJbEe0fTbR/NC38W+iyT6kuU4W/HrA5bOSW5ZJVkkVWaRZZJVmcLjpNWnEaqUWpJBckk1SYdM7BUr5mX1oFtqJtUFvaBLWhQ3AHYgJirqifPikpiQ4dOtC1a1fWrFnjsl0+Zzikg4UHF/L+zvcJcA/gXwP+Rf8m+i3cUhV2h525++fy0e6PCPUK5bWBr9Er3PXnCkpISKBTp0706tWLVatWXdEWV7G1mANZBziUdYi4nDgOZx/mZN5JrI4/59X3d/cn0ieSxj6NCfcOJ9wrnDDvMEI8Qgj2DCbIIwhfN1+nji6rSxpk4c8uzabIUlRxWXLh93r29Wd+H7L8X/kF7ZKUFdefuZ9DOiouO6QDu7QjkdgddhzSgU3asDvs2Bw27FI7tzqs2nQFDm3KgjJ7GaW2UkrtpZTYSiixlVBkLaLIWkShpZB8Sz55ZXkUWgsvmN/f3Z9GXo2I9Ikk0jeSpr5NifKPIto/mhDPkBrp7jgzudY777zDjBkzqv18tSEuO46n1z/N8bzj3NX+Lh7r/hhuRtc/QCu1MJVnNjzDzvSdXBd1HS/0faFOtHQdDgfDhg0jNjaWvXv3ElUDS1HaHDaSCpI4kXeCU/mnSCxIPGf6j1J76V8eYxRG/Nz88Hf3x8fsg4+bD75uvniaPPE0eeJl8sLd5I67UTu5Gdy0o8iNZkwGE2ahnRsNRozCiMlgwiAMGIURIUTFuQEDQghE+WpZZ34WiD+vL//onbn+zP3OEFz4s3n29SFeIVXeymuQhf9fW/7FF3FfOCFRzRMIPEweeJm88DJ74WnyxMes/cH6uPkQ4B6Av5s//u7+hHhqrZtgj2DCvMLwMjt/1kIpJePGjePnn39mx44ddOzY0emvWRNKbaW8G/suS+OW0i6oHW8NfsulD/z6LeE3Xtr0EnaHnRf6vsCYmDF6R6q0mTNnMmPGDObOncu9997r9NeTUpJvySetOO2creC8sjztZNEaS4UW7XSmUVVsK/7LFCKu7OPhHzOwycAqPbZBFv69GXtJyE+o1H0v9C185pu74jrBX77RBQKDMGAQ2rf/mZ/PtBDOtBZMwoTJoJ3MBnPFPDVuRjc8jB6YDCaX3xGZnp5Op06diIiIYOvWrbjXoXnpV59azYubXsRit/B0r6e5qdVNLvX7LrYW89b2t/jm6Dd0CunEm4PepKlf3Zkgbv/+/fTo0YORI0eyfPlyl/rdXohDOiq2tiu2wO0WbA4bNmnTzsu30u0Ou7YlL2XFFv2ZrfuK3oDzegTOXA9UXA8X7l24nN7hvas8Sq1BFn6l5v3www/ccMMNPPPMM7z++ut6x7kiaUVpPP/H82xN3epSO373Zuzl2Q3PkliQyL0d7+Whbg+5zLj8yrBYLPTp04fk5GT2799PWFiY3pGUchcr/GoPiHJFxowZw5QpU3jzzTdZu3at3nGuSCPvRsweMZsnez7JhuQN3PT9Taw+tVq3PBa7hQ92fcBdP9+F1WFl3rXzmN5jep0q+gAvvPACu3fv5rPPPlNFv45QLX7lihUWFtKjRw+KiorYs2fPX5bPqwvisuN4buNzHMk5wnUtruOZPs/U2Pw2lbEvYx8vbnqRY7nHGBM9hmf6PFMnduCe79dff+Xaa6/l/vvv5+OPP9Y7jnIe1dWj1KidO3fSt29frr/+epYtW+byfboXYrVbmbt/Lp/u/RRfsy+PdX+MG1ve6NSjSPPK8vh4z8csObyEEM8QXur3EoMjBzvt9ZwpPT2dzp07ExISwvbt2/FUC527HNXVo9So7t278+abb/Ldd9/xySef6B2nSsxGM/d3uZ8vR39JC/8WvLz5ZW778Ta2n95e469lc9hYcngJo5aNYsnhJdzS6haWj11eZ4u+w+Hg7rvvJi8vj6VLl6qiX8fUeotfCNEU+B8QDjiA2VLK9y/1GNXid00Oh4PRo0ezevVqtm7dSpcuXZz7gnYr5CRAfjIUZUBxFlgKQTq0OeMNJvAKAq9g8AmHkFbgGVCpp5ZSsjJhJTNjZ5JalErv8N5M6jCJgU0GVmtrpthazPJjy/nfwf+RXJhMn/A+/L3X32kT1KbKz+kK3nnnHf7+97/z0Ucf8cADD+gdR7kIl+nqEUJEABFSyp1CCF9gB3CjlPLgxR6jCr/rSk9Pp1u3bnh5eREbG4v/Rda1vWKFGZCyE1J2Q+puyDisFX1pv9wjz+UbAWHtocVAiBoMEV3BePGjlkttpXwR9wX/O/g/0ovTaRnQkhtb3siwZsOI9I2s1Es6pIO9GXtZlbCK5ceXk1eWR5fQLkztNJXBkYPrZLfY2TZu3MjQoUO58cYb+eqrr+r8+6nPXKbw/yWAEN8BH0gpf7vYfVThd21nCsHYsWP5+uuvq1YICtPhxDpI2AgJmyDzSPkNQmu5h7WH4JbaKaApeIWAdwi4+YDBCMIAtjIoyYbibG2rIOMwZMRpXx4Zh7Sn8/CHtmOg400QNeSiXwJWu5Wf439m0cFFHMrWHtsmsA1dw7oSExBDjH8Mfu7azlgpJWnFaRzPPc7x3ONsSd1CRkkGJoOJwU0GM6njJLqFdbvy34kLctoXveIULln4hRAtgPVARyll/nm33QfcB9CsWbMeCQmVOxBL0ceZTf///Oc/TJ8+/fIPsJXBqc1wbBUcXwNp+7Xr3f2hWV9o3h8ie0FEZ3C/8ARyV6QwA+I3wNFf4fAKKMsH71DoMQl6Tga/iIs+NKkgid9P/c6axDXEZcdddOoMgDCvMLqEdmFYs2EMjhx80cnv6iK73c7IkSPZuHEjW7ZscX7XnlJtLlf4hRA+wDrg31LKby91X9Xid31SSm666SZ+/PFH1qxZw8CBFzjEPD9FK7xHfoUTa8FaBEY3rdBHD4XoqyCii9aCdyZrqZZj92I4slJ7vfY3wqAnoFGHSz5USkl6cTon8k5QbCuuuD7YI5iYgJh6VejP949//IN//etftTYlg1J9LlX4hRBm4EdgpZRy5uXurwp/3ZCbm0vv3r3Jz88nNjaWyMYRkBSrFdmjK+H0Pu2O/k2h1QhodQ20GATuOi6ZmH0Cts2Bnf8DS4H2BTD0GQhrp18mF/Ttt99y8803M3nyZObMmaN3HKWSXKbwC60DeAGQLaWcXpnHqMJfdxzZsZ6ZD41iXGd/rok2IUpzQBihaR9ofQ20ulYrqq62Q7A4GzZ/CFs/AUsRdJ0Aw/4BvuF6J9PdgQMH6NOnD506dWLt2rV1ao6mhs6VCv9AYAOwD204J8BzUsqfLvYYVfhdWEmu1ld/Yp3WfVO+EzW1wMFx0YIBk/6JiLkKPPWfE6dSirJg40zY+imY3LXun74PgdlD72S6yMnJoXfv3hQWFrJjxw4aN26sdyTlClys8Nf6qtlSyo1wkUmoFdcmJeQlQdJ2SNwKCX/A6f2ABJOn1lff5TZoOZyPP/qSV1/9FzObJPB4xzpS9AG8g+Haf0PPe+G3F+H3V2DXIhj1LsRcrXe6WmW1WrnttttISEhg7dq1qujXI7Ve+JU64kyRT9sPqXu0U/JOKDyt3W7yhKa9Yeiz2gicpr21FnK5l1/uwIEDB5kxYwYtW7ZkzJi6M688AMExcPtiOL4aVjwJC8dBx1tg5OvgU/8nIpNS8uijj/Lbb78xb948+vevWyuYKZem+zj+ylBdPU4iJZTmQu4pyD4J2cch67g2hj79sLazE9DG0reGxl21IZaRPaFRRzBeehbJ4uJiBg8ezOHDh/njjz/q7vA/ayls/I/WBWT2gmtfg653uN5+ihr03nvv8fjjj/P000/zxhtv6B1HqSKX6eOvClX4r4DdphXsklwoydFOxVlQlKlNc1CYBgWpUHBaa9GX5Z/7+DNTHYS1g9C22vDGRh2rPPImJSWFPn36ALB169a63V2QcQS+fwQSt2hDT8e8D4HN9U5V43744QfGjh3LuHHj+Oqrr1x+jWXl4hpm4U+KhcyjF7mxEu/7or8bedbtZ92n4v7yz9vOPj/nZ0f5yQ4Oe/n5mcs27WS3gd0CDqs2T42tDOxl2rmtVGuJWou1USjWYigr0M4vxmACn0baSBWfcPCPhIBm2imwOQTFOGVo5e7duxk0aBAxMTGsW7eubh/t6XBA7FxY9bL2/zj8Zeg1BepJcdy6dStXXXUV7du3Z926dXh7e+sdqXIc9j8/E3ZL+efkrM+O46zPlcP252fvnM/i2afy6ys+w+X+8hk/20XqRXVrbPQQ8Ktag6lhFv4VM2B7HRxzLIxakTaatZPBrPWfG83aAU8mDzB7audu3lr3g5uXdoSru582jYFngDaSxiNAm9rAO0T7Wafuid9++43rr7+eQYMG8fPPP9f9IYG5ifDDo9o+gGb9YewH2n6BOiwuLo4BAwYQEBDAH3/8QaNGVVvur0qk1LY+C05r03cUZWhbqSXZ5VuuudrtZQV/nluKtYaOtVgr5vXVhG+g1fAqPbRhFv6irL92ZZytUkXwIvepeKw473nEWbeLs84N514nDH+eDMbyYm886371z6JFi7jzzju57bbb+Pzzz+t+F4KU2tG/vzyntTCHPgv9Hr7kJHCuKjU1lf79+1NUVMTmzZuJiXHCl5ilWDtgLuuYtj8pJwHyErV9TPkpF99adffTGi0eftrP7r7alqnZq7zh46kNNjB7gNEdTG7audHtrMaT6axT+eet4vNXfo4oP2r8rM/o2Z9huMDP57noZ7can2mfRlrDrgpcZjhnrfIO1k6KS5g4cSKpqak89dRTBAcH88EHH9TtmR2FgG4TIWYY/PQkrHoJDnwLN8zSpp6oI7Kzsxk5ciQZGRmsXbu2+kVfSsg5WT4abC+kH9KO78hJ4JzuEK8QrZuxUQftwD7fcG02Vd9G5ZPwhWpbrXXwi9TVqd+oUquefPJJMjMzeeutt/Dx8eGNN96o28UftAneblsEB7/TvgBmD4U+D8BVz9bMBHNOlJ+fz8iRIzl8+DArVqygZ8+/NA4vrzQPEreVH9+xTRv2W5an3WYwQXAraNwNutyhDRwIjoGgaJf/3dRnqvArtUoIwRtvvEFhYWFF8f/HP/6hd6zqEwI63KjtiFv1T9jyIRxcrg39bD/WJbvviouLGT16NLt27eLbb79l+PBK9iOX5kH8Rji5/tyD+IQBwjpAx3FaoY/ook2nbarj+3PqIVX4lVonhGDWrFkUFRXx4osv4u7uzlNPPaV3rJrhGQhj3tPG+f/4OHx1NzQfCCNfc6nun+LiYsaOHcvGjRtZsmTJpQ+ws9sgORaO/a7tzE7ZqY14MXlC017apHbN+kGTHvpOuKdUmir8ii4MBgNz5syhrKyMp59+GpvNxnPPPad3rJrTtDfctw52LoA1/4ZPh2hfBkOegsAWukYrKirihhtuYM2aNcybN4/bbrvtr3cqTIejv2kzqx5fo3XdCINW3Ac9qW3ZRPZSrfk6ShV+RTcmk4mFCxdiNBp5/vnnsdlsvPjii3rHqjlGE/SaDB1vhvVvw7bPYO8X2syfg5/UdmzWssLCQkaPHs2GDRtYsGABd955p3aDwwGpu7S1Eo6uhJRd2vW+EdD+Bmg5XCv2dWWyPeWSVOFXdGUymViwYAFGo5GXXnqJkpISXnvttbq/w/dsngHaxG/9HoINM7WtgF2LoN0Y6PuANmV1LbzfnJwcxowZw+bNm1m0aBHjx14L+7/RWvbHVmlj5xFaS/7qF7SRNuGdXHL/hFI99Xscv1Jn2O12HnzwQWbPns3UqVP5+OOPMRqdvBKXXvKStGmfdy7QdpSGd4JOt2rrAPtXbkH3K5WamsrY60cQUHiM9x8bRzu30+Wteqm14lsOh5YjtHM1BLreaJgHcCl1ipSSF154gddee42bb76ZxYsX1/0jfC/FUgR7lmqt/5Sd2nVN+2jzAEUN1ibDq2ofupTawVKpu8k9uIYT65bQKciK2Si0g5cie0HMVdoxCE26O3+5S0UXqvArdcaZmSEHDx7Mt99+S3BwA2iBZh3XDv46vEI78Ek6tCNPQ1pDaBvt3DsUvIK0FrowaPdx2LTpDM5MwpdzUnuurOMVs6uW2WF3uqDZwNuI6HuL9uXi4afv+1VqhSr8Sp2yZMkSJk2aRPPmzVmxYgWtWrXSO1LtKcmFhE3aymYZcZBxGHITLv84YdDWMw5uCcExbE8s5ZHX5pLrFs7yH36ibdu2To+uuBZV+JU6548//uDGG2/E4XDwzTffMHToUL0j6cdm0SYsK87SvhjOHDAljNrOY68Q7dxgRErJa6+9xgsvvMCAAQNYvnw5ISEh+uZXdNEw5+pR6rQBAwawdetWRo0axfDhw3nnnXd47LHH6teIn8oyuZXPZXPpxd8LCgq45557+Oabb5gwYQJz586t3/tJlCqp49MjKvVddHQ0W7ZsYfTo0Tz++ONMmDCBoqIivWO5pMOHD9O7d2+WL1/OO++8w8KFC1XRVy5IFX7F5fn7+/Ptt9/y73//m6VLl9K7d2/27dundyyXsnDhQnr16kVWVha//fYbM2bMaJhbRkqlqMKv1AkGg4HnnnuOlStXkpWVRa9evfjggw+oC/uonCk/P5+JEydy11130a1bN3bs2MFVV12ldyzFxanCr9QpI0aMYO/evQwbNoxHHnmE0aNHk5SUpHcsXaxdu5auXbuydOlSXnnlFdasWUPTpk31jqXUAfV65+53u5PZciLLSc9+ZZvR5291i/NuE+XXaD9TsZkuBBiEdqvBIM69LAQGg8AoBEYDGA2GinOTQWAyCswGA2aTwGw0YDYacDcZcDMZcDcZ8TAb8HIz4Wk24ulmxNvNiMno+m2BsLAwfvzxRz744AOefvppOnTowLvvvsvkyZMbRPdGfn4+Tz/9NJ988gkxMTGsX7+e/v376x2rUqSUlNkcFFvslFjtlFi0U5nNTpnNQZnNjsUmsdgdWG0OrHYHNofEVn5ud0hsDonDIbFL7dwhwSG1cyklDimREu0ysmLJWyklkj+XwD37Nu3y2Tmr9O6q9ku5jEn9o2gTXrNrF9Trwh93uoDfD6XX+PNe6X/vX/+I5Dm3VSzffNYf5pmfKf+jPnO9vfwBjvI/cEcN/625mQz4upvw9TDh62HG39OMv5eZQC8zQV5uhPi6E+ztToiPG438PAj398DDXPtHfQoheOSRR7j++uuZPHkyU6dOZcmSJcyaNYv27dvXep7aIKVk2bJlTJ8+naSkJJ544gleffVVvLyqtixfTeTJL7GRVlBKWn4pmYVlZBVayCy0kFtsIafYQm6xlfxSG/klVgpKrRRb7Nhq+I/WIMBoEIizG0SCisuc1Zj6S8OK8xtlf16oShvCGc2OMV0aAzVb+NU4/jpOlhd/u0P7IrA5JHa7xObQWkhWuwObXWtBWWwOLHYHZVYHpTa7dm7VWl7FFjvFZTYKLTYKS20UlNooKLWSW2Ilr1g7zy22XPCLJsDLTGN/TyIDPYkM9KJFiBfNg72JCvamSaAnRoNzW+EOh4PZs2fz7LPPUlhYyKOPPspLL72En1/9OTr10KFDPProo6xatYqOHTsye/Zs+vXr5/TXLbXaOZlZREJWEfFZxSRkFZOUU0xybgkpuSWUWh1/eYyb0UCAl5lALzf8vcz4eZjx8zTh52HG292Il5sJbzft3MPNiIfJgKebETfjn1ujbmdtpWpbrwaMBoHJIDCeOZVv8SoXpw7gUqrN7pDkFFvIKrSQXlBKWn4ZafmlpOaVkJxTQnJuCUk5JRRb7BWPcTcZiA71oWWYD23DfWkX4Uu7CD/C/TxqvFsmIyOD5557jrlz5xIaGsrzzz/PtGnT6vSQxpSUFF599VXmzJmDj48Pr7zyCg888AAmU81urNvsDk5mFnEwNZ/DpwuIO13AsfRCEnOKz9liDfJ2IzLQkyYBnjQO8CTC34NGftopxEfbIvR1NzWILre6QBV+pVZIKckoLCMhq5iTGUUcyyjkWHohcacLSM4tqbhfiI8bnZr40ykygG5NA+jWLIAAL7cayRAbG8tTTz3FmjVraNasGS+99BITJ07Eza1mnr82pKen88477zBr1ixsNhtTp07ln//8J6GhodV+biklidkl7DyVw+7EXPYl53EgJa+i9W42CmLKv6xbhvkQE+pDVIg3zYO98PUwV/v1ldqjCr+iu/xSK3GnCziYks++5Dz2JeVxNL2govsoJtSbXi2C6B2lnSIDq953LaVk1apVPPfcc8TGxtKkSROmT5/Offfd59JdQEePHuXdd99lwYIFlJWVMXHiRF5++WWio6Or/Jx2h+RQaj5bT2az7WQWOxJyySwsA8DTbKRjEz86NQmgYxM/2kX4ERPqg5vJ9XfyK5enCr/ikorKbOxJymXXqVx2JOQQG59NfqkNgMhAT/rHBNM/JoT+McGE+Xlc8fNLKVm5ciVvv/02q1evxtfXlwkTJjB16lS6d+9e02+nSqxWKz/++COfffYZv/zyC25ubtx1113MmDGDNm3aXPHzSSk5klbIH8cy2XQ8i60nsygo/502DfKkV/MgujUPpHuzANo08q0TI7mUqlGFX6kTHA5JXFoBW09ksflEFltOZJNXYgWgTSNfBrYKYVCrEPpGB1/xaKIdO3bwf//3f3z55ZeUlpbSrVs3br/9dm6++WZiYmKc8XYuyuFwsGnTJr755huWLl3K6dOnady4MZMnT+bBBx8kPPzSc/KcL6OgjA1HM9hwNJONxzLJKNBa9M2DvegfE0yfqGB6RwXROMDTGW9HcVEuVfiFECOB9wEjMEdK+cal7q8Kf8N1ppti47FMNh7NZFt8NhabAzeTgT5RQQxpHcqQ1qG0DPOp9A7F3NxcFi9ezIIFC9i+fTsAXbt25dprr2X48OEMGDAAT8+aL5BpaWmsXr2a33//nZ9++onU1FTc3d257rrrmDx5MiNHjqz0TluLzcHOUzmsjctg/ZEMDqbmA9rO1wEtQxjUMoT+LYOr1V2m1H0uU/iFEEbgCDACSAK2A+OllAcv9hhV+JUzSq12tp7MZv2RDNYdyeBYeiEAjf09GNJG+xLo3zIEv0ruhIyPj+fbb79l2bJlbNmyBZvNhpubG126dKF79+50796d1q1bEx0dTZMmTSq1HGRpaSnx8fGcOHGCAwcOsHPnTnbu3MmRI0cAbe6h4cOHc/PNNzNq1KhK73NIzC5m/dEM1sVlsOl4FoVlNkwGQY/mgQwu/wJsH+GnhjgqFVyp8PcDXpZSXlt++VkAKeXrF3uMKvzKxSTnlrAuLoO1cekVxdBoEHRvFsDAlqEMah1C5yb+lerHLiwsZMOGDaxZs4bY2Fh27txJXl5exe1ms5mQkBACAwMJCAioGCYqpaSoqIicnBxycnLIyjr3aPFmzZrRo0cPevfuzbBhw+jevXulvkAKy2xsOZ7FxmOZrD+awYkMbVbSJgGef37JxQSrkTbKRblS4b8FGCmlnFJ++U6gj5Ty4fPudx9wH0CzZs16JCRUYgUipUGz2h3sOpXLuiPprD+Syf6UPKQEXw8TfaKC6BcTQr/oYNqE+1bqoDIpJfHx8Rw7doyTJ09y8uRJMjMzKwq81WqtuK+XlxeBgYEEBgYSHh5OdHQ0UVFRtGnTptKLoJRY7Ow6lcOm49r+jT2JudgcEg+zgT5RwVq3VptQokO81Th5pVJcqfD/Dbj2vMLfW0r5yMUeo1r8SlXkFFn443gmfxzLZPPxLOKzigHti6BH80B6tQiiS2QAnZr44+9Vu61mKSVJOSXsS85j16kctsfnsD85D5tDYjQIOjXxp19MMINahdCjeSDuJrUYunLlXGkFriTg7CkEI4EUHXIo9VygtxujOzdmdOfGAKTklrDlRBbb47Vho2vj4iru2yLYi7bhfrQO96VNI19ahHjRLKj6ByydOaDtVFYxJzKKiEvTjoo9lJpPVpEF0KY46NLUn6mDo+ndIoieLQJV943iVHq0+E1oO3eHAcloO3fvkFIeuNhjVItfcYa8Yiv7kvPYk5TL/uQ84k4XEJ9VdM58RIFeZhr5eRDq606Ijzt+Hia83E3nzGQqpTbKpthio7DMRm6xlYzCMjILy0jNLaXE+ucUFh5mA60baV8unZsG0CXSnzbhvqpFrziFy7T4pZQ2IcTDwEq04ZzzLlX0FcVZ/L3MDGwVwsBWf/bBl1rtHEsv5FR2ccUpPV8r4icyiigss1FUZrvgDJPuJgPe7ib8Pc2E+rrTLsKPq9qE0TzYi6ZBXrQI9qZZkJfTJ61TlMvRZVpmKeVPwE96vLaiXIqH2UjHJv50bOJ/yfuV2ew4zpqY0s1kUAVdqTPq9Xz8iuIsqmtGqcvUJB2KoigNjCr8iqIoDYwq/IqiKA2MKvyKoigNjCr8iqIoDYwq/IqiKA2MKvyKoigNjCr8iqIoDYwq/IqiKA2MKvyKoigNjCr8iqIoDYwq/IqiKA2MKvyKoigNTK0vxFIVQogMoC4uuhsCZOodohY1tPcL6j03FHX1PTeXUoaef2WdKPx1lRAi9kKr39RXDe39gnrPDUV9e8+qq0dRFKWBUYVfURSlgVGF37lm6x2gljW09wvqPTcU9eo9qz5+RVGUBka1+BVFURoYVfgVRVEaGFX4a4EQ4kkhhBRChOidxdmEEG8LIQ4LIfYKIZYJIQL0zuQsQoiRQog4IcQxIcQzeudxNiFEUyHEGiHEISHEASHEY3pnqg1CCKMQYpcQ4ke9s9QUVfidTAjRFBgBnNI7Sy35DegopewMHAGe1TmPUwghjMCHwHVAe2C8EKK9vqmczgbMkFK2A/oCDzWA9wzwGHBI7xA1SRV+5/sP8BTQIPaiSyl/lVLayi9uASL1zONEvYFjUsoTUkoLsBQYq3Mmp5JSpkopd5b/XIBWDJvom8q5hBCRwChgjt5ZapIq/E4khLgBSJZS7tE7i07uBX7WO4STNAESz7qcRD0vgmcTQrQAugFbdY7ibO+hNdwcOueoUSa9A9R1QohVQPgFbnoeeA64pnYTOd+l3rOU8rvy+zyP1jWwuDaz1SJxgesaxFadEMIH+AaYLqXM1zuPswghRgPpUsodQoihOsepUarwV5OUcviFrhdCdAKigD1CCNC6PHYKIXpLKU/XYsQad7H3fIYQ4m5gNDBM1t8DRZKApmddjgRSdMpSa4QQZrSiv1hK+a3eeZxsAHCDEOJ6wAPwE0IsklJO1DlXtakDuGqJECIe6CmlrIsz/FWaEGIkMBMYIqXM0DuPswghTGg7r4cBycB24A4p5QFdgzmR0FowC4BsKeV0nePUqvIW/5NSytE6R6kRqo9fqWkfAL7Ab0KI3UKIT/QO5AzlO7AfBlai7eT8sj4X/XIDgDuBq8v/b3eXt4aVOka1+BVFURoY1eJXFEVpYFThVxRFaWBU4VcURWlgVOFXFEVpYFThVxRFaWBU4VcURWlgVOFXFEVpYFThV5QqEEL0Kl9zwEMI4V0+P31HvXMpSmWoA7gUpYqEEP9Cm8PFE0iSUr6ucyRFqRRV+BWlioQQbmhz9JQC/aWUdp0jKUqlqK4eRam6IMAHbW4iD52zKEqlqRa/olSREOJ7tJW3ooAIKeXDOkdSlEpR8/ErShUIIe4CbFLKz8vX390khLhaSrla72yKcjmqxa8oitLAqD5+RVGUBkYVfkVRlAZGFX5FUZQGRhV+RVGUBkYVfkVRlAZGFX5FUZQGRhV+RVGUBub/AQrORQvD/TX+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "number = [0,1,2,3,4] #which wavefunctions to plot, starting counting from zero\n",
    "zoom = -3.0 # zoom factor for wavefunctions to make them more visible\n",
    "plt.plot(x,V(x),'-k',label=\"V(x)\") # plot the potential\n",
    "for num in number:\n",
    "    plt.plot(x,zoom*wavefunctions[:,num]+energies[num],label=num) #plot the num-th wavefunction\n",
    "plt.ylim(-1,10); # set limits of vertical axis for plot\n",
    "plt.legend();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"Energy\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}