function dark_sub, index, data, index_out, $
                udata_in=udata_in, udata_out=udata_out, $
		update_index=update_index, $
		dc_only=dc_only, $
		force_scalar=force_scalar, interpolate=interpolate, $
		dc_scalar=dc_scalar, dc_interpolate=dc_interpolate, $
		float=float, norm_exp=norm_exp, shift_floor=shift_floor, new_floor=new_floor, $
		qstop=qstop, qdebug=qdebug, match_qual=match_qual,$
		orbit_correct=orbit_correct, no_orbit_correct=no_orbit_correct
;
;+
;NAME:
;	dark_sub
;PURPOSE:
;	Subtract the dark current from input FFI or PFI images
;CALLING SEQUENCE:
;	data_out = dark_sub(index, data, index_out)
;	data_out = dark_sub(index, data, /update_index)
; 
;       udata_out = 1		;for uncertainties
;	data_out = dark_sub(index, data, index_out, udata_out=udata_out)
;INPUT:
;	index	- The index structure for each image
;	data	- The data array (NX x NY x N)
;OUTPUT:
;	data_out- The dark current subtracted image
;	index_out- The updated index struture (history records appended)
;OPTIONAL KEYWORD INPUT:
;       udata_in - the uncertainty in the input data array. You pass this in if
;                  the data has been uncompressed already, if the data is of byte type,
;                  dark_sub will obtain the uncertainty
;       udata_out - The uncertainty in the output array, this will only be calculated
;                   if a variable is set and passed in as the udata_out keyword.
;	dc_only	- If set, the output is only the dark current
;		  images that go with each index.
;	dc_scalar - If set, then do not call GET_DC_IMAGE, just use
;		  a scalar value for the background level.
;	force_scalar - Same as "dc_scalar"
;	dc_interpolate - If set, then find a short and long exposure dark image
;		  of the same temperature around the time of the data, and
;		  interpolate between those exposures to create the dark
;		  frame.
;	interpolate - Same as "dc_interpolate"
;	float	- If set, return the data as floating point if necessary
;		  (for the case where PFI startline correction needs to be made)
;	norm_exp- If set, exposure normalize the data to DN/sec and return it as
;		  floating point.
;	shift_floor- If set, check the average of a 5x5 box in the lower 
;		  left corner and shift the image by that amount in order to
;		  adjust the image floor to the true zero.  This was 
;		  introduced because for short exposure (DPE=0) a small
;		  error in the "floor" is amplified quite a bit.
;	new_floor - If shift_floor is set, it will adjust the floor to be
;		  0 DN.  It can be shifted to a non-zero value by passing
;		  that value in through the "new_floor" keyword.
;	update_index - If set, then update the history portion of the
;		  input index (not pass separate index variable out)
;       orbit_correct - Use the orbital corrections for the dark current.
;                       Must also specify the /interpolate keyword for this to
;                       work.
;OPTIONAL KEYWORD OUTPUT:
;	match_qual - A structure which holds information on the
;		  quality of the match for the dark image (number of
;		  seconds between dark image and input image, how close
;		  the exposure duration is, how close the temperature is)
;METHOD:
;	1. Call SXT_DECOMP if "data" is byte type
;	2. Call GET_DC_IMAGE to get dark images
;	3. For PFI, extract the relevant portion of the images
;	4. Subtract off the dark current
;HISTORY:
;	Written 13-Oct-92 by M.Morrison
;	15-Oct-92 (MDM) - Old method was to call GET_DC_IMAGE one time
;			  outside the loop and to use IMAP.  For FR images
;			  this can be a lot of images/memory, so it was 
;			  changed to call GET_DC_IMAGE inside the loop and
;			  to sort the order that images are subtracted so
;			  that they are done in increasing DPE order.
;	19-Nov-92 (GAL) - update to allow input of a single image.
;	12-Mar-93 (MDM) - Modified to work in the cases when the
;			  SDC database is not available
;	18-Mar-93 (MDM) - Added /FORCE_SCALAR keyword
;			- Added code to handle the case where the full
;			  resolution PFI/OR falls outside of the dark current
;			  image (since it is only 1024x512).
;	12-May-93 (MDM) - Patch for case where whole full res PFI is outside
;			  of the FFI dark frame
;			- Modification to insert SUB_IMG into DATA_OUT using
;			  (0,0,i) notation instead of (*,*,i) because of case
;			  where there is something like a 64x64 imbedded in
;			  a 128x128 data array.
;	13-May-93 (MDM) - Added /FLOAT and /NORM_EXP options
;	14-May-93 (MDM) - Correct DN_SHIFT calculation for FR FFI dark images
;			  because they don't start at line 0.
;			- Added /SHIFT_FLOOR option
;V2.0	 6-Jun-93 (MDM) - Started tracking the version number
;			  Added call to HIS_INDEX to track history.
;	 8-Jun-93 (MDM) - Handle FR case where whole PFI falls outside
;			  of the FFI dark frame
;	10-Jun-93 (MDM) - Made exposure normalization in seconds 
;			  instead of in msec.
;	23-Jun-93 (MDM) - Put on-line
;			- Added call to HIS_INDEX for when exposure normalized
;			- Modified to use SHAPE_SAV instead of SHAPE_CMD
;V2.01	 3-Aug-93 (MDM) - Modified to work with full frame strips
;	16-Aug-93 (MDM) - Modified so that when DC_INDEX comes in defined, it
;			  does not cause problems (since GET_DC_IMAGE will not
;			  re-read the image if DC_INDEX comes in defined and is
;			  the same image which is requested)
;V2.02	16-Aug-93 (MDM) - Corrected an error caused when the full PFI or OR
;			  fell outside of the reference dark frame.
;			- Modified GET_DC_IMAGE to try to get the right 
;			  FFI strip when in FR (top half of bottom half)
;V2.03	16-Aug-93 (MDM) - Corrected bug introduced in 2.02
;V3.00	25-Aug-93 (MDM) - Major changes to GET_DC_IMAGE and DARK_SUB
;			- Added INDEX_OUT and /UPDATE_INDEX options
;			- Added interpolation option
;			- Allowed this one routine to be used for warm dark
;			  image subtraction also
;	26-Aug-93 (MDM) - Added $ys_dark_interp to switch interpolation on
;			- Removed DC_INDEX and DC_IMAP outputs (since the
;			  information is now saved in the history index)
;	 6-Oct-93 (SLF) - Add output keyword MATCH_QUAL (from get_dc_image)
;	13-Oct-93 (MDM) - Minor bug fix in the /DC_SCALAR option
;V3.01  16-Dec-93 (MDM) - Modified SEL_DC_IMAGE to not chose dark images on
;                         the "other side" of 13-Nov-92
;V3.02	14-Jan-94 (MDM) - Added /DC_INTERPOLATE and /DC_SCALAR
;			- Corrected error with history records when no
;			  SDC datasets are found
;V3.03	26-Jan-94 (MDM) - Corrected program version number variable (was
;			  3.01 for 3.02 modification
;	 9-Jun-94 (MDM) - Modified header (added information)
;V3.04	30-Nov-94 (MDM) - Modified to use /ORIGINAL switch on GT_CORNER
;			  call (because of how the DN shift is applied
;			  and since RD_XDA will now allow the reading of
;			  a portion of the image).
;V3.05	12-Dec-94 (MDM) - Modified to pass FLOAT keyword through to
;			  GET_DC_IMAGE
;       08-Feb-95 (jmm) - Added the UDATA keywords
;V3.10	27-Feb-95 (MDM) - Put McTiernan changes on-line
;			- Check that the input is 2-D or 3-D
;			- Put check that ROI location adjustment is done
;			  properly (it was applying a negative adjustment)
;V3.15  03-Feb-98 (TRM) - Added orbital correction code
;       10-Feb-98 (SLF) - made ORBIT_CORRECT the default, added keyword
;                            /NO_ORBIT_CORRECT to override  
;       11-Feb-98 (SLF) - made NO_ORBIT_CORRECT the default 
;V3.20  23-Jul-98 (SLF) - implement bug fix for ROI/PFI dn_shift correction
;                         as suggested by shimojo san.
;V3.25  02-Apr-99 (TRM) - changed val to a string so that extra values can
;                         be included in the uniqueness test.  The new values
;                         are the week id, to make sure that a new dark is
;                         selected if a week boundary is passed, and tfms
;                         which is used to force a new dark frame when 
;                         orbit correction is used.
;V3.26  14-Jul-99 (BNH) - Unbroke for /DC_SCALAR case.
;V3.27  22-Sep-99 (SLF) - Use bit 15 in Version for ORBIT CORRECT history
;V3.28  29-sep-00 (SLF) - Add calls to 'fstring' (allow more than 
;                         1024 input 'index,data'
;	29-Dec-01 (LWA) - Commented out '%DARK_SUB: applying SDC orbital correction'
;-
;orbit_correct=1-keyword_set(no_orbit_correct)
orbit_correct=keyword_set(orbit_correct)
progverno = round(3.28*1000)                         
progverno = progverno + ([0,2])(orbit_correct)^15    ; set bit for -> .HISTORY

n = n_elements(index)
qerr = 0
;
quncert = keyword_set(udata_out) 		;jmm, 8-feb-95

if (keyword_set(dc_only)) then begin
    nx = max(index.sxt.shape_cmd(0))
    ny = max(index.sxt.shape_cmd(1))
    data_out = intarr(nx, ny, n)
    if (quncert) then udata_out = data_out
end else begin
    siz = size(data)
    nx = siz(1)
    ny = siz(2)
    if (siz(0) eq 3) then nz=siz(3) else nz=1	;MDM added 27-Feb-95
    if (siz(0) lt 2) then begin
	tbeep, 3
	print, 'DARK_SUB: Error - input must be 2-D or 3-D.  Returning...'
	return, 0b
    end
    if (n_elements(index) ne nz) then begin
	tbeep, 3
	print, 'DARK_SUB: INDEX size does not match data size.  Returning...'
	return, 0b
    end
    typ = siz( siz(0)+1 )
; New stuff, jmm, 8-feb-95
    if (typ eq 1) then begin
       if(quncert) then data_out = sxt_decomp(data, udata_out) $
         else data_out = sxt_decomp(data)
    end else begin
       data_out = data
       if (quncert) then begin
          if (n_elements(udata_in) ne 0) then begin
             if (n_elements(udata_in) ne siz(siz(0)+2)) then begin
                message, 'UDATA_IN AND DATA HAVE DIFFERENT SIZES, USING SCALAR FOR UDATA', /cont
                udata_out = make_array(size = size(data_out)) + udata_in(0)
             end else udata_out = udata_in
          end else udata_out = make_array(size = size(data_out)) + 0
       end
    end
end
;
qfloat = 0					;added 13-May-93
if (keyword_set(float)) then qfloat = 1
if (keyword_set(norm_exp)) then qfloat = 1
if (keyword_set(shift_floor)) then qfloat = 1
if (qfloat) then data_out = float(temporary(data_out))
if (qfloat) and (quncert) then udata_out = float(temporary(udata_out))	 ;jmm, 1/19/95, MDM 27-Feb-95 added "quncert"
qinterp = keyword_set(interpolate) or keyword_set(dc_interpolate)
if (getenv('ys_dark_interp') ne '') then qinterp = 1
;
;---------------------------------------- Deterimine the optimal order of the images to do the subtraction
xcorn = fix(gt_corner(index, /lower, /x)/4)     ;0-255
ycorn = fix(gt_corner(index, /lower, /y)/4)
res = gt_res(index)
dpe = gt_dpe(index)
dpe = dpe > 14	;anything under 0.118 sec is treated the same - otherwise too slow and not significant differences
week = strarr(n_elements(index))
for iweek=0,n_elements(week)-1 do week(iweek) = tim2weekid(index(iweek))

;;do all FR in order of DPE, then all HR, then all QR
;val =      xcorn         +  ycorn*256L               + dpe*256L*256L  + res*256L*256L*64
;;      low 8-bits 0:7      second 8-bits  8:15          third 8-bits
;;      fourth 8-bits

; TRM converted val to a string and added week,tfms 02-Apr-1999
val = fstring(xcorn,format='(i5)') + $
      fstring(ycorn,format='(i5)') + $
      fstring(dpe,format='(i2)')   + $
      fstring(res,format='(i1)')   + $
      week
if keyword_set(orbit_correct) then begin
   tfms = sxt_uvf_info(index,/tfms)
   ; 0.01 minute granularity essentially
   ; guarantees that all val's will be unique
   val = val + fstring(fix(tfms*100.+0.5),format='(i5)')
endif
ispaces = where(byte(val) EQ 32b,nspaces)  ; Convert leading spaces to zeroes so that the sort will work
if nspaces GT 0 then begin
   vtemp = byte(val)
   vtemp(ispaces) = 48b
   val = string(vtemp)
endif

ss = sort([val])
;
qprint_message2 = 1
last_val = -999
qno_sdc = 0
if (keyword_set(force_scalar)) then qno_sdc = 1		;MDM added 18-Mar-93
if (keyword_set(dc_scalar)) then qno_sdc = 1		;MDM added 14-Jan-94
if (keyword_set(qstop)) then stop
;
dn_shift = 0
index_out = index
for iimg=0,n-1 do begin
    i = ss(iimg)		;do the subtraction in order of DPE so that
    ;
    ;
    ; TRM 1998-Feb-3
    ; Force get_dc_image for every data image so that the orbital corrections
    ; are applied properly for each data image.
    ;
    if ((not qno_sdc) and ((val(i) ne last_val))) then begin ;jmm, 8-feb-95
       udc_data = quncert       ;you need to pass in a value to get the DC uncertainty
;       if orbit_correct  then message,/info,"applying SDC orbital correction"
       get_dc_image, index(i), dc_index0, dc_data, interpolate = qinterp, match_qual = match_qual, $
         qerr = qerr, float = float, unc_data = udc_data,orbit_correct=orbit_correct
    end

    if (not qno_sdc) then $     ; (BNH) to fix /DC_SCALAR
        his_index, dc_index0	;MDM added 14-Jan-94 (since there is no history record for no-SDC file case)
    ;Do not re-read the data if the selected data is the same as the last data read
    ;
    if ((n_elements(dc_data) lt 10) or (qerr)) then begin
	;cannot find SDC dataset - use a scalar value for subtraction
        dc_data = fix(get_res_offset(index(i))) ;returns a scalar
        if(quncert) then udc_data = 1            ;jmm, 8-feb-95, i guess
	if (not qno_sdc) then begin
	    print, 'DARK_SUB: Warning - could not find SDC files'
	    print, 'DARK_SUB: Warning - using a value of ', strtrim(dc_data,2), ' for background for an image'
	end
	qno_sdc = 1
    end else begin
	qffi = gt_pfi_ffi(index(i), /true_ffi)
	dn_shift = 0
	if (qffi eq 0) then begin				;PFI/OR or FFI strip
	    npix = 2^gt_res(index(i))                         ;number of pixels binned (1,2 or 4)
	    c_cmd_y = gt_corner(index(i), /y, /original)	;MDM added /ORIGINAL 30-Nov-94
	    nguard = mask(index(i).sxt.flush, 4, 4)*8		;MDM corrected 10-Jun-93

	    c_cmd_y_dc = dc_index0(0).sxt.corner_cmd(1)*4	;*4 because reformatter error in FFI data
	    nguard_dc = mask(dc_index0(0).sxt.flush, 4, 4)*8
	    c_cmd_y_dc0 = (c_cmd_y_dc-nguard_dc) > 0

	    nsec_less = (c_cmd_y-c_cmd_y_dc0)*.008/npix - nguard*.008	;MDM 10-Jun-93
	    nsec_less = nsec_less > 0					;MDM 27-Feb-95
	    dn_shift = nsec_less * get_dc_rate(index(i))* (npix^2)
	    if (qprint_message2) then print, dn_shift, ' DN shift because of ROI location'
	    if (not qfloat) then dn_shift = fix(dn_shift + 0.5)		;round off
	    qprint_message2 = 0
	end
    end

    if (keyword_set(dc_only)) then begin
	if (n_elements(dc_data) eq 1) then data_out(*,*,i) = dc_data $
                                      else data_out(0,0,i) = dc_data - dn_shift ;insert into data cube (doesn't work for scalar)
	if(quncert) then if (n_elements(unc_data) eq 1) then udata_out(*, *, i) = udc_data $
          else udata_out(0, 0, i) = udc_data ;jmm, 8-feb-95
        
    end else begin
	;data_out(*,*,i) = data_out(*,*,i) - sub_img
        ;data_out(0,0,i) = data_out(*,*,i) - sub_img	;does not work - does every other line for 64x64 imbedded in 128x128

	siz = size(dc_data)
	nx2 = siz(1)-1
	ny2 = siz(2)-1
	;-- new if-then statements, jmm 8-feb-95        
	if (qno_sdc) then begin 
           data_out(0, 0, i) = data_out(*, *, i) - dc_data 	;subtract scalar value
           if (quncert) then if (qfloat) then udata_out(0, 0, i) = sqrt(udata_out(*, *, i)^2.+udc_data^2.) $;jmm, 8-feb-95
             else udata_out(0, 0, i) = fix(0.5+sqrt(udata_out(*, *, i)^2.+udc_data^2.)) ;round, don't truncate
        end else begin
           data_out(0, 0, i) = data_out(0:nx2, 0:ny2, i) - (dc_data - dn_shift)
           if (quncert) then if (qfloat) then udata_out(0, 0, i) = sqrt(udata_out(0:nx2, 0:ny2, i)^2.+ udc_data^2.) $
             else udata_out(0, 0, i) = fix(0.5+sqrt(udata_out(0:nx2, 0:ny2, i)^2.+udc_data^2.))
        end


	his_index, index_out, i, 'q_dark_sub', progverno
	if (qno_sdc) then begin
	    his_index, index_out, i, 'time_dark', dc_data			;put the scalar value into the time field
	end else begin
	    his_index, index_out, i, 'time_dark', dc_index0.his.time_dark
	    his_index, index_out, i, 'day_dark',  dc_index0.his.day_dark
	end
    end
    ;thought about not making "sub_img = dc_data(*,*,imap(i))" for FFI case
    ;because did not want to repeatedly copy the FFI image to "sub_img", but
    ;the command "dc_data(*,*,imap(i))" itself makes a separate copy in memory (doesn't it?)

    if (keyword_set(shift_floor)) then begin
	if (iimg eq 0) then shift_arr = fltarr(n)
	temp = data_out(0:4,0:4,i)
	signal = total(temp)/n_elements(temp)
	if (n_elements(new_floor) eq 0) then new_floor = 0.0
	shift_arr(iimg) = signal
	data_out(*,*,i) = data_out(*,*,i) - signal + new_floor
	if (iimg eq n-1) then print, 'DARK_SUB: Floor shift min,max,avg=', min(shift_arr), max(shift_arr), total(shift_arr)/n
    end

    if (keyword_set(norm_exp)) then begin
	data_out(*,*,i) = data_out(*,*,i) / gt_expdur(index(i)) * 1000	;added 13-May-93
	if(quncert) then udata_out(*,*,i) = udata_out(*,*,i) / gt_expdur(index(i)) * 1000;jmm added 8-Feb-95
	his_index, index_out, i, 'expdur', 1000.   ;normalized to 1 sec = 1000 msec
    end

    last_val = val(i)
end
;
if (keyword_set(qstop)) then stop
if (keyword_set(update_index)) then index = index_out
return, data_out
end