{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Catching close encounters using exceptions\n", "Sometimes one is interested in catching a close encounter between two planets. This can easily be done with REBOUND. What you do when a close encounter happens is up to you.\n", "\n", "Some integrators are better suited to simulate close encounters than others. For example, the non-symplectic integrator IAS15 has an adaptive timestep scheme that resolves close encounters very well. Integrators that use a fixed timestep like WHFast are more likely to miss close encounters.\n", "\n", "Let's start by setting up a two-planet system that will go unstable on a short timescale:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import rebound\n", "import numpy as np\n", "def setupSimulation():\n", " sim = rebound.Simulation()\n", " sim.integrator = \"ias15\" # IAS15 is the default integrator, so we don't need this line\n", " sim.add(m=1.)\n", " sim.add(m=1e-3,a=1.)\n", " sim.add(m=5e-3,a=1.25)\n", " sim.move_to_com()\n", " return sim" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's integrate this system for 100 orbital periods." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "sim = setupSimulation()\n", "sim.integrate(100.*2.*np.pi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rebound exits the integration routine normally. We can now explore the final particle orbits:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n" ] } ], "source": [ "for o in sim.orbits():\n", " print(o)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the orbits of both planets changed significantly and we can already speculate that there was a close encounter.\n", "\n", "Let's redo the simulation, but this time set the `sim.exit_min_distance` flag for the simulation. If this flag is set, then REBOUND calculates the minimum distance between all particle pairs each timestep. If the distance is less than `sim.exit_min_distance`, then the integration is stopped and an exception thrown. Here, we'll use the [Hill radius](https://en.wikipedia.org/wiki/Hill_sphere) as the criteria for a close encounter. It is given by $r_{\\rm Hill} \\approx a \\sqrt{\\frac{m}{3M}}$, which is approximately 0.15 AU in our case. \n", "\n", "This setup allows us to catch the exception and deal with it in a customized way. As a first example, let's catch the exception with a `try`-`except` block, and simply print out the error message. Additionally, let's store the particles' separations while we're integrating:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Two particles had a close encounter (d" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "fig = plt.figure(figsize=(10,5))\n", "ax = plt.subplot(111)\n", "ax.set_xlabel(\"time [orbits]\")\n", "ax.set_xlim([0,sim.t/(2.*np.pi)])\n", "ax.set_ylabel(\"distance\")\n", "plt.plot(times/(2.*np.pi), distances);\n", "plt.plot([0.0,12],[0.2,0.2]); # Plot our close encounter criteria;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We did indeed find the close encounter correctly. We can now search for the two particles that collided and, for this example, merge them. To do that we'll first calculate our new merged planet coordinates, then remove the two particles that collided from REBOUND and finally add the new particle." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of particles at the beginning of the simulation: 3.\n", "Two particles had a close encounter (d