Viewing contents of file '../idllib/contrib/buie/avger.pro'
;+
; NAME:
; avger
; PURPOSE: (one line)
; Temporal averaging of time-series data.
; DESCRIPTION:
; This program was written to perform N point averaging of raw photometry
; data. It will work on any temporal data streams or data that has clumpy
; independent variable values. The data are grouped together into bin
; that are specified by the THRESH input. Thresh specifies the size of a
; gap that will cause the group to be broken. The value for thresh is taken
; to be a multiple of the 'normal' spacing between points. If THRESH=2,
; then any gap twice as long as the previous point spacing will cause a
; break. Any number equal to or less than 1 will prevent all averaging.
; To prevent too much binning for long uniform data runs, MAXBIN puts an
; upper limit on the number of points that can be grouped together and
; XSPREAD limits the xspan within a single group.
;
; The THRESH criterion is applied to the data first for grouping,
; then MAXBIN and XSPREAD are used simultaneously to break up long
; binning strings.
;
; The data are averaged together using a weighted average (see MEANERR).
; The uncertainty returned is the standard deviation of the mean.
;
; CATEGORY:
; Numerical
; CALLING SEQUENCE:
; pro avger,x,y,err,maxbin,thresh,avgx,avgy,sigy
; INPUTS:
; x - Independent variable.
; y - Dependent variable.
; err - Uncertainty on y in units of standard deviation.
; This can be a scalar or a vector but must not be zero.
; maxbin - Maximum number of points to average together.
; thresh - Gap that will break grouping of data as a fraction of normal.
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; DATAERR - Output vector of sigma of the mean computed directly
; from the scatter in the data. If the number of points
; in an output bin is one, then the output error is the
; same as the input error.
; FORCEIT - 2xN vector, First column is a point number and the second
; column is a flag, 0 means force this point to bin, 1 means
; force a break at this point. Point numbers outside of
; the valid data range are silently ignored. This info
; if supplied overrides the breaking controlled by THRESH,
; MAXBIN, and XSPREAD allowing a direct modification of
; binning for pathalogical cases.
; Example:
; forceit=[[13,0],[14,1],[19,0]]
; would force points 13 and 19 to NOT end the binning
; and would force point 14 to be the end of a bin.
; When using this option, VERBOSE is especially useful.
; XSPREAD - maximum range of x allowed in a single averaged point.
; (default = no limit).
; VERBOSE - Flag, if true will cause a complete printout of how the
; vector is being binned.
; OUTPUTS:
; avgx - X value after binning.
; avgy - Y value after binning.
; sigy - New uncertainty.
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; The input vectors must have equal length and should be greater in length
; than 3. If the vectors are of length 2 or 3, the program will return a
; straight average of all input values. Scalar inputs are not allowed and
; will generate an error.
; PROCEDURE:
; MODIFICATION HISTORY:
; 93/05/11 - Written by Marc W. Buie, Lowell Observatory.
; 94/03/21, MWB, modified to output double precision if input is double.
; 94/02/25, MWB, added XSPREAD control over binning.
; 98/01/16, MWB, added DATAERR keyword
;-
pro avger,x,y,in_err,maxbin,thresh,avgx,avgy,sigy, $
DATAERR=dsigy,FORCEIT=forceit,XSPREAD=xspread,VERBOSE=verbose
if badpar(x,[2,3,4,5],1,caller='AVGER: (x) ',npts=nx,type=xtype) then return
if badpar(y,[2,3,4,5],1,caller='AVGER: (y) ',npts=ny,type=ytype) then return
if badpar(in_err,[0,2,3,4,5],[0,1],caller='AVGER: (err) ', $
npts=nerr,default=1) then return
if badpar(maxbin,[2,3],0,caller='AVGER: (maxbin) ') then return
if badpar(thresh,[2,3,4,5],0,caller='AVGER: (thresh) ') then return
if badpar(xspread,[0,2,3,4,5],0,caller='AVGER: (XSPREAD) ') then return
if badpar(forceit,[0,2,3],2,caller='AVGER: (FORCEIT) ') then return
if nerr eq 1 then err=replicate(in_err,ny) else err = in_err
test_zero=where(err eq 0.,count_zero)
if count_zero ne 0 then begin
print,'AVGER: the input sigma array must not contain zeros.
return
endif
if nx ne ny or nx ne n_elements(err) then begin
print,'AVGER: Error, the input arrays must have the same lengths.'
return
endif
nx = n_elements(x)
if nx eq 1 then $
message,'Input must be vectors longer than one element.'
if not keyword_set(xspread) then xspread=x[nx-1]-x[0]
if nx lt 4 then begin
meanerr,x,err,avgx
meanerr,y,err,avgy,sigy
endif else begin
;Compute the pre and post dx values, doesn't include first and last point.
predx = x[1:nx-2] - x[0:nx-3]
postdx = x[2:nx-1] - x[1:nx-2]
;Find break points for averaging based on dx test.
; First never breaks, last always breaks.
breakit = [0, (postdx gt thresh*predx), 1]
if keyword_set(verbose) then begin
print,'AVGER: breakit after dx test'
print,breakit
endif
; Check length of averaging run and break if length hits maxbin
count=0
istart=0
for i=0,nx-2 do begin
if breakit[i] then begin
count=0
istart=i+1
endif else begin
count=count+1
endelse
if (count ge maxbin) or (x[i+1]-x[istart]) gt xspread then begin
breakit[i] = 1
count=0
istart=i+1
endif
endfor
if keyword_set(verbose) then begin
print,'AVGER: breakit after run length and xspread check'
print,breakit
endif
; Apply FORCEIT if specified.
oldbreakit=breakit
if keyword_set(forceit) then begin
for i=0,n_elements(forceit)/2-1 do begin
idx = forceit[0,i]
flag = forceit[1,i] eq 1
if idx ge 0 and idx lt nx then $
breakit[idx] = flag
endfor
endif
if keyword_set(verbose) then begin
for i=0,nx-1 do $
print,i,oldbreakit[i],breakit[i],x[i],y[i], $
format='(i3,1x,i1,1x,i1,1x,g,1x,g)'
endif
;Find length of output vectors and initialize.
z=where(breakit eq 1, nout)
if xtype ne 5 and ytype ne 5 then begin
avgx = fltarr(nout,/nozero)
avgy = fltarr(nout,/nozero)
sigy = fltarr(nout,/nozero)
dsigy = fltarr(nout,/nozero)
endif else begin
avgx = dblarr(nout,/nozero)
avgy = dblarr(nout,/nozero)
sigy = dblarr(nout,/nozero)
dsigy = dblarr(nout,/nozero)
endelse
; This is the major working loop. Thumb through breakit and average points
; whenever it is true. Store input if length of run is 1. Otherwise,
; average the points.
j=0
nstart=0
for i=0,nx-1 do begin
if breakit[i] then begin
if nstart eq i then begin
avgx[j] = x[i]
avgy[j] = y[i]
sigy[j] = err[i]
dsigy[j] = err[i]
endif else begin
meanerr,x[nstart:i]-x[nstart],err[nstart:i],avgxval
meanerr,y[nstart:i],err[nstart:i],avgyval,sigyval,sigdata
avgx[j] = avgxval+x[nstart]
avgy[j] = avgyval
sigy[j] = sigyval
dsigy[j] = sigdata/float(i-nstart)
endelse
nstart = i+1
j=j+1
endif
endfor
endelse
end