/*
XANDER'S MANDELBROT PROGRAM OF DOOM
REVISION DATE:  10 January 2012

Expected input file format:
	line 1:		exponent (degree of the Julia or Mandelbrot set to generate)
	line 2:		x center
	line 3:		y center
	line 4:		zoom level
	line 5:		aspect ratio (width:height)
	line 6:		horizontal size of the output file
	line 7:		maximum number of iterations
	line 8:		number of colors
	line 9:		name of the output file (include file extension)
	line 10:	name of the color file (include file extension)
	line 11:	fractal to generate
					m - Mandelbrot / multibrot set (z^d+c)
					j - Julia set (z^d + K)
					s - burning ship fractal
					u - Mandelbrot / multibrot inversion ((1/z)^d+1/c)
					c - Collatz fractal
				use capital letters for NIC coloring
	line 12+:	additional parameters
					Re(K) and Im(K) for j or J
*/

#include<cstdlib>
#include<iostream>
#include<fstream>
#include<string>
#include<cmath>

using namespace std;

const int MAXCOLORS = 12800;
const double PI = 3.14159265;

int fact(int);
int choose(int, int);
void z_pow(double [], int);
int mandelbrot(double [], int, int);
int julia(double [], double [], int, int);
int ship(double [], int, int);
int collatz(double [], int);
void initColors(char [][3], int, ifstream &);
void bmpHeaders(int, int, int, ofstream &);
void writeValue(int, int, ofstream &);

int main(int argc, char *argv[]) {
	srand(time(NULL));

	ofstream fout;
	ifstream fin;
	int padding, color, d, xsize, ysize, numIterations, numColors;
	string input, parFile, imgName, colorFile;
	double xcen, ycen, aspect, zoom, xmin, xmax, ymin, ymax, step, z[2], K[2], ff=0, tmp, D;
	char colors[MAXCOLORS+1][3], fractal;
	bool invalid=true;
	time_t start, end;

	if (argc > 1) {
		parFile = argv[1];

		//	READ IN THE PARAMETER FILE
		cout << "Parameter file is:  " << parFile  << endl;

		fin.open(("par_files/"+parFile+".txt").c_str());
		if(fin.fail()) {
			cerr << "ERROR: The parameter file  'par_files/" << parFile << ".txt' could not be opened.\nAttempting to use 'parameters-default.txt'.\n";
			
			fin.open("par_files/parameters-default.txt");
			if(fin.fail()) {
				cerr << "ERROR: 'parameters-default.txt' could not be opened.\nAborting.\n";
				exit(1);
			}
		}
		
		fin >> d >> xcen >> ycen >> zoom >> aspect >> xsize >> numIterations >> numColors >> imgName >> colorFile >> fractal;

		switch(fractal) {
			case 'j':
				fin >> K[0] >> K[1];
				break;
			case 'J':
				fin >> K[0] >> K[1];
				break;
		}

		if(numColors > MAXCOLORS) {
			numColors = MAXCOLORS;
		}
		
		ycen = -ycen;
		xmin = xcen - zoom;
		xmax = xcen + zoom;
		ymin = ycen - (zoom * 1.0/aspect);
		ymax = ycen + (zoom * 1.0/aspect);
		
		fin.close();
		
		step = (xmax-xmin) / xsize;
		ysize = (ymax - ymin) / step;
	}
	else {
		cerr << "ERROR: No parameter file was specified.\nAborting.\n";
		exit(0);
	}
	
	fout.open(("output/"+imgName).c_str());
	if(fout.fail()) {
		cerr << "ERROR: The output file 'output/" << imgName << "'  could not be opened.\nAttempting to use 'temp.bmp'.\n";
		
		fout.open("output/temp.bmp");
		if(fout.fail()) {
			cerr << "ERROR: 'temp.bmp' could not be opened.\nAborting.\n";
			exit(1);
		}
	}

	fin.open(("color_files/"+colorFile).c_str());
	if(fin.fail()) {
		cerr << "ERROR: The color file 'color_files/" << colorFile << "' could not be opened.\nAttempting to use 'colors-default.bmp'.\n";

		fin.open("color_files/colors-default.bmp");
		if(fin.fail()) {
			cerr << "ERROR: 'colors-default.bmp' could not be opened.\nAborting.\n";
			exit(1);
		}
	}

	start = time(NULL);

	initColors(colors, numColors, fin);
	fin.close();

	writeValue(54,0,fout);

	padding = (4-((xsize*3)%4))%4;	// end of line padding for the bitmap

	xmin += (step/2);
	ymin += (step/2);

	//	THE HEART OF THE PROGRAM
	//		for each pixel, we determine whether it is in the set or not
	//		then color it appropriately
	for(int y = ysize; y >0; --y) {
		for(int x = 0; x < xsize; ++x) {

			z[0] = xmin + (x * step);
			z[1] = ymin + (y * step);

			switch(fractal) {
				case 'j':
				case 'J':
					color = julia(z, K, d, numIterations);
					break;
				case 's':
					color = ship(z, d, numIterations);
					break;
				case 'u':
				case 'U':
					D=(z[0]*z[0])+(z[1]*z[1]);
					z[0] = z[0]/D;
					z[1] = -z[1]/D;
					color = mandelbrot(z, d, numIterations);
					break;
				case 'c':
				case 'C':
					color = collatz(z, numIterations);
					break;
				case 'm':
				case 'M':
				default:
					color = mandelbrot(z, d, numIterations);
					break;
			}

			// NIC coloring, if specified by a capital letter
			if(fractal >= 'A' && fractal <= 'Z' && color != 0) {
				tmp = color + (log(log(81)) - (log(log(sqrt(z[0]*z[0]+z[1]*z[1]))))/log(d));
				tmp = tmp * ((double) numColors / (double) numIterations);
				color = (int) tmp;
			}
			
			if(color <= 0) {
				color = 0;
			} else {
				color %= (numColors-1);
				++color;
			}

			fout << colors[color][0] << colors[color][1] << colors[color][2];
		}

		writeValue(padding,0,fout);
	}
	
	fout.seekp(ios::beg);
	bmpHeaders(xsize, ysize, padding, fout);

	fout.close();

	end = time(NULL);
	
	printf("Render Time: %5d seconds\n", (int) (end-start));

	return 0;
}

