Modeling Basics

Here are a few examples of how to get started modeling transit light curves.

batman for simulating transits

The batman package is a Python toolkit for generating model transit light curves for planets of given parameters. Laura Kreidberg wrote it to make transit-modeling with Python easier, and she did a great job documenting how it works.

We wrote some quick wrapper so you can immediately start playing, by looking at transit models by themselves:

import henrietta as hsl

We can make a plot of a model transit light curve for some set of parameters.

hsl.example_transit_model(radius=0.1, a=4, b=0.9,  period=0.47, t0=0.05);

If I want to compare two or more different models, I can keep track of the plotting axes in which the first is plotted, and keep adding new curves to the same plot, on the same scale. For example:

ax = hsl.example_transit_model(radius=0.1, a=4, b=0.9,  period=0.47, t0=0.05);
ax = hsl.example_transit_model(radius=0.1, a=4, ax=ax);

If we want to access the modeling tool directly, we can use the BATMAN function to make an array of model flux values. The inputs to this code are an array of times and the transit model parameters.

import numpy as np, matplotlib.pyplot as plt

# make some times that span 0 to 5 days, in 2-minute increments
times = np.arange(0, 1, 2.0/60.0/24.0)

# use the BATMAN function to make a transit model
modelflux = hsl.BATMAN(times, period=4.2, t0=0.3, radius=0.11, a=9.0, b=0.2)

# plot the times and model fluxes
plt.scatter(times, modelflux, s=3)
plt.xlabel('Time (days)')
plt.ylabel('Relative Flux');

All of the other visualization tools described on this page are calling this BATMAN function to make their models.

making a simulated transit dataset

Sometimes it can be really helpful to be able to make a simulated noisy dataset, starting from some model parameters. The simulate_transit_data function can do that for you!

lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)
ax = lc.scatter()

visually comparing a model to data

We want to be able to compare our transit model to our data. Here’s a little wrapper called plot_with_transit_model, which can plot a light curve, and then overplot a transit model on top. As inputs, this function takes a light curve, and the transit parameters period, t0, radius, a, b as keyword arguments. The top panel shows the light curve data; the bottom panel shows the “residuals”, which are defined as the data minus the model.

lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)
hsl.plot_with_transit_model(lc, period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)

If we plot the wrong model for our data, our residuals show big bumps in them!

lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)
hsl.plot_with_transit_model(lc, period=1.23, t0=0.1, radius=0.1, a=10, b=0.0)

If we feed this function a folded light curve, it will plot everything in terms of phase instead of time. The x-axis will span only one orbital period, with data repeating back on it many times.

lc = hsl.simulate_transit_data(period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)
folded = lc.fold(period=1.23, phase=0.1/1.23)
hsl.plot_with_transit_model(folded, period=1.23, t0=0.1, radius=0.1, a=10, b=0.5)

keeping times straight

One tricky thing might be how to match up a mid-transit time you take off the NASA Exoplanet Archive (in BJD) with the times in your light curves. Remember that most Kepler light curves store their time as BKJD (Barycentric Kepler Julian Date), which is offset from actual BJD (Barycentric Julian Date) by 2454833 days. We’ve added a couple of small tools to help with this conversion:

# print the original t0, in BJD
t0 = 2455678.1234
print(t0, "is the BJD.")

# convert to BKJD (to match Kepler light curves)
t0_bkjd = hsl.bjd2bkjd(t0)
print(t0_bkjd, "is the BKJD.")

# convert back again
t0_bjd = hsl.bkjd2bjd(t0_bkjd)
print(t0_bjd, "is the BJD, after converting back and forth to BKJD.")