{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c6805951",
   "metadata": {},
   "source": [
    "# Multidimensional differential equations\n",
    "\n",
    "- Börge Göbel "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "921259c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29eb17ec",
   "metadata": {},
   "source": [
    "## Ball in a bowl: 2 uncoupled harmonic oscillators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b461fb52",
   "metadata": {},
   "source": [
    "We solve the following differential equation:\n",
    "\\\\( m\\ddot{\\vec{r}} = -\\xi\\dot{\\vec{r}}-\\nabla U(\\vec{r}) + \\vec{F}_\\mathrm{ext}.\\\\)\n",
    "\n",
    "For the profile of the bowl we consider a quadratic function\n",
    "\\\\( U(\\vec{r}) = U_0 r^2 = U_0 (x^2 + y^2). \\\\)\n",
    "\n",
    "The gradient is \n",
    "\\\\( \\nabla U(\\vec{r}) = 2U_0 \\begin{pmatrix}x\\\\y\\end{pmatrix} = 2U_0\\vec{r}. \\\\)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "34940b2b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.contour.QuadContourSet at 0x24af2a064c0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAADyCAYAAABkv9hQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAB1MUlEQVR4nO2dd5gcV5X2f7e6e/Jo8ozCSBolK2fJNgZssDEYJ8nGOBINu/4A42WXZMKCF7BhCQvLLp8JS/DyESVHbDA2YGNwDkqjHGY0SZqce6ZD1f3+qLrVVdXd0z1Zod/n6empfKu63nvCPedcIaUkgwwyOLOhTXcDMsggg8lHhugZZHAWIEP0DDI4C5AhegYZnAXIED2DDM4CZIieQQZnAfwptmfG3jLIYPIhJvsCGYmeQQZnATJEzyCDswAZomeQwVmADNEzyOAsQIboGWRwFiBD9AwyOAuQIXoGGZwFyBA9gwzOAmSInkEGZwEyRM8gg7MAGaJnkMFZgAzRM8jgLECG6BlkcBYgQ/QMMjgLkCF6BhmcBcgQfRogpSQcDhONRsmU285gKpCq8EQGEwzDMAiHwwwPD9vrfD4fgUAAv9+Pz+dDiEmvQ5DBWQaRQqJkxM0EQUpJNBolGo0yPDzMsWPHKCgooLi4mOzsbKSUNsFDoRCFhYVkZWVliH92YNJ/4AzRpwBKVTcMg7a2No4ePcr8+fMJhUJ0d3czPDxMQUEBJSUlFBcXc+TIEWpqasjLywMyEv8sQIbopzui0ShNTU1EIhH6+/sJh8OsXLnSJcGllPT399PT00N3dze9vb2UlJRQUVFBcXExWVlZGIZh7+/3++1PhvhnBDJEP13hVNXr6upobGxk4cKFzJ07F4BwOJyUoHv27KGyspKhoSG6u7uJRCLMmDHDlviBQMDlxPP7/bbE1zQtQ/zTD5P+g2WccZMAwzCIRCLouk5LSwtNTU3MnDmTefPmAaT0tGuaRkFBAVVVVdTU1GAYBn19fXR3d9PS0kI0GqWoqIji4mKKi4sRQhCNRgEQQrgkfob4GUCG6BMKKSW6rhOJRIhEIuzfvx+/38+iRYsIh8MJ9xdCIIx6pFaT9LyaptmkBtB13SZ+U1MTuq7b24uKigCIRCJAhvgZmMgQfYIgpbSleG9vL/v27WPhwoXMmjWLkydPEgqF7H2FEAghTMluBAkM/R/CeQ+DyLa3jyT1fT4fJSUllJSUANjX7O7upqGhASllhvgZuJAh+gRAjY0bhkF9fT0dHR2sX7/e9pqPRFwZfRYpW5Ch7yByPj2m6/t8PkpLSyktLQVMB6Aifn19PUIIiouLKSkpYcaMGUQiETo7OxkYGGDOnDm2je/z+TLEP0ORIfo44HS4hcNhamtrKSwsZPPmzWhaLOhwJKIbkQcJyyh65Nf4sj+CEAUpJXoq+P1+ysrKKCsrA0xp3tvbS2dnJ8eOHUPTNLKzTe1h1qxZtqmh2qoce36/39Y+Mji9kSH6GOEcG+/s7OTQoUMsXbqU8vLyuH1HIm7EeA0NA8kQ0eBXyMn/mn3+iUIgEKC8vNxuWyQSoaGhga6uLnbs2OEyBQoLCwmHw7apoWkagUDAlvgZ4p+eyBB9DFAON13XOXLkCAMDA2zatMmWkl4kI7qunyQiu1GjmLr+CAH9U5NOpEAgQFFREUIIFi5cSDgcpru7m9bWVg4fPozf77eJX1BQYBNfCIGmaXGqfganPjJEHwWcqvrQ0BC1tbVUVFSwcePGEcmZjOhDkQeJSIkE6xOlJ/gZhPiXKU12ycrKoqqqiqqqKgA7Yq+lpYX+/n6ys7Nt556S+GoUIUP80wMZoqcJNTZuGAYnT56krq6OlStX2kNeIyEZ0QfDjyMlSIQdmST1Z0C7AZg5oe0fDbKzs5k5cyYzZ5ptGB4epru7m+bmZgYGBsjJybGde/n5+RninwbIED0F1Nj4kSNHqK6u5sCBA+i6zrnnnksgEEjrHImIbhhRQvI4KihKSvPbQMKM/0TK+yb0PrwYjcaQk5PDrFmzmDVrFlJKm/gNDQ0MDg6Sm5trR+3l5ua6iK/rOj6fj/z8/AzxpxEZoo8A59h4U1MTJ06cYN68ecyZM2dUdnQiog/pOwjLCODDAKQ0CSAB6T/GkPEaxVw8cTeTpF1jOSY3N5fc3Fxmz56NlNIO1a2vr2dwcJD8/Hyb+H19fQwPD9tRgcq5p+L0M8SfGmSIngTOsfGmpiaGhoa44IILKCgoGPW5EhG9LfhrIoZFcpQ0F2Cp8a3ym8yUbz7lPdxCCPLy8sjLy2POnDlIKRkcHKSnp4djx47R19dHIBCw7fzs7GxCoRChUAgppUvNV8N5GUw8MkT3wBnGGo1G2bt3L1lZWRQUFJCfnz+mcyYiem90BwaBmG1uX1+RvpUTQw8yO+/aMd7J9EAIQUFBAQUFBVRXV9PS0kIwGLTNH29KrqZpmSIcU4AM0R1wjo339PSwf/9+Fi1axMyZM3nhhRcwDAOfzzfq83qJHow2E5L9mNJbc3jdhU10ieBY8PtU5lyGX8ubmBt0YCq9+jk5OVRXVzN37lxXSu6hQ4fsIhuK+EKIDPEnARmiW3Cq6nV1dXR2drJhwwZyc3MB07YcKzm8RG8ZfJyQ9CHxIYnZ5gpSmuRHBjnQ/31WFf3LGO8qdbsmG868e3XNGTNmMGPGDObNm4dhGPT399Pd3c3+/fvjUnKFEAwNDdnnyBB/bDjrie4cGw+FQtTW1lJcXJwwjNUwjDFdw0v0xsE/EpHmo5c41XVhLWsgzG3Hhx9nft4NFAbmjPUWpxVSyhEdbpqmUVRURFFRUdopuRnijx5nNdGdY+MdHR0cPnyYZcuW2THiTmiaNiFED+sD9OknAT8SgXlGkwiGdNjqhlLhDV7q+QaXVHxnTNeebngleiqMNiXXSfyTJ08yZ86cDPET4KwkutPhZhgGhw8fJhgMsnnzZrKyshIeM55EE+exTcHnCEuJtCS6LcUtWx1n8Iw0/2+LHOT40AvMzz1/TNdPhKmy0UdLdC/STcktKSmhqamJWbNmuSR+pvqOibOO6E5VPRgMUltby8yZM1m2bNmIL8FESfTDA48TkgGrLcJ2xAEue91cp6H4+Gz3vczOXk9ASxxPP9Z2TTYMw5jQsfKRUnKHhobYsWOHKyVX13W7fr7KxT8biX9WEd2pqp84cYLjx4+zcuVKuzjDSJgIiR6ODnNi6CBS88VscbDH0p0e9zjyy35e6P0Fbyy5dUxtmC6MV6KngjMlt7u7mzVr1sSl5CriFxYW2sSHs6sIx1lBdKeqrus6+/fvB+Dcc8/F70/vEYzX6x6NRvnTrt8yXCLBCLiI7JLqtmMOsJ1z5rpdA39iRf5bKcuqHlM7nDhdVPfRXitRSm53dzft7e0cOXIkLiU3Go2eFdV3zniiqzDW1157jYULF7J3717mz5/PnDmj82KPx+ve2tpKf38/waVHiUY1h8ouMKSwI+PUOqRToqv7MD3x9zd9m62Fn6K0pCztTmqke5psTCXREyEQCFBZWUllZSVAypRcbxGOM6X6zhlNdOfYeE9PD/v27WPt2rVjinAbi41uGAb79+8nHA6TX5jPgehRwoY/Nk4OtiNO/Q84VHhcDjppCEKinZd7/0LV8SUIISgpKaG0tJQZM2acknHjU0n0dK6TbkpuMuJ74/RPF+KfkUR3OtwikQh79+5FSsl55503ZjKM1kYPBoPs3r2bWbNmMWfOHJ7c8RCD+jDqkUspYuq5naqq1rmTXNz3JtihPc0/rb2cHPJs6XTo0CGys7Nd0ulUeAmniuhjNUWSpeQ2NTUlTckNhUK0tLRQVVVFXl7eaVF264wjujOMtbu7mwMHDrBkyRKGh4fHJfFGI9FbW1s5cuQIq1atoqioCMMwOJa9i4jhi5HZUtm9Q2kuCe5R62P3aPD/Tn6fD839tEstHR4epquri4aGBgYGBuwsstLSUjvCz/mcpgJTRXTnTDbjQbopuW1tbVRWVrqq7yiJfyrm4p9RRFcSXIWxdnd3s3HjRnJycjhy5Mi4Xrp0JLphGBw6dIhgMOjKV5dIOnzN6NJnS3LlZU9EdoWE4bHWd/1QHc91/50LSt5gb8vJyWH27Nl2+ujg4CDd3d0cPnyY4eFhCgsLKS0ttcekzyQbPVUE3lgwUkru0NAQu3btcqXkOnPxb7/9dj7/+c+zbNmydK4zF/hfzGojBvBDKeV/CiFKgd8ANUA9cL2Usts65jPABwAduENK+ceRrnFGEN07U2ltbS2lpaVs3rzZfsmURB5LUorz+GQYGhpi9+7dVFZWsnTpUtfL3RA8SlBEwfC7yCwdaalO77pLstv7x0v2B9oeYEXBKooDxXHtcWaRzZ07144p7+rqorm5meHhYXJycsjNzaW4uHjcjr1kmEqJPtkS1JmS29LSwsaNGwkGg3ZKbjAYpKCggKNHj9LV1UVOTk66p44CH5dSviaEKAReFUI8CbwP+LOU8mtCiDuBO4FPCyFWADcCK4HZwJ+EEOdIKfVkFzjtie4cG1dDKMuXL7cDKhTGS/SRJHp7ezuHDh1ixYoVtrR04i9dTxCV5nUVoWP2eXKpbl8t0ToAIny/8afcufCfU7bfGVO+YMECmpub7Syy+vp6NE2z7fuJdOydSUT3Qk2dpVJypZQMDAzw9NNPc/DgQa666io2btzId77znRFLjkkpTwAnrP/7hRD7gTnAFuBN1m73AU8Dn7bW/1pKGQLqhBBHgHOB55Nd47QlundsXKmnycJYxxPZlux4wzA4cuQIfX19Sa8bNaLs7z9k2+cx29u6DxKo7NI9tu5ch+M4gGNDDfyh7SneXvnmUd2PEIL8/HzXpI89PT2cPHnS5dgrLS0lPz9/zGQ9Y4guhwA/CBXVGN/pCyEoLCzkIx/5CPfffz9//etf2b9/P4WFhWlfRghRA6wHXgSqrE4AKeUJIUSltdsc4AXHYU3WuqQ4LYnuLPE0ODhIbW0ts2fPZvny5UlfqvES3SvRQ6EQu3fvpqSkZMQqsIcHDlveds1NcAeZvcSWVupavMc9sc3+YNsTLC84h5q80ccGKGRlZbkce8oWPX78uO3YU/a917E3Es4EosvoMfyR/0TP+U7ax0QiEXJzc9m4cWPaxwghCoD7gY9JKftGeG6JNozoQDrtiO4cGz9x4gQNDQ2sWrWKGTNmjHjcREr0zs5ODhw4kDTTzYnH2v5E1PDFElckDoI7fy/hUOndEtxNbGGvUOsjGHy3/pfcs/QOcnwTEwvvdUIpx54qFqFyxktKSpImApn3cXoTXRohtOF3oWvngyevfiIhhAhgkvwXUsoHrNWtQohZljSfBbRZ65uAuY7Dq4GWkc5/2hDd6XDTdZ19+/ahaVraYawTQXRd1zl69CidnZ22N38kDOvDHBqot9R2dR9aHJmTet1lvAPOXBaulRJoCXXxP02P8uF516CJifc+ex17KnW0ubkZwzAoKiqitLSU4uJilx/kdCa6lJJw8Fp8sgfNfxmaY32q40ZZPFQAPwb2Syn/w7HpEeC9wNes74cd638phPgPTGfcEuClka5xWhDdOTbe19fHvn37qKmpYfbs2WmfY7xE13WdhoYGqqqq2LRpU1ov1SvdewnqEUDDcKnqsZfA64WXHmmdTLLb6ayOlU917GBtwTm8sWx1yraNZxzdmTO+YMECdF2np6eHrq4u6urqbMdeaWnplDnJJuM6A8FPoRlHkMJPwPc6e70qYT0SRkn21wPvBvYIIXZa6z6LSfDfCiE+ADQA77TOvVcI8VtgH6bH/iMjedzhNCC6ruvU1dVRUlJCe3s7ra2trFu3zp6pNF2Mh+jd3d3U1dVRVlbGOeeck/ZxD594iohhedtl4uCYVFJdrUPtH3cV4eo4vl33EAvzZzMnZ2STAiZO/fT5fK5JHVU8+YkTJ+zQ4/Lycju6bDIk/EQFzCj0D/2SocjDZAkwRDXZWqz6byqiR6PRUY3uSCn/TmK7G+CSJMfcDdyd7jVOWaI7VXWlIhYXF3PuueeOqeceC9GllBw/fpzW1lYWLlxoxzyng85QDw3DbehGTFVPJdWlYz+vM04mfQ8UzHNEpcHnDvyC7636R/L9aY/jTiic8eRDQ0MsXryY/v5+u+67qgI7WsfeSJhIiR4M76Qr9GUCmOzz+dxlt3VdH/Faqrb9qYRTkujOsfGuri7a29upqalh0aJFYz6nz+cbFdEjkQi1tbXk5OSwefNmOjo67NlH0sHjJ19kWJeWAy42G4tLkjvI7CRycq97MpvdjRNDvdx96GHuXn79tMdeSynJzc2lsLDQ5djr6upyOfaURz/d2W+8mCiih6IdHB+4lYAATZi/w4zA2+OuNZLEVp3ZqYRTiujeEk/Hjh2jt7eXWbNmpfSqp4JypqWD3t5e9u7dy8KFC+1kh9GkqUaMKL9vfYWoodmkNmzbO4l97g2KsSvBJrDN0zCvX+w+ws+O/53317wx4fapDE31VoFVjj1VBbavr4+uri6ampowDMNOIvE69lJdZ7xE1/Uw+3rfgUYUgYaGRIo8Ar5Vnv1GVt0zEn0EOMfGh4eH2bNnD+Xl5WzatIljx46lTdJkSEd1l1LS2NhIc3NzXDrraApPHBlooTM8iEHsZTASSHW8Ul26ieeV7LGrx5M9Uct+3vQc8/LKuaRyecJ2ngqx7t5ikNFo1HbsHTt2zC4UUVpaSmFhYVIyj1eiSyl5ufs6YJCA8OGTZniT37ceTXOTOhXRBwYGMhI9EdTYuJSS1tZWjh075gonHa3anQipiK5mZfH5fJx77rlxP+RoJPr/NjxNVFWKUZLaKdVHsM/jnXHJJLjb655oFx342qHHmZNTwrIZ0zM762g1B7/f76oQoxx7Kl88JyfHtu+djj3DMMYVr/9c5wcIGZ0EhA8NSVj4QEJF4Iq4fVPZ6CrA6FTCtBLdq6ofPHiQcDgcN1PpaNTuZBiJ6P39/ezZs2fEIbt0JXpvOMir3fVEdWvs3EFcw+k996jtCgnt8wSeeNdmmZxIA3qUj+26nx9vuoU5ucUp2z8ZGI/m4C0UMTQ0RFdXV5xjLxQKjRi4MxL+2vFx+qN1ZAuBED40S5qHhZ+CwJvi9k8l0VVyy6mEaSO6c2xchbFWV1dTXV0d92L4fL5RebwTQdO0hOdoamqioaGBNWvWjPjjpCvRH2p5leGojnRknjkz0KRDqidKWrG3xanmiR10jhYm3dIdHuYjr93PDzfeQGVOgXWNsY+jTydyc3OZM2eOPaHjwMCAXROura2Nnp4eW+Kn49j7S+eX6IjsJQsNTQQQMkJE+PBh4KeagC8+SSkdZ1xGouMOY21ubqa5uZlVq1YlDf73+Xyu+bjGAq9EV9F1Usq0ouvSkehhI8q2hpcdTjhHVddkansi4kNSKZ4+P9V1zAOagj38846H+e+N11KSZQ5pTbdHfrxQSSSFhYX2VE5+v5/u7m4aGxtx1nxP5Nh7suMbtIReIVsIIiKAjwgRaanu+CnLvijhM9J1fcT35az3ujvHxqPRKPv27cPv9ye0iZ0Yb1Sb9xwDAwPs2bMnqQaRCOlI9B1dDbSHB3DWY7eruhputd06q4u4cVI8IdnFSGI9Ucvt/w70tfGpHY/xvU1bR3OC0wLKRndO9pDIsaeG8f469D80R14gIHwIJAGihAwfCIkPHUSAyuy3JryWrutkZyfPKRgcHIxLk55uTBnRnWPjvb297Nu3z56pNBV8Pt+E2egnTpygrq4urUQY7/EjSXRDSv7v4b9iGFatdsPtZXfllLvI7SRzPImT2d8jdRAj4eXOJu549VE+VbX+1PDEThASed0TOfa6urq4v/3b9PgP4ddM6R3Bh6oK4EPgEwGEzKXAnzgKMp2AGZX+e6pg0n9rp8NNSkl9fT0dHR2sX78+7TDWiZDoYBaIUGWeRuuhTSXR6/s72N/Xim5obnvcM3zmtMtxrI/9n4C4SZ1tqjNJ+zYQwAvtDXw5GOLORecxK/1DT2mkM7yWlZXFY9GfcVKrIwsNIX2YIVBRBD40IYkCmoxQ4l+HT0ts42cCZjxwqurhcJg9e/ZQVFQUN1NpKoxXogeDQQ4ePEggEGDt2rVjsk1TSfTvH3mWsC5xquPKy56I1G5CeyR5IpU9hec9XXVeCnPfV/va+LdDL/CDmbMoDEzcNE/ThVSx7oZh8O3jX6Y72khA08zKbJopxaNoaMIMH/YLg7D0M9S8mIMdBxM69tIJmDlriK4cbi+99BILFy7k0KFDLF261FajRoPxEL2trY3Dhw9TU1NDT0/PuIpDJpPorUP9PNN2zBwzlziI7SC9faiD1CN43hNnsKQXFZcu9vR38OHnHuXb572d8pzRJQmdahgpMm4oOsTddV8iJHvwCx9CmppNxPChaRKBH02aUt1nGASEnwtWvJOh/qjLsadInypp5VQcXpvw/EElxUOhkD10dvz4cTZt2jQmksPYAmYMw+DAgQM0NjayefNmu+zyWDGSRP/p0RcZjESQhkAaZmy7rgt0HQxdYOgCDM3cZggMQzM/UiB162MAunM50QezEzEm7rOz8yQ3P7WdI31dY342pwKSqe4NwWY+ffiL9EYHiBoaUWl+ItKHbn2by2YqcUT6yPUtIsdfSElJCQsXLmTjxo2sW7eOoqIiOjo66OrqYv/+/dTX19PX1xf3XowmYObWW29FCNEmhKhV64QQdwkhmoUQO63P5Y5tnxFCHBFCHBRCvC3d5zOhEt0Zxjo0NMSePXsQQoxYaikdjDZgZnh4mF27dlFRUWFXZI1Go+MuJZXo+I7hQR5p2I9hKGmrJHaCITWHVHe9Gx613XPl5Gr5eKW7dakTwQE++PRDfHHNG3jT/CWn5bBbIqI/1fEivzh5P4IofuEzxZqtsmM64SRoQuKTGhE0hPBTkxtff8/p2AsGgyxevJiBgQGam5vp6+sjNzeX0tJSe2w/3Tpx73vf+/jpT396GWa5Zye+LaX8pnPFWKq/2u1PqzVpQEpJKBRCSsnJkyepr69nxYoV7Nu3b9wvzmhU946ODg4ePBhXCXYiKswkkui/PPoaPaFhnHnhccR2qeixdTF4bXTPRZIFyng7hVTE91zS2dt0hIb56Et/5vpDB7h29kJ7auJRlCyeVjiJbhgG367/Da/17cInwG+tF9b9CgNLZdeISokwfPg0Aw0fEalRk3vBiNfSdZ3c3FwKCgqYOXMmqt57V1cX99xzDzt37uSTn/wkl19+OTfccMOIQ3EXXnghQLrq1BZGWf1VYcKILoSwZyo1DMPl2U4a72z0QOQZpLYYEViR9NzpJqQcOXKEnp4eNm3aFPdwxxtGm6jcc11HK/cdes0KfBEOu9ohsb3EdpIz2RCZTCLFkxI58TnjdonrbUTcPtt6WujODvCRgnzaDhwgEonYZaJKSkrGXC57sqGIfmKoi88d/DF9eg9+TQDK8QYYEjQrS0D6EUQRhsSnGUSlhiYk+aKCQv/IJqZXe3DWe7/33nu56KKL+Od//meeeeaZ8Tyv24UQ7wFewaz53s0Yqr8qTKjqXltbS3l5OXPmzElz4oRcjKGfoMsIvqKH0bTEzUmlETgrsm7atCnh/hM1RKfQ0tLCD3b9nUE9ahHTU9rJS+xEpPYS2kXSeLKL+DS2BEj2rGTyfRxNlcCfWo5zsLeHz29+A+dVzqa3t9eOL9c0zZb2hYWFI/42Uxlmq+s6v2l+jl+e+Aug49c006ehWSy3voQBQgMhdcvbLk2SSzNjbUne61O+b6kSdZSgO++888Z6O/cCX8b8Ob4MfAu4lTFUf1WYUKKvX78+7sdVandComvZRLTF6JGHiPTfRl7Rj0d9TVWRNZVHf6KIrpx87cEBfh88iTRESnLbj8RIIHkdEl7Eedon0j5P8mI6JL1zj8b+Pv7lmSfZsnApn9n8etsMUkEnTU1N9Pf322WgE6n5U5Xz3hzs4lvR5+luHMCngU8TRA0Nn22Tmyw3ve3S/JYGGhq69W1Ic1htWf7rUlxtZMEzEZ2blLLVca0fAY9ai6Ou/qowoURPpN76/X6i0WjSzCJ//pcJ9TyIHvkLIrST3Ox1aV1LSsmxY8fSrsg62tlQE8EwDF5++WUqKyt5KniCUFRPTG5lUzsltofYLlKnsrXN4XnPeUbYPxW8tnqi7RIGwxF+eaCW3x09xA8uuYJlpWVkZWXZs486q8Uc8Kj5xcXFkz676LAe4TuH/8zj7bsBHZ+wVHVACumS4mgGwhCg+fAh0a2htaj04bOkei5FlGdVj7td471vVeLZWrwGUB75UVd/VZhwonuRypHm8+WiaxcSiT5FsO865lQcSXkdFXxTUFCQdkXW8b5w3d3dBINBNm7cSI8fHnx5X0xCywRS20vshOPojgt4/3csi5T14kiP7N7OIslx0n1x+sNhbvn9Q1w6fwE3Ll3JuqqZBKy5wZ3VYnRdd6n5QghCoRB9fX0p1fzRYEiP8EDjLv6n7jlCDKNpEk0IJBo+K+lfOlR2YUl2U2U3iOg+8Ek0qeGTEl0KBD7Oyd+UVqnsdMo9p4ubbroJTGdauRCiCfgi8CYhxDrMX6ceuM0676irvypMeghsOh7zoqIf09KxEIMo9V23U1P63wn3k1LaZZ6WLFlizyoymXAWiMzLy6O0tJQ7//oww1E9Jrld5DaJHJPYiaS649tWmxOQwOszG8HyEM79R9gpHX+e3Ral1lsHPVlfx5/q67htzXrWVM5k48xZ5DpCiVXSiFLzBwcH2bNnT1pqfjoY1iP8un4XP294hb7oIEKAEOZoiNQM7P7esJKKfIZ5E4ZmDqlJDSElAoFmaGiaRJcGutQQAtYUprapU0XghcPhUdW9+9WvfsWvfvUrbyRyUht2tNVfFSad6H6/PyXRNU0jkH0LwdAvCUceIxj5DHmBeGdifX09ra2to4qTHw9U1Rm/38/mzZt54YUX+GtLPX9tOo7NGofktsntdcClInYSCRtvs8PIbE2NhK+os5dwNiuBii+B7+/eAcC7l6/ijXPnsaaiirwEL3cgECAnJ4cVK1akVPOT5R5EDYPGgR7uq3uNx1sOEZJhhGba/kKTaGq+eUMzWycMS3vXEDoIHwgMooYPIXQ0QzOPkxqaNDCkhi4lucxgXs7ClM/vdMxFh1NAdVeoKrqbA62/QWKwu+MKzpv5GsJSoyKRCENDQwwODo653PNoMTg4yO7du5k7dy7V1abd1heN8N+vPRcjrJPcign22DmJiZ2A1G5CJyByMu97gm1pQ8S+pbP9jn+FaznWEagO4Of7avl/+2oRwMLiEi6snscbq+exurKSLJ/P5YxLpeY7vfnZeXk0Dvby9Ml6Hm44QONQD4YwEMI8nzRA0wQYEl1INOForKYhjRjZ7QAZAbpNcoEuBT5pjqNrUrI0f3VaavvpGOcOU6S6R6PRtPYtz7+bloF/JcAQL7V/lvMqv2ar6rm5uSxevHhKSK7i41etWkVRUZG9/k/dbRzq7TQJbkC8VI994oidiNRJxtQTSvIEhBYjbEsKh53utQziFiwC2R2BpcoLez9lD8PR7m6O9nTzmwP7iBg6586aw8XVc5mlR+kcGiI/EMCvafiENQwpBMUlJWQVFJBdVUFDbw+PtjTw0r4XOTrUywBR0yIShnVpYT4uTSCEgYFJepWfApjDaThGWDRMBxymNAdhkdwsF6VLDZ80x9DPLUpvKOx0rAALp4iNrlBecAN1/V9BojMY+gM7jr+Z4ZYZrF27lkOHDk3oOHgiqKCb3t7euGmQXznZzB+7Wy0Suj3sKcntiJBzbXdu8673LKdU2dMle7zgdi9IdyNsh7+nA5BqpcCsQmk1OKhHQINnmxp5trkRAZQc2sdAJMw5JaW0BYPkBPzkBvy0BgfpjQyjW04NqRFTVzTMOQ1ND5rpRceU4lLTEIZ5jAE22TVAV89T00AaRNEISMOW5rrhQxM6PikwpMBAkC1msCRvcVqP73ScvAGmQHX3+/2jmvhgaclvqO26Hr/PYG/0K1y78THysvInpPgEkHRsV3nyCwsL42Lzg5EIn3/2L4QN6Sar9VIJj8MtmSPOLc3VsY51iUJhveSHEUk9pogK12Udertwqu+xPaSIJ785o7OT+DFp3xkcAgF7Otod55D2NUxNQSCkRGrCPNYQSCEtsgvrI81rq/0sb7ouLbJrHrIjEBropmqAIYSZyyOFJc0lUcPPshkr0h4ROB1z0eEUU90BRGQuQ5FitMAAJf4hftXyUW6d/7MJrTLj/aH6+vqora1l0aJFdrVRJ+5+4RnqenpibDA80ltJdCe546S5tez4P04z8G5PsEyi5RQY8RV2quOOda4KNk5WW5wz/3c0Wo/tYkr8GOmFg/T2OaznJJEmO61rWBq2Jd2ta2tq2gvzrwSzk5WYkl8z04A1BIZm2NGuGBpCmJ53IYUtzXVD4BOmNNclvLU88SQXiZBOTfczXqInwmgI2tzczPHjx9mw/AFeGrwcTWhIo5Htzd9klXbluFV3le7q/KGam5tpaGiIm7BB4bf7a3no0AFbgnsJ7lLLk5HbS2xbEyBJB+BuQ0KyezEG1T3uOHVrTuILR/sdHYBLdbclvLS3OfsGU0qr5yJjBBaWFmgkILwhQQi7UIYUIDRznelok0jUsQLNsttNsuMIlhEITWDYXnZMyW6p7oW+AubmpF/zPuOMS3aBNIbXVDKMrut2MkxOaBP9kR2U+MMcGnwOgzIu1q8aV1ucYbAqlDUcDrN58+aEwzsvtjTxX6+8hK7LePVc5YarZZetjlt1j1vvJH/senGETkTeZE65URA94a7C0Tb7j+O8wt0BKKnuDgwUMYkvZIygSkrbarhHyickfGw/qZmNk4YwpbuhVHkDqZmdgqFZRDbAh6nGC2nW7hPCHB0xpFkLQAozSEaXGpuL1+JPkmORCKcr0SfUhZ1seG0k1X1wcJCXXnqJoqIi1qxZYxPuopnfZ1gGMBCU+gZ5zfgLu/p3jKt9KoNteHiYl19+mby8PNauXZuQ5N3DQ3zrxedoHwqaKxRZdRC6sCW40C3SG+a3UIUhdGGqrdZ2TW2393N8y9g5bNvfua/u+CRYhx47Bu91DPc29NTniztGOu/P0V4ZO14zQHMso5vPQRgidq+Oa2GYz1AYlrquzmdJ6dhHIHTzo5alapPUzGVpbjOsc+rWt2EocisJju2AkwgM6eMtZZtH9Q6l44xLNxd9KjGtqvvJkyc5duxY0oqsyws+xq6+75Hrj5IdGebBgW1U9lazpih5SutI0DSN7u5u6uvr4/LVnYjqOldv/xXtwaBbRXekoioJbktT7zZwO98goRSPU/e96737e7anI8kT2ugiyUaP5E7sbSem0jskv1NlF8TUc+mR9FLE7lda0lup9VKLCQwppSmeLe1IGE5FSZqhrUp9t2LaJSA07CFBKbHUdt0q9WV+DCkQ0QDDDT20l5nzu6VTMDTV1E9nhdc94QUSqO5KbQ6FQmzevDlpyODK0ut4qe9HCAyqsoL0RAv4j6P38rEFH2ZdSeKJA5NBSkkwGOT48eMjJsF0DQ3xuaf/TPtg0GGDi5hU85AZS5ooQjrJndDu9uzj3BYXCJPKKRd3kym2O0mbaF/rNl0kdzrYZIzI9mmcDjZ1sCK9tSysY5Xqbqv2AtuJFyO8jNnvmqkV2R56S02XmkVuiWm3W8Ezaj0CDEOiaQJDSluF1ywHnDkfnsbbq86jqqiKrq4ujh8/jqZp9oSOM2bMSDp5Q6qa7qei6j7lkXHBYJDdu3czc+ZMli9fnnJY4x2zfs5vWt6NT0CBCNJPLncf/hEfnHcjb5t5blrtikaj1NbWous6K1asSEry55oa+ecnH6cvHLJVU1tFV6ql0/Y2SCi9XZI7kcT22OuJnHPOfRVEgnXJ9o2DGGG7srO9ktyxv8vB5tQG4rztJCY9FuGdUl7DtMHB4W2PdQTqUQscmYdCIIgNw0lDWPa+OW4u7d/GPD9WwU6z1Hbsf133sWXWeZRkz7AnfHBO6Hjw4EG7PFRpaSm5uebsNqfjvGswBRLdWdlFRZytXLnSniY3FWZkV6IxlzAnKc8ZpiMaQZDNd49t50D/CT66+KoRQxdVKOu8efOS9sRSSn66ewffevF5DENa9qKH4BAjvrJV7XUpyO3tABIRO5mqnkiyO7eNB4rgzkWvSq+kvIO4rvUO4ns7gzjSq4upYyyPubACZSSW+m1Lb0uaW+q8tCrEYFhj7iiyY2WoCdNRYJgdhqFLNJ/ZSdi2ujBV98X5syjLctvSzgkdlQbY1dXFoUOHCIVCFBUVMTw87IqW9OJUVd0nPJ7UK6FVHvjBgwftiqzpklzhAwv+h5Dhw0CjKqsfTegIJI+ceImbXvgv+sNDCY9ra2tj586drFy5kjlz5iQsPhHWdf75T4/zjReew9CTkFxJCN10HCmnk9NB5XJYKYeTdC/b397jncc499Nj53c57BI4xtL+yNgnznmXoL14nXZG7Fi3cy3JvXkdhd72qv11yzHnPafDWWc76qwO13zGwnLIqd9NWMVAYva4lLHeSAK61LhiZuJKRApCCPLz85k7dy5r165l06ZNVFZWMjw8zNGjR3nttdeoq6ujt7fX9U6NRnW/9dZbqaysxFMBtlQI8aQQ4rD1XeLYNqYKsDAFEn14eJhgMEggEGDDhg1jzkle6buenZEHyfXr5GgDBEU2Qkiah7p46zPf4qa5r+eOpW8GSBrK6iS6lJJXT7Twny+/yCsnWxzeZeEmeBIb3Bvl5pTeiYfaPMuJVH3HvknVdM/yqJ+m9Ehw54mke9llk3sluVovHNstVd6W7A773Zb6ajgNxzFKAzAAKystXrpb65ySHYfNroP0SXPUA2mfS1rblVT3SUGAABeXj87Ho+z3goICampqCAQCrnnbc3NzOXr0KMFgcFQVYG+//XY2btzoXH0n8Gcp5deEEHday58eTwVYmGSiq4qs2dnZLFyYOgVwJJxbfDnPtz6KTxjMzgvTF4kCphMvakh+cvR57jv6Km8pX8YWo5yKsrK4UFYhBB3BQV47doQHDu7nmcbjMZIaCaS4U01PoqLbtjwOKYl7vziV3UPshOq6arN3W4J9nEhGfJlsH5F4h6Qk9673qvCJSO9Q5d3edneHpmx4KWPqvG3na8o2NyW7S4231mM47ATr97Gnm7aXNc4vXUy2L/2ccSeUjZ5Izf/b3/7GwYMHefvb384b3/hGvvnNb9q2fSJceOGF1NfXe1dvAd5k/X8f8DTwacZRARYmgeiq/rmzIuurr7467vphPp+PKyMf5+HAf6ABZYFBBiPZON/QsB7ldy37eETXyGnOZubBQubkFpHvy2YwHOFoRzudoTBhVQLKVl8Tq+qu8WJwSeK49YZ7H+d68HQQzmU8++Pez943CbG9+yZDQnIn2cllazvOHbfeS3y1TyLSq2tqSQjv9AVolnRX5MaS7iQhu4zZ82Z7zIZJCUKa0tz+6QzBjdXpOXETIZEzTqn5H/7wh/nNb37D3//+d1577bWxlsquUmWkpJQnhBCqusqYK8DCJBA9FAqxc+dOiouL7YqsI1eCTQ+apuHX/ZTkLqM1fISinDCB0DBaNAv1NsZ6b0FQj3As3MPR3h5X8AYuMps2oG1LWwRPKMWdBMexzuF5j9s3kcqeZNl1TidSqfKJ4N3HS+xE57BIa+/uPEZzkNBBfi/xnUNqXklvdwCG45yK8M7t1j72kJzmUeUTkR3MITbr95C+WCRjrKyTQEqNmdlFrJgx9qklUwXMGIZBbm4ub3jDG8Z8jSRI1D2n8zYAk+CMO3bsGAsWLGDJktiMH+mEwaaCilP/l8VfIKz7kGjMLRjEp+loFqPM9y02i6ktUW31HIuYHgluONR3p4PI0SEo55J9nHJSOToAr9PK3YHEzutyZunO6zv28TqrvI4zb4SbDpr6GJ6PtT7RMa5rJ3LKSdwOORWN591Xxrc3kdPRdY+6e52rg5SO43Xhep4uP4pDTbfNKIOYxiZjkXRSwjvmrBvXezjSrK0TVN66VQgxC8D6brPWj7kCLEwC0VesWBFXdnm0GWyJ4ByP/9i8LxIxNHwISgL9VoaSEjPqX+cUScL+4dWLoSLdYqQUcWRMSL6RCO4ktJcA6lyeMFMvsZUJ4LyWi6xOgiU4dqRPIs+7WufsIOLamID8cd75JMSOu38P4ROG7HrJbmlaLrLL2HJMM3L/1m4IsrQs3jl3/bjew1QVXieg8u0jwHut/98LPOxYf6MQIlsIsYBRVICFKfC6w/inPVbnUB7z+UU1GHoxutZPeX6Y1qEwmggQdRmUxFRr5wuRQCokInmcLZ7gBYyzy70qPO5j4uzzBCq78L6gTnVZEtOvLSR8pWSqHTw8UPq3cC/aJoaTN04122qblDhSV7HsZmLqvMdWj7Pj1TmUOi+xQmGJZbRJ7OAYCWbnZUXECS1mrytPvXQ8d2k17ILSBQS0yZtpRko56gqwTz/9NMBSRwXYrwG/FUJ8AGgA3mmde8wVYGGKiD4RqrsKvJFSUl9fz7tD1/ODrPvQhKSmaIADnTmYgdEQe4PAEYplq+yuuHSvZHV2EJ7tCW3wJKSOG35LRHw8xHZ2AGqeMOf62F2lb50l2U84/5HgLfEcI7aI9QGK/I7tpq2Mm7DSs00R3tkmRXjNXjQDX0RsX6nKNKsdNJC6QPikaYurTsH6baWUbglvufLUuT+0OPXkDKkwEpHD4fCI4bFe/OpXv1L/eocALkly7TFVgIVJ8rp7MRESXQXe7Nq1i6ysLDZv3ozWmct/N/4QnyYpzemjNVJmliVyeovALU0dEtrrXY8LfHEQ3iXF0yG4jL3kru3E2uKV0s5tXoLjWU7YQaQLh/RN6IBzSnLHDhLiiC+cc1goKayI6UhPjXPKKagQWiW51f9G7Jx209TPapgOOhVZZ7IdTycuHc9ZY3F+OQsLy0b5oEaHUzUqDqZQdR+vjT4wMMDg4CALFixg9uzZAGyuWEOkPh9EiLI8g/bBYRBWVFIcwWLS3Kmyp03yERxGCSU9seNs4jjahZT2yxtHaqe0H480TwbvOYVntUMZ8qrxUt2cg/iKxDbpvYT2SHCnSm9Lb8Mh3Q1cGoCdOKPaYBNa2s9YFaiIg3UDHzrn/DE9KidGcsTBqVtdBibBGZcI45Xora2t7Nq1i9zcXJvkCr/c/O+EdA0DQU1pEE3zXMdJQIet7iKy2i8ByeMcTqkcTepYy7nlHKZDlwjD/Ni5287repxvXqdbfMisjPtoST6J9vU68eIcdE7nnCVhnfuY92DeD1b4sPNY2zvvPK8ee/5xH68zzuugc347n60yzRweePtHlVDky+WNlTVjfv8UTteiEzAFse4wdqJLKTl06BBNTU2ce+65SR/yP89/PxFDIBBUFPbGGGa2yHFC5ydWk935cnlJnsg+t789w0xegmMQI7eenNya4xxxhFYENsyPcHQWCUmS5JPY4251Our8nk7AG5+uRWPDdDVlXbxucX0c6VFJQdZxIxLe00nGEVu62+/tIBLVB3D94qYKwv8551x7jvTx4HSdvAGmSKKriRZHg3A4zKuvvooQgg0bNtg564mcIW+YuZpopIioFBTkSLICg7hYbb8EIuEwGB6iJCN53JCWh/SaY3+b4GmQO7beIp7hIbNHY3BL9ZjETevj3T8BoVQnoDoAL/kDWpRvvO8hvnTTH+OG18z78kh5L+ETSGa8v4VT8nskepxGptYT+7lNCHK0AO9csHJU714ypFNd5qyR6IkwWone29vLyy+/zLx581yBNyPNiPq7N3yBobCPqNSYVRYGLeK+O6dKJ4XbUaZelJFIrtYlkESJCO4sqxRHbvvYmLSOI7Xjhbdf/Kj5cQe9JFbLk6rr9sc6Xp0zmqBT8HSIqhP611seobykD39giEs37HKr6s7OUEl5B+HjtApnR+jRmlxS3Ety6Vl2am7KqScFNy9YM2FDaqdrBVg4BVX3pqYm9u3bx/r16+MmUUx1nu+tuR3dKjAwq6ofhGvQKO7l8DrQXA45r6T3vHheLcBL8ET2fGydRUQvmdR2J6EdnYFmTR8kvBI/AXGSqusO0sb8BdZ5PZ2Aq1OxzrP1ohfZuLyOMDohdP7hqr9YKrsluZOktTr38UbX2YFEybSsBOSP/Q4e08zxjwZsmTGLYDCY1ruXCqezjX7KjKMbhsG+ffvsSrCJHqgzaCYRlpbMIU+voo92EBozinvoay/FOZ7jlOReae71rntfOq+0ih0n3fsm6EScXvY4r7raD8d+4HhxHcc416cTPBN/CscBwj5QquNttzZ2vLlaWLb4ONdf9nciUhKyVGQRMHjjhl38/dW15m1oDm3JOVwmzRuWEitZhdhYuNPbDqg6ItL6o56TmXrquFdn++x2YoVTCK6bu4wsn58jR44wPDxMcXExpaWllJSUjCnv4nStLgOnyPDa0NCQXV5q3rx5SUMIndVqkuHhi/6JzX/8ArlZgtw8GMobItrrUKck8WRPJuG95HdKZ4kZleU4X5zEVwR3vIwurQFIGBjj6ggc2537eJFsvQNxj9XRUdjlm+2LCdfwWmlJNx/9wEPoQhKUscdnILjh6qd59pU1QCxvXAphD6UhsYNnzKE0aW9XQ2t4yI4z2MYiuP2MrA7C7kycJBfSes6COzdeSMDno7q6GsMw6Onpoauri7q6Ovx+P2VlZZSVlZGXl5dW2Go6w2veUaFTBdMeAtvZ2cmBAwdYsWKFXbtrLOdx4umLPsubn/kaAT/MKAvRFcwCPcu1j+tndZLPS/hUJPfuhyJ8jOBu6Q8C9TI6rm/vl0KiJ5DmaY+ve8fLwWa40oDtgBekTX5/VoTbP/wg/qwow1JYjyZWaJEAXPLml/jzn8+zx8cFEqnHJLywyBkLnnFId3UDDrLb4/SOjlQtO38v581LpD011HuXrCHgkL7O2VrBzLLs7Oykrq6OYDDIjBkz7O3JqrymI9FPVRt9SiLjEqnuKpS1vb19xKqsTiQqBZUI+Tk5vL38XH7f8Qo+DWbMHqSvzodQt+uRsKkkuncIDWna4l5V31z2ENwl5d3qe1Jiu7QOxzHObyfSkOZOiLheziHRRWwnKUDTDN5322MUlwwwJAN2wJlEsxz1ZkXVC9+8k6f+vAlD98VSUAUxCa/FAmvs0FiwYtRVuWfcZFfBNWA+R5/1nUCoSstEEEAAjX9eO3KATHZ2NrNnz2b27NkYhkF/fz+dnZ00NDTYnUJZWRkFBQX2O52OM+6sUt293nGv6h6NRtmzZw85OTls2rQp7amQR+PU+7cNl/Pww3uQOWGEEORVDzBUPwMptJiq7ZQOToKDi6Bup5lMKOHNjkK6SWw4yBynvscT2LWvl7xxywlU/nRgSWz7EJfO7uwEJEKTXHHz35i3sI1hAkjDlOReiS4R6Jpgyy1/5aH73mxWiLGIF8tDl6jplIRukl0R1kl2O1LO+ZsQ+99uv/OelR/AkuYfX3c+vlGMm2uaRlFREUVFRSxcuJBwOExXVxeNjY309/dTWFhIaWkp4XB4RIF01tvoTkk8MDDA7t27qampGbU9M9phup1bPs2Kh+4hJ9t8kQKzBtEbCl2O2jhJ6VXhXWROQnJDugiflOBGmuR2qewJnHjeGx1tHrQQjnPEyiqDW6Kf+5ZaVp5bR1j6MaQiuTXTiTSluk16KZi77CSFJQMMdBfYZZWFg9AgwRD2ZAsQ22aH1BuWxHeo5rbK7lLVVTstk0AzT1DgC3DLsjWjex4eZGVlMXPmTGbOnImUkoGBATo7Ozlx4gRCCIaGhigrK6OwsNAlpMYTMCOEqAf6MaP/o1LKTUKIUuA3QA1QD1wvpewey/mnhOhK9VEzs6xevXpM09ak8ro7IaXk2LFj3FN8Hp/rfYlAALRsiVERhBPWj5FMhU80ZmukQXJ7DNxBTu9yKnJ7O4fYDTluzvojPe99Kr677FlrWQhX56N+q3MuOMp5V9YSkT7C+C1iu4lu4LMUHzXXuMZl//g3tn/1MjMvG7PuuisF1XKWSY+6jjVxgz0Hm63bY0vxhPcppElyS83/3psvS/EQRgchBIWFhRQWFqLruv3enjhxgoMHD5KXl0dpaSn5+fkTMR3Tm6WUHY7lhIUix3LiKVHdDcNgeHiY5ubmEWdmSYV0vO4Qm7AhOzubqy98M4/9pZtnBw7j9wu0IgM9NIxsz4mp8Na3S2q6VPYEJDdAU0kVErMKi1d9N+IJLgyZnNx2W5ykdnQIiST6WKqaJJPoVmOrN7RwwU27CEs/Yem3iJxYonuJ7i/QWXLhMY78dVGM2CIm3W373KrOmkiyK0+8UFlstrvew3eBOXOq9RYvLS5jQ9Xkeb11XScrK4vi4mIqKyvtopCdnZ188pOf5OWXX+bb3/421113HRdeeGFaUzylwBbgTdb/9xErFDlqTHpkXLJQ1rEgHdU9GAzy8ssvU1FRwfLly9E0jR+95TrkUBaRsPmCauVR9MKwnWmV0B72qOEuh52T5AZxJDcjxdydg4pUswmrAlZ0GQsblRIMw14votJRVUZapoPEmlAMDMf5RvMxnMdb57WuWbm0g9e9fzcRfAzJLIIyi6CRxZC0Pob5CcpsgkYWQSPb8TH3WXxJPYG8YXcgjZSxSD9Pp+cyd6wacK4RiUSaj1LZ/WbPIKTgV1dcM7oXapTwOuNUUch58+bxwx/+kOrqai6++GIeeeSRsUzxLYEnhBCvCiH+0VrnKhQJVCY9OgUmVXXv7e2ltraWc845h8OHD4/7fJqmjTger4bqVq1aFTebxv6bPs45v/gWUupoPpCzIkRDgqw+R8fjkOIx9d2jmnsIbZLcuY8j/dRIIMFdwTXWBqftb7dFxo5xtA/cqr1z31QjwfapHIEygD0nWsnSLjZ/aDeGJhiWWaZJ7ZDkEmuCQixvu1TSXXngNUulF6z98B5e/fpGSzpbdrmlh0sdhOKLYZaGUBMjCglSlwifSGyJWDa7FKD7DLtkw5cvuMg1nDYZSBXrHg6Hueaaa7juuuvGcvrXSylbrKqvTwohDoy1nYkwaRK9sbHRDmWtqKgYlX2dDMnOIaXk+PHjHDlyhI0bNyadMue5q/8RY8iPHjXFgZwfJVIYcQqJOMSR3Dl+rmPHaLjsdQM7QcUlwR3SGyVFldTWPdLaKfUNiYgaaLphhc5KEmWxpZLkTknqPIemS0qXdLHuI/sxfBA0Apa09nz0LAYtKT6kW5LeyCJoBMxlPbafyJNUX1fniul3jVokSld1xCIgY9Lfy3ipgaFJyDN/tFJ/gC3nLE3nFRoXUmWvSSnTHkFKcGyL9d0GPIhZsz1ZochRY1Ik+pEjRxgYGHCFsiq1ezwlnxOp7oZhsHfvXoQQZtWZER50eWEhV+dV80iwCSNLInwQnRtFhAXZff5YuWJbnXR7012ESiTJVYfgOC5OgjtVVoiT3LHOIVngjHR/O48d4dk5T+EsBYUQzFjZy9Lbj2D4BENGtnV64egjrMAYwJCaNXaOPbRmzznu+S5YNkDBuk4GdpbF7HKsSRHVjKia+RzVkJwt1ZPcg7TG4Y08aw8DvrNu/Qh3PnEY6f0dz7wFQoh8QJNS9lv/vxX4ErFCkV/DXShy1JgUos+fPx9N01w3rsbS1fRIY4E3YGZ4eJhdu3alDJ114oa5C3hxXx+toX5kwEBogsjCKNRBTrfncSgCOUjrUu8TkNwpxWM2qKVaK5uV2LkAR8fgHX5zkNruGByNS6jbJkaCUSmT5Ju7qflAI7oQDOsBDGtPaTvcMJ1uNq80lwdeStCtoBkpzXnNVCdgSI3it3cw3JFDtCnfFeEmpJoJ1SK3CotVQ2qJCG+RPJpngN+8h/sufTu+np70H8Q4kEpQjYPsVcCD1rF+4JdSyseFEC+ToFDkWDApRM/KyoqTvBNVCVado6enh71797Js2TLKytKvBaZpGk/e8C7W/vSHyDyBzDIQGkTm6yAlee2BGFHBIc0dktlw2unJSW5LccNyRIGL4KpQQxy5ncQG20FlwyntHeti/3u2OXdWEx0gKby4l5nvbEUXGlEjy3UaRWRznbC967YHXsY87uYtKoJrViCN9b/UKLyxjd77ZmO0Z5v3LDBJjUTqZoSMHSSjnKNqQbXH8rCHcw2wYlbeOncBi0tKON7Xl/C3nmiMROTx1HSXUh4D1iZY30mSQpGjxaQNr8VdaIImcdB1nebmZhoaGli/fj15eXmjOofSCg5+6A6W3vtfZrRXFggNwvMlBKLkNzoei/R8UpHcQWohnZ52B3kVwW2Pv7TtcmcsuxAgVUegosu8avxo3y/dvJGC63oouawXXQh0w+8gtWqu5rhtEbfecEhu97JmqfJWByEFutTw39xBaHspojk3FnmoqXsXSJ/1DIWwHwnqvjGJHiowwBqmzifANy69lP7+/jHbxWNBMqIPDw+POM/adGNKAmZgYgpECiHo6elB13U2b948pnFKp/r/6vv/gY0/+RHSMJBZEnwQnmlgBCLMOOyPDS3jUdVtwuNS0W21HsAw7DBaW1/1EtyQVlBJrAOIhYti5rg71HvzGSi11sNwlySMe3COGzHI/0gv+auDRA1BlIDFKWHb5WaLY8vJiG6TGUAKm+BOshuGqcLrUqBv7Uf/uyR7V57bLlek93DI7mg0GCozQAWdheG5//N+87GkyCibKpzKZaRgiok+HokeDofZt28fAGvXrh2z48NJ9IKcHP7vJZfx4T89bjqKsgzwCaIl0L0qSskeDS3qeImkQ5o77W6Xum6q3kJ3kFE6verYBMd6HPYxAohKV+cQJ9V1h9ruJPeIkt3aWGSQ8/E+/LMiRAyNiDXoouLXsYkek8jS+gZz2AyphtawHXH2kJtll5uqfMxRp6t9pCBy/jDDC3Xyn8jDN+RDqvLR1jO0Y9ytr4jPYGg+4LP2CcGuD6lh5lOH6KdydRmYQtV9PETv7+9nz5491NTU0NzcPO5ZWZ0OvUuWLuGdh46xvf6wqUIGTEePzBV0bZSUvKqTNaTFyKzsa1tlJ84mt0mulpW6bw2rxUSVg+BOb7tXqoP5f9TB5iTSW9g7O3eVyBVRArcFIccgrPtQk13Y0txaUBLc+b+U8cQ3x8udpI4R35Dm83IS3HD8r5dIuq4bRjT7KNgZwN+nIaTP1d5wniRYLZBKGzaAEHxz5VJ2795t55JPFdFT2eAZia4uNIYCkRCLj1+zZg05OTk0NjaOqx1qWmcnvnLV23jpu000DA+ZeqIukQGTND0bBflHJAXNICzJJpRzzBvw4iS1ktx22qolN5UKrxxqDnVeKtVdj5HcZKAwTQGr/VIaMZXeXIlt55uXsiGFJHK9ju/iCLoQGLrftYvTNlfrlaqOIrN5IpWhi+oYFHkVgVUnECO2RX7DQXTDQf4q6LxER0YlDEuIWMaSPTOqNJ+5DmJY8J23vJk3Lq5heHiY7u5u9u/fz/DwMFlZWfT09DBjxoxJI71hGCMKmFO5jBRMseoeCoXS3l9KyZEjR+jr67Pj46WUEzK1k5PoKtjmK+ev5kNPv0LQkKa9bkjz6QgYXATDFVC2S6KpW5BOaW6uEsqm9pLcEQQjlH7qXcbdQYB1DiGQhm4Ph6Hrsf1GsssBfRYMfUiiVRpIQxB1/NzS/uMYL1edigWlqivpbT4vi8DSQXbr27TTsckvHZLclvwW0aUB0tCQhqm6yyzAL+2ZW9Wc9URBGxa8Z9lyLjpnIVJKsrOzmTlzJlVVVbS3t9PT08PJkyc5ePAg+fn5trQfz1CuF+mUej7riD5e1T0ajbJ7927y8/PZsGGDqwrseOFU3VWdOikl69ev54X161n/7z/CsDzx0gDpk+ADvRDaLoCiAzCj0Wl/K8lsxAXFuIfXYjLUGWdue9cMx75Ou1watjeaqOP5KYnvhRAQgIErBJE3C3wBHWlo6I4gyHiSx1zbLqXA44yDGLFRNrhhnsdLcDVKiJLohuo0LHIb1jkMgZquGl0gdGF/iwiIkGBjSSWfvOwiu/2GYSClJBqN0tPTQ1FREZWVlXYKaXd3N7W1tRiGQWlpKeXl5RQWFo7r/UmnMGRGdSf94bXBwUF27drFggULmDVr7BPWJ4OS6OFwmF27dlFeXk51dTVgdiS7PnMb677yA4xcMLIw7XW/tAsb9C6DgfmCquckgSguld1FcnCTXEqHI0/t65DuhkSzs/4EwjAsCa60j5jkd6v15rLKGBxaJui7KQAlEr8wMKLulzPWNYg4ZcDpXVdj2E6bHNzS3HbSSctJZy3bUltik1lKpY6LGMEloIhtWCq6IczS01GBCAnmaNn87NZr435DKSVHjx4lJyfHrglnGAa5ubn2jD6GYdDd3U1TU5NdQKKsrIzS0tJRJ1edztVl4BQbXuvo6ODgwYOsXr2aGTNmTEo7NE0jGAxy7NgxFi1aRFlZmR0IoXr8nZ+/jfVf/oEZopkFUhdIv0T6zP2MHDjxZo2cVknlqw5vOlh2pc0QN8kNpaqDLdmtPkRKs+aZTWxlJEtpTgusG+ZxGhA17HNomoZuGIQrBB3vykOfr+HzGQgDwiQiudsmj/0vYvtIN+FjyovDaeeU7LY0V8tWDrpap8gtLcltWOsMN8HRTYJrEdCigtww/NeNb6C3t5cZM2bYv4/SxHJycli0aBFCCJuEivBKa1MloYQQDAwM0N3dTWNjo10uqry8nPz8/JTS/nSevAFOkeE1Z/24TZs2jWrq2dEiGAzS1tbGhg0byM3NJVkiwiuf/SCbv/w/ptrpV2q8sFV5hGB4pqDhckHOScGs56M4h99sVZ7YuLsir5mhhWWLW8a2RXABSMvxIw0zHlRGdbMD0UBGFMlBaILB2dB5RQHhpX40YXYohnNIMEEEfLxJ7yCz+iNwqfLESXVznXk+B4mtfW0Jn4jcYEpyw0FwXZijFRGBLwLZEXjmM7fGlXQqKyujtbWVoqIiFixYEHdvmqbZv6dS8RXx8/PzycvLo7q6Gl3X6e7upr6+nsHBQYqKiigrK6OkpCRhfEYqGz0YDNqFJ09FTLuNrus6e/fuxe/3j6p+3Fhw/PhxOjo6qK6uJjc31yXFvdA0jf93/Rt596//hsgx1VHDB9IP0mfZ7hqgCYZn+qjfqpHdZlC6N0puq25laBm2Cm+7vEzPlTtTTXnwDGkle1ik1g00TWDoum1TY0iMLMHAimw635pHtEJD85vH6qr6YjJI15eFmIpuXwOHQgEuQpvrLLVbkd3VGTi2OWxxYdnydtaaYRIbRfCoMLP4whAIwwtfug2AqqoqqqqqkFLS29vL3r17MQyDaDSKEGJEiazeJa+0V2p4eXk55eXlCCHo6+uziZ+oFHTGRk8Cb5WZRMNrw8PD7Ny5kzlz5jB37ty0zzvasVPDMDhw4ADRaJQlS5Zw7Ngx/H4/FRUVCYv9qc4nJyeHZz/7Xt7w5fswskFkm5FeJtExCxz6rBLDQhCq9HGiyge6JLdVp7g2RP5J3V2JxrLZpTTMITMV060kumFKdGk23CSVIZFZgv4aPz0bchlcmQNZmFlgGBhRc4LJVEhEcNd2KRw7xobYzHXObcRIDTFJbrhJLlRHYI9OeKS4YUpyTcd0vEUk/rDgpbtvS/gbHjt2jHnz5jF37lxCoRAdHR0cPXqUYDBIcXEx5eXllJaWJiWkkvZ+v98l7VWJqMLCQubPn084HKa7u5ujR4/aEz/4fL4R1fuBgYHxlpGaVEyb6t7d3c2+ffvSqufuPc9oiB6JRNi5cydlZWXMmzcPKSUFBQV0dHSwd+9edF2nvLyciooKCgoKCIfD7N69m1mzZtlOume/8D7e8MWfmR7mAOAHw2cWU5A+TOnrw5xS1JLyQ7P8DM3yA5JAr0Hh0TCFR0IE+gyz2kpUkdy022N9ooEhDIwcH0OVWQRrshlYkk14ph/8AjvGVjdTPmMVHCwkexe96rpr2UN6Gb/eJanVPvY64QgJdkp7K/ZABRYZIhZkpAhuzf+mhU2Sv/C1eJJHo1E7S3HOnDmAWa55zpw5zJkzx56cob29naNHj5KdnW1L62Tx505pHwgEbMJLKQkEAlRWVlJRUWFL+8bGRgYHB+nv77elvfPcE1EBVghxGfCfmMbh/0gpvzauEzowZUR3jl83NjbS3Nycdj1373l0XU8rzl158BcuXEh5ebntdMvLy2PevHnMmzePSCRCR0cHdXV19Pf3E4lE4irUZmdn8+yX3s8bP/9T9CyQAUzC+yyVXjPH0KUmYmWMNWmSTggiRT66NubStTHXLG0swD+oIzWBf1AHBFErN0dmCYxsgQxY5zLrMJnnMkA4pzDxIpVQT3CIS3VX+3hJ75TmFpETkdqMC4jZ5Gofe/JFPUZ0NceciIBmSfLnvxFP8kgkwq5du5gzZ07SURjv5AzBYJCOjg72799PJBKxnW5FRUVJBUQi217XdQzDYMaMGZSUlNge++7ubg4ePEgkEqGkpIRIJDLuEFghhA/4HnAp0AS8LIR4REq5b8wndZ4/RWjfiBtHQiQSiYtAe/bZZykpKSEajbJy5coxFaHYsWMHS5cuTZm15iwrpfYdSQvo6Ojg8OHDzJ07l/7+fnp6eigoKKCiooLy8nJ7ePD1n/wRRo7loFMfS41XHzRpTyhgk1Utg7uuj8D0tjuTOoR5jLSJqzoNT6OFY1u6SCCtXdscdrpNaHCTGmx13Sa2bZpYHYFHXdcU+Q2HFI+CiErydcFTX09M8p07dzJv3jyqqqpGcZMxRKNRurq66OjooLe3l/z8fFvapxtQ09fXR21tLStWrHAJJuU3+P73v8/Pf/5z1q5dy/XXX8/NN988WgEmhBCvA+6SUr7NWvEZ6xpfHc2JkmHKJHooFCIYDDJ79mxqamrGHLyQTuBNQ0MDJ06cYP369QQCgRGdbmBqGK2trWzcuNH+8aWU9Pf309bWxvHjxwkEAlRUVPD0Pe/jzZ/6GSIXDN2sf+a22TEluzCXhSbtxA0pTDXfRWRrDNyUqsqdDnalVvXH2Xx3/KtrW7KeWSTdIYHanmBZOMgfU8uJSXhnlKBDlXdNE+1U1XVTXa/0Z/HI12+Na284HGbnzp0sWLCAioqKJHeVGn6/n8rKSrtq68DAAB0dHezatQuAsrKyEQNqBgcH2bt3L2vWrKGgoCBO2hcXF/PpT3+aZ599ls9+9rO8/PLLY62iNAdwxnc3AeeN5USJMCVE7+vrY8+ePWRnZyccEhkNRpqWyTAMDh48SDgcZv369TbBk5HcMAwOHTpENBplw4YNLokvhGDGjBnMmDGDxYsXMzQ0RHt7O3v37uU/37WBf/nf19CzJUbAdM4ZvhjhbZtdOEgvsGYvUQS3UjSVpLa2uyS+8HyjOgMv8RP+mxwJ7HGbyMJa54jTSSTVhZfsDmLHpDvuKZQd0zBrYcnayjK+/5n4oimhUIidO3eyePHiURUVSQVnjfYFCxYQDofp7Ozk+PHjDAwMUFRUZDv0/H4/g4OD7N69m9WrV9v2dyJP/nPPPcexY8dYvnw5F1544Zibl2DdmDVqLybV6w5mofv6+nrWrVvHrl27bDt5rEgm0ZUtV1xczJIlS+KCYLxQ00IVFRWxdOnSlG3Kzc112fXb5s7lvXc/zEDEwMgSCL+IqfCObyGw7XZFcOkguLCKLzjJD44OAGwNwPxXOGZTGc2TSwCl9buI7xwuc5I8Rm6nne7N0XcR3qmqK4JHJVpIctMbVnD7jRfFNUmNxCxdunRUTtqxICsri1mzZjFr1iwMw6C3t9f216jAqhUrVozoZHvttdf41Kc+xQsvvDAuzQNTgjuHnqqBlvGc0IlJtdH379/P4OAga9aswe/38+KLL7Jp06ZxFYg8ePAgpaWlrocaDAZdal4qkqtpmufPn8/MmTPH3BaAj969nR3N7RhZAsMvTJtdc5Dbh63Gx2x2bFXe/h/3stP+ls7b8PwvvetSQUnwZF54F/lFzBZ3Ehon6d2kVsuuj27a4iIKWkjy3duvYMOqeXFNCwaD7N69m2XLllFcXDyKm5pYBINBduzYQVVVFf39/YRCIUpKSigvL6ekpMSW6jt37uRDH/oQDzzwAIsWLRrPJYUQwg8cwiwd1Qy8DNwspdw73vuBSZToR44cQdM0W4WGWBjseCvBOlX3rq4u9u/fbzvdUpXc7e3ttYf1kpWFHg3+63PX8cifdvPN3/wdI1uR3SHd9Ri5peb+X8W32LOXODsAQIl32+mdgPCJ+zIvixPsJD3fJCAxnuVEarvh+PaSXZWxVjb5kMHj//UB8vLiIx+Vmrxy5cpJC39OB6qzWbNmjT0urqLo2tvbOXToEDt27KCtrY3f/e53PPzww+MlOQBSyqgQ4nbgj5jDaz+ZKJLDJEr0cDiM99zpesxHQl1dnT3lbVNTE01NTaxevZqsrKyUTreTJ09y/Phx1qxZM+H1vbr7BrjmjvswsjWMgOWVtxxzKNIrz7sgoWT3SniXre6000eS4Mm2jfBLxk8IYX2rEldOyU0Cgju2xz7SSk4x1fWsEPzpZx9KeP2BgQH27NnjsoWnA0NDQ+zatYsVK1Yk7WyklDz++OPcfffd9jv3i1/8goULF47n0uM1wlJfYLKIrut6XCTc7t27WbBgwbgiiBoaGhBCMDg4yPDwMCtWrEjpdJNSUldXR29vL6tXr56IObGS4tL3/F8ifiyyO1R5n2Vb+9yEdxHdQWyXHe5V0RORPt1XJZHqnkhtl/HLLkebjCe4Pf2UISEi0axAmGXleXz501ckjFrr6+uzvdrTGUKqSL58+fIRNb1Dhw7xnve8h1/84hesXr2a3t5e8vLyxjXVGGca0fft28esWbPG5WQ5fvw4jY2NVFVVsWDBgpT2uMp08vv9nHPOOVNSdujfvvY7/lrbYNrtAWFL9th4u4gRPQHhXcROtGwhLcme7BeUjlOlQXSXTe6Q5GrCCs3ACms1CS4iEl/I4L+/eA0zK/Nob2+nq6uL7OxsOzYhFAqxf/9+1qxZMy4tb7xQDsBUJK+rq+Omm27ivvvuY/369RPZhEkn+pSNo8P4C0QGg0Hq6+spLi5Oi+QqnLWyspJ58+KdP5OFL955Fe3tvdz84f9Fz9bMaDeHKq+m+TVVeJFSjYfEkjsh0Uehusep7OC2z5NIb/Nbejzr0g5p1cI62VHB739zu31e1bkHg0Ha29vZuXMng4ODzJkzh2g0Ou7RmLFCkXzZsmUjkryhoYGbb76ZH//4xxNN8inBpA+vOTEeoqvY+Llz59LX10coFBox+mhwcJA9e/awePFiysvLx3TN8aCioognt32Ud7//R5wcHMYIaEi/wPBJU43XHKR32u5C2FFyTnLHDbcxQUR3rnfa5MRs8tiMstKOgjMluDSjc6NWQcyIgRaRvPeaDbzr5tcnbEZeXh4FBQUIITjvvPMYGBigoaGB/v5+ioqKqKioGDExZSLhJPlIXv7m5mZuvPFG7r33XjZv3jzp7ZoMTJrqbhgGkUjEta6+vp6srCxXHHk6aG5uprGxkdWrV+Pz+WhqaqKjowMhhJ184HSudXZ2cvjwYVatWjWtzp1QKMTu3bvRRB6f+/KfTUed37LX/QJDOeqUZPeo8i4Jj0gQOktiYnvXJfoVvWPotqrurHuHLc2dc9HZzjZr2EzTpZl5psPvH/zoiM+kvb2duro61q1b5wpBVePYiVT80eZDpINQKGQ7h0cyJU+ePMl1113Ht7/9bS66KH7cf4Jw+troiYje2NiIlDJtNVpKyaFDhwgGg6xcuTLO6TY8PEx7ezvt7e1Eo1HKy8sxDIOuri7Wrl07qQUsUkFpFEuWLLGju+78xK/ZdagVI0szHXV2BpxwJcTEAmqEh+y47HYXRumMs/9PYJ/H1HXpHi6T0rbDhW4SXYtItHCUz3/q7bz+wmUjXrqtrY36+no7NHkkKBW/o6MDXdcpKyujoqJi3LXfIH2St7W18Y53vIOvf/3rXHLJhMyMlAynL9GllITDYde6lpYWQqFQWmGwqkBkQUEBCxcutFNTR7LHa2trGRgYIBAIUFZWRmVlJUVFRVNu+/X09Nhj+94RBsMwuPaq7xISIG11XoBPYAjARXgR740HVHTdqMju/SUdEtwVAWc4JLgqf2VLcmmPjWtRiQjrrFxYyTf++5aUz+TEiRM0NTWxbt26UXuoI5EInZ2ddHR0jFvFD4fD7NixgyVLloxYEaajo4N3vOMdfPnLX+ayyy4b1TXGgDOL6G1tbfT19bF48eIRjx0aGnJlLaVyuum6zp49eygoKGDRokUYhkFnZyft7e309fVRXFxsvxiT7XVva2ujrq6OtWvXjqhydnf1857rfoCR7cPI8oFmVrBBM9NeFeETSXZ7GdzOOc+yCwm87M553mLedMdkkkqq66Yk13QJuoEWNsjTNLb94WNpPZOWlhZOnDjB2rVrxz20OR4VX5E8VQx9d3c31157LZ///Oe56qqrxtXeNHFmEV2Rb9my5CqemiV1xYoVFBYWpiT58PAwu3fvprq6OqHtr4oStLW10d3dTUFBAZWVlZSVlU34eHpjYyNtbW2sWbMmbanV0tTJbbf8BLI0DL9mRdVpJtEFVnKMiFPnnUS3yzWnIdHVpEveOHazgKXyqsvYxJFKTY9KtKhOYcDH/3vsn9LuMBsbG2lvb2ft2rWT4mBLV8VPl+S9vb284x3v4BOf+ATXXntt0v0mGKcv0YG4CRt6enpobm5m5cqVCfdvaWnh+PHjrF69muzs7JSRbirYYtmyZWmNzTtTTzs7O8nKyrKdeeMp9q8mmxgeHmblypVj0hoGB0O894rvEkYiAz5r3N36CGHWkvPa7Ir4EOeRj28kHvtceobOZGzYTJd2ySsRNdAiOgvml/Gdn39wVPd0/Phxuru7WbNmzZTELyRT8QsKCtizZ49dgCQZ+vv7ue6667j99tu54YYbJr29DpxZRO/v76euro41a9a4L2IRpb+/3yZKKpK3tbXZUzWNNdhicHDQduYJIaioqKCysnJU4bGGYdj15RYvXjwh/oBvfuZ+/v70QZPwPp9VpkpzeOZFHNGll/Be2Oq6jH3bY+PSLmYZk+A6Imrw6bu38LpLEnfMI0FV7Fm1atW0TIKoVPzW1laam5spKChgzpw5SVX8wcFBrr/+em699Vbe/e53T3VzT2+ie+Pdh4aGOHDggCvgQKWL5uXlsWjRopSquppCqauri9WrV4839NBGKBSira3N5cGvrKwcseZ3JBJh9+7dVFRUTEpATjQa5RM3/4jjdZ1Iv2aq9LZa7yA9xFR6BTvf1fETOiW3Q4JjmLa30A1EVOfqmzfz/n8ZmwNKTaygwpOnc6bTSCTCjh07WLBgAfn5+UlV/OHhYW644QZuuukmPvCBD0xHU88soqvZUVTQgQpYqK6uZubMmWmFs+7fvx8hBMuWLZu0l0jVkWtra2NoaIjS0tI4D/7w8LA9o0xlZeWktMOLZx7byXe/+DuzyolPAy1GejfRkwykqznWFbkNAxE1KCnL5Y6vX4UIRBgcHKSkpISKigpXSmY6kFJy+PBhotEoy5cvn5ZINwVVhqqmpiYuT9yp4t9xxx0MDAxw4YUX8o1vfGPC4u1vvfVWHn30USorK6mtrQXMTMsbbriB+vp6ampq+O1vf6tMzjOL6Lqu8/LLL3P++efT29tLbW0ty5cvZ8aMGSlJrqRneXk58+bNm7KXSNd1urq67BGDoqIiCgsLaWpqYvny5dOaNw3w1KPP88AP/k5HS5BwxKxiGitdAwhrthifID8/izXn1vC+T19O+azEQ0tqGqP29nbbeVlRUUFZWdmI2pOUkoMHDwKkVchjMhGNRtmxYwfz588fsRMOh8O8613vYvHixQQCAQ4dOsTDDz88IW145plnKCgo4D3veY9N9E996lOUlpZy55138rWvfY3u7m7+/d//HU53onsLREopef7551mwYAH19fWsWbPGdoKNJDlUjvDChQunTHomgmEYNDY2UldXRyAQYMaMGZPmwU8HyqOtCntMNJTzUqm8qm6etx6+lJL9+/cTCAQmzE8xViiSpyooGYlEeP/7388FF1zAxz/+8Ulpc319PVdeeaVN9KVLl/L0008za9YsTpw4wZve9CbVOU76A5vytzMUCtHS0mLXaEvldOvu7ubAgQPTXpAAoLW1ldbWVl73uteRlZVle/BVaO9EePDTgZSSY8eOMTAwwLp16ybNhHHWzVu0aJGrbp6yc8vLyzl+/Dj5+fksXLhw2kmeTtXYaDTKP/zDP7Bx48ZJI3kitLa22iWrZ82aRVtb25RcF6aQ6CqoRUrJ2rVrU6rqYA63NTU1sX79+kmJd04XygHY3d3Nhg0bbOnpLB6p5nTbtWvXmD346bbl4MGDGIbBmjVrppRY3rp57e3tdh3ArKwsurq6Rm3XTxQUyaurq0ckua7rfPjDH2b58uV89rOfndaOaSoxqUR3Oq7U1EuDg4OcPHmSioqKpOqm8twODg6ycePGKclkSgYnsdauXZv0Jc7Ly6OmpoaamhpCoRDt7e32BALKg6+ytsYK51DedNvBmqbR1tZGTU0N1dXVrlJL6dr1EwVd1+1JHkaqAajrOnfccQfV1dXcddddU/78qqqqOHHihK26T6UZOukSXTndVL5vYWGhXdIpNzeXyspKysvL7RdCzXuWm5s75RLLC13Xqa2ttePt021LdnY21dXVVFdXu2aCGRwcHHMMvq7r7N69m5KSEmpqasZ4RxMDRazKykp72io1TZHTrm9oaMDn89l2/URrN6otO3fuZPbs2UlncgGzk/z4xz9OSUkJ99xzz7S8V1dffTX33Xcfd955J/fddx9btmyZsmtPqjOuubmZQ4cOsXbt2jinm5SSwcFBWltb6ejoICsri9LSUk6ePEl1dbU9x9Z0QQ0Fzp49e8LaksiDX1lZmTIGX5Wynj179qhTfCcaah60WbNmpdUWZde3t7dPeBaa6nBmzpw5YlsMw+DOO+8E4Lvf/e6UmBY33XQTTz/9NB0dHVRVVfFv//ZvbN26leuvv56GhgbmzZvHtm3bVGLN6e117+zsxO/32zNRjvTDtre3s2/fPgKBANnZ2fbsGtORaqq8/JNZtEJKacfgd3V1kZ+fb2s3TpMmFAqxa9cuampqpnXEAWJj03Pnzh1TmWw1ft3W1jau8XqIkbyqqmrEjtgwDL74xS/S39/P97///WkN4BkBpzfR77vvPhYuXMi6detGtLPVLJirV68mPz/flgLKKzlZjq1EUOWgp9LL7x3GysrKsmO0Dxw4wDnnnDNiSuVUQE2RNFEdjne8XnV06dj16ZJcSslXvvIVTpw4wY9//ONp9fWkwOlN9AcffJBf/vKXHDx4kIsvvpgtW7awefNml/quxoJVyWYvVGhqW1sbuq5TUVFBVVXVpBQTbG9vt+Pnp6JTSYZgMEhjYyNNTU3k5+cza9asKevoEkFlfi1atGhSNBw1J5pKNhrJrjcMg127dlFRUWH7B5Kd8+tf/zpHjhzhvvvum5Y4h1Hg9Ca6wtDQEI8//jjbt29n165dXHTRRVxxxRU8+uij3HjjjXHzniVDOBy2JX04HKa8vJyqqqoR49HTRVNTEydPnmTt2rVT4ikeCT09PRw4cMAuTa3ueSI9+OlChfqmKtQwkRgaGqKjo4P29nb7nisqKsjPz7ejI+fOnZv0eCkl//mf/8mOHTv45S9/Oe2/Zxo4M4juRCgU4qGHHuITn/gElZWVrF+/nmuvvZbXv/71o/pBvPHoZWVlVFVVjdrJ4ww+WbVq1bSrdx0dHRw9ejRh4Qqvjati8IuLiyeF9KrW+VTMg5YMznvu6OigsLCQhQsXJrXrpZTce++9/P3vf+e3v/3tpAcvTRDOPKIDfOELX2Dt2rVcddVVPPXUU9x///08++yznHvuuWzdupWLLrpoVD+Qrus26QcGBigtLaWqqirlEJZKklE136c7eEKVW3KOUiRDIg++GrueCIeTckimqnU+FTAMgz179lBcXExBQYHLrlfVZQKBAFJKfvzjH/PEE09w//33T2vNwFHizCR6IkSjUf72t7+xbds2/vrXv7J+/Xq2bt3KxRdfPKqoOC8BSkpKbKnnJICqSVdWVsb8+fMn45ZGhfHErafrwU8Xp8o8aOAmufN38tr1//u//0soFKKxsZEnnnhi0iIpa2pqKCwsxOfz4ff7eeWVVybitGcP0Z3QdZ3nnnuO7du385e//IUVK1awdetWLr300lE54ZRnt7W1ld7eXnvcOj8/nz179jBv3rxxz6Y6XijTYXBwcEKKNDgJoBJRVAx+OhKuv7+f2traaZ8HDczfr7a2lhkzZqQMEvre977Hb37zG4qLixkYGODPf/7zpEzxVFNTwyuvvDLRTsmzk+hOGIbByy+/zLZt23jyySdZvHgxV199NZdddtmo5nBTUq+pqYm2tjZKSkqorq6mrKxs2uxyZ3jtZOVvqxj89vZ2IDZUmajDPFXmQYMYyQsLC1NWDd62bRs/+clPeOyxxygoKGBgYGDSOqkM0acAhmGwc+dOtm/fzh/+8Afmzp3L1VdfzeWXX55WXnh3dzcHDx5k5cqVGIZhq315eXnjUnXHAhW3npuby6JFi6bEP6Bi8BN58Ht7ezlw4ABr166d1qFFMDtAFXqciuQPPfQQ9957L48++uiU+BIWLFhASUkJQghuu+02/vEf/3EiTpshejKol2H79u089thjlJeXs3XrVq644oqEVT5bW1vtKZO9udQDAwN2KG5OTo6t6k7WsIyKWy8tLZ02/0A0GrUdmH19fXZVmMrKyml1Skop2bt3L3l5eSmnIn7sscf49re/zWOPPTZlowItLS3Mnj2btrY2Lr30Uv7rv/6LCy+8cLynzRA9HSgVePv27fzud79jxowZXH311Vx11VVUVFTw4osvkpWVlZajyxl/7/f77VDciRqmOZXi1iE2fdX8+fPp7u62fRkT6cFPF1JK9u3bR05ODosWLRpx3yeeeIKvfvWr/P73vx+xfPNk4q677qKgoIBPfOIT4z1VhuijhUpxvf/++3nooYfo7e1l9uzZ3HvvvcyePXtU0spp32qaZtu3Y/XohkIhdu7cOaV15kZConnQJtqDny4UybOzs1OaMk899RR33XUXjz322JQ+x8HBQQzDoLCwkMHBQS699FK+8IUvTMRMLhmijxVSSm655RZKS0tZsGABDz/8MIZhcNVVV7F161aqq6tHRfrh4WE7FFdKaYfipmvPqnHpUyFuHUxTpqGhYcQpksbrwU8XoylF9be//Y3PfvazPPbYY1M+YnLs2DGuueYawDR9br75Zj73uc9NxKkzRB8PamtrWbVqFWC+TCdOnOD+++/nwQcfZGhoiCuuuIItW7aMugRSOBy2SR+NRm1Jn8xTrYasToVxaTADc5qbm1m3bt2opLSaFaW9vd3u7JJ58NOFlJIDBw7g9/tTkvz555/nE5/4BI8++ui0pzFPMDJEnyy0tbXx4IMP8sADD9DV1cXll1/O1q1bRx0hp0oqtba2EgqF7JdfxaKrCRdPhSErMGsEqJj+8ajiqrNrb2+38w5GG4OvSO7z+ViyZMmIx73yyivccccdPPLII5NSQ3+akSH6VKCzs5OHH36Y+++/n5MnT/K2t72Na665huXLl4/KGeX0ZAeDQXJzcxkcHGT9+vXTPmQFZvRdR0cHa9asmdDYAed9pxuDrxyoQoiUnevOnTv50Ic+xIMPPpjSE3+aIkP0qUZPTw+/+93veOCBB6irq+PSSy9l69atI9aLS4Tm5mbq6+spKCggGAxOegJKKtTX19PT0zPp86Cp+enb2tro7e21S2I7pzhW895LKVPWvqutreWDH/wg27dv55xzzpm0dk8zMkSfTvT39/PYY49x//33c/DgQS655BK2bNnCpk2bRiRLQ0ODLTn9fr/98re2ttrTOFdWVk5ZxVRndt5UD5f19vbagUkqCaW3txdIPdHD/v37ef/738+vf/1rVqxYMVXNng5kiH6qwJlTv3v3bi666CK2bNnC+eef75JUqeLW1TTOra2t9PT0uCaBmGgSOudBW7ly5bQHwgwMDHDgwAGCwSCFhYUjevAPHTrEe97zHn7xi1+wevXqaWjxlOL0JPq2bdu466672L9/Py+99BKbNm2yt331q1+1y/p897vf5W1ve9tYLjGtGB4e5sknn2T79u28+uqrXHDBBVx99dU89thj3HDDDWzevDktUimJ19raSldXlz13e3l5+bhtaKUe67o+7fOgqfYcOXKESCTC8uXL7eHKRB78uro6br75Zn72s5+5JuScSDz++OP80z/9E7qu88EPftAuHjlNOD2Jvn//fjRN47bbbuOb3/ymTfR9+/Zx00038dJLL9HS0sJb3vIWDh06NO3FHsaDcDjMk08+ycc+9jHy8vLYsGED11xzDRdeeOGooumklPT19dlqriqFPVL9+5HOdeDAAYQQ017/XbXn6NGjhEIhVqxYEdceVTlo165dfO5znyMSifClL32JW265ZVLarus655xzDk8++STV1dVs3ryZX/3qV9NpHkz6DzQpBtvy5ctZunRp3PqHH36YG2+8kezsbBYsWMDixYt56aWXJqMJU4asrCyOHDnCRz7yEV599VXe/e5384c//IE3vOEN3HbbbfzhD39geHg45XmEEBQVFbFkyRLOO+88Fi1aRDAY5NVXX2XHjh20tLQQiURSnkdFmPn9/lOC5GD6CJKRHMxnOGfOHNauXUthYSHve9/7eOyxx/jkJz85Ke156aWXWLx4MQsXLiQrK4sbb7xxwiZXPFUxpRXzmpubOf/88+3l6upqmpubp7IJk4I77rjDfoEvvvhiLr74YnRd59lnn+X+++/nrrvuYuXKlWzdupW3vOUtKQNMhBAUFBRQUFDAokWLGBwcpK2tjR07dtjx94lsW5URpxJCThWSDw0NpfQRnDx5khtuuIHvfOc7E5EkMiKam5tdNeeqq6t58cUXJ/Wa040xE/0tb3kLJ0+ejFt/9913J52BIpGZcCq8jONFonvw+XxceOGFXHjhhRiGwUsvvcT27dv56le/yuLFi9m6dStve9vb0sqbzs/PZ8GCBSxYsIChoSHa2trYvXs3QghX0s2ePXsoKiqa9plcFNTsNKtWrRrxd25ra+Od73wn3/jGNyad5HDmvocjYcxE/9Of/jTqY6qrq2lsbLSXm5qaTokMrsmGpmmcf/75nH/++XZO/bZt2/iP//gP5s2bZ+fUp5NPnZuby/z585k/f75dCru2tpb+/n572O5UQF1dHf39/SlJ3tHRwTvf+U7uvvtuLrnkkilp29n4Hk7q8Nqb3vQmlzNu79693HzzzbYz7pJLLuHw4cOntTNuPFA59du2beP3v/89FRUVbNmyhSuvvDLtxBc1mUFpaSmBQIDW1lYikYgrFHeqUV9fT19fX8px++7ubq699lo+//nPc9VVV01Z+6LRKOeccw5//vOfmTNnDps3b+aXv/wlK1eunLI2eHB6et0ffPBBPvrRj9Le3k5xcTHr1q3jj3/8I2Cq9j/5yU/w+/185zvf4e1vf/tYLnHGQXnKt2/fbldLufrqq7nyyiupqKhIKBXVVMHe3HZVCru1tZXh4WE7Dn0i5jtLhePHj9PT08Pq1atHJHlvby/veMc7+MQnPsG11147qW1KhN///vd87GMfQ9d1br311onKQhsrTk+iZzA+OHPqH374YbKzs7nqqqvYsmULM2fORAiR9jxo0WiUzs5OWltbxzWbazpoaGigu7s7Jcn7+/u57rrruP3227nhhhsmtA2nKTJETxd33XUXP/rRj6ioqADgnnvu4fLLL5/mVo0fUkoaGhrs9FqASy65hCeeeIKf/OQno8rkUqWwW1tb6e/vt0thqxpo40FDQwNdXV0pY+kHBwe5/vrrufXWW3n3u989rmueQcgQPV1MYFmfUxZSSnbv3s3VV1/N/PnziUQiXHnllWzZsoUFCxaMiqyJSmFXVVWNKf5eZcWlSvwZGhri+uuv55ZbbuHWW28d1TXOcEw60U/pmecycEMIwfPPP89Pf/pT3vzmN9s59f/yL/9CT08Pl19+OVu2bEkrp17TNMrKyigrK0NKSXd3N21tbRw6dIjCwkKqqqpcGWfJ0NTUlBbJh4eHueWWW3jnO9/J+9///jHdfwZjxxkl0X/2s58xY8YMNm3axLe+9a1pmy9sOtDZ2clDDz3EAw88QGtrqyunfrRz0XkzzqqqqhLG36sa+WvXrh2xQwiHw7zrXe/irW99Kx/96EfP+DHrMSCjujsxUpDO+eefT3l5OUII/vVf/5UTJ07wk5/8ZBpaOf3o6enhkUce4YEHHuD48eN2Tv1oc9HVvO2qZlxOTo5N+ra2NlpbW1OSPBKJ8L73vY/Xv/71fPzjH8+QPDEyRB8L6uvrufLKK6mtrZ3upkw7nDn1hw4dsnPqN27cOGpbXBWKbGlpIRqNsmjRIqqqqpIm70SjUT7wgQ+wfv16PvOZz2RInhynZ1LLdODEiRP2/w8++KBdFPJsR2FhITfeeCPbtm3j+eef5/Wvfz0/+MEPuOCCC/j0pz/Nc889h67raZ2roKCA3NxccnNz2bhxI7qus3PnTl599VUaGxsJhUL2vrqu86EPfYgVK1ZMGcnvuusu5syZw7p161i3bh2///3vJ/2apwvOGIn+7ne/m507dyKEoKamhh/84AfMmjVrupt1ysKZU//aa69xwQUXcM0113DBBRckTYs9ceIELS0trFu3zqWuO0thd3d389xzz9HY2Mj8+fO55557pkySn8YjLxnVPYPJRzgc5i9/+Qv3338/zz//POeddx5bt27ljW98o62Wnzx5kqamppQlotva2vinf/on9uzZw8yZM/ngBz84ZUNpGaKPcIEM0dPDKVaRZNIQjUZ55pln2LZtG3/729/YsGEDVVVV9Pf38/Wvf31EkhuGYT+X7373u/T09FBfX8+GDRumpO2n8chLhuinAk7BiiRTAl3X+epXv8oPfvADysvLWbZsGVu2bEmYU28YBl/4whcYGBjg+9///qQVoTxDR14yATOnApwVSQC7IsmZTnRV7FLNU/7SSy+xbds2vvrVr7JkyRK2bt3KW9/6VvLz8/nKV75CV1cXP/7xjye10my66dH/8A//wJVXXjlp7TjdkCF6GjgbK5IA+P1+l0R05tTv2LGDbdu28a1vfYtwOMw555zD9u3bpzXl+MSJE7YDNjPy4kaG6GngbKxIMhI0TWPjxo1s3LiRe+65h0cffZSLL7542usKfOpTn4obecnARIboaeBsrEiSLjRN4+qrr57uZgDw85//fLqbcMrijAmYmUxs3ryZw4cPU1dXRzgc5te//vUp83JnkEE6yEj0NOD3+/nv//5v3va2t9kVSaax7FAGGYwameG1DDKYfmSG185G1NTUUFhYiM/nw+/388orr0x3kzI4zZEh+imKp556ivLy8uluRgZnCDLOuAwyOAuQIfopCCEEb33rW9m4cSM//OEPp7s5GZwByBD9FMSzzz7La6+9xh/+8Ae+973v8cwzz0x3kyYV27ZtY+XKlWiaFuePUFNYLV261J4bIIPRI0P0UxAqGKeyspJrrrnmtJ9xNhVWrVrFAw88EDfv2r59+/j1r3/N3r17efzxx/nwhz+cdpGMDNzIEP0Uw+DgIP39/fb/TzzxxBkfs302TbM9Xch43U8xtLa2cs011wBmbvjNN9/MZZddNs2tmh6cqdNsTwcyRD/FsHDhQnbt2jXdzZhwZKbZnl6kiozL4AyEEOInwJVAm5RylbWuFPgNUAPUA9dLKbunuF1PA5+QUr5iLX8GQEr5VWv5j8BdUsrnp7JdZwIyNvrZiZ8BXnvgTuDPUsolwJ+t5enGI8CNQohsIcQCYAmQMdLHgAzRz0JIKZ8BujyrtwD3Wf/fB2ydqvYIIa4RQjQBrwMesyQ3Usq9wG+BfcDjwEeklBm3+xiQUd3PUgghaoBHHap7j5Sy2LG9W0p5WlRWzCA1MhI9gwzOAmSInoFCqxBiFoD13TbN7clgApEhegYKjwDvtf5/L/DwNLYlgwlGxkY/CyGE+BXwJqAcaAW+CDyE6fiaBzQA75RSeh12GZymyBA9gwzOAmRU9wwyOAuQIXoGGZwFyBA9gwzOAmSInkEGZwEyRM8gg7MAGaJnkMFZgAzRM8jgLECG6BlkcBbg/wOskg38yF1IEAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x, y = np.meshgrid(np.linspace(-10,10,201),np.linspace(-10,10,201))\n",
    "z = x**2 + y**2\n",
    "\n",
    "plotproj = plt.axes(projection='3d')\n",
    "plotproj.contour3D(x,y,z,200)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef679d02",
   "metadata": {},
   "source": [
    "### Revisit 1d harmonic oscillator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1d3e711b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pendulum geometry\n",
    "length = 2\n",
    "c = 9.81/length\n",
    "\n",
    "# Damping\n",
    "b = 0.1\n",
    "d = -1.0\n",
    "omega = 1.0\n",
    "\n",
    "def f_ODE(t,theta):\n",
    "    return theta[1], -b*theta[1] - c*np.sin(theta[0]) - d*np.sin(omega*t)\n",
    "\n",
    "theta00 = 2.0\n",
    "theta10 = 0.0\n",
    "\n",
    "solution_RK45 = integrate.solve_ivp(f_ODE, [0,100], [theta00, theta10], method = 'RK45', t_eval = np.linspace(0,100,201))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "87463039",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x24af7d30790>]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/YElEQVR4nO2deZgcVdX/v2cmM9kmYbISIGFLQoCACRBDAsi+BYIBX9HwiCD8NIKgor4gLo++LviqIK8CSogimwiCyiYBAmEVCDBAWJOQAAnZM9mTIdtkzu+P09eurqmqrqreZvl+nqef6q6+fftWdfX91jnn3nNFVUEIIYSEUVXpBhBCCGnbUCgIIYREQqEghBASCYWCEEJIJBQKQgghkXSpdANKQf/+/XXvvfeudDMIIaTd8Oqrr65W1QFB73VIodh7773R0NBQ6WYQQki7QUQWhb1H1xMhhJBIKBSEEEIioVAQQgiJhEJBCCEkEgoFIYSQSComFCIyRESeEpE5IvKOiHwzoIyIyHUiskBE3hSRQyvRVkII6cxUcnhsM4DvqOprItILwKsi8riqvuspMwHA8MzjcAA3ZraEEELKRMUsClVdrqqvZZ5vAjAHwB6+YpMA3K7GLAD1IrJbyRr1s58Bjz1WsuoJIaQ90iZiFCKyN4BDALzke2sPAIs9r5egtZi4OqaISIOINDQ2NqZryNVXUygIIcRHxYVCROoA/APAZaq60f92wEcCV1pS1WmqOkZVxwwYEDgLPT+9ewMb/U0ghJDOTUWFQkRqYCJxp6r+M6DIEgBDPK8HA1hWsgb17g1s2FCy6gkhpD1SyVFPAuBmAHNU9dqQYg8COC8z+mkcgA2qurxkjaJFQQghrajkqKcjAXwRwFsiMjuz7/sA9gQAVZ0KYDqA0wAsAPAxgAtK2iIKBSGEtKJiQqGq/0ZwDMJbRgFcUp4WwYRi6dKyfR0hhLQHKh7MblPssgstCkII8UGh8ELXEyGEtIJC4aV3b2DTJqClpdItIYSQNgOFwkvv3oAq0NRU6ZYQQkibgULhpXdv29L9RAgh/4FC4YVCQQghraBQeHFCwdnZhBDyHygUXmhREEJIKygUXigUhBDSCgqFFwoFIYS0gkLhhUJBCCGtoFB46dXLthQKQgj5DxQKL126AD17UigIIcQDhcIP8z0RQkgOFAo/FApCCMmBQuGn2EIxezbw6KPFq48QQspMpdfM/rOIrBKRt0PeP1ZENojI7MzjRyVvVLHXzb7qKuDSS4tXHyGElJlKLoUKALcCuAHA7RFlnlPVieVpDkwoVq4sXn1r19KVRQhp11TUolDVZwGsrWQbWlFs19O6dbbGBSGEtFPaQ4xivIi8ISKPiMjIkn9bsYVi/Xpg61agubl4dRJCSBlp60LxGoC9VHUUgOsB3B9WUESmiEiDiDQ0Njam/0YnFKrp6/Cybp1taVUQQtopbVooVHWjqm7OPJ8OoEZE+oeUnaaqY1R1zIABA9J/ae/ethRqMVa5a2nJBsYpFISQdkqbFgoRGSQiknk+FtbeNSX90u7dbbt1a+F1eS0TCgUhpJ1S0VFPInIXgGMB9BeRJQB+DKAGAFR1KoDPArhYRJoBbAEwWbVYPqEQampsu2NH4XWtX599TqEghLRTKioUqnpOnvdvgA2fLR+1tbbdvj28zOzZwC9+Adx8czaRYBAuPgFQKAgh7ZY27XqqCPksivnzgVNOAe69F5g7N7ouWhSEkA4AhcJPPqGYMgVYtcqeb9sWXRctCkJIB4BC4SefUCxdCuy+uz2nUBBCOgEUCj/5hGL79mxcIiqOAdD1RAjpEFAo/OQLZm/blhWKOBZFVZUtiEShIIS0UyqdFLDtkc+i8ApFHIuivt7mUlAoCCHtFAqFnziup9697Xkci6JPH6uLQkEIaadQKPwksSjyCYWzKLZupVAQQtotFAo/UTGKlhbLAhvX9eQsis2bKRSEkHYLhcJPlEXhhCGJRbHHHoAIhYIQ0m7hqCc/UULhhCFpjKJXLwoFIaTdQovCTxyhSDrqaft2CgUhpN1CofATFaNw+7p3B6qroy2KrVvt0acPsGULhYIQ0m6h68lPHIuia1cTlCiLwqXvqK/Pup5KnCGdEEJKAYXCT5xgdteu9oiyKFz6DicUzc35YxqEENIGoVD4iWNR1NbmFwq3Ql737tmYBt1PhJB2CIXCT1SMIonryQlNTQ2FghDSrqmoUIjIn0VklYi8HfK+iMh1IrJARN4UkUNL3qhiuZ5c2dra4gjF9u3Am2+m/zwhhKSk0hbFrQBOjXh/AoDhmccUADeWvEXV1bbN53oqt0Vxzz3AqFHAbbelr4MQQlJQUaFQ1WcBrI0oMgnA7WrMAlAvIruVtFEi1rnnG/WUz6IotlCszZymKVOAl15KXw8hhCSk0hZFPvYAsNjzeklmXytEZIqINIhIQ2NjY2HfGiYUlXQ9bdmS/e4//CF9PYQQkpC2LhQSsC9wMoKqTlPVMao6ZsCAAYV9a5hbKa3rqWdPe97UlL5NTij23JNBcUJIWWnrQrEEwBDP68EAlpX8W4vhevJaFPlWzXv9deDXv45u05Yt9p11dcDHH0eXJYSQItLWheJBAOdlRj+NA7BBVZeX/FuL4XryWhRdu+Z+3s8ddwDf/a6lMQ9jyxabk9GjB4WCEFJWKprrSUTuAnAsgP4isgTAjwHUAICqTgUwHcBpABYA+BjABWVpWD6LIqnryVkUYcLi0n1s3pzNTOvHKxQrVuQ/BkIIKRIVFQpVPSfP+wrgkjI1J0u+GEWxXU8u3cemTfGEghYFIaSMtHXXU2XI53pKa1GElXcWRVSQmkJBCKkQFIog4rieksQo3NyMfK6njRvD66NQEEIqBIUiiCihqKkBqqqSuZ4AK1+oRdGtmwlFIcNsCSEkIRSKIMLcStu3Zzv+JK6nfOW9MYowtm7NtSi4tgUhpExQKIKIsijcUNe4rqcumfECtbXB5ZubswIRN0YBZNOYE0JIiaFQBBFHKGprgZ077RHE9u3Z+AQQ7npy1gQQTyjcLG/GKQghZYJCEUTUqCdvzMHtC2LHjqzbCQh3PSUVCmdRUCgIIWWCQhFE1DwKr+vJ7Qtix46sqLg6g8q6QDZAoSCEtEkoFEHEdT0B4RaFcz05wlxPlRCKlSu5fjchJDYUiiCSuJ6iLIo4rqc4QqEaTygWLgT23x/46KPgegCgsREYPBiorwcuvji8HCGEZKBQBBF31JPbF4RXVIBw15OLUeyyS7hQ7NhhCQPdPAogWCjeeAOYNw949dXgegDLE9XcbMd4zz3h5QghJAOFIogoiyKu68lvUeRzPe25Z/jMbLcWRT6LYsMG2y6LyMTuJuvtu68JE+djEELyQKEIIiqYXQrXU20tMHBguEXh5kx4hSJodrazTpZHZGJ3nxs0yNrIWAUhJA8UiiCKFcyOO+qpvt6WSw0TirgWhROKOBbFbpmlx6PySxFCCCgUwcRxPSW1KKIm3PXpUzmhCPrOlhYLit92W3g9hJBOA4UiiCiLohSupz59bB2KtmJRrFyZPyhOCOk0VFQoRORUEZknIgtE5MqA948VkQ0iMjvz+FFZGlZbax29P9BbKtdTEouiWzdLC5I2mL15s22jhGLxYtuuWhVeTzl5+eXssRFCyk7FhEJEqgH8HsAEAAcCOEdEDgwo+pyqjs48flqWxjlLoLk5d38h8yiiXE8uRhEWXPYKhUj4mhTOolizJrxdcVxPS5bYtrExuI4k7NwZvRZ4PrZvBz71KeCqqwpvCyEkFZW0KMYCWKCqH6jqdgB3A5hUwfZkcR283/0UNI8ianis36KIcj316mWvgzpur1AA+YUCCF9Xu6kJqK4G+vWz11EWRaFCsXUrcPTRwDmRK95Gs3ixnbdnny2sLQDwxBPAY48VXg8hnYxKCsUeABZ7Xi/J7PMzXkTeEJFHRGRkWGUiMkVEGkSkobHQDi6OUDgRiJpw549R+Mu2tOQGs4FooejWzbZRQtGnjz0Pcz81NVkGWrc2d9D3FcP1pAp89avACy8Ac+akr2fRItu+9lr2PKTl4ouByy8vrA5COiGVFAoJ2Oef/fUagL1UdRSA6wHcH1aZqk5T1TGqOmbAgAGFtSwo/qBaHNeTN+7R1GRiscsuxbMoDsx47+IKRZRFsXp1tNvoJz8B/va34Pdeegm4/Xagrg5Yuza8jny4dCQ7dgANDenr+eADYMECexTiCiOkE1JJoVgCYIjn9WAAOb2bqm5U1c2Z59MB1IhI/5K3LMiiaG62Tj7JzGy/68nV4/BOpHNCEdRxe8sBwUKhagHfAw6w12GT7pqarPPu2dPiHVFCsXNnrjvLz403AtdcE/yeq+Ooo3LzWSXFWRSAWSdpmTHDtlu2RE9IJIS0opJC8QqA4SKyj4jUApgM4EFvAREZJGIr/4jIWFh715S8ZUFC4SyHJLme/K4ntz+ozqQWhX9mdlOTdezDhtmqevksCpHwkVaLF2e/K8z9pGpB89deCxaTNZmfafhwE7W0M8A/+sgC7yNGAM8/n64OICsUADB/fvp6COmEVEwoVLUZwKUAHgMwB8A9qvqOiFwkIhdlin0WwNsi8gaA6wBMVi1DcqIgoXAdfCGuJ395Zyl061a468l11n37WscaJhSbN2dXyevdu7VF0dxsnx01yl6HxXs2b7ayLS3BgebVq207bJht01oVixYBe+0FHHGEWRRpfv7mZmDmTODkk+11IULR0gI8/jhw553p6yCknVHReRSqOl1V91PVoap6VWbfVFWdmnl+g6qOVNVRqjpOVQvwPSQg390/kBWBJPMo/OWTCIVI9rujhGKXXYDdd89vUQDBk/yWL7fO8NBD7XWYUHjjDk891fr9NWvMxbXrrq3LJ2HRIkuYOHKk1Zkm5chrr9nnLrzQfocFC9K1RRU47jgTnHPPzYohIR0czswOIo7rqarKyiWZmQ2kdz25iXZAsFC4CWn19fYIm6DmFYpevVp3vC62cNhhtg1zPTnXUpcu4ULRr59ZOEA6i6Klxdqz116F1ePiHCNHAkOHprco1qwx62nsWHs9b166eghpZ1AogojjegLCJ9G5zxbT9eTcTkC0RVFfb2XDhpL6LYowoRg92rb5LIrjj7d1MNb4QkdOKNxw3TQWxapVdr68QpGmnpUrbbvrruYKSysUTnDOPtu2FArSSaBQBBHHogDC03IA8VxP7rPdutl3VldnxcNLqYQiKJjthGLffa2uMIvCddjHHmvbhQtz31+zBujfPysUaSwBNzR2zz0LE5wVK7KTDIcPB95/P90QWdeeo4+234tCQToJFIog4gpF167BQtHSYo98ricnCt6RVGEpPNxkO8A6+qgYRY8e0UJRV2fPgyyKlSvtu3bZBRgwIL9FMXy4bf0WxerVua6nNB28u4Mv1PW0cqUdS1WVWRRbtkTnwwrDCcW++1o9hQjFjBnALbd0roWjli+PXheetFkoFEEEdepBrqewtBxOYJK4nlyZMKHwWxTNzblC5hWK7t2DJ+S1tNj+qGD2xo1Wh4gtphQmFN7hr97X3vf79cvWVUhsoVCLYuXKbFDdtTeN+2nRIjv3/frZcN1ChOKKKyy4/pnPBFuRbYnGxsIEbetW4Fvfst9xypTC2/PAA8BNNxVeD4kNhSKIQi2KMFHxvuetM59QbN3aWiiAXDHYsMHq6dYt3PXk9vmD2d5OYMOG7KztAQOiXU89e9oIKyBXKJqbTbj69bO7+Pr6dB38smV2rPX1hVsUgwbZ88GDbZtm0t1HH1lnJ2JC8f77rRNHxmHHDuDdd20W/f33t921y5ubgSuvtBuGQoYD33QT8Nvf2rl/6KH0qVhUgcsuA848E7joIuC999K3yVtnW6WlxeZGtQEoFEFEBbPjCEWQRVGo6ymfULgstICV3bq19Z/ApRj3WhTOynA4iwLI73rq2zd7p+8VCteZu8SDffqk6+BXrcpaAi7FetoYhatn4MBs3UlxQ3UBE4odO4APP0xez7x59tnvfc/O4dNPJ6+jHHznO8CvfmXX7l13pa/n0UftfE2daq7PmTPT1fPqq8Dvfgecd5616brr0rcJAP7+d7sebrmlsHoczc2FZSHw8v77wMEH22CRNrBccSyhEJHTReQKEfmRe5S6YRUlyqKIs8aE+5x/hJS3HqAw1xMQLRTe+h1uNrdXKIBc95PXohg4MDzfk3Mtdeli3+sVCvfcCUXfvtEd/M9+Bhx5ZOv9jY3Zjh0wwUkqFKq5rqf6emtzGqH46COLlwDW8QHp3E9vvmnb0aOBY45pm0KxdStw663AF74AfO1r1rm7G40kbNsGPPOMzT059lizYh94IF2bnnzStr/6lWUkvuWWwiZyfvnLdkwXXgj8tMAVDNassTk2AweakLmRdmlYsAA4/HBL9//ss8All1Tc8skrFCIyFcDnAXwdlsjvbAB7lbhdlSXOhDtXLmzJVCD+PIq0QuFN47F+fdYScGX9Jr5fKILyS3ktivp6M3396UKArEUBmCBECUU+i+KRR2zWtT+usmqVWTWOvn2TdwwbN9o5dUJRVRXtUgtjyxb7jNeiANIJxRtv2PUwYoR1nh9+mJvTqhg88gjwxS+m69wBYPp0O3fnnw9MmmTnME2K9hdesHN30kl2fZ92GvDgg+lcKk89ZbnMBg0yF9THHwN//WvyegDLbNzSYqL92c8C//u/6XOAbdlia6a8/DIwebK5Er/5zXR1AcBvfmO/2yuvAN//PnDzzSa2FSSORXGEqp4HYJ2q/gTAeOQm8+t4xJ1HEbZkaiVcTxs2JBeKfBaFK5dUKNyM5TgWxc6d1nECrTvdYlgU3jkUjoEDkwuFGzbsLIp+/eyRxk/+5psWn6ipyQ4vTtsRbNlid53z5mU73+uvByZOBP7yF+BPf0pX75132jk77jhL7Ni3r8VTkvL442bBueM84ww7986qisuOHcBzz1l7ALPGhg5NJ14bNli7vvlNG9zwy19a/VdfnbwuwNxyc+YA994L3HEH8O1vm1i8+27yutavt6zL55wD7LefCUXXrumtsCIRRyhcb/OxiOwOYAeAfUrXpDZA3BhFmFCETc4Dgl1PxRCKzZuzw17D1tV2Hb4rl8+iyCcUTgj69ctNZ+FEo38m0W+URTF/frad3nUrVItjUbgFnFwwG0gnFN45HY5hw8yXnJQ33wQ+8Ql7ftBBhcUpfvUrc1/tvz8wZgxwwQXAN75hHfIRR9jdadik0DDWrwceftjujrt0sceECbmJFeMyYwYwfnz2WvvkJ207e3ayel55xa5DJxSAWSlPPx38H4ziuefMmjjhBHs9dKi52KZOTb5Yl6oJ80EH2TkHTCh69jSXalJuvdX+D1//ur3u2dOO+eGHk9dVROIIxb9EpB7A1bD1IRbCVqPruJTLoti2zfZXZX6GqFFP3nkUQULgnUiX1KJwQtHSYs/zWRQuc2xc15OzKIL8rK+/nn0+d272+caNdh69FkW+WEcQxbIovEN1HcOGJc8btXq1jeZyQlFVZQHL6dOTd3iA3WkedpilfG9stI7m8suBf/4T+OEPzc+d1D1zxx12HZ53XnbfoYfaOUvSkW7aZL+vt3MfOtSuq6RC4dLEOMsEMKHYtMnWPknCk0/af23cuOy+yy+3/0vSc/X883YsX/96NsVO//42Kuuee5Jfr3/8owmry7UGAKefbjdUFcx6HEcofq2q61X1H7DYxP4Afl7aZlWYUsyjCHM9xRlF5V1ZDwgXCmcppHU9NTVZZ57PonCZY6OEoqYm254+fcwtEuQvnz3byu61V65F4TrytuJ6+ugj6wjc8FrAhOKjj5KNSnnnHdsefHB2nwt+Tp+erE1Lltj5+9znrGOaO9fuvH/9axOgU0+1jvmf/4xfpyrwhz9YMNXbWY3MLC6ZZLXCl1+2mw/vQIXqahPJpELxwgvmrnNWKmACVFVlbqQkPPWUWVvem6+DDgIOOcTcdUm45Rb7v3zhC7n7zzrLjt0F4OOwYIG5qz7/+dz9p59u2wpaFXGE4kX3RFW3qeoG774OSaEWRVzX07ZtuRdroUKRz6IIGh4LZBMIum0+oXCdtVcoNm3KHrcbEeXusKJmZ7/+uv1JP/GJXIvCdeR+11NTUzJXysqV1pk46wYwoWhqCp6UGMayZfY5r/gPG2Yd6wcfxK/HxTr23ju7b8IESw2fNJ7ghGXiRNvW1Zn7ySFi1oYTpzg8/bT9DhdfnLvfrZyYxO/+wgvWhsMPz90/erQJRZKRPLNn5woXYDcOY8YkE4q1ay0m5rVyHF/8oq2i6L0Oo1C1ddhPPDH7X3GMHWv/ryRte+gh2zoXlmOffSyIn/RGooiECkVm0aDDAHQXkUNE5NDM41gAPcrVwIpQzmB2PqFoabH6vELhLkrX0e3caaKQ1PXkBMHN6nYuKL/ryW8J+F1LbuuEwAmFIyzfk6oJxSGHmI/9vfeyE9ici8NvUQTVE8WKFSY21dXZfa7OJG6U5cutM/fiZnkncT8tWWLbPfbI7uvSxWIL06cDX/mKBVkfeyx/R/qvf2U7kTBGjrRRVXFEUdViHn37mpXiZfBgE6IkQvHii/b97jpzjB5t15o/P1gYq1aZULtElV5OOsksl7BsyX6efTabLt7P5Ml2U3HHHfHq+vBDsyiD6urSxVyKcX5Hx0MPmSDvu2/r9447zs5nhSbgRVkUpwC4BrZE6bUAfpN5fAvA90vftAoSJRTezr8YQpHP9RQURPdbFG6bVCi6drWyTij8FoVzHcWxKICsgLg8T44wi2LFCis7erR1dtu3ZyewBbme0szO9s7KdqSZdBckFG5RpqRCUV/f+g70K1+x8/7AA+anPvXU6EluO3fa3IbTTstabkGMHGkdVRyX0Z13Wsf2ox/lDp4A7DsOPDC+ULS0WMd2xBGt33OLYsV1P7lRcYcc0vq9k06ycxGU6j6Ihga7afBaXo7ddrMOOa6rzn3n8ccHv3/yyRbbinN9rF9vQXa/NeEYP95u2NKMpCoCoUKhqrep6nEAvqSqx3kek1Q1gdMzHBE5VUTmicgCEbky4H0Rkesy778pIocG1VN0qqvtj+EXii5dsoFnIJnrqbraHkldT0HzN9yf2AmEfzRTlFBUV+e2yzsiKcyiyCcUzm/shMI7ExrITgT03/W5+MEee2Tvip3Z7+72vT7pNPme/G0BiicUffvasSURiqVLc+Mcjr33tvO3apUd3267ZV0RQSxaZL9/UOfpxcUW8rmfVq82S2b8eODSS4PLJBGKuXOt8xs/vvV7Bx9s/yMnAPlw5ZzAeBk/3q7TuC6e2bPtWvP+77xMnGhtj2PtPPmk3YTsv3/w+yedZNs4bZsxw6zpMKFwgfcXK+P1jxOjeF5EbhaRRwBARA4Ukf9X6BeLSDWA3wOYAOBAAOeIyIG+YhMADM88pgC4sdDvjdk4EwF/MNvbwQLJLAqgdfA7juspSCiqq+21Xyhcxx41PLauLvcOtL6+tespX4zCCYvruL0Whaq5CbyuFf/oKodXcNwENicUq1ZZO7zHHWVRTJ5sifb8uPxMXpIKxc6dVtYvFCLJRz4tWZJ7bvz1AXZNHH+8dURhbgs3Asa5v8IYNsyuu7ffji534432e9x0U66bzsuBB5pgxrHo/v1v2wYJRY8eNkfgtdfy1wNY5z5kSK6V6qitteHBcYXi9deDXViOCRNs+8gj0fWomkVx3HHhFt3QofZbu3MRxRNP2P/EH8/x1tW/PzBrVv66SkAcobgFtq51Jvsb3gNwWRG+eyyABar6gapuhw25neQrMwnA7WrMAlAvIrv5KyoJfhEIEop8o57ylY/jegoSCsD+bK4D9wtFlEXhd3l4LQp3x5/PonAxCzc23isUGzda+aRC0aeP1eM63VWrct1Orpz3c46WFpsMdvXVucMbt241q8UvFC5AHlcoVq82sfC7sIDkQhFmUfg5/nhrX9gdfFyh6NLF7nijLIpt24Df/946Se9oLD8uoB3HjfXQQzaSbb/9gt//5CdthFYc//3s2dGd+0kn2fnIN7vdxTqirLD99rO4Tz6hWLDARNM7XNePC+THGb47c6YJXpcu4XWNG9emLYr+qnoPgBYAUNVmAMWIqOwBYLHn9ZLMvqRlSkMcocjnevJbFH4hSOt6AnIXL/KPZooSCmdtOKIsipoae4QJhavLCcXq1dYRArlCETSxD8gKlBOA4cOznW5jY+6IJyDc9bRkiZ2n7t0tjbVzW7nAsV8oeva0R1yhcKkd/BYFYEKxcGG8IbI7dpgrLMyi8OImg4Ul0Js/36zDIPHyM3JktFDcfbcJ6mWXRdcTd+TT5s12h3/mmeF322PH2rlw10sYW7aYlZlPKID8VoVzYUXVJWKCOXNmdPp31/kHWUxeDj/cRsVFDZxYuNDKnHhidF3jx9u5KFbiwQTEEYomEekHQAFARMYBiDnEIJKgK8h/exGnjBUUmSIiDSLS0Jh0dmUQ/rv/crmedu7MHdkQRyj8MQpXp18ovGtROPwWhUi2HsDKBwlFjx5ZF0WPHvada9ZkFwRy6cfdcXfr1jpG4Tp8JwDeZUqDLAoX6/D/UZy4XHGFtfWVV+x10GxqR5K5FFFCcdRRZtH84x92Xm67zSZaeWeqO1assDvoOBbFXnvZ6Jewcfjz59v5igpkO0aOtLvtsEWDbrnFrA7X4Ua1qXv3/EIxY4Zdt5P8DgIPbob2yy9H1/XGG/Z/iLICDjzQ2nbffdF1ueB5UKzDy4QJ9l95/vnwMi+/bP8NJ55hOFdS1HG6mwF3cxCGi1MknWBYBOIIxbcBPAhgqIg8D+B2WILAQlmC3JxRgwEsS1EGAKCq01R1jKqOGeC/E01DMYQijesJaG11eN9zBAmFE4GqKivvFwp/KhCgtUXRq1duwL5nz9bDYzdtyhUTIDvpLsiiAIJX01u71s6Js0yGD7d5Blu2BFsU1dWtM9UCWXH5zGds+9Zbti2HUJx0krkr/u//bOLcl75kE6a+HvAXCRoaG8Wxx4b7t+fPD3fr+HEB7aD5Adu2md873+gpwK6LAw7IHxi//34T/099KrzMqFH2/8knFO74g0ZPOUQssd/jj2ev5SBefz081uHlmGPsWosaSfXSSzZyKiye4zjsMDtvUZ37zJlmGeYTnU9+0uqqQJwir1Co6msAjgFwBICvAhipqgkzegXyCoDhIrKPiNQCmAwTJC8PAjgvM/ppHIANqpoyxWNCunaNJxSqrcc2R7me/Ck8/BaF2+8t433PESUUQPDiRR9/3Nr11KeP/blaWnITAjrCLArnTnIMGmQdoRMKr0UBBAvFunXmdnIdlBtu+v77rRMCOoLWyJg/387jQQdZJ+wCt0Gzqb3t9S6H2tJiKS9GjcpOYHNECUVVlY0WamiwO9qf/tQE65lnWvvf3bmJY1EAdjxr1rS2Ttw6GPniEw4nKEEJDF95xa6xqE7dS76RT5s2WXxi4sRwfztgv9eoUfmF4vnn7brwj1zzc/bZdl4e9HchHvLFOhy9epkIhAnFtm1WV1jg2Utdnf2OYUKhalbj8cfnF+pevayuCsQpYq1HAQs8jwJwKGx00nl5yuclE+u4FBYonwPgHlV9R0QuEpGLMsWmA/gAwAIAfwTwtUK/Nzb+tSbChAJobVVEuZ78SQGLYVH4YxRA8HKoYRaFqv3BvQkBHXV1wULhtygOPtiS3S1bZnX6BSnMonDxCSArFP/4h4mve+2lf//WHeeCBTYqpKrK2uG1KAYNan3uACv//vvZtTYaGoCrrjJRePjh3DjI8uV2TGFDKs87zwTsmGMs2+eJJ9pn/EMsk1oUbtilP6vuwoV2fuIKhTs3QULx3HO2PeqoeHUdeKAdh/+3dFx3nd14BFlUfsaOtfMeNolM1SyKOG0bO9ashXvvDX5/wwazqILmTwRx3HEmYkEJMd94w/qDsWPj1XX44dl0Jn7eecfiQ/ncTo7x4010guoqIXmFQkTugE28OwrAJzOPmGc7GlWdrqr7qepQVb0qs2+qqk7NPFdVvSTz/sGq2lCM741FHIsiaBKd93U+11MpLYoePeJbFIDd3SexKPxCMXq0XfCvvRbcEe6yS36hcB3fDTdYx+Zy3HgJsyicqBx8sI3K2bEjeGisY8QIOz+u83Z3ab/5jW29d7pBcyi81NWZFTNjhrkiXG4jv9toyRL7vb3HHIV/yLAj7ognR9euNk8jTCgOOCB3vkoUUSOfNmwArrnG5gK4GEQUY8faDUqYK+u99+ymII5QOPfTjBnB82xeesmEJ1/w2XHssTavIShO4a6NOBYFYG6z9euDz1nc+IRj3Dg7z4Ws156COBbFGABHqurXVPXrmcc3St2wilMKi8IvPkHBbCCeUPTsmdz1FGZRACYUQRZFkFAExSicST9rVrBQhFkUTqgAe963b7ZzCIo1+YWipcUsA9dpHnywnWM3XDJMKJw7xv3hZs0yl9CkSdbpeF0F+YQCMDeZuz5c2gp/J7N0qZ2bOAFowAK0XbsWLhSAHa+/c9m509oY1+0ERI98+u1vrUP8yU/i1XXKKSasYRlbndAGrX4YxHnn2W8fVN+LLwbnnQrjyCPNdRaU/n3WLLse4lqGRx9t26A1R2bONIvPrXOSj6iJdzNmANOmpctCnIc4QvE2gBhj8DoYcWMUQOVdT01NVrf3+5LEKAD7gxdiUbi02aqt4xNAdIzCi7MMzjyzdR2ACcXq1Vn/vxsa6zrNgw6y7Vtv5bcogOxd9qxZdrfZu7d1hkmFwkt1tdXltygWL47fubh69tuvtVDMm2cCH9cKAKye997LjZu89Zb9JkmEYp997Fr0C8W6dcC111rW1HyzxR2DBlkQ/bbbsjm+vDz1lAWe3W+Vj9Gj7RG0BvaLL9q14b++w6irM4vHPzzZTbQ75pj4gr/PPnYT4heK5mbbF9eaAOx37NMnOKB9zTW2CFNUbCglUUkBHxKRBwH0B/CuiDwmIg+6R9Fb0tYoxKLYvt0uIv+ICK/rySX7K5bryT/stZgWhX/UU1Awu0+f7F1RmEURNDzWLxSuww8Tiv797bw50XF3105gDjggO2Jl69ZwoRg0yDqDefNs2OrChdm7tXHjzL2gao+kQgGYRfTOO7lDeefPT2YFABan8AvF3Ll2nHE7KsA626am3OU+3V1p1IgiP126WF1+obj2WvtN/ud/4tcF2HrVK1YAjz6au3/xYhtmfPbZyY7zggvM/eldQa+lJXsjkIQJEyzY713/et48i8OF5XcKQsSExT/AoaHBzlkSoaiqsuvzqady61q61ETt3HOTna+4Xxvx3jWwZIA9AJwJ4BewpIDXAsgzBKEDEHd4LBBsUdTUtP7BvFaKf71s9773Pe/zqJnZmzfnFwrV0loUQNb9FOV6chf39u1Wj18oLrwQ+MEP7C4sCOeOcu4nt8KcE4pu3cz95O4qw0x6kexdtrs7c0Jx+OE22uj9963D2rYtNy14HJyL49VXbbtunQ3HjXt37BgxwiZjea+JuXPD8wuFETTy6aWX7HyGneswRo7MTQmybh3wu99Zp+4sy7icfrq57X7+81yL89e/tmvlylYp4KL5whfsf3qjJ9vPnDl2bScViokTrQ3eWdpuXksSoQBMKFauzD3/991nNzVJ6/rsZ20Ah9fq/etfTRC/+MVkdcUkKingM6r6NICazPNnPPu6h32uw1Co68lfFsi1UvzLoHqfxxWKbdvMz+xdtMjhH/W0Y4ddSEEzswHreD/+OF6MIkwo3ESmMKFobs4et39WtuP4463TCMMvFIsW2V2u9zvvvjs7Ft65ooIYMcLuEF94wX5Lt96B6+RffDHeOP4gnPvFreDn4gNJhWL//bNxGMA6vOXL0wuFN04xa5aJY9I70KOOMreeO7Y//MHiVj/8YbJ6ADvv119vgnr00bZC3zXXWAbdL30pvu/e0a+fxSpuuSVrCbjOPalQjBpl19W//pXd9+STZqUGpQKPwh+n2LnTsvVOmJDMhQiYUHTvbi47wMTs9tvtt0xqscYkyvV0sYi8BWBEJnOre3wIoBjzKNo2SVxPQaOe/PEJV2dSiyJIUIBsh79lSzzXkxMNv+vJTbBryAwo84/xr6uzNrghjNu32yNIKFxnGnSx+vM9+RMLxsUJhRsiu2iRtdnr5hsxwtJlr1tngcIw9tvPPv/HP9qQVvdbjBxp3zN9uglFr17J75T79bPhmn6hSNrBu/LO/ZS2nsGD7bd3d7Tr1lld3uVA4zJ5sl3Lt9xi19jvfmcdXtJz5Pjc52zuxYoV5jq6/HKbqPbjH6er74or7Br97W/tv3T11Sb+cScoOkTMqnjsMauvpcVcPieckFxc99vPrNJbb83GOZYuzV1uNi69e9tcnbvvtv/m3XebhXf++cnriklU1OOvAB4B8L8AvPbfJlVNuB5lOyTJ8Ngw11NQnX6LohDXE2ACECQU/uGxTij8FkVVlVkRLk+Ov+PwLpLUq1c2XhEkFCefbG6SIFeGVyh23bV1qvK4uLsvr0URFIcQaR1H8TNihP1pm5qyw2IBE51Pf9p85LvvbgKYbwZuEIcckisUXbokd/Psv7/9Rq+/bp2DG2KZVCiqquwzrj1Jh3h66dvXYkh33mnnr7ExuYvIz6mnmu9/zhz7XxVyZzx8uN11X3+9TUxcvBi4+eZ0vvuJEy2j7j332H9w7dr8OZmCELFzdNFF9l+79Vb734WlFc/Hl75k5/+MM8wyPPJI4MtfTldXDKJcTxtUdaGqnqOqizyPji8SQOHDY4NcT153UDFcT0C4UPgtCvfcb1EAdle/YYNduH7XiD+DrD9zrBeR8I7Qb1GkFYog11NS94TDDfX89rdbrxJ31lnmTpk3L9moIC+HHGKfb2oyi2Do0OAbiCh69jQXiAs8z51rdSQVHMCE/LnnLB41a5b9XnHmOwRx4YX2G95wgy2bmvYceamqMmuuGO6T3/zGYmZ/+5u1LU3nDtgQ3nHjgEsuAb72NbN0zj47XV0XXGBW5n/9ly1Kdf754ZM483HCCZY2pqHB+oa77irJaCdH3JnZnY9CYhRhrievvz/I9eTq9wuFf8EkIFcogmIGYa4nv0UBZOMULpeMv81AViCiLIooXOyjUKHo2dOObfVqi3ksXZpeKEaNsru7n/609XsnnJA9xrizlv0ccojdcb/1lglGUivAccQR1rE3N5tQDBuWXHAAmyPS3GwjjJ55xoQy7nBRPyeeaJ3nXXdZjKIEI20KYsgQW/b0738H/vKX9O2rqbFjrKqya/fWW9Ode8D+37/4hVmnV19tsZi0iFi234ULbXTdkCH5PlEQpZOg9k6hFkWYULgAdBLXU1AKijgWxdat1lGJ5LcogGB/td+icBlIkwpFWIwiqVAA2Ul3S5ea3zitUADhd5rdutkY//vvT3/X7QLaDQ02NNafQyouRx5p60W89ZYJRb7kcWGMHWsjjP7nf0y4rroqXT2AdXY33JD+8+Wgqsru3gtl771t6On69dGDI+Jw7rk2MqtYwrrLLq0HoJQAWhRhFGpRBLmevJ1uEtdTWqEAst8Tx6II8leHuZ7SCoWbS7F2rf1Z0lzk/fubULiFasLmShTKNddYQDvonMVhyBATwqlT7RopxKIAbC3nBQuSj5xyVFebWM2bZ1bJd76Trp7OyKGHJh/GGkZbs75iQKEIw41Q8o77L4ZFAVinm2QeRZBQeIPMUULhBCKORVEOofC6nurrW7u64uBmZzuhKMSiiGLIkGSTofyI2JyQDz6w12nvRvfc04Zp/vzn5oYsZKz85z5n7br++uDripAAKBRh1NaaSDQ3m6uopSVZUsC4FkWhrqeNG61MmFA4gYiyKM44w0ZjBOVWShLMjsIvFGvWpHM7AVnXU9R6E22Fb3/bRvM8/XR6F5ZI1qq47rrWgfcknHKKzS849dT0dZBOB2MUYbjOefv27F1vEtdTlBXgtSi85cKC2VFC4Ub/+O/wvfMsvNsgi+LTn7ZHEK7eQi2Krl3t+JxQpEmL4fC6ngYODD6mtkR9vU0ALITvftfEohhDIIuxsBfpVFAowvB22m4MfRKhCLpzz2dRVFVZnUmEwq3SVohFEUWxgtlAbmLAZcviLSITxIABJlivvNK2rYlicthh9iCkAtD1FIbXoghbX6IYwWz/OGrvpDwgvkWRTyjcNq1Q+IfHpgnwevM9uZTbaTjtNKtr9uzSxScIIf+BQhGG16IohVCETaQrlVCEpfDIR1CMokePdDOVnVBs3Gj1BaUjj8Po0TbDeOLE4gx/JIREUhHXk4j0BfA3AHsDWAjgc6q6LqDcQgCbAOwE0KyqRVlZLxZei8ItO1gpiyLo7r2mxjpr53oKi1G4Dn7LFguKJh3pUlNjD69QJA1kO9wqd26t6rQWBWBJ2R56KP3nCSGxqZRFcSWAmao6HMBM5OaS8nOcqo4uq0gAuSOaopY2dWW8lMP1JGJiEGZRuM7cuYo+/tisiTRjuOvrsxPkgla3S1LP6tXmdgLSWxSEkLJSKaGYBCCTIxe3wda7aFsU4nrati2e66mqqnV+lrhCAZhQzJlj7fCvl+AfjrplS/qJY97lR8NSjMdhxAjLXurmPxRiURBCykalhGJXVV0OAJntwJByCmCGiLwqIlOiKhSRKSLSICINjd41ldNSjmB2UEKwpEIBWDpm/925XyicRZGGYgnFwQfbuXLrENOiIKRdULIYhYg8geC1tn+QoJojVXWZiAwE8LiIzFXVZ4MKquo0ANMAYMyYMRpUJhFei8LNzi50HoXX39/UFHyHn0QoBg2ydQ+++93W77nOvFgWxVtv2fOgVeni4mYmz5hhbqi07SGElJWSCYWqhub1FZGVIrKbqi4Xkd0ArAqpY1lmu0pE7gMwFkCgUBQdr0Xh8AtFVZU94loUQDaDbNB60e574wrFQw/Ze0HphauqLE5RLIvCLRS0eXP6uQv7758NwI8cma4OQkjZqZTr6UEAbjmm8wE84C8gIj1FpJd7DuBkAG/7y5WMODEKwCyENEKxZo1ZA36SCEW/ftFuIO8Et6D1suMyYIAJ286dhQWzu3XLrjVAtxMh7YZKCcUvAZwkIvMBnJR5DRHZXUSmZ8rsCuDfIvIGgJcBPKyqj5athXFiFG6f1+rYudMe5bAo8uEVii1b0lsU/fub+231anukdT0BFqcAGMgmpB1RkXkUqroGQKu0nBlX02mZ5x8AGFXmpmXxDn11Q0rjWBTueRyLIiibqFcoVAsTil69sik3Pv442IKJg8sNNGeOCU4hs6EPOgi4915aFIS0I5jrKQxvJtckQhFlfQC5FkU+11Nzs4lFsSyKQlxPgC3AAxQmFLQoCGl3MIVHGHEm3AHphGL9ervTz+d6CkvzERd/jKKQYDZQHKEYP97WeQha+4IQ0iahRRGGN5gdlmYcaC0UrnOPEorFi+15PouimEJRTIvCP7kvCYMGZdeRIIS0CygUYXiD2WFpxoFwiyKsc+/ZMzvUtJxCUYhF0b+/bd9/3+p0S6cSQjoFdD2FkWR4rH9t7bCyQG5OpnK5nlpaCrMoamqya1szrTchnQ4KRRhxYxS1tcliFN7OOsyicGt1F0MoVC2hn2phK8E59xOFgpBOB4UiDJewrxTBbEeQUHTrZp36jh3FEQrA1kgGCkuZ4YSikPgEIaRdwhhFFM4N5FJkBKXKKEQoglxPbtbz5s2FC4VLNb5ihW1pURBCUkChiMLNuq6psedBazmkFYra2tZrSAC5WV/bokVBoSCk00HXUxS1tdlgdljHn1Yo+vYNFp5SCIVb/8G9ToMb+UShIKTTQaGIwgWW8wmFd9RTnHkUQHg6Decu2rSpeELx6qu2HTYsXT2ATZKrqrIlSAkhnQq6nqJwFkVtbXjHHzbqKWoeBRAuFKWwKBoarJMfOjRdPQBwwQXAYYdlLQtCSKeBQhFFXIsirespCK9QtLRk25EGV9fChWYJpK0HsPjGuHHpP08IabdQKKJwwexSCEUc15MbZVXoqCcA2G+/dHUQQjo9jFFE4YbHllMoglxPQWtrx6Fr16zIUCgIISmhUERRSosizPXk5lEUI5gNZK0KCgUhJCUVEQoROVtE3hGRFhEZE1HuVBGZJyILROTKcrYRQGksikGDgCOOAI46Kvj96moTk40bLR05UNiwVvdZCgUhJCWVilG8DeAzAG4KKyAi1QB+D1sqdQmAV0TkQVV9tzxNRDyLwr8Uaj6h6NoVeP756O91yfy2b7dMrTU1iZueUxdAoSCEpKZSS6HOAQAJmnCWZSyABZklUSEidwOYBKC8QlFsiyIOvXub62nz5sKHo/bubeI0ZEhh9RBCOi1tOUaxB4DFntdLMvsCEZEpItIgIg2NjY3FaYF3eGxYnCBo4SKXUDAtvXqZRbF6deFCMWiQrVNd1ZZ/akJIW6ZkFoWIPAFgUMBbP1DVB+JUEbBPwwqr6jQA0wBgzJgxoeUSETeY7da2FokuGxfnempqAgYPLqyu667LBsUJISQFJRMKVT2xwCqWAPD6SwYDWFZgncmIG8wGTCxcOo9iCMUHHwBr1wKjRxdW1667FvZ5Qkinpy37I14BMFxE9hGRWgCTATxY1hbEtSiArPupGELRqxewYQPQ2JjN2koIIRWiUsNjzxKRJQDGA3hYRB7L7N9dRKYDgKo2A7gUwGMA5gC4R1XfKWtDu3a1taaXLQuPFXhXwnPbYlgUK1aYNcPcSoSQClOpUU/3AbgvYP8yAKd5Xk8HML2MTculttbWmgaA8eODy5TCoujdOxtXoFAQQipMW3Y9VR5vh3/EEcFlSiUUDgoFIaTCUCiicENihwwJn4dQqhiFgzEKQkiFoVBE4Tr8I48MLxMkFIXkZgJoURBC2hQUiihchx/mdgJaC4Vb6KgQKBSEkDYEhSIKl947jlAUc9STcz1VVwO77FJYXYQQUiBcuCiKs86yVeYOPTS8jBMFr+vJpQpPi7Mo+vVj6g1CSMWhUETRvz/w1a9GlynlqCcGsgkhbQDerhaKi2O4eQ/FdD0xPkEIaQNQKArFrVjX1GTbYloUFApCSBuAQlEoLh6xebNtiyEUPXpYbIJCQQhpA1AoCiVIKAqdRyECfOUrwBlnFFYPIYQUAQazC8UvFMWYRwEAU6cWXgchhBQBWhSFUgrXEyGEtCEoFIVSW2tDZCkUhJAOCoWiGNTVUSgIIR0WCkUxcEKxc6fN5KZQEEI6EBSKYuCEwuV7olAQQjoQlVoK9WwReUdEWkRkTES5hSLylojMFpGGcrYxERQKQkgHplLDY98G8BkAN8Uoe5yqri5xewrDLxSFzqMghJA2RKXWzJ4DACJSia8vPnV1wOLFtCgIIR2Sth6jUAAzRORVEZkSVVBEpohIg4g0NDY2lql5GZxF4RIDUigIIR2IklkUIvIEgEEBb/1AVR+IWc2RqrpMRAYCeFxE5qrqs0EFVXUagGkAMGbMGE3V6LQwRkEI6cCUTChU9cQi1LEss10lIvcBGAsgUCgqCoWCENKBabOuJxHpKSK93HMAJ8OC4G2PujpLM751q72mUBBCOhCVGh57logsATAewMMi8lhm/+4iMj1TbFcA/xaRNwC8DOBhVX20Eu3NS10doAqsXJl9TQghHYRKjXq6D8B9AfuXATgt8/wDAKPK3LR0OGGYP9+2u+1WubYQQkiRabOup3aFE4oFC2w7KCiGTwgh7RMKRTHwWhTdu2eXMiWEkA4AhaIYeC2K3XazFeoIIaSDQKEoBk4oPvqIbidCSIeDQlEMeva0bUsLA9mEkA4HhaIYeIfDUigIIR0MCkUx8AoFXU+EkA4GhaIY0KIghHRgKBTFoHv37EgnCgUhpINBoSgGVVXZgDZdT4SQDgaFolg49xMtCkJIB4NCUSzq6syyGDCg0i0hhJCiQqEoFnV1wMCBQHV1pVtCCCFFhUJRLOrq6HYihHRIKpJmvEPy3/8N7NxZ6VYQQkjRoVAUi0mTKt0CQggpCZVa4e5qEZkrIm+KyH0iUh9S7lQRmSciC0TkyjI3kxBCCCoXo3gcwEGq+gkA7wH4nr+AiFQD+D2ACQAOBHCOiBxY1lYSQgipjFCo6gxVbc68nAVgcECxsQAWqOoHqrodwN0A6N8hhJAy0xZGPV0I4JGA/XsAWOx5vSSzLxARmSIiDSLS0NjYWOQmEkJI56VkwWwReQJAUD6LH6jqA5kyPwDQDODOoCoC9mnY96nqNADTAGDMmDGh5QghhCSjZEKhqidGvS8i5wOYCOAEVQ3q2JcAGOJ5PRjAsuK1kBBCSBwqNerpVADfBfBpVf04pNgrAIaLyD4iUgtgMoAHy9VGQgghRqViFDcA6AXgcRGZLSJTAUBEdheR6QCQCXZfCuAxAHMA3KOq71SovYQQ0mmRYK9P+0ZEGgEsSvnx/gBWF7E57QEec8ensx0vwGNOyl6qGpjVtEMKRSGISIOqjql0O8oJj7nj09mOF+AxF5O2MDyWEEJIG4ZCQQghJBIKRWumVboBFYDH3PHpbMcL8JiLBmMUhBBCIqFFQQghJBIKBSGEkEgoFBk6w9oXIjJERJ4SkTki8o6IfDOzv6+IPC4i8zPbPpVua7ERkWoReV1E/pV53aGPWUTqReTvmXVf5ojI+E5wzN/KXNdvi8hdItKtox2ziPxZRFaJyNuefaHHKCLfy/Rp80TklLTfS6FAp1r7ohnAd1T1AADjAFySOc4rAcxU1eEAZmZedzS+CZvh7+jox/w7AI+q6v4ARsGOvcMes4jsAeAbAMao6kEAqmFpfzraMd8K4FTfvsBjzPy3JwMYmfnMHzJ9XWIoFEanWPtCVZer6muZ55tgnccesGO9LVPsNgBnVqSBJUJEBgM4HcCfPLs77DGLSG8ARwO4GQBUdbuqrkcHPuYMXQB0F5EuAHrAkoh2qGNW1WcBrPXtDjvGSQDuVtVtqvohgAWwvi4xFAoj0doXHQER2RvAIQBeArCrqi4HTEwADKxg00rBbwFcAaDFs68jH/O+ABoB3JJxt/1JRHqiAx+zqi4FcA2AjwAsB7BBVWegAx+zh7BjLFq/RqEwEq190d4RkToA/wBwmapurHR7SomITASwSlVfrXRbykgXAIcCuFFVDwHQhPbvcokk45efBGAfALsD6Cki51a2VRWnaP0ahcLoNGtfiEgNTCTuVNV/ZnavFJHdMu/vBmBVpdpXAo4E8GkRWQhzKR4vIn9Bxz7mJQCWqOpLmdd/hwlHRz7mEwF8qKqNqroDwD8BHIGOfcyOsGMsWr9GoTA6xdoXIiIwv/UcVb3W89aDAM7PPD8fwAPlblupUNXvqepgVd0b9rs+qarnomMf8woAi0VkRGbXCQDeRQc+ZpjLaZyI9Mhc5yfAYnAd+ZgdYcf4IIDJItJVRPYBMBzAy2m+gDOzM4jIaTBfdjWAP6vqVZVtUfERkaMAPAfgLWT99d+HxSnuAbAn7A93tqr6A2btHhE5FsB/q+pEEemHDnzMIjIaFryvBfABgAtgN4Yd+Zh/AuDzsNF9rwP4MoA6dKBjFpG7ABwLSye+EsCPAdyPkGPMLDd9IeycXKaqj6T6XgoFIYSQKOh6IoQQEgmFghBCSCQUCkIIIZFQKAghhERCoSCEEBIJhYKQMpDJ5vq1SreDkDRQKAgpD/UAKBSkXUKhIKQ8/BLAUBGZLSJXV7oxhCSBE+4IKQOZbL3/yqyVQEi7ghYFIYSQSCgUhBBCIqFQEFIeNgHoVelGEJIGCgUhZUBV1wB4XkTeZjCbtDcYzCaEEBIJLQpCCCGRUCgIIYREQqEghBASCYWCEEJIJBQKQgghkVAoCCGEREKhIIQQEsn/B1B9MOZaEKFQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel('t')\n",
    "plt.ylabel('theta')\n",
    "\n",
    "plt.plot(solution_RK45.t, solution_RK45.y[0], 'red')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "665666cb",
   "metadata": {},
   "source": [
    "**The normal harmonic oscillator is (mathematically) already a two-dimensional example:**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2dde00d2",
   "metadata": {},
   "source": [
    "Second-order differential equation: \\\\( y''(t) = f\\left(t,y(t),y'(t)\\right)\\\\)\n",
    "\n",
    "Introduce: \\\\( z_0(t) = y(t)\\\\) and \\\\( z_1(t) = y'(t)\\\\)\n",
    "\n",
    "\\\\( \\begin{pmatrix}z_0'(t)\\\\z_1'(t)\\end{pmatrix}=\\begin{pmatrix}z_1(t)\\\\f\\left(t,z_0(t),z_1(t)\\right)\\end{pmatrix}\\\\)\n",
    "\n",
    "Therefore, we can describe the second-order differential equation by a set of two first-order differential equations. We can solve both with our Euler method\n",
    "\n",
    "\\\\( z_0^{(n+1)} = z_0^{(n)} + z_1^{(n)} h \\\\)\n",
    "\n",
    "\\\\( z_1^{(n+1)} = z_1^{(n)} + f\\left(t,z_0^{(n)},z_1^{(n)}\\right) h \\\\)\n",
    "\n",
    "Or, going back to our initial nomenclature:\n",
    "\n",
    "\\\\( y_{n+1} = y_{n} + y'_{n} h \\\\)\n",
    "\n",
    "\\\\( y'_{n+1} = y'_{n} + f\\left(t,y_{n},y'_{n}\\right) h \\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "412f37a4",
   "metadata": {},
   "source": [
    "### Ball in a bowl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bcfb25f",
   "metadata": {},
   "source": [
    "**Now we have 2 harmonic oscillators. Therefore, mathematically speaking, we have a four-dimensional system**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1574872",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "520c1cfa",
   "metadata": {},
   "source": [
    "### Linear trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aae0fc93",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a5232a1f",
   "metadata": {},
   "source": [
    "### Elliptical trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f066f6a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f36c2779",
   "metadata": {},
   "source": [
    "### External force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a02579d6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
