==== Six Point Composite Gauss Formula ====
==== With Weight Factors Included ====
area = gquad6(fun,xlow,xhigh,mparts)
This function determines the area under 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 a six 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 structured in terms of a parameter mquad = 6 which
can be changed along with bp and wf to accommodate a different
order formula. The composite algorithm is described by the
following 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 d = (b-a)/m, d1 = d/2, a1 = a-d1,
m = mparts, and n = nquad.