Last active
July 8, 2018 22:38
-
-
Save hannorein/9241f10e27a6f69af2b087a64f791fab to your computer and use it in GitHub Desktop.
Using the REBOUND Kepler Solver
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Using the Kepler Solver in REBOUND" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"First, let's import REBOUND and a few wrapper functions from the ctypes library." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 1, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"Date on which REBOUND was compiled: Jul 8 2018 23:01:09.\n" | |
] | |
} | |
], | |
"source": [ | |
"import rebound\n", | |
"import numpy as np\n", | |
"from ctypes import c_double, byref, c_int\n", | |
"print(\"Date on which REBOUND was compiled: {}.\".format(rebound.__build__))" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Next up, we create a `Simulation` object. This is needed because the Kepler solver needs to know the gravitational constant (which is a property of the `Simulation` object and defaults to 1). " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 2, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"sim = rebound.Simulation()" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"Next up, we create the particle we want to evolver with the Kepler solver. Only the x, y, z, and vx, vy, vz properties of this particle matter." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"p = rebound.Particle(x=1.,vx=1.)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The mass of the particle is not used by the Kepler solver. It expects a mass as a separate argument when calling it (this allows you to use different coordinates systems). Here, we simply use:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 4, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"M = 1." | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The last argument that we pass to the Kepler solver is the timestep. As a test, let's evolve for one full orbit. In units where $G=1$ and the mass in the system is $M=1$, an orbit takes $2\\pi$ time units." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 5, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"dt = 2.*np.pi" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"We can now finally call the Kepler solver with this command:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 6, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"rebound.clibrebound.reb_whfast_kepler_solver(byref(sim),byref(p),c_double(M),c_int(0),c_double(dt));" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"The function will not return anything, but it will change the Particle that we pass it. To see the difference, let us print out the coordinates of `p`:" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 7, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"<rebound.Particle object, m=0.0 x=1.0000000000000044 y=0.0 z=0.0 vx=0.9999999999999991 vy=0.0 vz=0.0>\n" | |
] | |
} | |
], | |
"source": [ | |
"print(p)" | |
] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"As you can see we almost get exactly back what we put in, except a small difference at the 15th or 16th digit. This is because of floating point precision, which limits us to a relative precision of $\\approx 10^{-16}$. " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"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.6.5" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 2 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment