// Program to calculate the box counting dimension space of any attractor
// that can be implemented into a function for various values of its constants

// VERY slow. Also probably lots of kludgey code, especially the file writing
// part.

#include <iostream>
#include <fstream>
#include <math.h>
#include "plotter.h"

#define WARMUP 200	// to get minimum and maximum x/y values
#define THROWAWAY 100	// number of initial iterations we do not use because
			// the system has not been stabilized

#define ESCAPE	1000000000 // if above this, we quit and assume the particular
			   // attractor is unbounded (approximates infinity)
#define X_START -0.7
#define X_END	2.2
#define Y_START -1.3
#define Y_END	1.3

#define ITERATIONS 2200		// iterations per constant combination

// Size of box counting dimension calculation array
#define XSIZE	650
#define YSIZE	450

// Size of final output image
#define OUTXSIZE 400
#define OUTYSIZE 400

// special error values. They will be interpreted before written to the output
// file.
#define ERROR_INFINITE -1
#define ERROR_NOT_A_NUMBER -2

void calculate(double & x, double & y, double a, double b) {
	double xnew, ynew;

	// put your attractor here. This is the henon attractor.

	xnew = 1 - a*(x*x) + y;
	ynew = b * x;

	x = xnew;
	y = ynew;
}


double getDimension(double a, double b, double * output, int unique_number) {

	// gets dimension of a particular attractor
	// unique_number must be unique if you use the same output twice.
	// it is used for speed gain in not having to reinit the output array
	// every time - instead we simply count the number of pixels we set
	// to the value unique_number that has not already got that value.

	double xmin = 100, ymin = 100;
	double xmax = -100, ymax = -100;
	int counter;
	double x = 0.5, y = 0.5;
	double result;

	// for on-the-fly dimension calculation
	//int totalPixels = XSIZE*YSIZE;
	int covered_pixels = 0;
	int linear;

	Plotter tool;

	// run warmup to get minimum and maximum x/y values
	for (counter = 0; counter < THROWAWAY; counter++)
		calculate(x, y, a, b);
	
	for (counter = 0; counter < WARMUP; counter++) {
		calculate(x, y, a, b);

		if (isinf(x) || isinf(y)) return(ERROR_INFINITE);
		if (isnan(x) || isnan(y)) return(ERROR_NOT_A_NUMBER);
		if ((fabs(x) > ESCAPE) || (fabs(y) > ESCAPE))
			return(ERROR_INFINITE); // unbounded attractor
		
		if (x < xmin) xmin = x;
		if (y < ymin) ymin = y;
		if (x > xmax) xmax = x;
		if (y > ymax) ymax = y;
	}

	// set plotter data
	tool.setSourceBoundaries(xmin, ymin, xmax, ymax);
	tool.setTargetBoundaries(0, 0, XSIZE, YSIZE);
	tool.setTargetArray(output);

	for (counter = 0; counter < ITERATIONS; counter++) {
		calculate(x,y, a, b);
		
		if (isinf(x) || isinf(y)) return(ERROR_INFINITE);
                if (isnan(x) || isnan(y)) return(ERROR_NOT_A_NUMBER);
		if ((fabs(x) > ESCAPE) || (fabs(y) > ESCAPE)) 
			return(ERROR_INFINITE);
		
		// compute point we are to mark if not already marked
		// (and increment covered pixels if not already marked)
		
		linear = tool.getPlotPosition(x, y);
		if (linear != -1) {
			if (output[linear] != unique_number) {
				output[linear] = unique_number;
				covered_pixels++;
			}
		}
	}

	// get the side length of a square with area equal to the rectangle
	// a pixel (single value in output) corresponds to in the source
	// coordinate system, then calculate the box counting dimension.

	result = log(covered_pixels)/log(1/sqrt(tool.getPixelArea()));
	
	// if dimension < 0, something strange is going on, so abort
	if ((result < 0) || isinf(result)) return(ERROR_NOT_A_NUMBER); 
	
	return(result);
}

main() {
	double * output = new double[XSIZE*YSIZE];
	double * finalout = new double[OUTXSIZE*OUTYSIZE];
	double x, y;
	double increment_x, increment_y;
	double max = -100;
	ofstream xf("outdim.raw");
	int counter = 1;
	unsigned char outch;

	cout.width(4);
	cout.precision(4);
	cout.setf(ios::fixed);

	// clear arrays
	
	for (int ix = 0; ix < OUTXSIZE*OUTYSIZE; ix++) finalout[ix] = 0;
	for (int ix = 0; ix < XSIZE*YSIZE; ix++) output[ix] = 0;
	
	// plot values
	
	Plotter finalOutput;
	finalOutput.setSourceBoundaries(X_START, Y_START, X_END, Y_END);
	finalOutput.setTargetBoundaries(0, 0, OUTXSIZE, OUTYSIZE);
	finalOutput.setTargetArray(finalout);

	increment_x = finalOutput.getStepsX(OUTXSIZE);
	increment_y = finalOutput.getStepsY(OUTYSIZE);
	
	for (y = Y_START; y < Y_END; y = y + increment_y) {
		// display percentage cause it is really quite slow
		cout << y << " (";
		cout << ((counter-1)/((double)OUTXSIZE*OUTYSIZE)*100);
		cout << " %)\n";
		for (x = X_START; x < X_END; x = x + increment_x) {
			counter++;
			finalOutput.plot(x, y, getDimension(x, y, output, counter));
		}
	}

	// get maximum value so that we can scale it up to 255 later on
	// Using an exposure function later on would be preferrable, but
	// that hasn't been implemented.

	for (counter = 0; counter < (OUTXSIZE*OUTYSIZE); counter++) {
		if (finalout[counter] > max) max = finalout[counter];
	}

	// dump to file.

	for (counter = 0; counter < (OUTXSIZE*OUTYSIZE);counter++) {
		switch ((int)rint(finalout[counter])) {
			default:
				outch = (unsigned char)(rint((finalout[counter]
								/max)*255));
				xf << outch << outch << outch; 
				break;
			case ERROR_INFINITE:
				xf << (char)0 << (char)0 << (char)128;
				break;
			case ERROR_NOT_A_NUMBER:
				xf << (char)128 << (char)0 << (char)0;
				break;
		}
	}

	xf.close();
}