Fitting the Light Curves

PACMAN can fit time series observations of exoplanet transits, eclipses, and phase curves. It includes models for both the astrophysical signal and instrument systematic noise (e.g. ramps and slopes).

Here we describe the currently implemented models. Initial guesses and priors for the model parameters are listed in the fit_par.txt file.

The full model, which will be fit to the dataset, is the product of the physical and systematic flux. This ensures that we fit for the physical and systematic components simultaneously. For an example see Kreidberg et al., 2018.

Models

Instrument Systematics

Note

c is in log10. An average flux of 10^7 photoelectrons therefore leads to approximately c = 7.

  • model_ramp.py

    Exponential ramps fit to each orbit (free parameters: r1, r2, r3)

  • polynomial1.py

    Linear slope (free parameters: v)

  • polynomial2.py

    Quadratic trend (free parameters: v, v2)

  • exponential_visit.py

    Exponential trend over the whole visit (free parameters: exp1, exp2)

  • logarithmic_visit.py

    Logarithmic trend over the whole visit (free parameters: log1, log2)

  • upstream_downstream.py

    Offset between scan directions due to the upstream-downstream effect (free parameters: scale)

  • divide_white.py

    Uses the divide-white method, which assumes that the systematic parameters for the spectroscopic light curves are the same (have the same shape) as for the white light curve. See equation 2 in Kreidberg et al. (2014) for a reference. This model does not have any additional free parameters (so nothing has to be added to the fit_par.txt file that used).

  • constants_cj.py

    Alternative for model_ramp which fits a constant to every j-th exposure in an orbit. See equation 1 in Kreidberg et al. (2019) and references within for an application of this model.

Note

c is in log10. An average flux of 10^7 photoelectrons therefore leads to approximately c = 7.

  • uncmulti.py

    Scales the errorbars at every iteration of the sampler (free parameter: uncmulti_val). Does not work for the least squares fit only for the samplers, mcmc and dynesty.

Warning

With the current version of PACMAN, uncmulti_val has to be entered into fit_par.txt as the last parameter!!

This function has been implemented as an alternative to the ‘rescale_uncert’ technique. This rescales the errorbars of the flux measurements after the least squares routine so that chi2_red = 1. This might be problematic however, if the least squares is having troubles finding a good solution. Then the errorbars would be overestimated.

Astrophysical

  • transit.py

    Planetary transit (free parameters: t0, per, rp, a, inc, ecc, w, u1, u2). The parameters follow the parametrization from batman.

  • t0: transit midtime in days. Note that toffset from the pcf will be subtracted from this parameter to avoid fitting for high numbers.

  • per: period in days

  • rp: planet to star ratio Rp/Rs. unitless.

  • a: semi-major-axis to star radius ratio a/Rs. unitless.

  • inc: inclination in degree

  • ecc: eccentricity. unitless

  • w: argument of periastron in degrees. typically set to 90° if ecc = 0.

  • u1, u2: limb darkening

  • eclipse.py

    Secondary eclipse (free parameters: t_secondary, per, rp, fp, a, inc, ecc, w) The parameters follow the parametrization from batman.

  • t_secondary: eclipse midtime in days.

  • per: period in days

  • rp: planet to star ratio Rp/Rs. unitless.

  • fp: planet to star flux Fp/Fs. unitless

  • a: semi-major-axis to star radius ratio a/Rs. unitless.

  • inc: inclination in degree

  • ecc: eccentricity. unitless

  • w: argument of periastron in degrees. typically set to 90° if ecc = 0.

  • sine1.py

    Sinusoid (free parameters: a1, omega1, phi1)

  • sine2.py

    Sum of three sinusoids (free parameters: a1, omega1, phi1, a2, omega2, phi2, a3, omega3, phi3, a12, omega12, phi12, a22, omega22, phi22, a32, omega32, phi32)

  • sine_curve.py

    Sum of two sinusoids (free parameters: amp1, theta1, per, amp2, theta2)

The fit_par.txt file

This file has to be set up when running Stage 30. Here’s an example:

