fORWARD_FUNCTION DIcal_darkFrame, insertIRDarkModelParams,getModelIRDark, getModelVISDark, getSpecificDark, getBias, $
	dark_frame_generate

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Ideally uses a dark frame stored in the database for the given observation.
;;	However, if there is no known frame, then a model is used to calculate
;;	the dark frame based on the temperatures and other parameters of the
;;	observation. If the USEFIT keyword is set, then the model is used 
;;	automatically without first checking for a good dark frame.
;;
;;	Another option is to manually supply the dark frame. This is done by
;;	setting the DARK to an array the same size as the original image.
;;
;; CALLING SEQUENCE:
;;	RESULT = DIcal_darkFrame( in, fitsHdr, flags, /USEFIT, DARKFN=darkfn)
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The FITS header for the image
;;	flags - an image map specifying the flags for the pixels
;;
;; OUTPUTS:
;;	RETURN - the image after going through this calibration step
;;
;; OPTIONAL INPUT KEYWORDS:
;;	USEFIT - Use the fit model parameter to create the interpolation
;;	DARKFN - Fileanem of the dark frame to subtract. If this is not set,
;;		a frame is interpolated
;; EXAMPLE:
;;      IDL> out=DIcal_darkFrame(imgIn, fitsHdr, flags, /USEFIT)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - Read the FITS header
;;	getModelIRDark - Calculate a dark frame based on a model. In this file.
;;	getModelVISDark - Calculate a dark frame based on a model. In this file.
;;	getInterpDark - Calculates a dark frame by interpolating. In this file.
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-04-12  M. Desnoyer    Reordered functions
;;
;;-----------------------------------------------------------------------------

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Puts IR model dark parameters into the calibrated heder
;;
;; CALLING SEQUENCE:
;;	insertIRDarkModelParams, calHdr, darkHdr, mstFn
;;
;; REQUIRED INPUTS:
;;	calHdr - The header for the calibrated file
;;	darkHdr - The header for the dark frame
;;	mstFn - Filename of the master dark used
;;
;; OUTPUTS:
;;	None
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> insertIRDarkModelParams, calHdr, darkHdr
;;
;; PROCEDURES USED (i.e. called directly!):
;;	sxaddpar
;;
;; MODIFICATION HISTORY:
;;   2005-11-07  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------
pro insertIRDarkModelParams, calHdr, darkHdr, mstFn


;; Convert from keywords in darkHdr to the ones in the final header
;;	darkHdr		final
hdrConvert = [ $
	['MDARKVER',	'MDARKVER'], $
	['PREVMODE', 	'DRKPMODE'], $
	['TIMEISG',	'DARKISG'], $
	['COEFF_A0', 	'DARKA0'], $
	['COEFF_A1', 	'DARKA1'], $
	['COEFF_B0', 	'DARKB0'], $
	['COEFF_B1', 	'DARKB1'], $
	['COEFF_C0', 	'DARKC0'], $
	['TIMESCAL',	'DRKTMSCL'] $
]

;; Cycle through all of the values to change
n = (size(hdrConvert, /dimensions))[1]
FOR i=0, n-1 DO BEGIN
	tmp = sxpar([darkHdr], hdrConvert[0,i],count=nmatch)
	if nmatch eq 1 then fxaddpar, calHdr, hdrConvert[1,i], tmp[0]
ENDFOR

;; Identify whether the dark plateau was reached
tmp = sxpar([darkHdr], 'FLAG-DRK', count=n)
tmp = (n eq 0)?'T':'F'
fxaddpar, calHdr, 'DARKPLAT', tmp

;; Specify the master dark used
fxaddpar, calHdr, 'DARKFN', $
	strupcase(strmid(mstFn, strpos(mstFn, path_sep(), /reverse_search)+1))
return

