Lightcurve fitting with juliet

We have already exemplified how to fit a basic transit lightcurve in the Getting started section with juliet. Here, however, we explore some interesting extra features of the lightcurve fitting process, including limb-darkening laws, parameter transformations and fitting of data from multiple-instruments simultaneously, along with useful details on the model evaluations with juliet.

Before going into the tutorial, it is useful to first understand the lightcurve model that juliet uses. In the absence of extra linear terms (which we will deal with in the Incorporating linear models tutorial), a juliet lightcurve model for a given instrument \(i\) is given by (see Section 2 of the juliet paper)

\(\mathcal{M}_{i}(t) + \epsilon_i(t)\),


\(\mathcal{M}_{i}(t) = [\mathcal{T}_{i}(t)D_{i} + (1-D_{i})]\left(\frac{1}{1+D_{i}M_{i}}\right)\)

is the photometric model composed of the dilution factor \(D_{i}\), the relative out-of-transit target flux \(M_{i}\), and the transit model for the instrument \(\mathcal{T}_{i}(t)\) (defined by the transit parameters and by the instrument-dependant limb-darkening coefficients — see the Models, priors and outputs section for details). Here, \(\epsilon_i(t)\) is a stochastic process that defines a “noise model” for the dataset. In this section we will assume that \(\epsilon_i(t)\) is white-gaussian noise, i.e., \(\epsilon_i(t)\sim \mathcal{N}(0,\sqrt{\sigma_i(t)^2 + \sigma_{w,i}^2})\), where \(\mathcal{N}(\mu,\sigma)\) denotes a normal distribution with mean \(\mu\) and standard-deviation \(\sigma\), and where \(\sigma_i(t)\) are the errors on the datapoint at time \(t\) and \(\sigma_{w,i}\) is a so-called “jitter” term. We deal with more general noise models in the Incorporating Gaussian Processes tutorial.

The juliet lightcurve model is a bit different than typical lightcurve models which typically only fit for an out-of-transit flux offset. The first difference is that our model includes a dilution factor \(D_{i}\) which allows the user to account for possible contaminating sources in the aperture that might produce a smaller transit depth than the real one. In fact, if there are \(n\) sources with fluxes \(F_n\) in the aperture and the target has a flux \(F_T\), then one can show (see Section 2 of the juliet paper) that the dilution factor can be interpreted as

\(D_i = \frac{1}{1 + \sum_n F_n/F_T}\).

A dilution factor of 1, thus, implies no external contaminating sources. The second difference is that the relative out-of-transit target flux, \(M_i\) — which from hereon we refer to as the “mean out-of-transit flux” — is a multiplicative term and not an additive offset. This is because input lightcurves are usually normalized (typically with respect to the mean or the median), and in theory a simple additive offset might thus still not account for this pre-normalization of the lightcurve. \(M_i\), in turn, has a well defined interpretation: if the real out-of-transit flux was \(F_T + \sum_n F_n + E\), where \(E\) is an offset flux given by light coming not from the sources \(F_n\) or from the target, \(F_T\) (e.g., background scattered light, a bias flux, etc.), then this term can be interpreted as \(E/F_T\). As can be seen, then, with \(D_i\) and \(M_i\), one can uniquely recover the real (relative to the target) fluxes.

Transit fits

To showcase the ability of juliet to fit transit lightcurves, we will play with the HATS-46b system (Brahm et al., 2017), as the TESS data for this system has interesting features that we will be using both in this tutorial and in the Incorporating Gaussian Processes tutorial. In this tutorial in particular, we will play with the data obtained in Sector 2, because it seems the level of variability/systematics in this particular dataset is much smaller than the one for Sector 1 (which we tackle in the Incorporating Gaussian Processes tutorial). First, let us download and plot the TESS data, taking the opportunity to also put the data in dictionaries so we can feed it to juliet:

import juliet
import numpy as np
import matplotlib.pyplot as plt

