#
# qtrap
#
# Simple interface for trapzd. Returns estimate of integral of f from a to b,
# using the trapezoidal rule.
#
# The (global) variable QTRAP_EPS can be set to desired fractional accuracy
# and QTRAP_JMAX so that 2^(QTRAP_JMAX-1) is maximum allowed number of steps.
#
# Error 100 indicates that QTRAP_JMAX is too low to get desired accuracy.
#
# This routine is from "Numerical Recipes in C" by Press et. al.
# It uses trapzd(), defined in trapzd.ic.
#
# mws, 7/92. 
#
silent
echo "Reading trapzd..."
read "trapzd.ic"
echo "Defining qtrap..."
#
QTRAP_EPS = 1.0e-5
QTRAP_JMAX = 20
#
func qtrap(~f,a,b) = {
	local j,s,olds

	olds = -1.0e30;
	for (j = 1; j <= QTRAP_JMAX; j += 1) {
		s = trapzd(f,a,b,j)
		if (abs(s-olds) < QTRAP_EPS*abs(olds)) return s
		olds = s
	}
	error(99)	# too many steps
}

echo "Done."
verbose