`photoeccentric` Tutorial

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

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

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
import pickle
import scipy
import random

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

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

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

# Using dynesty to do the same with nested sampling
import dynesty

# And importing `photoeccentric`
import photoeccentric as ph

%load_ext autoreload
%autoreload 2

# pandas display option
pd.set_option('display.float_format', lambda x: '%.5f' % x)

spectplanets = pd.read_csv('../datafiles/spectplanets.csv')
muirhead_comb = pd.read_csv('../datafiles/muirhead_comb.csv')
muirheadKOIs = pd.read_csv('../datafiles/MuirheadKOIs.csv')
lcpath = '../datafiles/sample_lcs'

plt.rcParams['figure.figsize'] = [20, 10]

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

I’ll define the conversions between solar mass -> kg and solar radius -> meters for convenience.

smass_kg = 1.9885e30  # Solar mass (kg)
srad_m = 696.34e6 # Solar radius (m)

The Sample

I’m using the sample of “cool KOIs” from Muirhead et al. 2013, and their properites from spectroscopy published here.

I’m reading in several .csv files containing data for this sample. The data includes spectroscopy data from Muirhead et al. (2013), stellar and planet parameters from the Kepler archive, and distances/luminosities from Gaia.

# ALL Kepler planets from exo archive
planets = pd.read_csv('../datafiles/cumulative_kois.csv')

# Take the Kepler planet archive entries for the planets in Muirhead et al. 2013 sample
spectplanets = pd.read_csv('../datafiles/spectplanets.csv')

# Kepler-Gaia Data
kpgaia = Table.read('../datafiles/kepler_dr2_4arcsec.fits', format='fits').to_pandas();

# Kepler-Gaia data for only the objects in our sample
muirhead_gaia = pd.read_csv("../datafiles/muirhead_gaia.csv")

# Combined spectroscopy data + Gaia/Kepler data for our sample
muirhead_comb = pd.read_csv('../datafiles/muirhead_comb.csv')

# Only targets from table above with published luminosities from Gaia
muirhead_comb_lums = pd.read_csv('../datafiles/muirhead_comb_lums.csv')

Defining a “test planet”

I’m going to pick a planet from our sample to test how well photoeccentric works. Here, I’m picking KOI 818.01 (Kepler-691 b), a super-Earth orbiting an M dwarf. Exoplanet Catalog Entry

It has an orbital period of about 8 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 and using the masses/radii from the matching isochrones to constrian the stellar density.

nkoi = 818.01

I’ll read in a file with the MESA stellar isochrones for low-mass stars. I’ll use ph.fit_isochrone_lum() to find the subset of stellar isochrones that are consistent with a certain stellar parameters form Kepler-691 (Teff, Mstar, Rstar, and Gaia luminosity).

# # Read in MESA isochrones
isochrones = pd.read_csv('../datafiles/isochrones_sdss_spitzer_lowmass.dat', sep='\s\s+', engine='python')

Using ph.fit_isochrone_lum() to match isochrones to stellar data:

koi818 = muirhead_comb.loc[muirhead_comb['KOI'] == '818']

iso_lums = ph.fit_isochrone_lum(koi818, isochrones)

100%|██████████| 738479/738479 [00:06<00:00, 112685.12it/s]

# Write to csv, then read back in (prevents notebook from lagging)
iso_lums.to_csv("iso_lums_" + str(nkoi) + ".csv")

isodf = pd.read_csv("iso_lums_" + str(nkoi) + ".csv")

Define a KeplerStar object, and use ph.get_stellar_params and the fit isochrones to get the stellar parameters.

SKOI = int(np.floor(float(nkoi)))
print('KOI', SKOI)

star = ph.KeplerStar(SKOI)
star.get_stellar_params(isodf)

KOI 818

print('Stellar Mass (Msol): ', star.mstar)
print('Stellar Radius (Rsol): ', star.rstar)
print('Average Stellar Density (kg m^-3): ', star.rho_star)

