{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Plot planetary motion \n", "$$\n", "m\\frac{ d}{dt}\\vec r(t) = -GmM \\frac{\\hat r}{r^2}\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAGDCAYAAAAMFuaDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XeYVPXZxvHvQ1eQIh1EEUTFBgoq\nEqOA3RjBFsUS0CgxRmMssUcTu8ao8Y0xGDVgA0uiEjUqwiI2iKCIoiIgFoooIgIKSHneP35nnNll\ndncWzsyZ3b0/13WuU2fOPbOz88xpv2PujoiIyKaqk3QAERGpGVRQREQkFiooIiISCxUUERGJhQqK\niIjEQgVFRERioYIiiTKzCWZ2etI5qgszO8nMXkg6R1WZ2WVmdk/SOSS/VFAk78zsYzNbaWYrzGyR\nmf3TzJokkOHAAq5vhJm5mR1ZZvrt0fShOTxH52jZeqlp7v6Qux+ch8ixMbN+ZjYvc5q7X+/u+uFQ\nw6mgSKH81N2bAHsAewJXJJynSsys7kY87ENgSMZz1AOOA+bElUukmKigSEG5+3zgv8AuZeeZWVcz\nG29mX5nZYjN7yMyaZ8z/2MwuNLPpZvaNmT1iZo0y5h9hZtPMbKmZvWZmu0XTHwC2Bv4TbSVdFE1/\nzMw+j55ropntnPFcI8zsLjN71sy+Bc6Ptq7qZSxzjJlNq+Dl/gf4kZm1iMYPBaYDn2c8Rx0zu8LM\nPjGzL8zsfjNrFs2eGPWXRrn3MbOhZvZKxuP7mtkb0Wt4w8z6ZsybYGbXmNmrZrbczF4ws1bZgqa2\nKszsoijHQjMbZGaHm9mHZrbEzC7LWL5htLW1IOpuj6Y1Jvx9O0SZV5hZBzP7g5k9mPH4I81sRvS3\nmmBm3XP9O0vxUkGRgjKzTsDhwFvZZgM3AB2A7kAn4A9llvkZ4Yt5W2A3YGj0vHsA9wG/BFoCw4Ex\nZtbQ3U8BPiXaSnL3m6Pn+i/QDWgDvAk8VGZdJwLXAVsA/wd8BRyUMf9k4IEKXu4qYAxwQjT+c+D+\nMssMjbr+QBegCfDXaN5+Ub95lPv1zAea2ZbAM8Ad0Wu+FXjGzFqWeQ2nRq+xAXBhBXnbAY2AjsCV\nwD+i19gL+DFwpZl1iZa9HOgD9AR6AHsBV7j7t8BhwIIocxN3X1Am9/bAKOC3QGvgWUKxb5CxWNa/\nsxS3WldQzOy+6BfYuzksu5+ZvWlma83s2DLznot+XT2dv7Q1ypNmthR4BXgJuL7sAu4+293Huvtq\nd/+S8AW5f5nF7nD3Be6+hLAF0DOafgYw3N0nu/s6dx8JrCZ86WXl7ve5+3J3X00oXD0ytg4AnnL3\nV919vbuvAkYSvmBTX+aHAA9X8rrvB34ePe/+wJNl5p8E3OruH7n7CuBS4ITMLaEK/ASY5e4PuPta\ndx8FfAD8NGOZf7r7h+6+EniU9PuVzRrgOndfA4wGWgF/id6jGcAMwpd7KvfV7v5F9Lf6I3BKDpkB\njgeeif7Wa4BbgM2AvhnLlPd3liJW6woKMILwyycXnxJ+GWX70vgTuf8DCQxy9+buvo27nxV9wZVi\nZm3MbLSZzTezZcCDhC+1TJ9nDH9H+EUPsA1wQVTkl0bFqxNha2cDZlbXzG40sznRuj6OZmWu77My\nD3sQ+KmFEwp+Brzs7gsretHu/grhV/gVwNNZXncH4JOM8U+AekDbip63nMemHt8xY7y89yubr9x9\nXTScyrkoY/7KjMdny531vc6i1GPdfT3hvd7Y3FIkal1BcfeJwJLMadG+++fMbKqZvWxmO0bLfuzu\n04H1WZ5nHLC8IKFrjxsAB3Zz96aErQHL8bGfEX5dN8/oNo9+tRM9b6YTgYHAgUAzoHM0PXN9pR4T\nHf95HTiK8GOiot1dmR4ELmDD3V0ACwjFMGVrYC3hi7yypsDLPjb1+Pk55toU2XKndm1VKbeZGaH4\nFyK35FGtKyjluBs4x917EfYx/y3hPLXVFsAKwkHojsDvqvDYfwBnmtneFjQ2s5+Y2RbR/EWEYxSZ\n61pNOC6yOVl2wZXjfuAiYFfgiRwfcwfh2MvELPNGAeeZ2bbRls/1wCPuvhb4kvBjpkuWx0E49rC9\nmZ1oZvXM7HhgJ6AQu2FHAVeYWevoQP+VhMIJ4b1uWWb3YaZHgZ+Y2QFmVp9QbFcDr+U7tORXrS8o\n0T9xX+Cx6Iyd4UD7ZFPVWn8knFb8DeFg879zfaC7TyEcR/kr8DUwm9IHcm8gfAEuNbMLCYXhE8Kv\n4veASTmu6gnCr+snogPQuWRb4u7jPPvNh+4jbOlMBOYSDuSfEz3uO8JJAa9GuUsdD3L3r4AjCF/I\nXxEK3RHuvjjH17IprgWmEM5ae4dwUsO1Ua4PCAXnoyh3qV1h7j6TsPX5f8BiwjGfn7r79wXILXlk\ntfEGW2bWmbA/exczawrMdPdyi4iZjYiWf7zM9H7Ahe5+RP7SSrExsznAL939xaSziBSTWr+F4u7L\ngLlmdhyE/blm1iPhWFKkzOwYwjGC8UlnESk2tW4LxcxGAf0IZ/MsAq4ifDncRdjVVR8Y7e5Xm9me\nhF0cLQi7Ij53952j53kZ2JFw9slXwC/c/fnCvhopJDObQDhGcYr+1iIbqnUFRURE8qPW7/ISEZF4\nqKCIiEgscmneocZo1aqVd+7cudz53377LY0bNy5coBgoc2Eoc2Eoc2FUNfPUqVMXu3vrShd091rT\n9erVyytSUlJS4fxipMyFocyFocyFUdXMwBTP4TtWu7xERCQWKigiIhILFRQREYmFCoqIiMRCBUVE\nRGKhgiIiIrFQQRERkViooIiISCxUUEREJBYqKCIiEgsVFBERiYUKioiIxEIFRUREYqGCIiIisVBB\nERGRWKigiIhILFRQREQkFiooIiISCxUUERGJhQqKiIjEQgVFRERioYIiIiKxSLSgmNmhZjbTzGab\n2SVZ5t9mZtOi7kMzW5oxb13GvDGFTS4iImXVS2rFZlYXuBM4CJgHvGFmY9z9vdQy7n5exvLnALtn\nPMVKd+9ZqLwiIlKxJLdQ9gJmu/tH7v49MBoYWMHyg4FRBUkmIiJVlmRB6Qh8ljE+L5q2ATPbBtgW\nGJ8xuZGZTTGzSWY2KH8xRUQkF+buyazY7DjgEHc/PRo/BdjL3c/JsuzFwFaZ88ysg7svMLMuhEJz\ngLvPyfLYYcAwgLZt2/YaPXp0uZlWrFhBkyZNNvGVFZYyF4YyF4YyF0ZVM/fv33+qu/eudEF3T6QD\n9gGezxi/FLi0nGXfAvpW8FwjgGMrW2evXr28IiUlJRXOL0bKXBjKXBjKXBhVzQxM8Ry+15Pc5fUG\n0M3MtjWzBsAJwAZna5nZDkAL4PWMaS3MrGE03Ar4EfBe2ceKiEjhJHaWl7uvNbOzgeeBusB97j7D\nzK4mVMNUcRkMjI6qZEp3YLiZrSccB7rRM84OExGRwkusoAC4+7PAs2WmXVlm/A9ZHvcasGtew4mI\nSJXoSnkREYmFCoqIiMRCBUVERGKhgiIiIrFQQRERkViooIiISCxUUEREJBYqKCIiEgsVFBERiYUK\nioiIxEIFRUREYqGCIiIisVBBERGRWKigiIhILFRQREQkFiooIiISCxUUERGJhQqKiIjEQgVFRERi\noYIiIiKxUEEREZFYqKCIiEgsVFBERCQWKigiIhILFRQREYmFCoqIiMRCBUVERGKhgiIiIrGol3QA\nkWrrm2/ggw/g/fdDf9Ys+PRTWLgwdOvX52/dTZtCu3bQqRNstx3suCPssEPottkG6tbN37pFyqGC\nIlLWypUweTK8/jpMmgSvv06/L79MOlVpy5aF7sMPYdy4rIv0KzuhVy/o0yd0P/4xbL01mOU7qdQi\nKihSO337LYwdC//9b+g++6zqz7H55tC9e9g66N4dunULWwcdOkCbNtCwYfy5AdzD1tHChWGLaPZs\nmDkzbCXNnBmmZTN1aujuvDP7/AMPhEMPhSOPDK9FpIpUUKRmW7cOSkrg4YdDt3p1bo/r0wf69g39\nffZhwuzZ9OvXL69Rc2YGzZuHrnt3OOSQrItNmDAhnXnVKnjrrbDFNWkSTJwIn39e+gEvvhi6Cy8s\nPX3vveGEE+BnPwvFUqQcKihSc6xbB88+G36BP/985csfeCAcfnjott++4t0/s2fHlzMJjRrBPvuE\nLpsVK2DChLC19tRTMH9+et7kyaE777z0tLZt4ayz4IwzoH37vEaX6kNneUn1NXcunHlmKARmUK9e\n2F1TtpjstRfcfjssWBB2F6W6sWPDl+QOO+hYQpMmcMQRoRjPm5d+j9asCe/TaaeV3oW3aBFcdVXY\nYkm9/716wb/+ld+TEaSoqaBI9TFjBgwenP4C69IFhg8vvUzfvmHX1qpV6S/FyZPh3HP1S3pj1KsX\ntuTuvbf0e/rOO2ELJbMQv/kmHHtsOMPMLBz0Hz48992MUu2poEjx+vJLOPvsdAHZZRcYPTo9v359\n+MMfwnKpL7pXXw1FJ18HxCXYZZewNbN+fXjf16+HZ54JBT3ls8/CFmSjRum/3wsvJJdZ8i7RgmJm\nh5rZTDObbWaXZJk/1My+NLNpUXd6xrwhZjYr6oYUNrnkhXvYumjbNnwBtWlT+oykLbeE224L+/vd\n4fvvw26XVq2SyyyBWTgW9eqr6eI+aVLYBZkyY0Y4gSD1A+HUUzc8MUCqtcQKipnVBe4EDgN2Agab\n2U5ZFn3E3XtG3T3RY7cErgL2BvYCrjKzFgWKLnFatQp+//vwBVOnDpx0EnzxRXr+RRfB11+HL6iv\nvoLf/hYaN04ur+Ru773DAf7MLZgdd0zPHzEi7IY0g969Ydq0xKJKPJLcQtkLmO3uH7n798BoYGCO\njz0EGOvuS9z9a2AscGieckrM6i1bFg7ymsFmm8G116Zn9ukTLihM/cq96aZweqxUb6ktmPffD3/X\n5cvDD4mUqVNh993p178/tGwJzz2XXFbZaEkWlI5A5tVk86JpZR1jZtPN7HEz61TFx0qxWLnyh+Mh\n+w4cCP/8Z3reySeHs4bcQzHp0ye5nFIYTZrA1Vent14efhgaNAjzliyBww4LRah9+7DrTKoFc/dk\nVmx2HHCIu58ejZ8C7OXu52Qs0xJY4e6rzexM4GfuPsDMfgc0dPdro+V+D3zn7n/Osp5hwDCAtm3b\n9hqdeVC3jBUrVtCkSZP4XmQBFHXm9evp9MgjdL377g1mfTp4MB///Oesb9QogWBVV9Tvczmqa+aO\nc+fS/ZpraFSmuZvl3brx3lVXsbJjcf12rK7vc1Uy9+/ff6q79650QXdPpAP2AZ7PGL8UuLSC5esC\n30TDg4HhGfOGA4MrW2evXr28IiUlJRXOL0ZFmfm119zr1s284iN0Z5zhvnx5cWauhDIXxgaZn33W\nvV69DT9L553n/v33iWQsq0a8z5UApngO3+tJ7vJ6A+hmZtuaWQPgBGBM5gJmlnnhwJHA+9Hw88DB\nZtYiOhh/cDRNkrJyJQwZEnZT9O0brloHOOig9AWFd98ddnWI5Oqww8LFlevXh2thUm67LewiM4OX\nX04un5SSWEFx97XA2YRC8D7wqLvPMLOrzSx1ruFvzGyGmb0N/AYYGj12CXANoSi9AVwdTZNCe/XV\n8E+9+eZw//3p6ePGhSLywgu6oFA2nVk4kSPVMOaxx6bn7bdfmH/22bB2bXIZJdnrUNz9WXff3t27\nuvt10bQr3X1MNHypu+/s7j3cvb+7f5Dx2Pvcfbuo+2d565A8cIebbw7/xPvum54+bFj6auoBA5LL\nJzVb06bw2GPhc/bSS+npd94ZLnbt1q10W2RSMLpSXnK3YkW4MK1OHbj44vT0iRPDP/fw4bpCXQpr\nv/3CZ2/FChgYXXUwezZstVX4wTN2bLL5ahkVFKncokXhzoBbbJFuOqNHj3ABonu4WZNIkho3hief\nDJ/Hv/0tPf3gg0NhGTkyuWy1iAqKlG/u3HDgs1270AIthCvV160LVzW3bp1sPpFsfvWrUFimTElP\nGzo0FJYbbgjzJC9UUGRDH3yQbs13zZow7c9/Dv+It90WdnmJFLtevcJndu7ccNIIwGWXhc/vlVcm\nm62G0jeDpM2bFwpJ9+7pafffH/4pzz8/uVwim6Jz53DL5y+/TN/a+Jprwmf99tsTjVbTqKBIaOqi\nXbtwnCTl8cdDITnllORyicSpVSv48MPQ2GjnzmHaeeeFwpJ5N0rZaCootdmaNeEixJYtw4F3CGdq\nucMxxySbTSRfmjcPu8Eym86//fZQWB55JLlcNYAKSm11zTXhgPvrr4fxVEN9w4Ylm0ukUNq2DZ/5\nzGbzTzghFJa5c5PLVY3VSzqAFNj48XDAAenxQYPCfcB1oF1qqx49QmF58MH0Lt4uXUJ/zZpwG2TJ\nib5FaoulS8Mvr1QxadAg3LDqiSdUTEQg3EbBHY46Kj2tfn244ILkMlUz+iapDS67DFpk3NBy0iRY\nvTrcUldESvv3v8PtpVNuvTX8GPvoo+QyVRMqKDXZe++lL+YCuOSS8Ats772TzSVS7OrXD/8r//tf\nelrXrtChQ3KZqgEVlJrIPTT7vfPO6WlLlqQLi4jkZs89w//TkVED6AsXhh9pr76abK4ipYJS07z3\nXjgmkrond6pV1sxdXiJSNU89Fa5fSdl3XzU9lIUKSk1y6qnprZLWrcNxksz7RojIxmvePPw4u+ii\nML54MZjRePbsZHMVERWUmiC1GT5iRBgfNSq0BNygQaKxRGqkm24Ku5Aje55xRritg6igVHdtn3uu\n9IHCFSvCxVkikj8tWoStlZ//PIy/8EL4UbdqVbK5EqaCUl25Q8+edL/ppjCeutK9ceNkc4nUJiNH\nMjnz1tebbQbPPJNcnoTpEtDqaP78cEe6lPfeK91CsIgUzMpOnWD9+vQFwkccAX36pJs1qkW0hVLd\nPPtsupi0aMGEF19UMRFJmlnYQ3DNNWF80qT0tFpEBaU6Oe88+MlPwvCVV4YDg3XrJptJRNKuuCI0\nkZ9Sp07Yo1BLaJdXdeAO7dunm5gvKYF+/RKNJCLl6NYtNCpZv34Y32qrcH+hWnBLCG2hFLuVK8Ov\nnFQxWbRIxUSk2NWrF34IduwYxo89tlbcxEsFpZgtXJi+FzbAunXQpk1yeUSkaubNC42zQriJ1267\nJZsnz1RQitW0aenrS/r0Cb921My8SPVz3XXw5JNh+J13wsH6GkrfUMXo+edh993D8K9/XStPPxSp\nUQYODKf3p9TQoqKCUmz+/W849NAw/Le/wV//mmweEYlH9+6lG5isgacVq6AUkwcfTJ8J8uij8Ktf\nJZtHROLVvHk40SalTp1wUWQNoYJSLP7xj/T9rJ9+Go47Ltk8IpIfjRqF04pT6tYNJ9zUALoOpRjc\nfz8MGxaGX3wxfd93EamZ6tULRSR1YXK9emFLpZofW1FBSdqYMTBkSBieMAH23z/ROCJSIKndXamz\nN+vUqfbHVLTLK0kTJ4azPyAUFhUTkdrFDNauLT1ejamgJGXGjHQBGTkSfvrTZPOISDLq1i19H5Vq\nXFRUUJKwZAnssksYvuWW9E16RKR2atgQli5Nj++0U3JZNoEKSqGtXQstW4bhU06BCy5INo+IFIdm\nzWDmzDD8/vvV8rtBBaXQUi2QdugQzu4SEUnZfnt44okwfOut8NRTyeapokQLipkdamYzzWy2mV2S\nZf75ZvaemU03s3Fmtk3GvHVmNi3qxhQ2+UY66KD08Lx5yeUQkeI1aBBcdFF6ePHiZPNUQWIFxczq\nAncChwE7AYPNrOyOw7eA3u6+G/A4cHPGvJXu3jPqjixI6E0xcmS4xgTgu++q9YE3Ecmzm26Cxo3D\ncOvWyWapgiS3UPYCZrv7R+7+PTAaGJi5gLuXuPt30egkYCuqo7lzYejQMDx9Omy2WaJxRKQaWLEi\nPVxNfoAmWVA6Ap9ljM+LppXnF8B/M8YbmdkUM5tkZoPyETAW69ZBly5h+M9/hl13TTaPiFQf332X\nHv7Nb5LLkSPzhK7MNLPjgEPc/fRo/BRgL3c/J8uyJwNnA/u7++poWgd3X2BmXYDxwAHuPifLY4cB\nwwDatm3ba/To0eVmWrFiBU2aNNn0F5dh78GD2ezzz/luq6343wMPxPrckJ/M+abMhaHMhZHvzM3e\neYfdo2Ly+qhRrG7XbpOfs6qZ+/fvP9Xde1e6oLsn0gH7AM9njF8KXJpluQOB94E2FTzXCODYytbZ\nq1cvr0hJSUmF86ts1Cj30JiC+/r18T53JPbMBaDMhaHMhVGQzH37pr9LYlDVzMAUz+F7PcldXm8A\n3cxsWzNrAJwAlDpby8x2B4YDR7r7FxnTW5hZw2i4FfAj4D2KybJlMHhwGH733WqzD1REitCrr6aH\nW7VKLkclEiso7r6WsBvrecIWyKPuPsPMrjaz1FlbfwKaAI+VOT24OzDFzN4GSoAb3b24CkqzZqF/\n1lmw887JZhGR6m/JktD/6qvSBaaIJNrasLs/CzxbZtqVGcMHlvO414DiPbp9993p4TvvTC6HiNQc\nLVrApZfCDTfAvvsWZcvEulI+bqtWwS9/GYY//jjRKCJSw1x/fXq4CM8YVUGJW9euoX/88bDNNhUv\nKyJSVYsWhf6778KnnyabpQwVlDhNngwLFoThUaOSzSIiNVObNulmnIrsR6sKSpz69An98eN1VpeI\n5M8LL6SHi+jHqwpKXDIPxPfvn1wOEakdHnoo9E88MdkcGVRQ4uCePhBfZPs0RaSGyiwkJ52UXI4M\nKihx+O1vQ3+XXaBTp2SziEjtMWtW6D/8cLI5IpUWFDOrY2a7m9lPzGyAmbUtRLBqwx3uuCMMT5qU\nbBYRqV222y49fMQRyeWIlHtho5l1BS4mtKU1C/gSaARsb2bfEZpEGenu6wsRtGilbtPZu3f6/gUi\nIoXyySfhbK9nnkk6SYVbKNcCDwJd3f0Qdz/Z3Y/1cLOrI4FmwCmFCFm03OG228JwSUmyWUSkdtp6\n6/Twsccml4MKtlDcfTBA1Ajj6jKzv3H32/MZrFq45ZbQ794dqlmT2yJSg8ydC9tuC//6V6Ixcjko\n/3qO02qf1H2fJ0xINIaI1HKdO6eHR4xIKkX5BcXM2plZL2Cz6KD8HlHXD9i8YAmLVeYB+DZtkssh\nIgLw2GOhf+qpiUWoqLXhQ4ChhPu4/xlIXfq9DLgsv7GqgX32CX0dOxGRYpB5/GTBAujQoeARKjqG\nMhIYaWbHuHuyO+aKzZo16eF+/RKLISJSyt57hzYFO3ZMpHn7inZ5nWxmVl4xMbOuZrZv/qIVsSuu\nCP3jjks2h4hIpokTE119Rbu8WgLTzGwqMJX0dSjbAfsDi4FL8p6wGN18c+hntt8lIpK0Bg3Sw//+\nNxx9dEFXX+4Wirv/BdgDGAW0Bg6IxucDp7j7Me4+qyApi8nixenh5s2TyyEiks1f/hL6xxxT8FVX\neAtgd19nZlPdfWyhAhW9Cy8M/ct0XoKIFKHf/AbOPTeRVedyHcpkM3vMzA43000+GDky9FPHUURE\nilWBz0LNpaBsD9xNaGZltpldb2bb5zdWkfr++/TwZpsll0NEpCKpLZQBAwq62koLigdjo6ZYTgeG\nAP8zs5fMbJ+8Jywmd90V+oMGJZtDRKQit96ayGorPIYCYGYtgZMJWyiLgHOAMUBP4DFg23wGLCpX\nXhn6N9yQbA4RkYrUydhWWLcO6tYtzGpzWOZ1oCkwyN1/4u7/dve17j4F+Ht+4xWZZctCf8cdk80h\nIpKrq64q2KpyKSg7uPs17j6v7Ax3vykPmYrTqlVJJxARyd0//xn6111XsFXmdAylEEGK3kMPhf7A\ngcnmEBHJxZAhBV+l7imfq/vuC/1f/jLZHCIiuUjgKg8VlFy99lroH3RQsjlERKrq008LsppczvJq\nDZwBdM5c3t1Py1+sIlav0rdMRKQ4bLEFLF8Ol18ODzyQ99XlsoXyFOH+8S8Cz2R0IiJSzFItejz4\nYEFWl8vP7c3d/eK8Jylm8+eHfsuWyeYQEamKU0+Fiwv39Z3LFsrTZnZ43pMUs1R7OAVuxkBEZJO0\nbl3Q1eVSUM4lFJVVZrY86pblO1hRGT8+9Pv3TzaHiEgRq3SXl7tvUYggRe2NN0J/772TzSEiUsRy\nOmXJzI4E9otGJ7j70/mLVITmzAn9bt2SzSEiUsQq3eVlZjcSdnu9F3XnRtNqj5UrQ38LbayJSDX1\n7bd5X0Uux1AOBw5y9/vc/T7g0GjaJjOzQ81sppnNNrMN7k9vZg3N7JFo/mQz65wx79Jo+kwzOySO\nPCIiNdaMGXlfRa5XymfePL1ZHCs2s7rAncBhwE7AYDPbqcxivwC+dvftgNuAm6LH7gScAOxMKHB/\ni56vVvj7S3N4bc7iUtNem7OYv780J6FEIlL0vvwy76vIpaDcALxlZiPMbCQwFbg+hnXvBcx294/c\n/XtgNFC25cWBQHTPXR4HDohuQzwQGO3uq919LjA7er5aYbetmnH2w2/9UFRem7OYsx9+i922iqXW\ni0hNtHRp3leRy1leo8xsArAnYMDF7v55DOvuCHyWMT4PKHsa1Q/LuPtaM/sGaBlNn1TmsR1jyFQt\n9O3air+euDtnP/wW+7ZzXnn5Lf564u707doq6WgiUqySLChmtqO7f2Bme0STUvdD6WBmHdz9zU1c\nd7amMMs2lV/eMrk8NjyB2TBgGEDbtm2ZMGFCuYFWrFiRdX6/qF/RY5OwbztnzJw1HNm1Pt9/9i4T\nPqv8McWgvPe5mClzYShz/Pq0bUujRYt4/7PPWBTlzFfmirZQzid8Ef85yzwHNvWy8XlAp4zxrYAF\n5Swzz8zqEY7fLMnxsSGo+93A3QC9e/f2fv36lRtowoQJVDS/onmF9tqcxbzy8lsc2bU+r3xunDBg\nl2qzhVLZ+1yMlLkwlDkPOneGRYvofuSRdO/bF8hf5nKPobj7sGjwMHfvn9kRz1lebwDdzGxbM2tA\nOMg+pswyY4DUXWKOBcZHN/waA5wQnQW2LdAN+F8MmaqF1DGTv564O0d3a/DD7q+yB+pFRH44GF+A\nZlhyOSj/Wo7TqsTd1wJnA88D7wOPuvsMM7s6upAS4F6gpZnNJmwxXRI9dgbwKOG6mOeAX7v7uk3N\nVF1Mn/dNqWMmqWMq0+d9k3AKIc36AAAgAElEQVQyESk6BSwoFR1DaUc40L2Zme1O+rhFU2DzOFbu\n7s8Cz5aZdmXG8CrguHIeex1QuJslh5Umche0ss7cv+sG0/p2bVVtdnmJSAEtXx76zfJ/FmhFx1AO\nAYYSjk/cmjF9OXBZHjMVn44dQxP2H38M226bdBoRkaorwI/hcguKu48ERprZMe7+r7wnKWZ77hkK\nyhtvqKCIiJSj3GMoZnZyNNjZzM4v2xUoX3HYK7pm8n+15ri/iEiVVbTLq3HUb1KIIEWtT5/Qf/XV\nZHOIiFSFZ708L28q2uU1POr/sXBxilR07jaTJlW8nIhIMXnrrdDfYYeCrC6X5utvNrOmZlbfzMaZ\n2eKM3WG1Q8OGSScQEam6J54I/UGDCrK6XK5DOdjdlwFHEK5Q3x74XV5TiYjIpksVlKOOKsjqciko\n9aP+4cAod1+SxzzFK7WVMn9+sjlERHKVugfKnnsWZHW5FJT/mNkHQG9gnJm1BlblN1YR+tWvQv9v\nf0s2h4hIVdXJ9dZXm7iayhZw90uAfYDe7r4G+JYN71tS8517bujfemvFy4mIFIOFCwu+ykrvh2Jm\n9YFTgP3Cva14Cfh7nnMVn86dQ39V7ds4E5Fq6K67Qv+ccwq2yly2g+4CegF/i7o9omm1l4qKiBS7\na64J/QsuKNgqcykoe7r7EHcfH3WnEu7eWPscf3zo33FHsjlERHK1zTYFW1UuBWWdmf3QvK2ZdQFq\nTVPxpdx4Y+hffHGyOUREKjJzZujXq/SoRqxyWdvvgBIz+4jQhP02wKl5TVWsUsdRRESK2YUXhv6f\ns91wN38qLSjuPs7MugE7EArKB+6+Ou/JilXz5rB0KbzyCuy7b9JpREQ29PTToZ+63KFAcml6pRHw\na+APwJXAr6JptdMjj4T+sccmm0NEJJsFC9LD9euXv1we5LLL637CTbX+LxofDDxAOXdSrPEOPjj0\nFy1KNoeISDa//GXo33BDwVedS0HZwd17ZIyXmNnb+QpULaTu4Pj003DEEUmnERFJS+3uSh1HKaBc\nzvJ6y8z6pEbMbG+gdt8YZMyY0P/pT5PNISKS6fXX08MFPsMLcttC2Rv4uZl9Go1vDbxvZu8A7u67\n5S1dsdpjj/Tw8uWwxRbJZRERSUndu2n8+ERWn0tBOTTvKaqj888P7Xqdcgo8+WTSaUSktvv++/Rw\n//6JRMjltOFPChGk2rnpplBQnnoq6SQiInDaaaF/cnL3PyxMm8Y1Ub160KVLGC7wxUMiIqW4w0MP\nheF7700shgrKpkgdAEvgbAoRkR/88Y+h37MnNGiQWIxcLmw828xaFCJMtdOmDYQm/WH06GSziEjt\nlSooCR2MT8llC6Ud8IaZPWpmh5qlvkEFSDfCNnhwsjlEpHa67bbQb9oUWiT72z+XOzZeAXQD7gWG\nArPM7PrMFohrtW7d0sP33JNcDhGpfdzDGacAs2Ylm4Ucj6G4uwOfR91aoAXwuJndnMds1cen0SU6\nZ5yRbA4RqV1SxaR797ALPmG5HEP5jZlNBW4mXCG/q7v/inAXx2PynK966NQJttoqDF96abJZRKR2\nWLMGbr89DP/vf8lmieSyhdIKONrdD3H3x9x9DYC7rwfUkFXKu++G/o03wsqVyWYRkZqvT9Qi1pFH\nQpMmyWaJ5HIM5cryLm509/fjj1RNNWsWrpoH2HHHZLOISM02cya8+WYYfuKJZLNk0HUocRo5MvQ/\n/RSmTk02i4jUXKkfraNHQ53i+RovniQ1gRk891wY7t072SwiUjPdemt6+Pjjk8uRhQpK3A45JN1s\n9JAhyWYRkZplyRK44IIwPH9+slmyUEHJh2++Cf3774f3dZhJRGLSsmXoX3YZdOiQbJYsEikoZral\nmY01s1lRf4PLO82sp5m9bmYzzGy6mR2fMW+Emc01s2lR17Owr6ASm28ODz8chnfaKVx8JCKyKW7O\nuOzvuuuSy1GBpLZQLgHGuXs3YFw0XtZ3wM/dfWfCPVluN7PmGfN/5+49o25a/iNX0eDBoSkESN+H\nXkRkY3zyCVx8cRheuDDZLBVIqqAMBKJTohgJDCq7gLt/6O6zouEFwBdA64IljMOXX4b+iy/qJlwi\nsnHcoXPnMHzrrdCuXaJxKpJUQWnr7gsBon6FbQaY2V5AA2BOxuTrol1ht5lZw/xF3QQNGsDkyWH4\nqKPgq6+SzSMi1c+uu4Z+69Zw3nnJZqmEeZ7275vZi4SWisu6HBjp7s0zlv3a3bM2k2lm7YEJwBB3\nn5Qx7XNCkbkbmOPuV5fz+GHAMIC2bdv2Gl1BM/MrVqygSR6uOO165510evxxACaUlMT63PnKnE/K\nXBjKXBj5zNz+mWfY4ZZbAJjw4otQt24sz1vVzP3795/q7pVfC+HuBe+AmUD7aLg9MLOc5ZoCbwLH\nVfBc/YCnc1lvr169vCIlJSUVzt8kYcPVvWfPWJ82r5nzRJkLQ5kLI2+Z3347/b3x/vuxPnVVMwNT\nPIfv2KR2eY0BUhdpDAE2uDG7mTUAngDud/fHysxrH/WNcPzl3bymjcPq1aE/bRpcdVWyWUSkuC1f\nDj16hOF77602zTklVVBuBA4ys1nAQdE4ZtbbzFI3FfkZsB8wNMvpwQ+Z2TvAO4TGK68tbPyN0KBB\nupn7q6+GsWOTzSMixck9fYboccfBaaclm6cK6iWxUnf/Cjggy/QpwOnR8IPAg+U8fkBeA+ZLp07w\n/PPhavqDD4Y5c6BLl6RTiUgxadUqPfzoo8nl2Ai6Ur7QDj4Y/vCHMNy1q878EpG0Aw4IzasArF2b\nbJaNoIKShKuugp/9LAy3agWrViWbR0SSd/bZMH58GF62LLYzugpJBSUpjzwCu+wShjfbDNavTzaP\niCTn9tvhzjvD8Pz5sMUWyebZSCooSXrnnfRw3boqKiK10V13pS9YnD69KBt9zJUKStIyi4iKikjt\nMnw4nHVWGH7llfRV8dWUCkrSzFRURGqjf/wDzjwzDE+cCD/6UbJ5YqCCUgyyFZVqeIaHiOToL3+B\nYcPC8EsvwY9/nGyemKigFIuyRaV+ffjuu+TyiEh+XHgh/Pa3YbikBPbbL9k8MUrkwkYpR6qoNG8e\nThts3Bi++CK0Mioi1d+xx8K//hWGp01LN69SQ2gLpdiYhVsI77VXGG/TBmbNSjaTiGy6XXdNF5OP\nP65xxQRUUIrX5Mlw0klhePvt4Zlnks0jIhvn++/DD8V3ozZsFy+GbbZJNlOeqKAUswcfhJtuCsNH\nHAG//32yeUSkahYsgIYZ9/9btQpatkwuT56poBS7iy6CcePC8LXXpneFiUhxe+UV6NgxDO+yS2hF\nuGFx3lw2Lioo1cGAAemm7994I2w+L1+ebCYRKd8VV6RPBb7wwtKtYtRgOsuruujUKWwuN2oUxps2\nDb+ARKR4rFsHW24ZztIEeOIJGDQo2UwFpC2U6qRhw7DZnDpYv+++dL3rrmQziUjw8cdQr166mCxc\nWKuKCaigVE8PPghPPw1Ap0cfDbvAVq5MOJRILXbHHbDttmG4W7dwPVm7dslmSoAKSnX1k5+Eix5T\nNt883A1SRApn9Wr2O/BAOPfcMH7LLfDhh+FHXi2kglKdtW7NhJISOProMH7ooaEZB/dkc4nUBiUl\n0KgRddatC+OffAIXXJBspoSpoNQE//oXvPpqGH75ZahTJ30RlYjEyz38cBswAICv+vQJ07beOuFg\nyVNBqSn69g1X5KYumtp1VzjsMG2tiMTppZfCD7aXXw7jEybwzg03JJupiKig1CT164dmHVLtBT33\nXPjwv/56srlEqru1a0MTSP36hfGdd4Y1a2D//RONVWxUUGqio48OWytduoTxvn2he/cwTUSqZvTo\n8GMt1Ujryy+HXcr1dBlfWSooNVX9+jBnTrrZlg8+CNex3HJLsrlEqotPPw1naw0eHMYPOCCcDrzv\nvsnmKmIqKDXdgAHhn+C448L4734X/kmmT082l0ixWrcuFI/MFoFnzoQXX6y1pwPnSgWlNjCDRx+F\nzz9PT+vRAzp0SF/VKyLwf/8XdmWNHx/Ghw8PJ7Zsv32yuaoJFZTapG3b8M/xn/+E8YULoVmzcJFk\n6lx6kdrouefCD6/f/CaM//jH4UB86r7vkhMVlNroiCNCYbnkkjD+7LPhV9nllyebS6TQZswIheSw\nw9LTFiyAiROhbt3kclVTKii12Q03hFMfU/9M118f/rn+/Odkc4nk20cfhc/6Lrukp731Vvih1b59\ncrmqORWU2q5evbCFsnRp+mZAF14Y/tl0RpjUNHPmhM92167paU8+GQpJz57J5aohVFAkaNYM5s2D\nRYvCsRZInxGWug2xSHU1e3b4LG+3XXraAw+EQjJwYHK5ahgVFCmtTZtwNtgXX6Sb377kkvDPOHRo\n2EUmUl2MHx8+u926pac9/HAoJCefnFyuGkoFRbJr3TqcBfbll+lfdSNHQoMGsPvusGRJsvlEKnLn\nnaGQHHBAetro0aGQpC5UlNipoEjFWrUKTU6sWgXHHBOmTZsWGqE0g0mTks0nkrJyJZx4Yvhcnn12\nevrkyaGQHH98ctlqCRUUyU3DhvD44+Gq+z/+MT19n33CP/BFF+laFknG5MnhM7j55jBqVJi29dYw\nf34oJHvtlWy+WkQFRarGDK68MvyjptoJA/jTn8IZY1ttFZqpEMmntWvTJ4306ZOePmwYrF4dbnbV\noUNy+WqpRAqKmW1pZmPNbFbUb1HOcuvMbFrUjcmYvq2ZTY4e/4iZNShcevnBgAGhsHzzDRx+eJg2\nfz7suGP4Rz/tNPj222QzSs3yn/+Ez1b9+qVPa3/xxfBZHD48HOeTRCS1hXIJMM7duwHjovFsVrp7\nz6g7MmP6TcBt0eO/Bn6R37hSoaZN4Zlnwj/0ww+np//zn9CkSfgCuOce3exLNs6cOaHtOTM4MuNr\nYNCg8GPGvfTBd0lMUgVlIDAyGh4JDMr1gWZmwADg8Y15vOTZ4MHhH3zlSjjrrPT0M84IN/syo824\ncSouUrFPP4VDDklfO5JqHbttW5g6NXx+nngi/JiRopFUQWnr7gsBon6bcpZrZGZTzGySmaWKRktg\nqbuvjcbnAR3zG1eqrFGjcOqme/iFmXEV8k7XXvtDcfnhmgCRjz+mx/nnh8/FNtvACy+k591zTzgh\n5PPPYY89kssoFTLP0z+zmb0ItMsy63JgpLs3z1j2a3ff4DiKmXVw9wVm1gUYDxwALANed/ftomU6\nAc+6+67l5BgGDANo27Ztr9GjR5ebecWKFTRp0iTXl1gUqlvmxh99xA7XXEPTjz/eYN4nJ57Ipyef\nzLrNNit8sEpUt/cZqkfmFlOm0O3229l8/vwN5s3+9a+Zd9RRRd9IY3V4n8uqaub+/ftPdffelS7o\n7gXvgJlA+2i4PTAzh8eMAI4FDFgM1Ium7wM8n8t6e/Xq5RUpKSmpcH4xqtaZp09379HDPWyjlO72\n28/97bcTzZmpWr/PxeTbb91///vsf3PwD885x33duqRTVklRvs+VqGpmYIrn8B2b1C6vMcCQaHgI\n8FTZBcyshZk1jIZbAT8C3oteXAmhuJT7eKkGdt01XCTpHpp6Of309LyJE9MHYs3gl78MbY1J9bJ2\nbWhhoVOn8Hds3BiuuSY9v2tXeP75H0rK/KOPDrtDpVpK6i93I3CQmc0CDorGMbPeZnZPtEx3YIqZ\nvU0oIDe6+3vRvIuB881sNuGYyr0FTS/xa90a/vGP8MWybh3ce2/pA653353+UjKDCy4IB26luKxZ\nEy4u3G239Om9Q4eW/jEwZEho1sc9NNp48MGJxZV4JVJQ3P0rdz/A3btF/SXR9Cnufno0/Jq77+ru\nPaL+vRmP/8jd93L37dz9OHdfncTrkDypUydcw5I6JfTrr+H3vy+9zK23hgO3qQIzYEBohn/9+mQy\n11affAKpA+lm4RqQE0+Ed95JL3PQQaGJntSOrREj0g2PSo2ibUspfs2bw9VXp7+Q5s8PLSDXq5de\npqQk3Mq4bt30l9vRR4frY9auLf+5JXdz5oTC3rFj+j3u3Bluu630cgceGG6pu359+Hu98ALsvXci\nkaWwVFCk+unQIX23SffQf+QR2HPP0ss98US43XH9+ukvwJ12Crc6Tl3LIBv66quwFTFoUPp9S10P\ncu214Ra5KQ0bwmWXpdvNcoexY9PXkEitUq/yRUSKXL168LOfhS5lxYrQmOV998HLL6env/9+6K6/\nvvRzbLEF9O+f7nbdtWYfHF68OJz4MGFC6DJ3UZWnZUs49dRwTGTnnfMcUKojFRSpmZo0CV98Q4em\np7mHLZN//ztsvXzwQXre8uUwZkzosugH4eLM3XdP93v0KL4rtdeuDY1zvv02XcaMCVtyb78d7sSZ\nq/r1w9bJoEGhjbbmzSt/jAgqKFKbmEHv3qEru4WyaBG89FK4w9+ECdlbTJ42LXRV0awZtGgRvpSb\nN08PN22abi0glS3VX78+FLhvvgnd0qXp/tKl4d40Odi6opmNG0O/fqHbf/9w9XmRX0AoxU8FRQRC\nG1Fld5tlmPj88+y35Zbw1luhmzYt9FdXcoJhqigU0g47QI8efLTFFnQ56qiwRdWhg45pSN6poIjk\nYH3DhuGgf9kD/xVZty5saXz9deiWLk33ly1Ln+KcOjkg8ySBpk3D1k3z5qHftGl6OMcmMz6dMIEu\n/frlnldkE6mgiORL3brpXV3bbpt0GpG8q8GnsYiISCGpoIiISCxUUEREJBYqKCIiEgsVFBERiYUK\nioiIxEIFRUREYqGCIiIisVBBERGRWKigiIhILFRQREQkFiooIiISCxUUERGJhQqKiIjEQgVFRERi\noYIiIiKxUEEREZFYqKCIiEgsVFBERCQWKigiIhILFRQREYmFCoqIiMRCBUVERGKhgiIiIrFQQRER\nkViooIiISCxUUEREJBYqKCIiEotECoqZbWlmY81sVtRvkWWZ/mY2LaNbZWaDonkjzGxuxryehX8V\nIiKSKaktlEuAce7eDRgXjZfi7iXu3tPdewIDgO+AFzIW+V1qvrtPK0hqEREpV1IFZSAwMhoeCQyq\nZPljgf+6+3d5TSUiIhstqYLS1t0XAkT9NpUsfwIwqsy068xsupndZmYN8xFSRERyZ+6enyc2exFo\nl2XW5cBId2+esezX7r7BcZRoXntgOtDB3ddkTPscaADcDcxx96vLefwwYBhA27Zte40ePbrczCtW\nrKBJkyY5vLriocyFocyFocyFUdXM/fv3n+ruvStd0N0L3gEzgfbRcHtgZgXLngvcXcH8fsDTuay3\nV69eXpGSkpIK5xcjZS4MZS4MZS6MqmYGpngO37FJ7fIaAwyJhocAT1Ww7GDK7O6KtlAwMyMcf3k3\nDxlFRKQKkiooNwIHmdks4KBoHDPrbWb3pBYys85AJ+ClMo9/yMzeAd4BWgHXFiCziIhUoF4SK3X3\nr4ADskyfApyeMf4x0DHLcgPymU9ERKpOV8qLiEgsVFBERCQWKigiIhILFRQREYmFCoqIiMRCBUVE\nRGKhgiIiIrFQQRERkViooIiISCxUUEREJBYqKCIiEgsVFBERiYUKioiIxEIFRUREYqGCIiIisVBB\nERGRWKigiIhILFRQREQkFiooIiISCxUUERGJhQqKiIjEQgVFRERioYIiIiKxMHdPOkPBmNmXwCcV\nLNIKWFygOHFR5sJQ5sJQ5sKoauZt3L11ZQvVqoJSGTOb4u69k85RFcpcGMpcGMpcGPnKrF1eIiIS\nCxUUERGJhQpKaXcnHWAjKHNhKHNhKHNh5CWzjqGIiEgstIUiIiKxqHUFxcy2NLOxZjYr6rfIskxP\nM3vdzGaY2XQzOz5j3ggzm2tm06KuZx6zHmpmM81stpldkmV+QzN7JJo/2cw6Z8y7NJo+08wOyVfG\njch8vpm9F72v48xsm4x56zLe1zFFlHmomX2Zke30jHlDos/SLDMbUkSZb8vI+6GZLc2YV/D32czu\nM7MvzOzdcuabmd0RvZ7pZrZHxryk3uPKMp8UZZ1uZq+ZWY+MeR+b2TvRezyliDL3M7NvMv7+V2bM\nq/AzlRN3r1UdcDNwSTR8CXBTlmW2B7pFwx2AhUDzaHwEcGwBctYF5gBdgAbA28BOZZY5C/h7NHwC\n8Eg0vFO0fENg2+h56hZJ5v7A5tHwr1KZo/EVCXwecsk8FPhrlsduCXwU9VtEwy2KIXOZ5c8B7kv4\nfd4P2AN4t5z5hwP/BQzoA0xO8j3OMXPfVBbgsFTmaPxjoFURvs/9gKc39TNVXlfrtlCAgcDIaHgk\nMKjsAu7+obvPioYXAF8AlV7UE7O9gNnu/pG7fw+MJmTPlPlaHgcOMDOLpo9299XuPheYHT1f4pnd\nvcTdv4tGJwFbFSBXRXJ5n8tzCDDW3Ze4+9fAWODQPOXMVNXMg4FRBchVLnefCCypYJGBwP0eTAKa\nm1l7knuPK83s7q9FmaA4Psu5vM/l2ZT/gx/UxoLS1t0XAkT9NhUtbGZ7ESr2nIzJ10WbubeZWcM8\n5ewIfJYxPi+alnUZd18LfAO0zPGx+VDV9f6C8Ks0pZGZTTGzSWa2QaHPk1wzHxP9zR83s05VfGzc\ncl5vtEtxW2B8xuQk3ufKlPeaknqPq6rsZ9mBF8xsqpkNSyhTefYxs7fN7L9mtnM0LZb3uV4c6YqN\nmb0ItMsy6/IqPk974AFgiLuvjyZfCnxOKDJ3AxcDV2982vJXn2Va2VPyylsml8fmQ87rNbOTgd7A\n/hmTt3b3BWbWBRhvZu+4+5xsj49RLpn/A4xy99VmdiZhq3BAjo/Nh6qs9wTgcXdflzEtife5MsX2\nWc6ZmfUnFJR9Myb/KHqP2wBjzeyDaOshaW8SmlFZYWaHA08C3Yjpfa6RWyjufqC775KlewpYFBWK\nVMH4IttzmFlT4BngimgTPPXcC6PN8tXAP8nfrqR5QKeM8a2ABeUtY2b1gGaEzd1cHpsPOa3XzA4k\nFPcjo/cR+GH3Iu7+ETAB2D2fYSOVZnb3rzJy/gPoletj86Qq6z2BMru7EnqfK1Pea0rqPc6Jme0G\n3AMMdPevUtMz3uMvgCcozC7nSrn7MndfEQ0/C9Q3s1bE9T4X+qBR0h3wJ0oflL85yzINgHHAb7PM\nax/1DbgduDFPOesRDkBuS/og2c5llvk1pQ/KPxoN70zpg/IfUZiD8rlk3p2w+7BbmektgIbRcCtg\nFhtxUDBPmdtnDB8FTIqGtwTmRtlbRMNbFkPmaLkdCAeHLen3OVpfZ8o/WPwTSh+U/1+S73GOmbcm\nHJ/sW2Z6Y2CLjOHXgEOLJHO71OeBUOQ+jd7znD5Tla67UC+yWDrCMYZx0T/SuNSHk7D75Z5o+GRg\nDTAto+sZzRsPvAO8CzwINMlj1sOBD6Mv4MujaVcTftkDNAIeiz7U/wO6ZDz28uhxM4HDCvj+Vpb5\nRWBRxvs6JpreN3pf3476vyiizDcAM6JsJcCOGY89LXr/ZwOnFkvmaPwPlPnBk9T7TNhKWhj9X80j\n7CI6Ezgzmm/AndHreQfoXQTvcWWZ7wG+zvgsT4mmd4ne37ejz83lRZT57IzP8iQyimG2z1RVO10p\nLyIisaiRx1BERKTwVFBERCQWKigiIhILFRQREYmFCoqISDVXWaOQZZbdz8zeNLO1ZnZsmXnPmdlS\nM3t6Y3KooIgUiJmdaWY/j4aHmlmHjHn3mNlOeV7/bytY/2gz65bP9UtejSD3Ns4+JTR4+nCWeX8C\nTtnYECooIgXi7n939/uj0aGElqxT80539/fyte6oJYXTSH+JlFo/cBdwUb7WL/nlWRqFNLOu0RbH\nVDN72cx2jJb92N2nA+uzPM84YPnG5lBBkVrNzPaMGn1sZGaNLdwDZ5cyy3Q2sw/MbGRGA5GbR/MO\nMLO3ontf3JdqLNTMbrT0fV9uiab9wcwujHYz9AYeiu5JsZmZTTCz3tFyg6Pne9fMbsrIscLMrosa\n9ptkZm2zvJ47Uve4MLNDzGyimdUhtD32pruvzbZ+4GXgwKjwSM1wN3COu/cCLgT+lu8VqqBIrebu\nbwBjgGsJ98p50N2z7YfeAbjb3XcDlgFnmVkjwq6G4919V0LzFb8ysy0JTbTsHC1/bZl1Pg5MAU5y\n957uvjI1L9oNdROhAPQE9sxoEbgxodmXHsBE4IwsOS8Bjo8aLLyDcGX5euBHwNTy1h8tMxvokeU5\npZoxsyaEVhEeM7NpwHCgfb7Xq4IiEposOYjwq/3mcpb5zN1fjYYfJLQsuwMw190/jKaPJNzgaBmw\nCrjHzI4Gviv7ZBXYE5jg7l96uCXBQ9FzAnwPpA6WTiW02VSKh3vNnEG4b8hfPd2KcHvgy0rW/QWl\nd4NJ9VUHWBr9YEh13QuxUpHabkugCbAFoX20bMq2UVRe0+pEhWAv4F+EG7g9V4UsWZ8zssbTbSWt\no/zbT+wKfEXp4rCS8l9bSqNoOanm3H0ZMNfMjoMfbrGc961PFRSRsK/594StgZvKWWZrM9snGh4M\nvAJ8AHQ2s+2i6acAL0W7G5p5aB78t4RdV2UtJxSwsiYD+5tZKzOrG63rpVxfSHRDrQsIrTofZmZ7\nR7PeB7bLWDTb+rcnNBwo1YyZjQJeB3Yws3lm9gvgJOAXZpZqpHJgtOyeZjYPOA4YbmYzMp7nZUKD\nswdEz3NIVXLoAJzUatFptGvd/eHoC/w1Mxvg7uPLLPo+MMTMhhNaqr7L3VeZ2amE/dT1gDeAvxO2\neJ6KjrEYcF6WVY8A/m5mK4FUocLdF5rZpYRWjQ141sN9fHJ5LQbcC1zo4eZOvwBGmNmehKbhH6hg\n/U2BlR7dzVSqF3cfXM6sDU4ljo4bZr1dsbv/eFNyqLVhkUqYWWfgaXffpZJFi5qZPQFc5O6zssw7\nD1jm7vcWPpnUFNrlJVJ7XEL5Z/osJZxUILLRtIUiIiKx0BaKiIjEQgVFRERioYIiIiKxUEEREZFY\nqKCIiEgsVFBERCQW/0bPZ3UAAAAFSURBVA97G4BNUX8AcgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "#work in MKS units\n", "Tmax = 3*10**7 #seconds\n", "Dt = 1000 #seconds\n", "Nmax=int(Tmax/Dt)\n", "G = 6.67384*10**(-11) #m/sec**2\n", "M = 1.99*10**30\n", "t = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the time\n", "x = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the x position\n", "y = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the y position\n", "vx = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the x velocity\n", "vy = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the y velocity\n", "ax = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the x acceleration\n", "ay = np.zeros(Nmax, ) # numpy array holding the Nmax+1 values of the y acceleration\n", "# implementing the leapfrog method we will interpret the positions x[n] and \n", "# accelerations a[n] as their values at the time Dt*n while the velocity v[n] will be \n", "# treated as the velocity at the time (n+1/2)*Dt.\n", "t[0]=0\n", "x[0]=1.48*10**11 #meters\n", "y[0]=0\n", "vx[0]=0\n", "vy[0]=3.0*10**4/2 #meters/sec\n", "for n in range(0,Nmax-1):\n", " t[n] = Dt*n\n", " x[n+1]=x[n]+vx[n]*Dt\n", " y[n+1]=y[n]+vy[n]*Dt\n", " ax[n+1]=-(G*M)/((x[n+1])**2+(y[n+1])**2)*(x[n+1]/np.sqrt((x[n+1])**2+(y[n+1])**2))\n", " ay[n+1]=-(G*M)/((x[n+1])**2+(y[n+1])**2)*(y[n+1]/np.sqrt((x[n+1])**2+(y[n+1])**2))\n", " vx[n+1]=ax[n+1]*Dt+vx[n]\n", " vy[n+1]=ay[n+1]*Dt+vy[n]\n", "Figure = plt.figure(figsize=(6,6)) #create a fig object named \"fig\"\n", "Axis = Figure.add_subplot(111) #create an axis object within \"Figure\" named \"Axis\"\n", "Axis.plot(x,y,'r-') #Add a curve described by the arrays x and y to Plot, choose a red solid curve.\n", "Axis.plot(0.0,0.0,'x') #Locate the position of the sun.\n", "Axis.set(xlabel='x position x(t)', ylabel='y position y(t)', # add a tile and axis labels.\n", " title='Planetary Motion') \n", "Axis.set_aspect('equal', 'datalim') #set the scale of the x and y axes to be the same\n", "Axis.grid() #superimpose a grid on Plot\n", "#fig.savefig(\"test.png\") # save the figure as a png file.\n", "plt.show(Figure) # Plot the figure here." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "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.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }