<jupyter_start><jupyter_text>** Import **<jupyter_code>import kplr
import numpy as np
import matplotlib.pyplot as plt
from exoparams import PlanetParams
import batman
from lmfit import minimize, Parameters, Model
from tqdm import tqdm_notebook
plt.ion()
%matplotlib inline<jupyter_output><empty_output><jupyter_text>** Kepler Data Load **<jupyter_code>client = kplr.API()
koi3 = client.koi(3.01)
lc3 = koi3.get_light_curves(short_cadence=True)<jupyter_output><empty_output><jupyter_text>** Planet Parameters **<jupyter_code>keplaunch = 2454833.0
# Data below from Huber et al. 2017
kep3period = 4.88782433
kep3t0 = 2454957.812464 - keplaunch
kep3aoR = 14.64
kep3RpRs = 0.05856
kep3inc = 88.99
kep3ecc = 0.26493
kep3omeg = -162.149
kep3u1 = 0.646
kep3u2 = 0.048
while kep3t0 > keplaunch:
    kep3t0 -= kep3period<jupyter_output><empty_output><jupyter_text>** Arrays for Graphing Short Cadence Without NaNs **<jupyter_code>times, fluxs, ferrs = [], [], []
for lc in lc3:
    with lc.open() as f:
        hdu_data = f[1].data
        
        baddata   = hdu_data['sap_quality'] != 0
        baddata  &= np.isnan(hdu_data["time"])
        baddata  &= np.isnan(hdu_data["pdcsap_flux"])
        baddata  &= np.isnan(hdu_data["pdcsap_flux_err"])
        
        baddata  &= ~np.isfinite(hdu_data["time"])
        baddata  &= ~np.isfinite(hdu_data["pdcsap_flux"])
        baddata  &= ~np.isfinite(hdu_data["pdcsap_flux_err"])
        baddata   = np.where(baddata)
        
        if f[0].header['OBSMODE'] == 'short cadence':
            times.append(np.delete(hdu_data["time"]           , baddata, None))
            fluxs.append(np.delete(hdu_data["pdcsap_flux"]    , baddata, None))
            ferrs.append(np.delete(hdu_data["pdcsap_flux_err"], baddata, None))
plt.figure(figsize=(9, 9))
for i in range(len(times)):
    plt.errorbar(times[i], fluxs[i], yerr = ferrs[i], fmt = '.')

plt.xlabel("Time (KJD)", fontsize = 20)
plt.ylabel("Flux (Photons)", fontsize = 20)<jupyter_output><empty_output><jupyter_text>** Spliting and Phasing Transits **<jupyter_code>plt.figure(figsize=(9,9))
for pl in range(len(times)):
    phasei = ((times[pl] - kep3t0) % kep3period)/kep3period
    phasei[phasei > .5] -= 1
    plt.errorbar(phasei, fluxs[pl], yerr = ferrs[pl], fmt = '.')
epochsNums = []
for pl in tqdm_notebook(range(len(times)), total = len(times)):
    for k in range(int(max(times[pl]) / kep3period)):
        if (kep3t0 + k*kep3period) >= min(times[pl]) and (kep3t0 + k*kep3period) <= max(times[pl]):
            epochsNums.append(k)<jupyter_output><empty_output><jupyter_text>**Slice Epochs Into Arrays**<jupyter_code>sliceWidth = 0.5
timeSlices , fluxSlices, ferrSlices, epochs = [], [], [], []
for pl in tqdm_notebook(range(len(times)), total = len(times)):
    for k in epochsNums:
        tranCenter    = kep3t0 + k * kep3period
        if tranCenter >= min(times[pl]) and tranCenter <= max(times[pl]):
            transitSlices = abs(times[pl] - tranCenter) < sliceWidth
            if len(transitSlices) > 0 and len(times[pl][transitSlices]) > 0:
                timeSlices.append(times[pl][transitSlices])
                fluxSlices.append(fluxs[pl][transitSlices])
                ferrSlices.append(ferrs[pl][transitSlices])
                epochs.append(tranCenter)
plt.figure(figsize=(9,9))
for eachEpoch in range(len(timeSlices)):
    plt.errorbar(timeSlices[eachEpoch] ,fluxSlices[eachEpoch], ferrSlices[eachEpoch], fmt = '.' )
plt.figure(figsize = (9,9))
for timeNow, fluxNow, ferrNow, kep3epochNow in zip(timeSlices, fluxSlices, ferrSlices, epochs):
    plt.errorbar(timeNow - kep3epochNow, fluxNow / np.median(fluxNow) \
                 , ferrNow / np.median(fluxNow), fmt='.')<jupyter_output><empty_output><jupyter_text>**Outliers**<jupyter_code>plt.figure(figsize = (9,9))
