Skip to content

Instantly share code, notes, and snippets.

@hannorein
Created April 8, 2021 00:50
Show Gist options
  • Save hannorein/6d38796e09b3d2485e354757161f80f9 to your computer and use it in GitHub Desktop.
Save hannorein/6d38796e09b3d2485e354757161f80f9 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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.16.0\n",
"Apr 7 2021 19:03:43\n"
]
}
],
"source": [
"'''Let's try calculating the Lyapunov timescale using Rebound!'''\n",
"import matplotlib.pyplot as plt\n",
"import rebound\n",
"print(rebound.__version__)\n",
"print(rebound.__build__)\n",
"import numpy as np\n",
"from numpy import radians as rd"
]
},
{
"cell_type": "code",
"execution_count": 2,
"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'... "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rein/git/rebound/rebound/horizons.py:145: RuntimeWarning: Warning: Mass cannot be retrieved from NASA HORIZONS. Set to 0.\n",
" warnings.warn(\"Warning: Mass cannot be retrieved from NASA HORIZONS. Set to 0.\", RuntimeWarning)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"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": null,
"id": "waiting-little",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<ipython-input-49-fcc8ef9988fa>:31: RuntimeWarning: divide by zero encountered in double_scalars\n",
" lyapunov_CN[i] = 1./time * np.log(d/1.) # divide by initial perturbation (here 1.)\n",
"<ipython-input-49-fcc8ef9988fa>:31: RuntimeWarning: invalid value encountered in double_scalars\n",
" lyapunov_CN[i] = 1./time * np.log(d/1.) # divide by initial perturbation (here 1.)\n"
]
}
],
"source": [
"\n",
"twopi = rd(360) # 6.283185307179586\n",
"dt_out = 100 * twopi # Years * 2 Pi, output interval\n",
"tmax = 1e7 * 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",
"sim.ri_whfast.safe_mode = 0\n",
"orbits = sim.calculate_orbits() # Calculate the orbits\n",
"sim.dt = 0.08 * 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.log(d/1.) # divide by initial perturbation (here 1.)\n"
]
},
{
"cell_type": "code",
"execution_count": 52,
"id": "reverse-monster",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzKElEQVR4nO3dd3hUZfbA8e9JpyVI7wQIIEVpEUFAAUUpAuraQCzoYsdecHXV/dl31VWsi65iRdC10K0giCCEJiBI6E2kdwKEnN8fMxMmYTK5SaYlcz7Pcx9z37lz51xDcnLv+77nFVXFGGOMKUhMuAMwxhgT2SxRGGOM8csShTHGGL8sURhjjPHLEoUxxhi/LFEYY4zxyxKFMcYYvyxRGGOM8Ssu3AE4ISIXAf2AZOC/qvpNeCMyxpjoEfQ7ChF5R0S2icjSfO29ReR3EVklIiP8nUNVv1TVYcDNwBXBjNcYY0xeEuwSHiJyNnAAeF9VW7vbYoGVQC9gEzAPGATEAs/kO8X1qrrN/b4XgI9UdUFQgzbGGJMr6I+eVHWGiKTma+4IrFLVNQAi8gkwUFWfAS7Mfw4REeBZYIqTJFGtWjVNTc3/kcYYY/yZP3/+DlWtnr89XH0UdYGNXvubgDP9HD8cOA9IEZE0VX0z/wEiciNwI0CDBg3IyMgIYLjGGFP2ich6X+2lojNbVUcCIws5ZhQwCiA9Pd1K4hpjTICEa3jsZqC+1349d1uJiEh/ERm1d+/ekp7KGGOMW7gSxTygqYg0EpEE4EpgfElPqqoTVPXGlJSUEgdojDHGJRTDY8cAs4HmIrJJRG5Q1WzgduBrYDkwTlWXBTsWY4wxRReKUU+DCmifDEwO5GeJSH+gf1paWiBPa4wxUa1MlfCwR0/GGBN4hd5RiEgM0AaoAxwGlnomwJU1q7cfYOeBo3RsVCXcoRhjTMQoMFGISBPgQVzzFzKB7UAS0ExEDgH/Ad5T1ZxQBBoK/5y6gu+Wb+P+C5pz09mNcc3zM8aY6Obv0dOTwAdAE1W9QFWHqOqlqno6MABIAa4ORZBOlXR47L8ua8P5LWvy7JQVDHt/PnsPHQtwhMYYU/r4rfXkfuzUSVV/Dl1IJZeenq7FnZmtqrw7ax1PT15O7cpJvD64A6fVsz4PY0zZJyLzVTU9f7vfzmz3Y6XXghZVBBIRru/aiHE3d+b4ceUvb/zMB3PWE+ziicYYE6mcjHr6XkT+IlH2wL59g1OYeEc3OjWpyt+/XMpdYxdx8Eh2uMMyxpiQc5IobgI+BY6IyD4R2S8i+4IcV7EEuoRHlQoJjL7uDO7t1YwJi7cw8LVZZP65PyDnNsaY0qLQRKGqlVQ1RlUTVDXZvZ8ciuCKKhjzKGJihOHnNuXDG85kz6GjDHh1Fl8s3BSw8xtjTKQrNFGIyP9EpK+7YztqnZVWjUl3dOO0uincPXYxD32+hKxjx8MdljHGBJ2TX/5vAFcBmSLyrIg0D3JMEatmchIfDzuTm89pwpi5G/jLGz+zYeehcIdljDFB5eTR03eqehXQHlgHfCciP4vIUBGJD3aARRGKMuNxsTGM6HMqb12TzsZdh+j3yky+XrY1aJ9njDHh5uhxkohUBa4D/gosBF7GlTi+DVpkxRDKWk+9WtZk0h3dSK1agZs+mM9Tk37j2PEyM0ndGGNyOemj+AKYCZQH+qvqAFUdq6rDgYrBDjCS1a9Sns9u6czVnRry1sy1DBo1h617s8IdljHGBFSBiUJE6ri/HKmqLVX1GVX9w/sYXzP4ok1iXCxPXNSal69sy29/7KPfyJn8lLkj3GEZY0zA+LujeFtE5gAXiEh3ESkV62uHy8C2dRl/exeqVEjg6nd+4aXvVnI8x2ZzG2NKvwIThar2BboD04GLgTki8rmI3CgiDUITXumSVqMSX93ehYvb1uWl7zK5fvQ89hw6Gu6wjDGmRPwWBTzpYJFGQB+gN1BLVTsGK7Di8FrhblhmZmbY4lBVPp67gcfHL6NmchJvDulA67pWWNAYE9kKKgpYpESR74QJqhqRfy6XpHpsIC3csJtbP1rAzoNHeWJgK644w27EjDGRq8jVY0XkJ/d/93vXePL8N1KTRCRp1+AUJg7vSsfUKjz4vyU8+NmvNpvbGFPq+Ouj6Or+byXvGk+RXOspElWtmMh713fkth5NGJuxkcvenM3GXTab2xhTejidcNdeRO4QkeEi0i7YQZU1sTHC/Re4ZnOv23mQ/q/+xPTfy+Sy48aYMsjJhLtHgfeAqkA1YLSIPBLswMqiXi1rMuH2rtRKTmLo6Hm8/F0mOTaE1hgT4QrtzBaR34E2qprl3i8HLFLViC0OGCmd2QU5fPQ4f/tiCV8s3EyP5tX59xVtqVw+IdxhGWOiXLGWQnXbAiR57ScCmwMVWCCFoihgIJRLiOXFy9vwxMBW/LRqB/1f/YmlmyM7ZmNM9HKSKPYCy0RktIi8CywF9ojISBEZGdzwiiaURQFLSkS4unMqY2/qTLZ7be5xGRvDHZYxxpzESVmOL9ybx/TghBKd2ruH0A4fs5AHPvuVhRv28Fj/liTFx4Y7NGOMARwkClV9LxSBRLOqFRN5//qOvPDtSt6YvpplW/by+lXtqXdK+XCHZowxfifcJYvIMyLygYgMzvfa68EPLbrExcbwYO9TGXV1B9ZuP0j/V35ixsrt4Q7LGGP89lG8CwjwP+BK99rZie7XOgU9sih1fqtajB/elRqVkrj23bm88r0NoTXGhJe/RNFEVUeo6peqOgBYAPzgXu3OBFGjahX44razGNimDi98u5Lr35vHroNWMcUYEx7+EkWiiOS+rqpPAW8BM3BNvjNBVD4hjn9f0ZYnBrbi51U76fvyTOat2xXusIwxUchfopgA9PRuUNXRwL2A/XkbAp4htJ/fehZJ8TFcOWoOr09fZY+ijDEhVewy45Es0mdmF8f+rGM89PkSJv76B+c0q86Ll7ehasXEwt9ojDEOFafM+BDvR08+Xm8iIl0DFaDxr1JSPK8MaseTF7Vm9pqd9B05k1/W7Ax3WMaYKOBvHkVVYKGIzAfmA9txlfJIA84BdgAjgh5hEXitcBfuUIJCRBjSqSHtGlTm9o8XMuitOdzTqxm3dk8jJkbCHZ4xpozy++hJRGJx9VN0AWoDh4HlwBRV3RCSCIuhLD56yu/AkWz+9vkSxi/eQrem1fj3FW2pZo+ijDElEPClUCNZNCQKcK3N/cm8jTw+fhkp5eJ5+cp2dG5iA9KMMcVTkuqxJkKJCIM6NuDL27pQMTGOq96ew4vf/E728Zxwh2aMKUMsUZQBLWonM2F4Vy5uV4+RP6zi8v/YcqvGmMDxmyhEJEZELg9VMKb4KiTG8cLlbRg5qB2Zfx6gz8sz+WpRRC4bYowpZfwmClXNAR4IUSwmAAa0qcPkO7txaq1K3PnJIu4eu4j9WcfCHZYxphRz8ujpOxG5T0Tqi0gVzxb0yEyx1a9Snk9u7MRd5zXlq0Wb6TtyJgs27A53WMaYUsrJmtlrfTSrqjYOTkglFy2jnpzIWLeLOz9ZxNZ9Wdx1blNu7ZFGrM25MMb4UOxRT6rayMcWsUnC5JWeWoUpd3Wj32m1eeHblQx6aw6b9xwOd1jGmFKk0EQhIuVF5BERGeXebyoiFwY/NBMoyUnxvHxlW168vA3LNu+lz0szmPTrH+EOyxhTSjjpo3gXV7XYs9z7m4EngxaRCQoR4ZL29Zh8ZzcaV6/IbR8v4IHPFnPwSHa4QzPGRDgniaKJqv4TOAagqodwrXxnSqGGVSvw6c2dub1HGp/O38SFr/zEkk17wx2WMSaCOUkUR0WkHKDgqhoLHAlqVCao4mNjuO+C5owZ1omsY8e55I1ZvPnjalvnwhjjk5NE8RgwFagvIh8B3xPCuRUi0kJE3hSRz0TkllB9bjTo1LgqU+7sRq+WNXl2ygqufucXtu7NCndYxpgI42TU07fAJcB1wBggXVWnOzm5iLwjIttEZGm+9t4i8ruIrBIRv6XKVXW5qt4MXI6riq0JoMrlE3htcHue+8tpLFi/hz4vz+CbZVvDHZYxJoI4rfV0DnAu0APoVoTzjwZ6eze4S5e/BvQBWgKDRKSliJwmIhPzbTXc7xkATAImF+GzjUMiwhVnNGDiHV2pe0o5bvxgPk9O/I3j9ijKGIOzCXev41qsaIy76Qpgtare5ugDRFKBiara2r3fGXhcVS9w7z8EoKrPODjXJFXtV8BrNwI3AjRo0KDD+vXrnYRn8jmancPTk5cz+ud19G5Vi5eubEtSfGy4wzLGhEBBE+78rXDn0RNooe6MIiLvActKEEtdYKPX/ibgzIIOFpHuuB59JeLnjkJVRwGjwDUzuwTxRbWEuBgeH9CK+lXK8+Sk37j6v7/wznVnUCkpPtyhGWPCxEmiWAU0ADx/otd3t4WEuz9keqg+z7jc0LURNZMTueuTRQx5+xfeu74jlcsnhDssY0wYFNhHISITRGQ8UAlYLiLTRWQarqVQK5XgMzfjSjYe9dxtJSYi/UVk1N69Ni8gEC48vQ5vDunA8j/2c+WoOew4YKOijYlGBfZRiMg5/t6oqj86+oCT+yjigJW4Osc3A/OAwapaksdZeVhRwMCambmdYe9nULdyOT4e1omayUnhDskYEwRFLgqoqj96b8BCYInX5uRDxwCzgeYisklEblDVbOB24GtcdyfjApkkTOB1a1qd94Z2ZOveLAa9NYft++3Owpho4mTU043A/wFZQA6u8h0RWWZcRPoD/dPS0oZlZmaGO5wyZ+7aXVzzzi+kVq3AJzd2sj4LY8qYYpcZB+4HWqtqqqo2juQy46o6QVVvTElJCXcoZVLHRlV4+5ozWLPjINe8M5d9tnKeMVHBSaJYDRwKdiCmdOjatBpvXNWe37bsY+i786z6rDFRwEmieAj4WUT+IyIjPVuwAysOG/UUGue2qMkrg9qxcMNuhr2fQdax4+EOyRgTRE4SxX+AH4A5wHyvLeLYo6fQ6XNabV64vA2z1+zk5g/ncyTbkoUxZZWTCXfxqnpP0CMxpc7F7eqRdSyHhz5fwh1jFvLq4PbExzotH2aMKS2c/FRPEZEbRaS2iFTxbEGPzJQKgzo24LH+Lfl62Z/cO26xFRI0pgxyckcxyP3fh7zaFIi4kU9ew2PDHUpUGdqlEVnHcnhu6goS42J47i+nExNjiyAaU1YUmihUtVEoAgkEVZ0ATEhPTx8W7liizS3dm3D42HFGfp+JAs9echpx9hjKmDKh0EQhItf4alfV9wMfjinN7j6vKTECL32XyZ/7snjh8jbUqGTlPowp7Zz8yXeG19YNeBwYEMSYTCklItx1XjOeveQ05q7dRZ+XZjJtxbZwh2WMKSEnj56Ge++LSGXgk2AFZEq/Kzs2oH3DU7hjzEKGjp7H9V0a8WCf5iTG2QJIxpRGxXmIfBCIyH4Lm3AXOZrVrMSXt3Xh2s4NeWfWWi567WcWbdwT7rCMMcXgpCjgBFyjnMCVWFriqvg6IsixFZuVGY8s3y//k799sYRt+48wuGMD7u7VjGoVE8MdljEmn5Ishfq819fZwHpV3RSwyEyZd26LmnzXqAovfLOSD+as58uFm7npnCb8tVsjyic4+SdojAmnQu8oSiO7o4hcq7cf4J9TV/D1sj+pXimRu89rxhVn1CfW5l0YE3YF3VH4W+FuLSceOeWnqtokgPEFlCWKyDd//S6embyCjPW76dWyJq8Obmed3caEWXHWo0gn79DYM4EXcC1ctCgIMZoo0qFhFT69uTOPXtiSb3/7kxH/W0JZvLs1piwo8AGxqu4EEJEY4GpcCxgtAvqp6m8hia6IrIRH6SIiXN+1EYeOZvP8NytpWLU8d53XLNxhGWPyKfCOQkTiReQm4DdcE+0uUtUhkZokwMqMl1a39UjjknZ1efn7TH7K3BHucIwx+fgbcrIW1yinl4ANwOkicrrnRVX9PLihmWghIjx5cWuWbN7LXWMX8tnNZ5FarUK4wzLGuPnro/gOmAa0Afrn2y4MfmgmmpRPiOONIR04nqNc9fYvbNlzONwhGWPcbHisiShLN+9l0Kg5VKuUyNibOllRQWNCqDijnowJudZ1Uxh9/Rls3ZvF1W/P5c99WeEOyZioZ4nCRJwODavw9rXprN91kPP/PYPxi7fY0FljwqhMJQorClh2dEmrxuQ7utG4egXuGLOQWz9aQMa6XZYwjAmDQhOFiNzmLi3u2T9FRG4NalTFZMNjy5bG1Svy6U2deaB3c35YsY1L35zN0NHz2LT7ULhDMyaqOLmjGKaqezw7qrobsKVGTUjExcZwa/c0Mh45j0f6tWDu2l2c/+8ZtiCSMSHkJFHEikhuxTYRiQUSgheSMSerlBTPX7s15pu7z6Zx9Qrc/OF85q/fHe6wjIkKThLFVGCsiJwrIucCY9xtxoRcvVPK8/71Z1IzOYlbP5rPtv02KsqYYHOSKB7ENfHuFvf2PfBAMIMyxp8qFRJ4c0gH9h4+xu0fLeRodk64QzKmTCs0Uahqjqq+oaqXurf/qOrxUARnTEFa1knmub+cztx1u3jkS6s8a0wwFVjrSUTGqerlIrIEH+tSqOrpPt5mTMgMbFuX1dsPMvL7TGolJ3HXec2IsQWQjAk4f0UB73T/1+o6mYh193lN2bz7MCN/WMWs1Tu59/xmtKqTQkq5+HCHZkyZUeCjJ1X9w/3lraq63nsDInIehYk+IsLzl53Oi5e3Yc32Awx+6xfa/OMbrv7vLxzJtiekxgSCk87sXj7a+gQ6kECwmdnRSUS4pH09Zo3oyairO3Br9ybMzNzBi9+uDHdoxpQJ/hYuusXdP3GqiPzqta0Ffg1diM7ZzOzoVj4hjvNb1eKB3qdy5Rn1eWvGGuav3xXusIwp9fzdUXyMa+2Jr8i7FkUHVR0SgtiMKbZHLmxJ7ZRy3DtuMYeOZoc7HGNKNX99FHtVdR3wCLDV3TfRCBjiXfvJmEhUMTGOf112Out2HuK+TxeTk2PDZ40pLid9FP8DjotIGjAKqI/rbsOYiHZWk2o80q8Fk5ds5c6xixi/eAtLN++1pGFMEfkbHuuRo6rZInIJ8IqqviIiC4MdmDGB8Ndujdl18Chv/riaCYu3ANC/TR1evqKtzbkwxiEnieKYiAwCrsHVRwFgg9RNqfFA71MZ3rMp63cdZPyiLbw+fTVt6qXw126Nwx2aMaWCk0dPQ4HOwFOqulZEGgEfBDcsYwKrXEIsp9ZK5v4LmnN2s+q8MX01WcdsnoUxTjip9fSbqt6hqmPc+2tV9bngh2ZM4IkIt5zThJ0Hj3Lvp4vZn3Us3CEZE/Gs1pOJOp0aV+H+C5rz4rcr+WH5NlKrVeD6Lqlcll4/3KEZE5Gs1pOJOiLCbT3SOKtJVT6YvZ4Zmdv5+1dL6ZJWjTqVy4U7PGMijt9aT+7V7Ebnr/XknlNhTKnWrsEpvHhFW764tQsA93262MqVG+OD3z4K97oTOSJiNTFMmVW/Snke7teSn1fv5G9fLLF+C2PycTI89gCwRES+BQ56GlX1jqBFZUyIDe7YgBV/7OOjXzYwZu5GOjWuwptDOlC5vC0Pb4yTRPG5ewsbEakA/Ag8rqoTwxmLKZtiY4QnBrbmzMZVWbB+N+/PXsdL32Xy+IBW4Q7NmLArNFGo6nsiUg5ooKq/F+XkIvIOrs7wbara2qu9N/AyEAu8rarPFnKqB4FxRflsY4oqJkYY0KYOA9rU4Uh2Dh/9sp7rzkoltVqFcIdmTFgVOo9CRPoDi4Cp7v22IjLe4flHA73znS8WeA3XmhYtgUEi0lJEThORifm2GiLSC/gN2Ob0oowpqTvPbUpCbAx/eeNnuv3zB/719Qrr6DZRy8nM7MeBjsAeAFVdBDiqfaCqM4D8CwJ0BFap6hpVPQp8AgxU1SWqemG+bRvQHegEDAaGiYiTmI0pkVopSYy9qTOVkuLYuOswr01bzfSV28MdljFh4eSX7jFVzb9kXE4JPrMusNFrf5O7zSdVfVhV78JVsfYtVfX52SJyo4hkiEjG9u32A21KrnXdFKbd152VT/ahbuVyPPS/Jew6eJTs4yX5529M6eMkUSwTkcFArIg0FZFXgJ+DHNdJVHW0v45sVR2lqumqml69evVQhmbKMBEhIS6GR/u3ZOu+LNo/8S1pD09h1bYD4Q7NmJBxkiiGA62AI7j+qt/LiVnbxbEZ15oWHvXcbSVma2abYDmvRU3+2rURSfGuH5n3Z68Lb0DGhJAU1kEnIpep6qeFtfl5fyow0TPqSUTigJXAubgSxDxgsKouK3r4vqWnp2tGRkagTmdMHjeMnsf3K7bRqFoFHunXgnNb1Ax3SMYEhIjMV9X0/O1O7igectjm60PHALOB5iKySURuUNVs4Hbga2A5MC6QScKYYHu0f0uS4mNYu+MgN7yXwa6DR8MdkjFBVeAdhYj0AfoClwNjvV5KBlqqasfgh1c07qG8/dPS0oZlZmaGOxxThh3PUT7N2MiIz5dwdaeGPHJhCxJiYxCxVfNM6VWcO4otQAaQBcz32sYDFwQjyJJS1QmqemNKipWmMsEVGyNc2bEB152Vygdz1tP8kak8+pXdGJuyqcCZ2aq6GFgsIl8AB90FAj0T5hJDFJ8xEe3+C5pz9HgOk379gw/mrOfas1JpXK2CrcdtyhQnfRTfAN5F+ssB3wUnHGNKlwqJcTx98WlMHN4VgPNe/JFOz3zP0Wyba2HKDieJIklVcweNu78uH7yQis+Gx5pwqV+lPPed3wyAbfuP8PEvtmSLKTucJIqDItLesyMiHYDDwQup+KyPwoTT7T2bsvaZvtRMTuTxCb+ROmIS+2xtC1MGOEkUdwGfishMEfkJ1wio24MalTGllIjwUJ8WuftfL90axmiMCYxCJ9wBiEg80Ny9+7uqRvSfSTbhzoSbqtLl2R/YsjeLjo2q8Pa16SQnxYc7LGP8KsmEO3AliZZAe1xlwa8JZHCBYn0UJlKICA/2ORWAuWt38d+Za9mw81CYozKmeJyU8HgMV6nvlsBkXOtI/KSqlwY9umKyOwoTKfZlHePWDxfw06odAPx4f3caVrWFkExkKskdxaW46jJtVdWhQBvAeouNcSA5KZ7/G9iKNvUrA3DtO3Ot8qwpdZwkisPuNSCyRSQZ10pz9Qt5jzHGrXH1inx1Wxf6nV6bdTsPMfitObZanilVnCSKDBGpDLyFq4THAlyF/iKO9VGYSPZ4/1Z0bFSFbfuPMOS/v7Bpt/VZmNLBbx+FuCqc1VPVje79VCBZVX8NTXjFY30UJlIdOJJN68e+BuDSDvV4/rI2YY7ImBOK1Uehriwy2Wt/XaQnCWMiWcXEON66Jp1KiXF8Nn8TF74yk6xjx8MdljF+OXn0tEBEzgh6JMZEiV4tazJyUDsAlm7exwvf/G7JwkQ0J4niTGC2iKwWkV9FZImI2F2FMSXQ49QaLPx7L0TgrZlreW3aqnCHZEyBnCSKC4AmQE+gP3Ch+7/GmBI4pUICH91wJgCv/LCKaSu2hTkiY3wrNFGo6npVXY+rEKB6bRHHRj2Z0uastGrceW5TAIaOnkf2cStPbiJPoYlCRAaISCawFvgRWAdMCXJcxWLVY01pNPjMBlRIiAUg7eEp7Dlka3CbyOLk0dMTQCdgpao2wjVLe05QozImitRMTuKnB3vm7v+4cnsYozHmZE4SxTFV3QnEiEiMqk4DThpna4wpvlMqJPDxMFd/xZ2fLGLtjoNhjsiYE5wkij0iUhGYAXwkIi8D9q/YmAA7q0k1ru7UEIBP5m0IczTGnOAkUQzE1ZF9NzAVWI2NejImKJ64qDXtGlRmwfrd4Q7FmFxORj0dVNXjqpqtqu+p6kj3oyhjTBC0q38K89bt5r8/rQ13KMYAzkY9XSIimSKyV0T2ich+EdkXiuCKyobHmrLghm6NAJi69I8wR2KMi5NHT/8EBqhqiqomq2olVU0OdmDFYcNjTVlQt3I5+p5Wix0HTgyTzfxzP8PHLGTou3NZ+ef+MEZnolGcg2P+VNXlQY/EGJMrrUYlpi7dyp/7shg+ZiFz1+7KfW3a79t5c0h7ereuHcYITTQpMFGIyCXuLzNEZCzwJXDE87qqfh7c0IyJXmk1KpKjcO4LP3LgSPZJr9/84QIAXrqiLRe1qxvq8EyU8XdH4T2y6RBwvte+ApYojAmSlrVdT3c9SeLSDvW4rEM91u08yIP/W5J73F1jF1ErJYlOjauGJU4THfwuXFRa2cJFpizo+NR3bNt/hFu6N+HB3qfmtu89fIwznvqOo9kn6kJ9cetZtGtwSjjCNGVIkRcuEpF/ichNPtpvEpFnAx2gMSav+FjXj2fNSol52lPKxbPyyT58d885uW0Xv/4zz0y2rkQTHP5GPfUERvlofwtXqXFjTBDtzzoGuGpB+ZJWoyKzHzpRI+o/M9YwM9PqRJnA85coEtXHcylVzQEkeCEZYwAOu1e9q5bvjsJb7ZRyeZLF1f+dy8Zdh4Iem4ku/hLFYRFpmr/R3XY4eCEZYwDE/fdYufhYv8fVTinHksdPjDXp9s9pHD5qS6uawPGXKB4FpojIdSJymnsbCkxyv2aMCSJ1rw+WFF/4vNhKSfF89Nczc/fvGbcoWGGZKFTgv0BVnQJcBPQARru37sBfVHVy8EMrOivhYcqS7BxXokiM839H4dElrRr39moGwJSlWzlmq+WZAPH7p4qqLlXVa1W1g3u7VlWX+HtPOFkJD1OWeHoIEx3cUXjc1iMt9+tbPpwf6JBMlHL+L9AYExZO7ygAYmKEcTd1BuC75dvY5x45ZUxJWKIwJsIlxhXtx7Rjoyq5X6+zlfJMADgpM149FIEYY3zzTLwritcGtwdcE/EAPl+widQRk1i0cU8gQzNRwsm/wFki8o2I3CAiViPAmBCLKcaspbObVQPgeI4y4NWfuGfcYgAuem0WmVam3BSRkxXumgGPAK2A+SIyUUSGBD0yYwwAIkXPFJWS4unfpg6VEuP4dVPeUYC9/j2DQaPmBCo8EwUc3dOq6lxVvQfoCOwC3gtqVMaYEquTksR+HyXKAWav2cmPK63ch3HGSR9FsohcKyJTgJ+BP3AlDGNMBGtUrULu1y9d0ZYVT/TO8/q178xlwYbdoQ7LlEJO7igWA22B/1PVZqr6oKraAG1jItzp9Srnfn1Bq1okxcey7tl+PHphy9z2S9yd3cb44yRRNFbVu4ElIlIx2AEZYwKjfMKJ+RflvL4e2iU1z3ETFm8JVUimlHKSKFqJyEJgGfCbiMwXkdZBjssYU0Jxsb47wUWENU/3zd0fPmZhqEIypZSTRDEKuEdVG6pqA+BefK9TYYyJIP7mX3jP4AbYe8hmcJuCOUkUFVR1mmdHVacDFQo+3BgTCBOHd83Tn1BUcYVMwPCewf3U5N/IyVHmrdvFf35cbQUFTR5xDo5ZIyJ/Bz5w7w8B1gQvpLxEpDvwBK5HX5+4E5UxZV7ruim0rlv8ApdxDmZ0D+2Syruz1jEuYxPjMjbltj8zZQXvDj2DHs1rFPvzTdnh5I7ieqA68Ll7q+5uK5SIvCMi20Rkab723iLyu4isEpERhZxGgQNAErCpkGONMW7xBfRReHu4b4sCXxv67jyyjtkCSMbZzOzdqnqHqrZ3b3eqqtPB16OBPIO3RSQWeA3oA7QEBolIS/fCSBPzbTWAmaraB3gQ+EdRLs6YaBYXU/jfgXGxMdRO8b0mN8Cpf58KuEqBTF36Bz5WRzZRoNBHTyLSDLgPSPU+XlV7FvQer2NmiEhqvuaOwCpVXeM+/yfAQFV9BrjQz+l2AwUuHiwiNwI3AjRo0KCw0Iwp8wrro/D4Y29Wnv3kpDj2ZeWd0d3kb661yp64qDVXd2oYmABNqeGkj+JT4E3gbSAQ96F1gY1e+5uAMws4FhG5BLgAqAy8WtBxqjoK92is9PR0+7PHRL0Yh4mif5s6uXMppt7VjVNrJbN44x4GvjYLgFmrduQe+/cvl3Jp+3p55mWYss9JoshW1TeCHkkBVNXTN2KMCYKL251IFKfWSgagTf3Kua9f9fYveY6/e+wi3ry6Q8jiM+HnpDN7gojcKiK1RaSKZyvBZ24G6nvt13O3lZitmW1M0SXE+r476Na0ms/2qcu2cuBINh//ssE6u6OEk0RxLXA/roKA891bRgk+cx7QVEQaiUgCcCUwvgTny2VrZhtTdAnuFfTqVi6Xp/3Oc5vm2b/Da7/1Y1/zty+WcOrfp/LkxN+CH6QJKyejnhr52Bo7ObmIjAFmA81FZJOI3KCq2cDtwNfAcmCcqi4ryUUYY4rPkyjyD6dNT8374OCeXs18vv/tn9YydenW4ARnIoKTUU/X+GpX1fcLe6+qDiqgfTIwudDoikhE+gP909LSAn1qY8osJ/MtChtBdfOH81n+f72tk7uMcvLo6QyvrRvwODAgiDEVmz16MqboEt13FL6GCl7cri4ADaqWB+Dmc5rkvpbxyHn869LTc/dbPDo1eEGasHLy6Gm41zYMaA9YuXFjygjPUqsxPpZcrV/FlSAGtKkDQPfm1XNfq1YxkcvS6+c5/niOMn7xFjbtPhSscE0YOFoKNZ+DQKNABxIINurJmKLzpAdfs649bZ4k4unPSIo/8atj7TMnSpY3+dtk7hizkK7PTeON6avznOvi12fR+6UZHM22goOljZOlUCeIyHj3Ngn4Hfgi+KEVnT16MqboPEnA16OnnNxE4dr3PKaK9br7EB93IgDPTV3B71v3A/DFwk0s3LCHFVv30+yRKScdu/fQMTbstLuQSOVkwt3zXl9nA+tV1YrzGVNGeH7P5/i4o8hRzzGugzyJIn9yeKRfC56ctPyk91/w0gzWPduPu8cuztN++Ojx3I7v1BGTctufurg1V51pJUIijZM+ih9x3UWkAFVwJQtjTBkh7odPvur9eZKHJy94Jucdz8l78EXuTm+Aj4edyWUd6uXuHzxy8q+MFo9OJXXEpDxJAuDhL5Zy9j+n2XoYEcbJo6e/AnOBS4BLgTki4qjMeKhZH4UxxVfAEyTXa+5kkujum8j/i9x7Nb2zmlTjX5e1yd1v9djXRYpjw65DNH14Cs0emcKvm/YU6b0mOJx0Zt8PtFPV61T1WqADrpLfEcf6KIwpOvXZO+Hi6b/w3FkkuBNCdr47Cs8jKW91fJQvf3OI7xpR3mt4exzNzmHAq7NIHTGJ37bsKzBGE3xO+ih2Avu99ve724wxpUDL2smOjhNOvqXwdGJ7Rj8l+EgI4Ht97hkP9CDt4bwd1229ig16ZD7Vh5gYoXXdZJZu9p0Q+o6c6bP8+bpn+/k83gSWk0SxCvhFRL7CNTBiIPCriNwDoKovBjE+Y0wJTL+vO1UrJvg9xt9aRCfuKFz7BSWKWB8zt/Mvxbr2mb4cPJq3iGDrusm5SaZu5XIFJgrgpCQB0PP56fxwX/cC32MCw8mjp9XAl5wYPfcVsBao5N6MMREqtVoFKiXF+z3GkwwqJp78d6Pke/Tka1KeE5Pv6IaInFQupGNq1dyv61Z2Te5LKeeKt42Pu4/81uw4eFKHuLenJy9n4q9bihGx8VboHYWqlprlR63WkzFFV79KOR7qcyr93bOvvcXkDp0t2Wek1XAVc4jPtzxrh4an5H7tuQHx5KLT66aweOOePMc/0q8FMSKcXi+FS9+cndu+Yeeh3DIjHkPfncu037cDcPvHC+nevDqjh3Ys2YVEKSdFAasDDwCtgNzeKSdLoYaaqk4AJqSnpw8LdyzGlBYiwk1eNZy85U7GK+Fa2Z5HVvlX3fPuSPe8ln3c1Za/38Nff8TZ/5oGwJe3deHpScsZc2On3CThMf337fR4fjrT7FFVkTl59PQRsAJX2Y5/AOtwrSlhjCnjzmriejTUuUnVQo4sHu/846lQ62vin1MXvTaLuet25a7xnd/aHQfZn3Ws2OePVk46s6uq6n9F5E735LsfRcQShTFRID21CplP9cn96754PRQF804K+YfiliRh5Fc+IZZD7o700x7/BvB/h/LZ/E3c9+liXr+qPbd+tCC3fe0zffPMSh+/eAvdm1cnuZB+IICR32dy8Eg2FRPjGJ5vUahI5yRReNLvHyLSD9iCa4a2MSYK+Br6GijeucDztaeOVP7Z3yWRnBTPx8M6cdFrs3Lb3pi+muemrmDpPy7I05Hv3TnunSQAGj3kulNZ+0xfNu0+zB1jFuZ53ZNI9mcdo2JiHCLis7P9hW9Xsvrpvj5Hi0UiJ4niSRFJAe4FXgGSgbuDGpUxJip4D6Ly3EGUT4zj4NHjnFIhgaFdUnl31rpinbtx9Qqs2X4QgOycHNrUyzsR97mpKwDXsq4A/U6vzcgr2zk6tydhOG33xfvxmJP5IKrKyj8PUD4hloz1u+jVshazV+9k2Pt5V6b+9fHzHd3hFIWTRDFbVfcCe4EeAf30ALNRT8aULt6PcTz3D1ed2QBwLZKUFB9b7ERRrUJibqJoWSelwCq3HpN+/YNJv/5RrM8qKc9dx+yHelI7Je/a5dN/38Z17/p62r/YRxv8sHxbntpbgeAkUcwRkUXAu8AULenwhyCyUU/GBFcxp1Hk8dVtXXhy0m/MW7c7T5+Hek3qu7X7iT/2hvdMY1G+YbIAN53dmP/MWFPg57wyuB2fZmzkrLRqNK0RGWut9W9Th4ZVyvPqtFU+X+/8zA+Aq9TJzR/OL9ZnnN2seuEHFZGTRNEMOA+4HhgpIuOA0aq6MuDRGGPKvDb1K1O9UiKQdwJf/kWSPO49v7nP8zzUtwVXndkwd2hsfjWTk7i9Z2A6je+/oDn7s7KpVjEhTzn1lrWTeeHyNvR5eabf9+fvBL+mc0OGvZ/BLd3TfCaEwpJE2/qVWbRxDw2qlOepi1vTvFYlqldMJDtHg9Kn5GTCnQLfAt+KSA/gQ+A2913GCFWd7e/9xhiTX467+Kx3ThjSqSHf/vZn7jrdTjSoWh4R/2VIvF2RXp+xGRsLPW5ol1Qm/voH2/cfOelxkHeimDi8KzEx4qpXJcL6nQeZtXonPZpX5/3Z6zl4JJvhPZue9NirRnISX93eFYB5D5/HGU99V2hMT198GleeUf+kuSje8s98DxQnE+6qAkOAq4E/geHAeKAt8CkRuiyqMSZ4fFWLBd8VY325oVsjpi7bSnrqiZnZ9auUL1bdphn392DLnsN8On8TTWtU5JUfVnHAxxoYAM9dejo1U5IY+X2m33PGx8Yw5c5ubN2bdVKfwc8jenLoaDZpNSrlOR6gcfWKNK7uesz1t74tHMVfvVJibmf26u0HOPeFH3Nfe7hvCwa2rUONZGf/X4NFCutyEJGVwAfAu56V7UTkLlV9SUQeVNXnQhBnkaSnp2tGRkbhBxpjiuyN6avpeWoNmtfKW+rt0NFsYkRIio/NbfN00oayyuv2/UfYefAIp9byXTU3+3gOP63aQVqNinR9zvdjq8WPnZ9bcyqaiMh8VU0/qd1BopD8HdgiskFVGwQ4xoCxRGFMZFi6eS/bDxyhR/Ma4Q7FJ08i++CGjsSKsGVvFpd6rc4XbQpKFE77KE46X0CiCjAbHmtMZGldN7IXEbv7vGZsP5BFt6aBHylUlhR6R+HzTXZHYYwxZU6R7yhEZD/4XCNRgHI+2o0xxpRBBSYKVbVFiYwxxjgqM26MMSaKWaIwxhjjlyUKY4wxflmiMMYY45clCmOMMX5ZojDGGONXsSbcRToR2Q6sB1JwLbjk4b1f0NfVgB0BCCP/Z5fk2IJe93d9vtqi4Zqdfs9LyzU7aSut1+z0e+yrza7Z9zWX9HobqurJ09RVtcxuwKiC9v18nRGMzy7JsQW97u/6ovWai/A9LxXX7KSttF6z0++xXbPzaw7U9ebfyvqjpwl+9gv6OlifXZJjC3rd3/X5aouGa3b6PQ+UYF+zk7bSes1Ov8e+2uyag3/Nucrko6eSEJEM9VHrpCyza44Ods1lX7Cut6zfURTHqHAHEAZ2zdHBrrnsC8r12h2FMcYYv+yOwhhjjF+WKIwxxvhlicIYY4xfligKISIXichbIjJWRM4PdzyhICItRORNEflMRG4JdzyhICIVRCRDRC4MdyyhICLdRWSm+/vcPdzxhIKIxIjIUyLyiohcG+54QkFEurm/x2+LyM/FPU9UJgoReUdEtonI0nztvUXkdxFZJSIjAFT1S1UdBtwMXBGOeAOhiNe8XFVvBi4HuoQj3pIqyvW6PQiMC22UgVXEa1bgAJAEbAp1rIFSxGseCNQDjhEl16yqM90/yxOB94r9ocGYxRfpG3A20B5Y6tUWC6wGGgMJwGKgpdfrLwDtwx17qK4ZGABMAQaHO/ZgXy/QC7gSuA64MNyxh+iaY9yv1wQ+CnfsIbrmEcBN7mM+C3fsobhmr9fHAZWK+5lReUehqjOAXfmaOwKrVHWNqh4FPgEGistzwBRVXRDqWAOlKNfsPn68qvYBrgptpIFRxOvtDnQCBgPDRKRU/lwU5ZpVNcf9+m4gMYRhBlQRv8+bcF0vwPHQRRlYRf1ZFpEGwF5V3V/czyxwzewoVBfY6LW/CTgTGA6cB6SISJqqvhmO4ILE5zW7n1lfgusXyOTQhxU0Pq9XVW8HEJHrgB1ev0TLgoK+x5cAFwCVgVfDEFcwFfSz/DLwioh0A2aEI7AgKuiaAW4A3i3JyS1RFEJVRwIjwx1HKKnqdGB6mMMIOVUdHe4YQkVVPwc+D3ccoaSqh3D90owqqvpYSc9RKm+xg2QzUN9rv567rSyLtmuOtusFu2away4xSxQnzAOaikgjEUnA1bk5PswxBVu0XXO0XS/YNds1B0BUJgoRGQPMBpqLyCYRuUFVs4Hbga+B5cA4VV0WzjgDKdquOdquF+ya7ZqDd81WFNAYY4xfUXlHYYwxxjlLFMYYY/yyRGGMMcYvSxTGGGP8skRhjDHGL0sUxhhj/LJEYYwPIlJVRBa5t60istn99QEReT1In3mXiFxTxPc8LyI9gxGPMR42j8KYQojI48ABVX0+iJ8RByzAVco+2+F7YnGVanhLVaNiUS0THnZHYUwRuFeGm+j++nERec+9Utx6EblERP4pIktEZKqIxLuP6yAiP4rIfBH5WkRq+zh1T2CBqmaLSBMRWeD1mU09+yKyTkSec+9fpqrrgaoiUiv4V2+ilSUKY0qmCa5f8gOAD4FpqnoacBjo504WrwCXqmoH4B3gKR/n6QLMB1DV1cBeEWnrfm0oectE71TV9qr6iXt/AaV0JUJTOliZcWNKZoqqHhORJbhWGZvqbl8CpALNgdbAtyKC+5g/fJynNq4aPR5vA0NF5B5cS/B29HptbL73bgPqlOwyjCmYJQpjSuYIgKrmiMgxPdHpl4Pr50uAZarauZDzHMa1frXH/4DHgB+A+aq60+u1g/nem+R+vzFBYY+ejAmu34HqItIZQETiRaSVj+OWA2meHVXNwlUJ9A0KX52sGbA0MOEaczJLFMYEkXv94kuB50RkMbAIOMvHoVOAs/O1fYTrzuSbgs7v7gNJAzICEa8xvtjwWGMihIh8ATygqpnu/fuAFFX9u5/3XIxrSG2BxxhTUtZHYUzkGIGrUzvTnTQ8I6r8iQNeCHZgJrrZHYUxxhi/rI/CGGOMX5YojDHG+GWJwhhjjF+WKIwxxvhlicIYY4xfliiMMcb49f+muXH0FH+CMwAAAABJRU5ErkJggg==\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": 53,
"id": "abandoned-register",
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAFECAYAAADIlyJZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA120lEQVR4nO3dd3hUdfbH8fcJCEivooAUDUXAAkZAsay9gW137etawC6WdVdXXXV1bdvV3Z+9i4BdUIS1YgOkS5PeERGQUBNSzu+PmYQhJJMbyJ2b8nk9zzyZW2bmXCVzcr/lfM3dERERAUiLOgAREak4lBRERKSQkoKIiBRSUhARkUJKCiIiUkhJQURECikpiIhIoZolHTCzngFen+Pu08sxHhERiZCVNHnNzDYCEwBL8voO7t4+hLhERCQCJd4pABPc/bhkLzazT8s5HhERiVCJdwqFJ5jt6+7LUhSPiIhEKEhH88jQoxARkQohSFKYbGaHhR6JiIhELkjz0fdAOrAE2Eys49nd/aDwwxMRkVQKkhTaFbff3ZeEEpGIiEQm2egjYPuXv5ntBdQJPSIREYlMqX0KZnaGmc0DFgFjgMXAhyHHJSIiEQjS0Xw/0AeY6+4dgOOBcaFGJSIikQiSFHLcfS2QZmZp7v4ZkBFyXCIiEoFS+xSA9WZWH/gSGGxmq4mNQhIRkSomyOijesBWYncVFwGNgMHxuwcREalCSk0KUDgstaO7f2xmdYEa7r4x9OhERCSlgow+Ggi8CTwV39UaeDfEmEREJCJBOpqvA/oCGwDcfR6wV5hBiYhINIIkhWx331awYWY1gdLbnEREpNIJkhTGmNkdwJ5mdiLwBjAi3LBERCQKQUYfpQFXACcRK4Y3GnjWg/RQi4hIpRIkKfQHPnD3/NSEJCIiUQnSfHQeMM/M/mpmXcIOSEREohN0nkJD4ALgMmKdzC8AQzRXQUSkaglyp4C7byA2V2EosA9wNrEV2W4IMTYREUmxIH0KZxC7Q0gHXgZecvfV8ZnNs9y9fehRiohISgQpiPdL4F/u/kXiTnffYmZXhBOWiIhEIVCfgoiIVA8l9imY2fulvTjIOSIiUnmUeKdgZuuBL4o9GD8F6Obu+4UQl4iIRCBZn8KZAV6/rfRTRESkslCfgoiIFAo0T0FERKoHJQURESmkpCAiIoWSTl4zs8OBi4GjiJW32ArMAD4AXnX3zNAjFBGRlEk2JPVDYCXwHjARWA3UAToBxwL9gX+6+/DUhCoiImFLlhSau/uapC8OcI6IiFQeJfYpuPsaM6thZp8lOyecsEREJApJO5rdPQ/IN7NGKYpHREQiFKRK6iZgupl9BGwu2Onug0KLSkREIhEkKbwdf4iISBUXdDnOPYG27j4n/JBERCQqpU5eM7P+wFRgVHz7EDPTMFQRkSooyIzme4FewHoAd58KqFy2iEgVFCQp5BQzczk/jGBERCRaQTqaZ5rZhUANM+sIDAK+CTesYJo3b+7t27ePOgwRkUpl0qRJa9y9RXHHgiSFG4A7gWxgCDAauL/8wtt17du3Z+LEiVGHISJSqZjZkpKOlZoU3H0LcKeZPRLb9I3lGZyIiFQcQUYfHWZm04HviE1im2Zmh4YfWtKY+pvZ05mZKtIqIlKegnQ0Pwdc6+7t3b09cB3wQqhRlcLdR7j7lY0aqfqGiEh5CpIU8tz9y4INd/8KyA0vJBERiUqJfQpm1jP+dIyZPUWsk9mB84DPww9NRERSLVlH8z+KbN+T8Lz02hghis+y7p+enh5lGCIiVU6g2kcVVUZGhmtIqohI2ZjZJHfPKO5YqUNSzawxcAnQPvF8lc4Wkcriszmr+XreGu7q1zXqUCq8IJPXRgLjgOmovIWIFGP1hixmr9rIMZ2KnSQbuctemADAUZ1a0LZpXSYsXse5GftGHFXFFCQp1HH3W0KPREQqlA+n/0Crxnty8L6Nd9jv7rhDWpoV7jvrv1+zMjOLxQ+fnrL4snLyWLl+K/u1qF/iOdOWradDi3qF2799/tvC53948zs+vuVo0vdqEGqclU2QIamvmNlAM9vHzJoWPEKPLAlNXhMJ3zWDJ3Pmf7/eaf+pj35J+p0jd9i3MjOr2PdYsymbIx/5lPmrN+12PJuyc5m05GeGTVjKs18u5JbXp3LcP8awZVvxI+THLljLmf/9moPu/V+J73nCP78gOzePuT/uXKhhc3Yui9dsLuZVu29zdi63vjGN9rd/QPvbP2Dot0vZnF0xRvoHuVPYBvyNWP2jgl5pJ8Ly2e4+AhiRkZExMKoYRKqr71cFr3QzeuYqlv+8lee+WshD5xwEQHZuHjcPm8ofTu5C++b1SnmH7a5/bTKfz/mpcLthndjXV06uQ60dz31nynJuHjYt0Pve8fYM3pq8nOHX9+WgNo0L95/71FhmrtwAwKS7TqBZ/dqFxzK35jBr5QZaNa7D96s28vX8Ndx8Qiea1KtV9O0LZeXkcePQKYye+eNOx25/ezq3vz0dgPvP6s7FvdtiZjudlwpBksLvgHR3XxN2MCJStaTFv9gSBzmOW7iOkdNXsTErl1eu6B3ofbJy8nZICIm8mBHyQRMCwFuTlwNwxn++5j8X9mDl+q3sUSOtMCEAnPLolwy/vi+HP/QptWumkZ27c/fqy2OXMO+BU9mjxvYGmNy8fO4dMZN+B7Xi/KfHBYrnT+/O4E/vzgDgt4e346WxS+jdoSmDju9IRvsm1KqRFmrCCJIU5gNbQotARKoMdy/8wlqxfiurN2TH928/p6Arouho+E9m/8hr45fy3KWH7bB/1IxVXP3qpBI/c1tePs99tYjju+xFk7q1aFR3j12O//rXphS7/6eN2Rz+0KcAxSaEAh3v/BCAm0/oROe9GxTG/eq4pcWe36NtY6YsXV/i+700NlbMdPyidVz07Pgdjo0cdBRdWzUs8bW7KkhS2AxMNbPPiJXPBjQkVUR25g4Ff8T2ffjT7fsT/povuHvId+f1ics48YCWNKlXiytemhh/D9/hL+F7h89M+plPjVnIc18t4v73Z9G68Z58fftx5XU5u+xfH89Nevz3J3fmumO3T751dz6bs5rLXww+7+ru92bw5jVH7HKMJQmSFN6NPyoMzWgWqZhKmgqbn3Cg4Pt+3upN/OHN7zi6UwtO7tZyh3NrJLSObMsr/i/zDVmxjtnETuwV67cGivOL3x/L0X/7LNC55Wna3ScVeydjZhzXpSWLHz4dd2fYhGXU3iON0TN+5LSD9mHD1hzuijcpFXj4lweGEmOQ9RReCuWTd4M6mkUqpnx3apC8vbvgTmHrtjwg1jRz5zvbv/By8/OpkVajcDunhKRQYMzcHfsa8vN3TE37NKrDD5lZDBnYh8P3b1b6RSSYdvdJ1K1do7BZqMCXfziWfZvWZdqy9ezVsDb7NNqTrnePYkv8mhJ9MOhIurUKXtHZzDi/V1sAzu7RpnD/xX3alSn2XRVkRvMiivkDwN0jG30kIhVTSVVzrJjnefEv76IpJL9IDkgrY6fqfnfsOFx27B+PL/a8xQ+fzrNfLuQvH8wmo10TJi75ufDYx7ccQ/pe2+c/PHzOgRzYptFOX+6Jczhm3XcKAHe8M50xc35ixfqtLHrotMhGEe2qIM1HifUx6gC/BiKdpyAi5evnzdvocf9H/PPcgzmnZ5vSX1CC4kYCwfYmI9jelFSwL63IbKm8IpllV+uz/alfV644skPScwYctR8Djtr+923R/owCBX+5B/Hg2eE066RKqZPX3H1twmOFu/8bSN20RREJ3aK1sUlaL48tceneQPLyS0gKCfcD+fEv+YI7ACtyr5Cds2MTTEHfQaLXBsSGss79y6klxlJaQig2zkr2V30YgjQf9UzYTCN25xDkDkNEKphRM1bRqWX9nUpDeOEX9e69/6ffr6bfQa122p/4XZsbTxw1a2wfhZTo0L98zDk9W3P1MftTo0hAj55/CGce0hqgsKTG1cfsz5NjFuxwXirLbVQ1Qb7cE9dVyAUWA+eGEo2IhKpg3HzRL82CP/DL2n5f1MTFPxebFN6duoKHfxmb0VzQEVwrPslr1g8bdjr/7ckreHvyip32FySERLef2oWV67eyNSePZy4pthq0lEGQ0UfHpiIQEdl167dso2aNNOrX3rWb+IIv6t1NCi9+s5h7+nfdqRkmK2d773HB5K/VG3ee2JZMsr/+H7ugRxkjlZIkW47zkmQvdPeXyz+cYDRPQWRHh9z3EQ1q12T6n0/epdcX7fzdHR3+OLLYL/DN2bnUq12TdZu3JX393g3rsGrDjgX2Fjx42u4HJoEk+7PisBL2nwG0BiJLCpqnILKzjbtRZdO97HcKBSN1Ppm9c4G39Vt2/uLvds9oFjx4Gt8sKLmM2vXHpnPryZ0DxyDlr8Sk4O43FDy32L3gRcBtxBbceSD80EQkVQr7FIIU04/bsi2PerVrFpanSFRSDaH9i8whKHDpEe353UmdaFBn1+sWSflI2gBpZjWBS4FbiSWDX7n7nBTEJSIRKGjfDzI34P73Z/HbI9oXbr90ea/CRWy+mh+sqPKnvzsm6SI5knrJ+hSuA24EPgFOcffFqQpKRFKraOXSnLzik8LqhLb+oROWMXTCssLtYzq1oGGdmjvMK3j+0gx+0Wkv0tKM9rd/ULi/RYPaTLjzhHK8Aikvye4UHgdWA0cCfRNGExjg7n5QyLGJSIqY7ThnoKQidOc+NbbY/Rf0iq13/Pa1fTnhn2MK9x/XZXuhu8UPn86StZuZ/cMGTum+T7nELeUvWVIo+3RAEamUit4pbCthzYDFa4tfWqVgVbX0veoXrhEw+U8n7nReu2b1aNcs+GprknrJOpp3b767iFQoyfoJCo4U3Clk5+5c7fO8hLuEhQ+exsOjvuegNo12mqz2zrV9dz9YiYzKVYhUE7kl1CWC2JrDQGGl0GXrtq9LMGvlBprWq8X4ResK96WlGXecdkBIkUqUlBREqomSmoQArnplx+Uub31j+xrHpz325Q7H3r/hyPINTCqUUkclm1k9M0tL2E4zs7rhhlVqTP3N7OnMzMwowxCpVJIlhaKWriu+76DL3g3o3jr4gjFS+QSZqvIJkJgE6gIfhxNOMO4+wt2vbNRI/zhFgippRNHPRcpOfDWv5DkGo246ulxjkoonSPNRHXcvXATV3TdFfacgImVX0p3C5KU/77B98XPjC5/Puu9k7npnBjed0Im2zfRrXx0ESQqbzaynu08GMLNDgWCrY4tIhVHSnUJxZSoK1K1Vk3+ed0hIEUlFFKT56CbgDTP70sy+AoYB14calYiUu+9/2Fj4/Ot4GYq1m7IL9716Re8dzv/81l+kJC6pWIKspzDBzLoABaUL57h7TrhhiUh5+89n8wufX/TseH51aBvenLS8cN+RHZsXPu/WqiHtm2uSWXWUrPbRce7+qZmdU+RQJzPD3d8OOTYRKUfdWjVkdsIqZ4kJocCih07jjYnL6XewylBUV8nuFI4BPgX6F3PMASUFkUqkuCRQ4IvfxxZYNDPOPWzfVIUkFVCyMhf3xJ/e5+6LEo+ZmeoiiVQiicNMWzWqw8rM7dVOT+2+t0YWSaEgo4/eAnoW2fcmcGj5hyNS/fy0MZv3v1vJpUe032lt4/KSOMz069uP46NZPzJm7k/ceHxH9mpYJ5TPlMopWZ9CF6Ab0KhIv0JDQP+KRMrJDUMmM27hOvqmN6dTywbl/v6JI4xeG9gbM+OkbntzUre9y/2zpPJLdqfQGegHNGbHfoWNgNZGFikn67fEBvPllDCPYHcd+pftBQh6d2gWymdI1ZGsT+E94D0zO9zdi19ZQ0Qil6wk9l3vTi98ntGuCTXSwmmekqojyOS1s82soZntYWafmNlPZnZx6JGJVBMF/QgBlkUuVlbO9juMdQl1jJau3cKr45YWbg+5ss+ufYBUK0GSwknuvoFYU9JiIB34fRjBxCuyTjSzfmG8v0hFtLt/u09fsb1a8OOfzgNiC+Ic/bfPCvffcmIn9qgR5Nddqrsg/0r2iP88HXjD3QPXqzaz581stZnNKLL/FDObY2bzzez2hEO3Aa8HfX8RgYU/Fdar5IWvF3PAn0btsCDOOT1bM+j4jlGEJpVQkKQwwsy+JzYE9RMzawFklfKaAi8CpyTuMLMawH+BU4GuwAVm1tXMTgRmAasDvreIACszs0gcybo1Z8elNP957iGpDUgqtSC1j243s78Cme6eZ2abgTODvLm7f2Fm7Yvs7gXMd/eFAGY2NP5+9YF6xBLFVjMb6e7hDMcQqSLcncc+iTUZHdu5BZ/N+anw2D39u3JZX80zlbIpU+2jIhNrdrXMRWtgWcL2cqC3u18f/4xLgTUlJQQzuxK4EqBt27a7GIJIxVHwa7UrHc2JM5NfuKwXo2as4upXJ3H/Wd35TZ925RShVCcVrvaRu79YyvGngacBMjIydnG8hkjFUTM+TDQnv+w3xv3i6yf/5azuAJzSfW8WP3x6+QUn1U7S2kfxtZk/dPfy7PxdASRW3GoT3xeYmfUH+qenp5djWCLRqFUz1rVXljWUAVZlZvFzfOJbt1YNyz0uqZ6SdjTHm3D+UM6fOQHoaGYdzKwWcD4wvCxvoDWapSopGCpa1hnNv0moZ9SjbZNyjUmqryCjjz42s1vNbF8za1rwCPLmZjYEGAt0NrPlZnaFu+cSW7ltNDAbeN3dZ+7yFYhUcgWzjPPyg7eGTli8jnmrY0NRv/zDsaHEJdVTkCqp58V/Xpewz4H9Snuhu19Qwv6RwMgAn10sNR9JVZJmZUsKqzdm8esnt1ee2bepyl5L+Sn1TsHdOxTzKDUhhEnNR1KVFJQjCtqn0OuBTwqff3Wb7hKkfJWaFMysrpndZWZPx7c7qgyFSPkpuEHYsi0v6Xk5eflk/OWjwu0hA/vQponuEqR8BelTeAHYBhwR314B/CW0iESqmbq1agCwdnN2iedk5+bR8c4PWbMpVvDuN33acfj+KoMt5S9IUtjf3f8K5AC4+xZ2v4bXbjGz/mb2dGZm4DJMIhVWQUfz6g3FJ4Wt2/LofNeowu0uezfg7v5dUxKbVD9BksI2M9uTWOcyZrY/UPKfNCmgPgWpStbEV0b7ceOOv1Z5+c7Nw6ZywN3bE8IDZ3dn1E1Hq+KphCbI6KN7gVHAvmY2GOgLXBZmUCLVxdZteUxZuh6AEdNW0rBOTX57RHty8vI5/bGvdjj3ml/sz0W9VbpCwhWkIN7/zGwS0IdYs9GN7r4m9MhEqoEFP20iOzefurVqsGVbHoPHL2Xw+KU7nXfmIa34/UmdI4hQqpsgo48+cfe17v6Bu7/v7mvM7JPSXhcm9SlIVTE/PgHtpct70b31zqUqrj82nQUPnsaj5/cgTUtpSgqUmBTMrE585nJzM2uSMJu5PbFKp5FRn4JUFfNWb6RmmnFwm8YMv+5IRt10FEfs34xaNdN48bLDuPXkzlpXWVIqWfPRVcBNQCtgEttHHG0A/hNuWCJVX1ZOHv/9bAG1aqQVFsXrsndDXhuotZQlOsmqpD4KPGpmN7j74ymMSaRaeGvycgBO6LpXxJGIbBdkXFu+mTUu2Ig3JV0bXkilU5+CVAWTl6ynYZ2aPH5Bz6hDESkUJCkMdPf1BRvu/jMwMLSIAlCfglR2r09YxluTl9OrQ1P1GUiFEiQp1LCEdTjNrAZQK7yQRKq2JWs389fR3wNwyeHtow1GpIggk9dGAcPM7Kn49lXxfSJSRu7OKf/+kq05eVx1zH4c3alF1CGJ7CBIUriNWCK4Jr79EfBsaBGJVFHZuXlc++pktubEqqGeeXCkI7tFihVkRnM+8ET8ISK7aPTMH/nk+9UAjP3jcezTaM+IIxLZWZAZzR3N7E0zm2VmCwseqQguSUwafSSVyuI1mxk0ZAoAf+rXVQlBKqyg6yk8AeQCxwIvA6+GGVRpNPpIKpPcvHwGDY0lhF90bsEVR3aIOCKRkgVJCnu6+yeAufsSd78XOD3csESqhoU/bSL9zg/5bnkmffZryl9/eVDUIYkkFaSjOdvM0oB5ZnY9sZXX6ocblkjlN2vlBu54ZzoA9WvX5G+/Opi9GtaJOCqR5IIkhRuBusAg4H5iTUi/DTMokcrO3TnvqbFszM6lc8sGfHjjUapyKpVC0qQQn6h2nrvfCmxCi+uIlGpbbj7//GguG7NzaVC7Jg+c3V0JQSqNpEnB3fPM7MhUBSNSFTzz5UKeHLMAgG/+eBwN6uwRcUQiwQVpPppiZsOBN4DNBTvd/e3QoiqFmfUH+qenp0cVgshOMrfm8NDI2QydsAyILZyjhCCVTZDRR3WAtcBxQP/4o1+YQZVGQ1KlosnLdx77ZF5hQnjhssM4RiUspBIKMqNZ/QgiSXy3fD0XPTuejVm5nNJtbwYc1YGM9k2jDktkl5SaFMysDnAF0I3YXQMA7n55iHGJVBp/Gz2HjVm57NOoDree3In0vRpEHZLILgvSp/AK8D1wMnAfcBEwO8ygRCqDGSsy6ff4VwBcdcx+/PHUAyKOSGT3BelTSHf3PwGb3f0lYrOZe4cblkjF9+gn8wqf//rQNhFGIlJ+gtwp5MR/rjez7sAqQIvKSrX1xdyfuOT5bwG4oNe+PHj2gSSsQyVSqQVJCk+bWRPgT8BwYiUu7g41KpEKKisnrzAhANxwXEclBKlSgow+KlhQZwywX7jhiFRM7s6TYxYyasYPANx/ZjdOP6gVTetpZVqpWoKMPqoN/BJon3i+u98XXlgiFcvrE5fxyKjYuso92zbm4j7tdIcgVVKQ5qP3gExgEpAdbjjBaEazpNKMFZk8OPJ7DmvfhGcvOYw6tdKUEKTKCpIU2rj7KaFHUgbuPgIYkZGRMTDqWKTq2rotj/OfHsu05bEV/u7u141GdVW2Qqq2IENSvzGzA0OPRKSC+fcncwsTwu2nduHANiqrIlVfiXcKZjYd8Pg5l8XXZc4GDHB31xJSUiUNm7CUF79ZwuwfNtCzbWPevrZv1CGJpEyy5qNIi96JRGHRms3c9tb0wu0bjusYYTQiqZcsKewFNHf3DxN3mtmpwGpgSZiBiUTh3uEzqZFmDB7QmwP2aUijPdWHINVLsqTwCMWvtDYLeIFYKW2RSm/tpmxufWMan835CYA/nNKZPvs1izgqkWgkSwoN3H2nuwF3X2JmzUOMSSRltuXmc+5TY1n281Ya1qlJrw5NueyIDlGHJRKZZEmhSZJjdcs7EJEoDPl2KQt+2swzl2RwwgF7af6BVHvJhqR+bGYPWMJvicXcB3wafmgi4flu+XqeHLOAf308l94dmiohiMQlu1P4HfAsMN/Mpsb3HQxMBAaEHJdIaIZ+u5Tb346NMGrVqA73ndldCUEkrsSk4O6bgQvMbD9iq64BzHT3hSmJTCQEk5as4/a3p9OjbWMeOudAOrdsoIQgkiBIldSFgBKBVHpZOXncNGwqrRrV4dUrelOvdpAqLyLVi34rpMpblZnF4PFL+GjWjyxbt5XXBighiJSkwvxmmNkBwI1Ac+ATd38i4pCkCvh41o/c+uY01m/JoXn9Wgw8qgNHpGtEtUhJgqyn8Bgw1N2/Keubm9nzxMplrHb37gn7TwEeBWoAz7r7w+4+G7jazNKAlwElBdktP27I4trBk+nYsj5vXXME+7eoH3VIIhVekCqpk4C7zGyBmf3dzDLK8P4vAjuU3TazGsB/gVOBrsQ6s7vGj50BfACMLMNniBTr+a8WkZufzxMXHaqEIBJQqUnB3V9y99OAw4A5wCNmNi/Im7v7F8C6Irt7AfPdfaG7bwOGAmfGzx/u7qcCF5XhGkQKuTvrNm9jxLSVPPvVIs7q0Zq2zTTXUiSosvQppANdgHbA7N34zNbAsoTt5UBvM/sFcA5QmyR3CmZ2JXAlQNu2bXcjDKlqRkxbyT3DZ7Ju8zYAuuzdgPvO7F7Kq0QkUZA+hb8CZwMLgGHA/e6+vrwDcffPgc8DnPc08DRARkaGl3ccUjl9NW8NNw+bSrfWjbju2HTaNq1L3/Rm1K1VYcZSiFQKQX5jFgCHu/uacvrMFcC+Cdtt4vtEdsnyn7dwzauTSN+rPi9f3kvlrkV2Q5DJa0+Z2RlmdnR815j4Gsm7agLQ0cw6EEsG5wMXluUNzKw/0D89PX03wpCqID/fue2t78hz55lLMpQQRHZTqR3NZvYQsfkDs+KPQWb2YJA3N7MhwFigs5ktN7Mr3D0XuB4YTaxv4nV3n1mWoN19hLtf2aiR1syt7gaPX8LX89dy1+ld2bepOpRFdleQ5qPTgUPcPR/AzF4CpgB3lPZCd7+ghP0j0bBT2UXuzvQVmcz+YQMPjvyeozu14IJe+5b+QhEpVdBeuMZsH1oa+Z/naj6qvqYuW88f357O7B82ANCmyZ488ssDVdROpJwESQoPAVPM7DPAgKOB20ONqhTxPo0RGRkZA6OMQ1Lrxa8Xcd/7s2jZsA5//eVB9GjbmA7N61GzRpA5mCISRJCO5iFm9jmxyWsAt7n7qlCjEili9MxV/Pn9WRzfpSX/Ou9gGtRRh7JIGIL+iZUGrAHWA50SRiJFwsz6m9nTmZmZUYYhKTJzZSY3DZ3KQW0a858LeyghiIQoyOS1R4DzgJlAfny3A1+EGFdSaj6qPtZv2cZVr0yi0Z578Mwlh1JnjxpRhyRSpQXpUzgL6Ozu2SHHIrKD3Lx8Bg2dyuoN2Qy7qg97NagTdUgiVV6Q5qOFgO7XJaW2bMvl5ten8cXcn7j/rG70aNsk6pBEqoUgdwpbgKlm9glQeLfg7oNCi6oUGpJadeXnO5/PXc09w2eybN1WbjulC+cdpsKHIqkSJCkMjz8qDPUpVE0fz/qRe4bPZMX6rezXoh7DruxD7/2aRR2WSLUSZEjqS6kIRKqvDVk53DdiFm9OWk6XvRvw7/MO4ZTue6tTWSQCQUYfLSI22mgH7r5fKBFJtbJ07RYufHYcP2Rmcf2x6Qw6viO1amoymkhUgjQfJS6/WQf4NdA0nHCkOlmVmcVFz41jU3Yub1x9OD3VmSwSuSDLca5NeKxw938TK5IXGU1eq/zWbMrmomfH8fPmHF6+vJcSgkgFEaT5qGfCZhqxO4dIl7NSR3PltiErh0ue+5YV67fy8uW9OahN46hDEpG4IF/u/0h4ngssBs4NJRqp8tyd378xjbk/buS5Sw+jVwe1RIpUJEFGHx2bikCkenjh68WMnvkjd51+AMd0ahF1OCJSRJCV15qZ2WNmNtnMJpnZo2amweNSZp99v5q/fDCLE7u25IojO0QdjogUI8jYv6HAT8AvgV/Fnw8LM6jSqKO58vlmwRque20yXVs15N/nHaJFcUQqKHPfaQrCjieYzXD37kX2TXf3A0ONLICMjAyfOHFi1GFIEj9kbuVfH83lzUnL2b9FfQYP6M1eDVXYTiRKZjbJ3TOKOxako/l/ZnY+8Hp8+1fA6PIKTqqmDVk5PPH5Ap7/ahHucHnfDgw6oSMNtRaCSIUWJCkMBG4CXo1vpwGbzewqwN29YUixSSWUn+8MmbCUv4+ew89bcjjrkFb87qTO7Nu0btShiUgAQUYfNUhFIFL5zV+9kVvf+I6py9bTZ7+m3HV6V7q3bhR1WCJSBoEmoZlZE6AjsTIXALh7ZCuvScXi7gz5dhn3vT+TurVq8q/zDuasQ1qrM1mkEgoyo3kAcCPQBpgK9AHGAseFGplUCuu3bOP2t6YzauYqjurYnH/8+mB1JItUYkHuFG4EDgPGufuxZtYFeDDcsKQymLB4HYOGTGHNpmzuOK0LA47cj7Q03R2IVGZBkkKWu2eZGWZW292/N7POoUeWhFZei97omau44bUptGpch7euOUL1i0SqiCCT15abWWPgXeAjM3sPWBJmUKVx9xHufmWjRurEjMKoGau4dnBsItq71/VVQhCpQoKMPjo7/vReM/sMaASMCjUqqbA+n7OaG4ZM5sDWjXh1QG/q1460YK6IlLMgtY/+YWZdAdx9jLsPd/dt4YcmFc0389dw1SuT6NSyAS9d3ksJQaQKCtJ8NBt4xszGm9nVZqY2m2po4uJ1DHh5Im2b1uWVK3rTaE/NTBapioKsvPasu/cFLgHaA9+Z2WtmppLa1cR3y9dz2QsTaNmwDoMH9KZpvVpRhyQiIQm0QrqZ1QC6xB9rgGnALWY2NMTYpAKY/cMGfvPctzSqu4eK2YlUA0Emr/0L6Ad8Cjzo7t/GDz1iZnPCDE6iNX/1Jn7z3Hj23KMGrw3oQ6vGe0YdkoiELEhP4XfAXe6+uWCHmS1197ZAr9Aik0gtWbuZi54dB8Dggb1p20wF7USqgyBDUl8oZrfFj2mVmypoxfqtXPjMeLJz8xl6ZR/2b1E/6pBEJEUC9SkUI/nKPCHTymvhiSWEcWzYmsMrl/emy96qjC5SnZR4p2Bmt5R0CIj0T0d3HwGMyMjIGBhlHFXNkrWbufCZ8WzIyuGlK3pxYBuNPhapbpI1HyVbR+HR8g5EojV/9UYufGY8OXn5DBnYR+sgiFRTJSYFd/9zKgOR6MxauYHfPDceM2PolYfTeW+tqyRSXZXYp2Bmd8UX1ynp+HFm1i+csCRVvpq3hvOeGkutmmm8flUfJQSRai5Z89F04H0zywImAz8RW3mtI3AI8DFaV6HScnden7iMO9+ZQfpe9Xn+0sM0D0FEkjYfvQe8Z2Ydgb7APsAG4FXgSnffmpoQpbxtyMrhrndmMHzaSo5Mb87/XdyThnVUy0hEgs1TmAfMS0EskgITF6/jxqFTWbUhi9+d2Ilrj02nhlZLE5E41T6uJnLz8nn80/k8/uk8WjfZkzeuPpyebUvsMhKRakpJoRpYtm4LNw2byqQlP3NOj9b8+cxuNFBzkYgUI9nktQuA/7n72hTGI+XsvakruOudGTjw7/MO4aweraMOSUQqsGR3Cm2BN8xsD+AT4EPgW3ePtMSFBLMpO5e735vB25NX0LNtYx49vwf7NlVROxFJLtnoo0eIlcduAJwAXA48aWazia3RPNrdf0xNmFIWU5et58ahU1i2bguDju/IoOPSqVljV8tciUh1EmT00UbgnfiD+HrNpwIvAyeHGp2USV6+8+SYBfzro7m0bFiHYVcdzmHtm0YdlohUImXuaHb3WcAs4B/lHYyZnQWcDjQEnnP3/5X3Z1RVK9dv5eZhUxm/aB39DtqHB84+UOsoi0iZhd6mYGbPm9lqM5tRZP8pZjbHzOab2e0A7v6uuw8ErgbOCzu2qmLk9B849dEvmb4ik7//+mAev6CHEoKI7JJUNDS/CJySuCO+5vN/iTVDdQUuiDdLFbgrflyS2LItl9ve/I5rB0+mXbO6fDDoKH51aBvMNBlNRHbNLs1TMLP67r4pyLnu/oWZtS+yuxcw390Xxt9vKHBmvBP7YeBDd5+8K7FVF9OXZ3Lj0CksWruZa36xPzef0IlaNdWZLCK7Z1cnr80iNmR1V7UGliVsLwd6AzcQG+nUyMzS3f3Joi80syuBKwHatt2dECqn/Hzn2a8W8rfRc2hWrzaDB/TmiP2bRx2WiFQRFWrlNXd/DHislHOeBp4GyMjIqFZzJtZt3sYtr0/l8zk/cXK3ljx8zkE0qVcr6rBEpApJdqfwIPA3ILeYY7vbTrEC2Ddhu018XyBm1h/on56evpthVB4TFq/jhtemsG7zNu4/qzsX926rvgMRKXfJksJk4F13n1T0gJkN2M3PnQB0NLMOxJLB+cCFQV9cndZozs93nhizgH9+NJd9m+zJ29ceoaUyRSQ0yZLCZUBJdY8ygn6AmQ0BfgE0N7PlwD3u/pyZXQ+MBmoAz7v7zKDvWV2s3ZTNza9P44u5P9HvoH146JwDVchOREKVrMzFnCTHApe3cPcLStg/EhgZ9H0SVYfmo/EL1zJo6BR+3pLDA2d358Jeai4SkfBVyjGM7j7C3a9s1KjqNaPk5zv/+XQeFzwzjrq1avLOtUdwUe92SggikhJaT6ECWbMpm5uHTeXLeWs44+BWPHjOgdSvrf9FIpI6+sapIMYtXMugIVPI3JrDQ+ccyPmH7au7AxFJuUDNR2Z2ceLPqJlZfzN7OjMzM+pQdpt7rLLphc+Mo37tmrx7XV8uUP+BiEQkaJ/CLUV+Rqqq9ClszMrhmlcn8/CH33Nq930YfsORHLBPw6jDEpFqrKzNR/rztZzM/XEjV78yiSXrtnDX6QdwxZEddHcgIpFTn0IEhk9byW1vfke92jV5bUBveu/XLOqQRESASpoUKus8hW25+Tw4cjYvfrOYjHZN+O9FPWnZsE7UYYmIFNI8hRT5cUMWFzwzjhe/WczlfTsw5Mo+SggiUuEEvVOYG/9Z4ixnKdm4hWu5/rUpbNmWy+MX9KD/wa2iDklEpFiBkoK7n5/4U4Jxd575ciGPjJpDu2Z1GTKwNx1bNog6LBGREqlPISSbsnP5/RvT+HDGKk7ptjd/+/VBKmYnIhWe+hRCMO/HjZzxn6/436wfueO0LjxxcU8lBBGpFCrlnUJFNmLaSm576zvq1qrB4AG96aPhpiJSiQRKCmZ2JNDR3V8wsxZAfXdfFG5olUtOXj4Pjfye579exKHtmvDfC3uydyONLhKRyqXUpGBm9xBbVKcz8AKwB/Aq0Dfc0CqP1RuyuO61yUxY/DOXHtGeO047gFo1K2XLnIhUc0HuFM4GehBbnhN3X2lmGkITN37hWq57bQqbs3N59PxDOPOQ1lGHJCKyy4IkhW3u7mbmAGZWL+SYSlURRh/tMNy0aV0GD+hN572VK0WkcgvSxvG6mT0FNDazgcDHwDPhhpVc1KOPMrfmcOUrk3hw5Pec1LUl717fVwlBRKqEUu8U3P3vZnYisIFYv8Ld7v5R6JFVUNOXZ3Lta5P4YX0Wd/frymV926u6qYhUGUFnNH9kZuMLzjezpu6+LtTIKhh3Z/D4pdw3YhbN6tdi2FWHc2i7JlGHJSJSroKMProK+DOQBeQTW1PBgf3CDa3i2Jydy53vTOfdqSs5ulML/n3eITStVyvqsEREyl2QO4Vbge7uvibsYCqieT9u5JrBk1n40yZ+d2Inrjs2nbQ0NReJSNUUJCksALaEHUhF9O6UFfzx7enUq12DV6/ozRHpzaMOSUQkVEGSwh+Bb+J9CtkFO919UGhRlSLsIalZOXnc9/4sXhu/lF7tm/L4hT209oGIVAtBksJTwKfAdGJ9CpFz9xHAiIyMjIHl/d5L127h2tcmMWPFBq46Zj9+f1JnatbQ7GQRqR6CJIU93P2W0COpAP43cxW/e2MaBjxzSQYndm0ZdUgiIikVJCl8aGZXAiPYsfmoygxJzcnL52+j5/D0Fws5sHUj/u+inuzbtG7UYYmIpFyQpHBB/OcfE/ZVmSGpqzKzuGFIrJjdb/q0465+B1C7Zo2owxIRiUSQGc0dUhFIFL6at4Ybh05ha06eitmJiJAkKZjZce7+qZmdU9xxd387vLDClZ/vPP7pfP79yVzSW9TniYt7kr6XaheJiCS7Uzia2Kij/sUcc6DSJoVHRn/PU2MWck6P1vzl7O7UraUF6EREIHlSqAXg7pelKJaUufSI9uzfvD6/zmijYnYiIgmSDcA/JWVRpNg+jfbk3MP2VUIQESki2Z1CDTNrQqwA3k6q0pBUERGJSZYUugCTKD4pRDoktSKsvCYiUhUlSwqz3L1HyiIpgzDLXIiIVGcq6iMiIoWSJYVHUxaFiIhUCCUmBXd/MYVxiIhIBaDmIxERKVRqUjCzFqkIREREomfunvwEs7nAYmAY8La7/5yCuAIxs5+AJUAjIDPhUOJ2Sc+bA+W17nTRz9/V80o6Xtz+KK856PUGOVfXXPL+smxXxmsu6//jotsV+ZrL69910e3yuuZ27l78H/zuXuoD6AX8E1gIvA9cHOR1qXoAT5e0neT5xLA+f1fPK+l4cfujvOag16tr3r1rLst2Zbzmsv4/rkzXXF7/rlNxzUUfgfoU3P1bj62+1gtYB7wU5HUpNCLJdknPw/z8XT2vpOPF7Y/ymsvynrrm4MeTXV9p25Xxmsv6/7jodkW+5vL6d110O6zvsEJBmo8aAmcD5wP7A+8Ar7v7pLCDC5OZTXT3jKjjSCVdc/Wga64ewrrmIDWjpwHvAve5+9jyDiBCT0cdQAR0zdWDrrl6COWag9wpmLu7mdUHcPdNYQQiIiLRC9Kn0M3MpgAzgVlmNsnMuoccl4iIRCBIUngauMXd27l7W+B3VM9bNRGRKi9IUqjn7p8VbLj750C90CISEZHIBEkKC83sT2bWPv64i9h8hSrFzM4ys2fMbJiZnRR1PKlgZgeY2ZNm9qaZXRN1PKliZvXMbKKZ9Ys6llQws1+Y2Zfx/9e/iDqesJlZmpk9YGaPm9lvo44nFczsqPj/32fN7Jvdea8gSeFyoAXwdvzRIr6vwjOz581stZnNKLL/FDObY2bzzex2AHd/190HAlcD50URb3ko4zXPdvergXOBvlHEWx7Kcs1xtwGvpzbK8lXGa3ZgE1AHWJ7qWMtDGa/3TKANkEMlvV4o8+/yl/Hf5ffZ3XlkYcyIqygP4GigJzAjYV8NYAGxleNqERty2zXh+D+AnlHHnqprBs4APgQujDr2VFwzcCKxOTeXAv2ijj1F15wWP94SGBx17Cm43tuBq+LnvBl17Km45oTjrwMNdudzS52nYGadgFuB9iTMa3D340p7bdTc/Qsza19kdy9gvrsvBDCzocCZZjYbeBj40N0npzbS8lOWaya2ut5wYLiZfQC8ltJgy0kZr7k+sT6xrsBWMxvp7vmpjLc8lOWa3X1W/PjPQO3URVl+yvj/eBmwLX5OXsqCLGdl/V02s7ZAprtv3J3PDTJ57Q3gSeBZKvF/4AStif2jKbAc6A3cAJwANDKzdHd/MorgQlLsNcfbl88h9kUxMvVhharYa3b36wHM7FJgTWVMCEmU9P/5HOBkoDHwnwjiCktJv8uPAo+b2VHAF1EEFqKSrhngCuCF3f2AIEkh192f2N0Pqujc/THgsajjSCWPjST7POIwIuHVaBEpdy/oD6wW3H0LsS/IasXd7ymP9wnS0TzCzK41s33MrGnBozw+PCIrgH0TttvE91VlumZdc1VU3a4XUnDNQZLCb4HfA98Ak+KPieUZRIpNADqaWQczq0Ws03F4xDGFTdesa66Kqtv1QiquOeoe9pB774cAP7B9aNoV8f2nAXOJ9eLfGXWcumZds65Z11tRrjlIQbxLSkgmLyd9oYiIVDpBOpoPS3heBzgemAwoKYiIVDGl3ins9AKzxsBQdz8llIhERCQygZbjLGIz0KG8AxERkegFmdE8gljtFIhNsT6ASl43RkREiheko/mYhM1cYIm7V9oiUyIiUrJSm4/cfQwwB2gENCWWGEREpAoqNSmY2QDgW2I1cn4FjDOzSlE6W6Q0ZtbMzKbGH6vMbEX8+SYz+7+QPvOmkoZ6J3nN382swhehlMovSPPRHOAId18b324GfOPunVMQn0jKmNm9wCZ3/3uIn1GT2JDunu4e6K7bzGoQK2fwjLtXiwWgJDpBRh+tBRJLsW6M7xOpsuKrlb0ff36vmb0UX71siZmdY2Z/NbPpZjbKzPaIn3eomY0xs0lmNtrM9inmrY8DJrt7rpntb2aTEz6zY8G2mS02s0fi27929yVAMzPbO/yrl+osSFKYD4yP/2LcA4wD5prZLWZ2S7jhiVQY+xP7Qj8DeBX4zN0PBLYCp8cTw+PAr9z9UOB54IFi3qcvsfphuPsCINPMDokfu4wdSx+vdfee7j40vj2ZSrxCnlQOQWY0L4g/CrwX/9mg/MMRqbA+dPccM5tObGj2qPj+6cQWoOoMdAc+MjPi5/xQzPvsA8xO2H4WuCz+B9Z5xBZRKTCsyGtXA6127zJEkis1Kbj7n1MRiEgFlw3g7vlmluPbO+Pyif0eGTDT3Q8v5X22EisXU+At4B7gU2BSQd9d3OYir60Tf71IaIJMXmsB/AHoRsI/Zq8Ey3GKpNAcoIWZHe7uY+PNSZ3cfWaR82YD6QUb7p5lZqOBJyh9YZhOxFZCFAlNkD6FwcD3xEpb/BlYTKymt4jEufs2YkO2HzGzacBU4IhiTv2Q2ILsiQYTu+P4X0nvH08y6VTutUykEggyJHWSux9qZt+5+0HxfRPc/bCkLxSRYpnZO8Af3H1efPtWoJG7/ynJa84mNoy1xHNEykOQjuac+M8fzOx0YCWxmc0ismtuJ9bhPC+eIApGNiVTE/hH2IGJBLlT6Ad8SWxd0MeBhsCf3b2qL3snIlLtBEkKzYqMiBARkSoqSEfzODN7w8xOs/gAbBERqZqCJIVOwNPAb4i1gT5oZp3CDUtERKJQpuU4zexYYlP86xMbcne7u48NJzQREUm1QH0KwMXE7hR+BJ4DhgOHAG+4u5bmFBGpIoI0H40lNuLoLHc/3d3fBq5394nAk6FGJyIiKRXkTsG8yElmttTd24YamYiIpFyQ5TiLyxoahSQiUgUFaT4qTvDeaRERqTRKLHNhZhsp/svfgD1Di0hERCJTpiGpIiJSte1q85GIiFRBSgoiIlJISUFERAopKYiISCElBRERKaSkICIihf4fVcLuvcdLojoAAAAASUVORK5CYII=\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./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": "injured-willow",
"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