# First, get arrays of times, normalized-fluxes and errors for HATS-46
# from Sector 1 from MAST:
t, f, ferr  = juliet.get_TESS_data(''+\

# Put data arrays into dictionaries so we can fit it with juliet:
times, fluxes, fluxes_error = {},{},{}
times['TESS'], fluxes['TESS'], fluxes_error['TESS'] = t,f,ferr

# Plot data:
plt.errorbar(t, f, yerr=ferr, fmt='.')
Sector 2 data for HATS-46b.

Pretty nice dataset! The transits can be clearly seen by eye. The period seems to be about \(P \sim 4.7\) days, in agreement with the Brahm et al., 2017 study, and the time-of-transit center seems to be about \(t_0 \sim 1358.4\) days. Let us now fit this lightcurve using these timing constraints as priors. We will use the same “non-informative” priors for the rest of the transit parameters as was already done for TOI-141b in the Getting started tutorial:

priors = {}

# Name of the parameters to be fit:
params = ['P_p1','t0_p1','r1_p1','r2_p1','q1_TESS','q2_TESS','ecc_p1','omega_p1',\
              'rho', 'mdilution_TESS', 'mflux_TESS', 'sigma_w_TESS']

# Distributions:
dists = ['normal','normal','uniform','uniform','uniform','uniform','fixed','fixed',\
                 'loguniform', 'fixed', 'normal', 'loguniform']

# Hyperparameters
hyperps = [[4.7,0.1], [1358.4,0.1], [0.,1], [0.,1.], [0., 1.], [0., 1.], 0.0, 90.,\
                   [100., 10000.], 1.0, [0.,0.1], [0.1, 1000.]]

# Populate the priors dictionary:
for param, dist, hyperp in zip(params, dists, hyperps):
    priors[param] = {}
    priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

Now let’s fit the dataset with juliet, saving the results to a folder called hats46:

# Load and fit dataset with juliet:
dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, \
                   yerr_lc = fluxes_error, out_folder = 'hats46')

results =

As was already shown in the Getting started tutorial, it is easy to plot the juliet fit results using the function. In the background, this function extracts by default nsamples=1000 random samples from the joint posterior distribution of the parameters and evaluates the model using them — by default, a call to this function given an instrument name returns the median of all of those models. However, one can also retrieve the models that are about “1-sigma away” from this median model — i.e., the 68% credibility band of these models — by setting return_err=True. One can actually select the percentile credibility band with the alpha parameter (by default, alpha=0.68). Let us extract and plot the median model and the corresponding 68% credibility band around it using this function. We will create two plots: one of time versus flux, and another one with the phased transit lightcurve:

# Extract median model and the ones that cover the 68% credibility band around it:
transit_model, transit_up68, transit_low68  ='TESS', return_err=True)

# To plot the phased lighcurve we need the median period and time-of-transit center:
P, t0 = np.median(results.posteriors['posterior_samples']['P_p1']),\

# Get phases:
phases = juliet.get_phases(dataset.times_lc['TESS'], P, t0)

import matplotlib.gridspec as gridspec

# Plot the data. First, time versus flux --- plot only the median model here:
fig = plt.figure(figsize=(12,4))
gs = gridspec.GridSpec(1, 2, width_ratios=[2,1])
ax1 = plt.subplot(gs[0])

ax1.errorbar(dataset.times_lc['TESS'], dataset.data_lc['TESS'], \
             yerr = dataset.errors_lc['TESS'], fmt = '.' , alpha = 0.1)

# Plot the median model:
ax1.plot(dataset.times_lc['TESS'], transit_model, color='black',zorder=10)

# Plot portion of the lightcurve, axes, etc.:
ax1.set_xlabel('Time (BJD - 2457000)')
ax1.set_ylabel('Relative flux')

# Now plot phased model; plot the error band of the best-fit model here:
ax2 = plt.subplot(gs[1])
ax2.errorbar(phases, dataset.data_lc['TESS'], \
             yerr = dataset.errors_lc['TESS'], fmt = '.', alpha = 0.3)
