;+ ; NAME: ; eebls ; ; PURPOSE: ; >>>>>>>>>>>> This routine computes BLS spectrum <<<<<<<<<<<<<< ; ; [ see Kovacs, Zucker & Mazeh 2002, A&A, Vol. 391, 369 ] ; ; This is the slightly modified version of the original BLS routine ; by considering Edge Effect (EE) as suggested by ; Peter R. McCullough [ pmcc@stsci.edu ]. ; ; This modification was motivated by considering the cases when ; the low state (the transit event) happened to be devided between ; the first and last bins. In these rare cases the original BLS ; yields lower detection efficiency because of the lower number of ; data points in the bin(s) covering the low state. ; ; For further comments/tests see www.konkoly.hu/staff/kovacs.html ; ; N.B. This is an IDL version of the code originally created by ; Kovacs et al. I (Brian Jackson - decaelus@gmail.com) take no credit at all ; for developing this routine. ; ; CATEGORY: ; Astrophysical ; ; CALLING SEQUENCE: ; ; eebls,n,t,x,u,v,nf,fmin,df,nb,qmi,qma,$ ; p,bper,bpow,depth,qtran,in1,in2 ; ; INPUTS: ; n = number of data points ; t = array {t(i)}, containing the time values of the time series ; x = array {x(i)}, containing the data values of the time series ; u = temporal/work/dummy array, must be dimensioned in the ; calling program in the same way as {t(i)} ; v = the same as {u(i)} ; nf = number of frequency points in which the spectrum is computed ; fmin = minimum frequency (MUST be > 0) ; df = frequency step ; nb = number of bins in the folded time series at any test period ; qmi = minimum fractional transit length to be tested ; qma = maximum fractional transit length to be tested ; ; OUTPUTS: ; p = array {p(i)}, containing the values of the BLS spectrum ; at the i-th frequency value -- the frequency values are ; computed as f = fmin + (i-1)*df ; bper = period at the highest peak in the frequency spectrum ; bpow = value of {p(i)} at the highest peak ; depth= depth of the transit at *bper* ; qtran= fractional transit length [ T_transit/bper ] ; in1 = bin index at the start of the transit [ 0 < in1 < nb+1 ] ; in2 = bin index at the end of the transit [ 0 < in2 < nb+1 ] ; ; RESTRICTIONS: ; -- *fmin* MUST be greater than *1/total time span* ; -- *nb* MUST be lower than *nbmax* ; -- Dimensions of arrays {y(i)} and {ibi(i)} MUST be greater than ; or equal to *nbmax*. ; -- The lowest number of points allowed in a single bin is equal ; to MAX(minbin,qmi*N), where *qmi* is the minimum transit ; length/trial period, *N* is the total number of data points, ; *minbin* is the preset minimum number of the data points per ; bin. ; ; EXAMPLE: ; ;This is about 1 month of Kepler data ; n = 1440. ; ;Orbital period of ~5 hrs ; per = 20. + randomu(seed) - 0.5 ; t = dindgen(n) ; x = replicate(1.0, n) ; u = dblarr(n) ; v = dblarr(n) ; ;According to Kovacs+ (2002), this is the relationship between nf, meant to ; ; represent the nber of independent frequencies to test, and n ; nf = fix(n^(0.83)) ; fmin = 2./(max(t)-min(t)) ; fmax = 2./(t[1]-t[0]) ; df = (1./(0.05*per)-1/(5.*per))/(nf*10) ; ;must be less than 2000 ; nb = 50 ; folded_time = t mod per ; mid = randomu(seed)*per ; ;fake transit ; width = 0.03*per ; depth = 100d-6 ; ind = where(abs(folded_time-mid) le width) ; fake_trans = replicate(1.0, n_elements(folded_time)) ; fake_trans[ind] = 1.0-depth ; x = fake_trans ; qma = 0.05 ; qmi = 0.01 ; eebls,n,t,x,u,v,nf,fmin,df,nb,qmi,qma,$ ; p,bper,bpow,depth,qtran,in1,in2 ; freq = dindgen(nf)*df + fmin ; plot, 1./freq, p, /xlog, psym=2 ; print, bper, per ; ; MODIFICATION HISTORY: ; Original fortran code written by G. Kovacs, S. Zucker, & T. Mazeh ; Ported to IDL by Brian Jackson (decaelus@gmail.com), 2012 May ;- pro eebls,n,t,x,u,v,nf,fmin,df,nb,qmi,qma,$ p,bper,bpow,depth,qtran,in1,in2 y = dblarr(2000) ibi = dblarr(2000) p = dblarr(nf) minbin = 5 nbmax = 2000 ;Check that input parameters are in the right ranges if(nb gt nbmax) then begin print, "NB > NBMAX" stop end;if tot = t[n-1]-t[0] if(fmin lt 1./tot) then begin print, "fmin < 1/T!!" stop end;if if(qmi ge qma) then begin print, "qmi >= qma!" stop end;if rn = double(n) kmi = fix(qmi*double(nb)) if(kmi lt 1) then kmi = 1 kma = fix(qma*double(nb))+1 kkmi = fix(rn*qmi) if(kkmi lt minbin) then kkmi = minbin bpow = 0. ; ; The following variables are defined for the extension ; of arrays ibi() and y() [ see below ] ; nb1 = nb+1 nbkma = nb+kma ; ;================================= ; Set temporal time series ;================================= ; s=0.0d0 t1=t[0] u = t - t1 s = total(x) ;sum divided by total number of points s=s/rn ;v is the data minus the sum divided by the number of points v = x - s ; ;****************************** ; Start period search * ;****************************** ; for jf=0,nf-1 do begin ;trial frequency f0=fmin+df*double(jf) ;corresponding trial period p0=1.0d0/f0 ; ;====================================================== ; Compute folded time series with *p0* period ;====================================================== ; ;zero out arrays y = dblarr(2000) ibi = dblarr(2000) j = nb ; ;iterate through all the data points for i=0,n-1 do begin ;ph is orbital phase ph = u[i]*f0 ;remove integer part of ph to fold data ph = ph-fix(ph) j = 1 + fix(nb*ph) ibi[j] = ibi[j] + 1. y[j] = y[j] + v[i] end;for i ; ;----------------------------------------------- ; Extend the arrays ibi() and y() beyond ; nb by wrapping ; for j=nb1-1,nbkma-1 do begin jnb = j-nb ibi[j] = ibi[jnb] y[j] = y[jnb] end;for j ;----------------------------------------------- ; ;=============================================== ; Compute BLS statistics for this period ;=============================================== ; power=0.0d0 ; for i=0,nb-1 do begin s = 0.0d0 k = 0 kk = 0 nb2 = i+kma for j=i,nb2-1 do begin k = k+1 kk = kk+ibi[j] s = s+y[j] if(k lt kmi) then continue if(kk lt kkmi) then continue rn1 = double(kk) pow = s*s/(rn1*(rn-rn1)) if(pow lt power) then continue power = pow jn1 = i jn2 = j rn3 = rn1 s3 = s end;for j end;for i ; power = sqrt(power) p[jf] = power ; if(power lt bpow) then continue bpow = power in1 = jn1 in2 = jn2 qtran = rn3/rn depth = -s3*rn/(rn3*(rn-rn3)) bper = p0 ; end;for jf ; ; Edge correction of transit end index ; if(in2 gt nb) then in2 = in2-nb ; return end;pro