Stellar Mass (Msol):  0.5691952380952375
Stellar Radius (Rsol):  0.5488321428571431
Average Stellar Density (kg m^-3):  4849.937558132769

Define a KOI object.

koi = ph.KOI(nkoi, SKOI, isodf)
koi.get_KIC(muirhead_comb)

print('KIC', koi.KIC)

KIC 4913852

Creating a fake light curve based on a real planet

I’m pulling the planet parameters of Kepler-691 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.

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.

Let’s force the inclination to be 90 degrees for this example.

koi.planet_params_from_archive(spectplanets)
koi.calc_a(koi.mstar, koi.rstar)

print('Stellar mass (Msun): ', koi.mstar, 'Stellar radius (Rsun): ', koi.rstar)
print('Period (Days): ', koi.period, 'Rp/Rs: ', koi.rprs)
print('a/Rs: ', koi.a_rs)

koi.i = 90.
print('i (deg): ', koi.i)

Stellar mass (Msun):  0.5691952380952375 Stellar radius (Rsun):  0.5488321428571431
Period (Days):  8.11437482 Rp/Rs:  0.037204
a/Rs:  25.63672846518325
i (deg):  90.0

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’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.

Test Case 1: \(e = 0.0\), \(\omega = 90.0\)

I define a cadence length (~30 minutes, in days) that matches the Kepler long-cadence integration time, so I can create a fake light curve that integrates over the same time as real Kepler light curves.

I want to replicate the real Kepler light curve as closely as possible. So I am taking these parameters fromm the light curves

# Define the working directory
direct = 'tutorial01/' + str(nkoi) + '/e_0.0_w_90.0/'

First, reading in the light curves that I have saved for this planet.

KICs = np.sort(np.unique(np.array(muirhead_comb['KIC'])))
KOIs = np.sort(np.unique(np.array(muirhead_comb['KOI'])))

files = ph.get_lc_files(koi.KIC, KICs, lcpath)

# Stitching the light curves together, preserving the time stamps
koi.get_stitched_lcs(files)

# Getting the midpoint times
koi.get_midpoints()

starttime = koi.time[0]
endtime = koi.time[-1]

# 30 minute cadence
cadence = 0.02142857142857143

time = np.arange(starttime, endtime, cadence)

# Define e and w, calculate flux from transit model
e = 0.0
w = 90.0

params = batman.TransitParams()       #object to store transit parameters
params.t0 = koi.epoch                 #time of inferior conjunction
params.per = koi.period               #orbital period
params.rp = koi.rprs                  #planet radius (in units of stellar radii)
params.a = koi.a_rs                   #semi-major axis (in units of stellar radii)
params.inc = koi.i                    #orbital inclination (in degrees)
params.ecc = e                        #eccentricity
params.w = w                          #longitude of periastron (in degrees)
params.limb_dark = "nonlinear"        #limb darkening model
params.u = [0.5, 0.1, 0.1, -0.1]      #limb darkening coefficients [u1, u2, u3, u4]

t = time
m = batman.TransitModel(params, t, supersample_factor = 29, exp_time = 0.0201389)

flux = m.light_curve(params)

time = time
flux = flux

plt.plot(time-2454900, flux)
plt.xlim(109,109.7)
plt.xlabel('Time (BJD-2454900)')
plt.ylabel('Normalized Flux')
plt.title('Transit LC Model')

Text(0.5, 1.0, 'Transit LC Model')

To create a light curve with a target signal to noise ratio, we need the transit duration, number of transits, and the number of points in each transit, and the transit depth.

tduration = koi.dur/24.0

N = round(ph.get_N_intransit(tduration, cadence))
ntransits = len(koi.midpoints)
depth = koi.rprs**2

The magnitude of each individual error bar:

errbar = ph.get_sigma_individual(60, N, ntransits, depth)

Adding gaussian noise to produce a light curve with the target SNR:

