;+
; NAME:
;   cmb_autobad
;            
; EXAMPLE USAGE:
;   i=cmb_autobad(data)
;
; PURPOSE:   
; To filter data returning an array of same size as data of 1's and 0's,
; 0 indicates bad data element and 1 indicates a good data element.
; Based on D. A. Roberts, "An algorithm for finding spurious points in turbulent signals.", 
; COMPUTERS IN PHYSICS, JOURNAL SECTION, SEPT/OCT 1993.
;
; CATEGORY:
; Data fitering.
;
; CALLING SEQUENCE:                                   
; result= cmb_autobad(datain,sigmul,fillval=fillval,nsum=nsum)                    
;                                                     
; INPUTS:                                             
;   datain - one or two dimensional data array of n n_elements: datain[2,n] or datain[n].
;   sigmul - multiplicative factor of standard devidation in each sub array: 5 (default),  4 (less agressive), 6 (more agressive).
;            Absolute values of the (data - mean) > sigmul*standard_deviation 
;            will be flagged as bad (0), if less than then flagged as good (1).
; Keyword Inputs:                                                      
;   fillval - fill values: -1e31 (default).
;   nsum - size of sub-array for each filter step: 100 (default).
;
; OUTPUTS:
;   Interger array of same size as the data array:
;     1 is a good data element, 
;     0 is a bad data element.
; COMMON BLOCKS:
;   None.
;
; SIDE EFFECTS:
;   Unknown.
;
; RESTRICTIONS:
;   Unknown.
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
;   Code developed by Aaron Roberts and Scott Boardsen at GSFC.
;-

function cmb_autobad0,datain,sigmul,fillval=fillval,nsum=nsum, nfill = nfill, filter_by_fill=filter_by_fill
;input
; datain n elements array of one dimension
; sigmul -multiplicative factor of the running standard deviation
; nsum - size of sub interval
; fillval = data fill value
;output
; igood n elements array of one dimension, 0 bad data point, 1 good data point

if n_elements(fillval) eq 0 then fillval=-1e31
if n_elements(sigmul) eq 0 then sigmul = 5
if n_elements(nsum) eq 0 then nsum=100l

y = datain
npts = n_elements(y)
igood = intarr(npts)


if keyword_set( filter_by_fill ) then begin
  ii=where(y ne fillval)
  nfill = n_elements(where(y eq fillval))
  if (where(y eq fillval))(0) eq -1 then nfill =0
endif else begin
  ii=where(abs(y) lt 0.99*abs( fillval))
  nfill = n_elements(where(abs(y) ge 0.99*abs( fillval)))
  if (where(abs(y) ge 0.99*abs( fillval)))(0) eq -1 then nfill =0
endelse
if ii[0] ne -1 then igood[ii] = 1
ig=where(igood eq 1)

if n_elements(ig) lt nsum  then begin
   print,'Not enough data points in series, only filtered for fill values, returning.'
   return,igood
endif

ng = n_elements(ig)
igmax = max(ig)

; made j0 long (Ron Yurow 5/5/21)
;j0 = nsum-1
j0 = LONG (nsum-1)
ii = ig[0:j0] ;initial sub-interval

icnt=0l ;counter, not used
; made iend long (Ron Yurow 5/5/21)
; iend=0 ;end of data indicator
iend=0l ;end of data indicator

while iend lt 2 do begin

      if min(y[ii]) eq fillval then stop
      ymean = total(y[ii])/nsum
      yvar = sqrt( total((y[ii]-ymean)^2)/(nsum-1) )
      jj=where( abs(y[ii]-ymean) gt sigmul*yvar) ;bad points in sub-interval ii

      ;set igood to 0 for bad data points, and reform sub-interval indices array for next step
      if jj[0] ne -1 then begin
         kk = ii[jj]
         ibad = n_elements(jj)
         if ibad eq nsum then begin
            ja = j0+1 < ng-1
            jb = j0+100 < ng-1
	    ii = ig[ja:jb]
	    j0 = j0 + 100   
         endif else begin
	    if jj[0] eq 0 then ibad=ibad-1
	    igood[kk]=0  ;set bad data point indices to 0
	    ii[jj] = -1  ;set sub internval indices of bad data points to -1
	    ii=ii[1:*]   ;step to next data point
	    ii=ii[where(ii ne -1)] ;remove indices of bad data points
	    j1 = j0 +ibad +1
	    if j1 ge ng then ii = ig[ng-100:ng-1] else  ii = [ii, ig[j0+1:j1]]
	    if n_elements(ii) ne nsum then stop
	    j0 = j1
         endelse
      endif else begin
         j0 = j0+1
         if j0 ge ng then ii = ig[ng-100:ng-1] else  ii = [ii[1:*], ig[j0]] ;reform sub-interval indices array for next step
      endelse
      if max(ii) ge igmax then  iend =iend+1 ;check if end of data is reached.
      icnt= icnt + 1l
endwhile

return,igood
end


function cmb_autobad,datain,sigmul,fillval=fillval,nsum=nsum,ibad=ibad,filter_by_fill=filter_by_fill
;datain[ic,itime] or datain [ic] either a 1 or 2 dimensional array

if n_elements(sigmul) eq 0 then sigmul = 5
if n_elements(nsum) eq 0 then nsum=100
print,'Filtering data using cmb_autobad: sigmul=',sigmul,' nsum=',nsum
si=size(datain,/str)
if si.n_dimensions eq 1 then begin
   ig= cmb_autobad0(datain,sigmul,fillval=fillval,nsum=nsum, nfill=nfill,filter_by_fill=filter_by_fill) 
   ibad=where(ig[*] eq 0)
   if ibad[0] eq -1 then ibad=0 else ibad=n_elements(ibad)-nfill
   print,' Number of fill points = ',nfill,' Percent fill=',nfill*1.0/(n_elements(ig))*100,'%' 
   print,' Number of bad points = ',ibad,' Percent bad=', ibad*1./(n_elements(ig)-nfill)*100,'%' 
endif else begin
ig = fix(datain)*0
for ic = 0l, si.dimensions[0]-1 do begin
    ig[ic,*] = cmb_autobad0(reform(datain[ic,*]),sigmul,fillval=fillval,nsum=nsum, nfill=nfill, filter_by_fill=filter_by_fill)
    ibad=where(ig[ic,*] eq 0)
    if ibad[0] eq -1 then ibad=0 else ibad=n_elements(ibad)-nfill
    print,' Number of fill points = ',nfill,' Percent fill=',nfill*1.0/(n_elements(ig[ic,*]))*100,'%'     
    print,'Component=',ic,' Number of bad points = ',ibad,' Percent bad=', ibad*1./(n_elements(ig[ic,*])-nfill)*100,'%' 
endfor
endelse
return,ig
end