idx = np.argsort(phases)
ax2.plot(phases[idx],transit_model[idx], color='black',zorder=10)
Sector 2 data for HATS-46b along with the best-fit juliet model.

As can be seen, the lightcurve model is quite precise! In the code above we also made use of a function and a dictionary which we have not introduced in their entirety yet. The first is the juliet.get_phases(t, P, t0) function — this gives you back the phases at the times t given a period P and a time-of-transit center t0. The second is a very important dictionary: it was already briefly introduced in the Models, priors and outputs section, but this introduction did not pay justice to its importance. This is the results.posteriors dictionary. The posterior_samples key of this dictionary stores the posterior distribution of the fitted parameters — we make use of this dictionary in detail in the next part of the tutorial.

Transit parameter transformations

In the fit done in the previous section we fitted the Sector 2 TESS lightcurve of HATS-46b. There, however, we fitted for the transformed parameters r1_p1 and r2_p1 which parametrize the planet-to-star radius ratio, \(p = R_p/R_*\), and the impact parameter, in our case given by \(b = (a/R_*)\cos i\), and the limb-darkening parametrization q1_TESS and q2_TESS, which in our case parametrize the coefficients \(u_1\) and \(u_2\) of the quadratic limb-darkening law. How do we transform the posterior distributions of those parametrizations, stored in the results.posteriors['posterior_samples'] dictionary back to their physical parameters? juliet has built-in functions to do just this.

To transform from the \((r_1,r_2)\) plane to the \((b,p)\) plane, we have implemented the transformations described in Espinoza (2018). These require one defines the minimum and maximum allowed planet-to-star radius ratio — by default, within juliet the parametrization allows to search for all planet-to-star radius ratios from \(p_l = 0\) to \(p_u = 1\) (and these can be modified in the fit object — e.g.,,pl= 0.0, pu = 0.2)). The values used for each fit are always stored in results.posteriors['pl'] and results.posteriors['pu']. In our case, then, to obtain the posterior distribution of \(b\) and \(p\), we can use the juliet.utils.reverse_bp(r1,r2,pl,pu) function which takes samples from the \((r_1,r_2)\) plane and converts them back to the \((b,p)\) plane. Let us do this transformation for the HATS-46b fit done above and compare with the results obtained in Brahm et al., 2017:

fig = plt.figure(figsize=(5,5))
# Store posterior samples for r1 and r2:
r1, r2 = results.posteriors['posterior_samples']['r1_p1'],\

# Transform back to (b,p):
b,p = juliet.utils.reverse_bp(r1, r2, 0., 1.)

# Plot posterior distribution:

# Extract median and 1-sigma errors for b and p from
# the posterior distribution:
bm,bu,bl = juliet.utils.get_quantiles(b)
pm,pu,pl = juliet.utils.get_quantiles(p)

# Plot them:
             xerr = np.array([[bu-bm,bm-bl]]),\
             yerr = np.array([[pu-pm,pm-pl]]),\
             fmt = 'o', mfc = 'white', mec = 'black',\
             ecolor='black', ms = 15, elinewidth = 3, \
             zorder = 5, label = 'This work')

# Plot values in Brahm et al. (2017):
             xerr = np.array([[0.042,0.034]]), \
             yerr = np.array([0.0027]),zorder = 5,\
             label = 'Brahm et al. (2017)', fmt='o', \
             mfc = 'white', elinewidth = 3, ms = 15)

plt.ylabel('$p = R_p/R_*$')
Posterior distribution of the impact parameter and planet-to-star radius ratio.

The agreement with Brahm et al., 2017 is pretty good! The planet-to-star radius ratios are consistent within one-sigma, and the (uncertain for TESS) impact parameter is consistent at less thant 2-sigma with the one published in that work.

