SIP demo.¶

This demo shows how to compute an SIP for your favourite K2 target.

In [38]:

import numpy as np
import matplotlib.pyplot as plt
import wget
import h5py
import fitsio
from SIP import SIP, eval_freq
%matplotlib inline

Download the Campaign 1 Eigen Light Curves (ELCs) and a light curve for example star, EPIC 201133037.

In [3]:

wget.download("http://bbq.dfm.io/ketu/elcs/c1.h5")
wget.download("http://bbq.dfm.io/ketu/lightcurves/c1/201100000/33000/ktwo201133037-c01_lpd-lc.fits")

Out[3]:

'ktwo201133037-c01_lpd-lc.fits'

Load the light curve.

In [19]:

data = fitsio.read("ktwo201133037-c01_lpd-lc.fits")
aps = fitsio.read("ktwo201133037-c01_lpd-lc.fits", 2)
y = data["flux"][:, np.argmin(aps["cdpp6"])]
x = data["time"]
q = data["quality"]
m = np.isfinite(y) * np.isfinite(x) * (q==0)  # remove nans and bad data points
y, x = y[m], x[m]
y = y / np.median(y) - 1 # median normalise

Load the top 150 ELCs.

In [20]:

nELC = 150
with h5py.File("c1.h5", "r") as f:
    basis = f["basis"][:nELC, m]

Plot the raw light curve.

In [21]:

plt.plot(x, y, "k")
plt.xlabel("BJD - 2454833 (days)")
plt.ylabel("Normalised Flux")

Out[21]:

<matplotlib.text.Text at 0x10a58fed0>

Compute the SIP.

In [22]:

periods = np.arange(1., 70., .1)  # test periods ranging from 1 to 70 days
freqs = 1./periods
s2n, amp2s, w = SIP(x, y, basis, freqs)

plt.plot(periods, s2n, "k")
plt.xlabel("Period (days)")
plt.ylabel("Relative (S/N)^2")

Out[22]:

<matplotlib.text.Text at 0x10aad3690>

Now to compute the conditioned light curve. Firstly, locate the highest peak in the SIP.

In [25]:

peaks = np.array([i for i in range(1, len(periods)-1) if s2n[i-1] < s2n[i] and s2n[i+1] < s2n[i]])
l = s2n[peaks] == max(s2n[peaks])
Pmax = periods[peaks][l][0] 
print "Pmax = ", Pmax, "days"

Pmax =  24.8 days

Construct arrays.

In [39]:

AT = np.concatenate((basis, np.ones((3, len(y)))), axis=0)
ATA = np.dot(AT, AT.T)

Compute trends.

In [40]:

_, _, trends = eval_freq(x, y, 1./Pmax, AT, ATA, compute_trends=True)

In [43]:

plt.plot(x, y, ".5", label="raw")
plt.plot(x, y-trends, "k", label="conditioned")
plt.xlabel("BJD - 2454833 (days)")
plt.ylabel("Normalised Flux")
plt.legend(loc="best")

Out[43]:

<matplotlib.legend.Legend at 0x10aa142d0>

In [ ]: