ktransit
A simple exoplanet transit modeling tool in python
This python package contains routines to create and/or fit a transiting planet model. The underlying model is a Fortran implementation of the [Mandel & Agol (2002)] (http://iopscience.iop.org/1538-4357/580/2/L171/fulltext/16756.text.html) limb darkened transit model. The code will calculate a full orbital model and eccentricity can be allowed to vary.
With version 0.1, radial velocity data can now be calcaulted via the model and included in the fit
Installation
Install via pip
pip install ktransit
or via the git repository
git clone https://github.com/mrtommyb/ktransit.git
cd ktransit
python setup.py install
Making a model
The basic module, ktransit, will create a transit model. The code uses sensible defaults and only changes parameters explicitly stated. For example, to create a simple model just run
import ktransit
import matplotlib.pyplot as plt
M = ktransit.LCModel()
M.add_star()
M.add_planet()
M.add_data()
tmod = M.transitmodel
plt.plot(M.time,tmod)of if you'd just prefer to see an Earth-like transit model then run
time,earthlc = ktransit.give_me_earth()
plt.plot(time,earthlc)Now, the code allows the user to input their own parameters.
import ktransit
import matplotlib.pyplot as plt
import numpy
M = ktransit.LCModel()
M.add_star(
rho=1.5, # mean stellar density in cgs units
ld1=0.2, # ld1--4 are limb darkening coefficients
ld2=0.4, # if only ld1 and ld2 are non-zero then a quadratic limb darkening law is used
ld3=0.0, # if all four parameters are non-zero we use non-linear flavour limb darkening
ld4=0.0,
dil=0.0, # a dilution factor: 0.0 -> transit not diluted, 0.5 -> transit 50% diluted
zpt=0.0 # a photometric zeropoint, incase the normalisation was wonky
)
M.add_planet(
T0=1.0, # a transit mid-time
period=1.0, # an orbital period in days
impact=0.1, # an impact parameter
rprs=0.1, # planet stellar radius ratio
ecosw=0.0, # eccentricity vector
esinw=0.0,
occ=0.0) # a secondary eclipse depth in ppm
M.add_planet() # you can add as many planets as you like (up to 10)
M.add_data(
time=numpy.arange(0,10,0.0188), # timestamps to evaluate the model on
itime=numpy.zeros_like(numpy.arange(0,10,0.0188))+0.0188 ) # integration time of each timestamp
tmod = M.transitmodel # the out of transit data will be 0.0 unless you specify zpt
plt.plot(M.time,tmod)Fitting
This is the fun part!! Fit a transit model to data. I use the scipy Levenberg–Marquardt least-squares minimization You need to give the model an initial guess and tell it what parameters you want to allow to vary, everything else is kept fixed.
The code assumes that the observed data has been normalised and detrended. The default zeropoint is 0.0 so the out of transit data should have a mean of 0.0 unless you want to change the zpt parameter.
A simple fitting example, there is no transit in the fake data I create here.
from ktransit import FitTransit
time = np.arange(0,10,0.0188) # you need a time and a flux
flux = np.zeros_like(time) # there are no transits here :(
ferr = np.ones_like(time) * 0.00001 # uncertainty on the data
fitT = FitTransit()
fitT.add_guess_star(rho=7.0)
fitT.add_guess_planet(
period=365.25, impact=0.0,
T0=0.0, rprs=0.009155)
fitT.add_data(time=time, flux=flux, ferr=ferr)
vary_star = ['rho', 'zpt'] # free stellar parameters
vary_planet = (['period', # free planetary parameters
'T0', 'impact',
'rprs']) # free planet parameters are the same for every planet you model
fitT.free_parameters(vary_star, vary_planet)
fitT.do_fit() # run the fitting
fitT.print_results() # print some resultsBest-fitting stellar parameters
rho: 5.02680657221
Best-fitting planet parameters for planet 0
impact: 0.0803742216248
period: 1.00000143821
T0: 1.49999092263
rprs: 0.0999846857972
Best-fitting planet parameters for planet 1
impact: 0.387103707996
period: 1.2000614885
T0: 0.700001081846
rprs: 0.0198908445261
The above comes from fitting a model with two planets in.
You can make a pretty plot of the output of the fitting by running
fig = ktransit.plot_results(time,flux,fitT.transitmodel)
fig.savefig('my_beautiful_fit.png') # or fig.show()The plot below shows a three planet model fit to the Kepler-37 data. This plot can be created using the code in examples/
Including radial velocity data
New in version 0.1, you can now include radial velocity data.
import ktransit
import matplotlib.pyplot as plt
import numpy
M = ktransit.LCModel()
M.add_star(
rho=1.5,
ld1=0.2,
ld2=0.4,
ld3=0.0,
ld4=0.0,
dil=0.0,
zpt=0.0,
veloffset=10 # new keyword, the radial velocity zero-point offset in m/s
)
M.add_planet(
T0=1.0,
period=1.0,
impact=0.1,
rprs=0.1,
ecosw=0.0,
esinw=0.0,
occ=0.0,
rvamp=100.) # radial velocity semi-amplitude in m/s
M.add_data(
time=numpy.arange(0,10,0.0188),
itime=numpy.zeros_like(numpy.arange(0,10,0.0188))+0.0188 )
M.add_rv(rvtime=numpy.arange(0,10,2.), # radial velocity observation timestamps
rvitime=numpy.zeros_like(numpy.arange(0,10,2.))+0.02 ) # integration time of each timestamp
tmod = M.transitmodel
rvmodel = M.rvmodel
plt.plot(M.time,tmod)
plt.plot(M.rvtime,rvmodel)
