Skip to content

Instantly share code, notes, and snippets.

@hannorein
Created January 4, 2019 15:44
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save hannorein/284df7266ace9eb40a4027b682416c72 to your computer and use it in GitHub Desktop.
Save hannorein/284df7266ace9eb40a4027b682416c72 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"import rebound\n",
"import reboundx\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"sim = rebound.Simulation()\n",
"sim.add(m=1) # sun\n",
"sim.add(m=1e-6,a=1) # earth\n",
"sim.add(a=1,f=np.pi) # test particle\n",
"sim.move_to_com()"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"rebx = reboundx.Extras(sim)\n",
"rf = rebx.add(\"radiation_forces\")\n",
"rf.params[\"c\"] = 10064.915 # speed of light in AU/(year/(2pi))\n",
"# Set the beta parameter for the test particles\n",
"sim.particles[2].params[\"beta\"] = 0.1"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [],
"source": [
"# Logging and visualizing \n",
"sim.automateSimulationArchive(\"test.bin\",interval=0.1,deletefile=True) # a lot of outputs (just for testing)\n",
"sim.integrate(2.*np.pi) # one year"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rein/git/rebound/rebound/simulationarchive.py:127: RuntimeWarning: You have to reset function pointers after creating a reb_simulation struct with a binary file.\n",
" warnings.warn(message, RuntimeWarning)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXlY1NX+x19nZmAGhh1BQFBcchfJNNvMyiWLSsutbLHF9nv1ttyya9dLq3Zvv4qyzdtmdU3NFk3MMk3NzFxzxVxwAQVEQHaY7fz+mAEBQUC+w3pezzMPM2fOnHMYmPd8zvksXyGlRKFQKM6FrqkXoFAomj9KKBQKRa0ooVAoFLWihEKhUNSKEgqFQlErSigUCkWtKKFQKBS1ooRCoVDUihIKhUJRK4amXkBNtGvXTkZHRzf1MhSKVsvWrVtPSSlD6tK32QpFdHQ0W7ZsaeplKBStFiHE0br2VVsPhUJRK0ooFApFrSihUCgUtaKEQqFQ1IoSCoVCUStKKBQKRa0ooVAoFLWihEKhUNSKEgrFOUlMTmTk4pHEzIth5OKRJCYnNqifomXSbCMzFY1DYnIiCdsSSC9MJ8wcxrQB04jrElf+XPyGeErsJQCkFaYRvyEeoLxPffrVNp+i+SKaaxXugQMHShXC7V6qfsABTHoT8ZfFE9cljhGLRpCanoo1x4ot34Y934690I5XqRejQkdRVFSEw+Fg5ZGVFFmLABBCIDycN19vX+7qfxdms5nw8HBSRApfnPgCh68DnbcOIUSl+apbnxIV9yGE2CqlHFiXvsqiaMMkbEsoF4mCpAJKU0spTSvlzll34lfox7Hjx9B56fAI8EDvp8fga0Bv1mP3tdOrVy+8vb3R6XSs/nU1PsLHOagEaZNIq6TUWorBYODUqVPs3r2bxO2JFGYVYjttQzoknqGeeIZ68uiiRzl23TG6detGbGwsISEh9bJSFO5HCUUbJr0wvfx+zroc9CY9xggjxgFGfnroJx74/QFOWk+e9bpwczh/GfeX8sfzfeaTVphWbb+Z42aWP46ZF4PEacHai+1YTlqwnLRQerKU7du3s2jRIrZv305QUBBF4UWIjgKvrl54dfJC56mjxF5CwrYEJRRNgBKKNkyYOaz8Ax71YFR5e7g5nK5du/K4eLzarcm0AdMqjTNtwLQ69as4n95Lj1cnpwiEm8OZO24uAA6Hg/379zPitREUJReR+3supSdK8ershbmnmcKehZTcUILJZCofV21R3I8SijZMbR/wsg9bbR/Cuvari6DodDp69uxJr5G9ykXFXmyn6EARhfsKyVqcRcicEC6//HLGjRuHd6w3r+97XW1R3Iw6zGzjNPa3cV3nO9dB65B2Q1ixYgWLFy/mm2XfYOxixH+QP34D/DD4Ob/7ws3h/DjuR7f9Hq2B+hxmKqFQNFvqIip95/Ylb2ceeZvzyN+Vj7mnmcArAvHt78vu+3Y30cpbBo3u9RBCfATcAJyUUvat5nkBJADXA0XA3VLKbVrMrWi9xHWJq9W6iQiKQAwS+A/yx15sJ29zHqd+PEXaJ2lM3TGVe++9l9jY2EZacetFq8jMT4BR53j+OuAC1+0B4F2N5lU0hJ2L4PW+EB/g/LlzkTZ9G5FpA6Zh0jsPNvVeegKvDKT3s715d+m7BAcHc9NNNzFs2DB++OEHlh1apqJHzxNNLAop5TohRPQ5uowGPpXOfc5GIUSAECJcSnm2T03ROOxcBN9NBWux83FuivMxQMyE8+u7cxGseh5yU8E/EobNPHssjTnnQepV8I9//IMFCxbw4NQHOVl8kuDrgvG/2F8detYTzc4oXEKxrIatxzJgtpRyvevxKuBpKWWNhxDqjKIB1OUD+3pfCjKPkZbvIKNQcrpEklsCuboAci/6KwUFBeh0Ojw8PPDc9DYellw89QJfI/gbBf4mgX9wGIF/WUlUVBQeSd9UFhMADy+48U23i0VdGPHlCA5sPMCp5aewZFoIuSmEwMsDifCPaLOHni02MlMI8QDOrQkdO3Zs4tW0UKp8++dmHGPnWw+xo90KdmYZOHDgACdOnCDtyH5sDgj3FbQ36wj0Evgbwd+YRcAFefj4+GC327FarRTmnMbqAIvdQV4p5JZKckskuaVHyP7qatLT04nwkXT1d9A1UMcFwTr6hurpF1pIxE/PIaqzUBrZ8sgoysA3xhffGF8KDxRy8uuTZC7LJHd0LrYxNgyGZvVRaHY0lkXxPrBGSvmF6/GfwFXn2nooi6L+lJaWsu2pXmxIOsFvqXa2ptnJLJT0DdXTP8qPmDteoGfPnkRERBD+zS34lx7Hec5cAf8oeKyKt+D1vs7tRlVcfS0WC0efCOZQjoND2Q72ZznYnWlnV4YDmwP6XTyECy+8kKuuuoqh7bIJXPuPRrc8Ri4eeVb0aEFSAaeXnCbUEUp8fDwTJ05Ep2s7CdVN4h6tRSjigL/g9HoMBt6UUl58rvGUUFShmm9he5+x/PbbbyxdupT169ezY8cOeviVcFmUnsuiDAyM0NE1UIdeJwAB8acrj1fXrUJd+tYgJhm6cHZdNoetW7fy888/s2HNSi4IgqujDQzvYuCqaD0mg6heoDSkpriMf136L0yHTTzzzDOYTCbeffdd+vbt2yaiPRtdKIQQXwBXAe2ADOBfgAeAlPI9l3t0Dk7PSBFwz7nOJ0AJRSUqfFAtdsnqw3a++VOy5LAnoREdGTNmDFdffTWDBg3C57+XnPPb/6xx67oFqK1vHYXH8k9/Nh+3sfqwnR8O2dh10s41nQ3ccIEH179/iPDw8PqvrY6c68Nvt9uZO3cuM2fOZOgtQ0kelIzVw1r+2nNlubZUVMBVK8P+f31Y/ccR5u+2smSflV4hem7uaeDmQR3p+tL+yp3rYyloTR0PUSsK2akiBysO2lh22JMfjuiIjY3l3uF9GGv5Em9x5tu/sX6H9PR0LpxwIVl7swi/Ixy/C/3Kn2tt0Z5KKFoJx48fZ+7cuXzw+gtE+Aom9fVgYl8PInzL9tFVthNlNMFhYZ05h5CV9hjNd999x8cz72HD4QLG9/bgvgs9uLiD3nmW4ubtSRkx82LI35vPiXkn8Ir2IvzOcAw+BgSCnZN3un3+xqLFej3aJFU+1PKaf7I2pz1vv/02q1at4tZbb+WHh7rQ1zvz7Nf6R1Y/ZsyE5iMMVSlbVzVCZgTGjRvHuN1TOJ7nw2c7rUz6uphAk+CvF3sysW8KpnMOrg1h5jBkb0m357uR8VUGB589SMTkCHpc0aMRZm+eKIuiKanw7Wq1S+bvsvKf36xIn3AeeXw6d955J35+fk27nWgKKmxPHFLy/QEbczZb2JoOU/76FH/7298IDQ11m+VU9eCz8M9CTnx4gssuu4yvPvqKwMDABs/RHKiPRdF2fEHNkVXPYykp4v0tFrrPKeDTnVYSrjWy+1E/Hn30UadIgPOf/8Y3naY3LhO8tYoEOD/wHl4A6IQgrrsH39/djl/n/x+5ubn06tWLZ+4bQ9bCv7gERZ6JFtUgtDyuSxzxl8UTbg5HIOg2oBv/W/U/+kT2ISYmhg0bNjR4jpaGsiiaCJvNxqdj/XhhXQk9gnXMHGrksqiynWANZw9tiXNYC8eOHeOlibEs3nGaRwd58ORlRvyMrngQN59jLFu2jHvvvZcXX3yRDsM7tGgXqjrMbOb8/PPPTJ06laCiZF4aCld0rHJU1EiHdi2a+AAO59h5bm0pKw7amDnUyP0DPPDQ69wusvv372d43HCKI4sJvT0UnYfTMG9pLlS19WimHD16lPHjx3PPPfcQHx/Pmq8/4oquvpU7eXg5vz0V58Y/ks6BOj4Z48WKO7z5OslKv3cL+S7V3+1Td+/enW4zu2HJt3D4lcNYc5zxFmU1PVsjSijcRYW07OJXehP/8HgGDBhAv379SEpKYuzYsYj+E9vW2YOWVDjHiA3Ts/JOb16P8+Pvq2yMHTuWtFXvuzUt/pTjFFGPRuEb48uhFw5RfMR50FyxYHFrQrlH3UEFL8WGFBv3LPmTfmHJbJ//Gh2vfaRy3+bsymzOVHGzioAorntiJle/fxMvTp1E/9EP88owI3fHeiDOlUJ/npQVCg69KRRjuJEj/3eEDvd0oMeQ1ulCVRaFO1j1PMVFRTz+QwnjFhXz8jVGFo8z0nHvO029stZFzATnWU78aefPmAmYTCZe7LWfH+/wZs5mC6P+V8TxPIfTtbzqec2mrlgwx3+QP50e60Tap2n0OthLszmaE0oo3MD2fUcZ+N9CUvMc7HzYzNjeHs4nclObdmFthdxUYsP0/D7FzBVRBgbMLWTRHqum739VF2rXfl2Zs3gO3777LQkJre+cQm09NERKyTvvvEP8/4p5Y6SRSf0MldO4a4qkVGiLfyTkpmDQCf451Miobgbu+KaY7w578vbjeWfiUxpIdTU9r113LcOHDyc/P58ZM2acncbfQlEWhUYUFhZy1113MXfuXH5b8Bq3X+RX+Z9EeTMajwoHnQCDOujZ/mg7vKMHMnjwYP5c+obbDjo7derEunXrWLhwIdOnT6e5hh/UFyUUGnDgwAEuvfRShBD89ttvdLthqvJmNCXVRLJ6j32L9xf/xJMThzJk0uN8tykZrSM6ywgPD2fNmjWsXr2aMfeOYcSXI1p8QV8VcHU+VIgaTDzuzz1fnea5l2bz0EMPtRpTs9Xyel827jnM+C+LefAiT2YM8XRbZurCPxZy75h7MceYaT+2PdC8grJUwJU7KXN95qbwzuZS7l+QwpLxeh6+PFiJREsgN5VLIg1smmJmyZ9W7l5SgsUu3XLQ/OHBD+n4REfytuaRudyZ/dtSg7KUUNSXVc/jsBTx1MoS3vzdwvp7zVwabtfU9aZwI64D5XBfHWsmm8krlVz7eRHZhjDNp0ovTMfgZyD679Fkr8om55ec8vaWhhKKemLNTuHeJSWsP2bn13u96RLoeguV67NlUOGg0+wpWDzeiwEdjFz+UT7Hjx/XdKows1N8PAI9iH4imowvM8j7I6+8vSWhhKIelJaWMvYbyclCyco7vQn2rvD2Kddny6DKQac+sCP/987H3PPQNK666ipSfnxHM49IxaAsY4SRjtM6cuLDE4w2j9bol2k81GFmHbFYLIwbNw6PogwWXHkED0fj13NUuJfXnryLOR/N5+e7vOgU4PoSaODftmpB35iUGL5+62s2bdpEu3btNFx9/VGl8DTGarVy6623IoTgi+W/4LHv2+Zbk1Jx3jzeYRuGwQaumlfI2rvNdPTXnQn9Ps+/b3VBWbo0HRMmTOCHH37Aw8NDi6W7HWVR1ILNZuOOO+4gPz+fr7/+GqPR2NRLUriL+ABA8sbGUt7ebGXd3d6E++rQupCQ3W7nxhtvpFu3brz55puajVtflEXRUFxxEo7TKdyTKMjx7MSS1b8rkWjtuEK//3aJkQILjPisiLV3exMc3knTafR6PfPnz2fw4MFMe2UaSV2Tmn2VLHWYWZUKcRLTfyrhSGYB34xIx7R/aVOvTOFuKnhEZgzxJO4CA6Pml1Bw6VOaTxUQEMBjbz3GOy+9w+H9h5HI8iusN8foTSUUVVn1PFiLeWezhaV/2lhyq7fzQjQqTqL1U8EjIoSO2WO70W/gFdw1Yw6O/+ujeW7I13lfEzoulJR3UnBYHEDzDchSQlGV3FRWHLTxwrpSlt/uTZCXKG9XtAEq1LgQj+/h3WenkHnoD2YuOYjWuSHphekEXhmIqYOJ9IXpldqbG0ooqpBUEsJd3xTz5XivM8FUoOIk2ijGX2bx1Xgj/9tlZf4u17VINSqCE2YOQwhBxOQI8rfnk78rv7y9uaGEogJ5eXmM/qKQV671rVwZW6WIt11yUwk161h6qzfTVpSwI91e3t5QygKy9GY9HaZ04PiHxzEUG5g2YFqDx9YaJRQupJQ89NBDXHPdTdwT/1+VIq5w4rIk+7XX8+YoE2MXFXG6RGpiYVaskuXb25fwweGE/hzaLL0eyj3q4uOPP2bXrl1s2rQJvLyUMCicDJtZXij5tn4e/JZq564lpXz7xbOafMtWDMjKHZNL7969+eWXXxgyZIgGo2uHsiiApKQknn76aRYsWICXl1ftL1C0Harkhrw6riuZnp14Y/UJzafy9/fnjTfe4KGHHsJisWg+fkNo85GZxcXFDB48mKlTpzJlyhS3z6do+SQnJzN4YCzrHmhPL1OmpmH8Ukri4uIYMmQIzzzzjAarrRlVuKY2Klyc5+mRkfSO8OG+++5r6lUpWghdCrbw/BUw+bNkbA6Hpi5TIQRvv/02r776KsnJyRqsVhvanlBUiLz89ZiVxX/k8M5FBxG7vmzqlSlaCque56ELIcAk+Pevri2ChtcN6dy5M3F3x3HppEubTa1NTYRCCDFKCPGnEOKgEGJ6Nc/fLYTIFEL84bo1nY3virwssUmmfFfCW9eZCPIoVZGXirqTm4oQgg9u8uK13yz8eUo7lyk4U9P39NlDzoEcCg8UNovQ7gYLhRBCD7wNXAf0Bm4TQvSuputCKWWs6/ZBQ+c9b1x/zFfWW+jZTqcuzqOoPy7XaEd/Hf+80pMHl5U4y/JrFJSXsC0Bi95C6C2hpC9MR0rZ5KHdWlgUFwMHpZTJUkoLsABoviV8/CM5mO3grU0W3hxlqtSuUNSJCsljf7nYk0Kr5ONdQrOgvLIQ7oDLArCX2MnfkV+pvSnQQig6ACkVHqe62qoyVgixUwixWAgRpcG854W85p/8ZYWF6Vd4EuVfoYqRirxU1JUKLlO9Tsd/b4vmmbWC7MjhmgxfFsItdIL2t7Tn5FcnkQ7ZpKHdjXWY+R0QLaWMAVYC86rrJIR4QAixRQixJTMz0y0LWZkRyGFLMNNGdkVFXirOmwrJY7GvHmLM1QOZNba75rU2fS/0RXgIircVN2lod4PjKIQQlwLxUsprXY+fAZBSzqqhvx7IllL6n2tcd8RRSCkZNGgQTz/9NOPHj9d0bEUbZuciTvzvUfq9lcX2B10l9DSstanbo6PohyL+/ONPTa8d09hxFJuBC4QQnYUQnsCtQKUqL0KI8AoPbwKSNJi33nz11VdIKRk7dmxTTK9orax6nggvCw8P9OCfP5c62xroLo3rEseP435k5+SdbJu9DVEiWLdunUYLrj8NFgoppQ34C/ADTgFYJKXcI4R4Xghxk6vbVCHEHiHEDmAqcHdD560vNpuNZ599lpdffhmdru2FjyjciMtj9tTlRr4/YCMpU1t3qU6n44knnuA///mPJuOd1xq0GERKuVxK2V1K2VVK+ZKrbaaUcqnr/jNSyj5Syv5SyqullPu0mLc+fPrpp4SFhTFy5MjGnlrR2nF5zPyMgscv9eT5daWV2rXgrrvuYsuWLezdu1ezMetDm/hqLS0t5bnnnmPWrFnq+qAK7angLn10kCerku3szfHQ1JNmMpl49NFHefXVVzUbsz60bqFw5XR8MTGQHqYsLjWn1P4ahaK+VHCX+hp1PH51CM//2V1zT1q367vx2cLP6P1270YP62699ShcOR3SUsRrv5XynxEmZ44HKFeoQntiJpT/Xz2an090dDRHjx6lUydtSv0nJifyWtJrmGPM5GzIQT9ST/yGeIBGKXTTei0KV07HT8l2HBJGdtVrmrijUNSEr68vd143mPfuHahZ5e6EbQmU2EsIHBpIztqcRg/rbr1C4TpxnrPZwrTBnmfOJlROh8Ld7FzEIyGb+HBjFiU2bdLQy8K3zT3NOCwOSo6UVGp3N61XKPwjOXrawa/H7Ezq51GpXaFwK6uep7u/lQHhOhbu1qZyd3lYtxAEXh5IzvqcSu3upvUKxbCZvLfdwZ0xHpg9XdaEyulQNAYuq/WhgZ58sN16Vvv5UDGsO+DyAHI35eIpPRstrLvVCoWj7zg+TTJy/1BVTVvRyLis1rgLDBzIcrA/y16p/XyoWLHbGGLE3N7MLfpbGq1id6v1eqxfv5524VH0fmVHUy9F0dZwVe72oJg7Yzz4eLuVWdf5NNiarVix++XjL5O2KQ3u1GLBtdNqLYqFCxcyceLEpl6Goi1SIa7i3gs9mbfLge261zW1ZseMGcOSJUtorOLYrVIobDYbixcvVkKhaDpcaei93s4n4oL+rM8Lr/019aBXr14YjUa2b9+u6bg10SqFYu3atURFRdG1a9emXopCwejBXVnyr7GaXg1dCFFuVTQGrVIo1LZD0WzYuYibxE8s2ZmDlNqW9h89erQSivPFbrfzzTffMGGC8m4omgGrnicm2IpDwp5Mh7NNowjhrJAs9uzfQ+857s/9aF1CsXMRW5+4gHBdNp2+jtNEtRWKBuEq7X9jdwOJ+22V2htCYnIiL25+EdMFJgr2Fbi9pH/rEQpXEthPu04wvItBUxNPoThvXLETI7saWJlsO6v9fCnL/TD3NFO4rxDArbkfrUcoXElgK5NtDO+id7apJDBFU+OqVXFVtIHfj9spskpNIoTLcz+6myk6UHRWu9a0HqHITaXEJtl03M6QjoZK7QpFk+GKqfAN7ciFYXp+ORWkSYRwWY6HKdqEJcOCvdheqV1rWo9Q+Eey5YSdXu10+BpFpXaFoklxxVRcc8+z/Bx4qyaBV2W5HzqDDlMnE8XJxZj0JrflfrQeoRg2k19SdZWtCZUEpmhGXHzxxWzdulWTsSrmfnh19MIzw5P4y+LdlvvReoQiZgIbSrtzRc/2qCQwRXPkIq9Utv76M/Jf/poEXpWV9J81YRaX6S5za4JYq0oK230sm/4froZu3Zp6KQpFZXYuov2GmXgbHBw+7aCLSNGsNGO/fv14//33NVhkzbQai6KwsJCMjAw6d+7c1EtRKM7G5ZW7KELP1hOutHONvHJ9+vQhKSkJu93e4LFqotUIRVJSEt27d0ev1zf1UhSKs3F53/qF6tiTaT+rvSH4+vrSvn17Dh061OCxaqLVCMXevXvp3bt3Uy9Doagel/etT4ie3ScdZ7U3lL59+7Jnzx5NxqoOJRQKRWPgCrzqE6o7k/OhkVcuMTmRHfYdPLLoEbflfCihUCgaA1fgVY8uHTmU7cDuG6mJVy4xOZH4DfFY/C1Ys6xuy/loNUKRlJREr169mnoZCkXNxEzA+Pe9+AUGk3Xnz5q47styPjyCPbCcsgDuyflo+UKxcxHytT6kHjlI1JJbVBKYotnTvn17MjIyNBmrLLfDw98DW67trHataNlC4coYzTuZgkEHPiUnVMaootkTFhZGero2H+Sy3A69vx5bnu2sdq1o2ULh8k2nFzgI83H9KipjVNHM0dKiKMv5MPgZsOXZkFK6JeejZUdmunzQJwsloWZxVrtC0Rxp3769ZhZFWdj2G1vfYB/7CNWH8sRlT2gezt2yLQqXDzqnRBLkpTJGFS2DsLAwzSwKcIrFyvEr6RTRibmXznVLzkfLFgqXbzqnWBJgUpcNVLQMtNx6VCQkJIRTp05pPi5oJBRCiFFCiD+FEAeFENOred4ohFjoev53IUS0FvOW+aZzCCDQpDJGFS2DoKAgcnJyNB/X09MTm81We8fzoMFnFEIIPfA2MAJIBTYLIZZKKfdW6HYfkCOl7CaEuBV4BdCknn6ij5l3fH3JKbWxLyqCaT5mGudqjK2b6OlnB+wcma3e2Qaz7HE8Fv0X275ieC4ILrobbnhNk6F1Oh0Oh6P2jucztgZjXAwclFImSyktwAJgdJU+o4F5rvuLgWFCCEEDKYtKy7fmI5Fur0TcVqhOJM7Vrqgjyx6HLR9iEA5sDkDaYcuHznYNEEI0a6HoAKRUeJzqaqu2j5TSBuQCwQ2duCwqDQDXJRjdWYlYoWgQWz8BwKDDKRRV2htCYnIiu7J2ce+Ke92S79Gs3KNCiAeABwA6duxYa/+y6LPCg4XY8+1ntSsUzQp55n+00rWFZcPqSJTnezgsSHnGsgY084BoYVEcB6IqPI50tVXbRwhhAPyBrKoDSSnnSikHSikHhoSE1DpxWfSZtEkcFsdZ7QpFs0I4a6XYpdOqqNp+vjSGZa2FUGwGLhBCdBZCeAK3Akur9FkKTHbdHweslhpcr70sKk0IAS6dcGclYoWiQVx0NwBWu8RDf3b7+VJmQUurRHiIs9q1oMFC4Tpz+AvwA5AELJJS7hFCPC+EuMnV7UMgWAhxEHgcOMuFej6UVSI2eZiQUhJuDndrJeK2Qk3eDeX1aCA3vAYD78MmdRh0wmlJDLyvwV6PMgvaYXWg89Sd1a4FmpxRSCmXA8urtM2scL8EGK/FXFWJ6xLHiAtGsDVvKz+O+9EdU7RJlCi4iRteo6T4UoyWhfCvxZoMOW3ANOI3xCMtZywKrS3rZnWYeV7sXET4iRXkHc9xlkAfNlMFXCmaNadOnSI4uMFOv3LKLOhxtnHoPfWEm8OZNmCappZ1yxYKV5p5lLGQQqs8c2FiUGKhaJ7sXETmkn8RWpQHr/+q2RdbXJc4/IQfa+5YQ1iY9of5LTvXw5Vm3jlQR0lZ5KpKM1c0V1xfbJnZpwkxc+aLTYP6KVJKCgoKMJvNDV9nNbRsoXClk18QpMNqP7tdoWhWuL7YThZJQry1rZ9y+vRpDAYDvr6+DR6rOlq2ULjSyS8IFkggr8RRqV2haFa4vsCO50k6+GlbP+Xo0aN1ClI8X1q2ULjSzA06HToBezMdKs1c0XxxfYGl5jmI9NOd1X6+JCYncvcXd3NUf1SV668WV5o5/lGYDLC3IEClmSuaL8NmYtOZSC+QdPDVpn5KWfj2yeMn8Qj2UOX6ayRmAjy2m6D2kWwOukmJhKL5EjOBg/2nExngidGg06R+Sln4tiXLgkewB+CexMiW7R51kZicSKFfIfNWzuPQ4kOa+5AVCq3YY+tIn8tHQfwSTcYrC9O2nrLiFe11VrtWtHiLosz00kXpsGRaVE0KRbNmz5499OnTR7PxysK0S9NLMYYZz2rXihYvFGWml7mHGXu+HemQqiaFotmyZ88eTS99OW3ANIzCiOWkBc/2noB7EiNbvFCUmVje3bwBKM0srdSuUDQntLYo4rrE8XDHh/H08cTgZXBbYmSLP6MIM4eRVpiGwc+A0AuKkoowtTepmhSK5sXORRR9H8/h/Un0Wnk76OM1O3j3PenLqCGj+Hbyt5qEJ/qJAAAftElEQVSMVx0t3qIoq0kBYAg0UJBUoGpSKJoXrtDtLfuO0i9Uj6nouKaXvtyyZQsDBw7UZKyaaPFCUVaTItwcjinShPWYVdWkUDQvXKHbG1LsXBrpqlijYU7S5s2bGTRokCZj1USLFwpwisWP437k9ZsmYc0oZuTHk5wp5+pixYrmgCtE+5djdq7oqD+rvSE4HA62bdvGRRdd1OCxzkWrEAoAdi7CyDJsQA9HECN97ST+9HclFoqmxz8Sm0Oy/piNodH6Su0NZf/+/QQHB9OuXbsGj3UuWo1QJP7yPAmd/NGZdOT+nkeah4H4QB8Sf1Ep54omZthMtmV6EOWno11Z1qgGodsjF49k+KvDKQovcnvcUKsRigSjnRKdDu+u3uTvzAegRKcjwdiwUugKRYOJmUCivJpRfQKBhl/6sizIMK0wjYK9Bei7690eZNhqhCLd4DTp/C/2x5plxVZgq9SuUDQl3/52gDEvfAvxp+Gx3Zrkd0gpKUwqxNzb7PYgw1YjFGGeAQD49ncW7sjfkV+pXaFoKg4fPkxaWhqXXnqpJuOVBROWnihF6AWeIZ6V2t1BqxGKaZc8g0l4YPAz4BHgwen1pzEJD6Zd8kxTL03Rxlm6dCk33ngjer021m1ZMGHBzgJ8Y3wpu4yvO4MMW41QxHWJI/6KFwg3h+N3kR9FB4r4e//pKp5C0eR8++23jBkzRrPxyoIM83fn49PPB3D/ha9ajVDAmXiKr2Z8hZ/Zj/wt+U29JEVbZecieL0vGU/6sX3jOoa3z9Vs6LgucTzR+wlKk0vx6e3TKBe+avG5HtVxxRVXYJM2Hnv5MRIMCYSZw1SNCkXj4QrZxlrM/F0WxvTQ47XySTB5apbfUbi9kNHXj2bxg9pcRKg2WpVFUcaKoyvwvtibgtQCSk+WqhoVisbFFbItpeSTHVYm9/fQ/DISCxcuZOLEiZqNVxutUigStiXgO9gXnYeO0xtOA+4pD6ZQVIsrNHtrmoO8UnkmGlOjy0hkZmayadMm4uIaz0JulUKRXpiOV1cvhEGQvTabsgunqxoVikbBFZr9/hYL9w/wROfySmh1GYmvvvqK6667Dm9vb03GqwutUijCzGEInSBgSADS4gxKKWtXKNzOsJnk2Y0sTrJyT6yz4K2Wl5FYsGBBo247oJUKRZn7KPCKQKRdkvl9JiaHg2kZJ1SSmML9xExgbtG1jOrlR7ivXpNq2+AM3b78rcv5dfuvzCme06hnbq3S61Hm3UjYOIvjHVIpPVjMPbsziPO1qosYK9xOaWkpry9YTWLiWoiN1WTMsvyOI98fIeDyADIsGcRviAdoFG9eq7QowBVTkXGat/pIunjB5tV5zifURYwVbuazzz4jJiaGWI1EApwH9EUlReT8mkPg0ECgcQ/oW61QAJCbym19PUgvlXz+p42ePu0ZGRlBoi27qVemaKXY7Xb+/e9/M336dE3HTS9MJ3dTLl4dvSqV5W+sA/rWLRT+kawOMGO+JhhDqJHM5VnOOhUhQSqmQuEWvvnmG4KDg7nyyis1Hbe9d3uyf8omaFhQpfbGOqBv3UIxbCYJQYH4DQvGctLC6d9PY822UiKEiqlQaIcrXNs+058Xpt7BM7cNKU/U0ooRcgSOQkd5djS4P7+jIg0SCiFEkBBipRDigOtnYA397EKIP1y3pQ2Zs17ETCDdoMcjwAO/gX6YOpjI+CYDUDEVCo0oC9fOTWH+LgtmvY0b8+Zp7l1b/dFq7p96PxG+EQhEo+R3VKShXo/pwCop5WwhxHTX46er6VcspdTuZKcehJnDSStMI+T6EJJfSqY0o5Si5CK69uvaFMtRtDZc4dolNsk/fy7ls5u9ELYSZ7tGnrVNmzaxd+9elixZgtForP0FbqChW4/RwDzX/XmAdrm0GlEWU2EMN2LuZcbc00zG/Az+GvvXpl6aojXgCst+Y6OF2DA9QzoZKrVrwQsvvMBTTz3VZCIBDReK9lLKNNf9dKB9Df1MQogtQoiNQogaxUQI8YCr35bMzMwGLs1Jxet+hN4QSumBUsKMYZz+7bQm4yvaOP6RpOU7eHWDhVdHmiq1a8H27dvZunUrU6ZM0WS880WU5UHU2EGIn4DqjlZnAPOklAEV+uZIKc86pxBCdJBSHhdCdAFWA8OklIfONe/AgQPlli1b6vI71ItbbrmFiIgIvvnmG/bt24evr2/tL1IoamLnIiZPnkyol4P/lAmFh5cmkZjg/H8dMmQIjz32WIPHqooQYquUsk6XGKv1jEJKOfwcE2UIIcKllGlCiHDgZA1jHHf9TBZCrAEuBM4pFO5i1qxZXHzZxXj19KLrrV2JvSdW1apQnDcrTvix9qQvu/8aDCVpTkti2MwGiURiciIJ2xI4tOMQx9ce5/YXbtdwxedHQw8zlwKTgdmun0uqdnB5QoqklKVCiHbA5cC/GzjveXPQ4yDGC41YrBby1uVx5IojxJfEA40TCqtoPeTl5fHggw/ywSf/w2fECE3GLAvVLrYVk7YgjXY3t2P29tmYvExN+v/Z0DOK2cAIIcQBYLjrMUKIgUKID1x9egFbhBA7gJ+B2VLKvQ2c97xJ2JZA8Ohg8rblETQsiNT/plJsKSZh9ZPqMoSKejF9+nSGDx/OCI1EAs6U4s/bnIe0SAIuD2gWtVQaZFFIKbOAYdW0bwGmuO5vAPo1ZB4tSS9Mx+BnICQuhIKkAnTeOjKXZaK7KQRyU1TSmKJOrF27lqVLl7J7925Nx00vTMdR6iB9YTqR90cidKK8vSlp3ZGZ1VAW8ho0IgjLSQsBlwSQ9VMWXgcLnB1U0piiJlwRmAX/8GfKuJG889QkAgK0vW5MmDmMzOWZeHf1xtzTXKm9KWlzQlEWV6Ez6AifFE7md5lETWxPynup5Je6PEAa+sAVrYQKEZiPLi/mikjJTfmfar5VvcXvFnJW59B+4plIg8YM1a6JNicUFeMq/Pr5EBDpycBDuVwXpeMv35c4O2nkA1e0IlwRmPP+sLDlhJ0515k0tz7tdjufxn/K/U/eT6eOnZokVLsmWmXhmtqI6xLnfON3LiJj5KP0nJNNt4c6sOt/GWw9ZuaVu25G+T8UlchNZWeGnSdXlvLzZG/MnqK8XSveeustPDw8mPPsHHS65vUd3iaFopyYCWzJ3kXQ6U/YNf8kHR6MZN//HeXJyMXQ8ZImV3FF8+G0ZzhjF+3njWtN9A2tcGlAjazPQ4cO8eKLL7Jx48ZmJxLQBrceVUk49Tteg/0xRhjJ3ZhL5P2RHHzrILO/n93US1M0E+x2O3f+4M113Y3cHuNx5gmNCuY6HA6mTJnCP/7xD7p169bg8dxBmxeK9MJ0hBBETI4gd2MuQi8IiQth88ubyc3V7jJwipaJlJKpU6dS7BHIq3M+cBbKRWhSMDcxOZGRi0fS4fYObD+xnW43Nk+RgLa+9cDpdkorTMPga6DDlA6kfpBKtxe6YThlYOLEiSxbtgyDoc2/TW2W2bNns379etatW4envz9cNEmTccsiMLMPZJOZmEnXmV154fcX0Ov1zXLL2+YtijJ3KYBvX1/8L/Yn/ZN03kp4C4C//e1vTbk8RRPy2Wef8f777/P999/j7++v6dgJ2xIozCsk5d0UIu6KwDPEs1lEYNZEmxeKiu5SgSB2ciwhpSF8+9G3lN5WygfffkD326JJfLsvxAeoMO82wsqVK3nyySdZvnw5ERERmo+flp9G6txU/GL98B90RoSaOgKzJpRNTQV3qYtPoj7hgRsfIOK+CKKfjObwrMM8RCDvDfAiToV5t052LnLGROSmsul0ELd/msVXS5bRu3dvt0xX/EMx9iI7YRMrR1w2dQRmTbR5i6I65qfPJ/LhSFL/m4rD4iD66WgyVmXzxHZX5KYK825dVIi63Jhq5Yb/HuGjOB1D/NNqf+15sGLFCnLW5NBtajeE4UwR3uYQgVkTSiiqIb0wHXMPM6FjQjmWcAy9SU/npzqTvDKbNzaWOjupMO/Wgyvq8rcUGzd9UcwnY0zc0NXhli+Dw4cPM3nyZL798ltevv7l8i1vc4nArAm19aiGMk9I8DXBlJ4o5dg7x4h+LJrBT3bk7X8fptgKz9zQfF1ZinqSm8pPyTZu+6qYz272YlQ37eteAmRlZXHDDTcwY8YMhgwZArScGijKoqiGip6Q8Enh6Aw6Mj49ztOigLV3e/PpLhvP7O6Gw+Fo4pUqtGBJih+TvirmqwkVRAI0zfkpKCggLi6OG264galTp2o2bmOhhKIaKnpCdDodA//aG99DJWz//jQRkZ1Yt+g91u3LZMKECRQWFjb1chXniZSShIQEHlqaz/LJgVzZqYJIaBB1WRZQ1feDvnS6vBP+nfyZPbtlRvyqrUcNVPWEpMelM3ToUI6euoRjpxaRe08uhz8/TP/B/VmzYg2RkSrjtCVhsVh45JFH2LRpE79t3k503qZyr4dWdS/jN8RTbC0mZW4KeEDW9VksP7y8xWw3KlJrFe6mwl1VuBvCJ+s/4cGbHyT4+mCCrwlGSsnpFacpWVvC8iXLufjii5t6iYo6cPLkScaOHUtwcDCff/45Pj4+ms8xcvFITuSf4MQnJ7BkWuj0WCd0njrCzeH8OO5Hzec7H+pThVttPerB/PT5RD8VTeayTLLXZiOEIPC6QCLviiQuLo4FCxY09RIVtbBz504GDx7M0KFD+frrr90iEuAMqDrxyQlK00vpOK0jOk/nR625BlTVhtp61IP0wnQ8Qz3p/FRnjvz7CDgg6OogHL0d/PTTT4wePZotW7bw8m0D8Fz3smZmrOI8qRBEhX8k80uv4W+vLeDNN9/k1ltvddu0DoeD7M+zKU0vpdPjndCbzqSlN9eAqtpQFkU9KPsjG8OMRE93WhZZP2URZg6jf//+bNmyhQNb13LpTZPZn3wEkGcK9qqw78alQhBVXqmDuz7ZzwuvvcOPc550u0g88MADBOYH0uPvPSqJRHMOqKoNJRT1oKLb1BhqpPP0zmStyCLstzCWHVrGpDWTSL6riKKrgxn0cRHvb7EgpVSRnE2BK4hq/TEbF75fgJcBtkzxJvb4p26b0mKxcM8997B//342rt7IC9e80GICqmpDHWbWk7KrOKUXphNmDuOODnfw3JTnyAvNo/2d7RF6Z0iuTCnG8v5Roj0lH97kRZS/HuLV9U4bi+IZ/vzz5xLm77LybpyJ0T3LCs4It/wdsrOzGTt2LH5+fsyfPx+z2Vz7i5qY+hxmKqHQgGs+u4bN/94MOoh6OAq9l9PcbF9i5eqFR3jjdwsvXh/GlAUp6PX6WkZTNJTVq1fz6G3XERPi4O3rTbTzrmA4+0fBYw2/FkfFLwzfPF+Ovn6UW2+5lVdeeaXF/I2V16OROeU4Rae/dcIjwIPDLx/GkmkB4KTRwIwrjay6N5DPD/lz0UUXsWbNmqZdbCvm8OHD3HLLLdx33328PP2vLLwtuLJIaFS6rixGIq0wjYI/C9j07CYMVxm4+pGrW4xI1BclFBoQZg5DGAQR90QQMCSA5BeTKfyzkDCbHfyjiLn/XdZt3cuMGTO4++67ueWWW0hOTm7qZbcaCgoKePbZZxk0aBAXXXQRSUlJ3PzYq85SdRqWrisjYVsCxbZisldnc2zOMSIfiMR3qG+zLTqjBco9qgHTBkwjfkM8JfYS2o1shzHcSMqcFOKefgj5jzdZfng5CV9dS3pROl1f6orvdl8GDRrElClTmDFjBn5+fme58pRLtRqqvEeOq5/liz12pk+fztChQ9mxYwcdOnQ40z9mglvewxM5Jzj+2XGKk4vpMqMLxjAj0HJjJOqCEgoNKDvJLtuzdr+kO8+Meoa3H3+boRuGkj86H5vRBsBJ20nyYvN4fdnr/Dz3Z3r06MHfbx/Bg+YfMAvXBYhUcZyzKXN3WotxSMlXG5J58d93YgzuxMKFC7nssssaZRlHjhzh2Kxj6EJ0dPlnl1YRI1EX1GGmGykuLqb7zd3J2JZB1INReHfzLn+uLJT3jz/+4KW7h7L2QD5/vdiTv1zsSaCXq5iJRgdvrYLX+2LNPsbCPVZmrbfg4wkzrzRy/cDOiMf3NMoSEhMTue+++7hxyo1s77GdUkdp+XMmvanFuT/rc5ipLAo34uXlReBtgei66zj65lGChwUTEheCMIhyMzU2NpYvxwj2nfJm9noLXd/M544YT6YN9qQrqjgOQH5+Pv9dcYg3NpbSLUjHayNNjOyqRwgBecc1n6+qC/y+bvexImEFa9asYdGiRVx55ZVn9Zk2YFqLEon6ooTCzYSZw5AXSby6eHH8w+Mkv5RM5AORRHeLPtPJP5KepPDJGC+O5xmZs8nCJR8WcnkXH+4dsJRRo0bh6enZZL9DUyCl5JdffmHevHl8/fXXjIz25OuJegZGVPEqaHyd2DKPRonduQ3c/9t+7nzwTkZcN4Jdu3aV54ZUzS5u7Sivh5spi+b0CPSg0xOdCLgigMMvHabTH52w2WzOmgXtA4iJjmJkZAR/hPswa7iJI0+GEHfLbbz66qt06NCBhx9+mF9//ZWztoo7Fzkrg7ekCuHnWPPhw4d57rnn6NatG4888gi9evVi7969LJz3XwZ2qpLApZG7syIJ2xIosZdgK7CR+mEqJz49QYcpHbCOsbotgawloM4oGoGqZuq4wHEsmrWIQ8cPYRxvxKPrmcvUmRwO4osEcUPOeD2OHj3K/Pnz+fzzzykqKmLSpEnccccd9LLuKj/gK8fDSzM3oFuocChZRnqJJ4mmsXy+9k92797NbbfdxuTJkxkwYIBze1HxtW72DPX7pB85G3LIWJSB30A/2o9rj95Lj0Cwc/JOTedqalRkZgtASkn/x/uTNC8Jn94+hE0Iw+Dv3AnWVLNASsmOHTv43//+5wwTtp5iZLTk2q4Gru5swMezBRyCug4lN5+w88NBG4kHbBzKcTCihx+3zvyYuLg4jEZjkyxt3759XDbuMooLiomYHIF3l7MPn1sTjSYUQojxQDzQC7hYSlntJ1sIMQpIAPTAB1LKWuuBtXahAIiZF4Ot2Ebmkkxy1ufQblQ7gkcGo/fU1/rt5XA42PmIPz8esvLDIRu/p9rp117PlR31DI02cMlbKQQFBZ17AQ35hq7Ha3NyctixYwfbtm1j9XtP88sxG10CdQzvbCCuu4HLo/R46HWNkgtT3SHkJf6XMHv2bD7++GPGPjKWrV23UkrL9mjUhcb0euwGbgHeP8di9MDbwAggFdgshFgqpdzbwLlbPGHmMNJII+zWMAKvDiRjUQYHnjlAj9t74LjTgU6nq/F0XafTEdujE7FhKTx1uZFiq2Rjqp21R23853fYGh2Nj48Pffr0Kb/17duX3r17Oy+PV3ULUJ/YjRpea7c7OOI3iB07dvDHH3+U/8zOziYmJobY2FjuujSMj0efJsRc5XhM40PJ6qh6UJl6KpUHnnyAgjUFTJo4iZ07dxIREdHmPBp1QZOthxBiDfBkdRaFEOJSIF5Kea3r8TMAUspZ5xqzLVgUVf9xAawHrDiWOvASXoy6bxSrAlZRKmv4dqtmv192RiH7jSclJYXdu3ezZ88e9uzZw+7du0lKSiIwMJAO+mzCvayE+QhCzYJAkyDQSxAYHILXxA+w2+04HI6zbjabjewvHyMzK5tTRZKMQgepeZLUPAdpBdA+vAOxsbHExsbSv39/YmNj6dKlCzqdSxjOsWZ3n6uMXDyStMI0HKUOslZlcer7U/j286XPpD6sf3S9W+dujjS3OIoOQEqFx6nA4Oo6CiEeAB4A6Nixo/tX1sRUjegMM4cx7d5pXP/C9SxfvpzbH7ud4vxiQuJCCLg0AGEQ5ReyjesSd+aDVc0WQOB8Dzt27Mj1119fPqfD4eDYsWOkvdCX9AIP0gucH/ajuZI/MhzkJKVRfOw/6PV69Ho9Op2u0k2v1xN4OIsQM3QNFFwS6UGknyDKT0cHPx2eL6RU85tW4BxrdjfHTx7n1KpTZK3MwtzDTOfpnTF1MJFHntvnbunUKhRCiJ+A6mJTZ0gpl2i5GCnlXGAuOC0KLcdurtTkj4+LiyMqM4qCpAIyv8sk45sMgocHEzQ0iHQq5BTUM59Bp9MRHR1NdN9o55ahKv5R8NjKcw/yet+aX1sX3JSDUdOWISMjgzfeeIMDcw5g7m+m8zOdMUWYyl/XmkOvtaLWOAop5XApZd9qbnUVieNAxf+gSFebohbCfcLx6e1D56c702lqJ0pSStj/1H5OLzjN/v37y/uVXT8iZl4MIxePJDE5sfbBh810mvwVqWtcQkNe6yYqpn5LJCfyT/DY+48xJG4IPXr0IDc3l/cS36PbQ90qiURLLk/XmDTG1mMzcIEQojNOgbgVmNQI87Z4KmalekV7EfVgFLpcHT329eCKK67gggsuYMD1A1jfbj02kzPpLK0wjfgN8UAtl6tryBagCbcPNVEWKGXNtpKzPoecX3LQe+nRj9Rz5MgRAgICAGif3F4dVJ4HDXWP3gy8BYQAp4E/pJTXCiEicLpBr3f1ux54A6d79CMp5Uu1jd0WDjPrQk3mtNVqZcWKFUx5aQqndjoP5fwH++PTz6fZXT9CC87liTh9+jR9H+/L6d9OU5JSgt8gP4KGBmGKNqETulYXKKUVKuCqDREzLwZrgZW8TXnkbs6l+Egxvv198R/kT9LsJLy9vVu8u68675AuT8fVxVdzbOMxVq1ahamXCeMgI74xvuXX0IDWGSilFc3N66FwI2WxGEHXBBF0TRC2XBu5W3PJX51PaGgoXft1JSsqC2MvI17RXjVuTZpKTOoyb8K2BIqKiig6VETBrgIKdhdgzbaSEZPBrPtn8dFHH/Fr9q9niYk6f9AOZVG0cKr7ti2Ltbgy5EqufPlKUralULC3AGuWFe9u3nhf4E1UTBTr/r6u3OKoaYyqH9q6Ckpd+p1r3iHthvDbb7+xdu1aEr5MoPhYMaYoEz59ffDt54tXZy90+srbipZuOTU2auvRxjjXByRmXgwS59/Ylmej6EARhQcKKTpQhEgTdO3alZN+J5GhEmMHI8YII57tPKs956iroNS138jFI0k9lYrlpIXS46WUHC+hNLUU6wkrukIdgwYNctZ+IJGSyBJ0xspOOrWtaBhKKBTllEUjViXcHM6SuCUkJSUx+r3RlJwoofREKaUnSrFmWdH76vEM8WT8peOJiooiKCiIDw9+SL4hH71Zj86kQ+ehQ3gIwgLC+G7CdzgcDkpKShj/1XjS89KRFom92I690HkzW82MixrHsWPHOHToEJv3bEZaJB4hHhgjjJg6mDBFOm97n9hbXtG6PhaPou6oMwpFORVdrGWU7d29vLwYMGAAvUb2qiQm0i6x5ljxzfdlSNQQUlJSOHLkCMf/OO780BfZcZQ4cFgdSKvkkOUQHR7rgF6vx2QykW3LRufpFBGdlw69WY/erMditmDuaWbEiBE8/PDD/GPvP8jyyKqcSo5TxCqWva82glVtKxoVJRStnLp8yKqKidAL/Nr7EX9z5W/spMVJNVonFbcA57JiZoybUf74qfCn6nwA2dYqSjU3lFC0AWr7kNX1G/tc1sn59FOWQstBnVEo6oWWXg9F06IOMxUKRa2oa48qFApNUUKhUChqRQmFQqGoFSUUCoWiVpRQKBSKWlFCoVAoakUJhUKhqBUlFAqFolaabcCVECITOFqPl7QDTrlpOVqj1uoe1FrrRycpZUhdOjZboagvQogtdY0ya2rUWt2DWqv7UFsPhUJRK0ooFApFrbQmoZjb1AuoB2qt7kGt1U20mjMKhULhPlqTRaFQKNxEixUKIcR4IcQeIYRDCFHj6bEQYpQQ4k8hxEEhxPTGXGOFNQQJIVYKIQ64fgbW0M8uhPjDdVvayGs85/skhDAKIRa6nv9dCBHdmOurspba1nq3ECKzwns5pYnW+ZEQ4qQQYncNzwshxJuu32OnEGJAY6+xzkgpW+QN6AX0ANYAA2voowcOAV0AT2AH0LsJ1vpvYLrr/nTglRr6FTTRe1nr+wQ8Arznun8rsLAZr/VuYE5TrK/KOq4EBgC7a3j+euB7QACXAL839ZprurVYi0JKmSSl/LOWbhcDB6WUyVJKC7AAGO3+1Z3FaGCe6/48YEwTrOFc1OV9qvg7LAaGiarlsxuH5vI3rRUp5Tog+xxdRgOfSicbgQAhRHjjrK5+tFihqCMdgJQKj1NdbY1NeyllWVnqdKB9Df1MQogtQoiNQojGFJO6vE/lfaSUNiAXCG6U1dWwDhc1/U3Husz5xUKIqMZZWr1pLv+ftdKsq3ALIX4Cwqp5aoaUckljr+dcnGutFR9IKaUQoiZXUycp5XEhRBdgtRBil5TykNZrbQN8B3whpSwVQjyI0xK6ponX1KJp1kIhpRzewCGOAxW/TSJdbZpzrrUKITKEEOFSyjSXaXmyhjGOu34mCyHWABfi3I+7m7q8T2V9UoUQBsAfyGqEtVWl1rVKKSuu6wOcZ0TNkUb7/2worX3rsRm4QAjRWQjhifMQrlG9CS6WApNd9ycDZ1lDQohAIYTRdb8dcDmwt5HWV5f3qeLvMA5YLV0nco1MrWutss+/CUhqxPXVh6XAXS7vxyVAboUtavOiqU9TG3CifDPOPV0pkAH84GqPAJZXOVnej/ObeUYTrTUYWAUcAH4CglztA4EPXPcvA3bhPMXfBdzXyGs8630Cngduct03AV8CB4FNQJcm/NvXttZZwB7Xe/kz0LOJ1vkFkAZYXf+r9wEPAQ+5nhfA267fYxc1eO+aw01FZioUilpp7VsPhUKhAUooFApFrSihUCgUtaKEQqFQ1IoSCoVCUStKKBQKRa0ooVAoFLWihEKhUNTK/wP0P6fq6st1IwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Plotting\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.path import Path\n",
"import matplotlib.patches as patches\n",
"\n",
"sa = rebound.SimulationArchive(\"test.bin\")\n",
"verts, codes = sa.getBezierPaths(origin=0)\n",
"fig, ax = plt.subplots()\n",
"for j in range(sa[0].N):\n",
" path = Path(verts[:,j,:], codes)\n",
" patch = patches.PathPatch(path, facecolor='none')\n",
" ax.add_patch(patch)\n",
" ax.scatter(verts[::3,j,0],verts[::3,j,1])\n",
"ax.set_aspect('equal')\n",
"ax.autoscale_view()"
]
},
{
"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.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment