// plots lyapunov space for the logistic equation using differential calculus

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

// plot parameters
#define X_START -4.6
#define X_END 0
#define Y_START -4.6
#define Y_END -0.65

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

// restriction parameters (max and min values of iteration)
#define CHAOTIC_LIMIT 1.65
#define PERIODIC_LIMIT 1.43

// 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	200

// now for some functions we require

double getLyapunov (double x, double y, double delta_limit, int iterations) {
	// runs an iterative series to get the lyapunov exponent for
	// the logistic equation with r alternating between x and y.

	// doesn't employ a linked list, so we need number of iterations
	// future improvement: use a LL

	// first, run the logistic equation the number of iteration times
	int counter;
	double storage[iterations+1], r;
	double sum = 0;
	double intermediate;

	// initial values
	storage[0] = 0.5;

	for (counter = 1; counter < iterations; counter++) {
		// alternate r
		if ( (counter & 1) == 0 ) 
			r = x; 
		else	r = y; 
		
		storage[counter] = r * storage[counter-1] * (1-storage[counter-1]);
		
		// if it is infinite, there's nothing we can do -- abort.
		switch (isinf(storage[counter])) {
			default: break;
			case 1: return(0); break;
			case -1: return(0); break;
		}
	}

	// calculate the exponent
	for (counter = 1; counter < iterations; counter++) {
		if ( (counter & 1) == 0) r = x; 
		else r = y;

		// first calculate the inner part (that which we are to find
		// the logarithm of)
		// if it is zero, get out since log(0) = -inf
		intermediate = r-2*r*storage[counter];
		if (intermediate == 0) return(-PERIODIC_LIMIT);

		sum += log(fabs(intermediate));
	}
	// now divide
	sum = (sum / log(2)) / iterations;

	if (sum < -PERIODIC_LIMIT) return(-PERIODIC_LIMIT);
	if (sum > CHAOTIC_LIMIT) return(0);

	return(sum);

}

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 (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, and positive/negative for
	// hue (red for positive, green for negative).

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

			// and dump!

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