% L6Bmontecarlo
% Demonstration of simple Monte Carlo integration of an area
%
function p=L6Bmontecarlo(n, showpoints, showbounds)
  if nargin<1, n=100; end
  if nargin<2, showpoints=false; end
  if nargin<3, showbounds=false; end
  
  x = rand(1,n)*2 - 1;
  y = rand(1,n)*2 - 1;
  in = testpoint(x,y);
  
  p = sum(in)/n;
  
  if showpoints
    clf
    plot(x(in),y(in),'k.');
    axis equal
    axis([-1.2 1.2 -1.2 1.2])
    hold on
    plot(x(not(in)),y(not(in)),'r.');
  end
  
  if showbounds
    phi = linspace(0, 2*pi, 1000);
    plot(cos(phi), sin(phi), 'b-', 'linewidth', 2);
    hold on
    plot([-1 1 1 -1 -1], [-1 -1 1 1 -1], 'k-', 'linewidth', 2);
  end
end


% testpoint: return 1 if point inside figure, 0 if outside.
function in=testpoint(x,y)
  dist = sqrt(x.^2 + y.^2);
  in   = (dist < 1);
end