function V = VolumeCVL(u,v,f,fVals,tol,Level)
% u and v are row 3-vectors and (u(1),v(1)), (u(2),v(2)), and 
%    (u(3),v(3)) are the vertices of a triangle T.
% f is a function of two variable and f(x,y) is defined on T.
% fvals is a row 3-vector where fvals(i) = f(x(i),y(i)), i=1:3
% tol is a positive real number
% V is an approximation of the integral of f over T.
% Level is the number of allowable levels of subdivision.

% Centroid...
xc = mean(u); yc = mean(v); fc = f(xc,yc);
% Compute the volume over triangle with vertices 
%     (xc,yc)), (u(2),v(2)), (u(3),v(3))...
u1 = [xc u(2) u(3)]; 
v1 = [yc v(2) v(3)]; 
fVals1 = [fc fVals(2) fVals(3)];
Vol1 = Area(u1,v1)*sum(fVals1)/3;
% Compute the volume over triangle with vertices 
%     (u(1),v(1)), (xc,y(c) (u(3),v(3)),...
u2 = [u(1) xc u(3)]; 
v2 = [v(1) yc v(3)]; 
fVals2 = [fVals(1) fc fVals(3)];
Vol2 = Area(u2,v2)*sum(fVals2)/3;
% Compute the volume over triangle with vertices 
%     (u(1),v(1)), (u(2),v(2)), (xc,yc)...
u3 = [u(1) u(2) xc]; 
v3 = [v(1) v(2) yc]; 
fVals3 = [fVals(1) fVals(2) fc];
Vol3 = Area(u3,v3)*sum(fVals3)/3;
% Compare V_better to V...
A = Area(u,v);
V = A*(sum(fVals)/3);
V_better = Vol1+Vol2+Vol3;
    
if Level==0 || abs(V-V_better)<=tol*A
   % V_better is good enough...
   V = V_better;
else
   % Recur on three subtriangles..
   Vol1 = VolumeCVL(u1,v1,f,fVals1,tol,Level-1);
   Vol2 = VolumeCVL(u2,v2,f,fVals2,tol,Level-1);
   Vol3 = VolumeCVL(u3,v3,f,fVals3,tol,Level-1);
   V = Vol1+Vol2+Vol3;
end

function A = Area(u,v)
% u and v are row 3-vectors
% A is the area of the triangle with vertices (u(1),v(1)), (u(2),v(2)),
% (u(3),v(3));
d12 = sqrt((u(1)-u(2))^2+(v(1)-v(2))^2);
d23 = sqrt((u(2)-u(3))^2+(v(2)-v(3))^2);
d31 = sqrt((u(3)-u(1))^2+(v(3)-v(1))^2);
s = (d12+d23+d31)/2;
A = sqrt(s*(s-d12)*(s-d23)*(s-d31));
