function mandelbrot_area ( n, it_max )

%*****************************************************************************80
%
%% mandelbrot_area() estimates the area of the Mandelbrot set.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    03 November 2025
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    John D Cook,
%    Mandelbrot area and escape times.
%    https://www.johndcook.com/blog/2025/09/08/mandelbrot-area-and-escape-times/
%    Posted: 08 September 2025
%
%  Input:
%
%    integer N: the program will check N points.
%
%    integer IT_MAX: the maximum number of iterations.
%
%
%  Define the range.
%
  x_min = -2.00;
  x_max =  2.00;
  y_min = -2.00;
  y_max =  2.00;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'mandelbrot_area():\n' );
  fprintf ( 1, '  Estimate the area of the Mandelbrot set.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For %d random values C = X + i*Y\n', n );
  fprintf ( 1, '  in the complex circle of radius 2,');
  fprintf ( 1, '  carry out up to %d iterations of the map\n', it_max );
  fprintf ( 1, '  Z(n+1) = Z(n)^2 + C.\n' );
  fprintf ( 1, '  If the iterates stay in the circle of radius 2,\n' );
  fprintf ( 1, '  then C is a member of the Mandelbrot set.\n' );
  fprintf ( 1, '\n' );
%
%  Create an array of complex sample points in [x_min,x_max] + [y_min,y_max]*i.
%
  I = ( 1 : n );
  J = ( 1 : n );
  X = ( ( I - 1 ) * x_max + ( n - I ) * x_min ) / ( n - 1 );
  Y = ( ( J - 1 ) * y_max + ( n - J ) * y_min ) / ( n - 1 );
  [ Zr, Zi ] = meshgrid ( X, Y );
  C = complex ( Zr, Zi );
%
%  Carry out the iteration.
%
  d = zeros ( n * n, 1 );

  for i = 0 : it_max
    if ( i == 0 ) 
      Z = C;
      neps = numel ( C );
      ieps = 1 : numel ( C );
      d(ieps) = it_max + 1;
    else
      Z(ieps) = Z(ieps) .* Z(ieps) + C(ieps);
      U(ieps) = abs ( Z(ieps) );
      neps = ieps ( find ( 2.0 < U(ieps) ) );
      d(neps) = i;
      ieps = ieps ( find ( U(ieps) <= 2.0 ) );
    end
    fprintf ( 1, '%d  %d\n', i, numel ( ieps ) );
  end

  fprintf ( 1, '\n')
  fprintf ( 1, '  Estimated area = %g\n', 16.0 * numel ( ieps ) / n^2);
  fprintf ( 1, '  Cook estimate was 1.5129\n' );

  clf ( );
%
%  Display the data.
%
  t = delaunay ( Zr, Zi );
%
%  Set a nonsmooth color map.
%  colorcube, flag, lines, prism.
%
  colormap ( 'prism');
%
%  Make a color contour plot.
%
  h = trisurf ( t, Zr, Zi, d, 'FaceColor', 'interp', 'EdgeColor', 'interp' );

  view ( 2 );
  axis ( 'equal' );
  axis ( 'tight' );
  xtitle = sprintf ( '%g <---X---> %g', x_min, x_max );
  xlabel ( xtitle );
  ytitle = sprintf ( '%g <---Y---> %g', y_min, y_max );
  ylabel ( ytitle );
  title_string = sprintf ( 'Mandelbrot set, %d x %d pixels, %d iterations', ...
    n, n, it_max );
  title ( title_string );
  set ( gca, 'xticklabel', [] );
  set ( gca, 'yticklabel', [] );
  filename = 'mandelbrot_area.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end

