function mypi = silly_pi(num_iters)
%
% Calculate pi by approximating the area of a circle by a series
% of smaller and smaller boxes.  Iterate the given number of times.
% Print a table showing value of pi at each iteration.
%
% Arguments:
%                 num_iters    (input)       number of times we double the
%                                                 number of boxes
%        
%                 mypi        (output)       calculated approximation to pi
%
% MWR 4/22/2002
%

debug = 0;

% sanity check
if (num_iters <= 0)
  error('num_iters must be greater than zero');
end

% prepare to enter main loop
radius = 1.0;
num_boxes = 1;
old_sum = 0;
for (iter = 1 : num_iters)
  
  box_length = (2.0*radius)/num_boxes;
  box_rad = box_length*0.5;
  box_area = box_length*box_length;

  area_sum = 0;
  for (xi = 1 : num_boxes)
    xc = -1 + box_rad + box_length*(xi - 1);
    if (debug > 0)
      fprintf(1, 'iter %5d  box %5d  xc %9.4e \n', iter, xi, xc);
    end
    
    for (yi = 1 : num_boxes)
      yc = -1 + box_rad + box_length*(yi - 1);
   
      % now, find the center of this box, and decide whether it lies
      %  inside a circle of radius 1 or not.  If it does, we add
      %  its area to the running total
      dist = sqrt((xc*xc) + (yc*yc));
      if (dist < radius)
        area_sum = area_sum + box_area;
        if (debug > 1)
          fprintf(1, '    box %9.4e %9.4e  dist %9.4e inside  sum %9.4e\n', ...
                         xc, yc, dist, area_sum);
        end
      end

    end 
  end

  % check to see how much the sum has changed in the past iteration
  fractional_change = abs((area_sum - old_sum)/area_sum);
  if (debug > 0)
    fprintf(1, '  area = %9.4e  old = %9.4e   change %9.4e \n', ...
           area_sum, old_sum, fractional_change);
  end

  % print out a nice line summarizing this iteration's result
  mypi = area_sum/(radius*radius);
  fprintf(1, '  iter %5d  boxes %8d  pi %10.6f  frac_change %9.4e \n', ...
              iter, num_boxes*num_boxes, mypi, fractional_change);

  % prepare for next iteration
  old_sum = area_sum;
  num_boxes = num_boxes*2;
  iter = iter + 1;

end


% now calculate pi
mypi = area_sum/(radius*radius);


