function menger_sponge_chaos ( n )

%*****************************************************************************80
%
%% menger_sponge_chaos() draws a Menger sponge 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.
%

%
%  To make a scatter3 plot, I create all the data first, then
%  make a single call.
%
%  To carry out the affine map, it is easiest to deal with
%  iterates that are column vectors.  This means, sadly, that
%  the b matrix of translations is a 3 x 20 array (ugh).
%
%  Define the linear map.
%
  A = [ ...
      0.333, 0.000, 0.000; ...
      0.000, 0.333, 0.000;
      0.000, 0.000, 0.333 ];
%
%  Define the translations as a 20 x 3 matrix and transpose it.
%
  b = [ ...
     0.000,  0.000, 0.000; ...
     0.000,  0.000, 0.333; ...
     0.000,  0.000, 0.666; ...
     0.000,  0.333, 0.000; ...
     0.000,  0.333, 0.666; ...
     0.000,  0.666, 0.000; ...
     0.000,  0.666, 0.333; ...
     0.000,  0.666, 0.666; ...
     0.333,  0.000, 0.000; ...
     0.333,  0.000, 0.666; ...
     0.333,  0.666, 0.000; ...
     0.333,  0.666, 0.666; ...
     0.666,  0.000, 0.000; ...
     0.666,  0.000, 0.333; ...
     0.666,  0.000, 0.666; ...
     0.666,  0.333, 0.000; ...
     0.666,  0.333, 0.666; ...
     0.666,  0.666, 0.000; ...
     0.666,  0.666, 0.333; ...
     0.666,  0.666, 0.666 ]';
%
%  Random starting point in the unit square.
%
  x = zeros ( 3, n + 1 );
  x(1:3,1) = rand ( 3, 1 );
%
%  Iterate the map n times.
%
  for i = 1 : n
    j = randi ( [ 1, 20 ] );
    x(1:3,i+1) = A * x(1:3,i) + b(1:3,j);
  end
%
%  Plot the points.
%
  clf ( );
  markersize = 3.0;
  scatter3 ( x(1,:), x(2,:), x(3,:), markersize, 'b', 'filled' )
  axis ( 'equal' );
  azimuth = 90.0;
  elevation = 0.0;
  view ( azimuth, elevation );
  title ( 'Menger 3D sponge (essentially a 2D view)' );
  filename = 'menger_sponge_chaos.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"', filename );
  close ( );

  return
end

