Examples: Transit fit for a single planet
Note: This example shows you how to fit a single transiting planet for data taken in short cadence mode (i.e. short integration time). If you are dealing with long cadence data (e.g. Kepler/K2 data), you need to resample the model. Check here the corresponding example.
The simplest transiting planet model consists of a single planet, where the data was observed as short cadence (i.e., short integration time). In this example, we will model the light curve data of Pi Mensae (HD 39091) to estimate the transit parameters of Pi Mensae c (Gandolfi et al., 2018, Huang et al., 2018).
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/pimen
Now it is time to fill pimen directory with data and an input file. We will use the data used in Gandolfi et al., 2018. For details on how this data was obtained, we refer to that paper. You can download such data in this link.
Note how you have to 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
#---------------- -------- ------
8325.29698330 0.99919324 0.00021488
8325.29837220 0.99935260 0.00021488
8325.29976100 0.99939747 0.00021488
8325.30114990 0.99949504 0.00021488
...
...
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 Pi Men c
#Created by Oscar Barragan, December 2018.
#Light curve data file
fname_tr = ['pimen_lc_tess.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'
#We want the planet values in Earth units
unit_mass = 'earth'
#Stellar parameters taken from Gandolfi et al., (2018), A&A, 619, L10
mstar_mean = 1.02
mstar_sigma = 0.03
rstar_mean = 1.10
rstar_sigma = 0.01
tstar_mean = 5870.
tstar_sigma = 50.0
#Since we want to fit transit data, we need to include the next line
fit_tr = [True]
#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 fit eccentricity e, and angle of periastron w.
fit_t0 = ['u'] #uniform prior for T0
fit_P = ['u'] #uniform prior for P
fit_e = ['f'] #fixing e
fit_w = ['f'] #fixing w
fit_b = ['u'] #uniform prior for b
fit_a = ['u'] #uniform prior for a/R*
fit_rp = ['u'] #uniform prior for Rp/R*
fit_q1 = 'g' #gaussian prior for q1
fit_q2 = 'g' #gaussian prior for q2
#Set the prior ranges
#priors coming from Gandolfi et al., (2018)
#In practice, priors on T0 and P may come from detection algorithms, RV detections, etc.
#Minimum and maximum limits for T0
min_t0 = [8325.493694] #BJD
max_t0 = [8325.513694] #BJD
#Minimum and maximum limits for P
min_P = [6.256571] #days
max_P = [6.276571] #days
#Since we are fixing e, min_e = [0.0] will fix the eccentricity to zero, max_e is depreciated
min_e = [0.0] #this is fixed to zero in this case
max_e = [1.0] #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] #this is fixed to zero in this case
max_w = [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]
max_b = [1.0]
#Minimum and maximum limits for a/R*
min_a = [2.0]
max_a = [30.]
#Minimum and maximum limits for Rp/R*
min_rp = [0.0]
max_rp = [0.03]
#Gaussian priors on q1
min_q1 = 0.36 #Mean
max_q1 = 0.10 #Standard deviation
#Gaussian priors on q2
min_q2 = 0.25 #Mean
max_q2 = 0.10 #Standard deviation
#End of the input file for Pi Men c If you copy the input file example in the file inpy/pimen/input_file.py, you are ready to fit Pi Men c. In the main pyaneti directory type
./pyaneti.py pimen
The pyaneti execution will start and it should run in minutes in your laptop. Once the run finishes, open the outpy/pimen_out/pimenb_tr.pdf file. Does it look like this plot? If yes, congratulations, you have modelled Pi Men c!
How do your fitted parameters compare with the values reported by Gandolfi et al., 2018?
The last example corresponds to the simplest transiting planet model for a single planet. However, If you are dealing with long cadence data (e.g. Kepler/K2 data), you need to resample the model (see Kipping, 2010). This can be done easily with pyaneti by adding some extra lines in the input file. In this example, we will model the transits of K2-98b (Barragán et al., 2016).
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/K2-98
Now it is time to fill K2-98 directory with data and an input file. We will use the data used in Barragán et al., 2016. For details on how this data was obtained, we refer to that paper. You can download such data in this link.
Note how you have to 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
#---------------- -------- ------
7145.71976263 1.00001585 0.00013689
7145.74019450 1.00000581 0.00013689
7145.76062656 0.99986159 0.00013689
...
...
...
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 K2-98b
#Created by Oscar Barragan, October 2019.
#Light curve data file
fname_tr = ['K2-98-K2-flattened.dat']
#MCMC controls
thin_factor = 10
niter = 500
nchains = 100
#method to use 'mcmc' sample the posteriors, 'plot' creates the plots using posteriors from a previous run
method = 'mcmc'
#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 Barragán et al., (2016), AJ, 152, 193
mstar_mean = 1.07
mstar_sigma = 0.04
rstar_mean = 1.31
rstar_sigma = 0.08
tstar_mean = 6120.
tstar_sigma = 80.0
#Since we want to fit transit data, we need to include the next line
fit_tr = [True]
#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 fit eccentricity e, and angle of periastron w.
fit_t0 = ['u'] #uniform prior for T0
fit_P = ['u'] #uniform prior for P
fit_e = ['f'] #fixing e
fit_w = ['f'] #fixing w
fit_b = ['u'] #uniform prior for b
fit_a = ['u'] #uniform prior for a/R*
fit_rp = ['u'] #uniform prior for Rp/R*
fit_q1 = 'g' #gaussian prior for q1
fit_q2 = 'g' #gaussian prior for q2
#Set the prior ranges
#priors coming from Barragán et al., (2016)
#In practice, priors on T0 and P may come from detection algorithms, RV detections, etc.
#Minimum and maximum limits for T0
min_t0 = [7145.97] #BJD
max_t0 = [7145.99] #BJD
#Minimum and maximum limits for P
min_P = [10.130] #days
max_P = [10.140] #days
#Since we are fixing e, min_e = [0.0] will fix the eccentricity to zero, max_e is depreciated
min_e = [0.0] #this is fixed to zero in this case
max_e = [1.0] #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] #this is fixed to zero in this case
max_w = [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]
max_b = [1.0]
#Minimum and maximum limits for a/R*
min_a = [2.0]
max_a = [30.]
#Minimum and maximum limits for Rp/R*
min_rp = [0.0]
max_rp = [0.1]
#Gaussian priors on q1
min_q1 = 0.40 #Mean
max_q1 = 0.10 #Standard deviation
#Gaussian priors on q2
min_q2 = 0.26 #Mean
max_q2 = 0.10 #Standard deviation
#End of the input file for K2-98bNote: you can see that this file looks similar to the Pi Men c example, but with some extra lines where you have to specify the data cadence and the number of steps to perform the integration.
If you copy the input file example in the file inpy/K2-98/input_file.py, you are ready to fit K2-98b. In the main pyaneti directory type
./pyaneti.py K2-98
The pyaneti execution will start and it should run in minutes in your laptop. Once the run finishes, open the outpy/K2-98_out/K2-98b_tr.pdf file. Does it look like this plot? If yes, congratulations, you have modelled K2-98b!
How do your fitted parameters compare with the values reported in Barragán et al., 2016?
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 (
min_T0,max_T0) - Priors on Period (
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?