(NB: the noise is gaussian and uncorrelated, unlike the noise in real Kepler light curves)

noise = np.random.normal(0,errbar,len(time))
nflux = flux+noise

flux_err = np.array([errbar]*len(nflux))

plt.errorbar(time-2454900, nflux, yerr=flux_err, fmt='o')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.xlim(109,109.7)
plt.title('Transit LC Model with Noise')

Text(0.5, 1.0, 'Transit LC Model with Noise')

Fitting the transit

photoeccentric includes functionality to fit using juliet with multinest.

First, I’ll fit the transit shape with juliet. \(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 is not 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 close to the Kepler archive parameters.

koi.time = time
koi.flux = nflux
koi.flux_err = flux_err

Let’s just do the first 7 transits.

koi.midpoints = koi.midpoints[0:7]

koi.remove_oot_data(7, 6)

plt.errorbar(koi.time_intransit, koi.flux_intransit, yerr=koi.fluxerr_intransit, fmt='o')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.title('Transit LC Model with Noise')

Text(0.5, 1.0, 'Transit LC Model with Noise')

nlive=1000
nsupersample=29
exptimesupersample=0.0201389

dataset, results = koi.do_tfit_juliet(direct, nsupersample=nsupersample, exptimesupersample=exptimesupersample, nlive=nlive)

analysing data from tutorial01/818.01/e_0.0_w_90.0/jomnest_.txt

res = pd.read_table(direct + 'posteriors.dat')

# Print transit fit results from Juliet
res

	# Parameter Name	Median	Upper 68 CI	Lower 68 CI
0	P_p1	8.11419	0.00069	0.00072
1	t0_p1	2455009.34043	0.00091	0.00091
2	p_p1	0.03726	0.00090	0.00087
3	b_p1	0.34351	0.21283	0.21383
4	q1_KEPLER	0.41629	0.33855	0.25182
5	q2_KEPLER	0.45065	0.34201	0.28467
6	a_p1	23.96618	2.29076	2.99624
7	mflux_KEPLER	-0.00018	0.00007	0.00007
8	sigma_w_KEPLER	3.47350	33.94203	3.13059

# Save fit planet parameters to variables for convenience
per_f = res.iloc[0][2]
t0_f = res.iloc[1][2]
rprs_f = res.iloc[2][2]
b_f = res.iloc[3][2]
a_f = res.iloc[6][2]

i_f = np.arccos(b_f*(1./a_f))*(180./np.pi)

Below, I plot the transit fit corner plot.

Because I input \(e = 0.0\), the transit fitter should return close to the same parameters I input (because the transit fitter always requires \(e = 0.0\)).

# Plot the transit fit corner plot
p = results.posteriors['posterior_samples']['P_p1']
t0 = results.posteriors['posterior_samples']['t0_p1']
rprs = results.posteriors['posterior_samples']['p_p1']
b = results.posteriors['posterior_samples']['b_p1']
a = results.posteriors['posterior_samples']['a_p1']
inc = np.arccos(b*(1./a))*(180./np.pi)

params = ['Period', 't0', 'rprs', 'inc', 'a']

fs = np.vstack((p, t0, rprs, inc, a))
fs = fs.T

figure = corner.corner(fs, labels=params)

# Plot the data:
plt.errorbar(dataset.times_lc['KEPLER']-2454900, dataset.data_lc['KEPLER'], \
             yerr = dataset.errors_lc['KEPLER'], fmt = '.', alpha = 0.1)

# Plot the model:
plt.plot(dataset.times_lc['KEPLER']-2454900, results.lc.evaluate('KEPLER'), c='r')

plt.xlabel('Time (BJD)-2454900')
plt.xlim(np.min(dataset.times_lc['KEPLER'])-2454900, np.max(dataset.times_lc['KEPLER'])-2454900)
plt.ylim(0.998, 1.002)
plt.ylabel('Relative flux')
plt.title('KOI 818.01 Kepler LC with Transit Fit')
plt.show()

Determining T14 and T23

A crucial step to determining the \((e, w)\) distribution from the transit is calculating the total and full transit durations. T14 is the total transit duration (the time between first and fourth contact). T23 is the full transit duration (i.e. the time during which the entire planet disk is in front of the star, the time between second and third contact.)

Here, I’m using equations 14 and 15 from this textbook. We calculate T14 and T23 assuming the orbit must be circular, and using the fit parameters assuming the orbit is circular. (If the orbit is not circular, T14 and T23 will not be correct – but this is what we want, because they will differ from the true T14 and T23 in a way that reveals the eccentricity of the orbit.)

koi.calc_durations()

print('Total Transit Duration: ', np.mean(koi.T14_dist), '-/+', koi.T14_errs, 'hours')
print('Full Transit Duration: ', np.mean(koi.T23_dist), '-/+', koi.T23_errs, 'hours')

Total Transit Duration:  0.10485137003852651 -/+ (0.006525117642141129, 0.006687207067543222) hours
Full Transit Duration:  0.09561957555025782 -/+ (0.005883030475635945, 0.006018819870806616) hours

Get \(g\)

Finally, we can use all the values above to determine \(\rho_{circ}\). \(\rho_{circ}\) is what we would calculate the stellar density to be if we knew that the orbit was definitely perfectly circular. We will compare \(\rho_{circ}\) to \(\rho_{star}\) (the true, observed stellar density we calculated from spectroscopy/Gaia), and get \(g(e, w)\):

\[\rho_{\star}(e, \omega) = g(e, \omega)^{-3} \rho_{circ}\]

which is also defined as

\[g(e, \omega) = \frac{1 + e sin(\omega)}{\sqrt{1-e^2}}\]

Thus, if the orbit is circular \((e = 0)\), then \(g\) should equal 1. If the orbit is not circular \((e != 0)\), then \(\rho_{circ}\) should differ from \(\rho_{star}\), and \(g\) should be something other than 1. We can draw a \((e, w)\) distribution based on the value we calcaulte for \(g(e,w)\)!

ph.get_g_distribution() will help us determine the value of g. This function takes the observed \(\rho_{star}\) as well as the fit (circular) transit parameters and calculated transit durations, and calculates \(\rho_{circ}\) and \(g(e,w)\) based on equations 6 and 7 in Dawson & Johnson 2012.

Print \(g\) and \(\sigma_{g}\):

koi.get_gs()

g_mean = koi.g_mean
g_sigma = koi.g_sigma

g_mean

0.9898376567084125

g_sigma

0.14832653472072282

The mean of \(g\) is about 1.0, which means that \(\rho_{circ}\) agrees with \(\rho_{star}\) and the eccentricity of this transit must be zero, which is exactly what we input! We can take \(g\) and \(\sigma_{g}\) and use MCMC (emcee) to determine the surface of most likely \((e,w)\).

photoeccentric has the probability function for \((e,w)\) from \(g\) built in to ph.log_probability().

koi.do_eccfit(direct)

18743it [01:53, 164.43it/s, batch: 15 | bound: 0 | nc: 1 | ncall: 88101 | eff(%): 21.274 | loglstar:   -inf <  1.908 <  1.837 | logz:  0.870 +/-  0.052 | stop:  0.954]

with open(direct + '/kepewdres.pickle', 'rb') as f:
    ewdres = pickle.load(f)

labels = ["w", "e"]

fig = corner.corner(ewdres.samples, labels=labels, title_kwargs={"fontsize": 12}, truths=[w, e], plot_contours=True)

And here is the corner plot for the most likely values of \((e, w)\) that correspond to \(g = 1\). The \(e\) distribution peaks at 0!

Test Case 2: \(e = 0.3\), \(\omega = 90.0\)

Now let’s repeat this example with an eccentricity of 0.3 at periapse.

# Define the working directory
direct = 'tutorial01/' + str(nkoi) + '/e_0.3_w_90.0/'

# 30 minute cadence
cadence = 0.02142857142857143

time = np.arange(starttime, endtime, cadence)

# Define e and w, calculate flux from transit model
e = 0.3
w = 90.0

params = batman.TransitParams()       #object to store transit parameters
params.t0 = koi.epoch                 #time of inferior conjunction
params.per = koi.period               #orbital period
params.rp = koi.rprs                  #planet radius (in units of stellar radii)
params.a = koi.a_rs                   #semi-major axis (in units of stellar radii)
params.inc = koi.i                    #orbital inclination (in degrees)
params.ecc = e                        #eccentricity
params.w = w                          #longitude of periastron (in degrees)
params.limb_dark = "nonlinear"        #limb darkening model
params.u = [0.5, 0.1, 0.1, -0.1]      #limb darkening coefficients [u1, u2, u3, u4]

t = time
m = batman.TransitModel(params, t, supersample_factor = 29, exp_time = 0.0201389)

flux = m.light_curve(params)

time = time
flux = flux

plt.plot(time-2454900, flux)
plt.xlim(109,109.7)
plt.xlabel('Time (BJD-2454900)')
plt.ylabel('Normalized Flux')
plt.title('Transit LC Model')

Text(0.5, 1.0, 'Transit LC Model')

To create a light curve with a target signal to noise ratio, we need the transit duration, number of transits, and the number of points in each transit, and the transit depth.

tduration = koi.dur/24.0

N = round(ph.get_N_intransit(tduration, cadence))
ntransits = len(koi.midpoints)
depth = koi.rprs**2

The magnitude of each individual error bar:

errbar = ph.get_sigma_individual(60, N, ntransits, depth)

Adding gaussian noise to produce a light curve with the target SNR:

noise = np.random.normal(0,errbar,len(time))
nflux = flux+noise

flux_err = np.array([errbar]*len(nflux))

plt.errorbar(time-2454900, nflux, yerr=flux_err, fmt='o')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.xlim(109,109.7)
plt.title('Transit LC Model with Noise')

Text(0.5, 1.0, 'Transit LC Model with Noise')

Fitting the transit

photoeccentric includes functionality to fit using juliet with multinest.

First, I’ll fit the transit shape with juliet. \(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 is not 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 close to the Kepler archive parameters.

koi.time = time
koi.flux = nflux
koi.flux_err = flux_err

Let’s just do the first 7 transits.

koi.midpoints = koi.midpoints[0:7]

koi.remove_oot_data(7, 6)

plt.errorbar(koi.time_intransit, koi.flux_intransit, yerr=koi.fluxerr_intransit, fmt='o')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.title('Transit LC Model with Noise')

Text(0.5, 1.0, 'Transit LC Model with Noise')

nlive=1000
nsupersample=29
exptimesupersample=0.0201389

dataset, results = koi.do_tfit_juliet(direct, nsupersample=nsupersample, exptimesupersample=exptimesupersample, nlive=nlive)

analysing data from tutorial01/818.01/e_0.3_w_90.0/jomnest_.txt

res = pd.read_table(direct + 'posteriors.dat')

# Print transit fit results from Juliet
res

	# Parameter Name	Median	Upper 68 CI	Lower 68 CI
0	P_p1	8.11447	0.00020	0.00020
1	t0_p1	2455009.34045	0.00050	0.00053
2	p_p1	0.03779	0.00074	0.00081
3	b_p1	0.25005	0.20836	0.16863
4	q1_KEPLER	0.24137	0.33555	0.16504
5	q2_KEPLER	0.20592	0.33467	0.14925
6	a_p1	34.01496	1.53172	2.81683
7	mflux_KEPLER	-0.00005	0.00001	0.00001
8	sigma_w_KEPLER	1.76799	10.52057	1.49476

# Save fit planet parameters to variables for convenience
per_f = res.iloc[0][2]
t0_f = res.iloc[1][2]
rprs_f = res.iloc[2][2]
b_f = res.iloc[3][2]
a_f = res.iloc[6][2]

i_f = np.arccos(b_f*(1./a_f))*(180./np.pi)

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\)).

