<jupyter_start><jupyter_text>Short Cadence - fitting real Kepler light curves<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# from tqdm import tqdm
# from astropy.table import Table
# import astropy.units as u
# import scipy
# from scipy import signal


from astropy.io import fits

# Using astropy BLS and scipy curve_fit to fit transit
from astropy.timeseries import BoxLeastSquares

import dynesty

# And importing `photoeccentric`
import photoeccentric as ph

# # Random stuff
import scipy.constants as c
import os

import random

%load_ext autoreload
%autoreload 2

%matplotlib inline

plt.rcParams['figure.figsize'] = [12, 8]
nwalk = 64
nsteps_d = 1500
ndiscard_d = 1000
arrlen = (nsteps_d-ndiscard_d)*nwalk
smass_kg = 1.9885e30  # Solar mass (kg)
srad_m = 696.34e6 # Solar radius (m)<jupyter_output><empty_output><jupyter_text>Reading in Sample Data<jupyter_code>spectplanets = pd.read_csv('spectplanets.csv')
muirhead_comb = pd.read_csv('muirhead_comb.csv')
muirhead_comb_lums = pd.read_csv('muirhead_comb_lums.csv')
muirhead_lamost = pd.read_csv('muirhead_lamost.csv')
muirheadKOIs = pd.read_csv('MuirheadKOIs.csv')
def get_lc_files(KIC, KICs, lcpath):

    lclist = []

    for i in range(len(KICs)):
        templst = []
        for subdir, dirs, files in os.walk(lcpath):
            for file in files:
                if str(KICs[i]) in file:
                    templst.append(os.path.join(subdir, file))
        lclist.append(templst)


    files = lclist[int(np.argwhere(KICs==KIC))]

    return files



def sigma_clip(time, flux, fluxerr, sig=4):
    for i in tqdm(range(len(flux))):
        if flux[i] > np.nanmean(flux) + sig*np.nanstd(flux) or flux[i] < np.nanmean(flux) - sig*np.nanstd(flux):
            time[i] = np.nan
            flux[i] = np.nan
            fluxerr[i] = np.nan
            
    return time, flux, fluxerr
            
lcpath = '/Users/ssagear/Dropbox (UFL)/Research/MetallicityProject/photoeccentric/notebooks/datafiles/sample_SCs/'
KICs = np.sort(np.unique(np.array(muirhead_comb['KIC'])))
KIC = 4139816
kepid = KIC
NKOI = 812
KOI = 'KOI 812'

files = get_lc_files(KIC, KICs, lcpath)
hdus, t, f, fe, starts, stops = ph.get_stitched_lcs(files, KIC)
alltime_noclip = []
allflux_noclip = []
allfluxerr_noclip = []

for sublist in t:
    for item in sublist:
        alltime_noclip.append(item)
        
for sublist in f:
    for item in sublist:
        allflux_noclip.append(item)
        
for sublist in fe:
    for item in sublist:
        allfluxerr_noclip.append(item)
# alltime, allflux, allfluxerr = sigma_clip(alltime_noclip, allflux_noclip, allfluxerr_noclip, sig=4)

time, flux, flux_err = np.array(alltime_noclip), np.array(allflux_noclip), np.array(allfluxerr_noclip)

plt.scatter(time, flux, s=1)
plt.title(KOI)

plt.ticklabel_format(useOffset=False, style='plain')<jupyter_output><empty_output><jupyter_text>KOI 812 SCs<jupyter_code>kepname = spectplanets.loc[spectplanets['kepid'] == kepid].kepler_name.values[0]
isodf = pd.read_csv("datafiles/isochrones/iso_lums_" + str(kepid) + ".csv")

mstar = isodf["mstar"].mean()
mstar_err = isodf["mstar"].std()

rstar = isodf["radius"].mean()
rstar_err = isodf["radius"].std()
period, period_uerr, period_lerr, rprs, rprs_uerr, rprs_lerr, a_arc, a_uerr_arc, a_lerr_arc, i, e_arc, w_arc = ph.planet_params_from_archive(spectplanets, kepname)

# We calculate a_rs to ensure that it's consistent with the spec/Gaia stellar density.
a_rs = ph.calc_a(period*86400.0, mstar*smass_kg, rstar*srad_m)
a_rs_err = np.mean((a_uerr_arc, a_lerr_arc))

print('Stellar mass (Msun): ', mstar, 'Stellar radius (Rsun): ', rstar)
print('Period (Days): ', period, 'Rp/Rs: ', rprs)
print('a/Rs: ', a_rs)
print('i (deg): ', i)
# TODO: Take midpoints from csv files

