FUNCTION DIcal_flagSat, in, fitsHdr, FLAGS=flags

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Adds a saturation flag to the flag array for every pixel whose DN 
;;	value is >= satVal. There are three types of saturation levels that are
;;	flagged: 
;;		1. Some pixels fullwell
;;		2. Most pixels fullwell
;;              3. ADC Saturation
;;      And by definition, the levels are 1 < 2 < 3
;;
;; CALLING SEQUENCE:
;;	out = DIcal_flagSat(in, fitsHdr, FLAGS=flags)
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The FITS header for the image
;;
;; OUTPUTS:
;;	RETURN - the image after going through this calibration step
;;
;; OPTIONAL INPUT KEYWORDS:
;;	FLAGS - Array the same size as 'in' that contains flags for pixels. If
;;		This is not specified, a new array is created and populated
;;		with saturation flags. This new array can be assigned to
;;		a variable
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_flagSat(imgIn, fitsHdr, FLAGS=flags)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - Reads the FITS header
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-04-05  B. Carcich     Converted satVal array to LONG from STRING
;;
;;-----------------------------------------------------------------------------

;; Includes
@dical_flags

;; Do error handling 
;CATCH, error
error=0

IF error NE 0 THEN BEGIN
	CATCH, /CANCEL
	message, 'Flag Saturated Pixels - ' + !ERROR_STATE.MSG, /NONAME
ENDIF

dim = size(in, /dimensions)

;; Get the instrument
inst = SXPAR(fitsHdr, 'INSTRUME', COUNT=c1)

;; Get the date
date = SXPAR(fitsHdr, 'OBSDATE', COUNT=c2)

;; Make sure all items were in the header
IF c1 EQ 0 or c2 EQ 0 THEN message, 'Invalid FITS Header', /NONAME

;; Check the flag array
IF n_elements(flags) EQ 0 THEN flags = intarr(dim[0], dim[1])

;; Get the values for the various saturation levels
db = 'di_calib'
tbl = 'SATVALS'
sel = ['SomeFull', 'MostFull', 'ADC']
cond = "Date <= '"+date+"' AND Instrument='"+inst+$
	"' order by Date desc, Version desc limit 1"

satVals = di_sqlquery(getSQLServer(), db, tbl, sel, cond, dim=dim)
IF min(dim) EQ 0 THEN $
	message, 'Could not find valid saturation levels', /noname

;; Convert to integer
satVals = long(satVals)

;; Flag the saturated pixels
index = where(in GE satVals[0],cnt)
IF cnt NE 0 THEN flags[index] = flags[index] OR FLAG_SOMEFULLWELL
index = where(in GE satVals[1],cnt)
IF cnt NE 0 THEN flags[index] = flags[index] OR FLAG_MOSTFULLWELL
index = where(in GE satVals[2],cnt)
IF cnt NE 0 THEN flags[index] = flags[index] OR FLAG_ADCSAT

;; Update the FITS header
fxaddpar, fitsHdr, 'SATPIX', 'T'
fxaddpar, fitsHdr, 'SMSATVL', satVals[0]
fxaddpar, fitsHdr, 'MSTSATVL', satVals[1]
fxaddpar, fitsHdr, 'ADCSATVL', satVals[2]

RETURN, in
END