What about the limb-darkening coefficients? juliet also has a built-in function to perform the inverse transformation in order to obtain them — this is the juliet.utils.reverse_ld_coeffs() function — given a limb-darkening law and the parameters \(q_1\) and \(q_2\), this function gives back the limb-darkening coefficients \(u_1\) and \(u_2\). Let us plot the posterior distribution of the limb-darkening coefficients; let’s compare them to theoretical limb-darkening coefficients using limb-darkening (Espinoza & Jordan, 2015):

fig = plt.figure(figsize=(5,5))
# Store posterior samples for q1 and q2:
q1, q2 = results.posteriors['posterior_samples']['q1_TESS'],\

# Transform back to (u1,u2):
u1, u2 = juliet.utils.reverse_ld_coeffs('quadratic', q1, q2)

# Plot posterior distribution:

# Plot medians and errors implied by the posterior:
u1m,u1u,u1l = juliet.utils.get_quantiles(u1)
u2m,u2u,u2l = juliet.utils.get_quantiles(u2)
             xerr = np.array([[u1u-u1m,u1m-u1l]]),\
             yerr = np.array([[u2u-u2m,u2m-u2l]]),\
             fmt = 'o', mfc = 'white', mec = 'black',\
             ecolor='black', ms = 13, elinewidth = 3, \
             zorder = 5, label = 'This work')



Posterior distribution of the limb-darkening coefficients for HATS-46.

The agreement with the theory is pretty good in this case! It was kind of expected — HATS-46 is a solar-type star after all. Notice the triangular shape of the parameter spaced explored? This is what the \((q_1,q_2)\) sampling is expected to sample — the triangle englobes all the physically plausible parameter space for the limb-darkening coefficients (positive, decreasing-to-the-limb limb-darkening profiles). For details, see Kipping (2013).

Fitting multiple datasets

In the previous sections we have been fitting the TESS data only. What if we want to add extra datasets and fit all of them jointly in order to extract the posterior distribution of the transit parameters? As it was already mentioned, this is very easy to do with juliet: you simply add new elements/keys to the dictionary one gives as inputs to it. Of course, you also have to add some extra priors for the extra instruments: in particular, one has to define a jitter (\(\sigma_{w,i}\)), dilution factor (\(D_i\)), mean out-of-transit flux (\(M_i\)) and limb-darkening parametrization (\(q_1\) if a linear law wants to be assumed, or also give \(q_2\) if a quadratic law wants to be used). Let us fit the TESS data together with the follow-up lightcurves obtained by Brahm et al., 2017 from the Las Cumbres Observatory Global Telescope Network (LCOGT) and the 1m Swope Telescope. These can be obtained from CDS following the paper link, but we have uploaded them here and here so it is easier to follow this tutorial. Once that data is downloaded, we can load this data in juliet as follows:

# Add LCOGT and SWOPE data to the times, fluxes and fluxes_error dictionary.
# Fill also the priors for these instruments:
for instrument in ['LCOGT','SWOPE']:
    # Open dataset files, extract times, fluxes and errors to arrays:
    t2,f2,ferr2 = np.loadtxt('hats-46_data_'+instrument+'.txt',\
    # Add them to the data dictionaries which already contain the TESS data (see above):
    times[instrument], fluxes[instrument], fluxes_error[instrument] = \
                                                t2-2457000, f2, ferr2

    # Add priors to the already defined ones above for TESS, but for the other instruments:
    params = ['sigma_w_','mflux_','mdilution_','q1_','q2_']
    dists = ['loguniform', 'normal', 'fixed', 'uniform', 'uniform']
    hyperps = [[0.1,1e5], [0.0,0.1], 1.0, [0.0,1.0], [0.0,1.0]]

    for param, dist, hyperp in zip(params, dists, hyperps):
        priors[param+instrument] = {}
        priors[param+instrument]['distribution'], \
        priors[param+instrument]['hyperparameters'] = dist, hyperp

And with this one can simply run a juliet fit again:

dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, \
                      yerr_lc = fluxes_error, out_folder = 'hats46-extra')

results =

This can actually take a little bit longer than just fitting the TESS data (a couple of extra minutes) — it is a 17-dimensional problem after all. Let us plot the results of the joint instrument fit:

