;+
; NAME:
;        FASPER
;
; PURPOSE:
;        Given 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 creates an array wk1 with a
;        sequence of nout increasing frequencies (not angular frequencies)
;        up to hifac times the "average" Nyquist frequency, and creates
;        an array wk2 with the values of the Lomb normalized periodogram at
;        those frequencies.  The arrays x and y are not altered.  This
;        routine also returns jmax such that wk2(jmax) is the maximum
;        element in wk2, 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.
;
; CATEGORY:
;        Numerical Recipes routines.
;
; CALLING SEQUENCE:
;
;        FASPER, X, Y, Ofac, Hifac, Wk1, Wk2 [, Nout, Jmax, Prob]
;
; INPUTS:
;        X:       Abscissas array, (e.g. an array of times).
;
;        Y:       Ordinates array, (e.g. an array of corresponding counts).
;
;        Ofac:     Oversampling factor.
;
;        Hifac:    Hifac * "average" Nyquist frequency = highest frequency
;                  for which values of the Lomb normalized periodogram will
;                  be calculated.
;
; OUTPUTS:
;        Wk1:      An array of Lomb periodogram frequencies.
;
;        Wk2:      An array of corresponding values of the Lomb periodogram.
;
; OPTIONAL OUTPUTS:
;        Nout:     The dimension of Wk1 and Wk2, i.e. the number of frequencies
;                  calculated.
;
;        Jmax:     The array index corresponding to the MAX( Wk2 ).
;
;        Prob:     The False Alarm Probability (FAP) of the largest value of
;                  of Lomb normalized periodogram.
;
; MODIFICATION HISTORY:
;        Written by:    Han Wen, December 1994. (Adapted from a Numerical
;                       Recipes routine with the same name).
;        15-AUG-1996    Andrew Lee, Modified spread and fasper to use arrays
;                       starting at 0 instead of 1, and fixed some bugs where
;                       int was used instead of long.
;-


pro SPREAD, y, yy, n, x, m
;
;   Given an array yy(0:n-1), extirpolate (spread) a value y into
;   m actual array elements that best approximate the "fictional"
;   (i.e., possible noninteger) array element number x.  The weights
;   used are coefficients of the Lagrange interpolating polynomial

         nfac=[0,1,1,2,6,24,120,720,5040,40320,362880]

         if (m gt 10) then message,'factorial table too small in spread'
         ix=long(x)
         if x eq float(ix) then yy(ix)=yy(ix)+y $
         else begin

              ilo= ( long(x-0.5*m+1.0) > 0 ) < (n-m)
              ihi=ilo+m-1
              nden=nfac(m)
              fac=x-ilo
              for j=ilo+1,ihi do fac = fac*(x-j)
              yy(ihi) = yy(ihi) + y*fac/(nden*(x-ihi))
              for j=ihi-1,ilo,-1 do begin
                   nden=(nden/(j+1-ilo))*(j-ihi)
                   yy(j) = yy(j) + y*fac/(nden*(x-j))
              endfor
         endelse
end
; (C) Copr. 1986-92 Numerical Recipes Software


pro FASPER, x, y, ofac, hifac, wk1, wk2, nout, jmax, prob

         MACC = 4  ;Number of interpolation points per 1/4 cycle
                   ;of highest frequency

;   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
         nfreqt=long(ofac*hifac*n*MACC)     ;Size the FFT as next power
         nfreq=64L                          ;of 2 above nfreqt.

         while (nfreq lt nfreqt) do nfreq = ISHFT( nfreq,1 )
         ndim=long(ISHFT( nfreq,1 ))

         var=(stdev(y,ave))^2.              ;Compute the mean, variance
                                            ;and range of the data.
         xmin=MIN( x, MAX=xmax )
         xdif=xmax-xmin

         wk1 =complexarr(ndim)            ;Extirpolate the data into
         wk2 =complexarr(ndim)            ;the workspaces.
         fac=ndim/(xdif*ofac)
         fndim=ndim
         ck  =(x-xmin)*fac MOD fndim
         ckk =2.0*ck MOD fndim

         for j=0L,n - 1 do begin
              SPREAD, y(j)-ave,wk1,ndim,ck(j),MACC
              SPREAD, 1.0,wk2,ndim,ckk(j),MACC
         endfor

         wk1  = FFT( wk1,1, /OVERWRITE )    ;Take the Fast Fourier Transforms.
         wk2  = FFT( wk2,1, /OVERWRITE )

         wk1  = wk1(1:nout)
         wk2  = wk2(1:nout)
         rwk1 = double( wk1 ) & iwk1 = imaginary( wk1 )
         rwk2 = double( wk2 ) & iwk2 = imaginary( wk2 )
         df=1.0/(xdif*ofac)

         hypo2 = 2.0 * abs( wk2 )                  ;Compute the Lomb value for each
         hc2wt= rwk2/hypo2               ;frequency.
         hs2wt= iwk2/hypo2

         cwt  = sqrt(0.5+hc2wt)
         swt  = SIGN(sqrt(0.5-hc2wt),hs2wt)
         den  = 0.5*n+hc2wt*rwk2+hs2wt*iwk2
         cterm= (cwt*rwk1+swt*iwk1)^2./den
         sterm= (cwt*iwk1-swt*rwk1)^2./(n-den)

         wk1  = df*(findgen(nout)+1.)
         wk2  = (cterm+sterm)/(2.0*var)
         pmax = MAX( wk2, jmax )

         expy =exp(-pmax)                   ;Estimate significance of largest
         effm =2.0*(nout)/ofac              ;peak value.
         prob =effm*expy
         if (prob gt 0.01) then prob=1.0-(1.0-expy)^effm

end
; (C) Copr. 1986-92 Numerical Recipes Software