end

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Creates a dark frame based on a previously calculated model.
;;
;; CALLING SEQUENCE:
;;	RESULT = getModelVISDark(inst, date, dTemp, bTemp, intT, mode, hdr)
;;
;; REQUIRED INPUTS:
;;	inst - string: the instrument for the dark needed (e.g. HRIVIS)
;;	date - string: the latest date in ISO format of the model to use
;;	dTemp - Detector temperature
;;	bTemp - Bench temperature
;;	intT - Integration time
;;	mode - Image mode
;;	hdr - The FITS header for the image
;;
;; OUTPUTS:
;;	RETURN - the dark frame to subtract from the image
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> dark=getModelVisDark(det, inst, date, dTemp, bTemp, intT, mode, hdr)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	DI_sqlQuery - Query the SQL server
;;	getSQLServer - Retrives info on how to acess the SQL server
;;	getLocFn - Converts a remote filename into a local one
;;
;; MODIFICATION HISTORY:
;;   2004-07-19  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------
FUNCTION getModelVISDark, inst, date, dTemp, bTemp, intT, mode, hdr

;; Check the inputs
IF dTemp LT 0 OR bTemp LT 0 THEN $
	message, 'Invalid Temperature: Cannot compute dark frame'

;; Parameters to acess SQL database
serv = getSQLServer()
db = 'di_calib'

;; Get the filename for the model parameters from the database
sel = ['Filepath']
cond = "(Instrument='" + inst + "') AND (Date<='" + date + $
	"') AND (Mode="+strtrim(mode,2)+") order by Date DESC, Version desc limit 1"
tblNm = 'DARK_MODEL'
fn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)
IF max(dim) EQ 0 THEN $
	message, 'Database contained no dark model',/NONAME

;; Convert that filename into a local one
fn = getLocFn(fn, subdir='DARK_MODEL')

;; Load the model parameter coefficients into IDL. Loads a cube called
;; A. First two dimensions are the image dimensions. Last dimension
;; is the number of parameters in the model
message, /reset
A = readfits(fn[0],/silent)
IF !error_state.code LT 0 THEN message, 'Invalid dark parameters', /nonam

;; Create the dark frame based on the model
eg = A[*,*,1] - (7.021d-4*dTemp^2)/(1108.+dTemp)
dark = intT*1.12455d10*dTemp^1.5*exp(-5.802252562d3/dTemp*eg)*A[*,*,0]

;; Update the FITS header
fxaddpar, hdr, 'DARKFN', $
	strupcase(strmid(fn[0], strpos(fn[0], path_sep(), /reverse_search)+1))

RETURN, dark
END


;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Creates an IR dark frame based on a previously calculated model.
;;
;; CALLING SEQUENCE:
;;	RESULT = getModelIRDark(fitsHdr,img)
;;
;; REQUIRED INPUTS:
;;	fitsHdr - FITS header for the image. Assumes it is valid
;;	img - the original image
;;
;; OUTPUTS:
;;	RETURN - the dark frame to subtract from the image
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> dark=getModelIRDark(fitsHdr,img)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	DI_sqlQuery - Query the SQL server
;;	getSQLServer - Retrives info on how to acess the SQL server
;;	getLocFn - Converts a remote filename into a local one
;;	di_file_copy - Copies a file
;;
;; MODIFICATION HISTORY:
;;   2004-12-08  M. Desnoyer    Created
;;   2005-06-16  M. Desnoyer    Added img
;;
;;-----------------------------------------------------------------------------
FUNCTION getModelIRDark, fitsHdr, img

;; Get the routine to generate the darks

;; Parameters to access database
serv = getSQLServer()
db = 'di_calib'

;; Get the filename of the procedure to download
sel = ['Filepath']
cond = "Date <= '"+strtrim(sxpar(fitsHdr,'OBSDATE'),2)+"' order by Date DESC, Version desc limit 1"
tblNm = 'DARKIR_PRO'
fn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)
IF max(dim) EQ 0 THEN $
	message, 'Database contained no valid dark model procedure',/NONAME

;; Get the local filename of the code
fn = getLocFn(fn, subdir='DARKIR_PRO')

;; Copy the procedure into the working directory
di_file_copy, fn[0], "dark_frame_generate.pro", /f

;; Compile the procedure
resolve_routine, 'dark_frame_generate', /is_function, /compile_full_file

