%   ==== Six Point Composite Gauss Formula ====
%   ====   With Weight Factors Included    ====


function area = gquad6(fun,xlow,xhigh,mparts)
%  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.
%

%     by Howard B. Wilson, U. of Alabama, Spring 1990

%  The weight factors are
wf = [ 1.71324492379170d-01;   3.60761573048139d-01;...
        4.67913934572691d-01]; wf=[wf;wf([3 2 1])];

%  The base points are
bp = [-9.32469514203152d-01;  -6.61209386466265d-01;...
      -2.38619186083197d-01]; bp=[bp;-bp([3 2 1])];

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(:));