{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Tolman–Oppenheimer–Volkoff equations (Lindblom's form)\n",
    "\n",
    "We study the log-entalpy formulation of the Oppenheimer–Volkoff equations (also known as Lindblom's form). This formulation (here taken from [Martin Jakob Steil's thesis](https://theorie.ikp.physik.tu-darmstadt.de/nhq/downloads/thesis/master.steil.pdf) in Section 3.2.3) \n",
    "has several advantages w.r.t. the original formulation by Tolman–Oppenheimer–Volkoff:\n",
    "\n",
    "* It directly describes the radius and mass of the star.\n",
    "* It has a known integration domain.\n",
    "* It is not stiff.\n",
    "\n",
    "In this formulation, the state is:\n",
    "\n",
    "$$\n",
    "\\mathbf x=[r^2, m/r]=[x_0, x_1],\n",
    "$$\n",
    "\n",
    "and the parameters are:\n",
    "\n",
    "$$\n",
    "\\mathbf p = [\\eta_K, \\eta_\\Gamma],\n",
    "$$\n",
    "\n",
    "where $K = \\overline K(1+\\eta_K)$ and $\\Gamma = \\overline \\Gamma(1+\\eta_\\Gamma)$ will parameterize the equation of state of a polytrope (a gaseous sphere in hydrodynamic equilibrium).\n",
    "\n",
    "The equations are then written as:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\frac{dx_0}{dh} &= -\\frac{2x_0 (1-2x_1)}{4\\pi x_0 P(h) + x_1}\\\\\n",
    "\\frac{dx_1}{dh} &= \\left(2\\pi \\epsilon(h) - \\frac{x_1}{2x_0}\\right) \\frac{dx_0}{dh} \\\\\n",
    "\\end{cases},\n",
    "$$\n",
    "\n",
    "where an equation of state is needed to define the pressure $P(h)$ and energy density $\\epsilon(h)$ expressions. \n",
    "\n",
    "In the case of a relativistic polytrope here studied, the log-entalpy is introduced and defined as:\n",
    "\n",
    "$$\n",
    "h = \\log{[(\\epsilon+P) / \\rho]},\n",
    "$$\n",
    "\n",
    "where $\\rho$ is the (baryonic) mass density (following a [polytropic process](https://en.wikipedia.org/wiki/Polytropic_process)):\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\rho &= \\left(\\frac{P}{K}\\right)^{1/\\Gamma} \\\\\n",
    "\\epsilon &= \\rho + \\frac{P}{\\Gamma-1}\n",
    "\\end{cases}.\n",
    "$$\n",
    "\n",
    "Hence, in this case the log-entalpy can be written as:\n",
    "\n",
    "$$\n",
    "h = \\log\\left(1 + \\frac{\\Gamma}{\\Gamma-1} K^{\\frac 1\\Gamma} P^{1-\\frac 1\\Gamma}\\right),\n",
    "$$\n",
    "\n",
    "which can be inverted to yield the $P(h)$ form:\n",
    "\n",
    "$$\n",
    "P(h) = K^{-\\frac{1}{\\Gamma-1}}\\left[\\frac{\\Gamma-1}{\\Gamma}\\left(e^h-1\\right)\\right]^{\\frac{\\Gamma}{\\Gamma-1}}.\n",
    "$$\n",
    "\n",
    "\n",
    "**Note**: the dynamics is singular at $h=0$ (of a removable kind), so we will need to define the initial conditions at some small $\\delta h$ where the singularity is not present.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import heyoka as hy\n",
    "import numpy as np\n",
    "from scipy import optimize\n",
    "from copy import deepcopy\n",
    "\n",
    "import time\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0 - Preamble \n",
    "Let us start by introducing common definitions and the reference numerical values we will be using for the parameters $K, \\Gamma$ of the polytrope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# State variables\n",
    "x0, x1 = hy.make_vars(\"x0\", \"x1\")\n",
    "\n",
    "# The independent variable is log-enthalpy not time, thus we rename it as follows:\n",
    "h = hy.time\n",
    "\n",
    "# Parameters values\n",
    "Gamma_v = 2.0\n",
    "K_v = 100.0\n",
    "\n",
    "# Parameters expressions\n",
    "Gamma = Gamma_v * (1.0 + hy.par[0])\n",
    "K = K_v * (1.0 + hy.par[1])\n",
    "\n",
    "# Other expressions\n",
    "P = K ** (-1.0 / (Gamma - 1.0)) * ((Gamma - 1.0) / Gamma * (hy.exp(h) - 1.0)) ** (\n",
    "    Gamma / (Gamma - 1.0)\n",
    ")\n",
    "rho = (P / K) ** (1.0 / Gamma)\n",
    "eps = rho + P / (Gamma - 1.0)\n",
    "\n",
    "# The dynamics\n",
    "dx0dh = -2.0 * x0 * (1.0 - 2.0 * x1) / (4.0 * np.pi * x0 * P + x1)\n",
    "dyn = [\n",
    "    (x0, dx0dh),\n",
    "    (x1, (2.0 * np.pi * eps - x1 / 2.0 / x0) * dx0dh),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 - Standard Numerical Solution\n",
    "The Tolman–Oppenheimer–Volkoff equations in the Lindblom form are numerically well behaved. \n",
    "\n",
    "We here show a straightforward numerical integration using the Taylor adaptive scheme of *heyoka*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We instantiate the numerical ODE solver\n",
    "ta = hy.taylor_adaptive(dyn, state=[1.0, 1.0], time=0.1, tol=1e-16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the conditions at the star center fixing the pressure $P_c$, and thus computing the corresponding density $\\rho_c$, energy density $\\epsilon_c$ and log-entalpy $h_c$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Conditions at the star center:\n",
      "\n",
      "Pressure:  0.01\n",
      "Energy density:  0.02\n",
      "Density:  0.01\n",
      "log-entalpy:  1.0986122886681098\n"
     ]
    }
   ],
   "source": [
    "# Pressure (change this value to then approximate different part of the M-R curve)\n",
    "P_c = 1e-2\n",
    "\n",
    "# Density\n",
    "rho_c = (P_c / K_v) ** (1.0 / Gamma_v)\n",
    "\n",
    "# Energy density\n",
    "eps_c = rho_c + P_c / (Gamma_v - 1.0)\n",
    "\n",
    "# log-entalpy\n",
    "h_c = np.log(\n",
    "    1 + Gamma_v / (Gamma_v - 1) * K_v ** (1 / Gamma_v) * P_c ** (1.0 - 1.0 / Gamma_v)\n",
    ")\n",
    "\n",
    "print(\"Conditions at the star center:\")\n",
    "print(\"\\nPressure: \", P_c)\n",
    "print(\"Energy density: \", eps_c)\n",
    "print(\"Density: \", rho_c)\n",
    "print(\"log-entalpy: \", h_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the conditions at the star center (c) we derive the initial conditions (0) for the numerical integration offsetting the entalpy (the independent variable) by a small amount, thus avoiding the numerical (removable) singularity in the dynamics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since the dynamics is singular at h_c, we define the initial conditions at h_0 = h_c - dh\n",
    "dh = 1e-6  # note that if this gets too small the variational equations may have numerical issues at h_0.\n",
    "h_0 = h_c - dh\n",
    "\n",
    "# Initial conditions at h_0 (from Steil thesis). Could be improved at second order.\n",
    "x0_0 = 3.0 / (2.0 * np.pi * (3.0 * P_c + eps_c)) * dh\n",
    "x1_0 = 2.0 * eps_c / (3 * P_c + eps_c) * dh\n",
    "y_0 = [x0_0, x1_0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to perform a numerical integration, and extract from the result the star radius and mass:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total time to propagate: 0.00038123130798339844\n",
      "Outcome (should be time_limit if all went well):  time_limit\n",
      "Final value for the stellar radius is:  5.511134830596981\n",
      "Final value for the stellar mass is:  1.3684613613885503\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "ta.time = h_0\n",
    "ta.state[:] = y_0\n",
    "ta.pars[:] = [0.0, 0.0]\n",
    "time_grid = np.linspace(h_0, 0.0, 100)\n",
    "sol = ta.propagate_grid(time_grid)\n",
    "time_cost = time.time() - start\n",
    "print(\"Total time to propagate:\", time_cost)\n",
    "print(\"Outcome (should be time_limit if all went well): \", str(sol[0]).split(\".\")[1])\n",
    "\n",
    "# The mass and radius are easily computed from the last value reached\n",
    "rf = np.sqrt(sol[-1][:, 0][-1])\n",
    "Mf = sol[-1][:, 1][-1] * rf\n",
    "\n",
    "print(\"Final value for the stellar radius is: \", rf)\n",
    "print(\"Final value for the stellar mass is: \", Mf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we now plot the results of the numerical integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHVCAYAAACXAw0nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB9FElEQVR4nO3dd1hT59sH8G/YG2WDbEURcAHubStqnbWto9bRYWtdVevPau1Q26pddrpaZ4ejVu3SWrFVHDgRFQU3AiJDUDaEkDzvH9S8RVAZCScJ3891cV3m5Elyn5zk9s45z5AJIQSIiIiISK8ZSR0AEREREdUdizoiIiIiA8CijoiIiMgAsKgjIiIiMgAs6oiIiIgMAIs6IiIiIgPAoo6IiIjIALCoIyIiIjIAJlIHICWVSoVbt27B1tYWMplM6nCIqBqEEMjPz4eHhweMjPi79EGY34j0U11yXIMu6m7dugUvLy+pwyCiWkhJSYGnp6fUYegs5jci/VabHNegizpbW1sA5W+cnZ2dxNEQUXXk5eXBy8tL/f3VFytWrMDHH3+MtLQ0BAcH4/PPP0f37t2rbHv48GG88cYbuHjxIoqKiuDj44NXXnkFM2fOrPbrMb8R6ae65LgGXdTduyRhZ2fHpEekZ/TpkuLWrVsxY8YMrFixAl27dsXq1asxYMAAxMfHw9vbu1J7a2trTJ06Fa1bt4a1tTUOHz6MV155BdbW1nj55Zer9ZrMb0T6rTY5TiaEEFqIRS/k5eXB3t4eubm5THpEekIfv7cdO3ZEaGgoVq5cqd7WsmVLDBs2DEuWLKnWcwwfPhzW1tb4/vvvq9VeH98nIqrbd5e9jImItKi0tBQxMTGIiIiosD0iIgLR0dHVeo7Y2FhER0ejZ8+eD2wjl8uRl5dX4Y+IGhYWdUREWpSVlQWlUglXV9cK211dXZGenv7Qx3p6esLc3Bzh4eGYMmUKXnrppQe2XbJkCezt7dV/HCRB1PDobJ+6lStXYuXKlbhx4wYAIDg4GO+88w4GDBgAoHzI78KFC/HNN9/g7t276NixI5YvX47g4GAJoyai2oi7mYtfz6RicBsPtPFqJHU4WnF//xghxCP7zBw6dAgFBQU4duwY5s6di2bNmmH06NFVtp03bx5mzZqlvn2vszWRrhJCoEBehtv5ctwpLEVOkQJ5JQoUlipRUqqEQqWCUilgZCSDkUwGMxMjWJoaw8bCBI0sTeFgbQYXO3M4WpvD2Eh/+thqk84WdZ6enli6dCmaNWsGANi4cSOGDh2K2NhYBAcH46OPPsKyZcuwYcMGNG/eHO+//z769u2LS5cu6d2oOKKG7tczqVhzOBHZhaX4bGRbqcPRKCcnJxgbG1c6K5eZmVnp7N39/Pz8AACtWrVCRkYGFixY8MCiztzcHObm5poJmkiD7hSW4lJ6Pq5m5uN6ViGSsotw824RUu8Wo7BUWefnNzGSoUljS3g7WKGpsw0CXG0Q7GGPQDdbWJgaa2AP9IfOFnWDBw+ucPuDDz7AypUrcezYMQQFBeHzzz/H/PnzMXz4cADlRZ+rqys2bdqEV155pcrnlMvlkMvl6tvsc0IkPSEE9lwoL3j6BbtJHI3mmZmZISwsDJGRkXjyySfV2yMjIzF06NBqP48QokL+ItJFOUWliE3JwZnkHMSl5uLCrVxk5D38c2tjbgJHGzM0sjKDnYUJrM1MYGlmDBMjGYyNZBACUAqB0jIVikqVKJArkFOkQHZhKbIL5ChTCSRlFyEpuwiHrmSpn9fESIYgDzu093VAJ39HdG7qCBtznS17NEIv9k6pVGLbtm0oLCxE586dkZiYiPT09Aodj83NzdGzZ09ER0c/sKhbsmQJFi5cWF9hE1E1XLiVh5t3i2FhaoSezZ2lDkcrZs2ahbFjxyI8PBydO3fGN998g+TkZEyaNAlA+aXT1NRUfPfddwCA5cuXw9vbG4GBgQDK56375JNPMG3aNMn2gagqOUWlOHotG9HXsnEi8Q4uZeRX2c7bwQrNXW3g72wDX0dreDlYwqORJdztLWBlVvtSpEypQma+HCl3inAjuxBXMwtwKaMA51NzcaewFOdu5uLczVysPZwIEyMZwn0bIyLIDf1D3ODRyLLWr6urdLqoi4uLQ+fOnVFSUgIbGxvs3LkTQUFB6hFjVXU8TkpKeuDzsc8Jke7Zc778LF2v5i6wNDPMSyUjR45EdnY2Fi1ahLS0NISEhGD37t3w8fEBAKSlpSE5OVndXqVSYd68eUhMTISJiQmaNm2KpUuXPvAHK1F9EUIgLjUXfydk4sDl2zh3Mwf3T4zm52SNdl6N0NrTHq087RHoZgdrLZ0hMzE2gkej8gKxo79jhThTc4oRk3QXxxPvIPpqFm5kF+HY9Ts4dv0OFv0Rjw6+DngytAkGtXaHrYWpVuKrbzo9T11paSmSk5ORk5OD7du3Y82aNYiKikJOTg66du2KW7duwd3dXd1+4sSJSElJwZ49e6r1/JzHiUh6jy+LwtXMAnwxqi2Gtm3yyPb83lYP3yfSFKVK4ETiHfx5Pg17L2QgPa+kwv0BLjbo2swJnfwdEO7rACcb3ezbmZRdiL8TMrHnQjpO3rijLkatzIwxtG0TjO/ig0A36b8rdfnu6vSZOjMzM/VAifDwcJw8eRJffPEF3njjDQBAenp6haKuOh2PiUh3XM3Mx9XMApgay9A70EXqcIjoX0IIxKbk4Lczt/DHuTRkFfx/vzgrM2P0CHBG70Bn9GzuAjd7CwkjrT4fR2u80M0PL3TzQ3puCX49k4ptMTdxNbMAm08kY/OJZPRq4YxXezatcNZPn+h0UXe/ex2F/fz84ObmhsjISLRr1w5A+Vm9qKgofPjhhxJHSUTVde/Sa9dmTrAzkMsfRPosLbcYO06nYnvMTVzPKlRvb2RlioggVwwIcUfnpo56P6rUzd4Cr/Rsipd7+ON44h18d/QG9pxPx4FLt3Hg0m10beaI1yNaINS7sdSh1ojOFnVvvvkmBgwYAC8vL+Tn52PLli04cOAA9uzZA5lMhhkzZmDx4sUICAhAQEAAFi9eDCsrKzz77LNSh05E1bQrrryo62+Ao16J9IVSJXDgUiY2HU/G/kuZUP17WdLS1Bj9gl0xtG0TdAtwgqmx4a1XIJPJ0MnfEZ38HZGUXYhvDl7HT6dScORqNo5cjcag1u6Y90RLNNGTQRU6W9RlZGRg7NixSEtLg729PVq3bo09e/agb9++AIA5c+aguLgYkydPVk8+vHfvXs5RR6QnrmTkIyEtD6bGMoOcyoRI1+UWKbD5ZDJ+OJaEm3eL1ds7+DrgmXBPDGjlbvBTgPyXj6M1PniyFSb1bIov/76Cn0/fxB/n0hAZn4HpjwXg5R7+Ol/Y6vRACW1jR2Ii6Xy69xK++ucqHm/pgjXj21f7cfzeVg/fJ3qQG1mFWHs4ET/H3ESxonzyX3tLUzwT5onRHb3R1NlG4gh1w4VbuVj0ezyOJ94BAAS62eLjp9uglae9Vl/XYAdKEJFhEkLg1zO3AACD23hIHA1Rw3A+NRcrDlzFn+fT1SM/A91s8UJXPwxp66H3/eQ0LdjDHlte7oSdsal47494XEzPx5MrjmBm3+aY1LOpTi5NxqKOiOrd2Zu5SL5TBEtTY/QN4oh1Im2KSbqLL/++gqjLt9XberVwxsvd/dG5qeMj1yBuyGQyGYaHeqJnc2e89ct5/Hk+HR//dQmHr2Thy9Ht4GyrW9O3sKgjonr365lUAEDfINc6zSZPRA92JiUHyyIv4+C/xZyRDBjSxgOv9mqGFm7sf14TjjbmWDEmFD/H3MSC3y7g6PVsDP7qMFY8F6pTI2SZTYmoXpUpVfj9bBqA8v9giEizLqXn45O9lxAZnwEAMDaS4anQJpjSuxl8HK0ljk5/yWQyPBPuhXbejfDK9zG4drsQo1Yfw8fPtK7WxOn1gUUdEdWrqMu3kVUgh6O1GXq2MMy1XomkkJZbjGV7L2P76ZtQifIzc8NDPTG9TwC8Ha2kDs9gNHOxxa9Tu+H1n87grwsZeG3LGaTcKcKU3s0kv5TNoo6I6tW2UzcBAMPaNdH56QGI9EGhvAyro67hm0PXUaJQAQAGhLjh9YgWaObCkazaYGNugpVjwrDkzwR8eygRn+y9jDuFCrw9qKWkhR2LOiKqN3cKS/H3xfJLQs+Ee0ocDZF+uzeKfMmfCcjIK1/Gq71vY8x7oqVO9fMyVEZGMswfGIQmjSyx4Pd4rDuSiGJFGd4f1kqykbEs6oio3vwSmwqFUqBVE3udWDibSF/F38rDO7+ex6mkuwAALwdLvDmgJfqHuEl+CbChmdDVD9bmJnhj+zlsPpECAFj8ZCtJjgOLOiKqN9tiyi+98iwdUe0UyMvwWeRlbIi+AaVKwNLUGFP7NMOL3fw4z5yEngn3grmpMWZsicXmEymwNjPB/IH1fymWRR0R1YtzN3OQkJYHM2MjjnolqoW9F9Lx7m8XkJZbAgAY2Modbw1qCXd7/ViX1NANaeOBklIl5mw/hzWHE9HY2gxTejer1xhY1BFRvfjhWBIAYGBrdzSyMpM4GiL9kZlfgnd/vYA/z6cDALwdrPDesBD0bM7R47pmRHsvFMjLsOiPeHz81yV4O1jV66o5LOqISOtyixX47Wz5smBjOnpLHA2RfhBCYGdsKhb+Ho/cYgWMjWR4uYc/pvcJgKUZL7Xqqhe6+eHm3WKsO5KI17edhUcjS4T51M/AFRZ1RKR1O07fRIlChUA323pLbkT6LDO/BPN3nldPIBzsYYePnm6NYA/tLiZPmjF/YEsk3ynCvoQMTPohBrumd4OLrYXWX5eTRBGRVgkh8OPxZADlZ+k4Mo/o4facT0e/zw4iMj4DpsYyzI5ojl+mdGVBp0eMjWT4cnRbNHe1we18OaZtikWZUqX112VRR0RadfR6Nq5mFsDKzBjD2unGUjpEuqhAXob/bTuLST/E4G6RAi3d7fDb1G6Y2ieAE3XrISszE6x8LgzWZsY4nngHn+y9rPXX5KeEiLRq/ZEbAIAn2zWBrYWptMEQ6aizKTkY9OUhbIu5CZkMmNSzKX6d0hUt3Tmfoz5r6myDj55uAwD45uA1XL9doNXXY586ItKapOxC7Eso7xP0fFc/iaMh0j0qlcDaw4n4cM9FlKkEPOwt8NnItujo7yh1aKQhA1u741JGAMJ8GsPfWbvLtrGoIyKt2RB9A0IAvVo4cw1KovvcKSzF7G1n8c/FTADAE63csOTJ1rC34hltQzOrb/N6eR0WdUSkFfklCmw7Vb6CxAs8S0dUQWzyXUz58TRu5ZbAzMQI7w4OwrMdOJCI6oZFHRFpxdaTKSiQlyHAxQbdA5ykDodIJwgh8P2xJLz3RzwUSgE/J2ssfzYUQR7sO0d1x6KOiDROXqbEmkOJAIAXu/nx7AMRgBKFEm/ujMOO06kAyi+3fvhUaw4gIo1hUUdEGvdr7C2k55XA1c4cT4ZyGhOi1JxivPL9KZxPzYOxkQxz+wfipe78wUOaxaKOiDRKqRJYdfAaAOClbv4wN+FyRtSwnUi8g1d/iEF2YSkcrM3w9bPt0KUpuySQ5rGoIyKNioxPx/XbhbCzMMForvNKDdzWk8mYv/M8ylQCQe52+GZcGDwbW0kdFhkoFnVEpDFCCHz1z1UAwLjOvrAxZ4qhhkmpElj6ZwK+/bdv6cDW7vj46dawMuN3grSHny4i0ph9CZm4cCsP1mbGeKEbpzGhhqmotAyvbTmDyPjyibdfeywAMx4PYP850joWdUSkEUIIfL6vfG3D8V184WBtJnFERPXvdr4cL208ibM3c2FmYoSPn26NoW05WIjqB4s6ItKI/56le6m7v9ThENW7a7cLMH7dCdy8W4zGVqb4dlw4wn0dpA6LGhAWdURUZyqVwGeRPEtHDdfp5Lt4ccNJ3C1SwMfRChue7wA/J2upw6IGhkUdEdXZH3FpiE/Lg625Cc/SUYOz/2ImXv0xBiUKFdp42mPthPZwsjGXOixqgFjUEVGdKJQqfLr3EgDg5R7+PEtHDcrO2JuYve0clCqBXi2csWJMKEe4kmT4ySOiOtlyMgVJ2UVwsjHjiFdqUDZG38C7v10AADzZrgk+ero1TI2NJI6KGjIWdURUa4XyMnz59xUAwLQ+AbDmvHTUAAghsHz/VXyyt7wf6fNdffH2wCAYGXHKEpIWMzAR1drqqGu4nS+Ht4MVRnfg6hFk+IQQ+OivS1h5oHwpPM5BR7qERR0R1cqtnGJ8c+g6AGDegECYmfCyExk2IQQW/RGP9UduAADeGtiSA4NIp7CoI6Ja+fivSyhRqNDB1wH9Q9ykDodIq1QqgXd+O48fjiUDAN4bFoKxnXwkjoqoIhZ1RFRjscl3sTM2FQDw1qCWvPREBk2lEnj71/P48XgyZDLgw6daY0S4l9RhEVXCoo6IakT5739wAPBUqCdaezaSNiAiLRJC4K1fz2PTvwXdJ0+3wVNhnlKHRVQldoIhohrZfCIZ51PzYGthgrkDAqUOR2+sWLECfn5+sLCwQFhYGA4dOvTAtjt27EDfvn3h7OwMOzs7dO7cGX/99Vc9RktAeUH37m8X1AXdp8+woCPdxqKOiKrtTmEpPv6rfKLhWX2bw9mWs+ZXx9atWzFjxgzMnz8fsbGx6N69OwYMGIDk5OQq2x88eBB9+/bF7t27ERMTg969e2Pw4MGIjY2t58gbLiEE3t+VgO+OJkEmAz5+ug2Gh7KgI90mE0IIqYOQSl5eHuzt7ZGbmws7OzupwyHSef/bdhbbYm4i0M0Wf0zrBhMJJlrVx+9tx44dERoaipUrV6q3tWzZEsOGDcOSJUuq9RzBwcEYOXIk3nnnnSrvl8vlkMvl6tt5eXnw8vLSq/dJl3zy1yV8vf8qAODDp1phZHtO2UP1oy45TmfP1C1ZsgTt27eHra0tXFxcMGzYMFy6dKlCGyEEFixYAA8PD1haWqJXr164cOGCRBETGbboa1nYFnMTAPDBkyGSFHT6qLS0FDExMYiIiKiwPSIiAtHR0dV6DpVKhfz8fDg4ODywzZIlS2Bvb6/+8/JiR/7aWr7/qrqge29oMAs60hs6m5WjoqIwZcoUHDt2DJGRkSgrK0NERAQKCwvVbT766CMsW7YMX3/9NU6ePAk3Nzf07dsX+fn5EkZOZHhKFErM31k+OOK5Tt4I83lwcUEVZWVlQalUwtXVtcJ2V1dXpKenV+s5Pv30UxQWFmLEiBEPbDNv3jzk5uaq/1JSUuoUd0P13dEb6i4G8wYEYmxnX2kDIqoBnR39umfPngq3169fDxcXF8TExKBHjx4QQuDzzz/H/PnzMXz4cADAxo0b4erqik2bNuGVV16RImwig/TVP1eQmFUIF1tzzOnPwRG1cf+0L0KIak0Fs3nzZixYsAC//vorXFxcHtjO3Nwc5ubs41gXv8Sm4p1fy6/2TO/TDK/0bCpxREQ1o7Nn6u6Xm5sLAOrLD4mJiUhPT69wScPc3Bw9e/Z84CUNuVyOvLy8Cn9E9HBxN3OxKqp85YiFQ4JhZ2EqcUTao1QqsX37do2e7XdycoKxsXGls3KZmZmVzt7db+vWrXjxxRfx008/4fHHH9dYTFTZ/ouZmL3tLABgQhdfzOzbXOKIiGpOL4o6IQRmzZqFbt26ISQkBADUCbImlzTY54SoZuRlSszedhZKlcDA1u4Y0Mpd6pC0ytjYGM899xxu376tsec0MzNDWFgYIiMjK2yPjIxEly5dHvi4zZs3Y8KECdi0aRMGDhyosXiospiku3j1xxiUqQSebNcE7wwK4oTapJf0oqibOnUqzp07h82bN1e6ryaXNNjnhKhmvvr7Ki5l5MPR2gyLhgRLHU696NChAxITEzX6nLNmzcKaNWuwbt06JCQkYObMmUhOTsakSZMAlOemcePGqdtv3rwZ48aNw6effopOnTohPT0d6enp6isWpDlXM/PxwoaTKFGo0KuFMz56ujWMjFjQkX7S+aJu2rRp+O2337B//354ev7/HEFubuVrTdbkkoa5uTns7Owq/BFR1WKS7mLFgfIRgIuGhsDRpmH015o+fTrefPNNjf7oGzlyJD7//HMsWrQIbdu2xcGDB7F79274+JSvHZqWllZhzrrVq1ejrKwMU6ZMgbu7u/rvtdde01hMBGTklWD8upPILVagrVcjrBgTClOO6iY9prPz1AkhMG3aNOzcuRMHDhxAQEBApfs9PDwwc+ZMzJkzB0D51AEuLi748MMPqzVQQh/nuyKqD4XyMjzx5SEkZRdhWFsPfD6qndQhqWn7e2tkVP6fuo2NDYYMGYJevXqhXbt2aNWqFczMzDT+etrC/PZwBfIyjFh1FPFpefBzssb2V7vAwVp/ji8Zrrp8d2s8+rW4uBh37txBkyZNKmy/cOECgoM1d3lmypQp2LRpE3799VfY2tqqz8jZ29vD0tISMpkMM2bMwOLFixEQEICAgAAsXrwYVlZWePbZZzUWB1FD9N4f8UjKLoKHvQUWDg2ROpx6df36dZw9exZnz57FmTNnsGTJEty4cQPGxsYIDAzEuXPnpA6R6kihVGHyj6cRn5YHJxszbHy+Aws6Mgg1Kup+/vlnzJw5Ew4ODhBC4Ntvv0XHjh0BAGPHjsXp06c1Fti9mdd79epVYfv69esxYcIEAMCcOXNQXFyMyZMn4+7du+jYsSP27t0LW1tbjcVB1NDsOpeGLSdTyte6HNEW9paGO9r1v958800MGzYMHTp0gK+vL4YOHaq+Lz8/H2fOnGFBZwCEEHjn1/M4ePk2LE2NsW5Ce3g7WkkdFpFG1Kioe//993H69Gk4Ozvj1KlTGD9+PObPn49nn30Wmr6KW53nk8lkWLBgARYsWKDR1yZqqFLuFGHujvLCZXKvpujc1FHiiOpPWloaBg0aBGNjYwwePBjDhg3DY489BnNzc9ja2qJ79+7o3r271GFSHa2Kuo7NJ8p/tHw5uh1aezaSOiQijalRUadQKODs7AwACA8Px8GDBzF8+HBcvXqVw7+J9JxCqcL0LbHILylDmE9jzHi8Yc3TtX79egghcPjwYfz++++YNWsWUlNT0bdvXwwZMgSDBg2Ck5OT1GFSHew5n4YP91wEALwzKAh9gx4+TyCRvqnRMB8XF5cKlx8cHR0RGRmJhIQEXpYg0nNLdl9EbHIO7CxM8MWotg1yFKBMJkP37t3x0Ucf4eLFizhx4gQ6deqEb7/9Fk2aNEGPHj3wySefIDU1VepQqYbibuZixtYzAIDxnX3wfFc/aQMi0oIaZe3vv/++wjI1SqUSv//+O1avXo2oqCiNB0dE9eOPc7ew7kj53GyfjmgLz8bsYwQALVu2xJw5c3DkyBHcvHkT48ePx6FDh6qcM5N0V2ZeCSZ+dwolChV6NHfG24OCpA6JSCvqPKWJpaUlLly4AH9/f03FVG845J8IuJpZgKFfH0ZhqRKTejbF3AG6vbartr63CxYsQGhoKMLCwiqN7tdHzG/lShRKjPzmGM6m5KCZiw12TO5i0Evdkf6r1ylN7ndv9nV9LOqIGrq8EgVe/u4UCkuV6OjngNkRDasf3X8tWrRI3TfYyckJYWFhCA0NVRd69yYKJv0hhMCbO+JwNiUH9pamWDMunAUdGbQ6d5rRxuzrRKR9KpXAzC1ncD2rEO72Flg+JhQmDbAf3T3t27dHkyZN8NZbb2HBggVo0qQJdu/ejdGjR8Pf3x9OTk6IiIiQOkyqgbWHE7EjNhXGRjKsHBMKXydrqUMi0qo6n6l75plnAADBwcF6Pfs6UUPzaeQl/H0xE+YmRvhmbDicGsgyYA9y/PhxbNiwAW+++SbatWuHzz77DM2bN4dCocC5c+dw+vRpxMbGSh0mVdPhK1lYvDsBAPDWwJbo0owjl8nw1flneWJiInbu3InZs2ejqKgIS5YsQYcOHWBjY4PWrVtrIkYi0rAdp29i+f5rAIAlw1uhlae9xBHphgkTJuDy5csIDg5GeHg4/ve//0EulyMsLAwTJ07EihUrpA6RqiHlThGmbj4NlQCeCvXEhC6+UodEVC/qfKbOx8cHPj4+nH2dSE+cunEHc7fHASifYHh4qKfEEekWGxsbfPTRR5g4cSJmzZqFZs2aYfHixXjhhRekDo2qoUShxKs/xiCnSIHWnvb44MkQzqNKDYZWOtDcm319ypQp2nh6Iqql67cLMPG7UyhVqtA/2A2zI1pIHZJOUigUKC4uxqhRo+Dt7Y2JEyfizp07UodFjyCEwNu/nMf51Dw4WJth5XNhsDA1ljosonpT5zN1RKQfsgrkmLD+JO4WKdDG0x7LRraBkRHPYNzzwQcfIC4uDnFxcbh8+TKsra3RunVrdOzYEa+88grs7XmJWtdtOZmCbTE3YSQDvhrdDk0aWUodElG9YlFH1AAUysvw4oaTSL5TBC8HS6yd0B5WZvz6/9fbb78NX19fTJgwAaNHj0ZAQIDUIVENnE/Nxbu/XQAAzO7XAl05MIIaoIY7fwFRA1FapsKkH2Jw9mYuGlmZYsPzHRr8SNeqdOvWDdnZ2ViwYAHatm2Lzp07Y+rUqVi3bh3Onj0LpVIpdYj0ALlFCrz6YwxKy1R4vKUrJvVoKnVIRJLQSFF35MgRyOXySv8mImkpVQKvbzuLQ1eyYGVmjA3Pd0BTZxupw9JJBw8eRG5uLi5duoS1a9eie/fuSEhIwOzZs9GuXTvY2NigQ4cOUodJ9xFC4H8/n0XKnWJ4OVji02fYrYAaLo1cfxkwYADOnDkDf3//Cv8mIukIITB/Zxx+P3sLpsYyrHouDG29Gkkdls4LCAhAQEAARo0apd6WmJiIU6dOcZ46HbT+yA3sjc+AmbERVjwbBnsrrhhBDZdGirr/Lh9bx6VkiUgDhBB4f1cCtpxMgZEM+GxkW/Ro7ix1WHrLz88Pfn5+6snWSTecTcnBkj/LJxieP7Al51ukBo996ogMjBACH+65hLWHEwEAHz7VGoNae0gcFZFm5ZUoMHXzaSiUAgNC3DCuM9fmJWJRR2RAhBD4ZO8lrIoqXy3ivaHBeCbcS+KoiDSrvGvBeaTcKYZnY0ssfao1JxgmAos6IoNxr6C7t/zXgsFBGNvZV9qgiLRgW8xN/H72FoyNZPhydDvYW7IfHRHAoo7IIAghsHh3grqge3tQECZ09ZM4Kv3FEf266/rtArz7a/l8dK9HNEeod2OJIyLSHSzqiPScSiXw7m8X8O2h8j50i4YG48VuLOjqYsCAAUhNTa30b5JWaZkKr205g2KFEl2aOnI+OqL7cEp5Ij1WplRhzs/nsCM2FTIZsPjJVhjdwVvqsPQeR/TrpmWRlxGXWj6J9rIRbTkfHdF9NFLUvfnmm3BwcKj0byLSnhKFEtM2xyIyPgPGRjIsG9EGQ9s2kTosIq04dj0bqw+Wdy9YOrw13OwtJI6ISPdopKibN29elf8mIu3ILVZg4sZTOHHjDsxMjLD82VD0DXKVOiwircgrUeD1n85CCGBkuBf6h7hJHRKRTuLlVyI9k5ZbjAnrTuJSRj5szU3w7fhwdPJ3lDosIq1Z8NsFpOYUw9vBCm8PDpI6HCKdxaKOSI/E38rDCxtOIj2vBC625tj4Qge0dLeTOiwirdlzPg07Tqf+uzJKG9iY878togfht4NIT+y/lIlpm2JRIC9DgIsN1j/fHp6NraQOi0hrsgrkeHPneQDAq72aIsyH/bWJHkbjRV1OTg7++usvpKamQiaTwd3dHf369UPjxpxLiKi2NhxJxKI/4qESQCd/B6x+LpwLl5NBE0LgzR1xuFNYikA3W7z2WHOpQyLSeRqdp27t2rXo0KEDjh07BpVKBaVSiWPHjqFTp05Yu3atJl+KqEFQKFWYvzMOC34vL+ieCfPEdy90ZEGnZRzRL71fzqRib3wGTI1lWDaiLcxMOK0q0aPIhAYnYWrRogViYmJgY2NTYXt+fj7CwsJw+fJlTb2URuTl5cHe3h65ubmws2O/JNItdwpL8eoPMTieeAcyGTCnXyAm9fRv8Gtc8ntbPfr8PmXmleDxZVHIKynD7IjmmNonQOqQiOpNXb67Gr38KpPJUFBQUKmoKygoaPD/ERHVxPnUXLzyfQxSc4phY26CL0a1xWMtOWUJGT4hBN7cGYe8kjK0amKPST25agRRdWm0qPvkk0/Qs2dPhISEoEmT8klQb968iQsXLuDTTz/V5EsRGayfY25i/s44yMtU8HW0wrfjwhHgait1WET14pczqdiXkAlTYxk+eaYNTIx52ZWouupc1OXn58PWtvw/nEGDBmHAgAE4ceIEbt26BSEEmjRpgg4dOsDY2LjOwRIZshKFEgt/j8fmE8kAgD6BLvhsZFvYW7L/HDUMt/PlWPh7PADgtccC0MKNP2aIaqLORV337t2xZ88euLmVz/BtbGyMzp071zkwooYkKbsQUzadxvnUPMhk5f+hTe8TwLUtdVRMTAzCwsKkDsPgLPjtAnKKFAhyt8MrvOxKVGN1Pq8dHh6Ojh074uLFixW2x8bG4oknnqjr0xMZvF3n0jDoy8M4n5qHxlam2PB8B8x4vDkLOh325JNPSh2CwdlzPh274tJgbCTDR0+3hikvuxLVWJ3P1K1ZswYLFy5Et27d8Msvv8DFxQVvvfUWtm/fjiFDhmgiRiKDVFyqxKI/LmDziRQAQLhPY3z1bDu421tKHBkBwIgRI6rcLoTAnTt36jkaw5ZbrMA7v5ZPMvxKD3+ENLGXOCIi/aSRgRLvvvsuzMzM0LdvXyiVSvTr1w8nT55EaGioJp6eyOBcuJWL17acwdXMAshkwKs9m2JW3+bsFK5D9u3bh++//77SaH4hBA4ePChRVIbpwz0XkZkvh7+TNaY/xulLiGqrzkVdWloalixZgjVr1iAoKAgXL17EqFGjWNARVUGlEvj20HV8svcSFEoBF1tzfDayLbo2c5I6NLpPr169YGNjg549e1a6r127dhJEZJhO3riDTcfLBwctHt4KFqYcVEdUW3Uu6vz9/REYGIht27Zh4MCB+OuvvzBixAjcvHkTb7zxhiZiJDIIKXeKMHvbWRxPLL901zfIFUuHt4KjjbnEkdF/3RvRv2PHjge22bNnTz1GZLjkZUrM2xEHABgZ7oVO/o4SR0Sk3+p8rWf9+vWIjY3FwIEDAQD9+vXD/v378cUXX2Dy5Ml1DpBI3wkhsPVkMgZ8cQjHE+/AyswYS4a3wjdjw1jQ6aDu3bsjPT1d48+7YsUK+Pn5wcLCAmFhYTh06NAD26alpeHZZ59FixYtYGRkhBkzZmg8Hl3wTdR1XM0sgJONOd58oqXU4RDpvToXdaNGjaq0LTQ0FNHR0Thw4EBdn55Ir93KKcaE9SfxxvY4FMjL0N63Mf58rTtGd/DmKis6Shsj+rdu3YoZM2Zg/vz5iI2NRffu3TFgwAAkJydX2V4ul8PZ2Rnz589HmzZtavWauu5GViG+2n8VAPD2oJZcz5hIA7TWK9vX1xdHjhyp9eMPHjyIwYMHw8PDAzKZDL/88kuF+4UQWLBgATw8PGBpaYlevXrhwoULdYyaSDNUKoEfjych4rODiLp8G2YmRnjziUBsebkzfBytpQ6PHmLNmjV44YUX0K1bNxw+fBiXL1/GiBEjEB4eDnPz2p1ZXbZsGV588UW89NJLaNmyJT7//HN4eXlh5cqVVbb39fXFF198gXHjxsHe3vBGggoh8Pav51FapkL3ACcMaeMhdUhEBkGjy4Tdr3HjxrV+bGFhIdq0aYPnn38eTz31VKX7P/roIyxbtgwbNmxA8+bN8f7776Nv3764dOmSeoULIilcv12AeTvi1H3n2nk3wsdPt0EzF5tHPJJ0hSZH9JeWliImJgZz586tsD0iIgLR0dGaChlyuRxyuVx9Oy8vT2PPrWl/nEvDoStZMDMxwntDQ3jWmkhDtFrU1cWAAQMwYMCAKu8TQuDzzz/H/PnzMXz4cADAxo0b4erqik2bNuGVV16pz1CJAAClZSp8c/AavvznKkrLVLA0Ncb/+rXA+C6+MOZEwnpD0yP6s7KyoFQq4erqWmG7q6urRvvuLVmyBAsXLtTY82lLfokC7/1RvhTYlF7N4OvEM9dEmqKXk2IlJiYiPT0dERER6m3m5ubo2bPnQ3/5yuVy5OXlVfgj0oSj17LxxJeH8MneyygtU6FHc2fsndkDL3TzY0GnZ/z9/XHo0CFs27YNMTEx2LFjByZPnowPP/ywTs97/9koIYRGz1DNmzcPubm56r+UlBSNPbcmfb7vCjLz5fB1tMIrPf2lDofIoOjsmbqHuffrtqpfvklJSQ98nL78kiX9kZlfgqW7L2JHbCoAwMnGDG8PCsKQNh68pKSn1q9fX2EA2L0R/YMGDUJSUhJWrFhRo+dzcnKCsbFxpbNymZmZlXJYXZibm9e6z199uZiehw3RNwAAC4eGcE46Ig3TyzN199T0l6++/JIl3adQqrDucCIe+yQKO2JTIZMBz3b0xt+zemFo2yYs6PSYpkf0m5mZISwsDJGRkRW2R0ZGokuXLrUNU+8IIfDOLxegVAkMCHFDz+bOUodEZHD08kydm5sbgPIzdu7u7urtj/rlqw+/ZEn3RV/NwoLfL+ByRgEAoLWnPRYNDUFbr0bSBkZaVZcR/bNmzcLYsWMRHh6Ozp0745tvvkFycjImTZoEoPwHZ2pqKr777jv1Y86cOQMAKCgowO3bt3HmzBmYmZkhKCiozvsihd/O3sKJG3dgYWqEtwbp5z4Q6Tq9LOr8/Pzg5uaGyMhI9XI9paWliIqKqnO/F6IHScouxAe7ErA3PgMA0NjKFLP7tcCo9t7sN9dA1HZE/8iRI5GdnY1FixYhLS0NISEh2L17N3x8fACUD864f866/y5FFhMTg02bNsHHxwc3btyodfxSKZCXYfHuBADA1N7N0KSRpcQRERkmnS3qCgoKcPXqVfXtxMREnDlzBg4ODvD29saMGTOwePFiBAQEICAgAIsXL4aVlRWeffZZCaMmQ5RbrMDy/Vex4cgNlCpVMDaS4bmO3pjZtzkaWZlJHR7picmTJz9wlZ0NGzZU2iaE0HJE9efrf64iI08OH0crvNSdgyOItEVni7pTp06hd+/e6tuzZs0CAIwfPx4bNmzAnDlzUFxcjMmTJ+Pu3bvo2LEj9u7dyznqSGNKy1T48XgSvvz7Cu4WKQAA3QOc8PagIDR35eeMqDoSswqx9vB1AMA7g4I4OIJIi2TCkH4O1lBeXh7s7e2Rm5sLOzs7qcMhHaFSCfwRl4ZP/rqE5DtFAIAAFxvMH9gSvVq4SBwd8XtbPbryPr208ST2JWSiZ3NnbHi+PQcRET1CXb67Onumjqi+CSFw6EoWPvrrIs6nls9h6Gxrjll9m+OZME+YGOv1YHGiehd1+Tb2JWTCxEiGtwcFsaAj0jIWdUQAYpLu4KM9l9RLe9mYm+DlHv54sZsfrM35NSGqKYVSpV45YnwXXy6TR1QP+L8VNWjnbuZgWeRlHLh0GwBgZmKE5zr6YErvpnC04fQ3RLW16XgyrmYWwMHaDNMfC5A6HKIGgUUdNUjnU3Px+b4r2JdQPj2JsZEMz4R5YvpjAfDgdAtEdZJbpMBn+y4DAGY+HgB7S1OJIyJqGFjUUYNy7mYOvvz7qrqYM5IBQ9s2wWuPBXBhcSIN+fKfK8gpUiDAxQajO3hLHQ5Rg8GijhqEmKQ7+Oqfq+rLrEYyYHAbD0x/LABNndnXh0hTbmQV4rujNwAA8we25AAjonrEoo4MlhACB69kYcX+q+oBEPfOzE3t04zFHJEWfLjnIhRKgR7NnTkFEFE9Y1FHBkepEvjzfBpWRV1TT01iaizD02GemNSzKXwceZmVSBtiku7gz/PpMJIB859oKXU4RA0OizoyGMWlSvwck4I1hxORlF0+abClqTFGd/DGxB5+cLfnAAgibRFC4P1d5eu7jgj3Qgs3rrpCVN9Y1JHeu50vx/dHb+D7Y0nq5bwaWZlifGdfjO/iCwdrrs9KpG1/nk9HbHIOLE2NMatvc6nDIWqQWNSR3rqYnod1hxPxy5lbKC1TAQC8HCzxUjd/PBPuCSszfryJ6oNCqcLHf10CAEzs4Q8XOwuJIyJqmPi/HukVpUrgn4uZWH8kEdHXstXb23k3wsTu/ogIcuVoO6J6tuVEMhKzCuFobYaXe/hLHQ5Rg8WijvRCbpECP51KwffHkpB8p7y/nJEMGBDijhe6+SHMp7HEERI1TIXyMnzx9xUAwPTHAmDDZfWIJMNvH+m0C7dy8f3RJPxyJhUlivJLrHYWJhjd0RvjOvuiCVd/IJLUmkOJyCoohY+jFScaJpIYizrSOSUKJf48n4YfjiUjJumuenugmy3Gd/HFsLZNYGlmLGGERAQAdwpL8e2h6wCA1yNawMyEXR+IpMSijnTGtdsF2Hw8GT+fvomcf0exmhjJ0C/EDeM7+6K9b2PIZDKJoySie1bsv4oCeRmCPewwqJW71OEQNXgs6khSJQol9pxPx+YTyepVHwDAw94Cozt4Y2QHL7jYciQdka5JzSnGd8eSAABz+gfCyIg/uIikxqKOJJGQloetJ1OwMzYVucXlZ+WMZEDvFi4Y08kbPZu7wJj/SRDprC/3XUFpmQqd/B3QI8BJ6nCICCzqqB7lFinw27lb2HYqBedu5qq3N2lkiWfCPTEi3AseHPhApPMSswrx8+mbAID/9QtktwgiHcGijrRKqRI4fDUL206lYG98hnqSYFNjGfoGuWJEuBe6BzjzrByRHvks8jKUKoHHAl04nRCRDmFRR1pxKT0fO07fxM7YVGTmy9XbW7jaYkR7Lwxr6wFHG3MJIySi2riYnoffz90CAMyK4HJgRLqERR1pTGZeCX47ews7Y1Nx4VaeentjK1MMaeOBZ8K9EOxhx0s1RHps2d7LEAIY2ModwR72UodDRP/Boo7qpEBehj3n0/HrmVQcuZoFlSjfbmosQ+8WLngqzBO9W7hw/ioiAxB3Mxd74zNgJANm9g2QOhwiug+LOqqxEoUSUZdv47czt7AvIQPyf/vJAUCYT2MMa+uBQa090NjaTMIoiUjTPt93GQAwtG0TNHOxlTgaIrofizqqFoVShehr2fj97C38dT4d+fIy9X3+ztYY1rYJhrb1gI+jtYRREpG2nE3Jwd8XM2EkA6b1aSZ1OERUBRZ19EBKlcDx69n4/Vwa9pxPw91/V3kAADc7Cwxq7Y5h7ZqwnxxRA/DZv2fpnmznCX9nG4mjIaKqsKijCsqUKpxIvINdcWnYcz4d2YWl6vscrc0woJUbBrf2QHtfB84gT9RAxCbfxYFLt2FsJMP0x3iWjkhXsagjlJapcPR6NvacT8NfFzJw5z+FXGMrU/QPccPAVh7o5O8AE2MOeCBqaL74+woA4Ml2TdjFgkiHsahroEoUShy8fBt7zqdjX0IG8kr+v49cYytT9At2wxOt3NG5qSNMWcgRNVhnUnLUZ+mm9uZZOiJdxqKuAcktUuCfSxn463wGoi7fRrFCqb7PycYc/YJdMSDEnWfkiEjty3/P0g1t6wFfJ56lI9JlLOoM3K2cYkTGZ2BvfDqOX7+DsnsTyaF8zdV+wW7oF+yKcF8HLtVFRBXE3czFP+oRr5yXjkjXsagzMEIIXLiVh8j4DOxLyKiwsgMANHe1QUSQG/qHuHHUKhE91Nf7y8/SDWnjAT+epSPSeSzqDECJQonoa1n4OyET/1zMRFpuifo+I1n5hMB9g1zRN8iNiZmIquVSej7+upABmQyYwr50RHqBRZ2eSsstxv6Lt/HPxQwcvpqFEsX/r+pgZWaM7gFOeLylK/oEusDRxlzCSIlIHy3ffxUAMCDEDQGuXD2CSB+wqNMTZUoVYlNysP9i+dm4i+n5Fe53t7fAYy1d8FhLV3T2d4SFqbFEkRKRvrt+uwB/nLsFgGfpiPQJizodlplfgqhLt3Hg8m0cuny7wrQjMhnQzqsR+gSWF3KBbrbsH0dEGrEq6hpUAngs0AXBHvZSh0NE1cSiTocolCqcTrqLqMu3EXX5dqVBDvaWpujZ3Bm9A53Rs7kLHKzNJIqUiAzVrZxi7DidCgCYwjVeifQKizqJJWcX4eCV2zh4+Tair2WjQF5W4f5WTezRq4UzerVwQVuvRpx2hIi06ttD11GmEujk74BQ78ZSh0NENcCirp7llyhw9Fo2Dl3JwsErt5GUXVThfgdrM3Rr5oReLZzRPcAZzrYc5EBE9SO7QI7NJ5IBsC8dkT5iUadlCqUKZ1JycPhKFg5fzcKZlBwo/zMBsImRDKHejdGjuRO6BzijVRN7GPFsHBFJYEP0DZQoVGjtaY9uzZykDoeIaohFnYYJIXA5owCHr2Yh+moWjl3PRmGpskIbPydrdGvmhB7NndHJ3wG2FqYSRUtEVK5QXobvjiYBACb3asqBV0R6SO8X+FyxYgX8/PxgYWGBsLAwHDp0qN5jSLlThK0nkzF9cyzaf/A3+n1+EO/9EY+/L2aisFQJB2szDGztjqXDW+HQnN7YP7sX3hsWgr5BrizoiBqImuaqqKgohIWFwcLCAv7+/li1apVW49t8Ihm5xQr4O1mjb5CbVl+LiLRDr8/Ubd26FTNmzMCKFSvQtWtXrF69GgMGDEB8fDy8vb219roZeSU4ei0b0deyEH0tGzfvFle438LUCO19HdCtmRO6NnNCkLsdL6kSNWA1zVWJiYl44oknMHHiRPzwww84cuQIJk+eDGdnZzz11FMaj6+0TIW1hxMBAC/38OeALCI9JRNCiEc3000dO3ZEaGgoVq5cqd7WsmVLDBs2DEuWLHnk4/Py8mBvb4/c3FzY2dk9sF1ukQIHr9zG0evZOHYtG9ezCivcb2IkQ1uvRujS1BFdmjmhnXcjmJtw8l8ibaju91aX1DRXvfHGG/jtt9+QkJCg3jZp0iScPXsWR48erdZr1uR9+jnmJmZvOwtnW3McfqM38xeRhOqS4/T2TF1paSliYmIwd+7cCtsjIiIQHR1d5WPkcjnkcrn6dl5eXpXt7nchLRfTNseqb8tkQIiHPTo3dUSXpo5o7+sAa3O9fSuJSItqk6uOHj2KiIiICtv69euHtWvXQqFQwNS0creN2uY3lUpgddQ1AMALXf1Y0BHpMb2tRLKysqBUKuHq6lphu6urK9LT06t8zJIlS7Bw4cIav1aod2O08bRHmI8DOjd1RAc/B9hbsi8cET1abXJVenp6le3LysqQlZUFd3f3So+pbX4rUwmMbO+FrSdTMKaT9rqtEJH26f1AiftHaAkhHjhqa968ecjNzVX/paSkVOs1LEyN8evUbnhncBD6BrmyoCOiGqtJrnpQ+6q231Pb/GZmYoSXuvtj78wesOPALSK9prdn6pycnGBsbFzpl25mZmalX7j3mJubw9yck/kSUf2pTa5yc3Orsr2JiQkcHR2rfExd8xunMCHSf3p7ps7MzAxhYWGIjIyssD0yMhJdunSRKCoioopqk6s6d+5cqf3evXsRHh5eZX86IiJAj4s6AJg1axbWrFmDdevWISEhATNnzkRycjImTZokdWhERGqPylXz5s3DuHHj1O0nTZqEpKQkzJo1CwkJCVi3bh3Wrl2L2bNnS7ULRKQH9PbyKwCMHDkS2dnZWLRoEdLS0hASEoLdu3fDx8dH6tCIiNQelavS0tKQnJysbu/n54fdu3dj5syZWL58OTw8PPDll19qZY46IjIcej1PXV3p43xXRA0dv7fVw/eJSD81yHnqNOFePVvd+ZyISHr3vq8N+PdotTC/EemnuuS4Bl3U5efnAwC8vLwkjoSIaio/Px/29vZSh6GzmN+I9FttclyDvvyqUqlw69Yt2Nra1utw/ry8PHh5eSElJUWvL4sYwn5wH3RHdfdDCIH8/Hx4eHjAyEivx3pplVT5DTCMzyT3QXcYwn7UZB/qkuMa9Jk6IyMjeHp6Svb6dnZ2evsB/S9D2A/ug+6ozn7wDN2jSZ3fAMP4THIfdIch7Ed196G2OY4/c4mIiIgMAIs6IiIiIgPAok4C5ubmePfdd/V+yTJD2A/ug+4wlP0gwziW3AfdYQj7UV/70KAHShAREREZCp6pIyIiIjIALOqIiIiIDACLOiIiIiIDwKKOiIiIyACwqNOCDz74AF26dIGVlRUaNWpUrccIIbBgwQJ4eHjA0tISvXr1woULFyq0kcvlmDZtGpycnGBtbY0hQ4bg5s2bWtiDcnfv3sXYsWNhb28Pe3t7jB07Fjk5OQ99jEwmq/Lv448/Vrfp1atXpftHjRqlM/swYcKESvF16tSpQpv6PBY13QeFQoE33ngDrVq1grW1NTw8PDBu3DjcunWrQjttH4cVK1bAz88PFhYWCAsLw6FDhx7aPioqCmFhYbCwsIC/vz9WrVpVqc327dsRFBQEc3NzBAUFYefOnRqLl6rPEHIc85tu5Lfa7Icu5DidzW+CNO6dd94Ry5YtE7NmzRL29vbVeszSpUuFra2t2L59u4iLixMjR44U7u7uIi8vT91m0qRJokmTJiIyMlKcPn1a9O7dW7Rp00aUlZVpZT/69+8vQkJCRHR0tIiOjhYhISFi0KBBD31MWlpahb9169YJmUwmrl27pm7Ts2dPMXHixArtcnJydGYfxo8fL/r3718hvuzs7Apt6vNY1HQfcnJyxOOPPy62bt0qLl68KI4ePSo6duwowsLCKrTT5nHYsmWLMDU1Fd9++62Ij48Xr732mrC2thZJSUlVtr9+/bqwsrISr732moiPjxfffvutMDU1FT///LO6TXR0tDA2NhaLFy8WCQkJYvHixcLExEQcO3ZMIzFT9RlCjmN+0438Vpv9kDrH6XJ+Y1GnRevXr69WwlOpVMLNzU0sXbpUva2kpETY29uLVatWCSHKP8SmpqZiy5Yt6japqanCyMhI7NmzR+Oxx8fHCwAVPlBHjx4VAMTFixer/TxDhw4Vffr0qbCtZ8+e4rXXXtNUqA9U230YP368GDp06APvr89joanjcOLECQGgQtLR5nHo0KGDmDRpUoVtgYGBYu7cuVW2nzNnjggMDKyw7ZVXXhGdOnVS3x4xYoTo379/hTb9+vUTo0aN0lDUVFP6muOY34Y+8H59/b+mPnOcLuc3Xn7VAYmJiUhPT0dERIR6m7m5OXr27Ino6GgAQExMDBQKRYU2Hh4eCAkJUbfRpKNHj8Le3h4dO3ZUb+vUqRPs7e2r/XoZGRnYtWsXXnzxxUr3/fjjj3ByckJwcDBmz56N/Px8jcV+T1324cCBA3BxcUHz5s0xceJEZGZmqu+rz2OhieMAALm5uZDJZJUulWnjOJSWliImJqbC+wMAERERD4z56NGjldr369cPp06dgkKheGgbbXz+SbN0Lccxv+lGfqvrfvxXfeU4Xc9vJjVqTVqRnp4OAHB1da2w3dXVFUlJSeo2ZmZmaNy4caU29x6v6ZhcXFwqbXdxcan2623cuBG2trYYPnx4he1jxoyBn58f3NzccP78ecybNw9nz55FZGSkRmK/p7b7MGDAADzzzDPw8fFBYmIi3n77bfTp0wcxMTEwNzev12OhieNQUlKCuXPn4tlnn62wkLS2jkNWVhaUSmWVn+cHxZyenl5l+7KyMmRlZcHd3f2BbbTx+SfN0rUcx/ymG/mtLvvxX/WZ43Q9v/FMXTUtWLDggZ1k7/2dOnWqTq8hk8kq3BZCVNp2v+q0+a+a7EdVz1uT11u3bh3GjBkDCwuLCtsnTpyIxx9/HCEhIRg1ahR+/vln7Nu3D6dPn9aJfRg5ciQGDhyIkJAQDB48GH/++ScuX76MXbt2PTSumrw39XUcFAoFRo0aBZVKhRUrVlS4r67H4VFq+nmuqv3922vzHaHqMYQcx/ymG/mtPvbjHqlynK7mN56pq6apU6c+ctSMr69vrZ7bzc0NQHk17+7urt6emZmprtzd3NxQWlqKu3fvVvgFlZmZiS5dulT7taq7H+fOnUNGRkal+27fvl3p10RVDh06hEuXLmHr1q2PbBsaGgpTU1NcuXIFoaGhj2xfX/twj7u7O3x8fHDlyhUAmjkW9bEPCoUCI0aMQGJiIv75558Kv2CrUtPj8CBOTk4wNjau9Avzv5/n+7m5uVXZ3sTEBI6Ojg9tU5NjSQ9mCDmO+U038htguDlO5/NbjXrgUY3UtBPxhx9+qN4ml8ur7ES8detWdZtbt25pvfPq8ePH1duOHTtW7c6r48ePrzQS6UHi4uIEABEVFVXreKtS1324JysrS5ibm4uNGzcKIer3WNR2H0pLS8WwYcNEcHCwyMzMrNZrafI4dOjQQbz66qsVtrVs2fKhHYlbtmxZYdukSZMqdSQeMGBAhTb9+/fnQAkJ6WuOY377f1LmNyH0M8fpcn5jUacFSUlJIjY2VixcuFDY2NiI2NhYERsbK/Lz89VtWrRoIXbs2KG+vXTpUmFvby927Ngh4uLixOjRo6sc7u/p6Sn27dsnTp8+Lfr06aP1YeatW7cWR48eFUePHhWtWrWqNMz8/v0QQojc3FxhZWUlVq5cWek5r169KhYuXChOnjwpEhMTxa5du0RgYKBo166d1qYtqMk+5Ofni9dff11ER0eLxMREsX//ftG5c2fRpEkTyY5FTfdBoVCIIUOGCE9PT3HmzJkKw/nlcrkQQvvH4d6Q/7Vr14r4+HgxY8YMYW1tLW7cuCGEEGLu3Lli7Nix6vb3hvzPnDlTxMfHi7Vr11Ya8n/kyBFhbGwsli5dKhISEsTSpUs5pYlEDCHHMb/pRn6rzX5IneN0Ob+xqNOC8ePHCwCV/vbv369uA0CsX79efVulUol3331XuLm5CXNzc9GjRw8RFxdX4XmLi4vF1KlThYODg7C0tBSDBg0SycnJWtuP7OxsMWbMGGFraytsbW3FmDFjxN27dyu0uX8/hBBi9erVwtLSssr5gJKTk0WPHj2Eg4ODMDMzE02bNhXTp0+vNE+SVPtQVFQkIiIihLOzszA1NRXe3t5i/Pjxld7n+jwWNd2HxMTEKj9///0M1sdxWL58ufDx8RFmZmYiNDS0wq/j8ePHi549e1Zof+DAAdGuXTthZmYmfH19q/xPc9u2baJFixbC1NRUBAYGiu3bt2ssXqo+Q8hxzG+6kd9qsx+6kON0Nb/JhPi3tx4RERER6S2OfiUiIiIyACzqiIiIiAwAizoiIiIiA8CijoiIiMgAsKgjIiIiMgAs6oiIiIgMAIs6IiIiIgPAoo6IiIjIALCoI4PSq1cvzJgxQ+owiIg0jvmNHoVFHREREZEBYFFHREREZABY1JHBUalUmDNnDhwcHODm5oYFCxZIHRIRkUYwv9HDsKgjg7Nx40ZYW1vj+PHj+Oijj7Bo0SJERkZKHRYRUZ0xv9HDyIQQQuogiDSlV69eUCqVOHTokHpbhw4d0KdPHyxdulTCyIiI6ob5jR6FZ+rI4LRu3brCbXd3d2RmZkoUDRGR5jC/0cOwqCODY2pqWuG2TCaDSqWSKBoiIs1hfqOHYVFHREREZABY1BEREREZABZ1RERERAaAo1+JiIiIDADP1BEREREZABZ1RERERAaARR0RERGRAWBRR0RERGQAWNQRERERGQAWdUREREQGgEUdERERkQFgUUdERERkAFjUERERERkAFnVEREREBoBFHREREZEBYFFHREREZABY1BEREREZABZ1RERERAaARR0RERGRATCROgApqVQq3Lp1C7a2tpDJZFKHQ0TVIIRAfn4+PDw8YGTE36UPwvxGpJ/qkuMadFF369YteHl5SR0GEdVCSkoKPD09pQ5DZzG/Eem32uS4Bl3U2draAih/4+zs7CSOhoiqIy8vD15eXurvL1WN+Y1IP9UlxzXoou7eJQk7OzsmPSI9w0uKD8f8RqTfapPj2CGFiIiIyACwqCMiIiIyACzqiIiIiAwAizoiklyhvAzL919FiUIpdShERBqhVAncyCpETNId7L+YiT3n07X+mg16oAQR6YZP9l7C+iM3cOx6Nr5/saPU4RAR1ZhCqcLJxDs4eCULR69n41J6HkoUKvX9NuYm6B/iptUYWNQRkaTOpORgQ/QNAMDE7v7SBkNEVANCCMSm5GB7zE3sjkvD3SJFhfvNTYzgYmcOe0tT2FuaQqUSMDLS3sh9FnVEJBmFUoW5289BCODJdk3Qo7mz1CERET1SaZkKv55JxcajN3A+NU+93cnGDD2aO6NHgDNae9rDx9Eaxlos4u7Hoo6IJPPtoeu4mJ6PxlameGtgS6nDISJ6qBKFEptPJOPbg9dxK7cEAGBmYoRBrdzxZGgTdPZ3hImxdMMVWNQRkSQupefj88grAIC3BgbB0cZc4oiIiKpWplTh55ib+OLvK0j7t5hztjXH8119Maq9NxyszSSOsByLOiKqdwqlCq9vO4NSpQqPBbpgeGgTqUMiIqpS1OXb+GBXPC5nFAAA3OwsMKVPMzwT5gkLU2OJo6uIRR0R1bsV+6/hfGoe7C1NsWR4Ky75RUQ6Jzm7CAt/v4C/L2YCABpZmWJq72Z4rpOPzhVz97CoI6J6FZt8F1/+U37ZddHQYLjYWUgcERHR/5OXKbHywDWsOHANpWUqmBrLML6zL6b1CYC9lanU4T0Uizoiqjf5JQq8tuUMlCqBwW08MKSNh9QhERGpHb+ejXk743D9diEAoFszJywYEoxmLjYSR1Y9LOqIqN68++sFJN8pQpNGlvjgyRBediUinVAgL8PSPxPww7FkAICTjTneHRyEQa3d9SpPsagjonrx06kU7IhNhZEM+GJUW9hZ6PZlDCJqGI5czcKcn88hNacYADC6gxfm9m+p85daq8Kijoi07mJ6Ht759TwAYObjzRHu6yBxRETU0BWVlmHpnxfx3dEkAICXgyWWDm+Nrs2cJI6s9ljUEZFW5ZcoMPmH0yhRqNCzuTOm9G4mdUhE1MCdScnBrK1ncD2rvO/c2E4+mDsgENbm+l0W6Xf0RKTTVCqB1386i+tZhXC3t8BnI9tqdd1DIqKHUaoEVuy/is//vgKlSsDNzgIfP9Ma3QMMY4lCFnVEpDXL91/F3vgMmBkbYcWYUJ2ZdZ2IGp6bd4swc+sZnLxxFwAwuI0H3h8aopd95x6ERR0RacXfCRlYtu8yAOC9YcFo591Y4oiIqKHaHZeGN7afQ35JGWzMTfDesGA82c5T6rA0jkUdEWlcQloepm+OhRDAmI7eGNneW+qQiKgBKlEo8d4f8fjxePlUJW28GuHLUW3h42gtcWTawaKOiDTqdr4cL208hcJSJTr7O2LBkGCpQyKiBuhqZgGmbjqNi+n5kMmAST2bYlbf5jA1NpI6NK1hUUdEGlNcqsTE704hNacYfk7WWPlcqEEnUCLSTb/EpuLNnXEoKlXCycYMn41sazCDIR6GRR0RaYRSJTB9SyzOpOSgkZUp1owPRyMrDowgovpTolBi4e/x2Hyi/HJrZ39HfDGqbYNZY5pFHRHVmRACC367gMj4DJiZGGHNuHA0ddaPtRKJyDAkZxfh1R9jcOFWHmQyYFqfALz2WACMG9A0SizqiKjOvvj7Cr4/lgSZDPh8ZFuuGEFE9SoyPgOzfjqD/JIyOFib4fORbdGjueFfbr0fizoiqpPvj97A5/uuAAAWDgnGE63cJY6IiBqKMqUKyyIvY8WBawCAUO9GWD4mFO72lhJHJg0WdURUa9tjbuLtXy8AAF57LADjOvtKGxARNRjZBXJM3xKLI1ezAQAvdPXD3AGBMDNpuIOzWNQRUa3sjkvD/34+CwCY0MUXMx4PkDgiImoozqTk4NUfYpCWWwIrM2N8+FRrDG7jIXVYkmNRR0Q1tvdCOqZvjoVKACPDvfDOoCDIZA2nMzIRSWfziWS8++sFlCpV8HeyxuqxYQhwtZU6LJ3Aoo6IamTvhXRM2XQaZSqBoW09sHh4Kxg1oNFlRCQNeZkSC367gM0nUgAAEUGu+HREG9haGM7arXXFoo6Iqu2vC+mYuuk0FEqBwW088OkzbRrUdAFEJI303BJM+iEGZ1JyIJMBsyNaYHKvprxCcB8WdURULb+dvYWZW89AqSov6D4b0QYmXC2CiLTs5I07ePWH08gqkMPe0hRfjm6Hng1wupLqYFFHRI+07VQK3th+DioBDG/XBB893ZoFHRFplRACPx5PxoLfLqBMJRDoZovVY8Pg42gtdWg6i0UdET3U2sOJeO+PeADAqPZeWPwk+9ARkXbd339uYGt3fPx0a1iZsWx5GL47RFQlIQQ+i7yML/+5CgB4qZsf5g9syT4sRKRVmfklePWH04hJuguZDJjTLxCTevoz91QDizoiqqRMqcLbv55X/0p+vW9zTO3TjEmViLTqbEoOXvk+Bul5JbC1MMFXo9uhVwsXqcPSGyzqiKiCotIyTN98BvsSMmAkA94bFoIxHX2kDouIDNzO2Jt4Y3scSstUaOpsjW/HhcPf2UbqsPQKizoiUrudL8dLG0/i7M1cmJsY4cvR7dAv2E3qsIjIgClVAh/uuYhvDl4HADze0gWfjWzL+edqgUUdEQEALmfk44UNJ3HzbjEaW5lizfhwhPk4SB0WERmw3GIFpm+ORdTl2wCAqb2bYVbf5hyMVUss6ogIBy5lYtqmWOTLy+DjaIUNz3eAnxOnDSAi7bl2uwATN57C9axCWJga4ZNn2mBQa67fWhcs6ogaMCEE1h25gQ92xUMlgA5+Dlj9XBgaW5tJHRoRGbCoy7cxddNp5JeUwcPeAt+MC0dIE3upw9J7LOqIGih5mRJv7TyPbTE3AQDPhHnigydbwcyEkwoTkXbc/0MyzKcxVj0XBmdbc6lDMwgs6ogaoP+uo2gkA+YPDMILXX05ZQkRaY28TIm3fzmPn06V/5AcEe6J94aFwNzEWOLIDIfB/CRfsmQJZDIZZsyYIXUoRDrtROIdDPrqMM6k5MDe0hQbnu+AF7v5saDTooMHD2Lw4MHw8PCATCbDL7/88sjHREVFISwsDBYWFvD398eqVau0HyiRlmQVyPHcmuP46dRNGMmAtwcF4cOnWrOg0zCDKOpOnjyJb775Bq1bt5Y6FCKdJYTAmkPXMfrbY8gqkCPQzRa/T+2GHlwYW+sKCwvRpk0bfP3119Vqn5iYiCeeeALdu3dHbGws3nzzTUyfPh3bt2/XcqREmpeQloehXx/ByRt3YWthgvX8Iak1en/5taCgAGPGjMG3336L999/X+pwiHRSXokCb/x8Dn+eTwcADG3rgSXDW3EdxXoyYMAADBgwoNrtV61aBW9vb3z++ecAgJYtW+LUqVP45JNP8NRTT1X5GLlcDrlcrr6dl5dXp5iJNCEyPgOvbYlFUakSvo5WWDO+PZq5cEJhbdH7M3VTpkzBwIED8fjjjz+yrVwuR15eXoU/IkN3PjUXg786jD/Pp8PUWIZFQ4Px+ci2LOh02NGjRxEREVFhW79+/XDq1CkoFIoqH7NkyRLY29ur/7y8vOojVKIqCSGw4sBVvPz9KRSVKtGlqSN+mdKVBZ2W6XVRt2XLFpw+fRpLliypVnsmPWpIhBDYcCQRw1dEIym7CE0aWWLbpC4Y15kDInRdeno6XF1dK2xzdXVFWVkZsrKyqnzMvHnzkJubq/5LSUmpj1CJKilRKPH6T2fx0Z5LEAJ4rpM3Nr7QAY2sOFWStuntT/WUlBS89tpr2Lt3LywsLKr1mHnz5mHWrFnq23l5eSzsyCDdKSzFnJ/PYV9CBgDg8Zau+OSZ1kyqeuT+wlsIUeX2e8zNzWFuzmkhSFq38+V45ftTOJ2cA2MjGRYMDsLYzr5Sh9Vg6G1RFxMTg8zMTISFham3KZVKHDx4EF9//TXkcjmMjSuOqmHSo4bgyNUszNx6Bpn5cpgZG2H+wJYY19mHZ+f0iJubG9LT0ytsy8zMhImJCRwdHSWKiujhEtLy8NLGU0jNKYadhQlWjAlDtwAnqcNqUPS2qHvssccQFxdXYdvzzz+PwMBAvPHGG5UKOiJDJy9T4pO/LmHN4UQIATR1tsYXo9pxlnY91LlzZ/z+++8Vtu3duxfh4eEwNeUi56R79v07IKKwVAk/J2usGR+Ops7sP1ff9Laos7W1RUhISIVt1tbWcHR0rLSdyNAlpOVh5tYzuJieDwB4tqM33h4YBEsz/rjRBQUFBbh69ar6dmJiIs6cOQMHBwd4e3tj3rx5SE1NxXfffQcAmDRpEr7++mvMmjULEydOxNGjR7F27Vps3rxZql0gqlL5VEmJWPxnAoQAujR1xIoxoezqIRG9LeqICFCqBFYfvIbPIi9DoRRwtDbDh0+1xuNBro9+MNWbU6dOoXfv3urb9/r2jh8/Hhs2bEBaWhqSk5PV9/v5+WH37t2YOXMmli9fDg8PD3z55ZcPnM6ESAqlZSq89UuceoWI0R28sWhoMEyN9XoMpl6TiXu9bxugvLw82NvbIzc3F3Z2dlKHQ1Qj124X4H/bzuJ0cg6A8sEQS4a3Mvg1FPm9rR6+T6RNdwtLMemHGBxPvKNeIWJCF46s14S6fHd5po5IzyhVAusOJ+KTvZcgL1PB1twE7wwOwtNhnkyoRKR1VzML8OLGk0jKLoKNuQm+erYderdwkTosAos6Ir1yJSMf//v5HM6k5AAAugc44cOnWsOjkaW0gRFRg3D4ShZe/TEG+SVl8GxsibXj26OFm63UYdG/WNQR6YHSMhVWRV3D1/9cRamy/OzcmwNbYlR7L56dI6J68ePxJLzz6wUoVQJhPo2xemwYnGwMu7uHvmFRR6TjTiffxbztcbiUUT6ytU+gCz54MgTu9jw7R0Tap1QJfLArAeuOJAIAnmzXBEuGt4KFKUfX6xoWdUQ6Kq9EgU/+uoTvjyVBCMDB2gzvDg7CkDYePDtHRPWiQF6G6Ztj8c/FTADA7IjmmNK7GXOQjmJRR6RjhBDYFZeGRb/HIzNfDgB4KtQT8we2hIM1534iovpx824RXtp4ChfT82FuYoTPRrbFE63cpQ6LHoJFHZEOuZFViHd/u4Coy7cBAH5O1vhgWAi6NONSO0RUf2KT72LidzHIKpDD2dYca8aFo41XI6nDokdgUUekA0oUSqw8cA0ro66htEwFM2MjvNqrKV7t1ZT9VoioXv1x7hZe/+ks5GUqtHS3w9rx4RxhrydY1BFJSAiBfQmZWPTHBaTcKQZQPk3JoqEh8HOyljg6ImpIhBD4+p+r+DTyMgDgsUAXfDm6HazNWSroCx4pIolcv12Ahb/Hqy+1uttb4K2BQXiilRs7IRNRvZKXKTFvexx2xKYCAF7s5oc3n2gJYyPmIn3Coo6onuWXKPD1P1ex7kgiFEoBU2MZXuruj2l9msHKjF9JIqpfdwpLMen7GJy4cQfGRjIsGhqMMR19pA6LaoH/gxDVE6VK4OeYFHz81yVkFZQCAHq3cMY7g4N5qZWIJPHfJb9szU2w4rlQdA9wljosqiUWdUT14Nj1bLz3Rzwu3MoDAPg7WePtQUHoHcj1EolIGtFXszDphxjklZTBy8ES68a3R4Arl/zSZyzqiLToRlYhlvyZgL8uZAAAbC1M8NpjARjX2RdmJkYSR0dEDdVPJ1Pw5s44lKkEQr0b4Ztx4VzyywCwqCPSgruFpfjynyv44VgSFEoBYyMZRnfwwszHm8ORiZOIJKJSCXz01yWsiroGABjcxgMfP92aUycZCBZ1RBpUolBiY/QNLN9/FXklZQCAXi2cMf+JlrysQUSSKi5VYtZPZ/Dn+XQAwPQ+zTCzb3OOtjcgLOqINECpEvglNhXLIi8jNad8vrlAN1u8NTAI3QK4GgQRSSszrwQvfXcK527mwszYCB8+3QpPtvOUOizSMBZ1RHUghMCBS7fx4Z6LuJieD6B8vrnXI1rgyXZNOMcTEUnuYnoeXlh/ErdyS9DYyhSrx4ajg5+D1GGRFrCoI6ql08l38eGfF3E88Q6A8kEQk3s1w/Ndfdk/hYh0wv5LmZi2KRYF8jL4O1lj3YT28OUUSgaLRR1RDV3OyMcnf13C3vjyEa1mJkYY39kHk3s1Q2NrM4mjIyIq9/3RG3j3twtQCaCTvwNWPReGRlbMUYaMRR1RNaXcKcJnkZex80wqhACMZMBToZ6Y0bc5mnCxayLSEUqVwPu74rH+yA0AwNNhnlj8ZCtOo9QAsKgjeoT03BJ89c8VbD2ZgjKVAAD0D3bD7H7N0cyFI1qJSHcUysvw2pZY7EvIBAD8r18LTO7VlCNcGwgWdUQPkFUgx6oD1/D9sSTIy1QAgO4BTvhfvxZo7dlI2uCIiO6TlluMFzecQnxaHsxNjPDpiDYY1NpD6rCoHrGoI7rP3cJSfHvoOjZE30BRqRIA0N63MWZHtEBHf0eJoyMiqux8ai5e3HgSGXlyONmY4Ztx4Qj1bix1WFTPWNQR/Su3WIG1hxOx7nAiCuTlEwe38bTH6xEt0D3AiZcviEgn7YvPwPQtsSgqVSLAxQbrJrSHl4OV1GGRBFjUUYOXV6LA+sM3sObwdeT/uwpES3c7zOrbHI+3dGExR0Q6a/2RRLz3RzxUAujWzAnLx4TC3tJU6rBIIizqqMHKK1Fgw5EbWHPounpJr+auNpjxeHP0D3aDEScOJiIdVaZU4f1dCdgQfQMAMKq9F94bFgJTY45wbchY1FGDk1tcXsytPfz/xVyAiw2mPxaAga3cWcwRXnjhhWq1W7dunZYjIaqsQF6GaZtOY/+l2wCAeQMC8XIPf15VoPov6hQKBSIiIrB69Wo0b968vl+eGrDcIgXWHknE+iOJ6susAS42mPZvMcclveieDRs2wMfHB+3atYMQQupwiNTScovxwoZTSPh3hOvnI9tiQCt3qcMiHVHvRZ2pqSnOnz/PXxRUb+4WlmLt4URsiL6hHgDR3NUG0/rwzBxVbdKkSdiyZQuuX7+OF154Ac899xwcHLhWJkmr4ghXc6wZH462Xo2kDot0iCQX38eNG4e1a9dK8dLUgGQVyLHkzwR0+/AffL3/KgrkZQh0s8XyZ0Ox57UeGNzGgwUdVWnFihVIS0vDG2+8gd9//x1eXl4YMWIE/vrrL565I0nsi8/AM6uOIiNPjuauNtg5uQsLOqpEkj51paWlWLNmDSIjIxEeHg5r64qLCy9btkyKsMhAZOSVYHXUdWw6kYQSRfmkwUHudpj+WDNEBHEABFWPubk5Ro8ejdGjRyMpKQkbNmzA5MmToVAoEB8fDxsbG6lDpAZACIH1R27gvV3xEKJ8AvTlY0JhZ8ERrlSZJEXd+fPnERoaCgC4fPlyhft4WZZq6+bdIqyKuoafTt5EqbK8mGvtaY/pfQLwGKcmoTqQyWSQyWQQQkClUkkdDjUQZUoV3vsjHhuPJgEARnfwwqKhHOFKDyZJUbd//34pXpYMVGJWIVYeuIodp1PVa7OG+zTGtMcC0IOTBlMtyeVy7NixA+vWrcPhw4cxaNAgfP311+jfvz+MjPifKmlXgbwMUzedxoFLtyGTlY9wndidI1zp4TilCemtS+n5WL7/Kv44dwv/1nLo0tQR0/oEoJO/A5Mf1drkyZOxZcsWeHt74/nnn8eWLVvg6Mgl4qh+3MopxgsbTuJiej4sTMtHuPYP4QhXejQWdaR3zqbkYPn+q9gbn6He1ifQBVP7NONah6QRq1atgre3N/z8/BAVFYWoqKgq2+3YsaOeIyNDdz41Fy9sOInM/PIRrmvHh6MNB0RQNbGoI70ghMDxxDtYvv8qDl3JAgDIZMATIe6Y3Lspgj3sJY6QDMm4ceN4ppfqXWR8BqZvjkWxQokWrrZYOyEcno25hitVH4s60mlCCBy4dBvL91/FqaS7AABjIxmGtvXA5F7N0MyFIxBJ8zZs2CB1CNSACCGw7sgNvM8RrlRHLOpIJylVAnvOp2P5/quIT8sDAJiZGGFEuCde6dEUXg789UpE+q9MqcKiP+Lx3b8jXJ/t6I2FQ4I5wpVqhUUd6RSFUoVfYlOxMuoart8uBABYmRnjuU4+eKmbH1zsLCSOkIhIMzjClTSNRR3phBKFEj+dSsHqqOtIzSkGANhbmmJCF19M6OKLxtZmEkdIRKQ5lUe4tkP/EDepwyI9x6KOJJVfosAPx5Kx9vB1ZBWUAgCcbMwxsbsfxnTygY05P6JEZFjibpav4XpvhCvXcCVN4f+YJIk7haXYcCQRG6JvIK+kDADg2dgSr/RsimfCPGFhaixxhEREmvfXhXTM2HJGPcJ13fPt0aSRpdRhkYFgUUf1KiOvBN8evI4fjyejWKEEADR1tsbkXs0wpK0HOwcTkUESQmDNoUQs/jMBQgA9mjtj+bPtYMsRrqRBLOqoXiRnF2Fl1DVsj/n/dVlDmthhSq9m6BfsBiMjdgwmIsOkUKrw7m8XsOl4MgBgzL8jXE34I5Y0jEUdadWVjHysOHANv529BeW/a3l18HXA5N5N0bO5M0d5EZFByytRYMqPp3HoShZkMmD+Ey3xYjc/5j7SCv5MIK2Iu5mLSd/HoO9nB7EzNhVKlUDP5s746ZXO+GlSZ/Rq4cKkRg3KihUr4OfnBwsLC4SFheHQoUMPbHvgwAHIZLJKfxcvXqzHiKmuUu4U4emV0Th0JQuWpsZY9VwYXuKUJaRFenumbsmSJdixYwcuXrwIS0tLdOnSBR9++CFatGghdWgNWkzSHXz1z1UcuHRbva1/sBum9G6GVp5cyosapq1bt2LGjBlYsWIFunbtitWrV2PAgAGIj4+Ht7f3Ax936dIl2NnZqW87OzvXR7ikAWdScvDSxpPIKiiFi6051o5vzxxIWqe3RV1UVBSmTJmC9u3bo6ysDPPnz0dERATi4+NhbW0tdXgNihAC0dey8dU/V3Ds+h0AgJEMGNLGA5N7N0NzV1uJIySS1rJly/Diiy/ipZdeAgB8/vnn+Ouvv7By5UosWbLkgY9zcXFBo0aN6ilK0pTdcWmYufUM5GUqtHS3w7oJ4XC35whX0j69Ler27NlT4fb69evh4uKCmJgY9OjRQ6KoGpZ767J+9c8VnE7OAQCYGsvwVKgnJvVsCl8nFtdEpaWliImJwdy5cytsj4iIQHR09EMf265dO5SUlCAoKAhvvfUWevfu/cC2crkccrlcfTsvL69ugVONCSGwMuoaPtpzCQDQJ9AFX45ux/k2qd4YzCctNzcXAODg4PDANkx6mqFSCUQmZODrf64iLrX8fTczMcLo9l54pWdTeHDOJSK1rKwsKJVKuLq6Vtju6uqK9PT0Kh/j7u6Ob775BmFhYZDL5fj+++/x2GOP4cCBAw/80bpkyRIsXLhQ4/FT9ZSWqfDmzjj8HHMTADChiy/eHhQEY47sp3pkEEWdEAKzZs1Ct27dEBIS8sB2THp1o1IJ7LmQji//voKL6fkAAEtTY4zp6I2Xe/hzXVaih7i/c7wQ4oEd5lu0aFGhf3Dnzp2RkpKCTz755IFF3bx58zBr1iz17by8PHh5eWkgcnqUnKJSTPohBseu34GRDHh3cDDGd/GVOixqgAyiqJs6dSrOnTuHw4cPP7Qdk17tqFQCu8+n4cu/r+ByRgEAwMbcBOM6++DFbn5wtDGXOEIi3eXk5ARjY+NKZ+UyMzMrnb17mE6dOuGHH3544P3m5uYwN+d3sb4lZhXihQ0nkZhVCBtzE3z1bDv0buEidVjUQOl9UTdt2jT89ttvOHjwIDw9PR/alkmvZpQqgd1x5cXclczyYs7W3ATPd/XFC9380MjKTOIIiXSfmZkZwsLCEBkZiSeffFK9PTIyEkOHDq3288TGxsLd3V0bIVItHbuejUk/xCCnSIEmjSyxdkI4At3sHv1AIi3R26JOCIFp06Zh586dOHDgAPz8/KQOyWCoVAK77i/mLEzwYjc/PN/VD/aWXNaGqCZmzZqFsWPHIjw8HJ07d8Y333yD5ORkTJo0CUD5VYTU1FR89913AMpHx/r6+iI4OBilpaX44YcfsH37dmzfvl3K3aD/2HYqBW/ujINCKdDGqxG+HRcGF1t2QSFp6W1RN2XKFGzatAm//vorbG1t1Zc27O3tYWnJjvq1ca/P3Of7Lqsvs9pZmOCl7v6Y0NUXdlyjkKhWRo4ciezsbCxatAhpaWkICQnB7t274ePjAwBIS0tDcnKyun1paSlmz56N1NRUWFpaIjg4GLt27cITTzwh1S7Qv1QqgY/3XsLKA9cAAANbu+PTZ9rAwtRY4siIAJkQQkgdRG08qIPx+vXrMWHChGo9R15eHuzt7ZGbm1thgs+GRgiByPgMLIu8rB4AYWthgpe6+eP5bizmSLfwe1s9fJ80r7hUiVk/ncGf58tPIkzr0wwzH2/OtatJo+ry3dXbM3V6WovqFCEEoi7fxrLIyzh3s3xqEltzE7zQzQ8vdONlViKiezLySvDSxlOIS82FmbERlj7VCsNDH96Pm6i+6W1RR3VzIvEOPv7rIk7euAsAsDIzxoQuvni5hz8HQBAR/cf51Fy8tPEU0vNK4GBthtVjw9De98FzohJJhUVdA3M+NRcf/3UJUZfL12Y1MzHCuE4+mNSrKZw4NQkRUQV7zqdj5tYzKFYoEeBig7Xj28Pb0UrqsIiqxKKugbh+uwCfRl7GrnNpAAATIxlGtPfC9D4BcLPniC0iov+6f8mv7gFOWD4mlH2MSaexqDNwmXkl+OLvK9hyMgVKlYBMBgxt44GZfZvDx5FrsxIR3U9epsS8HXHYcToVADC+sw/eHhQEE2MjiSMjejgWdQaqQF6Gbw5ex7cHr6NYoQRQvrj0//q1QEt3joQjIqpKVoEck76PwamkuzA2kmHB4CCM7ewrdVhE1cKizsCUKVXYcjIFn++7jKyCUgBAO+9GmNs/EB39HSWOjohIdyWk5eGljaeQmlMMWwsTrBgTiu4BzlKHRVRtLOoMhBACBy7dxge7E3D131Ug/JysMadfC/QPcXvgvH5ERATsi8/Aa1tiUViqhK+jFdaMb49mLjZSh0VUIyzqDMDljHy890c8Dl3JAgA0tjLFjMeb49mO3jBlHxAiogcSQmD1wev4cM9FCAF0aeqIFWNCObUT6SUWdXosp6gUn0Vexg/Hk6FUCZgZG2FCV19M6d2MEwcTET1CiUKJN3fEYUds+YCIMR29sWBIMH8Mk95iUaeHlCqBTSeS8eneS8gpUgAA+ge7Yd4TgRzRSkRUDZn5JXjl+xjEJufA2EiGdwYFYXwXX6nDIqoTFnV6JibpDt7+5QLi0/IAAIFutnhncBC6NHWSODIiIv1wPjUXE787hbTcEthZmGDFmDB0C2AOJf3Hok5P3CksxZLdCdgWcxMAYGdhgtn9WuDZDt6cO4mIqJr+OHcLs7edRYlCBX9na6wZFw5/Zw6IIMPAok7HqVQCP8fcxOI/E9SXWkeEe+KN/oFw5LJeRETVolIJLIu8jK/3XwUA9GzujK+ebccVIsigsKjTYdduF2DejjicSLwDoPxS6wdPhiDMhwtJExFVV36JAjO3nsW+hAwAwMs9/PFG/0AYG3GqJzIsLOp0kEKpwjcHr+OLv6+gtEwFS1NjzOrbHM939eWlViKiGriRVYiJ353ClcwCmJkYYenwVhge6il1WERawaJOx8TfysOc7WdxPrV8IESP5s74YFgIvBysJI6MiEi/RF2+jembY5FbrICrnTm+GRuONl6NpA6LSGtY1OmIMqUKq6Ku4Yu/r0ChFLC3NMU7g4IwPLQJV4MgIqoBIQS++XdCYZUoXypx9XNhcLGzkDo0Iq1iUacDErMKMWPrGZxNyQEA9A1yxQfDQpiAiIhqqKi0DG9sj8PvZ28BAEaGe2HRsGCYmxhLHBmR9rGok5AQAltPpmDRH/EoKlXC1sIEC4cE48l2PDtHRFRTKXeK8PL3MUhIy4OJkQzvDA7C2E4+zKfUYLCok0husQJv7ojDrrg0AEBnf0d8OqINPBpZShwZEZH+OXj5Nqb923/OycYMK8aEoYMfZwqghoVFnQTO3czB5B9P4+bdYpgYyTC7Xwu83N0fRhxeT0RUI0IIrDhwDZ/uvQSVANp4NcKq50Lhbs8fyNTwsKirR0KUr9m68Ld4lCpV8Gxsia9Gt0M778ZSh0ZEpHfySxR4/aez2BtfPv/cqPZeWDiU/eeo4WJRV0/kZUq8/ct5/HSqfJmvx1u64tMRbWBvydnMiYhq6nJGPiZ9H4PrWYUwMzbCwqHBGN3BW+qwiCTFoq4eZOaXYNL3MTidnAMjGTC7XwtM6tGUl1uJiGrh1zOpmLs9DsUKJdztLbDyuTC05fxzRCzqtO1yRj6eX38SqTnFsLMwwVfPhqJnc2epwyIi0julZSp8sCseG48mAQC6NXPCF6Pach1son+xqNOi6GtZeOX7GOSXlMHfyRprxofD39lG6rCIiPROak4xpvx4Gmf+nc9zSu+mmNW3BddvJfoPFnVasud8OqZvjkWpUoX2vo3xzdhwNLY2kzosIiK988/FDMz66SxyihSwszDBZyPb4rGWrlKHRaRzWNRpwfaYm5iz/RyUKoEBIW74bGRbWJhyNBYRUU0olCp8svcSVkddBwC09rTH8mdDuRY20QOwqNOwn2NuYva2swCAp8M8sXR4K5gYG0kcFRGRfknNKcb0zbGISboLABjX2QfzB7bkdCVED8GiToN+P3sLc34uL+jGdfbBgsHBHOFKRFRDf11Ix5yfzyG3WAFbcxN8+HRrPNHKXeqwiHQeizoN2X8xEzO2noFKAKM7eGHhkGCuN0hEVAMlCiWW7E5Qj25t42mPr0aHwtuRl1uJqoNFnQbE38rD1E2noVQJPNmuCT4Y1ooFHRFRDVzJyMe0zbG4mJ4PAHi5hz9mR7SAmQm7rxBVF4u6OsrIK8GLG0+isFSJLk0d8dHTrXnJlYiomoQQ+P5YEj7YlQB5mQpONmb4+Jk26N3CRerQiPQOi7o6KFOqMHXTaaTllqCpszVWjgmDKQdFEBFVS2Z+Cd74+Rz2X7oNAOjR3BmfPNMaLrYWEkdGpJ9Y1NXBV/9cxckbd2FjboK149vD3orruBIRVcee8+mYt+Mc7hYpYGZihHkDAjG+sy+vdBDVAYu6WjqReAdf/XMFAPDBkyHwdbKWOCIiIt2XW6TAgt8vYGdsKgCgpbsdPh/ZFi3cbCWOjEj/sairhdIyFf7381moBPBUqCeGtm0idUhERDpv/8VMzN1xDhl5chjJgFd6NsWMxwM49xyRhrCoq4Xvjt5AUnYRnG3NsXBosNThEBHptJyiUiz6Ix47TpefnfNzssYnz7RBmE9jiSMjMiws6mrobmEpvvy7/LLr7IjmsDHnW0hEVBUhBP44l4aFv19AVkEpZDLgxa5+eD2iBSzNeHaOSNNYkdTQl/9cQV5JGQLdbPF0mJfU4RAR6aSUO0V497cL+OdiJgCgmYsNPnq6NUK9eXaOSFtY1NVAcakSPx5PBgDMH9gSxhylRURUgbxMiTWHEvHVP1dQolDB1FiGKb2b4dVeTdl3jkjLWNTVwJmUHJSWqeBqZ45uzZykDoeISKf8czEDi36Px43sIgBAJ38HvD+sFZq52EgcGVHDwKKuBk7duAMACPd14DJgRET/upSejw92J+Dg5fJJhF1szTHviUAMa9uEuZKoHrGoq4ET/xZ1HXwdJI6EiEh6qTnF+GLfZfwccxMqAZgay/BCVz9MeyyAg8iIJMBvXTWVKVU4nXQXABDuy46+RNRwpeeWYFXUNWw6kYzSMhUAYECIG+YOCISPIydiJ5IKi7pqupiej8JSJWzNTRDoZid1OERE9S4xqxBrDl3Htpib6mKuk78D5vQP5KhWIh2g96vPr1ixAn5+frCwsEBYWBgOHTqkldc5+e+l11Cfxhz1SkQ1VtNcFRUVhbCwMFhYWMDf3x+rVq2qp0grUqkEoi7fxsTvTqHPpwfw4/Hys3MdfB3w40sdsXliJxZ0RDpCr8/Ubd26FTNmzMCKFSvQtWtXrF69GgMGDEB8fDy8vb01+lqnbpRfeu3gx/50RFQzNc1ViYmJeOKJJzBx4kT88MMPOHLkCCZPngxnZ2c89dRTWo9XCIFLGfn47cwt/Hb2Fm7eLVbf91igCyb28EdHPw4YI9I1MiGEkDqI2urYsSNCQ0OxcuVK9baWLVti2LBhWLJkySMfn5eXB3t7e+Tm5sLO7sGXVIUQ6Lj4b2Tmy7H15U7o6O+okfiJqOaq+73VJTXNVW+88QZ+++03JCQkqLdNmjQJZ8+exdGjR6v1mjV5n27eLcK124W4frsA8bfycPhqFtJyS9T321mY4KkwT4zp6I1mLrbVen0iqp265Di9PVNXWlqKmJgYzJ07t8L2iIgIREdHV/kYuVwOuVyuvp2Xl1et10q+U4TMfDlMjWVo49Wo1jETUcNTm1x19OhRREREVNjWr18/rF27FgqFAqamppUeU9v8BgDPrTmunlvuHjMTI/Rs7owhbTzweEtXLutFpAf0tqjLysqCUqmEq6trhe2urq5IT0+v8jFLlizBwoULa/xad4sUCGliBytTE1iYMrERUfXVJlelp6dX2b6srAxZWVlwd3ev9Jja5jcAaOFmC1NjI/g7W6Opsw06+juig68DCzkiPaO3Rd099/fpEEI8sJ/HvHnzMGvWLPXtvLw8eHk9ev3Wtl6N8Me07lCq9PZKNRFJrCa56kHtq9p+T23zGwCsei6M/eOIDIDeFnVOTk4wNjau9Es3MzOz0i/ce8zNzWFubl7r1+SoVyKqqdrkKjc3tyrbm5iYwNGx6j69dclvLOiIDIPeTmliZmaGsLAwREZGVtgeGRmJLl26SBQVEVFFtclVnTt3rtR+7969CA8Pr7I/HRERoMdFHQDMmjULa9aswbp165CQkICZM2ciOTkZkyZNkjo0IiK1R+WqefPmYdy4cer2kyZNQlJSEmbNmoWEhASsW7cOa9euxezZs6XaBSLSA3p7+RUARo4ciezsbCxatAhpaWkICQnB7t274ePjI3VoRERqj8pVaWlpSE5OVrf38/PD7t27MXPmTCxfvhweHh748ssv62WOOiLSX3o9T11d6eN8V0QNHb+31cP3iUg/Nch56jThXj1bk/mciEha976vDfj3aLUwvxHpp7rkuAZd1OXn5wNAtYf9E5HuyM/Ph729vdRh6CzmNyL9Vpsc16Avv6pUKty6dQu2trb1OqT/3vxRKSkpen1ZxBD2g/ugO6q7H0II5Ofnw8PDA0ZGej3WS6ukym+AYXwmuQ+6wxD2oyb7UJcc16DP1BkZGcHT01Oy17ezs9PbD+h/GcJ+cB90R3X2g2foHk3q/AYYxmeS+6A7DGE/qrsPtc1x/JlLREREZABY1BEREREZABZ1EjA3N8e7775bpyXLdIEh7Af3QXcYyn6QYRxL7oPuMIT9qK99aNADJYiIiIgMBc/UERERERkAFnVEREREBoBFHREREZEBYFFHREREZABY1GnBBx98gC5dusDKygqNGjWq1mOEEFiwYAE8PDxgaWmJXr164cKFCxXayOVyTJs2DU5OTrC2tsaQIUNw8+ZNLexBubt372Ls2LGwt7eHvb09xo4di5ycnIc+RiaTVfn38ccfq9v06tWr0v2jRo3SmX2YMGFCpfg6depUoU19Houa7oNCocAbb7yBVq1awdraGh4eHhg3bhxu3bpVoZ22j8OKFSvg5+cHCwsLhIWF4dChQw9tHxUVhbCwMFhYWMDf3x+rVq2q1Gb79u0ICgqCubk5goKCsHPnTo3FS9VnCDmO+U038ltt9kMXcpzO5jdBGvfOO++IZcuWiVmzZgl7e/tqPWbp0qXC1tZWbN++XcTFxYmRI0cKd3d3kZeXp24zadIk0aRJExEZGSlOnz4tevfuLdq0aSPKysq0sh/9+/cXISEhIjo6WkRHR4uQkBAxaNCghz4mLS2twt+6deuETCYT165dU7fp2bOnmDhxYoV2OTk5OrMP48ePF/37968QX3Z2doU29XksaroPOTk54vHHHxdbt24VFy9eFEePHhUdO3YUYWFhFdpp8zhs2bJFmJqaim+//VbEx8eL1157TVhbW4ukpKQq21+/fl1YWVmJ1157TcTHx4tvv/1WmJqaip9//lndJjo6WhgbG4vFixeLhIQEsXjxYmFiYiKOHTumkZip+gwhxzG/6UZ+q81+SJ3jdDm/sajTovXr11cr4alUKuHm5iaWLl2q3lZSUiLs7e3FqlWrhBDlH2JTU1OxZcsWdZvU1FRhZGQk9uzZo/HY4+PjBYAKH6ijR48KAOLixYvVfp6hQ4eKPn36VNjWs2dP8dprr2kq1Aeq7T6MHz9eDB069IH31+ex0NRxOHHihABQIelo8zh06NBBTJo0qcK2wMBAMXfu3Crbz5kzRwQGBlbY9sorr4hOnTqpb48YMUL079+/Qpt+/fqJUaNGaShqqil9zXHMb0MfeL++/l9TnzlOl/MbL7/qgMTERKSnpyMiIkK9zdzcHD179kR0dDQAICYmBgqFokIbDw8PhISEqNto0tGjR2Fvb4+OHTuqt3Xq1An29vbVfr2MjAzs2rULL774YqX7fvzxRzg5OSE4OBizZ89Gfn6+xmK/py77cODAAbi4uKB58+aYOHEiMjMz1ffV57HQxHEAgNzcXMhkskqXyrRxHEpLSxETE1Ph/QGAiIiIB8Z89OjRSu379euHU6dOQaFQPLSNNj7/pFm6luOY33Qjv9V1P/6rvnKcruc3kxq1Jq1IT08HALi6ulbY7urqiqSkJHUbMzMzNG7cuFKbe4/XdEwuLi6Vtru4uFT79TZu3AhbW1sMHz68wvYxY8bAz88Pbm5uOH/+PObNm4ezZ88iMjJSI7HfU9t9GDBgAJ555hn4+PggMTERb7/9Nvr06YOYmBiYm5vX67HQxHEoKSnB3Llz8eyzz1ZYSFpbxyErKwtKpbLKz/ODYk5PT6+yfVlZGbKysuDu7v7ANtr4/JNm6VqOY37TjfxWl/34r/rMcbqe33imrpoWLFjwwE6y9/5OnTpVp9eQyWQVbgshKm27X3Xa/FdN9qOq563J661btw5jxoyBhYVFhe0TJ07E448/jpCQEIwaNQo///wz9u3bh9OnT+vEPowcORIDBw5ESEgIBg8ejD///BOXL1/Grl27HhpXTd6b+joOCoUCo0aNgkqlwooVKyrcV9fj8Cg1/TxX1f7+7bX5jlD1GEKOY37TjfxWH/txj1Q5TlfzG8/UVdPUqVMfOWrG19e3Vs/t5uYGoLyad3d3V2/PzMxUV+5ubm4oLS3F3bt3K/yCyszMRJcuXar9WtXdj3PnziEjI6PSfbdv3670a6Iqhw4dwqVLl7B169ZHtg0NDYWpqSmuXLmC0NDQR7avr324x93dHT4+Prhy5QoAzRyL+tgHhUKBESNGIDExEf/880+FX7BVqelxeBAnJycYGxtX+oX538/z/dzc3Kpsb2JiAkdHx4e2qcmxpAczhBzH/KYb+Q0w3Byn8/mtRj3wqEZq2on4ww8/VG+Ty+VVdiLeunWrus2tW7e03nn1+PHj6m3Hjh2rdufV8ePHVxqJ9CBxcXECgIiKiqp1vFWp6z7ck5WVJczNzcXGjRuFEPV7LGq7D6WlpWLYsGEiODhYZGZmVuu1NHkcOnToIF599dUK21q2bPnQjsQtW7assG3SpEmVOhIPGDCgQpv+/ftzoISE9DXHMb/9PynzmxD6meN0Ob+xqNOCpKQkERsbKxYuXChsbGxEbGysiI2NFfn5+eo2LVq0EDt27FDfXrp0qbC3txc7duwQcXFxYvTo0VUO9/f09BT79u0Tp0+fFn369NH6MPPWrVuLo0ePiqNHj4pWrVpVGmZ+/34IIURubq6wsrISK1eurPScV69eFQsXLhQnT54UiYmJYteuXSIwMFC0a9dOa9MW1GQf8vPzxeuvvy6io6NFYmKi2L9/v+jcubNo0qSJZMeipvugUCjEkCFDhKenpzhz5kyF4fxyuVwIof3jcG/I/9q1a0V8fLyYMWOGsLa2Fjdu3BBCCDF37lwxduxYdft7Q/5nzpwp4uPjxdq1aysN+T9y5IgwNjYWS5cuFQkJCWLp0qWc0kQihpDjmN90I7/VZj+kznG6nN9Y1GnB+PHjBYBKf/v371e3ASDWr1+vvq1SqcS7774r3NzchLm5uejRo4eIi4ur8LzFxcVi6tSpwsHBQVhaWopBgwaJ5ORkre1Hdna2GDNmjLC1tRW2trZizJgx4u7duxXa3L8fQgixevVqYWlpWeV8QMnJyaJHjx7CwcFBmJmZiaZNm4rp06dXmidJqn0oKioSERERwtnZWZiamgpvb28xfvz4Su9zfR6Lmu5DYmJilZ+//34G6+M4LF++XPj4+AgzMzMRGhpa4dfx+PHjRc+ePSu0P3DggGjXrp0wMzMTvr6+Vf6nuW3bNtGiRQthamoqAgMDxfbt2zUWL1WfIeQ45jfdyG+12Q9dyHG6mt9kQvzbW4+IiIiI9BZHvxIREREZABZ1RERERAaARR0RERGRAWBRR0RERGQAWNQRERERGQAWdUREREQGgEUdERERkQFgUUdERERkAFjUkUHp1asXZsyYIXUYREQax/xGj8KijoiIiMgAsKgjIiIiMgAs6sjgqFQqzJkzBw4ODnBzc8OCBQukDomISCOY3+hhWNSRwdm4cSOsra1x/PhxfPTRR1i0aBEiIyOlDouIqM6Y3+hhZEIIIXUQRJrSq1cvKJVKHDp0SL2tQ4cO6NOnD5YuXSphZEREdcP8Ro/CM3VkcFq3bl3htru7OzIzMyWKhohIc5jf6GFY1JHBMTU1rXBbJpNBpVJJFA0RkeYwv9HDsKgjIiIiMgAs6oiIiIgMAIs6IiIiIgPA0a9EREREBoBn6oiIiIgMAIs6IiIiIgPAoo6IiIjIALCoIyIiIjIALOqIiIiIDACLOiIiIiIDwKKOiIiIyACwqCMiIiIyACzqiIiIiAwAizoiIiIiA8CijoiIiMgA/B/gpafX9FjmzwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "r2_plot = sol[-1][:, 0]\n",
    "mr_plot = sol[-1][:, 1]\n",
    "\n",
    "# Plots\n",
    "ax = plt.subplot(2, 2, 1)\n",
    "ax.plot(-time_grid, r2_plot)\n",
    "ax.set_xlabel(\"h\")\n",
    "ax.set_ylabel(\"$x_0 = r^2$\")\n",
    "\n",
    "ax = plt.subplot(2, 2, 2)\n",
    "ax.plot(-time_grid, mr_plot)\n",
    "ax.set_xlabel(\"h\")\n",
    "ax.set_ylabel(\"$x_1 = M / r$\")\n",
    "\n",
    "ax = plt.subplot(2, 2, 3)\n",
    "ax.plot(-time_grid, np.sqrt(r2_plot))\n",
    "ax.set_xlabel(\"h\")\n",
    "ax.set_ylabel(\"r\")\n",
    "\n",
    "ax = plt.subplot(2, 2, 4)\n",
    "ax.plot(-time_grid, mr_plot * np.sqrt(r2_plot))\n",
    "ax.set_xlabel(\"h\")\n",
    "ax.set_ylabel(\"M\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 - Variational equations\n",
    "We seek an analytical expression for the stellar mass and radius, as a function of the variation of the conditions at the stellar center $\\delta h_c$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to obtain the expressions seeked, a variational equation for the initial value of the entalpy (the independent variable) $h_0$ is also needed. We thus introduce a last variable change $\\overline h_0 s = h$ so that the differential equations become:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\frac{dx_0}{ds} &= - \\overline h_0 \\frac{2x_0 (1-2x_1)}{4\\pi x_0 P(s) + x_1}\\\\\n",
    "\\frac{dx_1}{ds} &= \\left(2\\pi \\epsilon(s) - \\frac{x_1}{2x_0}\\right) \\frac{dx_0}{ds} \\\\\n",
    "\\end{cases},\n",
    "$$\n",
    "\n",
    "with:\n",
    "\n",
    "$$\n",
    "P(s) = K^{-\\frac{1}{\\Gamma-1}}\\left[\\frac{\\Gamma-1}{\\Gamma}\\left(e^{\\overline h_0 s}-1\\right)\\right]^{\\frac{\\Gamma}{\\Gamma-1}},\n",
    "$$\n",
    "\n",
    "and:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\rho(s) &= \\left(\\frac{P(s)}{K}\\right)^{1/\\Gamma} \\\\\n",
    "\\epsilon(s) &= \\rho + \\frac{P(s)}{\\Gamma-1}\n",
    "\\end{cases},\n",
    "$$\n",
    "\n",
    "and must be integrated for $s \\in [1, 0]$, which will correspond to $h \\in [\\overline h_0, 0]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The independent variable is s, not h = log(entalpy) nor t.\n",
    "s = hy.time\n",
    "\n",
    "# Important expressions\n",
    "P = K ** (-1.0 / (Gamma - 1.0)) * (\n",
    "    (Gamma - 1.0) / Gamma * (hy.exp(hy.par[2] * s) - 1.0)\n",
    ") ** (Gamma / (Gamma - 1.0))\n",
    "rho = (P / K) ** (1.0 / Gamma)\n",
    "eps = rho + P / (Gamma - 1.0)\n",
    "\n",
    "# We define the dynamics\n",
    "dx0ds = (-2.0 * x0 * (1.0 - 2.0 * x1) / (4.0 * np.pi * x0 * P + x1)) * hy.par[2]\n",
    "dyn2 = [\n",
    "    (x0, dx0ds),\n",
    "    (x1, (2.0 * np.pi * eps - x1 / 2.0 / x0) * dx0ds),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We augment the dynamics with the variational equations (we use 8th order but lower orders are good already)\n",
    "var_sys = hy.var_ode_sys(dyn2, args=[x0, x1, hy.par[2]], order=8)\n",
    "\n",
    "# We instantiate the Taylor adaptive integrator for the system of equations augmented with the variational ones\n",
    "ta_var = hy.taylor_adaptive(\n",
    "    var_sys,\n",
    "    state=[1.0, 1.0],\n",
    "    time=0.1,\n",
    "    tol=1e-18,\n",
    "    compact_mode=True,\n",
    ")\n",
    "# We copy, for future reference, the initial conditions on the variational state\n",
    "ic_var = list(deepcopy(ta_var.state[2:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total time to propagate: 13.353347539901733\n",
      "Outcome:  taylor_outcome.time_limit\n",
      "Final value for the stellar radius is:  5.51113483059698\n",
      "Final value for the stellar mass is:  1.3684613613885492\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "ta_var.time = 1.0\n",
    "\n",
    "ta_var.state[:] = y_0 + ic_var\n",
    "ta_var.pars[:] = [0.0, 0.0, h_0]\n",
    "\n",
    "s_grid = np.linspace(1.0, 0.0, 100)\n",
    "sol_var = ta_var.propagate_grid(s_grid)\n",
    "print(\"Total time to propagate:\", time.time() - start)\n",
    "print(\"Outcome: \", sol_var[0])\n",
    "\n",
    "rf = np.sqrt(sol_var[-1][:, 0][-1])\n",
    "Mf = sol_var[-1][:, 1][-1] * rf\n",
    "\n",
    "print(\"Final value for the stellar radius is: \", rf)\n",
    "print(\"Final value for the stellar mass is: \", Mf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have thus computed the following polynomial approximations (degree $k$):\n",
    "\n",
    "$$\n",
    "\\begin{array}{l}\n",
    "R_{\\odot} = \\mathcal P^k_R(\\delta x_{0_0}, \\delta x_{1_0},\\delta h_0) \\\\\n",
    "M_{\\odot} = \\mathcal P^k_M(\\delta x_{0_0}, \\delta x_{1_0}, \\delta h_0) \n",
    "\\end{array}.\n",
    "$$\n",
    "\n",
    "Note the underscript $0$ refers to the initial values, not to the star center which would here be indicated using $c$ rather than $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Conditions at the star center:\n",
      "\n",
      "Pressure:  0.04\n",
      "Energy density:  0.06\n",
      "Density:  0.02\n",
      "log-entalpy:  1.6094379124341003\n"
     ]
    }
   ],
   "source": [
    "# We repeat the computation changing the pressure at the center.\n",
    "\n",
    "# Pressure\n",
    "P_c_new = P_c * (1.0 + 3.0)\n",
    "\n",
    "# Density\n",
    "rho_c_new = (P_c_new / K_v) ** (1.0 / Gamma_v)\n",
    "\n",
    "# Energy density\n",
    "eps_c_new = rho_c_new + P_c_new / (Gamma_v - 1.0)\n",
    "\n",
    "# log-entalpy\n",
    "h_c_new = np.log(\n",
    "    1\n",
    "    + Gamma_v / (Gamma_v - 1) * K_v ** (1 / Gamma_v) * P_c_new ** (1.0 - 1.0 / Gamma_v)\n",
    ")\n",
    "\n",
    "print(\"Conditions at the star center:\")\n",
    "print(\"\\nPressure: \", P_c_new)\n",
    "print(\"Energy density: \", eps_c_new)\n",
    "print(\"Density: \", rho_c_new)\n",
    "print(\"log-entalpy: \", h_c_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We may now find the new starting conditions at $h_c-dh$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_0_new = h_c_new - dh\n",
    "\n",
    "# Initial conditions at h_0_new (from Steil thesis). Can be improved at second order.\n",
    "x0_0_new = 3.0 / (2.0 * np.pi * (3.0 * P_c_new + eps_c_new)) * dh\n",
    "x1_0_new = 2.0 * eps_c_new / (3 * P_c_new + eps_c_new) * dh\n",
    "y_0_new = [x0_0_new, x1_0_new]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the Taylor approximation to compute the stellar radius and mass:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final value for the stellar radius is:  4.902302695454508\n",
      "Final value for the stellar mass is:  1.1364855355000758\n"
     ]
    }
   ],
   "source": [
    "taylor_approx = ta_var.eval_taylor_map(\n",
    "    [x0_0_new - x0_0, x1_0_new - x1_0, h_0_new - h_0]\n",
    ")\n",
    "rf_new_taylor = np.sqrt(taylor_approx[0])\n",
    "Mf_new_taylor = taylor_approx[1] * rf_new_taylor\n",
    "\n",
    "print(\"Final value for the stellar radius is: \", rf_new_taylor)\n",
    "print(\"Final value for the stellar mass is: \", Mf_new_taylor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute these values by numerical integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total time to propagate: 0.0002846717834472656\n",
      "Outcome (should be time_limit if all went well):  time_limit\n",
      "Final value for the stellar radius is:  4.9020827509255644\n",
      "Final value for the stellar mass is:  1.13645940407093\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "ta.time = h_0_new\n",
    "ta.state[:] = y_0_new\n",
    "ta.pars[:] = [0.0, 0.0]\n",
    "time_grid = np.linspace(h_0_new, 0.0, 100)\n",
    "sol_new = ta.propagate_grid(time_grid)\n",
    "time_cost = time.time() - start\n",
    "print(\"Total time to propagate:\", time_cost)\n",
    "print(\n",
    "    \"Outcome (should be time_limit if all went well): \", str(sol_new[0]).split(\".\")[1]\n",
    ")\n",
    "\n",
    "# The mass and radius are easily computed from the last value reached\n",
    "rf_new = np.sqrt(sol_new[-1][:, 0][-1])\n",
    "Mf_new = sol_new[-1][:, 1][-1] * rf_new\n",
    "\n",
    "print(\"Final value for the stellar radius is: \", rf_new)\n",
    "print(\"Final value for the stellar mass is: \", Mf_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets see the error introduced by the Taylor approximation in this specific instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Absolute error on mass: -2.61314291458703e-05\n",
      "Absolute error on radius: -0.00021994452894347205\n",
      "\n",
      "Relative error on mass: -2.299371983923445e-05\n",
      "Relative error on radius: -4.486756754604851e-05\n"
     ]
    }
   ],
   "source": [
    "print(\"Absolute error on mass:\", Mf_new - Mf_new_taylor)\n",
    "print(\"Absolute error on radius:\", rf_new - rf_new_taylor)\n",
    "\n",
    "print(\"\\nRelative error on mass:\", (Mf_new - Mf_new_taylor) / Mf_new)\n",
    "print(\"Relative error on radius:\", (rf_new - rf_new_taylor) / rf_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The $M_\\odot$ vs $R_\\odot$ plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now loop over different star conditions (at the core) and plot the relation between mass and radius using the Taylor approximation and the numerical integrator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total computation time using the Taylor Map: 0.0064754486083984375\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "Rstar = []\n",
    "Mstar = []\n",
    "P_grid = np.linspace(P_c * (1.0 - 0.95), P_c * (1 + 20.0), 1000)\n",
    "for P_c_new in P_grid:\n",
    "    # Density\n",
    "    rho_c_new = (P_c_new / K_v) ** (1.0 / Gamma_v)\n",
    "\n",
    "    # Energy density\n",
    "    eps_c_new = rho_c_new + P_c_new / (Gamma_v - 1.0)\n",
    "\n",
    "    # log-entalpy\n",
    "    h_c_new = np.log(\n",
    "        1\n",
    "        + Gamma_v\n",
    "        / (Gamma_v - 1)\n",
    "        * K_v ** (1 / Gamma_v)\n",
    "        * P_c_new ** (1.0 - 1.0 / Gamma_v)\n",
    "    )\n",
    "\n",
    "    # Since the dynamics is singular at h_c, we define the initial conditions at h_0 = h_c - dh\n",
    "    h_0_new = h_c_new - dh\n",
    "\n",
    "    # Initial conditions at h_0 (from Steil thesis). Can be improved at second order.\n",
    "    x0_0_new = 3.0 / (2.0 * np.pi * (3.0 * P_c_new + eps_c_new)) * dh\n",
    "    x1_0_new = 2.0 * eps_c_new / (3 * P_c_new + eps_c_new) * dh\n",
    "    taylor_approx = ta_var.eval_taylor_map(\n",
    "        [x0_0_new - x0_0, x1_0_new - x1_0, h_0_new - h_0]\n",
    "    )\n",
    "    Rstar.append(np.sqrt(taylor_approx[0]))\n",
    "    Mstar.append(taylor_approx[1] * np.sqrt(taylor_approx[0]))\n",
    "time_cost = time.time() - start\n",
    "print(\"Total computation time using the Taylor Map:\", time_cost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total time using the numerical integrator: 0.10077452659606934\n"
     ]
    }
   ],
   "source": [
    "start = time.time()\n",
    "Rstar_gt = []\n",
    "Mstar_gt = []\n",
    "P_grid = P_c * np.logspace(0.01, 12, 1000) / 40.0\n",
    "for P_c_new in P_grid:\n",
    "    # Density\n",
    "    rho_c_new = (P_c_new / K_v) ** (1.0 / Gamma_v)\n",
    "\n",
    "    # Energy density\n",
    "    eps_c_new = rho_c_new + P_c_new / (Gamma_v - 1.0)\n",
    "\n",
    "    # log-entalpy\n",
    "    h_c_new = np.log(\n",
    "        1\n",
    "        + Gamma_v\n",
    "        / (Gamma_v - 1)\n",
    "        * K_v ** (1 / Gamma_v)\n",
    "        * P_c_new ** (1.0 - 1.0 / Gamma_v)\n",
    "    )\n",
    "\n",
    "    # Since the dynamics is singular at h_c, we define the initial conditions at h_0 = h_c - dh\n",
    "    h_0_new = h_c_new - dh\n",
    "\n",
    "    # Since the dynamics is singular at h_c, we define the initial conditions at h_0 = h_c - dh\n",
    "    h_0_new = h_c_new - dh\n",
    "\n",
    "    # Initial conditions at h_0 (from Steil thesis). Can be improved at second order.\n",
    "    x0_0_new = 3.0 / (2.0 * np.pi * (3.0 * P_c_new + eps_c_new)) * dh\n",
    "    x1_0_new = 2.0 * eps_c_new / (3 * P_c_new + eps_c_new) * dh\n",
    "\n",
    "    ta.time = h_0_new\n",
    "    ta.state[:] = y_0_new\n",
    "    outcome = ta.propagate_until(0.0)\n",
    "\n",
    "    # The mass and radius are easily computed from the last value reached\n",
    "    Rstar_gt.append(np.sqrt(ta.state[0]))\n",
    "    Mstar_gt.append(ta.state[1] * np.sqrt(ta.state[0]))\n",
    "time_cost = time.time() - start\n",
    "print(\"Total time using the numerical integrator:\", time_cost)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And lets have a look at the difference in the computed values for this range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa4d7f1b0b0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAGyCAYAAAAGdNXrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABd30lEQVR4nO3dd3gU5cLG4d/upleSkEACCYTeIYAgggoqchCx8YmISrNh53BAxQY2EAuix4J6BEQUG0oRBZRepBN6C52Q0NP77nx/BFcjxQSTzGbz3Ne1F5nJ7O6zScg+mXnnHYthGAYiIiIibsZqdgARERGRsqCSIyIiIm5JJUdERETckkqOiIiIuCWVHBEREXFLKjkiIiLillRyRERExC15mB3ALA6Hg6NHjxIYGIjFYjE7joiIiBSDYRikp6cTFRWF1XrxfTWVtuQcPXqU6Ohos2OIiIjIJTh8+DA1a9a86DaVtuQEBgYChV+koKAgk9OIiIhIcaSlpREdHe18H7+YSltyfj9EFRQUpJIjIiJSwRRnqIkGHouIiIhbUskRERERt6SSIyIiIm6p0o7JKS673U5+fr7ZMUTKnKenJzabzewYIiKlRiXnAgzDIDk5mZSUFLOjiJSbKlWqUL16dc0dJSJuQSXnAn4vOBEREfj5+emXvrg1wzDIysri+PHjAERGRpqcSETkn1PJOQ+73e4sOGFhYWbHESkXvr6+ABw/fpyIiAgduhKRCk8Dj8/j9zE4fn5+JicRKV+//8xrHJqIuAOVnIvQISqpbPQzLyLuRCVHRERE3JJKjriEUaNG0apVK7Nj0LlzZ4YMGWJ2DBERKQUqOW4mOTmZJ554gnr16uHj40O1atXo1KkTEyZMICsry+x4l2zx4sVYLJZSO6W/tB9PRERcj86uciP79u2jY8eOVKlShdGjR9O8eXMKCgrYvXs3EydOJCoqiptuuum8983Pz8fT07OcE5e+vLw8vLy8zI4hIiIuQCXHjTz88MN4eHiwbt06/P39neubN29Or169MAzDuc5isfDhhx/y888/8+uvvzJs2DBefPFFPvzwQ958800OHz5MbGwszz33HPfccw8ABw4cIDY2lo0bNzoPLaWkpBASEsKiRYvo3LkzixcvpkuXLvz666889dRTbN++nVatWjFp0iQaNmzofP7XXnuNt99+m6ysLHr37k14ePgFX9eBAwfo0qULACEhIQD079+fyZMn07lzZ5o1a4aXlxdTpkyhadOmfPbZZxfNWbt27Qs+HoDD4eDJJ5/kf//7H15eXgwePJhRo0Zd+jdG5B8wDIP0nHxOJZ8iJfEEGcdOkpOeQVZuAVm5BWTnFpCVm092XgH2AjsFdgd2uwPH2X8L7A4wDGxWsAE2C3hYzn5sBR8PKz4eVnx//9fTii0wAO/wqgRUDycoKpyQyKqE+nvjYdPOf6lYVHKKyTAMsvPt5f68vp62Yp3xcurUKebPn8/o0aOLFJw/++vjjBw5kjFjxvD2229js9n44YcfeOKJJxg/fjzXXXcdP/74IwMHDqRmzZrOUlBczz77LG+99Rbh4eEMHjyYQYMGsWLFCgC++eYbRo4cyfvvv8+VV17J559/zrvvvkudOnXO+1jR0dFMnz6dXr16sWvXLoKCgpxzugB89tlnPPTQQ6xYsaJIkbuQ4jze0KFDWb16Nb/99hsDBgygY8eOdO3atURfA5HiyMm3c+R0FolHT3LyUDKnE0+QmnSC3JOnsJ85g+VMCr7ZGXg4Lv77x/PsraTsQObZ2wW3sVhJ8w2gIKgKlpBQfKuFERgZTpWoaoRHVyOqWhVqhPgS5u+lM/TEpajkFFN2vp0mL8wr9+fd/lI3/Lz+/tuUkJCAYRhF9pYAVK1alZycHAAeeeQRxo4d6/xc3759GTRoUJHlAQMG8PDDDwMwdOhQVq1axZtvvlnikvPqq69y9dVXA/D000/To0cPcnJy8PHxYfz48QwaNIj77rsPgFdeeYVff/3VmfOvbDYboaGhAERERFClSpUin69Xrx6vv/66c/nAgQMXzfZ3j9eiRQtGjhwJQP369XnvvfdYsGCBSo5cOsPg9MlUDiQcIfFAMscPJXEm8QTpySdxpJwhKCcTT0eBc3MPzv3l7OlhwwgOhuBgrIEB+Hh54ONpK/zXywNfLw88PKx42KzYbFZsHjZsVmvhpI4WCw7AAdiNwmLjMKDAAdl2Bzn5DrILHGTlO8jNK8BIT8c4fQbOnMGalkpObj4hWWmQlQbJh2BH4WOdPnvb6OlDqk8AeYHBBEaGE1KjGtVqVadGnShia0dSI8wfm1XlR8qfSo6b+etfUWvWrMHhcHDXXXeRm5tb5HNt27Ytsrxjxw4eeOCBIus6duzIO++8U+IcLVq0cH78+yUCjh8/TkxMDDt27GDw4MFFtu/QoQOLFi0q8fPAua/jn/pzdijM//vlDkT+jmEYHE5KYe+mXSRt2U3q7n0UHDyMPTu7yHYeQMjZj71sVgIDfPAODca7aigB1aoSGBlOQLVwqkSFExIVgW/VEDBjFmqHA3tKKqlHj5F29ASpScdJO3qcrGMnyDl+ipzUDNJz8vFLPwnpJ+HoXlgPxym8rbNYyfEPJKB6BFVjqhFVO4pa9WoQU7cGHlXDwNu7/F+TVBoqOcXk62lj+0vdTHne4qhXrx4Wi4WdO3cWWf/7IaA/H4753fkOa/21JBmG4VxntVqd6353oZlx/zyI+ff7OxyOv30dl+Kvr6MkOc/nrwOwLRZLmWWXii81K4/Nm/ezd+MOTm5PIGf/IfxTTmHh3EOnnsFB+FWrSnBkOGE1IqhWqzqRMdUJiorAEhwMHi74K9lqxRYaQmhoCKHNGp37+ZwcSEkh9/hJkg8mcexgEqePHCM96QTZJ06RmpGDLSMVIyGVEwl7OAFsAmxWC+GB3lSNCKF6rUhi6tckqG5tiI6G6tXNKXTidlzwf5RrslgsxTpsZJawsDC6du3Ke++9x2OPPXbBcTkX07hxY5YvX06/fv2c61auXEnjxo0BnIODk5KSiIuLAyA+Pv6SnmfVqlVFnmfVqlUXvc/vZ0zZ7X8/Lqo4OUvyeCJ/dvREGlvW7ODgpp2c2rkPjhzGN6/wUKsHEADYLBZ8wqviW7c24Y3rEd2iPrUb18bP/9w/Nio8Hx+oXh3v6tWp1aIZtf78OYcDe2oahxMOs2/PERL3JXH8UBJpR4/jk5mGPTWH5NQktu5Jgl83EOLnRVQVX6pXDSSqSV3CGteDmBioWRPO84eayN9x3XdtKbEPPviAjh070rZtW0aNGkWLFi2wWq2sXbuWnTt30qZNm4vef/jw4fTu3ZvWrVtz7bXXMnv2bL7//nt+/fVXoHBv0OWXX85rr71G7dq1OXnyJM8991yJcz7xxBP079+ftm3b0qlTJ7744gu2bdt2wYHHALVq1cJisfDjjz9yww034OvrS0BAwHm3LU7OkjyeVG6pp1LZsjyevet3cHJHApbkZKxG4Z693992A/19CIiNIbxJPaKb16dWq0Z4hQSbF9pVWK3YQqpQ+7Iq1L6suXO1w2Gw/1QmW/cdY/u2Q+zdfZjUw0lEpp8k6fhJfI6mwuYjBPuupFaYHzGh/tRoWAvfOrULS090NISEgAY5y99QyXEjdevWZePGjYwePZoRI0Zw5MgRvL29adKkCcOGDXMOKL6QW265hXfeeYc33niDxx9/nNjYWCZNmkTnzp2d20ycOJFBgwbRtm1bGjZsyOuvv871119fopx33HEHe/fu5amnniInJ4devXrx0EMPMW/ehQd216hRgxdffJGnn36agQMH0q9fP+cp3+fzdzlL+nhSeTgcBtu2HWTrglWcWL8F+4GDWM6WGhtgAYLCQwhqUIcazetTv00TwurXds1DTS7KarVQNzyAuuEB3Ny+LgApWXmsP3iGtftPs2nrPk7v3Eu1lOMcOHWSkCNHsWxJonpQPLXC/KhT1Z/wyDAsvxeemBiIjNQhLjmHxSjOObduKC0tjeDgYFJTUwkKCiryuZycHPbv309sbCw+Pj4mJRQpf5X1Zz87t4C1v21j1+I1nN6wBe/TJ4tuEBFBWNMG1GndmGaXNyO4elXtRShj6Tn5rNp3muV7TrB222Fy9x8kMu0kUWknqJZxiireNupU9aduRAA1qvhi8/KEGjX+KD01a4Kfn9kvQ8rAxd6//0p/eohIpZSakcPKX9eye8k6srZuxze7cKYYb8DDw4OghnWIah9Hk2suI6p2lLlhK6FAH0+6NqlG1ybV4OZmHDmTxfI9J1m06zg/7TxG0JnjRKadpMaBE8RmJ9I42Ea9oynU3rf/j0kLq1b94/BWTAyEhqqcVjIqOSJSaaSlpLN63ip2L11P9vadeOTnAYVja3wDfAlt0YS6nVrTrEs7vAJLPnhfyk7NED/6tIuhT7sYcvLtLN9zkl+2H+PXHceYlZFLSHYaUeknqbvnNJ39cmntX0CMw8B68iRs2FD4IP7+hYXnz4e4dJjRrem7KyJuLScnj9/mrWLP3KVkb9uFcXbmYA/AP7QK1S5rSdNr21HnsmZY3OD6bZWBj6eN65pU47om1bA7DDYeOsO8bcnM2ZzErNQcZgE++Tk0SU2lZ6idzr451M5NwZKZCTt3Ft6gcAxPVNQfe3uiowuLkLgNlRwRcTuGYbBp/W42zFrMmd/WYsvOcn7OFhFB1OWtaN31cmJbNdThiwrOZrXQtnYobWuHMqJ7Y9YdPMOsTYn8tCWZDZk+bMiGF7OhXogP/WK9uDG0gNCTx+DwYcjMLPz38OE/HjAsrGjpqarxVxWZBh5r4LGIU0X/2T+afIYl0xdxZNEKPJOTnOs9goKIuPIy2tx0DXWbxpqYUMpLvt3BioSTzNp0lLlbk8nKK9yDZ7Na6NIwgt5tatAlwgPPo4mFJefQIThx4twH8vUtLDt160L79uX8KuR8NPBYRCoNu93BbwvWsWnWYnI2b8HmsOMJ2DxsVGnZjKY3Xk2rzm2xeurXXWXiabPSuWEEnRtG8NLNBczZfJRv1h1h/cEz/LqjcCxPeKA3vdvW5O6rriMy2Beys//Ys3PoECQmFq7bvRvsdpWcCkj/60WkQko+lMyyb38lcdFKSEkBCuexCYqOos71nWh3c2f8QquYGVFcRIC3B3dcFsMdl8WQcDydb9Yd4fsNRziRnsv7i/YyYck+ujWtRr8OtWlfvz6WBg0K72i3Q3JyYeH5mz0G4pp0uEqHq0ScXP1n3zAM4ldtY+20OWRu3PLHJH0+PlS7og3tel1H7RYNNIZC/la+3cGv24/x2W8HWLXvtHN9o+qB9L+iNje3inLpS/lUZiU5XGUtp0wXtXTpUnr27ElUVBQWi4UZM2b87X1yc3N59tlnqVWrFt7e3tStW5eJEyeWfVgRKXe5efnM/24R7941giXDR5O1YRMWw4Fn3Tq0fHwg98/6gN4vPULtlhpILMXjabPSvXkkXz3QgblDrqRv+xh8PW3sTE5nxPdb6DBmIW/N38XJjFyzo8o/4BI1NTMzk5YtWzJw4EB69epVrPv07t2bY8eO8emnn1KvXj2OHz9OQUFBGSetvEaNGsWMGTMu6YKcFcmnn37K119/zfz580v9sQcMGEBKSkqxSnxp27JlC927d2fXrl2XdPFWs5w6nc4vX8zlyNxFeKSmAGCx2Qi9vDUd7upJvRb1zA0obqFR9SBG39qcp7o14tv1h/l81UEOnsrivwsT+HjpPm5vW5P7r6xDrbCK839HCrlEyenevTvdu3cv9vZz585lyZIl7Nu3j9DQUABq165dRukqBsvf/PXav39/l702U+3atTl48CDTpk2jT58+RT7XtGlTtm/fzqRJkxgwYECZ5sjNzeWFF17gq6++KtPnKQu7d+9m+PDhrFixgry8PJo3b84rr7xCly5dAGjevDnt2rXj7bffvqSLqpa3IweSWfDZLE4u+Q1bXi4egEeAH9Fdr+LKu3oQWj3M7IjihoL9PLnvyjoM7BjL/G3JTFiyl01HUpm66hBfrj5E9+aRDL6qLs1r6uKrFYVLlJySmjVrFm3btuX111/n888/x9/fn5tuuomXX34ZX1/f894nNzeX3Nw/djumpaWVV9xykZT0x+myX3/9NS+88AK7du1yrrvQ16U85efn43mBydaio6OZNGlSkZKzatUqkpOTy23Pw/Tp0wkICODKK6/8R49zsdf5TxiGgd1ux+M8M7T26NGDBg0asHDhQnx9fRk/fjw33ngje/fupXr16gAMHDiQwYMHM2LECGwueiHDPZv3sOyzWaSu3YjF4cAG+FaPoPGt19Phtmvw9HW9cULifmxWC92bR/KvZtVZvf80E5bsZfGuE8zZnMSczUlc3SCcIdfVJy4mxOyo8jdcYkxOSe3bt4/ly5ezdetWfvjhB8aPH893333HI488csH7jBkzhuDgYOctOjq6HBOXverVqztvwcHBWCwW57KnpyeDBw+mZs2a+Pn50bx5c6ZNm+a875QpUwgLCytSAgF69epFv379zvt8DoeDl156iZo1a+Lt7U2rVq2YO3eu8/MHDhzAYrHwzTff0LlzZ3x8fJg6deoF8991110sWbKEw3+alGvixIncdddd57ypjxs3jubNm+Pv7090dDQPP/wwGRkZzs9PnjyZKlWqMGPGDBo0aICPjw9du3Yt8tjn89VXX3HTTTeVyuu02+0MHTqUKlWqEBYWxpNPPslfx/gbhsHrr79OnTp18PX1pWXLlnz33XfOzy9evBiLxcK8efNo27Yt3t7eLFu27JzcJ0+eJCEhgaeffpoWLVpQv359XnvtNbKysti2bZtzu27dunHq1CmWLFly0a+DGbas28mEh15lzsMjSVu9HovDgV+Dulz1wuM88M04rrrrBhUcKXcWi4XL64QxeWA7fn7iSm6Lq4HNamHJ7hPc+sFKBkxaw8ZDZ8yOKRdjuBjA+OGHHy66TdeuXQ0fHx8jJSXFuW769OmGxWIxsrKyznufnJwcIzU11Xk7fPiwARipqannbJudnW1s377dyM7O/mOlw2EYubnlf3M4Svw1nDRpkhEcHOxcPnLkiPHGG28YGzduNPbu3Wu8++67hs1mM1atWmUYhmFkZWUZwcHBxjfffOO8z4kTJwwvLy9j4cKFhmEYxsiRI42WLVs6Pz9u3DgjKCjImDZtmrFz507jySefNDw9PY3du3cbhmEY+/fvNwCjdu3axvTp0419+/YZiYmJ581bq1Yt4+233zZuuukm4+WXXzYMwzAyMzONoKAgY+PGjUZwcLAxadIk5/Zvv/22sXDhQmPfvn3GggULjIYNGxoPPfRQkdfv6elptG3b1li5cqWxbt06o127dsYVV1xx0a9blSpVjK+++qrIukt9nWPHjjWCg4ON7777zti+fbtx7733GoGBgcbNN9/sfOxnnnnGaNSokTF37lxj7969xqRJkwxvb29j8eLFhmEYxqJFiwzAaNGihTF//nwjISHBOHny5Dm5HQ6H0bhxY+O+++4zMjIyjPz8fOONN94wqlWrZpw5c6bItu3atTNGjRp1wa/BeX/2y9DmdTuN9we/Yrzd8U7j7Y53GuM69TU+eXyssX3tjnJ5fpGSOnAywxj2TbxRZ8Qco9ZTPxq1nvrR6D9xtbHh4Gmzo1UaqampF3z//qsKebgqMjKSGjVqEBz8x3HRxo0bYxgGR44coX79+ufcx9vbG29v70t/0vx8GD360u9/qZ55Bry8/tFD1KhRg2HDhjmXH3vsMebOncu3335L+/bt8fX1pW/fvkyaNInbb78dgC+++IKaNWvSuXPn8z7mm2++yVNPPeU8vDR27FgWLVrE+PHjef/9953bDRkyhNtuu61YOQcNGsR//vMfnn32Wb777jvq1q1Lq1atztluyJAhzo9jY2N5+eWXeeihh/jggw+c6/Pz83nvvfdof3byrs8++4zGjRuzZs0a2rVrd85jpqSkkJKSQlRU0atNX+rrHD9+PCNGjHAOpJ8wYQLz5s1zfj4zM5Nx48axcOFCOnToAECdOnVYvnw5H330EVdffbVz25deeomuXbte8OtmsVj45ZdfuPnmmwkMDMRqtVKtWjXmzp1LlSpVimxbo0YNDhw4cMHHKi9bN+xmySffkb9la+EKi5XQdq24+v7/o1aj2qZmE7mYWmH+vHF7Sx69ph7vLUzg+42JLN51gsW7TnBNowiGd2tI40jNqeMqKmTJ6dixI99++y0ZGRkEBAQAhQMvrVYrNWvWNDmd67Hb7bz22mt8/fXXJCYmOscn/Xmsy/33389ll11GYmIiNWrUcA70Pd+A5rS0NI4ePUrHjh2LrO/YsSObNm0qsq5t27bFztmjRw8efPBBli5dysSJExk0aNB5t1u0aBGjR49m+/btpKWlUVBQQE5ODpmZmc7X5OHhUeS5GzVqRJUqVdixY8d5S052djZAkblhLvV1pqamkpSU5Cwvf85jnD1ktX37dnJycs4pL3l5ecTFxV3wsc/HMAwefvhhIiIiWLZsGb6+vvzvf//jxhtvZO3atURGRjq39fX1JSsr6yKPVrZ2bEpg0YRvyPtLuelyfy+iG+lyC1JxnK/sLNx5nEW7jnNrXA2Gdm1AzRA/s2NWei5RcjIyMkhISHAu79+/n/j4eEJDQ4mJiWHEiBEkJiYyZcoUAPr27cvLL7/MwIEDefHFFzl58iTDhw9n0KBBZTfA1tOzcK9KeSuFAaxvvfUWb7/9NuPHj3eOZRkyZAh5eXnObeLi4mjZsiVTpkyhW7dubNmyhdmzZ1/0cf9agAzDOGddSQYNe3h4cM899zBy5EhWr17NDz/8cM42Bw8e5IYbbmDw4MG8/PLLhIaGsnz5cu69917y8/Mvmu9C6wDCwsKwWCycOXPu8fXSfp1QONYHYM6cOdSoUaPI5/66x/HvHnvhwoX8+OOPnDlzxjkx1gcffMAvv/zCZ599xtNPP+3c9vTp09StW7dEWUvDoV0H+eWDr0lbvwkLhV+/0Mta0fmB/1O5kQrt97LzcJd6vDlvF3O2JPH9hkR+3JREvw61eKRLPUL8/9neeLl0LjHweN26dcTFxTn/gh06dChxcXG88MILQOGZQ4cOHXJuHxAQwC+//EJKSgpt27blrrvuomfPnrz77rtlF9JiKTxsVN63UpjYbNmyZdx8883cfffdtGzZkjp16rBnz55ztrvvvvuYNGkSEydO5Lrrrrvg4OygoCCioqJYvnx5kfUrV66kcePG/yjroEGDWLJkCTfffDMhIeeeubBu3ToKCgp46623uPzyy2nQoAFHjx49Z7uCggLWrVvnXN61axcpKSk0atTovM/r5eVFkyZN2L59u3Pdpb7O4OBgIiMjWbVqVZE869evdy43adIEb29vDh06RL169YrcSjoo/vc9M1Zr0f/OVqvVWaZ+t3Xr1nP2FJWl40mnmPLc+3x//7Okr4/HgkFwm5bc+tHL3DNuuAqOuI3Yqv68f1drZjzSkcvrhJJnd/C/5fu56o1FfLA4gZx8u9kRKyWX2JPTuXPnc848+bPzze/SqFEjfvnllzJM5T7q1avH9OnTWblyJSEhIYwbN47k5ORz3qjvuusuhg0bxieffOLca3Yhw4cPZ+TIkc5xM5MmTSI+Pp4vvvjiH2Vt3LgxJ0+exM/v/Lt569atS0FBAf/973/p2bMnK1asYMKECeds5+npyWOPPca7776Lp6cnjz76KJdffvl5D1X9rlu3bixfvrzImJ9LfZ1PPPEEr732GvXr16dx48aMGzeOlLPXVwIIDAxk2LBh/Pvf/8bhcNCpUyfS0tJYuXIlAQEB9O/f/+JfqD/p0KEDISEh9O/fnxdeeAFfX18++eQT9u/fT48ePZzbHThwgMTERK677rpiP/alysjIZvZH35M051eseYVn7fk0bUyXh+6gYasGZf78ImZpFV2FafdfzpLdJ3jt553sTE7n9bm7+HL1IZ69oTH/alb9b+c1k9LjEiVHytbzzz/P/v376datG35+fjzwwAPccsstpKamFtkuKCiIXr16MWfOHG655ZaLPubjjz9OWloa//nPfzh+/DhNmjRh1qxZ5x30XVJhYRee6K1Vq1aMGzeOsWPHMmLECK666irGjBlzzqnufn5+PPXUU/Tt25cjR47QqVOnv73sx/3330/r1q1JTU11Dmq/1Nf5n//8h6SkJAYMGIDVamXQoEHceuutRb7mL7/8MhEREYwZM4Z9+/ZRpUoVWrduzTMlPCxatWpV5s6dy7PPPss111xDfn4+TZs2ZebMmbRs2dK53bRp07j++uupVatWiR6/JOx2B79Mm8+OL2ZgSU/DCnjH1OTywXcSd1X57UESMZPFYqFzwwiuqh/OjPhEXp+7iyNnsnnoiw1cXieUkT2banByOdEFOnWBziK6du1K48aNy/bQXxmbPHkyQ4YMKbLnpLh69+5NXFwcI0aMKP1gJsrNzaV+/fpMmzbtnIHUf/ZPfvbXLlzPbxOm4Th7+NAzNIRWA3pxxS2dsVhd4si4iCmy8gqYsHgvHy3dR26BA6sF7mwXw3+ub0ioxuuUWEku0Kk9OQIUDkidP38+Cxcu5L333jM7jmneeOMNZs2aZXaMUnfw4EGeffbZixacS7V/10F+GTeFrG07ALD4+FDvtn/RdeDNePn+g2kbRNyEn5cHQ69vSO/Lohnz807mbE7ii9WHmL3pKMO6NeSu9rWwWXUIqyyo5AgArVu35syZM4wdO5aGDRuaHcc0tWrV4rHHHjM7Rqlr0KABDRqU7liYzPQsZr73FcfnLgK7HaxWwjt3pMdjfagSrunuRf6qZogf7/dtTb/LT/Hi7O1sT0rjhZnb+G79EV69pbmuiVUGdLhKh6tEnIrzs284HCyeuYxN//sKzo4x8mvamK7/vodYTeQnUix2h8GXqw/y+rxdpOcUYLVAvw61GXp9A4J8Sv/ad+5Eh6tEpEwk7DjAvDcmk797NwC2kBDaDr6TDjd0LJXpDkQqC5vVwj0datOtWXVenbODmfFHmbzyAD9tSeKFnk3o0TxSZ2GVApUcEflb2Zk5/PDuNI79vBCLw47FZqPmjV3p8UhvfPy0t1PkUkUE+vBOnzhubxPN8zO3sv9kJo9+uZFvGxxh9G3NqVGljCa4rSRUci7irxOpibi78/3Mr5m3ilXvfY7jzBksQEDTRnQfPpAa9Uo2aaGIXFin+lX5+YkrmbBkLx8s2suS3Se4ftwSRtzQmL7tYrBqYPIl0Zic8xzTczgc7NmzB5vNRnh4OF5eXtptKG7NMAzy8vI4ceIEdrud+vXrc+p4CrPGTiR97QYALMHBXDa4L1fc2EmHpkTK0N4TGTz53WbWHyy8zEz72FDG9mpB7aolu3yMuyrJmByVnAt8kfLy8khKSjL1YoYi5c3Pz4/q1aqxcPoSdk7+Fkt2FlgsVO3amVv/fRf+gbrgoEh5sDsMpvx2gNfn7iI7346Pp5Vh1zdkYMfYSn+6uUpOMRTni2QYBgUFBdjtuuaIuD+bzcaxI8eZ9er/yNu5CwDPqEiuefJ+Grc9/zW/RKRsHTqVxdPfb2bl3lMAtI6pwpu3t6ROeIDJycyjklMMJfkiibg7w25n/uQf2fnFDIy8XCw2G3X+70a6P3gbHl46nVXETIZh8NXaw7w6ZwcZuQX4etp4pkdj7m4fUymHUqjkFINKjkiho/sTmfPSh2Tu2QeAV2xtejz3ILUalt01rkSk5I6mZDP8u02sSCjcq9O5YTiv92pBRFDlOsNRJacYVHKksjMcDuZPmcOuKdNx5OVh9/Si7p230PPem7DZdK0pEVfkcBhMXnmAsXN3klvgoIqfJ6Nvbc4NzSPNjlZuVHKKQSVHKrOkQ8eY+dKH5Ow8O6lfbCw3jnyI2Ho1TU4mIsWx51g6//4mnq2JaQDcGleDF29uWilmS1bJKQaVHKmMDIeDX79ZwLZPpkFuDnh4EHvHzdz4wK3aeyNSweQVOHh3wR4+WJyAw4DoUF/e7RNHXIx7XztOJacYVHKkskk9lcp3L35I+obNAHjViqb78w/pelMiFdz6g2d44quNHDmTjYfVwvBuDbn/yjpuO4GgSk4xqORIZbJpyQaWvv4x9tQ0sFqJue0Gbn70DmweNrOjiUgpSM3O55kftjBncxIAVzUI563bWxIe6G1ystKnklMMKjlSGeTn5TPzrc858tMCMAyM8HC6PvsQzTTvjYjbMQyDr9ceZtTsbeTkO6ga4M34O1rRqX5Vs6OVKpWcYlDJEXd3eM8hZr/wHnmHjwAQ2KkDdzx3HwEBuuCfiDvbfSydR7/cwO5jGVgs8HDnugzt2tBtZkpWySkGlRxxV4bDwYKvfmHr/6ZBXh6Gjy+tHulHl1uvNjuaiJSTnHw7L/+4nS9WHwLgirphvHtnHFUDKv7hK5WcYlDJEXeUmZbJ1yM/JO3sRTU969bhlpcepUat6iYnExEzzN50lKembyYrz061IG/e79uatrVDzY71j5Tk/VvnjIq4ib1b9zKx/wjS1m7AsFqJ/r8bGfzpiyo4IpVYz5ZRzHq0I/UiAjiWlkufj1fx6fL9VJb9G9qToz054gYWfDmXrR9PwyjIxxIUTJdnH6JFxxZmxxIRF5GZW8BT0zfz49mzr3o0j2Ts/7UgwNvD5GQlV5L374r36kTEKSczm29f+ohTK9YA4N2oAb1HP0FYhHtPBiYiJePv7cF/74yjba0QXpmzgzlbktiRnMYn/dpS142vaK7DVSIV1OHdB/l0wLOcWrEGw2Ih6rYbeeCjF1RwROS8LBYLAzrG8vWDHage5MO+E5nc8t4KFu08bna0MqOSI1IBLZu1lO8fGkV+UjKOgECueHEovYf21aUZRORvtakVwqzHOtK2VgjpuQUM+mwtHy7e65bjdPQbUaQCKcgv4ItXP2X96xMwcnPxqBNL309fof01bcyOJiIVSESgD1/efzl3tovBMGDs3J08/lU82Xl2s6OVKo3JEakgzpxI4aunxpG7OwGAiOu7cPuIgXh66r+xiJScl4eVMbc1p2lUEKNmbWP2pqPsO5HBx/3aUqOKe0waqj05IhVAQvwupg56prDgeHnReuj99H3hfhUcEfnH7r68Fl/c154wfy+2HU3jpv8uZ/3B02bHKhUqOSIubtn3i5jz79HYz6RAWCg93n2eq27rYnYsEXEj7euEMeuxTjSNCuJUZh53frKamfGJZsf6x1RyRFyUw+7gm9c/Y/24TzDy8/Fs3JD+E0dTv1lds6OJiBuqUcWXbwd34Pom1cgrcPDEV/GM/3V3hR6QrJIj4oLS0zL5+PHXODprHgBVu13Dgx88S0iYJq4UkbLj5+XBhLvb8OBVdQAY/+sehnwdT05+xRyQrAP6Ii5m/76jzHjqLSxJSWCz0ezBu7mubzezY4lIJWG1WhhxQ2Niq/rz3IytzIw/ypEz2Xx8TxvCKtgFPrUnR8SFrFu+me8Hj8SSlIQtMIDrxz6pgiMipujTLobPBrUj0MeD9QfPcMsHK0g4nmF2rBJRyRFxEfO++oWlz76JLSsTn5pR9Pn4JZpc3tzsWCJSiXWsV5UfHu5ITKgfh09n838TVlaoM69UckRM5rA7+HrsZHa8NwmrvYDgVs0Y8L+XCI/W1cNFxHz1IgL44eEraBldhZSsfPp+spq5W5PNjlUsKjkiJsrJzGbykNdJmj0fgBo9rqP/+KfwCfAzOZmIyB/CAryZdn97rm0UQW6Bg4e+WM+U3w6YHetvqeSImORk8ikm3fcCaRs3g81Gk4f6cfuIQVg9bGZHExE5h5+XBx/d08Z5KYgXZm5j7NydLn2Kuc6uEjHBvh0H+PHJN3CcOYPh68dVLzxGmytbmh1LROSiPGxWRt/ajKhgH976ZTcfLt5LcmoOY3u1wMvD9fabqOSIlLN1yzax7MV3seRkYwkL4+Y3hhPbIMbsWCIixWKxWHjs2vpUD/bh6e+38MPGRE5n5vHh3a3x83KtWuF6tUvEjS2asZRlz72FJScbn1ox3PPJSyo4IlIh3d42mk/7t8XX08aS3Se459M1pGblmx2rCJcoOUuXLqVnz55ERUVhsViYMWPGRbdfvHgxFovlnNvOnTvLJ7DIJZj9vxnEv/URFnsBQS2bMuDjkYRGhJgdS0TkknVuGMHU+9oTdHYunTs+/o3jaTlmx3JyiZKTmZlJy5Ytee+990p0v127dpGUlOS81a9fv4wSilw6w+HgqzGT2Dv5GyyGQfjVV9D/nafx8fc1O5qIyD/WplYI3wzuQHigNzuT0/m/Cb9x6FSW2bEAFxmT0717d7p3717i+0VERFClSpXSDyRSSvLz8pn6zHukrloLQK3/68Etj9+JxeoSf1+IiJSKRtWDmD74Cu7+dDWHTmfxfxNWMuXedjSqbu719ir0b9q4uDgiIyO59tprWbRo0UW3zc3NJS0trchNpCxlpWcy6ZFXSV21FsNqpclD93DrkLtUcETELcWE+fHd4A40rBbI8fRcek/4jfUHz5iaqUL+to2MjOTjjz9m+vTpfP/99zRs2JBrr72WpUuXXvA+Y8aMITg42HmLjo4ux8RS2Zw+kcLkB14ka8duDE9P2o94hOvvKvneShGRiiQiyIdvHuxA65gqpOUU8NrPO0ydR8diuNgsPhaLhR9++IFbbrmlRPfr2bMnFouFWbNmnffzubm55ObmOpfT0tKIjo4mNTWVoCBzd6eJezl66BjfDhmNcfwE+Plxzcv/pkX7pmbHEhEpN1l5BbwyZwdDuzagailfuTwtLY3g4OBivX+7xJic0nD55ZczderUC37e29sbb++KdYl4qXgOJhzhh3+PgTNnsAYHc+MbT1KnSazZsUREypWflwejbzX/AsNuU3I2btxIZGSk2TGkEtuzdS8/Dn8dS3o6tqph3P7OM1SvpZ9JERGzuETJycjIICEhwbm8f/9+4uPjCQ0NJSYmhhEjRpCYmMiUKVMAGD9+PLVr16Zp06bk5eUxdepUpk+fzvTp0816CVLJbVu7nV+eeQtLdjYekdW5891nCIusanYsEZFKzSVKzrp16+jSpYtzeejQoQD079+fyZMnk5SUxKFDh5yfz8vLY9iwYSQmJuLr60vTpk2ZM2cON9xwQ7lnF4lfsoElL/4X8nLxrBXD3e8+Q3CYxnmJiJjN5QYel5eSDFwSuZDVP6/kt7EToKAA7/p16ffuCPwD/cyOJSLitirlwGOR8rbsh8VsePt/4HDg07wpA8YNw8dXg9tFRFyFSo7IJVjwzS9s/u9nWAwHfpe1ZsBrT+Dl7Wl2LBER+ROVHJES+uWrX9j6/mQshkHIFZdx9+jHsXnYzI4lIiJ/oZIjUgLzv5jLtgmfYzEMql51OX1ffhSrrUJOHC4i4vZUckSKaf7nP7Ht4y+wGAZhV3fkrpcf0nWoRERcmEqOSDHM/+xHtn/yJRagapdO3PXiYBUcEREXp5Ij8jfmT5rF9k+/AqDqdVdz1wv3q+CIiFQAKjkiFzFv8mx2nC04YV07c9fz96ngiIhUECo5Ihfw82dz2PW/aQCEXd+Fu5+7VwVHRKQCUckROY/5X85j59kxOBFdr+ZOFRwRkQpHv7VF/mLR9IVs/3AKFgzCr+nEnc9rDI6ISEWk39wif7Ji1lLi35kIhkHolZfTd5TOohIRqaj021vkrDXzVrH2rU+wOBwEt2/L3a88qoIjIlKB6Te4CLBh4Tp+G/0+2O34t2nJPa89rpmMRUQqOP0Wl0pv829bWPbyfzHsdnyaN6P/60Px8NSYfBGRik4lRyq1HfF7WPD8eIz8fLwbNWDA28N0NXERETehkiOV1r7dh/j5qdex5GTjHVuL/uOfxMfHy+xYIiJSSlRypFI6evg4M4eOxZqZiXeNSO5+dwR+AX5mxxIRkVKkkiOVzumTqXzzxGgsKWfwqBrGHe88Q2BIkNmxRESklKnkSKWSkZ7Fl4+PhuPHsQUF0Wv8CEKrh5kdS0REyoBKjlQauTm5fD7kdQoOHcbi68uNbzxFZO0os2OJiEgZUcmRSsFeYGfq8LfJ3bUbvLy47tWhxDaNNTuWiIiUIZUccXuGw8GXz79P+sbNYLPR8dlHadquqdmxRESkjKnkiNv79q2pnFq2CsNiIW7IfVx2bVuzI4mISDlQyRG39tOUORydOReARvf24epbrzY5kYiIlBeVHHFbK3/+jV2fTAOg5k3d6D6gp8mJRESkPKnkiFvatnYHq1//CIvhIKTDZfQado/ZkUREpJyp5IjbObw3kV+efxtLfh6+jRrQ99VHsVj1oy4iUtnoN7+4lTMnUvhh2OuQkYFHZCR3jRuOp5cuuCkiUhmp5IjbyMnKYdq/x+I4cQJrlSr0Hv80AUH+ZscSERGTqOSIW7AX2Jn65NvkHTgIPj70eG0YETXCzY4lIiImUskRt/D16IlkxG8Bm42rX3iMus3qmB1JRERMppIjFd7Pn/3I8fmLAGjxaH/iroozOZGIiLgClRyp0NYuWMfOT78CIPq2Hlxz+3UmJxIREVehkiMV1r4dB1gx5gMsDgdBl7XmtiF3mh1JRERciEqOVEhnTqTw49NvQU4OXrVr0Xf0Y5oLR0REitC7glQ4eTl5TBv2Bo5Tp7CEhND7reH4+HqbHUtERFyMSo5UKIbDwRfPv0/e3v3g7c2NY/5D1WqhZscSEREXpJIjFcqM978l9be1GBYrHZ4erFPFRUTkglRypMJY8dNKDnw9C4AGA26nfdf2JicSERFXppIjFcLe7QdY+9YnWDAI69yJHvfebHYkERFxcSo54vJST6fz4zNvQW4u3nVj6fPC/WZHEhGRCsAlSs7SpUvp2bMnUVFRWCwWZsyYUez7rlixAg8PD1q1alVm+cQ89gI7055+G+PkqcKLbr7+H11VXEREisUlSk5mZiYtW7bkvffeK9H9UlNT6devH9dee20ZJROzff3GZ+Rs3wkeHnR76XHCdCaViIgUk4fZAQC6d+9O9+7dS3y/Bx98kL59+2Kz2Uq090cqhoXTF3J8zq8AtHqkHw1bNzI5kYiIVCQusSfnUkyaNIm9e/cycuTIYm2fm5tLWlpakZu4rp0bd7Ppv58BENn9WjrrmlQiIlJCFbLk7Nmzh6effpovvvgCD4/i7YwaM2YMwcHBzlt0dHQZp5RLdfL4GeY+Px5LQT5+jRvwf08NMDuSiIhUQBWu5Njtdvr27cuLL75IgwYNin2/ESNGkJqa6rwdPny4DFPKpSrIL+Cbp8ZBSgoeYaHcMXYoNg+b2bFERKQCcokxOSWRnp7OunXr2LhxI48++igADocDwzDw8PBg/vz5XHPNNefcz9vbG29vXd/I1X03dhJ5e/aClxc9Rg8lODTI7EgiIlJBVbiSExQUxJYtW4qs++CDD1i4cCHfffcdsbGxJiWTf2rZzCUkz10EQNwTg4htqks2iIjIpXOJkpORkUFCQoJzef/+/cTHxxMaGkpMTAwjRowgMTGRKVOmYLVaadasWZH7R0RE4OPjc856qTj27zzAhncmAxDe7RquvvkqcwOJiEiF5xIlZ926dXTp0sW5PHToUAD69+/P5MmTSUpK4tChQ2bFkzKWmZ7Fj8+9g5GXi1fdWHo/PcDsSCIi4gYshmEYZocwQ1paGsHBwaSmphIUpHEfZjEcDj79zzgy1m7AEhhI34mvEh5Z1exYIiLiokry/l3hzq4S9/LjxFlkrN0AViudn31YBUdEREqNSo6YZtOqrSRMmQ5AvTtvoWWnliYnEhERd6KSI6Y4fSKFRS+/j8VhJyiuBT0evM3sSCIi4mZUcqTcOQrsfPvsO5Caiq1qGL1feRSLVT+KIiJSuvTOIuVu9gffkL19F4aHJ91eeoKA4ACzI4mIiBtSyZFytWXFJvZ/OweAhvf2oUGLeiYnEhERd6WSI+XmzIkUFo6eAIaDgMta0/2ubmZHEhERN6aSI+XCYXfwzXP/xUhNxVq1Kn1efEjjcEREpEzpXUbKxayPvyd72w4MDw+6jXqUgCB/syOJiIibU8mRMrdl7Xb2T5sJQKN7etGwVQOTE4mISGWgkiNlKi01gwUvf4DFYScwrgX/GtjT7EgiIlJJqORI2TEMpo/6EE6fxhISwu0vPaxxOCIiUm70jiNlZuFX80hduxHDaqXLM4MJCtGFUEVEpPyo5EiZOJxwmC2ffA1A9G030KJDc5MTiYhIZaOSI6UuPy+fmSPfx8jLxatOLLc+eofZkUREpBJSyZFS9/34Lyk4eAh8fLjl5cewedjMjiQiIpWQSo6Uqg3LN3N09nwA4h66m6ha1U1OJCIilZVKjpSa1DMZLHntIyyGQcjlbbm61zVmRxIRkUpMJUdKzbcvf4Ql5QweoSH83/MPmB1HREQqOZUcKRWLpi8iY816sFjpPGIw/sEBZkcSEZFKTiVH/rHjh4+x+cOpANTo2ZVmOl1cRERcgEqO/COGw8GMUR9g5GRji67JLUP6mh1JREQEUMmRf+jnSbPJ2rUHw8OTHiMfxtPL0+xIIiIigEqO/AMHdx9i99TvAajb91bqNKptbiAREZE/UcmRS1KQX8DsFz+E/Hx8GtTlxntvMjuSiIhIESo5cklmTJhOwcGD4O3NzSMfxmrTj5KIiLgWvTNJiSVs3c/h7+YA0GJQbyJrRZqcSERE5FylWnL27t3L3r17S/MhxcUU5Bcw95UPsdgL8G/amC53djM7koiIyHl5lMaDnD59mn79+lGvXj0MwyAhIYGpU6cSEhJSGg8vLmTW+99QcORI4cU3XxiMxaqdgSIi4ppKpeSMGDGCoUOHcs01hdcqWrBgAcOGDePTTz8tjYcXF7F3SwKHv/8ZgGb330l4jXCTE4mIiFxYqfwZvm3bNq655hpeeuklTp8+zbXXXsuOHTtK46HFRRTkFzDv1Y8wHHZ8mjfj2tuvNTuSiIjIRZXqsQbDMJwfOxyO0nxoMdnsCdPJO5KI4ePLbc8/oMNUIiLi8krlcNUVV1zBl19+yciRIwGYNm0aHTt2LI2HFhewd8cBDnw3BwvQbFBvIqKqmh1JRETkb5VKyXnxxRd59NFHmT59OgDBwcF88MEHpfHQYjK73cGcVz/GYi/Ar3EDruvT1exIIiIixVKikjNq1Chat25NmzZtqFGjhnO9r68vn376KTk5ORiGga+vb6kHFXPMmTwbx4EDWLy8uem5B3WYSkREKowSlZyXXnoJi8UCQNWqVWnTpg2tW7emdevWtG3blpiYmDIJKeY4ejCZvV/MwAI0uOtmqmvSPxERqUBKVHIuu+wykpKSGDhwINWrV2fDhg389NNPvPHGGxQUFBASEkLr1q2ZP39+WeWVcmI4HMwe/QmWvFy8Y2vRrX9PsyOJiIiUSIlKzurVq5k8eTLPPPMMcXFxvP322zRo0ID8/Hw2b97Mhg0b2LhxY1lllXK05PtFZG/bATYb3UY8gNXDZnYkERGREinxAIsBAwawe/dumjZtStu2bRk+fDi5ubm0adOG+++/XwOO3UDKyRQ2f/IVADVu6kadJrEmJxIRESm5SxpFGhAQwOuvv8769evZuXMn9erVY+LEiaWdTUwy8/VJODIzsVSvzk2P9jY7joiIyCW55FNl8vPzyc7Opk+fPsTExHD//fdz+vTp0swmJli/eANnVq7FsFjoMmwg3t5eZkcSERG5JCUak/Pqq6+yZcsWtmzZwu7du/H396dFixa0b9+eBx98kODg4LLKKeUgJyePZe98BkDVq66gxeXNTU4kIiJy6Uq0J+f5559nzZo13HHHHWzfvp2UlBSWLl3Kf//7X+69915stksbnLp06VJ69uxJVFQUFouFGTNmXHT75cuX07FjR8LCwvD19aVRo0a8/fbbl/Tc8odZH34HJ05gCwzglmH3mB1HRETkHynRnpxOnTqxadMmRo0axdixY2nRooVzrpw2bdrQrFmzSyo6mZmZtGzZkoEDB9KrV6+/3d7f359HH32UFi1a4O/vz/Lly3nwwQfx9/fngQceKPHzCxxKOELizLlYgBb39iYwJMjsSCIiIv+IxfjzVTWLac+ePaxfv54NGzawfv16Nm7cSEpKCt7e3jRv3pw1a9ZceiCLhR9++IFbbrmlRPe77bbb8Pf35/PPPy/W9mlpaQQHB5OamkpQUOV+QzccDiYMfoXc7TvxbVCPB/43SjMbi4iISyrJ+/clXbuqfv361K9fnz59+jjX7d+/n3Xr1pkyT87GjRtZuXIlr7zyygW3yc3NJTc317mclpZWHtEqhKUzl5G7fSfYbPzrqXtVcERExC2UygU6AWJjY4mNjeX2228vrYf8WzVr1uTEiRMUFBQwatQo7rvvvgtuO2bMGF588cVyy1ZRZKRlsvGTaQDUuLErtRrWMjmRiIhI6ajQf7IvW7aMdevWMWHCBMaPH8+0adMuuO2IESNITU113g4fPlyOSV3XzPFfQFoa1rAwej6iOXFERMR9lNqeHDPExhbOxNu8eXOOHTvGqFGjuPPOO8+7rbe3N97e3uUZz+Xtit/DiV+WAtD+kbvx8fMxOZGIiEjpqdB7cv7MMIwiY27k4gyHg1/HTQbDgX/rlrS/vr3ZkUREREqVS+zJycjIICEhwbm8f/9+4uPjCQ0NJSYmhhEjRpCYmMiUKVMAeP/994mJiaFRo0ZA4bw5b775Jo899pgp+SuiRd8vIn/ffvD0pOew/mbHERERKXUuUXLWrVtHly5dnMtDhw4FoH///kyePJmkpCQOHTrk/LzD4WDEiBHs378fDw8P6taty2uvvcaDDz5Y7tkrosz0LDZP/BaAmJu7UT2musmJRERESt8lzZPjDirzPDlfvTaZ5B/nYw0N5YGv38LHV2OVRESkYijJ+7fbjMmR4jm05zDJPy8AoM0DfVRwRETEbankVDJzx30Gdjs+jRpwxQ1XmB1HRESkzKjkVCKr5q8ma8t2sFq5fmh/zWwsIiJuTe9ylURebj5rPvwSgGrXXUWdJrEmJxIRESlbKjmVxLxJs3CcOIElIIAej59/wkQRERF3opJTCaScSmXv9J8AaNT3ZoKqBJqcSEREpOyp5FQCc977GrKzsUVW57q+/zI7joiISLlQyXFzB/cc5viCZQB0GHwnNg+byYlERETKh0qOm5v3zlQsDjt+TRrS9trLzI4jIiJSblRy3NiG5ZvIit8CFivXPX6P2XFERETKlUqOm3LYHax4/wsAwq5qT51mdUxOJCIiUr5UctzUoukLsB8+gsXLm+6P9TU7joiISLlTyXFDuTm5bJ0yA4CYm7tStXqYuYFERERMoJLjhuZPmo2RcgaCgvjXfbeZHUdERMQUKjluJj0lnX0/zAOg0Z034evvY3IiERERc6jkuJm5H3+PkZWJNTyc6/pcb3YcERER06jkuJETyadJnLsIgLj+t+Hh6WFyIhEREfOo5LiRnz/8BvLy8IyuSaebrjQ7joiIiKlUctzEoX2JnFq0AoDL778di1XfWhERqdz0Tugmfnn/68LLNzSsT5subc2OIyIiYjqVHDewb/t+0tasB+DKwXeAxWJyIhEREfOp5LiBRR9+jcUw8G/RjMaXNTE7joiIiEtQyang9sTvJn3jZgwsXP3QHWbHERERcRkqORXc4k++AyCgbSsaNK9rchoRERHXoZJTge3cuJvMTVsxLBauvf//zI4jIiLiUlRyKrDFHxfuxQluG0edprEmpxEREXEtKjkV1Nb1u8jZshUsFq55oJfZcURERFyOSk4Ftfx/0wEIuSyO2o21F0dEROSvVHIqoJ2b9jj34lx9v/biiIiInI9KTgW04uxenKA2rbQXR0RE5AJUciqY/Vv3OufFueo+7cURERG5EJWcCmbJxB8A8ItrTr1mdUxOIyIi4rpUciqQw3uPkLI2HoCrBtxqbhgREREXp5JTgSyYOAMMB75NGtK4TUOz44iIiLg0lZwK4ljSKc4sXwNAh3tuMjmNiIiI61PJqSAWTp6NxV6AT+0YWnRqZXYcERERl6eSUwGkp2Zw7NelALS4owdYLCYnEhERcX0qORXAwi9+htwcrOHhXH5DR7PjiIiIVAgqOS4uLyePgz8uAKDRbddjtelbJiIiUhx6x3RxS6cvxJGWBoGBXH37dWbHERERqTBUclyY4XCwY8Z8AKJvuAZvH2+TE4mIiFQcKjkubO2iDdiTkrF4enJd3+5mxxEREalQXKLkLF26lJ49exIVFYXFYmHGjBkX3f7777+na9euhIeHExQURIcOHZg3b175hC1H67+dC0DVK9sTHBZkchoREZGKxSVKTmZmJi1btuS9994r1vZLly6la9eu/PTTT6xfv54uXbrQs2dPNm7cWMZJy8/BPYfJ2boDgE539TA5jYiISMXjYXYAgO7du9O9e/EPx4wfP77I8ujRo5k5cyazZ88mLi7uvPfJzc0lNzfXuZyWlnZJWcvLsi9+woKBX5NG1GpYy+w4IiIiFY5L7Mn5pxwOB+np6YSGhl5wmzFjxhAcHOy8RUdHl2PCkklPzeDkslUAtLq9m8lpREREKia3KDlvvfUWmZmZ9O7d+4LbjBgxgtTUVOft8OHD5ZiwZBZ//Qvk5mIND+eyay8zO46IiEiF5BKHq/6JadOmMWrUKGbOnElERMQFt/P29sbb2/VPwXbYHeyfsxCAujddh8XqFj1URESk3FXokvP1119z77338u2333Ldde4xUd6aBWtxnDoFPj50/r9rzY4jIiJSYVXY3QTTpk1jwIABfPnll/To4T5nH22aUXgJh/BO7fAP9DM5jYiISMXlEntyMjIySEhIcC7v37+f+Ph4QkNDiYmJYcSIESQmJjJlyhSgsOD069ePd955h8svv5zk5GQAfH19CQ4ONuU1lIbkI8fJ2rIdC9ChtwYci4iI/BMusSdn3bp1xMXFOU//Hjp0KHFxcbzwwgsAJCUlcejQIef2H330EQUFBTzyyCNERkY6b0888YQp+UvL0q/nYzEceNeNpU6TWLPjiIiIVGgusSenc+fOGIZxwc9Pnjy5yPLixYvLNpAJCvILOLpwBQCNenQ2N4yIiIgbcIk9OQKrf1kDqalY/Py54qarzI4jIiJS4ankuIitMwsHHFe7sp2uNi4iIlIKVHJcQNLBZLK27wTgit7Xm5xGRETEPajkuIDffliIxTDwqluHGF2nSkREpFSo5JjMcDg4vPg3AOp3u9LkNCIiIu5DJcdkm1dtwzh5Cry8uKKnSo6IiEhpUckxWfzsRQBUadNSMxyLiIiUIpUcE2Vn5XJ6TTwALW+82twwIiIibkYlx0Qrf1qBJTcHW0gVWnZqZXYcERERt6KSY6Ld85cDEHlle6w2fStERERKk95ZTXIq+RQ5O3cDcNnNnc0NIyIi4oZUckyy5qcVWBwObNE1qaW5cUREREqdSo5J9i1ZA0DMlZeZnERERMQ9qeSYIDnxBHn79gPQrkcnk9OIiIi4J5UcE6z+cTkWw8C7dgyRtSLNjiMiIuKWVHJMcHBZ4aGq2le1MzmJiIiI+1LJKWdJh45hP3AIgMtu0KEqERGRsqKSU842zFuJBQPv2FpUrRlhdhwRERG3pZJTzg7/tgGAqA5tTE4iIiLi3lRyylHq6XSyEwrPqmp9fXuT04iIiLg3lZxytPbX1VgcDjyrVyO6XrTZcURERNyaSk45Sli+EYBqbZqbnERERMT9qeSUk7y8AjK37QCg2TU6dVxERKSsqeSUky2/bcaWm4PV358GbRqbHUdERMTtqeSUk91nD1UFNm+M1cNmchoRERH3p5JTTk7EbwegVvsWJicRERGpHFRyysGxpFM4kpIAaHVVa5PTiIiIVA4qOeVg45LCQ1XeNaMIrRZqchoREZHKQSWnHBxctxWAiBaNTE4iIiJSeajklDHDMMjYsRuAOhqPIyIiUm5UcsrYob2J2FJTsFhtNL5ckwCKiIiUF5WcMrZj1RYAvGpF4+Pva3IaERGRykMlp4wd2Vx4qCq8ST2Tk4iIiFQuKjllyDAMUnbtBSA2TrMci4iIlCeVnDJ09FgKtlOnsFosNG6nkiMiIlKeVHLK0J71O7Bg4BVRFb/QKmbHERERqVRUcspQ4tY9AATUizU5iYiISOWjklOGTu85AEC1xnXMDSIiIlIJqeSUEbvDIPtQIgCxzeqanEZERKTyUckpI/sPHscrKwMPq5XaTbQnR0REpLyp5JSR/dv3AeBXrSoefpoEUEREpLy5RMlZunQpPXv2JCoqCovFwowZMy66fVJSEn379qVhw4ZYrVaGDBlSLjlL4vjuAwD4xtQwN4iIiEgl5RIlJzMzk5YtW/Lee+8Va/vc3FzCw8N59tlnadmyZRmnuzRnDh4FIKR2TZOTiIiIVE4eZgcA6N69O927dy/29rVr1+add94BYOLEiWUV6x/JPHoMb6B6XZUcERERM7hEySkPubm55ObmOpfT0tLK7Lmy8+zYT54CIKa+So6IiIgZXOJwVXkYM2YMwcHBzlt0dHSZPde+xNP45WXj62kjNDqyzJ6nIsgrcHD4dBYn0nMxDMPsOCIiUolUmj05I0aMYOjQoc7ltLS0Mis6yQcKx+P4VAkEH58yeQ5XZncYzNqUyEdL9rEzOd25PtjXk3s7xfLAVXXw8bSZmFBERCqDSlNyvL298fb2LpfnOn34GABeVcPK5flcya7kdHq8u4wCx7l7bVKz8xn3y25+3prMtPvbU8XPy4SEIiJSWVSaw1Xl6cyxwvE4fuGhJicpX79uP0a38UuLFBxPm4WO9cKoGvBHwdyRlMYjX27Q4SsRESlTLrEnJyMjg4SEBOfy/v37iY+PJzQ0lJiYGEaMGEFiYiJTpkxxbhMfH++874kTJ4iPj8fLy4smTZqUd/xzpJ44gycQHB5idpRys2zPCe6bss653DqmCq/1akGDaoEAGIbBzPijDPk6HoAVCaeYty2ZfzWr3GOWRESk7LhEyVm3bh1dunRxLv8+dqZ///5MnjyZpKQkDh06VOQ+cXFxzo/Xr1/Pl19+Sa1atThw4EC5ZL6YjFOphAAhEZVjT05yag73fLrGufx090Y8eFUdLBaLc53FYuGWuBqczszjpR+3AzBx+QGVHBERKTMuUXI6d+580UMXkydPPmedKx/qyD97enpl2ZMz/LtNzo8fu6Yeg6++8AVJ7768lrPkrDlwmozcAgK8XeLHUERE3IzG5JQyu8PAkZkJQFBYsMlpyl5iSjbL9px0Lg+5rsFFt/fysHJH2z/OatuRVHbzFYmISOWmklPKUrPz8bDbAQgM8jM5TdlbuPN4keUfNx/92/s0igx0fnzkTFapZxIREQGVnFJ3OjMPT3s+3h5WPH3K55R1Mx1Pyymy/MXqQxfY8g/Vgv6YOygjp6DUM4mIiIBKTqk7nZmHh92Or6cNPD3NjlPmGkcGFVnu0jDib+/jYf1jQDJ/GpwsIiJSmjTis5Rl5hXgYdjx9LCCh/t/ef/VtDr/vq4BK/eepGuTagzqGPu398nM+2PvTZCP+3+NRETEHHqHKWV5BQ4MLNislWMPhdVq4Ynr6vPEdfWLfZ8T6X9cKDU80P0P6YmIiDl0uKqU5RU4cFgs2CwWcDjMjuOSth3944yq+hGBF9lSRETk0qnklDJnybGq5FzIzPg/zsCqGqDrV4mISNlQySll+XYHDoutsOQU6Myhv0rLyXd+3KJmcJFZkUVEREqTSk4ps1ktZHt64zCALM0B81fLdv8xceDtf5oUUEREpLSp5JQyb08b2Z7e2B0GnJ35WP7w1do/5tG5un64iUlERMTdqeSUMi+b9WzJcWhPzl9k5BYUuQRETJj7zwgtIiLmUckpZd4eVtK8/Qv35Jw5Y3YclzJ70x8Djl+8qamJSUREpDJQySll/t4epPgGklvggJMn//4OlYRhGIz4fotzuWfLKBPTiIhIZaCSU8rCA7057RtEVp5dJedPVu495fz4hubVCfXXqeMiIlK2VHJKWdUAL077BZPrgPzUNEhNNTuSS3j0yw3Oj4d3a2RiEhERqSxUckpZgLcHNh9vTviHkJVrh8OHzY5kuuV7TnImq3B+nNYxVYit6m9yIhERqQxUckqZxWKhWpAPR4PCCye+27/f7EimKrA7uPvT1c7lcb1bmRdGREQqFZWcMlA3PID9IVGczsqDnTsr9eUdvlzzx7w4t8XVoLb24oiISDlRySkD9SICOBIcwYk8S+GEgIcO/f2d3FByag4vzNzmXH7+xiYmphERkcpGJacM1IsIwGG1sSno7GnS69aZG8gEhmEUOUw15rbmhOiMKhERKUcqOWWgUfVAAOb4RGMYBmzfDmlpJqcqX1+vPUzC8QwAIgK96XOZrlMlIiLlSyWnDDSODMLPy8Y+j0CSQ6sXjslZvNjsWOVm97F0nv7TxH/fDb5CVxsXEZFyp5JTBjxtVtrUCgFgRa2WhSs3boSkJBNTlY+svAKuf3upc/mdPq10jSoRETGFSk4ZaVc7FIAFGd7QrBkYBnz/PRQUmJys7DgcBr0/+s25fG2jCG5uVcPERCIiUpmp5JSRqxqEA7B41wmyr70eAgLgxAmYPbuw8LihUbO3sTXxj7FH7/VtbWIaERGp7FRyykiLmsHUDPElO9/O4iOZcOutYLXCpk2wYIHbFZ1paw4x5beDzuVVI67F18tmYiIREansVHLKiMVioUfzSAB+3JIEdetCjx6Fn1y+HH76yW0mCcwrcBS5wviMRzpSPdjHxEQiIiIqOWWqZ8vCeXLmb0vmeFoOtGkDN9wAFgusXQuTJ0NKiqkZ/ynDMBg1+48J/yYOaEur6CrmBRIRETlLJacMNasRTJtaIeTbDaauOnsop1076N0bvL0LZ0J+/31YuhTy880NewkcDoPnZ27ly9WHsFjgv3fGcU2jambHEhERAVRyyty9nWIBmLr6ENl59sKVjRvD4MEQE1NYbhYuhLffhkWLKsykgfl2B0O/iWfqqsKCM/a2Fs49VyIiIq7AYhhuNgK2mNLS0ggODiY1NZWgoKAye54Cu4POby7myJlshnZtwOPX1v/jk4YB27bBr78WPWwVE1NYhGrXhmrVCgcsm8kwICen8DpcWVlkn0nlzRnxbNl9lICCXB5qF8llET6Fea+80tysIiLi1kry/u1RTpkqLQ+blSf/1YjHp23kw8V7ueOyaKoFnR2Ua7EUzqHTpAns2AGrVxcewvr9BuDjU1h0IiIgPByCgyEwsPDm6ws2W+Hj/B3DKNxr9PstJweys/+4/XU5O9tZasjKcg6STsvJZ3b8UQIzcrnSauHGFpHEpjsgHfDStalERMR1aE9OGe/JgcLBub0+XMmGQyn0bBnFf++Mu/DGqamF17ratw8OHoS8vIs/uNUKnp6FBeP3PT6/f0sNA+z2P4rNP5SY5eDLbac46bBhCwhgUNcm1KtTHfz8wN8fwsKgZs1//DwiIiIXUpL3b5Wccig5AJuPpHDrByuxOwzG9W7Jba2LUQYcDjh2DI4fL5xI8ORJSE8vvGVkXPop6J6ehQOffX2L3nx8ii6fLS8OH18mbjjG2AUJ5NsNGkcG8Um/NtQM0eUaRESkfOlwlQtqUbMKQ66tz1u/7Ob5GVuJiwkhtqr/xe9ktUJkZOHtrxyOwr0zeXl/3ByOPw5d/f6vh0dhqfnzrQQXyzyelsN/vt3Esj0nAbiheXXe+L+W+HvrR0dERFyb3qnK0cNd6rEs4SRr9p9m0OS1TH/oCkL9L3Eci9VauDfG27t0Q55lGAazNh3lxdnbOZ2Zh4+nlZE9m9LnsmhdUVxERCoEnUJejmxWC+/1jaNGFV/2n8yk/8Q1pGT9zZgbE+w7kcHdn67mia/iOZ2ZR+PIIH58rBN3totRwRERkQpDJaecRQT6MHngZYT6e7ElMZU+H6/iyJkss2MBkJqdz+tzd/Kv8ctYkXAKbw8rw65vwMxHOlIvItDseCIiIiWigcflNPD4r3YfS6fvJ6s5mZFLqL8X7/WN44q6Vcs9B0BmbgFTVx3kg8V7Sc0uPAvrqgbhvHxzU2qF/c24IRERkXKks6uKweySA5CYks2Dn69ja2IaVgv061Cb/1zfgEAfz3J5/pMZuXy28gBTfjvoLDf1IwIY3q0hXZtU06EpERFxOSo5xeAKJQcgJ9/O8zO28u36IwBUC/Lm8Wvr06t1TXw8baX+fHaHwbI9J/hm3WF+2X6MfHvht792mB+PdKnHba1rYrOq3IiIiGsqyfu3S4zJWbp0KT179iQqKgqLxcKMGTP+9j5LliyhTZs2+Pj4UKdOHSZMmFD2QcuAj6eNN25vyef3tqN2mB/H0nJ59oetdBq7iLd/2U3C8fR//BzZeXYW7TzOiO830370AgZMWstPW5LJtxu0iq7ChLtbs+A/nbm9bbQKjoiIuA2XOIU8MzOTli1bMnDgQHr16vW32+/fv58bbriB+++/n6lTp7JixQoefvhhwsPDi3V/V3Rl/XDmDrmKL1cf4n/L9nE0NYd3FuzhnQV7aFAtgA51wmgVU4XmNYKJquKLn9e53zrDMEjPLeDAyUwSjmewMzmdtQdOszUx1bnHBiDY15Nb42rQu200TaLM24slIiJSllzucJXFYuGHH37glltuueA2Tz31FLNmzWLHjh3OdYMHD2bTpk389ttvxXoeVzlcdT75dgdzNicxMz6R5QknixSU3wX5eBDs54nNYsFqtZCZW8CZzHzy7OefBTky2IdrG0fQrWl12seG4eXhEjvxRERESsTtZzz+7bffuP7664us69atG59++in5+fl4ep47cDc3N5fc3FznclpaWpnnvFSeNiu3xNXglrgapGbls2TPCTYeOsOmwynsSk4nM89OWk4BaTkF571/1QBv6ob7U79aAHHRIbSLDaVmiK8GEouISKVSIUtOcnIy1apVK7KuWrVqFBQUcPLkSSLPcxmEMWPG8OKLL5ZXxFIT7OfJTS2juKlllHNdek4+yak5pOcW4HAY2B0Gfl4ehPh7EuLnpUsuiIiIUEFLDnDOXonfj7pdaG/FiBEjGDp0qHM5LS2N6OjosgtYhgJ9PMvtNHMREZGKqkKWnOrVq5OcnFxk3fHjx/Hw8CAsLOy89/H29sa7jK7zJCIiIq6nQo4+7dChA7/88kuRdfPnz6dt27bnHY8jIiIilY9LlJyMjAzi4+OJj48HCk8Rj4+P59ChQ0DhoaZ+/fo5tx88eDAHDx5k6NCh7Nixg4kTJ/Lpp58ybNgwM+KLiIiIC3KJw1Xr1q2jS5cuzuXfx87079+fyZMnk5SU5Cw8ALGxsfz000/8+9//5v333ycqKop33323ws6RIyIiIqXP5ebJKS+uPE+OiIiInF+Fu6yDiIiISGlTyRERERG3pJIjIiIibkklR0RERNySSo6IiIi4JZUcERERcUsqOSIiIuKWVHJERETELbnEjMdm+H0OxLS0NJOTiIiISHH9/r5dnLmMK23JSU9PByA6OtrkJCIiIlJS6enpBAcHX3SbSntZB4fDwdGjRwkMDMRisZgdxyktLY3o6GgOHz6sy024MH2fXJ++R65P3yPX54rfI8MwSE9PJyoqCqv14qNuKu2eHKvVSs2aNc2OcUFBQUEu8wMlF6bvk+vT98j16Xvk+lzte/R3e3B+p4HHIiIi4pZUckRERMQtqeS4GG9vb0aOHIm3t7fZUeQi9H1yffoeuT59j1xfRf8eVdqBxyIiIuLetCdHRERE3JJKjoiIiLgllRwRERFxSyo5IiIi4pZUclzEqFGjsFgsRW7Vq1c3O5b8RWJiInfffTdhYWH4+fnRqlUr1q9fb3Ys+ZPatWuf83/JYrHwyCOPmB1NziooKOC5554jNjYWX19f6tSpw0svvYTD4TA7mvxJeno6Q4YMoVatWvj6+nLFFVewdu1as2OVSKWd8dgVNW3alF9//dW5bLPZTEwjf3XmzBk6duxIly5d+Pnnn4mIiGDv3r1UqVLF7GjyJ2vXrsVutzuXt27dSteuXbn99ttNTCV/NnbsWCZMmMBnn31G06ZNWbduHQMHDiQ4OJgnnnjC7Hhy1n333cfWrVv5/PPPiYqKYurUqVx33XVs376dGjVqmB2vWHQKuYsYNWoUM2bMID4+3uwocgFPP/00K1asYNmyZWZHkRIYMmQIP/74I3v27HGp69RVZjfeeCPVqlXj008/da7r1asXfn5+fP755yYmk99lZ2cTGBjIzJkz6dGjh3N9q1atuPHGG3nllVdMTFd8OlzlQvbs2UNUVBSxsbH06dOHffv2mR1J/mTWrFm0bduW22+/nYiICOLi4vjkk0/MjiUXkZeXx9SpUxk0aJAKjgvp1KkTCxYsYPfu3QBs2rSJ5cuXc8MNN5icTH5XUFCA3W7Hx8enyHpfX1+WL19uUqqSU8lxEe3bt2fKlCnMmzePTz75hOTkZK644gpOnTpldjQ5a9++fXz44YfUr1+fefPmMXjwYB5//HGmTJlidjS5gBkzZpCSksKAAQPMjiJ/8tRTT3HnnXfSqFEjPD09iYuLY8iQIdx5551mR5OzAgMD6dChAy+//DJHjx7FbrczdepUVq9eTVJSktnxik2Hq1xUZmYmdevW5cknn2To0KFmxxHAy8uLtm3bsnLlSue6xx9/nLVr1/Lbb7+ZmEwupFu3bnh5eTF79myzo8iffPXVVwwfPpw33niDpk2bEh8fz5AhQxg3bhz9+/c3O56ctXfvXgYNGsTSpUux2Wy0bt2aBg0asGHDBrZv3252vGLRwGMX5e/vT/PmzdmzZ4/ZUeSsyMhImjRpUmRd48aNmT59ukmJ5GIOHjzIr7/+yvfff292FPmL4cOH8/TTT9OnTx8AmjdvzsGDBxkzZoxKjgupW7cuS5YsITMzk7S0NCIjI7njjjuIjY01O1qx6XCVi8rNzWXHjh1ERkaaHUXO6tixI7t27Sqybvfu3dSqVcukRHIxkyZNIiIiosigSXENWVlZWK1F335sNptOIXdR/v7+REZGcubMGebNm8fNN99sdqRi054cFzFs2DB69uxJTEwMx48f55VXXiEtLU1/1biQf//731xxxRWMHj2a3r17s2bNGj7++GM+/vhjs6PJXzgcDiZNmkT//v3x8NCvOVfTs2dPXn31VWJiYmjatCkbN25k3LhxDBo0yOxo8ifz5s3DMAwaNmxIQkICw4cPp2HDhgwcONDsaMVniEu44447jMjISMPT09OIiooybrvtNmPbtm1mx5K/mD17ttGsWTPD29vbaNSokfHxxx+bHUnOY968eQZg7Nq1y+woch5paWnGE088YcTExBg+Pj5GnTp1jGeffdbIzc01O5r8yddff23UqVPH8PLyMqpXr2488sgjRkpKitmxSkQDj0VERMQtaUyOiIiIuCWVHBEREXFLKjkiIiLillRyRERExC2p5IiIiIhbUskRERERt6SSIyIiIm5JJUdERETckkqOiIiIuCWVHBEREXFLKjkiUqFdddVVWCwWLBYLXl5eNG7cmC+//NLsWCLiAlRyRKTCMgyD+Ph43nzzTZKSkti1axf/+te/6NevH/v37zc7noiYTCVHRCqsPXv2kJ6ezr/+9S+qV69ObGws9957L3a7nV27dpkdT0RMppIjIhXW+vXrCQkJoUmTJgAcOXKEZ599Fm9vb5o3b25yOhExm4fZAURELtWGDRtITU0lMDAQh8NBdnY2vr6+TJgwgRo1apgdT0RMZjEMwzA7hIjIpbjmmmto1qwZjz/+OCkpKQwbNowOHTowZswYs6OJiAtQyRGRCiskJIQPP/yQPn36ALB9+3aaN29OQkICsbGxAGRlZTFy5EhWrFgBQNOmTXnttdcICwszLbeIlA+NyRGRCmnfvn2kpKTQrFkz57omTZpQr149pk2b5lz36KOP0rJlS1auXMnKlSvp06cP/fr1Q3/fibg/lRwRqZDWr1+Ph4cHDRo0KLK+a9eu/PDDDwBkZ2dz5swZ7r77bkaNGsWoUaO49tprqVu3LgkJCWbEFpFypIHHIlIhbdiwgQYNGuDl5VVkfdeuXfnggw84cuQIoaGhzvWPPvpoeUcUEZNpTI6IuLUBAwbQtWtX7rrrLgAWLFjAm2++yU8//YTFYjE5nYiUJZUcEXFrWVlZPPfcc6xYsQKLxULjxo158803NfBYpBJQyRERERG3pIHHIiIi4pZUckRERMQtqeSIiIiIW1LJEREREbekkiMiIiJuSSVHRERE3JJKjoiIiLgllRwRERFxSyo5IiIi4pZUckRERMQt/T8MYvtEZV4BNgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Rstar_gt, Mstar_gt, label=\"Ground truth\")\n",
    "plt.plot(Rstar, Mstar, c=\"r\", alpha=0.5, label=\"Taylor Map (order 8)\")\n",
    "\n",
    "plt.xlabel(r\"$R_{\\odot}$\")\n",
    "plt.ylabel(r\"$M_{\\odot}$\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusions**: \n",
    "* the use of the high order variational equations to solve the Tolman–Oppenheimer–Volkoff equations result in a relatively large convergence radius (when using the Lindblom form) of the resulting Taylor map.\n",
    "* the resulting speedup in the evaluation of the mass/radius curve seems to be of the order of one order of magnitude.\n",
    "* the Taylor map allows to capture well the maximum allowed stellar mass.\n",
    "\n",
    "We are not sure this is useful, but damn its cool!"
   ]
  }
 ],
 "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": 4
}