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

# include "leaf_chaos.h"

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

void leaf_chaos ( int n )

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

    leaf_chaos() draws a leaf using an iterated function system.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 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 A0[2*2] = { 
    0.800,  0.000,
    0.000,  0.800 };
  double A1[2*2] = { 
    0.500,  0.000,
    0.000,  0.500 };
  double A2[2*2] = { 
    0.355,  0.355,
   -0.355,  0.355 };
  double A3[2*2] = { 
    0.355, -0.355,
    0.355,  0.355 };
  double b0[2] = { 0.100, 0.040 };
  double b1[2] = { 0.250, 0.400 };
  double b2[2] = { 0.266, 0.078 };
  double b3[2] = { 0.378, 0.434 };
  FILE *command;
  char command_filename[] = "leaf_chaos_commands.txt";
  FILE *data;
  char data_filename[] = "leaf_chaos_data.txt";
  int i;
  int j;
  char plot_filename[] = "leaf_chaos.png";
  double *x;
/*
  Allocate space.
*/
  x = ( double * ) malloc ( 2 * ( n + 1 ) * sizeof ( double ) );
/*
  Initialize random number generator.
*/
  srand48 ( time ( NULL ) );
/*
  Random starting point in the unit square.
*/
  x[0+2*0] = drand48 ( );
  x[1+2*0] = drand48 ( );
/*
  Repeatedly compute and store A*x+b
*/
  for ( i = 0; i < n; i++ )
  {
    j = lrand48 ( ) % 4;
    if ( j == 0 )
    {
      r8mat_mv ( 2, 2, A0, x+(2*i), x+(2*(i+1)) );
      x[0+2*(i+1)] = x[0+2*(i+1)] + b0[0];
      x[1+2*(i+1)] = x[1+2*(i+1)] + b0[1];
    }
    else if ( j == 1 )
    {
      r8mat_mv ( 2, 2, A1, x+(2*i), x+(2*(i+1)) );
      x[0+2*(i+1)] = x[0+2*(i+1)] + b1[0];
      x[1+2*(i+1)] = x[1+2*(i+1)] + b1[1];
    }
    else if ( j == 2 )
    {
      r8mat_mv ( 2, 2, A2, x+(2*i), x+(2*(i+1)) );
      x[0+2*(i+1)] = x[0+2*(i+1)] + b2[0];
      x[1+2*(i+1)] = x[1+2*(i+1)] + b2[1];
    }
    else if ( j == 3 )
    {
      r8mat_mv ( 2, 2, A3, x+(2*i), x+(2*(i+1)) );
      x[0+2*(i+1)] = x[0+2*(i+1)] + b3[0];
      x[1+2*(i+1)] = x[1+2*(i+1)] + b3[1];
    }
  }
/*
  Create the data file.
*/
  data = fopen ( data_filename, "wt" );
  for ( i = 1; i <= n; i++ )
  {  
    fprintf ( data, "  %14.6g  %14.6g\n", x[0+2*i], x[1+2*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 output '%s'\n", plot_filename );
  fprintf ( command, "set title 'leaf_chaos'\n" );
  fprintf ( command, "set grid\n" );
  fprintf ( command, "set size square\n" );
  fprintf ( command, "unset key\n" );
  fprintf ( command, "set style data lines\n" );
  fprintf ( command, "plot '%s' using 1:2 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;
}