%  "Method 1":  We compute the volume of the unit sphere
%  as 8 times the volume of the piece of it in the 1st orthant
%  (the points where each of the four coordinates are nonnegative).
%   Call this piece Omega.
%
%  We can estimate the volume of Omega by uniformly sampling points 
%  in the unit box and counting the proportion of them that are
%  also in Omega.
% 
%  Let N be the number of sample points.
% challenge2.m  Dianne P. O'Leary

N = 1000     % as an example

%  Construct N samples, 3 coordinates each.

x = rand(N,3);

%  Compute} $x_1^2 + x_2^2 + x_3^2 - 1$ {\tt for each of them.
%  (This uses the vector capabilities of Matlab and is faster
%   than doing it component by component.)

y = x(:,1).^2 + x(:,2).^2 + x(:,3).^2 - 1;

%  Count the number of points in the region, and divide by N
%  to get the estimate of the volume of Omega.

Omegavol = sum(y < 0) / N;

%  Multiply by 8 to estimate the volume of the sphere.

Spherevol = 8 * Omegavol


%  "Method 2":  We can compute the volume of the unit sphere
%  as 2 times the integral from y=-1 to y=1 of 
%  the integral from x=-sqrt(1-y^2) to x= sqrt(1-y^2)
%  of sqrt(1-x^2-y^2).  
%  (see pointer)
% 
%  In other words, we need the integral of sqrt(1-x^2-y^2)
%  over the domain that is the unit circle.
% 
%  Let N be the desired number of sample points.

N = 1000 % as an example

%  Construct (approx.) N samples, uniformly distributed in the unit
%  quarter circle.
%  (Since the volume of the unit box is 1 and the volume of the quarter
%  circle is pi/4, we need to sample about 4N/pi points in the box
%  to get about N in the quarter circle.)

x = rand(ceil(4*N/pi),2);  % These are uniformly distributed in the square.
test = x(:,1).^2 + x(:,2).^2 -1;
[good_indices,z] = find(test < 0);
x = x(good_indices,:);    % These are uniformly distributed in the quarter circle.

%  Find the average function value.

z = sqrt(1-x(:,1).^2 -x(:,2).^2);
averagez = sum(z)/length(z);

%  The estimate of the integral is 8 times averagez times the area of the
%  quarter circle.

Spherevol = 2 * averagez * pi


