% Monte Carlo estimation of the area of a circle using N points, with plot.

function a = MCCircle1(N)
  x = [0.6765    0.1694    0.8962   -0.8779    0.1693   -0.4298    0.6555 ...
   -0.6180   -0.1149   -0.2132    0.6531    0.3537   -0.5848   -0.3638 ...
   -0.7324    0.3429    0.1420   -0.6605   -0.7047   -0.0478    0.8162 ...
    0.1044   -0.9341   -0.8923    0.6101   -0.0973   -0.2347    0.5793 ...
   -0.2714    0.0647    0.4233    0.7430   -0.3426    0.3002    0.9497 ...
   -0.8481    0.1740   -0.1722   -0.3817   -0.4723    0.5175    0.9904 ...
   -0.6269    0.5623   -0.6084    0.9847    0.6045   -0.1515    0.4577 ...
   -0.0033    0.6180   -0.2870   -0.8535    0.1820    0.8204   -0.6125 ...
   -0.1353    0.4983   -0.9216    0.8926    0.5273    0.1176   -0.6323 ...
   -0.0041    0.0357    0.9885    0.7097    0.9248    0.3579   -0.1930 ...
    0.8700   -0.0410   -0.5364   -0.2074    0.4102    0.1171    0.5133 ...
    0.9910    0.9249    0.0701    0.9277   -0.7687   -0.8971   -0.3913 ...
    0.1604    0.0619    0.8024    0.0811   -0.1360    0.0853    0.4248 ...
   -0.9667    0.6018   -0.7150   -0.0431   -0.4863   -0.2618    0.3235 ...
   -0.6608   -0.4424];
  y =  [-0.6036   -0.6099   -0.3463    0.7607   -0.0578   -0.1921   -0.6415 ...
    0.9378   -0.1851    0.6890    0.2307   -0.2468    0.7544    0.5697 ...
   -0.0701    0.6280    0.7969   -0.1415   -0.3313    0.1933    0.8040 ...
    0.4041   -0.2451    0.4699    0.9082    0.0856    0.0802   -0.3778 ...
   -0.8575   -0.6360   -0.8140   -0.0730   -0.9813    0.8301    0.2855 ...
   -0.9972   -0.9392   -0.5831   -0.0901   -0.7455   -0.9827    0.4542 ...
   -0.2918    0.5609   -0.1267   -0.1269   -0.9016   -0.9007   -0.8178 ...
    0.1881   -0.5178    0.6827    0.7144    0.9272   -0.0222   -0.5594 ...
   -0.5476    0.0736    0.5242   -0.3049   -0.0775    0.2786    0.8347 ...
   -0.6769    0.4313    0.1555   -0.1334    0.7685   -0.2139   -0.6420 ...
    0.2667    0.2480   -0.3441    0.6059    0.9990    0.9620   -0.7459 ...
   -0.5355   -0.9527    0.2149   -0.7784   -0.1851    0.7682    0.0963 ...
   -0.2620   -0.5833   -0.1181    0.9124   -0.7519   -0.0585    0.7138 ...
   -0.9132    0.3833    0.9580   -0.4335   -0.7324    0.3706    0.8189 ...
    0.2217    0.8000];
  d2 = x.^2 + y.^2;
  inside = d2 < 1;
  count = sum(inside)
  a = 4*count/N
  axis equal;
  axis([-1.2,1.2,-1.2,1.2]);
  plot([-1,1,1,-1,-1],[-1,-1,1,1,-1],'k-');
  plot(sin(2*pi*[0:100]/100),cos(2*pi*(0:100)/100),'k-');
  plot(x,y,'k*');
end