# Plot the transit fit corner plot
p = results.posteriors['posterior_samples']['P_p1']
t0 = results.posteriors['posterior_samples']['t0_p1']
rprs = results.posteriors['posterior_samples']['p_p1']
b = results.posteriors['posterior_samples']['b_p1']
a = results.posteriors['posterior_samples']['a_p1']
inc = np.arccos(b*(1./a))*(180./np.pi)

params = ['Period', 't0', 'rprs', 'inc', 'a']

fs = np.vstack((p, t0, rprs, inc, a))
fs = fs.T

figure = corner.corner(fs, labels=params)

# Plot the data:
plt.errorbar(dataset.times_lc['KEPLER']-2454900, dataset.data_lc['KEPLER'], \
             yerr = dataset.errors_lc['KEPLER'], fmt = '.', alpha = 0.1)

# Plot the model:
plt.plot(dataset.times_lc['KEPLER']-2454900, results.lc.evaluate('KEPLER'), c='r')

# Plot portion of the lightcurve, axes, etc.:
plt.xlabel('Time (BJD)-2454900')
plt.xlim(np.min(dataset.times_lc['KEPLER'])-2454900, np.max(dataset.times_lc['KEPLER'])-2454900)
plt.ylim(0.998, 1.002)
plt.ylabel('Relative flux')
plt.title('KOI 818.01 Kepler LC with Transit Fit')
plt.show()

