;+
; NAME:
;     RANGE_DETECTOR
; PURPOSE:
;     Routine meant to automatically inspect spectrum data and divide
;     the frequency range in subranges obeying the same trend.
;     The algorithm computes an intrinsic fluctuation level and cut the data range in
;     subranges in wich data points change below the fluctuation level
; CATEGORY:
;     OVSA APC CALIBRATION ANALYSIS
; CALLING SEQUENCE:
;     range_detector,freq,flux [,lf=lf,uf=uf,pf=pf,pd=pd,$
;                  width=width,range=range,i_lf=i_lf,i_max=i_max,i_uf=i_uf,nosmooth=nosmooth,$
;                  goright=goright,goleft=goleft,first=first]
; INPUTS:freq, flux
; OUTPUTS: lf:low frequency array
;          uf:high frequency array
;          pf:peak frequency array
;          pd: peak flux array
;          each index correspond to a data subrange
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
;          width: minimum number of points in a data subrange
;          the rest of the parameters are meant to be used internally
;
; ROUTINES CALLED:
;
; OUTPUTS:
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 14-Aug-2002 by Gelu M. Nita (gn2@njit.edu)
;     04-Jun-2006  DG
;       IUF was always one less than maximum range due to FOR loop
;       only going to IEND-1.  Changed to go to IEND.
;
;-

pro range_detector,freq,flux,range=range,width=width,lf=lf,uf=uf,pf=pf,$
                  i_lf=i_lf,i_max=i_max,i_uf=i_uf,nosmooth=nosmooth,$
                  goright=goright,goleft=goleft,first=first,fluct=fluct



if n_elements(width) eq 0 then width=3


if not keyword_set(first) then begin
 if (not keyword_set(goleft)) and (not keyword_set(goright)) then begin
  goleft=1
  goright=1
 end
end

sd=flux
fd=freq
if n_elements(sd)lt width then return

if not keyword_set(nosmooth) then sd=smooth(sd,2)

if n_elements(range) ne 2 then begin
 istart=0
 iend=n_elements(sd)-1
 endif else begin
  istart=range[0]<range[1]
  iend=range[1]>range[0]
 end



;Compute alowed fluctuation
if n_elements(fluct) eq 0 then begin
ds=sd-shift(sd,1)
fluct=median(abs(ds/sd))
end


find_max:

ssmax=max(sd[istart:iend],imax,/nan)
imax=imax+istart

find_ilf:
smin=1

ILF=IMAX
if IMAX gt istart then begin
 CMIN =SSMAX
 for i= IMAX-1,istart,-1 do begin
  if (sd[i] lt smin) or sd[i] gt (1.0+fluct)*cmin then $
  goto, find_iuf else begin
   ILF=I
   if sd[i] lt cmin then cmin=sd[i]
  endelse
 endfor
endif

find_iuf:

IUF=IMAX
IF IMAX lt iend then begin
 cmin=ssmax
 for i=imax+1,iend do begin
  if sd[i] lt smin or  sd[i] gt (1.0+fluct)*cmin  then $
  goto, start_fit else begin
   IUF=i
   if sd[i] lt cmin then cmin=sd[i]
  endelse
 endfor
endif

start_fit:

if iuf-ilf ge width then begin
 if n_elements(i_max) eq 0 then begin
  i_max=[imax]
  i_lf=[ilf]
  i_uf=[iuf]
 endif else begin
  i_max=[i_max,imax]
  i_lf=[i_lf,ilf]
  i_uf=[i_uf,iuf]
 endelse

end

 ;Recursive calling
 if (ilf ge width) and keyword_set(goleft) then $
 range_detector,freq,flux,range=[0,ilf],width=width,i_lf=i_lf,i_max=i_max,i_uf=i_uf,/nosmooth,/goleft
 if (iuf+width le n_elements(flux)-1) and keyword_set(goright) then $
 range_detector,freq,flux,range=[iuf,n_elements(flux)-1],width=width,i_lf=i_lf,i_max=i_max,i_uf=i_uf,/nosmooth,/goright



if n_elements(i_lf) ne 0 then begin
 ;sort ranges in ascending order
 sorted=sort(i_lf)
 i_lf=i_lf(sorted)
 i_uf=i_uf(sorted)
 i_max=i_max(sorted)

 lf=fd[i_lf]
 uf=fd[i_uf]
 pf=fd[i_max]
end

end