;; Now that we have the procedure to run, get the data it needs from the header
curMode = sxpar(fitsHdr, 'IMGMODE')
fpaTemp = sxpar(fitsHdr, 'IRFPAT')
covTemp = sxpar(fitsHdr, 'COVERT')
priMirTemp = sxpar(fitsHdr, 'PRIMIRT')
secMirTemp = sxpar(fitsHdr, 'SECMIRT')
intTime = sxpar(fitsHdr, 'INTTIME')
curSlot = sxpar(fitsHdr, 'FNIMGNMB')
simTemp = sxpar(fitsHdr, 'SMOBENT',count=c)
if (simTemp LT 0) THEN simTemp = sxpar(fitsHdr, 'OPTBENT')

;; Check the inputs
IF fpaTemp LE 0 OR simTemp LE 0 OR covTemp LE 0 OR priMirTemp LE 0 OR $
	secMirTemp LE 0 THEN $
	Message, 'Invalid temperatures in header: Dark cannot be calculated"

;; Now get the data from the last sequence.
;; Get the observation time of the first image in this sequence
serv = getSQLServer()
db = 'disdc'
sel = ['OBSJD']
cond = '(EXPID = '+strtrim(sxpar(fitsHdr, 'EXPID'),2)+$
	") AND (INSTRUME='HRIIR') AND (FNIMGNMB=1) AND (MSNCONFG='"+$
	sxpar(fitsHdr, 'MSNCONFG')+"') AND (OBSJD<="+$
	string(sxpar(fitsHdr,'OBSJD'),format='(D13.7)')+$
	") order by OBSJD DESC limit 1"
tblNm = 'HRIIR'
seqStart = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)


IF max(dim) EQ 0 THEN BEGIN
	;; Counldn't find a valid image so just set the parameters to dummy values
	seqGap = -999
	lastMode = -999
ENDIF ELSE BEGIN
	
	seqStart = seqStart[0]

	;; Figure out when the last sequence ended and what mode it was 
	serv = getSQLServer()
	db = 'disdc'
	sel = ['OBSENDJD','IMGMODE']
	cond = "(MSNCONFG = '"+sxpar(fitsHdr, 'MSNCONFG')+$
		"') AND (INSTRUME='HRIIR') AND (OBSENDJD < "+$
		seqStart+") order by OBSENDJD DESC limit 1"
	tblNm = 'HRIIR'
	data = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)

	;; Check results from the SQL server
	IF max(dim) EQ 0 THEN BEGIN
		;; Bad header so just insert dummy values
		seqGap = -999
		lastMode = -999
	ENDIF ELSE BEGIN
		;; Calculate the inter-sequence gap
		seqGap = double(seqStart)-double(data[0,0])

		;; Convert the sequence gap from days to miliseconds
		seqGap = seqGap * 86400000

		;; Get the last mode
		lastMode = fix(data[1,0])
	ENDELSE
ENDELSE

;; Get the master dark frame
serv = getSQLServer()
db = 'di_calib'
sel = ['Filepath']
cond = "Date <= '"+strtrim(sxpar(fitsHdr,'OBSDATE'),2)+"'" $
      +" AND instrument='" + strtrim(sxpar(fitsHdr,'INSTRUME'),2) + "'" $
      +" AND Mode = "+strtrim(curMode,2)+$
	" order by Date DESC, Version desc limit 1"
tblNm = 'DARK_MODEL'
fn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)
IF max(dim) EQ 0 THEN $
	message, 'Database contained no valid master dark',/NONAME

;; Get the local filename of the master dark
fn = getLocFn(fn, subdir='DARK_MODEL')

;; Generate the dark frame
out = dark_frame_generate(curMode,fpaTemp,simTemp,covTemp,priMirTemp,secMirTemp, $
	seqGap,intTime,lastMode,curSlot, img, darkHdr, MDARK_FN=fn[0])

;; Update the FITS header
insertIRDarkModelParams, fitsHdr, darkHdr, fn[0]

RETURN, out