Determining T14 and T23

koi.calc_durations()

print('Total Transit Duration: ', np.mean(koi.T14_dist), '-/+', koi.T14_errs, 'hours')
print('Full Transit Duration: ', np.mean(koi.T23_dist), '-/+', koi.T23_errs, 'hours')

Total Transit Duration:  0.07624977166072124 -/+ (0.001938522023045397, 0.0022380238235224087) hours
Full Transit Duration:  0.0700019311019035 -/+ (0.0018442657065528278, 0.00209380377961145) hours

Get \(g\)

Print \(g\) and \(\sigma_{g}\):

koi.get_gs()

g_mean = koi.g_mean
g_sigma = koi.g_sigma

g_mean

1.3096409093931336

g_sigma

0.09639732695175818

The mean of \(g\) is about 1.3, which means that \(\rho_{circ}\) agrees with \(\rho_{star}\) and the eccentricity of this transit must be non-zero. We can take \(g\) and \(\sigma_{g}\) and use nested sampling (dynesty) to determine the surface of most likely \((e,w)\).

photoeccentric has the probability function for \((e,w)\) from \(g\) built in to ph.log_probability().

koi.do_eccfit(direct)

18966it [01:42, 184.48it/s, batch: 13 | bound: 0 | nc: 1 | ncall: 167069 | eff(%): 11.352 | loglstar:   -inf <  2.339 <  2.218 | logz:  0.239 +/-  0.065 | stop:  0.978]

