Calculating transit timing variations#
In this example we will be using heyoka.py’s event detection system to find the transit times in a two-planet system. This example is directly inspired by a notebook in the documentation of the REBOUND software package.
We begin as usual with importing heyoka.py and a couple of extra useful modules:
import pykep as pk
import heyoka as hy
import numpy as np
Next, we define the initial conditions for the two planets in the system. The inner planet is on a low-eccentricity orbit, while the outer planet is on a circular orbit. Both orbits have zero inclination, so that the motion is constrained to the
# Initial cartesian conditions for the two planets.
# We use pykep to convert from Keplerian elements.
r1, v1 = pk.par2ic([1.0, 0.1, 0.0, 0.0, 0.25, 0.0])
r2, v2 = pk.par2ic([1.757, 0.0, 0.0, 0.0, 0.0, 0.0])
# Masses of the 3 bodies in the system.
masses = np.array([1, 1e-5, 1e-5])
We can now proceed to the definition of the system of differential equations via the model.nbody()
helper:
sys = hy.model.nbody(3, masses=masses)
The next step is the definition of the event for the detection of the planetary transit. In this setup, we assume that the observer of the system is in the direction of the positive
Because
Finally, we will log the trigger time of each event in a global list for further analysis.
Let us see the definition of the event:
# Global list of transit times.
tt_list = []
# The callback that will be invoked
# when a transit is detected.
def cb(ta, time, d_sgn):
# Append the transit time to the
# global list.
tt_list.append(time)
# Define the event object.
tt_event = hy.nt_event(
hy.expression("y_1"), # The event equation.
cb, # The callback.
direction=hy.event_direction.positive, # The event direction.
)
We can now proceed to the creation of the integrator object. After creation, we will adjust the system state so that the centre of mass of the system is in the origin with zero velocity:
ta = hy.taylor_adaptive(
sys, [0.0] * 6 + list(r1) + list(v1) + list(r2) + list(v2), nt_events=[tt_event]
)
# Adjust the COM.
com_x = np.sum(ta.state[0::6] * masses) / np.sum(masses)
com_y = np.sum(ta.state[1::6] * masses) / np.sum(masses)
com_z = np.sum(ta.state[2::6] * masses) / np.sum(masses)
com_vx = np.sum(ta.state[3::6] * masses) / np.sum(masses)
com_vy = np.sum(ta.state[4::6] * masses) / np.sum(masses)
com_vz = np.sum(ta.state[5::6] * masses) / np.sum(masses)
ta.state[0::6] -= com_x
ta.state[1::6] -= com_y
ta.state[2::6] -= com_z
ta.state[3::6] -= com_vx
ta.state[4::6] -= com_vy
ta.state[5::6] -= com_vz
Let us now propagate the system state for a few time units. Whenever a transit is detected, the callback will be invoked and the global list of transit times will be updated:
ta.propagate_until(1100.0)
(<taylor_outcome.time_limit: -4294967299>,
0.2296119265519109,
0.4928525973772395,
3344,
None,
None)
We can now proceed to the analysis of the transit times. First we will run a linear fit to remove the linear trend from the transit times:
N = len(tt_list)
A = np.vstack([np.ones(N), range(N)]).T
c, m = np.linalg.lstsq(A, tt_list, rcond=-1)[0]
We are now ready to visualize the transit time variations (TTVs):
%matplotlib inline
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(10, 5))
ax = plt.subplot(111)
ax.set_xlim([0, N])
ax.set_ylim([-0.8, 1.0])
ax.set_xlabel("Transit number")
ax.set_ylabel("TTV [hours]")
plt.scatter(
range(N),
(np.array(tt_list) - m * np.array(range(N)) - c) * (24.0 * 365.0 / 2.0 / np.pi),
);

The transit times deviate from a purely linear behaviour because of the perturbations induced onto the inner planet by the outer planet.