// plots lyapunov space for the henon attractor by a numerical difference
// method. It is not technically accurate, but produces nice images.

#include <iostream>
#include <math.h>
#include <stdio.h>

// plot parameters
#define X_START -1.2//-5.3
#define X_END 1.2//5.3
#define Y_START -1.2//-2.5
#define Y_END 1.2//2.5

// output parameters
#define XSIZE 2048
#define YSIZE 1760

// restriction parameters (max and min values of iteration)
#define CHAOTIC_LIMIT 7.5
#define PERIODIC_LIMIT 10.1

#define ITERATIONS 800

// now for some functions we require

inline double avg(double a, double b) { return ((a+b)/2); }
// alternately
inline double radius(double a, double b) {return (sqrt((a*a)+(b*b)));}

void calchenon(double & x, double & y, double a, double b) {
	// iterate the henon attractor once
	double oldx = x, oldy = y;

	x = 1 - (a*tan(oldx)) + oldy;
	y = b * oldx;
}

double getdifference(double a, double b, double & x, double & y, 
		int iterations) {
	// run the iterations for two parallel cases iterations number of
	// times, then dump it
	double xt, yt;

	xt = x - 0.001;
	yt = y - 0.001;

	for (int counter = 0; counter < iterations; counter++) {
		calchenon(x, y, a, b);
		calchenon(xt, yt, a, b);
		if (isinf(x)) return(x);
	}

	return(radius(x, y)-radius(xt, yt));
}

	

double getLyapunov (double a, double b, int iterations) {
	
	double x, y, xt, yt;
	double distance, old_distance;
	double l_e = 0;
	double sum;
	int counter;

	x = 0.5; y = 0.5;

	for (counter = 0; counter < (iterations/10); counter++) {

		// get the difference. The more apart two initially close points
		// are, the more bits we tend to lose for any given point, and
		// thus the more chaotic the attractor.

		l_e += fabs(getdifference(a, b, x, y, 10));

		switch (isinf(l_e)) {
			default: break;
			// UNBOUNDED
			case 1: return(l_e); break;
			case 2: return(-PERIODIC_LIMIT); break;
		}
	}

	sum = log(l_e/((counter-1)/10));
	
	if (sum < -PERIODIC_LIMIT) return(-PERIODIC_LIMIT);

	return(sum);

}

// the functions below are used because plotter was made after this app
// yes, it would be possible to cut them out and use plotter.h instead - that
// is left to the reader.

inline int linear(int x, int y, int xsize) { return(y*xsize+x); }
inline int convert(double min, double max, double current, double destmin, 
		double destmax) {

	// converts a number in the range [min...max] to a number in the
	// range [destmin...destmax]

	return( (int)rint(destmin+((current-min)/(max-min)*(destmax-destmin))));
}

inline double dconvert(double min, double max, double current, double destmin,
		double destmax) {
	return (destmin+((current-min)/(max-min)*(destmax-destmin)));
}
	
main() {
	// init arrays
	double * parameters = new double [ XSIZE * YSIZE ];
	int counter, sec;
	double x, y, res;
	double hsv[3];
	unsigned char rgb[3];		// a singular RGB triplet

	FILE * output;

	// the following are record numbers (storing minimum and maximum
	// values), and set to ridiculous values so that they are guaranteed
	// to be set at least once.
	double max_res = -1000, min_res = 1000;

	for (counter = 0; counter < YSIZE; counter++) {
		y = dconvert(0, YSIZE, counter, Y_START, Y_END);
		cout << counter << " " <<y << endl;
		
		for (sec = 0; sec < XSIZE; sec++) {
			x = dconvert(0, XSIZE, sec, X_START, X_END);
			
			res = getLyapunov(x, y, ITERATIONS);

			// res > 0 means the result is chaotic, otherwise
			// it is periodic.

			// Here we get the linear number corresponding to
			// the current coordinate, but do not translate
			// the scalar point (res) to color yet -- we need
			// minimum and maximum values for that.
			if (!isinf(res)) {
			if (res < min_res) min_res = res;
			else if (res > max_res) max_res = res;
			}

			parameters[linear(sec, counter, XSIZE)] = res;
		}
	}

	// open the file to which we are going to write
	// no error checking - proceed at own risk
	output = fopen("out.raw", "wb");

	// Now we want to color the chaotic regions with one hue, periodic
	// with another and vary the intensity according to the 
	// lyapunov exponents.

	// since all chaotic results are positive (see above), we use the
	// absolute value for intensity calculation, and positive/negative for
	// hue (red if chaotic, green if periodic).

	for (counter = 0; counter < YSIZE; counter++) {
		for (sec = 0; sec < XSIZE; sec++) {
			res = parameters[linear(sec, counter, XSIZE)];
			
			rgb[0] = 0;
			rgb[1] = 0;
			rgb[2] = 0;
			
			if (res < 0) {
				// periodic
				rgb[1] = convert(0, fabs(min_res), fabs(res),
						0, 255);
			} else {
				if (isinf(res)) {
					// unbounded
					rgb[0] = 0;
					rgb[1] = 0;
					rgb[2] = 32;
				} else {
					// chaotic
					rgb[2] = convert(0, max_res, 
							res, 0, 255);
				}
			}

			// and dump!

			fwrite(rgb, 3, 1, output);
		}
	}
			
	fclose(output);
}