``````

In [1]:

from __future__ import division, print_function

%matplotlib inline
%config InlineBackend.figure_format = "retina"

from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["font.size"] = 20

``````
``````

In [2]:

import h5py
import fitsio
import numpy as np
import matplotlib.pyplot as pl
from scipy.linalg import cho_solve, cho_factor

from astropy.io import fits
from astropy.wcs import WCS

``````

## Target pixel files

First we'll read in the data:

``````

In [3]:

``````

Let's look at one of the images (the last one in the time series). It's in the `FLUX` column of the TPF file:

``````

In [4]:

img = tpf["FLUX"][-1]
pl.imshow(img.T, cmap="gray", interpolation="nearest");

``````
``````

``````

## Centroiding

Here's some code for measuring the centroid using the 3x3 patch and modeling it as a 2D quadratic. I'll also calcuate the second derivatives of the flux in the 3x3 patch because this should be related to the shape of the PSF.

``````

In [5]:

# These are some useful things to pre-compute and use later.
_x, _y = np.meshgrid(range(-1, 2), range(-1, 2), indexing="ij")
_x, _y = _x.flatten(), _y.flatten()
_AT = np.vstack((_x*_x, _y*_y, _x*_y, _x, _y, np.ones_like(_x)))
_ATA = np.dot(_AT, _AT.T)
factor = cho_factor(_ATA, overwrite_a=True)

# This function finds the centroid and second derivatives in a 3x3 patch.
def fit_3x3(img):
a, b, c, d, e, f = cho_solve(factor, np.dot(_AT, img.flatten()))
m = 1. / (4 * a * b - c*c)
x = (c * e - 2 * b * d) * m
y = (c * d - 2 * a * e) * m
dx2, dy2, dxdy = 2 * a, 2 * b, c
return [x, y, dx2, dy2, dxdy]

# This function finds the centroid in an image.
# You can provide an estimate of the centroid using WCS.
def find_centroid(img, init=None):
if init is None:
xi, yi = np.unravel_index(np.argmax(img), img.shape)
else:
xi, yi = map(int, map(np.round, init))
ox, oy = np.unravel_index(np.argmax(img[xi-1:xi+2, yi-1:yi+2]), (3, 3))
xi += ox - 1
yi += oy - 1
assert (xi >= 1 and xi < img.shape[0] - 1), "effed, x"
assert (yi >= 1 and yi < img.shape[1] - 1), "effed, y"
pos = fit_3x3(img[xi-1:xi+2, yi-1:yi+2])
pos[0] += xi
pos[1] += yi
return pos

``````

We'll use astropy to get the WCS and estimate the location of the star.

``````

In [6]:

with fits.open("ktwo201466124-c01_lpd-targ.fits.gz") as hdus:
wcs = WCS(hdr)
# The order is the opposite of what I normally use...
init = wcs.wcs_world2pix(hdr["RA_OBJ"], hdr["DEC_OBJ"], 0.0)[::-1]

``````

Plot the initial centroid estimate:

``````

In [7]:

img = tpf["FLUX"][-1]
pl.imshow(img.T, cmap="gray", interpolation="nearest")
pl.plot(init[0], init[1], "or");

``````
``````

``````

And then refine it using the code above:

``````

In [8]:

cx, cy = find_centroid(img, init)[:2]

``````
``````

In [9]:

img = tpf["FLUX"][-1]
pl.imshow(img.T, cmap="gray", interpolation="nearest")
pl.plot(init[0], init[1], "or")
pl.plot(cx, cy, ".b");

``````
``````

``````

Now let's measure it at each time by looping over frames.

``````

In [10]:

coords = np.nan + np.zeros((len(tpf), 5))
for i, img in enumerate(tpf["FLUX"]):
msk = np.isfinite(img)
if np.any(msk) and np.any(img[msk] > 0.0):
coords[i] = find_centroid(img, init=init)

``````
``````

In [11]:

print(coords.shape)

``````
``````

(4022, 5)

``````

The elements of `coords` at each time are `x, y, d^2f/dx^2, d^2f/dy^2, d^2f/dxdy`.

## Aperture photometry

Let's just use the photometry that I measured already...

``````

In [12]:

``````

Then we'll choose the "best" aperture size, mask the missing data and normalize the fluxes:

``````

In [13]:

time = data["time"]
flux = data["flux"][:, np.argmin(aps["cdpp6"])]  # Lowest CDPP6.
qual = data["quality"]

good_mask = np.isfinite(time) & np.isfinite(flux) & (qual == 0)

``````
``````

In [15]:

pl.plot(good_time, good_flux, "k")
pl.xlim(1980, 1990);

``````
``````

``````

## Going forward...

I think that what we want to do is build a non-linear model of `coords` (the centroids that we measured) that best predicts `good_flux`. Crossfield et al. use a Gaussian Process and I think that that would be a good idea but there are lots of other things to try.

Probably the best first bet is try plotting all the things against each other and see where we get...

``````

In [16]:

pl.xlabel("x [pix]")
pl.ylabel("y [pix]");

``````
``````

``````
``````

In [17]:

import george
from george.kernels import Matern32Kernel, ExpSquaredKernel

``````
``````

In [18]:

X = np.concatenate((coords[good_mask, :2], np.atleast_2d(good_time).T), axis=1)
X -= np.mean(X, axis=0)[None, :]
X /= np.std(X, axis=0)[None, :]

``````
``````

In [36]:

d = 1000.0 * np.ones(X.shape[1])
d[-1] = 100.0
kernel = np.var(good_flux) * ExpSquaredKernel(d, ndim=X.shape[1])
gp = george.GP(kernel, mean=1.0)

``````
``````

In [37]:

gp.compute(X, 1e-8)

``````
``````

In [38]:

mu = gp.predict(good_flux, X, mean_only=True)

``````
``````

In [39]:

mu.shape

``````
``````

Out[39]:

(3631,)

``````
``````

In [40]:

pl.plot(good_time, good_flux, "k")
pl.plot(good_time, mu, "r")
# pl.xlim(2040, 2050);

``````
``````

``````
``````

In [43]:

pl.plot(good_time, good_flux - mu, "k");
# pl.xlim(2040, 2050);

``````
``````

``````
``````

In [131]:

np.sum((X[500] - X[1000]) ** 2)

``````
``````

Out[131]:

4.1494646870178595

``````
``````

In [ ]:

``````