;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; GETIRSPECTRALMAPSA.PRO from Dennis Wellnitz ca. 2011-03-16
;;; - N.B. GETIRSPECTRALMAPSA.PRO is not part of the operational pipeline
;;; - GETIRSPECTRALMAPSA.PRO shifts spectral map (smile) -15.0 pixels
;;; - GETIRSPECTRALMAP.PRO implements same shift with minimal changes to
;;;       previous repository version
;;; - Sylvia Protopapa tested GETIRSPECTRALMAPSA.PRO
;;; - See COMPAREGSM.PRO for comparison of the two routines' results
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION getIRSpectralMapSA, fitsHdr
;  , DOGEOM=dogeom, FN=fn, GEOMFN=geomfn
;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Calculatess the IR spectral map wavelengths and delta wavelengths
;;
;; 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).
;;
;; EXAMPLE:
;;      IDL> maps = getIRSpectralMap(fitsHdr)
;;
;;
;; 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-02-13  D. Wellnitz    Added FM's -15 pixel smile shift
;;-----------------------------------------------------------------------------

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 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

;; Use the local function to generate the map
	;; Create map of the smile offset
	;;smile = -102.298 + 102.298 * sqrt(1.0 + ((indgen(512) - 261.688)/740.984)^2)
	smile = -102.298 + 102.298 * sqrt(1.0 + ((indgen(512) - 15.0 - 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
	
END