#parameter     fixed   tied   value       lo_lim   lo_val    hi_lim   hi_val   prior   p1        p2       step_size  
per            True    -1     1.5804046   True     1.55      True     1.6      X       2.25314   2e-05    1e-07      
t0             False   -1     0.18        True     0.16      True     0.2      U       0.0       0.3      0.001      
t_secondary    True    -1     0.0         False    0.0       False    0.0      X       0.0       0.0      0.0        
w              True    -1     90.0        False    0.0       False    0.0      X       0.0       0.0      0.0        
a              True    -1     15.23       True     14.0      True     16.0     X       4.98      0.05     0.02       
inc            True    -1     89.1        True     88.0      True     89.9     X       85.3      0.2      0.02       
rp             False   -1     0.116       True     0.05      True     0.2      U       0.01      0.3      0.001      
fp             True    -1     0.0         False    0.0       False    0.002    X       0.0       1.0      0.0        
u1             False   -1     0.29        True     0.0       True     0.6      U       0.0       1.0      0.01       
u2             True    -1     0.0         False    0.0       False    0.0      X       0.0       0.0      0.0        
ecc            True    -1     0.0         False    0.0       False    0.0      X       0.0       0.0      0.0        
c              False   0      8.37        True     6.0       True     8.9      X       6.7       7.0      0.001      
c              False   1      8.37        True     6.0       True     8.9      X       6.7       7.0      0.001      
v              False   0      -1e-06      False    -0.0001   False    -1e-08   X       -7e-06    -1e-06   1e-05      
v              False   1      -1e-06      False    -0.0001   False    -1e-08   X       -7e-06    -1e-06   1e-05      
v2             True    -1     0.0         True     0.0       True     1.0      X       0.0       0.0      1e-11      
r1             False   0      0.1         False    0.0       False    1.0      U       -10.0     10.0     0.01       
r1             False   1      0.1         False    0.0       False    1.0      U       -10.0     10.0     0.01       
r2             False   0      0.0         False    5.5       False    7.5      U       -100.0    100.0    0.1        
r2             False   1      0.0         False    5.5       False    7.5      U       -100.0    100.0    0.1        
r3             True    0      0.0         False    -1.0      False    1.0      U       -10.0     10.0     0.001      
r3             True    1      0.0         False    -1.0      False    1.0      U       -10.0     10.0     0.001      
scale          False   0      0.0         False    -0.01     False    0.01     U       -0.1      0.1      0.002      
scale          False   1      0.0         False    -0.01     False    0.01     U       -0.1      0.1      0.002      
uncmulti_val   True    -1     1.0         False    1.0       False    7.0      U       0.1       10.0     0.01       

Let’s have a look at each column:

  • parameter

    A list of the different parameters in each model is above.

  • fixed

    If set to False, the parameter will be a free parameter.

  • tied

    If the user wants to tie a parameter over all visits, set -1.

    If the user does not want to tie a certain parameter, he or she has to duplicate the line as often as they have visits.

    Example: c in the template above. The code assumes that the user sorted the rows in the correct order.

  • value

    If fixed was set to True, this will be the used value for the parameter.

    If fixed was set to False, this is the initial guess for the parameter.

  • lo_lim

    Use lower bounds for the least squares routine?

  • lo_val

    lower bound value for the least squares routine.

  • hi_lim

    Use upper bounds for the least squares routine?

  • hi_val

    upper bound value for the least squares routine.

  • prior

    Prior for the sampling?

    • X: No prior

    • U: uniform prior

    • N: Gaussian prior

  • p1 & p2

    If prior = U -> lower and upper bounds for the uniform prior

    if prior = N -> mean and 1 sigma for the gaussian prior

  • step_size

    Sets a step_size for the least squares and sampling.

Add a new model

What do I need?

  • A model name

  • The parameters of your fitting model

Which files do I have to change?

  • src/pacman/lib/formatter.py

  • src/pacman/lib/functions.py

  • src/pacman/lib/models/MODEL_NAME.py

  • src/pacman/data/fit_par.py

What do I have to do?

As an example, we will create a new model which will fit a polynomial of third order to the whole visit. We will call this model polynomial3.

  1. We create a new py file in src/pacman/lib/models/ called polynomial3.py

  2. We create the code for this model which combines a fitting for the scale and a linear polynomial over the visit:

import sys
sys.path.insert(0,'..')
import numpy as np

def polynomial3(t, data, params, visit = 0):
    v, v2, v3 = params
    v = v[visit]
    v2 = v2[visit]
    v3 = v3[visit]

    idx = data.vis_idx[visit]

    return 1. + v*data.t_vis[idx] + v2*(data.t_vis[idx])**2 + v3*(data.t_vis[idx])**3

Note that we use the same function name (def polynomial3(...)) as file name. Our parameters for this model are therefore: v, v2, and v3.

  1. Add the model parameters into src/pacman/lib/formatter.py. Each parameter needs its own row. Below an example on how to do it with our new polynomial3.py model.

        elif 'polynomial2' in data.s30_myfuncs:
            self.v = params[data.par_order['v']*data.nvisit:(1 + data.par_order['v'])*data.nvisit]
            self.v2 = params[data.par_order['v2']*data.nvisit:(1 + data.par_order['v2'])*data.nvisit]
        elif 'polynomial3' in data.s30_myfuncs:
            self.v = params[data.par_order['v']*data.nvisit:(1 + data.par_order['v'])*data.nvisit]
            self.v2 = params[data.par_order['v2']*data.nvisit:(1 + data.par_order['v2'])*data.nvisit]
            self.v3 = params[data.par_order['v3']*data.nvisit:(1 + data.par_order['v3'])*data.nvisit]
        elif 'logarithmic_visit' in data.s30_myfuncs:
            self.log1 = params[data.par_order['log1']*data.nvisit:(1 + data.par_order['log1'])*data.nvisit]
            self.log2 = params[data.par_order['log2']*data.nvisit:(1 + data.par_order['log2'])*data.nvisit]
  1. Next we add the model to src/pacman/lib/functions.py. First we have to import the function at the very top:

