Kepler Data¶

The Data¶

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()

Filename: kplr007198959-2009259160929_llc.fits
No.    Name         Type      Cards   Dimensions   Format
  0  PRIMARY     PrimaryHDU      58   ()      
  1  LIGHTCURVE  BinTableHDU    155   4354R x 20C   [D, E, J, E, E, E, E, E, E, J, D, E, D, E, D, E, D, E, E, E]   
  2  APERTURE    ImageHDU        48   (12, 41)   int32

In [3]:

hdulist[1].header

Out[3]:

XTENSION= 'BINTABLE'           / marks the beginning of a new HDU               
BITPIX  =                    8 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  100 / length of first array dimension                
NAXIS2  =                 4354 / length of second array dimension               
PCOUNT  =                    0 / group parameter count (not used)               
GCOUNT  =                    1 / group count (not used)                         
TFIELDS =                   20 / number of table fields                         
TTYPE1  = 'TIME    '           / column title: data time stamps                 
TFORM1  = 'D       '           / column format: 64-bit floating point           
TUNIT1  = 'BJD - 2454833'      / column units: barycenter corrected JD          
TDISP1  = 'D14.7   '           / column display format                          
TTYPE2  = 'TIMECORR'           / column title: barycenter - timeslice correction
TFORM2  = 'E       '           / column format: 32-bit floating point           
TUNIT2  = 'd       '           / column units: day                              
TDISP2  = 'E13.6   '           / column display format                          
TTYPE3  = 'CADENCENO'          / column title: unique cadence number            
TFORM3  = 'J       '           / column format: signed 32-bit integer           
TDISP3  = 'I10     '           / column display format                          
TTYPE4  = 'SAP_FLUX'           / column title: aperture photometry flux         
TFORM4  = 'E       '           / column format: 32-bit floating point           
TUNIT4  = 'e-/s    '           / column units: electrons per second             
TDISP4  = 'E14.7   '           / column display format                          
TTYPE5  = 'SAP_FLUX_ERR'       / column title: aperture phot. flux error        
TFORM5  = 'E       '           / column format: 32-bit floating point           
TUNIT5  = 'e-/s    '           / column units: electrons per second (1-sigma)   
TDISP5  = 'E14.7   '           / column display format                          
TTYPE6  = 'SAP_BKG '           / column title: aperture phot. background flux   
TFORM6  = 'E       '           / column format: 32-bit floating point           
TUNIT6  = 'e-/s    '           / column units: electrons per second             
TDISP6  = 'E14.7   '           / column display format                          
TTYPE7  = 'SAP_BKG_ERR'        / column title: ap. phot. background flux error  
TFORM7  = 'E       '           / column format: 32-bit floating point           
TUNIT7  = 'e-/s    '           / column units: electrons per second (1-sigma)   
TDISP7  = 'E14.7   '           / column display format                          
TTYPE8  = 'PDCSAP_FLUX'        / column title: aperture phot. PDC flux          
TFORM8  = 'E       '           / column format: 32-bit floating point           
TUNIT8  = 'e-/s    '           / column units: electrons per second             
TDISP8  = 'E14.7   '           / column display format                          
TTYPE9  = 'PDCSAP_FLUX_ERR'    / column title: ap. phot. PDC flux error         
TFORM9  = 'E       '           / column format: 32-bit floating point           
TUNIT9  = 'e-/s    '           / column units: electrons per second (1-sigma)   
TDISP9  = 'E14.7   '           / column display format                          
TTYPE10 = 'SAP_QUALITY'        / column title: aperture photometry quality flag 
TFORM10 = 'J       '           / column format: signed 32-bit integer           
TDISP10 = 'B16.16  '           / column display format                          
TTYPE11 = 'PSF_CENTR1'         / column title: PSF-fitted column centroid       
TFORM11 = 'D       '           / column format: 64-bit floating point           
TUNIT11 = 'pixel   '           / column units: pixel                            
TDISP11 = 'F10.5   '           / column display format                          
TTYPE12 = 'PSF_CENTR1_ERR'     / column title: PSF-fitted column error          
TFORM12 = 'E       '           / column format: 32-bit floating point           
TUNIT12 = 'pixel   '           / column units: pixel (1-sigma)                  
TDISP12 = 'E14.7   '           / column display format                          
TTYPE13 = 'PSF_CENTR2'         / column title: PSF-fitted row centroid          
TFORM13 = 'D       '           / column format: 64-bit floating point           
TUNIT13 = 'pixel   '           / column units: pixel                            
TDISP13 = 'F10.5   '           / column display format                          
TTYPE14 = 'PSF_CENTR2_ERR'     / column title: PSF-fitted row error             
TFORM14 = 'E       '           / column format: 32-bit floating point           
TUNIT14 = 'pixel   '           / column units: pixel (1-sigma)                  
TDISP14 = 'E14.7   '           / column display format                          
TTYPE15 = 'MOM_CENTR1'         / column title: moment-derived column centroid   
TFORM15 = 'D       '           / column format: 64-bit floating point           
TUNIT15 = 'pixel   '           / column units: pixel                            
TDISP15 = 'F10.5   '           / column display format                          
TTYPE16 = 'MOM_CENTR1_ERR'     / column title: moment-derived column error      
TFORM16 = 'E       '           / column format: 32-bit floating point           
TUNIT16 = 'pixel   '           / column units: pixel (1-sigma)                  
TDISP16 = 'E14.7   '           / column display format                          
TTYPE17 = 'MOM_CENTR2'         / column title: moment-derived row centroid      
TFORM17 = 'D       '           / column format: 64-bit floating point           
TUNIT17 = 'pixel   '           / column units: pixel                            
TDISP17 = 'F10.5   '           / column display format                          
TTYPE18 = 'MOM_CENTR2_ERR'     / column title: moment-derived row error         
TFORM18 = 'E       '           / column format: 32-bit floating point           
TUNIT18 = 'pixel   '           / column units: pixel (1-sigma)                  
TDISP18 = 'E14.7   '           / column display format                          
TTYPE19 = 'POS_CORR1'          / column title: column position correction       
TFORM19 = 'E       '           / column format: 32-bit floating point           
TUNIT19 = 'pixels  '           / column units: pixel                            
TDISP19 = 'E14.7   '           / column display format                          
TTYPE20 = 'POS_CORR2'          / column title: row position correction          
TFORM20 = 'E       '           / column format: 32-bit floating point           
TUNIT20 = 'pixels  '           / column units: pixel                            
TDISP20 = 'E14.7   '           / column display format                          
INHERIT =                    T / inherit the primary header                     
EXTNAME = 'LIGHTCURVE'         / name of extension                              
EXTVER  =                    1 / extension version number (not format version)  
TELESCOP= 'Kepler  '           / telescope                                      
INSTRUME= 'Kepler Photometer'  / detector type                                  
OBJECT  = 'KIC 7198959'        / string version of target id                    
KEPLERID=              7198959 / unique Kepler target identifier                
RADESYS = 'ICRS    '           / reference frame of celestial coordinates       
RA_OBJ  =           291.366301 / [deg] right ascension                          
DEC_OBJ =            42.784367 / [deg] declination                              
EQUINOX =               2000.0 / equinox of celestial coordinate system         
EXPOSURE=          81.90687481 / [d] time on source                             
TIMEREF = 'SOLARSYSTEM'        / barycentric correction applied to times        
TASSIGN = 'SPACECRAFT'         / where time is assigned                         
TIMESYS = 'TDB     '           / time system is barycentric JD                  
BJDREFI =              2454833 / integer part of BJD reference date             
BJDREFF =           0.00000000 / fraction of the day in BJD reference date      
TIMEUNIT= 'd       '           / time unit for TIME, TSTART and TSTOP           
TELAPSE =          88.96781229 / [d] TSTOP - TSTART                             
LIVETIME=          81.90687481 / [d] TELAPSE multiplied by DEADC                
TSTART  =         169.51002764 / observation start time in BJD-BJDREF           
TSTOP   =         258.47783993 / observation stop time in BJD-BJDREF            
LC_START=       55002.01747542 / mid point of first cadence in MJD              
LC_END  =       55090.96492052 / mid point of last cadence in MJD               
DEADC   =           0.92063492 / deadtime correction                            
TIMEPIXR=                  0.5 / bin time beginning=0 middle=0.5 end=1          
TIERRELA=             5.78E-07 / [d] relative time error                        
TIERABSO=                      / [d] absolute time error                        
INT_TIME=       6.019802903270 / [s] photon accumulation time per frame         
READTIME=       0.518948526144 / [s] readout time per frame                     
FRAMETIM=       6.538751429414 / [s] frame time (INT_TIME + READTIME)           
NUM_FRM =                  270 / number of frames per time stamp                
TIMEDEL =     0.02043359821692 / [d] time resolution of data                    
DATE-OBS= '2009-06-20T00:10:27.146Z' / TSTART as UTC calendar date              
DATE-END= '2009-09-16T23:24:11.862Z' / TSTOP as UTC calendar date               
BACKAPP =                    T / background is subtracted                       
DEADAPP =                    T / deadtime applied                               
VIGNAPP =                    T / vignetting or collimator correction applied    
GAIN    =                111.9 / [electrons/count] channel gain                 
READNOIS=            80.176350 / [electrons] read noise                         
NREADOUT=                  270 / number of read per cadence                     
TIMSLICE=                    2 / time-slice readout sequence section            
MEANBLCK=                  728 / [count] FSW mean black level                   
LCFXDOFF=               419400 / long cadence fixed offset                      
SCFXDOFF=               219400 / short cadence fixed offset                     
CDPP3_0 =     4018.13134765625 / [ppm] RMS CDPP on 3.0-hr time scales           
CDPP6_0 =    3200.025146484375 / [ppm] RMS CDPP on 6.0-hr time scales           
CDPP12_0=     2240.04345703125 / [ppm] RMS CDPP on 12.0-hr time scales          
CROWDSAP=               0.9991 / Ratio of target flux to total flux in op. ap.  
FLFRCSAP=               0.9981 / Frac. of target flux w/in the op. aperture     
NSPSDDET=                    0 / Number of SPSDs detected                       
NSPSDCOR=                    0 / Number of SPSDs corrected                      
PDCVAR  =      8497.4873046875 / Target variability                             
PDCMETHD= 'regularMap'         / PDC algorithm used for target                  
NUMBAND =                    1 / Number of scale bands                          
FITTYPE1= 'prior   '           / Fit type used for band 1                       
PR_GOOD1=   0.9999999403953552 / Prior goodness for band 1                      
PR_WGHT1=           7.222428E7 / Prior weight for band 1                        
PDC_TOT =   0.9671307802200317 / PDC total goodness metric for target           
PDC_TOTP=    47.21249771118164 / PDC_TOT percentile compared to mod/out         
PDC_COR =   0.9740492701530457 / PDC correlation goodness metric for target     
PDC_CORP=    16.63991355895996 / PDC_COR percentile compared to mod/out         
PDC_VAR =   0.9556516408920288 / PDC variability goodness metric for target     
PDC_VARP=   25.864707946777344 / PDC_VAR percentile compared to mod/out         
PDC_NOI =   0.9995011687278748 / PDC noise goodness metric for target           
PDC_NOIP=     99.8350601196289 / PDC_NOI percentile compared to mod/out         
PDC_EPT =                  1.0 / PDC earth point goodness metric for target     
PDC_EPTP=    76.93851470947266 / PDC_EPT percentile compared to mod/out         
CHECKSUM= '5dE25ZB05bB05ZB0'   / HDU checksum updated 2015-09-07T00:10:06Z

