# include <cmath>
# include <cstdlib>
# include <cstring>
# include <fstream>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "cross_chaos.hpp"

//****************************************************************************80

void cross_chaos ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    cross_chaos() draws a cross using an iterated function system.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 August 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of iterations.
//
{
  double A[2*2] = { 
    0.333, 0.000,
    0.000, 0.333 };
  double b[2*5] = {
    0.333, 0.000,
    0.000, 0.333,
    0.333, 0.333,
    0.666, 0.333,
    0.333, 0.666 };
  ofstream command;
  string command_filename = "cross_chaos_commands.txt";
  ofstream data;
  string data_filename = "cross_chaos_data.txt";
  int i;
  int j;
  string plot_filename = "cross_chaos.png";
  double *x;
//
//  Allocate space.
//
  x = new double[ 2 * ( n + 1 ) ];
//
//  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 ( ) % 5;
    r8mat_mv ( 2, 2, A, x+(2*i), x+(2*(i+1)) );
    x[0+2*(i+1)] = x[0+2*(i+1)] + b[0+2*j];
    x[1+2*(i+1)] = x[1+2*(i+1)] + b[1+2*j];
  }
//
//  Create the data file.
//
  data.open ( data_filename.c_str ( ) );
  for ( i = 1; i <= n; i++ )
  {  
    data << "  " << x[0+2*i] << "  " << x[1+2*i] << "\n";
  }
  data.close ( );
  cout << "  Created data file '" << data_filename << "'\n";
//
//  Create the command file.
//
  command.open ( command_filename.c_str ( ) );

  command << "# " << command_filename << "\n";
  command << "#\n";
  command << "# Usage:\n";
  command << "#  gnuplot < " << command_filename << "\n";
  command << "#\n";
  command << "set term png\n";
  command << "set output '" << plot_filename << "'\n";
  command << "set title 'cross\\_chaos'\n";
  command << "set grid\n";
  command << "set size square\n";
  command << "unset key\n";
  command << "set style data lines\n";
  command << "plot '" << data_filename << "' using 1:2 with points pt 7 ps 0.5\n";
  command << "quit\n";

  command.close ( );

  cout << "  Created command file '" << command_filename << "'\n";

  delete [] x;

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

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

//****************************************************************************80
//
//  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:
//
//    11 April 2007
//
//  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;

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

  return;
}