procedure wstat (tabfil)

string	tabfil	=""	{prompt="file listing data values and errors"}
char	delim	=" "	{prompt="column separation character"}
int	nval	=0	{prompt="number of data values"}
real	wmean	=0.	{prompt="weighted mean of data values"}
real	wsigma	=0.	{prompt="sigma on weighted mean of data values"}
real	wrms	=0.	{prompt="weighted rms of data values"}
real	werr 	=0.	{prompt="weighted mean of the errors"}
bool	verbose	=no	{prompt="print results to screen?"}

# @(#) task wstat                     Author: R.A. Jansen -- Sep  1 2000
# @(#) 
# @(#) task to calculate the weighted statistics (mean, sigma, rms, and
# @(#) mean of errors) for a set of measurements and the errors thereon.
# @(#) Data and errors should be supplied as a 2-column ASCII text file.
# @(#) Any additional columns are ignored.  A warning is issued after
# @(#) 1024 lines are read.

begin
	string	infil, sep, cstr
	struct	line
	real	rval, rerr, vsum, wsum, dsum, ravg, rsig, rrms
        int	ival

	infil	= tabfil
	sep	= delim
	if ( sep == "" ) { sep = " " }
	
#	Check existence of input text file...
	chkimg (infil, "access", imtype="file", verbose=yes)
	if (!chkimg.ok) { return }

#	Input OK; loop over all lines in the file (skipping commented
#	lines) and calculate the weighted statistics.

#	First pass...
	ival = 0
	vsum = 0.
	wsum = 0.
	dsum = 0.
	list = infil
	while ( fscan(list,line) != EOF ) {
	   cstr = line
	   if ( substr(cstr,1,1) != "#" ) {
	      ival += 1
	      cparse (cstr, delim=sep)
	      rval = real(cparse.field1)
	      rerr = real(cparse.field2)
	      if (rerr == 0.) { rerr = 0.01*rval }
	      if (rerr == 0.) { rerr = 0.01 }
	      vsum = vsum + rval/rerr
	      wsum = wsum + 1.0/rerr
	      dsum = dsum + (1.0/rerr)**2
	      if ( ival>1024 ) { print("WARNING: File > 1024 lines! Continuing...") }
	   }
	}
	list = ""

#	Weighted mean of the data values...
	ravg = vsum/wsum

#	Second pass...
	rsig = 0.
	rrms = 0.
	list = infil
	while ( fscan(list,line) != EOF ) {
	   cstr = line
	   if ( substr(cstr,1,1) != "#" ) {
	      cparse (cstr, delim=sep)
	      rval = real(cparse.field1)
	      rerr = real(cparse.field2)
	      if (rerr == 0.) { rerr = 0.01*rval }
	      if (rerr == 0.) { rerr = 0.01 }
	      rsig = rsig + (1.0/rerr)*(rval-ravg)**2
	      rrms = rrms + ((1.0/rerr)*(rval-ravg))**2
	   }
	}
	list = ""

#	Standard deviation on weighted mean of data values ...
	rsig = sqrt( rsig/(wsum - dsum/wsum) )
#	Weighted rms of data values ...
	rrms = sqrt( rrms/(dsum*(ival-1)) )
#	Weighted mean of error values ...
	rerr = real(ival)/wsum

#	Print results ...
	if (verbose) {
	   print ("  N    wmean   wsigma   werr")
	   printf('%5g  %.4f  %.4f  %.4f\n',ival,ravg,rsig,rerr)
	}

#       update task parameters ...
	nval	= ival
	wmean	= ravg
	wsigma	= rsig
	wrms	= rrms
	werr	= rerr

end