psi_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical tests of ground state energy eigenfunction for the infinite square well potential\n",
    " This is an example of taking a python script and converting it to a jupyter notebook.\n",
    " \n",
    " Documentation for jupyter can be found at http://jupyter.org/index.html, including an [introduction to jupyter notebooks](https://jupyter-notebook.readthedocs.io/en/stable/examples/Notebook/What%20is%20the%20Jupyter%20Notebook.html).\n",
    " \n",
    " \n",
    "The original script (at the end of the notebook) was given to PHY-20006 students as an example on which to base their answers to a problem sheet question.\n",
    "\n",
    "This cell of the notebook is written in _markdown_ format. Double click on the cell to see the underlying rich-text markup \n",
    "language.\n",
    "\n",
    "Cells can be executed using the commands in the menu at the top of the notebook, or by using the keyboard short-cut *Shift-Enter*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following [\"magic\"](http://ipython.readthedocs.io/en/stable/interactive/magics.htm) handles the import commands for numpy and matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Physical constants in SI units from Math's handbook\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hbar = 1.055e-34  # Planck's constant / 2.pi\n",
    "a_0  = 5.292e-11  # Bohr radius\n",
    "m_e  = 9.109e-31  # Electron mass\n",
    "eV  =  1.602e-19  # Electron volts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysing $\\psi_1$, $n=1$ ground state for a particle trapped in a region of $a$.\n",
    "\n",
    "\n",
    "The 1-D time-independent Schrodinger equation (*TISE*) is  \n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{-\\hbar^2}{2m}\\frac{d^2\\psi}{dx^2} + V(x) \\psi = E \\psi\n",
    "\\end{equation}\n",
    "\n",
    "In this case we take $V(x)= 0$ in the region  $0<x<a$, in which case the ground state solution of the TISE is \n",
    "\n",
    "\\begin{equation}\n",
    " \\psi_1 = \\sqrt{\\frac{2}{a}}\\sin\\left(\\frac{\\pi x}{a}\\right)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters for the calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a =  4.32 a_0\n",
      "m =  2.33 m_e\n",
      "Grid size:  9999\n"
     ]
    }
   ],
   "source": [
    "# N.B. all numpy functions have were imported by the \"pylab\" magic statement at the top\n",
    "n_grid = 9999  # Do calculation on a grid of n_grid values.\n",
    "a = 4.32 * a_0  # Random value of well width a relative to Bohr radius\n",
    "m = 2.33 * m_e  # Random values of particle mass relative to electron mass\n",
    "n = 1\n",
    "print(\"a = \",a/a_0, \"a_0\")\n",
    "print(\"m = \",m/m_e, \"m_e\")\n",
    "print(\"Grid size: \",n_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a function for $\\psi_1$\n",
    " Note that all numpy functions are imported by the pylab magic above, so we can use \"sqrt\" instead of \"np.sqrt\", etc.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def psi(x, a, n=1):\n",
    "  norm = sqrt(2/a) # Normalisation factor\n",
    "  v = array(norm * sin(n*pi*x/a)) \n",
    "  return v\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid of $x$ values\n",
    "Exclude points at the boundaries where some functions may not be defined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid_step = a/(n_grid+1)\n",
    "x_grid = linspace(1,n_grid,n_grid) * grid_step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid of  $\\psi_1$ values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "psi_grid = psi(x_grid, a, n=n)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot of $\\psi_1$\n",
    "Not asked for in the original question, but very easy to do in a notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 102885.90731550295)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEOCAYAAABSLcpPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVXX+x/HXhx1RwAVREQVFRdyVzKzMMtcWq7GyZbK9\nxqa9SW2arKa9psV+ZctkWdliaZPllpq2mQtqqSAIiKK4iwKC7N/fHxxnkFQQ4X4v8Hk+HvfBueee\nc+77HpYP55zv+X7FGINSSinlCh62AyillGo4tOgopZRyGS06SimlXEaLjlJKKZfRoqOUUspltOgo\npZRyGS06SimlXEaLjlJKKZfRoqOUUspltOgopZRyGS/bAdxNixYtTEREhO0YSilVp6xZs2a/MSak\nsuW06FQQERFBXFyc7RhKKVWniMi2qixX66fXRGSaiOwVkY3l5jUTkUUikux8bVrutUkikiIiSSIy\nvNz8fiKywXltioiIM99XRD535q8UkYhy64xz3iNZRMbV9mdVSil1cq64pvMBMKLCvInAEmNMJ2CJ\n8xwRiQHGAt2cdd4UEU9nnanAbUAn53F0m7cAB40xUcArwPPOtpoBk4Ezgf7A5PLFTSmllOvVetEx\nxvwIZFaYPRqY7kxPBy4rN/8zY0yBMSYNSAH6i0hrINAYs8KUjcXwYYV1jm7rS2CIcxQ0HFhkjMk0\nxhwEFvHH4qeUUsqFbLVeCzXG7HKmdwOhznQYsL3ccjuceWHOdMX5x6xjjCkGsoDmJ9mWUkopS6w3\nmXaOXKyOJCcit4tInIjE7du3z2YUpZSq12wVnT3OKTOcr3ud+RlAeLnl2jrzMpzpivOPWUdEvIAg\n4MBJtvUHxph3jDGxxpjYkJBKW/wppZSqJltFZw5wtDXZOODrcvPHOi3SIilrMLDKORWXLSIDnOs1\nN1RY5+i2xgDfO0dPC4FhItLUaUAwzJmnlFLKklq/T0dEPgUGAy1EZAdlLcqeA2aKyC3ANuAqAGNM\nvIjMBBKAYuAuY0yJs6nxlLWE8wfmOw+A94CPRCSFsgYLY51tZYrIP4HVznJPGmMqNmhQSinlQlJ2\nUKCOio2NNXpzqFJKnRoRWWOMia1sOesNCZRSSjUcWnSUUkq5jBYdpZRSLqNFRymllMto0VFKKeUy\nWnSUUkq5jBYdpZRSLqODuCnlZnLyi9h5KJ+DeYUczi/mcEExuYXFCIIIeAj4eXsS6O9NsL83wY18\naB3kh5+3Z+UbV8oyLTpKWbIr6wgJO7NJ3J1Dwq5sUvceJuPQEXLyi6u1vdBAX8KbNiKyRQDd2gTS\nLSyIrq0Daeyrv+bKfehPo1Iusj0zj+Wp+1mZlsnKLZlkHDry39fCm/nTqWUT+kc2IyzYnzbB/jQP\n8KGxnxeNfb0IcApHqTGUGjhSWELWkSKyjxSRmVtIxqEjpGfmsT0zj+8T9/LFmrKRQEQgulUgZ3Vo\nzsCOzenfoRmBft5WPr9SoEVHqVpjjGFDRhaLEvawKGEPibtzAGge4EP/yGbcem4kPcKC6NKqCU1q\nsBAYY9iTXUD8ziw2ZGSxKi2Tj1duY9ovaXh6CAM6NGNEt1YM69aK0EC/GntfpapC+16rQPteU6dr\ne2Yes9dmMHvdDrYdyMND4IyIZgyNCWVwlxA6hjSmrLN018kvKmFd+iF+TN7HwvjdbNmXC8AZEU25\nsl84o3q21tNw6rRUte81LToVaNFR1VFcUsrC+D18tGIrK7ZkIgJndWjO5X3CuLBrKE0DfGxHPEbK\n3hwWbNzNV+sySN2XSyMfT0b1aM31A9rTOzzYdjxVB2nRqSYtOupUHMor5NNV2/no163szMqnXbNG\nXBXblsv7tiUs2N92vEoZY1ibfogv4rbzze87yS0sIbZ9U249N5KhMa3w9HDtEZmqu7ToVJMWHVUV\nBw4X8M5PW/hw+TaOFJUwsGNzbjo7kguiW9bZP9SHC4qZuXo77y9PY3vmEcKb+TN+cBRj+rXF21Nv\n6VMnp0WnmrToqJMpX2zyi0u4tFcb7jyvI11bB9qOVmNKSg2LEnYz9Yct/L79EOHN/Lnngk5c3icM\nLy0+6gS06FSTFh11PPlFJbz3cxpvLk0hr6is2Nx9QSeiWja2Ha3WGGNYmrSXlxdtZmNGNpEtApg0\nMpqhMaEubwih3J8WnWrSoqPKM8bw7fpdPDc/kYxDRxgaE8qEEV2IatnEdjSXMcawKGEPLyxMImXv\nYc6Oas4/Lo4hulX9ObpTp0+LTjVp0VFHbd6TwyOzNxC37SBdWwfyj4u6MjCqhe1Y1hSVlDJjxTZe\nWZxMTn4R157Zjr8NjybIX282VVp0qk2LjsovKuGNpSm89UMqjX29mDAimitjw+tsA4GadiivkFcX\nJ/Phr1tp3tiXJy/txojurfSUWwOnRaeatOg0bL+mHuCRrzaQtj+XK/qG8ehFMTRzs3ts3MWGHVlM\nnL2e+J3ZDI0J5Z+ju9MqSHs4aKi06FSTFp2GKb+ohOcXJPL+L1tp37wRT1/Wg3M6NdxTaVVVXFLK\nez+n8fKizfh4evDkZd24rHeYHvU0QFUtOtrvhWrw4ndmcd9nv5G89zA3Doxgwoho/H10mICq8PL0\n4I7zOjKieysenPk793/+O98n7uOp0d0JaqTXetQfadFRDVZpqeHtH7fw8qIkmjbyYfrN/Tmvc4jt\nWHVS++YBfH7HWbz1QyqvLNpM3NZM/nVVLwZ21KNFdSy900s1SAdzC7l5+mqeX5DI0JhQFt43SAvO\nafL0EO46P4pZfxmIv7cn1/17Ja8u3kxpqZ7CV/+jRUc1OOvSD3LRlJ9YnnKApy/vzhvX9nW7Djnr\nsl7hwXx7zzlc3ieMVxcnc+MHq8nMLbQdS7kJLTqqwTDGMH35Vq56+1c8PIRZfxnIdWe214vetaCR\njxf/urIXz1zegxWpB7h4yk/8tv2Q7VjKDWjRUQ1CQXEJD32xnslz4hnUKYS5d59Lj7ZBtmPVayLC\ntWe2Y9ZfBuLhIVz51nI+X51uO5ayTIuOqvf2Hy7g2ndXMmvtDu6/sDPv3hCrLatcqEfbIL69+xwG\ndGjOhFkbeHpuAiV6nafB0tZrql7btCubW6fHcSC3gDev68uoHq1tR2qQghv58P6NZ/DU3E28+1Ma\nqftyeW1s7xodplvVDXqko+qtxQl7+NPU5ZSUGr64Y6AWHMu8PD14/NJu/POy7vyweR9jpv7K9sw8\n27GUi2nRUfXSp6vSuf2jOKJaNubrv56t12/cyJ8HtGf6Tf3ZlXWEK6YuJ2Fntu1IyoW06Kh6xRjD\na4uTmTR7A4M6h/DZ7QMIDdT+wNzNOZ1aMOsvA/HyEK5++1d+TT1gO5JyEatFR0TuF5F4EdkoIp+K\niJ+INBORRSKS7HxtWm75SSKSIiJJIjK83Px+IrLBeW2KOG1gRcRXRD535q8UkQjXf0rlKiWlhkf/\ns5FXFm/mT33b8u4NsTTy0cuW7qpTaBNm/WUgrYL8GDdtFfM27LIdSbmAtaIjImHAPUCsMaY74AmM\nBSYCS4wxnYAlznNEJMZ5vRswAnhTRI52kDUVuA3o5DxGOPNvAQ4aY6KAV4DnXfDRlAX5RSWMn7GG\nGSvTGT+4Iy9d2RNvHVrZ7bUJ9ueLO8+iZ9sg7vpkLR/9utV2JFXLbP9WegH+IuIFNAJ2AqOB6c7r\n04HLnOnRwGfGmAJjTBqQAvQXkdZAoDFmhSnrMvvDCusc3daXwJCjR0Gq/jhSWMJtH8axMH4Pky+J\n4eER0XrDZx0S3MiHj289kyHRofzj63je+iHVdiRVi6wVHWNMBvASkA7sArKMMd8BocaYo8fZu4FQ\nZzoM2F5uEzuceWHOdMX5x6xjjCkGsoDmNf5hlDW5BcXc9MEqfk7Zz4tjenLT2ZG2I6lq8PP25K3r\n+3JprzY8Nz+R1xYno8Ou1E/WTng712pGA5HAIeALEbm+/DLGGCMitf6TJyK3A7cDtGvXrrbfTtWQ\n7Pwibnp/Nb9tP8SrV/dmdO+wyldSbsvL04NXru6Nj5cHryzeTEFxCX8b3kWPWusZm1dZLwTSjDH7\nAERkNjAQ2CMirY0xu5xTZ3ud5TOA8HLrt3XmZTjTFeeXX2eHcwovCPhDMxljzDvAO1A2iFvNfDxV\nmw7lFXLDtFUk7Mzm/67pw0i9B6de8PQQXvhT2fW4N5elUlBcyqMXddXCU4/YvKaTDgwQkUbOdZYh\nwCZgDjDOWWYc8LUzPQcY67RIi6SswcAq51RctogMcLZzQ4V1jm5rDPC90WP2Oi8rr4jr/r2SxF05\nvHV9Py049YyHh/DM5d25cWAE7/2cxhPfJOiptnrE2pGOMWaliHwJrAWKgXWUHW00BmaKyC3ANuAq\nZ/l4EZkJJDjL32WMKXE2Nx74APAH5jsPgPeAj0QkBcikrPWbqsNy8ou44f1VJO85zNs39OP8Li1t\nR1K1QESYfEkMHiJM+yUNP29PJozQU231geh/EMeKjY01cXFxtmOo48grLObGaatZm36QN6/ry7Bu\nrWxHUrXMmLJ7r2asTOe+Cztx34WdbUdSJyAia4wxsZUtp3fOqTohv6iE2z9cQ9y2TF4b20cLTgMh\nIvxzdHfyi0p5dXEyft6e3HleR9ux1GnQoqPcXmFxKeNnrOXnlP28dGUvLunVxnYk5UIeHsILY3pS\nUFzCc/MT8fXy0KbxdZgWHeXWSksN98/8je8T9/L05d0Z069t5SupesfTQ3jl6t4UFpfyxDcJBPp5\n8yf9WaiTbPdIoNQJGWN48tsE5q7fxaSR0Vx3ZnvbkZRF3p4evH5tH86Oas7Ds9azNHFv5Sspt6NF\nR7mtN5el8sHyrdx6TiR36Hl8Bfh6efLW9f2IbtWE8TPWsjb9oO1I6hRp0VFuaebq7by4MInLerfh\nkVFdbcdRbqSJnzcf3NSfloG+3PzBalL2HrYdSZ0CLTrK7SxO2MOkrzZwbqcWvDCmFx4eem+GOlZI\nE18+vLk/Xh7CuGmr2J2VbzuSqiItOsqtrNl2kLs+WUu3NoG8dX0/fLz0R1QdX/vmAXxwU3+yjhQx\nbtoqsvOLbEdSVaC/0cptpB/I4/YP42gV5Me0G88gwFcbV6qT6x4WxNt/7kfqvsPcNWMtxSWltiOp\nSmjRUW4h60gRN09fTXGpYdqNZ9Cisa/tSKqOODuqBU9f3p2fkvczeU689tPm5vRfSWVdUUkpf/1k\nLVv35/LRLWfSMaSx7Uiqjrn6jHak7c/jrR9SiWwRwK3ndrAdSZ2AFh1llTGGyXPi+Sl5Py+M6clZ\nHXWMPVU9Dw/vwrYDuTw9bxPtmwcwNCa08pWUy+npNWXVez+n8cnKdMYP7shVseGVr6DUCXh4CC9f\n1ZuebYO559N1bMzIsh1JHYcWHWXNooQ9PD1vE6N6tOKhYV1sx1H1gL+PJ+/e0I9mAT7cMn01e7O1\nKbW70aKjrEjek8N9n62jR1gQ/7qyt96Lo2pMyyZ+vHdjLDn5xdz58RoKiksqX0m5jBYd5XJZR4q4\n/aM1+Pt48faf++Hv42k7kqpnolsF8tKVvVibfojJX2uLNneiRUe5VEmp4d7P1rE9M4+p1/eldZC/\n7UiqnhrVozV/PT+Kz1Zv5+OV6bbjKIcWHeVSLy9KYlnSPh6/tBtnRDSzHUfVcw8M7cyQ6JY8MSee\nlVsO2I6j0KKjXGju+l28sTSVa/qHc92Z7WzHUQ2Ah4fwytjetGveiPEz1pJx6IjtSA2eFh3lEpt2\nZfPQF7/Tt10wj1/aDRFtOKBcI9DPm3dviKWwuJQ7Poojv0gbFtikRUfVuqwjRdzx0Rqa+Hnx1vX9\n8PXShgPKtTqGNObVsb3ZmJHN43Pibcdp0LToqFpljOFvX/zOzkNHmHp9X1oG+tmOpBqoIV1D/9uw\n4Iu47bbjNFhadFSt+vdPaXyXsIeJI6Pp114bDii77h/amYEdm/OPrzeyaVe27TgNkhYdVWvitmby\n3IJEhncL5ZZzIm3HUQpPD+G1sX0I9PNm/Iy1OgaPBVp0VK04cLiAv36yjrBgf14Y00sbDii3EdLE\nl/+7ti/pmXlM+HK93jjqYlp0VI0rKTXc9/lvZOYV8uZ1fQny97YdSalj9I9sxsQR0czfuJtpv2y1\nHadB0aKjatyUJcn8lLyfxy/pRvewINtxlDquW8+NZFhMKM/O28SabZm24zQYWnRUjfo5eT9Tvk/m\nij5hXNNfhypQ7ktEePHKXrQJ9ueeT38jK0+v77iCFh1VYw4cLuD+mb/RoUUAT13eXa/jKLcX5O/N\nlGv6sCc7n4mz9fqOK2jRUTXCGMNDX/xO1pEiXr+mL418dFBaVTf0Dg/moeFdmL9xN5+u0vt3apsW\nHVUj3v9lK0uT9vHIyGhi2gTajqPUKbn93A6c26kFT3wTz+Y9Obbj1GtadNRpi9+ZxXPzExkS3ZJx\nAyNsx1HqlHl4CP+6qhdN/Lz46ydrtX+2WqRFR52WvMJi7v50HcGNvHnxSr0fR9VdLZv48a+rerN5\nz2GemptgO069ZbXoiEiwiHwpIokisklEzhKRZiKySESSna9Nyy0/SURSRCRJRIaXm99PRDY4r00R\n5y+fiPiKyOfO/JUiEuH6T1m/PflNAmn7c3n16t40C/CxHUep03Je5xBuH9SBj1eks2DjLttx6iXb\nRzqvAQuMMdFAL2ATMBFYYozpBCxxniMiMcBYoBswAnhTRI52VzwVuA3o5DxGOPNvAQ4aY6KAV4Dn\nXfGhGoq563fx2ert/OW8jgyMamE7jlI14qFhXejZNoiHv1zPriwdf6emWSs6IhIEDALeAzDGFBpj\nDgGjgenOYtOBy5zp0cBnxpgCY0wakAL0F5HWQKAxZoUpa+/4YYV1jm7rS2CI6PmfGpFx6AgTZ6+n\nd3gw9w/tbDuOUjXGx8uDKWP7UFRS1iKztFSbUdckm0c6kcA+4H0RWSci/xaRACDUGHP0uHY3EOpM\nhwHl2zPucOaFOdMV5x+zjjGmGMgCmtfCZ2lQSksND80s+2WcMrYP3p62D5iVqlkRLQL4x8Ux/JJy\ngA+Wb7Udp16x+dfCC+gLTDXG9AFycU6lHeUcudT6vxkicruIxIlI3L59+2r77eq8ab+k8euWAzx2\nSQztmjeyHUepWnFN/3CGRLfkuQWJJGsz6hpjs+jsAHYYY1Y6z7+krAjtcU6Z4Xzd67yeAZTvV6Wt\nMy/Dma44/5h1RMQLCAIOVAxijHnHGBNrjIkNCQmpgY9Wf23ek8MLC5O4sGtLrorVbm5U/SUiPPen\nnjT29eK+z3+jsLjUdqR6wVrRMcbsBraLSBdn1hAgAZgDjHPmjQO+dqbnAGOdFmmRlDUYWOWcissW\nkQHO9ZobKqxzdFtjgO+N9nNRbYXFpdz/+W808fXi2St6avNoVe+FNPHl2St6EL8zm9eWbLYdp16w\n3VfJ3cAMEfEBtgA3UVYIZ4rILcA24CoAY0y8iMykrDAVA3cZY47ewTUe+ADwB+Y7DyhrpPCRiKQA\nmZS1flPVNGVJMvE7s3n7z/0IaeJrO45SLjG8Wyuuim3L1GWpnN+lJbEROgLu6RD9x/9YsbGxJi4u\nznYMt7Nm20GufGs5V/Rty0tX9rIdRymXOlxQzMjXfgRg/r2DaOxr+/919yMia4wxsZUtp82OVKXy\nCot5cOZvtA7yZ/IlMbbjKOVyjX29eOWq3mQcPMKT38TbjlOnadFRlXp67ia2ZeY5fVPpKKCqYYqN\naMYd53VkZtwOliburXwFdVxadNRJ/bh5HzNWpnPrOZEM6KC3OKmG7b4LO9E5tDETZ68n64gO+lYd\nWnTUCeXkFzFx1nqiWjbmwWFdKl9BqXrO18uTl67sxf7DhTz1rXYKWh1adNQJPTMvkd3Z+bw4pid+\n3p6Vr6BUA9CzbTB3DOrAF2v0NFt1aNFRx/VLyn4+XZXOred2oE+7ppWvoFQDcq+eZqs2LTrqD3IL\nipkwaz0dWgTwgHbmqdQf6Gm26tOio/7g+QWJZBw6wgt6Wk2pE9LTbNWjRUcdY8WWA3z46zZuGhip\nd14rVQk9zXbqtOio/8orLObhL9fTvnkj/jZcW6spVRk9zXbqql10RGRCTQZR9r24MIn0zDye/1NP\n/H30tJpSVVH+NNuPm3VolMpUueiIyMxyjy+AW2sxl3KxuK2ZfLB8K+POaq83gSp1iu4Z0okOIQE8\n8tUG8gqLbcdxa6dypJNtjLnKeVwJLK6tUMq18otKeHjWesKC/Xl4RLTtOErVOX7enjx3RU92HDzC\ny9/pEAgnU2nREZFPnVNpa0SkRbmX/l57sZQrvbkslS37cnnm8h4EaO+5SlVL/8hmXHdmO6b9ksbv\n2w/ZjuO2qnKkMxXIA/oBi0XkKxFpYozJrN1oyhU278lh6rIULu8TxqDOOmqqUqdjwshoQpr4MmHW\neopKdKTR46lK0YkAfgDuNMb0pmxY6SdqM5RyjdJSw6TZG2js68WjF3W1HUepOi/Qz5snR3cncXcO\n7/y4xXYct1SVotMJeApIFJF1wDDgChEZIiL6r3EdNmNVOmu2HeTRi2Jo3lhHAlWqJgzv1oqR3Vvx\n2pJk0vbn2o7jdqpSdF43xlxqjIkCLgA+BhoB1wHzajOcqj27s/J5fn4i50S14Iq+YbbjKFWvPHFp\nN3y9PJg4az2lpTo6c3lVKToLRCRNROYDk4EHgU+NMTcbY86o3Xiqtkyes5Hi0lKevrw7ImI7jlL1\nSstAP/4+qisr0zKZGbfddhy3UmlTJWNMXxHxBKKBLpQ1lV5Q28FU7VmwcTcL4/cwcWQ07ZsH2I6j\nVL109Rnh/Oe3DJ6et4kLolvSMtDPdiS3UKX7dIwxJcaYeGPMbGPMt8YYvfupjsrOL2LynI10bR3I\nLedE2o6jVL0lIjx7RU8Kikt5/Jt423Hchva91sC8sCCRfTkFPHdFD7w99duvVG2KbBHAPRdEMW/D\nbpYmaU/UoEWnQVmzLZOPV6Rz48BIeoUH246jVINw26AOdAwJYPLX8eQXldiOY50WnQaiqKSUR2Zv\npE2QHw8O04HZlHIVXy9P/nlZd9Iz83hjaYrtONZp0WkgPvhlK0l7cph8aTft6kYpFxvYsQVX9Anj\nrR9SSdl72HYcq7ToNAC7so7wyuLNDIluybCYUNtxlGqQHrmoK/7enjz6nw0Y03Dv3dGi0wD889sE\nSkoNj1/aTe/JUcqSFo19mTiyKyu2ZPLVugzbcazRolPPLUvay7wNu7n7gijCmzWyHUepBm3sGeH0\naRfM03M3kZXXMIe31qJTj+UXlfDY1/F0CAngtkEdbMdRqsHz8BCevqwHh44U8fzCRNtxrNCiU49N\nXZZKemYeT43ujq+XDj+tlDuIaRPITQMj+GRlWYe7DY0WnXoqbX8uU5elMrp3GwZGtah8BaWUy9w3\ntDOtg/x49D8bKW5g4+5o0amHjDE89vVGfL08+PsoHSdHKXfT2NeLyZd0Y9OubD5YvtV2HJfSolMP\nzd2wi5+S9/PgsM7ayaBSbmp4t1DO7xLCq4uT2ZudbzuOy2jRqWdy8ot48psEurUJ5M9nRdiOo5Q6\nARFh8iXdKCwu5dn5DadRgfWiIyKeIrJORL51njcTkUUikux8bVpu2UkikiIiSSIyvNz8fiKywXlt\nijg3o4iIr4h87sxfKSIRrv58rvbq4mT2HS7gqcu64+mh9+Qo5c4iWgRw+6AOfLUug1VpmbbjuIT1\nogPcC2wq93wisMQY0wlY4jxHRGKAsUA3YATwpjPOD8BU4DbKhtbu5LwOcAtw0Bn19BXg+dr9KHYl\n7c7hg+VbGXtGO/q0a1r5Ckop6+46P4qwYH8e+7phNCqwWnREpC1wEfDvcrNHA9Od6enAZeXmf2aM\nKTDGpAEpQH8RaQ0EGmNWmLK+JT6ssM7RbX0JDJF6eku+MYbH58TT2NeLh4d3sR1HKVVF/j6ePHpR\nVxJ35zBjZbrtOLXO9pHOq8DDQPnyHmqM2eVM7waOdhYWBpQf93WHMy/Mma44/5h1nIHnsoDmFUOI\nyO0iEicicfv27TutD2TLvA27+XXLAR4a3oWmAT624yilTsGI7q04J6oF//ouif2HC2zHqVXWio6I\nXAzsNcasOdEyzpFLrfeMZ4x5xxgTa4yJDQkJqe23q3F5hcU8PTeBmNaBXNu/ne04SqlTJCI8fmk3\n8gpLeGFB/W5UYPNI52zgUhHZCnwGXCAiHwN7nFNmOF+PDreXAYSXW7+tMy/Dma44/5h1RMQLCAIO\n1MaHsWnqslR2ZuXzxOhu2nhAqToqqmVjbjknkplxO1iXXn97KrBWdIwxk4wxbY0xEZQ1EPjeGHM9\nMAcY5yw2DvjamZ4DjHVapEVS1mBglXMqLltEBjjXa26osM7RbY1x3qNe9SmefiCPt3/cwmW923BG\nRDPbcZRSp+HuIZ0IDfTlsa/jKSmtV3+q/sv2NZ3jeQ4YKiLJwIXOc4wx8cBMIAFYANxljDk69ut4\nyhojpACpwHxn/ntAcxFJAR7AaQlXnzz5bQLeHsIk7XlAqTqvsa8Xj4zqyoaMLD5fvb3yFeogqWf/\n+J+22NhYExcXZztGlSxL2suN769m4sho7jyvo+04SqkaYIxh7DsrSNqTw9IHB9eZhkEissYYE1vZ\ncu54pKOqoLC4lCe/SaBDiwBuPjvSdhylVA0REZ4c3Z2c/GJe+i7Jdpwap0Wnjpr2Sxpb9ufy2CUx\n+Hjpt1Gp+qRLqyaMOyuCT1alszEjy3acGqV/reqgPdn5vL4kmQu7hjK4S0vbcZRSteDeCzvRrJEP\nT36TQH26DKJFpw56dt4mikoNj10cYzuKUqqWBPl78+CwLqzamsncDbsqX6GO0KJTx6zemsl/ftvJ\nHYM60K55I9txlFK16OozwunaOpBn5yWSX1RS+Qp1gBadOqSk1DD563jaBPkxfnCU7ThKqVrm6SFM\nviSGjENHeOfHLbbj1AgtOnXIF3HbSdiVzaRRXfH38ax8BaVUnTegQ3NG9WjF1GWp7Mo6YjvOadOi\nU0fk5Bfx0ndJxLZvysU9W9uOo5RyoUkju1JiDM/Xg8HetOjUEW8sTWX/4UIeuySGejo6g1LqBMKb\nNeKOQR0S5BWvAAARO0lEQVT4z287WbOtbg/2pkWnDth2IJdpP6cxpl9berYNth1HKWXBned1JDTQ\nlye+SaC0DvfLpkWnDnh2XiJensLfdHA2pRqsAF8vJo6MZv2OLGat3VH5Cm5Ki46b+zX1AAvidzN+\ncEdCA/1sx1FKWTS6Vxh92gXzwsIkDhcU245TLVp03FhJqeGf3yYQFuzPred2sB1HKWWZh4cw+ZJu\n7Msp4I2lKbbjVIsWHTf2vybS0fh5axNppRT0Dg/mir5hvPdTGtsO5NqOc8q06Lip8k2kL+qhTaSV\nUv8zYUQ0Xp7CM/M22Y5yyrTouCltIq2UOpHQQD/uOj+KhfF7+CVlv+04p0SLjhtKP5CnTaSVUid1\nyzmRtG3qz1NzN9Wpoa216LihZ+Zt0ibSSqmT8vP2ZMKIaDbtymbWmrrThFqLjpvRJtJKqaq6uGdr\n+rQL5sXvksitI02otei4EW0irZQ6FSLCoxfFsC+ngLfrSC/UWnTcyJdrtIm0UurU9HM6AX7nx7rR\nC7UWHTdxuKCYFxdu1ibSSqlTNmFENKUGXlyYZDtKpbTouIm3f0hl/+EC/nGxNpFWSp2a8GaNuOns\nCGavzWDDjizbcU5Ki44b2JV1hHd/2sLo3m3oFa5NpJVSp+6u86NoFuDDU3MTMMZ9m1Br0XEDLy3c\nTKmBh4ZpE2mlVPUE+nlz/9DOrEzL5LuEPbbjnJAWHcs2ZmQxe90Objo7gvBmjWzHUUrVYdecEU5U\ny8Y8Nz+RwuJS23GOS4uORcYYnpm3iWB/b8YPjrIdRylVx3l5evD3UV1J25/Lxyu22Y5zXFp0LFqa\ntJflqQe4d0gngvy9bcdRStUDg7uEcG6nFry2JJlDeYW24/yBFh1LiktKeWZeIpEtArhuQHvbcZRS\n9YSI8MiormTnF/H69+435o4WHUs+j9tOyt7DTBwZjbenfhuUUjWna+tAro4N58Nft7J1v3uNuaN/\n7Sw4XFDMK4s20z+iGcNiQm3HUUrVQw8M64y3pwfPzU+0HeUYWnQseGtZ2Vg5f7+oq94IqpSqFS2b\n+PGX8zqyIH43K7ccsB3nv7TouJjeCKqUcpVbz+1A6yA/npq7iVI3GXPHWtERkXARWSoiCSISLyL3\nOvObicgiEUl2vjYtt84kEUkRkSQRGV5ufj8R2eC8NkWcwwcR8RWRz535K0UkwtWfs6KXFm7GoDeC\nKqVqn7+PJw8N68KGjCy+Wb/TdhzA7pFOMfCgMSYGGADcJSIxwERgiTGmE7DEeY7z2ligGzACeFNE\njnbFPBW4DejkPEY4828BDhpjooBXgOdd8cFORG8EVUq52uV9wujaOpAXFyZRUFxiO469omOM2WWM\nWetM5wCbgDBgNDDdWWw6cJkzPRr4zBhTYIxJA1KA/iLSGgg0xqwwZR0OfVhhnaPb+hIYIpYuohhj\neHqu3giqlHItDw/hkVHR7Dh4hI9+tX/DqFtc03FOe/UBVgKhxphdzku7gaPNu8KA7eVW2+HMC3Om\nK84/Zh1jTDGQBTSv8Q9QBUuT9vLrFr0RVCnleud2CmFQ5xBe/z6FrLwiq1msFx0RaQzMAu4zxmSX\nf805cqn1q18icruIxIlI3L59+2p8+3ojqFLKtkkjo8nOL+KNZXZvGLVadETEm7KCM8MYM9uZvcc5\nZYbzda8zPwMIL7d6W2dehjNdcf4x64iIFxAE/KHtoDHmHWNMrDEmNiQkpCY+2jE+W603giql7Ora\nOpA/9W3LB79sZXtmnrUcNluvCfAesMkY83K5l+YA45zpccDX5eaPdVqkRVLWYGCVcyouW0QGONu8\nocI6R7c1BvjeuHigiZz8Il5drDeCKqXse2BoZ0Tg5UWbrWWw+W/32cCfgQtE5DfnMQp4DhgqIsnA\nhc5zjDHxwEwgAVgA3GWMOdoUYzzwb8oaF6QC85357wHNRSQFeACnJZwrvf3DFr0RVCnlFtoE+3Pz\nOZF8tS6DjRl2RhgVdx5hzobY2FgTFxdXI9vanZXP4JeWMiymFVOu6VMj21RKqdORnV/EeS8spWvr\nQGbcemaN/TMsImuMMbGVLacXGGrRK4s2U1Jq+NtwvRFUKeUeAv28uWdIJ5anHmDZ5ppvOFUZLTq1\nJHlPDl+s2c71A9rrjaBKKbdy3Zntad+8Ec/NS6TExd3jaNGpJc8vSCLAx4u7L+hkO4pSSh3Dx8uD\nh4dHk7Qnh1lrd1S+Qg3SolMLVm/NZPGmPdw5uCPNAnxsx1FKqT8Y1aMVvcOD+dd3SRwpdF33OFp0\napgxhmfmbSI00Jebz460HUcppY7r6Aije7ILmPZLmsveV4tODVsYv5t16Ye4/8LO+Pt4Vr6CUkpZ\n0j+yGUNjQpm6LJX9hwtc8p5adGpQUUkpLyxIIqplY8b0a1v5CkopZdmEEdEcKSrh9SXJLnk/LTo1\n6PPV29myP5cJI6Lx0u5ulFJ1QFTLxow9I5wZK9NJ259b6++nfxlrSG5BMa8uTia2fVMu7NrSdhyl\nlKqyey/shI+XBy8sSKz199KiU0Pe+zmN/YcLmDQqWru7UUrVKS2b+HHHoI7M37ibNdsO1up7adGp\nAfsPF/D2D6kM7xZKv/bNbMdRSqlTduu5kYQ08eWZeZuoze7RtOjUgNeXJJNfXMrDI6JtR1FKqWoJ\n8PXigaGdWbPtIAvj99Ta+2jROU1b9+cyY2U6V58RTseQxrbjKKVUtV3Zry1RLRvzwsJEiktKa+U9\ntOicphe/S8Lb04P7hmh3N0qpus3L04OHh3dhy75cZsbVTvc4WnROw+/bDzF3/S5uOzeSloF+tuMo\npdRpGxoTSr/2TXl18WbyCotrfPtadKrJGMOz8zfRLMCH2wZ1sB1HKaVqhIgwaWQ0e3MKeP+XrTW+\nfS061bRs8z5WbMnknguiaOLnbTuOUkrVmNiIsu5x3lqWSmZuYY1uW4tONZSUGp6fn0j75o249sz2\ntuMopVSNe3h4F3ILi/m/71NqdLtadKrhq3UZJO7O4aFhXfDx0l2olKp/OoU24arYcD5asZXtmXk1\ntl39i3mK8otKePm7JHq2DeKiHq1tx1FKqVpz34Wd8fQQ/vVdUo1tU4vOKZq+fCs7s/KZODIaDw/t\n7kYpVX+1CvLj5rMj+c9vO9mYkVUj29SicwoO5RXyxtIUBncJYWDHFrbjKKVUrbvjvI4EN/Lm+Rrq\nDFSLzil4c1kqOQXFTNDubpRSDUSQvzd/PT+Kn5L383Py/tPenhadKso4dIQPlm/l8j5hdG0daDuO\nUkq5zJ/Pak9YsD/PLdhEaenpdQaqRaeKXv5uMwAPDutiOYlSSrmWr5cnDw7rzMaMbL7dsOu0tqVF\npwo27cpm9rod3DgwgrBgf9txlFLK5S7rXXaW56WFSRQWV78zUC06VfD8gkSa+HoxfnBH21GUUsoK\nDw9hwogupGfm8cnKbdXfTg1mqpeWp+5nWdI+7jo/iuBGPrbjKKWUNed1DuGsDs2Z8n0KOflF1dqG\nFp2TKC01PDc/kTZBfowbGGE7jlJKWSUiTBwZTWZuIe/+uKVa29CicxJzN+xi/Y4sHhjWBT9vT9tx\nlFLKul7hwVzUszXv/pTG3pz8U15fi84JFBaX8tJ3SUS3asLlfcJsx1FKKbfxt2FdKCopZcqS5FNe\nV4vOCXy6Kp1tB/KYMCIaT+3uRiml/iuiRQDXntmOT1dtZ8u+w6e0rhad48jJL2LKkmQGdGjG4C4h\ntuMopZTbufuCTvh5efDSKXYG2iCKjoiMEJEkEUkRkYmVLf/uj1s4kFvIpJFdEdGjHKWUqiikiS+3\nDerAvA27WZd+sMrr1fuiIyKewBvASCAGuEZEYk60fHGJ4d2f0rioZ2t6hQe7KqZSStU5t57bgRaN\nfXh2ftU7A633RQfoD6QYY7YYYwqBz4DRJ1p4T04+RSWl/E27u1FKqZNq7OvFvUM6sSots8rrNISi\nEwZsL/d8hzPvuDJzC7nuzHZEtAio9WBKKVXXje3fjojmjaq8vFctZqkzROR24HbnacGTl/XY+KTN\nQCfWAjj9vsVrh2Y7de6aCzRbdblrNlfkal+VhRpC0ckAwss9b+vM+y9jzDvAOwAiEmeMiXVdvKrT\nbNXjrtncNRdotupy12zulKshnF5bDXQSkUgR8QHGAnMsZ1JKqQap3h/pGGOKReSvwELAE5hmjIm3\nHEsppRqkel90AIwx84B5VVz8ndrMcpo0W/W4azZ3zQWarbrcNZvb5BJjTm/oUaWUUqqqGsI1HaWU\nUm6iwRadyrrGkTJTnNfXi0hfN8o2WESyROQ35/GYi3JNE5G9IrLxBK/b3GeVZbO1z8JFZKmIJIhI\nvIjce5xlrOy3Kmaztd/8RGSViPzuZHviOMu4fL9VMZeVfVbu/T1FZJ2IfHuc16z9jv6XMabBPShr\nUJAKdAB8gN+BmArLjALmAwIMAFa6UbbBwLcW9tsgoC+w8QSvW9lnVcxma5+1Bvo6002AzW70s1aV\nbLb2mwCNnWlvYCUwwPZ+q2IuK/us3Ps/AHxyvAw2f0ePPhrqkU5VusYZDXxoyqwAgkWktZtks8IY\n8yNwsv4ubO2zqmSzwhizyxiz1pnOATbxxx4xrOy3KmazwtkXR/vM93YeFS9Au3y/VTGXNSLSFrgI\n+PcJFrH2O3pUQy06Veka55S6z6lBVX3fgc7h8XwR6eaCXFVha59VldV9JiIRQB/K/jsuz/p+O0k2\nsLTfnNNEvwF7gUXGGLfYb1XIBfZ+1l4FHgZKT/C69Z+1hlp06rq1QDtjTE/gdeA/lvPUBVb3mYg0\nBmYB9xljsl353pWpJJu1/WaMKTHG9KasF5H+ItLdVe99MlXIZWWficjFwF5jzBpXvF91NdSiU2nX\nOFVcpjZUpdue7KOH+KbsHiRvEWnhgmyVsbXPKmVzn4mIN2V/1GcYY2YfZxFr+62ybO7ws2aMOQQs\nBUZUeMnqz9uJclncZ2cDl4rIVspOy18gIh9XWMb672hDLTpV6RpnDnCD09pjAJBljNnlDtlEpJVI\n2ehyItKfsu/jARdkq4ytfVYpW/vMec/3gE3GmJdPsJiV/VaVbBb3W4iIBDvT/sBQoOKgLS7fb1XJ\nZWufGWMmGWPaGmMiKPu78b0x5voKi1n/HW0QPRJUZE7QNY6I3Om8/hZlPRiMAlKAPOAmN8o2BviL\niBQDR4CxxmmaUptE5FPKWua0EJEdwGTKLqRa3WdVzGZln1H23+efgQ3OdQCAR4B25bLZ2m9VyWZr\nv7UGpkvZIIwewExjzLdu8DtalVy29tlxucE+OzaPxX2hlFKqgWmop9eUUkpZoEVHKaWUy2jRUUop\n5TJadJRSSrmMFh2llFIuo0VHKaWUy2jRUUop5TJadJSqA0TkLRE523YOpU6X3hyqVB3g9BjQzxhT\nYjuLUqdDj3SUchNSNornUGf6KRF53ZnuCmw2xpSIyBgRWSFlI1f+LCIhVkMrdYq06CjlPiYDfxeR\n6ygb2+Y+Z/5IYIEzvdQYM8AY0wtYBFzl+phKVZ8WHaXchDP6qVA23PDYcqfShvO/onOjiKwSkd+B\n8UC+iASIyHQRedcpWEq5LS06SrkJEelBWS/Ghc7w0YhIIyDYGLNTRG6gbDjzC5wjnSQgHrgC+NIY\ncxtwqZ30SlWNFh2l3ICUjVM/g7Ix7A+LyNGBwc6nbKAwgB7AcmPMYRH5EzAQ2EDZQFxHhyDWhgbK\nrWnRUcoy52hmNvCgMWYT8E/Kru/AsddzPgDGi8gqyq75bDHG5FI2zn1bZxn9nVZuTZtMK+XGRGQt\ncKYxpugkywQA/wfkAz8bY2a4Kp9Sp0qLjlJKKZfRQ3GllFIuo0VHKaWUy2jRUUop5TJadJRSSrmM\nFh2llFIuo0VHKaWUy2jRUUop5TJadJRSSrmMFh2llFIu8/8b9aKKK6YyrgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a62bcf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(x_grid/a_0,psi_grid)\n",
    "xlabel('$x/a_0$')\n",
    "ylabel('$\\psi_1$')\n",
    "xlim(0,a/a_0)\n",
    "ylim(0,1.1*max(psi_grid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test that $\\psi_1$ is normalized\n",
    "$\\psi_1$ is real so $|\\psi_1^*\\psi_1| = \\psi_1^2$.\n",
    "\n",
    " The function used to do the numerical integration is *trapz*. You can see how to use this function by opening the cell below (double-click), placing the cursor inside the between \"trapz(\" and \")\" and using the keyboard short-cut *SHIFT-TAB*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normalisation test = 1.0000 (should be close to 1)\n"
     ]
    }
   ],
   "source": [
    "integrand = psi_grid**2\n",
    "test1 = trapz(integrand,x_grid)  # To see the how this function works, open t\n",
    "print('Normalisation test = {:0.4f} (should be close to 1)'.format(test1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate probability to observe in first quarter of the well\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P(x < a/4) = 0.0907 \n",
      "Expected value of P(x < a/4) = 0.0908 \n"
     ]
    }
   ],
   "source": [
    "first_quarter = (x_grid < 0.25*a) \n",
    "test2 = np.trapz(integrand[first_quarter],x_grid[first_quarter])\n",
    "print('P(x < a/4) = {:0.4f} '.format(test2))\n",
    "# Here we can check this against the calculated value\n",
    "truth2 = 0.25 -0.5/pi \n",
    "print('Expected value of P(x < a/4) = {:0.4f} '.format(truth2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate expectation value for $x$, $\\langle x \\rangle$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expectation value <x> = 2.1600 a_0 \n"
     ]
    }
   ],
   "source": [
    "integrand = psi_grid*x_grid*psi_grid\n",
    "test3 = np.trapz(integrand,x_grid)  # integral y(x)dx = trapz(y,x)\n",
    "print('Expectation value <x> = {:0.4f} a_0 '.format(test3/a_0))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot of $\\psi^2$\n",
    " ... so we can see if the answers above looks reasonable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 9623190840.122057)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEVCAYAAADzUNLBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclWX+//HXh10UUARxAUURBNx3syyz3CtbdLJdm6kx\nq29NTes07VnmtE2bma3TamVl5pJlZZq57wqIqOCKiCL7ev3+4NgPTfOgwHXO4fN8PHjIOecG3t0B\nb+77vu7rEmMMSiml1Kl42Q6glFLKPWhhKKWUcooWhlJKKadoYSillHKKFoZSSimnaGEopZRyiscV\nhoi8LSKZIrLRiW3PFZHVIlImIqOPe+0GEdnqeLuh9hIrpZR78LjCAN4Fhjm5bTowDvio6pMiEgo8\nAvQF+gCPiEiTmouolFLux+MKwxizCMiu+pyIxIjIPBFZJSK/iEi8Y9sdxpj1QMVxn2YosMAYk22M\nOQQswPkSUkopj+RjO0AdmQZMMMZsFZG+wGvAoD/ZvhWQUeXxLsdzSilVb3l8YYhII6A/8JmIHH3a\n314ipZRyTx5fGFSedjtsjOlWjY/ZDQys8jgS+KkGMymllNvxuGsYxzPGHAG2i8gYAKnU9RQfNh8Y\nIiJNHBe7hzieU0qpesvjCkNEPgaWAh1EZJeI/BW4BviriKwDNgGjHNv2FpFdwBjgDRHZBGCMyQae\nAFY43h53PKeUUvWW6PTmSimlnOFxRxhKKaVqhxaGUkopp3jUKKmwsDATHR1tO4ZnyNpa+W9YrN0c\nHmDHkR0ARAdHW82h1MmsWrUqyxgTfqrtPKowoqOjWblype0YnuGdkZX/jv/Wbg4PMH7eeADeGfaO\n5SRKnZiI7HRmOz0lpZRSyilaGEoppZyihaGUUsopWhhKKaWcooWhlFLKKVoYSimlnKKFoZRSyila\nGEoppZyihaGUUsopWhhKKaWcooWhlFLKKVoYSimlnKKFoZRSyilaGEoppZyihaGUUsopWhhKKaWc\n4lELKCnlCkrLK9h9qJD07AIOF5aSmVuMMYYZKzMI9PMmrJE/rUMDiQgOwNtLbMdVymlaGEqdoUP5\nJfyccoAVO7JZnX6YlP25lFeY319v0DoPgHtXrz/m4/x9vOjUKoTuUY3p164p58SGEeDrXafZlaoO\nLQylTkNecRmz1u5h1rrdrNhxiPIKQ5C/D91aN2ZQfDuimzYkKjSQsEZ+PLx8Bl4CT489n4KScjJz\ni0jPLmBbZj7rdh3m/d92Mn3xdgJ8vTg3NpzLe0RyQUIzfL31jLFyLVoYSlXDtgN5TP9lO7PW7ia/\npJz2zRpxy3kxDE6MoHOrELxOcIrJ36fyF39UaCAAHZoHHfN6cVk5y7dns2DzfuZv2sd3m/cTHuTP\nlb2iGH92NE0b+df+f5hSTtDCUMoJqZm5vLwwlW/W7cHX24uLu7bk6r6t6R7VGJEzuw7h7+PNgNhw\nBsSG8/BFifyccoCPlqXz6k+pvLV4O1f3bc3fz2tHs6CAGvqvUer0aGEo9ScO5Zfw3IJkPlqWToCv\nNzed246bBrQjrJb+6vfx9uKChAguSIggNTOX137cxru/7uCT5encNiiWG8+Jxt9Hr3MoO7QwlDqB\nigrDxyvSeXZeMnnFZVzXrw13XBhHaEO/OsvQvlkQz1/ZjdsGtWfSnCQmz0vikxXpPHlpJwbEhtdZ\nDqWO0sJQ6ji7Dxdy7+frWJJ6kH7tQnn0ko7ENw+2lqddeCOm39CLRSkHePSbTVz31nKu6tOaB0fE\nExTgay2Xqn+0MJSqYubqXTzy9SbKjWHSZZ25qk/UGV+jqCnnxoUz5/8G8MKCFKb9ksailAO8NLYb\nvaJDbUdT9YSO21MKKCot54GZ67lrxjoSWgQz745zubpva5cpi6MCfL15YEQCn0/oj7eXMHbab0z/\nJQ1jzKk/WKkzpIWh6r2M7ALGTF3Kx8szmDgwho9u6kvrpoG2Y/2pnm2a8M3t5zAovhlPfruFCR+s\nIq+4zHYs5eG0MFS9tmrnIUa9uoQdB/N58/pe3DssHh83uWEupIEvb1zXk4dGJvD9lkzGTF3KnsOF\ntmMpD+YePxlK1YK5G/Zy9Zu/ERTgw9e3ns3gxAjbkapNRPjbgHa8Pa43GdkFXPrqEjbsyrEdS3ko\nLQxVL03/JY2JH62mY8tgZt7Sn3bhjWxHOiPnxYXzxS398fX24i9vLOXnlAO2IykPpIWh6hVjDC8s\nSOHJb7cwrGNzPrqpn8dMvdGheRBf3tqf6LCG3PTeSuZt3Gc7kvIw1gpDRIaJSLKIpIrI/Sd4PURE\nvhGRdSKySUTG28ipPIcxhsnzknnph62M6RnJK1f38LjZYZsFBfDJTf1IbBnMrR+t5ss1u2xHUh7E\nSmGIiDfwKjAcSASuEpHE4za7FdhsjOkKDASeE5G6u81WeRRjDI/P3szUn7dxTd/WTL6ii8euRRES\n6MsHf+tLn+hQ7pqxjs9WZtiOpDyErSOMPkCqMSbNGFMCfAKMOm4bAwRJ5UD4RkA2oOMGVbUZY5g0\nZwvvLNnBjWe35clLO51wVllP0sjfh3fG9+ac9mHc98V6vlm3x3Yk5QFsFUYroOqfPbscz1X1CpAA\n7AE2AHcYYyqO/0QicrOIrBSRlQcO6IU+9UevLEzlzV+2c8NZbfj3RQkudzNebQnw9Wbadb3oFR3K\nPz5dy3eb9JqGOjOufNF7KLAWaAl0A14RkT9M6GOMmWaM6WWM6RUerhOyqWO9s2Q7zy1I4fLurXjk\n4o71piyOauDnzdvjetOxVQi3fbSGRTp6Sp0BW4WxG4iq8jjS8VxV44GZplIqsB2Ir6N8ygPMXL2L\nx77ZzJDECJ4d3cXjT0OdTCN/H94f34d24Q255YNVbNyt92mo02OrMFYAsSLS1nEheyww67ht0oEL\nAEQkAugApNVpSuW2fk3N4t7P19M/pin/vaq729y9XVtCAn1578Y+hDTwZfy7K9h1qMB2JOWGrPwU\nGWPKgNuA+cAWYIYxZpOITBCRCY7NngD6i8gG4AfgPmNMlo28yr2kZuYy4YNVtA1ryOvX9vS4obOn\nKyI4gHdv7ENRaTnj3llBTkGp7UjKzVj7s8sYM8cYE2eMiTHGPOV4bqoxZqrj/T3GmCHGmM7GmE7G\nmA9sZVXuIyuvmPHvrsDPp/LcfUgDXS+iqriIIN64ric7D+Zz8/9WUlL2h3EkSp1U/T5OVx6lqLSc\nv723kgO5xUy/oRdRoa4946wt/WPCmDK6K8u2Z/PoN5tsx1FuRBdQUh7BGMODX25gbcZhpl7bg25R\njW1HcmmXdm9F8v5cXv9pG4ktgrm2XxvbkZQb0CMM5RHeX7qTmat3c8cFsQzr1MJ2HLfwzyEdGNgh\nnEdnbWJZ2kHbcZQb0MJQbm/59myemL2ZC+KbcccFsbbjuA1vL+Glsd1pHRrIxA9Xs1vX0lCnoIWh\n3Nq+nCImfriaqNBAXhjbrd7ea3G6Qhr4Mu36XpSUVTDhf6soLiu3HUm5MC0M5bZKyiqY+OEqCkrK\neOO6ngQH6Iio09G+WSP+85eubNidw6Rvt9iOo1yYFoZyW//5LpnV6YeZMrorcRFBtuO4taEdm/PX\nc9ry3tKdfLt+r+04ykVpYSi39GNyJtMWpXFdvzaM7KIXuWvCfcPi6RbVmPu+WM+OrHzbcZQL0sJQ\nbmf/kSLunrGO+OZB/Gtkgu04HsPPx4tXr+mBj7cw8cPVFJXq9Qx1LC0M5VbKKwx3frKWwpJyj1wx\nz7ZWjRvw/F+6snnvEZ6eo9cz1LG0MJRbeWVhKkvTDvLEpZ1o36yR7TgeaVB8xO/XM35MzrQdR7kQ\nLQzlNlbtzOalHyrXthjdM9J2HI92z9AOxDcP4p7P1nMwr9h2HOUitDCUW8gvLuOuGeto2bgBj1/a\nyXYcjxfg682LY7txpLCU+77YgDHGdiTlArQwlFuYNGcL6dkFPDemK438dQq0uhDfPJh7h3Xg+y37\n+WRFxqk/QHk8LQzl8n5MzuTDZencNKAdfds1tR2nXrnx7Lac0z6Mx7/ZTNqBPNtxlGVaGMqlHcov\n4b7P19MhIoi7BsfZjlPveHkJ/xnTFT8fL+7+bB3lFXpqqj7TwlAu7d9fb+RQQQnPX9lVh9Ba0jwk\ngMcu6cia9MO8vXi77TjKIi0M5bK+WbeH2ev3cueFcXRsGWI7Tr02qltLLkyI4D/fJeupqXpMC0O5\npOz8Eh6ZtYmukSH8/dx2tuPUeyLCpMs64e/jxb2fr9dTU/WUFoZySY9/s4ncolKeHd0VH2/9NnUF\nzYIDePSSjqzceYh3f91hO46yQH8SlctZmLSfr9buYeLA9nRorrPQupLLurdiUHwzpsxP0gkK6yEt\nDOVScotKeejLjcQ2a8TE82Nsx1HHqTw11Rlf78pTUxV6aqpe0cJQLuXZecnsPVLE5NFd8PfRUVGu\nqHlIAP8emcjyHdl8ulJv6KtPtDCUy1i+PZv//baTcf2j6dG6ie046k+M6RVJ37ahPD1nC5m5Rbbj\nqDqihaFcQlFpOfd/sZ7IJg3455AOtuOoUxARJl3emaLSCp6crdOg1xdaGMolvP7TNtKy8pl0WWca\n6lxRbiEmvPI606x1e/g55YDtOKoOaGEo67Zn5fP6z9u4uGtLzo0Ltx1HVcMtA2NoF9aQh77aQGGJ\nrtDn6bQwlFXGGB7+eiP+3l78W5dbdTv+Pt48dVlnMrIL+e/CrbbjqFqmhaGs+nbDXn7ZmsXdQ+Jo\nFhxgO446DWfFNGVMz0jeXJRG0r4jtuOoWqSFoazJLSrl8W8207FlMNf2a2M7jjoDD45IILiBL//6\ncqPem+HBtDCUNc8vSOFAXjFPXdZZp/9wc00a+nH/8HhW7TzEzDW7bcdRtUR/SpUVm/bk8N6vO7i6\nT2u6RTW2HUfVgNE9IukW1Zhn5m7hSFGp7TiqFmhhqDpXUWF46KuNhDb0496h8bbjqBri5SU8MaoT\nB/NLeGFBiu04qhZoYag6N2NlBmvSD/PgiARCAn1tx1E1qHNkCFf3ac37S3fqBXAPpIWh6lROYSnP\nzk+md3QTLuveynYcVQvuGdqB4AAfHv56E8boBXBPooWh6tRL32/lUEEJj1zcERGxHUfVgsaBftwz\nNJ7l27OZtW6P7TiqBlkrDBEZJiLJIpIqIvefZJuBIrJWRDaJyM91nVHVrNTMXN5fuoOxvVvTqZUu\nuerJruwdRZfIEJ76dgu5egHcY1gpDBHxBl4FhgOJwFUiknjcNo2B14BLjDEdgTF1HlTVGGMMj32z\nmQZ+3vxzSJztOKqWeXsJj4/qRGZuMf/9Qe8A9xS2jjD6AKnGmDRjTAnwCTDquG2uBmYaY9IBjDGZ\ndZxR1aAFm/fzy9Ys/nFhHE0b+duOo+pAt6jGXNkrineW7GDbgTzbcVQNsFUYrYCqK6/scjxXVRzQ\nRER+EpFVInL9iT6RiNwsIitFZOWBAzpjpisqKi3nyW+3ENusEdedpXd01yf3DOtAgK83T32rU6B7\nAle+6O0D9ARGAkOBf4vIH85lGGOmGWN6GWN6hYfrTKeu6K3F20nPLuDhixPx1Tu665WwRv7cPqg9\nC5MydQp0D2Drp3c3EFXlcaTjuap2AfONMfnGmCxgEdC1jvKpGrIvp4hXf0xlSGIEA2K10OujcWdH\n06ZpIE/O3kxZeYXtOOoM2CqMFUCsiLQVET9gLDDruG2+Bs4RER8RCQT6Anpc62Ymz0uirMLw0MjE\nU2+sPJK/jzcPjkhga2YeHy9Ptx1HnQErhWGMKQNuA+ZTWQIzjDGbRGSCiExwbLMFmAesB5YD040x\nG23kVadn1c5DfLlmNzcNaEvrpoG24yiLhiRGcFa7pjy/IIWcAh1m666snVA2xswxxsQZY2KMMU85\nnptqjJlaZZspxphEY0wnY8yLtrKq6jPGMGnOFsKD/Jk4sL3tOMoyEeHfFyVyuLCUl3SYrdvSK5Cq\nVszftJ9VOw9x1+A4XaNbAZDYMpixvaN4f6kOs3VXWhiqxpWWV/DsvCTaN2vEmJ6RtuMoF3LX4Mph\ntpN0mK1b0sJQNe6TFRmkZeVz/7B4XRhJHSM8qHKY7Q9JmSzSYbZuR3+aVY3KKy7jpe9T6NM2lAsS\nmtmOo1zQuLOjaR0ayJPf6jBbd6OFoWrUm4vSyMor4cERCTobrTqhymG28aTsz2PGyl2246hq0MJQ\nNSbzSBFv/pLGyC4tdNlV9aeGdmxOrzZNeH5BCvnFZbbjKCdpYaga8+IPWyktr+DeoR1sR1EuTkR4\ncGQCWXnFTFuUZjuOcpIWhqoRqZm5fLoig2v6tqFN04a24yg30KN1E0Z2bsG0RWlkHimyHUc5QQtD\n1YjJ85IJ9PXm9kF6k55y3r3DOlBWUcHzC1JsR1FO0MJQZ2z59mwWbN7PhIExutaFqpY2TRtybb82\nzFiZQcr+XNtx1CloYagzcnQKkObBAdx4dlvbcZQb+r9BsTT09+HpOXozn6vTwlBnZO7GfazNOMxd\ng+No4OdtO45yQ00a+nHr+e35MfkAv6Zm2Y6j/oQWhjptJWWVU4DERTTiCp0CRJ2Bcf2jadW4AU/N\n2UJFhbEdR52EFoY6bR8vT2fHwQIeGJ6At5fepKdOX4CvN/cM7cCmPUf4et3xa6kpV6GFoU5LblHl\nNNVntWvKwA66kp46c5d0bUmnVsH8Z34KRaXltuOoE9DCUKfljZ/TyM4v4YER8ToFiKoRXl7CgyMS\n2H24kHeW7LAdR52AFoaqtn05RUxfnMYlXVvSJVKnAFE1p39MGIPim/Haj6lk55fYjqOOc8rCEJHB\nIvKmiHRzPL659mMpV/bCghTKKwz36BQgqhY8MDye/JIy/qsr87kcZ44wbgTuAa4VkUFAt5NtKCL3\n1VQw5ZqS9+Xy2aoMrj8rmqhQXadb1bzYiCCu7B3Fh8t2svNgvu04qgpnCiPXGHPYGPNPYAjQ++gL\nIjKjyttnwN9qK6hyDZPnJdHQ34fbztcpQFTtufPCOHy8vJgyP9l2FFWFM4Xx7dF3jDH3A+9Xee2I\nMeYvjrcxwPc1HVC5jl+3ZbEwKZNbz29Pk4Z+tuMoDxYRHMBNA9oye/1e1mYcth1HOZyyMIwxXx99\nX0Q+BgJFZLiIhAFPHrf5v2o4n3IRFRWGZ+Ym0TIkgHH9o23HUfXAzefF0LShH0/P2YIxejOfK6ju\nKKnXgQLgCiqPJl4QkaCjLxpjsmswm3IhszfsZf2uHO4e0oEAX50CRNW+Rv4+3HlhLMu2Z7MwKdN2\nHEX1CyMa+BmYYIzpBnwOPFbToZRrKS4rZ8r8JBJaBHNp91a246h6ZGyf1rQNa8gzc5N0/W8XUN3C\niKXyNFSSiKyh8iL45SJygYjo7b4e6sPf0snILuSB4fE6BYiqU77eXtw3rANbM/P4fJWu/21bdQvj\nZWPMJcaY9sAg4AMgELgGmFPT4ZR9OYWlvLxwKwNiwzg3Tv8mUHVvaMfm9GjdmOcXpFBQout/21Td\nwpgnIttFZC7wCHA38LEx5kZjTO9TfKxyQ1N/3sbhwlLuGxZvO4qqp0QqpwzJzC3mrV+2245Tr/lU\nZ2NjTA8R8QbigQ5UXvieVxvBlH3FZeW8vXg7l3VrRadWIbbjqHqsV3QoQztGMPXnbVzVtzVhurKj\nFdWeS8oYU26M2WSMmWmMmW2M0WNED7XrUCEGuGtInO0oSnHvsHiKyip0yhCLdPJBdUL5JWUcyCtm\nfP9oIpvoFCDKvpjwRlzVJ4qPlqWTdiDPdpx6SQtDnVB6dgHeXsLEgToFiHIdd1wQh5+PThliixaG\n+oPFW7PIKSwlsnEDQgJ9bcdR6nfhQf78/dwY5m7cx6qdh2zHqXe0MNQxKioMT8/dgr+PFxHBAbbj\nKPUHfxvQlvAgf50yxAItDHWMWev2sGnPEaKaBOKlK+kpF9TQ34d/XBjHyp2H+G7zfttx6hUtDPW7\notJypsxPplOrYJo20tlolev6S69IYsIbMnluEqU6ZUid0cJQv/vf0p3sPlzIg8MTEPToQrkuH28v\n7h+eQFpWPp+uyLAdp96wVhgiMkxEkkUkVUTu/5PteotImYiMrst89c3hghJeXriV8+LC6d8+zHYc\npU7pwoRm9IkO5cXvU8gr1tvB6oKVwnDcLf4qMBxIBK4SkcSTbDcZ+K5uE9Y/r/20jdziMu4frlOA\nKPcgIjwwIp6svBLeXJRmO069YOsIow+QaoxJM8aUAJ8Ao06w3e3AF4BOhl+LMrILeHfJDq7oEUlC\ni2DbcZRyWvfWTRjZuQVv/pJG5pEi23E8nq3CaAVUPfG4y/Hc70SkFXAZlYs2qVr0/IIUROCuwToF\niHI/9wztQElZBS/qlCG1zpUver8I3GeM+dMhECJys4isFJGVBw4cqKNonmPj7hy+XLObG89pS8vG\nDWzHUaraosMacm2/Nny6IoPUzFzbcTyarcLYDURVeRzpeK6qXsAnIrIDGA28JiKXHv+JjDHTjDG9\njDG9wsN1vYbqmjwviSaBvtwyMMZ2FKVO2+2D2tPA15vJ83TKkNpkqzBWALEi0lZE/ICxwKyqGxhj\n2hpjoo0x0VQuBTvRGPNV3Uf1XItSDvDL1ixuHxRLcIBOAaLcV9NG/twyMIYFm/ezfHu27Tgey0ph\nOKZEvw2YD2wBZhhjNonIBBGZYCNTfVNeYXh6bhKtQwO5tl8b23GUOmM3nt2WiGB/JumUIbWmWgso\n1SRjzByOW9bVGDP1JNuOq4tM9cmXa3azZe8RXr6qO34+rnwpSynnNPDz5u7BHbj3i/XM3biPEZ1b\n2I7kcfQ3RT1UWFLOf+Yn0zWqMRd10R8q5Tmu6BlJXEQjnp2XREmZThlS07Qw6qG3l2xn35Ei/jUi\nAdEJBpUH8fYSHhiewI6DBXy8PN12HI+jhVHPZOUV8/pP2xicGEGftqG24yhV4wZ2COesdk156Yet\n5BaV2o7jUbQw6pn//rCVwtJynQJEeayjU4Zk55fwxs86ZUhN0sKoR9IO5PHRsnSu6hNFTHgj23GU\nqjVdIhtzSdeWTF+cxr4cnTKkpmhh1COT5yXh7+PFHRfoFCDK890ztAPlFYYXFqTYjuIxtDDqiRU7\nspm/aT8TzoshPMjfdhylal1UaCDXnxXNZ6sySN6nU4bUBC2MesAYw6Q5W4gI9udvA9rZjqNUnbnt\n/PY09Pdh8rwk21E8ghZGPTBnwz7WpB/m7sEdaODnbTuOUnWmSUM/bj2/PQuTMvl1W5btOG5PC8PD\nlZRV8Oz8JOKbB3FFz0jbcZSqc+P6R9MyJIBn5iZRUaFThpwJLQwP98FvO9l5sID7h8fj7aU36an6\nJ8DXm7uHdGD9rhxmb9hrO45b08LwYDmFpfx34VbOaR/GeXE69buqvy7t3oqEFsFMmZ9EcVm57Thu\nSwvDg732Yyo5haU8MCJepwBR9VrllCHxZGQX8sFvOmXI6dLC8FDpBwt4Z8kOLu8eSceWIbbjKGXd\nuXHhDIgN4+WFW8kp1ClDTocWhod6eu4WvL2Ee4d1sB1FKZdx37B4cgpLef2nbbajuCUtDA/0W9pB\n5m7cx8SBMUQEB9iOo5TL6NQqhMu6teLtJdvZfbjQdhy3o4XhYcorDI9/s5mWIQHcdK7epKfU8e4a\nUjk1zvPf6ZQh1aWF4WE+X5XB5r1HuH9EAgG+epOeUseLbBLI+P7RzFyzi817jtiO41a0MDxIblEp\nU+an0LNNEy7WlfSUOqmJA9sTHOCr639XkxaGB3ntp21k5RXz8EWJOoxWqT8REujLnRfGsjg1iwWb\n99uO4za0MDxERnYBb/2ynct7tKJrVGPbcZRyedf2a0NcRCOe+HYzRaV6M58ztDA8xKQ5jmG0Q3Ul\nPaWc4evtxSMXdyQju5C3Fm+3HcctaGF4gKPDaG8ZGEPzEB1Gq5Szzm4fxtCOEbyyMJW9OTrM9lS0\nMNzcMcNoda0LpartoZGJlBvDM3N1zYxT0cJwcx8t28nmvUd4YESCrnWh1GmICg3k7+e24+u1e1i5\nI9t2HJemheHGDuYVM2V+Mv1jmnKRDqNV6rTdMjCGFiEBPPrNJsp1zYyT0sJwY1PmJ1NQUs5jl3TU\nYbRKnYFAPx8eGJHAxt1HmLEyw3Ycl6WF4abWZhzm05UZjOsfTWxEkO04Srm9i7u0oE90KFPmJ+ts\ntiehheGGKioMj3y9kbBG/txxYaztOEp5BBHh4YsTOVRQwgsLdJ6pE9HCcEMzVmawblcOD46IJyjA\n13YcpTxGp1YhXNO3Ne8v3cGmPTm247gcLQw3c7ighMnzkugd3YRLu7WyHUcpj3PPkHhCG/rx0Fcb\nqdAL4MfQwnAzz32XQk5hKY9d0kkvdCtVC0ICfXlwRAJr0g/zyQq9AF6VFoYb2bArhw+X7eS6fm1I\nbBlsO45SHuuy7q3o1y6UyfOSyMorth3HZWhhuImy8goe+HI9TRv5c/dQXXZVqdokIjx5aScKSsp4\neo7eAX6UFoabeG/pTjbuPsKjF3ckWC90K1Xr2jcL4qYB7fhi9S6WpR20HcclaGG4gd2HC3nuu2TO\n7xDOiM7NbcdRqt64fVAskU0a8NBXGykpq7AdxzprhSEiw0QkWURSReT+E7x+jYisF5ENIvKriHS1\nkdM2YyrvuTAGHh+lF7qVqksN/Lx57JKObM3M0ynQsVQYIuINvAoMBxKBq0Qk8bjNtgPnGWM6A08A\n0+o2pWuYv2k/32/J5B+DY4kKDbQdR6l654KECIZ2jOClH1LYeTDfdhyrbB1h9AFSjTFpxpgS4BNg\nVNUNjDG/GmMOOR7+BkTWcUbrcotKeXTWJhJaBHPj2W1tx1Gq3np8VCd8vb24/4sN9XoNcFuF0Qqo\nOsB5l+O5k/krMLdWE7mg575LYX9uEU9f3hkfb73cpJQtEcEBPDgigaVpB+v1vRku/1tIRM6nsjDu\nO8nrN4vIShFZeeDAgboNV4tW7sjmvaU7uOGsaLrpGt1KWTe2dxRntWvKpG+3sC+nyHYcK2wVxm4g\nqsrjSMdzxxCRLsB0YJQx5oTj2owx04wxvYwxvcLDw2slbF0rKi3nns/XE9mkAffoPRdKuQQR4enL\nO1NaUcGOpLINAAAOr0lEQVRDX9XPU1O2CmMFECsibUXEDxgLzKq6gYi0BmYC1xlj6tXUkc99l8z2\nrHwmX96Fhv4+tuMopRyiwxpy9+AOfL8lk9nr99qOU+esFIYxpgy4DZgPbAFmGGM2icgEEZng2Oxh\noCnwmoisFZGVNrLWtVU7DzF98Xau6dua/u3DbMdRSh1n/NnRdI0M4dFZm8jOL7Edp05Zu4ZhjJlj\njIkzxsQYY55yPDfVGDPV8f7fjDFNjDHdHG+9bGWtK5WnotbRMqQBD4xIsB1HKXUCPt5eTB7dhZzC\nUh7+eqPtOHXK5S961ycvfJ9C2oF8nrmiM430VJRSLiu+eTB3XBDL7PV7mbVuj+04dUYLw0WsTj/E\nm4vSGNs7igGxnnHxXilPdsvAGLpFNebfX22sN6OmtDBcQF5xGf/4dC0tQhrw4Eg9FaWUO/Dx9uL5\nv3SluKyce79YXy9GTWlhuIAnvtlMenYBL1zZTWeiVcqNtAtvxIMjEliUcoAPlqXbjlPrtDAsm7dx\nH5+uzOCW82Lo0zbUdhylVDVd168NA2LDmPTtFrZnefZcU1oYFu0/UsT9M9fTuVUId14YZzuOUuo0\niAhTRnfF11v4x6drKS333GnQtTAsqagw/POzdRSVlvPCld3w89H/FUq5q+YhATx9eRfWZhzmue88\n9z5j/S1lyVuLt/PL1iweGplI+2aNbMdRSp2hkV1acFWf1kz9eRs/p3jOvHZVaWFYsGrnISbPS2Jo\nxwiu6dvadhylVA155OJEOkQEcdena8k84nlDbbUw6tih/BJu/2g1LRoH8OzorrqCnlIeJMDXm1eu\n7k5+SRl3frqW8grPGmqrhVGHKioMd3+2jqy8El69ugchDXQIrVKeJjYiiMcu6civ2w7y+k+ptuPU\nKC2MOjTtlzQWJmXyr5EJdInUNS6U8lR/6RXFqG4teX5BCr9s9ZzrGVoYdWTptoNMmZ/MyM4tuP6s\nNrbjKKVqkYgw6bLOxDYL4vaP15CRXWA7Uo3QwqgDGdkF3PrRaqKbBvLMFZ31uoVS9UBDfx/euK4n\nFRWGv/9vFYUl5bYjnTEtjFpWUFLGzf9bRWl5BW9e34sgnfpDqXojOqwhL13VnS37jvDATPefb0oL\noxYZY7j38/Uk7TvCy1d1p1243m+hVH1zfodm3D04jq/W7uGtxdttxzkjWhi16LWftjF7/V7uHRrP\nwA7NbMdRSlkycWB7hnVszqQ5W/h+837bcU6bFkYt+XrtbqbMT+aSri2ZcF4723GUUhZ5eQnPX9mV\nTq1CuP3jNWzcnWM70mnRwqgFv6Ud5J7P1tOnbSjPju6iF7mVUgT6+TD9hl6ENvTjxndXsOdwoe1I\n1aaFUcO27s/l5vdXEhXagGnX9STA19t2JKWUi2gWFMA743tTWFLOje+uILeo1HakatHCqEH7cooY\n984K/Hy8eXd8HxoH+tmOpJRyMXERQbx+bU9SM/P463sr3Wq4rRZGDcnKK+aa6b+RU1jKO+N6ExUa\naDuSUspFnRMbxgtXdmPFjmxu+XAVJWXusYaGFkYNOFxQwrXTl7H7cCFvj+tN58gQ25GUUi7u4q4t\nmXRZZ35KPsA/3GSiQh/bAdxdblEpN7y9nLQD+bw1rpcus6qUctpVfVqTV1TGU3O20MDPm8lXdMHb\ny3UHyWhhnIGcglJueGc5m/Yc4fVrezIgNtx2JKWUm7np3Hbkl5Tx4vdbKS2v4LkxXfHxds2TP1oY\npykrr5jr3lrOtsw8Xr2mB4MTI2xHUkq5qTsvjMPX24sp85MpKavgpbHdXXLZZi2M07A3p5Brpi9j\nz+FC3hrXS48slFJn7Nbz2xPg680TszdT/MEqXr26Bw38XGtYvutVmIvbvOcIl7/2K5lHinn/xr5a\nFkqpGvPXc9ry1GWd+Ck5k7HTlnIgt9h2pGNoYVTDj0mZjJn6K8bAp3/vpxe4lVI17pq+bXjjul6k\n7M/jsteWkJqZazvS77QwnGCM4e3F2/nreyuIDmvIV7eeTceWOnRWKVU7BidG8Onf+1FUWsHlr/3K\nj0mZtiMBWhinlFtUym0freHx2Zu5ICGCGX8/i+YhAbZjKaU8XJfIxnw5sT+RTQIZ/+4Knv8u2fq9\nGloYf2LzniNc8soS5m3ax/3D43nj2p409NdxAkqpuhEVGsjMif0Z0zOS/y5M5Ya3l5N5pMhaHi2M\nEygtr+DlH7Yy6tXF5BWX8eHf+jLhvBi8XPiGGqWUZwrw9WbKmK5MvqIzK3ZkM/iFRXy9dreV1fv0\nz+XjbNqTw/1fbGDD7hwu6tKCx0d1IrShTiKolLLryt6t6RUdyj8/W8cdn6xl7oZ9PHxxIi0bN6iz\nDFoYDgfzivnPdyl8siKd0EA/XrumByM6t7AdSymlfhcT3ojPJ/Rn2qI0Xvw+hZ9TDjBxYAw3nduu\nTpZSqPeFcbighLeX7OCdJdsr56g/uy3/d0EsIQ18bUdTSqk/8PYSbhkYw0VdWjBpzhaeW5DCJysy\nuGVgDGN6ReLvU3vFUW8LY0dWPh8tT+fD33aSX1LO0I4R3DO0A+2bBdmOppRSpxQVGsjr1/bk19Qs\npnyXzENfbeSVhanceE40o3tG1cqpdGuFISLDgJcAb2C6MeaZ414Xx+sjgAJgnDFm9Zl8zYN5xfyQ\nlMnXa3ezJPUg3l7C8E7NuW1Qe+KbB5/Jp1ZKKSv6tw9jZkxTlqQe5OWFW5k0J4n/zE9hSMcILurS\nknPjwgj0q5lf9VYKQ0S8gVeBwcAuYIWIzDLGbK6y2XAg1vHWF3jd8a9TSssr2Hu4iC37jrA6/RAr\ndxxiTfohKgxENmnAP4fEMaZXFBHBek+FUsq9iQjnxIZxTmwYKftz+Xh5Ol+u2c3s9Xvx9/GiX7um\n9GzThO6tGxPbLIhmQf6nNerT1hFGHyDVGJMGICKfAKOAqoUxCnjfVI4d+01EGotIC2PM3pN90m0H\n8hj24iJyi8rYd6To95tcfL2FxJYh3DYoliGJEXRsGUzlAYxSSnmWuIggHrm4Iw+OSGDF9my+27yf\nJalZLNp6gKMjcf18vGgREkCgnw8NqzHBoa3CaAVkVHm8iz8ePZxom1bAMYUhIjcDNzseFs//x3kb\nT/QFU4FZZxC4BoQBWXYjnNTJs91otVjdc5+dxLu8WztJ/sij9lsdcdVcUEvZth77sI0zH+P2F72N\nMdOAaQAistIY08typBPSbNXnqrlAs50uV83mqrnAtbLZutN7NxBV5XGk47nqbqOUUqqO2CqMFUCs\niLQVET9gLH88YzQLuF4q9QNy/uz6hVJKqdpl5ZSUMaZMRG4D5lM5rPZtY8wmEZngeH0qMIfKIbWp\nVA6rHe/Ep55WS5FrgmarPlfNBZrtdLlqNlfNBS6UTWxMYKWUUsr96Gy1SimlnKKFoZRSyiluWRgi\nMkxEkkUkVUTuP8HrIiL/dby+XkR6uEiugSKSIyJrHW8P10Uux9d+W0QyReSE96lY3GenymVzn0WJ\nyI8isllENonIHSfYxtZ+cyZbne87EQkQkeUiss6R67ETbGNrnzmTzeb3m7eIrBGR2Sd4zco++wNj\njFu9UXmRfBvQDvAD1gGJx20zApgLCNAPWOYiuQYCsy3tt3OBHsDGk7xe5/vMyVw291kLoIfj/SAg\nxRW+16qRrc73nWM/NHK87wssA/q5yD5zJpvN77e7gI9O9PVt7bPj39zxCOP3aUWMMSXA0WlFqvp9\nWhFjzG9AYxGp7cUtnMlljTFmEZD9J5vY2GfO5LLGGLPXOCa8NMbkAluonG2gKlv7zZlsdc6xH/Ic\nD30db8ePrLG1z5zJZoWIRAIjgekn2cTKPjueOxbGyaYMqe42NnIB9HccUs4VkY61nKk6bOwzZ1nf\nZyISDXSn8q/Sqqzvtz/JBhb2nePUylogE1hgjHGZfeZENrDz/fYicC9QcZLXrX+fgXsWhjtbDbQ2\nxnQBXga+spzHHVjfZyLSCPgCuNMYc6Suv/6fOUU2K/vOGFNujOlG5ewMfUSkU118XWc4ka3O95mI\nXARkGmNW1fbXOlPuWBiuOq3IKb+mMebI0UNiY8wcwFdEwmo5l7NccioW2/tMRHyp/IX8oTFm5gk2\nsbbfTpXN9r4zxhwGfgSGHfeS9e+1k2WztM/OBi4RkR1UnsoeJCIfHLeN9X0G7lkYrjqtyClziUhz\nkcp51UWkD5X7/2At53KWS07FYnOfOb7uW8AWY8zzJ9nMyn5zJpuNfSci4SLS2PF+AyrXvEk6bjNb\n++yU2WzsM2PMA8aYSGNMNJW/NxYaY649bjOX+Pl0u9lqTe1NK1IXuUYDt4hIGVAIjDWOIRC1TUQ+\npnIESJiI7AIeofKin7V95mQua/uMyr/8rgM2OM57AzwItK6Sz8p+czKbjX3XAnhPKhdJ8wJmGGNm\n2/75rEY2m99vx3CRfXZsJkv7QimllJtxx1NSSimlLNDCUEop5RQtDKWUUk7RwlBKKeUULQyllFJO\n0cJQSinlFC0MpWqZiEwVkbNt51DqTOl9GErVMseNdT2NMeW2syh1JvQIQ6kaIJWLGQ12vP+kiLzs\neD8BSDHGlIvIaBH5TSoX8FksIuFWQytVTVoYStWMR4B/icg1VE41fqfj+eHAPMf7Pxpj+hljugIL\ngL/UfUylTp8WhlI1wLEQlFC5atrYKqefhvL/C2OcOJYIBSYCRSLSUETeE5E3HWWjlMvSwlCqBohI\nZyontytxrICHiAQCjY0xe0TkeipXZRzkOMJIBjYBlwOfG2NuAi6xk14p52hhKHWGpHKpzA+pXEYz\nT0SOrrFwPpVrLgB0Bn41xuSJyBVAf2ADlesaHF1JTS+KK5emhaHUGXAcRcwE7jbGbAGeoPJ6Bhx7\n/eJdYKKILKfyGkeaMSafyqU2Ix3b6M+jcmk6rFapWiIiq4G+xpjSP9mmIfAKUAQsNsZ8WFf5lKou\nLQyllFJO0UNgpZRSTtHCUEop5RQtDKWUUk7RwlBKKeUULQyllFJO0cJQSinlFC0MpZRSTtHCUEop\n5RQtDKWUUk75f3Ho7nQVeENDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10db84160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(x_grid/a_0,psi_grid**2)\n",
    "ytop = 1.1*max(psi_grid**2)\n",
    "plot(array([0.25*a/a_0,0.25*a/a_0]),array([0,ytop]))\n",
    "plot(array([test3/a_0,test3/a_0]),array([0,ytop]))\n",
    "\n",
    "xlabel('$x/a_0$')\n",
    "ylabel('$\\psi_{}^2$'.format(n))\n",
    "xlim(0,a/a_0)\n",
    "ylim(0,ytop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test time-independent Schrodinger equation and hence derive energy\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "del2psi = np.gradient(np.gradient(psi_grid,grid_step),grid_step) \n",
    "\n",
    "# Explicitly using potential V(x)=0 here so we know how to put in non-zero potentials later\n",
    "V = np.zeros(n_grid) \n",
    "\n",
    "# Left-hand side of TISE\n",
    "lhs = -hbar**2/(2*m)*del2psi + V*psi_grid\n",
    "\n",
    "# Divide LHS by u to get E\n",
    "# CAREFUL - must avoid points where psi is close to zero\n",
    "where_not_zero = (abs(psi_grid) >  1e-6 )  # List of true/false values\n",
    "E_grid = lhs[where_not_zero]/psi_grid[where_not_zero]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every point in E_grid should give the same value for E (+/- some numerical noise) but there may be a few points where the calculation is affected by numerical problems, so use the median value to get a robust estimate for the  typical value of E. Also print the mean absolute deviation to check that the energy derived is consistent for all valis values in the grid.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numerical calculation: E = 3.0909 +/- 0.0000 eV\n",
      "Expected value: E = 3.0909 eV\n"
     ]
    }
   ],
   "source": [
    "E = np.median(E_grid)              # Take median - more robust than average\n",
    "E_err = np.median(np.abs((E_grid-E))) # Should be close to 0 for valid solution\n",
    "\n",
    "#Report results in eV\n",
    "print(\"Numerical calculation: E = {:0.4f} +/- {:0.4f} eV\".format(E/eV,E_err/eV))\n",
    "\n",
    "# Expected value of energy\n",
    "E_expected  = n**2*pi**2*hbar**2/(2*m*a**2)\n",
    "print(\"Expected value: E = {:0.4f} eV\".format(E_expected/eV))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The original script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Test of ground state for infinite square well of width a\n",
      "a =  4.35 a_0\n",
      "m =  2.32 m_e\n",
      "Grid size:  9999\n",
      "Normalisation test = 1.0000 (should be close to  1.0)\n",
      "P(x < a/4) = 0.0907 \n",
      "Expected value of P(x < a/4) = 0.0908 \n",
      "Expectation value <x> = 2.1750 a_0 \n",
      "Energy E = 3.0615 +/- 0.0000 eV\n",
      "Expected value: E = 3.0615 eV\n",
      "--\n"
     ]
    }
   ],
   "source": [
    "# Numerical tests of solutions to T.I.S.E.\n",
    "# By Pierre Maxted, Nov. 2017\n",
    "from __future__ import (print_function) # So print() works in python 2.7\n",
    "from math import sqrt, pi\n",
    "import numpy as np\n",
    "\n",
    "# Physical constants in SI units from Math's handbook\n",
    "hbar = 1.055e-34  # Planck's constant / 2.pi\n",
    "a_0  = 5.292e-11  # Bohr radius\n",
    "m_e  = 9.109e-31  # Electron mass\n",
    "eV  =  1.602e-19  # Electron volts\n",
    "\n",
    "# Ground state of infinite square well potential, 0<x<a\n",
    "def psi_1(x,a):\n",
    "  norm = sqrt(2/a) # Normalisation factor\n",
    "  v = np.array(norm * np.sin(pi*x/a))  # x is an array, so use np.sin\n",
    "  return v\n",
    "\n",
    "# Start\n",
    "n_grid = 9999  # Do calculation on a grid of n_grid values.\n",
    "a = 4.35 * a_0  # Random value of well width \"a\" relative to Bohr radius\n",
    "m = 2.32 * m_e  # Random value of particle mass relative to electron mass\n",
    "\n",
    "print(\"\\nTest of ground state for infinite square well of width a\")      \n",
    "print(\"a = \",a/a_0, \"a_0\")\n",
    "print(\"m = \",m/m_e, \"m_e\")\n",
    "print(\"Grid size: \",n_grid)\n",
    "\n",
    "# Grid of x values\n",
    "# Exclude points at the boundaries where some functions may not be defined\n",
    "grid_step = a/(n_grid+1) \n",
    "x_grid = np.linspace(1,n_grid,n_grid) * grid_step\n",
    "\n",
    "# Test that psi_1 is normalised\n",
    "# psi_1 is real so |psi_1^* psi_1| = psi_1^2\n",
    "psi_1_grid = psi_1(x_grid, a)\n",
    "integrand = psi_1_grid**2\n",
    "test1 = np.trapz(integrand,x_grid)  # integral y(x)dx = trapz(y,x)\n",
    "print('Normalisation test = {:0.4f} (should be close to  1.0)'.format(test1))\n",
    "\n",
    "# Calculate probability to observe in first quarter of the well\n",
    "first_quarter = (x_grid < 0.25*a) \n",
    "test2 = np.trapz(integrand[first_quarter],x_grid[first_quarter])\n",
    "print('P(x < a/4) = {:0.4f} '.format(test2))\n",
    "# Here we can check this against the calculated value\n",
    "truth2 = 0.25 -0.5/pi \n",
    "print('Expected value of P(x < a/4) = {:0.4f} '.format(truth2))\n",
    "\n",
    "# Test expectation value for x\n",
    "integrand = psi_1_grid*x_grid*psi_1_grid\n",
    "test3 = np.trapz(integrand,x_grid)  # integral y(x)dx = trapz(y,x)\n",
    "print('Expectation value <x> = {:0.4f} a_0 '.format(test3/a_0))\n",
    "\n",
    "# Test time-independent Schrodinger equation and hence derive energy\n",
    "# TISE:  -hbar^2/2m d^2psi/dx^2 + V(x)*psi = E*psi\n",
    "# df/dx = gradient(f,grid_step) for regular grid of f(x) values\n",
    "del2psi = np.gradient(np.gradient(psi_1_grid,grid_step),grid_step) \n",
    "\n",
    "# Potential (V(x)=0)\n",
    "V = np.zeros(n_grid) \n",
    "\n",
    "# Left-hand side of TISE\n",
    "lhs = -hbar**2/(2*m)*del2psi + V*psi_1_grid\n",
    "\n",
    "# Divide LHS by psi to get E avoiding points where psi_1 is close to zero\n",
    "where_psi_not_zero = (abs(psi_1_grid) > 1e-6)  # List of true/false values\n",
    "E_grid = lhs[where_psi_not_zero]/psi_1_grid[where_psi_not_zero]\n",
    "E = np.median(E_grid)              # Take median - more robust than average\n",
    "E_err = np.median(np.abs((E_grid-E))) # median absolute deviation\n",
    "\n",
    "# Report results in eV\n",
    "print(\"Energy E = {:0.4f} +/- {:0.4f} eV\".format(E/eV,E_err/eV))\n",
    "E_expected  = pi**2*hbar**2/(2*m*a**2)\n",
    "print(\"Expected value: E = {:0.4f} eV\".format(E_expected/eV))\n",
    "print(\"--\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "&copy; Pierre Maxted, Keele University, 2017"
   ]
  }
 ],
 "metadata": {
  "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}