stevengiacalone / triceratops

Tool for Rating Interesting Candidate Exoplanets and Reliability Analysis of Transits Originating from Proximate Stars
https://triceratops.readthedocs.io
MIT License
20 stars 7 forks source link

issue with the plot of the different scenarios #29

Open astroMama opened 4 months ago

astroMama commented 4 months ago

I used the Juliet package to calculate the phase and I have already checked out the hyper parameters. I've run into a problem when I wanted to plot all the different scenarios( target.plot_fits(time=time, flux_0=flux, flux_err_0=flux_err). I don't obtain a good transit fit for the transit planet scenario even though I got 48% of probability for this scenario and a false positive probability of 0.36. Another strange thing is that there is no transit for some scenarios and weird one for others. Does someone have an idea about how to solve this issue ? Here is the code :

import numpy as np import pandas as pd import time from lightkurve import TessLightCurve import matplotlib.pyplot as plt

import triceratops.triceratops as tr

ID = 232967440 sectors = np.array([14]) target = tr.target(ID=ID, sectors=sectors)

ap = np.array([[1457, 1018], [1458, 1018], [1457, 1017], [1458, 1017]])

target.plot_field(sector=14, ap_pixels=ap)

print(target.stars)

apertures = np.array([ap]) target.calc_depths(tdepth=0.006551069280966382, all_ap_pixels=apertures)

print(target.stars)

import juliet import numpy as np import matplotlib.pyplot as plt import pandas as pd

target_ID = "232967440" save_path = r"/Users/maelravaux/Documents/QLP_lc_for_Mael/" lc_filename = save_path + '{}_lc_data.csv'.format(target_ID) lc_data = pd.read_csv(lc_filename, header=None,skiprows =1) lc_data = lc_data.T t = lc_data[0].values f = lc_data[1].values ferr = lc_data[2].values P_orb = 7.069946973323354

Put data arrays into dictionaries so we can fit it with juliet:

times, flux, flux_err= {},{},{} times['TESS'], flux['TESS'], flux_err['TESS'] = t,f,ferr

Plot data:

plt.figure() plt.errorbar(t, f, ferr) plt.xlim([np.min(t),np.max(t)])

priors = {}

Name of the parameters to be fit:

params = ['P_p1','t0_p1','r1_p1','r2_p1','q1_TESS','q2_TESS','ecc_p1','omega_p1',\ 'rho', 'mdilution_TESS', 'mflux_TESS', 'sigma_w_TESS']

Distributions:

dists = ['normal','normal','uniform','uniform','uniform','uniform','fixed','fixed',\ 'loguniform', 'fixed', 'normal', 'loguniform'] Rstar= 0.9554

Hyperparameters

hyperps = [[7.069946973323354,0.1],[1688.7192202723054,0.1], [0.,1],[0.08093867605147975,0.1],[0.1, 0.3], [0.1, 0.3], Rstar/15.749572358139467 , 90.,\ [100., 10000], 1.0, [0.,0.1], [0.1, 1000]]

for param, dist, hyperp in zip(params, dists, hyperps): priors[param] = {} priors[param]['distribution'], priors[param]['hyperparameters'] = dist, hyperp

Load and fit dataset with juliet:

dataset = juliet.load(priors=priors, t_lc = times, y_lc = flux, \ yerr_lc = flux_err, out_folder = 'TIC232967440')

results = dataset.fit()

Extract median model and the ones that cover the 68% credibility band around it:

transit_model, transit_up68, transit_low68 = results.lc.evaluate('TESS', return_err=True)

To plot the phased lighcurve we need the median period and time-of-transit center:

P, t0 = np.median(results.posteriors['posterior_samples']['P_p1']),\ np.median(results.posteriors['posterior_samples']['t0_p1'])

Get phases:

phases = juliet.get_phases(dataset.times_lc['TESS'], P, t0)# avec max P et max To

flux= dataset.data_lc["TESS"] flux_err=dataset.errors_lc["TESS"] time=phases*P_orb

lc_binsize = (time.max()-time.min())/100 lc = TessLightCurve(time=time, flux=flux, flux_err=flux_err).bin(time_bin_size=lc_binsize)

target.calc_probs(time=lc.time.value, flux_0=lc.flux.value, flux_err_0=np.mean(lc.flux_err.value), P_orb=P_orb)

df_results = target.probs print("FPP =", np.round(target.FPP, 14)) print("NFPP =", np.round(target.NFPP, 14)) print(df_results)

target.plot_fits(time=time, flux_0=flux, flux_err_0=flux_err)

FPPs = np.zeros(20) for i in range(20): target.calc_probs(time=lc.time.value, flux_0=lc.flux.value, flux_err_0=np.mean(lc.flux_err.value), P_orb=P_orb, parallel=True, verbose=0) FPPs[i] = target.FPP

meanFPP = np.round(np.mean(FPPs), 14) stdvFPP = np.round(np.std(FPPs), 14) print("FPP =", meanFPP, "+/-", stdvFPP)

stevengiacalone commented 3 months ago

It's difficult to find the issue without seeing your light curve. Can you please send me a zipped folder of your notebook and the related data? Please email it to steven_giacalone@berkeley.edu