# Extract new period and time-of-transit center:
P,t0 =  np.median(results.posteriors['posterior_samples']['P_p1']),\

# Generate arrays to super-sample the models:
model_phases = np.linspace(-0.04,0.04,1000)
model_times = model_phases*P + t0

# Plot figure:
fig = plt.figure(figsize=(10,3))
instruments = ['TESS','LCOGT','SWOPE']
alphas = [0.1, 0.5, 0.5]
for i in range(3):
    instrument = instruments[i]
    # Plot phase-folded data:
    phases = juliet.utils.get_phases(dataset.times_lc[instrument], P, t0)
    plt.errorbar(phases, dataset.data_lc[instrument], \
                 yerr = dataset.errors_lc[instrument], fmt = '.' , alpha = alphas[i])
    # Evaluate model in the supersampled times, plot on top of data:
    model_lc =, t = model_times)
    if i == 0:
        plt.ylabel('Relative flux')
Final fit involving TESS, Swope and LCOGT data.

Pretty nice fit! The Swope data actually shows a little bit more scatter — indeed, the \(\sigma_{w,SWOPE} = 1269^{+185}_{-155}\) ppm, which indicates that there seems to be some extra process happening in the lightcurve (e.g., systematics), which are being modelled in our fit with a simple jitter term. So, how does the posteriors of our parameters compare with that of the TESS-only fit? We can repeat the plot made above for the planet-to-star radius ratio and impact parameter to check:

Posterior distribution of a TESS-only fit and a TESS+LCOGT+Swope fit; compared with Brahm et al. (2017)

Interesting! The transit depth is consistent between fits and with the work of Brahm et al., 2017. Interestingly, the impact parameter is practically the same as the TESS-only fit, and just shrinked a little bit. It is still consistent at 2-sigma with the work of Brahm et al., 2017, however.

A word on limb-darkening and model selection

Throughout the tutorial, we have not explicitly defined what limb-darkening laws we wanted to use for each dataset. By default, juliet assumes that if the user defines \(q_1\) and \(q_2\), then a quadratic law wants to be used, whereas if the user only gives \(q_1\), a linear-law is assumed. In general, the limb-darkening law to use depends on the system under study (see, e.g., Espinoza & Jordan, 2016.), and thus the user might want to use laws other than the ones that are pre-defined by juliet. This can be easily done when loading a dataset via juliet.load using the ld_laws flag. This flag receives a string with the name of the law to use — currently supported laws are the linear, the quadratic, the logarithmic and the squareroot laws. We don’t include the exponential law in this list as it has been shown to be a non-physical law in Espinoza & Jordan, 2016.

Let us test how the different laws do on the TESS dataset of HATS-46b. For this, let us fit the dataset with all the available limb-darkening laws and check the log-evidences, \(\ln \mathcal{Z} = \ln \mathcal{P}(D | \textrm{Model})\) each model gives. Assuming all the models are equally likely, the different log-evidences can be transformed to odds ratios (i.e., the ratio of the probabilities of the models given the data, \(\mathcal{P}(\textrm{Model}_i|D)/\mathcal{P}(\textrm{Model}_j|D)\)) by simply substracting the log-evidences of the different models, i.e.,

\(\ln \frac{\mathcal{P}(\textrm{Model}_i|D)}{\mathcal{P}(\textrm{Model}_j|D)} = \ln \frac{\mathcal{P}(D | \textrm{Model}_i)}{\mathcal{P}(D|\textrm{Model}_j)} = \ln \frac{Z_i}{Z_j}\),

if \(P(\textrm{Model}_i)/P(\textrm{Model}_j) = 1\). juliet also extracts the model evidences in the results.posteriors dictionary under the lnZ key; errors on this log-evidence calculation are under lnZerr. Let us compute the log-evidences for each limb-darkening law and compare them to see which one is the “best” in terms of this model comparison tool:

