diff --git a/exercises/ex09_solution/dft.ipynb b/exercises/ex09_solution/dft.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..4d37f4427fea2fe335a2b3f72498e7efa0ebb6e3
--- /dev/null
+++ b/exercises/ex09_solution/dft.ipynb
@@ -0,0 +1,473 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.cm as cm\n",
+    "from numpy import sqrt, exp, sinh, conj, absolute, pi\n",
+    "from scipy import constants, special\n",
+    "plt.rcParams['figure.figsize'] = 16, 9"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We use atomic units; all wave functions are represented such that $u[i]$ is the value at the $i$-th bin; normalization then corresponds to $\\lVert u \\rVert_{L^2}^2 = dx \\sum_i u[i]^2 = 1$."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# (a) Schrödinger Solver"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Numerov Integration (from solution to exercise 3.2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Integrate $-\\psi''(x) = k(x) \\psi(x)$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "def numerov_step(k0, k1, k2, psi0, psi1, dx):\n",
+    "    \"\"\"compute psi2 via a single numerov step\"\"\"\n",
+    "    dx_sqr = dx**2\n",
+    "    c0 = (1 + 1/12. * dx_sqr * k0)\n",
+    "    c1 = 2 * (1 - 5/12. * dx_sqr * k1)\n",
+    "    c2 = (1 + 1/12. * dx_sqr * k2)\n",
+    "    return (c1 * psi1 - c0 * psi0) / c2\n",
+    "\n",
+    "def numerov(ks, psi0, psi1, dx):\n",
+    "    \"\"\"compute psi = [psi0, psi1, psi2, ...] for ks = [k0, k1, k2, ...] via iterated numerov steps\"\"\"\n",
+    "    psi = np.zeros_like(ks)\n",
+    "    psi[0] = psi0\n",
+    "    psi[1] = psi1\n",
+    "    for i in range(2, len(ks)):\n",
+    "        psi[i] = numerov_step(ks[i-2], ks[i-1], ks[i], psi[i-2], psi[i-1], dx)\n",
+    "    return psi"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Radial Schrödinger Solver (adapted from solution to exercise 2.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Find ground state for radial Schrödinger equation $-\\frac 1 2 u''(r) + V(r) u(r) = E u(r)$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "01: trying E = -5.000000 in [-10.000000, 0.000000]\n",
+      ", u(0) = 12884645614632673280.000000, #nodes = 0\n",
+      "02: trying E = -2.500000 in [-5.000000, 0.000000]\n",
+      ", u(0) = 65058210171.040329, #nodes = 0\n",
+      "03: trying E = -1.250000 in [-2.500000, 0.000000]\n",
+      ", u(0) = 54164.193034, #nodes = 0\n",
+      "04: trying E = -0.625000 in [-1.250000, 0.000000]\n",
+      ", u(0) = 0.750139, #nodes = 0\n",
+      "05: trying E = -0.312500 in [-0.625000, 0.000000]\n",
+      ", u(0) = -0.000956, #nodes = 1\n",
+      "06: trying E = -0.468750 in [-0.625000, -0.312500]\n",
+      ", u(0) = -0.007805, #nodes = 1\n",
+      "07: trying E = -0.546875 in [-0.625000, -0.468750]\n",
+      ", u(0) = 0.063054, #nodes = 0\n",
+      "08: trying E = -0.507812 in [-0.546875, -0.468750]\n",
+      ", u(0) = 0.005109, #nodes = 0\n",
+      "09: trying E = -0.488281 in [-0.507812, -0.468750]\n",
+      ", u(0) = -0.004283, #nodes = 1\n",
+      "10: trying E = -0.498047 in [-0.507812, -0.488281]\n",
+      ", u(0) = -0.000560, #nodes = 1\n",
+      "11: trying E = -0.502930 in [-0.507812, -0.498047]\n",
+      ", u(0) = 0.001995, #nodes = 0\n",
+      "12: trying E = -0.500488 in [-0.502930, -0.498047]\n",
+      ", u(0) = 0.000652, #nodes = 0\n",
+      "13: trying E = -0.499268 in [-0.500488, -0.498047]\n",
+      ", u(0) = 0.000031, #nodes = 0\n",
+      "  -> converged after 13 iterations: E_0 = -0.499267578125\n"
+     ]
+    }
+   ],
+   "source": [
+    "r_min = 0\n",
+    "r_max = 20\n",
+    "r_steps = 50000\n",
+    "rs, dr = np.linspace(r_min, r_max, r_steps + 1, retstep=True)\n",
+    "r_min, rs = rs[1], rs[1:]  # skip r=0\n",
+    "\n",
+    "def solve_schrodinger(Vs, E_min=-10, E_max=0, maxiter=100, stoptol = 0.0001, verbose=False, retnumiter=False):\n",
+    "    \"\"\"returns ground state energy, ground state wave function, and -- optionally -- the number of numerov iterations that were required\"\"\"\n",
+    "    for i in range(maxiter):\n",
+    "        E = (E_min + E_max) / 2.\n",
+    "        if verbose:\n",
+    "            print(\"%02d: trying E = %f in [%f, %f]\" % (i + 1, E, E_min, E_max),)\n",
+    "        \n",
+    "        # numerov integration starting from r*exp(-r) at r_max to r_min\n",
+    "        ks = 2 * (E - Vs)\n",
+    "        psi_last = rs[-1] * exp(-rs[-1])\n",
+    "        psi_2nd_last = rs[-2] * exp(-rs[-2])\n",
+    "        u = numerov(ks[::-1], psi_last, psi_2nd_last, dr)\n",
+    "        u = u[::-1]\n",
+    "        num_nodes = np.sum(u[1:] * u[:-1] < 0)\n",
+    "        if verbose:\n",
+    "            print(', u(0) = %f, #nodes = %d' % (u[0], num_nodes))\n",
+    "        \n",
+    "        # bisect       \n",
+    "        if num_nodes == 0 and absolute(u[0]) <= stoptol:\n",
+    "            if verbose:\n",
+    "                print(\"  -> converged after\", i + 1, \"iterations: E_0 =\", E)\n",
+    "            return (E, u, i+1) if retnumiter else (E, u)\n",
+    "        if num_nodes == 0 and u[0] > 0:\n",
+    "            E_min = E\n",
+    "        else:\n",
+    "            assert num_nodes > 0, \"expect #nodes > 0 since u[0] < 0 while u[infty] > 0\"\n",
+    "            E_max = E\n",
+    "    \n",
+    "    raise Exception(\"Not converged after %d iterations; u(0) = %f\" % (maxiter, u[0]))\n",
+    "\n",
+    "E, u = solve_schrodinger(-1./rs, verbose=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compare with exact solution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "E_exact = 0.5\n",
+    "u_exact = rs * exp(-rs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7QAAAIyCAYAAAD7WYbfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl0VfW9///XOyNJCBAIYZ4RFFQQFcUBAwRQFOcBtLa1\nxdpBq/5stWp7RWu9rb1+r1rtoLe2OFJn0aIMapwVUEFFJgWBkDCEzEDmz++Pc6CHmEAScrLP2ef5\nWCurZ9jn83mdZHW5Xnw+e29zzgkAAAAAgGgT53UAAAAAAABag0ILAAAAAIhKFFoAAAAAQFSi0AIA\nAAAAohKFFgAAAAAQlSi0AAAAAICoRKEFAAAAAEQlCi0AAAAAICpRaAEAEcHM/mFmd3idI9KZ2TAz\n+9TMSs3s6nac9wszGx+mse8ys5+HY+xG5vrIzI5oj7kAAOGX4HUAAEB0MrMNkn7onHsj5LXvSZrl\nnDvVu2S+d6OkN5xzx4RzkoZ/X+fckWGaJ1PS5ZKGtvG450gaKalOUr5z7rHgW3+U9FtJF7blfAAA\nb7BCCwBoa64tBzOz+LYczwcGSFrpdYg29H1J851zVXtfMLOZZrajtQOaWSdJ/+Wcu8s59wdJPzWz\nbsG3X5Y0wcyyDiU0ACAyUGgBAGFhZr8ws2cbvHa/mf1v8PExZvZxcOvsXEkdQo7bYGY3mtkKSRVm\nFmdmR5jZm2ZWbGafm9n0BmOPMbNPguM9bWZz925hNrNeZvasmW03s6/N7JoGn91gZjeY2Yrg+E+Z\nWVIj3+n7ZjYv5Pk6M/tXyPNNZnZ08PFNZvaVmZUFt+ueG3z9RjN7psG495nZvQfLamavS5og6cHg\nuIeZWb2ZDQ45Zt/W7YN9LzPra2bPBefaYWb3B19/VFJ/Sa8E5/llcKyJIZ89vKm/R3N/n0FnSHqr\nwWurJOU2cXxzjNf+pX+FAr83BYvzx5KmHsL4AIAIQaEFALQlC3n8uKSpwdWyvSutl0iaY2aJkl6Q\nNEdSV0nPSLqgwVgzFCg7XRT479U8Sa9J6i7p55KeMLPDgmMnSnpe0iPB8Z6SdF7wPVNgVe5TSb0k\nTZJ0rZlNbjDfRZKmSBokaZQCK4cNvSXplOC4vSQlShoXfD5YUppz7rPgsV9JOtk510nS7ZIeN7Me\nkuZKOsPM0oKfiwvO/cTBsjrnJkl6R9LPnHOdnHPrdPAV8Ua/V3DeVyRtUKC89glmk3Puu5I2SToz\nOM8fQwc0s4Rgzkb/Hi34fUrSUZLWNHhtoqQ3Ql8ws8Fm9t/B823/u8Hju8zs7JDD+0oqCXleIik0\n26pgJgBAlOMcWgDAoXjRzGpDnicrsPol59xWM3tHgWLzdwXK6Q7n3HILXFwowTl3f/Bzz5nZ0gZj\n3+ecy5ckMztFgbL4h+B7b5rZK5JmSrpD0omS4p1zDwTff8HMlgQfHy8p0zn3u+Dzb8zs/4KfXdRg\nvm3B+V6WNLrhl3XObTCzcjMbLWm4pAWSRpnZMEknKVA29x77XMjjZ8zsFkljnXMvm9knChTuxxUo\nrbucc0vN7IRmZg1lTbx+sO91ggKl+UbnXH3wtfebOfaJOvDf40DzNtRFUnmD1yZI+kXoC8659ZJu\nbmKMhjIkVYY8r5bUMeR5uaSezRwLABDBWKEFAByKc5xzXff+SPppg/cflfSd4OPLJO29ME8vSVsa\nHLuxwfO8kMe9JW1u5Pg+Ie83HG/v8QMk9TGzouBPsQLFqHuD47eFPN6t/QtQqLcUKFzjFdgWmysp\nW9JpCtk6a2bftcDViIuDc46UlBl8+ykFyp+C//tk8HH/ZmZtiaa+V19JG0PKbEsc7O9xoHkbKpaU\nvvdJcOV4iHOu4aptS5Rr/zKeIqko5Hm69l/BBQBEKVZoAQCH4mCrgy8qcL7nSElnSfpl8PUC7V9+\npECZ+yrkeehW2nxJ/Ro5fm/paWy8fsHxNkta75wbfpCszfW2pOmSBkr6naRSBcr6iZL+JElm1l/S\nQ5ImOOc+CL72qf7z+3pG0v+YWR8FVmpPDL7emqy7JaWGPO+pb5fNxmyW1N/M4pootQfaynywv0dL\nfCZpmIIr+wqsqC9reFBwS/eVjeSy4GsfOuf2nt/8taTjQo7pJumTkOdH6D//uAIAiGKs0AIAwsY5\nV6nAua1PSvrIObd31fUDSbVmdo2ZJZjZ+ZLGHmCojyTtDl5QKcHMshUoyHNDxqszs5+ZWbwFbtmy\nd7wlksqDn+0QfH+kmR33rVmaZ+8KbUpwS/Q7kk5XoDR9GjwmTVK9pEILXNDqCkn7bnvjnCsMjvMP\nBQrs3iLYmqyfSro0OM/pCqwUN8cSBf4h4PdmlmpmyWZ2Usj7WyUNbvyjTf49nmrm3KHmK7DCvdd4\nBbYwnxt6kHNuvXPuZufcLQ1+9r42L+TwtySNCXk+RtLrkmRmyZKOVdNbuAEAUYRCCwBorebenmeO\nAhf+eXTfB52rkXS+pCsk7VTgPNvnQj6z39jB46dLmiapUNIDki53zq1tMN4sBbawXqrARYuqgquP\nZylwDucGSdslPSypUyu+i4IXYipXYKVWzrlyBVYE33XOueBrqyTdI+lDBYrhSEnvNhjqSQXOn30i\nZOzWZL1O0tnB7z1TgYttHfR7BeearsDFkjYpsGJ7ccghv5f0m+DW5xtCxzrA32PdweZtxKMKXCQr\nOfh8rQLboYtbMMZ+nHO7Jd1tZr82s99I+qNzbnvw7bMlvemc29ra8QEAkcOC/+0N3wSBfy2+V4Hy\n/PeQC0jsfb+bAhfF6CUpXtI9zrl/hjUUAKDdmFlfSasl9XTOVbTjvB9K+otzbk57zYnWMbM7JW0P\nuUhYOOf6QNIPnXNfhnsuAED4hbXQBi/ssFaBf4HOl7RU0gzn3OqQY26T1ME5d7OZZSpw/k0P51xt\nY2MCAKJH8L8D/09SR+fcrDDPNV6B/4YUKnAhqj9LGrz3SrsAAMB/wn1RqLGS1jnnNkqSmc2VdI4C\n/1K/11YFtqJJgasO7qTMAkD0M7NUBa50u0GBW/aE23BJTytwgaT1ki6gzAIA4G/hLrR9tP+VFvP0\n7Yt+PCzpdTPLV+CS/peEORMAoB0Ez2NMP+iBbTffwwr8NwUAAMSISLgo1M2SVjjneks6RoHbOzR1\nrzoAAAAAACSFf4V2iwL3pdurr7594/uTFbiPn5xzX5vZBkmHq8E96MwsvFevAgAAAAB4yjl3sHvc\n7yfchXappKFmNkCBe93NUOCWAqFWScqR9J6Z9VDg5urrGxss3FdkBvBts2fP1uzZs72OAcQk/v8H\neIP/7wHeMGtRl5UU5kLrnKszs6slLdR/btuzysyuCrztHpL035L+YWYrJJmkG51zReHMBQAAAACI\nfuFeoZVz7jUFrjwZ+trfQh4XKnBzdgAAAAAAmi0SLgoFIIJlZ2d7HQGIWfz/D/AG/98DoodFy3mp\nZuaiJSsAAAAAoGXMLOIuCgUAAAAAEWPgwIHauHGj1zFi2oABA/TNN9+0yVis0AIAAACIGcFVQK9j\nxLSm/gatWaHlHFoAAAAAQFSi0AIAAAAAohKFFgAAAAAQlSi0AAAAAICoRKEFAAAAAEQlCi0AAAAA\noEl1dXVeR2gShRYAAAAAIsSgQYN0zz33aNSoUcrIyNDMmTNVVVWlOXPm6NRTT93v2Li4OK1fv16S\ndMUVV+hnP/uZpk2bpvT0dI0fP15bt27Vddddp4yMDI0YMUIrVqzY99mCggJdeOGFysrK0pAhQ/Sn\nP/1p33u33367LrroIl1++eXq0qWL5syZo+rqal133XXq06eP+vbtq+uvv141NTWSpBEjRmj+/Pn7\nPl9XV6esrCwtX748nL+qwO8g7DMAAAAAAJrtmWee0cKFC7VhwwatWLFCc+bMkRS4T2uohs+feeYZ\n3XXXXdq5c6cSExN14okn6vjjj1dRUZEuuOACXX/99ZIk55ymT5+uY445RgUFBXr99dd13333adGi\nRfvGmjdvni6++GKVlJTo0ksv1Z133qklS5bos88+04oVK7RkyRLdeeedkqSZM2fqySef3PfZ1157\nTd27d9fo0aPD8vsJRaEFAAAAgCCztvk5FNdee6169OihLl26aPr06U2udDrn9nt+3nnnafTo0UpK\nStJ5552ntLQ0XXbZZTIzXXLJJfvGWbJkiQoLC3XrrbcqPj5eAwcO1KxZszR37tx9Y40bN07Tp0+X\nJHXo0EFPPvmkbrvtNnXr1k3dunXTbbfdpkcffVSSdOmll2revHmqrKyUJD311FOaOXPmof0Smimh\nXWYBAAAAgCjQoCN6okePHvsep6amqqCgoMWfS0lJ+dbziooKSdKmTZu0ZcsWde3aVVKgGNfX12v8\n+PH7ju/Xr99+Y+fn56t///77ng8YMGBfriFDhmjEiBF6+eWXddZZZ2nevHm64447mvt1DwmFFgAA\nAAAiXFpamnbt2rXv+datW1s9Vr9+/TR48GCtWbOmyWMabmfu06ePNm7cqCOOOEKStHHjRvXu3Xvf\n+zNmzNCTTz6puro6jRw5UoMHD251vpZgyzEAAAAARLhRo0bpyy+/1Geffaaqqirdfvvt3yqdB7N3\ni/LYsWOVnp6uu+++W5WVlaqrq9PKlSu1bNmyJj87Y8YM3XnnnSosLFRhYaF++9vf6vLLL9/v/YUL\nF+ovf/mLLr300tZ9yVag0AIAAABAhGiqpB522GH6zW9+o0mTJmnYsGHfuuJxS8aOi4vTK6+8ouXL\nl2vQoEHKysrSlVdeqbKysiY/++tf/1rHHXecjj76aI0aNUrHHXecbr311n3v9+zZU+PGjdOHH36o\nSy65pMXZWssankgcqczMRUtWAAAAAJHJzL51MSW0r6b+BsHXW7TszAotAAAAACAqUWgBAAAAAFGJ\nQgsAAAAAiEoUWgAAAABAVKLQAgAAAACiEoUWAAAAABCVKLQAAAAAgKhEoQUAAAAARCUKLQAAAAAg\nKlFoAQAAACAGzZkzR6eeeqrXMQ4JhRYAAAAAYpBzTmbmdYxDQqEFAAAAgAhRUFCgCy+8UFlZWRoy\nZIgeeOABSdKZZ56pX/ziF/uOmzFjhmbNmiVJWr9+vSZNmqTMzExlZWXpO9/5jsrKyvYdm5eXpwsu\nuEBZWVnq3r27fv7zn2v16tX6yU9+og8++EDp6enq2rVr+37RNkKhBQAAAIAI4JzT9OnTdcwxx6ig\noECvv/667r33Xi1atEiPPPKIHn/8ceXm5uqJJ57QsmXLdP/99+/73C233KKtW7dq1apVysvL0+zZ\nsyVJ9fX1OuusszRo0CBt2rRJW7Zs0YwZM3T44Yfrr3/9q8aNG6fy8nIVFRV5+M1bL8HrAAAAAAAQ\nKez2ttmC625zLf7M0qVLVVhYqFtvvVWSNHDgQM2aNUtz587V5MmT9Ze//EXf/e53VVlZqZdeekmp\nqamSpCFDhmjIkCGSpG7duun666/XHXfcIUn66KOPVFBQoLvvvltxcYH1zJNOOqktvmJEoNACAAAA\nQFBrimhb2bhxo7Zs2bJv+69zTvX19Ro/frwk6ayzztLVV1+t4cOHa9y4cfs+t337dl177bV65513\nVFFRobq6un1j5OXlacCAAfvKrN/481sBAAAAQJTp16+fBg8erKKiIhUVFam4uFilpaV6+eWXJUm3\n3HKLRowYoYKCAs2dO3ff52655RbFxcVp5cqVKikp0eOPPy7n3L4xN23apPr6+m/NF+0XhJIotAAA\nAAAQEcaOHav09HTdfffdqqysVF1dnVauXKlly5bp7bff1pw5c/TYY4/pn//8p6655hoVFBRIksrL\ny9WxY0elp6dry5Yt+uMf/7jfmL169dKvfvUr7d69W1VVVXr//fclST169FBeXp5qamo8+b5tgULb\nxiqqK/TgO0/ooodu0uWP3KbnPl2kevftfw0BAAAAgFBxcXF65ZVXtHz5cg0aNEhZWVm68sorVVBQ\noO9///t68MEH1bNnT51yyimaNWuWrrjiCknSbbfdpo8//lhdunTR9OnTdcEFF+w35ssvv6x169ap\nf//+6tevn55++mlJ0sSJEzVy5Ej17NlTWVlZnnznQ2V7l6IjnZm5SM/68LvP6+ev/Uy1m4/VgLiT\nVGO7lddhvjp1qdGjFz+k6aPHHXwQAAAAAGFjZor0XuF3Tf0Ngq+3aB80hbaN3PDM/+q+j+7T91Kf\n1p9+NVbBC46puNjpB//zrObVXqNfjLlLf7jkB94GBQAAAGIYhdZ7FNoIc+/iubrhtZv01+Pf05WX\n9G30mH/MW6tZb03Vdcf8Rvd8h1ILAAAAeIFC6722LLTctucQrdm+Qb984xrd2G9xk2VWkq44e5iS\nOyzQ5YuzNfDV3rrmjNPbMSUAAAAA+A8rtIfAOafBs3OUkn+6Vj70SzXnqtd3zHlbt6+6WMt/ulRH\n9e8X/pAAAAAA9mGF1nttuULLVY4Pwd/ffUWbi7fqtduub1aZlaT/+t54neCu1cQ/Xa66Ru4FBQAA\nAABoHgptK9XV1+mGV2/SRRl3q3/flu3cXjz7RlVU7dHP/u/vYUoHAAAAAP7HObSt9Oc3X9Tu4s56\n6O5pLf5sakq8/jrtYf0gd5Kuy5uuw/v2DENCAAAAAA0NGDBA1tztlQiLAQMGtNlYnEPbSj1vPVmn\nxF+vZ++4sNVjjPrlL1SfWKbP73qoDZMBAAAAQPThHNp2suCLj7RjT4Hu+/F5hzTO09fcqpV1L+rV\nj1e2UTIAAAAAiB0U2lb4r3kPa3Ttj9Wnd/whjTO8f4ZOT7tZP3zqpjZKBgAAAACxg0LbQruqd+vj\nXc/r5jO/0ybjPXHdT7Vdn+uxN5a0yXgAAAAAECsotC1034IXlbT9BJ0/uXebjJfRKVnndL9RN/37\nd20yHgAAAADECgptCz300aM6o/d3FdeGv7m/XfUDbYtfouffX9F2gwIAAACAz1FoW6B4T4k2ufd1\n60Vnt+m4mV1SNLXTDbr++bvadFwAAAAA8DMKbQs8uHC+Urdna8yRaW0+9sM/vkqbExbr/S83tvnY\nAAAAAOBHFNoWePKTl3Rq1jlhGbtPZrpG2/d0/ZMPhmV8AAAAAPAbCm0zVdVWaU3tAl0zZXrY5rjn\nkqu1tPYR7SjZFbY5AAAAAMAvwl5ozex0M1ttZmvN7Fs3XDWzX5jZp2b2iZl9bma1ZtYl3Lla6qkP\n3lJ80UidfmpW2OaYMHqwsqpO1vX/fDxscwAAAACAX4S10JpZnKQHJE2VNFLSTDM7PPQY59z/OOeO\ncc6NkXSzpFznXEk4c7XGEx8u0sjkqW16dePG3HDytXp20/2qr3fhnQgAAAAAoly4V2jHSlrnnNvo\nnKuRNFfSgU5CnSnpqTBnapWlOxdr+shJYZ/nhvMnyFmd/jr//bDPBQAAAADRLNyFto+kzSHP84Kv\nfYuZpUg6XdJzYc7UYtvKC1Uat16zpo0N+1xxcabJmT/U/+b+PexzAQAAAEA0S/A6QIjpkt490Hbj\n2bNn73ucnZ2t7Ozs8KeS9Pc331DHnePVv09iu8z3+5nf1VEPDdfm7feqX1andpkTAAAAANpTbm6u\ncnNzD2kMcy5852qa2YmSZjvnTg8+/5Uk55z7QyPHPi/paefc3CbGcuHMeiBjf/sjpVSM0Ft/uK7d\n5uz7/12g7H5T9fj1P2q3OQEAAADAK2Ym55y15DPh3nK8VNJQMxtgZkmSZkia1/AgM+ss6TRJL4U5\nT6us3JWrC8ZMbNc5fzx2ll7Y9H/tOicAAAAARJOwFlrnXJ2kqyUtlLRS0lzn3Cozu8rMQpcez5W0\nwDm3J5x5WqOgdId223bNnHRku8574wVTVJVQoGfe/rxd5wUAAACAaBHWLcdtyastx/e+Ok+/fulB\nVfx1QbvPfcrsW1RZXatld93d7nMDAAAAQHuKxC3HUW/+5+9rWMpJnsx985mX6dOaJ1VTW+/J/AAA\nAAAQySi0B7G86H1lD/Wm0J55/Egl12Xq3hff8mR+AAAAAIhkFNoDqK6rVmHCJ7p0/AmeZcjpcZke\n/vAJz+YHAAAAgEhFoT2AxV8sV1zJEI0Z6d29YGdfMFNfJT6vorJKzzIAAAAAQCSi0B7A80s+VG93\nouI8/C2NGdpXXapG6XdPz/cuBAAAAABEIArtASzJ+1ijuh/ndQydO+QyPfkF244BAAAAIBSF9gA2\nVH6iCYeP8TqGbr/4Qm1NXaRvtpZ5HQUAAAAAIgaFtgm7qnerIukrnTPuSK+jqF/3LupRdap+/9y/\nvY4CAAAAABGDQtuEN1Z+roTiwzW4f7LXUSRJ5xx2oV5a+6zXMQAAAAAgYlBomzD/00/Uw42RmddJ\nAm4+/xxtTVusvO0VXkcBAAAAgIhAoW3Cks2faGSG9+fP7jUwq6u6V47T75/jascAAAAAIFFom/T1\n7k80fljkFFpJOnvohXphDduOAQAAAECi0Daquq5aZUmrdO6Jo7yOsp9bzj9X+akLtHXnbq+jAAAA\nAIDnKLSN+OirNbKy/hpxWKrXUfYzuGemulUerz8895rXUQAAAADAcxTaRiz+bKW6VB8ZMReECnXm\noIv07JfPeB0DAAAAADxHoW3E0m9WakDqSK9jNOqW889VXsqrKiqt8joKAAAAAHiKQtuI1cVf6Oie\nkVloh/fpoU5VI3T/vLe8jgIAAAAAnqLQNmJr3UqdOvxIr2M06bSe5+hfy1/yOgYAAAAAeIpC28Du\n6j3ak7RZU449zOsoTfr5lLO11uapttZ5HQUAAAAAPEOhbeC9tasVXzJUfXsneh2lSZNGHa4EpWjO\nwk+9jgIAAAAAnqHQNvDG5yvVtW5kRF7heC8z05i0s/X39+Z5HQUAAAAAPEOhbWDpxi80MC0yLwgV\n6ocnn6NPKjiPFgAAAEDsotA2sK5kpUb3jtwLQu31vUnjVJOSpzc+3uR1FAAAAADwBIW2ge1ulU4e\ndoTXMQ4qMT5Bh+lM3fca244BAAAAxCYKbYjqumpVJuXptKMHex2lWS4adbbe2kqhBQAAABCbKLQh\nVmzaICvvqwF9k7yO0izXnTVFpekfat2mUq+jAAAAAEC7o9CGeGflWqVXD4voKxyH6pbeUT2rT9H/\nm7fA6ygAAAAA0O4otCGWfbNWPRKHeR2jRXIGnKlX173qdQwAAAAAaHcU2hBrdqzV0C7RVWivnnqG\nNiW/qqrqeq+jAAAAAEC7otCGyNuzVqP6RlehPWHYYCXXZ+ifCz71OgoAAAAAtCsKbYjiuLU6aXh0\nFVpJGt1xmh77cL7XMQAAAACgXVFog8qrKlSTUKyTj+rrdZQWu2zsNH1STqEFAAAAEFsotEEfrFmn\n+LKh6poRfb+SH0w6RZUdv9Snawq9jgIAAAAA7Sb62luYvLt6rbrURt92Y0lKTU5W35qJun8+t+8B\nAAAAEDsotEErNq9Vnw7RWWglaeqQaVq4gW3HAAAAAGIHhTbo6+KvNLTrUK9jtNrPp52hgtQF2rW7\nzusoAAAAANAuKLRB26o3aGTvQV7HaLWj+vdVSm0fPTT/I6+jAAAAAEC7oNAGlcat15jB0VtoJen4\nzmfqqWVsOwYAAAAQGyi0kqpqq1STuEMnjoi+W/aE+t5J07Ri93w553USAAAAAAg/Cq2kzzdvlFX0\nVY/uCV5HOSTfyT5RNWnf6J0V+V5HAQAAAICwo9BK+mjtBqVWD5KZ10kOTWJ8ggbU5eihRYu8jgIA\nAAAAYUehlbRi43plxg32OkabmDJkit7YzP1oAQAAAPgfhVbS2sIN6pMW3ReE2uvq06dqa8oi7ams\n9zoKAAAAAIQVhVZSXvkGHZbpj0J71IB+Sq7vrjkLPvU6CgAAAACEFYVW0o7aDTqqnz8KrSQdnTpF\nTy5h2zEAAAAAf6PQSqpIXK+xh/njHFpJumjMVH1SutDrGAAAAAAQVjFfaIt3l6reqjVmeKbXUdrM\nlZPHa1fnj/XVpnKvowAAAABA2MR8oV361QYllA9SWlqU37MnROfUNHWvOkF/np/rdRQAAAAACJuY\nL7TLvtqg9Dr/nD+71ym9pujfaziPFgAAAIB/xXyhXZn/jbonDvQ6RpublT1VX2uh6rl7DwAAAACf\nivlCu7F4s3qn9fc6Rps7fcxRUnKZ5r+/wesoAAAAABAWMV9oC3Zv1qBu/byO0ebiLE5DbYr+/hbb\njgEAAAD4U8wX2qLaTRre03+FVpLOPHyK3i3g9j0AAAAA/CnmC21F/GaNGui/LceS9LPTJ6sw/Q0V\nl9Z4HQUAAAAA2lzYC62ZnW5mq81srZnd1MQx2Wb2qZl9YWZvhjvTXtV11apNKtSYw3q115TtanCP\nHupYM0gPv7rE6ygAAAAA0ObCWmjNLE7SA5KmShopaaaZHd7gmM6SHpR0lnPuSEkXhTNTqLUF+bJd\nPdU9M769pmx3x3aZqmc/5TxaAAAAAP4T7hXasZLWOec2OudqJM2VdE6DYy6V9JxzboskOecKw5xp\nn0++3qQOVf1l1l4ztr+ZY6fo892cRwsAAADAf8JdaPtI2hzyPC/4Wqhhkrqa2ZtmttTMLg9zpn1W\n5m1WJ/nzglB7fTf7ZFV1+lLL1xR5HQUAAAAA2lSC1wEUyDBG0kRJaZI+MLMPnHNfNTxw9uzZ+x5n\nZ2crOzv7kCZet22zuif5u9CmJCWrZ80p+tvCN/WX4Rd4HQcAAAAAJEm5ubnKzc09pDHCXWi3SAq9\nhHDf4Guh8iQVOucqJVWa2duSRkk6YKFtC5tKN6lfpxFtOmYkOq1vjhZ9tVgShRYAAABAZGi4SHn7\n7be3eIxwbzleKmmomQ0wsyRJMyTNa3DMS5JOMbN4M0uVdIKkVWHOJUnaVrlZgzP9vUIrST84bbI2\nxC1Wfb3XSQAAAACg7YS10Drn6iRdLWmhpJWS5jrnVpnZVWb2o+AxqyUtkPSZpA8lPeSc+zKcufYq\nqd+sI3r7v9DmHH2kLLlcr37wjddRAAAAAKDNmHPO6wzNYmaurbPG39JV716yVuNGZbbpuJFo2M2X\n6ehOE/TszbO8jgIAAAAA32Jmcs616B404d5yHLEqqnapPq5SRw3p5nWUdjH1sBy9l7/Y6xgAAAAA\n0GZittB9zxyxAAAgAElEQVR+tnGz4ir6qmNHH9+ENsRVkydpa+rr2r2HE2kBAAAA+EPMFtrlGzYr\ntcb/58/udWS//urguurxRZ95HQUAAAAA2kTMFtp1W/PVyfp4HaNdjUzJ0b+Wsu0YAAAAgD/EbKH9\nZme+MpN7ex2jXZ17VI6WFVFoAQAAAPhDzBbaLWX56tUxtgrtlZMnqKzz+9q6o8rrKAAAAABwyGK2\n0O6ozFf/jNgqtD06d1HnmhH626vvex0FAAAAAA5ZzBbakrp8DcmKrUIrScd3zdFLn7PtGAAAAED0\ni9lCuysuX4f37eV1jHY3Y2yOvqyk0AIAAACIfjFZaJ1zqkneqqMGxl6hvfTUcaru9KVWrCn2OgoA\nAAAAHJKYLLRby3ZKVR3Vv3cHr6O0u5SkZPWoOVkPLcr1OgoAAAAAHJKYLLRfbMxXwp7eSkjwOok3\nTu2do4Vfse0YAAAAQHSLyUL75eZ8pdbF3gWh9vpB9mRt0GLV13udBAAAAABaLyYL7bqt+eocF7uF\ndsqoo6SUYi34cJPXUQAAAACg1WKy0G4syldmh9gttHEWp0GapH+8xbZjAAAAANErJgttfnm+eqfH\nbqGVpClDc/ROPoUWAAAAQPSKyUJbWJWvAV1ju9BeNTlH21Je155KTqQFAAAAEJ1istCW1OdraI/Y\nuwdtqKP7D1CSOunJxV94HQUAAAAAWiUmC+3uuHwd0Te2V2glaWSHHD31EduOAQAAAESnmCu0dfV1\nqk3epqMG9fQ6iufOPipHy4ootAAAAACiU8wV2rziHVJlhnplJXkdxXNXTZmo0s7valthtddRAAAA\nAKDFYq7QfvFNgRIreyku5r75t/XsnKFONYfrofkfeB0FAAAAAFos5mrd2vytSqlju/Fex2XkaN4X\nr3sdAwAAAABaLOYK7YYd25Qe18PrGBHjkrGT9MVuzqMFAAAAEH1irtDmlWxT1yQK7V7fGX+yqjp/\nri++KvU6CgAAAAC0SMwV2q3l29Q9lUK7V2pSB2VVn6iHFrzldRQAAAAAaJGYK7SFldvUuxOFNtTJ\nvXO0YB3bjgEAAABEl5grtKW129SvG4U21PfH5+hrt1jOeZ0EAAAAAJov5gpthbZpcBaFNtS0Y0bL\npW3T60u3eB0FAAAAAJot5gptVcI2DetNoQ0VHxevgfUT9Y833/A6CgAAAAA0W0wV2tr6WtUlFuvw\n/pleR4k4OYNz9FYe59ECAAAAiB4xVWjzigqlygx175bgdZSIc2VOjvKTF6uqihNpAQAAAESHmCq0\nq/O2KaGqh8y8ThJ5jh00WInxSXr6zdVeRwEAAACAZompQrsuf5tS6jh/tjFmpsOTJunJD9h2DAAA\nACA6xFSh3bBjmzoahbYpZ43I0UeFFFoAAAAA0SGmCm1e8TZlJFJom/KjyRNV3OktFZXUeh0FAAAA\nAA4qpgptQdk2ZaVSaJsyIDNLHWsH6u+vLvU6CgAAAAAcVEwV2h2V29SrE4X2QI7pnKPnl7PtGAAA\nAEDki6lCW1KzTf26UmgP5MIxOfqsgkILAAAAIPLFVKGtcNs0KItCeyDfn3Cqdnf5RF9t3OV1FAAA\nAAA4oJgqtJUJ2zSsN4X2QDqlpCmz+lj97bV3vI4CAAAAAAcUM4W23tWrLqlQR/TL8jpKxBvXI0fz\nV7PtGAAAAEBki5lCW1C6U6rqpB7dE72OEvEuP3mS1tYulnNeJwEAAACApsVMoV2Tt13xVVmKi5lv\n3Hrnjj1edenf6L3l272OAgAAAABNipl6t37rDnWoZbtxcyTGJ6hf/Wn6v9ff8DoKAAAAADQpZgrt\npsJCpSrT6xhRY8KAHL25kfNoAQAAAESumCm0ecU7lJ5AoW2uWRNztDlpkWpqOJEWAAAAQGSKmUK7\nrbxQGUndvY4RNU4edrjiE2v14ttfex0FAAAAABoVM4V2x64dykxlhba5zEzDEnL02HtsOwYAAAAQ\nmWKm0BZXFapnOiu0LTFteI4+2Pa61zEAAAAAoFExU2jL6naodxdWaFviqimTVJj+hkrL6ryOAgAA\nAADfEjOFdrcrVP/urNC2xNAevZVa31NzFi73OgoAAAAAfEvMFNrK+EINymKFtqWO7jhJz3zMebQA\nAAAAIk9MFFrnnGqTdmhoHwptS50/OkfLSym0AAAAACJP2AutmZ1uZqvNbK2Z3dTI+6eZWYmZfRL8\n+XVbZyiv2iXVx6tfj9S2Htr3fjDxNFVkfKhvtuzxOgoAAAAA7CeshdbM4iQ9IGmqpJGSZprZ4Y0c\n+rZzbkzw5862zrF+6w5ZZaaSktp6ZP/r1rGzMqqP0sOvve91FAAAAADYT7hXaMdKWuec2+icq5E0\nV9I5jRxn4QzxVUGhkmq4IFRrndA9R6+sZNsxAAAAgMgS7kLbR9LmkOd5wdcaGmdmy83s32Y2oq1D\nfLN9hzrUc/5sa116Yo5WVS+Wc14nAQAAAID/SPA6gKSPJfV3zu02szMkvShpWGMHzp49e9/j7Oxs\nZWdnN2uCvKJCdYxjhba1Lj7pRH1vwRotW1mk44/s6nUcAAAAAD6Qm5ur3NzcQxrDXBiX3czsREmz\nnXOnB5//SpJzzv3hAJ/ZIOlY51xRg9dda7NefO89Wrk5Tyvv+d9WfR5S35umaVrPWXro+vO9jgIA\nAADAh8xMzrkWnY4a7i3HSyUNNbMBZpYkaYakeaEHmFmPkMdjFSjZRWpD2ysK1a0DW44PxWn9crR4\nPefRAgAAAIgcYS20zrk6SVdLWihppaS5zrlVZnaVmf0oeNiFZvaFmX0q6V5Jl7R1jp17dqh7R7Yc\nH4orsidpY/xi1dV5nQQAAAAAAsJ+Dq1z7jVJwxu89reQxw9KejCcGUprCtWrMyu0h2LiyKNkKSX6\n93sbdfb4AV7HAQAAAICwbzmOCBX1O9QvgxXaQxFncRoSN0lz3n7d6ygAAAAAIClGCu0eK9SALFZo\nD9XUw3L0XgHn0QIAAACIDDFRaKsTd2hIT1ZoD9VVk3O0Le11Veyq9zoKAAAAAPi/0NbW16o+sUxD\n+2R4HSXqjewzQB2sk55Y9IXXUQAAAADA/4U2v7hIquyiLp3jvY7iCyNTcvSvpWw7BgAAAOA93xfa\nrwp2KL4qU9ai2/OiKecenaOPiym0AAAAALzn+0K7aXuRkuu4IFRbuTJngsq6vKf8bdVeRwEAAAAQ\n43xfaPOKdqqD6+p1DN/o0amrOtcO0/+9+pHXUQAAAADEON8X2oKSInWMp9C2pWO7TtJLX7DtGAAA\nAIC3fF9ot5cXKT2RQtuWZhyfo5V7KLQAAAAAvOX7Qlu4q0gZyd28juErl516sqq7fKYVq8u8jgIA\nAAAghvm+0BZXFikzlRXatpSalKIetSfo4UVveR0FAAAAQAzzfaEtrdmp7ukU2rZ2Su8cLVzHtmMA\nAAAA3vF9od1VV6SenSm0be1743O0XotVX+91EgAAAACxyveFdo+K1Lcr59C2tTNGHSPXcasWf5Tv\ndRQAAAAAMcr3hbYqrkj9MlmhbWvxcfEapAl6JPd1r6MAAAAAiFG+L7S1iUUa2INCGw45g3P0Vt4i\nr2MAAAAAiFG+LrTVddVy8ZXq3yPd6yi+9NOpU7UtbaEqdnEiLQAAAID25+tCu6WoSKrsqrQ08zqK\nLx3db5CSLV1zXvvM6ygAAAAAYpCvC+2GrUVKqGa7cTiNSpuquUsXeh0DAAAAQAzydaHdtKNISXUU\n2nC6eMxUfVK2wOsYAAAAAGKQrwvtlqIidRCFNpx+OGmCdmcs0Zr1u7yOAgAAACDG+LrQ5pfsVMc4\n7kEbTp1TOiqr5jg9OD/X6ygAAAAAYoyvC+328iKlJ7JCG26n9pmiV9ey7RgAAABA+/J1od25q0gZ\nyRTacLtywlStt4Wqq/M6CQAAAIBY4utCW1RZpG6pFNpwm3zUaMWlFmve2xu9jgIAAAAghvi60JbV\n7FRWOufQhlucxemw+Ml65G22HQMAAABoP74utBV1RerVmRXa9jB9xFS9v41CCwAAAKD9+LrQ7lGR\nenel0LaHn0yZrKLOb2jHzlqvowAAAACIEb4utFVxReqXSaFtDwMzeyq9boAe+vcSr6MAAAAAiBG+\nLrS1iUUa3JNzaNvL8RlT9dwKth0DAAAAaB++LbTVddVy8XvUv0e611FixnfGTdXKygVyzuskAAAA\nAGKBbwvtlqJiqTJDqanmdZSYMfPkk1XT5Ut9uKLI6ygAAAAAYoBvC+0324qUUMP5s+2pQ2Ky+taN\n118XLvY6CgAAAIAY4NtCu2lHkRLrKLTtLWfQFL2xcaHXMQAAAADEAN8W2oLiYnVwGV7HiDk/njxV\nW1IWaM8eTqQFAAAAEF6+LbTbSkuUal28jhFzjh80TEkJ8Xpy0SqvowAAAADwOd8W2sKKEnVMYIW2\nvZmZjkyZqsc/5PY9AAAAAMLLt4W2aHeJOiWxQuuFC0dP1bJiCi0AAACA8PJvod1TrM7JFFovXJkz\nURVd39fXm/Z4HQUAAACAj/m20JZVl6hbGluOvdAtrYu6147WAy+/5XUUAAAAAD7m20JbUVuizDRW\naL1yWp8z9PLqV72OAQAAAMDHfFtod9cXK6sThdYrP5k0TRvi56umxuskAAAAAPzKt4W2UiXqlcGW\nY69MGHG04lN26+nF67yOAgAAAMCnfFtoq+NK1LsrK7ReMTMdmXyGHnl3vtdRAAAAAPiUbwttTUKx\n+mZSaL108THTtKSI82gBAAAAhIcvC229q5dLLFf/rM5eR4lpV03O0a6u72n117u8jgIAAADAh3xZ\naIt3lUvVaercKd7rKDEtI7WTsmqP059eedPrKAAAAAB8yJeFdtOOYsVVd1GcL79ddJnYf5rmr2Pb\nMQAAAIC258vKl1dYooRazp+NBD+bPE0bk+arstJ5HQUAAACAz/iy0OYXlSipnlv2RIKTDhuhxKR6\nPbFwtddRAAAAAPhMswqtmXU2s9PN7MdmdlXwccRecWlrSYk6iBXaSGBmGpV6hh59n23HAAAAANrW\nAQutmZ1iZvMkvS1ppqQBkgYFH79jZi+Z2Snhj9ky28uKlRZPoY0UM4+bpmWl3I8WAAAAQNs62Art\n+ZJucM6Ncs59zzl3s3PuV8HHR0v6ZfCYJgVXc1eb2Vozu+kAxx1vZjVmdsDxmmNHRYk6JrDlOFL8\ncOJE7en6kT5bU+51FAAAAAA+csBC65z7/yR9bWYXN/H+2uAxjTKzOEkPSJoqaaSkmWZ2eBPH/V7S\nghZkb1Lx7hJ1TmKFNlJ06tBRPetO1IP/fsPrKAAAAAB85KDn0Drn6iXd2Mrxx0pa55zb6JyrkTRX\n0jmNHHeNpGclbW/lPPsprixWlxQKbSSZPHCaXv2abccAAAAA2k5zr3K82Mx+YWb9zKzr3p9mfK6P\npM0hz/OCr+1jZr0lneuc+4ska2aeAyqvKVFmGluOI8nVU6cpL2W+du/m9j0AAAAA2kZzC+0lkn6m\nwMWhPg7+LGujDPdKCj239pBL7a66EmV2ZIU2khw3cJiSExL1z1e/8DoKAAAAAJ9IaM5BzrlBrRx/\ni6T+Ic/7Bl8LdZykuWZmkjIlnWFmNc65eQ0Hmz179r7H2dnZys7ObnTS3fXF6tGZQhtJzEzHdJym\nJz56VT+94Civ4wAAAADwWG5urnJzcw9pDHOu6S2gZnaac+6tAw5glu2cazSFmcVLWiNpkqQCSUsk\nzXTOrWri+H9Ietk593wj77kDZQ2V+ouj9NCUJ/WdKRSnSPLnRfN1w/N/0O4/vyVrk83lAAAAAPzC\nzOSca1FTONiW47PM7CMzu8vMzjezcWZ2UvDxf5vZUklnNPVh51ydpKslLZS0UtJc59wqM7vKzH7U\n2EdaEr4pNfEl6t2NFdpI8/3TJqiq63J9sKLI6ygAAAAAfOCAK7SSZGYdFbgy8Sn6z/bhjZLelTTP\nOVcR1oT/ydHsFdq4X3fUmisLdNiA9DCnQksNuuUcndDxYs295TKvowAAAACIIOFYoVWwsA6VtFWB\nLcNLgo+HtleZbYmauhq5uCr1zerodRQ04pwjpuv1zS97HQMAAACADzT3KscVIT+1CmwzHhimTIdk\na0mpVNVZKSmcpBmJrpt2pgq7LFD+1hqvowAAAACIcs0qtM65e0J+ficpW9LgsCZrpc07ShRfw/mz\nkWpgt17KqD9M9734jtdRAAAAAES55q7QNpSqwC14Ik7ezmIl1lFoI1l2r+l6fiXbjgEAAAAcmmYV\nWjP73Mw+C/6sVOBWPPeGN1rr5BeVKLk+w+sYOICrp0zX+oSXVVnZJhe1BgAAABCjEpp53Fkhj2sl\nbXPO1YYhzyHbXlaiDtbZ6xg4gAlHjFJCcrUef221Zp17hNdxAAAAAESp5p5DuzHkZ0uklllJKiwv\nVWochTaSmZmOSTtL/3ifbccAAAAAWq+159BGrOLdZUpLpNBGuu+Nm66Py19WM28tDAAAAADf4rtC\nW7KnVOmJnbyOgYP4/mkTVN31M72zbKfXUQAAAABEKd8V2tKqUnXuwAptpEtJ7KBBmqj7X5vvdRQA\nAAAAUcp3hbaipkwZKRTaaHDeiOl6cwvn0QIAAABoHd8V2l11pcpIZctxNLh22pkqzlikzfnVXkcB\nAAAAEIV8V2j31JUpM50V2mjQL6OHMuqH674X3/Y6CgAAAIAo5LtCW6lS9ehMoY0WE/tM1wtfsu0Y\nAAAAQMv5rtBWx5UqqzNbjqPFNVOna0Piy9qzh/v3AAAAAGgZ3xXa2rgy9erKCm20OHXYUUpKdnrk\n3yu9jgIAAAAgyviu0NYllqp3N1Zoo4WZ6YRO5+qR91/wOgoAAACAKOOrQltTVyPF1ahn11Svo6AF\nfjLhXH1W/YLq6rxOAgAAACCa+KrQbi8rk6o6KSXFvI6CFrjohFOkTnl67vVvvI4CAAAAIIr4qtBu\n2VmquFq2G0eb+Lh4HZU8XX9+4yWvowAAAACIIr4qtAVFZUqo5YJQ0ej7J56rD0tfkONixwAAAACa\nyVeFdmtxqZIcK7TRaNaEHFV3/VRvLtnhdRQAAAAAUcJXhXZHWamSHSu00Sg1KUVDbbLue/UVr6MA\nAAAAiBK+KrSF5WVKiaPQRqsZo8/Tm1u5fQ8AAACA5vFVod25q1Sp8Ww5jlbXnnGmKrrl6vM1FV5H\nAQAAABAFfFVoS3aXqWMiK7TRqltaF/WuP1H/8+ICr6MAAAAAiAL+KrSVpeqUzAptNDt3+Hl6dQPb\njgEAAAAcnK8KbXl1qbp0YIU2mt1w1jkqzJivzfnVXkcBAAAAEOF8VWgrasuUkUqhjWaDMnsrww3T\nPc/leh0FAAAAQITzVaHdXVeqbh3ZchztpvY/Ty+setHrGAAAAAAinK8KbaUrU/d0Vmij3Q3TztXm\ntBdVUlrvdRQAAAAAEcxXhbbKSpXVmUIb7Y4dOFxp8Rm677kPvY4CAAAAIIL5qtDWxJWqZwZbjv1g\nUs+L9Ngnz3odAwAAAEAE81WhrY0vU++urND6wS/PulDrk59VWTnbjgEAAAA0zjeF1jknl1Sm3t3S\nvY6CNnDSkJFKTUzTfc8u8ToKAAAAgAjlm0JbXrlbqk1W1y6JXkdBGzAzTexxkR77mG3HAAAAABrn\nm0K7ZWeprLqz4nzzjXDjmRfpq6RnVVHhvI4CAAAAIAL5pv7lF5UqvpYLQvnJyYcdqZSkZN3/3FKv\nowAAAACIQL4ptFuLy5RYxwWh/MTMlJ11keYsfcbrKAAAAAAikG8K7baSUiU5Vmj95pfTLtI6th0D\nAAAAaIRvCm1heZk6GCu0fnPa8KOVkpSoB57/2OsoAAAAACKMfwptRalS4ym0fmNmGp95of7JtmMA\nAAAADfim0BbvLlNaPFuO/eiX0y7S2oRntWsX244BAAAA/IdvCm1pZZnSktK9joEwmHDEaHVINv35\nhU+9jgIAAAAggvim0JZXlyudQutLgW3HF+kfH7HtGAAAAMB/+KbQ7qouV+cObDn2qxvOuFCrE55h\n2zEAAACAffxTaGvL1SWFFVq/yhkxRh2STX96jqsdAwAAAAjwTaHdU1+uLqkUWr8yM+X0mKm/f/SU\n11EAAAAARAjfFNpKV6ZuHSm0fnbrOTP1dcq/VFxS73UUAAAAABHAN4W2WuXK7ESh9bMTBh+h9PhM\n/WHuO15HAQAAABABfFNoa+LK1Z1C63tnDpipx1ew7RgAAACAjwptXXy5srpQaP3uv86bofzOz2pz\nfrXXUQAAAAB4zD+FNqFcPTMotH53eK8B6qbh+t3cRV5HAQAAAOAxXxTauvo6Kb5SPTLSvI6CdnDB\n8Jl6bg3bjgEAAIBYF/ZCa2anm9lqM1trZjc18v7ZZrbCzD41s2VmNrGlcxTvqpCqOyo11domNCLa\nredepJ3dXtGqr3Z7HQUAAACAh8JaaM0sTtIDkqZKGilpppkd3uCwxc65Uc65YyRdIemhls5TUFwu\nq0mX0WdjQr+MHuqjsfrtv17xOgoAAAAAD4V7hXaspHXOuY3OuRpJcyWdE3qAcy50ma2jpMKWTrK9\nuFzxdZw/G0suO3qm/r2JbccAAABALAt3oe0jaXPI87zga/sxs3PNbJWk+ZJ+3tJJtpWWKYFCG1Nu\nPPs8lWe+oQ+Xl3gdBQAAAIBHErwOIEnOuRclvWhmp0h6TNLwxo6bPXv2vsfZ2dnKzs6WJO0oLVeS\no9DGkq6pXTTYJup3z7+gl0df4XUcAAAAAC2Um5ur3NzcQxoj3IV2i6T+Ic/7Bl9rlHPuXTNLMLNu\nzrmdDd8PLbShdlaUK9kotLHmyhMu1ex//1XOXcH50wAAAECUCV2klKTbb7+9xWOEe8vxUklDzWyA\nmSVJmiFpXugBZjYk5PEYSWqszB5I8a5ydaDQxpyfnz5dNV2X64U3Nh/8YAAAAAC+E9ZC65yrk3S1\npIWSVkqa65xbZWZXmdmPgoddYGZfmNknku6TdElL5yneXa7UBAptrElJ7KBjUy7UH+Y/4XUUAAAA\nAB4I+zm0zrnX1OCcWOfc30Ie3y3p7kOZo7SyXGkJnQ5lCESpm874ri5+9Efas+cmpaSw7xgAAACI\nJeHectwuyivL1TGJFdpYdN6xJykptVL3P/OJ11EAAAAAtDN/FNrqcnVKptDGIjPT1F6X668fPup1\nFAAAAADtzBeFdldtuTp3oNDGqtsvuFwb05/SloIar6MAAAAAaEe+KLS768rUJZVCG6uO7jtEmTZM\ns594zesoAAAAANqRLwptZX25unak0MayS4+8XM+ufczrGAAAAADakS8KbZXK1Y1CG9N+fe7FKu2+\nQO99Uux1FAAAAADtxBeFtsbKldmJQhvLMjtmaFj8FN3x7DNeRwEAAADQTvxRaOPKldWZQhvrrjn1\nu8otflR1dV4nAQAAANAefFFo6+LL1TOjk9cx4LEfTTxd9Rlr9egrX3kdBQAAAEA7iPpC65yTSyxX\nzwxWaGNdYnyiTul0mf5n8T+8jgIAAACgHUR9oa2srZJcvLp2SfQ6CiLAnRf8UKuS52hHIfuOAQAA\nAL+L+kK7tbhcqk5XfLzXSRAJTj7sSGUk9NFv5izwOgoA4P9v787Doyrv94/fzyxZJxtJ2DUo4vKl\nVcRqXb8ibmhVrPoVsVoVN0C0LnXBFWoVKa1WsWBVtO5brdWf1q1KFOuOO4osFoSyBkjILElmeX5/\nJGqkIAQyPHNm3q/rOtfMnJwc78vLYbj9POcMAABplgWFdq38CZYb4zu//PGZeuTLabLWdRIAAAAA\n6eT5Qrt8TaP8SQotvnPdcScpXP2KXn5zhesoAAAAANLI84V25dpGBVMUWnynvLBUu+Ydq3F/e8B1\nFAAAAABp5PlCW9fYqDxRaPF9Vx15pt6JT1MkwrpjAAAAIFt5vtCuDjeqwFBo8X3H/2R/FRQlNPHh\nt11HAQAAAJAmni+0ayKNKvSVuo6BDGOM0bE1I3TX+9NcRwEAAACQJp4vtA2xRhUFmNDiv9047DQt\n7/KkPv4i7DoKAAAAgDTwfqFtalQoSKHFf6vp0kPb+Q/QlQ897joKAAAAgDTwfKFtbGlUST6FFut3\n0aCz9PLquxSPu04CAAAAoLN5vtBG4hRabNjIwUfKlC3WbY997DoKAAAAgE7m+UIbTYRVVhByHQMZ\nKuAL6OieZ+uWN+5wHQUAAABAJ/N8oW1KhVVeRKHFhv1u+JlaUvGoPvy80XUUAAAAAJ3I84W22VJo\n8cO2r+qlfoHBuvT+h1xHAQAAANCJPF9oWxRWlxCFFj/sqsNHqjYyVdGodR0FAAAAQCfxfKGNK0Kh\nxUadst/Byi+O6qYH33YdBQAAAEAn8XyhTfjCqiwtdh0DGc5nfDqx77ma+v5U11EAAAAAdBLPF9qk\nP6zqMia02LgJw07XqqpnVPvuKtdRAAAAAHQCzxfaVDCsruUUWmxc99Iq7Zp/tMY+9hfXUQAAAAB0\nAk8X2kQqIflaVFVW4DoKPGL8MSP1bvLPqm9IuY4CAAAAYAt5utCuCUeklpDy843rKPCIYwbsq5KC\nQl017WXXUQAAAABsIU8X2pUNEZlEsQx9FpvIGKNzd79A982+TSmGtAAAAICnebzQhuVPcv0sOua6\nn5+s5qr3dM/Tc11HAQAAALAFPF1o69ZSaNFxRXmFOrz6LF3/4mTXUQAAAABsAW8X2sawgpZCi467\n5RejtKjiQb3/6VrXUQAAAABsJk8X2tXhsPJEoUXH9eu6jXbJO1QX3fcX11EAAAAAbCZPF9qGaER5\npth1DHjUjUMv0JvJyVq9hrtDAQAAAF7k6UJbHw2rwMeEFpvnmN33VXlBmS6763nXUQAAAABsBk8X\n2oZYWEV+Ci02jzFGY/a8QA/Nu1XJpOs0AAAAADrK04W2sSmsogCFFpvvymOGKVn1iab+9XPXUQAA\nAAB0kLcLbUtYRUGuocXmyw/k6+geI3XjK7e6jgIAAACggzxdaCPxiErymdBiy9x26mgtq3pcL7yx\n3HaHofsAAB8DSURBVHUUAAAAAB3g6UIbS4QptNhivcq7au/QMF386O2uowAAAADoAE8X2mgyrLJC\nCi223J9OuUSzQ3fosy8jrqMAAAAA2ESeLrRNybDKiyi02HK7b9tP/fIO0Oi77nEdBQAAAMAm8nSh\nbVZY5cXcFAqd4/fHXap/pW7WshUJ11EAAAAAbAJPF9q4jahLiAktOsfRA/ZRdUEvXTD1SddRAAAA\nAGwCbxdaX1iVFFp0oqsHX6qnVkxSNGpdRwEAAACwEZ4utElfWFWlFFp0ntEHH6380rCuurvWdRQA\nAAAAG+HtQhsIq7qcQovO4zM+jd7tEt352SQlk67TAAAAAPghni201lrZQFhdy7kpFDrX+ONPVaLy\nY0168EPXUQAAAAD8AM8W2uZEiySjitI811GQZQqDBTptx19rwhs3KJVynQYAAADAhni20NatjUgt\nIfn9rpMgG9188jmKVs/Q7Y/Pch0FAAAAwAakvdAaY4YYY2YbY+YYYy5fz89PNsZ83La9YYz58aac\nd0V9WL4E188iPUL5xTqp5kKNf2WCLDc8BgAAADJSWgutMcYn6XZJh0vqL2m4MWbndQ77StL/Wmt3\nk/RbSXdtyrnr1oblT1JokT6Tf3meGqpf0L1Pz3MdBQAAAMB6pHtCu5ekudbahdbauKRHJQ1tf4C1\n9m1rbUPby7cl9dqUE9etDSuQ4oZQSJ/ywlIN7Xmexj53E1NaAAAAIAOlu9D2krSo3evF+uHCepak\n5zflxKsawwpaJrRIr6mnXaC66qf02Atfu44CAAAAYB0B1wG+YYw5SNIZkvbf0DHjxo379vnCVJ7y\nDIUW6dW1pFKHV5+lS5+apJOOmOw6DgAAAJA1amtrVVtbu0XnMDaNaymNMXtLGmetHdL2+gpJ1lo7\ncZ3jdpX0pKQh1tr5GziXbZ/14nse0eOfPKPFf3wkbfkBSVq8ZrlqJu2ivx3ymYYO7uk6DgAAAJCV\njDGy1pqO/E66lxy/J2kHY0yNMSZP0kmSnml/gDFmW7WW2VM3VGbXpz4WVoGPa2iRfr0ruungLmdo\nzGMTXEcBAAAA0E5aC621NilpjKSXJM2S9Ki19gtjzLnGmHPaDrtGUhdJU4wxHxpj3t2Uczc2hVUU\nYMkxto57z7pcSyof5lpaAAAAIIOkdclxZ1p3yfHhN16vtZFmvXXDbx2mQi45+taxev+z1Vpy559l\nOrQQAgAAAMDGZOKS47SJtERUHGRCi63nnjMv1crqJ3XfM1+5jgIAAABAXi608bBK8im02HqqQ130\n815j9Ov/9xu+lxYAAADIAJ4ttNFEWCUF3BQKW9edZ1yk+q7PaeoTs11HAQAAAHKeZwttLBlWWSET\nWmxdFUVlGt7nYl350nilUq7TAAAAALnNs4W22YZVXkShxdY35bTzFamerpsf/NR1FAAAACCnebbQ\nttiIKooptNj6SvJDOnPnKzRuxpWKx12nAQAAAHKXdwutCauimGto4catvxilVNUsXTal1nUUAAAA\nIGd5ttAmTERdQhRauJEfyNc1+96gP825TGvXcstjAAAAwAXPFtqkP6KKUJHrGMhhlx81TCWlKY34\nwxOuowAAAAA5ybOFNuWPqrKUCS3c8Rmfbh86SU81jtXCxS2u4wAAAAA5x7OF1gYiqi6j0MKt4Xsf\npJrinXTyzXe4jgIAAADkHE8W2pZEXDJJlZfkuY4C6P5TJ+qt4A169+MG11EAAACAnOLJQtsQjUrx\nYuXlGddRAO2/44/1k9Kf6ZQ7JrqOAgAAAOQUTxbalQ0RmQTLjZE5Hjnnes0vv1MPPPtv11EAAACA\nnOHJQruqMSpfkjscI3P0re6lk2ou0nlP/1qJhOs0AAAAQG7wZKFd0xiRP8mEFpnlrhEXK171gS6a\n/KrrKAAAAEBO8GShXR2OKmCZ0CKzFOUV6sZBv9fUf/9Ky1cypgUAAADSzZOFtj4SUdAyoUXmufCw\n49S9pFonTrzTdRQAAAAg63mz0EYjChoKLTKPMUaPnPZHzfCP04z3V7uOAwAAAGQ1TxbatdGo8g1L\njpGZDthxV+3f5QQN//N1stZ1GgAAACB7ebPQNkWU72NCi8z1+KjfaHnV45p4/4euowAAAABZy5OF\nNtwcVYGfCS0yV/fSKl068EZd++5IrV6TdB0HAAAAyEqeLLSNzREVBZjQIrP99oQz1KUsqONuuMt1\nFAAAACArebLQRlootMh8PuPTX0+fqtcD1+iFN5a7jgMAAABkHU8W2mg8quI8lhwj8+2/4491ZPcR\nOvkvlyjBV9MCAAAAncqThTaWiCiUz4QW3vDoqGsVrZ6hX936iusoAAAAQFbxZqFNUmjhHaH8Yt18\nyGTdsWiU5i9sch0HAAAAyBqeLLTNqahKC1hyDO8YffAx2rH8Rzripuv5bloAAACgk3iz0NqIyoqY\n0MJbnj//T/qq/C5N+MsHrqMAAAAAWcGThbbFRlVayIQW3tKnsoeu3uv3uu6DM7R4aYvrOAAAAIDn\nebLQxhVRRYgJLbznumNP1bYVvTXktxNdRwEAAAA8z5OFNuGLqAuFFh5kjNELY/6s2aW36ZaHPnMd\nBwAAAPA0TxbapC+qihBLjuFN/br21iW736jL3xyh5Sv5cloAAABgc3my0Kb8EVWWMqGFd930f2ep\ne0WpDhs/yXUUAAAAwLM8WWhtIKqqMia08C5jjF4ec49mldyiCfd+6DoOAAAA4EmeK7SJZFLyxdWl\npMB1FGCL7NR9W43f54+65qNfaM5XMddxAAAAAM/xXKFdE4lI8SLl5RnXUYAtdtUxJ6t/1W4aPOFy\npVKu0wAAAADe4rlCu2ptVCbBcmNkj1cumqK6qr9r5KSXXEcBAAAAPMVzhbZubUS+JDeEQvaoClVo\n2jH3alrdCL323irXcQAAAADP8FyhXd0YkT/FhBbZ5Rf7HKzDeg7T0XeerUjEuo4DAAAAeILnCu2a\ncFSBFBNaZJ+nxtyoYNXXOuSq211HAQAAADzBe4U2ElFAFFpkn4Jgvqaf95jeK7xev7n7fddxAAAA\ngIznuULbEI0qTyw5RnbatXdfTfzfKRr/+TC990mD6zgAAABARvNgoY0o3zChRfa65IgT9L89huiQ\nyWdxPS0AAADwAzxXaNc2RZTvo9Aiu/3jwj/IXzVfh101xXUUAAAAIGN5rtCGm6Iq8LPkGNmtMFig\nV0c9rncKx+vqO952HQcAAADISJ4rtI3NERUGmNAi+w3Ydgfddsg0TZh/gp6fsdR1HAAAACDjeK7Q\nRuNRFQWY0CI3jD74aJ3Y92wd+/D/adGSFtdxAAAAgIziuUIbiUdUnMeEFrnjoXOvUZ9ulfrpdReq\nhU4LAAAAfMtzhTaWiCiUT6FF7vAZn965/AGtrXxVR1w5zXUcAAAAIGN4rtA2JaMK5bPkGLmlvLBU\nr57zd70WHKuxU950HQcAAADICN4rtKmISguY0CL37LX9zppy2F/0uwXH64HnvnIdBwAAAHDOc4W2\nORVRaRETWuSmcw46UmN2u1pnvPQzvf3xGtdxAAAAAKc8V2jjiqq8kAktctetvzhPg2sO10FTT9Di\npdwlCgAAALkr7YXWGDPEGDPbGDPHGHP5en6+kzHmTWNMkzHm4o2dL66IyosptMhtz1/4B/XuVqQ9\nrh2tWMy6jgMAAAA4kdZCa4zxSbpd0uGS+ksabozZeZ3DVkk6X9KkTTlnwkRVEWLJMXKb3+fXzLGP\nqKVqpva6ZIKSSdeJAAAAgK0v3RPavSTNtdYutNbGJT0qaWj7A6y1ddbamZISm3LCpD+iLiVMaIHS\ngpBmXvyc5pXerUMuu1uWQS0AAAByTLoLbS9Ji9q9Xty2b7Ml/RFVlVJoAUnavrqn/jXyRb0RvFYn\njX/KdRwAAABgqwq4DtAR48aNk/1Xg+5N3qKjjzxUgwYNch0JcG5gn376xynP6oiHhuhXt1To1osG\nuY4EAAAAbFRtba1qa2u36BzGpnGdojFmb0njrLVD2l5fIclaayeu59jrJDVaa2/ewLlsIplQ4DdB\nNY1NKj/fpC034EUPvzVdpz49TL/d+UWNPX1313EAAACADjHGyFrboaKX7iXH70nawRhTY4zJk3SS\npGd+4PgfDL82FpPihZRZYD1O3ucg/X7QHbr6859pyhOfu44DAAAApF1alxxba5PGmDGSXlJreZ5m\nrf3CGHNu64/tncaYbpLel1QiKWWM+ZWk/7HWhtc936q1UZkEdzgGNuSiIcepsSmq8/91qAL+V3TO\nceveVBwAAADIHmm/htZa+4KkndbZ9+d2z5dL2mZTzrWmMSaTLOzcgECWufbYUxRPJjXqrUMU9E/X\nGUP7uY4EAAAApIWnbgq1JhKVL8WEFtiY648/TYlkQmfNOFgB/3SdelRf15EAAACATuepQlsfjilA\noQU2yYQTz1Q8mdDp0wfLmOk65Wfbu44EAAAAdCpPFdqGaFR+y5JjYFP9fvi5Slmr06YfqHDsJY08\nYRfXkQAAAIBO46lCuzYaU9AyoQU64uaTR6q0oFjnvT1Y9ZHndMVpA11HAgAAADqFpwptQyyqoGFC\nC3TUuONOVWlhsS59bYgapj6lCaP2cx0JAAAA2GLp/h7aTtXYFFWeYUILbI6LjzhOdx7xgH638FiN\nmvSy6zgAAADAFvNUoQ03xZTvY0ILbK4zDzxcjx33N9216hQddeWDSqVcJwIAAAA2n7cKbXNU+X4m\ntMCWOGGvA/TaiFf1cuoqDbjgRjU1WdeRAAAAgM3iqUIbbYmpgEILbLH9duyvzy9+S1+HntB2F4zU\nylUJ15EAAACADvNUoY3EoyoMsOQY6Ax9u/bUwuteV37XBdr+ymM1a27YdSQAAACgQzxVaKPxqIqC\nTGiBzlJWWKK545/Vj/p014DJ++qxF//tOhIAAACwyTxVaJsSMRUFmdACnSnoD+rNK+7S2T85Sye/\nuo8uum2660gAAADAJvFUoY0loyrOZ0ILdDZjjKb88gLde/RDmrxkuPa/+E9qbuZmUQAAAMhsniq0\nLcmYQnkUWiBdfrn/wZo5+k19kjdVfc4/RwsWN7mOBAAAAGyQpwptcyqqUAFLjoF02m3b7bXourfU\npVeD+k3cRw88N891JAAAAGC9PFVoW2xMpYVMaIF0Kyss0WfXPqaRe52l02fsqxPHP6Fk0nUqAAAA\n4Ps8VWjjiqqECS2wVRhjNPnU8/Tc8Of1/2JXqM/oMfp6SbPrWAAAAMC3PFdoy4qY0AJb05Dd9tCi\nq2equPsS9Z2wt6b8dZbrSAAAAIAkjxXapImprJgJLbC1VYXK9cW4J3Xhvufp/JmDtN8lf1RjOOU6\nFgAAAHKcpwptwhdVeTETWsAFY4wmDT9L7498S3PzHle3Xx+qZ2oXuY4FAACAHOapQpvyxVQRotAC\nLu1es4OW3PC6hu42WMc+v4eOueY+xWJ8Zy0AAAC2Pk8VWhuIqiLEkmPAtYAvoEdGXaXnf/GCXo//\nUVUXHa7HXvq361gAAADIMZ4rtJWlTGiBTHH4rgO18oZ3dfzuh+jk6Xtq30tu1qo1CdexAAAAkCM8\nVWglqTQUdB0BQDtBf1D3n3uZ3jv3bS3Mf1Y9rtlHv3twpiyrkAEAAJBm3iq08SL5vJUYyBkD++yg\nxTe8oosPGK2rPjtKvUedoxkfrHQdCwAAAFnMU/XQJFhuDGQyY4xuGnaGll79hXbsU6wDH+uvgy6f\nrLrVLEMGAABA5/NUofUluSEU4AVVoXJNv+IWvTZiuub4/q4e43bXhZNfVoJeCwAAgE7kqULrTzGh\nBbzkgJ36a/GN/9QNh43TnYvPU9n5h+r3D3F9LQAAADqHtwqtZUILeI0xRpcddbzqb5ilEXsfr7Gf\nHq3qUcP1yIvzXUcDAACAx3mq0AYsE1rAq/ICQU0+baRWjZujQf3765Tan6pm9Cj9vXah62gAAADw\nKG8VWjGhBbyutCCkv55/tRZe9oV27lOu418cqN6jz9bjL3/lOhoAAAA8xlOFNs8woQWyRe+Kar14\n2QQtvnyOBvTtruGv7qkeo07X/c/N4RpbAAAAbBIKLQCnepRX6tlLrteyK+dr752214g39lOX0T/X\nuHtmqKWFZgsAAIAN81Shzfex5BjIVtUl5Xrqwmu1+toFOvEnh2rilyNUcvFPdcpNj/I9tgAAAFgv\nTxXaAj8TWiDblRYW689njlZ4wmzdOOQqvbRmirrd0Ff7XnaTXpu53HU8AAAAZBBvFdoAE1ogV/h9\nfl1y1FCtmPi6nj7lSUUL5+qgJ3dS9ehhuuqu6YrFWI4MAACQ6zxVaAsDTGiBXHTU7j/RR+OnafnY\nBRo64ADdNvcClYzdWYPG/kEzPmRqCwAAkKs8VWiLghRaIJdVl5Tr7nPGaO3ET/TQifeoLvCJDnxi\nJ5WN/plOm/SoFi6JuY4IAACArchbhTaPJccAJGOMhu27nz67/j6tuWaxRu4/XC/X3aPtJvfSNmPO\n1LX31KoxnHQdEwAAAGlmrEe+8NEYY0+9eZruv2iE6ygAMtTcZUt07V8f1nNfP6iwXa7tmn+u4bsd\nr0tOOFAVZQHX8QAAAPADjDGy1poO/Y6XCu3IPz2sqaOHu44CwAPemTdXk559Uv/8z5Na61ugbaJD\n9X/9T9ClJwxWt6o81/EAAACwjqwvtJdN+7smjhjqOgoAj/nk6wW66ekn9cLXT2pN4HNVrh2sA3se\nqVGHHKmD9+op06E/NgEAAJAOWV9of/PQi7rm5MNcRwHgYQvrVmjy8y/qmS/+ofl6Uf7GGvXPP1LD\ndj9SZx+5lyorgq4jAgAA5KSsL7S3PjVDFxy7v+soALJEPJnQ42++rXvf+IfeWf28wsGvVNawn/bo\nMlgn7DFYpxyym0pCftcxAQAAckLWF9q/vDhTpx020HUUAFlq8eo6Tfvna3p21nTNir6qmH+ZqsKD\ntGf1QTpmwL4aNmhXVZQxwQUAAEiHrC+0T834XMfuv4vrKAByxLxlS3X3K7V6fvZ0zYu9pWj+AoUa\nfqIdi/bRoL77avgBe2uPXaq4BhcAAKATZH2hrf1wgQ4cUOM6CoActWJtvR55/R3949M39dHqN7Uy\n7135It3VPbWn+lcO1KB+e+jYfQZo5z5llFwAAIAOyvpC+9Hc5dpth66uowCAJCmRTOqVT2fp6fff\n17tff6D50Q9Un/+J/NEeqk4O1C5lA7X/DgN1xMAfa89duikQoOUCAABsSNYX2q+XNWqbbiHXUQBg\ngxLJpF7//Es9O/MDvbXgA80Jf6D64GdKWaviSH/1Cv5IO1f210+3669DB/TXHjtXy+dznRoAAMC9\nrC+0kWhCRYXccRSAt1hrNX/5cr3wwSy9OXeWPlsxS4uaZqkhf5ZsIqhQ0y7q6u+nmpId1L9HP+3V\nt58O+FFf1fQoZukyAADIGVlfaFMpy1/uAGQNa62+XLJUr378pd5fMFezV8zVoshcrbLzFCv4Sqa5\nQsXN/dTVv4NqSvqqb1WN/qfXthq4fY322LGnQsX8Dz4AAJA9sr7QeiUrAGypZCqlz75erBmz5mnm\ngrmas/IrLY18rVXJhYoEFyqZVydftKeK49uqwlejHkU12q6iRjt07aV+PXroRzU99T99qlSQz3pm\nAADgDRRaAMgRTfEWfTh/sWbOX6hZixdqXt1CLW5cqFUtSxQ2S9QcXKpUXoN8sa7Kb+mhkHqqItBD\nXYt6qHdpT/Wp6qE+1dXq061KO/aq1jZdS7hpFQAAcCojC60xZoikP0rySZpmrZ24nmNuk3SEpIik\n0621H63nGAotAHRAU7xFn3+9TLMWLtWcZUv175VLtLhhqVZEl2p1fIkitk5N/pVK5NVJvhb5mqoU\njFepIFWtkK9KZcFqdSmoUrdQtbqVVqp7WYW6l5erZ5dybVNVoW26lqm8JI9LQQAAQKfIuEJrjPFJ\nmiPpYElLJL0n6SRr7ex2xxwhaYy19mfGmJ9KutVau/d6zkWhBRyora3VoEGDXMdAmjXGYpq/dJXm\nL12pBSvqtGj1Si1pqNOK8EqtbqpTfUudoql6NWmNWnz1SgTqZfPrpUS+fC3lCiTKFUyVq0DlKvKV\nKxQoV2leuUrzS1WaH1JpYUjlhSGVF4dUGSpRZUlI1WWtW7eKElWWFjAhXg/ef4AbvPcANzan0AbS\nFabNXpLmWmsXSpIx5lFJQyXNbnfMUEn3S5K19h1jTJkxppu1dnmaswHYBHyo54aSwkIN2L63Bmzf\ne5N/x1qr1eGIFq2s1+K6ei1ZXa9lDfVaubZedeF6rYquUWPTai2LLVR0VVhNybCabFgtNqy4r1EJ\nX1hJf1g2EJb8LVI8JF88JH8ypECqRAFbpKAKFTQFCppC5fsKle9v3QoDhSoItD4W5RWqOK9QxfkF\nChUUqqRtKy1q3UIF+SrKz1OoME/FBa1b6/Og/P7MLtG8/wA3eO8B3pHuQttL0qJ2rxerteT+0DH/\nadtHoQWADGaMUWVJSJUloQ4V4fVpSSS0siGsFfVhrWwIa2VDo1aHowo3xdTYFFO4KaZwc0zRliZF\nW2KKxmOKJSKqj9apuTGm5mRMzamYWmxMCdukuGJKmJiSJqaUaVHKtMj6vtvkb5H8cSkZlJJ5Mqk8\nmWR+66PNk6/d5td3W0B5Cpg8+U1QPuOX3wS+3QImIL8vIL/xK+Br2+drv/kV8Lc+D/oCCvgDCn67\n+ds9Dygv0Pr4r1lf6Q9/my6/z/ftFmj33O9vfQz6/a0/87c7zr/+x7yAX35/6+tA2zkC7V5/8zwY\n8CvgN/L5Mrv0AwByW7oLLQAAG5UXCKhXZbl6VZZvtX9mKmUVa2lRuKlFkViLIk1tW3OLou22WHOL\nYi1tW7xFTfFmNccTSqQSiieTiie/eZ5QIplQPJVQMpVse0yoJdWkaCKhhG19nUwllbAJpWxCye9t\nSaWUUEqtr1OKK7Jont6fsUhWqQ1sSVmTktrtk/n+Y+vPvztOJvXfz7957Uu2Pv92s5I1kvW1Psq0\nPard83aP69ln1j32ez9bZ78k2bb96zu27fU3Z/3mmHXPs6Fj25/3u/1qd65196yzz5r1HPb9HWZ9\nx6yzeu6//9nrOc8m5PvvozbvmM36d7HRvG3HtNu9oX83mXgdfuStuZrS8K7rGFtmvf+9onOs/794\nuJHuQvsfSdu2e927bd+6x2yzkWMktU4DAGx948ePdx0ByFlN7yxwnMBKSm7Rb2/OzwDXou/Mcx0B\nwCZId6F9T9IOxpgaSUslnSRp+DrHPCPpPEmPGWP2llS/vutnO3pxMAAAAAAgu6W10Fprk8aYMZJe\n0ndf2/OFMebc1h/bO621/zDGHGmMmafWr+05I52ZAAAAAADZIe3fQwsAAAAAQDr4XAfYFMaYIcaY\n2caYOcaYy13nAXKFMWaBMeZjY8yHxhiP3x0DyFzGmGnGmOXGmE/a7aswxrxkjPnSGPOiMabMZUYg\nW23g/XedMWaxMeaDtm2Iy4xAtjHG9DbGvGqMmWWM+dQYc0Hb/g5/9mV8oTXG+CTdLulwSf0lDTfG\n7Ow2FZAzUpIGWWt3t9au+5VbADrPvWr9nGvvCkn/tNbuJOlVSWO3eiogN6zv/SdJN1trB7ZtL2zt\nUECWS0i62FrbX9I+ks5r63gd/uzL+EKr1u+tnWutXWitjUt6VNJQx5mAXGHkjT8nAE+z1r4hac06\nu4dKuq/t+X2Sjt2qoYAcsYH3n8SX3gBpY61dZq39qO15WNIXav22mw5/9nnhL6q9JC1q93px2z4A\n6WclvWyMec8Yc7brMECO6frNXf+ttcskdXWcB8g1Y4wxHxlj7mbJP5A+xpg+kgZIeltSt45+9nmh\n0AJwZz9r7UBJR6p1Kcj+rgMBOYy7OAJbzxRJ21trB0haJulmx3mArGSMCUn6q6RftU1q1/2s2+hn\nnxcK7X8kbdvude+2fQDSzFq7tO1xpaSn1HoJAICtY7kxppskGWO6S1rhOA+QM6y1K+13XwVyl6Q9\nXeYBspExJqDWMvuAtfbptt0d/uzzQqF9T9IOxpgaY0yepJMkPeM4E5D1jDFFbf/XTMaYYkmHSfrM\nbSogqxl9/5q9ZySd3vb8NElPr/sLADrN995/bX+R/sZx4vMPSId7JH1urb213b4Of/Z54nto226V\nfqtaC/g0a+1NjiMBWc8Ys51ap7JWUkDSQ7z3gPQwxjwsaZCkSknLJV0n6e+SnpC0jaSFkk601ta7\nyghkqw28/w5S6zV9KUkLJJ37zXV9ALacMWY/Sa9L+lStf9e0kq6U9K6kx9WBzz5PFFoAAAAAANbl\nhSXHAAAAAAD8FwotAAAAAMCTKLQAAAAAAE+i0AIAAAAAPIlCCwAAAADwJAotAAAAAMCTKLQAAAAA\nAE+i0AIAkOGMMcZ1BgAAMhGFFgCADGOMqTHGzDbG3GeM+VRSb9eZAADIRMZa6zoDAABoxxhTI2m+\npH2ste+5zgMAQKZiQgsAQGZaSJkFAOCHUWgBAMhMEdcBAADIdBRaAAAyEzeCAgBgIyi0AABkJm5y\nAQDARnBTKAAAAACAJzGhBQAAAAB4EoUWAAAAAOBJFFoAAAAAgCdRaAEAAAAAnkShBQAAAAB4EoUW\nAAAAAOBJFFoAAAAAgCdRaAEAAAAAnvT/AevuR+tUu+dfAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x10762b3c8>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# normalize u & u_exact are the values of a wave function normalized to norm 1\n",
+    "u /= np.linalg.norm(u) * sqrt(dr)\n",
+    "u_exact /= np.linalg.norm(u_exact) * sqrt(dr)\n",
+    "\n",
+    "# plot\n",
+    "plt.figure()\n",
+    "plt.title('Hydrogen wavefunction ($l=0$)')\n",
+    "plt.plot(rs, u, label='numerov')\n",
+    "plt.plot(rs, u_exact, label='exact')\n",
+    "plt.xlabel('r')\n",
+    "plt.ylabel('u(r)')\n",
+    "plt.legend(loc='best')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# (b) Poisson Solver"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Verlet Integration"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Integrate $U''(x) = f(x)$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "def verlet(fs, U0, U1, dx):\n",
+    "    \"\"\"compute U = [U0, U1, U2, ...] for fs = [f0, f1, ...] via verlet integration\"\"\"\n",
+    "    dx_sqr = dx**2\n",
+    "    U = np.zeros_like(fs)\n",
+    "    U[0] = U0\n",
+    "    U[1] = U1\n",
+    "    for i in range(2, len(fs)):\n",
+    "        U[i] = 2 * U[i-1] - U[i-2] +  fs[i-1] * dx_sqr\n",
+    "    return U"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Integrate $U''(r) = -\\frac{u^2(r)} {r}$ subject to boundary conditions $U(0) = 0$ and $U(r_{\\max}) = 1$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "def solve_poisson(u):\n",
+    "    # start verlet integration from U(0) = 0, U(dx) = dx (can always be extended to a solution, but does not necessarily satisfy the right-hand side boundary condition)\n",
+    "    fs = -u**2 / rs\n",
+    "    U0, U1 = rs[0], rs[1]\n",
+    "    #U0, U1 = 0, dr\n",
+    "    U = verlet(fs, U0, U1, dr)\n",
+    "    \n",
+    "    # fix right-hand side boundary condition\n",
+    "    U = U - (U[-1]-1) / rs[-1] * rs\n",
+    "    return U\n",
+    "\n",
+    "U = solve_poisson(u)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Exact solution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "U_exact = -(rs + 1) * exp(-2 * rs) + 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7QAAAIwCAYAAAC2kSfUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VOW9x/HvLwlJgCRAgJAAWdg3cd83gtalFavVahGt\nS6u3taW19drq1bait71a23vrgtVqtS5YUdvrvqCiEVRAQUVZwhaTkI0t7CH7c//IwI1pCEmYmTNz\n5vN+vXgxmXPmzPfMwEu/PM95jjnnBAAAAABAtInzOgAAAAAAAN1BoQUAAAAARCUKLQAAAAAgKlFo\nAQAAAABRiUILAAAAAIhKFFoAAAAAQFSi0AIAEOHM7AEzu6WT+75rZt/rYPt4M/s4eOn+5fi3mtmT\noTp+d5nZRDP7wOscAIDgotACAA6KmX1pZqe1ee4KM5t/EMc8qNeHk5nlmlmzmQXlv6ntnbtz7lrn\n3O+CcXxJt0u6q9X7Bf37kxRxN7l3zn0haauZneN1FgBA8FBoAQCh0q1SY2bxkuxArw9WgQyCvVkt\nyMcLOjPLlJQv6cVO7B70DIHv1kt/l/RDjzMAAIIoUv5nAADgY2Z2o5mtNbMdZrbMzM5vte0KM3vf\nzP7HzDZLmi3pAUknmNlOM6sO7Pc3M/uzmb1qZjsl5ZtZopn90cxKzKwysD2p1bGnmNmnZrY18B4T\nO8jYbGY/MbN1ZrbRzFqPYpqZ/crMis2sysweM7PUwOb3Ar9vC5zfcYHXfM/MVpjZFjN73cxy2rzX\nD8xstZlVm9nMwPNjOzj32wOP+5rZy4GMWwKPh3TyqzhD0ifOufpO7i8zu8HM/tHmuXvN7E+Bx3lm\nVmBm281sjqQBrfbbO3r9PTMrkTQ38Pw3A38Oqs3sncB5733NkWb2SeB4z5rZ7L3nHti+3+80MNr8\n72a2NLD9aTNLbBW9QNLpZtajs+cPAIhsFFoAQCi0Ha1cK+kk51yapNskzTKzQa22HxfYJ0PSZWoZ\nRVvgnEt1zqW32u8SSf/pnEuV9IGk30saKenQwO9DJP1GkszsCEmPSLpGUrqkv0h66QBl5nxJRwZ+\nnWf/fy3qVZIulzRJ0nBJqZLuD2w7NfB7mnMuzTm3yMzOk3RT4HgDJc2X9HSb9zpH0lGSDpN0sZmd\n6Zwr7ODc94qT9KikbEk5kmokzezgnFqbKGlVJ/Zr/f3NknSWmaVJ+0ZZvyPp8cD2v0v6WC1F9reS\nrmjneKdKGhs4zqjAa36qls/mdUkvm1lC4Lv538D5pavlM/vWvlCd+04vknSmpGFq+Wyv3LvBOVch\nqUHSmE58BgCAKEChBQAEwwuB0bbqwKji/a03Ouf+6ZzbEHj8nKQ1ko5ttUu5c+7Pzrlm51xdB+/z\nonNuYeA4dWopNj93zm13zu2WdKdaSq8C2x50zi12LZ6UVCfp+A6Of2fgWGWS7m51rGmS/sc5V+Kc\nq5H0H5KmBqY97y1/rUvgDyTd4Zxb7ZxrDuQ63MyyW+1zh3Nup3NuvaR3JR3eQa59nHPVzrnnnXN1\ngXO+Q/9fqg+kr6Sd7Ty/3+/POVellkJ+UeCpr0va5Jz7LDDqfLSk3zjnGpxz8yW93DaypFudc3sC\n39l3JL3inHvHOdck6Y+SkiWdqJbvJt45N9M51+Sce17SR62O1Znv9B7n3Abn3LZAlraf687A5wAA\n8AEKLQAgGM5zzqXv/SXpR603mtnlraaJbpU0Qa2mpkpa38n32befmQ2U1EvSklZF7HVJ/QO75Er6\n91ZFbaukoZIGd3D8slaPS1rtOzjwc+ttCZIGqf1rTXMl3dMq15bAfq2nBm9o9bhGUkoHufYxs55m\n9pfA9Odtapny3NfMOnMN71a1jC631eH3J+kJtYycS9KlkvauYpwlaatzbk+rfVt/Tnu1/ly/8lk6\n51xg+5DAtvI2r239Z6Mz3+mBPtdUSdvayQgAiEIUWgBAMOy3TAVG8R6S9CPnXD/nXD9Jy9u8pm0p\n3N+CRK2f36yWwjKhVRnr65zrE9i+XtLvWm3r55xLcc4908F5tB5BzZVUEXhcEfi59bYGtZSn9rKW\nSvpBO++9sIP3bu8c23ODpFGSjnHO9dX/j852ptB+Lml0O88f6LUvSDrUzCZImiLpqcDzlZL6mVnP\nVvvmtH2xvnpObT9LqeVzLw8cb2g72/bqzne6j5kNltRDnZt2DQCIAhRaAECo9ZbULGmzmcWZ2VWS\nDjnAazZIGtrR9a6Bkb2HJd0dGK2VmQ0xszMDuzws6YdmdmxgW28z+4aZ9e7gfX8RWHQpWy3XeM4O\nPP+0pJ8HFkBKkfQ7SbMD04k3Bc5vRKvj/EXSzWY2PvDefczs2wc4586ee4qkPZJ2mFm6pBmdPK4k\nvSXpyDYLJR2Qc65WLde2/l3SosCUbDnnSiUtlnSbmfUws5Mlndvm5W3L8rOSzjGzyYHrZm+QVCvp\nQ0kLJDWa2Y/NLD5wLXLrqend+U5bmyTpHedcQ2fPHQAQ2Si0AICD1eGIonNupaT/lrRQUpVaphu/\nf4BjvqOWUdwqM9vYwX43qmUxqYWB6bdvKjAC6ZxbopZrLmcGpv2uVvsLFrX2oqQlkj5Ry/WXjwae\nf1Qt02znSVqnlpHhnwbeZ49aCu4HgWmwxzrnXlDLdbOzA7k+l3R2q/fpaET6QOd+t1qmWm9WSwl8\nrYNjfXWDcxsDxz+/9dP727+Nx9WyqNQTbZ6fppZrWLdI+rX+f7Godo/vnFutlunLM9XyjwHnSDrX\nOdcYKJoXSLpaLdOjp6nle6gLvPZA3+mBzuVSSQ8e6EQBANHDWv6BO0QHN3tELVOTNjjnDt3PPveq\nZYGJ3ZKudM59FrJAAADsh5k1SxrpnCvyOksomdk4SY85547r4uuGSiqUlOmc2xWScO2/70JJDzjn\n2hblrh5noloWlDopOMkAAJEg1CO0f5N01v42mtnXJY1wzo1Sy4qQ/KspAAAh5Jxb2Y0yG6eWa3dn\nh7rMmtmpZjYoMOX4CrWMCr9xsMd1zn1BmQUA/0kI5cGdc++bWduFH1o7T4GpS4H79vUxs0F7b+0A\nAEAYhW7KUhQzs15qua73S7XMqAq1MWq5zraXpCJJF/L/BQCA/Qlpoe2EIfrqcvzlgef4DxcAIKyc\nc/FeZ4hEgfvutnern1C938NqWfwJAIADYlEoAAAAAEBU8nqEtlxfvb/cUP3rDdUlSWbGVDAAAAAA\n8DHnXGfuq75POAqtaf83bH9J0o8lPWNmx0va1tF1MqFckRmxwzmnsi3b9NGqEq2uqNKXm6pUvm2D\nNuyu0ua6Ku1oqlJNXJXqe2ySS9omNfRSfH0/9WhMV5Lrp17WTykJ6Urr0U9pyWlKS0pRWnLLr769\nUtSvd2/1S0nRgNQUpaf2VkbfFGX07a203omKi+vS38+IMGPGDM2YMcPrGEBM4u8f4A3+7gHeMOv6\n/yuHtNCa2d8l5Uvqb2alkm6VlCjJOececs69Frgh+lq13LbnqlDmQeyoqavXe8vWav6KQi2vKFLJ\n9hJtqC3WNpWotmex5ExJe3LUu3mw+sRnqn9yprJScnTMkGM0bGCmRmZmavSQgcrL7KuUXj28Ph0A\nAAAA7Qj1KsfTOrHP9FBmgL/VNzbqrU9X6rVPP9Gn5StUvKtQW2yl6nuWKmF3rvo1jdXg5BHK7TtS\nZ438mg7NzdWxo/M0cmhfdeMfgAAAAABEEK+voQU6rbnZae7nK/XMhx/q47JP9GXtEu3stUwJu7OV\n6Y7QqD4TNXXsFTppzFidfsRI9U1N9DqyL+Tn53sdAYhZ/P0DvMHfPSB6WLRcl2pmLlqyIjgam5r0\njw8/0bOL5umjqvmq7PG+rK6PhjSfpMMzjtbksUfq/BMOU15W2O4mAQAAACBEzKzLi0JRaBFRiqo2\n6Z5X5+jV1a/ry7g5it8zSCMS8jUp7xRdduopOmniEKYKAwAAICLk5eWppKTE6xhRJzc3V8XFxf/y\nPIUWUWlZSbl++/xzemP9s9qeuEIZuycrf+g3dO3XztakI7IpsAAAAIhIgQLmdYyos7/PjUKLqFG2\npVo3P/20Xvlytrb1WKHc2vM07bCLdcOFp6lfGte+AgAAIPJRaLuHQouo5JzTI2/P011zH9bauFeU\nuevr+u5hl+nGb5+h9D6UWAAAAEQXCm33BLPQssoxQm5Xba2uf/wJPbn2T2pqiNPp/a7RP6bdo0NH\n9fc6GgAAAIAoRqFFyKzfXK1rH/2z3qieqbRdR+uW4x/QTVMnKSGBi2IBAACAaPH444/rr3/9q+bP\nn+91lH9BoUXQbdqxQ5c/8CfN2Xqfsvd8U09MeUfTzhjvdSwAAAAA3WSdXKn1qquuUnZ2tm6//fYQ\nJ2pBoUXQ7K6r1VUP3K9/Vt2lrN1n6aWpH2nKScO9jgUAAACgm5qamryO0KE4rwPAH/744qtK/9Uh\nmrt2nh4/7R2V3fcEZRYAAADwyF133aWLLrroK89dd911+tnPfqYdO3bo+9//vgYPHqzs7Gz9+te/\n3rdI0+OPP66TTz5Z119/vQYMGKDbbrvtX45dWFioM888U/3799e4ceP03HPPSZIefvhhPfXUU7rr\nrruUlpam8847L+TnyQgtDspHa4p04V+vU2X9Kt0w/n7dcfVZ3DcWAAAA8NjUqVN1++23a/fu3erd\nu7eam5v13HPP6YUXXtCVV16pzMxMFRUVadeuXZoyZYpycnJ0zTXXSJIWLVqkadOmaePGjWpoaNDs\n2bP3HbempkZnnnmmfvvb32rOnDn6/PPP9bWvfU0TJ07UNddcow8//DCsU44ZoUW3NDU366oH7tPx\njxyrvPiTVPWbL3TnNZRZAAAAYC+z4PzqjpycHB155JF6/vnnJUlz585V7969lZeXp9dee01/+tOf\nlJycrAEDBuhnP/uZnn766X2vHTJkiH70ox8pLi5OSUlJXznuK6+8omHDhunyyy+Xmemwww7ThRde\nuG+UNtwYoUWXLV73pb7x4Pe0q7ZOz17wob49ebTXkQAAAICI4/Utai+55BI9/fTTuuyyy/T0009r\n2rRpKikpUUNDg7KysgIZnZxzysnJ2fe67Ozs/R6zpKRECxcuVHp6+r7XNzU16fLLLw/tyewHhRZd\n8utnntF/fTpdJ/f4pV677Xr17hXvdSQAAAAA7bjooot0ww03qLy8XM8//7wWLVqktLQ0JScna8uW\nLftdubijFY2zs7OVn5+vOXPmdPm1ocCUY3TKnvo6HXfbdN2x6Bbdf/ybeu+/fkGZBQAAACLYgAED\nNGnSJF111VUaPny4Ro8erczMTJ155pn6+c9/rp07d8o5p6KiIs2bN69Tx5wyZYpWr16tWbNmqbGx\nUQ0NDVq8eLFWrVolSRo0aJCKiopCeVpfQaHFAX1Rul6Zt5ysoo2VWvbTJfrh+Ud4HQkAAABAJ0yb\nNk1z587VpZdeuu+5J554QvX19Ro/frzS09N10UUXqaqqqlPHS0lJ0ZtvvqnZs2dr8ODBGjx4sG66\n6SbV1dVJkr7//e9r+fLlSk9P1wUXXBCSc2rNnNcTuzvJzFy0ZPWT5xct1sX/e76O13V657c3qEcP\nVn0CAAAApJbptXSUrtvf5xZ4vkuFgxFa7Ndv//m/uvD5r+v7WTM1//e/oMwCAAAAiCgsCoV2Xf3Q\nTP1tzZ2674Q39OMLjvI6DgAAAAD8Cwot/sVF996h54v/qhcvmK8pJw/zOg4AAAAAtItCi32cczr7\nD7fonYoX9PZ35yn/qCFeRwIAAACA/aLQYp+v//EWFZS9rgXXvqejxw30Og4AAAAAdIhCC0nSxffe\nobnlL2rhte/pqHEDvI4DAAAAAAdEoYWuefh+/W/xI5p7+TzKLAAAAICoQaGNcTOee1aPrrpTL14w\nT5OOHOx1HAAAAADoNAptDHvmwwW6ffF0PXjKW6xmDAAAACDqxHkdAN5YUvSlLn3xQk3P/pv+7bzD\nvI4DAAAAIIo9/vjjOuWUU8L+vhTaGLR55w6d+sAU5Sf8h+6dfo7XcQAAAABEOeeczCzs70uhjTHO\nOZ1w5/fVv+YkvfmfP/E6DgAAAIAQqays1Le//W1lZGRoxIgRmjlzpiTpnHPO0Q033LBvv6lTp+rq\nq6+WJBUVFen000/XgAEDlJGRocsuu0w7duzYt29ZWZkuvPBCZWRkaODAgfrpT3+qwsJCXXvttVqw\nYIFSU1OVnp4etnOk0MaYKx68W6U7vtSiGfcqjm8fAAAA8CXnnM4991wdccQRqqys1Ny5c3X33Xfr\nrbfe0qOPPqpZs2apoKBATz31lBYvXqx777133+tuvvlmVVVVaeXKlSorK9OMGTMkSc3NzZoyZYqG\nDRum0tJSlZeXa+rUqRo7dqwefPBBnXDCCdq5c6eqq6vDdp4sChVDnpr/gWYV36mXvrNIWQOTvY4D\nAAAA+JrdFpwpuO5W1+XXfPzxx9q8ebNuueUWSVJeXp6uvvpqzZ49W2eccYYeeOABXX755aqtrdWL\nL76oXr16SZJGjBihESNGSJL69++vn//857r99tslSYsWLVJlZaXuuusuxQVGx0488cRgnGK3UWhj\nRFn1Fl35ylRNz35UU07O8zoOAAAA4HvdKaLBUlJSovLy8n3Tf51zam5u1qmnnipJmjJliqZPn64x\nY8bohBNO2Pe6jRs36rrrrtP8+fO1a9cuNTU17TtGWVmZcnNz95XZSBA5SRAyzjlN/u8fatieb+ue\nH7MIFAAAAOB32dnZGj58uKqrq1VdXa2tW7dq+/btevnllyVJN998s8aPH6/KykrNnj173+tuvvlm\nxcXFafny5dq2bZtmzZol59y+Y5aWlqq5uflf3s+LBaEkCm1MuOnvs1S8a6UKfnWHPPpzBgAAACCM\njj32WKWmpuquu+5SbW2tmpqatHz5ci1evFjz5s3T448/rieffFKPPfaYfvKTn6iyslKStHPnTqWk\npCg1NVXl5eX6wx/+8JVjZmVl6aabblJNTY3q6ur04YcfSpIGDRqksrIyNTQ0hPU8KbQ+t7S4RH/8\n4nrdN3mWBmdw3SwAAAAQC+Li4vTKK6/os88+07Bhw5SRkaFrrrlGlZWVuvLKK3X//fcrMzNTJ598\nsq6++mpdddVVkqRbb71VS5YsUd++fXXuuefqwgsv/MoxX375Za1Zs0Y5OTnKzs7Ws88+K0k67bTT\nNGHCBGVmZiojIyNs52l7h48jnZm5aMkaKZxzyrrpdI3QWfrg9zd6HQcAAADwFTMTHaXr9ve5BZ7v\n0pxSFoXysX+f9Zi21uzU6/91w4F3BgAAAIAoQ6H1qS83btQ9y2/SvZPfUFpqvNdxAAAAACDomHLs\nU+N+dals5xCtuOcur6MAAAAAvsSU4+5hyjE6NPONOVpds1BFN37hdRQAAAAACBlWOfaZusZ6/eKd\nn2r68PuUO7iX13EAAAAAIGQotD5z9UMzlbhrhP7n2m94HQUAAAAAQoopxz5SsnmT/l56hx7/5nzF\nsw4UAAAAEFK5ubky69Iln1DL5xYsLArlI0ff9kPtqE7W6nvu9joKAAAAAHQJi0LFsLe/+EKf1Dyv\npdMLvY4CAAAAAGHBNbQ+cfVTv9LkhP/QxFH9vI4CAAAAAGHBCK0PPLdggdY3fqoPr3/G6ygAAAAA\nEDaM0EY555ymv3Czvpl2qwZnJHsdBwAAAADChkIb5R56e6621Jfrbz+7wusoAAAAABBWFNoo5pzT\nTW/doksH/6f6pjF7HAAAAEBsodBGsQfeelO76nbrgekXeR0FAAAAAMKOQhvFbn37d7o462b16snX\nCAAAACD20ISi1FPz56u6sVz3/+hir6MAAAAAgCcotFHqxlfv0JS+N3LtLAAAAICYRRuKQq9+8qkq\nmj7X4muf9zoKAAAAAHiGEdoo9PN/3Kn8xOuVOTDJ6ygAAAAA4BlGaKPM56UlWtv8tl754V+9jgIA\nAAAAnmKENspcN+t+ja69QqPzUr2OAgAAAACeYoQ2iuzYs1vzdj6qf5z/kddRAAAAAMBzjNBGkV88\n9YT6bD9F38of7nUUAAAAAPAchTZKNLtmPbn6Xk0/+jqvowAAAABARKDQRokH3nxLjbWJ+tVlk7yO\nAgAAAAARgUIbJf7w7gP6xoDpSkw0r6MAAAAAQERgUagosLqyQqU2T29dMcvrKAAAAAAQMRihjQK/\nfPpvyt11sUblpngdBQAAAAAiBiO0Ea7ZNev1DQ/rD5P+6XUUAAAAAIgojNBGuD/PeUuupr9+/K2j\nvI4CAAAAABGFQhvh/ue9h3T2wH9TfLzXSQAAAAAgsjDlOIIVbdygYr2j1777N6+jAAAAAEDEYYQ2\ngv3mub9r8M7zNHZYmtdRAAAAACDiUGgj2EvFT+qKI77rdQwAAAAAiEgU2ghVsHy5drmNuuk7+V5H\nAQAAAICIRKGNULe/+KQmNF+q1BRWgwIAAACA9rAoVARqds16f/tTevC0172OAgAAAAARixHaCPTI\nOwXSngG68huHeB0FAAAAACIWhTYC3fPukzol7buK49sBAAAAgP2iMkWYmvparWh+Qb867xKvowAA\nAABARKPQRpj7XpujntsP0+RjsryOAgAAAAARjUIbYR77+DlNGniR1zEAAAAAIOJRaCNITX2tVrlX\ndeO5F3odBQAAAAAiHoU2gtz32pvquf1QTToq0+soAAAAABDxKLQRhOnGAAAAANB5FNoIsae+Tqvc\nK/rlFKYbAwAAAEBnUGgjxMzX31TPHROVfzSrGwMAAABAZ1BoI8RjH/1Dp/ZnujEAAAAAdBaFNgI0\nNDWqsPlV/eTM87yOAgAAAABRg0IbAWZ/sEDxu7L19RNzvI4CAAAAAFGDQhsB/jrvJR3W85sy8zoJ\nAAAAAEQPCm0E+GjHS7riuG96HQMAAAAAogqF1mMfrlqlOrdLV59zpNdRAAAAACCqUGg9dvcbL2lY\n/TeVnMx8YwAAAADoipAXWjM728wKzWy1md3Yzvb+Zva6mX1mZl+Y2ZWhzhRJ3il7SReMZ7oxAAAA\nAHRVSAutmcVJminpLEkTJF1iZmPb7DZd0mfOucMlTZb032aWEMpckaJ0yyZtSfhC119wmtdRAAAA\nACDqhHqE9lhJa5xzJc65BkmzJbW92WqVpNTA41RJW5xzjSHOFRHuefV1DdhxurIGJnkdBQAAAACi\nTqhHQodIWt/q5zK1lNzWHpY018wqJKVI+k6IM0WMVwrf0KmDv+51DAAAAACISpEwtfc/JC11zk02\nsxGS3jKzQ51zu9ruOGPGjH2P8/PzlZ+fH7aQwdbU3KR17i3dc9rvvY4CAAAAAGFXUFCggoKCgzqG\nOeeCk6a9g5sdL2mGc+7swM83SXLOud+32uc1Sb9zzn0Q+HmupBudc4vbHMuFMmu4vfDRx/r2U1eq\n/k/LFcda0wAAAABinJnJOdel27+Eukp9LGmkmeWaWaKkqZJearPPSklfkyQzGyRptKSiEOfy3KPz\n5mh03FmUWQAAAADoppDWKedck1pWMX5T0nJJs51zK83sB2b2b4Hd7pB0tJktlfSWpF8656pDmSsS\nfLDhDZ034WyvYwAAAABA1ArplONg8tOU4407tmnQ77NV+pONys7s6XUcAAAAAPBcJE45RjsemDNX\nadtOoswCAAAAwEGg0Hrg+c/n6PgBTDcGAAAAgINBoQ0z55xW1L+hy086y+soAAAAABDVKLRhtmD1\najU2ShdPHut1FAAAAACIahTaMHu04B0NqT9NPXp06VpnAAAAAEAbFNowe/fLd3TK0NO8jgEAAAAA\nUY9CG0bNrlklVqArTp3sdRQAAAAAiHoU2jCa+8Uyudq+OuO4bK+jAAAAAEDUo9CG0WPz3lFu82mK\n41MHAAAAgINGtQqj+evf1eQ8phsDAAAAQDBQaMOkoalRZfHv6arJFFoAAAAACAYKbZi8uuRTxe8e\nqpMOG+R1FAAAAADwBQptmMz64F0Nj5ss4/azAAAAABAUFNow+aDyHZ0xgvvPAgAAAECwUGjDoKGp\nURt6fKgrJp/qdRQAAAAA8A0KbRi88dlSxe/K1tHj+3sdBQAAAAB8g0IbBs8seF85OoXrZwEAAAAg\niCi0YfBh2fs6cejJXscAAAAAAF+h0IaYc06lmq/vnEChBQAAAIBgotCG2Mfr1qm5sYfOOi7X6ygA\nAAAA4CsU2hB7av77GlR7snr04AJaAAAAAAgmCm2IvbvufR2dwXRjAAAAAAg2Cm2Iral7X+cfeYrX\nMQAAAADAdyi0IVS8aaNqE6p0cf4Er6MAAAAAgO9QaEPoyfc+UN+dJyo1Jd7rKAAAAADgOxTaEHpj\nxfs6JI3rZwEAAAAgFCi0IbRixwKdMfZEr2MAAAAAgC9RaEOkrrFe25KX6pJJR3sdBQAAAAB8iUIb\nIq8uWaoeO0ZqVG6K11EAAAAAwJcotCHy4uJFGmrHex0DAAAAAHyLQhsiC8sW6ujM47yOAQAAAAC+\nRaENkZKmRTr3CAotAAAAAIQKhTYE1m/ZrLqEjfrWKeO8jgIAAAAAvkWhDYFn3v9IqTuOUUpvPl4A\nAAAACBUaVwi8tWKRRvVkQSgAAAAACCUKbQgsrV6oU0dw/SwAAAAAhBKFNsiaXbM29fhIF51AoQUA\nAACAUKLQBtkHhaul2n46fmKG11EAAAAAwNcotEH2j4WLlNFwnOL4ZAEAAAAgpKhdQfZh8WId0u8Y\nr2MAAAAAgO9RaINs7e4lmjT6KK9jAAAAAIDvUWiDqKGpUduSl+qCE47wOgoAAAAA+B6FNogKlq9U\n/K6hGjc8zesoAAAAAOB7FNogennxEmU0HSUzr5MAAAAAgP9RaINoYckSje/H9bMAAAAAEA4U2iBa\ns3uJ8lkQCgAAAADCgkIbJA1Njdqe/Lm+xYJQAAAAABAWFNogeW95oeJ2D9b44X28jgIAAAAAMYFC\nGyQvL1migY0sCAUAAAAA4UKhDZIFxUs0gQWhAAAAACBsKLRBsnb3Ep06ikILAAAAAOFCoQ2CxqYm\nbUteqgtYEAoAAAAAwoZCGwTzVhbKajI1YURfr6MAAAAAQMyg0AbBq0s+1cCGI1kQCgAAAADCiEIb\nBB+VLNWeGGmFAAAgAElEQVSotMO8jgEAAAAAMYVCGwRrdi7VcbkUWgAAAAAIJwptEGyOX6qzD6fQ\nAgAAAEA4UWgPUtHGKjWpQacePtTrKAAAAAAQUyi0B+nlj5YqZddhSkxkRSgAAAAACCcK7UGat3qp\nspOYbgwAAAAA4UahPUhfbFqqwwZRaAEAAAAg3Ci0B6m8cakmj6PQAgAAAEC4UWgPwp76OtUkrdO5\nx4/3OgoAAAAAxBwK7UF4a+kKJewaoayByV5HAQAAAICYQ6E9CG9+vlQZzUw3BgAAAAAvUGgPwuL1\nSzW2L4UWAAAAALxAoT0I63Yt1fHDKLQAAAAA4AUKbTc551TdY6m+cSSFFgAAAAC8QKHtpsKKcjU3\nJej4QzK9jgIAAAAAMYlC201vfLJMqTUTFR/vdRIAAAAAiE0U2m5asHa5hiZO8DoGAAAAAMQsCm03\nLd+0XOMHUmgBAAAAwCsU2m4qq1+mE0Yc4nUMAAAAAIhZFNpuaHbN2pm8QmcdOd7rKAAAAAAQsyi0\n3fBZcYlU21cTRvT1OgoAAAAAxCwKbTe8+elypdUdIjOvkwAAAABA7KLQdsOiouXKTmJBKAAAAADw\nEoW2G1ZsWaZDMii0AAAAAOAlCm03VDQs10mjWOEYAAAAALxEoe2ixqYm7Uou1NePZoVjAAAAAPAS\nhbaLPl5bpLg9GRo+NMXrKAAAAAAQ0yi0XfTW0uXqW88KxwAAAADgNQptF31UvEw5PVkQCgAAAAC8\nRqHtosItyzVxEIUWAAAAALxGoe2iyqblOmUMKxwDAAAAgNcotF1Q39igmp5rdNZRY72OAgAAAAAx\nj0LbBR+tLVLc7izlZPXyOgoAAAAAxDwKbRcULCtUn4ZxXscAAAAAAIhC2yVLSgo1JInpxgAAAAAQ\nCSi0XbBqS6HGDaDQAgAAAEAkoNB2QUV9oY4eRqEFAAAAgEgQ8kJrZmebWaGZrTazG/ezT76ZfWpm\ny8zs3VBn6g7nnHYkFmryRAotAAAAAESChFAe3MziJM2UdLqkCkkfm9mLzrnCVvv0kXS/pDOdc+Vm\nNiCUmbrry00b5ZrjdMToiIwHAAAAADEn1CO0x0pa45wrcc41SJot6bw2+0yT9E/nXLkkOec2hzhT\nt7yztFC9do9VQoJ5HQUAAAAAoNAX2iGS1rf6uSzwXGujJaWb2btm9rGZfTfEmbpl4dpCZcQz3RgA\nAAAAIkVIpxx3UoKkIyWdJqm3pAVmtsA5t7btjjNmzNj3OD8/X/n5+WGKKC3bUKgRfSi0AAAAABAM\nBQUFKigoOKhjhLrQlkvKafXz0MBzrZVJ2uycq5VUa2bzJB0mqcNCG24luwt16ejTPHt/AAAAAPCT\ntoOUt912W5ePEeopxx9LGmlmuWaWKGmqpJfa7POipJPNLN7Mekk6TtLKEOfqsi1WqFPGMkILAAAA\nAJEipCO0zrkmM5su6U21lOdHnHMrzewHLZvdQ865QjObI+lzSU2SHnLOrQhlrq7aWVujhqRK5R8+\nzOsoAAAAAIAAc855naFTzMx5lfXVxUv1rVnTVH/3ck/eHwAAAAD8zszknOvSbWVCPeXYF+avLFS/\nZqYbAwAAAEAkodB2wqdlhcrtNc7rGAAAAACAVii0nbBuW6HGD2KEFgAAAAAiCYW2E6qaCnXCSAot\nAAAAAEQSCu0BNLtm7U5ercmHjvY6CgAAAACgFQrtAawoK5fVpWlUTprXUQAAAAAArVBoD2DesjXq\nVTtK1qXFowEAAAAAoUahPYAlxWuUET/K6xgAAAAAgDYotAewcuMa5aVRaAEAAAAg0lBoD2D9rjUa\nn0mhBQAAAIBIQ6E9gM1ujY4dQaEFAAAAgEhDoe1AU3OTant+qUkTR3odBQAAAADQBoW2A5+XrJft\nGaCcrF5eRwEAAAAAtEGh7cD85WuUUsctewAAAAAgElFoO7CkeI0yErh+FgAAAAAiEYW2A4Wb1mgY\nt+wBAAAAgIhEoe3A+t1rNHEwhRYAAAAAIhGFtgPVtkbHjqTQAgAAAEAkotDuR0NTo+qSSzTp0OFe\nRwEAAAAAtINCux+frCtRXE2msgYmex0FAAAAANAOCu1+vL9yjVLrmW4MAAAAAJGKQrsfn5Ss0aAe\nFFoAAAAAiFQU2v1YtWmNhvel0AIAAABApKLQ7kfZHm7ZAwAAAACRjEK7H1ttjY4bRaEFAAAAgEhF\noW1HfWOD6pPLdOrEYV5HAQAAAADsB4W2HZ8WlSquJksD0xO9jgIAAAAA2A8KbTsWri5SSv1wr2MA\nAAAAADpAoW3HZyXrNCBhhNcxAAAAAAAdoNC2Y83mIuWmMkILAAAAAJGMQtuOst1FGpNBoQUAAACA\nSEahbcfm5nU6Io8pxwAAAAAQySi0bTjnVJNUpBPGMkILAAAAAJGMQttGWfUWueY4TRjez+soAAAA\nAIAOUGjb+GBlkZJ2j1AcnwwAAAAARDRqWxufFBWpn5huDAAAAACRjkLbxvLKdRrck0ILAAAAAJGO\nQttG8fYijUin0AIAAABApKPQtlFVV6SJQ7llDwAAAABEuoTO7GRmfSSdIClPkpNUImmBc2576KJ5\nY0f8Oh0zkhFaAAAAAIh0HY7QmtnJZvaSpHmSLpGUK2lY4PF8M3vRzE4Ofczw2FNfp8akDTphfLbX\nUQAAAAAAB3CgEdoLJP27c25NexvNbLSkH0p6P9jBvLB4bYnid2erT2qnBq4BAAAAAB7qcITWOXe9\npHVmdvF+tq8O7OMLC1evU0oD040BAAAAIBoccFEo51yzpF+GIYvnPl9fpIweFFoAAAAAiAadXeX4\nbTO7wcyyzSx976+QJvPA2s1Fyk1jhWMAAAAAiAadvVj0O4Hff9zqOSfJV8OZZTXr9K3hJ3odAwAA\nAADQCZ0qtM65YaEOEgmqXZGOyPNVRwcAAAAA3zrQbXsmHegAZpYftDQecs5pT3KRThxPoQUAAACA\naHCgEdopZnaXpLmSFkuqlGSSMiUdI+lrkt6RVBDCjGFRvGmTXGOSRuf08ToKAAAAAKATOiy0zrlf\nmFmKpPMknSEpJ7CpVC33nv2dc25XaCOGx6LVxUreM0xmXicBAAAAAHRGh4XWzFrfY3atpDWSNkl6\n3zn3ZSiDhdvS4mL1cXlexwAAAAAAdNKBbtuT2upXSuD3oyW9bmZTQ5wtrAqripWZnOd1DAAAAABA\nJx1oyvFt7T0fuAft25JmhyKUF4q3FSunzzivYwAAAAAAOulAI7Ttcs5Vq2VxKN+oqi3W2EF5XscA\nAAAAAHRStwqtmU2WtDXIWTy1TcU6LC/P6xgAAAAAgE460KJQX0hybZ5Ol1Qh6fJQhQo355xqk0t0\n3Jhcr6MAAAAAADrpgPehbfOzk7TFObc7RHk8UbJ5s9SYrBFD07yOAgAAAADopAMtClUSriBeWrSq\nWEl78rgHLQAAAABEkW5dQ+s3n3EPWgAAAACIOhRatdyDdlAS188CAAAAQDSh0KrlHrS5ffK8jgEA\nAAAA6AIKrVruQTsmM8/rGAAAAACALqDQquUetIfn5nkdAwAAAADQBTFfaPfeg/bY0VxDCwAAAADR\nJOYLbenmLVJjokbl9PE6CgAAAACgC2K+0C5aVayk2lzuQQsAAAAAUSbmC+3SEu5BCwAAAADRKOYL\n7cqqYg1KzPM6BgAAAACgi2K+0BZv5R60AAAAABCNYr7QVtUWa8ygPK9jAAAAAAC6KOYL7TYV61Du\nQQsAAAAAUSemC61zTrVJ3IMWAAAAAKJRTBfasupqueZ4jc7p63UUAAAAAEAXxXShXVRYrMSaPMXF\n9KcAAAAAANEppqvc0tISpbkcr2MAAAAAALohpgvtmqr1GpjI9bMAAAAAEI1iutAWbyvVkNRsr2MA\nAAAAALohpgtt1Z5SjejPlGMAAAAAiEYxXWi3Nq3XuMEUWgAAAACIRjFdaHcnlOrwYUw5BgAAAIBo\nFLOFtrahXk1Jm3Xk6CyvowAAAAAAuiFmC+3nxeWKq8lSau8Er6MAAAAAALohZgvtp+tK1bOe6cYA\nAAAAEK1ittAuLytVvzgWhAIAAACAaBWzhXbd5vUalEyhBQAAAIBoFbOFdv2OUmX3YcoxAAAAAESr\nmC20m+rWa3QGI7QAAAAAEK1ittBuV6kOyabQAgAAAEC0itlCW5tUqiNHMuUYAAAAAKJVTBbajdu3\ny1mTxuT08zoKAAAAAKCbYrLQLlm7Xj1qcpSQYF5HAQAAAAB0U0wW2s+KS5XSxHRjAAAAAIhmIS+0\nZna2mRWa2Wozu7GD/Y4xswYzuyDUmVZXrlf/BBaEAgAAAIBoFtJCa2ZxkmZKOkvSBEmXmNnY/ex3\np6Q5ocyzV1F1qbJ6U2gBAAAAIJqFeoT2WElrnHMlzrkGSbMlndfOfj+R9A9JG0OcR5JUsXu9hqUz\n5RgAAAAAolmoC+0QSetb/VwWeG4fMxss6Xzn3AOSwrJK05bGUo3NYoQWAAAAAKJZgtcBJN0tqfW1\ntfsttTNmzNj3OD8/X/n5+d16w11xpTosj0ILAAAAAF4pKChQQUHBQR3DnHPBSdPewc2OlzTDOXd2\n4OebJDnn3O9b7VO096GkAZJ2S/o359xLbY7lgpG1qblZCbf2VNXPtmtQ/+SDPh4AAAAA4OCZmZxz\nXZq1G+oR2o8ljTSzXEmVkqZKuqT1Ds654Xsfm9nfJL3ctswG0+qKDbK6vpRZAAAAAIhyIS20zrkm\nM5su6U21XK/7iHNupZn9oGWze6jtS0KZR5KWrC1Vch3TjQEAAAAg2oX8Glrn3BuSxrR57i/72fd7\noc6zfP16pYkVjgEAAAAg2oV6leOIs3pDqQYmMUILAAAAANEu5gpt6fb1GprKCC0AAAAARLuYK7Qb\na8s0rP9Qr2MAAAAAAA5SzBXabc1lGpNFoQUAAACAaBdzhbYmvlwTcym0AAAAABDtYqrQNjY1qTG5\nSoePyPI6CgAAAADgIMVUoV1dsUFWm670PoleRwEAAAAAHKSYKrRLi8qVVMd0YwAAAADwg5gqtCvK\nypTqKLQAAAAA4AcxVWjXbixT/8QhXscAAAAAAARBTBXa9dvLNbg3I7QAAAAA4AcxVWirasqUm06h\nBQAAAAA/iKlCW91YplGDmHIMAAAAAH4QU4V2d1y5JmQzQgsAAAAAfhAzhdY5p/qeZTpiJCO0AAAA\nAOAHMVNoSzdXS43JGprR2+soAAAAAIAgiJlCu7SoXIl7hsrM6yQAAAAAgGCImUK7rLRMKc1cPwsA\nAAAAfhEzhXbNhjKlJ1BoAQAAAMAvYqbQlmwt06BeLAgFAAAAAH4RM4W2cne5cvoyQgsAAAAAfhEz\nhXZLfZlGZlBoAQAAAMAvYqbQ7rQyjRvClGMAAAAA8IuYKbR1SeU6fAQjtAAAAADgFzFRaDfv2Cln\nDRo1tK/XUQAAAAAAQRIThfazonIl7BmihATzOgoAAAAAIEhiotAuKy1Xr0amGwMAAACAn8REoV1V\nWaZ+8RRaAAAAAPCTmCi0xVvKlJHMCscAAAAA4CcxUWjLd5ZraB9GaAEAAADAT2Ki0G6uL9OIARRa\nAAAAAPCTmCi0212Zxg5hyjEAAAAA+ElMFNraxHIdmscILQAAAAD4ie8L7Z76ejUnbtXE4RleRwEA\nAAAABJHvC+2ykkrF1WQqOcn3pwoAAAAAMcX3LW95SYWSGwZ7HQMAAAAAEGS+L7SrKiuUKgotAAAA\nAPiN7wvtl5srlN6DQgsAAAAAfuP7Qlu2vUKZvSm0AAAAAOA3vi+0G/dUKLsvhRYAAAAA/Mb3hXZr\nQ6WGD6TQAgAAAIDf+L7Q7rIKjRmc5XUMAAAAAECQ+b7Q1iVW6JA8RmgBAAAAwG98XWh31OyRS6jR\nmOx0r6MAAAAAAILM14X2i+JKxddkqUcP8zoKAAAAACDIfF1oV6yvUHIj040BAAAAwI98XWhXV1Yo\nzSi0AAAAAOBHvi60xVsq1D+RQgsAAAAAfuTrQlu+o0KZvSm0AAAAAOBHvi60G/dUKKcfhRYAAAAA\n/MjXhXZrY4VGZFBoAQAAAMCPfF1od8dVaOwQCi0AAAAA+JGvC219UoUOzaPQAgAAAIAf+bbQbt6x\nU05NGjY4zesoAAAAAIAQ8G2h/aK4Ugm1WYqPN6+jAAAAAABCwLeFdsX6CvVsZLoxAAAAAPiVbwvt\nmqoKpRmFFgAAAAD8yreFtnhLhQYkUWgBAAAAwK98W2grdlYoK4VCCwAAAAB+5dtCu6m2UjnpFFoA\nAAAA8CvfFtptTRUamUGhBQAAAAC/8m2hrYmr0LihFFoAAAAA8CtfFlrnnOqTKzRxWJbXUQAAAAAA\nIeLLQlu1bYfUHK+cQaleRwEAAAAAhIgvC+3nX1aoR+1gmXmdBAAAAAAQKr4stCvXV6hXE9fPAgAA\nAICf+bLQrtlQoT5xFFoAAAAA8DNfFtqS6goNTKbQAgAAAICf+bLQVu6sUFYqKxwDAAAAgJ/5stBu\nqqtQbjojtAAAAADgZ74stNubKzR8ICO0AAAAAOBnviy0e+KrNHYIhRYAAAAA/Mx3hdY5p4akSh2S\nR6EFAAAAAD/zXaHduH2n5OKUPSjF6ygAAAAAgBDyXaFdVlKlhD1ZMvM6CQAAAAAglHxXaFeVVSq5\nKdPrGAAAAACAEPNdoV23oUppxvWzAAAAAOB3viu0JdWVSk9khBYAAAAA/M53hbZyZ5UyelNoAQAA\nAMDvfFdoN+2p1NA0phwDAAAAgN/5rtBua6xS3gBGaAEAAADA73xXaHepSqOyGKEFAAAAAL/zXaGt\nS6zUuGxGaAEAAADA73xVaOsaGtScuFXjcwd6HQUAAAAAEGK+KrSryjbKageoZ3K811EAAAAAACHm\nq0K7orRKSfVMNwYAAACAWOCrQru6qlK9HQtCAQAAAEAs8FWhLd5Upb4JjNACAAAAQCzwVaEt216l\nAcmM0AIAAABALPBVod2wq1JZqYzQAgAAAEAs8FWhra6vUk4/Ci0AAAAAxIKQF1ozO9vMCs1stZnd\n2M72aWa2NPDrfTOb2N332uEqNWIQU44BAAAAIBaEtNCaWZykmZLOkjRB0iVmNrbNbkWSTnXOHSbp\nt5Ie7u777Ymv0ujBjNACAAAAQCwI9QjtsZLWOOdKnHMNkmZLOq/1Ds65hc657YEfF0oa0p03cs6p\nIalKh+RRaAEAAAAgFoS60A6RtL7Vz2XquLBeLen17rxR5dYdUnO8hgxM6c7LAQAAAABRJsHrAHuZ\n2WRJV0k6eX/7zJgxY9/j/Px85efn7/t5eXGVEmqzZBa6jAAAAACA4CgoKFBBQcFBHcOcc8FJ097B\nzY6XNMM5d3bg55skOefc79vsd6ikf0o62zm3bj/Hch1lve/lAt3y9m+04555QcsPAAAAAAgPM5Nz\nrktDlKGecvyxpJFmlmtmiZKmSnqp9Q5mlqOWMvvd/ZXZzijaWKW0OK6fBQAAAIBYEdIpx865JjOb\nLulNtZTnR5xzK83sBy2b3UOSfi0pXdKfzcwkNTjnju3qe5VUVyo9kVv2AAAAAECsCPk1tM65NySN\nafPcX1o9vkbSNQf7PpU7qzSoNyO0AAAAABArQj3lOGw276nS0D6M0AIAAABArPBNod3WVKm8AYzQ\nAgAAAECs8E2h3aUqjcpihBYAAAAAYoVvCm19YqXGZzNCCwAAAACxwheFtra+Qc2J2zQ2Z4DXUQAA\nAAAAYeKLQrty/QbF1Q5UclK811EAAAAAAGHii0K7orRKifVMNwYAAACAWOKLQru2qkopYkEoAAAA\nAIglvii0X26uVN94RmgBAAAAIJb4otCWb6/SwJ6M0AIAAABALPFFod2wu1JZqYzQAgAAAEAs8UWh\nra6vUk46hRYAAAAAYokvCu2O5iqNGMSUYwAAAACIJb4otHsSKjVmMCO0AAAAABBLor7QOufUmFSl\nCXkUWgAAAACIJVFfaCu37pCaE5TVv7fXUQAAAAAAYRT1hXZl6QYl1A2SmddJAAAAAADhFPWFdnXF\nBiU3DvI6BgAAAAAgzKK+0H65cYNSjEILAAAAALEm6gttaXWV+vVgQSgAAAAAiDVRX2grd27QwJ6M\n0AIAAABArIn6QrupZoOy0ii0AAAAABBror7QbmvYoOx+FFoAAAAAiDVRX2h3ug0alkGhBQAAAIBY\nE/WFtjZhg0ZlUWgBAAAAINZEdaF1zqkxaYPG5VBoAQAAACDWRHWh3bh9l+RMQwameB0FAAAAABBm\nUV1oV5RsUHzdIJl5nQQAAAAAEG5RXWhXV2zQ/7V3PzF2nXcdh7+/mbGdEkdTK5N4koxrp3WKo2xC\nJSKqIlE2NHQTxAI1K2BRZdEKJDYtbMoSNpWKKlQBAQUEagEJmgWCgkqEWEAjQSCQv6Sq7TjxcaiT\n2I3tely/LO6NcAa7kYHrM++c55FGvvf4+Oo3izPv/fice+amSy43BgAAmKKug/abp05mb9bHHgMA\nAIARdB20x08P2bfLGVoAAIAp6jpoXz0zZO09ghYAAGCKug7aU+eH3HGLoAUAAJiiroP2jYtDDuwT\ntAAAAFPUddCebUMO3SZoAQAApqjroL2wPOSeOwUtAADAFHUdtJt7htx7QNACAABMUbdB+59nvpPU\n93Lg9lvGHgUAAIARdBu0zxwdsnxhPUtLNfYoAAAAjKDboH3+xJCbLrncGAAAYKq6DdpvnhqyN4IW\nAABgqroN2uOnh6zuErQAAABT1W3QvnJmyNpNghYAAGCqug3a184NueMWQQsAADBV3Qbt65tDNvYJ\nWgAAgKnqNmjPXh5y922CFgAAYKq6Ddrzy0PuuUPQAgAATFW3Qbu5e8i97xO0AAAAU9Vl0J4+ey5Z\nvpiD+1fHHgUAAICRdBm0zxwbsnx+f5aWauxRAAAAGEmXQfvCiSF7LrncGAAAYMq6DNqXhiF7I2gB\nAACmrMugPX56yOqKoAUAAJiyLoP2lTND1t4jaAEAAKasy6A9dW7I+l5BCwAAMGVdBu3rF4ds7BO0\nAAAAU9Zl0J69POTQmqAFAACYsi6D9vzSkA/euT72GAAAAIyoy6C9uOdk7j3gDC0AAMCUdRe0b751\nIVm+kEPr7x17FAAAAEbUXdA+c2zI0oXbs7xcY48CAADAiLoL2hdODNlzyeXGAAAAU9dd0L40DNnb\nBC0AAMDUdRe0x749ZHVF0AIAAExdd0F74s0hazcJWgAAgKnrLmhfOzdk/15BCwAAMHXdBe3pi0M2\n9glaAACAqesuaM9cPpm719bHHgMAAICRdRe055eGHL7DGVoAAICp6y5oN3cPOXJA0AIAAExdV0F7\n9tx303a9lQ/cuW/sUQAAABhZV0H77LFTWbpwW1aWuxobAACABeiqDJ8/MWTPpsuNAQAA6CxoXxqG\n3NwELQAAAJ0F7bFvD1ldEbQAAAB0FrSvvDnk1j2CFgAAgM6Cdjg3ZP2W9bHHAAAAYBvoKmhPf/dk\nNt7rDC0AAACdBe2Zy0MOrglaAAAAOgvaczXk8LqgBQAAoLOg3dw95MgBQQsAAEBnQdt2n8k9d906\n9hgAAABsA10F7dKFtexa6WpkAAAAFqSrOty96XJjAAAAZroK2puboAUAAGCmq6BdXRG0AAAAzHQV\ntLfuWR97BAAAALaJroJ2/15naAEAAJjpKmjvWhW0AAAAzHQVtIduE7QAAADMdBW0h9cFLQAAADML\nD9qqerCqnquqF6rqM9fY5zeq6sWqeqqq7r/Wax3ZELQAAADMLDRoq2opyReTfCzJfUkerqojW/b5\nySQfaK3dk+SRJF+61ut9cGNtgdMCV/PEE0+MPQJMluMPxuHYg34s+gztA0lebK0dba1tJvlykoe2\n7PNQkt9PktbaPyZZraqrnordvWt5kbMCV2FRh/E4/mAcjj3ox6KD9q4kx694/vJ82/fb58RV9gEA\nAIB36OqmUAAAAPC2aq0t7sWrfiTJr7bWHpw//2yS1lr79Sv2+VKSv22tfWX+/LkkP9ZaG7a81uIG\nBQAAYHSttbqe/VcWNcjck0kOV9XBJK8m+USSh7fs83iSTyX5yjyA39gas8n1f2MAAADsbAsN2tba\n96rq00m+ltnlzY+21p6tqkdmf91+q7X2F1X18ar6jyRvJfn5Rc4EAADAzrDQS44BAABgUbq4KVRV\nPVhVz1XVC1X1mbHngamoqm9V1b9U1T9X1TfGngd2qqp6tKqGqvrXK7btq6qvVdXzVfVXVbU65oyw\nU13j+PtcVb1cVf80/3pwzBlhp6mqjar6elX9e1U9XVW/MN9+3Wvftg/aqlpK8sUkH0tyX5KHq+rI\nuFPBZFxO8tHW2g+11h4YexjYwX4vs3XuSp9N8jettR9M8vUkv3zDp4JpuNrxlySfb619aP71lzd6\nKNjhLiX5pdbafUk+nORT88a77rVv2wdtkgeSvNhaO9pa20zy5SQPjTwTTEWlj58T0LXW2t8neX3L\n5oeSPDZ//FiSn7qhQ8FEXOP4S2ZrILAArbWTrbWn5o+/k+TZJBv5X6x9PbxRvSvJ8SuevzzfBixe\nS/LXVfVkVX1y7GFgYm5/+67/rbWTSW4feR6Ymk9X1VNV9Tsu+YfFqapDSe5P8g9J9l/v2tdD0ALj\n+Uhr7UNJPp7ZpSA/OvZAMGHu4gg3zm8meX9r7f4kJ5N8fuR5YEeqqr1J/jTJL87P1G5d69517esh\naE8ked8Vzzfm24AFa629Ov/ztSR/ltlHAIAbY6iq/UlSVetJTo08D0xGa+219t+/CuS3k/zwmPPA\nTlRVK5nF7B+01r4633zda18PQftkksNVdbCqdif5RJLHR54Jdryq+oH5/5qlqm5O8hNJ/m3cqWBH\nq7zzM3uPJ/m5+eOfTfLVrf8A+H/zjuNv/kb6bT8d6x8swu8meaa19oUrtl332tfF76Gd3yr9C5kF\n+KOttV8beSTY8arq7szOyrYkK0n+0LEHi1FVf5Tko0luTTIk+VySP0/yJ0kOJDma5Gdaa2+MNSPs\nVOjL7fMAAAEKSURBVNc4/n48s8/0XU7yrSSPvP25PuD/rqo+kuTvkjyd2XvNluRXknwjyR/nOta+\nLoIWAAAAturhkmMAAAD4HwQtAAAAXRK0AAAAdEnQAgAA0CVBCwAAQJcELQAAAF0StAAAAHRJ0ALA\nNldVNfYMALAdCVoA2Gaq6mBVPVdVj1XV00k2xp4JALajaq2NPQMAcIWqOpjkpSQfbq09OfY8ALBd\nOUMLANvTUTELAN+foAWA7emtsQcAgO1O0ALA9uRGUADwLgQtAGxPbnIBAO/CTaEAAADokjO0AAAA\ndEnQAgAA0CVBCwAAQJcELQAAAF0StAAAAHRJ0AIAANAlQQsAAECXBC0AAABd+i+7RhbW9sd4DgAA\nAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x109da3320>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.title('Hartree potential (Hydrogen)')\n",
+    "plt.plot(rs, U, label='verlet')\n",
+    "plt.plot(rs, U_exact, label='exact')\n",
+    "plt.xlabel('r')\n",
+    "plt.ylabel('U(r)')\n",
+    "plt.legend(loc='best')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# (c) Helium"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "DFT for two electrons with given nuclear potential."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "00: eps = -1.993684, #iterations = 41\n",
+      " eps_old is none \n",
+      "01: eps = -0.393066, #iterations = 12\n",
+      ", |eps - eps_old| = 1.600618\n",
+      "02: eps = -0.672200, #iterations = 22\n",
+      ", |eps - eps_old| = 0.279133\n",
+      "03: eps = -0.545101, #iterations = 19\n",
+      ", |eps - eps_old| = 0.127099\n",
+      "04: eps = -0.596247, #iterations = 20\n",
+      ", |eps - eps_old| = 0.051146\n",
+      "05: eps = -0.574245, #iterations = 19\n",
+      ", |eps - eps_old| = 0.022001\n",
+      "06: eps = -0.583467, #iterations = 20\n",
+      ", |eps - eps_old| = 0.009222\n",
+      "07: eps = -0.579557, #iterations = 20\n",
+      ", |eps - eps_old| = 0.003910\n",
+      "08: eps = -0.581207, #iterations = 16\n",
+      ", |eps - eps_old| = 0.001650\n",
+      "09: eps = -0.580502, #iterations = 19\n",
+      ", |eps - eps_old| = 0.000706\n",
+      "E =  -2.83096625252\n",
+      "CPU times: user 30 s, sys: 53.5 ms, total: 30.1 s\n",
+      "Wall time: 30.1 s\n"
+     ]
+    }
+   ],
+   "source": [
+    "%%time\n",
+    "def two_electron_dft(V_nucl, maxiter=100, stoptol=0.001, verbose=False):\n",
+    "    # initial estimate for u(r) corresponding to V_h = V_xc = 0\n",
+    "    V_hartree = np.zeros_like(V_nucl)\n",
+    "    V_xc = np.zeros_like(V_nucl)\n",
+    "    eps = None\n",
+    "    \n",
+    "    for i in range(maxiter):\n",
+    "        # compute new u(r)\n",
+    "        eps_old = eps\n",
+    "        eps, u, numiter = solve_schrodinger(V_nucl + V_hartree + V_xc, retnumiter=True)\n",
+    "        u /= np.linalg.norm(u) * sqrt(dr)\n",
+    "        if verbose:\n",
+    "            print(\"%02d: eps = %f, #iterations = %d\" % (i, eps, numiter),)\n",
+    "\n",
+    "        # converged?        \n",
+    "        if eps_old is not None:\n",
+    "            if verbose:\n",
+    "                print(\", |eps - eps_old| = %f\" % abs(eps - eps_old))\n",
+    "            if abs(eps - eps_old) < stoptol:\n",
+    "                return 2*eps - dr * np.dot(V_hartree, u**2) - 0.5 * dr * np.dot(V_xc, u**2)\n",
+    "        elif verbose:\n",
+    "            print(\" eps_old is none \")\n",
+    "\n",
+    "        # update Hartree potential by integration the electron density\n",
+    "        U = solve_poisson(u)\n",
+    "        V_hartree = 2. * U / rs\n",
+    "        \n",
+    "        # update exchange-correlation potential\n",
+    "        #V_xc = -(3. / pi * 2. * u**2 / (rs**2 * 4 * pi)) ** (1./3)\n",
+    "        \n",
+    "        # exchange potential as in equation 8.47 in the notes:\n",
+    "        rhos = (1./3.*2.*u**2/rs**2)**(-1./3.)\n",
+    "        V_xc = -(3./(2.*pi))**(2./3.)*1./rhos*(1.+0.0545*rhos*np.log(1.+11.4/rhos))\n",
+    "\n",
+    "    raise Exception(\"Not converged after %d iterations; |eps - eps_old| = %f\" % (maxiter, abs(eps - eps_old)))\n",
+    "\n",
+    "print(\"E = \", two_electron_dft(-2./rs, verbose=True))\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We should obtain $ \\epsilon \\approx 0.52 \\text{ a.u.}$ and $ E \\approx 2.72 \\text{ a.u.} $  for the exchange potenial used in exercise sheet 10.\n",
+    "\n",
+    "The exact value is $ E = -2.903 \\text{ a.u.}$  (Computational Physics, J. Thijssen, chapter 4)\n",
+    "\n",
+    "Remember from the previous exercise that the Hartree ground state energy is -2.85516035589. \n",
+    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.5.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}