Bug

[svaverbe]

While testing the grazing template in TLS, an issue appears with a fit that has a clear underestimation of the transit duration. It is not clear if this is due to a bug in the code.

import numpy as np import transitleastsquares from transitleastsquares import transitleastsquares,cleaned_array,catalog_info,transit_mask from matplotlib.backends.backend_pdf import PdfPages import matplotlib.pyplot as plt from astroquery.mast import Catalogs

def TLS(tic,Tmag,ra,dec,time,flux,detrended_flux,flux_err,Pmin,Pmax,duration_grid_step,ntransits_min):

Applies the TLS algorithm to the detrended data to search for transits

ab, mass, mass_min, mass_max, radius, radius_min, radius_max = catalog_info(TIC_ID=np.int(tic))

if np.isnan(mass):
   mass_param=1.0
else:
   mass_param=mass

if np.isnan(radius):
   radius_param=1.0
else:
   radius_param=radius

model = transitleastsquares(time, detrended_flux,flux_err)

results = model.power(u=ab,period_min=Pmin,period_max=Pmax,duration_grid_step=duration_grid_step,M_star=mass_param, \

R_star=radius_param,M_star_max=1.5,transit_template='grazing')

results = model.power(u=ab,period_min=Pmin,period_max=Pmax,duration_grid_step=duration_grid_step,transit_template='grazing')
ntransits=results.transit_count

transit_depth=1.0-results.depth

if ~np.isnan(ntransits): 

   with PdfPages('detection_TLS_'+str(tic)+'.pdf') as pdf:

        filename='detection_output_TLS_'+str(tic)+'.dat'
        file1=open(filename,"w")

        plt.figure(figsize=(8, 6))
        plt.scatter(time, flux, color='blue',s=2, alpha=0.5, zorder=0,label='Normalized PDCSAP flux')
        plt.scatter(time, detrended_flux, color='red',s=2, alpha=0.5, zorder=0,label='Detrended PDCSAP flux')
        plt.legend(loc='best')
        plt.xlabel('BKJD')
        plt.ylabel('Normalized flux')
        plt.title('Detrended PDCSAP lightcurve:TIC '+str(tic))
        pdf.savefig()
        plt.close()

        plt.figure(figsize=(8, 6))
        ax = plt.gca()
        ax.axvline(results.period, alpha=0.4, lw=3)
        plt.xlim(np.min(results.periods), np.max(results.periods))

        for n in range(2, 10):
            ax.axvline(n*results.period, alpha=0.4, lw=1, linestyle="dashed")
            ax.axvline(results.period / n, alpha=0.4, lw=1, linestyle="dashed")

        plt.ylabel(r'SDE')
        plt.xlabel('Period (days)')
        plt.semilogx(results.periods, results.power, color='black', lw=0.5)
        plt.title('P='+str("{:.4f}".format(results.period))+' d')
        pdf.savefig()
        plt.close()

        plt.figure(figsize=(8, 6))
        plt.plot((results.model_folded_phase-0.5)*results.period, results.model_folded_model, color='red',label='transit model')
        plt.scatter((results.folded_phase-0.5)*results.period, results.folded_y, color='blue', s=10, alpha=0.5, zorder=2,label='folded data') 
        plt.xlim(-5*results.duration,5*results.duration)
        plt.legend(loc='best')
        plt.xlabel('Time from central transit(days)')
        plt.ylabel('Normalized flux')
        pdf.savefig()
        plt.close()

        plt.figure(figsize=(8, 6))
        in_transit = transit_mask(time, results.period, results.duration, results.T0)
        plt.scatter(time[in_transit], detrended_flux[in_transit], color='red', s=2, zorder=0,label='in transit data')
        plt.scatter(time[~in_transit], detrended_flux[~in_transit], color='blue', alpha=0.5, s=2, zorder=0,label='out of transit data')
        plt.plot(results.model_lightcurve_time, results.model_lightcurve_model, alpha=0.5, color='red', zorder=1)
        plt.xlim(time.min(), time.max())
        plt.legend(loc='best')
        plt.xlabel('BKJD')
        plt.ylabel('Normalized flux')
        pdf.savefig()
        plt.close()

        print('TIC id:',tic,file=file1)
        print('Period', format(results.period, '.5f'), 'd at T0=', results.T0,file=file1)
        print(len(results.transit_times), 'transit times in time series:', ['{0:0.5f}'.format(i) for i in results.transit_times],file=file1)
        print('Number of data points during each unique transit', results.per_transit_count,file=file1)            
        print('Transit duration (days)', format(results.duration, '.5f'),file=file1)
        print('Transit depths (mean)', results.transit_depths,file=file1)

        print('The number of transits with intransit data points', results.distinct_transit_count,file=file1)
        print('The number of transits with no intransit data points', results.empty_transit_count,file=file1)
        print('Transit depth', format(results.depth, '.5f'), '(at the transit bottom)',file=file1)
        print('Transit depth uncertainties', results.transit_depths_uncertainties,file=file1)
        print('Overall transit SNR',results.snr,file=file1)
        print('Transit SNRs', results.snr_per_transit,file=file1)
        print('Odd even mismatch',results.odd_even_mismatch,file=file1)

        plt.figure(figsize=(8, 6))
        plt.errorbar(results.transit_times,results.transit_depths,yerr=results.transit_depths_uncertainties,
        fmt='o', color='red')
        plt.plot((time.min(), time.max()),(np.mean(results.transit_depths), np.mean(results.transit_depths)),
        color='black', linestyle='dashed')
        plt.plot((time.min(), time.max()), (1, 1), color='black')
        plt.xlabel('BKJD')
        plt.ylabel('Normalized flux')
        pdf.savefig()
        plt.close()

        file1.close()

return results

tic=236887394 Pmin=0.5 Pmax=100.0 ntransits_min=2 duration_grid_step=1.1 # default=1.1

ticinfo = Catalogs.query_criteria(catalog="Tic", ID=tic) # stellar data from TIC

Vmag=ticinfo['Vmag'][0] Tmag=ticinfo['Tmag'][0] ra=ticinfo[0]['ra'] dec=ticinfo[0]['dec']

fname='detrendedlightcurve'+str(tic)+'.dat'

data=np.genfromtxt(fname,dtype='float',delimiter=" ") # detrended flux time series time=data[:,0] flux=data[:,1] detrended_flux=data[:,2] nobs=len(time)

flux_err=np.ones((nobs))

results=TLS(tic,Tmag,ra,dec,time,flux,detrended_flux,flux_err,Pmin,Pmax,duration_grid_step,ntransits_min)

print(results)

detrended_lightcurve_236887394.txt detection_TLS_236887394.pdf