TISE-Numerov-var.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical solution to the 1-dimensional Time Independent Schroedinger Equation\n",
    "Based on the paper \"Matrix Numerov method for solving Schroedinger's equation\" by Mohandas Pillai, Joshua Goglio, and Thad G. Walker, _American Journal of Physics_ **80** (11), 1017 (2012).  [doi:10.1119/1.4748813](http://dx.doi.org/10.1119/1.4748813)\n",
    "\n",
    "## Variational method\n",
    "Using the Numerov approximation for the kinetic energy operator to calculate the expectation value of energy of a trial wavefunction.  From here, the energy may be minimized to find a variational solution.  Plot the ground state of this solution and the \"true\" solution, and compare their energies.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.379434Z",
     "start_time": "2019-01-18T21:41:04.370438Z"
    }
   },
   "outputs": [],
   "source": [
    "# import some needed libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.optimize as opt\n",
    "import scipy.integrate as integrate\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#import timeit # for timing the execution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.395425Z",
     "start_time": "2019-01-18T21:41:04.384430Z"
    }
   },
   "outputs": [],
   "source": [
    "autoscale = False # set this equal to true to use Pillai's recommended step sizes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is the finite square well.\n",
    "\n",
    "Using units such that &hbar;, *m*, and *a* are one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.429402Z",
     "start_time": "2019-01-18T21:41:04.399421Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# values of constants\n",
    "hbar = 1.0\n",
    "mass = 1.0 # changing the mass will also change the energy scale\n",
    "a = 1.0 # half width of the well (called a in the homework)\n",
    "V0 = 18.0*np.sqrt(hbar/mass)/a**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Exact solution\n",
    "Following McIntyre section 5.5, the exact energy must be found by solving a transcendental equation. We first define the dimensionless parameter\n",
    "\n",
    "$$\n",
    "z_0 = \\sqrt{\\frac{2mV_0a^2}{\\hbar^2}}.\n",
    "$$\n",
    "\n",
    "Then, we solve the equations\n",
    "$$\n",
    "z \\tan z = \\sqrt{z_0^2 -z^2}\\qquad\\text{and}\\qquad -z \\cot z = \\sqrt{z_0^2 -z^2},\n",
    "$$\n",
    "for all of the possible values of $z$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this problem, we only need the ground state, which will be given by the lowest $z$ value of the tangent equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.0\n"
     ]
    }
   ],
   "source": [
    "z0 = np.sqrt(2*mass*V0)*a/hbar\n",
    "print(z0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.344751045375789\n"
     ]
    }
   ],
   "source": [
    "#solve for z.  We know that the first solution must be in the range [0,π/2] by comparing to the \"graphical method\" of McIntyre.\n",
    "z_solution = opt.root_scalar(lambda z: z*np.tan(z)-np.sqrt(z0**2-z**2), bracket=(0,np.pi/2))\n",
    "\n",
    "z=z_solution.root\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know $z$, we can calculate energy and the exact wavefunction.\n",
    "\n",
    "$$\n",
    "z = \\sqrt{\\frac{2mEa^2}{\\hbar^2}} \\qquad \\rightarrow \\qquad E=\\frac{\\hbar^2}{2ma^2}\\,z^2.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9041776870196385\n"
     ]
    }
   ],
   "source": [
    "E_exact = (z*hbar/a)**2/2/mass\n",
    "print(E_exact)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need four parameters for the wavefunction:\n",
    "$$\n",
    "k = \\frac{z}{a},\\qquad q=\\frac{\\sqrt{z_0^2-z^2}}{a}. %k^2+q^2 = (z_0/a)^2\n",
    "$$\n",
    "Putting off normalization:\n",
    "$$ D = 1,\\qquad A=e^{qa}\\cos(ka) =e^{\\sqrt{z_0^2-z^2}} \\cos z $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=1.34475, q=5.84736, D=1.00000, A=77.619\n"
     ]
    }
   ],
   "source": [
    "k_exact = z/a\n",
    "q_exact = np.sqrt(z0**2-z**2)/a\n",
    "D_exact = 1.\n",
    "A_exact = np.exp(np.sqrt(z0**2-z**2))*np.cos(z)\n",
    "print(f\"k={k_exact:.5f}, q={q_exact:.5f}, D={D_exact:.5f}, A={A_exact:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def psi_exact(x):\n",
    "    # This is a long way to write the function, but it is clearer than the standard ways\n",
    "    # normalized (integrated by hand)\n",
    "    norm = 1/np.sqrt(a*(1+np.sin(2*z)/(2*z)+np.cos(z)**2/np.sqrt(z0**2-z**2)))\n",
    "    # The booleans are zero except when true\n",
    "    return norm* (((x<-1) * (A_exact*np.exp(q_exact*x))+\n",
    "           ((x>=-1) * (x<=1)) * (D_exact*np.cos(k_exact*x))+\n",
    "           (x>1) * (A_exact*np.exp(-q_exact*x))))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#check normalization.  This should be True.\n",
    "np.isclose(integrate.quad(lambda x: np.abs(psi_exact(x))**2, -10,10)[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fed39a34a00>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+j0lEQVR4nO3de3xU9Z3/8feZmWRyT0hCQgIhBJA7ogZUUFqxFaVda29K19+KbXF/IrVdpa2t9bdtdbfLbnfrsr1AdRW7F2rpzdW21Eq7XlDqhfstIPcASQghZCYJSSYzc35/TM5AIIFMksmZy+v5eMwDcnIm88kQZt75fj/n+zVM0zQFAABgE4fdBQAAgORGGAEAALYijAAAAFsRRgAAgK0IIwAAwFaEEQAAYCvCCAAAsBVhBAAA2MpldwF9EQwGVVNTo+zsbBmGYXc5AACgD0zTVHNzs0pLS+Vw9D7+ERdhpKamRmVlZXaXAQAA+uHYsWMaNWpUr5+PizCSnZ0tKfTN5OTk2FwNAADoC6/Xq7KysvD7eG/iIoxYUzM5OTmEEQAA4szlWixoYAUAALYijAAAAFsRRgAAgK0IIwAAwFaEEQAAYCvCCAAAsBVhBAAA2IowAgAAbEUYAQAAtiKMAAAAWxFGAACArQgjAADAVnGxUR6AxNPa4deW6jPacrRJWWkuXVeRr8klOXI6Lr2hFoDEQxgBMGQON7Tq+Xer9c7hRu064VEgaHb7fLbbpcoxwzRnXIH+z3XlynTzEgUkA8M0TfPyp9nL6/UqNzdXHo9HOTk5dpcDIEKmaernm47p2y/tUVtnIHx8ZF66Zo0ZJk9bpzYdOaPmDn/4c2OHZ+r7n7la00bm2lEygEHQ1/dvfu0AEFWetk5944Wd+t2OWknS9WPztXBWmWaNydeoYRnh8wJBU1W1Xr1zuFH//sYhHTrVqk+u3KivLZikz98wRobB9A2QqBgZARA1m4826kvPb9OJpja5HIa+PH+i7v/AWDku0xdyptWnR361Q+v3nJQkzZs4XP985wwVZrmHomwAg6Sv799cTQMgKt7c36CFT72tE01tKstP1y+WzNYDN427bBCRpGGZqXr6nkr93R1Tlepy6NV9p/SpVRvV2OobgsoBDDXCCIBBd7ihVUvXbJY/aOqWKcVa96W5unr0sIi+hmEYumf2GL304A0aNSxdR0+f1QP/vVk+fzBKVQOwC2EEwKDytHVq8X+8J2+7X1ePztMP/vJqZael9PvrTRqRo9WfnaXMVKfeOdyob720W3EwuwwgAoQRAIPGHwjqwZ9u0aFTrSrNTdNT91QqLcU54K87oThb3//Lq2UY0vPvVus/Nh4ZeLEAYgZhBMCg+c66Km3Y36D0FKf+/d6ZKspOG7Sv/aHJxXp0wSRJ0hO/3aM33j81aF8bgL0IIwAGxc83HdNzbx2RJP3rwhmaWjr464P89dyx+nTlKAVN6Qs/3aIjDa2D/hgAhh5hBMCA1Te36+9+s0eStOyWCbptWklUHscwDH3nE9M0s3yYmtv9+tsXd9E/AiQAwgiAAfvH3+9Vc4dfM0bl6sF546P6WG6XU9+7a4ZSnQ5t2N+gP+yui+rjAYg+wgiAAdl0pFG/3nJChiE9fse0Pq0jMlDlBZm6/4NjJUl/99sqtfkCl7kHgFhGGAHQb4Ggqb99cbckaeHMMl1Vljdkj730pvEamZeuE01tWvnagSF7XACDjzACoN/WvHNUVbVe5aan6JHbJg3pY6enOvW3fzFFkvTU64doZgXiGGEEQL+cbunQv/xhnyTpK7dOVH5m6pDXcOvUYs29olC+QFCP/4bF0IB4RRgB0C//9PJeedv9mlqao7uvHW1LDYZh6PGPTVWK09Cr+07pj1X1ttQBYGAIIwAitrvGo59vOi5JeuKOaXIOQdNqb8YOz9J9c61m1j3yB9i7Bog3hBEAEVv12kFJ0u0zSlVZHtkGeNHw4Lzxys9MVXXjWf1uZ63d5QCIEGEEQESONLRqXdcb/tKbxtlcTUim26XPzRkjKRSU6B0B4gthBEBEnt5wSEFTmjdxuCaX5NhdTtg9s8uVkerU3rpmvca+NUBcIYwA6LN6b7t+2dUr8sBN0V1pNVJ5GanhRlprGglAfCCMAOiz1W8dkS8QVGX5MM0aY3+vyIUWz61QitPQu4cbtfloo93lAOgjwgiAPvG2d2rN20clSQ98cJwMw74raHpTkpuuT1w9UpK06rVDNlcDoK8IIwD65L/fPqrmDr8mFGfp5klFdpfTq//7gXEyDOmPVSf1/slmu8sB0AeEEQCX1d4Z0Oo3j0iSlnxw3JBshtdf44uydOuUEZKkH79O7wgQDwgjAC7rV1uOq6GlQyPz0nX7jFK7y7msJV2XHL+0rUYnmtpsrgbA5RBGAFySaZp67q0jkqTFN1YoxRn7LxtXleVpzrgC+YOm/rurzwVA7Ir9VxUAttpS3aQD9S1KT3Hqzpmj7C6nzxbNLpck/XLzcZaIB2IcYQTAJa19r1qS9JHpJcpOS7G5mr67eVKxCjJTdaq5Q6/uYxE0IJYRRgD0qqXDr9/uCC39vnBWmc3VRCbV5dAnrwld5rv2vWM2VwPgUggjAHr1ux01OusLaGxhZkwucnY5VoB6dV+96r3tNlcDoDeEEQC9skYU7pxZFpOLnF3O+KJsXTM6T4GgqV9tOWF3OQB6QRgB0KP9J5u1pbpJToehT1WOtLucfvvMrNB+NT/fdIzdfIEYRRgB0CNrVOTmSUUqyk6zuZr+++iVJcpMdepwQ6vePcx+NUAsIowAuIjPH9Svt4amNRbOjK/G1Qtlul36iytDC7Wt3UQjKxCLCCMALvKnqpNqbPWpKNutmyYOt7ucAburq5F13c5aeds7ba4GwIUIIwAuYo0gfLpylFxxsOLq5VwzOk9XFGWpvTOo32yvsbscABeI/1cZAIOq1tOm198PLRJ2V5xP0VgMwwhf5suaI0DsIYwA6OZ3O2plmtKsMcM0pjDT7nIGzSeuHimnw9CO4x4dPd1qdzkAzkMYAdDN73aGVly1mj4TRUGWW7PHFkg69z0CiA39CiMrV65URUWF0tLSVFlZqQ0bNlzy/DVr1mjGjBnKyMhQSUmJPve5z+n06dP9KhhA9NQ0tWlrdZMMQ1owbYTd5Qy6j0wvkRRqZAUQOyIOI2vXrtVDDz2kxx57TFu3btXcuXO1YMECVVdX93j+m2++qUWLFmnx4sXavXu3fvGLX+i9997TfffdN+DiAQwu6016Vnm+inLid22R3tw6tVhOh6FdJ7yqPn3W7nIAdIk4jDz55JNavHix7rvvPk2ePFkrVqxQWVmZVq1a1eP5b7/9tsaMGaMvfelLqqio0I033qj7779fmzZtGnDxAAaXFUY+Mj3xRkWk0FTN9WPzJTFVA8SSiMKIz+fT5s2bNX/+/G7H58+fr40bN/Z4nzlz5uj48eNat26dTNPUyZMn9ctf/lIf/ehH+181gEFX09SmLdYUTdd0RiJiqgaIPRGFkYaGBgUCARUXF3c7XlxcrLq6uh7vM2fOHK1Zs0YLFy5UamqqRowYoby8PP3gBz/o9XE6Ojrk9Xq73QBE1+93hf4PzyrPV3ECTtFYbps6Qg5D2nnCw1QNECP61cB64e6dpmn2uqPnnj179KUvfUnf/OY3tXnzZr388ss6fPiwlixZ0uvXX758uXJzc8O3srLEWOsAiGW/2xFaDCxRp2gsBVluzR7HVTVALIkojBQWFsrpdF40ClJfX3/RaIll+fLluuGGG/TVr35VV155pW699VatXLlSq1evVm1tzy8Ejz76qDweT/h27BiLFAHRlCxTNBamaoDYElEYSU1NVWVlpdavX9/t+Pr16zVnzpwe73P27Fk5HN0fxul0SlKv23m73W7l5OR0uwGIHmuKZmb5sISeorHcylQNEFMinqZZtmyZnnnmGa1evVpVVVV6+OGHVV1dHZ52efTRR7Vo0aLw+bfffrt+/etfa9WqVTp06JDeeustfelLX9K1116r0tLEWlQJiFfnrqJJ/FERSSrMcuv6rgXQ1u1idASwmyvSOyxcuFCnT5/WE088odraWk2bNk3r1q1TeXm5JKm2trbbmiOf/exn1dzcrB/+8If68pe/rLy8PN188836p3/6p8H7LgD0W62nTZuPnpEkLZiWHGFECgWvjQdPa93OWi354Di7ywGSmmH2NlcSQ7xer3Jzc+XxeJiyAQbZ6jcP64nf7tHM8mH65QM9T7cmooaWDl37nT8qaEobHpmnsvwMu0sCEk5f37/ZmwZIci/vDvWLJMsUjaUwy63rKkJTNS/v6nlpAgBDgzACJDHP2c7wFM38qT1fEZfIrO/5f/fW21wJkNwII0ASe33/KQWCpiYUZ2nUsOSbprh5UpEk6b0jjfK2d9pcDZC8CCNAEnu1a0RgXtebcrIpL8jU2OGZ8gdNvbm/we5ygKRFGAGSVCBo6rV9XWFkYnKGEenc985UDWAfwgiQpLYda9KZs53KTnOpsnyY3eXYxpqqeW3fKQWDMX9xIZCQCCNAkrKmaD4wYbhSnMn7UjBrTL6y3C41tHRoV43H7nKApJS8r0BAknu1a4rm5iSeopGkVJdDN44vlMRUDWAXwgiQhE5627W7xivDkG6aONzucmxnTdW8ShgBbEEYAZKQ9aY7Y1SeCrLcNldjv5smhQLZ9uMenWrusLkaIPkQRoAkZE1H3Jykl/ReqCg7TdNH5kqSXn//lM3VAMmHMAIkmQ5/QG8eCK2pQRg5Zx5TNYBtCCNAknn3cKPO+gIanu3WlBI2nrTM6+qdeeP9U+oMBG2uBkguhBEgyVhTNPMmDpfDYdhcTeyYMSpPBZmpau7wa9ORM3aXAyQVwgiQZF6lX6RHDoehD3aNjliXPQMYGoQRIIkcbmjVkdNnleI0dOMVXNJ7ISugsd4IMLQII0AS2bA/dKXIzPLQqqPobu4Vw+UwpAP1LarztNtdDpA0CCNAEnmr6yqaG68otLmS2JSbnqLpo/IknXuuAEQfYQRIEoGgqT8fPC1JmjOuwOZqYtcNXc/NWwcJI8BQIYwASWLXCY+87X5lu13hBb5wsRu69ql560CDTJNdfIGhQBgBkoT1m/714wrkSuJdei+nsnyY3C6HTno7dPBUq93lAEmBVyQgSWw8EJqiuYEpmktKS3Fq5phhkqSNTNUAQ4IwAiSB9s6A3jvSKOncNAR6N2fcuakaANFHGAGSwJajZ9ThD6oo263xRVl2lxPzrMD254OnFQjSNwJEG2EESAJWv8iccQUyDJaAv5zpI3OVneaSt92vXSc8dpcDJDzCCJAE3rL6RZii6ROnw9DssVziCwwVwgiQ4LztndpxvEkSYSQS1nNlNf4CiB7CCJDg3j54WkFTqijMVGleut3lxI0bxodGRt470qj2zoDN1QCJjTACJLiNrLraL+OGZ6ko260Of1Bbjp6xuxwgoRFGgAQX3o+GKZqIGIYRfs7oGwGiizACJLCT3nbtr2+RYUizGRmJ2JyuMPImfSNAVBFGgARmrSA6tTRHeRmpNlcTf6y+kZ3Hm+Rp67S5GiBxEUaABBa+pHccUzT9UZKbrrGFmQqa0juHGB0BooUwAiSwt7veQJmi6b85XaMjbx9qtLkSIHERRoAEVdPUpuNn2uQwpJlj8u0uJ25dWxEKI+8eYWQEiBbCCJCgrI3xppbmKsvtsrma+HVtV5DbU+NVczt9I0A0EEaABPXu4VAYubaCUZGBGJGbptH5GQqa0pbqJrvLARISYQRIUFYYmcUUzYBZz+G7h5mqAaKBMAIkoDOtPu2vb5EkzRozzOZq4t+1FaHn8L3DrMQKRANhBEhAVr/I+KIsFWS5ba4m/llNrNuONbFPDRAFhBEgAdEvMrjGFGSoMMstXyCoHcc9dpcDJBzCCJCArJGRa+kXGRSGYei6CvpGgGghjAAJprXDr101XknSLEZGBo3Ve/PuEfpGgMFGGAESzJbqMwoETY3MS9fIvHS7y0kYVrDbcvSM/IGgzdUAiYUwAiSY9+gXiYpJI3KUneZSS4dfVbXNdpcDJBTCCJBg3iGMRIXTYWhmuTVVwz41wGAijAAJpMMf0LZjTZJY7CwawvvU0MQKDCrCCJBAdp3wqMMfVEFmqsYNz7S7nIQTXvzsyBmZpmlzNUDiIIwACeSd85aANwzD5moSz/SReXK7HGps9engqRa7ywESBmEESCA0r0ZXqsuhq0fnSZLeZWl4YNAQRoAEEQia2tS1BgZhJHqsheTeo4kVGDSEESBB7K3zqrnDryy3S5NLcuwuJ2HNCq/EShgBBgthBEgQW46GRkWuHp0np4N+kWi5ZvQwOQzpRFOb6jztdpcDJATCCJAgtlQ3SZIqu9bCQHRkul2aNCI08rSlmr4RYDAQRoAEYb0xXj2aMBJtVhOrNRoFYGAII0ACaGjp0NHTZyVJV5Xl2VtMErimK/AxMgIMDsIIkAC2dk3RXFGUpdz0FHuLSQLXdE2F7TrhVYc/YHM1QPwjjAAJYGvXb+jXMEUzJMYUZCg/M1W+QFB7arx2lwPEPcIIkACs6YJryvPsLSRJGIahq7umw6zGYQD9RxgB4pw/ENT2Yx5JNK8OJWuqhr4RYOAII0Cc21vXrLbOgLLTXBo/PMvucpKGdUXNVq6oAQaMMALEOatf5KqyPDlY7GzIzBiVJ4ch1XjaWfwMGKB+hZGVK1eqoqJCaWlpqqys1IYNGy55fkdHhx577DGVl5fL7XZr3LhxWr16db8KBtCddSUNzatDK9Pt0sSuxc+2MlUDDEjEYWTt2rV66KGH9Nhjj2nr1q2aO3euFixYoOrq6l7vc9ddd+lPf/qTnn32We3bt0/PP/+8Jk2aNKDCAYSca14ljAy1a6zFzwgjwIC4Ir3Dk08+qcWLF+u+++6TJK1YsUJ/+MMftGrVKi1fvvyi819++WW9/vrrOnTokPLzQxtMjRkzZmBVA5AknW7p0BEWO7PNNaOHac071VxRAwxQRCMjPp9Pmzdv1vz587sdnz9/vjZu3NjjfV566SXNnDlT3/3udzVy5EhNmDBBX/nKV9TW1tbr43R0dMjr9Xa7AbiYNUUznsXObGGNRu084ZHPH7S5GiB+RTQy0tDQoEAgoOLi4m7Hi4uLVVdX1+N9Dh06pDfffFNpaWl64YUX1NDQoKVLl6qxsbHXvpHly5fr8ccfj6Q0ICmFp2i6pgswtKzFzxpbfdpd4+HSaqCf+tXAahjdO/ZN07zomCUYDMowDK1Zs0bXXnutPvKRj+jJJ5/UT37yk15HRx599FF5PJ7w7dixY/0pE0h4W1h51VYsfgYMjojCSGFhoZxO50WjIPX19ReNllhKSko0cuRI5ebmho9NnjxZpmnq+PHjPd7H7XYrJyen2w1Ad/5AUDuOhxY7o3nVPtZzzxU1QP9FFEZSU1NVWVmp9evXdzu+fv16zZkzp8f73HDDDaqpqVFLS0v42Pvvvy+Hw6FRo0b1o2QAkrTvZLPO+gLKdrPYmZ2skZGtjIwA/RbxNM2yZcv0zDPPaPXq1aqqqtLDDz+s6upqLVmyRFJoimXRokXh8++++24VFBToc5/7nPbs2aM33nhDX/3qV/X5z39e6enpg/edAEnGmha4ajSLndlpRllo8bMTTW066WXxM6A/Ir60d+HChTp9+rSeeOIJ1dbWatq0aVq3bp3Ky8slSbW1td3WHMnKytL69ev1xS9+UTNnzlRBQYHuuusu/f3f//3gfRdAErKWIadp0l7W4mdVtV5tOXpGC6aX2F0SEHciDiOStHTpUi1durTHz/3kJz+56NikSZMumtoBMDDbjjVJOrdHCuxzzeg8VdV6tfVYE2EE6Af2pgHikKetU4caWiVJV43Ks7cYaEZX38j2roAIIDKEESAO7ey6imZ0foaGZabaXA1mdAXCXSc8CgRNe4sB4hBhBIhD2483SZKuHJV76RMxJMYXZSkj1alWX0CHTrVc/g4AuiGMAHFoR1cYmcEUTUxwOgxNKw0Fw+1do1YA+o4wAsQha7GzGWyOFzOsUSorKALoO8IIEGfqve2q9bTLYUjTRrI6cay4kiZWoN8II0CcsaYBrijKVkZqv67ORxRYVzVV1Tazgy8QIcIIEGd20Lwak8ry0zUsI0W+QFB767x2lwPEFcIIEGeskZEr6ReJKYZhaHrX6AhNrEBkCCNAHDFNMzwywmJnsWdG12gVfSNAZAgjQBypbjyrprOdSnU6NHFEtt3l4AJXdgVErqgBIkMYAeKINfw/uTRHqS7++8Yaa2TkQH2LWjv8NlcDxA9ezYA4sqNr+H8GzasxqSgnTSW5aQqaoaXhAfQNYQSII+HFzugXiVnnFj8jjAB9RRgB4oQ/ENTOE9bKq4yMxCqrb2QbfSNAnxFGgDhx4FSL2joDynK7NLYwy+5y0IsZNLECESOMAHFix7HQqMi0kTlyOAybq0FvpndN0xxrbFNjq8/maoD4QBgB4sR2a6deFjuLabnpKRpbmCmJ0RGgrwgjQJygeTV+0MQKRIYwAsSB9s6AqmpD+52wJ03ss5pYWYkV6BvCCBAHqmq98gdNFWSmamReut3l4DKsq522H/fINE2bqwFiH2EEiAPWAlrTR+XKMGhejXVTS3PldBhqaOnQSW+H3eUAMY8wAsSBXSdCUzTTRzJFEw/SUpwaPzx0+TUrsQKXRxgB4oC12Nk0wkjcsP6tdhJGgMsijAAxrsMf0PsnmyURRuLJtJE5kqTdNYQR4HIII0CMe7+uRf6gqWEZKSrNTbO7HPSRFRytKTYAvSOMADHu/Ckamlfjx5SSHBmGVOdt16lmmliBSyGMADFuVw39IvEo0+0Kr8S6i6ka4JIII0CM222NjJQSRuKNFSB308QKXBJhBIhhnYGgquqs5tUcm6tBpKZzRQ3QJ4QRIIbtP9kinz+o7DSXRudn2F0OIjS1lCZWoC8II0AMC/eLlNK8Go+mdo1mnWhq05lWn83VALGLMALEsF3hK2mYoolHOWkpGlMQGtGiiRXoHWEEiGG7WHk17k1lvRHgsggjQIzyB4LaUxt6AyOMxK9p4b4RRkaA3hBGgBh1qKFV7Z1BZaY6VVGQaXc56CfrihqmaYDeEUaAGGX9Jj2lNEcOB82r8Wpqaajf5+jps/K0ddpcDRCbCCNAjGKn3sQwLDNVI/PSJbFpHtAbwggQo3Z3NTyy8mr8mx5eiZUmVqAnhBEgBgWDZvi3aEZG4p91aTZ9I0DPCCNADDp8ulWtvoDSUhwaN5zm1Xg3lWXhgUsijAAxyGpenVySI5eT/6bxzppqO9zQqpYOv83VALGHVzkgBu2uoV8kkQzPdmtETppMU6qqpW8EuBBhBIhBO4+HRkam0y+SMKzeH+vfFsA5hBEgxpimGV55dUope9IkCmu9kT2MjAAXIYwAMabG0y5PW6dcDkNXFGfZXQ4GiRUs99QQRoALEUaAGGO9WY0vypLb5bS5GgyWKSWhMLK/vlk+f9DmaoDYQhgBYowVRpiiSSyjhqUrO82lzoCpA/UtdpcDxBTCCBBjrMXOrN+kkRgMwwj/m7IsPNAdYQSIMVaD41Qu60041r8pTaxAd4QRIIZ42jp1/EybJEZGEhFNrEDPCCNADLEWxBqZl67cjBSbq8FgswLmnlqvTNO0uRogdhBGgBhC82piG1+UpRSnoeZ2f3gEDABhBIgp4cXOmKJJSKkuh64oypZE3whwPsIIEEMYGUl89I0AFyOMADHC5w9qf32zJEZGEtn5fSMAQggjQIw4UN+izoCp7DSXRg1Lt7scRAkjI8DFCCNAjDi/X8QwDJurQbRM7hoZOdHUJs/ZTpurAWIDYQSIEfSLJIfc9JTwyBdTNUAIYQSIEXtqWQY+WdA3AnRHGAFigGmajIwkEfpGgO4II0AMONHUJm+7XylOI7wOBRIXIyNAd/0KIytXrlRFRYXS0tJUWVmpDRs29Ol+b731llwul6666qr+PCyQsKzfkMcXZSvVxe8Iic4aGTlQ3yyfP2hzNYD9In7VW7t2rR566CE99thj2rp1q+bOnasFCxaourr6kvfzeDxatGiRPvShD/W7WCBRsfJqchmZl66cNJc6A2Z4bRkgmUUcRp588kktXrxY9913nyZPnqwVK1aorKxMq1atuuT97r//ft19992aPXt2v4sFEtVu+kWSimEY4X/r3fSNAJGFEZ/Pp82bN2v+/Pndjs+fP18bN27s9X7PPfecDh48qG9961t9epyOjg55vd5uNyCRhZtXGRlJGlNKciXRxApIEYaRhoYGBQIBFRcXdzteXFysurq6Hu+zf/9+ff3rX9eaNWvkcrn69DjLly9Xbm5u+FZWVhZJmUBc8Zzt1Imm0A6ujIwkj6mlNLECln51yl24OqRpmj2uGBkIBHT33Xfr8ccf14QJE/r89R999FF5PJ7w7dixY/0pE4gL1pvRqGHpyk1PsbkaDBUreFbVeGWaps3VAPbq21BFl8LCQjmdzotGQerr6y8aLZGk5uZmbdq0SVu3btWDDz4oSQoGgzJNUy6XS6+88opuvvnmi+7ndrvldrsjKQ2IW1VdYWTSCEZFksm44VlKcRpq7vDr+Jk2leVn2F0SYJuIRkZSU1NVWVmp9evXdzu+fv16zZkz56Lzc3JytHPnTm3bti18W7JkiSZOnKht27bpuuuuG1j1QALYW2f1i7C+SDJJdTk0bniWJGlvHVfUILlFNDIiScuWLdM999yjmTNnavbs2Xr66adVXV2tJUuWSApNsZw4cUL/+Z//KYfDoWnTpnW7f1FRkdLS0i46DiQr641oMs2rSWdKSY721jVrb61Xt0y5eHQZSBYRh5GFCxfq9OnTeuKJJ1RbW6tp06Zp3bp1Ki8vlyTV1tZeds0RACGBoKl9XWFkEmEk6UwqyZa2MjICGGYcdE55vV7l5ubK4/EoJ4cXbCSOA/Ut+vCTrys9xaldj98qp+PiRnAkrjfeP6VFq9/V2MJM/e9XbrK7HGDQ9fX9m3WnARtZ/SITRmQTRJLQpK4+ocOnW9XmC9hcDWAfwghgo721oeF5mleTU1F2mgqzUmWa0vsnmapB8iKMADbisl5Y//bWKBmQjAgjgI2sxsVJIxgZSVbWv31VLSMjSF6EEcAmnrZzy8AzMpK8rKuoqlgWHkmMMALYxLqkd2ReunIzWAY+WU3u6hfaW9fMsvBIWoQRwCbn+kWYoklm44uy5HQY8rR1qs7bbnc5gC0II4BNrIbFSVxJk9TcLqfGDc+UdO7qKiDZEEYAm1gNi/SLwPoZ2EPfCJIUYQSwQfC8ZeDZkwbWzwDLwiNZEUYAGxxtPKu2zoDcLofGFLB1fLKzpur2MjKCJEUYAWxgvelMKM6Wy8l/w2Q3uWua5lBDq9o7WRYeyYdXQcAGVSx2hvMU57iVl5GiQNDUgfoWu8sBhhxhBLCBdVkv/SKQJMMwwqMjLH6GZEQYAWzAZb240KTzFj8Dkg1hBBhize2dOtbIMvDobjIb5iGJEUaAIWZtFV+c41Z+ZqrN1SBWWCMjVbUsC4/kQxgBhtieWtYXwcUmFGfLYUiNrT6dau6wuxxgSBFGgCG2N7wnDWEE56SlOFVRGFoWvoq+ESQZwggwxPaGV16leRXdTSrhihokJ8IIMIRM89wy8BNZYwQXmNz1M7GPkREkGcIIMISOn2lTS4dfLoehsYVZdpeDGDOhmMt7kZwII8AQsq6kGTc8S6ku/vuhO6uP6GB9izoDQZurAYYOr4bAENrLFA0uYdSwdGWkOuULBHX0dKvd5QBDhjACDCH6RXApDofBVA2SEmEEGEL72CAPlzGJJlYkIcIIMER8/qAOngrtyMrICHpj/WwwMoJkQhgBhsihhhb5g6ay3S6NzEu3uxzEqImMjCAJEUaAIWK9uUwYkS3DMGyuBrFqYlfPSHXjWbV2+G2uBhgahBFgiNC8ir4oyHKrMMst6dyl4ECiI4wAQ4TmVfQVTaxINoQRYIhYDYnWpZtAb2hiRbIhjABDoLm9Uyea2iQxMoLLs8II0zRIFoQRYAhYbyrFOW7lZaTaXA1iHdM0SDaEEWAInFsGPsfmShAPrijKlmFIp1t9OtXcYXc5QNQRRoAhQPMqIpGe6lR5foYkRkeQHAgjwBAIX9ZL8yr66FwTq9fmSoDoI4wAUWaapvadZI0RRMaa0mNkBMmAMAJEWX1zh5rOdsrpMDS+KMvuchAnwk2sXFGDJEAYAaLMal4dU5ChtBSnzdUgXpx/eW8gaNpcDRBdhBEgyvZ1zflP4koaRGBMQabcLofaO4OqbjxrdzlAVBFGgChj5VX0h9Nh6Iri0LTePppYkeAII0CUsUEe+mtisdXE2mJzJUB0EUaAKAoETe2vD72RsMYIIjVxRNfIyElGRpDYCCNAFB053SqfP6j0FKdGdy1iBfSVdXkvG+Yh0RFGgCjaF+4XyZLDYdhcDeKNNZp2pKFV7Z0Bm6sBoocwAkTRXvpFMABF2W7lZaQoaEoH6ukbQeIijABR9D4b5GEADMMIbyHASqxIZIQRIIrCy8BzWS/6aSIrsSIJEEaAKGnzBXTkdKskpmnQf+c2zCOMIHERRoAo2V/fLNOUCjJTNTzbbXc5iFPhPWpY+AwJjDACRAmLnWEwWCv3nvR2qOmsz+ZqgOggjABRso9l4DEIstNSNDIvXRJNrEhchBEgSqyGQ1ZexUDRxIpERxgBooQ1RjBYaGJFoiOMAFHQ2OrTqeYOSdIVTNNggKzRtfcJI0hQhBEgCqy5/bL8dGW5XTZXg3hn9R3tO9ks0zRtrgYYfIQRIAqsyzCtLeCBgRg3PEsuh6Hmdr9qPO12lwMMOsIIEAX7Tob2EaF5FYMh1eXQ2OGZkpiqQWIijABRYI2MTCCMYJBY+xvRxIpERBgBBplpmnqfkREMsonFWZJYiRWJiTACDLLjZ9rU0uFXitNQRWGm3eUgQTAygkTWrzCycuVKVVRUKC0tTZWVldqwYUOv5/7617/WLbfcouHDhysnJ0ezZ8/WH/7wh34XDMS697sWpho3PEspTvI+Boc1ynboVKs6A0GbqwEGV8SvlGvXrtVDDz2kxx57TFu3btXcuXO1YMECVVdX93j+G2+8oVtuuUXr1q3T5s2bNW/ePN1+++3aunXrgIsHYhGLnSEaRualKzPVKV8gqCMNrXaXAwyqiMPIk08+qcWLF+u+++7T5MmTtWLFCpWVlWnVqlU9nr9ixQo98sgjmjVrlq644gr9wz/8g6644gr95je/GXDxQCxigzxEg8NhhBfQY6oGiSaiMOLz+bR582bNnz+/2/H58+dr48aNffoawWBQzc3Nys/P7/Wcjo4Oeb3ebjcgXljTNBNZeRWDLLwSK3vUIMFEFEYaGhoUCARUXFzc7XhxcbHq6ur69DW+973vqbW1VXfddVev5yxfvly5ubnhW1lZWSRlArbpDAR18FToShpGRjDY2KMGiapf3XWGYXT72DTNi4715Pnnn9e3v/1trV27VkVFRb2e9+ijj8rj8YRvx44d60+ZwJALNReayna7wtu+A4MlvHsvYQQJJqJNMwoLC+V0Oi8aBamvr79otORCa9eu1eLFi/WLX/xCH/7why95rtvtltvtjqQ0ICbsPW+xs74EdCAS1tRfdeNZtXb4lcm+R0gQEY2MpKamqrKyUuvXr+92fP369ZozZ06v93v++ef12c9+Vj/96U/10Y9+tH+VAnHAmsufQL8IoqAgy63CrNAvavvrW2yuBhg8EU/TLFu2TM8884xWr16tqqoqPfzww6qurtaSJUskhaZYFi1aFD7/+eef16JFi/S9731P119/verq6lRXVyePxzN43wUQI6zhc1ZeRbRMCk/V0NiPxBFxGFm4cKFWrFihJ554QldddZXeeOMNrVu3TuXl5ZKk2trabmuOPPXUU/L7/frCF76gkpKS8O1v/uZvBu+7AGLEnpqu3XoJI4gSK4zsriGMIHH0a8Jx6dKlWrp0aY+f+8lPftLt49dee60/DwHEnZqmNtV42uV0GJo+MtfucpCgrhqdJ0naUn3G3kKAQcRa1cAgsd4cJpdk01iIqKksHyZJqqptVmuH3+ZqgMFBGAEGyeajoTBSOXqYzZUgkZXkpqs0N02BoKntx5vsLgcYFIQRYJBs6Qoj15QTRhBd1s+Y9TMHxDvCCDAI2nyBcEPhzDG9b3UADIaZXWFkE2EECYIwAgyC7ceb5A+aGpGTptLcNLvLQYKrLA8F3i1HzygYNG2uBhg4wggwCML9IuXDWHkVUTepJFvpKU552/3hvZCAeEYYAQYB/SIYSilOh2aUhS4f38xUDRIAYQQYoGDQ1Oauy3pnEkYwRGZ2TdXQN4JEQBgBBuhQQ6uaznYqLcWhKaU5dpeDJFHJFTVIIIQRYICsN4MrR+Upxcl/KQyNq7tWYj3U0KrGVp+9xQADxCsnMEDnN68CQyUvI1Xji7IkMTqC+EcYAQbI6hdh5VUMNetnbjP71CDOEUaAAWg669OB+tCllVxJg6FWOaYrjBwhjCC+EUaAAbA2xxs7PFP5mak2V4NkY00Nbj/eJJ8/aHM1QP8RRoABYHM82GlsYabyMlLU4Q9qT63X7nKAfiOMAANA8yrsZBjGub4RmlgRxwgjQD91BoLafswjSZo5hjACe1h9I1xRg3hGGAH6qarWq7bOgHLTUzS2MMvucpCkrJGRTUcbZZpsmof4RBgB+mnD/gZJoSkah4PN8WCPK0flKdXp0Elvh/bXs2ke4hNhBOinl3fVSZJumVJscyVIZumpTt14RaGkcz+TQLwhjAD9cKzxrHae8MhhSPMJI7DZbdNGSJJ+TxhBnCKMAP3wh92hF/1rK/JVkOW2uRoku1smF8vpMFRV69XR0612lwNEjDAC9IM1HH7b1BE2VwJIwzJTdf3YfElM1SA+EUaACNV728N7gdw2rcTmaoAQ62eRqRrEI8IIEKE/7K6TaYa2cB+Rm2Z3OYAk6dYpxTIMaduxJtV62uwuB4gIYQSIkPWb54JpTNEgdhTlpIXXHGGqBvGGMAJEoLHVp3cON0qSbpvKFA1ii3VVDWEE8YYwAkTgj3tOKhA0NaUkR6MLMuwuB+jGCiPvHWlUQ0uHzdUAfUcYASLw+121kpiiQWwaNSxD00fmKmhKr+w+aXc5QJ8RRoA+8rZ36q0DpyVJC6YTRhCbwlM1u5mqQfwgjAB99OreevkCQY0vytL4omy7ywF6ZI3abTzQIM/ZTpurAfqGMAL00e93stAZYt/Y4VmaUJwlf9DUH6uYqkF8IIwAfVDnadef9oZe2D8ynatoENs+Or1UkrTmnaM2VwL0DWEE6IPn3jqszoCpayvyNaU0x+5ygEv6y+vKlOp0aEt1kzYdabS7HOCyCCPAZXjbO7XmnWpJ0pIPjrW5GuDyirLT9KnKkZKkH79+yOZqgMsjjACX8fw71Wrp8OuKoizdNKHI7nKAPrlv7lgZhvTHqpM6UN9idznAJRFGgEvw+YNa/dZhSdJff2CsHA7D5oqAvhk3PEsfnlwsSfr3NxgdQWwjjACX8OK2Ezrp7VBxjlt3XFVqdzlARKxpxRe2nlC9t93maoDeEUaAXgSDpp7u+o3yczdUyO1y2lwREJnK8nxVlg+TLxDUcxuP2F0O0CvCCNCL196v1/76FmW5Xbr7utF2lwP0y/0fCI2O/PfbR9XS4be5GqBnhBGgF9ZVCHdfN1o5aSk2VwP0z4cnF2vs8Ew1t/v1s3er7S4H6BFhBOjBpiONevdwo1Kchj53wxi7ywH6zeEw9H/nhkZHnn3zsNo7AzZXBFyMMAJcwOcP6rEXdkmSPnn1KJXkpttcETAwH796pEbkpKnW065/Xf++3eUAFyGMABdY+doB7TvZrILMVH1twSS7ywEGLC3Fqb//+DRJ0r9vOKQdx5vsLQi4AGEEOM++umb96NUDkqRvf2yq8jNTba4IGBwfnlKs22eUKmhKj/xyh3z+oN0lAWGEEaCLPxDUI7/crs6AqQ9PLtZfXMmGeEgs3759ivIzU7W3rlmrXjtodzlAGGEE6LL6rcPaftyj7DSXvvOJaTIMVltFYinIcutbt0+RJP3w1f3aV9dsc0VACGEEkHS4oVXfeyXU2Pf/PjpZxTlpNlcERMfHZpTqw5OL1Bkw9civdigQNO0uCSCMAO2dAS37+TZ1+IO6cXyh7ppZZndJQNQYhqG///h0Zbtd2n6sSf/2p/12lwQQRpDcOgNBLV2zRVurm5Ttdmn5J6czPYOENyI3Td/smq75/p/267/ePmpzRUh2hBEkrWDQ1Fd/sV3/u7debpdDz352lsryM+wuCxgSd84s0xdvHi9J+uaLu/TithM2V4RkRhhBUjJNU0/8do/+Z1uNnA5DK//PNbq2It/usoAhteyWCbrn+nKZpvTln2/Xq/vq7S4JSYowgqT0b3/ar5907WL6vTtn6EOTi+0tCLCBYRh6/GNTdfuMUvmDph74783adKTR7rKQhAgjSCo+f1B/99s9WvHHUNPe4x+bqo9fPdLmqgD7OByGvnfnDH1wwnC1dwb12efe0+931tpdFpIMYQRJo6apTZ95+s969s3DkqSv3jpR984ZY29RQAxIdTn047+q1PVj89XS4dcDa7boid/sYZVWDBnCCJLCa/vq9dHvb9CW6iZlp7n09D2V+sK88XaXBcSM9FSn/mvxdbr/A6Edfle/dVgLn/6zaprabK4MycAwTTPmV7zxer3Kzc2Vx+NRTk6O3eUgjnjaOvXD/92vZ948LNOUpo3M0cq7KzW6gKtmgN6s33NSX/75Nnnb/RqWkaK//YspuuOqkXI6uOwdkenr+zdhBAmpvTOg/9h4RCtfOyhPW6ck6a+uH63/99EpSktx2lwdEPuONZ7V0jVbtPOER5I0aUS2vnbbJN00cThr8aDPCCNISu2dAb247YRW/HG/aj3tkqQJxVn62m2TuGIGiFCHP6Bn3zysVa8dVHO7X5J0bUW+vjJ/omaNGUYowWURRpA0TNPUjuMe/XLzcb20vSY8ElKam6aHb5mgT14ziuFlYACazvq06rWDem7jkXBT67jhmfp0ZZk+ec1I9nJCrwgjSGj+QFA7Tnj05v4G/XZHjd4/2RL+XElumhbfWKG/ur6cKRlgENV62vT9P+3XC1tPqL0zFEochjT3iuGaP7VYN44v1Oj8DEZMEBbVMLJy5Ur98z//s2prazV16lStWLFCc+fO7fX8119/XcuWLdPu3btVWlqqRx55REuWLOnz4xFG0Nrh1946r3Ye9+itg6f19sHTau7whz/vdjl027QR+nTlKM0ZV8hICBBFze2dWrezVr/cfFzvHTnT7XOjhqXrxvGFum5svqaW5mpsYaZcTi7cTFZRCyNr167VPffco5UrV+qGG27QU089pWeeeUZ79uzR6NGjLzr/8OHDmjZtmv76r/9a999/v9566y0tXbpUzz//vD71qU8N6jeD+BYMmqpv7lB149nw7WB9i/bUenXkdKsu/EnNTU/RnHEF+uCE4VowvUS56Sn2FA4kscMNrfrN9hq9eaBBW6vPqDPQ/T9qqsuhicXZmlySrbHDs1Q2LEOj80O33Az+zya6qIWR6667Ttdcc41WrVoVPjZ58mR9/OMf1/Llyy86/2tf+5peeuklVVVVhY8tWbJE27dv15///Oc+PSZhJH4EgqZ8/qDaOgNq7fDrrC+gVp9fZzsC8rZ3ytMWunnbOnXmrE+nmjvCt4YWn3yB3hdZKs5xa0pJjmZV5OvG8YWaWprLCAgQQ876/Hr3cKPe3N+g7cebtKfGq1ZfoNfzs9wuDc92qzArtetPt/IyUpWT5lJueopy01OUnZaiTLdTGamu0J8pLrlTHHK7HEwHxYG+vn+7IvmiPp9Pmzdv1te//vVux+fPn6+NGzf2eJ8///nPmj9/frdjt956q5599ll1dnYqJeXiZNzR0aGOjo5u30w0/Grzce2q8XQ7NpQdND3lwAuPnH+K2fVZ65h50Tnmuc+ZofNDf17wsRn6SkHr76YUNE0FTVOBoLr+NMN/+oPn/xmUP2DKFwiqMxBUp99UZyCoDn9QHf7ARb8VRcrpMDQyL12j8zNUlp+hisIMTSnJ1eSSbBVkuQf0tQFEV0aqSzdNLNJNE4skhUY7j505qz01XlXVNav6dKuqG8/q2Jk2nWruUEuHXy0dfh1uaO3X47ldoVCS6nIq1WkoxeVQitMhl8NQitMhp8OQy2HI5TTkdBhyGKE/nYYhh8OQw5AcRui40fV3w5AMhf4uQzJ07phx/sfhHGSE/35+NDp3zOj28fkuPDSU4aqnh/rUNaM0bWTukNVwvojCSENDgwKBgIqLu18iWVxcrLq6uh7vU1dX1+P5fr9fDQ0NKikpueg+y5cv1+OPPx5Jaf3y+vun9NL2mqg/TrJKS3EoM9WlDLdTmaku5aSlKCfdpZz0FOWkpWhYRui3IetWmJWqETlpzC8DCcLhMFRekKnygkwtmN79tf6sz686T3t4VPRUc7tOtXR0jZz6Q3+2h0ZR23wBtfoCOuvzd/uFJ/RLUFCSXxi4q0cPi48wYrkwvZmmeclE19P5PR23PProo1q2bFn4Y6/Xq7Kysv6Uekm3TClWWX56v+5rXJRp+3i/C+7W41e54KSezuktdYfTu3HuuHFecndayd8wutJ/6AXDMM79lmD91uB0hD52Oc79huF0hv5MdTqU4nIo1emQy2kozeXsGjp1yu1yKC3FyRQKgF5lpLo0dniWxg7Piuh+vq5RWCuItHcG1NEZlD8YGq31dY3WdhvN7RrdtW7WKLApMzxCHAyaCpjn3p9Co8XnRpWl80alze6jztK5c84/L/xxT9/IBSf1dE5/RurNnh+tT64oiuzfYjBFFEYKCwvldDovGgWpr6+/aPTDMmLEiB7Pd7lcKigo6PE+brdbbnf0h+Rvn1Gq22eURv1xAACDI9XlUKrLoWy7C8Ggimg8PDU1VZWVlVq/fn234+vXr9ecOXN6vM/s2bMvOv+VV17RzJkze+wXAQAAySXiyflly5bpmWee0erVq1VVVaWHH35Y1dXV4XVDHn30US1atCh8/pIlS3T06FEtW7ZMVVVVWr16tZ599ll95StfGbzvAgAAxK2Ie0YWLlyo06dP64knnlBtba2mTZumdevWqby8XJJUW1ur6urq8PkVFRVat26dHn74Yf3oRz9SaWmpvv/97/d5jREAAJDYWA4eAABERV/fv7mGEgAA2IowAgAAbEUYAQAAtiKMAAAAWxFGAACArQgjAADAVoQRAABgK8IIAACwFWEEAADYKuLl4O1gLRLr9XptrgQAAPSV9b59ucXe4yKMNDc3S5LKyspsrgQAAESqublZubm5vX4+LvamCQaDqqmpUXZ2tgzDGLSv6/V6VVZWpmPHjrHnTR/wfPUdz1Xf8Vz1Hc9V3/Fc9V00nyvTNNXc3KzS0lI5HL13hsTFyIjD4dCoUaOi9vVzcnL4YY0Az1ff8Vz1Hc9V3/Fc9R3PVd9F67m61IiIhQZWAABgK8IIAACwVVKHEbfbrW9961tyu912lxIXeL76jueq73iu+o7nqu94rvouFp6ruGhgBQAAiSupR0YAAID9CCMAAMBWhBEAAGArwggAALAVYeQ8H/vYxzR69GilpaWppKRE99xzj2pqauwuK+YcOXJEixcvVkVFhdLT0zVu3Dh961vfks/ns7u0mPSd73xHc+bMUUZGhvLy8uwuJ6asXLlSFRUVSktLU2VlpTZs2GB3STHpjTfe0O23367S0lIZhqH/+Z//sbukmLV8+XLNmjVL2dnZKioq0sc//nHt27fP7rJi0qpVq3TllVeGFzubPXu2fv/739tSC2HkPPPmzdPPf/5z7du3T7/61a908OBBffrTn7a7rJizd+9eBYNBPfXUU9q9e7f+9V//VT/+8Y/1jW98w+7SYpLP59Odd96pBx54wO5SYsratWv10EMP6bHHHtPWrVs1d+5cLViwQNXV1XaXFnNaW1s1Y8YM/fCHP7S7lJj3+uuv6wtf+ILefvttrV+/Xn6/X/Pnz1dra6vdpcWcUaNG6R//8R+1adMmbdq0STfffLPuuOMO7d69e+iLMdGrF1980TQMw/T5fHaXEvO++93vmhUVFXaXEdOee+45Mzc31+4yYsa1115rLlmypNuxSZMmmV//+tdtqig+SDJfeOEFu8uIG/X19aYk8/XXX7e7lLgwbNgw85lnnhnyx2VkpBeNjY1as2aN5syZo5SUFLvLiXkej0f5+fl2l4E44fP5tHnzZs2fP7/b8fnz52vjxo02VYVE5PF4JInXp8sIBAL62c9+ptbWVs2ePXvIH58wcoGvfe1ryszMVEFBgaqrq/Xiiy/aXVLMO3jwoH7wgx9oyZIldpeCONHQ0KBAIKDi4uJux4uLi1VXV2dTVUg0pmlq2bJluvHGGzVt2jS7y4lJO3fuVFZWltxut5YsWaIXXnhBU6ZMGfI6Ej6MfPvb35ZhGJe8bdq0KXz+V7/6VW3dulWvvPKKnE6nFi1aJDNJFqmN9LmSpJqaGt1222268847dd9999lU+dDrz3OFixmG0e1j0zQvOgb014MPPqgdO3bo+eeft7uUmDVx4kRt27ZNb7/9th544AHde++92rNnz5DX4RryRxxiDz74oD7zmc9c8pwxY8aE/15YWKjCwkJNmDBBkydPVllZmd5++21bhq2GWqTPVU1NjebNm6fZs2fr6aefjnJ1sSXS5wrdFRYWyul0XjQKUl9ff9FoCdAfX/ziF/XSSy/pjTfe0KhRo+wuJ2alpqZq/PjxkqSZM2fqvffe07/927/pqaeeGtI6Ej6MWOGiP6wRkY6OjsEsKWZF8lydOHFC8+bNU2VlpZ577jk5HAk/yNbNQH6uEHoBrKys1Pr16/WJT3wifHz9+vW64447bKwM8c40TX3xi1/UCy+8oNdee00VFRV2lxRXTNO05T0v4cNIX7377rt69913deONN2rYsGE6dOiQvvnNb2rcuHFJMSoSiZqaGt10000aPXq0/uVf/kWnTp0Kf27EiBE2Vhabqqur1djYqOrqagUCAW3btk2SNH78eGVlZdlbnI2WLVume+65RzNnzgyPrlVXV9N71IOWlhYdOHAg/PHhw4e1bds25efna/To0TZWFnu+8IUv6Kc//alefPFFZWdnh0ffcnNzlZ6ebnN1seUb3/iGFixYoLKyMjU3N+tnP/uZXnvtNb388stDX8yQX78To3bs2GHOmzfPzM/PN91utzlmzBhzyZIl5vHjx+0uLeY899xzpqQeb7jYvffe2+Nz9eqrr9pdmu1+9KMfmeXl5WZqaqp5zTXXcPllL1599dUef4buvfdeu0uLOb29Nj333HN2lxZzPv/5z4f//w0fPtz80Ic+ZL7yyiu21GKYZpJ0ZwIAgJiUXBP9AAAg5hBGAACArQgjAADAVoQRAABgK8IIAACwFWEEAADYijACAABsRRgBAAC2IowAAABbEUYAAICtCCMAAMBWhBEAAGCr/w9nGSNeB/4KEAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.linspace(-3,3,101),psi_exact(np.linspace(-3,3,101)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerov method\n",
    "The second solution is a discrete position basis matrix approximation using Numerov's approximation for the kinetic energy.  (See paper cited above.)  After calculating the approximate Hamiltonian as a matrix, the approximate eigenstates and energies are found by diagonalizing the Hamiltonian.\n",
    "\n",
    "We'll also reuse the Hamiltonian from the Numerov method to calculate our variational energies in the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.453387Z",
     "start_time": "2019-01-18T21:41:04.437398Z"
    }
   },
   "outputs": [],
   "source": [
    "# bounds (These are overwritten if autoscale=True)\n",
    "xmin = -5.*a # lower bound of position\n",
    "xmax = 5*a # upper bound of position\n",
    "n = 200 # number of steps (may be overwritten if autoscale == True)\n",
    "dx = (xmax-xmin)/(n-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.478372Z",
     "start_time": "2019-01-18T21:41:04.460383Z"
    }
   },
   "outputs": [],
   "source": [
    "# the function V is the potential energy function\n",
    "def V(x):\n",
    "    # make sure there is no division by zero\n",
    "    # this also needs to be a \"vectorizable\" function\n",
    "    # uncomment one of the examples below, or write your own.\n",
    "    #return 0.5*mass*omega**2*x*x # harmonic oscillator\n",
    "    return np.piecewise(x, [np.abs(x)<= a, np.abs(x)>a],[0.0,V0])\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.504357Z",
     "start_time": "2019-01-18T21:41:04.481371Z"
    }
   },
   "outputs": [],
   "source": [
    "if (autoscale): \n",
    "    # set autoscale=True to use the scalings recommended by Pillai et al.\n",
    "    #Emax is the maximum energy for which to check for eigenvalues\n",
    "    Emax = 20.0\n",
    "    #The next lines make some reasonable choices for the position grid size and spacing\n",
    "    xt = opt.brentq(lambda x: V(x)-Emax ,0,5*Emax) #classical turning point\n",
    "    dx = 1.0/np.sqrt(2*Emax) #step size\n",
    "    # bounds and number of steps\n",
    "    n = np.int(0.5+2*(xt/dx + 4.0*np.pi)) #number of steps\n",
    "    xmin = -dx*(n+1)/2 \n",
    "    xmax = dx*(n+1)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.534338Z",
     "start_time": "2019-01-18T21:41:04.510353Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-5.0, 5.0, 200)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xmin, xmax, n #show the limits and number of steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.553327Z",
     "start_time": "2019-01-18T21:41:04.539335Z"
    }
   },
   "outputs": [],
   "source": [
    "#define the x coordinates\n",
    "x = np.linspace(xmin,xmax,n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.594301Z",
     "start_time": "2019-01-18T21:41:04.557326Z"
    }
   },
   "outputs": [],
   "source": [
    "#define the numerov matrices\n",
    "B = np.matrix((np.eye(n,k=-1)+10.0*np.eye(n,k=0)+np.eye(n,k=1))/12.0)\n",
    "A = np.matrix((np.eye(n,k=-1)-2.0*np.eye(n,k=0)+np.eye(n,k=1))/(dx**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.626280Z",
     "start_time": "2019-01-18T21:41:04.597299Z"
    }
   },
   "outputs": [],
   "source": [
    "#calculate kinetic energy operator using Numerov's approximation\n",
    "KE = -0.5*hbar**2/mass*B.I*A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.640278Z",
     "start_time": "2019-01-18T21:41:04.632279Z"
    }
   },
   "outputs": [],
   "source": [
    "#calculate hamiltonian operator approximation\n",
    "H = KE + np.diag(V(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.700235Z",
     "start_time": "2019-01-18T21:41:04.647268Z"
    }
   },
   "outputs": [],
   "source": [
    "#Calculate eigenvalues and eigenvectors of H\n",
    "energies, wavefunctions = np.linalg.eigh(H) # \"wavefunctions\" is a matrix with one eigenvector in each column."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, the eigenvalue calculation took 3.3s for _n_ = 500."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:04.715226Z",
     "start_time": "2019-01-18T21:41:04.704233Z"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8945998246109116\n"
     ]
    }
   ],
   "source": [
    "true_gs_energy = energies[0]\n",
    "print(true_gs_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational approximation\n",
    "Here we reuse the Hamiltonian matrix from above, but we use it to calculate the approximate energy expectation values for our trial wavefunctions.  Then we adjust the trial wavefunctions to minimize the energy expectation value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rectangle function for making it easier to write our trial wavefunctions\n",
    "def rectangle(x,a):\n",
    "    '''rectangle is one for -a<x<a and zero otherwise'''\n",
    "    return np.heaviside(x+a,0)-np.heaviside(x-a,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#trial wavefunction\n",
    "def psi(x,b=1.0):\n",
    "    '''Wavefunction as a function of position x.\n",
    "    The variable b is the parameter that we tune to minimize the energy.'''\n",
    "    # Note, normalization is calculated below.\n",
    "    #psi_unnorm = np.exp(-x**2*b) # gaussian\n",
    "    #psi_unnorm = np.cos(b*np.pi*x/2.0)*rectangle(x,1.0/b)# cosine piece\n",
    "    psi_unnorm = np.cos(b*np.pi*x/2.0)**2*rectangle(x,1.0/b)# cosine-squared piece\n",
    "    norm = 1.0/np.sqrt(np.dot(psi_unnorm,psi_unnorm)*(x[1]-x[0]))\n",
    "    return norm*psi_unnorm\n",
    "\n",
    "#limits on b (to help the minimization function)\n",
    "bmin = 0.1\n",
    "bmax = 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ev(psi): #ev of some wavefunction psi\n",
    "    '''Calculate expectation value of energy'''\n",
    "    return np.dot(np.dot(psi,H),psi)[0,0]*(x[1]-x[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "lines_to_next_cell": 2,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# minimize expectation value with respect to b\n",
    "result=opt.minimize_scalar(lambda b: ev(psi(x,b)),(bmin,bmax))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, the variational calculation took 0.3s with _n_=500."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of b is 0.6881, with an energy of 0.9234, compared to the energy from the Numerov diagonalization of 0.8946.\n"
     ]
    }
   ],
   "source": [
    "print(f\"The value of b is {result.x:.4f}, with an energy of {result.fun:.4f}, compared to the energy from the Numerov diagonalization of {true_gs_energy:.4f}.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-01-18T21:41:42.346083Z",
     "start_time": "2019-01-18T21:41:41.669500Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG2CAYAAACZEEfAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhJUlEQVR4nO3dd3xUVf7/8dfMpBKSQIAUegDpvWmQKkpTiqKiouCK+ENFV/laFlzr6qIruqwN1KXYK2JZFQEpgoAQJHQQaaEkhARISEidmd8fk0QiLQmT3Ds37+djZ2Xu3HvnM8Mh+ZzPOfdcm9vtdiMiIiJiUXajAxARERGpSEp2RERExNKU7IiIiIilKdkRERERS1OyIyIiIpamZEdEREQsTcmOiIiIWJqSHREREbE0JTsiIiJiaUp2RERExNIMTXZmzJhB+/btCQsLIywsjLi4OL7//vvzHrN8+XK6dOlCUFAQTZo0YebMmZUUrYiIiPgiQ5Od+vXr8/zzzxMfH098fDxXXHEFw4cPZ+vWrWfdf+/evQwZMoRevXqxYcMGpkyZwv3338+8efMqOXIRERHxFTaz3Qg0IiKCF198kXHjxp3x2qOPPsrXX3/N9u3bi7dNmDCBjRs3snr16soMU0RERHyEn9EBFHE6nXz22WdkZWURFxd31n1Wr17NgAEDSmwbOHAgs2bNIj8/H39//zOOyc3NJTc3t/i5y+Xi2LFj1KpVC5vN5t0PISIiIhXC7XZz8uRJ6tati91etoEpw5OdzZs3ExcXR05ODtWrV2f+/Pm0bt36rPsmJycTFRVVYltUVBQFBQWkpqYSExNzxjFTp07l6aefrpDYRUREpHIdOHCA+vXrl+kYw5OdFi1akJCQwIkTJ5g3bx5jx45l+fLl50x4/lyNKRqFO1eVZvLkyUyaNKn4eXp6Og0bNuTAgQOEhYV56VOIXJwvvviCv/zlL/Ts2ZNvv/3W6HCkiqtXrx6ZmZls2LCBJk2aGB2OCAAZGRk0aNCA0NDQMh9reLITEBBAs2bNAOjatSvr1q3jP//5D2+++eYZ+0ZHR5OcnFxiW0pKCn5+ftSqVeus5w8MDCQwMPCM7UVXgImYQbVq1QBwOBxql2IaoaGhao9iOuWZgmK6dXbcbneJOTani4uLY9GiRSW2LVy4kK5du551vo6Ir9D8MTEjtUuxCkOTnSlTprBixQr27dvH5s2beeyxx1i2bBmjR48GPENQY8aMKd5/woQJ7N+/n0mTJrF9+3Zmz57NrFmzeOihh4z6CCJeZbKLI6WKUjsUqzF0GOvIkSPcdtttJCUlER4eTvv27VmwYAFXXXUVAElJSSQmJhbvHxsby3fffceDDz7I66+/Tt26dXnllVcYOXKkUR9BxCvUgxYzUrsUqzA02Zk1a9Z5X587d+4Z2/r06cOvv/5aQRGJGEs9ajEDs7ZDp9NJfn6+0WFIBQoICCjzZeWlYfgEZRFRD1rMySzt0u12k5yczIkTJ4wORSqY3W4nNjaWgIAAr55XyY6IiZi1Ry1Vi9naYVGiExkZSbVq1UyThIl3uVwuDh8+TFJSEg0bNvTq37OSHRET0A9vMSMztEun01mc6JxriRGxjjp16nD48GEKCgq8epW16S49F6nKzNajlqrJTO2waI5O0VpUYm1Fw1dOp9Or51WyI2ICZuhBi/yZmdqlmWKRilNRf89KdkRMxEw9aqm61A7FapTsiJiAeq1iRmqXlePxxx/nrrvuKvX+r732GsOGDavAiKxHyY6IiahHLWagdnjxhg4dypVXXnnW11avXo3NZuPXX3/lyJEj/Oc//2HKlCmlPvf48eNZt24dK1eu9Fa4lqdkR8QE1IMWM1K7LL9x48axZMkS9u/ff8Zrs2fPpmPHjnTu3JlZs2YRFxdH48aNS33uwMBAbrnlFl599VUvRmxtSnZETEQ9ajEDtcOLd8011xAZGXnGnQBOnTrFJ598wrhx4wD4+OOPSwxJHT16lOjoaP75z38Wb/vll18ICAhg4cKFxduGDRvGl19+SXZ2dsV+EIvQOjsiJqAetJiRWdul2+3m1KlThrx3aRc19PPzY8yYMcydO5cnnnii+JjPPvuMvLw8Ro8ezfHjx9myZQtdu3YtPq5OnTrMnj2bESNGMGDAAFq2bMmtt97KPffcw4ABA4r369q1K/n5+axdu5Y+ffp4/4NajJIdERNRj7occjLgtx/g8AZwu/7Y7vCHxr2gSV/w8+7S81Zn9nZ46tQpqlevbsh7Z2ZmEhISUqp977jjDl588UWWLVtGv379AM8Q1nXXXUfNmjVJSEjA7XZTt27dEscNGTKE8ePHM3r0aLp160ZQUBDPP/98iX1CQkKoUaMG+/btU7JTCkp2REzArD1o08pJh+3fwLavYc9ScOadfb9Vr0BgGDQfBK2HwSUDlfiUgdrlxWnZsiU9evRg9uzZ9OvXj927d7NixYri4aiiIaigoKAzjp02bRpt27bl008/JT4+/qz7BAcHG1bh8jVKdkRMxOw9asO53bBlHnz/KJxK/WN7rUug2ZXgH/zHtlNpnopPZjJs/tTzqN0Chr0KDS+t/Nh9iNnbYbVq1cjMzDTsvcti3LhxTJw4kddff505c+bQqFEj+vfvD0Dt2rUBOH78OHXq1Clx3J49ezh8+DAul4v9+/fTvn37M8597NixM46Ts1OyI2IC6kGXQvoh+HYS/LbA8zyiKXS4CVoNg8iWZz/G5YKDaz1VoI0fQ+pOmD0Quo+H/k9AYGjlxe+DzNoubTZbqYeSjHbjjTfy17/+lQ8//JB33nmH8ePHF3+vTZs2JSwsjG3bttG8efPiY4rm9IwaNYqWLVsybtw4Nm/eTFRUVPE+u3fvJicnh06dOlX6Z/JFuhpLxETM3qM2zPp34PVLPYmO3R/6ToF71kCfR86d6ADY7dDwMhj4HExcBx1HA25Y+xa8EQd7llfaR/AlaofeU716dUaNGsWUKVM4fPgwt99+e/FrdrudK6+88oz1ch577DHS09N55ZVXeOSRR2jVqlXx1VtFVqxYQZMmTWjatGllfAyfp2RHxATM2oM2nNsNy56Hb+6HvJNQvztMWAl9Hy373JtqETDiDbhtPtRoCOkH4P2RsO2riondAtQuvWPcuHEcP36cK6+8koYNG5Z47a677uLjjz/G5fJMrl+2bBnTp0/nvffeIywsDLvdznvvvcfKlSuZMWNG8XEfffQR48ePr9TP4cs0jCViIupRn8bthkVPeCYZA/SdDL0fBrvj4s7b9ApPVeiribD1C/jsLzBiBnQYdfExW4TaoXfFxcWd8zsdMGAA9erV45NPPuHmm2+mb9++xXd6L9KwYUNOnDhR/HzLli0kJCTw6aefVmTYlqLKjogJqAf9Jy4XfPfQH4nOoOeh798uPtEpEhACI/8LHW8FtxPm/z+In+Odc1uI2mXFs9lsvPXWWxQUFJT6mMOHD/Puu+8SHh5egZFZiyo7IiaiHjWeROfriZDwAWCDof+BLmO9/z52h+fKrIBqnjk8/3sACnLgsru9/14+Ru2wcnXo0IEOHTqUev/TFxeU0lFlR8QE1IM+zfIXPImOzQHXvV0xiU4Rux0G/wsuf8DzfMHfYMe3Ffd+PkbtUqxClR0RE6nyPeod38HywpVih70C7W8o1WHJ6Tn8npKJmz++Pz+7ndYxYYRX8z//wTYbXPkU5J/yVHi++H8wfgnUaX7+40TEZyjZETEB9aCB1F2euTMA3e+CTreee9fMXJbuSGHt3mOs3XeM/WlnX0XWZoOW0WFcGhtB99gI+raoQ7WAs/zYs9lg4D8heQskroJPRsOdP0JQmDc+mc9SuxSrULIjYiJVtrKTexI+Hg25GdCwhyfxOItDJ7J5a/luPl53gNyCP+6DZbdBbO0QAvz+mMCclVtA4rFTbE/KYHtSBnNX7aNmNX/G9YzltrjGhAf/qeLj8Icb34E3+0Dqb/Dl3XDje56hriqkyrZBsTQlOyImUKV70C4XzJ/gWd04NAZumOtJPE6zPy2LN5buZt6vBylweX4Zt60XRu9L6tA9NoIujWoSGnTmcFXKyRzi9x1n7d5j/LjjCAeOZTNt4W+8uXwPY3s0ZlzPWGqGnLZeT/VIGPU+zBkEO/4HK1/yXO5eRVXpdimWomRHxESqZK96/WxPYuEI8FRSQv9YEt/tdvP+mv08++324kpOj6a1mHhFM+Ka1LrgL+PI0CCGtIthSLsY/n51K77dnMTrS3/ntyOZvLb0dz5ed4CXb+xA7+an3V+ofhe4+iX4+j5Y8hzE9oEG3Svko5tRlWyDYnlVqz4rYlJVtgd94gAsetLz56v+AQ26Fb90LCuPu95bz+NfbSW3wEWPprWYd3cPPhx/GT2a1i7zd+bnsDO8Yz0W/LU3M2/tQrPI6qRm5jJm9lqe+3YbeacNi9F5DHS4GXB7Fh8syPXCh/U9VbZdiuUo2RExkSrVq3a74X8PQl4mNLjMMym50KrdqQz+z08s2naEAIedx69pzfvjLqVLo5oX/bZ2u41BbaP53309GRPXCIC3V+zluhk/s/voaXfSHvhPCIn0DK/99OJFv6+vqFJtUKoMJTsiJlAle9CbPoXfF3mGr4a9WjwReN76g9z63184kpFLkzohfHFPD8b1jMVu9+53FOTv4JnhbXl7TFdqVvNny6EMRrz+M+v3H/fsUC0Crp7m+fPKf0PyZq++vy+oku3Si26//XZsNtsZj0GDBlXK+z/11FN07NixUt7L7JTsiJhIlelVZx6FBY96/tzn0eI1bd5fs5//+2wjLjdc26ke/7uvJ23rVeyS+Fe1juL7v/amS6OanMwp4LZZv7Bqd6rnxdbDodVQcBV4hrOcpV/S31dVmTZYSQYNGkRSUlKJx0cffWR0WFWOkh0RE6hyPejvH4Hs4xDVDi7/KwD/XbGHv3+5BYDbezTmpRs6nH1NnAoQHR7Ee+O607NZbU7lOfnLnHUs25nieXHINAgKh6QEWPN6pcRjFlWuXVaAwMBAoqOjSzxq1qzJsmXLCAgIYMWKFcX7vvTSS9SuXZukpCQAFixYQM+ePalRowa1atXimmuuYffu3SXOf/DgQW666SYiIiIICQmha9eu/PLLL8ydO5enn36ajRs3FleU5s6dW5kf3VR0NZaIiVSJXvXOBZ67jdscMPw13HY/Xv1xFy8v+g2Ae/o25eGBLSr9F221AD/+O7YrEz/8lcXbUxj/bjyv3tyZQW2jPfN3vroXlv7TU+mJaFKpsVUmX2iDbreb7HynIe8d7O/wStvs27cvDzzwALfddhsbN25k3759PPbYY3z00UfExMQAkJWVxaRJk2jXrh1ZWVk88cQTXHvttSQkJGC328nMzKRPnz7Uq1ePr7/+mujoaH799VdcLhejRo1iy5YtLFiwgMWLFwNU6RuHKtkRMYEq04N25sPCxzx/jrsX6nZk9sq9xYnOQwOaM/GKSwwLL8jfwYxbu/DAJwl8uymJez/8lXfv6M7lHUd75hjtXQ6Ln/YsPlgFmLVdZuc7af3ED4a897ZnBpap4vi///2P6tWrl9j26KOP8vjjj/Pss8+yePFi7rrrLrZu3cptt93GtddeW7zfyJEjSxw3a9YsIiMj2bZtG23btuXDDz/k6NGjrFu3joiICACaNWtWvH/16tXx8/MjOjq6PB/VUjSMJWIivtCrvijr50La71CtNvR+mCU7jvDst9sAeHRQS0MTnSL+Djuv3NSJoR3q4nS5ufv99exOzSpc1dkG276EA2uNDrPCWL4NVrJ+/fqRkJBQ4nHvvfcCEBAQwPvvv8+8efPIzs5m+vTpJY7dvXs3t9xyC02aNCEsLIzY2FgAEhMTAUhISKBTp07FiY6cmyo7IiZg1h60V+VkwLLCm3z2/Rs7TsB9H27A7YabuzdgQh/zDA057DZevL49h09ks37/ccbNXcf8ey6nZqfRsOF9WPh3uOMHzz21LMys7TLY38G2ZwYa9t5lERISUqLa8merVq0C4NixYxw7doyQkJDi14YOHUqDBg14++23qVu3Li6Xi7Zt25KXl+eJJTi4HJ+galJlR8RELN2r/nk6nEqFWs042vxmxs2NJyvPSY+mtXhmeFvT/WIN8nfw5m1dqF8zmH1pp7j7g/Xk9ZoMfsFw4BfY/o3RIVYIX2iDNpuNagF+hjy82U53797Ngw8+yNtvv81ll13GmDFjcLk8i1umpaWxfft2/v73v9O/f39atWrF8ePHSxzfvn17EhISOHbs2FnPHxAQgNNpzNwms1GyI2ICZvtF73XpB2G150qmvH5PcdeHGzl0IpvY2iG8Mboz/g5z/iiqXT2QWWO7UT3QjzV7jvH4kmO4e0z0vLj4SSjIMzbACmb5dlkJcnNzSU5OLvFITU3F6XRy2223MWDAAP7yl78wZ84ctmzZwksvvQRAzZo1qVWrFm+99Ra///47S5YsYdKkSSXOffPNNxMdHc2IESP4+eef2bNnD/PmzWP16tUANG7cmL1795KQkEBqaiq5uVVzJXBQsiNiKr7Qqy6XJc9BQQ407MEzuxqxIfEE4cH+zBrblRrVAi58vIFaRIfy6s2dsNvgk/gDzAsaCSF14NgeWD/H6PC8zrJt0CALFiwgJiamxKNnz54899xz7Nu3j7feeguA6Oho/vvf//L3v/+9+Gqrjz/+mPXr19O2bVsefPBBXnyx5EreAQEBLFy4kMjISIYMGUK7du14/vnncTg8Q20jR45k0KBB9OvXjzp16lTp9X1s7irWsjMyMggPDyc9PZ2wsDCjwxEBYNmyZfTr14/WrVuzdetWo8PxrqRN8GZvwM3a/p9x47f5ALx7R/eSN+A0uZnLd/P89zsI9newov8+ai97FIIj4K8JnnV4LCInJ6d4LkhGRgahoaGGx7N3715iY2MJCgoyNBapeOf7+76Y39+q7IiYiCX7HkueBdzktLiWu5d5hkXu7BnrU4kOwF29mhDXpBbZ+U7Gb26Fu3YLyD4Gq14zOjSvsmQblCpPyY6ICVh2bsThBNj1A26bnSdPDictK4+W0aE8PKiF0ZGVmd1u46UbOxAW5MeGQ5l8U+svnhd+eRNy0o0NroJYtl1KlaNkR8RELNerXuG5kebe6EF8sieAQD87r9zciUC/sl2+axZ1awQz9br2ADy4qT7ZNS6B3HRY+7bBkXmP5dqgCEp2REzBkj3oI9tg+ze4sTHxwBUATB7ckuZRxs4BuVhXt49hZOf6ON12/pV1jWfj6tchN9PYwCqAJdulVElKdkRMxFK96hWeS2hXB17OtoK69Gleh7E9Ghsbk5c8Naw1DSKCefdkZ1ID6nvm7ljkyixLtUGRQkp2REzAcj3o1N89N/sEns0YQvVAP14Y2d4ynzM0yJ9/jeyAEwf/yhri2fjzK5CfbWxgXmaVvy8RJTsiJmKZXvXKl8HtYhld2OZuzKSrmhMdbq3LhuOa1uK6TvX4wtmTFHsdyEqBX98zOqyLZpk2KHIaJTsiJmCpHvTxfbDxYwD+nTuc1jFhjIlrZGxMFWTykFZUCwrildzCuTs/T4cC66xSa6l2KVWakh0RE7FEr3rldHA7+cnZjk0049lr2+Jn0ttBXKw6oYE8PKglnzn7cMRdEzIOwUbfXqXWEm1Q5E+s+RNIxMdYpgedlYo74UMAXisYwU3dGtK5YU2Dg6pYt3RvSMv6dXir4GrPhtWvQ+HNHH2dZdqlhezbtw+bzUZCQoIpzlNWc+fOpUaNGpX6nmBwsjN16lS6detGaGgokZGRjBgxgp07d573mGXLlmGz2c547Nixo5KiFqk4Pt+rXj8HmzOXja4m7A5uz6M+uHhgWTnsNp67th2fuvqR4Q6G1N9g9xKjwyo3n2+DFtegQQOSkpJo27ZtqY+5/fbbGTFixEWfx5cZmuwsX76ce++9lzVr1rBo0SIKCgoYMGAAWVlZFzx2586dJCUlFT8uueSSSohYpGJYogddkIfrF8/ierMLBvG3Ia1Mf5NPb2lbL5yRca341NkXAHfhHd59nSXapYXk5eXhcDiIjo7Gz8/vos7lrfP4CkOTnQULFnD77bfTpk0bOnTowJw5c0hMTGT9+vUXPDYyMpLo6OjiR9FdXkV8mU/3qrfOx551hCPuGiRGD2Bk5/pGR1SpHrjyEub5XY3TbcO2ZwmkbDc6pHLx6TZoIm+++Sb16tXD9achzWHDhjF27Fh2797N8OHDiYqKonr16nTr1o3FixeX2Ldx48Y8++yz3H777YSHhzN+/Pgzhp+cTifjxo0jNjaW4OBgWrRowX/+85/iczz11FO88847fPXVV8UjIcuWLTvrMNby5cvp3r07gYGBxMTE8Le//Y2CgoLi1/v27cv999/PI488QkREBNHR0Tz11FMlYn755Zdp164dISEhNGjQgHvuuYfMTOMX3DTVnJ30dM/9ZSIiIi64b6dOnYiJiaF///4sXbr0nPvl5uaSkZFR4iFiNj7fg3a7yVv5KgDvFgzgoSHtsNt9/DOVUY1qAQzv14OFrq4AFKyeYXBEF8+07dLthrwsYx6lTAZvuOEGUlNTS/x+On78OD/88AOjR48mMzOTIUOGsHjxYjZs2MDAgQMZOnQoiYmJJc7z4osv0rZtW9avX8/jjz9+xvu4XC7q16/Pp59+yrZt23jiiSeYMmUKn376KQAPPfQQN954I4MGDSoeCenRo8cZ5zl06BBDhgyhW7dubNy4kRkzZjBr1iyeffbZEvu98847hISE8Msvv/Cvf/2LZ555hkWLFhW/brfbeeWVV9iyZQvvvPMOS5Ys4ZFHHinVd1aRTFO/crvdTJo0iZ49e553DDEmJoa33nqLLl26kJuby3vvvUf//v1ZtmwZvXv3PmP/qVOn8vTTT1dk6CJe47O96sTVBBzdTI7bnz2NbuDhZrWNjsgQt/dozAMrhzM4f53n8vsrn4SQWkaHVSY+0QbzT8E/6xrz3lMOQ0DIBXeLiIhg0KBBfPjhh/Tv3x+Azz77jIiICPr374/D4aBDhw7F+z/77LPMnz+fr7/+mokTJxZvv+KKK3jooYeKn+/bt6/E+/j7+5f4HRcbG8uqVav49NNPufHGG6levTrBwcHk5uYSHR19znjfeOMNGjRowGuvvYbNZqNly5YcPnyYRx99lCeeeAK73VMbad++PU8++SQAl1xyCa+99ho//vgjV111FQAPPPBAiVj+8Y9/cPfdd/PGG29c8DurSKap7EycOJFNmzbx0Ufnv2yzRYsWjB8/ns6dOxMXF8cbb7zB1VdfzbRp0866/+TJk0lPTy9+HDhwoCLCF7kopu1Bl1L60lcAmO/qycRrLjU4GuME+Tu4YsBwNrli8XPlkrNmltEhXRRfb5dGGz16NPPmzSM317P20gcffMBNN92Ew+EgKyuLRx55hNatW1OjRg2qV6/Ojh07zqjsdO3a9YLvM3PmTLp27UqdOnWoXr06b7/99hnnuZDt27cTFxdX4u/88ssvJzMzk4MHDxZva9++fYnjYmJiSElJKX6+dOlSrrrqKurVq0doaChjxowhLS2tVHNxK5IpKjv33XcfX3/9NT/99BP165d9nP+yyy7j/fffP+trgYGBBAYGXmyIIpXCJ3rVf+I+tpfQfT8AkHjJWG6uG25wRMYa2aUB05ZcS/vsl8lf8xZBfR4AP9+ZqO0TbdC/mqfCYtR7l9LQoUNxuVx8++23dOvWjRUrVvDyyy8D8PDDD/PDDz8wbdo0mjVrRnBwMNdffz15eXklzhEScv4q0qeffsqDDz7ISy+9RFxcHKGhobz44ov88ssvZfpYbrf7jOS2qC2cvt3f37/EPjabrXhe0v79+xkyZAgTJkzgH//4BxEREaxcuZJx48aRn59fpni8zdBkx+12c9999zF//nyWLVtGbGxsuc6zYcMGYmJivBydSOXx5R70wYWv0AAXK13tuOWagUaHYziH3Ub3a+7gyKezico/yrH4z4i4bLTRYZWLadulzVaqoSSjBQcHc9111/HBBx/w+++/07x5c7p06QLAihUruP3227n22msByMzMPGOIqjRWrFhBjx49uOeee4q37d69u8Q+AQEBOJ3O856ndevWzJs3r0TSs2rVKkJDQ6lXr16pYomPj6egoICXXnqpeNiraO6Q0Qwdxrr33nt5//33+fDDDwkNDSU5OZnk5GSys/+4md7kyZMZM2ZM8fPp06fz5ZdfsmvXLrZu3crkyZOZN29eiTFOEV/lE73q07hyMonY4bk1xIEWf6FBROl7vVbWt3V9loYOByBz+asGR1M2vtYGzW706NF8++23zJ49m1tvvbV4e7Nmzfjiiy9ISEhg48aN3HLLLWdcuVUazZo1Iz4+nh9++IHffvuNxx9/nHXr1pXYp3HjxmzatImdO3eSmpp61irLPffcw4EDB7jvvvvYsWMHX331FU8++SSTJk0qTlwupGnTphQUFPDqq6+yZ88e3nvvPWbOnFnmz1QRDE12ZsyYQXp6On379iUmJqb48cknnxTvk5SUVGLsMS8vj4ceeoj27dvTq1cvVq5cybfffst1111nxEcQ8QrT9qAvYPviOYRwikR3FIOG+2b1oiLYbDZaD72PPLeDhtnbSdqx1uiQysVX26WZXHHFFURERLBz505uueWW4u3//ve/qVmzJj169GDo0KEMHDiQzp07l/n8EyZM4LrrrmPUqFFceumlpKWllajyAIwfP54WLVoUz+v5+eefzzhPvXr1+O6771i7di0dOnRgwoQJjBs3jr///e+ljqVjx468/PLLvPDCC7Rt25YPPviAqVOnlvkzVQSbu4ql8RkZGYSHh5Oenk5YWJjR4YgAsHbtWi699FIaNWpUrlK2EdxuN7ue7U5z528saziRvnc8Z3RIprPmhWFclr2cNRHDuez+d40Op1TS0tKoXdtzNV1BQYHha5jl5OSwd+9eYmNjCQoKMjQWqXjn+/u+mN/fprkaS6Qq88UedPyaZTR3/kae20H7a+658AFVUI3eEwBom/YDh1OOGhxN2fliuxQ5GyU7IibiK4VWt9vN8eVvArCr1hVERJZuAmNV0/KywRzyq091Ww7rvnnL6HBKxVfaoEhZKNkRMQFf60H/sn0fl2d7bnZZ98q7DY7GxGw28tp7LrBolvgpKenZFzjAXHytXYqci5IdERPxlV715gWzCLHlkhrYkJqtrjA6HFNr3P9O8vCnjW0fXy/41uhwLshX2qBIWSjZETEBX+pBr92TxuUnvgYg4LJxnjVP5JxsIbU41ngIADW2vkdqZq7BEZWemdqlkrCqoaL+npXsiJiIL/xA/27BN7S27yffFkDYpWMufIAQ1c8z1DfEtpr3lm42OJrzM1sbLFqx99SpUwZHIpWhaAVpb18FaIrbRYhUdWbqQZ/PlkPptD78BfhBXoth+FeLMDokn2BreBmZYZdQPWMXWfEfkHFVB8KC/C98oOBwOKhRo0bx/ZeqVavmM/9epGxcLhdHjx6lWrVq+Pl5Nz1RsiNiImbrVf/Ze8s28ZRjNQAhPe4yOBofYrMRcvl4+P4Rrncv4tO1idzZu6nRUZ2VGdtg0d26T7/hpFiT3W6nYcOGXk9oleyImIAv9FQPn8jGf/t8gv3yyKlxCUENuhsdkk+xtR9FwQ+P05IDvLjyR26/PBY/h3lnEpipTdpsNmJiYoiMjDT8hpJSsQICAkp9e4qyULIjYiJm7FUXeWf1PkbalwMQ1G2MJiaXVXANbK2uga3z6HNqId9vuYahHeoaHdUZzNwGHQ6H4Ss6i28yb7dCpAoxUy/6bLJyC1j9yyo62X/HZXNA+1FGh+STHJ09N4Ic7viZd37aYerEwuxtUqQslOyImIhZf/l9Fn+AIQWeRQRtl1wFoVEGR+SjYvvgDK1HuO0U0UlLWL//uNERncGsbVDkYijZETEBM/einS4376zczXWOlQDYOt1qcEQ+zO7A0clz5+sbHMv574q9Bgd0bmZukyJlpWRHxETM2KtetC2Z2PTVRNpO4A6uBZcMNDok39bRk+z0sm9m47at7E/LMjigkszYBkUulpIdERMwcy/6vyv2coPDMzHZ1mEU+AUYHJGPi2gCjS7HbnNzrX0Fc37eZ3REZ2XmNilSVkp2REzEbL3qjQdOsHv/fvrbf/Vs6Dja2ICsovB7vN7xE5/GJ5KRY57Lqc3WBkW8QcmOiAmYtRf93pr9DHesIsDmhJgOEN3W6JCsofVw3P4hNLEn0yp/O/N/PWR0RGcwa5sUKQ8lOyImYqZe9YlTeXyz8XDxEBYdNTHZawKrY2szAvBMVH5vzX7T/N2bJQ4Rb1KyI2ICZuxFf77+IE2ce2lj34/bEQDtrjc6JGspHMoa6ljDoZRU1uw5ZnBAJZmxTYqUl5IdERMxS6/a5XLz/pr9XFt0uXnzQaCbfnpXox5QszEhthyusv/K+7/sNzoiwDxtUMSblOyImIDZetE/704lMS2TEY5Vng0dbjI2ICuy2aDdjQCMcKzkhy3JpGTkGBzUH8zWJkUuhpIdERMxS6/6vdX76WHfSqTtOATXhGZXGR2SNRXedqOPYxM1XCf4eN0BgwMyTxsU8SYlOyImYKZe9OET2SzefqR4CIs212ltnYpSuxnU64IDF9c4VvPR2kQKnC6jowLM1SZFLpaSHRETMUOv+uO1iQS6cxjiF+/ZoJt+VqzC7/d6/1Ukpefw444UQ8MxQxsU8TYlOyImYJZedL7TxUfrDnCV/VeC3dlQoxE06G50WNbW5jqwOWjL7zSxHeb9NeaYqGyWNiniDUp2REzE6F71om1HOHoylxsDCycmtx/lmUgrFad6HWh2JQAjHD+zYlcq+1KNu1+W0W1QpCIo2RExAbP0oj9ed4DapNPDvdGzQUNYlaO956qsmwJXA24+jTd+orJZ2qSINyjZETERI3vVh09ks2LXUa5xrMaOE+p18UyglYrXYggEVCfSmUwX22/M+/WgYROVVdkRK1KyIyKAZ8VktxturbbGs0FVncoTUA1aDQPgpsBVHMnI5addRw0OSsQ6lOyImEDRkIFRvWqXyzN00sR2mGb5v4HN4Zk4K5WncCjrGsca/Cng03UHDQmjqA1qGEusRMmOiLBmTxoHj2dzQ0BhVafpFZ6Js1J5YntDaAzBzpP0tSewePsRUjNzjY5KxBKU7IiYgNGVnU/iDwBubggsGsK60ZA4qjT7H9W0W0PXU+By8+WGQ5Uehio7YkVKdkSquPRT+Xy/JZm2tr3UzjsIfsGeCbNS+dqOBKBHwVqCyOWTdQc0YVjEC5TsiJiAkZWdrzceIq/Axe1hhSsmtxgEgdUrPQ4B6nWGGo3wd2Yz0H8ju1IySThwolJDUGVHrEjJjkgV90n8AWy4GMRqz4a21xsbUFVms0Fbz1DWHTV+BTDFmjsivk7JjogJGFXZ2Xo4nS2HMojz+43quUcgMBwu0R3ODVU4lNUuaw3VOcU3G5M4lVdQaW+vyo5YkZIdkSrs8/Wey5sn1Nrg2dBqKPgFGhiRENUWajfH7srjprDNZOYW8P3mZKOjEvFpSnZETMCIyk6+08XXCYfxo4DLsld4NrYbWWnvL+dgsxVXd0aHeOZRza/Eq7JU2RErUrIjUkWt3JVKWlYeg6vtICDvBITUgca9jQ5LoPgS9Mbpv1CDk/y8O5UjGTkGByXiu5TsiJiAEZWdLwqrBXcWToSlzbXg8Ku095fzqNMcotthcxUwoc5W3G74KqFyqjuq7IgVKdkRqYJO5uSzcGsygeTR9uRKz0ZdhWUuhUNZIwpXtf7i18pfYFDEKpTsiJhAZVd2vt+STG6Bi5tqbMeRnwnhDaFB90p5bymlNtcCEHUsnrqODHYkn2R7UkaFv60qO2JFSnZEqqD5hVWC26oXLiTY9lrPxFgxj5qNoV5XbG4X98VsBSp3orKIlSjZETGByqzsHD6RzZq9aVQjhyYnfvZsbKursEyp8O9lkNvz9/RVwiGcroptI6rsiBUp2RGpYr5KOIzbDeOjd2EvyIGIJhDd3uiw5GzajABs1Ez7lRZBJziSkcvq3WlGRyXic5TsiJhAZVV23G438zd4FhK8PnCdZ2MbDWGZVlhdaHQ5AA/EbAHgi8K/v4qiyo5YkZIdkSpk6+EMfjuSSU2/XOqnFV6FVTgRVkyq8F5ZvfM8Cz8u2JJcqbePELECJTsiJlBZlZ0vCye43lfvd2wFOVDrEs/tCcS8Wg8Hm4OQtM30qHGCU3lOFm49UmFvp8qOWJGhyc7UqVPp1q0boaGhREZGMmLECHbu3HnB45YvX06XLl0ICgqiSZMmzJw5sxKiFfFtTpebrzceBuAah2ftFg1h+YCQ2tCkLwD3Rm4CKP57FJHSMTTZWb58Offeey9r1qxh0aJFFBQUMGDAALKyss55zN69exkyZAi9evViw4YNTJkyhfvvv5958+ZVYuQi3lUZlZ11+46RcjKXmKA86hwpvBeWhrB8Q+FVWV0zlwLw029HOXEqr0LeSpUdsSJD14ZfsGBBiedz5swhMjKS9evX07v32e/RM3PmTBo2bMj06dMBaNWqFfHx8UybNo2RI3X5rMi5fFNYDfhrvV3YDuVB7RYQ2crgqKRUWl4N/wsg8NhOBkUeY0FKBAu2JHNT94ZGRybiE0w1Zyc9PR2AiIiIc+6zevVqBgwYUGLbwIEDiY+PJz8//4z9c3NzycjIKPEQMZuKruzkO118vyUZgKtY5dmoISzfEVwDml0JwJ01NwDwzaaKGcpSZUesyDTJjtvtZtKkSfTs2ZO2bc89YTI5OZmoqKgS26KioigoKCA1NfWM/adOnUp4eHjxo0GDBl6PXcTsVu1O41hWHo2q5RORpKuwfFLhUFaHE0sAN6t3p5FyUndCFykN0yQ7EydOZNOmTXz00UcX3PfPPY7z9UQmT55Menp68ePAgQPeCVjEiyq6slM0hHV//d+wufIhsjVEtqyQ95IK0nwQ+AXjn76XkdGpuNzw/eZkr7+NKjtiRaZIdu677z6+/vprli5dSv369c+7b3R0NMnJJf+Bp6Sk4OfnR61atc7YPzAwkLCwsBIPkaokt8DJD4VDWP2dhbeHUFXH9wRWhxaDABgbvh74I4kVkfMzNNlxu91MnDiRL774giVLlhAbG3vBY+Li4li0aFGJbQsXLqRr1674+/tXVKgiFaoiKzvLdx7lZG4BzUPzCS8awmo9wuvvI5WgjWeBwTbHlmCzuYnff5xDJ7K9+haq7IgVGZrs3Hvvvbz//vt8+OGHhIaGkpycTHJyMtnZf/zjnTx5MmPGjCl+PmHCBPbv38+kSZPYvn07s2fPZtasWTz00ENGfAQR0ytak+W+ejuxuQo8iwjWaW5wVFIul1wF/iE4Th7klrpHAfifqjsiF2RosjNjxgzS09Pp27cvMTExxY9PPvmkeJ+kpCQSExOLn8fGxvLdd9+xbNkyOnbsyD/+8Q9eeeUVXXYuPq2iKjun8gr4cXsKAH3yNYTl8/yDocVgAG6tXjiU5eWrslTZESsydJ2d0vxgnzt37hnb+vTpw6+//loBEYlYy+LtKWTnO2lbs4DQw7oKyxLaXgdbPqf5sR/xsw9iy6EM9hzNpEmd6kZHJmJappigLFLVVVRlp2gC630x27G5nRDdHmo19ep7SCVr2h8CQnGcPMyYBp6hrG82Jnnt9KrsiBUp2RGxqJM5+Sz/zfPL8PJcVXUswz8IWg4B4KZq8QB8v8V7yY6IFSnZETGBiqjsLNmRQl6Bi061Cgg5rPk6llJ4VVazo4vxt7vZkXyS3UczvXJqVXbEipTsiFjUd5s9vf27o7Zhc7ugbieIuPDyDuIDmvaDwHDsmcncXv8IAN9vVnVH5FyU7IiYgLd70Vm5BSzb6RnC6pHzk2ejqjrW4RfouTkocGPwOgC+8/JqyqrsiJUo2RGxoKU7U8gtcNGpZi4hSWs8G7WQoLUUJq9NU3/E3+5mW1IG+1KzLvq0FXXLEhEjKdkRMQFv96KL7pk0oWgIq14XqNnIq+8hBmvSF4JqYM9K4fa6nqvuiu5s7w2q7IiVKNkRMZmL7Vln5zlZssOzkODlOSs8GwsntIqF+AVAq2sAuMGLV2WpsiNWpGRHxAS82Yte/ptnIcEONbIJSf7Fs7H1cK+dX0ykcCirWeqP+NucbDqYzoFjp7xyalV2xEqU7IiYzMX2rIsmqt4buRUbbqjfHWo08EZoYjaxfSA4AvupVMYWDmUtuMihLFV2xIqU7IiYgLd60Tn5Tn7c7rkUuUdu0RCWrsKyLIc/tBoKwI3BawH41kuXoKuyI1aiZEfEZC6mZ/3Tb0fJynPSPiyL6kc8lyRrCMviioeyluJnc5Jw4ASHTmSX+3Sq7IgVKdkRMQFv9aKLrsa5N2qrZ0PDOAiv55Vzi0k17gXVamHPOcbYmAPAxQ9lgSo7Yi1KdkRMprw967wCF4u3eYawLtdCglWHww9aDQNgVOECgxezmrIqO2JFSnZETMAbvehVu1M5mVtAu+oZVD/6K2Ar/iUoFte28F5Zacvwo4D1icdJOZlzUadUZUesRMmOiMmUt2f9w1ZPVefuyC2eDY16QFiMt8ISM2t0OYTUwZ5znFuj9uF2w+JtKeU6lSo7YkVKdkRM4GJ70U6Xm0WFQ1g9dRVW1WN3FE9EHxXsWWDwh60XN29HlR2xEiU7IiZTnp71hsTjpGbm0iLoGGFpG8Fm1xBWVVO4Snbz48vxp4BVu1PJyMkv82lU2RErUrIjYgIX24su6sXfU6doCOtyCI262LDElzS8DKpH48hN54aau8h3ulm6o3xDWaDKjliLkh0Rkylrz9rtdhfP1+lTsNKzUUNYVc9pQ1m3hKwHYGFhuygLVXbEipTsiJjAxfSidySfJPHYKZr5pVDj+BYNYVVlhUluq/QVBJDPsp0p5OQ7y3UqVXbESpTsiJhMWXvWfwxhbfJsiO0D1et4OyzxBQ0uhdC6OPJPMrz6DrLynPz8e2qZTqHKjliRkh0RE7iYXnTRENYVBYVXYbUd6Y2QxBfZ7dBmBAC3hnqGssp7VZYqO2IlSnZETKYsPesDx06xPSmDlo6D1Di5C+z+0OqaCoxOTK9wKKvtyZ8JJI/F21MocLpKfbgqO2JFSnZETKC8veiiXvv4mgmeDc2uhOCaXopKfFK9rhBWH0dBFoODtnIsK4/4/cfLfBpVdsRKlOyImExZetaeq23cXOkqvApLQ1hy2lDWmLBfgbINZamyI1akZEfEBMrTi07NzGXd/mO0se0j/FQi+AVDi8EVEJ34nMIFBttnrSaIXBZuPVLmJEaVHbESJTsiJlPaX0pLtqfgdsMd4Z6JqDQfCIHVKzAy8Rn1OkONhvg5TzHAfxOHTmSzLSmjVIeqsiNWpGRHxATK04tetP0INlwMcP/s2aAhLClisxVPVL4tbANQ9huDqrIjVqJkR8RkStOzzsl3smLXUTrbdhGaewQCQuGSqyohOvEZhclOp5xfCCaHxdtLt5qyKjtiRUp2REygrL3olbtSycl3cXO1dZ4Nra4B/+AKiEx8VkxHqNkYP2c2/R0JbD6UTlJ6dqkPV2VHrETJjojJlKZnvXj7ERw4GWRb7dmgISz5s9OGskZX98zrWrz9wkNZquyIFSnZETGBsvSiXS43i7encKl9O9ULjnvW1WnSt+KCE99VeFVWt/x4Qshm8bbS3xhUlR2xEiU7IiZzoZ51wsETpGbmcp3/Gs+G1sPB4V8JkYnPiW4HtZrh58rlKvt6Vu9OIzO34LyHqLIjVqRkR8QEytKLXrztCP4UMNhROF9HQ1hyLjYbtL0egJuC15DndPHTb0dLeagqO2IdSnZETOZCPetF247Q076ZENdJqB4FjS6vpMjEJ7XzJDvdXBupScYFh7JU2RErUrIjYgKl7UXvS81iV0omwx2FE5PbXAt2RwVGJj6v9iUQ3R6H28kQx1qW7CzdjUFV2RErUbIjYjLn61kv3n6EQPIY6Fe4arKGsKQ02t0AwLX+qzlxKv+8NwZVZUesqFzJTlZWlrfjEKnSStuLXrTtCP3sCQS7syG8IdTvVsGRiSW09VyV1ZkdxJBWqquyVNkRKylXshMVFcUdd9zBypUrvR2PSJV3rp718aw84vcfZ6hjlWdD22s9E1BFLiS8PjTsgR031zhWs2j7uW8MqsqOWFG5kp2PPvqI9PR0+vfvT/PmzXn++ec5fPiwt2MTqTJK04te9lsKQa5TXOlI8GzQEJaURTtPexnut5r9aaf4PSXzvLursiNWUq5kZ+jQocybN4/Dhw9z991389FHH9GoUSOuueYavvjiCwoKzr+Og4ic27l61kt2HOVK+3oCyYNazSC6fSVHJj6t9QiwOWhr20usLYklO86+mrIqO2JFFzVBuVatWjz44INs3LiRl19+mcWLF3P99ddTt25dnnjiCU6dOuWtOEUs7UK96AKni+U7UxjqOO32EOp5S1mE1Iam/QAYZl/Fj+dIdoqosiNWclHJTnJyMv/6179o1aoVf/vb37j++uv58ccf+fe//838+fMZMWKEl8IUqTrO1rNev/84tpwT9HFs8mwovA2ASJkULjA4zLGK9fuPceJU3hm7qLIjVuRXnoO++OIL5syZww8//EDr1q259957ufXWW6lRo0bxPh07dqRTp07eilPE0i7Ui16yM4XBjrX444SothDZspIiE0tpeTX4BdG0IIlW7j0s/+0owzvWO+uuquyIlZSrsvOXv/yFunXr8vPPP5OQkMDEiRNLJDoATZo04bHHHvNGjCJVytl61ku2p3Cto/Dqx8I1U0TKLCgMWgwG4FrHz2edt6PKjlhRuSo7SUlJVKtW7bz7BAcH8+STT5YrKJGq5ny96APHTpF9dC+XBu7AjQ2bkh25GO1Hwdb5DHOs4o2dyRQ4Xfg5zuz3qrIjVlKuZKegoICMjIwztttsNgIDAwkICLjowESqqj/3rJfuTGG4/WcAbLG9IPzsww4ipdLsStzBEdTJPkbb3A1sOHAp3RpHFL+syo5YUbmGsWrUqEHNmjXPeNSoUYPg4GAaNWrEk08+ict14fuviMj5e9E/bjvyxxBW+5sqKSKxLIc/tsI1mkY4VvLj9rNflaXKjlhJuZKduXPnUrduXaZMmcKXX37J/PnzmTJlCvXq1WPGjBncddddvPLKKzz//PPnPc9PP/3E0KFDqVu3LjabjS+//PK8+y9btgybzXbGY8eOHeX5GCKmdHrP+lReASf3rqOZ/TAuRyC0GmpgZGIZHTxJ80B7PKu27yvxkio7YkXlGsZ65513eOmll7jxxhuLtw0bNox27drx5ptv8uOPP9KwYUOee+45pkyZcs7zZGVl0aFDB/7yl78wcmTpV4PduXMnYWFhxc/r1KlTno8hYnqrfk/jGlYAYGt1jWeCqcjFqtcFZ80mVDu+h6apyzhwrDcNIkrOw1RlR6ykXMnO6tWrmTlz5hnbO3XqxOrVnkXPevbsSWJi4nnPM3jwYAYPHlzm94+MjDzj6i8Rqzi9Z71sx2H+WngvLFv7UUaFJFZjs+HocBMs+yfXOlayZEcKY3s0BlTZEWsq1zBW/fr1mTVr1hnbZ82aRYMGDQBIS0ujZs2aFxfdOXTq1ImYmBj69+/P0qVLz7tvbm4uGRkZJR4iZvTnnrTb7SZr22Lq2DLIC4yAplcYFJlYUntPZf5y+xbWb9l2xsuq7IiVlKuyM23aNG644Qa+//57unXrhs1mY926dezYsYPPP/8cgHXr1jFqlHd7ojExMbz11lt06dKF3Nxc3nvvPfr378+yZcvo3bv3WY+ZOnUqTz/9tFfjEKlIRT3r7Ukn6Zu7BBxgb3c9OPwNjkwsJSKW7OiuBCfHU/fAt2TlDiAk0E+VHbGkciU7w4YN47fffmPmzJns3LkTt9vN4MGD+fLLL2ncuDEAd999tzfjBKBFixa0aNGi+HlcXBwHDhxg2rRp50x2Jk+ezKRJk4qfZ2RkFFefRMzEZrOV+EWzcutebrPHA+DXSVdhifcFdbkFvo1nqG0FK39PZWCb6OLXVNkRKylzspOfn8+AAQN48803mTp1akXEVCaXXXYZ77///jlfDwwMJDAwsBIjErk4RQlPzqYvCbblkRHSiLC6nQ2OSqzI1uZanN89Qhv7fhZuWM3ANteqsiOWVOY5O/7+/mzZssU0Wf+GDRuIiYkxOgyRi3b6v6ljWXl0T//e86T9TbrDuVSMahEcr++ZCxa5Zx4u1x+Jjll+xot4Q7kmKI8ZM+asE5TLKjMzk4SEBBISEgDYu3cvCQkJxVdxTZ48mTFjxhTvP336dL788kt27drF1q1bmTx5MvPmzWPixIkXHYuIWbjdbtZtWM9l9u24sBF22ZgLHyRSTuFxtwMwwLmcLQdTVdkRSyrXnJ28vDz++9//smjRIrp27UpISEiJ119++eVSnSc+Pp5+/foVPy+aWzN27Fjmzp1LUlJSicvX8/LyeOihhzh06BDBwcG0adOGb7/9liFDhpTnY4iYyuk9adevHwCwP7w7seH1jQpJqgD/FgNJd9SkjvM4q1Z9iT3EM29HlR2xknIlO1u2bKFzZ88cgt9++63Ea2X5B9K3b9/z9iLmzp1b4vkjjzzCI488UvpARXxQfn4+HY99B4C742iDoxHLc/iR1Gg44XvmErX7M1LaqVIu1lOuZOdCa9uISNkVdRSStyznUtLIIIRGl994gaNELl5k73GwZy5d89YxP+sEoMqOWEu55uwU+f333/nhhx/Izs4GtPKmiDcE7fwCgE01r8IREGxwNFIVRDRuz06/FvjZXFRL+tnocES8rlzJTlpaGv3796d58+YMGTKEpKQkAO68807+7//+z6sBilQVNpuNGkHQKsNzh3NnBw1hSeU51Pg6ANqdLLwXmyo7YiHlSnYefPBB/P39SUxMpFq1P24eN2rUKBYsWOC14ESqmtHdahFAPjtcDejYra/R4UgVEn35aHLc/jR0H6JbgwCjwxHxqnIlOwsXLuSFF16gfv2SV4lccskl7N+/3yuBiVRFf+kUBMAv4YMJD9EvHKk8rRrXZ6kjDoDxl0caHI2Id5Ur2cnKyipR0SmSmpqq1YpFyqlNHTtdap4k3+3ArTucSyWz2WwcbDwSgFFNswlyaA6mWEe5kp3evXvz7rvvFj+32Wy4XC5efPHFEuvmiEjpje3oqeT86OpMj/YtLrC3iPfFdhnIAVcdwvzy6ROZYXQ4Il5TrkvPX3zxRfr27Ut8fDx5eXk88sgjbN26lWPHjvHzz5rJL1Jm+Tnc1s4BOFkU0J9pkdWNjkiqoB6X1OEtV28esM9jWMOTRocj4jXlquy0bt2aTZs20b17d6666iqysrK47rrr2LBhA02bNvV2jCLWt+0rIgKdHHTXxt24l66EEUNUC/BjXbXeFLjtdApPh5TtRock4hXlquwAREdH8/TTT3szFpEqyx0/GxvwUcEVXNoo3OhwpAqLqlmDxYe7MMixDuLnwJB/GR2SyEUrd7Jz4sQJ1q5dS0pKCi6Xq8Rrp9+8U0Qu4MhWbAfWkO928EnOpXxQN+TCx4hUkFbhTj440J9BjnW4Nn6E/conIUBtUnxbuZKdb775htGjR5OVlUVoaGiJkrvNZlOyI1IW8XMAWOjqQuKevQQ4NIQlxqkZ4GZJcgj760fSKDcFtnwBnW8zOiyRi1KuOTv/93//xx133MHJkyc5ceIEx48fL34cO3bM2zGKWFduJmz8GIAPnFdy6vdfDA5IBE7tXs+Hzv6eJ/GzjQ1GxAvKlewcOnSI+++//6xr7YhIGWz5HPJOsscVzSpnS7L3xOsec2Iot9tN9u61fObsQz4OOPwrHN5gdFgiF6Vcyc7AgQOJj4/3diwiVYvbDetmAfChsz/5ybtxnUo3OCgRyD38Gyfy/fjOealnQ+FQq4ivKtecnauvvpqHH36Ybdu20a5dO/z9/Uu8PmzYMK8EJ2Jph3+F5E3k4888Zy/y9n0LoMqOGMrtdoPbRfCJvXwQ0J/hjlWw+XMY8A8I0pWC4pvKleyMHz8egGeeeeaM12w2G06n8+KiEqkKCudCfOe6lOOEkb9vvcEBifyh2vHdrK0zlP32+jTKPwibPoXu440OS6RcyjWM5XK5zvlQoiNSClmpsOkzAN7N70/9msG4jh8CVNkRYxW1v+CM/fjZ7czJvcLzwi9vwp+WGRHxFWVKdoYMGUJ6+h9zCp577jlOnDhR/DwtLY3WrVt7LTgRy4qfDc5cDgS3ZL27OVe2ijI6IpESHM5cujWO4DNnH/L8qkPaLvh9sdFhiZRLmZKdH374gdzc3OLnL7zwQolLzQsKCti5c6f3ohOxooJcWPs2ADNzBwE2rmwVVbxelSo7YqSi9mez2ejfKpIsglkcNNDz4prXDYxMpPzKlOz8+YewfiiLlMOWeZCVQl61aD451ZnQQD+6x0YYHZXIGa5oGQnAC8f74LbZYc8yOLLV2KBEyqFcc3ZEpJzcbljt6R2vrn0dBfjRu3kdAvzsquyIKZxe2WlSpzqxtUPY76xNcr0Bnh1Wv2FgdCLlU6Zkx2aznXE3Zt2dWaQM9v4ER7bg9q/Gy8d6ANC/VaTBQYmcW1F15yP7UM+GzZ9CZoqBEYmUXZkuPXe73dx+++0EBgYCkJOTw4QJEwgJ8dwk7vT5PCJyFms8veKUJtexcaOdYH8HV7X2TE5WZUfM4PTKDsCwDnWZtXIvM/fU4v4GXfBLWu9ZDLPfZCPDFCmTMiU7Y8eOLfH81ltvPWMf3QRU5BxSf4ffFgAwO38QAEM7xBAa5H++o0QM1b5+OK1iwtielMGK2jfSL2k9rPsv9HwQ/IOMDk+kVMqU7MyZoyXDRcrtlxkA5DcdwNydfoCLm7o3LH5ZlR0xgz9Xdmw2Gzd3b8ATX23lhf3N6RteH1v6Qdj8me6GLj5DE5RFKkPmUdjwAQBLa15PboGLFlGhdGpQw9i4REpheMd6BPnb2ZGSzaHmhdX7Va9qkUHxGUp2RCrD6tegIBt33c78+/cYAG7q3qDEBH9VdsQM/lzZAQgP9mdIO0+7fTOzl+ceWak7YftXhsQoUlZKdkQqWlZa8SKC+9rcy/bkkwT42bm2Uz2DAxMpvZsLh1w/35JBbtf/59m4/EVVd8QnKNkRqWhr3oD8LIhux1vJlwAwuG00NaoFlNhNlR0xg7NVdgC6NqpJ0zohZOc7+TpoGASEQspW2PmtEWGKlImSHZGKlH0c1r4FQE6P/+PrjUkA3NSt4fmOEjEdm81W3G7f3ZAOlxZVd/7lWSxTxMSU7IhUpF/ehNwMiGzNl9mdyMpzEls7hMuanHl7CFV2xAzOVdkBuK5zPfwdNjYfSmd741vBPwSSN8FvP1R2mCJlomRHpKLkZBQvIuju/TAfrTsIwKhuDbTyuPikWtUDGdAmGoD3NmVC9zs9L/yk6o6Ym5IdkYqy9i3ISYfazVnmiGPjwXQC/OyM7Fz/rLursiNmcL7KDsCtlzYC4PP4gyS1uRP8guHQetj9Y6XFKFJWSnZEKkJORvENP129HuLFhb8DMDauEXVCA42MTOSiXNYkgrgmtchzuvj3quPQ9Q7PC8teUHVHTEvJjkhFWPlvyD4GtZrxrSuObUkZhAb6cU/fZuc8RJUdMYMLVXZsNhsPD2oBwOfrD7K3xR2e6s7BtbD960qLU6QslOyIeNuJA8VzdQr6P81Li3cDML53E2qGBJzvSBGf0LlhTQa0jsLlhn/9nA497vO8sOhJKMgzNjiRs1CyI+JtPz4DBTnQqCefZrRjX9opaoUEcEfP2PMepsqOmMGFKjtFHhrYApsNvt+SzObGt0P1KDi+F9a9XQlRipSNkh0Rbzq0HjZ/CkBu/2f4z5JdAEy8ohnVA8t0310RU2seFVq8CvgLSw5Cv8c8Lyz/F5w6ZmBkImdSsiPiLW43/PB3z5/b38TcvTU4kpFLvRrB3HLphRcRVGVHzKC0lR2AB69sjr/DxsrfU1kVOggi20DOCfjpxQqOUqRslOyIeMuO/0HiKvAL4thljzJjuWeuzgNXXkKgn8Pg4ES8r0FENUYXXor+zx9+o+DKZzwvrH0b0nYbGJlISUp2RLyhIA8WPQGAO24iDy1M48SpfFpGh3LdOdbV+TNVdsQMylLZAbi3XzNCg/zYciiDV/Y3hGZXgisfFj9ZkWGKlImSHRFv+Hk6HNsDIZF8HHAdS3akEOBnZ/pNHXHYtVqyWFed0ECeu7YdAK8t2cWWNg+DzQ7bv9FtJMQ0lOyIXKwj2zyTMoHkuMd56odEAP42qCUto8NKfRpVdsQMylrZARjWoS7XdaqHyw0TFp4it9sEzwvfPOBZRVzEYEp2RC6GswC+uhdc+TibD+aO+MbkFrjodUltbu/R2OjoRCrN08Pb0CAimIPHs3kifThENIGTh2Hh40aHJqJkR+SirH4NDv8KQeG8Ue0etiWfpGY1f166oQP2Mg5fqbIjZlCeyg5AaJA//76xI3YbfLIxjZVtnvK88Os7sHupl6MUKRslOyLllboLlv4TgIQ2j/LyLycBeGFkeyLDgoyMTMQQXRtHMPGKSwC4+6cg0lrf7nnh6/shN9O4wKTKU7IjUh4up2f4yplLanQvrl/dGLcbbu/RmAFtost1SlV2xAzKW9kpcv8Vzbg0NoKTuQUM29Gf/NAGkJ4Ii5/yYpQiZWNosvPTTz8xdOhQ6tati81m48svv7zgMcuXL6dLly4EBQXRpEkTZs6cWfGBivzZipfhwC8U+Fdn5MFRFLhgRMe6PH5Na6MjEzGUn8PO22O70r5+OIdOOXjgVOFd0de9DTsXGBucVFmGJjtZWVl06NCB1157rVT77927lyFDhtCrVy82bNjAlClTuP/++5k3b14FRypymh3fwdJnAXgi91b2F0QwoHUU027ocFGXmauyI2ZwsZUdgLAgf975S3daRIXybVYLPncM8bzwxXg4+ps3whQpE0Nv1jN48GAGDx5c6v1nzpxJw4YNmT59OgCtWrUiPj6eadOmMXLkyAqKUuQ0R3fCF3cB8JF7IB/m9abXJbV59ZZO+Dk0KixSpGZIAO/d2Z0bZ65mctpNXFJtPx1yt8LHN8OdP0JwDaNDlCrEp346r169mgEDBpTYNnDgQOLj48nPzzcoKqkyso/j/uhmyDvJGlcrHs8dTffGEbx1W1ev3A5ClR0xA29UdopEhgbxwfjLiKwRyh2nJpJELUj73VPhcTkv+vwipeVTyU5ycjJRUVEltkVFRVFQUEBqaupZj8nNzSUjI6PEQ6TMXE4KPhuH7dhuDrprc0/eXxnaqRHv3NGd4ADd90rkXOrVCOaLe3rQqGEj7sydRLY7AHYtxPXjP4wOTaoQn0p24MzexoV6IVOnTiU8PLz40aBBgwqPUSzG5eT4J/fgt+dHst0B3F3wf/x1WBwv39jBq4mOKjtiBt6s7BSJCgvi47vi6HpZXx7N9wwD23/+NxnLSjdfU+Ri+VSyEx0dTXJycoltKSkp+Pn5UatWrbMeM3nyZNLT04sfBw4cqIxQxSKOpZ9k83+uo+bOj3G6bTzjdx9P3XUTY3s09uovAxGrC/Cz8/TwtvS74W5muK4FIGzZY/wy5xFy8goMjk6sztAJymUVFxfHN998U2LbwoUL6dq1K/7+/mc9JjAwkMDAwMoITyykwOni01U7aPzjBHqwkVy3H+/UfZwHb5lAZGjFLBioyo6YQUVUdk53baf67Ij+Dx+9H8LNWe9z6f43+exfh6hx7TSubB2tToRUCEMrO5mZmSQkJJCQkAB4Li1PSEggMdFzI8XJkyczZsyY4v0nTJjA/v37mTRpEtu3b2f27NnMmjWLhx56yIjwxYLSs/N5+6c9XPOvb2ixaAw92Eg2Qey+ajZ3/b8HKizREalKWsaEc9NDr7Gx3RQAbij4Hxkfj2fEq8uYt/4guQWavCzeZWhlJz4+nn79+hU/nzRpEgBjx45l7ty5JCUlFSc+ALGxsXz33Xc8+OCDvP7669StW5dXXnlFl53LRXG73exIPslHaxP5fP0B+hSsZq7/u0Tbj5PrH4b/rZ/TutGlFR6HKjtiBhVd2Slis9noMPJRchrVJ+B/9zHSsYKWqYlM+XwcU79vw62XNeSGrg2oVyO4QuOQqsHQZKdv377n/cE+d+7cM7b16dOHX3/9tQKjkqrA7Xaz+VA6329JZsGWZPamZlHflsJrfnO5IiABAFfNWAJv+gCi2hgbrIiFBXUdDaERuOZPoE3OfuYHPskHOf15cfEopi/eRYf64QxqG8PgttE0rh1idLjio3xqzo5IeblcbnalZLJ2bxq/7D3G2r3HSDmZiwMnXW2/MSZgPaMdPxLgzsVt98fW80HsvSaBf+X1KlXZETOorMpOCS0GY58YD4sex77xI27zW8zQgPXMyruSHw914oWDJ3hhwQ6a1A6he2xE8aN+zWqVF6P4tKqb7Gz7GqrrH4rVuIGTOQWkZuZw+EQOh05kc/h4DofST1GQl0eALZ8aFDCUPDoF7uMKx0aquTL/OLhxL2xXvwx1mhv5MUSqnup14NqZ0PEW+N8kaqTt4v/8PuP//D4jzVGbBXkd2X08hrzjfqz+1Z9lbn9CggOpWyOY+jWCqVsjmJjwIGqHBlLNX2tfWVLmqXIfWnWTna/ugUDN+rcaGxBW+Gjy5xcDznKACwiOgEsGQOth0GIIGHQ1iCo7YgaGVHZOF9sb7v4ZNn0CO7+H3UupVZDKaMdi+HMO4wTSCh9ifbnl/9lYZZOd/LqXkh9c8uMb/SvGfa4IvBhYeU/ldhfG5/kfRb+P3bj/+LP7j8/gLtwP9x/7uNxuXIX/dZ/2Z5fbM8x0+jan243T6abA5Xm4ypgABDjsBAc4qBbgoFqAH9UCHAQHBWH3C4SiR3gDaD4Q6ncDu3qCIqbhFwidx3ge+dmwdwXs/hGyUqEgB5x5OPOyOZWbx6k8p+eRW0BOvpN8V+l/Vths4GezYbfbcNhtOP70Z4fdhs0GdpvnvzZs2G2eRLB4e+F5Tt+nKE+0Ff+fZ3vRe57+mu2PZyVf4/wJ5wVT0QvsUBGpbEUnyM7sfGBxuY6tsslOpz13YQ/UMJavCQlwEB7sT1iwP1FhQUSGBhIVFkRUWCD1I6rRMKIa9WsGe+VeVZVNlR0xA8MrO3/mHwzNB3gep3EAoYWP02XlFnDg+Cn2p53i8Ils0jLzSM3MJTUzj7SsXFIzc0nL9CRJ4ltcuadQsuMl5/r3fb5/9uf6oXCuY873M8R2rqPKFVfZ3uN8n/2Mnk9hr6foYbdR+N/Tt/1xTICfnUA/O0H+DgL97AT62wn0K/yzn51A/9P+7OcguDCpKUpswoP9CQvy053FReS8QgL9aBkdRsvosPPudyqvgLTMPLLyCsjKLSAr1+n5b17Rfz3bc/Nd5Dtd5Dnd5DtdxY+8gj89d7rJK/D82e0+WzW76LlnW4nqtstdXDE/1z6nU3+o7KpssrPl6YGEhZ3/H4NIZVJlR8zAdJWdClItwI9qEdb8FfjnnyF//pFytp8wZxxzxut/PkfZEjBv/FjLyMggenr5jrXm37SIiEgV9edEtXR5q/mT2/yLuPGyxgRETEKVHTGDqlLZkapFyY6IiIhYmpIdEZNQZUfMQJUdsSIlOyIiImJpSnZETEKVHTEDVXbEipTsiIiIiKUp2RExCfWkxUzUHsVKlOyImIyGscRIan9iRUp2RExCPWkxE7VHsRIlOyImo561GEntT6xIyY6ISagnLWai9ihWomRHxGTUsxYjqf2JFSnZETEJ9aTFTNQexUqU7IiYjHrWYiS1P7EiJTsiJqGetJiJ2qNYiZIdEZNRz1qMpPYnVqRkR8Qk1JMWM1F7FCtRsiNiMupZi5HU/sSKlOyImIR60mImao9iJUp2RExGPWsxktqfWJGSHRGTUE9azETtUaxEyY6IyahnLUZS+xMrUrIjYhLqSYuZqD2KlSjZETEZ9azFSGp/YkVKdkRMQj1pMRO1R7ESJTsiJqOetRhJ7U+sSMmOiEmoJy1movYoVqJkR8Rk1LMWI6n9iRUp2RExCfWkxUzUHsVKlOyImIx61mIktT+xIiU7IiIiYmlKdkRMomjYQD1rMVJR+9MwlliJkh0RERGxNCU7Iiahyo6YgSo7YkVKdkRERMTSlOyImIQqO2IGquyIFSnZEREREUtTsiNiEqrsiBmosiNWpGRHRERELE3JjohJqLIjZqDKjliRkh0RERGxNMOTnTfeeIPY2FiCgoLo0qULK1asOOe+y5Ytw2aznfHYsWNHJUYsUjFU2REzUGVHrMjQZOeTTz7hgQce4LHHHmPDhg306tWLwYMHk5iYeN7jdu7cSVJSUvHjkksuqaSIRURExNcYmuy8/PLLjBs3jjvvvJNWrVoxffp0GjRowIwZM857XGRkJNHR0cUPh8NRSRGLVBxVdsQMVNkRKzIs2cnLy2P9+vUMGDCgxPYBAwawatWq8x7bqVMnYmJi6N+/P0uXLj3vvrm5uWRkZJR4iIiISNVhWLKTmpqK0+kkKiqqxPaoqCiSk5PPekxMTAxvvfUW8+bN44svvqBFixb079+fn3766ZzvM3XqVMLDw4sfDRo08OrnEPEWVXbEDFTZESvyMzqAP/+Dcrvd5/xH1qJFC1q0aFH8PC4ujgMHDjBt2jR69+591mMmT57MpEmTip9nZGQo4REREalCDKvs1K5dG4fDcUYVJyUl5Yxqz/lcdtll7Nq165yvBwYGEhYWVuIhYkaq7IgZqLIjVmRYshMQEECXLl1YtGhRie2LFi2iR48epT7Phg0biImJ8XZ4IiIiYhGGDmNNmjSJ2267ja5duxIXF8dbb71FYmIiEyZMADxDUIcOHeLdd98FYPr06TRu3Jg2bdqQl5fH+++/z7x585g3b56RH0PEK1TZETNQZUesyNBkZ9SoUaSlpfHMM8+QlJRE27Zt+e6772jUqBEASUlJJdbcycvL46GHHuLQoUMEBwfTpk0bvv32W4YMGWLURxARERGTM3yC8j333MM999xz1tfmzp1b4vkjjzzCI488UglRiVQ+VXbEDFTZESsy/HYRIiIiIhVJyY6ISaiyI2agyo5YkZIdERERsTQlOyImocqOmIEqO2JFSnZERETE0pTsiJiEKjtiBqrsiBUp2RERERFLU7IjYhKq7IgZqLIjVqRkR0RERCxNyY6ISaiyI2agyo5YkZIdERERsTQlOyImocqOmIEqO2JFSnZERETE0pTsiJiEKjtiBqrsiBUp2RERERFLU7IjYhKq7IgZqLIjVqRkR0RERCxNyY6ISaiyI2agyo5YkZIdERERsTQlOyImocqOmIEqO2JFSnZERETE0pTsiJiEKjtiBqrsiBUp2RERERFLU7IjYhKq7IgZqLIjVqRkR0RERCxNyY6ISaiyI2agyo5YkZIdERERsTQlOyImocqOmIEqO2JFSnZERETE0pTsiJiEKjtiBqrsiBUp2RERERFLU7IjYhKq7IgZqLIjVqRkR0RERCxNyY6ISaiyI2agyo5YkZIdERERsTQlOyImocqOmIEqO2JFSnZERETE0pTsiJiEKjtiBqrsiBUp2RERERFLU7IjYhKq7IgZqLIjVqRkR0RERCxNyY6ISaiyI2agyo5YkZIdERERsTQlOyImocqOmIEqO2JFSnZERETE0pTsiJiEKjtiBqrsiBUp2RERERFLU7IjYhKq7IgZqLIjVqRkR0RERCzN8GTnjTfeIDY2lqCgILp06cKKFSvOu//y5cvp0qULQUFBNGnShJkzZ1ZSpCIVS5UdMQNVdsSKDE12PvnkEx544AEee+wxNmzYQK9evRg8eDCJiYln3X/v3r0MGTKEXr16sWHDBqZMmcL999/PvHnzKjlyERER8RV+Rr75yy+/zLhx47jzzjsBmD59Oj/88AMzZsxg6tSpZ+w/c+ZMGjZsyPTp0wFo1aoV8fHxTJs2jZEjR1Zm6CJeV9ST/v333/npp58MjkaqqiNHjgCq7Ii1GJbs5OXlsX79ev72t7+V2D5gwABWrVp11mNWr17NgAEDSmwbOHAgs2bNIj8/H39//zOOyc3NJTc3t/h5RkaGF6IX8T673VNofe2113jttdcMjkaquqL2KGIFhiU7qampOJ1OoqKiSmyPiooiOTn5rMckJyefdf+CggJSU1OJiYk545ipU6fy9NNPey9wkQpy++23s3XrVnJycowORaq4iIgIhg4danQYIl5j6DAWnFkqdbvd5y2fnm3/s20vMnnyZCZNmlT8PCMjgwYNGpQ3XJEKc+WVV7JhwwajwxARsRzDkp3atWvjcDjOqOKkpKScUb0pEh0dfdb9/fz8qFWr1lmPCQwMJDAw0DtBi4iIiM8xbFA2ICCALl26sGjRohLbFy1aRI8ePc56TFxc3Bn7L1y4kK5du551vo6IiIiIoTPQJk2axH//+19mz57N9u3befDBB0lMTGTChAmAZwhqzJgxxftPmDCB/fv3M2nSJLZv387s2bOZNWsWDz30kFEfQUREREzO0Dk7o0aNIi0tjWeeeYakpCTatm3Ld999R6NGjQBISkoqseZObGws3333HQ8++CCvv/46devW5ZVXXtFl5yIiInJONncVW641IyOD8PBw0tPTCQsLMzocERERKYWL+f2thRRERETE0pTsiIiIiKUp2RERERFLU7IjIiIilmb4CsqVrWg+tu6RJSIi4juKfm+X57qqKpfspKWlAeiWESIiIj4oLS2N8PDwMh1T5ZKdiIgIABITE8v8ZUlJRfcZO3DggC7jv0j6Lr1D36P36Lv0Hn2X3pGenk7Dhg2Lf4+XRZVLdux2zzSl8PBwNTovCQsL03fpJfouvUPfo/fou/QefZfeUfR7vEzHVEAcIiIiIqahZEdEREQsrcolO4GBgTz55JMEBgYaHYrP03fpPfouvUPfo/fou/QefZfecTHfY5W7N5aIiIhULVWusiMiIiJVi5IdERERsTQlOyIiImJpSnZERETE0qp8sjNs2DAaNmxIUFAQMTEx3HbbbRw+fNjosHzKvn37GDduHLGxsQQHB9O0aVOefPJJ8vLyjA7NJz333HP06NGDatWqUaNGDaPD8SlvvPEGsbGxBAUF0aVLF1asWGF0SD7np59+YujQodStWxebzcaXX35pdEg+aerUqXTr1o3Q0FAiIyMZMWIEO3fuNDosnzRjxgzat29fvChjXFwc33//fZnOUeWTnX79+vHpp5+yc+dO5s2bx+7du7n++uuNDsun7NixA5fLxZtvvsnWrVv597//zcyZM5kyZYrRofmkvLw8brjhBu6++26jQ/Epn3zyCQ888ACPPfYYGzZsoFevXgwePJjExESjQ/MpWVlZdOjQgddee83oUHza8uXLuffee1mzZg2LFi2ioKCAAQMGkJWVZXRoPqd+/fo8//zzxMfHEx8fzxVXXMHw4cPZunVrqc+hS8//5Ouvv2bEiBHk5ubi7+9vdDg+68UXX2TGjBns2bPH6FB81ty5c3nggQc4ceKE0aH4hEsvvZTOnTszY8aM4m2tWrVixIgRTJ061cDIfJfNZmP+/PmMGDHC6FB83tGjR4mMjGT58uX07t3b6HB8XkREBC+++CLjxo0r1f5VvrJzumPHjvHBBx/Qo0cPJToXKT09vVw3axMpj7y8PNavX8+AAQNKbB8wYACrVq0yKCqRP6SnpwPo5+JFcjqdfPzxx2RlZREXF1fq45TsAI8++ighISHUqlWLxMREvvrqK6ND8mm7d+/m1VdfZcKECUaHIlVEamoqTqeTqKioEtujoqJITk42KCoRD7fbzaRJk+jZsydt27Y1OhyftHnzZqpXr05gYCATJkxg/vz5tG7dutTHWzLZeeqpp7DZbOd9xMfHF+//8MMPs2HDBhYuXIjD4WDMmDFodK/s3yPA4cOHGTRoEDfccAN33nmnQZGbT3m+Syk7m81W4rnb7T5jm0hlmzhxIps2beKjjz4yOhSf1aJFCxISElizZg133303Y8eOZdu2baU+3q8CYzPMxIkTuemmm867T+PGjYv/XLt2bWrXrk3z5s1p1aoVDRo0YM2aNWUqkVlRWb/Hw4cP069fP+Li4njrrbcqODrfUtbvUsqmdu3aOByOM6o4KSkpZ1R7RCrTfffdx9dff81PP/1E/fr1jQ7HZwUEBNCsWTMAunbtyrp16/jPf/7Dm2++WarjLZnsFCUv5VFU0cnNzfVmSD6pLN/joUOH6NevH126dGHOnDnY7ZYsGpbbxbRJubCAgAC6dOnCokWLuPbaa4u3L1q0iOHDhxsYmVRVbreb++67j/nz57Ns2TJiY2ONDslS3G53mX5PWzLZKa21a9eydu1aevbsSc2aNdmzZw9PPPEETZs2rfJVnbI4fPgwffv2pWHDhkybNo2jR48WvxYdHW1gZL4pMTGRY8eOkZiYiNPpJCEhAYBmzZpRvXp1Y4MzsUmTJnHbbbfRtWvX4upiYmKi5o6VUWZmJr///nvx871795KQkEBERAQNGzY0MDLfcu+99/Lhhx/y1VdfERoaWlx1DA8PJzg42ODofMuUKVMYPHgwDRo04OTJk3z88ccsW7aMBQsWlP4k7ips06ZN7n79+rkjIiLcgYGB7saNG7snTJjgPnjwoNGh+ZQ5c+a4gbM+pOzGjh171u9y6dKlRodmeq+//rq7UaNG7oCAAHfnzp3dy5cvNzokn7N06dKztr+xY8caHZpPOdfPxDlz5hgdms+54447iv9d16lTx92/f3/3woULy3QOrbMjIiIilqaJFSIiImJpSnZERETE0pTsiIiIiKUp2RERERFLU7IjIiIilqZkR0RERCxNyY6IiIhYmpIdERERsTQlOyIiImJpSnZERETE0pTsiIjPO3r0KNHR0fzzn/8s3vbLL78QEBDAwoULDYxMRMxA98YSEUv47rvvGDFiBKtWraJly5Z06tSJq6++munTpxsdmogYTMmOiFjGvffey+LFi+nWrRsbN25k3bp1BAUFGR2WiBhMyY6IWEZ2djZt27blwIEDxMfH0759e6NDEhET0JwdEbGMPXv2cPjwYVwuF/v37zc6HBExCVV2RMQS8vLy6N69Ox07dqRly5a8/PLLbN68maioKKNDExGDKdkREUt4+OGH+fzzz9m4cSPVq1enX79+hIaG8r///c/o0ETEYBrGEhGft2zZMqZPn857771HWFgYdrud9957j5UrVzJjxgyjwxMRg6myIyIiIpamyo6IiIhYmpIdERERsTQlOyIiImJpSnZERETE0pTsiIiIiKUp2RERERFLU7IjIiIilqZkR0RERCxNyY6IiIhYmpIdERERsTQlOyIiImJpSnZERETE0v4/cqOC/O1edrsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "zoom = 2.0 # zoom factor for wavefunctions to make them more visible\n",
    "plt.plot(x,V(x),'-k',label=\"V(x)\") # plot the potential\n",
    "plt.plot(x,zoom*psi_exact(x)+E_exact,label=\"Exact\")\n",
    "#plt.plot(x,zoom*np.sign(wavefunctions[n//2,0])*wavefunctions[:,0]+energies[0],label=\"Numerov\") #plot the num-th wavefunction\n",
    "plt.plot(x,zoom*psi(x,result.x)+result.fun,label=\"variational\") # plot the variational wavefunction\n",
    "plt.ylim(-0.1,3); # set limits of vertical axis for plot\n",
    "plt.xlim(-3,3); # set limits of horizontal axis for plot\n",
    "plt.legend();\n",
    "plt.xlabel(\"x\");\n",
    "plt.ylabel(\"Energy\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "Python (scipy)",
   "language": "python",
   "name": "scipy"
  },
  "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.10.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}