mpts = pd.read_csv('tmidpoints_new.csv', comment='#')
transitmpt = mpts.loc[mpts['KOI (star)'] == NKOI]['Transit Epoch (BJD)'].values[0]
midpoints = np.concatenate((np.arange(transitmpt, starts[0], -period), np.arange(transitmpt, stops[-1], period)))
midpoints
timez = []
fluxez = []
ferrz = []

for i in range(len(midpoints)):
    
    try:
        m, b, t1bjd, t1, fnorm, fe1 = ph.do_linfit(time, flux, flux_err, midpoints[i], 100, 5)
        timez.append(t1bjd)
        fluxez.append(fnorm)
        ferrz.append(fe1)

    except TypeError:
        continue
        
ttimes = np.array(timez).flatten()
tfluxes = np.array(fluxez).flatten()
tferr = np.array(ferrz).flatten()

tfluxes = np.nan_to_num(tfluxes, nan=1.0)
tferr = np.nan_to_num(tferr, nan=np.nanmedian(tferr))

def get_KIC(KOI, muirhead_comb):
    return muirhead_comb[muirhead_comb['KOI'] == str(KOI)].KIC.item()
priortransform = [2., 19., 0.1, 0.3, 20., 40., 2., 88., 0.1, transitmpt]

dres, pdist, rdist, adist, idist = ph.fit_keplc_dynesty(812, midpoints, ttimes, tfluxes, tferr, priortransform, arrlen, 150, spectplanets, muirhead_comb)
np.savetxt('812scperiods.csv', pdist, delimiter=',')
np.savetxt('812scrprs.csv', rdist, delimiter=',')
np.savetxt('812scars.csv', adist, delimiter=',')
np.savetxt('812scinc.csv', idist, delimiter=',')
np.savetxt('812sct0.csv', t0dist, delimiter=',')
rho_star, mass, radius = ph.find_density_dist_symmetric(mstar, mstar_err, rstar, rstar_err, arrlen)
T14dist = ph.get_T14(pdist, rdist, adist, idist)
T14errs = ph.get_sigmas(T14dist)

T23dist = ph.get_T23(pdist, rdist, adist, idist)
T23errs = ph.get_sigmas(T23dist)<jupyter_output><empty_output><jupyter_text>Get $g$<jupyter_code>gs, rho_c = ph.get_g_distribution(rho_star, pdist, rdist, T14dist, T23dist)

g_mean = ph.mode(gs)
g_sigma = np.mean(np.abs(ph.get_sigmas(gs)))
g_mean
g_sigma
g = g_mean
gerr = g_sigma<jupyter_output><empty_output><jupyter_text>Dynesty (e,w) Fitting LIkelihood Functions<jupyter_code>ndim = 2

def loglike(theta):
    """The log-likelihood function."""

    w, e = theta
    
    model = (1+e*np.sin(w*(np.pi/180.)))/np.sqrt(1-e**2)
    sigma2 = gerr ** 2
    
    return -0.5 * np.sum((g - model) ** 2 / sigma2 + np.log(sigma2))

def unife_prior_transform(utheta):
    """Uniform eccentricity prior"""
    
    uw, ue = utheta
    w = 360.*uw-90.
    e = 1. * ue
    
    return w, e


def betae_prior_transform(utheta):
    """Uniform eccentricity prior"""
    
    uw, ue = utheta
    w = 360.*uw-90.
    
    a, b = 0.867, 3.03
    e = scipy.stats.beta.ppf(ue, a, b)
    
    return w, e<jupyter_output><empty_output><jupyter_text>Uniform Prior<jupyter_code>dsampler = dynesty.DynamicNestedSampler(loglike, unife_prior_transform, ndim=2, bound='multi', sample='rstagger')
dsampler.run_nested()
dres = dsampler.results
from dynesty import plotting as dyplot

fig, axes = dyplot.cornerplot(dres, show_titles=True, title_kwargs={'y': 1.04}, labels=["w", "e"],
                              fig=plt.subplots(2, 2, figsize=(8, 8)))<jupyter_output><empty_output><jupyter_text>Beta-distribution Prior<jupyter_code>dsampler = dynesty.DynamicNestedSampler(loglike, betae_prior_transform, ndim=2, bound='multi', sample='rstagger')
dsampler.run_nested()
dres = dsampler.results
fig, axes = dyplot.cornerplot(dres, show_titles=True, title_kwargs={'y': 1.04}, labels=["w", "e"],
                              fig=plt.subplots(2, 2, figsize=(8, 8)))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14214601769911506, 'fraction_numerical': 0.024059734513274336, 'mean_word_length': 4.334070796460177, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '234825', 'license': 'MIT', 'n_tokens_mistral': 3116, 'n_tokens_neox': 2895, 'n_words': 569, 'path': 'notebooks/.ipynb_checkpoints/SC_keplerlc_fitting_dynesty-checkpoint.ipynb', 'repo_name': 'ssagear/metallicity', 'size': '7232'}