{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ODEs with parameters\n",
    "==================\n",
    "\n",
    "The values of numerical constants in heyoka.py can either be specified when constructing an ODE system, or they can be loaded at a later stage when the ODE system is being integrated. The latter type of numerical constant is known as a *parameter*.\n",
    "\n",
    "Let's start by importing heyoka.py and NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import heyoka as hy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will be integrating the pendulum ODE:\n",
    "\n",
    "$$\n",
    "   \\begin{cases}\n",
    "   x^\\prime = v \\\\\n",
    "   v^\\prime = -g \\sin x\n",
    "   \\end{cases},\n",
    "$$\n",
    "\n",
    "where $g$ is the value of the gravitational acceleration ($9.8\\,\\mathrm{m}/\\mathrm{s}^2$ on Earth). Let's first create the symbolic state variables $x$ and $v$, which represent respectively the pendulum's angle and its time derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, v = hy.make_vars(\"x\", \"v\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we don't want to fix $g$ to a specific numerical value, when writing down the ODE system we will be implementing $g$ as a parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_sys = [(x, v),\n",
    "           (v, -hy.par[0] * hy.sin(x))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The syntax ``par[0]`` means that the actual value of $g$ will be the first value (index 0) loaded from the parameter array when integrating the ODE system.\n",
    "\n",
    "Let's now create the integrator object, using as initial conditions $\\left( \\pi/2, 0\\right)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tolerance               : 2.2204460492503131e-16\n",
      "High accuracy           : false\n",
      "Compact mode            : false\n",
      "Taylor order            : 20\n",
      "Dimension               : 2\n",
      "Time                    : 0.0000000000000000\n",
      "State                   : [1.5707963267948966, 0.0000000000000000]\n",
      "Parameters              : [0.0000000000000000]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ta = hy.taylor_adaptive(ode_sys,\n",
    "                       [np.pi/2, 0.])\n",
    "\n",
    "print(ta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the screen output ``Parameters = ...``, heyoka.py detected that ``ode_sys`` contains one parameter, and set its value to zero. We can change the value of the parameter by directly accessing the parameters array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tolerance               : 2.2204460492503131e-16\n",
      "High accuracy           : false\n",
      "Compact mode            : false\n",
      "Taylor order            : 20\n",
      "Dimension               : 2\n",
      "Time                    : 0.0000000000000000\n",
      "State                   : [1.5707963267948966, 0.0000000000000000]\n",
      "Parameters              : [9.8000000000000007]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ta.pars[0] = 9.8\n",
    "\n",
    "print(ta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that it is also possible to directly set the value of the parameters on construction via the ``pars`` keyword argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tolerance               : 2.2204460492503131e-16\n",
      "High accuracy           : false\n",
      "Compact mode            : false\n",
      "Taylor order            : 20\n",
      "Dimension               : 2\n",
      "Time                    : 0.0000000000000000\n",
      "State                   : [1.5707963267948966, 0.0000000000000000]\n",
      "Parameters              : [9.8000000000000007]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "ta = hy.taylor_adaptive(ode_sys,\n",
    "                       [np.pi/2, 0.],\n",
    "                       pars=[9.8])\n",
    "\n",
    "print(ta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now integrate the ODE system for a few time units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_grid = np.linspace(0,10,1000)\n",
    "e_hist = ta.propagate_grid(t_grid)[4][:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now move to Mars, where the gravitational acceleration on the surface is $3.71\\,\\mathrm{m}/\\mathrm{s}^2$ (instead of good ole Earth's $9.8\\,\\mathrm{m}/\\mathrm{s}^2$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reset time and state.\n",
    "ta.time = 0\n",
    "ta.state[:] = [np.pi/2, 0]\n",
    "\n",
    "# Change gravity.\n",
    "ta.pars[0] = 3.71\n",
    "\n",
    "# Integrate again.\n",
    "m_hist = ta.propagate_grid(t_grid)[4][:,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the evolution of $x$ in time to show how a pendulum swings more slowly on Mars:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABtKklEQVR4nO29d5hcV33w/zlTdme2V21f7Upa9W5JLnKRC7hgMPg12E5CCXkxJPACyUsSSAIBQn7pDZLACySYBEI1GIML2NjGNm6SJVla9bIrbe+zdWZ3yvn9ceburlazU2+ZXd/P8+yz0syde8/Oved8z7cLKSU2NjY2NjaL4bB6ADY2NjY22Y0tKGxsbGxs4mILChsbGxubuNiCwsbGxsYmLragsLGxsbGJi8vqARhBRUWFbGpqsnoYNjY2NkuG1157bVBKWRnrvWUpKJqamjhw4IDVw7CxsbFZMgghLiz2nm16srGxsbGJiy0obGxsbGziYgsKGxsbG5u42ILCxsbGxiYutqCwsbGxsYmLLShsbGxsbOJiCwobGxsbm7gsyzwKGxsA/D448wsYuQD55bD6ZihdafWobDLF74PTPwffRcivgDU3Q0mj1aNa1tiCwiYpWrtG+fKz5xgLBLl/TyN3bKmxekiLEwnDq1+FX/4FBCfnXhcO2PFuuPUvIbdQ10tOzYT4t2fOsr99hB0NJXz4pjUUedy6XmOp88sTffzXSxfwuB186IbV7GgsTe0EkTC8/GV45i8hODX3unDAzvfAm/8Scgv0HbSJ9I0F+OenTtM2OMktG6r47b3NOB3C6mEBIJZj46Jdu3bJdDKzX/u7t+LyFNK861aKdt0Lbo8Bo1t6HLo4wv1fe5m8HBeFHhcXhqb4s7ds4H9ft8rqoV1OOAg/egCO/QhaboUb/giqNsNYF+z/OrzyFajcAO/5CRTErFaQMjOhCL/19VfYf2GYzbXFHOseZVNtMd/74FXk5dh7MYDv7b/IHz90lPpSL4FgmFF/kAd/ew9711Qkd4JwEB76HTj+E1h7O1z/h1C9GUY71abglf+n7vN7HlZaxhKjfzzAO/7tRQYnpllVWcCJnjHu3lHHP7xrG0KYIyyEEK9JKXfFes/2UUSZCYbInxmkdvAFin7+UUJf2gUXXrJ6WJYzOR3io989RHl+Lj//+PX88g9u4LZN1fz14yc52Ttm9fAuRUp4+PeUkLjlc/Ab34P6XUrgl6+G2/4KfushGD4P/3UXTI/rctl/ffoMr7YP88/3buen/+davvruXbR2j/I3j5/U5fxLnbP9E3zmJ8e4rqWCX/7fG/jlH+yjqTyfP/j+YcYCwcQnkBJ+/EElJN78Bbj/O1B/Bbhy1X29/W/gN38IQ2ei93XC+D9KR6SUfOIHRxienOEHH7qaxz92HR+/pYUfHerioYNdVg8PsFhQCCH+UwjRL4RoXeT9fUKIUSHE4ejPZ4waS47bxfo/fYnBDx7lAf6M/skw8pt3wslHjbrkkuDBF9vpGPbzz/dtp7IwF5fTwV/dvYUir5vPPnLM6uFdystfhqPfhxv/FK79OMTaia2+Ce7/Hxg4AY98VC1CGdDl8/Pvz57jHTvquGt7HQC3bKzit65cybdeuUjb4GSCMyx//uaJk+S4HPzDO7eR63JSnOfmH961jYHxaf716bOJT/DSv0LrQ3DzZ+Ca/xP7vrbcAvd+G/qPw08/lvF9NZNnTw3w3OkBPnHrOrbWlwDw0Zta2N5Qwt/9/CT+mbC1A8R6jeJB4LYExzwvpdwe/fm80QPaUFvMrW+9n9smP8dI8Ub4wfug/ddGXzYrmZwO8dXnznPLhhXsbiqbfb00P4cP37iGl88Pc+jiiIUjnMfAKXjqz2HdW+C6T8Q/dvVNcNOnlebx+nczuuzXnjsPwCduXXfJ6x+9uYVcl4MvPX0mo/MvdY53j/Hk8T5+59pmVhTNmXK31pfw1m21fPvlC4xOxdEq+o7DU5+FDW+Fa/8g/sVablGbhNYfKsGyRPjXZ87SWJbHe66eC7RwOAR/cscG+sam+f6BDgtHFx2PlReXUj4HDFs5hli8Y0cdtdXVvG/mj5Aljco2Ojlo9bBM52dHuhn1B/ndfasve+/e3Q0Uelz856/bzR/YQiIRpR3k5MNb/wUcSTzWez8O9bvhF38KU+k9gpPTIb5/oIO3ba+lrsR7yXuVhbncvbOOnx3pwTc1k9b5lwPfekU5r3/7mubL3vvg9auZnAkvvhBGIko7yC2CO/8ltiaxkGt/H2p3whOfTPu+mklr1yivXRjhvdc04XZe+tzuaS5jW30x33r5Alb7kq3WKJLhaiHE60KIx4UQmxY7SAjxgBDigBDiwMDAQEYXdDgEv7tvNUeGHBze80/qgfvZ72d0zqXId/d30LKigJ0xolMKcl3cvaOOXxzrZTwZO7ORHH8YOl6GN/1F8g5qh0MJFb8Pnv2rtC77eGsvUzNh7t8TOzTzN/asZCYU4UdZYmc2G/9MmJ8e7uaOzTUU510eAbaxtojtDSU8dLAz9gmO/Qg6X1VRavnlyV3U4YS3fRGmhuC5v89g9Obw/QMdeNwO7rmiPub7v3nVSs70T/DaBWs192wXFAeBlVLKbcCXgIcXO1BK+VUp5S4p5a7KysyjWW7dVE1pnpuvnc1XkTMnHoEzT2V83qVCx/AUhy76uOeK+kWjLu7aUcd0KMITrb0mj24e4ZAKl6zcANt/I7XPVm2Cne+GA99QuRYp8uNDnawsz2PXythhnhtri9hYU8SjR3tSPvdy4OmT/YxPh7hnV+xFEODunXWc7B3nRM+CwIhwUN3Xqs2w9b7ULly9RT0L+78GPuvNNosRiUgeb+3lxnUrKPbGDqW+Y0sNOS4HPzti7TOU1YJCSjkmpZyI/vsxwC2EMCX2zeN2ctf2On55op/JK34XytfAE3+sFqY3AE8e7wPgts3Vix6zo6GEhjIvj1m5EL7+HRg6Czf9mdpNpsr1f6Ti8J/725Q+NuoP8vL5Ye7YUhM3fPH2zdW8dmGE3tFA6mNb4jx1oo+y/ByubF5cG3jLlhqEgJ8fW7DZOPw/Kjrtpj9LzpS4kBs+CQh4Pnu1itcujjAwPs3tcXKSCnJd3LC2ksdbe4hErDM/ZbWgEEJUi+gsFELsQY13yKzr37a5mulQhGfOjcItn1ULUusPzbq8pTx1oo+WFQWsLM9f9BghBDevr+Kl80MEghZEZkQi8Ot/gZrtsP4t6Z2juE5pFUe+D+N9SX/sudMDhCOSWzasiHvc7VuUoH3yuIValwWEwhGePtnPvnWVcZPGygty2d5QwjOn5pmLIxF44Z+Ur2FtoliXRShpgO33w+HvwERmpmijeOxoDzkuBzetj/8M3bGlmr6xaY50jZo0ssuxOjz2O8BLwDohRKcQ4neEEB8SQnwoesg9QKsQ4nXgi8B90kSvzu6mMioKcpRpZf2dSqX91d8ue61idCrIK23DvGljVcJjb1hXSSAY4ZU2CxyHZ59SsfNXfyQ5R+diXPV7ytSx/2tJf+SZk/2U5rnZ3hA/u3h1ZQH1pV6eP/PGCoZ47cIIo/4gt2xI/AzduG4FRzp9DE5MqxfO/BxG2uCaTO/rhyE8rRIts5BfnRrgmtXlFOTGT8q8rkWZ0l84Y53Aszrq6X4pZY2U0i2lrJdS/oeU8itSyq9E3/9XKeUmKeU2KeVVUsoXzRyf0yHYt24FL5wdJCJR6uzwOTj2YzOHYTovnhskHJEJdzoAV68qJ9fl4JmT/SaMbAEv/zsU1sDGuzI7T/lqpZHs/w+YmUp4eCQiefb0APvWrUhYYkEIwXUtFbx0fohQOJLZOJcQz5wawOVQf3siblq/AinVwgmofJiiOtjwtswGUblWZXHv/xoE/ZmdS2d6Rv2cH5zk2iQy0ysKctlYU2TpZiOrTU/ZwN415fimghzvGYN1dyhfxStfsXpYhvJK2zBet5NtDSUJj/W4nVy1qpznzN7tDJ6B88/A7v8NrpzMz3fV74J/WAUtJOBk7zjDkzNJLYIAe9dUMB4IcdRC04HZvNI2xLaGEgqTqHe1saaIioIcXjg7CP0noe1XsOcD4NShVtZVH1IRUCd+lvm5dOTFs8qCfs3q5J6h61oqOHhxhMlpa6wZtqBIwN7ojXzh7KByqu3+AHQdgK7XLB6Zcbx8fogrVpZeFte9GFetKuf8wCRDmunADA7/T7TI32/pc76Ve6FsFRz874SH7m9XZrY9zWUJjlRoi8Gvz74xzE9TMyGOdo4m/f04HILdTWXqez38LXC4VPFGPWi6HkpWwqHE99VMfn1ukLL8HNZXJ1ec8tqWCoJhyatWmHixBUVCVhR5WFtVMDfJt/8G5BTAq8nbs5cSvqkZTvWNc9Wq5CY5wO4mZac/YFasdyQMR76nyoYXLh6VlRJCKKFz4QUYOhf30Ffbhqkt9lBfmpfUqcvyc9hQU8TL57M/AUwPDl7wEYpIrkxSUIDyB/aMTBA+/F3lwNarsJ8juplo+xWMtOtzTh146dwQV68qx5FkddhdK8twOcTsJsVsbEGRBNesruDVtmGC4Qh4imDrvdD6I5Wstcx4tW0YKeHKVUkmOAFb6ovJcTnMSwpqe05Vg001byIR2+5XWsrhby96iJSSV9uH2Z3CIgiws7GE1zt8loY4msWrbUM4BFyxSH5JLHY3lXG94wjOqQFj7isCDi1+X82kZ9RPz2hgdoOVDN4cJxtqijhoUckcW1AkwRUrS5kORTjZE602uv03VTTF8Z9YOzADeO3iCDlOB1vri5P+TK7Lyda6YvN2O4f/BzzFymekJ0W1sOpGOPqDRYvKXRiaYmB8OmmzisbOxlLGp0OcHVhalU3TYX/7CJtqi5PyT2hsqCnkXvfzTLhKoOXN+g6opAGar1eh7VlQLPDwRR8A21Psx6E2G6OWBEXYgiIJdjSWAHCoIyrN63Yqp/aR71k3KIM42jnK+ppCcl2pJa/taiqjtWvU+HyKmSk4+TPY9A5j+oVseofqnNZzOObbhzt8QGq7ZZh7hg5aXIrBaCIRydGuUbY1JL/RAHCF/exzHOJZ5zX6OLEXsvlulcDXe0T/c6fIoQ4fOU4HG2pSa561c2Up/mCYU336lMdPBVtQJEFdiZcVhblzk1wI2HYfXPh1WqUfshVtkm+pS22SA2xvKCEYlpeXYtCbc0+r7mab3mHM+de/RTlTjz0c8+2jXaPkuhysqUytk1pzRT7FXjeHorvJ5Ur70CQT0yG21pWk9sGzT+GR03x3Yocxm431bwXhVCZjizl80cemuqKUN2NazbWDFjxDtqBIAiEEOxpLOBTdTQLKTwHLKlP7wvAU44FQSmYnjS3Rz7QaHQJ64qfgLVVRSkaQVwbNN6hCgzHMFK1do2yoKcKVZESYhvYMWWVjNgstBHhzqpuN448wk1PKS+H1nOw1YMecXw6r9qkcKAvNT8FwhCNdPrYnEXq+kPpSLxUFubOmKzOxBUWS7Gws5cLQ1FwIaEkj1F2hFq5lwpFOHwBbUt0NArXFHkrz3MbmCoRm4NTjyjdhhHlCY9PbVYRMz+uXvByJSI51j6WlcQFsrSvm3MBEVjSiMYqjnUrjaqlKQeMKBuD0E8y03E4Yp3HP0KZ3gO/CZffVTM70TRAIRtISFEIINtUWcazb/HwcW1AkibZjPtY9z7Sy4a3QfSirK1SmQlqTPIoQgs11xbR2GWh6an8Opkczz9hNxPo7lZliQfKdZlZJV1BsrC0iIrHExmwWR7tG2VhblHQODqASJ2cmyN/+vyjNc9PaadBCuPY2QMDpnxtz/iTQTLMba4rS+vym2iLO9k8wHTJ3s2ELiiTRbuzx+Tb49W9Vv5dJu9Qj6UzyeWypK+Z037hxDu0TP1U5LKv2GXN+jbwyWHkNnP7FJS9rO91NdelOciVgjncb7MexiLQ1rpM/g9xiRPP1bK4rNk6jKKhUPdRPP2HM+ZPgRM8YOS4HzRWLF9uMx6baYkIRyZk+c6PnbEGRJCV5OdQWey511lasUX0QloH5KRKRHEvTka2xpU49xKeMsDFLqXaCa242JtppIS1vgr6jMDrXdOhYt5rka6tSi1bRqC/1UpjrssR0YAZt6WhcUqo+L6tvBFfO7GbDsB3z2tug+yCMW1PN92TvOOuqClP2cWlsrFWbFLOfIVtQpMDG2qLLo3o2vBUuvpi1pYyTpXPEz+RMmA1pqsTA7GdP9hqwY+4/DuM9+sfYL4Z2nbNPzr7U2jXK+urCtDUuIQQbaosu1UqXEZqmpC1mSdF7FCZ6Z7/v9TVFhCKStsFJI4Y4V7b8zC/iH2cAUqqowFTDYuezsiyPglzXpSZwE7AFRQpsqCni3MDkpaaV9W8BGVElr5cwp6N283R3ywANZXl43A5OG6EWn4ku2Ktv1v/csahcD8WNc9cFTvWOJ12bZzE21hRxsmec8DLM0D7TN45DqNLqSaMJ4jW3ALAu+vwZopWC6mpYVG+Jn2JgfJqhyZmMNmMOh2BDTaHp5ktbUKTAxpoiwgvtg9VbIX/FkhcUp2YFReqObA2nQ7BmRcGs0NGVs0+ptphFi3cD0xUhlPnp3DMQmmZoQk3yTAQpqGfIHwxzYcigHbOFnO6boKkiH487hfyAM0+pOVSo+lY0V+TjcghjniFQ93XtrSofJ2RiEUvm/Jvrq9MXFADrqgs53TeOia15bEGRChtmHdrz7IMOh7Kbn3taFatbopzuG6e22JNS2YVYrK0q1H+ST4/DxZfV92wmLW+G4CRc+PWsltSSoaDQIsrO9i+/Uh6n+8ZZuyKF78fvg45XlECOojl6T/Ua+P2suVklbXbuN+4aMTgRLQGUbsSTRsuKQsYCIQbGzRN0tqBIgcayPPJznLM3fJY1t6heBt2HLRmXHpzqHWdthmYVUIKib2ya0amgDqOK0vY8RIKz5gnTaL4OHG44/yxn+jPXuADWrIgKimVW8ykQDNM+NJna93P+WZBhWPOmS15eW104+30bQtO1qvjj+V8Zd40YnOwdo6bYQ3FeZpuxlugzdMbEzYYtKFLA4RCsrS683Fm76kZALFnzUygc4fzA5Kx9OBO0c5zWc6Kf+yW486HhKv3OmQw5+VC/G9qe40zfBIW5LqqLMou4KvS4qS7ycNbk8EajOT8wSUSS2maj7Vcq3Ll+1yUvr6sq5OLwFFMzBjXp8RSrZNnzzxpz/kU42z+RsUYKsCYqjA0zz8XAFhQpsqaygHMDC+zL+eWqUOASFRTtQ1PMhCMZ299hzrSi60Pc9rzKa9Cjk12qrLoBel6nq6eblqoCRCY9nKO0VBWYuhs0g7SCIdqeh8arL8uyX1tViJQGm+dW7VPNxwLmhJlGIpLzA5Mp1wiLRWVBLsVet61RZDNrVhQwMD7NqH+BaWXNLarznX/p1fLRI+JJo67ES36Ok9N6Ra1M9MPgKWUusILm60FGKBnYr8v3A+oZOjcwsax6U5zuG8flEDSVJ5lINtYDQ2fU97uAddUGRz6BEhQyDO2/Nu4a8+gZC+APhlm9Ir1Eu/kIIWhZUWCqVmoLihTRQv8u2+0036DCZC+8aMGoMuN03zhCzNnPM0EIQUtVoX4hsu0vqN/N1+lzvlSp24V0edky87ou3w8oZ+TUTJjuUb8u58sGTveNs6oynxxXkktKnPvaWJZHjsth7I65fje480wzP52L/i0phQ7HoaWqgNP95kU+2YIiRbTF4txCZ2T9LnB55ibAEqJtcJLaYi/enNTKHi/GmhUFnB/US1A8DzmFUL1Nn/OliisHX+UurnEc002j0Mxzy8n8dG5gMjVB2v6c8hVUb73sLadD0Fyez/mFJl49ceUqs1fbc8ZdYx7aeqGboFhRiG8qyNDkjC7nS4SlgkII8Z9CiH4hROsi7wshxBeFEGeFEEeEEDvNHuNCGsryyHE6ZncIs7hyoWGPWtiWGG2Dk6yqzFwl1miuyKdvbJrJaR2cke0vwMqrwenK/FzpDqHoCtY5OmnJn9LlfJqd+rJnaIkSDEfoGJ5KrX5R2/OqVLwj9uakuSKfNr02G4ux8hoYOAFTxndmPDcwQZHHRUWBPn42Q3yBcbBao3gQuC3O+7cDLdGfB4AvmzCmuDgdguaK/NiOtqbroLfVlAdPL6RU5RKSti0nwarogpFxGYbxPhg8bZ1/Isphx2YAqkYO6nK+0vwcKgpyTC/sZhSdI35CEZn8MzTaCSNtar4sQnNlPheHp4xt+7nyGvX74svGXSPKuf5JVq/QJxgC5jQTQ7WueVgqKKSUzwHxVtW7gP+SipeBEiGESam5i6M5Iy+j6VpALik/xdDkDOOBUNrVLGPRXKmToLgQNePFWVDM4GV/AwFycHTot6A0V+TTvkyys9uj9znpZ6gtqnXH8Ts1V+QTDEs6Rwz049TuBGeOqtVmMGcHJnSJeNKoLvKQ63LMfvdGY7VGkYg6YH6zh87oa5chhHhACHFACHFgYMDYAn2ro7udy8pp110BLu+S8lOkPMmTQNtZZiwo2l+A3KKYdmwzOTc8TbtnI1x8SbdzrixfPoKiLdVn6OKL4CmBFZsWPUQ3rTQebo+asxf0u6+xGPUHGRifZrVOwRCgcrpWlueZ9gxlu6CIpafFdPNLKb8qpdwlpdxVWVlp6KBWryggIrn8Ji1BP8V5AwSFx+2krsSrg0bxEjRcaal/IhyRXByaYrBsp6p0Oq2PTbipPI++sWnjkspMpH1okkKPi7L8JO3vHa+qeeJYfPnRnsfzRu+YG6+GnsMwY9x1tHmwSsc5BmpD1j6kj98sEdkuKDqBhnn/rwe6LRrLLNpD3D4Y4yY1XQd9S8dP0T44icshqC/16nre5or8zCa536ccjQ1X6jamdOj2+ZkJRwjWXanCnzsP6HLepugzdHHYnIluJG2DkzRX5Cdnf/ePwMBJJSjiUJafQ5HHZY5DOxLS7b7GQisA2aSzoGiuyOfi0JQplYizXVA8ArwnGv10FTAqpeyxelAry7RJHmMhbNqrfutopjCStsFJGsvy0m6kshjNFfm0DUykH+fdFZ24CRYUo9G0xvxVV6n6QDo5PjXznFk2ZiPRBEVSaAtygg2AEILmygJjTU8Qfb6EofP1QnTX31iWp+t5myrymQlH6PYZn49jdXjsd4CXgHVCiE4hxO8IIT4khPhQ9JDHgPPAWeBrwO9ZNNRLKM5zU5rnjq321e5QheQ6XjV/YGmQ0iRPgeaKfMYCIYbTjfPueFUtzHVX6DuwFNEW8sbaalXmXKcFpbFcLRpmmQ6MYjoUpsvnTz7iqeMV1Y88ifu6qsLgXAqI5nJsNjQA5cLQFFVFuamVX0+C2c2GCX4K64y/gJTy/gTvS+DDJg0nJRrL82P3FHB7oWbbkhAUkYikfWiSa9dU6H7u+ZFP5QW5qZ+g41XVZCZXPwdgOrQNTuF1O6kqyoXGq+DQtyEcvKw+UaoUedyU5+cs+b4UF4emkDIFH1fHK1C9RRVcTEBzRT4/PtSFfyasWzJoTBqvgUP/DeGQIf6wi8OTrNQx/FxjzgQ+yXUtxvpls930lLU0lefNqpSX0XCl6ssbMidrMl16xwIEghHdbacAzZlEPkXCykRRb63ZCdRurUmzvzdepfpT9B7R5dxNFfnGm1YMJqWIp3AIOl9L2u+0Mqp1dYwYrHXV71b9KQZOGHL6C0NTrNTZ7ARQVZSL1+2kLZavVGdsQZEmK8vylKMzFCMhqGEPhAIqSiaLaTcoGgOgtsSLENCRThz8wEmYGbfckQ3qO2quiE5yrcy5To7PlfE2G0uE9lQctf3HlKBN0u/UEF1cO4x2+NdHzWAGOLSnZkL0j0/PCj09EcK8EFlbUKTJyvJ8IhI6Y+12tAWu4xVzB5UiF6ITsNGAhzjH5aC22EtnOpNc+94adus7qBQJhSNcHJ6as78X1UJBtSpPrQNN5fn0jAYuz8dZQrQNTlGWn0OxNwlTnGaOTXIDoDl/DY8MK22GvHJDBMXF2Tmm/2YMoombJmiltqBIE22HEHNHWFQDxY1ZLyg6hqdwOQQ1xfqGxmo0lHnTm+Qd+yG/Uk1gC+n2BS4tTSGEKv6oo0YBizxDS4SO4anZnX/ig1+Bwloork/q8PL8HPJynMYLCiGgbtdcpJ2OaPe2yYDNGChh2jniN7xkvS0o0kRzTi3qjGzYk/UO7Y4RP7UlXpwOferPLKShNC9NQfGK8k/oVBcnXTTb+CULYd1OGD6nS56MmVErRtE5MkVDsjk4Ha+oeZHkfRVC0FiWZ7zpCdQGYOAUBMYSH5sCF6OCQgup15v6sjxmwhH6De6fbQuKNKkoULudRcMbG66E8W5VAC1LUbtBY7QJULud/vHp1Ewrk0NqIbY4fwLmbOOXJCPWRdt2dmdeILAp0WYjywlHJF0+f3IaxcQA+C5e1vY0EQ1laW42UqXuCkDqcl/n0z40SbHXnXGf7MXQhLTRDn9bUKSJciTlL/4QawtdFpuf1G7QGJUY5nwfMf04i6FNVIvzJ0BVRXU6BDXF8/pk1+4ABHRlvqAU57kp9LiMLXxnIH1jAYJhmVxWf/ch9bs2tU4BSqPwG9+gp84Yh/bF4SlDHNkaZjn8bUGRAU3xIg6qNqlGRjosKEYwNRNicGImeftyGtSXpuGM7D4ECJWLYjEdI1PUFHsuzVr3FEHlOt0c2vWleUtWUGiLU1Kbje6DqPuaWoHHhlIv/mCYwQmDQ829JVDeott91bgwNGVIDoVGXUlUoxg29hmyBUUGNJbn0Tnsj11rxelWiUVZKii6oouT3jWe5tM4u9tJ4SHuPgQVLWpBtpjOEX/sRbDuCrXz1GGXW1/qTU3jyiI0AZfUZqP7kBKwual1CdS0UlPMT/W7oHO/LvcVVNRcl89vSA6FhieaDGqbnrKYxllHUiD2AbU7oed1lUCWZcR01OpMRUEOXneKUSvdh6LmHevpGJ6KLUjrdsLUoLK5Z0hDVKMwq/exnnSMTCEE1JZ44h8opdowpXFfG83KpQC1AZgc0OW+AvSNTxOOSOoM3IyBeoZs01MWo5lWFjUd1O1UCUYDp0wcVXJou3wjNQohRGohsmM9MN6Tsh3bCALBMP3j07EFqebQ1sFMUV/qZWomzMhUMONzmU3HsD/aQCdBeY2xbpjsT+u+pmW+TBfN0a5TmGxPtFjfJT4uA2goM958aQuKDNDsg4uaDjQHmc6RFHrQMTyFx+2gMp06TCmQUnhjz2H1Ows0ii5fHEE663/SR1BAig7/LKFjZBGNayHa85/GfZ01rZghKFZsUh3vug/rcjrtGdLWCaNoKPXSM+onaGDbWFtQZIA2SboWk+Zlq1WHNp0dZHqgJnmebj18F6M+qhYnZVrpPqQqxlZvMXRMyTDrqI2lUTjdytmui6BIoJVmMZ3DSUbNdR8Ch0tVaU2DRrNCZF05ahPQ87oup+sZVSbpGoMFRX1ZHhGJoeXGbUGRAR63k4qC3MUnucMBtduz0qHdMexPPlEqAxrL8phM1rTSdRAqN0COcX6TZOlM5Oyv2Q49RyCS2S6ubolqFDOhCL1jAeqT8XF1HYQVG1Rl5TQwwwY/S802JSh08Bl1+/wUelwU5BpbpFsT1kZGPtmCIkNU1EqcG1S7E/qOQcjYzMlU6RhJofRCBmgLbcKJLmV2ObJHpnA7BVWFi9iXa7Yp/9PwuYyuU+x1U7QEcyl6Rv1EJIk3G7P3NX2/U12pl96xACEDTSuz1GyHgA98FzI+VbcvYLjZCZhNmjUy8skWFBlSlyi8sW4nRILQ22reoBIw6g8yHggZmmynoe2YE6rFo50qkqh2u+FjSobOET91JV4ci5U30fI8dDBTLMVcirlgiATP0EibWngz2ADUlXiJSFUW33C0+6qDn6Jn1G+4IxugptiLyyEM1bpsQZEh9aVeun2BxYtyaTupLPJTxCxNYRDajqorkaBIM3PXKDoTFburXAfO3DkHfAYsxVyKufDqBM+Qdl/r0r+vtSUJfIF6UrVJ+VN02AB0+/yzYzcSp0NQW+JNr6R/ktiCIkPqS1UuxcDEIqal4npVCTWLIp86Tcih0Cj2usnPcdLtS7Ab1ByeVZsMH1MydI744wtSp1s3x2f9Esyl6BxJsvJw10ElUFdsTPtas1rpqAmCwpWr/CkZbgD8Ub+cGYIC1GajyzY9ZS8JwxuFULtkbWeVBWhmAzNMT0Ko3U6XL8FD3H1ILSZu41X1RExOhxianElsVtHJ8bkUcyk6hpOsPNzzuhKoGbSOrTNTo4BooEJm97Vn1JwcCo26Em9irT0DbEGRIfWzuRRxblLNNhg8DTPZYV7o8vkpyHVR5DWnZXptiTe+RiGl2sFlkX8CkjDN1WyDwGjGjs+lmEvROTKV2FErJfS1Zhzu7HE7Kc/PoSuRVqoXNdtgaiijys/a826WRlFb4qV/fDp2x00dsAVFhsyFN8YTFFtBRqD/uEmjio/mZDM6h0KjrtQb35k91g3+EahOrWCcUWjaT1IaBWRsflqKuRQ9o4HEi+Dsfc08L6au1Ngd8yVojvcM7qtmJqs1qCnYQupKvUgJvaPGCFNbUGRIXo6L8vyc+JNcWwB1SuTJlG5fEpNcR+pKvAxNzuCfWaTmVV80IqwqvYQsvZnbDSYwG6zYqIvjs75saWkUoXCEvrEAdYm+Hx3va21xgs2GnlRtAuHMyE/R7fMjBFQVG1v5QGO2SkQiE2+aWCoohBC3CSFOCSHOCiE+GeP9fUKIUSHE4ejPZ6wYZyISRq2UNIKnJIsEhTnRGBragruoM7L3qPpdlb7DU096RlUfihWL5VBouD0qQTDD+1rkUbkURpeK1ou+8WkiMomMYx3va12ply6zHP5uL1Suz+i+9vgCVBTkJq6DpROaoEgYNJImlgkKIYQT+DfgdmAjcL8QItYT9byUcnv05/OmDjJJEqrFIlqHv/eIeYNahEAwzNDkDLUmOdkA6kqUaWXRHWHfsagwLTZtTPHo8QWoKsxNrkVszTYVc5/hAlZb4p11gGY72n1MuNnQ8b7Wlai+FKY5/DO8r92jflPnWHX0WkY5/K3UKPYAZ6WU56WUM8B3gbssHE/a1JfmJd7tVG+FvuMQtjayRas/Y4VGsehD3NcKVdbXd9LoHvUnX5+nZptKFBzvyeiaCR3+WcSsoEi0EOp4X2tnd8xm+Sm2q4q3ad5Xs7V2j9tJZWFu4ujCNLFSUNQBHfP+3xl9bSFXCyFeF0I8LoRYNMheCPGAEOKAEOLAwMCA3mONS32pl+lQnFwKUCF34WnLS47Plj5OZF/WkaoiDw6xyCQP+mHobNoF44ygZzQwu0NLiE4O7Zpiz5LRKJIqdqfzfa1PJmhETzQHfBoVFaSU9IwGEueY6IyRmw0rBUUsvX7hlvwgsFJKuQ34EvDwYieTUn5VSrlLSrmrsrJSv1EmgRbZEFft01pAWmx+Mqv08XzcTgfVRZ7Y4Y39J1REWJYk2mmTPGmzQfVmQGQsKGpLvIxMBRd3+GcR3T4/RYmK3fUf1/W+mq5RaOPuO5ryR8f8IaZmwomDIXSm3sBcCisFRSfQMO//9UD3/AOklGNSyonovx8D3EKICvOGmBzaQ9wTLzStfA2481TFUQvRdhxJ75h1YtGkuyyLeBqenGEmFEl+N5iTD+Wr5xy3aZLQ4Z9FJBU113dM/dbpvpbmufG6neaFyHqKlX8lDY2iK1kfjs5srS9mVYUx/bmtFBT7gRYhRLMQIge4D3hk/gFCiGoRDfYXQuxBjXfI9JEmYHaSx3uIHU41aSzWKHpG/aZGY2gsqhb3toI7H0qbTR3PYsz5cFIQpFWb5wRemmiCqWcJ+CmSsr/3tkJOgW73VQiROB9Hb6q2pHVfzc7K1vjgDav5j/ftNuTclgkKKWUI+Ajwc+AE8H0p5TEhxIeEEB+KHnYP0CqEeB34InCfzMKCOMVetduJq1GAMj/p0MMgE7p8/sTx7wZQF+3CdVnxxL5jKnzSkR0pPdpClJJ9uXozjLRDYCzt62rmy6WgUSRVFbWvNZpnot99rTW4TMVlVG9WfpYUKyp0W2DeNRpLZ6eU8jEp5Vop5Wop5V9GX/uKlPIr0X//q5Ryk5Rym5TyKinli1aOdzGEENSUJOGMrN4KM+Oq9LJFmB2NoVFb4iUYlgzOd/hLqWzAWWJ2gvmO2lQ0iqjjM4PMey0xK9s1iqSK3UmpNAqd/U51JV7z6j2Bei5lRPnRUqB7NIDbKagwuM2wmWTHNm4ZoDJHE2kU0QgZi8xPVkVjALNaTOf8HeFop6qVlEURT92jfjXJ81OY5Nr4M/BT5LpUeKOpppU0mC1NEU+QjnbCtP73ta7Ew9DkDIGgSQ5/LfIpRYd2j8+vIv2SycNZItiCQieSCm9csUG3WvfpMOoPWhKNAYsk3ens8NSDHp8KjU1pkhfVqcz7DP0UtcWerDc9aRpP3M3GbICCvrkx2jWNqmd0GSUrIacwZYe22SVyzMAWFDpRk0z1RldutOSDNRqFVdEY6poxHP7aTi1LQmNBs7+n+P0IoXafGXYxrCn2JvZzWUxS9nfte9C5JIvmFzFNmDoc6tlMcQNgdla2GdiCQifqSjzK5J6oXWPN1oxDKdOlx+TSx/Mp9Lgp9LgutTH3tkJpE+QWmj6exej2pZBDMZ+qTcpHEUnfLFJT4qHHl90NjLpHo8XuiuJ8R31HDbmvWoKfaRoFKPNZb2vSASjhiKRvzNYobBZhNrwx0UNctVmVBpjoN2FUl5KUfdlA6kq8dM//fvqOZZXZSZvkSZfvmE/VZghOwXD6gQp1JV4mZ8KMBUJpn8Noun0qvDrHFWfpMOi+VkeFk6laV9VmFYCSZM+RwYlpgmGZ3jOUxdiCQie0xTdx5FN0AmVoz06HLl8ajlodqSn2zJmeZqZg+FxWCYqhiWlCEZle/LsO91XbbGSzQzthH4qZSRg6p0sPioV4c5yU5rnNLXUy69BO7r4mXQdriWELCp2Ym+RJaBSQsT07HXp8KuLJqmiMmpJ5NnitdEdWRTwl4ahdjMoNqodBJoIi2c2GhXT5Etjf+08A0rANQHWx19wQ4hUbQTiSnq9md7YzC1tQ6ER+rosijyvxbjCvDAprLdEoVA6FdTuduhIvw1oDo2x0ZM8m26XxHbk9UNGS0QagNtnNhkVIKelJFNHTa+x9rS32mGt6ysmDstVJz9cekzvbmYUtKHQk6Z4C1ZvnQkNNRBW7s+4B1hbgnlG/+vtzCqCkybLxLKQ70xLsGZbyqCzMxeUQWatRjPqD+IPh+IK075gKKS1ZacgYqq2oslu9OekAlG5fgLwcp2n96M1ief01FqNs8Ensdqo2wblnIDQDrhzjB4ZqX9lrcTTG/OKJq7TM3Swp3QFKo8h1OSjNc6d3gurN0PpD1SfaW5ryx50OQVVRks+QBSQVXt1n7H3VquwGgmE87szrlQWDQTo7OwkE4nzna34X6n1w/JgyQ8Xhpuog17+lipMnT2Y8NqPweDzU19fjdif/nNuCQkdqS7wc7vAlPrBqM0SCMHjKEKdfLPrHpwlHpKl9KBYyV459Su08t/wvy8YSC81RG61DmTpaglnfMWi6Nq1T1JZ4staZnTC8Wkr1t299l2FjmB/51KxDpdTOzk4KCwtpampa/L4HRmH4PJSvhNyCuOc72z+OQwhWVcY/ziqklAwNDdHZ2Ulzc/IFG5MS+0KIXUKI3xdC/J0Q4vNCiHcJIcrSHu0yJemeAtXzFhSTyAbbqVbPaKKvXZV4yKKIJ4h2tsskWqU680CFbE66mw2vXuw78l2A6TFD/U56O/wDgQDl5eXxNweu6JwJJr7mTFjGDx22GCEE5eXl8TWoGMT9i4QQ7xNCHAQ+BXiBU0A/cC3wpBDim0KIxjTHvOy4xAYfj7LV4PKYmninmTOs1Ci0ekbOgehCapI2lSxaVFjaFFRBXkVazW40aku89I4GLq+ymwV0+xIUu5styWLcfTWiHHtCDdLpVhFtofjzOiIloXAEtzN7BQUk8ffGIJHpKR/YK6WM+Q0JIbYDLcDFlK+8DJkfIhtX9XS6oHK9qZFPvZmEfupIbbGH/JFoNc4V+pZ4yIRQOEL/eCCzqDAh1G46k8inEg8z4QhDkzNUFmZX9dFunz9+HazeVkDoXrpjPtpmrDdRBQQ9EQLc3oQaRTCssrdTFRROp5MtW+aE63333ccnP/nJpD//8MMPs3btWjZuVN/7vn37+Pu//3t27dqV0jjiEVdQSCn/LcH7h3UbyTIgpS5l1Zvh9M8NHtEc3aN+FY3hsdYtVVvipbL9rGpok8DeayZ949NEpA6CtHoLvPo1CIfUhiBF5ifdZZugSFgHq+8olK1SXf8MwuNWSXem+3HcXpgcUn6YRXbkwZDSAnOcqe3YvV4vhw8fTmtYoVCIhx9+mDvvvHNWUBhBItPTF+P9GDaqJYrWXjQptbhqC0wOwHifwaNS9I4GqCn2pO+o1YmaYi8rg23ILEq0A+jVupJlapqr2gzhadXwJg2SNl9aQMI6WAb0oIiFJX4ctxeIQGh60UPS1SgW4/Of/zy7d+9m8+bNPPDAA7M1wPbt28ef/MmfcMMNN/A3f/M3PPLII/zhH/4h27dv59y5cwD84Ac/YM+ePaxdu5bnn38+47Ek2vK8Fv29F9gIfC/6/3fOe88mSq7LSUVBTnKTfH7z9sIqYweGyhHIhmzRxkJJI71Ml20km4oczPpwMi29ML+Ux4r1KX9cu0fZFiIbSVQHa3pCNeTa/huGj6Wm2HNpzTCd+NxPj3G8e5EuhTKianm59qtWATEIhiPMhCLk5869v7G2iD9/a3zh6ff72b59++z/P/WpT3HvvffykY98hM985jMAvPvd7+ZnP/sZb33rWwHw+Xz86le/AuDMmTPceeed3HPPPbPnCIVCvPrqqzz22GN87nOf46mnnkr498cjkenpm6Cc2sCNUspg9P9fAX6R0ZWXKbULC98txuyCcgzW3GLsoFA75rUrKg2/TiJaRAcOIRnIX0OD1YOZx1yf4wyFacU6tZD0tcKWexIfv4DSPDe5LkfWaRSDiepgad39TIhkqynx8NrFEcOvcwla/oRcvIpsRKbnKF7M9PTMM8/wt3/7t0xNTTE8PMymTZtmBcW9994b95x33303AFdccQXt7e0pj2khyRpRa4FCYDj6/4LoazYLqCn2cH5gMvGB3lIoqjel5lMwHKF/fDorKlo2TCvV+KJ7VVYJim5fgHw9fDiuHCUs0gx9FkIkv9kwkZ5EwRBaBJ8JJsWaYi++aBi6NyfzpDuNRDt/+k+AMwfKV8d8u21wklA4QktV5uXVA4EAv/d7v8eBAwdoaGjgs5/97CUhrfn58f1AubnKv+V0OgmFMq9GnKwx7a+BQ0KIB4UQDwIHgf8v46svQ2qKvXQn21MgjaYo6dA/Po2UOphVdKBi8gzj0sv5UIXVQ7mEnlE/NZkk280nwxItNcWe2bpT2cKcxrXIM9TXCrnFUGy8+LfMj+PyQmhxAR7UMTRWEwoVFRVMTEzwwx/+cNFjCwsLGR8f1+W6i5HUXyWl/AZwJfDj6M/VmlnK5lJqSzzJ9xSo3gyDp+M6yPQgo2J3OuMdPskp2ZiVO2bdvp+qTTDWBVPDiY+NQTYm3c1pFIsJimPq7zYhWEILGjG1gREoh3Z4BiKx53YwHMGdRrKd5qPQfj75yU9SUlLCBz7wAbZs2cLb3/52du/evejn77vvPv7u7/6OHTt2zDqz9SYVPXsa6AE8wFohxFop5XOGjGoJM9fAyE+xN0EtlarN6qEbOKU63xlERuWz9URKRP8xLrqvybodc7cvwIbqIn1OVjXP/9R8Xcofry3x0DcWIBSO4MqS5K2e0QA5Lgdl+TFqk0Ui6m81wZEN86rsWhL5BAQDl4V2hyMRwhGJO8XQWIBwOHYlhy984Qt84QtfuOz1Z5999pL/7927l+PHj8d8v6KiQhcfRbIlPP438Bzwc+Bz0d+fzfTiQojbhBCnhBBnhRCXZZgIxRej7x8RQuzM9JpGM9vAKJmolRSboqSLbqGfmRIt8TBQsDarNIrpUJjBiWn9vp+qeZFPaVBT7CUilckwW+j2+RcPr/a1w8yEaSVZ5jQKs3Mpos9HjMS7YFjLocgOwa43yf5VHwN2AxeklDcCO4CBTC4shHAC/wbcjgq9vV8IsTBj5HZU5ncL8ADw5UyuaQazCVPJPMRlq5Td02CHtuaoLcy1NtlO+zunStZlVeG7vlG1IOtWB6tgRbSUR5qCYraeUfYI0954pjnt+TVJUHjcTsryc8zfbDgWL+Whdw5FtpHsXxWQUgYAhBC5UsqTwLoMr70HOCulPC+lnAG+C9y14Ji7gP+SipeBEiFETYbXNZQVhbk4HSI5jcLhhBUbTNAoAvo5ajOh7xggkCs20TcWIJwl9Yy69da4hIj2MEjvvtbOM19mC8qHs4gg7YuW316xwbTxVBd5zPdRxCnlMWMLCgA6hRAlwMOoYoA/AbozvHYd0DH/GtHXUj0GACHEA0KIA0KIAwMDGSk7GeFyOqgqzE1Oo4BohEyrKg1gED2ZVkXVi2iJh8ryMoJhyeBEdphWDKmDVbUZBk6qUh4pklKGvwmEtWS7eBFPZatVNziTsKwcu9urfBQL5mswJBGQlo9iKZBs1NM7pJQ+KeVngU8D/wG8PcNrx/pGF66WyRyjXpTyq1LKXVLKXZWV1iaW1ZR4k3+IqzbD1BCM9xo2Hl0jejIhWuJBKwORLean7kShn+lQtVmFUg6nHoVS5HGRn+NMfrNhMAmT7XqPmt77vLrYY25hQI1FSnkEo4EHlmvtBpFQUAghHEKIWR1aSvkrKeUjUXNRJnTCJTlX9VyupSRzTNZRk0pf3/kRMgYwE4owMDFNtdURT1qJh+otl3S6ywZ6fAG1OOvpw5kt0ZK6+UkIQU2JN2s0irjJdoExFaRgcu/z+Ul3pqL1pljgp9AzhyIbSfiXSSkjwOsG9J3YD7QIIZqFEDnAfcAjC455BHhPNPrpKmBUStmj8zh0R/XODiSfdAcZ9TCIR/94ACnjNJsxi3klHmrnVUjNBnpG/frXwaqMlvJI009RY0Vv6EXQQpmrYz1DJvSgiIV1SXexI5+C4UjKVWM1hBC8+93vnv1/KBSisrKSO++8M+1h6k2yW6ga4JgQ4lVgtj6FlPJt6V5YShkSQnwEFWrrBP5TSnlMCPGh6PtfAR4D7gDOAlPAb6d7PTOpKfYwE1I9BRZt8qLhLVHZrAZFPmm7wZiT3Ey0Eg9VmyjyusjLcWZN4btunwGmOVcuVKxNW1OsLfZystfYbNtk0Z6hmMJU05hMNj3N5Ssl6P2iNw6HEhbzBIWUkpmwpMibnkaRn59Pa2srfr8fr9fLk08+SV1dTFfsooRCIVwu46Iakz3z54y4uJTyMZQwmP/aV+b9WwIfNuLaRjK/C1dCQQHK/GRQ5FPcSW4mfcdUiYeSxrl6RlmkUWxvLNH/xFWb4cKv0/poTYmHwYlpZkIRy1tr9oz6yXU5KM2LkUDa1wqeEihKbWHLlBor/VwuLwTn6rmFIhIpZUamp9tvv51HH32Ue+65h+985zvcf//9s+XBX331VT7+8Y/PCpJvfOMbrFu3jgcffJBHH32UQCDA5OQk3/72t7n33nsZGxsjFArx5S9/meuuSz3hMxZxBYUQQkRDU3+V6BhdRrNMqCuZy6XYUl+c+APVW+DML1Q0hVvfnW1cs4GZ9LVeUuIhW0wr/pkwI1NBY0xzVZvg6PdVKY+81FrM1xZ7kRL6xgI0lJkXTRSL7ni9THpblUA02YlrSBmPxz+ZXHvi8IzqOZJTAAgcUrJqJozH7VAaxyUD3QK3/3XCU9533318/vOf58477+TIkSO8//3vnxUU69ev57nnnsPlcvHUU0/xJ3/yJzz00EMAvPTSSxw5coSysjL+4R/+gVtvvZU//dM/JRwOMzU1leo3sCiJNIpnhBAPAT+RUs62O436FK4F3gs8Azyo24iWAbMJU8nudqo3gwzDwAmo3aHrWHpGAxTkuijyJCgnYiRaiYdt98++VFeSHaYV3cqLx6I6/VIe2jPU7fNbLih6F8uhiISV72nne0wfk2VJd3BpyXHhnPVFZhLxtHXrVtrb2/nOd77DHXfcccl7o6OjvPe97+XMmTMIIQgGg7PvvelNb6KsTG1Cdu/ezfvf/36CwSBvf/vbL+lxkSmJBMVtwPuB7wghmgEfqtaTE9WP4p/sdqiXU56fQ47LkXrkU+9RAwSF33ptwndBlXiYZ8euKfYyMD7NdChMrku/UtGpMhvRY0R5kwxqPs23wVtNj8/PVavKL39jpF018zEpI3shNcUefct4JLHzB5RG0XdMtQkoqMQ3MU23z8/GmiLIwPz0tre9jU984hM8++yzDA0Nzb7+6U9/mhtvvJEf//jHtLe3s2/fvtn35pcbv/7663nuued49NFHefe7380f/uEf8p736CPEEzUuCgD/Dvy7EMINVAB+KaVPl6svU4QQqXXhKm1WaqwBDu24pRfMQvO/zIuM0RbmvtFpGsut2zFrNm7dynfMp6AqWsoj9Yi2lPqvG0g4IukbX6QOlok9KGJRU+yhc8SC72dBKY9gOIJDCJyOzMxv73//+ykuLmbLli2XFPYbHR2ddW4/+OCDi37+woUL1NXV8YEPfIDJyUkOHjyom6BIWvxJKYNSyh5bSCRHTXEKmaMOh2G9KbpHA8YsgqnQ2wqIS0o8aH6cLosd2oZGhQkRva+pRz7l5bgo9rotz6UYGJ8mHJGx83D6WpUZpjL1lq96oPV+MZ0FpTxmQiqHItNku/r6ej72sY9d9vof/dEf8alPfYq9e/cuWmkWVNXY7du3s2PHDh566KGY50oXi6vELV9qi728fH4o8YEaVZvh6A9VaQCdHIMzoQiDE9PWm576WlVXsHklHiyLg19Az6if8vwcPG6DzF/VW2D/11UpD2dq0y0bHP7a9WM6+3tbobxlrvy2ydSUeBgLhJicDumbLJkMbi9MDoGUBMPplRfXmJiYuOy1ffv2zZqYrr76ak6fPj373l/8xV8A8L73vY/3ve99s6+/973v5b3vfW/a44jH8k0ltJiaEg990d1YUlRvhulR8F1MfGyS9I1Fk+2sLi/e13qZHTtbsrO7fQFjy69XbYqW8jif8kdVCLG1309cjavvmGVmJ7C4eOK8Uh7LPSsbku9HsbD8N0KIfXoPZjlRU+wlHJH0jyfr0Na/N8XcJLfQ9BQYU07PBQuKFrVivenJb6xpbtahnbqfosaqekbzmM3DWfgd+X0wetEyRzbMbTYsEabRUh4y6CdkC4pZvi+E+ONoKQ2vEOJLwF8ZObClTl2qD3HVRkDo6tCOazYwi/4T6neMBSUbekN3+wLGJiNqpTzS8FPUFHsYnpwhEDS5ntE8enwq2a5kYbLdbOkOKwWFlUl36tqR4BQScLuWZzFAjWQFxZWo4nwvomo0dQN7jRrUcmCu+UySD3FOvrLj9x7RbQxZUb5D20nHWFCsNq2MBYJMTIeMjQrTSnmksQHIhhDZnjElSC9z1FpUumM+VUUehMi8JWpa+cLRUh5yRs3vpdTZLp2/N9m/Lgj4AS8qj6ItWizQZhHml/FIGp1LefSOBijMdVFoZbJd71FV4qG4/rK3aos9loZ/avfG8PImaUY+pZy4aQA9vkV6mfQehbxyKLSuj5jb6aCqMLO+FB6Ph6GhofSEhcuLCAVmx7IUkFIyNDSEx5Pa5ijZUIH9wE9Q7VDLgf8nhLhHSnlPasN845BWT4HqzXD8YWXX9xRlPIZuXxYk2/W2qsifGJFcNSVexgMhxgNBS4TZbA6F0c7+qs1w9AfgHwFvadIfm62ya6VGMRrg6tUxku16j1pSumMhNSWZRYbV19fT2dlJWs3OAmMQ8DEgAzjH8nAskV4UHo+H+vrLN27xSFZQ/I6U8kD0373AXUKId8f7wBsdradASrud6q3qd98xWHl1xmPoHVMtUC0jElZ/y67YRX/nRz5ZIiiMLN8xn/kZ2k3XJv2xuU531mgUoXCE/vHpyx3Z4ZDyPe35gCXjmk9tiZfj3WNpf97tdtPc3Jzeh888CT95F3/LZ/nGZ38/7TEsBZLtcHcgxmv/rf9wlhcpNTCCeQuKPuanbl+AmiILNYqhcyp7dRGHp9Wd7np8ARxC9Tk3lPk1n1LA43ZSblU9I2BgQku2W/AMDZ1RRfG0jY2F1EYTWy2pSxp9rnd5s75FTsYsDcPaEqW2OEVnbVGtMk0kU8EyAVqynaE5AonQHNnVsZvaWBreiNIoqos8uIy2LxdUKXt+Gvc1U9NKJsyVqF/wDPVa78jWqC3xMh2KMDyZacPNNCisZlQUsdmpX+5TtmILCgOpLfEyOKEK3yWFELo5tPvGtPaVFgqK3qMqNLRyXcy3VxTm4hDWZWd3+/zmmOZm72s6IbLWtUTVrltdtOA76j0CzhwVzWUxlkaGCcFJ2cjqcLv51zYZW1AYyPzCd0lTvRX6jiv7fgbE7XNsFr2tqg6QK7Zpx+V0UFXksUyj6DGzYGLVZmXXT/G+WhkZNpuHc5lGcVTdV6eF0XRRrKwZNjUTojXUQNX0+Yzna7ZjCwoDmYtaSTHyKeRX9v0MmOuzYLFGsYjZScOqTndSSnpGDU62m492X1Ms5aFFhk1Mhwwa2OL0jAbwuB0Ue+cJBCmj99V6/wRYG0Lc7QtwQjbijkynVaJlKWELCgNJOekOMir5MB9tl25ZeOzEAEz0Jszctarw3dDkDDOhiHlZ61Wb1O8U/RQ1FkY+9UYrD1+SbDfRB1ODCTcAZqH1frHC4d/l83Mi0qj+Y1Ar42zBFhQGMqtRpGJa0Uo+ZFjKo9vnp8hjYbJdAke2Rl2Jl+7RgOlRK5oWY1r4cOV61cMgRT/FrMPfooXwso1GFjmyQYWh16ZS0l9Hun1+zso6pHAa0ksmm7AFhYF4c5yU5LlT2zG7ctWikmHkU7fPT12phS00ZxeU+IKiptjDTCjCkMlRK5rwNq1Xh1bKI8Wdp5UaRbfPP+sDmEUrMWNhjaeFWGW+7Pb5CTlyVan1NAIVlhK2oDCYmlRDZEGXyKcun586K0Nje49CUR3klcU9TNvRmx3ZM+vDMfM7qk498kmvekapMh0K0z8+TV3pQkFxFEoawVti6njiUVPstSTqqWtEhVeLan1L72QjlggKIUSZEOJJIcSZ6O+YdQ2EEO1CiKNCiMNCiMuS/pYCaanF1ZthvAcmB9O+bles3aCZxOhBEQurola6fX5yXA7K83PMu2jVJhjtUKU8ksTtdFBZkKtvb+gk0AT3Zc9QX+slLW2zgboSD31jAUJhc8vPdfn8KiKsenP0vvpMvb6ZWKVRfBL4pZSyBfhl9P+LcaOUcruUcpc5Q9OX2pI0djuauSZN89NYIMh4IGReRM9CggEYOJWUw9OqTneqRawn4/aVKTHbc+R4Sh+rSecZyhBtc3OJoJiZgqGzWePI1qgp8RKR0DeeQhi6DnSPRjdjVell3i8lrBIUdwHfjP77m8DbLRqH4dSUeBj1B5maSSG8McMmRrOTfKHZwCwGToAMJ7WglOXnkOtymG5jVlVRTf5+tMinFO+rFc7azljPUP8JkJGsExRzGf7mfUfhiKRH62ViCwrDqJJS9gBEf69Y5DgJ/EII8ZoQ4oF4JxRCPCCEOCCEOJBWJUiDSCvyKT9avjnNSIquES1RyiJBkaQjG6JRK9HIJzPpGTW4BWosCqtVKY+UHdpKozAzMqxrxI8QC8KrNUd2lkQ8aVhRM2xgfJpQRKo5VlgN3rKMQ9qzGcM6kgshngKqY7z1pymcZq+UslsIsQJ4UghxUkr5XKwDpZRfBb4KsGvXLgsqhMVmvmllzYqC5D+YgUNbmzD1lgmKo+DOh9LkqnKa3ekuFI7QNxYw34cjhNIqUjQp1pZ4mJoJM+YPUbyw05xBdPv8VBbkkutyzr3YexRyi6BkpSljSJYaC2qGdc3XuIRIK1BhKWGYRiGlvEVKuTnGz0+APiFEDUD0d/8i5+iO/u4HfgzsMWq8RlGbblRP9RYYOAmh1O2unT4/OU4HFQUGV0VdjN4jakF0JPd4md3prm98moi0qLxJzTblowgHk/9IOhn+GdLl88eOeMqCHhQLKch1UeRxmern6lrow6narEvpnWzFKtPTI8B7o/9+L6op0iUIIfKFEIXav4E3A0suBm0uvDGNyKdISAmLFOn2KbOKw2HBhI6EoecI1G5P+iO1xR76xwMETYpaMa1hUSxqtqsS3Vov8WQ+kk6Gf4ZcFjUXDikNt2abaWNIBbNzKeaeoeh3VLNNlWgZPG3aGMzEKkHx18CbhBBngDdF/48QolYI8Vj0mCrgBSHE68CrwKNSyicsGW0G5LjUzj7lh7hmu/rdfTjla3aNTJmXSLaQobMQnJwbfxLMRq2MmaNVdI5MAdBQZkFCYu0O9bvncPIfScfPlQGRqKP2Eo1i8DQEp1LaAJiJ2Vpp14ifYq+bgtyo9T6D+boUsERQSCmHpJQ3Sylbor+Ho693SynviP77vJRyW/Rnk5TyL60Yqx6kFSJb2gy5xSktKBrdCye5mWgTJRWNYl6nOzPoHI4R+mkWpc3Kzp/CglJZmIvLIUzbMQ9OTDMTjlz6/WjPYQobADOpLTG3ym63z39psEhFi/LLpTFflwJ2ZrYJ1JV4ZiORksbhgJqtKe9QZkIR+sYtcNRq9BwGlxcqYvegiIXZUSsdI1NUFubicTsTH6w3DoeqvJrCguJ0qMiwzlSfoTTpjJVD0X1YLYQVLaaMIVVqS7z4poJMmlRl97LKBw6n8ivaGoVNujSU5tE54icSSTEYq3aHsguHkq+D1DcWQEqLdsugJkr1ZnAmH1BndtRK54ifeqs0LlDaVm9rSg7t+lIvHVGTmdHEzMPpOaw2Lg4LhGsS1EfrmpklTGNWPqjdrgI5lqFD2xYUJlBflsdMOMLARIoRTLXbITyjEtiSRJsolpieIhE1UVI0T5gdtaIEhYUFEzWH9sCppD/SUJpHx7BJi+DCPJxIWEU8ZanZCaAh+rx3DBsvTBetfFC7Q/lxlqFD2xYUJlCf7kOsOT67DyX9kcuiMcxk6CzMTKTl8DQraiUckXT7/LMLiyVo308K5qeGMtVWNxA0frfa5fNT6HFRpJWoz3JHNswFJpihdWmh7pfNsWXs0LYFhQk0lKb5EJc2g6c4pQdPi++2pLNdBg5Ps6JW+sYChCLSWo2ibDXkFKZ0X+dMK8YvhJeVF9fGmcUaRXl+Dl630xTTU5dP3YPLtPZl7NC2BYUJzGkUKT7EQqjJmcKD1zXip6LAIkdt92FweVQ/jRQxq9OdtpBY6qPQAhVS1CgAOkxYCDtHFgqKQ1ntyAZVCqa+1GuK6UkzzV3mo5h1aCdvAVgq2ILCBDxuJysKc9PbDdZuV6UBknRod4xM0VhmYcRTVWqObI3aEi8jU0H8M8aaVrR7YKmgALUB6G1ViWxJMKtRGLwQSinpWujs7zmsFsAsdWRrNJTlmSJILw5PketysKIwRuWD2u3Kn7PMHNq2oDAJtdtJ4yGu3aEc2v3Jlaa+ODxFoxWJZJFIyhnZ89GypI2Ohe+0umCiRu32aCZvcg7tyoJcclwOwxdC31SQ8enQXDKi5sjOYv+ERkOpl87hKcOLJ14cnqKhLC92ifpl6tC2BYVJqN1OGrvBWQdZYnU2GI7Q7fNbIyiGz8HMeNp2bK2ekdGd7jqGp6gqssg0Nx+tFEaSfgqHQ5lWjPZRXIxqLLPPkObIzmL/hEZDWR7j0yFG/cmHHadDx3CcObZMHdq2oDCJhtI8ekbT6MJV2gSekqTs2d0+PxGpwnFNJ42M7PnUzeZSGK9RWOrI1ihfAzkFKdmz600IkdU2M43l0e9o9r7uMPS6emBGLoWUko7hqcWj5papQ9sWFCbRUOZVzU5SLVMhhFp8k9ihXLYbNJPugyojOw1HNswVTzS6JWqnb8p6/wQoe3/tDuhKvsNvgwlJd9ozpEXqLQVHtsasw99AP85lprmFOJwqUKHroGFjsAJbUJhEfbohsgC1O1WGdjD+IqrtNi0RFJ371cLnTK9fQo7LQU2Rx9BJHgpH6PEFskNQANTvUvb/BPd19vDSPLVQBYwzrXQMT1FRkEO+Vuyucz/U7cx6RzZkOMeSZFbjijfH6q6AntdTqqiQ7diCwiQaZqNW0tgx1+9WJcd7jsQ97OLwFDlOB1VFJudQhKbVxKjPrK15Y3keFwwUFH3RrmRZYXqCpO+rhrZjNtK0ojlqASXAeo9kfF/NotjrpsjjMtQ8N6txxRMU9btV5v0y6nhnCwqTqCnx4BBp7na0idr5atzDOoaVWcVpdh+K3qMqMqt+d0anWVmWz4WhSZ0GdTmdC80qVlOn3df9SR0+m7hpoDC9JGqu54gSZBneVzNpKMsz1OGftKAA6Ejuvi4FbEFhEm6ng5riNBOCClao9pMJFpRLdoNmoo0rU0FRkcfgxAwTBlUA1bSVrDE9FVZBcWPSgkIbt1EahYqaC8wJCm1cdUtDo4BoTSwDNa6OYT/l+TlzfShiUVwHhbVJ39elgC0oTERVAE3zIW7Yk3CHogSFBYtg534oqoeimoxOs7IsH4CLQ8bsCC8MTeKKhplmDfVXQNdrSR1alp9DXo7TMBt8jy9AOCLnNhud+6GkUQm0JUJDmQohNiqXoiPZzVj9LltQ2KRHQ1le+maD+t0w3g2jXTHfHvUHGfUHrXNk62DHXhkNybw4bIz5qX1QTXKXM4se+/rdMNoB470JDxVC0FiWZ5ggvSxqrvPAkjI7gXJoB4JpVGpOkqS19oY94LsAE/2GjMNssmjGLH+ayvPoH59Or7lKAj9Fh1WhseN94Luoy4Kixe63G7QQtg9NzgqjrEH73jqTC5NtKs+nzSA/ziX297FuGOtccoJCu78XDHiGQrMJrUlopCne12zHFhQm0lxRAKgFK2WqtqiCe4s8eB3JONmMQMsD0GFBKfK4KcvPMWSSSylpH5ykqTxf93NnRPVWcLiTNlM0VeTTMTyVeuJmElwcnsLtFFQXeeb5nfbofh0jaa5Q97dtUH9h2jOqKg8ntRmr2QYO17IxP9mCwkSaKqI75sE0FkJXjioP0BFbo7hglaDoPKAWupqtupyusSzPENPT4MQMkzNhmrJNo3B7VMG9JHeezRV5BMPSkJLsKmouT0XNde4HZ64a2xKirsSLyyEMERRJRTxpuL3R+2oLCpsU0XazaWkUAA27o4k8l9tf2wYmqSjInWs2Yxad+1XrU7c+DuKm8jxDNArtO19ZkWUaBSizYvfBpFqjalrp+cEJ3YfRNjh5qX+iZpvaoCwhXE4HjeV5tBsgKM5Hz7m6siC5D9TvVhnay6CSrC0oTCQ/10VVUS7nB9J8iLVEnt7LE3naBidZZfYiGA6qiB0dzRON5fl0+/zMhPQ1rWgLR3O2mZ4AGq9Shfd6EyfezWml+i6EUkr1DFXmq4zi7kNLJtFuIc3l+YZoFOcHJsjLccYuLx6L+t0QnFRtApY4lggKIcQ7hRDHhBARIcSiT6MQ4jYhxCkhxFkhxCfNHKNRNFfkp69RaAvyxZcve+v84OSsfdY0el5XC9zKa3Q75cqyPCJS/05uF4amcDqENb3EE9EY/f4uvJjw0MqCXPJznLo7/PvGpvEHw6yqjBYqDAV0va9m0lyhBEUkom+IbFt0jsUsLx6LxqvV74sv6ToOK7BKo2gF7gaeW+wAIYQT+DfgdmAjcL8QYqM5wzMO7SFOi6Ia1R71wq8veXksEGRwYprmSpMFhTYOHRcUbcesdymPtqFJGkq9uLMpNFZj9r4mXlCEEDRl8gwtwvkBZcpaVZE/d1+1hW6J0VSRz3QoQu+Yvn6ctlQ3YyUNKqFywXxdilgya6SUJ6SUiTq27AHOSinPSylngO8Cdxk/OmNpKs9neHIm/Zr5TXvVzjMyZ5qZNauYrVFceBHKW1TmuE5ofpy2dM1zi3BhaJKV2Wh20li5Fy5eel8XIyOtdBHOz3+GLvxaVQHOr9D1GmaxyoDIp5lQhI7hqdTNuyuvUfPE4GZKRpOF26tZ6oCOef/vjL4WEyHEA0KIA0KIAwMDA4YPLl20xTxtG/PKayHgu6TjnTYhTPVRRMJqB6yzeaK8IJfSPDdnB/Rz1qrQ2Knsi3iaz8qrwT+SVMe75op8Okf09eO0DU7icTuoLnDBxVeU4FqiNBkgKC4OTxGRpK61N+2FyQEYPKPbWKzAMEEhhHhKCNEa4ydZrSCWIXBRsSyl/KqUcpeUcldlZWV6gzaBjOO8m6ITeJ46e35gEiHmNZsxg/7jMD1qyIKyurKAc/36CYqhSVU/qikbI540NIGbhJmiqTyfcETqWsrj/MAEzRUFOPpbVafCJeqfAKgu8uBxO3QVFG2zGleSEU8aKy+fr0sRwwSFlPIWKeXmGD8/SfIUnUDDvP/XA936j9RcGsvzECIDQVHSqOye7S/MvtQ2OEl9qZdcl4k9AzTHqwELyurKAs7paHrShI7pprlUKG2Gwpqk/BRNmWqlMZiNmmvX3+9kNg6HoKk8X+fvJ/oMpWq+LFsFBVVJBSpkM9lsetoPtAghmoUQOcB9wCMWjyljcl1O6kq8me12ND9F1O6pnGwp7nQypf0FJbBKGhIfmyKrV+QzODHN6JQ+DXrORAXF2qpCXc5nCEIo53ES9my9bfAzoQgdI/6of+JFJbSKanU5t1VkFDQSg7bBScrzcyjOSzFPSQilVVz49ZL2U1gVHvsOIUQncDXwqBDi59HXa4UQjwFIKUPAR4CfAyeA70spl35AMrCqsoBzmdjgV+6FqUEYODUX/27mbllKFfJn0K5TS2jSy09xtn+C/BwnNcUmN3RKlZXXqMKPw+fjHlaan0N5fg5n+vT5fi4OTxGOSFZVeJVDvWnp+ic0mivyuTg8pZsf51x/BuHnK6+BsS5VJHCJYlXU04+llPVSylwpZZWU8tbo691SyjvmHfeYlHKtlHK1lPIvrRirEaxdUcDZ/gnC6cZ5z/opXqB7NMDEdIiWKhM1iv7jykHXdK0hp1+zQv0tGQnTeZzpH2dNVWHy8e9WsepG9fv8swkPbakq4HT/uC6XPRvVuDY6O5VDfQk7sjXWVhUSikhdosOklJzuH6clXY1Umydti2YDZD3ZbHpatqytLmQ6Gm6XFqXNUNwA55/ldK9aLEw1q5x7Rv1efaMhp68vzSPH6dDNoX2mb4KWFSab5tKhfLW6r+eeTnjo2qpCzvZN6NJ34XSfeoZWjb2iXlhlzH01E23jdKo3c2E6MDGNbyrI2nQ3Y5Xrlf9JmzdLEFtQWIC2qJ/qS/MhFkIt0uef43TPiDrnCjMFxdNQsRaK6w05vdMhaK7I10WjGJ0K0j8+vTQEhXZf256HcPxS9C1VhYxPh+gZzTyp7FTfOI1leeS0/woqN2TcgCobWF1ZgEPMCcFMON2rnsN16W7GhIDVN8H5Z5Zs3SdbUFiAtmidyeQhXn0TTI8yfWE/1UWe1J1s6RIMKIenwbvO1SvyZ53QmXB2QH3HpprmMmHVjSrsuPtg3MPWRp8hfRbCcTZV5qj7uvqmjM+XDXjcTpoq8vX5fvq0ZyiDzdjqm5RZr+dwxuOxAltQWEB+rov6Ui+nMnFGNt8AwkFF3wusrTZRm+h4BUJ+wxeUdVVFXByeSq/J0zw0h2+LmRpXJqzaB4iEZgpNK83UoT0dCnN+cJIb886qgpPLRFCA0rJP6+DwP903Tll+DhUFGVTSXbVP/U7CrJiN2ILCItZVFWamUeSVIWt3sn7qwOzu0hTOPa0ashgcGbOhphAp4WSGNuYz/RN43A7qSrKwGGAs8sqgdrsyU8ShND+HioLcjHfM5wcmCUckO4KHwJmzpPMnFrK2upALQ5MEgpmZe071jbO2qiCzYIj8ClW2/awtKGxSoKWqkHMDEwQz6FQ2Wnst2zjLpjL9u50tyrmnVRXbXGN36BtqigA40TOW0XlO9Y7TsqIQhyPLI57ms+pG1aAqEP9vX1tVwOkMzXOaoGkYeUWVO8/J4jInKbKuqpCInIvqSgcpJWf6JvQJFll9k2plnOC+ZiO2oLCIddUFBMOSCxmE750q2I1TSHYED+s3sHiMdqmeCWvfbPil6ku9FHpcGQkKKSXHe8bYGBU6S4Y1N4MMJwyTXRvVSjMpp32qd5wGxzCeoeOw+ua0z5ONaFFKmWhdWvi5PoLiZoiEoP35zM9lMragsIh1VWrxOtad/kL4anA1w7KA+n6Twu5OP6F+r7sj/nE6IIRgQ3VRRoKidyzA8OQMm+qWmKBouAo8JXDqsbiHbawpYmomTFsGm43TfeO8qyjaCGv9W9I+TzbSVJFPjtORkfnyZPT5W6eHH7DhSsgpnJtHSwhbUFhES1UBOS4HrV2jaZ/jSPckr7r34Dr3ZFJtNDPm1GOqdk3FWuOvhfJTnOxNf8d8rEtN8iWnUThdsPZWtaDECZPdXFcMkNEzdKx7jFscB1W5+IqWtM+TjbidDjbUFHK0M/3v52jXKELo9Ay5cqDlTXDq8SUXJmsLCotwOx1srCniSAYPcWvXKJ1VN0Jg1PiiY9PjKrN03R0qLtwENkR3zOlWST3eM4YQsH6pCQpQ37N/BDou72ao0VJVQK7LkfZC2D8eYGJ0mLX+Q7DeeC3RCrbUF9PaNZr2ZqO1a5RVFfnk57r0GdD6t6iqBp379TmfSdiCwkK21hdzrHssrYd4YHyantEArpabweWBk48aMMJ5nHsawjOw7nZjrzMPzaF9PE3z3LHuUZrK8ynQa5KbyZqbVRTSycXNT2rHXMTRNDWK1q5RbnC8jlOGTDEnWsHWuhLGp0Npl/Jo7RpjS1Rz04WWN6v7euKn+p3TBGxBYSGb64qZmA6lZWPWzA3rG6tVNMWpx4ytTnnyMfCWKvu5SayrLsTtFBxJeyEcY2PtEtQmQEWVNd8AJ38W975uqUt/s3G0c4w3OV8jklcB9bszGW3WsqVeLfLpCNOB8Wl6xwKzJj5d8BRF7+ujS6qarC0oLGSr9hCnYTrQHvxNtUVKnR3tgO5Duo5vlqBfPdjr36Ls5ybhcTvZUFPEoYsjKX+2fzxAl8/P9voS/QdmFhvuVBVH42Tzbslgs3Gqo5c3OQ/iWH8HOEzsZWIiLSvSN89pmzFdBQWoeTTSdkmXymzHFhQWsqayAI/bkZaf4kinsp0Wetyw/k6lzh79gQGjBM78QnU923yPMeePw46GEo50jqZcaffQRR8AO1eW6D8os9jwNnC44cji9zVdh7aUkpLOp8kjAFvemdEwsxmX08Gm2qK0tNLXO30IEd2M6cn6O0E4jZuvBmALCgtxzTq0fSl9TkrJaxeG2bmyVL3gLYG1t8HRHyYsJpcWR38I+ZXQdJ3+507AjsZSpmbCKcfCH7w4gtsp2FSr827QTPLKlE279aFFo2RaqtRmQxOMydI54ueGmeeYyqlcFmXF47G1voRjXalvNl67MMK6qkK1GdOTgkpYcwsc+T5ETEyWzQBbUFjMFStLOdI5mlKZgXMDk4xMBdnTVDb34tZ3wWQ/tD2r7wADY0qj2PQOU81OGjsaSwBSXggPXfSxsbYYj3uJm1S2vgsmehftZeB2OtjRUMr+9uGUTnv4TDv7HIfxr33bsjU7aexoLGFyJpxSTk4oHOHghRF2z59jerLtXtXMaIkk39mCwmKubC5nJhzhcIcv6c9oi8KuptK5F1veDJ5itUvRk2M/glAAtrxL3/MmSWNZHmX5Obx2IXk/RTAc4Uinj51RIbOkWXsb5BbFNVPsaS7jRM8YY4Hkc2lCR35IrghRcuVv6jHKrObK5nIAXj4/lPRnTvaOMzkTvnSO6cm6O9R9PfI9Y86vM7agsJjdTWUIAa+cT35HuL99mIqCnEtbM7py1a7/xE9VXoVevPYgrNgE9bv0O2cKCCG4alUZL50bTLpJj9LQIuxaadBu0EzcHtj4Njj28KI1gq5sLiMiSV6YSsnmnh9x0b0KZ/1O/caapVQXe1hZnscrbanNMVBC2BDcXnVfj/8EZvTr7W0UtqCwmOI8N+uri3i1PbndjpSSV9uGuWJl6eXVLK/4bQhOweHv6DO47sMqkuqK95mWZBeLq1dX0D0aoH0oucS7F84MIgRcs7rc4JGZxK73Q3ASXv9uzLd3NJbicgheTXIhHD37Mi2RNs43vcvS+2omVzaXsb99OOkw4lfbhqkr8VJTbGDV4R3vhpkJ/bSKsW6Y6NfnXAuwBUUWcGVzGa9dGEnKT9E2OEnniJ9r11Rc/mbtdlXZdf/X9HGSvfagSubbao3ZSWNvdMH/9dnBpI5/4ewAW+qKKc3PoH9ANlF3BdTuhP1fjxl7781xsq2hhBeT/H58z3+VSZlL8Z7lb3bSuLK5HN9UMKm6T8FwhBfODsaeY3rScCVUb4VXvqpPTsXTX4B/26Oai+mMLSiygBvWVhIIRpKyoT5zagCAfetWxD5gzwdg6Cycz7Du/cQAvP4dFTrpLcnsXBnSXJFPTbEnKUExMR3i0EUfe42e5Gaz5wMweArafhXz7X1rK3m9c5SB8en45xnvo7bjpzzhuJ6tqxsMGGh2cl2Leh6eOZV4x33wwgjjgRA3rq80dlBCwJUfhIETmTu1x7qVf3LLO5W5UmdsQZEFXL26HK/byVMn+hIe++ypftasKKChbJG+ARvvgoJqeOGfMxvUK1+B0DTs/Vhm59EBIQTXtVTwwplBpkPxta4Xzw4SikiuW26CYtPdkL8Cnv/HmG/fvKEKgGdOxl8IIy/9Ow4Z4tSq9+FcSj06MmRFkYdt9cU8eTzxHHvm1AAuhzBns7H5Hsgrhxe/lNl5Xv6yKk1/9Yf1GdcCLBEUQoh3CiGOCSEiQohFvaRCiHYhxFEhxGEhxAEzx2gmHreT69dW8PSJ/rgO28npEK+cH+am9YtoE6Cc2td+XO1Q2l9Ib0DT48p8teGtWVNR9PbNNYxPhxJqFU8c66XY62a3UU5Iq3B7lNBu+xVcvLxQ4IaaQmqKPfzyZJyF0O9D7v86j4evZMvW5e/EXsgtG6o43OGjfzy+aebZU/3sairVP38iFm4PXP0RFYLekWahQP8IHPiGCmYpbdJ1eBpWaRStwN1A7ODwS7lRSrldSmlN2I1J3Lyhiu7RQNyaNL843stMOMIt0d3jolzxPiiogmf/Oj3b54tfUpFT1/5+6p81iL1rKijyuHj0SO+ix8yEIjx5vI83bazC7VyGyvKu34a8CnVfFyCE4OYNK3ju9ODifcZf/BLO4ARfi7yN69cabFbJQm7ZqOZNPK3iTN84J3vHefPGarOGBXseUFrFs/9fep9//h+VU/zaP9B3XPOwZDZJKU9IKU9Zce1s5daN1eS4HDz0Wueix/z4UDf1pV52rUwQ2+32qoem/XlV+z4Vxrrh119Upo667Nl15rgcvGljNU8e713U6f/MqX7GAyHu2GLiJDeTnHwlvM8/A6cub37z9u11+INhHm+NIUxHO5Ev/SuPi2upWX8VxV4TdstZxvrqQlZX5sedYw8f7sLpELx1W615A8stgL0fVxWaz/4ytc/6OuCV/wfb7oPqzYYMD7LfRyGBXwghXhNCPBDvQCHEA0KIA0KIAwMDAyYNTz+K89zctqmahw93x1wIO4aneOHMAG/fXpdc/+fdvwOV6+GJP4aZFPo5PPFJZeu8+TMpjN4c7t5Zx1ggxM+O9MR8/9uvXKSm2MP1Lct4t7znAdVk6IlPqmKN87hiZSlN5Xl8f3/HpZ+REh7/YyKRCF/wv5P/dUW9iQPOHoQQ3Lu7gYMXfZyJURJmJhThh691cu2aCioLc80d3JUfhLLV8Ngnko9akhIe/b8qs/7GPzV0eIYJCiHEU0KI1hg/d6Vwmr1Syp3A7cCHhRDXL3aglPKrUspdUspdlZVLc6H4jSsbGfUH+d7CiQ78xwttOITgN69qTO5kTje85R/AdxF+8WfJfab1IZUAtO9TUNacwsjN4ZrV5bSsKOAbv267zJdzpm+c504PcN/uRlzL0eyk4cqBt/y9qj76i09f8pYQgt+8ciWvtg9fmnx39Adw8md8y/tbUNLIvnVLc37owd0768lxOfjqc+cve++nr3fTNzbNb+9tMn9grlw1X4fPwy8/l9xnXv8unPk53PRpKDE2gs2wGSWlvEVKuTnGz09SOEd39Hc/8GNgj1HjzQaubC5jT1MZX3723CV25otDU3x3/0Xu2l6XWgJQ07Vwzf+BA/+RuLRHbys88lGo2wXXfDTNv8BYhBA8cP0qjnWP8dMFWsXfPHGKwlwX7756pUWjM5FV+5QDdP/XVMHGefzGlY2U5rn5xydPKWHacwR++nFGK3byuaGb+OANq5an/yZJKgpy+Y09jfzoUBfnBiZmXw8Ew3zx6TOsqyrkBqv8N6tvhD0fhJf/XWXix6PndfjZ70PjNUobMZisfWKEEPlCiELt38CbUU7wZYsQgj++fR194wH+/JFjSCkJhiP88UNHcDkcfOLWNHpV3/QZVR304d+FEz+LfUz/Cfj2O1WznHv/25Lif8ly9856NtcV8fmfHqcz2iL14UNdPHWij9+9cTVlyyXJLhE3f0YtEj/+0CVd8PJzXfzBm9by67ND/PTJp+B/3kXYU8xvjf0eq1cUce/uN07uxGL83o2rKch18fvfOzxr5v3bJ05xYWiKP3/rxssrHpjJm/9CJc3+6ANw+uexj+lthW/do6oLv+ubphR1FMnWz9H1okK8A/gSUAn4gMNSyluFELXA16WUdwghVqG0CAAX8D9Syr9M5vy7du2SBw4s3Wjaf/zFKb749Fn2rilnPBDiSOcof//ObdyTrm05MAb/dRd0H1QaxrV/oB6y0LRKqvvFZ1SY3m/9yFCHmF6c6Rvn7i+/SF6Ok91NZTze2ssVjaV8+wNXvrF2y36fuq89r8PejyqHaF4ZkRk///2Vv+IdQ18l4srjo65P8+pUNd974Gq2NZRYPOjs4InWXj70rdfYUFNEVVEuz54a4H3XNPHZt22yemgq3PW/7lLa4N6PqXB3b6nyXRz+Fjz1OcgpgPc8DJXrdLusEOK1xaJLLREURrPUBYWUkq8/38aDL7bjcTv46M0t3LW9LrOTBv3w+B/DwW+CcEBxg6oLE/JD49Xwjv8HpUvHbNPaNcrnfnqMtsFJbtlQxZ/duXFp9sbOlJkpePwP4dC3VDOc4vrZ+9pesJ2PBD6Eu6yRP3vLBq5YDkUSdeSJ1l7+8clTTARC3LenkQ/fuCZ7khCnJ9R8PTz/vvapSs4rr4V3fEV3v4QtKGzm6DsGx34MIxdU7PbaN0PzPnC8gXbiy5Heo8qu7bsI+RXQ8iZYdeMbpujfsqX3qJqvvg51X9feqnpuG3Bf4wmKN+AW7A1O1Sb1Y7O8qN6ifmyWF1lyX+1tpI2NjY1NXGxBYWNjY2MTF1tQ2NjY2NjExRYUNjY2NjZxsQWFjY2NjU1cbEFhY2NjYxMXW1DY2NjY2MTFFhQ2NjY2NnFZlpnZQogB4EKaH68A4vfbXH7Yf/Py543294L9N6fKSillzNK5y1JQZIIQ4sByb7u6EPtvXv680f5esP9mPbFNTzY2NjY2cbEFhY2NjY1NXGxBcTlftXoAFmD/zcufN9rfC/bfrBu2j8LGxsbGJi62RmFjY2NjExdbUNjY2NjYxMUWFFGEELcJIU4JIc4KIT5p9XiMRgjRIIR4RghxQghxTAjxMavHZBZCCKcQ4pAQ4mdWj8UMhBAlQogfCiFORu/31VaPyWiEEL8ffa5bhRDfEUJ4rB6T3ggh/lMI0S+EaJ33WpkQ4kkhxJno71I9rmULCtTCAfwbcDuwEbhfCLHR2lEZTgj4v1LKDcBVwIffAH+zxseAE1YPwkT+BXhCSrke2MYy/9uFEHXAR4FdUsrNgBO4z9pRGcKDwG0LXvsk8EspZQvwy+j/M8YWFIo9wFkp5Xkp5QzwXeAui8dkKFLKHinlwei/x1GLR521ozIeIUQ98Bbg61aPxQyEEEXA9cB/AEgpZ6SUPksHZQ4uwCuEcAF5QLfF49EdKeVzwPCCl+8Cvhn99zeBt+txLVtQKOqAjnn/7+QNsGhqCCGagB3AKxYPxQz+GfgjIGLxOMxiFTAAfCNqbvu6ECLf6kEZiZSyC/h74CLQA4xKKX9h7ahMo0pK2QNqMwis0OOktqBQiBivvSHihoUQBcBDwMellGNWj8dIhBB3Av1SytesHouJuICdwJellDuASXQyR2QrUbv8XUAzUAvkCyF+y9pRLW1sQaHoBBrm/b+eZaiqLkQI4UYJiW9LKX9k9XhMYC/wNiFEO8q8eJMQ4lvWDslwOoFOKaWmLf4QJTiWM7cAbVLKASllEPgRcI3FYzKLPiFEDUD0d78eJ7UFhWI/0CKEaBZC5KAcX49YPCZDEUIIlN36hJTyH60ejxlIKT8lpayXUjah7vHTUsplvdOUUvYCHUKIddGXbgaOWzgkM7gIXCWEyIs+5zezzB3483gEeG/03+8FfqLHSV16nGSpI6UMCSE+AvwcFSHxn1LKYxYPy2j2Au8GjgohDkdf+xMp5WPWDcnGIP4P8O3oJug88NsWj8dQpJSvCCF+CBxERfcdYhmW8xBCfAfYB1QIITqBPwf+Gvi+EOJ3UALznbpcyy7hYWNjY2MTD9v0ZGNjY2MTF1tQ2NjY2NjExRYUNjY2NjZxsQWFjY2NjU1cbEFhY2NjYxMXW1DY2MRBCFEuhDgc/ekVQnRF/z0hhPh3g675cSHEe+K8f6cQ4nNGXNvGJhZ2eKyNTZIIIT4LTEgp/97Aa7hQ8f87pZShRY4R0WP2SimnjBqLjY2GrVHY2KSBEGKf1s9CCPFZIcQ3hRC/EEK0CyHuFkL8rRDiqBDiiWipFIQQVwghfiWEeE0I8XOt1MICbgIOakJCCPFRIcRxIcQRIcR3AaTa3T0L3GnKH2vzhscWFDY2+rAaVb78LuBbwDNSyi2AH3hLVFh8CbhHSnkF8J/AX8Y4z15gftHCTwI7pJRbgQ/Ne/0AcJ3uf4WNTQzsEh42NvrwuJQyKIQ4iioD80T09aNAE7AO2Aw8qSxHOFElsBdSw6V1iY6gym88DDw87/V+VGVUGxvDsQWFjY0+TANIKSNCiKCcc/5FUPNMAMeklInakPqB+W0734JqPPQ24NNCiE1Rs5QneqyNjeHYpicbG3M4BVRq/aqFEG4hxKYYx50A1kSPcQANUspnUM2WSoCC6HFrgdYYn7ex0R1bUNjYmEC0xe49wN8IIV4HDhO7R8LjKA0ClHnqW1Fz1iHgn+a1Mb0ReNTIMdvYaNjhsTY2WYYQ4sfAH0kpzyzyfhXwP1LKm80dmc0bFVtQ2NhkGdEmQ1VSyucWeX83EJRSHjZ1YDZvWGxBYWNjY2MTF9tHYWNjY2MTF1tQ2NjY2NjExRYUNjY2NjZxsQWFjY2NjU1cbEFhY2NjYxOX/x+FDzm5vXsjMQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.pylab import plt\n",
    "plt.rcParams[\"figure.figsize\"] = (12,6)\n",
    "\n",
    "plt.plot(t_grid, e_hist, label=\"Earth\")\n",
    "plt.plot(t_grid, m_hist, label=\"Mars\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"x (rad)\")\n",
    "plt.legend();"
   ]
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}