function vol = volMonteCarlo(region,x0,x1,y0,y1,z0,z1,n)
%VOLMONTECARLO estimates the volume of a region in 3-space.
% VOLMONTECARLO(region,x0,x1,y0,y1,z0,z1,n) estimates
% the volume of region between the given limits for x, y, and z
% with a Monte Carlo method with n trials.  The region should
% be specified by a function that takes a matrix with 3 columns
% and outputs a column vector with a 0 if a certain row
% is not in the region and a 1 (or 1 multiplied by an integrand)
% if the row is in the region.  For example one could take
% sphere = @(X) heaviside(ones(size(X,1),1)-X(:,1).^2-X(:,2).^2-X(:,3).^2)

X = rand(n,3);  % random numbers between 0 and 1, now rescale
X(:,1)= x0+(x1-x0)*X(:,1);
X(:,2)= y0+(y1-y0)*X(:,2);
X(:,3)= z0+(z1-z0)*X(:,3);
% At this point each row of X is in the desired range.
vol = sum(region(X))*(x1-x0)*(y1-y0)*(z1-z0)/n;
end