; MODIFICATION HISTORY:
; Version Date: 
;   v2. 2009 Nov 10  AJS  Update to call ra_darkrate.pro

PRO ra_darksubtract, calibration_dir, primary_header, image, error_image, wavelength_indices

CSPICE_UTC2ET, sxpar(primary_header, 'STRTSCET'), et_start

exptime = float(sxpar(primary_header, 'EXPTIME'))
hvsetc = sxpar(primary_header, 'HVSETC')

CASE 1 OF 
   hvsetc LE -3900 AND hvsetc GT -4000: darkfile = calibration_dir + 'dark/RA_DARK_3.9KV_V003.FIT'
   hvsetc LE -3800 AND hvsetc GT -3900: darkfile = calibration_dir + 'dark/RA_DARK_3.8KV_V003.FIT'
   hvsetc LE -3700 AND hvsetc GT -3800: darkfile = calibration_dir + 'dark/RA_DARK_3.7KV_V003.FIT'
   ELSE: GOTO, skip
ENDCASE

rdfits_struct, darkfile, temp, /silent
dark = temp.im0
darkerr = temp.im1

; The region in which the dark is measured is [50:950, *]
darktotal = total(dark[50:950, *])
restofdark = total(dark[0:49,1:31]) + total(dark[951:*,*])

darkrate = ra_darkrate(et_start)

darknorm = darkrate / darktotal * exptime

IF n_params() GT 3 THEN begin
; Interpolate dark to data image
   yindex = (fltarr(1024) +1.) # findgen(32)
   dark = interpolate(dark, wavelength_indices, yindex, cubic = -0.5)
   darkerr = interpolate(darkerr, wavelength_indices, yindex, cubic = -0.5)

   image = adderr(image, error_image, dark * darknorm, darkerr * darknorm, $
                  err = error_image, /subtract)
ENDIF ELSE image -= (darkrate * exptime) / (901. * 32.)

sxaddpar, primary_header, 'DARKFILE', /savecomment, file_basename(darkfile)
sxaddpar, primary_header, 'DARKNORM', darknorm/exptime, /savecomment
sxaddpar, primary_header, 'DARKFLAG', 'T', /savecomment
skip:foo = 0. ; No dark subtraction

END