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


// the program now produces something almost, but not entirely unlike shadows


/* Twiddle these constants to get different results */

#define LOW_LYAP_LIMIT 0.1	// all attractors with lyap. exponent less than
				// this are set as not viable
#define LOW_FDIM_LIMIT  1.0	// same, but with correlation dimension
#define HIGH_FDIM_LIMIT 3.2	// same, but greater than this

/* Twiddle these constants to affect output */
#define XZ 1024   // x size of raw file output
#define YZ 768   // y size

#define ZZ_L 0x20 // lowest grey intensity, to not interfere with future shadows
#define ZZ 0xFF   // highest grey intensity
#define SZ 0x00	  // shadow intensity (when shadows are implemented)

/* These define the range of the random function p_PopulateConstants */
#define MAX 1.8 //2.1
#define STEPS 32 //42
#define STEPSIZE 0.1

/* Colors */
#define HUE 29
#define SAT_BRIGHT 255
#define SAT_SHADOW 200

/* Iteration settings */
#define WARMUP 70000
#define LIMIT  850000

//#define SHADOWS

// y axis goes straight up, x axis to the right, and z axis outwards
// all the angles below are in degrees

#define r_phi	3	// rotation that leaves x axis unchanged
#define r_theta	3	// rotation that leaves y axis unchanged
#define r_eta	3	// rotation that leaves z axis unchanged

#define setsin sin(r_phi), sin(r_theta), sin(r_eta)
#define setcos cos(r_phi), cos(r_theta), cos(r_eta)

inline double degToRadian(double x) {
	return (x/180)*3.1415926;
}

double ipow(double mantissa, unsigned int exponent);


//////////////////////////////////////////////////////////////////////////////

class Coordinate {
	// Q&D i know
	public:
		double x;
		double y;
		double z;
		
		Coordinate(); 
		Coordinate(double xz, double yz, double zz);
		Coordinate operator *= (double b);
		Coordinate operator * (Coordinate);
		Coordinate operator + (Coordinate);
		Coordinate operator / (Coordinate);
		Coordinate operator - (Coordinate);
		void zero();
		bool is_greater(int than);
		void adjust(bool minima, Coordinate source);
		double getRadius();
		double getBiasedRadius(Coordinate bias);
		double getAverage();
		void set(double xz, double yz, double zz);
		bool isNaN();
		bool isZero();
		
};

Coordinate::Coordinate() { x = 0; y = 0; }

bool Coordinate::is_greater(int than) {
	if ( (x > than) || (y > than) || (z > than) ) return(true);
	return (false);
}

bool Coordinate::isNaN() {
	return ( isnan(x) || isnan(y));
}

void Coordinate::adjust(bool minima, Coordinate source) {
	if (minima) {
		if (x > source.x) x = source.x;
		if (y > source.y) y = source.y;
		if (z > source.z) z = source.z;
	} else {
		if (x < source.x) x = source.x;
		if (y < source.y) y = source.y;
		if (z < source.z) z = source.z;
	}
}

double Coordinate::getRadius() {
	return(	sqrt(ipow(x, 2) + ipow(y, 2) + ipow(z, 2)) );
}

double Coordinate::getBiasedRadius(Coordinate bias) {
	// calculate radius when relative (0,0) is absolute location (bias)
	return(sqrt(ipow(x-bias.x, 2) + ipow(y-bias.y, 2) + ipow(z-bias.z, 2)));
}

double Coordinate::getAverage() {
	return ( (x+y+z)/3);
}

Coordinate::Coordinate(double xz, double yz, double zz) {
	x = xz;
	y = yz;
	z = zz;
}

void Coordinate::set(double xz, double yz, double zz) {
	x = xz;
	y = yz;
	z = zz;
}

Coordinate Coordinate::operator*= (double b) {
	x = x * b;
	y = y * b;
	z = z * b;
}

Coordinate Coordinate::operator+ (Coordinate p) {
	Coordinate temp;
	temp.x = x + p.x;
	temp.y = y + p.y;
	temp.z = z + p.z;
	return(temp);
}

Coordinate Coordinate::operator* (Coordinate p) {
	Coordinate t;
	t.x = x * p.x;
	t.y = y * p.y;
	t.z = z * p.z;
	return(t);
}

