Viewing contents of file '../idllib/iuedac/iuelib/pro/bins.pro'
;**************************************************************************
;+
;*NAME:
;
; BINS (General IDL Library 01) 6-JAN-83
;
;*CLASS:
;
; Resampling
;
;*CATEGORY
;
;*PURPOSE:
;
; To bin flux data on a specified wavelength grid with or without weights.
;
;*CALLING SEQUENCE:
;
; BINS,WAVE,FLUX,WEIGHT,WCENTER,WIDTHS,WMEAN,WSIGMA,WGT,$
; file=file,/silent,/overwrite
;
;*PARAMETERS:
;
; WAVE (REQ) (I) (1) (I L F D)
; Required input vector giving the wavelength scale for the
; flux data which are to be binned.
;
; FLUX (REQ) (I) (1 2) (I L F D)
; Required input vector or array giving the flux data to be binned.
; If flux is an array (e.g. an ELBL file) the binning is only
; performed in the wavelength direction
;
; WEIGHT (OPT) (I) (1 2) (F D)
; Optional input vector or array giving point by point weights
; for the flux vector. If specified as a scalar, WEIGHT =
; FLUX * 0 + WEIGHT; otherwise, WEIGHT must have same dimensions
; as FLUX.
;
; WCENTER (REQ) (I/O) (0 1) (I L F D)
; Required input scalar or vector giving the centers for
; the wavelength bin(s). The units MUST be the same as for
; the WAVE vector. This vector is recomputed for output to the
; actual value used (i.e. to describe any edge effects).
;
; WIDTHS (REQ) (I/O) (0 1) (I L F D)
; Required input scalar or vector giving the full width in the
; same units as the WAVE array for each bin. This is
; recomputed to the actual value used (i.e. to describe bins
; truncated at ends of array).
;
; WMEAN (REQ) (O) (0 1 2) (F D)
; Required scalar, vector or array output variable containing the
; (weighted) mean value of the FLUX array in each bin.
;
; WSIGMA (REQ) (O) (0 1 2) (F D)
; Required scalar, vector or array output variable containing the
; (weighted) rms standard deviation from the WMEAN in each bin.
;
; WGT (REQ) (O) (0 1) (F D)
; Required scalar or vector output variable containing the
; new weight values in each bin. If weight is not specified,
; WGT equals the number of points in each bin.
;
; FILE (KEY) (i) (0) (IS)
; This keyword allows you to write your results to a file. You
; set it equal to a filename (e.g., file='mybins.txt') or use a
; default filename BINS.TXT by setting this keyword to 1
; (i.e., file=1 or /file).
;
; silent (key) (i) (0) (i)
; If set, results will not be displayed on the screen.
;
; overwrite (key) (i) (0) (i)
; If set, and the results are to be written to a file, the program
; will not check for prior existence of the file before writing.
;
;*EXAMPLES:
;
; with WEIGHT assigned by procedure WEIGHT,
; BINS,WAVE,FLUX,WEIGHT,WCENTER,WIDTHS,WMEAN,WSIGMA,WGT
;
; For uniform weighting:
; BINS,WAVE,FLUX,WCENTER,WIDTHS,WMEAN,WSIGMA,WGT
;
; For weighted binning of a 2-D line-by-line file using 50 A
; bins from 1150 to 1600 angstroms:
; READFILE,IMAGET,LABEL,H,W,FIMAGE,EIMAGE ;create image array
; WEIGHT,EIMAGE,WEIGHT ;use epsilons for weight
; WIDTHS = 50 ;specify constant width
; WCENTER = 1150 + INDGEN(10)*50 ;generate wavelength centers
; BINS,W,FIMAGE,WEIGHT,WCENTER,WIDTHS,WMEAN,WSIGMA,WGT ;run bins
;
;*SYSTEM VARIABLES USED:
;
; !NOPRINT if=0, the results are printed in tabular form
; >0, the results are not printed.
;
;*INTERACTIVE INPUT:
;
;*SUBROUTINES CALLED:
;
; TABINV
; PARSHIFT
; cpychk (Checks for prior existence of file and queries user if nec.)
;
;*FILES USED:
;
;*SIDE EFFECTS:
;
;*RESTRICTIONS:
;
; modified for sun/unix idl version 1.1
;
;*NOTES:
;
; Actually uses (WEIGHT>0) for weights.
; All inputs may be an array, vector or scalar.
; If WEIGHT is a scalar then all points have this weight.
; If WIDTHS is a scalar then all bins have this width.
; The results are printed out if !NOPRINT = 0.
; WCENTER and WIDTHS are recomputed to be the actual values used.
; Typing BINS without parameters will display the procedure call
; statement.
; FLUX and WEIGHT may be specified as arrays so that line-by-line
; type files can binned using weights.
;
; tested with IDL Version 2.1.0 (sunos sparc) 19 Jun 91
; tested with IDL Version 2.1.0 (ultrix mispel) N/A
; tested with IDL Version 2.1.0 (vms vax) 19 Jun 91
;
;*PROCEDURE:
;
; Treats original data rigorously as binned data.
;
;*MODIFICATION HISTORY:
;
; Jan 6 1983 RJP GSFC Modified BINS to allow weights
; Apr 17 1987 RWT GSFC Added PARCHECK, made WEIGHT optional
; May 8 1987 RWT GSFC Fix PARSHIFT
; Jul 23 1987 RWT GSFC Use !NOPRINT to suppress tabular printout
; Aug 25 1987 RWT GSFC correct calculation of N for scalar WCEN,
; and add listing of procedure call statement
; Mar 7 1988 CAG GSFC Add VAX RDAF-style prolog
; 5-23-88 RWT fix scaling factor error in printout format
; 7-28-88 RWT allow FLUX to be an array
; 1-17-89 RWT remove recalculation of total weight (IWGT)
; 5-16-89 RWT & jtb remove LOOKUP commands, use N_ELEMENTS,
; and use vector subscripts in TOTAL commands
; sep 13 1989 jtb @gsfc correct calculation of npoints and weights
; for 'end effects'
; sep 18 1989 jtb @gsfc modifications for unix/sun idl
; apr 23 1990 rwt & jtb correct error for epsilon weighting and
; improve iwsigma, iwmean, and pwgt calculation
; Jun 19 1991 PJL GSFC cleaned up; tested on SUN and VAX; updated prolog
; 10 Mar 1992 LLT print out j+1 instead of j (so line numbers for
; line-by-line files are correct)
; 6 Feb 1995 LLT add file, silent, and overwrite keywords.
;-
;*****************************************************************************
pro bins,wave,flux,weight,wcenter,widths,wmean,wsigma,wgt,file=file,$
silent=silent,overwrite=overwrite
;
; check parameters & redefine if optional parameters are missing
;
nparm = n_params(0)
if nparm eq 0 then begin
print,'BINS,WAVE,FLUX,weight,WCENTER,WIDTHS,WMEAN,WSIGMA,WGT,$'
print,' file=file,/silent,/overwrite'
retall
endif ; nparm
parcheck,nparm,[7,8],'BINS'
wei = weight
wcen = wcenter
if nparm eq 8 then wid = widths
if nparm eq 7 then begin
parshift,1,wei,wcen,wid ; shift input parameters
wei = 1.0
endif ; npar eq 7
;
s =size(wcen)
if s(0) lt 1 then wcen= fltarr(1) + wcen ;convert scalar to array
n = fix(n_elements(wcen))
swid=size(wid)
if swid(0) lt 1 then wid=wid+fltarr(n)
;
; rescale wavelength and flux to avoid trouble with numerical precision
;
wscale = total(wcen) / n ; mean wavelength
wcent = wcen-wscale
w = wave-wscale
scale = total(abs(flux)) / n_elements(flux) ; normalization value
scale = scale+(scale eq 0)
;
; weight = weight>0
;
swei = size(wei)
if (swei(0) lt 1) then wei = (flux*0.0 + wei)
wei = wei > 0.0
;
; determine parameter sizes using dimensions of flux
;
sf = size(flux)
if sf(0) gt 1 then begin
nmax = sf(2) - 1 ; number of rows
wmean = fltarr(n,sf(2)) ; output array of binned fluxes
endif else begin
nmax = 0 ; set nmax = 0 for vectors
wmean = fltarr(n) ; output vector of binned fluxes
endelse ; sf(0)
wsigma = wmean
wgt = wmean
iwmean = fltarr(n)
iwsigma = iwmean
iwgt = iwmean
;
; begin loop for arrays
;
if !noprint or (not keyword_set(silent)) then unit=[-1]
cpy=keyword_set(file)
if cpy then begin
sfile=size(file)
if sfile(n_elements(sfile)-2) ne 7 then file='bins.txt'
if not keyword_set(overwrite) then cpychk,file,cpy
if cpy then begin
openw,ufile,file,/get_lun
if n_elements(unit) eq 1 then unit=[unit,ufile] else unit=[ufile]
endif ;cpy
endif ;keyword_set(file)
for j=0,nmax do begin
if nmax gt 0 then begin
iflux = flux(*,j)
iwei = wei(*,j)
endif else begin
iflux = flux
iwei = wei
endelse ; nmax
;
; setup quantities for integration
;
ind = indgen(n_elements(w))
wdiff = (w(ind+1) - w(ind-1)) / 2.0
area = iflux / scale * wdiff * iwei
sqrs = area * iflux / scale
;
; compute integration limits
;
; low side
;
tabinv,w,wcent-wid/2.,lin
ind = fix(lin) & r= lin-ind
lw = w(ind)*(1-r) + w(ind+1)*r
lind = fix(lin+0.5)
rl = lw - (w(lind) + w(lind-1)) /2.
fl = iflux(lind)/scale
wgtl = iwei(lind)
;
; high side
;
tabinv,w,wcent+wid/2.,uin
ind = fix(uin) & r= uin-ind
uw = w(ind)*(1-r) + w(ind+1)*r
uind = fix(uin+0.5)
ru = (w(uind) + w(uind+1)) /2. - uw
fu = iflux(uind)/scale
wgtu = iwei(uind)
;
; recompute actual areas to be averaged
;
wcent = (uw+lw) /2.
wid = uw - lw
npoints= uind - lind + 1
itemp = iwei * wdiff
;
; integrate by summation
;
for i=0,n-1 do begin
iwmean(i) = total( area(lind(i):uind(i)) )
iwsigma(i) = total( sqrs(lind(i):uind(i)) )
iwgt(i) = total( itemp(lind(i):uind(i)) )
endfor ; integration by summation
;
; correct for end effects
;
lends = wgtl * rl
uends = wgtu * ru
npoints = npoints - rl/wdiff(lind) - ru/wdiff(uind)
iwgt = iwgt - lends - uends
iwmean = iwmean - fl * lends - fu * uends
iwsigma = iwsigma - fl * fl * lends - fu * fu * uends
;
; compute statistics
;
iwsigma = iwsigma * iwgt - iwmean * iwmean
pwgt = (iwgt gt 0) / (iwgt + (iwgt eq 0)) ;ensure non-zero division
iwsigma = sqrt(iwsigma *(iwsigma gt 0)) * pwgt * scale
;
; note that we avoid error sqrt(neg or zero) can occur when npoints
; is small due to the numerical approximations
;
iwmean = iwmean * pwgt * scale
iwgt = iwgt/(wid+(wid le 0)) * (wid gt 0) * npoints ; avg wt * npts
;
; printout option
;
for u=0,n_elements(unit)-1 do begin ; printout
printf,unit(u),j+1,' Weighted Values'
fm = "(1x,f10.2,f7.2,2e12.3,f10.2)"
if nparm eq 8 then begin
printf,unit(u),' Wcenter Width Mean Sigma Weight'
for i=0,n-1 do $
printf,unit(u),form=fm,wcent(i)+wscale,wid(i),iwmean(i),iwsigma(i),iwgt(i)
endif else begin
printf,unit(u),' Wcenter Width Mean Sigma Npnts'
for i=0,n-1 do $
printf,unit(u),form=fm,wcent(i)+wscale,wid(i),iwmean(i),iwsigma(i),npoints(i)
endelse ; nparm
if nmax gt 0 then if j ne nmax then printf,unit(u),' '
endfor ; printout
;
; define output parameters
;
if nmax gt 0 then begin
wmean(0,j) = iwmean
wsigma(0,j) = iwsigma
wgt(0,j) = iwgt
endif else begin
wmean = iwmean
wsigma = iwsigma
wgt = iwgt
endelse ; nmax
endfor ; array loop
if cpy then free_lun,ufile
;
; convert output to scalar if necessary
;
if s(0) lt 1 then begin ; convert output to scalar
wcent= wcent(0) & wid= wid(0)
wmean= wmean(0) & wsigma= wsigma(0) & wgt=wgt(0)
endif ; s(0)
;
; rescale wavelengh scale
;
wcen=wcent + wscale
;
; reset input parameters if weight was not specified
;
if nparm eq 7 then begin
parshift,-1,wei,wcen,wid,wmean,wsigma,wgt
wgt = npoints
endif ; nparm eq 7
wcenter = wcen
widths = wid
;
return ; bins
end ; bins