PRO fitcont, flux, nanValue=nanValue, continuum=continuum, rms=rms, filter=filter

y = flux
; Checking for NaN's and later on doing interpolation to
; replace them with linearly interpolated values taken over median
; values. All NaNs are treated by default if values are set to NaN. If
; not, the keyword nanValue shoud be set to the new nan value:
; ex. new=fitcont(old,0.0) ; where missing data are 0.0 in the input spec 
if (n_elements(nanValue)) then begin
   addNaNs=where(flux EQ nanValue,count)
   add=count
   nans=where(flux EQ nanValue OR (finite(flux,/NaN) NE 0),count)
   if (add ne 0) then y[addNaNs]=replicate(!Values.F_NAN, add)
endif else nans=where(finite(y,/NaN) NE 0,count) 
totNans=count

if (n_elements(filter) EQ 0) then filter=300 
m = n_elements(flux) 
if (m LT filter) then message, 'Filter larger then the spectrum!'
; Ycont will hold continuum values at the end
Ycont = flux
; yy is tmp array
yy = y
; Filter = 300px wide calculating median
; The loop consisting of two steps:
; 1) Filter will move from the first to the last pixel and in each step 
; median will be taken as continuum value (= median(300 points)
; / for the beginning and end of spectra where there are no fiter/2
; values on the side where there are enough points filter/2 will be
; taken and on the other side what is left. This is quite good
; solution, for it handles properly 'ending' of the spectrum. 

; 2) Additionaly 3 sigma clipping will be performed on these new values
; to exclude strong emission and absorption. So all values that are
; larger than 3*rms (rms=sqrt((new_continuum-prev_continuum)^2/m)) are
; replaced by linear interpolation of the nearest points that should
; be kept. 

; Both steps are repeated until one of the conditions are met: rms
; stopped decreasing or number of interation (=10) is excedeed.

niter=10 & iter=0
prevRMS=1. & rms=.5
x=indgen(m)
count = 0 ; to prevent iteration when all points should be excluded
while ((prevRMS/rms gt 1.) AND (iter le niter) AND (count ne m)) do begin 
 prevRMS = rms
 for i=0,m-2 do begin
 if (i ge (0+filter/2)) then begin 
        if (i le (m-1 - filter/2)) then $ 
           Ycont[i] = median(yy[(i-filter/2):(i+filter/2)]) $ 
           else Ycont[i]=median(yy[i-1-filter/2:*])
        endif else Ycont[i]=median(yy[0:i+filter/2])
    endfor ; END filter slidding 
  rms=sqrt(total((Ycont-yy)^2/m,/NaN))
  clip = where(yy ge (Ycont+3*rms) OR yy le (Ycont-3*rms),count,complement=keep)
; print,'iter=',iter,'  rms = ',rms,prevRMS/rms, ' eliminate', count
if (count ne 0 and count ne m) then yy=nneighbor(Ycont,clip,keep)
iter++ 
endwhile ; END 3 sigma-clipping

; Missing values y[nans] are simply
if (totNans ne 0) then begin 
   y[nans] = Ycont[nans]
   if (Finite(y[m-1],/NaN) EQ 0) then y[m-1]=Ycont[m-2] & Ycont[m-1]=Ycont[m-2]
endif

continuum=Ycont


END