Coordinate Coordinate::operator-(Coordinate p) {
	Coordinate t;
	t.x = x - p.x;
	t.y = y - p.y;
	t.z = z - p.z;
	return(t);
}

Coordinate Coordinate::operator/(Coordinate p) {
	Coordinate t;
	t.x = x / p.x;
	t.y = y / p.y;
	t.z = z / p.z;
	return(t);
}

void Coordinate::zero() {
	x = 0;
	y = 0;
	z = 0;
}

Coordinate real_coord(Coordinate minimum, Coordinate maximum, 
		      Coordinate position, Coordinate dimensions) {

	// scale a coordinate (position) in area [minimum..maximum] to 
	// its correct representation in area [0..dimensions]

	return((position - minimum)/(maximum-minimum)*dimensions);
}

bool Coordinate::isZero() {
	return( (x == 0) && (y == 0) && (z == 0) );
}


///////////////////////////////////////////////////////////////////////////////


			// Once and Only Once
class Rotation {

	private:
		Coordinate sine, cosine;
		Coordinate center;
		int number;	// for finding center of mass on-the-fly
		void rotateOneAxis(double sine_angle, double cosine_angle, 
				double & first, double & second,
				double & third);

	public:
		Rotation(double phi, double theta, double eta);
		void setCenter(Coordinate inCenter);
		void setCenter(Coordinate * points, int N);
		Coordinate getCenter();
		void setNewRotation(double phi, double theta, double eta);
		Coordinate rotate(Coordinate point);

		// figuring the center from center of mass, on the fly
		void setNumberofPoints(unsigned int inNumber);
		bool addPoint(Coordinate point);
};

void Rotation::rotateOneAxis(double sine_angle, double cosine_angle,
		double & first, double & second, double & third) {

	// internal function to rotate about one axis.
	// alters first, second, and third.

	double tempFirst, tempSecond;

	tempFirst = (sine_angle * first) + (cosine_angle * second);
	tempSecond = (cosine_angle * first) - (sine_angle * second);

	first = tempFirst;
	second = tempSecond;
}

void Rotation::setNewRotation(double phi, double theta, double eta) {
	sine.set(sin(phi), sin(theta), sin(eta));
	cosine.set(cos(phi), cos(theta), cos(eta));
}

Rotation::Rotation(double phi, double theta, double eta) {
	setNewRotation(phi, theta, eta);
	center.set(0,0,0);
	number = 0;
}

void Rotation::setCenter(Coordinate inCenter) {
	center = inCenter;
}

void Rotation::setCenter(Coordinate * points, int N) {
	// sets center based on N points around the center
	// does this by summing up coordinates, then dividing by N

	Coordinate tempCenter(0,0,0);

	for (int counter = 0; counter < N; counter++) 
		tempCenter = tempCenter + points[counter];

	center.set(tempCenter.x/N, tempCenter.y/N, tempCenter.z/N);
}

Coordinate Rotation::getCenter() {
	Coordinate xyz = center;
	return(xyz);
}

Coordinate Rotation::rotate(Coordinate point) {
	// rotate point around axes phi, theta, eta
	// (coordinates sine and cosine refer to sine and cosine of these)
	// and spit out the rotated coordinate

	Coordinate result, working = point;

	// .x = phi, .y = theta, .z = eta

	// first normalize center
	working = working - center;

	// then rotate:


	/* rotate about theta (keeps y value intact)
	result.x = (cosine.y * working.x) - (sine.y * working.z);
	result.y = working.y;
	result.z = (sine.y * working.x) + (cosine.y * working.z);

	// rotate about phi (keeps x value intact)
	working.x = result.x;
	working.y = (cosine.x * result.y) + (sine.x * result.z);
	working.z = (sine.x * result.y) - (cosine.x * result.z);

	// rotate about eta (keeps z value intact)
	result.x = (cosine.z * working.x) - (sine.z * working.y);
	result.y = (sine.z * working.x) + (cosine.z * working.y);
	result.z = working.z;*/

	rotateOneAxis(sine.x, cosine.x, working.x, working.y, working.z);
	//rotateOneAxis(sine.z, cosine.z, working.z, working.x, working.y);
	//rotateOneAxis(sine.y, cosine.y, working.y, working.z, working.x);

	// normalize back

	return(working + center);
}

