/*
XANDER'S MANDELBROT PROGRAM OF DOOM
REVISION DATE:  12 August 2011

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:		name of the output file (include file extension)
	line 9:		name of the color file (include file extension)
	line 10:	fractal to generate
					m - Mandelbrot / multibrot set (z^d+c)
					M - Mandelbrot / multibrot set with NIC colors
					j - Julia set (z^d + K)
					J - Julia set with NIC colors
					s - burning ship fractal
					u - Mandelbrot / multibrot inversion ((1/z)^d+1/c)
					U - Mandelbrot / multibrot inversion with NIC colors
	line 12+:	additional parameters
					Re(K) and Im(K) for j or J
					number of colors for M and J
*/

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

using namespace std;

const int MAXCOLORS = 12800;

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);
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 '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 RETURNS A PIXEL COLOR USING THE ESCAPE TIME ALGORITHM.
	THE ALGORITHM GIVES UP AFTER 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 A PIXEL COLOR USING THE ESCAPE TIME ALGORITHM.  THE
	ALGORITHM GIVES UP AFTER 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;
	} else

	return iteration;
}

/*	DETERMINES WHETHER A GIVEN POINT IN THE COMPLEX PLANE IS IN THE dTH DEGREE
	BURNING SHIP FRACTAL, AND RETURNS A PIXEL COLOR USING THE ESCAPE TIME
	ALGORITHM.  THE ALGORITHM GIVES UP AFTER 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;
}

/*	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;
}