; NAME:
;     subspecsky
;
; PURPOSE:
;     To subtract the sky from a spectral image
;
; CATEGORY:
;     Data Reduction
;
; CALLING SEQUENCE:
;     subspecsky,images,var,edgecoeffs,norders,ybuffer,xranges,$
;                UPDATE=update,WIDGET_ID=widget_id,PLOT=plot,CANCEL=cancel
;
; INPUTS:
;     image      - An array [ncols,nrows ,nimages] of sky images
;     var        - The array [ncols,nrows ,nimages] of variance images
;     edgecoeffs - Array [degree+1,2,norders] of polynomial coefficients 
;                  which define the edges of the orders.  array[*,0,0]
;                  are the coefficients of the bottom edge of the
;                  first order and array[*,1,0] are the coefficients 
;                  of the top edge of the first order.
;     norders    - The number of orders 
;     ybuffer    - Number of pixels to buffer the edge of the orders
;     xranges    - An array [2,norders] of pixel positions where the
;                  orders are completely on the array
; 
; OUTUTS:
;     The images and var variables are modified in the program
;
; KEYWORD PARAMETERS:    
;     UPDATE    - If set, the program will launch the Fanning
;                 showprogress widget.
;     WIDGET_ID - If given, a cancel button is added to the Fanning
;                 showprogress routine.  The widget blocks the
;                 WIDGET_ID, and checks for for a user cancel
;                 command.
;     PLOT      - Set to plot the results as the program runs.
;     CANCEL    - Will be set on return if there is a problem.
;
; PROCEDURES CALLED:
;     Requires the Astronomy User's Library
;
;  PROCEDURE: 
;     The program loops over every order and uses the edgecoeffs to 
;     determine the top and bottom of the slit at every column.  The 
;     median sky level is then subtracted.
;
;REVISION HISTORY:
;     2000-11-19 - Written by M. Cushing, Institute for Astronomy, UH
;     2001-05-21 - Modified to propagate errors
;     2001-10-04 - Removed start and stop input and added xranges input 

pro subspecsky,images,var,edgecoeffs,norders,ybuffer,xranges,$
               UPDATE=update,WIDGET_ID=widget_id,PLOT=plot,CANCEL=cancel

cancel  = 0

;  Check parameters

if n_params() lt 6 then begin
    
    cancel = 1
    print, 'Syntax - subspecsky,images,var,edgecoeffs,norders,ybuffer,$'
    print, '                    xranges,UPDATE=update,$'
    print, '                    WIDGET_ID=widget_id,PLOT=plot,CANCEL=cancel'
    return

endif

zparcheck, 'subspecsky', images, 1, [2,3,4,5], [2,3], 'Data Cube'
zparcheck, 'subspecsky', var, 2, [2,3,4,5], [2,3], 'Variance cube'
zparcheck, 'subspecsky', edgecoeffs, 3, [2,3,4,5], [2,3], 'Edge coefficients'
zparcheck, 'subspecsky', norders, 4, [2,3,4,5], 0, 'Norders'
zparcheck, 'subspecsky', ybuffer, 5, [2,3,4,5], 0, 'Ybuffer'
zparcheck, 'subspecsky', xranges, 6, [2,3,4,5], [1,2], 'Xranges'


;  Get size of image and initialize x array.

images  = float(temporary(images))
s       = size(images)
ncols   = s[1]
nrows   = s[2]
nframes = (s[0] eq 3) ? s[3]:1

;  Set up the Fanning showprogress object if requested.

if keyword_set(UPDATE) then begin

    cancelbutton = (n_elements(WIDGET_ID) ne 0) ? 1:0
    progressbar = obj_new('SHOWPROGRESS',widget_id,color=2,$
                          CANCELBUTTON=cancelbutton,title='Spextool',$
                          MESSAGE='Subtracting the Sky Level...')
    progressbar -> start

endif

idx = 0
if keyword_set(PLOT) then begin

    window, /free
    win_idx = !d.window
    re = ' '

endif

for i = 0, nframes-1 do begin

;  Loop over the orders

    for j = 0,norders-1 do begin

        start = xranges[0,j]
        stop  = xranges[1,j]  
        x     = findgen(stop-start+1)+start
        
        botedge = round( poly(x,edgecoeffs[*,0,j]) )
        topedge = round( poly(x,edgecoeffs[*,1,j]) )
        for k = 0,n_elements(x)-1 do begin
            
;  Check to make sure the order is on the array and the column has
;  more than 2 data points.
            
            if botedge[k] gt ncols-1 then goto, cont1
            if topedge[k] lt 0       then goto, cont1

;  Clip the values to 0 < value < ncols-1
            
            bot = 0 > botedge[k]
            top = (ncols-1) < topedge[k]

            bot = bot+ybuffer
            top = top-ybuffer
            
            if top-bot lt 2 then goto, cont1
            
;  Compute background of the pixels in the slit and subtract.
            
            slit_flux = images[x[k],bot:top,i]
            slit_var  = var[x[k],bot:top,i]
            z         = where(findoutliers(slit_flux,3) eq 1,count)

            bg_flux  = (count eq 0) ? total(slit_flux/slit_var)/$
              total(1./slit_var):total(slit_flux[z]*1./slit_var[z])/$
              total(1./slit_var[z])
            bg_var   = (count eq 0) ? 1./total(1./slit_var):$
              1./total(1./slit_var[z])

            if keyword_set(PLOT) then begin
                
                xx = findgen(top-bot+1)+bot
                title='Frame '+string(i+1,format='(i2.2)')+', '+$
                  'Order '+string(j+1,format='(i2.2)')+', '+$
                  'Column '+strtrim(x[k],2)
                
                ptitle=title+' -- BG LEVEL = '+strcompress(bg_flux,/RE)+$
                  ' +- '+strcompress(sqrt(bg_var),/RE)
                plot, xx,images[x[k],bot:top,i],/xsty,/ysty,title=ptitle,$
                  psym=10
                plots, !x.crange,[bg_flux,bg_flux],color=2
                read, re
                
                images[x[k],bot:top,i] = images[x[k],bot:top,i]-bg_flux
                var[x[k],bot:top,i] = var[x[k],bot:top,i]+bg_var
                plot, xx,images[x[k],bot:top,i],/xsty,/ysty,title=ptitle,$
                  psym=10
                plots, !x.crange,[0,0],color=3
                plots, !x.crange,[bg_flux,bg_flux],color=2
                read, re
                
            endif else begin

                images[x[k],bot:top,i] = images[x[k],bot:top,i]-bg_flux
                var[x[k],bot:top,i] = var[x[k],bot:top,i]+bg_var
            
            endelse

        endfor
        cont1:
       
        idx = idx + 1
        if keyword_set(UPDATE) then begin
            
            if cancelbutton then begin
                
                cancel = progressBar->CheckCancel()
                if cancel then begin
                    
                    progressBar->Destroy
                    obj_destroy, progressbar
                    return
                    
                endif
                
            endif
            percent = (idx+1)*(100./float(norders*nframes))
            progressbar->update, percent
            
        endif

    endfor
    
endfor

if keyword_set(UPDATE) then begin

    progressbar-> destroy
    obj_destroy, progressbar
    
endif

if keyword_set(PLOT) then wdelete, win_idx

end