Viewing contents of file '../idllib/iuedac/iuelib/pro/bfit.pro'
;+******************************************************************************
;
;*NAME:
;
; bfit
;
;*CALLING SEQUENCE:
;
; BFIT,W,F,Q,B,lines=lines,nfirst=nfirst
;
;*PURPOSE:
;
; To create a background vector from the SILO image using the same method
; as SWET (the extraction step of NEWSIPS).
;
;*PARAMETERS:
;
; W (req) (i) (1) (f)
; Wavelengths from SILO file.
;
; F (req) (i) (2) (f)
; Fluxes from SILO file.
;
; Q (req) (i) (2) (i)
; Quality flags from SILO file.
;
; B (req) (o) (1) (f)
; Output background vector.
;
; lines (key) (i) (1) (i)
; Indices of lines to be used for background.
;
; nfirst (key) (i) (0) (i)
; The first NFIRST elements of LINES will be considered part of the
; FIRST swath; remaining LINES will be part of the SECOND swath.
; Default is seven.
;
;*SUBROUTINES CALLED:
;
; parcheck
;
;*FILES USED:
;
; None.
;
;*NOTES:
;
; Although this code has been optimized in the sense that unused
; variables and oversized arrays appearing in the FORTRAN code
; have been dispensed with (and the capabilities of IDL are used),
; the logic is the same as in the original. One place this might
; lead to problems is in the patching of bad points in column 1 of
; the background array. If the very first point in either swath
; is bad, it is replaced with its neighbor in the previous row---
; without checking that point first. Bad points in other columns
; are replaced with the point in the preceding column. This could
; theoretically propagate a bad pixel along the row if there were
; many bad pixels next to each other.
;
;*PROCEDURE:
;
; This is an IDL version of the subroutine IUEBKG (which is a subroutine
; of SWET---Signal Weighted Extraction Technique, used by NEWSIPS to
; create the MXLO file). It first replaces points in the two background
; swaths with the value of the closest (shortward) good point. The two
; regions of patched fluxes are summed separately and then individual
; points in the difference spectrum (region 1 minus region 2) are assigned
; tiny weights if they differ from the mean value (of the d. spectrum) by
; more than 2*sigma. (Points outside the target ring---i.e., all of the
; contiguous bad region at the edge of the spectral lines, are also given
; low weights.) Then a Chebyshev fit is done to the average of all the
; patched background lines. There are comments in this code showing
; which sections are analogous to which of IUEBKG's subroutines.
;
;*MODIFICATION HISTORY:
; Written by LLT.
; 2 Dec 94 LLT fix smoothing, background patching, weighting.
;
;-******************************************************************************
pro bfit,w,f,q,cb,lines=lines,nfirst=nfirst
npar=n_params(0)
if npar eq 0 then begin
print,'bfit,w,f,q,b,lines=lines,nfirst=nfirst'
retall
endif ;npar eq 0
w=float(w)
s=size(f)
npts=s(1)
nlines=s(2)
;Get line numbers of background region lines. These are for large aperture.
if keyword_set(lines) then r=lines else begin
r=indgen(14)+31
r(7)=indgen(7)+63
endelse ;keyword_set(lines)
nr=n_elements(r)-1
if not keyword_set(nfirst) then nfirst=7
;The following steps are taken care of by subroutine SUMBKG when IUEBKG is
;run (SWET).
;Find the indices of points with unacceptable quality flags, and make a copy
;of the flux array with those bad points set to zero.
bad=(q le -256) or ((abs(q) or (not 64)) eq -1)
;If there are bad points in the background region lines, replace them with the
;most recent previous good point value. Start with the points in the first
;column, replacing bad points with the closest previous good point in that
;column, only for background lines (as is done in SUMBKG).
t=f
mgood=0
for ii=0,nr do begin
i=r(ii)
bittest,abs(q(*,i)),6,res
res=res or (q(*,i) le -256)
if res(0) then t(0,i)=t(0,i-1)
for j=1,npts-1 do if res(j) then t(j,i)=t(j-1,i)
wres=where(res eq 0,nwres)
nlast=max(wres)
mgood=nlast > mgood
wres=wres(nwres-20:nwres-1)
pave=total(t(wres,i))/20.
for j=nlast+1,npts-1 do t(j,i)=pave
endfor ; i=0,nr
;SUMBKG is called twice by IUEBKG (once for each background region). In IDL
;all the rows are considered together (for discarding bad points), then the
;two summed background spectra are calculated.
b1=t(*,r(0)) ;First background region
for i=1,nfirst-1 do b1=b1+t(*,r(i)) ;Sum up lines
b2=t(*,r(nfirst)) ;Second background region
for i=nfirst+1,nr do b2=b2+t(*,r(i)) ;Sum up lines
;The smoothed difference between the sums of the two background regions will be
;calculated. The first 320 elements of this difference spectrum will be
;considered separately from the remainder. The average and standard deviation
;is found, and any element where the difference minus the average is greater
;than 2*sigma has its inverse weight set to 1e10. This step is performed by
;AVEBKG.
iw=replicate(1.,npts) ;Inverse weight array
d=b1-b2 ;Difference spectrum
;This is from the subroutine SMOOTH, called by AVEBKG:
for i=0,npts-1 do $ ;Smooth (not like IDL smooth)
d(i)=total(d((0>(i-3)):((npts-1)<(i+3))))/(7+(0<(i-3))+(0<(npts-4-i)))
;Back to AVEBKG:
half=npts/2
av1=total(d(0:half-1))/half ;Average of first half points
sg1=(total((d(0:half-1)-av1)^2)/half)^0.5 ;Standard deviation thereof
av2=total(d(half:mgood))/(mgood-half+1) ;Average of remaining points
sg2=(total((d(half:mgood)-av2)^2)/(mgood-half+1))^0.5 ;Standard deviation
bd1=where(abs(d(0:half-1)-av1) gt 2.0*sg1,nb1) ;"Bad" points of first half
bd2=where(abs(d(half:mgood)-av2) gt 2.0*sg2,nb2) ;"Bad" points of remainder
if nb1 gt 0 then iw(bd1)=1.e10 ;Set inverse weights of bad
if nb2 gt 0 then iw(bd2+half)=1.e10 ;points to 1.e10, discard
iw=iw(0:mgood) ;points past target ring.
;Back from AVEBKG to IUEBKG briefly.
b=(b1+b2)/(nr+1) ;Average background spectrum
;Chebyshev polynomials expect input to be in the range -1 to +1. Therefore,
;create appropriate arrays and rescale the background vector. This is BASISFIT.
x=findgen(mgood+1)/(mgood+1)*2-1
bmax=max(abs(b))
b=b/bmax
;Get Chebyshev coefficents. Here, BASISFIT calls LFIT.
a=fltarr(7) ;Solution vector
covar=fltarr(7,7) ;Covariance matrix
cheb=fltarr(mgood+1,7)
;This is where LFIT calls FUNCS to calculate the Chebyshev coefficients.
cheb(*,0)=1
cheb(*,1)=x
for j=2,6 do cheb(*,j)=2.*x*cheb(*,j-1)-cheb(*,j-2)
;Back to LFIT. Initial values for covariance matrix and solution vector.
bm=b(0:mgood)
for j=0,6 do begin
for k=0,j do covar(j,k)=total(cheb(*,j)*cheb(*,k)/iw^2)
a(j)=total(bm*cheb(*,j)/iw^2)
endfor ;j=0,6
for i=1,6 do covar(0:i-1,i)=covar(i,0:i-1) ;Above diagonal = below diagonal
;Gauss-Jordan elimination, done by subroutine GAUSSJ.
iss1=intarr(7) ;Indices for first subscript (columns in IDL)
iss2=intarr(7) ;Indices for second subscript (rows in IDL)
for i=0,6 do begin
m=where(covar eq max(covar),nm) ;one-dimensional index of maximum value
m=m(nm-1) ;Use last value (if duplicates exist)
ss1=fix(m/7) ;Subscript 1 (column number)
ss2=m mod 7 ;Subscript 2 (row number)
if ss1 ne ss2 then begin ;If not on diagonal, switch first
dummy=covar(ss1,*) ;subscript (columns). This will be the
covar(ss1,*)=covar(ss2,*) ;"pivot" column.
covar(ss2,*)=dummy
dummy=a(ss1) ;Switch around values in A array
a(ss1)=a(ss2) ;so they match COVAR
a(ss2)=dummy
endif ; ss1 ne ss2
iss1(i)=ss1 ;Store first coordinate of pivot elt.
iss2(i)=ss2 ;Store second coordinate of pivot elt.
piv=1./covar(ss2,ss2) ;Divide "pivot" column by max. value
covar(ss2,ss2)=1. ;Preserve maximum value---don't divide
covar(ss2,*)=covar(ss2,*)*piv ;it by itself.
a(ss2)=a(ss2)*piv ;Scale corresponding A accordingly
others=where(indgen(7) ne ss2) ;All columns EXCEPT "pivot" column
prow=covar(others,ss2) ;will be reduced. Same for A.
covar(others,ss2)=0.
for j=0,5 do $
covar(others(j),*)=covar(others(j),*)-prow(j)*covar(ss2,*)
a(others)=a(others)-a(ss2)*prow
endfor ;i=0,6
;Skipping things that don't affect the results (COVSRT etc).
;Now, BASISFIT calls a function CHEBEV to generate the actual background fit.
d=x*0.
dd=d
for j=6,1,-1 do begin
sv=d
d=x*2.*d-dd+a(j)
dd=sv
endfor ;j=6,0,-1
cb=[(x*d-dd+a(0))*bmax,fltarr(npts-mgood)]
;Back to IUEBKG.
if (mgood+1) lt npts then cb(mgood+1:*)=cb(mgood)
;SWET itself: cb=cb*13 (scaling to object fluxes) not done in BFIT.
return
end ;bfit