H2 Morse Potential.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# H2 Morse Potential"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the whole program at https://github.com/basnijholt/VMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit, brentq\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E = np.loadtxt('VMC H2/E.dat')\n",
    "E_all = np.loadtxt('VMC H2/E_all.dat')\n",
    "N, steps, eq_steps, N_s, minimization_steps = np.loadtxt('VMC H2/data_int.dat') \n",
    "s_min, s_max = np.loadtxt('VMC H2/data_float.dat') \n",
    "s_vector = np.loadtxt('VMC H2/s.dat') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/bnijholt/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:6: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n"
     ]
    }
   ],
   "source": [
    "#avg = np.mean(E, axis = 0)\n",
    "avg = np.min(E_all, axis=0)\n",
    "var = np.var(E, axis = 0)\n",
    "x = np.linspace(s_min,s_max,N_s)\n",
    "y_min = np.min(avg)\n",
    "y_max = np.max(avg[len(avg)/2:])\n",
    "s_min = s_vector[np.argmin(avg)]\n",
    "D = (y_min-y_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH01JREFUeJzt3Xt0lPW97/H3lwQQcBsRA7nZugEBL/XG8XK2LU5Re0Er\noq42pRfsObuuI9WWntKNgKsG3VY5a596wbNjtR5LT9FtS1GsSytYTN29eEFBhVB0WawQIF7KRYUC\nE77nj5mESZhJJnlm5snM83mtlZXJwy/zfEX95Jfv83t+j7k7IiJS+gaEXYCIiBSGAl9EJCIU+CIi\nEaHAFxGJCAW+iEhEKPBFRCKiPMg3m9lw4GHg48BbwBfdfVeacW8Bu4CDwAF3PzvIeUVEpPeCzvCv\nB5529/HAKmBuhnEHgZi7n6GwFxEJR9DAnwosTr5eDFyWYZzl4FwiIhJA0BAe6e6tAO6+HRiZYZwD\nK83sRTP7ZsBziohIH/TYwzezlcCo1EMkAvyGNMMz7dNwnrtvM7NKEsG/wd1/3+tqRUSkz3oMfHe/\nKNOfmVmrmY1y91YzqwLeyfAe25Kf3zWzR4CzgbSBb2ba3EdEpJfc3XoaE7Sl8xhwVfL1DGB51wFm\nNtTMjky+HgZ8BljX3Zu6e84+5s69C9hC4peP9o8tzJt3V07P4+7ceOONOX/PQn6oftWv+sOvoy8f\n2Qoa+AuBi8xsI3ABcBuAmVWb2ePJMaOA35vZGuA54NfuviLgebM2c+Y0qqoe6XSsqmoZM2deXqgS\nRET6hUDr8N39b8CFaY5vAy5Jvt4EnB7kPEHU1dUxefJ7tLQ08Pzz8IlPwAknQG1tbVgliYiEIlDg\nF4slSxoAOPtsuPvuxOd8iMVi+XnjAlH94VL94Sr2+rNhven/FIKZeb5q+sxnYPbsxGcRkVJhZngB\nLtoWlaOPhp07w65CRCQcCnwRkYhQ4IuIRIQCX0QkIhT4IiIRocAXEYkIBb6ISEQo8EVEIkKBLyIS\nEZEK/IoKBb6IRFekAl8zfBGJskgF/tChcOAA7N8fdiUiIoUXqcA3S8zyd+0KuxIRkcKLVOCD2joi\nEl0KfBGRiFDgi4hEhAJfRCQiFPgiIhGhwBcRiQgFvohIRCjwRUQiIlDgm9lwM1thZhvN7Ckzq8gw\nrsLMfmlmG8xsvZmdE+S8QSjwRSSqgs7wrweedvfxwCpgboZxdwJPuPuJwGnAhoDn7TMFvohEVdDA\nnwosTr5eDFzWdYCZHQV8yt0fAHD3uLvvDnjePlPgi0hUBQ38ke7eCuDu24GRacb8I/CemT1gZi+b\n2b1mNiTgeftMgS8iUVXe0wAzWwmMSj0EOHBDmuGe4RxnAt9y99VmdgeJVtCNmc7Z0NDQ8ToWixGL\nxXoqM2sKfBEpdk1NTTQ1NfX6+8w9XUZn+c1mG4CYu7eaWRXwTLJPnzpmFPAndx+d/PqTwBx3/0KG\n9/QgNfXko49g5MjEZxGRUmBmuLv1NC5oS+cx4Krk6xnA8q4Dki2fzWY2LnnoAqA54Hn7bOjQxH74\n2hNfRKIm6Az/GOAXwHHAX4EvuvtOM6sG7nP3S5LjTgN+AgwE/gJ8w93T7kqf7xk+wLHHwoYNUFmZ\n19OIiBREtjP8QIGfD4UI/LFj4ckn4YQT8noaEZGCKFRLpyjpwq2IRFFkA1+PORSRqIls4GuGLyJR\no8AXEYkIBb6ISEQo8EVEIkKBLyISEQp8EZGIUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhERDLw\nhw2DffvgwIGwKxERKZxIBr4ZVFRoAzURiZZIBj4kAl9tHRGJksgGvvr4IhI1CnwRkYhQ4IuIRIQC\nX0QkIiId+FqlIyJREunA1wxfRKJEgS8iEhGBAt/MhpvZCjPbaGZPmVlFmjHjzGyNmb2c/LzLzL4d\n5Ly5oMAXkagJOsO/Hnja3ccDq4C5XQe4++vufoa7nwlMBD4CHgl43sAU+CISNUEDfyqwOPl6MXBZ\nD+MvBN50980BzxuYAl9EoiZo4I9091YAd98OjOxh/JeAhwKeMycU+CISNeU9DTCzlcCo1EOAAzek\nGe7dvM9A4FISbaBuNTQ0dLyOxWLEYrGevqXXFPgiUqyamppoamrq9feZe8aM7vmbzTYAMXdvNbMq\n4Bl3PzHD2EuBme7+uR7e04PUlK0PPoCamsRnEZFiZma4u/U0LmhL5zHgquTrGcDybsZ+mX7SzgE4\n8kjYuxfi8bArEREpjKAz/GOAXwDHAX8FvujuO82sGrjP3S9Jjhua/PPR7t7tnLpQM3yAESPg9dcT\nn0VEilW2M/xAgZ8PhQz8MWNgxYrEZxGRYlWolk5R04VbEYmSSAe+nnolIlES6cDXDF9EokSBr8AX\nkYhQ4CvwRSQiFPgKfBGJCAW+Al9EIiLyga/HHIpIVEQ68P/2t/U89dQsXnmlOexSRETyLpKBH4/H\nmT37dm65ZSXvvLOQKVNW8P3v30FcG+uISAmL5NYKV145h+XLZxCPn9RxrLy8malTF7N06cK8nltE\nJNe0tUI3xo2rIx7v/PjdeLyC8ePrQqpIRCT/Ihn4M2dOo6qq82N1q6qWMXPm5SFVJCKSfz0+8aoU\n1dXVMXnye7S0NPDHP8LEiTB6NNTW1oZdmohI3kSyh59q4kS45x4466yCnVJEJKfUw89STQ1s2xZ2\nFSIi+Rf5wK+uhq1bw65CRCT/Ih/4NTUKfBGJBgW+WjoiEhGRD3y1dEQkKiIf+GrpiEhUKPDV0hGR\niIj8Ovy2NjjiCNi7F8ojeRuaiBS7gqzDN7PhZrbCzDaa2VNmVpFh3HfNbJ2ZvWpmS8xsUJDz5lJZ\nGVRWQmtr2JWIiORX0JbO9cDT7j4eWAXM7TrAzGqA64Az3f1UEts51Ac8b06pjy8iURA08KcCi5Ov\nFwOXZRhXBgwzs3JgKNCv4lUrdUQkCoIG/kh3bwVw9+3AyK4D3H0r8L+Bt4EWYKe7Px3wvDmlC7ci\nEgU9XqY0s5XAqNRDgAM3pBl+2NVWMzuaxG8CHwd2AUvNbLq7P5jpnA0NDR2vY7EYsVispzIDUUtH\nRIpJU1MTTU1Nvf6+QKt0zGwDEHP3VjOrAp5x9xO7jLkS+Ky7fzP59deAc9z92gzvWdBVOgD33gsv\nvAA/+UlBTysikhOF2i3zMeCq5OsZwPI0Y94GzjWzI8zMgAuADQHPm1Nq6YhIFAQN/IXARWa2kUSQ\n3wZgZtVm9jiAu78ALAXWAK+QaAndG/C8OaWWjohEQeRvvILE7P7007UWX0SKU7YtHQU+ibtthwyB\njz6CgQMLemoRkcD0xKteKCuDY4/VDF9ESpsCP0l9fBEpdQr8JK3UEZFSp8BP0vYKIlLqFPhJaumI\nSKlT4CeppSMipU6Bn6SWjoiUOgV+klo6IlLqFPhJaumISKnTnbZJ7c+23bNHd9uKSHHRnba9VFYG\nI0fqblsRKV0K/BTq44tIKVPgp9BKHREpZQr8FLpwKyKlTIGfQjN8ESllCvwU6uGLSClT4KdQS0dE\nSpkCP4VaOiJSyhT4KdTSEZFSpjttU+jZtiJSjHSnbR+UlUFlJWzfHnYlIiK5p8DvQm0dESlVgQLf\nzIab2Qoz22hmT5lZRYZx3zGz15If3w5yznzTSh0RKVVBZ/jXA0+7+3hgFTC36wAzOxn478B/AU4H\nLjGz0QHPmzdaqSMipSpo4E8FFidfLwYuSzPmROB5d9/n7m3As8DlAc+bN2br+dGPZvHKK81hlyIi\nklNBA3+ku7cCuPt2YGSaMeuATyXbP0OBKcBxAc+bc/F4nNmzb+ehh1by5psLmTJlBd///h3E4/Gw\nSxMRyYnyngaY2UpgVOohwIEb0gw/bD2lu//ZzBYCK4EPgTVAW3fnbGho6Hgdi8WIxWI9lRlYff18\nli+fQTx+EgBbt87ijjua2bRpPkuXLsz7+UVEstXU1ERTU1Ovv6/HwHf3izL9mZm1mtkod281syrg\nnQzv8QDwQPJ7bgE2d3fO1MAvlHHj6ojHO19zjscrGD++ruC1iIh0p+tEeMGCBVl9X9CWzmPAVcnX\nM4Dl6QaZWWXy88eAacCDAc+bczNnTqOq6pFOx6qqljFzZr+93CAi0itBA38hcJGZbQQuAG4DMLNq\nM3s8ZdyvzGwdiR8IM919d8Dz5lxdXR2TJ7/H+ec3MGRIAxMnNjB58vvU1taGXZqIZGHt2vVcfPGh\nBRddvy7Usf5MWyukccUV8MUvwpe+FGoZIpKFeDzO9dcv4qGHjK1br6G6+m5qa9fQ0jKRbdtmUlPT\nSH19G+7Oww+Xs3XrNXk5Nn063Hrrtaxbt5H58+/jhz+8mtNOO4m1a9d3+jofst1aQYGfxvz5MHgw\n/OAHoZYhIlm48so5nRZcwBzg68DJKaOuSX6cmrdjZWWvMnr0D/joo1jGHzztPxTKy3u8fNor2QZ+\nbs9aIiZMgCefDLsKEcnG2LFdF1zUAUd3GVUHjMjrsba2n/DGGzcDnwBg27Z32LZtLu0/ePrDyj/t\npZPGhAnw5z+HXYWIQPe98z/9qZk1a6YxcGDqgotpDBjw807vUVnZRmXl0rweGzp0F3BMypHDf/CE\nvfJPLZ00du9ObLHwwQcwQD8SRULRtTefrndeXt7IhAlw0kl/o7X10P+sb731Kscff6jd0r72oqWF\nvB0bPnw3zz03mu3br00e2cKAAUs4eHBOx5iqqkWsXn15zheDqIcfUG0t/OlP8LGPhV2JSDQd3puH\ndL3z8vJmpk5d3C9ukPzKVxo6/RBI94NnyZKGnJ9XPfyA2ts6CnyRwkld0TJmzOE3Q6brnYfdJkmV\njzDPJTUsMlAfXyQ7uVj/3r6X1cUXr+SJJxby6U+v4J573mfQoGWdzpWud64bJLOnlk4GixbBhg3w\n7/8ediUi4Ui3fjz12Mknj8vJ+vfq6kYGD36WzZtvpq3t0FLKsrJmqqvnMmbMGR3HMvXT+/vMOt/U\nww9o5Uq49VZYtSrsSkTyr7sgzxTcw4Y9y6ZNNxOPt4d0kPXvNwJXA6kXM1uYN28Zt9xyXc7/eUtN\ntoGPu/erj0RJ4Xv7bffq6rCrEAlmzZp1PmXKd3zt2vVpjx04cMC/970feU3N7Q5/95qa233s2Gle\nXr7OwVM+/ofDK12O/cBhS8rXd3X52h3+NctjL/iwYbd1OlZVdZdv2bIlxL+94pHMzR7zVRdtM6it\nTSzP3LULKtI+uFGkf8k8S1/I2rWN1Nc/mTJLX8iLLzZSXn4Dra03c/DgoZuDYAfZ3bh0CUOH/pw9\ne9qXHSbWv6cuQ6ysbAOW8u673+n2WFXVc5x77vvs2NHQcay2Fu1llWNq6XRj4kRobISzzw67EpHO\nemrBDBnyLG+91bknnn0r5UWGDVvFRx+lBvdNQEWXkF7EueduZseOoR3Hgqx/j3ofPggty8yB8eMT\nK3UU+BK29oC/6ab/xkMP/bZj5r5mTSNlZTewdWvfZ+nDhv28U7hnmm3DDlpaOh9bsuR/5eyfUfJP\ngd+NCRNg48awq5Co6a418/TTX2b//gUc2q8lU7gfHuS9aaUoyEuTAr8bEybAww+HXYWUsu7Cfc2a\nRI+9peXQ7H3//k/Teb8W0CxdsqXA74ZuvpJ86bxPTOKi6tChN7Bp06G+e/rZ+zRgKaBZuvSeLtp2\nY+9eOOaYxCZqOd6+WiIodTZ/882L0+wTk90F1COO+DKnn348gwcPBnTBU3TRNieGDEnsmrlpE5xw\nQtjVSDHprlXzhz80MmDA1jT7xGTbmhmvgJc+0Qy/B1OmwDXXwBe+EHYlUgwOf9xeI2Vlz9LScjPu\nh5ZImq1i4MDn2L9/XsexdMscNXuXbGiGnyPtfXwFvmTSXasm0yoa9/GMHXsflZUNHcfUd5d8U+D3\nYMIEeP75sKuQ/qK7Vs0f/9gIZNuqWcaKFf+mO0mloBT4PZgwARYvDrsKCVumVTV/+cuhJZM7d87C\nbBWDBy9m377UVo22DZD+QYHfgwkTEtsku4P1vBedlJj2Gf2+ffv43e+u62jVZLqb1X08Y8aoVSP9\nU6DAN7MrgQbgROAsd385w7jPAXeQeODK/e4e/rPIslRZmfj83nuHXkvp6zqjHzLkn9WqkaIXdIb/\nGok7QX6caYCZDQDuBi4AtgIvmtlydy+KW5rMoK5uPZdeeh/33HPoQRBSerq7+Lp3763AL4DvdoxX\nq0aKTaDAd/eNAGbdNjvOBt5w978mx/4HMBXo94HfPst7801jz56FTJnSyPTpK7j11msp151YJaPr\nbP6FF9JdfK0DtnHccf/M6NGJ56eqVSNFJ5tN83v6AJ4BzszwZ1cA96Z8/VXgrm7eK1fPBAjsiiv+\nxcvL13d6KEN5+Xq/4op/Cbs0CSj1ISDp/j2b/dYHDbpFD+SQokCuHoBiZiuBUamHAAfmu/uvc/nD\np11DQ0PH61gsRiwWy8dpejRuXN1hfdt4vILx4+tCqUeC6zqbf+ml9He9Zlonr3aN9AdNTU00NTX1\n/huz+anQ0wfdz/DPBX6T8vX1wJxu3itvPwV7a/PmzV5VtUizvCLX02weNJuX4kYIjzjM1Md/ERhr\nZh8HtgH1wJdzeN68qaurY/Lk92hpaWDNGjj+eDjlFM3yikW6tfPDhqW7MUqzeYmGQHvpmNllwCLg\nWGAnsNbdP29m1cB97n5JctzngDs5tCzztm7e04PUlC/z5sHAgbBgQdiVSLauvHJOmh0pE3vYHDjQ\neQ+b1asvV8BL0SrIXjru/ijwaJrj24BLUr7+DTA+yLnCNmkSLCyauweiK3Vp5fHHH34NBsZzwgma\nzUs0abfMLO3eDTU18P77kNyGXPqRrrtU/sM/NHLgwE7cR7Bv33Ud4zSbl1Kk3TJz7KijEtssrF4N\n550XdjUC3d8o9cEHsygra6a6ei5jxrzf8T2azUuUKfB7YdIkePZZBX7Yul6MXbOmkbKywy/GtrVV\n8PWvX8gtt1yX4Z1EokUtnV549FH48Y/hySfDriTaMl2MHTTo8AeKqH0jUaCWTh588pMwYwbE43rG\nbaGltm/Gjk1/MVZLK0W6pxl+L51ySmJ//IkTw64kGrpejB0+vBH3nXz00QgOHNDFWBHQDD9vJk2C\n//xPBX6h1NfP79S+2bFDF2NF+kqB30uTJsEvfgGzZoVdSelKbd+MGnV4+0YXY0X6Ri2dXmppgdNP\nh3fe0ROwcq1r++bIIxvZv38ngwaN4MMP1b4RyUQtnTyprU2syd+wAU7Ss1ByItNjBD/8cBbl5c0c\nffRcJk5U+0YkKAV+H7Svx1fgB5PNYwTjcbVvRHJlQNgFFKP2wJdg6uvnc+edn00+EHxwymMED6mq\nWsbMmZeHUp9IqVHg98GkSfD00+uZMmUWr7zSHHY5RWXt2vVcfHHi723MmK4XZA89RvD88xs4//wG\nJk9+X+0bkRzRRdteisfjzJmziDvuMA4evIaamkamT0fPue1B1wuyI0Y00ta2kz17RrB/vy7IigSR\n7UVbBX4vpbutv7y8malTF7N0qfZPziTd39uh9fRndBxLPBi8IYQKRYqXVunkiZ5zm73U9fQnnKD1\n9CJh0wy/l7Zs2cJZZz3K9u3XdhxTG6Kzru2byspGYCc7d2o7BJF80Aw/T1Kfc/v664mHofzTP2ld\neKqu2yG8++4sBgxopqZG2yGIhEmB3wftPebnn4evfhV+/vNw6+kPUts36dpeBw+qfSMSNrV0AnCH\nk0+Ge+9NbJ0cRV3bN6NGNTJkyE42bx5BW5vaNyKFoJZOAZjBVVfBT38a3cDv2r5pbZ2FWTO1tWrf\niPQ3CvyAvva1xBYLd94Jw4aFXU3hpWvfuKt9I9If6U7bgKqrExdtly0Lu5LCab9bdvXqZvbvn4bZ\nI53+XNshiPRPgXr4ZnYl0ACcCJzl7i9nGHc/cAnQ6u6n9vCeRdPDb/fLX0JjI6xaFXYl+dW1X19e\n3shxx8Epp/yN3bsPzR1085RIYRXkTlszGw8cBH4MzO4m8D8JfAj8rBQDf9++RMitXg3HHx92Nfmj\nu4xF+qdsAz9QS8fdN7r7G0C3J3L33wM7gpyrPxs8GOrr4bbbDm0MVipSNzsrK9NdxiLFTBdtcyAe\nj7NrV6LV0da2kLVrG5k+fUVRb6jWda/6Z55ppKxsJ8OHL2PHjtTllurXixSLHtPIzFYCo1IPAQ7M\nd/df56OohoaGjtexWIxYLJaP0+RM+9LEtrZEq2Pr1lnccUczmzbNL9pWR9fllnv3Hnr61Kmnarml\nSJiamppoamrq9ff1GPjuflFfCgoiNfCLQalsqNbTw8P19CmR/qHrRHjBggVZfV8ul2X2dMHAshhT\nlGbOnEZVVfEuTYzH48yefTsXX7ySJ55YyKRJK7j//vc58sjOa02L6Z9JRA4XqMFsZpcBi4BjgcfN\nbK27f97MqoH73P2S5LgHgRgwwszeBm509weCld5/pG6otncvvPwynHNO8bQ6urZvdu9OtG8qK/Xw\ncJFSor108mDBAnj1VfjVr8KupHvtLZzhwytYsuRqIDXMW5g3b5naNyJFQE+8CtHf/w6nnAJ33QVT\npoRdzeG63kAFt1JePoR4fE7HGG12JlI8tHlaiI44Au6+G771LRgxYj033ZS4EHraaSf1/M0FUF8/\nn0cfPbSqCBqIx7/Lscd+jZNPHgOofSNSihT4eXLhhXHKyhZx4YXGhx+Gvza/vX1z3XVX09xcR1tb\nRZcRs7n6arVwREqZWjp50l+2IWhv3zz4oLFt2zWYNXLGGTtpaRlBa6v2qxcpBWrphCzMtfmp6+kX\nLFjM8uUzOHgw8YPHfRavvtpMVdVczj9fK3BEokQz/DxJ97Dz4cMX8dprl/Puuzs7AjmXff2uF2OP\nPrqRPXteYv/+29AKHJHSpVU6/cBXvtJAS0vi9a5d0Nx8kM9+djgvvZQI5JqaRqZPJ1BfP3U2f/PN\niw9rIw0YsIqBA59j3755HcfUvhEpLQr8fujCC+fw29/OAIL39dM9S7as7CW2bj18Nn/SSbOprBzf\ncUT71YuUloJsjyy9c/bZdcDhff2KioGHbaucui1xumP19fO5887PsnXrLGAwra2z2Lr1GwwatLjT\n+1dVLWPFin+jqamh40NhLxJNmuEX0OF9/Thm3+Coo05j167rqKlppL6+DXfn4YfLO9o+XY8dfXQj\nZi+xY4dm8yKilk6/ldrXX7fuJd5//4fAJ1JGXJP8OLXbY2aJ3vz+/erNi0SdlmX2U6mz7XnzFnHr\nrcd0GVEHjOjxmPt4xo69j8rKQ++npZUi0h3N8EOUbulmZeVNQAXvvvudbo9pNi8i7dTSKRKpLR5I\nzNKBrI6pNy8ioMAXEYkMLcsUEZFOFPgiIhGhwBcRiQgFvohIRCjwRUQiQoEvIhIRCnwRkYgIFPhm\ndqWZrTOzNjM7M8OYOjNbZWbrzew1M/t2kHOKiEjfBJ3hvwZMA37XzZg48D/d/WTgvwLfMrMJAc/b\nLzU1NYVdQiCqP1yqP1zFXn82AgW+u2909zeAjHd4uft2d1+bfP0hsIHOe/qWjGL/D0b1h0v1h6vY\n689GQXv4ZnY8cDrwfCHPKyIiWWyPbGYrgVGphwAH5rv7r7M9kZkdCSwFvpOc6YuISAHlZPM0M3sG\n+J67v5zhz8uBx4En3f3OHt5LO6eJiPRSoR+A0t3J/i/Q3FPYQ3ZFi4hI7wVdlnmZmW0GzgUeN7Mn\nk8erzezx5OvzgK8Ak81sjZm9bGafC1q4iIj0Tr/bD19ERPKj39xpa2afM7M/m9nrZjYn7Hp6w8zu\nN7NWM3s17Fr6othvjjOzwWb2fPI3yNfM7Mawa+otMxuQ/O33sbBr6S0ze8vMXkn+/b8Qdj29ZWYV\nZvZLM9uQ/H/gnLBrypaZjUvpnKwxs13d/f/bL2b4ZjYAeB24ANgKvAjUu/ufQy0sS2b2SeBD4Gfu\nfmrY9fSWmVUBVe6+Nrma6iVgarH8/QOY2VB332NmZcAfgG+7e9GEj5l9F5gIHOXul4ZdT2+Y2V+A\nie6+I+xa+sLMfgr8zt0fSC4wGeruu0Muq9eSOboFOMfdN6cb019m+GcDb7j7X939APAfwNSQa8qa\nu/8eKMr/2KE0bo5z9z3Jl4NJLEYIfyaTJTOrA6YAPwm7lj4y+k+W9IqZHQV8yt0fAHD3eDGGfdKF\nwJuZwh76z7+kWiC1yC0UWeCUimK9OS7ZElkDbAdWuvuLYdfUC7cD36eIfkh14cBKM3vRzL4ZdjG9\n9I/Ae2b2QLItcq+ZDQm7qD76EvBQdwP6S+BLP1DMN8e5+0F3PwOoA84xs5PCrikbZnYx0Jr8Dcvo\nfnlzf3Weu59J4reUbyVbnMWiHDgT+D/Jf4Y9wPXhltR7ZjYQuBT4ZXfj+kvgtwAfS/m6LnlMCiTZ\nu1wK/D93Xx52PX2V/HX8GaBYlv6eB1ya7IM/BHzazH4Wck294u7bkp/fBR4h0aItFluAze6+Ovn1\nUhI/AIrN54GXkv8OMuovgf8iMNbMPm5mg4B6oNhWKxTr7Kxd1jfH9TdmdqyZVSRfDwEuAorigrO7\nz3P3j7n7aBL/3a9y96+HXVe2zGxo8jdDzGwY8BlgXbhVZc/dW4HNZjYueegCoDnEkvrqy/TQzoHc\n3mnbZ+7eZmbXAitI/BC63903hFxW1szsQSAGjDCzt4Eb2y8CFYOUm+NeS/bBHZjn7r8Jt7KsVQOL\nk6sUBgAPu/sTIdcUFaOAR5JbopQDS9x9Rcg19da3gSXJtshfgG+EXE+vmNlQEhdsr+5xbH9Yliki\nIvnXX1o6IiKSZwp8EZGIUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhEKPBFRCLi/wME4JqUZJTS\ndgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6e8dc6a898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fast non fancy plot\n",
    "plt.plot(s_vector,avg, '-p')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# datablock\n",
    "data_block_size = 100\n",
    "def data_block(x):\n",
    "    correct_length = len(x)-len(x)%data_block_size # make sure the list is an integer multiple of data_block_size\n",
    "    xx = x[:correct_length].reshape(-1,data_block_size)\n",
    "    avgs = np.mean(xx, axis = 0)\n",
    "    var = np.std(xx,axis = 0)\n",
    "    avg = np.mean(avgs, axis = 0)\n",
    "    error = np.std(avgs)/np.sqrt(np.size(avgs))\n",
    "    return avgs, var, avg, error\n",
    "\n",
    "avg = np.zeros(int(N_s))\n",
    "error = np.zeros(int(N_s))\n",
    "\n",
    "for i in range(int(N_s)):\n",
    "    avg[i] = data_block(E[:,i])[2]\n",
    "    error[i] = data_block(E[:,i])[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting the data to the Morse potential and find the spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D_e: 0.349708450087 a: 0.735882020509 r_e 1.61182703519 minimum energy: -1.16979204679\n",
      "[ 0.0172305   0.03819248  0.04071281  0.01051653]\n"
     ]
    }
   ],
   "source": [
    "r_e = s_min\n",
    "def V(r, D_e, a, r_e, const):\n",
    "    return D_e*(1.0-np.exp(-a*(r-r_e)))**2 + const\n",
    "\n",
    "popt, pcov = curve_fit(V, s_vector, avg, maxfev=10000)\n",
    "\n",
    "[D_e, a, r_e, const] = popt\n",
    "print('D_e:', D_e, 'a:', a, 'r_e', r_e, 'minimum energy:', const)\n",
    "print(np.sqrt(pcov.diagonal()))\n",
    "\n",
    "# the spectrum\n",
    "mu = 916.57\n",
    "w_0 = a*np.sqrt(2*D_e/mu)\n",
    "delta = D_e/w_0\n",
    "v_max = int(2*D_e/w_0)\n",
    "def e_v(v):\n",
    "    return w_0*(v+0.5)-w_0**2/(4*D_e)*(v+0.5)**2+const\n",
    "\n",
    "# find where to draw the lines for the spectrum\n",
    "def roots(x):\n",
    "    return V(x, D_e, a, r_e, const)-e_v(i)\n",
    "def find_left():\n",
    "    return brentq(roots, 0.0, r_e)\n",
    "def find_right():\n",
    "    return brentq(roots, r_e,1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The constants in units of cm$^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40005.8008431 510.717525502\n"
     ]
    }
   ],
   "source": [
    "from scipy import constants\n",
    "#36117*6.63e-36*3e10/1.6e-19*27.2-1\n",
    "Har = constants.physical_constants['Hartree energy in eV'][0]\n",
    "print((D_e+1)/(constants.h*constants.c/constants.e*Har), \n",
    "      np.sqrt(pcov.diagonal()[0])/(constants.h*constants.c/constants.e*Har))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Morse potential with spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEnCAYAAACKbmVGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnWl4VFXSgN9iX8O+KwHjzIgoIIsfuAZGVBRZFHEJIG4o\nKuOCg6iDCTCM4iCKjI4bIAwgyqKIqIhCENGICyioIEsSkE2QNREhIfX9uJ3QabqT7qQ73Z3U+zz9\nJPfcc8+pu/StPqfqVImqYhiGYRjBoly4BTAMwzBKF6ZYDMMwjKBiisUwDMMIKqZYDMMwjKBiisUw\nDMMIKqZYDMMwjKBiisUwDMMIKqZYDMMwjKBiisUIGBFZLyKX+Fk3VUS6+dj3ZxFZIyKHROS+4EoZ\negq7DgWde3H6Cma7/vQXwj5Cfh5GeCizikVE0kTkDxGp61G+RkRyRKR5uGQrLkF+oZ3Slqqeo6qf\nBqH5EcAyVa2lqv8pbmPeZBWRW0RkZXHb9obndQjlizKI1/wUQnyPjRASqcq5zCoWQIFU4KbcAhE5\nB6jq2hcQIlI+eKKVGWKBHwI9qAjXuszGLbLnsmBK+/UJ1/mVZcUC8D/gFrftW4Dp7hVE5CwRWS4i\nB0RknYhc47YvVURGiMh3QIaIlBORR0TkFxE5LCI/iUhXV90mIjJPRH4VkS0iMqwgwVxtjxSRH0Tk\nNxGZIiKV/JBpBtAcWOSS4eHC+nb1NVxEvhORgyLyhohU9taWW/1urv8fEZHNrv3rRaSPPxdeRD4B\nugIvuI49M9Br7U8/gSIig0XkXbftTSLyptv2NhFp4yZT7nXIvVbvuc7n7zgK7TzXdT3guq6VfPQ7\nQkTmepRNEpHnPPtycb63Z8Otrrfn8pT75Oc9blXIfcl9dk45x6I+H65jfT63fvRb2LGe16e9iHwr\nzrTsWyIyR0TGur4/8zzkel5EnvUhs9fvv1u/vr7ThX1HTxOR+a79e10yFHTv3M+vvDizMGe4tTdN\nRMa41X/YdS2PiMirItJQRN53tfuRiNTy974BoKpl8oMzWukG/AT8BUfJbgNOB3JcN6wCsAl4xPV/\nV+Aw8Ce3Nr4FmgKVgT+72mjk2t8caAkI8DXwOFAeaAFsBroXIt/3rrZrA58BYwqTye3Yrq7/C+3b\nVT8FaOTq60dgiGdbntfO9f91bud7PZCRu+1Z18s5Lgduc/0f0LUu6J56lA0GPg3guWgJ7Hf93wRI\nA7a5ts8AfvPVn+e1Kui6eum3uevaVXdtlwN2Ap28XHOvz4ZHv/muVUH3qaB77Od98XmOfvTr69ko\n8LktqF8/j3X/3lZ03ef7XPX7Asdwvm+NgSNAjOvY8sAeoJ0Xmb1+//34ThcmbzlgLTABqAJUAi4o\n5N553v8TwBludabhemZc9T8H6uM883tc8rRx9fUJMCqQ92tZH7HAyVFLdxwls9NtX2ecL/p4Vc1W\n1eXAe7hNnwGTVHWnqh7DuXmVgHNEpIKqblPVVKATUF9Vx6nqCVVNA14DbixEtsmutg8C41z9+iMT\nOA8rwPk++vasP0lV97j6WgS089LWKajqfFXd4/p/Ls5L6PxCzssbXfw4L/dr7Yt3RGR/7gd4IRAh\nXPfriIi0Ay4BlgA7ReTPru3C7DWe16qg6+re7zacl0FfV9FfgUxV/cpHP57Pxs1e+s27Vn7cJ1/3\n2N/vgNdzLMbz4c93xle/vp55z2Nzr09noLyq/sdV/21gtUvm3Tj3/HrXcT2Avaq61ovMvr7/7nj7\nThd2rv+H88Ifoap/qOpxVf3crU1v987zu+LzO+wm1z5V3eU63y9V9XtVPQ68DZxXyPH5qBBI5VLK\nTOBTnF+qMzz2NQW2e5SlA83ctn/J/UdVt4jIA0AS0FpEPgSG49gSmrledODc5HKufgviF7f/013y\nNPFDJnea+9n3Hrf/f3f1UygiMgh4EOdXFkB1nF8+geLPef1C4fR2vfxy5bsFuN1t+0zgXNfnPVX9\n1ksbK3B+mZ8JJAMHgHgc5bfCDxncCeS6voHzopnp+ju7gLqez4Znu/muVTHukz/fAZ/nWIx+/fnO\n+OrXn2fe/fo0BXZ49O9+ztOBu4EpQALOj9FT8Pj+ny0iS4Dhrpe1t35zv9OFyXsakK6qOd769YE/\n3xV33K/lUS/bNQJprMyPWFy/FFNxfoks8Ni9E2dqzJ3m5H8I8xmGVXWOql7sqgfwFM5DulVV67o+\nddTxhLqGgnHvO9Ylz063tv2Rqah9e2srH+J4zr0C3ONqtw6OMb6wX0feCPha+xKrkP3XuNp8FnjY\nR51PcRTJRTiK5FPgUpwRS0GKpbhOAnOBeBFphjNyKUixeHs2vMrix30qSG5/njevFPP5KM5z68+x\n7ue8i1N/mLlf33eANiLSGugJzPLVsdv3P9ZV9FQB7ebet8Lk3Q40F+92RV/3zrP8d6Ca23ZjX+cQ\nDMq8YnFxG85c71GP8i+B312GsAoiEo/zYL3hrRFx1mV0dRnkjuNo+hM4w+ojrnaquIxprUWkYyFy\n3SsizcRxiX4MmONqK9OLTHPcjtuNYw+gGH3nssetLU+q49ij9rkMoLcC5/jZricBXeuioqrPqupq\nnF+BntMUueSOWKqq6k6cqYErgXrAmgKad7/uRZFtn6vvaTgvmo0FVPf2bPiisPtU0D3+Eu/Pmz/3\npTjPh6/ntkMIjv0COCEi97rq9sZtus41nTQfR9F/qapeRwM+vv+eowxv9+0rH/LmfkdX4yi/p0Sk\nmjiONRe49hV079xZA9zsug9X4vxQChllWbHkaXRVTfWYElFXeRbOL9yrgH3Af4CBqrrJsw0XlXF+\noezF+SXSAHjMNYTtiTMHnAr8CrwKxBQi42zgIxxD3iZgXAEy/ex23FPAKNfQ+gE/+i7oF+uTuW2J\nyEPu9VX1J+AZHCPqbqA1jkHSnYLadr8HgV7rAtvzgz44c9ynNuL0eQTXVISqHgG2AJ+pqnsfnv3l\nXXcRGR6gPLnMxrGveP4q9uz3lGfDl1x+3KeC7nGR74sf/RZ0rK/vTK53UtCOdZ3jtcAdONOeN+PY\nbNxtedNxpk89p8vd8fb9f9SjjrfvdIHvB9f+a4A/4TgHbAf6u9rzee88eADo5Tq/m3DsJvioX9yR\nN5L/exJeRKQO8CbOEDEN6K+qh7zUexBn3jwHWAfcqqrHRSQRuBPnxoDzUv+wJGQPNiKSCtyuqsvC\nLUtpQxx32WSgsdsL0jDyEJEU4L+qOt21fTqOc09jVc0oYptl5jsdaSOWkcDHqvoXYBmnantEpCkw\nDGivqm1wHBDcvT0mqmp71ycqlYoROkSkLzAKZ2qjfyHVjTKCiFwiIo1c01C34IxOPnTtK4fjhDOn\nqEqlrBFpXmG9OTn3Nx3nV+VIL/XKA9VFJAfHIOVuuCyK4TgSiZyhZCnC5Ur6dqEVjbLGX4C3cN4n\nW4HrVHWPiFTDsWPkOvgUhzLznY60qbD9qlrX17Zb+d9w5pR/Bz5S1YGu8kScBXGHcBb4DPc2lWYY\nhmGEjhKfChORpSLyvdtnnetvLy/VT9F6IlIbZ2QTi+MDXkNEcheHvYizurQdjrFwYohOwzAMw/BB\niU+FqWp3X/tEZI+INHINQRtz0gjvzmU4rpj7XccsAC4AZqvqXrd6r+J4dnjrJ3KGaYZhGFGEqhZq\nbog04/27OFNZ4IRZWeilzjags8vfW3BcM38CcCmjXK4F1vvqSIsYYyzkn9at0e++K1YbiYmJ4T+P\nEH7s/KL7U5rPL2LP7dgxNDkZfewxtFMntGZNtEcP9N//Rr/5Bs3O9qsdf4k04/144C0RuQ0n3EF/\ncCJ/Aq+qak9VXS1OtNE1QJbr7yuu458WJ8ZTDo678l0lLH/xadQI9uwpvJ5hGEZBpKfDhx/CBx/A\n8uXwl7/A5ZfDv/8NXbpAJa+BtoNCRBnvSwoR0Yg975tv5qeWLXmzYsUiN5GcnEx8fHzwZIow7Pyi\nm9J8fuE8N8nJ4fRffuHPP//MnzZtokZGBkcuvJDGgwdD9+7QoEHx+xBB/ZgKM8USaTz4IJx2Ggwf\nXuQmSvMXF+z8op3SfH4lfm6ZmbB0KSxcCO+9B82aQc+ezqdTJygf3DxfplgKIKIVy1NPwf798PTT\n4ZbEMIxI5OBBR4ksWACffOIokN69oVcviI0t/PhiYIqlACJasUybBsnJMH16oVUNIxy0aNGC9PT0\ncIthhJDY2FjS0tJOKfdXsUSa8d4w470R4aSnpwfkIWREH47DbdExxRJpNGrEkS1beCYpKdySGIZh\nFAmbCos0fvkFzj8fdnrmbTKMyMA1HRJuMYwQ4usem42lACJasRw/DtWrw7FjUC7S1q8ahimWskBx\nFYu9uSKNSpWgZk347bdwS2IYZYYnn3ySIUOGBL1uYZQrV46tW7f6VXf06NEMHDgwKP2GGrOxRCKN\nG/PVe++x2DxvDCNgXn/9dSZOnMiWLVuoVasWffr04cknn6RWrVo+j3n00VNSPwWlbmEEaiT3t/6t\nt97K6aefzpgxY4oiVrExxRKJNGpEp+bN6XTrreGWxDBOYfTo0V7LU1PTGTXqdXbsyKFZs3KMHTuY\nli0DW1dR3DaeeeYZJkyYwIwZM+jWrRs7duxg6NChdO/enc8//5wKFU595Z04cYLyQV5I6C+ldkox\n7MHRwvBxTjuCueEG1Vmzwi2FYXjF2/dn69Y0jYsbrpChoAoZGhc3XLduTfO73eK2cfjwYa1Ro4bO\nmzcvX3lGRoY2aNBAp02bpqqqSUlJ2q9fPx0wYIDWqlVLp0yZoklJSTpgwIC8Y6ZPn66xsbFav359\nHTt2rLZo0UI/+eSTvONz66alpamI6PTp07V58+baoEEDHTduXF47q1ev1i5dumjt2rW1adOmet99\n92lWVlbefhHRLVu2eD2f1NRUvfTSSzUmJkYvv/xyve+++3TgwIF5+6+//npt3Lix1q5dWy+99FL9\n8ccfVVX1lVde0YoVK2rlypW1Zs2a2qtXL1VVfeqppzQuLk5r1qyprVu31rffftvntfT1jnSVF/qO\nNRtLJGJrWYwoY9So19myZTRQ3VVSnS1bRjNq1Osl1sbnn3/OsWPH6Nu3b77y6tWrc9VVV7F06dK8\nsnfffZf+/ftz8OBBbr7ZSeeUO830448/cu+99/LGG2+wa9cuDh06xE4PL03PKalVq1axadMmPv74\nY8aMGcPGjRsBKF++PM899xz79+/niy++YNmyZbz44ot+nc/NN99Mp06d2LdvH//4xz+Y7rFo+qqr\nrmLLli38+uuvtG/fPu887rzzThISEhgxYgSHDx9m4UInSPyZZ57JqlWrOHz4MImJiQwYMIA9IXrP\nmGKJRBo3hl27wi2FYfjNjh05nFQIuVRn586cEmtj37591K9fn3JevCmbNGnCvn378ra7dOnCNddc\nA0CVKlXy1Z0/fz69evWiS5cuVKhQoVA7hYiQlJREpUqVaNOmDW3btuW7774DoH379px//vmICM2b\nN2fIkCGsWLGi0HPZvn07X3/9NWPGjKFixYpcfPHFefLmMnjwYKpVq0bFihV54okn+O677zhy5IjP\nNq+77joaNWoEwPXXX8+f/vQnVq9eXagsRcFsLJFIs2bsWbaM/9oiSSNKaNasHJBJfsWQSdOm/v92\nLW4b9evXZ9++feTk5JyiXHbt2kX9+vXztk8//XSf7ezcuTPf/qpVq1KvXr0C+859YQNUq1aNjIwM\nADZt2sRDDz3E119/zdGjR8nOzqZDhw6FnsvOnTupU6cOVatWzSuLjY3ll19+ASAnJ4fHHnuMefPm\nsW/fPkQEEWHfvn3UrFnTa5szZszg2WefzQvVkpmZmU/ZBhNTLJFI06Y0OnGCJFMsRgTizXg/duxg\nUlIS3aayMomLS2Ts2GF+t1vcNrp06ULlypVZsGAB/fr1yyvPyMjggw8+4KmnnsorK8i7qkmTJvz8\n889520ePHuW3Irr/Dx06lPbt2/Pmm29SrVo1Jk2axPz58ws9rkmTJhw4cICjR4/mKZdt27blKcxZ\ns2axaNEili1bRvPmzTl06BB16tTJcwbwPL9t27YxZMgQli9fTpcuXQA477zzQuY8YFNhkUizZrBj\nR7ilMAy/adkylqVLh5GQMIGuXRNJSJjA0qXDAvLoKm4bMTExPPHEEwwbNowlS5aQnZ1NWloaN9xw\nA82bN2fAgAF+tdOvXz8WLVpESkoKWVlZhf7AK+jlfOTIEWJiYqhWrRobNmzgv//9r18yNG/enI4d\nO5KYmEhWVhafffYZixadzLSekZFB5cqVqVOnDpmZmTz66KP5lEmjRo3yrY/JzMykXLly1K9fn5yc\nHKZNm8b69T4T7BYbG7FEIk2bWkgXI+po2TKWmTMTw9rG3//+d+rXr8/DDz/M1q1biYmJoW/fvsye\nPZuKfibPO/vss5k8eTI33HADv//+Ow888AANGzakcuXKXut7jg7ctydMmMCQIUN4+umnOe+887jx\nxhtZtmyZz2PdmT17NoMGDaJevXp06dKFW265hYMHDwIwaNAglixZQrNmzahXrx5jx47l5Zdfzjv2\n9ttv5/rrr6du3brEx8ezYMECHnroITp37kz58uUZNGgQF110kV/XoyhYSJdIRNVZfb9zJ8TEhFsa\nw8hHWQvpkpmZSe3atdm8eTOxIc53EilYSJfSiIiNWgwjjLz33nscPXqUzMxMhg8fTps2bcqMUgkG\nNhUWqTRrxtrFi3lnzpxwS2IYZY6FCxfmxeXq2LEjc+x7GBA2FRapJCTAlVdClASdM8oOZW0qrCxi\nU2GlFfMMMwwjSjHFEqmYjcUwjCglohSLiNQRkY9EZKOILBERr3GuReR+EVnn+vwt0OOjAhuxGIYR\npUSa8X4k8LGqPi0ijwCPusryEJHWwO1ARyAb+FBE3lPVrf4cHzU0bcqhDRt41lbfG4YRZUSU8V5E\nNgCXquoeEWkMJKvqWR51+gFXqOqdru1/AH+o6gR/jncdE/nG+7Q0uOQS2LYt3JIYRj7MeF/6KW3G\n+4aqugdAVXcDDb3UWQ9c7Jr2qgZcBeRGjGvkx/ERS2pqOgMGjKZr10QGPzoF3b0bcvyPDmsYhne2\nb99OTExMSBXirbfeyhNPPOFX3fT0dMqVK0dOKf1+l7hiEZGlIvK922ed628vL9VPeQpUdQMwHlgK\nvA+sAU746C5qflalpqbTvftkZs16mOTk0UyfM5IDJyqQ/vW34RbNMKKGFi1aUK1aNWJiYqhZsyYx\nMTHs3r2b008/ncOHD+eFUOnatStTp07Nd2wg+eeDgb9phlesWFFgNOZIpMRtLKra3dc+EdkjIo3c\nprJ+9dHGNGCa65hxwHbXrt3+HA/kCywXHx9PfHx8oKcSVLwlOdqeE8dbSa8x7v2O4RTNMKIGEWHx\n4sV07dq1SMdGIqoaNtmSk5NJTk4O/EB/0kyW1AdnJPKI6/9HgKd81Gvg+tsc+BGICfB4r2k3w0l8\n/BOudKwnP4vpoY+ee2O4RTOMfETi9ycX9xTC7uSmED5x4oQ+/vjjWr58ea1atarWrFlThw0bppdc\ncomKiFavXl1r1qypb731lqqqLlq0SNu1a6e1a9fWCy+8UL///vu8Nr/99ltt3769xsTE6A033KA3\n3nijjho1yqtcJ06c0OHDh2v9+vU1Li5OX3jhBS1XrpyeOHFCVVWnTZumrVq10po1a2pcXJy+/PLL\nqqqamZmpVatW1fLly2uNGjW0Zs2aumvXrkJTHhcXX/cYP1MTh12Z5BMG6gIfAxuBj4DarvImwHtu\n9T7FsbWsAeILO95LP0W93iEjISHJLde383mVW/S183uGWzTDyEckfn9yKUixuL/I4+PjdcqUKfnq\niIhu3bo1b/vbb7/Vhg0b6ldffaU5OTk6Y8YMbdGihR4/flyPHz+usbGxOmnSJM3OztZ58+ZpxYoV\nfSqW//73v9qqVSvdsWOHHjhwQLt27ZpPnvfff19TU1NVVfXTTz/VatWq6Zo1a1RVNTk5WU8//fR8\n7X3zzTf65Zdfak5Ojqanp+vZZ5+tkyZNKtpF80JxFUtEuRur6n7gMi/lu4CebtuXBHJ8NOAtydHR\nOhuJP62pJfwyootgTdsU0dDep08fKlRwXm25IeP97/Jkn6+++ip33303HTs6U9EDBw5k3LhxpKSk\nAJCdnc3f/uYso7vuuuvo1KmTz3bnzp3LAw88QNOmTQF49NFH86Uo7tGjR97/F198MZdffjkrV66k\nXbt2Xttr37593v/uKY9z5Qk3EaVYyjK5SY5GjZrAzp05NG1ajptaXUP9tFTONcViRBDeMkjmI4Se\nV/6wcOHCItlYPElPT2fGjBlMnjwZcJROVlYWO10RMZo1a5avfkHRjz3THXvW/eCDDxgzZgw///wz\nOTk5HD16lDZt2vhsr6gpj0uKSHM3LtPkJjlatmw0M2cmUr9tG1t9bxgBon4oNn+M4aeffjqPP/44\n+/fvZ//+/Rw4cICMjAxuuOEGmjRpwg6P7+a2AtacNWnShO3bt+dtp6en5/1//Phx+vXrx4gRI9i7\ndy8HDhygR48eeefhTdahQ4fSqlUrtmzZwsGDBxk3bpxf511SmGKJZCxemGEEDfcXr2fqXoDGjRvn\nK7vzzjt56aWXWL16NeAk/Hr//ffJzMykS5cuVKhQgcmTJ5Odnc2CBQvy6nmjf//+PP/88+zYsYMD\nBw4wfvz4vH3Hjx/n+PHj1K9fn3LlyvHBBx/w0Ucf5ZP1t99+4/Dhw3llRU15XFKYYolkmjUzxWIY\nAVDQSMR93/3338/cuXOpV68eDzzwAACJiYkMGjSIunXrMm/ePDp06MCrr77KfffdR926dfnzn//M\n9OnTAahYsSILFixg2rRp1KtXj7lz53Ldddf57PvOO+/kiiuuoG3btnTs2DFf3Ro1avD888/npRKe\nM2cOvXv3ztv/l7/8hZtuuokzzjiDunXrsnv3biZMmMCsWbOIiYnhrrvu4sYbbyzyNQsFERXSpaSI\nipAu4Ky6r1IFjhwBH/m2DaOksZAupZ/ihnQx430kU64cNG7MF2+/zZING8ItjWEYhl/YiCXS6dwZ\nJk6ECy4ItySGAdiIpSxQ2oJQGp6YncUwjCjDFEuk07SpuRwbhhFVmGKJdGzEYhhGlGGKJdKxEYth\nGFGGeYVFOjZiMSKM2NjYiA0xbwSHgsLT+IN5hUU6Gzbw++WX8/Rtt4VbEsMIGg327uWmN94g87LL\nOH32bKgQ2b9xBwwYzaxZD3MyXxJAJgkJExg7djCjRr3Ojh05NGtWjrFjB9Oy5akv5oLamDkzMaTy\nBwt/vcLCHio/HB8iOOz3KRw5olq1qmpOTrglMYzgsHixaoMGqq+/Hm5J/MZbviRQ7dz5fo2LG+6W\n8iJD4+KG69ataae0sXVrmt91IxX8DJtvNpZIp0YNqFoV9u0LtySGUTxUYcIEuOMOWLgQbrkl3BL5\nTbNm5YBMj9JMdu/efkrm1y1bRjNq1OuntJEbwTwhYQJduyaSkDCBpUuHeR3dRDuRPf40HGJjIT0d\nGjQItySGUTSysmDoUPj6a/jyS4iyHO7e8iXFxSXSoMGZpKVV96hdnZ07c7y2kxvBvLRjI5ZooHlz\nKCAkt2FENAcPQo8esGcPfPZZ1CkV8D3aiIurhreRTNOmZfvVasb7aOD++6FFC3jwwXBLYhiBkZYG\nV10F3bs7oYnKlw+LGKmp6V4N7L7KA2m3e/fJp4xkSusUlwWhLE3YiMWIRr75Bnr1gpEjYdiwsInh\n7eWfkpLI1Kl9ue22t08pD0QpeMv8OnZs6VQqgWAjlmhg3jz2Pv88L3TrFm5JDMMvzty8mb5vv03q\nyJGcM2pUWGXx5ebbosUg0tJmnFIeTe6/JY2NWEoTsbE0yMwkKSkp3JIYRuG8/jq89BIsXco5ERCV\ne8eOHPIrD4DqHDxY3Wu5L8O74T+mWKIBmwozooWnn4YXX4TkZDjrrHBLA7i7CucfmdSuncnBg6eW\nl3XDezCwqbBoQBWqVXPWslT3/IVlGBFATg488gh88AEsWeKEIooQfBnYvdlYcg3vQLGM+qUVf6fC\nIkqxiEgd4E0gFkgD+qvqIS/17gfucG2+pqqTXOWJwJ3Ar659j6nqh16Ojy7FAvDnP8O770bMr0DD\nyCM721n0+PPP8N57ULduuCU6hVzvr5MG9vxeYe7lQJny9AqEaFUs44HfVPVpEXkEqKOqIz3qtAbe\nADoB2cCHwF2qutWlWI6o6sRC+ok+xdK9Ozz8MFxxRbglMYyTHDsGN90ER4/C/PnOyDrKKQ0xvUJF\ntGaQ7A1Md/0/HejjpU4r4EtVPaaqJ4AVwLVu+0tn2NXmzZ3V94YRKfz+u+NOXK6cE6KlFCgV8G3s\nN6O+/0SaYmmoqnsAVHU30NBLnfXAxSJSR0SqAVcB7kt57xORtSLymojUCr3IJURsrBnwjcjh0CFn\n9NykCcyZA5UqhVuioOErLpgZ9f2nxKfCRGQp0Mi9CFDgH8DrqlrXre5vqlrPSxu3AvcCGcAPwDFV\nfUhEGgD7VFVF5J9AE1W93cvx0TcVNn06u2fN4qUIcN80yjZVjh5l4MyZcP75NJs/3xmxRBC2mj50\nRKuN5ScgXlX3iEhjYLmqtirkmHHAdlV9yaM8Flikqm28HKOJiSfnSuPj44mPjw/GKYSO5cshKQlW\nrAi3JEZZ5rffHHtffDw88wyEMeGXNwUCwTG8+zL2lzWSk5NJTk7O2x49enRUKpbxwH5VHe/LeO+q\n10BV94pIcxzjfWdVPSwijV1TaIjIg0AnVb3Zy/HRN2LZuhW6dXNiLxlGONi7Fy67zAko+eSTYVcq\n3hRI69bCu+8mYYb30BCtK+/HA2+JyG1AOtAfQESaAK+qak9XvfkiUhfIAu5R1cOu8qdFpB2Qg+Ou\nfFdJCh9STjvNSVGcnR3x2faMUsivvzo/bPr2hTFjSlSpeBuZjBr1utc8KEeODMUM7+Enot5Qqrof\nuMxL+S6gp9v2JT6OHxQ66cJMpUpOPpZdu6Iy7LgRxeSOVMKkVLwFkGzQoDzeFIhIBt5W2ZvhvWSJ\nKMViFEKLfn/YAAAgAElEQVSuy7EpFqOk2LfPUSo9e5aIUvEcnRw5kuF1ZHLixCC8KZDOnWNZv/7U\nhFxjx4YvunJZxBRLNGEux0ZJsn+/Y6jv0QPGjSsRpeI5OqlSZQjeRiaNG59O+fKnKpBnn30AwMLY\nhxlTLNFEbopiwwg1hw/DlVc6dpUQGOr9tZv88ccZeBuZxMXVYfbswT4ViBnqw0tEeYWVFFHpFQbw\nwgvsXLKEV9q3D7ckRimm4vHjDJg1iwrt2tHsnXdColS8eXQ1aFCelJTxHrXTqVp1NEePTsbWlISf\naPUKMwoiNpam2dmWl8UIHX/8AddcA127wtSpIZn+8uXR5d1uUp/LL69FjRo2tRVNmGKJJixemBFK\nsrLghhuc6MRTpgRlRb23KS9fsbgKspuYIokuTLFEE7nGe9WwLk4zSiE5OXD77Y5ymTsXypcvdpO+\nXIVbtxaKYjcxogezsUQbtWpBampE5rwwohRVuP9++PZb+OijIkcp9uYq7G0VfO/e/2D9+vIWiysK\nMRtLaSXXM8wUixEsxoyBTz910gkXQ6n46yp8+HAMS5feZiOTUowplmijRQtnxHLeeeGWxCgNvPQS\n/O9/sGoV1K5d5GYCcRVu2rQcLVvGmktwKaZIikVEqgN/uBJtGSXJmWfCli3hlsIoDSxYAGPHOqOV\nRo0Kr+/Cf4P8HVStOuwUV2FbBV/68UuxiEg54EYgAScl8DGgsojsAxYDL6vq5pBJaZwkLg7WrQu3\nFEa0s2IF3H03LFniPFN+EphB3lyFyyp+Ge9FZAXwMbAQWK+qOa7yukBX4GbgbVWdGUJZg0ZUG++X\nLOHAY48x6Zprwi2JEaU03LOHQf/7H5vHjKHdQw8FdKyvfPBmkC8bBNt4f5mqZnkWuqIRz8cJY18x\nQBmNohAXR539+22RpFE0fvkFLrgAXn2VdjfdFPDhvtagmEHecMcvxeKuVESkDvAnoIrb/k+9KR4j\nBMTGOnlZjh8vVXnGjRLg8GG4+mq47z7wQ6l4s6WczAdvBnnDNwEZ70XkDuB+4DRgLdAZ+ALoFnzR\nDK9UrOgk/UpLgz//OdzSGNFCVhb06wcXXgh//3uh1X3ZUqZO7UtKioWlNwomUK+w+3GM9ymq2lVE\nzgL+FXyxjAKJi3M8w0yxGP6g6hjqK1eG55/3K2qDr3her7wygaVLh9mUl1EggSqWP1T1DxFBRCqr\n6gYR+UtIJDN8c+aZsNmc8Aw/GT8e1qxx3Ip9pLX2nPbasuV3fKX4tSkvozACVSy/iEht4B1gqYgc\nwMlNb5QkuSMWwyiMefPgxRfhiy+gRg2vVbxNe9WocROW4tcoKgE9JaraV1UPqmoSMAqYAvQJhWBG\nAdiIxfCH1ath6FBYuBCaNfNZzdu0V0bGeGrUGIajXOCkLWVwSEU2SgdFDumiqiuCKYgRAHFxZK5b\nx7/N5djwQa1Dh7h9yhTSR47k3ELC/3h3IW7FOefEEBdnthQjcPxdeT8Z8LmiUFX/FjSJjMI54wyq\n79lD0qhRQQlvbpQyMjLgoosgMZFzhw/PtysQF+K4uDpmSzGKhL8jlq/d/h8N2NMWTqpVg3r1YMcO\nJ/mXYeSSkwODBkH79uCxqt5ciI2Swt8FktNz/xeRB9y3g4mI9AOSgFZAJ1X91ke9K4HncGxEU1R1\nvKu8DvAmEAukAf1V9VAoZA07uXYWUyyGO4mJ8Ouv8MYbp7gVmwuxUVIUxcYSyiBb64C+wMu+KrgC\nYv4H+CuwE/hKRBaq6gZgJPCxqj4tIo8Aj7rKohZvUxctW8ae9AzrZmtTDRdz5jgh8FevJnXnbr9T\nApsLsRFsIiofi6puBBApcAXX+cAmVU131Z0D9AY2uP5e6qo3HUgmihWLr6mLpUuH0dLC5xvurFkD\nw4bBxx+Tmnk0oJTA5kJsBBu/nigROSIih0XkMNAm9//c8hDL6EkzYLvb9i+uMoBGqroHQFV3Aw1L\nWLag4mvqYtSo150Ri7kcGwB790LfvvDCC9C2rc/nRiSbuLhEzIXYCDX+2lhqBqtDEVkKuGcVEpzp\ntcdVdVGw+nERpbHxHQqaurBFkgbgxAC7/nq4+Wbo3x+wCMRG+PHX3bjQBCb+1AFQ1e7+CueDHYC7\nxfo0VxnAbhFppKp7RKQx8KuvRtzDzsfHxxMfH19MsYJPQZFk80Ysqn7FfjJKKcOHO16CY8fmFVkE\nYiNYJCcnk5ycHPiBqlroB8dWMQxo7lFeCSey8XRgsD9t+dnfcqCDj33lgc04nl+VcKIst3LtGw88\n4vr/EeApH21oNLB1a5rGxQ1XyFBHg2RoXNxw3bo1zalQt67qnj3hFdIIH9Onq555pqat/U4TEpI0\nPv4JTUhI0hUrPiv4uTGMIuJ6dxb6Dvc3g2QV4Dac1MQtgYM4+VjKAx8BL6rqmsDV2in99AEmA/Vd\nfaxV1R4i0gR4VVV7uupdCUzipLvxU67yusBbwOk4Mcz6q+pBL/2oP+cdCeR6hZ2cuhh8curi//6P\nbwcO5N19+8Iqo1HyNN61i4EzZ7IiMYm/vbLnlDUoU6f25ZVXPvb+3BhGEfE3g6RfisWj4Yo4L/6j\n3l7a0UA0KZYCuflm6NEDBg4MtyRGSbJ/P3TsCE8+yYBFG7ymCk5ImGBTXkbQ8VexBOxnqKpZqror\nWpVKqcKCUZY9TpxwflD07Qs33FCwg4dhhAlzYI9mTLGUKVJT03m7XTfWpmxg0K7qpKamuxnq3bG1\nKUZ48XcdyxQR6ScizVzbtVx2FyOcnHUW/PRTuKUwSoDU1HTGXvAQndZv5spDKfzvjUfo3n0yQ4Zc\nZmtTjIjD35X3u1V1ntv278ClItIKSFfVd4MvmlEoZ50FGzc6gQfL2S/U0sxzDz7Hk7tX0o/57KEx\ngMX5MiIWfxVLKoCIXI0TIHI1sFxVPxaRTwBTLOEgJgbq1oVt26BFi3BLY4SKY8cYuuwtnuYRPuNi\ntx0W58uITALNILkYaAxsxFkxD/B2sIUyAqBVK/jxx3BLYYSShx7icO1aTOQujx1mSzEiE3+fyidF\nZJaIDMVZX7JXVbNd+34PjWiGX5x9ttlZShmpqekMGDCarl0TmXzR9WQtXkzDRW8SF5eE2VKMaMDf\nqbDHcaa//g9ngeSXIpINfAfUA6aGRjyjUFq1YtfChbx85Ei4JTGCwIEDh/jf/w5x4MBkzmQnN3AB\nPet35aWYGLOlGFFDwAsk8w4UqYETwv5+Ve0dVKlCTKlZIAnw6acwciR8/nm4JTGCwIABo5k162Gq\nUI4v6MLL3MVLDLIFj0ZE4O8CySLnY1HVDGCZiNhP5XBy9tmOjcWCUZYKchc8PsvdbOAsXuJuQGzB\noxFVFDvRl6p+FQxBjCJSvz5UrAi7d0OTJuGWxigmzZqVoz/T+Suf0IFvcHxkzEhvRBcRlUHSKCK5\nBnxTLFGHZ+rpv13dirg3BtM9ZylHiOGkkX5YuEU1DL8xxVIayHU57tYt3JIYAeCZeroiBxgxrxUH\n7h3C2fs/pvbOj8xIb0QlASkWERkGzFTVAyGSxygKrVqZy3EU4plC+CnGsvVYJ+b9VpuZs8xQb0Qv\ngU7cNgK+EpG3RORKEbMWRwS5BnwjqnCPTHwVi7mO+dzGdHbuKiUei0aZJdCV9/8A/gRMAQYDm0Tk\nXyISFwLZDH+xEUtUkhuZuAk7mcLtDGAmB6hshnoj6inSOhYRaQvcClyJk0a4M7BUVUcEV7zQUKrW\nsYDjalyrFqSlObHDjIjD00ifu2L+8sue579bv2Ul8YzhYeLiElm61GwqRmQSkgySInI/MAjYB7wG\nvKOqWSJSDtikqlExcil1igWcNMUDBvDub7+FWxLDA/fV9Lnpg+vUGcbAgbX469df0fKHLTzY7g4a\nn1bBUggbEU2oFkjWBa5V1XT3QlXNEZGeAbZlBJNWrWhftSrtk5LCLYnhwYABozlw4J+czPRYnQMH\nJlNr4wP02vQzfP81HzdvHk4RDSOoBKRYVNWnq4qq2iR/ODEDfsTiLX1wDXK469P5MOMVMKVilDIC\ndTd+yEvxIeAbVV0bHJGMItGqFSxfHm4pDC+cTB98Urk8zz1satqCZv36hU0uwwgVgbqfdATuBpq5\nPnfhGPBfFZGoMNyXWix8fsQyduzgfOmDr2cG8RUX0+KdmWGVyzBCRaDG+0+Bq1wBKHMjHC/GUS7f\nqOrZIZEyyJRK4/2JE05Gyd27oWbNcEtTZvHm/dWyZWxeefbWA7z87atkzHmDZn16hVtcwwgIf433\nqKrfH2ADUNFtuzKwwfX/mkDa8tF+P2A9cAJoX0C9K12y/Aw84laeCPwCfOv6XOnjeC2VnH++6sqV\n4ZaizLJ1a5rGxQ1XyFDHBzxD4+KG69ataU6F7GzVSy9V/de/wiqnYRQV17uz0Hd5oFNhs3CSfCWK\nSCKwCpgtItWBYFiO1wF9gRW+Krhcm/8DXAG0Bm4SkbPcqkxU1fauz4dBkCl6aNcO1pqpK1x4hmiB\n6mzZMppRo153NidMcPTNCJs1Nko3fhvvXeFbXgc+AC50Fd+tql+7/k8orjCqutGtL1+cj7NmJt1V\ndw7QG2cEA06c8bLJeefBV5bFIFx48/6C6k4ulTVr4JlnnPtTvnw4xDOMEsPvEYtrGPS+qn6tqpNc\nn68LPTD4NAO2u23/4irL5T4RWSsir4lIrZIVLcycd57zAjPCwknvL3cyiW14AhIS4NlnIdYWPxql\nn0AXSH4rIp20GMm9RGQpTjDLvCJAgcdVdVFR23XxIjBGVVVE/glMBG73VjHJbSFhfHw88fHxxew6\nAjj3XNiwAY4fh0qVwi1NqcabkX7s2MGkpCS6TYc5uVSeq/IbtG0LN98cXqENI0CSk5NJTk4O+LhA\nvcI24AShTMP5aSY4g5k2AfdccD/LgeGq+q2XfZ2BJFW90rU90iXDeI96scAib7KVSq8wF5ktWjDj\nyivZ07hxuEUptRQUogVg+fJtHDlSnZo1Mxl6Rg43L1/Gupkzueiaa8Iqt2EUl1CFdLmiiPIUBV/C\nfwWc6VIcu4AbgZsARKSxqu521bsWx8OsTFH9wgsZ2rkzDB4cblFKLb5CtPz22wRmznQLTvGba6Ty\n9ttcZEnYjDJEoF5h24CLgVtcxnMl/7RWsRCRPiKyHSda8nsi8oGrvImIvAegqieA+4CPgB+AOXoy\nnMzTIvK9iKwFLgUeDJZsUYPZWUJOgUb6XFRh6FDo398yexpljkBHLC8COUA3YAxwBJgPdAqGMKr6\nDvCOl/JdQE+37Q+Bv3ipNygYckQ1550H774bbilKNd5CtEBm/jwqs2fDDz/AjBklLJ1hhJ9AbSzf\nqmp7EVmjque5yr5T1bYhkzAElAYbi68V3vz2G5xxBhw4AOUsYVQo8MxVn2ukz8ujsn07dOgAS5Y4\nit4wSgmhsrFkiUh5nCkwRKQBzgjGKEG8vdhSUtxebLVqQWoqxEVFepyIx5sSX7p0GKNGTWDnzhya\nNi3H2LGua5+T49i3HnjAlIpRdvFneb6eDIWSALwL7ADGARuB6wNpIxI+RHlIl4SEJLewIZoXPiQh\nIcmp0KuX6ty54RWylFBomBZPJk1S7dJFNSurZAU1jBKAUIR0UdVZwAjgX8BOoI+qzg2emjP8oVDj\nsRnwg0ahYVrc2bABxoyB6dOhQqCTAYZReghIsYhIZaA9UAuoB1wvIk+EQjDDN75WeOcZj887z2KG\nBQm/PMAAsrJg4EAYOxb+9KeSEs8wIpJArbsLceJyZeO82XI/Rgnimd8j13g8duxgZ7NdOxuxBIlC\nlXguTz4J9erB3XeXlGiGEbEE6hW2XlXPCaE8JUJp8go7aTx2eYUBqJJVqxbPDRlCZo0aYZUz2ilo\nlX2dOs5K+yY7dzJg1iy+e/11Lrj++rDKaxihxF+vsEAVyyvAZFVdVxzhwk1pUCyF8te/wt//Dlde\nGW5JoobCknR5VeJ//OG4Fj/2mBNo0jBKMaFSLD8CZwKpwDFCFCss1JQJxTJihJNJctSocEsSFRS6\nNsUXf/+749o9dy4UmO0hNFx99dW89tprNGnSpMT7NsoeoVrH0qOI8hglzQUXwCuvhFuKqMG395dH\n/C93Vq6EWbPgu+/ColTeeustPvzwQ7Kyskq8b8MoCL+M9yIyAkCd+GDnq2p67ge4K5QCGkWkSxf4\n4gtnwZ5RKH57f+WSkeEshHzpJWjQIMTSncqRI0dYtWpVifdrGP7gr1fYjW7/P+qxzybxI5FGjaBu\nXWdthVEofnt/5fLww3DxxdCrV6hF88pLL73EkCFDKPVTukZU4q9iER//e9s2IoULLoDPPw+3FFFB\noS7c7ixZAu+/D889l6946tSpTJgwgTvvvDOvrG3btrzzzilxVYvFN998Q6tWrahhHn9GhOKvjUV9\n/O9t24gUchXLHXeEW5KIwpf3l8/4X+4cPOhcz2nToHbtvOLvv/+eWrVqcfHFF9OlSxdeffVVAC65\n5BLS0tLy6mVnZ3PPPfeQnZ0NcMqIw2UcRUS48cYbufzyy/Ptz8nJ4b333iMxMZH09HQkDLYdwygM\nv7zCROQEJzNGVgV+z90FVFHViiGTMASUCa8wcIzKN9xg02FuFNn7K5dbbnG87f7zn3zFn3/+OV26\ndCEpKYkdO3bw2muvAbBs2TKys7NPURBFZcqUKVx22WXExsaSnp7OGWecQWpqKs2bNw9K+4ZREP56\nhfk1Faaq5VU1RlVrqmoF1/+521GlVMoU55wDO3c6ofQNIMDYX568844zAhw//pRdF1xwASLClClT\nGDToZFqgtLQ0Lr300qDIvmvXLv744w9iYx0FWCZ+HBlRiUXKK82ULw//93+QkgJXXx1uaSKCgL2/\nctm7F+65x1mvUt3zeIeffvqJ3377jYsuuiiv7NChQ1SuXDlvOysri3vvvTdvKswbvqbCPvnkE1JT\nU3nssccA2LNnDwDjx4+nc+fODBw4sOBzMIwSIqAFkqWFMjMVBpCYSPqWLUw788xwSxIRLFiwjnXr\nZuCZ/fHccwdx7bXnej9Ilf5z51KjTRuaz5njs+0ffviBbt265b3wFyxYQNu2bYkLUV6cFStW0LVr\nV9LS0mwqzCgRQrLyvrRQphTLkiXw1FOwfHm4JYkIimRjmT0b/vUv+PprqFKlwPZHjBhBlSpVaNq0\nKa1atQraNJgn//73v1m8eDErV67k6quvpm/fvtx6660h6cswcjHFUgBlSrEcPAinnw7790NFM4dB\nIQE8Pdm504kW/cEHTkywELNlyxZee+01qlevzqFDhwAYN24clSpVCnnfhlEYplgKoEwpFnCM+NOn\nl8iLMdLw5VrsF6pw1VXQuTMk+gjrEmQ6duzIiBEj6N+/PwCdOnWiffv2vPzyyyXSv2EURKhihRnR\nSO56ljKmWLxNe6WkBOBa/NprjtHeZSwvCVSVlJSUPMUSFxfHypUrS6x/wwgGgSb6MqKRMroCv1iu\nxampjkKZMaNEpxC/+eYbJk6cmLe9efNmOpSxHwRG9BNRikVE+onIehE5ISLtC6g3RUT2iMj3HuV1\nROQjEdkoIktEpFbopY4CLr4YkpOdqZ0yRJFdi3NynACTI0fC2WfnFZdUyJZcUlJS2Lt3bz5FYxjR\nQKRNha0D+gKFTShPAyYDMzzKRwIfq+rTIvIITsDMkUGXMtqIi4OqVWH9ejjXh0ttKeRkYMn8rsU+\nA0vm8txzjhJ+4IG8opIK2QKwdetW3nnnHebNm8esWbNoEIboyYZRLFQ14j7AcqB9IXVige89yjYA\njVz/NwY2+DhWyxx33aX6zDPhlqJE2bo1TePihitkqKMpMjQubrhu3Zrm+6AfflCtX191y5Z8xatW\nrdKcnBx94okn9Pbbb88r/+STT3TJkiUhkf/w4cPasmVLfeKJJ0LSvmEEiuvdWeg7PKKmwoJAQ1Xd\nA6Cqu4GGYZanxElNTWfAgNF07ZrIgAGjSU1Nd3Zcfjl89FF4hQsh3s47N7BkQsIEunZNJCFhQsGG\n+6wsGDQIxo2DM87It6skQrZ4UrNmTXr16sW4cePYtm1bSPowjFBQ4lNhIrIUaORehBMh+XFVXRTk\n7nwaFZKSkvL+j4+PJz4+PshdlzwFekF16+YEUPzjj0IX+UUbhXl/+cwA6cnYsU4eGzcbijuhDtny\n+eef06dPH6ZMmcI111wDQKVKlVDVvDUthlGSJCcnk5ycHPiB/gxrSvpD0afCfiL/VNhPPo4t3ngw\nQklISHKb9tG86Z+EhCSnQufOqh9/HF4hQ0Ch5+0PKSmqjRqp7tzps8r69eu1YcOGedvz58/XzZs3\nF0f0fKxdu1ZPO+00TUlJUVXVnJwc7dSpk1500UVB68MwigN+ToVFmvHencIW4YiXOu8Cg4HxwC3A\nwuCLFbkU6gV1+eVse+01ppaydRGrVm3F23mvWrUl38jUFxWPH+ful19m99/+RusmTXzWa926Nbfc\ncgtPPPFEXsiWYMYBa9u2LVOnTiU5OZnFixeTmppKhw4deOqpp4LWh2GUCP5on5L6AH2A7cBRYBfw\ngau8CfCeW73ZwE7gGLANuNVVXhf4GNgIfATU9tFPMJV4xFDoL/fPPlNt1y68QoaAYo9Y7rlHdcCA\n0AppGKUA/ByxWEiXUkShARazsqBBA9i40bEllBKKlbzrww9hyBD4/vt8GSFDwd69e5k0aRInTpxg\n7dq1dO7cmccff5wKFSJ54sAwTmKxwgqgtCoW8CPAYp8+0L8/3Hxz2GQsDr5ifwUUWDKXffugbVuY\nORO6dg257EOGDGHSpElUrVqVP/74gw4dOnDZZZcxadKkkPdtGMHAFEsBlGbFUigvvgirV8Prr4db\nkoApdlphd1Thuusct+IJE0Ihbj42b95M7969mT17Nm3btgXgqaeeIikpicOHD1v0YiMqCGpqYqMU\nkbueJQoVa7Fif3ny+uuwebOzZqUEqFSpEnv37uXnn3/OK6tRowZZWVkcPny4RGQwjJLCJnfLGnFx\nzjqWdeugTZtwSxMQRY795cnWrTBiBCxbBq41KFOnTmX//v1s3LgxL1xL27ZtGT16NH369Cm27M2b\nN+fXX3/NV7Z69Wpat25N/fr1i92+YUQSpljKGiLQty/Mmxd1iqXIsb/cyc6GAQOcyMWuuGklGQcs\nl7S0NBYsWMBHpTgaglF2MRtLWSQlBW69FX780VE0UUJQbCxPPOHYmN5/H8o5Cunzzz+nS5cuJCUl\nsWPHDl577TUAli1bRnZ2doEKoihkZ2fTvXt37rzzTm6OUicKo2xixvsCKPOKRRViY52X6znnhFsa\nn3jzAAMC9/7K5dNP4YYbYM0aaNz4lN2nnXYas2fP5pJLLgGc6bGEhIR8IVuCwb333ssVV1xBr169\ngtquYYQayyBp+EYErr8e5s6NWMVSUPwvv2N/uXPgAAwc6GSF9KJUQh0HLJdnnnmGnj170qNHDwDe\nfPNNevbsSfXqnrYjw4hebMRSVklJgdtuc6bDIpABA0Yza9bDeNpTEhImBK5YVJ2RSuPG8PzzXqv8\n8MMPdOvWjT179gCwYMEC2rZtG9SQLbNnz2b16tV5SkVVmTNnDq9Hoeu3UTaxqbACMMUCqPJH48ZM\nue469jaMvOwC06dvJS3NM48btGgxiFtuOcPLEb7p8M03dPrqK36eMYNLCrCXjBgxgipVquTFAQtm\nOPxNmzZx7rnnkpWVla/84osvLlr0WMMIA6ZYCqCsKRZfq9V56CGIiQE/AjWWNEEbsaxbB926wcqV\ncNZZwRbTMMoUplgKoCwplgI9qXbvhNtvD/t0mC8jfbE9wDIzoWNHePRRJ4FXCfD+++8zbdo05s6d\nWyL9GUZJYoqlAMqSYinwl/+MUY532IcfQuvWYZGvIMUHxfAAA8elWrVEwtcsXLiQTz/9lHXr1pGd\nnc2yZctC3qdhlDTmFWYAhaxWL1cO+vWDN9+EMWPCIV4BYVqcKa8ieYABzJgBX3wBX38dLFELpHfv\n3vTu3ZvRo0ezYsWKEunTMCIVixVWyjm5Wt0dt9Xqt94KU6c6IfXDQNDCtLjz/fcwfLgTXaBGjeKI\nZxhGETDFUsoZO3YwcXGJnFQuzlRTrh2DNm2gZUtYtCgs8hWq+ALl0CEnavFzz+Wt0Zk6dSoTJkzg\nTrdc9m3btuWdd94pWh+GYRSITYWVclq2jGXp0mGMGjXBzVbhYQC/5x4nnP6114ZUFm9G+rFjB5OS\nkniKjWXs2GGBd6DqjMC6d4eEBCA8ccAMo6xjxnsDjh1zjPjJySFzyQ2pkT6Xf//biSawcmVe1OKS\njgOWa2Mx471RGjHjveE/lSvDHXfASy85U0ghIGRG+lw++ggmToQvv8xTKgAXXHABAFOmTGH27Nl5\n5WlpaSS4RjWGYQQXUyxlGPepqTa1M5mYPJ3y48ZBCOJWhcRIn8umTU4csHnzoHnzU3aXVBwwwzAc\nTLGUUTynppLJpGe1RZz3/GTqPzoy6P0FJZeKNw4fht69HXfpiy/2WiUnJ4eYmBjKucLkL1iw4JTI\nwhUrVuSVV14pniycapsxjLKI2VjKKN4WTnbnXV6scRszH7q32HlaDhw4xPLl2zhypDo1a2bSvn0d\n3n03hwMHJpNrY6lTZxgDB9aiTp1aRepDVLlxzhwO16xJ5oQJxMfH+6wbyjhgAB999BHz5s3j/fff\nZ//+/Vx77bVceOGFDB06NKj9GEY4sZX3BWCKBbp2TSQ5eXS+MiGHjTWa8acZLzpZJouIL0P91Kl9\neeWVj4tvpM/loYec3CpLlkClSkVvxzAMv4hK472I9AOSgFZAJ1X91ke9KUBPYI+qtnErTwTuBHKT\niz+mqh+GVOgoxdvUlHKUtzv+lRGPPw69ekH58kVq25eh/pVXihDy3hfPP++Eolm1ypSKYUQYEaVY\ngHVAX+DlQupNAyYDp8ZVh4mqOjHYgpU2fK0fuX7KPzl60w3MuOg65lRpmz8ashe8rU0JqaEe4J13\nYNQXSg4AAA9cSURBVPx4R6nUqROcNj3Iyspi5MiRNGzYkOzsbPbt28eECRMoX0RlaxhliYhSLKq6\nEUCk4Al+Vf1MRHzNoURPEvcw4mvhJMCgHWcwYccq7mcOxzhBSsrJaazCIhCnpCTSurUQEkM9OAnK\n7rzTGa20aFH89nwwatQojh8/ziOPPALAgw8+yKOPPsrTTz8dsj4No9SgqhH3AZYD7QupEwt871GW\nCKQCa4HXgFo+jlXDOwkJSQoZupBr9H6eVWc5+49ao8atChmu7QyNixuuvXo97Fameft6935A4+KG\nn1J/69a04gm3Zo1qw4aqixcH52R9cOzYMY2JidGVK1fmla1atUrr1asX0n4NI9JxvTsLfYeX+IhF\nRJYCjdyLAAUeV9XiBqx6ERijqioi/wQmArd7q5jkltwqPj6+QI+iskTuNNbjjONjLmMKt5PBW2Rk\n5HpzQa7N5MiRoXib8jp8OIalS28rOIxMoPz0E/ToAS+8AFddVfR2/OC7774jIyMjX1riFi1asH//\nftasWcN5550X0v4NI1JITk4uUobTElcsqto9hG3vddt8FfCpqJIiMGtiJJBr1F/PuSzmasbzCPdS\nF28KRCQDX1NeLVvGBs9Qv2ULXH65Y1fp14+pU6eyf/9+Nm7cmBf7q23btowePZo+ffoUu7vt27cD\nUN1toWjNmjUB2LFjhykWo8zg+aN79OjRviu7EVE2Fg8Ks5WIZx0Raayqu12b1wLrQyFYacbdqP8Q\nE1lLW66r0pj5fzyKpwLp3DmW9euDFEDSF5s2OUElH38cBg0qkaCSR48eBaBKlSp5Zbmr9I8cORK8\nczOMUkpEKRYR6YPj7VUfeE9E1qpqDxFpAryqqj1d9WYD8UA9EdkGJKrqNOBpEWkH5ABpwF1hOI2o\nxtOo/7/y3Zi1ZhHbGv2dr9L/jbsCefbZBwCCO+Xlztq1zrTXmDFOLDMgIyODa6+9lqSkpHyjk759\n++YLx1KhQoUir6SvXbv2KWUZGRlAfmVjGIZ3bIGkUTiPPsrvX65mSJOL2blLg7O4sTBWrnTyqrz4\nopPl0oPTTjuN2bNnc8kllwBOzpWEhIR88b+Kytq1a+nQoQO7d++mQYMGAKSnp9OyZUu++uorOnTo\nUOw+DCMaicoFkkaEMno01bp0YeZ1DeDee0Pf37x5MHQozJ7tTIN5EOqgkm3atKFevXqkpqbmKZYf\nfviBmJgYzj333GCcoWGUbvxxHSttH8zdOHA2btQ/6tXTuddeq4mJiSH5jB41Sj+74AI92qiR6jff\n+BRl/fr12rBhw7zt+fPn6+bNm4NymomJiYrjpZjv06ZNGx0xYoTf9RMTEwNq3+pb/Wioj5/uxjYV\nZvjP+vXOCGLiRLjppuC2vW8f3HijE/zyjTegfv0Cq4c6qGRmZiYPPvggsbGxeSvvn3nmGSpZ+Bij\nDGNBKAvAFEsxWL/ecf2dONFRBMVFFd580wkoecst8M9/FjlGmWEYocUUSwGYYikm69fDFVc4imXM\nmKInBktNdWwpO3fCyy9Dly7BldMwjKDir2IJQvAmo8xxzjmOK/Devc7/H3wQ2PHr1zvxvjp2hK5d\n4ZtvTKkYRinCRixG8Vi6FO6+G5o2dXK49O4NbqFQAMjJcUKypKQ4017r18M998Bdd4HL68owjMjH\npsIKwBRLkDl2DD75xAlnv3AhVKsGNWtC5cpQoYKjVOrXd0YlV1wB11/v7DMMI6rwV7GE3fU3HB+8\nuNMRwS5+0VS/HOifQNuAvnLHHaqffaa6Z0/UyG/1rb7V910fzN3YJzZiMQzDCBwz3huGYRhhwRSL\nYRiGEVRMsRiGYRhBxRSLYRiGEVRMsRiGYRhBxRSLYRiGEVRMsRiGYRhBxRSLYRiGEVRMsRiGYRhB\nxRSLYRiGEVRMsRiGYRhBxRSLYRiGEVQiSrGISD8RWS8iJ0SkvY86p4nIMhH5QUTWicjf3PbVEZGP\nRGSjiCwRkVolJ71hGIYBEaZYgHVAX2BFAXWygYdUtTXQBbhXRM5y7RsJfKyqfwGWAY+GUthIJTk5\nOdwihBQ7v+imNJ9faT63QIgoxaKqG1V1E+AzLLOq7lbVta7/M4CfgGau3b2B6a7/pwN9QihuxFLa\nH247v+imNJ9faT63QIgoxRIoItICaAekuIoaquoecBQQ0DA8khmGYZRdKpR0hyKyFGjkXoSTsexx\nVV0UQDs1gHnA/aqa6aOaZfMyDMMoYSIyg6SILAeGq+q3PvZXAN4DPlDVSW7lPwHxqrpHRBoDy1W1\nlZfjI++kDcMwogB/MkiW+IglAAoSfirwo7tScfEuMBgYD9wCLPR2sD8XxjAMwygaETViEZE+wGSg\nPnAQWKuqPUSkCfCqqvYUkQuBT3E8yNT1eUxVPxSRusBbwOlAOtBfVQ+G41wMwzDKKhGlWAzDMIzo\nJ6q9woqCiFwpIhtE5GcReSTc8gQTEZkiIntE5PtwyxIKClocG+2ISGUR+VJE/r+9+4+1uq7jOP58\ngQKC8qNy5rpyC805GmqQwiZkRb+mzmzLxdSx6axVNCup2bR0q7Xyr344ZjGNH2lqumjVysGCodAv\nEuSHgNVIhEqmSQkygQuv/vh+DvtyOOfcc/Vz7/dy7vuxnZ3v+X4/5/N5fz/n3vM+35+fDWnd7qo6\npv4gaZik9ZJ+WXUsuUl6TtLG9Bn+uep4cpM0TtKjkral/8HpTcsOpS0WScOAvwKzgX8B64A5trdX\nGlgmkmYC+4Glti+sOp7c0gkZb7X9dDor8CngYx30+Y22fUDScGAtcIvtjvqCkvQlYBow1vbVVceT\nk6QdwDTbe6uOpT9IWgystr0onUA12vYrjcoOtS2WS4G/2d5p+zDwMMVFlR3B9hqgI/+oodeLY096\ntg+kyZEUJ9Z01K8+SV3AFcB9VcfST0SHfqdKGgvMsr0IwHZPs6QCHdoJLbwN2FV6vZsO+mIaSkoX\nx/6p2kjySbuJNgAvACtsr6s6psy+C3yFDkuYJQZWSFon6VNVB5PZO4CXJC1KuzIXSjqtWeGhllhC\nB6i7OHZ/1fHkYvuo7XcDXcB0SZOrjikXSVcCe9IWp2h9OcHJ6jLbUym2yualXdOd4hRgKrAgreMB\ninszNjTUEss/gYml111pXjhJpH27jwE/sd3wOqWTXdrFsAr4aNWxZHQZcHU6DvEQ8H5JSyuOKSvb\n/07PLwLLKHa9d4rdwC7bf0mvH6NINA0NtcSyDjhPUrekEcAciosqO0mn/hqsaXZx7ElN0ltqwzyk\nXQwfAjripAQA27fbnmh7EsX/3Urbc6uOKxdJo9OWNJLGAB8GtlQbVT7pHoy7JJ2fZs0GtjYrP5iv\nvM/O9hFJnweWUyTV+21vqzisbCT9FHgf8GZJzwN31Q62dYJ0cez1wOZ0LOLYxbHVRpbF2cCSdObi\nMOAR27+pOKbQvrOAZel2UacAD9peXnFMud0CPCjpVGAHcGOzgkPqdOMQQgj9b6jtCgshhNDPIrGE\nEELIKhJLCCGErCKxhBBCyCoSSwghhKwisYQQQsgqEksIIYSsIrGEEELIKhJLCKEykkZWHUPILxJL\nGHCS9vWyfJykzw5UPE1iaBnjG6j3Dklb0kiD6yVd0h/ttBnLCf0sac0Atn8VcHqbZadU2VehbyKx\nhCr0dh+hCcDn+lqppJw338xyr6NyTJJmUNxS/WLbFwEf5PjxgbLrpU9O6GfbA3Kr9zQa6Bm2/9NO\nedubgVnp7tZhkIvEEiohaaKkrWnAoC2SHi/tFvk2MCn9or87lb8+jQm/XtK9KnRL2i5piaTaF88J\ndaZym0ttz5d0Z+n13NJY5Utqs0vLT2g7zV+WBnXaLOnmNK8+pq7Sap8NvGS7B8D2y7ZfaNVOqm+b\npAfSuv1M0qg+9ElXs1ib9PO+9HxrKrtJ0hdKMXa3+Nz64ibgF318zwrg2tfRVhhotuMRjwF9AK8A\n3cBhYEqa9whwXZruBjaVyl9AMbzB8PR6AXBDKncEuKT0vkP1dTaobz5wZ5qeTHF7+gnp9fhajK3a\nris7CthMsQXQDfTUYqpb7zHAhtTeAuC9ba7jUWBGmn8/cGu7fVKqv1msmxp8NlOBjansGIrbv1/U\nqo/r6pgJfAf4IcXYKzMb9MX36l7PAeYCtwM3t/jb+UHVf7/x6P0Rm5WhKgZ2uNjFAfAU8PYmZWdT\nfNmtS1sLo4A9wJPAcz5+CN9/NKhzbYs4PgA8ansvgO3/ttk2wBclXZOmu4B3pmU73WBYYduvSpoK\nzErtPizpq7aX9rKOz9v+Y6rmAYrblx8EprXZJ61irSeKQbmW2X4NQNLPU8wbU5lGfVz2IrAfWAms\ntn2wQTujjjVYjPHxEds3SroN2CHpUmAscK7tHzV4fxjEIrGEKpW/cI5Q+rKpI2CJ7TuOmyl1A6+2\nUWcPMLw0v1k7fWn7corkMN32QUmrSvXWx3SMi5/dTwBPpF1Vc4GlLdrpblDN0fS8uJ0+6SXWZlod\nm2n5udl+VtJ7gLttH25Sx6ml6RuAX6XpqcA9wDdsf1nSZEnn2K4dixrTS9xhEIhjLKEKqnuutw84\no/T6d8AnJJ0JIGmCpNoQ0/V1NKpzD3Bmet9I4KrSspXAtZLeVKu7rp5mbY8D9qYv6guAGb3EgKTz\nJZ1XmnUxsLONdZwoaXqavg5Yk+Jut0+axVrfzzVPAtdIGqViNMSPp3kt16+0ngJGtEgqUCT7mvHA\ndhUDSI2l2BI7LS3bTzGIVs2RVm2HwSG2WEIVXPd8/EL7ZUlrJW0Cfmv7NklfB5arGGHxEDCPImHU\n13FCnbZ7JH2TYmjq3cC20rKtkr4FrJbUQ3EM5KZaPba3Sfpag7YfBz4j6RngWeAPrWJITgfuUTEE\ncQ/wd+DTvbSzJ9U/T9Ii4BngXtuvtShf337DWFM//77cz8VsPy1pceovAwttbyzV19sZcxMpdpG1\ncqA0vZRiKN93pfjOAv6Xlo3n+F125feFQSpGkAxhEEu7tn5te0rVseQkaT7F0OD1x7Rqy2dS7C67\n0Pb307xzgctt/3jgIg2vR2yxhDD4deKvv/uATwILGy20XbtQc1Vp9pXNyofBJbZYQgiVSFslO0sH\n5luVnQScY3t1/0cW3qhILCGEQU/SCNuHqo4jtCcSSwghhKzidOMQQghZRWIJIYSQVSSWEEIIWUVi\nCSGEkFUklhBCCFlFYgkhhJBVJJYQQghZRWIJIYSQVSSWEEIIWUViCSGEkFUklhBCCFlFYgkhhJBV\nJJYQQghZ/R/O3P5CaW2wugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6e8de7da20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "s_m = 6\n",
    "plt.plot(s_vector,avg, 'bo', label='Original data')\n",
    "long_s = np.linspace(0,s_m,100)\n",
    "plt.plot(long_s, V(long_s, D_e, a, r_e, const), 'r-', label='Fitted data')\n",
    "plt.title(r'Morse potential for H$_2$ with vibrational energy spectrum')\n",
    "plt.ylim(-1.2,-0.85)\n",
    "plt.xlim(0,s_m)\n",
    "plt.xlabel(r'Internuclear Seperation $r$ $(a_0)$')\n",
    "plt.ylabel(r'Energy $E$ (Ha)')\n",
    "text_grote = 16\n",
    "plt.hlines(const, 0, s_m, linestyles='dashed')\n",
    "i = 0\n",
    "plt.hlines(e_v(0), find_right()+0.55, s_m, linestyles='dashed')\n",
    "plt.hlines(const+D_e, 0, s_m, linestyles='dashed')\n",
    "\n",
    "for i in range(v_max):\n",
    "    right = find_right()\n",
    "    left = find_left()\n",
    "    plt.hlines(e_v(i), left, right, linestyles='solid', alpha=0.5)\n",
    "    if i < 5:\n",
    "        plt.annotate(\n",
    "            r'$v = %s$' %i, xy=(right+0.1, V(right, D_e, a, r_e, const)-.002), xycoords = 'data',\n",
    "            xytext = (0, 0), textcoords = 'offset points', size=text_grote)\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(1.2*s_max, const+D_e), xycoords = 'data',\n",
    "    xytext = (1.2*s_max, const), textcoords = 'data',\n",
    "    arrowprops = {'arrowstyle':'<->'})\n",
    "\n",
    "plt.annotate(\n",
    "    r'$D_e$', xy=(1.2*s_max, (y_min-y_max)/2+y_max), xycoords = 'data',\n",
    "    xytext = (0, 0), textcoords = 'offset points', size=text_grote)\n",
    "\n",
    "plt.annotate(\n",
    "    '', xy=(1.1*s_max, const+D_e), xycoords = 'data',\n",
    "    xytext = (1.1*s_max, e_v(0)), textcoords = 'data',\n",
    "    arrowprops = {'arrowstyle':'<->'})\n",
    "\n",
    "plt.annotate(\n",
    "    r'$D_0$', xy=(1.1*s_max, (y_min-y_max)/2+y_max), xycoords = 'data',\n",
    "    xytext = (0, 0), textcoords = 'offset points', size=text_grote)\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}