{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6882e778",
   "metadata": {},
   "source": [
    "# Chaotic Hyperion\n",
    "In this example, we simulate the spin of Hyperion. The spin evolution is governed by an ordinary differential equation that is coupled to the moon's orbit. \n",
    "\n",
    "We start by importing REBOUND, numpy and matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e350dbbb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import rebound\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc568b2",
   "metadata": {},
   "source": [
    "The right hand side of ODEs can be implemented in either python or in C. Although not absolutely necessary for this example, we here show how to implement the RHS in C. This is often significantly faster than using a python callback function. \n",
    "\n",
    "We use a simple spin model which is one second order ODE, or a set of two coupled first order ODEs. For more details on the physics behind this model, see  Danby (1962), Goldreich and Peale (1966), and Wisdom and Peale (1983). The RHS of this set of ODEs implemented in C is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "79a48333",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting rhs.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile rhs.c\n",
    "#include \"rebound.h\"\n",
    "void derivatives(struct reb_ode* const ode, double* const yDot, const double* const y, const double t){\n",
    "    struct reb_orbit o = reb_orbit_from_particle(ode->r->G, ode->r->particles[1], ode->r->particles[0]);\n",
    "    \n",
    "    double omega2 = 3.*0.26; \n",
    "    yDot[0] = y[1];\n",
    "    yDot[1] = -omega2/(2.*o.d*o.d*o.d)*sin(2.*(y[0]-o.f));\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33cd1c9c",
   "metadata": {},
   "source": [
    "We now compile this into a shared library. We need the REBOUND headers and library for this. The following is a bit of hack: we just copy the files into the current folder. This works if you've installed REBOUND from the git repository. Otherwise, you'll need to find these files manually (which might depend on your python environment). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7d3c176a",
   "metadata": {},
   "outputs": [],
   "source": [
    "!cp ../src/librebound.so .\n",
    "!cp ../src/rebound.h .\n",
    "!gcc -c -O3 -fPIC rhs.c -o rhs.o\n",
    "!gcc -L. -shared rhs.o -o rhs.so -lrebound "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2421c0e7",
   "metadata": {},
   "source": [
    "Using ctypes, we can load the library into python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11736a91",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctypes import cdll\n",
    "clibrhs = cdll.LoadLibrary(\"rhs.so\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37ac5371",
   "metadata": {},
   "source": [
    "The following function is setting up the N-body simulation as well as the ODE system that governs the spin evolution. Note that we set the `derivatives` function pointer to the C function we've just compiled. You could also set this function pointer to a python function and avoid all the C complications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "18043d1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup():\n",
    "    sim = rebound.Simulation()\n",
    "    sim.add(m=1)             # Saturn\n",
    "    sim.add(a=1, e=0.123233) # Hyperion, massless, semi-major axis of 1\n",
    "    sim.integrator = \"BS\"\n",
    "    sim.ri_bs.eps_rel = 1e-12  # tolerance\n",
    "    sim.ri_bs.eps_abs = 1e-12\n",
    "    \n",
    "    ode_spin = sim.create_ode(length=2, needs_nbody=True)\n",
    "    ode_spin.y[0] = 0.01  # initial conditions that lead to chaos\n",
    "    ode_spin.y[1] = 1\n",
    "    ode_spin.derivatives = clibrhs.derivatives\n",
    "    \n",
    "    return sim, ode_spin"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85c097cf",
   "metadata": {},
   "source": [
    "We will create two simulations that are slightly offset from each other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8182b908",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim, ode_spin = setup()\n",
    "sim2, ode_spin2 = setup()\n",
    "ode_spin2.y[0] += 1e-8 # small perturbation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c597773f",
   "metadata": {},
   "source": [
    "With these two simulations, we can measure the growing divergence of nearby trajectories, a key feature of chaos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "id": "4eb25927",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = 2.*np.pi*np.linspace(0,30,100) # a couple of orbits\n",
    "obliq = np.zeros((len(times)))\n",
    "obliq2 = obliq.copy()\n",
    "\n",
    "for i, t in enumerate(times):\n",
    "    sim.integrate(t, exact_finish_time=1)\n",
    "    sim2.integrate(t, exact_finish_time=1)        \n",
    "    obliq[i] = ode_spin.y[0]\n",
    "    obliq2[i] = ode_spin2.y[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98a92abf",
   "metadata": {},
   "source": [
    "Finally, let us plot the divergence as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "id": "13ac8f3a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEMCAYAAADal/HVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2UElEQVR4nO3deXxU5dXA8d9JiBAXDChuEQStxVqtLBFUqkWqFVd4SYAAUjZBcMWdKgqoCIK+KBb1pewqi6yiolRFBEHZAlao4oJaiRYVCbJESMJ5/5gZmAwzd+5kZpKZzPl+PvmQubnLMwTumfuc5zmPqCrGGGNMKGlV3QBjjDGJzQKFMcYYRxYojDHGOLJAYYwxxpEFCmOMMY4sUBhjjHFkgcIYY4wjCxTGGGMc1ajqBoQjIu2Bq4HawERV/WfVtsgYY1JLlTxRiMgkEflBRDYGbG8rIptF5AsRGQSgqgtUtS/QH+hcFe01xphUVlVdT1OAtv4bRCQdGAdcCZwNdBGRs/12Gez9uTHGmEpUJV1PqrpMRBoGbG4BfKGqWwBEZCbQTkQ+AUYCb6hqQahzikg/oB/AUUcd1fyss86KS9uNMZVr586dfPXVVwA0atSIY4891nH/or0l/PeXXykpO0C6CAiUHQhd0+7cbOfzuVW0t4TComIO+NXPSxMh68gMivaWHLb9gEOdPV+bQp1TwrynwPO4tW7dup9UtV7g9kTKUWQD3/q93gq0BG4FLgOOFZHfqOrzwQ5W1fHAeICcnBxdu3ZtnJtrjImnsrIyHnroIR577DGaNWvGnDlzaNSoUYXO1WrkEgqLig/bnp2VyYpBbaJt6kEL1hcyevFmvisq5tjMDERgx94STgyyb7oIZSGCRc2sTO65ojHtm2aXO+cp3u13zNpAuDBRkfcmIt8E255IgSIoVR0LjK3qdhhjKs+PP/5Ily5deOedd+jbty9jx46lVq1aIfcPdjNt3zT74M/vuaIxf5v3McUlZQe3ZWakc88VjV23KdQ1gm0HDrteoDJVMjPSg+5TWFTM3+Z9DED7ptnl3gvA6MWbgwa+ir63cBIpUBQC9f1en+rdZoxJIatWrSIvL48ff/yRSZMm0atXL8f9F6wvLHdTDrzJ+v/pFEwqco213/zM3HWFh22vlZHmGCTA84n/nisah7zpF5eUMXrx5qBtDBb4BFC/87p9b24kUqBYA5wpIo3wBIh8oGvVNskYU1lUleeee46BAwdy6qmn8sEHH9C0adOwx41evPmwm3Kwm2ywT+ZuhbrGjFXfHtZ9VFxSFjZI+D7x+9rUaNDrQbuSvgvx1BBt4ItUlQQKEZkBtAaOF5GtwBBVnSgitwCLgXRgkqpuqor2GWMq1549e7jxxht56aWXuPrqq3nhhReoU6eOq2ND3UxDba+IUOcKlWNwEuwT/ylZmUGfKk7Jyiz3OlwXW7xU1ainLiG2LwIWVXJzjDFV6LPPPiM3N5dNmzbxyCOPcP/995OW5n7kvtubbDRCXSNUQjorM4N9pQcOy4mM6HBuRF1JhUXFtBq5JGjeI1gXW7xYCQ9jTJWZN28eOTk5fP/99yxevJjBgwdHFCTAc5PNzEgvty3WydxQ1+jSsn7Q7UOv+z0jOpxLdlYmgucpIlSQAM+N3rc/HMo3wKGAMOzVTSG72OItkXIUxpgUUVpayv3338/o0aNp0aIFs2fPpkGDBhU6V0X669104QTuk9s8m3c//fGwY3JOqxvyXJF80vflK4IN5XXKe8Syiy0U0Qr0sSU6m0dhTOLatm0b+fn5LF26lAEDBjBmzBhq1qxZadcPHMEEh7qF4NDQU/9P9f77xLubJ1RiO5RYzgURkXWqmhO43Z4ojDGVZsWKFXTs2JGioiKmTZtG9+7dK70NoUYwDV24qVxeIfBm7TRcNZZC5UNC5T1i2cUWiuUojDFxp6o89dRTtG7dmqOOOooPP/ywSoIEhO6qKSouCTustTK6eULlQyLNe8SSPVEYY+Jq165d3HDDDbz88su0a9eOKVOmkJWVFfV5KzpUNNQndjdiOZIqlHA5l8oIDIEsUBhj4uaTTz4hNzeXzZs3M2LECO69996IRzUF42Y2diihynnUykhjx96SkMdVVjcPRDc5MB6s68kYExcvv/wyLVq0YPv27bz99tsMGjQoJkECnGdjh+M/FNW/C2fItb8/rMtHvH9WZjdPIrInCmNMTJWUlHDvvffy1FNPceGFFzJ79myys2N7g412NrbTJ/aqmPmc6CxQGGNi5rvvvqNTp06sWLGC2267jdGjR3PEEUe4OjaS6qzxmo2daF0+icIChTEmJpYuXUrnzp3Zs2cPM2bMID8/3/WxkVZnzW2eXW47VDyHUFX1k5KJ5SiMMVFRVUaPHs1ll11GnTp1WLVqVURBApyrswbb/uKH/6FmjTTqHJkR1VBRX4AqLCpGORSIFqy3FQ782ROFMabCdu7cSa9evZg/fz55eXlMnDiR2rVrR3yeilRnLSouITMjnTGdm8S8fHhlTKxLJvZEYYypkI0bN3L++eezcOFCnnzySV5++eUKBQkInVtIFwm63SfaoniVUaK8OrBAYYyJ2PTp02nZsiW7du3i3Xff5c4770TC3NSdRFKdNVA0N/VQAaoyJtYlEwsUxhjX9u/fzy233EK3bt1o3rw5BQUFXHzxxa6PX7C+kFYjl9Bo0Ou0GrnkYC4g1NyGR9ufW678djDR3NQro0R5dWDVY40xrnz77bd07NiRVatWcddddzFixAgyMjJcH+9UtdVNPiDa453Oa6OePEJVj7VAYYwJ65133iE/P59ff/2VyZMnk5eXF/E5gq2z4BNsedBg7KYeX1Zm3BgTsQMHDjBy5EgefPBBzjrrLObNm0fjxhXrlnHKJbit1WQT4qqG5SiMMUHt2LGD9u3b88ADD9CpUydWrVpV4SAB4XMJlbWsp4mcBQpjzGE2bNhATk4Ob7zxBmPHjmX69OkcffTRUZ0zWOI4kA1LTUzW9WSMKWfq1Kn079+f4447jmXLlnHhhRfG5Lz+6yyEylXYsNTEZE8UxhgA9u3bR//+/enZsycXXnghBQUFMQsSPu2bZrNiUBue6tzEhqUmEXuiMMbwzTffkJeXx9q1a7nvvvt49NFHqVEjfreHcKu4mcRigcKYFLd48WK6du1KaWkp8+fPp3379pVyXacRTDYMNrFY15MxKerAgQMMGzaMK6+8kuzsbNauXVtpQcKJVXRNPBYojElB27dv55prrmHo0KFcf/31fPjhh5x55plV3SwgumVOTXxY15MxKWbdunXk5uby3Xff8dxzz3HjjTdGVdAvnEi7kayia+KxQGFMCpkwYQK33HILJ5xwAu+//z4tWrSI2bmDBQQg6Mp1EHoGdryWOTUVZ11PxqSA4uJi+vTpQ9++fbnkkksoKCiIeZAIllcY9uqmiLuRrKJr4gn5RCEiY10c/4uqDo5he4wxMbZlyxZyc3PZsGEDDz74IEOGDCE93XmGdKRC5RUCt/k4dSPZ0NnE49T11A54KMzxgwALFMYkqNdee43u3bsf/P7qq6+Oy3UizR+E60ay4n+JxSlQjFHVqU4Hi0idGLfHGBMDZWVlDBkyhOHDh9OkSRPmzp3L6aefHrfrhcorAAjgv5iBdSMln5A5ClV9KtzBbvaJlogcJSJTReQfItIt3tczJtn99NNPXHnllQwfPpzevXuzcuXKuAYJcC74p3iCBRxauc6eFpJL2GS2iIwSkdoikiEi74jIjyJyfTQXFZFJIvKDiGwM2N5WRDaLyBciMsi7uQMwR1X7AtdFc11jqrvVq1fTrFkzli1bxoQJE5g4cSKZmfEfLeS/lGkwiidIrBjUxoJEEnIz6ukvqvoLcA3wNfAb4J4orzsFaOu/QUTSgXHAlcDZQBcRORs4FfjWu1vwzJgxKU5Vee655/jjH/9Ieno6K1asoE+fPpXaBl/Bv1AzMmweRPJyEyh8eYyrgdmqujPai6rqMuDngM0tgC9UdYuq7gdm4kmob8UTLNy215iUsnfvXnr06MFNN93E5Zdfzrp162jevHmVtSdUotrmQSQvNxPuXhORT4FiYICI1AN+jUNbsjn05ACeANESGAv8XUSuBl4NdbCI9AP6ATRo0CAOzTMm8Xz++efk5uayceNGhg0bxuDBg0lLi+zzVDQF+IIde88VjctNsgNLYCe7sP+iVHUQcBGQo6olwF48n/QrharuUdVeqjpAVV9y2G+8quaoak69evUqq3nGVJkFCxaQk5NDYWEhb7zxBg899FCFgkRFC/CFOhY4mK8QLIFdHYR9ohCRI4GbgAZ4PrGfAjQGXotxWwqB+n6vT/VuM8b4KS0t5YEHHmDUqFHk5OQwZ84cTjvttAqdy6kAX7gbu9OxlrSuXtx8/JgM7MfzVAGem/ejcWjLGuBMEWkkIkcA+cDCOFzHmKS1bds2/vKXvzBq1ChuvPFGli9fXuEgAdEV4LPifanDTY7iDFXtLCJdAFR1r0RZalJEZgCtgeNFZCswRFUnisgtwGIgHZikqpuiuY4x1cnKlSvp2LEjP//8M1OmTKFHjx5Rn9NtAT7/XMSxmRmIlJ9E53SsSX5uAsV+EcnE++9CRM4A9kVzUVXtEmL7ImBRNOc2prpRVZ555hnuuusuTjvtND788EPOO++8qM7pu/EXFhWHnTnty0X4upmKiktCnteS1tWTm0AxBHgTqC8iLwGtgJ7xbJQxqcz/0/uJmZDxwT94f/ErXHvttUybNo2srKyoz+9/4/fNnPZNigsc9RQsFxFMsGNN9RA2UKjqWyJSAFyA59/T7ar6U9xbZkwK8r+Jl2z/loL5j1HycyFZl/yVny7oztKv9tC+aVZU1wh24/cPEqMXb+aOWRsODnd1k3MQYMWgNlG1yyQuUQ3V0+jdwZOP6AacrqoPi0gD4CRVXV0ZDayInJwcXbt2bVU3w5iItRq5hMKiYvZ8+j7b33gaqXEEx197D5kNmwDOn/zdajTo9ZD5hcyM9MPmP9TKSGPH3tDdTXCoPIdJbiKyTlVzAre7GfX0LHAh4Msr7MJTasMYE2OF23fx8zv/4KdXRpJxXANO7vH0wSABh3IJhUXF3DFrAw0HvU6rkUtczXvwCZVsThcJOtxVlZAF/8DyEqnATaBoqao3452Nrao7gCPi2ipjUtD333/PjjkPsmvtKxzT7BpO6jaSGrWPD7m/f9BwO0kOQq8gVxaid2FncUm5CXRZmRnUOTLDJtOlEDfJ7BJvwT7fqKd6wIG4tsqYFLNs2TI6derEvp2/cHL7ezmi8SURHV9cUsbAWRsYvXhz2C6pUCvI+UZBBTolK9MWEkpxbgLFWGA+cIKIDAfysFXtjIkJVeXJJ59k0KBBnHHGGbz99tt8UVIn5NDVcPzLaIQLFsF+bjWa3ImmPlYyckxmi0gantFOPwN/xpNLe0dVP6mc5lWMJbNNMvjll1/o3bs3c+fOpUOHDkyePJnatWuX28dpvoOTiiaXg02sK9pbkhI3Q7cChxeDJ6BWhy64UMlsN6Oe1qtq07i1LA4sUJhEt2nTJjp06MCXX37JiBEjuPvuuwlX8CCSoCHAVyMrvj52db4ZRss3Mi1QdRj5FSpQuOl6ekdEcoF5Gi6qGGPCmjFjBjfccAPHHHMM77zzDn/6059cHeffXeQfNIIJNrIpku6SaIoFVnepWOPKzainG4HZwD4R+UVEdonIL3FulzHVzv79+7n99tvp2rUrTZs2paCgwHWQCORbTe6pzk2CjmAKzCtEWk48FW+GbqXiwkxu1qM4RlXTVPUIVa3tfV073HHGmEO2bt1K69atGTt2LHfccQfvvvsup5xyStTn9V+r2jd0tVZGGnfM2lBufoXTE0IwqXgzdCvU8OLqnPR3sx5FsyCbdwLfqGpp7JtkTPWyZMkS8vPzKS4uZtasWXTq1Cmm5/d1SQXmFfxHQEX6hGCr1IUWanhxde6Sc5OjeBZoBnzsfX0usBE4VkQGqOo/49U4Y5KZqjJq1Cjuv/9+GjduzNy5c/nd734Xt+s5PTW4LSfuk4o3w0ik2rwSN4HiO6CPb20IETkbeBi4F5gHWKAwJsDOnTvp0aMHr7zyCp07d2bChAkcffTRIfePxbh8p6eGMZ2bRPyEkGo3QxOam0DxW/8FhFT13yJylqpuiXL9ImOqpX/961/k5uby9ddf89RTT3Hbbbc5Dn116jKK5Ebt9NRgTwgmGm4CxSYReQ6Y6X3dGfi3iNQEnEtKGpNipk2bRv/+/cnKymLp0qW0atUq7DGxGooaLq9gTwimotwMj+0JfAEM9H5t8W4rAS6NT7OMSS779u1jwIAB9OjRgxYtWlBQUOAqSEDshqIGjoCygn0mVtwsXFQsIs8Cr6lq4Fi63fFpljHJ4z//+Q95eXmsWbOGe+65h8cee4waNdw8rHtEmmh2Yk8NJh7CPlGIyHXABjzLoSIiTURkYZzbZUxS+Oc//0mzZs349NNPmTt3LqNGjYooSEBqjss3ycVN19MQoAVQBKCqG4BG8WuSMYnvwIEDPPLII7Rt25aTTz6ZtWvX0qFDhwqdy7qMTKJztR6Fqu4MGLVhNZ9MytqxYwfdu3fn9ddfp1u3bvzf//0fRx11VFTntC4jk8jcjnrqCqSLyJnAbcDK+DbLmMRUUFBAbm4uhYWFjBs3jgEDBoSt+mpMsnPT9XQr8HtgHzAD+AXP6CdjUsqkSZO46KKLKC0tZfny5dx0000WJExKcDPqaS/wgPfLmJTz66+/cuuttzJhwgQuu+wypk+fTr169aq6WcZUmpCBQkRexSEXoarXxaVFxiSQr776iry8PAoKCrj//vt5+OGHSU9PD7pvqi2PaVKH0xPFE94/OwAnAS96X3cBtsWzUcZUtQXrC7n/6WlsnvUYgnL/mMkMH9jTcf9YlOEwJhGFDBSq+h6AiDwZsDTeqyJi64yaamvu2v9w4x2D2P7+DDJOaES99vczb/tJnO+3tkPgU0O49R6CHWNPICZZuFkz+xPgalXd4n3dCFikqvGrlxwlWzPbVNRPP/3EmRddRdHnazjqnMuo+5cBpGXUBDyLAu0rPRB0Hek7Zm0I2U+bmZF+2DG5zbOZu67Q1qQ2CSXUmtluAkVbYDyeGk8CnAb0S+R1KCxQmIpYs2YNeXl5/KfwO+pe1p+jz7vC1aimbG+pjWBlONJFKAvyfyzU9uysTFYMalOB1hsTvVCBws2opze98yfO8m76VFX3xbqBxlQVVWX8+PHcdtttnHTSSfyh/1h2Ht3A9fG+ACGUH/0hEDQY4LDd1qQ2iSjkPAr/JVBVdZ+qfuT92hdsH2OS0d69e+nVqxf9+/enTZs2FBQUMKzPdUFrL9U5MsPxXIonOMDhQcMtW5PaJCKnJ4rJItKaQ//2g5kINI1lg4ypLF988QW5ubl8/PHHDB06lAcffJC0tDTaH+f5eWCiGThsvYdASuhupXCsEKBJVE6B4lhgHc6B4sfYNseYyrFw4UL++te/kp6ezqJFi2jbtm25nzvVXvIFkFChoCJBIttGPZkE5jQ8tmEltsORiLQHrgZqAxMTOZFuEltpaSkPPfQQI0aMoHnz5syZM4eGDRsC7ibM+QeQViOXRJTADkXAEtgmobmp9RQVEZkkIj+IyMaA7W1FZLOIfCEig5zOoaoLVLUv0B/PUqzGROyHH36gbdu2jBgxgn79+vH++++XCxJ/m/cxhd4nBd+EuQXeuRPBhFpHokvL+hHlOCwvYRJdZCusVMwU4O/ANN8GEUkHxgGXA1uBNd7FkNKBEQHH91bVH7zfD/YeZ0xEPvzwQ/Ly8ti+fTuTJ0+mZ8+e5X5ekXWrfduDPYXknFbXVY7D8hImGcQ9UKjqMhFpGLC5BfCF3yS+mUA7VR0BXBN4DvEMZh8JvKGqBcGuIyL9gH4ADRq4H9poqjdV5dlnn+WOO+6gfv36fPDBBzRp0uSw/Sq6bnWoXIabHIfNxjbJImyg8N6kuwGnq+rDItIAOElVV0dx3WzgW7/XW4GWDvvfClwGHCsiv1HV5wN3UNXxeCYGkpOTYwsrGfbs2UO/fv2YPn0611xzDdOmTaNOnTpB943lutVObIEik4zc5CieBS7EUwwQYBeV3P2jqmNVtbmq9g8WJIwJ9Nlnn9GyZUtmzpzJ8OHDeeWVV0IGCbB1q41x4qbrqaWqNhOR9QCqukNEjojyuoVAfb/Xp3q3GRO1efPm0bNnT2rWrMmbb77J5ZdfHvYYp3yDManO1ZrZ3uSzAohIPeBAlNddA5zpLTBYCOQDXaM8p0lxpaWl/O1vf+OJJ56gZcuWzJ49m/r164c/0Mu6hYwJzk3X01hgPnCCiAwH3gcec3sBEZkBfAA0FpGtItJHVUuBW4DFwCfAy6q6KeLWG+P13//+lz//+c888cQT3HTTTbz33nsRBQljTGhuigK+JCLrgD/jmRvUXlU/cXsBVe0SYvsiYJHb8xgTyvLly+ncuTNFRUW88MILXH/99a6PtTUhjAnPzainC4BNqjrO+7q2iLRU1VVxb50xISxYX8ioNz/l07dnsmPpJE7ObsCqVav4srQurUYucXXjt1XpjHHHTdfTc8Buv9e7vduMqRIL1hdy74xVfDRtGDuWTCDzNy05uvMTzPqSiGZXh1uVzhjj4SZQiPqtbqSqB6icGd3GBDVs2mK+mng7ezevIKt1T+q1v5/96bV48cP/RHTjr+gkO2NSjZsb/hYRuY1DTxE34VntzphKN2vWLD4adzOSUYsT8x+lVoM/hD2msKiYJsP+iQgU7S052CVVWZPsjEl2bp4o+gMX4RnG6ptB3S+ejTIm0P79+xk4cCD5+fkcfcoZnNzzKVdBwqeouIQde0vKdUldelY9m2RnjAthA4Wq/qCq+ap6gqqeqKpd/Yr0GRN33333HZdeeilPP/00t99+OxNffpVj6p4Y1TmLS8p499MfGdHhXLKzMhE8a0KM6HCuJbKNCeBm1FM9oC/Q0H9/Ve0dv2YZ47F06VLy8/PZvXs3M2fOpHNnT5X5jIwjGL14c9CuI7e+Kyq2SXbGuOCm6+kVPKvdvQ287vdlTNyoKqNGjeLPf/4zderUYfXq1QeDBHiGr64Y1IanOjep0PrWYLkIY9xyk8w+UlXvi3tLjPHauXMnvXr1Yv78+XTs2JGJEydyzDHHBN03VI0mcF7f2nIRxrgnGmbJRhF5FFjpnUmdFHJycnTt2rVV3QxTAR9//DG5ubls2bKF0aNHM3DgQDyV7iPnP+v62MyMg6Oe/L+32djGHCIi61Q157DtLgLFLuAoYL/3SwBV1drxaGgsWKBITi+99BJ9+/YlKyuLWbNmcfHFF8f8GoGzscHzdGFJbGNCBwo3o56OUdU0Va2lqrW9rxM2SJjks2/fPm655Rauv/56WrRoQUFBQVyCBNhsbGMqImygEI/rReRB7+v6ItIi/k0zqeDbb7/lT3/6E+PGjeOuu+7i7bff5qSTTorb9Ww2tjGRc5PMfhbP+hNtgEfw1HoaB5wfx3aZFPDWW2/RpUsX9u/fz5w5c8jNzXXcPxaVXm02tjGRczM8tqWq3gz8Cp4V7oBoV7gzKezAgQMMHz6cK664ghNPPJE1a9a4ChKRFPwLxZY8NSZybgJFPFa4Mylqx44dtGvXjsGDB5Ofn8+qVato3Dj8TTpWuYX2TbNtNrYxEXLT9RS4wl0eMDiurTLV0vr168nNzWXr1q0888wz3Hzzza6HvsYyt2CzsY2JjGOgEJE04CvgXiq4wp0xAJMnT+amm27iuOOO47333uPCCy+M6HjLLRhTdRy7nrxrT4xT1U9VdZyq/t2ChInEr7/+Sr9+/ejduzcXXXQRBQUFEQcJsNyCMVXJTdfTOyKSC8zTcLPzjPHz9ddfk5eXx7p16xg0aBCPPPIINWp4/smFmjV9SlYml55Vj3c//THo6CZb39qYyhfJzOwyoBibmW1cePPNN+nWrRtlZWVMnTqVdu3aHfxZsNnRTqp65nQshuUakwxiMTM7w2Zmm3AOHDjA0KFDueqqqzj11FNZu3ZtuSABwUcwOanKmdOxGpZrTDKzmdkmZrZv387VV1/NsGHD6N69Ox988AG/+c1vDtuvIiOVqmrmtJX8MMbdPIpngQuBrt7XvpnZxhy0du1amjdvzpIlS3j++eeZMmUKRx55ZNB9KzJSqapGN1nJD2NsZraJkqoyfvx4WrVqhaqyfPlybrzxRsf5EcFGMDmpytFNoQKUDcs1qcRmZpsKKy4upk+fPtx44420bt2adevW0aJFCxasL6TVyCU0GvQ6rUYuOaw/P3B2dFZmBnWOzDg4U/r6CxqU+1mtjDTumLUh6LnizYblGuNu1FM3oDPQDJiKd2a2qs6Of/MqxkY9xd+XX35Jbm4uH330EQ8++CBDhgwhPT096IgmwfMpIzvCEUOJsnaEjXoyqSLUqKeQ8yhEpJGqfqWqL4nIOmxmtvF69dVX6d69O2lpabz++utcddVVB38WLPnr+yjiGzEEuLrROiWSK/NGbSU/TKpz6nqaAyAi79jMbANQVlbGAw88wHXXXcfpp5/OunXrygUJCJ/kjWTEkCWSjUkMTjOz00TkfuC3InJn4A9V9X/j1yyTaH788Ue6du3K22+/TZ8+fXjmmWfIzDw8oRuqJpM/tzd6q+9kTGJweqLIxzMbuwZwTJAvkyJWrVpFs2bNWL58ORMmTGDChAlBgwS4G9Hk9kZviWRjEoPTE0VbVX1cRGqq6sOV1iKTMFSV5557joEDB5Kdnc3KlStp1qyZ4zH+NZkKi4oPJrJ9IrnRW30nYxJDyFFPIrJBVZuISIGqOt8dEoyNeorenj176N+/Py+++CJXXXUVL7zwAnXr1o34PDZiyJjkEfGoJ+ATEfkcOEVE/uV/LjxFAf8Q60aGIiJHAe8BQ1X1tcq6bqr67LPPyM3NZdOmTQwbNozBgweTluZmys3hbMSQMckvZKBQ1S4ichKwGLiuIicXkUnANcAPqnqO3/a2wNNAOjBBVUeGOdV9wMsVaYOJzPz58+nZsycZGRm8+eab/OUvf6nqJhljqpjjehSq+l/gvCjOPwX4OzDNt8E7y3sccDmwFVgjIgvxBI0RAcf39l7/30CtKNphwigtLeWBBx5g1KhRnH/++cyePZvTTjutqptljEkAThPuXlbVTiLyMeXzka67nlR1mYg0DNjcAvhCVbd4rzMTaKeqI/A8fQS2ozWe9TDOBopFZJF35T0TI9u2bSM/P5+lS5cyYMAAxowZQ82aNau6WcaYBOH0RHG798/Dbt5Ryga+9Xu9FWgZamdVfQBARHoCP4UKEiLSD+gH0KBBg1i1tdpbsWIFnTp1YseOHUybNo3u3btXdZOMMQnGKUfxvffPbyqvOaGp6pQwPx8PjAfPqKfKaFMyU1XGjh3L3XffTcOGDXnjjTf4wx8qbXyCMSaJOHU97aJ8l1M5UaxyVwjU93t9qnebqSS7d++mb9++zJw5k3bt2jF16lSOPfbYqm6WMSZBOT1RHAMgIo8A3wMv4MlPdANOjuKaa4AzRaQRngCRz6FFkUycffrpp3To0IHNmzczYsQI7r333goPfY2UzamofPZ3bmLBcdST13Wq6j/y6TkR+Qh4KNyBIjIDaA0cLyJbgSGqOlFEbsEz7DYdmKSqmyJvuonU7Nmz6d27N0ceeSRvv/02l156aczOHe6GFFgyPNJKsiZy9nduYsXNR8k9ItJNRNJFJM27PsUeNydX1S6qerKqZqjqqao60bt9kar+VlXPUNXh0bwBE15JSQl33nknnTp14txzz6WgoCDmQeJv8z6msKgY5dANyX+RIVt7uvLZ37mJFTeBoivQCdjm/eqIdRUlje+//542bdowZswYbr31VpYuXUp2dmw/Tbq5IVnJ8Mpnf+cmVsJ2Panq10C7+DfFxNqyZcvo1KkTu3btYvr06XTp0iUu13FzQ7KS4ZXP/s5NrFROFtNUKlXliSeeoE2bNhx77LGsXr06bkECQt94/LdbyfDKZ3/nJlbcJLNNEvnll1/o3bs3c+fOJTc3l0mTJlG7dkVHMjvzJbDdlBO3kuGVz/7OTayELDN+cAeRdFUtc9wpwaRqmfFNmzbRoUMHvvzySx5//HHuvPNORCQu1wocUQMcDBbZdkMyJilVpMy4z+ciMheYrKr/jn3TTCxMnz6dvn37Urt2bZYsWcIll1wSl+v4P0UE8gWJFYPaxOXaxpiq4SZHcR7wGTBBRD4UkX4iEp++DBOx2au/4uSL/odu3bqRfsLpDJ/6elyDhG8YbCg2osaY6idsoFDVXar6D1W9CM+6EEOA70Vkqoj8Ju4tNCFNeHMNf+1wFf/9YAHHnN+eOnmPMnr5D+XmL8RSsGGwgWxEjTHVT9iuJ+/6EVcDvYCGwJPAS8DFwCLgt3FsnwlhyZIlDMjLpaxkH8e3G8RRZ/0RODR/IXBWtH9C89Kz6vHupz86JjiDzbQO97RgI2qMqZ7cJLO3AO8CE1V1ZcDPxqrqbXFsX4VU52T2gQMHePzxxxk8eDDpdU6l3v/8jYzj6h+2ny+hDByWdA6UmZHOiA7nHgwWwRLVmRnp1MpIY8fekqDnsAS2MckvVDLbTaD4o6q+H7CtlaquiHEbY6a6BoqioiJ69OjBwoULyc/P56uzrue/Dh/yw93cA/lu9qGS1VmZGewrPXBYAPEPMsaY5BUqULhJZo8Nsu2Z6JtkIvHRRx+Rk5PDokWLePrpp5k+fTqDrmty2IQqf8UlZa6DBByq0RQqWV1UXELNGmnUOTIDwRNYLEgYU/05rUdxIXARUE9E7vT7UW08VV9NJZk6dSr9+/enbt26LF26lFatWgHlJ1Q5jUSKRHFJGekilIV40iwqLiEzI50xnZtYgDAmRTg9URwBHI0nmBzj9/ULkBf/ppl9+/bRv39/evbsyQUXXEBBQcHBIOHTvmk2Kwa1ITvEaKOszAzHp45gylTDPqlYBVJjUofTwkXvAe+JyJREWQ41lXzzzTd07NiRNWvWcN999/Hoo49So0boQWr3XNE4aAJ66HW/Bwg66snpKaRmjTTH/IbNlzAmdTh1PT2lqgOBv4vIYf0QqnpdPBuWyhYvXkzXrl0pLS1l/vz5tG/fPuwx4er6BOsmCja6ycfXxVTnyIygwcLmSxiTOpzmUbzg/fOJymiI8Qx9HT58OEOGDOGcc85h7ty5nHnmma6Pb980O6K8QbgcR3FJGTVrpJGZkX7Yk4rNlzAmdYQdHpuMknF47M8//8z111/PG2+8Qffu3Xn++ec58sgjK+36jQa9TrB/CQKM6dzEKpAakwIqXBRQRL6Cw+8hqnp6jNqW8goKCsjNzaWwsJBnn32W/v37IyJh16GOJadFbiJ9UjHGVC9uqsf6R5daeJZCrRuf5qSeiRMncvPNN3PCCSewfPlyWrZsCRyePygsKuaOWRsYOGsDWZkZiEDR3pKYBZBQyXDrYjLGuCkKuN3vq1BVn8JT+8lEobi4mD59+nDDDTdwySWXUFBQcDBIQPACfL7HuqLiEnbsLUE5NEku2kKA7ZtmM6LDuWRnZdpkOmNMOW66npr5vUzD84RhK+NFYcuWLeTl5bF+/XoeeOABhg0bRnp6+XkLkQw/DVYIsCKsi8kYE4ybG/6Tft+XAl8DneLSmhTw2muv0b17dwBeffVVrrnmmqD7hcoZhGLzGowx8eKm6+lSv6/LVbWvqtq03AiVlZXx4IMPcu2119KwYUPWrVsXMkiAJ2cQyYxqm9dgjIkXN11Pdzr9XFX/N3bNqZ5++uknunbtyltvvUWvXr0YN24cmZnON/bAOQ6+9aiDsaSzMSae3I56Oh9Y6H19LbAa+DxejapOVq9eTV5eHj/88AP/+Mc/uOGGG1wf658z8B8qe2wcRj0ZY0wobgLFqUAzVd0FICJDgddV9fp4NizZqSrPP/88t99+O6eccgrvv/8+OTmHzWM5TKi5E5ZoNsZUFTeB4kRgv9/r/d5tJoS9e/fSv39/XnjhBa688kpefPFF6tYNP/Uk2NyJv837GAheq8kYYyqDm4WLpgGrRWSo92liFTAlno1KZp9//jkXXHABL774IsOGDeO1115zFSQg+NwJK+ltjKlqYZ8oVHW4iLwBXOzd1EtV18e3WclpwYIF9OjRgxo1arBo0SLatm0b0fGhhrja0FdjTFVyKjNeW1V/EZG6eOZOfO33s7qq+nP8m5ccSktLGTx4MI8//jg5OTnMmTOH0047zfGYYLkIp3pLkarMOlHGmOrN6YliOnANsI7yIzN9IzWtKCCwbds2unTpwrvvvkvfvn0ZO3YstWrVOvjzYDdsIGguIrd5NnPXFUZdb8lyHcaYWLIy41FYuXIlHTt25Oeff6bvoMdYV/M8x4AAnht/qJXjsr3HRfsk0GrkkqBPJtlZmawY1CbCd2mMSRURlxkPqPF0GFUtiEXDkpGq8swzz3DXXXfRoEEDHpuykPEbyygu9tycfZ/ga2WkBU1OB1tRDjy5iFgMg7VchzEmlpy6np50+JkClfLRVETSgEeA2sBaVZ1aGdcNZffu3fTr148ZM2Zw7bXXMm3aNK5+voDikv3l9nMKCKEonqeBaPMJscx1GGNMyEChqpdGe3IRmYQnz/GDqp7jt70t8DSQDkxQ1ZEOp2mHZ9LfdmBrtG2KxubNm8nNzeWTTz6h2833seWUy2g6ckXI0hqhZGVmsK/0QNBAEot8gq0tYYyJpbDzKESklojcKSLzRGSuiAwUkVrhjvOaApQbIyoi6cA44ErgbKCLiJwtIueKyGsBXycAjYGVqnonMCCSNxdLc+bMIScnh23btvHQuJdYV6c13/2yzzFIZGVmHFbYLzMjnaHX/f7g2g/BRDt3wtaWMMbEkpuZ2dOAXcAz3tddgRfwrHTnSFWXiUjDgM0tgC9UdQuAiMwE2qnqCDxPH+WIyFYOzQwP2ZcjIv2AfgANGjQI1zTX5qz+mtvvvpfvls/m6Pq/Y8T/TWXyR7soLnHu7/cFBCBkcrp90+yQa1VHm0+wkh/GmFhxEyjOUdWz/V6/KyL/juKa2cC3fq+3Ai1D7AswD3hGRC4GloXaSVXHA+PBM+opivYdNOmfBdxyQw+Kv93IMc2upk6bG3hy5XbH3IPAwaJ9d8zaEHbkkuUTjDGJzk0JjwIRucD3QkRaAvEfe+qlqntVtY+q3qqq4yrrusuXL6d/3uX8+v3nHHfNXdS9fACSnkFxSRnpIkGPyc7KZEznJuwrPeB6qdJg605YPsEYk0ichsd+jGcgTgawUkT+4/1RA+DTKK5ZCNT3e32qd1uV8k2MK9yxF9n4Ov9Z/A/Sap/ISX99mCPqNSy3b5kqmRnpQZPFTvWagj1V+K87YbOojTGJyKnrKfTya9FZA5wpIo3wBIh8PHmPKuObybxn9y62v/E0ezev4OjGF9Gow938cuCIw/Z3mhh3x6wNQa/hlHOwfIIxJpE5DY/9xve9iJzHoaKAy1X1IzcnF5EZQGvgeG9SeoiqThSRW4DFeIbHTlLVTRVsf0z4ngJ+WjSG4s9XkdW6F7VbdCCt5hFkBgxj9T05hLq5W87BGFPduBkeezvwEnCC9+tFEbnVzclVtYuqnqyqGap6qqpO9G5fpKq/VdUzVHV4NG8gFnyf9utc8ldOzH+UY1vmIiLsLC6JeJip5RyMMdWNm1FPfYCWqroHQEQeBz7g0HDZpOd7Csg4rj4Zx9Uvtz3YcqROo5ks52CMqW7cBAqh/PyFMu+2asPNTOZIKrJazsEYU524CRSTgVUiMt/7uj0wMW4tqgJungIiHc1kjDHVhZsV7v5XRJYCf/RuqpYr3IV7CrCKrMaYVOXmicJXUjxly4qDjWYyxqQuV4Ei1QRblc4qshpjUpWbEh4pxZe0LiwqLleCA7CKrMaYlGRPFAGcktYrBrWxwGCMSTn2RBHAktbGGFOeBYoAoZLTlrQ2xqQqCxQBrASHMcaUZzmKAFaCwxhjyrNAEYSV4DDGmEOs68kYY4wje6KoRMEm8tmTizEm0VmgqCSRVJ81xphEYl1PlcRpIp8xxiQyCxSVxCbyGWOSlQWKSmIT+YwxycoCRSWxiXzGmGRlyexKYhP5jDHJygJFJbKJfMaYZGRdT8YYYxxZoDDGGOPIAoUxxhhHFiiMMcY4skBhjDHGkQUKY4wxjixQGGOMcWSBwhhjjCMLFMYYYxxZoDDGGOPIAoUxxhhHCV/rSUQaAGOBn4HPVHVkFTfJGGNSSlyfKERkkoj8ICIbA7a3FZHNIvKFiAwKc5pzgTmq2htoGrfGGmOMCSreTxRTgL8D03wbRCQdGAdcDmwF1ojIQiAdGBFwfG/gQ2COiPQGXohze40xxgSIa6BQ1WUi0jBgcwvgC1XdAiAiM4F2qjoCuCbwHCJyNzDEe645wOR4ttkYY0x5VZGjyAa+9Xu9FWjpsP+bwFAR6Qp8HWonEekH9PO+3C0imyvYvuOBnyp4bKKx95J4qsv7AHsviSja93FasI0Jn8xW1Y1Anov9xgPjo72eiKxV1Zxoz5MI7L0knuryPsDeSyKK1/uoiuGxhUB9v9enercZY4xJQFURKNYAZ4pIIxE5AsgHFlZBO4wxxrgQ7+GxM4APgMYislVE+qhqKXALsBj4BHhZVTfFsx0Rirr7KoHYe0k81eV9gL2XRBSX9yGqGo/zGmOMqSashIcxxhhHFiiMMcY4skDhJ8LSIglLRL4WkY9FZIOIrK3q9kQiWNkXEakrIm+JyOfeP+tUZRvdCvFehopIofd3s0FErqrKNrohIvVF5F0R+beIbBKR273bk+734vBekvH3UktEVovIR973Msy7vZGIrPLex2Z5Bw1Fdy3LUXh4S4t8hl9pEaCLqv67ShtWASLyNZCjqkk3gUhELgF2A9NU9RzvtlHAz6o60hvA66jqfVXZTjdCvJehwG5VfaIq2xYJETkZOFlVC0TkGGAd0B7oSZL9XhzeSyeS7/ciwFGqultEMoD3gduBO4F5qjpTRJ4HPlLV56K5lj1RHHKwtIiq7gdmAu2quE0pR1WX4akU7K8dMNX7/VQ8/7ETXoj3knRU9XtVLfB+vwvPaMVskvD34vBeko567Pa+zPB+KdAGmOPdHpPfiwWKQ4KVFknKf0B4/rH8U0TWeUubJLsTVfV77/f/BU6sysbEwC0i8i9v11TCd9f489ZuawqsIsl/LwHvBZLw9yIi6SKyAfgBeAv4EijyTkOAGN3HLFBUT39U1WbAlcDN3i6QakE9faXJ3F/6HHAG0AT4HniySlsTARE5GpgLDFTVX/x/lmy/lyDvJSl/L6papqpN8FS4aAGcFY/rWKA4pNqUFlHVQu+fPwDz8fwDSmbbvH3Lvj7mH6q4PRWmqtu8/7kPAP8gSX433j7wucBLqjrPuzkpfy/B3kuy/l58VLUIeBe4EMgSEV8dv5jcxyxQHFItSouIyFHeJB0ichTwF2Cj81EJbyHQw/t9D+CVKmxLVHw3Vq//IQl+N96k6UTgE1X9X78fJd3vJdR7SdLfSz0RyfJ+n4lnIM4neAKGr5BqTH4vNurJj3dI3FN4FlGapKrDq7ZFkROR0/E8RYCnOvD0ZHof3rIvrfGUS94GDAEWAC8DDYBvgE6qmvBJ4hDvpTWe7g3FUzb/Rr9+/oQkIn8ElgMfAwe8m+/H07efVL8Xh/fSheT7vfwBT7I6Hc+H/pdV9WHvPWAmUBdYD1yvqvuiupYFCmOMMU6s68kYY4wjCxTGGGMcWaAwxhjjyAKFMcYYRxYojDHGOLJAYVKOiGSJyE1+r08RkTlOx1TwOr6KpA9HeFxPEfl7iJ+t9P7ZUES6hjlPprcS6n4ROT6SNhjjzwKFSUVZwMFAoarfqWpe6N2jMkZVH3K7s9+M2qBU9SLvtw0Bx0ChqsXe8g7fub2+McFYoDCpaCRwhvfT9mjvp/ONcPDT/ALv+gpfi8gtInKniKwXkQ9FpK53vzNE5E1v4cXlIhK2xo53/YYF3sJzH3onTPmePF4QkRXAC97d64vIUvGs9TDE7xy+aqEjgYu97+EOEfm9d22CDd7znxnDvy+T4hw/vRhTTQ0CzvF+2vZVEfV3Dp6qorWAL4D7VLWpiIwB/opn9v54oL+qfi4iLYFn8ZR3djIMWK+q7UWkDTANz2xggLPxFHMsFpGeeGoNnQPsBdaIyOuq6r8I1SDgblW9xvsengGeVtWXvCVo0iP4+zDGkQUKYw73rnetgl0ishN41bv9Y+AP3sqjFwGzPaWDAKjp4rx/BHIBVHWJiBwnIrW9P1uoqsV++76lqtsBRGSe91in1Qo/AB4QkVPxLFrzuYv2GOOKdT0Zczj/ujgH/F4fwPPhKg1Pzf8mfl+/i/KaewJeB9bWcay1o6rTgeuAYmCR94nFmJiwQGFS0S7gmIoe7F2/4CsR6QieiqQicp6LQ5cD3bzHtAZ+ClzXwc/l3pxGJp4VylYE/Lzce/AWgtuiqmPxVAv9g+s3ZEwYFihMyvF26awQkY0iMrqCp+kG9BGRj4BNuFs2dyjQXET+hScZ3cNh39V41kz4FzA3ID+Bd3uZiHwkInfgWfN5o3e1s3Pw5D+MiQmrHmtMnIjIUGC3qj5Rxe34GshR1Z+qsh0medkThTHxsxvoF+mEu1jxTbgDMji09oIxEbMnCmOMMY7sicIYY4wjCxTGGGMcWaAwxhjjyAKFMcYYRxYojDHGOPp/stT5eQpjkGUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,1)\n",
    "ax.set_xlabel(\"time [orbits]\")\n",
    "ax.set_ylabel(\"obliquity difference [degrees]\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_ylim([1e-8,1e2])\n",
    "o1 = np.remainder((obliq-obliq2)*180/np.pi,360.)\n",
    "o2 = np.remainder((obliq2-obliq)*180/np.pi,360.)\n",
    "ax.scatter(times/np.pi/2.,np.minimum(o1,o2))\n",
    "ax.plot(times/np.pi/2.,1e-8/np.pi*180.*np.exp(times/(np.pi*2.)/1.2),color=\"black\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2bc381f",
   "metadata": {},
   "source": [
    "On a log-linear scale, we see that the divergence follows a straight line, indicating exponential growth. The Lyapunov timescale is approximately 1.2 orbits. About 25 days! (Wikipedia says ~30 days which is close enough given the very simplistic model)"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}