SIMBAD info: http://simbad.u-strasbg.fr/simbad/sim-id?Ident=KIC7198959
Lightcurve data from: https://archive.stsci.edu/kepler/publiclightcurves.html
In [1]:
# !curl -O http://archive.stsci.edu/pub/kepler/lightcurves/0071/007198959/kplr007198959-2009259160929_llc.fits
In [2]:
from astropy.io import fits
hdulist = fits.open('kplr007198959-2009259160929_llc.fits')
hdulist.info()
In [3]:
hdulist[1].header
Out[3]:
In [4]:
from astropy.table import Table
data = Table(hdulist[1].data)
data
Out[4]:
In [5]:
df = data.to_pandas()[['TIME', 'SAP_FLUX', 'SAP_FLUX_ERR']]
df.shape
Out[5]:
In [6]:
df = df.dropna()
df.shape
Out[6]:
In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy.stats import LombScargle
plt.style.use('seaborn-whitegrid')
In [8]:
from astropy.stats import LombScargle
In [9]:
ls = LombScargle(df['TIME'], 1, center_data=False, fit_mean=False)
freqW, powerW = ls.autopower(minimum_frequency=0,
maximum_frequency=200)
In [10]:
# Find the maximum near 2 hours
f, p = ls.autopower(minimum_frequency=1.95*24,
maximum_frequency=2.05*24,
samples_per_peak=100)
f_ny = f[np.argmax(p)]
In [11]:
t_sorted = np.sort(df['TIME'])
p_ny = 24 * 60 * 60 / f_ny
delta_t = (t_sorted[1:] - t_sorted[:-1]) * 24 * 60 * 60
In [12]:
ls = LombScargle(df['TIME'], df['SAP_FLUX'], df['SAP_FLUX_ERR'])
freq, power = ls.autopower(minimum_frequency=0,
maximum_frequency=200)
fmax = freq[np.argmax(power)] / 24
In [13]:
fig, ax = plt.subplots(2, 2, figsize=(12, 5))
fig.suptitle('Kepler object ID 007198959', size=14)
fig.subplots_adjust(hspace=0.35, wspace=0.15, left=0.07, right=0.97)
# upper left
ax[0, 0].plot(df['TIME'], df['SAP_FLUX'] / 1E6, 'ok', markersize=2, rasterized=True)
ax[0, 0].set(ylabel='SAP flux ($10^6 e^-/s$)',
title='Observed light curve',
xlim=(168, 260))
# bottom left
left, bottom, width, height = ax[1, 0].get_position().bounds
ax[1, 0].set_position([left, bottom + 0.15, width, height-0.15])
inset = fig.add_axes([left, bottom, width, 0.1])
ax[1, 0].plot(t_sorted[:-1], delta_t / 60, 'ok', markersize=2, rasterized=True)
ax[1, 0].axhline(p_ny / 60, color='gray', linestyle='--')
ax[1, 0].set(xlim=ax[0, 0].get_xlim(),
ylim=(10, 10000),
yscale='log',
ylabel='$\Delta t$ (min)',
title='Time between observations')
ax[1, 0].xaxis.set_major_formatter(plt.NullFormatter())
inset.plot(t_sorted[:-1], 1000 * (delta_t - p_ny), 'ok', markersize=2, rasterized=True)
inset.axhline(0, color='gray', linestyle='--')
inset.set(xlim=ax[0, 0].get_xlim(),
ylim=(-100, 100),
xlabel='Observation time (days)',
ylabel='$\Delta t - p_{W}$ (ms)')
inset.yaxis.set_major_locator(plt.MaxNLocator(3));
# Upper right
ax[0, 1].plot(freqW / 24, powerW, '-k', rasterized=True);
ax[0, 1].set(xlim=(0, 6.5),
ylim=(-0.1, 1.1),
ylabel='Lomb-Scargle power',
title='Window Power Spectrum');
ax[0, 1].annotate('', (0, 0.6), (f_ny / 24, 0.6),
arrowprops=dict(arrowstyle='<->', color='gray'));
ax[0, 1].text(f_ny / 48, 0.6, r'$({0:.1f}\ {{\rm minutes}})^{{-1}}$'.format(24 * 60 / f_ny),
size=12, ha='center', va='bottom');
# Lower right
ax[1, 1].plot(freq / 24, power, '-k', rasterized=True)
ax[1, 1].fill_between([0.5 * f_ny / 24, 1.5 * f_ny / 24], -0.05, 1,
color='gray', alpha=0.3)
ax[1, 1].text(2.25, 0.95, r"(Aliased Region)", size=14, color='gray', ha='right', va='top')
ax[1, 1].text(fmax, 0.85, r"$f_{{max}}=({0:.2f}\ {{\rm hours}})^{{-1}}$".format(1 / fmax),
size=12)
ax[1, 1].set(xlim=(0, 2.3),
ylim=(-0.05, 1.0),
xlabel='frequency (hours$^{-1}$)',
ylabel='Lomb-Scargle power',
title='Light Curve Power Spectrum');
fig.savefig('fig16_kepler_data.pdf')
In [14]:
n_o = 5
T = df['TIME'].max() - df['TIME'].min()
delta_f = 1. / n_o / T
print("f_ny =", f_ny)
print("T =", T)
print("n_grid =", f_ny / delta_f)