FORWARD_FUNCTION DICAL_vismodulemask, DICAL_vismodulestripes, dical_darkframe, 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.
;;;  DICAL_vismodulemask - Computes mask frames
;;;  DICAL_vismoduledark - Shift the Dark image to well remove the stripe feautres
;;;  DICAL_vismodulestripes - Removes stripe features
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-04-12  M. Desnoyer    Reordered functions
;;   2009-09-01  F. Merlin      Add DI_modulemask, DI_moduledark and DI_modulestripes procedures
;;   2010-08     B. Carcich     Modified for use in pipeline  
;;   2011-05     B. Carcich     Added call to DICALIR_VARINTTIME for line-dependent integration time offset
;;   2011-05     B. Carcich     Added return of column profile to call to dical_vismodulestripes
;;
;;-----------------------------------------------------------------------------

;;-----------------------------------------------------------------------------
;; 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
;;   2012-11-30  B. Carcich     Added keyword argument commentStruct
;;
;;-----------------------------------------------------------------------------
pro insertIRDarkModelParams, calHdr, darkHdr, mstFn, commentStruct=commentStructArg

;; Copy commentStruct keyword argument, if it exists and is a structure,
;;   to variable CS, else copy a dummy default structure (csDefault)

csDefault = { _:0b }   ;;; Only value is non-string, so nothing will be used

cs = n_elements(commentStructArg) eq 1L $
     and size(commentStructArg,/type) eq size(csDefault,/type) $
     ? commentStructArg[0] $
     : csDefault

;; CS tag names, if present, are keywords from darkHdr; corresponding
;;   string values are to be used as comment for converted header keywords

tns = tag_names( cs )

;; 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'] $
  , ['DRKMSCP',   'DRKMSCP'] $
  , ['DRKMSCPP',  'DRKMSCPP'] $
  , ['DRKIRPRO',  'DRKIRPRO'] $
  , ['TEMPFPA',   'TEMPFPA'] $
  , ['TEMPSIM',   'TEMPSIM'] $
  , ['TEMPSCAL',  'TEMPSCAL'] $
  , ['DRKTOOHI',  'DRKTOOHI'] $
  , ['DARK99  ',  'DARK99  '] $
  ]

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

    ;; Look for darkHdr keyword in tag names (TNS) of comment structure (CS)

    iw = where( tns eq hdrConvert[0,i], iwct)

    if iwct eq 1 and size(cs.(iw[0]>0),/type) eq size('',/type) then begin
      fxaddpar, calHdr, hdrConvert[1,i], tmp[0], cs.(iw[0])
    endif else begin
      fxaddpar, calHdr, hdrConvert[1,i], tmp[0]
    endelse
  endif
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

commentStruct = { _ : 0b }  ;;; dummy commentStruct which may be used later

;; 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')
fnPro=fn[0]

;; Copy the procedure into the working directory
di_file_copy, fnPro, "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', count=c)
if c eq 1 then fpaTemp = sxpar(fitsHdr, 'IRFPAT', count=c)
if c eq 1 then covTemp = sxpar(fitsHdr, 'COVERT', count=c)
if c eq 1 then priMirTemp = sxpar(fitsHdr, 'PRIMIRT', count=c)
if c eq 1 then secMirTemp = sxpar(fitsHdr, 'SECMIRT', count=c)

;; - For OPTical BENch Temperature, try to get a SMOothed BENch Temp. first:
simTemp = sxpar(fitsHdr, 'SMOBENT',count=csimTemp)
if simTemp LT 0 or csimTemp ne 1 THEN simTemp = sxpar(fitsHdr, 'OPTBENT', count=csimTemp)

;; Check the temperature inputs
IF c ne 1 or csimTemp ne 1 $
  or 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'

intTime = sxpar(fitsHdr, 'INTTIME', count=cInttime) $
          + dicalir_varinttime(mode=curMode)     ;;; Add any mode-dependent offset(s)
IF cInttime ne 1 then $
  Message, 'Invalid integration time in header: Dark cannot be calculated'

curSlot = sxpar(fitsHdr, 'FNIMGNMB', count=cSlot)
IF cSlot ne 1 then $
  Message, 'Invalid slot number 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 = (di_scu_missionstruct()).mrgdb
sel = ['OBSJD']
cond = '(EXPID = '+strtrim(sxpar(fitsHdr, 'EXPID'),2)+$
  ") AND (INSTRUME='HRIIR') AND (FNIMGNMB=1) AND (MSNCONFG='"+$
  sxpar(fitsHdr, 'MSNCONFG')+"') AND (OBSJD<="+$
  strtrim(string(sxpar(fitsHdr,'OBSJD'),format='(D20.8)'),2)+$
  ") order by OBSJD DESC limit 1"
tblNm = 'HRIIR'
;seqStart = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)
dim=[0,0]

IF max(dim) EQ 0 THEN BEGIN
  ;; Couldn'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 = (di_scu_missionstruct()).mrgdb
  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])

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; 2012-11-30 BTCarcich
;;; Look for dark scale factor in MANUAL_SCALE_FACTOR table.
;;; Initial values provided for Hartley 2 encounter by Silvia Protopapa.