nsigma = 5
for itran in range(len(epochs)):
    fluxdiff = np.diff(fluxSlices[itran])
    plt.plot(fluxdiff)
    plt.axhline(np.nanmedian(fluxdiff) + nsigma * np.nanstd(fluxdiff), \
                color = 'k', linestyle = '-', linewidth = 1)
    plt.axhline(np.nanmedian(fluxdiff) - nsigma * np.nanstd(fluxdiff), \
                color = 'k', linestyle = '-', linewidth = 1)
fluxSlicesTemp = fluxSlices.copy()

nsigma = 5
for itran in range(len(epochs)):
    fluxdiff = np.diff(fluxSlicesTemp[itran])
    outliers = np.where(abs(fluxdiff - np.nanmedian(fluxdiff)) > nsigma * np.nanstd(fluxdiff))[0]
    for o in outliers:
        beforeOutlier = list(fluxSlicesTemp[itran][o - 10: o + 0])
        afterOutlier  = list(fluxSlicesTemp[itran][o + 1: o + 11])
        fluxSlicesTemp[itran][o] = np.median(beforeOutlier + afterOutlier) #+ fluxdiff[o]
    fluxdiff = np.diff(fluxSlicesTemp[itran])
    print(itran, np.where(abs(fluxdiff - np.nanmedian(fluxdiff)) > nsigma * np.nanstd(fluxdiff))[0]) 
plt.figure(figsize = (9,9))
nsigma = 5
for itran in range(len(epochs)):
    fluxdiff = np.diff(fluxSlicesTemp[itran])
    plt.plot(fluxdiff)
    plt.axhline(np.nanmedian(fluxdiff) + nsigma * np.nanstd(fluxdiff), \
                color = 'k', linestyle = '-', linewidth = 1)
    plt.axhline(np.nanmedian(fluxdiff) - nsigma * np.nanstd(fluxdiff), \
                color = 'k', linestyle = '-', linewidth = 1)
for itran in range(len(epochs)):
    ferrNans = np.where(np.isnan(ferrSlices[itran]))
    for n in ferrNans:
        ferrNan = np.nanmedian(ferrSlices[itran])
        ferrSlices[itran][n] = ferrNan
    
    print(itran, np.where(np.isnan(ferrSlices[itran])))<jupyter_output>0 (array([], dtype=int64),)