void Rotation::setNumberofPoints(unsigned int inNumber) {
	number = inNumber;
	if (inNumber != 0) center.set(0,0,0); // so addpoints will work
}

bool Rotation::addPoint(Coordinate point) {
	if (number < 1) return(false);
	Coordinate toAdd(point.x/number, point.y/number, point.z/number);
	center = center + toAdd;
	return(true);
}
				
	

///////////////////////////////////////////////////////////////////////////////

class FractalTypes {
	public:
		// organization functions
		void p_populateConstants(double * what, int number);
		void revPop(double * what, int number);
		void invRevPop(double * what, char * input, int number);
		double getLyap(Coordinate * iterates, double * a);
		
		// Strange attractors
		Coordinate henon(Coordinate old, double alpha, double beta);
		Coordinate henonGeneralized(Coordinate old, double * a);
			// polynomial attractor
};

Coordinate FractalTypes::henon(Coordinate old, double alpha, double beta) {

	Coordinate newC;

	newC.x = 1 - (alpha*old.x*old.x)+old.y;
	newC.y = beta * old.x;

	return(newC);
}

Coordinate FractalTypes::henonGeneralized(Coordinate old, double * a) {
	/*//Xn+1 = a1 + a2*Xn + a3*(Xn)^2 + a4*Xn*Yn + a5*Yn + a6*(Yn)^2
	//Yn+1 = a7 + a8*Xn + a9*(Xn)^2 + a10*Xn*Yn + a11*Yn + a12*(Yn)^2

	Coordinate newC;

	double sP_oldX = old.x * old.x;
	double sP_oldY = old.y * old.y;
	double sP_XY = old.y * old.x;

	newC.x = a[0] + a[1]*old.x + a[2]*sP_oldX + (a[3]*sP_XY) +
		(a[4]*old.y) + (a[5]*sP_oldY);
	
	newC.y = a[6] + (a[7]*old.x) + (a[8]*sP_oldX) + (a[9]*sP_XY) +
		(a[10]*old.y) + (a[11]*sP_oldY);

	return(newC);*/

	// 	S	L	O	W. 'nuff said.
	// fixed it a bit by using i_pow, but still not as fast as hardcoding
	
	Coordinate newC;
	unsigned char degree = 2;
	unsigned char counter, sec, tet, tri = 0;
	/*unsigned char y_split_point = ((degree+1)>>1)*(degree+2);
	unsigned char z_split_point = y_split_point * 2;*/
	double coeff;

	newC.zero();

	for (counter = 0; counter <= degree; counter++) {
		for (sec = 0; sec <= degree-counter; sec++) {
			for (tet = 0; tet <= sec; tet++) {
				coeff = ipow(old.x, counter) * ipow(old.y, sec) * ipow(old.z, tet);
				newC.x = newC.x + (a[tri++] * coeff);
				newC.y = newC.y + (a[tri++] * coeff); 
				newC.z = newC.z + (a[tri++] * coeff);
				//cout << (int)tri << " - ";
			}
		}
	}
	//cout << "\n";
	return(newC);
	//*/
}

double ipow(double mantissa, unsigned int exponent) {
	/* Integer POWer-of function for polynomials (positive integer 
	   exponents only)
	 ya ya, 0^0 is incorrect here.. bwah*/

	if (exponent == 0) return(1);

	double result = mantissa;

	while (--exponent > 0) result *= mantissa;

	return(result);
}

double FractalTypes::getLyap(Coordinate * iterates, double * a) {

	// gets Lyapunov exponent estimate

	// iterates: to conserve calculations (cuts time in half)

	// output: a number that's greater the more sensitive the system
	// is to noise, and thus indirectly a degree of its chaos

	// sets that give negative values are very likely to be periodic
	// e.g ellipses and such
	
	Coordinate g(0.2,0.2, 0.2);

	double old_distance = 0;
	double distance = 0;
	double l_e = 0;
	int counter;
	
	g = henonGeneralized(g, a);
	old_distance = fabs(iterates[0].getAverage() - g.getAverage());

	for (counter = 0; counter < 35; counter++) {
		g = henonGeneralized(g, a);
		distance = fabs(iterates[counter].getAverage() - 
				g.getAverage());
		if (old_distance > 0) l_e = l_e + (distance/old_distance);
		old_distance = distance;
	}

	return(log(l_e/(counter-1))/log(2));
}
		