serv = getSQLServer()
db = 'di_calib'
tick="'"
sel = ['scalfact','provenance']

cond = " scstart < " + tick + strtrim(sxpar(fitsHdr,'SCSTOP'))  + tick $
     + " AND scstop  > " + tick + strtrim(sxpar(fitsHdr,'SCSTART')) + tick $
     + " AND expid=" + tick + strtrim(sxpar(fitsHdr,'EXPID')) + tick $
     + " ORDER BY scstart DESC LIMIT 1"

tblNm = "MANUAL_SCALE_FACTOR"
scalfactRtn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)

if max(dim) gt 0 then begin
  ;;; if SCALFACT was found, update dark and darkHdr as follows:
  ;;; - re-scale dark (variable out) by SCALFACT / TIMESCAL
  ;;;   - Use 1 for TIMESCAL if it is not present or it is LE 0d0
  ;;; - set TIMESCAL to -999
  ;;; - Provide TIMESCAL comment indicating overridden value
  ;;; - Update DRKMSCP and DRKMSCPP keywords in dark header, which
  ;;;   will be added to fitsHdr by procedure insertirdarkmodelparams

  scalfact = double(scalfactRtn[0])
  provenance = strtrim(scalfactRtn[1],2)
  oldTIMESCAL = double(strtrim(sxpar(darkHdr,'TIMESCAL'),2))
  out *= oldTIMESCAL le 0d0 ? scalfact : (scalfact / oldTIMESCAL)
  sxaddpar,darkHdr,'TIMESCAL',-999d0

  commentStruct = create_struct( commentStruct $
    , 'TIMESCAL' $
    ,  ' Manual dark scaling(' $
       + strtrim(string( f='(F20.5)', oldTIMESCAL),2) $
       + ') vetoed by DRKMSCP')

  sxaddpar,darkHdr,'DRKMSCP',scalfact
  sxaddpar,darkHdr,'DRKMSCPP',provenance
endif

iSlash=strpos(fnPro,'/',/reverse_search)+1L
darkHdr = [ darkHdr $
          , "DRKIRPRO= '" $
            +strmid(fnPro,iSlash) $
            +"' / DARK_IR IDL PRO used" $
          , '' $
          ]
;; Update the FITS header
insertIRDarkModelParams, fitsHdr, darkHdr, fn[0], commentStruct=commentStruct

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 = (di_scu_missionstruct()).mrgdb
sel = ['Filepath']
cond = '(EXPID = '+strtrim(sxpar(fitsHdr, 'EXPID'),2)+$
  ") AND (INSTRUME='HRIIR') AND (FNIMGNMB=3) AND (MSNCONFG='"+$
  sxpar(fitsHdr, 'MSNCONFG')+"') AND (OBSJD<="+$
  strtrim(string(sxpar(fitsHdr,'OBSJD')+1,format='(D20.8)'),2)+$
  ") 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'
tblNm = 'Darks'
sel = ['Filepath','Filename']

;;; - Search for calibration files which match the Exposure ID,
;;;   instrument, and have either the same date as the day of,
;;;   or the day before, the observation

obsdate = strtrim( sxpar(fitsHdr,'OBSDATE'), 2 )
ymd = long( strsplit( strmid(obsdate,0,10), '-', /extract) )       ;;; Year, month, day of obs
caldat, julday( ymd[1], ymd[2], ymd[0], 12 ) - 1d0, mm1, dm1, ym1  ;;; Y,M,D of day before
obsdatem1 = string( f='(i4.4,2(1h-,i2.2))', ym1, mm1, dm1 )        ;;; 'YYYY-MM-DD' of day before

cond = 'StartExpID <= '+string(expid,format='(F11.3)')+$
  " AND Instrument='"+strtrim(sxpar(fitsHdr,'INSTRUME'),2)+$
  "' AND EndExpID >= "+string(expid,format='(F11.3)')+$
  " AND Date <= '"+obsdate+"'"+ $
  " AND Date >= '"+obsdatem1+"' ORDER BY Date DESC, Version DESC"

;;; Make the MySQL SELECT call

webFn = DI_sqlQuery(serv, db, tblNm, sel, cond, DIM=dim)

IF max(dim) EQ 0 THEN return, -1   ;;; No specific calibration file found

if n_elements(dim) eq 1L then dim = [dim,1L]

;;; Compare RECSCLKD (encoded SCLK) from fitsHdr to STSCLKDP and
;;; EDSCLKDP in calibration files

;;; Ensure fn and noStartEndFn do not exist
;;; - fn will be the calibration file to use if the RECSCLKD from
;;;     the data file is between STSCLKDP and EDSCLKDP of the
;;;     calibration file
;;; - if a file is the first in which STSCLKDP is not present in a
;;;     calibration file, put it into noStartStopFn

ptr_free,ptr_new(fn,/no_copy)
ptr_free,ptr_new(noStartEndFn,/no_copy)

