pro sigclip, data_in, threshold=threshold, fwindow, order, degree, help=help, maxiter=maxiter, verbose=verbose, stddevin=stddevin, plot=plot

        ;NAME:
        ;     pro sigclip.pro (version 2.1)
        ;
        ; PURPOSE:
        ;     This program fits a smooth and continuous function
        ;     to input data using a Savitsky-Golay filter and sigma
        ;     clipping. Upon exit, the input data are overwritten with
        ;     the smoothed data.
        ;
        ; CALLING SEQUENCE:
        ;     sigclip, data_in [, threshold=threshold, /help]
        ;
        ; INPUTS:
        ;     data_in   - a 1-dimensional fltarr to be smoothed
        ;     threshold - The sigma clipping threshold
        ;     fwindow   - Window size for SAVGOL filter
        ;     order     - Order for SAVGOL filter
        ;     degree    - Degree for SAVGOL filter
        ;     maxiter   - THIS DOES NOTHING. IT IS HERE FOR COMPATIBILITY WITH OLDER SOFTWARE
        ;     verbose   - THIS DOES NOTHING. IT IS HERE FOR COMPATIBILITY WITH OLDER SOFTWARE
        ;     stddevin  - THIS DOES NOTHING. IT IS HERE FOR COMPATIBILITY WITH OLDER SOFTWARE
        ;
        ; OUTPUTS
        ;
        ; KEYWORD PARAMETERS:
        ;      help      - Print a help message
        ;      plot      - Make a plot to show the fitting
        ;
        ; EXAMPLE
        ;
        ;
        ;
        ; REFERENCE:
        ;
        ;       This program has been approved for release by
        ;       B.J. Rauscher. It contains no proprietary information.
        ;
        ;
        ; MODIFICATION HISTORY:
        ;       Written by:
        ;       Modified:    B.J. Rauscher, March 18, 2003
        ;                    Initial release
        ;
        ;       4/8/03 -- B.J. Rauscher, Use interpolation to fill
        ;                 "gaps" of outliers. Use only one iterations
        ;                 since this seems to suffice.
        ;  -


       ;; Handle request for help
       if (not keyword_set(help)) then begin
           help = 0
       endif else begin
           stop, 'Usage: sigclip, data_in [, threshold=threshold, fwindow=fwindow, order=order, degree=degree, maxiter=maxiter, /verbose, /plot, /help]'
       endelse

       ;; Set defaults
       if (not(keyword_set(threshold))) then threshold = 5.0; Default sigma clipping threshold =3 sigma
       if (not(keyword_set(fwindow))) then fwindow = 64
       if (not(keyword_set(order))) then order = 0
       if (not(keyword_set(degree))) then degree = 3
       if (not(keyword_set(plot))) then plot = 0; Default does not make a plot

       ;; Create the Savitsky-Golay kernel
       savgol_filter = savgol(fwindow, fwindow, order, degree)

       ;; Do the sigma clipping. It seems like only 1 iteration is
       ;; needed to get essentially everything.
       for niter = 0, 1 do begin

           ;; For the 0th iteration, use simple sigma clipping to
           ;; identify outliers. This works better in practice because
           ;; the first or last row in an SCA can occasionally be
           ;; anomalous.
           if (niter eq 0) then begin

               djs_iterstat, data_in, mean=mean, sigrej=3, mask=mask

           endif else begin

               ;; For subsequent iterations, apply a Savitsky-Golay filter and calculate residuals against this.

               ;; Smooth the data
               data_out = convol(reform(data_in), savgol_filter, /EDGE_TRUNCATE)

               ;; Calculate a residual array
               residual = sqrt((data_in - data_out)^2)

               ;; Identify statistical outliers. Use 3 sigma clipping
               djs_iterstat, residual, sigma=rms_error, sigrej=3, mask=mask

           endelse

           ;; The first and last pixels must be good. If they
           ;; are not, replace them with an average of neighboring pixels

           ;; Handle pixel 0
           if (mask[0] eq 0) then begin
               djs_iterstat, data_in[0:fwindow/2-1], mean = meanpx, sigrej=3
               data_in[0] = meanpx
               mask[0] = 1
           endif

           ;; Handle last pixel
           if (mask[n_elements(data_in)-1] eq 0) then begin
               djs_iterstat, data_in[n_elements(data_in)-fwindow/2:n_elements(data_in)-1], mean = meanpx
               data_in[n_elements(data_in)-1] = meanpx
               mask[n_elements(data_in)-1] = 1
           endif

           ;; Use linear interpolation to fill in outliers.
           ;; First look for statistical outliers
           for j = 0L, n_elements(data_in) - 2 do begin
               if (mask[j] eq 0) then begin
                   ;; Find the next "good" pixel
                   for k = j + 1, n_elements(data_in) - 1 do begin
                       if (mask[k] eq 1) then begin
                           slope = (data_in[k] - data_in[j-1])/(k - j + 1)
                           for u = j, k - 1 do data_in[u] = data_in[j-1] + slope * (u - j + 1)
                           j = k
                           break
                       endif
                   endfor
               endif
           endfor

       endfor

       if (plot) then begin

           plot, data_in, yrange=[-10,10], psym=1
           oplot, data_out, thick=3, color=2

       endif

       ;; Overwrite the input data with the result
       data_in = data_out

   end