function sierpinski_carpet_chaos ( n )

%*****************************************************************************80
%
%% sierpinski_carpet_chaos() draws a Sierpinski carpet using an iterated function system.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    20 August 2025
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    Scott Bailey, Theodore Kim, Robert Strichartz,
%    Inside the Levy dragon,
%    American Mathematical Monthly,
%    Volume 109, Number 8, October 2002, pages 689-703.
%
%    Michael Barnsley, Alan Sloan,
%    A Better Way to Compress Images,
%    Byte Magazine,
%    Volume 13, Number 1, January 1988, pages 215-224.
%
%    Michael Barnsley,
%    Fractals Everywhere,
%    Academic Press, 1988,
%    ISBN: 0120790696,
%    LC: QA614.86.B37.
%
%    Michael Barnsley, Lyman Hurd,
%    Fractal Image Compression,
%    Peters, 1993,
%    ISBN: 1568810008,
%    LC: TA1632.B353
%
%    Alexander Dewdney,
%    Mathematical Recreations,
%    Scientific American,
%    Volume 262, Number 5, May 1990, pages 126-129.
%
%    Bernt Wahl, Peter VanRoy, Michael Larsen, Eric Kampman,
%    Exploring Fractals on the Mac,
%    Addison Wesley, 1995,
%    ISBN: 0201626306,
%    LC: QA614.86.W34.
%
%  Input:
%
%    integer n: the number of iterations.
%

%
%  Define the linear map.
%
  A = [ ...
      0.333,  0.000;
      0.000,  0.333 ];
%
%  Define the translations.
%
  b0 = [ 0.000;  0.666 ];
  b1 = [ 0.333;  0.666 ];
  b2 = [ 0.666;  0.666 ];
  b3 = [ 0.000;  0.333 ];
  b4 = [ 0.666;  0.333 ];
  b5 = [ 0.000;  0.000 ];
  b6 = [ 0.333;  0.000 ];
  b7 = [ 0.666;  0.000 ];
%
%  Random starting point in the unit square.
%
  x = rand ( 2, 1 );
%
%  Iterate the map n times.
%
  clf ( );

  hold ( 'on' );

    for i = 1 : n
      j = randi ( [ 0, 7 ] );
      if ( j == 0 )
        x = A * x + b0;
      elseif ( j == 1 )
        x = A * x + b1;
      elseif ( j == 2 )
        x = A * x + b2;
      elseif ( j == 3 )
        x = A * x + b3;
      elseif ( j == 4 )
        x = A * x + b4;
      elseif ( j == 5 )
        x = A * x + b5;
      elseif ( j == 6 )
        x = A * x + b6;
      elseif ( j == 7 )
        x = A * x + b7;
      end
      plot ( x(1), x(2), 'bo', 'markersize', 1 );
    end

  hold ( 'off' );

  axis ( 'equal' );
  filename = 'sierpinski_carpet_chaos.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"', filename );
  close ( );

  return
end

