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]:
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]:
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]:
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"
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]:
In [ ]: