animation-1d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IPython Notebook to animate the results provided by our solver of the time-dependent Schroedinger equation for a quantum mechanical point particle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation\n",
    "import numpy as np\n",
    "from IPython.display import display, clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set the path/ filename of the file the contains the evolution of the wave function and the file that contains the potential here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputfilePathReal = 'sample_outputs/hydr_wavefunction_output_real.txt'\n",
    "inputfilePathImag = 'sample_outputs/hydr_wavefunction_output_imag.txt'\n",
    "potentialfilePath = 'sample_outputs/hydr_potential.txt'\n",
    "observablesfilePath = 'sample_outputs/hydr_observables_output.txt'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read the data from the files and save it in lists of floats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = 1.602176634*1e-19\n",
    "hbar = 1.054571817*1e-34\n",
    "\n",
    "potData = [] #Array to save the potential data to be plotted as a reference (units do not match up with wave function)\n",
    "with open(potentialfilePath,'r') as potentialfile:\n",
    "    for line in potentialfile:\n",
    "        potData.append(line.strip())\n",
    "    a = float(potData[0].split(\"\\t\")[1]) # retrieve lattice spacing a\n",
    "    epsilon = float(potData[0].split(\"\\t\")[3]) # retrieve energy scale epsilon\n",
    "    tauhat = float(potData[0].split(\"\\t\")[5]) # retrieve time step tauhat\n",
    "    omega = float(potData[0].split(\"\\t\")[7]) # retrieve omega\n",
    "    del potData[0] # remove first element of the list\n",
    "    potData = [float(potData[i])*epsilon/e*1e3 for i in range(len(potData))] # cast list entries to float\n",
    "\n",
    "plotDataReal = [] #Array to save the real part of the wave functions\n",
    "with open(inputfilePathReal,'r') as inputfileReal:\n",
    "    for line in inputfileReal:\n",
    "        plotDataReal.append(line.strip().split('\\t')) #Tab is used as seperator between values\n",
    "    for j in range(len(plotDataReal)): #Cast list items from string to float type\n",
    "        plotDataReal[j] = [float(plotDataReal[j][i])*a**(-1/2)*1e-5 for i in range(len(plotDataReal[0]))] # cast entries to float\n",
    "    plotDataReal = np.array(plotDataReal) # cast to numpy array\n",
    "        \n",
    "plotDataImag = [] #Array to save the imaginary part of the wave functions\n",
    "with open(inputfilePathImag,'r') as inputfileImag:\n",
    "    for line in inputfileImag:\n",
    "        plotDataImag.append(line.strip().split('\\t'))\n",
    "    for j in range(len(plotDataImag)):\n",
    "        plotDataImag[j] = [float(plotDataImag[j][i])*a**(-1/2)*1e-5 for i in range(len(plotDataImag[0]))]\n",
    "    plotDataImag = np.array(plotDataImag)\n",
    "    \n",
    "obsData = [] # Array to save the observable data\n",
    "with open(observablesfilePath,'r') as obsfile:\n",
    "    for line in obsfile:\n",
    "        obsData.append(line.strip().split('\\t'))\n",
    "    for j in range(len(obsData)):\n",
    "        obsData[j] = [float(obsData[j][i]) for i in range(len(obsData[0]))]\n",
    "    obsData = np.array(obsData).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find the absolute maximum of the wave function over the whole time to use it as the y- maximum of the plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "maximum = 0.0 #Initialise a variable to store the absolute maximum of all wave functions\n",
    "\n",
    "for i in range(len(plotDataReal)): #Iterate over all time steps\n",
    "    if max(plotDataReal[i]) > maximum: #If the maximum in the current time step is bigger than the maximum up to here...\n",
    "        maximum = max(plotDataReal[i]) #...set it as the new maximum\n",
    "        \n",
    "for i in range(len(plotDataImag)): # Check for a bigger value in the imaginary part\n",
    "    if max(plotDataImag[i]) > maximum:\n",
    "        maximum = max(plotDataImag[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate a series of plots that is shown over time. Each plot contains the wave function at consecutive points in the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEMCAYAAADAqxFbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOXZ//HPlX0PISC7JmyRJRAwiBtCXakVFVGUapGq1da6Um3t8lO0i/XRVovLY7XuWsQN2/rUDQpFFIWoQVkjSyQssgTIAllnrt8fZzLMZJ1AJjNJrvfrFTJzzsmZe4ZkvnMv575FVTHGGGPqRIS6AMYYY8KLBYMxxhg/FgzGGGP8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8WPBYIwxxk9UqAtwJHr06KEZGRmhLoYxxnQon3322V5V7dnScR0yGDIyMsjLywt1MYwxpkMRkW8COc6akowxxvixYDDGGOPHgsEYY4wfCwZjjDF+LBiMMcb4sWAwxhjjx4LBGGOMn5BfxyAiccBSIBanPK+r6t3N/UxxcTHPPfec37YRI0Ywbtw4ampqePnllxv8TE5ODjk5ORw6dIhXX321wf7c3FxGjhxJSUkJCxYsaLD/5JNPJisri7179/L222832H/66aczcOBAvv32W959990G+88880wGDBhAUVERixYtarB/8uTJ9O7dm82bN7N06dIG+88//3x69OjBhg0bWL58eYP9U6dOJTU1ldWrVzd6jcf06dNJSEggPz+f/Pz8BvuvuOIKoqOjWblyJWvWrGmwf9asWQB8/PHHFBQU+O2Ljo7miiuuAOC///0vW7Zs8dufkJDA9OnTAVi4cCHbtm3z25+SksLFF18MwLvvvsu3337rtz89PZ0pU6YA8K9//Yvi4mK//b1792by5MkAvPnmm5SWlvrt79+/P2eddRYAr776KocOHfLbn5mZycSJEwF4+eWXqamp8ds/dOhQTjnlFIAGv3dgv3v2u9dxf/eaEvJgAKqAM1S1XESigWUi8o6qfhLqghljTFckqhrqMniJSAKwDPiJqn7a1HG5ublqVz4bY0zriMhnqprb0nFh0ccgIpEikg/sBj5oLhSMMcYEV1gEg6q6VDUH6A+cKCIj6x8jIteJSJ6I5O3Zs6f9C2mMMV1EWARDHVU9ACwBJjey70lVzVXV3J49W5wc0BhjzBEKeTCISE8R6ea5HQ+cBawPbamMMabrCodRSX2A50UkEieoXlXVhmPyjDHGtIuQB4OqfgmMCXU5jDHGOELelGSMMSa8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8WPBYIwxxo8FgzHGGD8WDMYYY/xYMBhjjPFjwWCMMcaPBYMxxhg/FgzGGGP8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8WPBYIwxxo8FgzHGGD8WDMZ0EAUFBfzqV79ix44dAR3/9ttv8/DDD+NyuYJcMtPZRIW6AMaYwMyaNYvly5ezadMm5s+f3+yxO3fuZMqUKQCkpaVx1VVXtUcRTSdhNQZjOoCCggKWL18OwFtvvUVFRUWzx7/xxhve2y+++GJQy2Y6n5AHg4gMEJHFIrJORNaIyC2hLpMx4Wbx4sUAPPzww1RXV7Ny5cpmj//kk0/o168fs2fPZtmyZS0GiTG+Qh4MQC3wM1UdBpwE/FREhoe4TMaElU8//ZQePXowY8YMAD777LNmj8/LyyM3N5dJkyZRVVXV4vHG+Ap5MKjqTlX93HO7DFgH9AttqYwJL5988gnjx4+nZ8+epKWlsWHDhiaPrayspKCggJycHMaMGQPAl19+2V5FNZ1AyIPBl4hkAGOAT0NbEmPCR1VVFRs2bGDMmDGICFlZWc0GQ2FhIarK4MGD6devH2lpaRYMplXCJhhEJAl4A7hVVUsb2X+diOSJSN6ePXvav4DGhMimTZtwu90cf/zxAExJTSVz1aomj9+yZQsAmZmZiAijRo1iVTPHG1NfWASDiETjhMLLqvpmY8eo6pOqmququT179mzfAhoTQnW1g6ysLADufP99ntm/n+rKykaP9w0GgOHDhzdbwzCmvpAHg4gI8DSwTlX/HOryGBNu6t7Uhw4dCkCEKgC7P/+80eO3bNlCbGwsvXv3BmDQoEHs37+f/fv3t0NpTWcQ8mAATgV+AJwhIvmer/NCXShjwsWGDRvo27cvKSkp4DPsdP9XXzV6/NbNm7kpLY2IvXsBJxjAaZIyJhAhDwZVXaaqoqqjVDXH8/XvUJfLmHCxceNGhgwZ4tzx6V87UFjY6PFjv/ySB779Fm68EbBgMK0X8mAwxjRv69atHHfccc6d8nLv9vJt2xo9fvi33zo3Fi0Ct5uBAwcCTsAYEwgLBmPCWG1tLdu3b+fYY491NvgEQ0UTk+kNPnTIubFvH2zcSGJiIr1797YagwmYBYMxYWznzp24XK5Gg8FdXNzg+IpDhxjodlPoqSWwfj0AgwcPtmAwAbNgMCaMbd26FaDRYJCSkgbH71m/nligODvb2VBQAEBGRgbffPNNUMtqOg8LBmPCWN2beWPBEOFzu86BtWudG8OHQ8+e3mAYMGAA27dvt7UZTEAsGIwJY3U1hgEDBjgbPGFQGRVFRF1fgo+Dng7mhMGDYfBgDn25kaws2Lz5u9TW1vJtXce0Mc2wYDAmjG3dupXu3buTlJTkbPAEw8HERCKrqhocX715MwCpw4bBgAGUryuioADeeusUINYbNMY0x4LBmDC2devWw81IcLjGkJJCnMvVYJ0F986dAKQPHw4DBpBSto2oSKWqKhKYQFFRUXsV3XRgFgzGhDG/axgAKishOhpNSiIB2Ou5urmOHjhAJRCbmkplzwHEaSX/74ZiRBQ42WoMJiAWDMaEsQY1hupqiIkhMjmZRBoGA2VllEc4f9bbI5x+idOOK2LoUCEqapwFgwmIBYMxYaqkpISSkpLGgyE1lUSg/hT0EQcPUhEVBcA3rv4AZEQWMWYMiIy1piQTEAsGY8JUg2sYwAmG6GhiunVrtMYQfegQlTExABRUODWG3jVFjBgBNTX9KCzc3S5lNx2bBYMxYarJYIiJITYtjQSguN7VzzGVlVTHxgKwtrgX1USTUFxE3Rx833wT1R5FNx2cBYMxYWrLliKgb6PBEJOWRiJwoN4aC3E1NdQmJABQtD2C3dH9oOhwMOzfn95gJJMx9VkwGBOmFiwYAGzn0KHehzfW1Dh9DCkpRAFl9WoMCbW1uDzXPOzeDfsSBvgFAwy1fgbTIgsGY8LUF1+MAOCrr3z+TD01Bjy1ggqfYKisrCRJFZKTAWfphtJuA2DbNpKTIS2tGhhiwWBaZMFgTJiqqIgEwK9/uS4YEhMBqPQJhv3795MCRKSmAk6NoSLdCQZnXQY3MNgm0zMtsmAwJkzV1lYC4Nda5BmVVBcM1T59DPt37SIeiExLo7oaSkqgtu8Ap/lp926GDo0GjrNrGUyLLBiMCUO1tbXU1grQRI0hPt656zP1dun27QBEp6d7VwCV/s61DBQVkZkZCfT3dGob0zQLBmPC0I4dOwCnSchvdu16fQy1paXeXeWeFd3ijjnGGwwxgzyzshYV4cysEcXGjQ1nZTXGlwWDMWHom2+2UhcMfrNre0Yl1QWDq6zMu+ugZwK9uF692O25ji0h63AwZGQ4NwsLNYglN52BBYMxYWjTpu2AcwWzXzDUqzG4y8txu93OLk+bU0KvXt5+idRBPSAuzqfGALt2xXt/xpjGWDAYE4a+/vrw1BXNBUOcKuWetqZaT0d0Qq9e1HU9dEsT6N8fioqou07O5erHrl27gv4cTMfVYjCISGkLX2UiUnA0hRCRZ0Rkt4isPprzGNNZbNlyuMf54EGfHfWCIQE4cOAAAC7P96ju3b3BkJoKDHAucouPh27dqoDjbMiqaVYgNYZNqprSzFcycLDFszTvOWDyUZ7DmE5j27bDY1Qb1Biio32CIdkbDN40SEmhpASiojyDlzzBANC/vwsLBtOSQIJhWhsd0yRVXQrsO5pzGNOZ7NjhNAvFxzfflJTAb7zBIHUd0SkpHDjg1BZEcIJhxw6orWXwYLuWwbSsxWBQ1c2+90UkUUQimzvGGHN0du503uzT0xsfleSKca5jiKeCrVudAyLKy3EDJCZSUuJpRgInGNxu+PZbbzAUFlowmKYF0scQISLfF5H/E5HdwHpgp4isEZEHRGRIS+doCyJynYjkiUhe/cVJjOlMSkpKOHTIGTXUIBg8NYY9xRFUEEcChygocI6NOnSIQ5GREBHRMBgAtm71jEyKo6DgQHs9HdMBBdKUtBgYBPwS6K2qA1T1GGAC8AnwRxG5MohlBEBVn1TVXFXN7dmzZ7AfzpiQKdq0ib9QwAhWNxkMO3bAIRJI4BC7djnXJURVVlIZHQ3gHwyDBjnfN270DlndssWGq5qmBbJqx1mqWlN/o6ruA94A3hCR6DYvmTFdVOkHH3Azm8niZ/yt+3t4l09wu6G2FmJi2LcPjiGeBA6xd68zdUZsZSVVnkV6Skpg4EDPzw0c6PREb9hAxuXOph07Ytr3SZkOJZA+hgahICKxInKpiLzV1DGtISLzgOVAlohsE5FrjuZ8xnRk5WvXApDJFrp1A5fLyQRqPH9m0dGUlR2uMezb5/wZx9XUUB0XB9SrMURHO+GwYYP36ueKimMo8ZlnyRhfAV/gJo4zReRZYAfwPZxhpkdNVWeoah9VjVbV/qr6dFuc15iOqHrTJu/tujf36uq6f4CYGL9g2L/fqbAn1NZS65lczy8YALKyYP16kpMhObkKyLAhq6ZJTQaDiFzl6XD+fyLyILAJ+CnwDlCuqrNU9a32KqgxXYXUzZJKDd26OduqqvALhvLyumAop7w8hqqqKpJVcSUm4nZDaWkjwbBxI7hc9OtXC2SyZcuW9nxapgNprsYwG7gKuAfIBSap6sWq+irgao/CGdMl7XMu6UmhlJQUZ1N1NYebknxqDIlSyqFDUZSWlpICaFIS5eWgWi8Yjj/eSZfCQs+6DBls3LixHZ+U6Uia63x+BHgeJxiSgP+IyCpgHjbHkjFB4Xa7Ec/cR8mU1S274NQY3P5NSRUk0FO2U1UVQ2lpKelAeWqq/3QYdbKzne/5+QwdOgjIYMOGr9vhGZmOqMk3eFX9m6p+T1XvUdU7gCHA48B5QLKIvCgiF7dXQY3pCnbu3EmiZ+bTGGqIj6gCGu9jqI5OIFEOUVUVS2lJiXdZz0aDYfRopxN6xQoyMwHiWLPGdwUgYw4L+JO/Ohap6tVAP+BtnKYmY0wb2bRpEyk+9xNdzkI8fn0MnlFJNdEJJFBJTU085bt2EYGzrGejwRAbCzk58Omn3pFJBQVHNZjQdGKBXPmcVH+bqlaq6nxVvbCpY4wxrbdp0yaSfe4nuJz5jxqrMbhiEojXSmprE6j0rMwTlZ7eeDAAjB8PeXlk9ncCYe/eRA4ePNr5L01nFEiN4R8i8icROV1EEus2ishAEblaRN7DZkY1pk2sX7+eFOBA7DEAJNQ2UmPwBIM7LoF4rcLtTqLKM01MbI8eTQfDGWfAwYNk7vjIsyHTOqBNowK5wO1MYBFwPbBGREpEpBh4CegLXKWqrwe3mMZ0DWvWrCFFhOLYvgDEuZyO6PqjksrLnWCI0yqERA596wRDXM+eTQfD2WdDTAxxC98mPb0GyODrr60D2jQUyJQYqOq/gX8HuSzGdHkFa9YQp8q+qGMYBMS6KwFPjSHicI2hshLvDKtxVLL/G2dFtvhevSj5wjmsQTAkJTm1htdfZ8jAP1BcnMn69Z8E/0mZDqc1Vz4vCmSbMebIHDx4kOLCQgCKI5ympFh1gqF+H0N1Nbhi69ZkOETpdufah1hPjSEy0rtkg79rroFvvmFq/AdERQ3lyy+/DOZTMh1UIJ3PcSLSHeghImki0t3zlYHTlGSMaQPr16/3djzvEScYYnxrDD6jkqqqnKYkcIKhcrfTfiSe4areRXrqu/BC6NOHi7Y9gsvVjy++WBPEZ2Q6qkBqDNcDnwHHe77Xff0DeCx4RTOma1m9erU3GHapEwzRrsZrDPWDoXaf0xdRt6xng2akOtHR8NOfMnTze4zUNWzc6Kbcc0GdMXUC6Xz+CzAY+J2qDlTVTM/XaFV9NPhFNKZrWLlyJb09lzrvcjtrjjRaY2gkGKTMOY7k5OaDAeAnP8EVn8jP+R8gi6+++ioIz8Z0ZAH1MaiqC+eKZ2NMkKxcuZIxgwcDsKPWqTFE1frUGHxGJVVVgcYfDob42lqqRSA2tuVg6N6d2quv43Je4Vh6k5+fH6ynZDqo1sx59L6ITBNptOXSGHMUqquryc/PJ9uzxNqOGqfGUNeU1FiNwTcYUoigIsoZZNhiMACxv7gNRbgjspCVK1e2/RMyHVpAw1U9ZgOJQK2IVAKCM1NGSvM/ZoxpyVdffUV1dTXH93XGcxRVNVJjwH9UUt2wIycYlIqYGFIJLBgYMID/9L6CH377Ks8v6d3mz8d0bK2ZKylZVSNUNUZVUzz3LRSMaQNLly4FYHBv5016rzsNd0SkNxh8awzuqBhqa0ESDwdDN1xUNbZ6WzNWTLyDRCo4b8sWduzY0cbPyHRkrbmOYWxzX8EspDGd3cKFC8nKyiItMhKAcpKojYojsqbhqKQqt7Ni2+Fg2Ec3aqhKSMDtDjwY0k8fwb84n5uI4KP332/7J2U6rNb0MTwOfAI8CTwFfArMBf4EPNj2RTOma6iurmbp0qWceeaZUFaGxsfjIgpXdBwR1ZVERPjXGKo0BoCIpLpgKCGVGmp9FumpW/mtOSecAPfzC3rgpuKJJ4L19EwH1JpgKAROUNVcVT0BGANsVNXvqOoZQSmdMV3AokWLKC8vZ/LkyVBejjvBmazYHRMHlZXExHgGJHlGJVW5nK7BuhpDkhygG5W4k5KaniepEaNGwScRp7ImfjDjV66kpq5z23R5rQmG41XVO+BZVVcDOW1fJGO6lt/+dgdRUW8zYcI5UFaGK8G5zK0uGKKjfZqSYmKoqnYGBkYlxUF0NKkUk0oF2tQiPU2Ij4eR2cL8njeT5XaT/6c/BekZmo6mNcGwTkT+JiKTRGSiiDwFrAtWwYzpCoqL97N8+TXU1n6P//43FsrKqI13gkFj4701Bt9gqPtgHxsnkJbGMewimQpoapGeZpxwAvyt7Fr2isAjj7T9EzQdUmuC4YfAGuAW4FZgrWebMeYI3XvvG97bX3wBlJVRE+cJhrh6TUl1NQZntU9iY4Hu3clgG+Cs3nbggLMvkD4GgIkTYef+eJZnn8UJO3fyzZIlbfK8TMfWmuGqlar6kKpOBa4BFql6pn48SiIyWUQ2iMhGEbmzLc5pTLjbWVDAzEd+w5P8CICvv8YJhljPgog+weCtMXgm0ANPMKSlkcG3AEQ3t0hPE84+2/m+ZsJfqAEKZs5EVdvi6ZkOrDXDVZeISIpnptV84FkR+fPRFkBEInEm4/suMByYISLDj/a8xoSziooK5p93HifoLn7E35hx2lY2bgTKyqiOcWoMEhsLFRUNmpLq1xgy1VnWM7pv31YHQ58+Tif02/nD+Or00zmzqIgXbr21TZ+r6Xha05SUqqqlwMXAs56RSWe1QRlOxBndtFlVq4FXgAvb4LzGhB2Xy8XChQs55eSTOXfTJqpxrkk4NfITvv0WKCujsi4Y4uOgqqph57MnGGJigO7dvefW3r1bHQwAl18OH30Eqf/zBgfi4jhl7lxuvfxyNm3adPRP2HRIrZkSI0pE+gDTgV+3YRn6AUU+97cB45v7gX35m3glfWrDHc1WgZveJ777VJreV+8scqSP18yP1X+8hudsfKoqUbf3DM2ds/7Zm3sOzU+K1cTPaUvPoWmCNnnall+Xw0f6bm/wHNT3yMbPqU3uk8bP6fODzZ0z1nWIqpokHuJmBiZmMIxVXMvj/DXyBoZVfsGuXdPRmDIqo5xgiIiPhX1Vh/sYImoa1hjS0ryP0S1rFCV5zszanklaAzJrFtx9N/zhiR489c47JJ1zDr+YP59n5r9JeeKJJMRkkkAk8dSSSDXRWuv3Sxzo5Gk2y1rH0ZpguBd4D1imqitFZCDQFgvGNvbr0uCvS0SuA64DGE0kZ+77sNGTNfNn67Ov4R9vYD/XXvuaFh7l1KA9XnNvM+3+/9fMO9mR/B9VSRx9Irey0DWVqoMx7EsawIvlP2DugD8z4OAGqqoUqsup8A2GqipiEjw1Bqk3KikWSE/3nr8mLYMDB5pZpKcJffrA7Nlw//3Qp88kRtz8KcOevo07D3xI5MGP4OBHANQSyUESqCGqmdfN3v3DVWs+sLUYDCIyA3hfVV8DXqvbrqqbgWlHUsB6tgEDfO73BxpM3KKqT+JcdU1ubq72zMtrg4c2pp1VVFB593188cgybiv/A/0yY4k/PpPuawuddRVUORThdD5HJHiCIcY/GPxqDBkZ3lOXVkQHPB1GfffcA1u2wH33AYxh8OAlTJtZwZQzDjL+1CiiUhOJio7mCE5twkmAnxgCqTEcB7wmItHAIuAdYIW23dCFlcAQEckEtgOXA99vo3MbE17i44n7n3tJuxqG3Q8PXQfyUgYpH35MCqUAHIx05qaM8gSDt49BGhmVNGIEADvpTWmpM09SoENVfcXGwvz58NBDzntHnz4A8Z4v09UEsoLbHz1TXpwHrAKuBj4Xkb+LyEwR6XU0BVDVWuBGnGaqdcCrqmoL0ZpO7fjj4dln4eSTgcxMossPMJDNAJREdkfEpympuVFJY8eyZebdXMZ8yspg3z6/bodW69u3LhRMVxZwH4OqlgELPF94hpR+F3gBOPdoCqGq/wb+fTTnMKbDyswE4AQ+A+BAZDrx8SBxsf4XuEk1pKb6j0oSYd/Nc/jwBSgthb17vacz5oi1pvPZj6quFRG3qtoEK8YcDc87eS5Ov9k+cYKBWE+NIVqprhaQRkYlcbhPYf9+2LMHevRo5/KbTqc11zE0xiZxN+ZoDRwIwPgIZ4nNYu3uBENcHKgSH13b+FxJnmDo6awCyvbtTh9D3X1jjlQgo5LmNrULOIJuLmOMn27doHt3svatB+DbmnSSk/G+88dHVFFdHd34qCQgJcW5vXatc9+CwRytQJqSfgj8DKhqZN+Mti2OMV3UwIGwbx+HIpPYcyjRLxgSIquoqUnyG5UUEQFRnr9eEejVC9Z4hmxYU5I5WoEEw0pgtap+XH+HiMxp8xIZ0xUNGgR5eeyMPpayciElBb9gcJqPDtcY6moLdY45Blatcm5bjcEcrUCC4RKg0VlUVdXGPxjTFjzXI+zQvpSVQe/eHA6GiMoWg6FXL+8Cbxx3XLuV2nRSgVzHsE9VD7VHYYzpsqZOpTI6ib+5r6asjEb6GHDe+T3BEBPj/+N11x5ERUH//u1actMJtRgMIvJkWxxjjGnGyJHMvecAL9TMYNcu/2CIpcpvoZ7q6oY1hlGjnO9RUYf7How5UoH8Cl0kIs0tyCPAd9qoPMZ0WenHRAJQWekJhrg4wKkxuFyg1dVIE01J4z3zEU+Z0o4FNp1WIMFwRwDHND7VqTEmYL6dxqmpeN/946SKCFyIy9VkH8O4cbBgAZx6avuV13ReLQaDqj7fHgUxpqvzDYaePfFrSorG07PsGa5aPxhE4KKL2qecpvM72iufjTFtxPf6g2OOwS8YYvBc7txEjcGYtmTBYEyY8J3V1He4aqxWHq4xNDEqyZi2ZMFgTJhISjp8e9gwmqwxNDYqyZi2FNDANhHpj7OAzgSgL1ABrAb+D3hH1bvgsDHmKDz9NBQWekLC8+4fY01Jpp0FMones0A/4G3gfmA3EAcMBSYDvxaRO1V1aTALakxXcPXVPnfqgkEtGEz7CqTG8CdVXd3I9tXAmyISAxzbtsUyxtRdxxDj9gmGJkYlGdOWApkSYzWAiBxTf5+IZKlqtapuDEbhjOnSPO/+0W6rMZj21ZrO5w9FZHrdHRH5GZ5lPo0xQeANBhuVZNpXa2ZVmQQ8KSKXAr2AdcCJwSiUMQbvogvRLhuVZNpXwDUGVd0JvAucDGQAL6hqeZDKZYwBiI0lylVFbN06WdbHYNpBwDUGEfkA2AmMBPoDz4jIUlW9PViFM6bL8wRDXY1BY2ItGEzQtaaP4TFVnamqBzwd0qcAJUEqlzEGGgRDbWRs3WZjgiaQ9RgEQFXf8t2uqrWq+lvfY1pLRC4VkTUi4haR3CM5hzGdWmwskT5NSdXE1G02JmgCqTEsFpGbRMTvWgURiRGRM0TkeeCqI3z81cDFgF0cZ0xj4uKIrD1cY6jCc9GbjUoyQRRIH8Nk4GpgnohkAgdwrnyOBN4HHlLV/CN5cFVdB3CEFQ5jOr/YWCKrK701hhqxGoMJvkDWY6gEHgceF5FooAdQoaoHgl04Y7q82FgifGoMlWp9DCb4ApkraabnZjmwubW1AxFZCPRuZNevVfUfrTjPdcB1AMceazNwmC4iNpaImsN9DFVqNQYTfIE0JdWt55wAHC8iClyrqnmBPICqnnWkhat3nieBJwFyc3O1Lc5pTNiLjSWi7KDVGEy7CqQp6Ye+90Xke8BbIjJaVYuDVjJjjBMMxfu8NYZKt9UYTPAF0pQ0u5HNB4D7RWStqv75SB9cRKYCjwA9gf8TkXxVPfdIz2dMpxMbC9VOH4Mboao2ErBRSSa4AmlKSm5kWwRO01JSI/sCpqoLsIn4jGlaXBxSVUVCRBWuiFiqa5wRfFZjMMEUSFPSPb73ReRsnE7gW1R1T7AKZozBSYDKSuIiq6kVZ2bVus3GBEsgTUnPeG7GA8d7fuZSCwVj2kFsLHhqDLURsRYMpl0E0pT0X8/3g8Bm4AtVtVFBxrQHTzDESjW1EVZjMO0jkKak59ujIMaYRniCIS6yihqfGoN1Pptgas3sqsaY9uYJhniposanj8GzHLQxQWHBYEw487QZJVFOjcRaMJh2YcFgTDjzBEOyllIjMVRW+m02JigsGIwJZ56qQaKWUSWx3mCwGoMJJgsGY8JZXVOSu5QanD6GyEiICnhRXmNaz4LBmHDmCYZEVylVODUGa0YywWbBYEw486RASu1+KiSBykr1Pt9iAAAY4UlEQVRrRjLBZ8FgTDjzqR4cIoGqKgsGE3wWDMaEM59gOEiiNSWZdmHBYEw48w0GtaYk0z4sGIwJZ/WCwZqSTHuwYDAmnPmkQLk7wZqSTLuwYDAmnPmkQJk70ZqSTLuwYDAmnPkGg8uakkz7sGAwJpwlJnpvlrmsKcm0j05zYX1NTQ3btm2jsm4yGRP24uLi6N+/P9HR0aEuSvhKS/PeLKm1piTTPjpNMGzbto3k5GQyMjIQkVAXx7RAVSkuLmbbtm1kZmaGujjhKz7ee3OXuwfl5RYMJvg6TVNSZWUl6enpFgodhIiQnp5uNbxW2ENPioutKckEX6cJBsBCoYOx/6/W2c0xVFX5dTsYExQhDQYReUBE1ovIlyKyQES6hbI8bWHBggWICOvXrwdgyZIlnH/++Ud93lmzZvH66683e8ySJUv4+OOPj/qxTJi5/372DBpPCc6fR3JyiMtjOr1Q1xg+AEaq6iigAPhliMtz1ObNm8dpp53GK6+80u6PbcHQSf385yx74BPv3aSkEJbFdAkhDQZVfV9Vaz13PwH6h7I8R6u8vJyPPvqIp59+2i8YSktLmTp1KsOHD+fHP/4xbrcbl8vFrFmzGDlyJNnZ2Tz00EMA5Ofnc9JJJzFq1CimTp3K/v37GzxORkYGe/fuBSAvL49JkyZRWFjIE088wUMPPUROTg4ffvghe/bsYdq0aYwbN45x48bx0Ucftc8LYdqcb/ORBYMJtnAalXQ1ML8tTnTrrbeSn5/fFqfyysnJ4eGHH272mLfeeovJkyczdOhQunfvzueffw7AihUrWLt2LccddxyTJ0/mzTffJDMzk+3bt7N69WoADhw4AMDMmTN55JFHmDhxInfddRf33HNPi48LTlj8+Mc/Jikpidtvvx2A73//+9x2222cdtppbN26lXPPPZd169YdzctgQiQh4fBtCwYTbEGvMYjIQhFZ3cjXhT7H/BqoBV5u5jzXiUieiOTt2bMn2MU+IvPmzePyyy8H4PLLL2fevHkAnHjiiQwcOJDIyEhmzJjBsmXLGDhwIJs3b+amm27i3XffJSUlhZKSEg4cOMDEiRMBuOqqq1i6dOkRl2fhwoXceOON5OTkcMEFF1BaWkpZWdnRP1HT7qzGYNpT0GsMqnpWc/tF5CrgfOBMVdVmzvMk8CRAbm5uk8cBAX3CbmvFxcX85z//YfXq1YgILpcLEeG8885rMPpGREhLS2PVqlW89957PPbYY7z66qve5qSWREVF4Xa7AZod7ul2u1m+fDnxPmPhTcfkW2OwzmcTbKEelTQZ+AVwgaoeCmVZjtbrr7/OzJkz+eabbygsLKSoqIjMzEyWLVvGihUr2LJlC263m/nz53Paaaexd+9e3G4306ZN47e//S2ff/45qamppKWl8eGHHwLw4osvemsPvjIyMvjss88AeOONN7zbk5OT/WoE55xzDo8++qj3fls3r5n2k5Jy+LbVGEywhXpU0qNAMvCBiOSLyBMhLs8RmzdvHlOnTvXbNm3aNP7+979z8sknc+eddzJy5EgyMzOZOnUq27dvZ9KkSeTk5DBr1izuu+8+AJ5//nnuuOMORo0aRX5+PnfddVeDx7r77ru55ZZbmDBhApGRkd7tU6ZMYcGCBd7O57lz55KXl8eoUaMYPnw4TzzRYV/eLq9Hj8O3u3X4Qd0m3EkzrTdhKzc3V/Py8vy2rVu3jmHDhoWoROZI2f9b4OpaJEtLrTnJHBkR+UxVc1s6LtQ1BmNMK1komGALp+Gqxphm/PvfUFIS6lKYrsCCwZgO4rvfDXUJTFdhTUnGGGP8WDAYY4zxY8FgjDHGjwVDG4qMjCQnJ4eRI0cyZcoU7/xHR8J3ory2ZDOwGmNaYsHQhuLj48nPz2f16tV0796dxx57LNRF8lNbW2vBYIxpkQVDkJx88sls377de/+BBx5g3LhxjBo1irvvvtu7/aKLLuKEE05gxIgRPPnkky2eNykpiZ/97GeMHTuWM888k7oJBZ966inGjRvH6NGjmTZtGocOOTOMzJo1i9mzZ/Od73yHyy67rMHU3MYYU1+nHK56663Q1tMC5eRAoHPzuVwuFi1axDXXXAPA+++/z9dff82KFStQVS644AKWLl3K6aefzjPPPEP37t2pqKhg3LhxTJs2jfT09CbPffDgQcaOHcuf/vQn7r33Xu655x4effRRLr74Yn70ox8B8Jvf/Iann36am266CYCCggIWLlxIZGQkc+bM8Zua2xhj6rMaQxuqqKggJyeH9PR09u3bx9lnnw04wfD+++8zZswYxo4dy/r16/n6668BmDt3LqNHj+akk06iqKjIu70pERERXHbZZQBceeWVLFu2DIDVq1czYcIEsrOzefnll1mzZo33Zy699FK/OZWMMaY5nbLGEIJZt4HDfQwlJSWcf/75PPbYY9x8882oKr/85S+5/vrr/Y5fsmQJCxcuZPny5SQkJDBp0qRmp9FuTN2U3rNmzeKtt95i9OjRPPfccyxZssR7TKKtHm+MaQWrMQRBamoqc+fO5cEHH6SmpoZzzz2XZ555hvLycgC2b9/O7t27KSkpIS0tjYSEBNavX88nn3zSwpmdNRZef/11AP7+979z2mmnAVBWVkafPn2oqanh5ZebXO+owdTcxhhTX6esMYSDMWPGMHr0aF555RV+8IMfsG7dOk4++WTA6UB+6aWXmDx5Mk888QSjRo0iKyuLk046qcXzJiYmsmbNGk444QRSU1OZP99ZDfW3v/0t48eP57jjjiM7O7vJN/8pU6ZwySWX8I9//INHHnmECRMmtN2TNsZ0CjbtdgeTlJTkrXl0Bl3l/82YcGDTbhtjjDkiFgwdTGeqLRhjwpMFgzHGGD8WDMYYY/xYMBhjjPFjwWCMMcaPBUMbSkpKCvpjnHLKKUF/jLZ04MABHn/88VAXwxjTChYMHUxbTZldW1vbJudpjsvlsmAwbcr396nuGqw5c+b43e8ovvnmG84880xGjRrFpEmT2LZtW6PHzZ8/n1GjRjFixAh+/vOfe7dXVVVx2WWXMXjwYMaPH09hYaF333333cfgwYPJysrivffea3XZLBiCYMmSJUycOJHp06czdOhQ7rzzTl5++WVOPPFEsrOz2bRpEwD/+te/GD9+PGPGjOGss85i165dAOzZs4ezzz6bsWPHcv3113Pcccd5F+2pq5UsWbKESZMmcckll3D88cdzxRVXeP8w7r33XsaNG8fIkSO57rrrvNsnTZrEr371KyZOnMjvf/97MjMzqampAaC0tJSMjAzv/TqzZs3ixz/+MRMmTGDo0KG8/fbbABQWFjJhwgTGjh3L2LFjvYG1ZMkSvvOd7/D973+f7Oxs7rzzTjZt2kROTg533HFHMF920wX4BsNDDz3E3/72Nw4ePMivf/1rPvjggxCXrnVuv/12Zs6cyZdffsldd93FL3/5ywbHFBcXc8cdd7Bo0SLWrFnDrl27WLRoEQBPP/00aWlpbNy4kdtuu41f/OIXAKxdu5ZXXnmFNWvW8O6773LDDTfgcrlaVzhVDdkX8FvgSyAfeB/oG8jPnXDCCVrf2rVrD9+55RbViRPb9uuWWxo8Zn2JiYmqqrp48WJNTU3VHTt2aGVlpfbt21fvuusuVVV9+OGH9RbPufbt26dut1tVVZ966imdPXu2qqr+9Kc/1T/84Q+qqvrOO+8ooHv27GnwGCkpKVpUVKQul0tPOukk/fDDD1VVtbi42FumK6+8Uv/5z3+qqurEiRP1Jz/5iXffrFmzdMGCBaqq+te//tX7+L6uuuoqPffcc9XlcmlBQYH269dPKyoq9ODBg1pRUaGqqgUFBVr3f7J48WJNSEjQzZs3q6rqli1bdMSIEU2+Zn7/b8a04LLLLtO4uDgdPXq03n777XrfffdpTEyMLl26VFVV3W633n777TpixAgdOXKkvvLKK6rq/F6efvrpeumll+qQIUP0F7/4hb700ks6btw4HTlypG7cuFFVVXfv3q0XX3yx5ubmam5uri5btsy7/ayzztIxY8boddddp8cee6z3b/LCCy/UsWPH6vDhw/Wvf/1rwM9l+PDhWlRU5C13cnJyg2NWrFihZ555pvf+Cy+84P0bPuecc/Tjjz9WVdWamhpNT09Xt9utf/jDH7zvH/WPA/I0gPfYUNcYHlDVUaqaA7wN3BXi8rSZcePG0adPH2JjYxk0aBDnnHMOANnZ2d4q37Zt2zj33HPJzs7mgQce8E6VvWzZMi6//HIAJk+eTFpaWqOPceKJJ9K/f38iIiLIycnxnnfx4sWMHz+e7Oxs/vOf//hNwV03ZTfAtddey7PPPgvAs88+yw9/+MNGH2f69OlEREQwZMgQBg4cyPr166mpqeFHP/oR2dnZXHrppaxdu9avXJmZmUfwqhnTvD/+8Y8MGjSI/Px8+vXrR48ePbj55pt59913+eCDD3jzzTfJz89n1apVLFy4kDvuuIOdO3cCsGrVKv7yl7/w1Vdf8eKLL1JQUMCKFSu49tpreeSRRwC45ZZbuO2221i5ciVvvPEG1157LQD33HMPZ5xxBp9//jlTp05l69at3jI988wzfPbZZ+Tl5TF37lyKi4sB528tJyenwdcLL7wAwOjRo3njjTcAWLBgAWVlZd6frTN48GDWr19PYWEhtbW1vPXWWxQVFQHOZJwDBgwAICoqitTUVIqLi/22A/Tv399v0bBAhHQSPVUt9bmbCLRNI2Go5t32ERsb670dERHhvR8REeFt37/pppuYPXs2F1xwAUuWLGl1W6nvY0RGRlJbW0tlZSU33HADeXl5DBgwgDlz5vhN5e07Bfepp55KYWEh//3vf3G5XIwcObLRx6mb2tv3/kMPPUSvXr1YtWoVbrebuLi4Rh/DmGC55ZZbEBHmzJnDnDlzUFVmz57NjBkziIyMpFevXkycOJGVK1eSkpLi/bAGNPiwtnjxYgAWLlzo9yGntLSUsrIyli1bxoIFC4CGH9bmzp3r3Ve3pkp6erp3gsumPPjgg9x4440899xznH766fTr14+oKP+35LS0NP73f/+Xyy67jIiICE455RQ2b94MNP4+ISJNbm+NUNcYEJHfi0gRcAXN1BhE5DoRyRORvLrlLDu6kpIS+vXrB8Dzzz/v3X7aaafx6quvAs4iP/v37w/4nHUh0KNHD8rLy71TdDdl5syZzJgxo8naAsBrr72G2+1m06ZNbN68maysLEpKSujTpw8RERG8+OKLTbZh2jTfJljq3uzqPlA19aZYJ5APa263m+XLl5Ofn09+fj7bt28nOTm5yfP6rqmyatUqxowZ4/0bbKnG0LdvX958802++OILfv/73wPOlP31TZkyhU8//ZTly5eTlZXFkCFDAKcmUFd7qK2tpaSkhO7du/ttB6dlom/fvi28mv6CHgwislBEVjfydSGAqv5aVQcALwM3NnUeVX1SVXNVNbdnz57BLna7mDNnDpdeeikTJkygR48e3u13330377//PmPHjuWdd96hT58+JCcnB3TObt26eZt4LrroIsaNG9fs8VdccQX79+9nxowZTR6TlZXFxIkT+e53v8sTTzxBXFwcN9xwA88//zwnnXQSBQUFTdYS0tPTOfXUUxk5cqR1Ppuj1tIHjdNPP5358+fjcrnYs2cPS5cu5cQTTwz4/Oeccw6PPvqo936+Z43gpj6sNbemyvz5870B4/s1c+ZMAPbu3Yvb7QacUURXX311o2XavXs3APv37+fxxx/3Nm9dcMEF3g+Ur7/+OmeccQYiwgUXXMArr7xCVVUVW7Zs4euvv27VawCEtvPZ9ws4DlgdyLEtdj53cJWVlVpTU6Oqqh9//LGOHj06aI/12muv6ZVXXtnk/quuukpfe+21oD1+Z/p/M+1jxowZOmLECL399tsb7Guu8/l73/ue97iJEyfqypUrG+zbs2ePTp8+XbOzs3XYsGF6/fXXq6rqrl279IwzztAxY8borbfeqn369NHKykqtrKzUyZMna3Z2tl5yySU6ceJEXbx4cUDP47XXXtPBgwfrkCFD9JprrtHKykrvPt+/+csvv1yHDRumw4YN03nz5nm3V1RU6CWXXKKDBg3ScePG6aZNm7z7fve73+nAgQN16NCh+u9//9u7nQA7n0MdBkN8bt8EvB7Iz3X2YCgoKNCcnBwdNWqU5ubm6ooVK4LyODfeeKMOGjRIN2zY0OQxFgzGtO+HtWAKNBhCvYLbH0UkC3AD3wA/DnF5wsKQIUP44osvgv44dSMxmvPcc88FvRzGhLutW7cyffp03G43MTExPPXUU6EuUlCFelTStFA+vjHGBKK9PqyFi5CPSmpLGuAwTxMe7P/LmPDUaYIhLi6O4uJie7PpIFSV4uJiv+sfjDHhIdR9DG2mf//+bNu2jc5yjUNXEBcXR//+/UNdDGNMPZ0mGKKjo20aBmOMaQOdpinJGGNM27BgMMYY48eCwRhjjB/piKN4RKQM2BDqcoSJHsDeUBciTNhrcZi9FofZa3FYlqq2OPFaR+183qCquaEuRDgQkTx7LRz2Whxmr8Vh9locJiJ5gRxnTUnGGGP8WDAYY4zx01GD4clQFyCM2GtxmL0Wh9lrcZi9FocF9Fp0yM5nY4wxwdNRawzGGGOCpEMFg4hMFpENIrJRRO4MdXlCSUSeEZHdIrI61GUJJREZICKLRWSdiKwRkVtCXaZQEZE4EVkhIqs8r8U9oS5TqIlIpIh8ISJvh7osoSQihSLylYjkBzIyqcM0JYlIJFAAnA1sA1YCM1R1bUgLFiIicjpQDrygqiNDXZ5QEZE+QB9V/VxEkoHPgIu64u+FiAiQqKrlIhINLANuUdVPWvjRTktEZgO5QIqqnh/q8oSKiBQCuaoa0PUcHanGcCKwUVU3q2o18ApwYYjLFDKquhTYF+pyhJqq7lTVzz23y4B1QL/Qlio0PKs3lnvuRnu+OsYnvyAQkf7A94C/hbosHU1HCoZ+QJHP/W100TcA0zgRyQDGAJ+GtiSh42k6yQd2Ax+oapd9LYCHgZ/jLB3c1Snwvoh8JiLXtXRwRwoGaWRbl/00ZPyJSBLwBnCrqpaGujyhoqouVc0B+gMnikiXbGYUkfOB3ar6WajLEiZOVdWxwHeBn3qaopvUkYJhGzDA535/YEeIymLCiKc9/Q3gZVV9M9TlCQeqegBYAkwOcVFC5VTgAk/b+ivAGSLyUmiLFDqqusPzfTewAKdpvkkdKRhWAkNEJFNEYoDLgX+GuEwmxDwdrk8D61T1z6EuTyiJSE8R6ea5HQ+cBawPbalCQ1V/qar9VTUD573iP6p6ZYiLFRIikugZmIGIJALnAM2OZuwwwaCqtcCNwHs4HYyvquqa0JYqdERkHrAcyBKRbSJyTajLFCKnAj/A+USY7/k6L9SFCpE+wGIR+RLng9QHqtqlh2kaAHoBy0RkFbAC+D9Vfbe5H+gww1WNMca0jw5TYzDGGNM+LBiMMcb4sWAwxhjjx4LBGGOMHwsGY4wxfiwYjAlTnmt2PhaR10Wko67PbjogCwZjwtc04H5gKzAqxGUxXYgFgzFhQESWeNYaucBn8z+B14Dz8LlSVUQeEJFvReT29i6n6RqsempM+LhCVX0XUYkGvgR64DNDqKreISIH27twpuuwGoPp1ERknIh86VndLNGzstnIesckisj/eVY+Wy0il3m2F4rI/Z5V0VaIyGDP9iki8qlnZbCFItLLsz1JRJ71rJT1pYhM82w/R0SWi8jnIvKaZybYQPwAZy2BT3HmtzGmXVgwmE5NVVfiNMn8Dvgf4CVVrT+B2GRgh6qO9qyG5zuPTKmqngg8ijO/Pzgro52kqmNwZu78uWf7/wNKVDVbVUcB/xGRHsBvgLM80x7nAbNbKrdncsBLgdeBecAVrXzqxhwxa0oyXcG9OJPKVQI3N7L/K+BBEbkfeFtVP/TZN8/n+0Oe2/2B+Z5lRWOALZ7tZ+HM5AmAqu73rAswHPjIea8nBmfyw5acDmxQ1b0i8g7wqIgkqqo1IZmgsxqD6Qq6A0lAMhBXf6eqFgAn4ATEfSJyl+/uRm4/AjyqqtnA9T7nFBouHiU4s5zmeL6Gq2ogM+FeCYzzrCewAUgFLgrg54w5ahYMpit4EqeZ52Wc4Z9+RKQvcEhVXwIeBMb67L7M53vdJ/1UYLvn9lU+x76PMzV83XnTgE+AU336JxJEZGhzhRWRWOB8YKCqZnjWFLgQa04y7cSCwXRqIjITqFXVvwN/xPkUfka9w7KBFZ61kn+N0x9RJ1ZEPgVuAW7zbJsDvCYiHwJ7fY79HZDm6cBeBXxHVfcAs4B5nnUSPgGOb6HY3wM+U9Uyn21LgbEickwgz9uYo2HrMRjTBE8zTq6q7m3p2DZ4rCXA7fWGqzZ3/BygXFUfDGa5TNdkNQZjwsM+4Ll6F7g1SkQewOmDsI5oExRWYzDGGOPHagzGGGP8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8fP/ATpAflQY9GngAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEMCAYAAADAqxFbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOXZ//HPlX0PISC7JmyRJRAwiBtCXakVFVGUapGq1da6Um3t8lO0i/XRVovLY7XuWsQN2/rUDQpFFIWoQVkjSyQssgTIAllnrt8fZzLMZJ1AJjNJrvfrFTJzzsmZe4ZkvnMv575FVTHGGGPqRIS6AMYYY8KLBYMxxhg/FgzGGGP8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8WPBYIwxxk9UqAtwJHr06KEZGRmhLoYxxnQon3322V5V7dnScR0yGDIyMsjLywt1MYwxpkMRkW8COc6akowxxvixYDDGGOPHgsEYY4wfCwZjjDF+LBiMMcb4sWAwxhjjx4LBGGOMn5BfxyAiccBSIBanPK+r6t3N/UxxcTHPPfec37YRI0Ywbtw4ampqePnllxv8TE5ODjk5ORw6dIhXX321wf7c3FxGjhxJSUkJCxYsaLD/5JNPJisri7179/L222832H/66aczcOBAvv32W959990G+88880wGDBhAUVERixYtarB/8uTJ9O7dm82bN7N06dIG+88//3x69OjBhg0bWL58eYP9U6dOJTU1ldWrVzd6jcf06dNJSEggPz+f/Pz8BvuvuOIKoqOjWblyJWvWrGmwf9asWQB8/PHHFBQU+O2Ljo7miiuuAOC///0vW7Zs8dufkJDA9OnTAVi4cCHbtm3z25+SksLFF18MwLvvvsu3337rtz89PZ0pU6YA8K9//Yvi4mK//b1792by5MkAvPnmm5SWlvrt79+/P2eddRYAr776KocOHfLbn5mZycSJEwF4+eWXqamp8ds/dOhQTjnlFIAGv3dgv3v2u9dxf/eaEvJgAKqAM1S1XESigWUi8o6qfhLqghljTFckqhrqMniJSAKwDPiJqn7a1HG5ublqVz4bY0zriMhnqprb0nFh0ccgIpEikg/sBj5oLhSMMcYEV1gEg6q6VDUH6A+cKCIj6x8jIteJSJ6I5O3Zs6f9C2mMMV1EWARDHVU9ACwBJjey70lVzVXV3J49W5wc0BhjzBEKeTCISE8R6ea5HQ+cBawPbamMMabrCodRSX2A50UkEieoXlXVhmPyjDHGtIuQB4OqfgmMCXU5jDHGOELelGSMMSa8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8WPBYIwxxo8FgzHGGD8WDMYYY/xYMBhjjPFjwWCMMcaPBYMxxhg/FgzGGGP8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8WPBYIwxxo8FgzHGGD8WDMZ0EAUFBfzqV79ix44dAR3/9ttv8/DDD+NyuYJcMtPZRIW6AMaYwMyaNYvly5ezadMm5s+f3+yxO3fuZMqUKQCkpaVx1VVXtUcRTSdhNQZjOoCCggKWL18OwFtvvUVFRUWzx7/xxhve2y+++GJQy2Y6n5AHg4gMEJHFIrJORNaIyC2hLpMx4Wbx4sUAPPzww1RXV7Ny5cpmj//kk0/o168fs2fPZtmyZS0GiTG+Qh4MQC3wM1UdBpwE/FREhoe4TMaElU8//ZQePXowY8YMAD777LNmj8/LyyM3N5dJkyZRVVXV4vHG+Ap5MKjqTlX93HO7DFgH9AttqYwJL5988gnjx4+nZ8+epKWlsWHDhiaPrayspKCggJycHMaMGQPAl19+2V5FNZ1AyIPBl4hkAGOAT0NbEmPCR1VVFRs2bGDMmDGICFlZWc0GQ2FhIarK4MGD6devH2lpaRYMplXCJhhEJAl4A7hVVUsb2X+diOSJSN6ePXvav4DGhMimTZtwu90cf/zxAExJTSVz1aomj9+yZQsAmZmZiAijRo1iVTPHG1NfWASDiETjhMLLqvpmY8eo6pOqmququT179mzfAhoTQnW1g6ysLADufP99ntm/n+rKykaP9w0GgOHDhzdbwzCmvpAHg4gI8DSwTlX/HOryGBNu6t7Uhw4dCkCEKgC7P/+80eO3bNlCbGwsvXv3BmDQoEHs37+f/fv3t0NpTWcQ8mAATgV+AJwhIvmer/NCXShjwsWGDRvo27cvKSkp4DPsdP9XXzV6/NbNm7kpLY2IvXsBJxjAaZIyJhAhDwZVXaaqoqqjVDXH8/XvUJfLmHCxceNGhgwZ4tzx6V87UFjY6PFjv/ySB779Fm68EbBgMK0X8mAwxjRv69atHHfccc6d8nLv9vJt2xo9fvi33zo3Fi0Ct5uBAwcCTsAYEwgLBmPCWG1tLdu3b+fYY491NvgEQ0UTk+kNPnTIubFvH2zcSGJiIr1797YagwmYBYMxYWznzp24XK5Gg8FdXNzg+IpDhxjodlPoqSWwfj0AgwcPtmAwAbNgMCaMbd26FaDRYJCSkgbH71m/nligODvb2VBQAEBGRgbffPNNUMtqOg8LBmPCWN2beWPBEOFzu86BtWudG8OHQ8+e3mAYMGAA27dvt7UZTEAsGIwJY3U1hgEDBjgbPGFQGRVFRF1fgo+Dng7mhMGDYfBgDn25kaws2Lz5u9TW1vJtXce0Mc2wYDAmjG3dupXu3buTlJTkbPAEw8HERCKrqhocX715MwCpw4bBgAGUryuioADeeusUINYbNMY0x4LBmDC2devWw81IcLjGkJJCnMvVYJ0F986dAKQPHw4DBpBSto2oSKWqKhKYQFFRUXsV3XRgFgzGhDG/axgAKishOhpNSiIB2Ou5urmOHjhAJRCbmkplzwHEaSX/74ZiRBQ42WoMJiAWDMaEsQY1hupqiIkhMjmZRBoGA2VllEc4f9bbI5x+idOOK2LoUCEqapwFgwmIBYMxYaqkpISSkpLGgyE1lUSg/hT0EQcPUhEVBcA3rv4AZEQWMWYMiIy1piQTEAsGY8JUg2sYwAmG6GhiunVrtMYQfegQlTExABRUODWG3jVFjBgBNTX9KCzc3S5lNx2bBYMxYarJYIiJITYtjQSguN7VzzGVlVTHxgKwtrgX1USTUFxE3Rx833wT1R5FNx2cBYMxYWrLliKgb6PBEJOWRiJwoN4aC3E1NdQmJABQtD2C3dH9oOhwMOzfn95gJJMx9VkwGBOmFiwYAGzn0KHehzfW1Dh9DCkpRAFl9WoMCbW1uDzXPOzeDfsSBvgFAwy1fgbTIgsGY8LUF1+MAOCrr3z+TD01Bjy1ggqfYKisrCRJFZKTAWfphtJuA2DbNpKTIS2tGhhiwWBaZMFgTJiqqIgEwK9/uS4YEhMBqPQJhv3795MCRKSmAk6NoSLdCQZnXQY3MNgm0zMtsmAwJkzV1lYC4Nda5BmVVBcM1T59DPt37SIeiExLo7oaSkqgtu8Ap/lp926GDo0GjrNrGUyLLBiMCUO1tbXU1grQRI0hPt656zP1dun27QBEp6d7VwCV/s61DBQVkZkZCfT3dGob0zQLBmPC0I4dOwCnSchvdu16fQy1paXeXeWeFd3ijjnGGwwxgzyzshYV4cysEcXGjQ1nZTXGlwWDMWHom2+2UhcMfrNre0Yl1QWDq6zMu+ugZwK9uF692O25ji0h63AwZGQ4NwsLNYglN52BBYMxYWjTpu2AcwWzXzDUqzG4y8txu93OLk+bU0KvXt5+idRBPSAuzqfGALt2xXt/xpjGWDAYE4a+/vrw1BXNBUOcKuWetqZaT0d0Qq9e1HU9dEsT6N8fioqou07O5erHrl27gv4cTMfVYjCISGkLX2UiUnA0hRCRZ0Rkt4isPprzGNNZbNlyuMf54EGfHfWCIQE4cOAAAC7P96ju3b3BkJoKDHAucouPh27dqoDjbMiqaVYgNYZNqprSzFcycLDFszTvOWDyUZ7DmE5j27bDY1Qb1Biio32CIdkbDN40SEmhpASiojyDlzzBANC/vwsLBtOSQIJhWhsd0yRVXQrsO5pzGNOZ7NjhNAvFxzfflJTAb7zBIHUd0SkpHDjg1BZEcIJhxw6orWXwYLuWwbSsxWBQ1c2+90UkUUQimzvGGHN0du503uzT0xsfleSKca5jiKeCrVudAyLKy3EDJCZSUuJpRgInGNxu+PZbbzAUFlowmKYF0scQISLfF5H/E5HdwHpgp4isEZEHRGRIS+doCyJynYjkiUhe/cVJjOlMSkpKOHTIGTXUIBg8NYY9xRFUEEcChygocI6NOnSIQ5GREBHRMBgAtm71jEyKo6DgQHs9HdMBBdKUtBgYBPwS6K2qA1T1GGAC8AnwRxG5MohlBEBVn1TVXFXN7dmzZ7AfzpiQKdq0ib9QwAhWNxkMO3bAIRJI4BC7djnXJURVVlIZHQ3gHwyDBjnfN270DlndssWGq5qmBbJqx1mqWlN/o6ruA94A3hCR6DYvmTFdVOkHH3Azm8niZ/yt+3t4l09wu6G2FmJi2LcPjiGeBA6xd68zdUZsZSVVnkV6Skpg4EDPzw0c6PREb9hAxuXOph07Ytr3SZkOJZA+hgahICKxInKpiLzV1DGtISLzgOVAlohsE5FrjuZ8xnRk5WvXApDJFrp1A5fLyQRqPH9m0dGUlR2uMezb5/wZx9XUUB0XB9SrMURHO+GwYYP36ueKimMo8ZlnyRhfAV/gJo4zReRZYAfwPZxhpkdNVWeoah9VjVbV/qr6dFuc15iOqHrTJu/tujf36uq6f4CYGL9g2L/fqbAn1NZS65lczy8YALKyYP16kpMhObkKyLAhq6ZJTQaDiFzl6XD+fyLyILAJ+CnwDlCuqrNU9a32KqgxXYXUzZJKDd26OduqqvALhvLyumAop7w8hqqqKpJVcSUm4nZDaWkjwbBxI7hc9OtXC2SyZcuW9nxapgNprsYwG7gKuAfIBSap6sWq+irgao/CGdMl7XMu6UmhlJQUZ1N1NYebknxqDIlSyqFDUZSWlpICaFIS5eWgWi8Yjj/eSZfCQs+6DBls3LixHZ+U6Uia63x+BHgeJxiSgP+IyCpgHjbHkjFB4Xa7Ec/cR8mU1S274NQY3P5NSRUk0FO2U1UVQ2lpKelAeWqq/3QYdbKzne/5+QwdOgjIYMOGr9vhGZmOqMk3eFX9m6p+T1XvUdU7gCHA48B5QLKIvCgiF7dXQY3pCnbu3EmiZ+bTGGqIj6gCGu9jqI5OIFEOUVUVS2lJiXdZz0aDYfRopxN6xQoyMwHiWLPGdwUgYw4L+JO/Ohap6tVAP+BtnKYmY0wb2bRpEyk+9xNdzkI8fn0MnlFJNdEJJFBJTU085bt2EYGzrGejwRAbCzk58Omn3pFJBQVHNZjQdGKBXPmcVH+bqlaq6nxVvbCpY4wxrbdp0yaSfe4nuJz5jxqrMbhiEojXSmprE6j0rMwTlZ7eeDAAjB8PeXlk9ncCYe/eRA4ePNr5L01nFEiN4R8i8icROV1EEus2ishAEblaRN7DZkY1pk2sX7+eFOBA7DEAJNQ2UmPwBIM7LoF4rcLtTqLKM01MbI8eTQfDGWfAwYNk7vjIsyHTOqBNowK5wO1MYBFwPbBGREpEpBh4CegLXKWqrwe3mMZ0DWvWrCFFhOLYvgDEuZyO6PqjksrLnWCI0yqERA596wRDXM+eTQfD2WdDTAxxC98mPb0GyODrr60D2jQUyJQYqOq/gX8HuSzGdHkFa9YQp8q+qGMYBMS6KwFPjSHicI2hshLvDKtxVLL/G2dFtvhevSj5wjmsQTAkJTm1htdfZ8jAP1BcnMn69Z8E/0mZDqc1Vz4vCmSbMebIHDx4kOLCQgCKI5ympFh1gqF+H0N1Nbhi69ZkOETpdufah1hPjSEy0rtkg79rroFvvmFq/AdERQ3lyy+/DOZTMh1UIJ3PcSLSHeghImki0t3zlYHTlGSMaQPr16/3djzvEScYYnxrDD6jkqqqnKYkcIKhcrfTfiSe4areRXrqu/BC6NOHi7Y9gsvVjy++WBPEZ2Q6qkBqDNcDnwHHe77Xff0DeCx4RTOma1m9erU3GHapEwzRrsZrDPWDoXaf0xdRt6xng2akOtHR8NOfMnTze4zUNWzc6Kbcc0GdMXUC6Xz+CzAY+J2qDlTVTM/XaFV9NPhFNKZrWLlyJb09lzrvcjtrjjRaY2gkGKTMOY7k5OaDAeAnP8EVn8jP+R8gi6+++ioIz8Z0ZAH1MaiqC+eKZ2NMkKxcuZIxgwcDsKPWqTFE1frUGHxGJVVVgcYfDob42lqqRSA2tuVg6N6d2quv43Je4Vh6k5+fH6ynZDqo1sx59L6ITBNptOXSGHMUqquryc/PJ9uzxNqOGqfGUNeU1FiNwTcYUoigIsoZZNhiMACxv7gNRbgjspCVK1e2/RMyHVpAw1U9ZgOJQK2IVAKCM1NGSvM/ZoxpyVdffUV1dTXH93XGcxRVNVJjwH9UUt2wIycYlIqYGFIJLBgYMID/9L6CH377Ks8v6d3mz8d0bK2ZKylZVSNUNUZVUzz3LRSMaQNLly4FYHBv5016rzsNd0SkNxh8awzuqBhqa0ESDwdDN1xUNbZ6WzNWTLyDRCo4b8sWduzY0cbPyHRkrbmOYWxzX8EspDGd3cKFC8nKyiItMhKAcpKojYojsqbhqKQqt7Ni2+Fg2Ec3aqhKSMDtDjwY0k8fwb84n5uI4KP332/7J2U6rNb0MTwOfAI8CTwFfArMBf4EPNj2RTOma6iurmbp0qWceeaZUFaGxsfjIgpXdBwR1ZVERPjXGKo0BoCIpLpgKCGVGmp9FumpW/mtOSecAPfzC3rgpuKJJ4L19EwH1JpgKAROUNVcVT0BGANsVNXvqOoZQSmdMV3AokWLKC8vZ/LkyVBejjvBmazYHRMHlZXExHgGJHlGJVW5nK7BuhpDkhygG5W4k5KaniepEaNGwScRp7ImfjDjV66kpq5z23R5rQmG41XVO+BZVVcDOW1fJGO6lt/+dgdRUW8zYcI5UFaGK8G5zK0uGKKjfZqSYmKoqnYGBkYlxUF0NKkUk0oF2tQiPU2Ij4eR2cL8njeT5XaT/6c/BekZmo6mNcGwTkT+JiKTRGSiiDwFrAtWwYzpCoqL97N8+TXU1n6P//43FsrKqI13gkFj4701Bt9gqPtgHxsnkJbGMewimQpoapGeZpxwAvyt7Fr2isAjj7T9EzQdUmuC4YfAGuAW4FZgrWebMeYI3XvvG97bX3wBlJVRE+cJhrh6TUl1NQZntU9iY4Hu3clgG+Cs3nbggLMvkD4GgIkTYef+eJZnn8UJO3fyzZIlbfK8TMfWmuGqlar6kKpOBa4BFql6pn48SiIyWUQ2iMhGEbmzLc5pTLjbWVDAzEd+w5P8CICvv8YJhljPgog+weCtMXgm0ANPMKSlkcG3AEQ3t0hPE84+2/m+ZsJfqAEKZs5EVdvi6ZkOrDXDVZeISIpnptV84FkR+fPRFkBEInEm4/suMByYISLDj/a8xoSziooK5p93HifoLn7E35hx2lY2bgTKyqiOcWoMEhsLFRUNmpLq1xgy1VnWM7pv31YHQ58+Tif02/nD+Or00zmzqIgXbr21TZ+r6Xha05SUqqqlwMXAs56RSWe1QRlOxBndtFlVq4FXgAvb4LzGhB2Xy8XChQs55eSTOXfTJqpxrkk4NfITvv0WKCujsi4Y4uOgqqph57MnGGJigO7dvefW3r1bHQwAl18OH30Eqf/zBgfi4jhl7lxuvfxyNm3adPRP2HRIrZkSI0pE+gDTgV+3YRn6AUU+97cB45v7gX35m3glfWrDHc1WgZveJ777VJreV+8scqSP18yP1X+8hudsfKoqUbf3DM2ds/7Zm3sOzU+K1cTPaUvPoWmCNnnall+Xw0f6bm/wHNT3yMbPqU3uk8bP6fODzZ0z1nWIqpokHuJmBiZmMIxVXMvj/DXyBoZVfsGuXdPRmDIqo5xgiIiPhX1Vh/sYImoa1hjS0ryP0S1rFCV5zszanklaAzJrFtx9N/zhiR489c47JJ1zDr+YP59n5r9JeeKJJMRkkkAk8dSSSDXRWuv3Sxzo5Gk2y1rH0ZpguBd4D1imqitFZCDQFgvGNvbr0uCvS0SuA64DGE0kZ+77sNGTNfNn67Ov4R9vYD/XXvuaFh7l1KA9XnNvM+3+/9fMO9mR/B9VSRx9Irey0DWVqoMx7EsawIvlP2DugD8z4OAGqqoUqsup8A2GqipiEjw1Bqk3KikWSE/3nr8mLYMDB5pZpKcJffrA7Nlw//3Qp88kRtz8KcOevo07D3xI5MGP4OBHANQSyUESqCGqmdfN3v3DVWs+sLUYDCIyA3hfVV8DXqvbrqqbgWlHUsB6tgEDfO73BxpM3KKqT+JcdU1ubq72zMtrg4c2pp1VVFB593188cgybiv/A/0yY4k/PpPuawuddRVUORThdD5HJHiCIcY/GPxqDBkZ3lOXVkQHPB1GfffcA1u2wH33AYxh8OAlTJtZwZQzDjL+1CiiUhOJio7mCE5twkmAnxgCqTEcB7wmItHAIuAdYIW23dCFlcAQEckEtgOXA99vo3MbE17i44n7n3tJuxqG3Q8PXQfyUgYpH35MCqUAHIx05qaM8gSDt49BGhmVNGIEADvpTWmpM09SoENVfcXGwvz58NBDzntHnz4A8Z4v09UEsoLbHz1TXpwHrAKuBj4Xkb+LyEwR6XU0BVDVWuBGnGaqdcCrqmoL0ZpO7fjj4dln4eSTgcxMossPMJDNAJREdkfEpympuVFJY8eyZebdXMZ8yspg3z6/bodW69u3LhRMVxZwH4OqlgELPF94hpR+F3gBOPdoCqGq/wb+fTTnMKbDyswE4AQ+A+BAZDrx8SBxsf4XuEk1pKb6j0oSYd/Nc/jwBSgthb17vacz5oi1pvPZj6quFRG3qtoEK8YcDc87eS5Ov9k+cYKBWE+NIVqprhaQRkYlcbhPYf9+2LMHevRo5/KbTqc11zE0xiZxN+ZoDRwIwPgIZ4nNYu3uBENcHKgSH13b+FxJnmDo6awCyvbtTh9D3X1jjlQgo5LmNrULOIJuLmOMn27doHt3svatB+DbmnSSk/G+88dHVFFdHd34qCQgJcW5vXatc9+CwRytQJqSfgj8DKhqZN+Mti2OMV3UwIGwbx+HIpPYcyjRLxgSIquoqUnyG5UUEQFRnr9eEejVC9Z4hmxYU5I5WoEEw0pgtap+XH+HiMxp8xIZ0xUNGgR5eeyMPpayciElBb9gcJqPDtcY6moLdY45Blatcm5bjcEcrUCC4RKg0VlUVdXGPxjTFjzXI+zQvpSVQe/eHA6GiMoWg6FXL+8Cbxx3XLuV2nRSgVzHsE9VD7VHYYzpsqZOpTI6ib+5r6asjEb6GHDe+T3BEBPj/+N11x5ERUH//u1actMJtRgMIvJkWxxjjGnGyJHMvecAL9TMYNcu/2CIpcpvoZ7q6oY1hlGjnO9RUYf7How5UoH8Cl0kIs0tyCPAd9qoPMZ0WenHRAJQWekJhrg4wKkxuFyg1dVIE01J4z3zEU+Z0o4FNp1WIMFwRwDHND7VqTEmYL6dxqmpeN/946SKCFyIy9VkH8O4cbBgAZx6avuV13ReLQaDqj7fHgUxpqvzDYaePfFrSorG07PsGa5aPxhE4KKL2qecpvM72iufjTFtxPf6g2OOwS8YYvBc7txEjcGYtmTBYEyY8J3V1He4aqxWHq4xNDEqyZi2ZMFgTJhISjp8e9gwmqwxNDYqyZi2FNDANhHpj7OAzgSgL1ABrAb+D3hH1bvgsDHmKDz9NBQWekLC8+4fY01Jpp0FMones0A/4G3gfmA3EAcMBSYDvxaRO1V1aTALakxXcPXVPnfqgkEtGEz7CqTG8CdVXd3I9tXAmyISAxzbtsUyxtRdxxDj9gmGJkYlGdOWApkSYzWAiBxTf5+IZKlqtapuDEbhjOnSPO/+0W6rMZj21ZrO5w9FZHrdHRH5GZ5lPo0xQeANBhuVZNpXa2ZVmQQ8KSKXAr2AdcCJwSiUMQbvogvRLhuVZNpXwDUGVd0JvAucDGQAL6hqeZDKZYwBiI0lylVFbN06WdbHYNpBwDUGEfkA2AmMBPoDz4jIUlW9PViFM6bL8wRDXY1BY2ItGEzQtaaP4TFVnamqBzwd0qcAJUEqlzEGGgRDbWRs3WZjgiaQ9RgEQFXf8t2uqrWq+lvfY1pLRC4VkTUi4haR3CM5hzGdWmwskT5NSdXE1G02JmgCqTEsFpGbRMTvWgURiRGRM0TkeeCqI3z81cDFgF0cZ0xj4uKIrD1cY6jCc9GbjUoyQRRIH8Nk4GpgnohkAgdwrnyOBN4HHlLV/CN5cFVdB3CEFQ5jOr/YWCKrK701hhqxGoMJvkDWY6gEHgceF5FooAdQoaoHgl04Y7q82FgifGoMlWp9DCb4ApkraabnZjmwubW1AxFZCPRuZNevVfUfrTjPdcB1AMceazNwmC4iNpaImsN9DFVqNQYTfIE0JdWt55wAHC8iClyrqnmBPICqnnWkhat3nieBJwFyc3O1Lc5pTNiLjSWi7KDVGEy7CqQp6Ye+90Xke8BbIjJaVYuDVjJjjBMMxfu8NYZKt9UYTPAF0pQ0u5HNB4D7RWStqv75SB9cRKYCjwA9gf8TkXxVPfdIz2dMpxMbC9VOH4Mboao2ErBRSSa4AmlKSm5kWwRO01JSI/sCpqoLsIn4jGlaXBxSVUVCRBWuiFiqa5wRfFZjMMEUSFPSPb73ReRsnE7gW1R1T7AKZozBSYDKSuIiq6kVZ2bVus3GBEsgTUnPeG7GA8d7fuZSCwVj2kFsLHhqDLURsRYMpl0E0pT0X8/3g8Bm4AtVtVFBxrQHTzDESjW1EVZjMO0jkKak59ujIMaYRniCIS6yihqfGoN1Pptgas3sqsaY9uYJhniposanj8GzHLQxQWHBYEw487QZJVFOjcRaMJh2YcFgTDjzBEOyllIjMVRW+m02JigsGIwJZ56qQaKWUSWx3mCwGoMJJgsGY8JZXVOSu5QanD6GyEiICnhRXmNaz4LBmHDmCYZEVylVODUGa0YywWbBYEw486RASu1+KiSBykr1Pt9iAAAY4UlEQVRrRjLBZ8FgTDjzqR4cIoGqKgsGE3wWDMaEM59gOEiiNSWZdmHBYEw48w0GtaYk0z4sGIwJZ/WCwZqSTHuwYDAmnPmkQLk7wZqSTLuwYDAmnPmkQJk70ZqSTLuwYDAmnPkGg8uakkz7sGAwJpwlJnpvlrmsKcm0j05zYX1NTQ3btm2jsm4yGRP24uLi6N+/P9HR0aEuSvhKS/PeLKm1piTTPjpNMGzbto3k5GQyMjIQkVAXx7RAVSkuLmbbtm1kZmaGujjhKz7ee3OXuwfl5RYMJvg6TVNSZWUl6enpFgodhIiQnp5uNbxW2ENPioutKckEX6cJBsBCoYOx/6/W2c0xVFX5dTsYExQhDQYReUBE1ovIlyKyQES6hbI8bWHBggWICOvXrwdgyZIlnH/++Ud93lmzZvH66683e8ySJUv4+OOPj/qxTJi5/372DBpPCc6fR3JyiMtjOr1Q1xg+AEaq6iigAPhliMtz1ObNm8dpp53GK6+80u6PbcHQSf385yx74BPv3aSkEJbFdAkhDQZVfV9Vaz13PwH6h7I8R6u8vJyPPvqIp59+2i8YSktLmTp1KsOHD+fHP/4xbrcbl8vFrFmzGDlyJNnZ2Tz00EMA5Ofnc9JJJzFq1CimTp3K/v37GzxORkYGe/fuBSAvL49JkyZRWFjIE088wUMPPUROTg4ffvghe/bsYdq0aYwbN45x48bx0Ucftc8LYdqcb/ORBYMJtnAalXQ1ML8tTnTrrbeSn5/fFqfyysnJ4eGHH272mLfeeovJkyczdOhQunfvzueffw7AihUrWLt2LccddxyTJ0/mzTffJDMzk+3bt7N69WoADhw4AMDMmTN55JFHmDhxInfddRf33HNPi48LTlj8+Mc/Jikpidtvvx2A73//+9x2222cdtppbN26lXPPPZd169YdzctgQiQh4fBtCwYTbEGvMYjIQhFZ3cjXhT7H/BqoBV5u5jzXiUieiOTt2bMn2MU+IvPmzePyyy8H4PLLL2fevHkAnHjiiQwcOJDIyEhmzJjBsmXLGDhwIJs3b+amm27i3XffJSUlhZKSEg4cOMDEiRMBuOqqq1i6dOkRl2fhwoXceOON5OTkcMEFF1BaWkpZWdnRP1HT7qzGYNpT0GsMqnpWc/tF5CrgfOBMVdVmzvMk8CRAbm5uk8cBAX3CbmvFxcX85z//YfXq1YgILpcLEeG8885rMPpGREhLS2PVqlW89957PPbYY7z66qve5qSWREVF4Xa7AZod7ul2u1m+fDnxPmPhTcfkW2OwzmcTbKEelTQZ+AVwgaoeCmVZjtbrr7/OzJkz+eabbygsLKSoqIjMzEyWLVvGihUr2LJlC263m/nz53Paaaexd+9e3G4306ZN47e//S2ff/45qamppKWl8eGHHwLw4osvemsPvjIyMvjss88AeOONN7zbk5OT/WoE55xzDo8++qj3fls3r5n2k5Jy+LbVGEywhXpU0qNAMvCBiOSLyBMhLs8RmzdvHlOnTvXbNm3aNP7+979z8sknc+eddzJy5EgyMzOZOnUq27dvZ9KkSeTk5DBr1izuu+8+AJ5//nnuuOMORo0aRX5+PnfddVeDx7r77ru55ZZbmDBhApGRkd7tU6ZMYcGCBd7O57lz55KXl8eoUaMYPnw4TzzRYV/eLq9Hj8O3u3X4Qd0m3EkzrTdhKzc3V/Py8vy2rVu3jmHDhoWoROZI2f9b4OpaJEtLrTnJHBkR+UxVc1s6LtQ1BmNMK1komGALp+Gqxphm/PvfUFIS6lKYrsCCwZgO4rvfDXUJTFdhTUnGGGP8WDAYY4zxY8FgjDHGjwVDG4qMjCQnJ4eRI0cyZcoU7/xHR8J3ory2ZDOwGmNaYsHQhuLj48nPz2f16tV0796dxx57LNRF8lNbW2vBYIxpkQVDkJx88sls377de/+BBx5g3LhxjBo1irvvvtu7/aKLLuKEE05gxIgRPPnkky2eNykpiZ/97GeMHTuWM888k7oJBZ966inGjRvH6NGjmTZtGocOOTOMzJo1i9mzZ/Od73yHyy67rMHU3MYYU1+nHK56663Q1tMC5eRAoHPzuVwuFi1axDXXXAPA+++/z9dff82KFStQVS644AKWLl3K6aefzjPPPEP37t2pqKhg3LhxTJs2jfT09CbPffDgQcaOHcuf/vQn7r33Xu655x4effRRLr74Yn70ox8B8Jvf/Iann36am266CYCCggIWLlxIZGQkc+bM8Zua2xhj6rMaQxuqqKggJyeH9PR09u3bx9lnnw04wfD+++8zZswYxo4dy/r16/n6668BmDt3LqNHj+akk06iqKjIu70pERERXHbZZQBceeWVLFu2DIDVq1czYcIEsrOzefnll1mzZo33Zy699FK/OZWMMaY5nbLGEIJZt4HDfQwlJSWcf/75PPbYY9x8882oKr/85S+5/vrr/Y5fsmQJCxcuZPny5SQkJDBp0qRmp9FuTN2U3rNmzeKtt95i9OjRPPfccyxZssR7TKKtHm+MaQWrMQRBamoqc+fO5cEHH6SmpoZzzz2XZ555hvLycgC2b9/O7t27KSkpIS0tjYSEBNavX88nn3zSwpmdNRZef/11AP7+979z2mmnAVBWVkafPn2oqanh5ZebXO+owdTcxhhTX6esMYSDMWPGMHr0aF555RV+8IMfsG7dOk4++WTA6UB+6aWXmDx5Mk888QSjRo0iKyuLk046qcXzJiYmsmbNGk444QRSU1OZP99ZDfW3v/0t48eP57jjjiM7O7vJN/8pU6ZwySWX8I9//INHHnmECRMmtN2TNsZ0CjbtdgeTlJTkrXl0Bl3l/82YcGDTbhtjjDkiFgwdTGeqLRhjwpMFgzHGGD8WDMYYY/xYMBhjjPFjwWCMMcaPBUMbSkpKCvpjnHLKKUF/jLZ04MABHn/88VAXwxjTChYMHUxbTZldW1vbJudpjsvlsmAwbcr396nuGqw5c+b43e8ovvnmG84880xGjRrFpEmT2LZtW6PHzZ8/n1GjRjFixAh+/vOfe7dXVVVx2WWXMXjwYMaPH09hYaF333333cfgwYPJysrivffea3XZLBiCYMmSJUycOJHp06czdOhQ7rzzTl5++WVOPPFEsrOz2bRpEwD/+te/GD9+PGPGjOGss85i165dAOzZs4ezzz6bsWPHcv3113Pcccd5F+2pq5UsWbKESZMmcckll3D88cdzxRVXeP8w7r33XsaNG8fIkSO57rrrvNsnTZrEr371KyZOnMjvf/97MjMzqampAaC0tJSMjAzv/TqzZs3ixz/+MRMmTGDo0KG8/fbbABQWFjJhwgTGjh3L2LFjvYG1ZMkSvvOd7/D973+f7Oxs7rzzTjZt2kROTg533HFHMF920wX4BsNDDz3E3/72Nw4ePMivf/1rPvjggxCXrnVuv/12Zs6cyZdffsldd93FL3/5ywbHFBcXc8cdd7Bo0SLWrFnDrl27WLRoEQBPP/00aWlpbNy4kdtuu41f/OIXAKxdu5ZXXnmFNWvW8O6773LDDTfgcrlaVzhVDdkX8FvgSyAfeB/oG8jPnXDCCVrf2rVrD9+55RbViRPb9uuWWxo8Zn2JiYmqqrp48WJNTU3VHTt2aGVlpfbt21fvuusuVVV9+OGH9RbPufbt26dut1tVVZ966imdPXu2qqr+9Kc/1T/84Q+qqvrOO+8ooHv27GnwGCkpKVpUVKQul0tPOukk/fDDD1VVtbi42FumK6+8Uv/5z3+qqurEiRP1Jz/5iXffrFmzdMGCBaqq+te//tX7+L6uuuoqPffcc9XlcmlBQYH269dPKyoq9ODBg1pRUaGqqgUFBVr3f7J48WJNSEjQzZs3q6rqli1bdMSIEU2+Zn7/b8a04LLLLtO4uDgdPXq03n777XrfffdpTEyMLl26VFVV3W633n777TpixAgdOXKkvvLKK6rq/F6efvrpeumll+qQIUP0F7/4hb700ks6btw4HTlypG7cuFFVVXfv3q0XX3yx5ubmam5uri5btsy7/ayzztIxY8boddddp8cee6z3b/LCCy/UsWPH6vDhw/Wvf/1rwM9l+PDhWlRU5C13cnJyg2NWrFihZ555pvf+Cy+84P0bPuecc/Tjjz9WVdWamhpNT09Xt9utf/jDH7zvH/WPA/I0gPfYUNcYHlDVUaqaA7wN3BXi8rSZcePG0adPH2JjYxk0aBDnnHMOANnZ2d4q37Zt2zj33HPJzs7mgQce8E6VvWzZMi6//HIAJk+eTFpaWqOPceKJJ9K/f38iIiLIycnxnnfx4sWMHz+e7Oxs/vOf//hNwV03ZTfAtddey7PPPgvAs88+yw9/+MNGH2f69OlEREQwZMgQBg4cyPr166mpqeFHP/oR2dnZXHrppaxdu9avXJmZmUfwqhnTvD/+8Y8MGjSI/Px8+vXrR48ePbj55pt59913+eCDD3jzzTfJz89n1apVLFy4kDvuuIOdO3cCsGrVKv7yl7/w1Vdf8eKLL1JQUMCKFSu49tpreeSRRwC45ZZbuO2221i5ciVvvPEG1157LQD33HMPZ5xxBp9//jlTp05l69at3jI988wzfPbZZ+Tl5TF37lyKi4sB528tJyenwdcLL7wAwOjRo3njjTcAWLBgAWVlZd6frTN48GDWr19PYWEhtbW1vPXWWxQVFQHOZJwDBgwAICoqitTUVIqLi/22A/Tv399v0bBAhHQSPVUt9bmbCLRNI2Go5t32ERsb670dERHhvR8REeFt37/pppuYPXs2F1xwAUuWLGl1W6nvY0RGRlJbW0tlZSU33HADeXl5DBgwgDlz5vhN5e07Bfepp55KYWEh//3vf3G5XIwcObLRx6mb2tv3/kMPPUSvXr1YtWoVbrebuLi4Rh/DmGC55ZZbEBHmzJnDnDlzUFVmz57NjBkziIyMpFevXkycOJGVK1eSkpLi/bAGNPiwtnjxYgAWLlzo9yGntLSUsrIyli1bxoIFC4CGH9bmzp3r3Ve3pkp6erp3gsumPPjgg9x4440899xznH766fTr14+oKP+35LS0NP73f/+Xyy67jIiICE455RQ2b94MNP4+ISJNbm+NUNcYEJHfi0gRcAXN1BhE5DoRyRORvLrlLDu6kpIS+vXrB8Dzzz/v3X7aaafx6quvAs4iP/v37w/4nHUh0KNHD8rLy71TdDdl5syZzJgxo8naAsBrr72G2+1m06ZNbN68maysLEpKSujTpw8RERG8+OKLTbZh2jTfJljq3uzqPlA19aZYJ5APa263m+XLl5Ofn09+fj7bt28nOTm5yfP6rqmyatUqxowZ4/0bbKnG0LdvX958802++OILfv/73wPOlP31TZkyhU8//ZTly5eTlZXFkCFDAKcmUFd7qK2tpaSkhO7du/ttB6dlom/fvi28mv6CHgwislBEVjfydSGAqv5aVQcALwM3NnUeVX1SVXNVNbdnz57BLna7mDNnDpdeeikTJkygR48e3u13330377//PmPHjuWdd96hT58+JCcnB3TObt26eZt4LrroIsaNG9fs8VdccQX79+9nxowZTR6TlZXFxIkT+e53v8sTTzxBXFwcN9xwA88//zwnnXQSBQUFTdYS0tPTOfXUUxk5cqR1Ppuj1tIHjdNPP5358+fjcrnYs2cPS5cu5cQTTwz4/Oeccw6PPvqo936+Z43gpj6sNbemyvz5870B4/s1c+ZMAPbu3Yvb7QacUURXX311o2XavXs3APv37+fxxx/3Nm9dcMEF3g+Ur7/+OmeccQYiwgUXXMArr7xCVVUVW7Zs4euvv27VawCEtvPZ9ws4DlgdyLEtdj53cJWVlVpTU6Oqqh9//LGOHj06aI/12muv6ZVXXtnk/quuukpfe+21oD1+Z/p/M+1jxowZOmLECL399tsb7Guu8/l73/ue97iJEyfqypUrG+zbs2ePTp8+XbOzs3XYsGF6/fXXq6rqrl279IwzztAxY8borbfeqn369NHKykqtrKzUyZMna3Z2tl5yySU6ceJEXbx4cUDP47XXXtPBgwfrkCFD9JprrtHKykrvPt+/+csvv1yHDRumw4YN03nz5nm3V1RU6CWXXKKDBg3ScePG6aZNm7z7fve73+nAgQN16NCh+u9//9u7nQA7n0MdBkN8bt8EvB7Iz3X2YCgoKNCcnBwdNWqU5ubm6ooVK4LyODfeeKMOGjRIN2zY0OQxFgzGtO+HtWAKNBhCvYLbH0UkC3AD3wA/DnF5wsKQIUP44osvgv44dSMxmvPcc88FvRzGhLutW7cyffp03G43MTExPPXUU6EuUlCFelTStFA+vjHGBKK9PqyFi5CPSmpLGuAwTxMe7P/LmPDUaYIhLi6O4uJie7PpIFSV4uJiv+sfjDHhIdR9DG2mf//+bNu2jc5yjUNXEBcXR//+/UNdDGNMPZ0mGKKjo20aBmOMaQOdpinJGGNM27BgMMYY48eCwRhjjB/piKN4RKQM2BDqcoSJHsDeUBciTNhrcZi9FofZa3FYlqq2OPFaR+183qCquaEuRDgQkTx7LRz2Whxmr8Vh9locJiJ5gRxnTUnGGGP8WDAYY4zx01GD4clQFyCM2GtxmL0Wh9lrcZi9FocF9Fp0yM5nY4wxwdNRawzGGGOCpEMFg4hMFpENIrJRRO4MdXlCSUSeEZHdIrI61GUJJREZICKLRWSdiKwRkVtCXaZQEZE4EVkhIqs8r8U9oS5TqIlIpIh8ISJvh7osoSQihSLylYjkBzIyqcM0JYlIJFAAnA1sA1YCM1R1bUgLFiIicjpQDrygqiNDXZ5QEZE+QB9V/VxEkoHPgIu64u+FiAiQqKrlIhINLANuUdVPWvjRTktEZgO5QIqqnh/q8oSKiBQCuaoa0PUcHanGcCKwUVU3q2o18ApwYYjLFDKquhTYF+pyhJqq7lTVzz23y4B1QL/Qlio0PKs3lnvuRnu+OsYnvyAQkf7A94C/hbosHU1HCoZ+QJHP/W100TcA0zgRyQDGAJ+GtiSh42k6yQd2Ax+oapd9LYCHgZ/jLB3c1Snwvoh8JiLXtXRwRwoGaWRbl/00ZPyJSBLwBnCrqpaGujyhoqouVc0B+gMnikiXbGYUkfOB3ar6WajLEiZOVdWxwHeBn3qaopvUkYJhGzDA535/YEeIymLCiKc9/Q3gZVV9M9TlCQeqegBYAkwOcVFC5VTgAk/b+ivAGSLyUmiLFDqqusPzfTewAKdpvkkdKRhWAkNEJFNEYoDLgX+GuEwmxDwdrk8D61T1z6EuTyiJSE8R6ea5HQ+cBawPbalCQ1V/qar9VTUD573iP6p6ZYiLFRIikugZmIGIJALnAM2OZuwwwaCqtcCNwHs4HYyvquqa0JYqdERkHrAcyBKRbSJyTajLFCKnAj/A+USY7/k6L9SFCpE+wGIR+RLng9QHqtqlh2kaAHoBy0RkFbAC+D9Vfbe5H+gww1WNMca0jw5TYzDGGNM+LBiMMcb4sWAwxhjjx4LBGGOMHwsGY4wxfiwYjAlTnmt2PhaR10Wko67PbjogCwZjwtc04H5gKzAqxGUxXYgFgzFhQESWeNYaucBn8z+B14Dz8LlSVUQeEJFvReT29i6n6RqsempM+LhCVX0XUYkGvgR64DNDqKreISIH27twpuuwGoPp1ERknIh86VndLNGzstnIesckisj/eVY+Wy0il3m2F4rI/Z5V0VaIyGDP9iki8qlnZbCFItLLsz1JRJ71rJT1pYhM82w/R0SWi8jnIvKaZybYQPwAZy2BT3HmtzGmXVgwmE5NVVfiNMn8Dvgf4CVVrT+B2GRgh6qO9qyG5zuPTKmqngg8ijO/Pzgro52kqmNwZu78uWf7/wNKVDVbVUcB/xGRHsBvgLM80x7nAbNbKrdncsBLgdeBecAVrXzqxhwxa0oyXcG9OJPKVQI3N7L/K+BBEbkfeFtVP/TZN8/n+0Oe2/2B+Z5lRWOALZ7tZ+HM5AmAqu73rAswHPjIea8nBmfyw5acDmxQ1b0i8g7wqIgkqqo1IZmgsxqD6Qq6A0lAMhBXf6eqFgAn4ATEfSJyl+/uRm4/AjyqqtnA9T7nFBouHiU4s5zmeL6Gq2ogM+FeCYzzrCewAUgFLgrg54w5ahYMpit4EqeZ52Wc4Z9+RKQvcEhVXwIeBMb67L7M53vdJ/1UYLvn9lU+x76PMzV83XnTgE+AU336JxJEZGhzhRWRWOB8YKCqZnjWFLgQa04y7cSCwXRqIjITqFXVvwN/xPkUfka9w7KBFZ61kn+N0x9RJ1ZEPgVuAW7zbJsDvCYiHwJ7fY79HZDm6cBeBXxHVfcAs4B5nnUSPgGOb6HY3wM+U9Uyn21LgbEickwgz9uYo2HrMRjTBE8zTq6q7m3p2DZ4rCXA7fWGqzZ3/BygXFUfDGa5TNdkNQZjwsM+4Ll6F7g1SkQewOmDsI5oExRWYzDGGOPHagzGGGP8WDAYY4zxY8FgjDHGjwWDMcYYPxYMxhhj/FgwGGOM8fP/ATpAflQY9GngAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(0,len(plotDataReal[0]))*a*1e10 # Generate x-axis here\n",
    "\n",
    "fig, ax = plt.subplots() # Generate an empty plot\n",
    "l, = ax.plot([min(x),max(x)],[-maximum*1.1,maximum*1.1]) # Initialize plot range\n",
    "\n",
    "def animate(i):\n",
    "    plt.clf()\n",
    "    plt.plot(x,np.sqrt(plotDataReal[i]**2+plotDataImag[i]**2), color=\"black\", label =\"Absolute\")\n",
    "    plt.plot(x,plotDataReal[i], color=\"blue\", label=\"Real part\")\n",
    "    plt.plot(x,plotDataImag[i], color=\"red\", label=\"Imaginary part\")\n",
    "    plt.axhline(max(np.sqrt(plotDataReal[0]**2+plotDataImag[0]**2)),linestyle = 'dashed', color = 'grey') #Add a dashed grey line as reference for the maximum of the start wave function\n",
    "    plt.xlim(min(x),max(x)) #Set the overwritten plot range again\n",
    "    plt.ylim(-maximum*1.1,maximum*1.1)\n",
    "    plt.xlabel(\"x space [Å]\")\n",
    "    plt.ylabel(\"Ψ(x) [1/sqrt(Å)]\")\n",
    "    plt.text(max(x),-maximum,\"t*omega=%02.3f\"%(obsData[0][i]),horizontalalignment='right')\n",
    "    plt.legend()\n",
    "\n",
    "for i in range(len(plotDataReal)): #Iterate over all plots\n",
    "    animate(i)\n",
    "    clear_output(wait=True)\n",
    "    display(fig)\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graphical output of the observables\n",
    "### Valid values for obs are:\n",
    "###    - \"norm-1\"  to display the time evolution of the normalisation minus 1\n",
    "###    - \"energy\"  to display the time evolution of the energy expectation value\n",
    "###    - \"x\"  to display the time evolution of the spacial coordinate expectation value\n",
    "###    - \"p\"  to display the time evoltuion of the momentum expectation value\n",
    "###    - \"delta-x\"  to display the width of the wave in real space\n",
    "###    - \"delta-p\"  to display the width of the wave in momentum space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = \"energy\" # Set the requested observable to plot here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFi9JREFUeJzt3X2UJXV95/H3R0aQSFCQkTOIyaABFHMMkitKTHwgsotkFWUlamICLh7W1dW4WZMlx+wmZh8OrpuE42GPLgIZNDqKgBFFRRw1kxh86EGUwRERJDoyh2mXFcVn4Lt/VA20bT/8Zqar+3b3+3XOPbdu3V9VfX/T3fO5Vb9bVakqJEmaz4OWugBJ0vJgYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJarJmqQtYSIccckitX79+qcuQpGVjy5Yt36qqtS1tV1RgrF+/nomJiaUuQ5KWjST/3NrWQ1KSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpyWCBkeTiJDuTbJ0y7/QkNya5L8loluUeneQTSbb1bf9gqBolSe2G3MPYAJw8bd5W4DRg8xzL3QP8x6p6PPBU4FVJjhmkQklSs8FuoFRVm5OsnzZvG0CSuZbbAezop7+bZBvwKOBLQ9UqSZrfWI9h9IHzJOAzS1uJJGlsAyPJAcDlwGur6jtztDs7yUSSicnJycUrUJJWmbEMjCQPpguLd1bVFXO1raoLqmpUVaO1a5vuYy5J2gNjFxjpBjguArZV1V8tdT2SpM6QX6vdCFwLHJ1ke5KzkrwgyXbgBOCqJFf3bQ9L8qF+0acBvwecmOT6/nHKUHVKktoM+S2pl8zy1vtmaHs7cEo//Y/A7F+jkiQtibE7JCVJGk8GhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmgwWGEkuTrIzydYp805PcmOS+5KM5lj25CQ3JflqknOGqlGS1G7IPYwNwMnT5m0FTgM2z7ZQkn2A/w08BzgGeEmSYwaqUZLUaLDAqKrNwJ3T5m2rqpvmWfR44KtVdWtV/Rh4N3DqQGVKkhqN4xjGo4BvTHm9vZ8nSVpC4xgYmWFezdo4OTvJRJKJycnJAcuSpNVtHANjO/DoKa8PB26frXFVXVBVo6oarV27dvDiJGm1GsfA+BxwZJIjkuwLvBi4colrkqRVb8iv1W4ErgWOTrI9yVlJXpBkO3ACcFWSq/u2hyX5EEBV3QP8e+BqYBtwaVXdOFSdkqQ2qZp1eGDZGY1GNTExsdRlSNKykWRLVc16XtxU43hISpI0hgwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktRkzWxvJDmtYfkfVtWHFrAeSdKYmjUwgLcB7wcyR5unAwaGJK0CcwXGh6vq38y1cJK/XeB6JEljaq4xjJfNt3BVvXQBa5EkjbG5AuObSd6W5MQkcx2WkiStAnMFxuOBCeA/A99Icl6SpyxOWZKkcTNrYFTV/62q/1NVzwKOB74GnJfkliT/fdEqlCSNhabzMKrqduAi4C3Ad4GXD1mUJGn8zBkYSR6S5PQkVwC3AL8J/Alw2GIUJ0kaH7MGRpJ3AV8HXgS8C/jFqjqjqj5cVfe2rDzJxUl2Jtk6Zd7BSa5JcnP/fNAsy/7PJDcm2ZbkzQ68S9LSmmsP42rgsVX1wqq6rKp+uAfr3wCcPG3eOcCmqjoS2NS//ilJfg14GvBE4JeBJwPP2IPtS5IWyFyD3pdU1XeTHJrkoiQfAUhyTJKzWlZeVZuBO6fNPhW4pJ++BHj+TIsCDwH2BfYDHgzc0bJNSdIwWga9N9DtbazrX38FeO1ebPPQqtoB0D8/cnqDqroW+ASwo39cXVXb9mKbkqS91BIYh1TVpcB9AFV1D9A0hrGnkvwS3XkghwOPAk5M8vRZ2p6dZCLJxOTk5JBlSdKq1hIY30vyCLrDRCR5KnDXXmzzjiTr+nWtA3bO0OYFwKer6u6quhv4MPDUmVZWVRdU1aiqRmvXrt2LsiRJc2kJjD8ErgQem+RTwNuBV+/FNq8Ezuinz6C7Iu50XweekWRNkgfTDXh7SEqSltBcV6sFoKquS/IM4Gi6S53fVFU/aVl5ko3AM4FDkmwH/gw4F7i0Hzj/OnB633YEvKKqXg5cBpwI3EC3Z/ORqvrAbvZNkrSAUlUzv5EcV1XXzblwQ5vFNBqNamJiYqnLkKRlI8mWqhq1tJ1rD+NvkjyTuW+gdBHwpN2oTZK0TM0VGA8DtjB3YPi1JElaJWYNjKpav4h1SJLGXNPVaiVJMjAkSU0MDElSk3kDI8nlSX4rieEiSatYSwi8Bfgd4OYk5yZ53MA1SZLG0LyBUVUfq6rfBY4DbgOuSfJPSV7WX7ZDkrQKzHtpEID+4oMvBX4P+DzwTuDX6a4F9cyhilssb/jAjXzp9u8sdRmStEeOOexA/uy5Txh8O/MGRn8/78cB7wCeu+teFsB7kngdDklaJVr2MM6vqo/P9Ebr9UfG3WIksyQtdy2B8fAkp02bdxdwQ1XNdC8LSdIK1BIYZwEn0N0yFboxi08DRyX5i6p6x0C1SZLGSEtg3Ac8vqruAEhyKN1XbZ8CbKYb25AkrXAt52Gs3xUWvZ3AUVV1J9B0IyVJ0vLXsofxD0k+CLy3f/2vgc1JHgp8e7DKJEljpSUwXgWcRnfeReju6X15dbfqe9aAtUmSxsicgZFkH+Dqqno2cPnilCRJGkdzjmFU1b3A95M8bJHqkSSNqZZDUj8EbkhyDfC9XTOr6jWDVSVJGjstgXFV/5AkrWLzBkZVXZJkf+AXquqmRahJkjSGWm6g9FzgeuAj/etjk1w5dGGSpPHScuLenwPH059zUVXXA0cMWJMkaQy1BMY9VXXXtHk1RDGSpPHVMui9NcnvAPskORJ4DfBPw5YlSRo3LXsYrwaeAPwI2Ah8B3jtkEVJksZPyz29v19Vr6+qJ1fVqJ/+4XzLJbk4yc4kW6fMOzjJNUlu7p8PmmXZX0jy0STbknwpyfrd6ZQkaeG1fEvqqCQX9P+Bf3zXo2HdG4CTp807B9hUVUcCm/rXM3k78KaqejzdgLs3apKkJdYyhvFe4K3AhcC9rSuuqs0z7BmcSncDJoBLgE8C/2lqgyTHAGuq6pp+PXe3blOSNJyWwLinqt6yQNs7tKp2AFTVjiSPnKHNUcC3k1xB9/XdjwHn9Ne1kiQtkZZB7w8keWWSdf0YxMFJDh6wpjXAbwCvA54MPAY4c7bGSc5OMpFkYnJycsCyJGl1awmMM4A/ovsq7Zb+MbGH27sjyTqA/nmmsYntwOer6taqugf4O+C42VZYVRf0g/GjtWvX7mFZkqT5tFxLaiHP6r6SLoDO7Z/fP0ObzwEHJVlbVZPAiex5QEmSFsisexhJ/njK9OnT3vsf8604yUbgWuDoJNuTnEUXFCcluRk4qX9NklGSC+H+e3C8DtiU5Aa6u/y9bXc7JklaWOnutDrDG8l1VXXc9OmZXo+L0WhUExPujEhSqyRbqmrU0nauMYzMMj3Ta0nSCjdXYNQs0zO9liStcHMNev9Kku/Q7U3s30/Tv37I4JVJksbKrIFRVfssZiGSpPHWch6GJEkGhiSpjYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpiYEhSWpiYEiSmhgYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaDBoYSS5OsjPJ1inzDk5yTZKb++eD5lj+wCTfTHL+kHVKkuY39B7GBuDkafPOATZV1ZHApv71bP4r8PfDlCZJ2h2DBkZVbQbunDb7VOCSfvoS4PkzLZvkV4FDgY8OVqAkqdlSjGEcWlU7APrnR05vkORBwF8CfzTfypKcnWQiycTk5OSCFytJ6ozroPcrgQ9V1Tfma1hVF1TVqKpGa9euXYTSJGl1WrME27wjybqq2pFkHbBzhjYnAL+R5JXAAcC+Se6uqrnGOyRJA1qKwLgSOAM4t39+//QGVfW7u6aTnAmMDAtJWlpDf612I3AtcHSS7UnOoguKk5LcDJzUvybJKMmFQ9YjSdpzqaqlrmHBjEajmpiYWOoyJGnZSLKlqkYtbcd10FuSNGYMDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNTEwJElNDAxJUhMDQ5LUxMCQJDUxMCRJTQwMSVITA0OS1MTAkCQ1MTAkSU0MDElSEwNDktTEwJAkNRksMJJcnGRnkq1T5h2c5JokN/fPB82w3LFJrk1yY5IvJnnRUDVKktoNuYexATh52rxzgE1VdSSwqX893feB36+qJ/TLn5fk4QPWKUlqMFhgVNVm4M5ps08FLumnLwGeP8NyX6mqm/vp24GdwNqh6pQktVnsMYxDq2oHQP/8yLkaJzke2Be4ZY42ZyeZSDIxOTm5oMVKkh4wtoPeSdYB7wBeVlX3zdauqi6oqlFVjdaudUdEkoay2IFxRx8EuwJh50yNkhwIXAX8aVV9ehHrkyTNYrED40rgjH76DOD90xsk2Rd4H/D2qnrvItYmSZrDkF+r3QhcCxydZHuSs4BzgZOS3Ayc1L8mySjJhf2ivw08HTgzyfX949ih6pQktUlVLXUNC2Y0GtXExMRSlyFJy0aSLVU1amk7toPekqTxYmBIkpoYGJKkJgaGJKmJgSFJamJgSJKaGBiSpCYGhiSpyYo6cS/JJPDPe7j4IcC3FrCcpWRfxpN9GU+rvS+/WFVNV25dUYGxN5JMtJ7tOO7sy3iyL+PJvrTzkJQkqYmBIUlqYmA84IKlLmAB2ZfxZF/Gk31p5BiGJKmJexiSpCYrMjCSnJzkpiRfTXLODO/vl+Q9/fufSbJ+ynt/0s+/Kcm/bF3nUBa6L0keneQTSbYluTHJHyzXvkx5b58kn0/yweF7cf82h/gde3iSy5J8uf/5nLCM+/If+t+vrUk2JnnIOPclySP6v4u7k5w/bZlfTXJDv8ybk2Q59iXJzyW5qv/9ujHJubtdVFWtqAewD3AL8BhgX+ALwDHT2rwSeGs//WLgPf30MX37/YAj+vXs07LOZdSXdcBxfZufB76yXPsyZbk/BN4FfHC5/o71710CvLyf3hd4+HLsC/Ao4GvA/n27S4Ezx7wvDwV+HXgFcP60ZT4LnAAE+DDwnOXYF+DngGdN+f36h93ty0rcwzge+GpV3VpVPwbeDZw6rc2pdH+cAJcBv9l/ajgVeHdV/aiqvgZ8tV9fyzqXRV+qakdVXQdQVd8FttH9gS+7vgAkORz4LeBCFs+C9yXJgXS3Jr4IoKp+XFXfXo596dutAfZPsobuP6rbB+4H7EVfqup7VfWPwA+nNk6yDjiwqq6t7n/atwPPH7QXnQXvS1V9v6o+0U//GLgOOHx3ilqJgfEo4BtTXm/nZ/9DvL9NVd0D3AU8Yo5lW9Y5hCH6cr9+F/ZJwGcWsObZDNWX84A/Bu5b+JJnNURfHgNMAn/TH167MMlDhyl/5jqn1TNjm5a+VNU3gf8FfB3YAdxVVR8dpPpZ6pxaz2xtpvVlrnVun2edQxiiL/dL8nDgucCm3SlqJQbGTMcXp38VbLY2uzt/aEP0pVsoOQC4HHhtVX1njytst+B9SfKvgJ1VtWVvi9tNQ/xc1gDHAW+pqicB3wMWY6xsiJ/LQXSffo8ADgMemuSle1Vlm73py96scwhD9KVbqNvr2wi8uapu3Z2iVmJgbAcePeX14fzs7vD9bfp/vIcBd86xbMs6hzBEX0jyYLqweGdVXTFI5T9riL48DXhektvodtlPTPK3QxQ/W53T6pmxzW78jm2vql17e5fRBcjQhujLs4GvVdVkVf0EuAL4tUGqn6XOafXM2GZaX+Za59TDNsvhb38+FwA3V9V5u13V0IM3i/2g+6R2K92nm12DRU+Y1uZV/PRg0aX99BP46UG8W+kGn+Zd5zLqS+iOw5633H8u05Z9Jos36D1IX+gGIY/up/8ceNNy7AvwFOBGurGL0B1nf/U492XK+2fys4PenwOeygOD3qcs4778N7oPiw/ao7qG7vhSPIBT6L79cwvw+n7eXwDP66cfAryXbpDus8Bjpiz7+n65m5jyDYKZ1rkc+0L37YkCvghc3z8G/wMY6ucy5f1nskiBMeDv2LHARP+z+TvgoGXclzcAXwa2Au8A9lsGfbmN7hP63XSf3o/p54/6ftwCnE9/wvNy6wvdXkrRfdFl19/+y3enJs/0liQ1WYljGJKkARgYkqQmBoYkqYmBIUlqYmBIkpoYGFr2ktyb5Popj3P6+Z9Msuj3ak6yIckLZ5h/ZpLDpry+MMkxC7zt9Ul+kOT6KfNe01/99p2zLPPY/t/t7oWsRSvPmqUuQFoAP6iqY5e6iAZn0n2f/3aAqnr5QNu5Zdq/xyvpzpH42kyNq+oW4FgDQ/NxD0OrQpK3JJno7wPwhinzb0vyxiSf7R+/1M8/vb+XwxeSbO7n7ZPkTUk+l+SLSf5tPz9Jzk/ypSRXAY+cYfsvpDsB7J39p/n9p+4B9fcueGOSLUk+luT4/v1bkzxvru3P0++30l3Y8Mp096h4xpQ9sc8n+fm9/sfV6rEYZyz68DHkA7iXB85cvR54UT//k8Conz64f96nn//E/vVtPHAW7e/Tny0O3EB35VXo70sBnA38aT+9H91Z2UcApwHX9Os+DPg28MIZ6ry/nhnqKx44G/99wEeBBwO/Alw/1/anbWM9sHXavNuAQ/rpDwBP66cPANZMaXf3Uv8sfYz3w0NSWglaDkn9dpKz6Q7DrqO7VMIX+/c2Tnn+6376U8CGJJfSXTwP4F8AT5wyPvEw4Ei6+1hsrKp7gduTfHwP+vBj4CP99A3Aj6rqJ0luoAuBubY/46GmWXwK+Kt+POOKqto+3wLSLgaGVrwkRwCvA55cVf8vyQa66/DsUtOnq+oVSZ5Cd3Om65McS3fxuVdX1dXT1n8Ke3/J659U1a513Af8qK/jvv5KpMy2/d1RVef2h81OAT6d5NlV9eW9KVyrh2MYWg0OpLu/xF1JDgWeM+39F015vha6bw5V1Weq6r8A36K7jPTVwL/rLw9PkqP6mxxtBl7cjzGsA541Sx3fpbst7p6abfvN+n7dUFVvpDuk9bi9qEerjHsYWgn2n/o1UuAjVXX/zYeq6gtJPk93ye1b6Q7LTLVfks/QfYB6ST/vTUmOpPtUv4nu8tJfpDs8dF2S0N0h7/l0Yw4n0h1K+grw97PUuQF4a5If0N0jenddOMv2d8drkzyLbtznS3SX65aaeLVarWr9zZdGVfWtpa5lIaS77e4Hq+qX92DZu6vqgAUvSiuGh6SkleVe4GHT9rjmtOvEPeCO4crSSuAehiSpiXsYkqQmBoYkqYmBIUlqYmBIkpoYGJKkJgaGJKnJ/wdAnU/QVa38nQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if(obs == \"norm-1\"):\n",
    "    plt.plot(obsData[0]*tauhat/omega*1e15,obsData[1])\n",
    "    plt.xlabel(\"Elapsed time [fs]\")\n",
    "    plt.ylabel(\"Absolute deviation from 1\")\n",
    "    plt.show()\n",
    "    \n",
    "elif(obs == \"energy\"):\n",
    "    plt.plot(obsData[0]*tauhat/omega*1e15,obsData[2]*epsilon/e)\n",
    "    plt.xlabel(\"Elapsed time [fs]\")\n",
    "    plt.ylabel(\"Energy [eV]\")\n",
    "    plt.show()\n",
    "    \n",
    "elif(obs == \"x\"):\n",
    "    plt.plot(obsData[0]*tauhat/omega*1e15,obsData[3]*a*1e10)\n",
    "    plt.xlabel(\"Elapsed time [fs]\")\n",
    "    plt.ylabel(\"x expectation value [Å]\")\n",
    "    plt.show()\n",
    "    \n",
    "elif(obs == \"p\"):\n",
    "    plt.plot(obsData[0]*tauhat/omega*1e15,obsData[4]*hbar/a)\n",
    "    plt.xlabel(\"Elapsed time [fs]\")\n",
    "    plt.ylabel(\"p expectation value [kg*m/s]\")\n",
    "    plt.show()\n",
    "    \n",
    "elif(obs == \"delta-x\"):\n",
    "    plt.plot(obsData[0]*tauhat/omega*1e15,obsData[5]*a*1e10)\n",
    "    plt.xlabel(\"Elapsed time [fs]\")\n",
    "    plt.ylabel(\"x width of the wave function [Å]\")\n",
    "    plt.show()\n",
    "    \n",
    "elif(obs == \"delta-p\"):\n",
    "    plt.plot(obsData[0]*tauhat/omega*1e15,obsData[6]*hbar/a)\n",
    "    plt.xlabel(\"Elapsed time [fs]\")\n",
    "    plt.ylabel(\"momentum width of the wave function [kg*m/s]\")\n",
    "    plt.show()\n",
    "    \n",
    "else:\n",
    "    print(\"Invalid value!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}