Skip to content

Instantly share code, notes, and snippets.

@hannorein
Last active April 2, 2021 20:32
Show Gist options
  • Save hannorein/83ccbc13d9d5292d80c9802e9d1d6250 to your computer and use it in GitHub Desktop.
Save hannorein/83ccbc13d9d5292d80c9802e9d1d6250 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "educational-contrast",
"metadata": {},
"outputs": [],
"source": [
"'''Let's try calculating the Lyapunov timescale using Rebound!'''\n",
"import matplotlib.pyplot as plt\n",
"import rebound\n",
"import numpy as np\n",
"from numpy import radians as rd"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "partial-adams",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Searching NASA Horizons for 'Sun'... Found: Target body name: Sun (10).\n",
"Searching NASA Horizons for 'Mercury'... Found: Target body name: Mercury Barycenter (199).\n",
"Searching NASA Horizons for 'Venus'... Found: Target body name: Venus Barycenter (299).\n",
"Searching NASA Horizons for 'Earth'... Found: Target body name: Earth-Moon Barycenter (3).\n",
"Searching NASA Horizons for 'Mars'... Found: Target body name: Mars Barycenter (4).\n",
"Searching NASA Horizons for 'Atalante'... Found: Target body name: 36 Atalante (A855 TA).\n",
"Searching NASA Horizons for 'Jupiter'... Found: Target body name: Jupiter Barycenter (5).\n",
"Searching NASA Horizons for 'Saturn'... Found: Target body name: Saturn Barycenter (6).\n",
"Searching NASA Horizons for 'Uranus'... Found: Target body name: Uranus Barycenter (7).\n",
"Searching NASA Horizons for 'Neptune'... Found: Target body name: Neptune Barycenter (8).\n",
"Searching NASA Horizons for 'Pluto'... Found: Target body name: Pluto Barycenter (9).\n"
]
}
],
"source": [
"\n",
"# Note: Atalante is added after Mars. This is because by default WHFast uses Jacobi coordinates.\n",
"\n",
"bodies = ['Sun', 'Mercury', 'Venus', 'Earth', 'Mars', \"Atalante\", # Use these bodies\n",
" 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto']\n",
"date = '2013-08-09 12:00' # Start on this date\n",
"\n",
"# Set up simulation.\n",
"sim = rebound.Simulation()\n",
"sim.add(bodies, date=date) # Add the Sun and all the planets\n",
"sim.move_to_com() # Move to the center of momentum frame\n",
"sim.save(\"ss.bin\")"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "waiting-little",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-45-b3a9471059ff>:21: RuntimeWarning: divide by zero encountered in double_scalars\n",
" lyapunov_CN[i] = 1./time * np.abs(d)\n"
]
}
],
"source": [
"\n",
"twopi = rd(360) # 6.283185307179586\n",
"dt_out = 20 * twopi # Years * 2 Pi, output interval\n",
"tmax = 1e5 * twopi # Years * 2 Pi, end time\n",
"times = np.arange(0., tmax, dt_out) # Output times, in units yrs * 2 Pi\n",
"timeyr = times / twopi # The times in human-friendly years\n",
"megno = np.zeros(len(times)) # An empty array to store megno values\n",
"lyapunov_CN = megno.copy() # An empty array to store LCN values\n",
"\n",
"# Load in cached simulation\n",
"sim = rebound.Simulation(\"ss.bin\")\n",
"\n",
"sim.integrator = 'whfast' # Use WHfast\n",
"orbits = sim.calculate_orbits() # Calculate the orbits\n",
"sim.dt = 0.05 * orbits[0].P # Use a timestep = 5% of Mercury's period\n",
"\n",
"# Create a set of variational/shadow particles\n",
"var = sim.add_variation()\n",
"\n",
"# Perturb the orbit of the asteroid (particle number 5).\n",
"# Ideally, perturb each position and velocity coordinates\n",
"# by a random amount (here only x is perturbed).\n",
"# Note that the initial perturbation has magnitude 1. \n",
"var.particles[5].x = 1. \n",
"\n",
"# Now integrate, output, repeat.\n",
"for i, time in enumerate(times):\n",
" sim.integrate(time, exact_finish_time=0) # Integrate\n",
" # Calculate divergence of shadow particle\n",
" d = np.sqrt(np.sum(np.square(var.particles[5].xyz + var.particles[5].vxyz) ))\n",
" lyapunov_CN[i] = 1./time * np.abs(d)/1. # divide by initial perturbation (here 1.)\n"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "reverse-monster",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAENCAYAAAARyyJwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3Z0lEQVR4nO3dd3hUZfbA8e9JJ5TQewlSpZfQBBQUFUUEsXewYK9r3dXfuusq1l3XLq6C7qpgVxRQUZEuTRCQJr1Ih1BDy/n9MZMwSSYzdzJzM5PkfJ5nHnLv3LlzuJA5c99yXlFVjDHGmMLERTsAY4wxsc0ShTHGmIAsURhjjAnIEoUxxpiALFEYY4wJyBKFMcaYgCxRGGOMCcgShTHGmIASoh2AEyIyGBgAVALeUtVvoxuRMcaUHeL2zGwReRs4D9imqm189vcH/g3EA/9R1accnKsK8JyqXh/ouOrVq2t6enpYcRtjTFkzb968HapaI//+4rijGA28DLybs0NE4oFXgDOBjcAcEfkST9IYke/116nqNu/Pj3hfF1B6ejpz584NP3JjjClDRGSdv/2uJwpVnSIi6fl2dwV+V9XVACIyBhikqiPw3H3kISICPAVMUNX5LodsjDHGR7Q6s+sBG3y2N3r3FeYOoB9wkYjc7O8AERkuInNFZO727dsjF6kxxpRxJaIzW1VfBF4McsxIYCRARkaGlcQ1xpgIidYdxSaggc92fe8+Y4wxMSZaiWIO0ExEGotIEnAZ8GW4JxWRgSIyMjMzM+wAjTHGeLieKETkA2Am0EJENorI9ap6DLgd+AZYCnyoqkvCfS9VHaeqw9PS0sI9lTHGGK+gfRQiEge0B+oCh4DFPsNVg1LVywvZPx4Y7/Q8xWHV9v3s3H+Ero2rRjsUY4yJGYUmChFpAjyIZ7TRSmA7kAI0F5GDwBvAO6qaXRyBOiEiA4GBTZs2LdLrn56wjO+XbePB/i24sfdJeEblGmNM2VbozGxvk9GrwDTNd5CI1ASuAHar6juuRxmijIwMLcqEu71ZR3nw41+ZsHgLZ7WqxbMXtyetXKILERpjTOwRkXmqmlFgf6ASHt5mp+6qOsPN4CKtqIkCQFV5e/paRoxfSt3K5Xj1yk60qWd9HsaY0q+wRBGwM9vbrBS0ZEZpIiJc36sxY2/qzpFj2Qx5bQZj56zH7ZpYxhgTq5yMevpeRC6UEtBgH8nhsZ0bVeXrO3vRrXFVHvxkEfd//CuHjhyPQJTGGFOyBK0eKyL7gPLAMSALEEBVtZL74RVNOE1P+R3PVl78fiUv/rCSFrUq8uqVnTipRoWInNsYY2JJkZqeAFS1oqrGqWqSqlbybsdskoi0+DjhnjObM2poF7buzeL8l6czftEf0Q7LGGOKTdBEISKfiMi53o7tMqtPi5p8fWdvmtWqwK3vzefv437jyLGYGRlsjDGucfLh/xpwJbBSRJ4SkRYux1RkbpfwqFu5HGOH92BYz3Tenr6Gy0bO5I/MQ668lzHGxAonTU+TVPVKoBOwFpgkIjNEZJiIxNQkg+Io4ZGUEMdfB7bm5Ss6snzLPga8OI2pK62suTGm9HLUnCQi1YChwA3AL3iWMO0EfOdaZDHuvHZ1+fKOXlSvkMQ1b8/mhUkrOJ5tQ2iNMaWPkz6Kz4CpQCowUFXPV9WxqnoHUKaH/zSpUYHPb+vJBR3q8cKklQwdNZtdB45EOyxjjImoQhOFiNT1/viiqrZS1RGqmme4j79hVGVNalICz1/SnhFD2vLzml0MeHEq89fvjnZYxhgTMYHuKP4jIrOAs0Wkj4jE/Gp40VqPQkS4vGtDPr3lFBLihUten8mo6WtsNrcxplQIVuspBegDnAP0BNYDE4GJqrq+OAIsikhOuAtV5sGj/OmjhUxaupUBbevw1IVtqZgSU33+xhjjV5GKAvo5SWM8SaM/UFtVu0YuxMiJZqIAT2HBN6as5tlvltOoaiqvXtWJlrXLzBxFY0wJVeSZ2b5UdY2qvqqq5wO9IhZdKSMi3HxaE967oRv7Dh9j8CvT+WTexmiHZYwxRRKoM3ua9899IrI3/5+qasN7guh+UjW+vrMXHRpU5k8fLeThT38l66gVFjTGlCyFJgpV7eX9s6JvjaeyVuspXDUrpvC/67txa58mfDB7Axe+NoP1Ow9GOyxjjHHM6YS7TiJyp4jcISId3Q6qqKI16imYhPg4HujfkreuzWDDroOc99JUflzmeNlxY4yJKicT7v4PeAeoBlQHRovII24HVhTFUcIjHGecXIuv7+xN/SqpXPfOHP49aSXZNpvbGBPjnKxHsRxor6pZ3u1ywAJVjdnigNEe9RTMoSPH+ctni/j0l02c0bIm/7y0g63NbYyJunBGPW0GUny2k4FNkQqsLCqXFM/zl7Tn74Na89OK7Zz/8jSWbdkb7bCMMcYvJ4kiE1giIqNFZBSwGNgjIi+KyIvuhld6iQjX9EhnzPDuHDpynAtemcEXCyz/GmNij5OyHJ95HzkmuxNK2ZSRXpWv7uzFbe/N564xC1i4IZOHz21JYnyZXifKGBNDgiYKVX2nOAIpy2pWTOH9G7vzxNdLeXv6GhZvzuTlKzpSs2JK8BcbY4zLAk24qyQiI0TkvyJyRb7nXnU/tLIlMT6Ox85vzQuXduDXjXsY+NI05q2zKrTGmOgL1L4xChDgE+Ay79rZyd7nurseWRHE6jyKUAzuWI/Pbu1JckI8l42cyX9nrbMqtMaYqAqUKJqo6kOq+rm3ttN84AfvancxKdbnUTh1cp1KjLu9F72aVufRzxdz30dW+sMYEz2BEkWyiOQ+r6pPAG8CU/BMvjMuSktN5K1ru3B3v2Z8Mn8jF742gw27rPSHMab4BUoU44DTfXeo6mjgT4AVBCwGcXHC3f2a89a1GazfdZCBL09jyort0Q7LGFPGBCoK+ICqTvKzf6KqNnM3LOPrjJNrMe72XtSulMK1o2bz8g9W+sMYU3wCjXq6yrfpyc/zTUTE1qQoJunVy/Ppradwfvu6PPftCm763zz2Zh2NdljGmDIg0DyKasAvIjIPmAdsx1PKoylwGrADeMj1CE2u1KQEXri0A+3rV+bJ8UsZ/PJ0Xr+6M81rVYx2aMaYUixQ09O/gU7AB0AN4Azv9ibgalW9UFVXFkuUJpeIcF2vxrx/Y3f2ZnlWz/vq183RDssYU4qFtGZ2SRHr1WMjZeveLG753zzmr9/D8FNP4oGzW5BgpT+MMUUUkTWzTWypVSmFMcN7cE2PRoycspqr35rNjv2Hox2WMaaUsURRwiUlxPH3QW14/uL2zF+/m4EvTeOX9Vb6wxgTOQEThYjEicglxRWMKboLO9fnk1tOIT5OuPSNWfx35lor/WGMiYiAiUJVs4EHiimWsJWGWk/haFMvja/u6EXPptV49Isl3Pb+fDIP2RBaY0x4nDQ9TRKR+0SkgYhUzXm4HlkRlJZaT+GonJrEW9d24eFzWvLtkq2c99JUFm7YE+2wjDElmJNEcSlwG54aTzlzKkr/kKISLC5OuOm0Joy9qQfZ2XDR6zP4z9TV1hRljCmSoIlCVRv7eZxUHMGZ8HRuVIXxd/amb4ua/OPrpdz47lz2HLQyXcaY0ARNFCKSKiKPiMhI73YzETnP/dBMJKSlJvLG1Z3568BW/LRiO+f+eyrz1u2KdljGmBLESdPTKDzVYk/xbm8C/uFaRCbiRIRhPRvzyS2nkBAfxyVvzOK1yaussKAxxhEniaKJqj4DHAVQ1YN4Vr4zJUy7+pX56s5e9G9Tm6cnLmPo6Dk2Qc8YE5STRHFERMoBCp6qsYB9upRQlVISefnyjjxxQRtmrd7Juf+eyqzVO6MdljEmhjlJFH8FJgINROQ94HtK0NwKU5CIcGW3RnxxW08qpCRwxZuz+PeklRy3pihjjB+OigJ618nujqfJaZaq7nA7sHCUlaKAkXDg8DEe/Xwxn/6yiVOaVOOFSztQs1JKtMMyxkRBuEUBT8NTZrwv0DuSgZnoKp+cwPOXtOfZi9rxy/o9nPfSNH7fti/aYRljYoiT4bGvAjcDi4DFwE0i8orbgZniIyJcnNGAz247hWyFy0bOYvkWSxbGGA8ndxSnA2er6ihVHQWc691nSpmWtSsx9qbuxMcJl42cyZLNZbNmljEmLyeJ4negoc92A+8+Uwo1qVGBscN7UC4xnive/JlfN+6JdkjGmCgrNFGIyDgR+RKoCCwVkcki8iOw1LvPlFLp1csz9qYeVExJ4Mo3f2beOlvfwpiyrNBRTyJyWqAXqupPrkRUMI6TgbuA6sD3qvpasNfYqKfI2LznEFe8OYvt+w7zt0FtuKBjPeLjbK6lMaVVYaOeHK+ZLSKVgIScbVUNWjBIRN4GzgO2qWobn/39gX8D8cB/VPUpB+eKA95V1auCHWuJInK27s1i+LtzWbgxk6Y1K3DfWc05u3VtRCxhGFPaFHl4rIgMF5EtwK94youHUmZ8NNA/3/nigVeAc4BWwOUi0kpE2orIV/keNb2vOR/4Ghjv8H1NhNSqlMLnt/XktSs7oarc/L/5nP/ydKau3B7t0IwxxSToHYWIrAR6FHWSnYikA1/l3FGISA/gMVU927v9MICqjnBwrq9VdUCw4+yOwh3Hjmfz2S+beGHSSjbtOcRNp57EA/1bWnOUMaVEYXcUCf4OzmcVcDCCsdQDNvhsbwS6FXawiPQBhgDJBLijEJHhwHCAhg0bFnaYCUNCfBwXZzTg/A51+cdXS3ljympW7zjAC5d2oHyyk/9KxpiSyMlv98PADBH5GZ9igKp6p2tR+VDVycBkB8eNBEaC547C3ajKtuSEeB4f3IYmNcrz969+4+LXZ/LW0AzqpJWLdmjGGBc4mUfxBvADMIsTS6HOC+M9N+GZi5GjvnefKWGG9mzMW0O7sH7XQQa9PJ3Fm2yCnjGlkZNEkaiq93pnZr+T8wjjPecAzUSksYgkAZcBX4ZxvlwiMlBERmZm2gdWcenboiaf3HIKifFxXD5yFj9byXJjSh0niWKCd+RTHRGpmvNwcnIR+QCYCbQQkY0icr2qHgNuB77BM3nvQ1VdUuS/gQ9VHaeqw9PS0iJxOuNQi9oV+ejmHtSslMw1b89m0m9box2SMSaCnIx6WuNnt6rqSe6EFD4b9RQduw4cYeio2SzZvJfnLm7HBR3rRzskY0wIijyPQlUb+3nEZJKwpqfoqlo+ifdv7E63xlW5Z+xC/j7uNw4fOx7tsIwxYXJyR3GNv/2q+q4rEUWA3VFEV9bR44wYv5R3Zq7j5DqVeOnyDjStaeXBjIl14Sxc1MXn0Rt4DDg/otGZUiUlMZ6/DWrDW9dmsHVvFue9NI3nv13Oxt2RnI5jjCkujms95b5ApDIwRlX7Bzs2WuyOInZs25fFo58v5ltvB3fvZjW498zmdGhQObqBGWMKCLsooM+JEoHFqtoiUsFFiogMBAY2bdr0xpUrV0Y7HONj4+6DfDh3Ix/MXo8qfH/vaaSlJkY7LGOMj3CKAo4TkS+9j6+A5cBnbgQZLhseG7vqV0nl3jObM2poF3YfPMKT45dGOyRjjENOSng85/PzMWCdqm50KR5TyrWpl8YNvRrzxpTVDO5Yjx5NqkU7JGNMEE6Gx/7k85huScKE6+5+zWlYNZU/f7aIrKM2fNaYWBdohbs1QGEdGKqqTVyLqoisj6LkmLpyO1e/NZtalZJpVacSbetX5vpejUkrZ/0WxkRLyJ3ZIpK/TSAOuAS4D5ivqhdGPMoIsVFPJcMXCzbx47JtLNuyjxVb91G7UgrPX9LBmqOMiZIij3ryLkF6NXA/sAB4UlV/cyPISLFEUfIs2LCHu8f8wrpdB/nzOSdz46kxOfnfmFIt5FFPIpIoIjcBv+GZaDdYVa+K9SRhSqYODSrz9Z296d+6Nk+MX8q3S7ZEOyRjjFegUU9r8IxyegFYD7QTkXY5T6rqp+6GZsqa8skJ/OvSDmzeM5N7xi7gs9t60ryWlf4wJtoCjXqaBPwItAcG5nuc535oobOigCVfSmI8b1ydQWpyAte/M4dNew4BoKrsPnAkytEZUzaFPDO7JLA+ipJv4YY9XPXWz6SVS+S5i9vz6uRVzPh9B5/f1pM29WxCpTFuCKcooDHFrn2Dyrx/Q3f2ZR3jspGzmLt2FwnxwtvT/S2PYoxxkyUKE7Pa1k9j7E3dGXpKOt/cfSqXZDTgq4V/sGP/4WiHZkyZYonCxLSWtSvx2PmtaVA1lWt6pHPkeDYf/Lw+2mEZU6Y4KQp4m7e0eM52FRG51dWojPGjac0K9G5WnXdnrWPV9v0AbNh1kMyDR6McmTGlm5MJdwtUtUO+fb+oakc3AysKK+FR+s1du4uho+Zw6OhxmtQoz4qt+2lRqyKf39aTcknx0Q7PmBItnM7seBERnxPFA0mRDC5SrMx46ZeRXpXJ9/fhym4NqZyaxI29G7Ni2z4e+3JJtEMzptRyUmZ8IjBWRN7wbt/k3WdMVFSvkMzfB7XJ3U5KiOOVH1fRv21t+jSvwX9nreOMk2tRr3K5KEZpTOnh5I7iQTwT727xPr4HHnAzKGNCcXe/5tSsmMx/Z65j1upd/N8XSxhhCyMZEzFB7yhUNRt4zfswJuYkxsdxSUYDXp38O3sOemZvT1i8hY27D1K/SmqUozOm5AtUFPBD75+LROTX/I/iC9GY4C7t0oBshfnr9zCkYz0ARk9fG92gjCklAt1R3OX9MybrOhnjq0HVVHo3q86MVTu5v38LjmUrY+Zs4K5+zZj++04qlUvglCbVox2mMSVSoYlCVf/w/nirqj7o+5yIPI2n78KYmPHE4Las2rGfOmnluKF3Y75cuJmnJixj7JwNKPDMhe24sHP9aIdpTInjpDP7TD/7zol0IJFg1WPLtobVUunboiYA7epXpmvjqrz383pSEuPp1LAyf/5sEYeP2RrdxoQqUB/FLSKyCGiZr39iDRCTfRQ2j8L4urG3Z5W8W/o0YVjPxhw+ls28dbvp98+fePH7lZTGysnGuCFQH8X7wARgBPCQz/59qrrL1aiMiYB+J9dkzPDudEmvymbvuhZvT1vD79v288/vVlC1fBJXdW8U5SiNiX2F3lGoaqaqrgUeAbao6jqgMXCVb+0nY2KViND9pGrExwn1q5Sjcmoik5ZuIyk+jsbVyzNp6VYA9h8+xsbdB6McrTGxy0kfxSfAcRFpCowEGuC52zCmxBAR2noXPOrQsDLt66exYss+jh7PptPj39Hr6R+jHKExsctJoshW1WPAEOAlVb0fqONuWMZEXrv6nkTR46RqNK9dkc2ZWTw1YRlHjmUDsDfrKIs3ZfLJvI3RDNOYmOOk1tNREbkcuAbPetkAie6FZIw7OjeqAkDvZtXZm+UpTT56xtrc51dvP8DgV6YDMKRTPXxqYRpTpjm5oxgG9ACeUNU1ItIY+K+7YRkTeX1b1OSrO3qRkV6V5rUqAnA8W7n5tCYATF2xPffY7fsPs+fgEfZl2VoXxgRNFKr6m6reqaofeLfXqOrT7odmTGSJCG28/RT1KpejfFI8cQLDeqYTHyc8/92K3GM37j7EsNFzeODjX9mw6yAnPzqRJZttfo4pmwptehKRD1X1Eu9cigIDzlW1nauRGeMiEaFDw8rEiVCrUgr1Kpdj/a6DVCufxM4DR/ht814WbNhDnUopLP1jL4eOHmfRxkxa17U5OqbssVpPpsx67arO5PRCrN/lGR770Dktuf/jX/liwSZUYXNmFqu2HwDInYthTFkTaB7FH97V7Ear6rr8j2KM0TEr4WFCUSklkYopnnEZd/drRpzA4I71qFo+iTlrd+ceN3Wlp+9ic2YW3y7ZwtSV2zl6PDu3Q9yY0s7JmtnfA0NUtcR8+mZkZOjcuXOjHYYpYbKzlbg44bRnf2TdzoM0rJqae6cBcEqTasxYtROA/q1rM3HJFtaMONdGR5lSI5w1s/cDi0TkLRF5MecR+RCNia64OM8HfnKC59fi5Ss6Ui4xPvf5jbtPND1NXLIFgE3WHGXKACfzKD71PowpE16+ohObdh/KrUD7k3fYrO/dRY4VW/fx+S+bSIiPyx1ma0xp42Qp1HdEpBzQUFWXF0NMxkRV81oVc+dZnNmqVm6i8Ofn1bt4Y8pqgNxEcd9HC+nVtDqDvSvtGVPSBW16EpGBwAJgone7g4h86XJcxsSE01t61rfo3ezE6nhxPl0SOUkC4MuFm1FVPp63kbvHLiiuEI1xnZM+iseArsAeAFVdAJzkWkTGxJC6lcsx9YG+/OfaE/17OZP28rvzg1949pu8N937Dx9j9hqrym9KNieJ4qifEU/ZbgRjTCxqUDWV5IR4epxUDYDalVIKPfaT+XkLCt49ZgGXvDGTXQeOuBqjMW5ykiiWiMgVQLyINBORl4AZLsdlTMx557qu/PrYWaSVK7wm5ta9h3N/vvbt2blrXhw4fMz1+Ixxi5NEcQfQGjiMZx2KTE7M2jamzEhKiKNSSiIJ8Z5fm3qVywU83rcTfNeBI2QdPc74RX/YEqymxHEyPHaAqv4F+EvODhG5GPjItaiMiWHePEHDqqmO51EM8pYvB/jXpe25oGN9N0IzxhVO7igedrjPmDIhzjsTO1ATVCArt+6PZDjGuC5Q9dhzgHOBevlmYlcCrMHVlFk5M7eb16rAxCVwavMaTAkw1yK/nfuPsHzLPlrUruhWiMZEVKA7is3AXCALmOfz+BI42/3QjIlNt/dtxpXdGnJznyasfWoA13RvFNLrx87dwNkvTGH5ln0uRWhMZBV6R6GqC4GFIvIZcEBVjwN4K8omF1N8xsSctNREnrigbe52anJ8gKMLd/YLU/JsT76vD+nVy4cVmzFucNJH8S3gO7yjHDDJnXD8E5HyIjJXRGxtDBNzUpOcjAkJbsrK7Yyds57vftsakfMZEylOEkWKqub2vnl/TnVychF5W0S2icjifPv7i8hyEfldRB5ycKoHgQ+dvKcxxa18kueOom5aCved1ZwruzUEoH6VwMNn81u74yAPfrKIG9+1EvkmtjhJFAdEpFPOhoh0BpzWVh4N9Pfd4W26egU4B2gFXC4irUSkrYh8le9RU0TOBH4Dtjl8T2OKVc68ivh44fbTm5Ho3W5Vp1JI53l7+prcn48ez2a3zeY2McLJPfPdwEcishkQoDZwqZOTq+oUEUnPt7sr8LuqrgYQkTHAIFUdgZ9lV0WkD1AeT1I5JCLjVdVKiJiYUTU1CYDz29cFTiSIAe3q8G0Rm5Hu/2ghny/YbAsjmZjgpMz4HBFpCbTw7lququGsAVkP2OCzvRHoFuD9/wIgIkOBHYUlCREZDgwHaNiwYRjhGROatNREFv7fWVRM8fw6XZxRn4z0KqRXK89z3y7ngo71WbIpk++XOb8p/nzBZgC+WLCZ05rXoEr5JFdiN8YJp71wLfB8o08BOokIqvque2EVpKqjgzw/EhgJnqVQiyMmY3KkpZ6YfCcinFSjAgBTHzgdgOPZSpM/jw/5vHePXUDvZtW55bQmNK1ZgZoBChIa45agiUJE/gr0wZMoxuPpW5gGFDVRbAIa+GzX9+4zptSKjyt689HUlTuYunIHDaqWy008xhQnJ53ZFwFnAFtUdRjQHvBfkN+ZOUAzEWksIknAZXgm8YVNRAaKyMjMzPxV0Y0p+TbssvW5TXQ4SRSHvP0Cx0SkEp7RRw2CvAYAEfkAmAm0EJGNInK9qh4Dbge+AZYCH6rqkqKFn5eqjlPV4Wlp4eQxY2LXsi172R+kZLmVNDeR5iRRzBWRysCbeEp4zMfz4R+Uql6uqnVUNVFV66vqW97941W1uao2UdUnihq8MSXRZ7eeUuTX9n9hKpeNnEnnx79j7tqCK+f9uGwbrf/6jd/nfKkqr/+0ih37Dwc8zhgIkijEMy5vhKruUdXXgTOBa71NUMaYEHRNr8rA9nXp2LBKWOdZvGkvOw8c4d/fryzw3NSVOwC4/f1f+GjuhgLP51iwYQ9PTVjGfR8tDCsWUzYETBTqWWFlvM/2WlX91fWoisj6KEws+/DmHrx0eceInW//4WN8PG9jnoWQsr0/b9mbxf0f/8pFr82g+V8mFEgah495RplPXr6dOWt3kXkwnBHvprRz0vQ0X0S6uB5JBFgfhSlLflm/h/s+WsgTXy8FPM1Ja3YcyHPM3HW7OXI8m8e+XMKBw8cY9PI0/jtzLXsOnpj1ffHrMxn48rRijd2ULBJsWUYRWQY0BdYBB/DMzlZVbed+eEWTkZGhc+davRwTu3buP8yHczfy9MRl/P7EOYz7dTN108qxY/8ROjWqTI8RP4R0vrb10li0Kbw76bVPDQjr9YePHScpPs5mkpdgIjJPVTPy73cy4a7ErD0hIgOBgU2bNo12KMYEVK1CMrf0acItfZoA5Fka9Xh26PNFw00S4LkjKeqH/L6so7R97Fvu7teMu/s1DzsWE1uCNj2p6jpVXYenEKD6PGKONT2Z0iCcyXnhGDpqDvPW7SK7CIlqj7eP4+N5GyMdlokBQROFiJwvIiuBNcBPwFpggstxGWOK2U8rtnPhazMZNWMtAL9t3kuwpukc2/ZluRiZiTYnndmPA92BFaraGM8s7VmuRmWMiZoJi/7gx2XbOPfFqXwYYIitrwtf80ytsu6J0slJojiqqjuBOBGJU9UfgQKdHcaY0mHuut3cNeYXAFZs3V/g+acnLiP9oa8BzwS/xT79Ixt2HXJ8F3I8W5m8fJvj4/3ZsOsgU1duL/LrjTNOEsUeEakATAHeE5F/4xn9FHNsHoUpLWpWTM5dKc/XWa1qeYsD9nX1/fdmecqArNxWMFG8NnkV4On8HjZ6Due9lHdo7dqdBx29R5M/j2foqDkBl349nq30fuYHvly42e/zpz77I1e/NdvR+5VGe7OOssu7wNXR49ls2+tOE6CTRDEIT0f2PcBEYBUw0JVowmSd2aa0mP2XfjxxQdsC+0dek8HUB06nQVVHqxGHbcqK7cxb578cyJHj/tcPu+GdOcxYtYOso8eZsOgPv8f4juzaGuDD7eCRY2zYdYiHP/E/zzeMm5Ggjh3PZm9W8ImI63YeQFXZujeLVdsLJlZ/tmRmhTy6bc2OAzz+1W+oKnPW7uL29+fT7rFv6fT4dwCc9syPdH3ye1dqfTlZuMj37uGdiEdgjCmSiskJ7CuGAoAbdx+icyPPz1lHj+fu7/T37/wev2r7Aa5482eqlU9i54EjjB3enW4nVctzzKfzQxsddeDI8eAH+Th87Dgrt+6nTb0TXxq37ztMWrlEkhKcfD+Gez5cyLiFm3Pnlyz9Yy/Na1XMMyrtrWlrePyr37isSwPGzPH057xwaQc+mL2esTf1YPmWfaSVS+Sqt34mOSGOJZv3cmlGA8bO3cCtfZrwQP+WvP/zeiqnJvLa5FW0rF2R9bsO8vOaE8l5SKd6fDr/xEoM2/YdZly+O6xuT05i615P3a41Ow7k+XtHgpP1KIYATwM18Uy2y5lwF9qCwMaYkE17sC9ZR7P5cuHmAoX+pj14Oht2HyzQ9BNpd41ZQP82tUlOiOfNKatz9wf78N7pbRLJPHTiW/mSzZlkHT3O/R+fuEN49IslXJzRgJTE+ALnCFYpN8fmPYeoW7lc7naLRyYCMOOh03P3d3liEj2bVuPd67rlfti3fHQCWUez+f2Jc0iIj+OIt7RJUkJcng/jxZsyOe+laVzYqT4jhrSl+SMTuKVPEz6a60l4OUkCPItNAazcuo+zX5hSINax3gECr05eRZXUJJ4YvzT3OX/zYXyTBFAgSQC5ScItTibcPQMMVNWlQY80xkRU/SqeJqZ7zyw4iS0tNZHEhPLFEseeg0epVSk+5G/2vo5nKwNe9J/Ubn9/Pv+59kSloG17s/hmyRYe/cL/CgT5S5V8v2wbV3f33PYc9WkSyzx0NE8Cmf77Th75fBEjhrRj1PQ1ZB31HLs36xg79x/mzH95PthXPXlunvPnNCl9Mn9jbpNbTl9NYXLOFYhvkogUN5rjnNyDbS0pScI6s01Zk5JQ8Fu4W44dz+b1nwJ/OPqTM9s7O8An2KSledcT7/rk9wWSxOnPT879ue9zk+n73GT86e/zLf79n9cDMHn5ifN/MNvzjf5v437L87ob3z1R9se3r6H3Mz9w15gFudv+vtHHEnVhPnShiUJEhnibneaKyFgRuTxnn3d/zLHObFPWxBXTLO7lW/axYMOesM4x7fcdAZ9//KvfOFZIBznA6u3BB1vOWLWDVT7HzfE2163flXck1pbMvB3oqlroaK2StrKgEPn/E4GannxHNh0EzvLZVuDTiEdjjAnZ+zd044r//Ozqe1zzdvhDUIeNmhPw+bemreHgkeOMGFJwtFcwj36+mD7Na3DFm3mvw6bdng/5Q/mazPKPZso/gitnyGlJ5MYdRaGJwhYnMqZkOKVpdd68JiNP00ksCeX77Qez1wdMFKc9+yN/H9TG73O9n/mxwL6cUWEjJizLs39hvruj7fvydgZfNrLkFp8o1j4KEXlWRG7ys/8mEXkq8qEYY4rqzFa1oh1CoR79YnGeYbXhWLfzIHPWBF7mNb8VW/cV2PfZL3lHEp3/8vSw4oolbkwtCdSZfTow0s/+N4HzXIjFGFMK/ZGZVejkPH8mBZipDYE7xf05y8/oo2VbCiaP0iLYeulFEShRJKufIiyqmk1od5PFxkY9GRObQvlsX/rH3oDPvxpkWKoToSabksTfHVS4AiWKQyLSLP9O776YHAZgo56MiVEhfC4//90K9+Lw2lOK1wivXiE54ucMlCj+D5ggIkNFpK33MQz42vucMSbGnda8RrRDAOD92eujHUKZ4W+Ge7gKTRSqOgEYDPQFRnsffYALVXV8xCMxxkRMckIca58aQKu6sVFp5+mJy4IfZCKiTb3I/5sHLOGhqouBayP+rsaYYhGTnYnGVdEq4WGMKQG6plfN/TlnpTlbca7ssURhjCnUO9d1ZdK9p+XZl7+cQ+dGVYozJBMFxT2PAgARiY3eMGNMQOWS4kmv5qk2e13PxgD0alY9zzGDO9Qt9rhM8ZqyIvJLwzopMz5dRNYCY4FPVXV3xKOIEBEZCAxs2rRptEMxJioS4uNyF9oB6J5vwSCxtqhS7/CxyMyC9xX0jkJVmwOPAK2BeSLylYhcFfFIIsDmURgTWGpS8ZUlN9ER70JFYUd9FKo6W1XvBboCu7AlUY0pkQZ1qBftEIzL6qaVC35QiJz0UVQSkWtFZAIwA/gDT8IwxpQgFZMTivRts1fT6rx3QzcXIjJuaOjtp4okJ30UC4HPgb+r6syIR2CMKRbf3nsq4JmMd/hY4UX6qpZPyrMeQ1pqImnlEl2Pz8QuJ01PJ6nqPcAiEangdkDGmMgaNawLd57RjDreJomv7+yV5/nG1fOuu12zYt5aQXHWAV6iuDFgwUmiaC0ivwBLgN9EZJ6I+F85xBgTc/q2qMm9ZzbP3U6Iy/tr//ltPfNsvz20S55tpx87V3VvWKT4TGT5KfodNieJYiRwr6o2UtWGwJ/wv06FMaYEyP+FM3+zUp20lDzbZ5xc09EM71Z13Blt6JvkTHQ4SRTlVTV3jUFVnQyUL/xwY0wsC/SFc+UT5+RpuqhWPinqI6Ws4Ss0bjQ9OenMXi0ijwL/9W5fBayOeCTGmKjL/xETSslqdaV4hIkFTu4orgNqAJ96HzW8+4wxJVCgL5xF/TZ6T7/mIReju7tfgXXRTIxyMjN7t6reqaqdvI+7YrWMhy2FakxwgeZSFPZM+aTCGx8aVy/PXf2ahXw/cXc/Z30P57arE+KZnWtVJ/DaDZVTbVgwOJtw11xERorItyLyQ86jOIILlZXwMCa4+lVSeWxgK7/PFXZDkV69PK9f1cn/a3J+cGG0zbLH+9OkRuRG5fdqmrdI4rvX5507fEHHvP0x395zasTeu7i40afjpOnpI+AXPPWe7vd5GGNKqKHe6rL55TQ9vXF15wLP9W8T+Jt9OGkiKcH/R1Gkl/WskW+OSP4P1XL5amFVSrE7CnCWKI6p6mveek/zch6uR2aMiZqTa4ewnGYEvsLOfaRf+CcJolfT6vxtUOs8+/Int3qV89ZJSooveUv25J9AGQlOrsI4EblVROqISNWcR8QjMcaUGA/0b5H7c06eCKflyY1v7g2q5v3Qv6tfs6DvMzhf01OcC5VY3damXuSb3p0kimvxNDXNAOZ5H3MjHokxJmYEGvyUXi2VW/sUXPMl0jOC7z+7RdBjWtct/M5n6gOn59n2F17+fUnxcVbXyo+g8yhU1X9jpjGmzJn58OlUzPetPFITvDo2rMwv6/cAMPzUk7itb/AFyOpXKceSzXsdnb95rYKd4hVT8n4E1qiY7EoJjJIuaKIQkWv87VfVdyMfjjEmltXxs9bBMxe1A/K29z8y4GT+8fXSkM7tm24ePqdlnufeva4rVcsncd5L00I6p6/KqUkF9vnrLA81T3RuVIV562JyxkDEOGl66uLz6A08BpzvYkzGmCgL5SahU8MqAJRPPvG984beJ3HXGUWfUJf/LuXU5jUKbXt/fFBrPrv1lDz72tRz3hmfv5ZU7Xy1rnIkFNJfcWW30l8M0cmEuzt8HjcCnQArN25MKXJR5/p5tovS+nJhp7zn8J3YVyU1kakP9KV/69oAPNg/7x1DOK7ukU5Hb7IC+Nel7Xnv+u6OX39Kk7zril9RyAd/YcmzeoUTQ26b1iydH41FGft1ALB+C2NKkWcvaseqJ88N6xz5Z3xf36txbgK6pkc6Daqm5taDSvezClsofR0D29f1vMbP2NzBHeqRFsKM6vxvG2qPi+/f+/SWNUN8dcngpI9iHCeaH+OBk4EP3QzKGFO8RIT4CI8ELZ+cwHMXt+epIW0LJJGcD+dv7j6VQ0ePA6GNmrqia0PGLdxM/SoF+0x8E86wnumMmr42pLjdqL5aXGpVSg5+UBE4qR77nM/Px4B1qrrRlWiMMaVOQoBJay1qVyzSOXs0qcbIqzvTp0Xgb/DDTmkccqIozGnNazBp6baInMstpzar4cp5nfRR/AQsB9KAqniShTHGRNVZrWsXWvojNM7uIIKVMAGoXcl/R7ibzmlTO/dnt26GnBQFvAGYDQwBLgJmiYiVGTfGhCxQ61L+OkzFpUAfhXfb3yinq7s3AuC+s5rTrv6JUVjf3nMqP//5DK49JZ3Xr+pEI28fTIcGlXOPGTWsS+5Q4kjyN/kx0pyk4/uBjqo6VFWvBToDD7ob1gki0kdEporI6yLSp7je1xjjpoIfws9c1N7RK9vXd16iwt9iSu/f2C3Pdly+TJGzdWmXBnn2x8dBTu4on5zAw+ecTHJCHG3qptG8VkVqVUohPk7o36ZO7jmqVzgxd6Nvi5p5tiPhk1t60LZ+Go8PbhPR8+bnJFHsBPb5bO/z7gtKRN4WkW0isjjf/v4islxEfheRh4KcRoH9QApgfSPGlGCBuqtzSmf0bVF4O/vyf/Tnk1tOKfT5Au/noH+8Xf75GYW035zeslaeju4eTaqx/B/n+B1h1bqu55w5dxTdT/KUx8s/qz0cL13ekc6NPOfNSZ498g31jRQnndm/Az+LyBd4/p0HAb+KyL0AqvrPAK8dDbwM5M7iFpF44BXgTDwf/HNE5Es8I6pG5Hv9dcBUVf1JRGoB/wSudBCzMSaGFdaWPuvhMwIuFpScUHjZ8VFDu/DjcgedzfmSR/7Cfz28H+oD2ubtk4iPE8cjs569uB3DeqZz4IhnRFfOOhjp1SJX2dW3zlW7+pWZ90g/qlWI3qinVd5Hji+8fwYdrqCqU0QkPd/ursDvqroaQETGAINUdQRwXoDT7Qai04hpjImIYJ+zhc2KdqJvy5r0jcA8hqY1K7L2qQG526lJ8Rz0fuDnCNZnnJqUQEa6J+H87/puud/0a1RMZu1TA7j49RnMWRvZsh9uJQlwVhTwbxF+z3rABp/tjUC3Qo5FRIYAZwOV8dydFHbccGA4QMOGpX9KvTElWXHNVIhkeT8p4vl6Nase/CAHHhvYisfG/XYinmKc7+Fkwl0N4AGgNZ5+AgBU9fRCXxRBqvop8KmD40YCIwEyMjKs/KMxQfzzkvacHGTNaOP/LsitD+lXrujEAx8vzG2y8tU4gkvChspJZ/Z7wDI8ZTv+BqwF5oTxnpsA3+EE9b37jDHFaEin+qU+UfjrUyjqt0jf3BBuKfKrvMNs8xvQrg4zHjqjwMgsgBa1ijY5MRKc9FFUU9W3ROQu7+S7n0QknEQxB2gmIo3xJIjLgCvCOF8uERkIDGza1P1xxcaYoijem32n7/bS5R2pkOz/4/B/N3RlzOwNlIvg+t2DOtQjI70qPZ/6ocBzaamJBZZk9e0zyVGchUacJIqj3j//EJEBwGY8M7SDEpEPgD5AdRHZCPzVm3RuB77BM9LpbVVdEnLkfqjqOGBcRkbGjZE4nzHGHbFQT6lTw8q5P+cUGfSnc6OqucNQc0Qi/vzJIFTFOUHRSaL4h4ikAX8CXgIqAfc4ObmqXl7I/vHAeKdBGmNKh+JePK6On1FUOTH4W7SopPB3h+EmJ4lipqpmAplAX5fjCYs1PRkTWYG+ODerWSHkb7UNqnpKWwSaKxFJqUmFf8QV5aaguBKdbxI7rbk7hf5C4SRRzBKRBcAoYILG8IKy1vRkTGQF+m3/7t7TQj7fw+e2pGfT6nRJd9R67Qp/ZT1C5XbLWa1KKYwa1oXOjapQKYKzuYvKyain5niGnV4NrBSRJ0WkeZDXGGPKoJTEwB8pyQnxnNmqVjFFE5i/RY+CyZnEXRw9LH1b1IyJJAHOJtwp8B3wnYj0Bf4H3Oa9y3hIVWe6G6IxJlpC/eb80/192ZKZ5U4wEdKpYRVa1q5YpOVY7z2zBUezlYszGgQ/2AWJ8cLR48XfqONkwl014Co8dxRbgTuAL4EOwEfE0LKo1kdhTHTVqpRCrSisyRCK8skJTLz71CK9Ni01kScvaBvhiJyb85d+ZB3NLvb3ddL0NBPPSKfBqjrAO1P6dlWdC7zuanQhUtVxqjo8Lc15GWJjTOmXFGCVvWi6sXdo37MrpyaFVQ+rqJx0Zrfw04F9L/CCqj7tQkzGmCjLGXXTsohLlcaSKff3pXxybA6F/cuAVvxlQCvSH/o62qEE5LSPIr/oz5YxxrimRsVkxgzvTtv8azWUQA29q83Fsqu6N6Ra+dgtju3kjsKfmBwia30UxkRO95PcWQTHFPSPwdHr93Ci0IY7EdknInv9PPYBhc93jyLrozDGmMgr9I5CVUt+46QxxpiwxeZQAGOMMTHDEoUxxpiASlWiEJGBIjIyMzMz2qEYY0ypUaoShXVmG2NM5JWqRGGMMSbyLFEYY4wJSGJ4eYkiE5HtwLoivjwNzyJNbr422HGFPR/K/vz78m9XB3YEjbToYvk6FvZcLF7HwuKK5OvsOkbuddH+3Q73OjZS1YIrJamqPXwewEi3XxvsuMKeD2V//n1+tueW1evo9JrFwnUM51radSze6xjOtYzU77Zb19GangoaVwyvDXZcYc+Hsj//vnD+XkURy9exsOdi8TqG8552HSPznqG8rlT+bpfKpicTnIjMVdWMaMdR0tl1jAy7jpHh1nW0O4qya2S0Aygl7DpGhl3HyHDlOtodhTHGmIDsjsIYY0xAliiMMcYEZInCGGNMQJYoDCIyWETeFJGxInJWtOMpyUTkZBF5XUQ+FpFboh1PSSYi5UVkroicF+1YSioR6SMiU73/J/sU9TyWKEopEXlbRLaJyOJ8+/uLyHIR+V1EHgJQ1c9V9UbgZuDSaMQby0K8lktV9WbgEqBnNOKNVaFcR68HgQ+LN8rYF+J1VGA/kAJsLOp7WqIovUYD/X13iEg88ApwDtAKuFxEWvkc8oj3eZPXaEK4liJyPvA1ML54w4x5o3F4HUXkTOA3YFtxB1kCjMb5/8epqnoOnqT7t6K+oSWKUkpVpwC78u3uCvyuqqtV9QgwBhgkHk8DE1R1fnHHGutCuZbe47/0/nJeWbyRxrYQr2MfoDtwBXCjiNhnlVco11FVs73P7waSi/qeha6ZbUqlesAGn+2NQDfgDqAfkCYiTVX19WgEV8L4vZbeduAheH4p7Y4iOL/XUVVvBxCRocAOnw88419h/x+HAGcDlYGXi3pySxQGVX0ReDHacZQGqjoZmBzlMEoNVR0d7RhKMlX9FPg03PPY7VzZsglo4LNd37vPhM6uZWTYdYwMV6+jJYqyZQ7QTEQai0gScBnwZZRjKqnsWkaGXcfIcPU6WqIopUTkA2Am0EJENorI9ap6DLgd+AZYCnyoqkuiGWdJYNcyMuw6RkY0rqMVBTTGGBOQ3VEYY4wJyBKFMcaYgCxRGGOMCcgShTHGmIAsURhjjAnIEoUxxpiALFEY44eIVBORBd7HFhHZ5P15v4i86tJ73i0i14T4mudE5HQ34jEmh82jMCYIEXkM2K+qz7n4HgnAfKCTd/KUk9fE4ynV8Kaq2oJTxjV2R2FMCLwrhn3l/fkxEXnHu4LYOhEZIiLPiMgiEZkoIone4zqLyE8iMk9EvhGROn5OfTowX1WPiUgTEZnv857NcrZFZK2IPO3dvlhV1wHVRKS2+397U1ZZojAmPE3wfMifD/wP+FFV2wKHgAHeZPEScJGqdgbeBp7wc56ewDwAVV0FZIpIB+9zw4BRPsfuVNVOqjrGuz0fW03PuMjKjBsTngmqelREFgHxwETv/kVAOtACaAN8JyJ4j/nDz3nq4KnRk+M/wDARuRfP8rRdfZ4bm++124C64f01jCmcJQpjwnMYQFWzReSonuj0y8bz+yXAElXtEeQ8h/Csa5zjE+CvwA/APFXd6fPcgXyvTfG+3hhXWNOTMe5aDtQQkR4AIpIoIq39HLcUaJqzoapZeCqBvkbeZid/mgOLIxOuMQVZojDGRd71iy8CnhaRhcAC4BQ/h04ATs237z08dybfFnZ+bx9IU2BuJOI1xh8bHmtMjBCRz4AHVHWld/s+IE1VHw3wmgvwDKkt9BhjwmV9FMbEjofwdGqv9CaNnBFVgSQAz7sdmCnb7I7CGGNMQNZHYYwxJiBLFMYYYwKyRGGMMSYgSxTGGGMCskRhjDEmIEsUxhhjAvp/1RJFIvs6s84AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Make some plots\n",
"# First, the LCN vs time.\n",
"plt.plot(timeyr, lyapunov_CN)\n",
"plt.xlabel('Time (yr)')\n",
"plt.ylabel('Lyapunov Characteristic Number (2pi/yr)')\n",
"plt.xscale('log')\n",
"plt.yscale('log')"
]
},
{
"cell_type": "code",
"execution_count": 57,
"id": "abandoned-register",
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAFECAYAAADIlyJZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7S0lEQVR4nO3dd3hUddbA8e9JKKGGLp0gIIgIAgFEwF5wFevada3YVtFVV33Xrrv23lZZu2KvICgoKIiiQOgdpHcCISQhPef9YyZhkkxmbpK5U5LzeZ55ktsP13jP3F8VVcUYY4wBiIt0AMYYY6KHJQVjjDElLCkYY4wpYUnBGGNMCUsKxhhjSlhSMMYYU8KSgjHGmBJ1KtogIgMcHJ+vqotDGI8xxpgIkoo6r4lIBjAHkADHd1XVJBfiMsYYEwEVvikAc1T1+EAHi8i0EMdjjDEmgip8UyjZQaSTqm4KUzzGGGMiyElF8yTXozDGGBMVnCSFeSIyyPVIjDHGRJyT4qMVQHdgA5CFp+JZVbWv++EZY4wJJydJoYu/9aq6wZWIjDHGREyg1kfAgYe/iLQBElyPyBhjTMQErVMQkTNEZDWwDpgOrAe+czkuY4wxEeCkovkR4Ehglap2BU4Afnc1KmOMMRHhJCnkq+puIE5E4lT1JyDZ5biMMcZEQNA6BWCviDQGfgHGichOPK2QjDHG1DBOWh81ArLxvFVcAiQC47xvD8YYY2qQoEkBSpql9lDVH0WkIRCvqhmuR2eMMSasnLQ+Gg18DrzuXdUB+NrFmIwxxkSIk4rmvwPDgH0AqroaaONmUMYYYyLDSVLIVdW84gURqQMEL3MyxhgTc5wkheki8i+ggYicBHwGTHA3LGOMMZHgpPVRHHA1cDKewfAmA2+okxpqY4wxMcVJUhgFTFTVovCEZIwxJlKcFB9dAKwWkSdFpJfbARljjIkcp/0UmgIXAVfiqWR+G/jI+ioYY0zN4uRNAVXdh6evwsdAO+BsPDOy3exibMYYY8LMSZ3CGXjeELoD7wHvqupOb8/mZaqa5HqUxhhjwsLJgHjnAs+p6gzflaq6X0SudicsY4wxkeCoTsEYY0ztUGGdgoh8G+xgJ/sYY4yJHRW+KYjIXmCG343eXYDDVPVgF+IyxhgTAYHqFM50cHxe8F2MMcbECqtTMMYYU8JRPwVjjDG1gyUFY4wxJSwpGGOMKRGw85qIDAUuBUbgGd4iG1gCTAQ+UNV01yM0xhgTNoGapH4HbAW+AeYCO4EE4BDgOGAU8Kyqjg9PqMYYY9wWKCm0UtXUgAc72McYY0zsqLBOQVVTRSReRH4KtI87YRljjImEgBXNqloIFIlIYpjiMcYYE0FORknNBBaLyA9AVvFKVR3jWlTGGGMiwklS+NL7McYYU8M5nY6zAdBZVVe6H5IxxphICdp5TURGAQuA773LR4iINUM1xpgayEmP5geBwcBeAFVdANhw2cYYUwM5SQr5fnouF7kRjDHGmMhyUtG8VEQuBuJFpAcwBvjN3bCcadWqlSYlJUU6DGOMiSkpKSmpqtra3zYnSeFm4B4gF/gImAw8Errwqi4pKYm5c+dGOgxjjIkpIrKhom1Bk4Kq7gfuEZEnPIuaEcrgjDHGRA8nrY8GichiYBGeTmwLRWSg+6EZY4wJNyfFR28CN6rqLwAiMhx4G+jrZmDGGGPCz0nro8LihACgqjOBAvdCCk5ERonI2PR0m87BGGNCqcKkICIDRGQAMF1EXheRY0XkGBF5Ffg5bBH6oaoTVPXaxEQbp88YY0IpUPHRM2WWH/D5PfjYGMYYY2JOhUlBVY8LZyDGGGMiL2hFs4g0A/4GJPnub0NnG1MzrN2VScvG9UlsUDfSodR4X83fzOY92dx8Qo9Ih1IhJ62PJgG/A4ux4S2MqXGOf2Y6XVs14qc7jo10KNWSnp1PUZHSvFG9SIdSoX98shAg5pNCgqre5nokxpiIWZeaFXynKNfvoSkArH/8tAhHUt7kpdu57v2USIfhiJMmqe+LyGgRaSciLYo/rkdmjDE1xJsz11Vq/7yCIvZk5bkUTWBOkkIe8BQwC0jxfmzAIWNqsey8QpZv2xfpMIL611eLmb1uT6WPKyxSnv1hFXv3R+bBfP0HKQx45IeIXNtJUrgd6K6qSara1fux+RSMqcX+8ckCTn3hF/bl5Ec6lIA+/GMj578+q9LH/bxyJy9OXc393yx1Iargpq3YGZHrgrOksAbY73YgxpjYMXeD59t3Tn5hhfvMXrcHJ9P9RqP8Qk/cgf59lSEBti3fti9iRUX+OKlozgIWiMhPeIbPBqxJqjG1xYPjl7JkSzqf33CUz9pAjzmYsnQ7176fwsNnHsbfhia5Gp+bwpHSTn3BM4rQyb0PYuzfksNwxcCcJIWvvZ+o4Z03elT37t0jHYoxNd47v62veGMFT83NadkArN0VuVZN1XlLkcA5zxVTlu0I/0X9cDKfwrvhCKQyVHUCMCE5OXl0pGMxpjYK9tCMxEO1Jkq6eyK/3HkcnVo0DNs1ncynsE5E1pb9hCM4Y0x0C/ZdPNrqFOZvTOOHCr6R78rI5cHxS8kvPNBHNxTh788r4I8qtIAqNn3VruoHUQlOio98C7kSgPMA66dgTC0W7EUgWl8Uzn7VM728vw5uD4xfwqTF2xnarWXI4k/fn0+/h6f43ZZXUER6dunWW0l3T+Tbm4eXWhfutOqk+Gh3mVXPi0gKcL87IRljQmn2uj0k1I2jb8dmkQ4lrCr7Lb/A2+JIVZEQlX+t3ul/9uL5G9N4cepqflpZ/i1g/MKtIbl2VTkZEG+Az2IcnjcHJ28YxpgoUNxO343hH6KsdKhawlkPUvzG4k9ahJunOnm4+86rUACsB853JRpjTI0SizlDlSqVf63ZmUHH5g1JqBtfret/lrLZT0DhE7SiWVWP8/mcpKqjVXVlOIIzxsSmYMUv2XmF5BVUbdDl9alZPDppebUrsbPzCkm6eyLfLNgCwOSl/iqgA19DVVm2dR/p2fmc+OwM/vn5opJtSXdP5KVpa6oVo78ICovU1Qr8QNNx/i3Qx7WIjDFRz2lRS0XPrkPv/57TXvzF/8Ygrv8ghbEz1rJmZ2bgawc5z/Z9OQA898OqUutzCgod1ymM+2Mjf3nxF6Ys3Q7AhDL1AaFoObQlLbskCWTlFtDtX5NCkmwqEuhNYVAFn0eAt1yLyBgTM7SCR2+gZ2pRkeeY1QEe6iOfn8GIJ6f5P977gCyq5pflir5tP//jap99YOe+HE56djqb08qP9lM8KOCG3e6NBPT6jLV8OncTAGneAfo+mbPJtesFmo7z5uLfxZM2LwHuwjPhzn9ci8gYE/XEYaG7v6Tx4eyNJb+nZ+f7nfFtxXb/rXZ8r10UpAjl6/lbgsTmsb7MAz0jp6DkX5eamcunczexemcm7/++gdEjDqaoSGnTNIGUDXtKWiylZpaMAMSK7ftCPovdXV8splOLhlz8vz9Cel5/AlY0i0gd4ArgDjzJ4K9Wn2CMKVbRczlQytiZceABmltQCFTuAVr8FlLRtYublN7+2cIKz5FXUBSw/vYrb0JZuDmdhZvTAU/RTfK/fwRg0pgRnPvfA6OvfuzzzX3k81UrFgvm+jBN0hOoTuHvwDJgIDBSVa+whGCMgeo13yx1aJAioP15BSzZkk5mbkHJujgJ/KbgpA420FvGnqw8Ji7e5ieWAyOm7s7KLbfdbftyDtyDiFQ0Ay8BTYHhwHgRWeT9LBaRRQGOM8bUEk4fTamZuSz2fuOuzPF/7szi9JdmlvqW7O9NYe76A8NIxGIz2Mramp7j2rkDFR91de2qxpiY5vttX1V59ec/eWrySi4a3JnHzjncZ5vn56kv/MKujNxyHeiCfeH9cr6nzf78jWkl65Zu9VTu+n7bv9OnKag67mhQG9JH5VX4pqCqGwJ93AhGRBqJyFwROd2N8xtT2302dxNJd08kfX/1Zkwr/qaqqizcnM5Tkz0lyx/N3sjuzNySr/Prd3uGzt7lU4/gW/RUUeulYm//ur7CbV/N31IyXeba1ANDdDt91Fe2BCbbp/hoi3do8JrIycxrVSYib4nIThFZUmb9SBFZKSJrRORun013AZ+6GZMxtVnxQ3aTn+aVTpVtmllYVLoT2sB//8gr3nb0v64pO3Sa85ZLZc1YtavUfMvv/LaeYY/7b7bqz7rU0nM7FM+u5pTvFJl3f7m4UsfGErfHMHoHeBl4r3iFiMQDrwAnAZuBOSIyHuiAp2I7weWYjDFVMHvdHgZ3bUFO/oEkMGHhNj6Zs7HcvsUdw4Ip29egog5pCvztrdnl1mfllZ8us6I3gL+8cKBVUHZeIYWV7OiQW8Ue2LHGyXwKjUQkzmc5TkQczfigqjOAsgOJDwbWqOpaVc0DPgbOBI4FjgQuBkb7XtMYE1778wrYn1dQal3xwHq+nvh+Rbl2/sH4Fh/d/OG8UttOfHY6Jz47vVLnK6u4yMrXqh0ZZPvMtzzPp47ClObkwTsV8E0CDYEfq3HNDoBvd7zNQAdVvUdVbwU+BP6nqn7Tsohc6613mLtrV3gnnzAm1jltStr7/sn0fdD/PADVHU10xfZ9Jb/P27iXsTP+DHrMfj9vBBX5ouyAcnj6GPiavHR70PqM2spJUkhQ1ZJ3Ou/vrs0Np6rvqOq3AbaPVdVkVU1u3bq1W2EYU+sVFCnv+pmfubojTE9avL3U8qOTVlTrfP/6qnT5vr9HfdmxjLLyCjknwPDVtZmTpJDlO6eCiAwEqlP1vgXo5LPc0bvOGOOyyra4eWD80pBev6Aw9OXyH/5Ruk5j7Iy1QesLtu7NpqC6gydF2LeL3JmMx0lSuBX4TER+EZGZwCfATdW45hygh4h0FZF6wIXA+MqcQERGicjY9PTynWGMMc7tzKh8J6jqzEp2XZiGanhoQulkVjbi3PzSySna5pJ2YubqVFfO62Q+hTlAL+AG4HrgUFV19F9WRD4CZgE9RWSziFytqgV4kspkYDnwqapW6uuIqk5Q1WsTExMrc5gxtc60Ff4nqS82+D9Ty61buGlvwGOqU3w01adZp6/q9psoa08lZy8LNGJrbVNhk1QROV5Vp4nIOWU2HSIiqOqXwU6uqhdVsH4SMKlyoRpjKuuqd+aWWnbyJf/MV351KZqK9Xt4Cl/eeFTIzvftotJjF71epjJ72bZ9pZZPfm5GyK4dLm693ATqp3AMMA0Y5S8eIGhScIuIjAJGde/ePVIhGFMrPTppOZcM6VypY57/cVXwnYBHvl1WlZAcKVu5bSoWaD6FB7y/Pqyq63y3iUhEx0VS1QnAhOTk5NGRjMOY2mbsjLWM6tu+Usf4Tlpjop+TiuYv/Kz7PNSBGGPCpzpFD26174/But6IWpvqTj1IoPkUeonIuUCiiJzj87kCG4rCmJhUXKcw6uWZpeYoAPjvz8E7kUHVxy4KxnJC5WTkFATfqQoCvSn0BE4HmuGpVyj+DAAiWmxjTVKNqbyCwiK2+4zDv2xr6cpWJz2LAZZsdef/u2Ctnkx4BKpT+Ab4RkSGqmr5QU8iyOoUjCnvz12ZTFi4lVtO6OG3L8FDE5aRmnmgqeazP5SeSLFInbXX/78aPEJoLHGruM1JncLZItJUROqKyFQR2SUil7oTjjGmqi594w+e/3E1S7Z43gAenlC6Nc+Py0v3Wdixr/SUkqpq5foxZEcVOh464SQpnKyq+/AUJa0HugP/dCUaY0yVFQ/tPOrlmQC89WupRoPlHvj+5hfIzHOnnNqEXmWH/nbKSVKo6/15GvCZqka8IN/qFIwJLNvfPANBqnKz8wsrHBnVRJ9IVDQXmyAiK4CBwFQRaQ24N2u0AzbMhTHl+dYHPDppuZ/t4YzGxConYx/dDRwFJKtqPpCFZ1IcY0yUqspAd8ZAJcc+KtOiIWLDXBhjoN9DUzi1T1seP7cvqZm5pPkMKuevuNleFIwTgd4UjvH+HOXnc7rLcRljKlBUpKRl5ZGenc/HczyTGK7eUbp369z1ZWfBteKjmiahrjszFgcc+8g7T/J3qvqpK1evIhsQz9RmL0xdzQtTA48nlOZnKOrUzFw/e5pY1bBeoPFMqy5gqvHOk3ynK1euBqtoNrXZ5KXlR/y0+YZrH7cmBnLy/vGjiNwhIp1EpEXxx5VojDFVYzmh1nHrP7mTpHAB8HdgBpDi/cwNeIQxJmyS7p7Igs17Ix0GAIOSmjvet2Wjei5GUvPtDfFsdcWcNEnt6udzsCvRGGOq5MnvVwbfKQw+vW6o430HdHGeQEz4BE0KItJQRO4VkbHe5R4iYq2PjDHl+BuIryJ/Py5wQ5EmCdWvSG0agnPUNk6Kj94G8vB0YAPYAvzbtYgcsGEujIl9zRvWDbj9iqOSqn2Nlo3rV/sctY2TpNBNVZ8E8gFUdT+4NMuGQ9b6yJjY16Vlo4DbK/uQufn48m8ewRJPLEtq2dCV8zpJCnki0gBvZbeIdAOswbMxEbBkS7prA6FV1+VDu4T0fCP7tKvU/pcMKX/9Q9s1LbV8zfCITi8fUpUpqqsMJ0nhQeB7oJOIjAOmAne5Eo0xpkKFRcrpL81ky97sSIfiV3ETycq0QKrI0odOoXf7pn63PX/BESQ2KP8G0DYxgW9vHl5q3TUjSreJ6dm2Sanlvw7sWM1IIydi/RRUdQpwDnAF8BGegfF+ciUaY0yFuv1rkuvXuP6Ybnxw9ZBqneOukb0Cbj/mkNbMuefEgPsUVfDA696mMWf178Dn1/tv5VS2crpsEUtcmW/Xzfwkl9rOSeujqaq6W1Unquq3qpoqIlPDEZwxJrx6tGnM8B6tqnWO5KTAfVvfvWowrZsErgBuUDfe7/r6dTyPrB4HNeH+03uXP65e6ePKFrHEx0mpCuy6dUo/AnuVeZOIZmHvvCYiCd6ey61EpLlPb+YkoINL8RhjosCH11T+baEqpRm/3X08L13Uv9S69Y+fRp14z6Np6u3H8P2tI/wee+WwpJLf/3lKTwDaNEkoWVdcxPTh6AP/FhG49cQeJctXDevKQG9/ic+uH8qHo48s2ea734hqJko3tEtMCL5TFQR6U7gOT+/lXhzoyZwCfAO87Eo0xtRSGTn5bNidxbQVO7ji7dmltuUVFLEnKy8scfTv3AzwfKOujm/+Pqxc+T7A7H+dUGq5fbMGjOrXnpl3Hef3PN1aN6ZX2wN1C5ceeaAy2fctoHOL8i1xfrv7eACO6taKFt7e00O6tqRZw3osf3gk8+47idZN6pe0clKFxvU9xU914oQxx3uSQvc2jTmjX/ty52/kfSt5/bKBfmP3ddYR5Y+vrrrx4R8l9QXgBRG5WVVfcuXqVWSjpJqa5rzXZrFiewYinodTfmFRyf/0t3w8n++WlB8EL+QxDOzIwa0bAzAoSBGQrz4dmrJkyz6a+TT/7NepGQBT/nE0Jz83A4CTex9Em6b+v912bO6seWXZb+wr/z2Sn1fu4uTeB5Ws+/Xu42lYN55G9Q883ubdd1Kp4xrUiy8pahrarSVzN6TRpkl96tWJ46sbj6Jbm8bExQmr/n0qcQKrd5YemhzgqfP68fcP53FUt5Z8f+sICgqV//tyMYu3pHv/TQ3YnJZNv07NeP7C/gzs0pyNe/bzy+pU0rPzySsoYrdPsh/YpTkpG9Ic3Qc3OenuVyQizVR1L4CINAcuUtVXXY0sAFWdAExITk4eHakYjAmlFdszAE9FaKEquQUHkkI4EgJA19YH+g3ExQnDurfk1zW7efGi/oz5aL7fYx4+8zAuHdKFcX9s4PxBncpt79Gmccnv1WlB2aFZA7+trurXieeUw9qW27cybj3xEP46sGNJv4n+nQ+0nqrnrXM4tF1Txt80jO5tGpOVW8ianZkM7daSdY+dBkCvtp6EOKHM29G61KyS+pPLhiaVu3ZxC6L9eYUlSWx7eg63f7aAVy8eyKqdGTw7ZRWz1u4GPEVrN45LYdLi7X6b4IaCk6QwWlVfKV5Q1TQRGQ1ELCkYU1Pc/unCUsvxcUJhkfLqT2uYuHgbX984LGyxXHd0t1LL4i1YqagD2CNnHlbyoPP3wANPEU/zhnVJ259Pu8TAD+v2iQmc5PON31eLRvXYsjebOnGhLzKJj5OgHekA+nZsBnjmMQhWUV6sa6sgHfS8mdL3raZtYgLjrvHUbQxKasFH1x7Jvpx8du7zdA97+aIBrDguo8Imu9XlJCnEi4ioN6WJSDxgwxsaEwJfzNvsd/2rP/8JQP9HfnD1+g3rxbM/r5CeBzUpV49wSp+2zFyTSlLLRrx31WCy8wu57v2Uku3nJZd/M/Bn3n0nMX7h1nLf6Mv67f9OqHDbm5cn8+PynbR1qXI12jVNqEvTBE9yjosT1xICOEsK3wOfiMjr3uXrvOuMMSGWV1AU1utNGjOCY5/+uVynLoBLh3Tm3AEdaFivDp28FblNE+qwz9ujOqGCZqNliQhnHlG9BottmiZw8ZDO1TqHccZJUrgLTyK4wbv8A/CGaxEZY8LihF5tSGrViM+vH8ph7cuPIyYi5aZ8/PH2Yxj8H+umVJMFTQreKTn/6/0YY2JcfJzw7pWD6dfJkwiCdTbz5dsPwNRMQZOCiPQAHgN6AyV/ETbRjjGxaelDpzgu+qnI4EokEhNbnBQfvQ08ADwHHAdcibOB9Iyp9XILCpm/cS+j35vLz3ccy9rULPILiqhfN46ebd2rLAykuglh4f0nk1DPHgE1lZOk0EBVp3pbIG0AHhSRFOB+l2MzJuYNfWxaSW/kgf/+sVrnalQvniUPnYKI8MyUlbw0bY3f/Vb/51TOfPlX2iUm8OAZhzHiydCOX5lYg+coMM6SQq6IxAGrReQmPDOvNQ5yjDEGQjY8xZ0je3LVsK4l7dr7edvMF2tSvw4ZuQVcOKgTdePjmHSLZ7wg3+GVX764P4d3sImpTGBOksItQENgDPAIniKky90MKhgb5sJEs2enrGTNrkxevST4mDhO3Xhs6b/1+PgDfQouGtyJx87p6/c43/GBTu8b+vF3TM0TMCl4O6pdoKp3AJl46hMizoa5MNHsRW+xzobdWSE5331+hoju4u03cNfIXtxwbLdy233dObInR/doHZJYTM0nwWbvEZHfVfXIgDtFSHJyss6dOzfSYRhTStLdE6t9jm9vHk7zRvUCjuOzLT2btk0TXJuW0dRcIpKiqsn+tjkpPpovIuOBz4CSrz6q+mWI4jPGlNHHQdl/sLGEjKkKJ0khAdgNHO+zTgFLCsYYU8M46dEcFfUIxtQG464ZUmpeAmPCzUmP5gTgauAwSvdovsrFuIyJSkVFisiBVj17svJo0aged36+kEMOasL6KlYuXzioE1cP70qPg2JnjmBTMznplvg+0BY4BZgOdAQy3AzKmGh18L8mcftnnjkQVu3IYMAjP/DJnI18Oncz/564nA9+31ip8910nKep6TUjDraEYKKCk6TQXVXvA7JU9V3gNKDys3obU0N8OW8LABt27wfg/d83VPlc/zjpEH687Ri6t7H+oCY6OEkK+d6fe0WkD5AItHEvJGOik+9cB+MXbiVtv6e38rKt+yp1noN9ZuOKjxNLCCaqOGl9NNY7L/N9wHg8Q1zYuEem1tmfV1Dyu++cxUWBu/qU0+Ogxtx4XHd2ZuSEKjRjQsZJ66PiCXWmAzZctqm19ucVhuQ8Z/fvwMg+7UJyLmNCzUnro/rAuUCS7/6q+rB7YRkTeYVFysvT1nDZ0C60aFSv1JtCZdx/em8uPyqJ/MKiag9bbYzbnBQffQOkAylArrvhGBM9Zq5J5bkfV7E2NZPJS7dzQq+DKnV8h2YNePq8fgzu2oL4OCE+zhKCiX5OkkJHVR3peiTGRJmcfE9x0da92eTkFzFx8bZKHd+iUT2GdmvpRmjGuMZJ66PfRORw1yMxJsoUFHpqkHfsq9oL8iHW78DEoArfFERkMZ4xjuoAV4rIWjzFRwKoqvofwL2KRORQPHM3tAKmqup/Q3l+YyorPdvTGnvHPuethO75y6Gcl9yRSYu3c1Z/m7/AxJ5AxUenV/fkIvKW9zw7VbWPz/qRwAtAPPCGqj6uqsuB672zvL0HWFIwEaGqFBRpST+EXJ/+Cf70bteUDbuzyMorpGmDOjRrWI+Lh3QOR6jGhFyg4qM2QG9V3eD7AXrj+TbvxDtAqfoI78Q9rwCnes91kYj09m47A5gITKrUv8KYEHpx6hp63POd4zeESbeM4OrhXQFo3rCem6EZ47pASeEJYJmf9cuAp5ycXFVnAHvKrB4MrFHVtaqaB3wMnOndf7yqngpc4uT8xoTSvpx88gqKeO7HVQAs3pJe4b79OzcrtXzT8T14/oIjOKl35VooGRNtAhUfNfG+GZSiqhtExOmbgj8dgE0+y5uBISJyLHAOUJ8Abwoici1wLUDnzvaKbkKjoLCIvg9O4YReB0ZwWbrF//AV9552KBcN7sxhD0wuWVevThxn9e/gepzGuC1QUmgeYFvDUAeiqj8DPzvYbywwFjzTcYY6DlM7pWxIA2Dqip0l6/IKD9QlnHVEe75esBXwjGgKcE7/DjROcNKq25jYEegv+kcR+Q9wr3onchbPIPIPAdOqcc0tQCef5Y7edY6JyChgVPfu3asRhqntdmfmsmzbPkb0aM2UZTsC7tunQyI9DmrC6h0HRo1/9oIjXI7QmPALlBRuB94A1ojIAu+6fsBc4JpqXHMO0ENEuuJJBhcCF1fmBKo6AZiQnJw8uhpxmFruqckr+XjOJsZdM4Rf16SWrI8TaFA3niyfsY6SWjbiRKsvMLVAhUlBVbPwtAw6GM+sawBLVXWt05OLyEfAsUArEdkMPKCqb4rITcBkPE1S31LVpVX9BxhTVXuyPE1O7/16CVvSskvWt22awNZ0T8uj0/u249tF2+jZ1jqimdrBySipawHHiaDMsRdVsH4S1uzURMDU5Tu464tF/HjbMaRmenoqr0v1TKE5KKk5c9an0b5ZA/5zzuFsTsvmokGduGtkLzq1CHk1mjFRyckwF1FHREaJyNj09IqbDBpTrKCwqGQcozdnriM1M48v5m1hc1o2/Tomlux3dI/WALRv1oDjerbhsiO7UCc+zhKCqVViMimo6gRVvTYxMTH4zqbW+/fE5Zz1yq+oKvXqeP7k35q5jp0ZuRzXqw2JDeoiAsN6eFpat2/WIJLhGhNRQZOCiLwoIkeFIxhj3DBn/R5WbM9g5Y6MkqKiLXs9dQhdWjbkyINb0L11Y3q0aUzj+nXo3b5pJMM1JqKcNLJOAe4VkZ7AV8DHqjrX3bACsyapJpgte7Npn5hAkcLqnZkATFy0jU179nPVsK689es6ADo2b8hj5/QlO7+QJgl1+eNfJ9Cwns17YGqvoG8Kqvquqv4FGASsBJ4QkdWuRxY4Jis+MhVan5rFiCemMWnxdjbsziKvoAgReH3GWooUDu/YlPreYqQOzRrQolE9OniLjBrVr4OnO44xtVNl6hS6A72ALsAKd8IxpvpSNqRRpDB91U5WeTubnT+wE3ne0U4PbtWYiWOG8/fjutEuMSGSoRoTdZzUKTzpfTN4GFgCJKvqKNcjM6aKigeym7V2Nyu3ZyIC94/qTesm9QE4uHUjurdpwj9P6WVvBcaU4aRO4U9gqKqmBt3TmCiwaPNeADbtyWbaih10btGQRvXr8PMdx7JqRwZNEupGNkBjopiTOoXXgaNE5GnvJ+JvCdZPwZT1565M8guLKCgsYunWfQzr7pkbeeHmdEZ4m5o2ql+H/p0DjfNojHFSfPQYnmkyl3k/Y0TkUbcDC8Qqmo2vdalZnPzcDB6btILVOzPJLSjirwM70r1NY07t05b7Tu8d6RCNiRlOio9OA45Q1SIAEXkXmA/8y83AjHHqrZnrKCxS3pu1vmS2tCM6NWfyrUcTH2d1BsZUhtPWR818frev5yaiMnMLSkY1TcvK47OUTZx46EHUrxPHxMXbuOKoJLq2amQJwZgqcPKm8BgwX0R+AgQ4Grjb1aiMCeD5H1bxxsx1fHvzcKav2kVOfhH/PKUnG3Z3ZFdmLhcPthn5jKkqJ6OkfiQiP+PpvAZwl6pudzWqIKxHc+2Vk1/IF/M2AzB2xlpmrd3N0Ye0pmfbJja8tTEh4LT4KA5IBfYCh4jI0a5F5IBVNNdek5duJ21/Pr3bNWX8wq3sysjlmuFdIx2WMTVG0DcFEXkCuABYChRPWqvADBfjMsavj2dvonOLhrxw4RGc9NwMeh7UpKTJqTGm+pzUKZwF9FTVXJdjMSag7xZvY9ba3fzfqb3ocVATHhzVm8M6JFqvZGNCyElSWAvUBSwpmLDKKyjig9838NufuxnctTkvTV3DEZ2aceUwT3HRFcOs2MiYUHOSFPYDC0RkKj6JQVXHuBaVqfW2pWdz0djfWb97P62b1OfH5TtomlCHly7qXzJRjjEm9JwkhfHeT9Sw1kc135u/rGNzWjZvXzGIY3u2ZtWOTOrVsakxjXGbkyap74YjkMpQ1QnAhOTk5NGRjsWEXlZuAZ/M3cSph7fjuF5tAKy5qTFh4qT10To8rY1KUdWDXYnI1Ep7svJ4b9Z6zh3QkemrdpGRU8AVR3WJdFjG1DpOio+SfX5PAM4DWrgTjqmNdmXkcukbf7ByRwavT19L44Q6HN4hkQE2oqkxYedk6OzdPp8tqvo8nkHyjKm2nRk5XDh2Fhv37OeFC49gyMEt2JWRy9XDu1pTU2MiwEnx0QCfxTg8bw5O3jCMCSgnv5DR76WwLT2Hd68azOCuLTijX3vWpWbRtVWjSIdnTK3k5OH+jM/vBcB64HxXojG1hqpyz1dLWLhpL69dOpDBXT0lkiLCwa0bRzg6Y2ovJ62PjgtHIKb2UFVemraGL+Zt5pYTejCyT9tIh2SM8XJSfNQSeAAYjqcV0kzgYVXd7XJsgWKyfgoxaue+HO74fBEzVu3ijH7tueWEHpEOyRjjQ1TLtTYtvYPID3gGv/vAu+oS4FhVPdHl2IJKTk7WuXPnRjoME8T8jWl8t2Q7K7ZnMH9DGvlFRdx7Wm8uGdLZKpONiQARSVHVZH/bnNQptFPVR3yW/y0iF4QmNFPTbdqzn4v/9weFqvRo05iRfdpy3THd6N7G6g2MiUZOksIUEbkQ+NS7/FdgsnshmZpCVbn36yXECUy9/VjaN2sQ6ZCMMUE4GVlsNPAhkOf9fAxcJyIZIrLPzeBMbBu/cCvTV+3in6f0tIRgTIxw0vrIBp0xlZaWlcfDE5ZxRKdmXDY0KdLhGGMcctQJTUSaAz3wDHMBgKrazGumnOy8QiYu3sZbM9eRnp3PuHMPJz7OKpONiRVOmqReA9wCdAQWAEcCs4DjXY3MxJTcgkJe+HE17/++gYycAg5u1Yhnzu9Hr7ZNIx2aMaYSnLwp3AIMAn5X1eNEpBfwqLthmViyZmcGYz5awLJt+zi9bzsuO7ILg7u2sOamxsQgJ0khR1VzRAQRqa+qK0Skp+uRmainqnw4eyOPfLuMhvXq8OblyZxw6EGRDssYUw1OksJmEWkGfA38ICJpwAY3gzLRr6hIefjbZbzz23pG9GjFM+f1o03ThOAHGmOimpPWR2d7f31QRH4CEoHvXY0qCBvmIrLyC4u46/NFfDl/C1cP78o9fzmUOKtMNqZGCNpPQUSeEZHeAKo6XVXHq2qe+6FVTFUnqOq1iYmJkQyjVsrJL+SGD1L4cv4W7jj5EO49zRKCMTWJk+Kj5cD/RKQO8DbwkaqmuxuWiUZ79+dx7fspzFm/h0fO6sNlR9p0mcbUNE5mXntDVYcBfwOSgEUi8qGI2JDatcjaXZmc/epvLNi4lxcu7G8JwZgayskwF4hIPNDL+0kFFgK3icjHLsZmosRvf6Zy9qu/kZ6dz4ejh3BGv/aRDskY4xInndeeA04HpgGPqups76YnRGSlm8GZyFJVPpmziXu/XkLXVo148/JBdG7ZMNJhGWNc5KROYRFwr6pmFa8QkY2q2hkY7FpkJqLmrt/DU5NX8se6PYzo0YpXLhlA04S6kQ7LGOMyJ01S3/azWrzbrMK5hlm6NZ1npqxi2oqdtGpcn4fOOIxLhnSmTryjkkZjTIxzNCCeH4GnazMxJye/kHu+WsIX8zaT2KAud43sxeVHdaFhvar+iRhjYlGF/8eLyG0VbQJs2qwaJCe/kNHvzeWX1anccGw3rj+mG4kNrKjImNoo0NfAQPMovBDqQExk7M8r4Jp35zJr7W6ePLcv5w/qFOmQjDERVGFSUNWHwhmICb/M3AKuensOczfs4dnz+3F2/46RDskYE2EV1h6KyL3eyXUq2n68iJzuTljGbfty8vnbm3+QsjGNFy7sbwnBGAMELj5aDHwrIjnAPGAXnpnXegBHAD9i8yrEpLSsPC5/ezbLt+3jlYsHMLJP20iHZIyJEoGKj74BvhGRHsAwoB2wD/gAuFZVs8MTogmllA1pjPloPrsycvnvJQM5sbfNf2CMOcBJP4XVwOowxIKInAWcBjQF3lTVKeG4bm1QVKSM/WUtT01eSbvEBD69fihHdGoW6bCMMVHG9R5JIvKWiOwUkSVl1o8UkZUiskZE7gZQ1a9VdTRwPXCB27HVFrszc7nq3Tk8/t0KTjr0ICaOGWEJwRjjVzh6Jr0DvAy8V7zCO8DeK8BJwGZgjoiMV9Vl3l3u9W431fTH2t2M+Xg+aVn5PHLmYVx6ZBebO9kYU6FArY8uEpGW1b2Aqs4A9pRZPRhYo6prvRP2fAycKR5PAN+p6rzqXrs2KyxSXpq6mov+9zsN6sbz5Y1HcdnQJEsIxpiAAr0pdAY+E5G6wFTgO2C2qoZiiIsOwCaf5c3AEOBm4EQgUUS6q+prZQ8UkWuBawE6d+4cglBqnp0ZOfzjkwX8umY3Z/Rrz6PnHE7j+jZchTEmuECtj57AMzx2EzwP6quA10RkOZ45mier6o5QBqOqLwIvBtlnLDAWIDk52cZgKuPXNanc8vECMnLyefycw7lgUCd7OzDGOOak9VEG8JX3g3e+5lPx1BGcUsXrbgF8x1Po6F1nqqigsIgXp67mpZ/W0K11Y8ZdM4SebQONVGKMMeVVukzBWxm8DHimGtedA/QQka54ksGFwMVODxaRUcCo7t27VyOEmmN7eg5jPp7P7HV7OG9gRx468zAb3dQYUyXhaJL6ETAL6Ckim0XkalUtAG4CJgPLgU9VdanTc6rqBFW9NjEx0Z2gY8hPK3fylxd/YcmWdJ49vx9PndfPEoIxpspcf3qo6kUVrJ8ETHL7+jVVQWERT09ZxWvT/6RX2ya8fPEAurexEc2NMdVTpaQgIo1VNTPUwVTi+rW6+Cg1M5ebP5zPrLW7uWhwZx4Y1ZuEuvGRDssYUwNUtfhoWfBd3FObi48WbNrLqJdmMm9jGk+f14/HzjncEoIxJmRs5rUY8tHsjTzwzVLaNK3PFzccRZ8OtS8pGmPcFaj46FHgKaDAzzabxT2McvILeeCbpXwydxMjerTixQv707xRvUiHZYypgQIlhXnA16qaUnaDiFzjXkjB1aY6hS17s7nhgxQWbU7npuO684+TDiE+zjqjGWPcEegb/5XAhgq2JbsQi2O1pU7h1zWpjHppJmt3ZTH2soHccUpPSwjGGFcFGuZiZYBtIR3ewpSmqrw2fS1PTV5Bt9aNee2ygXRrbdU4xhj3WS+nKJOZW8A/P1vId0u2c9rh7Xjyr31pZIPZGWPCJCafNjW1TmHNzkyue38u63fv556/HMo1I7raYHbGmLCKyVZENbFO4fsl2zjrlV/Zuz+f968ezOijD7aEYIwJO0dJQUQu9f1pQqewSHni+xVc/8E8urVpzISbh3NUt1aRDssYU0s5fVO4rcxPEwJ7svK4/K3Z/PfnP7l4SGc+ve5I2jdrEOmwjDG1WGXrFKw8I0QWbd7LDR/MY1dmLk+e25fzB3UKfpAxxrjMKpoj4NM5m7j3myW0blyfz68fSt+OzSIdkjHGAFbRHFa5BYX835eLufOLRQxKas6Em4dbQjDGRJWYfFOIRdvSs7n+g3ks3LSX64/pxh0nH0Kd+JjMycaYGsxpUljl/VlhL2dTsd/+TOXmD+eTk1/Ia5cOYGSfdpEOyRhj/HKUFFT1Qt+fxhlV5Y1f1vH49ytIatmQ1y87ku5tmkQ6LGOMqZAVH7kkK7eAO79YxMRF2xh5WFuePr8fjW24CmNMlLOnlAvW7Mzghg/m8eeuTO4a2Yvrj7HeycaY2BCTSSGam6SOX7iVu79YRIO68bx/9RCGdbfeycaY2OF0mIvhInKl9/fWItLV3bACi8YmqbkFhdz/zRLGfDSf3u2aMnHMCEsIxpiYE/RNQUQewDOpTk/gbaAu8AEwzN3QYsfmtP38fdw8Fm5OZ/SIrtw5shd1rbmpMSYGOSk+Ohvoj2d6TlR1q4hYExqvn1bs5NZPFlBUpNbc1BgT85wkhTxVVRFRABFp5HJMMaGwSHnuh1W8/NMaDm3XlP9eMoCkVnZrjDGxzUlS+FREXgeaicho4Crgf+6GFd12ZeRyy8fz+e3P3VyQ3ImHzjyMhLrxkQ7LGGOqLWhSUNWnReQkYB+eeoX7VfUH1yOLUrPX7eGmD+eRnp3Pk3/ty/nJNrqpMabmcNqj+QcR+aN4fxFpoap7XI0sgEg0SVVV/vfLWp74fiWdmjfgnSsH07t907Bd3xhjwsFJ66PrgIeAHKAIz5wKChzsbmgVU9UJwITk5OTR4bheenY+//xsIVOW7WDkYW158ry+NE2oG45LG2NMWDl5U7gD6KOqqW4HE42Wbk3nxnHz2JKWzb2nHcrVw7ta72RjTI3lJCn8Cex3O5Boo6p8OncT932zlBYN6/HJdUcysEuLSIdljDGucpIU/g/4zVunkFu8UlXHuBZVhGXnFXLfN0v4PGUzw7u34vkLj6BV4/qRDssYY1znJCm8DkwDFuOpU6jR1u7K5MZx81i5I4MxJ/TglhN6EB9nxUXGmNrBSVKoq6q3uR5JFJi0eBt3fr6IOvHC21cM4tiebSIdkjHGhJWTpPCdiFwLTKB08VHEmqSGWn5hEY9/t4I3Z67jiE7NeOWSAXRo1iDSYRljTNg5SQoXeX/+n8+6iDZJDaVt6dnc9OF8UjakccVRSfzrL4dSr44NZmeMqZ2c9GiO6DDZbvpl9S5u+XgBufmFvHxxf07v2z7SIRljTERVmBRE5HhVnSYi5/jbrqpfuheWu4qKlJemreH5qavo0aYxr14ykO5tGkc6LGOMibhAbwpH42l1NMrPNgViMimkZeVx6ycLmL5qF2f378B/zu5Dw3oxOQGdMcaEXKCnYT0AVb0yTLE4VtWxjxZu2suN4+axKyOX/5zdh4sHd7beycYY4yNQjerIsEVRSZWdjlNVGffHBs57bRYAn98wlEuGdLGEYIwxZQR6U4gXkeZ4BsArJ1aapGbnFXLP14v5ct4WjjmkNc9fcATNG9WLdFjGGBOVAiWFXkAK/pNCTDRJXZeaxQ0fpLByRwb/OPEQbj6+O3HWO9kYYyoUKCksU9X+YYskxCYv3c4dny4k3nonG2OMYzWu2U1BYRFPTVnJ69PX0rdjIq9eMoCOzRtGOixjjIkJgZLCC2GLIkR2ZuQw5qP5/L52D5cM6cz9o3pTv47NnWyMMU5VmBRU9Z0wxlFtc9bv4e/j5rEvJ59nzuvHuQM7RjokY4yJOTFffKSqvDlzHY99t4JOzRvw7lWDObSdzZ1sjDFV4WSO5taquiscwVRWkSo3fTifiYu3cXLvg3j6/H42d7IxxlSDqGrgHURWAeuBT4AvVTUtDHE5IiK7gA1VPDwRSHf52GD7Bdrub5uTdWWXWwFuz69d1Xtp9zF4XKE+rqr30un6YPc2mu9jZY51+z76WxfKv8kuqtra7xZVDfoBBgPPAmuBb4FLnRwXzR9grNvHBtsv0HZ/25ys87M8N1rvpd3H8N7H6txLp+uD3dtovo+VOdbt++jk3rp1Lx1NHKCqs9Uz+9pgYA/wrpPjotyEMBwbbL9A2/1tc7KuOv+uqqrqNe0+huaalTmuqvfS6Xqn99tN0fz/dmXWR+Rv0knxUVPgbOBCoBvwFfCpqqa4H56pLhGZq6rJkY4j1tl9DA27j6Hj1r100vpoIfA18LCqzgp1AMZ1YyMdQA1h9zE07D6Gjiv30smbgqiqikhjAFXNdCMQY4wxkeekTuEwEZkPLAWWiUiKiPRxOS5jjDER4CQpjAVuU9UuqtoZuB17BTTGmBrJSVJopKo/FS+o6s9AI9ciMsYYEzFOksJaEblPRJK8n3vx9FcwMUhEzhKR/4nIJyJycqTjiVUicqiIvCYin4vIDZGOJ5aJSCMRmSsip0c6llglIseKyC/ev8ljq3MuJ0nhKqA18KX309q7zkQJEXlLRHaKyJIy60eKyEoRWSMidwOo6teqOhq4HrggEvFGq0rex+Wqej1wPjAsEvFGq8rcR6+7gE/DG2X0q+R9VCATSAA2V+u6wVofmegnIkfj+YN4T1X7eNfFA6uAk/D8kcwBLlLVZd7tzwDjVHVeZKKOPpW9jyJyBnAD8L6qfhihsKNOZe4j0AFoiedhlqqq30Yk6ChUyfu4QlWLROQg4FlVvaSq13UyIN4hwB1Aku/+qnp8VS9qQktVZ4hIUpnVg4E1qroWQEQ+Bs4UkeXA48B3lhBKq8x9xDMz4XhgvIhMBCwpeFXyPjbGU0fZG8gWkUmqWhTOeKNVZe5j8Zc9IA2oX53rOum89hnwGvAGUFidi5mw6gBs8lneDAwBbgZOBBJFpLuqvhaJ4GKI3/voLbc9B8//gJPCH1bM8XsfVfUmABG5As+bgiWEwCr6ezwHOAVoBrxcnQs4SQoFqvrf6lzERA9VfRF4MdJxxDpvK7yfIxxGjaExNqlXtFHV4jrfanNS0TxBRG4UkXYi0qL4E4qLG1dtATr5LHf0rjOVY/cxNOw+hobr99HJm8Ll3p//9FmnwMGhDMSE3Bygh4h0xfNHcyFwcWRDikl2H0PD7mNouH4fg74pqGpXPx9LCFFERD4CZgE9RWSziFytqgXATcBkYDmekW2XRjLOaGf3MTTsPoZGpO6jkwHx/uZvvaq+F8pAjDHGRJ6T4qNBPr8nACcA8wBLCsYYU8NUuvOaiDQDPlbVka5EZIwxJmIcTcdZRhbQNdSBGGOMiTwnPZon4GltBBAPHIqNU2KMMTWSk4rmY3wWC4ANqlqtAZeMMcZEJydNUqcDK4FEoAWexGCMMaYGCpoUROQaYDaecV7+CvwuIjZ0tqkRRKSliCzwfraLyBbv75ki8qpL17y1oqbeAY55WkRsEErjOifFRyuBo1R1t3e5JfCbqvYMQ3zGhI2IPAhkqurTLl6jDp4m3QO8HZGcHBOPZziD/6mqTYxkXOWk9dFuIMNnOcO7zpgayzuT1bfe3x8UkXe9M1ttEJFzRORJEVksIt+LSF3vfgNFZLqIpIjIZBFp5+fUxwPzVLVARLqJyDyfa/YoXhaR9SLyhHf5PFXdALQUkbbu/+tNbeYkKawB/vD+j/EA8DuwSkRuE5Hb3A3PmKjRDc8D/QzgA+AnVT0cyAZO8yaGl4C/qupA4C3gP37OMwxIAVDVP4F0ETnCu+1K4G2ffXer6gBV/di7PA+b5c24zEmP5j+9n2LfeH82CX04xkSt71Q1X0QW42ma/b13/WI8E1D1BPoAP4gI3n22+TlPOzxj1hR7A7jS+wXrAjyTqBT7pMyxO4H21ftnGBNY0KSgqg+FIxBjolwugHfKw3w9UBlXhOf/IwGWqurQIOfJxjNcTLEvgAeAaUBKcd2dV1aZYxO8xxvjGied11oDdwKH4fPHbNNxGlPKSqC1iAxV1Vne4qRD/IxguRzoXrygqjkiMhn4L3B1kGscgmcmRGNc46ROYRywAs/QFg8B6/GM6W2M8VLVPDxNtp8QkYXAAuAoP7t+BxxdZt04PG8cUyo6vzfJdAfmhiJeYyripElqiqoOFJFFqtrXu26Oqg4KeKAxxi8R+Qq4U1VXe5fvABJV9b4Ax5yNpxlrhfsYEwpOKprzvT+3ichpwFY8PZuNMVVzN54K59XeBFHcsimQOsAzbgdmjJM3hdOBX/DMC/oS0BR4SFXHux+eMcaYcHKSFFqWaRFhjDGmhnJS0fy7iHwmIn8RbwNsY4wxNZOTpHAIMBa4DE8Z6KMicoi7YRljjImESk3HKSLH4eni3xhPk7u7VXWWO6EZY4wJN0d1CsCleN4UdgBvAuOBI4DPVNWm5jTGmBrCSfHRLDwtjs5S1dNU9UvgJlWdC7zmanTGGGPCysmbgmiZnURko6p2djUyY4wxYedkOk5/WcNaIRljTA3kpPjIH+e108YYY2JGhcNciEgG/h/+AjRwLSJjjDERU6kmqcYYY2q2qhYfGWOMqYEsKRhjjClhScEYY0wJSwrGGGNKWFIwxhhTwpKCMcaYEv8PALjV1pwbEikAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Make some plots\n",
"plt.plot(timeyr, 1 / np.array(lyapunov_CN)/twopi)\n",
"plt.xlabel('Time (yr)')\n",
"plt.ylabel('Lyapunov Time = 1 / (Lyapunov Characteristic Number) [year]')\n",
"plt.yscale('log')\n",
"plt.xscale('log')\n",
"plt.show()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fewer-charger",
"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.9.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment