% MonteCarlo calculation of the volume of the unit sphere
% in 3 dimensions. Uses STRATIFIED SAMPLING.

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

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

% xyz offsets for all the points
Ninside = 0; % initialize total of "in" points

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

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

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