FORWARD_FUNCTION DIcalIR_DNtoRad, getIRCalConstMap

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Converts the IR image from raw values to radience.
;;
;; CALLING SEQUENCE:
;;	out = DIcalIR_DNtoRad(in, fitsHdr)
;;
;; 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 in radience units
;;
;; OPTIONAL INPUT KEYWORDS:
;;	CONSTFN - Filename of the const map to use. If not specified, then
;;		this is found in a database
;;	SPECFN - Filename of the spectral maps to use. If not specified, then
;;		this is found in a database
;;	GEOMFN - Filename of the calibration data needed for geometric
;;		calibration of the maps. If not specified, then this is
;;		found in a database
;;	SPECMAP - Variable to contain the spectral map for this image
;;
;; EXAMPLE:
;;      IDL> outImg = DIcalIR_DNtoRad(imgIn, fitsHdr)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - get a FITS header item
;;	getIRCalConstMap - determine the calibration constant map
;;	getIRSpectralMap - determine the IR spectral map
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-04-12  M. Desnoyer    Reordered functions
;;   2005-04-15  B. Carcich     Preserve FITS header comments and identify 0 areas in
;;				calibration factors
;;   2005-05-10  M. Desnoyer    Force an output of 0 if calibration factor is 0
;;   2005-06-14  M. Desnoyer    Added ability to handle flight calibration constants
;;				which is a wavelength dependant lookup table
;;
;;-----------------------------------------------------------------------------

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Determines the IR absolute calibration to use. In flight, the format
;;	of this input changed. So, if MSNCONFIG == FLIGHT in the MAPHDR, then
;;	the output is a lookup table for a multiplicative factor. The key to the
;;	table is the wavelength.
;;
;; CALLING SEQUENCE:
;;	const = getIRCalConstMap(fitsHdr, mode, temp, date, FN=fn, MAPHDR=maphdr)
;;
;; REQUIRED INPUTS:
;;	fitsHdr - FITS header for the image
;;	mode - The image mode
;;	temp - The bench temperature of the image
;;	date - String specifying the date the image was taken
;;
;; OUTPUTS:
;;	RETURN - A map of ST calibration constants for the given binning and
;;		bench temperature
;;
;; OPTIONAL INPUT KEYWORDS:
;;	FN - Filename of the file containing the constant map. It must be
;;		stored as a FITS image.
;;	MAPHDR - Optional keyword to store the header of the constant map 
;;
;; EXAMPLE:
;;      IDL> calConst = getIRCalConstMap(fitsHdr, mode, temp, date)
;;
;; PROCEDURES USED (i.e. called directly!):
;;
;; MODIFICATION HISTORY:
;;   2004-07-21  M. Desnoyer    Created
;;   2005-04-19  M. Desnoyer    Changed temperature retrieval to not interpolate
;;				if the temperature is outside the range of valid maps
;;				found
;;   2005-06-14  M. Desnoyer    Added MAPHDR
;;
;;-----------------------------------------------------------------------------
FUNCTION getIRCalConstMap, fitsHdr, mode, temp, date, FN=fn, MAPHDR=maphdr

ON_ERROR, 2

;; Get the constant map
IF not keyword_set(fn) THEN BEGIN
	;; Get the known maps with bracketing bench temperatures
	db = 'di_calib'
	tbl = 'IR_ABSCAL'
	sel = ['Filepath', 'Temperature']
	condlo = 'Mode = '+strtrim(mode,2)+' AND Date <= "'+date+$
		'" AND Temperature <= '+strtrim(temp,2)+ $
		' order by Date desc, Version desc, Temperature desc  limit 1'
	condhi = 'Mode = '+strtrim(mode,2)+' AND Date <= "'+date+$
		'" AND Temperature >= '+strtrim(temp,2)+ $
		' order by Date desc, Version desc, Temperature limit 1'

	lo = di_sqlquery(getSQLserver(), db, tbl, sel, condlo, dim=dimlo)
	hi = di_sqlquery(getSQLserver(), db, tbl, sel, condhi, dim=dimhi)

	foundlo = min(dimlo) GT 0
	foundhi = min(dimhi) GT 0

	IF (NOT foundlo) AND (NOT foundhi) THEN $
		message, 'Could not Find Calibration File', /NONAME

	IF foundlo AND foundhi THEN info = [[lo],[hi]] $
	ELSE IF foundlo THEN info = lo $
	ELSE IF foundhi THEN info = hi

	fn1 = getLocFn(info[0,*], subdir='IR_ABSCAL')

	;; Open the files
	message, /reset
	IF foundlo THEN BEGIN
		lo = readFits(fn1[0], calHdr, /silent)
		lo = DICal_rotCalImg(lo,'HRIIR',calHdr)
	ENDIF
	IF foundhi THEN BEGIN
		hi = readFits(fn1[(foundlo)?1:0], calHdr, /silent)
		hi = DICal_rotCalImg(hi,'HRIIR',calHdr)
	ENDIF
	IF !ERROR_STATE.CODE LT 0 THEN message, !ERROR_STATE.MSG, /NONAME

	mapHdr = calHdr

	;; Interpolate linearily between the known maps
	IF foundlo AND foundhi THEN BEGIN
		info = info[*,uniq(info[0,*])]
		out = (hi-lo)/(fix(info[1,1])-fix(info[1,0]))*(temp-fix(info[1,0]))+lo
	ENDIF ELSE IF foundlo THEN out = lo $
	ELSE IF foundhi THEN out = hi

	;; Update the FITS header
	fxaddpar, fitsHdr, 'CALMAP1', $
		strupcase(strmid(fn1[0], strpos(fn1[0], path_sep(), /reverse_search)+1))
	IF foundlo AND foundhi THEN $
		fxaddpar, fitsHdr, 'CALMAP2', $
			strupcase(strmid(fn1[1], strpos(fn1[1], path_sep(), /reverse_search)+1))

ENDIF ELSE BEGIN
	out = readFITS(fn, mapHdr, /silent)
	IF !ERROR_STATE.CODE LT 0 THEN message, !ERROR_STATE.MSG, /NONAME

	;; Update the FITS header
	fxaddpar, fitsHdr, 'CALMAP1', $
		strupcase(strmid(fn, strpos(fn, path_sep(), /reverse_search)+1))
ENDELSE

RETURN, out

END

;;-----------------------------------------------------------------------------
;; Main function
;;-----------------------------------------------------------------------------
FUNCTION DIcalIR_DNtoRad, in, fitsHdr, CONSTFN=constfn, SPECFN=specfn, $
	GEOMFN=geomfn, SPECMAP=specmap

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

IF error NE 0 THEN BEGIN
	CATCH, /CANCEL
	message, 'Convert to Radience Units - ' + !ERROR_STATE.MSG, /NONAME
ENDIF

;; Get the integration time in milliseconds
intTime = SXPAR(fitsHdr, 'INTTIME', COUNT=c1)

;; Get the mode
mode = SXPAR(fitsHdr, 'IMGMODE', COUNT=c2)

;; Get the bench temperature
temp = sxpar(fitsHdr, 'SMOBENT', COUNT=c3)
if (temp EQ -1) THEN temp = sxpar(fitsHdr, 'OPTBENT', COUNT=c3)

;; Get the image's units
units = SXPAR(fitsHdr, 'BUNIT', COUNT=c4)

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

;; Make sure all the header items exist
IF (c1 EQ 0) OR (c2 EQ 0) OR (c3 EQ 0) OR (c4 EQ 0) OR (c5 EQ 0) THEN $
	message, 'Invalid FITS Header', /NONAME

;; Check the units of the file to see if this can be applied
IF strtrim(units,2) EQ 'W/[m^2 sr um]' THEN RETURN, in ;; This has already been applied
IF strtrim(units,2) NE 'DATA NUMBER' THEN $
	message, 'Image must be in DN units to convert to Radiance'

out = double(in)

;; Determine if the image has been geometrically calibrated and if so,
;; apply the same calibration to the spectral map
geom = SXPAR(fitsHdr, 'CALGEOM', COUNT=c1)
IF c1 EQ 0 THEN geom = 0

;; Get the spectral maps. First dimension is wavelength [um], second dimension is delta.wavelength [um]
specMap = getIRSpectralMap(fitsHdr, DOGEOM=geom, FN=specfn, GEOMFN=geomfn)

;; Convert the integration time to seconds
intTime = double(intTime)*1d-3

;; Identify the proper calibration constant (QET) based on the mode. (e/photons)
ST = getIRCalConstMap(fitsHdr, mode, temp, date, FN=constfn, MAPHDR=maphdr)

;; The in flight calibration constant map is different, so handle the change
IF strtrim(sxpar(maphdr,'MSNCONFG'),2) EQ 'FLIGHT' THEN BEGIN
	;; In this case, the constant map is actually a lookup table with three rows
	;; the first row is the wavelength in um. The second row is the constant for
	;; areas not under the anti-saturation filter. The third row is the constant for
	;; areas under the anti-saturation filter. The location of the anti-saturation
	;; filter is specified in the data's header. In all cases, the constants are in 
	;; units of [W/(m^2 sr)]/(DN/s)

	;; Get the location of the anti-saturation filter
	antiSatLoc = intarr(2)
	antiSatLoc[0] = sxpar(mapHdr,'ASATST',count=c1)
	antiSatLoc[1] = sxpar(mapHdr,'ASATEND',count=c2)
	IF c1 EQ 0 OR c2 EQ 0 THEN message, 'Invalid calibration constant table'

	;; Make the calibration constant map
	dim = size(out,/dimensions)
	fac = dblarr(dim[0],dim[1])
	IF antiSatLoc[0] GT 0 THEN fac[0,0] = $
		interpol(ST[*,1],ST[*,0],specMap[*,0:antiSatLoc[0]-1,0])
	fac[0,antiSatLoc[0]] = $
		interpol(ST[*,2],ST[*,0],specMap[*,antiSatLoc[0]:antiSatLoc[1],0])
	IF antiSatLoc[1] LT dim[1]-1 THEN fac[0,antiSatLoc[1]+1] = $
		interpol(ST[*,1],ST[*,0],specMap[*,antiSatLoc[1]+1:*,0])

	;; The deltawavelength in binned data is twice as large, but the data is averaged
	;; so, we need have a factor of two to counteract this in binned data
	binFac = sxpar(fitsHdr, 'IMGMODE')
	binFac = (binFac LE 3 OR binFac EQ 5)+1

	;; Fill in the rest of the calibration factor. Note that what was just calculated
	;; is a multiplicative factor whereas the value needed below is a divisional
	;; factor 
	fac = specMap[*,*,1]*(intTime/binFac)/fac
ENDIF ELSE BEGIN
		
	;; Identify the second constant we need (gOmegaA) (m^2 sr DN / electron)
	const = 1.177d-13

	;; Get the conversion factor from W-s to photons for all wavelengths (photons/[Ws])
	P = specMap[*,*,0]/1.98648d-19

	fac=(intTime*const*P*ST*specMap[*,*,1])

	fxaddpar, fitsHdr, 'CALCONST', const
ENDELSE

;; Change the Physical unit in the FITS header

fxaddpar, fitsHdr, 'BUNIT', 'W/[m^2 sr um]'
FXADDPAR, fitsHdr, 'CALCORR', 'T'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;BTCarcich2005-04-13 START
;;;
;;; Calculate calibration factors for entire array
;;; - Check for ZERO calibration factors
;;; - Put count and envelope of ZERO calibration factors into FITS header
;;;   if there's anything bad

iw=where(fac eq 0d0,iwcount,complement=good)
if iwcount gt 0 then begin
  dims=size(st,/dim)
  dim0=dims[0]
  xmax = max(iw mod dim0, min=xmin)
  ymax = max(iw / dim0, min=ymin)
  xextent = strjoin( strtrim([xmin,xmax],2), ':')
  yextent = strjoin( strtrim([ymin,ymax],2), ':')
  zeroextent = strjoin( [ '[', strjoin([xextent,yextent],','), ']'], '')
  out[iw] = 0
endif


;; Identify bad locations in the FITS header
fxaddpar, fitsHdr, 'CALNUMB0', iwcount

;; Perform the unit conversion
out[good] = out[good]/fac[good]
RETURN, out
;;;BTCarcich2005-04-13 END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

END