;+
; NAME:
;        LOGREBIN
;
; PURPOSE:
;        Rebins a XYdY distribution into equal logarithmic (base 10)
;        intervals in X.
;
; CATEGORY:
;        Math.
;
; CALLING SEQUENCE:
;
;        LOGREBIN, X, Y, Dy, Ny, Nbins, DlogX
;
; INPUTS:
;        X:        Array of abcissae.
;
;        Y:        Array of Y values.
;
;        Dy:       Array of 1 sigma error bars along Y.
;
;        Ny:       Number of points used to determine each Dy, LONG.
;
;        Nbins:    Number of bins of the rebinned XYdY distribution.
;
; OUTPUTS:
;        X,Y,Dy:   The X-Y-dY distribution rebinned into Nbins equal
;                  logarithmic intervals in X, Array(Nbins)
;
;        Ny:       Number of points used to determine each rebinned Dy,
;                  LONARR(Nbins).
;
;        DlogX:    The width or bin size of each equal logarithmic
;                  interval in X.
;
; RESTRICTIONS:
;        Each input data point, X, Y, Dy is assumed to be derived from
;        averaging Ny distributions of XY.
;
; EXAMPLE:
;        ; Create some simple data from a SIN wave
;        X = FINDGEN(100)+50
;        Y = ABS(SIN(X/10))
;        dY= SQRT(Y/10)
;        ny= 5
;        nbins = 30
;
;        ; Plot original and rebinned XYdY distributions as log-scalar plots
;        ploterr,  X, Y, Dy, PSYM=3, /XTYPE, XRANGE=[MIN(X,MAX=mx),mx]
;        LOGREBIN, X, Y, Dy, Ny, Nbins, DlogX
;        oploterr, X, Y, Dy
;
; MODIFICATION HISTORY:
;        Written by:    Han Wen, August 1996.
;        19-AUG-1996    Corrected numerical round off errors and the case
;                       when Ny=1, leading to negative or NaNQ variances.
;                       Location of the last bin edge defined to be
;                       consistent with [min,max) convention.
;-
pro LOGREBIN, X, Y, Dy, Ny, Nbins, DlogX

;   Check integrity of input parameters

         NP   = N_PARAMS()
         if (NP lt 5) then message, $
              'Must be called with 5-6 parameters: X, Y, Dy, Ny, Nbins [,DlogX]'

         if (N_ELEMENTS(Ny) gt 1) then message, 'Ny parameter must be a scalar.'

         X0   = MIN(X,MAX=X1)
         if (X0 lt 0) then message, 'X parameter must be > 0.'

;   Define logarithmic interval, dlogX such that the last logarithmic bin
;   boundaries are [ logXmax+eps - dlogX, logXmax+eps ).

         logX = ALOG10(TEMPORARY(X))
         logX1= ALOG10(DOUBLE(X1))
         eps  = 10.^(FLOOR(ALOG10(logX1))-7)
         dlogX= ((logX1+eps) - ALOG10(DOUBLE(X0)))/Nbins

;   Rebin distribution with considerations for divide by 0s or round off
;   errors leading to negative variances, (>0)

         Y2   = Dy^2*(Ny-1) + Y^2
         hkeys= { BINSIZE: dlogX, MAX: logX1 }
         Ysum = HIST1D( logX, TEMPORARY(Y) , _EXTRA=hkeys, $
                        OBIN=ologX, DENSITY=nX )
         Y2sum= HIST1D( logX, TEMPORARY(Y2), _EXTRA=hkeys )

         nX_  = (nX > 1)               ; > 1 to avoid divide by 0s
         Yavg = Ysum /nX_
         Y2avg= Y2sum/nX_

         Ny   = nX*LONG(Ny)
         Ny_  = Ny > 2
         var  = (Ny/((Ny_-1.)))*( Y2avg - Yavg^2 ) > 0
         Dy   = SQRT( var/Ny_ )

         X    = 10L^ologX
         Y    = Yavg
end