# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "menger_sponge_chaos.h"

/******************************************************************************/

void menger_sponge_chaos ( int n )

/******************************************************************************/
/*
  Purpose:

    menger_sponge_chaos() draws a Menger sponge using an iterated function system.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 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.
*/
{
  double A[3*3] = { 
    0.333, 0.000, 0.000,
    0.000, 0.333, 0.000,
    0.000, 0.000, 0.333 };
  double b[3*20] = {
     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 };
  FILE *command;
  char command_filename[] = "menger_sponge_chaos_commands.txt";
  FILE *data;
  char data_filename[] = "menger_sponge_chaos_data.txt";
  int i;
  int j;
  char plot_filename[] = "menger_sponge_chaos.png";
  double *x;
/*
  Allocate space.
*/
  x = ( double * ) malloc ( 3 * ( n + 1 ) * sizeof ( double ) );
/*
  Initialize random number generator.
*/
  srand48 ( time ( NULL ) );
/*
  Random starting point in the unit square.
*/
  x[0+3*0] = drand48 ( );
  x[1+3*0] = drand48 ( );
  x[2+3*0] = drand48 ( );
/*
  Repeatedly compute and store A*x+b
*/
  for ( i = 0; i < n; i++ )
  {
    j = lrand48 ( ) % 20;
    r8mat_mv ( 3, 3, A, x+(3*i), x+(3*(i+1)) );
    x[0+3*(i+1)] = x[0+3*(i+1)] + b[0+3*j];
    x[1+3*(i+1)] = x[1+3*(i+1)] + b[1+3*j];
    x[2+3*(i+1)] = x[2+3*(i+1)] + b[2+3*j];
  }
/*
  Create the data file.
*/
  data = fopen ( data_filename, "wt" );
  for ( i = 1; i <= n; i++ )
  {  
    fprintf ( data, "  %14.6g  %14.6g  %14.6g\n", x[0+3*i], x[1+3*i], x[2+3*i] );
  }
  fclose ( data );
  printf ( "  Created data file '%s'\n", data_filename );
/*
  Create the command file.
*/
  command = fopen ( command_filename, "wt" );

  fprintf ( command, "# %s\n", command_filename );
  fprintf ( command, "#\n" );
  fprintf ( command, "# Usage:\n" );
  fprintf ( command, "#  gnuplot < %s\n", command_filename );
  fprintf ( command, "#\n" );
  fprintf ( command, "set term png\n" );
  fprintf ( command, "set style line 1 lt 3 pt 4\n" );
  fprintf ( command, "set output '%s'\n", plot_filename );
  fprintf ( command, "set title 'menger\\_sponge\\_chaos'\n" );
  fprintf ( command, "set xlabel '<-- X -->'\n" );
  fprintf ( command, "set ylabel '<-- Y -->'\n" );
  fprintf ( command, "set zlabel '<-- Z -->'\n" );
  fprintf ( command, "set grid\n" );
  fprintf ( command, "set pointsize\n" );
  fprintf ( command, "set view equal xyz\n" );
  fprintf ( command, "unset key\n" );
  fprintf ( command, "set style data lines\n" );
  fprintf ( command, "set view 0, 0, 1, 1\n" );
  fprintf ( command, "splot '%s' with points pt 7 ps 0.5\n", data_filename );
  fprintf ( command, "quit\n" );

  fclose ( command );

  printf ( "  Created command file '%s'\n", command_filename );

  free ( x );

  return;
}
/******************************************************************************/

void r8mat_mv ( int m, int n, double a[], double x[], double ax[] )

/******************************************************************************/
/*
  Purpose:

    r8mat_mv() multiplies a matrix times a vector.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

    For this routine, the result is returned as an argument.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 April 2014

  Author:

    John Burkardt

  Input:

    int M, N, the number of rows and columns of the matrix.

    double A[M,N], the M by N matrix.

    double X[N], the vector to be multiplied by A.

  Output:

    double AX[M], the product A*X.
*/
{
  int i;
  int j;
  double *y;

  y = ( double * ) malloc ( m * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    y[i] = 0.0;
    for ( j = 0; j < n; j++ )
    {
      y[i] = y[i] + a[i+j*m] * x[j];
    }
  }

  for ( i = 0; i < m; i++ )
  {
    ax[i] = y[i];
  }

  free ( y );

  return;
}