END

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Retrieves the third frame in the scan to use as a dark. Also linearizes it
;;
;; CALLING SEQUENCE:
;;	RESULT = getThirdIRDark(fitsHdr,img)
;;
;; REQUIRED INPUTS:
;;	fitsHdr - FITS header for the image. Assumes it is valid
;;	img - the original image
;;
;; OUTPUTS:
;;	RETURN - the dark frame to subtract from the image
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> dark=getThirdIRDark(fitsHdr, img)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	getLocFn - Converts a remote filename into a local one
;;
;; MODIFICATION HISTORY:
;;   2005-07-03 M. Desnoyer	Created
;;
;;-----------------------------------------------------------------------------
FUNCTION getThirdIRDark, fitsHdr, img

;; Get the filename for the third frame in the sequence
serv = getSQLServer()
db = 'disdc'
sel = ['Filepath']
cond = '(EXPID = '+strtrim(sxpar(fitsHdr, 'EXPID'),2)+$
	") AND (INSTRUME='HRIIR') AND (FNIMGNMB=3) AND (MSNCONFG='"+$
	sxpar(fitsHdr, 'MSNCONFG')+"') AND (OBSJD<="+$
	string(sxpar(fitsHdr,'OBSJD')+1,format='(D13.7)')+$
	") order by OBSJD DESC limit 1"
tblNm = 'HRIIR'
webFn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)

IF max(dim) EQ 0 THEN return, getModelIRDark(fitsHdr,img)

;; Get the third frame locally
fn = getLocFn(webFn, subdir='IRFRAME3')

;; Open up the third frame
dark = readfits(fn[0])
IF n_elements(dark) EQ 1 THEN return, getModelIRDark(fitsHdr, img)

;; Lineairze it
dark = dical_lineardn(dark, fitsHdr)

;; Update the header
fxaddpar, fitsHdr, 'DARKFN', $
	strupcase(strmid(fn[0], strpos(fn[0], path_sep(), /reverse_search)+1))

return, dark

end

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Retrieves a dark frame from the database for the given observation.
;;	If no frame is available, then -1 is returned. The dark must be
;;	linearized and must also include any bias.
;;
;; CALLING SEQUENCE:
;;	RESULT = getSpecificDark(fitsHdr)
;;
;; REQUIRED INPUTS:
;;	fitsHdr - header for the image
;;
;; OUTPUTS:
;;	RETURN - the dark frame to subtract from the image or -1 if no valid
;;		frame was found
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> dark=getSpecificDark(fitsHdr)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	DI_sqlQuery - Query the SQL server
;;	getSQLServer - Retrives info on how to acess the SQL server
;;	getLocFn - Converts a remote filename into a local one
;;
;; MODIFICATION HISTORY:
;;   2005-07-21  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------
FUNCTION getSpecificDark, fitsHdr

;; Create an exposure ID string with the part after the decimal being the image number
expid = double(strtrim(sxpar(fitsHdr, 'EXPID'),2)+'.'+string(sxpar(fitsHdr, 'FNIMGNMB'), format='(I3.3)'))

;; Get the filename for the specific frame
serv = getSQLServer()
db = 'di_calib'
sel = ['Filepath']
cond = '(StartExpID <= '+string(expid,format='(F11.3)')+$
	") AND (Instrument='"+strtrim(sxpar(fitsHdr,'INSTRUME'),2)+$
	"') AND (EndExpID >= "+string(expid,format='(F11.3)')+$
	") AND Date <= '"+strtrim(sxpar(fitsHdr,'OBSDATE'),2)+"' order by Date DESC, Version desc limit 1"
tblNm = 'Darks'
webFn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)

IF max(dim) EQ 0 THEN return, -1

;; Found a dark to use, so go get it
fn = getLocFn(webFn, SUBDIR='DARKS')

out = readfits(fn,darkhdr)

;; Update the fits header
fxaddpar, fitsHdr, 'DARKFN', $
	strupcase(strmid(fn, strpos(fn, path_sep(), /reverse_search)+1))
fxaddpar, fitsHdr, 'DARKALG', 'SPECIFIC'

