Viewing contents of file '../idllib/ghrs/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
;
;*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.
;
;*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:
; None.
;
;*SUBROUTINES CALLED:
; TABINV
; PARSHIFT
;
;*FILES USED:
; None.
;
;*SIDE EFFECTS:
; None.
;
;*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.
;
;*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
; Mar 4 1991 JKF/ACC - moved to GHRS DAF (IDL Version 2)...added
; longword for GHRS spectra.
;-
;*****************************************************************************
pro bins,wave,flux,weight,wcenter,widths,wmean,wsigma,wgt
;
; 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'
retall & end
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
end
;
s =size(wcen)
if s(0) lt 1 then wcen= fltarr(1) + wcen ;convert scalar to array
n = long(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
end else begin
nmax = 0 ; set nmax = 0 for vectors
wmean = fltarr(n) ; output vector of binned fluxes
end
wsigma = wmean
wgt = wmean
iwmean = fltarr(n)
iwsigma = iwmean
iwgt = iwmean
;
; begin loop for arrays
;
for j=0L,nmax do begin
if nmax gt 0 then begin
iflux = flux(*,j)
iwei = wei(*,j)
end else begin
iflux = flux
iwei = wei
end
;
; setup quantities for integration
;
ind = lindgen(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 = long(lin) & r= lin-ind
lw = w(ind)*(1-r) + w(ind+1)*r
lind = long(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 = long(uin) & r= uin-ind
uw = w(ind)*(1-r) + w(ind+1)*r
uind = long(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)) )
end ; 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
;
if !noprint eq 0 then begin ; printout
print,j,' Weighted Values'
fm = "(1x,f10.2,f7.2,2e12.3,f10.2)"
if nparm eq 8 then begin
print,' Wcenter Widths Mean Sigma Weight'
for i=0,n-1 do $
print,form=fm,wcent(i)+wscale,wid(i),iwmean(i),iwsigma(i),iwgt(i)
end else begin
print,' Wcenter Widths Mean Sigma Npnts'
for i=0,n-1 do $
print,form=fm,wcent(i)+wscale,wid(i),iwmean(i),iwsigma(i),npoints(i)
end
print,' '
end ; printout
;
; define output parameters
;
if nmax gt 0 then begin
wmean(0,j) = iwmean
wsigma(0,j) = iwsigma
wgt(0,j) = iwgt
end else begin
wmean = iwmean
wsigma = iwsigma
wgt = iwgt
end
end ; array loop
;
; 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)
end
;
; 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
end
wcenter = wcen
widths = wid
;
return ; bins
end