pro refsub_h1rg_stamp, input_data, output_data, fwindow=fwindow, order=order, degree=degree, plot=plot, simple=simple, verbose=verbose

       ; Name:
       ; refsub_H1RG.pro
       ;
       ; PURPOSE: 
       ; This program performs a reference pixel correction
       ; on Rockwell HAWAII-1RG detector data when read out using 2
       ; outputs in a 25x25 pixel subarray.
       ;
       ; CALLING SEQUENCE:
       ;       refsub_H1RG, input_data, output_data
       ;
       ; INPUTS:
       ;       input_data    - A 50x50 array containing the image data.
       ;       output_data   - A 50x50 fltarr for the result
       ;
       ; KEYWORD PARAMETERS:
       ;       fwindow       - The size of the filter window. 
       ;       order         - Filter order, see the SAVGOL help page
       ;       degree        - Filter degree, see the SAVGOL help page
       ;       plot          - Set this =1 to plot the fits
       ;       simple        - Set this =1 to subtract a median reference pixel
       ;                     
       ;       verbose       - Set this =1 for verbose mode.
       ;
       ;
       ; EXAMPLE
       ;        IDL> fits_read, 'theImage.fits', myData, myHeader
       ;        IDL> myNewDate = fltarr(1024, 1024)
       ;	IDL> refsub_H1RG, myData, myNewData
       ;
       ;
       ; REFERENCE:
       ;
       ;
       ; MODIFICATION HISTORY:
       ;       Written by:
       ;       Modified:    B.J. Rauscher, IDTL, August 4, 2002
       ;                    Initial release
       ;                    B.J. Rauscher, IDTL, August 16, 2002
       ;                    Modified to exclude outtermost reference pixels
       ;                    from statistics. These appear to be anomalously
       ;                    bright.
       ;                    B.J. Rauscher, IDTL, September 20, 2002
       ;                    Modified refsub_h1rg to get this program.
       ;
       ;-
	

       ; The reference pixels on the Rockwell HAWAII-1RG SCA are the
       ; outermost 4 pixels all around the chip. This program makes a
       ; reference pixel correction using data in the regions [1:3,*]
       ; and [46:48,*].

       ;; Set defaults
       if (not(keyword_set(plot))) then plot = 0 ; Default is no plots
       if (not(keyword_set(simple))) then simple = 0; Default is Savitzky-Golay filter
       if (not(keyword_set(verbose))) then verbose = 0; Default is quiet mode
       
       ;; We will use a Savitzky-Golay filter to smooth the reference pixel data.
       ;; See the SAVGOL procedure's help page for more information on
       ;; the following.
       if (not(keyword_set(fwindow))) then fwindow = 8   ; Default uses 8 pixels on
                                                         ; either side of current pixel
                                                         ; for smoothing
       if (not(keyword_set(order))) then order = 0; Use zero to smooth
       if (not(keyword_set(degree))) then degree=3; Three seems reasonable

       ;; Extract and average reference pixels
       row = intarr(50)            ; Y values vector
       mean_refcol_l = fltarr(50)  ; Signals vector on left
       mean_refcol_r = fltarr(50)  ; Signals vector on right
       for y=0, 49 do begin
           row[y] = y
           
           ;; Do not use the first reference pixel on each output. It
           ;; appears to be anomalous (contains signal, this may be memory).
           mean_refcol_l[y] = mean(input_data[1:3,y])       ; Average left side reference pixels
           mean_refcol_r[y] = mean(input_data[46:48,y]) ; Average right side reference pixels
       endfor

       ;; Sigma clip to get rid of outliers
       sigclip, mean_refcol_l, threshold=3.0, maxiter=5
       sigclip, mean_refcol_r, threshold=3.0, maxiter=5
       if (verbose) then begin
           print, 'Mean reference pixel on left = ', mean(input_data[0:3,*])
           print, 'Mean reference pixel on right = ', mean(input_data[46:48,*])
       endif

       ;; Apply reference pixel correction
       smooth_refcol_l = fltarr(50)
       smooth_refcol_r = fltarr(50)
       if (not simple) then begin
           ;; Apply Savitzky-Golay filter. See the SAVGOL help page for
           ;; more information.
           savgol_filter = savgol(fwindow, fwindow, order, degree)
           smooth_refcol_l = convol(mean_refcol_l, savgol_filter, /EDGE_TRUNCATE)
           smooth_refcol_r = convol(mean_refcol_r, savgol_filter, /EDGE_TRUNCATE)
       endif else begin
           ;; Apply simple reference pixel correction.
           for y=0, 49 do begin
               smooth_refcol_l[y] = mean(mean_refcol_l)
               smooth_refcol_r[y] = mean(mean_refcol_r)
           endfor
       endelse
               


       ;; Make plots if requested
       ;;
       if (plot) then begin

           set_device = 'X'
           !P.MULTI = [0,1,2]   ; Plot 1 plot over the other

           ;; First plot left hand reference pixels
           plot, row, mean_refcol_l, $
             title = 'Mean Reference Column [1:3,*]', $
             xtitle = 'Row Number', $
             ytitle = 'Signal (ADU)', $
             xrange = [0,49], $
             xstyle = 1

           ;; Draw a smooth curve through pixels
           oplot, row, smooth_refcol_l, thick=3

           ;; Now plot right hand reference pixels
           plot, row, mean_refcol_r, $
             title = 'Mean Reference Column [46:48,*]', $
             xtitle = 'Row Number', $
             ytitle = 'Signal (ADU)', $
             xrange = [0,49], $
             xstyle = 1

           ;; Draw a smooth curve through points.
           oplot, row, smooth_refcol_r, thick=3

           !P.MULTI = 0 ; Return to default plotting mode

       endif

       ;; Perform reference pixel correction
       for x=0, 24 do begin
           output_data[x,*] = input_data[x,*] - smooth_refcol_l
           output_data[49-x,*] = input_data[49-x,*] - smooth_refcol_r
       endfor
       
end