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)
For models using
constant, PACMAN checks whether the initial value and prior bounds ofcare plausible for the light curve being fitted. Theconstantmodel usesflux = 10**c, socshould be close tolog10of the raw flux level. This helps catch cases where a white-light normalization is accidentally used for spectroscopic light curves, or vice versa.
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.
-
Linear trend with shift of the spectrum on the detector in spatial direction (rows). Independent application for both scan directions. Implemented similarly to the method explained in Tsiaras et al. (2016). (free parameters: vf, vr)
Note
Requires that stage 20 was executed with calculate_rowshift True.
-
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 prior p1 p2
per True -1 1.5804046 X 2.25314 2e-05
t0 False -1 0.18 U 0.18 0.19
t_secondary True -1 0.0 X 0.0 0.0
w True -1 90.0 X 0.0 0.0
a True -1 15.23 X 4.98 0.05
inc True -1 89.1 X 85.3 0.2
rp False -1 0.116 U 0.01 0.2
fp True -1 0.0 X 0.0 1.0
u1 False -1 0.29 U 0.0 1.0
u2 True -1 0.0 X 0.0 0.0
ecc True -1 0.0 X 0.0 0.0
c False 0 8.45 U 8.2 8.5
c False 1 8.45 U 8.2 8.5
v False 0 0 U -1e-05 1e-05
v False 1 0 U -1e-05 1e-05
v2 True -1 0.0 X 0.0 0.0
r1 False 0 0.1 U 0 1.0
r1 False 1 0.1 U 0 1.0
r2 False 0 0.1 U 0.0 0.1
r2 False 1 0.1 U 0.0 0.1
r3 True -1 0.0 U -10.0 10.0
scale False 0 0.0 U -0.1 0.1
scale False 1 0.0 U -0.1 0.1
uncmulti_val False 0 2.0 U 0.1 8.0
uncmulti_val False 1 2.0 U 0.1 8.0
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, they have 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.
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
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
sysandsys_porder. If you are adding a new model fitting for an astrophysical effect you use insteadastroandastro_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 prior p1 p2
per True -1 1.5804046 X 2.25314 2e-05
t0 False -1 0.18 U 0.18 0.19
t_secondary True -1 0.0 X 0.0 0.0
w True -1 90.0 X 0.0 0.0
a True -1 15.23 X 4.98 0.05
inc True -1 89.1 X 85.3 0.2
rp False -1 0.116 U 0.01 0.2
fp True -1 0.0 X 0.0 1.0
u1 False -1 0.29 U 0.0 1.0
u2 True -1 0.0 X 0.0 0.0
ecc True -1 0.0 X 0.0 0.0
c False 0 8.45 U 8.2 8.5
c False 1 8.45 U 8.2 8.5
v False 0 0 U -1e-05 1e-05
v False 1 0 U -1e-05 1e-05
rowshift_vr True -1 0.0 X -1.0 1.0
rowshift_vf True -1 0.0 X -1.0 1.0
v2 True -1 0.0 X 0.0 0.0
v3 False -1 0.0 U -1.0 1.0
r1 False 0 0.1 U 0 1.0
r1 False 1 0.1 U 0 1.0
r2 False 0 0.1 U 0.0 0.1
r2 False 1 0.1 U 0.0 0.1
r3 True -1 0.0 U -10.0 10.0
scale False 0 0.0 U -0.1 0.1
scale False 1 0.0 U -0.1 0.1
uncmulti_val False 0 2.0 U 0.1 8.0
uncmulti_val False 1 2.0 U 0.1 8.0
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!