{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.sparse as sp\n",
      "import scipy.sparse.linalg as alg\n",
      "import math\n",
      "from time import time\n",
      "\n",
      "import pylab as pl\n",
      "\n",
      "%pylab inline --no-import-all"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 44
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# x - bezwymiarowa siatka polozen\n",
      "\n",
      "# Stale 'rho' zadeklerowane dla ponizszych typow potencjalow to bezwymiarowe zmienne\n",
      "# z bezwymiarowego rownania Schrodingera.\n",
      "\n",
      "\n",
      "def rho_Morse(m, V_c, C, h):\n",
      "    # 1/C - stala o wymiarze dlugosci\n",
      "    # m - masa czastki\n",
      "    # V_c - charakterystyczna energia\n",
      "    # h - stala Diraca\n",
      "    return 2.0 * m * V_c / (C * C * h * h)\n",
      "def rho_Morse2():\n",
      "    return 25.0\n",
      "def Morse(x):\n",
      "    # potencjal Morse'a: v(x) = V_c * ( exp(-2*c*x) - 2*exp(-c*x) )   \n",
      "    # bezwymiarowy potencjal Morse'a: V(x) = exp(-2*x) - 2*exp(-x)   \n",
      "    return np.exp(-2.0 * x) - 2.0 * np.exp(-x)\n",
      "    #return (1.0 - np.exp(-x)) ** 2\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(-3.0, 9.0, 300)\n",
      "pl.ylim(-1.5, 1.5)\n",
      "pl.plot(x, np.exp(-2.0 * x) - 2.0 * np.exp(-x))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 46,
       "text": [
        "[<matplotlib.lines.Line2D at 0xa342c8c>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD9CAYAAABQvqc9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHHxJREFUeJzt3X1UVHX+B/D3xIwVqDz4MOAMRgY6oAgYRe35kVM6ehTl\naNqKW0cWzeVY6va429lz2uCcIk3b9KzVsTZLN1codxXqIEfUHY+rIXkwNWUTXbEBhHx+QuPB+/vj\nJoo8zdw7cPnOfb/OuYeZ4c69HxXffv3M936vQZIkCUREpAt3aV0AERH1HIY+EZGOMPSJiHSEoU9E\npCMMfSIiHWHoExHpiKrQnzt3LsxmM2JjY9v9vtPpRGBgIBISEpCQkIA333xTzemIiEglo5o3Z2Rk\nYNGiRZgzZ06H+4wdOxYFBQVqTkNERF6iaqSfnJyM4ODgTvfhtV9ERL2HqpF+VwwGA/bs2YO4uDhY\nLBYsX74cMTEx7e5HRESe83Rg3a0f5I4ZMwYulwsHDhzAokWLMG3atA73lSRJ2O2NN97o9nM8/riE\n7dvFrF3033vWz/p766ZEt4Z+v3794O/vDwCYNGkSGhsbce7cue48pc/y9wfq67WugohE162hX1dX\n1/KvUWlpKSRJQkhISHee0mcx9InIG1T19GfPno2dO3fizJkzCA8PR3Z2NhobGwEAmZmZ2LhxIz78\n8EMYjUb4+/sjNzfXK0X3Nna7vdvP4e8PXL3q/eP2RO3difVri/WLxyApbQx5swiDQXF/Si+efx6I\niZG/EhEByrKTV+QKgu0dIvIGhr4gGPpE5A0MfUF0V0+fiPSFoS+IgACO9IlIPYa+INjeISJvYOgL\ngqFPRN7A0BcEe/pE5A0MfUFwpE9E3sDQFwQ/yCUib2DoC4IjfSLyBoa+INjTJyJvYOgLgiN9IvIG\nhr4g2NMnIm9g6AuCI30i8gaGviBMJuDGDeCX2xUQESnC0BeEwSC3ePhhLhGpwdAXSN++wJUrWldB\nRCJj6AukXz/g8mWtqyAikTH0BcKRPhGpxdAXCEOfiNRi6AuE7R0iUouhLxCO9IlILYa+QBj6RKQW\nQ18g/fox9IlIHYa+QPr2ZU+fiNRh6AuE7R0iUouhLxCGPhGpxdAXCKdsEpFaqkJ/7ty5MJvNiI2N\n7XCfxYsXIyoqCnFxcdi/f7+a0+keR/pEpJaq0M/IyEBRUVGH3y8sLMSxY8dQUVGBjz76CAsWLFBz\nOt1j6BORWqpCPzk5GcHBwR1+v6CgAOnp6QCApKQkXLhwAXV1dWpOqWts7xCRWsbuPHh1dTXCw8Nb\nnlutVlRVVcFsNrfZNysrq+Wx3W6H3W7vztKExJE+kb45nU44nU5Vx+jW0AcASZJaPTcYDO3ud3vo\nU/sY+kT6dueAODs72+NjdOvsHYvFApfL1fK8qqoKFoulO0/p03hFLhGp1a2hn5qainXr1gEASkpK\nEBQU1G5rh9zDK3KJSC1V7Z3Zs2dj586dOHPmDMLDw5GdnY3GX+7cnZmZicmTJ6OwsBCRkZEICAjA\np59+6pWi9apPH0CSgIYG+TERkacM0p1Ndy2KMBja9P6pfcHBwPHjQEiI1pUQkdaUZCevyBUMp20S\nkRoMfcFwBg8RqcHQFwxDn4jUYOgLhu0dIlKDoS8YjvSJSA2GvmAY+kSkBkNfMLwql4jUYOgLhlfl\nEpEaDH3BsL1DRGow9AXD0CciNRj6guGUTSJSg6EvGI70iUgNhr5gGPpEpAZDXzBs7xCRGgx9wXCk\nT0RqMPQFw9AnIjUY+oLhFblEpAZDXzC8IpeI1ODtEgXT1ATccw/Q2AgYDFpXQ0Ra4u0SdcBolEOf\nLR4iUoKhL6DAQODCBa2rICIRMfQFFBQEXLyodRVEJCKGvoA40icipRj6AuJIn4iUYugLKDCQoU9E\nyjD0BcT2DhEpxdAXENs7RKQUQ19AbO8QkVKqQ7+oqAg2mw1RUVFYunRpm+87nU4EBgYiISEBCQkJ\nePPNN9WeUveCgtjeISJljGre3NzcjIULF2Lbtm2wWCx46KGHkJqaiujo6Fb7jR07FgUFBaoKpVs4\n0icipVSN9EtLSxEZGYmIiAiYTCakpaUhPz+/zX5cV8e7GPpEpJSqkX51dTXCw8NbnlutVuzdu7fV\nPgaDAXv27EFcXBwsFguWL1+OmJiYNsfKyspqeWy322G329WU5tPY3iHSJ6fTCafTqeoYqkLf4MYy\nj2PGjIHL5YK/vz+2bNmCadOm4ejRo232uz30qXMc6RPp050D4uzsbI+Poaq9Y7FY4HK5Wp67XC5Y\nrdZW+/Tr1w/+/v4AgEmTJqGxsRHnzp1Tc1rd4zx9IlJKVegnJiaioqIClZWVaGhoQF5eHlJTU1vt\nU1dX19LTLy0thSRJCAkJUXNa3eM8fSJSSlV7x2g0YtWqVZg4cSKam5sxb948REdHY/Xq1QCAzMxM\nbNy4ER9++CGMRiP8/f2Rm5vrlcL1rG9f4Pp1+YYqRlV/gkSkN7xzlqBCQoCKCmDAAK0rISKt8M5Z\nOsIPc4lICYa+oBj6RKQEQ19QnKtPREow9AXFkT4RKcHQFxRH+kSkBENfUBzpE5ESDH1BMfSJSAmG\nvqDY3iEiJRj6guJIn4iUYOgLiuvvEJESDH1BcaVNIlKCoS8otneISAmGvqD4QS4RKcHQF1RwMHD+\nvNZVEJFouLSyoJqagHvuAX7+GfDz07oaItICl1bWEaMR6N+fLR4i8gxDX2ADBgBnzmhdBRGJhKEv\nsIEDgbNnta6CiETC0BfYgAEMfSLyDENfYAMHsr1DRJ5h6AuMI30i8hRDX2AMfSLyFENfYGzvEJGn\nGPoC40ifiDzF0BcY5+kTkacY+gLjPH0i8hRDX2Bs7xCRp7jgmsAaGoCAAPmrwaB1NUTU0zRZcK2o\nqAg2mw1RUVFYunRpu/ssXrwYUVFRiIuLw/79+9Wekn7Rpw9w773ApUtaV0JEolAV+s3NzVi4cCGK\niopw5MgRbNiwAeXl5a32KSwsxLFjx1BRUYGPPvoICxYsUFUwtcYWDxF5QlXol5aWIjIyEhERETCZ\nTEhLS0N+fn6rfQoKCpCeng4ASEpKwoULF1BXV6fmtHQbztUnIk8Y1by5uroa4eHhLc+tViv27t3b\n5T5VVVUwm82t9svKymp5bLfbYbfb1ZSmGxzpU28lScCNG0Bzs7zdfOzOazdu3DpGe1tH31PyHneP\nd/uv685fZ2fP7Xb5hkfe4HQ64XQ6VR1DVegb3Pz08M4PGtp73+2hT+7jXH266fp14OJF+TOeq1eB\n+np5u3bt1uM7n1+7Jm8NDUBjo/z19sedvXbz650BfvOxJAF33SVvfn7ydvNxZ6/dfM/NmDAY2t86\n+p6S97h7vJvujLDOnj/0kPdC/84BcXZ2tsfHUBX6FosFLper5bnL5YLVau10n6qqKlgsFjWnpdtw\nrr5vkSTgyhXg9Gl5++mn1o/Pnr0V7Bcvtt4AIDBQ3gICAH9/ebv33vYfBwYCYWFyIPXpI28mU9vH\nnb1mMt0K7zvD/M6gpN5BVegnJiaioqIClZWVGDJkCPLy8rBhw4ZW+6SmpmLVqlVIS0tDSUkJgoKC\n2rR2SDm2d8RSXw8cPw6cPAm4XEBVlfz15uPqajkoBw0CBg+Wv97cQkOBkSNvBfudm7dGk+TbVIW+\n0WjEqlWrMHHiRDQ3N2PevHmIjo7G6tWrAQCZmZmYPHkyCgsLERkZiYCAAHz66adeKZxkAwcChw5p\nXQXdTpLkED9wACgvByoq5O3YMfkf6PvvByIiAKsVCA8Hnnji1mOLBejbV+tfAfkyXpwluLw84J//\nBL74QutK9Km5GThyBCgtlUP+wAHg4EF51D16tDwyj4q6tVmtcuuDyBuUZKeqkT5pj+2dnnX1KrB7\nN7Bnj7zt3QuYzUBSEhAfD0ydCsTFya0Zot6IoS+4QYPkD/moe9y4IY/et26Vt717gYQE4P/+D1i0\nCFi/Xv4zIBIFQ19wZjPAa928q7ERcDqBjRuB/Hz5Q9IJE4AXXpDnXPfrp3WFRMox9AU3cCBw7hzQ\n1AQY+aepmCQB//kPsHYtsHkzEBkJPPmk/FpkpNbVEXkPY0JwRiMQHCxfoBUaqnU14qmqAj77TN7u\nvhv47W+B/fvlmTREvoih7wNCQ+UWD0PffSUlwIoVcp8+LQ3YsAFITOTFROT7eBMVH2A2A7W1WlfR\n+924AWzaBDzyCPD008CjjwKVlcAHH8iXyjPwSQ840vcB/DC3c5IEFBQAWVnyMgGvvy5PreR8edIj\nhr4PuNneobZ27ABefVW+iCo7G0hN5Yie9I2h7wPY3mmrshJ45RVg3z5g2TJgxgx5lE+kd/xr4APY\n3rnl+nXgjTeABx+Ur4wtLweeeoqBT3QTR/o+IDSUI30A+PZbID0dsNmA777jtEui9jD0fUBYGHDq\nlNZVaOfnn+V+/SefACtXArNmsW9P1BGGvg8YMkReh12PysvlkH/gAXmNHF6rQNQ5djp9QEiI3Mu+\nelXrSnrW558Djz0GLF4M/OtfDHwid3Ck7wMMBnm0f+qUPtaJaWiQV7j897+B7dvldeuJyD0c6fsI\ni0UfLZ4zZwCHQ/7get8+Bj6Rpxj6PmLIEKCmRusquteRI/LNSh59VF5OoX9/rSsiEg/bOz7C10f6\nu3fLSx2/8448LZOIlGHo+whfHukXFspB//nnwMSJWldDJDa2d3yEr07b3LABmDsX+OorBj6RN3Ck\n7yOsVsDl0roK7/rHP+T1c7ZtA0aN0roaIt/A0PcRQ4f6Vujn5gIvvwwUFzPwibzJIEmSpHkRBgN6\nQRlCa2wEAgLkC7RMJq2rUWfjRnkePgOfqHNKspM9fR9hMsmrbYr+Ye6OHcBzzwFbtjDwiboDQ9+H\nDB0K/Pij1lUoV1Ym36/2yy+B+HitqyHyTQx9HyJy6B8/DkyZAqxeDYwdq3U1RL6Loe9DRA39Cxfk\nwH/9dWD6dK2rIfJtimfvnDt3DrNmzcLJkycRERGBL774AkFBQW32i4iIQP/+/eHn5weTyYTS0lJV\nBVPHhg4Fvv9e6yo809QkL408YQKwYIHW1RD5PsUj/SVLlsDhcODo0aMYN24clixZ0u5+BoMBTqcT\n+/fvZ+B3s/vuk+8NK5IXXpBXCX33Xa0rIdIHxaFfUFCA9F8WQUlPT8fmzZs73JfTMXvGsGHA//6n\ndRXu++QTebZOXh5g5BUjRD1C8V+1uro6mM1mAIDZbEZdB3fmNhgMGD9+PPz8/JCZmYn58+e3u19W\nVlbLY7vdDrvdrrQ03br/fuDkSaC5GfDz07qazpWVAa+9BuzaBQQGal0NkRicTiecTqeqY3R6cZbD\n4UBtO3fcfuutt5Ceno7z58+3vBYSEoJz58612ffUqVMICwvD6dOn4XA48Ne//hXJycmti+DFWV5j\nsQDffCP393ur8+eBxETg7beBX/9a62qIxKUkOzsd6RcXF3f4PbPZjNraWoSGhuLUqVMYPHhwu/uF\nhYUBAAYNGoTp06ejtLS0TeiT99xs8fTW0L9xQ14xc8oUBj6RFhT39FNTU7F27VoAwNq1azFt2rQ2\n+9TX1+Py5csAgKtXr2Lr1q2IjY1VekpywwMP9O6+/jvvyHe/WrZM60qI9Elx6L/22msoLi7G8OHD\nsWPHDrz22msAgJqaGqSkpAAAamtrkZycjPj4eCQlJWHKlCmYMGGCdyqndg0bJl/o1Bvt3g2sWAF8\n8QXQp4/W1RDpExdc8zHr1wMFBfKMmN7k0iV5aYUVK4DUVK2rIfINXHCNMGIE8MMPWlfR1qJF8g3N\nGfhE2uLsaB8zYgRQUSF/YHpXL/knPS8PKCmRp2kSkbZ6SSyQt/TrBwQH954bqrhc8ih//Xp5vX8i\n0hZD3weNGAH8979aV3FreuaLL8rz8olIewx9H2Sz9Y7QX70aqK8H/vAHrSshopvY0/dBNhtw+LC2\nNfz4I/DnPwM7d/b+JSGI9IQjfR80apS2SyxLEvC738kraMbEaFcHEbXF0PdBsbHAwYNy+Gph3Tqg\nro5tHaLeiKHvgwYOlGfxnDzZ8+c+dQp49VVgzRr5Zu1E1Lsw9H3U6NHyaL8nSRLw/PNyaychoWfP\nTUTuYej7KC1Cf+NGedbQ66/37HmJyH0MfR+VkNCzV8CeOQMsXizfDevuu3vuvETkGS645qOOHwfG\njgWqqnrmfM88AwweDPzlLz1zPiLqhpuokLiGDQOuXQNqaoAhQ7r3XF9/Ld+t69Ch7j0PEanH9o6P\nMhiAhx8Gvv22e89z8SKwYAHwt78B/v7dey4iUo+h78Meflhe3bI7vfoqkJICPP54956HiLyD7R0f\nlpwMvPFG9x1/+3agqEjbq3+JyDP8INeHXb0qf7h6+rT3Wy+XL8vTQj/4AJg0ybvHJiL38M5Z1EpA\nABAXJ3/I6m1//KPc0mHgE4mF7R0fZ7cDO3YA48Z575g7dgBffcXZOkQi4kjfx02aBBQWeu94V64A\nzz4rr5UfFOS94xJRz2BP38c1NQFms7wkg8Wi/njPPy/fGOXTT9Ufi4jUYU+f2jAagYkT5XaMWtu3\nAwUFwHvvqT8WEWmDoa8DTz0F5OaqO8ZPP8n3u12zhm0dIpGxvaMDP/8sL8Xw3XdAeLjn779xQ74A\nKyEByMnxfn1EpAzbO9Suu+8GZs4EPvtM2fvffRe4dAnIzvZqWUSkAY70deLQIbm3X1kJ9Onj/vt2\n7QJmzJDX8Lnvvm4rj4gU4EifOhQbC4wcCfz97+6/59gx+fOAzz9n4BP5Co70daSkRG7zHD3a9bIM\nP/0kr93z0ktAZmbP1EdEnunRkf6XX36JkSNHws/PD2Wd3KKpqKgINpsNUVFRWLp0qdLTkRc88oh8\nY5U//anz/aqr5f1mz2bgE/kaxaEfGxuLTZs24bHHHutwn+bmZixcuBBFRUU4cuQINmzYgPLycqWn\nJC9YtQrYvLnjNs/Bg/IIPyMDyMrq0dKIqAcoXnvHZrN1uU9paSkiIyMREREBAEhLS0N+fj6io6OV\nnpZUCg6Wl2VwOIAffgBefBEYMEC+reL77wMffwysXAk8/bTWlRJRd+jWBdeqq6sRftvEcKvVir17\n97a7b9Ztw0q73Q673d6dpelaTAywb5+8UmZEBHDXXYCfHzBrFnDggHeWayAi73M6nXA6naqO0Wno\nOxwO1NbWtnk9JycHU6dO7fLgBoPB7UKy2EvoUWFhwLp1QGOjvDZ+cLB8i0Ui6r3uHBBnK7h4ptPQ\nLy4u9viAt7NYLHC5XC3PXS4XrFarqmOSd5lMQEiI1lUQUU/xyjz9jqYMJSYmoqKiApWVlWhoaEBe\nXh5SU1O9cUoiIlJAcehv2rQJ4eHhKCkpQUpKCib9cgulmpoapKSkAACMRiNWrVqFiRMnIiYmBrNm\nzeKHuEREGuLFWUREguIyDERE1CmGPhGRjjD0iYh0hKFPRKQjDH0iIh1h6BMR6QhDn4hIRxj6REQ6\nwtAnItIRhj4RkY4w9ImIdIShT0SkIwx9IiIdYegTEekIQ5+ISEcY+kREOsLQJyLSEYY+EZGOMPSJ\niHSEoU9EpCMMfSIiHWHoExHpCEOfiEhHGPpERDrC0Cci0hGGPhGRjjD0iYh0hKHvBU6nU+sSFBO5\ndoD1a431i0dx6H/55ZcYOXIk/Pz8UFZW1uF+ERERGD16NBISEvDwww8rPV2vJvIPjsi1A6xfa6xf\nPEalb4yNjcWmTZuQmZnZ6X4GgwFOpxMhISFKT0VERF6iOPRtNpvb+0qSpPQ0RETkRQZJZSI//vjj\nePfddzFmzJh2vz9s2DAEBgbCz88PmZmZmD9/ftsiDAY1JRAR6ZanEd7pSN/hcKC2trbN6zk5OZg6\ndapbJ9i9ezfCwsJw+vRpOBwO2Gw2JCcnt9qH/xMgIuoZnYZ+cXGx6hOEhYUBAAYNGoTp06ejtLS0\nTegTEVHP8MqUzY5G6vX19bh8+TIA4OrVq9i6dStiY2O9cUoiIlJAcehv2rQJ4eHhKCkpQUpKCiZN\nmgQAqKmpQUpKCgCgtrYWycnJiI+PR1JSEqZMmYIJEyZ4p3IiIvKc1MssX75cMhgM0tmzZ7UuxSOv\nvPKKZLPZpNGjR0vTp0+XLly4oHVJXdqyZYs0YsQIKTIyUlqyZInW5Xjkxx9/lOx2uxQTEyONHDlS\nWrlypdYlKdLU1CTFx8dLU6ZM0boUj50/f16aMWOGZLPZpOjoaOmbb77RuiSP5OTkSDExMdKoUaOk\n2bNnS9evX9e6pE5lZGRIgwcPlkaNGtXy2tmzZ6Xx48dLUVFRksPhkM6fP9/lcXrVFbkulwvFxcW4\n7777tC7FYxMmTMDhw4dx4MABDB8+HG+//bbWJXWqubkZCxcuRFFREY4cOYINGzagvLxc67LcZjKZ\n8N577+Hw4cMoKSnB+++/L1T9N61cuRIxMTFCzmD7/e9/j8mTJ6O8vBwHDx5EdHS01iW5rbKyEh9/\n/DHKyspw6NAhNDc3Izc3V+uyOpWRkYGioqJWry1ZsgQOhwNHjx7FuHHjsGTJki6P06tC/6WXXsI7\n77yjdRmKOBwO3HWX/NuZlJSEqqoqjSvqXGlpKSIjIxEREQGTyYS0tDTk5+drXZbbQkNDER8fDwDo\n27cvoqOjUVNTo3FVnqmqqkJhYSGeffZZ4WawXbx4Ebt27cLcuXMBAEajEYGBgRpX5b7+/fvDZDKh\nvr4eTU1NqK+vh8Vi0bqsTiUnJyM4OLjVawUFBUhPTwcApKenY/PmzV0ep9eEfn5+PqxWK0aPHq11\nKaqtWbMGkydP1rqMTlVXVyM8PLzludVqRXV1tYYVKVdZWYn9+/cjKSlJ61I88uKLL2LZsmUtgwWR\nnDhxAoMGDUJGRgbGjBmD+fPno76+Xuuy3BYSEoKXX34ZQ4cOxZAhQxAUFITx48drXZbH6urqYDab\nAQBmsxl1dXVdvqdHf9ocDgdiY2PbbAUFBXj77beRnZ3dsm9vHPl0VP9XX33Vss9bb72FPn364De/\n+Y2GlXZNxHZCe65cuYKZM2di5cqV6Nu3r9bluO3rr7/G4MGDkZCQ0Ct/1rvS1NSEsrIyPPfccygr\nK0NAQIBbrYXe4vjx41ixYgUqKytRU1ODK1euYP369VqXpYrBYHDr77XiZRiU6Gje//fff48TJ04g\nLi4OgPzf3gcffBClpaUYPHhwT5bYqa6uW/jss89QWFiI7du391BFylksFrhcrpbnLpcLVqtVw4o8\n19jYiBkzZuCZZ57BtGnTtC7HI3v27EFBQQEKCwtx/fp1XLp0CXPmzMG6deu0Ls0tVqsVVqsVDz30\nEABg5syZQoX+vn378Ktf/QoDBgwAADz55JPYs2cPnn76aY0r84zZbEZtbS1CQ0Nx6tQpt/KyV/y/\nctSoUairq8OJEydw4sQJWK1WlJWV9arA70pRURGWLVuG/Px83HPPPVqX06XExERUVFSgsrISDQ0N\nyMvLQ2pqqtZluU2SJMybNw8xMTF44YUXtC7HYzk5OXC5XDhx4gRyc3PxxBNPCBP4gPyZSnh4OI4e\nPQoA2LZtG0aOHKlxVe6z2WwoKSnBtWvXIEkStm3bhpiYGK3L8lhqairWrl0LAFi7dq17g59uml2k\nyv333y/clM3IyEhp6NChUnx8vBQfHy8tWLBA65K6VFhYKA0fPlx64IEHpJycHK3L8ciuXbskg8Eg\nxcXFtfyeb9myReuyFHE6ndLUqVO1LsNj3333nZSYmCjUNOXbLV26tGXK5pw5c6SGhgatS+pUWlqa\nFBYWJplMJslqtUpr1qyRzp49K40bN86jKZuqF1wjIiJx9Ir2DhER9QyGPhGRjjD0iYh0hKFPRKQj\nDH0iIh1h6BMR6cj/A2+88d+FbafAAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xa3423cc>"
       ]
      }
     ],
     "prompt_number": 46
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# rozmiar macierzy pieciodiagonalnej:\n",
      "N = 1500\n",
      "    \n",
      "# konce bezwymiarowego przedzialu obustronnie domknietego [x_min,x_max]:\n",
      "x_min = -3.0\n",
      "x_max = 9.0\n",
      "\n",
      "# ilosc wartosci wlasnych i wektorow wlasnych, ktore chcemy wyznaczyc:\n",
      "N_eigvals = 3\n",
      "\n",
      "# wybor potencjalu:\n",
      "#potential = Konwent\n",
      "potential = Morse\n",
      "# opcjonalne argumenty wybranego potencjalu:\n",
      "args_potential = ()\n",
      "\n",
      "# wybor parametru rho (odpowiedni dla danego potencjalu):\n",
      "#rho = rho_Konwent\n",
      "rho = rho_Morse2\n",
      "# opcjonalne argumety wejsciowe funkcji zwracajacej wartosc bezwymiarowej zmiennej rho\n",
      "args_rho = ()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 47
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def EigenvalEigenvecProblem5(N, x_min, x_max, N_eigvals, potential, *args_potential):\n",
      "    \n",
      "    \"\"\"\n",
      "    funkcja 'EigenvalEigenvecProblem5' rozwiazuje algebraiczne zagadnienie wlasne \n",
      "    wykorzystujac pieciopunktowa metode centralnych roznic skonczonych do przyblizenia \n",
      "    drugiej pochodnej.\n",
      "    \"\"\"\n",
      "    \n",
      "    # dyskretna siatka polozen, gdzie N - ilosc punktow siatki:\n",
      "    x = np.linspace(x_min, x_max, N)\n",
      "    \n",
      "    # krok dyskretnej siatki:\n",
      "    dx = (x_max - x_min) / (N - 1)\n",
      "    \n",
      "    # elementy na glownej diagonali macierzy pentadiagonalnej:\n",
      "    diagonal = 30.0 * np.ones(N) + 12.0 * dx * dx * rho(*args_rho) * potential(x, *args_potential)\n",
      "    # elementy na pierwszych diagonalach macierzy pentadiagonalnej:\n",
      "    diag_edge1 = -16.0 * np.ones(N - 1)\n",
      "    # elementy na drugich diagonalach macierzy pentadiagonalnej:\n",
      "    diag_edge2 = 1.0 * np.ones(N - 2)\n",
      "    # tablica niezerowych element\u00f3w macierzy pentadiagonalnej:\n",
      "    pentadiagonal = np.array([diag_edge2, diag_edge1, diagonal, diag_edge1, diag_edge2])\n",
      "    \n",
      "    # numery diagonali, gdzie 0 oznacza diagonale glowna, \n",
      "    # -1 -pierwsza diagonale dolna, 1 -pierwsza diagonale gorna,\n",
      "    # -2 -druga diagonale dolna, 2 - druga diagonale gorna:\n",
      "    offsets = [-2, -1, 0, 1, 2]\n",
      "    # macierz trojdiagonalna (hamiltonian) zapisana w postaci macierzy rzadkiej:\n",
      "    hamilton = sp.diags(pentadiagonal, offsets)\n",
      "        \n",
      "    start = time()\n",
      "    \n",
      "    # wyznaczenie N_eigvals wartosci wlasnych i odpowiadajacych im wektorow wlasnych\n",
      "    eigvals, eigvecs = alg.eigsh(hamilton, k = N_eigvals, which = \"SM\", tol = 1.0E-4)\n",
      "    # argument \"SM\" oznacza, \u017ce wyznaczone zostanie N_eigvals najmniejszych wartosci\n",
      "    # wlasnych i odpowiadajacych im wektorow wlasnych, a argument \"tol\" oznacza wzgledna \n",
      "    # dokladnosc, z jaka wyznaczone zostana wartosci wlasne\n",
      "    \n",
      "    # !UWAGA! Argument 'eigvecs' stanowi macierz o 'N' wierszach i 'k' kolumnach,\n",
      "    # wiec i-ty wektor wlasny ma 'i' kolumn o dlugosci 'N'.\n",
      "    \n",
      "    #print \"'alg.eigsh(...)' time: %.2fs\" % (time() - start)\n",
      "    \n",
      "    eigvals = eigvals / (12.0 * dx * dx * rho(*args_rho))\n",
      "    \n",
      "    return eigvals, eigvecs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# DYSKUSJA BLEDOW\n",
      "\n",
      "eigvals_a, eigvecs_a = EigenvalEigenvecProblem5(N, x_min, x_max, N_eigvals, potential, *args_potential)\n",
      "\n",
      "print eigvals_a\n",
      "\n",
      "N_list = range(100, 800, 200)\n",
      "measures = 5\n",
      "\n",
      "# blad wzgledny:\n",
      "time_list = []\n",
      "digits_list = []\n",
      "\n",
      "for N in N_list:\n",
      "    eigvals_num, eigvecs_num = EigenvalEigenvecProblem5(N, x_min, x_max, N_eigvals, potential, *args_potential)\n",
      "    \n",
      "    time_start = time()\n",
      "    for m in xrange(measures):\n",
      "        EigenvalEigenvecProblem5(N, x_min, x_max, N_eigvals, potential, *args_potential)\n",
      "    exec_time = (time() - time_start) / measures\n",
      "    \n",
      "    time_list.append( exec_time )\n",
      "    digits = [-math.log10( abs( (e_num - e_a) / e_a ) )\n",
      "                  for e_num, e_a in zip(eigvals_num, eigvals_a) ]\n",
      "    digits_list.append(digits)\n",
      "    \n",
      "for idx in xrange(N_eigvals):\n",
      "    \n",
      "    digit_data = [digit[idx] for digit in digits_list]\n",
      "    pl.plot(N_list, digit_data, \".-\", label = \"$\\widetilde{E}_{%d}$\" % idx)\n",
      "    \n",
      "pl.legend(bbox_to_anchor=(1.16, 0.7))\n",
      "pl.xlabel(\"$N^2$\", {\"fontsize\": \"xx-large\"})\n",
      "pl.ylabel(\"$l_n$\", {\"fontsize\": \"xx-large\", \"rotation\": \"horizontal\"})\n",
      "pl.xlim(None,1000)\n",
      "pl.show()\n",
      "\n",
      "for idx in xrange(N_eigvals):\n",
      "    \n",
      "    digit_data = [digit[idx] for digit in digits_list]\n",
      "    pl.plot(time_list, digit_data, \".-\", label = \"$\\widetilde{E}_{%d}$\" % idx)\n",
      "    \n",
      "pl.legend(bbox_to_anchor=(1.16, 0.7))\n",
      "pl.xlabel(\"$t$\", {\"fontsize\": \"xx-large\"})\n",
      "pl.ylabel(\"$l_n$\", {\"fontsize\": \"xx-large\", \"rotation\": \"horizontal\"})\n",
      "#pl.xlim(None,1000)\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# WYKRES ZALEZNOSCI ENERGII DLA KOLEJNYCH POZIOMOW ENERGETYCZNYCH:\n",
      "\n",
      "eigvals, eigvecs = EigenvalEigenvecProblem5(N, x_min, x_max, N_eigvals, potential, *args_potential)\n",
      "\n",
      "\n",
      "\n",
      "fig, ax1 = pl.subplots()\n",
      "ax1.set_xticklabels([])\n",
      "ax1.set_yticklabels([])\n",
      "ax1.set_xlabel(\"\")\n",
      "ax1.set_ylabel(\"$E_n$\", {\"fontsize\": \"xx-large\", \"rotation\": \"horizontal\"})\n",
      "ax1.set_xlim(x_min, x_max)\n",
      "\n",
      "ax2 = ax1.twinx()\n",
      "ax2.set_xticklabels([])\n",
      "#ax2.plot((eigvals), \"o-\")\n",
      "for eigval in eigvals:\n",
      "    ax2.plot([eigvals[0],eigvals[-1]], [eigval,eigval], \"b\")\n",
      "ax2.set_yticks(eigvals)\n",
      "ax2.set_yticklabels([\"$n = %d$\" % i for i in range(0,len(eigvals))])\n",
      "\n",
      "#pl.plot(x, x ** 2)\n",
      "\n",
      "#pl.savefig(\"num_har_E.png\")\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(x_min, x_max, N)\n",
      "\n",
      "# WYKRES FUNKCJI FALOWEJ:\n",
      "\n",
      "pl.figure(figsize = (10, 6))\n",
      "\n",
      "shift = 0.24\n",
      "\n",
      "for idx in range(N_eigvals)[0:3]:\n",
      "    \n",
      "    psi = eigvecs[:,idx]\n",
      "    pl.plot(x, psi + idx * shift, label = \"$\\psi_{%d}(x)$\" % idx)\n",
      "    \n",
      "pl.legend(bbox_to_anchor=(1.16, 0.7))\n",
      "pl.xlabel(\"$x$\", {\"fontsize\": \"xx-large\"})\n",
      "pl.ylabel(\"$\\psi_i(x)$\", {\"fontsize\": \"xx-large\", \"rotation\": \"horizontal\"})\n",
      "pl.xlim(x_min, x_max)\n",
      "#pl.ylim(-shift, shift * N_eigvals)\n",
      "pl.xticks([x_min, 0.0, x_max], [\"$x_{min}$\", \"$0.0$\", \"$x_{max}$\"], size = 15)\n",
      "pl.yticks([])\n",
      "#plt.savefig(\"num_harr_psi3.png\")\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": "*"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# WYKRES GESTOSCI PRAWDOPODOBIENSTWA:\n",
      "\n",
      "pl.figure(figsize = (10, 6))\n",
      "\n",
      "#shift2 = shift**2\n",
      "shift2 = 0.012\n",
      "\n",
      "for idx in range(N_eigvals)[0:4]:\n",
      "    \n",
      "    psi = eigvecs[:,idx]\n",
      "    pl.plot(x, psi**2 + idx * shift2, label = \"$P_{%d}(x)$\" % idx)\n",
      "    \n",
      "pl.legend(bbox_to_anchor=(1.16, 0.7))\n",
      "pl.xlabel(\"$x$\", {\"fontsize\": \"xx-large\"})\n",
      "pl.ylabel(\"$P_i(x)$\", {\"fontsize\": \"xx-large\", \"rotation\": \"horizontal\"})\n",
      "#pl.xlim(x_min, x_max)\n",
      "#pl.ylim(-shift2, shift2 * N_eigvals)\n",
      "pl.xticks([x_min, 0.0, x_max], [\"$x_{min}$\", \"$0.0$\", \"$x_{max}$\"], size = 15)\n",
      "pl.yticks([])\n",
      "#pl.savefig(\"num_harr_P3.png\")\n",
      "pl.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAr0AAAGBCAYAAABrdAmTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8leX9//H3SUIYIZAEwkwIJJBAGGFvWS5UXGhxUrSt\nYm2rrVZ/8NVaO7RfR23r+jpqW7UOXFBEFKpMkT0SCIQwAgkQAiEQAhmE5Pz+uD2QnbPvM17Px+M8\nIOdc57o/zLxz5XNfl8VqtVoFAAAABLAQswsAAAAAPI3QCwAAgIBH6AUAAEDAI/QCAAAg4BF6AQAA\nEPDCzC7A3SwWi9klAAAA2I2NtLwj4EKv5Nt/eZ588kk9+eSTZpcBD+HPFwB8i6//v8xinffQ3gAA\nAICAR+gFAABAwCP0etmkSZPMLgEexJ8vAPgW/l+GjSXQjiG2WCw+3dMLAABgQ27xHlZ6AQAAEPAI\nvQAAAAh4hF4AAAAEPEIvAAAAAh6hFwAAAAGP0AsAAICAR+gFAABAwCP0AgAAIOARegEAABDwCL0A\nAAAIeIReAAAABDxCLwAAAAIeoRcAAAABj9ALAACAgEfoBQAAQMAj9AIAACDgEXoBAAAQ8Ai9AAAA\nCHiEXgAAAAQ8Qi8AAAACHqEXAAAAAY/QCwAAgIBH6AUAAEDACzO7AAB+6NQp6ZNPpPJy6cYbpe7d\nza4IAIAmsdILwDE7dkgDB0pLl0pbtkhpadKCBWZXBQBAk1jpBWC/kyela66R/vd/pTvuMJ7bvFm6\n+mopJkaaMMHc+gAAaITFarVazS7CnSwWiwLslwT4jnvvlVq0kF55pfbzixdLP/uZsQocEWFObQDg\nh8gt3kPoBWCf3bulSy6RsrOlqKj6r99xh9Stm/Tcc96vDQD8FLnFewi9AOzzi19I0dHS73/f8OvH\njkn9+hl9vgkJ3q0NAPwUucV7CL0Amnf2rNSjh7RtmxQf3/i4uXOl4mLp1Ve9VxsA+DFyi/cQegE0\n7+9/lz7/XPrPf5oed/y4lJJi9PZ26+ad2gDAj5FbvIctywA07913pZ/8pPlxsbFGb+9rr3m+JgAA\nHMBKL4CmFRZKSUlSQYHUqlXz4zMzpcsvlw4eNHZ6AAA0itziPaz0AmjaokVGiLUn8EpS//5S797N\nt0IAAOBFhF4ATVuwQLr+esfec9990ptveqYeAACcQHsDgMaVl0udOkkHDhgnrtmrtNS4kW33bqlz\nZ4+VBwD+jtziPaz0AmjcunXG3ruOBF5JatNGuvZa6aOPPFMXAAAOIvQCaNzy5dLkyc699/bbpfff\nd289AAA4idALoHErVkiTJjn33ssuk/btk/bvd2dFAAA4hdALoGFlZdLmzdL48c69v0ULacYM6YMP\n3FsXAABOIPQCaNjatdLAgVLbts7Pceut0rx57qsJAAAnEXoBNOzbb6VLLnFtjjFjjEMtaHEAAJiM\n0AugYevXS6NGuTZHaKh03XUcVAEAMB2hF0B9Vqu0YYProVeSbrjBOOACAAATEXoB1JeTI4WHS3Fx\nrs916aXStm3S8eOuzwUAgJMIvQDqc0drg02rVtIVV0iLFrlnPgAAnEDoBVCfO0OvZLQ40NcLADAR\noRdAfRs3SiNHum++q682TncrLXXfnAAAOIDQC6C26mopI0MaPNh9c0ZHGyF6yRL3zQkAgAMIvQBq\n279fiokxgqo7XXed9Pnn7p0TAAA7EXoB1JaeLqWluX/eadOkL74wVpIBAPAyQi+A2jwVenv1kmJj\njf1/AQDwMkIvgNo8FXol6dpr2boMAGAKQi+A2rZt82zopa8XAGACi9VqtZpdhDtZLBYF2C8J8J6T\nJ6UePaTiYinEA18TV1VJXbpImzcb1wGAIEdu8R5WegFclJEhDRjgmcArSaGhxp69tDgAALyM0Avg\nIk/289pMm0aLAwDA6wi9AC7yRui98kppzRrpzBnPXgcAgBoIvQAuysjwfOht104aNUr6+mvPXgcA\ngBoIvQAM1dXSrl1Saqrnr0WLAwDAywi9AAyHDkmRkVJUlOevde21nM4GAPAqQi8AQ1aW1K+fd66V\nmCh16CBt2uSd6wEAgh6hF4Bh1y6pb1/vXY8WBwCAFxF6ARi8udIrcSQxAMCrCL0ADLt2eTf0jh4t\n5eUZDwAAPIzQC8Dg7faGsDDpqqtY7QUAeAWhF4BUVCSVlUndu3v3urQ4AAC8hNALwOjn7dtXsli8\ne90rr5RWr5bOnvXudQEAQYfQC8D7rQ027dtLI0ZwOhsAwOMIvQC8v3NDTbQ4AAC8gNALwPs7N9Q0\nbZoRejmdDQDgQYRewJfs2yctXChVVHj3uraeXjP07m0cfbx5sznXBwAEBUIvUNf8+VJqqjRjhnTq\nlPeu+9130qhR0h/+IN1wg1RV5Z3rlpdLhw9LSUneuV5DbKu9AAB4CKEXqGnjRmn2bOmVV6SICOPn\n3lBVJf34x9Kbbxrh9/Rp6e23vXPt7GypVy+pRQvvXK8h117LkcQAAI8i9AI2Vqv0q19Jzz0nTZ4s\nvfqqtG6dtGGD56/92WdSTIyxwtuihVHDH//onT5XM1sbbMaOlQ4eNFacAQDwAEIvYPPtt9Lx49Kd\ndxoft24tPfSQ9MILnr/2P/8p/fSnF/fJHTPG6HP1xlZeZt7EZhMWJk2dSosDAMBjCL2AzWuvSfff\nL4WGXnzuzjulL7802g08pahIWrNGuvHGi89ZLNKsWdIHH3juujZmbldWEy0OAAAPIvQCklRaaqwy\n3nFH7ec7dJAmTPBsGFu2TBo3zughrumGG4zrnj/vuWtL5h1MUdeVV0qrVhl/FgAAuBmhF5CkpUul\n4cOljh3rv3bttdLixZ679n//K11+ef3nExKkuDhp/XrPXbuqyriRzRdCb3S0NGwYp7MBADyC0AtI\nxjZlNdsLarrqKmnJEs9tIdZY6JWkKVOkFSs8c11Jys01gn7btp67hiNuukn6+GOzqwAABCBCL1BZ\nabQ23HBDw6/Hxxs7K+zc6f5r79tn7JPbv3/Dr0+eLC1f7v7r2vhKa4PNTTcZfxbl5WZXAgAIMIRe\nYO1aqWdPo5WgMWPHGjebuduyZdKll17ctaGuSy4x2hs8dUKbL+zcUFPXrtLgwdJXX5ldCQAgwBB6\ngVWrjBXVpowb55nQu369NHp0469HRUkpKZ7bK9hXdm6o6ZZbpHnzzK4CABBgCL3AypXSxIlNjxk3\nzjgpzd02bpRGjmx6zKRJnuvr9bX2BslocfjyS3ZxAAC4FaEXwa2y0jh1bfz4psf17SudPCkdPeq+\na5eWSnv3SoMGNT1uzJjgWumNjZVGjPDsjhkAgKBD6EVw27xZSkoytstqSkiIET7d2eKwdauUmiq1\nbNn0uBEjjNBrtbrv2pJx+lx1tdSpk3vndQdaHAAAbkboRXBbudI4fMIeI0caIdldNm40Am1z4uON\nG93y8tx3belia0NjN9GZafp0Y+/kkhKzKwEABAhCL4LbqlXN9/PaDB5srM66i72h12K5uNrrTr7Y\n2mATE2P8uXz6qdmVAAACBKEXwauqymhXsHeld8gQ94berVuloUPtGztypBGS3cnXtiur6667pLff\nNrsKAECAIPQieKWnS926GTdO2SM+3rjxzR03s1VUSDk59u+cMGKE+0NvVpbv7dxQ0zXXSDt2GL9P\nAAC4iNCL4GXPVmU1WSzua3HYtUtKTGz+Jjab4cONfuLqatevXbMGX17pbdlSuvVW6Z13zK4EABAA\nCL0IXo7cxGbjrhaHHTukgQPtH9+xo9S+vbR/v+vXlqSzZ6Vjx4yT6HzZrFlGi4M7wz4AICgRehGc\nqqul1asdW+mVjNC7bZvr19++3bHQKxn9v+7qKc7Olnr3lkJD3TOfpwwbJkVEeO5wDgBA0CD0Ijhl\nZho7BHTr5tj7Bg92X+gdMMCx97jzRrqdO327tcHGYpFmz5b+7//MrgQA4OcIvQhOq1Y53togScnJ\nxn65rh6R6+xK75Ytrl3XZudOqX9/98zlaTNnSl9/LeXnm10JAMCPEXoRnBy9ic2mRQspJcUIjc46\ndcp4ONpPO2SIEXrdcTJbZqZxGpw/aN9emjFD+vvfza4EAODHCL0IPlarczex2QwcKGVkOH/9HTuM\nVdYQB//5de9u1H7kiPPXtsnM9J+VXkn66U+lN96Qzp83uxIAgJ8i9CL4ZGdLrVo5v3PBoEFGe4Kz\nnGltkIz+VnfczFZWJh06ZNzI5i8GDzb2Sf7Pf8yuBADgpwi9CD7OtjbYuLrS68xNbDbuuJlt924p\nKclo1fAnv/619Nxz7mnvAAAEHUIvgo+zN7HZDBpkhF5nw5ezK73Sxb5eV+zc6T/9vDVdf7104oT0\n7bdmVwIA8EOEXgQXWz+vKyu9Xbsa+/wWFDh3fUcPpqjJHe0N/tbPaxMaKj38sLHaCwBBKCYmRhaL\nhUcjj5iYmCZ//wi9CC779hmB1ZV+VovF+b7ew4eN43VjY527dlKSVFRkPBpgtWf12V9XeiXjhLb1\n643gDgBB5uTJk7JarTwaeZw8ebLJ3z9CL4LLsmXSpZcawdUVzvb1utLPKxk7PgweXG+112q16onl\nT6jtn9pq6OtDtbtwd+Nz+NN2ZXW1bm2s9v7ud2ZXAgDwM4ReBJdly6QpU1yfx9mVXlf6eW0auJnt\n/zb9nxbuXqjsn2dr9rDZuuLfV+hU+an67y0vl3JzpT59XKvBTD/7mXGEtDtOxgMABA1CL4JHdbUR\neidPdn0uV1Z6Bw1y7dp1Qu/JspP67Yrf6sObP1T3dt01e/hsXd3naj227LH67929W0pMlMLDXavB\nTBER0ty50hNPmF0JAMCPEHoRPDIzpXbtpIQE1+fq31/KynL8sISMDNdDb53jiN/c8qau7nO1+nbs\ne+G530/6vT7c8aHyivNqv9d2MIa/u/deI/ivXWt2JQAAP0HoRfBwV2uDJLVtK3XrJu3ZY/97zp0z\nDsZwtZ+2Xz/p4EHp7FlZrVa9ueVN/XT4T2sNiY2I1W0DbtNbW9+q/d70dKMn2N+1aiU99ZT04IPG\nCj4AAM0g9CJ4uDP0So739e7ebZwC17q1a9dt0cIIzunp2py/WZI0qvuoesN+MvQn+ue2f6qquuri\nk9u2BUbolaQ77zS2MXv7bbMrAQD4AUIvgsP588b+vO7o57VxtK/XHa0NNt/v1/vZrs90U7+bZGlg\nN4rBXQarY5uO+nr/18YTVmtghd6QEOnFF6X/+R+puNjsagAATsrJyWl2TH5+vkpLS126DqEXwWHN\nGuMGrs6d3Tenoyu97gy939/MtnTfUl3d5+pGh909+G69m/Gu8UF+vhF8u3VzTw2+YMQI6dprpTlz\nzK4EAOCE/fv3a926dc2Oi42N1bPPPuvStQi9CA6LFhnhyJ1MXumt2rxRWYVZDbY22Fyfcr2+3Pul\nzlefv7jK6+oexb7m2WelL76Qvv7a7EoAAA56/fXXddtttzU7LiwsTNdcc43eeecdp69F6EVw+Pxz\n94fepCTp2DHp9Gn7xrsz9A4cKGXt1tjOw9UyrGWjw+LbxyuhfYK+y/susFobaoqKkt58U/rxj+3/\nswAAeMymTZt09dVXa8KECXrrrbf0+uuv6/7779eKFStqjUtPT1dcXJzd844YMUJfu7DAQehF4Nuz\nxwhDQ4e6d97QUOOGsh07mh9bWCidOSP16OGea7dpo+Od2+om9Wt26HUp12nh7oXGFl9pae65vq+5\n8krjcd99RgsHAMA0w4cPV5s2bXTPPffoxz/+sWbPnq2f/exnuuWWW2qNW7RokaY4eIN5bGys9u7d\n61RdhF4EvgULjFXeEA/8dbe3r9d2KIUbWwu2dJEmFbVvdtx1Kdfp8+zPpQ0bpJEj3XZ9n/PXv0o7\nd0qvvGJ2JQAQ1KxWq1auXKlx48ZdeG7//v2KjIysNW7jxo1KdXAbz7S0NG3evNmpusKcehfgT957\nzwhEnmBvX687WxsknT13VqtiTuuKnOZ3LRjSZYjCjxaqqqRaof58/HBz2rSRPv1UGjPGWNUfO9bs\nigDA69yxtuLqN8wyMjIUFhamxMRESVJZWZneeOMNvfzyy7XGlZaW1tt9aOHChQoNDdXq1as1cOBA\nffXVV3rsscfUt69xAFN0dLSys7OdqovQi8CWmSmdOCFNmOCZ+QcNkj77rPlxW7ZI48e77bIZBRk6\nPqCXwpZtbHasxWLRXRV9ld//jOIC7Sa2upKSjH17p083tqhLSTG7IgDwKl/o8Fq+fLl69OihefPm\nqbKyUiUlJXr55ZeVUOdE1Kqqqlof5+bmKjU1Vb1799YTTzyhOXPmqH379upRozWwdevWOnfunFN1\nEXoR2N57T7rtNs+0NkhGj+y2bVJVldHj25h166Rf/9ptl916dKtajBorvfqxVFIi1fmWUV2XFrTV\nuh6lutltFfiwq66S/vQnaepUY6u6QNqiDQD8wPLlyzVr1qxaPbxFRUV65pln1KlTJw0aNEjDhg1T\nWFjtGGoLtwUFBYqMjFRUVJSmTZtWa0xxcbFiYmKcqoueXgSuykrpnXekmTM9d40OHYxQlZnZ+Jii\nImOPXFePH65h69GtGhg/XBo2TPruu2bHJ+8t0qftDsnqC0sA3nD33dLs2cZhJHl5ZlcDAEGjurpa\nq1ev1tSpU2s9//bbb2vy5MmaOXOmXnjhBUlSly5ddObMmQtjsrKylJ6ersWLF2vC99+hXbRoUa15\n8vPz1bt3b6dqI/QicH3yidSnj9F360ljxzYdPDdskIYPb3ol2EHbjm7TkK5DjLaN1aubHlxZqdYZ\nu7Q1voX2Fjl3x6tfmjPHCL6XXGLs4AEA8Kj09HTNnTtXFRUV9bYn279/v7p27aqwsDAVFRVJkiZO\nnKgNGzZcGLN06VItWrRIVqtV5eXlmj9/vjp16lRrnm3bttW6Qc4RtDcgcP3tb945qWvsWGnVKmO7\nrIasWyeNHu22y1VWVWrn8Z0a1HmQdMkZ6Y9/bPoN69bJkpKi4X37aeXBlerTIYBvZqvroYekdu2M\nfur33pMuu8zsigAgYKWlpSktLU3PPPNMvdeqq6sV+v3ij+3mtenTp+v555+/sG3ZAw880OT85eXl\nateunVq1auVUfaz0IjAtX27sjevuAykaMmZM0yu9bg69WYVZim8Xr7bhbY1rb9kiNXUe+VdfSVdd\npbHxY7X20Fq31eE3fvITad48o83lmWeM/msAgFelpKSooKDgQnCVpKioKHXs2FGFhYV2zfHhhx9q\n9uzZTtdA6EXgqa42bhp76im3thQ0ql8/Y4eIgoL6r50/L61fL41q/KhgR207uk2Du3x/slpkpLH3\nblMn1Hz1lTR1qsbEjdHavCAMvZI0aZLxxceXX0oTJ0pObmwOAHDOHXfcoa+//lrvvPOOHnzwwQvP\nP/jgg5o/f36z78/Ly1N0dLRSXNiVh/YGBJ5335VatJBmzPDO9UJCpHHjjC2y6l5z82YpLk7q3Nlt\nl9tZuFP9Y/tffOKGG6T586Xrrqs/uKBA2r9fGj1aA0OkvNN5KiorUkxr5+589WsJCdKyZUbby+jR\n0j33SHPnGu0PAACP6tChgx555JF6z1ssFt1zzz3Nvj8+Pl7x8fEu1cBKLwJLTo70yCPSq6+69fSz\nZk2daqwi1vXf/0qXX+7WS2UVZqlfbI3jh2+4QVq0yFhVrmvJEmnKFCksTGEhYRrZfaTWHVrn1nr8\nSkiI9KtfGYeF5OdLvXtLTz5ptMIAAAIaoReBo6REuuUWY/Vu6FDvXvuqq4zQW7dfdMEC6Zpr3Hqp\nrMIs9e3Y9+ITPXoYj2+/rT/4gw+kmy/uzhu0fb11desm/etfxs4Xhw9LyclG7+/KlUZ7DAAg4BB6\nERiKi41wmZYm1egV8pqkJKl7d+MGOpv9+6XcXKOH1E0qqyqVczJHfWLq7MBw663SP/9Z+7m8PKOP\n9frrLzw1Jm6Mvstrfl/foJGSIr35prHPckqK9ItfSL16GT9+8UXTNwgCAPwKoRf+b9kyafBg4/H6\n6547fa05P/yh9Pe/X/z49delO+6QwtzXOr/v5D7Ft49Xy7CWtV/40Y+khQulQ4cuPvfCC8YhDW3a\nXHhqdNxobTy8UeerG2iFCGZduxptMRkZRqtIXJz0/PNSbKyxx/L99xtfVGzaZHyBBQDwOxZrgB3R\nZLFYgufUqWBVVWUcNrBihfTvf0tHjxoBr6EbubyppMRY8V261AhLaWnSxo3GyqGbzN81X//Y9g99\nftvn9V/8zW+kXbukjz+Wduwwenm3bTNWoGtIfSVV79/0/sUdINC4sjLj93DDBuOxc6fxd69NG6Ml\nIiHBuEmxSxfj0bmzFB1t3BwXGWn8GBFh3hdiAHyeI7mFjNO05n5/2L3B2775RvrsM+Pndf9gan5s\nz8+deY8Z13THeyorjeN8jx83VjM7dzbuwH/oISPsunE11WmRkdJf/mK0WUREGDdMuTHwSt/fxNax\nX8Mvzp0rXXqpEXazsqS//rVe4JWM1d61eWsJvfZo3drYC3nMmIvPWa3GTXB79hgtJEePGo+MDGO3\njFOnpNOnjS+CTp82gnNEhBGUw8ObfoSFGTdghoQYj5o/r/txQ6/Ze/OmPePcNcaMuYC6Onc2FgYQ\n1HwgKQSZDh2k1NSLH9f9D7zmx/b83Jn3mHFNV98TGip17Gj8/sXF+e42U3fcYewIUFJiBFA3yzqR\npYkJjfQIt2ljtHp8/rlx/PKQIQ0OG9FthDblb3J7bUHDYjFuhOvWzb7xVVXSmTNSeblUUSGdO1f/\nUVFhPKqrjYfV2vDPm3vNHvasErlrjBlzAQ2JijK7AvgA2hsAPzL676P15yv+rHE9nDt3XJI2Ht6o\nn3z+E6Xfl+7GygAAzqC9wX2a+/2h0QzwI3uL9qp3TG+X5hjUeZD2nNij0kp2JgAABA9CL+AnisuL\nVVFVoU4RnVyap2VYS6XGpir9KCu9AIDgQegF/MS+k/uUGJ0oixtu5Bnebbg2HtnohqoAAHBNTk5O\ns2Py8/NV6uLe6YRewE/sK9qnpOgkt8w1otsIbTrCzWwAAHPt379f69ata3ZcbGysnn32WZeuRegF\n/MS+k/uUFOOe0Du823BCLwDAdK+//rpuu+22ZseFhYXpmmuu0TvvvOP0tQi9gJ/Yd9J9K72psak6\nWHxQJRUlbpkPAACbTZs26eqrr9aECRP01ltv6fXXX9f999+vFStW1BqXnp6uuLg4u+cdMWKEvv76\na6frYp9ewE/sK9qnGakz3DJXi9AWGtR5kLbkb9HEno3s+wsAgBOGDx+uNm3a6LbbbtPMmTMlSZmZ\nmZoyZYoKCgoujFu0aJFuuOEGh+aOjY3V3r171bu34zsZsdIL+Al3tjdI9PUCADzDarVq5cqVGjfu\n4p7y+/fvV2RkZK1xGzduVGrNA7vskJaWps2bNztVFyu9gB0+2fmJnlzxpK5NuVZPTXlKIRbvfr1Y\ncb5CR88cVY/2Pdw25/Buw7V4z2K3zQcAMJ/ld67v8GP9rWsHYGRkZCgsLEyJiYmSpLKyMr3xxht6\n+eWXa40rLS2ttyPRwoULFRoaqtWrV2vgwIH66quv9Nhjj6lv376SpOjoaGVnZztVF6EXaEZWYZbu\nW3Sf/j393/qfb/5HKR1SdNfgu7xaw8Hig4prF6ewEPf9kx3ebbh+v/L3bpsPAGA+VwOrOyxfvlw9\nevTQvHnzVFlZqZKSEr388stKSEioNa6qqqrWx7m5uUpNTVXv3r31xBNPaM6cOWrfvr169Li44NO6\ndWudO3fOqboIvUAzfrP8N5ozfo6m9p6q6FbR+sHHP9DMQTMVGhLqtRpyi3OV0D6h+YEOSOmQomNn\nj+lk2UlFt45269wAgOC1fPlyzZo1S7fccsuF54qLi/XZZ59p9+7dmjt3riRjR4aabOG2oKBAkZGR\nioqK0rRp02qNKS4uVkxMjFN10dMLNCG/JF9f7/9a9w67V5I0Km6UurTtoiX7lni1jrziPLe2NkhS\naEiohnQdQl8vAMBtqqurtXr1ak2dOrXW8+3bt9ewYcNqrdJ26dJFZ86cufBxVlaW0tPTtXjxYk2Y\nMEGScbNbTfn5+U7dxCYReoEmzcucp+tTrle7lu0uPDdz0Ex9vPNjr9aRW5zr9tArsV8vAMB90tPT\nNXfuXFVUVNTbnqwhEydO1IYNGy58vHTpUi1atEhWq1Xl5eWaP3++OnXqVOs927Ztq3WDnCNobwCa\nsCBrgR4e83Ct565JvkZ/XP1HVVurvXZDW+7pXI2JG+P2eYd3He71AA8ACExpaWlKS0vTM888Y9f4\n6dOn6/nnn9eUKVMkSQ888ECT48vLy9WuXTu1atXKqfpY6QUaUVxerM35m3VZ4mW1nk+MTlRM6xht\nO7rNa7V4aqV3RHe2LQMAeIfVWvsmu6ioKHXs2FGFhYV2vf/DDz/U7Nmznb4+oRdoxJq8NRrRbYRa\nt2hd77UJCRO0+uBqr9XiqdCbFJ2kknMlKjhT0PxgAACcdObMGX366afavHmzduzYceH5Bx98UPPn\nz2/2/Xl5eYqOjlZKSorTNVisdWO3n7NYLPW+kgCcMefrOWoV1kpPTnqy3mv/zvi3FmQt0CczPvF4\nHVarVW2ebqPCRwoVER7h9vkvf/dy/XLUL3VN8jVunxsA0DRHcgsZp2nN/f6w0gs0YnXuak1ImNDg\na2Pjx2rtobVeqaOwtFARLSI8EnglbmYDAAQHQi/QgNLKUqUfTdfouNENvt4rqpfOnjur42ePe7yW\n3OJcxbeP99j8w7sO18YjGz02PwAAvoDQCzRg/aH1Gth5oNq0aNPg6xaLRWld0pRekO7xWjzVz2tj\nu5mNb5kBAAIZoRdowLpD6zQ2fmyTY4Z0GaKt+Vs9XounQ298u3hVW6t1uOSwx64BAIDZCL1AA7Yc\n3aJhXYc1OWZwl8HaVuD5bctyT+eqRzvPhV6LxaIR3Udo42FaHAAAgcuU0JuTk6O3337bofc8/fTT\ntY6uAzxpS/4WDe06tMkxg7sM9spevZ5e6ZW+v5ktn5vZAACBy+uh98SJE3r88cd15513OvS+m2++\nWffdd581lNF/AAAgAElEQVSHqgIuOll2UsfOHlOfmD5NjkuNTVXOyRyVVZZ5tB5vhN4R3VjpBQAE\nNodC74oVKzRjxgyFh4crNjZWd955p2bOnKmbb75Zw4YN0yOPPKLy8vIm55gzZ44eeeQRhYaGOlRo\ncnKy+vfvr3feeceh9wGO2nZ0m9I6pyk0pOm/o+Gh4UrukKzM45kercdrK73czAYACGBhjgyeNGmS\nxo8fr6ioKP385z/Xb3/72wuvlZaWasiQITp+/Lj+9a9/Nfj+Xbt26cCBAxo8eLBTxd57770aP368\nbr/9doWFOVQ6YDd7Whts+nbsq6zCLA3vNtwjtVScr1BRWZG6tO3ikflturTtoojwCO0/uV9JMUke\nvRYAAGZwuL1h48aNKi0t1SWXXFLr+TZt2mjIkCH67LPPGn3vSy+9pLvuusvhIm0iIyM1ZMgQ/ec/\n/3F6DqA5W47aH3pTOqZo94ndHqvl0OlD6hbZrdlVZ3fgkAoAQCBzOPSuWrVKoaGhGj26/qb92dnZ\n6tSpU6PvXbJkicaMGePoJWsZP368FixY4NIcQFPSj6YrrXOaXWNTOqRod2HjobeiQqqqcr4Wb7Q2\n2HBIBQDADDk5Oc2Oyc/PV2lpqUvXcTj0rly5UkOGDFGbNrU37c/MzNS2bdv0wAMPNPi+AwcO6OTJ\nk0pMTKz32t/+9jf96U9/0g9+8APt2rVLTz31lP74xz/q4Ycfrjd2+PDhWrdunaNlA3aprKrUvpP7\n1LdjX7vGp3RofKX35ZelmBipd29p+3bn6vFm6LUdUgEAgLfs37/frlwXGxurZ5991qVrORR6q6qq\ntGbNGo0fP77W89u3b9ftt9+uxx9/vNHQm5ub2+Aq8IsvvqipU6dq7ty5uuKKKzRp0iTdfffdqqio\n0L///e964zt27Kj9+/fr/PnzjpQO2GVP0R7FtYtT6xat7Rqf3CFZe07sUbW1utbz334r/elPUmam\n9Mc/StdeKznzBao3Q++wrsO0JX+LqqpdWJoGAMABr7/+um677bZmx4WFhemaa65xaUMDh+4G27p1\nq0pKSpSTk6O5c+dKkk6dOqXo6Gh99NFHSklJafS9x44dU1RUVL3nz58/f+F9hw4d0pAhQ9StWzfN\nnj1bM2fOrDc+JiZGVqtVJ0+eVGxsrCPlA83aeXyn+sf2t3t8ZMtIRbeOVl5xnhKiEi48/8gj0vPP\nSz17Go+FC6W//U36/p+N3XJP52p4V8/cJFdXhzYdFBsRq+wT2eoX288r1wQABJ5NmzbpiSee0Jkz\nZzRr1iydP39e6enpmjFjhiZNmnRhXHp6uuLi4uyed8SIEXrppZf0wx/+0Km6HAq9K1eulGSszsbH\nxzt0oerqaoWE1F9Yfuihhy78fNWqVbryyislqdHfBNscFovFoesD9sg8lqn+newPvdLFFgdb6F2/\nXiookGbMuDjm8celK6+UHn5YCg+3f+7c4lxN7zvdoXpcMbyb0ddL6AUAOGv48OFq06aNbrvttgsL\nmJmZmZoyZYoKCgoujFu0aJFuuOEGh+aOjY3V3r171bt3b4frcqi9YeXKlUpKSmo28D7++OP1dliI\njY1VUVFRo++pqKjQ+vXrNXHixCbnLioqUmhoqDp27Gh/4YCdMo9nOrTSK32/g0ONm9nefVf68Y+l\nmltRDxxo9PYuXuxYPd5sb5CMQyro6wUAuMJqtWrlypUaN27chef279+vyMjIWuM2btyo1NRUh+ZO\nS0vT5s2bnarL7tBbXV2tb7/9VpdeemmzY//whz/o+uuvr/Vc9+7ddeLEiVrPVVZWatmyZZKkdevW\nyWq1asSIEZKMk9saOqr4xIkT6ty5s71lAw7JPJ6p1FjH/gEmxyQruyhbkrFTwyefSLfcUn/cHXdI\n779v/7xWq1W5xbmKb+/Yd1VcYVvpBQD4KYvF9YeLMjIyFBYWdmHzgrKyMr3xxht6+eWXa40rLS2t\n9537hQsX6osvvtCcOXP03nvvaebMmcrKyrrwenR0tA4dOuRUXXa3N2RkZOjUqVOaPHlys2Mbaj1I\nTk5WeHi4Dh8+rO7du0uS3njjDT388MM6efKkFi5cqI4dO144dOKVV17Rz3/+83rzbN26VaNGjbK3\nbMBu56rOaV+R/Ts32CTFJOnrnK8lSenpUlSUsapb1/Tp0qOPSpWVUosWzc9bVFak8NBwtWvZzqF6\nXDGs6zBlFGToXNU5hYc60IcBAPANPnCy5vLly9WjRw/NmzdPlZWVKikp0csvv6yEhIRa46rq7OmZ\nm5ur1NRU9e7dW0888YTmzJmj9u3bq0ePi9/xbN26tc6dO+dUXc2G3gMHDmjOnDnaunWrLBaLXn31\nVWVnZ+s3v/lNvbGFhYVaunSpFixYoI8++qje61dccYVWr16tW2+9VZI0ceJETZ8+XU8//bSmT5+u\niIgIPfzww2rdurVuvvlmxcTE1Jtj9erV9VaRAXfYW7RXPdr3UKuwVg69LzE6UftP7pckLVsmNfbN\nkNhYqU8fad06qc7ZLg3ydmuDZNyY1yemj7bmb9WoOL64BAA4bvny5Zo1a5ZuqfFtzz179mj+/PnK\nyMjQtddeq6FDh9Y7XdcWbgsKChQZGamoqChNmzat1pji4uIG86E9mg29PXv21IcffmjXZFu3btXV\nV1+tv/zlLw2+PmvWLL344osXQu+AAQP0fo3v99bs/WhIeXm5tmzZUm95HHAHZ25ik6ReUb104NQB\nVVurtXx5iH70o8bHXnGFtHSp74ZeSRobP1bf5X1H6AUAOKy6ulqrV6+ulwUXLVqkcePG6bLLLtPs\n2bP1/vvvq0uXLjpz5ozatm0rScrKylJFRYW2bNmiCRMmXHhfzeCbn5+vfv2cu9na4cMpmnL55Zfr\nH//4R6NHDU+aNEllZWXKzs52av63335b99xzj1q2bOlClUDDnLmJTZIiwiPUrmU75Z08qm+/lWrs\nxlLPFVdIS5bYN69ZoXdM3Bh9d+g7r18XAODf0tPTNXfuXFVUVGjFihW1XvvVr36lkSNHKi8vT716\n9ZJkfMd/w4YNF8YsXbpUixYtktVqVXl5uebPn1/vjIdt27Y1u0jaGIe2LLPHe++9p2XLltVL5jav\nvfaaHn30UX3wwQcNbmHWmPz8fP33v//VvHnz3FkucEHm8UyntwdLjE7Ul+v2KzGxmzp0aHzc2LHS\n7t3SiRNqcpxk7NHbo505K71zvpkjq9XK1oAAALulpaUpLS1NzzzzTKNj5s+fr8cee0ySNH36dD3/\n/POaMmWKJDV6wJlNeXm52rVrp1atHGtDtHHrSq8k9e7dW59//rlGjx7d4OsJCQn69a9/rVdeecWh\neZ977jm99tprCq25DxTgRpnHHN+5wSYxOlEr0vc327YQHi6NGSOtWdP8nGat9CZGJ+p89Xnlnc7z\n+rUBAIFr4cKFeuCBB3T48GFJUlRUlDp27KjCwkK73v/hhx9q9uzZTl/f7aF33rx5uvPOO5vcR3fE\niBH6xS9+4dC8L7zwAnvzwmPOVZ1TzqkcpXRs/FTBpiRGJyrzSI6G23F42ujRxgEWzckrzjMl9Fos\nFqPFIY8WBwCAe8yfP19/+MMfNH369FqbHTz44IOaP39+s+/Py8tTdHR0k6f/Nsft7Q2AP9pzYo9T\nOzfYJEYlKvf0Cg0d2vzYUaOMI4qbY9ZKr2S0OKw9tFa3DrjVlOsDAALLjTfeqBtvvLHe8xaLRffc\nc0+z74+Pj3f4NOC63L7SC/gjZ29is+nSMlFnWuxXXzu2+B05Utq40TjIojGVVZU6XnpcXSO7Ol2T\nK2w7OAAAECgIvYC+D71ObFdmU3okUaGx+xVmx/dOOnSQunaVdu5sfMzhksPq0raLwkLM+WbMsK7D\ntPP4TpVWlppyfQAA3I3QC+j7PXpdWOnNy+ym6pYnVH6+3K7xo0YZh1Q0xszWBklq3aK1BnQaoE1H\nNplWAwAA7kToBWSs9Dq7c4Mkbd0Sqo4teujAqQN2jR85UqqxNWE9ZodeiRYHAEBgIfQi6J2rOqec\nkzlK6eD8HaHp6VKvqIvHETdnyBDjPY3xhdA7Ln6cVueuNrUGAADchdCLoJd9Ils9o3qqZZhzJ/1V\nVUnZ2dKA7vaH3kGDpMxM6fz5hl/PLTbnYIqaJiRM0JrcNTpf3UiRAAD4EUIvgl7mMdduYjt4UOrY\nUUrpZH/ojYw0bmbbs6fh131hpbdTRCfFtYvTtqPbTK0DAGCIjo6WxWLh0cgjOjq6yd8/Qi+Cnqvb\nle3cKaWmSr2ieinnVI7d70tLa7zFIbc4V/HtXduP0B0m9ZykFQdWmF0GAEBSUVGRrFYrj0YeRUVF\nTf7+EXoR9FwNvbt2Sf36Sb2ie9m90isZoXdbA4uoVqtVB4sPmr7SKxF6AQCBg9CLoLfz+E6Xdm6w\nhd7E6ETlnMyR1Wq1632DBze80ltcUSyLLGrfsr3TNbnLxISJ+jb3W/p6AQB+j9CLoFZ+vlwHTh1Q\nSkfnd26wtTdEtYpSaEioisqa/vaKTWMrvbZ+XovF4nRN7hIbEav49vHamr/V7FIAAHAJoRdBLasw\nS0nRSQoPDXfq/VbrxZVeyejrtbfFoUcP6exZqW4Lki/cxFYTLQ4AgEBA6EVQ216wXQM7D3T6/fn5\nUsuWxtHC0vctDnbezGaxGGF5167az/ta6J3cc7K+yfnG7DIAAHAJoRdBbfux7RrYyfnQu2ePlJx8\n8eNe0b2Uc9L+HRz69TPaI2rytdA7pdcUrclbo7LKMrNLAQDAaYReBLUdx3ZoQKcBTr9/3z6pd++L\nHzu6bVlqav2VXl/ZucEmqlWUBncZrFUHV5ldCgAATiP0Iqi5utK7d2/t0JsYbf8BFZIReuuu9B44\ndUC9ono5XZMnTE2aqq/2fWV2GQAAOI3Qi6B1qvyUTpWfUkJUgtNz1A29zqz0NhR6e0b1dLomT5ja\ne6q+2kvoBQD4L0IvgtaOYzvUP7a/QizO/zOoG3oTohKUV5ynquoqu96fkCCdOCGVlBgfl58v18my\nk+oa2dXpmjxhSNchOlF6QgdOHTC7FAAAnELoRdDacWyHSzs3WK1G6E1Kuvhcq7BW6timow6XHLZr\njtBQ40a4rCzjY9vxw64EcU8IsYToyt5XasneJWaXAgCAU3zrMyvgRduPbdeAWOdvYjt+XAoPl6Kj\naz/v6A4ONVscDpw6oIT2zrdbeNLUpKlavHex2WUAAOAUQi+Clqt79NZtbbBxpa/XF/t5ba7qc5WW\n5yzX2XNnzS4FAACHEXoRlKxWq7HS68J2ZY2G3mj7T2WTam9b5suhN6Z1jEbHjeaGNgCAXyL0Iigd\nOHVAES0i1Cmik9Nz7NsnJSbWfz4xyv5T2aTaB1T4cuiVpJv63aRPd31qdhkAADiM0IugtCV/i4Z2\nHerSHAcPSj171n/e0Z7epCTp0CGpvNz3Q+/1fa/Xl3u/VMX5CrNLAQDAIYReBKUtR10Pvbm5xpZj\ndfWKcqy9oUULIzzv2+f7obdL2y4a0GmAvsn5xuxSAABwCKEXQckdK725uVKPBk4L7hbZTUVlRSqr\nLLN7rj59pJ27K1RUVqSubX1rj966pvedTosDAMDvEHoRdKxWqzYf2exS6K2uNloS4uPrvxYaEqoe\n7XvoYPFBu+dLTpY2Zucqrl2cQkNCna7LG25OvVkLshao/Hy52aUAAGA3Qi+CzpGSI5Kk7pHdnZ6j\noECKipJatWr4dUd3cEhOlrbn+XZrg018+3gN6TJEC3cvNLsUAADsRuhF0Nl0ZJOGdh0qi8Xi9BwH\nDzbc2mCTGJ3o0M1sycnS3hP+EXolaVbaLL2d/rbZZQAAYDdCL4LO2kNrNTputEtzNHYTm42jB1Qk\nJ0v5pf4Teqf3m67v8r5Tfkm+2aUAAGAXQi+CztpDazU2fqxLczS30uvoDg7duknlrXPUKdw3jyCu\nKyI8Qjf0vUHvbX/P7FIAALALoRdBpbKqUpuPbNao7qNcmqe5ld7EaMcOqLBYpPAuexV+po9LdXnT\nXWl36R9b/yGr1Wp2KQAANIvQi6CSXpCunlE91b5Ve5fmaXal9/sb2ewNhFarVefbZ+vcUf8JvRMS\nJijEEsKevQAAv0DoRVBZm7dWY+LHuDxPcyu90a2iFWIJ0YmyE3bNV1haqLCQUOXv6+Bybd5isVj0\ny9G/1F/W/cXsUgAAaBahF0Hlu0PfaWyca/28UvMrvRaLRSkdUrS7cLdd82WfyFa3Vn2Une1yaV51\nx8A7tPHwRrt/nQAAmIXQi6BhtVq1LGeZJvea7NI8p09LlZVSTEzT4/p27Kuswiy75txTtEd9Yvpo\nzx6XSvO61i1aa/bw2Xpxw4tmlwIAQJMIvQgaO47tUGR4pMvbgtmOH25um9++Hfsq64T9oXdwj2Rl\nZ0v+dl/Y/cPv1wfbP2D7MgCATyP0Imh8k/ONLk281OV5muvntXG0vSEtro9CQ6Xjx10s0Mu6RnbV\n3UPu1lOrnzK7FAAAGkXoRdD4JucbXdrL9dDbXD+vjUPtDSf2KLlDspKT5Xd9vZI0Z9wcfbDjA4dO\noQMAwJsIvQgK56rOafXB1Zrc07V+Xsn+ld6kmCTlFueq4nxFk+OsVqv2Fu1Vn5g+fht6YyNi9bMR\nP9PvVv7O7FIAAGgQoRdBYcWBFUqNTVVsRKzLc9m70hseGq6EqATtO7mvyXEHTh1QVKsoRbaMVHKy\n/O5mNpuHxzysJfuWaP2h9WaXAgBAPYReBIXPdn2m6f2mu2Uue1d6JftaHDIKMpTWJU2SlJws7fbT\n3b/at2qvP1/xZ9276F5VVlWaXQ4AALUQehHwqqqrtCBrgW7se6Nb5rN3pVeS+sf21/aC7U2OySjI\n0KDOgyRJffr470qvJN024DZ1aduFAysAAD6H0IuAt/bQWnVu21lJMUkuz1VZKR07JnXvbt/4oV2H\nauvRrU2OSS9IV1pnY6W3d29p3z6putrVSs1hsVj0f9f8n55d82yzYR8AAG8i9CLg/Wvbv3Rr/1vd\nMtfhw1LnzlJYmH3jh3Ydqi35W5ocU3Olt21bKTpaOnTI1UrNkxidqL9c+Rfd/PHNOl1x2uxyAACQ\nROhFgCsuL9YnOz/R3UPudst8jvTzSlKvqF4qOVei42cb3nz37LmzOnT6kJI7JF94zt9bHCRpZtpM\nTe45WT/6z49k9bfTNgAAAYnQi4D23vb3dFniZerStotb5nOkn1cyvt0/uMvgRlscdhzbob4d+yos\n5OLSsb9uW1bXX6f+VUdKjujhpQ8TfAEApiP0ImBVVVfp5Q0v66fDf+q2OW1HEDuiqRaHbUe3Xdi5\nwSYQVnolqVVYK31x+xf6Jucb/X7l780uBwAQ5Ai9CFjvb39fMa1jNKXXFLfNefCgY+0NkjS0y1Bt\nzt/c4GtrD63VmLgxtZ4LlJVeSYpuHa2ldy7VvMx5emjJQ6qqrjK7JABAkCL0IiCdqzqnJ1c+qaem\nPCWLxeK2eZ1Z6R0VN0pr89Y2+C3+NXlrNDZ+bK3nAmWl16Zz285a86M12np0q6Z/NF3F5cVmlwQA\nCEKEXgSk3638nQZ0GqCJPSe6dV5nVnqToo2t0vYW7a31fH5JvgpLC5Uam1p7fJJxnfPnXSrVp0S3\njtaSO5corl2cBr02SN/s/8bskgAAQYbQi4Cz+uBq/WPrP/TGtDfcOq/V6txKr8Vi0eRek/VNTu2g\nt2TfEl2WeJlCLLX/GbZsKXXtKh044GLBPiY8NFyvXP2K3pj2hu7+z9269ZNbta+o6SOaAQBwF0Iv\nAkpGQYZu/vhm/fP6f6pz285unfvECSOQRkY6/t5pfabpP7v/U+u5xXsW66reVzU4PpD6euu6sveV\n2vWzXRrQaYBG/X2UZs6fqQ2HN5hdFgAgwBF6ETAWZS/S5e9erpeueklTe091+/zOrPLaXN3naq3J\nXaMTpSckGfsHL923VNcmX9vg+EDr660rIjxCj094XNm/yFZa5zTd+smt6vdKPz227DFtPLyRG94A\nAG5n57lSgG+yWq1af3i9nvvuOW06skmfzfhM43qM88i1nOnntYlsGanp/abrtU2v6bEJj+nt9Ld1\nWeJlio2IbXB8crK0e7cLxfqJmNYx+vXYX+uhMQ9p05FNmp81X7MWzNLhksMa2X2kRnYfqb4d+qpv\nx75K6Ziidi3bmV0yAMBPWawBtmu8xWLx6Y3wT5adVGFpoSTJqot11q25sdccfT7Q5jpVfkrHzh5T\nfkm+0gvStf7werVp0Uazh83WfcPvU5sWbeQpf/ubtHev9NJLzr1/1/FdmvCvCXpx6ov65ZJf6r8z\n/3vh+OG6vvxS+stfpKVLXSjYj50oPaG1h9Zq05FN2n1it7IKs5R9IluhllB1advlwiO6dbTahrdV\nZHjkhR8jwiPUIqSFWoS2uPBjWEhYrecsFossslzY2cP287o/NvaapFo/B4KNO3fF8YYWIS2UEOXk\nqoWH+XpuCSSEXi/7x9Z/6E/f/unCxzU/adb9T6Sx1xx9PpDmat+yvWIjYtWlbRcN6jRIw7oNU5+Y\nPl75D/ihh6Ru3aRf/9r5OT7O/Fh/Xf9X3T/8ft0x6I5Gx+3bJ116aeDdzOYKq9Wq0xWndfTM0QuP\nU+WnVHKuRGfOnbnw49lzZ1VZXanKqspGf7TKeuH/CdvP6/7Y2Gu2Wmp+cQYEE1/+HNuYxOhELZ3p\nm6sIvp5bAgmhF7DTzTdLM2YYD087f15q21Y6dUpq1crz1wMAmIPc4j3cyAbYyZWeXkeFhUk9exor\nvgAAwHWEXsBOruze4IxA38EBAABvIvQCdigrk4qLpc7u3fq3SYG8Vy8AAN5G6AXskJcnxcVJIV78\nF8NKLwAA7kPoBezgzX5emz59WOkFAMBdCL2AHbzdzysZ7Q2s9AIA4B6EXsAOZqz0du9ubFlWUuLd\n6wIAEIgIvYAdzFjpDQmRevc2ToEDAACuIfQCdjBjpVeirxcAAHch9AJ2MGOlV6KvFwAAdyH0As2o\nrpYOHZLi471/bbYtAwDAPQi9QDOOHpWio6VWrbx/bQ6oAADAPQi9QDPMam2QWOkFAMBdCL1AM8y6\niU2SOnWSKiuloiJzrg8AQKAg9ALNMHOl12JhtRcAAHcg9ALNMHOlV6KvFwAAdyD0As0wc6VXIvQC\nAOAOhF6gGWav9PbtK+3aZd71AQAIBIReoBlmr/Smpko7d5p3fQAAAoHFarVazS7CnSwWiwLslwQT\nnT4tdesmlZQYN5WZobxciooyamjRwpwaAACeQW7xHlZ6gSbYVnnNCryScShGfLy0d695NQAA4O8I\nvUATDh40t7XBpl8/+noBAHAFoRdowsGDUs+eZldhhF76egEAcB6hF2jCgQO+EXpTU1npBQDAFYRe\noAm+EnpZ6QUAwDWEXqAJvhR6d++WqqrMrgQAAP9E6AWa4CuhNzJS6tDB6DEGAACOI/QCjTh71tgb\nt3Nnsysx0NcLAIDzCL1AI2zHD5u5R29N9PUCAOA8Qi/QCF9pbbBhpRcAAOcReoFGHDhgrPT6ClZ6\nAQBwHqEXaISvHExhY1vp5Yh2AAAcR+gFGuFr7Q0dOkgtW0pHjphdCQAA/ofQCzTC10KvZKz2Zmaa\nXQUAAP6H0As0whdD76BB0vbtZlcBAID/IfQCDSgtlU6f9p09em3S0qT0dLOrAADA/xB6gQYcPCj1\n6CGF+Ni/kEGDCL0AADjDxz6lA77BF1sbJKl/fyk7Wzp3zuxKAADwL4ReoAE5Ob4Zetu0MfYOzsoy\nuxIAAPwLoRdowN69Uu/eZlfRsLQ0KSPD7CoAAPAvhF6gAfv2+Xbopa8XAADHEHqBBuzdKyUlmV1F\nw7iZDQAAx1ms1sA61NRisSjAfknwsupqKSJCKiw0fvQ1eXnSiBHS0aNmVwIAcBW5xXtY6QXqOHxY\nio72zcArSXFxUkWFVFBgdiUAAPgPQi9Qhy/fxCZJFgt9vQAAOIrQC9Thyzex2QwbJm3ebHYVAAD4\nD0IvUIcv38RmM2KEtHGj2VUAAOA/CL1AHb7e3iBJw4cTegEAcAShF6jDH0JvUpJ09iw7OAAAYC9C\nL1CD1Wr09Pp6e4PFYqz2btpkdiUAAPgHQi9Qw7FjUqtWUlSU2ZU0j75eAADsR+gFavCHm9hsCL0A\nANiP0AvU4A/9vDa20MtBPgAANI/QC9Swe7fUp4/ZVdine3epRQspN9fsSgAA8H2EXqCGrCypXz+z\nq7DfiBHShg1mVwEAgO8j9AI17NrlX6F33Djp22/NrgIAAN9H6AW+V1kp5eRIyclmV2K/Sy6RVq82\nuwoAAHwfoRf43t69Uny81LKl2ZXYb9gwKTtbKi42uxIAAHwboRf4nr+1NkhSeLjR1/vdd2ZXAgCA\nbyP0At/zx9ArSRMm0OIAAEBzCL3A93btkvr2NbsKx9HXCwBA8wi9wPf8daV39Ghp61apvNzsSgAA\n8F2EXkDS+fNG6B0wwOxKHNe2rZSayn69AAA0hdALyNgBoXt3I0D6o0mTpGXLzK4CAADfRegFJKWn\nS4MGmV2F8668UvrqK7OrAADAdxF6AUkZGVJamtlVOG/8eKM948QJsysBAMA3EXoBGSu9/hx6W7aU\nJk6U/vtfsysBAMA3EXoBGSu9/tzeIElTp0pffml2FQAA+CaL1Wq1ml2EO1ksFgXYLwkeduKElJgo\nnTolWSxmV+O8/fulsWOlI0ekEL6cBQC/QG7xHj41Iuht2SINHuzfgVcygnv79karBgAAqI3Qi6C3\nfr00apTZVbjHVVdJX3xhdhUAAPgeQi+C3vr10siRZlfhHjfeKH36qdlVAADgewi9CGpWq3GSWaCs\n9I4fL+XnS3v3ml0JAAC+hdCLoHbwoHHTV1yc2ZW4R2ioNH269MknZlcCAIBvIfQiqNlWef39Jraa\nbr5Z+ugjs6sAAMC3EHoR1NauDZx+XpuJE6Vjx6Tt282uBAAA30HoRVBbscIIiYEkNFT64Q+lt982\nu6IpxNYAAAmbSURBVBIAAHwHh1MgaBUVST17SoWFUni42dW41+7dRpg/dEgKCzO7GgBAY8gt3sNK\nL4LWypXGCWaBFnglKSVFSkqSFi40uxIAAHwDoRdBa9kyafJks6vwnF/8QnrpJbOrAADANxB6EbSW\nLw/s0HvTTVJ2tpSRYXYlAACYj9CLoHTggFRQIA0danYlntOihfTTn0ovvGB2JQAAmI9bXBCUFi6U\npk0L/Ju8fvYzqU8fad8+o8cXAIBgxUovgtKCBdINN5hdhedFR0v33y89/bTZlQAAYC62LEPQsW1V\ndvSo1KaN2dV4XlGRsZvDihVS//5mVwMAqInc4j2s9CLozJ8vXXZZcAReSYqJkR5/XPrVryT+XwUA\nBCtCL4LOW29Jd99tdhXedf/9xkEV8+ebXQkAAOagvQFBZdcu6dJLpdzcwL+Jra5vv5V+8ANp+3ap\nY0ezqwEASOQWb2KlF0HlzTelu+4KvsArSePHS7ffLt17L20OAIDgw0ovgkZhoXFD19atUo8eZldj\njooKacIEafp06f/9P7OrAQCQW7wnCNe7EKxeeMH49n6wBl5JatlS+vRTadQoqVcvacYMsysCAMA7\nCL0ICkeOSK+/Lm3ZYnYl5ouLk778Urr8cqlVK+m668yuCAAAzyP0IuBZrcZxvD//uZSQYHY1vmHQ\nIGnRIiPwHj0q3XOPZLGYXRUAAJ5DTy8C3r/+JT37rNHL27Kl2dX4luxso783LU169VWpfXuzKwKA\n4EJu8R52b/CyFStWmF1CUFm+3Lhh6+OPvRN4/e3PNzlZ2rjRCLsDB0rvvitVV5tdFQC4j7/9vwzP\nIfR6Gf/4vOeTT4wbtT780HvH7/rjn2/r1sYq7/vvGz8OHmxs7Xb2rNmVAYDr/PH/ZXgGoRcB59Ah\nadYs6eGHpSVLpMmTza7IP4wfL333nfTcc0a/b/fuxm4Xb78t7d3L3r4AAP/GjWzweyUlUlaWtGmT\ntHixtGaNcQDDjh1SZKTZ1fkXi0W68krjceyY9MUXRgB+/HGprEwaMkRKSpISE41Hp05Shw4XHy1a\nmP0rAACgYQF5IxsAAIC/CLAo5rMCbqWXvzgAAACoi55eAAAABDxCLwAAAAIeoRcAAAABj9ALAACA\ngBdwN7J5W1FRkf7+97+rurpa5eXlSklJ0YEDB1RYWKg///nPZpcHNzt79qweffRRpaSkaM+ePRo3\nbpxuvfVWl8cCAJrH51y4gtDrojfeeEOPPvqoQkJC1LlzZ7322mvq3LmzFi9ebHZp8IBHH31UycnJ\neuCBB1RVVaV+/fppwIABGjBggEtjAQDN43MuXBFw+/R6U1VVlbZt26Zhw4bp+PHj6tmzp4qLixUW\n5vjXEldddZW+/PJLD1QJdyktLVVsbKwyMzPVs2dPSdKDDz6oiIgIPf30006PBQA0z52fcxGc6Ol1\nQWhoqIYNGyZJWrVqlUaNGuX0Pz4Cr+/bvXu3ysrKlJCQcOG5+Ph4rVq1yqWxAIDmufNzLoIToddN\nVq1apXHjxl34eP369SZWA08oKChQy5Yta536FxERofz8fJfGAgAcw+dcOIMvkVzw0Ucf6amnntLm\nzZv1+eef68knn5QkzZ8/X8nJyVq7dq0+//xz9evXT23bttUXX3yh22+/XaWlpfr444/11ltvKS8v\nT998843OnDmjX/7yl/ruu++0aNEipaamKjIyUosXL9b//u//Kjo62txfLFRUVKRWrVrVei48PFyn\nT592aSwAoHnu+JwbFhamBQsW6NixY6qurtaePXv05z//WVarVR9++OGF58+dO6fExET94Ac/MPcX\nDbdipdcFcXFxGjhwoJ577jn985//1PLly/Xqq68qNDRU/fv3V0VFhWJiYiRJN954oyIiIlRVVaVp\n06apqKhIpaWlOnTokBISEnTu3DlJ0rlz5xQTEyOr1arrr79ebdq0UV5enpm/THwvMjKy3nOlpaX1\nwq2jYwEAzXP1c25ZWZn27t2rd955R/fee686d+6spKQkSVJ6erqmTZum0aNH6/+3cz+hsL1xHMc/\nxp8aRYrNSEoWajaGzKgpzZgFEykLKTtK2VlYsZ2FBUmUjd2U7VCzkSkUFmNDCYuRUEwWGKkZM/6d\n3+JX95dyf7iue26n92t3znme+n5Xz6en5zyXl5fq7OyU2+02s118A3Z6v8Dr9crr9f549vl8r777\n/X7NzMwoHA5Lkg4ODjQ9Pa3n52fd3d2ptLRULS0tGhgYUCgUenPO3t6enE7nH+oI/8fhcCidTr96\nl06nVV1d/aWxAID3fXXNLSkp0eTkpLq6uiRJ8XhcfX19kiSXyyVJ2trakt/v55Ydi2Kn9xu9vLzo\n9vZWpaWlymQyKioqks1m08bGhgKBgDY3N/Xw8KDLy0tdXV0pk8m8mnN8fKzq6mrF43GzW4Gk+vp6\n2e12XVxc/Hh3dHQkj8fzpbEAgK/7yJqbSqXkdDplGIbi8bjKysp0cnKi3d1dJZNJra+vy+PxKJFI\n8A+GBRF6v9Hp6amam5slSfv7+woEApKkgoIC5eXlKT8/X09PT6qqqlIymVRxcfGrOYZhqLy8XNls\n1rQe8J/CwkJ1d3draWlJkpTNZhWLxTQ4OChJCofDmpiY+NBYAMDv9d6aa7PZ1N/fr9XVVUUiEblc\nLsViMdXU1Gh5eVnRaFRNTU2KxWI6PDyUw+Ewsx18A+7pBT7h5uZGw8PDamxsVCKRUCAQUG9vryRp\nZGREZ2dnikQi744FAAB/FqEXAAAAlsfxBgAAAFgeoRcAAACWR+gFAACA5RF6AQAAYHmEXgAAAFge\noRcAAACWR+gFAACA5RF6AQAAYHkFZhcAAH+zmZkZZTIZ7ezsKBQKaXFxUYZhKJVKaWpqyuzyAAAf\nxE4vAPzE7OysgsGgxsbG1NbWJr/fr4GBAeVyOS0sLJhdHgDgEwi9APATT09PqqurkySdn5+roaFB\nlZWVGhoa0ubmpsnVAQA+I88wDMPsIgDgb9fa2qr29naNjo6aXQoA4Bew0wsA78jlctre3pbP5zO7\nFADALyL0AsAbHh8ftba2JkmKx+MyDENut1uSdH19rXA4bGZ5AIBPIvQCwBvm5+fV0dGh+/t7RaNR\nVVRUqKDg3wtv5ubm1NXVZXKFAIDP4EwvALxhf39f4+Pjqq2tVTAY1MrKitLptOx2u3p6euRyucwu\nEQDwCYReAAAAWB7HGwAAAGB5hF4AAABYHqEXAAAAlkfoBQAAgOURegEAAGB5hF4AAABYHqEXAAAA\nlkfoBQAAgOURegEAAGB5/wCT+t4GlHZ4CwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xa8ce52c>"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}