#
# qromb
#
# Performs numerical integration of a function, using Romberg's method.
#
# Error 102 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, and polint(), defined in polint.ic.
#
# mws, 7/92. 
#
silent
echo "Reading polint..."
read "polint.ic"
echo "Reading trapzd..."
read "trapzd.ic"
echo "Defining qromb..."

QROMB_EPS = 1e-6
QROMB_JMAX = 20
QROMB_JMAXP = QROMB_JMAX+1
QROMB_K = 5
array qromb_s[QROMB_JMAXP+1]
array qromb_h[QROMB_JMAXP+1]
array qromb_c[QROMB_K]
array qromb_d[QROMB_K]

func qromb(~f,a,b) = {
	local ss,dss,k,j

	qromb_h[1]=1
	for (j=1;j<=QROMB_JMAX;j+=1) {
		qromb_s[j]=trapzd(f,a,b,j);
		if (j >= QROMB_K) {
			for (k = 1; k <= QROMB_K; k+=1) {
				qromb_c[k] = qromb_h[j-QROMB_K+k]
				qromb_d[k] = qromb_s[j-QROMB_K+k]
			}
			polint(qromb_c,qromb_d,QROMB_K,0,ss,dss)
			if (abs(dss) < QROMB_EPS*abs(ss)) return ss
		}
		qromb_s[j+1]=qromb_s[j]
		qromb_h[j+1]=0.25*qromb_h[j]
	}
	error(102)	# too many iterations required
}

echo "Done."
verbose