examples/2_advanced_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/stammler/simframe/HEAD?labpath=examples%2F2_advanced_integration.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This notebook teaches how to:**  \n",
    "create groups, create immutable fields, display the table of contents, and set up complex integration instructions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Advanced Integration\n",
    "\n",
    "Let's have a look at orbital integration.\\\n",
    "Consider a planet Earth orbiting the Sun. How does its orbit look like?\n",
    "\n",
    "The equations of motion are determined by a set of two differential equations for the position $\\vec{r}$ and velocity $\\vec{v}$ of Earth and Sun..\n",
    "\n",
    "$\\frac{\\mathrm{d}}{\\mathrm{d}t} \\vec{r} = \\vec{v}$\n",
    "\n",
    "$m\\frac{\\mathrm{d}}{\\mathrm{d}t} \\vec{v} = \\vec{F}_\\mathrm{G}$\n",
    "\n",
    "The gravitational force $F_\\mathrm{G}$ of a body of mass $M$ and position $\\vec{R}$ acting on a body of mass $m$ at position $\\vec{r}$ is given by\n",
    "\n",
    "$\\vec{F}_\\mathrm{G} = -GmM\\frac{\\vec{r}-\\vec{R}}{\\left| \\vec{r}-\\vec{R} \\right|^3}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating groups\n",
    "\n",
    "First we need to add fields for postion, velocity, and mass for Earth and Sun.\\\n",
    "To have more structure we group the fields by Earth and Sun. This can be done with `addgroup(name)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simframe import Frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = Frame(description=\"Earth-Sun system\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.addgroup(\"Sun\")\n",
    "sim.addgroup(\"Earth\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The frame object has now two groups for Earth and Sun, that can be adressed just as fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Frame (Earth-Sun system)\n",
       "------------------------\n",
       "    Earth        : Group\n",
       "    Sun          : Group\n",
       "  -----\n",
       "    Integrator   : \u001b[93mnot specified\u001b[0m\n",
       "    Writer       : \u001b[93mnot specified\u001b[0m"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will keep the problem general in three dimensions, but we will only use two of them for simplicity.\\\n",
    "Although spherical coordinates would be more suitable for the problem, we will use Carteesian coordinates here.\n",
    "\n",
    "First, we need to define a few constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "AU      = 1.495978707e11        # Astronomical unit [m]\n",
    "day     = 8.64e4                # Day [s]\n",
    "G       = 6.6743e-11            # Gravitational constant [m³/kg/s²]\n",
    "year    = 3.15576e7             # Year [s]\n",
    "\n",
    "M_earth = 5.972167867791379e24  # Mass of the Earth [kg]\n",
    "M_sun   = 1.988409870698051e30  # Mass of the Sun [kg]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating immutable fields\n",
    "\n",
    "We can now fill the groups with the variables we need, starting with the masses.  \n",
    "Again: the initial values must have the correct shape and data types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.Earth.addfield(\"M\", M_earth, description=\"Mass [kg]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The groups and the fields within can be easily accessed with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group\n",
       "-----\n",
       "    M            : Field (Mass [kg])\n",
       "  -----"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Earth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.972167867791379e+24"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Earth.M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mass of the Earth shall be constant throughout the simulation. We therefore set a flag so we cannot accidentally change its value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.Earth.M.constant = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group\n",
       "-----\n",
       "    M            : Field (Mass [kg]), \u001b[95mconstant\u001b[0m\n",
       "  -----"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Earth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we add the mass of the Sun, which we directly set to constant when adding the field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.Sun.addfield(\"M\", M_sun, constant=True, description=\"Mass [kg]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group\n",
       "-----\n",
       "    M            : Field (Mass [kg]), \u001b[95mconstant\u001b[0m\n",
       "  -----"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Sun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do all sorts of operations with those fields, such as calculating the mass ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0034893488507934e-06"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Earth.M / sim.Sun.M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to add fields for position and velocity. We initialize them with zeros and think later about the initial conditions.\\\n",
    "The important thing is that the have the right shape upon initialization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.Earth.addfield(\"r\", np.zeros(3), description=\"Position [m]\")\n",
    "sim.Earth.addfield(\"v\", np.zeros(3), description=\"Velocity [m/s]\")\n",
    "sim.Sun.addfield(\"r\", np.zeros(3), description=\"Position [m]\")\n",
    "sim.Sun.addfield(\"v\", np.zeros(3), description=\"Velocity [m/s]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group\n",
       "-----\n",
       "    M            : Field (Mass [kg]), \u001b[95mconstant\u001b[0m\n",
       "    r            : Field (Position [m])\n",
       "    v            : Field (Velocity [m/s])\n",
       "  -----"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Earth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group\n",
       "-----\n",
       "    M            : Field (Mass [kg]), \u001b[95mconstant\u001b[0m\n",
       "    r            : Field (Position [m])\n",
       "    v            : Field (Velocity [m/s])\n",
       "  -----"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.Sun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setting the integration variable**\n",
    "\n",
    "For our simulation we need an integration variable. In our case this is the time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.addintegrationvariable(\"t\", 0., description=\"Time [s]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll set the step size to a constant value of one day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 1.*day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_dt(frame):\n",
    "    return dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.t.updater = f_dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to integrate for two years and want to have a snapshot every ten days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "snapwidth = 10.*day\n",
    "tmax = 2.*year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.t.snapshots = np.arange(snapwidth, tmax, snapwidth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** If the initial value of the integration variable is smaller the first snapshot, `simframe` automatically writes an output with initial conditions, if a writer is set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Frame (Earth-Sun system)\n",
       "------------------------\n",
       "    Earth        : Group\n",
       "    Sun          : Group\n",
       "  -----\n",
       "    t            : IntVar (Time [s]), \u001b[95mIntegration variable\u001b[0m\n",
       "  -----\n",
       "    Integrator   : \u001b[93mnot specified\u001b[0m\n",
       "    Writer       : \u001b[93mnot specified\u001b[0m"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Setting the writer**\n",
    "\n",
    "As a writer we use the namespacewriter, which does not write the data into files, but stores them within a buffer in the writer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simframe import writers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.writer = writers.namespacewriter()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The namespacewriter is by default not writing dump files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Writer (Temporary namespace writer)\n",
       "-----------------------------------\n",
       "    Data directory : data\n",
       "    Dumping        : \u001b[93mFalse\u001b[0m\n",
       "    Verbosity      : 1"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.writer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Adding differential equations**\n",
    "\n",
    "As a next step we'll add differential equations to the quantities. The differential equations for the positions are simple. We simply return the velocities, which we can adress with the `frame` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dr_Earth(frame, x, Y):\n",
    "    return frame.Earth.v\n",
    "\n",
    "def dr_Sun(frame, x, Y):\n",
    "    return frame.Sun.v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the differential equations of the velocities we'll write a little helper function that computes the gravitational accelleration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gravitational acceleration\n",
    "def ag(M, r, R):\n",
    "    direction = r-R\n",
    "    distance = np.linalg.norm(direction)\n",
    "    return -G * M * direction / distance**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dv_Earth(frame, x, Y):\n",
    "    return ag(frame.Sun.M, frame.Earth.r, frame.Sun.r)\n",
    "\n",
    "def dv_Sun(frame, x, Y):\n",
    "    return ag(frame.Earth.M, frame.Sun.r, frame.Earth.r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to add the differential equations to their fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.Earth.v.differentiator = dv_Earth\n",
    "sim.Earth.r.differentiator = dr_Earth\n",
    "sim.Sun.v.differentiator   = dv_Sun\n",
    "sim.Sun.r.differentiator   = dr_Sun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Displaying the table of contents\n",
    "\n",
    "You can also display the complete tree structure of your frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[94mFrame (Earth-Sun system)\u001b[0m\n",
      "    - \u001b[94mEarth\u001b[0m: \u001b[94mGroup\u001b[0m\n",
      "        - \u001b[94mM\u001b[0m: Field (Mass [kg]), \u001b[95mconstant\u001b[0m\n",
      "        - \u001b[94mr\u001b[0m: Field (Position [m])\n",
      "        - \u001b[94mv\u001b[0m: Field (Velocity [m/s])\n",
      "    - \u001b[94mSun\u001b[0m: \u001b[94mGroup\u001b[0m\n",
      "        - \u001b[94mM\u001b[0m: Field (Mass [kg]), \u001b[95mconstant\u001b[0m\n",
      "        - \u001b[94mr\u001b[0m: Field (Position [m])\n",
      "        - \u001b[94mv\u001b[0m: Field (Velocity [m/s])\n",
      "    - \u001b[94mt\u001b[0m: IntVar (Time [s]), \u001b[95mIntegration variable\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sim.toc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting memory usage\n",
    "\n",
    "It is possible to print out the memory usage in bytes of a group with the command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "341.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.memory_usage()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further details can be displayed with the keyword argument `print_output`, while `skip_hidden` will ommit attributes starting with underscore `_`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    - \u001b[94mEarth\u001b[0m:                                   total:  56   B\n",
      "        - \u001b[94mM\u001b[0m:                                     (1,)   8   B\n",
      "        - \u001b[94mr\u001b[0m:                                     (3,)  24   B\n",
      "        - \u001b[94mv\u001b[0m:                                     (3,)  24   B\n",
      "    - \u001b[94mSun\u001b[0m:                                     total:  56   B\n",
      "        - \u001b[94mM\u001b[0m:                                     (1,)   8   B\n",
      "        - \u001b[94mr\u001b[0m:                                     (3,)  24   B\n",
      "        - \u001b[94mv\u001b[0m:                                     (3,)  24   B\n",
      "    - \u001b[94mt\u001b[0m:                                         (1,)   8   B\n",
      "\n",
      "Total: 120   B\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "120.0"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.memory_usage(print_output=True, skip_hidden=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up complex integration instructions\n",
    "\n",
    "Next we need to set up the integrator. We integrate all quantities with the explicit Euler 1st-order scheme as in the previous tutorial, but this time he have to integrate four fields. This can be easily achieved by adding four instructions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simframe import Integrator\n",
    "from simframe.integration import Instruction\n",
    "from simframe import schemes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator = Integrator(sim.t, description=\"Euler 1st-order\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "instructions_euler = [Instruction(schemes.expl_1_euler, sim.Earth.r),\n",
    "                      Instruction(schemes.expl_1_euler, sim.Earth.v),\n",
    "                      Instruction(schemes.expl_1_euler, sim.Sun.r  ),\n",
    "                      Instruction(schemes.expl_1_euler, sim.Sun.v  ),\n",
    "                     ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator.instructions = instructions_euler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initial conditions**\n",
    "\n",
    "Before we can start the simulation, we have to think about initial conditions.\n",
    "\n",
    "If we simply set the Sun at rest in the center of our simulation and only set the position and velocity of the Earth, then the center of mass would have a non-zero momentum and would slowly drift away from its initial position.\n",
    "\n",
    "So what we do first: we set the Sun's positon to zero (i.e., don't do anything) and the Earth's position to a distance of 1 AU in positive x-direction. Then we'll offset their positions to center the system on the center of mass instead onto the Sun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_Earth_ini = np.array([AU, 0., 0.])\n",
    "r_Sun_ini = np.zeros(3)\n",
    "# Center of mass\n",
    "COM_ini = (M_earth*r_Earth_ini + M_sun*r_Sun_ini) / (M_earth+M_sun)\n",
    "# Offset both positions\n",
    "r_Earth_ini -= COM_ini\n",
    "r_Sun_ini -= COM_ini"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We save them in a separate variable instead of assigning them directly for later use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inital orbital velocities of the Earth shall be in positive y-direction, the velocity of the Sun in negative y-direction.\\\n",
    "For calculating the value we use the reduced mass $mu$ of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = M_earth*M_sun / (M_earth+M_sun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "v_Earth_ini = np.array([0., np.sqrt(G*M_sun/M_earth*mu/AU), 0.])\n",
    "v_Sun_ini   = np.array([0., -np.sqrt(G*M_earth/M_sun*mu/AU), 0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we assign them to their fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.Earth.r = r_Earth_ini\n",
    "sim.Earth.v = v_Earth_ini\n",
    "sim.Sun.r   = r_Sun_ini\n",
    "sim.Sun.v   = v_Sun_ini"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Starting the simulation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving frame 0000\n",
      "Saving frame 0001\n",
      "Saving frame 0002\n",
      "Saving frame 0003\n",
      "Saving frame 0004\n",
      "Saving frame 0005\n",
      "Saving frame 0006\n",
      "Saving frame 0007\n",
      "Saving frame 0008\n",
      "Saving frame 0009\n",
      "Saving frame 0010\n",
      "Saving frame 0011\n",
      "Saving frame 0012\n",
      "Saving frame 0013\n",
      "Saving frame 0014\n",
      "Saving frame 0015\n",
      "Saving frame 0016\n",
      "Saving frame 0017\n",
      "Saving frame 0018\n",
      "Saving frame 0019\n",
      "Saving frame 0020\n",
      "Saving frame 0021\n",
      "Saving frame 0022\n",
      "Saving frame 0023\n",
      "Saving frame 0024\n",
      "Saving frame 0025\n",
      "Saving frame 0026\n",
      "Saving frame 0027\n",
      "Saving frame 0028\n",
      "Saving frame 0029\n",
      "Saving frame 0030\n",
      "Saving frame 0031\n",
      "Saving frame 0032\n",
      "Saving frame 0033\n",
      "Saving frame 0034\n",
      "Saving frame 0035\n",
      "Saving frame 0036\n",
      "Saving frame 0037\n",
      "Saving frame 0038\n",
      "Saving frame 0039\n",
      "Saving frame 0040\n",
      "Saving frame 0041\n",
      "Saving frame 0042\n",
      "Saving frame 0043\n",
      "Saving frame 0044\n",
      "Saving frame 0045\n",
      "Saving frame 0046\n",
      "Saving frame 0047\n",
      "Saving frame 0048\n",
      "Saving frame 0049\n",
      "Saving frame 0050\n",
      "Saving frame 0051\n",
      "Saving frame 0052\n",
      "Saving frame 0053\n",
      "Saving frame 0054\n",
      "Saving frame 0055\n",
      "Saving frame 0056\n",
      "Saving frame 0057\n",
      "Saving frame 0058\n",
      "Saving frame 0059\n",
      "Saving frame 0060\n",
      "Saving frame 0061\n",
      "Saving frame 0062\n",
      "Saving frame 0063\n",
      "Saving frame 0064\n",
      "Saving frame 0065\n",
      "Saving frame 0066\n",
      "Saving frame 0067\n",
      "Saving frame 0068\n",
      "Saving frame 0069\n",
      "Saving frame 0070\n",
      "Saving frame 0071\n",
      "Saving frame 0072\n",
      "Saving frame 0073\n",
      "Execution time: \u001b[94m0:00:00\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sim.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Reading and plotting**\n",
    "\n",
    "Reading data from the namespace writer works identical to the writer discussed earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_euler = sim.writer.read.all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_orbits(data):\n",
    "    fig, ax = plt.subplots(dpi=150)\n",
    "    ax.set_aspect(1)\n",
    "    ax.axis(\"off\")\n",
    "    fig.set_facecolor(\"#000000\")\n",
    "    imax = data.t.shape[0]\n",
    "    for i in range(imax):\n",
    "        alpha = np.maximum(i/imax-0.1, 0.5)\n",
    "        ax.plot(data.Sun.r[i, 0], data.Sun.r[i, 1], \"o\", c=\"#FFFF00\", markersize=4, alpha=alpha)\n",
    "        ax.plot(data.Earth.r[i, 0], data.Earth.r[i, 1], \"o\", c=\"#0000FF\", markersize=1, alpha=alpha)\n",
    "    ax.plot(data.Sun.r[-1, 0], data.Sun.r[-1, 1], \"o\", c=\"#FFFF00\", markersize=16)\n",
    "    ax.plot(data.Earth.r[-1, 0], data.Earth.r[-1, 1], \"o\", c=\"#0000FF\", markersize=4)\n",
    "    ax.set_xlim(-1.5*AU, 1.5*AU)\n",
    "    ax.set_ylim(-1.5*AU, 1.5*AU)\n",
    "    fig.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAJJCAYAAAC+gKM0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAABcSAAAXEgFnn9JSAAAacklEQVR4nO3dX4ydd53f8Y9RBDfthm5C1jHpYk52SEISQtb5n2WB3iz4kqutkmoHrZBWoW3UdqtuK5KqNWq36j/1D9FK0YpByt72EtKbwrYFYoObACaBjHLoRWS7bLIl297ARdyL50w9k3zteDznnN/z5/W6eTQQfL4m9sz7/H6/5zmHklwIAAB7vKv1AAAAfSSSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAK17QeAJiajya5PckPk7zQdBKAyxFJwBp9NMkfJDmc5HySP4xQAvrKdhtM2t1J/sbiug63pwukX11cb1/T696V5K8vrgBXxkoSTNbdSf5RLq7q/LMkz6/4NX+4eK0srj9c8eslXRj9g1z8ff6LJN9bw+sCQyeSYLLuSBcOf3XX16uOpBfSbbGt80zSh7P39/nhiCTgSogkmKwz2buqc2ZNr/tC1nsO6cXs/X2+uMbXBobsUJILrYcAWrk73QrSmax+Famlu9KtIL0Yq0jAlRJJ0Bu/nuTOJD9I8j8az8LBfSTJbUleSvL9xrMAV8N2G/TCryd5IhcPF5+IUBqyjyT5/Vz89/mvIpRgeDwCAHrhznQ/UN+/uN7ZdhwO6LbsfdTBbW3HAa6KlSTohR9k7+HiHzSchYN7KXv/fb7UcBbgajmTBL3hTNK4OJMEQyeSAAAKziQBABREEgBAQSQBABTc3QaXdCzd4dvvJzndeBZ4J3cmuTXJj+LuSFgOkQSlY0n+cZJfSfK/kvyTCCX6684kfzfJjUnOJfk3EUpwcLbboPSRdIH0/sX1I23Hgcu6NV0g/eriemvbcWAkrCRB6fvpVpCyuHrODX32o3QrSFlcf9RwFhgPz0mCS3ImiSFxJgmWTSQBABScSQIAKIgkAICCSAIAKIgkAICCSAIAKIgkAICCSAIAKHjiNgN0T5K7knwvyXcbzwJTc0eSW5L8OMmZxrPAaokkBuaeJP80yQ1JfprkyQglWJc7kjyeix+k++8ilBgz220MzF3pAun9i+tdbceBSbklXSB9YHG9pe04sGJWkhiY76VbQcri+r2Gs8DU/Dh7P0j3xw1ngdXz2W0MkDNJ0I4zSUyHSAIAKDiTBABQcCYJuAIbSWZJ5km2J/j6wBSJJOAdbCTZTHIkydkkW1lvqLR+fWCqbLcB72CWLlCOLq6zib0+MFVWkoB3ME+3gpPFdT6x1wemyt1twBVofSao5eu3/r0DrYgkgEtyHgqmzJkkGIyNJL+1uLIezkPBlDmTBINgRaMN56FgykQSK3Bvko8meSHJd5pOMh67VzR2vhZJq7edLkidSVqv25N8KMnLSX7YeBamTCSxZPcm+WKSG9J9AO0XIpSWwYpGO9sRR+t0e5K/leTGdB+i+x8ilGhFJLFkH00XSEd2fS2SDs6KBlPxoXSB9IFdX4sk2hBJLNkL6VaQsri+0GyS8bGiwRS8nG4FKYvryw1nYeo8AoAVcCYJOAhnkugHkQQAULDdBtBbnvYNLYkkgF7ybCxozRO3Yd88+Zp18LRvaM1KEuyLd/esi2djQWsiCfbFk69ZF8/GgtZEEuyLd/esk2djQUseAQD75o4jgCkQSQAABXe3AQAURBIAQEEkAQAURBIAQMEjACbrviR3J3k+yanGswBA/4ikSbovyT9P8r4kf5bkH0YoAeP34XSP8NhO8mLjWRgCkTRJd6cLpBt3fS2SgDH7cJLH0n3fO5fkqQgl3olImqTn060gZXF9vuEsQH+M+UGpG+kC6eiur0USlyeSJulUui02Z5KAHWP/8ObtdCtIWVzH9HtjVUTSZJ3KuOJozO+AYR3G/uHNL6bbYnMmiSsnkhiBsb8DhnWYwoc3vxhxxH6IJEZg7O+AYR22073BsCILO0QSIzCFd8CwDtsRR3DRoSQXWg8BB+dMEgDLJZIAAAo+uw0AoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKPuB2UO5PcneS55OcbDwLAIybSBqM+5P8YZL3JfmzJH8QoQQMkw+kZhhE0mDcnS6QDu/6WiQBQ7ORZDPJkSRnk2xluKF0W7rfz3aSlxrPwiqIpMF4Pt0KUhbX5xvOsl/eNQI7ZukC6eiur4f4feG2JL+Xi7H3RxFK4yOSBuNkui22oZ1JGtO7RuDg5um+F2RxnTec5SA20n1f+8Cur0XS2IikQTmZ4cTRjrG8awSWYzvdm6Whry5vZ2/sDfX3weWIJFZsLO8ageXZzvCj4qV0W2zOJI3ZoSQXWg/B2DmTBMDwiCQAgIInbgMAFEQSAEBBJAEAFEQSAEBBJAEAFEQSAEBBJAEAFEQSAEBBJAEAFEQSAEBBJAEAFEQSAEBBJAEAFEQSAEDhmtYDTMMDSY4lOZ3kucazAABXQiSt3ANJ/mWS65O8luTvRygBLNtGklmSeZLtxrMwFrbbVu5YukA6vLgeazsOwOhsJNlM8tuL60bLYRgRK0krdzrdClIW19MNZ3kn3okBQzRLciTJ0V1fD+V72K1Jbk7ySpIfNZ6FtxJJK/dcui22vp9J2nkndiTJ2SRbGc43GWDa5um+b2VxnTecZT9uTfK5XPy++3SEUr+IpLV4Lv2Nox1DficGTNt2ujd2Q1sJvznd990P7PpaJPWJSGJhqO/EAJIujIYSRzteyd7vu680nIXKoSQXWg9BXziTBLBeziT1mUgCACh4BAAAQEEkAQAURBIAQEEkAQAURBIAQEEkAQAURBIAQEEkAQAURBIAQEEkAQAURBIAQEEkAQAURBIAQEEkAQAURBIAQOGa1gMMy4NJjiU5neTbjWcBAFZJJF2xB5P86yTXJXk9yd+LUAKA8bLddsWOpQukX1lcj7UdBwBYKStJV+x0uhWkLK6nG84CAKzaoSQXWg8xHH0+k7SRZJZknmS78SwAQzWU76W3JLk5yStJftx4lvESSaOwkWQzyZEkZ5Nspd9/uQH6aCjfS29J8ru5OOcfRyithjNJozBL95fl6OI6azoNwDAN5Xvpzdk7581NpxkzZ5JGYZ7u3UQW13nDWQCGaijfS1/J3jlfaTjLuNluG42h7KMD9NlQvpc6k7QOIgkAoOBMEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABSuaT1Aew8luSfJd5N8q/EsAEBfTDySHkryb5P8cpI/T/J3IpQAgGTy2233pAukGxbXe9qOAwD0xsRXkr6bbgUpi+t3G84CAPTJoSQXWg/RljNJAMDbiSQAgMLEt9v6aCPJLMk8yXbjWQAYtg/l4s+UlxvPMjwiqVc2kmwmOZLkbJKtCCWAvuvrm9sPJflsLv5M+XKE0v6IpF6ZpfvDfHTX1336CwfAXn1+c1v9TBFJ+yGSemWe7i9ZFtd5w1kAeGd9fnPrZ8pBiaRe2U73LqSPy7YAvF2fQ+TldFtsziRdLXe3AcCB9PVMEgclkgAAChP/WBIAgJpIAgAoiCQAgIJIAgAoiCQAgIJIAgAoiCQAgIJIAgAoiCQAgIJIAgAoiCQAgIJIAgAoiCQAgIJIAgAoXNN6gNV7OMm9Sb6T5JuNZwEAhmLkkfRwkn+f5K8k+d9J/naEEgBwJUa+3XZvukB63+J6b9txAIDBGPlK0nfSrSBlcf1Ow1kAgCE5lORC6yFWy5kkAGD/JhBJAAD7N/IzSQAAV2fkZ5IAgLfbSDJLMk+y3XiW/hJJa+UPJQCtbSTZTHIkydkkW/EzqSaS1sYfSgD6YJbuZ9HRXV/7eVQRSWvjDyUAfTBP92Y9i+u84Sz9JpLWxh9KgGnq21GL7XS7GX2aqZ88AmCt+vYXBYDVctRiyKwkrdV2/OUAmBJHLYZMJAHAyjhqMWS22wBgpRy1GCqRBABQ8LEkAAAFkQQAUBBJAAAFkQQAUBBJAAAFkQQAUBBJAAAFkQQAUBBJAAAFn90G9MLhw8n11yfvfnfyi18kr72WnD/feipgykQS0MRsljzySHL//cmxY10kvdX588np08nJk8kzzyQ/+cn65wSma2Sf3fYbSe5LcirJf288C1A5fjz5/OeTT30qedeuDf8LF5JDhy799ZtvJs8+m3zpS8lXv7q+eYHpGlEk/UaS/5jkvUl+luRvRihBf9x0U/L0010cHdSzzyaf+1zy6qsH/7UALmVEB7fvSxdI1y+u97UcBthlczM5c6YLpAsHfFt24UL365w5k/zO7yxlPIDSiCLpVLoVpNcW11MthwEWTpxIvvzl5Npru693b6FdjZ3//bXXJltb3a8PsAoj2m5LnEmCfjlxIvnCF9bzOk8+ufrXAaZlZJEE9MXmZreCtM7X+8pX1vd6wPiJJGDpbrqpOzO0s8W2Dm+8kdxxh8PcwPKM6EwS0BdPP73eQEq613v66fW+JjBuIglYquPHl3MX237t3PV2/Ph6XxeGbyPJby2u7CaSgKX6/Oe760HvYtuvndd77LH1vi4M20aSzSS/vbgKpd1E0kqpc6ZlNlvOwyIP4tOfTj74wbYzwHDMkhxJcnRxnTWdpm9E0sqoc6bnkUf2ftRIC+96V/Loo21ngOGYJzmb5H8urvOm0/SND7hdmd11vvP1drNpYB3uv7/1BJ37PHAfrtB2kq10P6Pm8XNqL5G0Mjt1nqhzpuLYsdYTdPoyBwzDdsRRzXOSVmoj6pypOHw4OXeuu8ts3Ye2d9t5/RtvTM6fbzcHMHxWklZKnTMd11/fXVsG0u7Xv+46kUTfeOM8NCIJWIp3v7v1BHu95z2tJ4Dddm7mOZLuCMZWhFL/ubsNWIpf/KL1BHv9/OetJ4Dd3Go/RFaSgKV47bXu2pczSa+/3m4GeDs38wyRg9vA0pw71x3gbu3cueTIkdZTwFs5kzQ0ttuApTl9uvUEnb7MAXttJ/nPEUjDIZKApTl5svUEnVOnWk8AjIHtNmBpZrNke7vtR5O8+Wbya7+W/OQn7WYAxsFKErA083ny7LNtZ/ja1wQSsBwiCViqL32pu15Y8xr1zus99dR6XxcYL9ttwNJ97WvJpz7V5nWPH1//6wLjJJKApbvppuTMmeTaa9f3mm+8kdxxR/Lqq+t7TWDcbLcBS/fqq8njj6/3NR9/XCAByzWSSPpYkt9fXIE++MpXki9+cT2vdeJE93oAyzSCjyX5WJKnkrw3yc+SPJbkvzWcB9jxxBPdgeonnljda5w4kTz55Op+fWC6RrCSdH+6QLpucb2/5TDAWzz5ZPLZz3ZnhpKD3/W2879/441kc1MgAaszgkg6mW4F6fXFtSeP/AX+v62t7lD1s88e/MNvDx3qfp077rDFBqzWSO5u+1i6FaSTsdUG/Xb8ePLYY8mnP733ydwXLuwNqLd+/eab3S3+Tz2VfPWr65sXmK6RRBIwNB/8YPLoo8l99yXHjiU33vj2f+bcue7Dak+dSp55xpO0gfUSSUAvHD6cXHdd8p73JD//efL668n5862nAqZMJAEAFEZwcBsAYPlG8JwkAODqbCSZJZkn2W48S/+IJACYpI0km0mOJDmbZCtCaS/bbQAwSbN0gXR0cZ01naaPrCStjCVMAPpsnm4FKYvrvOEs/eTutpWwhAnAEHhDfzlWklZi9xLmztf+8AHQN9vx8+nSRNJKWMIEgKGz3bYyljABYMhEEgBAwSMAAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKI/vstt9M8kCS55L818azAABDNqKPJfnNJH+U5NokbyT5vQglAOBqjWi77YF0gfTLi+sDbccBAAZtRJH0XLoVpD9fXJ9rOw4AMGgj2m5LnEkCAJZlZJEEALAcI9puAwBYHpEEAFAQSQAABZEEAFAY2RO3AYBL20gySzJPst14lv4TSQAwCRtJNpMcSXI2yVaE0uXZbgOASZilC6Sji+us6TRDYCUJACZhnm4FKYvrvOEsw+BhkgAwGc4k7YdIWit/OAFgKGy3rY0DcwAwJA5ur40DcwAwJFaS1saBOQAYEmeS1sqZJAAYCpEEAFBwJgkAoCCSAAAKIgkAoCCSAAAKIgkAoCCSAAAKE3iY5MeTPJjk20n+tPEsAMBQjPw5SR9P8nSSX0ryF0k+F6EEAFyJkW+3PZgukN67uD7YdBoAYDhGHknfTreC9LPF9dtNpwEAhmPk222JM0kAwNWYQCQBAOzfyLfbAACujkgCACiIJACAgkgCACiIJACAwgQ+lgQApmAjySzJPMl241nGQSQBwOBtJNlMciTJ2SRbEUoHZ7sNAAZvli6Qji6us6bTjIWVJAAYvHm6FaQsrvOGs4yHJ24DwCg4k7RsIqlX/AEHgL6w3dYbDt0BQJ84uN0bDt0BQJ9YSeoNh+4AoE+cSeoVZ5IAoC9EEgBAwZkkAICCSAIAKIgkAICCSAIAKIgkAIDCxJ+T9IkkDyX5VpJvNJ0EAOiXCT8C4BNJ/jjJX07yf5L8boQSALBjwtttD6ULpGsX14fajgMA9MqEI+lb6VaQ3lhcv9V2HACgVya83ZY4kwQAXMrEIwkAoDbh7TYAgEsTSQAABZEEAFAQSQAAhYk/cRsA+mwjySzJPMl241mmRyQBQC9tJNlMciTJ2SRbEUrrZbsNAHppli6Qji6us6bTTJGVJADopXm6FaQsrvOGs0yTh0kCQG85k9SSSAIAKNhuGzTvMABgVUTSYLnrAQBWyd1tg+WuBwBYJStJg+WuBwBYJQe3B82ZJABYFZEEAFBwJgkAoCCSAAAKIgkAoODutiv2ySQPJ/lmkq83ngUAWDUHt6/IJ5N8OclfSvJ/k3w2QgkAxs122xV5OF0g/dLi+nDbcQCAlRNJV+Sb6VaQ/mJx/WbbcQCAlbPddsWcSQKAKRFJAAAF220AAAWRBABQEEkAAAWRBABQ8MRtAFiZjSSzJPMk241nYb9EEgCsxEaSzSRHkpxNshWhNCy22wBgJWbpAuno4jprOg37ZyUJAFZinm4FKYvrvOEsXA0PkwSAlXEmachEEgBAwXbbpHmHAwCXIpImy10XAHA57m6bLHddAMDlWEmaLHddAMDlOLg9ac4kAcCliKSV+USSR5LckOSnSf4kyTcazgMA7IdIWrr3JvlPST5Z/HdfT/KZJD9b4zwAwNUQSUv3X9IF0oV0//fu2Pn660n+WoO5AID9cHfbUn0idSAlF3v0k0k+vt6xAIB9E0lL9cji+tZAylv+80fXMAsAcBAiaaluWPI/BwC0IpKW6qdL/ucAgFZE0lL9yeJ6qbPwO//5M2uYBQA4CJG0VN9Id/daddPg7rvb/nS9YwEA+yaSlu4zuRhKu+0E0mfWPhEAsH+ek7QyH093F9vOE7efiRUkABgOkQQAb+OzLUmuaT0AAPTLRpLNJEeSnE2yFaE0Tc4kAcAes3SBdHRxnTWdhnasJAHAHvN0K0hZXOcNZ6ElZ5IA4G2cSUIkAQCUbLexJN51ATAuIoklcCcIAOPj7jaWwJ0gAIyPlSSWwJ0gAIyPg9ssiTNJAIyLSAIAKDiTBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAUfSwLAwHniP6shkgAYsI0km+k+XPtskq0IJZbFdhsAAzZLF0hHF9dZ02kYFytJAAzYPN0KUhbXecNZGBsfcAvAwDmTxGqIJHrONz8A2rDdRo85kAlAOw5u02MOZALQjpUkesyBTADacSaJnnMmCYA2RBIAQMGZJACAgkgCACiIJACAgkgCACiIJACAguckAbAGHufB8IgkAFbMRwwxTLbbAFgxHzHEMFlJAmDFfMQQw+SJ2wCsgTNJDI9IYuJ84wagZruNCXOYFIBLc3CbCXOYFIBLs5LEhDlMCsClOZPExDmTBEBNJAEAFJxJAgAoiCQAgIJIAgAouLsNYPLcwAAVkQQwaR6qCpdiuw1g0jxUFS7FShLApHmoKlyK5yTBWjjzQZ/58wkVkQQr58wHwBA5kwQr58wHwBA5kwQr58wHwBDZboO1cOYDYGhEEgBAwXYbQC9ZfYTWRBJA77gjEvrA3W0AveOOSOgDK0kwSrZqhs0dkdAHDm7D6NiqGQehC61ZSYLR2b1Vs/O1H7LDsx3/3qAtkQSjY6tmNazswNTYboNR8gN9uWxhwhRZSYJRWvdWzdijzBYmTJFIAg6o1SrLOsPMFiZMkUgCDqjFKsu6w2x78RpjXi0D3kokAQfUYpWlRZi52wymRiQBB9RilcX2F7B67m4DBmrsh8WB1kQSAEDBB9wCABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBAQSQBABREEgBA4f8BSwHCDiCB6XQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_orbits(data_euler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Leapfrog integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the Earth is not on a circular orbit. Its orbital distance is increasing and because of that the Earth could not achieve two full orbital cycles.\n",
    "\n",
    "The problem is the simple integration scheme used here.\\\n",
    "Every integration scheme has numerical error. Some have larger ones than others. Euler's 1st order method is simply not suited for orbital integration.\n",
    "\n",
    "But there is a way out: [Symplectic integration](https://en.wikipedia.org/wiki/Symplectic_integrator)\\\n",
    "Symplectic integrators conserve the energy of the system.\n",
    "\n",
    "One of these methods is the Leapfrog method.\\\n",
    "Leapfrogging means the velocity and the position are not updated synchronous, but in between each other. They are leapfroggin each other.\n",
    "\n",
    "In our case we first update the velocities by using half of the step size. Then we update the positions for a full time step. And finally, we update the velocities for another semi time step. We can do this by using the `fstep` keyword argument for the instructions, which is the fraction of the time for which this instruction is applied.\n",
    "\n",
    "But there is one caveat: usually the fields that are integrated in an instructions set are only updated after all instructions have been executed, such that the second instruction still sees the original value of the first variable and not already the new value.\n",
    "\n",
    "But for leapfrog integration, the instruction for the positions already need to know the new values of the velocities. And the second set of velocity instructions already need to know the new positions. We therefore have to add `update` instructions in between. \n",
    "\n",
    "We do not have to add update instructions after the second set of velocity operations, because at the end of an instruction set all fields contained in the instruction set will be updated.\n",
    "\n",
    "The new leapfrog instruction set now looks like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "instructions_leapfrog = [Instruction(schemes.expl_1_euler, sim.Sun.v,   fstep=0.5),\n",
    "                         Instruction(schemes.expl_1_euler, sim.Earth.v, fstep=0.5),\n",
    "                         Instruction(schemes.update,       sim.Sun.v             ),\n",
    "                         Instruction(schemes.update,       sim.Earth.v           ),\n",
    "                         Instruction(schemes.expl_1_euler, sim.Sun.r,   fstep=1.0),\n",
    "                         Instruction(schemes.expl_1_euler, sim.Earth.r, fstep=1.0),\n",
    "                         Instruction(schemes.update,       sim.Sun.r             ),\n",
    "                         Instruction(schemes.update,       sim.Earth.r           ),\n",
    "                         Instruction(schemes.expl_1_euler, sim.Sun.v,   fstep=0.5),\n",
    "                         Instruction(schemes.expl_1_euler, sim.Earth.v, fstep=0.5),\n",
    "                        ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can assign this new instruction set to our frame. To rerun the simulation we have to reset the initial conditions, that we've saved earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.integrator.instructions = instructions_leapfrog\n",
    "sim.integrator.description = \"Leapfrog integrator\"\n",
    "sim.Earth.r = r_Earth_ini\n",
    "sim.Earth.v = v_Earth_ini\n",
    "sim.Sun.r   = r_Sun_ini\n",
    "sim.Sun.v   = v_Sun_ini\n",
    "sim.t       = 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally we have to reset the buffer of the namespace writer. Otherwise, we would simply add snapshots to the old dataset. Furthermore, we descrease the `verbosity` of the writer to prevent it from writing information on screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.writer.reset()\n",
    "sim.writer.verbosity = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Execution time: \u001b[94m0:00:00\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sim.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_leapfrog = sim.writer.read.all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAJJCAYAAAC+gKM0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAABcSAAAXEgFnn9JSAAAZWUlEQVR4nO3dS6yf913n8U+rqt0MJJAQbDe05sDJrc4Nt3YuQNPZTPGSFVIqYYQioRSIuImL2iBwNcNouM7QCMlUdaWMNCsWs2jDqi2QNElrenPbJFYPLCLblAQSYNMuYhbP/8jH7vfY5397rq+XVD06aY6f3yLJeZ/f7/s8/zcluRgAAC7z5q4XAADQRyIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACiIJAKAgkgAACm/pegHA2N2a5EeSfDPJi3N+72aSjSRbSc6ueF0AVyeSgDW6NcnPJzmQ5FySj2XvobSZ5PiO7z2V+UPpllyKrJfm/F5g6kQScBXL7uT8SJrIObjj671G0sYV37sx5xpuSfJzuRRZH49QAuYhkoBdrGIn55uz783s+s05vnfriu/dmvPeVWTNE0mO+mDqRBKM3m25NBP0whzft+xOTtLsGn0si80knU0TZouGyjKRtYpABIZOJMGo3ZbkkVz6YX8yew+lZXdytr2Y+Qe2t53N4nHyUpojtkVmklYRiMDQiSQYte2ZoHfu+HqvkbTsTk4fvJTF5pBWFYjAkIkkGLVlZoKS5XZyhmxVgWiuCYbsTUkudr0IYJ0WnUliOeaaYOjsJEHvLbsb8ULEURfMNcHQiSToNbsRw2WuCYZOJEGv2Y0YrjEMvsO0iSToNbsRw7bo4LuBb+gDg9vQe35gTosjVugLO0nQmtvT/AA8m+Qbc3zfVB/DnypHrNAXIglacXuSX8il3YG/yHyhxHQ4YoW+EEnQis1c/ubrzYgkaga+oS9EErTibC7fHfCDj6txxAp9YHAbWrPoTBIAXRBJAKPiaUhYFcdtAKPh9QGwSm/uegEArMrO1wccmH0NLMpOEszFUQZ95vUBsEpmkmDPHGUwBEIeVsVOEuyZNyEzBF4fAKsikmDPHGUATInjNpiLowzGyD/XUBFJAJNm1g524xUAAJPmtQGwGzNJTNgdufQxIV/veC3QFbN2sBvHbUzUHUkeTbI/yfkkT0QoMV1mkqBiJ4mJ2kwTSAd3fC2SmCqvDYCKSGKizqbZQcrs6gcEAJdz3MaEmUkCYHciCQCg4LgNgAUZ+GbcRBIAC/ASSsbPyyQBWICXUDJ+dpIAWICXUDJ+BrcZAXMR0A3/7jFuIomBMxcBwHqYSWLgzEUAsB5mkhg4cxEArIfjNkbAXAQAqyeSAGiJX2gYFsdtALTAQxYMj8FtAFrgIQuGx04SAC3wkAXDYyYJgJaYSWJYRBI99a4ktyR5KcnXOl4LAFPkuI0eeleSX0qyP8n5JP8nQgmAthncpoduSRNI75xdb+l2OQBMkp0keuilNDtImV1f6nAtAEyVmSR6ykwSAN0SSQAABTNJAAAFM0kA9JT3KtEtkQRAD/msN7rnuA2AHvJZb3TPThIAPeSz3uiep9sA6CkzSXRLJNES/7EDYFgct9ECA5gADI/BbVpgABOA4bGTRAsMYAIwPGaSaImZJACGRSQBABQctwEwcHaqWQ+RBMCAeXqW9fF0GwAD5ulZ1sdOEgAD5ulZ1sfgNgADZyaJ9RBJLOlQkluTvJjkTMdrAYDVcdzGEg4leSzJ/iTnk/xZhBIAY2FwmyXcmiaQ3jm73trtcgBghewksYQX0+wgZXZ9scO1AMBqmUliSWaSABgnkQQAUDCTBABQEEkAAAWD2wBMhJdOMh+RBMAE+CBc5ue4DYAJ8EG4zM9OEgAT4INwmZ9XAAAwEWaSmI9IAgAoOG7jKvzWBcB0iSR24UkQAKbN023swpMgAEybnSR24UkQAKbN4DZXYSYJgOkSSQAABTNJAAAFM0kAkMSIAVcSSQDgtScUHLcBgNeeULCTBABee0LB020AkMRMElcSSQAABcdtk3ZnktuSvJDkqx2vBQD6RSRN1p1JfjXJ/iTnk/xxhBIAXOLptsm6LU0gvWN2va3b5QBAz9hJmqwX0uwgZXZ9ocO1AED/GNyeNDNJALAbkQQAUDCTBABQEEkAAAWRBABQ8HQbAOyJjy2ZGpEEANe0meR4kgNpPgD3VITS+DluA4Br2kgTSAdn141OV0M77CQBwDVtpdlByuy61eFaaIv3JAHAnphJmhqRBABQcNw2en7zAYBFiKRR8zQGACzK022j5mkMAFiUnaRR8zQGACzK4PbomUkCgEWIJACAgpkkAICCSAIAKBjcBoClmP0cK5EEAAvzProxc9wGAAvzProxs5MEAAvzProx8woAAFiKmaSxEkkAAAXHbaNyV5Lbk3wjyVc6XgsADJtIGo27kvx6kn1JLiT5wwglAFicp9tG4/Y0gfSO2fX2bpcDAANnJ2k0vpFmBymz6zc6XAsADJ/B7VExkwQAqyKSAAAKZpIAAAoiCQCgIJIAAAoiCQCg4BUAALByPs9tDEQSAKzUZpLjSQ4kOZfkVITSMDluA4CV2kgTSAdn141OV8Pi7CQBwEptpdlByuy61eFaWIaXSQLAyplJGgORBABQMJMEAFAwkzRotnMBYF1E0mB5xBQA1slx22B5xBQA1slO0mB5xBQA1snTbYNmJgkA1kUkAQAUzCQBABREEgBAweA2AKyV+dGhEkkAsDbeaTdkjtsAYG28027I7CQBwNp4p92QeQUAAKyVmaShEkkAAAXHbYNzd5I7knw9yZc7XgsAjJdIGpS7k/xmkn1JLiT5nxFKALAenm4blDvSBNIPza53dLscABgxO0mD8vU0O0iZXb/e4VoAYNwMbg+OmSQAaINIAgAomEkCACiIJACAgkgCACh4ug3o3L59yY03Jm99a/Kd7ySvvJJcuHDt7wNYJ5EEtG5jI3n44eTo0eTw4SaSrnThQnL6dPLcc8mTTyb/8A/trxOYNk+3Aa05diz54AeT978/efOOw/6LF5M3vWn3r994I3nqqeSjH00++cn21gur5YNuh0YkAWt3883JyZNNHC3rqaeSRx5JXn55+T8L2rOZ5HiSA0nOJTkVodR/BreBtTp+PDlzpgmki0v+SnbxYvPnnDmT/OzPrmR50JKNNIF0cHbd6HQ17I1IAtbmxInk4x9Prruu+XrnEdoitr//uuuSU6eaPx+GYSvNDtI/zq5bna6GvXHcBqzFiRPJhz7Uzn0ef3z994HlmUkaGpEErNzx480OUpv3+8Qn2rsfMA0iaVD8FkL/3XxzMzO0fcTWhtdfTw4dMswNrJaZpMHYfjLiZ2bXzS4XA7s6ebLdQEqa+5082e49gfETSYPhyQj679ix1TzFNq/tp96OHWv3vsC4iaTB8GQE/ffBDzbXZZ9im9f2/R59tN37AuNmJmlQzCTRXxsbydmzl79Ju21vvJH86I/6CBNgNewkDcrZJH8dgUQfPfxwt4GUNPf/wAe6XQMwHiIJWImjR7teQePIka5XAIyF4zZgJc6fT/bt63oVzToOHOh6FXClh5I8nOSmJN9K8n+TfKbD9bAXIglY2r59TZxcvNj+0PZO2/ffvz+5cKG7dcAl1yf5qyTvK/6/Tyf56SSvtbge5uG4DVjajTc21y4Daef9b7ih23XAJduBdOV+xMXZX/+r1lfE3okkYGlvfWvXK7jc297W9QogaY7YtgPpyt8gtg9y3pfkve0uiz0TScDSvvOdrldwuW9/u+sVQNLMICXfHUi54q97JLOvRBKwtFdeaa5tv2n7Stv3f/XVbtcBjZtW/PfRNpEELO3CheZ/fZhJOn/e0DZ98a0V/320TSQBK3H6dNcraPRlHdA85p/s/hD59l9/soW1sAiRNCj3pDnjvqfbZUDhuee6XkHj+ee7XgFs+0yax/yrt+1sD3N/Osln210We+Y9SYNxT5LfSrIvyYUkf5DkSx2uBy7ns9ugcn28J2m47CQNxrvSBNI7Ztd3dbscuMLWVvLUU92u4VOfEkj0zWtJ/mua1wH8ZZL/P7s+NPvrr3WzLPbkLV0vgL36WpodpMyuX+twLVD76EeTY8faf/P29v2eeKK9e8J8PhvHasPjuG1Q7kmzg/S1OGqjrz71qeT97+/mvseOtX9fYLxEErBSN9+cnDmTXHdde/d8/fXk0KHk5ZfbuycwfmaSgJV6+eXkscfavedjjwkkYPVEErByn/hE8pGPtHOvEyea+wGsmkgC1uLDH24CZp1OnEgef3y99wCmy0wSsFbHjyd/+qfNjNKyT71tf//rrzdHbHaQgHUSScDa3XxzcvLkap56e+qp5JFHzCAxRJtJNpJsJTnb8VrYC5EEtObYseTRR5Of+qnL38x95Q7TlV+/8UbziP8TTySf/GR764XV2UxyPMmBJOeSnIpQ6j+RBLTuh384+cAHkiNHksOHk/37v/vvOX+++bDa559PnnzSm7QZuv+W5GeSHEzyj0n+X5K/7nA97IVIAjq3b19yww3J296WfPvbyauvJhcuXPv7YDjsJA2RSAKAVphJGhqRBABQ8J4kAIDCW7peAPOyXQsAbRBJg2LwDwDa4rhtUDbSBNLB2XWj09UAwJjZSRqUrTQ7SJldtzpcCwCMm6fbBsdMEgC0QSQBABTMJAEAFEQSAEBBJAEAFEQSAEDBKwAAYG3el+TBJE8n+XTHa2Fenm4DgLV4X5KPJ/kvSf4jyc9FKA2L4zYAWIsH0wTS986uD3a7HOYmkgBgLZ5Os4P0b7Pr090uh7k5bgOAtTGTNGQiabDuTXIoyZkkX+x4LQAwPp5uG6R7k/xOkn1JLiT57xFKALBaZpIG6VCaQPqh2fVQt8sBgBGykzRIZ9LsIGV2PdPhWgBgnMwkDZaZJABYJ5EEAFAwkwQAUBBJAAAFkQQAUBBJAAAFrwAAgJXaTLKRZCvJ2Y7XwjJEEgCszGaS40kOJDmX5FSE0nA5bgOAldlIE0gHZ9eNTlfDcuwkAcDKbKXZQcrsutXhWliWl0kCwEqZSRoLkQQAUHDcNip+ewGAVRFJo+GJCgBYJU+3jYYnKgBglewkjYYnKgBglQxuj4qZJABYFZEEAFAwkwQAUBBJAAAFkQQAUBBJAAAFrwAAgIU9lOSBJM8k+UynK2H1PN0GAAt5KMnHknxPkn9P8vMRSuPiuA0AFvJAmkC6bnZ9oNvlsHIiCQAW8kyaHaTXZ9dnul0OK+e4bfR+LMmdSb6a5O87XgvA2DwUM0njJZJG7ceSfDjJviQXkpyIUAKAvXHcNmp3pgmkt8+ud3a7HAAYEK8AGLWvptlByuz61Q7XAgDD4rht9MwkAcAiRBIAQMFMEgBAQSQBABREEgBAwdNtAHBNm0k2kmwlOdvxWmiLSAKAq9pMcjzJgSTnkpyKUJoGx20AcFUbaQLp4Oy60elqaI+dJAC4qq00O0iZXbc6XAtt8p4kALgmM0lTJJIAAAqO2ybNb0YAsBuRNFme1gCAq/F022R5WgMArsZO0mR5WgMArsbg9qSZSQKA3YgkAICCmSQAgIJIAgAoiCQAgIKn2wAg701yf5LPJflsx2uhLwxuAzBx701yMsn3Jvm3JI9EKJE4bgNg8u5PE0jXz673d7oa+kMkATBxn0uzg/Ta7Pq5TldDfzhuYxeHk9yV5CtJTne8FoB1M5PEdxNJFA4n+d0kP5jkn5L8XoQSAFPjuI3CXWkC6e2z613dLgcAOuAVABS+kmYHKbPrVzpcCwB0w3EbuzCTBMC0iSQAgIKZJACAgkgCACiIJACAgqfbAJiAzSQbSbaSnO14LQyFSAJg5DaTHE9yIMm5JKcilNgLx20AjNxGmkA6OLtudLoahsNOEgAjt5VmBymz61aHa2FIvCcJgAkwk8T8RBJL8B8dAMbLcRsLMggJwLgZ3GZBBiEBGDc7SSzIICQA42YmiSWYSQJgvEQSAEDBTBIAQEEkAQAUDG4DMGA/meS+JM8m+ZuO18LYmEkCYKB+MslfJLkuyetJfiFCiVVy3AbAQN2XJpC+f3a9r9vlMDoiCYCBejbNDtK/zK7PdrscRsdxGy14d5K7k3w5yRc6XgswLmaSWB+RxJq9O8nvJ7kpybeSPB6hBMAQOG5jze5OE0hvn13v7nY5ALBHXgHAmn05zQ5SZtcvd7gWANg7x220wEwSAMMjkgAACmaSAAAKZpIA6KHNJBtJtpKc7XgtTJVIAqBnNpMcT3IgybkkpyKU6ILjNgB6ZiNNIB2cXTc6XQ3TZScJgJ7ZSrODlNl1q8O1MGWebgOgh8wk0T2RRA/5jyMA3XPcRs8Y2ASgHwxu0zMGNgHoBztJ9IyBTQD6wUwSPWQmCYDuiSQAgIKZJACAgpkkAFrwE0mOJnkuyd92vBbYG5EEwJr9RJInklyf5LUkj0YoMQSO2wBYs6NpAumG2fVol4uBPbOTxMC9J8k9Sb6U5POdrgTYzXNpdpAyuz7X2UpgHp5uY8Dek+QjSW5K8q0kH4pQgr4yk8Tw2EliwO5JE0gHdnwtkqCf/jbiiKERSQzYl9LsIGV2/VJnKwFgfBy3MXBmkgBYD5EEAFBw3AbAAnzGIuMnkgCY02aS42kemjiX5FSEEmPkZZIAzGkjTSAdnF03Ol0NrIudJADmtJVmBymz61aHa4H1MbjNBJmlgOX594jxE0lMjFkKAPbGTBITY5YCgL0xk8TEmKUAYG8ctzFBZikAuDaRBDBpP57kSJLnk/xdx2uBfnHcBjBZP57kz5Ncn+S1JL8YoQSXGNwGmKwjaQLpxtn1SJeLgd6xkwR7ciTJvUm+mOZYAsbg+TQ7SJld/bMNO5lJgms6kuR/JPmBJP+c5LfjhwnjYSYJdmMnCa7p3jSBtH/H1yKJsfi7iCOoiSS4pi+m2UHK7PrFDtcCQFsct8GemEliCLwDDFZJJAGMgs8lhFXzCgCAUfC5hLBqZpJg7RyB0AafSwir5rgN1soRCG0S5LBKdpJgrXYegWx/7YcX63I2/vmC1RFJsFaOQACGynEbrJ0jEPbqwSTvSfL5JE93vBZAJEFvHc2ldzM91/FaWL8Hk/zvJN+X5F+T/HKEEnTLcRv00tEkf5BLnxf3WxFKY/eeNIH0Azu+FknQJe9Jgl7a/ry4fbPrvd0uhxZ8Ps0O0j/Prp/vdjmAnSToJ58XNz1PpzliM5MEfWEmCXrLTNJwGdaHMRBJMFp+UHfDC0RhLBy3wSj5Qd0dLxCFsRBJMEp+UC/ugSTvTvKFJM8s8P1eIApjIZJglJb9QX1fksNJTid5doXr6rsHkvxJku9P8i9JfiXzh9LZNDt3jjph6EQSjNIyP6jvS/K/ktyY5JUkv5HphNK70wTSTTu+XmQ3yWeowRiIJBitRX9QH04TSPt2fL1IJC0zOH5/Lu1kfW6Bey96/y+k2UHK7PqFBe8NjIFIAq5wOs0OUmbX0wv8GcsMjt+f5I+S3JDk1SS/lvlDadH7P5PmiG2ZmSRgLEQScIVn0xyxLTOTtMzg+OE0gfSDO76eN5KWuf8zEUdAIpKA0rNZbg5pmcHx02l2kDK7LrKT5QkzYHleJgmsyRBnkgAuEUkAAIU3d70AAIA+EkkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAAWRBABQEEkAAIX/BDCi1OpRd4ieAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_orbits(data_leapfrog)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}