/****************************
 * MANDELBROT SET
 *
 * Command line tool to produce images of the Mandelbrot set.
 *
 * Copyright (c) 2008, Per Liedman (per at liedman.net)
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *     * Redistributions of source code must retain the above copyright
 *       notice, this list of conditions and the following disclaimer.
 *     * Redistributions in binary form must reproduce the above copyright
 *       notice, this list of conditions and the following disclaimer in the
 *       documentation and/or other materials provided with the distribution.
 *     * Neither the name of the <organization> nor the
 *       names of its contributors may be used to endorse or promote products
 *       derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY PER LIEDMAN ``AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL PER LIEDMAN BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <stdio.h>
#include <stdlib.h>

float *mandelbrot(float x1, float xi1, float x2, float xi2, int w, int h, int max_iter, float *buffer) {
	float step_x = (x2 - x1) / w;
	float step_xi = (xi2 -xi1) / h;	
	float x = x1;
	float xi = xi1;
	int i;
	int j;
	int iter;
	float x_c;
	float xi_c;
	float x_temp;
	float xi_temp;
	
	for (j = 0; j < h; j++) {
		x = x1;

		for (i = 0; i < w; i++) {
			x_c = x;
			xi_c = xi;
		
			for (iter = 0; x_c * x_c + xi_c * xi_c < 4 && iter < max_iter; iter++) {
				x_temp = x_c * x_c - xi_c * xi_c + x;
				xi_temp = 2 * x_c * xi_c + xi;
				
				x_c = x_temp;
				xi_c = xi_temp;
			}
			
			*(buffer++) = iter;
			
			x = x + step_x;		
		}
		
		xi = xi+ step_xi;
	}
	
	return buffer;
}

void pgm(int *buffer, int w, int h, int maxval, FILE *f) {
	int x, y;
	int c_count = 0;

	fprintf(f, "P2\n");
	fprintf(f, "%d %d %d\n", w, h, maxval);
	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			c_count += fprintf(f, "%d ", *(buffer++));
			
			if (c_count > 64) {
				fprintf(f, "\n");
				c_count = 0;
			}
		}
	}
	
	fprintf(f, "\n");
}

void ppm(float *buffer, int w, int h, void (color_function(float, unsigned char*)), FILE *f) {
	int x, y;
	unsigned char o_buffer[16384 * 3];
	int o_buf_c = 0;
	unsigned char *o_buf_p = o_buffer;

	fprintf(f, "P6\n");
	fprintf(f, "%d %d 255\n", w, h);
	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			color_function(*(buffer++), o_buf_p);
			o_buf_p += 3;
			o_buf_c++;
			if (o_buf_c >= 16384) {
				fwrite(o_buffer, sizeof(unsigned char) * 3, o_buf_c, f);
				o_buf_c = 0;
				o_buf_p = o_buffer;
			}
		}
	}
	
	fwrite(o_buffer, sizeof(int) * 3, o_buf_c, f);
}

const float R_CENTER = 0.5f;
const float G_CENTER = 0.6f;
const float B_CENTER = 0.05f;
const float R_COFF = 5.0f;
const float G_COFF = 8.0f;
const float B_COFF = 4.0f;

void color_function(float v, unsigned char *result) {
	if (v > 255.0) {
		*result++ = 0;
		*result++ = 0;
		*result++ = 0;
	} else {
		float f = v / 255.0;

		float r_val = f - R_CENTER;
		float g_val = f - G_CENTER;
		float b_val = f - B_CENTER;
	
		float r = ((1 - R_COFF * r_val * r_val) * 255.0f);
		float g = ((1 - G_COFF * g_val * g_val) * 255.0f);
		float b = ((1 - B_COFF * b_val * b_val) * 255.0f);
	
		r = r < 0 ? 0 : r;
		g = g < 0 ? 0 : g;
		b = b < 0 ? 0 : b;
		r = r > 255 ? 255 : r;
		g = g > 255 ? 255 : g;
		b = b > 255 ? 255 : b;
	
		*result++ = (unsigned char)r;
		*result++ = (unsigned char)g;
		*result++ = (unsigned char)b;
	}
}

int main(int argc, char **argv) {
	if (argc < 3) {
		fprintf(stderr, "Missing arguments.\n");
		fprintf(stderr, "SYNOPSIS:\n\n");
		fprintf(stderr, "%s [X] [Xi] [DIAMETERX] [MAXITER] [WIDTH] [HEIGHT].\n", argv[0]);
		return 1;
	}

	float x_c, xi_c, diameter_x;
	float x1, xi1, x2, xi2;
	float aspect_ratio;
	int max_iter, w, h;
	FILE *output = stdout;

	x_c = atof(argv[1]);
	xi_c = atof(argv[2]);
	diameter_x = atof(argv[3]);
	max_iter = atoi(argv[4]);
	w = atoi(argv[5]);
	h = atoi(argv[6]);
	aspect_ratio = (float)w / (float)h;
	
	if (argc > 7) {
		output = fopen(argv[7], "wb");
		
		if (output == NULL) {
			fprintf(stderr, "Couldn't open file \"%s\".\n", argv[7]);
			return 1;
		}
	}
	
	x1 = x_c - diameter_x / 2.0f;
	xi1 = xi_c - diameter_x / 2.0f / aspect_ratio;
	x2 = x_c + diameter_x / 2.0f;
	xi2 = xi_c + diameter_x / 2.0f / aspect_ratio;

	float *buffer = malloc(w * h * sizeof(float));
	
	if (buffer != NULL) {
		mandelbrot(x1, xi1, x2, xi2, w, h, max_iter, buffer);
		ppm(buffer, w, h, color_function, output);
		
		return 0;
	} else {
		fprintf(stderr, "Couldn't allocate image buffer of size %d bytes\n", w * h * sizeof(int));
		
		return 1;
	}
}