Hatom.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c43dd52e",
   "metadata": {},
   "source": [
    "$$ \n",
    "\\begin{aligned}\n",
    "&\\frac{1}{r^2} \\frac{d}{dr} r^2 \\frac{d}{dr} f(r) \\\\\n",
    "= &\\frac{1}{r^2} \\frac{d}{dr} r^2 \\frac{df}{dr} \\\\\n",
    "= &\\frac{1}{r^2} \\left( 2 r  \\frac{df}{dr} +    r^2 \\frac{d^2f}{dr^2}    \\right) \\\\\n",
    "= &\\frac{2}{r}  \\frac{df}{dr} + \\frac{d^2f}{dr^2}\n",
    "\\end{aligned}\n",
    "$$\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\frac{1}{r} \\frac{d^2}{dr^2} r f(r) \\\\\n",
    "= &\\frac{1}{r} \\frac{d}{dr} \\left( f +  r \\frac{df}{dr}    \\right) \\\\\n",
    "= &\\frac{1}{r} \\left( 2 \\frac{df}{dr} +  r \\frac{d^2f}{dr^2}   \\right) \\\\\n",
    "= &\\frac{2}{r} \\frac{df}{dr} + \\frac{d^2f}{dr^2}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6216e13f",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\frac{1}{2} \\frac{1}{r} \\frac{d^2}{dr^2} r \\psi(r) + \\frac{1}{2 r^2} l (l+1) \\psi(r) - \\frac{Z}{r} \\psi(r) = E \\psi(r)\n",
    "$$\n",
    "$$\n",
    " \\phi(r) = r \\psi(r)\n",
    " $$\n",
    " $$\n",
    "-\\frac{1}{2} \\frac{d^2}{dr^2} \\phi(r) + \\frac{1}{2 r^2} l (l+1) \\phi(r) - \\frac{Z}{r} \\phi(r) = E \\phi(r)\n",
    "$$\n",
    "$$\n",
    "-\\frac{1}{2} \\frac{d^2}{dr^2} \\phi(r) - \\frac{Z}{r} \\phi(r) = E \\phi(r)\n",
    "$$\n",
    "$$\n",
    " r \\in [0,\\infty]\n",
    "$$\n",
    "$$\n",
    "  \\psi_{0} (r) = \\exp(-Z r)\n",
    "  $$\n",
    "$$\n",
    " \\psi(r) = 10^{-10} \\implies r > \\log(10^{-10}) / Z \\approx 2.3 \\times 10 / Z \n",
    " $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9c82b644",
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting up everything we need for the matrices\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "from scipy.special import roots_legendre, eval_legendre\n",
    "import math\n",
    "\n",
    "Z = 1\n",
    "qn_l = 0 # quantum number l\n",
    "L = math.log(10) * 10 / Z\n",
    "\n",
    "n = 50\n",
    "epsilon = 1e-6\n",
    "roots, weights = roots_legendre(n)\n",
    "P = np.array([eval_legendre(i, roots) for i in range(n+1)])\n",
    "dP = np.zeros((n+1, n))\n",
    "dP[1] = np.ones(n)\n",
    "for i in range(2,n+1):\n",
    "    dP[i] = (2*i-1) * P[i-1] + dP[i-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "47763d77",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[23.02585093  0.          0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          7.67528364  0.         ...  0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          4.60517019 ...  0.          0.\n",
      "   0.        ]\n",
      " ...\n",
      " [ 0.          0.          0.         ...  0.24237738  0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.23737991\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         ...  0.          0.\n",
      "   0.23258435]]\n"
     ]
    }
   ],
   "source": [
    "overlap_matrix = np.zeros((n,n))\n",
    "for l in range (n):\n",
    "    for m in range (n):\n",
    "        approx_overlap = np.dot(P[l]*P[m], weights) * L * 0.5\n",
    "        if abs(approx_overlap) < epsilon:\n",
    "            approx_overlap = 0\n",
    "        overlap_matrix[l, m] = approx_overlap\n",
    "print(overlap_matrix)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b816251",
   "metadata": {},
   "source": [
    "PE: $$\\int_{-1}^{1} P_l \\frac{Z}{\\frac{L}{2}(y + 1)} P_m \\frac{L}{2} dy$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c908961a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-8.99841068  6.99841068 -5.99841068 ...  0.12248299 -0.08081633\n",
      "   0.04      ]\n",
      " [ 6.99841068 -6.99841068  5.99841068 ... -0.12248299  0.08081633\n",
      "  -0.04      ]\n",
      " [-5.99841068  5.99841068 -5.99841068 ...  0.12248299 -0.08081633\n",
      "   0.04      ]\n",
      " ...\n",
      " [ 0.12248299 -0.12248299  0.12248299 ... -0.12248299  0.08081633\n",
      "  -0.04      ]\n",
      " [-0.08081633  0.08081633 -0.08081633 ...  0.08081633 -0.08081633\n",
      "   0.04      ]\n",
      " [ 0.04       -0.04        0.04       ... -0.04        0.04\n",
      "  -0.04      ]]\n"
     ]
    }
   ],
   "source": [
    "pe_matrix = np.zeros((n, n))\n",
    "for l in range (n):\n",
    "    for m in range (n):\n",
    "        approx_pe = np.dot(P[l] * P[m] / (roots + 1), weights) * -Z\n",
    "        if abs(approx_pe) < epsilon:\n",
    "            approx_pe = 0\n",
    "        pe_matrix[l, m] = approx_pe\n",
    "print(pe_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "96accbf8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00\n",
      "  0.00000000e+00 0.00000000e+00]\n",
      " [0.00000000e+00 8.68588964e-02 0.00000000e+00 ... 8.68588964e-02\n",
      "  0.00000000e+00 8.68588964e-02]\n",
      " [0.00000000e+00 0.00000000e+00 2.60576689e-01 ... 0.00000000e+00\n",
      "  2.60576689e-01 0.00000000e+00]\n",
      " ...\n",
      " [0.00000000e+00 8.68588964e-02 0.00000000e+00 ... 9.79768351e+01\n",
      "  0.00000000e+00 9.79768351e+01]\n",
      " [0.00000000e+00 0.00000000e+00 2.60576689e-01 ... 0.00000000e+00\n",
      "  1.02146062e+02 0.00000000e+00]\n",
      " [0.00000000e+00 8.68588964e-02 0.00000000e+00 ... 9.79768351e+01\n",
      "  0.00000000e+00 1.06402148e+02]]\n"
     ]
    }
   ],
   "source": [
    "ke_matrix = np.zeros((n, n))\n",
    "for l in range (n):\n",
    "    for m in range (n):\n",
    "        approx_ke = np.dot(dP[l] * dP[m], weights) / L\n",
    "        if abs(approx_ke) < epsilon:\n",
    "            approx_ke = 0\n",
    "        ke_matrix[l, m] = approx_ke\n",
    "print(ke_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5dc93f9a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "angular_ke_matrix = np.zeros((n, n))\n",
    "for l in range (n):\n",
    "    for m in range (n):\n",
    "        approx = np.dot(P[l] * P[m] / (1 + roots)**2, weights) * qn_l * (qn_l + 1) / L\n",
    "        if abs(approx) < epsilon:\n",
    "            approx = 0\n",
    "        angular_ke_matrix[l, m] = approx\n",
    "print(angular_ke_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bdf820ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.34294481e-02 -1.07654941e-01  1.18724566e-01 -7.83959027e-02\n",
      "  1.40633078e-02  4.26417611e-02 -7.35007348e-02  7.76013945e-02\n",
      " -6.45539199e-02  4.55903067e-02 -2.83051658e-02  1.57664128e-02\n",
      " -7.98647935e-03  3.71540711e-03 -1.59949154e-03  6.41112296e-04\n",
      " -2.40471634e-04  8.47696630e-05 -2.81890099e-05  8.87152698e-06\n",
      " -2.65003480e-06  7.53290279e-07 -2.04240932e-07  5.29310855e-08\n",
      " -1.31372245e-08  3.12815551e-09 -7.15764019e-10  1.57613336e-10\n",
      " -3.34441802e-11  6.84378063e-12 -1.34891477e-12  2.53616855e-13\n",
      " -4.35231576e-14  5.66588983e-15 -2.98743443e-16  8.55171753e-16\n",
      " -2.22500604e-15  3.38975605e-15 -3.92760978e-15  3.64052199e-15\n",
      " -2.66362646e-15  1.01747111e-15  1.03481114e-15 -2.96192342e-15\n",
      "  4.48003041e-15 -5.10866900e-15  4.67789309e-15 -3.17399505e-15\n",
      "  9.32396113e-16  1.51111486e-15]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.014872513059998168, 0.014872513059998156)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#obtaining coefficients\n",
    "phi = lambda r: r * (math.e)**(-Z * r)\n",
    "c = np.zeros(n)\n",
    "for l in range(n):\n",
    "    c[l] = ((2*l+1)/2) * np.dot(phi(L * 0.5 * (roots + 1)) * P[l], weights)\n",
    "print(c)\n",
    "fit = lambda x : np.dot(eval_legendre(range(n), 2 * x/L -1), c)\n",
    "fit(6),phi(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3cf0eb43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2334eb86c50>,\n",
       " <matplotlib.lines.Line2D at 0x2334eb86c80>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQd0lEQVR4nO3deXhU5eH28e9MtpmQBUIgCwQS9i0QZAlR3Gok4IpbEW3B1OqvuFSbqhUroNU2StUXrVRaW1uXqtRaqVWL0khwC4sgsgiUJRC2BAhmZSbLzHn/GDIYTSCT7UyS+3NdcyUzeXLOPQlxbp855zkWwzAMRERERPyY1ewAIiIiImeiwiIiIiJ+T4VFRERE/J4Ki4iIiPg9FRYRERHxeyosIiIi4vdUWERERMTvqbCIiIiI3ws0O0BrcLvdHDp0iPDwcCwWi9lxREREpAkMw6C8vJz4+His1tPPoXSKwnLo0CESEhLMjiEiIiLNsH//fvr27XvaMZ2isISHhwOeJxwREWFyGhEREWmKsrIyEhISvK/jp9MpCkvd20AREREqLCIiIh1MUw7n0EG3IiIi4vdUWERERMTvqbCIiIiI31NhEREREb+nwiIiIiJ+T4VFRERE/J4Ki4iIiPg9FRYRERHxeyosIiIi4vdUWERERMTvqbCIiIiI31NhEREREb+nwtJO3v/15/z9Z5+ZHUNERKRD6hRXa/Z3NSdqmPrgeAB6DdjChXeOMjmRiIhIx6IZlnawf+1h7+e/fqjaxCQiIiIdkwpLO8hfd8z7+fzjd8Mnn5gXRkREpANSYWkH+VsqAZjGe5zHx/DwwyYnEhER6VhUWNrBnp0uAJL6uyEwEP77X4yPNcsiIiLSVCos7cAoKSWSEpJGhXH8hju4n2yuutpidiwREZEOQ4WlHWRHZFNCD342q5jyW3/OU2Txr2Pn8PGzX5odTUREpENQYWkP+fkABAxKov85ffnR8DwAHprvMjOViIhIh6HC0tYqK+HIEc/nSUkAzH1+IEFU8+HXZ/HR7zTLIiIiciYqLG1s+4eHGMFWZga9AT16AJycZVkNwMMLNMsiIiJyJiosbWzXuq/Zxgi2B9Zf3faBP2uWRUREpKlUWNpY/lcOAJJ6lNR7vF9aH24e4Zll+fUCrX4rIiJyOs0qLIsXLyYxMRGbzUZqaipr165tdOw///lPxo8fT/fu3enWrRspKSm8/PLL9cbcdNNNWCyWerepU6c2J5rfyc83AEiKq/rO1+b+aSC3WZ/jT19fDatWtXc0ERGRDsPnwrJ06VKysrJYsGABGzZsYMyYMWRkZHCk7sDSb4mKiuKXv/wleXl5bNq0iczMTDIzM3n//ffrjZs6dSqHDx/23l577bXmPSM/k3/YBkDSgO+uu9IvrQ+Lb91EAge0+q2IiMhp+FxYnnrqKW655RYyMzMZMWIES5YsITQ0lBdeeKHB8RdccAFXXXUVw4cPZ+DAgdx1112MHj2aT751PZ2QkBBiY2O9tx4nD1Dt6PK/9jyPpBH2hgc88AAEBcHKlTg/+Kgdk4mIiHQcPhWW6upq1q9fT3p6+qkNWK2kp6eTl5d3xu83DIOcnBx27NjBeeedV+9rubm59O7dm6FDhzJnzhyKi4sb3U5VVRVlZWX1bv7IcBvkO2MBSBoX1fCghAT2fP8XTOctLpkR1o7pREREOg6fCsuxY8dwuVzExMTUezwmJobCwsJGv6+0tJSwsDCCg4O59NJL+d3vfsfFF1/s/frUqVN56aWXyMnJ4fHHH2fVqlVMmzYNl6vhU36zs7OJjIz03hISEnx5Gu3GcfA4SeQTSQmJZ8c3Oi7ozjm8xyWsLDmL3EUb2y+giIhIB9EuZwmFh4ezceNG1q1bx69//WuysrLIzc31fv3666/niiuuIDk5menTp/POO++wbt26emO+ae7cuZSWlnpv+/fvb4+n4bPQonw2MpaS2OGE9mzkLSEgITWeH488uS7Lw0Z7xRMREekwfCos0dHRBAQEUFRUVO/xoqIiYmNjG9+J1cqgQYNISUnh5z//Oddeey3Z2dmNjh8wYADR0dHs2rWrwa+HhIQQERFR7+aXTi7JX7fC7enM/fMgAqkht2Qs+R/5ZwETERExi0+FJTg4mHHjxpGTk+N9zO12k5OTQ1paWpO343a7qar67mm+dQ4cOEBxcTFxcXG+xPM/PhSWhNR4ku2egvblfw61ZSoREZEOx+e3hLKysnj++ed58cUX2bZtG3PmzKGyspLMzEwAZs2axdy5c73js7OzWbFiBXv27GHbtm08+eSTvPzyy/zgBz8AoKKignvvvZfVq1ezd+9ecnJyuPLKKxk0aBAZGRmt9DTN8fO/ncVItvByxVVNGj863nOg8aY1jraMJSIi0uEE+voNM2bM4OjRo8yfP5/CwkJSUlJYvny590DcgoICrNZTPaiyspLbbruNAwcOYLfbGTZsGK+88gozZswAICAggE2bNvHiiy9SUlJCfHw8U6ZM4ZFHHiEkJKSVnqY5th7szleMpCbqeJPGJw+vhd2w+X8d+3mLiIi0NothGB3+KM+ysjIiIyMpLS31q+NZhgbn87+aJD588gsuzBp7xvGrn1nLs3f9jwtjtnFz4a/bIaGIiIh5fHn99nmGRZrGXetmb43nVOakib2a9D2TZvRn0l2pcMQCJ34JoaFtGVFERKTD0MUP28ihL4qoJoQAauk7vvEzqOqJiYHevcEwYOvWtg0oIiLSgaiwtJH8NZ5rK/ULPESgrekTWa5RY9jBEPauzG+raCIiIh2OCksbyd9UDkBS+DGfvu8XpXMZxg4WvRZz5sEiIiJdhApLGwn5upAxbCQ5vvFrIjVkVEoQAJv2+s/BwyIiImZTYWkjM7q9w0bGsmjmWp++L/nCngBsKu2H4e7wJ3CJiIi0ChWWtuLDKrffNGJaIlZcFBs9Kdx0pA2CiYiIdDwqLG3E2HOysAwY4NP32aPsDAneC8Cmd3VNIREREVBhaRPVFdVEHtjCKDZT2tO3wgKQ3Nszs7Lps4rWjiYiItIhqbC0gYI1hykngj0MIGJg0xaN+6bRQz0Xhty8Xev6iYiIgFa6bRP5644B/UkMPozFOtDn758y1UptzkOcx0FgcqvnExER6Wg0w9IG8rdUApAU2bSLHn7bxGv78RAP8739L0JNTWtGExER6ZBUWNpA/m43AEmxjuZtoH9/CA/3lJUdO1oxmYiISMekwtIG8g8GAzAgyd28DVgsHB56Ae9yCVv/s68Vk4mIiHRMKixtIP+4Z5XapGG2Zm/jVyfu4TLe5ZV/dmutWCIiIh2WDrptAyNqN+PAzaDx3Zu9jdGjga9g024VFhERERWW1lZRwV+qbvB8PqWk2ZtJPrc7vA6bj8e3SiwREZGOTG8Jtba6JfmjoiAystmbSb6sPwD7XX34Or+kFYKJiIh0XCosrax2114M8PkaQt8W2S+S/gEHANj8jg68FRGRrk2FpZUtfjmCCMrIKnuoxdtKjj4EwOaPS1q8LRERkY5MhaWV7cmHCsIJ7B7W4m2NHngCgE2bLS3eloiISEemg25bWf5hOwBJA1veBa+7oorhn/2AiSHVwHkt3p6IiEhHpcLSyvJLegCQNDK0xdtKuaIfKff/DXZ1A7cbrJoQExGRrkmvgK3IcBvkV8UBkDQuquUbHDwYQkKgsvLU2UciIiJdkApLKzq2o5hKPMeu9E9rhfVTAgP5PPEanuanrF9W0PLtiYiIdFAqLK0of3URAPHWw9i6N39Z/m9aXPsT7uZp3nlHB96KiEjXpWNYWlFQ0QEuZzdR0YFAXKtsM3mkC3bDpv+1TgESERHpiDTD0orG8gVvcyV/zXi91bY5+hzPhRQ3H41ptW2KiIh0NCosranuwNgWrnL7TaMv6QvArpr+VB490WrbFRER6UhUWFpR+c7CVlmW/5t6j+pNjPUIBla2vru31bYrIiLSkaiwtKKzPnmaCMpYVzW6Vbeb3H0/AJtzi1t1uyIiIh2FCksrcVW72FcTTwXhxIxu3eNNRieVA7Bpo7tVtysiItJR6CyhVnLoiyJqiCeIavqMi23Vbf/ftce5bv0kRkVEAee36rZFREQ6As2wtJL8tUcB6Bd4mIDggFbd9pApiUxiDWFfrQXDaNVti4iIdATNKiyLFy8mMTERm81Gamoqa9eubXTsP//5T8aPH0/37t3p1q0bKSkpvPzyy/XGGIbB/PnziYuLw263k56ezs6dO5sTzTR7vvS8bZMUcaz1Nz5ihOc6QsXFcPhw629fRETEz/lcWJYuXUpWVhYLFixgw4YNjBkzhoyMDI4cOdLg+KioKH75y1+Sl5fHpk2byMzMJDMzk/fff987ZuHChTzzzDMsWbKENWvW0K1bNzIyMnA6nc1/Zu0sf2ctAEm9K1t/4zYbb8bezp08w5q/72v97YuIiPg5nwvLU089xS233EJmZiYjRoxgyZIlhIaG8sILLzQ4/oILLuCqq65i+PDhDBw4kLvuuovRo0fzySefAJ7ZlUWLFvHggw9y5ZVXMnr0aF566SUOHTrEsmXLWvTk2lP+Ac/hQEn92+bA2DcCrudZ7mTViuo22b6IiIg/86mwVFdXs379etLT009twGolPT2dvLy8M36/YRjk5OSwY8cOzjvvPADy8/MpLCyst83IyEhSU1ObtE1/MZ7PuZy3GTu2ba75kzzUU1Q2faXjpEVEpOvx6dXv2LFjuFwuYmLqn7YbExPD9u3bG/2+0tJS+vTpQ1VVFQEBAfz+97/n4osvBqCwsNC7jW9vs+5r31ZVVUVVVZX3fllZmS9Po038tOoJfspBuGpNm2x/9KRQ+C9sLoxuk+2LiIj4s3Y5Syg8PJyNGzeybt06fv3rX5OVlUVubm6zt5ednU1kZKT3lpCQ0Hphm6OqCg4d8nzeiqvcftPoDM/FFLc5k6iurGmTfYiIiPgrnwpLdHQ0AQEBFBUV1Xu8qKiI2NjG1x6xWq0MGjSIlJQUfv7zn3PttdeSnZ0N4P0+X7Y5d+5cSktLvbf9+/f78jRanfN/BRQbPTBCu0F028yA9Du7LxGUUkMwO97f2yb7EBER8Vc+FZbg4GDGjRtHTk6O9zG3201OTg5paWlN3o7b7fa+pZOUlERsbGy9bZaVlbFmzZpGtxkSEkJERES9m5ly3y4jmmLONj4FS9scw2KxWkgO95whtDmn4TOyREREOiufj+DMyspi9uzZjB8/nokTJ7Jo0SIqKyvJzMwEYNasWfTp08c7g5Kdnc348eMZOHAgVVVVvPfee7z88ss899xzAFgsFu6++24effRRBg8eTFJSEvPmzSM+Pp7p06e33jNtQ/lbPKcy94po29OwR/cr4dOtsGuTrtosIiJdi8+FZcaMGRw9epT58+dTWFhISkoKy5cv9x40W1BQgNV6auKmsrKS2267jQMHDmC32xk2bBivvPIKM2bM8I657777qKys5NZbb6WkpITJkyezfPlybDZbKzzFtpe/x3Mqc1Js2xaWB2/M5+EHrqFXxETg4jbdl4iIiD+xGEbHX+u9rKyMyMhISktLTXl76LqEPP5xII3/N30Vd7/Vhtf6+fRTmDwZ+vYFk4/bERERaSlfXr91LaFWkH88EoCk4W08IzRqlOfjgQNw/Hjb7ktERMSPqLC0gnyH52ymAeN6tO2OIiPJjnyMS3iXz9/SDIuIiHQdKiwtVHOihuNGFAB9Unq1+f5yAy/iP1zCF6vMXyxPRESkvWid9xZyn3DyI17CiY1u0de0+f6GJ1TyQTFs29o21ywSERHxRyosLRTidvBnfuy5E3FDm+9v+AgLbITt+0PbfF8iIiL+Qm8JtZTD4flos7XZonHfNHyS5wDfbV83vrKwiIhIZ6PC0kK1FU7KCKfGFt4u+xue3geAfbV9OHFMC8iJiEjXoMLSQlu2QCRl9C/b1C776zWsJz0txRhY2bGioF32KSIiYjYVlhZylFYDYLNWt88OLRaGhx8gglKKvixsn32KiIiYTIWlhZwVtQDYA9qpsADvTX+eErozNfjDdtuniIiImXSWUAs5ytu/sISPGeD5ZNu2dtuniIiImTTD0kKOkzMstsDa9tvp8OGejyosIiLSRaiwtJCz0rOAm70dC0tpnxFM5y1GbV1KrbMdi5KIiIhJVFhayFFXWIJq2m2f4SMS+IApbGUke1bpmkIiItL5qbC0UGJ4Md9nKZNjdrXbPq2BVobaPac0b//oSLvtV0RExCwqLC2U3n8nS7me+8auaNf9Du9dDMC2DY523a+IiIgZVFha6ptL87ej4QM9b0Ft26kTvUREpPNTYWmhmooq3FjAbm/X/Q5PCQFgW2H3dt2viIiIGVRYWujeFVMIwM2Da69o1/0OP68XANsq+2G4jXbdt4iISHtTYWkhZ5XnCs3BIW1/peZvGnRhAhGUMpidlO043K77FhERaW8qLC3kqPL8CNv5HSFCIkIoGTSB9Ywn8pAWkBMRkc5NhaWFHNUBANhC2/9HaRmhFW9FRKRrUGFpIWfNyRmW0PZ9SwjwLtFvfKXCIiIinZsKSws5ajynFdu7tf+PcpXlAkaxmUuWzmr3fYuIiLQnLeLRQo7aIABs3QLafd/2IQlsZSTFJUXtvm8REZH2pMLSQufYvyCybD9947q3+76HXZwAQKE7hpJ9pXTvH9nuGURERNqD3hJqocd7LuRdLiN1bHW77zuibwR9rJ5Tmrd/UNDu+xcREWkvKiwt5XR6Prb3ec0nDe9+CIBteSWm7F9ERKQ9qLC0kHHCnGsJ1RnetxyAbVtcpuxfRESkPaiwtFB00RZCqWRnUYQp+x8+3HM69bZ9oabsX0REpD2osLRQpRGKg1BCIkJM2f+otHBGsZmBVVtN2b+IiEh70FlCLeB2GVTheSvIFhFsSoZzZ8Sz+e44KLeCc6Zpb02JiIi0Jc2wtEBVqdP7ub2HSUUhJga6dwe3G/73P3MyiIiItDEVlhZwlFR5P7dHmXOWEBYLDB+OGwvOTSosIiLSOamwtIDja88MSwC1BNrMe3ftIccvCKecx/4UbVoGERGRttSswrJ48WISExOx2Wykpqaydu3aRsc+//zznHvuufTo0YMePXqQnp7+nfE33XQTFoul3m3q1KnNidaunKWeGRY7DlNzdIuP5ATd2LbHnAN/RURE2prPhWXp0qVkZWWxYMECNmzYwJgxY8jIyODIkSMNjs/NzWXmzJmsXLmSvLw8EhISmDJlCgcPHqw3burUqRw+fNh7e+2115r3jNpRkLuKKbzPhcGfmppj+DjPKc3bj0aZmkNERKStWAzDMHz5htTUVCZMmMCzzz4LgNvtJiEhgTvvvJP777//jN/vcrno0aMHzz77LLNmea4yfNNNN1FSUsKyZct8fwZAWVkZkZGRlJaWEhHRjuuhfP45TJgACQlQYN7S+Ls/3Megi/oTgpPKqiACgtv/QowiIiK+8uX126cZlurqatavX096evqpDVitpKenk5eX16RtnDhxgpqaGqKi6s8G5Obm0rt3b4YOHcqcOXMoLi5udBtVVVWUlZXVu5nCcfKtIJOW5a+TOLkvITipwsbeTw6YmkVERKQt+FRYjh07hsvlIiYmpt7jMTExFBYWNmkbv/jFL4iPj69XeqZOncpLL71ETk4Ojz/+OKtWrWLatGm4XA0vN5+dnU1kZKT3lpCQ4MvTaD0mX0eoTkBwAENt+wDYlltkahYREZG20K5nCT322GO8/vrrvPXWW9i+scDZ9ddfzxVXXEFycjLTp0/nnXfeYd26deTm5ja4nblz51JaWuq97d+/v52eQX3vfBRBBKVcsvf3puz/m4b38sxIbdtg7gHAIiIibcGnwhIdHU1AQABFRfX/L76oqIjY2NjTfu8TTzzBY489xgcffMDo0aNPO3bAgAFER0eza9euBr8eEhJCREREvZsZKstclBOBw2LuDAvA+cnHmc5bDHBoiX4REel8fCoswcHBjBs3jpycHO9jbrebnJwc0tLSGv2+hQsX8sgjj7B8+XLGjx9/xv0cOHCA4uJi4uLifInX7hyVbgDsQbUmJ4E5syp5i6u5xvGK2VFERERanc9vCWVlZfH888/z4osvsm3bNubMmUNlZSWZmZkAzJo1i7lz53rHP/7448ybN48XXniBxMRECgsLKSwspKKiAoCKigruvfdeVq9ezd69e8nJyeHKK69k0KBBZGRktNLTbBt1hcUW2PCxNu1q+HDPx23bwLcTv0RERPyez8uzzpgxg6NHjzJ//nwKCwtJSUlh+fLl3gNxCwoKsFpP9aDnnnuO6upqrr322nrbWbBgAQ899BABAQFs2rSJF198kZKSEuLj45kyZQqPPPIIISH+vRCa03FyhiXYDwrLkCEYFiuHS+xE7SvClnj6t+hEREQ6kmatJ3/HHXdwxx13NPi1bx8ou3fv3tNuy2638/777zcnhunqzmq2hbjNDQJgs3FW0GY2Vo/gw39+wYVZKiwiItJ56FpCLeCsW4bFHwoL0Kd7JQDb15WbnERERKR1qbC0QF/bMc7mUwZFl5gdBYDh/TyFZdtXOoZFREQ6FxWWFrh1SC6fMpm7J683OwoAw0d5luTfdiDM5CQiIiKtS4WlJepWuv3GInhmGp7WHYBtJfHmBhEREWllKiwt4SfXEqozbEo/AA664yg7YNL1lURERNqACksLzP74x8RzkL9tPv3Kve2lR2IkMdYjAGxfYc7lCkRERNpCs05rFo+jzjAOE0+NdY/ZUbwy+/6XmoJD9DgyABhpdhwREZFWoRmWFnDUePqevZv//BizL/2EJ7iXwSXrzI4iIiLSavznlbYDctQEAWAL86OJqm8u0S8iItJJqLC0gNPlKSz2sACTk3zD8OEcJZq1G/yoRImIiLSQXtVawOEKBsAe7j8/xkM9RtKHowTsr+VERTXBYcFmRxIREWkxzbC0gMPtKQO28CCTk5wSNzaWCEpxEcj/VuwzO46IiEirUGFpgZHW7YzmSyJ7+s8Mi8VqYWRYAQCbPzxqchoREZHWocLSAv8Jmc6XpDBkuB8dwwIk9/0agM3rq01OIiIi0jpUWFrCz1a6rZM8ynPxw827Q01OIiIi0jpUWJqrthZcLs/nfnItoTrJkyMB2FysawqJiEjnoMLSTBVHHfRjH0PYQZXVz2ZYLusPwD5XX11TSEREOgX/OVq0gzlx3Ml+PBcbDI4wTE5TX9TAHtwf9iyJFZuxbMuEvpPMjiQiItIiKizN5Cj1HNAaghOL1b/eEgLInvwuLF8Oe84CVFhERKRj01tCzeQsrQLAbnGanKQRycmej5s3m5tDRESkFaiwNJOjrAbw38LiHDqGNUzk37nhZkcRERFpMb0l1Ex1hcVm9c+1TraHjWcSa+jx1dcUuw0sVovZkURERJpNMyzN5J1hCfDPwjJ8WiIB1PK10YNDGwrNjiMiItIiKizNFOx2MoKtDLIdNDtKg0IiQhgS7LmW0Obl/plRRESkqVRYmmnywMNsZRT/GvVLs6M0KjmmCIDNeRUmJxEREWkZFZbm8tNl+b8peYjn7arN23WokoiIdGwqLM3VEQrLRE+2zYejTU4iIiLSMioszfTSqv6MYCv3/+9HZkdpVPKUOAC2ORKpddaanEZERKT5VFiaqeh4ENsYweEa/529SJzcl0XB9/Eul8LOnWbHERERaTYVlmbyviMU4jY3yGlYA63cNfYjLuJDArdpxVsREem4VFiayXFygVtbiH9d+PA76pbo37TJ3BwiIiItoMLSTE6nZ+VYu82/C8uRxIn8jRt46T3/fetKRETkTHS+azM56gqL/54kBMBXoeP5AbeQtLmAWWaHERERaSbNsDSTs9rzo7PZTA5yBsmX9gMgv7Yf5YfKTU4jIiLSPCoszdTDWkoCBURFusyOclo9h/Qkzuq5ltDW9/aZnEZERKR5mlVYFi9eTGJiIjabjdTUVNauXdvo2Oeff55zzz2XHj160KNHD9LT078z3jAM5s+fT1xcHHa7nfT0dHb6+Wm4/2/YHyigPz++KN/sKGeU3OMAAJs/Om5yEhERkebxubAsXbqUrKwsFixYwIYNGxgzZgwZGRkcOXKkwfG5ubnMnDmTlStXkpeXR0JCAlOmTOHgwVMX5Fu4cCHPPPMMS5YsYc2aNXTr1o2MjAycTmfzn1lbq8vm7+8JAcmJnmsJbf7Svw8QFhERaYzPheWpp57illtuITMzkxEjRrBkyRJCQ0N54YUXGhz/t7/9jdtuu42UlBSGDRvGn/70J9xuNzk5OYBndmXRokU8+OCDXHnllYwePZqXXnqJQ4cOsWzZshY9uTbVAZbmr5M81nNs9eaCCJOTiIiINI9PhaW6upr169eTnp5+agNWK+np6eTl5TVpGydOnKCmpoaoqCgA8vPzKSwsrLfNyMhIUlNTG91mVVUVZWVl9W7tbebWB5nIGvL2xrX7vn2VfEFPADaX9sNwa5ZFREQ6Hp8Ky7Fjx3C5XMTExNR7PCYmhsLCwiZt4xe/+AXx8fHeglL3fb5sMzs7m8jISO8tISHBl6fRKjafGMA6JuLA/2dYRlySyNuWK1lvnAVN/D2JiIj4k3Y9S+ixxx7j9ddf56233sLWgmM/5s6dS2lpqfe2f//+VkzZNA5XMAD2cP9fysbWw87lQ3bQnwIsW7REv4iIdDw+FZbo6GgCAgIoKiqq93hRURGxsbGn/d4nnniCxx57jA8++IDRo0d7H6/7Pl+2GRISQkRERL1be3O4QwCwRQS3+76bpW6J/s0qLCIi0vH4VFiCg4MZN26c94BZwHsAbVpaWqPft3DhQh555BGWL1/O+PHj630tKSmJ2NjYetssKytjzZo1p92m2ZyGp7DYI4JMTtI022Iv5GHms/iN3mZHERER8ZnP72dkZWUxe/Zsxo8fz8SJE1m0aBGVlZVkZmYCMGvWLPr06UN2djYAjz/+OPPnz+fVV18lMTHRe1xKWFgYYWFhWCwW7r77bh599FEGDx5MUlIS8+bNIz4+nunTp7feM21lDsPzlpY9smPMsGwPPYuHmMRZm7Zxu9lhREREfORzYZkxYwZHjx5l/vz5FBYWkpKSwvLly70HzRYUFGC1npq4ee6556iurubaa6+tt50FCxbw0EMPAXDfffdRWVnJrbfeSklJCZMnT2b58uUtOs6lLRluA+fJg23t3UNMTtM0yVPiYCF85UjEVe0iIDjA7EgiIiJNZjEMo8Of51pWVkZkZCSlpaXtcjxLVamTpO7HcWJj7/5AIvr6//om7lo34UEOTtCN7e/tYei0AWZHEhGRLs6X129dS6gZQtwODtGH4/QkIsb/T2sGsAZaGdltLwCb/1t0+sEiIiJ+RoWlOepWuQ0IgKCOcdAtQHK851pCm9ZVmZxERETENyoszVF3HaEOsCz/NyWPdAOweVfHyi0iIqLC0gy7d9RyNp9yTfVrZkfxSfLZ4QDsONbT5CQiIiK+UWFphq+P1pLH2axzjTU7ik/SruvLVwxnU81wqKw0O46IiEiTqbA0g6O8FgCbtdrkJL4JTezN8JivCaQWtm41O46IiEiTqbA0g7O8BgB7QMcqLICW6BcRkQ5JhaUZ6mZY7AE1Jifx3cruV/FDXuLJF3Uci4iIdBwqLM3gqPScbWMP6niFZW94Mq/wQ97b0s/sKCIiIk2mwtIMjgoXALbAWpOT+C75As/MyuaSBJOTiIiINJ0KS3PU1BBOGeEhHe8YlhGXJGLBzVGjF0VbjpodR0REpElUWJrhprFfUkYkfz//92ZH8VlodCiDgvYBsPm9/SanERERaRoVluaoW5rfT68mfSbJvTzXEtr8WZnJSURERJpGhaU56gpLB1uav07yYM+lBTZ/FWByEhERkaZRYWmG5/JSmML7/GXXuWZHaZbkCSFYcVF+TBdBFBGRjkGFpRm2FkWzginsqYwxO0qzXDY7mkq68YbzCnC5zI4jIiJyRioszeCssgAd9h0hQoYPwGa3et7a2rXL7DgiIiJnpMLSDI5qz4/NHmpykOYKCICUFM/n69aZGkVERKQpVFiawVntOVjVZu+4P763ev6YSeRx71NxZkcRERE5o0CzA3REjhpPYbF367iFxZk4jDVMwvjfFrOjiIiInFHHfcU1kfNkYbF14MKSek1fADZWDqaqTGcLiYiIf+u4r7gmMtwGFtzYu3XcdUySzksg2nKMakLY+MZOs+OIiIiclgpLM6xMuhkXAVx+XqnZUZrNYrUwMTofgDX/KTY5jYiIyOmpsDSHw4EFsHbroOc1n5Q6qhKAtet1KJOIiPg3FZbmcHqWtu+o1xKqk3pRGABrDvYxOYmIiMjpqbA0ww1FT3EN/2B/aYTZUVpkwsxB9GU/o2vWU1Oot4VERMR/qbA0wztVF/NPrqHK0rFnWKIGdGf/kHTe5FqCNmoBORER8V8qLM3gwHPsij0y2OQkrSA11fNxzRpzc4iIiJyGCouPap211BIEgL1Hx55hAWDiRAzgyEfbzU4iIiLSKJ0e4iPHcQcQDoCte8cvLDvjzuMcinCvDOCo28BitZgdSURE5Ds0w+IjZ+mpVWFtkSEmJmkd/TOGUUokxUZP9uQWmB1HRESkQSosPnKUeApLCE6sgR3/xxccFszYMM9Kt2vePGByGhERkYZ1/FfcdlY3w2Kj81x/J3Wg55TmNZ/WmpxERESkYSosPhoSV04tARyKGWt2lFaTeo7nUKY1u3qanERERKRhKiy+cjgIwE1oqNlBWk/dlZu/qBxMdXnnmTkSEZHOo1mFZfHixSQmJmKz2UhNTWXt2rWNjt26dSvXXHMNiYmJWCwWFi1a9J0xDz30EBaLpd5t2LBhzYnW9hwOz0d7x76O0DcNuKAfPS3FVBPCl//cbXYcERGR7/C5sCxdupSsrCwWLFjAhg0bGDNmDBkZGRw5cqTB8SdOnGDAgAE89thjxMbGNrrdkSNHcvjwYe/tk08+8TVau1i/OZgZvM5Dx39qdpRWY7Fa+HFSDveykB571psdR0RE5Dt8LixPPfUUt9xyC5mZmYwYMYIlS5YQGhrKCy+80OD4CRMm8Nvf/pbrr7+ekJDGTwMODAwkNjbWe4uOjvY1WrvYdyCAvzOD/1ammR2lVT02axsL+QWD9nxgdhQREZHv8KmwVFdXs379etLT009twGolPT2dvLy8FgXZuXMn8fHxDBgwgBtvvJGCgsbXBKmqqqKsrKzerb04KlwA2INq2m2f7UJL9IuIiB/zqbAcO3YMl8tFTExMvcdjYmIoLCxsdojU1FT++te/snz5cp577jny8/M599xzKS8vb3B8dnY2kZGR3ltCQkKz9+0rR6UbAFugq9322S4mTqSESFbs7E/5vuNmpxEREanHL84SmjZtGtdddx2jR48mIyOD9957j5KSEv7+9783OH7u3LmUlpZ6b/v372+3rM4TnsJiD+5khSUqitSgDUxhBZ+9ssfsNCIiIvX4VFiio6MJCAigqKio3uNFRUWnPaDWV927d2fIkCHs2rWrwa+HhIQQERFR79ZeHCcMoBMWFmBCn8MArMmpMDmJiIhIfT4VluDgYMaNG0dOTo73MbfbTU5ODmlprXcQakVFBbt37yYuLq7Vttla6s5qtgW7zQ3SBiaO86x0u2ZLN5OTiIiI1OfzW0JZWVk8//zzvPjii2zbto05c+ZQWVlJZmYmALNmzWLu3Lne8dXV1WzcuJGNGzdSXV3NwYMH2bhxY73Zk3vuuYdVq1axd+9ePvvsM6666ioCAgKYOXNmKzzF1uV0ej7abYa5QdpA6iWeM7PWHkvCcHe+5yciIh1XoK/fMGPGDI4ePcr8+fMpLCwkJSWF5cuXew/ELSgowGo91YMOHTrE2LGnlrF/4okneOKJJzj//PPJzc0F4MCBA8ycOZPi4mJ69erF5MmTWb16Nb169Wrh02t9v5rwb+7PuwLr934GnG92nFaVcu0ggm+u4pgRTf5HBQy4oJ/ZkURERIBmFBaAO+64gzvuuKPBr9WVkDqJiYkYxun/b/31119vTgxTBFafIIJyiPCL45VbVUhECCndtrC2chRr3jygwiIiIn6j873qtrVOuDT/N9VduXntJ9UmJxERETmlWTMsXdnvNp3PRs7lhwd6cIHZYdrAzKucjNj0Ey6gFDrlMxQRkY5IMyw+WnFoJC9wMztL/PPSAS2VduMAfsIfGLbtLajWLIuIiPgHFRYfOWsCALB3CzA5SRsZNAiioqCqCjZtMjuNiIgIoMLiM0et5100e1gn/dFZLOwcOZ0l/B8fvnTA7DQiIiKACovPnLVBANi6dd7Df16x/JA5LOGv7/rfaeUiItI1qbD4yOHyFBZ7eOctLKnf86x0u2Z/vMlJREREPFRYfOR0BQNg68SFZeL1AwD4X00SX+eXmBtGREQEFRafOdyewmIPDzI5SduJHtqTAYEFAKx7fbfJaURERFRYfLY1PI2DxDNypNlJ2lZqH88Bt2v+W25yEhERERUWn3WvKiKewwRHds6VbuuknlUD6MrNIiLiH1RYfGEYpy7XbLOZm6WNpV7SE4B1RxN15WYRETGdCosPqsur+D+WcBeLqLJ27hmWsdcN4oOAaewwBmPZm292HBER6eJUWHxQWezkj/wfz3AXltDOXVhCIm1cPKmc7pRCTo7ZcUREpItTYfGBs7QKACsuguyd97RmrylTPB8/+MDcHCIi0uV1gVfd1uMo8RQWOw4s1jCT07S9islT+Q2B5C07lxXOWgJt+uciIiLm0AyLDxylnqsX2yxVJidpH/bJ43jOchu5tefy+SvbzY4jIiJdmAqLD5zlnlN97dauUVgCggNIj98GwPt/O2ZyGhER6cpUWHzgKOtahQUg46JaAD5YH2VyEhER6cpUWHxQV1hsATUmJ2k/U/4vCYA15SMo2VdqchoREemqVFh8cN6QQvbRj3eG3mN2lHbT7+y+DAvejYtAPly8zew4IiLSRamw+CCktpJ+7Kdf9zKzo7SrKcM91xX64J1qk5OIiEhXpcLiC4fD89HeuReN+7Yp00OJ5ihhB3d4Lk8gIiLSzrSwhg8+3tydZTzBuKMWbjA7TDuamjWCol9HYS2rhl0XwODBZkcSEZEuRjMsPli/uztP8XPeOZZqdpR2FRDRDeu553juaNVbERExgQqLDxwnPG+H2IPdJicxwZQpGMD+ZevNTiIiIl2QCosP6g5hsYV0vcKyf8xl9OUAI/77NNUVOvhWRETalwqLD5xOz0d7SNc78LTPxSOosQRTQTirX/jK7DgiItLFqLD4wOG0AF3uJCEArIFWLu7/PwDeX1pibhgREelyVFh84Kjy/LhsNpODmCTjYs/M0gcbe5ucREREuhoVFh84qz0/LnuoyUFMcvFtntOZ158YxrEdxSanERGRrkSFxQdPjvgzWxnBD8/JNzuKKeJSYki2/Q8DK/9dvMPsOCIi0oWosPgghiJGsI1evcxOYp6MUYcA+GC5y+QkIiLSlWilW1900aX5v+nqH3bD/fmTTD/+CRiTwWIxO5KIiHQBKiw+WJx/CUc4nx98HU1XXZw+7dZk0n5xHhQ74auvYORIsyOJiEgX0Ky3hBYvXkxiYiI2m43U1FTWrl3b6NitW7dyzTXXkJiYiMViYdGiRS3eplmeL7ycX7GA/K+7mx3FPDYbnH++5/P33zc3i4iIdBk+F5alS5eSlZXFggUL2LBhA2PGjCEjI4MjR440OP7EiRMMGDCAxx57jNjY2FbZplkcrmAAbOFBJicxl/N7l7CCdP745wCzo4iISBfhc2F56qmnuOWWW8jMzGTEiBEsWbKE0NBQXnjhhQbHT5gwgd/+9rdcf/31hISEtMo2zeJ0e4qKPaJrF5bdQy9hCiu466tbcXztNDuOiIh0AT4VlurqatavX096evqpDVitpKenk5eX16wAzdlmVVUVZWVl9W7tweH2rBhnC+vah/6MuHwgfayHcWLnkz9sNTuOiIh0AT4VlmPHjuFyuYiJian3eExMDIWFhc0K0JxtZmdnExkZ6b0lJCQ0a9++chieGSJ794ZniroKi9XClIG7APjgzXKT04iISFfQIddhmTt3LqWlpd7b/v3722W/Tk7OsEQEt8v+/NmUqZ7jV97fEm9yEhER6Qp8KizR0dEEBARQVFRU7/GioqJGD6hti22GhIQQERFR79bWap211HLyGJYuPsMCkH77UCy42ewcwuGNRWf+BhERkRbwqbAEBwczbtw4cnJyvI+53W5ycnJIS0trVoC22GZbsFY72UQya5hIZEwXvfrhN0QP7cm40O0AfPD7nSanERGRzs7no0ezsrKYPXs248ePZ+LEiSxatIjKykoyMzMBmDVrFn369CE7OxvwHFT71VdfeT8/ePAgGzduJCwsjEGDBjVpm/7AWuUgmS2eO2EqLABTUo7w+WcjyFtZxWyzw4iISKfmc2GZMWMGR48eZf78+RQWFpKSksLy5cu9B80WFBRgtZ6auDl06BBjx4713n/iiSd44oknOP/888nNzW3SNv1C3bL8wcFg7ZCH/rS6/7vbzvWfJTOqpAjchfq5iIhIm7EYhmGYHaKlysrKiIyMpLS0tM2OZzm2ehfPpb1Id1sVdzoWtsk+OpzqaujZEyoqYPVqSE01O5GIiHQgvrx+63+Jm+hggYv5PMKvq+8xO4r/CA6GK64AwPXS30wOIyIinZkKSxM5y2sAsFurTE7iX45P/xE38RcGLfk51RXVZscREZFOSoWliRzltQDYrXpR/qaIy8/nfes09rr7szz7C7PjiIhIJ6XC0kSOMs8Miy2gxuQk/iXQFsiNY7cB8NKLbpPTiIhIZ6XC0kTOShcA9kDNsHzbD3/hWe323wfP4vjur01OIyIinZEKSxM5KuoKS63JSfzPmOuGMNq2g2pC+Pu8zWbHERGRTkiFpYkclZ63O2xBLpOT+KdZF3suVPnyO93NDSIiIp2SCksTXTF0B59yNtljXjc7il+64VfDsOLis/LR7MrZZ3YcERHpZHxe6bar6hVUQi/yIHaQ2VH8UlxKDLcmvEuf/XlE/CcCLrrP7EgiItKJaIalqeqW5rfbzc3hx557vIwH+TW933wO3DpjSEREWo8KSxPlbo/l/3E3nx4fbnYU/3XllRAeDnv3wqefmp1GREQ6ERWWJnp7+2Cy+H/8e3+K2VH8V2gozqtm8gbXsvj+ArPTiIhIJ6LC0kQOp+dHpXeETm/12Dl8nzd44LPLcBx3mB1HREQ6CRWWJnJUeX5UNpvJQfzceXeMpl/AAcqI5O1fbTQ7joiIdBIqLE3krD45wxJqMTmJf7MGWvnhpF0AvPxagMlpRESks1BhaSJHjefF195NheVMfvjLfgAsP3IWRVuOmpxGREQ6AxWWJnJUe5assYXqR3YmQ6cNILXbFlwE8tq8r8yOIyIinYBefZvIWespLPZuepujKWZdVgzAS+/3NjmJiIh0BiosTfR034W8zxQmj600O0qHMOORZIKoJthRStlqzbKIiEjLqLA00VjLl0xhBbFxOoalKXoOjmLf1J+wmjQi3nrR7DgiItLBqbA0ldPp+aiFWJos7tbLPZ+88gq4dJVrERFpPhWWJnqh+Aqe58ccrw4zO0rHccklEBXF8UMO9r+upfpFRKT5VFia6J7SedzK8xw5ocLSZCEhPD/6GeI4zAPz9E9NRESaT68iTeQwPEvc2iKCTU7SsSTfOIZqQvhn/ljKD5WbHUdERDooFZYmMNwGTjzHrti7h5icpmNJ/dFIhgbv4QTd+MOt682OIyIiHZQKSxNUlVV5P1dh8Y3FauH+Hx4CYOF7I6k8otPCRUTEdyosTeAoOVVYbN119UNf/eDZSQwM3MdRoxe/z1xndhwREemAVFiawFniOaXZioug0CCT03Q8gbZAHpy9H4CF/xmlWRYREfGZCksT1M2w2HFgsWrhuOaom2WpNEJZN+9ts+OIiEgHo8LSBDFhlbzHNF4Pv9XsKB1WoC2Q1x7YzF4SueCfP4WKCrMjiYhIB6LC0gTdLCeYxnIui/zY7Cgd2oR5U+k9MAKOHYPf/97sOCIi0oGosDSFw+H5qGX5WyYwEObNA+DT36yiolCzLCIi0jQqLE1wYL/BX7iJ/9Skmx2l47vxRm6NeI3Jpe+yOPNzs9OIiEgHocLSBF98FcKP+AsLjt5udpSOLzCQyTP6AvDb95M1yyIiIk3SrMKyePFiEhMTsdlspKamsnbt2tOOf+ONNxg2bBg2m43k5GTee++9el+/6aabsFgs9W5Tp05tTrQ24ayoBcAeWGNyks7hhmcmMTgon2KjJ89qlkVERJrA58KydOlSsrKyWLBgARs2bGDMmDFkZGRw5MiRBsd/9tlnzJw5k5tvvpkvvviC6dOnM336dLZs2VJv3NSpUzl8+LD39tprrzXvGbUBR6UbAHtgrclJOodAWyDzfnQQgCfeT9Y1hkRE5Ix8LixPPfUUt9xyC5mZmYwYMYIlS5YQGhrKCy+80OD4p59+mqlTp3LvvfcyfPhwHnnkEc466yyeffbZeuNCQkKIjY313nr06NG8Z9QG6gqLLchlcpLOY+aiU7Msi3+kawyJiMjp+VRYqqurWb9+Penppw4+tVqtpKenk5eX1+D35OXl1RsPkJGR8Z3xubm59O7dm6FDhzJnzhyKi4t9idamnCdOzrAEq7C0lnqzLB9olkVERE7Pp8Jy7NgxXC4XMTEx9R6PiYmhsLCwwe8pLCw84/ipU6fy0ksvkZOTw+OPP86qVauYNm0aLlfDBaGqqoqysrJ6t7bkPas5RIWlNc1cNIkhQfmEG2Xszv672XFERMSPBZodAOD666/3fp6cnMzo0aMZOHAgubm5XHTRRd8Zn52dzcMPP9xu+U4VFne77bMrCLQF8m72l/S/5zqCXo2A33wfwsPNjiUiIn7IpxmW6OhoAgICKCoqqvd4UVERsbGxDX5PbGysT+MBBgwYQHR0NLt27Wrw63PnzqW0tNR7279/vy9Pw2czBn7OUr5P5pgv2nQ/XdGguy4jaMgAOH4cvnVck4iISB2fCktwcDDjxo0jJyfH+5jb7SYnJ4e0tLQGvyctLa3eeIAVK1Y0Oh7gwIEDFBcXExcX1+DXQ0JCiIiIqHdrSyPCCvg+bzA+8Vib7qdLOrn6bQ2BvPDrQ5QWlJqdSERE/JDPZwllZWXx/PPP8+KLL7Jt2zbmzJlDZWUlmZmZAMyaNYu5c+d6x991110sX76cJ598ku3bt/PQQw/x+eefc8cddwBQUVHBvffey+rVq9m7dy85OTlceeWVDBo0iIyMjFZ6mi2kpfnb1syZXNkth5srf8dd39tsdhoREfFDPheWGTNm8MQTTzB//nxSUlLYuHEjy5cv9x5YW1BQwOHDh73jzz77bF599VX++Mc/MmbMGP7xj3+wbNkyRo0aBUBAQACbNm3iiiuuYMiQIdx8882MGzeOjz/+mJCQkFZ6mi3zyb4E3uBa9pxo/G0saYGAAB5Y2B0rLl7cPZk37234jDMREem6LIZhGGaHaKmysjIiIyMpLS1tk7eHru2bx5sH01g8YxW3vX5+q29fPB44O5fsvAuIshxn8+fVxJ+lgigi0pn58vqtawk1gaPaczKVPVQ/rrb00Adnc5Z9G8eNKH405QCGu8N3aRERaSV6BW4CR42nsNi6BZicpHMLDgvmlaVB2HDwfvF4Fl//sdmRRETET6iwNIGz9uQMS5gKS1sbfvkgFl6zBoAFb4ykYv0OkxOJiIg/UGFpAocrCFBhaS+3v34eP+v3Jp8wmbD/uxGqq82OJCIiJlNhaQKHKxgAe0SQyUm6Bmuglac+m8TwHkWwfj386ldmRxIREZOpsDRBXWGxhfnFlQy6hj594A9/AOCz3+Sy9k+bTA4kIiJmUmFpgkVhD/JnfsSAgRazo3Qt113Hm+c9zbnGKmbOidQVnUVEujAVlia4yv1PfsRfiI4PNjtKl5P+8mwSAg6xp7Y/P7two9lxRETEJCosTeF0ej7abObm6IIi+0Xy4pPFWHDz5/+dy7K5a8yOJCIiJlBhOQN3jYu3ai7lP0ylJlDXEjLD+XelcM+EjwC45fGB7P3kgMmJRESkvamwnMGJYgdX8xaX8B8VFhM98t80UuzbOWZEc+GFBvs+VWkREelKVFjOwFFS5f3c1l1vCZklJCKEdz+OZHBQPntrE3jrir/A/v1mxxIRkXaiwnIGzlJPYQmmCmugflxmih8Xx8pPQ3i656+4+/h8uPBCOKCZFhGRrkCvwGdQN8Nix2lyEgHoMyGen36RCUlJsHs3Jy64hMKNhWbHEhGRNqbCcgaOUs+y8HarCovfSEiAlSup7DecS3c/zQUTT3B4Y5HZqUREpA2psJyBo7wWAJtV17PxK/37c/y199kTMJgdNQP4XmoFhZuOmJ1KRETaiArLGTjLawCwq7D4nYSzE1j5XzcJAQfZXj2QiyaWcWTrUbNjiYhIG1BhOYPBPY/zB25lXtyfzI4iDRhwQT9WflBLH+thvqoaxPfGlXB02zGzY4mISCtTYTmDPmGl3MrzzOzzkdlRpBEDv9eflcuriLceZmvVYC466zgFn+qUZxGRzkSF5UwcDs9HuxaN82eDL07kw3ccxFqLOOjsScDUi+GVV8AwzI4mIiKtQIXlDAoOWFlBOltrhpgdRc5g6LQBfLaymjdGPkSfih3wwx/CDTfgOFxidjQREWkhFZYzePvzeKawggW7Z5kdRZog6bwEvrfx/8GvfgUBAfzr9RMMS6hg1dMbzY4mIiItoMJyBs4TbgDswS6Tk0iTBQbCvHkYn3zK4yELKHD15cK7RzM3LZfqCp3tJSLSEamwnIH3EJYQFZaOxjIplQ/2DuHmIR9jYOWx1ReQ1ns329/bY3Y0ERHxkQrLGThPLnBrD9HBmx1RWGwYf9pxLm/eu5ooy3E2OIZz1qWxPDfzI1zVKqEiIh2FCssZOE4WFptNhaUju3rhJDZ/Xs3FUetxEMptr5/H1iFXwcsvQ22t2fFEROQMVFjOwOm0AGC3mRxEWiz+rFiWF43l6atXMTv4VUbv+zfMmgXDhvHp3HeoOVFjdkQREWmECssZOKo9PyItw9I5WAOt/PTN8/nr0csgOxuio9m/u4oLH5vC4IgiltzwEVVlVWbHFBGRb1FhOYMb+nzEk2TxvRGFZkeR1hQRAfffD3v3svMnTxFlLWGfqy9zXjuPgT2O87trV+E47jA7pYiInKTCcgYXR64li//HhKFlZkeRttCtG9977jryj4bz9DWriLcWctAdx0/fPJ/E6AruHL2KY//OA7fb7KQiIl2aCsuZeE8T0ntCnZk9ys5P/3E+u4/34PfXf0S/gAMcMXrxl83j6HbF9yAhAe66i4P/+hx3rcqLiEh7U2E5gw1HE8hjEiXuCLOjSDuwRYYw57Xz2FUWw7sPreOJSW9ijwyBQ4fgmWeYMt1Of1shWeNyWf2nLRhunT0mItIeLIbR8a8OV1ZWRmRkJKWlpUREtG6xmNhtK+tOjOTf89Zy2a8mtuq2pYOoqoIVKzj8l+UM/edvKOfUv7G+AYc4r28+Z0+s5aKZvRl2xRAICDAxrIhIx+HL67dmWM7A4QoGwB4eaHISMU1ICFx2GXFvPsuRr4NZNncNN/T/lDDKOeCK59V953DHG+fz+6tXQI8eMGUKzgcfZcVj6yk/VG52ehGRTkGvwmfgLSwRQSYnEX9g627jyt+kcuVvwHHcwad/2sBn/ynl083hXFT5GZSXw4oVrF3hZAoPYp3rYpRtByN6FzNsQDXDRgcz7JyeDLkoAXvPULOfjohIh9GsGZbFixeTmJiIzWYjNTWVtWvXnnb8G2+8wbBhw7DZbCQnJ/Pee+/V+7phGMyfP5+4uDjsdjvp6ens3LmzOdFancPtKSy2MHU7qc8eZSf9vrOYv/JC3j82nisr/gYbN8LixZScdyWJgftxE8Am51BeLzibh3Iv4PpnziZlxlBejf4pJCbC1KnsuulRfn/9Kt7+5RrWv7KNwk1HdGCviMi3+PwqvHTpUrKysliyZAmpqaksWrSIjIwMduzYQe/evb8z/rPPPmPmzJlkZ2dz2WWX8eqrrzJ9+nQ2bNjAqFGjAFi4cCHPPPMML774IklJScybN4+MjAy++uorbDZzl5h1GiEA2CODTc0hHUBAAIwZA2PGcMVtcAVwaEMhn7+5jx0bHWzfGcD2wki2V/RlqLEd9u2DfftYSQK382C9TQVSQ1zAYeJtX7Nw3FLOG1MK8fHscSfyxdf9iYqzEdU3lKj+4UQlRRIaHYrFajHneYuItAOfD7pNTU1lwoQJPPvsswC43W4SEhK48847uf/++78zfsaMGVRWVvLOO+94H5s0aRIpKSksWbIEwzCIj4/n5z//Offccw8ApaWlxMTE8Ne//pXrr7/+jJna8qDbUMsJHISS//EBEif3bdVtS9dlHD2GZcd22L6df79r5c+rR3KwLJyDziiK3L1wc+rA3Q+5kAvJBeAP3MpP+MN3thdMFVHWEv7adx4ZsV9CeDirnKk8f/ASwuwuwkLdhIVBWBh0C7di72blwjHH6Z/gBrudo85wdh+LxBYWiD0iCFtEMEH2QILsgQR3CyI0Moig0CAIDASLipGItA5fXr99mmGprq5m/fr1zJ071/uY1WolPT2dvLy8Br8nLy+PrKyseo9lZGSwbNkyAPLz8yksLCQ9Pd379cjISFJTU8nLy2uwsFRVVVFVdWr59LKytlnUzXAbOPAcZ2CL0AyLtB5Lr2joNRkmT+byH8Pl3/harbOWws2HObS5mIM7KhjTYyaUpEJhIT23JDF555ccd3ajuDaC4+7u1BBMNSEUumMILtgJBZ63aLczkL9xTqMZ3mI6/fkXAB9wAz/gb42OfYUbuZFXAXg74Cpmul4miFoCLbUEWlwE4PJ8tLh5LPpJvt/9AwgM5DPnWdx58H6sFoMAixurxcBqPfX5nQn/4urYz8BqZUtFIr/4381YLGC1uE9+NLBawGIxuLHfJ1zVdx1YLORX9uaBzTOxWAwseDqUxWKc7FIWpies5+rEDQAcdnTnl19c4/m5U79vWSwGU/p8xXUDPWOPV3XjgXXTT46sG3Tq3rlxu7hhyOcAVNTa+EXelY3+zFJ772XW8HUA1LisZH16TYPjLBYYE32Qm0eu9j5210dXYxgNF8PhUUXMSf7Ee/++T6+gyhVYP/NJAyOP8dOUj7z3H8y7hPLqhmet+4aVcO+4D733f7Umg2JntwbHxoSW88CEFd77j32ezuHKhl9seoSc4KFJy733n9pwAQUVPU7eq585NLCa35zzrvf+4i8ns7OkV4PbDbK6+O25b3vvP79lEluL4xocC/DUecuwWjz/f/7Xryay8WifRsc+ds6/sQV6Lor62o6zWFPYv9Gxj6S9R3iw5/XozZ2j+fjQwEbHPjjxA6LtlQC8vWckH+4f0ujY+8blEB/meW17f98w/rN3eKNjfzY2l/4RXwOwcv8g/rUnudGxt43+hCE9jgLwycEk/rErpdGxPx6Zx6hozyrvW4tjGRlzDJ54otHxbc7wwcGDBw3A+Oyzz+o9fu+99xoTJ05s8HuCgoKMV199td5jixcvNnr37m0YhmF8+umnBmAcOnSo3pjrrrvO+P73v9/gNhcsWGAA37mVlpb68nTOyHXCaSzkHuNh5hknDpe06rZFWoPb5TbKD5cb+z47YHzx+naj7J1VhvH224bx6qvGl/PeMJ68YqXx8IUrjXsnrDTmjFxl/HDAx8ZVcXnGtF5rjbXjfmIYaWmGMXas8Y8+PzWSAvcZcdbDRg/LccPGCSOAGgMMAwxjKdcZdXeWcp338YZuL3CT9867TDvt2N9xu/dOLuedduxj3Oe9s45xpx07j4e9d75i2GnH/pzfeu/sI+G0Y3/C7713jhF12rE/5EXvnRPYTjv2av5R7wELrkbHTuW9eg+EUdbo2HNZVe+B3hQ2OvYsPq/3QBK7Gx07jK/qPTCKTY2OTWBfvQcmsrrRsVEcq/fAheQ0OtbGiXoPXMbbp/0Zu7B473yf1087toww751M/nzasYX09t65g2dOO3Y3Sd479/HYacduYpT3zsPMO+3Y1Uz03nmCrNOOzeFC753f85PTjn2by7x33uESwwgJafX/hpWWlhpNff3ukEeSzp07t96sTVlZGQkJCa2+H6sV7p0f6lntVmd0iB+yWC2ExYYRFhtGvzSAod6vjT55a9wE72fXnLx9m7vWTa2zFmvtC2D8AWpquLSkhj2F+6murMFV48ZV46a22nNz1bgZEPMTCJ8NLhcTjll4b8c63C4Dt8vA5eLk525ctXBW/wzoPQkMg6HHg/nL5k9wuwwMt4HbAMPtuSqC220wKfEs6Ps7MAwSSm0s2rDq5H9WDc9Ht+e/rABp/QdB/8cB6FVhI3tdLhinvm584/PUvnEw6FEwDCKdIfwqL/fUD8AwMOpmAQyDcfGRMOQhMAzsNUEs+DT3m0PrGdM7CIbPAyDQZWXex7l8W933jIh2wahTxzE9mLvKu99vb3dw1AkY/Uvv/Xs/Xku1q+G1f/pHlsDYU2Pv/vRLKmsani2OCyuHcQ9479+2egdfOwsaHNsrtBImnhp7y7p9FFUWNzg2MsQJaafG3rS+kPTy3AbH2gNrYPKpsTduLCGtpOGxgVY3nHdq7Pc3Oxld3PBYAM5/AE7OsFy11WDQ0cbHBk/+GQS6ALh0u424wsbHdjv7NgiuBmDKzu5EHGx8bPfUH4Hdc42yC3fHELi/8bG9xt8AYRUATN7bnwf2Nj42/qxrIMLzLkXq/v48sLvxsf3HXAo90gA462BfHtjZ+NhBoy6CaM9/RQYVR0PM+EbHtgefjmGprq4mNDSUf/zjH0yfPt37+OzZsykpKeFf//rXd76nX79+ZGVlcffdd3sfW7BgAcuWLePLL79kz549DBw4kC+++IKUlBTvmPPPP5+UlBSefvrpM+Zqy2NYREREpG202cJxwcHBjBs3jpycHO9jbrebnJwc0tLSGvyetLS0euMBVqxY4R2flJREbGxsvTFlZWWsWbOm0W2KiIhI1+LzW0JZWVnMnj2b8ePHM3HiRBYtWkRlZSWZmZkAzJo1iz59+pCdnQ3AXXfdxfnnn8+TTz7JpZdeyuuvv87nn3/OH//4RwAsFgt33303jz76KIMHD/ae1hwfH19vFkdERES6Lp8Ly4wZMzh69Cjz58+nsLCQlJQUli9fTkxMDAAFBQVYracmbs4++2xeffVVHnzwQR544AEGDx7MsmXLvGuwANx3331UVlZy6623UlJSwuTJk1m+fLnpa7CIiIiIf9DFD0VERMQUuvihiIiIdCoqLCIiIuL3VFhERETE76mwiIiIiN9TYRERERG/p8IiIiIifk+FRURERPyeCouIiIj4PRUWERER8Xs+L83vj+oW6y0rKzM5iYiIiDRV3et2Uxbd7xSFpby8HICEhASTk4iIiIivysvLiYyMPO2YTnEtIbfbzaFDhwgPD8disbTqtsvKykhISGD//v26TpGJ9HvwD/o9+Af9HvyDfg8tZxgG5eXlxMfH17twckM6xQyL1Wqlb9++bbqPiIgI/YP0A/o9+Af9HvyDfg/+Qb+HljnTzEodHXQrIiIifk+FRURERPyeCssZhISEsGDBAkJCQsyO0qXp9+Af9HvwD/o9+Af9HtpXpzjoVkRERDo3zbCIiIiI31NhEREREb+nwiIiIiJ+T4VFRERE/J4KyxksXryYxMREbDYbqamprF271uxIXcpDDz2ExWKpdxs2bJjZsTq9jz76iMsvv5z4+HgsFgvLli2r93XDMJg/fz5xcXHY7XbS09PZuXOnOWE7sTP9Hm666abv/H1MnTrVnLCdVHZ2NhMmTCA8PJzevXszffp0duzYUW+M0+nk9ttvp2fPnoSFhXHNNddQVFRkUuLOS4XlNJYuXUpWVhYLFixgw4YNjBkzhoyMDI4cOWJ2tC5l5MiRHD582Hv75JNPzI7U6VVWVjJmzBgWL17c4NcXLlzIM888w5IlS1izZg3dunUjIyMDp9PZzkk7tzP9HgCmTp1a7+/jtddea8eEnd+qVau4/fbbWb16NStWrKCmpoYpU6ZQWVnpHfOzn/2Mf//737zxxhusWrWKQ4cOcfXVV5uYupMypFETJ040br/9du99l8tlxMfHG9nZ2Sam6loWLFhgjBkzxuwYXRpgvPXWW977brfbiI2NNX772996HyspKTFCQkKM1157zYSEXcO3fw+GYRizZ882rrzySlPydFVHjhwxAGPVqlWGYXj+7QcFBRlvvPGGd8y2bdsMwMjLyzMrZqekGZZGVFdXs379etLT072PWa1W0tPTycvLMzFZ17Nz507i4+MZMGAAN954IwUFBWZH6tLy8/MpLCys97cRGRlJamqq/jZMkJubS+/evRk6dChz5syhuLjY7EidWmlpKQBRUVEArF+/npqamnp/D8OGDaNfv376e2hlKiyNOHbsGC6Xi5iYmHqPx8TEUFhYaFKqric1NZW//vWvLF++nOeee478/HzOPfdcysvLzY7WZdX9+9ffhvmmTp3KSy+9RE5ODo8//jirVq1i2rRpuFwus6N1Sm63m7vvvptzzjmHUaNGAZ6/h+DgYLp3715vrP4eWl+nuFqzdF7Tpk3zfj569GhSU1Pp378/f//737n55ptNTCZivuuvv977eXJyMqNHj2bgwIHk5uZy0UUXmZisc7r99tvZsmWLjqMziWZYGhEdHU1AQMB3jvQuKioiNjbWpFTSvXt3hgwZwq5du8yO0mXV/fvX34b/GTBgANHR0fr7aAN33HEH77zzDitXrqRv377ex2NjY6murqakpKTeeP09tD4VlkYEBwczbtw4cnJyvI+53W5ycnJIS0szMVnXVlFRwe7du4mLizM7SpeVlJREbGxsvb+NsrIy1qxZo78Nkx04cIDi4mL9fbQiwzC44447eOutt/jwww9JSkqq9/Vx48YRFBRU7+9hx44dFBQU6O+hlektodPIyspi9uzZjB8/nokTJ7Jo0SIqKyvJzMw0O1qXcc8993D55ZfTv39/Dh06xIIFCwgICGDmzJlmR+vUKioq6v1fen5+Phs3biQqKop+/fpx99138+ijjzJ48GCSkpKYN28e8fHxTJ8+3bzQndDpfg9RUVE8/PDDXHPNNcTGxrJ7927uu+8+Bg0aREZGhompO5fbb7+dV199lX/961+Eh4d7j0uJjIzEbrcTGRnJzTffTFZWFlFRUURERHDnnXeSlpbGpEmTTE7fyZh9mpK/+93vfmf069fPCA4ONiZOnGisXr3a7EhdyowZM4y4uDgjODjY6NOnjzFjxgxj165dZsfq9FauXGkA37nNnj3bMAzPqc3z5s0zYmJijJCQEOOiiy4yduzYYW7oTuh0v4cTJ04YU6ZMMXr16mUEBQUZ/fv3N2655RajsLDQ7NidSkM/f8D4y1/+4h3jcDiM2267zejRo4cRGhpqXHXVVcbhw4fNC91JWQzDMNq/JomIiIg0nY5hEREREb+nwiIiIiJ+T4VFRERE/J4Ki4iIiPg9FRYRERHxeyosIiIi4vdUWERERMTvqbCIiIiI31NhEREREb+nwiIiIiJ+T4VFRERE/J4Ki4iIiPi9/w8wwr7WoZPKzQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# phi(X) plotted with fit, showing very small error\n",
    "X = np.linspace(start=0,stop=L,num=51)\n",
    "F = [fit(x) for x in X]\n",
    "plt.plot(X,F,\"r\",X,phi(X),\"b--\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a1040fc4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.2499999999999986 0.12499999999999258 0.0 -0.25000000000001144\n",
      "0.4999999999999731 0.0 -1.0000000000000513 -0.5000000000000782\n"
     ]
    }
   ],
   "source": [
    "# numerical overlap computation\n",
    "overlap_estimate = np.dot(c, np.dot(overlap_matrix, c))\n",
    "ke_integral = np.dot(c, np.dot(ke_matrix, c))\n",
    "angular_ke_integral = np.dot(c, np.dot(angular_ke_matrix, c))\n",
    "pe_integral = np.dot(c, np.dot(pe_matrix, c))\n",
    "\n",
    "ke_approx = ke_integral / overlap_estimate\n",
    "pe_approx = pe_integral / overlap_estimate\n",
    "angular_ke_approx = angular_ke_integral / overlap_estimate\n",
    "\n",
    "energy_approx = ke_approx + angular_ke_approx + pe_approx\n",
    "print(overlap_estimate, ke_integral, angular_ke_integral, pe_integral)\n",
    "print(ke_approx, angular_ke_approx, pe_approx, energy_approx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "09bbea5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def energy(c):\n",
    "    overlap_integral = np.dot(c, np.dot(overlap_matrix, c))\n",
    "    ke_integral = np.dot(c, np.dot(ke_matrix, c))\n",
    "    pe_integral = np.dot(c, np.dot(pe_matrix, c))\n",
    "#    angular_ke_integral = np.dot(c, np.dot(angular_ke_matrix, c))\n",
    "    return (ke_integral + pe_integral) / overlap_integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "df0f7479",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient(c):\n",
    "    h_matrix = ke_matrix + pe_matrix # + angular_ke_matrix\n",
    "    integral = np.dot(h_matrix - energy(c) * overlap_matrix, c)\n",
    "    overlap_integral = np.dot(c, np.dot(overlap_matrix, c))\n",
    "#    print(energy(c) * overlap_matrix, overlap_integral, integral, h_matrix)\n",
    "    return 2 * integral / overlap_integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "21e4141a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.87806476e+01, -3.41475424e-01, -1.02248964e-01, -5.14551408e-02,\n",
       "       -1.91514721e-02,  4.32656753e-02,  1.29072282e-01,  2.35982122e-01,\n",
       "        3.63160039e-01,  5.10162899e-01,  6.76723277e-01,  8.62666998e-01,\n",
       "        1.06787398e+00,  1.29225750e+00,  1.53575239e+00,  1.79830797e+00,\n",
       "        2.07988363e+00,  2.38044598e+00,  2.69996693e+00,  3.03842237e+00,\n",
       "        3.39579129e+00,  3.77205503e+00,  4.16719688e+00,  4.58120225e+00,\n",
       "        5.01405792e+00,  5.46574218e+00,  5.93624286e+00,  6.42581753e+00,\n",
       "        6.93483429e+00,  7.46690521e+00,  8.02722967e+00,  8.65385072e+00,\n",
       "        9.35708987e+00,  1.02053260e+01,  1.12273142e+01,  1.24711509e+01,\n",
       "        1.41645041e+01,  1.59441173e+01,  1.88479265e+01,  2.15479047e+01,\n",
       "        2.73277810e+01,  3.13736469e+01,  4.37506719e+01,  5.08106882e+01,\n",
       "        8.57022232e+01,  9.82902826e+01,  2.34646669e+02,  2.70343178e+02,\n",
       "        2.28077232e+03,  2.47497011e+03])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# this finds all eigenstates possible... first one is ground state\n",
    "scipy.linalg.eigh(ke_matrix + pe_matrix, overlap_matrix, eigvals_only=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "90d98e2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-4.00000001,  3.99999999, -4.00000001,  3.99999999, -4.00000001,\n",
       "        3.99999999, -4.00000001,  3.99999999, -4.00000001,  3.99999999,\n",
       "       -4.00000001,  3.99999999, -4.00000001,  3.99999999, -4.00000001,\n",
       "        3.99999999, -4.00000001,  3.99999999, -4.00000001,  3.99999999,\n",
       "       -4.00000001,  3.99999999, -4.00000001,  3.99999999, -4.00000001,\n",
       "        3.99999999, -4.00000001,  3.99999999, -4.00000001,  3.99999999,\n",
       "       -4.00000001,  3.99999999, -4.00000001,  3.99999999, -4.00000001,\n",
       "        3.99999999, -4.00000001,  3.99999999, -4.00000001,  3.99999999,\n",
       "       -4.00000001,  3.99999999, -4.00000001,  3.99999999, -4.00000001,\n",
       "        3.99999999, -4.00000001,  3.99999999, -4.00000001,  3.99999999])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gradient(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dda33758",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     fun: -18.780647545715542\n",
      "     jac: array([ 8.35419399e-06, -3.30017360e-06,  4.37349622e-07,  1.77614019e-06,\n",
      "        1.32416608e-06, -3.15098753e-06,  2.44236515e-06,  3.31868423e-06,\n",
      "       -3.88975596e-06, -7.58150173e-07,  1.70489338e-06, -3.72398065e-07,\n",
      "       -1.93723074e-06,  1.92673290e-06, -2.71227875e-06,  8.27662229e-08,\n",
      "       -5.51671292e-07,  6.09276595e-07,  1.16130635e-06, -1.72882288e-06,\n",
      "       -1.24296784e-07, -6.13646622e-07, -8.87545411e-07,  1.00212487e-07,\n",
      "       -2.21195074e-06,  5.02587398e-07, -3.22059264e-06,  1.06256780e-06,\n",
      "       -4.34988609e-06, -1.81294921e-07, -3.08877318e-06,  1.14450755e-06,\n",
      "       -4.01077733e-06, -9.85629571e-07, -2.21890494e-06, -6.68939471e-08,\n",
      "       -3.24653304e-06,  5.73144241e-07, -3.60616169e-06,  1.39805396e-06,\n",
      "       -8.40929026e-06,  2.30050795e-06, -2.16717844e-06,  8.12820306e-07,\n",
      "       -8.12766603e-06,  1.25587569e-06, -3.20283838e-06, -3.25177096e-07,\n",
      "       -1.99043968e-06, -2.43183438e-06])\n",
      " message: 'Optimization terminated successfully.'\n",
      "    nfev: 136\n",
      "     nit: 45\n",
      "    njev: 136\n",
      "  status: 0\n",
      " success: True\n",
      "       x: array([  2.2665747 ,  -6.69282839,  10.80676772, -14.42629573,\n",
      "        17.40592741, -19.64490436,  21.09290718, -21.74973847,\n",
      "        21.6631395 , -20.91815545,  19.62665682, -17.91833138,\n",
      "        15.92878207, -13.7887467 ,  11.61534826,  -9.5076213 ,\n",
      "         7.54093922,  -5.76898316,   4.22451181,  -2.91829207,\n",
      "         1.84765784,  -0.99804661,   0.34662708,   0.13397438,\n",
      "        -0.47195821,   0.69519944,  -0.82926855,   0.89633685,\n",
      "        -0.91522716,   0.89997757,  -0.86412109,   0.81596066,\n",
      "        -0.76103067,   0.70331859,  -0.64779012,   0.59495228,\n",
      "        -0.54650079,   0.50113069,  -0.46257908,   0.42555413,\n",
      "        -0.40294872,   0.37054273,  -0.37110127,   0.34514928,\n",
      "        -0.42445366,   0.39496922,  -0.72464015,   0.68588032,\n",
      "        -1.94892585,   1.87363183])\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(1000)\n",
    "c0 = np.random.random(n)\n",
    "result = scipy.optimize.minimize(energy, c0, jac=gradient, method=\"CG\")\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "6da2773f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nnp.random.seed(2023)\\nc0 = np.random.random(n)\\nconjugate_gradient(energy, gradient, c0)\\n'"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def conjugate_gradient(function, gradient, c0):\n",
    "    p = 1\n",
    "    epsilon = 1e-3 # threshold for error\n",
    "    position = c0\n",
    "    r = -gradient(position)\n",
    "    d = r\n",
    "    old_r = r\n",
    "\n",
    "    # calculating step size\n",
    "    def diff(f, x, h=1e-3):\n",
    "        return (f(x+h) - f(x-h)) / (2*h)\n",
    "    \n",
    "    def diff_2nd(f, x, h=1e-3):\n",
    "        return (diff(f, x+h) - diff(f, x-h)) / (2*h)\n",
    "    \n",
    "    while np.linalg.norm(r) > epsilon:\n",
    "        line_function = lambda s: function(position + s * d)\n",
    "        diff_line_function = diff(line_function, 0)\n",
    "        diff_2nd_line_function = diff_2nd(line_function, 0)\n",
    "        step_size = - diff_line_function / diff_2nd_line_function\n",
    "\n",
    "        position += step_size * d\n",
    "        old_r = r\n",
    "        r = -gradient(position)\n",
    "        beta = np.dot(r, r - old_r)/np.dot(old_r, old_r)\n",
    "        d = r + beta * d\n",
    "#        print(p, function(position), np.linalg.norm(r), step_size, diff_line_function, diff_2nd_line_function)\n",
    "        p += 1\n",
    "        \n",
    "    return position, function(position), gradient(position)\n",
    "\n",
    "'''\n",
    "np.random.seed(2023)\n",
    "c0 = np.random.random(n)\n",
    "conjugate_gradient(energy, gradient, c0)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "13f817a9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7346.952545953678 74145.58529088588 0.0 -212126.11159153937\n",
      "10.092019082350197 0.0 -28.872666628065737 -18.78064754571554\n"
     ]
    }
   ],
   "source": [
    "cc = result.x\n",
    "\n",
    "overlap_estimate = np.dot(cc, np.dot(overlap_matrix, cc))\n",
    "ke_integral = np.dot(cc, np.dot(ke_matrix, cc))\n",
    "angular_ke_integral = np.dot(cc, np.dot(angular_ke_matrix, cc))\n",
    "pe_integral = np.dot(cc, np.dot(pe_matrix, cc))\n",
    "\n",
    "ke_approx = ke_integral / overlap_estimate\n",
    "pe_approx = pe_integral / overlap_estimate\n",
    "angular_ke_approx = angular_ke_integral / overlap_estimate\n",
    "\n",
    "energy_approx = ke_approx + angular_ke_approx + pe_approx\n",
    "print(overlap_estimate, ke_integral, angular_ke_integral, pe_integral)\n",
    "print(ke_approx, angular_ke_approx, pe_approx, energy_approx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3c9ee827",
   "metadata": {},
   "outputs": [],
   "source": [
    "# integrated function that does all the calculations together\n",
    "def setup(n, Z=1, epsilon=1e-6):\n",
    "\n",
    "    roots, weights = roots_legendre(n)\n",
    "    P = np.array([eval_legendre(i, roots) for i in range(n + 1)])\n",
    "    dP = np.zeros((n+1, n))\n",
    "    dP[1] = np.ones(n)\n",
    "    for i in range(2,n+1):\n",
    "        dP[i] = (2*i-1) * P[i-1] + dP[i-2]\n",
    "    \n",
    "    overlap_matrix = np.zeros((n,n))\n",
    "    for l in range (n):\n",
    "        for m in range (n):\n",
    "            approx_overlap = np.dot(P[l]*P[m], weights) * L * 0.5\n",
    "            if abs(approx_overlap) < epsilon:\n",
    "                approx_overlap = 0\n",
    "            overlap_matrix[l, m] = approx_overlap\n",
    "\n",
    "    pe_matrix = np.zeros((n, n))\n",
    "    for l in range (n):\n",
    "        for m in range (n):\n",
    "            approx_pe = -Z * np.dot(P[l] * P[m] / (roots + 1), weights)\n",
    "            if abs(approx_pe) < epsilon:\n",
    "                approx_pe = 0\n",
    "            pe_matrix[l, m] = approx_pe\n",
    "\n",
    "    ke_matrix = np.zeros((n, n))\n",
    "    for l in range (n):\n",
    "        for m in range (n):\n",
    "            approx_ke = np.dot(dP[l] * dP[m], weights) / L\n",
    "            if abs(approx_ke) < epsilon:\n",
    "                approx_ke = 0\n",
    "            ke_matrix[l, m] = approx_ke\n",
    "    \n",
    "    # from exact phi\n",
    "    phi = lambda r: r * (math.e)**(-Z * r)\n",
    "    c_exact = np.zeros(n)\n",
    "    for l in range(n):\n",
    "        c_exact[l] = ((2*l+1)/2) * np.dot(phi(L * 0.5 * (roots + 1)) * P[l], weights)\n",
    "    \n",
    "    def energy(c):\n",
    "        overlap_integral = np.dot(c, np.dot(overlap_matrix, c))\n",
    "        ke_integral = np.dot(c, np.dot(ke_matrix, c))\n",
    "        pe_integral = np.dot(c, np.dot(pe_matrix, c))\n",
    "        return (ke_integral + pe_integral) / overlap_integral\n",
    "\n",
    "    def gradient(c):\n",
    "        h_matrix = ke_matrix + pe_matrix\n",
    "        integral = np.dot(h_matrix - energy(c) * overlap_matrix, c)\n",
    "        overlap_integral = np.dot(c, np.dot(overlap_matrix, c))\n",
    "        return 2 * integral / overlap_integral\n",
    "\n",
    "    np.random.seed(1000)\n",
    "    c0 = np.random.random(n)\n",
    "    result = scipy.optimize.minimize(energy, c0, jac=gradient, method=\"CG\")\n",
    "\n",
    "    return energy(c_exact), result.fun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "729a266a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.5000000000000783, -18.780647545715542)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "setup(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "881ede9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAckAAAEmCAYAAAAJLWsmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABA0UlEQVR4nO3dd3hTZcMG8PskadKd7qalm5YWKFM2KkUQUGQICiqgRUQZIkNlvH6or4psFBeKFgFREZWhICiyRQQFCyJQ2tJJ995pk5zvj0JeCqS00PYk7f27rlxNTk6S21p69znjOYIoiiKIiIjoBjKpAxAREZkrliQREZEJLEkiIiITWJJEREQmsCSJiIhMYEkSERGZwJIkIiIygSVJRERkgkLqAE3JYDAgLS0NDg4OEARB6jhERCQBURRRXFwMb29vyGS1jxVbVEmmpaXB19dX6hhERGQGUlJS4OPjU+s6LaokHRwcAFR/YxwdHSVOQ0REUigqKoKvr6+xE2pjcSX54YcfYvny5cjIyECnTp3w/vvvo0ePHnV67dVNrI6OjixJIqIWri673SzqwJ1vvvkGc+bMwWuvvYZTp06hU6dOGDx4MLKysqSORkREzZBFleSqVaswefJkTJw4Ee3atcPHH38MW1tbrFu3TupoRETUDFlMSVZWVuLkyZMYOHCgcZlMJsPAgQNx7NgxCZMREVFzZTH7JHNycqDX6+Hp6VljuaenJy5cuHDT12i1Wmi1WuPjoqKiRs1IRETNi8WMJG/H4sWLoVarjTee/kFERPVhMSXp5uYGuVyOzMzMGsszMzOh0Whu+poFCxagsLDQeEtJSWmKqERE1ExYTEkqlUrcdddd2Ldvn3GZwWDAvn370Lt375u+RqVSGU/3aKjTPs7lxmPZqfewPf4XZJbm3PH7ERGR+bKYfZIAMGfOHDz11FPo1q0bevTogXfffRelpaWYOHFik2U4lxeDHBxBYelB/BkPlGttIdd7wEMVgDbqdujl2QWedq5NloeIiBqPRZXk2LFjkZ2djVdffRUZGRno3Lkz9uzZc8PBPI3pkZAHMdIwGGdyYvB3zmkkVcSiBJeRgyMoKDmIEyVAeYUtFAZPuKsCEapui96aLnC3dWmyjERE1DAEURRFqUM0laKiIqjVahQWFjb4jDs6gw6nsy/g75x/kFx2ESX6y7BS5UEhNwAAyivsoDB4wkMVgDB1O/TUdIG7rXODZiAiolurTxewJBtRlb4Kp3Mu4O+cM0gpi0OJIRVWyvwbitNTFYQwp7boqekCNxunRs9FRNSSsSRNaOqSvJlKfRWis8/jdO4ZJJfGolRMg1KZD7ncAFEEKrT2sBb90NGpJx4IuA9qlZ0kOYmImiuWpAnmUJI3U6mvwt/Z/+J07j9ILo1FuXAJ1qoy6PQy6Cq84G/bGUN8B6KNi7/UUYmILB5L0gRzLcnriaKIv7PP4UDaAaRp/4HKOhuCAJSWqeGiaIs+HvegX6vusJJb1HFXRERmgSVpgqWU5PUyy3KxJ+lXnC/6CwarJFgpdKioVEKhC0K4ugceDLgPrjaW899DRCQllqQJllqS16rUV+JA6u84nn0UhWIMrFUl0Otl0FZo4GvdEYN8ByDcrbXUMYmIzBZL0oTmUJLX+zf3IvZe3ofU8jOwss6ETABKyxyhloWhl/vduM+3J1QKK6ljEhGZDZakCc2xJK+VW5GPPUn78G/hCegUicbNsrIqf7R17IGhAf3hacdJDYioZWNJmtDcS/JaOoMOh9P+wLHMI8g3XIBKVQy9QUBFmQZhDn0xJngYXGwcpI5JRNTkWJImtKSSvF5swSX8kvIrEstPQmmdjcoqBYTKENzr+QCG+PeGQi6XOiIRUZNgSZrQkkvyWhfyYrEjaRtyDH9DaaVFSZk9POTdMCpoJNq6+kkdj4ioUbEkTWBJ1qQz6LA35QB+y/oZBqsEiADKSr3RQR2BMcEPQG1tK3VEIqIGx5I0gSVpWl5FLr67tB2xpUehVBWgXKuEoioMA72HYqB/N8gEQeqIREQNgiVpAkvy1kRRxJncs9iVsgOF4j9QKKpQWOKEVlY98EjrEQhx8ZI6IhGZMVEUoRdFGEQRoihChAhRRPVjVH8Fqp+vXgcwQASM64jGdUXx2teKV14LuNnawVGluu2MLEkTWJL1U6nXYnfyzzieuw+wSoHeIEN5iR+6uPTHo8GDYH8HP6REVD+iKKJSr4dWr0OF7sabVqe76XPXL7v6WKc3QC8aoDOIMIgG6AwG6A0idKIBesPVxwboRbH68dWvhurXVd8XjfevvtbQBJWyatADGBnW7rZfz5I0gSV5+zLK0rE1YRsSK/6AlVUJSsutodKHY0irh9DPrwM3xxLdgkEUUaStQH5FBQorKpBfUX7lawUKKspRUFFx5VZ9v0irhVanQ4X+fyVYn1/WCpkM1nIFVAoFrBUKqBRyWMuv3reClUwGuUwGhUyAXJBBYXwsg1wQIL/y9epyuUwGhSCDXHZl2ZX7N3ut7JobIEAmAIIgQAYBwtX7ggABuPJVgCBUPye75r6A69a7sjzQyQVutrd/zARL0gSW5J0TRRF/Zf+Fny//iBKch0ymR0GxG/yte2Fs8HAEOLlJHZGo0ekNBmSUliCvvBwF5eXIv7bktBUoKC+/8vXK4yuFeLNfttYKBZxU1nCysan+al19c1SpjAVnrbCq/npN6V1ffNYKqxrPKWSyJv++WAqWpAksyYZVoSvDj0m78Hf+AQjKDFRWyVFZFoJhvo9igH9nCBxdkgUrraxEclEhUgoLkFxYiOTCAqQUFSKpsBCXiwpRZTDUWN9KJoPa2hrO1jbGonO62f0rheh8Zbk1p41scixJE1iSjSelJAlbE75HWtVxyORVKCxshT5uwzA2dACsOFEBmSGDKCKrtKRGAV69n1xYiNzyMuO6NgoF/NRO8FOr4et45ataDXdbOzhb20BtbQ07Kyv+YWghWJImsCQbX4WuDN8nbMWZol+gsCpBQbEzQmwGYGK7h6G2tpE6HrUwFboqpBQWVRfflVFhUmH115SiImj1OuO6HnZ21UXoqP5fIaqr77vZ2LIAmxGWpAksyaZjEPX49fKvOJC5HTKrTJSU28BZ7IWJYY/Dn/stqRFU6fWIyc1BdEZ69S0zHZfy843PK+Vy+Dmq4at2gr+6+uvVQvRxdISNFTd7thQsSRNYktKIzv0b25O+RqU8DtoqBWTaDhgT9Bi6e4VIHY0sWHpxMaIz042l+E9WJip0OihkMoS5uaOLxgvt3T0Q4OQMP7UaHnb2PAqbALAkTWJJSiu1JBlfx29CruFvQDCgrDgQA71G4aHWvfnLi2pVVlWFs1mZiM5Ix98Z6TidkY6M0hIAgJe9A7povND5yi3cw4MHw1CtWJImsCTNQ3FlIb6O24z4ikNQKCpQUOSBruoHMC7sQVhzk1eLZxBFJOTnITozw1iIF3KyoRdF2CgU6OipQSeNV3UxenrB095e6shkYViSJrAkzUuVvhI/JP2IE3m7oFAWoKjUAb5W92Jiu0fhbsv/Py1Ffnk5Tmdm4O+MNJzOyEB0ZjqKtFoAQLCzCzp7VZdhF40XQlzdeP4f3TGWpAksSfMkiiJ+y/gNey5/ByhTUKZVwU7XFRNCn0CoSyup41EDq9LrcTI9DfsT4nEgMQHx+XkAAGdra3TWeKOzRoPOGi908tTAUWUtcVpqjliSJrAkzV9MwQVsSfgSZcI56A0y6MpCMcJ/DO716cBD8C1YQUU5DiUlYl9CPA4nJaJIq4W7rR36BwSil48fumi84KdW8/8xNQmWpAksScuRU56FTbGbkK6rnpyguMgX93qMwOiQCMi5uc3siaKIS/l52JdwCfsTLuFk+mXoRRHt3T1wX2AQBgS2RriHJw/YIkmwJE1gSVqecl0ZtsR/h39L9sLKqhR5hR7o4TQST7QdCJVCIXU8ukalXo8/01Kx/0oxJhUWQCVXoK+vHwYEtUb/gEBo7B2kjknEkjSFJWm5DKIeO5N+wu+530NuVYTsfA90chiGyPD7YadUSh2vxcorL8OhxCubUZMTUVJZCY2dPfoHBuG+wCD08fHjSfpkdliSJrAkLZ9e1GNv6k84nP0dBEUx0nM90Nb2AUzqeD9cbG7/0jlUN6IoIjYvF/sS4rE/4RL+zkiHQRTR0VOD+wKCMCAwCO3cPbhvkcxasyzJRYsWYdeuXYiOjoZSqURBQUG934Ml2XzoRR32p+/GwczvIMqKkZrtgSDrQXi200B4O/D/bUPS6nQ4cTkV+xOrN6OmFBXCRqHA3X7+uC+wejOqhx3PVSTL0SxL8rXXXoOTkxNSU1MRFRXFkiQAgM5QhSNZe7A/41vohRIkZXmglaI/pnS+H0HOLlLHs1illZU4mJiAPfEXcTAxAaVVVfCyd8CAwCDcF9gavX18uU+YLFazLMmr1q9fj1mzZrEkqYYqQyWOZu/B/ozvUIVSJGa4w1W8B892GYCOnhqp41mEIq0W+xPisTvuIg4nJUGr16G9uwcGtw7BgKDWCHN142ZUahbq0wXN+k9BrVYL7ZWZO4Dqbww1T1YyJSI8h6Ov+2D8nr0H+4XvUCluxaK/f4NNVU882/k+9Pbx5S/56+SVl2HvpXjsiYvF7ylJqDIY0NnTC3N698Hg1iHwUztJHZFIUs26JBcvXoz//ve/UsegJmQlU6Gf5wj0dh+MY9m7sV/4HpWGnXjn7Am8c6IrnuncD/cHBbfo8/OyS0vxc3ws9sTH4nhqCgyiiO7ePlhwdz8Mah3MfbpE15B0c+v8+fOxdOnSWtc5f/48wsLCjI/rs7n1ZiNJX19fbm5tQbT6chzL+QkHMrei0qBFbJonygrD8XSnezAitC2s5HKpIzaJtOIi/Bwfh91xF3Ey7TJkgoBePr4YEtwGg4KC4W5nJ3VEoiZjMfsks7OzkZubW+s6QUFBUF5zHhz3SdLtqNCX4vecXTictQOV+kpcSPVEXl4YIjv2wZj2HWDbDM/lSyoowJ74i9gTF4vTmRlQyuTo6+ePIcEhGBjYGs42NlJHJJKExeyTdHd3h7u7u5QRqIWwltvhPs8x6O36II7m/Ail/EdU+hzGD6kX8cmpIDwR3hMjQtvC19Gy5w+Ny8vFnrhY7Im7iHM52bBWKHCvfwBWdXoQ9wUGwVGlkjoikUWxmH2SycnJyMvLQ3JyMvR6PaKjowEAwcHBsOf15KiObBT2GKh5HH3cHsKR7B2wku9CW59MHEhNxAd/ecFB4WC8VmEnjcbsr0RRXlWFi7k52J94CbtjLyIuPw92VlboHxiEad17ISIgsFmOkomaisWcAhIZGYkNGzbcsPzAgQOIiIio03twcytdr0RXgCNZO/BH7m4IUECp7YWYVA9EZ2ShuLJ6f3ZrZxfjVe87e2oQ6ube5Nc0rNLrkVCQj4u5ObiYm4uLuTmIyc1BcmEBRACOKhUGBrbGkOAQ3OMXwHMYiWphMfskmxpLkkwpqsrDrxlf41T+friqvDFIMwEqXRBOZ2UiOiMd0RnpuJCTDb0owkahQLiH55VrHlaPOr0cGmbiboMoIrWo8EoJ/q8ME/LzUGUwAAA87OzQxsUNbVzd0MbVFW1c3dDO3QPKFnIQEtGdYkmawJKkW0kvT8Du9A2ILzmDILtwPOAdCW+bIADVmzbPZmfidEYGojPS8XdGOtJLigEAnnb26KTRoIvGC509vdDBU1PrZk5RFJFdVoqY3BzE5OTgYl71CDE2NwflOh2A6tFhdRG6oY2LK0Kv3OcBN0R3hiVpAkuS6kIURcQUn8Se9I3I0V5GZ+cI3K95Amor1xvWzSotuTLSrC7OM1kZKKuqgkwQ0MbVrXrfpqcGPo5qXMrP+9/m0rwcFFRUAACsFQqEuFSPCK8WYRtXV3ja2Vv0QURE5oolaQJLkupDL+rxV95e/JqxGVWGCtztPhL3uI+ASm56JKc3GBCbl4vTV0aa0ZkZiM3NgQhAIZMhyMm5ugzd3IybTH0cHXkhaaImxJI0gSVJt6NCX4pDWVvxe85O2Mirj47t6twfMqFu+wBLKiuRWVIMX7UT9xsSmQGWpAksSboT+ZVZ+CXjS5wpOAKNtT8e8IpEsEMnqWMRUT2xJE1gSVJDSCm7iJ/SPkdyWQzaOHTFA15PwcPaV+pYRFRHLEkTWJLUUERRxL+Fx/BzxhcoqMxGN9f7McBzLOwVTlJHI6JbsJhp6YgslSAICHfqgzDH7vgj9yccyPoOp/MPo5/HaPRxewhWMuWt34SIzB5HkkQNoExXjP2ZW3A8dw8crVwwSDMeHZ3u5ikcRGaIm1tNYElSY8vRpmFP+kacLzoBH5sQPOgdCX+7tlLHIqJrsCRNYElSU0ko+Rc/pX+OtPJLCFf3xiDNBLiqNFLHIiKwJE1iSVJTMogGnC44jL0ZX6JEV4h73Eein8doKGW8XBWRlHjgDpEZkAkydHGOQHt1bxzO2orD2dtwOv8whraahLaO3aWOR0R1wLmwiBqZUqbCQM3jeKHNu3BTeWNT4mJ8kfA28iozpY5GRLfAza1ETUgURZwrOo5daetQqitCP49RuMd9JE8ZIWpC3NxKZKYEQUB7dS+EOHTGgcxvcSDzW0TnH8JDrSahjUNXqeMR0XW4uZVIAkqZNQZ7TcCMNu9ArXTDhoS38FXiMhRUZksdjYiuwc2tRBITRRH/FB7FT2mfo0Jfhv6ej6Kv2zAoZKYv2kxEt4+bW4ksiCAI6Oh0N9o4dMX+zG/wa8ZX+Dv/AIa1ehat7TtIHY+oRePmViIzYS23xYPeEzE9ZCVs5Y5Yd+k1fJO0CkVVeVJHI2qxuLmVyAyJoojogoPYnb4RVQYtBng+ht5uD0IucOMP0Z3i5lYiCycIAro490eYYw/szfgKe9I34lT+fgxv9SwC7NpJHY+oxeDmViIzZiO3w/BWkzE1eCmUMhU+jf8/fJfyHkqqCqSORtQicHMrkYUwiAaczNuHXzI2wSDqMVDzBHq6DoZMkEsdjciicIJzE1iS1ByU6Yrxc8YmnMz7FRrrAAxv9Sz87EKljkVkMerTBdzcSmRhbBUOeNhnKp4LXgxBEPBJ/AJsT/0YFfpSqaMRNTscSRJZMIOox/Hcn/FLxiZYy20x3PtZtFX3kDoWkVnjSJKohZAJcvR2exAz26yGxjoAm5KWYHPSCh7YQ9RAWJJEzYCT0h1PBryCMb6zEF9yFu9efAF/5x9AC9pQRNQoWJJEzYQgCOjkfC9mha5GG4cu+C7lfaxPeBP5lVlSRyOyWCxJombGTqHGGL/ZmBDwH2RrU/HexVn4PWcnDKJe6mhEFsciSjIxMRGTJk1CYGAgbGxs0Lp1a7z22muorKyUOhqR2Qpz7IaZbVaji3N/7Epbh7XxryCzIlnqWEQWxSKmpbtw4QIMBgM++eQTBAcH4+zZs5g8eTJKS0uxYsUKqeMRmS2V3AbDW01GR6e7sS31I3wY+xL6eYxGP/dRvBQXUR1Y7Ckgy5cvx5o1a3Dp0qU6v4angFBLVmWoxMGs73A4axvcVN4Y5TsdvrZtpI5F1ORaxCkghYWFcHFxqXUdrVaLoqKiGjeilspKpsT9micwLWQZrGRKfBK3ALvS1qHSUCF1NCKzZZElGRcXh/fffx/PPfdcrestXrwYarXaePP19W2ihETmy8smEM8FL8EQryfxZ+4veC9mFmKLo6WORWSWJN3cOn/+fCxdurTWdc6fP4+wsDDj48uXL6Nfv36IiIjAZ599VutrtVottFqt8XFRURF8fX25uZXoilxtBranfoRLpWfRxbk/HvSKhK3CQepYRI3KYiY4z87ORm5ubq3rBAUFQalUAgDS0tIQERGBXr16Yf369ZDJ6jcQ5j5JohuJooiT+fuwO2095DIrDPOejHB1bwiCIHU0okZhMSVZH5cvX0b//v1x1113YdOmTZDL6395IJYkkWlFVXnYefkz/Fv0B9o69sCwVpOhtnKVOhZRg2t2JXn58mVERETA398fGzZsqFGQGo2mzu/DkiS6tbOFx/Dj5U9RZajEEK8n0c1lIGSCRR6+QHRT9ekCizhPcu/evYiLi0NcXBx8fHxqPGcBHU9kUcLVvdHargN2p2/Ajssf40zBbxjlOx0uSk+poxE1OYsYSTYUjiSJ6ieu+DS2pX6EMn0xHvSaiG4uA7mvkixeizhPkogaX7BDJ8xo8w46qPti++U12Ji4CEVVeVLHImoyLEkiqpW13BajfKdjQsB/kFZ+Ce9dnIXT+Ue4q4NaBJYkEdXJ1QnTQ+w7Y0vKO9icvBKlOs5iRc0bS5KI6sxW4YCx/nPwmN+LiC85g/cuzsL5oj+ljkXUaFiSRFRvHZz6Ymab1WhlE4xNiYvxfcr7qNCXSh2LqMGxJInotjhYOWNCwAI87DMd/xb+gfcuzkZ88RmpYxE1KJYkEd02QRDQzWUAZrR5B65KDdYlvI4fL3+KSoP21i8msgAsSSK6Y85KD0wMeh0PeU/Cybx9+ODiHCSXxkgdi+iOsSSJqEHIBBl6uw3F821WwVbhgLXxr+Dn9E3QGaqkjkZ021iSRNSg3FTemNx6Ee7XPI6jOT/go7i5SCtPkDoW0W1hSRJRg5MLcvTzGI1pwcsgAPg4bh4OZH4HvaiXOhpRvbAkiajRaGwCMDV4Ge5xH4l9mZuxNu4/yK5IlToWUZ2xJImoUSlkVrhf8wSeC34bFfpSfBD7En7P2QmDaJA6GtEtsSSJqEn42rbB9DYr0cP1fuxKW4d1l15HfmWW1LGIasWSJKImo5SpMNR7EiYF/Rf5lZl47+IsnMzbx8nSyWyxJImoyQXZd8CMNu8gXN0HW1M/xObklSjXlUgdi+gGLEkikoS13BajfZ/HY34vIa7kNN6PnYOEkn+ljkVUQ71K8sEHH0RhYaHx8ZIlS1BQUGB8nJubi3bt2jVYOCJq/jo49cGMkHfgovRE1KVX8UvGl9CLOqljEQEABLEeOwPkcjnS09Ph4eEBAHB0dER0dDSCgoIAAJmZmfD29oZeb57nQhUVFUGtVqOwsBCOjo5SxyGiaxhEPY5kb8evGZvhZROIsX6z4arykjoWNUP16YJ6jSSv71PubCeihiK7MgHB/04VeZEH9ZDkuE+SiMyKj20IpoesQAd1Xx7UQ5KrV0kKggBBEG5YRkTUkFRyG4zync6DekhyivqsLIoiIiMjoVKpAAAVFRWYMmUK7OzsAABaLa8hR0QNp4NTH/jatsF3KasRdelV3OsxCgM8x0Iu1OtXF9Ftq9eBOxMnTqzTep9//vltB2pMPHCHyDIZRD0OZ2/DvoxveFAP3bH6dEG9StLSsSSJLFtKWSy+TX4Xxbp8POQ9CV2d7+MuH6q3Rju6lYhISr5XDurhTD3UVFiSRGRRVHIbztRDTYYlSUQWqXqmnlVwVnpwph5qNCxJIrJYTkp3TAr6LwZqHseRrO1YG/cKcrXpUseiZoQlSUQWTSbIEeHxCJ4Nfhtl+mJ8EPsiTuXt50w91CBYkkTULPjahuD5kJUIV/fG96kf8KAeahAWU5LDhw+Hn58frK2t4eXlhQkTJiAtLU3qWERkRqoP6pmBx/xeRFxxNN6PnYOk0gtSxyILZjEl2b9/f2zZsgUxMTH4/vvvER8fj0ceeUTqWERkhjo49cWMNu/AycoNn8X/Hw5nbYVBNEgdiyyQxU4m8MMPP2DkyJHQarWwsrKq02s4mQBRy6IXddiXsRmHsrcixL4zHvF7AfYKJ6ljkcSa/WQCeXl5+PLLL9GnT59aC1Kr1aKoqKjGjYhaDrmgwCCv8YgMXIi08kv44OKLuFRyVupYZEEsqiTnzZsHOzs7uLq6Ijk5GTt27Kh1/cWLF0OtVhtvvr6+TZSUiMxJiEMXPN9mFdxVrbDu0uvYl/kNDKJ5XhyezIukm1vnz5+PpUuX1rrO+fPnERYWBgDIyclBXl4ekpKS8N///hdqtRo7d+40OXejVqutcWWSoqIi+Pr6cnMrUQtlEPU4mPU99mduQaBdOzzqNwuOVi5Sx6ImZjETnGdnZyM3N7fWdYKCgqBUKm9YnpqaCl9fX/z+++/o3bt3nT6P+ySJCAAulZzFluR3YBD1eNRvJkIcukgdiZpQfbpA0ouyubu7w93d/bZeazBUH6nGa1gSUX0F2Yfj+TYr8V3K+1if8CbudR+FgZrHeJ1KuoFF/EQcP34cf/75J+6++244OzsjPj4eCxcuROvWres8iiQiupa9wglPBryC37J3YG/Gl0gs/Rdj/V6Ek9JN6mhkRiziwB1bW1ts3boVAwYMQGhoKCZNmoSOHTvi0KFDUKlUUscjIgslE2S41+NhPNP6LRRV5eKD2Dk4X/Sn1LHIjFjseZK3g/skiciUMl0xvk/9ABeK/kRft2EYpBkPhaxu52CTZWn250kSETU0W4UDxvvPx1Dvp/FH7m6sjf8P8rQZUsciibEkiYiuEAQBfdwewrOt30a5vgQfxL6Efwp+lzoWSYglSUR0HR/bYEwPWYE2Dp2xOXkFdlz+BFWGSqljkQRYkkREN2Ett8NYvxcxotVzOJW3Hx/HzUd2xWWpY1ETY0kSEZkgCAJ6uA7G1OCl0IlV+CjuZUTnH5I6FjUhliQR0S1obAIwLXgZ2jn2xLcpq/F9ygeoNFRIHYuaAEuSiKgOVHIbPOo3E6N9ZuCfgqP4KHYuMiuSpY5FjYwlSURUD11d+mNayHLIBBnWxM7j5tdmjiVJRFRPHtY+mBK8FOFOvfFtympsT/2YR782UxYxdysRkblRylQY7TMD/rZtsTPtM1wuj8fj/i/BRekpdTRqQBxJEhHdJkEQ0N31fuPkAx/GvsS5X5sZliQR0R1qZdsa00NWINCuPTYlLsbP6V9AL+qljkUNgCVJRNQAbOR2GOc/D0M0T+K37B34/NLrKK7KlzoW3SGWJBFRAxEEAfd4jMTTQf9FjjYNH8S+iISSf6WORXeAJUlE1MAC7dtjesgKeKh8EHXpNRzO2gqDaJA6Ft0GliQRUSNwsHJGZNBruNfjYfycsQlfJi5Bua5E6lhUTyxJIqJGIhfkGKQZhwkB/0FS2QV8GPsSLpfFSx2L6oElSUTUyMIcu2F6yArYKhzwSfwCnMj9GaIoSh2L6oAlSUTUBJyVHni29dvo5jIQOy5/gu9S3uMk6RaAJUlE1EQUMisMb/UsxvjOwrmi41gTOw/ZFalSx6JasCSJiJpYJ+d7MTV4KUSI+ChuLs4U/CZ1JDKBJUlEJAEPa19MDV6KMMfu+CZ5FX68/Cl0hiqpY9F1OME5EZFEVHIbjPGdhQC7ttiVtg6pZXF4zP9FOCs9pI5GV3AkSUQkIUEQ0NN1CJ5tvQglugJ8GPsSYopOSR2LrmBJEhGZAR/bEEwPWQE/21BsTHwLezO+goGTpEuOJUlEZCZsFQ4YH7AA92uewKGsrdiQsAhlumKpY7VoLEkiIjMiE2SI8HgEkYELkVYez1l6JMaSJCIyQ8EOnTA9ZAXsFI5YG/8fnMzbJ3WkFoklSURkppyU7pjcehE6O/fD1tQPsT11DU8TaWI8BYSIyIxZyZR42GcafGxCsDPtM6SXJ+Jx/5fhpHSTOlqLYHEjSa1Wi86dO0MQBERHR0sdh4ioSXR3vR+TWy9CiS4fH8W+hPiSf6SO1CJYXEnOnTsX3t7eUscgImpyPrbBmBayAhqbAHx+6b84krWdVxNpZBZVkrt378Yvv/yCFStWSB2FiEgSdgpHRAYuxL3uI7EnYyO+TlqOCn2Z1LGaLYspyczMTEyePBlffPEFbG1t6/QarVaLoqKiGjciIksnE+QY5DUeT/jPRVzJaXwcNw9ZvJpIo7CIkhRFEZGRkZgyZQq6detW59ctXrwYarXaePP19W3ElERETau9uhemBi+DABnWxM3F2YLfpY7U7EhakvPnz4cgCLXeLly4gPfffx/FxcVYsGBBvd5/wYIFKCwsNN5SUlIa6b+EiEga7tatMCV4CUIduuLr5BXYk74Rek5n12AEUcK9vtnZ2cjNza11naCgIIwZMwY//vgjBEEwLtfr9ZDL5Rg3bhw2bNhQp88rKiqCWq1GYWEhHB0d7yg7EZE5EUURv+fsxJ70DQiwa4ex/nNgr3CSOpZZqk8XSFqSdZWcnFxjf2JaWhoGDx6M7777Dj179oSPj0+d3oclSUTNXULJv/g6eQUUggKP+78MX9s2UkcyO/XpAouYTMDPz6/GY3t7ewBA69at61yQ9aHX61FVxVkt6PZYWVlBLpdLHYNaqED79pgesgKbk1bg0/j/w0Pek9DdZVCNLXFUdxZRkk1FFEVkZGSgoKBA6ihk4ZycnKDRaPiLiSShtnLFpKA38FP6euy4/AlSyi5ieKtnYSVTSR3N4lhkSQYEBDTKCbRXC9LDwwO2trb8BUf1JooiysrKkJWVBQDw8vKSOBG1VAqZFYa3mgxf22BsT/0EGRVJeMJ/LpyVHlJHsygWWZKNQa/XGwvS1dVV6jhkwWxsbAAAWVlZ8PDw4KZXklQX5/7wtA7AV0nL8GHsSxjrNxshDl2kjmUxLOI8yaZwdR9kXScqIKrN1Z8j7tsmc+BtE4jpwcvha9sGGxLewoHMb2EQDVLHsggsyetwEys1BP4ckbmxUdhjQsB/0N9zDPZlbsaXSUs5nV0dsCSJiFoImSDDAM+xGB+wAAkl/2JN3Dxkczq7WrEkqU4CAgLw7rvv3tF7HDx4EIIgNNjRw4mJibxkGtFtCHPshmkhyyCDDGvi5uFc4QmpI5ktlmQzcezYMcjlcgwdOlTqKACAiIgIzJo1q8ayPn36ID09HWq1WppQRGTkpvLGlOAlCLbvhC+TluDXjM3cT3kTLMlmIioqCjNmzMDhw4eRlpYmdZybUiqVPHeQyIyo5DZ43P9l3K95AgezvsWmxCUo15dKHcussCSbgZKSEnzzzTeYOnUqhg4divXr1xufu7qJc9++fejWrRtsbW3Rp08fxMTEGNeJj4/HiBEj4OnpCXt7e3Tv3h2//vqryc97+umn8dBDD9VYVlVVBQ8PD0RFRSEyMhKHDh3C6tWrjRPVJyYm3nRz69GjRxEREQFbW1s4Oztj8ODByM/PBwDs2bMHd999N5ycnODq6oqHHnoI8fHxDfNNIyIA1QeZRXg8ggkBryCp7DzWxM5FVgUvBnEVz5OsRXlVFeLz85r8c1s7u8DGyqrO62/ZsgVhYWEIDQ3F+PHjMWvWLCxYsKDGiO2VV17BypUr4e7ujilTpuDpp5/G0aNHAVSX7IMPPohFixZBpVJh48aNGDZsGGJiYm6YEhAAnnnmGdx7771IT083niy/c+dOlJWVYezYsXjkkUdw8eJFhIeH44033gAAuLu7IzExscb7REdHY8CAAXj66aexevVqKBQKHDhwAHp99RUMSktLMWfOHHTs2BElJSV49dVX8fDDDyM6OhoyGf++I2pIoY5dMS14Gb5MWoo1cfPwiO8LaK/uJXUsybEkaxGfn4fhmzc1+ef+8Nh4hHt41nn9qKgojB8/HgAwZMgQFBYW4tChQ4iIiDCus2jRIvTr1w9A9SXKhg4dioqKClhbW6NTp07o1KmTcd0333wT27Ztww8//IDnn3/+hs/r06cPQkND8cUXX2Du3LkAgM8//xyPPvqocV5dpVIJW1tbaDQak7mXLVuGbt264aOPPjIua9++vfH+6NGja6y/bt06uLu749y5cwgPD6/rt4eI6shV5YXnWi/G1tQP8VXSMkR4PIIBnmMhE1ruhBgsyVq0dnbBD4+Nl+Rz6yomJgYnTpzAtm3bAAAKhQJjx45FVFRUjZLs2LGj8f7V0V9WVhb8/PxQUlKC119/Hbt27UJ6ejp0Oh3Ky8uRnJxs8nOfeeYZrF27FnPnzkVmZiZ2796N/fv31+u/Mzo6Go8++qjJ52NjY/Hqq6/i+PHjyMnJgcFQfVBBcnIyS5KokajkNnjM70UcyW6NXzK+RFr5JYzxnQUbhb3U0STBkqyFjZVVvUZ0UoiKioJOp4O3t7dxmSiKUKlU+OCDD4zLrK7ZfHt1M+zV0nnppZewd+9erFixAsHBwbCxscEjjzyCyspKk5/75JNPYv78+Th27Bh+//13BAYG4p577qlX9qvTt5kybNgw+Pv749NPP4W3tzcMBgPCw8NrzUVEd04QBNzr8TC8bALxTfIqfBQ3F+MD5sHT2l/qaE2OO3YsmE6nw8aNG7Fy5UpER0cbb6dPn4a3tze+/vrrOr3P0aNHERkZiYcffhgdOnSARqO5Yf/h9VxdXTFy5Eh8/vnnWL9+PSZOnFjjeaVSady3aErHjh2xb9++mz6Xm5uLmJgY/N///R8GDBiAtm3bGg/oIaKmEeLQGdOCl0EpU+HjuAX4p+B3qSM1OZakBdu5cyfy8/MxadIkhIeH17iNHj0aUVFRdXqfkJAQbN261ViwTzzxhHGUWZtnnnkGGzZswPnz5/HUU0/VeC4gIADHjx9HYmJijU2l11qwYAH+/PNPTJs2DWfOnMGFCxewZs0a5OTkwNnZGa6urli7di3i4uKwf/9+zJkzp27fGCJqMC4qDZ4LXowwh27YnLwCP6dvgkGs/Q/g5oQlacGioqIwcODAm56cP3r0aPz11184c+bMLd9n1apVcHZ2Rp8+fTBs2DAMHjwYXbt2veXrBg4cCC8vLwwePLjG5l6gehOuXC5Hu3bt4O7uftP9m23atMEvv/yC06dPo0ePHujduzd27NgBhUIBmUyGzZs34+TJkwgPD8fs2bOxfPnyW2YiooanlFljjN9sPOD1FI5kb8eGhEUo0xVLHatJCGJjXJjRTBUVFUGtVqOwsBCOjo41nquoqEBCQgICAwNhbW0tUULLUlJSglatWuHzzz/HqFGjpI5jVvjzRM1VXPFpfJO8Ciq5Lcb7z4PGJkDqSPVWWxdcjyNJqjeDwYCsrCy8+eabcHJywvDhw6WORERNJNihE6aFLIe1zBYfxy3AmYLfpI7UqFiSVG/Jycnw9PTEV199hXXr1kGh4EHSRC2Js9IDzwa/jXbqnvgmeRV2p22Avpnup+RvN6q3gIAAtKCt9ER0E0qZCo/6zkQrm9bYk74B6RUJeMzvRdgqHKSO1qA4kiQiotsiCAL6ug/DxKDXkF6eiI9iX0ZaeYLUsRoUS5KIiO5IkH0HTA9ZDhuFPdbGLUB0/iGpIzUYliQREd0xJ6U7nm29COFOffBtymr8lPZ5s9hPyX2SRETUIKxkKoz2mQFvmyDsTluPzIpkjPWbY9H7KTmSJCKiBiMIAvq4PYTIoNeQVn4Ja+LmIqM8SepYt40lSUREDa61fQdMDVkGpcwGn8QvwNnCY1JHui0sSZJUZGQkRo4cKXWMerHEzERScFF64rngtxHqcBe+TlqOvRlfwSDeel5oc8KSbCaOHTsGuVyOoUOHSh2lXlavXo3169c3+uew2IikoZRZY6zfHAzWjMehrO+xKXExKvSlUseqM5ZkMxEVFYUZM2bg8OHDSEtLa/TPa6hrOqrVajg5OTXIexGReaq+PuUoTAh4BUml57Embj6yKy5LHatOWJLNQElJCb755htMnToVQ4cOrTEyO3jwIARBwK5du9CxY0dYW1ujV69eOHv2rHGd9evXw8nJCdu3b0dISAisra0xePBgpKSkGNd5/fXX0blzZ3z22Wc1Ju1OTk7GiBEjYG9vD0dHR4wZMwaZmZkAgAsXLsDW1hZfffWV8X22bNkCGxsbnDt3DsCNI7yIiAjMmDEDs2bNgrOzMzw9PfHpp5+itLQUEydOhIODA4KDg7F7927ja/R6PSZNmoTAwEDY2NggNDQUq1evrpF9w4YN2LFjBwRBgCAIOHjwIAAgJSUFY8aMgZOTE1xcXDBixIga19LU6/WYM2cOnJyc4Orqirlz53K2IaLbFOrYFVNDlkGAgDVx83Ch6C+pI90SS7IWlQYtLpfFN/mt0qCtV84tW7YgLCwMoaGhGD9+PNatW3fDL/KXX34ZK1euxJ9//gl3d3cMGzYMVVVVxufLysqwaNEibNy4EUePHkVBQQEee+yxGu8RFxeH77//3njtSYPBgBEjRiAvLw+HDh3C3r17cenSJYwdOxYAEBYWhhUrVmDatGlITk5GamoqpkyZgqVLl6Jdu3Ym/3s2bNgANzc3nDhxAjNmzMDUqVPx6KOPok+fPjh16hQGDRqECRMmoKysDED1hOs+Pj749ttvce7cObz66qv4z3/+gy1btgCovmzXmDFjMGTIEKSnpyM9PR19+vRBVVUVBg8eDAcHBxw5cgRHjx6Fvb09hgwZYhwpr1y5EuvXr8e6devw22+/IS8vD9u2bavX/x8i+h83lTemBC9BkH04NiUuxoHMb836D09eKuuKm13a6HJZPD6Ke7nJc04LXo5Wtq3rvH7fvn0xZswYzJw5EzqdDl5eXvj2228RERGBgwcPon///ti8ebOxvPLy8uDj44P169djzJgxWL9+PSZOnIg//vgDPXv2BFA9Cmzbti2OHz+OHj164PXXX8fbb7+Ny5cvw93dHQCwd+9ePPDAA0hISICvry8A4Ny5c2jfvj1OnDiB7t27AwAeeughFBUVQalUQi6XY8+ePRAEAUD1SLKgoADbt28HUD2S1Ov1OHLkCIDqkZxarcaoUaOwceNGAEBGRga8vLxw7Ngx9OrV66bfk+effx4ZGRn47rvvbvo5ALBp0ya89dZbOH/+vDFPZWWlcVQ9aNAgeHt7Y/bs2Xj55eqfA51Oh8DAQNx111013ut6vFQWUe0MogEHsr7F/sxv0F7dC6N9ZkAlt2mSz67PpbIsZjKBgIAAJCXVPNdm8eLFmD9/fqN9pru1D6YFN/2Fft2tfeq8bkxMDE6cOGEc3SgUCowdOxZRUVGIiIgwrte7d2/jfRcXF4SGhuL8+fPGZQqFwlhqQPUo0MnJCefPn0ePHj0AAP7+/saCBIDz58/D19fXWJAA0K5dO+Prrr7funXr0KZNG8hkMvz777/GQjKlY8eOxvtyuRyurq7o0KGDcZmnpycAICsry7jsww8/xLp165CcnIzy8nJUVlaic+fOtX7O6dOnERcXBweHmic6V1RUID4+HoWFhUhPTzf+4XD1+9StWzez/suXyBLIBBkGeI6Fl3UAvk1ZjU/iF2Cc/3y4qjRSR6vBYkoSAN544w1MnjzZ+Pj6X24NTSlT1WtEJ4WoqCjodDp4e3sbl4miCJVKhQ8++KBBP8vOzu62Xnf69GmUlpZCJpMhPT0dXl5eta5vZWVV47EgCDWWXS1Zg6H6UPLNmzfjpZdewsqVK9G7d284ODhg+fLlOH78eK2fU1JSgrvuugtffvnlDc9d+8cAETWeduqemKJagi8Tl2JN3Fw85jcHwQ6dpY5lZFH7JB0cHKDRaIy32/2l3VzodDps3LgRK1euRHR0tPF2+vRpeHt74+uvvzau+8cffxjv5+fn4+LFi2jbtm2N9/rrr//tRI+JiUFBQUGNda7Xtm1bpKSk1DjA59y5cygoKDDuc8zLy0NkZCReeeUVREZGYty4cSgvL2+Q//6rjh49ij59+mDatGno0qULgoODER8fX2MdpVIJvb7mPJJdu3ZFbGwsPDw8EBwcXOOmVquhVqvh5eVVo2x1Oh1OnjzZoPmJWjpPaz9MDV4KH9sQrE94C79l7zCbrTUWVZJLliyBq6srunTpguXLl0On09W6vlarRVFRUY1bc7Jz507k5+dj0qRJCA8Pr3EbPXo0oqKijOu+8cYb2LdvH86ePYvIyEi4ubnVOKrUysoKM2bMwPHjx3Hy5ElERkaiV69exk2tNzNw4EB06NAB48aNw6lTp3DixAk8+eST6NevH7p16wYAmDJlCnx9ffF///d/WLVqFfR6PV566aUG/T6EhITgr7/+ws8//4yLFy9i4cKF+PPPP2usExAQgDNnziAmJgY5OTmoqqrCuHHj4ObmhhEjRuDIkSNISEjAwYMH8cILLyA1NRUAMHPmTCxZsgTbt2/HhQsXMG3aNBQUFDRofiICbBT2eDLgP7jHfQR2p2/AtymrUVXPgxgbg8WU5AsvvIDNmzfjwIEDeO655/D2229j7ty5tb5m8eLFxhGBWq2use+sOYiKisLAgQOhVqtveG706NH466+/cObMGQDVf2DMnDkTd911FzIyMvDjjz9CqVQa17e1tcW8efPwxBNPoG/fvrC3t8c333xT6+cLgoAdO3bA2dkZ9957LwYOHIigoCDj6zZu3IiffvoJX3zxBRQKBezs7LBp0yZ8+umnNU7huFPPPfccRo0ahbFjx6Jnz57Izc3FtGnTaqwzefJkhIaGolu3bnB3d8fRo0dha2uLw4cPw8/PD6NGjULbtm0xadIkVFRUGHfmv/jii5gwYQKeeuop46bchx9+uMGyE9H/yAQ5BntNwFi/OThX+AfWxr+CgspsSTNJenTr/PnzsXTp0lrXOX/+PMLCwm5Yvm7dOjz33HMoKSmBSqW66Wu1Wi202v/9JVJUVARfX986H93aHFw9ujU/P9/kSfvr16/HrFmzOEJqQM3154moqaSVJ+DLxCWoMmjxuP/LCLRv32DvbTFHt7744ouIjIysdZ2goKCbLu/Zsyd0Oh0SExMRGhp603VUKpXJAiUiIvPlbROIaSHLsTlpBdZdeh1DvZ9GT9chtzw6vqFJWpLu7u63fRRhdHQ0ZDIZPDw8GjgVERGZAzuFIyKDXsWe9A34Me1TpJVfwvBWz0Ihs7r1ixuIRZwCcuzYMRw/fhz9+/eHg4MDjh07htmzZ2P8+PFwdnaWOp5Zi4iIuOVRYpGRkbcc0RMRSUEuKDDUexK8rAOx4/InyNKm4An/uXC0cmmSz7eIklSpVNi8eTNef/11aLVaBAYGYvbs2ZgzZ47U0YiIqAl0dbkP7ta++CppKZLLYhCu7n3rFzUAiyjJrl271jjPrzGZy7k5ZNn4c0TU8HxtQzA79AMoZU13MJzFnALS2K7O6HJ10myiO3H15+j62YOI6M40ZUECFjKSbApyuRxOTk7G+UBtbW2b/CgqsnyiKKKsrAxZWVlwcnKCXC6XOhIR3QGW5DU0muqJda+dOJvodjg5ORl/nojIcrEkryEIAry8vODh4VHjWotE9WFlZcURJFEzwZK8Cblczl9yRETEA3eIiIhMYUkSERGZwJIkIiIyoUXtk7x6gndzu64kERHV3dUOqMukHy2qJIuLiwGg2V1XkoiI6q+4uPim1+O9lqTXk2xqBoMBaWlpcHBwMNuJAq5e8zIlJeWW1zkzN8wuDWaXBrNLoyGyi6KI4uJieHt7Qyarfa9jixpJymQy+Pj4SB2jThwdHS3uh/cqZpcGs0uD2aVxp9lvNYK8igfuEBERmcCSJCIiMoElaWZUKhVee+01qFQqqaPUG7NLg9mlwezSaOrsLerAHSIiovrgSJKIiMgEliQREZEJLEkiIiITWJJEREQmsCTNxOLFi9G9e3c4ODjAw8MDI0eORExMjNSx6m3JkiUQBAGzZs2SOkqdXL58GePHj4erqytsbGzQoUMH/PXXX1LHuiW9Xo+FCxciMDAQNjY2aN26Nd588806zUUphcOHD2PYsGHw9vaGIAjYvn17jedFUcSrr74KLy8v2NjYYODAgYiNjZUm7DVqy11VVYV58+ahQ4cOsLOzg7e3N5588kmkpaVJF/gat/qeX2vKlCkQBAHvvvtuk+WrTV2ynz9/HsOHD4darYadnR26d++O5OTkBs/CkjQThw4dwvTp0/HHH39g7969qKqqwqBBg1BaWip1tDr7888/8cknn6Bjx45SR6mT/Px89O3bF1ZWVti9ezfOnTuHlStXwtnZWepot7R06VKsWbMGH3zwAc6fP4+lS5di2bJleP/996WOdlOlpaXo1KkTPvzww5s+v2zZMrz33nv4+OOPcfz4cdjZ2WHw4MGoqKho4qQ11Za7rKwMp06dwsKFC3Hq1Cls3boVMTExGD58uARJb3Sr7/lV27Ztwx9//AFvb+8mSnZrt8oeHx+Pu+++G2FhYTh48CDOnDmDhQsXwtrauuHDiGSWsrKyRADioUOHpI5SJ8XFxWJISIi4d+9esV+/fuLMmTOljnRL8+bNE++++26pY9yWoUOHik8//XSNZaNGjRLHjRsnUaK6AyBu27bN+NhgMIgajUZcvny5cVlBQYGoUqnEr7/+WoKEN3d97ps5ceKECEBMSkpqmlB1ZCp7amqq2KpVK/Hs2bOiv7+/+M477zR5tlu5WfaxY8eK48ePb5LP50jSTBUWFgIAXFxcJE5SN9OnT8fQoUMxcOBAqaPU2Q8//IBu3brh0UcfhYeHB7p06YJPP/1U6lh10qdPH+zbtw8XL14EAJw+fRq//fYbHnjgAYmT1V9CQgIyMjJq/Oyo1Wr07NkTx44dkzBZ/RUWFkIQBDg5OUkd5ZYMBgMmTJiAl19+Ge3bt5c6Tp0ZDAbs2rULbdq0weDBg+Hh4YGePXvWujn5TrAkzZDBYMCsWbPQt29fhIeHSx3nljZv3oxTp05h8eLFUkepl0uXLmHNmjUICQnBzz//jKlTp+KFF17Ahg0bpI52S/Pnz8djjz2GsLAwWFlZoUuXLpg1axbGjRsndbR6y8jIAAB4enrWWO7p6Wl8zhJUVFRg3rx5ePzxxy1i0vClS5dCoVDghRdekDpKvWRlZaGkpARLlizBkCFD8Msvv+Dhhx/GqFGjcOjQoQb/vBZ1FRBLMX36dJw9exa//fab1FFuKSUlBTNnzsTevXsbZ39AIzIYDOjWrRvefvttAECXLl1w9uxZfPzxx3jqqackTle7LVu24Msvv8RXX32F9u3bIzo6GrNmzYK3t7fZZ2+OqqqqMGbMGIiiiDVr1kgd55ZOnjyJ1atX49SpU2Z72UBTDAYDAGDEiBGYPXs2AKBz5874/fff8fHHH6Nfv34N+nkcSZqZ559/Hjt37sSBAwcs4rJeJ0+eRFZWFrp27QqFQgGFQoFDhw7hvffeg0KhgF6vlzqiSV5eXmjXrl2NZW3btm2UI+Qa2ssvv2wcTXbo0AETJkzA7NmzLW40DwAajQYAkJmZWWN5Zmam8TlzdrUgk5KSsHfvXosYRR45cgRZWVnw8/Mz/rtNSkrCiy++iICAAKnj1crNzQ0KhaLJ/u1yJGkmRFHEjBkzsG3bNhw8eBCBgYFSR6qTAQMG4J9//qmxbOLEiQgLC8O8efMgl8slSnZrffv2veE0m4sXL8Lf31+iRHVXVlZ2w8Vi5XK58a9sSxIYGAiNRoN9+/ahc+fOAKovrHv8+HFMnTpV2nC3cLUgY2NjceDAAbi6ukodqU4mTJhww/EDgwcPxoQJEzBx4kSJUtWNUqlE9+7dm+zfLkvSTEyfPh1fffUVduzYAQcHB+O+GLVaDRsbG4nTmebg4HDDflM7Ozu4urqa/f7U2bNno0+fPnj77bcxZswYnDhxAmvXrsXatWuljnZLw4YNw6JFi+Dn54f27dvj77//xqpVq/D0009LHe2mSkpKEBcXZ3yckJCA6OhouLi4wM/PD7NmzcJbb72FkJAQBAYGYuHChfD29sbIkSOlC43ac3t5eeGRRx7BqVOnsHPnTuj1euO/WxcXFyiVSqliA7j19/z6QreysoJGo0FoaGhTR73BrbK//PLLGDt2LO699170798fe/bswY8//oiDBw82fJgmOYaWbgnATW+ff/651NHqzVJOARFFUfzxxx/F8PBwUaVSiWFhYeLatWuljlQnRUVF4syZM0U/Pz/R2tpaDAoKEl955RVRq9VKHe2mDhw4cNOf76eeekoUxerTQBYuXCh6enqKKpVKHDBggBgTEyNtaLH23AkJCSb/3R44cEDq6Lf8nl/PnE4BqUv2qKgoMTg4WLS2thY7deokbt++vVGy8FJZREREJvDAHSIiIhNYkkRERCawJImIiExgSRIREZnAkiQiIjKBJUlERGQCS5KIiMgEliQREZEJLEkiIiITWJJEREQmcIJzohYgIiICHTt2hLW1NT777DMolUpMmTIFr7/+utTRiMwaR5JELcSGDRtgZ2eH48ePY9myZXjjjTewd+9eqWMRmTVOcE7UAkRERECv1+PIkSPGZT169MB9992HJUuWSJiMyLxxJEnUQnTs2LHGYy8vL2RlZUmUhsgysCSJWggrK6sajwVBgMFgkCgNkWVgSRIREZnAkiQiIjKBJUlERGQCj24lIiIygSNJIiIiE1iSREREJrAkiYiITGBJEhERmcCSJCIiMoElSUREZAJLkoiIyASWJBERkQksSSIiIhNYkkRERCawJImIiExgSRIREZnw/1+viFcwxxoxAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = []\n",
    "e_exact = []\n",
    "e_approx = []\n",
    "\n",
    "for i in range(1,17,1):\n",
    "    exact, approx = setup(i)\n",
    "    X.append(i)\n",
    "    e_exact.append(exact)\n",
    "    e_approx.append(approx)\n",
    "\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(5, 3))\n",
    "ax1.set_xlabel(r'n')\n",
    "ax1.set_ylabel(r'E')\n",
    "\n",
    "colors = plt.cm.viridis([.5, .8])\n",
    "\n",
    "p1 = ax1.plot(X, e_exact, color=colors[0], ls='-', lw=1, label='Analytical')\n",
    "p2 = ax1.plot(X, e_approx, color=colors[1], ls='-', lw=1, label='Approximated')\n",
    "\n",
    "ax1.legend(handles=p1+p2, loc='best')\n",
    "\n",
    "plt.savefig('fig/hatom-broken.png', bbox_inches='tight', dpi=300)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}