void FractalTypes::p_populateConstants(double * what, int number) {
	// set what[0] .. what[number-1] inclusive to some random value
	// from -MAX to MAX in steps of STEPSIZE

	// this we do by generating (for each what[x]) a random number
	// 0..STEPS-1 inc, multiplying by STEPSIZE, and subtracting MAX.

	int counter;
	
	for (counter = 0; counter < number; counter++)
		what[counter] = MAX-((random() % STEPS)*STEPSIZE);

}

void FractalTypes::revPop(double * what, int number) {
	// generate ascii representation of constants

	int counter, temp;
	for (counter = 0; counter < number; counter++){
		temp = (int)(rint((MAX + STEPSIZE + what[counter])*10));
		if (temp < 10) {
			cout << (char)(temp+'0');
		} else {
			cout << (char)(temp-10+'A');
		}
	}
	cout << "\n";
}

void FractalTypes::invRevPop(double * what, char * input, int number) {
	int counter, temp;

	for (counter = 0; counter < number; counter++) {
		temp = input[counter];
		if (temp < 'A') temp = temp - '0';
		else temp = temp+10-'A';
		
		what[counter] = (temp-(10*MAX)-(10*STEPSIZE))/10;
	}
}
		
///////////////////////////////////////////////////////////////////////////////

unsigned int getPos(int x, int y, int xsize) {
	return((y*xsize)+x);
}

double getCorrelationDimension (Coordinate * iterates, int iterations,
		Coordinate att_min, Coordinate att_max) {

	// gets correlation dimension (F)
	// does this by counting number of points within two virtual
	// circles, based on iterates[0..iterations-1]

	// the circles are scaled to 6 and 0.6 % of the attractor size, as
	// given by upper left and lower right corners att_min and att_max

	int within_a = 0, within_b = 0, counter;
	double a, b, radius;

	Coordinate bias;

	// set radius of circles. b is the greater one

	b = sqrt( pow(att_max.y-att_min.y,2) + pow(att_max.x-att_min.x, 2) );
	a = b * 0.006;
	b = b * 0.06;

	// set center to a place we know will be populated with at least
	// one dot
	
	bias = iterates[iterations-1];

	// count points

	for (counter = 0; counter < iterations; counter++) {
		radius = fabs(iterates[counter].getBiasedRadius(bias));

		if (radius < b) {
			within_b++;
			if (radius < a) within_a++;
		}
	}

	// calculate F
	if (within_b == 0) return(0);	// avoid division by zero
	if (within_a == 0) return(0);	// ditto

	return(log((double)within_b/(double)within_a)/log(10));
}

