Examples: Multi planet transit fitting
Note: This example assumes that you have seen the example fitting one single transiting planet and know to deal with short and long cadence data.
The system GJ 9827 contains (at least) three transiting planets. They were discovered by K2 on its Campaign 12 ( see Niraula et al., 2017, Prieto-Arranz et al., 2018, and Rodriguez et al., 2017 for more details).
First, let us make a directory to create our setup. If you have pyaneti installed, you only need to create a target directory in the inpy directory. If you are inside the pyaneti main directory, you only need to type
mkdir inpy/gj9827
Now it is time to fill gj9827 directory with data and an input file. We will use the data used in Niraula et al., 2017. For details on how this data was obtained, we refer to that paper. You can download such data in this link.
Note how you must create light curve data files for pyaneti. The first column corresponds to time (in days), the second column to the flux (the flux has to be normalized to 1), and the third column to the error of each measurement. A light curve file should look like this:
#---------------- -------- ------
#BJD - 2450000(d) Flux e_flux
#---------------- -------- ------
7738.39776200 1.00000711 0.00004331
7738.41819434 0.99993498 0.00004331
7738.43862668 1.00005847 0.00004331
7738.45905902 1.00001860 0.00004331
...
...
The next step is to create the input file to run pyaneti. This file has to include information about the light curve data file, MCMC controls, units to use, kind of priors, prior limits, etc.
The input file that you should use must look like this
#Input file for GJ9827
#Created by Oscar Barragan, October 2019.
#Light curve data file
fname_tr = ['GJ9827_lc.dat']
#MCMC controls
thin_factor = 1
niter = 500
nchains = 100
#method to use 'mcmc' sample the posteriors, 'plot' creates the plots using posteriors from a previous run
method = 'mcmc'
#method = 'plot'
#In order to re-sample the model, we just need to specify the cadence time (t_cad) in units of days and the number
#of steps for the integration (n_cad), default values are t_cad = 2./60./24. (2 min) and n_cad = 1 (i.e., no resampling)
#For K2 long cadence (30 min) data we will integrate the model for 30min with 10 steps
n_cad = 10
t_cad = 30./60./24.
#We want the planet values in Earth units
unit_mass = 'earth'
#Stellar parameters taken from Niraula et al., 2017, AJ, 154, 266
#https://iopscience.iop.org/article/10.3847/1538-3881/aa957c/meta
mstar_mean = 0.659
mstar_sigma = 0.060
rstar_mean = 0.651
rstar_sigma = 0.065
tstar_mean = 4255.
tstar_sigma = 110.
#Now we have to specify the number of planet signals that we are going to fit
#In this case we are dealing with 3 planets
nplanets = 3
#Since we want to fit transit data, we need to include the next line
#Now we are dealing with a multi-planet fit, we need to especify True for each fitted parameter
#In this case, we are fitting 3 planets, so let us write
fit_tr = [True,True,True]
#Given that now we are fitting 3 planets, we have to say to the code that we are NOT fitting RV data in this case, so we have to specify "False" for each planet
fit_rv = [False,False,False]
#Specify which kind of priors we want
#For a transit fit of a single planet with circular obit
#we fit only for time of minimum conjunction T0, period P, impact parameter b, scaled semi-major axis (a)
#scaled planet radius rp, and limb darkening coefficients q1 and q2
#We do not fit eccentricity e, and angle of periastron w.
#For the multi-planet case, we need to add an element to the list for each planet, in this case, we are dealing with 3 planets, so we write
fit_t0 = ['u','u','u'] #uniform prior for T0
fit_P = ['u','u','u'] #uniform prior for P
fit_e = ['f','f','f'] #fixing e
fit_w = ['f','f','f'] #fixing w
fit_b = ['u','u','u'] #uniform prior for b
fit_a = ['u','u','u'] #uniform prior for a/R*
fit_rp = ['u','u','u'] #uniform prior for Rp/R*
fit_k = ['f','f','f'] #No fit for semi-amplitude
#Note how we only indicate one element for the LD coefficients since they do not depend on the number of planets
fit_q1 = 'u' #gaussian prior for q1
fit_q2 = 'u' #gaussian prior for q2
#Set the prior ranges
#priors coming from Niraula et al., (2017)
#In practice, priors on T0 and P may come from detection algorithms, RV detections, etc.
#Minimum and maximum limits for T0
min_t0 = [7738.82,7738.50,7740.960] #BJD
max_t0 = [7738.83,7738.60,7740.970] #BJD
#Minimum and maximum limits for P
min_P = [1.207,3.645,6.201] #days
max_P = [1.210,3.655,6.203] #days
#Since we are fixing e, min_e = [0.0] will fix the eccentricity to zero, max_e is depreciated
min_e = [0.,0.,0.] #this is fixed to zero in this case
max_e = [1.,1.,1.] #the upper limit is not important if we are fixing a variable
#Since we are fixing w, min_w = [0.0] will fix the periastron to zero, max_w is depreciated
min_w = [0.,0.,0.] #this is fixed to zero in this case
max_w = [2*np.pi,2*np.pi,2*np.pi] #the upper limit is not important if we are fixing a variable
#Minimum and maximum limits for b
min_b = [0.,0.,0.]
max_b = [1.,1.,1.]
#Minimum and maximum limits for a/R*
min_a = [2.,2.,2.]
max_a = [30.,30.,50.]
#Minimum and maximum limits for Rp/R*
min_rp = [0.,0.,0.]
max_rp = [0.1,0.1,0.1]
#We are fixing k
min_k = [0.,0.,0.]
max_k = [1.,1.,1.]
#Uniform priors for LD following Kipping 2013
min_q1 = 0.
max_q1 = 1.
#Gaussian priors on q2
min_q2 = 0.
max_q2 = 1.
#End of the input file for GJ9827If you copy the input file example in the file inpy/gj9827/input_file.py, you are ready to fit the three planets of GJ9827. In the main pyaneti directory type
./pyaneti.py gj9817
The pyaneti execution will start and it should take some minutes in your laptop (now we are fitting 3 planets, this takes a bit more). Once the run finishes, you are ready to check if your three planets are there open the files: outpy/gj9827_out/gj9827b_tr.pdf, outpy/gj9827_out/gj9827c_tr.pdf, and outpy/gj9827_out/gj9827d_tr.pdf (note how pyaneti labelled each planet with a different letter). Do they look like this plot? If yes, congratulations, you have modelled the multi-planet transiting system GJ9827!
How do your fitted parameters compare with the values reported in Niraula et al., 2017?
You can use the same input file presented in this example to fit your own data for all the planets that you want. You just have to be sure that you are indicating to the code the number of planets to fit. The main parameters you should change from the example file to perform a fit to your own transit data are system dependents. They are
- Priors on Time of minimum conjunction for each planet (
min_T0,max_T0) - Priors on Period for each planet (
min_P,max_P) - priors on Limb Darkening Coefficients
- Name of your light curve file (
fname_tr) - Stellar parameters (They do not affect the MCMC run, but the physical parameter estimations)
- Be sure that your input light curve is normalized to one.
- How are you dealing with the limb darkening coefficients? (priors: uniform or Gaussian, do not fix them!).
- Are you dealing with long-cadence data?
TRAPPIST-1 is small star hosting 7 transiting planets. This has been one of the most exciting exoplanet discoveries in the last years.
In this link, you can find the K2 short-cadence data of EPIC 246199087 (also known as TRAPPIST-1). A great analysis of TRAPPIST-1's K2 data is presented in Luger et al., 2017. Can you create an input file to fit all the planets in the system?
Here you can see the pyaneti fit for the seven planets.
