Skip to content

Instantly share code, notes, and snippets.

@hannorein
Created November 10, 2019 21:39
Show Gist options
  • Save hannorein/400555b439abfc0a62ef7e20a8c26221 to your computer and use it in GitHub Desktop.
Save hannorein/400555b439abfc0a62ef7e20a8c26221 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"import rebound\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Searching NASA Horizons for 'Sun'... Found: Sun (10).\n",
"Searching NASA Horizons for 'Mercury'... Found: Mercury Barycenter (199).\n",
"Searching NASA Horizons for 'Earth'... Found: Earth-Moon Barycenter (3).\n"
]
}
],
"source": [
"sim = rebound.Simulation()\n",
"sim.add([\"Sun\",\"Mercury\",\"Earth\"])"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADoCAYAAABM+DfFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3xU5b3v8c8vk4RAbiTEEAgESAgECDfJFqnQVkUB6xa3t6NV8dLKUbd7W+05FVu3FltPa7vVbq3WetSttd6xFxSQCojijVtAIIRLLgQSkgBJSCAht8lz/siEEyEhCVkzz2Tm93698mJmzZpnfbMm/OaZZ9ZajxhjUEop5RshtgMopVQw0aKrlFI+pEVXKaV8SIuuUkr5kBZdpZTyIS26SinlQ1aLrojMFZHdIpInIos6eDxFRD4WkS0isk1ELrORUymlnCK2jtMVERewB7gEKAY2AjcYY3a2W+cFYIsx5g8iMh5YbowZeaZ2ExISzMiRZ1wlYDQ1NdHU1ER4eDihoaG24ygV0DZv3nzEGHNOb9ux+T/1PCDPGFMAICJvAfOBne3WMUCM53YscLCrRkeOHMmmTZscjupfmpqaKCgowO12k5aWRr9+/WxHUkGqpaUFEUFEbEfxOhEpcqIdm0U3GTjQ7n4xMP2UdX4O/ENE/g2IBGZ31JCILAQWAqSkpDge1J8cPXqUwsJCEhMTGTp0aFD8sSv/FRISgjEmqIpvb/n7F2k3AK8YY4YBlwGvichpmY0xLxhjsowxWeec0+vev18yxlBUVMS+ffsYPXo0ycnJ+geu/IKInCy+elmBrtns6ZYAw9vdH+ZZ1t4PgLkAxpgvRSQCSAAO+SShn2hqamLPnj24XC4yMzN1/Fb5pbbCq87MZk93I5AuIqNEJBy4Hlh6yjr7gYsBRGQcEAEc9mlKy44fP8727duJjY0lIyNDC67ya/rpq2vW/gcbY5pF5B5gJeACXjbG5IjIo8AmY8xS4MfA/xWR+2j9Uu1WE0RvpYcOHeLAgQOkpqYSFxdnO45SPeZ2u1mxYgVbtmxh6tSpzJs3D5fLZTuWVVa7TcaY5cDyU5Y93O72TuACX+eyzRjDvn37qKmpYfz48fTv3992JKV6zO12M2fOHNavX09tbS2RkZFMnz6dlStXBnXh9fcv0oKO2+0mNzeXxsZGMjMzteCqPmvFihWsX7+e48ePY4zh+PHjrF+/nhUrVtiOZpUWXT/S2NjIjh07iIiIYMyYMUHdG1B935YtW6itrf3GstraWrZu3WopkX/Qousn6urq2L59OwkJCaSmpuoXEqrPmzp1KpGRkactj4+Pt5DGf2jR9QM1NTXk5OSQkpJCcnKy7ThKOWLevHlMnz6dqKgoRISoqCjGjRvHz3/+cz777DPb8azR448sq6iooLCwkDFjxhAbG2s7jlKOcblcrFy5khUrVrB161amTJnCvHnzWLNmDVdddRUvvPACV155pe2YPmftgjfekpWVZfrKtReOHDlCQUEBEyZM6PBjmFKBKjs7m4ULF3L//ffz/e9/33acbhGRzcaYrN62oz1dSw4dOsS+ffu04KqgdO655/Lqq68yZ84cmpubWbBgge1IPqNF14KysjL279/PxIkT9ZAwFbQmTJjARx99xJw5c3C73dx22222I/mEFl0fKy0tpbi4WAuuUsC4ceNYtWoVl1xyCc3Nzdxxxx22I3mdFl0fKi0tpaSkhIkTJxIREWE7jlJ+YcyYMaxevZpLLrmEfv36BfxQgxZdHzl8+DBFRUVMnjxZC65Spxg9ejSrVq1i9uzZhIWFccMNN9iO5DVadH2gsrKSvLw8Jk2apEMKSnUiLS2NDz74gEsvvZS4uDjmzp1rO5JX6MkRXlZdXc3u3bv1KAWlumHChAm899573H777Xz55Ze243iFFl0vOn78ODt37iQjI4OYmJiun6CU4vzzz+eVV17hmmuuYceOHbbjOE6LrpfU19ezY8cORo8erdfCVaqHLr30Up566imuvfZaioocmQ/Sb+iYrhe43W527NhBcnIygTpnm1Ledt1111FVVcWVV17JJ598EjCfFrWn6zBjDDt37iQmJobhw4d3/QSlVKcWLlzIzJkzufnmm3G73bbjOEKLrsPy8/NpaWkhPT3ddhSl+jwR4cknn6S+vp4HH3zQdhxHaNF10MGDB6msrGTChAl6PVylHBIWFsabb77JsmXLePnll23H6TWrRVdE5orIbhHJE5FFnaxznYjsFJEcEXnD1xm76+jRoyevp6Az9irlrPj4eP72t7/x8MMPs27dOttxesVa0RURF/AsMA8YD9wgIuNPWScdeBC4wBgzAfiRz4N2Q2NjIzk5OYwZM0ZPflDKS9LT03nttddYuHAhpaWltuOcNZs93fOAPGNMgTGmEXgLmH/KOncAzxpjqgCMMYd8nLFLxhhycnIYOnRo0E9DopS3XXjhhSxYsICbb76Z5uZm23HOis2imwwcaHe/2LOsvTHAGBH5XES+EpEOzwsUkYUisklENh0+fNhLcTtWUFCAiDBy5EifblepYPXAAw/Qv39/Hn30UdtRzoq/f5EWCqQD3wVuAP6viAw8dSVjzAvGmCxjTJYvj4s9cuQI5eXl+sWZUj4UEhLCf//3f/PGG2/w4Ycf2o7TYzaLbgnQ/kDWYZ5l7RUDS40xTcaYQmAPrUXYuvr6+pPXVAgLC7MdR6mgkpCQwJ/+9CcWLlzIgQMHun6CH7FZdDcC6SIySkTCgeuBpaes8zdae7mISAKtww0FvgzZkbZx3OTkZJ1MUilLvvWtb3HfffexePHiPnXihLWia4xpBu4BVgK5wDvGmBwReVRErvCsthKoEJGdwMfA/zbGVNhJ/P/t378fYwwjRoywHUWpoHbvvfdSXFzM008/bTtKt+lswD10/PhxNm/ezHnnnaeHhynlB4qKipg1axYrVqxgwoQJXtuOU7MB+/sXaX6lpaWFnJwc0tPTteAq5SdGjBjBL37xC374wx/S2NhoO06XtOj2QEFBAREREQwdOtR2FKVUOwsWLGDo0KH8+te/th2lS1p0u6m6upqysjLGjRtnO4pS6hQiwrPPPsvLL7+MN4cXnaBFtxvaLteYlpZGeHi47ThKqQ4kJSXxu9/9jh//+Mc0NTXZjtMpLbrdsH//fsLCwhgyZIjtKEqpM5g/fz4DBw7kueeesx2lU1p0u1BfX09BQQHjx4/vemWllFVt19998sknKS4uth2nQ1p0u7Br1y5SUlIYMGCA7ShKqW5IS0vjzjvv5Cc/+YntKB3SonsGhw4dora2llGjRtmOopTqgfvuu4+cnBxWrlxpO8pptOh2oqWlhV27dpGRkUFIiO4mpfqSiIgInnrqKX72s59x4sQJ23G+QatJJ4qKioiKimLQoEG2oyilzsJFF13E1KlTef75521H+QYtuh1obGyksLCQsWPH2o6ilOqFRYsW8fTTT3PkyBHbUU7SotuB/Px8kpKSiIyMtB1FKdULaWlpXHPNNfz2t7+1HeUkLbqnqK2t5eDBg4wePdp2FKWUAx544AHeeecdCgqsXxUW0KJ7mj179jBq1Cg980ypAJGQkMA999zD4sWLbUcBtOh+Q3V1NUePHtXr5CoVYO6++242bNjAli1bbEfRotve3r17SU1NxeVy2Y6ilHJQ//79+elPf8pjjz1mO4oW3TbV1dXU1NQwfPjwrldWSvU51113Hbm5udZ7u1p0PfLy8khNTdUTIZQKUP369ePee+/lP//zP63m0AoD1NTUcPToUe3lKhXgFixYwKZNm8jJybGWwWrRFZG5IrJbRPJEZNEZ1rtaRIyI9Hp+oo7oWK5SwSEiIoJ77rmHJ554wloGa0VXRFzAs8A8YDxwg4icdv1EEYkG7gXWeyNHbW0tVVVVpKSkeKN5pZSfuf3221m3bh179+61sn2bPd3zgDxjTIExphF4C5jfwXq/AB4H6r0RorCwkGHDhmkvV6kgERkZyb333ssf/vAHK9u3WXSTgQPt7hd7lp0kIucCw40xy87UkIgsFJFNIrLp8OHD3Q7Q1NRESUkJI0eO7H5qpVSfd9111/G3v/2Nqqoqn2/bb79IE5EQ4Engx12ta4x5wRiTZYzJOuecc7q9jQMHDpCYmEhEREQvkiql+prExEQuvfRSXn/9dZ9v22bRLQHaHy4wzLOsTTSQCawVkX3A+cBSp75MM8ZQWFioFyhXKkgtXLiQl156Cbfb7dPt2iy6G4F0ERklIuHA9cDStgeNMdXGmARjzEhjzEjgK+AKY4wj8yuXl5cTERHBwIEDnWhOKdXHnHvuuSQkJPDRRx/5dLvWiq4xphm4B1gJ5ALvGGNyRORREbnC29svKirSsVylgtzChQv54x//6NNtWh3TNcYsN8aMMcakGWMe8yx72BiztIN1v+tUL7euro6qqiqdUl2pIDd//nxCQ0PJy8vz2Tb99os0bzpw4ADJycl6yq9SQS48PJyMjAzeeecdn20z6KqOMYb9+/fryRBKKQCuv/563n33XVpaWnyyvaAruhUVFYSFhREbG2s7ilLKD0yYMIHY2Fg+//xzn2wv6Iqu9nKVUqe6/vrrefvtt32yraAqus3NzRw5coRhw4bZjqKU8iNXX301K1asoLa21uvbCvX6FvxIWVkZ0dHROv9ZgDPG0NTURHNzM83NzbjdbowxJ8fsjDEnv0R1uVyICC6X67QfFTzOOeccLrzwQj766COuvPJKr24rqIruwYMHGTp0qO0YymENDQ3U1dXR0NBAQ0MDbrcbl8tFeHg4oaGhuFwuQkNDCQkJISQkBBEBWosvQEtLC8YYGhsbMcbQ3NwMQGhoKGFhYYSGhp5sRwWuOXPm8Je//EWLrlPcbjeHDh1iypQptqOoXnK73Rw7dozjx49TW1tLaGgokZGR9O/fn7i4OMLDw08W1rPVVnzdbjdNTU00NDQAEBYWdrIQq8ByySWXsGjRIurq6hgwYIDXthM0fznl5eUn/0OqvscYQ3V1NVVVVZw4cYKoqCiio6NJSkrySgEUkZMFto3b7cbtdtPQ0EB9ff3J4qs94MAwcOBAzj33XNasWcPll1/ute0ETdHVoYW+qaGhgYqKCiorKxkwYADx8fGMGDHCyoktbWO94eHhJ8eNGxsbEZGTQxCqb7vssstYvny5V4tuUBy90NLSQnl5uZ7224ecOHGCwsJCdu3aBcCYMWNITU1l4MCBfnEmoYgQHh5O//79CQ8Pp6Wl5eR4suq75s6dy+rVq2lqavLaNuz/9frAkSNHiI6O1uvm9gH19fXk5+ezZ88eBgwYQGZmJkOHDvXrYaGQkBDCw8NP9oAbGxu1+PZRQ4YMITU1lS+++MJr2wiKz0OHDh2iJxc3V77X3NxMSUkJFRUVDB06lFGjRvlFj7Yn2oYZ2g5Pc7vdhISE0NLSwooVK9iyZQtTp05l3rx5Og7sx66++mq++OILvvOd73il/aAouuXl5XrUgh87cuQIRUVFDBo0iEmTJvX5sdG2436hdUqouXPnsmHDBmpra4mMjGT69OmsXLlSC6+fysrK4t577+XBBx/0Svt9qytxFtqO4YyPj7cdRZ2isbGRnTt3cvDgQTIyMhg5cmSfL7inWrlyJRs2bOD48eMYYzh+/Djr169nxYoVtqOpTkycOJGysjLKy8u90n7AF922oYXeHrepnFVRUcHXX39NdHQ0EydOJDIy0nYkr9iyZctpp5bW1taydetWS4lUV1wuFxdccAGfffaZV9oP+KJbXl5OYmKi7RjKwxhDQUEBhYWFZGRkMHz48IB+Q5w6deppbyghISHExMRYSqS6Y9asWaxbt84rbQd80a2oqNCi6yeamprYtm0b9fX1TJkyhejoaNuRvG7evHlMnz6dqKgoRISoqCjGjRvHk08+yd///nfb8VQn2opu26niTgqsAbRT1NXV0dTURFRUlO0oQa+uro6dO3eSkJAQVHPTuVwuVq5cyYoVK9i6dStTpkxh3rx5bNu2jX/5l3+hqKiIf//3f7cdU50iLS2N0NBQ9u3b5/iM4VZ7uiIyV0R2i0ieiCzq4PH7RWSniGwTkdUiMqIn7VdUVOgXaH6gpqaGrVu3MmzYsKAquG1cLheXX345Dz30EJdffjkul4upU6fy6aef8sknn/Af//EfXulRqbMnIkybNo1NmxyZlvEbrBVdEXEBzwLzgPHADSIy/pTVtgBZxphJwBLgNz3ZRmVlJYMGDXIirjpLVVVVbN++nYyMDJKSkmzH8SspKSn88Y9/ZNWqVfzqV7/Swutnpk6dypYtWxxv12ZP9zwgzxhTYIxpBN4C5rdfwRjzsTGmznP3K6BHVx+vqKjQomtRVVUVO3bsYPz48fqJoxMJCQksW7aM999/X3u8fmbKlCleOcrEZtFNBg60u1/sWdaZHwAdHtwoIgtFZJOIbDp8+DDQekWo6upq4uLinMqreqC6uprt27czYcIEfQ26EB8fz7Jly9i3bx9PP/207TjKIzMzk927d5+8rKdT+sTRCyJyE5AF/Lajx40xLxhjsowxWW2n+x49epTo6Gg968eC2tpatm3bpj3cHoiPj+fxxx/nueee44MPPrAdRwH9+/cnLS2NnTt3OtquzaJbAgxvd3+YZ9k3iMhs4GfAFcaYbr/lVFZWkpCQ0OuQqmcaGxvJzs4mNTVV938PJScn895773HPPfeQnZ1tO47CO+O6NovuRiBdREaJSDhwPbC0/QoiMhX4I60F91BPGq+urtZDxXzMGMO2bdtISkoiOflMI0WqM5mZmTzzzDM89thjHDlyxHacoJeVlUVJyWl9wV6xVnSNMc3APcBKIBd4xxiTIyKPisgVntV+C0QB74rIVhFZ2klzp6muriY2Ntbx3Kpzu3btQkQYPXq07Sh92j//8z+TmZnJnXfeeXIyTWVHWloaX375paNtWh3TNcYsN8aMMcakGWMe8yx72Biz1HN7tjFmsDFmiufnijO3+P/V1NRo0fWhsrIyKioqmDRpUkCf1usrDz30EDExMbz44ou2owS19PR08vPzHX3z67Loisi/iUif+vq5rq7u5LQqyvvq6+vZsWMHmZmZ35hTTJ29sLAwHnroIR5//HH27NljO07Qio6OJioqitLSUsfa7E5PdzCwUUTe8ZxB5vfdGO3l+taOHTsYMWIEAwcOtB0loKSmpvLzn/+cRx55RIcZLEpPTycvL8+x9rosusaYh4B04CXgVmCviPwfEUlzLIXDqqur9SpOPlJcXMyJEyd0HNdLbrzxRiorK/nzn/9sO0rQGj16NHv37nWsvW6N6ZrW02TKPD/NQBywRER6dFqurxw7dkyPXPCBxsZG9uzZo+O4XhQSEsITTzzBq6++ytGjR23HCUo+7+mKyL0ispnW6x58Dkw0xtwFTAOudiyJg9qmRVHetWvXLhITE3Uox8syMzOZMmUKTz31lO0oQWn06NHU19c71l53errxwFXGmDnGmHeNMU0AxpgWwHuTw/dCbW0tAwYMsB0joNXU1FBWVsbYsWNtRwkKDzzwAF999RUHDx60HSXoDBkyxNGTVbozpvuIMaaok8dyHUviIO3pet/OnTtJTU3VoxV8JDExkZkzZ/K73/3OdpSgk5SU5POjF/oUYwwul0uLgRdVVlZSXV3t+MWd1Zndfffd7N69m7KyMttRgsqAAQPo37+/Y+0FXNFtaWnRoQUvy83NZcyYMXoxIR8bNGgQEydO5OWXX7YdJegMGTLEsbYCrui63W49csGLampqqK6uJiUlxXaUoHTHHXewdetWTpw4YTtKUHHyAvwBV3SNMTq04EX5+fmkpqZqL9eSESNG0K9fP5YtW2Y7SlDRnu4ZGGPo16+f7RgBqampif379+tYrmU33XQTn332me0YQWXEiB5Nz3hGAVd0W1pa9JoLXlJUVMTgwYP1Tc2yiy66iDVr1lBcXGw7StBwsqYEXNE1xmjR9ZJ9+/YF5Wy+/iYsLIzLLruMlStX2o4SNJy8rEBAFl3tiTnv2LFj1NXVMXjwYNtRFHDZZZfx/vvv244RNKKjox1rK+CKrg4veEdJSQnDhw/Xayz4iRkzZpCXl6ezS/iI9nS7oEXXeQcOHGDo0KG2YyiPsLAwLr/8cr744gvbUYKC9nTPoKWlRQ9nclh9fT3Hjx+nbaZl5R/GjRunRddHtKfbhZCQgPy1rCkvLyc+Pj6g9uuxY8f4p3/6J44fP247ylk7//zzqaiosB0jKERERDjWltX/RZ6ZKHaLSJ6ILOrg8X4i8rbn8fUiMrI77QZScfAH5eXlAfcF2vLly9m0aRPLly+3HeWsjR49mvXr1+t1dn3AyU/P1qqTiLiAZ4F5wHjgBhEZf8pqPwCqjDGjgaeAx7tq1xijRddhVVVVJCYm2o7hiDvuuIOYmBhuvvlmoPVEg5iYGBYuXGg5Wc+5XC4mTJhATk6O7SgBLyCKLnAekGeMKTDGNAJvAfNPWWc+8Krn9hLg4u7M0aZF1znNzc1UVVUFzPxnP/3pTxk5ciShoaEAhIaGMnLkSB588EHLyc5ORkYGu3fvth0j4DlZU2xWp2TgQLv7xZ5lHa5jjGkGqoFBpzYkIgtFZJOIbKqpqdGi66CqqipiY2MDZp+OGjWKxYsX09TURFRUFE1NTSxevLjPntqcnp5OQUGB7RgBL1B6uo4xxrxgjMkyxmTFxMTosaQOOnr0KHFxcbZjOGrJkiVER0fz2GOPERUVxXvvvWc70llLS0tzdCoZ1TEnOx2hjrXUcyXA8Hb3h3mWdbROsYiEArHAGb+uNcbgdrv1sDGHVFdXB9ylMhctWsR//dd/kZCQwPe///0+fVHwpKQkvv76a9sxAl6gFN2NQLqIjKK1uF4PfP+UdZYCtwBfAtcAazwzE3dKRGhpafFC3OAUiMfnTpw48eTthIQEEhISLKbpncGDB1NfX09LS0vADAH5IydrirVXyTNGew+wEsgF3jHG5IjIoyJyhWe1l4BBIpIH3A+cdlhZR9xutzciB6WampqA6+kGkoiICOrq6qirq7MdJaA5edF4mz1djDHLgeWnLHu43e164NqetKk9XWedOHHC0fmhlPNiYmL0zdHLfD0Fe5+jPV1nGGNoaWlx9Gwc5bwJEyZoT9fLnOzpBlzR1Z6uc5qbm2lqatKxQj9XVlZGQ0OD7RgBTYvuGYgITU1NtmMEhObm5pMnESj/FRoaSnNzs+0YAU2HF85ARPRd3yHNzc06nttH6KcR79Ke7hlo0XWW7kv/19jYqLOleFl1dbVjbQVc0Q0JCdFC4RARoYvDopUf0MlCva+ystKxtgKu6IqInhbpkJCQEB3T7QNycnKIjY21HSOgOTktUkAWXe3pOiM8PJxjx47ZjqHOoKWlhWPHjjk6s4E6nZMXiw+4oqvDC84JDQ3F7XbrIXh+rKKigrFjx+oXaV6mRfcMQkJC+vQULP4mKirK0W9ulbPKysoICwuzHSPgadE9Ay26zgoPD6e2ttZ2DNWJ0tJS0tLSbMcIeE6eYh2QRbepqUlPkHBIVFSUjuv6sfz8/ICbv87fNDU1kZ2d7Vh7AVd0obVQaG/XGdHR0Y4eo6iclZeXx+jRo23HCGglJSUkJSU51l5AFt3o6Ggtug6Ji4tz9BhF5az6+nrS09Ntxwho+/fvZ/jw4V2v2E0BWXS1p+ucQYMGaU/XT1VWVpKbm0tKSortKAFt//79ju7jgCy62tN1Tnx8PIcPH9YLqvihr7/+mszMTD1czMsOHDigRbcrMTEx+pHYIaGhoScLr/IvW7ZsYfr06bZjBLyioiItul2Jj4/Xouug5OTkPj15Y6Bav34906ZNsx0j4IkII0eOdKy9gCy6UVFRNDY20tjYaDtKQEhKSqKoqMh2DNXOwYMHiYiIYMyYMbajBLQTJ06wYcMGRowY4VibVoquiMSLyEcistfzb1wH60wRkS9FJEdEtonI/+hB+/qtu4OGDx9OaWmpjuv6kbVr1+rpvz6wZ88eUlNTHb3wk61XbBGw2hiTDqym41l+64AFxpgJwFzgdyIysLsb0KLrnPDwcIYNG8aBAwdsR1EeH3/8MRdffLHtGAFv165dZGRkONqmraI7H3jVc/tV4MpTVzDG7DHG7PXcPggcAs7p7gbi4+OpqqpyIKqC1t7unj17bMdQtA4tuFwuHc/1gdzcXMaNG+dom7aK7mBjTKnndhlwxvMYReQ8IBzI7+TxhSKySUQ2tX3LrseXOis9PZ2CggIdYvADH3zwAePHj9drHftAn+rpisgqEdnRwc/89uuZ1qkJOp2eQESGAK8BtxljOrzGoDHmBWNMljEm65xzWjvDgwYNoqysTKdjd0hkZCTDhg0jP7/D9z3lI263mzVr1nD55ZfbjhLw3G43e/fudfzLSq+9VRpjZnf2mIiUi8gQY0ypp6ge6mS9GGAZ8DNjzFc92X5YWBixsbFUVFSQmJjYo+yqY2PGjGH79u2MHTvWdpSg9emnn3LOOeeQmppqO0rAy8vLIysry9ErjIG94YWlwC2e27cAfz91BREJB/4K/MkYs+RsNpKUlKTHlzpo9OjRHDlyRMfKLTHG8NZbb3HttdfajhIUsrOzvXIFN1tF99fAJSKyF5jtuY+IZInIi551rgO+DdwqIls9P1N6spHBgwdTXl7uZO6g5nK5mDx5Mps3b7YdJSh9/fXXREVFMXPmTNtRgkJ2drZXvqy0UnSNMRXGmIuNMenGmNnGmErP8k3GmB96bv/ZGBNmjJnS7mdrT7ajRdd5kydPJjc3l7q6OttRgooxhpdffplLLrlEj831AWMMmzdvDpyi6yvR0dEAehFuBw0YMIBx48Y5elFn1bXs7GwqKiq45JJLbEcJCsXFxYgIycnJjrcd0EUXYNiwYZSUlNiOEVDOO+88srOz9UpuPtLS0sIf/vAH7rrrLlwul+04QaGtlysijrcd8EU3KSlJz6RyWExMDJMmTeLzzz+3HSUoLF++nMjISC644ALbUYLGli1bvHbyScAX3eHDh1NSUkLr4cDKKeeffz4lJSV6yUcvq66uZtWqVfzrv/6rV3pd6nTGGL744gvOP/98r7Qf8EU3MjKSAQMGaHFwWHWiC3cAAA/3SURBVEREBOeeey4ffvihvqF50XPPPcfYsWP1amI+tGfPHiIiIhy9slh7AV90obW3q0MMzps8eTLGGLZu7dFBJaqbNmzYwI4dO7jttttsRwkqn376KbNmzfJa+0FTdPfv3287RsAREebOncu6dev0OhcOq6ys5IknnuDBBx8kIiLCdpygsm7dOi26vTVkyBCqq6upr6+3HSXgJCYmMn36dJYuXarDDA4xxvCb3/yGefPmkZmZaTtOUDl69Ch5eXlevYJbUBRdl8vF4MGDdfYDLznvvPMIDw/Xoxkc8vrrrxMdHc3NN99sO0rQ+eyzz5g+fTrh4eFe20ZQFF1ovW5AXl6e7RgBSUT43ve+x+bNm3Uf99Lnn3/OihUr9JhcS7Zv385FF13k1W0ETdEdOXIkBw8epKGhwXaUgBQVFcXVV1/N0qVLqaiosB2nT8rLy+Pdd9/lkUceIT4+3nacoFNXV8eyZcu8Op4LQVR0w8LCGDZsGPv27bMdJWANGzaM2bNn88Ybb1BbW2s7Tp9SWlrKI488wrXXXquHh1nyySefMGXKFAYO7PasYGclaIou6BCDL0yaNIlJkybxxhtv6GzM3VRRUcFvfvMbbrzxRj3rzKIPP/yQefPmeX07QVV0dYjBN77zne+QkpKihbcbjhw5wqJFi5g5cyaXXXaZ7ThBq7q6ms2bN/Pd737X69sKqqIbFhbGqFGjtLfrA5deeilxcXG89tpr+ibXifLych5++GHmzZvH1VdfbTtOUFu9ejUzZswgMjLS69sKqqILrUMMubm5tmMEPBHhiiuuYPDgwbz11ls6xnuKgoICfvKTnzB//nyuuuoq23GC3ooVK5g7d65PthV0RTclJYVjx45RWVlpO0rAazuULDk5mRdeeEGPavDYvHkzDz/8MHfeeSdz5syxHSfolZSUUFlZ6bMZOYKu6IaEhDBu3Dh27txpO0pQEBFmz57NrFmzeOutt4J6NmFjDG+//Ta///3vefjhh5kxY4btSAr4y1/+wgUXXEC/fv18sr2gK7oA48aNY/fu3To9uw9lZWUxb9483n33XdauXRt0pwzX1NTwzDPPsGXLFp588kk9LMxPNDU1sXTpUp8O8VgpuiISLyIfichez79xZ1g3RkSKReT3Tm0/NjaWQYMGUVhY6FSTqhtSU1O5++67yc/P59VXXw2ai+RkZ2fzox/9iIEDB/LLX/6SuLhO/9yVj61du5ZRo0YxcuRIn23TVk93EbDaGJMOrPbc78wvgE+dDjB+/Hh2797tdLOqCzExMdx6662kpKTwzDPPkJ2dHbC93mPHjvHcc8/x+uuvc99997FgwQJCQ0Ntx1LtvPfeez4/csRW0Z0PvOq5/SpwZUcricg0YDDwD6cDpKWlUVpaql/uWOByubjooou4/fbbyc7O5qWXXgqoi8wbY1i7di0/+tGP6NevH7/85S+ZOHGi7VjqFEVFReTn53PhhRf6dLtio5chIkeNMQM9twWoarvfbp0QYA1wEzAbyDLG3NNJewuBhQApKSnTuns1sQ0bNnDs2DEuvvjis/5dVO+43W6+/PJLVq9ezfTp0/n2t7/NgAEDbMc6a9u2bePPf/4zcXFxXHfddaSlpdmOpDrx/PPPExoayg9/+MNurS8im40xWb3drtc+64jIKiCpg4d+1v6OMcaISEeV/25guTGmuKu5oYwxLwAvAGRlZXX7XWTixIn86U9/YsaMGX36P3pf5nK5mDlzJpMnT2b16tU8/vjjfOtb32LWrFl95jUxxpCTk8M//vEPCgsLufHGG5k+fbrOaebHampqePfdd3nzzTd9vm1bPd3dwHeNMaUiMgRYa4wZe8o6rwOzgBYgCggHnjPGnGn8l6ysLLNp06ZuZ1mzZg2RkZFMnz69p7+G8oLKykpWrVpFbm4ukydPZtasWQwaNMh2rA41NTWxceNG3n//ferr65k/fz4zZ87Ucds+4KWXXqK4uJhHHnmk289xqqdrq+j+FqgwxvxaRBYB8caYn5xh/Vs5w/BCez0tulVVVSxZsoTbbrtN/7P4kerqaj777DM2btxIcnIyWVlZZGZmEhYWZjWXMYbi4mI++eQT1q1bx/Dhw5k7dy7Tpk3Tnm0f0dDQwBVXXMHzzz/PqFGjuv08vx9e6MKvgXdE5AdAEXAdgIhkAXcaY7o3yOKAuLg4Bg8eTG5urn7Z4UdiY2P53ve+x6WXXsq2bdtYv349S5YsISMjg8zMTDIyMujfv79PsrS0tFBUVER2djZfffUVDQ0NfPvb32bx4sUkJXU0gqb82fvvv8+ECRN6VHCdZKWn60097ekCHDx4kJUrV3LLLbcQEhKU54v0CdXV1Wzfvp2cnBwKCgoYMmQIqamppKamMmLECKKjox3ZTmNjI0VFRRQWFrJ792527txJbGwsWVlZTJkyhbS0NO3V9lEtLS1cddVVLF68mMmTJ/fouX29p+tXhg4dSkxMDLm5uUyYMMF2HNWJ2NhYZs6cycyZM2lsbKSwsJD8/HzWrl1LaWkpzc3NDBkyhKFDhxIZGcnAgQOJjY2lX79+hIeHnxw+crvdNDc3U1dXx7Fjxzh+/DhHjhyhvLyc8vJyKisrSUpKYtSoUWRlZXHTTTfpCQ0BYu3atSQlJfW44DpJi67HjBkz+PDDD8nIyNC5qfqA8PBwxo4dy9ixrd+/GmOoqamhtLSUI0eOUFlZSX5+Ps3NzVRUVNDQ0EBkZCRHjx7F5XIRHR1NS0sL0dHRxMbGEhcXR0ZGBomJiSQmJlofO1bOc7vdPP/889x///1Wc2jR9Rg6dChxcXHs3LlTx3b7IBEhNjaW2NhY21GUn1q+fDmDBg2yfqSSDmC2M2PGDDZs2KAXwlEqwDQ2NvLiiy9y1113WR+P16LbTlJSEgkJCWzfvt12FKWUg/76178yevRoJk2aZDuKFt1TzZgxg40bN+rcXkoFiBMnTvDKK69w55132o4CaNE9TWJiIikpKWzYsMF2FKWUA5YsWUJWVhbp6em2owBadDt0wQUXsG3btqC53qtSgerw4cO88cYbftPLBS26HYqKimLatGmsW7fOdhSlVC8888wzXHnllSQnJ9uOcpIW3U5MmzaNsrIyiouLbUdRSp2FrVu3sm3bNhYsWGA7yjcE3GnAInIM8JcpIRKAI7ZDeGiWjmmWjmmW0401xvT6XPNAPDlitxPnRztBRDZpltNplo5plo75SxYR6dlFXTqhwwtKKeVDWnSVUsqHArHovmA7QDuapWOapWOapWP+ksWRHAH3RZpSSvmzQOzpKqWU39Kiq5RSPtQni66IXCsiOSLS4plXrbP15orIbhHJ80yA2bZ8lIis9yx/W0TCe5ElXkQ+EpG9nn9Pm2JARC4Uka3tfupF5ErPY6+ISGG7x6Z4M4tnPXe77S1tt9zX+2WKiHzpeS23icj/aPdYr/dLZ69/u8f7eX7PPM/vPbLdYw96lu8WkTk93XYPc9wvIjs9+2C1iIxo91iHr5UXs9wqIofbbfOH7R67xfN67hWRW3yQ5al2OfaIyNF2jzm9X14WkUMisqOTx0VEnvZk3SYi57Z7rGf7xRjT536AccBYYC2tswR3tI4LyAdSaZ2+/WtgvOexd4DrPbefB+7qRZbfAIs8txcBj3exfjxQCQzw3H8FuMah/dKtLMDxTpb7dL8AY4B0z+2hQCkw0In9cqbXv906dwPPe25fD7ztuT3es34/YJSnHZcXc1zY7u/hrrYcZ3qtvJjlVuD3nfzdFnj+jfPcjvNmllPW/zfgZW/sF0973wbOBXZ08vhlwApAgPOB9We7X/pkT9cYk2uM6eqss/OAPGNMgTGmEXgLmC8iAlwELPGs9ypwZS/izPe00d22rgFWGGPqerFNp7KcZGO/GGP2GGP2em4fBA4B5/Rim+11+PqfIeMS4GLPfpgPvGWMaTDGFAJ5nva8ksMY83G7v4evgGFnua1eZzmDOcBHxphKY0wV8BEw14dZbgDe7MX2zsgY8ymtnaHOzAf+ZFp9BQwUkSGcxX7pk0W3m5KBA+3uF3uWDQKOGmOaT1l+tgYbY0o9t8uAwV2sfz2n//E85vnI8pSI9PNBlggR2SQiX7UNc2B5v4jIebT2ePLbLe7Nfuns9e9wHc/vXU3rfujOc53M0d4PaO1RtenotTpb3c1ytWe/LxGR4T18rtNZ8Ay3jALWtFvs5H7pjs7y9ni/+O1pwCKyCkjq4KGfGWP+7i9Z2t8xxhgR6fQYPM8740RgZbvFD9JalMJpPQ7wAeBRL2cZYYwpEZFUYI2IbKe14PSIw/vlNeAWY0yLZ3GP9ksgEJGbgCzgO+0Wn/ZaGWPyO27BEe8DbxpjGkTkf9L6SeAiL26vO64Hlhhj2s+j5ev94hi/LbrGmNm9bKIEGN7u/jDPsgpaPxqEeno3bcvPKouIlIvIEGNMqad4HDpDU9cBfzXGNLVru6032CAi/w38L29nMcaUeP4tEJG1wFTgPSzsFxGJAZbR+mb6Vbu2e7RfOtDZ69/ROsUiEgrE0vr30Z3nOpkDEZlN65vVd4wxDW3LO3mtzra4dJnFGFPR7u6LtI7Ntz33u6c8d+1Z5uhWlnauB/71lJxO7pfu6Cxvj/dLIA8vbATSpfUb+XBaX7ilpnX0+2Nax1YBbgF603Ne6mmjO22dNi7lKUhtY6pXAh1+e+pUFhGJa/uoLiIJwAXAThv7xfO6/JXWsbIlpzzW2/3S4et/hozXAGs8+2EpcL20Ht0wCkgHznYqkS5ziMhU4I/AFcaYQ+2Wd/hanWWO7mYZ0u7uFUCu5/ZK4FJPpjjgUr75ic3xLJ48GbR+QfVlu2VO75fuWAos8BzFcD5Q7ekY9Hy/OPkNoK9+gH+hdeykASgHVnqWDwWWt1vvMmAPre+AP2u3PJXW/0R5wLtAv15kGQSsBvYCq4B4z/Is4MV2642k9V0x5JTnrwG201pU/gxEeTML8C3P9r72/PsDW/sFuAloAra2+5ni1H7p6PWndYjiCs/tCM/vmef5vVPbPfdnnuftBub18u+1qxyrPH/HbftgaVevlRez/ArI8WzzYyCj3XNv9+yrPOA2b2fx3P858OtTnueN/fImrUfPNNFaW34A3Anc6XlcgGc9WbfT7qipnu4XPQ1YKaV8KJCHF5RSyu9o0VVKKR/SoquUUj6kRVcppXxIi65SSvmQFl2llPIhLbpKKeVDWnRV0BGRf/Jc0CVCRCKl9Xq+mbZzqeCgJ0eooCQiv6T1jLT+QLEx5leWI6kgoUVXBSXP+f4bgXrgW+abV7BSymt0eEEFq0FAFBBNa49XKZ/Qnq4KStI6r9ZbtF4ce4gx5h7LkVSQ8Nvr6SrlLSKyAGgyxrwhIi7gCxG5yBizpqvnKtVb2tNVSikf0jFdpZTyIS26SinlQ1p0lVLKh7ToKqWUD2nRVUopH9Kiq5RSPqRFVymlfOj/AYFB46UVGQCRAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 360x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"rebound.OrbitPlot(sim,xlim=[-1,1]);"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [],
"source": [
"def calculate_angle(sim):\n",
" v1 = np.array(sim.particles[1].xyz) - np.array(sim.particles[0].xyz)\n",
" v2 = np.array(sim.particles[2].xyz) - np.array(sim.particles[0].xyz)\n",
" v1 /= np.linalg.norm(v1)\n",
" v2 /= np.linalg.norm(v2)\n",
" return np.arccos(np.dot(v1,v2))/np.pi*180."
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.2867051844859263"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"calculate_angle(sim) # [deg]"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [],
"source": [
"Nsamples = 1000\n",
"times = np.linspace(0.,1.,Nsamples) # [days]\n",
"angles = np.zeros(Nsamples)\n",
"for i,t in enumerate(times):\n",
" sim.integrate(t*2.*np.pi/365.25)\n",
" angles[i] = calculate_angle(sim)"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZfbA8e9JCCSUkAhBepMeSoDQgvpDsbBWelFQQEHA7rquZVVQV921rYUiKE2pUhTL2kGE0EKvIr33DgmQ5Pz+mAsbkWQuJJPJZM7nee6Tuf1cSs7ce8/7vqKqGGOMCV4h/g7AGGOMf1kiMMaYIGeJwBhjgpwlAmOMCXKWCIwxJsgV8HcAl6pkyZJauXJlf4dhjDEBZfHixQdUNeZi6wIuEVSuXJmkpCR/h2GMMQFFRLZmts4eDRljTJCzRGCMMUHOEoExxgQ5SwTGGBPkLBEYY0yQs0RgjDFBzhKBMcYEuaBJBIdOnmHQl6tJOZvm71CMMSZPCZpEMHfDAUYnbqHbiPkcPHHa3+EYY0yeETSJ4PYGZRlyVyPW7DpG+6GJbD5w0t8hGWNMnhA0iQDgL/XKML5Pc44ln6X9kLks3nrI3yEZY4zfBVUiAGhcKZppA1pSPCKMbiMW8N+Vu/0dkjHG+FXQJQKAKiWLMLV/ArFlIxkwfgkf/boJG7vZGBOsgjIRAJQoWogJfZpzc53SvPL1WgZ9uYa0dEsGxpjgE7SJACA8LJTBdzfivqurMDpxC/0+XUzyGSsvNcYEl6BOBAChIcLzt9Xhxdvr8OPavXQdMZ8DVl5qjAkimSYCETnmZTouIutzM1hf6tWyCkPvbsy63cdoPySRjftP+DskY4zJFVndEWxU1cgspmJAvirGb1O3NBP6NufE6VQ6DE1k0RYrLzXG5H9ZJYIOLvZ3s01AaVQxmukDEoguXJC7P1rA1yusvNQYk79lmghUdZO3nbPaRkTCRWShiCwXkdUiMugi2/QUkf0issyZ7ncfuu9UKlGEaf0TqF+uOA+OX8Lw2RutvNQYk295fVnsvAu48P3AdhGZLiJVs9j1NHC9qjYA4oA2ItL8IttNUtU4Z/roMq8jx0UXKcin9zfj1nplePWbdbw4Y7WVlxpj8qUCLrb5D7ADGA8I0BW4ClgCjARaXWwn9XyFPvfGNcyZAuo3aXhYKO93a0i56AiGz97EriPJvNetIYULuvljM8aYwOCmfPQOVf1QVY+r6jFVHQ7crKqTgOisdhSRUBFZBuwDflDVBRfZrIOIrBCRKSJS4dIvwbdCQoRnb6nNS3fG8vO6fXQdPp/9x6281BiTf7hJBKdEpLOIhDhTZyDFWZflN3xVTVPVOKA80FRE6l6wyZdAZVWtD/wAjLnYcUSkr4gkiUjS/v37XYSc8+5pUZkPe8Szfu9x2g2Zy4Z9x/0ShzHG5DQ3ieBuoAeeb/V7nc/dRSQCeMjNSVT1CDATaHPB8oOqeu7r9UdA40z2H66q8aoaHxMT4+aUPnFjnSuZ1LcFKWfTaD8kkQWbDvotFmOMySleE4GqblLV21W1pKrGOJ83qGqyqs7JbD8RiRGRKOdzBHAjsO6CbcpkmL0DWHt5l5F7GlSIYvqAlpQsVogeHy/ki2U7/R2SMcZki5uqoRoi8pOIrHLm64vIP1wcuwwwU0RWAIvwvCP4SkReEpE7nG0ecUpLlwOPAD0v7zJyV4UrCjOtfwJxFaJ4dOIyhs6y8lJjTOASb7/AROQX4G/Ah6ra0Fm2SlUvfN6fK+Lj4zUpKckfp/6TlLNp/G3KCr5cvou7mlXkpTtiKRAa9N03GWPyIBFZrKrxF1vnpg6ysKouFJGMy1JzJLIAFx4Wyrtd4igXFcGwXzay+0gyH9zViCKFrLzUGBM43Hx9PSAiV+FUCIlIR8D6XXCEhAhP/6UWr7Styy/r99Nl+Dz2HUvxvqMxxuQRbhLBg8CHQC0R2Qk8BvT3aVQBqHvzSoy4J56N+07Sbkgiv++18lJjTGBwWzV0AxAD1FLVq1V1i88jC0Cta1/JpAeaczo1nfZDE5m30cpLjTF5X6Yvi0Xkiax2VNW3fRKRF3npZXFmth86Ra/Ri9h68CRvdGxA24bl/B2SMSbIZfWyOKs7gmLOFI/nUVA5Z+oHNMrpIPOTClcUZmq/BBpXiuaxScsYPHODlZcaY/KsTMtbVHUQgIjMBhqp6nFnfiDwda5EF8CKFw5jTO+mPDVlBW989xs7Dp/i5TvrWnmpMSbPcVPneCVwJsP8GWeZ8aJQgVD+0yWOCtGF+WDmBnYdSWHw3Y0oauWlxpg8xM3X07HAQhEZ6NwNLABG+zKo/EREePLmmrzWvh5zNhyg87B57LXyUmNMHuKmauifQC/gsDP1UtXXfB1YftOtaUU+ujeerQdP0m7wXH7bY+Wlxpi8wdUDa1VdoqrvOtNSXweVX11XsxSTHmhBarrScWgiiRsO+DskY4zJPBGIyBJvO7vZxvxR3XLFmf5gS8pEhXPvqIVMW7LD3yEZY4JcVm8tazs9h2ZGgOI5HE9QKBcVwWf9Euj3yWKemLycHYeTefj6alzQn5MxxuSKrBJBLRf7p+VUIMGmeISnvPTpqSt4+4f17Dh8in+2q0eYlZcaY3JZVu0ItuZmIMGoYIEQ3urcgPLREbz38wZ2H01hyN2NKBYe5u/QjDFBxL5++pmI8MRNNfl3h/okbjxIp2Hz2H002d9hGWOCiCWCPKJzkwqM7NmE7YdO0W5wImt3H/N3SMaYIGGJIA/5vxoxTO7XAkXpNGwev/6+398hGWOCgJsxi9uLyO8iclREjonIcRHx+nVVRMJFZKGILHfGJR50kW0KicgkEdkgIgtEpPLlXUb+EVu2ONMHtKR8dAS9Ri1ictJ2f4dkjMnn3NwR/Bu4Q1WLq2qkqhZT1UgX+50GrlfVBkAc0EZEml+wzX3AYVWtBrwD/OtSgs+vykZFMLlfC5pXLcFTU1bwzg/rrfdSY4zPuEkEe1V17aUeWD1OOLNhznThb7M7gTHO5ylAa7FiegAiw8MY1asJHRuX592ffufJz1ZwJjXd32EZY/IhN91gJonIJOBzPN/yAVDVad52FJFQYDFQDRisqgsu2KQcsN05XqqIHAVKAAcuOE5foC9AxYoVXYScP4SFhvBGx/pUiC7MOz+uZ8+xZIZ2b0yklZcaY3KQmzuCSOAUcBNwuzPd5ubgqpqmqnFAeaCpiNS9nCBVdbiqxqtqfExMzOUcImCJCI/eUJ03OzVgwaZDdBo6j11HrLzUGJNzvN4RqGqv7J5EVY+IyEygDbAqw6qdQAVgh4gUwNNlhQ30exEdG5endGQ4/T9dTLshcxnZswmxZa2HD2NM9rmpGholIiMvnFzsFyMiUc7nCOBGYN0Fm80A7nU+dwR+Vnsrmqmrq5fks/4tCBGh87B5/LLeykuNMdnn5tHQV3iGpvwa+AnPo6ITWe7hUQaY6XRctwj4QVW/EpGXROQOZ5uPgRIisgF4Anj6Ui8g2NQqHcn0AS2pWKIIvUcvYtKibf4OyRgT4ORSv4CLSAgwR1UTfBNS1uLj4zUpKckfp85Tjqec5cHxS5m9fj8PX1+NJ26sYb2XGmMyJSKLVTX+Yusup2VxdaBU9kIy2VUsPIyP742nS3wF3v95A3+dvNzKS40xl8Xry2IROY6n/l+cn3uAv/s4LuNCWGgIr3eoR/noCN76YT27j6YwrEdjikdYeakxxj03YxYXy9iiWFVrqOrU3AjOeCciPNy6Om93bkDS1kN0GpbITisvNcZcgqyGqqwkIsUzzF8nIu+KyOMiUjB3wjNutW9UnjG9mrL7aArtBs9l1c6j/g7JGBMgsrojmAwUARCROOAzYBuefoOG+D40c6kSqpVkav8ECoQInT+cx8zf9vk7JGNMAMgqEUSo6i7nc3dgpKq+BfQCmvo8MnNZalxZjOkPtqRKySLcPyaJ8QusvNQYk7WsEkHGWsTr8bQhQFWtNCWPuzIynEkPtOCa6iV5dvpK3vhunfVeaozJVFaJYKaITBaRd4Fo4GcAESkDnMmN4MzlK1qoAB/dE0+3phUZPHMjj01axunUNH+HZYzJg7IavP4REemKp4Xw1ap61llVGnguN4Iz2VMgNIRX29WlwhUR/Pvb39hzNIXhPeIpXtjKS40x/5NV1dB3eJLAf1V157nlqrpUVb/LjeBM9okIA1pV492ucSzddoQOwxLZfuiUv8MyxuQhWT0auhc4DAwUkSUiMlRE7hSRIrkUm8lBd8aVY+x9Tdl3LIV2QxJZucPKS40xHpkmAlXdo6qjVbUrEA+MBRoD34vIjyLyVG4FaXJG86olmDYggUIFQjzlpeusvNQY464b6hKqmq6q81T1BVVtCXTFM5aACTDVShVj+oAEripVhPvHWnmpMcZdp3PzReQzEbnl3HjCqnpAVcf5ODbjI6Uiw5nU18pLjTEebhJBDWA40AP4XUReFZEavg3L+FoRp7y0a5MKDJ65kccnLbPeS40JUm46nVNV/UFVuwF98LxEXigiv4hIC59HaHymQGgIr7Wvx99ursnny3Zx78iFHE0+631HY0y+4uodgYg8KiJJwJPAw0BJ4K/AeB/HZ3xMRHjwumq83bkBi7Z4ei/dZb2XGhNU3DwamodneMq2qnqrqk5T1VRVTQKG+TY8k1vaNyrPmN5N2X0khXZD5rJ6l5WXGhMsskwEIhIKfKmqL6vqjgvXq+q/sti3gojMFJE1IrJaRB69yDatROSoiCxzphcu6ypMjmhZrSSf9W9BiAhdPpzP7PX7/R2SMSYXZJkIVDUNuNyxiVOBv6pqHaA58KCI1LnIdr+qapwzvXSZ5zI5pFbpSKYPaEn56Ah6j17EZ0nb/R2SMcbH3DwaWiYiM0Skh4i0Pzd520lVd6vqEufzcWAtUC6b8ZpcULp4OJ/1a0GLq0rwtykr+M+P66281Jh8zE0iCAcO4umK+nZnuu1STiIilYGGwIKLrG4hIstF5L8iEpvJ/n1FJElEkvbvt8cVuaFYeBgjezahY+Py/OfH33lqygrOpll5qTH5kfj6m56IFAV+Af6pqtMuWBcJpKvqCRG5BXhXVatndbz4+HhNSkryXcDmD1SV//z4O+/+9DvXVC/JkLsbUSzcei81JtCIyGJVjb/Yuky7oc6w8yjgT9lCVXu72DcMmAqMuzAJOMc4luHzNyIyRERKquoBb8c2uUNEePzGGpSLiuCZ6Svp/OF8RvdqwpWR4f4OzRiTQ9w8GvoK+NqZfsJTSnrC205OdxQfA2tV9e1Mtil9rtsKEWnqxHPQXegmN3VuUoGRPZuw7eBJ2g2ey297jvs7JGNMDrnkR0MiEgLMUdUsq4lE5GrgV2AlcO7h8rNARQBVHSYiDwH98VQYJQNPqGpiVse1R0P+tWrnUXqPXkTy2TQ+7NGYhKtK+jskY4wLWT0aupxEUBP4WlWr5URwl8oSgf/tPJJMr1EL2XzgJG90bEDbhlYMZkxel1UicNPFxHEROXZuAr4E/p7TQZrAUS4qgs/6JdC4UjSPTVrG4JkbrLzUmADm9WWxqhbLjUBMYCkeEcaY3k15asoK3vjuN3YeSealO2IpEOrmtZMxJi9xc0fQTkSKZ5iPEpG2vg3LBIJCBUJ5p3McA1pdxfgF2+gzNomTp1P9HZYx5hK5+fr2oqqe74FMVY8AL/ouJBNIQkKEp9rU4p/t6vLL+v10HT6ffcdT/B2WMeYSuEkEF9vG6yMlE1zublaJEffEs2HfCdoPSWTDPq8VxsaYPMJNIkgSkbdF5CpnehtY7OvATOBpXftKJj3QnJSzaXQYmsjCzYf8HZIxxgU3ieBh4AwwCZgIpAAP+jIoE7jql49i+oCWlChakO4fLeCrFbv8HZIxxossH/E44xEMUtUncykekw9UuKIwU/sl0PeTJB4av5RdR5Lpc01VnEbkxpg8xs14BFfnUiwmH4kuUpBP7mvGrfXK8Oo36xg4YzVp6dbWwJi8yM1L36UiMgP4DDh5buHFOpEzJqPwsFDe79aQslHhjPh1M7uOpvBe14ZEFAz1d2jGmAxyZTwCE7xCQoTnbq3DwNvr8OPavXQbMZ+DJ077OyxjTAZuWhb3yo1ATP7Ws2UVSheP4NGJS2k/NJHRvZpSpWQRf4dljMFdy+IaIvKTiKxy5uuLyD98H5rJb9rULc2Evs05npJKh6GJLN562N8hGWNw92hoBPAMcBZAVVcAXX0ZlMm/GlWMZlr/BIqFF+CuEfP5dtUef4dkTNBzkwgKq+rCC5ZZhzLmslUuWYRp/ROoXSaS/uMWM2ruZn+HZExQc5MIDojIVTjDVYpIR2C3T6My+V6JooWY0Kc5N9a+kkFfruGVr9aQbuWlxviFm0TwIPAhUEtEdgKP4RlVzJhsiSgYytDujemZUJmP5mzmoQlLSDmb5u+wjAk6XhOBqm5S1RuAGKCWql6tqlu87SciFURkpoisEZHVIvLoRbYREXlPRDaIyAoRaXRZV2ECVmiI8OLtdXjultp8s3IP3T9awOGTZ/wdljFBJdPyURF5IpPlAGQ2IH0GqcBfVXWJiBQDFovID6q6JsM2fwGqO1MzYKjz0wQREaHPtVUpGxXB45OX0WFYIqN7NqViicL+Ds2YoJDVHcGbQHegBFAUKHbBlCVV3a2qS5zPx4G1wIWD294JjFWP+UCUiJS55Ksw+cKt9csw7v5mHDxxhvZD57J8+xF/h2RMUMgqETQEvgduBSoBc4GXVHWQqg66lJOISGXneAsuWFUO2J5hfgd/ThaISF8RSRKRpP3791/KqU2AaVL5Cqb2TyA8LJSuw+fz09q9/g7JmHwv00SgqstV9WlVjQM+xvPtfY2I3HEpJxCRosBU4DFVPXY5QarqcFWNV9X4mJiYyzmECSDVShVl+oCWVL+yKH3GJvHp/K3+DsmYfM1Ny+IYPN/m6+H5xr7P7cFFJAxPEhiXSSd1O4EKGebLO8tMkIspVoiJfZvTqmYp/vH5Kv717TorLzXGRzJNBCLSW0S+xdPrqACdVfVG51m+V+J5q/wxsDaLF8szgHuc6qHmwFFVtTYKBoDCBQswvEdj7mpWkaGzNvL45GWcTrXyUmNyWladzn0ErAK2AjcDN2UcWERVvT0iagn0AFaKyDJn2bNARWf/YcA3wC3ABuAUYB3cmT8oEBrCP9vWpXx0BP/+9jf2HE1heI94ihcO83doxuQbWSWC67JzYFWdg+dOIqttFBv20nghIgxoVY2yxSP425TldByWyKheTSgfbeWlxuSETBOBqv6Sm4EY403bhuUoFVmIBz5ZTPshiYzs2YS65Yr7OyxjAp6bLiaMyTMSrirJ1P4JFAgRunw4j1/WWzmxMdllicAEnBpXFmP6gy2pVKIIvUcvYtKibf4OyZiA5joRiIg9kDV5xpWR4Uzu14KW1Ury96krefv73/C8cjLGXCo37QgSRGQNsM6ZbyAiQ3wemTFeFC1UgI/vjadzfHne+3kDT362gjOp6f4Oy5iA43XMYuAdPOWjM8DT4lhErvVpVMa4FBYawr861KdcVGHe+XE9e4+lMLR7I4qFW3mpMW65ejSkqtsvWGStekyeISI8ekN13uhYn/mbDtJp2Dx2H032d1jGBAw3iWC7iCQAKiJhIvIknp5EjclTOsVXYGTPJuw4nEy7wYms23NZXVsZE3TcJIJ+eBp9lcPTD1Ac1gjM5FHX1ohh8gMtUJROQ+cxd8MBf4dkTJ7nZoSyA6p6t6peqaqlVLW7qh7MjeCMuRx1ykYyfUBLykZFcO/IhUxbssPfIRmTp7mpGhojIlEZ5qNFZKRvwzIme8pGRTC5XwuaVrmCJyYv54Off7fyUmMy4ebRUH1VPT9UlKoextMttTF5WvGIMEb3akr7huV48/v1PDNtJWfTrLzUmAu5KR8NEZFoJwEgIle43M8YvytYIIS3OjegbFQEH8zcwO6jKQy5uxFFCtk/YWPOcXNH8BYwT0ReFpFXgETg374Ny5icIyI8eXNNXm1XjzkbDtBl+Dz2HUvxd1jG5BluXhaPBToAe4E9QHtV/cTXgRmT0+5qVpGP7oln0/6TtBuSyO97j/s7JGPyBLd9Da0DpuFpXXxCRCr6LiRjfOe6WqWY1LcFp1PT6TA0kfmbrADOGDdVQw/juRv4AfgK+Nr5aUxAqle+ONMHJBBTrBD3fLyQGct3+TskY/zKzR3Bo0BNVY1V1fqqWk9V6/s6MGN8qcIVhZnWvyVxFaN4ZMJShv2y0cpLTdBy1cUEcPRSDywiI0Vkn4isymR9KxE5KiLLnOmFSz2HMdlRvHAYY3s35bb6ZXj9v+t44YvVpKVbMjDBx00N3SZgloh8DZw+t1BV3/ay32jgA2BsFtv8qqq3uYjBGJ8IDwvlva4NKRcVwYezN7H7aArvd2tIRMFQf4dmTK5xc0ewDc/7gYJAsQxTllR1NnAoW9EZkwtCQoRnbqnNS3fG8vO6vXQdMZ8DJ05739GYfMLrHYGqDvLh+VuIyHJgF/Ckqq6+2EYi0hfoC1CxohUsGd+4p0VlSkeG88jEpbQfksjoXk2oGlPU32EZ43NuqoZiROQNEflGRH4+N+XAuZcAlVS1AfA+8HlmG6rqcFWNV9X4mJiYHDi1MRd3U2xpJvRpzsnTqXQYmsjirXZTa/I/N4+GxuFpR1AFGARsARZl98SqekxVTzifvwHCRKRkdo9rTHY1rBjNtAEJFI8I464RC/h21W5/h2SMT7lJBCVU9WPgrKr+oqq9geuze2IRKS0i4nxu6sRirXtMnlCpRBGmDWhJbNlI+o9bwsg5m/0dkjE+46Zq6Kzzc7eI3Irnef4V3nYSkQlAK6CkiOwAXgTCAFR1GNAR6C8iqUAy0FWtkNvkIVcUKcj4Ps15dOJSXvpqDTuPJPPcLbUJCRF/h2ZMjhJvv3tF5DbgV6ACnmf5kcBAVf3S9+H9WXx8vCYlJfnj1CZIpaUrL3+1htGJW7ilXmne7hxHeJiVl5rAIiKLVTX+Yuvc3BEcVtWjeBqVXeccsGUOxmdMnhYaIgy8I5by0RG88vVa9h1bwIh74okuUtDfoRmTI9y8I3jf5TJj8rX7r6nK4LsasWLnUToMTWTbwVP+DsmYHJHpHYGItAASgBgReSLDqkjA7otNULq1fhlKRRaiz9gk2g+dy8f3NqFBhSjvOxqTh2V1R1AQKIonWWRsUXwMz4teY4JSk8pXMLV/AuFhoXQdPp+f1u71d0gmCPiylsbNy+JKqrrV+RwCFFXVYz6LyAt7WWzyiv3HT3PfmEWs2nmUl+6sS/fmlfwdksmnTp1J5dGJy+jQqDxt6pa+rGNk9bLYzTuC10QkUkSKAKuANSLyt8uKxJh8JKZYISb2bU6rmqX4x+er+Ne360i33ktNDjt44jTdRizgp7V7OXTyjE/O4SYR1HHuANoC/8XTwriHT6IxJsAULliA4T0ac1ezigydtZHHJy/jdGqav8My+cTWgyfpMDSRdbuPMay759+ZL7gpHw0TkTA8ieADVT0rIva1xxhHgdAQ/tm2LuWjI/j3t7+x91gKH/aIp3hEmL9DMwFs+fYj9B69iHRVxvdpTuNK0T47l5s7gg/x9C9UBJgtIpXwvDA2xjhEhAGtqvGfLnEs3nqYjkMT2Xkk2d9hmQA1c90+ug6fT0TBUKb0T/BpEgAXiUBV31PVcqp6i3psxWlYZoz5o7YNyzGmd1P2HEuh3eC5rN51yYP7mSA3edF27h+bRNWYIkwbkMBVudAVeqaJQES6Oz+fuHACHvF5ZMYEqISrSjKlXwIFQoTOw+bxy/r9/g7JBABV5d0ff+epqStIuKoEkx5oQali4bly7qzuCIo4P4tlMhljMlGzdDGmP9iSiiWK0Hv0IiYv2u7vkEwelpqWzrPTV/LOj+tp36gcI3s2oWghN69wc4bXdgR5jbUjMIHkeMpZBoxbwq+/H+DR1tV57IbqOL2vGwN42gg8PH4pP63bx0PXVeOvN9Xwyb+Ry+p0TkTey+qgqmqPh4zxolh4GCN7NuHZaSt596ff2XUkmVfb1yMs1E2dhsnvDp44Te8xSazccYRX2vqvUWJW9x6LnZ8tgTrAJGe+E7DGl0EZk5+EhYbw7471KRcdwX9+/J09x1IYcncjioVbeWkw23rwJPeOXMjuoykM696Ym2Ivr8VwTnDTxcR84GpVTXXmw4BfVbV5LsT3J/ZoyASyyUnbeXbaSqpfWYxRPZtQunjuvAw0eUvGNgIf3dvE5+WhkP0uJqLx9Dh6TlFnmTHmEnWOr8DInk3YdvAk7YfM5bc9x/0dksllud1GwA03ieB1YKmIjBaRMcAS4FVvO4nISBHZJyKrMlkvIvKeiGwQkRUi0ujSQjcmMF1bI4bJ/VqQmq50HJZI4sYD/g7J5BJ/tBFww02DslFAM2A6MA1ooapjXBx7NNAmi/V/Aao7U19gqItjGpMvxJYtzvQHW1KmeDj3jlzI50t3+jsk40P+bCPghqvSBVXdo6pfONMel/vMBg5lscmdwFintfJ8IEpEyrg5tjH5QbmoCD7r53k08NikZQyeucGnfc4b//B3GwE3/FnDVg7I2Mpmh7PsT0Skr4gkiUjS/v3WStPkH8UjwhjTuylt48ryxne/8dznq0hNS/d3WCaHnDqTygOfLGbCwu08dF013urUIE+WDuettJQJVR0ODAdP1ZCfwzEmRxUqEMrbneMoGxXBkFkb2XM0hfe7NaRIHvvWaC5NXmkj4IbX1CQin7hZdhl2AhUyzJd3lhkTdEJChKfa1OKVtnWZ9ZunqmT/8dP+DstcpgvHEcjLSQDcPRqKzTgjIqFA4xw49wzgHqd6qDlwVFV358BxjQlY3ZtXYsQ98WzYd4J2Q+ayYd8Jf4dkLtHy7UdoPySRo8lnGd+nuV8birmVVe+jz4jIcaC+iBxzpuPAPuALbwcWkQnAPKCmiOwQkftEpJ+I9HM2+QbYBGwARgADsnsxxuQHrWtfycS+zUk5m0aHoYks2pJVzYXJS/JiGwE33LQsfk1Vn8mleLyylsUmWGw7eIqeoxay40gy73SO49b6VlSXl01etJ1npq+kVulijOrVJE+Vh0I2WycJJ8QAABruSURBVBar6jMiUk5EEkTk2nNTzodpjMmoYonCTO2fQP1yxXlw/BI++nWTlZfmQXm9jYAbXssSROR1oCuejubOjcqtwGwfxmWMAaKLFOTT+5vxxORlvPL1WnYcTub52+oQGmJdWecFqWnpPP/FKiYs3E77RuX4V4f6ebI81Bs39WntgJqqaiUMxvhBeFgoH3RrxKvF1/LRnM3sPprMu10bEh4W6u/QglpujSOQG9ykrk2A9ZdrjB+FhAj/uK0OL95eh+/X7KXbiPkcOnnG32EFrYMnTtNtxAJm/raPV9rW5cmbawZsEgB3dwSngGUi8hNw/q7ABqYxJvf1almFMsXDeXTiMtoPmcvoXk2pXLKI9x1NjslL4wjkFDd3BDOAl4FEPIPVnJuMMX7Qpm4ZxvdpztHks7QfmsjSbYf9HVLQCMQ2Am64GrNYRCKAiqr6m+9DypqVjxrjsWn/CXqOWsS+4ym817VhvvmllFfNXLePAeOWUKJoQcb0bppnupB2K1vloyJyO7AM+NaZjxORGTkbojHmUlWNKcq0AQnULB3JA58uZuy8Lf4OKd/Kq+MI5BQ3j4YGAk2BIwCqugyo6sOYfO7zzz9HRFi3bp2/Q8nUq696HfvnD1544QV+/PFHn8Ry+vRpbrjhBuLi4pg0adIf1s2fP59mzZoRFxdH7dq1GThwoE9iMBdXsmghJvZpTutaV/LCF6t57Zu1pKdbW4Ockh/aCLiiqllOwHzn59IMy1Z4289XU+PGjTW7OnfurFdffbW+8MIL2T6WN2fPnr2s/YoUKZLDkVy+efPmaevWrS+6rkaNGrps2TJVVU1NTdXVq1fnZmjGkZqWrs9/vlIr/f0rfXDcYk0+k+rvkALe2dQ0fXrqcq3096/08UlL9Uxqmr9DyhYgSTP5vermjmC1iNwFhIpIdRF5H8+L44B04sQJ5syZw8cff8zEiRPPL581axbXXnstt956KzVr1qRfv36kp3v6hS9atCiPP/44sbGxtG7dmnNjIowYMYImTZrQoEEDOnTowKlTpwDo2bMn/fr1o1mzZjz11FOcPHmS3r1707RpUxo2bMgXX3i6aho9ejTt27enTZs2VK9enaeeegqAp59+muTkZOLi4rj77rv/EH9aWho9e/akbt261KtXj3feeef8OadMmUJSUhJxcXHExcVRr1698yVtGzdupE2bNjRu3JhrrrnmondDhw4dom3bttSvX5/mzZuzYsUK9u3bR/fu3Vm0aBFxcXFs3LjxD/vs27ePMmU8XR+EhoZSp04dAAYOHMibb755fru6deuyZcsWtmzZQu3atenTpw+xsbHcdNNNJCcnX85fpckgNEQYdEcsz/ylFl+t2M09Ixdy5JSVl16uk6dT6TM2iQkLt/PgdVfl2XEEckxmGeLcBBQG/gkscqZXgHBv+/lqyu4dwaeffqq9e/dWVdUWLVpoUlKSqqrOnDlTCxUqpBs3btTU1FS94YYb9LPPPjuXSfXTTz9VVdVBgwbpgw8+qKqqBw4cOH/c5557Tt977z1VVb333nv11ltv1dRUz7eyZ555Rj/55BNVVT18+LBWr15dT5w4oaNGjdIqVarokSNHNDk5WStWrKjbtm1T1czvCJKSkvSGG244P3/48OHz5zwX7zlPPvmkPvnkk6qqev311+v69etVVXX+/Pl63XXX/enYDz30kA4cOFBVVX/66Sdt0KDB+T+bW2+99aLxDBo0SKOiorRt27Y6bNgwTU5OVlXVF198Ud94443z28XGxurmzZt18+bNGhoaqkuXLlVV1U6dOp3/szE544tlO7X6s9/o9W/O1G0HT/o7nICz71iK3v7+r1rl6a/00/lb/B1OjuFy7wicLqdfUtXnVLWJM/1DVVN8m558Z8KECXTt2hWArl27MmHChPPrmjZtStWqVQkNDaVbt27MmTMHgJCQELp06QJA9+7dzy9ftWoV11xzDfXq1WPcuHGsXr36/LE6depEaKin5ef333/P66+/TlxcHK1atSIlJYVt27YB0Lp1a4oXL054eDh16tRh69atWcZftWpVNm3axMMPP8y3335LZGTkRbebNGkSS5Ys4fXXX+fEiRMkJibSqVMn4uLieOCBB9i9+889fs+ZM4cePXoAcP3113Pw4EGOHTuWZTwvvPACSUlJ3HTTTYwfP542bbIaptqjSpUqxMXFAdC4cWO2bNnidR/j3h0NyvLJfU3Zf/w07YcmsmrnUX+HFDA27T9Bh6GJrN97nOE94rm7Wd4eRyCnZNmgTFXTROTq3ArG1w4dOsTPP//MypUrERHS0tIQEd544w2AP7UMzKyl4LnlPXv25PPPP6dBgwaMHj2aWbNmnd+mSJH/NfJRVaZOnUrNmjX/cJwFCxZQqFCh8/OhoaGkpqZmeQ3R0dEsX76c7777jmHDhjF58mRGjhz5h21WrVrFwIEDmT17NqGhoaSnpxMVFcWyZcuyPPbluuqqq+jfvz99+vQhJiaGgwcPUqBAgfOP1gBSUv733eHCa7ZHQzmvWdUSTBuQwL0jF9H5w3kMvrsR19Us5e+w8rTFWw9z/5hFhIgwsW8L4ipE+TukXOPmoddSEZkhIj1EpP25yeeR+cCUKVPo0aMHW7duZcuWLWzfvp0qVarw66+/ArBw4UI2b95Meno6kyZN4uqrPTkwPT2dKVOmADB+/Pjzy48fP06ZMmU4e/Ys48aNy/S8N998M++///75niOXLl3qNdawsDDOnj37p+UHDhwgPT2dDh068Morr7BkyZI/rD9y5AjdunVj7NixxMTEABAZGUmVKlX47LPPAE9iWr58+Z+Ofc0115y/jlmzZlGyZMlM7zjO+frrr89f1++//05oaChRUVFUrlz5fGxLlixh8+bNXq/Z5KxqpYoxfUACVWOKcP+YJCYs3ObvkPKs71bv4a4R8ykeEcbU/glBlQTAXSIIBw4C1wO3O9NtvgzKVyZMmEC7du3+sKxDhw7nHw81adKEhx56iNq1a1OlSpXz2xYpUoSFCxdSt25dfv75Z1544QUAXn75ZZo1a0bLli2pVatWpud9/vnnOXv2LPXr1yc2Npbnn3/ea6x9+/alfv36f3pZvHPnTlq1akVcXBzdu3fntdde+8P6L774gq1bt9KnT5/zL40Bxo0bx8cff0yDBg2IjY09/8I6o4EDB7J48WLq16/P008/zZgxY7zG+cknn1CzZk3i4uLo0aMH48aNIzQ0lA4dOnDo0CFiY2P54IMPqFGjhtdjmZxXKjKcSX1bcE31kjwzbSVvff+bdWV9gbHzttD/08XULhPJ1P4JQdllh6uWxXmJr1oWz5o1izfffJOvvvrqT+uKFi3KiRM2ZKAJXKlp6fzj81VMXLSd9g3L8XqH+hQskI+rYFxIT1f+/d1vDPtlIzfULsX73RoRUTD/9uiaVctiN+MRhAP34Rm7+HxLClXtnWMRGmN8qkBoCK+1r0e5qAje+mE9e4+nMLR7YyLDg7Nj4dOpaTw1ZQVfLNvF3c0qMuiOWArk5/JQL9xc+SdAaeBm4BegPHDczcFFpI2I/CYiG0Tk6Yus7yki+0VkmTPdfynB56RWrVpd9G4AsLsBky+ICA+3rs5bnRqwYNMhOg+bx+6jwfei/mjyWXqOXMQXy3bxVJuavNK2blAnAXCXCKqp6vPASVUdA9wKNPO2k1N6Ohj4C1AH6CYidS6y6SRVjXOmjy4hdmPMZejQuDyjezVlx+Fk2g1OZO3urEuE85PdR5PpPGweSVsP8U6XBgxoVS2gxxHIKW7GIzhXunJEROoCewA3dWhNgQ2quglARCYCd+IZ8tIvWrVq5a9TG5PnRBaOYW/NDtzy1o+UWv8FEceybsMS6M5ElGRvrY6khxak1PrPeXfuv3nX30Fdoowl6jnJzR3BcBGJBp7HMzbBGuBfLvYrB2zPML/DWXahDiKyQkSmiEiFix1IRPqKSJKIJJ3r3sEYkz0FT+2nzOpPKXD6GHtrdeBEyVh/h+QzyZEV2R17FwiUWTOBiGNWSpuRz6qGRKQj0EZV73fmewDNVPWhDNuUAE6o6mkReQDooqrXZ3VcG4/AmJx1LOUs/T9dzNwNB/nrjTV46Pr89bjki2U7efKz5VQuUYTRvZtSLirC3yH5RXbHIyghIu+LyBIRWSwi/3F+gXuzE8j4Db+8s+w8VT2oqueGv/wIaOziuMaYHBQZHsaonk1p36gcb/2wnmemreRsWrr3HfM4VWXYLxt5dOIyGlWMZkq/hKBNAt64eTQ0EdgHdAA6AgeASVnu4bEIqC4iVUSkINAVz6Ol80SkTIbZO4C1boI2xuSsggVCeKtTAx65vhoTF23n/jFJnDiddXcneVlauvLijNW8/t913N6gLGPva0rxwsFZKuuGm0RQRlVfVtXNzvQKcKW3nVQ1FXgI+A7PL/jJqrpaRF4SkTuczR4RkdUishx4BOh5eZdhjMkuEeGJm2ryevt6zNlwgC4fzmPfscDrXzLlbBr9P13M2Hlb6XttVd7tEkehAvm3oVhO8PqOQETeBhYCk51FHYGmqvqkj2O7KHtHYIzvzfxtHw+OW0J04YKM7tWE6lcW83dIrhw6eYb7xyxi6fYjvHBbHXq1rOLvkPKMrN4RuEkEx4EiwLmHhiHASeezqmrWvZLlMEsExuSOVTuP0mv0Ik6fTWP4PfE0r+rm1aD/bDt4intHLWTnkWTe7RLHX+qV8b5TEMnWy2JVLaaqIapawJlCnGXFcjsJGGNyT91yxZnWP4FSkeHc8/FCZizf5e+QMrVixxHaD53L4VNnGH9/M0sCl8hVu2qn6+m3ReQtEWnr66CMMXlDhSsKM7VfAnEVo3hkwlI+/GVjnuu9dOa6fXT5cD7hYaFM6ZdAfOUr/B1SwHFTPjoE6AesBFYB/URksK8DM8bkDcULh/HJfU25rX4ZXvvvOl74YjVp6XkjGUxcuI37xyZxVakiTBuQQLVSRf0dUkBy08XE9UBtZ8xLRGQMsDrrXYwx+UmhAqG817Uh5aIj+PCXTew+msL73Rr6rdvm9HTlrR9+Y/DMjfxfjRgG392IooXc/DozF+Pm0dAGoGKG+QrOMmNMEAkJEZ75S21eujOWn9ftpeuI+Rw4cdr7jjnsdGoaj01axuCZG+nWtAIf3RtvSSCb3CSCYsBaEZklIjPx9DUU6QxfOcPLvsaYfOaeFpUZ1r0xv+05RvshiWzan3vdtB8+eYYeH3leXP+9TS1ebVePsCDvQjonuCkf/b+s1qvqLzkakRdWPmpM3rB022HuH5NEuiof3RtP40q+fUm79eBJeo1axI7DybzVuQG3Nyjr0/PlN9ktH/0l4wSkAZ0zzBtjglDDitFMG5BAVOGC3DViAd+u2u2zcy3Zdpj2QxI5dOoM4/o0sySQw9yWjzYUkTdEZAvwMtYnkDEGqFSiCFP7JxBbNpL+45Ywcs7mHD/Ht6t20234fIqGF2Ba/wSaWHlojsv0DYuI1AC6OdO5juZEVa/LpdiMMQHgiiIFGd+nOY9OXMpLX61h55FknrulNiEh2evKWlX5eM5m/vnNWhpWiGLEPfGUKFooh6I2GWX1qn0d8Ctwm6puABCRx3MlKmNMQAkPC2XI3Y155es1fDxnM7uPJvN25zjCwy6vvPRsWjovfbmGT+Zv5ZZ6pbN1LONdVomgPZ6uo2eKyLd4uqPOP6NVGGNyVGiI8OLtsZSLiuCVr9ey79gCRtwTT3SRgpd0nMMnzzBg3BLmbTrIA/9Xlb/fXCvbdxcma5m+I1DVz1W1K1ALmAk8BpQSkaEiclNuBWiMCSz3X1OVwXc1YsXOo7QdMpe1u4+53nfDvuO0HTKXxVsP81anBjzzl+w/YjLeuakaOqmq41X1djyjjC0F/u7zyIwxAevW+mWY0Kc5KWfTaDdkLp8v3el1nx/X7KXd4EROnk5jQt/mdGhcPhciNeCyaugcVT2sqsNVtbWvAjLG5A+NK0Xz5cNXU79cFI9NWsaLX6wi5Wzan7Y7nZrGS1+u4f6xSVQqWZgZD7WkcaVoP0QcvKxdtjHGZ0oVC2dcn2a89s06Rs7dzOzfD/DsLbVpXasUIjD79wP88+s1rN97gp4JlXnmllo2mpgfeG1ZnK2Di7QB3gVCgY9U9fUL1hcCxuIZtP4g0EVVt2R1TGtZbExgmvP7AZ6dvpJth04RVTiMEBEOnTxD+egIXrozlutreR0B12RDVi2LfXZHICKhwGDgRmAHsEhEZqjqmgyb3QccVtVqItIV+BfQxVcxGWP85+rqJfnpr//Ht6v2kLjxIKpK0ypXcEu9MlYa6me+fDTUFNigqpsARGQicCeeTuvOuRMY6HyeAnwgIqJ5beQLY0yOCAsN4fYGZa2LiDzGl932lQO2Z5jf4Sy76DaqmgocBf40MKqI9BWRJBFJ2r9/v4/CNcaY4BQQ/bc6lUrxqhofExPj73CMMSZf8WUi2IlnEJtzyjvLLrqNiBQAiuN5aWyMMSaX+DIRLAKqi0gVESmIp7uKCweymQHc63zuCPxs7weMMSZ3+exlsaqmishDwHd4ykdHqupqEXkJSFLVGcDHwCcisgE4hCdZGGOMyUU+bVCmqt8A31yw7IUMn1OATr6MwRhjTNYC4mWxMcYY37FEYIwxQc6nXUz4gojsB7Ze5u4l8Yy2FkzsmoODXXNwyM41V1LVi9bfB1wiyA4RScqsr438yq45ONg1BwdfXbM9GjLGmCBnicAYY4JcsCWC4f4OwA/smoODXXNw8Mk1B9U7AmOMMX8WbHcExhhjLmCJwBhjgly+TAQi0kZEfhORDSLy9EXWFxKRSc76BSJSOfejzFkurvkJEVkjIitE5CcRqeSPOHOSt2vOsF0HEVERCfhSQzfXLCKdnb/r1SIyPrdjzGku/m1XFJGZIrLU+fd9iz/izCkiMlJE9onIqkzWi4i85/x5rBCRRtk+qarmqwlPB3cbgapAQWA5UOeCbQYAw5zPXYFJ/o47F675OqCw87l/MFyzs10xYDYwH4j3d9y58PdcHVgKRDvzpfwddy5c83Cgv/O5DrDF33Fn85qvBRoBqzJZfwvwX0CA5sCC7J4zP94RnB8iU1XPAOeGyMzoTmCM83kK0FpEJBdjzGler1lVZ6rqKWd2Pp7xIQKZm79ngJfxjIWdkpvB+Yiba+4DDFbVwwCqui+XY8xpbq5ZgUjnc3FgVy7Gl+NUdTae3pgzcycwVj3mA1EiUiY758yPiSDHhsgMIG6uOaP78HyjCGRer9m5Za6gql/nZmA+5ObvuQZQQ0Tmish8EWmTa9H5hptrHgh0F5EdeHo7fjh3QvObS/3/7pVPu6E2eY+IdAfigf/zdyy+JCIhwNtATz+HktsK4Hk81ArPXd9sEamnqkf8GpVvdQNGq+pbItICzxgndVU13d+BBYr8eEcQjENkurlmROQG4DngDlU9nUux+Yq3ay4G1AVmicgWPM9SZwT4C2M3f887gBmqelZVNwPr8SSGQOXmmu8DJgOo6jwgHE/nbPmVq//vlyI/JoJgHCLT6zWLSEPgQzxJINCfG4OXa1bVo6paUlUrq2plPO9F7lDVJP+EmyPc/Nv+HM/dACJSEs+jok25GWQOc3PN24DWACJSG08i2J+rUeauGcA9TvVQc+Coqu7OzgHz3aMhDcIhMl1e8xtAUeAz5734NlW9w29BZ5PLa85XXF7zd8BNIrIGSAP+pqoBe7fr8pr/CowQkcfxvDjuGchf7ERkAp5kXtJ57/EiEAagqsPwvAe5BdgAnAJ6ZfucAfznZYwxJgfkx0dDxhhjLoElAmOMCXKWCIwxJshZIjDGmCBnicAYY4KcJQKT60SkhIgsc6Y9IrIzw3yij845wemp8XFfHD83iEgrETkqIt9ksn60iHTMwfM9LiLbROSDnDqmyZvyXTsCk/c5de1xACIyEDihqm/66nwiUhpooqrVLrKugNPfVKD4VVVvy40Tqeo7InIYT5ckJh+zOwKTp4jICednKxH5RUS+EJFNIvK6iNwtIgtFZKWIXOVsFyMiU0VkkTO1vMhhvwfKOXcc14jILBH5j4gkAY+KSGUR+Vn+N1ZDRefYo0VkqNN52yYnppEislZERmcS/xYRGSQiS5w4aznLrxCRz51zzBeR+s7ylSIS5bQSPSgi9zjLx4rIjV7+rEREPhBPX/0/AqUyrHvB+fNYJSLDnW2vEpElGbapfm7e+fM9N16Fz5KyyZssEZi8rAHQD6gN9ABqqGpT4CP+18Pku8A7qtoE6OCsu9AdwEZVjVPVX51lBVU1XlXfAt4HxqhqfWAc8F6GfaOBFsDjeJr2vwPEAvVEJC6TuA+oaiNgKPCks2wQsNQ5x7PAWGf5XKClc8xNwDXO8haAt8dk7YCaePrgvwdIyLDuA1Vtoqp1gQjgNlXdCBzNEHcvYJSIlHCOFevE94qX85p8xhKBycsWqepup4O8jXi+2QOsBCo7n28APhCRZXh+UUeKSFEXx56U4XML4NxIXp8AV2dY96XTXcFKYK+qrnR6tVydIYYLTXN+Ls6wzdXOsVHVn4ESIhIJ/IpnIJJr8SSOeiJSDjisqie9XMO1wARVTVPVXcDPGdZdJ57R91YC1+NJNOBJlL1EJBTo4lz3UTzjNXwsIu3xdFtggoglApOXZewhNT3DfDr/e78VAjR3vu3HqWo5VT3h4tjefsleGEPG818YQ2b7pGWxzTmz8dwFXAPMwtNZWkc8CeKyiEg4MAToqKr1gBF4OmIDmAr8BbgNWKyqB513JE3xDNJ0G/Dt5Z7bBCZLBCbQfU+GgUiyeFyTlUT+1/Hg3WTjl3AWfnWOjYi0wvP46JiqbsfTZXJ1Vd0EzMHzOGm2i2POBrqISKh4Rqi6zll+7pf+Aefu6Hwlkaqm4OnAbSgwyomnKFBcVb/B8wisQXYu1AQeqxoyge4RYLCIrMDz73k2nvcKl+JhPM/K/4bnG3m2e3O8iIHASCfOU/yvG3SABXh61gRPwngNT0LwZjqexz5r8HTFPA9AVY+IyAhgFbAHT1fOGY3D807g3KO2YsAXzp2EAE9cyoWZwGe9jxoTIJw7iSezWz4qIk/iuQN43sW2PYF4VX0oO+c0eZs9GjImcJwB6mbWoMwNEZmOp8LoXRfbPg48Axy73POZwGB3BMYYE+TsjsAYY4KcJQJjjAlylgiMMSbIWSIwxpggZ4nAGGOC3P8DZfaTbtvGfJoAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig,ax = plt.subplots(1,1)\n",
"ax.set_xlabel(\"Time from now [days]\")\n",
"ax.set_ylabel(\"Apparent distance Mercury/Sun [deg]\")\n",
"ax.plot(times,angles)\n",
"ax.hlines([0.5],xmin=times[0],xmax=times[-1])\n",
"ax.text(times[0],0.6,\"Apparent size of Sun\");"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Mercury is closest to the Sun's center in 17.945946 hours.\n"
]
}
],
"source": [
"print(\"Mercury is closest to the Sun's center in %f hours.\" % (times[np.argmin(angles)]*24.))"
]
},
{
"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.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment