FUNCTION getIRSpectralMap, fitsHdr, DOGEOM=dogeom, FN=fn, GEOMFN=geomfn
;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Determines the IR spectral map to use
;;
;; CALLING SEQUENCE:
;;	maps = getIRSpectralMap(fitsHdr, /GEOM, FN=fn, GEOMFN=geomfn)
;;
;; REQUIRED INPUTS:
;;	fitsHdr - The FITS header for the image being calibrated
;;
;; OUTPUTS:
;;	RETURN - Spectral maps the same size as the original image. maps[*,*,0]
;;		is the wavelength (um) and maps[*,*,1] is the delta wavelength (um).
;;
;; OPTIONAL INPUT KEYWORDS:
;;	DOGEOM - Flagged if the spectral map needs to be run through the 
;;		geometric calibration
;;	GEOMFN - Filename of the file used for geometric calibration
;;	FN - The filename of the spectral map to use. If this is not specified,
;;		one is interpolated using spectral maps in a database.
;;
;; EXAMPLE:
;;      IDL> maps = getIRSpectralMap(fitsHdr, /GEOM)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	dical_geom - Geometrically calibrates the spectral map if needed
;;
;; MODIFICATION HISTORY:
;;   2004-07-21  M. Desnoyer    Created
;;   2005-04-18  M. Desnoyer	Forced uniqueness on temperatures when retrieving
;;				maps from the server
;;   2005-06-17  M. Desnoyer    Added option to use Dennis' routine to generate maps
;;
;;-----------------------------------------------------------------------------

ON_ERROR, 2

;;;;;;;;;;;
;; This is a hack. I'd prefer having the routine be dynamically downloaded but there's
;; no time. So just call it directly and return. Leave the old sql code in for now
doFunc = 1	;; Flag to force the local function to be called instead of querying the database
	 
;; Get the mode
mode = sxpar(fitsHdr, 'IMGH030', COUNT=c1)

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

;; Get the date
date = sxpar(fitsHdr, 'OBSDATE', COUNT=c3)

;; Make sure we could get everything from the header
IF (c1 EQ 0) OR (c2 EQ 0) OR (c3 EQ 0) THEN $
	message, 'Invalid FITS header', /NONAME

;; See if we should use the local function to generate the map
IF doFunc THEN BEGIN
	;; Create map of the smile offset
	smile = -102.298 + 102.298 * sqrt(1.0 + ((indgen(512) - 261.688)/740.984)^2)
	smile = transpose(rebin(smile,512,1024))

	shift = 5
	map = rebin(indgen(1024),1024,512)	;; Routine can only produce the 1024 case
	map = map-shift+smile
	irrega, temp, map, wl, dwl

	dim = [sxpar(fitsHdr, 'NAXIS1'),sxpar(fitsHdr,'NAXIS2')]
	out = dblarr(dim[0],dim[1],2)

	;; Crop/bin as appropriate
	IF dim[0] EQ 1024 THEN BEGIN
		out[0,0,0] = wl
		out[0,0,1] = dwl	;; unbinned
	ENDIF ELSE BEGIN
		wl = rebin(wl,512,256)

		;; Binned, so summing is required in the delta wavelength
		even = indgen(512)*2
		dwl = dwl[even,*]+dwl[even+1,*]
		dwl = rebin(dwl,512,256)

		;; Crop as appropriate
		out[0,0,0] = wl[*,128-dim[1]/2:128+dim[1]/2-1]
		out[0,0,1] = dwl[*,128-dim[1]/2:128+dim[1]/2-1]
	ENDELSE
	RETURN, out
ENDIF
	
;; Get the spectral maps
IF not keyword_set(fn) THEN BEGIN
	;; Get the known maps with bracketing bench temperatures
	db = 'di_calib'
	tbl = 'SPEC_MAPS'
	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

	fnm = getLocFn(info[0,*], subdir='SPEC_MAPS')

	;; Open the files
	message, /reset
	IF foundlo THEN BEGIN
		lo = readFits(fnm[0], calHdr,/silent)
		lo = DICal_rotCalImg(lo,'HRIIR',calHdr)
	ENDIF

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

	;; 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, 'SPECMAP1', $
		strupcase(strmid(fnm[0], strpos(fnm[0], path_sep(), /reverse_search)+1))
	IF foundlo AND foundhi THEN $
		fxaddpar, fitsHdr, 'SPECMAP2', $
			strupcase(strmid(fnm[1], strpos(fnm[1], path_sep(), /reverse_search)+1))
	
ENDIF ELSE BEGIN
	out = readFITS(fn, /silent)
	IF !ERROR_STATE.CODE LT 0 THEN message, !ERROR_STATE.MSG, /NONAME

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

;; Geometrically calibrate the two maps
IF keyword_set(DOGEOM) THEN BEGIN
	out[*,*,0] = dical_geom(out[*,*,0], fitsHdr, fn=geomfn)
	out[*,*,1] = dical_geom(out[*,*,1], fitsHdr, fn=geomfn)
ENDIF
 
RETURN, out

END