{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c6ee499b-1d71-4ecd-9e58-0d5e40c411ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Markdown\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# Electron & proton masses, Boltzmann & Planck constants\n",
    "from astropy.constants import m_e, m_p, k_B, h\n",
    "# Useful units\n",
    "from astropy.units import eV, kg, m, Kelvin\n",
    "\n",
    "# Ionization energies\n",
    "I1 = 13.598 * eV\n",
    "I2 = 24.587 * eV\n",
    "I3 = 54.416 * eV\n",
    "\n",
    "# Partition function ratios\n",
    "omega1 = 1\n",
    "omega2 = 4\n",
    "omega3 = 1\n",
    "\n",
    "# Species Masses\n",
    "m_H = m_p + m_e\n",
    "m_He = 4 * m_H\n",
    "\n",
    "# Number abundances\n",
    "nu_1 = 0.9\n",
    "nu_2 = 0.1\n",
    "\n",
    "# Density and temperature of the medium\n",
    "rho = 1e-7 * kg / m**3\n",
    "T = 12500 * Kelvin\n",
    "\n",
    "# And finally, beta (=1 since we neglect radiation pressure)\n",
    "beta = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9d10423a-9460-4d55-bc7c-db4c7806c866",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "The value of $\\mu_0$ is approximately 1.3000"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# We give mu0 in atomic mass units\n",
    "mu_0 = (nu_1 * m_H + nu_2 * m_He) / m_H\n",
    "Markdown(rf\"The value of $\\mu_0$ is approximately {mu_0:.04f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3f91a9c4-acf8-4ca1-bd17-65ff5c0a2d51",
   "metadata": {},
   "outputs": [],
   "source": [
    "def P(E):\n",
    "    \"\"\" Eq (8) of the question sheet \"\"\"\n",
    "    return k_B / m_H * (1 + E) / mu_0 * rho * T\n",
    "\n",
    "def K(E, I_i, omega_i):\n",
    "    \"\"\" Eq (6) of the question sheet \"\"\"\n",
    "    numerator = (2 * np.pi * m_e)**(3/2) * (k_B * T)**(5/2)\n",
    "    denominator = beta * P(E) * h**3\n",
    "    return numerator / denominator * omega_i * np.exp(-I_i / (k_B*T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c0a85377-e623-4ad9-a499-a996632fa15a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def getE(X1, X2, X3):\n",
    "    \"\"\" Eq (1) of the solution sheet\"\"\"\n",
    "    return nu_1 * X1 + nu_2 * (X2 + 2 * X3)\n",
    "    \n",
    "def getalpha(E):\n",
    "    \"\"\" alpha defined as in the solution sheet \"\"\"\n",
    "    return (1 + E) / E\n",
    "    \n",
    "def getXs(X1, X2, X3):\n",
    "    \"\"\" Here, we update X1, X2 and X3 based on the current guess \"\"\"\n",
    "    # Compute E:\n",
    "    E = getE(X1, X2, X3)\n",
    "    \n",
    "    # Obtain the {Ki}:\n",
    "    K1 = K(E, I1, omega1)\n",
    "    K2 = K(E, I2, omega2)\n",
    "    K3 = K(E, I3, omega3)\n",
    "    alpha = getalpha(E)\n",
    "    X1new = alpha * K1 / (1 + alpha*K1)\n",
    "    X2new = alpha * K2 / (1 + alpha*K2*(1 + alpha*K3))\n",
    "    X3new = alpha**2 * K2 * K3 / (1 + alpha*K2*(1 + alpha*K3))\n",
    "    \n",
    "    return X1new.value, X2new.value, X3new.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5f678de3-8906-499b-9f5f-1623dbe0854d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Starting guess\n",
    "X1 = 1\n",
    "X2 = 0.0\n",
    "X3 = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "887f71ab-fcf0-4417-a766-8af9e134b2c0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "After one iteration: $X_1=0.99629$, $X_2=0.03834$, $X_3=0.00000$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X1, X2, X3 = getXs(X1, X2, X3)\n",
    "Markdown(f\"After one iteration: $X_1={X1:.05f}$, $X_2={X2:.05f}$, $X_3={X3:.05f}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ada62ed3-844b-494c-863c-23b4f55759e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "After two iterations: $X_1=0.99629$, $X_2=0.03832$, $X_3=0.00000$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X1, X2, X3 = getXs(X1, X2, X3)\n",
    "Markdown(f\"After two iterations: $X_1={X1:.05f}$, $X_2={X2:.05f}$, $X_3={X3:.05f}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b1998dae-14e3-452b-a2bb-5d40e4ac557e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "After three iterations: $X_1=0.99629$, $X_2=0.03832$, $X_3=0.00000$"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X1, X2, X3 = getXs(X1, X2, X3)\n",
    "Markdown(f\"After three iterations: $X_1={X1:.05f}$, $X_2={X2:.05f}$, $X_3={X3:.05f}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fff0fb2e-215a-418c-86fb-890c5345d77e",
   "metadata": {},
   "source": [
    "We see that in 3 iterations, we already converge to a level of precision reaching $10^{-5}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4e12e477-e717-4155-9c35-45497aec3a6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can automate this process using a simple function\n",
    "def doIterations(niteration=5):\n",
    "    # Starting guess\n",
    "    X1 = 1\n",
    "    X2 = X3 = 0\n",
    "    E = 0.9\n",
    "    \n",
    "    # Iterate:\n",
    "    for i in range(niteration):\n",
    "        X1, X2, X3 = getXs(X1, X2, X3)\n",
    "        \n",
    "    return X1, X2, X3, E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f07bc5c7-d691-49a0-a56a-92fba33e066c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We are interested in particular in helium, so we\n",
    "# will store its ionization fraction (both HeII and HeIII)\n",
    "# over a range of temperatures\n",
    "X2s, X3s = [], []\n",
    "Ts = np.linspace(5000, 40000)\n",
    "\n",
    "# Compute the result for each temperature:\n",
    "for T in Ts:\n",
    "    T = T * Kelvin\n",
    "    X1, X2, X3, E = doIterations()\n",
    "    X2s.append(X2)\n",
    "    X3s.append(X3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a5d3823e-56d5-47d0-983b-a961f392504e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x109fcb230>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcVFJREFUeJzt3Xd8U/X6wPFPkjbpXnSXlhZa9l61CAJSRfQievXKxQW4/blxgQvH1brFqyheUXGDoiLKUECm7L1HaaFldFG6R9rk/P5IGyi00ELak/G8X+bVk5NvznlOEpOH79QoiqIghBBCCOEktGoHIIQQQghhS5LcCCGEEMKpSHIjhBBCCKciyY0QQgghnIokN0IIIYRwKpLcCCGEEMKpSHIjhBBCCKfipnYALc1sNnPs2DF8fX3RaDRqhyOEEEKIRlAUheLiYiIjI9Fqz10343LJzbFjx4iOjlY7DCGEEEJcgMzMTFq3bn3OMi6X3Pj6+gKWF8fPz0/laIQQQgjRGEVFRURHR1t/x8/F5ZKb2qYoPz8/SW6EEEIIB9OYLiXSoVgIIYQQTkWSGyGEEEI4FUluhBBCCOFUXK7PTWOZTCaqqqrUDkM0E3d3d3Q6ndphCCGEaAaS3JxBURSysrIoKChQOxTRzAICAggPD5f5joQQwslIcnOG2sQmNDQULy8v+eFzQoqiUFZWRk5ODgAREREqRySEEMKWJLk5jclksiY2rVq1Ujsc0Yw8PT0ByMnJITQ0VJqohBDCiUiH4tPU9rHx8vJSORLREmrfZ+lbJYQQzkWSm3pIU5RrkPdZCCGckyQ3QgghhHAqqiY3K1asYOTIkURGRqLRaJgzZ855n7Ns2TJ69+6NwWAgPj6eGTNmNHucQgghhHAcqiY3paWl9OjRg6lTpzaqfHp6Otdccw1Dhw5l69atPProo9x111388ccfzRypEEIIIRyFqqOlRowYwYgRIxpdftq0acTFxfHOO+8A0KlTJ1atWsV7773H8OHDmytMAdYasnHjxqkah3A9ZrNCtVlBqwE3neXfYyazQmF5FWZFwawoYPkPy6ZS89cy7L+WRqMhxMeA3s1yjNLKasqMJjzctfh6uENN+cpqM25aDTqtRvplCeGgHGoo+Jo1a0hOTq6zb/jw4Tz66KMNPqeyspLKykrr/aKiouYKTwiXYjIr6LSnfvy3ZJzkQE4JXSP96RzpB8ChvFI+WZFGRZWJMmM15VVmKowmKqpNGKvNVJsVqkxmqqrNVNVsV5sUFk8YTLi/BwCvzt/DZ6vSuX9IO56+qiMAmfllDHl7WZNjnvvgpXRvHQDAl2sO8ebCffyrT2ve+lcPAMqrTHR+4VRNsJtWg5tOg5tWi5tOg7tOi6e7Di+9Dk+9zrr90OUJ9Ii2HDc1p5jl+/NICPXhsvYh1mOZzQparSRLQrQEh0pusrKyCAsLq7MvLCyMoqIiysvLrXOXnC4lJYWXXnqppUJU1bhx4ygoKDir79KyZcsYOnQoJ0+eJCAgoEVjGj9+PFFRUfznP/856zF7jFdYKIpCXomR9LxSDuWVkpZXSnpeCbnFlZwsq+JkmZGi8ir2vjLCWhPy5epDzNl6jEkjOlqTm4LyKr5fn9Hk81eZzNbt2nzAbD5VC3N6UqXRgIZTo9801n2WBzRYanFQQFtPTczpu6pPO0ftfcs+M+dye1KsdXvDoZO88vtuhnUMrZPc9HjpTwzuWoJ9DLTy0RPbypsO4b60D/OlQ5gvgd76c55DCNF4DpXcXIhJkyYxYcIE6/2ioiKio6NVjMhxGI1G+vfvD0B+fj4AU6ZMAWD9+vXo9ef+MjaZTPz+++/MmzevWeMUF+9oQTk/bswkPa/Ucsstpbiy+rzPKyg3EuprqWHpGuVPQXkVMUGn5omKDPBgwhXt8XQ/VdPhqddhcNOid9PiptWid7PUjLjrtLjX1I7U1toAPHZFex4cmoDB/VQXwdaBnqS9dvVF1YT835B4/m9IfJ19vgY3dr40nGqTmSqTgqm2NsmsYDKbqaw219RCWW612+3DfK3HiArw5B/dI+gW5W/dV2aspriymuJKyCsxAvB36ok65w7xNdAhrCbZCfehfZgvnSL88HCXCSaFaCqHSm7Cw8PJzs6usy87Oxs/P796a20ADAYDBoOhJcJzGGazmTfeeIP//e9/ZGVl0b59e55//nluvPHGOuX0ej1bt24FLqzPzerVq3F3d6dfv34tFrNonIO5JWg1GuKCvQEoKDMyZfGBOmU0GssPdVywN22DvYkN9ibC35NAL3eCvPUEeOlpdVptw12D2nLXoLZ1jhHq68HDwxIuKlYvvRuckUdrNBqaozuMRqPBx3BxX4uXtQ+pU2MD4OmuY/2zwzhRYiSvpJKcokoO5pawP7uYfdnFZOaXk1tcSW5xJatS86zP+/DmXvyje+RFxSOEK3Ko5CYpKYn58+fX2bdo0SKSkpKa/dxlxvP/K/ZMep3W2gGy2mTGaDKj1Wjq/EusoeN66ZvvrUlJSeGbb75h2rRpJCQksGLFCm699VZCQkIYPHiwzc4zd+5cRo4caZNOmS0VsyuYsng/UxYf4Ka+rXnzRktfk7hgb0b3jSYuxNuazEQHeUmtgY1oNBpCfT2stVxnKq2s5kBOCfuzLMnO/uxi0nJLGdoh1Frm42UHWZt2gjsGxjH4jORJCFGXqslNSUkJqamp1vvp6els3bqVoKAgYmJimDRpEkePHuWrr74C4L777uPDDz/kqaee4o477uCvv/7ihx9+aJFmj9M7GTbW1Jt7c013y6KMf+zK5oHvNpMYF8Sse08lYwPfWEp+qfGs5x56/ZoLivP333/Hx8enzj6TyWTdrqys5LXXXmPx4sXWpLBt27asWrWKTz75xKaJwq+//sp77713UfG2dMzOZvexIj5ZcZDbLmlD39ggAPq0CcRdp8FYfaofiZfejTdu7K5WmC7P2+BGz+gAetZ0Sq7Pb9uOsft4kfU7BeBkqZGKahMR/vXXXAuhClM16NStO1H17Bs3bmTo0KHW+7V9Y8aOHcuMGTM4fvw4GRmnOiPGxcUxb948HnvsMd5//31at27N9OnTZRj4aYYOHcrHH39cZ9+6deu49dZbAUhNTaWsrIwrrriiThmj0UivXr0aPG5Th4Dv2bOHY8eOMWzYsIuK92JidmUVVSbeX3KA/61Iw2RWCPYxWJObS9q2YvPzV1iHPwvH8N8xPfljVzbJnU4Nqvhp8xHeWLiX25NiefjyBPy95D0VKijIgMOr4fDflr+t+8P1H5//ec1I1eRmyJAhdeahOFN9sw8PGTKELVu2NGNU9dv9ctMTKL3uVAfI4V3C2P3y8LNGa6x6euiZT7so3t7exMfX7SR55MgR63ZJSQkA8+bNIyoqqk45W/ZNmjt3LldccQUeHvVXw9c6X7zQcjE7i9UH83jm5x0cOlEGwIiu4YzqearfhqXjrqy84mjiQ32JD/Wts+/wiTKqTAqfrUrn581HeDS5PTcnxsj7K5qPokD5SfCy/GMJYxn8txeYT+tiYW56Nw5bc6g+N2q62D4wbqf1v7HlcZuqc+fOGAwGMjIymrU559dff+Wee+6xybFaKmZHV1hWxavzd/PDRktyGOZn4JVRXbmyS7jKkYnm8sp1XUnuHMZ/ft/NgZwSJs/dxVdrDvHcNZ0Z0iFEJiEUtmGqhkMrYdfPsG8B+LeGe5ZZHtN7WWpqTEZoMwDaXAoxiaqGC5LcuBxfX1+eeOIJHnvsMcxmMwMHDqSwsJC///4bPz8/xo4de9HnyMnJYePGjcydO9cGEbdMzI5MURTm78hi8txd5JVYJqy89ZIYnrqqI37S9OT0BrcP4dJHBvH9hkzeW7Sfg7mljJ+xgUEJwTx3TWc6hPue/yBCnMlssjQx7foZds+FslOj+Kgqh6oKcK+pmR83D7T2VVsoyY0LeuWVVwgJCSElJYW0tDQCAgLo3bs3zzzzjE2O/9tvv9G/f3+Cg4Ntcjxo/pgd1fHCcp6fs5PFe3IAaBfizes3dKdfTf8a4RrcdFpuu6QN1/aI5KOlqXzx9yFWHshjxPsrGNM/hseuaE+wjzThiiaY83+wfeap+55B0Pla6HI9xCSB22mfJztLbAA0yrk6vTihoqIi/P39KSwsxM/Pr85jFRUVpKenExcXd96+IqJh1157LQMHDuSpp55SO5RzcvT3+/ftx5j40w5KKqtx12m4f0g8Dwxth8FNhm+7usMnSnl9wV4W7MwCwNfDjWm39uHSeNv9g0M4kbJ8WPcJ9BwDgbGWfdt/gPlPQqd/QJd/QtxloFO3Jvhcv99nkpobYXMDBw5kzJgxaofh1A7mljDhh20Yq830jgng9Ru615klV7i2Nq28+fjWPqxLO8Er83az82gR42dsYMb4fgxoJwmOqFF0HNZ8CBu/gKpSSzPTwMcsj3W+znJzc8xlQSS5ETZn7zU2ziC2lTcPDo1n+5FCPrmtT521loSoldi2FT/dP4CHv9/C4RNldI449792hYvIT4e/34et31o6AgOEd7PcajloUlNLkhshHJBOq+HhYQmy0rQ4L4Objqk396a4opoAL8f+wRIXSVHg1wdh23eg1EziGX0JXPYExCfTLGuaqESSGyEcyO5jRbQN8bYuiyCJjWgMN522zqrjX605hKLA2AGx6gUlWp5GY7kpZksyM+hxy/BtJyTJjRAOIruoglumryXMz4MvxveTKffFBdmScZIXft0FQOdIPxlZ58yqK2HF29DrVghsY9k3+CnodxdE9lQ1tOYmyY0QDuLIyXK0Gg1ajYYgb2leEBemZ3QAD10eT7nRRN82gWqHI5rLsS2W4dw5uyFzLdw+11JrExBjuTk5SW6EcBB92gSy4NFBlBtNMtxbXDCNRsPjV3ZAURTrDMYVVSYMblqZ0dgZVBthxVuw8h1QTOAVDP3udqr+NI0hyY0Qdu70H6FQX8ebj0fYp9rPVGW1ibu/2kjrQE/+c103GXnnyI5vhzn3Q/ZOy/0u18PVb4O36w3/l+RGCDtWWlnNzdPXcc+gtlzTPULtcIQT2pB+kr9T8zArEOSt58nhHdUOSVyItGXwzQ2WRSs9g+Cad6DrP9WOSjX2N2eyEMLqpd92sS2zgFfn7abcaFI7HOGEBiYE89aNPQCYtjyNHUcKVY5IXJDoS6BVPHT8BzywzqUTG5DkRgi79fv2Y/yw8QgaDbw3uieeeulnI5rHDX1a84/uEZjMCk/Otsx8LRxA9i7L3DVgmV14/AIY/Q34hKoblx2Q5EYIO3TkZBmTft4BwIND40ls20rliISze+naLrTy1rM3q5gPl6aqHY44n50/wyeXwZKXTu3zCnK5jsMNkeRGnNOMGTOYMWOGzcqJxnnix20UV1TTKyaAh4clqB2OcAGtfAy8NKoLAB8tTWXXMWmeslubv4af7rT0rynIBLM0WZ9Jkhsh7Mz+7GLWpuXjrtPw/uheuOvkf1PRMq7pFsFVXcKpNis8+eN2qkzSPGV31nwEcx+0zDLceyz883+glSbrM8m3phMZN24c11133Vn7ly1bhkajoaCgoMVjOp/x48fz3HPPNfi4I17TxZq3/TgAlyWEENPKS+VohCvRaDS8cl1XArzc2X28iI+XHVQ7JFFLUWD5m/DHJMv9pAdh5PuS2DRAhoKLsxiNRvr37w9Afn4+AFOmTAFg/fr16PX6JpVriMlk4vfff2fevHm2vgSHNm+HJbmRod9CDSG+Bl4c2YVHZ23lg78OcGWXMDqGy2riqlv8Ivw9xbI95BnLMgrSv6ZBUnPjgsxmMykpKcTFxeHp6UmPHj2YPXu29XG9Xs/WrVvZunUrL7/8Mi+//LL1/ukJS2PLNWT16tW4u7vTr1+/Zr8mR7E/u5jUnBL0Oi3JncPUDke4qFE9I0nuFEaVydI8VS3NU+oL7QRoYPhrMORpSWzOQ2puGstY2vBjGp1lGF6jymrB3fP8ZfXeTYuvCVJSUvjmm2+YNm0aCQkJrFixgltvvZWQkBAGDx7cbOc909y5cxk5cqRNpny3l2u6WL/XNkm1D8bPw13laISr0mg0vHZ9V9ann2DH0UK++PsQd1/WVu2wXFuPf0Nkbwhpr3YkDkGSm8Z6LbLhxxKuhFt+PHX/rXioKqu/bJuBMP60Zpgp3aDsxNnlXrywkQq///47Pj4+dfaZTKd60ldWVvLaa6+xePFikpKSAGjbti2rVq3ik08+adFE4Ndff+W99947bzlHuqaLoSgK82uapK7uJk1SQl2hfh68MLILS/fmcH3vKLXDcT3VRvjzWRj0BPjW1OJKYtNoktw4maFDh/Lxxx/X2bdu3TpuvfVWAFJTUykrK+OKK66oU8ZoNNKrV6+zjjdu3LhGnbex5Wrt2bOHY8eOMWzYsPOWtfU12av92SXSJCXsyg29o7ixT2u1w3BNfz4H6/8HRzbCXUtAK71ImkKSm8Z65ljDj2nO6K3+5DkmwNKc8QF9dMeFx1QPb29v4uPj6+w7cuSIdbukpASAefPmERVV919jBoPBprGcy9y5c7niiivw8Dj/QpCOck0XK9zfg5R/diOrsEKapIRdOLPJ+PCJUtq0ar4mc1Fj1y+w/hPL9uCnJLG5AJLcNFZT+sA0V1kb6Ny5MwaDgYyMDFWba3799VfuuecemxzLXq7pYvl7ujOmf4zaYQhxlnKjiSd+3MbiPdkseGQQbUN8zv8kcWFOHIRfH7JsD3wMOoxQNx4HJcmNi/H19eWJJ57gsccew2w2M3DgQAoLC/n777/x8/Nj7NixzR5DTk4OGzduZO7cuTY5nj1ckxDOzMNdS2F5FVUmM4t2Z3PvYElumkVVOfwwFozFEDMAhjY8B5g4N0luXNArr7xCSEgIKSkppKWlERAQQO/evXnmmWda5Py//fYb/fv3Jzg42GbHVPuaLtaPGzMpM5q4ulsEIb6O05QmXINGo+Htf/WgymQmOkgmlmw2C56G7B3gFQw3fg46+Ym+UBpFqV1S1DUUFRXh7+9PYWEhfn51J6aqqKggPT2duLi4RvUFERfm2muvZeDAgTz11FOqxmEv77eiKAx7dzlpuaVMGd2T63rJyBQhXE55AXx6OeSnwW2/QLuhakdkd871+30mSQtFixs4cCBjxoxROwy7YTIr3JLYhj92ZTGsU6ja4QhxXul5pfh7uhPkff7JOkUjeQbAPcsgbakkNjYgNTensZd/yYuWIe+3EE339h/7mLoslYeGxjPhyg5qhyNcSFNqbmR8mRBCiEbrHOmHosDXaw9TUWU6/xNEwxQFfn0A1k6zbAubkeRGCBWl55Xyw4ZMCsqMaociRKNc2TmM1oGenCyr4qfNR87/BNGwzV/Clm/gj2cgZ4/a0TgVSW6EUNHPm4/w1E/befqn7WqHIkSjuOm0jL80DoDPVqVjNkuNwwU5vh3m1wyqGPY8hHVWNx4nI8mNECpRFIV522UtKeF4RveLxtfgRlpuKUv35agdjuOpKIQfx4KpEhKGw4BH1I7I6UhyUw8X62PtstR+n/ccLyYtrxS9m5ZhnWQtKeE4fAxujEm0zKY9fWW6ytE4oCUvW4Z8+0fD9dNkeYVmIK/oadzdLev5lJU1sKK3cCq173Pt+97SalcAH9I+BB+DzMogHMu4AbHotBrWpJ1g59FCtcNxHLn7YeMXlu3rPgKvIHXjcVLyjXoanU5HQEAAOTmWalYvL6+zFo4Tjk9RFMrKysjJySEgIACdTnf+JzVDDPNqkptrukuTlHA8kQGeXNMtgrnbjvHZqnTeG91T7ZAcQ+Y60GggYQTEXaZ2NE5LkpszhIeHA1gTHOG8AgICrO93S9t9vIj0vFIM0iQlHNhdg+KYu+0Yv207xlNXdSDC31PtkOxf79sgJgm0Lf+PKlciyc0ZNBoNERERhIaGUlVVpXY4opm4u7urUmNTy9ok1UGapITj6t46gP5xQaxPz+fL1YeZOKKj2iE5huB4tSNwevKt2gCdTqfqj59wXqePkrqme6TK0Qhxce4e1Jb16fl8t+4wD10ej7ck6/VLXwGegRDeTe1IXIJ0KBaihe06VsShE2WWJqmOspaUcGzDOoYS28qLoopqFu3OVjsc+1RdaZmJeNog2Dtf7WhcgqTYQrSw2iapoR1C5V+5wuFptRpevLYLXno3+sUGqh2OfVr/KRRkgG8EtB2idjQuQb5ZhWhBMkpKOKMhHaQGskFl+bDiTcv20GdB76VuPC5CmqWEaEGpOSUcrmmSulyapIQTKqmsVjsE+7LyHcuMxKFdoOfNakfjMqTmRogWFOxj4O1/9aC0slqapITTeX3BXr5ec4iv70qkd4w0UZGfDuv/Z9m+8mUZ/t2C5NtViBYU6K3nxj6t1Q5DiGaRV1JJqdHEvO3HJbkByzILJiO0HQrxyWpH41IkuRFCCGET/zekHSN7RHJZQrDaoahPUSyT9R1aCVe+onY0LkejqL16YAsrKirC39+fwsJC/Pz81A5HuJgfN2bSOtCL3m0CMLhJFbUQTq+qAtw91I7CKTTl91tqboRoIWXGap6cvR2ATc8lY/CR5EY4L0VRZG0+kMRGJTJaSogWUlxRTXKnMHq09qeVj0HtcIRoFsZqM8/P2cnAN5ZSWO6CS9iYquDbm2D3XEvTlFCF1NwI0ULC/DyYPrav2mEI0az0blrWpZ/gaEE5i3dnc4OrdaDfNAMO/AFHN0G7y8Hgo3ZELklqboQQQtjU1d0sE1TWzsbtMswm+Pu/lu0hEyWxUZEkN0K0kMKyKlys/75wUbXJzcoDeRRVuFDT1IFFUJgBHgHQ61a1o3FpktwI0UKunLKcXq8sYl9WsdqhCNGs2of5Eh/qg9FkZskeF1pMc8N0y99et4K7p7qxuDhJboRoAQVlRrKLKikoqyIyQEZPCOd3dddwAObvyFI5khaSnwapiy3bfe9QNxahfnIzdepUYmNj8fDwIDExkfXr15+z/JQpU+jQoQOenp5ER0fz2GOPUVFR0ULRCnFh9meXABAV4Imvh7vK0QjR/K6uWRh2+f5cil2haWrjF4BimYm4VTu1o3F5qiY3s2bNYsKECUyePJnNmzfTo0cPhg8fTk5OTr3lv/vuOyZOnMjkyZPZs2cPn332GbNmzeKZZ55p4ciFaJp92ZamqPZh0sFQuIYOYb60DfbGWG3mr731f6c7ldiBEHcZ9LtL7UgEKic37777LnfffTfjx4+nc+fOTJs2DS8vLz7//PN6y69evZpLL72Um2++mdjYWK688krGjBlz3toeIdS2v6afTftwX5UjEaJlaDQa1xo11X44jP0NOoxQOxKBismN0Whk06ZNJCefWkxMq9WSnJzMmjVr6n3OgAED2LRpkzWZSUtLY/78+Vx99dUNnqeyspKioqI6NyFaWm3NTYcwSW6E6xjRzdLvZtm+XEorq1WORrgS1ZKbvLw8TCYTYWFhdfaHhYWRlVV/B7Sbb76Zl19+mYEDB+Lu7k67du0YMmTIOZulUlJS8Pf3t96io6Nteh1CnI+iKOy3NktJciNcR+cIP2JbeVHpzE1TWTtg+VtQ7EKjwhyA6h2Km2LZsmW89tprfPTRR2zevJmff/6ZefPm8corDa+4OmnSJAoLC623zMzMFoxYCMgttoyS0mogPlT63AjXodFoGFHTNLVgp5M2Ta37BJb+B/58Vu1IxGlUW34hODgYnU5HdnbdbDc7O5vw8PB6n/P8889z2223cdddlg5b3bp1o7S0lHvuuYdnn30WrfbsXM1gMGAwyDo+Qj21TVKxrbzxcJfFMoVruaZbBB8vO0h2UaXzLaZZfhJ2zLZs971T3VhEHarV3Oj1evr06cOSJUus+8xmM0uWLCEpKane55SVlZ2VwOh0lh8LmflV2KvaSfukSUq4oi6Rfix/cgg/3T/AuRIbgK3fQ3U5hHWFmEvUjkacRtWFMydMmMDYsWPp27cv/fv3Z8qUKZSWljJ+/HgAbr/9dqKiokhJSQFg5MiRvPvuu/Tq1YvExERSU1N5/vnnGTlypDXJEcLeWPvbyEgp4YI0Gg1tWnmrHYbtmc2nZiTudyc4W+Lm4FRNbkaPHk1ubi4vvPACWVlZ9OzZk4ULF1o7GWdkZNSpqXnuuefQaDQ899xzHD16lJCQEEaOHMmrr76q1iUIcV77aibwk5FSwtWVVlbjptNgcHOCf4ymL4f8g6D3hW43qR2NOINGcbH2nKKiIvz9/SksLMTPz0/tcISTM5sVur74B2VGE4snXEZ8qCQ4wjVN/nUnMzdk8v6/e3JV1wi1w7l4M2+Bvb9D/3vg6rfUjsYlNOX3W9WaGyGcnVarYfmTQ9mfXeycVfNCNJLeTUtltZm1afmOn9woCni1AjdP6Uhsp6TmRgghRLPLzC+juKKaThG+ztOxuLIYDFIb21Kk5kYIIYRdiQ7yUjsE25PExm5JciNEM/p0RRonSo1c3yuKDjJaSgjA0hdNq3XQ2psjm0DnDhHd1Y5EnINDzVAshKP5ectRpi0/SGZ+mdqhCKG6oooqHpu1lYFv/EVltUntcC7Mohfgk0Gw4TO1IxHnIMmNEM1obFIbbrukDV2ipH+XED56N9YcPMGxwgpW7s9TO5ymy9kDh1eBRgftr1I7GnEOktwI0Yz+3T+GV67rSoS/p9qhCKE6rVbDVV0ty+vMd8S1pmprazpeDf5R6sYizkmSGyGEEC3mmu6WYeCLdmc7VtNUdSVsn2XZluHfdk+SGyGaya5jhWzLLKDMWK12KELYjT4xgYT4GiiuqGbToZNqh9N4acugsgh8IyBusNrRiPOQ5EaIZvLhX6mMmvo3363LUDsUIeyGVqthUHwwACtTHajfze65lr8d/wFa+em0d/IOCdFM9tUsmClDwIWo69Ka5GbVAQdJbhQFMtZYtjtfq24solFknhshmkFFlYlDeaWALJgpxJkGJliSm53HCjlZaiTQW69yROeh0cAD6+DQKogZoHY0ohGk5kaIZnAwtwSzAgFe7oT4GtQORwi7EubnQfswHxQFVh88oXY4jaNzh3ZDQSd1Ao5AkhshmsH+miap9mFOtI6OEDY0MD4EgFWpuSpHch6KYrkJhyLJjRDNYF9WCSBNUkI0ZFBN09TKA3nY9frNGWtgSjdY9rrakYgmkORGiGZgrbmRzsRC1CuxbRDuOg1HTpZz+IQdL0+yey4UZsLJQ2pHIppAkhshmsG+rJqRUlJzI0S9vPRu9I4JBGDH0UKVo2mA2Qx7aoaAd5JRUo5EekYJYWMlldUcLSgHoH2Yj8rRCGG//nNdV4K89bTysdNO98c2Q9FR0PtAu8vVjkY0gSQ3QtjYgZomqTA/AwFedj7EVQgVJdh7zebuXy1/E64Edw91YxFNIs1SQtjY6SOlhBAOSlFONUnJxH0OR5IbIWxMRkoJ0XgLdx5n9CdrmLb8oNqh1JW1w9KJ2M0D4q9QOxrRRNIsJYSNVVSb0LtpZaSUEI2QV2JkXXo+CnDf4HZqh3OK3hv63Q2KCQzSd87RaBS7nmDA9oqKivD396ewsBA/Pz+1wxFOqtpkxqQoGNx0aocihF07XljO4t3ZDEoIITbYW+1whB1ryu+31NwI0QzcdFr5n0uIRojw9+S2pFi1wxBORvrcCCGEEKdLXWxZJNNsUjsScYEkuRHChr5fn8GI91cyfWWa2qEI4TCKK6r4eu1hXpy7S+1QLBa9CDOuge2z1I5EXCBJboSwoZ1HC9lzvIj8UqPaoQjhMKpNCi/8upMZqw+RXVShbjD5aZC9AzQ6aH+VurGICybJjRA29MDQeD69vS+jekapHYoQDiPQW0/XSH8AVh3IUzeY3TVz28QNAq8gdWMRF0ySGyFsKDLAkys6h9FBhoEL0SQDa1YJX5WqcnIja0k5BUluhBBCqG5Q/KnkRrUZSgqPwNFNgAY6/kOdGIRNSHIjhI3sOlbIB0sOsFrtf3kK4YD6xAbi4a4lt7iS/dkl6gSx5zfL35gk8A1TJwZhE5LcCGEj69PzeWfRfr5cc0jtUIRwOAY3Hf3jWgGw8kCuOkEcXm35K2tJOTxJboSwkeOFllEeUQFeKkcihGM6vWlKFTd9BXf9Bd3+pc75hc3IJKpC2MixgnIAIgM8VI5ECMdU26l4XVo+ldWmll++RKOB1n1a9pyiWUjNjRA2UltzE+HvqXIkQjimjuG+BPsYKK8ysflwQcue3Gxu2fOJZiXJjRA2crym5iZCam6EuCAajYaB8ZZ+N6tSW7DfTVk+vNMBfn0QTFUtd17RbCS5EcIGqk1msmpmVo2UmhshLtiltf1uWnIyv/TlUJoDRzaAzr3lziuajSQ3QthATnElZgXctBpCfA1qhyOEwxqUEALA9qOFFJS10DImB/+y/G03rGXOJ5qddCgWwgaOF1qapML8PNBpNSpHI4TjCvf34OHL4+kU4YeHewt0KFYUOLjMst1uaPOfT7QISW6EsIFjBTVNUtLfRoiLNuHKDi13shMHoTADdHpoM6DlziualTRLCWEDtTU3MlJKCAeTttTyNzoR9N7qxiJsRpIbIWygtuZGRkoJYRtbMk7y/uIDHK0ZhdhsrP1tpEnKmUizlBA2UFtzExUgNTdC2ELK/L2sP5RPqJ+BMf1jmu9EET2gIAPaXd585xAtTpIbIWwgwt+T+FAfYoJk6QUhbOHKLmGE+BqIDmzm/6eGTLTchFPRKKqtLa+OoqIi/P39KSwsxM/PT+1whBBCCNEITfn9lj43QgghXFPmBjCWqR2FaAaS3AghhLBLiqJwMLeE9LxS2x+8ohA+Hw5vxEKpSquQi2YjyY0QF2nNwRP0fPlP7vlqo9qhCOFU3lt8gGHvLGfasoO2P3j6SlBMEBAN3sG2P75QlSQ3QlykYwXlFJRVUWqsVjsUIZxKz2h/ANYfyrf9wa1DwGWUlDOS0VJCXKSru0XQNcofs2v1zRei2fVpE4RGA+l5peQUVRDqZ8N5pGon72sr89s4I6m5EeIieep1dAj3pVOEjL4Twpb8Pd3pXPP/1bp0G9benDwE+Wmg0UHsQNsdV9gNSW6EEELYrf5xQQCst2Vyc7B2yYX+4CH/KHFGktwIcZH+u+QAU5emklNUoXYoQjidxLhWAKxLP2G7g0qTlNOTPjdCXKRPV6ZRXFHNlZ3DbNsnQAhBv9hAAPZnl5BfaiTIW3/xBx36rGWhTOlM7LSk5kaIi1BSWU1xhWWUVISsKyWEzbXyMZAQ6gPABluNmgrpAEkPQGgn2xxP2B1JboS4CMdrViz29XDDxyAVoUI0h8S2ln4369KaYUi4cEqS3AhxEY4VWvrZyGrgQjSf/jX9btYfskG/m+VvwpZvoLzg4o8l7Jbqyc3UqVOJjY3Fw8ODxMRE1q9ff87yBQUFPPDAA0RERGAwGGjfvj3z589voWiFqKu25ibCX/raCNFcEmtGTO0+VkRRRdWFH6iyxJLc/PoAlNmwg7KwOxeU3BQUFDB9+nQmTZpEfr6lmnDz5s0cPXq0SceZNWsWEyZMYPLkyWzevJkePXowfPhwcnJy6i1vNBq54oorOHToELNnz2bfvn18+umnREVFXchlCHHRamtupL+NEM0nzM+D2FZemBXYdOjkhR/o8N9groKANhDU1nYBCrvT5E4C27dvJzk5GX9/fw4dOsTdd99NUFAQP//8MxkZGXz11VeNPta7777L3Xffzfjx4wGYNm0a8+bN4/PPP2fixIlnlf/888/Jz89n9erVuLu7AxAbG3vOc1RWVlJZWWm9X1RU1Oj4hDif2pqbSKm5EaJZJca1QqfVUFltuvCDWJdcGAoajW0CE3apyTU3EyZMYNy4cRw4cAAPj1Nf6FdffTUrVqxo9HGMRiObNm0iOTn5VDBaLcnJyaxZs6be58ydO5ekpCQeeOABwsLC6Nq1K6+99homU8Mf9pSUFPz9/a236OjoRscoxPkcK6xtlpKaGyGa06vXd2XJ40O4qmvEhR+kdvI+GQLu9Jqc3GzYsIF77733rP1RUVFkZWU1+jh5eXmYTCbCwsLq7A8LC2vwOGlpacyePRuTycT8+fN5/vnneeedd/jPf/7T4HkmTZpEYWGh9ZaZmdnoGIU4n+MFlmapSGmWEqJZuekusoto4VHI2wcaLcRdZpughN1qcrOUwWCot2ln//79hISE2CSohpjNZkJDQ/nf//6HTqejT58+HD16lLfeeovJkyc3GK/BYGjWuIRrUhTFWnMTGSDNUkK0hCqTmcpqc9OnXqidlTiyF3gG2j4wYVeanApfe+21vPzyy1RVWXqsazQaMjIyePrpp7nhhhsafZzg4GB0Oh3Z2dl19mdnZxMeHl7vcyIiImjfvj06nc66r1OnTmRlZWE0Gpt6KUJclIKyKiqqzACES58bIZrdR8tS6f7in/xv+cGmP/nEQUAjTVIuosnJzTvvvENJSQmhoaGUl5czePBg4uPj8fX15dVXX230cfR6PX369GHJkiXWfWazmSVLlpCUlFTvcy699FJSU1Mxm83Wffv37yciIgK93gZTcgvRBLW1NsE+egxuuvOUFkJcrEAvPeVVJnYeu4CBIcmT4cmD0P/sbhXC+TS5Wcrf359Fixbx999/s23bNkpKSujdu3edjsGNNWHCBMaOHUvfvn3p378/U6ZMobS01Dp66vbbbycqKoqUlBQA7r//fj788EMeeeQRHnroIQ4cOMBrr73Gww8/3ORzC3GxNGgYlBCMn4e72qEI4RKGdwmnV0wA7UN9L+wA3q1sG5CwWxc8X/yll17KpZdeelEnHz16NLm5ubzwwgtkZWXRs2dPFi5caO1knJGRgVZ7qnIpOjqaP/74g8cee4zu3bsTFRXFI488wtNPP31RcQhxITpH+vH1nYlqhyGEywjy1l/YwpmKIkO/XYxGURSlKU94+OGHiY+PP6u25MMPPyQ1NZUpU6bYMj6bKyoqwt/fn8LCQvz8/NQORwghRHP79l9gLIMrX4aoPmpHIy5QU36/m9zn5qeffqq3xmbAgAHMnj27qYcTwmFVm8znLySEsKn92cU8OnMLE37Y2rgnVJVD2nI4vAr0Ps0am7AfTU5uTpw4gb+//1n7/fz8yMvLs0lQQjiCmz9dR8+X/2TZvvqXCxFC2J5ZUZiz9RgLd2ZR1Zh/YGSsAVMl+EVBcPvmD1DYhSYnN/Hx8SxcuPCs/QsWLKBtW1mrQ7iOowXlFJRV4etxwV3XhBBN1D7UF39Pd8qMJnY1ZtTUoVWWv3GXSb8bF9Lkb+UJEybw4IMPkpuby+WXW+YLWLJkCe+8847d97cRwpYWPDqI4wUVxAR5qR2KEC5Dq9XQLzaIxXuyWZd2gp7RAed+wqG/LX/bXNwAGOFYmpzc3HHHHVRWVvLqq6/yyiuvAJbFKz/++GNuv/12mwcohL3y83DHL1yGgQvR0i5pa0lu1qfnc+/gdg0XrCqHo5ss27GS3LiSC6pPv//++7n//vvJzc3F09MTHx/ppCWEEKJl9I8LAmD9oXxMZgWdtoHmpiMbwFwFvhEQGNeCEQq1XVRngeZeS0oIe7XxUD6/bDlKnzaB/LN3a7XDEcKldI7ww8fgRnFFNXuziugSefYgFwB0Bmg/Avwipb+Ni2lyh+Ls7Gxuu+02IiMjcXNzQ6fT1bkJ4Qq2Zhbw7boM/torI6WEaGluOi192lgWv1yfnt9wwZhEuHkm/OPdFopM2Ism19yMGzeOjIwMnn/+eSIiItBINixc0LGCCgAiAzxVjkQI19Q/Lojl+3NZn57P+EulyUnU1eTkZtWqVaxcuZKePXs2QzhCOIbjNYtmRshq4EKoIrG23016PoqinP0P7ZIcqK6AgBgVohNqa3KzVHR0NE1csUEIp3Os0FJzE+EvNTdCqKF76wAMblpOlBo5mFtydoHNX8KUbvDbIy0fnFBdk5ObKVOmMHHiRA4dOtQM4QjhGI4XWGpuIgOk5kYINejdtPSOsfS7WZtWT7+bw6stf0M7t2BUwl40uVlq9OjRlJWV0a5dO7y8vHB3rzvPR37+OTp3CeEEjNVmcksqAelzI4SaEtsGsSbtBOvT87n1kjanHjBVQcY6y7ZM3ueSmpzcyCzEwtVlF1WgKJZ/Obby1qsdjhAu6x/dI4gL9iYxrlXdB45vg6pS8AyUmhsX1eTkZuzYsc0RhxAO41jBqc7EMlpQCPXEh/oSH+p79gO160nFDABtk3tfCCdwUZP4VVRUYDQa6+zz8/O7qICEsHfHrZ2Jpb+NEHaptr+NLLngspqc0paWlvLggw8SGhqKt7c3gYGBdW5COLtjhbWdiaW/jRBqO3KyjGnLDzJ9ZZplh9kEGWss220GqBeYUFWTk5unnnqKv/76i48//hiDwcD06dN56aWXiIyM5KuvvmqOGIWwK8drJ/CTYeBCqC4tt5TXF+zli78PWXYoZrj2A7jk/yC8u6qxCfU0uVnqt99+46uvvmLIkCGMHz+eQYMGER8fT5s2bfj222+55ZZbmiNOIeyGtc+NDAMXQnV92gRyZecw+scFWRbR1LlDl+ssN+Gympzc5Ofn07ZtW8DSv6Z26PfAgQO5//77bRudEHaoR3QAVWaFhPo6MgohWpS3wY3/3d5X7TCEnWlyctO2bVvS09OJiYmhY8eO/PDDD/Tv35/ffvuNgICAZghRCPvy8LAEtUMQQtTHbIa/p0B0IsRcAlpZzNlVNbnPzfjx49m2bRsAEydOZOrUqXh4ePDYY4/x5JNP2jxAIYQQ4lwURSEzv4zlf6+AJS/Bt/+y9L0RLqvJNTePPfaYdTs5OZm9e/eyadMm4uPj6d5dOm8J52asNlNlMuNtuKhZFIQQNlRcWc1lby3lVu2fDHYHovuDzv28zxPO66K/odu0aUObNm3OX1AIJ7A+PZ9bP1tH75gAfv4/mUNDCHvg5+FO5wg/EnP3WHbIkgsur1HJzX//+99GH/Dhhx++4GCEsHdZRZZh4D4e8q9CIexJYmwQiSdqkhuZvM/laRRFUc5XKC4urs793NxcysrKrB2ICwoK8PLyIjQ0lLS0tGYJ1FaKiorw9/ensLBQZlMWF6Sooopyo4kwPxkKLoS9WLlmNYP+GIERd/TPHQU3g9ohCRtryu93ozoUp6enW2+vvvoqPXv2ZM+ePeTn55Ofn8+ePXvo3bs3r7zyik0uQAh75ufhLomNEHamt7IbgM3meE5UqByMUF2TR0s9//zzfPDBB3To0MG6r0OHDrz33ns899xzNg1OCCGEaAzvnM0ArDN3Yn16vsrRCLU1Obk5fvw41dXVZ+03mUxkZ2fbJCgh7NUzv+zg2V92cORkmdqhCCFON/K/TGs/nR9Ng1knyY3La3JyM2zYMO699142b95s3bdp0ybuv/9+kpOTbRqcEPZEURTmbDnKt+syMFbLHBpC2BWdG627XsoRJUSSG9H05Obzzz8nPDycvn37YjAYMBgM9O/fn7CwMKZPn94cMQphF4rKqykzmgCIkEUzhbA7/eOCANibVURhWZXK0Qg1NXmem5CQEObPn8/+/fvZu3cvAB07dqR9+/Y2D04Ie3Ks0LJgZqCXO556mdZdCLux5BUoySa03120DfEmLbeUDYfySe4cpnZkQiUXPIlf+/btJaERLuV4TXITGSC1NkLYlR0/QsFh6HwdiXGhpOWWsi79hCQ3LqzJyY3JZGLGjBksWbKEnJwczOa6fQ/++usvmwUnhD05WmAZXypNUkLYkcIjlsRGo4OYRBJLivh+fab0u3FxTU5uHnnkEWbMmME111xD165d0Wg0zRGXEHbneEFtzY3McSOE3Ti82vI3ogcYfEls60ZMkBedI/xQFEV+o1xUk5ObmTNn8sMPP3D11Vc3RzxC2K3jhVJzI4TdObTK8rfNAMDy/+eKp4aqGJCwB00eLaXX64mPj2+OWISwa8ek5kYI+1NbcxM7UN04hF1pcnLz+OOP8/7779OIJamEcCpScyOEnSnOhhMHAA3EXFLnIZNZITWnRJ24hOqa3Cy1atUqli5dyoIFC+jSpQvu7nVXR/75559tFpwQ9sJsVsiyJjdScyOEXSg+BsHtLYtkegZad+eVVDL0rWWUV5nY8eJwmbrBBTU5uQkICOD6669vjliEsFsnSo0YTWY0GgiX5EYI+xDZCx7cAFXldXa38tbj42H5eUvLK6FLpL8a0QkVNTm5+eKLL5ojDiHsWu0cN6G+Btx1TW7NFUI0J/e6TcUajYYf70siwt8TnVZGS7miC57ETwhX4qV349/9oqV6Wwh7UW0EjQZ07vU+3DrQq4UDEvakyclNXFzcOecNSEtLu6iAhLBH8aE+vH5Dd7XDEELUOvAn/Hw3dLsRrv3gnEVlvhvX0+Tk5tFHH61zv6qqii1btrBw4UKefPJJW8UlhBBCNOzwaqgqs8xM3ICU+XtYsDOLt27sTmLbVi0YnFDbBc1QXJ+pU6eycePGiw5ICHuUXVSBj8ENb4O05AphFw7XTN53jvltjhSUk5Ffxvr0fEluXIzNekaOGDGCn376yVaHE8Ku3PfNJrpM/oM/d2WpHYoQoqIQsnZYtmOSGix2SVwQgKwz5YJsltzMnj2boKAgWx1OCLtSVF4FyDBwIexC5npQzBAYC/5RDRbrH2eprdl0+CRVJnOD5YTzaXIde69evep0zFIUhaysLHJzc/noo49sGpwQ9mLJ40MoqqjC011GSwmhOut6UudeciEh1IdAL3dOllWx42ghvWMCz1leOI8mJzfXXXddnftarZaQkBCGDBlCx44dbRWXEHbHz6P+IadCiBZWu55UzWKZDdFqNfSPC+KPXdmsS8uX5MaFNDm5mTx5cnPEIYQQQjROp5Fg8IXYS89btH9cK/7Ylc3atBPcP6RdCwQn7MEFDf0wmUz88ssv7NmzB4DOnTszatQo3NxkJIlwPn/tzebrNYcZlBDCHQPj1A5HCHHpw5ZbIwxoZ+l3sz49n8pqEwY3aVp2BU3ORnbt2sXIkSPJzs6mQ4cOALzxxhuEhITw22+/0bVrV5sHKYSa9hwvZum+XFr5GNQORQjRRB3DfQn2MZBXUsnmwwUktZMh4a6gyaOl7rrrLrp27cqRI0fYvHkzmzdvJjMzk+7du3PPPfc0R4xCqOpogWVdqcgAz/OUFEI0u7RlUNz4KRk0Gg2DEoIBWHkgt5mCEvamycnN1q1bSUlJITDwVMeswMBAXn31VbZs2WLT4ISwB8dqkpuoABkGLoSqqivhu9HwTgc4cbDRTxsYb0luVqXmNVdkws40Oblp37492dnZZ+3PyckhPj7eJkEJYU9qk5sIf6m5EUJVx7ZAdQV4h0BQ20Y/bWBNzc2Oo4WcLDU2V3TCjjQquSkqKrLeUlJSePjhh5k9ezZHjhzhyJEjzJ49m0cffZQ33nijueMVosUdK6gApFlKCNVZ57cZYFkRvJHC/DzoEOaLosDGwyebKThhTxrVoTggIOCsiftuuukm6z5FUQAYOXIkJpOpyUFMnTqVt956i6ysLHr06MEHH3xA//79z/u8mTNnMmbMGEaNGsWcOXOafF4hzqeoooqSymoAIqVZSgh1Wee3Of8Q8DO99s9uhPoaiA7ysnFQwh41KrlZunRpow62Y8eOJgcwa9YsJkyYwLRp00hMTGTKlCkMHz6cffv2ERoa2uDzDh06xBNPPMGgQYOafE4hGqu2SSrQyx0vvUx1IIRqTNWQuc6yfZ7J++rTp41M4OdKGvVtPXjw4AYfKy4u5vvvv2f69Ols2rSJBx98sEkBvPvuu9x9992MHz8egGnTpjFv3jw+//xzJk6cWO9zTCYTt9xyCy+99BIrV66koKCgweNXVlZSWVlpvV9UVNSk+IRrOyYjpYSwD1nbwFgCHv4Q2kXtaISdu+CFM1esWMHYsWOJiIjg7bff5vLLL2ft2rVNOobRaGTTpk0kJyefCkirJTk5mTVr1jT4vJdffpnQ0FDuvPPO854jJSUFf39/6y06OrpJMQrXdlT62whhH2qbpGIGgPbCfrr+3JXFXV9u5IeNmTYMTNijJtWzZ2VlMWPGDD777DOKioq46aabqKysZM6cOXTu3LnJJ8/Ly8NkMhEWFlZnf1hYGHv37q33OatWreKzzz5j69atjTrHpEmTmDBhgvV+UVGRJDii0aw1N7IauBDq6nYT+ISD94VPwncgp4TFe7LRaOCmvvI74MwandyMHDmSFStWcM011zBlyhSuuuoqdDod06ZNa8746iguLua2227j008/JTg4uFHPMRgMGAwys6y4MNIsJYSd8A2D7v+6qEMM7xKGRgND2jfcn1M4h0YnNwsWLODhhx/m/vvvJyEhwSYnDw4ORqfTnTVvTnZ2NuHh4WeVP3jwIIcOHWLkyJHWfWazGQA3Nzf27dtHu3ayMJqwnePSLCWE04gP9SU+1FftMEQLaHTD5apVqyguLqZPnz4kJiby4Ycfkpd3cbM96vV6+vTpw5IlS6z7zGYzS5YsISkp6azyHTt2ZMeOHWzdutV6u/baaxk6dChbt26V5iZhc//oEcGY/tF0DJcvRCFUs28hrHoPcvaoHYlwEI2uubnkkku45JJLmDJlCrNmzeLzzz9nwoQJmM1mFi1aRHR0NL6+Tf8BmDBhAmPHjqVv377079+fKVOmUFpaah09dfvttxMVFUVKSgoeHh5nLcwZEBAAIAt2imZxe1Ks2iEIIbZ9B7t/BbMJQjtd1KFKK6tZvCebjBNlPDTMNq0Qwv40ucu5t7c3d9xxB6tWrWLHjh08/vjjvP7664SGhnLttdc2OYDRo0fz9ttv88ILL9CzZ0+2bt3KwoULrZ2MMzIyOH78eJOPK4QQwgkoykVN3nemgvIqHpm5lSlLDlBUUXXRxxP2SaPUTi98EUwmE7/99huff/45c+fOtUVczaaoqAh/f38KCwvx8/NTOxxhxwrLq8gtriDC3xNvg0zgJ4QqcvfD1H7g5gETM8Dt4geIXP72MtLySvnktj4M73J2/05hn5ry+33B89ycTqfTcd1119l9YiNEU6zYn0vyuysY/8UGtUMRwnUd/tvyt3U/myQ2AINqFtJcdUBWCXdWNkluhHBG5UYTvh5usqaUEGqyNkk1fcmFhgxMCAFg5YFcmx1T2BepaxeiATf1i+amftFUmcxqhyKEa1KUUzU3NuhvU+uStkHotBoOnSgjM79MFtN0QlJzI8R5uOvkfxMhVFGSA2X5oHWzNEvZiK+HO71jAgBYKU1TTkm+tYUQQtgn3zBLJ+J7loPetrUrA+MtTVOrUqVpyhlJciNEA/41bTXjv1hPdlGF2qEI4brc9BBu+3nMBtZ0Kv479QQm80UPGhZ2RvrcCFGPMmM1Gw6dBMBTr1M5GiGErfVo7Y+vhxuF5VXsOFpIz+gAtUMSNiQ1N0LU41jNmlK+Bjf8PNxVjkYIF1R0DKYmwvwnLR2LbcxNp2VAO8sK46tk1JTTkeRGiHrIauBCqOzwasjdC5nrQaNpllOcGhIunYqdjSQ3QtSjNrmJkDluhFBHMwwBP9NlNf1uNmecpLSyutnOI1qe9LkRoh7HCi3NUlJzI4RK0ldY/sY2X3LTppU3Y5Pa0CXSv7kqh4RKJLkRoh61NTdRktwI0fIKMuFEKmh0EDuwWU/10ijbj8QS6pNmKSHqcarPjTRLCdHi0pZa/kb1AQ9/dWMRDkmSGyHqYU1u/KXmRogWd7AmuWk3tEVOdyC7mM9WpZMjc1o5DWmWEuIMiqJInxsh1BQQY7m1bZnk5onZ29mWWYCvhxs39Y1ukXOK5iXJjRBnOFFqxFhtRqOBcH9plhKixV3xEiS/2HKn6xSKn4cbwT76FjunaF6S3AhxhtomqVBfgyyaKYRaWnD40oOXJ/Bgi51NtAT55hbiDLWzE0uTlBAqyN0HJplzRlwcqbkR4gzxod48fkV7gqSKWoiWVVUO0waBmwc8sA78Ilr09NlFFRzKKyWxbasWPa+wPUluhDhDfKgvDw3zVTsMIVxPxhowVYJXK/ANb9FTr0/P56ZP1hDmZ2DNxGFotTKrnyOTZikhhBD24fQh4C08ZXD31v5463VkF1WyJbOgRc8tbE+SGyHOsDnjJKk5xRirzWqHIoRrqZ28r4WGgJ/Ow13HsE5hACzcebzFzy9sS5IbIc5w79ebSH53Bfuzi9UORQjXUZILWTss222HqBLCiK6WprAFO7NQFEWVGIRtSHIjxGlMZoVW3np8PdyIkDluhGg56cstf8O6gU+IKiEM7hCCh7uWIyfL2XWsSJUYhG1Ih2IhTqPTalj46GVqhyGE67H2txmiWgheejeGtA9l4a4sFuw8TtcoWdfKUUnNjRBCCPX1uxMGT4TO16kaxohu0jTlDKTmRgghhPqieltuKru8Yyh6nZa03FIO5JTQPkymhXBEUnMjxGmmr0xj+Hsr+GxVutqhCCFU4OvhzqCEYAAW7MhSORpxoSS5EeI0B3NL2JddTHFFldqhCOE6NkyH3XOh0j5GKF5lHTUlQ8IdlSQ3QpzmqKwrJUTLMlXBohfhh9sg74Da0QBwRecw3LQa9mYVk55XqnY44gJIciPEaWpXBI/0l+RGiBZxdBMYi8EzECJ6qB0NAAFeepLaWdaXktobxyQdioWooSgKx2uTmwCZ40aIFlE7BDxuMGh16sZymht6tybE10CfmEC1QxEXQJIbIWoUlVdTajQB0iwlRItJW2b5267ll1w4l+t6RXFdryi1wxAXSJqlhKhxtKbWppW3Hg93+/kXpBBOq6IIjmywbKuwnpRwXpLcCFHD2t9Gam2EaBmHVoFigqC2ENhG7WjOoigKO48W8s3aw2qHIppImqWEqHG8UPrbCNGijm2x/LXTWpuc4kr+8cEqNBq4snMYoX7y3eAoJLkRokbtMPAIGSklRMu4/FnodYvaUTQozM+DS+Nb4Wtwp6SymlC1AxKNJsmNEDVqm6WipFlKiJYTGKt2BOf0zZ2JaDQatcMQTSR9boSoIX1uhBBnksTGMUlyI0SN44W1sxNLu7oQze63R+G7f0PmBrUjaZS03BJ2HStUOwzRSJLcCFHjyeEdePyK9sQFe6sdihDOzWyGPb/B/gVgMqodzXl9vfYwl7+znDcX7lM7FNFI0udGiBoyYZcQLSR7J5Tlgbs3tO6ndjTnldTWshTD6oN5FJZX4e/prnJE4nyk5kYIIUTLSqtZciF2ILjp1Y2lEeJDfUgI9aHKpLBkT7ba4YhGkORGCODwiVKW788lM79M7VCEcH6160m1HaJqGE0xoms4AAt2ZqkciWgMSW6EAObvyGLs5+t5d9F+tUMRwrlVFsPh1ZZtO1tP6lxGdIsAYPn+XEoqq1WORpyPJDdCAD4GHR3DfYkP9VE7FCGc2/4/wFQJQe0gpKPa0TRax3Bf2oZ4Y6w28+vWo2qHI85DkhshgNuSYln46GU8MDRe7VCEcG4e/hAzALr+ExxoDhmNRsMtiZb1r75ecxhFUVSOSJyLJDdCCCFaTsIVcMcCGPqs2pE02Y29W+PhrmVvVjGbDp9UOxxxDpLcCCGEaHkOVGtTy9/LnVE9LFNGfC0rhds1SW6EyyutrKbb5D8Y/t4KKqpMaocjhPM6uBRKT6gdxUW5LcnSNDV/x3HySipVjkY0RJIb4fKOF5ZTXFnN8cJyPNx1aocjhHMylsLMm+HtBMhLVTuaC9Y1yp+e0QFUmRRmbchUOxzRAEluhMs7WlC7ppQsmClEszmwCKrKwL81tGqndjQX5bZLLLU3363LwGyWjsX2SJZfEC5PVgMXogXs/tXyt/Moh+xvc7prukew42gh/+rbGq3Wsa/FWUlyI1zeqeRGVgMXollUlVvmtwHocp2qodiCh7uOF6/tonYY4hykWUq4vKM1yU2Ev9TcCNEsDiyCqlLwj4bI3mpHI1yA1NwIl1dbcxMlzVJCNA8napI63YHsYv63Io1wfw8ev7KD2uGI09hFzc3UqVOJjY3Fw8ODxMRE1q9f32DZTz/9lEGDBhEYGEhgYCDJycnnLC/EuSiKwv7sEgBig71VjkYIJ2SqgtRFlu3O16kaiq1l5Jfx46YjfL32sEwjYWdUT25mzZrFhAkTmDx5Mps3b6ZHjx4MHz6cnJycessvW7aMMWPGsHTpUtasWUN0dDRXXnklR4/KWh+i6Y4XVpBfakSn1dAx3FftcIRwPjp3eGgzjPoIovqoHY1NDekQyrgBsXw2th8GN9V/TsVpNIrKC2QkJibSr18/PvzwQwDMZjPR0dE89NBDTJw48bzPN5lMBAYG8uGHH3L77beft3xRURH+/v4UFhbi5+d30fELx/bHrizu/XoTHcN9WfjoZWqHI4QQogFN+f1WNdU0Go1s2rSJ5ORk6z6tVktycjJr1qxp1DHKysqoqqoiKCio3scrKyspKiqqcxOi1q6jhQB0i/JXORIhhBC2ompyk5eXh8lkIiwsrM7+sLAwsrKyGnWMp59+msjIyDoJ0ulSUlLw9/e33qKjoy86buE8dtQkN10luRHC9vb/CZ+PgG0z1Y6kWWWcKOO5OTt4ce4utUMRNRy6kfD1119n5syZ/PLLL3h41D9HyaRJkygsLLTeMjNlumxxys5jlpo8SW6EaAa7foaM1XBsi9qRNKvckgq+WZvB9+szKCyrUjscgcpDwYODg9HpdGRnZ9fZn52dTXh4+Dmf+/bbb/P666+zePFiunfv3mA5g8GAwWCwSbzCuSiKwhfj+rHrWCGdI6T/lRA2VW2EvfMt251HqRtLM+sdE0inCD/2HC/ix02Z3DWordohuTxVa270ej19+vRhyZIl1n1ms5klS5aQlJTU4PPefPNNXnnlFRYuXEjfvn1bIlThhDQaDV2j/BndLwZPvSyYKYRNpS2DykLwCYfoS9SOpllpNBrrelPfynpTdkH1ZqkJEybw6aef8uWXX7Jnzx7uv/9+SktLGT9+PAC33347kyZNspZ/4403eP755/n888+JjY0lKyuLrKwsSkpK1LoEIYQQZ6qduK/TSNCq/lPT7Eb1jMTX4EZ6Xil/H8xTOxyXp/onbvTo0bz99tu88MIL9OzZk61bt7Jw4UJrJ+OMjAyOHz9uLf/xxx9jNBq58cYbiYiIsN7efvtttS5BOKgvVx/i6zWHyC6qUDsUIZyLqQr2/m7ZdoK1pBrD2+DGDX1aA/D1msMqRyNUn+empck8N6LWJa8tIauogh/vS6JfbP1TCQghLkDqYvjmBvAOgcf3gdY1mn1Tc4pJfncFWg2sfPpyWdLFxhxmnhsh1FJtMvOvvq0Z0iFEOhMLYWt6H0gYDl1vdJnEBiA+1JcB7VphVuDDvw6oHY5Lk5obIYQQzUNRnGqhzMbYeCifG6etQauBhY9eRvswWdbFVqTmRgghhPpcLLEB6BsbxFVdwjEr8PqCvWqH47IkuREuaefRQnKLK9UOQwjnc/AvOOnaHWqfHtERN62Gv/bmsDpVRk6pQZIb4ZLu/XoT/V5dzIZD+WqHIoTzMFXDT3fD+90hc73a0agmLtibWxJjAHh1/h6Z90YFqs5QLIQa8kuNHC0oB6BDuLSHC2Ez++ZBWR54tYLIXmpHo6qHhyWw+3gRD16egFbres1zapPkRricnTWLZcYFe+Pn4a5yNEI4kTVTLX/73gE61/5/q5WPgR/vG6B2GC5LmqWEy9l5zJLcdImU0XJC2EzmBshcBzo99Ltb7WjsTpXJrHYILkWSG+FyamtuuslK4ELYzpoPLX+73QS+YerGYkeqTGamLk1l8JtLyS81qh2Oy5DkRricnUeLAOgqyY0QtnHyMOyZa9lO+j91Y7EzOo2GeduPc6ywgpkbMtQOx2VInxvhUgrLqsjILwOga6QkN0LYRPYuy6zErftCWBe1o7ErWq2GF6/twtGCMkb1iFI7HJchyY1wKbtq+ttEB3ni7+XaHR6FsJmOV8OE3VB2Qu1I7FL/uCBA1q9rSdIsJVzKDulvI0TzMPhCYKzaUdi9MmM1mTW1x6L5SHIjXMrOY5b+Nl2kSUqIi2eqgkOrLGtIifPacCifIW8t45GZW3CxZR1bnCQ3wqXISCkhbGj3rzDjGvj2RrUjcQhtgrworqhmc0YBC3ZmqR2OU5PkRriM4ooq0vNKARkpJcRFUxRY/YFlOzpR3VgcRKifB/dc1haANxbuxVgtc980F0luhMs4WlBOsI+eqABPgrz1aocjhGM7vBqObwU3D+h7p9rROIx7LmtLiK+BwyfKmL4qTe1wnJYkN8JldAz3Y8Ozycx/eJDaoQjh+Gon7esxBrxbqRuLA/E2uPHk8A4AvPPnfjbK4r3NQpIb4VI0Go0MARfiYuWlwr4Flu1LZNK+pvpXn9aM7BGJyazw4HdbOFFSqXZITkeSGyGEEE2z9iNAgfZXQUh7taNxOBqNhpR/dqNtsDdZRRU89sM2zGYZPWVLktwIl1BaWc2AlCXc9eUGKqpMaocjhOMym+HIBst20gPqxuLAfAxufHRrbzzctazYn8tHy1LVDsmpSHIjXMKe40UcK6xgx9FCPNx1aocjhOPSauGe5XD7XIiV/msXo2O4H6+M6grAu4v2s/pgnsoROQ9JboRL6Brlz4/3JfGf67qpHYoQjk+rhbaDQaNROxKH96++0fyrT2vMCjz8/VZyiivUDskpSHIjXIKHu45+sUFc0TlM7VCEcFx5qVBVrnYUTuflUV3pEOZLXkklszcdUTscpyALZwohhDg/RYEfboeSbPj3txBzidoROQ1PvY6Pbu3N2rQT3Nw/Ru1wnIIkN8LplRtNvLFwL12j/Lm+VxQ6rVSlC9Fk++ZDzi5w94aQDmpH43TahfjQLsRH7TCchjRLCae3J6uIGasP8fqCvUheI8QFKC+AeY9btvvfDZ6Bqobj7ArLqnjix20cL5QmwAslyY1werWLZXaN8kMjHSCFaLo/n4Xi4xDUDoZMVDsap/f0T9uZvekIj/+wTVYPv0DSLCWcnqwELsRFSF0MW74BNDBqKrh7qh2R05t0dUcyT5YxaUQn+QfZBZLkRji9HUeLAOgSKcmNEE1SUQRzH7FsJ94LbZLUjcdFtGnlze8PDZTE5iJIs5RwahVVJg5kFwPQrbUkN0I0WbshEBgLw15QOxKXcnpi83dqHhNmbcVYbVYxIsciNTfCqe3LKqbarBDo5U6kv4fa4QjhWDz8LE1RFYWg91Y7GpdUXFHF/327mcLyKgrKq/jolt4yy3ojSM2NcGo7j9V2JvaXKl4hGqvaaJnXppaH1HqqxdfDnff/3RMPdy1/7c1h3BfrKamsVjssuyfJjXBqp0ZKyZezEI3257Pw7b+gUGbLtQdDOoTy1R2J+BjcWJuWzy3T11FQZlQ7LLsmyY1wajtrOhPLSCkhGunQ37D+f5C6CPL2qx2NqNE/Lojv7k4kwMudbZkF/Pt/a2UdqnOQ5EY4razCCvZmSXIjRKMZy+DXByzbvW+HdperG4+oo3vrAGbdk0SIr4G9WcWM/mQtRwtkor/6SHIjnNaHSw9QZVLoHxtEdJCX2uEIYf/++g+cTAe/KLjyP2pHI+rRIdyXH+9NIirAk/S8Uv718WrS80rVDsvuSHIjnFJ2UQWzNmQC8PiV7VWORggHkLEO1n5k2R75vnQitmOxwd7Mvj+JtiHeHCusYNSHq/h58xGZzfg0ktwIpxTqa+CLcf25a2AciW1bqR2OEPatqrymOUqBHjdDwhVqRyTOI8Lfkx/uTaJHdABFFdVM+GEb932zCbNZEhyQ5EY4KY1Gw8CEYJ77R2e1QxHC/lVXgpsH+ITDVa+pHY1opGAfAz/dl8STwzvgrtMQ28obrawODIBGcbF6rKKiIvz9/SksLMTPz0/tcEQzKDea8NTLJFdCNElJLuTtg9iBakciLsC+rGLatPKyTvCXmV+Gp15HsI9B5chspym/31JzI5zKvqxi+r+2mPcW7Zf2ZyHOxVgKe347dd8nRBIbB9Yh3Nea2JjMCg99v4Ur31vB2rQTKkemDkluhFOZvSmT4opqUnNKZEZiIRpSUQTf3ACzbq1Z8Vs4kxMllVRUmaiqNtOmlWuOFJW1pYRTmTSiE71iAmkf5qt2KELYp7J8+OafcGwLGPwhWEYTOptQPw/mPjiQfVnFRPh7AqAoCt+vz+SabhH4e7mrHGHzkz43QgjhKkpy4KvrIGcXeLWC236BiB5qRyVawOrUPG6evg5vvY5bk9pw58A4Qn0dazFh6XMjXM6Rk2WUGWUxOSEaVHgEvhhhSWx8wmHcfElsXIkGOob7Umo08cnyNAa+sZTn5+wkM79M7ciahdTcCKdw86dr2Z9dzJTRvRiYEKx2OELYl/IC+GQQFGSAfzTc/iu0aqd2VKKFKYrCX3tz+HBpKlsyCgDQaTWM6hnJ/w1pR3yofTfnN+X3W/rcCIe3OjWP1QdPoNdpiQvxVjscIeyPZwB0Hw07f4Lb50JAtNoRCRVoNBqGdQrj8o6hrE3L56Nlqaw8kMfPm4/yy5ajXNEpjBv7tOay9iHWkVeOSmpuhENTFIUbp61h0+GTjE1qw0ujuqodkhD2oSwfio9DWBfLfUWByiJZVkHUsS2zgI+WpfLHrmzrPl+DG1d0DuOxK9rb1bp80udGuIzl+3PZdPgkBjctDwyNVzscIdRXbYQ1U+G/PWH2HWCq6Yum0UhiI87SIzqAT27ry6LHLuPOgXGE+3lQXFnNnK1HMbifShGyCiuoMplVjLRppFlKOCxFUXjnz/0A3J7UhlA/x+r5L4RNKQrsWwB/Pgf5By37/KOhJAv8W6sbm7B7CWG+PP+Pzjx7dSc2Z5xkx9HCOqOpHp21hX1ZxUz5dy8Gtw9RMdLGkeRGOKw/d2ez42ghXnod9w2WzpHChWXtgD+egfQVlvveIXD589DrVtA6dt8J0bK0Wg19Y4PoGxtk3VduNJGWW8rJsiraBp/q1/jz5iOsT8+nf1wQ/eOCaB1oP01YktwIh1RUUcW7NbU2d1waRysnWj9FiCY5vg0+GQwooDNA0v/BwAngIX0KhW146nWsnng5O44W1umDs3BnFn/uzmbmhkwAIv096B8XxIB2wdzUT91O65LcCIdirDbzzdrDfPDXAU6WVeHr4cbdg9qqHZYQLefkYTi+FTqPstwP725ZE8o7GJJfgsA2qoYnnJObTkuvmMA6+25PiiUuxJv16fnsOFLIscIK5mw9RvqJMkluhGgMRVGYvyOLN//Yy+ETlkmn4kN9eO36bi4xlbhwYYoCufssi1zu/c1SU6N1hycHW4Z4azRwy2xwlz5nomUNTAi2zitWZqxmS0YB69LzCbeD/o+S3AiHYFbgw6WpHD5RRoivgQlXtOdffVrjppMBf8JJ5afB5q8tSc2JA6f2a7QQcwmU5lqSG5DERqjOS+/GpfHBXBpvH5Oo2sUvw9SpU4mNjcXDw4PExETWr19/zvI//vgjHTt2xMPDg27dujF//vwWilS0pIO5JdYlFXRaDc9f04lHkxNY9sQQxvSPkcRGOL6qCktNzLaZsOgFyNp56rGMtbDqXUtio9NDwpVw7QfwxAEY9zsEJ6gXtxB2TvWam1mzZjFhwgSmTZtGYmIiU6ZMYfjw4ezbt4/Q0NCzyq9evZoxY8aQkpLCP/7xD7777juuu+46Nm/eTNeuMoGbo6qdP8G9JmH5YMkBpiw5wGPJCTx4ueVLfEB8MAPs5F8FQpyTooCxxLJQpUcAeLey7M87AGs/guJsyNtvGbKtnDZ3iF9rCK/5HovoYelX0+laSLhC5qgRoglUn6E4MTGRfv368eGHHwJgNpuJjo7moYceYuLEiWeVHz16NKWlpfz+++/WfZdccgk9e/Zk2rRp5z1fc85QbKys4ETW4QYfN+t9UTwCLHdMVehKswAI9jFYf9SLKqoprazCyycA/1ZhlqLV1eQcPdjgcRV3L8yeNV+eihld8VEAAr0NeLhZjltSWU1xRRUGTx+CQqNqiprJyjxAQ58Axc0Ts9epZEJXlImiQJCPHs+aqblLKqvJLzGi9/AiPCrGWnbv3j0Yq6uoMilUmRWqq81UmcxUKwqlVVrSKnw4crKcIyfLMZ7MJK+ojNf/2d3afvvn7ixemruby7tE8cptV5wKqugYmBtYIFPrBn6Rp+4XZ4HJWH9Zjbbu3B/F2WCqrL8smrrT1ZfkQHVFA2WBgFOvA6V5UHWOhen8oy19JgBKT0BVacNl/aJODesty7f8eDbENxJ0bo0sGwG6mn5L5QWWWWwb4hMGbjUj0yoKLbeGeIeeai6pKIKKglOPnfmh8wkFd0/LdmUxlJ04o8xp296hYPA5VbY4q6ascvZfv0jwrOkEWX4STqRZkgnrzXRqO7j9qc9PcRakLQNTleUzZK62/K2utLyfCcMhJtFStnYYtrHM8lhlseV9ry63PH7V63DJ/ZbtzPXw2WmfZ7AkP2FdILQTdLne0jlYCHEWh1lbymg0smnTJiZNmmTdp9VqSU5OZs2aNfU+Z82aNUyYMKHOvuHDhzNnzpx6y1dWVlJZeepHq6joHF/cF+lo6jbifryywcc/rh7JG9VjAIjRZLPC8NhZZfxqbuuCbyDxwc8BOJl3jIgv+jd43FnVQ3i6+h4AvClnl8edZ5Xxqblt8h1K0ONzAEsiea7jLjb14q6qJ6339xnGYtBU1XvcnYaehE9abt0fPjOZAOr/Qd1ijudR48vW+6sNTxNpyId5p8pcCVzpARR3Bk77Mfjy2rr9D04XGAuPbDt1/7ubLFX+9fEOhSdPO86PYyGj/s8ceh945uip+7/cBweX1F9Wo4XJJ0/d//1RS5+JhjybfSoB+GMSbJ/VcNmn0sGrZu6JJS/Dpi8aLvvYrlPJ28p3YM2HDZd9YAOEtLdsr/kQVrzVcNl7lkFkL8v2hs9gyUsNlx2/ANoMsGxv+x4WPNVw2Vt+goRky/auX2DuQw2Xvelr6HytZXv/H/DT2Z93q+s/gR7/tmwfXg0zb2647D+mQN/xlu2cPfDLvQ2X9Qw6ldxUlZ+aX+ZM7t51E+zAOLjsKUsyF9QWQjuDb/ipBFcIYROqJjd5eXmYTCbCwsLq7A8LC2Pv3r31PicrK6ve8llZWfWWT0lJ4aWXzvEFbEsaLRXKOUbuaN0w1NSkuHOqrN5Nh7bmu63arFBtMqPo6h7nXMdVTjuuoc5xtWhrvjStx9XUfcvPdVyz1h2P06bfrsQdRbE0HelqAjaZFarMZhRt3eOYNHoqFX2dfbVf356entzcK4bWgZ60DvQiaLEvSnlZnTJWurrHwM0Abg10ntQZzr7fUNkz9+vcz1H2zOPqGy6r0Ta+7Jm054jhTOeK1xLIacfVnbvs6T+sWrcmHPdcZTV1XwutDtw8z3He08u6WZKCM8vUnlt72mdY5w4G/5qHNDXlT/t7+ufH3aumpkxruWl1p7Y1OjCctiKydzC0u9zyfK2b5Tw6veWvuzeEdztVtlU83PCZ5fh6L0sy7B1smUhPf8ZCrj4hcPmzDbxmQghbUbVZ6tixY0RFRbF69WqSkpKs+5966imWL1/OunXrznqOXq/nyy+/ZMyYMdZ9H330ES+99BLZ2dlnla+v5iY6OloWzhRCCCEciMM0SwUHB6PT6c5KSrKzswkPD6/3OeHh4U0qbzAYMBhk9lohhBDCVag6llav19OnTx+WLDnVf8FsNrNkyZI6NTmnS0pKqlMeYNGiRQ2WF0IIIYRrUX0o+IQJExg7dix9+/alf//+TJkyhdLSUsaPt3Tsu/3224mKiiIlJQWARx55hMGDB/POO+9wzTXXMHPmTDZu3Mj//vc/NS9DCCGEEHZC9eRm9OjR5Obm8sILL5CVlUXPnj1ZuHChtdNwRkYG2tM6Gw4YMIDvvvuO5557jmeeeYaEhATmzJkjc9wIIYQQArCDeW5aWnPOcyOEEEKI5tGU32+Zv14IIYQQTkWSGyGEEEI4FUluhBBCCOFUJLkRQgghhFOR5EYIIYQQTkWSGyGEEEI4FUluhBBCCOFUJLkRQgghhFOR5EYIIYQQTkX15RdaWu2EzEVFRSpHIoQQQojGqv3dbszCCi6X3BQXFwMQHR2tciRCCCGEaKri4mL8/f3PWcbl1pYym80cO3YMX19fNBqNTY9dVFREdHQ0mZmZLrlulatfP8hrINfv2tcP8hq4+vVD870GiqJQXFxMZGRknQW16+NyNTdarZbWrVs36zn8/Pxc9kMNcv0gr4Fcv2tfP8hr4OrXD83zGpyvxqaWdCgWQgghhFOR5EYIIYQQTkWSGxsyGAxMnjwZg8GgdiiqcPXrB3kN5Ppd+/pBXgNXv36wj9fA5ToUCyGEEMK5Sc2NEEIIIZyKJDdCCCGEcCqS3AghhBDCqUhyI4QQQginIsnNaV588UU0Gk2dW8eOHa2PV1RU8MADD9CqVSt8fHy44YYbyM7OrnOMjIwMrrnmGry8vAgNDeXJJ5+kurq6Tplly5bRu3dvDAYD8fHxzJgxoyUur14rVqxg5MiRREZGotFomDNnTp3HFUXhhRdeICIiAk9PT5KTkzlw4ECdMvn5+dxyyy34+fkREBDAnXfeSUlJSZ0y27dvZ9CgQXh4eBAdHc2bb755Viw//vgjHTt2xMPDg27dujF//nybX++Zznf948aNO+szcdVVV9Up48jXn5KSQr9+/fD19SU0NJTrrruOffv21SnTkp/7qVOnEhsbi4eHB4mJiaxfv97m13ymxrwGQ4YMOetzcN9999Up46ivwccff0z37t2tE64lJSWxYMEC6+PO/v7D+V8DZ37/6/P666+j0Wh49NFHrfsc7nOgCKvJkycrXbp0UY4fP2695ebmWh+/7777lOjoaGXJkiXKxo0blUsuuUQZMGCA9fHq6mqla9euSnJysrJlyxZl/vz5SnBwsDJp0iRrmbS0NMXLy0uZMGGCsnv3buWDDz5QdDqdsnDhwha91lrz589Xnn32WeXnn39WAOWXX36p8/jrr7+u+Pv7K3PmzFG2bdumXHvttUpcXJxSXl5uLXPVVVcpPXr0UNauXausXLlSiY+PV8aMGWN9vLCwUAkLC1NuueUWZefOncr333+veHp6Kp988om1zN9//63odDrlzTffVHbv3q0899xziru7u7Jjxw5Vr3/s2LHKVVddVeczkZ+fX6eMI1//8OHDlS+++ELZuXOnsnXrVuXqq69WYmJilJKSEmuZlvrcz5w5U9Hr9crnn3+u7Nq1S7n77ruVgIAAJTs7W/XXYPDgwcrdd99d53NQWFjoFK/B3LlzlXnz5in79+9X9u3bpzzzzDOKu7u7snPnTkVRnP/9b8xr4Mzv/5nWr1+vxMbGKt27d1ceeeQR635H+xxIcnOayZMnKz169Kj3sYKCAsXd3V358ccfrfv27NmjAMqaNWsURbH8UGq1WiUrK8ta5uOPP1b8/PyUyspKRVEU5amnnlK6dOlS59ijR49Whg8fbuOrabozf9zNZrMSHh6uvPXWW9Z9BQUFisFgUL7//ntFURRl9+7dCqBs2LDBWmbBggWKRqNRjh49qiiKonz00UdKYGCg9TVQFEV5+umnlQ4dOljv33TTTco111xTJ57ExETl3nvvtek1nktDyc2oUaMafI4zXb+iKEpOTo4CKMuXL1cUpWU/9/3791ceeOAB632TyaRERkYqKSkptr/QczjzNVAUy4/b6V/0Z3K21yAwMFCZPn26S77/tWpfA0Vxnfe/uLhYSUhIUBYtWlTnmh3xcyDNUmc4cOAAkZGRtG3blltuuYWMjAwANm3aRFVVFcnJydayHTt2JCYmhjVr1gCwZs0aunXrRlhYmLXM8OHDKSoqYteuXdYypx+jtkztMexJeno6WVlZdeL19/cnMTGxzjUHBATQt29fa5nk5GS0Wi3r1q2zlrnsssvQ6/XWMsOHD2ffvn2cPHnSWsZeX5dly5YRGhpKhw4duP/++zlx4oT1MWe7/sLCQgCCgoKAlvvcG41GNm3aVKeMVqslOTlZ9deg1rfffktwcDBdu3Zl0qRJlJWVWR9zltfAZDIxc+ZMSktLSUpKcsn3/8zXoJYrvP8PPPAA11xzzVlxOuLnwOUWzjyXxMREZsyYQYcOHTh+/DgvvfQSgwYNYufOnWRlZaHX6wkICKjznLCwMLKysgDIysqq88bWPl772LnKFBUVUV5ejqenZzNdXdPVxlxfvKdfT2hoaJ3H3dzcCAoKqlMmLi7urGPUPhYYGNjg61J7DLVcddVV/POf/yQuLo6DBw/yzDPPMGLECNasWYNOp3Oq6zebzTz66KNceumldO3a1RpfS3zuT548iclkqrfM3r17bXaN51PfawBw880306ZNGyIjI9m+fTtPP/00+/bt4+effwYc/zXYsWMHSUlJVFRU4OPjwy+//ELnzp3ZunWry7z/Db0G4PzvP8DMmTPZvHkzGzZsOOsxR/wekOTmNCNGjLBud+/encTERNq0acMPP/xgV0mHaDn//ve/rdvdunWje/futGvXjmXLljFs2DAVI7O9Bx54gJ07d7Jq1Sq1Q1FNQ6/BPffcY93u1q0bERERDBs2jIMHD9KuXbuWDtPmOnTowNatWyksLGT27NmMHTuW5cuXqx1Wi2roNejcubPTv/+ZmZk88sgjLFq0CA8PD7XDsQlpljqHgIAA2rdvT2pqKuHh4RiNRgoKCuqUyc7OJjw8HIDw8PCzeo/X3j9fGT8/P7tLoGpjri/e068nJyenzuPV1dXk5+fb5HWpfdxetG3bluDgYFJTUwHnuf4HH3yQ33//naVLl9K6dWvr/pb63AcHB6PT6ezyNahPYmIiQJ3PgSO/Bnq9nvj4ePr06UNKSgo9evTg/fffd6n3v6HXoD7O9v5v2rSJnJwcevfujZubG25ubixfvpz//ve/uLm5ERYW5nCfA0luzqGkpISDBw8SERFBnz59cHd3Z8mSJdbH9+3bR0ZGhrVdNikpiR07dtT5sVu0aBF+fn7W6s2kpKQ6x6gtc3rbrr2Ii4sjPDy8TrxFRUWsW7euzjUXFBSwadMma5m//voLs9ls/QJISkpixYoVVFVVWcssWrSIDh06EBgYaC3jCK/LkSNHOHHiBBEREYDjX7+iKDz44IP88ssv/PXXX2c1n7XU516v19OnT586ZcxmM0uWLFH9NajP1q1bAep8Dhz5NTiT2WymsrLSJd7/htS+BvVxtvd/2LBh7Nixg61bt1pvffv25ZZbbrFuO9znoEndj53c448/rixbtkxJT09X/v77byU5OVkJDg5WcnJyFEWxDIWLiYlR/vrrL2Xjxo1KUlKSkpSUZH1+7VC4K6+8Utm6dauycOFCJSQkpN6hcE8++aSyZ88eZerUqaoOBS8uLla2bNmibNmyRQGUd999V9myZYty+PBhRVEsQ8EDAgKUX3/9Vdm+fbsyatSoeoeC9+rVS1m3bp2yatUqJSEhoc5Q6IKCAiUsLEy57bbblJ07dyozZ85UvLy8zhoK7ebmprz99tvKnj17lMmTJ7fIUOhzXX9xcbHyxBNPKGvWrFHS09OVxYsXK71791YSEhKUiooKp7j++++/X/H391eWLVtWZ5hrWVmZtUxLfe5nzpypGAwGZcaMGcru3buVe+65RwkICKgz+kKN1yA1NVV5+eWXlY0bNyrp6enKr7/+qrRt21a57LLLnOI1mDhxorJ8+XIlPT1d2b59uzJx4kRFo9Eof/75p6Iozv/+n+81cPb3vyFnjhBztM+BJDenGT16tBIREaHo9XolKipKGT16tJKammp9vLy8XPm///s/JTAwUPHy8lKuv/565fjx43WOcejQIWXEiBGKp6enEhwcrDz++ONKVVVVnTJLly5Vevbsqej1eqVt27bKF1980RKXV6+lS5cqwFm3sWPHKopiGQ7+/PPPK2FhYYrBYFCGDRum7Nu3r84xTpw4oYwZM0bx8fFR/Pz8lPHjxyvFxcV1ymzbtk0ZOHCgYjAYlKioKOX1118/K5YffvhBad++vaLX65UuXboo8+bNa7brrnWu6y8rK1OuvPJKJSQkRHF3d1fatGmj3H333Wf9T+bI11/ftQN1PpMt+bn/4IMPlJiYGEWv1yv9+/dX1q5d2xyXXcf5XoOMjAzlsssuU4KCghSDwaDEx8crTz75ZJ15ThTFcV+DO+64Q2nTpo2i1+uVkJAQZdiwYdbERlGc//1XlHO/Bs7+/jfkzOTG0T4HGkVRlKbV9QghhBBC2C/pcyOEEEIIpyLJjRBCCCGciiQ3QgghhHAqktwIIYQQwqlIciOEEEIIpyLJjRBCCCGciiQ3QgghhHAqktwIIYQQwqlIciOEEM1s2bJlaDQaNBoN1113XZOeO27cOOtz58yZ0yzxCeFsJLkRwsXU/lA2dHvxxRfVDtHmYmNjmTJlitphsG/fPmbMmGG9P27cuLOSndmzZ+Ph4cE777wDwPvvv8/x48dbMEohHJ+b2gEIIVrW6T+Us2bN4oUXXmDfvn3WfT4+PmqE1WSKomAymXBza7mvMaPRiF6vv+Dnh4aGEhAQ0ODj06dP54EHHmDatGmMHz8eAH9/f/z9/S/4nEK4Iqm5EcLFhIeHW2/+/v5oNJo6+2bOnEmnTp3w8PCgY8eOfPTRR9bnHjp0CI1Gww8//MCgQYPw9PSkX79+7N+/nw0bNtC3b198fHwYMWIEubm51ufV1lC89NJLhISE4Ofnx3333YfRaLSWMZvNpKSkEBcXh6enJz169GD27NnWx2ubdhYsWECfPn0wGAysWrWKgwcPMmrUKMLCwvDx8aFfv34sXrzY+rwhQ4Zw+PBhHnvsMWvtFMCLL75Iz54967w2U6ZMITY29qy4X331VSIjI+nQoQMAmZmZ3HTTTQQEBBAUFMSoUaM4dOjQRb0vb775Jg899BAzZ860JjZCiAsjNTdCCKtvv/2WF154gQ8//JBevXqxZcsW7r77bry9vRk7dqy13OTJk5kyZQoxMTHccccd3Hzzzfj6+vL+++/j5eXFTTfdxAsvvMDHH39sfc6SJUvw8PBg2bJlHDp0iPHjx9OqVSteffVVAFJSUvjmm2+YNm0aCQkJrFixgltvvZWQkBAGDx5sPc7EiRN5++23adu2LYGBgWRmZnL11Vfz6quvYjAY+Oqrrxg5ciT79u0jJiaGn3/+mR49enDPPfdw9913N/k1WbJkCX5+fixatAiAqqoqhg8fTlJSEitXrsTNzY3//Oc/XHXVVWzfvv2CanaefvppPvroI37//XeGDRvW5OcLIc7Q5HXEhRBO44svvlD8/f2t99u1a6d89913dcq88sorSlJSkqIoipKenq4AyvTp062Pf//99wqgLFmyxLovJSVF6dChg/X+2LFjlaCgIKW0tNS67+OPP1Z8fHwUk8mkVFRUKF5eXsrq1avrnPvOO+9UxowZoyiKoixdulQBlDlz5pz3urp06aJ88MEH1vtt2rRR3nvvvTplJk+erPTo0aPOvvfee09p06ZNnbjDwsKUyspK676vv/5a6dChg2I2m637KisrFU9PT+WPP/6oN57a2E+ePFln/9ixYxW9Xn/W61cfQPnll1/OWUYIYSE1N0IIAEpLSzl48CB33nlnnRqO6urqs/p8dO/e3bodFhYGQLdu3ersy8nJqfOcHj164OXlZb2flJRESUkJmZmZlJSUUFZWxhVXXFHnOUajkV69etXZ17dv3zr3S0pKePHFF5k3bx7Hjx+nurqa8vJyMjIymnL5DerWrVud2pht27aRmpqKr69vnXIVFRUcPHiwycfv3r07eXl5TJ48mf79+ztMnych7JkkN0IIwJIkAHz66ackJibWeUyn09W57+7ubt2u7cNy5j6z2dzkc8+bN4+oqKg6jxkMhjr3vb2969x/4oknWLRoEW+//Tbx8fF4enpy44031unPUx+tVouiKHX2VVVVnVXuzPOVlJTQp08fvv3227PKhoSEnPOc9YmKimL27NkMHTqUq666igULFpyVOAkhmkaSGyEEYKltiYyMJC0tjVtuucXmx9+2bRvl5eV4enoCsHbtWnx8fIiOjiYoKAiDwUBGRkad/jWN8ffffzNu3Diuv/56wJJ8nNm5V6/XYzKZ6uwLCQkhKysLRVGsCdrWrVvPe77evXsza9YsQkND8fPza1KsDWnTpg3Lly+3JjgLFy6UBEeIiyCjpYQQVi+99BIpKSn897//Zf/+/ezYsYMvvviCd99996KPbTQaufPOO9m9ezfz589n8uTJPPjgg2i1Wnx9fXniiSd47LHH+PLLLzl48CCbN2/mgw8+4MsvvzzncRMSEvj555/ZunUr27Zt4+abbz6r1ig2NpYVK1Zw9OhR8vLyAMsoqtzcXN58800OHjzI1KlTWbBgwXmv45ZbbiE4OJhRo0axcuVK0tPTWbZsGQ8//DBHjhy54NcnOjqaZcuWkZOTw/DhwykqKrrgYwnh6iS5EUJY3XXXXUyfPp0vvviCbt26MXjwYGbMmEFcXNxFH3vYsGEkJCRw2WWXMXr0aK699to6Ewa+8sorPP/886SkpNCpUyeuuuoq5s2bd95zv/vuuwQGBjJgwABGjhzJ8OHD6d27d50yL7/8MocOHaJdu3bWpqNOnTrx0UcfMXXqVHr06MH69et54oknznsdXl5erFixgpiYGP75z3/SqVMn7rzzTioqKi66Jqd169YsW7aMvLw8SXCEuAga5cxGZyGEsLFx48ZRUFDgsssHLFu2jKFDh3Ly5MlzTuJ3LhqNhl9++aXJyzcI4Yqk5kYIIVpI69atGTNmTJOec99998kIKiGaSGpuhBDNztVrbsrLyzl69ChgWd4iPDy80c/NycmxNk9FREScNXpLCHE2SW6EEEII4VSkWUoIIYQQTkWSGyGEEEI4FUluhBBCCOFUJLkRQgghhFOR5EYIIYQQTkWSGyGEEEI4FUluhBBCCOFUJLkRQgghhFP5fy+wQVlQBC/7AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Ts, X2s, '-.', label=r\"He$^+$ / He\")\n",
    "plt.plot(Ts, X3s, '--', label=r\"He$^{++}$ / He\")\n",
    "plt.xlabel('Temperature [K]')\n",
    "plt.ylabel('Abundance')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2822683-337b-4ad2-8d4e-59bd87ce7439",
   "metadata": {},
   "source": [
    "Note that it really needs to get hot for helium to be fully ionized !"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "astro_base",
   "language": "python",
   "name": "astro_base"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}