Viewing contents of file '../idllib/ghrs/pro/centerflux.pro'
pro centerflux,idlist,params,table
;
;+
; centerflux
;
; Routine to compute average flux values for the selected portion of
; the diode arrays.
;
; CALLING SEQUENCE:
; centerflux,idlist,params,table
;
; INPUTS:
; idlist - integer list of observation id numbers or
; string array of file names.
;
; OPTIONAL INPUTS:
; params - standard GHRS parameter specification.
; 0 - use all defaults (default if params not supplied)
; 1 - interactive change parameters
; string - 'parname=parvalue,parname=parvalue,...'
; filename - name of a text file containing non defaulted
; parameters, one per line in the form parname=parvalue
;
; The following parameters are available with the specified
; defaults:
; D1=200 ;first diode in range to use
; D2=300 ;last diode to range to use
; DTYPE=NET ;type of data to average
; WAVE=1 ;tabulated versus wavelenghts
; 0 = no 1 = yes
; if 0 then wavelength vectors
; do not need to be present in the
; input data files.
; MAXEPS=0 ;maximum allowed data quality
; flag to accept
; USEEPS=1 ;use epsilon values. If 0 then
; no epsilon values need be present
; in the data files
; IDENT=_R ;file identification concatenated
; to the file names
; MAXNUM=1000 ;number number of flux vectors that
; can be processed.
; ORDER=0 ;order of polynomial to fit to
; data. 0 = use average
; table - output stsdas table name. If not supplied then the
; output tablename will be CENTERFLUX.
;
; HISTORY:
; version 1 D. Lindler April 1989
; 27-JAN-1992 JKF/ACC - converted to IDL V2.
;-
;------------------------------------------------------------------------------
nparms = n_params(0)
if nparms lt 1 then $
message,' Calling Sequence: centerflux,idlist,params,table'
if nparms lt 2 then params=0
if nparms lt 3 then table='centerflux'
;
; change idlist to vector if scalar supplied
;
s=size(idlist) & ndim=s(0) & datatype= vmscode(s(ndim+1))
nfiles=n_elements(idlist)
if ndim gt 0 then begin ;is it already a vector
files=idlist
end else begin
if datatype eq 1 then begin ;is it a string
files=strarr(80)
files(0)=idlist
end else begin ;integer scalar
files=indgen(nfiles)+idlist
endelse
endelse
;
; get default parameters
;
getpar,params,'centerflux',par
d1=fix(par(0))>0
d2=fix(par(1))<499>d1
type=strtrim(par(2))
wave_comp=fix(par(3))
maxeps=float(par(4))
eps_comp=fix(par(5))
ident=strtrim(par(6))
maxnum=fix(par(7))
order = fix(par(8))
;
; create output table arrays
;
if wave_comp then begin
wave=dblarr(maxnum)
wlow=dblarr(maxnum)
whigh=dblarr(maxnum)
endif
npoints=intarr(maxnum)
carpos=lonarr(maxnum)
orders=intarr(maxnum)
flux=fltarr(maxnum)
num=0 ;number of spectra processed
;
; loop on files
;
first=1 ;flag for first line of data
for ifile=0,nfiles-1 do begin
if !dump then print,' Observation ',files(ifile)
;
; open file and read data
;
hrs_open,files(ifile),in,'read',ident
hrs_read,in,type,ih,data
if !err lt 0 then begin
print,'CENTERFLUX-- No data of type ',type,' found'+ $
' for observation ',files(ifile)
retall
endif
s=size(data) & ns=s(1)
if s(0) gt 1 then nspec=s(2) else nspec = 1
if wave_comp then begin
hrs_read,in,'wave',ih2,wave
if !err lt 0 then begin
print,'CENTERFLUX-- No wavelenghts found '+ $
'for observation ',files(ifile)
retall
endif
s=size(wave)
if (ns ne s(1)) then $
message,'CENTERFLUX-- wavelengths/data not '+ $
'compatible for file '+files(ifile)
;
; If 2D, check compatibility...
;
if (s(0) gt 1) and (nspec ne s(2)) then $
message,'CENTERFLUX-- wavelengths/data not '+ $
'compatible for file '+files(ifile)
endif
if eps_comp then begin
hrs_read,in,type+'/eps',ih2,eps
if !err lt 0 then begin
print,'CENTERFLUX -- no epsilons found for'+ $
'observation ',files(ifile)
retall
endif
endif
hrs_close,in
;
; if first file then get grating mode and aperture, else verify
; that the grating mode and aperture are the same
;
gmode=ih(48)
binid=ih(53)
sclamp=ih(39)*2+ih(38)
det=ih(31)
case sclamp of
1: aper='SC1'
2: aper='SC2'
else: if binid eq 2 then aper='LSA' else aper='SSA'
endcase
if first then begin
detector=det
grating=gmode
aperture=aper
first=0
end else begin
if (grating ne gmode) or (aperture ne aper) or $
(det ne detector) then begin
print,'CENTERFLUX-- all observations must be for'+ $
'then same grating, aperture and detector'
retall
endif
endelse
;
; compute starting and ending data point position for region to average
;
; on main array
case ns of
512: xsteps=1
500: xsteps=1
1000: xsteps=2
2000: xsteps=4
else: begin
print,'CENTERFLUX--data length must be '+ $
' 512, 500, 1000 or 2000 points'
retall
end
endcase
i1 = d1*xsteps
i2 = d2*xsteps + (xsteps-1)
;
; loop on spectra within the file
;
for ispec=0,nspec-1 do begin
d = data(i1:i2,ispec)
xx = findgen(i2-i1+1)
;
; find good data points
;
if eps_comp then begin
e=eps(i1:i2,ispec)
good=where(e le maxeps) & ngood=!err
if ngood lt 1 then begin
print,'CENTERFLUX-- all data bad for'+ $
'file ',files(ifile),' spectra',ispec+1
goto,next_ispec
endif
d=d(good)
xx = xx(good)
endif else ngood=i2-i1+1
;
; determine wavelength range
;
if wave_comp then begin
if eps_comp then begin
i1w=i1+good(0) ;first good point
i2w=i1+good(ngood-1) ;last good point
end else begin
i1w=i1
i2w=i2
endelse
wlow(num)=wave(i1w,ispec)
whigh(num)=wave(i2w,ispec)
endif
;
; compute average flux or midpoint of the polynomial
;
npoints(num)=ngood
if order eq 0 then begin
flux(num)=total(d)/ngood
end else begin
coef = poly_fit(xx,d,order)
mid = (i2-i1)/2.0
ff = coef(0)
for ii=1,order do ff = ff + coef(ii)*mid^ii
flux(num) = ff
end
;
; extract other useful parameters from header
;
carpos(num)=ih(43,ispec)
orders(num)=ih(50,ispec)
num=num+1
next_ispec:
endfor; ispec
endfor; ifile
;
; create output table
;
if num lt 1 then begin
print,'CENTERFLUX-- all data determined as unusable'
retall
endif
tab_create,tcb,tab
num=num-1
tab_put,'ORDER',orders(0:num),tcb,tab
tab_put,'CARPOS',carpos(0:num)+(carpos lt 0)*65536L,tcb,tab
tab_put,'WAVELENGTH',(wlow(0:num)+whigh(0:num))/2.0,tcb,tab
tab_put,'FLUX',flux(0:num),tcb,tab
tab_put,'WLOW',wlow(0:num),tcb,tab
tab_put,'WHIGH',whigh(0:num),tcb,tab
tab_put,'NPOINTS',npoints(0:num),tcb,tab
;
; create header
;
h=strarr(80) & h(0)='END'
sxaddpar,h,'DETECTOR',detector
gratings=['NDF','G-1','G-2','G-3','G-4','G-5','E-A','E-B']
grating=gratings(grating)
sxaddpar,h,'GRATING',grating
sxaddpar,h,'APERTURE',aperture
sxaddpar,h,'D1',d1
sxaddpar,h,'D2',d2
if eps_comp then sxaddpar,h,'MAXEPS',maxeps
;
; write the table
;
tab_write,table,tcb,tab,h
return
end