function area = gquad (fun,xlow,xhigh,mparts,bp,wf)
% area = gquad (fun,xlow,xhigh,mparts,bp,wf)
% This function evaluates the integral of an externally
% defined function fun(x) between limits xlow and xhigh. The
% numerical integration is performed using a composite Gauss
% integration rule. The whole interval is divided into mparts
% subintervals and the integration over each subinterval
% is done with an nquad point Gauss formula which involves base
% points bp and weight factors wf. The normalized interval
% of integration for the bp and wf constants is -1 to +1. The
% algorithm is described by the summation relation
% x=b j=n k=m
% integral( f(x)*dx ) = d1*sum sum( wf(j)*fun(a1+d*k+d1*bp(j)) )
% x=a j=1 k=1
% where bp are base points, wf are weight factors
% m = mparts, and n = length(bp) and
% d = (b-a)/m, d1 = d/2, a1 = a-d1
% The base points and weight factors must first be generated
% by a call to grule of the form [bp,wf] = grule(nquad)
% by Howard Wilson, U. of Alabama, Spring 1990
bp=bp(:); wf=wf(:);
d = (xhigh - xlow)/mparts; d2 = d/2; nquad = length(bp);
x = (d2*bp)*ones(1,mparts) + (d*ones(nquad,1))*(1:mparts);
x = x(:) + (xlow-d2); fv = feval(fun,x); wv = wf*ones(1,mparts);
area=d2*(wv(:)'*fv(:));