<jupyter_start><jupyter_text>`photoeccentric` Tutorial (Only Transit Fit) Full tutorial at https://github.com/ssagear/photoeccentric/blob/master/notebooks/photoeccentric_tutorial.ipynb In this tutorial, I will create a simulated transit based on a Kepler planet and demonstrate how to use `photoeccentric` to recover the planet's eccentricity using the photoeccentric effect [(Dawson & Johnson 2012)](https://arxiv.org/pdf/1203.5537.pdf). The code I'm using to implement the photoeccentric effect is compiled into a package called `photoeccentric`, and can be viewed/downloaded here: https://github.com/ssagear/photoeccentric<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 os

# Using `batman` to create & fit fake transit
import batman

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

import scipy.optimize as opt

# Using emcee & corner to find and plot (e, w) distribution
import emcee
import corner

# And importing `photoeccentric`
import photoeccentric as ph

%load_ext autoreload
%autoreload 2
nwalk = 64
nsteps = 3000
ndiscard = 1000
arrlen = (nsteps-ndiscard)*nwalk<jupyter_output><empty_output><jupyter_text>I'll define the conversions between solar mass -> kg and solar radius -> meters for convenience.<jupyter_code>smass_kg = 1.9885e30  # Solar mass (kg)
srad_m = 696.34e6 # Solar radius (m)<jupyter_output><empty_output><jupyter_text>The Sample I'm using the sample of "cool KOIs" from [Muirhead et al. 2013](https://iopscience.iop.org/article/10.1088/0067-0049/213/1/5), and their properites from spectroscopy published here.<jupyter_code>muirhead_data = pd.read_csv("datafiles/Muirhead2013_isochrones/muirhead_data_no_missing_data.txt", sep=" ")<jupyter_output><empty_output><jupyter_text>I'm reading in a file containing data for all Kepler planets from the Exoplanet Archive (`planets`), then only taking these data for planets in the Muirhead et al. 2013 sample (`spectplanets`).<jupyter_code># ALL Kepler planets from exo archive
planets = pd.read_csv('datafiles/exoplanetarchive/cumulative_kois.csv')

# Take the Kepler planet archive entries for the planets in Muirhead et al. 2013 sample
spectplanets = planets[planets['kepid'].isin(list(muirhead_data['KIC']))]
spectplanets = spectplanets.reset_index()<jupyter_output><empty_output><jupyter_text>Now, I'm reading in the entire Kepler/Gaia dataset from [gaia-kepler.fun]. I'm again matching these data with the objects in our sample (`muirhead_gaia`). I'm using the DR2 data with a 4 arcsecond search radius.Then I'll combine the spectroscopy data with Kepler/Gaia data for our sample.<jupyter_code># Kepler-Gaia Data

# This file is too big to upload to github -- download from https://gaia-kepler.fun (4" search radius)
kpgaia = Table.read('datafiles/Kepler-Gaia/kepler_dr2_4arcsec.fits', format='fits').to_pandas();

# Kepler-Gaia data for only the objects in our sample
muirhead_gaia = kpgaia[kpgaia['kepid'].isin(list(muirhead_data.KIC))]
muirhead_gaia = muirhead_gaia.reset_index()
muirhead_gaia.rename(columns={"index": "KIC"})

# Combined spectroscopy data + Gaia/Kepler data for our sample
muirhead_comb = pd.concat([muirhead_data, muirhead_gaia], axis=1)
muirhead_comb_nn = muirhead_comb[muirhead_comb.KOI.notnull()]

