usage: vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy)
or
vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,nquadx,nquady)
This function evaluates the integral of an externally
defined function fun(x,y) between limits xlow and xhigh
and ylow and yhigh. The numerical integration is performed
using a Gauss integration rule. The integration
is done with an nquadx by nquady Gauss formula which involves base
point matrices bpx and bpy and weight factor matrix wfxy. The normalized
interval of integration for the bpx, bpy and wfxy constants is -1 to +1
(in x) and -1 to +1 (in y). The algorithm is described by the
summation relation
x=b j=nx k=ny
integral( f(x)*dx ) = J*sum sum( wfxy(j,k)*fun( x(j), y(k) ) )
x=a j=1 k=1
where wfxy are weight factors,
nx = nquadx = number of Gauss points in the x-direction,
ny = nquady = number of Gauss points in the y-direction,
x = (xhigh-xlow)/2 * bpx + (xhigh+xlow)/2 = mapping function in x,
y = (yhigh-ylow)/2 * bpy + (yhigh+ylow)/2 = mapping function in y,
and J = (xhigh-xlow)*(yhigh-ylow)/4 = Jacobian of the mapping.
The base points and weight factors must first be generated
by a call to grule of the form [bpx,bpy,wfxy] = grule2d(nquadx,nquady)
The first form of gquad2d is faster when used several times, because
the points and weights are only calculated once.
The second form of gquad2d is usefull if it is only called once (or a
few times).
