;+
; NAME:
;	O_SKYSUBCORR
;
; PURPOSE:
;	fix incomplete sky subtraction from reduced cubes
;
; EXPLANATION:
;	if the sky subtraction isn't perfect, telluric correction is poor
;
; CALLING SEQUENCE:
;
;
; INPUTS:
;	xreg/yreg - 2 element vector describing sky region ex. [0,10]
;
; OUTPUTS:
;	fits file - does not overwrite original file
;
; OPTIONAL OUPUTS:
;	
;
; OPTIONAL INPUT KEYWORDS:
;       dir, filename, avgfunc
;
; EXAMPLE:
;	o_skysubcorr, xreg=[3,7], yreg=[10,54]
;
; ERROR HANDLING:
;	none
;
; RESTRICTIONS:
;	one file at a time
;
; NOTES:
;       
;
; PROCEDURES USED:
;	Functions:	o_make_filename
;                       fxaddpar
;
; MODIFICATION HISTORY:
;	2007 Jun 28 jlyke based upon o_telluric.pro
;	2007 Sep.  dlm: minor change of filedir to dir
;-
PRO O_SKYSUBCORR, dir=dir, file=file, $
              xreg=xreg, yreg=yreg, avgfunc=avgfunc

; set default name addition
add = 'ssc'

if ( not keyword_set(dir) ) then begin
    dir = './'
endif
if ( not keyword_set(file) ) then begin
    file = dialog_pickfile( title='File to fix sky subtraction', $
                                     filter='*.fits', $
                                     get_path=new_path, $
                                     /MUST_EXIST )
    ; separate the directory from the filename
    dir = new_path
    temp = strsplit(file, '/', /extract)
    filename = temp[n_elements(temp)-1]
endif else filename=file

; change logic here...only do something if the regions are defined
if ( keyword_set(xreg) or keyword_set(yreg) )  then begin
    if ( not keyword_set(avgfunc) ) then begin
        ; I think we should only do means...
        avgfunc='mean'
    endif else begin
        avgfunc=strlowcase(avgfunc)
    endelse

print, avgfunc
; region to calc dc level offset
;xreg = [0,10]
;yreg = [57,67]

    data = readfits(dir+filename, data_hdr)
    noise = readfits(dir+filename, exten_no = 1)
    qual = readfits(dir+filename, exten_no = 2)

    zeroth = n_elements(data[*,0,0])
    first  = n_elements(data[0,*,0])
    second = n_elements(data[0,0,*])

    dc = fltarr(zeroth)
    dcplot = fltarr(zeroth)

; define dc based upon avgfunc, but make the default "mean"

    if ( avgfunc eq 'median') then begin
        for i=0, zeroth-1 do begin
            dc[i] = median(data[i, yreg[0]:yreg[1], xreg[0]:xreg[1]])
        endfor
        plot, dc
        dclevel = median(dc)
        dcstdev = meanabsdev(dc, /median)
;dclevel = 0.0
    endif else begin
        avgfunc = 'mean'
        for i=0, zeroth-1 do begin
            dc[i] = mean(data[i, yreg[0]:yreg[1], xreg[0]:xreg[1]])
        endfor
        plot, dc
        dclevel = mean(dc)
        dcstdev = meanabsdev(dc)
;dclevel = 0.0
     endelse
        
     for i=0, zeroth-1 do begin
        dcplot[i] = dclevel
        data[i,*,*] = data[i,*,*] - dclevel
     endfor

; provide some output for the user
     print, "dclevel", dclevel
     print, "dcstdev", dcstdev
     
     plot, dc
     oplot, dcplot
     oplot, dcplot+dcstdev, linestyle=2
     oplot, dcplot-dcstdev, linestyle=2

; update the header info
    fxaddpar, data_hdr, 'DCLEVEL', dclevel, ' Bkg Sky Level Subtracted', $
	after='CUNIT1'
    fxaddpar, data_hdr, 'DCSTDEV', dcstdev, ' Bkg Sky Std Dev', $
	after='DCLEVEL'
    fxaddpar, data_hdr, 'AVGFUNC', avgfunc, ' Bkg Sky Avg Function', $
	after='DCSTDEV'
    ; make an xreg string
    xregion = '[' + strcompress(strjoin(string(xreg),','),/remove_all) + ']'
    fxaddpar, data_hdr, 'AVGXREG', xregion, ' Bkg Sky Avg X Region', $
	after='AVGFUNC'
    yregion = '[' + strcompress(strjoin(string(yreg),','),/remove_all) + ']'
    fxaddpar, data_hdr, 'AVGYREG', yregion, ' Bkg Sky Avg Y Region', $
	after='AVGXREG'

; create a new filename
    sscfile = o_make_filename(filename, add)

    writefits, dir+sscfile, data, data_hdr
    writefits, dir+sscfile, noise, /append
    writefits, dir+sscfile, qual, /append

; endif if for xreg/yreg keywords set
endif
END