1 (array([], dtype=int64),)
2 (array([], dtype=int64),)
3 (array([], dtype=int64),)
4 (array([], dtype=int64),)
5 (array([], dtype=int64),)
6 (array([], dtype=int64),)
7 (array([], dtype=int64),)
8 (array([], dtype=int64),)
9 (array([], dtype=int64),)
10 (array([], dtype=int64),)
11 (array([], dtype=int64),)
12 (array([], dtype=int64),)
13 (array([], dtype=int64),)
14 (array([], dtype=int64),)
15 (array([], dtype=int64),)
16 (array([], dtype=int64),)
17 (array([], dtype=int64),)
18 (array([], dtype=int64),)
19 (array([], dtype=int64),)
20 (array([], dtype=int64),)
21 (array([], dtype=int64),)
22 (array([], dtype=int64),)
23 (array([], dtype=int64),)
24 (array([], dtype=int64),)
25 (array([], dtype=int64),)
26 (array([], dtype=int64),)
27 (array([], dtype=int64),)
28 (array([], dtype=int64),)
29 (array([], dtype=int64),)
30 (array([], dtype=int64),)
31 (array([], dtype=int64),)
32 (array([], dtype=int64),)
33 (array([], dtype=int64),)
34 (array([], dtype=int6[...]<jupyter_text>** Fitting the transit with BATMAN **<jupyter_code>timeSlicesMod = []
for itran in range(len(epochs)):
    timeSlicesMod += [np.linspace(np.nanmin(timeSlices[itran]), np.nanmax(timeSlices[itran]), \
                                  len(timeSlices[itran]))]
plt.figure(figsize=(4,4))
for itran in range(len(epochs)):
    plt.plot(timeSlices[itran])
    plt.plot(timeSlicesMod[itran], '-')
tparams = batman.TransitParams()
tparams.t0 = kep3t0                     #time of inferior conjunction
tparams.per = kep3period                #orbital period 
tparams.rp = kep3RpRs                   #planet radius (in units of stellar radii)
tparams.a = kep3aoR                     #semi-major axis (in units of stellar radii)
tparams.inc = kep3inc                   #orbital inclination (in degrees)
tparams.ecc = kep3ecc                   #eccentricity
tparams.w = kep3omeg                    #longitude of periastron (in degrees)
tparams.limb_dark = "quadratic"         #limb darkening model
tparams.u = [kep3u1, kep3u2]            #limb darkening coefficients [u1, u2, u3, u4]

sliceLcModel = []
for itran in range(len(epochs)):
    sliceLcModel += [batman.TransitModel(tparams, timeSlicesMod[itran], transittype ='primary')]
lcKepModel = []
for itran in range(len(epochs)):
    lcKepModel += [sliceLcModel[itran].light_curve(tparams)]
plt.figure(figsize = (9,9))

for timeSlice, epoch, BatModel in zip(timeSlicesMod, epochs, lcKepModel):
    plt.plot(timeSlice - epoch, BatModel)
    plt.xlabel("Normalized Time", fontsize = 20)
    plt.ylabel("Relative Flux (Photons)", fontsize = 20)<jupyter_output><empty_output><jupyter_text>** LMFIT Model fitting **<jupyter_code>def batman_wrapper_lmfit(period, tcenter, inc, aprs, rprs, ecc, omega, u1, u2, 
                         offset, slope, curvature,
                         times, ldtype='quadratic', transitType='primary'):

    # period, tcenter, inc, aprs, rprs, ecc, omega, u1, u2 = init_params
    
    if offset == 1.0 and slope == 0.0 and curvature == 0.0:
        OoT_curvature = 1.0 # OoT == Out of transit
    else:
        OoT_curvature = offset + slope*(times-times.mean()) + curvature*(times-times.mean())**2
    
    bm_params           = batman.TransitParams() # object to store transit parameters

    bm_params.per       = period   # orbital period
    bm_params.t0        = tcenter  # time of inferior conjunction
    bm_params.inc       = inc      # inclunaition in degrees
    bm_params.a         = aprs     # semi-major axis (in units of stellar radii)
    bm_params.rp        = rprs     # planet radius (in units of stellar radii)
    bm_params.ecc       = ecc      # eccentricity
    bm_params.w         = omega    # longitude of periastron (in degrees)
    bm_params.limb_dark = ldtype   # limb darkening model # NEED TO FIX THIS
    bm_params.u         = [u1, u2] # limb darkening coefficients # NEED TO FIX THIS

    m_eclipse = batman.TransitModel(bm_params, times, transittype=transitType)# initializes model

    return m_eclipse.light_curve(bm_params) * OoT_curvature
initialParams = Parameters()

initialParams.add_many(
    ('period' , kep3period, False),
    ('tcenter' , kep3t0, True, kep3t0 - 0.1, kep3t0 + 0.1),
    ('inc' , kep3inc, True, 80., 90.),
    ('aprs' , kep3aoR, True, 10, 20),
    ('rprs' , kep3RpRs, True, 0.01, 0.1),
    ('ecc' , kep3ecc, False),
    ('omega' , kep3omeg, False),
    ('u1' , kep3u1, True, 0.6, 0.7),
    ('u2' , kep3u2, True, 0.0, 0.1),
    ('offset', 1.0, True, 0.0),
    ('slope', 0.0, True),
    ('curvature', 0.0, True))
initialParams
lc = Model(batman_wrapper_lmfit, 
              independent_vars = ['times', 'ldtype', 'transitType'])
fitResult = []
for itran in tqdm_notebook(range(len(epochs)), total = len(epochs), desc = 'fitResult'):
    fitResult += [lc.fit(fluxSlices[itran],
                   weights     = 1 / ferrSlices[itran] ** 2.,
                   params      = initialParams,
                   method      = 'powell',
                   times       = timeSlicesMod[itran],
                   ldtype      = 'quadratic',
                   transitType = 'primary')]
plt.figure(figsize=(9,9))
for timeSlice, epoch, lmModel, fluxSlice in zip(timeSlicesMod, epochs, fitResult, fluxSlices):
    plt.plot(timeSlice - epoch, fluxSlice / fluxSlice.mean(),'.')
    plt.plot(timeSlice - epoch, lmModel.best_fit / lmModel.best_fit.mean())
plt.figure(figsize=(9,9))
for timeSlice, epoch, lmModel in zip(timeSlicesMod, epochs, fitResult):
    plt.plot(timeSlice - epoch,lmModel.best_fit / lmModel.best_fit.mean())
n = 3
plt.figure(figsize=(9,9))
plt.plot(timeSlices[n], fitResult[n].best_fit )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15607040229885058, 'fraction_numerical': 0.03394396551724138, 'mean_word_length': 3.5034371209057826, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '347697', 'license': 'MIT', 'n_tokens_mistral': 4488, 'n_tokens_neox': 4177, 'n_words': 840, 'path': 'Kepler Data Analysis For HAT-P-11 b all Kepler (with Graphs).ipynb', 'repo_name': 'exowanderer/Fitting-Exoplanet-Transits', 'size': '11136'}