; flattopi.pro - find indices of the longest flattop
;		 Uses +/- 3% as a default
;
;   to test:	
;	IDL> data = [10,20,30,30,30,30,30,30,50, 30,30,30,20,10]
;	IDL> data = [10,20,30,30,50,30,30,30,30,50,30,30,30,30,30,30,20,10]
;	IDL> data = [10,20,30,30,30,30,30,30,50,30,30,30,30,50,30,30,20,10]
;	IDL> data = [10,20,30,30,30,50,30,30,30,30,30,30,30,50,30,30,20,10]
;	IDL> flattopi, data, imin, imax
;	IDL> print, data[imin:imax]

pro flattopi, data, imin, imax,  percent=percent

if n_elements( percent ) eq 0 then percent = 3.0

nbins = nint( 100./percent) < n_elements( data )

minData = MIN( data, MAX=maxData)
yTolerance = (maxData-minData)*percent/100.

bins = FINDGEN(N_ELEMENTS(data))/(N_ELEMENTS(data)-1)*(maxData-minData) + minData
BinSize = bins[1]-bins[0]
   ; compute x & y histogram values:
plothist, data, xhist, yhist, bin=binsize, /noplot

bins = xhist	; (already the same)
hist = yhist	; counts in corresponding bin

   ; find the bin with the most values
maxNum = MAX( hist, index )

   ; find the points in the band considered to be a flattop
flatInds = where( data gt bins[index]-yTolerance and  $
                  data le bins[index]+yTolerance )

   ; find where there are breaks in the flattop
deltaInds =  flatinds[1:*] - flatinds[0:*]
breakInds = where( deltaInds gt 1, breakCount )
if breakCount eq 1 then begin	; just one break, so take longest segment
   if breakInds[0] lt n_elements(deltaInds)/2 then begin
      iMin = flatinds[ breakInds[0]+1 ]
      iMax = flatinds[ n_elements(flatinds)-1 ]
   endif else begin
      iMin = flatInds[ 0 ]
      iMax = flatinds[ breakInds[0] ]
   endelse
   
endif else if breakCount gt 0 then begin

      ; several breaks in flattop found; find the longest without a break
   
   breakinds = [ breakInds, n_elements( deltaInds) ]	; add last segment

   distance = breakInds[1:*] - breakInds[0:*]
   distance = [ (breakInds[0]-0), distance ]		; add first segment

       ; find index of longest segment
   dum = max( distance, maxBrInd )

   if maxBrind eq 0 then iMin = flatinds[ 0 ] $
                    else iMin = flatinds[ breakinds[ maxbrind-1 ]+1]
   iMax = flatinds[ breakinds[ maxbrind ]]

endif else begin	; no breaks in flattop, so take the whole thing
   iMin = flatInds[0]
   iMax = flatInds[ n_elements(flatInds-1) ]
endelse

end