with open(direct + '/kepewdres.pickle', 'rb') as f:
    ewdres = pickle.load(f)

labels = ["w", "e"]

fig = corner.corner(ewdres.samples, labels=labels, title_kwargs={"fontsize": 12}, truths=[w, e], plot_contours=True)

And here is the corner plot for the most likely values of \((e, w)\) that correspond to \(g = 1.3\). The \(e\) distribution peaks at \(e = 0.3\)!

Test Case 3: \(e = 0.3\), \(\omega = 270.0\)

Now let’s repeat this example with an eccentricity of 0.3 at apoapse.

# Define the working directory
direct = 'tutorial01/' + str(nkoi) + '/e_0.3_w_270.0/'

# 30 minute cadence
cadence = 0.02142857142857143

time = np.arange(starttime, endtime, cadence)

# Define e and w, calculate flux from transit model
e = 0.3
w = 270.0

params = batman.TransitParams()       #object to store transit parameters
params.t0 = koi.epoch                        #time of inferior conjunction
params.per = koi.period                       #orbital period
params.rp = koi.rprs                       #planet radius (in units of stellar radii)
params.a = koi.a_rs                        #semi-major axis (in units of stellar radii)
params.inc = koi.i                      #orbital inclination (in degrees)
params.ecc = e                       #eccentricity
params.w = w                        #longitude of periastron (in degrees)
params.limb_dark = "nonlinear"        #limb darkening model
params.u = [0.5, 0.1, 0.1, -0.1]      #limb darkening coefficients [u1, u2, u3, u4]

