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