# Only targets from table above with published luminosities from Gaia
muirhead_comb_lums = muirhead_comb_nn[muirhead_comb_nn.lum_val.notnull()]<jupyter_output><empty_output><jupyter_text>Defining a "test planet" I'm going to pick a random planet from our sample to test how well `photoeccentric` works. Here, I'm picking Kepler-1582 b, a super-Earth orbiting an M dwarf [Exoplanet Catalog Entry](https://exoplanets.nasa.gov/exoplanet-catalog/2457/kepler-1582-b/). It has an orbital period of about 5 days.First, I'll use the spectroscopy data from Muirhead et al. 2013 and Gaia luminosities to constrain the mass and radius of the host star beyond the constraint published in the Exoplanet Archive. I'll do this by matching these data with stellar isochrones [MESA](https://iopscience.iop.org/article/10.3847/0004-637X/823/2/102) (check this ciation) and using the masses/radii from the matching isochrones to constrian the stellar density.<jupyter_code># Kepler ID for Kepler-1582 b
kepid = 5868793
kepname = spectplanets.loc[spectplanets['kepid'] == kepid].kepler_name.values[0]

kp1582b = muirhead_comb.loc[muirhead_comb['KIC'] == kepid]

# Read in MESA isochrones
# File is not in Github repo
isochrones = pd.read_csv('datafiles/Muirhead2013_isochrones/isochrones_sdss_spitzer_lowmass.dat', sep='\s\s+', engine='python')<jupyter_output><empty_output><jupyter_text>Using `ph.fit_isochrone_lum()` to match isochrones to stellar data<jupyter_code>iso_lums = ph.fit_isochrone_lum(kp1582b, muirhead_comb, isochrones, gaia_lum=True)
# Write to csv, then read back in (prevents python notebook from lagging)
iso_lums.to_csv("datafiles/isochrones/iso_lums_" + str(kepid) + ".csv")
isodf = pd.read_csv("datafiles/isochrones/iso_lums_" + str(kepid) + ".csv")<jupyter_output><empty_output><jupyter_text>I'm determining the mass and radius constraints of this star based on the isochrones that were consistent with the observational data.<jupyter_code>mstar = isodf["mstar"].mean()
mstar_err = isodf["mstar"].std()

rstar = isodf["radius"].mean()
rstar_err = isodf["radius"].std()<jupyter_output><empty_output><jupyter_text>Using `ph.find_density_dist_symmetric()` to create a stellar density distribution from symmetric (Gaussian) distributions based on mstar and rstar (from isochrones). Note: this does not necessarily mean the resulting density distribution will appear symmetric.<jupyter_code>rho_star, mass, radius = ph.find_density_dist_symmetric(mstar, mstar_err, rstar, rstar_err, arrlen)
plt.hist(rho_star, bins=20)
plt.xlabel('Stellar Density Histogram (kg m^-3)', fontsize=20)<jupyter_output><empty_output><jupyter_text>Creating a fake light curve based on a real planet I'm pulling the planet parameters of Kepler-1582 b from the exoplanet archive using `ph.planet_params_from_archive()`. This will give me the published period, Rp/Rs, and inclination constraints of this planet. (It will also return some other parameters, but we don't need those right now).  I'm calculating a/Rs using `ph.calc_a()`, instead of using the a/Rs constraint from the Exoplanet Archive. The reason is because a/Rs must be consistent with the density calculated above from spectroscopy/Gaia for the photoeccentric effect to work correctly, and the published a/Rs is often inconsistent. a/Rs depends on the orbital period, Mstar, and Rstar.<jupyter_code>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)<jupyter_output>Stellar mass (Msun):  0.18620280162489727 Stellar radius (Rsun):  0.21038847941706276
Period (Days):  4.83809469 Rp/Rs:  0.036066
a/Rs:  32.64384979582086
i (deg):  89.98<jupyter_text>Now, I'll create a fake transit using `batman`.I'm creating a model with the period, Rp/Rs, a/Rs, and inclination specified by the Kepler catalog entry and the density constraints I determinedI'll create the transit model with an $e$ and $w$ of my choice. This will allow me to test whether `photoeccentric` accurately recovers the $(e,w)$ combination I have input. I'll start with $e = 0.0$ and $w = 90.0$ degrees.  $e = 0.0$, $\omega = 90.0$<jupyter_code>time = np.linspace(-24, 24, 10000)
def planetlc(time, per, rp, a, inc, w):
    params = batman.TransitParams()       #object to store transit parameters
    params.t0 = 0.                        #time of inferior conjunction
    params.per = per                      #orbital period
    params.rp = rp                        #planet radius (in units of stellar radii)
    params.a = a                          #semi-major axis (in units of stellar radii)
    params.inc = inc                      #orbital inclination (in degrees)
    params.ecc = 0.0
    params.w = w                          #longitude of periastron (in degrees)
    #params.limb_dark = "linear"
    #params.u = [0.3]
    params.limb_dark = "quadratic"
    params.u = [0.1, 0.3]
    #params.limb_dark = "uniform"
    #params.u = []

    #times to calculate light curve
    m = batman.TransitModel(params, time)
    
    flux = m.light_curve(params)
    
    return flux
