FUNCTION getIRSpectralMap, fitsHdr, DOGEOM=dogeom, FN=fn, GEOMFN=geomfn $
       , skipFunction=skipFunction, smileShift=smileShiftArg $
       , cacheSize=cacheSizeArg, cacheMiss=cacheMiss $
       , debug=debugArg
FORWARD_FUNCTION di_sqlquery, getSQLserver, getLocFn, dical_geom
;;-----------------------------------------------------------------------------
;; 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:
;;	/SKIPFUNCTION    - If set, do not use smile function, use spectral map file instead
;;                         - N.B. Default (SKIPFUNCTION not set) is to use smile function
;;                         - see comments for keyword FN below if set
;;	SMILESHIFT=-15.0 - Offset to use with smile function
;;                         - N.B. defaults to -15.0
;;	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
;;
;;   2011-05-12  B. Carcich     Added FM's smile shift per DDW 16.Mar, 2011 email and
;;                              code (GETIRSPECTRALMAPSA.PRO) and SP testing 2011-05-12
;;                              - smile shift defaults to -15.0
;;
;;-----------------------------------------------------------------------------

doDebug = keyword_set(debugArg) ? 1b : 0b
noDebug = 1b - doDebug
if noDebug then ON_ERROR, 2

common gismCache_cmn, tempArr, naxis1Arr, naxis2Arr, smileShiftArr, ptrMapArr, oldestIdx

;;;;;;;;;;;
;; 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 = keyword_set(skipFunction) ? 0b : 1b
	 
;; Get the mode
mode = sxpar(fitsHdr, 'IMGH030', COUNT=c1)

;; Get the temperature
temp = sxpar(fitsHdr, 'SMOBENT', COUNT=c2)
if (temp EQ -1 or c2 ne 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
	smileShift = n_elements(smileShiftArg) eq 1L ? float(smileShiftArg[0]) : -15.0
	dim = [sxpar(fitsHdr, 'NAXIS1'),sxpar(fitsHdr,'NAXIS2')]

	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;; Check cached results

	cacheSize = n_elements(cacheSizeArg) eq 1L ? long(cacheSizeArg[0]) : 0L

	if cacheSize gt 0L and n_elements(tempArr) gt 0L then begin
	  iw = where( tempArr eq temp $
	              and naxis1Arr eq dim[0] $
	              and naxis2Arr eq dim[1] $
	              and smileShiftArr eq smileShift $
	            , iwCount)
	  if iwCount gt 0L then begin
	    cacheMiss = 0L
	    return, *(ptrMapArr[iw[0]])
	  endif
	endif
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

	smile = -102.298 + 102.298 * sqrt(1.0 + ((indgen(512) + smileShift - 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

	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

	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;;; Save result (pointer to OUT) to cache ...
	if cacheSize gt 0L then begin
	  cacheMiss = 1L
	  ptrMap = ptr_new(out)
	  if n_elements(tempArr) eq 0L then begin
	    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	    ;;; ... cache is empty, save new ptr as first entry:
	    tempArr       = [temp]
	    naxis1Arr     = [dim[0]]
	    naxis2Arr     = [dim[1]]
	    smileShiftArr = [smileShift]
	    ptrMapArr     = [ptrMap]
	    oldestIdx = 0L
	  endif else if n_elements(tempArr) lt cacheSize then begin
	    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	    ;;; ... cache is shorter than allowed length, save new ptr to end
	    ;;; - recycle oldest ptr to beginning of array, reset index
	    if oldestIdx gt 0L then begin
	      tempArr       = shift( tempArr, -oldestIdx)
	      naxis1Arr     = shift( naxis1Arr, -oldestIdx)
	      naxis2Arr     = shift( naxis2Arr, -oldestIdx)
	      smileShiftArr = shift( smileShiftArr, -oldestIdx)
	      ptrMapArr     = shift( ptrMapArr, -oldestIdx)
	      oldestIdx = 0L
	    endif
	    tempArr       = [tempArr,temp]
	    naxis1Arr     = [naxis1Arr, dim[0]]
	    naxis2Arr     = [naxis2Arr, dim[1]]
	    smileShiftArr = [smileShiftArr, smileShift]
	    ptrMapArr     = [ptrMapArr, ptrMap]
	  endif else begin
	    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	    ;;; ... cache is equal or longer than allowed, save new ptr
	    ;;;       over oldest value
	    ;;; - increment oldest index, wrap modulo length-of-cache, 
	    ;;;     not by allowed length i.e. cache cannot shrink
	    tempArr[oldestIdx]       = temp
	    naxis1Arr[oldestIdx]     = dim[0]
	    naxis2Arr[oldestIdx]     = dim[1]
	    smileShiftArr[oldestIdx] = smileShift
	    ptrMapArr[oldestIdx]     = ptrMap
	    oldestIdx = (oldestIdx+1L) MOD n_elements(tempArr)
	  endelse
	  ;;;
	endif
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	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