02_KeplerProblem.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Keplerjev problem\n",
    "\n",
    "<img src=\"imgs/solarSystem.png\" style=\"width:600px; margin:0px;\">\n",
    "\n",
    "**Enačba gibanja**\n",
    "\n",
    "$ m \\,\\mathbf{a} = \\mathbf{F_g} = - G \\frac{m M}{r^2} \\frac{\\mathbf{r}}{r}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dt = 5e-06\n"
     ]
    }
   ],
   "source": [
    "k = 1 # G*M\n",
    "\n",
    "# Začetni pogoji:\n",
    "t = [0] # [s]\n",
    "r = [[1, 0]] # začetna lega\n",
    "v = [[0, .4]] # začetna hitrost\n",
    "\n",
    "# Račun s končnimi diferencami:\n",
    "t_final = 5\n",
    "N = 1000000\n",
    "\n",
    "dt = t_final/N # velikost časovnega koraka [s]\n",
    "print(f'dt = {dt}')\n",
    "\n",
    "def distance(p1,p2):\n",
    "    \"\"\" Izračuna razdaljo med točkama p1 in p2 \"\"\"\n",
    "    \n",
    "    return np.sqrt( (p1[0]-p2[0])**2 + (p1[1]-p2[1])**2 )\n",
    "\n",
    "def iterationStep(t1, r1, v1, dt):\n",
    "    \"\"\" Naredi korak iteracije \"\"\"\n",
    "    \n",
    "    x1, y1 = r1\n",
    "    vx1, vy1 = v1\n",
    "    \n",
    "    R = distance([0,0], r1)\n",
    "    a = [-k*x1/R**3, -k*y1/R**3]\n",
    "    \n",
    "    vx2, vy2 = vx1 + a[0]*dt, vy1 + a[1]*dt\n",
    "    x2, y2 = x1 + vx1*dt, y1 + vy1*dt\n",
    "    \n",
    "    return t1+dt, [x2, y2], [vx2, vy2]\n",
    "\n",
    "for i in range(N):\n",
    "    \n",
    "    iterationStep\n",
    "    \n",
    "    t2, r2, v2 = iterationStep(t[-1], r[-1], v[-1], dt)\n",
    "\n",
    "    t.append(t2)\n",
    "    r.append(r2)\n",
    "    v.append(v2)\n",
    "\n",
    "r = np.array(r)\n",
    "v = np.array(v)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.14176367374566184,\n",
       " 1.062142347169399,\n",
       " -0.3261189760950101,\n",
       " 0.3253084539082648)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAHwCAYAAABOjq0vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUVf7H8fdJIyGNHjqhtwChN4EAUgRBxb42VMSOa1kLixXburZ1dX+2VVxFgUVBV0FUJID0IiK9lxBKCBDS25zfHxPGAIEESHKTyef1PPM4t8y93zmJ+XDbOcZai4iIiHgnH6cLEBERkZKjoBcREfFiCnoREREvpqAXERHxYgp6ERERL6agFxER8WIKehEvYoyJNMZYY4zfeX7+BmPMD8VdV3EzxsQaY8acYdkFtYGIt1HQizjAGLPLGHNxvunrjDFHjTH9nKzLWjvZWjvYyRpEpHgp6EUcZoy5BXgHGG6tne9gHWXiCLis1CHiLRT0Ig4yxowFXgOGWGsX55vfwxiz2BhzzBjzmzEmJt+yWGPMS8aY5caYJGPM18aYamfYfrgx5t/GmP3GmH3GmOeNMb55y0YbYxYZY94wxhwBnsmb90u+z7cyxvxojDlijNlsjLkmX30HTmwrb94Vxpi1Z6hjkjHm3bxtJRtj5htjGuVbbo0x9xpjtgJb8+b1MsasyPuOK4wxvU7ZbNNiboNjxpgdefsdbYzZa4w5lPcPMZFyS0Ev4py7gYnAQGvtyhMzjTH1gO+A54FqwCPAl8aYmvk+ezNwG1AXyAHeOsM+Pslb3gzoCAwG8l/b7g7sAGoBL+T/oDEmGPgR+Dxv+fXAv4wxba21S4FUYEC+j/wpb90zuSHv+9YA1gCTT1l+eV49bfJC+7u871UdeB34zhhTvYTaYG3efj4HpgBd89a/EXjbGBNylu8lUrZZa/XSS69SfgG7gOPA14DPKcseAz49Zd4c4Ja897HAy/mWtQGyAF8gErCAHxABZAJB+da9HpiX9340sOeU/YwGfsl7fy2w8JTl7wFP571/Hvgo730o7uBvdIbvOwmYkm86BMgFGuRNW2BAvuU3ActP2cYSYHQJtMHWfMva5X02It+8RCDa6d8ZvfQ635euhYk45y7gSeBDY8zt1toTI0w1Aq42xozIt64/MC/f9N5873fnLa9xyvYb5c3fb4w5Mc/nlM/u5cwaAd2NMcfyzfMDPs17/zmw2BhzNzAKWG2t3X2W7Xn2Za1NybtcUDff/Py11M37XvntBuqdofYLaYOD+d6n59V36jwd0Uu5paAXcc4hYCAwH/gX7lP54A6hT621d5zlsw3yvW8IZAOHT5m/F/fRbA1rbc4ZtnO24Sv3AvOttYMK/KC1G4wxu4FLKPy0/Uk1550KrwbEn6GWeNwhnV9D4PuCtseFtYGIV9M1ehEHWWvjcV/nHmqMeSNv9mfACGPMEGOMrzEm0BgTY4ypn++jNxpj2hhjKgPPAdOttbmnbHs/8APwmjEmzBjjY4xpeg6P8H0LtDDG3GSM8c97dTXGtM63zufAOKAv8N9CtjfMGHORMSYA97X6ZdbaM51RmJW37z8ZY/yMMdfiPj3/bSm3gUi5p6AXcVhe2A0ArjLGvJQ3fRkwHkjAfVT6F07+//VT3Ne9DwCBuMO2IDcDAcAG4CgwHahTxLqScd+4dh3uI+wDwN+ASvlW+wKIAX621h4uZJOfA08DR4DOuG/OO9O+E4FLgYdxXyN/FLj0lH2UeBuIeAPzx2VBESkPjDGxwGfW2g9LYNu3ATdaawcUuvK5bXcSEGetnVCc2xWRwumIXkTyawvsdLoIESk+uhlPRAAwxswEmgNXO12LiBQfnboXERHxYjp1LyIi4sUU9CIiIl7M667R16hRw0ZGRjpdRrFKTU0lODjY6TLKPLVT0aidikbtVHRqq6IpyXZatWrVYWttzYKWeV3QR0ZGsnLlysJXLEdiY2OJiYlxuowyT+1UNGqnolE7FZ3aqmhKsp3yeqkskE7di4iIeDEFvYiIiBdT0IuIiHgxr7tGX5Ds7Gzi4uLIyMhwupTzEh4ezsaNG50uo8yraO0UGBhI/fr18ff3d7oUESnDKkTQx8XFERoaSmRkJPnGpC43kpOTCQ0NdbqMMq8itZO1lsTEROLi4mjcuLHT5YhIGVYhTt1nZGRQvXr1chnyIgUxxlC9evVye5ZKREpPhQh6QCEvXke/0yJSFBXi1L0UzaJFi2jYsCHHjx/HWktUVJTTJYmIyAWqMEf05UVISEiB85966il++uknAN58803S0tLOabujR49m+vTpZ12nRYsW/PnPf2b8+PE0a9YMgHfffZf//Oc/57Svc/Xiiy8Wab3IyEgOHz5cYnXMnDmT5557DoAhQ4YQHR3tedWtW5fu3bsXy37yf49evXqddd2srCz69u1LTk5OsexbRCoeHdGXEdZazjaS4IkAAnfQ33jjjVSuXLlYa6hZsyZjxowhPDycwMBAAO66665i3UdBXnzxRcaPH1/i+ynMK6+8wjfffAPAnDlzPPNTU1Pp3Lkzzz//fLHvc/HixWddHhAQwMCBA5k6dSo33HBDse9fRLyfjuhLyeuvv05UVBRRUVG8+eabAOzatYvWrVtzzz330KlTJ/bu3QvAww8/TKdOnRg4cCAJCQnAH0fkb731FvHx8fTv35/+/fsDcPfdd9OlSxfatm3L008/XWgtzz33HF27diUqKoqxY8dirSU+Pp7o6GieeOIJ7rnnHnx9fdm9ezfPPPMMr776KgDbt29n6NChdO7cmT59+rBp0ybP/B49etC1a1eeeuopz1mJlJQUBg4cSKdOnWjXrh1ff/31abU8/vjjpKenEx0d7Qmyzz77jG7duhEdHc2dd95Jbm7uaZ8raJ3c3FxGjx5NVFQU7dq144033gDggw8+oGvXrnTo0IErr7yywLMhW7ZsoVKlStSoUeO0ZQ888ADDhg1j0KBBZ22H0aNHc/fdd9O/f3+aNGnC/Pnzue2222jdujWjR48u8GdRlLa6/PLLmTx5coGfFxEpTIU7on/2f+vZEH+8WLfZpm4YT49oe8blq1at4uOPP2bZsmVYa+nevTv9+vWjatWqbN68mY8//ph//etfgPvosVOnTrz22ms899xzPPvss7z00kuebY0bN47XX3+defPmeULphRdeoFq1auTm5jJw4EDWrl1L+/btz1jPfffdx1NPPQXATTfdxLfffsuIESNYs2YNAO+88w7z58+nUaNGJ31u7NixvPvuuzRv3pxly5Zxzz338PPPP/PAAw/wwAMPcP311/Puu+961g8MDGTGjBmEhYVx+PBhevTowciRI0+6iezll1/m7bff9ux748aNTJ06lUWLFuHv788999zD5MmTufnmmz2fOdM6kZGR7Nu3j3Xr1gFw7NgxAEaNGsUdd9wBwIQJE/j3v//N/ffff9J3W7RoEZ06dTqtrWbMmMHKlStZtmxZoe0AcPToUX7++We++eYbRowYwaJFi/jwww/p2rUra9asITo6usCfydnaKioqihUrVpzx5ykicjYVLuid8Msvv3DFFVd4Ri0aNWoUCxcuZOTIkTRq1IgePXp41vXx8eHaa68F4MYbb2TUqFGFbn/atGm8//775OTksH//fjZs2HDWoJ83bx6vvPIKaWlpHDlyhLZt2zJixAgATzAtXLjwpM+kpKSwePFirr76as+8zMxMAJYsWcLMmTMB+NOf/sQjjzwCuC9HjB8/ngULFuDj48O+ffs4ePAgtWvXPmNtc+fOZdWqVXTt2hWA9PR0atWqVaR1YmJi2LFjB/fffz/Dhw9n8ODBAKxbt44JEyZw7NgxUlJSGDJkyGn73b9/PzVrnjzw0759+xg3bhxz5syhUqVKhbYDwIgRIzDG0K5dOyIiImjXrh0Abdu2ZdeuXWcM+rO1la+vLwEBARWqnwARKT4VLujPduRdUs527b2wIQsLe4Rq586dvPrqq6xYsYKqVasyevTosz5bnZGRwT333MPKlStp0KABzzzzjGf9/fv3c/vtt/PNN9+cdlOgy+WiSpUqniPvopg8eTIJCQmsWrUKf39/IiMjC33u21rLLbfcctJZjKKuk5yczG+//cacOXN45513mDZtGh999BGjR49m5syZdOjQgUmTJhEbG3vaNoOCgkhKSjptH48//jht2rTxzC+sHU78g8DHx8fz/sT02W6oK6ytMjMzPfdNiIicC12jLwV9+/Zl5syZpKWlkZqayowZM+jTp0+B67pcLs/d8Z9//jkXXXTRaeuEhoaSnJwMwPHjxwkODiY8PJyDBw8ye/bss9ZyIjxq1KhBSkqKZ1/Z2dlcc801/O1vf6NFixanfS4sLIzGjRvz3//+F3AH4W+//QZAjx49+PLLLwGYMmWK5zNJSUnUqlULf39/5s2bx+7dBY+i6O/vT3Z2NgADBw5k+vTpHDp0CIAjR46c9rkzrZOYmIjL5eLKK69k4sSJrF69GnD/A6BOnTpkZ2ef8Vp369at2bZtm2f61VdfJTAwkHvvvbfI7XAhztZWiYmJ1KxZU13dish5qXBH9E7o1KkTo0ePplu3bgCMGTOGjh07smvXrtPWDQ4OZv369XTu3Jnw8HCmTp162jpjx47lkksuoU6dOsybN4+OHTvStm1bmjRpQu/evc9aS5UqVbjjjjto164dkZGRntPfixcvZsWKFTz99NOeG/pmzZoF/HFWYfLkydx99908//zzZGdnc91119GhQwfPUwCvvfYaw4cPJzw8HIAbbriBESNG0KVLF6Kjo2nVqlWBNY0dO5b27dvTqVMnJk+ezPPPP8/gwYNxuVz4+/vzzjvvnHS/QJs2bQpcx+VyMWrUKFwuF4DniH/ixIl0796dRo0a0a5dO88/kvLr27cvDz/8MNZajDFMmDCB+vXrn3SqvWrVqsybN++M7XA+TrTt2dpq3rx5DBs27Ly2LyLieazLW16dO3e2p9qwYcNp88qT48ePO7bv++67z3700UdnXSc1NdW6XC5rrbVffPGFHTlyZGmUdpoLbadx48bZH3/8sZiqKdzhw4dtw4YNC13viiuusJs2bSpw2fn8bs+bN++cP1MRqZ2KTm1VNCXZTsBKe4Zc1Kl7OaMnn3ySZcuWMXLkyLOut2rVKqKjo2nfvj3/+te/eO2110qpwuI1fvz4c+6I6HzFx8fTs2dPz42LZ5KVlcXll19Oy5YtS6UuEfE+OnUvZzRx4kQmTpxY6Hp9+vQpluvUTouIiCj0HzXFpW7dumzZsqXQ9QICAk56tFBE5FzpiF5ERMSLKehFRES8mIJeRETEi+ka/amSt8PG12DXZ5CTAn4hEHkjtH4YQps6XZ2IiMg50RF9fvGzYVZ72P4h5CQD1v3f7R+658efvTOa0lbYEKcXqqjD4cbExLBy5cpz2vZVV13Fjh07APjoo49o164d7du3JyoqqsDBb0rCI4884umjXkTEWynoT0jeDguvgtw0sNknL7PZ7vkLr3KvV0YUNsTphTqfce+LYv369eTm5tKkSRPi4uJ44YUX+OWXX1i7di1Lly49az/9xen+++/n5ZdfLpV9iYg4RUF/wsbXwJV99nVc2bDpjXPe9K5du2jVqhVjxowhKiqKG264gZ9++onevXvTvHlzli9fDsDy5cvp1asXHTt2pFevXmzevNld2saNniFZ27dvz9atW4E/hji99tprPb3YgXu41C+//JJJkyYxatQohg4dSvPmzXn00Uc963zxxRe0a9eOqKgoHnvssdNqLmg43B9++IGePXvSqVMnrr76alJSUk77XGHbBXcPe5dddhkAhw4dIjQ01PNdQkJCaNy4MQBr1qyhR48etG/fniuuuIKjR48C7jMIjz32GN26daNFixaeAXhyc3N55JFHPGcH/vnPfwLu5/z79etH586dGTJkCPv37wegUaNGJCYmcuDAgUJ+giIi5diZetIpr6/z7hlvaqi1kyn8NS2s8G2dYufOndbX19euXbvW5ubm2k6dOtlbb73VulwuO3PmTHvZZZdZa61NSkqy2dnZ1lprf/zxRztq1ChrrbVjx461n332mbXW2szMTJuWlmattTY4ONhaa+1XX31lb775Zs/y+vXr27S0NPvxxx/bxo0b22PHjtn09HTbsGFDu2fPHrtv3z7boEEDe+jQIZudnW379+9vZ8yYcVrdjRo1sgkJCdZaaxMSEmyfPn1sSkqKtdbal19+2T777LPWWmv79etnV6xYUeTt9u3b165du9Zaa21OTo4dPHiwbdCggR09erT95ptvPOu1a9fOxsbGWmutffLJJ+0DDzzg2d9DDz1krbX2u+++swMHDrTWWvv666/bUaNGedowMTHRZmVl2Z49e9pDhw5Za62dMmWKvfXWWz37GDNmjJ0+fXoRfoplk3rGKzlqp6JTWxWNUz3j6Wa8E3JOPzotUHYR1ztF48aNTxqydODAgZ7hTE/0eZ+UlMQtt9zC1q1bMcZ4Bnrp1q0bL774InFxcYwaNYrmzZuftO1LLrmEcePGkZmZyffff0/fvn0JCgoC3APAnOh7vk2bNp7BX2JiYjzDst5www0sWLCAyy+//Iz1L126lA0bNnj60s/KyqJnz54nrbNixYoibTf/kLC+vr58//33rFixgrlz5/Lggw+yatUqHnzwQY4dO0a/fv0AuOWWW04aGvbE8L2dO3f2tF9sbCz33Xcffn7uX+tq1aqxbt061q1bx6BBgwD3UX+dOnU826lVqxbx8fFn/N4iIuWdgv4Ev5C8G/AK4R9S+DoFOHXI0vzDmZ4YvvTJJ5+kf//+zJgxg127dhETEwPANddcQ0xMDN999x1Dhgzhww8/ZMCAAZ7tBQYGEhMTw5w5c5g6dSrXX399gfv19fUlJycHe5Zhc8/EWsugQYP44osvzrpOUQQFBZ00BKsxhm7dutGtWzcGDRrErbfeyoMPPnjWbZz4Xie+04n9nzqsr7WWtm3bsmTJkgK3k5GR4flHkYiIN9I1+hMibwRTyDCgxh8ibyqxEpKSkqhXrx4AkyZN8szfuXMnTZo0Ydy4cYwcOZK1a9ee9tnrrruOjz/+mIULFzJkyJCz7qd79+7Mnz+fw4cPk5ubyxdffOE5cs4v/3C4PXr0YNGiRZ6hXNPS0k7rwrWo280/JGx8fLxnOFlwX5dv1KgR4eHhVK1a1XP9/dNPPy1wW/kNGDCAd9991xP8R44coWXLliQkJHiCPjs7m/Xr13s+s2XLFqKios66XRGR8kxBf0Lrh8GnkKD38YdWZz/SvBCPPvooTzzxBL179yY3N9cz/6uvviIqKoro6Gg2bdpUYN/ngwcPZsGCBVx88cUEBAScdT916tThpZdeon///nTo0IFOnTp5bo7L78RwuP3796dmzZpMmjSJ66+/nvbt29OjRw82bdp0XtsdPnw4sbGxgDt4H3nkEVq1akV0dDRTp07lH//4BwCffPIJf/nLX2jfvj1r1qzhqaeeOuv3uuWWW2jYsCHt27enQ4cOfP755wQEBDB9+nQee+wxOnToQHR0tOdphezsbLZt20aXLl3Oul0RkfLMnM9p3LKsS5cu9tRnujdu3Ejr1q0L/3D8bPcjdK7skx+xM/7ukO8zHepeUswVFy45OZnQ0NBS329JSU9Pp3///ixatAhfX99i2+65ttOMGTNYvXp1kQbuKauK/LudT2xsrOeykJyZ2qno1FZFU5LtZIxZZa0t8KhFR/T51b0Ehq2FZmPBPwzwcf+32Vj3fAdC3hsFBQXx7LPPsm/fPkfryMnJ4eGHH3a0BhGRkqab8U4V2hS6vu1+SYkp7D6C0pD/Ln4REW9VYY7ove0ShYh+p0WkKCpE0AcGBpKYmKg/jOI1rLUkJiYSGBjodCkiUsZViFP39evXJy4ujoSEBKdLOS8ZGRn6g14EFa2dAgMDqV+/vtNliEgZVyGC3t/f39N/enkUGxtLx44dnS6jzFM7iYicrkKcuhcREamoFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFHA16Y8xQY8xmY8w2Y8zjZ1nvKmOMNcZ0Kc36REREyjvHgt4Y4wu8A1wCtAGuN8a0KWC9UGAcsKx0KxQRESn/nDyi7wZss9busNZmAVOAywpYbyLwCpBRmsWJiIh4AyeDvh6wN990XN48D2NMR6CBtfbb0ixMRETEW/g5uG9TwDzrWWiMD/AGMLrQDRkzFhgLEBERQWxsbPFUWEakpKR43XcqCWqnolE7FY3aqejUVkXjVDs5GfRxQIN80/WB+HzToUAUEGuMAagNfGOMGWmtXZl/Q9ba94H3Abp06WJjYmJKsOzSFxsbi7d9p5KgdioatVPRqJ2KTm1VNE61k5On7lcAzY0xjY0xAcB1wDcnFlprk6y1Nay1kdbaSGApcFrIi4iIyJk5FvTW2hzgPmAOsBGYZq1db4x5zhgz0qm6REREvImTp+6x1s4CZp0y76kzrBtTGjWJiIh4E/WMJyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl5MQS8iIuLFFPQiIiJeTEEvIiLixRT0IiIiXkxBLyIi4sUU9CIiIl7Mz+kCRKRsstaSmeMix2XJzbXkuFzkuqx7Ou+/Obnu5caAjzH4GDDGeN77GIOPj/u9v68Plfx8CPT3xc/HYIxx+iuKVAgKehEv43JZjqVnk5T3Ong8g7ij6cQdTSPuaDqb9qZzLHYOyRk5Tpd6Rn4+hnpVg6gbHuT+b5Ug6lUJJDzIn7Agf6oHV6JacABVK/vj56sTkyJno6AXKSdSM3PYn5TOjoRU1u1LYuOBZDYfSGbPkbTz2JqrWGry8zH4+/rg62PIdVmy847wL1SOy7I7MY3dief23YL8fWkREUKbumG0rRtOvSpB1A4PpHZYIFUq++ssglRICnqRMsBaS0JKJtsOprB05xFW7T7Cqt1Hycg+/0CuHOBLw2qVqV+1MvWrBnlee7duoF+vboQF+hMS6EeQvy++Ps4G4InLBBnZuaRk5nA8PYfjGdkcT8/maFoW8ccy2HcsnX1H09lzJI19x9IL3E56di6/xSXxW1wSsPeM+2tUvTJdGlWjc6OqRFavTINqlakTHqizA+KVFPQipcTlshxMzuC3vUks2JrAgi0JxB0tOLDOJCzQj44Nq9Khfjit64TRsHplaoRUomrlAAL8ihZSsYc30yIi9Hy+QokxxhDo70ugvy9VKgdA1aJ/1lpLSmYOiSlZJKZmkpCcxYGkdH7fd5xf9x5lR0LqaZ85cbbgy9VxBW6zV9Pq1PbJIrPmAZrVCqFhtcr46x8BUk4p6EWKmctliTuazoKtCcxet59F2xKL/NnoBlXo1bQ6XSOr0bRmCLXCKhHo71uC1ZZ/xhhCA/0JDfQnskbwGdez1nIsLZv4pHT2H8tg79E0lu88wi9bD5OcefL9Cou3u39mX21dddL8uuGBDG5bm/b1w2laM4SmtUIIqaQ/o1K26TdU5AIkpWezZPth/rsyjrmbDhXpMzEta9K3eU26N6lGZPVgghUUpcIYQ9XgAKoGB9C2bjgAt/Zu7FlureVQcia7E9PYlZjKD8s3sCnZ/6SzLvFJGUxavOu0bTeuEcyI9nVoV78KreuEUq9KkO4HkDJDf2FEiuhQcgZz1h9k+sq9edeAz6ySnw/D29VhSFRtOjeqSvXgAP3hL+OMMUSEBRIRFki3xtWolbKdmJgYAHJyXew+ksbWgyls3H+cuZsOsm7fcc9ndx5O5a2ft520PT8fw6hO9ejRpDodGlShcfVgfBy+F0IqJgW9SAEysnNZtO0wkxbvYuHWw2ddt1fT6lzdpT59mtekRkilUqpQSpOfr4/7VH3NEIZG1ebBQS0AyHVZdh5OYX38cVbsOsKM1ftIzcoF3E8OTFsZx7SVf9wHUC04gKu71Kdzw6pEN6hCrbBAR76PVCwKehFg75E0Plu6m3//svOsj4eN7FCXqzrXp2tkNYICdO28ovP1MTSrFUqzWqFcFl2P5y9vh7WWA8czWL/vOL/vS2LGr/s8j0AeSc3ivfk7TtpG18iqXNq+Ll0jq9GydqjjT0CI91HQS4W0IyGFDxbu4IvlZ34EK7pBFW7tHcmAVrUIDfQvxeqkPDPGUCc8iDrhQVzcJoIHB7XAWvcNmr/uPcaS7Yl8sXyPZ/0Vu46yYtdRz3TLiFAu71iPbo2r0q5elSI/TSFyJgp6qRB2JKTwr9jtTF9V8ONUAHf2bcIN3RvRoJpupJLiZYyhQTX38/ojO9TlpVHtyMjOZcP+46zefZSvVu9jw373Nf/NB5P52/ebPJ+tXzWIP3VvSJ9mNWlbN0zX+eWcKejFK6Vm5jB52W5enLWpwOU+Bh4Y2II/dW9IzVBdV5fSF+jvS6eGVenUsCpj+jTB5bJsPpjMsh2JzFp3gOU7jwAQdzSdV77fzCtsBqBP8xpcElWHPs1r0KBaZSe/gpQTCnrxCtZaNiTm8uo/F550N3R+jwxuwY09Grk7ZBEpY3x8DK3rhNG6ThijezfGWsu2Qyks2ZHIV6v3sWbvMQAWbj180g2iN3RvyKA2EfRoUl19LkiBFPRSbmXm5PLZ0j1M/HZDvrkZnnfD29fhwYtb0KxWSOkXJ3KBjDE0jwileUQoN/eMJDvXxa97jrFgSwKTFu8iJa+Tn8nL9jB5mfuaf6+m1RnZoS4DWtXSHf3ioaCXcuVwSiavztnMlBWn30RXLdDw4lUdubh1hPosFylH7xgAACAASURBVK/j7+tDt8bV6Na4Go8MacnR1Cx+2XaY79cd4Lvf9wPuHv1O9OpXIySAG3s04pKoOrSICNF9JxWYgl7KvPhj6fx1xu/M25xw2rKBrWox4dI2NK4RTGxsLDFRdRyoUKT0VQ0OYESHuozoUJd/5LpYvecYP244wAcLdwJwOCWLN3/ayps/bcUYuDemGZe0q02bOmEK/QpGQS9l0sHjGUyYuY4fNxw8bdmd/Zpwb/9mhOmRNxHA3aHPiaP9vw5vw/aEFH7acJDPl+9hd2Ia1sLb87bx9jx373139WvKsHa1aVcvXKFfASjopcw4kprFU1+v49u1+09b9tjQVtx2USSV/HSzkUhhmtYMoWm/EO7s15T9SenM/v0Any7dzc7D7pH83p2/nXfnb8fXx3D/gGZcHl3vrAMCSfmmoBdH5eS6+GDhzpOeGz7h4UEtGNuvicJd5ALUCQ/itosac9tFjYk/ls6s3/fz+bI97DicSq7Lek7vN6sVws09G3Fp+7pUC9aTKd5EQS+OWLztMDd/tPy07mbv6teUP1/cXI8JiZSAulWCGNOnCWP6NCHuaBpfr4nnnXnbSMvKZduhFJ76ej1Pfb2emJY1ubZLAy5uE4G/bmwt9xT0UmqOpWXxwJQ1zN9y8k11vZtV5/VroonQ40AipaZ+1crc278Z98Q05be4JL5cFcenS3cDELs5gdi8m1/v7NeE67o2pLFO7ZdbCnopcT+sP8DYT1edNC/Q34dPbu1G9ybVHapKRMD9vH50gypEN6jChEtbM29TAv9duZe5mw4B8N78Hbw3fwedG1Xlph6NGBpVW2fcyhkFvZSIlMwcHpn2G9+vP3DS/Pv6N+OBi5vrdKBIGVTJz5ehUbUZGlWbg8czmLZiL/+Yu5Ucl2XV7qOs2n0UpsIdfRpzc89IdcFbTijopVhtiD/OsLcWnjQvtJIfU+7sQdu64Q5VJSLnKiIskPsHNuee/s2Yv+UQk5fu8Rzlf7BwJx8s3MmgNhHc2jsSa888tLM4T0EvxeKr1XE8NO23k+bd2juSJy5prWE2RcoxXx/DgFYRDGgVQfyxdKYs38NbP7ufx/9xw0F+3HCQiMqGPwfv4fLoegQF6LR+WaOgl/OWk+vimf+t57Ole06a/+HNXbi4TYRDVYlISalbJYiHBrfknv7N+Oa3eP4vdjs7D6dyMM3yxFe/88RXvzNuYHNG94rUI3pliIJezll6Vi5j/rOCRdsSPfOqVvZnxj291emGSAUQ6O/LNV0acHXn+izZnsjfvl7Jbwm5ALw1dytvzd3K6F6RjOnTmPpVdR3faQp6KbKk9GyufW8Jmw4ke+YNbFWLt//USafrRCogYwy9mtXgwc6BNIrqygcLd/B53kh6kxbvYtLiXYzqWI87+zWlZe1Qh6utuBT0UqhDyRmM+OcvHDye6Zl3Q/eGPDuyrUaJExEAGtcI5sUr2vHAwOZ8uHCHZ3Cdr37dx1e/7uOSqNr8+eIWCnwHKOjljI6mZjH4zQUkJP8R8PcPaMZDg1poIAwRKVBEWCB/Hd6Ge2KaMWnxLv4xdysAs9cdYPa6A1zavg5/vrg5zWop8EuLgl5Ok5qZw5X/t/ikU/QThrfm9osaK+BFpEiqBgfw4KAW3NG3CZ8s3sXf52wG4Nu1+/l27X4ui67LuIHNaVozxOFKvZ+CXjyyclzcOmn5STfZPTSoBfcPaKaAF5HzElLJj3v7N+PGHo3498Idnkfzvl4Tz9dr4rmuawMeGtSCWuoCu8Q4eoHVGDPUGLPZGLPNGPN4AcsfMsZsMMasNcbMNcY0cqJOb2et5Zlv1tNiwmxPyI/uFcnOl4YxbmBzhbyIXLDwIH8eGtyS1U8OYmzfJp75U1bspduLc3n9xy2kZuY4WKH3cizojTG+wDvAJUAb4HpjTJtTVvsV6GKtbQ9MB14p3Sq937dr42n8xCwmLd4FwPB2ddj+4jCeGdlWAS8ixa5acADjh7Vm2fiBXN+tgWf+W3O30uvln/ls6W5ycl0OVuh9nDyi7wZss9busNZmAVOAy/KvYK2dZ61Ny5tcCtQv5Rq91paDyUQ+/h33ff4rAHXCA1n37BDeuaETvj4KeBEpWRFhgbw0qj3f/7kPfZrXANyP8E6YuY6h/1jIom2HHa7Qezh5jb4esDffdBzQ/Szr3w7MLtGKKoCM7FwGvjaffcfSPfN+eqgfzWrphhgRKX2taofxn9u6EbslgRe+28i2QylsO5TCDR8uY3i7Ovx1eGvqVglyusxyzTg1GIEx5mpgiLV2TN70TUA3a+39Bax7I3Af0M9am1nA8rHAWICIiIjOU6ZMKdHaS1tKSgohIRcexLN3ZjN1c5Zn+r7oSnSp7T33YxZXO3k7tVPRqJ2KrrjaKtdlmbc3h6mbs8jOd/b+yub+DG3sj385P9tYkr9T/fv3X2Wt7VLQMieDvifwjLV2SN70EwDW2pdOWe9i4J+4Q/5QYdvt0qWLXblyZQlU7JzY2FhiYmLO+/M7ElIY8Np8z/Tl0XV549por7sGf6HtVFGonYpG7VR0xd1Wh1MyeXn2JqavivPMi6xemReuaEfvZjWKbT+lrSR/p4wxZwx6J6/RrwCaG2MaG2MCgOuAb/KvYIzpCLwHjCxKyMvJXC7Lnz5YelLILx8/kDev6+h1IS8i3qNGSCVevboD/72rJy0j3B3r7EpM44YPl/Ho9N9ISst2uMLyxbGgt9bm4D4dPwfYCEyz1q43xjxnjBmZt9rfgRDgv8aYNcaYb86wOTnFqt1HaDJ+Fou3ux+Xe/PaaHa9PFzPqopIudE1shrfjruICcNbc+Ks/bSVcQx8fT6zf9/vbHHliKMXaK21s4BZp8x7Kt/7i0u9qHLO5bKM+r/FrNl7DIAmNYP54c991Se9iJRL/r4+jOnThMFtavPYl2tZsiORwymZ3D15NUPaRjDx8ihqheoA5mz019+LrNzlPoo/EfKf3d6dnx+OUciLSLnXsHplJo/pzgtXRBHo7/6bNmf9QYa8sYA56w84XF3ZpgTwAtZabpu0gqveXQJAy4hQtr84jIual9+bVkRETuXjY7iheyPmPhxD3xY1ATials2dn67iL//9jRT1rFcgBX05d/B4Bo2fmMXPm9z3Kn52e3fmPNhXnd6IiNeqVyWIT27tysTL2nr+1v13VRyX/GMBK3cdcbi6skdBX459sngX3V+c65neNHGojuJFpEIwxnBTz0hmP9CHVnlj3O89ks5V7y7hrblbcbmceXS8LFLQl0O5LkuX53/k6W/WA/CXIS3Z9fJwAv19Ha5MRKR0tYgI5ev7enP7RY09817/cQu3fLycxJTT+lerkBT05cz+pHSajp/F4RR3D3cLH+3Pvf2bOVyViIhzKvn58uSlbZh0a1eqVvYHYOHWwwx7ayHLd+pUvoK+HPl2bTw9X/oZcF+j2vHiMBpUq+xwVSIiZUNMy1p8O64PHRpUAeDg8UyueW8J783fjlO9wJYFCvpy4s5PV3pGmnt4UAsWPT4AH91wJyJyknpVgph2Zw9u6tHIM++l2Zt4cOoaMrJzHazMOd4zoomXynVZRn+fCqQC8M19vWlfv4qzRYmIlGGV/HyZeHkUnRpV4fEvfyczx8XMNfHsPJzKezd1oXZ4xepgR0f0ZdixtCyajv+j48D1zw5RyIuIFNEVHesz7c6e1AqtBMBvcUmMfPsXT6diFYWCvozaEH+c6Od+BCAsAHa+NIzgSjoBIyJyLjo0qMLX9/Wmbd0wAA4lZ3Lte0sqVG96Cvoy6Pt1+xn21kIARnWqx1sDgjXanIjIeaoTHsR/7+rJkLYRAGTmuLjz01V8unS3w5WVDgV9GfPJ4l3c9dlqAF68oh2vXxPtcEUiIuVf5QA//u+GztzW+4/n7Z+cuY6/z9nk9Xfk61xwGfK37zfxf7HbAfj09m70aV7T4YpERLyHj4/hyUtbUyM0gFe+3wzAO/O2c/B4Jn+7sr3Xdh2uoC8jHpjyK1+viQfg2/svIqpeuMMViYh4H2MM98Q0o3pwAE989TsuC9NXxZGZ4+L1azrg74WjfSroy4DbJ61gbt6gNAv+0p+G1dUJjohISbq2a0OqVA7g/s9/JSvXxf9+iycrJ5d/Xt+JAD/vCnvv+jbl0B3/WekJ+RV/vVghLyJSSoa0rc17N3UmwPeP8e3v+myV13Wso6B30N2freLHDQcBWDnhYmrmPespIiKlo3+rWrx/c2fPUfzPmw5xz+TVZOW4HK6s+CjoHfLnKb8ye537Oc4Vf72YGiEKeRERJ8S0rMUHN3c5KewfmraGXC8Z6lZB74C/z9nEzLwb75aPH6gjeRERh/VrUZP3buyMX96d99+u3c9fZ/zuFY/eKehL2dQVe3hnnvsRuoWP9qdWWMXqc1lEpKzq36oWr17dwTM9ZcVeXpq9ycGKioeCvhQt3JrAY1/+DsDMe3triFkRkTLm8o71eHpEG8/0+wt28MniXc4VVAwU9KVk5+FUbvr3cgDeu6kz0Q00OI2ISFl0a+/G3Ne/mWf62f+t5+dNBx2s6MIo6EtBRnYu/V+NBWD8sFYMaVvb2YJEROSsHh7cgkvb1wHAZeG+z39lfXySw1WdHwV9CbPW0urJ7wHo3aw6Y/s2dbgiEREpjDGGV6/uQIf67l5K07JyuX3SShKSMx2u7Nwp6EvYyLcXed5PHtPDwUpERORcBPr78sHNXaidd9P0geMZ3P/FanJyy9cz9gr6EvTp0t38vs99qmfL85c4XI2IiJyrWmGB7mfs83rPW7rjCH+fs9nhqs6Ngr6ExB1N48mZ6wCIfSTG6/pOFhGpKNrVD+fJfHfiv7dgB7N/3+9gRedG6VMCrLVc9Ld5AEwY3prIGsEOVyQiIhfixu4NGdmhrmf60S/Xsu9YuoMVFZ2CvgQMf+sXAKoFBzCmTxOHqxERkQtljOHFUe1oUtN94JackcNDU8tHN7kK+mK2YEsCG/YfB9x92IuIiHcIqeTHm9dG45vXTe6ynUf4YOEOh6sqnIK+GOXkurj5I3enOF/d08vzyyAiIt6hff0qJ3Wm89oPm9kQf9zBigqnoC9GA1+fD0C3xtXo1LCqw9WIiEhJuG9AM9rVC8cYyHFZnvhqbZk+ha+gLyZLtieyOzENgKlj9by8iIi38vf14bVrOuBrDP6+PvwWl8R/luxyuqwzUtAXk+s/WArA9Lt6YoxO2YuIeLMWEaHc0bcJWTkuggN8eXXOZuLL6F34Cvpi8Nz/NgAQWsmPLpHVHK5GRERKw/0DmlGvShAAWbkunv9ug8MVFUxBf4FSM3P4aNFOAJaMH+hwNSIiUloqB/jxzMi2pGblElzJj1m/H2DlriNOl3UaBf0FGvqPBQCM7hVJSCU/h6sREZHSNKhNBBc1q0FGdi7BAb5M/G4jrjJ2Y56C/gIcOp7B3iPuazJP5+seUUREKo7HhrYiI9tFaKA/v+09xv/Wxjtd0kkU9BfgRDe3z4xooxvwREQqqHb1wxnZoS5H0rKoHhzAmz9tLVMj3Cnoz9OexDSy8n6Qo3s3drgaERFx0iODW5LrsgT6+7LzcCrfri07g94o6M/ToDfcneO8cW0HhysRERGnNaxemcs61CUxNZNaoZX4589by0wnOgr685CUnk1mjvto/oqO9R2uRkREyoK7YprmXav3Y3tCKt+vO+B0SYCC/rzc9O9lAIwb2NzhSkREpKxoERHK4DYRHErOpGZoJT7Oe/TaaQr6c+RyWdbGJQHw4MUKehER+cOYPk1IzsihamV/Vu4+yrp9SU6XpKA/V/+K3QZAjybVdKe9iIicpGtkVVpEhJCV4yLI37dM9IGvoD9Hr/6wBYB3/tTJ4UpERKSsMcZwQ/dG7EpMo1mtEL5eE09yRrajNSnoz0Hc0TTP++ohlRysREREyqorOtUj0N+HXJclM8fl+E15CvpzMO6LXwGYeFlbhysREZGyKizQn4GtIzhwPIP6VYOYuWafo/Uo6M/B6j3HALixRyOHKxERkbJsRPu6HEnNIiIskMXbEzmQlOFYLQr6Itp75I/T9roJT0REziamZU3PQGfWwo8bDzpWi4K+iMbP+B3Q4DUiIlK4QH9fBrSqxe7EVBpWq8zPCvqyb+HWwwDc3DPS2UJERKRc6NeiJodTsqgdHsii7Ylk5jjTJa6Cvgjyjy3s66PT9iIiUrg+LWoA4OdjyMpxsfFIriN1KOiLYNY69yhEXRpVdbgSEREpL2qFBtK6Thjp2bn4+xq2HHVm6FoFfRG88aO7k5yHB7d0uBIRESlPujeuxuYDybSpE8bWozqiL7O2J6QC7m5vRUREiiq6QRXSsnIJC/JnZ5KLjOzSD3sF/TnQY3UiInIuohtUAcDHGHIsbDqQXOo1KOgLkZiS6XQJIiJSTjWqXpnwIH9yXO7r81sOKujLnNl5fRRf0bGew5WIiEh5Y4yhac1gsnMs/j6wRUf0Zc+CLQmAu5cjERGRc9WkZgi7ElOpG+LDlkMppb5/BX0hTnSUc1GzGg5XIiIi5VGTmsEcSs4kLMAQl6879dKioC9Eet4dkhqWVkREzkf9qpUBsMD+pAysLd0e8hT0IiIiJahm3oGir3EfPB5PzynV/Svoy6rk7bD8HpgWRr/4ATAtzD2dvN3pykRE5BzUDHUH/Yke1A8cL90hax0NemPMUGPMZmPMNmPM4wUsr2SMmZq3fJkxJrL0q3RA/GyY1R62fwg5yRgs5CS7p2e1dy8XEZFyoUZIAAA5eT3gpmRWkCN6Y4wv8A5wCdAGuN4Yc+oYsLcDR621zYA3gL+VbpUOSN4OC6+C3DSw2Scvs9nu+Quv0pG9iEg5ERTgC/xxRJ+WVUGCHugGbLPW7rDWZgFTgMtOWecy4JO899OBgcbbu6fb+Bq4ss++jisbNr1ROvWIiMgFCfB1R+2JoE/NLN1ucAsNemPMfcaYkhi2rR6wN990XN68Atex1uYASUD1Eqil7Nj12elH8qey2bDr09KpR0RELogxhkp+f8RtZk7pBr1fEdapDawwxqwGPgLm2OJ5NqCgI/NTt1uUdTDGjAXGAkRERBAbG3vBxZ2qJLZZkH45KQV+6VPZ7GTml1JN5UVKSkqp/ZzKM7VT0aidik5tVbhKPpas7FzAsHnjRmKPbS21fRca9NbaCcaYJ4HBwK3A28aYacC/rbUXcqE4DmiQb7o+EH+GdeKMMX5AOHCkgBrfB94H6NKli42JibmAsk7x/XcAFOs2z2ZaiPvGu0IY/9DSq6mciI2NVZsUgdqpaNRORae2KtzaGPj0fz/z5KJ02rdrS0xUnVLbd5Gu0ecdwR/Ie+UAVYHpxphXLmDfK4DmxpjGxpgA4Drgm1PW+Qa4Je/9VcDPxXQ2oeyKvBGM/9nXMf4QeVPp1CMiIsUi1+WOL59SvtWsKNfoxxljVgGvAIuAdtbau4HOwJXnu+O8a+73AXOAjcA0a+16Y8xzxpiReav9G6hujNkGPASc9gie12n9MPgUEvQ+/tDqwdKpR0REisWJfnJCKhXlqnnxKcreagCjrLW788+01rqMMZdeyM6ttbOAWafMeyrf+wzg6gvZR7kT2hT6THc/QufKPvnGPOPvDvk+093riYhIuZGS7T6irxocUKr7LfSI3lr71Kkhn2/ZxuIvSah7CQxbC83Ggn8YFgP+Ye7pYWvdy0VEpFxJycoL+sqlG/Sle/5Aii60KXR9G7q+zXzd6CIiUu4dzwv6KpULuTxbzNTXfRG5XN59D6CIiJSshHRLrdBKBPr7lup+FfSFaFozGIAN+487XImIiJRnh9JcNKpeudT3q6AvxOC2tQGYvW6/w5WIiEh5dijNUjO0Ep8u2UX8sfRS26+CvhCXtnd3avDtWgW9iIicn6S0bI5lWtKzcnny6/VsLMWzxAr6QrSpEwbA7sQ0hysREZHyal180knT7eqFl9q+FfSF8PbB8kREpOT9vs8d9BaICKtErbDAUtu3gv4cHDye4XQJIiJSDv0el0T1QMOeI2mlejQPCvoiuaNPYwA+WrTT4UpERKS8sdaybGcitSobdiSk0jWyWqnuX0FfBLdf1ASA9+bvcLgSEREpbzYfTOZwShY+eVeCezerUar7V9AXQe3w0ruWIiIi3mXxtkTP+/Agf1rn3eRdWhT0RRQc4O7JaN2+pELWFBER+UPslgQiq1cmPsXSq2l1fH3K2DC14vb8FVEAPPX1OocrERGR8iIpLZvF2w5TKzSQo5mWwW0jSr0GBX0RXR5dD4DVe445XImIiJQXP208SE7eWCm+Bga2VtCXWfmfp99yMNnBSkREpLyYve4AdcIDiU9Kp20NX8ICS3fkOlDQn5O3ru8IwF2frXK4EhERKesOp2QSu/kQtUIrEXc0ne61S3fUuhMU9OdgZIe6AOxISHW4EhERKetmrN5HjstigdBKfnSp7edIHQr6c9SgWhAA01bsdbgSEREpq6y1TFmxhxYRIWw5mMzI6LpU8nWmS3UF/Tn6fEwPAB79cq3DlYiISFm1fOcRtiekUsnPl4xsF9d1behYLQr6c9SgWmXP++0JKQ5WIiIiZdUHC3cQHuTPweMZRDeoQlS90u0kJz8F/Xl450+dALjkzYUOVyIiImXNtkPJ/LTxEBFhlTiUnMmdfZs4OhKqgv48DG9fB4CsXBeHNKKdiIjk88GCnVTy8yE9O5fI6pUZ3La2o/Uo6M/TsyPbAjD4zQUOVyIiImXF3iNpfPVrHHXCA9l7JJ0xfZqUepe3p1LQn6dbekUCcCwtm4TkTGeLERGRMuEfc7diMKRm5VKvShBXd6nvdEkK+gsxYXhrAHq+NNfhSkRExGlbDybz1eo4aoVVIiE5kwcHtaCSnzOd5OSnoL8AY/q4x6nPcVmNaiciUsG99sMWAvx8SM/KpVmtEK7oWM/pkgAF/QX75LZuAFz6z18crkRERJyyaNthvl9/gOAAPxJTs/jLkJaOX5s/QUF/gfq1qOl5P3XFHgcrERERJ2Tnunj6m/VUDvAlOTOHPs1rMLhN6Y9SdyYK+mKw8NH+ADz25e9k5bgcrkZERErTpEW72HYohUB/X6y1PDOyraPPzZ9KQV8MGlSrzMV5YwwPemO+w9WIiEhp2XcsnTd/2kJooB9HUrO4/aImNK0Z4nRZJ1HQF5MPbu4MwO7ENFbtPuJwNSIiUtKstTz+5VrSsnNJy8qlQbUg7h/QzOmyTqOgLybGGL64wz3gzZX/t4Rcl3W4IhERKUlTVuxl4dbD+Pv4kOuyvHJlB4IrOTMU7dko6ItRz6bVaVU7FICBr8U6W4yIiJSYfcfSeeG7jfj7GrJyXdzSsxE9m1Z3uqwCKeiL2axxfQDYlZjGd2v3O1yNiIgUt5xcF3+e8ispmTlk51oaVqvMY5e0crqsM1LQFzMfH8P3f3aH/b2fryYpPdvhikREpDi98dMWVuw6CoC/r+Gf13ekckDZO2V/goK+BLSqHcZNPRoB0OHZH7BW1+tFRLzBgi0J/Ct2u2f6saGt6NCgioMVFU5BX0ImXh7leX/Vu0scrERERIrDgaQMHpq2hhPHbgNb1eL2ixo7W1QRKOhL0JbnLwFg1e6jfL5MveaJiJRXGdm5jP10JYdTsgCoVyWIV6/uUKY6xjkTBX0JCvDz4aeH+gEwfsbvbIg/7nBFIiJyrqy1PDp9LWvj3IOXVQ7w5cNbulA1OMDhyopGQV/CmtUK4eVR7QAY9tZCjqRmOVyRiIici/+bv51vfov3TL9+TTSt64Q5WNG5UdCXguu6NeSaLvUB6DTxR3Jy1R++iEh58N3a/fx9zmbP9EODWjA0qraDFZ07BX0peeWqDtSrEgRAs7/O1p34IiJl3OJth3lw6h83313RsV6Z7OK2MAr6UvTLY/097ztO/NHBSkRE5GzW7Uti7KeryMo7A9uvRU1euap9ubj57lQK+lJkjGHrC+478Y+lZTP0zQUOVyQiIqfanZjK6I9XkJKZA0CH+uH864ZO+PuWz8gsn1WXY/6+Pmx8bigAmw4kM/rj5Q5XJCIiJ+w9ksb17y/lcEomAE1qBPPR6K5lcrCaolLQOyAowJffnhoMQOzmBO79fLXDFYmIyN4jaVz3/lLikzIAaFS9Mp/f0YPqIZUcruzCKOgdEl7Zn+XjBwLuuzpvn7TC4YpERCquuKNpXP/BUvYdSwegYbXKfHFHD2qHBzpc2YVT0DuoVlggy//qDvu5mw5x3fvqKldEpLTtSXSHfNxRd8g3qBbEF2N7UDfvSanyTkHvsFqhgax+chAAS3ccYfhbCx2uSESk4ti4/zhXvruYvUfcId+ouvtIvp6XhDwo6MuEasEB/Pa0+5r9+vjjtJyg5+xFRErayl1HuPa9JSQku2+8a10njOl39aJ+1coOV1a8FPRlRHiQPxueGwJAZo6Lxk/MIitHPeiJiJSEeZsPceO/l3E8w/0IXdfIqkwZ24OaoeX7xruCKOjLkMoBfmx/cZhnusWE2RxV3/giIsXqs6W7uX3SCjKy3QdT/VvW5D+3dSc8yN/hykqGgr6M8fUx7Hp5OE1qBgPuHvT2HM91uCoRkfIv12V57n8bmDBzHa68q6PXd2vI+zd3ISjA19niSpCCvoz6+eEYRnWsB8BTizOYukLj2YuInK/UzBzG/mclHy3a6Zk3YXhrXrwiqtz2eFdU3v3tyrnXr43mlavaA/DYl79z92erHK5IRKT82ZOYxlXvLmHupkMABPn78v5NnRnTp0m57Lv+XCnoy7hrujRgYm/3Yx6z1x0g8vHvyNYwtyIiRTJv0yEu/edCNu4/DkDd8ED+e1dPBrctX0PNXggFfTnQINSHdc8O8Uw3/+tsthxMdrAiEZGyzeWyvPHjFm6dtMJzZ33vZtX53/0XEVUv3OHqSpeCvpwIqeTHzpeG0aFBFQAGv7GAd+Ztc7gqEZGy51haFrd/soJ/zN3qmXd3TFM+ubVbuC/FeAAAGWJJREFUue+3/nwo6MsRYwxf39ubV6/uAMDf52ymw7M/kOtS5zoiIgBLticy9M2FzNucALgPkt69sTOPDW2Fn5ffdHcmFfNbl3NXda7PkicGAJCUnk3T8bP4PS7J4apERJyTnevi73M2cf0HSzlw3D36XPv64Xx7/0UMjao41+MLoqAvp+qEB7HzpWG0qRMGwIi3f+GBKb86XJWISOnbk5jG1e8u4Z152z3z7urXlOl39SKyRrCDlZUNCvpyzBjDrAf68P5NnQH4ek08kY9/x/6kdIcrExEpedZavli+hyFvLmDN3mMA/H97dx4fVXnvcfzzZA9kI4QMgUBYEpawCURBBIUiSEWrVm210qq19ra96r2t7atWvb1ebSvWW221tvZWvWptpZsoFS+WLSgKKLIIsiXshJANEsgekuf+MWMkEGCUJGfy5Pt+veaVmXMOmR+/JPOdc84zz0mNj+aP35jAPZ8fRlSEIg4U9E6YOaI3Wx+c1fz4woeX8Yt/bvewIhGR9nXgSDVfe+49fvTKJmoa/LOHXjbCx6J/v5iLMlM8ri60KOgdERsVzp65s7nv8uEAPLksnwH3LGT/4WqPKxMRaTvWWv60Zh+XPraCt/NKAUjqFsmTN47l6TnjSe4e5XGFoUdB75jbLx7UfH17gCk/X87df9moy96KSKe3p7SKrz77HvfO39R8QZpZI3qz+LuXcOWYPl1ilrvPwpOgN8YkG2MWG2PyAl97tLLNecaYVcaYj4wxHxpjvuxFrZ1Rcvco9sydzU+vGQnA39cdYOCPNDJfRDqnuuON/GpJHlP/O5eV+f69+B6Bvfjfzhnn5KVl25JXe/T3AEuttVnA0sDjk1UDX7PWjgBmAb80xiR1YI2d3k0TMtj64CzioiMA/8j8Cx9eSk29roYnIp3DyrxSZv3ybR5fsqN52Q3n92PZ3VO1Fx8kr4L+KuCFwP0XgKtP3sBau8Namxe4fxAoBnp1WIWOiI0KZ/N/Xcbzt54PQGFFLcN/vIjHNFhPREJY0dFa7nx5PXOeXcPu0ioAhvWO5+/fvpC5146mh87FB814ce7WGFNurU064fERa+0ph+9PWH8B/jcEI6y1p1zRxRjzTeCbAD6fb/y8efPaoWrvVFZWEhcXd87fx1rL/35Uz1sHjjcvu39CDJk93LgOc1v1yXXqU3DUp+C1Za/qGi2LdjcwP7+heVm4geuHRDEjI4LwsM67B9+ev1PTpk37wFqb09q6dgt6Y8wSoLXpiO4DXgg26I0xaUAucLO1dvXZnjcnJ8euXbv2sxUdonJzc5k6dWqbfb/y6nrGPrSYE3/0a+6dji8hps2ewwtt3SdXqU/BUZ+C1xa9amqyvLaxgIff2Ebxsbrm5V8c25cfzBpKWmLsOVbpvfb8nTLGnDboI9rlGQFr7aVnKKjIGJNmrS0MBHnxabZLABYC9wcT8hKcpG5R7H54Nqt2lnHj7/1tnfCzpQxPS2D+dyYRE+nGHr6IdA4f7D3Mg69vZWNg0huAnIwe/McV2c0X8pLPzqtz9AuAmwP3bwZeO3kDY0wUMB940Vr71w6srcu4cHBP9sydzX9emQ3A1sKjDPuPRfzwbx/SpAvliEg723boKLe/uJZrf7uqOeT7JsXy1FfG8ddvXaiQbyNeBf1cYIYxJg+YEXiMMSbHGPNMYJsvARcDtxhjNgRu53lTrttuvWggu352OdePTwfgz2v3M+jeN3j0zW36/L2ItLk9pVX827z1zPrl2yzeUgRAVEQYP/r8MJbefQmzR6dpNH0bardD92dirS0DpreyfC3wjcD9l4CXOri0LisszPDo9WN46OqRXPf0u2wuOMpTy3fy1PKdfPfSIdw1PVN/eCJyTg6W1/Dksjxefm9/i+V3Tc/iG1MGkhAT6VFlbvMk6CV0xUSG8/qdU6ioaeCap95hV2kVjy/ZweNLdvDDWcP41iWDFPgi8qnsP1zN0yt28sc1+1osv33KQL51yWB6xmnCm/akoJdWJcZGsuz7UzlSVc8VT66koLyGRxZt45FF2/iXSwbxw8uGEdaJP+YiIu0vv7iS3+Tm88q6ghbLvzoxg3+dlknvxM79SZ/OQkEvZ9SjexTv3PM5yirrmP3ESg4dreV3K3bxuxW7mD0qjce+PIboCI3SF5FPbC6o4De5+byx6VCL5bdPGcjtUwaR2sk/ytvZKOglKD3joll973SO1jZw83PvsX5fOQs3FbJwUyGj+iby0m0TSOym82siXVVTk2XFjhKeWbmLd/LLWqy763OZ3HrRQM1m5xEFvXwqCTGRzP/ORdQfb+L7f93Igo0H2VRQwZgH/wnA63dOZmTfRI+rFJGOUlPfyLJ9Ddz38+UUlNc0L4+NDOfO6ZnMmZihQXYeU9DLZxIVEcYTN47ll18+j8cW7+DXy/MBuOLJlQA8cGU2N08aoIF7Io4qOlrLi6v28NTynS2WD/XF862pg5g9qg9REboSeihQ0Ms5CQszfP+yoXz/sqEs31bMrc+/D8AD/9jCA//YwrShvfjVjWP1jl7EAU1Nlnd2lvLS6r28+VFRi3Uzsn3cNnkgEwYm6w1+iFHQS5uZNiyVPXNnU1hRw03PrGFXSRXLt5cw+gH/Yf3ffXU8M7N9ehEQ6WTKKuv42wcH+E3uTipqGlqsu/WiAQwLK+LLs1udZl1CgIJe2lxaYizL7p7K8cYmHlm0jd+/vRuAf/nDBwBMyUrhiRvGamCOSAhrarKs3l3Gy+/t5x8bD7ZYN7JvArdMGsgVo9OIiQwnN7fEoyolGAp6aTcR4WHcNzub+2Zns7mggjnPrqG8uoG380oZ+9BiAH5y9UhuvKB/p770pIhL9pRW8fd1B3hqeT4nX/JizsT+fOWCDLL7JHhTnHwmCnrpECP7JrLhxzM53tjEE0vzeGKZf/De/a9u5v5XNxMbGc6zt+QwaXCKx5WKdD0VNQ0s/LCQP6zey9bCoy3WjUlP5KaJGVw+Ko24aEVGZ6SfmnSoiPAwvjdzKN+bOZS9ZVXc9fJ6Nh6ooKahka/8fg0Ao9MTefLGsWT07O5xtSLuqqlvJHd7Ma9uKDhlYF1URBjfmTqYa8el0y+5m0cVSltR0ItnMnp257U7JgPwbn4pt72wlpqGRj48UMElj+YCMDkzhZ9dM4r+PfViI3KuahsaWbGjhAUbDrJwU+Ep67+Uk8514/tx/oAeGjTrEAW9hIRJmSlsfWgWTU2Wl9/fx33zNwOwMr+Uix9dDsAFA5J5+NpRDO4V52WpIp1KbUMjK/NKWbDxIAtOGlQHcMXoNK46ry9TslKIidR01i5S0EtICQsz3DQhg5smZFB/vImXVu/lwde3APDensNM/8UKAIanJfCTq0cwrr/2PEROVlZZx7JtxSzYeJC380pPWX/ZCB9Xn9eXacNSFe5dgIJeQlZURBhfnzyQr08eyPHGJua9v5/7X/Xv6W8tPMq1v13VvO2DV40g7eQhwiJdyK6SShZvKWLe+/vZXVp1yvrZo9L4/KjefG5YKt2i9NLfleinLZ1CRHgYcyZmMGdiBo1Nljc2FXL/q5ubJ+/48Wsf+Tf850KuGduX7146ROf1xWlVdcdZtbOM5duLT7nO+8dumTSAmdk+zh+YTGS4pqPtqhT00umEhxmuHNOHK8f0ASCv6BgPLdzKWzv8k3bMX1/A/PWfXP/6h7OGMWdif+I1Da90Yk1Nli2FR1mxo4S/rN3P3rLqU7ZJjY9mzsQMZmT7GNY7Xqe1BFDQiwOyfPG8+PULyM3NJefCyTy3cjePLd7RvP6RRdt4ZNE2ADJ6duMHlw3l0uE+nZuUkGatJb+4ktW7D/Pm5kOszD/1XDvA1ef14dJsH5MzU0jqptkm5VQKenFKXHQEd03P4q7pWQDsLKnkiaV5vLbBP9p4b1k1d/xpffP2fZNi+ddpmXzhvD6aDEQ81dRk2VF8jDW7DvPGpkLW7D7c6nbD0xK4dlxfLh7Si6zUOO21y1nplU2cNrhXHL+6YSy/umGsf+7uXWU8vmQH7+85AkBBeQ33zt/EvfM3ARATGca3L8nk+px0+iTFelm6OK6y7jgfHihn/b5y5q8vIL+4stXtMnp24/rx6Uwc1JPR6Um69Kt8agp66TLCwgyTMlOYlOmfZrexybJ+3xGeXrGLJVv9M4PVNjTx+JIdPL7kk0P/4/oncfOkAUwf7tNev3wmTU2WXaWVrNtXzls7Snj9w1Mnq/nYoF7duW58OhMG9mR0eqIG0ck506uWdFnhYYacAck8MyAZ+OSc6PPv7mkxinndvnLW7dvQ4t/OyPZxXWAvKzFWg/zkEw2NTeQVVfLRwQrW7D7M3z44cMbtZ2b7mDYslXH9e5CZGqcLPEmbU9CLBBhjyPLF89NrRvHTa0YB/vnAV+aX8uKqPS0mHlm8pYjFW1rOD56VGscXx6Xz+ZG96Z/cjTC9YDuvrLKOvOJKdhQdY/m2YpZvP/PlWvslx/LFsemMz+jBmH5JepMoHUJBL3IGsVHhzMj2MSPb17zscFU9y7YV88q6A7y7s6x5eV5xZYsR/h8b2z+Jy0b0ZurQXgxKidM51k7GWktJZR35Rf5AX5lf1nyq50yG+uKZPTqNUemJjOiTQGp8TAdUK3IqBb3Ip5TcPYrrxqdz3fj05mW1DY1sLqjg1Q0FvLKugOr6xuZ16/f5B1zN/b+WbwBS4qK5OCuFizJTGNMviX7JsURH6CN/XrDWUnysjr1l1ewpq2LTgQoWbaymZNHCoP799GGpTMlKYWjvBLL7JGhPXUKKgl6kDcREhpMzIJmcAcn85OpRzcvLq+vZVFDB4i1FvPnRIYqO1jWvK62s45X1BbxywuQ+JxrTL4kLBvTgvH49yPLF0TsxhvjoCH2c6jOoP95E0dFaDpbXUFhRy46iY3yw9wjr9h2hoTH4qZNnZvu4KDOFLF8cQ3zx9OwepZ+HhDwFvUg7SuoWxZSsXkzJ6sWDV41sXl7b0MjesmrW7TvCyvxS3skvpby6ocW/3bi/nI37y4Hdp/3+g3t1Z3haAkN88WSlxpHeoxvJcVH07B7l/IRATU2WipoGyqrqKDlWT2llHQeO1LD90FG2HTrGtkPHPtP3nZKVwqTBKdQU7eaKqRPon9zN+V6K2xT0Ih6IiQxnaO94hvaO58YL+rdY9/E54V0lVXyw9wgb9pezYX85JcfqTvk+O0uq2FlSBZz+41onCg8zJEfDoO2r6JsUS1pSDH2SYkmNjyEuOoKYyDCiI8KJiQwjJjI8cPMvO5fR4NZamqz/I431jU3U1DdS2+C/1TQ0UlP/ydejtQ2UVtZTdLSWgiM17D9Szf7DNdQ0NJ79ic4iNjKc8Rk9GJfRg8G9ugd6EIsvPpqIkz7Glpu7nyG++HN+ThGvKehFQowxhtT4GFLjY5g4qGer2zQ2Wcqq6jhUUcvOkkryiyvZUeT/2tqVy078dyU1UHKaWdc6q5S4aIb1jmdY73gyU+NITYgmJS6annHRXeLohsiZKOhFOqHwsE/eDIxOTzrr9rUNjVTWHaey9ji576wmM3sMhwLnrA+W11ByrI7SqnqO1jRQUdNAeXU9HXXV34SYCJK7R5HULYrk7v5bz+5R+BJiSImPJjE2svmWEBNBQmykJpER+RQU9CJdwMeH4VPiohmQGM7krBSvSxKRDqK3xSIiIg5T0IuIiDhMQS8iIuIwBb2IiIjDFPQiIiIOU9CLiIg4TEEvIiLiMAW9iIiIwxT0IiIiDlPQi4iIOExBLyIi4jAFvYiIiMMU9CIiIg5T0IuIiDhMQS8iIuIwBb2IiIjDFPQiIiIOU9CLiIg4TEEvIiLiMAW9iIiIwxT0IiIiDlPQi4iIOExBLyIi4jAFvYiIiMMU9CIiIg5T0IuIiDhMQS8iIuIwBb2IiIjDFPQiIiIOU9CLiIg4TEEvIiLiMAW9iIiIwxT0IiIiDvMk6I0xycaYxcaYvMDXHmfYNsEYU2CM+XVH1igiIuICr/bo7wGWWmuzgKWBx6fzELCiQ6oSERFxjFdBfxXwQuD+C8DVrW1kjBkP+IB/dlBdIiIiTvEq6H3W2kKAwNfUkzcwxoQBvwB+0MG1iYiIOMNYa9vnGxuzBOjdyqr7gBestUknbHvEWtviPL0x5g6gm7X258aYW4Aca+0dp3mubwLfBPD5fOPnzZvXRv+L0FBZWUlcXJzXZYQ89Sk46lNw1KfgqVfBac8+TZs27QNrbU5r69ot6M/EGLMdmGqtLTTGpAG51tqhJ23zR2AK0ATEAVHAb6y1ZzqfT05Ojl27dm07Ve6N3Nxcpk6d6nUZIU99Co76FBz1KXjqVXDas0/GmNMGfUS7POPZLQBuBuYGvr528gbW2ps+vn/CHv0ZQ15ERERa8uoc/VxghjEmD5gReIwxJscY84xHNYmIiDjHkz16a20ZML2V5WuBb7Sy/Hng+XYvTERExDGaGU9ERMRhCnoRERGHKehFREQcpqAXERFxmIJeRETEYQp6ERERhynoRUREHKagFxERcZiCXkRExGEKehEREYcp6EVERBymoBcREXGYgl5ERMRhCnoRERGHKehFREQcpqAXERFxmIJeRETEYQp6ERERhynoRUREHKagFxERcZiCXkRExGEKehEREYcp6EVERBymoBcREXGYgl5ERMRhCnoRERGHKehFREQcpqAXERFxmIJeRETEYQp6ERERhynoRUREHKagFxERcZiCXkRExGEKehEREYcp6EVERBymoBcREXGYgl5ERMRhCnoRERGHKehFREQcpqAXERFxmIJeRETEYQp6ERERhynoRUREHKagFxERcZiCXkRExGEKehEREYcp6EVERBymoBcREXGYgl5ERMRhCnoRERGHKehFREQcpqAXERFxmIJeRETEYQp6ERERhynoRUREHKagFxERcZiCXkRExGEKehEREYcp6EVERBymoBcREXGYgl5ERMRhCnoRERGHKehFREQcpqAXERFxmLHWel1DmzLGlAB7va6jjaUApV4X0QmoT8FRn4KjPgVPvQpOe/Ypw1rbq7UVzgW9i4wxa621OV7XEerUp+CoT8FRn4KnXgXHqz7p0L2IiIjDFPQiIiIOU9B3Dv/jdQGdhPoUHPUpOOpT8NSr4HjSJ52jFxERcZj26EVERBymoA8hxphZxpjtxph8Y8w9rayPNsb8ObB+jTFmQMdX6b0g+vQ9Y8wWY8yHxpilxpgML+r02tn6dMJ21xljrDGmS46aDqZPxpgvBX6nPjLG/KmjawwFQfzd9TfGLDfGrA/87V3uRZ1eM8Y8Z4wpNsZsPs16Y4x5ItDHD40x49q9KGutbiFwA8KBncAgIArYCGSftM13gKcD928A/ux13SHap2lAt8D9b6tPrfcpsF088BawGsjxuu5Q7BOQBawHegQep3pdd4j26X+AbwfuZwN7vK7bo15dDIwDNp9m/eXA/wEGmAisae+atEcfOi4A8q21u6y19cA84KqTtrkKeCFw/2/AdGOM6cAaQ8FZ+2StXW6trQ48XA2kd3CNoSCY3yeAh4CfA7UdWVwICaZPtwNPWWuPAFhrizu4xlAQTJ8skBC4nwgc7MD6Qoa19i3g8Bk2uQp40fqtBpKMMWntWZOCPnT0Bfaf8PhAYFmr21hrjwMVQM8OqS50BNOnE92G/91zV3PWPhljxgL9rLWvd2RhISaY36chwBBjzDvGmNXGmFkdVl3oCKZPDwBzjDEHgDeAOzumtE7n076GnbOI9vzm8qm0tmd+8kcigtnGdUH3wBgzB8gBLmnXikLTGftkjAkDHgdu6aiCQlQwv08R+A/fT8V/dOhtY8xIa215O9cWSoLp043A89baXxhjLgT+EOhTU/uX16l0+Ou49uhDxwGg3wmP0zn10FfzNsaYCPyHx850iMhFwfQJY8ylwH3AF6y1dR1UWyg5W5/igZFArjFmD/5zhQu64IC8YP/uXrPWNlhrdwPb8Qd/VxJMn24D/gJgrV0FxOCf211aCuo1rC0p6EPH+0CWMWagMSYK/2C7BSdtswC4OXD/OmCZDYzu6ELO2qfAIenf4Q/5rng+Fc7SJ2tthbU2xVo7wFo7AP9Yhi9Ya9d6U65ngvm7exX/AE+MMSn4D+Xv6tAqvRdMn/YB0wGMMcPxB31Jh1bZOSwAvhYYfT8RqLDWFrbnE+rQfYiw1h43xtwBvIl/hOtz1tqPjDEPAmuttQuAZ/EfDsvHvyd/g3cVeyPIPj0KxAF/DYxV3Get/YJnRXsgyD51eUH26U1gpjFmC9AI/MBaW+Zd1R0vyD7dDfzeGPNd/Ieib+mCOyIYY17Gf5onJTBe4T+BSABr7dP4xy9cDuQD1cCt7V5TF/w5iIiIdBk6dC8iIuIwBb2IiIjDFPQiIiIOU9CLiIg4TEEvIiLiMAW9iIiIwxT0IiIiDlPQi8g5McacH7iudowxpnvgmu0jva5LRPw0YY6InDNjzE/wT3kaCxyw1j7scUkiEqCgF5FzFpj//H3817WfZK1t9LgkEQnQoXsRaQvJ+K8vEI9/z15EQoT26EXknBljFgDzgIFAmrX2Do9LEpEAXb1ORM6JMeZrwHFr7Z+MMeHAu8aYz1lrl3ldm4hoj15ERMRpOkcvIiLiMAW9iIiIwxT0IiIiDlPQi4iIOExBLyIi4jAFvYiIiMMU9CIiIg5T0IuIiDjs/wEW+S8Whv167wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Opazuj vezane in nevezane orbite\n",
    "\n",
    "plt.figure(figsize=(8,8))\n",
    "\n",
    "plt.plot(r[:,0], r[:,1], label='orbita lažjega telesa (Zemlja)')\n",
    "plt.plot(0,0, '.', color='orange', ms=20, label='masivno telo (Sonce)')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.title('Keplerjev problem')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.grid()\n",
    "\n",
    "plt.axis('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kepler(y, t):\n",
    "    '''\n",
    "    Parametri:\n",
    "        y : zajema podatka o legi in hitrosti [x, y, vx, vy].\n",
    "        t : čas\n",
    "    Izhod:\n",
    "        dydt : časovni odvodi (\"diference\") [dx/dt, dy/dt, dvx/dt, dvy/dt]\n",
    "    '''\n",
    "    \n",
    "    x1, y1 = y[0], y[1]\n",
    "    vx1, vy1 = y[2], y[3]\n",
    "    \n",
    "    R = distance([0, 0], [x1, y1])\n",
    "    dydt = [vx1, vy1, -k*x1/R**3, -k*y1/R**3]\n",
    "    \n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.linspace(0, 5, 1001)\n",
    "y0 = [1, 0, 0, .4]\n",
    "\n",
    "sol = odeint(kepler, y0, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "\n",
    "plt.plot(sol[:,0], sol[:,1], label='orbita lažjega telesa (Zemlja)')\n",
    "plt.plot(0,0, '.', color='orange', ms=20, label='masivno telo (Sonce)')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.title('Keplerjev problem')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.grid()\n",
    "\n",
    "plt.axis('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Klasični model vodikovega iona $H_2^+$\n",
    "\n",
    "Obravnavajmo problem treh teles, pri čemer bomo predpostavili, da sta protona $p$ nepremična, elektron $e$ pa se giblje v njunem potencialu.\n",
    "\n",
    "<img src=\"imgs/hydrogenIon.png\" style=\"width:200px; margin:0px;\">\n",
    "\n",
    "$$F_j = - \\frac{e_0^2}{4\\pi\\epsilon_0 r_j^2} \\frac{\\mathbf{r}_j}{r_j},\\qquad j=1,2$$\n",
    "\n",
    "$\\mathbf{F}_j$ ... elektrostatična sila $j$-tega protona $p$ na elektron $e$<br>\n",
    "$\\mathbf{r}_j$ ... vektor od lege $j$-tega protona $p$ do lege elektrona $e$\n",
    "\n",
    "Poskusi najti kakšno zanimivo orbito, npr. podobno tisti, ki si jo videl/a na [spletni strani](http://zimskasola.fmf.uni-lj.si/#delavnice)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}