Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

issue with the plot of the different scenarios #29

Open
astroMama opened this issue Jun 14, 2024 · 1 comment
Open

issue with the plot of the different scenarios #29

astroMama opened this issue Jun 14, 2024 · 1 comment

Comments

@astroMama
Copy link

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
Copy link
Owner

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 [email protected]

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants