# Eclipse of WASP-189 b¶

## Introduction¶

Warning

This tutorial does **not** use real CHEOPS data, only a realistic example
data set. Don’t interpret any science results from this example!

The strength of some CHEOPS observations will come from the spacecraft’s ability to return to targets and repeatedly measure low signal-to-noise events until we can confidently detect small photometric features, like the depth of secondary eclipse of a hot exoplanet. This is the case for WASP-189 b, a very-hot Jupiter which orbits a bright A star (Anderson et al. 2018).

## Sigma clipping¶

First, let’s import some packages we’ll need, including `linea`

:

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binned_statistic
from batman import TransitModel
from linea import JointLightCurve, Planet
# Load the transit parameters for WASP-189 b
p = Planet.from_name('WASP-189 b')
# Load the example light curve of WASP-189 b (built into the package)
lcs = JointLightCurve.from_example()
```

In the last line, we initialized a `JointLightCurve`

object using the
`from_example`

method, which loads example (fake) data
of four visits of WASP-189 b.

Then we sigma-clip the light curve – this removes outliers.

```
# Sigma clip the light curve to remove outliers
for lc in lcs:
lc.sigma_clip_flux(sigma_upper=3, sigma_lower=3, plot=True)
plt.show()
```

(Source code, png, hires.png, pdf)

## Regression analysis¶

We assemble each indidiual visit’s design matrix by calling
`combined_design_matrix`

, which combines them into one
big design matrix. We’ll also concatenate the design matrix of basis vectors
with an eclipse model, computed by `batman`

:

```
all_lcs = lcs.concatenate()
eclipse_model = TransitModel(p, all_lcs.bjd_time[~all_lcs.mask],
supersample_factor=3,
transittype='secondary',
exp_time=all_lcs.bjd_time[1]-all_lcs.bjd_time[0]
).light_curve(p) - 1
X = np.hstack([
lcs.combined_design_matrix(),
eclipse_model[:, None]
])
```

Now `X`

contains our “final” design matrix, consisting of the all of the
basis vectors we can think of, and the eclipse model.

To do the linear regression, simply call the `regress`

method:

```
r = lcs.regress(X)
```

The solution to the linear regression is stored in `r`

. One neat measurement
we can pull directly from the `r`

object is the best-fit eclipse depth and
its uncertainty. The transit model we intiailized in `eclipse_model`

has a
depth of 1 ppm, so the least squares weight for the eclipse (last) basis vector
is the amplitude of the secondary eclipse in units of ppm. The eclipse depth
and uncertainty are:

```
>>> print(f"Eclipse depth = {r.betas[-1] :.0f} ± {np.sqrt(np.diag(r.cov))[-1] :.0f} ppm")
Eclipse depth = 77 ± 7 ppm
```

Finally, let’s plot the best-fit detrended light curve and eclipse model:

```
# Compute orbital phase for every time
phases = np.concatenate([lc.phase(p)[~lc.mask] for lc in lcs])
sort = np.argsort(phases)
# Compute the best-fit systematics model, without removing the eclipse
obs_eclipse = all_lcs.flux[~all_lcs.mask] / (X[:, :-1] @ r.betas[:-1])
# Compute the best-fit eclipse model
eclipse_model = X[:, -1] * r.betas[-1]
# Binned light curve:
bs = binned_statistic(phases[sort], obs_eclipse[sort], bins=30)
bincenters = 0.5 * (bs.bin_edges[1:] + bs.bin_edges[:-1])
# Create plot:
fig, ax = plt.subplots(figsize=(4, 3), sharex=True)
ax.plot(phases, obs_eclipse, '.', color='silver')
ax.plot(phases[sort], eclipse_model[sort] + 1, 'r', lw=2)
ax.plot(bincenters, bs.statistic, 's', color='k')
for sp in ['right', 'top']:
ax.spines[sp].set_visible(False)
ax.ticklabel_format(useOffset=False)
ax.set(xlabel='Phase', ylabel='Flux',
ylim=[0.99965, 1.0004])
plt.show()
```

(Source code, png, hires.png, pdf)

Note

The above plot is a simulated example light curve, **not** real
CHEOPS observations. Do not make any conclusions about the planet from
this fake dataset.

We can clearly see the ~80 ppm secondary eclipse which occurs when the planet is occulted by the star.