clear all
  close all
 
%
% Input number of random points
%

  N = input('N = ')

%
%  Compute uniformly distributed points and Sobol points.
%

  U = rand(N,2);
  p = sobolset(2,'Skip',1e3,'Leap',1e2);
  S = net(p,N);

%
% Now let's approximate pi by comparing the area of a circle inscribed in a square.  The square has area 1 and is centered at (0.5, 0.5).
%

% Compute distance from (0.5, 0.5)

  r_U = sqrt((U(:,1)-0.5).^2+(U(:,2)-0.5).^2); 
  r_S = sqrt((S(:,1)-0.5).^2+(S(:,2)-0.5).^2);

% Count fraction of points within 1/2 unit of (0.5, 0.5)

  frac_U = sum(r_U <= 0.5) / N; 
  frac_S = sum(r_S <= 0.5) / N;

% Since square has side 1, circle has radius (1/2) and should have area of pi*(1/2)^2
% frac is approximately pi/4 so pi is approximately 4*frac

  mypi_U = 4*frac_U
  mypi_S = 4*frac_S

  error_U = abs(mypi_U - pi)
  error_S = abs(mypi_S - pi)