function hexagon_chaos ( n )

%*****************************************************************************80
%
%% hexagon_chaos() draws a fractal by applying randomly chosen maps.
%
%  Discussion:
%
%    An initial point p inside the unit hexagon is randomly chosen.
%    On each iteration, a side of the hexagon is chosen at random,
%    and p plus this side forms a triangle.  The next p is the barycenter
%    of this triangle - which is simply the average of p and the two
%    hexagon vertices.
%
%    The set of points p will gradually display a fractal of hexagons.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 April 2026
%
%  Author:
%
%    John Burkardt.
%
%  Reference:
%
%    John D Cook,
%    Random Hexagon Fractal,
%    https://www.johndcook.com/blog/2026/04/09/random-hexagon-fractal/
%    Posted 09 April 2026.
%
%  Input:
%
%    integer n: the number of times the mapping is to be applied.
%

%
%  Define vertices of unit hexagon.
%
  v = hexagon01_vertices ( );
%
%  Pick a random starting point in the hexagon.
%
  p = hexagon01_sample ( );
%
%  Iterate the map n times.
%
  clf ( );

  hold ( 'on' );
%
%  Hexagon edges.
%
    i = 6;
    for ip1 = 1 : 6
      plot ( [v(i,1),v(ip1,1)], [v(i,2),v(ip1,2)], 'm-', 'linewidth', 2 ); 
      i = ip1;
    end
%
%  Hexagon vertices.
%
    for i = 1 : 6
      plot ( v(i,1), v(i,2), 'b.', 'markersize', 15 ); 
    end
%
%  Chaos points.
%
    for i = 1 : n
%
%  Randomly pick the hexagon side with vertices v(j) and v(j+1).
%
      j = randi ( [ 1, 6 ] );
      k = mod ( j, 6 ) + 1;
%
%  Average the current point and the two vertices.
%  We are computing the barycenter of the triangle we have formed.
%
      p = ( p + v(j,1:2) + v(k,1:2) ) / 3.0;
%
%  Display the new point.
%
      plot ( p(1,1), p(1,2), 'ro', 'markersize', 2 );

    end

  hold ( 'off' );

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

  return
end
function p = hexagon01_sample ( )

%*****************************************************************************80
%
%% hexagon01_sample() samples uniformly from the regular unit hexagon.
%
%  Discussion:
%
%    The unit hexagon has center (0,0) and "radius" 1.
%
%    Instead of six triangles, we think of the hexagon as three rhombuses.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    09 April 2022
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    real P(1,2), the point.
%

%
%  Set basis vectors for each rhombus.
%
  v = [ ...
    -1.0, 0.0; ...
     0.5, sqrt(3.0)/2.0; ...
     0.5, -sqrt(3.0)/2.0; ...
    -1.0, 0.0 ];
%
%  Assign the point randomly to one of the rhombuses.
%
  rh = randi ( [ 1, 3 ], 1, 1 );
%
%  Set barycentric coordinates of the point in its rhombus.
%
  r = rand ( 1, 2 );

  if ( rh == 1 )
    p(1,1:2) = r(1,1:2) * v(1:2,1:2);
  elseif ( rh == 2 )
    p(1,1:2) = r(1,1:2) * v(2:3,1:2);
  else
    p(1,1:2) = r(1,1:2) * v(3:4,1:2);
  end

  return
end
function v = hexagon01_vertices ( )

%*****************************************************************************80
%
%% hexagon01_vertices() returns the vertices of the unit hexagon.
%
%  Discussion:
%
%    The unit hexagon has maximum radius 1, and is the hull of the points
%
%      (   1,              0 ),
%      (   0.5,   sqrt (3)/2 ),
%      ( - 0.5,   sqrt (3)/2 ),
%      ( - 1,              0 ),
%      ( - 0.5, - sqrt (3)/2 ),
%      (   0.5, - sqrt (3)/2 ).
%
%      120_____60
%        /     \
%    180/       \0
%       \       /
%        \_____/
%      240     300
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 April 2026
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    real V(2,6), the coordinates of the vertices.
%
  a = sqrt ( 3.0 ) / 2.0;

  v = [ ...
    1.0, 0.0; ...
    0.5,   a; ...
   -0.5,   a; ...
   -1.0, 0.0; ... 
   -0.5,  -a; ...
    0.5,  -a ];
 
  return
end