# Load Sector 1 data for HATS-46b again:
t, f, ferr  = juliet.get_TESS_data(''+\

# Put data arrays into dictionaries so we can fit it with juliet:
times, fluxes, fluxes_error = {},{},{}
times['TESS'], fluxes['TESS'], fluxes_error['TESS'] = t,f,ferr

# Define limb-darkening laws to test:
ld_laws = ['linear','quadratic','logarithmic','squareroot']

for ld_law in ld_laws:
     priors = {}
     # If law is not the linear, set priors for q1 and q2. If linear, set only for q1:
     if ld_law != 'linear':
         params = ['P_p1','t0_p1','r1_p1','r2_p1','q1_TESS','q2_TESS','ecc_p1','omega_p1',\
                   'rho', 'mdilution_TESS', 'mflux_TESS', 'sigma_w_TESS']

         dists = ['normal','normal','uniform','uniform','uniform','uniform','fixed','fixed',\
                      'loguniform', 'fixed', 'normal', 'loguniform']

         hyperps = [[4.7,0.1], [1358.4,0.1], [0.,1], [0.,1.], [0., 1.], [0., 1.], 0.0, 90.,\
                        [100., 10000.], 1.0, [0.,0.1], [0.1, 1000.]]
         params = ['P_p1','t0_p1','r1_p1','r2_p1','q1_TESS','ecc_p1','omega_p1',\
                   'rho', 'mdilution_TESS', 'mflux_TESS', 'sigma_w_TESS']

         dists = ['normal','normal','uniform','uniform','uniform','fixed','fixed',\
                      'loguniform', 'fixed', 'normal', 'loguniform']

         hyperps = [[4.7,0.1], [1358.4,0.1], [0.,1], [0.,1.], [0., 1.], 0.0, 90.,\
                        [100., 10000.], 1.0, [0.,0.1], [0.1, 1000.]]

     for param, dist, hyperp in zip(params, dists, hyperps):
         priors[param] = {}
         priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

     dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, \
                           yerr_lc = fluxes_error, out_folder = 'hats46-'+ld_law, \
                           ld_laws = ld_law)

     results =
     print("lnZ for "+ld_law+" limb-darkening law is: ",results.posteriors['lnZ']\

In our runs this gave:

lnZ for linear limb-darkening law is:      64202.653 +- 0.040
lnZ for quadratic limb-darkening law is:   64202.182 +- 0.018
lnZ for logarithmic limb-darkening law is: 64202.652 +- 0.077
lnZ for squareroot limb-darkening law is:  64202.786 +- 0.041

At face value, the model with the largest log-evidence is the square-root law, whereas the one with the lowest log-evidence is the quadratic law. However, the difference between those two log-evidences is very small: only \(\Delta \ln Z = 0.60\) in favor of the square-root law, or an odds ratio between those laws of \(\exp\left(\Delta \ln Z\right) \approx 2\) — given the data, the square-root law model is only about two times more likely than the quadratic law. Not much, to be honest — I wouldn’t bet my money on the quadratic law being wrong, so our assumption of a quadratic limb-darkening law in our analyses above seems to be very good. It is unlikely more complex limb-darkening laws would have given better results, by the way: note how the simpler linear law is basically equally likely to the square-root law (\(\exp\left(\Delta \ln Z\right) \approx 1\)).

What if more than one instrument is being fit; how do we define limb-darkening laws for each instrument? The ld_laws flag can also take as input a comma-separated string where one indicates the law to be used for each instrument in the form instrument-ldlaw. For example, if we wanted to fit the TESS, LCOGT and Swope data and define a square-root law for the former and logarithmic law for the other instruments, we would do (assuming we have already loaded the data and priors to the priors, times, fluxes and fluxes_error dictionaries):

dataset = juliet.load(priors=priors, t_lc = times, y_lc = fluxes, \
                      yerr_lc = fluxes_error,\
                      ld_laws = 'TESS-squareroot,LCOGT-logarithmic,SWOPE-logarithmic')

results =