#
# qsimp
#
# Simple interface for trapzd. Returns estimate of integral of f from a to b,
# using Simpson's rule.
#
# The (global) variable QSIMP_EPS can be set to desired fractional accuracy
# and QSIMP_JMAX so that 2^(QSIMP_JMAX-1) is maximum allowed number of steps.
#
# Error 101 indicates that QSIMP_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 qsimp..."
#
QSIMP_EPS = 1e-6
QSIMP_JMAX = 20
#
func qsimp(~f,a,b) = {
	local j,s,st,ost,os

	ost = os = -1e30
	for (j=1;j<=QSIMP_JMAX;j+=1) {
		st=trapzd(f,a,b,j)
		s=(4*st-ost)/3
		if (abs(s-os) < QSIMP_EPS*abs(os)) return s
		os=s
		ost=st
	}
	error(101)
}

echo "Done."
verbose