t = time
m = batman.TransitModel(params, t, supersample_factor = 29, exp_time = 0.0201389)

flux = m.light_curve(params)

time = time
flux = flux

plt.plot(time-2454900, flux)
plt.xlim(109,109.7)
plt.xlabel('Time (BJD-2454900)')
plt.ylabel('Normalized Flux')
plt.title('Transit LC Model')

Text(0.5, 1.0, 'Transit LC Model')

To create a light curve with a target signal to noise ratio, we need the transit duration, number of transits, and the number of points in each transit, and the transit depth.

tduration = koi.dur/24.0

N = round(ph.get_N_intransit(tduration, cadence))
ntransits = len(koi.midpoints)
depth = koi.rprs**2

The magnitude of each individual error bar:

errbar = ph.get_sigma_individual(60, N, ntransits, depth)

Adding gaussian noise to produce a light curve with the target SNR:

noise = np.random.normal(0,errbar,len(time))
nflux = flux+noise

flux_err = np.array([errbar]*len(nflux))

plt.errorbar(time-2454900, nflux, yerr=flux_err, fmt='o')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.xlim(109,109.7)
plt.title('Transit LC Model with Noise')

Text(0.5, 1.0, 'Transit LC Model with Noise')

Fitting the transit

koi.time = time
koi.flux = nflux
koi.flux_err = flux_err

Let’s just do the first 7 transits.

koi.midpoints = koi.midpoints[0:7]

koi.remove_oot_data(7, 6)

plt.errorbar(koi.time_intransit, koi.flux_intransit, yerr=koi.fluxerr_intransit, fmt='o')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.title('Transit LC Model with Noise')

Text(0.5, 1.0, 'Transit LC Model with Noise')

nlive=1000
nsupersample=29
exptimesupersample=0.0201389

dataset, results = koi.do_tfit_juliet(direct, nsupersample=nsupersample, exptimesupersample=exptimesupersample, nlive=nlive)

analysing data from tutorial01/818.01/e_0.3_w_270.0/jomnest_.txt

res = pd.read_table(direct + 'posteriors.dat')

# Print transit fit results from Juliet
res

	# Parameter Name	Median	Upper 68 CI	Lower 68 CI
0	P_p1	8.11407	0.00024	0.00023
1	t0_p1	2455009.34021	0.00057	0.00057
2	p_p1	0.03688	0.00069	0.00065
3	b_p1	0.36570	0.17872	0.22039
4	q1_KEPLER	0.35974	0.26895	0.17441
5	q2_KEPLER	0.29702	0.31601	0.18794
6	a_p1	17.29780	1.10237	1.78141
7	mflux_KEPLER	-0.00037	0.00002	0.00002
8	sigma_w_KEPLER	3.34227	24.49567	3.00753

# Save fit planet parameters to variables for convenience
per_f = res.iloc[0][2]
t0_f = res.iloc[1][2]
rprs_f = res.iloc[2][2]
b_f = res.iloc[3][2]
a_f = res.iloc[6][2]

i_f = np.arccos(b_f*(1./a_f))*(180./np.pi)

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\)).

