;+
; NAME:
;        PSDREBIN
;
; PURPOSE:
;        Calculates the average power spectra from a given array of event times.
;
; CATEGORY:
;        Math.
;
; CALLING SEQUENCE:
;
;        PSDREBIN, Ts, Bin, Nbin, Timedel, PSDs
;
; INPUTS:
;        Ts:       Array of event times.
;
;        Bin:      The length of one time bin in the same UNITS as the Ts array.
;
;        Nbin:     The number of time bins in each time segment.
;
;        Timedel:  Time resolution or UNITS of Ts in SECONDS.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;
;        VERBOSE:  Setting the keyword will display messages updating the user
;                  on which time segment the routine is currently processing.
;
;        NLOGBIN:  Number of equal logarithmic bins to divide the frequencies
;                  into, (Default=100).
;
; OUTPUTS:
;        PSDs      Structure holding the average power spectras of the data
;                  and of expected Poisson statistics. Tags are defined as
;                  follows:
;
;                  freq      : Frequencies, [Hz]
;                  dfreq     : Width of the corresponding equal logarithmic
;                              binning of the frequencies and powers, [Hz]
;                  pwr       : average power spectra
;                  dpwr      : square root of variance of average
;                  params    : Structure holding the parameters used for
;                              analyzing the data and the poisson statistics:
;                              bin     : see Input definition above
;                              nbin    : "     "      "         "
;                              timedel : "     "      "         "
;                  mean      : Average number of counts/bin over all times FFTed.
;                  sigma:    : Sqrt of the variance of the counts/bin about
;                              the mean over all times FFTed.
;                  Leahy_coeff: Multiply this value with the average powers
;                              to obtain the Leahy normalized power spectra
;                  rms_coeff : Multiply this value with the average powers
;                              to obtain the fractional rms amplitude normalized
;                              power spectra
;
; EXAMPLE:
;        Consider an array of times, ts from the 4 microsec XTE data:
;
;        Timedel   = 3.814697265625d-6      ; seconds
;        Bin       = 13                     ; 1/(nbin*bin)   = Minimum frequency
;        Nbin      = 2L^19                  ; 1/(2*bin)      = Nyquist frequency
;
;        PSDREBIN, Ts, Bin, Nbin, Timedel, PSDs, /VERB
;        Freq      = PSDs.Freq
;        dFreq     = PSDs.dFreq
;        Pwr       = PSDs.pwr
;        dPwr      = PSDs.dpwr
;        PLeahy    = Pwr *PSDs.Leahy_coeff
;        DPLeahy   = dPwr*PSDs.Leahy_coeff
;        PLOTPSD, Freq, PLeahy, DFreq, DPleahy, PSYM=3
;
; MODIFICATION HISTORY:
;        Written by:    Han Wen, September 1996.
;-
pro PSDREBIN, Ts, Bin, Nbin, Timedel, PSDs, VERBOSE=Verbose, NLOGBIN=Nlogbin

         if (N_ELEMENTS(Nlogbin) eq 0) then Nlogbin = 100L

         if (KEYWORD_SET(Verbose)) then $
              print,'Calculating FFT Power Spectras of data...'
         rpsd = MULTIPSD( Ts, Bin, Nbin, Npsd, Mu, Sigma_mu, VERBOSE=Verbose )
         tpsd = Nbin*Bin*Timedel
         pwr  = rpsd(*,0)
         dpwr = rpsd(*,1)

         allp = TOTAL(pwr)             ; total unnormalized power
         avgI = Mu/(Bin*Timedel)       ; average intensity
         rms  = Sigma_mu/Mu            ; fractional rms amplitude
         norm1= rms^2*tpsd             ; rms   normalization
         Nrms = norm1/allp             ; rms   "coefficient"
         Nlhy = 2/(Nbin*Mu)            ; Leahy "coefficient"
         norm2= Nlhy*allp              ; Leahy normalization
         rerr = (norm1*avgI - norm2)/norm2
         if (rerr gt 0.01) then begin
         ;    print,'Normalizations rms: ',norm1,', Leahy: ',norm2,' avgI:',avgI
         ;    message,'Error in normalizations.'
         endif

         if (KEYWORD_SET(Verbose)) then $
              print,'Rebinning to equal logarithmic intervals...'
         freq =(FINDGEN(Nbin/2)+1)/tpsd
         LOGREBIN, freq, pwr, dpwr, Npsd, Nlogbin, dlogF

         readme    = [ 'Date: '+SYSTIME(),                                              $
                       'freq      : Frequencies, [Hz]                                 ',$
                       'dfreq     : Width of the corresponding equal logarithmic      ',$
                       '            binning of the frequencies and powers, [Hz]       ',$
                       'pwr       : average power spectra                             ',$
                       'dpwr      : square root of variance of average                ',$
                       'params    : Structure holding the parameters used for         ',$
                       '            analyzing the data:                               ',$
                       '            bin     : The length of one time bin in the same  ',$
                       '                      UNITS as the Ts array.                  ',$
                       '            nbin    : The number of time bins in each time    ',$
                       '                      segment.                                ',$
                       '            timedel : Time resolution or UNITS of Ts in       ',$
                       '                      SECONDS.                                ',$
                       'mean      : Average number of counts/bin over all times FFTed.',$
                       'sigma:    : Sqrt of the variance of the counts/bin about      ',$
                       '            the mean over all times FFTed.                    ',$
                       'Leahy_coeff: Multiply this value with the average powers      ',$
                       '            to obtain the Leahy normalized power spectra      ',$
                       'rms_coeff : Multiply this value with the average powers       ',$
                       '            to obtain the fractional rms amplitude normalized ',$
                       '            power spectra                                     ']
         readme    = TRANSPOSE(readme)
         dfreq     = 10^dlogf
         PSDs      = { $
                   readme    : readme, $
                   freq      : freq,   $
                   dfreq     : dfreq,  $
                   pwr       : pwr,    $
                   dpwr      : dpwr,   $
                   npts      : Npsd,   $
                   params    : { Bin   : Bin, $
                                 Nbin  : Nbin, $
                                 Timedel : Timedel }, $
                   Mean      : Mu, $
                   Sigma     : Sigma_Mu, $
                   Leahy_coeff : Nlhy, $
                   rms_coeff : Nrms  }
end