/*	COMPUTES n! */
int fact(int n) {
	int res=1;
	
	for(int i=2; i<=n; ++i) {
		res *= i;
	}
	
	return res;
}

/*	COMPUTES C(n,k) */
int choose(int n, int k) {
	int res = 1;
	
	if(k>n/2) {		//  C(n,k) = C(n,n-k), so work with the smaller of the two
		k = n-k;
	}
	
	for(int i=n; i>=n-k+1; --i) {
			res *= i;
	}
	
	return res/fact(k);
}

/*	COMPUTES z^d WHERE z IS COMPLEX AND d IS AN INTEGER */
void z_pow(double z[2], int d) {
	double res[2], delta;
	
	res[0]=0; res[1]=0;		// the real and imaginary parts of the result
	
	for(int i=0; i<=d; ++i) {
		delta = choose(d,i)*pow(z[0],d-i)*pow(z[1],i);
		switch(i%4) {
			case 0:
				res[0] += delta;
				break;
			case 1:
				res[1] += delta;
				break;
			case 2:
				res[0] -= delta;
				break;
			case 3:
				res[1] -= delta;
				break;
		}
	}
	
	z[0] = res[0];
	z[1] = res[1];

	return;
}

/*	DETERMINES WHETHER A GIVEN POINT IN THE COMPLEX PLANE IS IN THE dTH DEGREE
	MULTIBROT SET, AND EITHER RETURNS THE NUMBER OF ITERATIONS UNTIL THE ESCAPE
	CRITERION IS MET OR numIterations IF THE ESCAPE CRITERION IS NOT MET IN
	FEWER THAN numIterations ITERATIONS. */
int mandelbrot(double z[2], int d, int numIterations) {
	int iteration = 1;
	bool inSet = true;
	double c[2];
	
	c[0]=z[0];
	c[1]=z[1];
	
	if((z[0]*z[0])+(z[1]*z[1]) > 4) {
		inSet = false;
	}

	while(inSet && iteration < numIterations) {
		z_pow(z,d);
		z[0] += c[0];
		z[1] += c[1];
		++iteration;
	
		if((z[0]*z[0])+(z[1]*z[1]) > 4) {
			inSet = false;
		}
	}
	
	if(inSet) {
		iteration = 0;
	}

	return iteration;
}

/*	DETERMINES WHETHER A GIVEN POINT IN THE COMPLEX PLANE IS IN THE dTH DEGREE
	JULIA SET, AND RETURNS THE NUMBER OF ITERATIONS UNTIL THE ESCAPE
	CRITERION IS MET OR numIterations IF THE ESCAPE CRITERION IS NOT MET IN
	FEWER THAN numIterations ITERATIONS. */
int julia(double z[2], double K[2], int d, int numIterations) {
	int iteration = 1;
	bool inSet = true;
	double c[2];
	
	c[0]=z[0];
	c[1]=z[1];
	
	if((z[0]*z[0])+(z[1]*z[1]) > 4) {
		inSet = false;
	}

	while(inSet && iteration < numIterations) {
		z_pow(z,d);
		z[0] += K[0];
		z[1] += K[1];
		++iteration;
	
		if((z[0]*z[0])+(z[1]*z[1]) > 4) {
			inSet = false;
		}
	}
	
	if(inSet) {
		iteration = 0;
	}

	return iteration;
}

