# 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.

### 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
])
```

```#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!