// renders the well known Mu Molecule (aka Mandelbrot set)

// uses Plotter.

// TODO: replace with
//	Plotter			[ DONE ]
//	Coordinate / BasicCoordinate

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

#define PI	3.1415926535

#define NUMBER_OF_ITERATIONS 1000	// if it hasn't escaped by this number
					// of iterations, get out
#define MAGNITUDE_LIMIT      2		// mandelbrot limit

// output file data
#define XSIZE 400
#define YSIZE 400

// input data
#define X_START -2
#define X_END 0.9
#define Y_START -1.4
#define Y_END 1.4

#define SATURATION 1
#define VALUE 1

// program specific optimization - used to cut away a square root and thus
// be faster.
// (if sqrt(x^2+y^2) > limit) becomes (if x^2+y^2 > limit*limit)
#define realMagLimit	     MAGNITUDE_LIMIT*MAGNITUDE_LIMIT


// color space transformations

void HSV_to_RGB_internal(double h, double s, double v, double & R, double & G,
		double & B) {

	// H is given on [0, 6] or. S and V are given on [0, 1]. 
	// RGB are each returned on [0, 1]. 
	double m, n, f; 
	int i; 

	if (s == 0) { R = v; G = v; B = v; return; } 

	i = (int)floor(h); 
	f = h - i; 
	if ( !(i&1) ) f = 1 - f; // if i is even 
	m = v * (1 - s); 
	n = v * (1 - s * f); 

	switch (i) { 
		case 6: 
		case 0: R = v; G = n; B = m; break; 
		case 1: R = n; G = v; B = m; break;
		case 2: R = m; G = v; B = n; break;
		case 3: R = m; G = n; B = v; break;
		case 4: R = n; G = m; B = v; break;
		case 5: R = v; G = m; B = n; break;
	} 

}

void HSV_to_RGB(double * hsv, int * rgb, double scaling_factor) {
	// NO ERROR CHECKING.
	
	//hsv[0] = hue, hsv[1] = saturation, hsv[2] = value
	//rgb[0] = red, rgb[1] = green, rgb[3] = blue
	//scaling factor is max value for each of the rgb values

	// hsv[0] = [0 ... 6], rest are [0 ... 1]
	
	double normalized_R, normalized_G, normalized_B;

	HSV_to_RGB_internal(hsv[0], hsv[1], hsv[2], normalized_R, normalized_G,
			normalized_B);
	
	/*cout << hsv[0] << "\t" << hsv[1] << "\t" << hsv[2] << "\n";
	cout << normalized_R << "\t" << normalized_G << "\t" << normalized_B;
	cout << "\n";*/
	
	rgb[0] = (int)(rint(normalized_R * scaling_factor));
	rgb[1] = (int)(rint(normalized_G * scaling_factor));
	rgb[2] = (int)(rint(normalized_B * scaling_factor));
}

//----------------------------------------------------------------------------//

void calculate(double & x, double & y, double C_real, double C_imaginary) {

	double old_x = x;

	// x is real part, y is imaginary
	// (denote the complex number x + yi = z
	//  then we are iterating z^2 + C

	x = (old_x*old_x) - (y*y) + C_real;
	y = 2*old_x*y + C_imaginary;
}

double get_magnitude(double x, double y) {
	        // radius
	        return(sqrt((x*x)+(y*y)));
}

bool magnitude_greater(double x, double y, double limit_squared) {
	        // optimized trick, but makes code unreadable
	        return ( ((x*x)+(y*y)) > limit_squared);
}

double getCount(double x, double y, double max) {
	// this does the hard part
	// set C = (x,y), then start iterating with x0,y0 = (0,0)
	// if we escape, return the number of iterations if less than max
	// (otherwise return max). if we don't, return 0.

	// returns normalized doubles [0..1]

	// See descriptions of escape-time on the web for information on the
	// algorithm.

	double xn = x, yn = y, oldxn, oldyn, logabs;
	double angle;

	for (int counter = 0; counter < NUMBER_OF_ITERATIONS; counter++) {

		if (magnitude_greater(xn, yn, realMagLimit)) {
			if (counter < max) {
				return (counter/max);
				// uncomment the two lines below to add angles
				// to the information conveyed by color
				// this would result in a sort of "true-color
				// mandelbrot"
				/*angle = (atan(yn/xn)+PI/2)/PI;
				return((counter+angle)/max);*/
			}
			else			return(1);
		}

		oldxn = xn; oldyn = yn;
		calculate(xn, yn, x, y);
	}

	// magnitude was not greater. return zero.
	return(0);
}

double processHue(double input) {

	// Originally, mandelbrot renders were mapped to the standard
	// VGA set which was basically logarithmic (each color was high-
	// contrast wrt the previous) for the first 15 colors, then linear,
	// but since we're now mapping to true-color, the high-contrast
	// effect must be approximated somehow - here by sines.

	// the advantage of doing this is that the function will be continuous
	// if we add "distance estimation" algorithms later on.

	return((1+sin(input*2*PI))/2);
		
}


main() {

	double incrementX; // X amount in source space that corresponds to one
			   // pixel to the right in output space
	double incrementY; // ditto, but with Y and down, respectively

	double xpos, ypos; // counters

	// clean up array
	double * output = new double [XSIZE*YSIZE];
	for (int ix = 0; ix < (XSIZE*YSIZE); ix++) output[ix] = 0;

	Plotter tool;

	// set required plotter parameters
	tool.setSourceBoundaries(X_START, Y_START, X_END, Y_END);
	tool.setTargetBoundaries(0, 0, XSIZE, YSIZE);
	tool.setTargetArray(output);
	
	incrementX = tool.getStepsX(XSIZE);	// get X increment rate
	incrementY = tool.getStepsY(YSIZE);	// get Y increment rate

	// iterate for all (x,y)

	for (ypos = Y_START; ypos < Y_END; ypos += incrementY) {
		
		cout << ypos << "\n";
		
		for (xpos = X_START; xpos < X_END; xpos += incrementX) 
			tool.plot(xpos, ypos, getCount(xpos, ypos, 
						256-1));
	}

	// dump plot
	
	ofstream outputFile("qnd.raw");
	double hsv[3];
	int rgb[3];

	// set HSV values
	
	hsv[1] = SATURATION;
	hsv[2] = VALUE;

	for (int counter = 0; counter < XSIZE*YSIZE; counter++) {
		// set hue
		hsv[0] = processHue(output[counter])*6;

		// convert to rgb and dump
		HSV_to_RGB(hsv, rgb, 255);
		outputFile << (char)rgb[0] << (char)rgb[1] << (char)rgb[2];
	}

	outputFile.close();
			
}