/*	DETERMINES WHETHER A GIVEN POINT IN THE COMPLEX PLANE IS IN THE dTH DEGREE
	BURNING SHIP FRACTAL, AND RETURNS THE NUMBER OF ITERATIONS UNTIL THE ESCAPE
	CRITERION IS MET OR numIterations IF THE ESCAPE CRITERION IS NOT MET IN
	FEWER THAN numIterations ITERATIONS. */
int ship(double z[2], int d, int numIterations) {
	int iteration = 1;
	bool inSet = true;
	double c[2];
	
	c[0]=z[0];
	c[1]=z[1];
	
	if((z[0]*z[0])+(z[1]*z[1]) > 4) {
		inSet = false;
	}

	while(inSet && iteration < numIterations) {
		z[0] = abs(z[0]);
		z[1] = abs(z[1]);
		z_pow(z,d);
		z[0] += c[0];
		z[1] += c[1];
		++iteration;
	
		if((z[0]*z[0])+(z[1]*z[1]) > 4) {
			inSet = false;
		}
	}
	
	if(inSet) {
		iteration = 0;
	}

	return iteration;
}

/*	DETERMINES WHETHER A GIVEN POINT IN THE COMPLEX PLANE IS IN THE COLLATZ
	FRACTAL SET, AND RETURNS THE NUMBER OF ITERATIONS UNTIL THE ESCAPE
	CRITERION IS MET OR numIterations IF THE ESCAPE CRITERION IS NOT MET IN
	FEWER THAN numIterations ITERATIONS. */
int collatz(double z[2], int numIterations) {
	int iteration = 1;
	bool inSet = true;
	double cosz[2],tmp[2];

	if((z[0]*z[0])+(z[1]*z[1]) > 16384) {
		inSet = false;
	}

	while(inSet && iteration < numIterations) {
		cosz[0] = cos(PI*z[0])*cosh(PI*z[1]);
		cosz[1] = -sin(PI*z[0])*sinh(PI*z[1]);

		tmp[0] = (z[0]*cosz[0])-(z[1]*cosz[1]);
		tmp[1] = (z[0]*cosz[1])+(z[1]*cosz[0]);
		
		z[0] = (2+7*z[0]-2*cosz[0]-5*tmp[0])/4;
		z[1] = (7*z[1]-2*cosz[1]-5*tmp[1])/4;

		++iteration;
	
		if((z[0]*z[0])+(z[1]*z[1]) > 16384) {
			inSet = false;
		}
	}
	
	if(inSet) {
		iteration = 0;
	}

	return iteration;
}

/*	READS A BITMAP FILE INTO AN ARRAY THAT WILL BE USED TO COLOR THE FINAL IMAGE
	*KNOWN BUG*
		IF THE INPUT BITMAP'S WIDTH IS NOT A MULTIPLE OF 4, THE RESULTS WILL BE
		UNPREDICTABLE */
void initColors(char colors[][3], int numIterations, ifstream &fin) {
	fin.seekg(54);
	
	for(int i=0; i<=numIterations; ++i) {
		for(int j=0; j<=2; ++j) {
			fin.get(colors[i][j]);
		}
	}

	return;
}

/*	WRITES THE BITMAP HEADERS TO THE OUTPUT FILE */
void bmpHeaders(int xsize, int ysize, int padding, ofstream &fout) {
	int bmpSize, numPixels;
	numPixels = (xsize+padding)*ysize;
	bmpSize = numPixels*3+54;
	
	// Magic Number
	fout << "BM";
	
	// Size of BMP File
	writeValue(4,bmpSize,fout);
	
	// Unused
	writeValue(4,0,fout);
	
	// Offset
	writeValue(4,54,fout);
	
	// Header Size
	writeValue(4,40,fout);
	
	// Height and Width
	writeValue(4,xsize,fout);
	writeValue(4,ysize,fout);
	
	// Color Planes
	writeValue(2,1,fout);
	
	// Bits per Pixel
	writeValue(2,24,fout);
	
	// Compression
	writeValue(4,0,fout);
	
	// Raw BMP Size
	writeValue(4,bmpSize-54,fout);
	
	// Horizontal and Vertical Resolution
	writeValue(4,2835,fout);
	writeValue(4,2835,fout);

	// Colors in Palette
	writeValue(4,0,fout);
	
	// Important Colors
	writeValue(4,0,fout);
	
	return;
}

/*	WRITES A VALUE TO A FILE 
	(I honestly don't remember exactly what this does anymore. :\ ) */
void writeValue(int num, int value, ofstream &fout) {
	int tmp;
	
	for(int i=0; i<num; i++) {
		if(value > 255) {
			tmp = value / 256;
			tmp = value * 256;
			fout << (char) (value-tmp);
			value /= 256;
		}
		else if(value >= 0 && value <=255) {
			fout << (char) value;
			value = -1;
		}
		else {
			fout << (char) 0;
		}
	}
	return;
}