PRO scargle,t,c,om,px,fmin=fmin,fmax=fmax,numf=numf, $ nu=nu,period=period,omega=omega, $ fap=fap,signi=signi,simsigni=simsigni, $ pmin=pmin,pmax=pmax,old=old, $ psdpeaksort=psdpeaksort,multiple=multiple,noise=noise, $ debug=debug,slow=slow ;+ ; NAME: ; scargle ; ; ; PURPOSE: ; Compute the lomb-scargle periodogram of an unevenly sampled ; lightcurve ; ; ; CATEGORY: ; time series analysis ; ; ; CALLING SEQUENCE: ; scargle,t,c,om,px,fmin=fmin,fmax=fmax,numf=numf,pmin=pmin,pmax=pmax, ; nu=nu,period=period,omega=omega,fap=fap,signi=signi ; ; ; INPUTS: ; time: The times at which the time series was measured ; rate: the corresponding count rates ; ; ; OPTIONAL INPUTS: ; fmin,fmax: minimum and maximum frequency (NOT ANGULAR FREQ!) ; to be used (has precedence over pmin,pmax) ; pmin,pmax: minimum and maximum PERIOD to be used ; omega: angular frequencies for which the PSD values are ; desired ; fap : false alarm probability desired ; (see Scargle et al., p. 840, and signi ; keyword). Default equal to 0.01 (99% significance) ; noise: for the normalization of the periodogram and the ; compute of the white noise simulations. If not set, equal to ; the variance of the original lc. ; multiple: number of white noise simulations for the FAP ; power level. Default equal to 0 (i.e., no simulations). ; numf: number of independent frequencies ; ; ; KEYWORD PARAMETERS: ; old : if set computing the periodogram according to J.D.Scargle ; 1982, ApJ 263, 835. If not set, computing the periodogram ; with the fast algorithm of W.H. Press and G.B. Rybicki, ; 1989, ApJ 338, 277. ; debug: print out debugging information if set ; slow: if set, a much slower but less memory intensive way to ; perform the white noise simulations is used. ; ; OUTPUTS: ; om : angular frequency of PSD ; psd : the psd-values corresponding to omega ; ; ; OPTIONAL OUTPUTS: ; nu : normal frequency (nu=omega/(2*!DPI)) ; period: period corresponding to each omega ; signi : power threshold corresponding to the given ; false alarm probabilities fap and according to the ; desired number of independent frequencies ; simsigni : power threshold corresponding to the given ; false alarm probabilities fap according to white ; noise simulations ; psdpeaksort : array with the maximum peak pro each simulation ; ; PROCEDURE: ; The Lomb Scargle PSD is computed according to the ; definitions given by Scargle, 1982, ApJ, 263, 835, and Horne ; and Baliunas, 1986, MNRAS, 302, 757. Beware of patterns and ; clustered data points as the Horne results break down in ; this case! Read and understand the papers and this ; code before using it! For the fast algorithm read W.H. Press ; and G.B. Rybicki 1989, ApJ 338, 277. ; ; The code is still stupid in the sense that it wants normal ; frequencies, but returns angular frequency... ; ; EXAMPLE: ; ; ; ; MODIFICATION HISTORY: ; Version 1.0, 1997, Joern Wilms IAAT ; Version 1.1, 1998.09.23, JW: Do not normalize if variance is 0 ; (for computation of LSP of window function...) ; Version 1.2, 1999.01.07, JW: force numf to be int ; Version 1.3, 1999.08.05, JW: added omega keyword ; Version 1.4, 1999.08 ; KP: significance levels ; JW: pmin,pmax keywords ; Version 1.5, 1999.08.27, JW: compute the significance levels ; from the horne number of independent frequencies, and not from ; numf ; Version 1.6, 2000.07.27, SS and SB: added fast algorithm and FAP ; according to white noise lc simulations. ; Version 1.7, 2000.07.28 JW: added debug keyword, sped up ; simulations by factor of four (use /slow to get old ; behavior of the simulations) ;- ;; defaults IF n_elements(noise) EQ 0 THEN noise = double(sqrt((moment(c))[1])) IF n_elements(multiple) EQ 0 THEN multiple = 0 IF n_elements(fap) EQ 0 THEN fap = 0.01 ;; make times manageable (Scargle periodogram is time-shift invariant) time = t-t[0] ;; number of independent frequencies ;; (Horne and Baliunas, eq. 13) n0 = n_elements(time) horne = long(-6.362+1.193*n0+0.00098*n0^2.) IF (horne LT 0) THEN horne=5 IF (n_elements(numf) EQ 0) THEN numf = horne ELSE horne = numf ;; min.freq is 1/T IF (n_elements(fmin) EQ 0) THEN BEGIN IF (n_elements(pmax) EQ 0) THEN BEGIN fmin = 1.D0 / max(time) END ELSE BEGIN fmin = 1.D0 / pmax END ENDIF ;; max. freq: approx. to Nyquist frequency IF (n_elements(fmax) EQ 0) THEN BEGIN IF (n_elements(pmin) EQ 0) THEN BEGIN fmax = n0 / (2.D0*max(time)) END ELSE BEGIN fmax = 1.D0 / pmin END ENDIF ;; if omega is not given, compute it IF (n_elements(omega) EQ 0) THEN BEGIN om = 2.D0 * !DPI* (fmin+(fmax-fmin)*findgen(numf)/(numf-1.D0)) END ELSE BEGIN om = omega END ;; False Alarm Probability according to Numf signi = -alog( 1.D0 - ((1.D0-fap)^(1./horne)) ) ;; Periodogram IF keyword_set(old) THEN BEGIN ;; slow version ;; Subtract mean from data cn = c-mean(c) ;; computing the periodogram px = fltarr(numf) FOR i=0L,numf-1L DO BEGIN tau = atan(total(sin(2.D0*om[i]*time))/total(cos(2.D0*om[i]*time))) tau = tau/(2.*om[i]) co = cos(om[i]*(time-tau)) si = sin(om[i]*(time-tau)) px[i]=0.5D0*(total(cn*co)^2/total(co^2)+total(cn*si)^2/total(si^2)) ENDFOR ;; correct normalization var = (moment(cn))[1] IF var NE 0. THEN BEGIN px = px/var END ELSE BEGIN print,'Scargle Warning: Variance is zero' END ;; some other nice helpers ;; computed here due to memory usage reasons nu = om/(2.D0*!dpi) period = 1.D0/nu return ENDIF ;; Ref.: W.H. Press and G.B. Rybicki, 1989, ApJ 338, 277 ;; Eq. (6); s2, c2 s2 = dblarr(numf) c2 = dblarr(numf) FOR i=0L,numf-1L DO BEGIN s2[i] = total( sin(2.D0*om[i]*time) ) c2[i] = total( cos(2.D0*om[i]*time) ) ENDFOR ;; Eq. (2): Definition -> tan(2omtau) ;; --- tan(2omtau) = s2 / c2 omtau = atan(s2/c2) / (2.D0) ;; cos(tau), sin(tau) cosomtau= cos(omtau) sinomtau= sin(omtau) ;; Eq. (7); total(cos(t-tau)^2) and total(sin(t-tau)^2) tmp = c2*cos(2.D0*omtau) + s2*sin(2.D0*omtau) tc2 = 0.5D0*(n0+tmp) ; total(cos(t-tau)^2) ts2 = 0.5D0*(n0-tmp) ; total(sin(t-tau)^2) ;; clean up tmp = 0. & omtau= 0. s2 = 0. & t2 = 0. ;; computing the periodogram for the original lc ;; Subtract mean from data cn = c - mean(c) ;; Eq. (5); sh and ch sh = dblarr(numf) ch = dblarr(numf) IF (multiple GT 0 AND NOT keyword_set(slow)) THEN BEGIN sisi=dblarr(n0,numf) coco=dblarr(n0,numf) FOR i=0,numf-1L DO BEGIN sisi[*,i]=sin(om[i]*time) coco[*,i]=cos(om[i]*time) sh[i]=total(cn*sisi[*,i]) ch[i]=total(cn*coco[*,i]) END END ELSE BEGIN FOR i=0L,numf-1L DO BEGIN sh[i] = total( cn*sin(om[i]*time) ) ch[i] = total( cn*cos(om[i]*time) ) ENDFOR END ;; Eq. (3) px = (ch*cosomtau + sh*sinomtau)^2 / tc2 + $ (sh*cosomtau - ch*sinomtau)^2 / ts2 ;; correct normalization px = 0.5D0*px/(noise^2) ;; --- RUN SIMULATIONS for multiple > 0 IF multiple GT 0 THEN BEGIN IF (multiple*min(fap) LT 10) THEN BEGIN message,'WARNING',/informational message,'Number of iterations (multiple keyword)',/informational message,'not large enough for false alarm probability',/informational message,'requested (need multiple*FAP > 10 )',/informational ENDIF IF (keyword_set(debug)) THEN BEGIN t0=systime(1) END psdpeak = dblarr(multiple) FOR m=0L,multiple-1L DO BEGIN IF (keyword_set(debug)) THEN BEGIN IF ((m+1) MOD 100 EQ 0) THEN BEGIN message,'...working on '+strtrim(m+1,2)+'th simul. ('+$ strtrim(string(format='(F8.3)',100.*m/multiple),2)+ $ '% done)',/informational exec=systime(1)-t0 message,' time expired : '+strtrim(exec,2)+'s',/info message,' time per step: '+ $ string(format='(F10.6)',exec/(m+1))+'s',/info message,' time remaining: '+ $ strtrim(exec/(m+1)*(multiple-m),2)+'s',/info END END ;; white noise simulation cn = randomn(anyseed,n0)*noise cn = cn-mean(cn) ;; .. force OBSERVED count rate to zero ;; Eq. (5); sh and ch IF (NOT keyword_set(slow)) THEN BEGIN FOR i=0L,numf-1L DO BEGIN sh[i]=total(cn*sisi[*,i]) ch[i]=total(cn*coco[*,i]) END END ELSE BEGIN FOR i=0L,numf-1L DO BEGIN sh[i] = total( cn * sin(om[i]*time) ) ch[i] = total( cn * cos(om[i]*time) ) ENDFOR END ;; Eq. (3) ; computing the periodogram for each simulation psdpeak[m] = max ( (ch*cosomtau + sh*sinomtau)^2 / tc2 + $ (sh*cosomtau - ch*sinomtau)^2 / ts2 ) ENDFOR ;; False Alarm Probability according to simulations IF n_elements(psdpeak) NE 0 THEN BEGIN idx = sort(psdpeak) ;; correct normalization psdpeaksort = 0.5D0 * psdpeak[idx]/(noise^2) simsigni = psdpeaksort[long((1.-fap)*(multiple-1))] ENDIF END ;; some other nice helpers ;; computed here due to memory usage reasons nu = om/(2.D0*!dpi) period = 1.D0/nu END