void findAttractors(int space_limit, double lower_difference_limit, 
		double upper_difference_limit, 
		double lower_F_limit, double upper_F_limit,
		double lyap_lower_limit,
		double * constants, Coordinate & min, Coordinate & max,
		double & F, double & L) {

	#define ITERATIONS 200
	#define ITERATIONS_THOROUGH 450
	
	Coordinate C, CBack[ITERATIONS_THOROUGH];
	FractalTypes X;
	double diff, startpoint;
	int counter, sec;

	// avoid a warning

	min.zero();
	max.zero();

	diff = lower_difference_limit - 1; // to make it work first time

	while ( C.is_greater(space_limit) || diff < lower_difference_limit 
			|| F > upper_F_limit || F < lower_F_limit 
			|| diff > upper_difference_limit 
			|| L < lyap_lower_limit) {

		// reset parameters so the above doesn't fall through if
		// unwarranted
		F = lower_F_limit - 1;
		L = lyap_lower_limit - 1;

		// reset coordinates and difference
		C.set(0.1, 0.1, 0.1);
		min.set(100, 100, 100);
		max.set(-100, -100, -100);
		diff = 0;

		// set random constants
		X.p_populateConstants(constants, 36);
		/*X.invRevPop(constants, "EFI8PEBKLNB3H8D84N43E24FK6DM2J", 30);
		X.revPop(constants, 30);*/
		
		// "warm up" by running 15 times
		CBack[0] = C;
		for (counter = 0; counter < 15; counter++)
			CBack[0] = X.henonGeneralized(CBack[0], constants);

		// then start computing. get attractor size while we're at it
		for (counter = 1; counter < ITERATIONS; counter++) {
			CBack[counter] = X.henonGeneralized(CBack[counter-1],
					constants);

			min.adjust(true, CBack[counter]);
			max.adjust(false, CBack[counter]);
			
			if ( (CBack[counter].is_greater(space_limit)) ||
					CBack[counter].isNaN()) break;
		}

		// error checking - remove those that are unbounded
		if (CBack[counter].is_greater(space_limit)) continue;

		// now calculate delta_radius ("diff") for the 10 last points
		// if too small, the attractor lacks detail and thus is
		// uninteresting.

		startpoint = CBack[counter-1].getRadius();

		for (sec = counter-2; sec > (counter-10); sec--) 
			diff += (startpoint - CBack[sec].getRadius());

		if ( (diff < lower_difference_limit) || (diff > 
					upper_difference_limit)
				|| isnan(diff)) continue;

		// compute L

		L = X.getLyap(CBack, constants);

		cout << "L: " << L;

		if ( (L < lyap_lower_limit) || isnan(L) ) continue;

		// compute F (rough approximation)

		F = getCorrelationDimension(CBack, ITERATIONS, min, max);
		cout << " F: " << F << " " << lower_F_limit <<"\n";

		// if it falls within limits, calculate more thoroughly
		if ( (F < lower_F_limit) || (F > upper_F_limit) ) continue;
		
		for (counter = ITERATIONS-1; counter < ITERATIONS_THOROUGH;
				counter++)
			CBack[counter] = X.henonGeneralized(CBack[counter-1],
					constants);

		F = getCorrelationDimension(CBack, ITERATIONS_THOROUGH-1, 
				min, max);
	}
}

void convolve(double * array, double & min, double & max, 
		int xsize, int ysize ) {

	// run a convolution kernel to make the output look more 3d-ish

	double process[4];
	double result;

	min = 100;
	max = -100;

	int counter, sec;

	for (sec = 0; sec < ysize; sec++) {
		for (counter = 0; counter < xsize; counter++) {
			process[0] = array[getPos(counter, sec, xsize)];
			process[1] = array[getPos(counter+1, sec, xsize)];
			process[2] = array[getPos(counter, sec+1, xsize)];
			process[3] = array[getPos(counter+1, sec+1, xsize)];

			result = 0;

			// convolution kernel below

			result =  (process[0]*4) + (process[1]*4) +
				  (process[2]) +   (-process[3]*2);

			
			if ((sec < ysize-2) && (counter < xsize-2)) {
				if (result < min) min = result;
				if (result > max) max = result;
			}

			array[getPos(counter, sec, xsize)] = result;
		}
	}
}

// -------------------------- rendering functs ----------------------------- //

void warm_up(int N, Coordinate & renderCoord, Coordinate & min, 
		Coordinate & max, double * constants, Rotation & center) {
	
	// gets viewport size coordinates by running a "warm up" rendering 
	// of N times, recording the minimum and maximum values.

	// rotation & center is used to find the center of mass.

	int counter;
	FractalTypes renderer;

	// before starting the warm up, get initial within bounds by rendering
	// 100 times.

	for (counter = 0; counter < 100; counter++)
		renderCoord = renderer.henonGeneralized(renderCoord, constants);

	// now perform the warmup

	center.setNumberofPoints(N);
	
	for (counter = 0; counter < N; counter++) {
		renderCoord = renderer.henonGeneralized(renderCoord, constants);

		center.addPoint(renderCoord);

		min.adjust(true, renderCoord);
		max.adjust(false, renderCoord);
	}

	// and we're done.
	center.setNumberofPoints(0);


}

/*bool zbuffer_adjust(Coordinate position, Coordinate sourceMin, 
		Coordinate sourceMax,
		Coordinate destTrans, double * array){

	Coordinate translated = real_coord(SourceMin, SourceMax, position, 
			destTrans);

	int linearPos = getPos((int)rint(translated.x), (int)rint(translated.y),
			(int)rint(destTrans.x));

	if ((linearPos < 0) || (linearPos > destTrans.x*destTrans.y) )
		return(false);

	if (translated.z > array[linearPos]) {
		/* more later */


