PRO DI_writefits_flagmap, fn, flag

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Writes a map of bit flags to a FITS file as an extension. Since the 
;;	flag array should be sparse, the data is transformed into an array
;;	of long integers where the most significant 12 bits are the flags
;;	and the least significant 20 bits specify the the location of that
;;	flag in the array as a 1 dimensional index.
;;	As of 5/24/05, the array is no longer compressed and only 8 bits of flags are
;;	allowed
;;
;; CALLING SEQUENCE:
;;	DI_writefits_flagmap, fn, flag
;;
;; REQUIRED INPUTS:
;;	fn - filename of the FITS file to extend the flag map onto
;;	flag - A MxN array of flag values, one for each pixel, to store
;;
;; OPTIONAL INPUT KEYWORDS:
;;	none
;;
;; EXAMPLE:
;;      IDL> DI_writefits_flagmap, fn, flag
;;
;; PROCEDURES USED (i.e. called directly!):
;;    writeFITS, sxaddhist, sxaddpar, mkhdr
;;
;; MODIFICATION HISTORY:
;;   2004-08-18  M. Desnoyer    Created
;;   2005-05-24  M. Desnoyer	Hey, it's my birthday. Yay! Anyway, I removed
;;				the space saving compression because it's very 
;;				complicated to get the headers right for PDS etc.
;;				so now, it's just a byte array the same size as
;;				the original image.
;;   2005-06-02  B. Carcich	grammar (The word "data" is plural)
;;
;;-----------------------------------------------------------------------------

on_error, 2

;; Get the dimensions of the flag map
dim = size(flag, /dimensions)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;; M. Desnoyer 05/24/05
;; Compress the map
;;w = where(flag NE 0, cnt)
;;IF cnt NE 0 THEN comp = long(w) or ishft(long(flag[w]),20) ELSE comp = [0]
;;host_to_ieee, comp
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
flag = byte(flag)

;; create the header for the extension
mkhdr, hdr, flag, /image
sxdelpar, hdr, 'O_BZERO'
sxaddpar, hdr, 'BZERO', 0, ' Data zero location'
sxaddpar, hdr, 'BSCALE', 1, ' Data scaling factor'
sxaddpar, hdr, 'BUNIT', 'N/A', ' Data are unitless'

;; determine the comment

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;; 05/24/05 M. Desnoyer
;;com = 'This extension contains a compressed format of a flag map. '+$
;;	'It is an array where each element represents a non-zero entry of the map.'+$
;;	' The 12 most significant bits specify the flag, while the 20 least '+$
;;	'significant bits specify the one dimensional index of that flag '+$
;;	'in the map.'
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

com = 'This extension contains a flag map for the data. Attached to every pixel '+$
	'is a set of flags, one per bit. These flags contain information about the '+$
	'quality of the data in the pixel.'
wrds = strsplit(com)
wrds = [wrds, max(wrds)+strlen(strmid(com,max(wrds)))+1]
cur = 0
lines = ''
WHILE max(cur LT wrds) EQ 1 DO BEGIN
	next = wrds[max(where(wrds LT cur+67))]
	lines = [lines, '** '+strmid(com, cur, next-cur-1)]
	cur = next
ENDWHILE
com = [$
	'**********************************************************************', $
	lines[1:*], $
	'**********************************************************************'$
	]
sxaddhist, com, hdr, /comment
sxaddpar, hdr, 'EXTNAME', 'FLAGS', ' Unique name for this image extension'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;; 05/24/05 M. Desnoyer
;;sxaddpar, hdr, 'MAPAXIS1', dim[0], ' Horizontal size of the map'
;;sxaddpar, hdr, 'MAPAXIS2', dim[1], ' Verical size of the map'
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Write the extension
writefits, fn, flag, hdr, /append

END