# Plot the transit fit corner plot
p = results.posteriors['posterior_samples']['P_p1']
t0 = results.posteriors['posterior_samples']['t0_p1']
rprs = results.posteriors['posterior_samples']['p_p1']
b = results.posteriors['posterior_samples']['b_p1']
a = results.posteriors['posterior_samples']['a_p1']
inc = np.arccos(b*(1./a))*(180./np.pi)

params = ['Period', 't0', 'rprs', 'inc', 'a']

fs = np.vstack((p, t0, rprs, inc, a))
fs = fs.T

figure = corner.corner(fs, labels=params)

# Plot the data:
plt.errorbar(dataset.times_lc['KEPLER']-2454900, dataset.data_lc['KEPLER'], \
             yerr = dataset.errors_lc['KEPLER'], fmt = '.', alpha = 0.1)

# Plot the model:
plt.plot(dataset.times_lc['KEPLER']-2454900, results.lc.evaluate('KEPLER'), c='r')

# Plot portion of the lightcurve, axes, etc.:
plt.xlabel('Time (BJD)-2454900')
plt.xlim(np.min(dataset.times_lc['KEPLER'])-2454900, np.max(dataset.times_lc['KEPLER'])-2454900)
plt.ylim(0.998, 1.002)
plt.ylabel('Relative flux')
plt.title('KOI 818.01 Kepler LC with Transit Fit')
plt.show()

Determining T14 and T23

koi.calc_durations()

print('Total Transit Duration: ', np.mean(koi.T14_dist), '-/+', koi.T14_errs, 'hours')
print('Full Transit Duration: ', np.mean(koi.T23_dist), '-/+', koi.T23_errs, 'hours')

Total Transit Duration:  0.1453917067458673 -/+ (0.0024099101024557534, 0.0029148854669333313) hours
Full Transit Duration:  0.1328974496125958 -/+ (0.002286337605429317, 0.002373614271698049) hours

Get \(g\)

Print \(g\) and \(\sigma_{g}\):

koi.get_gs()

g_mean = koi.g_mean
g_sigma = koi.g_sigma

g_mean

0.7138312721212587

g_sigma

0.056832865957797074

The mean of \(g\) is about 1.0, which means that \(\rho_{circ}\) agrees with \(\rho_{star}\) and the eccentricity of this transit must be zero, which is exactly what we input! We can take \(g\) and \(\sigma_{g}\) and use nested sampling (dynesty) to determine the surface of most likely \((e,w)\).

photoeccentric has the probability function for \((e,w)\) from \(g\) built in to ph.log_probability().

koi.do_eccfit(direct)

22184it [02:18, 160.69it/s, batch: 15 | bound: 0 | nc: 1 | ncall: 262862 | eff(%):  8.439 | loglstar:   -inf <  2.868 <  2.745 | logz:  0.281 +/-  0.074 | stop:  0.795]

with open(direct + '/kepewdres.pickle', 'rb') as f:
    ewdres = pickle.load(f)

labels = ["w", "e"]

fig = corner.corner(ewdres.samples, labels=labels, title_kwargs={"fontsize": 12}, truths=[w, e], plot_contours=True)

And here is the corner plot for the most likely values of \((e, w)\) that correspond to \(g = 0.7\). This \(e\) distribution peaks at 0.3 too!

photoeccentric Tutorial

The Sample

Defining a “test planet”

Creating a fake light curve based on a real planet

Test Case 1: \(e = 0.0\), \(\omega = 90.0\)

Fitting the transit

Determining T14 and T23

Get \(g\)

Test Case 2: \(e = 0.3\), \(\omega = 90.0\)

Fitting the transit

Determining T14 and T23

Get \(g\)

Test Case 3: \(e = 0.3\), \(\omega = 270.0\)

Fitting the transit

Determining T14 and T23

Get \(g\)

`photoeccentric` Tutorial