{ "cells": [ { "cell_type": "markdown", "id": "eac3271e", "metadata": {}, "source": [ "# The wavy ramp\n", "\n", "In this example, we will integrate the dynamics of a point mass bouncing down a non-linear ramp. The point mass behaves like a rubber ball, meaning that at every bounce it will dissipate part of its kinetic energy due to the [inelastic](https://en.wikipedia.org/wiki/Inelastic_collision) character of the collisions.\n", "\n", "Let us begin, as usual, with the definition of the dynamics:" ] }, { "cell_type": "code", "execution_count": 1, "id": "553b9d9f", "metadata": {}, "outputs": [], "source": [ "import heyoka as hy\n", "\n", "x, y, vx, vy = hy.make_vars(\"x\", \"y\", \"vx\", \"vy\")\n", "\n", "eqns = [(x, vx),\n", " (y, vy),\n", " (vx, hy.expression(0.)),\n", " # Downwards constant acceleration.\n", " (vy, hy.expression(-9.8))]" ] }, { "cell_type": "markdown", "id": "cbe82079", "metadata": {}, "source": [ "In other words, we are setting up a 2D dynamical system in which the only acceleration is due to the (constant) gravitational field.\n", "\n", "The contour of the ramp is defined by the non-linear equation\n", "\n", "$$\n", "y = 1 - x + 0.05 \\cos\\left(11 \\pi x \\right).\n", "$$\n", "\n", "Let us take a look:" ] }, { "cell_type": "code", "execution_count": 2, "id": "7faf59f9", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiTUlEQVR4nO3deXRUZZ7/8fc3K0RkCwmEzRASlICJQpAltmxDFpDNhUUBO0QDIorjzGk36N/hqMPY2tPIiIjYLDbN0t3gCK1C2yraA6IEsd0RRFQaRlFAh7AmeX5/JDOThkAKqMqt5fM6h3Oouk/u/Twn4cPNrarnmnMOEREJfVFeBxAREf9QoYuIhAkVuohImFChi4iECRW6iEiYiPHqwC1atHCpqaleHV5EJCRt3br1O+dcUm3bPCv01NRUSktLvTq8iEhIMrMvz7RNl1xERMKECl1EJEyo0EVEwoQKXUQkTKjQRUTChApdRCRMqNBFRMKECl1EJEyo0EVEwoQKXUQkTKjQRUTChApdRCRMqNBFRMKECl1EJEyo0EVEwoQKXUQkTKjQRUTCRMgV+t69e8nIyODPf/6z11FERIJKyBV648aN6dq1K8OHD2fGjBlexxERCRqe3VP0fJkZOTk5tG3blqeeeoqvvvqKRYsWERUVcv83iYj4Vci2YIsWLSgqKmLdunWUlJR4HUdExHMhW+gAF198MePGjWP16tU88MADXscREfFUSBc6VF1THzduHHPmzGHFihVexxER8UzIFzpAYmIiI0eOpKSkhO3bt3sdR0TEE3UWupktNLNvzezDM2w3M5tjZjvN7H0z6+b/mHVLT0+nR48eDB06lPLyci8iiIh4ypcz9MVAwVm2FwIZ1X9KgHkXHuv8XH311Zw8eZK77rrLqwgiIp6ps9Cdc28CB84yZDjwnKuyGWhqZin+CnguoqKiGD58OEuWLOH111/3IoKIiGf8cQ29DfB1jcd7qp87jZmVmFmpmZXu37/fD4c+XbNmzfiHf/gHJkyYwPHjxwNyDBGRYOSPQrdannO1DXTOPeOcy3HO5SQlJfnh0LW78soriY+P55577gnYMUREgo0/Cn0P0K7G47bAXj/s97yZGYMHD2bx4sX89a9/9TKKiEi98UehrwEmVL/bpRfwg3Nunx/2e0ESExPp1asXEyZMoLKy0us4IiIB58vbFpcDbwGXmtkeMys2s8lmNrl6yEvALmAnsACYErC056hPnz7s27ePJ554wusoIiIBV+fiXM65sXVsd8AdfkvkRzExMQwZMoSZM2cyYcIEEhMTvY4kIhIwYfFJ0bNJTU0lLS1NC3iJSNgL+0IHGDhwIOvXr+fNN9/0OoqISMBERKE3atSI/v37U1xcTEVFhddxREQCIiIKHaBbt24cO3aMRx991OsoIiIBETGFHhUVxeDBg5k1axb/9V//5XUcERG/i5hCB2jdujWXX345xcXFXkcREfG7iCp0gL59+/KXv/yF9evXex1FRMSvIq7QGzRoQH5+PrfddhsnT570Oo6IiN9EXKEDdOnShbi4OGbMmOF1FBERv4nIQjczCgsL+fd//3c+//xzr+OIiPhFRBY6QIsWLejRowcTJkzwOoqIiF9EbKFD1S3rPvnkExYuXOh1FBGRCxbRhR4bG8vQoUP5x3/8R7766iuv44iIXJCILnSAtLQ0srKyGDFihNZNF5GQFvGFDtCvXz+++eYb/vmf/9nrKCIi502FTtW66TfccAMLFixg2bJlXscRETkvKvRqzZs35/rrr2fSpEls3rw5oMeqqKjgiy++4OOPP+bAgQMBPZaIRA4Veg0dOnRg4MCBFBYWsm3bNr/u+4033qCoqIiOHTty0UUX0aVLF3r16kXLli1p1qwZBQUFrF271q/HFJHIUuct6CLNlVdeSWVlJQMGDODFF1+kT58+572vsrIyZs+ezbPPPsvBgwfp3Lkzubm5pKSkkJCQAEBlZSXff/89n376KePHj6d9+/YsWrSI7t27+2tKIhIhVOi16N69O9HR0eTl5fGLX/yCKVPO7b7XO3fuZObMmbzwwgskJyfTs2dPOnXqRHR09Gljo6KiSEpKIikpid69e/P2229zzTXXcP/99zN9+nR/TUlEIoAK/QyuuOIKEhMTmT59OsuWLWPu3LlkZ2efcXxFRQXPP/88s2fPZuvWrWRlZTF+/HiSk5N9PmZMTAy5ubmkp6fzb//2b2zfvp0lS5YQFaUrYyJSN3POeXLgnJwcV1paes5f9z+XMcrLywOQ6nQnTpxg48aNbNmyhUsvvZT8/Hy6d+9OUlIShw4d4v3332fTpk1s2rSJhg0bcsUVV3DFFVfQsGHDCzru4cOHWbZsGX369GHlypUqdREBwMy2Oudyat2mQvfN8ePH+fTTT9m9ezcHDhzg6NGjxMfH07x5c1q1akV6ejpJSUl+PeaRI0dYsmQJQ4cOZcGCBX7dt4iEprMVui65+Cg+Pp7s7OyzXnbxt4SEBG6++WaeffZZunbtyrRp0+rt2CISevR7fJBr3Lgxo0eP5oEHHmDDhg1exxGRIKZCDwFt2rRh0KBBjBkzhh9++MHrOCISpFToIeKKK64gOTmZsWPHeh1FRIKUCj1EmBmDBw/mrbfeYv78+V7HEZEgpEIPIQ0bNmTEiBHce++9fPPNN17HEZEg41Ohm1mBmW03s51mdl8t25uY2Voz+6uZfWRmRf6PKgCpqalkZGRwyy23eB1FRIJMnYVuZtHAXKAQyATGmlnmKcPuAD52zmUD/YBfmlmcn7NKtYEDB7Jp0yb+8Ic/eB1FRIKIL2foVwE7nXO7nHMngBXA8FPGOOBiMzOgEXAAqN9P/kSQhg0bUlBQwNSpUykrK/M6jogECV8KvQ3wdY3He6qfq+lJoDOwF/gAmOacO+1+bmZWYmalZla6f//+84wsAJmZmTRr1oy77rrL6ygiEiR8KXSr5blT1wvIB94DWgNXAE+aWePTvsi5Z5xzOc65HH9/TD7SmBkFBQWsWLHC72u3i0ho8qXQ9wDtajxuS9WZeE1FwGpXZSfwBXCZfyLKmTRt2pSrr76an/70p3i1Jo+IBA9fCn0LkGFmHapf6BwDrDllzFfAQAAzawlcCuzyZ1CpXc+ePfnmm2948sknvY4iIh6rs9Cdc+XAVGA98AnwO+fcR2Y22cwmVw97COhjZh8ArwL3Oue+C1Ro+T/R0dEMGTKEGTNmcPDgQa/jiIiHfFpt0Tn3EvDSKc89XePve4E8/0YTX7Vv356MjAxuu+02vZVRJILpk6JhYsCAAaxfv5433njD6ygi4hEVephISEhg0KBBTJw4kYqKCq/jiIgHVOhhJCsri8rKSh5++GGvo4iIB1ToYeR/VmR87LHH+Nvf/uZ1HBGpZyr0MJOcnExOTo4W7xKJQCr0MJSbm8vWrVtZtWqV11FEpB6p0MNQXFwcQ4YMYcqUKRw9etTrOCJST1ToYSojI4OkpCQt3iUSQVToYSwvL49ly5bxzjvveB1FROqBCj2MNWnShL59+zJq1ChdehGJACr0MHfVVVeRkJDATTfd5HUUEQkwFXqYMzOGDRvGhg0bmDNnjtdxRCSAVOgRoGHDhowaNYoHHniAl19+2es4IhIgPq22KKGvdevWDBs2jNGjR/Pqq6/So0cPv+5/x44d/Pa3v2Xjxo18+eWXHD16lNjYWFq3bk2/fv0oLi6mQ4cOfj2miPw9naFHkEsvvZT+/fszaNAgNm7ceMH7O3r0KI8//jidO3cmKyuLVatWERsbS+/evSksLKRfv34kJiayZs0aunTpQmFhIV988YUfZiIitdEZeoTp1q0bMTExFBQU8OSTT57XEgG7du1i5syZPP/88yQmJtK9e3duuOEGYmJO/3G65JJL6NatG4cPH2bjxo1kZWXx0EMPcffdd/thNiJSkwo9AmVlZdG4cWPuvPNOXnjhBebPn48vN+1+5ZVXmDVrFm+99RZZWVmMHz+e5ORkn47ZqFEj8vPzyczMZObMmWzbto1FixYRFaVfEkX8Rf+aIlRqaiqTJk3iyy+/JC0tjZtuuonXXnuNkydP/u+YY8eO8eKLL1JSUkKbNm0YNWoU0dHR3HXXXQwePNjnMq+pXbt2FBcX8+qrrzJs2DAqKyv9OS2RiGZe3S0+JyfHlZaWnvPXlZWVMXv2bMrLywOQKjLt37+fbdu2sWPHDg4dOkRCQgKVlZUcOXKE5ORk2rdvT2ZmJu3atfPbGfXRo0f5zW9+Q+/evXXbPJFzYGZbnXM5tW5ToUtNx48f59ixYwBcdNFFtV4X95cjR46wcOFCbrvtNv7lX/4lYMcRCSdnK3RdcpG/Ex8fT5MmTWjSpElAyxyqbps3ZswY5syZw8qVKwN6LJFIoEIXT7Vo0YKRI0dSUlKitzSKXCAVunguPT2drKwsRo4cqRdJRS6ACl2CQr9+/di/fz8/+9nPvI4iErJU6BIUYmJiuO6665g3bx7vvvuu13FEQpIKXYJGixYtyM3NZdy4cbr0InIeVOgSVHr16sWhQ4eYOXOm11FEQo4KXYJKdHQ0w4YN45e//CW7du3yOo5ISPGp0M2swMy2m9lOM7vvDGP6mdl7ZvaRmb3h35gSSVJSUsjOzmbixIleRxEJKXUWuplFA3OBQiATGGtmmaeMaQo8BQxzznUBbvR/VIkkffv2Zdu2baxatcrrKCIhw5cz9KuAnc65Xc65E8AKYPgpY24CVjvnvgJwzn3r35gSaeLj4xk0aBDTpk3jxIkTXscRCQm+FHob4Osaj/dUP1dTJ6CZmW0ws61mNsFfASVyde3albi4OO6//36vo4iEBF8K3Wp57tQVvWKA7sAQIB+YYWadTtuRWYmZlZpZ6f79+885rEQWM6OgoID58+drWQARH/hS6HuAdjUetwX21jJmnXOuzDn3HfAmkH3qjpxzzzjncpxzOb7cUEGkZcuWZGdnU1xc7HUUkaDnS6FvATLMrIOZxQFjgDWnjHkB+ImZxZhZAtAT+MS/USVSXXPNNZSWlrJmzak/diJSU52F7pwrB6YC66kq6d855z4ys8lmNrl6zCfAOuB94B3gWefch4GLLZGkQYMGDBo0iDvuuEMvkIqchU/vQ3fOveSc6+Sc6+ice6T6uaedc0/XGPOYcy7TOdfVOTc7QHklQl1++eXEx8dz3321fgxCRNAnRSVEmBmFhYXMnz+fHTt2eB1HJCip0CVktGjRgh49ejBhgt4VK1IbFbqElKuvvprPPvuMJUuWeB1FJOio0CWkxMbGUlhYyD/90z/x448/eh1HJKio0CXkZGRk0KZNG116ETmFCl1CUkFBAa+//jorV670OopI0FChS0hKSEhg6NCh3H777Xz7rdaCEwEVuoSwjIwM0tPTGTFihG5ZJ4IKXUJcXl4eX375Jffcc4/XUUQ8p0KXkBYbG8uoUaP49a9/zfLly72OI+KpGK8DiFyopk2bct1113HbbbeRnJzMwIED/X6MsrIyXn75Zd5991327t1LVFQUbdu2JTc3l0GDBhEVpXMj8Z5+CiUspKWlMWTIEEaOHMnGjRv9ss/y8nKWLFlCz549SUxMZNq0aaxfv55du3axc+dO1q5dy80330xSUhJTpkzhhx9+8MtxRc6XztAlbGRmZlJeXk5BQQELFy7kxhvP79a233zzDbNmzWLFihVUVlaSk5PDtGnTSEhIOG2sc459+/bxxhtvkJqayhNPPKH3x4tnVOgSVrKysmjUqBFFRUX86U9/Yu7cucTFxfn0tf/5n//JrFmzeP311+nQoQP5+fl06NABs9pu2lXFzGjdujWjRo1i586dTJ06lddee42FCxfqMozUO/3ESdhJS0ujpKSEDRs20L59ex599FEOHz5c69jdu3fz85//nI4dOzJ48GAOHz7M5MmTGTVqFGlpaWct81Olp6dTUlLCK6+8QkFBAeXl5f6akohPzLlTbw9aP3Jyclxpaek5f11ZWRmzZ8/WPxapk3OOzz//nLfeeos9e/aQlpZG27ZtiYuL48CBA+zevZuDBw/SsWNHLr/8cjIyMoiOjr7g4x4/fpzly5fTuXNn/vjHP+pMXfzKzLY653Jq26ZLLhK2zIz09HTS09M5cuQIf/vb3zh06BDOOVq1akXXrl1p1aqVX0q8pvj4eMaMGcPixYuZOHEiixcv9uv+Rc5EhS4RISEhgYyMjHo7XoMGDRg3bhwLFizgiSeeYNq0afV2bIlc+l1QJEAaNWrEjTfeyIMPPsjmzZu9jiMRQIUuEkBt27alf//+XH/99Wd8YVbEX1ToIgHWvXt3mjRpQlFRkddRJMyp0EUCzMwYMmQI69at4/e//73XcSSMqdBF6kGjRo249tpruf3227VEgASMCl2knlx22WWkpKQwadIkr6NImFKhi9SjvLw81q5dy4YNG7yOImFIhS5Sjxo3bky/fv0oLi6moqLC6zgSZlToIvUsJyeHkydPMmPGDK+jSJhRoYvUs6ioKIYMGcKcOXP48ssvvY4jYUSFLuKBVq1akZ2drfemi1/5VOhmVmBm281sp5ndd5ZxPcyswsxu8F9EkfDUt29ftm7dyvPPP+91FAkTdRa6mUUDc4FCIBMYa2aZZxj3KLDe3yFFwlF8fDz5+flMnTqVEydOeB1HwoAvZ+hXATudc7uccyeAFcDwWsbdCawCvvVjPpGwlpmZSUJCAvfee6/XUSQM+FLobYCvazzeU/3c/zKzNsBI4Omz7cjMSsys1MxK9+/ff65ZRcKOmVFYWMgzzzzDjh07vI4jIc6XQq/tHlyn3uZoNnCvc+6sb6x1zj3jnMtxzuUkJSX5GFEkvCUmJtKzZ09uueUWr6NIiPOl0PcA7Wo8bgvsPWVMDrDCzHYDNwBPmdkIfwQUiQR9+vThs88+Y+nSpV5HkRDmS6FvATLMrIOZxQFjgDU1BzjnOjjnUp1zqcAfgCnOuf/wd1iRcBUbG0thYSF33303ZWVlXseREFVnoTvnyoGpVL175RPgd865j8xssplNDnRAkUiRnp5OSkoKd911l9dRJET59D5059xLzrlOzrmOzrlHqp972jl32ougzrmfOuf+4O+gIpEgLy+PFStW8O6773odRUKQPikqEkQaN27MNddcw4QJE6isrPQ6joQYFbpIkOnRowfff/89jz76qNdRJMSo0EWCTHR0NMOGDeORRx7hgw8+8DqOhBAVukgQat26Nbm5uYwcOVLLAojPVOgiQapXr16YmVZkFJ+p0EWCVFRUFCNGjOCll17iV7/6lddxJATEeB1ARM6sUaNGjBkzhhkzZpCens7QoUMDdqzvvvuODz/8kGPHjpGcnEzXrl2Ji4sL2PHE/1ToIkGuVatWDB8+nLFjx7Jq1Sry8/P9st/y8nJWrlzJ8uXLefvttzl06BDNmjUjJiaGI0eOcPToUbp27UpRURG333470dHRfjmuBI45d+o6W/UjJyfHlZaWnvPXlZWVMXv2bMrLywOQSiR4ffrpp/zxj3/kueeeY+TIkee9n3379jFr1iyWL19OTEwMnTt3plOnTrRs2ZKoqP+7CnvkyBG2b9/OO++8Q3R0NPPmzePaa6/1x1TkApjZVudcTm3bdIYuEiIuu+wyoqOjGT9+PHfeeSePPPLI3xVwXTZt2sRDDz3Ehg0byMjIYNiwYbRr1w6z2hZUhYSEBK688kqys7P58MMPGTNmDDfffDPz5s07p+NK/dF3RSSEZGRkUFRUxG9+8xsyMzNZt27dWccfPnyYOXPm0KVLF/Lz8zl27BhTpkzh+uuvp3379mcs85qioqLIysqipKSEtWvXkpeXp7dSBildchEJQZWVlWzdupWNGzfSpEkT+vXrR3Z2NklJSfzwww98/PHHbNmyhQ8++IDWrVuTnZ1Nly5diIm5sF/Kjx8/zsqVK2nXrh1//vOfL3h/cu7OdslFhS4SwioqKti9eze7d+/mwIEDHDt2jLi4OJo1a0arVq1IS0ujUaNGfj3miRMnWLp0KdnZ2bzwwgs+neWL/+gaukiYio6OpmPHjnTs2LHejhkXF8fYsWNZtGgR9913n9acCSK6hi4i56xhw4aMHj2auXPnsnr1aq/jSDUVuoicl8TERIYPH87EiRPZs2eP13EEFbqIXIBOnTqRmZnJDTfcoPXbg4AKXUQuyIABA/jiiy+YNWuW11EingpdRC5IbGwsI0eOZNasWXz22Wdex4loKnQRuWApKSl0796dcePG4dVboUWFLiJ+8pOf/IRdu3bx9NOn3Tte6okKXUT8IjY2liFDhvDAAw/w/fffex0nIqnQRcRv0tLSSE1N5dZbb/U6SkRSoYuIXw0aNIhXXnmF119/3esoEUeFLiJ+ddFFF9G/f39uvfVWKioqvI4TUVToIuJ33bp14+TJkzz88MNeR4koKnQR8buoqCgGDx7M448/zt69e72OEzFU6CISECkpKVx++eUUFRV5HSViqNBFJGD69u3L5s2befHFF72OEhF8KnQzKzCz7Wa208zuq2X7zWb2fvWfTWaW7f+oIhJqGjRoQH5+PpMmTdJt6+pBnYVuZtHAXKAQyATGmlnmKcO+APo657KAh4Bn/B1UREJTZmYmDRs2ZPr06V5HCXu+nKFfBex0zu1yzp0AVgDDaw5wzm1yzh2sfrgZaOvfmCISqsyMwsJC5s6dy+eff+51nLDmS6G3Ab6u8XhP9XNnUgy8XNsGMysxs1IzK92/f7/vKUUkpCUmJtKzZ08mTJjgdZSw5kuh13YH2FqXUzOz/lQV+r21bXfOPeOcy3HO5SQlJfmeUkRCXp8+fdi+fTtLly71OkrY8qXQ9wDtajxuC5z2xlIzywKeBYY757Qyj4j8ndjYWAoLC7n77rs5fPiw13HCki+FvgXIMLMOZhYHjAHW1BxgZu2B1cB455xWuBeRWqWnp5OSksKUKVO8jhKW6ix051w5MBVYD3wC/M4595GZTTazydXDfg4kAk+Z2XtmVhqwxCIS0vLy8li9ejXr1q3zOkrYifFlkHPuJeClU557usbfbwW0XqaI1Klx48YMHjyY8ePHs337dpo3b+51pLChT4qKSL3r0qULl1xyCdddd53XUcKKCl1EPJGfn88nn3zCz372M6+jhA0Vuoh4Ii4ujjFjxjBv3jyee+45r+OEBRW6iHimefPm3HjjjUyZMoU1a9bU/QXnoLKykjfffJMpU6bQrVs3WrZsSdOmTUlMTKRz584UFRWxceNGvx7Taz69KCoiEiiXXHIJI0aM4KabbmLp0qWMGDHigva3f/9+fvGLX7B06VKOHj1Kp06d6NixI7169aJBgwaUl5fz/fff8+mnn1JQUECnTp146qmn6Nmzp38m5CEVuoh4LiMjg5EjRzJu3Djuv/9+HnzwwXPexzvvvMNDDz3Eq6++SocOHcjLy6NDhw6Ynf5h98TERDp16sSAAQPYsmULAwYMYPLkyTz22GNERYXuhYvQTS4iYSU9PZ1bbrmFX/3qV/Tp04cdO3bU+TXl5eUsWrSI7OxsBgwYwI8//sjkyZMZNWoUaWlptZZ5TXFxceTm5lJcXMyyZcvIy8vj2LFj/ppSvVOhi0jQSE5OZtKkScTGxpKVlcXAgQP59a9/zb59+3DOUVlZya5du1iwYAHDhw+nefPmTJ8+ndTUVO6++24GDBhAkyZNzvm4iYmJFBUV8dVXX9GvXz+OHz8egNkFnjlX6zpbAZeTk+NKS8/9A6VlZWXMnj2b8vLyAKQSkWBRVlbGe++9x+eff86ePXuoqKjAOUd8fDxt2rShXbt2ZGZm4s+F/k6ePMnKlStJSUlhw4YNQXn5xcy2Oudyatuma+giEpQuuugicnNzyc3NBarKNioqiujo6IAdMzY2ltGjR7NkyRKKi4tZtGhRwI4VCMH334+ISC1iY2MDWuY1jzN69GhWr17NE088EfDj+ZMKXUTkFBdffDGjR4/mwQcf5L333vM6js9U6CIitWjTpg25ublcf/31IfMiqQpdROQMevXqRVRUFCUlJV5H8YkKXUTkDKKiohg+fDirVq3ilVde8TpOnVToIiJncfHFFzNw4ECKi4uD/tKLCl1EpA5XXnklsbGx3HPPPV5HOSsVuohIHcyMa6+9lsWLF7Nt2zav45yRCl1ExAfNmzend+/eTJw4Ea8+YV8XFbqIiI969+7N119/zfz5872OUisVuoiIj2JiYigsLGT69On893//t9dxTqNCFxE5B+np6bRs2ZI777zT6yinUaGLiJyjvLw8fv/73wfdsgAqdBGRc9S0aVN69epFcXGx11H+jgpdROQ89O7dm927d/Pss896HeV/qdBFRM5DbGwshYWF3HfffUHzAqkKXUTkPGVkZJCSksIdd9zhdRRAhS4ickHy8vJYtWoVW7Zs8TqKCl1E5EI0btyYa665hltuuYXKykpPs/hU6GZWYGbbzWynmd1Xy3YzsznV2983s27+jyoiEpx69OjBwYMHefzxxz3NUWehm1k0MBcoBDKBsWaWecqwQiCj+k8JMM/POUVEglZ0dDRDhgzh4YcfZt++fZ7l8OUM/Spgp3Nul3PuBLACGH7KmOHAc67KZqCpmaX4OauISNBq164dXbp0YfTo0Z5liPFhTBvg6xqP9wA9fRjTBvi7/6rMrISqM3iAw2a2/ZzS/p8WwHfn+bWhSnOODJpzGDCzuoZcyJwvOdMGXwq9tmSnrh3pyxicc88Az/hwzLMHMit1zuVc6H5CieYcGTTnyBCoOftyyWUP0K7G47bA3vMYIyIiAeRLoW8BMsysg5nFAWOANaeMWQNMqH63Sy/gB+ecd68MiIhEoDovuTjnys1sKrAeiAYWOuc+MrPJ1dufBl4CBgM7gSNAUeAiA364bBOCNOfIoDlHhoDM2YL1VkoiInJu9ElREZEwoUIXEQkTQV3okbjkgA9zvrl6ru+b2SYzy/Yipz/VNeca43qYWYWZ3VCf+QLBlzmbWT8ze8/MPjKzN+o7o7/58LPdxMzWmtlfq+cc6NfiAsrMFprZt2b24Rm2+7+/nHNB+YeqF2A/B9KAOOCvQOYpYwYDL1P1PvhewNte566HOfcBmlX/vTAS5lxj3GtUvQB/g9e56+H73BT4GGhf/TjZ69z1MOcHgEer/54EHADivM5+AXO+BugGfHiG7X7vr2A+Q4/EJQfqnLNzbpNz7mD1w81Uvec/lPnyfQa4E1gFfFuf4QLElznfBKx2zn0F4JwL9Xn7MmcHXGxVH7NsRFWhl9dvTP9xzr1J1RzOxO/9FcyFfqblBM51TCg51/kUU/U/fCirc85m1gYYCTxdj7kCyZfvcyegmZltMLOtZjah3tIFhi9zfhLoTNWHEj8ApjnnvF2PNrD83l++fPTfK35bciCE+DwfM+tPVaFfHdBEgefLnGcD9zrnKnxYIyMU+DLnGKA7MBBoCLxlZpudc58FOlyA+DLnfOA9YADQEXjFzP7inPsxwNm84vf+CuZCj8QlB3yaj5llAc8Chc657+spW6D4MuccYEV1mbcABptZuXPuP+olof/5+rP9nXOuDCgzszeBbCBUC92XORcB/+qqLjDvNLMvgMuAd+onYr3ze38F8yWXSFxyoM45m1l7YDUwPoTP1mqqc87OuQ7OuVTnXCrwB2BKCJc5+Paz/QLwEzOLMbMEqlY4/aSec/qTL3P+iqrfSDCzlsClwK56TVm//N5fQXuG7oJzyYGA8nHOPwcSgaeqz1jLXQivVOfjnMOKL3N2zn1iZuuA94FK4FnnXK1vfwsFPn6fHwIWm9kHVF2OuNc5F7LL6prZcqAf0MLM9gD/D4iFwPWXPvovIhImgvmSi4iInAMVuohImFChi4iECRW6iEiYUKGLiIQJFbqISJhQoYuIhIn/D97MEUjFlxanAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from matplotlib.pylab import plt\n", "import numpy as np\n", "plt.rcParams[\"figure.figsize\"] = (9,6)\n", "\n", "fig = plt.figure()\n", "x_grid = np.linspace(0, 1., 1000)\n", "plt.plot(x_grid, (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), 'k-', linewidth=1)\n", "plt.fill_between(x_grid, -1., (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), color='gray')\n", "plt.ylim(0, None);" ] }, { "cell_type": "markdown", "id": "a02d74cf", "metadata": {}, "source": [ "We want the particle to bounce off both the wavy ramp and the ground (i.e., $x$ axis). Thus, we will have to define two different event equations and two different callback functions to implement the bouncing behaviour. Let us begin with the event equations:" ] }, { "cell_type": "code", "execution_count": 3, "id": "04c9f0e8", "metadata": {}, "outputs": [], "source": [ "# Event equation for the ramp.\n", "eq_w_curve = y - (1. - x + 0.05 * hy.cos(11 * np.pi * x))\n", "\n", "# Event equation for the ground.\n", "eq_bottom = y" ] }, { "cell_type": "markdown", "id": "4479fba2", "metadata": {}, "source": [ "Next, we need to implement the callbacks. Because we are dealing with inelastic collisions, when the particle hits a surface it will lose part of its kinetic energy. The [coefficient of restitution](https://en.wikipedia.org/wiki/Coefficient_of_restitution) $C_R \\in \\left[0, 1\\right]$ regulates how much energy is lost at every bounce.\n", "\n", "Additionally, because of the loss of energy, the particle will experience increasingly frequent collisions if it gets stuck in one of the local minima of the curve or when it starts bouncing off the ground. In practice, we want to stop the simulation once the collisions become too frequent.\n", "\n", "Let us see the code:" ] }, { "cell_type": "code", "execution_count": 4, "id": "c6a7ab81", "metadata": {}, "outputs": [], "source": [ "# The coefficient of restitution.\n", "CR = .8\n", "\n", "# Global variables to track\n", "# the time of the last collision\n", "# and the collision points.\n", "last_t = 0.\n", "bounce_points = []\n", "\n", "# Callback for bouncing against the curve.\n", "def cb_w_curve(ta, mr, d_sgn):\n", " global last_t\n", "\n", " # If the last collision happened\n", " # too recently, return False to\n", " # stop the simulation.\n", " if ta.time - last_t < 1e-6:\n", " return False\n", "\n", " # Update last_t.\n", " last_t = ta.time\n", "\n", " # Update bounce_points.\n", " x, y = ta.state[0:2]\n", " bounce_points.append((x, y))\n", " \n", " # Compute the normal unit vector\n", " # using the gradient of the event\n", " # equation.\n", " grad = np.array([1+0.05*11*np.pi*np.sin(11*np.pi*x), 1])\n", " grad_uvec = grad / np.linalg.norm(grad)\n", " \n", " # Compute the component of the velocity\n", " # across the normal vector.\n", " xy_vel = ta.state[2:4]\n", " vproj = np.dot(xy_vel, grad_uvec)\n", " \n", " # Flip it and rescale it according\n", " # to the coefficient of restitution.\n", " Dv = -vproj*grad_uvec\n", " xy_vel += (1. + CR) * Dv\n", "\n", " return True\n", "\n", "# Callback for bouncing off the ground.\n", "def cb_bottom(ta, mr, d_sgn):\n", " global last_t\n", " \n", " # If the last collision happened\n", " # too recently, return False to\n", " # stop the simulation.\n", " if ta.time - last_t < 1e-6:\n", " return False\n", "\n", " # Update last_t.\n", " last_t = ta.time\n", "\n", " # Update bounce_points.\n", " x, y = ta.state[0:2]\n", " bounce_points.append((x, y))\n", " \n", " # Flip the y component of the velocity,\n", " # and rescale it according\n", " # to the coefficient of restitution.\n", " ta.state[3] = -CR*ta.state[3]\n", " \n", " return True" ] }, { "cell_type": "markdown", "id": "ba77fb48", "metadata": {}, "source": [ "We are now ready to create the integrator object. As initial conditions, we will put the particle at rest above the ramp on the $y$ axis:" ] }, { "cell_type": "code", "execution_count": 5, "id": "285b4776", "metadata": {}, "outputs": [], "source": [ "ta = hy.taylor_adaptive(eqns, [0., 1.2, 0., 0.], t_events = [hy.t_event(eq_w_curve, callback = cb_w_curve, direction = hy.event_direction.negative),\n", " hy.t_event(eq_bottom, callback = cb_bottom, direction = hy.event_direction.negative)])" ] }, { "cell_type": "markdown", "id": "e7b95f4e", "metadata": {}, "source": [ "Note how, as explained in the [Keplerian billiard](<./The Keplerian billiard.ipynb>) example, we assigned specific directions to the collision events. This is done in order to avoid triggering spurious double-bounce events that would lead to the particle penetrating through the ramp or the ground.\n", "\n", "We can now integrate the system for a few time units:" ] }, { "cell_type": "code", "execution_count": 6, "id": "a1fab993", "metadata": {}, "outputs": [], "source": [ "t_grid = np.linspace(0, 10, 10000)\n", "oc, _, _, _, res = ta.propagate_grid(t_grid)\n", "\n", "# Transform bounce_points in a NumPy array\n", "# for ease of use.\n", "b_points = np.array(bounce_points)" ] }, { "cell_type": "markdown", "id": "1c503814", "metadata": {}, "source": [ "Let us take a look at the trajectory:" ] }, { "cell_type": "code", "execution_count": 7, "id": "cff0eff9", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsIAAAIICAYAAABkYYgLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAB400lEQVR4nO3dd3gc1dk28Pts06qXberNknvvxgVsY4wBx/QYQglJ6CXkSwLhTQK8SUh9U0iAhB4goRdjwIDBGNyLZFuyLNmWbEtW771tm++PlWSVlbSSZou19++6uOKdnTlztFnM7eNnniMkSQIRERERkb9ReHsCRERERETewCBMRERERH6JQZiIiIiI/BKDMBERERH5JQZhIiIiIvJLDMJERERE5JdU3rqxXq+XkpOTvXV7IiIiIvITmZmZNZIkGfof91oQTk5ORkZGhrduT0RERER+QghR5Ow4SyOIiIiIyC8xCBMRERGRX2IQJiIiIiK/xCBMRERERH6JQZiIiIiI/BKDMBERERH5JQZhIiIiIvJLDMJERERE5JcYhImIiIjILzEIExEREZFfYhAmIiIiIr/EIExEREREfolBmIiIiIj8EoMwEREREfklBmEiIiIi8ksMwkRERETklxiEiYiIiMgvMQgTERERkV9iECYiIiIiv8QgTERERER+iUGYiIiIiPwSgzARERER+SUGYSIiIiLySwzCREREROSXGISJiIiIyC8xCBMRERGRX2IQJiIiIiK/xCBMRERERH6JQZiIiIiI/BKDMBERERH5JQZhIiIiIvJLDMJERERE5JcYhImIiIjILzEIExEREZFfYhAmIiIiIr80bBAWQrwkhKgSQuQM8v53hBDZXf/sEULMkn+aRERERETycmVF+N8ALh3i/TMALpQkaSaAXwN4ToZ5ERERERG51bBBWJKkHQDqhnh/jyRJ9V0v9wGIl2lusrPa7ChvbEeb2ertqRARERGRl8ldI/x9AJ/KPKZsqpo7seR3X+GjrDJvT4WIiIiIvEwl10BCiJVwBOFlQ5xzB4A7ACAxMVGuWxMRERERjZgsK8JCiJkAXgCwQZKk2sHOkyTpOUmS5kuSNN9gMMhxayIiIiKiURlzEBZCJAJ4H8DNkiSdHPuUiIiIiIjcb9jSCCHEGwAuAqAXQpQAeAyAGgAkSfoXgEcB6AA8I4QAAKskSfPdNWEiIiIiIjkMG4QlSbphmPd/AOAHss2IiIiIiMgDuLMcEREREfklBmEiIiIi8ksMwkRERETklxiEiYiIiMgvMQgTERERkV9iECYiIiIiv8QgTERERER+iUGYiIiIiPwSgzARERER+SUGYSIiIiLySwzCREREROSXGISJiIiIyC/5ZRCWJG/PgIiIiIi8za+CsBDengERERER+Qq/CsJERERERN0YhImIiIjILzEIExEREZFfYhAmIiIiIr/EIExEREREfolBmIiIiIj8EoMwEREREfklBmEiIiIi8ksMwkRERETklxiEiYiIiMgvMQgTERERkV9iECYiIiIiv8QgTERERER+iUGYiIiIiPwSgzARERER+SUGYSIiIiLySwzCREREROSX/DIIS96eABERERF5nV8FYQHh7SkQERERkY/wqyBMRERERNSNQZiIiIiI/BKDMBERERH5JQZhIiIiIvJLDMJERERE5JcYhImIiIjILzEIExEREZFfYhAmIiIiIr/EIExEREREfolBmIiIiIj8EoMwEREREfklBmEiIiIi8ksMwkRERETklxiEiYiIiMgvMQgTERERkV/yqyCszXsPuzQPYOOWmcBfpwPZb3t7SkRERETkJSpvT8Bjst9G2Jc/RoSi3fG6sRj46AHHr2de7715EREREZFX+M+K8LZfQWFt73vM0g5s+5V35kNEREREXuU/QbixxOlhe2MJduXXeHgyRERERORt/hOEw+OdHq6EHje9uB83v7gf+ZXNHp4UEREREXmL/wTh1Y8C6sC+x9SB0F/5G/zi8inILmnE5X/fhae+yofFZvfOHImIiIjIY/wnCM+8Hlj/dzSKcNgloMYaCKz/O9SzN+IHy1Ox7ccXYs00E/5v60nc8Nw+VDZ1eHvGRERERORG/hOEAWDm9Xgu6H7cVngl5v9H26dbhD4kAE/fOBdPbpyN3PImXP73Xdh/utaLkyUiIiIid/KvINwlMTER5eXlKCsrG/Dehtlx2HTvUoRpVbj5xQPYcrTcCzMkIiIiInfzyyCsUqmQkJCADz/80On7E02h+OCepZgZH457Xz+E1/YVeXiGRERERORufhmEASAhIQHbtm0b9P3wIDVe+/4irJ5sxC835eCVPYWemxwRERERuZ3fBuGkpCRkZmYOeU6gRol/3jQPa6aa8NjmY3jr4FkPzY6IiIiI3M1vg3BcXBxKS0vR0NAw5HlqpQJP3TgHF0404GfvH2XNMBEREdE44bdBWKPRwGg0Dlke0S1ApcSzN8/D3MRI/OitIzhS3OD+CRIRERGRW/ltEAaA2NhY7Ny506VztWolnrt5HoxhAfjBKxkobWh38+yIiIiIyJ38OghHR0cPWyfcmy4kAC/dugCdFht+8EoG2s02N86OiIiIiNzJr4NwbGws8vPzR3RNuikU/7hxDo5XNOHxzcfcNDMiIiIicje/DsIGgwF1dXWoqakZ0XUXTTLi3ovS8FZGMd7NLHHT7IiIiIjInYYNwkKIl4QQVUKInEHeF0KIvwshCoQQ2UKIufJP0z1UKhVMJhO2b98+4msfvDgdi1Oj8ItNR3GystkNsyMiIiIid3JlRfjfAC4d4v11ANK7/rkDwD/HPi3PiYmJwZ49e0Z8nUqpwN83zkGwRoUfvXUEZqvdDbMjIiIiIncZNghLkrQDQN0Qp2wA8KrksA9AhBAiRq4JupvRaER2dvborg3T4rdXz8CxsiY8tb1A5pkRERERkTvJUSMcB6C41+uSrmMDCCHuEEJkCCEyqqurZbj12BmNRpw6dWrU16+dFo2r58bh6e0FyGJ/YSIiIqLzhhxBWDg5Jjk7UZKk5yRJmi9J0nyDwSDDrcfOaDSivLwcVqt11GM8tn4ajKEB+PE7Wei0sqUaERER0flAjiBcAiCh1+t4AGUyjOsRgYGBCAwMxOHDh0c9RnigGr+7egYKqlrw3DenZZwdEREREbmLHEF4M4BburpHLAbQKElSuQzjeozJZMK+ffvGNMZFk4y4fGYM/rG9AIU1rTLNjIiIiIjcxZX2aW8A2AtgkhCiRAjxfSHEXUKIu7pO2QLgNIACAM8DuMdts3UTg8EwphXhbo9eMRUBSgV++WEOJMlpdQgRERER+QjVcCdIknTDMO9LAO6VbUZeYDQacezY2HeJM4Vp8ZO1k/DY5mP4KLsc35oVK8PsiIiIiMgd/HpnuW4mkwlFRUWyjHXT4iTMiAvHE5/kos08+gfwiIiIiMi9GIQB6PV61NbWorV17LW9SoXAY+unorKpE8/t4INzRERERL6KQRiOrZYjIyNx4MABWcabnxyFy2fE4NlvTqOisUOWMYmIiIhIXgzCXUwmEw4ePCjbeD9bNxk2u4T/23pCtjGJiIiISD4Mwl30ej2ysrJkGy8hKgi3LUvGe4dKkFPaKNu4RERERCQPBuEuBoMBx48fl3XMe1emITJIgz98Ju+4RERERDR2DMJdjEYjzp49K+uYYVo17rloAnbm12Df6VpZxyYiIiKisWEQ7hIVFYWGhgY0NzfLOu5Ni5NgDA3An7ee4CYbRERERD6EQbiLUqmEXq+XrXNEN61aiftXpeFgYT125NfIOjYRERERjR6DcC9GoxEZGRmyj/vtBYmIiwjkqjARERGRD2EQ7kXuzhHdNCoFfnhxOrJLGvFlXpXs4xMRERHRyDEI9+KOzhHdrp4Th4SoQDy1vYCrwkREREQ+gEG4F6PRiKKiIreMrVIqcNeFE5BV3IC9p9hBgoiIiMjbGIR7iYyMRHNzs+ydI7pdMzcextAAPP11gVvGJyIiIiLXMQj3olAooNfrsX//freMr1UrcfvyVOwuqMXhs/VuuQcRERERuYZBuJ/o6GgcPHjQbePfuCgR4YFqPPP1Kbfdg4iIiIiGxyDcj16vR3Z2ttvGDw5Q4balyfgitxInKtxTgkFEREREw2MQ7kev1yMvL8+t9/juBckIVCvx4q7Tbr0PEREREQ2OQbgfd3aO6BYRpMG18+Kx6UgZalo63XovIiIiInKOQbif8PBwtLW1ob7evQ+zfXdpMsxWO/6zz72hm4iIiIicYxDuR6FQwGAwuGWr5d4mGEKwarIR/9lXhE6rza33IiIiIqKBGISdiI6OxoEDB9x+n+8tTUFNixmbj5S5/V5ERERE1BeDsBMGgwFZWVluv8/SNB0mmULx0u5CbrtMRERE5GEMwk7odDocP37c7fcRQuB7y5KRV96Evae57TIRERGRJzEIO+GJzhHdNsyOQ2SQmg/NEREREXkYg7ATYWFhMJvNqKmpcfu9tGolrpufgK3HKlHV1OH2+xERERGRA4OwE0IIGI1G7Nu3zyP3u2FhIqx2CW8dLPbI/YiIiIiIQXhQJpPJI50jACBFH4zl6Xq8ceAsbHY+NEdERETkCQzCg9DpdMjOzvbY/b6zKBFljR34+kSVx+5JRERE5M8YhAdhMBhw4sQJj91v9RQTjKEBfGiOiIiIyEMYhAdhNBpx9uxZj/X3VSsV2LgwEV+frEZxXZtH7jmcotpWvJtZAjvLNYiIiGgcYhAeREhICCRJwtmzZz12z40LEiAAvHHAc/ccjCRJ+P4rGfjJO1n45Gi5t6dDREREJDsG4UEIIRAdHY2dO3d67J6xEYFYNdmIdzNLYLXZPXZfZ4rr2lFQ1QIA+CK30qtzISIiInIHBuEhREdHY8+ePR6957XzElDV3Imd+e7vYTyUIyUNAIBUfTAyi+q9OhciIiIid2AQHkJsbCwyMjI8es9Vk42ICtbgnUzv9hQ+VtoIjUqBK2bGoLShHW1mq1fnQ0RERCQ3BuEhxMXF4eTJkx69p0alwJWz4/BlbhXqW80evXdvZ+vakBAZiMkxYQCA09WtXpsLERERkTswCA8hMjISFovF42H4uvnxMNvs+PBIqUfv21tJfTviI4OQrAsG4AjGREREROMJg/AQhBBISEjAZ5995tH7TokJw7TYMLyTWeLR+/ZW2tCOuMhAxIRrAQDljR1emwsRERGROzAIDyM2Nha7du3y+H2vmxePY2VNyC1r8vi9WzutqGs1Iz4yEBFBagSoFKhobPf4PIiIiIjciUF4GImJiTh48OCIr7PZbMjKykJHx+hWUjfMjoNGqfDKQ3PlXaE3LiIQQgjEhGu5IkxERETjDoPwMOLj41FWVoaqqiqXr9mxYwfi4uKwbNkyxMTEYOvWrSO+b2SwBhdPNeLDI2UwWz3bU7imxfGQnj4kAAAQHa5FBYMwERERjTMMwsNQq9VISkrCu+++69L5R48exfr167Fo0SL85Cc/wcUXX4zrr78eJSUjr/e9ek486lrN2JlfPeJrx6Kuq1tFVLAGAGAM1aKqudOjcyAiIiJyNwZhFyQnJ2PLli3Dnmez2bBx40bMmzcPc+fOBQBMnz4d6enpeOCBB0Z83xUTDYgIUuPDI2UjvnYsaruCsK4rCEcFa1Df5r1WbkRERETuwCDsggkTJmDv3r2w24cuUXj88cdRV1eHpUuX9jm+YsUKfPrppygvLx/RfTUqBS6fEYMvcivR2um5DS1qWxyrv5FdQTgySIPmDissXt72mYiIiEhODMIuMJlMUCgUQ9b65ufn469//SvWr18PlUrV572wsDCkpaXhL3/5y4jvvWF2HNotNnyRWznia0errtWM8EA11ErH1yMqWA0AXBUmIiKicYVB2AVCCEybNg3PP//8oOfcfPPNmDNnDuLi4py+P3v2bLzzzjsjvvf8pEjERQR6dHON2lZzT1kEcG5luL7V4rE5EBEREbkbg7CLZsyYga1bt6KlpWXAe08//TQKCgpw4YUXDnp9SkoKKisrcfz48RHdV6EQWD8rFjvya3pKFtytrsXc86AccO6huTovbvlMREREJDcGYRfp9XrExcXhiSee6HO8qKgIjzzyCDZs2AC1Wj3o9UqlEunp6XjllVdGfO8r58TCZpew5ejIaoxHq7a1E7qQgUGYpRFEREQ0njAIj8CFF16Ip556CmfOnAEAtLe3Y926dZg9ezYSExOHvX7ixIn49NNPR3zfydFhmGQKxSYPdY9oaLMgIrBXEA7iijARERGNPwzCIxAbG4sFCxZg2bJlePzxxzFz5kwoFAqsXLnSpetTUlJw4sQJmM0jD5Qb5sQis6gexXVtI752pJo7rAgLPPfAX0RQd40wgzARERGNHwzCI7RixQosWrQImzZtwqRJk3DNNddAoXDtYwwJCUFYWBi+/PLLEd93/cxYAMDH2e4tj7DY7Gi32BCqPVfmoVEpEKBSoNmDLdyIiIiI3I1BeISEEJg1axauuuoqLFiwwOUQ3C0pKWlU5REJUUGYFR+OT3PcG4SbOxxhN0zbtwVcqFaN5g52jSAiIqLxg0HYwxITE7F79+5RXXvZjBhklzS6tTyiqd0RdnuvCAOOYNzUwRVhIiIiGj8YhD0sKSkJJ06cGHaXOmcumxEDAG7tHtGzIhzYNwiHBqp73iMiIiIaDxiEPSwsLAxKpRJHjx4d8bUJUUGYEReOLTkVbpiZQ1NH94pwv93xtCqWRhAREdG4wiDsBXFxcdi2bduorr1sRgyyihtQUu+e8ojusBvWrzQiVKvqKZsgIiIiGg8YhL0gOjoae/fuHdW1l82IBgB8etQ9q8JN7Y7yh/4rwqEBLI0gIiKi8YVB2AtiY2ORk5MzqmuTdMGYFhuGT9xUJ9w0xIowgzARERGNJwzCXhAbG4vCwsJRPTAHOMojjhQ3oLShXeaZnXtYLsRJ+7R2iw0W2+jmTERERORrGIS9ICQkBFqtFpmZmaO6vrt7xKduWBVu6rAgJEAFpUL0Od5dKtHCVWEiIiIaJxiEvSQ6OnrU/YRT9MGYEhPmljZqzR3WAZtpAOfaqTWxcwQRERGNEwzCXmIwGHDo0KFRX79uejQOFzegqrlDxlk5NtTov5kGcG5FmHXCRERENF4wCHuJyWRCbm7uqK9fM9UESQK25VXJOCug1WxFcIBywPGQAEcQbu1kECYiIqLxgUHYS4xGI4qKikZ9/eToUCREBWLrMXnbqLV22hAcMLA0IkjjCMdtZpus9yMiIiLyFgZhL9HpdGhsbERtbe2orhdCYM2UaOw+VSvrKm2b2doTenvrDsetZq4IExER0fjAIOwlSqUSer0ee/bsGfUYa6aaYLbaseNktWzzajPbEKwZYkW4kyvCREREND64FISFEJcKIU4IIQqEED9z8n64EOIjIUSWEOKYEOI2+ac6/phMJhw8eHDU1y9IjkREkBpbcytlm1Ob2YYgJzXC3eGYK8JEREQ0XgwbhIUQSgBPA1gHYCqAG4QQU/uddi+AXEmSZgG4CMCfhRAamec67hgMBmRlZY36epVSgVWTjfjqeJVsG104SiOcrAgHsEaYiIiIxhdXVoQXAiiQJOm0JElmAG8C2NDvHAlAqBBCAAgBUAeAS4fDMBqNOHHixJjGuGSqCY3tFhwsrBvzfGx2CR0Wu9MaYY1SAZVCsGsEERERjRuuBOE4AMW9Xpd0HevtKQBTAJQBOArgh5IkDViiFELcIYTIEEJkVFfLV9d6vjKZTCgpKRnTGCsmGhCgUmDrsbGXR7R1lT04qxEWQiBQo+SKMBEREY0brgRh4eSY1O/1WgBHAMQCmA3gKSFE2ICLJOk5SZLmS5I032AwjHCq409oaChsNtuY2qgFaVRYlqbHF7mVkKT+/7eMTHtXyA10siIMOAJyG2uEiYiIaJxwJQiXAEjo9ToejpXf3m4D8L7kUADgDIDJ8kxx/BJCwGQyYe/evWMaZ81UE0ob2pFb3jSmcVq7grCzDTUAR51wK1eEiYiIaJxwJQgfBJAuhEjpegBuI4DN/c45C2A1AAghTAAmATgt50THK6PROKatlgFg9RQThAC+zB3bLnPdq73OHpYDulaEWSNMRERE48SwQViSJCuA+wB8DiAPwNuSJB0TQtwlhLir67RfA7hACHEUwDYAD0uSVOOuSY8ner0eOTk5YxrDEBqAmfER2H5irEHYsdrr7GG57uNcESYiIqLxwvnSXz+SJG0BsKXfsX/1+nUZgEvknZp/MBqN2L1795jHWTnJgCe35aO2pRO6kIBRjXEuCA+yIhygQlVzx6jnSERERORLuLOclxmNRpSUlIz5QbdVk42QJGBH/ui7cXSXPQy1Isyd5YiIiGi8YBD2suDgYCiVSpw5c2ZM40yPDYc+JABfHR99EO55WG6oGmGWRhAREdE4wSDsA+ToHKFQCFw0yYAdJ6thHeUuc+3dD8sN0jUiUKPkFstEREQ0bjAI+wA5OkcAwMpJRjS2W3C4uGFU17cO87BccIBjQ42xlnEQERER+QIGYR8gR+cIAFg+UQ+lQmD78dF1j2gz2yAEoFUNViOsgs0uodM6uhVnIiIiIl/CIOwDjEYj8vPzxzxOmFaN+UmR+Gq0QbjTiiC1EgqFs80Ez60Ut7NOmIiIiMYBBmEfYDAYUFZWBrt97CutqyYbcbyiGeWN7SO+ttVsQ+AgD8oBgFbtCMIdVgZhIiIiOv8xCPuAwMBABAQE4NSpU2Mea+VkIwDg6xMj7x7RbrYOur0yAASquSJMRERE4weDsI8wmUzYt2/fmMdJN4YgLiJwVOURbWZbT9h1Rqt2fF06LKwRJiIiovMfg7CPMJlMsnSOEEJg1WQjdhfUoHOEJQztFltP+YMz3e+1W7giTEREROc/BmEfodPpcPToUVnGWjnZgDazDQfP1I/ouk6LvWfV15nuINzJIExERETjAIOwjzAYDCgoKJBlrMWpOqiVAjtHuN1yh3XoFeFArggTERHROMIg7CO6O0fYbGMPmUEaFeYnRWFHfs2Iruuw2AbtIQz06hrBGmEiIiIaBxiEfURAQABCQkJw8uRJWcZbPlGPvPImVDV3uHxNh8WOwEF2lQO4IkxERETjC4OwD5GrcwQArEg3AAB2F7i+KtxhsQ1TI6zoOY+IiIjofMcg7ENMJhMyMzNlGWtqTBiigjXYeXJkQThgqNIIjbLnPCIiIqLzHYOwD4mKisKxY8dkGUuhEFiWpseO/BpIkuTSNR1W+9Dt01QMwkRERDR+MAj7EKPRiPz8fNnGWzHRgJqWThyvaB72XJtdgtk6dPs0tVJAqRCsESYiIqJxgUHYh+j1elRWVsJsNssy3vJ0PQC41Eate/ONoVaEhRDQqhTsGkFERETjAoOwD1Gr1QgPD0dubq4s45nCtJhkCsUOF+qEu8OtVjX0VyJQo+SKMBEREY0LDMI+Jjo6Gnv37pVtvOXpehworEO7eejw2l33O9SKMAAEqJSsESYiIqJxgUHYxxiNRhw6dEi28ZZPNMBsteNAYd2Q57kahAM1DMJEREQ0PjAI+xidToecnBzZxluYHAWNSoGdJ4euE+4pjRjiYbnu91kjTEREROMBg7CPMRqNOH36tGzjBWqUWJgchZ3DbLfc0fWwXMBwK8Jq5bBlFkRERETnA5W3J0B96XQ61NXVob29HYGBgbKMuSxdj99/ehxVTR0whmmdntNTGjHEhhqAo3SipdM67D0lScKXeVV488BZVDR1YHJ0GG5fkYLJ0WEj/wGIiIiI3IArwj5GqVRCp9PJttUyACyd4Gijtvd07aDndAfhQM3wQXi4FWGbXcLP3juK21/NwPGKZhhCA7D1WAWu+PsuvJNRPMLZExEREbkHg7APMhqNsnaOmBobhvBANfYUDBWEXa0RVqLTOnSN8K8+Ooa3Mopxz0UT8PVPL8K/b1uIHQ+txOJUHX76bjY+y6kY+Q9BREREJDMGYR9kMpmwf/9+2cZTKgQWp0Zhz+nB64RdLY0IVCuGXBH+JLscr+wtwg+WpeChSydDrXR8xSKDNXjh1vmYnRCBn7yThbKG9lH8JERERETyYRD2QQkJCcjKypJ1zAsm6FFc147iujan759bER7+YbnuB+v6a+6w4PGPjmFmfDgeXjd5wPtatRL/uGEObHYJv9wkX2cMIiIiotFgEPZBsbGxKCsrQ0tLi2xjXjBBBwDYc8r5qvC5PsLDl0YMtiL81PYC1LR04tcbpvesBPeXEBWEBy9Ox7bjVdh7avBSDSIiIiJ3YxD2QWq1GkajEV9++aVsY6YZQ6APCcCeQcJn9yrvsDvLddUI2+1Sn+P1rWa8trcIG2bFYlZCxJBj3HpBMkxhAfi/rScgSdKQ5xIRERG5C4Owj0pISJA1CAshcMEEHfacqnUaPrtLIwJUQ38lut832/o+MPfynkK0mW24Z2XasHPRqpW4d2UaMovqkVlU7+qPQERERCQrBmEfFRsbO6oWalVVVfjggw/Q0NAw4L0LJuhQ3dyJgqqBJRedFhsCVAoIIYYcvzsI9+4c0Wm14T/7inDxFBMmmkJdmue18+IRqlXh33sKXTqfiIiISG4Mwj4qOTkZubm5sFgsLl/z9NNPIyUlBffeey+Sk5Oxffv2Pu8vTXP0E3ZWHtFhsQ1bFgGc23mus9cDc1/kVqKu1YybFie6PNcgjQobFyTg05wKVDR2uHwdERERkVwYhH1UWFgYwsPD8dlnn7l0/ltvvYWHH34YN954I+68806sWrUK1113Herq6nrOSYgKQnxkoNMH5tottmEflAN6lUb0WhF+62Ax4iICsTzd4NJcu31nURJsdgmbjpSO6DoiIiIiOTAI+7C0tDS8//77w55XVlaGO+64A1dddRXi4+MBALNmzUJ0dDQeeeSRPudeMEGHfafrYOv3sFuHxe7ainC/0ojShnbszK/B9fMToFQMXVbRX7I+GHMSI7DpMIMwEREReR6DsA9LTU3F119/Pex5t956KyZNmoS0tL4Pqq1YsQKvv/56nzZsF0zQo7Hdgrzypj7ndlhsCBxJEO56uO7To+UAgCvnxA57rTNXzYnD8YrmAfMhIiIicjcGYR+WmJiIqqoqHDt2bNBz3nrrLRw4cACrV68e8J7JZILRaMSzzz7bc2xJVz/h3QV9yyM6rPae+t+hBKj61ghvOVqOabFhSNIFD/8DOXH5jBioFAKbs8pGdT0RERHRaDEI+zCVSoXp06fjySefdPp+a2sr7r//fqxbtw5ardbpOTNnzsRrr73W89oUpsUEQzB293tgrsNig3aY1mlA3xrhsoZ2HDrbgMtmxLj6Iw2gCwnA4lQdvsitHPUYRERERKPBIOzjZsyYgQ8++ABWq3XAe3fffTeMRiOmTJky6PWTJ09GXl4eamvPBd8lE3TILKyD1da7BZqLK8LqczXCn+VUAADWTY92+edx5uIpRhRUteBMTeuYxiEiIiIaCQZhHxcfH4+goCD87W9/63N869ateO+997B27dohr9dqtYiPj8cbb7zRc2xRig6tZhuOlZ2ry+3uIzwcjbK7NMKO7SeqkGYMQaohZAQ/0UCrp5gAANvyuCpMREREnsMg7OOEELjooovwxBNPoLzc8WDaqVOnsHHjRlx++eUICwsbdoy0tDRs3ry55/WilCgAwP4z51aJzTY7NK6URnStCNe3mbH/dB0umjiylmnOJEQFYXJ0KLayPIKIiIg8iEH4PJCamorp06dj3rx5uP322zF//nzMnz8f06ZNc+n6tLQ0ZGRk9GytbAzTIkUfjANnzvUYNlvtCFC6XiO842Q1zDY7Lpw09iAMABdPMSGjsA5NHa5vIEJEREQ0FgzC54nVq1dj2bJlyMvLw4YNG7B06VKXrzUYDOjs7OzTfWJRShQOnDnXT9hRI+xKEHaURmzNrUSgWokFyVEj/EmcW5qmh10C9p+uG/5kIiIiIhkwCJ8nhBCYNm0a1qxZg5SUlBFfm5SUhI8//rjn2MKUKDR1WHGiohmAY0VYM4IVYbPVjiUTdC5twuGKuUkR0KoVA9q6EREREbkLg7CfSEhIwDfffNPzelGqo59wd52w2epajXDvc5an62WbX4DKsbrMIExERESewiDsJxITE5GVldXzOi4iEPGRgT2lCJ1WW0/Zw1B6d5ZY3BWm5bI0TY/8qhZUNXXIOi4RERGRMwzCfiI6Oho1NTWoqTm34rowJQoHCutgsdlhl+DSirCqV/nEJFOorHNcluZYYd59iqvCRERE5H4Mwn5CpVLBZDLhq6++6jm2OEWHulYzcrv6CbsShHtTKISsc5waE4ZQrQoHztTLOi4RERGRMwzCfiQmJga7d+/ueb2wq5/wzvxqAHBpQ42qZveVLSgUAnMTI5FZxM4RRERE5H4Mwn4kOjoahw4d6nmdpAuCKSwAO/MdpQiurAgf7Fqtlbssotv8pEicrGxBYxv7CRMREZF7MQj7kdjYWOTn5/e8FkJgUYoO+7s21nClfdqBri4Taaaxbas8mHlJkQCAQ8UsjyAiIiL3YhD2IwaDAXV1dairO1d60F0eAQABLvQEzihyBNTuXerkNishAkqFQGYhgzARERG5F4OwH1GpVDAajfj66697ji1OPReEh1sR7rDYcLxrA45Oi90tcwwOUGFKTCgyWCdMREREbsYg7Geio6Nx4MCBntcTDOdKHIZ7WO5YWWOfLZndZX5SFLKKG2Gxue8eRERERAzCfsZgMCA7O7vntRACumANgOGD8JHiRgDABEMwOq02t81xTmIE2i02nKxsdts9iIiIiBiE/YzRaOzzwBwAzIwPBwDUtZmHvDaruAEx4VokRAW5dUV4ZnwEACCntNFt9yAiIiJiEPYzRqMRpaWlsNvPBdkZXcHz6DDBM7ukATPjwxGgUsDsxiCcFBWE0AAVsksYhImIiMh9GIT9THBwMIQQOHHiRM+xCYZgAMDRIYJnQ5sZhbVtmJUQgQCV0q0rwgqFwPS48GGDOREREdFYMAj7GSEEoqOjsXfv3p5j9q5WaEMF4ayu92bHR0CjUqDT4r4aYcBRrnG8vNmtK89ERETk3xiE/ZDBYOizw1x32GzutKK10+r0muziBggBTO8qjXDnijAATI8Lh9lm5wNzRERE5DYMwn7IYDAgJyen53XvVdcjxQ1Orzla2ogUXTDCtGoEqJRuX6ntfoCPdcJERETkLgzCfshoNOL06dM9r3uv7mYMsqNbXkUTpsSEAQAC1O5fEU6MCkKYVsU6YSIiInIbBmE/ZDQaUVFRAavVUQbRHWpT9cFOd3Rr7rCguK4dU2JCATh2oDPb7LDb3bPNMuCoZZ4RH46jpQ1uuwcRERH5NwZhP6TVahEUFISsrCwA50ojlkzQ4fDZhp7d47p1b6vce0UYAMxu3vltakwYTla2wMod5oiIiMgNXArCQohLhRAnhBAFQoifDXLORUKII0KIY0KIb+SdJsnNZDL1bLXcabVDo1RgQXIUWjqtOF7R1OfcvHLH654grFI6rrO4N6BOig6D2WpHYW2bW+9DRERE/mnYICyEUAJ4GsA6AFMB3CCEmNrvnAgAzwD4liRJ0wBcJ/9USU56vR5HjhwB4FgR1qgUmJcUCQDILOpbJ5xX3oTwQDViwrUAAI1SOK5z80rt5GhHKUb/YE5EREQkB1dWhBcCKJAk6bQkSWYAbwLY0O+cGwG8L0nSWQCQJKlK3mmS3AwGA3JzcwEAZpsNASoF4iMDYQoLGPDAXG55M6bEhEIIRwDWqDxTGpFmDIFSIXC8nC3UiIiISH6uBOE4AMW9Xpd0HettIoBIIcTXQohMIcQtck2Q3MNgMKCwsBCAo8RBo1JACIH5yVHIKDz3wJzNLuFEr44RAKBWOr42Fjd3jtCqlUjVB3NFmIiIiNzClSAsnBzr3y5ABWAegMsBrAXwSyHExAEDCXGHECJDCJFRXV094smSfAwGAyorK2G1WmG22XtWeecnRaKssQNlDe0AgMLaVnRY7H2CcPe5Fg88xDYpOrTnYT0iIiIiObkShEsAJPR6HQ+gzMk5n0mS1CpJUg2AHQBm9R9IkqTnJEmaL0nSfIPBMNo5kwwCAgIQHByM7OxsmK12BPQE4SgAQEZXnXD3g3JTnawIu7uXMOB4QK+kvh1NHRa334uIiIj8iytB+CCAdCFEihBCA2AjgM39zvkQwHIhhEoIEQRgEYA8eadKcjOZTNi/f7+ja0RXEJ4SE4ogjRKZXeURx8uboVQIpBlDeq7z5Ipw9wNzJ7kqTERERDIbNghLkmQFcB+Az+EIt29LknRMCHGXEOKurnPyAHwGIBvAAQAvSJKUM9iY5BsMBgOysrIcXSO6VnlVSgVmJ0TgYNcDc/lVzUjSBUGrVvZc132uu7dZBoDJXSvRLI8gIiIiualcOUmSpC0AtvQ79q9+r/8E4E/yTY3cTa/XIzc3F9Gz7T29gQFHnfBT2wvQ2mlFflUL0nutBgO9HpazuW9nuW6x4VqEalV8YI6IiIhkx53l/JjBYMDp06fRabX1lDsAwJykSNglRz/hoto2pBtD+1znydIIIQTSjSHIr2xx+72IiIjIvzAI+7HuzhHtneY+QXh2fAQA4IPDpbDZpT71wQCg7tpQwxMPywGOfsKnqhmEiYiISF4Mwn5Mo9EgNDQUjRXFPV0jACAyWINUfTA+OFwKAAOCcIAHV4S771/TYkZDm9kj9yMiIiL/wCDs50wmE+pKTvVZEQaA2YkRAAAhgAkG5zXCnnhYDjgXxAuquCpMRERE8mEQ9nMGgwEtFWf6rAgDwJzESACAJAGBGmWf9zxZIwwAaQZHjTKDMBEREcmJQdjP6fV6tFUV9ekaAQBzEiIGveZc1wjPBOG4yEBoVArWCRMREZGsGIT93GUJrcheexSPHVoK/HU6kP02gIF1wb15cmc5AFAqBFL1wVwRJiIiIlm51EeYxqfpUh7WR2RAI7oCbWMx8NEDAIDy2MsHve7cw3Lu7yPcLc0YgqySBo/dj4iIiMY/rgj7sdXYBY2w9j1oaQe2/Qr5led2cuuw2Pqc4umH5QBHEC6pbx8wFyIiIqLRYhD2Y+EYZNvixhLk9ypDOFbWd1c3pUJAqRAeqxEGHEFYksA6YSIiIpINg7Afa0So8zfC43G6uhXCsW8GDp+tH3CKWun5IAywcwQRERHJh0HYj23DMpj7l4mrA4HVj6KwthULk6MQFxGIw8UNA65VKxUee1gOAFL0wVAI4BSDMBEREcmEQdiP5Ygp+AhrUGsNgl0CSux6lCz7AzDzehTWtCJFH4w5iRE4XDRwRThApfDoinCASom4yECcqW3z2D2JiIhofGMQ9nM5Ygr+iu9D89s2rJWexOOFU9HUYUFtqxnJ+mDMSYxEWWMHKho7+lynVio8+rAcACTrglFY0+rRexIREdH4xSBMUKvViAgPx7pYK77Mq8Kmw6UAHMFzTtdWy0eK+64Kazy8Igw4yiMKa1ohSZ5r20ZERETjF4MwAQBMJhOMbWegC9bg0Q+PAXAEz2mxYdAoFTh8tqHP+WqlwqN9hAEgSReM5k4r6lrNHr0vERERjU8MwgTAsdVy7tEs3H3RhJ5jSbogBKiUmBobNiAIazz8sBwApOiDAACFtSyPICIiorFjECYAgMFgwLFjx/CdRUk9x7RqJQBgTmIEsksb+pRCqL1QGpGsCwYAnKnhA3NEREQ0dgzCBACIjo7GqVOnEKhR9hzr7h88NzESHRY7jpef24BDoxQef1guPjIISoXgA3NEREQkCwZhAgDodDq0tbWhuLgYmq4tlP/xVQEAYHZCBADgSElDz/neeFhOo1IgLiKQpRFEREQkCwZhAgAoFArExcXhwy2fw2yzIypYg6+OV+FoSSPiIwMRFaxBdq+NNRwPy3k2CANAsj6YQZiIiIhkwSBMPWJjY7Ftxx4AwKNXTEWYVoV/fJUPIQRmxYcjq/eKsBcelgOAFF0QCmva2EKNiIiIxoxBmHrExMQgOysLADA9LhzfW5aCrbmVyCtvwsz4CORXtaCl0wrAOw/LAY4Wai2dVtS0sIUaERERjQ2DMPVYn9yBr9YU4HTAjZjw30W4IyIToQEqPPVVAWYnRECSgJzSRgBAgFIBsxeCcIre0TmiiOURRERENEYMwgQAmC7l4dtBu5EUDigEIBpLEPT5j/C79OPYklOO4AAVACCrq05YrVTAYvV8eUKyvruFGoMwERERjQ2DMAEAVmMXNLD2PWhpx7qq56BRKvBeZgniIwORXeJYEVarhFdKI+IjAx0t1LgiTERERGPEIEwAgHA0Oz2uaCrFt2bF4oPDpYgJ1+JI14qwRqns00d4z549eOqpp1BeXu7WeaqVjhZqZ+va3XofIiIiGv8YhAkA0IhQp8dL7Tq8k1kCs82Og4X1KG1oR01LJ9QqAbPNDrvdjttuuw2XXHIJnnzySUyZMgV79uxx61wTogJRXMfd5YiIiGhsGIQJALANy2CGqs8xSR2ItuX/g0UpUX2Ol+98FfceuRJ5yo1o+t9EBJz4EHfddRduuukmrFixAtdccw3a2923YpsQGYSSegZhIiIiGhsGYQIA5Igp+Ahr0IBQ2CWgQQqFWP93TFzzfbx5x2L8dO0kAMC3FLswYd//IKyzAgoBRIhmPHkxsCSkBAAwd+5chISE4LHHHnPbXBOiglDTYkZrp3X4k4mIiIgGwSBMPXLEFDwpbsf1x9bgog/1wMzrAQBCCNy7Mg1LUnV4SPU2gkTfHr4BChtWY1fPucuWLcO///1v2Gw2t8wzPjIQAFBSzzphIiIiGj0GYRogPT0dx48fR3V1dZ/jd180AXGixuk1vR+2S0pKglKpxOuvv+6W+SVGBQEA64SJiIhoTBiEaQCtVovU1FT885//7HN8ecd2DNY5uPfDdkIIzJgxA6+88opb5pfQHYRZJ0xERERjwCBMTs2ZMwcvvvgiJOlc9BXbfgWFGHiuBMfDdr1NnToVe/fuhdks/1bIumANAtVKFLOFGhEREY0BgzA5lZaWhtbWVrz22mvnDjaWDHp+jpjS53VkZCTCw8Px8ccfyz43IYSjhdoYVoQLa1rxSbZ7ex4TERGRb2MQJqcUCgUuueQS/OQnP+mpFW5TRzk9d7AexCkpKdi0aZNb5pcQGTSmGuEH3zqCe18/hFPVLTLOioiIiM4nDMI0qEmTJmHChAmYOXMm1q5di3vfr0KnXdnnHDNUA8oiuqWmprptc42EKEcQ7l26MRLdO+QdOdsg36SIiIjovMIgTENau3YtVq9eDYVCgfDlP8DHikvQgFBIABoQio+wZkBZRLfExEQUFxejpsZ5p4mxSIgKQqvZhvo2y4ivbeo4d81Zdp4gIiLyW6rhTyF/JoRAeno60tPTAQA5iEQOnAff/jQaDWJiYvDpp5/i5ptvlnVeCV29hIvr2hAVrBnRtaW9+g9XNnXIOi8iIiI6f3BFmNwqLi4OX331lezjjqWFWnlje69fMwgTERH5KwZhcqv4+HhkZmbKPm5PEB5FC7WGrnKKFH0w6lrlb+9GRERE5wcGYXKr+Ph4nDp1Cna7XdZxQwJUiArWjKrGt7HdEYQTooJ6fk1ERET+h0GY3CosLAwajQYHDx6Ufez4yECUjKI0ojv8xkcGMggTERH5MQZhcru4uDh8/fXX8o8bETiqGt/GdgtCAlTQBWvQ1GGB3T66FmxERER0fmMQJrczGo04cOCA7OPGRgSitL59xL2Em9qtCNOqEB6ohiQBzZ1W2edGREREvo9BmNwuJiYGeXl5so8bGxGIdout5+E3VzW2WxAWqEZ4oBoA0MTyCCIiIr/EIExuFx0djbNnz8r+wFxchBYAUNowss4RTe0WhAeqEdYVhFknTERE5J8YhMntQkNDIUmS7KvCsRGOTTXKRhqEOxwrwkEax3bR7RabrPMiIiKi8wODMLmdEAKxsbHYuXOnrOPGjTIIt5qtCAlQIUjj2FixlTXCREREfolBmDzCZDLJ3kItKliDAJUCZSPsHNFutkGrViI4QNnzmoiIiPwPgzB5hMlkwtGjR2UdUwiBuIjAEdcIt5ltCNIoEaTuWhFmECYiIvJLDMLkESaTCYWFhbKPGxsROKLSCEmS0G7pCsJdK8JtZpZGEBER+SMGYfIInU6H+vp6tLS0yDpubIQWpfWuB+FOqx2SBEdpRFeNcBtXhImIiPwSgzB5hEqlQmRkJPbv3y/ruLERgahq7kSn1bUw210PHKRRQqtWQAigjQ/LERER+SUGYfIYo9GIjIwMWcfsbqFW2djp0vltlnNBWAiBILWSNcJERER+ikGYPEan0yErK0vWMeO7grCrD8y1d9UDa9WO+uCgABVLI4iIiPwUgzB5jNFoxIkTJ2Qdc6SbarSbHbvbdfcQDtIo+bAcERGRn2IQJo8xGAwoLi6WdczocMc2y64G4e7QG9i9IqxRobWTK8JERET+iEGYPKa7c0Rzc7NsY2rVSuhDAlDW6OKKcFeNcGDX9spatcLlB+2IiIhofGEQJo9RqVSIioqSvXNEXIQWpQ2u7S7Xu2sEAASoFOi02GWdDxEREZ0fGITJowwGg1s6R5TWt7l0bveDcd2lEVq1kivCREREfopBmDxKr9fLvtVybEQgyhs7IEnSsOe2WwauCHdwRZiIiMgvMQiTRxkMBtk7R0SHadFmtqHZhY0xuksjAnuCMFeEiYiI/BWDMHmU0WjE2bNnZR2zu3NERePwdcI9D8upzz0sxxVhIiIi/8QgTB4VFRWF+vp6tLa2yjbmSIJwp9UGpUJApXR89bkiTERE5L8YhMmjVCoVIiMjceDAAdnGjA4bQRC22BGgOve154owERGR/3IpCAshLhVCnBBCFAghfjbEeQuEEDYhxLXyTZHGG5PJJGvnCFNXEC53IQibbXZoegXh7hVhVx60IyIiovFl2CAshFACeBrAOgBTAdwghJg6yHl/APC53JOk8UWv1yMrK0u28TQqBfQhGlQ0uRCErXZolH1XhO0SYLUzCBMREfkbV1aEFwIokCTptCRJZgBvAtjg5Lz7AbwHoErG+dE45I7OEaYwLSpc2F2u02pHgLrvijAAdFhYJ0xERORvXAnCcQCKe70u6TrWQwgRB+AqAP8aaiAhxB1CiAwhREZ1dfVI50rjhMFgQFFRkaxjxoRrUdHUOex5/VeEu0Nxp5V1wkRERP7GlSAsnBzr//fIfwPwsCRJQy6rSZL0nCRJ8yVJmm8wGFycIo033Z0j2tpc2w3OFSNZEdZ0rQIDgJYrwkRERH7LlSBcAiCh1+t4AGX9zpkP4E0hRCGAawE8I4S4Uo4J0vijUqkQFRWFgwcPyjZmTLgW9W2WYQNtp9XWp2sEV4SJiIj8lytB+CCAdCFEihBCA2AjgM29T5AkKUWSpGRJkpIBvAvgHkmSNsk9WRo/TCaTrEG4u3NE5TAPzJmtA7tGAFwRJiIi8kfDBmFJkqwA7oOjG0QegLclSTomhLhLCHGXuydI45NOp0N2drZs48WEBwIYvpew2WbnijAREREBAFSunCRJ0hYAW/odc/pgnCRJ3x37tGi8MxgMOH78uGzjRYcHAMCwLdQ6LXbognu1T+taEe7kphpERER+hzvLkVcYDAYUFhbKNl70CFaENU5WhDu4zTIREZHfYRAmr5C7c0RIgAohAaphd5cb0D6tKxR3skaYiIjI7zAIk1d0d46Qc6vl6HDtsCvCjq4R59qndQdhs407yxEREfkbBmHyGpPJhAMHDsg2XnSYdtga4f5dIzRKRyi28GE5IiIiv8MgTF6j0+mQlZUl23iurAj3D8JqlWO/GLONQZiIiMjfMAiT1xgMBuTl5ck2XnSYFtUtnbAOEWo7rX3bp6m76oUtDMJERER+h0GYvMZgMKCoqEi28aLDtbDZJdS0mJ2+b7dLsNqlvqUR3TXCLI0gIiLyOwzC5DU6nQ4NDQ1obm6WZbzort3lyhvbnb7fXf7Qt0ZY0ec9IiIi8h8MwuQ1SqUSOp0Oe/bskWW86PDubZY7nb7fvWlG764RPaURVnaNICIi8jcMwuRVsbGx2L17tyxjmcK6g7DzB+Y6bY5ewb1XhJUKAaVCsEaYiIjIDzEIk1eZTCbZegnrgjVQK8WgLdS664ADlH2/9mqlYGkEERGRH2IQJq+KiYmRrXOEQiFgDNWicpAWap3dQVjdPwgr+LAcERGRH2IQJq+Kjo5GaWkpzGbnnR5GyhQWMOyKsKbfinCASsHSCCIiIj/EIExeFRAQgPDwcOzbt0+W8aLDtYPWCPcEYRVXhImIiIhBmHxAbGwsduzYIctYxlDt4F0jrAO7RgCOIMwVYSIiIv/DIExeFx0djf3798szVrgWLZ1WtHRaB7w32IqwRqWAxcb2aURERP6GQZi8Li4uDtnZ2bKM1b2pRoWTB+a6V33VStHnuFqp6FktJiIiIv/BIExeFxcXh4qKCtTX1495rKF6CZt7gnC/FWEl+wgTERH5IwZh8jqVSoXY2Fhs2bJlzGOd211uYBC2dpU/DAjC7BpBRETklxiEySckJCTgiy++GNE1drsdL7zwAh555BEUFRUBcLRPA+C0hdpQpRHsGkFEROR/VN6eABHgCMJ79+51+fy2tjasWLECJSUlMBqN+Oc//4kvv/wS8+fPR6hW5XRTDcsgpRFqpQKtTh6uIyIiovGNK8LkE5KSklBYWIi6urphz5UkCVdccQVaW1tx++2345prrsHixYtx3XXXwWq1IjpM63RF2GofvDTCzK4RREREfodBmHxCQEAAEhMT8eabbw577pNPPomjR4/i6quvhkrl+EuNxYsXw2q14m9/+xuiw7WocNJLuHtFWNWvNEKjVMBstcnwU7hXS6cVmUX1+CK3El/kViKntBEdFt+fNxERka9iaQT5jAkTJmDTpk245557Bj2ntLQUjz32GK655hpoNJqe40IILFu2DE8//TSu+eMaFFTVDLi2u1ewWtG/NEL4bB9hs9WOTUdK8V5mCQ4U1kHqN81AtRJL0/S4aXEiLpxogBDC+UBEREQ0AIMw+YxJkybhxRdfREdHB7RardNzbrnlFkyePBlJSUkD3ktPT8cnn3yCttLjqGoOhc0uQak4Fwx7aoRV/VaEfbRrxNZjFfjVx7koqW9HqiEY961Mw6z4CJjCtLBLEkrq23HgTC225FTgy5crMSs+HL/aMB2zEiK8PXUiIqLzAksjyGdERUXBYDDgxRdfdPr+O++8g4yMDKxatcrp+0qlErNmzULWlx/AZpdQ29K3PMLaXRoxYEXYt7pGtJmt+NFbR3DHa5kICVDh5e8uwLb/dyF+fMkkXDzVhBnx4ZiVEIHLZ8bgfzdMx+6HV+GP185EWWMHrnpmN57eXgC73TdXuImIiHwJgzD5lFmzZuG5554bcLy1tRX33Xcf1q1bN+hqMQBMnDgRxzN2ARjYQq2nNMJZ+zQfWRGuaOzA9c/uxYdHSvHD1enYfN8yrJxsHLLkQaNS4Pr5Cfjqxxfiipmx+NPnJ3DnfzJZP0xERDQMBmHyKdOnT0dhYSE+/vjjPsdvvvlm6PV6TJkyZcjr4+Pj0dHSCEtDxYBtli02O1QKMSBUBvhIaURZQzuu/dceFNa04YVb5+NHayZCo3L9X9FQrRpPbpyNx9ZPxZd5lfjuywfQwrZwREREg2IQJp+iVquxatUq3HnnnT1bLj/66KPYvn07rrjiimGvVygUSJs4CR2nM1DZ3K80wi4N6BgB+EZpRGVTB258fh8a2y14/fZFWDXZNKpxhBC4bWkK/nr9bBwsrMcdr2ag8zzoiEFEROQNfFiOfM6sWbNQVlaGiRMnwmQyoaysDDfffDOCgoJcuj4lKQGn8o4O2FTDbLUP6CEMOIKwXcKAh+s8pc1sxW0vH0R1cyde+8EizIyPGPOYV86Jg12S8P/ezsJP3snGk9+eDYUXfjYiIiJfxiBMPkcIgXXr1uHs2bNobW3Fhg0b+rRKG87VaVY8mXwUiXsXAbnxwOpHgZnXw2p3HoS7yw/MVjsCNUrZfg5XSJKEn76TjbyKJrx06wLMTYyUbeyr58ajsqkTf/jsOKbEhOKei9JkG5uIiGg8YBAmnySEcNoibTjTpTysD94DR2yWgMZi4KMHAAAW6+QBD8oB5x6es9jtCIRng/DzO0/jk6Pl+Nm6yVg52Sj7+HddmIrc8ib83+cnMCchEksm6GS/BxER0fmKNcI0rqzGLmjQ7wExSzuw7Vew2O0DWqcBgKqrZMDq4U01csua8KfPT+CSqSbcuSLVLfcQQuB3V89Aij4Y979xGHWtZrfch4iI6HzEIEzjSjianb/RWAKrTXLahUHVVS5h9WDniA6LDQ++dRgRQRr8/pqZbt0RLiRAhae/MxeN7WY8+mGO2+5DRER0vmEQpnGlEaHO3wiP72mf1t+50gjPrQj/7ct8nKxswZ+unYmoYNfrn0drcnQYfrg6HR9nl+PTo+Vuvx8REdH5gEGYxpVtWAZz/9J3dSCw+lFYbFLP6m9v3eUSnloRzq9sxgs7T+PaefG4aJL8dcGDuevCCZgRF45ffpiDpg6Lx+5LRETkqxiEaVzJEVPwEdagAaGwS0CxJRzNa/4CzLweFpsdGicPy3X3FrZ4oEZYkiT8YlMOggNUeGTdZLffrzeVUoHfXT0Dta1m/P3LfI/em4iIyBcxCNO4kyOm4ElxO+ZsmYZZ+y/FJtsFAACr3T70irDd/SvCHx4pw/4zdXj40snQhQS4/X79TY8Lx8YFifj3nkIUVA1ST01EROQnGIRp3EqONUBRnY/NWWUAHCu+ztqnda8Iu7trRIfFhj9+dhzT48KwcUGCW+81lJ9cMhFBGiX+96Ncr82BiIjIFzAI07gVGxsLS9kJHCysR2lDOyy2wXaW6wrCbn5Y7j/7ilDW2IFH1k3x6i5vupAA/PDiidiZX4PdBTVemwcREZG3MQjTuKXX69HWUA3JasZHWWWw2iSnXSM88bBcU4cFT28vwPJ0PZam6d12H1fdtDgRseFa/OnzE5Akz/ZPJiIi8hUMwjRuqVQq6HQ6xEvV2HykbNAVYU88LPf8jtOob7Pg4Us9+4DcYAJUSjywOh1HihvwZV6Vt6dDRETkFQzCNK6ZTCYY2wqRW96E4xXNg5RGuPdhucY2C17adQaXz4jB9Lhwt9xjNK6ZF49kXRD+vPUE7B7soUxEROQrGIRpXDMYDOgszUP3xm1OH5Zz8xbLr+wtRKvZhvtWpbll/NFSKxX44cXpOF7RjG3HuSpMRET+h0GYxjWTyYSiUyexJFUHAE7bp3WvCFvcUCPcZrbi5d1nsHqyEVNiwmQff6zWz4xFXEQg/vXNKW9PhYiIyOMYhGlcMxqNKC4uxuUzYgAAp6pbBpyjcmPXiDcOFKO+zYJ7VvrWanA3lVKB25enILOoHgcL67w9HSIiIo9iEKZxLTw8HBaLBVNCzQCAw2cbBpzT3TVC7hVhs9WO53ecxuLUKMxLipR1bDl9e0EiooI1+NfXXBUmIiL/wiBM45oQAiaTCXlHDgx6jtpNG2p8mlOOiqYO3Lligqzjyi1Qo8StS5Kx7XgVTlZytzkiIvIfDMI07hmNRmRmZva8zu8X9lRu6hrxyp5CpOiDceFEg6zjusMtS5IQoFLg1b2F3p4KERGRxzAI07hnMBiQk5PT8/qTo+V93lcr5O8jfLSkEYfONuDmxUle3UXOVZHBGqyfFYv3D5WiqcPi7ekQERF5BIMwjXtGoxEFp87Vv356tKLP+8qe9mnyrQj/e08hgjRKXDs/XrYx3e3WJcloM9vwfmaJt6dCRETkEQzCNO4ZDAaUl5VBkiQYQwNworIZBVXnukecK42QZ0W4tqUTH2WX4Zq58QjTqmUZ0xNmxIdjdkIEXt1XxG2XiYjILzAI07gXHBwMpUoFW3Mt1s+KBQB8lnOuPEItc/u0tzKKYbbacesFSbKM50m3XpCE09Wt2FVQ4+2pEBERuR2DMPkFozEalpoixEcGYn5SJD7pVR7R3T5NjtIISZLw9sFiLEyJQpoxdMzjedplM2IQFazBmweKvT0VIiIit2MQJr+gMxphqSmCSqnAuhkxyCtvwpmaVgDnVoTleFjuYGE9CmvbcP38hDGP5Q0BKiWunB2HL3IrUd9q9vZ0iIiI3IpBmPxCpM4Ac81ZqBUC66ZHA3D0+QUcvYaVCiFL+7S3DhYjJECFy2ZEj3ksb7lufjzMNjs+PFLq7akQERG5FYMw+YX1yR04tiYD394yE7EvL8ADhsPY0quNmkohxryhRnOHBVuOlmP9rBgEaVRjnbLXTIkJw4y4cLydwe4RREQ0vjEI07g3XcrDHbqDSAq1Q0ACGotxf+s/kFq+BcV1bQAAtVIx5tKIT7LL0W6xnbdlEb1dNz8eueVNyClt9PZUiIiI3IZBmMa91diFAGHtc0xt78BDqrfxRW4lAEClHHtpxNsZxUg3hmB2QsSYxvEF35oVC41KgXfZU5iIiMYxBmEa98LR7PR4rKIWW3Md3SNUirGtCBfVtuLQ2QZcMy8eQvj+TnLDiQjS4JKpJmw6UgqzVd6tp4mIiHwFgzCNe41w3sasOcCEA2fqUN9qhlopxtQ+7aOsMgDo6VM8Hlw1Jw4NbRbsKqj29lSIiIjcgkGYxr1tWAYz+j68JqkD0bjkEdglYNvxqq7SiNGtCEuShA+PlGFBciTiIgLlmLJPWJ5uQHigGpuPlHl7KkRERG7BIEzjXo6Ygo+wBjXWINgloMSuR/7CJ5Bw4a2ICddi67EKqBUKWEa5Iny8ohn5VS341uw4mWfuXRqVApfNiMbW3Eq0m23eng4REZHsGITJL+SIKfiL/TZof9eODaqn8buSGRBCYM1UE3bkV8Nit4+6fdrmrDIoFQKXTT9/ewcPZv3MWLSZbdh2vNLbUyEiIpIdgzD5jYCAAAQHB2NtArD9RDVOVjbjkqnR6LDYUVzXPqrSCEmSsPlIGZal6aELCXDDrL1rUaoOxtCAnhpoIiKi8YRBmPxKdHQ0YjqLoVUr8PLuM1iUGoVQraN+eDTt0w6drUdpQzs2zB4/D8n1plQIXD4zBttPVKOpw+Lt6RAREcnKpe2vhBCXAngSgBLAC5Ik/b7f+98B8HDXyxYAd0uSlCXnRInkYDQacSInC1etWoT3D5XiobWTsXqyEZuOlKHTMvIg/FFWOTQqBdZMNblhtr7hW7Ni8fLuQmw9Volr58V7ezqysdrs+DKvCtvyKnHobD3KGztgtUmICFJjSkwYlqXpsWF2LIxhWm9PlYiI3GTYFWEhhBLA0wDWAZgK4AYhxNR+p50BcKEkSTMB/BrAc3JPlEgOOp0Ox44dw3cvSEGn1Y7XD5zFJdMctb37ztSOaCxJkrD1WAVWpBsQqlW7Y7o+YXZCBGLDtfj8WIW3pyILu13C2weLsewP23HXfzLxRV4lUg0h2LggEbctS8aKiQZUNHbgiS15WPL7r/DI+9mobOrw9rSJiMgNXFkRXgigQJKk0wAghHgTwAYAud0nSJK0p9f5+wCMn2UjGleMRiMOHDiASdGhWJamx2t7i/D5gysAANIIS4SPljairLED/++SSW6Yqe8QQuCSadF448BZtJmtCNK49BdJPqm4rg0PvHkYh882YE5iBH595XSsmmyEUjFwE5RT1S14ZU8h3jhwFh9nl+Px9dNw9dy4cbFhChERObhSIxwHoLjX65KuY4P5PoBPnb0hhLhDCJEhhMiormaTfvI8vV6PyspKWCwW3LY0GRVNHdiRf+67KI0gDX9+rAJKhcDFU4zumKpPWTstGp1WO745cf7+e7v9eBUu+/tOFFS14K/fnoX3774Aa6aanIZgAJhgCMGvNkzHFz+6EJOjQ/Hjd7Lwyw9zxrTxChER+RZXgrCz/0o4TQtCiJVwBOGHnb0vSdJzkiTNlyRpvsFgcH2WRDJRq9UIDw9Hbm4uVk4yIlkXhJd2n+l5/3iF8+2YnfkspwKLU6MQEaRxx1R9yoLkSEQGqfHZeVoeselwKX7wagYSo4Kw5YHluGqO61thJ+uD8eYdS3Dnhan4z76z+MGrGeiwsK8yEdF44EoQLgGQ0Ot1PIABvZSEEDMBvABggyRJIyu2JPIgk8mEffv2QaEQuPWCZBw+2wBV16rgV8erXBqjoKoZp6pbsXba+Osd7IxK6Xgg8Ku8Kpit59eK6AeHS/DgW0ewMDkKb925BAlRQSMeQ6kQeGTdFPz2qhn4+kQ17v5PJjqtDMNEROc7V4LwQQDpQogUIYQGwEYAm3ufIIRIBPA+gJslSTop/zSJ5GM0GnHo0CEAwHXzExAaoOrpIbwtz7WNIz4/5jjvkqn+EYQBR3lEc6cVe07VeHsqLttxsho/fScbS1J1ePm2BQgJGFt9842LEvHbq2Zg+4lq/OSd7BGV0hARke8ZNghLkmQFcB+AzwHkAXhbkqRjQoi7hBB3dZ32KAAdgGeEEEeEEBlumzHRGHV3jgCAkAAVrpt/7i88Dhc3oLalc9gxPj9WgdkJEYgO95/WWkvT9AjWKHv+EODrTlY24+7/ZCLdFIpnb5kHrVopy7g3LkrEw5dOxkdZZXh6e4EsYxIRkXe4tKGGJElbJEmaKEnSBEmSnug69i9Jkv7V9esfSJIUKUnS7K5/5rtz0kRjYTQacerUqZ7XtyxJ6vm1JAFfD/NAWFlDO7JLGv2mLKKbVq3ERZON+CK3AvZR7MLnSa2dVtzz30MI1Kjw8ncXIEzm9nZ3XZiKq+bE4f+2nsSXuefHHwyIiGgg7ixHfken06GmpgadnY6V32R9cM97+hDNsHXC20843l8zdfx3i+hvzRQTalrMyC5t9PZUBiVJEn7+wVGcrm7B32+Y7ZZVeyEEfnf1DEyLDcNP381CFfsMExGdlxiEye+oVCpERkbi8OHDPcf0Id2dHwR2nKyGZYgWWduPVyEhKhATDCFunqnvuXCiAQrh+kOF3vDJ0XJsOlKGBy+eiAsm6N12H61aiSc3zkG7xYYfv5Pl86vkREQ0EIMw+aXo6Gjs3bu353V3mUNNSyeaO604WFjn9LoOiw27C2qxcpLRLzdWiAzWYE5iJLb7aBCubenEYx8ew6z4cNxz0QS33y/NGIJfXD4VO/Nr8PqBs26/HxERyYtBmPySXq/vsyIc2O9Bqq/ynAe9fadr0W6xYeVk/yuL6LZqshFHSxt9shzg8Y9y0dRhwR+vnQWV0jO/vX1nUSIumKDDHz47jurm4R+0JCIi38EgTH5Jr9cjN7dnl3AolY7V3e5Nxgb7q/+vT1RDq1ZgSarO7XP0VSsnOf4QMNxDhZ62M78aH2WV4b6V6ZgUHeqx+woh8Osrp6PTYsdvPskd/gIiIvIZDMLkl4xGIwoLC3teqxQCKoXoCXmna1pxpqa1zzWSJOGr41VYOkEvWyuu89GUmFBEh2l9qk7YarPjVx/lIjEqCHddlOrx+08whOCuiybgwyNl2FNw/vRZJiLydwzC5Jd0Oh1aWlpQXl4OAFAqFLDaJdy4KLHnnP5B71R1K87Wtfl1WQTgWAFdOdmAXQU1PrPL3H/3n0V+VQt+fvkUBKi884eUey6agPjIQDyxJY8PzhERnScYhMkvKRQKxMXFYevWrQDQs8XysnQ94iICAQBfHe/bH/brrrZp/h6EAUd5REunFRmDPFToSY1tFvzli5NYmqbDJVNNXpuHVq3ET9dOwrGyJmzOGrALPRER+SAGYfJbcXFx+OabbwAAqq4aYUkCvr3AsdPc7oJaNHdYes7/6ngVJplCe4KyP1uapodGqfCJ8ojnd55GY7sFP79sqtc7eayfGYvpcWH40+cn0GGxeXUuREQ0PAZh8ltxcXHIyHDsBt69ImyzS7hmXnzPOTtOOuo9W7taql002eD5ifqg4AAVFqREYme+d+th61rNeHn3GVw+IwZTY8O8OhcAUCgE/mfdFJQ2tOO/+9lOjYjI1zEIk9+Ki4tDQUEBbDYblArHvwpWu4S4iMCerhDdK54HztTBYpOwIp1BuNuyNANOVDZ7tY3asztOoc1iw4MXp3ttDv1dkKbHklQdnv3mFFeFiYh8HIMw+a3Q0FAEBATg4MGDfVaEAWDjQkd5xHuHSiBJEnbm1yBApcC8pEivzdfXLE937Nq2y0tdEqqbO/HqniJsmBWLdJPn2qW54v7Vaahq7sQ7GcXengoREQ2BQZj8WnJyMjZv3gxlVxC2dm2t3L3THADkljdhV0E1FqZEQatWora2Fj/+8Y9x5513oqCgwCvz9gVTY8IQFazBLi+VR7y0+ww6rTY8sNp3VoO7LUnVYV5SJP71zWmf6axBREQDMQiTX0tOTsa2bdug7npYztq1IqxVK7GmqwPBi7vO4GRlC5al6ZGfn49p06bhs88+Q2ZmJhYsWICjR496bf7epFAIXDBBh10FNZAkz7YLa+m04j/7inDp9GikGkI8em9XCCFw/6o0lDa044PDJd6eDhERDYJBmPxaamqqI8jaHat2tl79X+9bmQYAeP9QKQBgbmwQ1qxZg8mTJ+P666/H+vXrMXfuXHz729/2eBD0FcvT9ahq7sTJyhaP3vfNA2fR3GHFHSsmePS+I3HhRANmxIXj2R2n2VeYiMhHMQiTXwsLC0NoaChyMvcCOLciDAAz48N7fq1WCvzt8Z8gMDAQF110Uc/xpUuXorKyEq+//rrH5uxLlnU9PLgz33PbLVtsdry06wwWpkRhdkKEx+47UkII/GB5Ck5Xt+IbD34+RETkOgZh8nvp6enY9flmAIDNfq6eUwiBnyccxS7NAzihugG/M36En39rYp9etUqlEsuWLcMf/vAHj8/bF8RFBCLVEOzRNmpbjpajrLEDdyz3/FbKI7VuegyMoQF4eXeht6dCREROMAiT3/vuvGC8MmELTgfciMRXFwHZbzveyH4b36v7K+IVNVAIIDFc4PrAXZgu5fW5fvr06cjPz8fJkye9MHvvW56mx/4ztei0eqZV2Eu7C5FqCMaq82CHP41KgVuWJGHHyWrkVzZ7ezpERNQPgzD5telSHm6JzEBiGKAQgKalFPjoAUcY/vRhKG19e+RqYMVq7Op7TKPBlClT8NRTT3ly6j5jWboBHRY7Movq3X6vnNJGZBU34ObFSVAovLuLnKtuWJiIAJUCL+8p9PZUXCZJEqqbO3HobD0OFtahoKqZ3S+IaFxSeXsCRN60GruggbXvQUs7zB//FGpzA5xFrXAMXNmbNGkSPv/8c/dM0sctTo2CUiGw91QtLpigd+u9/rOvCIFqJa6eGz/8yT5CFxKAK2fH4f1DJfjZuskI06q9PaVBnapuwWt7i7D1WAXKGvv9IVCpwJIJOlw1Jw6Xz4yBWsl1FCI6/zEIk19zFmoBQN3ZADHIgmMjBm7ekJqaivfffx+VlZUwmUxyTtHnhWrVmB4Xjn2na916n6YOCz48UoZvzYpFeKDvhklnvrM4EW9lFOPDI2W4eXGSt6czQHVzJ574JBcfZpVBrVRgRboB31+eihR9ENRKBWpaOnG0pAlf5lXiwbeO4P+2nsAvr5jap982EdH5iEGY/FojQhHhLAwPEoIlANuwbMBxjUaDpKQkvPHGG3jwwQdlneP5YHFKFF7afQbtZhsCNUq33OP9zBK0W2y4yQeD5HBmxIVjakwY3jxw1ueC8Gc5FXj4vWy0m22468IJ+MGyFOhCAgacd9Uc4BeXT8H2E1X40+cncOdrmbhiZgz+cM1MBAfwPyVEdH7i322RX9uGZTD3+/NgmwWoa3d+fhu0yBFTnL6XkpLix+UROlhsEg6fdU+dsCRJ+O/+s5gVH44ZvdranS+EELhhYQKOlTXhaEmjt6cDwPGZ/vWLk7jrP5lI1gVhyw+X4eFLJzsNwd0UCoHVU0z46P5l+OnaSdhytBxXPbMbxXVtHpw5EZF8GITJr+WIKfgIa9CAUEgAGhCKN9suwAcdiwcEZDNU+AwrBx0rMTERWVlZbp6xb5qfHAmFgNvKIw6dbUB+VQtuXJTolvE9YcOcOGjVCrxx8Ky3pwJJkvDY5mN4cls+rp4bh7fuXII048CSn8GolQrcuzINr31/ESqbOnH9s3txutqzm6oQEcmBQZj8Xo6YgifF7fiV+H94UtyO4ojFKI26YEBA/ghrBl0NBoCYmBjU1taipMT/ttQN1aoxIy4c+07XuWX89w+VQKtW4LIZMW4Z3xPCtGpcPiMWm4+UobXTOvwFbiJJEv73o1y8urcId6xIxZ+vmwWtenTlLEvT9Hjj9sUwW+349nP7uDJMROcdBmGiQfQPyEOFYMCxuUZ8fDw++eQTD83QtyxO1eFIcQPazfL2E+6w2PBRVhkunRaNUB/uuOCKGxYmoKXTii1Hy702hxd3ncG/9xTiB8tS8Mi6yX02iBmNqbFheOOOxei02HDrywfQ0GaWaaZERO7HIEwko4SEBHz11VfenoZXLE7VwWyzy14nvC2vCk0dVlwz7/xpmTaYeUmRSNIFYdORUq/c/8vcSjyxJQ/rpkfjfy6bMuYQ3G2iKRTP3zIfJXXtuOe/h2DrtVU5EZEvYxAmklF8fDzrhGWuE37vUAmiw7Ru71HsCUIIXDk7DntO1aKiX59edyupb8OP3j6C6bHh+Mv1s2XfkGRRqg6/uWo69pyqxZNf+ucui0R0/mEQJpJRTEwMioqKYLf73y5cPXXCZ+SrE65u7sQ3J6tx5Zw4KM+TneSGc+WcOEgS8FFWmcfuabXZ8cM3jwAS8Mx35rqtxd318xNwzdx4/GN7AXbmV7vlHkREcmIQJpJRSEgINBoNMjMzvT0Vr1icqsORsw3osMhTJ/zhkVLY7BKumRsny3i+IEUfjFkJEfjgsOfKI57aXoDMonr85qrpSIgKcuu9fn3lNKTqg/HQu9lo7rC49V5ERGPFIEwks7i4OHzzzTfenoZXdNcJHyqSp054c1YZZsSFI93kemuv88FVs2ORW96EExXOdzaUU35lM57eXoBvzYrFhtnu/wNFkEaF/7tuFiqbOvC7T4+7/X5ERGPBIEwkM5PJhAMHDnh7Gl4xLzkSQgAHCsdeHnG2tg3ZJY1YP+v8bZk2mCtmxUKpEG5/aM5ul/A/HxxFcIAKj66f6tZ79TYnMRLfX5aC1/efxd5T7t16m4hoLBiEiWQWGxuLnJwcb0/DK8K0akyODkOmDCvCn3S1GDufewcPRh8SgOXpenyUVQZJcl+HhbcyinGwsB4/v2wK9EPsGOcOP75kEhKiAvH45mOw2vyvZp6Izg8MwkQy8+cH5gBgflIkDhXVjzn8bDlajtkJEYiPdG9Nq7dcNiMGJfXtyCltcsv4TR0W/OnzE1iYEoVrvdB6TqtW4ueXTcWJyma8fsD7u+kRETnDIEwks9DQUAghkJ+f7+2peMX85Ei0mm04Pob616LaVhwtbcQVM8ffanC3S6aaoFIIbMlxz+Yaz2w/hfo2Mx69Yqps/YJHau00Ey6YoMOft55EfSs32iAi38MgTOQGJpMJe/bs8fY0vGJeUiQAjKk8orssYt04LIvoFhGkwZIJOnx6tFz28ojiuja8tPsMrpoTh+lx4bKOPRJCCDy2fhqaOyz4x1cFXpsHEdFgGISJ3MBgMPhtC7W4iEDEhGuRMZYgnF2OuYkRiIsIlHFmvmfd9BgU1rYhr1ze7hF/+eIkBICfXDJJ1nFHY1J0KK6ZG4//7C/y+CYizjS0mfFFbiVe3HUGT36Zjxd3ncHXJ6rQ2M5Wb0T+SOXtCRCNR0aj0W8fmBNCYF5SJDJG2TmisKYVx8qa8MsrPNflwFsumWbCLzYdxac55ZgaGybLmKeqW/DhkVLcvjwVsT7yB4kHVqdj05FSPLU9H7+5cobH7y9JEr46XoV/7ynEroIaOFuAVykEVk024ralKVgyQefxORKRdzAIE7mB0Wj0262WAccDcx9nl6O0oX3Eq7qfH6sAAFw6PdodU/Mp+pAALErR4ZOj5fh/aybKUsv79FcF0KgUuH1FqgwzlEdCVBC+vSABbx0sxp0rJrh9U4/e8iub8eiHx7D3dC1iwrW4f2Ualk80IN0YgpAAFZo6rDhe0YSvT1Tj/UOl2Jq7DxdNMuCJq2aM+7+RICKWRhC5hcFgQEVFBSwW//zr1vnJUQAwqlXhL3IrMS02zG9CyGUzonG6uhX5VS1jHquwphWbjpTipkVJHm+XNpz7VqZDCIFnvvZMrbAkSXg7oxjrn9qFvIom/HrDNOx8aCX+3yWTsCA5ChFBGqiUCkQFa3DBBD3+57Ip2PXwSvzPZZNx8EwdLv3rDnx61D0PMhKR72AQJnIDrVaLoKAgHD582NtT8YrJ0aEI1iiRUTiyOuHalk5knq3HmqkmN83M96yd5lj5/jynYsxjPbW9AGqlAndc6Durwd2iw7W4bl483jtUiqpm99YKS5KE327Jw0PvZmNuYiS2/mgFbl6SDJVy6P/kadVK3LFiAj794QqkmUJw938P4ZmvC9za65mIvItBmMhN/HmHOZVSgTmJkSN+YG7b8SpIEnDxFP8JwsYwLWYlRODL41VjGqe8sR2bDpfihoWJMIZqZZqdvG5fngqLzY5/7y502z3sdgk/e+8ont95BrcsScJr31804s8jUReEN25fjPWzYvHHz07gz1tPumm2RORtDMJEbmIwGHDo0CFvT8Nr5iVF4kRFE5o6XC8P+SK3ErHhWkyT6cGx88XFk43IKm4Y00rpK3uKYJckfH9Ziowzk1eyPhjrpkfjtX1FaOm0uuUev92Sh7cyinH/qjT877emQakYXd21Vq3Ek9+ejRsWJuCp7QX4xzb/7AtONN4xCBO5icFgQG5urren4TULkqNgl4DDZxtcOr/DYsPO/GpcPNXktQ0gvGV11wr49lGuCreZrXjjwFmsnRbt0QfRRuPOFRPQ3GHFm27Ybe6Fnafxwq4z+O4FybI8fKhQCDxx5QxcPScOf/7iJD48UirTTInIVzAIE7mJ0WhEUVGRt6fhNbMTIyAEcPisa+URu/Jr0GGx+1VZRLcpMaGIDdfiy7zRBeH3DpWisd3i06vB3WYlRGBRShRe3l0Im12+2tvdBTX47ZY8XDYjWtbd9BQKgd9fMxMLkiPx0LvZOFrSKMu4ROQbGISJ3ESv16OmpgadnZ3enopXhASoMNEY6vKK8Jd5lQgJUGFxqv/1cBVCYPUUU9cfBmwjutZul/DyrjOYFR/es6ufr7ttaTJKG9rx1RjrortVNHbgh28eRqohBH+6dhYUoyyHGIxGpcA/b5oHfUgA7vpP5ojKfYjItzEIE7mJWq1GeHg4MjIyvD0Vr5mdEIGskoZhn7q32yV8mVeFCycZoFH5529Lq6YY0W6xYe/p2hFd983JapyuacX3lqWcNyUlF08xITpMi1f3Fo55LLtdwo/eOoI2sw3/umkuggPc0x5fHxKAf9w4BxVNHXh0k39ulkM0Hvnnf3GIPMSfO0cAjvKIhjYLCmvbhjzvaGkjalo6cfEUo4dm5nuWpOoQpFFiW17liK777/6z0IcE4LIZMW6amfxUSgVuXJSInfk1OF09tv7J/z1wFntP1+KXV0xFmjFUphk6NzcxEvevSsOmI2XYnFXm1nv119hmQUl9G6qaO9jOjUhG3FmOyI30ej2ys7O9PQ2vmZ0QAQA4UlyPFH3woOd9c7IaQgAr0g0empnv0aqVWJamx7a8Kvx6g+TS6m5FYwe+Ol6JOy+cAPUwPXJ9zcaFCfj7tnz8d//ZUW+nXVLfht9vycPydD02LkiQeYbO3bcyDV+fqMbjm49hRboeEUEat9zHbndsC705qwx7T9eiuvlciVWgWon5yZFYOy0aV8+NQ5CG/yknGq3z63dOovOMwWBAXl6et6fhNRNNoQjSKHFkmDrhb05WY0ZcOHQ+thuap62eYkR5YwdOVrq2SvpORjHsEjwWAuVkDNXi0unReCejeMR10YBj04yff+AoUfjd1TM8VhaiUirw26tmoLHdgj98dtwt99h+vAqXPrkDP3g1Azvzq7E8TY+fXzYFf7x2Jv73W9Nw/fx4lDa04xebcrDkd1/h+R2nYbHZ3TIXovGOf4wkciODwYBdu3Z5expeo1QIzIgLx5HihkHPaWyz4PDZety7Ms1zE/NRKyY6VsR3nKzGpOih/5rfZpfw5sFiLE3TIUk3+Gq7L7txUSI+zi7H58cqsGF23Iiu3ZZXhW9OVuPRK6YiPtKzLeOmxobhe0uT8fzOM7h2XjzmJUXJMm5LpxWPbsrB+4dLkWoIxpMbZ+OyGTGDrvZnFtXh79sK8MSWPLx/uBTPfGfukH/zQkQDcUWYyI30ej1qa2v9tnMEAMxJjERuedOgq367Cmpgl4CLJvlvWUS3mPBApBtDsCO/ethzd+ZXo7ShHRsXJHpgZu6xOEWH+MhAvJtZMqLrOq02/OaTXKQZQ3DzkiQ3zW5oD148EbHhWjz64THYZWgDV1zXhmue2YMPs8pw/6o0fPrD5dgwO27Ikpd5SVH4920L8K+b5qG8sR3f+scufHV8ZDXmRP6OQZjIjVQqFSIiIpCZmentqXjN7IQIWGwSjpU1OX3/m5NVCNOqMCs+wrMT81ErJhqw/0wd2s1Dlwu8eaAYUcEaXDLt/O27rFAIXDM3HrsKalDa0O7ydf/eXYjC2jb88oqpXquNDg5Q4aFLJ+NYWRM+yh7bg3MFVS246pk9KG9sx79vW4AfXzIJASqlS9cKIXDp9Gh88sByJOuDcfurmdz4g2gEGISJ3MzfO0fMSYwAAKflEZIk4ZuT1VieboDqPHvYy12Wp+thttqx/8zgbdQa2szYdrwSV86Oczkw+apr58VDkoD3XVwVbmyz4KntBVg12YgLJ3r3bxG+NSsW02LD8KfPT6DTOvI6ZwA4U9OKG5/fB0DCe3dfgOWjfGA0LiIQb9yxGPOTIvHgW0fwkYe7WhCdr/hfHiI3MxgMyMrK8vY0vMYUpkVMuNZpED5R2YzKpk6vBxpfsihFB41KgR0nawY955Oj5bDYJFw9d2R1tb4oISoIS1J1ePdQiUttwV7YdRrNHVb8dO0kD8xuaAqFwCPrpqCkvh2v7R35LpK1LZ24+cX9sNolvH77YqSbxtb+LSRAhVe+txALkqLw47ezsPfUyHpSE/kjPixH5GZ6vd6vO0cAjvKII8UDt1r+5oSjFnYFg3CPQI0Si1KisHOIOuFNh0uRZgzBtNgwD87Mfa6dF48fv5OFg4X1WJgy+INnda1mvLTrDC6fEYMpMb7xsy9L12N5uh7//PoUvrMoCYEa11bozVY77v7vIVQ3d+LtO5dg4hhDcDetWonnb5mPa/+1B3e+loFPHliOhKixP0xot0vYfaoG35yoRnZpI2qaOyGEY6ORmfHhWDnJiMWpOtl39SNyN64IE7mZwWBAYWGht6fhVbMTIlBc147alr4PDX5zshqTo0MRHa710sx804p0A/KrWlDmpG62uK4NBwvrcdWcuPNmJ7nhrJsRjSCNEh8cHrq29dkdp9BmseHBi9M9NDPX/HB1OmpbzXj9wFmXr3nik1wcOFOHP147E7O6+m3LJTxIjZe+uwAAcPd/M0fVnq6bzS7h7YxirPjTdtz84gG8uq8INruEKTFhmBwdBovNjlf2FOHGF/Zj5Z+/xnuZJbDJ8PAgkacwCBO5mU6nQ21tLTo6Orw9Fa85t7FGQ8+xdrMNGYX1WJ6u986kfFj3CrmzVeHuB6G+NSvWo3NypyCNCmummvBpTvmg/XBrWzrx6p4ibJgVO+YSArnNT47CklQdnttxyqXQuf14FV7ZW4TvL0sZcds4VyVEBeHP189GTmkTfvNJ7qjGOFXdgiuf3o2H3s2GLiQA/7hhDrIfuwTv3X0Bnv7OXDz9nbl4/56lyH78Ejy5cTZCtSr8+J0sXPuvPSisaZX5JyJyDwZhIjdTqVSIjIzE4cOHvT0Vr5kRHw6lQuBwr401MovqYbbZcUEag3B/E00hiA7TDqgTliQJHxwuxYLkSFn+utuXrJ8Zi4Y2C3blO6+NfnVvEdotNty3yjf7Td+/Kg2VTZ3DtoKraenET9/NwuToUDx0qXvrnNdMNeEHy1Lwn31nsePk8C35evs4uwzr/7ELJfVteHLjbGy65wKsnxULrXpg6YdWrcSG2XHYfO8y/OX6WThV1YLL/r4TX+aylRv5PgZhIg8wmUzYv3+/t6fhNUEaFSaaQpFV0tBzbPepGqgUAguT5dmMYDwRQmBZuh67T9X06VF7rKwJp6pbceWc8/8huf5WTDQgTKty2u2g3WzDq3sLcfEUI9KMvrUa3G3JBB3mJkbgn1+fgnWIXd5+8UEOmjqs+NvG2R7p+PGTtZMwwRCMR94/iuYOi0vXvLz7DO57/TCmxoRhS1c/Y1fKcBQKgavnxuPzH61AmjEEd7yWgdf2jfwhwv5aOq3IK2/C0ZJGFNe1sfSCZMWH5Yg8wN87RwDArPhwfHasApIkQQiBPadqMTshAsEB/G3ImQsm6PBuZgnyKpowLTYcALA5qwwqhcDlM2K8PDv5aVQKrJseg0+OlqPDYuuz8vhuZjHq2yy4Y8UEL85waEII3H1RGm5/NQOfH6vE5TMH/n/01fFKfHasAj9dOwmToz3zsJ9WrcSfrpuFa/+5B3/47Dh+c+WMIc9/Yedp/OaTPKydZsKTG+c4XQEeTkx4IN68YzEeeOMwfrkpBwLATYtHtvFJaUM73th/Fp8dq0BBVd8txwPVSixL1+OqOXFYOy0aSj6gR2PAFWEiD9Dr9cjNHV2d3ngxMz4CDW0WnK1rQ1OHBUdLGnDBBJ23p+WzlnR9Nt0tsCRJwqc55ViapkdEkMabU3Ob9bNi0dJpxdcnqnqO2ewSXth1BrMTIrAgOdKLsxveqslGJEYF4aXdZwa812Gx4bHNx5BmDMHty1M9Oq+5iZG4ZUky/rv/LHJKGwc97/1DJfjNJ3m4fEYMnr5x7qhCcLcgjQrPfGceVk024hebclze5KOx3YJfbDqKC/+4Hf/85hSMoQH46dpJePrGuXj+lvn4/dUzcM28OOSUNuKe/x7C6j9/zd30aEwYhIk8gJ0jgFkJjlXNI8UN2H+6DnYJWDKB9cGDiQkPRKo+GHu6gvCxsiYU17XjshnRXp6Z+yxOjYI+RIPNvcojvsitRFFtG+5ckerzXTKUCoHvXpCMzKJ6ZPXrm/3M9gIU17Xj1xumQ6Py/H96f7RmIiKDNPjfj4457de873QtHno3GxdM0OEv354lywY3GpUCz3xnLhamROGn72Yju1dplDO78muw5i/f4PX9Z3HDwkTseGglXr99Me5dmYbLZ8ZgzVQTNi5MxG+unIFdD6/Cv26aC5VSge/9OwMPvnkYbWbriOZntdlRXNeGzKJ6ZJc0oKrZfx9o9mf8O0kiD9DpdKirq0N7ezsCAwO9PR2vmGgKRYBKgeySRtglCQEqBeYmRXh7Wj5tyQQdPjxSBqvNjk9zyqFUCKyZOn6DsErpKI94J7MYbWYrgjQq/GdfEWLDtbhk2vnxc183Px5/+eIkXt59Bn/bOAcAUNHYgWd3nMa3ZsX2rPR7WnigGg+tnYSfvX8Um7PK+nSrqGrqwH2vH0aiLgjP3jxP1tplrVqJf35nLr711G7c8WomNt+/FMbQvu0SJUnCy7sL8ZtPcpFmDMGLty7AjPjwIcdVKgQunR6DVZNNeObrAjy5LR/HK5rx/C3zh3yQ1G6X8E1+Nd4+WIxd+TVo7uwbnk1hAbhkajS+vSAB0+OGngOND1wRJvIApVKJqKgoZGZmensqXqNWKjA9LhzZJQ3Ye6oWC5Kjzvvtgd3tggl6tHRakV3aiE+PVmBxahSigsdnWUS3ddOj0WGxY8fJapyubsGughrcsDDxvKkDDdWqcd38eHxytByVTY4Vxr99eRJ2SfL6bnjXzU/A9Lgw/PGzc1tC2+wS7nv9MFo7rfjXTfMQqlXLfl9dSACev2U+GtrN+Mk72X1WpCVJwl+/OIlffZyLNVNN+OCepcOG4N40KgUevHgiXrltIcoa2vHtZ/fizCCt2zIK63D5P3bhtpcP4mBhPS6fGYPfXz0DL9+2AC/cMh+PXjEV85Ii8XZGMa74xy7c8WoGiuvahp1Du9mGQ2frsfVYBb7MrcSR4oYx9W4mz3JpRVgIcSmAJwEoAbwgSdLv+70vut6/DEAbgO9KknRI5rkSnde6O0csW7bM21Pxmpnx4fjv/rMwW+346drx0wfXXRanOjpqvLKnEKdrWvG9ZSlenpH7LUyJQmSQGp/lVOBgYT1UCoFvL0zw9rRG5NYlyXh5dyHePliMdTOi8XZGMW69INnrLe+UCoGH1k7GLS8dwJsHHHN6fudpHCisw1+unyXb7nbOTI0Nw88vm4JffngMr+4twq0XJAMA/vZlPv7+VQGunx+P3189c9Q7062YaMAbdyzGzS8ewLef3Yv37r6g5/O22Oz43ZbjeGn3GcSGa/Hn62Zh/axYpyUq31uWgsZ2C17ZU4hnvzmFtX/bgcfXT8P1C/p+B602O7bkVOCtg2dx4EwdLLa+5SYapQJLJuiwcUECLhnmgb7WTitaOq0ICVDx4WEvGPYTF0IoATwNYA2AEgAHhRCbJUnq/eTPOgDpXf8sAvDPrv8loi7sHAFcJu3C9xV/RGxADawZcYDucWDm9d6els/ShQTgrshM3JT7AP4aUANpdzwQ9Ni4/sxUSgV+GpuNC/P+iRjU4p4gI3RnfnNe/czJ+mD8NCYL1+z8IUw7a7BLo0Oo6dcApnl7aliersePTEdwydYHIG2txXpJB3XiHbhqzmVuv/dNi5Ow7XgVsj99HpY9H0DVUobr7DoYJ9yNG66+bMzbM0+LDcebdyzGtf/cg9ee+xN+pnkboqkUDQoDajquxa1LvoOH101GkGbo6BMeqMYDq9Nxzbx4PPRuFh56LxvKnHdwdcOLEI2l6AiKwV/sG/Fcw3wkRAXie8tSMC8xEtHhWkgSUN7YgcyiOnycXY67/3sIk6ND8fdpBZiY8xegsQRSeByOpP8Qz9bPw4HCOtS1mnvuHR2mxYKUKFw5OxYXdX4N5fZfA40lQHg8zBf9Ann6S1Hd3Am7JEEfGoAp0WEIPP4esO1XPedh9aN9/33Jfnvw9529B/Q9ln4JkL/V+evArgdY2+sACAAjaG33+OAPbnqScFY03+cEIZYAeFySpLVdrx8BAEmSftfrnGcBfC1J0htdr08AuEiSpPLBxp0/f76UkZEx9p9ghP70pz+hrW34v+ogkltubi5Onz6NgwcPensq3pH9NuybH4DC2mvbYHUgsP7v51XI8ajst2H+4D5opF5bU4/3zyz7bdg+vB9KW68Hl863nzn7bVg/vB8qX/wZst+G7cMHoLSd+/dQUgVCfMszc2s88F9otjyIQJwLf5I6EELGz6bgy5cQu/NhBIlz97AqtVBt+MeI72G12fHRf57E2tO/7TNeBwJwYuFvMOPS2wcN8Da7hE+OluPwx8/ip+Zn+lzfJmnwe9U9aJ10NSaaQhCiVaGx3YITFc3YlV+Dpe1f4Y+aF6FFZ59rfmb5ATbbz/2t4lWq3fi9+gUE9Po9olME4Kng+/GF6kJcjl24s/FvfX4PsSq02Drhf6BWKbDy5G/6fE/tQgUIAYX9XM9pCY6IO9jrMfFgGBZCZEqSNH/AcReC8LUALpUk6Qddr28GsEiSpPt6nfMxgN9LkrSr6/U2AA9LkjRo0mUQJn9TXV2Nt956C9XVI9vhadz463SgsXjg8fAE4Ec5np/P+cAfP7Px8DP78s/g7bl54v5y32OM49n/Mg2KpoE7DkrhCRBOrrfY7LD83zQEtTvZXCYoFidv2AuFEChvbMeiDy9EuLliwHk1SiMeSXoDTxTeAKO9asD7JXZHx554hfOdHD3mPAnC1wFY2y8IL5Qk6f5e53wC4Hf9gvBDkiRl9hvrDgB3dL2cBODE6H+kMdED8PL/++MCP0d5+MXnOC9GMW+w9zLL7XI8RTjuPkcPfGbOePVz9NLPLCtf/hm8PbdR3n9E30m5f8axjjea6129ZrjzhnrfF3j434ckSZIM/Q+6UpVdAqB3lXg8gP5/THHlHEiS9ByA51y4p1sJITKc/amARoafozz4OcqDn6M8+DnKg5+jfPhZyoOfo3OutE87CCBdCJEihNAA2Ahgc79zNgO4RTgsBtA4VH0wEREREZG3DbsiLEmSVQhxH4DP4Wif9pIkSceEEHd1vf8vAFvgaJ1WAEf7tNvcN2UiIiIiorFzqWGdJElb4Ai7vY/9q9evJQD3yjs1t/J6ecY4wc9RHvwc5cHPUR78HOXBz1E+/Czlwc/RiWEfliMiIiIiGo+4xTIRERER+aVxG4SFEJcKIU4IIQqEED9z8r4QQvy96/1sIcRcb8zT17nwOV4khGgUQhzp+udRb8zT1wkhXhJCVAkhnDad5PfRNS58jvw+ukAIkSCE2C6EyBNCHBNC/NDJOfxODsPFz5HfyWEIIbRCiANCiKyuz/F/nZzD76MLXPws+Z3sZVxuas1toeXh4ucIADslSbrC4xM8v/wbwFMAXh3kfX4fXfNvDP05Avw+usIK4MeSJB0SQoQCyBRCfMHfI0fMlc8R4HdyOJ0AVkmS1CKEUAPYJYT4VJKkfb3O4ffRNa58lgC/kz3G64rwQgAFkiSdliTJDOBNABv6nbMBwKuSwz4AEUKIGE9P1Me58jmSCyRJ2gGgbohT+H10gQufI7lAkqRySZIOdf26GUAegLh+p/E7OQwXP0caRtd3rKXrpbrrn/4PMPH76AIXP0vqZbwG4TgAvfdDLMHA35xcOcffufoZLen6a5hPhRDTPDO1cYffR/nw+zgCQohkAHMA7O/3Fr+TIzDE5wjwOzksIYRSCHEEQBWALyRJ4vdxlFz4LAF+J3uM1yAsnBzr/yciV87xd658Rofg2LZwFoB/ANjk7kmNU/w+yoPfxxEQQoQAeA/Ag5IkNfV/28kl/E46McznyO+kCyRJskmSNBuOnWkXCiGm9zuF30cXufBZ8jvZy3gNwrJtC+3nhv2MJElq6v5rmK5+02ohhN5zUxw3+H2UAb+PruuqH3wPwH8lSXrfySn8TrpguM+R38mRkSSpAcDXAC7t9xa/jyM02GfJ72Rf4zUIc1toeQz7OQohooUQouvXC+H4TtV6fKbnP34fZcDvo2u6PqMXAeRJkvSXQU7jd3IYrnyO/E4OTwhhEEJEdP06EMDFAI73O43fRxe48lnyO9nXuOwawW2h5eHi53gtgLuFEFYA7QA2StylZQAhxBsALgKgF0KUAHgMjocY+H0cARc+R34fXbMUwM0AjnbVEgLA/wBIBPidHAFXPkd+J4cXA+CVrk5FCgBvS5L0Mf+bPSqufJb8TvbCneWIiIiIyC+N19IIIiIiIqIhMQgTERERkV9iECYiIiIiv8QgTERERER+iUGYiIiIiPwSgzARERER+SUGYSIiIiLySwzCREREROSX/j+7wOGMOJB/AwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.rcParams[\"figure.figsize\"] = (12,9)\n", "fig = plt.figure()\n", "plt.plot(res[:,0],res[:,1])\n", "x_grid = np.linspace(0, 1., 1000)\n", "plt.plot(x_grid, (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), 'k-', linewidth=1)\n", "plt.plot(b_points[:, 0], b_points[:, 1], 'o')\n", "plt.fill_between(x_grid, -1., (1 - x_grid + 0.05 * np.cos(11 * np.pi * x_grid)), color='gray')\n", "plt.ylim(0, None);" ] }, { "cell_type": "markdown", "id": "d39bab91", "metadata": {}, "source": [ "We can clearly see the dissipative behaviour of the inelastic collisions at play. Indeed, despite the fact that we asked to integrate the system up to $t=10$, the integration was stopped at" ] }, { "cell_type": "code", "execution_count": 8, "id": "f75d9d74", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.7390191265468795" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ta.time" ] }, { "cell_type": "markdown", "id": "6a619f29", "metadata": {}, "source": [ "instead, and the outcome ``oc`` informs us that the integration was stopped because of a stopping terminal event:" ] }, { "cell_type": "code", "execution_count": 9, "id": "935109c2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "oc" ] }, { "cell_type": "markdown", "id": "09145c42", "metadata": {}, "source": [ "A value of ``oc`` negative and greater than ``taylor_outcome.success`` indicates a stopping terminal event. The index of the stopping terminal event is ``-oc - 1``:" ] }, { "cell_type": "code", "execution_count": 10, "id": "ba9dc815", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stopped by terminal event? True\n", "Index of the terminal event: 1\n" ] } ], "source": [ "print(\"Stopped by terminal event? {}\".format(oc > hy.taylor_outcome.success and int(oc) < 0))\n", "print(\"Index of the terminal event: {}\".format(-int(oc)-1))" ] } ], "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": 5 }