function mandelbrot_matlab (  )

%*****************************************************************************80
%
%% MANDELBROT_MATLAB computes an image of the Mandelbrot set.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    15 May 2017
%
%  Author:
%
%    John Burkardt
%
%  Local parameters:
%
%    Local, integer XNUM, the number of pixels in the X direction.
%
%    Local, integer YNUM, the number of pixels in the Y direction.
%
%    Local, integer COUNT_MAX, the number of iterations.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MANDELBROT_MATLAB\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Compute the Mandelbrot set and save it in a graphics file.\n' );

  xmin = - 1.00;
  xmax = - 0.60; 
  ymin =   0.00;
  ymax =   0.40;
  xnum = 1000;
  ynum = 1000;
  count_max = 200;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  For each point C = X + i*Y\n' );
  fprintf ( 1, '  with X range [%f,%f]\n', xmin, xmax );
  fprintf ( 1, '  and  Y range [%f,%f]\n', ymin, ymax );
  fprintf ( 1, '  carry out %d iterations of the map\n', count_max );
  fprintf ( 1, '  Z(n+1) = Z(n)^2 + C.\n' );
  fprintf ( 1, '  If the iterates stay bounded\n' );
  fprintf ( 1, '  then C is a member of the Mandelbrot set.\n' );

  mandelbrot_image ( xmin, xmax, ymin, ymax, xnum, ynum, count_max );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MANDELBROT_MATLAB\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function mandelbrot_image ( xmin, xmax, ymin, ymax, xnum, ynum, count_max )

%*****************************************************************************80
%
%% MANDELBROT_IMAGE creates an image of the Mandelbrot set.
%
%  Discussion:
%
%    Over the rectangle [XMIN,XMAX] x [YMIN,YMAX], determine the Mandelbrot
%    count for a grid of XNUMxYNUM points, using a particular value of
%    COUNT_MAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    15 May 2017
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XMIN, XMAX, YMIN, YMAX, the physical limits of the rectangle.
%
%    Input, integer XNUM, YNUM, the number of points in X and Y directions.
%
%    Input, integer COUNT_MAX, the maximum number of iterations.
%

%
%  Start the timer;
%
  tic;

  count = mandelbrot_set ( xmin, xmax, ymin, ymax, xnum, ynum, count_max );

  wall_time = toc;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Calculation took %g seconds.\n', wall_time );
%
%  Save the data in a graphics file.
%
  file_name = 'mandelbrot_matlab.pgm';
  comment = 'PGM gray map created by mandelbrot_matlab.m';

  pgma_write ( file_name, comment, xnum, ynum, count_max, count );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved in graphics file "%s"\n', file_name );

  return
end
function count = mandelbrot_set ( xmin, xmax, ymin, ymax, xnum, ynum, count_max )

%*****************************************************************************80
%
%% MANDELBROT_SET computes the Mandelbrot count for a grid of points.
%
%  Discussion:
%
%    Over the rectangle [XMIN,XMAX] x [YMIN,YMAX], determine the Mandelbrot
%    count for a grid of XNUMxYNUM points, using a particular value of
%    COUNT_MAX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    15 May 2017
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, float XMIN, XMAX, YMIN, YMAX, the physical limits of the rectangle.
%
%    Input, int XNUM, YNUM, the number of points in X and Y directions.
%
%    Input, int COUNT_MAX, the maximum number of iterations.
%
%    Output, int COUNT(YNUM,XNUM), the count for each point in the grid.
%
  count = zeros ( ynum, xnum );

  for i = 1 : ynum

    y = ( ( ynum - i     ) * ymax   ...
        + (        i - 1 ) * ymin ) ...
        / ( ynum     - 1 );

    for j = 1 : xnum
  
      x = ( ( xnum - j     ) * xmin   ...
          + (        j - 1 ) * xmax ) ...
          / ( xnum     - 1 );

      c = complex ( x, y );

      count(i,j) = mandelbrot_count ( c, count_max );

    end
  end

  return
end
function value = mandelbrot_count ( c, count_max )

%*****************************************************************************80
%
%% MANDELBROT_COUNT returns the Mandelbrot count for a single point.
%
%  Discussion:
%
%    Starting with the value 0, repeatedly square and add complex value C.
%
%    Return number of applications of this process at which the
%    value exceeds 2 in norm.
%
%    Repeat no more than COUNT_MAX times.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    15 May 2017
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex C, the value added at each step.
%
%    Input, integer COUNT_MAX, the maximum number of iterations.
%
%    Output, integer VALUE, the number of operations at which
%    the iterate first exceeded 2 in norm.  If this never happens,
%    return COUNT_MAX.
%
  value = count_max;

  z = complex ( 0.0, 0.0 );

  for k = 0 : count_max - 1

    if ( 2.0 <= abs ( z ) )
      value = k;
      break;
    end

    z = z * z + c;

  end

  return
end
function pgma_write ( file_name, comment, width, height, maxval, gray )

%*****************************************************************************80
%
%% PGMA_WRITE writes gray scale data into an ASCII PGM file.
%
%  The PGMA file will have the format:
%
%    P2
%    # file_name created by Matlab PGMA_WRITE.
%    width height
%    maxval, the maximum legal value for the data.
%    row 1
%    row 2
%    ...
%    row height
%
%  Example:
%
%    P2
%    # feep.pgm
%    24 7
%    15
%    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
%    0  3  3  3  3  0  0  7  7  7  7  0  0 11 11 11 11  0  0 15 15 15 15  0
%    0  3  0  0  0  0  0  7  0  0  0  0  0 11  0  0  0  0  0 15  0  0 15  0
%    0  3  3  3  0  0  0  7  7  7  0  0  0 11 11 11  0  0  0 15 15 15 15  0
%    0  3  0  0  0  0  0  7  0  0  0  0  0 11  0  0  0  0  0 15  0  0  0  0
%    0  3  0  0  0  0  0  7  7  7  7  0  0 11 11 11 11  0  0 15  0  0  0  0
%    0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 May 2017
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the name of the ASCII PGM file to be created.
%
%    Input, string COMMENT, a comment.
%
%    Input, integer WIDTH, HEIGHT, the number of columns and rows of data.
%
%    Input, integer MAXVAL, the maximum value allowed.
%
%    Input, integer GRAY(HEIGHT,WIDTH) is the gray scale data.
%
  fid = fopen ( file_name, 'wt' );

  if ( fid < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PGMA_WRITE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'PGMA_WRITE - Fatal error!' );
    return;
  end
%
%  First line identifies the file as an ASCII PGM file.
%
  fprintf ( fid, 'P2\n' );
%
%  Second line is a comment with the filename and creator.
%
  fprintf ( fid, '#%s\n', comment );
%
%  Third line holds HEIGHT and WIDTH.
%
  fprintf ( fid, '%d %d\n', height, width );
%
%  Fourth line contains the maximum gray value.
%
  fprintf ( fid, '%d\n', maxval );
%
%  Subsequent lines contain the gray values, no more than 10 to a line.
%
  for i = 1 : height
    for j = 1 : width
      fprintf ( fid, '%d ', gray(i,j) );
      if ( j == width || rem ( j, 10 ) == 0 )
        fprintf ( fid, '\n' );
      end
    end
  end

  fclose ( fid );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