void render(Coordinate renderCoord, Coordinate min, Coordinate max,
		int times, double * constants, double * output_array,
		Coordinate array_size, bool rotate, Rotation rotateBy,
		bool display_progress,
		bool * background_array, bool * shadow_array) {

	// renderCoord is the initial coordinate. min and max sets the 
	// source viewport, and (0,0)-array_size sets the destination
	// viewport. The z values are not changed, but dumped to output_array.

	// if the system is to be rotated before it is dumped (eg for
	// calculating z-buffer shadows), set rotate to true.
	// rotateBy should have its center and rotation coordinates set.

	// if display_progress is on, progress will be shown by a percentage
	// indicator.

	// background_array points to an array, inited to all true, that when
	// exiting will contain true for the pixels that were not visited.

	// ---------------------------------------------------------------- //
	
	// renders a given attractor "times" times, using the constants in
	// *constants, and dumps the result to "output_array"

	// remember to clean up output_array before using this function.

	int counter, location;
	FractalTypes renderer;
	Coordinate srcViewport, destViewport, oldRenderCoord; // normal
	Coordinate rotatedSrcViewport, rotatedDestViewport, rotatedMin,
		rotatedMax;
	int rotatedLocation;
	double fudge_factor = (max.z-min.z)/100;
	int xyzzy = (int)array_size.x*(int)array_size.y;

	rotatedMin = rotateBy.rotate(min);
	rotatedMax = rotateBy.rotate(max);

	double * lightPOV = new double[(int)array_size.x*(int)array_size.y];
	Coordinate * kludge = new Coordinate[xyzzy];
	
	for (counter = 0; counter < array_size.x*array_size.y; counter++) 
		lightPOV[counter] = -100;
	

	for (counter = 0; counter < times; counter++) {
		// display progress if requested

		if (display_progress)
			if (counter % (times/100) == 1) 
				cout << (counter*100)/times << "\t" << flush;
		
		// render
		oldRenderCoord = renderCoord;
		renderCoord = renderer.henonGeneralized(renderCoord, constants);
		
		// set coordinate so we don't mess up renderCoord if rotating
		srcViewport.set(renderCoord.x, renderCoord.y, renderCoord.z);
		
		// get the rotated version for shadowing
		rotatedSrcViewport = rotateBy.rotate(srcViewport);
		
		// translate from source viewport to target viewport, and
		// find where in the array to put it
	
		destViewport = real_coord(min, max, srcViewport, array_size);
		rotatedDestViewport = real_coord(rotatedMin, rotatedMax, 
				rotatedSrcViewport,array_size);

		/* meep meep */
		//destViewport = rotatedDestViewport;
		
		// now this should be cleaned up soon
		location = getPos((int)rint(destViewport.x), 
				(int)rint(destViewport.y), (int)array_size.x);
		rotatedLocation = getPos((int)rint(rotatedDestViewport.x),
					(int)rint(rotatedDestViewport.y),
					(int)array_size.x);
		

		// replace current point iff the point we have is closer to the
		// viewer than the point currently stored is. If so, also set
		// it not to be background.

		// clip off out-of-bounds results
		if (location >= (array_size.x * array_size.y)) continue;
		if (location < 0) continue;
		
		// if clause 1
		if (output_array[location] < destViewport.z) {
			output_array[location] = destViewport.z;
			shadow_array[location] = false; 
				// it's the only way to be sure
			if (background_array != NULL)
				background_array[location] = false;
			kludge[location] = srcViewport;

		}

		// shadow compensation
		if (rotatedLocation >= (array_size.x * array_size.y)) continue;
		if (rotatedLocation < 0) continue;
		
		if (lightPOV[rotatedLocation] < rotatedDestViewport.z) {
			lightPOV[location] = rotatedDestViewport.z;
		}

		if (output_array[location] == destViewport.z) {
			if (lightPOV[rotatedLocation] > 
					(output_array[location]+fudge_factor))
				shadow_array[location] = true;
		}

		// and continue
	}

	for (counter = 0; counter < (array_size.x * array_size.y); counter++){
		if (background_array[counter])
			if (lightPOV[counter] != -100)
				shadow_array[location] = true;
	}

	// QND
	
	for (counter = 0; counter < (array_size.x*array_size.y); counter++) {
		shadow_array[counter] = false;
	}

	for (counter = 0; counter < (array_size.x*array_size.y); counter++) {
		rotatedSrcViewport = rotateBy.rotate(kludge[counter]);

		rotatedDestViewport = real_coord(rotatedMin, rotatedMax,
				 rotatedSrcViewport,array_size);
		 
		rotatedLocation = getPos((int)rint(rotatedDestViewport.x),
				  (int)rint(rotatedDestViewport.y),
				  (int)array_size.x);

		if (output_array[rotatedLocation] < lightPOV[rotatedLocation]) {
				shadow_array[counter] = true;
		}
	}
}