;; See if it is an IR manually scaled dark and if so, fill in the appropriate areas of the header
garb = sxpar(darkHdr, 'MDARKVER',count=nmatch)
IF nmatch GT 0 THEN BEGIN
	;; Hack in the correct filename

	scalmdnm = strtrim(sxpar(darkHdr, 'SCALMDNM', count=c1),2)

	insertIRDarkModelParams, fitsHdr, darkHdr $
                               , (c1 gt 0) ? scalmdnm[0] : 'UNK'

	tmp = sxpar(darkHdr, 'SCALFAC')
	fxaddpar, fitsHdr, 'DARKMSCL', tmp

	fxaddpar, fitsHdr, 'DARKALG', 'MODEL*MANUAL'
ENDIF

return, out

END

;;-----------------------------------------------------------------------------
;; Main function
;;-----------------------------------------------------------------------------
FUNCTION DIcal_darkFrame, in, fitsHdr, flags, USEFIT=usefit, DARKFN=darkfn

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

IF error NE 0 THEN BEGIN
	CATCH, /CANCEL
	message, 'Subtract Dark Frame - ' + !ERROR_STATE.MSG, /NONAME
ENDIF

;; Get the instrument
inst = STRTRIM(SXPAR(fitsHdr, 'INSTRUME', COUNT=c1),2)

;; Get the detector ID
det = SXPAR(fitsHdr, 'IMGH021', COUNT=c2)

;; Get the date when the picture was taken
date = SXPAR(fitsHdr, 'OBSDATE', COUNT=c3)

;; Get the detector temperature
dTemp = SXPAR(fitsHdr, (det EQ 0)?'CCDT':'IRFPAT', COUNT=c4)

;; Get the bench tempearture
bTemp = sxpar(fitsHdr, 'SMOBENT', COUNT=c5)
if (bTemp EQ -1) THEN bTemp = sxpar(fitsHdr, 'OPTBENT', COUNT=c5)

;; Get integration time
intT = SXPAR(fitsHdr, 'INTTIME', COUNT=c6) * 1d-3     ;;; Seconds

c7=1

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

;;Make sure we were able to get all of the data from the header
IF (c1 EQ 0) OR (c2 EQ 0) OR (c3 EQ 0) OR (c4 EQ 0) OR $
	(c7 EQ 0) OR (c6 EQ 0) OR (c7 EQ 0) OR (c8 EQ 0) THEN $
	message, 'Invalid FITS file', /NONAME

out = double(in)

IF keyword_set(darkfn) EQ 0 THEN BEGIN
	dark = -1
	;; Try to get the specific dark frame for this observation
	IF NOT keyword_set(usefit) THEN BEGIN
		dark = getSpecificDark(fitsHdr)
	ENDIF

	;; Find calculate the interpolated dark frame
	IF n_elements(dark) EQ 1 THEN BEGIN  ;; Use model
		IF det EQ 1 THEN BEGIN ;; IR
			dark = getModelIRDark(fitsHdr,in)
			fxaddpar, fitsHdr, 'DARKALG', 'MODEL'
		ENDIF ELSE BEGIN
			dark = getModelVISDark(inst, date, dTemp, bTemp, intT, mode,fitsHdr)+$
				getBias(fitsHdr,inst,date,mode,flags,IMG=out,$
					OCLOCK = (mode NE 8) AND (mode NE 7))
			fxaddpar, fitsHdr, 'DARKALG', 'MODEL'
		ENDELSE
	ENDIF
ENDIF ELSE BEGIN
	dark = readFITS(darkfn, darkHdr)

	IF n_elements(dark) NE n_elements(in) THEN message, 'Dark frame is different size from the image'

	;; Add the filename to the FITS header
	fxaddpar, fitsHdr, 'DARKFN', $
		strupcase(strmid(darkfn, strpos(darkfn, path_sep(), /reverse_search)+1))
	fxaddpar, fitsHdr, 'DARKALG', 'USER FILE'
ENDELSE

;; Subtract the dark frame from the image but leave the bias rows
z = getBiasCols(fitsHdr, COMPLEMENT=imgCols)
out[imgCols[0]:imgCols[1],*] = $
	out[imgCols[0]:imgCols[1],*] - dark[imgCols[0]:imgCols[1],*]

;; Update the FITS header
fxaddpar, fitsHdr, 'DARKCORR', 'T'

RETURN, out

END