// plots lyapunov space for henon attractor. more technically accurate 
// algorithm, but slower.

#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

/*#define X_START 0.5
#define X_END 2
#define Y_START 0.2
#define Y_END 0.5*/

// output parameters
#define XSIZE 640
#define YSIZE 480

// restriction parameters (max and min values of iteration)
#define CHAOTIC_LIMIT /*4.5*/ 1
#define PERIODIC_LIMIT /*6.1*/ 0.90

// as in HSV system. Saturation is constant, and valence varies with lyap.
// exponent. All are out of 255 (max value of a byte)

#define CHAOTIC_HUE	30
#define PERIODIC_HUE	240

#define SATURATION	250

#define ITERATIONS	800
#define P_ITERATIONS	20
#define WARMUP		200

#define FIXED_DELTA 0.000000000000001	// d0

// 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, int iteration) {
	double oldx = x, oldy = y;
	//double r;

	
	x = 1 - (a*oldx*oldx) + oldy;			// HENON
	//x = 1 - a*fabs(oldx) + oldy;			// LOZI
	y = b * oldx;
}

double readjust_orbit(double bias, double old_coord, double old_delta, 
		double new_delta) {
	// returns new coordinate after having changed its distance from
	// center (center is set so that bias is at origin)

	// yeah I know, this could be done better using polar coordinates,
	// but whatever..

	return (bias + (old_delta/new_delta) * (old_coord - bias));
}

double getLyapunov (double a, double b, double delta_limit, int iterations) {
	
	double x = 0.5, y = 0.5;
	double xt, yt, sum, radiusA, average;
	int counter;

	// Iterate until orbit is on attractor
	for (counter = 0; counter < WARMUP; counter++) {
		calchenon(x, y, a, b, counter);
	}

	// Select a nearby point for the other orbit
	xt = x + FIXED_DELTA;
	yt = y + FIXED_DELTA;

	// loop for average
	for (counter = 0; counter < ITERATIONS; counter++) {
		
		// advance both points one iteration
		//for (int sec = 0; sec < 2; sec++) {
		calchenon(x, y, a, b, counter);
		calchenon(xt, yt, a, b, counter);
		//}

		// calculate new radius
		radiusA = radius(x, y) - radius(xt, yt);

		// add to the sum
		sum = sum + log(fabs(radiusA/FIXED_DELTA));

		if (isnan(sum)) return(0);

		// Readjust orbit
		xt = readjust_orbit(x, xt, FIXED_DELTA, radiusA);
		yt = readjust_orbit(y, yt, FIXED_DELTA, radiusA);

	}

	average = sum / counter;

	// band limiting
	if (average > CHAOTIC_LIMIT) return(CHAOTIC_LIMIT);
	if (average < -PERIODIC_LIMIT) return(-PERIODIC_LIMIT);

	return(average);
		
}

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, (double)1/255, 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");

	// see the other lyapunov programs for what this does

	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[0] = convert(0, max_res, 
							res, 0, 255);
				}
			}

			// and dump!

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