In this notebook we will replicate the K2SFF method from Vanderburg and Johnson 2014. The paper introduces a method for "Self Flat Fielding", by tracking how the lightcurve changes with motion of the spacecraft.
from pyke import KeplerTargetPixelFile
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
#from oktopus import UniformPrior, JointPrior
#from pyke import PRFPhotometry, SceneModel
#from pyke.utils import KeplerQualityFlags
#! wget https://www.cfa.harvard.edu/~avanderb/k2/ep60021426alldiagnostics.csv
import pandas as pd
df = pd.read_csv('ep60021426alldiagnostics.csv',index_col=False)
df.head()
BJD - 2454833 | Raw Flux | Corrected Flux | X-centroid | Y-centroid | arclength | Correction | Thrusters On | |
---|---|---|---|---|---|---|---|---|
0 | 1862.502368 | 0.995119 | 0.995985 | 25.135097 | 24.661074 | 2.327480 | 0.999130 | 1.0 |
1 | 1862.522801 | 0.997313 | 0.996767 | 25.289752 | 24.418689 | 1.175322 | 1.000548 | 1.0 |
2 | 1862.543235 | 0.996713 | 0.996136 | 25.288052 | 24.429406 | 1.214627 | 1.000580 | 0.0 |
3 | 1862.563668 | 0.996930 | 0.996277 | 25.275216 | 24.448405 | 1.306617 | 1.000656 | 0.0 |
4 | 1862.584102 | 0.996862 | 0.996228 | 25.253864 | 24.480184 | 1.460259 | 1.000636 | 0.0 |
Let's use the provided $x-y$ centroids, but we could compute these on our own too.
col = df[' X-centroid'].values
col = col - np.mean(col)
row = df[' Y-centroid'].values
row = row - np.mean(row)
def _get_eigen_vectors(centroid_col, centroid_row):
centroids = np.array([centroid_col, centroid_row])
eig_val, eig_vec = np.linalg.eigh(np.cov(centroids))
return eig_val, eig_vec
def _rotate(eig_vec, centroid_col, centroid_row):
centroids = np.array([centroid_col, centroid_row])
return np.dot(eig_vec, centroids)
eig_val, eig_vec = _get_eigen_vectors(col, row)
v1, v2 = eig_vec
The major axis is the last one.
plt.figure(figsize=(5, 6))
plt.plot(col*4.0, row*4.0, 'ko', ms=4)
plt.plot(col*4.0, row*4.0, 'ro', ms=1)
plt.xticks([-2, -1,0, 1, 2])
plt.yticks([-2, -1,0, 1, 2])
plt.xlabel('X position [arcseconds]')
plt.ylabel('Y position [arcseconds]')
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.plot([0, v1[0]], [0, v1[1]], color='blue', lw=3)
plt.plot([0, v2[0]], [0, v2[1]], color='blue', lw=3);
Following the form of Figure 2 of Vanderburg & Johsnon 2014.
rot_colp, rot_rowp = _rotate(eig_vec, col, row)
You can rotate into the new reference frame.
plt.figure(figsize=(5, 6))
plt.plot(rot_rowp*4.0, rot_colp*4.0, 'ko', ms=4)
plt.plot(rot_rowp*4.0, rot_colp*4.0, 'ro', ms=1)
plt.xticks([-2, -1,0, 1, 2])
plt.yticks([-2, -1,0, 1, 2])
plt.xlabel("X' position [arcseconds]")
plt.ylabel("Y' position [arcseconds]")
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.plot([0, 1], [0, 0], color='blue')
plt.plot([0, 0], [0, 1], color='blue');
We need to calculate the arclength using: $$s= \int_{x'_0}^{x'_1}\sqrt{1+\left( \frac{dy'_p}{dx'}\right)^2} dx'$$
where $x^\prime_0$ is the transformed $x$ coordinate of the point with the smallest $x^\prime$ position, and $y^\prime_p$ is the best--fit polynomial function.
z = np.polyfit(rot_rowp, rot_colp, 5)
p5 = np.poly1d(z)
p5_deriv = p5.deriv()
x0_prime = np.min(rot_rowp)
xmax_prime = np.max(rot_rowp)
x_dense = np.linspace(x0_prime, xmax_prime, 2000)
plt.plot(rot_rowp, rot_colp, '.')
plt.plot(x_dense, p5(x_dense));
@np.vectorize
def arclength(x):
'''Input x1_prime, get out arclength'''
gi = x_dense <x
s_integrand = np.sqrt(1 + p5_deriv(x_dense[gi]) ** 2)
s = np.trapz(s_integrand, x=x_dense[gi])
return s
plt.plot(df[' arclength'], arclength(rot_rowp)*4.0, '.')
plt.plot([0, 4], [0, 4], 'k--');
It works!
Now we apply a high-pass filter. We follow the original paper by using BSplines with 1.5 day breakpoints.
from scipy.interpolate import BSpline
from scipy import interpolate
tt, ff = df['BJD - 2454833'].values, df[' Raw Flux'].values
tt = tt - tt[0]
knots = np.arange(0, tt[-1], 1.5)
t,c,k = interpolate.splrep(tt, ff, s=0, task=-1, t=knots[1:])
bspl = BSpline(t,c,k)
plt.plot(tt, ff, '.')
plt.plot(tt, bspl(tt))
[<matplotlib.lines.Line2D at 0xa1fddacf8>]
Spline fit looks good, so normalize the flux by the long-term trend.
Plot the normalized flux versus arclength to see the position-dependent flux.
norm_ff = ff/bspl(tt)
Mask the data by keeping only the good samples.
bi = df[' Thrusters On'].values == 1.0
gi = df[' Thrusters On'].values == 0.0
al, gff = arclength(rot_rowp[gi])*4.0, norm_ff[gi]
sorted_inds = np.argsort(al)
We will follow the paper by interpolating 15 bins of means. This is a piecewise linear fit.
knots = np.array([np.min(al)]+
[np.median(splt) for splt in np.array_split(al[sorted_inds], 15)]+
[np.max(al)])
bin_means = np.array([gff[sorted_inds][0]]+
[np.mean(splt) for splt in np.array_split(gff[sorted_inds], 15)]+
[gff[sorted_inds][-1]])
zz = np.polyfit(al, gff,6)
sff = np.poly1d(zz)
al_dense = np.linspace(0, 4, 1000)
interp_func = interpolate.interp1d(knots, bin_means)
plt.figure(figsize=(5, 6))
plt.plot(arclength(rot_rowp)*4.0, norm_ff, 'ko', ms=4)
plt.plot(arclength(rot_rowp)*4.0, norm_ff, 'o', color='#3498db', ms=3)
plt.plot(arclength(rot_rowp[bi])*4.0, norm_ff[bi], 'o', color='r', ms=3)
#plt.plot(al_dense, sff(al_dense), '-', color='#e67e22')
#plt.plot(knots, bin_means, '-', color='#e67e22')
plt.plot(np.sort(al), interp_func(np.sort(al)), '-', color='#e67e22')
plt.xticks([0, 1,2, 3, 4])
plt.xlabel('Arclength [arcseconds]')
plt.ylabel('Relative Brightness')
plt.title('EPIC 60021426, Kp =10.3')
plt.xlim(0,4)
plt.ylim(0.997, 1.002);
Following Figure 4 of Vanderburg & Johnson 2014.
Apply the Self Flat Field (SFF) correction:
corr_flux = gff / interp_func(al)
plt.figure(figsize=(10,6))
dy = 0.004
plt.plot(df['BJD - 2454833'], df[' Raw Flux']+dy, 'ko', ms=4)
plt.plot(df['BJD - 2454833'], df[' Raw Flux']+dy, 'o', color='#3498db', ms=3)
plt.plot(df['BJD - 2454833'][bi], df[' Raw Flux'][bi]+dy, 'o', color='r', ms=3)
plt.plot(df['BJD - 2454833'][gi], corr_flux*bspl(tt[gi]), 'o', color='k', ms = 4)
plt.plot(df['BJD - 2454833'][gi], corr_flux*bspl(tt[gi]), 'o', color='#e67e22', ms = 3)
#plt.plot(df['BJD - 2454833'][gi], df[' Corrected Flux'][gi], 'o', color='#00ff00', ms = 4)
plt.xlabel('BJD - 2454833')
plt.ylabel('Relative Brightness')
plt.xlim(1862, 1870)
plt.ylim(0.994, 1.008);
Following Figure 5 of Vanderburg & Johnson 2015.
Let's compute the CDPP:
from pyke import LightCurve
#lc = LightCurve(time=df['BJD - 2454833'][gi], flux=corr_flux*bspl(tt[gi]))
lc = LightCurve(time=df['BJD - 2454833'][gi], flux=df[' Corrected Flux'][gi])
lc.cdpp(savgol_window=201)*1.4
38.872509891106382
The end.