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
-
A constant (free parameters: c)
Note
c is in log10. An average flux of 10^7 photoelectrons therefore leads to approximately c = 7.
-
Exponential ramps fit to each orbit (free parameters: r1, r2, r3)
-
Linear slope (free parameters: v)
-
Quadratic trend (free parameters: v, v2)
-
Exponential trend over the whole visit (free parameters: exp1, exp2)
-
Logarithmic trend over the whole visit (free parameters: log1, log2)
-
Offset between scan directions due to the upstream-downstream effect (free parameters: scale)
-
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).
-
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.
-
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
-
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
-
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.
-
Sinusoid (free parameters: a1, omega1, phi1)
-
Sum of three sinusoids (free parameters: a1, omega1, phi1, a2, omega2, phi2, a3, omega3, phi3, a12, omega12, phi12, a22, omega22, phi22, a32, omega32, phi32)
-
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
.
We create a new py file in src/pacman/lib/models/ called polynomial3.py
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
.
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]
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
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
])
Note, that the model, which we are currently adding is on to fit for systematics. Thats why we use
sys
andsys_porder
. If you are adding a new model fitting for an astrophysical effect you use insteadastro
andastro_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
])
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
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']
Run s30!