In [4]:

from astropy.table import Table
data = Table(hdulist[1].data)
data

Out[4]:

TIME	TIMECORR	CADENCENO	SAP_FLUX	SAP_FLUX_ERR	SAP_BKG	SAP_BKG_ERR	PDCSAP_FLUX	PDCSAP_FLUX_ERR	SAP_QUALITY	PSF_CENTR1	PSF_CENTR1_ERR	PSF_CENTR2	PSF_CENTR2_ERR	MOM_CENTR1	MOM_CENTR1_ERR	MOM_CENTR2	MOM_CENTR2_ERR	POS_CORR1	POS_CORR2
float64	float32	int32	float32	float32	float32	float32	float32	float32	int32	float64	float32	float64	float32	float64	float32	float64	float32	float32	float32
169.520244707	0.00276929	2965	1.04089e+07	80.6514	17445.8	6.81535	nan	nan	256	nan	nan	nan	nan	658.98490534	5.95852e-06	48.3103679836	6.3673e-05	-0.00062082	-0.0417439
169.540678885	0.00276986	2966	1.02251e+07	79.9479	17458.6	6.81904	nan	nan	256	nan	nan	nan	nan	658.984756055	6.02055e-06	48.3645953434	6.32755e-05	-0.000859906	-0.0414842
169.561113062	0.00277044	2967	1.00993e+07	79.4513	17447.1	6.81976	nan	nan	256	nan	nan	nan	nan	658.984589033	6.06324e-06	48.3997750666	6.29868e-05	-0.00104417	-0.0411753
169.581547239	0.00277102	2968	1.01336e+07	79.5974	17440.8	6.83005	nan	nan	8576	nan	nan	nan	nan	658.984422904	6.04838e-06	48.390495693	6.31236e-05	-0.00107995	-0.0406672
169.601981415	0.0027716	2969	1.02778e+07	80.1858	17444.1	6.82776	nan	nan	393600	nan	nan	nan	nan	658.984122755	5.99424e-06	48.3358836774	6.35298e-05	-0.00127112	-0.0405857
169.622415492	0.00277217	2970	1.02574e+07	80.072	17459.8	6.82085	nan	nan	256	nan	nan	nan	nan	658.984584229	5.99411e-06	48.3514164486	6.3458e-05	-0.00144995	-0.0403933
169.642849768	0.00277275	2971	9.98968e+06	79.0959	17483.7	6.82168	nan	nan	384	nan	nan	nan	nan	658.984344926	6.09163e-06	48.4532447707	6.3e-05	-0.00170547	-0.0402648
169.663283944	0.00277332	2972	9.77098e+06	78.2334	17483.0	6.82184	nan	nan	384	nan	nan	nan	nan	658.98438983	6.17017e-06	48.4774463806	6.24124e-05	-0.00160952	-0.0400706
169.68371802	0.0027739	2973	1.00747e+07	79.3498	17492.0	6.82596	nan	nan	8576	nan	nan	nan	nan	658.984164904	6.05106e-06	48.4005135707	6.30872e-05	-0.00186426	-0.0400735
169.704152195	0.00277448	2974	1.16034e+07	85.2008	17466.9	6.84612	nan	nan	384	nan	nan	nan	nan	658.98393699	5.57952e-06	48.1491895595	6.69151e-05	-0.00199162	-0.0393474
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
258.283726632	0.00270851	7309	1.31303e+07	90.7647	17822.3	6.88246	1.31696e+07	90.8242	0	nan	nan	nan	nan	658.976566253	5.28842e-06	48.556622252	6.65606e-05	-0.00338648	0.0333565
258.304159514	0.00270789	7310	1.29471e+07	90.1064	17815.1	6.87775	1.29858e+07	90.159	0	nan	nan	nan	nan	658.976785138	5.28611e-06	48.5489829917	6.71985e-05	-0.00362269	0.0343312
258.324592595	0.00270728	7311	1.27119e+07	89.2577	17829.2	6.86093	1.27434e+07	89.3091	0	nan	nan	nan	nan	658.97690318	5.31089e-06	48.4355448483	6.73667e-05	-0.00367731	0.0329374
258.345025577	0.00270666	7312	1.24702e+07	88.3857	17808.9	6.86379	1.25071e+07	88.4538	0	nan	nan	nan	nan	658.976988051	5.34708e-06	48.2808913315	6.72824e-05	-0.00375214	0.0333151
258.365458458	0.00270604	7313	1.22394e+07	87.5396	17795.8	6.86721	1.22863e+07	87.6056	8192	nan	nan	nan	nan	658.976720505	5.3801e-06	48.1381236156	6.72635e-05	-0.00439123	0.0352742
258.385891439	0.00270542	7314	nan	nan	nan	nan	nan	nan	32800	nan	nan	nan	nan	nan	nan	nan	nan	nan	nan
258.406324519	0.0027048	7315	1.16802e+07	85.472	17810.2	6.84704	1.17144e+07	85.5362	0	nan	nan	nan	nan	658.978009467	5.4965e-06	48.0239072802	6.65793e-05	-0.00313238	0.0327484
258.4267574	0.00270418	7316	1.13527e+07	84.2527	17809.3	6.85194	1.1391e+07	84.311	0	nan	nan	nan	nan	658.977539358	5.58907e-06	48.1019007076	6.57001e-05	-0.00367803	0.033705
258.44719038	0.00270356	7317	1.10874e+07	83.2584	17793.6	6.84324	1.11284e+07	83.3245	0	nan	nan	nan	nan	658.977533777	5.6692e-06	48.1753740325	6.50243e-05	-0.00364378	0.0347512
258.46762346	0.00270294	7318	1.09189e+07	82.6193	17790.1	6.85121	1.09536e+07	82.6836	0	nan	nan	nan	nan	658.977162406	5.72111e-06	48.2262843684	6.45673e-05	-0.0039766	0.0342977

In [5]:

df = data.to_pandas()[['TIME', 'SAP_FLUX', 'SAP_FLUX_ERR']]
df.shape

Out[5]:

(4354, 3)

In [6]:

df = df.dropna()
df.shape

Out[6]:

(4083, 3)

In [7]:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy.stats import LombScargle

plt.style.use('seaborn-whitegrid')

Data and Window¶

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)

/Users/jakevdp/anaconda/envs/python3.5/lib/python3.5/site-packages/astropy/stats/lombscargle/implementations/fast_impl.py:123: RuntimeWarning: invalid value encountered in true_divide
  power = (YC * YC / CC + YS * YS / SS)

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')

No description has been provided for this image

Size of required grid¶

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)

f_ny = 48.9393547811
T = 88.947378753
n_grid = 21765.1366282