Viewing contents of file '../idllib/ghrs/pro/calhrs_bck.pro'
;***************************************************************************
;+
;*NAME:
; calhrs_bck
;
;*PURPOSE:
; Routine to compute and subtract the background from GHRS spectra.
;
;*CALLING SEQUENCE:
; calhrs_bck,ih,data,eps,fwid,bcktab,bck_scale,first,nf,nmerge,ihf,flux,
; ihb,back,epsb,log
;
;*PARAMETERS:
; INPUTS:
; ih - header array for unmerged data bins (128 x N)
; data - data array for unmerged data bins (512 x N)
; eps - data quality array for unmerged data bins (512 x N)
; fwid - integer array of 4 elements
; fwid(0) - interorder median filter width
; fwid(1) - interorder mean filter width
; fwid(2) - sky median filter width
; fwid(3) - sky mean filter width
; fwid(4) - order of polynomial for smoothing the background
; fwid(5) - order of polynomial for smoothing the sky
; bcktab - table of background scaling fudge factors
; bck_scale - vector of fudge factors to use instead of table values
; bck_scale(0) = back_a
; bck_scale(1) = back_b
; bck_scale(2) = back_c
; bck_scale(3) = back_d
; first - first bin in unmerged data for each readout.
; integer vector of length M
; nf - number of bins found for each readout
; integer vector of length M
; nmerge - number of bins merged for gross spectra
; (0 means no merging)
;
; INPUT/OUTPUTS:
; ihf - headers for the gross spectra (128 x M)
; flux - flux vectors for each readout for which the background
; is to be subtracted.
;
; OUTPUTS:
; ihb - headers for the background spectra (128 x M)
; back - background spectra (500 x M)
; epsb - epsilon for the background spectra (500 x M)
;
; OPTIONAL INPUT/OUTPUTS:
; log - processing log (string array)
;
; METHOD:
; The background is determined using one of four possiblities.
; type 0 - sky
; type 1 - interorder observed width main diode array
; type 2 - background diodes from background data bins
; type 3 - background diodes from gross spectral bins
; If nmerge is 0 (no spectral merging) then type 3 background
; is used. If substep bin ids 5 or 6 are present then type 0
; is used. If substep bin ids 3 or 4 are present the type 1
; is used. if substep bin ids greater than 6 are present then
; type 0 is used. If none of the previous conditions are
; satisfied, type 3 is used. If multiple types can be used, the
; the one with the smaller type number is used.
; HISTORY:
; version 1 D. Lindler March 1989
; version 2.0 D. Lindler March 1990 - added polynomial fit
; to background and the background scale factors
; (ie the Cardelli fudge factors)
; version 2.1 D. Lindler Feb, 1993 added default scattered light
; table to CALHRS.
;-
;-------------------------------------------------------------------------
pro calhrs_bck,ih,data,eps,fwid,bcktab,bck_scale,first,nf,nmerge,ihf,$
flux,ihb,back,epsb,log
VERSION = 2.1
FILL = 800 ;epsilon value for fill data
NORM5 = 64. ;binid=5 normalization of sky from ssa to lsa
NORM6 = 1/64. ;binid=6 normalization of sky form lsa to ssa
;
; Determine type of background subtraction. Lower type takes precidence
;
binids_all=ih(53,*)
if nmerge eq 0 then begin
type = 3
end else begin
binids_present=intarr(16) ;up to sixteen possible ids
binids_present(binids_all)=1 ;which types are present
if total(binids_present(5:6)) gt 0 then type = 0 $
else if total(binids_present(3:4)) gt 0 then type = 1 $
else if total(binids_present(8:15)) gt 0 then type = 2 $
else type = 3
end
;
; filter widths --------------------------------------------------------------
;
if type eq 0 then begin
median_width=fwid(2)
mean_width=fwid(3)
poly_order = fwid(5)
endif
if type eq 1 then begin
median_width=fwid(0)
mean_width=fwid(1)
poly_order = fwid(4)
endif
;
; update log
;
case type of
0: stype='Sky spectrum'
1: stype='Interorder (main diode) array'
2: stype='Background diodes from background bins'
3: stype='Background diodes from gross spectral bins'
endcase
if type lt 2 then hist=strarr(5) else hist=strarr(2)
hist(0)='CALHRS_BCK version '+string(version,'(f5.2)')+ $
': Background subtraction'
hist(1)=' Type: '+stype
if type lt 2 then begin
hist(2)=' Background smoothing: Median ' + $
'filter width ='+string(median_width,'(i3)')
hist(3)=' Mean '+ $
'filter width ='+string(mean_width,'(i3)')
hist(4)=' Background fit by polynom' + $
'ial of order ='+string(poly_order,'(i3)')
endif
if n_elements(log) gt 0 then sxaddhist,hist,log
if !dump gt 0 then printf,!textunit,hist
;
; read background scale factors ------------------------------------------------
;
if type eq 1 then begin
if strupcase(strtrim(bcktab)) ne 'NONE' then begin
hist = ' Background scale factor table: '+strtrim(bcktab)
sxaddhist,hist,log
if !dump gt 0 then printf,!textunit,hist
;
; get from table (if not found in table use defaults a=1,b=1,c=0,d=0)
;
if ihf(53) eq 1 then aperture = 'SSA' else aperture = 'LSA'
bck_a = 1.0 & bck_b = 1.0 & bck_c = 0.0 & bck_d = 0.0
order = ihf(50)>1
grating = [' ','G-1','G-2','G-3','G-4','G-5','E-A','E-B']
grating = grating(ih(48))
table_ext,bcktab,'GRATING,ORDER,BACK_A,BACK_B,BACK_C,BACK_D'+ $
',APERTURE',gratings,orders,back_a,back_b, $
back_c,back_d,aper
n = n_elements(gratings)
good = where(orders eq order) & ngood = !err
if ngood gt 0 then begin
for i=0,ngood-1 do begin
ipos = good(i)
if (grating eq strtrim(gratings(ipos))) and $
(aperture eq strtrim(aper(ipos))) then begin
bck_a = back_a(ipos) & bck_b = back_b(ipos)
bck_c = back_c(ipos) & bck_d = back_d(ipos)
end
end
end
end else begin
;
; no table supplied, use input parameter values
;
bck_a = 1.0
bck_b = 1.0
bck_c = 0.0
bck_d = 0.0
end
;
; change any user supplied values
;
if strtrim(bck_scale(0)) ne '' then bck_a = float(bck_scale(0))
if strtrim(bck_scale(1)) ne '' then bck_b = float(bck_scale(1))
if strtrim(bck_scale(2)) ne '' then bck_c = float(bck_scale(2))
if strtrim(bck_scale(3)) ne '' then bck_d = float(bck_scale(3))
; update log
;
hist = strarr(3)
hist(0) = ' Interorder background scale factors:'
hist(1) = ' BACK_A ='+string(bck_a,'(F7.4)') + $
' BACK_C ='+string(bck_c,'(F7.4)')
hist(2) = ' BACK_B ='+string(bck_b,'(F7.4)') + $
' BACK_D ='+string(bck_d,'(F7.4)')
if n_elements(log) gt 0 then sxaddhist,hist,log
if !dump gt 0 then printf,!textunit,hist
;
; substep information needed to apply scale factors
;
s = size(flux) & ns=s(1)
nxsteps = ns/500 ;number of xsteps
step_index = indgen(500)*nxsteps ;every nxstep point
end
;
; set up output arrays -------------------------------------------------------
;
nreads = n_elements(first) ;number of readouts
ihb = intarr(128,nreads)
back = fltarr(500,nreads)
epsb = intarr(500,nreads)
not_found=intarr(nreads) ;flags for readouts with missing back.
;
; weights for computing each diode from left and right background diodes
;
if type ge 2 then begin
index=findgen(500)
wright=(index+16.0)/533. ;distance betwen bck diodes=534
wleft=(517-index)/533.
endif
;
; wieghts for expanding background for merged flux
;
if nmerge gt 1 then begin
index=indgen(nmerge*500);vector of data point numbers
index1=index/nmerge ;first diode to interpolate
index2=(index1+1)<499 ;seconde diode to interpolate
weight2=(index mod nmerge)/float(nmerge)
weight1=1.0-weight2
endif
;
; loop on readouts
;
for iread=0,nreads-1 do begin
pos1=first(iread) ;first bin position
pos2=pos1+nf(iread)-1 ;last bin positions
binids=binids_all(pos1:pos2)
notfound=0 ;flag for case where background bin
; missing
;
; Main diode array background bins --------------------------------------------
;
if type lt 2 then begin
if type eq 0 then good=where((binids eq 5) or (binids eq 6)) $
else good=where((binids eq 3) or (binids eq 4))
ngood=!err
if ngood gt 0 then begin
backpos=pos1+good ;background bin positions
;
; extract background bins
;
interorder_present = bytarr(2) ;upper/lower flags
for i=0,ngood-1 do begin
case binids(good(i)) of
3: begin ;upper inter.
interorder_present(0) = 1
bupper = data(6:505,backpos(i))
eupper = eps(6:505,backpos(i))
end
4: begin ;upper inter.
interorder_present(1) = 1
blower = data(6:505,backpos(i))
elower = eps(6:505,backpos(i))
end
else: begin
b = data(6:505,backpos(i))
e = eps(6:505,backpos(i))
end
endcase
endfor
;
; apply the Jason Cardelli fudge factors for type=1 (binids 3 and 4)
;
if type eq 1 then begin
;
; compute object rate per diode
;
object = flux(*,iread)
o = fltarr(500)
for i=0,nxsteps-1 do $
o = o + object(step_index+i)
o = o/nxsteps ;object for each diode
;
; compute average uncorrected background (UB)
;
if interorder_present(0) then begin
ub = bupper
e = eupper
end else begin
ub = fltarr(500)
e = replicate(10000L,500) ;no data
end
if interorder_present(1) then begin
same = where(elower eq e) ;same eps
if !err gt 0 then ub(same) = $
(ub(same)+blower(same))/2.0
better = where(elower lt e)
if !err gt 0 then begin ;better data
ub(better) = blower(better)
e(better) = elower(better)
end
end
;
; compute uncorrected net
;
net = o - ub
ave_net = total(net)/500.0
;
; compute corrected background
;
if interorder_present(0) then begin
b = bupper*bck_a
e = eupper
end else begin
b = fltarr(500)
e = replicate(10000L,500) ;no data
end
if interorder_present(1) then begin
same = where(elower eq e) ;same eps
if !err gt 0 then b(same) = $
(b(same)+blower(same)*bck_b)/2.0
better = where(elower lt e)
if !err gt 0 then begin ;better data
b(better) = blower(better)*bck_b
e(better) = elower(better)
end
end
b = b - bck_c*net + bck_d*ave_net
;
; compute epsilons
;
if interorder_present(0) then e = eupper $
else e = intarr(500)
if interorder_present(1) then e = e>elower
end; if type eq 1
;
; scale sky taken in opposite aperature
;
case binids(good(0)) of
5: b=b*norm5 ;normalize sky for dif-
6: b=b*norm6 ; ferent aperture size
else:
endcase
end else notfound=1
;
; BACKGROUND DIODES
;
end else begin
if type eq 2 then begin
good=where(binids gt 6)
ngood=!err
end else begin
good=where((binids gt 0) and (binids lt 3))
ngood=!err
end
if ngood gt 0 then begin
values=fltarr(ngood*4) ;background values
right=intarr(ngood*4) ;right hand side flag
left=intarr(ngood*4) ;left hand side flag
nv=0 ;number of values
eps_v=fltarr(ngood*4) ;their data quality
for i=0,ngood-1 do begin ;ext. background diodes
bin=pos1+good(i)
d=data(*,bin) ;data vector
e=eps(*,bin) ;epsilon vector
if type eq 3 then begin
;from gross spectral bin
diodes=[0,1,510,511] ;use all four
end else begin
binid=binids(good(i))
case binid of
8: diodes=[0,510] ;upper background
9: diodes=[1,511] ;lower background
10: diodes=0 ;upper left
11: diodes=1 ;lower left
12: diodes=510 ;upper right
13: diodes=511 ;lower right
14: diodes=[0,1] ;both left
15: diodes=[510,511] ;both right
endcase
endelse
values(nv)=d(diodes)
eps_v(nv)=eps(diodes)
right(nv)=diodes gt 256
left(nv)=diodes lt 256
nv=nv+n_elements(diodes)
end; for i loop on bins
;
; separate left and right background diodes
;
left=where(left) & nleft=!err
if nleft gt 0 then begin
goodleft=where(eps_v(left) eq 0)
nleft=!err
if nleft gt 0 then left=left(goodleft)
endif
right=where(right) & nright=!err
if nright gt 0 then begin
goodright=where(eps_v(right) eq 0)
nright=!err
if nright gt 0 then right=right(goodright)
endif
if nright gt 0 then $
right_val=total(values(right))/nright $
else notfound=2
if nleft gt 0 then $
left_val=total(values(left))/nleft $
else notfound=3
if (nleft gt 0) then begin
if (nright gt 0) then $
b=right_val*wright + left_val*wleft $
else b=replicate(left_val,500)
end else begin
if (nright gt 0) then $
b=replicate(right_val,500) $
else notfound=1
endelse
e = 0
end else notfound=1
end;if itype ne 2
;
; Place results in output arrays
;
ihb(0,iread)=ihf(*,iread)
if notfound eq 1 then begin
back(*,iread)=0
epsb(*,iread)=fill
end else begin
if type eq 1 then back(0,iread)=ub else back(0,iread)=b
epsb(0,iread)=e
end
not_found(iread)=notfound
;
; smooth background if type 0 or 1
;
if (notfound eq 0) and (type lt 2) then $
bck_smooth,median_width,mean_width,poly_order,b
;
; expand to same size as flux array for half or quarter stepped data
; and subtract from gross
;
if notfound ne 1 then begin
if nmerge gt 1 then b=b(index1)*weight1+b(index2)*weight2
flux(0,iread)=flux(*,iread)-b
ihf(66,iread)=ihf(66,iread) or 64 ;flag as done
endif
end; iread loop on readouts
;
; print errors to processing log
;
bad=where(not_found gt 0)
nbad=!err
if nbad gt 0 then begin
hist=strarr(nbad)
for i=0,nbad-1 do begin
iread=bad(i)
case not_found(iread) of
1: error=': No valid background data found'
2: error=': No valid right hand side background diodes'
3: error=': No valid left hand side background diodes'
endcase
hist(i)=' Readout'+string(iread+1,'(i5)')+error
endfor
if n_elements(log) gt 0 then sxaddhist,hist,log
if !dump gt 0 then printf,!textunit,hist else print,hist
endif
return
end