% MonteCarlo calculation of the volume of the unit 2-sphere.
% Uses STRATIFIED SAMPLING.

function volume = twosphere(Npoints)
% Npoints is the total number of points where the function will be
% evaluated. We need this to be an exact square; if not, we'll get
% as close as we can to Npoints with a perfect square.
N = round(Npoints^(1/2)); % Number of abcissas per dimension
Npoints = N^2; % Actual number of points that will be used.

h = 1/N; % step size for dicing the unit cube


% Create square matrices of x and y
x = ones(N,1) * (h*(0:N-1));
y = (h*(0:N-1))' * ones(1,N);

% Random offsets in x,y
xoff = h * rand(N,N);
yoff = h * rand(N,N);
   
% square of distance from the origin for each x,y,z
r2 = (x+xoff).^2 + (y+yoff).^2;
   
Ninside = sum(sum(r2<1)); % sum rows & columns

volume = (Ninside/Npoints)*(2^2);
   % 2^2 is the number quadrants in the 2-sphere. 