# Calculate flux from transit model
w = 90.0
e = 0.0
flux = planetlc(time, period, rprs, a_rs, i, w)

# Adding some gaussian noise
noise = np.random.normal(0,0.00005,len(time))
nflux = flux+noise
flux_err = np.array([0.00005]*len(nflux))
plt.errorbar(time, nflux, yerr=flux_err)
plt.xlabel('Time')
plt.ylabel('Flux')
#plt.xlim(-1, 1)<jupyter_output><empty_output><jupyter_text>Fitting the transit I'm using the Astropy BLS method to determine the period of the fake light curve, and use this as the period "guess" (`per_guess`).<jupyter_code>per_guess = ph.bls(time, nflux)<jupyter_output><empty_output><jupyter_text>Now, I'm fitting the transit shape with scipy.curve_fit. $p$, $Rp/Rs$, $a/Rs$, $i$, and $w$ are allowed to vary as free parameters. The transit fitter, `ph.planetlc_fitter`, fixes $e = 0.0$, even if the input eccentricity is not zero! This means that if e != 0, the transit fitter will fit the wrong values for $a/Rs$ and $i$ -- but they will be wrong in such a way that reveals the eccentricity of the orbit. More on that in the next section.I enter an initial guess based on what I estimate the fit parameters will be. For this one, I'll enter values pretty close to what I input.<jupyter_code>print('Period fit: ', per_guess)
# Inital guess: period, rprs, a/Rs, i, w
p0 = [per_guess, rprs, 35, 89]

dr = 'e_' + str(0) + '_w_' + str(w)
direct = '/Users/sheilasagear/Dropbox (UFL)/Research/MetallicityProject/photoeccentric/notebooks/plots_tutorial/' + dr + '/'

if not os.path.exists(direct):
    os.mkdir(direct)

# EMCEE Transit Model Fitting
res, errs, pdist, rdist, adist, idist = ph.mcmc_fitter(p0, time, nflux, flux_err, nwalk, nsteps, ndiscard, e, w, direct)

p_f, rprs_f, a_f, i_f = res
perr_f, rprserr_f, aerr_f, ierr_f = errs

# Create a light curve with the fit parameters
fit = ph.planetlc_fitter(time, p_f, rprs_f, a_f, i_f, 0.0)<jupyter_output>100%|██████████| 3000/3000 [06:20<00:00,  7.88it/s]<jupyter_text>Below, I print the original parameters and fit parameters, and overlay the fit light curve on the input light curve. Because I input $e = 0.0$, the transit fitter should return the exact same parameters I input (because the transit fitter always requires $e = 0.0$).<jupyter_code>plt.plot(time, nflux, c='blue', alpha=0.5, label='Original LC')
plt.plot(time, fit, c='red', alpha=1.0, label='Fit LC')
plt.xlim(-0.1, 0.1)
plt.legend()

print('Stellar mass (Msun): ', mstar, 'Stellar radius (Rsun): ', rstar)

print('\n')

print('Input params:')
print('Period (days): ', period, 'Rp/Rs: ', rprs)
print('a/Rs: ', a_rs)
print('i (deg): ', i)
print('w (deg): ', w)

print('\n')

print('Fit params:')
print('Period (days): ', p_f, 'Rp/Rs: ', rprs_f)
print('a/Rs: ', a_f)
print('i (deg): ', i_f)<jupyter_output>Stellar mass (Msun):  0.18620280162489727 Stellar radius (Rsun):  0.21038847941706276


Input params:
Period (days):  4.83809469 Rp/Rs:  0.036066
a/Rs:  32.64384979582086
i (deg):  89.98
w (deg):  90.0


Fit params:
Period (days):  4.838099106532112 Rp/Rs:  0.03611172346928493
a/Rs:  32.513671039509035
i (deg):  89.84574555803674
{'fraction_non_alphanumeric': 0.11059984941018991, 'fraction_numerical': 0.03865138458964277, 'mean_word_length': 4.982982982982983, 'pattern_counts': {'":': 1, '<': 55, '<?xml version=': 0, '>': 56, 'https://': 7, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '287929', 'license': 'MIT', 'n_tokens_mistral': 4591, 'n_tokens_neox': 4104, 'n_words': 1278, 'path': 'notebooks/.ipynb_checkpoints/transit_fit_only-checkpoint.ipynb', 'repo_name': 'ssagear/metallicity', 'size': '11953'}