recsclkd = sxpar(fitsHdr,'RECSCLKD', count=ct)

if ct eq 1L and dim[n_elements(dim)-1] gt 0L then begin

  ;;; Loop through calibration file Filepaths that match the SELECT call
  newV = ''

  for i=0L,dim[1]-1L do begin

    ;;; Get filename (HRIIR_YYMMDD_Version_Mode_ExpID*), parse to
    ;;; eliminate version (third underscore-separated token) as a unique
    ;;; part of filename, don't use calibration files that have the same
    ;;; basic filename but differ in version

    lastV = newV
    toks = strsplit(webFn[1,i],'_',/extract)
    toks[2] = 'V'                                ;;; Replace the version
    newV = strjoin( toks, '_')
    if newV eq lastV then continue

    ;;; Get Filepath and FITS header from calibration file [i]

    nextFn = getLocFn(webFn[0,i], SUBDIR='DARKS')
    calfileHdr = headfits( nextFn )

    ;;; Look for ST/EDSCLKDP in FITS header, go to next calibration file
    ;;; if either is not found

    stsclkdp = sxpar( calfileHdr, 'STSCLKDP', count=ct )
    if ct eq 1L then begin
      edsclkdp = sxpar( calfileHdr, 'EDSCLKDP', count=ct )
    endif else begin
      if n_elements(noStartEndFn) eq 0L then noStartEndFn = nextFn
    endelse
    if ct ne 1L then continue

    ;;; Compare RECSCLKD (encoded SCLK) from fitsHdr to STSCLKDP and
    ;;; EDSCLKDP in calibration file

    if stsclkdp le recsclkd and recsclkd le edsclkdp then begin

      ;;; If matched, this is the file; the ORDER BY above

      fn = nextFn
      break

    endif

  endfor
endif

;;; If no cal file satisfied the STSCLKDP/RECSCLKD/EDSCLKDP test, then
;;; check if there was a cal file that did not have STSCLKDP in its FITS
;;; header

if n_elements(fn) eq 0 then begin
  if n_elements(noStartEndFn) eq 0L then return, -1L
  fn = noStartEndFn
endif

;; Found a dark to use, so go get it

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 $
       , rtnPY=rtnPY

;; 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)
      darkalg = fxpar( fitsHdr,'DRKMSCP') le 0d0 ? 'MODEL' : 'MODEL*MANUAL(PROTOPAPA)'
      fxaddpar, fitsHdr, 'DARKALG', darkalg
    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

z = getBiasCols(fitsHdr, COMPLEMENT=imgCols)

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;;; Test if the stripes removal algorithm is to be attempted i.e. if 
  ;;; the signal in 75% or less of the the FOV is high (>500 DN above dark) 
  ;;; i.e. > 0.5% error

  out2=out-dark
  if det eq 1 then begin    ;;; IR
    nb_est = 2.0
  endif else begin
    out_est=where(out2[imgCols[0]:imgCols[1],imgCols[0]:imgCols[1]] gt 500, tot_est) 
    nb_est=tot_est/((imgCols[1]-imgCols[0]+1)^2.0)
  endelse

  didRmStripe = 'F'

  IF nb_est LE 0.75 and det EQ 0 THEN begin

    ;;; Compute the mask

    inimagem=out
    darkimagem=dark
    mask=DICAL_vismodulemask(inimagem, darkimagem, imgCols)

    out1=out-dark                           ;;; Subtract the dark frame

    ;;; Subtract the residual stripe features when it is possible but 
    ;;; leave the bias rows

    INSTRUME = SXPAR(fitsHdr, 'INSTRUME', count=res)
    if res eq 1 then begin
      TESTNAVi = strmid(strtrim(INSTRUME,2),0,3) eq 'NAV'
    endif else begin
      INSTN = SXPAR(fitsHdr, 'INSTRNAV', count=res)
      TESTNAVi = res eq 1
    endelse
    
    out2 = DICAL_vismodulestripes(out1, mask, imgCols, isNav=TESTNAVi, rtnProfileY=rtnPY) 

    ;;; If DICAL_vismodulestripes succeeds, out2 contains the stripe-removed out
    ;;; If DICAL_vismodulestripes fails, out2 = FIX([0])

    IF N_ELEMENTS(out2) EQ N_ELEMENTS(out) THEN BEGIN
      didRmStripe = 'T'                           ;;; Success, set keyword
    ENDIF else begin
      out2=out-dark                               ;;; Failure, subtract dark
    ENDELSE

    IF TESTNAVi THEN BEGIN
      poszero=where(out EQ 0,poszeroCount)
      if poszeroCount gt 0L THEN out2[poszero]=0
    ENDIF

  ENDIF

  fxaddpar, fitsHdr, 'RMSTRIPE', didRmStripe, ' Stripe removal algorithm applied (T/F)'

  out[imgCols[0]:imgCols[1],*]=out2[imgCols[0]:imgCols[1],*]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

RETURN, out

END