Created
June 6, 2024 13:31
-
-
Save hannorein/477d1820c949aa954fba391a41bd54ee to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "id": "9422750a-f104-4fa6-a6bc-f8eec2659781", | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import rebound\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "import numpy as np" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "id": "f4234cfd-8374-4271-967c-284cb405bc4b", | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "sim = rebound.Simulation()\n", | |
| "sim.add(m=1)\n", | |
| "sim.add(m=1e-3,a=1.0,e=0.1)\n", | |
| "sim.add(m=1e-3,a=3.1345,e=0.1)\n", | |
| "sim.init_megno()\n", | |
| "times = np.linspace(0,1e6,1000)\n", | |
| "megno = np.zeros(len(times))\n", | |
| "for i, t in enumerate(times):\n", | |
| " sim.integrate(t)\n", | |
| " megno[i] = sim.megno() " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "id": "c8c4e266-9baf-41c8-b8f3-9026727da871", | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "sim = rebound.Simulation()\n", | |
| "sim.add(m=1)\n", | |
| "sim.add(m=1e-3,a=1.0,e=0.1)\n", | |
| "sim.add(m=1e-3,a=3.1345,e=0.1)\n", | |
| "sim.integrator = \"whfast\"\n", | |
| "sim.dt = sim.particles[1].P/20.12345\n", | |
| "sim.init_megno()\n", | |
| "megno_whfast = np.zeros(len(times))\n", | |
| "for i, t in enumerate(times):\n", | |
| " sim.integrate(t)\n", | |
| " megno_whfast[i] = sim.megno() " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "id": "2703ab5b-320c-4a1b-ae7a-aa54485aa44c", | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x1198acac0>" | |
| ] | |
| }, | |
| "execution_count": 13, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAGwCAYAAABCV9SaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABG8ElEQVR4nO3deXxTZcL28SvpXqAtBdpSaAUE2XcEiwrMiCIwLuO86PCg1o15VBwEFJWZcX8cGMcFdFTEjXFFHRFHXBDRoigurIogggIVaFkEWkqhS3K/f5wmNFCwhBNy2vy+n09Mc87JyX3Hkly9t+MyxhgBAADUc+5wFwAAAOBEIPQAAICIQOgBAAARgdADAAAiAqEHAABEBEIPAACICIQeAAAQEaLDXYATzev1auvWrWrUqJFcLle4iwMAAGrBGKO9e/cqMzNTbndwbTYRF3q2bt2qrKyscBcDAAAE4eeff1bLli2Dem7EhZ5GjRpJst60pKSkMJcGAADURnFxsbKysvzf48GIuNDj69JKSkoi9AAAUMccz9AUBjIDAICIQOgBAAARgdADAAAiAqEHAABEBEIPAACICIQeAAAQEQg9AAAgIhB6AABARCD0AACAiEDoAQAAEYHQAwAAIgKhBwAARISIu+AoguT1SqU7JZdbik+RovjVAQDULXxz4cgq9ksrXpa+e1P6+UvJU25tj4qTMrpKpwyROl0oNTslrMUEAKA2XMYYE+5CnEjFxcVKTk5WUVGRkpKSwl0cZzJGWvWGNO8vUsm2ajtckmr4dWkzSMr5s9T2LMnlOkGFBABEEju+v2npQaADRdKc66Xv51qPk7OlvqOl9kOlxq2tULN7o7Tpc2nNf6X1C6Sf8qxbiz7SOf8nnZQTxgoAAFAzWnpw0O6N0st/lHaskdwx0oCbpTMmSNGxR3nOJunLJ6WlM6WKfda29sOls++RmrY9EaWuHU+FJJfkjqI1CgDqIDu+vwk9sPzyo/TcMKmkUGrUXPrjy1KLXrV//t5CKW+ytOx5yXgld7R06jXSwFulxNTQlftIZVn7npT/hVSwQiraIpXvtfZFx0vJLaXGraTMXlJWX+mk/lJsgxNbRgDAMSH0BIHQU4NdG6SZw6XiLVJaJ+nSN6SkzODOtf17af4d0rp51uP4ZGnARKnvn6ToOPvKfKiKA9J3s6XlL1pdbzWNPTqSqDhrXFKH4VKn86WExqEqJQAgSISeIBB6DlG6S3p6sLTrR6lpe+mKd6SGzY7/vD9+LH3wN2nbKutx41bSWXdas73cNi4Ptednackz0tJ/S/t3Hdzeoo81sLpFHym1tdSgqTVA+8AeqWiztPMHafMSadNn0p78g8+LirXGL3UfKbUdLEXF2FdWAEDQCD1BIPRU46mQXrxI2vCJlJwlXfOh1CjDvvN7PdaU94/uPTgLrOkp0uk3Sl0vPvpYoaMxRtr4qfTVDOn7d6zuNElKaimdepV17pSs2p9r+xpp7TvSqtnS9tUH9yU2lbr+P6nbxVZX2IkcC1RZbo2x2vWjtH+3VLbXej9jEqS4RlYXXcpJ1v8vxigBiACEniAQeqoYI80dLy19ToptKF01T8roEprXKiuRPn9U+uIJqazI2uYLFF1HSJk9rQHGv2b3RmnN21YX1o7vD25vPUDq+7/SKece36KJxkiF30orZ0nfvibt23FwX5O2UrdLrPKmtg7+NY7EUyFtWWq1kP2UJ21ZInkrf/15CalSZg+peQ8pO8eaORfXyP7yAUCYEXqCQOip8uWT0nu3SHJJI1+xunRC7UCxNcvri8elvQUHt8enSK3OsMYTNWkrJaRYIah0t9VCtG2VtHV5YNCJaSD1GGkNlk7raH9ZPZXSjx9Z4WfNXKly/8F9LfpYAavd2VJGt+C667xea5bchk+skLNxkVReEnhMbEMptY3UMM0KMi63NXbpwB6p6GdrgLbxBD7HHS216G0FwdYDpJZ9pZj4Yy8fADgMoScIhB5ZX+Yv/sHqFjr7Xun0sSf29T2V0k8fW11f6z+Uyopr9zyXWzrpdKnTBVaXU3xyaMvpU7bX6kb75lUroPi60ySpQTMrWLToZa1SndxSSmphlc3lslqPykukku3SjrVWcNu8RMr/3Oq2qi4h1RpQ3WaQFVgatzp611VlmbTtO2uG2palVnDavTHwmKg4Kbuf1HqgdcvsGbpLiFSWS7+ss8ZLFRdYMwFLd1mtWJ4ya7xUXCPrvUnKlFJaWXVs3IrLmgD4VYSeIER86PnlR+mp31qtBd1HShc+Ed4xIZ5KqxXn5y+kneukXT9ZIcFbac2iSmxiDbDO7GGFiwZNwldWqWo6/LvSug+tAORbm6gmUbHWOJxDW2N8YhpYgcQXdNK7Hv8g792brPFOPy20WpFKCgP3xzaypuiflGO1UjXvbg3yPhbGWO/DtlVWd+C276yxUDt/qF2X3KGi4qT0TlZ5Mrpa46cyuoR2th+AOofQE4SIDj0HiqSnz5Z2rpVanirlzqXr43hUlklblklbl1ktLTt/sGaGHdqCI0nRCdZijc06SOldrO685t1DOzvMGCtIbqgKQBs/rblsjZpb3Wi+VqqEFCkm0QptlQek8n3WzLiizdZt109S6S81v2ZcktXdmNTCGmSdkGqFl6hY69ptZXutwL3nZ2nPJqtlqqL08PO4Y6T0zlZXXYteVhBq1r52Y78A1EuEniBEbOjxeqRXRlrr5zTKlP70sb0ztXBQ+T7rVllmfUnHp1izrsI9y8rrtVpnNiy0wlrBSmt2WDBcbmv8VXoXq1UmvYs1Jiu55bHV0+uVdm+QCr+xWo0KVlotfzWFqpgGVotfZs+DYSjlpON/Xw8UV42R2mwtX1Cy3QpiFfutoBYdb/3/i0m0xlf5wmFKthQfQZ8hQJgReoIQsaHnw7ukRQ9bH+BXvndsqy2j/irba03Z35NvfekXb7G2VZRaoS0mwQob8cnWl31yS2s5gGYdrH2hYIxVni1Lq1rRllvjlg4d6C1Z3Z9pnawlF3zli0+yyhybaI2/qiyzWqxKd1ndcnsLrFvRZqvFyTejMBhJLa2uubSOUlpnKwA2PYX1nYAQIPQEISJDzzevS7OvsX6+6Gmp24jwlgc4Vl6P1X3o705cZrUMeSvsOX9C46rQlGV198U2qOrii7FaeypKraUX9hZKxZutmXPVF8OsLirOCkEZXa0uzIyuVktYXEN7ygpEKEJPECIu9GxaLD1/gTV75ozx0uC7wl0iwB6VZVZ33S8/Wt1Te36uaqkqsQaYl5da3YvR8da4ooTGUsN0K9Q0Sq9qHcqSklsEt7bR/t3WZVe2fydtW20N5i5cdfA6bwFc1rip5lWDtTO6WbdG6cf9NsjrscpSsl3at13at9PqHvRWWq1mxmu1esUlW61gCalWnRumM0YKdQqhJwgRFXp++VF6+izrA7HD76SLn+dDDgglr1fas9FqhSr8ViqoGqu0d2vNx8cnH+ya8w0ij2tk3aLjrcDi9VitTaW7rNal0l1S6U4r3JRst36uvoxCbbmjraUDkrOlpu2sbrmmp0jNTrG67ey8XAxgA0JPECIm9Oz7RXpmsDXTJrOXdU2t2MRwlwqITCU7pG2HBKFf1gUXVo4kIdVaN6phmtWqFRVrDTiXrO65A0XWmlilu6TirUdeSkE6ONuwSTtrwLr/drIVzIAwsOP7mxXB6qOyEumVP1qBJyVb+p9XCTxAODVsJjX8rXTybw9uKy89OIC8KN8KIgeKrYHkZcUHZ/+5oqzFGxNSrYHbiU2kxKqA4ws5iU2ObfC012ONTyrabC0bsPOHqts66Zf11grkvtaqQyU2tQJQSpbVRdYwTWqYUbXUQdUst+h4K3B5yq1b5QErdO3fXdVitbvqtkvav8eqc3mJ9dlVXmKVzx1VdYuxWr4SUqyZkAmNrRaqpBZV47BaWC1TDZrROoVfRUtPfVNeKr18sbUmS3yydPV8a30TAKgNT6W1htKOtVZr1C8/Vt3WH77YpZNExVXNLsyuumVZSxokZ1mPG2WErnu/svxg1+P+XdZyBy631YXojrLKltDYCqvxyQwzCBItPQhUcUB6dZQVeGIbSZfOJvAAODZR0VY3VpOTD99XttdqQf7lR2vQeMk2ae82KwwdKLa+7Cv3W398SVYXW3Ssde9rpUlMte6r3+IaWdeai2to3bujqlYz91qXMTlQZC1qeaCoqntuy8ElFoo2W61WnjJr3akjrT3ljrFahVKyrXFMSc2tABKXZA3wjkuyQoqPt6Kq1a3q5nvt6uGmtKq1qqblFI4mPtlqHfO1WCVlHv5zQuPwr+1VD9HSU19UlkmvXS798L7VvHzpbOtSAwBQ33kqrAC0J79qte/8qhl9+VarVU0X57Wd62CIi21gzZzzVlqv67tQcG2vMyhZXYSNMqzFZJOaW7MOkzKrVjpvfHA2XlySFRrd0YEhqWL/wYVSy6sFtwNFVpei7+cDewK3l+2VZKz6uNzWOWMbHBxg73u9Bk0PdrFW72pNaByylixaemAp2yvN+h/rUgPR8dLIWQQeAJEjKubgxWtr4qmsWpDSF4TyrVaqA8VWEPGNpaoejFxRVRfITTr4hZ+QWtVSdeh9Y6sl69fGFPlarUp/qVrzaasV1oq3Vt02W/elv1jjoHZvPPwiwk7nclvjvrL6Sn98KdylOQyhp64r3WVdMX3rMqtZ+I8vSW0GhrtUAOAcUdFVY3yyrAvuhq0cMVUtJE2PPvSgsswKacUF1nIHxVWriBdvtcLSgT1VLTdVoU1H6LBxRVnfC7ENrK7D+BSray2h6j4++fBtcY2s5xmvdV6v11r3yv96vq6+ndK+HdbMxH07rDWi9u+2nrdvu3WMAxF67LR1hfTjAqn/jdY/slDbud6apfXLOusvjkv/Y12TCABQd0XHHb3lqjpfKPH6WqmM1bUWk2id50SOC/JUWOtH7dvh2PFIYQ09kydP1uzZs/X9998rISFB/fv31z/+8Q+1b3/0wbevv/66br/9dm3cuFHt2rXTP/7xDw0bNuwElfoo/n2elbrL9oZ+5eP1C6T/XGml6aSW0mX1Y9Cyx2tU4fGq3ONVeaVXFR6vKiqNKrxeebzGf/Mao0qvkbfaNo+pts9j3Xu8ksdYx9nJ5ZLcLpei3C7/fZTb2nak7b5tR9zudinK5ZLbLUW5XIqJdis2yq2YKLei3M78AAEQZm53cCuKh0JUjDX+KKl5uEtyRGENPQsXLtSYMWN06qmnqrKyUn/5y190zjnnaPXq1WrQoEGNz/n88881cuRITZ48Wb/73e/08ssv68ILL9SyZcvUpUuXE1yDQ/gGqS16WBp4mxQTb/9reCqkvCnSpw9KMlLLqn7ThmlBn9IYo7JKb9XNo7IK6+dy32PfvgpPjcf5j6k4/Pjq4aXcYw4GmerbK62QU+GxQgsO53JJMVG+EORSdLWfY6qCUcDP0W7FRrkU7bZ+jolyKTbKreiqY3xhKjb64H1s1XP823zbazzukG1Vr+9y6F93ACA5bPbWjh07lJaWpoULF2rAgAE1HnPJJZdo3759mjt3rn/baaedph49emj69Om/+hohnb01Oetg8GnSVup6sdTlImuJ90P4g0aFVwcqPTpQ4dGBCm/VvUcHKq2ffeGjwuNVgz1rddqqu5VWbC0Y9k3aBXo/+ybtN9H+FpHyaq0k5dWCy+Hbqu6rtjuV78s6ym190VstI1ZLSFSUr2XEurdaTw5pZfG1pLgll379C7m239leX6tSVUuSr4XJ17rkrWp58lZrgTJG/hap6tv9x5qDz3POv8pj4ws/vnAUc0hwqh62AsNT1XOiohQT7VLcIcfVeJ6qe+v/tw62lPlayw5tWXO55HLJ/zvicsn/+xHwu1H9x5o3B4S7wO3Vjz/4wPe74TWyhklUe2xk/f/2PzaHP65+7622338v32PfsVXP9x58jveQ81Y/PvC1rN/rYzq+enm8NR/v+ZX9Rzzfoa/vrc3xh79nNdXPW+0fmksu//8/l8v6v+dyVd2q9rmq/U/27z/0eFkbAx7XcC7f60RVtRq7XNVai32fYa7qj1X1e+v7Pa96frXff3e1fe5qn4O+3/Xq5wo47yH7Dn0dt//f0tFeR9We41KDuCi1aWbvRXbr3eytoiJr4FNqauoRj1m8eLEmTJgQsG3IkCGaM2dOjceXlZWprKzM/7i4+BimDB4j06yDXJu/sh78sl7K+7uU93etd7fRh1Fn6D3laJOnqT/M1PaLLUn7NC76DY2I+kDRLq+KTaImVVyjd/JPk/I321oHl0uKi3YrLjrKuo+p9rNve4z7KMdY+2Oj3P77Q1sGAloHol3+L8K46OqtFJHdcuD1Wt13/laxqpawSn8rmbWv0nvw54O36s8zqqi0jqs4pKWtwmMF75pa3Q62vh0MyoceX+7xHtYyZ4VuaV95qKcHA3CyntkpevP608NdjMM4JvR4vV6NGzdOp59++lG7qQoLC5WeHnhl4vT0dBUW1rxS6OTJk3X33XfbWtYj2bK7VC0l3VR+rVwuo2HuL3Wm+1u19f6ktt6fdK2e11JvO/3X21/vmn7aoRRJVqqOj3YrPiZK8TG+IBGl5u7dOq/sHZ1d+o4aeq0rN69sOEBvN/+zkhOb64pqf1X7w4R/W1TAX9NxMW7FVf8LuyqgxPrDjPU4koOGk7jdLsW6rf+3TuYbg1VTeArcZlTu8VSFJRMYsqoFrYpDAleZ77mVnoPnOSSYeata2TxVf8EHtroFtqD5W9yMc7pS3f6/9K2/lN1VrQE1bfe1CFhjyqoeq2q/+5DH1cadVX/sPuR8bv/5qp2jtscf2gpwxOcePDbgub92/KEtEzXsr13Zazi+aptcssb+quqi9FUtb9ZjqxVN1bcfss9UHeDbriOc57DtVft8LVieaq1P1u+rqrZXa62q+tnjf87h+2r63Q9oaTuk5bl6y5jv8ZH2mWr/zqrvO/TfoNcYNWkQdyL/GdWaY0LPmDFjtGrVKi1atMjW806aNCmgZai4uFhZWVm2vobPvvJKSVKJq4HS+16kHZk36NOofcre9qHSNr2jRoVfqLd7nXq71+mu2BfkyT5drvZDFZXR2bp2TnmptY7EjrXST3nSlqUH141oeoo09B/qfvJv1T0kpQeOndVVZIX1usgXhqqr/tBUmwp8tJbZIz3Ht883EP3QL1/+wABOLEeEnhtuuEFz587VJ598opYtWx712IyMDG3bti1g27Zt25SRkVHj8XFxcYqLO7GJ838HnaxeZ1dvreoo6c/WWgur50jf/keuLUsUvelTadOnRz9Z1mlS/xuk9sO4XgtgM7fbJXctxnoBqB/CGnqMMfrzn/+sN998U3l5eWrduvWvPicnJ0cLFizQuHHj/Nvmz5+vnBwHrED8a63lSc2l066zbrs3St+9KW1eIu343prm7rtgXuOTrAW02vzGWkwLAAAct7CGnjFjxujll1/WW2+9pUaNGvnH5SQnJyshIUGSdPnll6tFixaaPHmyJOnGG2/UwIED9eCDD2r48OGaNWuWlixZohkzZoStHj6uqtRTq78bG7eSzhgfyuIAAIBqwjpK8oknnlBRUZEGDRqk5s2b+2+vvvqq/5j8/HwVFBT4H/fv318vv/yyZsyYoe7du+s///mP5syZE/41egLQXA4AgNOEvXvr1+Tl5R22bcSIERoxYkQISnR8nDEXBAAA1MTZ82HrGF/3llOvOQIAQCQj9IQEoQcAAKch9IQAkQcAAOch9IQAY3sAAHAeQg8AAIgIhB4b+dfpoX8LAADHIfSEBKkHAACnIfQAAICIQOixEd1bAAA4F6HHRsZ/T+oBAMBpCD02chnfBUcJPQAAOA2hJxTo3wIAwHEIPQAAICIQemxE+w4AAM5F6AkFurcAAHAcQo+NuOYWAADOReixkYvYAwCAYxF6AABARCD02OjgisyM6QEAwGkIPQAAICIQemx0cEQPLT0AADgNocdGDGQGAMC5CD2hwJgeAAAch9ADAAAiAqHHRq5D7gEAgHMQeuzkG9JD9xYAAI5D6AEAABGB0GMrZm8BAOBUhJ6QoHsLAACnIfTY6OBlKMJcEAAAcBhCT0iQegAAcBpCDwAAiAiEHhtxGQoAAJyL0GMjLjgKAIBzEXpCgcwDAIDjEHpsRNYBAMC5CD2hwJx1AAAch9BjI/86PbT5AADgOIQeGxn/PaEHAACnIfSEAL1bAAA4D6HHRqzTAwCAcxF6AABARCD0hAC9WwAAOA+hx0YuQ/cWAABORegJAWZvAQDgPISeEHAxfQsAAMch9IQCoQcAAMch9IQAkQcAAOch9NiIdXoAAHAuQk8o0L0FAIDjEHoAAEBEIPTYiO4tAACci9ATEnRvAQDgNIQeW1ktPQzpAQDAeQg9IUHqAQDAaQg9AAAgIhB6bET7DgAAzkXoCQHDoB4AAByH0BMCRB4AAJyH0GMj1ukBAMC5CD028kceurcAAHAcQo+NfC09RB4AAJyH0BMSxB4AAJyG0BMKdG8BAOA4hB4bEXUAAHAuQo+dDo5kDmcpAABADcIaej755BOdd955yszMlMvl0pw5c456fF5enlwu12G3wsLCE1PgWiLyAADgPGENPfv27VP37t312GOPHdPz1q5dq4KCAv8tLS0tRCU8VqzTAwCAU0WH88WHDh2qoUOHHvPz0tLSlJKSUqtjy8rKVFZW5n9cXFx8zK93zBjIDACA49TJMT09evRQ8+bNdfbZZ+uzzz476rGTJ09WcnKy/5aVlRXy8pF5AABwnjoVepo3b67p06frjTfe0BtvvKGsrCwNGjRIy5YtO+JzJk2apKKiIv/t559/Dln5uAwFAADOFdburWPVvn17tW/f3v+4f//++vHHH/Xwww/rhRdeqPE5cXFxiouLO1FFrEJTDwAATlOnWnpq0rdvX61fvz7cxZBESw8AAE5W50PPihUr1Lx583AX4xC09AAA4DRh7d4qKSkJaKXZsGGDVqxYodTUVGVnZ2vSpEnasmWLnn/+eUnS1KlT1bp1a3Xu3FkHDhzQ008/rY8++kgffPBBuKpQIwYyAwDgPGENPUuWLNFvfvMb/+MJEyZIknJzczVz5kwVFBQoPz/fv7+8vFw33XSTtmzZosTERHXr1k0ffvhhwDnCiawDAIBzuYwxETUQpbi4WMnJySoqKlJSUpKt595+V2ulaZfWX/Se2nbrb+u5AQCIZHZ8f9f5MT0AAAC1QeixEbO3AABwLkJPKDCSGQAAxyH02IiWHgAAnIvQExK09AAA4DSEnhBwEXoAAHAcQo+NiDoAADgXoScEGMcMAIDzEHoAAEBEIPTYitlbAAA4FaEnFNz0bwEA4DSEHhuxTg8AAM5F6AkJWnoAAHAaQk8IsE4PAADOQ+ixEVEHAADnIvSEgIuFegAAcBxCDwAAiAiEHlsxewsAAKci9IQC3VsAADgOocdGrNMDAIBzEXpCgpYeAACchtATAvRuAQDgPIQeG9G9BQCAcxF6QoKmHgAAnIbQAwAAIgKhx0a07wAA4FyEnhBwuYk/AAA4DaHHVgxkBgDAqQg9IUFLDwAATkPoAQAAEYHQYyPW6QEAwLkIPSHgcvG2AgDgNHw7AwCAiEDosZGL3i0AAByL0BMCLq44CgCA40Qfz5N37NihtWvXSpLat2+vZs2a2VKouoqBzAAAOFdQLT379u3TVVddpczMTA0YMEADBgxQZmamrr76apWWltpdxjqIlh4AAJwmqNAzYcIELVy4UP/973+1Z88e7dmzR2+99ZYWLlyom266ye4y1jn0bgEA4DxBdW+98cYb+s9//qNBgwb5tw0bNkwJCQm6+OKL9cQTT9hVvjqF7i0AAJwrqJae0tJSpaenH7Y9LS2N7i1JjA8HAMB5gvp2zsnJ0Z133qkDBw74t+3fv1933323cnJybCscAACAXYLq3po2bZqGDBmili1bqnv37pKklStXKj4+XvPmzbO1gAAAAHYIKvR06dJF69at00svvaTvv/9ekjRy5EiNGjVKCQkJthawLnK5GckMAIDTBL1OT2JiokaPHm1nWQAAAEIm6NCzbt06ffzxx9q+fbu8Xm/AvjvuuOO4C1YXMXsLAADnCir0PPXUU7ruuuvUtGlTZWRkBFx2weVyRWzoAQAAzhVU6Pm///s/3Xfffbr11lvtLk+dRksPAADOFdSU9d27d2vEiBF2l6Xe4IKjAAA4T1ChZ8SIEfrggw/sLgsAAEDIBNW91bZtW91+++364osv1LVrV8XExATsHzt2rC2Fq2to3wEAwLlcxphjHojSunXrI5/Q5dJPP/10XIUKpeLiYiUnJ6uoqEhJSUm2nrvszqaKc1Wo8Oqlyshqa+u5AQCIZHZ8fwfV0rNhw4agXgwAACBcuDKmjZi9BQCAcwXV0jNhwoQat7tcLsXHx6tt27a64IILlJqaelyFq7sY3QMAgNMEFXqWL1+uZcuWyePxqH379pKkH374QVFRUerQoYMef/xx3XTTTVq0aJE6depka4GdjZYeAACcKqjurQsuuECDBw/W1q1btXTpUi1dulSbN2/W2WefrZEjR2rLli0aMGCAxo8fb3d56wQuOAoAgPMENXurRYsWmj9//mGtON99953OOeccbdmyRcuWLdM555yjnTt32lZYO4Ry9lb5namKdXm0bfRypbdoY+u5AQCIZHZ8fwfV0lNUVKTt27cftn3Hjh0qLi6WJKWkpKi8vDyoQtVVtO8AAOBcQXdvXXXVVXrzzTe1efNmbd68WW+++aauvvpqXXjhhZKkr776SqeccoqdZa0zXEyKAwDAcYIayPzkk09q/Pjx+uMf/6jKykrrRNHRys3N1cMPPyxJ6tChg55++mn7SgoAAHAcghrT41NSUuJffblNmzZq2LChbQULlVCO6am8s7GiXV7t+NNKNctsZeu5AQCIZGEb0+NTWFiogoICtWvXTg0bNtRx5Kd6htE9AAA4TVCh55dfftFZZ52lU045RcOGDVNBQYEk6eqrr9ZNN91kawHrElZkBgDAuYIKPePHj1dMTIzy8/OVmJjo337JJZfo/ffft61wdZXLRUsPAABOE9RA5g8++EDz5s1Ty5YtA7a3a9dOmzZtsqVgAAAAdgqqpWffvn0BLTw+u3btUlxc3HEXqq6ifQcAAOcKKvSceeaZev755/2PXS6XvF6v7r//fv3mN7+p9Xk++eQTnXfeecrMzJTL5dKcOXN+9Tl5eXnq1auX4uLi1LZtW82cOTOIGoSWIf0AAOA4QXVv3X///TrrrLO0ZMkSlZeX65ZbbtF3332nXbt26bPPPqv1efbt26fu3bvrqquu0kUXXfSrx2/YsEHDhw/Xtddeq5deekkLFizQNddco+bNm2vIkCHBVAUAAESIoEJPly5dtHbtWj322GNq1KiRSkpKdNFFF2nMmDFq3rx5rc8zdOhQDR06tNbHT58+Xa1bt9aDDz4oSerYsaMWLVqkhx9++Iihp6ysTGVlZf7HvstkhILbZc3ectHRBQCA4wQVeiQpPj5eZ599trp37y6v1ytJ+vrrryVJ559/vj2lO8TixYs1ePDggG1DhgzRuHHjjvicyZMn6+677w5JeY6I2VsAADhOUKHn/fff12WXXaZdu3YdtiChy+WSx+OxpXCHKiwsVHp6esC29PR0FRcXa//+/UpISDjsOZMmTdKECRP8j4uLi5WVlRWS8gEAAOcKaiDzn//8Z1188cXaunWrvF5vwC1UgSdYcXFxSkpKCriFGuv0AADgPEGFnm3btmnChAmHtbqEWkZGhrZt23ZYWZKSkmps5QEAAPAJKvT8v//3/5SXl2dzUX5dTk6OFixYELBt/vz5ysnJOeFlOQzXHQMAwNGCGtPzr3/9SyNGjNCnn36qrl27KiYmJmD/2LFja3WekpISrV+/3v94w4YNWrFihVJTU5Wdna1JkyZpy5Yt/jWBrr32Wv3rX//SLbfcoquuukofffSRXnvtNb3zzjvBVCNkmL0FAIDzBBV6XnnlFX3wwQeKj49XXl5ewBgWl8tV69CzZMmSgMUMfQOOc3NzNXPmTBUUFCg/P9+/v3Xr1nrnnXc0fvx4TZs2TS1bttTTTz/tiDV6jDFEHQAAHMxlDp1+VQsZGRkaO3asbrvtNrndQfWQhU1xcbGSk5NVVFRk66Bm4/XIdU+qJGn3mO/VuFnt1ysCAABHZ8f3d1CJpby8XJdcckmdCzwnDLO3AABwnKBSS25url599VW7y1KnGS8DmQEAcLKgxvR4PB7df//9mjdvnrp163bYQOaHHnrIlsLVVQxkBgDAeYIKPd9++6169uwpSVq1alXAPhbmAwAAThRU6Pn444/tLkedZ0T3FgAATsZI5BBwuWntAgDAaQg9Ngli5j8AADiBCD12CQg9tPQAAOA0hJ6QIPQAAOA0hB6b0LkFAICzEXpsUy320NADAIDjEHpsEjikh9QDAIDTEHrsUi31EHkAAHAeQk8o0NIDAIDjEHpsworMAAA4G6HHNnRvAQDgZIQeAAAQEQg9NuEqFAAAOBuhxy7VZ28xkBkAAMch9IQCoQcAAMch9NiGgcwAADgZoSckiD0AADgNoccmDGQGAMDZCD12CRjIHMZyAACAGhF6QoHUAwCA4xB6bMJlKAAAcDZCj22qz96ipQcAAKch9NgkcCAzoQcAAKch9NiFgcwAADgaoScUSD0AADgOoccmDGMGAMDZCD22IfYAAOBkhJ6QoHsLAACnIfTYhMtQAADgbIQeuzB7CwAARyP02CRwmR5SDwAATkPosQ0rMgMA4GSEnpAg9AAA4DSEHpswkBkAAGcj9NiGgcwAADgZoScUSD0AADgOoccudG8BAOBohB6buKqv0xPGcgAAgJoRemxiApp6iD0AADgNocc21QcyE3oAAHAaQk8oEHoAAHAcQo9NWKcHAABnI/TYhYHMAAA4GqEnJIg9AAA4DaHHJvRuAQDgbIQe23AZCgAAnIzQYxNTfSQzqQcAAMch9NjFsE4PAABORugBAAARgdBjEwYyAwDgbIQe2xB7AABwMkKPTXxDeryG8TwAADgRoQcAAEQEQo9t6N4CAMDJCD02I/oAAOBMhB4AABARCD12MbTxAADgZIQeuxjfHbO3AABwIkKPbWjpAQDAyQg9NjH+e1p6AABwIkeEnscee0ytWrVSfHy8+vXrp6+++uqIx86cOVMulyvgFh8ffwJLCwAA6qKwh55XX31VEyZM0J133qlly5ape/fuGjJkiLZv337E5yQlJamgoMB/27Rp0wks8REwkBkAAEcLe+h56KGHNHr0aF155ZXq1KmTpk+frsTERD377LNHfI7L5VJGRob/lp6efgJLXDNzyD0AAHCWsIae8vJyLV26VIMHD/Zvc7vdGjx4sBYvXnzE55WUlOikk05SVlaWLrjgAn333XdHPLasrEzFxcUBNwAAEHnCGnp27twpj8dzWEtNenq6CgsLa3xO+/bt9eyzz+qtt97Siy++KK/Xq/79+2vz5s01Hj958mQlJyf7b1lZWbbXQxLdWwAAOFzYu7eOVU5Oji6//HL16NFDAwcO1OzZs9WsWTM9+eSTNR4/adIkFRUV+W8///xzSMpl/B1bzN4CAMCJosP54k2bNlVUVJS2bdsWsH3btm3KyMio1TliYmLUs2dPrV+/vsb9cXFxiouLO+6y/jpaegAAcLKwtvTExsaqd+/eWrBggX+b1+vVggULlJOTU6tzeDweffvtt2revHmoilkrhhWZAQBwtLC29EjShAkTlJubqz59+qhv376aOnWq9u3bpyuvvFKSdPnll6tFixaaPHmyJOmee+7RaaedprZt22rPnj365z//qU2bNumaa64JZzUAAIDDhT30XHLJJdqxY4fuuOMOFRYWqkePHnr//ff9g5vz8/Pldh9skNq9e7dGjx6twsJCNW7cWL1799bnn3+uTp06hasKkiQXA5kBAHA0lzGR9W1dXFys5ORkFRUVKSkpybbzFuavU8azfXTAxCj+7p22nRcAANjz/V3nZm8BAAAEg9Bjm4hqMAMAoM4h9NjE10vI7C0AAJyJ0AMAACICoccmrNMDAICzEXoAAEBEIPTYJbJm/gMAUOcQemxG9AEAwJkIPTZjRA8AAM5E6LENbTwAADgZoccmzN4CAMDZCD02cdHSAwCAoxF6bGLEiswAADgZoQcAAEQEQo9dWKcHAABHI/TY5OBAZgAA4ESEHgAAEBEIPbahjQcAACcj9NjE+O+ZvQUAgBMReuziG9RD5gEAwJEIPbbxtfWQegAAcCJCDwAAiAiEHruwTg8AAI5G6LGJOeQeAAA4C6EHAABEBEKPTejdAgDA2Qg9NmOdHgAAnInQYxeaegAAcDRCj21M1X9p6QEAwIkIPQAAICIQemxD9xYAAE5G6LELmQcAAEcj9NiOMT0AADhRdLgLUF8YmnoAoF7xeDyqqKgIdzEiSmxsrNzu0LXHEHpsRvQBgLrNGKPCwkLt2bMn3EWJOG63W61bt1ZsbGxIzk/oAQCgGl/gSUtLU2Jiolwuhi2cCF6vV1u3blVBQYGys7ND8r4TemxiWJwQAOo8j8fjDzxNmjQJd3EiTrNmzbR161ZVVlYqJibG9vMzkNlmLE4IAHWXbwxPYmJimEsSmXzdWh6PJyTnJ/TYhZYeAKg36NIKj1C/74Qe23AZCgAAnIzQAwBAPTBo0CCNGzcu3MVwNEKPTejdAgCE0+zZs3Xvvffacq777rtP/fv3V2JiolJSUmo8xuVyHXabNWuWLa8fKszeAgCgHkhNTbXtXOXl5RoxYoRycnL0zDPPHPG45557Tueee67/8ZECklPQ0gMAwFEYY1RaXnnCb8e6FEr17q0XXnhBffr0UaNGjZSRkaH/+Z//0fbt2/3H7t69W6NGjVKzZs2UkJCgdu3a6bnnnvPvv/vuuzV+/Hh17dr1qK+ZkpKijIwM/y0+Pv6Yynyi0dJjE9bpAYD6aX+FR53umHfCX3f1PUOUGBvc13RFRYXuvfdetW/fXtu3b9eECRN0xRVX6N1335Uk3X777Vq9erXee+89NW3aVOvXr9f+/fuP+XXGjBmja665Rm3atNG1116rK6+80tEz3wg9NmP2FgAg3K666ir/z23atNEjjzyiU089VSUlJWrYsKHy8/PVs2dP9enTR5LUqlWrY36Ne+65R7/97W+VmJioDz74QNdff71KSko0duxYu6phO0KPXWjpAYB6KSEmSqvvGRKW1w3W0qVLddddd2nlypXavXu3vF6vJCk/P1+dOnXSddddpz/84Q9atmyZzjnnHF144YXq37//Mb3G7bff7v+5Z8+e2rdvn/75z386OvQwpsc2rNMDAPWRy+VSYmz0Cb8F2020b98+DRkyRElJSXrppZf09ddf680335RkDVCWpKFDh2rTpk0aP368tm7dqrPOOks333zzcb1P/fr10+bNm1VWVnZc5wklQg8AAPXI999/r19++UVTpkzRmWeeqQ4dOgQMYvZp1qyZcnNz9eKLL2rq1KmaMWPGcb3uihUr1LhxY8XFxR3XeUKJ7i2bMJAZAOAE2dnZio2N1aOPPqprr71Wq1atOmz9njvuuEO9e/dW586dVVZWprlz56pjx47+/fn5+dq1a5fy8/Pl8Xi0YsUKSVLbtm3VsGFDvf3229q2bZtOO+00xcfHa/78+fr73/9+3K1FoUboAQCgHmnWrJlmzpypv/zlL3rkkUfUq1cvPfDAAzr//PP9x8TGxmrSpEnauHGjEhISdOaZZwYsLHjHHXfo3//+t/9xz549JUkff/yxBg0apJiYGD322GMaP368jDFq27atHnroIY0ePfrEVTQILhNhTRTFxcVKTk5WUVGRkpKSbDvv+pWfqe2bw7RdqUq7a4Nt5wUAnDgHDhzQhg0b1Lp1a8evOVMfHe39t+P7mzE9tomo7AgAQJ1D6LEZs7cAAHAmQo9dIquXEACAOofQYxPjv6elBwAAJyL0AACAiEDosQvdWwAAOBqhx3Z0bwEA4ESEHgAAEBEIPTaJsDUeAQCocwg9NiP6AACcbObMmUpJSfnV4+666y6lp6fL5XJpzpw5IS/XiUDosQ1xBwBQP6xZs0Z33323nnzySRUUFGjo0KHHfc5WrVpp6tSpx1+448AFR21mXAxkBgDUbT/++KMk6YILLpCrHn2v0dIDAMDRGCOV7zvxt2MYKzp37lylpKTI4/FIklasWCGXy6XbbrvNf8w111yjSy+91P943rx56tixoxo2bKhzzz1XBQUFkqxurfPOO0+S5Ha7/aHn66+/1tlnn62mTZsqOTlZAwcO1LJly6q9TUZ33XWXsrOzFRcXp8zMTI0dO1aSNGjQIG3atEnjx4+Xy+UKW5CipccmxusNdxEAAKFQUSr9PfPEv+5ftkqxDWp16Jlnnqm9e/dq+fLl6tOnjxYuXKimTZsqLy/Pf8zChQt16623SpJKS0v1wAMP6IUXXpDb7dall16qm2++WS+99JJuvvlmtWrVSldeeaU/CEnS3r17lZubq0cffVTGGD344IMaNmyY1q1bp0aNGumNN97Qww8/rFmzZqlz584qLCzUypUrJUmzZ89W9+7d9ac//UmjR4+27z06Ro5o6XnsscfUqlUrxcfHq1+/fvrqq6+Oevzrr7+uDh06KD4+Xl27dtW77757gkpaG/WnGRAAUDckJyerR48e/pCTl5en8ePHa/ny5SopKdGWLVu0fv16DRw4UJJUUVGh6dOnq0+fPurVq5duuOEGLViwQJLUsGFD/0DnjIwMZWRkSJJ++9vf6tJLL1WHDh3UsWNHzZgxQ6WlpVq4cKEkKT8/XxkZGRo8eLCys7PVt29ff8BJTU1VVFSUGjVqFHDOEy3sLT2vvvqqJkyYoOnTp6tfv36aOnWqhgwZorVr1yotLe2w4z///HONHDlSkydP1u9+9zu9/PLLuvDCC7Vs2TJ16dIlDDUAANRrMYlWq0s4XvcYDBw4UHl5ebrpppv06aefavLkyXrttde0aNEi7dq1S5mZmWrXrp0+++wzJSYm6uSTT/Y/t3nz5tq+fftRz79t2zb97W9/U15enrZv3y6Px6PS0lLl5+dLkkaMGKGpU6eqTZs2OvfcczVs2DCdd955io4Oe9TwC3tLz0MPPaTRo0fryiuvVKdOnTR9+nQlJibq2WefrfH4adOm6dxzz9XEiRPVsWNH3XvvverVq5f+9a9/neCSBzLM3gKA+snlsrqZTvTtGMe9DBo0SIsWLdLKlSsVExOjDh06aNCgQcrLy9PChQv9rTySFBMTc0gVXb+63lxubq5WrFihadOm6fPPP9eKFSvUpEkTlZeXS5KysrK0du1aPf7440pISND111+vAQMGqKKi4pjqEUphDT3l5eVaunSpBg8e7N/mdrs1ePBgLV68uMbnLF68OOB4SRoyZMgRjy8rK1NxcXHALZSIPgCAcPCN63n44Yf9AccXevLy8jRo0KDjOv9nn32msWPHatiwYercubPi4uK0c+fOgGMSEhJ03nnn6ZFHHlFeXp4WL16sb7/9VpIUGxvrH2gdLmENPTt37pTH41F6enrA9vT0dBUWFtb4nMLCwmM6fvLkyUpOTvbfsrKy7Cn8IdzuKO03sapwxYbk/AAAHE3jxo3VrVs3vfTSS/6AM2DAAC1btkw//PBDQEtPMNq1a6cXXnhBa9as0ZdffqlRo0YpISHBv3/mzJl65plntGrVKv3000968cUXlZCQoJNOOkmStU7PJ598oi1bthwWlk6UsHdvhdqkSZNUVFTkv/38888heZ1Teg1Swt071OqOVSE5PwAAv2bgwIHyeDz+0JOamqpOnTopIyND7du3P65zP/PMM9q9e7d69eqlyy67TGPHjg0Ye5uSkqKnnnpKp59+urp166YPP/xQb7/9tpo0aSJJuueee7Rx40adfPLJatas2XGVJVguE8aLRpWXlysxMVH/+c9/dOGFF/q35+bmas+ePXrrrbcOe052drYmTJigcePG+bfdeeedmjNnjn9q3NEUFxcrOTlZRUVFSkpKsqMaAIB64sCBA9qwYYNat26t+Pj4cBcn4hzt/bfj+zusLT2xsbHq3bu3f5qcJHm9Xi1YsEA5OTk1PicnJyfgeEmaP3/+EY8HAACQHDBlfcKECcrNzVWfPn3Ut29fTZ06Vfv27dOVV14pSbr88svVokULTZ48WZJ04403auDAgXrwwQc1fPhwzZo1S0uWLNGMGTPCWQ0AAOBwYQ89l1xyiXbs2KE77rhDhYWF6tGjh95//33/YOX8/Hy53QcbpPr376+XX35Zf/vb3/SXv/xF7dq105w5c1ijBwAAHFVYx/SEA2N6AABHwpie8KrXY3oAAHCiCGsPcIxQv++EHgAAqvhWKi4tLQ1zSSKTb3XnqKiokJw/7GN6AABwiqioKKWkpPivQ5WYmCjXMV4OAsHxer3asWOHEhMTQ3a9LkIPAADV+K4A/msX4IT93G63srOzQxY0CT0AAFTjcrnUvHlzpaWlOepimZEgNjY2YMa23Qg9AADUICoqKmRjSxAeDGQGAAARgdADAAAiAqEHAABEhIgb0+Nb+Ki4uDjMJQEAALXl+94+ngUMIy707N27V5KUlZUV5pIAAIBjtXfvXiUnJwf13Ii79pbX69XWrVvVqFEj29cBKC4uVlZWln7++ed6eV2v+l4/iTrWB/W9fhJ1rA/qe/0k++tojNHevXuVmZkZ9LT2iGvpcbvdatmyZUhfIykpqd7+Ekv1v34SdawP6nv9JOpYH9T3+kn21jHYFh4fBjIDAICIQOgBAAARgdBjo7i4ON15552Ki4sLd1FCor7XT6KO9UF9r59EHeuD+l4/yZl1jLiBzAAAIDLR0gMAACICoQcAAEQEQg8AAIgIhB4AABARCD02eeyxx9SqVSvFx8erX79++uqrr8JdJEnS5MmTdeqpp6pRo0ZKS0vThRdeqLVr1wYcc+DAAY0ZM0ZNmjRRw4YN9Yc//EHbtm0LOCY/P1/Dhw9XYmKi0tLSNHHiRFVWVgYck5eXp169eikuLk5t27bVzJkzDytPqN+nKVOmyOVyady4cfWqflu2bNGll16qJk2aKCEhQV27dtWSJUv8+40xuuOOO9S8eXMlJCRo8ODBWrduXcA5du3apVGjRikpKUkpKSm6+uqrVVJSEnDMN998ozPPPFPx8fHKysrS/ffff1hZXn/9dXXo0EHx8fHq2rWr3n333eOun8fj0e23367WrVsrISFBJ598su69996Aa+zUpTp+8sknOu+885SZmSmXy6U5c+YE7HdSXWpTlmOtY0VFhW699VZ17dpVDRo0UGZmpi6//HJt3bq13tTxUNdee61cLpemTp1aZ+pYm/qtWbNG559/vpKTk9WgQQOdeuqpys/P9++vc5+vBsdt1qxZJjY21jz77LPmu+++M6NHjzYpKSlm27Zt4S6aGTJkiHnuuefMqlWrzIoVK8ywYcNMdna2KSkp8R9z7bXXmqysLLNgwQKzZMkSc9ppp5n+/fv791dWVpouXbqYwYMHm+XLl5t3333XNG3a1EyaNMl/zE8//WQSExPNhAkTzOrVq82jjz5qoqKizPvvv+8/JtTv01dffWVatWplunXrZm688cZ6U79du3aZk046yVxxxRXmyy+/ND/99JOZN2+eWb9+vf+YKVOmmOTkZDNnzhyzcuVKc/7555vWrVub/fv3+48599xzTffu3c0XX3xhPv30U9O2bVszcuRI//6ioiKTnp5uRo0aZVatWmVeeeUVk5CQYJ588kn/MZ999pmJiooy999/v1m9erX529/+ZmJiYsy33357XHW87777TJMmTczcuXPNhg0bzOuvv24aNmxopk2bVifr+O6775q//vWvZvbs2UaSefPNNwP2O6kutSnLsdZxz549ZvDgwebVV18133//vVm8eLHp27ev6d27d8A56nIdq5s9e7bp3r27yczMNA8//HCdqeOv1W/9+vUmNTXVTJw40SxbtsysX7/evPXWWwGfaXXt85XQY4O+ffuaMWPG+B97PB6TmZlpJk+eHMZS1Wz79u1Gklm4cKExxvpwiomJMa+//rr/mDVr1hhJZvHixcYY6x+G2+02hYWF/mOeeOIJk5SUZMrKyowxxtxyyy2mc+fOAa91ySWXmCFDhvgfh/J92rt3r2nXrp2ZP3++GThwoD/01If63XrrreaMM8444n6v12syMjLMP//5T/+2PXv2mLi4OPPKK68YY4xZvXq1kWS+/vpr/zHvvfeecblcZsuWLcYYYx5//HHTuHFjf519r92+fXv/44svvtgMHz484PX79etn/vd///e46jh8+HBz1VVXBWy76KKLzKhRo+p8HQ/9MnFSXWpTlmDqWJOvvvrKSDKbNm2qV3XcvHmzadGihVm1apU56aSTAkJPXapjTfW75JJLzKWXXnrE59TFz1e6t45TeXm5li5dqsGDB/u3ud1uDR48WIsXLw5jyWpWVFQkSUpNTZUkLV26VBUVFQHl79Chg7Kzs/3lX7x4sbp27ar09HT/MUOGDFFxcbG+++47/zHVz+E7xneOUL9PY8aM0fDhww8rQ32o33//+1/16dNHI0aMUFpamnr27KmnnnrKv3/Dhg0qLCwMeO3k5GT169cvoI4pKSnq06eP/5jBgwfL7Xbryy+/9B8zYMAAxcbGBtRx7dq12r17d63eh2D1799fCxYs0A8//CBJWrlypRYtWqShQ4fWmzr6OKkutSmLXYqKiuRyuZSSklJv6uj1enXZZZdp4sSJ6ty582H763IdvV6v3nnnHZ1yyikaMmSI0tLS1K9fv4AusLr4+UroOU47d+6Ux+MJ+B8qSenp6SosLAxTqWrm9Xo1btw4nX766erSpYskqbCwULGxsf4PIp/q5S8sLKyxfr59RzumuLhY+/fvD+n7NGvWLC1btkyTJ08+bF99qN9PP/2kJ554Qu3atdO8efN03XXXaezYsfr3v/8dUMajvXZhYaHS0tIC9kdHRys1NdWW9+F463jbbbfpj3/8ozp06KCYmBj17NlT48aN06hRo+pNHX2cVJfalMUOBw4c0K233qqRI0f6LzxZH+r4j3/8Q9HR0Ro7dmyN++tyHbdv366SkhJNmTJF5557rj744AP9/ve/10UXXaSFCxf6X7eufb5G3FXWI9mYMWO0atUqLVq0KNxFsc3PP/+sG2+8UfPnz1d8fHy4ixMSXq9Xffr00d///ndJUs+ePbVq1SpNnz5dubm5YS6dPV577TW99NJLevnll9W5c2etWLFC48aNU2ZmZr2pY6SqqKjQxRdfLGOMnnjiiXAXxzZLly7VtGnTtGzZMrlcrnAXx3Zer1eSdMEFF2j8+PGSpB49eujzzz/X9OnTNXDgwHAWL2i09Bynpk2bKioq6rDR6tu2bVNGRkaYSnW4G264QXPnztXHH3+sli1b+rdnZGSovLxce/bsCTi+evkzMjJqrJ9v39GOSUpKUkJCQsjep6VLl2r79u3q1auXoqOjFR0drYULF+qRRx5RdHS00tPT63T9JKl58+bq1KlTwLaOHTv6Z1D4zn+0187IyND27dsD9ldWVmrXrl22vA/HW8eJEyf6W3u6du2qyy67TOPHj/e33tWHOvo4qS61Kcvx8AWeTZs2af78+f5WHt9r1+U6fvrpp9q+fbuys7P9nz2bNm3STTfdpFatWtX5OjZt2lTR0dG/+tlT1z5fCT3HKTY2Vr1799aCBQv827xerxYsWKCcnJwwlsxijNENN9ygN998Ux999JFat24dsL93796KiYkJKP/atWuVn5/vL39OTo6+/fbbgH+8vg8w3z+InJycgHP4jvGdI1Tv01lnnaVvv/1WK1as8N/69OmjUaNG+X+uy/WTpNNPP/2wZQZ++OEHnXTSSZKk1q1bKyMjI+C1i4uL9eWXXwbUcc+ePVq6dKn/mI8++kher1f9+vXzH/PJJ5+ooqIioI7t27dX48aNa/U+BKu0tFRud+DHUVRUlP+vzfpQRx8n1aU2ZQmWL/CsW7dOH374oZo0aRKwv67X8bLLLtM333wT8NmTmZmpiRMnat68eXW+jrGxsTr11FOP+tlTJ78/jmnYM2o0a9YsExcXZ2bOnGlWr15t/vSnP5mUlJSA0erhct1115nk5GSTl5dnCgoK/LfS0lL/Mddee63Jzs42H330kVmyZInJyckxOTk5/v2+KYfnnHOOWbFihXn//fdNs2bNapxyOHHiRLNmzRrz2GOP1Tjl8ES8T9Vnb9WH+n311VcmOjra3HfffWbdunXmpZdeMomJiebFF1/0HzNlyhSTkpJi3nrrLfPNN9+YCy64oMYp0D179jRffvmlWbRokWnXrl3A1Nk9e/aY9PR0c9lll5lVq1aZWbNmmcTExMOmzkZHR5sHHnjArFmzxtx55522TFnPzc01LVq08E9Znz17tmnatKm55ZZb6mQd9+7da5YvX26WL19uJJmHHnrILF++3D9zyUl1qU1ZjrWO5eXl5vzzzzctW7Y0K1asCPjsqT5LqS7XsSaHzt5yeh1/rX6zZ882MTExZsaMGWbdunX+qeSffvqp/xx17fOV0GOTRx991GRnZ5vY2FjTt29f88UXX4S7SMYYaxpiTbfnnnvOf8z+/fvN9ddfbxo3bmwSExPN73//e1NQUBBwno0bN5qhQ4eahIQE07RpU3PTTTeZioqKgGM+/vhj06NHDxMbG2vatGkT8Bo+J+J9OjT01If6vf3226ZLly4mLi7OdOjQwcyYMSNgv9frNbfffrtJT083cXFx5qyzzjJr164NOOaXX34xI0eONA0bNjRJSUnmyiuvNHv37g04ZuXKleaMM84wcXFxpkWLFmbKlCmHleW1114zp5xyiomNjTWdO3c277zzznHXr7i42Nx4440mOzvbxMfHmzZt2pi//vWvAV+QdamOH3/8cY3/7nJzcx1Xl9qU5VjruGHDhiN+9nz88cf1oo41qSn0OLmOtanfM888Y9q2bWvi4+NN9+7dzZw5cwLOUdc+X13GVFvyFAAAoJ5iTA8AAIgIhB4AABARCD0AACAiEHoAAEBEIPQAAICIQOgBAAARgdADAAAiAqEHAABEBEIPgDrrrrvuUo8ePY56zBVXXKELL7zwhJQHgLNFh7sAABBK06ZNU/WF5wcNGqQePXpo6tSp4SsUgLAg9ACoc4wx8ng8tTo2OTk5xKUBUFfQvQXAEcrKyjR27FilpaUpPj5eZ5xxhr7++mtJUl5enlwul9577z317t1bcXFxWrRokf+5Tz75pLKyspSYmKiLL75YRUVF/n3Vu7euuOIKLVy4UNOmTZPL5ZLL5dLGjRu1e/dujRo1Ss2aNVNCQoLatWun55577oTWH0DoEXoAOMItt9yiN954Q//+97+1bNkytW3bVkOGDNGuXbv8x9x2222aMmWK1qxZo27dukmS1q9fr9dee01vv/223n//fS1fvlzXX399ja8xbdo05eTkaPTo0SooKFBBQYGysrJ0++23a/Xq1Xrvvfe0Zs0aPfHEE2ratOkJqTeAE4fuLQBht2/fPj3xxBOaOXOmhg4dKkl66qmnNH/+fD3zzDM69dRTJUn33HOPzj777IDnHjhwQM8//7xatGghSXr00Uc1fPhwPfjgg8rIyAg4Njk5WbGxsUpMTAzYl5+fr549e6pPnz6SpFatWoWqqgDCiJYeAGH3448/qqKiQqeffrp/W0xMjPr27as1a9b4t/lCSXXZ2dn+wCNJOTk58nq9Wrt2ba1f/7rrrtOsWbPUo0cP3XLLLfr888+DrAkAJyP0AKgzGjRoEJLzDh06VJs2bdL48eO1detWnXXWWbr55ptD8loAwofQAyDsTj75ZMXGxuqzzz7zb6uoqNDXX3+tTp06HfW5+fn52rp1q//xF198Ibfbrfbt29d4fGxsbI0zv5o1a6bc3Fy9+OKLmjp1qmbMmBFkbQA4FWN6AIRdgwYNdN1112nixIlKTU1Vdna27r//fpWWlurqq6/WypUrj/jc+Ph45ebm6oEHHlBxcbHGjh2riy+++LDxPD6tWrXSl19+qY0bN6phw4ZKTU3VXXfdpd69e6tz584qKyvT3Llz1bFjx1BVF0CYEHoAOMKUKVPk9Xp12WWXae/everTp4/mzZunxo0bH/V5bdu21UUXXaRhw4Zp165d+t3vfqfHH3/8iMfffPPNys3NVadOnbR//35t2LBBsbGxmjRpkjZu3KiEhASdeeaZmjVrlt1VBBBmLlN9qVIAAIB6ijE9AAAgIhB6AABARCD0AACAiEDoAQAAEYHQAwAAIgKhBwAARARCDwAAiAiEHgAAEBEIPQAAICIQegAAQEQg9AAAgIjw/wE4REV40T1EkwAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<Figure size 640x480 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig, ax = plt.subplots()\n", | |
| "ax.set_xlabel(\"orbits\")\n", | |
| "ax.set_ylabel(\"megno\")\n", | |
| "ax.plot(times/(2*np.pi), megno, label=\"ias15\")\n", | |
| "ax.plot(times/(2*np.pi), megno_whfast, label=\"whfast\")\n", | |
| "ax.legend()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "id": "5796b66b-75b4-4d57-9653-3fa8e8174c7a", | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3 (ipykernel)", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.9.6" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 5 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment