FUNCTION DIcal_flagGaps, in, fitsHdr, FLAGS=flags

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Sets all pixels whose DN value is 0 to NaN. ie. flags the 
;;	gaps in data with NaN. This is because we should only get 0 DN value
;;	if the data got lost somwhere
;;	
;;	!! Warning !!! - This routine is obsolete because gaps are earlier
;;		on the SDC when the images are reconstructed from the original
;;		byte stream
;;
;; CALLING SEQUENCE:
;;	out = DIcal_flagGaps(in)
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The fits header for this 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 gap flags. This new array can be assigned to
;;		a variable
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_flagGaps(imgIn)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - Reads an item in the FITS header
;;
;; MODIFICATION HISTORY:
;;   2004-06-02  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------

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

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

;; Identify where the image is
out = double(in)

;; Get the detector
det = SXPAR(fitsHdr, 'IMGH021', COUNT=c1)

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

;; Set up the flag map
dim = size(out, /dimensions)
IF n_elements(flags) eq 0 THEN flags = intarr(dim[0],dim[1])

@dical_flags

;; Flag the gaps with NaN
;; The easiest type of gaps to find are where we have a string of 85's 1017 
;; bytes long. This occurs on a CFDP data gap. (I assume that means 
;; communications). Also, data comes down row by row with words rotating 
;; through the quadrants in order. Because of the difficulty of reforming
;; the original byte stream from the image, we're just going to break up
;; the image into quadrants and look for strings of 85's that are 254-255 bytes
;; long (VIS) and 508-509 (IR)
IF det EQ 0 THEN BEGIN ;; VIS
	quadCnt = 4
	err = getBIAScols(fitsHdr, COUNT=xshift)
	err = getPOCrows(fitsHdr, COUNT=yshift)
	xshift = xshift/2 & yshift = yshift/2
	rotates = [7,2,0,5]
	q = dblarr(dim[0]/2,dim[1]/2,quadCnt,/NOZERO)
	q[*,*,0] = out[0:dim[0]/2-1, dim[1]/2:*]
	q[*,*,1] = out[dim[0]/2:*, dim[1]/2:*]
	q[*,*,2] = out[0:dim[0]/2-1, 0:dim[1]/2-1]
	q[*,*,3] = out[dim[0]/2:*, 0:dim[1]/2-1]
ENDIF ELSE IF det EQ 1 THEN BEGIN
	quadCnt = 2
	xshift = 0
	yshift = 0
	rotates = [0,5]
	q = dblarr(dim[0]/2,dim[1],quadCnt,/NOZERO)
	q[*,*,0] = out[0:dim[0]/2-1,*]
	q[*,*,1] = out[dim[0]/2:*,*]
ENDIF ELSE message, 'Invalid detector', /NONAME

;; Rearrage the quadrants so they are in the order data was received
FOR i=0, quadCnt-1 DO $
	q[*,*,i] = shift(rotate(q[*,*,i],rotates[i]),-xshift,-yshift)
q = reform(q,dim[0]*dim[1]/quadCnt,quadCnt,/overwrite)

;; Identify the gaps
gaps = where(q EQ 21845,cnt)
IF cnt EQ 0 THEN goto, ret

;; Find the strings of 85 bytes
i=1
j=0
first = gaps[0]
last = first
WHILE i LT cnt DO BEGIN
	IF gaps[i] NE last[j] + 1 THEN BEGIN ;; Found the end of a gap
		first = [first,gaps[i]]
		last = [last,gaps[i]]
		j = j+1
	ENDIF ELSE last[j] = gaps[i]
	i = i+1
ENDWHILE

;; Make sure the strings are the right size
bigGaps = where((last-first LE 508/quadCnt) AND $
	(last-first GE 508/quadCnt-2),cnt)
IF cnt EQ 0 THEN GOTO, ret
first=first[bigGaps]
last=last[bigGaps]

;; Search around the big gaps for single 85 bytes in another number
n = (size(first,/dimensions))[0]-1
FOR i =0, n DO BEGIN
	low = first[i]-1
	IF low GE 0 THEN BEGIN
		IF uint(q[low]) MOD 256 EQ 85 OR uint(q[low])/256 EQ 85 THEN $
			first[i] = low
	ENDIF
	high = (last[i]+1)<(n_elements(q)-1)
	IF high LT n_elements(q) THEN BEGIN
		IF uint(q[high]) MOD 256 EQ 85 OR uint(q[high])/256 EQ 85 THEN $
			last[i] = high
	ENDIF
	
	;; Fill the gap with NaNs
	q[first[i]:last[i]] = !VALUES.D_NAN
ENDFOR

;; Rebuild the image
q = reform(q, dim[0]/2, dim[1]/(quadCnt/2), quadCnt,/OVERWRITE)
FOR i=0, quadCnt-1 DO $
	q[*,*,i] = rotate(shift(q[*,*,i],xshift,yshift),rotates[i])

IF det EQ 0 THEN BEGIN ;;VIS
	out[0:dim[0]/2-1, dim[1]/2:*] = q[*,*,0]
	out[dim[0]/2:*, dim[1]/2:*] = q[*,*,1]
	out[0:dim[0]/2-1, 0:dim[1]/2-1] = q[*,*,2]
	out[dim[0]/2:*, 0:dim[1]/2-1] = q[*,*,3]
ENDIF ELSE BEGIN;;IR
	out[0:dim[0]/2-1,*] = q[*,*,0]
	out[dim[0]/2:*,*] = q[*,*,1]
ENDELSE

;; Add the gaps to the flag map
w = where(finite(out) EQ 0, cnt)
IF cnt GT 0 THEN flags[w] = flags[w] OR FLAG_GAP

ret:

;; Update the FITS header
fxaddpar, fitsHdr, 'GAPFL', 'T'

RETURN, out
END