#! /usr/bin/env python3
#
def menger_sponge_chaos_test ( ):

#*****************************************************************************80
#
## menger_sponge_chaos_test() tests menger_sponge_chaos().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    19 August 2025
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'menger_sponge_chaos_test():' )
  print ( '  python version: ' + platform.python_version ( ) )
  print ( '  numpy version:  ' + np.version.version )
  print ( '  menger_sponge_chaos() uses an iterated map to plot' )
  print ( '  the 3D Menger sponge.' )
  print ( '' )
  print ( '  The image really has to be viewed interactively' )
  print ( '  in order to rotate the cube and see the holes.' )

  n = 10000
  print ( '  Apply the iteration map', n, 'times.' )
  menger_sponge_chaos ( n )
#
#  Terminate.
#
  print ( '' )
  print ( 'menger_sponge_chaos_test():' )
  print ( '  Normal end of execution.' )

  return

def menger_sponge_chaos ( n ):

#*****************************************************************************80
#
## menger_sponge_chaos() draws the Menger sponge using an iterated function system.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    19 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.
#
  from mpl_toolkits.mplot3d import Axes3D
  import matplotlib.pyplot as plt
  import numpy as np
#
#  Define the linear map.
#
  A = np.array ( [ \
    [  0.333,  0.000,  0.000 ], \
    [  0.000,  0.333,  0.000 ],
    [  0.000,  0.000,  0.333 ] ] )
#
#  Define the translations.
#
  b = np.array ( [ \
    [ 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.
#
  xplot = np.zeros ( [ n + 1, 3 ] )
  xplot[0,:] = np.random.random ( size = 3 )
#
#  Iterate the map n times.
#
  for k in range ( n ):
    i = np.random.choice ( 20 )
    xplot[k+1,:] = np.dot ( A, xplot[k,:] ) + b[i,:]
#
#  Plot the data.
#
  fig = plt.figure ( )
  ax = fig.add_subplot ( 111, projection = '3d' )
  ax.scatter ( xplot[:,0], xplot[:,1], xplot[:,2], s = 4 )
  ax.set_title ( 'Menger sponge', fontsize = 16 )
  ax.set_xlabel ( '<--- x --->' )
  ax.set_ylabel ( '<--- y --->' )
  ax.set_zlabel ( '<--- z --->' )
  filename = 'menger_sponge_chaos.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
#
#  Turn this to True to enable interactive viewing.
#
  if ( False ):
    plt.show ( )
  plt.close ( )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  menger_sponge_chaos_test ( )
  timestamp ( )