{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "510a865d-ebbf-4fb6-ae57-dfde70702520",
   "metadata": {},
   "source": [
    "<center><img src=\"headers/header_arial.png\" alt=\"Drawing\" style=\"width: 1000px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f7cd575-b794-4a1a-ac9f-e84c79e9bda1",
   "metadata": {},
   "source": [
    "<div style=\"border-top: 3px solid #000; text-align: center; font-size: 24px; font-weight: bold; padding: 20px 0;\">\n",
    "    Problem set 5: Thomas-Fermi model for many-electron atoms \n",
    "</div>\n",
    "<div style=\"text-align: center; font-size: 15px; font-weight: normal; color: black; margin: 20px;\">\n",
    "    For questions contact: tabea.buhler@epfl.ch\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd9a84d-0803-4c61-8b25-68afe5475670",
   "metadata": {},
   "source": [
    "<div style=\"text-align: center; font-size: 22px; font-weight: bold; color: black; margin: 20px;\">\n",
    "    Exercise 1: From central field to the Thomas-Fermi equation for screening\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea2cfd9c-ee0e-4afa-9d58-90ab70fa1a4a",
   "metadata": {},
   "source": [
    "### Solution exercise 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52cda5b0-6954-4845-a7a8-957d286064dd",
   "metadata": {},
   "source": [
    "a) It must be $F \\rightarrow 1$ for $r \\rightarrow 0$, which corresponds to no screening near the nucleus."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3bb8ce7-a43a-4501-bac3-b22b1d355913",
   "metadata": {},
   "source": [
    "b) We write the expressio for $\\mu$ as the energy cost to add one particle to the system:\n",
    "$$\n",
    "\\mu = \\frac{p^2}{2m} + V(r) = \\frac{\\hbar^2k_F^2(r)}{2m} + V(r)\n",
    "$$\n",
    "Inserting $k_F(r) = (3 \\pi^2 n(r))^{1/3}$ leads to:\n",
    "$$\n",
    "\\mu = \\frac{\\hbar^2  (3 \\pi^2 n(r))^{2/3}}{2m} + V(r).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2dc916e-1a67-4987-8a46-1a578657e538",
   "metadata": {},
   "source": [
    "c) Assuming $\\mu = 0$ leads to:\n",
    "$$\n",
    "n(r) = \\frac{1}{3 \\pi ^2} \\bigl( \\frac{2m}{\\hbar^2 4 \\pi \\epsilon_0} \\frac{Z e^2}{r} F(r) \\bigr)^{3/2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9894ae6-4f39-497d-815e-8bc5f8a6cd07",
   "metadata": {},
   "source": [
    "d) From the Poisson equation for $\\phi(r)$ it follows: \n",
    "$$\n",
    "\\partial^2_r (r \\phi(r)) = \\frac{Z e}{4 \\pi \\epsilon_0} \\partial^2_r F(r) = - \\frac{\\rho(r)}{\\epsilon_0}.\n",
    "$$\n",
    "Using the relation $\\rho = -e n(r)$, we arrive at the desired result: \n",
    "\n",
    "$$\n",
    "\\frac{Z}{r}\\partial^2_r F(r) = 4 \\pi n(r).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3f98e70-fe89-4bd8-a790-536708711a69",
   "metadata": {},
   "source": [
    " e) Using the results of question c) and d) we get a differential equation for $F$:\n",
    "\n",
    "$$\n",
    "\\partial^2_r F(r) = \\frac{4}{3 \\pi} \\frac{1}{Z} \\Bigl( \\frac{2 m}{\\hbar^2 4 \\pi \\epsilon_0} Ze^2  \\Bigr)^{3/2} r^{-1/2} F^{3/2}(r).\n",
    "$$\n",
    "\n",
    "We perform the change of variable $F(r) = \\tilde{F}(x)$ with $r = \\alpha x$ and $\\alpha = \\frac{a_0}{2} (\\frac{3 \\pi}{4})^{2/3} Z^{-1/3}$, where the definition of the Bohr radius $a_0 = \\frac{4 \\pi \\epsilon_0 \\hbar^2}{me}$ is used to simply the expression for $\\alpha$. This leads to the Thomas-Fermi equation for screening:\n",
    "\n",
    "$$\n",
    "\\partial^2_x \\tilde{F}(x) = x^{-1/2}\\tilde{F}^{3/2}(x)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3db8f9f1-57c5-438b-933e-cf9f03f8eeb8",
   "metadata": {},
   "source": [
    "<div style=\"text-align: center; font-size: 22px; font-weight: bold; color: black; margin: 20px;\">\n",
    "    Exercise 2: The universal scaling solution\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb008064-dc5a-446e-b280-f2c662a70913",
   "metadata": {},
   "source": [
    "### Solution exercise 2 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee77b580-e25d-4f77-8e86-bdc735ac3a45",
   "metadata": {},
   "source": [
    "a) Given the change of varible in exercise 1 we obtain the following scaling laws:  \n",
    "- $r \\propto$ $Z^{-1/3}$ (heavier atoms are smaller)\n",
    "- $V \\propto \\frac{Z}{r} F \\propto Z^{4/3}$\n",
    "- $n \\propto \\frac{Z}{r} \\partial^2_r F \\propto Z^{2}$\n",
    "- $E_{kin} \\propto k_F^2 \\propto Z^{4/3}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d512175-0b1f-4f43-b18e-2fff72e83160",
   "metadata": {},
   "source": [
    "b) We propose that for $x \\rightarrow \\infty$, $\\tilde{F} \\propto \\beta x^{-b}$, with $b > 0$. As a leading contribution to the Thomas-Fermi equation we get:\n",
    "$$\n",
    "b(b+1) \\beta x^{1/2 - b -2} + \\mathcal{O}(x^{-5/2 - b}) = \\beta^{3/2} x^{-3b/2} + \\mathcal{O}(x^{-(b+1)3/2}).\n",
    "$$\n",
    "Matching the behaviour on both sides we get:\n",
    "$$\n",
    "-(b+2) + \\frac{1}{2} = -\\frac{3}{2} b\n",
    "$$\n",
    "$$\n",
    "b(b+1)\\beta = \\beta^{3/2}\n",
    "$$\n",
    "which has the solutions $\\beta = 144$ and $b = 3$. From this we follow that $\\tilde{F}(x) \\propto 144x^{-3}$ is the leading behaviour at large $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac03ac1c-b0e1-4309-b2cc-533aba7c99d2",
   "metadata": {},
   "source": [
    "c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6e4ab501-735d-4f15-851b-9716ddc50028",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b75d6c08-09b5-444e-afcc-b53e2e5d696d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_differential_equation(func, x_span, y0, method='LSODA', atol=1e-8, rtol=1e-6, t_eval=None, plot_result=True):\n",
    "    \"\"\"\n",
    "    Solver for nonlinear differential equations.\n",
    "\n",
    "    Parameters:\n",
    "    - func: Callable. Defines the system of equations dy/dx = f(x, y).\n",
    "            It should take the form func(x, y), where y can be a scalar or a list of values.\n",
    "    - x_span: Tuple. The range of x values (x_start, x_end).\n",
    "    - y0: List or array-like. Initial conditions for the system.\n",
    "    - method: String. Integration method ('RK45', 'RK23', 'LSODA', etc.).\n",
    "    - t_eval: Array-like. Points at which to store the computed solution.\n",
    "    - plot_result: Boolean. Whether to plot the result.\n",
    "\n",
    "    Returns:\n",
    "    - solution: ODE solution object from scipy.integrate.solve_ivp.\n",
    "    \"\"\"\n",
    "\n",
    "    # Provide default evaluation points if not specified\n",
    "    if t_eval is None:\n",
    "        t_eval = np.linspace(x_span[0], x_span[1], 500)\n",
    "\n",
    "    # Solve differential equation\n",
    "    sol = solve_ivp(func, x_span, y0, method=method, t_eval=t_eval)\n",
    "\n",
    "    if plot_result:\n",
    "        plt.figure(figsize=(7, 5))\n",
    "        plt.plot(sol.t, sol.y[0])\n",
    "        plt.xlabel('x')\n",
    "        plt.ylabel('y(x)')\n",
    "        plt.show()\n",
    "    return sol\n",
    "    \n",
    "\n",
    "def DGL_system(x, y):\n",
    "    # Rewrite secon order differential equation as a system of first order differential equations:\n",
    "    # z = y'; z' = y''\n",
    "\n",
    "    dydx = y[1]\n",
    "    \n",
    "    f = (1+4/3*x**(3/2))\n",
    "    \n",
    "    if x == 0: \n",
    "        dzdx = - 4*x**(1/2)/f*y[1]\n",
    "    else:\n",
    "        dzdx = y[0]/(x**(1/2)*f)*(f**(3/2)*y[0]**(1/2) - 1) - 4*x**(1/2)/f*y[1]\n",
    "        \n",
    "    return [dydx, dzdx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0301038e-0f04-485f-8e10-997433862643",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAHACAYAAADjth/nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4cElEQVR4nO3de3RU1d3/8c/MJDO5J4SQSQKBIKCCXAUT8fJYH1OpWiy9UrVCsdpHF7VgelGqwK+1ErWKrCqFaqu1XUWxtlKrlBZTb2iUm1Gs3OQaCEkIgUxuZJKZ8/tjkkkiCSSSmTOZeb/WmpXJPvvMfIepzWftvc8+FsMwDAEAAMBUVrMLAAAAAKEMAAAgJBDKAAAAQgChDAAAIAQQygAAAEIAoQwAACAEEMoAAABCAKEMAAAgBESZXUCweb1elZWVKTExURaLxexyAABAmDMMQ7W1tcrKypLV2v14WMSFsrKyMmVnZ5tdBgAAiDClpaUaMmRIt8cjLpQlJiZK8v3DJCUlmVwNAAAIdy6XS9nZ2f4M0p2IC2VtU5ZJSUmEMgAAEDRnWjbFQn8AAIAQQCgDAAAIAYQyAACAEEAoAwAACAGEMgAAgBBAKAMAAAgBhDIAAIAQQCgDAAAIAYQyAACAEEAoAwAACAGmhrK33npL06dPV1ZWliwWi9asWXPGc9544w1deOGFcjgcGjlypP7whz8EvE4AAIBAMzWU1dfXa8KECVq+fHmP+u/bt0/XXXedrrzySpWUlGj+/Pm69dZb9a9//SvAlQIAAASWqTckv+aaa3TNNdf0uP/KlSs1fPhwPfroo5Kk0aNHa8OGDXrsscc0bdq0QJUJAAAQcP1qTVlxcbHy8/M7tU2bNk3FxcXdntPU1CSXy9XpEUiHjjfoK09s0HW/fjug7wMAAMJLvwpl5eXlcjqdndqcTqdcLpcaGxu7PKewsFDJycn+R3Z2dkBrdETZ9OGhGn1yxKUWjzeg7wUAAMJHvwpln8eCBQtUU1Pjf5SWlgb0/VLj7bJaJMOQquvdAX0vAAAQPkxdU9ZbGRkZqqio6NRWUVGhpKQkxcbGdnmOw+GQw+EIRnmSJJvVotR4h6rqmnS0rknpSTFBe28AANB/9auRsqlTp6qoqKhT2/r16zV16lSTKupaWoJdklRVx0gZAADoGVNDWV1dnUpKSlRSUiLJt+VFSUmJDh48KMk39Thr1ix//9tvv1179+7VT3/6U+3YsUO/+c1v9MILL+iuu+4yo/xuDUr0jcxV1TaZXAkAAOgvTA1lmzdv1qRJkzRp0iRJUkFBgSZNmqRFixZJko4cOeIPaJI0fPhwvfrqq1q/fr0mTJigRx99VL/73e9CbjuMQQm+UHa0jlAGAAB6xtQ1ZV/4whdkGEa3x7varf8LX/iCPvjggwBWdfbSGCkDAAC91K/WlPUX7WvKCGUAAKBnCGUBkNY6fclCfwAA0FOEsgBoC2VHmb4EAAA9RCgLAP/Vl0xfAgCAHiKUBUDbSFl1g5tbLQEAgB4hlAVAp1stNbCuDAAAnBmhLAB8t1pqvQKzllAGAADOjFAWIGlsIAsAAHqBUBYg3GoJAAD0BqEsQNr3KiOUAQCAMyOUBQi7+gMAgN4glAUIu/oDAIDeIJQFCLv6AwCA3iCUBQi7+gMAgN4glAUIC/0BAEBvEMoCJC3Rt9C/ut4tj9cwuRoAABDqCGUBkhpnl8UieQ3pWD2jZQAA4PQIZQESZbNqILdaAgAAPUQoCyDWlQEAgJ4ilAUQoQwAAPQUoSyA2NUfAAD0FKEsgNhAFgAA9BShLIDaN5BloT8AADg9QlkAsaYMAAD0FKEsgNISmb4EAAA9QygLoPaF/kxfAgCA0yOUBdCg1unL6vombrUEAABOi1AWQKnx7bdaqq5ntAwAAHSPUBZAUTarUuPYqwwAAJwZoSzAuAITAAD0BKEswNISfSNlXIEJAABOh1AWYIyUAQCAniCUBdigBHb1BwAAZ0YoC7C2DWSrmL4EAACnQSgLMP9NyZm+BAAAp0EoC7C2Xf1Z6A8AAE6HUBZggxJZUwYAAM6MUBZgbaGMWy0BAIDTIZQF2MB4h6ytt1o6xroyAADQDUJZgNmsFg1sXexfyboyAADQDUJZEKQntoWykyZXAgAAQhWhLAj8oczFSBkAAOgaoSwI0hNjJDF9CQAAukcoC4L0JKYvAQDA6RHKgqBt+pINZAEAQHcIZUEwiOlLAABwBoSyIPBPX7LQHwAAdINQFgQdpy8Ng139AQDAqQhlQZDWunms2+NVTWOzydUAAIBQRCgLgphom5JjoyWxrgwAAHSNUBYkbCALAABOh1AWJOxVBgAATodQFiRtu/qzVxkAAOgKoSxI2m9KTigDAACnIpQFySBCGQAAOA1CWZCkJ7Xu6u9iTRkAADgVoSxIuP8lAAA4HUJZkDB9CQAATodQFiRtI2V1TS1qcLeYXA0AAAg1poey5cuXKycnRzExMcrLy9PGjRtP23/ZsmU677zzFBsbq+zsbN111106eTL012klOKIUG22TxBQmAAA4lamhbPXq1SooKNDixYu1detWTZgwQdOmTVNlZWWX/VetWqV77rlHixcv1vbt2/X73/9eq1ev1s9+9rMgV957FoulwwayhDIAANCZqaFs6dKluu222zRnzhyNGTNGK1euVFxcnJ5++uku+7/77ru69NJLdeONNyonJ0dXX321brjhhjOOroUKbrUEAAC6Y1ooc7vd2rJli/Lz89uLsVqVn5+v4uLiLs+55JJLtGXLFn8I27t3r9auXatrr7222/dpamqSy+Xq9DBL267+3GoJAAB8VpRZb1xVVSWPxyOn09mp3el0aseOHV2ec+ONN6qqqkqXXXaZDMNQS0uLbr/99tNOXxYWFurnP/95n9b+eXEFJgAA6I7pC/1744033tCSJUv0m9/8Rlu3btXf/vY3vfrqq7r//vu7PWfBggWqqanxP0pLS4NYcWf+NWVMXwIAgM8wbaQsLS1NNptNFRUVndorKiqUkZHR5TkLFy7UzTffrFtvvVWSNG7cONXX1+v73/++7r33Xlmtp2ZMh8Mhh8PR9x/gc2D6EgAAdMe0kTK73a7JkyerqKjI3+b1elVUVKSpU6d2eU5DQ8Mpwctm820zYRhG4IrtI4PY1R8AAHTDtJEySSooKNDs2bM1ZcoU5ebmatmyZaqvr9ecOXMkSbNmzdLgwYNVWFgoSZo+fbqWLl2qSZMmKS8vT59++qkWLlyo6dOn+8NZKHOyJQYAAOiGqaFs5syZOnr0qBYtWqTy8nJNnDhR69at8y/+P3jwYKeRsfvuu08Wi0X33XefDh8+rEGDBmn69Ol64IEHzPoIveJsnb6srnerqcUjR1ToB0kAABAcFqM/zPv1IZfLpeTkZNXU1CgpKSmo720Yhs5buE7uFq/e/umVyk6NC+r7AwCA4Otp9uhXV1/2dxaLRRlJvtGycheL/QEAQDtCWZD5Q1kNoQwAALQjlAWZM9kXyioYKQMAAB0QyoIso/UKTEbKAABAR4SyIHOypgwAAHSBUBZkGUxfAgCALhDKgoyrLwEAQFcIZUHWNn1Z4WrqF7eGAgAAwUEoC7K2UOZu8ep4Q7PJ1QAAgFBBKAsye5RVA+PtkrgCEwAAtCOUmaB9CpNQBgAAfAhlJmi7ApPF/gAAoA2hzARObrUEAAA+g1BmggymLwEAwGcQykyQkdx6qyVCGQAAaEUoMwHTlwAA4LMIZSbg6ksAAPBZhDITtK0pO97QrKYWj8nVAACAUEAoM0FKXLTsUb5/+kpXk8nVAACAUEAoM4HFYuHG5AAAoBNCmUkyWOwPAAA6IJSZxJnMYn8AANCOUGaSjKTWvcoYKQMAACKUmcbJmjIAANABocwkGUxfAgCADghlJmlb6H+E6UsAACBCmWmyUmIl+UbKPF7D5GoAAIDZCGUmSU90yGqRmj2GqurYQBYAgEhHKDNJlM3qn8IsO9FocjUAAMBshDITtU1hlp1gXRkAAJGOUGaiTH8oY6QMAIBIRygzUVZK6/RlDaEMAIBIRygz0WBGygAAQCtCmYkyk1lTBgAAfAhlJvJPXzJSBgBAxCOUmaht+vJYvVsnmz0mVwMAAMxEKDNRcmy04uw2SdxuCQCASEcoM5HFYlFmMlOYAACAUGa6LK7ABAAAIpSZbjC7+gMAABHKTNe+LQYjZQAARDJCmcnY1R8AAEiEMtOxqz8AAJAIZabL7LCmzDAMk6sBAABmIZSZrG1LjMZmj040NJtcDQAAMAuhzGQx0TalJdglsa4MAIBIRigLAdyYHAAAEMpCADcmBwAAhLIQ4N/Vn+lLAAAiFqEsBGQxfQkAQMQjlIUA7n8JAAAIZSFg8ABfKDt0vMHkSgAAgFkIZSFgSGsoq3A16WSzx+RqAACAGQhlIWBgvF2x0TZJ0mGmMAEAiEiEshBgsViUneobLSutZgoTAIBIRCgLEdkD4iRJh44zUgYAQCQilIWI7FRfKCtlsT8AABGJUBYi2hb7H6pmpAwAgEhEKAsRjJQBABDZTA9ly5cvV05OjmJiYpSXl6eNGzeetv+JEyc0d+5cZWZmyuFw6Nxzz9XatWuDVG3gtK0pY6E/AACRKcrMN1+9erUKCgq0cuVK5eXladmyZZo2bZp27typ9PT0U/q73W598YtfVHp6ul588UUNHjxYBw4cUEpKSvCL72NtV18eb2hWXVOLEhymfjUAACDITP3Lv3TpUt12222aM2eOJGnlypV69dVX9fTTT+uee+45pf/TTz+t6upqvfvuu4qOjpYk5eTkBLPkgEmMiVZKXLRONDSrtLpBozOTzC4JAAAEkWnTl263W1u2bFF+fn57MVar8vPzVVxc3OU5L7/8sqZOnaq5c+fK6XRq7NixWrJkiTye7nfBb2pqksvl6vQIVUxhAgAQuUwLZVVVVfJ4PHI6nZ3anU6nysvLuzxn7969evHFF+XxeLR27VotXLhQjz76qH75y192+z6FhYVKTk72P7Kzs/v0c/Ql/way7FUGAEDEMX2hf294vV6lp6frySef1OTJkzVz5kzde++9WrlyZbfnLFiwQDU1Nf5HaWlpECvuHUbKAACIXKatKUtLS5PNZlNFRUWn9oqKCmVkZHR5TmZmpqKjo2Wz2fxto0ePVnl5udxut+x2+ynnOBwOORyOvi0+QIaktu3qTygDACDSmDZSZrfbNXnyZBUVFfnbvF6vioqKNHXq1C7PufTSS/Xpp5/K6/X623bt2qXMzMwuA1l/kz2g7f6XTF8CABBpTJ2+LCgo0FNPPaVnn31W27dv1x133KH6+nr/1ZizZs3SggUL/P3vuOMOVVdXa968edq1a5deffVVLVmyRHPnzjXrI/SpjhvIGoZhcjUAACCYTN0SY+bMmTp69KgWLVqk8vJyTZw4UevWrfMv/j948KCs1vbcmJ2drX/961+66667NH78eA0ePFjz5s3T3XffbdZH6FODU3wjZQ1uj6rr3RqY0D+mXQEAwNmzGBE2JONyuZScnKyamholJYXeXmB5S15ThatJa+ZeqonZKWaXAwAAzlJPs0e/uvoyEnAFJgAAkYlQFmK4MTkAAJGJUBZiuAITAIDIRCgLMexVBgBAZCKUhZhhraHswDFCGQAAkYRQFmJy0uIl+UbK3C3eM/QGAADhglAWYtITHYqz2+Q1WOwPAEAkIZSFGIvFomEDfaNl+6vqTa4GAAAEC6EsBA1P860r28+6MgAAIgahLATlMFIGAEDEIZSFoLbF/vuPEcoAAIgUhLIQNLw1lO1jpAwAgIhBKAtBbdOXZSca1dTiMbkaAAAQDISyEJSWYFeCI8q3LQY3JgcAICIQykKQb1sM3xWY+6oIZQAARAJCWYjyL/ZnXRkAABGBUBaihreuK9vHFZgAAESEqM9zUnNzs8rLy9XQ0KBBgwYpNTW1r+uKeIyUAQAQWXo8UlZbW6sVK1boiiuuUFJSknJycjR69GgNGjRIw4YN02233aZNmzYFstaI0rar/wF29QcAICL0KJQtXbpUOTk5euaZZ5Sfn681a9aopKREu3btUnFxsRYvXqyWlhZdffXV+tKXvqTdu3cHuu6w598Wo6ZRJ5vZFgMAgHDXo+nLTZs26a233tIFF1zQ5fHc3FzdcsstWrlypZ555hm9/fbbGjVqVJ8WGmlS4+1KjIlS7ckWHaxu0LnORLNLAgAAAdSjUPbcc8/16MUcDoduv/32syoIPhaLRcPT4vXRoRrtq6onlAEAEOZ6ffXl0aNHuz22bdu2syoGnXFjcgAAIkevQ9m4ceP06quvntL+yCOPKDc3t0+Kgk9O6way3JgcAIDw1+tQVlBQoK9//eu644471NjYqMOHD+uqq67Sww8/rFWrVgWixojVti3G3qOEMgAAwl2vQ9lPf/pTFRcX6+2339b48eM1fvx4ORwOffTRR/rqV78aiBoj1ohBCZKkPYQyAADC3ufa0X/kyJEaO3as9u/fL5fLpZkzZyojI6Ova4t4I9J9oayqrkknGtwmVwMAAAKp16HsnXfe0fjx47V792599NFHWrFihe68807NnDlTx48fD0SNESvBEaWs5BhJ0qeVdSZXAwAAAqnXoex///d/NXPmTL333nsaPXq0br31Vn3wwQc6ePCgxo0bF4gaI9rI1q0wCGUAAIS3Xt/78t///reuuOKKTm0jRozQO++8owceeKDPCoPPyEEJemvXUe0mlAEAENZ6PVL22UDmfyGrVQsXLjzrgtDZyNZ1ZYyUAQAQ3noUyp5//vkev2Bpaaneeeedz10QOhvlJJQBABAJehTKVqxYodGjR+vhhx/W9u3bTzleU1OjtWvX6sYbb9SFF16oY8eO9XmhkWpk67YYh080qr6pxeRqAABAoPRoTdmbb76pl19+WY8//rgWLFig+Ph4OZ1OxcTE6Pjx4yovL1daWpq++93v6uOPP5bT6Qx03RFjQLxdA+PtOlbv1t6j9Ro3JNnskgAAQAD0eKH/9ddfr+uvv15VVVXasGGDDhw4oMbGRqWlpWnSpEmaNGmSrNbPte0ZzmBkeoKO7avWp0drCWUAAISpXl99+aMf/Ujf+973NGPGjACUg66MTE/Q+/uqtbuCdWUAAISrXg9t1dTUKD8/X6NGjdKSJUtUVlYWiLrQAVdgAgAQ/nodytasWaPDhw/rjjvu0OrVqzVs2DBdc801evHFF9Xc3ByIGiPeqHQ2kAUAINx9rkVggwYNUkFBgT788EO9//77GjlypG6++WZlZWXprrvu0u7du/u6zojWNlJ2oLpB7havydUAAIBAOKuV+UeOHNH69eu1fv162Ww2XXvttdq2bZvGjBmjxx57rK9qjHjOJIcSHFHyeA3tP1ZvdjkAACAAeh3Kmpub9de//lVf/vKXNWzYMP3lL3/R/PnzVVZWpmeffVavvfaaXnjhBf3iF78IRL0RyWKx+EfLWOwPAEB46vXVl5mZmfJ6vbrhhhu0ceNGTZw48ZQ+V155pVJSUvqgPLQZmZ6gktITrCsDACBM9TqUPfbYY/rmN7+pmJiYbvukpKRo3759Z1UYOhvVNlJWWWtyJQAAIBB6HcpuvvnmQNSBMzjX6bsCc1cFoQwAgHDEFvz9xHkZvlC292i9mlo8JlcDAAD6GqGsn8hMjlFSTJRavIb2VHIFJgAA4YZQ1k9YLBadn5kkSdpR7jK5GgAA0NcIZf3I6NYpzB3lrCsDACDcEMr6kfMy2kbKCGUAAIQbQlk/cn5m60jZEaYvAQAIN4SyfuS81m0xKmubVF3vNrkaAADQlwhl/Ui8I0rDBsZJYrE/AADhhlDWz7SNlu04wroyAADCCaGsn2nbFmM768oAAAgrhLJ+ZkzrYv/tTF8CABBWCGX9zAVZyZKkneW1crd4Ta4GAAD0FUJZPzNkQKySYqLU7DG0u5J1ZQAAhAtCWT9jsVj8o2X/PcwUJgAA4YJQ1g+NHexb7P9xWY3JlQAAgL4SEqFs+fLlysnJUUxMjPLy8rRx48Yenff888/LYrFoxowZgS0wxIwd3DpSVsZIGQAA4cL0ULZ69WoVFBRo8eLF2rp1qyZMmKBp06apsrLytOft379fP/7xj3X55ZcHqdLQcUGWb6TskzKXPF7D5GoAAEBfMD2ULV26VLfddpvmzJmjMWPGaOXKlYqLi9PTTz/d7Tkej0c33XSTfv7zn+ucc84JYrWhYXhagmKjbWps9mhfVZ3Z5QAAgD5gaihzu93asmWL8vPz/W1Wq1X5+fkqLi7u9rxf/OIXSk9P1/e+970zvkdTU5NcLlenR39ns1o0unW/MqYwAQAID6aGsqqqKnk8Hjmdzk7tTqdT5eXlXZ6zYcMG/f73v9dTTz3Vo/coLCxUcnKy/5GdnX3WdYeCtnVlHx9msT8AAOHA9OnL3qitrdXNN9+sp556SmlpaT06Z8GCBaqpqfE/SktLA1xlcIxt3Rbjo0OEMgAAwkGUmW+elpYmm82mioqKTu0VFRXKyMg4pf+ePXu0f/9+TZ8+3d/m9fp2tY+KitLOnTs1YsSITuc4HA45HI4AVG+u8dntI2UeryGb1WJyRQAA4GyYOlJmt9s1efJkFRUV+du8Xq+Kioo0derUU/qff/752rZtm0pKSvyP66+/XldeeaVKSkrCZmqyJ0alJyrOblO926M9R1nsDwBAf2fqSJkkFRQUaPbs2ZoyZYpyc3O1bNky1dfXa86cOZKkWbNmafDgwSosLFRMTIzGjh3b6fyUlBRJOqU93NmsFo0dnKyN+6r1YekJnetMNLskAABwFkwPZTNnztTRo0e1aNEilZeXa+LEiVq3bp1/8f/BgwdltfarpW9BMzE7xRfKDp3QN6dEzighAADhyGIYRkTtPupyuZScnKyamholJSWZXc5ZeeWjMv1g1QcaNzhZ/7jzMrPLAQAAXehp9mAIqh+bMCRFkrSj3KWTzR5ziwEAAGeFUNaPDRkQq4HxdjV7DG0/wiayAAD0Z4SyfsxisWj8EN/WGB+WnjC3GAAAcFYIZf3chOwUSdKHbCILAEC/Rijr5y4cOkCStOXAcZMrAQAAZ4NQ1s9NHJoii0U6WN2gytqTZpcDAAA+J0JZP5cUE63zWjeO3XrghLnFAACAz41QFgYuHOabwtx6kClMAAD6K0JZGJjcuq5s8/5qkysBAACfF6EsDEzJ8YWyjw+ziSwAAP0VoSwMDE2NU1qCXW6PV/8tY2sMAAD6I0JZGLBYLGyNAQBAP0coCxOTWxf7b9xHKAMAoD8ilIWJ3OGpkqRN+6vl9RomVwMAAHqLUBYmxg5OVrzdpprGZu0orzW7HAAA0EuEsjARbbNqco5vtOz9fcdMrgYAAPQWoSyM5LVOYb6/l/3KAADobwhlYeTic3yhbOP+ahkG68oAAOhPCGVhZNzgFMVEW1Vd79buyjqzywEAAL1AKAsj9iirf2uM9/eyrgwAgP6EUBZmLh4+UJL07h5CGQAA/QmhLMxcMjJNki+UedivDACAfoNQFmYmDElWoiNKNY3N3AcTAIB+hFAWZqJsVl08wjeFueHTKpOrAQAAPUUoC0OXtU5hbthNKAMAoL8glIWhy0b5Qtnm/cfV6PaYXA0AAOgJQlkYOictXpnJMXJ7vNp8gN39AQDoDwhlYchisejS1inMt5nCBACgXyCUhakrzh0kSXp9R6XJlQAAgJ4glIWp/zl3kGxWi3ZX1qm0usHscgAAwBkQysJUcmy0prTecuk/jJYBABDyCGVh7H/PT5ckFRHKAAAIeYSyMHbVaF8oe2/PMdU3tZhcDQAAOB1CWRgbMShBQ1Pj5PZ49Q67+wMAENIIZWHMYrH4pzBZVwYAQGgjlIW5jqHMMAyTqwEAAN0hlIW5vHNSFWe3qbK2Sf8tc5ldDgAA6AahLMw5omz+G5S/tr3C5GoAAEB3CGUR4ItjnJKkdR+Xm1wJAADoDqEsAlw9JkPRNot2lNfq08o6s8sBAABdIJRFgOS4aP8NytduO2JyNQAAoCuEsghx7bhMSYQyAABCFaEsQlw9xqkoq28Kc89RpjABAAg1hLIIkRJnb5/C/IjRMgAAQg2hLIJc1zqF+SpTmAAAhBxCWQS5+gKmMAEACFWEsgjScQrzlQ8ZLQMAIJQQyiLMVyZmSZL+9sEh7oUJAEAIIZRFmC+NzVC83aYDxxq0+cBxs8sBAACtCGURJs4epevG+xb8v7j5kMnVAACANoSyCPSNydmSfFdhNrhbTK4GAABIhLKIdFHOAA1NjVNdU4v+9V9uUg4AQCgglEUgi8Wib0weIkl6cQtTmAAAhAJCWYT62oWDJUnv7jmmQ8cbTK4GAAAQyiLUkAFxumTEQBkGo2UAAIQCQlkE+3buUEnSqvcPqtnjNbkaAAAiG6Esgn3pggylJThUWdukf/+3wuxyAACIaISyCGaPsurGPN9o2bPF+80tBgCACEcoi3A35g6VzWrRxn3V2lHuMrscAAAiVkiEsuXLlysnJ0cxMTHKy8vTxo0bu+371FNP6fLLL9eAAQM0YMAA5efnn7Y/Ti8jOUZfuiBDkvTH4gMmVwMAQOQyPZStXr1aBQUFWrx4sbZu3aoJEyZo2rRpqqys7LL/G2+8oRtuuEGvv/66iouLlZ2drauvvlqHDx8OcuXh4+apwyRJL209rJrGZpOrAQAgMlkMwzDMLCAvL08XXXSRnnjiCUmS1+tVdna27rzzTt1zzz1nPN/j8WjAgAF64oknNGvWrDP2d7lcSk5OVk1NjZKSks66/nBgGIa+tOxt7ayo1b3XjtZt/3OO2SUBABA2epo9TB0pc7vd2rJli/Lz8/1tVqtV+fn5Ki4u7tFrNDQ0qLm5WampqV0eb2pqksvl6vRAZxaLRbdcliNJ+t2GvWpq8ZhbEAAAEcjUUFZVVSWPxyOn09mp3el0qry8Z/dkvPvuu5WVldUp2HVUWFio5ORk/yM7O/us6w5HX500RBlJMapwNemlrUwFAwAQbKavKTsbDz74oJ5//nm99NJLiomJ6bLPggULVFNT43+UlpYGucr+wR5l1a2XD5ck/fatvfJ4TZ3VBgAg4pgaytLS0mSz2VRR0Xnj0oqKCmVkZJz23EceeUQPPvig/v3vf2v8+PHd9nM4HEpKSur0QNduyB2qlLho7auq17qPezZSCQAA+oapocxut2vy5MkqKiryt3m9XhUVFWnq1Kndnvfwww/r/vvv17p16zRlypRglBoR4h1Rmj01R5L0mzc+lcnXgAAAEFFMn74sKCjQU089pWeffVbbt2/XHXfcofr6es2ZM0eSNGvWLC1YsMDf/6GHHtLChQv19NNPKycnR+Xl5SovL1ddXZ1ZHyGsfPeSHMVG2/TfMpfe2HXU7HIAAIgYpoeymTNn6pFHHtGiRYs0ceJElZSUaN26df7F/wcPHtSRI0f8/VesWCG3261vfOMbyszM9D8eeeQRsz5CWBkQb9d3Lvbdemnpv3cxWgYAQJCYvk9ZsLFP2Zkdq2vS/zz8uurdHv3mpgt17bhMs0sCAKDf6hf7lCE0DUxw6NbLfRvIPvLvnWrxeE2uCACA8EcoQ5duvXy4BsRFa+/Rev116yGzywEAIOwRytClxJhozb1ypCRp2Wu7dbKZXf4BAAgkQhm69Z2LhykzOUZHak7qmXf2m10OAABhjVCGbsVE2/Tjq8+TJD3+n92qcJ00uSIAAMIXoQyn9dVJgzVpaIoa3B4Vrt1udjkAAIQtQhlOy2q16BfXj5XFIq0pKdPm/dVmlwQAQFgilOGMxg1J1rcvypYkLfr7f7lZOQAAAUAoQ4/8+OrzlBQTpU+OuPTsu/vNLgcAgLBDKEOPDExw6KdfOl+S9Kt/7dTBYw0mVwQAQHghlKHHbswdqovPSVVjs0d3//Uj7osJAEAfIpShx6xWix76+njFRFtVvPeYnttYanZJAACEDUIZemXYwHj9ZJpvGnPJ2u06dJxpTAAA+gKhDL323UtyNHnYANU1tWj+8yXcsBwAgD5AKEOv2awWPfatiUp0RGnzgeP6ddFus0sCAKDfI5Thcxk6ME4PfG2cJOnx1z9V8Z5jJlcEAED/RijD53b9hCx9a8oQGYY0f/UHOlbXZHZJAAD0W4QynJX/d/0FGjEoXhWuJt3x561qZn0ZAACfC6EMZyXOHqXf3jxZCY4obdxXrftf+cTskgAA6JcIZThrI9MTtWzmRFks0h+LD+i5jQfNLgkAgH6HUIY+kT/GqR998VxJ0qK/f6x391SZXBEAAP0LoQx9Zu6VI3Xd+Ew1ewz93x+36JMyl9klAQDQbxDK0GcsFose/eYE5Q1PVW1Ti777zEaVVrPjPwAAPUEoQ5+KibbpyVlTdH5GoiprmzT76Y06WstWGQAAnAmhDH0uOTZaf5iTq8EpsdpbVa+bfveeqtjDDACA0yKUISAykmP051vz5ExyaFdFnW566n02lwUA4DQIZQiYnLR4Pf/9qUpPdGhnRa1u+t37qqw9aXZZAACEJEIZAmp4Wrye+/7FSk90aEd5rb6xolgHjtWbXRYAACGHUIaAGzEoQX+5faqGpsbpYHWDvr6iWP8tqzG7LAAAQgqhDEExbGC8Xrx9qkZnJqmqrknf/u17emNnpdllAQAQMghlCJr0pBit/r+L/fuY3fKHTfr9hn0yDMPs0gAAMB2hDEGVFBOtP30vTzOnZMtrSPe/8onu/utHOtnsMbs0AABMRShD0NmjrHrw6+N033WjZbVIL2w+pBnL39GnlXVmlwYAgGkIZTCFxWLRrZefoz/MydXAeLt2lNfq+ic26G9bD5ldGgAApiCUwVT/c+4g/XPe5bpkxEA1uD0qeOFDFbxQoprGZrNLAwAgqAhlMF16Uoz+9L08FXzxXFkt0t+2HtYXl76p9Z9UmF0aAABBQyhDSLBZLfrhVaP0wv9N1Tlp8aqsbdJtf9ysHz73AbdnAgBEBEIZQsqUnFStnXe5br9ihKwW6eUPy3TV0jf1x+L9avF4zS4PAICAsRgRtkmUy+VScnKyampqlJSUZHY5OI2PDp3QT1/8SDvKayVJ5zkTtXj6GF0yMs3kygAA6LmeZg9CGUJai8er5zYe1KPrd+lEg2/xf/5op3509bkancn3BwAIfYSybhDK+qcTDW4te223/vTeAXm8vv/Jfnl8pubnn6uR6QkmVwcAQPcIZd0glPVvn1bWadlru/TKR0ckSVaLNH1Clv7vf0ZoTBbfJwAg9BDKukEoCw+flLn02Gu7Om2bcfmoNN1+xQhdMmKgLBaLidUBANCOUNYNQll42XaoRr99a4/Wbjui1llNnZ+RqJvyhmrGpMFKjIk2t0AAQMQjlHWDUBaeSqsb9Lu39+qFzYfU2Hpz8zi7TV+ZmKUbc4dp3JBkkysEAEQqQlk3CGXhraahWX/74JD+/P7BTjc4H52ZpK9MzNL1E7KUlRJrYoUAgEhDKOsGoSwyGIahjfuqtWrjQf1zW7ncHTaezR2eqhkTB+tLYzOUGm83sUoAQCQglHWDUBZ5TjS4tXZbudaUHNbGfdX+dqtFmjIsVV8c49QXxziVkxZvYpUAgHBFKOsGoSyylZ1o1MsflukfH5bpv2WuTsfOdSboC+el67KRacodnqqYaJtJVQIAwgmhrBuEMrQ5dLxBr31SofXbK/T+3mq1eNv/U7BHWTVl2ABdNipNl48cpDFZSbJZ2WYDANB7hLJuEMrQlZqGZr2xq1Ibdldpw6dVOlJzstPxBEeUJg1N0ZRhqZqSM0ATs1MU74gyqVoAQH9CKOsGoQxnYhiG9lbV+wPae3uOqbappVMfm9Wi0ZmJGj8kRWOzkjVucLLOzUiQI4opTwBAZ4SybhDK0Fser6Gd5bXacqBamw8c1+b9x3X4ROMp/aJtFp3rTNTYrGSNHZykc52JOteZqAFc4QkAEY1Q1g1CGfrCkZpGbT1wQh+X1ejjwzXadrhGJxqau+yblmDXqPREjXImaJQzUaPSEzQqPUGp8XZuBwUAEYBQ1g1CGQLBMAwdPtHoD2iflLm0q6KuyxG1NomOKA0dGKdhA+M0NDVewwbGaVhqnIYOjFNmciwXFgBAmCCUdYNQhmCqb2rRnqN12lVRp92Vtdrd+rO0uvuwJkl2m1WDB8QqMzlGmcmxykqJUUZyjLKSY5WZEqPMpFglxUYx0gYA/UBPsweXjwEBFO+I0vghKRo/JKVT+8lmj0qrG3TgWIMOVDfo4LH61p8NKj3eILfHq31V9dpXVd/ta8fZbcpM9oW1tARHh4ddaYkODWr9fWCCXdE2a4A/KQDgbBHKABPERNt868uciacc83gNHalpVGl1o8pdjSo7cVJHahpVXnPS//x4Q7Ma3B7tOVqvPUe7D25tUuKilZbgUGq8XQPiopUSa1dKXLRS4lp/xnZ4HhetAXF2Ns8FgCAjlAEhxma1aMiAOA0ZENdtn5PNHh2pOakjJxpVUXtSVbVuVdU16Whdk6rq3KqqbVJVXZOO1bvl8Ro60dDc7YUI3XFEWZUSF62kmGglxEQpwRGlxNafCQ5fW6Ijyn8sISZKSTHtxxLsUYpz2BilA4AeColQtnz5cv3qV79SeXm5JkyYoMcff1y5ubnd9v/LX/6ihQsXav/+/Ro1apQeeughXXvttUGsGDBXTLRNw9PiNfwM9+v0eg0db3D7glpdk6rr3TrR2KyaBrdONDTreEOzahrbnrtV0+gLby1eQ00tXlW4mlThajqrWqOsFsXabYqz2xQbbVOsParDc1s3z9v7OKKsckRb5YhqfR5lkz3Kekq7PcqqKKuFdXYA+i3TQ9nq1atVUFCglStXKi8vT8uWLdO0adO0c+dOpaenn9L/3Xff1Q033KDCwkJ9+ctf1qpVqzRjxgxt3bpVY8eONeETAKHLarVoYIJDAxMcOk+nTpV2xTAM1bs9Ol7vC2u1Tc2qO9miuibfo7bt+ckW1Z5s7tzW1t7UIneLV5LU4jVUe9LXJ9CsFvlCWrTVH9Taw1x7oIu2WRRtsyrK1vrcalV0lEVR1vZwF23r3M9us7T279Du72dVlM1yyjGb1feaVqs6/bT5j/l+2iwWWbnaFoh4pl99mZeXp4suukhPPPGEJMnr9So7O1t33nmn7rnnnlP6z5w5U/X19XrllVf8bRdffLEmTpyolStXnvH9uPoSCI6mFo9Our1qaG5Rg9ujRrdHjc2eDs87tLs9amhuf+7v1+wLd00tXjU1e9XU4lFTi7e9rcWjZk94XEBusUg2S3tYs/pDm1W2rkJda9+OAc/aMei19rFYLLJaJKvFIqu19WeHNovFIpu1/Xlbu81qkaXtPP/5nzmv7ViHvr737HCe1dLFe/qm6a2d+vqeW1r/LSwWS+tzi7/N2tqhrd1qkSzqeN6pz31Zt729/X18P9WxrbXd2vr8dOd1rrH9PHXoY/3McbXW89nz2r7/jv9baH13/+8W/7HWto79GB0Oef3i6ku3260tW7ZowYIF/jar1ar8/HwVFxd3eU5xcbEKCgo6tU2bNk1r1qzpsn9TU5OamtqnX1wu19kXDuCMfCNUNiUrOqDv4/EarSHN0ymsnWxuf94pyDV71OI11OLxyu3x/Wz2eNXsMdTs8aql9fVavF61eAy5Pb6fnft41dxiqNnb8ZjveNvrery+1/J6Dd9PwzhtgDQMqcXw9T27CWNEuu5CXPsxf+rr9pjllGP+s7o4z9Jl/7bfLZ85sWPfrgJo+7HOYbMtDJ+pf1d1dvXv8NnXefLmKRo6sPu1vMFgaiirqqqSx+OR0+ns1O50OrVjx44uzykvL++yf3l5eZf9CwsL9fOf/7xvCgYQcmyta9Zi7f3jalGv15DHMOTx+h4t3vbnvt+98nrl+9ka0lo8hv+5t8tzOvxutAZCjyHDkLyGIa//p+/8tt8NQ/K0thuG/Mc8hiGjrb+//TPndWhr6+vxqvN5/vdUp76eU87z/S5DMuR7fUPt72VIUts5bcdb2w1/n/bP6+vffvyU81r7dD7Pd9Ip793h9byfOe90NZvJ/2/W9supPYJYTf/h9njMLsH8NWWBtmDBgk4jay6XS9nZ2SZWBCCSWa0WWWURO46EP6ND+PxsCGx7LnUOUW0ritr6+3/Rmfu3HzM+c96px4xTjp0aJo0O79tV/1Pq7OK1jS5q766Wju/bVZ1dfVb/q3fxWf2fq4vP07Fz2/GslFiZzdRQlpaWJpvNpoqKik7tFRUVysjI6PKcjIyMXvV3OBxyOBx9UzAAAD3kW3cndZhQA07L1A2E7Ha7Jk+erKKiIn+b1+tVUVGRpk6d2uU5U6dO7dRfktavX99tfwAAgP7A9OnLgoICzZ49W1OmTFFubq6WLVum+vp6zZkzR5I0a9YsDR48WIWFhZKkefPm6YorrtCjjz6q6667Ts8//7w2b96sJ5980syPAQAAcFZMD2UzZ87U0aNHtWjRIpWXl2vixIlat26dfzH/wYMHZbW2D+hdcsklWrVqle677z797Gc/06hRo7RmzRr2KAMAAP2a6fuUBRv7lAEAgGDqafbgpnQAAAAhgFAGAAAQAghlAAAAIYBQBgAAEAIIZQAAACGAUAYAABACCGUAAAAhgFAGAAAQAghlAAAAIYBQBgAAEAJMv/dlsLXdVcrlcplcCQAAiARtmeNMd7aMuFBWW1srScrOzja5EgAAEElqa2uVnJzc7fGIuyG51+tVWVmZEhMTZbFYAvIeLpdL2dnZKi0t5abnJuO7CA18D6GD7yI08D2EjmB8F4ZhqLa2VllZWbJau185FnEjZVarVUOGDAnKeyUlJfEfW4jguwgNfA+hg+8iNPA9hI5AfxenGyFrw0J/AACAEEAoAwAACAGEsgBwOBxavHixHA6H2aVEPL6L0MD3EDr4LkID30PoCKXvIuIW+gMAAIQiRsoAAABCAKEMAAAgBBDKAAAAQgChDAAAIAQQyvrY8uXLlZOTo5iYGOXl5Wnjxo1mlxRxCgsLddFFFykxMVHp6emaMWOGdu7caXZZkPTggw/KYrFo/vz5ZpcScQ4fPqzvfOc7GjhwoGJjYzVu3Dht3rzZ7LIijsfj0cKFCzV8+HDFxsZqxIgRuv/++894T0ScvbfeekvTp09XVlaWLBaL1qxZ0+m4YRhatGiRMjMzFRsbq/z8fO3evTuoNRLK+tDq1atVUFCgxYsXa+vWrZowYYKmTZumyspKs0uLKG+++abmzp2r9957T+vXr1dzc7Ouvvpq1dfXm11aRNu0aZN++9vfavz48WaXEnGOHz+uSy+9VNHR0frnP/+pTz75RI8++qgGDBhgdmkR56GHHtKKFSv0xBNPaPv27XrooYf08MMP6/HHHze7tLBXX1+vCRMmaPny5V0ef/jhh/XrX/9aK1eu1Pvvv6/4+HhNmzZNJ0+eDF6RBvpMbm6uMXfuXP/vHo/HyMrKMgoLC02sCpWVlYYk48033zS7lIhVW1trjBo1yli/fr1xxRVXGPPmzTO7pIhy9913G5dddpnZZcAwjOuuu8645ZZbOrV97WtfM2666SaTKopMkoyXXnrJ/7vX6zUyMjKMX/3qV/62EydOGA6Hw3juueeCVhcjZX3E7XZry5Ytys/P97dZrVbl5+eruLjYxMpQU1MjSUpNTTW5ksg1d+5cXXfddZ3++0DwvPzyy5oyZYq++c1vKj09XZMmTdJTTz1ldlkR6ZJLLlFRUZF27dolSfrwww+1YcMGXXPNNSZXFtn27dun8vLyTv8flZycrLy8vKD+DY+4G5IHSlVVlTwej5xOZ6d2p9OpHTt2mFQVvF6v5s+fr0svvVRjx441u5yI9Pzzz2vr1q3atGmT2aVErL1792rFihUqKCjQz372M23atEk//OEPZbfbNXv2bLPLiyj33HOPXC6Xzj//fNlsNnk8Hj3wwAO66aabzC4topWXl0tSl3/D244FA6EMYW3u3Ln6+OOPtWHDBrNLiUilpaWaN2+e1q9fr5iYGLPLiVher1dTpkzRkiVLJEmTJk3Sxx9/rJUrVxLKguyFF17Qn//8Z61atUoXXHCBSkpKNH/+fGVlZfFdgIX+fSUtLU02m00VFRWd2isqKpSRkWFSVZHtBz/4gV555RW9/vrrGjJkiNnlRKQtW7aosrJSF154oaKiohQVFaU333xTv/71rxUVFSWPx2N2iREhMzNTY8aM6dQ2evRoHTx40KSKItdPfvIT3XPPPfr2t7+tcePG6eabb9Zdd92lwsJCs0uLaG1/p83+G04o6yN2u12TJ09WUVGRv83r9aqoqEhTp041sbLIYxiGfvCDH+ill17Sf/7zHw0fPtzskiLWVVddpW3btqmkpMT/mDJlim666SaVlJTIZrOZXWJEuPTSS0/ZFmbXrl0aNmyYSRVFroaGBlmtnf/02mw2eb1ekyqCJA0fPlwZGRmd/oa7XC69//77Qf0bzvRlHyooKNDs2bM1ZcoU5ebmatmyZaqvr9ecOXPMLi2izJ07V6tWrdLf//53JSYm+tcDJCcnKzY21uTqIktiYuIpa/ni4+M1cOBA1vgF0V133aVLLrlES5Ys0be+9S1t3LhRTz75pJ588kmzS4s406dP1wMPPKChQ4fqggsu0AcffKClS5fqlltuMbu0sFdXV6dPP/3U//u+fftUUlKi1NRUDR06VPPnz9cvf/lLjRo1SsOHD9fChQuVlZWlGTNmBK/IoF3nGSEef/xxY+jQoYbdbjdyc3ON9957z+ySIo6kLh/PPPOM2aXBMNgSwyT/+Mc/jLFjxxoOh8M4//zzjSeffNLskiKSy+Uy5s2bZwwdOtSIiYkxzjnnHOPee+81mpqazC4t7L3++utd/m2YPXu2YRi+bTEWLlxoOJ1Ow+FwGFdddZWxc+fOoNZoMQy2EQYAADAba8oAAABCAKEMAAAgBBDKAAAAQgChDAAAIAQQygAAAEIAoQwAACAEEMoAAABCAKEMAAAgBBDKAAAAQgChDAAAIAQQygDgNI4ePaqMjAwtWbLE3/buu+/KbrerqKjIxMoAhBvufQkAZ7B27VrNmDFD7777rs477zxNnDhRX/nKV7R06VKzSwMQRghlANADc+fO1WuvvaYpU6Zo27Zt2rRpkxwOh9llAQgjhDIA6IHGxkaNHTtWpaWl2rJli8aNG2d2SQDCDGvKAKAH9uzZo7KyMnm9Xu3fv9/scgCEIUbKAOAM3G63cnNzNXHiRJ133nlatmyZtm3bpvT0dLNLAxBGCGUAcAY/+clP9OKLL+rDDz9UQkKCrrjiCiUnJ+uVV14xuzQAYYTpSwA4jTfeeEPLli3Tn/70JyUlJclqtepPf/qT3n77ba1YscLs8gCEEUbKAAAAQgAjZQAAACGAUAYAABACCGUAAAAhgFAGAAAQAghlAAAAIYBQBgAAEAIIZQAAACGAUAYAABACCGUAAAAhgFAGAAAQAghlAAAAIYBQBgAAEAL+P8+cIKr9txdDAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 700x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define x range and initial conditions\n",
    "x_range = (0, 10)\n",
    "initial_conditions = [1, -1.5886]  # [y(0), y'(0)]\n",
    "\n",
    "# Solve differential equation with initial conditions\n",
    "solution = solve_differential_equation(DGL_system, x_range, initial_conditions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67dad927-cead-4f9e-b25e-f2351389ea64",
   "metadata": {},
   "source": [
    "d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "aa13885e-a74f-446e-b104-f4e5b3738240",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8IUlEQVR4nO3de3yT9cH//3eStkkpPVJ6gkLLSUDOIFjwOKtMvdnY3OR2Ohhubvplm9jf7RQV2G6nqFPGnAiTqXNu3qJOnVPUsc7DkCpYQFHOx5ZC0xbapue0yfX7o22wApVCkitNX8/HIw/Iletq3s0meT8+1+f6XBbDMAwBAACECavZAQAAAPyJcgMAAMIK5QYAAIQVyg0AAAgrlBsAABBWKDcAACCsUG4AAEBYiTA7QLB5vV4dPnxYsbGxslgsZscBAACnwTAM1dTUKCMjQ1Zr52MzPa7cHD58WJmZmWbHAAAAZ6C4uFj9+/fvdJ8eV25iY2MltX44cXFxJqcBAACnw+VyKTMz0/c93pkeV27aT0XFxcVRbgAA6GZOZ0oJE4oBAEBYodwAAICwQrkBAABhhXIDAADCCuUGAACEFcoNAAAIK5QbAAAQVig3AAAgrFBuAABAWKHcAACAsGJquXn//fc1Y8YMZWRkyGKx6NVXX/3KY959911NmDBBdrtdQ4YM0Z/+9KeA5wQAAN2HqeWmrq5OY8eO1fLly09r//379+vqq6/WpZdeqi1btmj+/Pn60Y9+pLfffjvASQEAQHdh6o0zr7zySl155ZWnvf/KlSuVnZ2tRx55RJI0YsQIrVu3Tr/97W81ffr0QMU8LR6voaO1Tap3e5SVHGNqFgAAerJuNeemoKBAubm5HbZNnz5dBQUFpzymqalJLperwyMQ1u+t0OT78/XjZz8OyM8HAACnp1uVm9LSUqWmpnbYlpqaKpfLpYaGhpMes2TJEsXHx/semZmZAcmWEuuQJJXVNAXk5wMAgNPTrcrNmViwYIGqq6t9j+Li4oC8T2qcXZJUVd+sxmZPQN4DAAB8NVPn3HRVWlqanE5nh21Op1NxcXGKjo4+6TF2u112uz3g2eKjIxUVYZW7xavymiZlJvUK+HsCAIATdauRm5ycHOXn53fYtnbtWuXk5JiU6DiLxaK+vVtLFKemAAAwj6nlpra2Vlu2bNGWLVsktV7qvWXLFhUVFUlqPaU0e/Zs3/4333yz9u3bp1/84hfasWOHHn/8cb3wwgu67bbbzIh/gvZTU2WuRpOTAADQc5labj7++GONHz9e48ePlyTl5eVp/PjxWrRokSTpyJEjvqIjSdnZ2XrjjTe0du1ajR07Vo888oj++Mc/mn4ZeDsmFQMAYD5T59xccsklMgzjlK+fbPXhSy65RJs3bw5gqjOX0jZy42TkBgAA03SrOTehLjWOkRsAAMxGufGjvrFMKAYAwGyUGz9KiWVCMQAAZqPc+BGnpQAAMB/lxo/aR26O1bnlbvGanAYAgJ6JcuNHib2iFGG1SJLKaxm9AQDADJQbP7JaLcy7AQDAZJQbP+vLvBsAAExFufGzVEZuAAAwFeXGz9pXKWbkBgAAc1Bu/Mx3fykX5QYAADNQbvwsrW3OzeHqBpOTAADQM1Fu/CwjIVqSdKSaOTcAAJiBcuNnGQltIzdVDZ3e8RwAAAQG5cbP2kdu6t0eVTc0m5wGAICeh3LjZ45Im/rEREmSSqqYdwMAQLBRbgKgffTmcBXzbgAACDbKTQB8cd4NAAAILspNABwfuaHcAAAQbJSbAOjXVm6YcwMAQPBRbgKAkRsAAMxDuQkAJhQDAGAeyk0AtE8odtY0qtnjNTkNAAA9C+UmAJJj7IqyWWUYUim3YQAAIKgoNwFgtVqUzuXgAACYgnITIBnxbfNuuDs4AABBRbkJkPZJxSWVlBsAAIKJchMgmUmt5ab4GOUGAIBgotwESGZiL0lScWW9yUkAAOhZKDcBMqBPa7kpOka5AQAgmCg3ATIgqbXcHK5qYK0bAACCiHITIH1722WPsMprSEdYqRgAgKCh3ASI1WpRZhKnpgAACDbKTQBlJrZeMUW5AQAgeCg3AdQ+74YrpgAACB7KTQBxWgoAgOCj3ASQb+SGcgMAQNBQbgIok3IDAEDQUW4CqL3cVNY3y9XYbHIaAAB6BspNAPW2R6hPTJQkRm8AAAgWyk2A+SYVH6XcAAAQDJSbAMtOjpEk7T9aZ3ISAAB6BspNgGX1aS03ByooNwAABAPlJsCykltPS+2n3AAAEBSUmwDznZaqYM4NAADBQLkJsKy2clNR26QaLgcHACDgKDcBFueIVHLv1svBD3LFFAAAAUe5CYL2ScXMuwEAIPAoN0HQfmqKK6YAAAg8yk0QHJ9UTLkBACDQKDdB4DstxUJ+AAAEHOUmCNrXuuG0FAAAgUe5CYL201KV9c2qqnebnAYAgPBGuQmCXlERSo93SJL2lteanAYAgPBGuQmSISm9JUl7yig3AAAEEuUmSNrLzW4n5QYAgECi3ASJb+SG01IAAAQU5SZIhqbESuK0FAAAgWZ6uVm+fLmysrLkcDg0ZcoUbdiwodP9ly1bpnPOOUfR0dHKzMzUbbfdpsbGxiClPXPtIzeHKhtU724xOQ0AAOHL1HKzevVq5eXlafHixdq0aZPGjh2r6dOnq6ys7KT7P/fcc7rzzju1ePFibd++XU8++aRWr16tu+66K8jJuy4pJkpJMa030NxXzno3AAAEiqnlZunSpbrppps0d+5cjRw5UitXrlSvXr301FNPnXT/9evXa9q0afre976nrKwsXXHFFbruuuu+crQnVHDFFAAAgWdauXG73SosLFRubu7xMFarcnNzVVBQcNJjpk6dqsLCQl+Z2bdvn9asWaOrrrrqlO/T1NQkl8vV4WEW3xVTZTWmZQAAINxFmPXGFRUV8ng8Sk1N7bA9NTVVO3bsOOkx3/ve91RRUaELLrhAhmGopaVFN998c6enpZYsWaJf/epXfs1+poYycgMAQMCZPqG4K959913df//9evzxx7Vp0ya9/PLLeuONN3Tvvfee8pgFCxaourra9yguLg5i4o44LQUAQOCZNnKTnJwsm80mp9PZYbvT6VRaWtpJj1m4cKG+//3v60c/+pEkafTo0aqrq9OPf/xj3X333bJaT+xqdrtddrvd/7/AGWgvNweO1svd4lVURLfqlgAAdAumfbtGRUVp4sSJys/P923zer3Kz89XTk7OSY+pr68/ocDYbDZJkmEYgQvrJ2lxDvW2R8jjNXTwKFdMAQAQCKYOHeTl5WnVqlV65plntH37dt1yyy2qq6vT3LlzJUmzZ8/WggULfPvPmDFDK1as0PPPP6/9+/dr7dq1WrhwoWbMmOErOaHMYrFoMKemAAAIKNNOS0nSrFmzVF5erkWLFqm0tFTjxo3TW2+95ZtkXFRU1GGk5p577pHFYtE999yjkpIS9e3bVzNmzNB9991n1q/QZUP69tYnxVXaXVarK80OAwBAGLIY3eF8jh+5XC7Fx8erurpacXFxQX//le/t1QNv7tA3xmbo0evGB/39AQDojrry/c2M1iBrvxx8l5O1bgAACATKTZCdk9Z6A8295bVq9nhNTgMAQPih3ARZv4Roxdoj1OwxtLecScUAAPgb5SbILBaLhqe3jt7sOMKpKQAA/I1yY4IR6a0TobYfMe8+VwAAhCvKjQmGp7WVm1JGbgAA8DfKjQlGtJ2WYuQGAAD/o9yYYFhqrCwWqbymSRW1TWbHAQAgrFBuTBBjj9DApF6SmFQMAIC/UW5M0j6peEcpp6YAAPAnyo1JfJOKGbkBAMCvKDcmYVIxAACBQbkxSftpqT1l3IYBAAB/otyYpF9CtHrbI+T2eLWvvM7sOAAAhA3KjUmsVouGt91Ek0nFAAD4D+XGRO33mNp2mHIDAIC/UG5MNLpfvCRpa0m1yUkAAAgflBsTjfpCufF6DZPTAAAQHig3JhqWGquoCKtqGlt08Fi92XEAAAgLlBsTRdqsGtl2STinpgAA8A/Kjcl8824OVZkbBACAMEG5Mdno/kwqBgDAnyg3JmsfufmsxMWkYgAA/IByY7KhKb1lj7CqtqlFB46yUjEAAGeLcmOyCJtVIzOYVAwAgL9QbkLAmLZTU58eotwAAHC2KDchYHT/BEmM3AAA4A+UmxDQPqn4c1YqBgDgrFFuQsDgvjGKjrSpzu3Rvopas+MAANCtUW5CQITNqlH9WicVby6qMjcMAADdHOUmREwYkChJ2kS5AQDgrFBuQsT4tnKzuajS5CQAAHRvlJsQMWFAgiRpp7NGNY3N5oYBAKAbo9yEiJQ4h/olRMswpE+KuSQcAIAzRbkJIRMGcmoKAICzRbkJIe2npjZRbgAAOGOUmxDSfsXU5uIqGQaL+QEAcCYoNyFkRHqc7BFWVdU3a18FdwgHAOBMUG5CSFSEVWP6t96KYdNBTk0BAHAmKDchhsX8AAA4O5SbEMNifgAAnB3KTYiZMDBBUutiftUNLOYHAEBXUW5CTEqsQ4OSY2QYUuHBY2bHAQCg26HchKDJ2UmSpI/2U24AAOgqyk0I8pWbfZQbAAC6inITgtrLzWcl1aprajE5DQAA3QvlJgT1T+ylfgnRavEa2swl4QAAdAnlJkS1j95s2H/U5CQAAHQvlJsQxaRiAADODOUmRE1pKzebi6vU1OIxOQ0AAN0H5SZEZSfHKLm3Xe4Wrz4prjY7DgAA3QblJkRZLBbf6A3zbgAAOH2UmxDWPu+mYB/lBgCA00W5CWHThvSRJG08UKnGZubdAABwOig3IWxw395KjWudd1N4kLuEAwBwOig3IcxisWjakGRJ0ro9FSanAQCge6DchLgL2srNB5QbAABOC+UmxLWP3GwtqVZ1fbPJaQAACH2ml5vly5crKytLDodDU6ZM0YYNGzrdv6qqSvPmzVN6errsdruGDRumNWvWBClt8KXGOTQkpbcMQyrYx+gNAABfxdRys3r1auXl5Wnx4sXatGmTxo4dq+nTp6usrOyk+7vdbl1++eU6cOCAXnrpJe3cuVOrVq1Sv379gpw8uC5g3g0AAKctwsw3X7p0qW666SbNnTtXkrRy5Uq98cYbeuqpp3TnnXeesP9TTz2lY8eOaf369YqMjJQkZWVldfoeTU1Nampq8j13uVz++wWCZNqQZP1p/QGt38N6NwAAfBXTRm7cbrcKCwuVm5t7PIzVqtzcXBUUFJz0mNdee005OTmaN2+eUlNTNWrUKN1///3yeE69BsySJUsUHx/ve2RmZvr9dwm0KYOSZLVI+yrqVFLVYHYcAABCmmnlpqKiQh6PR6mpqR22p6amqrS09KTH7Nu3Ty+99JI8Ho/WrFmjhQsX6pFHHtGvf/3rU77PggULVF1d7XsUFxf79fcIhjhHpMZmJkiS1u0uNzcMAAAhztTTUl3l9XqVkpKiJ554QjabTRMnTlRJSYl+85vfaPHixSc9xm63y263Bzmp/108rK82F1Xp3Z3lmnXeALPjAAAQskwbuUlOTpbNZpPT6eyw3el0Ki0t7aTHpKena9iwYbLZbL5tI0aMUGlpqdxud0Dzmu2Sc1IkSet2V6jZ4zU5DQAAocu0chMVFaWJEycqPz/ft83r9So/P185OTknPWbatGnas2ePvN7jX+67du1Senq6oqKiAp7ZTGP6xatPTJRqmlq4FQMAAJ0w9VLwvLw8rVq1Ss8884y2b9+uW265RXV1db6rp2bPnq0FCxb49r/lllt07Ngx3Xrrrdq1a5feeOMN3X///Zo3b55Zv0LQWK0WXTSsryTpnZ0nv1QeAACYPOdm1qxZKi8v16JFi1RaWqpx48bprbfe8k0yLioqktV6vH9lZmbq7bff1m233aYxY8aoX79+uvXWW3XHHXeY9SsE1SXn9NUrm0v03s5yLbhyhNlxAAAISRbDMAyzQwSTy+VSfHy8qqurFRcXZ3acLqmsc2vir9fKa0jr7/yaMhKizY4EAEBQdOX72/TbL+D0JcZEaVzbJeHv7eKScAAAToZy0820XzX1zg7m3QAAcDKUm27m0rZy88GeCjW1nHplZgAAeirKTTdzbkac+sbaVef26MN9x8yOAwBAyKHcdDNWq0W5I1qvJlu77eS3qQAAoCej3HRDV4xsLzdOeb096mI3AAC+EuWmG8oZ3EcxUTY5XU3aWlJtdhwAAEIK5aYbckTadPE5rasVr93m/Iq9AQDoWSg33dTlbaem/sm8GwAAOqDcdFNfOydVNqtFu5y1OlBRZ3YcAABCBuWmm4rvFakp2UmSODUFAMAXUW66MU5NAQBwIspNNzb93DRJ0scHK+V0NZqcBgCA0EC56cYyEqI1YUCCDEN6c+sRs+MAABASzqjcbN++XYsXL9bXvvY1DR48WOnp6RozZozmzJmj5557Tk1NTf7OiVO4ekyGJOkNyg0AAJK6WG42bdqk3NxcjR8/XuvWrdOUKVM0f/583XvvvbrhhhtkGIbuvvtuZWRk6MEHH6TkBMFVo4+fmiqt5tQUAAARXdn5mmuu0e23366XXnpJCQkJp9yvoKBAv/vd7/TII4/orrvuOtuM6ER6fLQmDkxU4cFKvfnZEc2dlm12JAAATNWlcrNr1y5FRkZ+5X45OTnKyclRc3PzGQfD6btqdLoKD1ZqzVbKDQAAXTotdTrFRpLq6+u7tD/OTvupqY0HODUFAMAZXy112WWXqaSk5ITtGzZs0Lhx484mE7qo/dSUJK1hYjEAoIc743LjcDg0ZswYrV69WpLk9Xr1y1/+UhdccIGuuuoqvwXE6ZkxJl2S9PctJxZOAAB6ki7NufmiN954Q8uXL9eNN96ov//97zpw4IAOHjyo119/XVdccYU/M+I0/NfYDN37xnZ9cqha+8prNahvb7MjAQBgijMuN5I0b948HTp0SA8++KAiIiL07rvvaurUqf7Khi5I7m3XRUOT9c7Ocr265bDyLh9mdiQAAExxxqelKisrdc0112jFihX6wx/+oGuvvVZXXHGFHn/8cX/mQxfMHN9PkvTq5hIZhmFyGgAAzHHG5WbUqFFyOp3avHmzbrrpJv3lL3/Rk08+qYULF+rqq6/2Z0acpstHpqpXlE1Fx+q1qajK7DgAAJjijMvNzTffrPfff1/Z2cfXVZk1a5Y++eQTud1uv4RD1/SKitDX226m+epmJhYDAHqmMy43CxculNV64uH9+/fX2rVrzyoUzlz7qanXPz0sd4vX5DQAAARfl8pNUVFRl374ydbBQWBNHdxHfWPtqqxv1js7y8yOAwBA0HWp3Jx33nn6yU9+oo0bN55yn+rqaq1atUqjRo3S3/72t7MOiK6JsFn17bbRmxc2FpucBgCA4OvSpeDbtm3Tfffdp8svv1wOh0MTJ05URkaGHA6HKisrtW3bNn3++eeaMGGCHnroIRbzM8l3J2XqD+/v0zs7y+R0NSo1zmF2JAAAgqZLIzd9+vTR0qVLdeTIET322GMaOnSoKioqtHv3bknS9ddfr8LCQhUUFFBsTDQkpbcmDUyU15D+tumQ2XEAAAgqi9HFBVH27dun7OxsWSyWQGUKKJfLpfj4eFVXVysuLs7sOAHzwsZi/eJvnyo7OUb//v8u7rb/ewEAIHXt+7vLV0sNHTpU5eXlvuezZs2S0+nsekoE1NVj0hUTZdP+ijptPFBpdhwAAIKmy+XmywM9a9asUV1dnd8CwT9i7BH6rzEZkqTVTCwGAPQgZ7zODULftef1lyS9sfWwquubTU4DAEBwdLncWCyWE+ZvMJ8jNE0YkKjhabFqbPbqxUJGbwAAPUOX7wpuGIZ+8IMfyG63S5IaGxt18803KyYmpsN+L7/8sn8S4oxZLBZ9P2eg7n7lM/31oyLdOC1bVitFFAAQ3rpcbubMmdPh+Q033OC3MPC/meP66YE1O7S/ok4f7K3QhUP7mh0JAICA6nK5efrppwORAwESY4/QNRP760/rD+jPBQcpNwCAsMeE4h7ghvMHSJLytztVUtVgchoAAAKLctMDDEmJ1dTBfeQ1pOc+Omh2HAAAAopy00N8//yBklrXvGlq8ZicBgCAwKHc9BCXj0xVapxdFbVuvfVZqdlxAAAIGMpNDxFhs+r6Ka2jN3/8z/4TVpoGACBcUG56kBvOHyh7hFVbS6r10f5jZscBACAgKDc9SFJMlL47qfWWDKve32dyGgAAAoNy08P88IJBslik/B1l2lNWY3YcAAD8jnLTw2Qnx+jyEamSWufeAAAQbig3PdCPLxokSXp5c4nKa5pMTgMAgH9RbnqgiQMTNX5AgtwtXj1bcMDsOAAA+BXlpgeyWCz68YWtozd//vCg6t0tJicCAMB/KDc91BXnpmlgn16qqm/Wcx8VmR0HAAC/odz0UDarRfMuGSJJWvnePjU2c0sGAEB4oNz0YN+a0E/9EqJVUduk5zcwegMACA+Umx4s0mbV/7t0sCRpxXt7Gb0BAIQFyk0P952J/ZUe75DT1aQXCw+ZHQcAgLNGuenh7BE23Xxx6+jNynf3yt3iNTkRAABnJyTKzfLly5WVlSWHw6EpU6Zow4YNp3Xc888/L4vFopkzZwY2YJibdV6mUmLtKqlq0MubGL0BAHRvppeb1atXKy8vT4sXL9amTZs0duxYTZ8+XWVlZZ0ed+DAAf3P//yPLrzwwiAlDV+OSJtv1eLf/3uPmlqYewMA6L5MLzdLly7VTTfdpLlz52rkyJFauXKlevXqpaeeeuqUx3g8Hl1//fX61a9+pUGDBgUxbfi64fyBSo1rHb1h3RsAQHdmarlxu90qLCxUbm6ub5vValVubq4KCgpOedz//u//KiUlRT/84Q+/8j2amprkcrk6PHAiR6RNt142TJL02L/3qLaJVYsBAN2TqeWmoqJCHo9HqampHbanpqaqtLT0pMesW7dOTz75pFatWnVa77FkyRLFx8f7HpmZmWedO1x9d1J/ZSfH6GidW0+t447hAIDuyfTTUl1RU1Oj73//+1q1apWSk5NP65gFCxaourra9yguLg5wyu4r0mZV3uWtozdPvL9Px+rcJicCAKDrIsx88+TkZNlsNjmdzg7bnU6n0tLSTth/7969OnDggGbMmOHb5vW2XrocERGhnTt3avDgwR2OsdvtstvtAUgfnq4ena6V7+3V54ddevydPbrnv0aaHQkAgC4xdeQmKipKEydOVH5+vm+b1+tVfn6+cnJyTth/+PDh2rp1q7Zs2eJ7fOMb39Cll16qLVu2cMrJD6xWi37x9eGSWu8YXlLVYHIiAAC6xtSRG0nKy8vTnDlzNGnSJE2ePFnLli1TXV2d5s6dK0maPXu2+vXrpyVLlsjhcGjUqFEdjk9ISJCkE7bjzF00NFnnD0rSh/uO6eG3d+q3s8aZHQkAgNNmermZNWuWysvLtWjRIpWWlmrcuHF66623fJOMi4qKZLV2q6lB3Z7FYtHdV43UN5av0yubSzQ7Z6DGD0g0OxYAAKfFYhiGYXaIYHK5XIqPj1d1dbXi4uLMjhPSbn/xE71YeEjjByTo5VumymKxmB0JANBDdeX7myERnNLt089RryibNhdV6bVPDpsdBwCA00K5wSmlxDk079IhkqQH3tyhBje3ZQAAhD7KDTr1wwuy1S8hWkeqG/XE+/vMjgMAwFei3KBTjkibFlzVemn4ivf2qPhYvcmJAADoHOUGX+nq0emakp2kxmavfvna5+phc9ABAN0M5QZfyWKx6L5vjVKkzaL8HWX65zbnVx8EAIBJKDc4LUNSYnXThYMkSb967XPVcddwAECIotzgtP3sa0PVPzFah6sb9Wj+brPjAABwUpQbnLboKJv+95vnSpL+uG6/dpS6TE4EAMCJKDfokq8NT9XXz02Tx2vojr9tVYvHa3YkAAA6oNygy375jXMV64jQJ8VVenLdfrPjAADQAeUGXZYW79DCq0dKkh5Zu0t7y2tNTgQAwHGUG5yR707qr4uG9ZW7xas7XvpUHi9r3wAAQgPlBmfEYrFoybdHKybKpo8PVuqZ9QfMjgQAgCTKDc5Cv4RoLbhqhCTpobd3aH9FncmJAACg3OAsfW/yAE0b0keNzV7NX71FzVw9BQAwGeUGZ8Vqteg33xmruLarpx779x6zIwEAejjKDc5aRkK0fv2t0ZKkx97Zo8KDlSYnAgD0ZJQb+MU3xmZo5rgMebyG8l7Ywr2nAACmodzAb371zVHKiHfo4NF6/eofn5sdBwDQQ1Fu4Dfx0ZF65NpxslikFz4+pFc2HzI7EgCgB6LcwK9yBvfRz782VJJ09yufaU8ZqxcDAIKLcgO/+/llQ5UzqI/q3R7N++smNbg9ZkcCAPQglBv4nc1q0e+uG6fk3nbtdNYw/wYAEFSUGwRESqxDj/536/yb5zcW68WPi82OBADoISg3CJipQ5I1/7JhkqS7X/1Mm4tY/wYAEHiUGwTUz742RJePTJW7xaufPFsop6vR7EgAgDBHuUFAWa0W/XbWOA1L7a2ymib95NlCNTYzwRgAEDiUGwRcb3uEVs2epPjoSG0prtI9r34mwzDMjgUACFOUGwTFwD4xeux742W1SC8VHtLTHxwwOxIAIExRbhA0Fw7tq7uuGiFJ+vUb25S/3WlyIgBAOKLcIKh+eEG2rp3UX15D+ulzm/XpoSqzIwEAwgzlBkFlsVh037dG68KhyWpo9ujGP32s4mP1ZscCAIQRyg2CLtJm1ePXT9CI9DhV1DbpB09vUFW92+xYAIAwQbmBKWIdkXr6B+cpPd6hveV1uvFPG1XX1GJ2LABAGKDcwDRp8Q79ae5kxUdHalNRFWvgAAD8gnIDU52TFqs/zT1PvaJsWrenQj//v81q8XjNjgUA6MYoNzDd+AGJ+uOcSYqKsOqf25z6xUufyutlkT8AwJmh3CAkTB2crBXXT1CE1aKXN5do0WusYgwAODOUG4SMy0akaumscbJYpL98WKTFr31OwQEAdBnlBiHlG2Mz9OA1Y2SxSH8uOKh7Xv2MU1QAgC6h3CDkXDspUw9/Z6wsFumvHxXprle2UnAAAKeNcoOQdM3E/vrtteNktUjPbyzWL/72qTwUHADAaaDcIGTNHN9Pv501zncn8dtWb5G7hcvEAQCdo9wgpH1zXD89et14RVgteu2Tw7rpzx+r3s1KxgCAU6PcIOT915gMrZozSY5Iq97bVa4b/vgR96ICAJwS5QbdwqXnpOivPzrfd6uGa/9QoNLqRrNjAQBCEOUG3cbEgYl68eYcpcbZtctZq2tWrNduZ43ZsQAAIYZyg25lWGqsXrp5qrKTY1RS1aBvP75e63ZXmB0LABBCKDfodjKTeulvt0zV5Kwk1TS16AdPb9DzG4rMjgUACBGUG3RLSTFRevZHk/Wt8f3U4jV058tbteTN7Sz2BwCg3KD7skfYtPTasZqfO1SS9If39uknfylUTWOzyckAAGai3KBbs1gsmp87TMtmjVNUhFVrtzk1c/kH2lNWa3Y0AIBJKDcICzPH99OLP8lRWpxDe8vrNHP5B1q7zWl2LACACSg3CBtjMxP0j59doMnZSaptatFNf/5YS9fu4p5UANDDUG4QVvrG2vXXH03RD6ZmSZIezd+t2U99pLIaFvwDgJ6CcoOwE2mz6pffOFdLrx2r6EibPthzVFf97j/6z+5ys6MBAIKAcoOw9e0J/fWPn12g4Wmxqqh1a/ZTG/Sbt3eoxcOdxQEgnIVEuVm+fLmysrLkcDg0ZcoUbdiw4ZT7rlq1ShdeeKESExOVmJio3NzcTvdHzzYkpbdenTdN108ZIMOQlr+zV//9xIc6VFlvdjQAQICYXm5Wr16tvLw8LV68WJs2bdLYsWM1ffp0lZWVnXT/d999V9ddd53eeecdFRQUKDMzU1dccYVKSkqCnBzdhSPSpvu+NVrLvzdBsfYIfXywUl9f9h+98HGxDIPJxgAQbiyGyf+6T5kyReedd54ee+wxSZLX61VmZqZ+9rOf6c477/zK4z0ejxITE/XYY49p9uzZX7m/y+VSfHy8qqurFRcXd9b50b0UHa3X/NWbtamoSpKUOyJVS749Wn1j7eYGAwB0qivf36aO3LjdbhUWFio3N9e3zWq1Kjc3VwUFBaf1M+rr69Xc3KykpKSTvt7U1CSXy9XhgZ5rQJ9eevHmqfrF189RpM2if213avqy9/Xm1iNmRwMA+Imp5aaiokIej0epqakdtqempqq0tPS0fsYdd9yhjIyMDgXpi5YsWaL4+HjfIzMz86xzo3uzWS36f5cM0Ws/bZ1sfKzOrVv+ukm3Pr9ZR2ubzI4HADhLps+5ORsPPPCAnn/+eb3yyityOBwn3WfBggWqrq72PYqLi4OcEqFqRHqcXvvpBZp36WBZLdLftxzWZUvf00uFh5iLAwDdmKnlJjk5WTabTU5nx2XynU6n0tLSOj324Ycf1gMPPKB//vOfGjNmzCn3s9vtiouL6/AA2kVFWHX79OF6+f9N0/C0WFXVN+t/XvxENzz5kQ5U1JkdDwBwBkwtN1FRUZo4caLy8/N927xer/Lz85WTk3PK4x566CHde++9euuttzRp0qRgREWYG9d264Y7vj5c9girPthzVNOXva/H390jdwvr4gBAd2L6aam8vDytWrVKzzzzjLZv365bbrlFdXV1mjt3riRp9uzZWrBggW//Bx98UAsXLtRTTz2lrKwslZaWqrS0VLW13AUaZyfSZtUtlwzWP2+7SBcMSVZTi1cPvbVTV/7ufb2/i9WNAaC7ML3czJo1Sw8//LAWLVqkcePGacuWLXrrrbd8k4yLiop05MjxK1lWrFght9ut73znO0pPT/c9Hn74YbN+BYSZgX1i9OwPJ+uR745Vcu8o7S2v0+ynNujHf/5YxcdY/A8AQp3p69wEG+vcoCuqG5r1u3/t1jMFB+TxGoqKsOrmiwbplkuGKDrKZnY8AOgxuvL9TbkBTsNuZ41++Y/P9cGeo5KktDiH8i4fpmsm9pfNajE5HQCEP8pNJyg3OFOGYejtz0t17+vbVVLVIEkaltpbd3x9uL42PEUWCyUHAAKFctMJyg3OVmOzR3/58KB+/+89qm5oliRNzk7SgiuHa/yARJPTAUB4otx0gnIDf6luaNaKd/fq6Q/2q6ntcvErR6Vpfu4wnZMWa3I6AAgvlJtOUG7gb4erGvTbtbv00qZDMgzJYpGuGp2uWy8bqmGplBwA8AfKTScoNwiUnaU1+l3+Lq3Z2npfNEoOAPgP5aYTlBsE2o5Slx7N331CyZl3yRCNzOD/cwBwJig3naDcIFi+XHIk6aJhfXXzRYOUM7gPV1cBQBdQbjpBuUGw7Sh1afk7e/XGp4flbfuvbXS/eP3k4kH6+rlpirCZvlA4AIQ8yk0nKDcwS/Gxev3xP/u0+uNiNTa3Xl01IKmX5k7L0ncm9lesI9LkhAAQuig3naDcwGzH6tz6c8EBPbP+gCrrW9fJiYmy6ZqJ/TU7J0tDUnqbnBAAQg/lphOUG4SKBrdHLxUW65mCg9pTdvyu9hcOTdacnCxdOjyFWzsAQBvKTScoNwg1hmHogz1H9UzBAf1ru1Pt/0VmJkXr+ikDdc2E/uobazc3JACYjHLTCcoNQlnxsXr95cODen5jse/WDhFWi3JHpGrW5ExdNLQvozkAeiTKTScoN+gOGtwevfZJif5vQ7G2FFf5tqfHO/Tdif313UmZykzqZV5AAAgyyk0nKDfobnaUurR6Y7Fe2VyiqrYJyBaLNG1wsr41vp+mj0pTb3uEySkBILAoN52g3KC7amz2aO02p1ZvLNa6PRW+7Y5Iqy4fmaaZ4zJ00bC+imTdHABhiHLTCcoNwkHxsXq9srlEr24u0b6KOt/2pJgoXT06XTPHZ2jCgERWQQYQNig3naDcIJwYhqFPD1Xr1S0l+scnh1VR6/a9lhHv0NdHpeuq0WmaMCBRViYiA+jGKDedoNwgXLV4vPpg71G9urlEb39eqnq3x/daSqxdV45K05Wj03VeVhJXXAHodig3naDcoCdobPbo/V3levOzUv1rm1M1TS2+15J7R+mKc9N0xchUnT+ojxyRNhOTAsDpodx0gnKDnqapxaMP9lToza2l+uc2p2/9HEnqFWXTBUOSlTsiVZcOT2GxQAAhi3LTCcoNerJmj1cFe4/qrc9Llb/dKaerqcPrYzMTlDs8RZeNSNWI9FgmJAMIGZSbTlBugFaGYejzwy7lby9T/g6nPj1U3eH11Di7LhzaVxcOTda0IclK7s2oDgDzUG46QbkBTs7patQ7O8r0r+1lWrenXI3N3g6vn5sRpwuH9tVFQ5M1MStR9gjm6gAIHspNJyg3wFdrbPao8GCl3t9drv/sqtC2I64Or0dH2jRlUJKmDU7WlEFJGpkepwgWDwQQQJSbTlBugK4rr2nSB3sqWsvO7gqV13ScqxNrj9B52Umakp2k8wf10bkZlB0A/kW56QTlBjg7hmFop7NG63ZX6MN9R/XR/mOqaWzpsE9ve4QmZSXq/EF9dF5Wkkb1i+M0FoCzQrnpBOUG8C+P19D2Iy59uO+oPtx3TBv2H5XrS2UnymbV6P7xmjgwURMGJGrCwASlxDpMSgygO6LcdIJyAwRWe9n5aP8xfbjvqDYdrNTROvcJ+2UmRWvigMTWwjMwUeekxnIqC8ApUW46QbkBgsswDB08Wq/Cg5UqLKrUpoOV2ums0Zf/5XFEWnVuRrzG9G9/JCi7Twz3xAIgiXLTKcoNYD5XY7O2FFWp8GClNhVVanNRlWqbWk7Yr7c9QqP6xWls/wSN7h+vsf0T1D8xmsUFgR6IctMJyg0QerxeQ/sq6vTpoSp9eqhaW0uq9fnh6hPW2pGkhF6RGpkep5HpcRrR9hiS0ltREZzSAsIZ5aYTlBuge2jxeLW7rFZbD1Xr05LW0rP9iEvNnhP/yYq0WTS4b2+NzOhYepJiokxIDiAQKDedoNwA3VdTi0e7Smu1/YhL2464fH9++VL0dqlxdg1LjdXQlFgNTe2tYam9NSQlVvHRkUFODuBsUW46QbkBwothGCqpatD2IzWtZeewS9tLXTp4tP6Ux6TG2X2Fx1d8UmIV34vSA4Qqyk0nKDdAz1Db1KKdpTXaU1ajXc5a7S6r1W5njY5UN57ymL6xdg1N6a3s5JgOj8ykXorkMnXAVJSbTlBugJ6tprFZu8tqtcdZq13OGl/pOdxJ6bFZLRqQ1EvZyTHK6hOj7L4xGtRWfNLiHFyuDgRBV76/I4KUCQBCQqwjsnWV5AGJHbbXNDZrT1mt9pXXaX9F62NfRZ0OVNSpodnj2/ZljkirsvrEaEBSr9ZHn17KTOylzKRe6p8YLUckt50Ago2RGwDohGEYcrqatK+itrXglNfpwNHW4lN0tF4t3s7/CU2Ns2tAUmvZyUzsWIBSYu2M+gCnidNSnaDcAPCXFo9XhyobtL+iTsWV9So6Wq+iY/UqrmxQ8bH6ky5M+EVREVb1T4hWRkK0MhIcbX9GKyP++HNGfoBWnJYCgCCIsFmVlRyjrOSYE14zDEOV9c0qPtZaeIqO1etQZVv5OdagkqoGuVu82td2+utU+sREKT3B0VZ4otWvrQClJzjULyFayb3tsjH6A3RAuQGAALBYLEqKiVJSTJTGZiac8HqLx6sj1Y06VNmgw1Vtj+oGHa5q1OGq1vJT7/boaJ1bR+vc+qzEddL3sVpar/JKjXO0PexKjXUoNf7487Q4h+KjI7ltBXoMyg0AmCDCZm2dh5PU66SvG4YhV0OLSqoadKS6oa3wtBafI20lqNTVKI+3dU6Q09UkqfqU72ePsPrKTkqcQ2ntRSjOob697UqOtSu5t10J0ZHMA0K3R7kBgBBksVgU3ytS8b0iNTLj5PMLPF5DR2tbi02pq1FOV6PKXI1tf2+Ss21bZX2zmlq8vtNjnYmwto44JfsKT5T69rarb1v5ad3e+npiryhOiSEkUW4AoJuyWS1KiXMoJc6h0Yo/5X6NzR6V17SXndYiVNZWfEpdjaqodauitklV9c1q8Roqq2lSWU2TdKTz97dapKSYtgIUa1efmCglxkQpqVfbn194JPaKUmKvSEWwGCKCgHIDAGHOEWnr9BRYO3eLV8fq3CqvaVJFbZPKa1v/rKhpLT/HH25V1rvlNeTbtqO05rSyxDkiTig9STFfKENfKEYJ0ZGKi45kdAhdRrkBAEhqvTQ9Ld6htHjHV+7b4mkrQm1lp6KmSZX1rZOfK+vcOlbn7vC8qqFZhiG5GlvkamzRgU7u/fVlsY4IJfSKVHx0pBKio1pP10VHKiG67c+25/HRUcf36xWp6Egbk6h7KMoNAKDLImxW3ymx0+HxGqpuaNaxtuLTXn6OfaEMHatv+3u9W8dq3apzeyRJNY0tqmlsUbEaupQx0mbpWHjaylBcdKTiHBGKdUQqLrr1z9j2547jz1ljqPui3AAAAs5mPX5p/Olq9nhV3dCs6oZmVdU3y9XQrKoGt6rqv7yt/bnbt3+zx1Czx/CdNjsTUTarYh0RiotuLz8RirVHfmlb23NHhOIckb7nMfbW/e0RVkaPTEC5AQCEpEib1XeFVlcYhqF6t8dXgFoLj7vD89bRoGa52v5sHx1yNTSr1t0iw5DcHq9vnaEzZbNaFBNlU297a+GJsUeod9uj9e+2Dts7vn7iccw/Oj2UGwBAWLFYLL5CkJEQ3eXjvV5Dte7jZaemQwFqbps31LEQdShIjc2qbzul5vEavnlG/hAdaTtpKeoVZWt7HP97dNSJ26OjbIr5wt97RdnCcm4S5QYAgC+wWi2Kc0QqzhGpfmdQjqTWUlPvblFdk0e1TS2qa2pRbdujzvfc49t+4uttx7lbVNvY4rtBa0OzRw3NHlXU+u/3tVhaS9OXS9AJZSmydTQpOsqmXpFtr9lbX3NE2tp+RoSiI23q3XZVnFkoNwAA+JnNammbfxN51j/LMAw1tXhPLD1fKEp1Ta2lp97donq3R/VNHtU3e9TQ9rzOffzvrY8WNTZ7236+fNulMz8F90Vj+sfrtZ9e4JefdSYoNwAAhDCLxSJHZOvoSJ/e/vu5Xq+hhmaP6twtavCVnuMFqeFLz+vbClKd77Xj2xubj//Z0OxRTJS59YJyAwBAD2S1Hp+b5G+GYfj9Z3YF62ADAAC/MnuCMuUGAACElZAoN8uXL1dWVpYcDoemTJmiDRs2dLr/iy++qOHDh8vhcGj06NFas2ZNkJICAIBQZ3q5Wb16tfLy8rR48WJt2rRJY8eO1fTp01VWVnbS/devX6/rrrtOP/zhD7V582bNnDlTM2fO1GeffRbk5AAAIBRZDJNn/UyZMkXnnXeeHnvsMUmS1+tVZmamfvazn+nOO+88Yf9Zs2aprq5Or7/+um/b+eefr3HjxmnlypVf+X4ul0vx8fGqrq5WXFyc/34RAAAQMF35/jZ15MbtdquwsFC5ubm+bVarVbm5uSooKDjpMQUFBR32l6Tp06efcv+mpia5XK4ODwAAEL5MLTcVFRXyeDxKTU3tsD01NVWlpaUnPaa0tLRL+y9ZskTx8fG+R2Zmpn/CAwCAkGT6nJtAW7Bggaqrq32P4uJisyMBAIAAMnURv+TkZNlsNjmdzg7bnU6n0tLSTnpMWlpal/a32+2y27t2R1kAANB9mTpyExUVpYkTJyo/P9+3zev1Kj8/Xzk5OSc9Jicnp8P+krR27dpT7g8AAHoW02+/kJeXpzlz5mjSpEmaPHmyli1bprq6Os2dO1eSNHv2bPXr109LliyRJN166626+OKL9cgjj+jqq6/W888/r48//lhPPPGEmb8GAAAIEaaXm1mzZqm8vFyLFi1SaWmpxo0bp7feess3abioqEhW6/EBpqlTp+q5557TPffco7vuuktDhw7Vq6++qlGjRpn1KwAAgBBi+jo3wcY6NwAAdD/dZp0bAAAAfzP9tFSwtQ9UsZgfAADdR/v39umccOpx5aampkaSWMwPAIBuqKamRvHx8Z3u0+Pm3Hi9Xh0+fFixsbGyWCx+/dkul0uZmZkqLi5mPk8A8TkHB59zcPA5Bw+fdXAE6nM2DEM1NTXKyMjocKHRyfS4kRur1ar+/fsH9D3i4uL4DycI+JyDg885OPicg4fPOjgC8Tl/1YhNOyYUAwCAsEK5AQAAYYVy40d2u12LFy/mXlYBxuccHHzOwcHnHDx81sERCp9zj5tQDAAAwhsjNwAAIKxQbgAAQFih3AAAgLBCuQEAAGGFcuMny5cvV1ZWlhwOh6ZMmaINGzaYHSnsLFmyROedd55iY2OVkpKimTNnaufOnWbHCmsPPPCALBaL5s+fb3aUsFRSUqIbbrhBffr0UXR0tEaPHq2PP/7Y7FhhxePxaOHChcrOzlZ0dLQGDx6se++997TuT4RTe//99zVjxgxlZGTIYrHo1Vdf7fC6YRhatGiR0tPTFR0drdzcXO3evTto+Sg3frB69Wrl5eVp8eLF2rRpk8aOHavp06errKzM7Ghh5b333tO8efP04Ycfau3atWpubtYVV1yhuro6s6OFpY0bN+oPf/iDxowZY3aUsFRZWalp06YpMjJSb775prZt26ZHHnlEiYmJZkcLKw8++KBWrFihxx57TNu3b9eDDz6ohx56SL///e/Njtat1dXVaezYsVq+fPlJX3/ooYf06KOPauXKlfroo48UExOj6dOnq7GxMTgBDZy1yZMnG/PmzfM993g8RkZGhrFkyRITU4W/srIyQ5Lx3nvvmR0l7NTU1BhDhw411q5da1x88cXGrbfeanaksHPHHXcYF1xwgdkxwt7VV19t3HjjjR22ffvb3zauv/56kxKFH0nGK6+84nvu9XqNtLQ04ze/+Y1vW1VVlWG3243/+7//C0omRm7OktvtVmFhoXJzc33brFarcnNzVVBQYGKy8FddXS1JSkpKMjlJ+Jk3b56uvvrqDv+/hn+99tprmjRpkr773e8qJSVF48eP16pVq8yOFXamTp2q/Px87dq1S5L0ySefaN26dbryyitNTha+9u/fr9LS0g7/fsTHx2vKlClB+17scTfO9LeKigp5PB6lpqZ22J6amqodO3aYlCr8eb1ezZ8/X9OmTdOoUaPMjhNWnn/+eW3atEkbN240O0pY27dvn1asWKG8vDzddddd2rhxo37+858rKipKc+bMMTte2Ljzzjvlcrk0fPhw2Ww2eTwe3Xfffbr++uvNjha2SktLJemk34vtrwUa5Qbd0rx58/TZZ59p3bp1ZkcJK8XFxbr11lu1du1aORwOs+OENa/Xq0mTJun++++XJI0fP16fffaZVq5cSbnxoxdeeEF//etf9dxzz+ncc8/Vli1bNH/+fGVkZPA5hzFOS52l5ORk2Ww2OZ3ODtudTqfS0tJMShXefvrTn+r111/XO++8o/79+5sdJ6wUFhaqrKxMEyZMUEREhCIiIvTee+/p0UcfVUREhDwej9kRw0Z6erpGjhzZYduIESNUVFRkUqLwdPvtt+vOO+/Uf//3f2v06NH6/ve/r9tuu01LliwxO1rYav/uM/N7kXJzlqKiojRx4kTl5+f7tnm9XuXn5ysnJ8fEZOHHMAz99Kc/1SuvvKJ///vfys7ONjtS2Lnsssu0detWbdmyxfeYNGmSrr/+em3ZskU2m83siGFj2rRpJyxlsGvXLg0cONCkROGpvr5eVmvHrzqbzSav12tSovCXnZ2ttLS0Dt+LLpdLH330UdC+Fzkt5Qd5eXmaM2eOJk2apMmTJ2vZsmWqq6vT3LlzzY4WVubNm6fnnntOf//73xUbG+s7dxsfH6/o6GiT04WH2NjYE+YwxcTEqE+fPsxt8rPbbrtNU6dO1f33369rr71WGzZs0BNPPKEnnnjC7GhhZcaMGbrvvvs0YMAAnXvuudq8ebOWLl2qG2+80exo3Vptba327Nnje75//35t2bJFSUlJGjBggObPn69f//rXGjp0qLKzs7Vw4UJlZGRo5syZwQkYlGuyeoDf//73xoABA4yoqChj8uTJxocffmh2pLAj6aSPp59+2uxoYY1LwQPnH//4hzFq1CjDbrcbw4cPN5544gmzI4Udl8tl3HrrrcaAAQMMh8NhDBo0yLj77ruNpqYms6N1a++8885J/z2eM2eOYRitl4MvXLjQSE1NNex2u3HZZZcZO3fuDFo+i2GwTCMAAAgfzLkBAABhhXIDAADCCuUGAACEFcoNAAAIK5QbAAAQVig3AAAgrFBuAABAWKHcAACAsEK5AQAAYYVyAwAAwgrlBgAAhBXKDYBur7y8XGlpabr//vt929avX6+oqCjl5+ebmAyAGbhxJoCwsGbNGs2cOVPr16/XOeeco3Hjxumb3/ymli5danY0AEFGuQEQNubNm6d//etfmjRpkrZu3aqNGzfKbrebHQtAkFFuAISNhoYGjRo1SsXFxSosLNTo0aPNjgTABMy5ARA29u7dq8OHD8vr9erAgQNmxwFgEkZuAIQFt9utyZMna9y4cTrnnHO0bNkybd26VSkpKWZHAxBklBsAYeH222/XSy+9pE8++US9e/fWxRdfrPj4eL3++utmRwMQZJyWAtDtvfvuu1q2bJmeffZZxcXFyWq16tlnn9V//vMfrVixwux4AIKMkRsAABBWGLkBAABhhXIDAADCCuUGAACEFcoNAAAIK5QbAAAQVig3AAAgrFBuAABAWKHcAACAsEK5AQAAYYVyAwAAwgrlBgAAhJX/H4IRUhlkABAbAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "initial_conditions = [1, -1.5886] \n",
    "\n",
    "def f(x): \n",
    "    return (1+4/3*x**(3/2))\n",
    "\n",
    "def F(x_range):\n",
    "    solution = solve_differential_equation(DGL_system, x_range, initial_conditions, plot_result=False)\n",
    "    return solution.t, solution.y[0]*f(solution.t)\n",
    "\n",
    "x, Fx = F(x_range)\n",
    "\n",
    "plt.plot(x, Fx)\n",
    "plt.ylabel('F(x)')\n",
    "plt.xlabel('x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6e5a2d7-9f39-4a22-8486-d79c8069d7d2",
   "metadata": {},
   "source": [
    "We see $\\tilde{F}(x) \\rightarrow 1$ for $x \\rightarrow 0$, which indicates the region of no screening."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65eb49bb-4660-4dd6-bbbf-b96d6bb2ac9e",
   "metadata": {},
   "source": [
    "d) We consider the effective potential $V_{eff}(r) = -\\frac{Z}{r}F(r) + \\frac{\\ell(\\ell + 1)}{2 r^2}$. The potential has a bound state if it has a minimum with negative energy. We see that for a given $\\ell$, the potential can reach negative values if $Z$ gets large enough. We determine the minimal $Z$ in order for $V_{eff}$ to take negative values and obtian:\n",
    "$$\n",
    "Z(\\ell) = \\frac{4}{3 \\pi} \\Bigl ( \\frac{\\ell(\\ell + 1)}{\\tilde{F}(x) x} \\Bigr) ^{3/2}.\n",
    "$$\n",
    "To evaluate this expression, we numerically determine the maximum of $\\tilde{F}(x) x$, using the solution of $\\tilde{F}$ we found above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a4d7e0b9-e03b-4150-8186-21e9683701e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQ+ElEQVR4nO3deVjU1f4H8PcszAz7KiCIAoILoqAgihsuKJYt3rLMXy5hWd30ppK2p93bQpoaLS5pmTfLslXLjFIUVxQFNXdURBBklx1mYGZ+f4xMcUUTnZnvMPN+Pc88I19mhg+TNe/O+ZxzRFqtVgsiIiIiCyEWugAiIiIiQ2K4ISIiIovCcENEREQWheGGiIiILArDDREREVkUhhsiIiKyKAw3REREZFGkQhdgahqNBgUFBXB0dIRIJBK6HCIiIroFWq0W1dXV8PHxgVh887EZqws3BQUF8PPzE7oMIiIiug15eXno1KnTTR9jdeHG0dERgO7NcXJyErgaIiIiuhVVVVXw8/PTf47fjNWFm+apKCcnJ4YbIiKiduZWWkrYUExEREQWheGGiIiILArDDREREVkUq+u5ISIiaiu1Wo3Gxkahy7B4Mpnsb5d53wqGGyIiohvQarUoLCxERUWF0KVYBbFYjICAAMhksjt6HYYbIiKiG2gONp6enrCzs+Pmr0bUvMnulStX0Llz5zt6rxluiIiIWqFWq/XBxt3dXehyrEKHDh1QUFCApqYm2NjY3PbrsKGYiIioFc09NnZ2dgJXYj2ap6PUavUdvQ7DDRER0U1wKsp0DPVeM9wQERGRRTGLcLN8+XL4+/tDoVBgwIABSE9Pv+Fj161bB5FI1OKmUChMWC0RERGZM8HDzcaNG5GQkICFCxciMzMTYWFhiIuLQ3Fx8Q2f4+TkhCtXruhvly5dMmHFREREZM4EDzfLli3DjBkzEB8fj5CQEKxatQp2dnZYu3btDZ8jEong7e2tv3l5eZmwYiIiIjJngoYblUqFjIwMxMbG6q+JxWLExsYiLS3ths+rqalBly5d4Ofnh/vvvx8nT5684WOVSiWqqqpa3IiIiMhyCRpuSktLoVarrxt58fLyQmFhYavP6d69O9auXYvNmzfjiy++gEajwaBBg3D58uVWH5+YmAhnZ2f9zc/Pz+C/BxERkbkoKSmBt7c33n77bf21/fv3QyaTISUlpdXn5OTkQCQS4ZtvvsHQoUNha2uL/v37IysrC4cOHUJkZCQcHBxw1113oaSkRP+8Q4cOYfTo0fDw8ICzszNiYmKQmZmp/35qaipkMhn27Nmjv7Z48WJ4enqiqKjICL+9jkir1WqN9up/o6CgAL6+vti/fz+io6P1159//nns2rULBw8e/NvXaGxsRM+ePTFp0iS88cYb131fqVRCqVTqv66qqoKfnx8qKyvh5ORkmF+EiIgsTkNDAy5evIiAgIDrF67kbQIub/r7F3HoCvR+reW1428ANRf+/rmdxgN+42+t2P+xdetWjB8/Hvv370f37t0RHh6O+++/H8uWLWv18Tk5OQgICECPHj2QlJSEzp07Y/r06WhsbISjoyPefPNN2NnZ4eGHH0ZsbCxWrlwJANixYwcKCgoQGRkJrVaLpUuXYsuWLTh37hwcHR0B6D7Tv/nmGxw7dgzZ2dkYOHAgvv32W9x3333X1XGz97yqqgrOzs639Pkt6A7FHh4ekEgk16W3oqIieHt739Jr2NjYoG/fvjh//nyr35fL5ZDL5XdcK5m5rOVAYzVg4wh0myl0NURk6dR1gLLs7x8n97j+WmPlrT1XXdf2uq65++67MWPGDDz66KOIjIyEvb09EhMT//Z58+bNQ1xcHABg9uzZmDRpElJSUjB48GAAwOOPP45169bpHz9y5MgWz1+9ejVcXFywa9cu3HPPPQCAN998E9u2bcOTTz6JEydOYNq0aa0GG0MSdFpKJpMhIiKixTCZRqNBSkpKi5Gcm1Gr1Th+/Dg6duxorDKpPSg7BJTs090TERmbxA6Qu//9zcb5+ufaON/acyV3tjPykiVL0NTUhG+//RZffvnlLf2Pfp8+ffR/bm4Z6d27d4trf13NXFRUhBkzZiA4OBjOzs5wcnJCTU0NcnNz9Y+RyWT48ssv8f3336OhoQHvvffeHf1et0Lws6USEhIwbdo0REZGIioqCklJSaitrUV8fDwAYOrUqfD19dUnzv/85z8YOHAggoKCUFFRgXfffReXLl3CE088IeSvQURE1sRv/G1PGV03TWUkFy5cQEFBATQaDXJyclqElBv563lOzbsF/+81jUaj/3ratGkoKyvD+++/jy5dukAulyM6OhoqlarF6+7fvx8AUF5ejvLyctjb29/R7/Z3BA83EydORElJCRYsWIDCwkKEh4cjOTlZnxhzc3MhFv85wHT16lXMmDEDhYWFcHV1RUREBPbv34+QkBChfgUiIiKzolKpMHnyZEycOBHdu3fHE088gePHj8PT09OgP2ffvn1YsWIF7r77bgBAXl4eSktLWzzmwoULmDt3LtasWYONGzdi2rRp2L59e4vPdkMTPNwAwKxZszBr1qxWv5eamtri6/fee88kQ1pERETt1SuvvILKykp88MEHcHBwwNatWzF9+nRs2bLFoD8nODgY69evR2RkJKqqqjB//nzY2trqv69WqzF58mTExcUhPj4eY8eORe/evbF06VLMnz/foLX8leCb+BEREZHhpKamIikpCevXr4eTkxPEYjHWr1+PPXv26Fc5Gcqnn36Kq1evol+/fpgyZQqeffbZFqNDb731Fi5duoSPP/4YANCxY0esXr0ar776Ko4dO2bQWv5K0KXgQmjLUjIyH5X1jThXVI38inrkV9SjtFqFOlUTapRNqFOpMV3ybziJK1EncsG30jfhZm8DDwc53B3k8HZSwN/DDh2dbSER83RfIro1N10KTkZhEUvBiVqj1WpxqawOe8+XIu1CGf7Ir0Beef1Nn3N3pwY0SJWoaKrH9zfY0FEmEaOzux16eDsi1NcZvX2dEerjDGc7m1YfT0RE7RPDDZmNc0XV+PlYAX7+4woultZe930fZwU6u9vBx8UWno4KOMglsJNJYS+XIKLQFTZqoF7kghdCe6C8VomyGhVKapQoqKhHbnkdVGoNzhfX4HxxDbb8cUX/uj28HTEw0B0DA90QFeAON3uZKX9tIiIyMIYbElSTWoPkk4VYu/ciMnMr9NdtJCL06+yKIUEeiPB3RUhHJ7jY3SR0pNkDygZA7oDg6K7XfVut0aKgoh7ZpbU4WVCJk/lVOJ5fidzyOpwprMaZwmqs258DAAj1dcKoHl4YHeKFXj5O+uWQRETUPjDckCBqlU348uAlrNuXg4LKBgCAVCxCTLcOuC/cB6N6esFB3oa/np7DgKYaQOrQ6rclYhH83Ozg52aHmG4d9NdLqpVIv1iOgxfLcCC7DFlFNTiRX4UT+VV4P+UcvJ0UGNXTE/f08cGAADeI2bNDZHWsrDVVUIZ6r9lQTCbV0KjGhoO5WJF6HqU1uk2e3O1leHRgF0we2BmejsI27ZXWKLHzTDG2ny7CnnOlqFOp9d/zdlLgnj4dcX+4L0J9OaJDZOnUajWysrLg6ekJd3d3ocuxCpWVlSgoKEBQUFCLzQOBtn1+M9yQyew4U4QFm0/i8lVdc3AXdzvMHB6E+8J9oLCRCFzd9Roa1UjLLkPy8UJsPXEF1Q1N+u8FeNjjwX6+eCjSD15OXEVBZKmuXLmCiooKeHp6ws7Ojv9TY0QajQYFBQWwsbFB586dr3uvGW5uguHG9IqqGrBw80kknywEoBsBmR0bjAkRnWAjaR9bLSmb1Nh1tgSbjxVg+6kiKJt0249LxCKM6uGJSQM6Y1hwBy41J7IwWq0WhYWFqKioELoUqyAWixEQEACZ7PoeS4abm2C4Ma3fThbihe//QEVdIyRiER4fEoA5scGwk7Xfdq8aZROSTxTi6/RcHL50VX/d18UWE/v74f8GdIaHA0+iJ7IkarUajY2NQpdh8WQy2Q2PZWC4uQmGG9NoaFTjzV9O4YsDupNhQ32dsPjBMIT4GOk9T38aUJUDMjcgapVxfkYrzhVV46v0PHyfeRmV9br/8MkkYtwX7oP4wf7o5dPKicBERNRmDDc3wXBjfIWVDXhq/WEcu1wJAHhyWCDmjekOmdSIU1BpjwHKMkDuDkSvM97PuYGGRjV+PXEF/91/CUfzKvTXBwa6IX5wAGJ7enHKiojoDnCHYhLM0bwKPPn5YRRXK+FiZ4P3H+nbYum1pVLYSPCPvp3wj76dkJl7FZ/ty8HW41dwILscB7LLEehhj6djumJ8X1/jhjwiIuLIDRnOzrPFeHp9BpRNGgR7OuCTaZHo4m5vmh8u8MhNa65U1mN92iV8eTBXP2XV0VmBJ4YGYlKUX7vuOyIiMjVOS90Ew41xJJ8oxL++ykSjWosR3Tvgg0l94agw4ZlNZhhumtUom/DVwVys2ZON4molAMDFzgaPDfLHY4P8b77zMhERAWjb5zfHx+mO/XSsADM36ILNuN4dsXpqpGmDjZlzkEsxY1gg9rwwAokP9Ia/ux0q6hqRtP0chizaiWW/n0VlHVdhEBEZCsMN3ZHUs8WYu/Eo1BotHujni/cfCW83e9eYmlwqwaSozkh5bjg+nNQXPbwdUaNswgc7zmPI4h14f/s5VDUw5BAR3Sl+CtFt++NyBZ75MhNqjRbjw32wZEIYpAw2f0siFuHeMB9sfXYoVk3uh+5ejqhuaMJ727MwdNFOfLTjHGqUTX//QkRE1Cp+EtFtyS2rw/R1h1CnUmNosAcWTwjjoZJtJBaLMDa0I36dPRQf/V9fBHk6oLK+EUt+z8LQRTuwMvUC6v9ythUREd0ahhtqszpVE55cfxilNSr08nHCyskRXN58B8RiEe7p44Pf5gzD+4+EI9DDHlfrGrEo+QxGLEnFxkO5aFJrhC6TiKjd4GopahOtVou5G49i09ECeDjI8cuzQ8zj4MiyQ4BaCUjkgHt/oau5I01qDTYfLcB727P0h4wGeTrg+bjuGB3ixYP7iMgqcSn4TTDc3JnP03KwYPNJSMQibHhiAAYEugtdksVSNqnxxYFcfLTjHK5eW00V0cUVL93VA5H+bgJXR0RkWgw3N8Fwc/uO5VVgwqr9aFRr8eq4nnhiaKDQJVmFqoZGrN6VjU/2ZqOhUTc9FdvTC8+P7Y5uXo4CV0dEZBoMNzfBcHN76lVqjPtgD7JLa3FXqDdWPNqP0yMmVlTVgKTt5/DN4TyoNVqIRcAjUZ2RMLobTyEnIovHTfzI4BJ/PY3s0lp4OcmR+EBv8ws21eeByjO6ewvl5aRA4gO98ducYYjr5QWNFthwMBfD303FytQLaGjkyioiIoAjN0KX0y6kni3GY58dAgCsfzwKQ4PN8CBMMz5+wVgOZpfhzV9O43i+7vT1Tq62ePGuHhjXu6P5hU8iojvEkRsymKqGRrzw/R8AgMcG+ZtnsLFSAwLdsXnmYCx9KAxeTnJcvlqPWRuOYMKqNBzNqxC6PCIiwTDc0E0t+z0LRVVK+Lvb4YWxPYQuh/6HWCzCgxGdsHPecMyJDYatjQQZl65i/PJ9mPP1ERRU1AtdIhGRyTHc0A2dyK/E52k5AIA3x/eGrUwibEF0Q3YyKebEdsPOecPxYL9OAIBNRwswYkkqlv5+FrU8zoGIrAjDDbVKrdHilR+PQ6MF7g3zwZBgD6FLolvg7azA0ofD8POsIYgKcIOySYMPd5zH8CWp+OZwHjQaq2qxIyIrxXBDrfoqPRfHLlfCUS7Fa+N6Cl0OtVHvTs7Y+ORArJrcD53d7FBSrcTz3/2Bf6zYhyO5V4Uuj4jIqBhu6DpVDY1Yti0LAPDcmG7wNIfjFajNRCLdwZzbEobh5bt7wEEuxbHLlfjHiv2Y9+0xlFQrhS6RiMgoGG7oOqtSL6C8VoWuHewxeWAXocuhOySXSvDksK7YMS8GEyJ0/TjfZVzGyCWp+GRPNhp5KCcRWRiGG2rhSmU9Pt17EQDw4l09IZXwr4il8HRUYMlDYfjhmUHo08kZ1comvPnLaYxN2o0950qELo+IyGD4yUUtLPs9C8omDaL83RDb01PocsgI+nV2xaZnBmPxg33gbi/DhZJaTPk0HU9+fhi5ZXVCl0dEdMe4QzHpnSmswl3v74FWC/z4zCD07ewqdEm3rqkegBaACJDaCl1Nu1FZ34j3t5/Df9NyoNZoIZOK8dSwQPxzeFfYyaRCl0dEpMcdium2fJByDlotcHdv7/YVbABdoJHaMdi0kbOtDRbcG4Lk2UMxOMgdqmtLx2OX7sKWPwpgZf/vQ0QWguGGAABnC6ux9XghAGBObDeBqyFTC/ZyxBePD8Cqyf3g62KLgsoGzNpwBI+sPoAzhVVCl0dE1CYMNwQA+Gin7jTtu3t7o5uXo8DVkBCal46nPBeDubHdIJeKcfBiOcZ9sBdvbDmF6oZGoUskIrolDDeE88U12PJHAQBg1ohggau5TXmbgJwNunu6IwobCWbHBiPluRiM7eUNtUaLT/dexKilu7D5aD6nqojI7DHcEJbvPA+tFhgT4oUQn3baZH15E5Dzle6eDKKTqx1WTYnAuvj+8He3Q3G1ErO/PopJaw4gq6ha6PKIiG6I4cbK5ZXXYfPRfADAs6Pa6agNGdXw7p74be4wzBvTDQobMQ5kl+Pu9/fg7a2nUcMDOYnIDDHcWLl1+3Og0QJDgz0Q6ussdDlkpuRSCWaNDMa2uTEYE+KFJo0Wq3dnY9TSVPx8jKuqiMi8MNxYseqGRmw8lAcAeHxIgMDVUHvg52aH1VMj8dlj/dHF3Q5FVUr866sjePSTgzhfzKkqIjIPDDdW7LuMy6hRNqFrB3sMC+4gdDnUjozo4Ynf5gxDwmjdqqr9F8owNmkPEn89jVpOVRGRwBhurJRao8W6/TkAgPjBARCLRcIWRO2OwkaCZ0fppqpie3qiSaPFx7uyEbtsF7Yev8KpKiISDMONldpxphiXyurgbGuDB/r5Cl0OtWOd3e3wybT++GRqJPzcbHGlsgHPfJmJqWvTcaGkRujyiMgKMdxYqc/26U7+fiTKj2cIkUHEhnhh29wYzB4VDJlUjD3nSjE2aTcWJ59BnYpTVURkOgw3VuhiaS32XyiDSARMGdhF6HLIgihsJJg7uhu2zR2GEd07oFGtxYrUC4hdugvJJzhVRUSmwXBjhb4+lAsAiOnWAZ1c7QSuxkAcugJO3XX3JLgu7vZY+1h/rJ4SoT+r6ukvMjHts0O4WFordHlEZOFEWiv7X6m2HJluiRrVGkQnpqC0RoVVkyMwNtRb6JLIwtWr1FiReh4f78qGSq2BTCLGk8MCMXNEEGxlEqHLI6J2oi2f3xy5sTIpp4tQWqOCh4Mco3p6Cl0OWQFbmQTPjemO3+YOw7BuHaBSa/DRzvOIXbYLv50s5FQVERkcw42V+Spdt2nfQ5GdYCPhP34ynQAPe/w3vj9WTe4HXxdb5FfU46n1GZi+7hByOFVFRAbETzcrcvlqHXafKwEAPNLfT+BqyBqJRCKMDe2IbQnDMHNEV9hIRNh5tgRj3tuNZb+fRb1KLXSJRGQBGG6syHcZl6HVAoO6uqOLu73Q5RjW8TeAzHm6ezJ7djIp5sf1wG9zhmFosAdUag0+2HEeo9/bhW2nijhVRUR3hOHGSmi1Wmw6ojv9++FICxy1qbkAVJ3V3VO7EdjBAZ9Pj8LKR/vBx1mBy1frMePzw3j8v4dxqYxTVUR0exhurMSxy5XIKauDrY0Eo0O8hC6HSE8kEuGu3h2x/bkY/HO4bqpqx5lijH5vN5Zty0JDI6eqiKhtzCLcLF++HP7+/lAoFBgwYADS09Nv6Xlff/01RCIRxo8fb9wCLcDmo7pRm9EhXrCXc0diMj92MileGNsDyXOGYUiQB1RNGnyQcg6j39uF7aeKhC6PiNoRwcPNxo0bkZCQgIULFyIzMxNhYWGIi4tDcXHxTZ+Xk5ODefPmYejQoSaqtP1qUmvw87ErAIDxfX0Erobo5rp2cMD6x6Ow4tF+6OisQF55PZ74/DAeX3cIuWV1QpdHRO2A4OFm2bJlmDFjBuLj4xESEoJVq1bBzs4Oa9euveFz1Go1Hn30Ufz73/9GYGDgTV9fqVSiqqqqxc3apGWXobRGCVc7GwwN7iB0OUR/SyQS4e7eHbE9IQZPx3SFVCxCyplixL63C0nbOVVFRDcnaLhRqVTIyMhAbGys/ppYLEZsbCzS0tJu+Lz//Oc/8PT0xOOPP/63PyMxMRHOzs76m5+fBTbT/o1NRwoAAOP6dOTeNtSu2MulePGuHkieMxSDg9yhatIgafs5jHlvN3ac4VQVEbVO0E+60tJSqNVqeHm1bHD18vJCYWFhq8/Zu3cvPv30U6xZs+aWfsZLL72EyspK/S0vL++O625PGhrV+O2k7r0cH+4rcDVEtyfI0xFfPD4AH/1fX3g7KZBbXofp6w7jif8eRl45p6qIqKV21VlaXV2NKVOmYM2aNfDw8Lil58jlcsjlciNXZr5SzxajRtkEXxdb9OvsKnQ5RLdNJBLhnj4+GNHdEx+knMOney9i++ki7DlXgmeGB+GpmEAobHhWFREJHG48PDwgkUhQVNRyeLmoqAje3tcf6HjhwgXk5OTg3nvv1V/TaDQAAKlUirNnz6JrV54K/VfJJ3SjNnf39oZYLBK4GqI7Zy+X4qW7e2JCRCcs2HwSadlleG97Fn44chmv39sLI3rwzDQiayfotJRMJkNERARSUlL01zQaDVJSUhAdHX3d43v06IHjx4/j6NGj+tt9992HESNG4OjRo1bZT3MzqiYNUs7oVp3F9bLw0787jQf8J+nuySoEezliw4wB+HBSX3g5yXGprA7x6w5hxuecqiKydoJPSyUkJGDatGmIjIxEVFQUkpKSUFtbi/j4eADA1KlT4evri8TERCgUCoSGhrZ4vouLCwBcd510q6SqG5rQwVFu+VNSfuOFroAEIBKJcG+YD0b00E1Vrd17EdtOFWF3VglmjQjCjGGcqiKyRoKHm4kTJ6KkpAQLFixAYWEhwsPDkZycrG8yzs3NhVjMFT63o7mReHSIF6ekyKI5yKV4WT9VdQIHssuxdFsWvs+8jNfv64Xh3TlVRWRNRForO6GuqqoKzs7OqKyshJOTk9DlGI1ao8WAt1NQWqPEf6dHIaYb97ch66DVavHTsQK89ctpFFcrAQBxvbzw2j0h6ORqJ3B1RHS72vL5zSERC3Uk9ypKa5RwVEgRHegudDnG11QPNNXp7smqiUQi3B/ui5TnYjBjaAAkYhF+O1mE2GW78P72c9wAkMgKMNxYqOZVUqN6eEImtYJ/zIf+CeydqLsnAuCosMEr40Kw9dmhGBDghoZGDd7bnoVRS3ch+cQVWNmgNZFVsYJPPeuj1Wrx2ylduLH4VVJEf6O7tyO+fnIgPpzUFx2dFcivqMfTX2RiyqfpOFdULXR5RGQEDDcW6EJJDfLK6yGTihHTnb02RM2rqlKei8G/RgZBJhVj7/lSjH1/D/7z8ylU1jcKXSIRGRDDjQXaeaYEADAgwA12MsEXxBGZDTuZFM+N6Y7tc2MwJsQLao0Wa/ddxMglqfg6PRdqDaeqiCwBw40F2nlWt3HfCC5/JWpVZ3c7rJ4aic+nR6FrB3uU1arw4g/HMX75PmRcuip0eUR0hxhuLEyNsgmHcsoBgNvQE/2NYd06IHnOMLw6ricc5VIcz6/Egyv3I+GboyiuahC6PCK6TQw3Fmbf+VI0qrXo4m6HAA97ocshMns2EjGeGBqIHfOG4+HITgCAHzLzMWJJKj7edQGqJo3AFRJRWzHcWJhUTkkR3ZYOjnIsnhCGTTMHI9zPBbUqNRJ/PYOxSbv1U71E1D4w3FgQrVaL1LO6ZuLhXCVFdFvC/Vzwwz8HYclDYfBwkCO7tBbxnx3C4+sOIae0VujyiOgWMNxYkLNF1bhS2QC5VIyB1rArMZGRiMUiTIjohJ3zYvDksEBIxSKknCnGmPd2I3HraVQ3cOk4kTnjOmEL0rwEfFBXd+s7CTn0VUDTBIj5V5oMx1Fhg5fv7omHI/3wny2nsDurBB/vzsb3mZcxb0x3PBTpBwkPpSUyOxy5sSC7s5qnpKyw38YxCHDuobsnMrAgTwf8N74/1j4WiUAPe5TW6JaO3/vhXqRdKBO6PCL6Hww3FqKhUY2MXN3+HIODPASuhsjyiEQijOzhheQ5w/DaPSFwUkhx6koVJq05gKfXZyC3rE7oEonoGoYbC5F56SpUTRp4OcnRtQOXgBMZi0wqxuNDApA6fwSmDOwCsQhIPlmI2GW7kPgr+3GIzAHDjYXYd6EUADCoqwdEIivsASg7BBTv1d0TmYCbvQxvjA/Fr7OHYWiwB1RqDT7elY0RS3Zh4yEe5UAkJIYbC7H/2rx/dFcrXSWVtRw4tUh3T2RC3b0d8fn0KHw6LRIBHvYorVHihe91/TgHstmPQyQEhhsLUN3QiD8uVwLQrZQiItMSiUQY1dMLvzUf5XCtH+eR1Qfwzy/Yj0Nkagw3FuBQTjnUGt2RC51c7YQuh8hqyaS6oxxS5w3H5IGdIRYBv57Q9eMsSj7DfhwiE2G4sQD7zuuGvjlqQ2Qe3B3keHN8b/w6exiGBOn6cVamXmA/DpGJMNxYgOZ+m0FduQScyJx093bE+sej8MnUlv0493y4F3vPlQpdHpHFYrhp58prVTh9pQoAeOQCkRkSiUSIDdH147xyt64f5/SVKkz+9CDiP0tHVlG10CUSWRyGm3aueTVGdy9HdHCUC1wNEd2ITCrGjGGB2DV/BB4b5A+pWISdZ0swNmk3Xv7xOEqqlUKXSGQxGG7aufSL5QCAgYFuAldCRLfCzV6G1+/rhd/nDkNcLy9otMCGg7kY/u5OfLTjHOpVaqFLJGr3GG7aucOXdOEm0p/hhqg9CezggI+nRGLjkwPRp5MzalVqLPk9CyOXpuKHzMvQsOmY6LYx3LRjNcomnCrQ9dtE+rsKXI3AJApAaqu7J2pHBgS6Y9Mzg/H+I+HwdbHFlcoGJHxzDPct56GcRLdLKnQBdPuO5F6FRgt0crVFR2dbocsRVtQqoSsgum1isQj3h/sirpc3PtuXgxU7z+NEvu5Qztiennjxrp4I8nQQukyidoMjN+3Y4RzdKeCRXax81IbIQihsJPjn8K5InT8cU6O7QCIWYfvpYsQl7cZrm06grIZNx0S3guGmHWO/DZFlcneQ4z/3h+K3OcMQ29MLao0W6w9cwvB3U7Ei9TwaGtl0THQzDDftVKNagyO5FQCA/gw3RBYpyNMBn0yLxIYZAxDq64RqZRMWJ5/F8HdTudMx0U0w3LRTp69UoU6lhpNCimDOxQMX1gJnP9DdE1mYQV098NPMIXhvYhh8XWxRWNWAF74/jrFJu7H9VBG0WoYcor9iuGmnDjX32/i7QSwWCVyNGSjeDVzZprsnskBisQj/6NsJKc/F4NVxPeFiZ4NzxTV44vPDmPjxAWTmXhW6RCKzwXDTTmVc67eJYDMxkVVR2EjwxFDdTsdPx3SFXCpGek45HlixH0+vz8CFkhqhSyQSHMNNO6TVavUjN+y3IbJOzrY2ePGuHkidPxwPR3aCWAQknyzEmPd245Ufj6O4ukHoEokEw3DTDuWV16OkWgkbiQh9OjkLXQ4RCaijsy0WTwjDr7OHIbanJ9QaLb48mIvh76Zi2bYs1CibhC6RyOQYbtqho5crAAAhHZ2gsJEIWwwRmYXu3o74ZFp/bHxyIML9XFCnUuODlHOIWbwT/92fA1WTRugSiUyG4aYd+iOvAgAQ5uciaB1EZH4GBLrjx2cGYeWj/RDgYY+yWhUW/nQSo9/bhZ+PFXBlFVkFhpt26Ni1kZuwTi6C1kFE5kkkEuGu3h3x+9xheGN8KDwc5LhUVod/fXUE9y/fh73nSoUukcioGG7amSa1BsfzKwEAYX7styGiG7ORiDFlYBfsmj8cc2O7wV4mwR+XKzH504P4vzUHcPTaKDCRpWG4aWeyimrQ0KiBg1yKQA9u3kdEf89eLsXs2GDsen4E4gf7QyYRY/+FMoxfvg9Pfn4YWUXVQpdIZFAMN+3MH9empPp0cubmfX/l3h/oMFh3T0St8nCQY+G9vbBjXgwmROiWj/9+qghjk3bjuW+OIa+8TugSiQxCKnQB1Db6fhs2E7fUbabQFRC1G51c7bDkoTA8NSwQS3/PQvLJQnyfeRk/HcvHowO6YOaIIHRwlAtdJtFt48hNO3M071q/DZuJiegOBXs5YtWUCGyaORiDg9zRqNZi3f4cxLy7E0t/P4uqhkahSyS6LQw37Ui9Sq2fG2czMREZSrifC758YiC+fGIAwjo5o06lxoc7zmPY4p34eNcFNDSqhS6RqE0YbtqRkwWVUGu08HSUw9tJIXQ5RGRhBgd5YNPMwVg1OQJBng6oqGtE4q9nEPPuTnx58BIa1dwIkNoHhpt25OhfNu8TidhM3ELGXCDtMd09Ed02kUiEsaHe+G3OMCx5KAy+LrYoqlLilR9PYNTSXfg+4zLUGm4ESOaN4aYdOXZZ128Tzmbi66muAsoy3T0R3TGJWIQJEZ2wY14MFt4bAg8HGXLL6/Dct8cw+r1d+OlYATQMOWSmGG7akePXVkr19mW/DRGZhlwqQfzgAOx+fgReGNsDLnY2yC6pxbNfHcFd7+9B8olCHulAZofhpp2oUTYhp0y3B0Uoww0RmZidTIp/Du+KPc+PQMLobnBUSHG2qBpPf5GBez/aix1nihhyyGww3LQTp69UAQA6OivgZi8TuBoislaOChs8OyoYe58fiX+NDIK9TIIT+VWYvu4w/rFiP/acK2HIIcEx3LQTpwp04Sako5PAlRARAc52NnhuTHfseWEknooJhMJGjKN5FZjyaTomfnwAB7PLhC6RrBjDTTtxskDXTBziw3BDRObDzV6Gl+7qid3N51ZJxUjPKcfE1Qcw+ZODyMxlkz+ZHsNNO3Hq2rRUL4YbIjJDno4KLLy3F3bNH47JAzvDRiLC3vOleGDFfsR/lo7j11Z7EpkCw0070KjWIKuwBgAQ0pHNxERkvjo62+LN8b2x47nhmBjpB4lYhJ1nS3DvR3vx5OeH9f2DRMbEcNMOnC+ugUqtgaNcik6utkKXQ0T0t/zc7LBoQh+kJMTgH319Ibp2Avld7+/BU+sP66faiYyBp4K3A83NxD19nCAWc2fiVgXGAxolIOZJxkTmxN/DHu9NDMfMEV2RtP0cfjl+Bb+dLMJvJ4swJsQLz44K5vYWZHBmMXKzfPly+Pv7Q6FQYMCAAUhPT7/hY3/44QdERkbCxcUF9vb2CA8Px/r1601Yrek199twpdRNeMUAHcfo7onI7AR5OuKj/+uH3+cMw31hPvqRnHs+3Isn/nsYJ/I5kkOGI3i42bhxIxISErBw4UJkZmYiLCwMcXFxKC4ubvXxbm5ueOWVV5CWloY//vgD8fHxiI+Px2+//Wbiyk2HK6WIyFIEeznig0l9sW3uMNwf7gOxCNh+ujnkHGLjMRmESCvwbksDBgxA//798dFHHwEANBoN/Pz88K9//QsvvvjiLb1Gv379MG7cOLzxxhvXfU+pVEKpVOq/rqqqgp+fHyorK+HkZP5hQavVIuzfv6OqoQm/PDsEvXw4fEtEluNCSQ0+2nEem4/mo/moqpE9PDF7VDDCeI4e/UVVVRWcnZ1v6fNb0JEblUqFjIwMxMbG6q+JxWLExsYiLS3tb5+v1WqRkpKCs2fPYtiwYa0+JjExEc7Ozvqbn5+fweo3hfyKelQ1NMFGIkKwp6PQ5ZivunygNld3T0TtRtcODnhvYji2J8Tggb6+EIuAHWeKcf/yfYj/LB1H8yqELpHaIUHDTWlpKdRqNby8vFpc9/LyQmFh4Q2fV1lZCQcHB8hkMowbNw4ffvghRo8e3epjX3rpJVRWVupveXl5Bv0djO3ktWbiYE9HyKSCzyKar2OvAIdm6u6JqN0J7OCAZRPDkfLccDzYr5N+Cfn45fswbW06NwOkNmmXq6UcHR1x9OhR1NTUICUlBQkJCQgMDMTw4cOve6xcLodc3n5X0OhXSrGZmIisQICHPZY+HIZ/jQzCRzvP48cj+diVVYJdWSUY1q0DZo8KRkQXV6HLJDMnaLjx8PCARCJBUVFRi+tFRUXw9va+4fPEYjGCgoIAAOHh4Th9+jQSExNbDTftXVZRNQCgZ0dOSRGR9fD3sMeSh3QhZ/nO8/g+Mx+7s0qwO6sEQ4M9MHtUMCL93YQuk8yUoPMcMpkMERERSElJ0V/TaDRISUlBdHT0Lb+ORqNp0TRsSc5eCzfdvBhuiMj6dHG3x+IJYdh5bcdjqViEPedKMWFVGh5ZnYZ950t5CjldR/BpqYSEBEybNg2RkZGIiopCUlISamtrER8fDwCYOnUqfH19kZiYCEDXIBwZGYmuXbtCqVRi69atWL9+PVauXCnkr2EUDY1qXCqrA8BwQ0TWrbO7bsfjWSODsCL1PL7LuIwD2eU4kH0QfTu7YNaIIIzs4QmRiBudkhmEm4kTJ6KkpAQLFixAYWEhwsPDkZycrG8yzs3NhVj85wBTbW0tnnnmGVy+fBm2trbo0aMHvvjiC0ycOFGoX8FosktqodZo4aSQwsup/fYNEREZip+bHRIf6IN/jQzG6t3Z+Co9F0dyK/D4fw8jpKMTZo0Mwthe3tzN3coJvs+NqbVlnbzQNh/Nx+yvjyKyiyu+++cgocsxb2mPAcoyQO4ORK8TuhoiMpGSaiU+2ZuN9WmXUKdSAwCCPB0wc0RX3NvHB1IJV5lainazzw3dXHMzcTdvTkkREbWmg6McL93VE/teGIlnRwXDUSHF+eIazN14DCOX7sLX6blQNWmELpNMjOHGjJ0trAEAdPN0ELgSIiLz5movQ8Lobtj34kjMj+sON3sZcsvr8OIPxzH83Z347/4cNDSqhS6TTIThxoydK+bIDRFRWzgpbDBzRBD2vjACr47rCU9HOQoqG7Dwp5MYsmgnVu++gFplk9BlkpGx58ZM1avUCFmYDK0WOPxqLDwc2FB8U8pyQKsBRGJAzr0viEinoVGNbzMuY1XqBeRX1AMAXOxs8PjgAEwd5A9nWxuBK6RbxZ4bC3C+uAZaLeBmL2OwuRVyN0DhwWBDRC0obCSYMrALUucPx+IJfeDvboeKukYs3ZaFIe/swJLfzqK8ViV0mWRgDDdm6s/N+9hvQ0R0p2wkYjwc6YeU54bj/UfC0c3LAdXKJny08zwGvZOC1386iYJrIzvU/jHcmKlz3JmYiMjgJGIR7g/3RfLsYfh4SgR6+zqjoVGDdftzMGzxTsz79hjOF9cIXSbdIcE38aPW8diFNipIBtQNgEQB+IwVuhoiMnNisQhxvbwxJsQLe8+XYsXOC0jLLsN3GZfxfeZlxIV445kRXdGnk4vQpdJtMHi4qaurg52dnaFf1uqcK7q2DJzh5tZc+vrPTfwYbojoFolEIgwN7oChwR1wJPcqVqRewLZTRUg+WYjkk4UYEuSBZ4Z3RXRXdx7t0I7c1rTUqFGjkJ+ff9319PR0hIeH32lNVq+6oVHf1c+eGyIi0+jb2RVrpkbi97nD8EA/X0jEIuw9X4r/++Qgxq/Yj99OFkKjsaoFxu3WbYUbhUKBPn36YOPGjQB0p3K//vrrGDJkCO6++26DFmiNzl2b7/V0lMPFTiZwNURE1qWblyOWPRyO1HnDMTW6C+RSMY7lVeCp9RmIS9qN7zMuo1HNXY/N2W1NS/3yyy9Yvnw5pk+fjs2bNyMnJweXLl3Cli1bMGbMGEPXaHUuXAs3QdyZmIhIMH5udvjP/aH418hgfLbvItanXcK54ho89+0xLNuWhRlDAzCxf2fYyiRCl0r/47Z7bmbOnInLly9j0aJFkEqlSE1NxaBBPNzREC6W1gIAAjvYC1wJERF1cJTj+bE98PTwrvjyQC4+3XsR+RX1eP3nU/hwx3nED/bHlGhuCGhObmta6urVq3jwwQexcuVKfPzxx3j44YcxZswYrFixwtD1WaXskmvhxoMjN0RE5sJJYYN/Du+KvS+MwBvjQ9HJ1RZltSos+T0Lg9/ZgcRfT6O4ukHoMgm3GW5CQ0NRVFSEI0eOYMaMGfjiiy/w6aef4rXXXsO4ceMMXaPVyS7VTUtx5IaIyPzodz2eNxxJE8PR3csRNcomfLwrG0MW7cTLPx7Xj8CTMG4r3Dz99NPYvXs3AgIC9NcmTpyIY8eOQaXiNtZ3Qq3RIqesDgBHboiIzJlUIsb4vr74dfZQfDI1Ev06u0DVpMGGg7kYuTQVz3yZgWN5FUKXaZV4cKaZyS2rw7B3d0ImFeP0f8ZCIua+Crck7bE/97mJXid0NURkhbRaLdIvluPj3dnYcaZYfz060B1PxQQiplsH7pVzB9ry+c0dis3MhWtTUv7udgw2bWHnC0jtAJmr0JUQkZUSiUQYEOiOAYHuOFNYhdW7s/HT0QKkZZchLbsMPbwd8XRMV4zr0xE2Ep5+ZEwcuTEzn+69iDe2nMLYXt5YNSVC6HKIiOgOFFTU49O9F/FVei7qVGoAgK+LLZ4YGoCJ/f1gJ+MYw61qy+c3o6OZyS5hMzERkaXwcbHFa/eEIO3FUZgf1x0eDjLkV9Tj3z+fwqB3dmDZ72dRVqMUukyLw3BjZvTLwDuwmZiIyFI429lg5ogg7H1hJN76Ryj83e1QUdeID3acx6B3duC1TSeQe20xCd05hhszw2XgRESWS2EjwaMDuiDlueFY8Wg/9OnkDGWTBusPXMLwJTsxa0MmTuRXCl1mu9fmyb7Tp0/j66+/xp49e3Dp0iXU1dWhQ4cO6Nu3L+Li4vDggw9CLpcbo1aLV6NsQlGVbniyK5eBt83pJUBjFWDjBPScJ3Q1REQ3JRGLcHfvjrgr1Btp2WX4eFc2dmWVYMsfV7DljysYEuSBp2ICMSTIgyusbsMtNxRnZmbi+eefx969ezF48GBERUXBx8cHtra2KC8vx4kTJ7Bnzx5UVVXh+eefx5w5c8wy5JhzQ/Hxy5W496O9cLeXIeO10UKX075wKTgRtXOnCqqwevcF/PzHFaivnT7ey8cJT8V0xd2h3pBa+QoroywFf/DBBzF//nx89913cHFxueHj0tLS8P7772Pp0qV4+eWXb7lo4pQUEZE1C/FxQtIjfTEvrjs+2XMRGw/l4WRBFZ796ggWu9pixtBAPBzpx4M6b8Eth5usrCzY2Pz9oWDR0dGIjo5GY2PjHRVmjS7wTCkiIqvXydUOr9/XC7NHBWP9gUtYtz8Hl6/WY+FPJ5G0PQvTBvljysAucHcwv9kRc3HLY1y3EmwAoK6urk2Ppz9xGTgRETVztZfh2VHB2PfCSLxxfy/4udnial0jkrafw6B3duDVTTzD6kZuawJv1KhRyM/Pv+56eno6wsPD77Qmq8Vl4ERE9L9sZRJMifbHzueG48NJfdHbV7fC6osDujOsnlp/GBmXrgpdplm5rXCjUCjQp08fbNy4EQCg0Wjw+uuvY8iQIbj77rsNWqC10Gq1yCnThZsAD47cEBFRS1KJGPeG+eCnWYPx1YyBGNnDE1ot8NvJIjy4cj8eXLkfyScK9c3I1uy29n3+5ZdfsHz5ckyfPh2bN29GTk4OLl26hC1btmDMmDGGrtEqlNQoUadSQyQCOrvZCV0OERGZKZFIhOiu7oju6o5zRdVYsycbm44UIOPSVWRcykCAhz0eHxKACRGdoLCxzubj2z7UYubMmbh8+TIWLVoEqVSK1NRUDBo0yJC1WZXmnSl9nG0hk1r3cj8iIro1wV6OWDwhDPPGdMd/03LwxYFcXCytxaubTmDZtixMGdgFU6Otr/n4tj5Fr169igcffBArV67Exx9/jIcffhhjxozBihUrDF2f1bh0Ldxw1IaIiNrK00mB+XE9sP/FkVh4bwg6udqivFaF91N0zcev/Ghdzce3NXITGhqKgIAAHDlyBAEBAZgxYwY2btyIZ555Br/88gt++eUXQ9dp8S6V68JNF3eGm9vSMQ5oqgWk7FciIutlL5cifnAApgzsguSThVi9Oxt/XK7ElwdzsSE9F2NCvPDksEBEdHETulSjuq2Rm6effhq7d+9GQECA/trEiRNx7NgxqFQqgxVnTfKuhZvODDe3x38SEPSE7p6IyMpJJWLc08cHm2cOxtdPDsSoFs3HaXhgxT4kn7hisc3Ht3z8gqUw1+MXHlixD5m5FVj+f/0wrk9HocshIiILc764Gmt2X8SPR/KhUmsAAP7udnh8aCAm9Otk9jsft+Xz+5ZHbnJzc9tURGv74NCN5XJaioiIjCjI0xGLJvTB3hdHYNaIIDjb2iCnrA6vbTqBwYt24L1tWSirUQpdpkHccrjp378/nnrqKRw6dOiGj6msrMSaNWsQGhqK77//3iAFWoMaZRNKa3TTeZyWIiIiY/J0VGBeXHekvTQSr7fSfPzyj8f1O+a3V7fcUHzq1Cm89dZbGD16NBQKBSIiIuDj4wOFQoGrV6/i1KlTOHnyJPr164fFixdzM782aF4G7mpnAycFj624LTwVnIioTexkUjw2OACTB3bBbyeLsHr3BRy7XIkNB3PxVXouRvdsbj52hUgkErrcNmlzz019fT1++eUX7N27F5cuXUJ9fT08PDzQt29fxMXFITQ01Fi1GoQ59twkn7iCp7/IRJifCzbPHCx0Oe0Tww0R0R3RarVIv1iONXuysf10sf56384ueHJoIMb08oZELFzIacvnd5uXgtva2mLChAmYMGHCbRdILXGPGyIiEppIJMKAQHcMCHTH+eJqfLLnIn7IzMeR3Ar888tMdHazw+NDAvBQZCfYyW57D2CTaNNS8OzsbFjZ4iqT0O9xw3BDRERmIMjTEe88qGs+/tfIILjY2SC3vA4LfzqJ6MQdePe3MyiuahC6zBtqU7gJDg5GSUmJ/uuJEyeiqKjI4EVZm+aeGzYTExGROfF0VOC5Md2x/8WReOP+XvB3t0NlfSOW77yAIYt2Yv63x3C2sFroMq/TpnDzv6M2W7duRW2t9WznbCy5HLkhIiIzZieTYkq0P1KeG45VkyMQ2cUVKrUG32ZcRlzSbkxbm46950rNZnbHvCfNrECjWoP8inoAQBd3Hh1ARETmSyIWYWyoN8aGeiMz9yo+2ZON5BOF2JVVgl1ZJejZ0QlPDgvAPX18YCMR7hDoNv1kkUh03XKw9rY8zNwUVNRDrdFCLhXD09G6Tm0lIqL2q19nV6x4NAKp80bgsUH+sLWR4PSVKszdeAyjl+1C47VdkIXQppEbrVaLxx57DHK57kO4oaEBTz/9NOztW444/PDDD4ar0ML9daWUWMAldkRERLejs7sdXr+vF+bEBuPLg7lYtz8H0V09BB25aVO4mTZtWouvJ0+ebNBirBFPAyciIkvgYifDzBFBeGJoAOpVakFraVO4+eyzz4xVh9XKLdM1ZHd2Y7/NHen5HKBpBMTc4ZmISEhyqQRyqbCHcLKhWGDNK6U6u9kKXEk759Jb6AqIiMhMCDchRgCgXynVyZXTUkRERIbAcCOw/Ku6cOPrypEbIiIiQ+C0lIDqVE24WtcIgOHmjlUc/7PnhlNURERWjeFGQM2jNo4KKZwUbIS9I6eX8lRwIiICwGkpQV2+1m/j68JRGyIiIkNhuBHQ5avNzcQMN0RERIbCcCMgfTMxR26IiIgMxizCzfLly+Hv7w+FQoEBAwYgPT39ho9ds2YNhg4dCldXV7i6uiI2NvamjzdnzcvA2UxMRERkOIKHm40bNyIhIQELFy5EZmYmwsLCEBcXh+Li4lYfn5qaikmTJmHnzp1IS0uDn58fxowZg/z8fBNXfufyr+o28OMeN0RERIYjeLhZtmwZZsyYgfj4eISEhGDVqlWws7PD2rVrW338l19+iWeeeQbh4eHo0aMHPvnkE2g0GqSkpJi48juXz4ZiIiIigxM03KhUKmRkZCA2NlZ/TSwWIzY2Fmlpabf0GnV1dWhsbISbm1ur31cqlaiqqmpxMweqJg2Kq5UAOC1FRERkSIKGm9LSUqjVanh5ebW47uXlhcLCwlt6jRdeeAE+Pj4tAtJfJSYmwtnZWX/z8/O747oN4UplPbRaQGEjhru9TOhyiIiILIbg01J34p133sHXX3+NH3/8EQqFotXHvPTSS6isrNTf8vLyTFxl65pXSvm42EIkEglcDRERkeUQdIdiDw8PSCQSFBUVtbheVFQEb2/vmz53yZIleOedd7B9+3b06dPnho+Ty+WQy+UGqdeQuIGfgXFXYiIiukbQkRuZTIaIiIgWzcDNzcHR0dE3fN7ixYvxxhtvIDk5GZGRkaYo1eDyuYEfERGRUQh+tlRCQgKmTZuGyMhIREVFISkpCbW1tYiPjwcATJ06Fb6+vkhMTAQALFq0CAsWLMCGDRvg7++v781xcHCAg4ODYL9HW3GlFBERkXEIHm4mTpyIkpISLFiwAIWFhQgPD0dycrK+yTg3Nxdi8Z8DTCtXroRKpcKECRNavM7ChQvx+uuvm7L0O6LfnZgjN0RERAYl0mq1WqGLMKWqqio4OzujsrISTk5OgtUxdPEO5JXX45unohEV0PoydmqDnK+AplpAag/4TxK6GiIiMrC2fH4LPnJjjdQaLa5UNADgyI3BXPkNUJYBcneGGyIiK9eul4K3V8XVDWjSaCERi+DlaH4ruYiIiNozhhsBNPfbeDspIJXwHwEREZEh8ZNVADwNnIiIyHgYbgRwpfJavw2XgRMRERkcw40ACq+FG2/n1o+MICIiotvHcCMAfbhxYrghIiIyNIYbAVyp4sgNERGRsTDcCKDo2shNR4YbIiIig+MmfibWpNaguJrTUgbnEgo0VgE2wu06TURE5oHhxsRKapTQaAGpWAR3B27gZzA95wldARERmQlOS5lY8zJwLycFJGKRwNUQERFZHoYbEyvShxuO2hARERkDw42JXdE3E3MDPyIiImNgz42JFXIZuHEcewVQXQVkrkDYW0JXQ0REAmK4MTFu4GckdfmAsgxoqhO6EiIiEhinpUyMRy8QEREZF8ONiV2p0p0Izg38iIiIjIPhxoS0Wi2KKpUAdEvBiYiIyPAYbkyovFYFlVoDgOGGiIjIWBhuTKh5GbiHgxwyKd96IiIiY+AnrAkV6ZeBcwM/IiIiY2G4MaEr+mXg3MCPiIjIWBhuTOjPZeAcuSEiIjIWbuJnQs27E/PoBSPo8gigbgAkbNQmIrJ2DDcmxN2JjchnrNAVEBGRmeC0lAldqdRt4MfdiYmIiIyH4caEiqp0G/gx3BARERkPp6VMpLqhETXKJgCcljIKZTmg1QAiMSB3E7oaIiISEMONiRRX60ZtHOVS2Mv5thtcZoLuVHC5OxC9TuhqiIhIQJyWMpHia1NSHRy5DJyIiMiYGG5MpKSG4YaIiMgUGG5MpPjaHjcMN0RERMbFcGMizSM3no5sJiYiIjImhhsTKWHPDRERkUkw3JjInyM3DDdERETGxHBjIiXVHLkhIiIyBYYbEylmuCEiIjIJhhsTaFRrUF6rAsBpKSIiImPjVrkmUHqt30YqFsHVTiZwNRYq7C1AqwZEEqErISIigTHcmEBzv42HgxxisUjgaiyUna/QFRARkZngtJQJ8OgFIiIi02G4MQEuAyciIjIdTkuZAEduTKBoF6BRAmI54BUjdDVERCQghhsTKKnRnSvFkRsjyv4MUJYBcneGGyIiK8dpKRPgBn5ERESmw3BjAn9u4MdDM4mIiIyN4cYEOHJDRERkOgw3RqbVavUjN+y5ISIiMj6GGyOramiCqkkDgCM3REREpsBwY2Ql1bqVUo4KKRQ2PBqAiIjI2BhujIxTUkRERKbFcGNkbCYmIiIyLW7iZ2Ql+pEbLgM3Kplry3siIrJaDDdGxpEbE4l4T+gKiIjITHBaysjYc0NERGRagoeb5cuXw9/fHwqFAgMGDEB6evoNH3vy5Ek8+OCD8Pf3h0gkQlJSkukKvU0cuSEiIjItQcPNxo0bkZCQgIULFyIzMxNhYWGIi4tDcXFxq4+vq6tDYGAg3nnnHXh7e5u42tvDcENERGRagoabZcuWYcaMGYiPj0dISAhWrVoFOzs7rF27ttXH9+/fH++++y4eeeQRyOXtIyyU1erCjbt9+6i33cpaDpx8R3dPRERWTbCGYpVKhYyMDLz00kv6a2KxGLGxsUhLSzPYz1EqlVAqlfqvq6qqDPbaf0et0aK8VgUA8HCQmeznWqWyQ4CyDJC7C10JEREJTLCRm9LSUqjVanh5ebW47uXlhcLCQoP9nMTERDg7O+tvfn5+Bnvtv1NRp4JGq/uzqz3DDRERkSkI3lBsbC+99BIqKyv1t7y8PJP97LJrozYudjawkVj8W01ERGQWBJuW8vDwgEQiQVFRUYvrRUVFBm0WlsvlgvXnlNY099tw1IaIiMhUBBtOkMlkiIiIQEpKiv6aRqNBSkoKoqOjhSrLoMpqdCM37g5sJiYiIjIVQXcoTkhIwLRp0xAZGYmoqCgkJSWhtrYW8fHxAICpU6fC19cXiYmJAHRNyKdOndL/OT8/H0ePHoWDgwOCgoIE+z1upOzayA2biYmIiExH0HAzceJElJSUYMGCBSgsLER4eDiSk5P1Tca5ubkQi/8cXCooKEDfvn31Xy9ZsgRLlixBTEwMUlNTTV3+32peKcVl4ERERKYj+NlSs2bNwqxZs1r93v8GFn9/f2i1WhNUZRilzeGGIzdEREQmwyU8RtQ8LcWeGyIiItMRfOTGkjU3FHtwtZTxeQ4DmmoAqYPQlRARkcAYboyorJarpUym63ShKyAiIjPBaSkj0u9zw54bIiIik2G4MRJlkxrVDU0AAA+uliIiIjIZhhsjaV4GLhWL4GTL2T8iIiJT4aeukfy5O7EMIpFI4GqsQPrTgKockLkBUauEroaIiATEkRsj+fNcKU5JmYS6AWiq190TEZFVY7gxkr+O3BAREZHpMNwYSVktTwQnIiISAsONkfBEcCIiImEw3BhJ6bVw48aRGyIiIpNiuDGSq3XNJ4Iz3BAREZkSw42RNO9z48pwQ0REZFIMN0bSPHLDaSkiIiLTYrgxEv3IjR3DDRERkSlxh2IjaFRr9OdKceTGRLrNBNRKQMLVaURE1o7hxgiap6REIsDZ1kbgaqyEe3+hKyAiIjPBaSkjuFrbCEA3JSUR81wpIiIiU2K4MYI/+204akNERGRqnJYyAq6UEkD1eUDTBIilgGOQ0NUQEZGAGG6MgCulBHDiTUBZBsjdgeh1QldDREQC4rSUEVyt5cgNERGRUBhujKC8jrsTExERCYXhxgj0IzecliIiIjI5hhsjKK+7thScIzdEREQmx3BjBH/23HApOBERkakx3BgBV0sREREJh+HGCLjPDRERkXAYbgysoVGNOpUaAHtuiIiIhMBwY2DNozZSsQiOcu6RSEREZGr89DWwspo/97gRiXhopsn0XwlAC4DvORGRtWO4MTB9vw2biU1Lait0BUREZCY4LWVg+pVSXAZOREQkCIYbA6u4toEfV0oREREJg9NSBtY8LeXCaSnTytsEqOsAiR3gN17oaoiISEAMNwbWPHLjYstpKZO6vAlQlgFyd4YbIiIrx2kpA6usvxZu7BhuiIiIhMBwY2AVzdNStpyWIiIiEgLDjYFVXBu5cebIDRERkSAYbgyskj03REREgmK4MbAKfc8Np6WIiIiEwHBjQBqNVt9z48ppKSIiIkEw3BhQjaoJGq3uz06cliIiIhIEw40BVdTqpqRsbSRQ2EgEroaIiMg6cRM/A6qob96dmKM2JufQFZB7ADbOQldCREQCY7gxoObdiZ05JWV6vV8TugIiIjITnJYyoAruTkxERCQ4hhsDquTuxERERIJjuDEg/aGZHLkhIiISDHtuDIhHLwjo+BtAY6WuoZj9N0REVo3hxoCaR25cuTux6dVcAJRlgNxd6EqIiEhgnJYyoD9PBOfIDRERkVAYbgyIq6WIiIiEx3BjQM0jN85cLUVERCQYhhsDquTIDRERkeDMItwsX74c/v7+UCgUGDBgANLT02/6+G+//RY9evSAQqFA7969sXXrVhNVemNarZZLwYmIiMyA4OFm48aNSEhIwMKFC5GZmYmwsDDExcWhuLi41cfv378fkyZNwuOPP44jR45g/PjxGD9+PE6cOGHiyluqVanRdO1IcG7iR0REJBzBw82yZcswY8YMxMfHIyQkBKtWrYKdnR3Wrl3b6uPff/99jB07FvPnz0fPnj3xxhtvoF+/fvjoo49MXHlLzf02MqkYChvB31YiIiKrJeinsEqlQkZGBmJjY/XXxGIxYmNjkZaW1upz0tLSWjweAOLi4m74eKVSiaqqqhY3Y/hzjxsbiEQio/wMIiIi+nuChpvS0lKo1Wp4eXm1uO7l5YXCwsJWn1NYWNimxycmJsLZ2Vl/8/PzM0zx/6O+UQ1HuZRTUkLpNB7wn6S7JyIiq2bxOxS/9NJLSEhI0H9dVVVllIDT398Nx/8dB/W1vhsyMb/xQldARERmQtBw4+HhAYlEgqKiohbXi4qK4O3t3epzvL292/R4uVwOuVxumIJvgUTMKSkiIiIhCTotJZPJEBERgZSUFP01jUaDlJQUREdHt/qc6OjoFo8HgG3btt3w8URERGRdBJ+WSkhIwLRp0xAZGYmoqCgkJSWhtrYW8fHxAICpU6fC19cXiYmJAIDZs2cjJiYGS5cuxbhx4/D111/j8OHDWL16tZC/BhEREZkJwcPNxIkTUVJSggULFqCwsBDh4eFITk7WNw3n5uZCLP5zgGnQoEHYsGEDXn31Vbz88ssIDg7Gpk2bEBoaKtSvQERERGZEpNVqraoDtqqqCs7OzqisrISTk5PQ5RAREdEtaMvnN3ebIyIiIovCcENEREQWheGGiIiILArDDREREVkUhhsiIiKyKAw3REREZFEYboiIiMiiMNwQERGRRWG4ISIiIosi+PELpta8IXNVVZXAlRAREdGtav7cvpWDFawu3FRXVwMA/Pz8BK6EiIiI2qq6uhrOzs43fYzVnS2l0WhQUFAAR0dHiEQig752VVUV/Pz8kJeXx3OrjIjvs2nwfTYNvs+mw/faNIz1Pmu1WlRXV8PHx6fFgdqtsbqRG7FYjE6dOhn1Zzg5OfFfHBPg+2wafJ9Ng++z6fC9Ng1jvM9/N2LTjA3FREREZFEYboiIiMiiMNwYkFwux8KFCyGXy4UuxaLxfTYNvs+mwffZdPhem4Y5vM9W11BMRERElo0jN0RERGRRGG6IiIjIojDcEBERkUVhuCEiIiKLwnBjIMuXL4e/vz8UCgUGDBiA9PR0oUuyOImJiejfvz8cHR3h6emJ8ePH4+zZs0KXZdHeeecdiEQizJkzR+hSLFJ+fj4mT54Md3d32Nraonfv3jh8+LDQZVkUtVqN1157DQEBAbC1tUXXrl3xxhtv3NL5RHRju3fvxr333gsfHx+IRCJs2rSpxfe1Wi0WLFiAjh07wtbWFrGxsTh37pzJ6mO4MYCNGzciISEBCxcuRGZmJsLCwhAXF4fi4mKhS7Mou3btwsyZM3HgwAFs27YNjY2NGDNmDGpra4UuzSIdOnQIH3/8Mfr06SN0KRbp6tWrGDx4MGxsbPDrr7/i1KlTWLp0KVxdXYUuzaIsWrQIK1euxEcffYTTp09j0aJFWLx4MT788EOhS2vXamtrERYWhuXLl7f6/cWLF+ODDz7AqlWrcPDgQdjb2yMuLg4NDQ2mKVBLdywqKko7c+ZM/ddqtVrr4+OjTUxMFLAqy1dcXKwFoN21a5fQpVic6upqbXBwsHbbtm3amJgY7ezZs4UuyeK88MIL2iFDhghdhsUbN26cdvr06S2uPfDAA9pHH31UoIosDwDtjz/+qP9ao9Fovb29te+++67+WkVFhVYul2u/+uork9TEkZs7pFKpkJGRgdjYWP01sViM2NhYpKWlCViZ5ausrAQAuLm5CVyJ5Zk5cybGjRvX4u81GdZPP/2EyMhIPPTQQ/D09ETfvn2xZs0aocuyOIMGDUJKSgqysrIAAMeOHcPevXtx1113CVyZ5bp48SIKCwtb/PfD2dkZAwYMMNnnotUdnGlopaWlUKvV8PLyanHdy8sLZ86cEagqy6fRaDBnzhwMHjwYoaGhQpdjUb7++mtkZmbi0KFDQpdi0bKzs7Fy5UokJCTg5ZdfxqFDh/Dss89CJpNh2rRpQpdnMV588UVUVVWhR48ekEgkUKvVeOutt/Doo48KXZrFKiwsBIBWPxebv2dsDDfULs2cORMnTpzA3r17hS7FouTl5WH27NnYtm0bFAqF0OVYNI1Gg8jISLz99tsAgL59++LEiRNYtWoVw40BffPNN/jyyy+xYcMG9OrVC0ePHsWcOXPg4+PD99mCcVrqDnl4eEAikaCoqKjF9aKiInh7ewtUlWWbNWsWtmzZgp07d6JTp05Cl2NRMjIyUFxcjH79+kEqlUIqlWLXrl344IMPIJVKoVarhS7RYnTs2BEhISEtrvXs2RO5ubkCVWSZ5s+fjxdffBGPPPIIevfujSlTpmDu3LlITEwUujSL1fzZJ+TnIsPNHZLJZIiIiEBKSor+mkajQUpKCqKjowWszPJotVrMmjULP/74I3bs2IGAgAChS7I4o0aNwvHjx3H06FH9LTIyEo8++iiOHj0KiUQidIkWY/DgwddtZZCVlYUuXboIVJFlqqurg1jc8qNOIpFAo9EIVJHlCwgIgLe3d4vPxaqqKhw8eNBkn4ucljKAhIQETJs2DZGRkYiKikJSUhJqa2sRHx8vdGkWZebMmdiwYQM2b94MR0dH/dyts7MzbG1tBa7OMjg6Ol7Xw2Rvbw93d3f2NhnY3LlzMWjQILz99tt4+OGHkZ6ejtWrV2P16tVCl2ZR7r33Xrz11lvo3LkzevXqhSNHjmDZsmWYPn260KW1azU1NTh//rz+64sXL+Lo0aNwc3ND586dMWfOHLz55psIDg5GQEAAXnvtNfj4+GD8+PGmKdAka7KswIcffqjt3LmzViaTaaOiorQHDhwQuiSLA6DV22effSZ0aRaNS8GN5+eff9aGhoZq5XK5tkePHtrVq1cLXZLFqaqq0s6ePVvbuXNnrUKh0AYGBmpfeeUVrVKpFLq0dm3nzp2t/vd42rRpWq1Wtxz8tdde03p5eWnlcrl21KhR2rNnz5qsPpFWy20aiYiIyHKw54aIiIgsCsMNERERWRSGGyIiIrIoDDdERERkURhuiIiIyKIw3BAREZFFYbghIiIii8JwQ0RERBaF4YaIiIgsCsMNERERWRSGGyIiIrIoDDdE1O6VlJTA29sbb7/9tv7a/v37IZPJkJKSImBlRCQEHpxJRBZh69atGD9+PPbv34/u3bsjPDwc999/P5YtWyZ0aURkYgw3RGQxZs6cie3btyMyMhLHjx/HoUOHIJfLhS6LiEyM4YaILEZ9fT1CQ0ORl5eHjIwM9O7dW+iSiEgA7LkhIotx4cIFFBQUQKPRICcnR+hyiEggHLkhIougUqkQFRWF8PBwdO/eHUlJSTh+/Dg8PT2FLo2ITIzhhogswvz58/Hdd9/h2LFjcHBwQExMDJydnbFlyxahSyMiE+O0FBG1e6mpqUhKSsL69evh5OQEsViM9evXY8+ePVi5cqXQ5RGRiXHkhoiIiCwKR26IiIjIojDcEBERkUVhuCEiIiKLwnBDREREFoXhhoiIiCwKww0RERFZFIYbIiIisigMN0RERGRRGG6IiIjIojDcEBERkUVhuCEiIiKL8v/axZ0jP8PF/gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x max =  2.1\n"
     ]
    }
   ],
   "source": [
    "# Define x range and initial conditions\n",
    "initial_conditions = [1, -1.5886]  # [y(0), y'(0)]\n",
    "\n",
    "def f(x): \n",
    "    return (1+4/3*x**(3/2))\n",
    "\n",
    "def F(x_range):\n",
    "    solution = solve_differential_equation(DGL_system, x_range, initial_conditions, plot_result=False)\n",
    "    return solution.t, solution.y[0]*f(solution.t)\n",
    "\n",
    "x, Fx = F(x_range)\n",
    "\n",
    "max_index = np.argmax(Fx*x)\n",
    "x_max = x[max_index]\n",
    "Fx_max = Fx[max_index]\n",
    "\n",
    "plt.plot(x, Fx*x)\n",
    "plt.vlines(x_max,0,0.5, linestyle='dashed', linewidth=2, alpha = 0.7, color = 'orange', label='x max')\n",
    "plt.legend()\n",
    "plt.ylabel('F(x) x')\n",
    "plt.xlabel('x')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "print('x max = ', np.round(x_max,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fe341d48-d844-4dfa-ac53-a2ecfb523e89",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "l = 1, Z =  3.5\n",
      "l = 2, Z =  18.4\n",
      "l = 3, Z =  52.1\n"
     ]
    }
   ],
   "source": [
    "l = [1,2,3]\n",
    "\n",
    "for l in l: \n",
    "    print(r'l = {}, Z = '.format(l), np.round(4/(3*np.pi)*(l*(l+1)/(Fx_max*x_max))**(3/2), 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5169e24d-13ab-4d5b-8d7a-045ae339309c",
   "metadata": {},
   "source": [
    "From the periodic table we know experimentally that the values of $Z$ at which the first $p, d$ and $f$ atoms appear are $Z = 5, 21, 58$. Even though these values do not agree with the predictions from the model we obtained here, the Thomas-Fermi model still predicts the right trend of the dependence of $Z$ on $\\ell$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d07ee9c-0833-40f5-8e8f-59b3abbf710e",
   "metadata": {},
   "source": [
    "e) One drawback of the Thomas-Fermi model is that it includes self-screening and therefore overestimates the total screening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79ae5021-bf1c-4b6e-bba3-2496052aa7f8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}