Examples: RV fit for a single planet
One of the simplest cases to model consist of a planet detected in RV data with a single spectrograph. For the example case we will model the RV data of 51 Peg to infer the parameters of 51 Peg b, the first exoplanet detected around a Sun-like star.
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/51Peg
The next step is to get the RV data. Birkby et al., 2017 compiled the RV measurements of 51 Peg and you can find it here. That data file include data taken with the spectrograph ELODIE. Some of this data are part of the famous paper of 1995. You can find the data set that we will use here. Just download it and put it inside inpy/51Peg.
Note how you have to create RV data files for pyaneti. First column corresponds to time (in days), second column to the RV (in km/s), third column to the error of each RV measurement (in km/s) and a fourth column to specify the spectrograph.
#---------------- --------------- ----------- ------
#BJD (d) RV (km/s) e_RV (km/s) Spectrograph
#---------------- --------------- ----------- ------
2449610.53276 -33.258 0.009 ELODIE
2449612.47166 -33.225 0.009 ELODIE
2449655.31126 -33.272 0.007 ELODIE
...
...
The next step is to create the input file to run pyaneti. This file has to include information about the spectrograph label in the RV 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 51 Peg b
#Created by Oscar Barragan, December 2018.
#This contains the instrument name included in the fourth column of the RV data file
telescopes = ['ELODIE']
#This corresponds to the label that you want to put to each instrument
telescopes_labels = ['ELODIE Spectrograph']
#RV data file
fname_rv = ['51pegbrvkms.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'
#We want the planet values in Jupiter units
unit_mass = 'jupiter'
#Stellar data taken from http://exoplanet.eu/catalog/51_peg_b/
#Stellar mass
mstar_mean = 1.11
mstar_sigma = 0.06
#Stellar radius
rstar_mean = 1.266
rstar_sigma = 0.046
#Stellar temperature
tstar_mean = 5793.
tstar_sigma = 70.
#Since we want to fit RV data, we need to include the next line
fit_rv = [True]
#Specify which kind of priors we want
#For a RV fit of a single planet with circular obit
#we fit only for time of minimum conjunction T0, period P, Doppler semi-amplitude k and systemic velocities v0
#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'] #We fix e
fit_w = ['f'] #We fix w
fit_k = ['u'] #Uniform prior for k
fit_v0 = 'u' #Unform prior for v0
#Set the prior ranges
#For 51 Peg b I set priors based on previous published values
#In practice, priors on T0 and P may come from periodogram analysis, transit ephemeris, etc.
#Minimum and maximum limits for T0
min_t0 = [2449798.7] #BJD
max_t0 = [2449799.1] #BJD
#Minimum and maximum limits for P
min_P = [4.22] #days
max_P = [4.24] #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 k
#I set max_k = 1.0 km/s, but be aware that this limit depends on the signal that you will fit.
min_k = [0.0] #km/s
max_k = [1.0] #km/s
#End of input file for 51 Peg b If you look the input file again, you will notice that we did not specify priors for the systemic velocity. This is because pyaneti creates priors for this automatically.
If you copy the input file example in the file inpy/51Peg/input_fit.py, you are ready to fit the RV data for 51 Peg b. In the main pyaneti directory type
./pyaneti.py 51Peg
The pyaneti execution will start and it should run fast in your laptop. Once the run finishes, open the outpy/51Peg_out/51Pegb_rv.pdf file. Does it look like this plot? If yes, congratulations, you have modelled the signal of 51 Peg b.
How do your fitted parameters compare with the original values reported in the first 51 Peg b paper?
Given that 51 Peg host a famous exoplanet, it has been observed with more than one spectrograph. Fortunately, pyaneti can deal with this kind of data sets by changing only two lines in our previous input file. If you look inside the 51pegbrvkms.dat file, you may notice that there is more than one spectrograph. The available spectrographs are ELODIE, Lick6, Lick8, Lick13, HIRES, and HARPS.
Now re-open the input file for 51 Peg and change the telescopes and telescopes_labels to
...
#This contains the instrument name included in the fourth column of the RV data file
telescopes = ['ELODIE','Lick6','Lick8','Lick13','HIRES','HARPS']
#This corresponds to the label that you want to put to each instument
telescopes_labels = ['ELODIE','Lick6','Lick8','Lick13','HIRES','HARPS']
...What about the priors for each offset? In this case, you do not need to specify the priors for each spectrograph. Pyaneti creates an automatic prior for each offset (only if the variable fit_v0 is set to 'u'). Now just re-run pyaneti and it will fit an offset for each different instrument.
Now your RV curve should look like this:
where each symbol/color correspond to a different instrument.
The main parameters you should change from the example file to perform a fit to your own RV data are:
- Priors on Time of minimum conjunction (
min_T0,max_T0) - Priors on Period (
min_P,max_P) - Name of your RV file (
fname_rv) - Telescope labels (
telescopesandtelescopes_labelsaccording to your RV input file) - Stellar parameters (They do not affect the MCMC run, but the minimum mass calculation)
There are other parameters than you may want to change to optimise the MCMC run (e.g., priors on Doppler semi-amplitude).
