;+ ; NAME: ; PERIOD ; ; PURPOSE: ; Given data points with abscissas x (which need not be equally spaced) ; and ordinates y, and given a desired oversampling factor ofac (a ; typical value being 4 or larger), this routine fills array px ; with an increasing sequence of frequencies (not angular frequencies) ; up to hifac times the "average" Nyquist frequency, and fills array ; py with the values of the Lomb normalized periodogram at those ; frequencies. The arrays x and y are not altered. The routine also ; returns jmax such that py(jmax) is the maximum element in py, and ; prob, an estimate of the significance of that maximum against the ; hypothesis of random noise. A small value of prob indicates that ; a significant periodic signal is present. ; ; (Adapted from routine of the same routine in Numerical Recipes in C, ; Second edition). ; ; CATEGORY: ; Math. ; ; CALLING SEQUENCE: ; ; PERIOD, X, Y, Ofac, Hifac, Px, Py, Nout, Jmax, Prob ; ; INPUTS: ; X: The abscissas of the data points. ; ; Y: The ordinates of the data points. ; ; Ofac: Oversampling factor. ; ; Hifac: Maximum frequency of Px = Hifac * "average" Nyquist ; frequency,. ; ; OUTPUTS: ; Px: The frequencies of the Lomb normalized periodogram. ; ; Py: The spectral powers of the Lomb normalized periodogram. ; ; Nout: The number of elements of the Px or Py array. ; ; Jmax: Index to the maximum element in Py. ; ; Prob: Significance of Py(Jmax) against the hypothesis of ; random noise. ; ; MODIFICATION HISTORY: ; Written by: Han Wen, August 1996. ; 07-AUG-1996 Eliminated redundant variables. ;- pro PERIOD, x, y, ofac, hifac, px, py, nout, jmax, prob ; Check dimensions of input arrays n = N_ELEMENTS( x ) sz = SIZE( y ) if n ne sz(1) then message, 'Incompatible arrays.' nout=0.5*ofac*hifac*n var=(STDEV(y,ave))^2. ;Compute the mean, variance ;and range of the data. xmin=MIN( x, MAX=xmax ) xdif=xmax-xmin xave=0.5*(xmax+xmin) pymax=0.0 pnow =1.0/(xdif*ofac) ;Starting frequency arg = 2.0*!dpi*(x-xave)*pnow ;Initialize values for the wpr = -2.0*SIN(0.5*arg)^2 ;trigonometric recurrences at each wpi = SIN(arg) ;data point. The recurrences are wr = COS(arg) ;done in double precision wi = wpi px = pnow + FINDGEN(nout)/(ofac*xdif) py = FLTARR(nout) yy = y-ave for i=0L,nout-1 do begin ;Main loop over the frequencies to ;be evaluated sumsh= TOTAL(wi*wr) ;First, loop over the data to get sumc = TOTAL((wr-wi)*(wr+wi)) ;tau and related quantities. wtau =0.5*ATAN(2.0*sumsh,sumc) swtau=SIN(wtau) cwtau=COS(wtau) ss =wi*cwtau-wr*swtau ;Get the periodogram value. cc =wr*cwtau+wi*swtau sums =TOTAL(ss^2) sumc =TOTAL(cc^2) sumsy=TOTAL(yy*ss) sumcy=TOTAL(yy*cc) wtemp=wr wr =(wtemp*wpr-wi*wpi)+wr ;Update the trigonometric recurrences wi =wi*wpr+wtemp*wpi+wi py(i)=0.5*(sumcy^2/sumc+sumsy^2/sums)/var endfor pymax=MAX(py,jmax) expy =EXP(-pymax) effm =2.0*nout/ofac prob =effm*expy if prob gt 0.01 then prob=1.0-(1.0-expy)^effm end