%==========================================================================
%
% This function computes quadrature weights for evaluating the surface
% integral of a scalar function f(x,y,z) over the surface S. The surface S
% may be specified implicitly as a level surface h(x,y,z)=0 (along with the
% Quadrature_Nodes and Triangles) with boundary b(x,y,z)>=0.
%
% Inputs: Quadrature_Nodes  A set of points located exactly on the
% surface that can be a (Number_of_Quadrature_Nodes X 3) Array or
% a (Number_of_Quadrature_Nodes X 2) Array
%
% Triangles  A triangulation of the set Quadrature_Nodes. This
% should be an array where row k contains the indicdes in
% Quadrature_Nodes of the vertices of triangle k
%
% Boundary_Edge_Nodes  The set of boundary edges. This should
% be an array where rows 2l1 and 2l contain the two or
% threedimensional coordinates of the vertices of edge l. If
% this is given as an empty array, then the algorithm will try to
% determine the boundary edges by finding those edges for which
% both vertices have b approximately zero.
%
% h  For the surface S defined implicitly by h(x,y,z)=0, row i
% in the output of h should contain
% h(Quadrature_Nodes(i,:))
% h should take in Quadrature_Nodes as specified above
%
% gradh  The gradient of the function h. Row i in the output
% of gradh should contain
% [dh/dx(Quadrature_Nodes(i,:)),dh/dy(Quadrature_Nodes(i,:)),dh/dz(Quadrature_Nodes(i,:)]
% or
% [dh/dx(Quadrature_Nodes(i,:)),dh/dy(Quadrature_Nodes(i,:))]
% corresponding to the specification of Quadrature_Nodes. gradh
% should take in Quadrature_Nodes as specified above.
%
% b  For the boundary of S defined implicitly by b(x,y,z)=0,
% row i in the output of b should contain
% b(Quadrature_Nodes(i,:))
% b should take in Quadrature_Nodes as specified above.
%
% gradb  The gradient of the function b. Row i in the output
% of gradb should contain
% [db/dx(Quadrature_Nodes(i,:)),db/dy(Quadrature_Nodes(i,:)),db/dz(Quadrature_Nodes(i,:)]
% or
% [db/dx(Quadrature_Nodes(i,:)),db/dy(Quadrature_Nodes(i,:))]
% corresponding to the specification of Quadrature_Nodes. gradb
% should take in Quadrature_Nodes as specified above.
%
% options (optional)  an optional structure that is used to
% indicate the types of surfaces and boundaries that are present.
% There are four options that can be specified:
% Planar_Surface_Flag  0 for nonplanar surface (default)
% 1 for planar surface
% Planar_Surface_Normal  the normal to the planar surface
% given as a 1x3 vector. The
% default is [0 0 1]
% Planar_Boundary_Flag  0 for nonplanar boundary (default)
% 1 for planar boundary
% Planar_Boundary_Normal  the normal to the planar bounary
% given as a 1x3 vector. The
% default is [0 0 1]
% Example: To specify that the boundary is planar with normal
% [1/2,1,2] set
% options.Planar_Boundary_Flag=1;
% options.Planar_Boundary_Normal=[1/2,1,2];
% and pass options to this function.
%
% Output: Quadrature Weights  A set of quadrature weights
% (Number_of_Quadrature_Nodes X 1) vectorcorresponding
% to the set of points Quadrature_Nodes
%
% Examples of how to use this function with various options are available
% on Matlab Central's file exchange in Bounded_Smooth_Surface_Quadrature_RBF_Test.m
%
% This implementation uses the method and default settings discussed in:
%
% J. A. Reeger and B. Fornberg "Numerical quadrature over
% smooth surfaces with boundaries".
%
% NOTE: The main loop of this method (over each triangle) can be easily
% parallelized if you have access to the parallel toolbox. In such a case,
% change the for loop to a parfor loop.
%
%==========================================================================
1.3.0.0  Included missing .mat files for test code. 

1.2.0.0  Removed dependence on parallel toolbox. 

1.1.0.0  Updated description. 
Create scripts with code, output, and formatted text in a single executable document.