from ..lib.models.polynomial2 import polynomial2
from ..lib.models.polynomial3 import polynomial3
from ..lib.models.logarithmic_visit import logarithmic_visit
  1. Next, we have to add the model to the file. See an example here for our model.

            elif f == "polynomial2":
                self.sys.append(polynomial2)
                self.sys_porder.append([
                    data.par_order['v']*data.nvisit,
                    data.par_order['v2']*data.nvisit
                ])
            elif f == "polynomial3":
                self.sys.append(polynomial3)
                self.sys_porder.append([
                    data.par_order['v']*data.nvisit,
                    data.par_order['v2']*data.nvisit,
                    data.par_order['v3']*data.nvisit
                ])
            elif f == "logarithmic_visit":
                self.sys.append(logarithmic_visit)
                self.sys_porder.append([
                    data.par_order['log1']*data.nvisit,
                    data.par_order['log2']*data.nvisit
                ])
  1. Note, that the model, which we are currently adding is on to fit for systematics. Thats why we use sys and sys_porder. If you are adding a new model fitting for an astrophysical effect you use instead astro and astro_porder. See an example below for the transit model:

            elif f == "transit":
                self.astro.append(transit)
                self.astro_porder.append([
                    data.par_order['t0']*data.nvisit,
                    data.par_order['per']*data.nvisit, 
                    data.par_order['rp']*data.nvisit, 
                    data.par_order['a']*data.nvisit, 
                    data.par_order['inc']*data.nvisit, 
                    data.par_order['ecc']*data.nvisit, 
                    data.par_order['w']*data.nvisit, 
                    data.par_order['u1']*data.nvisit, 
                    data.par_order['u2']*data.nvisit
                    ])
  1. Add the new lines with your model parameters to your fit_par.txt:

#parameter     fixed   tied   value       lo_lim   lo_val    hi_lim   hi_val   prior   p1        p2       step_size  
per            True    -1     1.5804046   True     1.55      True     1.6      X       2.25314   2e-05    1e-07      
t0             False   -1     0.18        True     0.16      True     0.2      U       0.0       0.3      0.001      
t_secondary    True    -1     0.0         False    0.0       False    0.0      X       0.0       0.0      0.0        
w              True    -1     90.0        False    0.0       False    0.0      X       0.0       0.0      0.0        
a              True    -1     15.23       True     14.0      True     16.0     X       4.98      0.05     0.02       
inc            True    -1     89.1        True     88.0      True     89.9     X       85.3      0.2      0.02       
rp             False   -1     0.116       True     0.05      True     0.2      U       0.01      0.3      0.001      
fp             True    -1     0.0         False    0.0       False    0.002    X       0.0       1.0      0.0        
u1             False   -1     0.29        True     0.0       True     0.6      U       0.0       1.0      0.01       
u2             True    -1     0.0         False    0.0       False    0.0      X       0.0       0.0      0.0        
ecc            True    -1     0.0         False    0.0       False    0.0      X       0.0       0.0      0.0        
c              False   -1     8.37        True     6.0       True     8.9      U       6.7       7.0      0.001      
v              False   -1     0           True     -1        True     1.0      U       -7e-06    -1e-06   1e-05      
v2             False   -1     0.0         True     -1        True     1.0      U       0.0       0.0      1e-11      
v3             False   -1     0.0         True     -1        True     1.0      U       0.0       0.0      1e-11      
r1             False   -1     0.1         False    0.0       False    1.0      U       -10.0     10.0     0.01       
r2             False   -1     0.0         False    5.5       False    7.5      U       -100.0    100.0    0.1        
r3             True    -1     0.0         False    -1.0      False    1.0      X       -10.0     10.0     0.001      
scale          False   -1     0.0         False    -0.01     False    0.01     U       -0.1      0.1      0.002      
uncmulti_val   True    -1     1.0         False    1.0       False    7.0      X       0.1       10.0     0.01       
  1. Add the new model to the list of models in s30_myfuncs in the obs_par.pcf:

##30
s30_myfuncs                  ['constant','upstream_downstream','model_ramp','polynomial3','transit']
  1. Run s30!