void outcoord(Coordinate a, Coordinate b) {
	cout << rint(a.x) << "(" << rint(b.x) << ")" << "\t";
	cout << rint(a.y) << "(" << rint(b.y) << ")" << "\t";
	cout << a.z << "(" << b.z << ")" << endl;
}

void dispcoord(Coordinate a) {
	cout << "(" << a.x << ", " << a.y << ", " << a.z << ")" << endl;
}

int exposure_function(double value, Coordinate min, Coordinate max,
		int maximum_out){
	double range, output;

	// simple model of a camera: to span intensities

	value = value-min.z;
	//range = log(1/maximum_out)/(min.z-max.z);
	// something is wrong with log
	range = 0.02120459518219654415;

	output = 1-exp(value*range)*(maximum_out);
	
	return ( (int)rint(output) );
}

// ------------------------------ HSV to RGB ---------------------------------

void hsv_to_rgb(double h, double s, double v, unsigned char & r, 
		unsigned char & g, unsigned char & b) {
	// stolen from www.alvyray.com/Papers/hsv2rgb.htm

	Coordinate real(0, 0, 0); // H, S, V

	real.set(h*6/255, s/255, v/255);

	/*float realH = h*6/255;
	double realS = s/255;
	double realV = v/255;
	double blah = 1.3734;*/

	double m, n, f;
	int i;

	Coordinate out (0,0,0);

	i = (int)floor(real.x);
	f = real.x - i;

	if ( !(i&1) ) f = 1-f;

	m = real.z * (1-real.y);
	n = real.z * (1 - real.y * f);

	switch(i) {
		case 6:
		case 0: out.set(real.z, n, m); break;
		case 1: out.set(n, real.z, m); break;
		case 2: out.set(m, real.z, n); break;
		case 3: out.set(m, n, real.z); break;
		case 4: out.set(n, m, real.z); break;
		case 5: out.set(real.z, m, n); break;
	}

	r = (int)rint(out.x*255);
	g = (int)rint(out.y*255);
	b = (int)rint(out.z*255);
}
	
