Skip to content

Instantly share code, notes, and snippets.

@hannorein
Created April 26, 2021 00:30
Show Gist options
  • Save hannorein/29052d430ab724d6adea25a8fc4c4a81 to your computer and use it in GitHub Desktop.
Save hannorein/29052d430ab724d6adea25a8fc4c4a81 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 42,
"id": "f952f306",
"metadata": {},
"outputs": [],
"source": [
"import rebound\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 40,
"id": "d89d4b71",
"metadata": {},
"outputs": [],
"source": [
"sim = rebound.Simulation()\n",
"rebound.data.add_outer_solar_system(sim)\n",
"\n",
"sim.add(a=55, e=0.1) # add test particle\n",
"i = sim.N - 1 # index of test particle \n",
"\n",
"var = sim.add_variation(testparticle=i) # create a new variational particle linked to the testparticle\n",
"var.particles[0].x = 1. # initial perturbation to the variational particle\n",
" # we use index 0 here because there is only one variational particle\n",
"\n",
"times = np.linspace(100,1e6,100)\n",
"lyapunov_exponent = np.zeros(len(times))\n",
"for i, time in enumerate(times):\n",
" sim.integrate(time, exact_finish_time=0) \n",
" \n",
" # Calculate divergence of variational particle\n",
" d = np.sqrt(np.sum(np.square(var.particles[0].xyz + var.particles[0].vxyz)))\n",
" lyapunov_exponent[i] = 1./time * np.log(d/1.) # divide by initial perturbation (here 1.)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "89c0faea",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAHlCAYAAABmsfzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABPFUlEQVR4nO3dd3ykdbn//9eV3nuyJbvZkmR3WXapobcFRBFBFJGfCCrKga9w7Md2jp5jOyp6PBz7AfTYFQQBBQUVpfctdLZle7al92TTrt8fM1nDkuzOZmYyk5n38/GYR3LfM3PPlZslVz7t+pi7IyIikqxSYh2AiIhILCkRiohIUlMiFBGRpKZEKCIiSU2JUEREkpoSoYiIJDUlQhERSWpKhCIiktTiPhGa2RFmdpOZ/c7Mrot1PCIikliimgjN7Cdm1mhmLx9w/nwzW29m9Wb22YNdw93XuvsHgcuA06IZr4iIJB+LZok1MzsT6AZ+4e7LgudSgQ3AeUADsBK4HEgFvn7AJT7g7o1m9lbgOuCX7v6bqAUsIiJJJ6qJEMDM5gN/HJMITwG+6O5vCh7/K4C7H5gEx7vWn9z9LRM8dy1wLUBubu7xS5YsicwPICIiCWH16tXN7l5+4Pm0GMRSCewYc9wAnDTRi81sBXAJkAncN9Hr3P0W4BaAuro6X7VqVQRCFRGRRGFm28Y7H4tEeFjc/WHg4RiHISIiCSoWs0Z3AnPHHM8JnhMREZlysUiEK4FaM1tgZhnAu4B7InFhM7vIzG7p6OiIxOVERCQJRHv5xK3AU8BiM2sws6vdfQj4EPAXYC1wu7u/EonPc/d73f3awsLCSFxORESSQFTHCN398gnO38dBJr6IiIhMlbivLCMiIhJNSoTj2NvZz3k3PsI9L+yKdSgiIhJlCZUIIzVZpjgng83NPWzY0xWhyEREJF4lVCKM1GSZjLQU5pXksKmpO0KRiYhIvEqoRBhJC8vzlAhFRJKAEuEEairy2Nrcy9DwSKxDERGRKFIinEB1eS4DwyM0tPXFOhQREYkiJcIJVFfkAah7VEQkwSVUIoxkibXqMiVCEZFkkFCJMJIl1gpz0inLy2RTY08EIhMRkXiVUIkw0qrLc9UiFBFJcEqEB1FdkUd9UzfuHutQREQkSpQID6K6PI/23kFaewZiHYqIiESJEuFBVJfnArCpSeOEIiKJasJtmMzsuyG8v9PdPx/BeMJiZhcBF9XU1ETketXl/5g5euKCkohcU0RE4svBWoQXA6sP8XhHtAM8HJHemLeyKJvMtBQ2NWrCjIhIojrYxrz/4+4/P9ibzaw4wvHElZQUU81REZEEN2GL0N2/fag3h/Ka6a6mIk9jhCIiCeyQk2XM7JtmVmBm6Wb2dzNrMrMrpyK4eFBdnsuOtl76B4djHYqIiERBKLNG3+juncCFwFagBvhUNIOKJ9XlebjD1ha1CkVEElEoiTA9+PUtwB3uHn4hz2lk/8xRlVoTEUlIoSTCe81sHXA88HczKwf6oxvW5ESy6PaoBWW5mKn4tohIogolEX4BOBWoc/dBoBd4a1SjmqRIL58AyM5IpbIoW4lQRCRBhZIIn3L3VncfBnD3HuD+6IYVX6q1hEJEJGEdrLLMTKASyDazYwELPlUA5ExBbHGjujyPZ7e0MjLipKTYod8gIiLTxsEW1L8JuAqYA9w45nwX8G9RjCnu1FTk0Tc4zK6OPuYUJ9XfACIiCW/CRBisKvNzM3uHu985hTHFnZqKwMzRjY3dSoQiIgnmYF2jV7r7r4D5ZvaJA5939xvHeVtCqq0YXULRzdmLK2IcjYiIRNLBukZzg1/zpiKQeFacm0FZXgYb92rCjIhIojlY1+jNwa9fmrpw4ld1eR4bG7tiHYaIiETYwVqEAAQX0F8DzB/7enf/QPTCmpxI70c4Vu2MPO55fhfujplmjoqIJIpQ1hH+ASgE/gb8acwj7kRjQf2omvI8OvuHaOraF/Fri4hI7ByyRQjkuPtnoh5JnKudkQ9AfWM3FQVZMY5GREQiJZQW4R/N7IKoRxLnascsoRARkcQRSiL8KIFk2GdmnWbWZWad0Q4s3pTnZ5KflaYJMyIiCeaQXaPunj8VgcQ7M6O2Io96tQhFRBJKKLNGzxzvvLs/Gvlw4ltNRR4PrmuMdRgiIhJBoUyWGbsbfRZwIrAaOCcqEcWx2op8bl/VQFvPAMW5GbEOR0REIiCUrtGLxh6b2Vzg29EKKJ7VzAhMmKlv6uaE3JIYRyMiIpEQymSZAzUAR0Q6kOmgpjw4c1Sl1kREEkYoY4TfAzx4mAIcA6yJYkxxq7Iom+z0VE2YERFJIKGMEa4a8/0QcKu7PxGleMISzRJrACkpRnVFrpZQiIgkkFDGCH8+FYFEgrvfC9xbV1d3TbQ+o7Yin6c3t0Tr8iIiMsUmM0aY1Goq8tjd0U9X/2CsQxERkQhQIjxMo7vVb2rqiXEkIiISCSEnQjPLiWYg08X+mqN7NU4oIpIIDpkIzexUM3sVWBc8PtrMfhj1yOJUVUkOGakpKr4tIpIgQmkR/g/wJqAFwN1fAMYtu5YM0lJTOLKygJVbW2MdioiIREBIXaPuvuOAU8NRiGXaOL2mjBd2tNPRpwkzIiLTXSiJcIeZnQq4maWb2SeBtVGOK66dXlPGiKNlFCIiCSCURPhB4J+BSmAngcoy/xzFmOLesVXFZKen8kR9c6xDERGRMIWyoL4ZuGIKYpk2MtJSOGlhCY8rEYqITHsTJsIDaoy+jrt/JCoRTROn15Txn39ay672PmYXZcc6HBERmaSDdY2uIrDvYBZwHLAx+DgGSPrN+E6vLQNQq1BEZJqbsEU4WmPUzK4DTnf3oeDxTcBjUxNe/Fo8I5+yvAyeqG/msrq5sQ5HREQmKZTJMsVAwZjjvOC5pGZmnFZTxhP1zbhP2IMsIiJxLpREeAPwnJn9zMx+TmAvwq9FN6zp4fSaMpq7B1i3R+XWRESmq1Bmjf7UzO4HTiIweeYz7r4n6pFNA6fVBMYJn6hv5ohZBYd4tYiIxKNQi26fCJxBoLTaCdELJzxmdpGZ3dLR0TElnze7KJuF5bmaMCMiMo2FUnT7BuCjwKvBx0fMLC67Rt39Xne/trCwcMo+8/SaMp7Z3MrA0MiUfaaIiEROKC3CC4Dz3P0n7v4T4HzgwuiGNX2ctaicvsFhVZkREZmmQu0aLRrz/dQ1t6aBM2rLKcpJ567ndsY6FBERmYRDTpYBvk5g1uhDgBEYJ/xsVKOaRjLSUnjr0bP57coddPYPUpCVHuuQRETkMByyRejutwInA3cBdwKnuPtvox3YdHLJcXPYNzTC/S/tjnUoIiJymEKZLPN2oNfd73H3e4B+M3tb1CObRo6eU8jC8lzuXKPuURGR6SaUMcIvuPv+9Qju3g58IWoRTUNmxjuOm8OzW1rZ0dob63BEROQwhJIIx3tNKGOLSeVtx1YCcLcmzYiITCuhJMJVZnajmVUHHzcS2JVCxqgsyuaUhaXctaZBtUdFRKaRUBLhh4EB4LfAbUA/Sb5D/UQuOa6SrS29rNneHutQREQkRKHMGu1x98+6e527n+Du/+buPVMR3HTz5uWzyEpP4c41DbEORUREQjRhIjSzLx7qzaG8JpnkZaZx4VGzuXvNTtp6BmIdjoiIhOBgk17+ycw6D/K8Ae8CvhjRiKa5a89cyO9WN/Dzp7bysTcsinU4IiJyCAfrGv0RkH+QR17wNTLGohn5vOGIGfzsya30DgzFOhwRETmECVuE7v6lqQwkkVx/djWX/HAvtz67g6tPXxDrcERE5CBCLboth+G4qmJOWlDCjx/brO2ZRETinBJhlFy3oprdHf384XktsBcRiWeh1Bo9LZRz8lpnLSrniFkF3PTIJkZGtMBeRCRehdIi/F6I52QMM+O6FdVsaurhvpe1K4WISLyacLKMmZ0CnAqUm9knxjxVAKRGO7BE8Jbls/jhQ/XccP863nDEDLLSddtEROLNwVqEGQSWSKTx2mUTncCl0Q9t+ktNMf7jwqU0tPXxf49viXU4IiIyjoMtn3gEeMTMfubu26YwpoRyak0Z5y2dwQ8fquedx8+hoiAr1iGJiMgYoYwRZprZLWb2VzN7cPQR9cgSyOcuOIKB4RG+9df1sQ5FREQOEEoivAN4Dvg88KkxjyljZrlmtsrMLpzKz42U+WW5vP+0BdyxuoGXd3Yc+g0iIjJlQkmEQ+7+v+7+rLuvHn2EcnEz+4mZNZrZywecP9/M1ptZvZl9NoRLfQa4PZTPjFcfOqeGkpwMvnTvK1pOISISR0JJhPea2fVmNsvMSkYfIV7/Z8D5Y0+YWSrwA+DNwFLgcjNbambLzeyPBzwqzOw84FWgMfQfK/4UZKXzmfOXsHJrG796RkOuIiLx4mC7T4x6X/Dr2O5QBxYe6o3u/qiZzT/g9IlAvbtvBjCz24CL3f3rwOu6Ps1sBZBLIGn2mdl97v66umVmdi1wLUBVVdWhQouJd9bN4Y8v7ebr961jxaIKqkpzYh2SiEjSC2Vj3gXjPA6ZBA+iEtgx5rgheG6iz/+cu38M+A3wo/GSYPB1twQ3D64rLy8PI7zoMTNuuGQ5aSnGp373grpIRUTiQCgl1nLM7PNmdkvwuDYWk1bc/Wfu/sep/txIm12UzecvPIJntrTyy6fVRSoiEmuhjBH+FBggUGUGYCfwn2F85k5g7pjjOcFzSeOyurmctaicG+5fx7aWnliHIyKS1EJJhNXu/k1gEMDdewnsTj9ZK4FaM1tgZhkEdrm/J4zrTTtmxg3vWE5aqvGRW59j39BwrEMSEUlaoSTCATPLJjBBBjOrBvaFcnEzuxV4ClhsZg1mdrW7DwEfAv4CrAVud/dXJhX96z/vIjO7paMj/tfqzSrM5r8uPZoXGjr46p/WxjocEZGkZe4Hn7ARXL7weQKzNv8KnAZc5e4PRz26Saqrq/NVq1bFOoyQfPVPr/Kjx7bwnXcdw8XHTDhnSEREwmRmq9297sDzh1w+4e4PmNka4GQCXaIfdffmKMSYlD59/hKe297Ov971EktnFVA7Iz/WIYmIJJVQd6jPAtoI7Dyx1MzOjF5IySU9NYXvv/s4stNTue7Xa+jsH4x1SCIiSSWU5RPfAJ4APsc/6ox+MspxTcp0GiMca2ZhFt+7/Fi2Nvdw/a/WMDg87lJJERGJglDGCNcDR7l7SBNk4sF0GiMc6/ZVO/j0717k0uPn8F+XHoVZOJNzRURkrEmPEQKbgXRCnCkqk3dZ3Vx2tvXxnb9vpLIom4+ftyjWIYmIJLxQEmEv8LyZ/Z0xydDdPxK1qJLYx95Qy872QDKcXZTF/3dCfNZNFRFJFKEkwntIsgXvsWRmfP2S5TR27eOzd71EakoKlx4/J9ZhiYgkrFCWT/w8WAFmtJ9uvbtramMUpaemcMt7jueffr6KT/3uBdydd9bNPfQbRUTksIUya3QFsJHAHoI/BDbE6/KJ6TprdDxZ6an8+H11nF5TxqfvfJE7Vu049JtEROSwhbKO8L+BN7r7We5+JvAm4H+iG9bkuPu97n5tYWFhrEOJiKz0VH703n8kw58/uTXWIYmIJJxQEmG6u68fPXD3DQRmkcoUGE2G5y6ZwRfueYWv37dW+xiKiERQKIlwlZn92MxWBB8/AqbfIr1pLCs9lZvfczxXnlzFzY9u5qO/fV47VoiIREgos0avA/4ZGF0u8RiBsUKZQqkpxlcuXkZlUQ7f+PM69nb088Mrj6MsLzPWoYmITGuHrCwDEJw1egQwQmDW6EC0AwvHdK0sE6o/PL+TT//uRUpyM7jpyuM5em5RrEMSEYl7E1WWCWXW6FuATcB3gO8D9Wb25siHGL5EmjV6MBcfU8md151KihnvvPkpbl+pGaUiIpMVSq3RdcCF7l4fPK4G/uTuS6YgvklJ9BbhqNaeAT586xqeqG/hsro5fPGtR5KTEUpvt4hI8pl0ixDoGk2CQZuBrohFJpNWkpvBz99/Iv98djV3rG7gou89zqu7OmMdlojItBLqrNH7zOwqM3sfcC+w0swuMbNLohyfHEJaagqfetMSfnX1SXT1D/G2HzzBT5/YoiUWIiIhCiURZgF7gbOAFUATkA1cBFwYtcjksJxWU8b9Hz2D02vL+NK9r3LFj59hR2tvrMMSEYl7Ic0afd2bzDLieeZosowRjsfduW3lDr76p7WMuPOvFxzBFSdWkZKivQ1FJLmFM2v0YTObP+b4BGBlZMOTSDEzLj+xir98/EyOqyrm33//Mpf/6Gk27NWwrojIeELpGv068Gczu97MvgrcArw/umFNTrIsnwhFZVE2v7z6RG64ZDnr93ZxwXce46t/epXufUOxDk1EJK6EuqB+BfAA0Awc6+57ohtWeJK5a3Q8rT0DfPPP67ht5Q5mFGTyqTct4ZJjK9VdKiJJJZyu0X8HvgecCXwReDi4yF6miZLcDG54x1Hcdf2pzCjI4pN3vMCF33ucJ+ubYx2aiEjMhdI1Wgqc6O5PufvNBLZh+lhUo5KoOK6qmN9ffxrfedcxdPQN8u4fP8MHfraS9Xs0figiySvUrtFsoGrsdkzxTF2jh9Y/OMxPn9jKDx+up3vfEO84bg4fP28RlUXZsQ5NRCQqwukavQh4Hvhz8PgYM7sn4hHKlMpKT+W6FdU89umz+afTF3DP87s4+1sP87X71tLeG7crY0REIi6UrtEvAicC7QDu/jywMGoRyZQqysngc29ZykOfWsFFR83mR49t5sxvPsRNj2yif1B7HopI4gslEQ66+4HrEUaiEYzETmVRNv992dHc/9EzOH5eMTfcv46zv/Uwv125naFh/ecWkcQVSiJ8xczeDaSaWa2ZfQ94MspxSYwsmVnAT99/IrdeczIzCrL4zJ0v8cZvP8p9L+1mMlWIRETiXSiJ8MPAkcA+4DdAB3E6a1QL6iPnlOpS7r7+VG5+z/GkmnH9r9dw4fce52+v7lVCFJGEMqlao/FOs0Yja3jEufu5nXz37xvZ3trLUXMK+fgbFrFicTlmWpQvItPDRLNGlQglZIPDI9y9ZifffXAjDW19HD23iI+9oZYVi5QQRST+KRFKxAwMjXDnmga+/2A9O9v7OHpOIR9TC1FE4tykE6GZlbt7U9QiiwIlwqkxMDTCXWsa+P5D9TS09XFMsIV4llqIIhKHwkmEG4CtwG+Bu9y9LSoRRpAS4dQ6sIV4bFURHz1XCVFE4ktYXaNmdiLwLuBtwKvAbe7+q0gHGSlKhLExMDTCHat38IMH69nV0c/Rc4v46Lk1nL24QglRRGIuImOEZlYG3Ahc4e6pEYwvopQIY2u0hfiDYJfpssoC/t+Z1bx52UzSUkNZsSMiEnnh1BotMLP3mdn9BBbS7yZQck1kXBlpKVx+YhUPfXIF33zHUfTuG+bDtz7HOf/9CL98ait9AyrdJiLxI5Qxwi3A74Hb3f2pqQgqXGoRxpeREeeBtXu56ZFNPLe9nYKsNP6/E+by3lPmM7ckJ9bhiUiSCGeyjLm7m1kegLt3RynGiFEijE/uzuptbfz0ya38+eU9jLhz7pIZXHFyFWfWlpOaonFEEYmeiRJhWgjvPdLMfgmUBK5jTcD73P3lSAcpic3MqJtfQt38EvZ09PPrZ7Zx67Pb+dvavVQWZXP5iXO5rG4uFQVZsQ5VRJJIKC3CJ4HPuftDweMVwNfc/dSoR3eYgnsnXlRTU3PNxo0bYx2OhGBgaIQHXt3Lr5/ZxpObWkhNMc5eXM5ldXM5e0kF6ZpcIyIREk7X6AvufvShzsUTdY1OT1uae7h91Q7uXN1AY9c+yvIyeOvRlbz92EqWVRZoCYaIhCWcRHg3sAb4ZfDUlcDx7v72iEcZIUqE09vQ8AiPbGjijlUNPLiukYHhEWoq8rj46Nm8efksairyYh2iiExD4STCYuBLwOnBU48BX4znCjNKhImjvXeA+17aw93PNbBya+Cf3KIZebx52SzOWzqDI2erpSgioVHRbZn29nT08+eXd3Pfy3tYubUVd6jIz+TsxRWcvaScU2vKKMhKj3WYIhKnwmkRLgI+CcxnzCxTdz8nwjFGjBJh4mvu3sfD65t4aF0jj25oomvfEKkpxjFzizijtowzass4ak6RJtuIyH5hTZYBbgJWA/tLgrj76kgHGSlKhMllcHiE1dvaeHxjM4/VN/NiQzvukJeZxskLSzitpozTa8qoqchTN6pIEgsnEa529+OjFlkUKBEmt/beAZ7a1MLj9c08Xt/MtpZeINCNelpNGadUl3LyglLmlmQrMYokkXAS4ReBRuBuYN/oeXdvjXCMEaNEKGPtaO3lyU3NPF7fwpP1zbT0DAAwqzCLkxeWcsL8Ek6YX0x1eR4pqm4jkrDCSYRbxjnt7r4wUsFFmhKhTGRkxKlv6uaZzS08vbmVZ7a00twd+PuuOCeduvklnLywlJMXlnDEzAIlRpEEolmjIuNwd7a29LJyayurtgYS42hXamF2OqfVlHLWonLOXFTOrMLsGEcrIuGYdK1RM3vveOfd/ReRCEwklsyMBWW5LCjL5bK6uQDsau/j6c0tPLWphcc2NnPfS3sAWDwjnxVLyjl3yQyOqyrS3ooiCSKUrtHvjTnMAs4F1rj7pdEMLBxqEUqkuDsb9nbzyIZGHl7fxLNbWhkacQqz0zmjtowViys4a1E55fmZsQ5VRA4hYl2jZlYE3Obu50cotohTIpRo6ewf5PGNzTy4LpAYR8cXl1UWsGJRBSsWl3PMXLUWReJRJBNhOvCyuy+OVHCRpkQoU2FkxHl1dycPrw8kxTXb2xhxKMhK44zaclYsLuesxeVU5GtbKZF4EM4Y4b3AaLZMBY4Abo9seCLTT0qKsayykGWVhXzonFo6egd5YlPz/sT4p5d2A3DUnEJWLK7g7MXlHDWnSBsQi8SZUMYIzxpzOARsc/eGqEY1SdqPUOKF+2hrsYkH1zXyXLC1WJKbwZm1ZZxRW87ptWXM0CbEIlMmrK5RM5sJnEigZbjS3fdEPsTIUdeoxJu2ngEe3djEw+ubeHRD0/5F/Ytm5HFqdRknLSjhhAUllOVp0o1ItISzoP6fgP8AHgQMOAv4srv/JBqBRoISocSzkRFn7Z5OHt8YKAG3cmsr/YMjACwsz6VuXjHHVRVz3LxialTtRiRiwkmE64FT3b0leFwKPKnJMiKRMTA0wsu7Onh2SyvPbmllzfY22nsHAcjPTOOouYUcM7eIo+cUcfTcInWnikzSpCfLAC1A15jjruA5EYmAjLSUQAuwqpgPnlWNu7OluYfV29p4fkc7LzS0c/MjmxkaCfzRWp6fyVGVhRxZWcjSWfksnVXInOJstRxFJimURFgPPGNmfyAwRngx8KKZfQLA3W+MYnwiScfMWFiex8LyPN4ZrHbTPzjMyzs7eGlnBy81dPDizg4eXN/IaIdOXmYaNRV51FTkUV0e+LqgLIe5JTlkpqXG8KcRiX+hJMJNwceoPwS/5kc+HBEZT1Z6KnXzS6ibX7L/XO/AEOv3dLF2dxdrd3dS39jNIxua+N3qf0zqNoPZhdn7y8gtKMtlYXku1eV5VBapFSkCKrotknA6+gbZ3NTNtpZetjT3sLWlh63NPWxu6qFr39D+12Wnp7KwPJfaYEuypiKfmoo85pXmkK7KOJKAwllQXw58GjiSQK1RANz9nIhGKCIRUZidzrFVxRxbVfya8+5OS88Am5t6qG/sDjyaunl2Syu/f37X/telpxrzSnOpCXaxLp6ZzxGz8plfmqvScZKQQuka/TXwW+BC4IPA+4CmaAYlIpFnZpTlZVKWl8mJC0pe81zPviE2NXWzcW8gOdY3drOhsYsH1u5lODhJJyMthSUz8zl6ThFHzQnMZNVmxpIIQlk+sdrdjzezF939qOC5le5+wpREOAnqGhWJjH1Dw9Q3dgfHIjt5eWcnL+3soDvYxZqflcaxVcUcX1XMsVWB5R2F2ekxjlpkfOEsnxgMft1tZm8BdgElB3m9iCSIzLRUjpxdyJGzC/efGx5xNjd18/yOdtZsb2fNtja+/fcN+2ewLizP5Zhgq3H5nCKWziogO0MzVyV+hdIivBB4DJgLfA8oAL7k7vdEP7zJUYtQZGp19g/y4o4OXmho57nt7Ty/o33/FlWpKUZtRR7LKws5am4RR1UWsmRWvpZ1yJQLp7JM6WhVmelCiVAkttydvZ37eLGhnZd2dvBiQ2ANZGuwxmp6qrF4Zj7LKwMtx6WzClg8M5+sdCVHiZ5wEuFG4Hngp8D9Pg3WWygRisQfd2dne9/+pPhS8GtHX2D0JcWgujyPJbMKqK3ICzxm5FFVkktGmmarSvjCSYQGvAH4AHACgb0If+buG6IRaCQoEYpMD+7OjtY+Xt3dwau7OnllVyfr93bR0Na3/zUpBrOLsplfmsu80hyqSgIVc6pKcphbnENBdhqBX1MiBxeRHerN7GzgV0AegVbiZ939qUgFGSlKhCLTW+/AEJsae9jY2MXWll62tfTs/zpakHxUfmYalcXZzCnOYU5xNrOLsqgsymFmYSYV+VlUFGRqPFKA8BbUlwJXAu8B9gIfBu4BjgHuABZENFIRSXo5GWksn1PI8jmFr3uuq3+QHa19bG/tpaGtl4a2Pna09rKjtZenN7fsX9oxVlFOOmV5mZTmZlCWn0l5XiYVBZnMyM9iRkEWs4uymF2UrTHKJBXK8omngF8Cbxvdmd7MPubu3zazm6IanYjIAfKz0lk6O52lswvGfb6zf5CdbX3s6eynsbOfxs597O3qp6V7gJbuAdbu6uTRrn2vKTc3qiI/kznF2cwrzWVuSQ7zSnIC3bGlOZTnZaoLNkGFNEZ44AQZM9vu7lVRjSwM6hoVkUPpHRiisXMfezr72d3RR0NrHw1tgZbm9tZednX0MfY3X05GKlUlOSwsHy1gHihBV1uRR25mKG0KibVJd41OMEtUfxaJyLSWk5HG/LI05pfljvv8vqHhQGIMjk1ua+1lW0sv63Z38ddX9u7fHxKgqiQnWJO1gKWzCjhydgFzirPVgpwmJvtnTNwvoRARCUdmWirV5YH9HQ80ODzC9tbe/eXn1u/tYt3uTv6+di8jY/aIXDQjj8UzC1g8I49FM/KprsijIl9drPFmwq5RM+ti/IRnQLa7x21fgLpGRSQW+gaG2bC3i1d3d7J2d+f+JDl2pmt+ZhoLy3NZWJ63f3/I0b0iczLi9tdqQjjsrlF318a7IiKHITsjlaPnBoqPj3J3Grv2Ud/Yzabgzh6bmrp5ZnMLdz+38zXvn1mQxfyyHOaX5lIVXDNZVZLD7KJsSnMz1JKMkrj/88PMVgBfAV4BbnP3h2MZj4jI4TAzZhQElmmcVlP2mud6B4YCmyc397K1JbB58pbmbv62di/N3QOveW1GWgqzC7OoKMiiPC+T8vxMyvIyKMnNpCQ3g9K8DIpzMijLy6AgK13bYx2GqCZCM/sJgX0MG9192Zjz5wPfAVKBH7v7DQe5jAPdBDYFbohiuCIiUyonI+11u3uM6t43xPaWXna09bK7vY/dHf3sbO+jsWsfa3d38ujGfXT1v34JCAQKnY8mxZLcwKM0N4OinAyKctIpzsmgMDudgux0CrPTKMzOoDgnPWk3Xj6syjKHfXGzMwkksV+MJkIzSwU2AOcRSGwrgcsJJMWvH3CJDwDN7j5iZjOAG939ikN9rsYIRSQZ9A8O09YbWB/Z2hN4tPQM0NYzQEvPvv3nW4LPjdZ1nUhhdjqluRlUFGQyqzCbWYWBQgOjaypnF2VP67qv4exHOGnu/qiZzT/g9IlAvbtvDgZ2G3Cxu3+dQOtxIm1A5kRPmtm1wLUAVVVxu8RRRCRistJTgwkrO6TXD484HX2DtPUGkmJn3yCd/UO0945JpN0D7O3s59ktrezp7Gd4zDKRFINZhdnMLcneX+u1qjRnf6IsmabjmLEYI6wEdow5bgBOmujFZnYJ8CagCPj+RK9z91uAWyDQIoxEoCIiiSQ1xfZ3lYZieMRp7OrfX9Jue0sPO4Il7R5e30Rj177XvD4vM23/DNjRmbE15XksLM+N6/J1cT9Zxt3vAu6KdRwiIskmNcX2tzhPXFDyuuf7BoZpaAsUGtjeGpjws6W5hzXb27j3xV37K/OYQWVR9v4kObqTSGVxNpVF2eRnpU/xT/ZasUiEOwnsdj9qTvCciIhMI9kZqdTOyKd2xutX2/UPDrO1pYdNjT37l4xsbenh7ud2vm6ST2F2OrMKs5hVmMXM4NjkzILAziEzC7OinixjkQhXArVmtoBAAnwX8O4YxCEiIlGSlZ7KkpkFLJn52uLo7k5rzwDbW3vZ2R6o77qzLTArdndHYOPmlp7XLh3557Or+dSblkQt1mgvn7gVWAGUmVkD8AV3/z8z+xDwFwIzRX/i7q9E6PMuAi6qqamJxOVERCTCzIzSvExK8zI5tqp43NfsGxoO7BrS2c/ezn0smKAebMRiiubyiVjR8gkRETnQRMsnpu+CEBERkQhQIhQRkaSWUInQzC4ys1s6OjpiHYqIiEwTCZUI3f1ed7+2sPD1dftERETGk1CJUERE5HApEYqISFJTIhQRkaSWkOsIzawJ2BaBS5UBzRG4TiLSvZmY7s3EdG8mpnszsUjdm3nuXn7gyYRMhJFiZqvGW3wpujcHo3szMd2bieneTCza90ZdoyIiktSUCEVEJKkpER7cLbEOII7p3kxM92ZiujcT072ZWFTvjcYIRUQkqalFKCIiSU2JUEREkpoSoYiIJDUlQhERSWpKhCIiktSUCEVEJKkpEYqISFJTIhQRkaSmRCgiIklNiVBERJKaEqGIiCQ1JUIREUlqSoQiIpLUlAhFRCSpKRGKiEhSUyIUEZGkpkQoIiJJTYlQRESSmhKhiIgkNSVCERFJakqEIiKS1JQIRUQkqSkRiohIUlMiFBGRpKZEKCIiSU2JUEREkpoSoYiIJDUlQhERSWpKhCIiktSUCEVEJKkpEYqISFJTIhQRkaSmRCgiIkktLdYBRENZWZnPnz8/1mGIiEgcWb16dbO7lx94PiET4fz581m1alWswxARkThiZtvGO6+uURERSWpKhCIiktSUCEVEJKkpEYqISFJTIhQRkaSmRCgiIkktIZdPiIjI9DIy4rT2DrCno5/dHf3s6eynsbOfPR39nHtEBecvmxW1z1YiFBGRqHN32noH2dLcw9bmHra39rKzvY+GtsDXvR37GBgeec17UgzK8zNZMqsgqrHFfSI0sxXAV4BXgNvc/eFYxiMiIhMbGBphR1svm5t62NTUzabGbuqbutnc1ENH3+D+15nBjPwsKouzOXZuMbOXZzOrMIuZhVnMLAh8Lc3NIC01+iN4MUmEZvYT4EKg0d2XjTl/PvAdIBX4sbvfADjQDWQBDTEIV0RExujqH2RHax872nrZ0drLtpZetrf2sq2lhx1tfQyP+P7XludnUlOex4VHzWJBWS4LynKZX5bLnOJsMtNSY/hT/EOsWoQ/A74P/GL0hJmlAj8AziOQ8Faa2T3AY+7+iJnNAG4Erpj6cEVEkkfPviF2d/Sxq72fXe19bG8NJLodwa9tvYOveX1BVhrzSnM5cnYhFx09e3/CW1ieR2F2eox+itDFJBG6+6NmNv+A0ycC9e6+GcDMbgMudvdXg8+3AZkTXdPMrgWuBaiqqop4zCIi04270zMwTEv3Plp7BmjvG6S9d4C2nkE6+gbp7B+ks2+Ijr4BWnoGaO0ZoLV7gK59Q6+5TlqKUVmcTVVJDucvm8W80hzmFucwtySbeSW5FObEf7I7mHgaI6wEdow5bgBOMrNLgDcBRQRakeNy91uAWwDq6up8oteJiExnfQPDNHfvo7FrH83BBNfaM0BL9wCtPfto6RmgrTeQ0Jp7BhgYGpnwWvmZaRRkp1OQnU5pbgZzi3Moyc1gRkEWs4uymFUYGLebVZg1JWN1sRJPiXBc7n4XcFes4xARmQqd/YNsbwmMu21r7aGhrY9d7X3sDnZTHthaG5WbkUpxbgaluRmU52WyeEYBZXkZlOQGHqV5GRTlZFCUnU5xTgb5WWkJndwORzwlwp3A3DHHc4LnREQSyr6hYba39LKpqYctzT1sae4Ofu2huXvgNa8tyc1gdlEWVaU5nLywhBmFWZTlZVKen0l5Xub+RJeVHh8TT6ajeEqEK4FaM1tAIAG+C3j34VzAzC4CLqqpqYlCeCIih6d3YIgNe7tZv6eT+sZuNgWXFOxo7WXMxErK8jJZWJ7LG46YwfyyXOaX5lBVkktVaQ55mfH0azoxmfvUD6eZ2a3ACqAM2At8wd3/z8wuAL5NYPnET9z9q5O5fl1dnWtjXhGZKu7OzvY+1u7u4tVdnby6u4O1u7vY0dbL6K/YjLQUFpblUl2RR3V5HtXl/1hKUJA1vSebTBdmttrd6w48H6tZo5dPcP4+4L4pDkdEJCTDI87Otj42NwcWiG9sDLT2Nuztpjs4dmcGC0pzWVZZwKXHz2HRjHyWzMxnbkkOqSkW459AxnPYidDMSkJ42Yi7tx9+OCIisdU/OExDW3CyypiF4ttae2lo7XtNGbDC7HQWz8znkuMqWTQjn6WzC1g8I59cdWdOK5P5r7Ur+DjYnzapgBbziUhc6R0YorlrgKbufTR29rO3s5+9XfvY1d5HQ1ug7uXezn2veU9eZhpVJTksnpHPeUtnUF2Wx4Jgt2ZpbgZmauVNd5NJhGvd/diDvcDMnptkPGHRZBmR5DQ84jR17WNney87g8sM9nT009jVT2PnPvZ29dPcNUDf4PDr3puWYswqymJOUQ5n1pYzpziHeaU5VJXmMK8ksK5OyS6xHfZkGTPLcvf+cF8TTZosI5J42nsH2NoSKPM1WuNyR+s/di8YHH7t77L8zDQqCjKpyM+ioiCw1KA0L5OyvAzK8jOZkZ/FjIJMinMySNHYXVKI2GSZ0QRnZtVAg7vvC+4QcRTwC3dvj2USFJHpa2h4hG2tvWzc2019Yxf1jd1saella/Nrdy6AwPq6uSU5HFlZyJuWzWROcQ5zirKpLM5mdlG2lh1IyML5l3InUGdmNQRKm/0B+A1wQSQCE5HE1j84zNrdnby8s4NXdnXy6u5O1u/pYt+YkmCVRdksLM/loqNnMb80l6qSQJfl3OIcTUiRiAnnX9KIuw+Z2duB77n792I1Nigi8W1gaIT1e7p4oaGdlxo6eHFnBxv2du3frqcoJ50jZxfw3lPmsXhmAYtmBNbaKdnJVAjnX9mgmV0OvA+4KHhOq0JFktzQ8Ajr93bxUkMHL+0MPNbt7tq/7KAoJ53llYWcs2QhyyuLWD6nkNmFWZqQIjETTiJ8P/BB4KvuviVYGu2XkQlrcjRrVGTq7enoZ/W2Np7b3hZo8e3soH8wkPTys9JYXlnI+0+bz1FzijhqTiFzirOV9CSuhJMIz3P3j4weBJNhTCfJuPu9wL11dXXXxDIOkUTVOzDEyzs7ebGhned3tPPc9nZ2tvcBgRJiy2YX8O4T53H03EKOnlNEVUmOZmRK3AsnEb4P+M4B564a55yITDMjI4HamRsbu1i7u4v1e7pYFywcPVosenZhFsfOK+bq0xdw/LxijphVQEaatvWR6WcyJdYuJ7ArxAIzu2fMU/lAa6QCE5HoGxgaYVtLoGZmfWP3/q+bm7pfN3vziFn5nH/kTI6eW8RRc4ooz8+MYeQikTOZFuGTwG4CO0f895jzXcCLkQhKRCKrrWeAzc09bG3uYXNzINnVN3azraWXoTH7Ac0pzqamIo/TqkupqcijuiKPxTPztTuCJLTJLKjfBmwDTol8OCIyWf2Dw2wO7ne3qambrc09bGkJFIxu7/3HYvS0FGNeaQ41FXmcv2wmNRV51Fbks7A8l5wMLVeQ5DOZrtHH3f10M+sCxtY0MsDdvSBi0R0mzRqVZODu7Oro5+WdHazd3cna3Z2s29PF9tZ/7H1nBrMLs5lflsNblgcWoy8MFoqeW5JDeqrG8kRGxWRj3mhTrVFJJE1d+3hhRzsvNLTzYkMHL+/soKVnAPjH3ndLZuWzaEZ+oDuzPI8FZblkpafGOHKR+BKVjXnNLBWYMfY67r49nGuKJLOBoRFe3d3Jmm1trNne9prlCakpRm1FHucsqeCoOYUcWVnIkpn56s4UCdOk/w8ysw8DXwD2AqPTy5xA8W0RCUHfwDBrtrfxzJZWnt3SwnPb2/fP1pxdmMWxVcVcdep8jqkqYtnsQrIz1MoTibRw/pT8KLDY3VsiFYxIohsZcV7e1cFjG5t5fGMzq7e1MTA8QorBkbMLueKkedTNL+a4qmJmFmbFOlyRpBBOItwBdEQqEJFE1dE7yKMbm3hofSOPrG/aP753xKwCrjptPqdUl1I3r5h8LVEQiYlwEuFm4GEz+xOwb/Sku98YdlQi05i7s25PFw+tb+ThdU2s3t7G8IhTlJPOWYvKWbG4nNNryrUgXSROhJMItwcfGcFHzGn5hMRKR98gT9Q388j6Jh7Z0MSezkDZ3WWVBVx3VjVnL6ngmLlFpKrupkjc0fIJkUlwd17Z1ckjG5p4eH0ja7a3MzziFGSlcUZtOWcuKmPF4gpmFGicTyReRHz5hJk9xGsX1APg7udM9poi8ax73xCPb2zmoXWNPLi+kaauwIjAaKtvxeJyjplbRJoWq4tMK+F0jX5yzPdZwDuAofDCEYkf7s76vV08sr6JRzc2sXJLYIZnfmYaZy4u5+zFFZy5qIyKfLX6RKazSSdCd199wKknzOzZMOMRiRl3Z1tLL09tbuHpzS08tamFxmCrb8nMfK46bT5nL66gbn6xSpSJJJBwukZLxhymAMcDhWFHJDJFhkectbs7Wbm1lVVb21i5tXV/4ivPz+TkhaWcUVPGGYvKmFWYHeNoRSRawukaXU1gjNAIdIluAa6ORFAi0TA84ryyq4OnN7fwzOZWnt3aSld/oDe/siibU6pLOWF+CadUl7KwLBczzfAUSQbhdI0uiGQgIpE22tX5eH0zT9Q38+SmFjr6AtsRLSjL5S3LZ3HSwhJOXFBKZZFafCLJStV6JaG09gwEEt/GZh6vb95fsHp2YRZvXDqD02rKOKW6VMsaRGS/hEqEWlCffAaHR1i9rY1HNjTx2MYmXtnViTsUZKVxanUZH1xRzWnVpSxQV6eITEAL6mXa6egb5MF1e/n72kYe3dBEZ/8QaSnGcVXFnF5bxhm1ZSyvLNR6PhF5jWgsqDfgCmChu3/ZzKqAme6uJRQSce29A/z1lb3c9/JunqhvZnDYKcvL5E1HzuScJRWcXlumotUiMinhdI3+kMA+hOcAXwa6gDuBEyIQlwj7hoZ5cG0jdz+3k4fWNzI47Mwpzub9py3g/GUzOWZOESmq3SkiYQonEZ7k7seZ2XMA7t5mZnFRfFumL3fn5Z2d3L5qB394fied/UOU5WXy3lPmc/Exs1leWaixPhGJqHAS4aCZpRKsN2pm5fxjp3qRw9LRN8jdaxr47aoG1u7uJDMthfOXzeQdx83h1OpSjfeJSNSEkwi/C9wNVJjZV4FLgc9HJCpJCu7O8zva+c0z27n3xV30D46wvLKQr7xtGW89ejaF2RrzE5HoC2dB/a/NbDVwLoHqMm9z97URi0wSVv/gMPe+sIufP7WVl3d2kpORytuPncMVJ1WxrFJV+kRkah12IjygxmgjcOvY59y9NRKBSeJp7OrnF09u4zfPbqe1Z4Daijy+cvGRvO3YSs34FJGYmUyLcGyN0SqgLfh9EYEd61V6TV5jc1M3P3psC3euaWBweIQ3HDGDq06dz6nVpZr4IiIxd9iJcLTGqJn9CLjb3e8LHr8ZeFtEo5Npbe3uTr7/YD33vbyb9NQU3nn8HP7pjIUsKMuNdWgiIvuFM1nmZHe/ZvTA3e83s29GIKZJU4m1+PDKrg6++/eN/OWVveRlpnHdWdW8/7QFlOdnxjo0EZHXCScR7jKzzwO/Ch5fAewKP6TJc/d7gXvr6uquOeSLJeJe3dXJt/+2gb++upf8rDQ+em4tHzhtAYU5Gv8TkfgVTiK8HPgCgSUUAI8Gz0mSWbu7k+/8bSN/fmUP+VlpfOwNtbz/tAVa/iAi00I4yydagY+aWX7g0LsjF5ZMB6/u6uS7fw8mwMxgC/B0JUARmV7CKbq9HPgFUBI8bgbe5+4vRyg2iVOvaQFmpvGRc2u5Wl2gIjJNhdM1ejPwCXd/CMDMVgC3AKeGH5bEow17u/j23zZw30tjWoBKgCIyzYWTCHNHkyCAuz9sZpoXn4C2NPfwPw9s4N4Xd5GbkcZHzqnh6tMXKgGKSEIIJxFuNrN/B34ZPL4S2Bx+SBIv9nb2852/b+S3K3eQkZrCB8+q5tozFlKcq01GRCRxhJMIPwB8CbiLQKWZx4LnZJrr6Bvkfx/exE+f2MKIO1eeVMU/n1NDRX5WrEMTEYm4cGaNtgEfiWAsEmP9g8P88qltfP+hejr7B3nbMZV84rxFzC3JiXVoIiJRE86s0QeAd7p7e/C4GLjN3d8UodhkigyPOL9/bic3PrCBne19nLWonM+cv4SlswtiHZqISNSF0zVaNpoEYf8O9RXhhyRTxd15eH0T3/jzOtbt6WJ5ZSH/delRnFpTFuvQRESmTDiJcMTMqtx9O4CZzSO4W73Ev5d3dvDVP63lqc0tzCvN4fvvPpYLls0iJUW7QYhIcgknEX4OeNzMHiGwDdMZwLURiUqiZld7H9/6y3ruem4nJbkZfPniI3nXCVVkpKXEOjQRkZgIZ7LMn83sOODk4KmPuXtzZMKSSOsfHObmRzbzw4frceC6FdVct6KaAm2IKyJJbjI71M909z0AwcT3x4O9RmLvofWNfPGeV9jW0stbjprFv755CXOKNRNURAQm1yK8DzguAq+RKGvs6uc/fv8Kf35lDwvLc/nV1Sdxeq0mwoiIjDWZRHi0mXUe5HkDDvZ81Ghj3gB3554XdvGFe16hd2CYT71pMdecsVDjgCIi4zjsROjuqdEIJBK0MS80d+/jc3e/xF9e2csxc4v41juPpqYiL9ZhiYjErXBmjUqceWxjEx//7Qt09g/y2Tcv4ZozFpKq5RAiIgelRJgABodH+NZf13PzI5uprcjjV/90IktmqiqMiEgolAinuV3tfVz/6zU8v6Odd59Uxb+/ZSnZGXHbey0iEnfCSoRmdjpQ6+4/NbNyIM/dt0QmNDmUpza18KHfrGHf0Ag/vOI4Llg+K9YhiYhMO+EU3f4CUAcsBn4KpAO/Ak6LTGgyEXfnJ09s5Wv3rWVeaQ63vKdOE2JERCYpnBbh24FjgTUA7r7LzPIjEpVMaGBohH+7+yV+t7qBNy6dwX9fdjT5qg4jIjJp4STCAXd3M3MAM8uNUEwyga7+Qa7/9Roe29jMR86t5WPn1qpItohImMJJhLeb2c1AkZldQ2B3+h9FJiw50N7Ofq766Uo27O3im5cexWV1c2MdkohIQgin6Pa3zOw8AlVkFgP/4e4PRCwy2W9zUzdX/vgZOvoG+clVJ3DWovJYhyQikjDCmjUaTHxKflG0uamby3/0NEPDzm//3yksqyyMdUgiIgllMrtPdDH+BrwGuLtrJXeEjE2Cv7nmZBbP1FwkEZFIm0ytUf02ngJKgiIiUyPsyjJmVgFkjR67+/Zwr5ns9nT0c8WPn2Fo2Ln12pNZNENJUEQkWia9L4+ZvdXMNgJbgEeArcD9EYoraXX1D/L+n62kq3+IX159kpKgiEiUhbNB3VeAk4EN7r4AOBd4OiJRJanB4RGu//UaNu7t4odXHMfS2RpuFRGJtnAS4aC7twApZpbi7g8RKLkmk+Du/OtdL/HYxma+dslyztQSCRGRKRHOGGG7meUBjwK/NrNGoCcyYSWfmx/dzO9WN/DRc2u1WF5EZAqF0yK8GOgFPg78GdgEXBSJoJLN05tb+K+/rOcty2fxsTfUxjocEZGkEk6LsALY7e79wM/NLBuYAbREJLIk0djZz4dvfY55JTnc8I7lmKl2qIjIVAqnRXgHMDLmeDh4TkI0NDzCh299jq7+Qf73yuO1i4SISAyEkwjT3H1g9CD4fUb4ISWP/35gA89saeVrb1+uBfMiIjESTiJsMrO3jh6Y2cVAc/ghJYenN7fwvw9v4vIT53LJcXNiHY6ISNIKJxF+EPg3M9tuZjuAzwD/LzJhvZaZ5ZrZKjO7MBrXn2q9A0N8+ncvMq80h3+/cGmswxERSWqTToTuvsndTwaWAke4+6nuXh/Ke83sJ2bWaGYvH3D+fDNbb2b1ZvbZMU99Brh9srHGm2/cv47trb188x1HkZMRdpU7EREJQzgl1j5qZgUE1g5+28zWmNkbQ3z7z4DzD7heKvAD4M0EkuvlZrY0uOfhq0DjZGONJ09tauHnT23jqlPnc9LC0liHIyKS9MLpGv2Au3cCbwRKgfcAN4TyRnd/FGg94PSJQL27bw5OvLmNwFrFFQRKub0buMbMxo3ZzK4Ndp+uampqmszPE3U9+4b49J0vMK80h0+fvzjW4YiICOGtIxxd8HYB8At3f8XCWwRXCewYc9wAnOTuHwIws6uAZncfGee9uPstwC0AdXV14+2XGHM3PrCBhrY+fnvtKeoSFRGJE+H8Nl5tZn8FFgD/amb5vHZdYUS5+8+ide2pUN/Yzc+f3Mq7TpjLiQtKYh2OiIgEhZMIrwaOATa7e6+ZlQLvD+N6O4GxRTbnBM8lhK/dt5as9FT+5Y3qEhURiSfhjBE6gUktHwke5zJmg95JWAnUmtkCM8sA3gXcE8b14sYjG5p4cF0jHz6nhrK8zFiHIyIiY4STCH8InAJcHjzuIjDr85DM7FbgKWCxmTWY2dXuPgR8CPgLsBa43d1fOZyAzOwiM7ulo6PjcN4WVUPDI3zlj68yrzSHq06bH+twRETkAOF0jZ7k7seZ2XMA7t4WbMkdkrtfPsH5+4D7JhuQu98L3FtXV3fNZK8Rab95djv1jd3c/J7jyUxLjXU4IiJygLA25g2u/XMAMysnipNlpqOOvkFufGADp1aX8salM2IdjoiIjCOcRPhd4G6gwsy+CjwOfC0iUSWIXz29jfbeQf7tgiO0vZKISJyadNeou//azFYD5xJYU/g2d18bscgmwcwuAi6qqamJZRgA9A8O89MntnDWonKWVRbGOhwREZlAOC1CgL3AY8CTQLaZHRd+SJPn7ve6+7WFhbFPPHes2kFz9wDXr6iOdSgiInIQk24RmtlXgKuATQTHCYNfzwk/rOltaHiEmx/dzHFVRVo8LyIS58KZNXoZUD12c14J+NNLu2lo6+MLFx2psUERkTgXTtfoy0BRhOJIGO7O/z68idqKPM5dUhHrcERE5BDCaRF+HXguuKfgvtGT7v7Wid+S+B5a38i6PV3ceNnRpKSoNSgiEu/CSYQ/B74BvEScrB+Mh1mjP35sC5VF2Vx09OyYxSAiIqELJxH2uvt3IxZJBMS6skxDWy9PbmrhX85bRHpquBNyRURkKoSTCB8zs68TKIw9tmt0TdhRTVN/eH4XAG87tjLGkYiISKjCSYTHBr+ePOZc0i6fcHfuXNPASQtKmFuSE+twREQkROFUljk7koFMdy80dLC5qYcPnqkF9CIi08lhJ0Izu9Ldf2VmnxjveXe/Mfywpp+71jSQmZbCm5fPjHUoIiJyGCbTIswNfs0f5zkf59yUidWs0YGhEe55YRdvOnIm+VnpU/rZIiISnsNOhO5+c/Dbv7n7E2OfM7PTIhLVJMVq1uhD6xtp7x3kkuM0SUZEZLoJZ47/90I8l/DuWtNAeX4mp9eUxToUERE5TJMZIzwFOBUoP2CcsABIui3Y23oGeHBdI1edOp80rR0UEZl2JjNGmAHkBd87dpywE7g0EkFNJw+8upfBYefiY9QtKiIyHU1mjPAR4BEz+5m7b4tCTNPKY/XNlOdncuTsgliHIiIikzDpvjwlQRgZcZ6sb+b0mjJttyQiMk0l1KCWmV1kZrd0dHRMyeet29NFS88Ap2mSjIjItDXpRGhm5ZEMJBLc/V53v7awsHBKPu/x+iYAzRYVEZnGwmkRPmFmfzWzq82sOGIRTSOP17dQU5HHzMKsWIciIiKTFM4Y4SLg88CRwGoz+6OZXRmxyOLcvqFhnt3SotagiMg0F9YYobs/6+6fAE4EWgls1psU1mxrp39wROODIiLTXDhjhAVm9j4zux94EthNICEmhcfrm0hNMU5aWBLrUEREJAzh7Ef4AvB74Mvu/lRkwpk+Hq9v4Zi5RRSoyLaIyLQWTiJc6O5uZnlmlufu3RGLKs519A7yUkM7HzqnNtahiIhImMIZIzzSzJ4DXgFeNbPVZrYsQnHFtac2tzDiWjYhIpIIwkmEtwCfcPd57l4F/EvwXMxM1YL6x+ubyM1I5diqoqh+joiIRF84iTDX3R8aPXD3h/nHpr0xMVUL6p+sb+HEBSWka7cJEZFpL5zf5JvN7N/NbH7w8Xlgc6QCi1e9A0NsaenhmLlJWUNARCThhJMIPwCUA3cFH+XBcwltc1MP7lA7Iy/WoYiISARMetaou7cBH4lgLNPCxsYuAGorlAhFRBLBpBOhmS0CPgnMH3sddz8n/LDiV31jN6kpxrzSmA6HiohIhISzjvAO4Cbgx8BwZMKJfxv3djO/NIeMNE2UERFJBOEkwiF3/9+IRTJN1Dd1s6giP9ZhiIhIhITTrLnXzK43s1lmVjL6iFhkcWjf0DDbWnqp0figiEjCCKdF+L7g10+NOefAwjCuGde2NvcyPOKaMSoikkDCmTW6IJKBTAf1jYFyqtXlSoQiIokinFmj7x3vvLv/YvLhxLeNjV2YKRGKiCSScLpGTxjzfRZwLrAGiFkiNLOLgItqamqicv36xm7mFueQnZEaleuLiMjUC6dr9MNjj82sCLgt3IDC4e73AvfW1dVdE43r1zd2a6KMiEiCieRiuB4gYccNh4ZH2Nzco4oyIiIJJpwxwnsJzBIFSAWOAG6PRFDxaEdbHwNDI2oRiogkmHDGCL815vshYJu7N4QZT9zauDdQY1SJUEQksUy6a9TdHwHWA4VACYFkmLDqmwJLJ5QIRUQSy6QToZn9E/AscAlwKfC0mSXsNkz1e7uZVZhFflZ6rEMREZEICqdr9FPAse7eAmBmpcCTwE8iEVi8qW/SjFERkUQUzqzRFqBrzHFX8FzCGRlxLZ0QEUlQ4bQI64FnzOwPBGaPXgy8aGafAHD3GyMQX1zY1dFH78CwEqGISAIKJxFuCj5G/SH4NeH2KBqtMVqr7ZdERBJOOJVlvhTJQOLZPxKhWoQiIokmnAX15cCngSMJ1BoFwN3PiUBccWVTUzcluRkU52bEOhQREYmwcCbL/BpYR6Cs2peArcDKCMQUdzY19lCjHSdERBJSOImw1N3/Dxh090fc/QNAwrUGIdAirK7IjXUYIiISBeFMlhkMft1tZm8BdhGoMJNQ2noGaOkZ0B6EIiIJKpxE+J9mVgj8C/A9oAD4eESimqRo7Ee4uVm70ouIJLJwukafcvcOd3/Z3c929+Pd/Z6IRTYJ7n6vu19bWFgYsWuOzhhVIhQRSUzhJMKnzewOM7vAzCxiEcWZTU09ZKSlUFmcHetQREQkCsJJhIuAW4D3ABvN7GtmtigyYcWPTY3dLCzLJTUlYXO9iEhSC2cbJnf3B9z9cuAa4H3ASjN7xMxOiViEMbapqVvdoiIiCSycbZhKzeyjZrYK+CTwYaCUwOSZ30QovpjaNzTM9tZeqsu1dEJEJFGFNVmGwEzRt7n7W9z9LuBD7r4KuCki0cXYtpZeRhyqVVpNRCRhhbN8YrG7+wHnPgF8292/EcZ148YmzRgVEUl4YY0RjnM6oWaUbGoKJMKF6hoVEUlY4XSNjme85DhtbWrqobIom5yMcBrOIiISzw77N7yZdTF+wjMgoRbbbWrqVmtQRCTBHXYidPek2J3W3dnU2M076+bGOhQREYmiSHeNJoy9nfvoGRjWjFERkQSnRDiB0YkyWkMoIpLYlAgnMFpsWxvyiogkNiXCCWxq6iY/M43y/MxYhyIiIlGkRDiBTU3dLKzII4E31hAREZQIJ7SpsUfjgyIiSUCJcBzd+4bY09mv0moiIklAiXAcg0MjXHPGAk5eWBrrUEREJMpUO2wcxbkZfO4tS2MdhoiITAG1CEVEJKkpEYqISFJTIhQRkaQW94nQzI4ws5vM7Hdmdl2s4xERkcQSk0RoZj8xs0Yze/mA8+eb2XozqzezzwK4+1p3/yBwGXBaLOIVEZHEFasW4c+A88eeMLNU4AfAm4GlwOVmtjT43FuBPwH3TW2YIiKS6GKSCN39UaD1gNMnAvXuvtndB4DbgIuDr7/H3d8MXDHRNc3sWjNbZWarmpqaohW6iIgkmHhaR1gJ7Bhz3ACcZGYrgEuATA7SInT3W4BbAMysycy2RSCmMqA5AtdJRLo3E9O9mZjuzcR0byYWqXszb7yT8ZQIx+XuDwMPH+Z7yiPx2Wa2yt3rInGtRKN7MzHdm4np3kxM92Zi0b438TRrdCcwd8zxnOA5ERGRqImnRLgSqDWzBWaWAbwLuCfGMYmISIKL1fKJW4GngMVm1mBmV7v7EPAh4C/AWuB2d38lFvGNcUuMPz+e6d5MTPdmYro3E9O9mVhU7425ezSvLyIiEtfiqWtURERkyikRiohIUkv6RDheWbcDns80s98Gn3/GzObHIMyYCOHefMLMXjWzF83s72Y27hqdRHSoezPmde8wMzezpJkWH8q9MbPLgv92XjGz30x1jLESwv9TVWb2kJk9F/z/6oJYxBkLE5XeHPO8mdl3g/fuRTM7LmIf7u5J+wBSgU3AQiADeAFYesBrrgduCn7/LuC3sY47ju7N2UBO8PvrdG9e97p84FHgaaAu1nHHy70BaoHngOLgcUWs446je3MLcF3w+6XA1ljHPYX350zgOODlCZ6/ALgfMOBk4JlIfXaytwgnLOs2xsXAz4Pf/w4418xsCmOMlUPeG3d/yN17g4dPE1j7mQxC+XcD8BXgG0D/VAYXY6Hcm2uAH7h7G4C7N05xjLESyr1xoCD4fSGwawrjiykfv/TmWBcDv/CAp4EiM5sVic9O9kQ4Xlm3yole44ElHh1A6ZREF1uh3Juxribw11oyOOS9CXbbzHX3P01lYHEglH83i4BFZvaEmT1tZueTHEK5N18ErjSzBgIlJT88NaFNC4f7OylkcV9iTeKfmV0J1AFnxTqWeGBmKcCNwFUxDiVepRHoHl1BoBfhUTNb7u7tsQwqTlwO/Mzd/9vMTgF+aWbL3H0k1oElsmRvEYZS1m3/a8wsjUB3RcuURBdbIZW8M7M3AJ8D3uru+6Yotlg71L3JB5YBD5vZVgLjGfckyYSZUP7dNAD3uPugu28BNhBIjIkulHtzNXA7gLs/BWQRKDgtUSzDmeyJMJSybvcA7wt+fynwoAdHbhPcIe+NmR0L3EwgCSbLOA8c4t64e4e7l7n7fHefT2D89K3uvio24U6pUP6f+j2B1iBmVkagq3TzFMYYK6Hcm+3AuQBmdgSBRKh95QLuAd4bnD16MtDh7rsjceGk7hp19yEzGy3rlgr8xN1fMbMvA6vc/R7g/wh0T9QTGMh9V+winjoh3pv/AvKAO4Lzh7a7+1tjFvQUCfHeJKUQ781fgDea2avAMPApd0/4XpYQ782/AD8ys48TmDhzVZL84T1aenMFUBYcI/0CkA7g7jcRGDO9AKgHeoH3R+yzk+Qei4iIjCvZu0ZFRCTJKRGKiEhSUyIUEZGkpkQoIiJJTYlQRETi2qEKco/z+sMq6q5ZoyIiEtfM7Eygm0Ct0WWHeG0tgaIE57h7m5lVHGqds1qEIpNgZkVmdv2Y49lm9rsofM4XzWxncK1ZuNe6ysy+H4m4DvIZHzSz9475vNmHeP3HzWx7tOOS6W28gtxmVm1mfzaz1Wb2mJktCT512EXdlQhFJqeIwBZdALj7Lne/NEqf9T/u/h9RunZEuftN7v6L4OFVwEETobv/DzAtfjaJO7cAH3b344FPAj8Mnj/sou5KhCKTcwNQbWbPm9l/mdn80fGLYEvo92b2gJltNbMPWWAT4+eC/2OWBF830V+0EzKzPDP7qZm9FNyc9B3B85cHz71sZt8Y8/r3m9kGM3sWOG3M+XIzu9PMVgYfp43zWa9pQZrZH81sRfD7bjP7qpm9EPyZZgTPf9HMPmlmlxIoxP7r4D3KNrMb7B8bOX9rMjddBAL/HwCnEqhq9TyBUo+jWzKNLep+OYFKPUUHu54SocjkfBbY5O7HuPunxnl+GXAJcALwVaDX3Y8FngLeG3zNRH/RHsy/E6ixuNzdjwIeDHY/fgM4BzgGOMHM3maBvdq+RCABnk5go9dR3yHQ0jwBeAfw49B/dABygafd/WgCmw9fM/ZJd/8dsAq4wt2PAXKAtwNHBuP+z8P8PJGxUoD24P9/o48jgs8ddlH3pK41KhJFD7l7F9BlZh3AvcHzLwFHHfAX7eh7MkO47hsYU+82OBngTOBhd28CMLNfE9jtmwPO/5ZAt9HodZaO+ewCM8tz9+4Qf74B4I/B71cD5x3i9R0ENij+PzP745j3ihw2d+80sy1m9k53v8MC/5CPcvcXCBR1vxz4qYVY1F2JUCQ6xm5JNTLmeITA/3f7/6Kd4rhGpQAnu3v/QV4zxGt7jbLGfD84phj0MIf4XRIsOH0igZ0VLgU+RKAFK3JIExTkvgL4XzP7PIHi3LcBLzCJou7qGhWZnC4C+w5Oirt3AlvM7J0AFnB0CG99APjn0QMzKwaeBc4yszIzSyXw1/AjwDPB86Vmlg68c8x1/sqY3c/N7JhxPmsrcIyZpZjZXODEw/gRYcw9CraAC939PuDjQCg/qwgA7n65u89y93R3n+Pu/+fuW9z9fHc/2t2XuvuXg691d/9E8Nxyd7/tUNdXIhSZhOBfmE8EJ6f81yQvcwVwtZm9ALwCXBzCe/4TKA5+7gvA2cE92T4LPETgL+LV7v6H4PkvEhiXfAJYO+Y6HwHqghNXXgU+OM5nPQFsAV4FvgusOcyf72fATcHJDPnAH83sReBx4BOHeS2RqNGCepE4ZmZfBLrdPWFnWZrZVUCdu38o1rFIclKLUCS+dQPXRmJBfTyywAa0/wp0xjoWSV5qEYqISFJTi1BERJKaEqGIiCQ1JUIREUlqSoQiIpLU/n9Ewo8xuJHVOAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 504x576 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fix, ax = plt.subplots(2,1,sharex=True,figsize=(7,8))\n",
"ax[0].set_yscale(\"log\")\n",
"ax[0].set_ylabel(\"Lyapunov exponent [code units]\")\n",
"ax[0].plot(times,lyapunov_exponent)\n",
"ax[1].set_yscale(\"log\")\n",
"ax[1].set_xlabel(\"time [code units]\")\n",
"ax[1].set_ylabel(\"Lyapunov timescale [code units]\")\n",
"ax[1].plot(times,1./lyapunov_exponent);"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8ef02b9",
"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.8.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment