/* mandelbrot
* Copyright (C) 2008 Wikimedia Foundation
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

#include <stdio.h>
#include <iostream>
#include <cmath>
using namespace std;


/********************************** settings **********************************/
const int w = 4000; // image width
const int h = 1500; // image height
const double mx = 0.7, my = 0.0; // relative position of the origin 
const int u = 0.3 * w; // unit
const int iters = 2000; // maximum number of iterations
const bool smooth = true; // flattens the color steps
const int subpix = 16; // must be >= 1; decrease for much faster evaluation!
/******************************************************************************/

double fc(double x)
{	// makes the color slightly smoother
	if (x <= 0.0) return(0.0);
	if (x >= 1.0) return(1.0);
	return(2.0 * (1.0  + x * sqrt(x) - (1 - x) * sqrt(1 - x)) - 3.0 * x);
}

void color(double x, double rgb[])
{
	x += 0.52;	// arbitrary phase shift
	x -= (int)x;
	if (x < 0)
		x += 1.0;
	x *= 6;
	int xx = x;
	x -= xx;
	switch (xx)
	{
		case 0:	rgb[2] = 0; rgb[0] = 1.0; rgb[1] = fc(1.0 - x);
			break;
		case 1:	rgb[1] = 0; rgb[0] = 1.0; rgb[2] = fc(x);
			break;
		case 2:	rgb[1] = 0; rgb[2] = 1.0; rgb[0] = fc(1.0 - x);
			break;
		case 3:	rgb[0] = 0; rgb[2] = 1.0; rgb[1] = fc(x);
			break;
		case 4:	rgb[0] = 0; rgb[1] = 1.0; rgb[2] = fc(1.0 - x);
			break;
		case 5:	rgb[2] = 0; rgb[1] = 1.0; rgb[0] = fc(x);
			break;
	}
}

double flatten(double x0, double y0, double x, double y)
{	// creates a continuous gradient instead of color steps
	int steps = 0;
	double sqrxy = x * x + y * y;
	if (sqrxy >= 4.0)
	{
		while (sqrxy < 65536)	// high value >> 2
		{
			double x_old = x;
			x = x * x - y * y;
			y = 2 * x_old * y;
			x += x0;
			y += y0;
			steps++;
			sqrxy = x * x + y * y;
		}
		double m = steps + 1.5 - log2(log2(sqrxy));
		return(m);
	}
	return(0.0);
}

bool outcircle(double cx, double cy, double r, double x, double y)
{	// checks if (x,y) is outside the circle around (x0,y0) with radius r
	x -= cx;
	y -= cy;
	if (x * x + y * y > r * r)
		return(true);
	return(false);
}

int main()
{
	// open image file
	FILE *image_file = fopen("mandelbrot.ppm", "wb");
	fprintf(image_file, "P6\n%d %d\n255\n", w, h);
	
	// iterate pixels
	for(int row = 0; row < h; row++) {
	for(int column = 0; column < w; column++)
	{		
		double sum[3] = { 0.0, 0.0, 0.0 };
		// iterate subpixels
		for (int i = 0; i < subpix; i++) {
		for (int j = 0; j < subpix; j++)
		{
			// pixel -> coordinates
			double x0 = (-w * mx + column
					+ (j + 0.5) / subpix) / u;
			double y0 = (h * (1.0 - my)
					- row - (i + 0.5) / subpix) / u;
			double x = x0, y = y0;
			int n = -1;
			
			// skip values we know they are inside
			if ((outcircle(-0.11, 0.0, 0.63, x0, y0) || x0 > 0.1)
				&& outcircle(-1.0, 0.0, 0.25, x0, y0)
				&& outcircle(-0.125, 0.744, 0.092, x0, y0)
				&& outcircle(-1.308, 0.0, 0.058, x0, y0)
				&& outcircle(0.0, 0.25, 0.35, x0, y0))
			{
				for (int it = 0; it <= iters; it++)
				{	// iteration steps
					/********* fractal generation ********/
					double x_old = x;
					x = x * x - y * y;
					y = 2 * x_old * y;
					x += x0;
					y += y0;
					/*************************************/
					if (x * x + y * y > 4)
					{	// value is outside
						n = it;
						break;
					}
				}
			}
			
			double subpixel[3];
			if (n < 0)
			{	// black color inside
				for (int c = 0; c < 3; c++)
					subpixel[c] = 0;
			}
			else
			{
				double m = n;
				if (smooth)
					m += flatten(x0, y0, x, y);
				if (m <= -1.8)
					m = -0.8;
				m = 0.23 * log(1.8 + m);
				color(m, subpixel);
			}
			for (int c = 0; c < 3; c++)
				sum[c] += subpixel[c];
		}}
		
		// subpixels finished -> make arithmetic mean
		char pixel[3];
		for (int c = 0; c < 3; c++)
			pixel[c] = (int)(255.0 * sum[c] / (subpix * subpix)  + 0.5);
		fwrite(pixel, 1, 3, image_file);
		//pixel finished
	}}
	fclose(image_file);
	return(0);
}