int main(int argc, char * * argv) {

	Coordinate C, Q, real, realrot, rotmin, rotmax, min, max, center;
	Coordinate output(XZ, YZ-1, ZZ);
	Coordinate identity(0,0,0);
	Coordinate w_up[WARMUP];

	// rotation coordinates/functions
	//Rotation rotFunct((double)r_phi, (double)r_theta, (double)r_eta);
	Rotation rotFunct(degToRadian(r_phi), degToRadian(r_theta), 
			degToRadian(r_eta));
	
	unsigned int arraysize = XZ*YZ+1;

	if (argc < 2) {
		cout << "usage: "<<argv[0] << 
			" [number of viable attractors to skip]\n";
		return(-1);
	}

	srandom(atoi(argv[2]));

	double * temparray = new double [arraysize];
	bool * isbackground = new bool [arraysize];
	bool * isshadow = new bool [arraysize];
	
	FractalTypes X;
	unsigned int counter, p, q, oldq;
	double diff, F, L;
	double cMi, cMa, cMed;

	double constants[57];

	int sec = 0;
	int numerum = atoi(argv[1]);

	FILE * output_file;

	// prepare data

	for (counter = 0; counter < arraysize; counter++) {
		temparray[counter] = -100;
		isbackground[counter] = true;
		isshadow[counter] = false;
	}

	C.set(0.1, 0.1, 0.1);
	diff = 0;

	if (numerum < 1) numerum = 1;

	// find attractors

	for (sec = 0; sec < numerum; sec++)
		findAttractors(10, 0.005, 10, LOW_FDIM_LIMIT, HIGH_FDIM_LIMIT, 
				LOW_LYAP_LIMIT, constants, min, max, F, L);

	// print data about attractor

	cout << "Correlation dimension: " << F << " ( "
	   << (double)((F-LOW_FDIM_LIMIT)*100/(HIGH_FDIM_LIMIT-LOW_FDIM_LIMIT)) 
			<< "% )\n";

	cout << "Lyapunov Exponent: " << L << "\n";
	
	X.revPop(constants, 36);

	// warmup to get accurate min and max values
	// some attractors make this strategy fail because some very few points
	// (1% or less) are set to way off the main attractor... to fix later?
	
	C = X.henonGeneralized(C, constants);
	w_up[0] = C;

	// set values that we are sure is going to be changed
	rotmin.set(200, 200, 200);
	min.set(200,200,200);
	max.set(-200,-200,-200);
	rotmax.set(-200, -200, -200);
	
	warm_up(WARMUP, w_up[0], min, max, constants, rotFunct);

	/*for (counter = 1; counter < WARMUP; counter++) {
		w_up[counter] = X.henonGeneralized(w_up[counter-1], constants);
		
		if (w_up[counter].isNaN()) {
			cout << "Unbounded at iteration" << counter << ".\n";
			return(-1);
		}
		
		min.adjust(true, w_up[counter]);
		max.adjust(false, w_up[counter]);
	}*/

	// multiply min and max to prevent attractor from getting 
	// squashed to the edges

	/*center.x = min.x+((max.x - min.x)/2);
	center.y = min.y+((max.y - min.y)/2);
	center.z = min.z+((max.z - min.z)/2);*/

	min *= 1.05;
	max *= 1.05;

	
	/* get new max and min coords
	#ifdef SHADOWS
	for (counter = 1; counter < WARMUP; counter++) {
		C = rotateCoordinate(costheta, sintheta, cosphi, sinphi,
				coseta, sineta, center, w_up[counter]);
		
		rotmin.adjust(true, C);
		rotmax.adjust(false, C);
	}
	#endif*/
	
	rotmin = min;
	rotmax = max;

	render(Q, min, max, LIMIT, constants, temparray, output, false,
			rotFunct, true, isbackground, isshadow);

	for (counter = 0; counter <= arraysize; counter++) {
		if (temparray[counter] == -100) temparray[counter] 
			= (min.z*1.05);	// background is behind the entire 
					// object
	}

	//convolve(temparray, cMi, cMa, XZ, YZ); 

	cMed = (cMa-cMi)/(ZZ-ZZ_L);

	output_file = fopen("out.raw.col", "wb");
	
	unsigned char r[3];

	for (counter = 0; counter < arraysize; counter++) {

		if (isbackground[counter]) {
			r[0] = 64;			// R
			r[1] = 64;			// G	 blue background
			r[2] = 128;			// B
		} else {
			if (temparray[counter] > 255) cout << "." << flush;
			if (!isshadow[counter]) {
			hsv_to_rgb(HUE, SAT_BRIGHT, temparray[counter],
					r[0], r[1], r[2]);
			} else {
				hsv_to_rgb(HUE, SAT_SHADOW, temparray[counter],
						r[0], r[1], r[2]);
			}
			//}
			// variant of real_coordinate, but only z axis
			// so we render a height map
			//r[0] = ZZ_L + (int)rint((temparray[counter]-cMi)/cMed);

			/*r[0] = ZZ_L + exposure_function(temparray[counter], 
					identity, output, 255-ZZ_L);*/
				
			
			/*r[1] = r[0];
			r[2] = r[0];*/
		}
		

		/*if (isshadow[counter]) {
			r[0] = (int)rint((double)(r[0] * 0.30));
			r[1] = (int)rint((double)(r[1] * 0.30));
			r[2] = (int)rint((double)(r[2] * 0.30));
		}*/

		
		fwrite(r, 3, 1, output_file);
	}

	fclose(output_file);
	cout << "\n";
}